含剪切節理面巖質邊坡滑裂面位置及穩定性研究

陳東宇 劉文連 眭素剛 許漢華

摘 要：同向雙平面滑動是存在單一地質斷層面（剪切節理面）巖質邊坡的常見破壞模式之一，但對該種類型的滑裂面計算方法并不充足。為了能夠更加高效準確地尋找邊坡的滑裂面位置，判斷邊坡的穩定性，基于極限平衡理論和非線性數學規劃模型，假設滑體的滑動方式為同向雙平面滑動，再假設目標函數為該巖質邊坡的安全系數，運用MTALAB全局最優搜索法，計算含剪切節理面工程邊坡在天然工況作用下的滑裂面位置及穩定性，并與極限分析法、強度折減法和畢肖普法進行對比分析。研究結果表明：基于極限平衡理論和非線性數學規劃模型得出的滑裂面位置與安全系數基本一致，驗證此類方法的可行性，為存在單一地質斷層面巖質邊坡的滑裂面計算和穩定性分析提供了新依據。

關鍵詞：巖質邊坡；滑裂面計算；安全系數；極限平衡法；非線性理論；最優化原理

中圖分類號：TU452

文獻標志碼：A

巖質邊坡的地層是分布在我國西南地區的一種工程性質較差的軟巖，從上世紀60年代國家開始對西南地區進行開發建設，巖質邊坡就慢慢地出現在工程作業中。在工程作業中，滑坡等地質災害就出現在了我們的視線當中，其不僅制約了工程建設的發展，而且還對人類的生產生活產生了許多影響［1-2］。針對這一特殊的軟質巖土，預測邊坡穩定性及滑坡失穩面的位置，對土建、水利、交通等方面有非常重要的意義。

在結構面控制型邊坡中，主要的失穩模式有平面滑動、楔型滑動、傾倒破壞等［3-4］。目前常見的邊坡穩定性預測分析方法有：極限平衡法、強度折減法、系統分析法、數值模擬等［5-8］。李芬［9］、朱擎［10］、程小龍［11］等通過強度折減法和雙強度折減法解決了雙平面破壞模式下的臨界失穩問題，討論了巖質邊坡的破裂角和內摩擦角之間的關系，進而提出安全系數的簡化估算公式。陳建宏等［12］使用了極限平衡分析的上下限法，并把模糊化處理參數和置信水平等概念加入到了平面滑動的巖質邊坡工程中，此方法能夠有效地解決平面滑動巖質邊坡的參數不確定問題。CHENG［13］、ARDESTANI［14］等采用極限平衡法，通過二維平面應力分析求解邊坡的臨界滑動面及穩定性系數。SCHLOTFELDT等［15］以反傾巖質邊坡為研究對象，把數值分析與極限平衡法相結合，提出了這類巖質邊坡的研究方法。張崇波［16］對香港秀茂坪巖質邊坡進行了可靠性分析，利用拉丁超立方抽樣可靠性分析方法，分析了同時考慮參數不確定性和最危險滑面不確定性的平面滑動巖質邊坡可靠性，優化了可靠性分析方法的計算精度和效率，將研究成果應用到工程實踐，并將該方法應用到山東臨沂換流站站址邊坡的評價當中。陳志強等［17］提出了如果巖質邊坡的滑裂面位置和產狀無法確定，就可以用巖質邊坡平面滑動滑面極限傾角來計算邊坡的最小安全系數，并考慮拉裂縫深度對滑坡穩定性的影響。高丙麗等［18］基于坐標投影法，提出了三維單滑面和雙滑面型塊體的穩定系數計算方法，并基于Matlab開發出適于巖質邊坡工程中平面多面體塊體和曲面塊體穩定性分析的CPG程序，實現了結構面、臨空面及不穩定塊體的空間表示及可視化。但對存在剪切節理面巖質邊坡的滑裂面位置的計算除了工程軟件外，并沒有其他更加簡便的方法了。

針對以上研究的不足之處，本文以兩種不同地層的巖質邊坡為研究對象，將極限平衡理論和非線性數學規劃理論結合起來，通過求解最優化方法確定安全系數從而求解巖質邊坡滑裂面，建立一種直接計算同向雙平面滑動邊坡滑裂面位置的計算方法。

1 計算原理

對于穩定性受確定性結構面控制的巖質邊坡，應根據結構面的產狀與強度參數，采用極限平衡方法計算邊坡的穩定性。根據結構面的空間展布，邊坡失穩類型分別為平面滑動和空間滑動。嚴格意義上邊坡滑體都是空間塊體，但對于單一結構面控制的滑體，或由兩個及兩個以上平面構成的畫面，只要這些平面走向大致相同、與邊坡破面走向平行或接近平行，且滑體兩側不受約束或約束不大時，即可按平面滑動進行分析，即本文算法的基本原理。如圖1所示，本文研究對象主要是存在兩種不同地層的復合巖質邊坡，該類邊坡主要是由于不同紀元產生的不同類型巖石共同組成的一種特殊的復合地層巖質邊坡，地層Ⅱ的強度較大，其次是地層Ⅰ，該巖質邊坡強度較弱的部分即最容易發生滑裂的部位是兩種地層的交界面（剪切節理面），通常的破壞形式為同向雙平面滑動。傳統計算方法常見的有有限元法、強度折減法、極限平衡法，可通過工程軟件大量的計算來實現，但計算結果通常為圓弧型滑裂面或為不規則滑動面，不符合邊坡實際滑裂面形狀。

本文算法基于塑性力學的基本理論，將邊坡安全系數作為目標函數，將底滑面和后緣滑裂面的位置坐標作為決策變量，同時考慮滑體非線性約束，如極限平衡方程，底滑面、后緣滑裂面的屈服條件，即Mohr-Coulomb屈服條件，最后使用優化算法可同時求解得到邊坡的安全系數及底滑面和后緣滑裂面控制點的位置坐標。

2 邊坡滑塊計算模型

將邊坡滑塊離散出來成為一塊狀單元如圖2所示，該塊體單元為一四面體，每條邊為一結構面，除了該單元體的自重，底滑面和后緣滑裂面還受到法向力和剪力；由于已有極限平衡法里的剛性假設，該滑體的內部不會發生變形，故可以直接對該滑體進行受力分析。該滑塊單元的受力如圖2所示，滑塊在形心處受到自身重力G，θ為兩種地層交界線（底滑面）與水平方向的夾角，β為后緣滑裂面與水平方向的夾角，由于不施加外部荷載，故臨空面AD和頂面CD是不受力的，滑塊底滑面AB受到支持力NAB和剪切力SAB，后緣滑裂面也受到支持力NBC和剪切力SBC，此時該滑塊單元處于極限平衡狀態。

3 巖質邊坡的力學機制分析

3.1 滑塊的力學模型及幾何關系

根據含剪切節理面的巖質邊坡變形破壞機制，可將發生破裂的面分為底滑面AB和后緣滑裂面BC兩個部分，如圖1所示，滑裂面、臨空面和頂面共同組成滑塊。

假設滑塊為均質剛體，則滑塊的非線性約束重力方程為

式中：G為滑體的自重；ρ為第四系地層的密度；g為重力加速度；（xA，yA）、（xB，yB）、（xC，yC）、（xD，yD）分別為A、B、C、D 4點的坐標。

3.2 底滑面和后緣滑裂面的幾何關系

如圖1所示，以O點為原點來確定A、B、C、D各點的坐標。

（1）底滑面AB的坐標及長度描述

式中：θ是底滑面AB與水平方向的夾角，θ以逆時針方向為正；xA是A點的x坐標，yA是y點的y坐標，xB是B點的x坐標，yB是B點的y坐標，lAB是底滑面AB的長度。

（2）后緣滑裂面BC的坐標及長度描述

式中：β是后緣滑裂面與水平方向的夾角，β以逆時針為正；xC是C點的x坐標，yC是C點的y坐標，lBC是后緣滑裂面BC的長度。

（3）后緣滑裂面BC兩點的附加約束

式中：xD是D點的x坐標；L是巖質邊坡的寬度；d是邊坡基底的高度；H是巖質邊坡的高度。其中，xB在xD和L之間取間隔為1。

3.3 含剪切節理面巖質邊坡滑塊的非線性數學規劃模型

根據滑塊的平衡方程、屈服條件和幾何約束條件，可以建立約束非線性數學規劃模型，由于約束條件較多，為了確保求解方程的可行性，本文采取最優化方法中的全局最優解。

（1）建立目標函數

Maximize K

式中：K為巖質邊坡的安全系數；Maximize表示“使最大”。

（2）建立滑體的非線性約束平衡方程

式中：NAB是底滑面受到的法向力，NAB以受壓為正；SAB是底滑面AB的剪力，SAB以對滑體產生逆時針的轉動效果為正；NBC是后緣滑裂面BC的法向力，NBC以受壓為正；SBC是后緣滑裂面BC的剪力，SBC以對滑體產生逆時針的轉動效果為正。

（3）建立底滑面AB的非線性約束Mohr-Coulomb屈服條件

式中：K為巖質邊坡的安全系數；NAB是底滑面AB受到的法向力，NAB以受壓為正；SAB是底滑面AB的剪力，SAB以對滑體產生逆時針的轉動效果為正；φAB是底滑面AB的內摩擦角；cAB是底滑面AB的凝聚力；lAB是底滑面AB的長度。

（4）建立后緣滑裂面BC的非線性約束Mohr-Coulomb屈服條件

式中：K為巖質邊坡的安全系數；NBC是后緣滑裂面BC受到的法向力，NBC以受壓為正；SBC是后緣滑裂面BC的剪力，SBC以對滑體產生逆時針的轉動效果為正；φBC是后緣滑裂面BC的內摩擦角；cBC是后緣滑裂面BC的凝聚力；lBC是后緣滑裂面BC的長度。

（5）建立非線性數學規化模型

將上述目標函數安全系數K、滑塊受力的極限平衡方程式、底滑面和后緣滑裂面的屈服條件式以及各點坐標的附加幾何約束條件式聯立，即可得到復合地層巖質邊坡滑體的非線性數學規劃模型如下：

上式中，部分參數可以在邊坡中測量得知。

4 非線性數學規劃模型的求解

巖質邊坡后緣滑裂面位置穩定性計算方法的模型是一種比較典型的非線性數學規劃模型，本文的模型計算流程如下。

（1）定義目標函數：定義滑面的材料參數，抗剪斷凝聚力和內摩擦角度；定義幾何、荷載參數，包括用幾何法確定各點坐標和角度等確定部分并給予賦值，將未知部分作為未知數后續來進行求解。

（2）列出復合地層巖質邊坡滑體的非線性數學規劃的函數模型：先建立底滑面和后緣滑裂面的方程和滑體的重力計算公式；再根據上一步驟點的坐標確定底滑面lAB和后緣滑裂面lBC的長度；最后列出滑體的平衡方程和底滑面、后緣滑裂面的屈服條件以及部分額外幾何平面約束。

（3）利用MATLAB軟件中的功能函數進行求解，得到巖質邊坡的安全系數K和后緣滑裂面BC的具體位置，即得到復合地層巖質邊坡的后緣滑裂面位置的最優解。

5 基本算例分析

5.1 邊坡算例參數

本文選取的巖質邊坡具體參數為坡底寬L=100 m，坡高H=35 m，邊坡基底高度d=15 m，邊坡頂面寬度為67.19 m；臨空面AD的傾斜角度為68°，底滑面AB的傾角為30°，邊坡巖體的密度ρ=2 500 kg/m3，后緣滑裂面BC的傾角β為決策變量；對于滑體ABCD，其中，A點和D點坐標已知xA=19.920 8 m、yA=18.457 1 m、xD=32.806 2 m、yD=50 m，B點和C點坐標為決策變量。

在該算例中，由于滑塊長期只受重力荷載的影響，故底滑面AB和后緣滑裂面BC的凝聚力和內摩擦角均為固定值，具體參數見表1。

5.2 本文算法的邊坡算例計算結果

巖質邊坡滑體的非線性數學規劃模型是一個非線性數學規劃模型，其目標函數是安全系數，使用全局最優解求解巖質邊坡滑體的非線性數學規劃模型，可求解得到以下決策變量：K、xC、β。其中，xB取坐標間距為1，根據xB的確定坐標求解得到的決策變量如表2所示。

由表2可以看出，MTALAB最優化算法，在以xB取間隔為1且能得到最優解的情況下得到的安全系數最小為0.791 7，即文本文算法最終結果。

5.3 本文算法結果驗證

為了驗證最優化算法求解最優滑裂面和安全系數的可靠性與和合理性，采用OPTUM G2 2020軟件與Rocscience.Slide.v6007軟件進行對比驗算。本次驗算通過算例選取模型進行建模并施加標準邊界條件，選用OPTUM G2 2020軟件中極限分析的上下限法及強度折減的上下限法進行滑裂面搜索，但搜索的滑裂面并非平面雙折線型滑裂面，搜索得到的巖質邊坡總位移如圖3—6，圖中紅色實線為極限平衡最優搜索法結果，綠色多段線為軟件搜索滑裂面，為了方便對比參考，將軟件搜索的滑裂面取底滑面的右端點與滑裂面和頂面的交點的連線近似模擬同向雙平面型滑裂面，即藍色直線；選用Rocscience.Slide.v6007軟件中的畢肖普法可以完成非圓弧滑裂面（平面雙折線型滑裂面）的搜索，搜索得到巖質邊坡的安全系數K，底滑面和后緣滑裂面的A點和B點坐標及后緣滑裂面傾角β如表3所示。

圖3—6分別為G2軟件通過極限分析和強度折減法的上下限法搜索的滑裂面位置，可以看出，G2軟件搜索的滑裂面底滑面為直線形，但后緣滑裂面為近似圓弧形。在近似滑裂面與本文方法的結果滑裂面對比中，底滑面的長度誤差在0.5 m之內，后緣滑裂面誤差在1～2 m之間。從圖中的安全系數對比，基于極限平衡法的最優化搜索的安全系數低于G2軟件的各個方法所求的安全系數，最大差值為0.130 8，差距最小的是與極限平衡法的下限法的對比，誤差為8.36％，差距最大的是與強度折減法的上限法的對比，誤差為16.52％。通過對比近似滑裂面的傾角，極限平衡和強度折減法的上下限法與剛體平衡法的后緣滑裂面傾角相差3～4°，誤差為5.09％～6.78％。

表3為畢肖普法與全局最優搜索法結果，從表3可以看出，本文的全局最優搜索法與Slide軟件的畢肖普法搜索結果非常接近，誤差最大的依然為所計算的安全系數，誤差為11.04％。對于滑裂面位置，通過兩點坐標對比分析，誤差最大的為B點x坐標，誤差為5.88％，而C點坐標誤差僅為0.60％。后緣滑裂面傾角畢肖普法得到的結果為56.71°，本文算法得到的結果為58.91°，誤差為3.88％。

圖7—9為畢肖普法與全局最優搜索法的結果，從圖7—9可知，全局最優搜索法與Slide軟件的畢肖普法搜索結果非常接近。兩種方法的安全系數都隨后緣滑裂面凝聚力c的增大而增大，在后緣滑裂面凝聚力c=120 kPa時差值最大，差值為0.176 3，誤差為16.45％。對于滑裂面位置，通過兩點坐標對比分析，B點x坐標在后緣滑裂面凝聚力c=85 kPa時差值最大，差值為2.84 m，但對一個100 m寬的邊坡來說，誤差僅為2.84％，而C點x坐標在后緣滑裂面凝聚力c=115 kPa時差值最大，差值為1.81 m，誤差僅為1.81％。

5.4 結果分析與討論

通過算例驗證，我們不難發現本文算法的安全系數普遍偏低，相差最多為0.1左右，但在邊坡防護治理工作中，可以起到預防邊坡失穩的作用。由于計算方法和計算原理的不同，誤差必然是存在的，在圖3—7中我們可以發現極限平衡和強度折減的上下限法的滑裂面位置與本文算法搜索的滑裂面位置差距較大，這是因為OPTUM G2 軟件無法實現同向雙平面型滑裂面搜索，對于能夠實現同向雙平面型滑裂面搜索的Slide軟件，畢肖普法與全局最優搜索法（本文算法）結果非常接近。

通過觀察后緣滑裂面傾角，我們可以發現復合地層巖質邊坡的后緣滑裂面傾角接近60°，根據邊坡滑動破壞的普遍形式，平面滑動的邊坡巖體是沿單一地質斷層面（剪切節理面）的剪切位移，此時邊坡傾角β、地層交界面傾角θ及其內摩擦角φ之間的關系為β>θ>φ。本文算法的實施例以及Slide軟件的計算結果符合規定，進一步驗證本文算法的可行性與正確性。

6 結論

（1）本文以復合地層巖質邊坡為研究對象，基于非線性理論的非線性數學規劃模型計算復合地層巖質邊坡的穩定性及底滑面和后緣滑裂面位置，通過最優化方法中的全局最優搜索的計算結果，并與Slide軟件的結果進行對比，驗證了本文算法的可行性和準確性。

（2）傳統的極限平衡法計算邊坡的穩定性系數時，需要假定破壞面，然后再根據假定的每一個破壞面的安全系數，來確定邊坡最易發生破壞的滑裂面，強度折減法對c、tan φ簡單折減后得到的應力場已經不是原巖質邊坡的真實應力場，故而得出的結果也存在偏差，而本文的方法可以省略滑裂面位置的假設，直接通過全局優化算出邊坡的滑裂面位置，相比于極限平衡法，本文的編制程序更加簡單、工程應用簡便，可將其應用于含剪切節理面地層巖質邊坡的設計中。

（3）計算邊坡穩定性及邊坡滑裂面的工程軟件很多，主要以圓弧形滑裂面為主，但計算平面型滑裂面算法的軟件還需開發，本文算法程序簡單且計算精度較高，為同向雙平面滑動的邊坡提供了一種簡單高效的計算方法。在本文算法中，邊坡的滑裂面位置與Slide軟件的計算結果非常接近，安全系數比軟件的計算結果要普遍偏低，在工程應用中可以起到防范于未然的作用。

（4）本文提出了一種用于特殊巖質邊坡折線形滑裂面分析的計算方法，旨在解決現有軟件如OPTUM G2和Madis civil在處理此類問題時的不足。所開發的計算程序具有簡單直觀、原理明確且計算精度高的特點。本方法理論基礎扎實，且在工程實踐中操作簡便，適用于處理含有單個地質斷裂面（如斷層面或剪切節理面）的巖質邊坡的滑裂面位置分析。該算法不僅能夠準確預測滑裂面的位置，還能清晰地展示滑塊的自重、底部滑動面和后緣滑裂面的受力狀況，從而直觀地呈現邊坡的穩定性狀況。此外，算法提供的詳細信息對于邊坡支護結構的設計同樣具有重要價值。

參考文獻：

[1]黃潤秋. 20世紀以來中國的大型滑坡及其發生機制［J］. 巖石力學與工程學報， 2007（3）： 433-454.

［2］ 夏元友， 李梅. 邊坡穩定性評價方法研究及發展趨勢［J］. 巖石力學與工程學報， 2002（7）： 1087-1091.

［3］ 楊肖鋒， 魯祖德， 陳從新， 等. 板裂結構順層巖質邊坡滑移-彎曲破壞機制的力學模型研究［J］. 巖土力學， 2022， 43（增刊1）： 258-266.

［4］ 黃少平， 晏鄂川， 尹曉萌， 等. 不同臨空條件的層狀反傾巖質邊坡傾倒變形幾何特征參數影響規律［J］. 地質科技通報， 2021， 40（1）： 159-165.

［5］ LUO G， HU Q， TAN J. Process stability analysis for high slope based on limit equilibrium method and strength reduction finite element method［J］. Mining and Metallurgical Engineering， 2013， 33（2）： 14-17.

［6］ YAN L， SUN Y. The analysis on rock slope stability based on strength reduction technique［J］. Journal of Guangxi University（Natural Science Edition）， 2008， 33（3）： 235-238.

［7］ 張文蓮， 孫曉云， 陳勇， 等. 基于巖體抗壓強度折減的邊坡穩定性分析方法［J］. 巖土力學， 2022， 43（增刊1）： 607-615.

［8］ 李寧， 郭雙楓， 姚顯春. 再論巖質高邊坡穩定性分析方法［J］. 巖土力學， 2018， 39（2）： 397-406，416.

［9］ 李芬， 成濤. 基于雙折減系數法的巖質邊坡平面滑動分析［J］. 公路工程， 2018， 43（5）： 294-299.

［10］朱擘， 何光春. 強度折減法在平面滑動型巖質邊坡的應用［J］. 水運工程， 2012（9）： 65-69.

［11］程小龍， 張緒進， 朱擘. 雙強度折減法在平面滑動巖質邊坡中的應用［J］. 水運工程， 2013（3）： 72-76.

［12］陳建宏， 鐘福生， 陳定坤. 平面滑動型巖質邊坡極限平衡分析的上下限法［J］. 中南大學學報（自然科學版）， 2013， 44（8）： 3310-3315.

［13］CHENG Y， LANSIVAARA T， WEI W. Two-dimensional slope stability analysis by limit equilibrium and strength reduction methods［J］. Computers and Geotechnics， 2007， 34（3）： 137-150.

［14］ARDESTANI A， AMINI M， ESNAEILI K. A two-dimensional limit equilibrium computer code for analysis of complex toppling slope failures［J］. Journal of Rock Mechanics and Geotechnical Engineering， 2021， 13（1）： 114-130.

［15］SCHLOTFELDT P， ELME D， PANTON B. Overhanging rock slope by design： an integrated approach using rock mass strength characterisation， large-scale numerical modelling and limit equilibrium methods-ScienceDirect［J］. Journal of Rock Mechanics and Geotechnical Engineering， 2018， 10（1）： 72-90.

［16］張崇波. 基于LHS方法的巖質邊坡平面滑動可靠性分析［D］. 南京： 南京大學， 2015.

［17］陳志強， 王亮清， 劉順昌. 巖質邊坡平面滑動滑面極限傾角的研究［J］. 人民長江， 2013， 44（13）： 39-42.

［18］高丙麗， 李鐸， 李朗， 等. 基于坐標投影法巖質邊坡塊體穩定性分析及其可視化研究［J］. 巖土力學， 2022， 43（1）： 181-194.

Study on the Location and Stability of Slip Fracture Surface

of Rocky Slope with Shear Joint

Abstract：

Codirectional biplane sliding is one of the common failure modes of rock slopes with a single geological fault plane （shear joint face）， but the calculation methods for this type of slip surface are insufficient. In order to find the position of the slip surface of the slope more efficiently and accurately， and judge the stability of the slope， based on the limit equilibrium theory and nonlinear mathematical programming model， assuming that the slider has fallen off and the shape is an uncertain quadrilateral slide， and then assuming that the objective function is the safety factor of the rock slope， the position and stability of the slip fracture surface of the engineering slope with shear joint surface under natural working conditions are calculated by using MTALAB's global optimal search method. Then the results are analyzed and compared with those ofthe strength reduction method and the Bishop rigid body equilibrium method. The results show that the position of the slip fracture surface based on the limit equilibrium theory and the nonlinear mathematical programming model is basically consistent with the safety factor， which verifies the feasibility of such a method， providing a new basis for the calculation and stability analysis of the slip fracture surface of rock slope with a single geological fault layer.

Key words：

rock slope; calculation of slip surface; safety factor; limit equilibrium method; nonlinear theory; optimization principle