王述紅， 魏 崴， 韓文帥， 陳 浩

(東北大學 資源與土木工程學院， 遼寧 沈陽 110819)

邊坡最小安全系數(shù)是評估邊坡安全狀態(tài)和確定最危險滑動面的重要量化指標，在邊坡穩(wěn)定性分析中具有重要的研究意義.求解邊坡的最小安全系數(shù)可以轉(zhuǎn)換為搜索在一定約束條件下相應(yīng)優(yōu)化問題的全局最優(yōu)解.由于受到邊坡結(jié)構(gòu)的多樣性、滑裂面形式和組成邊坡巖土體材料賦存狀態(tài)等不定因素的影響，該優(yōu)化問題往往是典型的非凸優(yōu)化問題，表現(xiàn)出強烈的非線性和復雜性[1].

伴隨著計算機領(lǐng)域的快速發(fā)展，各種優(yōu)化算法例如遺傳算法[2-4]、粒子群算法[5-7]、模擬退火算法[8-9]等已被成功應(yīng)用于邊坡工程的多個領(lǐng)域，均不同程度地推進了邊坡穩(wěn)定性分析方法的完善和發(fā)展.針對不同類型的邊坡問題，每種算法均有其特定的優(yōu)勢屬性，但由于單一算法的局限性，同樣存在不足之處.例如，遺傳算法子代數(shù)量與父代個體適應(yīng)度呈現(xiàn)典型的正相關(guān)，容易發(fā)生早熟現(xiàn)象；粒子群算法較容易受到局部最優(yōu)解的牽引而趨于同一性；模擬退火算法的尋優(yōu)性能較大程度上依賴于退火程度和冷卻速率的參數(shù)影響.灰狼優(yōu)化(grey wolf optimization，GWO)算法的靈感源于自然界中狼群狩獵機制的元啟發(fā)式算法，具有結(jié)構(gòu)簡單、需調(diào)節(jié)參數(shù)少、易于實現(xiàn)等特點[10].函數(shù)優(yōu)化測試表明[11]，其算力與遺傳算法、粒子群優(yōu)化等主流優(yōu)化算法相當，目前已廣泛應(yīng)用于多個學科的各個領(lǐng)域[12-14].但是在處理相對復雜的尋優(yōu)問題時，基礎(chǔ)GWO算法同樣存在尋優(yōu)停滯、陷入局部最優(yōu)陷阱的情況.

針對傳統(tǒng)灰狼算法存在初始種群不均勻、早熟收斂等問題，本文根據(jù)NFL定理[15]從算法的群初始化、擾動策略、關(guān)鍵參數(shù)調(diào)節(jié)等方面進行改進，提出了一種混沌灰狼優(yōu)化(chaotic grey wolf optimization，CGWO)算法用于確定二維圓弧滑裂面邊坡的最小安全系數(shù)及其臨界滑裂面.改進后的算法具有較強的尋優(yōu)性能和魯棒性，克服了基本算法在應(yīng)對復雜問題時收斂效率低、容易陷入局部最優(yōu)解等缺陷.選取圓弧滑裂面與邊坡出滑點的橫坐標XL和XR以及滑裂面圓心的橫坐標X0作為設(shè)計變量,改善了以往搜索區(qū)域設(shè)置困難、主觀性強等情況.選取ACADS邊坡考核題對該方法的有效性進行驗證，并與已有的算法及有限元強度折減法進行對比，驗證了本文方法的準確性和可靠性.同時，該方法適用性強、搜索區(qū)域易于設(shè)置，具有一定的工程指導和推廣意義.

1 計算模型

1.1 邊坡穩(wěn)定性分析力學模型

簡化的Bishop方法是我國《公路路基設(shè)計規(guī)范(JTGD30—2004)》推薦使用的方法，已被嚴格證明具有較高的計算精度[16]，在實際工程中被廣泛應(yīng)用于評估邊坡的穩(wěn)定性.該方法假定邊坡潛在的滑裂面為圓弧滑裂面，各土條之間只存在水平推力.將滑裂體垂直分成若干塊，第i塊土體的受力分析如圖1所示.

考慮各土條之間的垂直平衡及對圓心的力矩平衡，當圓弧滑裂面的控制參數(shù)給定以后，邊坡的安全系數(shù)表達式為

(1)

(2)

式中：Fs為邊坡安全系數(shù)；mi為滑動面第i個土條的質(zhì)量；φi，ci分別為第i個土條在滑動面上土體的內(nèi)摩擦角和黏聚力；bi為第i個土條沿滑動面的長度；θi為第i個土條底邊傾角；pi為作用于第i個土條邊上的孔隙水壓力.

圖1 簡化的Bishop方法力學模型

1.2 滑裂面搜索模型

在構(gòu)建邊坡安全系數(shù)搜索模型時，多選取圓弧滑裂面的圓心坐標(X0,Y0)和半徑R作為設(shè)計變量[2-3,8]，如圖2所示，控制變量用紅色字體標出.

圖2 文獻[2-3,8]采用的優(yōu)化模型

該方法所得最優(yōu)解無需轉(zhuǎn)換，尋優(yōu)結(jié)果即為控制滑裂面位置形態(tài)的幾何參數(shù).但在實際工程中，邊坡面形狀往往比較復雜，滑裂面圓弧與邊坡面左、右交點的坐標不易準確求解.不管是采用近似解法還是精確解法，都會對后續(xù)條分法的精度和計算速度產(chǎn)生影響.另外，該設(shè)計變量的搜索區(qū)域的標定相對復雜.對于不同類型的邊坡，要求使用者憑借自身經(jīng)驗和查閱工程地質(zhì)資料來確定大致的搜索范圍[8]，并通過多次試算來壓縮搜索區(qū)域以減少算法的尋優(yōu)壓力.當搜索區(qū)域設(shè)置不當時，容易遺漏近坡面小半徑及遠坡面大半徑等可能潛在的滑裂面形式，如圖2藍線和紅線圓弧滑裂面所示.

針對以上問題，本文選取圓弧滑裂面與邊坡出滑點的橫坐標XL和XR以及滑裂面圓心的橫坐標X0作為設(shè)計變量，建立如圖3所示的邊坡穩(wěn)定性分析計算模型.對于不同的邊坡算例，只需更改相應(yīng)的邊界控制點即可完成設(shè)計變量搜索范圍的設(shè)置.圖中，PL，PR為邊坡的左、右出滑點；Xa，Xb，Xa,b分別為滑裂面搜索區(qū)域左、右邊界控制點及中點的橫坐標.

圖3 本文采用的優(yōu)化模型

在該搜索區(qū)域內(nèi)任意選擇滿足約束條件的XL，XR及X0參數(shù)，根據(jù)坡面曲線可直接求得XL和XR，由對應(yīng)的幾何關(guān)系可以推導出

(3)

(4)

坡體巖土材料參數(shù)已知，通過改變設(shè)計變量即可實現(xiàn)對任意圓弧滑裂面的標定，則安全系數(shù)可用設(shè)計變量表達為

F=Φ(X)

,

(5)

X=[XL,XR,X0]T

.

(6)

式中：F為任意圓弧滑裂面所對應(yīng)的安全系數(shù)；Φ為設(shè)計變量X關(guān)于安全系數(shù)F的函數(shù).

1.3 安全系數(shù)優(yōu)化模型

在構(gòu)建邊坡穩(wěn)定性分析的優(yōu)化模型時，需要過濾不符合實際情況的滑裂面形式，避免圓弧滑裂面圓心在邊坡內(nèi)部或滑裂曲線越界穿過最下層土層等錯誤情況的發(fā)生.也即是設(shè)計變量需要滿足相應(yīng)的約束條件以構(gòu)造符合真實情況的潛在滑裂面.滿足約束條件后的邊坡安全系數(shù)尋優(yōu)模型最終可表示為

(7)

式中：約束條件即為搜索區(qū)域；Xc為全局最小安全系數(shù)所對應(yīng)臨界滑裂面的設(shè)計變量；h′為底部土層下邊界到地面線的距離；Yd為X0所對應(yīng)邊坡面幾何曲線的縱坐標.

2 混沌GWO算法

2.1 基本GWO算法

灰狼優(yōu)化(grey wolf optimization，GWO)算法由Mirjalili等[10]于2014年提出，是一種基于灰狼狩獵過程和社會等級的元啟發(fā)式算法.在每一代狼群中，種群領(lǐng)導者被定義為α狼，負責協(xié)助α狼進行決策的下屬被定義為β狼，地位最低的狼被稱為δ狼.不在以上三個等級(優(yōu)勢狼)中任意一級的狼被稱為ω狼(劣勢狼).

灰狼群體逐漸接近并包圍獵物的行為，用數(shù)學模型可表達為

D=|C⊙Xp(t)-X(t)|,

(8)

Xi(t+1)=Xp(t)-A⊙D.

(9)

式中：D為灰狼與獵物間距離向量；t為當前迭代的次數(shù)；Xp(t)為當前灰狼種群的位置向量；A和C為系數(shù)向量，可定義為

A=2a⊙r1-a，

(10)

C=2r2

,

(11)

(12)

式中：tmax為最大迭代次數(shù)；r1，r2均為各元素在[0,1]產(chǎn)生的隨機向量；a為距離控制參數(shù)向量，初值a0=2.“⊙”符號的含義為兩向量對應(yīng)位置的元素做乘法.

狩獵過程中，ω狼根據(jù)α狼、β狼、δ狼的指導從各方向逼近獵物的數(shù)學描述如下：

(13)

ω狼相對于α狼、β狼、δ狼的前進方向與步長為

(14)

ω狼最終的位置向量為

(15)

根據(jù)上述模型，狼群在進化過程中通過不斷更新位置來包圍獵物，直至達到最大迭代次數(shù)或求解精度.

2.2 混沌初始種群

初始種群的質(zhì)量在一定程度上影響算法的全局尋優(yōu)性能和收斂速度.基本GWO算法通過計算機生成偽隨機數(shù)在決策空間內(nèi)的初始種群，該策略在處理較復雜的問題時，會導致算法種群多樣性降低，進而影響算法的全局尋優(yōu)性能和收斂速度.而混沌系統(tǒng)作為一類同時具備周期性和不收斂性的隨機有界系統(tǒng)，其所對應(yīng)的混沌映射所具備的規(guī)則性、隨機性、遍歷性和不可預測性使之成為可靠的隨機源之一.這些特點使得混沌系統(tǒng)能夠有效地維持算法種群的多樣性，增加算法跳出局部解的概率.其中Tent混沌序列具有較好的區(qū)間均勻分布特性和迭代效率，多被應(yīng)用于豐富和維持算法的種群多樣性[17-18].

Tent映射數(shù)學表達式：

(16)

式中：xt為第t次迭代所對應(yīng)的Tent混沌映射值；u為調(diào)節(jié)參數(shù).

當u=0.5時，即為典型的Tent映射形式.單梁等[17]證明此時的Tent映射與λ=4的logistic映射相比具有更快的迭代速度和[0，1]區(qū)間均勻分布特性.同時Hung等[18]也指出，由于典型的Tent映射受字節(jié)長度的限制和計算機精度的影響，存在例如{0.2,0.4,0.6,0.8}的小周期點以及{0,0.25,0.5,0.75}的不穩(wěn)定周期點，當?shù)趖次迭代的混沌序列為以上點集合的元素時，混沌迭代容易趨于周期性或收斂到穩(wěn)定值0.

為了避免以上情況的發(fā)生，本文取u=0.8，同時考慮到混沌映射對初值具有強烈的敏感依賴性，通過多次試算確定初值x0=0.456 7.如圖4所示，改進后的Tent映射能夠有效避免收斂至不動點0，并且具有較好的均勻遍歷特性和迭代速度.

2.3 混沌擾動策略

當目標函數(shù)在極小的求解域內(nèi)存在多個界限模糊的局部最優(yōu)解時，種群容易受到局部最優(yōu)解的吸引導致過早趨于同一化.采用變異擾動策略，能夠減少這一情況發(fā)生的可能性，增強算法的全局尋優(yōu)能力.為了維持種群的豐富度，降低GWO算法過早收斂的概率，采用上述改進的Tent混沌映射對種群進行混沌擾動，位置更新公式如下：

(17)

式中,Xω(t)為第t次迭代通過式(15)所計算的當前狼群位置.

圖4 改進Tent混沌序列及分布圖

2.4 混沌非線性調(diào)節(jié)策略

由于實際待優(yōu)化問題的多樣性，基本GWO算法中參數(shù)向量a的線性遞減策略往往不能很好地適應(yīng)復雜的尋優(yōu)搜索過程.Gandomi等[19]驗證了將混沌映射應(yīng)用于算法參數(shù)調(diào)節(jié)的有效性,這是由于混沌系統(tǒng)具有的非重復性和遍歷性，使混沌序列能夠取到當前迭代任意不穩(wěn)定周期軌道中臨近點的所有值.同時Yang等[20]指出組合優(yōu)化算法的效率受混沌序列的統(tǒng)計特性以及非線性函數(shù)全局最優(yōu)位置的影響.根據(jù)Kohli等[21]的相關(guān)研究，本文選取Chebyshev混沌映射對關(guān)鍵參數(shù)向量a進行混沌非線性調(diào)節(jié).提出改進的混沌非線性調(diào)節(jié)策略如下：

a(t)=(2-2t/tmax)·|zt|，

(18)

zt+1=cos(k·cos-1(zt)),zt∈[-1 1]

.

(19)

式中，zt為第t次迭代Chebyshev混沌映射的序列值，其中取系數(shù)k=4,初值z0=0.45，tmax=1 000時Chebyshev混沌映射及改進后參數(shù)a的迭代曲線如圖5所示.

圖5 Chebyshev混沌序列及a迭代曲線

2.5 CGWO算法流程

本文提出的混沌灰狼算法(CGWO)的具體實現(xiàn)步驟:

1) 設(shè)置種群的規(guī)模為N，維數(shù)為d，最大迭代次數(shù)tmax，通過式(10)～式(12)計算參數(shù)向量A，C及a；

2) 通過式(17)利用改進Tent映射序列產(chǎn)生種群初始個體Xi(i=1,2,…,N)，并計算所有個體的適應(yīng)度值fi(i=1,2,…,N);

3) 將適應(yīng)度值按照優(yōu)劣程度進行排序，將排名前三的灰狼個體依次設(shè)置為α狼、β狼、δ狼，其余個體均為ω狼，位置向量分別為Xα,Xβ,Xδ,Xω;

4) 根據(jù)式(18)計算混沌非線性變化參數(shù)a，然后根據(jù)式(10)和式(11)更新A和C;

5) 根據(jù)式(13)計算ω狼與α狼、β狼、δ狼的距離向量，然后根據(jù)式(14)計算狼的移動方向向量，最后利用式(15)對群體內(nèi)所有個體的位置向量進行更新.并重新計算適應(yīng)度值，更新α狼、β狼、δ狼;

6) 判斷t是否達到tmax值，如果達到則輸出最優(yōu)解，即α狼的適應(yīng)度值，否則返回3)繼續(xù)執(zhí)行，直至滿足規(guī)定的迭代次數(shù)或計算精度在設(shè)置的誤差范圍以內(nèi).

2.6 算法驗證

在MATLAB仿真平臺選取13個基準函數(shù)測試CGWO算法的性能，并與基本灰狼優(yōu)化(GWO)算法、粒子群優(yōu)化(PSO)算法、正弦余弦算法(SCA)、鯨魚優(yōu)化算法(WOA)的仿真結(jié)果進行對比分析.其中，F(xiàn)1～F7是連續(xù)單峰函數(shù)，F(xiàn)8～F13是連續(xù)多峰函數(shù)，分別對應(yīng)文獻[10]中的F1～F13.此次實驗，為了避免隨機性和偶然性對實驗結(jié)果的影響，每個基準函數(shù)分別采用不同的算法在30維獨立運行30次.五種算法的參數(shù)設(shè)置：種群數(shù)量為30，最大迭代次數(shù)為500，PSO算法的最大慣性質(zhì)量、最小慣性質(zhì)量、學習因子1和學習因子2分別為0.9，0.2，2和2.統(tǒng)計每次實驗各函數(shù)的適應(yīng)度平均值和標準差，具體實驗結(jié)果見表2.由于篇幅原因，僅展示典型多峰函數(shù)F9五種算法完整的仿真實驗過程及各算法隨機的一次收斂過程.

表2 五種算法基準測試函數(shù)尋優(yōu)結(jié)果對比

由表2測試結(jié)果可知，CGWO算法在理論最小值為0的單峰函數(shù)和多峰函數(shù)(F8函數(shù)除外)[10]中的收斂效率和求解精度明顯優(yōu)于其他四種算法.13個基準函數(shù)中，WOA算法僅在F8函數(shù)搜索到了最優(yōu)解，PSO算法在2個函數(shù)中具有較好的尋優(yōu)結(jié)果，而CGWO算法有10個函數(shù)搜索到了優(yōu)于其他四種算法的最優(yōu)解，并且其尋優(yōu)結(jié)果的標準差最小，表明該算法具有較好的尋優(yōu)性能和魯棒性.在F1～F4，F(xiàn)9和F11共6個基準函數(shù)中，CGWO算法能夠準確收斂到理論最優(yōu)值0，且尋優(yōu)成功率為100%，具有相當?shù)姆€(wěn)健性.由圖6和圖7各算法的迭代曲線可以清楚地看到，CGWO具有較好的收斂速度和全局尋優(yōu)能力，當其他算法在局部最優(yōu)解徘徊時該算法仍能夠有效地進行全局尋優(yōu).對比基本GWO算法，由表1知CGWO和基本GWO具有相同的設(shè)置參數(shù)，但不管是應(yīng)對單峰函數(shù)還是多峰函數(shù)，CGWO算法的尋優(yōu)精度和收斂速度都有顯著的提升，更能勝任復雜問題的尋優(yōu)工作.

圖6 F9各算法迭代曲線

圖7 F9各算法完整仿真實驗過程

3 實例分析

基于上述優(yōu)化模型，選取澳大利亞計算機應(yīng)用協(xié)會(ACADS)邊坡考核題[22]驗證CGWO算法的計算效率和準確性，并與基本GWO算法、PSO算法以及SCA的尋優(yōu)結(jié)果進行對比分析.各算法的參數(shù)設(shè)置與前文仿真實驗相同，分別獨立運行30次，取歷次結(jié)果的最小值作為該算法搜索得到的邊坡安全系數(shù)，考核題邊坡各土層材料見表3，尋優(yōu)結(jié)果見表4.CGWO動態(tài)尋優(yōu)過程(僅表示出歷代狼所搜索的滑動面位置)如圖8所示.

表3 考核題邊坡材料參數(shù)[22]

表4 推薦答案及各算法安全系數(shù)結(jié)果

圖8 CGWO算法動態(tài)尋優(yōu)過程

本算例中，邊界控制點Xa=20.0 m，Xb=70.0 m，Xa,b=45.0 m,則搜索區(qū)域設(shè)置如下：XL=20.0～45.0 m,XR=45.0～70.0 m，X0=XL-XR.圖8形象展示了CGWO算法搜索全局臨界滑裂面(最小安全系數(shù))的動態(tài)過程.算法通過改變設(shè)計變量能夠有效模擬各種形式的圓弧滑裂面，快速尋找并過濾可能的局部臨界滑裂面，最終確定全局臨界滑裂面如圖8紅線所示，其對應(yīng)的安全系數(shù)為1.389，略小于推薦答案[22]給出的1.390，表明算法具有較好的求解精度，能夠可靠地搜索到考核題的全局最小安全系數(shù).

與此同時，PSO,GWO,WOA及SCA的計算結(jié)果分別為1.405，1.406，1.408和1.411，均大于推薦答案.圖9直觀地給出了各算法不同安全系數(shù)所對應(yīng)的臨界滑裂面.可以看到，安全系數(shù)顯著影響著圓弧滑裂面的位置和幾何特征，在其值相差不大的情況下，滑裂面存在著一定的差異，左、右出滑點最大相差近1 m，從另一方面說明了精準求解邊坡安全系數(shù)對確定邊坡臨界滑裂面的必要性.

圖9 各算法所求邊坡臨界滑裂面

橫向?qū)Ρ雀魉惴ǖ牡屎褪諗克俣?由圖10算法安全系數(shù)的收斂過程可知，在迭代30次左右時，各算法均在局部較優(yōu)解1.41附近受到了一定的影響，此時除CGWO以外的四種算法均不同程度地陷入了局部最優(yōu)陷阱.以GWO算法為例，在迭代次數(shù)為67，254以及349時，算法均在局部最優(yōu)解徘徊搜索并停滯較長時間，最終在其附近收斂至確定值.而CGWO算法通過引入混沌擾動策略和參數(shù)非線性調(diào)節(jié)機制，能夠有效地跳出局部最優(yōu)解，持續(xù)進行全局開發(fā).以上結(jié)果表明，CGWO算法具有較高的計算精度和收斂速度，能夠準確地搜索到邊坡全局最小安全系數(shù).

圖10 各算法安全系數(shù)的收斂過程

進一步采用強度折減法對該算例進行穩(wěn)定性分析，建立如圖11所示的邊坡有限元模型.該模型采用6節(jié)點三角形單元，節(jié)點總數(shù)為662個，共劃分為1 261個單元，各土層采用Mohr-Coulomb材料模型，詳細的材料參數(shù)設(shè)置見表3，各土層剪脹角取值為0°.

圖11 邊坡模型網(wǎng)格劃分(單位：m)

本文以邊坡巖土體產(chǎn)生足夠大的塑性應(yīng)變并導致計算結(jié)果不收斂作為邊坡失穩(wěn)的判據(jù).當?shù)降?1次時，在坡腳處產(chǎn)生了超過設(shè)定閾值的塑性應(yīng)變，此時邊坡的安全系數(shù)為1.380，所對應(yīng)的邊坡累計位移云圖及等效塑性應(yīng)變云圖如圖12和圖13所示.該方法所求安全系數(shù)與本文方法相比稍小，分析可能是由于模型缺少設(shè)置各土層材料剪脹角屬性，導致計算得到的安全系數(shù)略小.當剪脹角參數(shù)設(shè)置為10°時，所得邊坡的安全系數(shù)為1.402，與極限平衡法的計算結(jié)果較為吻合.

圖12 強度折減累計位移云圖

圖13 等效塑性應(yīng)變云圖

同時可以看到，由于計算結(jié)果不收斂導致僅在邊坡坡腳處產(chǎn)生了明顯的塑性集中區(qū)域.并且當土層分層情況相對復雜時，邊坡內(nèi)部并未產(chǎn)生明顯的塑性貫通區(qū)域，不易直觀地判定邊坡臨界滑裂面的所在位置.本文方法所得的臨界滑裂面如圖13紅色虛線所示，與有限元強度折減法所求出的塑性發(fā)展區(qū)域契合度較高，再次驗證了本文方法的有效性.同時，與建模及材料設(shè)置過程相對復雜的有限元強度折減法相比，本文所提出的方法得益于CGWO算法出色的全局尋優(yōu)能力，僅通過設(shè)置邊界控制點即可實現(xiàn)搜索區(qū)域的全邊坡范圍覆蓋，在初始應(yīng)力未知的情況下就可以快速搜索得到邊坡的最小安全系數(shù)及其對應(yīng)的臨界滑裂面，具有一定的工程指導和推廣意義.

4 結(jié) 論

1) 基于混沌理論，從算法的群初始化、擾動策略、關(guān)鍵參數(shù)調(diào)節(jié)三個方面對基本GWO算法進行改進.仿真實驗結(jié)果表明，所提出的CGWO算法相比基本GWO算法，具有較強的尋優(yōu)能力和魯棒性，更適合處理復雜多峰函數(shù)的尋優(yōu)問題.

2) 基于簡化Bishop法，采用邊坡左、右出滑點及圓弧滑裂面圓心的橫坐標構(gòu)建安全系數(shù)的優(yōu)化模型，該方法具有操作簡單、變量范圍易于設(shè)置等優(yōu)點.選取ACADS邊坡考核題進行計算，CGWO算法具有較高的收斂效率，能夠避免局部最優(yōu)解的干擾，可靠地搜索到全局最小安全系數(shù).

3) 所提出的CGWO算法適用于確定二維圓弧滑裂面巖土邊坡的安全系數(shù)及其對應(yīng)的臨界滑裂面，同時對解決其他低高維優(yōu)化問題具有一定的借鑒意義.與建模過程較為復雜的有限元強度折減法相比，該方法能夠在初始應(yīng)力未知及缺少建模參數(shù)的情況下對邊坡的安全狀態(tài)進行準確評估，具有一定的工程實用價值.