本發(fā)明屬于巖體變形和受力狀態(tài)分析的巖土工程問題，具體涉及一種實(shí)現(xiàn)工程巖體力學(xué)參數(shù)反演的改進(jìn)粒子群優(yōu)化方法。

背景技術(shù)：

工程巖體變形和受力狀態(tài)分析的難題之一是如何恰當(dāng)?shù)毓烙?jì)巖土體的力學(xué)模型參數(shù)。最優(yōu)估計(jì)的巖土體力學(xué)模型參數(shù)是通過比較現(xiàn)場觀測到的信息數(shù)據(jù)與理論模型得到的。通過定義目標(biāo)函數(shù)，將參數(shù)識別反問題轉(zhuǎn)化為優(yōu)化問題處理。伴隨著計(jì)算機(jī)技術(shù)的發(fā)展，正分析的理論和計(jì)算方法逐漸成熟，觀測儀器的精度也逐步提高，根據(jù)現(xiàn)場觀測數(shù)據(jù)進(jìn)行巖土力學(xué)模型參數(shù)反演具有良好的應(yīng)用前景。目前粒子群算法中的適應(yīng)值在算法前期改善較快，而在中后期增加緩慢，迭代次數(shù)過多，占用較多運(yùn)行時(shí)間，效率較低。

技術(shù)實(shí)現(xiàn)要素：

發(fā)明目的：本發(fā)明的針對現(xiàn)有技術(shù)中存在的問題，提供一種實(shí)現(xiàn)工程巖體力學(xué)參數(shù)反演的改進(jìn)粒子群優(yōu)化方法。

技術(shù)方案：本發(fā)明所述的一種實(shí)現(xiàn)工程巖體力學(xué)參數(shù)反演的改進(jìn)粒子群優(yōu)化方法，包括下述步驟：

(1)根據(jù)邊坡工程巖體的已知性能參數(shù)設(shè)置改進(jìn)粒子群優(yōu)化算法參數(shù)；

(2)隨機(jī)初始化每個(gè)粒子的速度與位置，定義反演適應(yīng)度函數(shù)，計(jì)算各粒子的適應(yīng)值，根據(jù)適應(yīng)值評價(jià)各粒子位置的優(yōu)劣程度；

(3)進(jìn)入快速粒子群搜索階段，不斷更新當(dāng)前粒子的個(gè)體極值及全局極值，獲取粒子個(gè)體最優(yōu)位置和粒子種群最優(yōu)位置；

(4)檢驗(yàn)是否滿足直接搜索引入準(zhǔn)則，如滿足，進(jìn)入直接搜索階段，否則，返回到上一步驟繼續(xù)快速粒子群搜索迭代；

(5)進(jìn)入直接搜索階段，將快速粒子群搜索階段得到的最優(yōu)解p0作為直接搜索階段的初始解，粒子種群最優(yōu)適應(yīng)值pbest作為直接搜索階段的初始代最佳適應(yīng)值，進(jìn)行精度搜索；

(6)引入改進(jìn)算法終止判斷準(zhǔn)則，如滿足，迭代優(yōu)化終止，此時(shí)輸出的最優(yōu)解p0即為反演所得最優(yōu)的邊坡工程力學(xué)參數(shù)；如不滿足，返回上一步驟繼續(xù)搜索。

上述步驟(1)中，設(shè)置的改進(jìn)粒子群優(yōu)化算法參數(shù)包括粒子種群大小m、慣性權(quán)重ω、加速因子c1、c2、解的維數(shù)d、每一維最大速度vimax、每一維上的解上下限ximax、ximin、初級精度ε1、高級精度ε。

步驟(2)中，適應(yīng)度函數(shù)為式中，fi(x)為第i個(gè)量測點(diǎn)上發(fā)生的位移量的計(jì)算值，ui為相應(yīng)的位移量的實(shí)測值。f(x)的值越小，粒子適應(yīng)度越高，位置越好，解的品質(zhì)越高，也就是說該粒子對應(yīng)的一組力學(xué)參數(shù)更優(yōu)。

較優(yōu)的，在計(jì)算各粒子的適應(yīng)值前，先判斷粒子是否違反約束，對于違反約束的粒子采用隱匿墻處理；否則通過適應(yīng)度函數(shù)計(jì)算各粒子的適應(yīng)值。

優(yōu)選的，步驟(3)中，針對較差粒子施以“加時(shí)飛行策略”，根據(jù)下式更新粒子新位置：

xid^t+1＝xid^t+ti*νid^t+1；

νid^t+1＝ωνid^t+c1r1(pid^t-xid^t)+c2r2(p^tgd-xid^t)；

式中：

ti為第i個(gè)粒子的飛行時(shí)間，fi為較差粒子的適應(yīng)值，farv為粒子群體的平均適應(yīng)值；

vid^t+1為每個(gè)粒子的新速度，t為當(dāng)前進(jìn)化代數(shù)；d＝1，2，...，n，n為粒子空間搜索的維數(shù)；i＝1，2，...，m，m為種群規(guī)模；w為慣性權(quán)重；

c1、c2為加速因子，c1調(diào)節(jié)粒子向自身最優(yōu)位置飛行的步長，c2調(diào)節(jié)粒子向全局最優(yōu)位置飛行的步長；

r1、r2為分布在[0，1]之間的隨機(jī)數(shù)；

pid為第i個(gè)粒子d維上的個(gè)體最優(yōu)位置；xid為第i個(gè)候選解中第d個(gè)待優(yōu)化力學(xué)參數(shù)變量。

同時(shí)，針對最優(yōu)粒子采取“保持最優(yōu)位置策略”，群體中除較差粒子和最優(yōu)粒子以外的粒子采用正常的位置更新公式xid^t+1+1＝xid^t+νid^t+1計(jì)算其新位置。

進(jìn)一步的，本步驟中，更新粒子個(gè)體極值時(shí)，若粒子的當(dāng)前適應(yīng)值較其歷史最好個(gè)體極值小，則對其進(jìn)行更新操作，否則個(gè)體極值不變；更新當(dāng)前粒子的全局極值時(shí)，比較其當(dāng)前適應(yīng)值與全局極值，如前者小則用其更新全局極值，且全局極值對應(yīng)的參數(shù)向量x為粒子種群的當(dāng)前代最優(yōu)解。

上述步驟(4)中，直接搜索引入準(zhǔn)則包括如下兩種情形：

(1)在指定的迭代次數(shù)內(nèi)解的質(zhì)量無明顯改進(jìn)，即個(gè)體最優(yōu)位置幾乎不變，出現(xiàn)了解的收斂停滯現(xiàn)象；

(2)粒子的適應(yīng)值滿足時(shí)，ε1為設(shè)定的初級精度；

當(dāng)滿足上述兩個(gè)條件之一時(shí)，算法可進(jìn)入直接搜索階段。

具體的，步驟(5)中，在最優(yōu)解p0對應(yīng)的最優(yōu)反演參數(shù)組合的附近空間內(nèi)生成m-1個(gè)新的力學(xué)參數(shù)組，利用直接搜索法進(jìn)行精度搜索，range線性遞減，并計(jì)算每一個(gè)力學(xué)參數(shù)組合方案的適應(yīng)值，若該值比種群最優(yōu)適應(yīng)值pbest小，則用該值更新pbest，最終將適應(yīng)值最小的一組參數(shù)更新為最優(yōu)解p0；其中，range為優(yōu)化參數(shù)變化范圍矢量，線性遞減range＝(1-α)×range，其中α為優(yōu)化參數(shù)的領(lǐng)域搜索范圍線性遞減參數(shù)。

上述步驟(6)中，算法終止判斷準(zhǔn)則包括下述兩種情形：

(1)最大直接搜索代數(shù)k或在指定的代數(shù)內(nèi)解的質(zhì)量無明顯改進(jìn)；

(2)粒子適應(yīng)度函數(shù)值ε1為設(shè)定的初級精度；

當(dāng)滿足上述兩個(gè)條件之一時(shí)，算法即可終止。

有益效果：與現(xiàn)有技術(shù)相比，本發(fā)明的優(yōu)點(diǎn)在于：本發(fā)明對粒子群優(yōu)化算法進(jìn)行改進(jìn)，先通過快速粒子群搜索獲取粒子個(gè)體最優(yōu)位置和粒子種群最優(yōu)位置，針對不同位置的粒子實(shí)行不同的位置更新方式，然后進(jìn)一步通過精度搜索、range線性遞減得到最優(yōu)解，最終反演得到最優(yōu)的邊坡工程力學(xué)參數(shù)；將該方法應(yīng)用于工程巖體力學(xué)參數(shù)反演，可以提高優(yōu)化反演效率，減少迭代次數(shù)，減少正算計(jì)算量。

附圖說明

圖1為實(shí)施例中邊坡開挖前后的二維離散模型，其中，(a)為開挖前的二維離散模型，(b)為開挖后的二維離散模型；

圖2為實(shí)施例中應(yīng)用本發(fā)明的改進(jìn)粒子群優(yōu)化方法反演邊坡工程巖體力學(xué)參數(shù)過程中粒子的迭代收斂圖。

具體實(shí)施方式

下面結(jié)合附圖對本發(fā)明的技術(shù)方案作進(jìn)一步說明。

本發(fā)明的一種實(shí)現(xiàn)工程巖體力學(xué)參數(shù)反演的改進(jìn)粒子群優(yōu)化方法，包括下述步驟：

(1)設(shè)置改進(jìn)粒子群優(yōu)化算法參數(shù)；根據(jù)工程巖體的已知性能參數(shù)設(shè)置改進(jìn)粒子群優(yōu)化算法的參數(shù)，設(shè)置的參數(shù)包括粒子種群大小m、慣性權(quán)重ω、加速因子c1、c2、解的維數(shù)d、每一維最大速度vimax、每一維上的解上下限ximax、ximin、初級精度ε1和高級精度ε；

(2)隨機(jī)初始化每個(gè)粒子的速度與位置，定義反演適應(yīng)度函數(shù)，計(jì)算各粒子的適應(yīng)值，通過適應(yīng)值評價(jià)各粒子位置的優(yōu)劣程度；

反演適應(yīng)度函數(shù)為：fi(x)為第i個(gè)量測點(diǎn)上發(fā)生的位移量的計(jì)算值，ui為相應(yīng)的位移量的實(shí)測值。f(x)的值越小代表粒子適應(yīng)度越高，位置越好，解的品質(zhì)越高，也就是說該粒子對應(yīng)的一組力學(xué)參數(shù)更優(yōu)。

本步驟中，在計(jì)算各粒子的適應(yīng)值前，要判斷粒子是否違反約束，粒子越出解空間的邊界即判定粒子違反約束；對于違反約束的粒子采用隱匿墻處理，隱匿墻處理是指粒子允許越過物理限定區(qū)域，但是越界的粒子并不估計(jì)其適應(yīng)度；否則通過適應(yīng)度函數(shù)計(jì)算各粒子的適應(yīng)值。

(3)進(jìn)入快速粒子群搜索階段：不斷更新當(dāng)前粒子的個(gè)體極值及全局極值，獲取粒子個(gè)體最優(yōu)位置和粒子種群最優(yōu)位置；

具體而言，根據(jù)粒子的適應(yīng)值判斷粒子的優(yōu)差，適應(yīng)度數(shù)值越小，認(rèn)為其對應(yīng)的適應(yīng)度越好，所在位置就越好，所代表的一組反演參數(shù)就越優(yōu)。

針對較差粒子施以“加時(shí)飛行策略”，根據(jù)下式更新粒子新位置xid^t+1＝xid^t+ti*νid^t+1，

其中，第i個(gè)粒子的飛行時(shí)間fi為較差粒子適應(yīng)值，farv為粒子群體平均適應(yīng)值；

每個(gè)粒子的新速度νid^t+1＝ωνid^t+c1r1(pid^t-xid^t)+c2r2(p^tgd-xid^t)，式中，t為當(dāng)前進(jìn)化代數(shù)；d＝1，2，...，n，n為粒子空間搜索的維數(shù)；i＝1，2，...，m，m為種群規(guī)模；w為慣性權(quán)重；c1、c2為加速因子，c1調(diào)節(jié)粒子向自身最優(yōu)位置飛行的步長，c2調(diào)節(jié)粒子向全局最優(yōu)位置飛行的步長；r1、r2為分布在[0，1]之間的隨機(jī)數(shù)；pid為第i個(gè)粒子d維上的個(gè)體最優(yōu)位置；xid為第i個(gè)候選解中第d個(gè)待優(yōu)化力學(xué)參數(shù)變量。

針對最優(yōu)粒子采取“保持最優(yōu)位置策略”，群體中除此以外的粒子采用正常的位置更新公式xid^t+1+1＝xid^t+νid^t+1計(jì)算其新位置；

得到每個(gè)粒子的新速度和新位置后，再次計(jì)算其適應(yīng)度值，根據(jù)適應(yīng)度值更新當(dāng)前粒子的個(gè)體極值及全局極值；其中，更新粒子個(gè)體極值時(shí)，若粒子的當(dāng)前適應(yīng)值較其歷史最好個(gè)體極值小，則對其進(jìn)行更新操作，否則個(gè)體極值不變；更新當(dāng)前粒子的全局極值時(shí)，比較當(dāng)前適應(yīng)值與全局極值，如前者小則用其更新全局極值，且全局極值對應(yīng)的參數(shù)向量x為粒子種群的當(dāng)前代最優(yōu)解。

(4)檢驗(yàn)是否滿足直接搜索引入準(zhǔn)則，如滿足，進(jìn)入直接搜索階段，否則，返回到上一步驟繼續(xù)快速粒子群搜索迭代；

當(dāng)滿足下述兩條件之一時(shí)，代表滿足直接搜索引入準(zhǔn)則，算法就可以進(jìn)入直接搜索階段：第一種，在指定的迭代次數(shù)內(nèi)解的質(zhì)量無明顯改進(jìn)，即個(gè)體最優(yōu)位置幾乎不變，出現(xiàn)了解的收斂停滯現(xiàn)象；第二種，粒子適應(yīng)度值時(shí)，ε1為設(shè)定的初級精度，是一個(gè)較小的正數(shù)。

(5)進(jìn)入直接搜索階段：將快速粒子群搜索階段得到的最優(yōu)解p0作為直接搜索階段的初始解，粒子種群最優(yōu)適應(yīng)值pbest作為直接搜索階段的初始代最佳適應(yīng)值，進(jìn)行精度搜索；

具體而言，在最優(yōu)解p0周圍一定范圍range內(nèi)隨機(jī)生成m-1個(gè)解(即在前期優(yōu)化m迭代所得最優(yōu)反演參數(shù)組合的附近空間內(nèi)生成m-1個(gè)新的力學(xué)參數(shù)組開始利用直接搜索法進(jìn)行精度搜索，range線性遞減)，并計(jì)算每一個(gè)力學(xué)參數(shù)組合方案的適應(yīng)值，與種群最優(yōu)適應(yīng)值pbest比較起來較小，則用其值更新pbest，將適應(yīng)值最小的一組參數(shù)更新為最優(yōu)解p0。

每一組反演參數(shù)組合包括下述參數(shù)：粒子種群大小m、慣性權(quán)重ω、加速因子c1、c2、解的維數(shù)d、每一維最大速度vimax、每一維上的解上下限ximax、ximin、初級精度ε1、高級精度ε。range為優(yōu)化參數(shù)變化范圍矢量，線性遞減range＝(1-α)×range，其中α為優(yōu)化參數(shù)的領(lǐng)域搜索范圍線性遞減參數(shù)。

(6)引入改進(jìn)算法終止判斷準(zhǔn)則，如滿足，迭代優(yōu)化終止，此時(shí)輸出的最優(yōu)解p0即為反演所得最優(yōu)的邊坡工程力學(xué)參數(shù)；如不滿足，返回上一步驟繼續(xù)搜索。

當(dāng)滿足下述兩條件之一時(shí)，代表滿足算法終止判斷準(zhǔn)則，算法可以停止：第一種，最大直接搜索代數(shù)k或在指定的代數(shù)內(nèi)解的質(zhì)量無明顯改進(jìn)；第二種，粒子適應(yīng)度值

實(shí)施例

如圖1，邊坡分三步開挖(邊坡分上下兩種介質(zhì))，高為50m，上兩級坡比為1﹕1，開挖高度為25m，平臺(tái)為5m；第三級開挖坡比為1﹕0.5，開挖高度為20m。上層介質(zhì)30m，下層介質(zhì)20m。

考慮分步開挖卸荷對邊坡變形的影響，圖1中(a)圖和(b)圖分別給出了其開挖前后的離散模型，反分析的基礎(chǔ)數(shù)據(jù)為設(shè)在各坡角處觀測點(diǎn)1，2，3的豎向位移。本構(gòu)模型采用摩爾-庫侖彈塑性模型，邊坡材料參數(shù)值如表1所示。

表1材料參數(shù)表

先對該邊坡用flac進(jìn)行正演計(jì)算，以三個(gè)測點(diǎn)1、2、3的開挖完成時(shí)的豎向位移(向上為正)3.90mm，12.8mm，3.72mm作為真實(shí)觀測值，進(jìn)行參數(shù)反演計(jì)算，待反演的參數(shù)共2個(gè)，分別為上層介質(zhì)變形模量e1和下層介質(zhì)變形模量e2，其參數(shù)反演范圍分別為0.2～0.8gpa，4.0～8.0gpa。

以上表1中內(nèi)聚力c，內(nèi)摩擦角φ、泊松比μ和介質(zhì)容重γ作為初始參數(shù)，采用本發(fā)明的改進(jìn)粒子群優(yōu)化方法對邊坡工程的力學(xué)參數(shù)進(jìn)行反演，反演過程中粒子的迭代收斂圖如圖2，反演得到的一組最優(yōu)邊坡工程力學(xué)參數(shù)如下表2。

表2邊坡巖體變形模量反演結(jié)果

同時(shí)采用標(biāo)準(zhǔn)粒子群優(yōu)化算法和線性遞減慣性權(quán)重的粒子群優(yōu)化算法對實(shí)施例中的邊坡工程力學(xué)參數(shù)進(jìn)行反演，三種粒子群優(yōu)化算法的運(yùn)行情況如下表3，其中，fpso-ds代表本發(fā)明的改進(jìn)粒子群優(yōu)化算法，pso-ldw代表標(biāo)準(zhǔn)粒子群優(yōu)化算法；spso代表線性遞減慣性權(quán)重的粒子群優(yōu)化算法。

表3：三種粒子群優(yōu)化算法的運(yùn)行情況

由表3可知，與現(xiàn)有的標(biāo)準(zhǔn)粒子群優(yōu)化算法和線性遞減慣性權(quán)重粒子群優(yōu)化算法相比，本發(fā)明的實(shí)現(xiàn)邊坡工程力學(xué)參數(shù)反演的改進(jìn)粒子群優(yōu)化算法可以有效提高優(yōu)化反演效率，減少迭代次數(shù)，減少正算計(jì)算量。