基于改進狼群算法的邊坡極限平衡穩定性分析方法研究

馬金城 石 崇 張一平 張文浩 王樂榮

(1.河海大學 巖土力學與堤壩工程教育部重點實驗室,南京 210098;2.河海大學 巖土工程科學研究所,南京 210098)

邊坡穩定性分析是巖土工程十分重要的研究內容,涉及水利水電、公路、鐵路、礦山工程等諸多領域,是判斷邊坡穩定性,是否需要加固及采取何種防護措施的主要依據.在進行邊坡穩定分析時,關鍵是如何得出滑面并利用滑動力、抗滑力計算安全系數、優化得出最危險滑面.因此,對于危險滑面搜索方法進行理論和應用研究非常有價值.

在常規的條分法邊坡穩定計算中,如瑞典條分法和畢肖普法假定滑坡體為不變形的剛體,利用自重計算滑面上滑動力與抗滑力,再利用抗滑力矩與滑動力矩之比來定義穩定性系數[1].而自然巖土體處于彈塑性變形狀態,達到屈服時應力可向周圍傳遞,因此條分法得到的各條分底部抗滑力、滑動力往往與實際偏差很大.而數值模擬方法綜合考慮物理方程、平衡方程、幾何邊界條件,因此可以逼真地得到巖土體應力場,但不能直接給出穩定性系數.如果將這兩者有效結合起來,利用數值模擬得到真實的應力場,采用矢量和法可定義穩定性系數,對邊坡穩定性計算有很好參考價值.因此近些年很多學者致力于探索危險滑面搜索方法的工作,如遺傳算法[2],蟻群算法[3],粒子群算法[4],模擬退火法[5]等.這些算法在邊坡滑面搜索中得以運用,與早期使用的窮舉法、二分法相比,上述方法的精度和收斂性提高了很多,對邊坡穩定研究和實際工程應用起到了積極的作用,但蟻群算法容易出現停滯等現象[6],粒子群算法會過早陷入局部最優等缺點[7],模擬退火算法和遺傳算法對參數選擇要求較高,如若參數選擇不當會使搜索陷入局部極值或搜索結果不令人滿意[8].

狼群算法[9]是2011年提出的一種新型的智能優化算法,該算法在全局搜索和局部開發能力上具有良好的適用性,被廣泛地運用到路徑規劃等實際問題中,具有魯棒性好,全局尋優能力強,收斂速度快等優點,但目前尚未運用到邊坡穩定性分析之中.本文基于傳統的狼群算法,嘗試結合邊坡穩定性計算特點進行優化改進,建立狼群滑面位置搜索及穩定計算模型,并用已有成果驗證其可行性和準確性.

1 極限平衡穩定性分析基本原理

1.1 抗滑力和下滑力的計算

對一個二維邊坡中潛在的滑面,計算穩定時一般如圖1所示將滑面劃分為若干段,不同段相交點稱為一個狀態點,從起始點到終止點通過狀態點首尾相連可形成一條路徑,每條路徑都是一條潛在的滑動面.

圖1 滑面受力示意圖及單元體圖

設滑動方向向右,逆時針為下滑力正方向,則其中某一段上的抗剪強度τf和切向應力τ的計算公式如下：

第i段的平均滑動力：

第i段的平均抗滑力：

式中：σx、σy、τxy分別為單元形心沿x,y方向的正應力和切應力；φ和c分別為材料的內摩擦角和內聚力；σα為斜面的法向應力；α為斜面與水平面的夾角.

1.2 安全系數的定義

利用數值模擬的方法獲取邊坡真實應力場后,計算出每一段上的滑動力和抗滑力,然后采用矢量和法定義安全系數,安全系數定義為：

式中：A和B分別為邊坡起始點和終止點；n為滑動面所經過的段的數量；τ為沿滑動面滑動方向的剪應力；dL為沿滑動面每一段的長度增量；τi為從第i-1個節點到第i個節點路徑上的剪應力；τfi為從第i-1段到第i段路徑上的抗剪強度；Δd i為從第i-1階段到第i階段路徑的長度[10].

2 滑面搜索的改進狼群算法

對于已知滑面邊坡,采用上述極限平衡穩定性分析方法可直接求出安全系數,但對于滑面未知的邊坡,如何確定滑動面上點的位置,找到最危險滑面是非常困難的[11].從數學問題看,確定滑面實際上是圖形路徑的優化問題.近年來,狼群算法作為一種生物仿生算法被提出,可用于圖形最優化問題的求解.因此,本文針對二維邊坡的網格,將滑出滑入點當做頭狼與母狼的位置,某點的應力當作當前位置的獵物氣味濃度,網格節點當作猛狼,網格內部的某些點為探狼,通過不同種類狼之間的信息交互、“勝者為王”的頭狼產生機制和“強者生存”的狼群更新機制來尋求全局路徑的最優解,可以用于解決邊坡最危險滑面的搜索.

2.1 經典狼群算法基本原理

經典狼群算法分為如下幾個步驟：

1)頭狼生成：在某一待尋優區域內,將位置最優的狼視為頭狼.頭狼是狼群的核心,負責狼群的指揮工作[12].

2)游走行為：選取除頭狼外最佳的數匹狼作為探狼,執行游走行為,根據某種計算規則得到探狼處氣味濃度,然后和當前頭狼氣味濃度進行比較,選擇二者中較優者為頭狼.

3)召喚行為：當頭狼發出召喚行為后,猛狼以較大步長向頭狼位置逼近.如果在奔襲過程中,猛狼位置優于頭狼,則猛狼代替頭狼發出召喚行為,否則猛狼繼續奔襲,直至距離頭狼位置小于圍攻距離,此時轉化為圍攻行為[13].

4)圍攻行為：當猛狼奔襲到距頭狼較近位置時,猛狼聯合探狼對獵物進行圍攻.

5)“強者生存”的狼群更新機制：剔除掉較差的狼,只有適應性較好的狼才能生存下來.

2.2 狼群算法的改進

1)模型的建立：根據邊坡形狀及地層分布狀況,結合邊界條件,對邊坡進行合理的有限單元網格劃分,單元格數量由上述因素及所需結果精度共同確定.得到有限單元網格后,計算邊坡的應力場,得到每個節點處的應力和坐標[14].將所得應力場視為草原平面,也即狼群活動的平面,將節點處的應力值視為當前位置的獵物氣味濃度.

2)頭狼生成：在狼群中,頭狼和母狼處于支配地位,頭狼負責狼群的指揮,母狼負責配合頭狼工作.頭狼和母狼對整個狼群的優劣起決定作用.如圖2 所示,頭狼和母狼在可能滑入點與滑出點的某一區域內,初始頭狼與母狼分別為該區域內橫坐標值最小的點.算法中引入母狼的原因是為配合頭狼工作,兩者共同把握整個狼群的位置,避免陷入局部最優.

圖2 可能圓心位置與滑出滑入點切線斜率

3)狼群區域劃分：如圖3 所示,以滑入點(X0,Y0)為坐標原點建立搜索區域坐標系,用相同間距為x d=(X n-X0)/n的縱向線將待尋優區域劃分成n部分,其中X n為滑出點橫坐標.將每條縱向線鄰域視為一個搜索區域,分別記為LQ1,LQ2,…,LQn.將各縱向線橫坐標依次記為X1=(X n-X0)/n,X2=2(X n-X0)/n,…,X n=X n.在給定滑出滑入點位置后,需確定下一個點(X ML1,Y ML1)的搜索位置,潘家錚等人通過研究表明,對于簡單均質土坡,作邊坡中點的垂線和法線,以坡面中點為圓心,以1/4坡長和5/4坡長為半徑畫同心圓,最危險滑弧的圓心在如圖2所示陰影部分ABDC范圍內.據經驗可知,圓心縱坐標越小,滑出點和滑入點與圓弧相切切線斜率越大(滑出點和滑入點切線的斜率僅用于確定滑出滑入點的下一個點搜索范圍,并不對滑動面是否為圓弧面做限制).設滑出點處的切線最大斜率絕對值為k1,滑入點處最大斜率絕對值為k2.

圖3 狼群搜索區域示意圖

如圖4(a)所示,LQ1的搜索范圍為：

如圖4(b)所示,在LQ i中,橫坐標在i(X n-X0)/n的左右鄰域內波動.

圖4 狼群搜索區域及探索步長

由于滑動面點為一下凹曲線,因此點(X0,Y0)與點(X ML1,Y ML1)所在直線為LQ2搜索區域的下限,X=X1-δ,X=X1+δ是X=X1左右鄰域.點(X0,Y0)與點(X n,Y n)所在直線為LQ2搜索區域上限.也即LQ2的搜索范圍為：

根據上述規則,由前一搜索區域所得猛狼位置來推算下一狼群搜索區域,直至最后一匹猛狼產生.這樣做可以有效縮小搜索區域的范圍,提高搜尋效率,減少搜索時間.

4)游走行為：探狼從每個待尋優區域最上側的點開始以固定步長y d向下游走探索.處于頭狼和母狼附近的探狼分別以固定步長x d水平向右游走探索.

5)召喚行為：當頭狼和母狼及探狼已處于各自狩獵位置后,頭狼發出召喚行為,召喚探狼周圍某一范圍內猛狼.視探狼所處單元格內部為召喚區域.處于單元格節點處的猛狼在接收到召喚指令后,向附近區域探狼奔襲,通過召喚區域上猛狼處氣味濃度綜合確定探狼處氣味濃度,也即由已知多邊形節點的坐標和應力通過機器學習的方法確定待求點的應力狀態(具體確定過程見(8)).

6)圍攻行為：可能最優路徑上的頭狼和猛狼經上述方式產生后,對獵物展開圍攻,形成一狼群.

7)“強者生存”“優勝劣汰”的狼群更新機制：將適應度D作為目標函數,則狼群算法可將求最小安全系數的問題轉換為求D最大值問題,目標函數定義為：

適應度D為狼群對草原適應性的好壞,D越大,則狼群越優.在完成每次可能最優路徑適應度D i計算后,將其適應度D i與當前適應度最好的Dlead進行比較,直至最后一條路徑計算結束,將適應度最好的狼群視為優勝狼群,也即：

其中：Dmax為全局最優解,其對應的各種路徑為最適宜路徑.

8)機器學習方法確定猛狼處氣味濃度：在Python 中使用sklearn庫中的Gradient Boosting Re-gressor模塊,將給定數據(結點坐標和應力值)分為訓練集和測試集[15](將部分已知節點應力和坐標作為訓練集,其余節點坐標和應力作為測試集,在訓練集中找尋規律,在測試集中驗證規律[16].可采用L1范數損失函數(loss function)來估量模型的預測值f(x)與真實值Y的不一致程度,其值的大小可以用來評價預測的準確性,該函數把目標值(Y)與估計值(f(x))的絕對差值的總和最小化[17]：

相比于傳統的插值法求中間部分應力,機器學習方法可以很直觀地通過S值的大小來判斷對某點預測的準確性.

3 程序實現與算例驗證

根據以上原理,形成改進狼群算法邊坡滑面搜索方法,計算流程圖如圖5所示.

圖5 改進狼群算法流程圖

本文采用python語言編制程序,引入一經典算例來驗證狼群算法在邊坡滑面搜索以及計算安全系數的準確性.該算例為澳大利亞計算機協會于1987年委托BDonald和PGiam 設計ACADS邊坡穩定性分析程序的考核題,該題的參考滑動面如圖6(a)所示,推薦答案的安全系數為1.00[18],具體參數為：容重γ=20 k N/m3,泊松比μ=0.25,黏聚力c=3.0 kPa,內摩擦角φ=19.6°.由于該算例已經過大量學者的引用和對比分析驗證,因此其結果視為標準值.

基于改進狼群算法理論的邊坡臨界滑動面搜索方法搜索出最危險滑動面如圖6(b)所示,對應的安全系數為1.03675,計算所得安全系數與推薦值1.00誤差約3.6%,采用本文方法得到的臨界滑動面與推薦滑動面基本一致,這表明基于改進狼群算法的邊坡穩定性分析是可行的,其值偏高的原因是數值應力結果比剛體極限平衡協調,因此有更高的承載能力.

圖6 邊坡滑面(單位：m)

4 工程案例研究

溪洛渡水電站左岸谷肩堆積體的典型地質剖面如圖7(a)所示.本文采用flac-3D 軟件建立如圖7(b)所示數值模型,選用摩爾-庫侖的本構模型,以及表1中巖土物理力學參數進行數值計算,得到其應力場后,采用Python語言編寫程序,運用改進狼群算法搜索其最危險滑面,并計算其安全系數.

表1 巖土物理力學參數

如圖7(a)所示,IN13鉆孔傾斜儀與剖面較為接近,圖7(c)表明,溪洛渡水電站左岸谷肩堆積體的最危險滑面位于滑坡堆積體與宣威組上部交界處.這與歷史監測資料所揭示的滑動面位置一致.文獻[19]在自然工況下采用改進動態規劃算法,圓弧形滑面計算所得安全系數為1.100,而在相同工況下采用狼群算法得到的安全系數為1.044,如圖7(d)所示.經對比滑面可知,狼群算法搜索得到的滑面可保證更多比例滑面位于宣威組,而現場勘查表明該邊坡處于蠕滑狀態,自然工況安全系數為1.00～1.05之間,計算結果接近蠕滑狀態自然工況下安全系數的上限1.05,故采用本方法得到的結果更加接近于實際.

圖7 工程案例計算結果

5 結論

基于數值計算得到的邊坡應力場,提出一種基于改進狼群算法的邊坡臨界滑動面搜索方法,為邊坡穩定分析提供一條新的有效思路.得到主要結論如下：

1)改進狼群算法是利用數值模擬手段獲得真實的應力場后,采用矢量和法來定義安全系數.在邊坡穩定性分析過程中可得到較為準確的邊坡最危險滑動面和最小安全系數.

2)在運用改進狼群算法的過程中,用機器學習方法代替傳統的插值法來計算指定坐標點處的應力,可以通過損失函數直觀地看出計算結果的準確性.

3)通過對經典邊坡算例的分析,本文所提方法的計算結果與典型算例參考答案具有較好的一致性,證明了該算法的合理性.