基于FOA法搜索危險滑動面的三維極限平衡分析

常晨瑩，林 川，袁占軍，趙蘭浩

(1.河海大學 水利水電學院，江蘇 南京 210098；2.黃河水利委員會供水局，河南 鄭州 450003)

1 研究背景

在邊坡穩定分析領域，剛體極限平衡法是理論最成熟的方法之一，已經成功應用于許多工程實踐[1-3]之中。陳祖煜等[4]在二維Spencer法的基礎上進行了三維擴展，提出了一個理論基礎更為嚴密、計算步驟相對簡單、收斂性能更好的三維極限平衡分析法。極限平衡法對邊坡進行穩定分析時，只能對已確定的滑動面進行計算，這就需要先找到最不利滑動面，再進行計算。傳統的搜索滑動面的方法有二分法、坐標輪換法等，這些方法對于復雜邊坡，搜索極易陷入局部極小值，很難搜索出真正的最不利滑動面。隨著科技的發展，一大批智能優化算法隨之涌現，主要有遺傳算法、模擬退火法和仿生算法。謝桂華等[5]利用遺傳算法建立邊坡可靠度分析模型，該方法對目標函數要求比較寬松，可以得到全局最優解，但是遺傳算法對實際問題的搜索空間大小的變化適應能力較差，且計算量大，容易出現早熟現象；李守巨等[6]基于模擬退火法提出了一種搜索邊坡臨界滑動面的數值方法，該法極大地提高了計算速度，并可以保證全局收斂，但模擬退火法在搜索滑動面時，參數的選擇對計算的結果影響很大。仿生算法[7-8]是借助仿生智能思想和計算機技術而設計的智能算法。其中，果蠅優化算法[9]是一種過程簡單、全局收斂性高的仿生優化算法。一些學者[10-12]已成功將其應用到邊坡穩定分析中。吳成濱等[13]在用果蠅優化算法對彈性模量進行反演計算時，在算法中添加觸壁機制，當果蠅種群越界后就會自動返回，并且對果蠅種群的初始位置進行限定，彌補了算法后期收斂速度慢的缺陷。

本文在三維極限平衡法基礎上增加了地震慣性力和加固措施力的作用，采用改進的果蠅算法不斷地搜尋可能的滑動面，再通過改進的三維極限平衡法計算危巖體擬定滑動面的安全系數，最終得到危巖體最危險滑動面的穩定安全系數。

2 計算方法

2.1 改進三維極限平衡法

原有的三維極限平衡法只能用來分析天然狀態下的邊坡，本文對其進行改進，在力的平衡方程和力矩平衡方程中增添地震慣性力和加固措施作用力，使其既能進行危巖體動力穩定分析，又能進行加固措施貢獻度分析。

危巖體穩定性分析用的三維極限平衡法模型是建立在三維極限平衡分析方法[14]基礎之上的。將危巖體劃分成若干條塊，其中某一條塊受力模型見圖1。x和z的正方向分別為滑動方向和重力的反方向，y軸的正方向按右手法則確定。列界面法向應力Q，平行于y軸；行界面條間力E，平行于xoz面，且與x軸夾角為β；底滑面法向力為Ni；底滑面剪切力Ti，與xoz面的夾角為ρ。

根據三維極限平衡分析法計算之前要建立S-S′坐標系，見圖2。S軸與條柱體行界面上側向力E方向一致，與水平面夾角為β。S軸與S′軸垂直。在S-S′坐標系中，Ni′為滑動面法向力Ni在平面XOZ里分量，Ti′為切向力Ti在滑動面內平行于平面XOZ方向上的分量，與水平面夾角為αx。

圖1條柱受力模型

圖2條柱力在S方向上的投影

建立S′方向的靜力平衡方程:

(Wi+Fei+Ni+Ti+Ri)S′=0

(1)

摩爾-庫倫準則:

(2)

連理(1)、(2)兩式可得到Ni和Ti，其中Ni為:

Ni=

(3)

式中:tanφe=tanφ/F,ce=c/F。

S方向的整體靜力平衡方程式:

∑(Wi+Fei+Ni+Ti+Ri)S=0

(4)

y方向的整體靜力平衡方程式:

(Ni+Ti+Ri)y=0

(5)

繞y軸的整體力矩平衡方程式:

∑(Wi×C+Fei×C+Ni×P+Ti×P+

Ri×Q)y=0

(6)

式中:F為抗滑穩定安全系數；Ti為滑動底面法向力，其方向導數是關于ρ的函數；Ni是滑動底面法向力，其方向導數可由計算程序導出為已知量；Wi為條柱體的重力，作用在條塊體的質心；Fei為條柱體受到的地震作用力，作用在條塊體的質心，本文采用擬靜力法[15]考慮地震作用；Ri為危巖體所受的加固措施力；C為條塊的質心坐標；P為條塊滑動面的形心坐標；Q為加固措施作用力的作用點坐標。

平衡方程中共有3個未知數，分別為F、β、ρ，利用牛頓-勒普生法來迭代求解。

2.2 果蠅優化算法尋找最不利滑動面過程

果蠅優化算法的思路是果蠅利用嗅覺搜尋食物，飛近食物后利用視覺飛向食物。它具有全局優化的能力且算法概念非常簡單，已被成功地應用于多個領域。

當滑坡的入滑點A與出滑點B坐標為已知時，只要確定出滑點的B處切線的斜率便可確定一條圓弧滑動線，進一步由若干條滑動線確定整個滑動面。根據改進的果蠅算法，本文通過果蠅種群搜索出滑點處切線的垂線的斜率K值來確定滑動面，并對搜索范圍進行合理限制。如圖3所示，每條圓弧滑動線的圓心都在入、出滑點A、B連線的垂直平分線l1上，K值下界為直線l1的斜率。K值的上界要滿足出滑的運動容許條件，即圓弧滑動面在出滑點的切線斜率要大于出滑邊界的切線斜率，如圖4所示，滑動線應像圖4(a)一樣順滑，而避免出現如圖4(b)一樣的折拐現象，因此K值上界為出滑邊垂線l2的斜率。

(AB為圓弧滑動面，BC為滑坡的出滑邊界)

圖3K值上、下界示意圖

圖4邊坡滑動線示意圖

將果蠅優化算法尋找最不利滑動面過程以程序實現，程序流程圖見圖5，程序實現步驟如下：

步驟1：設置初始果蠅種群的數量以及計算時的迭代次數，并對果蠅位置(即K值)限制范圍。

步驟2：在搜索范圍內隨機確定每個果蠅種群的初始位置，提供滑動面給改進三維極限平衡法的程序，計算得到該種群所在位置的抗滑穩定安全系數值，并記錄最優種群的位置。

步驟3：所有果蠅種群飛到的最優位置，再從該位置向周圍繼續搜尋，并計算味道濃度(即安全系數值)。

步驟4：對新位置的味道濃度值與原來的最優值進行比較，如果前者優于后者，則新位置替換最優位置，反之最優位置不變。

步驟5：判斷是否達到最大的迭代步數，是則循環結束，否則回到步驟3繼續循環。

圖5計算流程圖

3 工程實例

3.1 工程概況

索風營水電站位于烏江干流六廣河河段，Dr2危巖體在水電站進水口右上方的T1m灰巖陡壁上,最低處高程870 m，最高處高程1 085 m。Dr2危巖體內夾層較發育，對危巖體穩定影響較大的夾層有4條，分別為發育于T1y3底部的J1、中部的J2、發育于T1y3頂部，與T1m1接觸的J3、發育于T1m1下部，出露于危巖體陡崖腳附近的J4。危巖體主要裂縫有3條，分別為位于危巖體后緣的L1；位于L1外側的L2；位于L2外側的L3，在高程980 m處露出陡壁外，形成危巖體外側倒三角體。

高程1 070 m以上的部分為危巖體頂部，是由灰巖T1m形成的緩坡平臺，地表有殘破積黃色黏土夾碎石層分布，坡度5°～10°；高程960～1 070 m的陡壁為中部危巖體，主要由T1m灰巖構成，坡度大于70°，局部呈倒懸狀；底部是高程960 m以下部分，外側為早期上部傾斜滑倒下來的崩塌堆積物(III號堆積體)形成的斜坡，坡度24°～40°。

總體來看，危巖體呈單斜構造，下寬上窄，平面長軸呈南北向，向西凸起，東側由弧狀拉裂縫L1構成后緣邊界，底部止于T1y3泥巖。危巖體剖面如圖6所示。

圖6危巖體剖面圖

水電站所處位置的地震基本烈度為6度，雖然地震烈度不高，但危巖體的安全對電站運行有極其重要的影響，電站運行期間該危巖體一旦失穩將影響整個大壩的安全和正常運行。為防止危巖體上部發生傾倒變形，主要以“拉”的方式對上部的進行加固；中部和下部主要采用阻滑和抗剪的措施來防止可能發生的座滑失穩。概括起來為“頂戴帽，腰系帶，腳穿鞋”。

“頂戴帽”主要是Dr2危巖體頂部裂縫外部設置截水溝，并在裂縫處設置地梁，在混凝土地梁中設47根1 000 kN有粘結預應力錨索；“腰系帶”主要是在L1外側1 040～1 060 m高程崖面采用176根傾角為15°的2 000 kN無粘結錨索進行加固處理；“腳穿鞋”主要是在危巖體下部T1y3泥巖內布置一排6@φ7.0 m的深層抗滑樁，930 m高程施工平臺布置7條錨固洞，并設置傾角為15°的2 000 kN有粘結預應力錨索20根，形成洞、樁、錨聯合加固；并在危巖體倒懸體采用貼坡混凝土進行支撐。

3.2 模型及計算參數

考慮到危巖體邊坡特點且保證計算效果，將其模型鉛直剖分成30個塊體(如圖7所示)，共5行6列。其中6列對應7條滑動邊界，每條滑動邊界由出滑點處切線的斜率K控制，入、出滑點分別為泥巖與危巖體和泥巖與堆積體接觸的邊界，見圖8。假定滑動邊界為圓弧。每個滑動面是由兩側的滑動邊界、入滑邊界和出滑邊界組成的曲面。Dr2危巖體計算參數見表1。

表1 危巖體計算參數表

條柱體受力情況要考慮地震作用，水平地震慣性力Fei的方向與危巖體滑動方向相同。考慮加固措施的作用時，錨索錨固力、錨固洞抗滑力和抗滑樁抗滑力的作用分解到加固措施所在的每個塊體上面，包括力的大小、方向和作用點，錨固洞和抗滑樁按鋼筋混凝土抗剪強度進行極限承載力計算，加固措施作用力見表2。泥巖灌漿加固通過修改灌漿后泥巖材料的c、φ值來實現。

調用三維極限平衡法進行穩定分析時，設置F、β、ρ初值分別為1.3、28.66°、0°，最大迭代次數為1 000。用果蠅優化算法搜索滑動面時，令果蠅種群數量為5，最大迭代次數為10；并對危巖體7條滑動邊界出滑點處的7個K值的搜索范圍進行限制，上界為出滑邊界垂線斜率，下界為入、出滑兩點連線的垂直平分線的斜率，其取值見表3，計算中為了降低計算量，將K1～K2、K3～K6、K7用3個K值等效，取包含K界限的較小區間。

4 計算結果與分析

本工程為乙類建筑，根據《水利水電工程邊坡設計規范》的規定，Dr2危巖體等級為A類II級，在短暫及偶然設計工況下抗滑穩定安全系數為F=1.10。

根據危巖體的加固情況，為反映單種加固措施的抗震安全穩定貢獻度，分6次單獨添加單個加固措施對危巖體進行抗震復核，計算結果見表4。

圖7危巖體計算模型

以添加了所有加固措施后的危巖體的穩定分析為例，危巖體最危險滑動面安全系數的搜尋過程如圖9、10所示。歷史最優值在迭代后期趨于穩定，最終得到最不利穩定安全系數為1.1138。危巖體最危險滑動面見圖11。

從表4可看出：

(1)地震狀態下，天然危巖體沿最危險滑動面的抗滑穩定安全系數為1.0124，略高于1，危巖體處于臨界破壞狀態，與實際工程情況相符，不滿足規范要求。而采取了一系列加固措施加固后，危巖體最不利滑動面的抗滑穩定安全系數提高到了1.1138，達到預期目標，滿足規范要求。

圖8危巖體出、入滑點示意圖

表2 加固措施合力作用表

表3 K值取值范圍表

表4 地震狀態中各種工況下危巖體的抗滑穩定安全系數及滑動參數

圖9果蠅種群尋優過程

圖10歷史最優值變化圖

圖11最危險滑動面圖

(2)地震工況下，單種加固措施安全穩定貢獻度由大到小排序為：泥巖灌漿加固>上部錨索>錨固洞>抗滑樁>頂部錨索。頂部錨索設置在危巖體最頂部，主要用于防治頂部危巖體傾倒破壞，對整體的穩定的貢獻幾乎為零。泥巖灌漿加固增加了巖體抗剪斷摩擦力，再聯同上部錨固索的錨固作用，危巖體與L1裂縫后的T1m灰巖貼合得更加緊密，限制危巖體上部的傾倒變形。下部930 m高程聯系大梁設置20根2 000 kN錨索，與錨固洞和抗滑樁形成聯合加固，既為危巖體整體增大了阻滑力，又增強了危巖體底部基座的剛度，極大地限制了底部危巖體的坐滑變形。

(3)在不同的加固措施作用下，最危險滑動面位置變化幅度不大，比較穩定，其中K3～K6一直穩定于4.27，K7穩定在3.24左右，兩者分別為K6和K7的上界。由此可見出滑點處最不穩定，危巖體最危險滑動面在Ⅲ號堆積體與泥巖接觸邊界。

5 結 論

本文以索風營Dr2危巖體為背景，采用果蠅優化算法搜索危巖體最不利滑動面，并以三維極限平衡法對危巖體進行了地震工況下的安全穩定分析，具體結論如下：

(1)對三維極限平衡法的平衡方程進行改進，增添地震慣性力和加固措施作用力，使方法不僅可以進行動力分析，還可以分析加固措施對危巖體穩定的貢獻度。

(2)建立了以改進果蠅智能算法為基本框架的三維極限平衡法，改善了傳統極限平衡法無法所搜滑動面的缺陷，從而降低人為因素對參數選取的干擾，提高了邊坡失穩的預測精度。

(3)加固前的危巖體在地震狀態下處于危險臨界狀態，加固后的危巖體抗滑穩定安全系數為1.1138，滿足規范要求，可以維持穩定。