基于混合變異策略差分進(jìn)化算法的邊坡滑裂面搜索研究

張子映， 柴軍瑞， 張書濱， 錢武文

(1.西安理工大學(xué) 西北旱區(qū)生態(tài)水利工程國家重點實驗室培育基地， 陜西 西安 710048；2.江西省河道湖泊管理局， 江西 南昌 330009)

1 研究背景

邊坡滑動面搜索是邊坡穩(wěn)定計算中的一個關(guān)鍵問題，其實質(zhì)為安全系數(shù)[1-4]與最危險滑動面[5-7]的搜索問題。在這一領(lǐng)域出現(xiàn)了較多的優(yōu)化搜索方法，模式搜索法、動態(tài)規(guī)劃法、數(shù)學(xué)歸納法、變分法、隨機(jī)搜索法等。隨著優(yōu)化方法研究的深入，人工智能搜索算法如遺傳算法[8-9]、模擬退火算法、粒子群優(yōu)化算法、蟻群算法[10]、神經(jīng)網(wǎng)絡(luò)以及蒙特卡洛法[11]等被用于邊坡最危險滑動面搜索，且目前已成為滑動面搜索研究中的主流方法。這些智能尋優(yōu)法各有其優(yōu)點但也都存在不足之處，例如：遺傳算法存在早熟收斂、算法性能對參數(shù)選取強(qiáng)依賴性而參數(shù)不易確定等缺陷[9]；模擬退火算法在實際應(yīng)用中由于退火不充分并不能保證一定能搜索到全局最優(yōu)值；粒子群優(yōu)化算法[12-13]易陷入局部最優(yōu)；蟻群算法[14-15]最初被設(shè)計用于組合優(yōu)化問題，其缺少用于連續(xù)數(shù)值優(yōu)化的操作算子，用于邊坡滑動面搜索問題時，需要對邊坡進(jìn)行一定精度的離散,操作過程復(fù)雜，同時也存在易于陷入局部最優(yōu)的問題，有待進(jìn)一步改進(jìn)。

差分進(jìn)化(differential evolution，簡稱 DE)算法是由 Storn等[16]于1995年提出的一種智能優(yōu)化算法，該方法具有操作簡單、收斂性好、全局尋優(yōu)能力強(qiáng)的優(yōu)點[17-19]。錢武文等[20]受Nelder-Mead方法啟發(fā)，提出了一種名為“反射變異”的新型變異策略，具有更明確的變異方向，可即使基向量擁有一定的隨機(jī)性，但仍為一種貪婪搜索策略，為進(jìn)一步增強(qiáng)其全局收斂性能，本文在此基礎(chǔ)上引入一個基本變異策略形成混合變異策略，并將該優(yōu)化方法用于邊坡臨界滑動面的搜索中，通過兩個經(jīng)典考核算例，驗證了該方法應(yīng)用到邊坡穩(wěn)定性研究是合理可行的。

2 改進(jìn)的差分進(jìn)化算法

算法的全局搜索能力和局部開采能力是相互對立的，變異策略高度影響著差分進(jìn)化算法的性能。參考文獻(xiàn)[20]提出了一種名為“反射變異”的新型變異策略，并與參數(shù)自適應(yīng)調(diào)整方法結(jié)合，形成了一種具有良好全局探測能力和局部開采能力的差分進(jìn)化算法(RMADE)。考慮到“反射變異”從某種程度上仍屬于貪婪變異策略，仍易陷入局部最優(yōu)，因此為進(jìn)一步改進(jìn)其全局收斂能力，本文在RMADE的基礎(chǔ)上引入一種基本變異策略組成混合變異策略共同控制子代個體的生成。

2.1 混合變異策略

受Nelder-Mead方法啟發(fā)，文獻(xiàn)[20]提出了一種名為反射變異策略的新型變異策略，該策略不論求解問題的維數(shù)，始終選擇種群中的4個任意個體組成單純形，每完成一次變異操作須引入一個單純形。該策略操作如下:

Vi,g=Xo,g+F(Xa,g-Xd,g)

(1)

式中：Vi,g為i個體在g代的試向量；Xo,g為單純形中心；Xa,g和Xd,g分別為單純形中的最好點和最壞點，F(xiàn)為變異縮放因子。

不同于Nelder-Mead方法，本文使用的單純形中心Xo,g是根據(jù)單純形頂點的函數(shù)值計算，并且除最壞點外的各個頂點均按其函數(shù)值權(quán)重分配到中心中。單純形中心Xo,g的計算公式如下：

(2)

式中：w1，w2和w3分別為除最壞點外的其余單純形頂點的權(quán)重；Xa,g，Xb,g和Xc,g為除最壞點外的剩余單純形頂點。

這種變異策略類似于Nelder-Mead方法中的反射操作，具有明確的搜索方向。即使基向量擁有一定的隨機(jī)性，但仍為一種貪婪搜索策略，為進(jìn)一步增強(qiáng)其全局收斂性能，本文還引入一種基本子代生成策略(DE/current-to-rand/1)，其變異向量生成如下。

Vi,g=Xi,g+F(Xc,g-Xi,g+Xa,g-Xb,g)

(3)

式中：Xi,g為當(dāng)前操作個體，其他同上。

引入DE/current-to-rand/1后，一個種群中將同時存在兩種變異策略，本文根據(jù)各策略的成功率(成功更新子代個體的比率)為各個個體輪盤賭選擇其對應(yīng)的變異策略。成功率的計算公式如下：

(4)

式中：Nk為種群中使用策略k的個體總數(shù)；srk為策略k的成功率，在初始種群時將其設(shè)置為0.5，以保證算法開始時各變異策略擁有總體相等的個體數(shù)。

2.2 控制參數(shù)自適應(yīng)調(diào)整機(jī)制

傳統(tǒng)的差分進(jìn)化算法采用固定的變異因子和交叉因子，而針對不同的優(yōu)化問題，使用固定的控制參數(shù)將嚴(yán)重妨礙算法的性能。

參數(shù)自適應(yīng)調(diào)整機(jī)制允許每一個體使用不同于其他個體的縮放因子F和交叉因子CR去生成子代個體，且能根據(jù)算法進(jìn)化過程中成功個體的經(jīng)驗反饋調(diào)節(jié)控制參數(shù)。本文采用SHADE[19]算法中的參數(shù)自適應(yīng)調(diào)整方法。

2.3 算法的基本思想和步驟

本文在文獻(xiàn)[20]中算法的基礎(chǔ)上，引入基本的變異策略DE/current-to-rand/1，從而進(jìn)一步改進(jìn)算法的全局搜索能力，并在降低早熟收斂可能性的同時保持原算法的收斂速度。使用自適應(yīng)參數(shù)調(diào)節(jié)方法進(jìn)一步提升了算法的整體性能。算法的具體步驟如下：

Step 1：初始化參數(shù)和初始種群；

Step 2：判定是否達(dá)到終止條件,若是轉(zhuǎn)Step 7，若否則轉(zhuǎn)Step 3；

Step 3：為每個個體分配變異因子和交叉因子；

Step 4：根據(jù)各變異策略的成功率使用輪盤賭方法為各個體分配變異策略；

Step 5：對每一靶個體，按公式(1)或(3)完成變異操作，然后完成交叉操作；

Step 6：計算試向量的函數(shù)值并與靶個體比較，選擇較優(yōu)個體替換靶個體；

Step 7：輸出最優(yōu)個體。

3 基于DE的邊坡穩(wěn)定分析

3.1 滑動面上任意點應(yīng)力張量計算

為克服極限平衡法分析結(jié)果不考慮土體應(yīng)力應(yīng)變狀態(tài)的缺陷，本文采用有限元軟件進(jìn)行邊坡極限平衡分析，即先計算獲得邊坡高斯應(yīng)力場，然后根據(jù)雙線性插值計算得到各單元節(jié)點的應(yīng)力值。

(5)

式中：Ni為雙線性插值函數(shù)；ξ，η為自然坐標(biāo)，本文采用笛卡爾坐標(biāo)系。

生成滑動面y=y(x)，對滑動面離散成多個線單元。依次對滑動面上每個線單元的每個高斯點查找其歸屬的有限元單元，由這個單元的各個節(jié)點的方向應(yīng)力形函數(shù)插值得到高斯點的方向應(yīng)力，即σx，σy，τxy。

(6)

3.2 安全系數(shù)計算

(8)

式中：c、φ分別為材料的凝聚力和內(nèi)摩擦角。σn、τ分別為曲線上任意一點土體方向應(yīng)力和這點位置沿曲線的剪應(yīng)力。

若將滑動面視為由一系列連續(xù)線段組成，將滑裂面離散化處理，應(yīng)用高斯積分得邊坡安全系數(shù)為:

(9)

3.3 任意滑裂面構(gòu)造方法

滑裂面由n個控制點{(x1,y1) , (x2,y2),…, (xn,yn) }組成的n-1個小片段近似構(gòu)造。由于滑入點和滑出點都位于邊坡輪廓線上，其縱坐標(biāo)y1和yn都可據(jù)已知的輪廓線求出，因此通過事先固定中間控制點的橫坐標(biāo)，進(jìn)一步減少搜索變量的個數(shù)，即由n個控制點近似的滑裂面的搜索變量為S={x1,y2,…,yn-1,xn}。

合理的滑裂面并非這些點在平面內(nèi)任意分布的組合，通常將這些片段組成的滑裂面向上凹時視為合理，如圖1所示，這意味著滑裂面通常需要滿足以下特征[21]，即：

α1≥α2≥…≥αn-1

(10)

式中：αi為線段i的傾角；n為滑裂面控制點的個數(shù)。

式(10)亦可用線段斜率ki表示，即：

k1≥k2≥…≥kn-1

(11)

使用點坐標(biāo)表示后，有：

(12)

則優(yōu)化問題變?yōu)閹в衝-2個約束的優(yōu)化問題，優(yōu)化變量為S={x1,y2,…,yn-1,xn}，即：

minF(S),S={x1,y2,…，yn-1,xn}

(13)

j=1,2,…,n-2;li≤Si≤ui；

i=1,2,…,n

在智能優(yōu)化領(lǐng)域，約束的存在易造成可行域的減少進(jìn)而使得解的搜索變得更加復(fù)雜。邊坡臨界滑裂面的搜索屬于有約束優(yōu)化問題，考慮到一般情況不合理的滑裂面對應(yīng)的安全系數(shù)比合理滑裂面更大的特性，可將有約束優(yōu)化直接按無約束優(yōu)化處理。

圖1 合理的邊坡滑動面

3.4 求解最危險滑裂面

求解最危險滑裂面的程序流程如下(程序框圖見圖2)：

(1)根據(jù)邊坡地層分布情況及邊坡剖面界限，運用有限元軟件進(jìn)行合理的有限元網(wǎng)格劃分，并計算邊坡應(yīng)力場。

(2)差分進(jìn)化算法應(yīng)用于滑裂面搜索，每個個體代表一個滑裂面，每個滑裂面由n個控制點組成的曲線，故每個個體就是一個n維的向量Xi=[x1,y2,…,yn-1,xn]。

(3)在算法搜索前，根據(jù)剖面的幾何模型對個體的每個維度設(shè)定合理的搜索區(qū)域。其控制方式為：第一個滑入點和最后一個滑出點的搜索范圍位于坡表，中間控制點縱坐標(biāo)為搜索區(qū)間隨機(jī)初始化值y2,y3,…yn-1。

(4)構(gòu)建初始種群{X1,X2,…,XNP}后，使用有限元極限平衡法計算個體Xi的安全系數(shù)Ki，用差分進(jìn)化算法搜索出最佳滑裂面Xbest。

圖2 基于差分進(jìn)化算法的有限元極限平衡法的流程圖

4 算例驗證

本算例使用8個滑裂面控制點(滑入點、滑出點以及6個中間控制點)，因此只需搜索滑入點和滑出點的X值以及3個中間控制點的Y坐標(biāo)值。為了使中間控制點更好地反映滑面的形狀，盡可能地將中間控制點的橫坐標(biāo)均勻分布在搜索域內(nèi)。種群規(guī)模設(shè)置為40，最大迭代次數(shù)設(shè)為150。

4.1 均質(zhì)邊坡

引用澳大利亞計算機(jī)應(yīng)用協(xié)會(ACADS)通用考題1，坡高為10 m，坡角30°，土容重γ=20 kN/m3，黏聚力c=3 kPa，土的彈性模量E=104kN/m2，內(nèi)摩擦角φ=19.6°，泊松比ν=0.25。推薦的安全系數(shù)為1.0。應(yīng)力計算模型如圖3所示，假定區(qū)域左右邊界為法向約束，底部邊界為固定端約束。

建立有限元模型，使用理想彈塑性本構(gòu)模型計算基于Mohr-Coulomb強(qiáng)度破壞準(zhǔn)則的邊坡土體的應(yīng)力場分布，進(jìn)而以本文方法確定最危險滑動面的位置和其對應(yīng)的安全系數(shù)。為了進(jìn)行比較，將結(jié)果與使用傳統(tǒng)極限平衡法spencer、Bishop和Janbu法穩(wěn)定性分析的結(jié)果比較，將各方法所得到的安全系數(shù)和最危險滑裂面匯總繪制于圖4中。由圖4可知，圓弧滑動搜索結(jié)果為KJanbu=0.966，KBishop=0.996，KSpencer=0.995，推薦安全系數(shù)為1.0。本文方法搜索的邊坡最危險滑面安全系數(shù)為1.055，僅與推薦的安全系數(shù)相差5.5%，在合理范圍，本文方法正確。

圖3 邊坡幾何模型

圖4 各方法搜索的最危險滑裂面對比

4.2 非均質(zhì)邊坡

引用澳大利亞計算機(jī)應(yīng)用協(xié)會(ACADS)一通用考題。坡體由3種土料構(gòu)成，坡高為 10 m，坡比為1∶2，邊坡幾何特征如圖5所示，其材料特性如表1所示，推薦的安全系數(shù)為1.0。

圖5 邊坡幾何模型

表1 非均質(zhì)邊坡材料性質(zhì)

建立有限元模型，使用理想彈塑性本構(gòu)模型計算基于Mohr-Coulomb強(qiáng)度破壞準(zhǔn)則的邊坡土體的應(yīng)力場分布，進(jìn)而以本文方法確定最危險滑動面的位置和其對應(yīng)的安全系數(shù)。同上，將各方法所得到的安全系數(shù)和最危險滑裂面匯總繪制于圖6中。

由圖6可知，圓弧滑動搜索結(jié)果為KJanbu=1.341，KBishop=1.442，KSpencer=1.419，推薦安全系數(shù)為1.39。本文方法搜索的邊坡最危險滑面安全系數(shù)為1.375，比推薦的安全系數(shù)小1.08%，在合理范圍。本文計算結(jié)果處于Janbu法與Bishop法和Spencer法之間，且略小于Bishop法和Spencer法，說明本文方法的準(zhǔn)確性?；谟邢拊獜椝苄宰冃蔚臉O限平衡法滑面起始部分較圓弧法略有“上翹”，尾部略有“下沉”，考慮到土體的彈塑性變形這種現(xiàn)象是可以解釋的。說明是否考慮塑性變形對邊坡的滑面位置和安全系數(shù)影響較大。

圖6 各方法搜索出的最危險滑裂面對比

4.3 收斂性能比較

使用上節(jié)非均質(zhì)邊坡作為計算模型，將本文改進(jìn)的差分進(jìn)化算法與單純形法、DE、RMADE進(jìn)行比較，研究不同算法對邊坡臨界滑動面搜索的影響。鑒于單純形法為無約束優(yōu)化算法，而滑面受到邊坡坡體區(qū)域的限制，因此本文對單純形法中的初始單純形生成方法進(jìn)行了略微調(diào)整，即引入差分進(jìn)化算法中的初始種群技術(shù)生成初始單純形。

文中使用的單純形法的反射系數(shù)、擴(kuò)張系數(shù)以及收縮系數(shù)分別為1.0、1.3和0.7，DE的縮放因子F=0.8，交叉因子CR=0.8。種群規(guī)模NP=40，最大迭代次數(shù)為200。為了研究算法搜索臨界滑面的穩(wěn)定性，各算法均獨立運行10次，記錄最大值、平均值、最小值以及標(biāo)準(zhǔn)差于表2中，各算法的典型收斂過程線繪于圖7中。

表2 獨立10次計算下各算法性能比較

圖7 各算法的收斂過程

由表2可知，單純形法結(jié)果最為發(fā)散，說明其已陷入局部最優(yōu)，本文算法的平均值和標(biāo)準(zhǔn)差最小，說明其性能較其他方法要好且更穩(wěn)定。由圖7可知，本文算法擁有較快的收斂速度。綜合可知，本文算法性能優(yōu)于其他幾種方法，說明本文對RMADE的改進(jìn)具有實際意義。

5 結(jié) 論

(1)本文在RMADE的基礎(chǔ)上引入一個基本的變異策略DE/current-to-rand/1而形成了一種改進(jìn)的差分進(jìn)化算法，該方法根據(jù)各變異策略的子代成功率使用輪盤賭選擇自動選擇個體變異策略，并與其他優(yōu)化方法進(jìn)行比較驗證了本文提出的算法具有更好的尋優(yōu)性能。

(2)將本文改進(jìn)的差分進(jìn)化算法與有限元極限平衡法結(jié)合，應(yīng)用于邊坡穩(wěn)定性分析中，通過兩個算例驗證了該法解決邊坡穩(wěn)定性問題的有效性。