含軟弱土層邊坡穩定性分析的全局滑面自動搜索技術

(武漢工程大學 土木工程與建筑學院，湖北 武漢 430074)

1 概述

現有商業軟件[1]在進行邊坡穩定性分析時，一般采用滑面窮盡搜索算法，以獲取全局最小安全系數及其臨界滑面。在滑面選擇上，有圓弧滑面、多段折線滑面及多次樣條曲線滑面[2]。通常，這些商業軟件能較好地處理邊坡穩定性計算問題。但是，在面對具有軟弱下臥層的邊坡時，圓弧滑面假設在搜索時有可能無法找到全局臨界滑面；多次樣條曲線滑面一般是在圓弧極限滑面的基礎上進行局部微調得到的；此外，除非完全指定折線段的空間組成，商業軟件通常只能采用三段線搜索臨界滑面，且折線滑面需要工程師進行事先的初步預估，大致選定臨界滑面折點的可能位置及首尾折線傾角范圍，整個計算過程存在明顯的人為因素影響，效率不高。

Malkawi于2001年提出了多段折線臨界滑面的Monte Carlo隨機搜索技術[3]，其采用逐步逼近策略，在坡體地質較為簡單時具備全自動全局搜尋功能。張魯渝[4-5]將之改良后進行了各種邊坡的隨機轉角及常值轉角的穩定性分析計算，效果較好。白桃[6]也采用Monte Carlo隨機搜索技術進行了邊坡等線段及等角度穩定性分析計算，結果表明此方法在不太復雜土質邊坡中應用效果良好。而一旦在復雜土質邊坡中使用時，基于其從頭至尾的掃描特性，掃描結果會依賴于初始滑面位置，易導致計算結果離散，使得滑面落入局部極值中。

為解決上述問題，筆者引入模擬退火法(SAA)這一全局極值解算方法，結合改進后的Monte Carlo隨機搜索技術，采用Morgenstern-Price方法進行含軟弱土層邊坡滑面的全局自動搜索。

2 邊坡安全系數計算方法

一般而言，在進行邊坡穩定性的極限平衡法分析時，滑面條塊都會在納入邊坡拐點坐標的情況下進行水平方向上的等分。但在進行多層邊坡穩定性分析或者土體參數存在空間變異性時，條塊等分在計算上反而存在不便的地方。白桃[7]進行了邊坡截面四面形劃分的穩定性計算分析，但是對于含不規則土層截面的邊坡無法進行網格劃分。三角形具備良好的網格適應能力，本文以三角形網格劃分為例，進行邊坡條塊劃分說明。圖1為澳大利亞計算機應用協會(ACADS)發布的標準考題截面之一[8]。圖1(b)為圖1(a)在邊坡坡腳的局部放大圖。

提取多段折線滑面G(x)與邊坡截面三角網格相交后的所有交點，其與多段折線拐點合并后按x坐標排序，就可以獲得多段折線滑面條塊劃分點位。如圖1(b)所示，G(x)上任意相鄰交點(含拐點)i與j，其與頂面線S(x)之間的部分就構成一個單獨的計算條塊。i、j的中點O與S(x)之間的垂直距離設定為條塊高度。對于分層土體，每個條塊高度范圍內的土體重量進行疊加計算。O點所在三角形網格的強度參數取定為條塊底面的強度參數，用于滑塊的穩定性計算。

圖1 折線滑面條分法示意圖Figure 1 Diagram of Slices Method using n-segment Polyline

對于圓弧/非圓弧滑面的安全系數計算，數學方程較為完備的是Morgenstern-Price方法。為利于程序編制，本文采用朱大勇提出的簡明Morgenstern-Price法[9]。

3 邊坡臨界滑面的掃描算法

3.1局部搜索掃描算法

圖2 局部臨界滑面搜索掃描算法[6]Figure 2 Searching Strategy for Local Critical Slip Surface

圖2為白桃[6]進行多段折線滑面的局部臨界滑面掃描算法，其核心思想就是采用Monte Carlo隨機方法進行等角/等邊逐段搜索最優解的方法進行邊坡穩定的分析計算。從步驟1至步驟8為一輪，經多輪掃描一般可獲取均值邊坡的極限滑面。在步驟2和步驟8中，采用等角度規避邊坡地面線拐點，其余旋轉步驟均采用等長線段旋轉方式。這種方法存在明顯的短板，即面對含多層軟弱層土質邊坡時，容易陷入局部極值的問題。文獻[6]中有專門的論述。

3.2全局搜索掃描算法

3.2.1多段折線模擬退火法

為解決上述掃描中的局部極值問題，本文通過引入模擬退火法[10]。模擬退火法用于描述系統在溫度t下從能量E(i)狀態i進入具有能量E(j)狀態j的機制。t隨算法進程逐步遞減，類似于固體退火過程中的溫度角色。計算持續進行Metropolis算法：“產生新解—判斷—接受/舍棄”的迭代過程，對應著固體在某一恒定溫度下趨于熱平衡的過程。模擬退火算法從某個初始解出發，經大量解的變換最終求得最優化問題的整體最優解。本文對每次局部掃描獲得的極限滑面施加一定大小的抖動，以概率接受的方式進行全局極限滑面的位置篩選，計算流程見圖3。

模擬退火法是一種全局近似算法，在溫度下降過程中，系統接受惡化解的概率不斷降低，直到退火溫度后達到收斂解。在模擬退火的過程中，需要注意以下3個方面的技術問題：

a.初始溫度選擇。

退火初始溫度T需能保證平穩分布中每個狀態的概率趨于相等，如式(1)，本文選取初始溫度為0.5 ℃。

exp[(F1-F2)/T]≈1

(1)

b.溫度下降方法。

研究表明，退火溫度下降方法的選取直接關系到計算量的大小和計算結果的準確性。本文選用線性下降方法，降溫系數α=0.9，見式(2)。

Ti+1=αTi

(2)

c.算法終止條件。

模擬退火算法從初始溫度開始，通過在每一溫度的迭代和溫度的下降，最后達到終止原則而停止。本文設定的終止溫度為1×10-5℃，當退火溫度達到終止溫度時，認定掃描結束，此時的安全系數為全局最小安全系數。

圖3 模擬退火全局掃描計算流程圖Figure 3 Global searching flowchart using SAA

3.2.2圓弧窮盡掃描算法

為對比說明，對窮盡掃描算法[7]進行移植，分析本文采用的三角網格邊坡截面，搜索圓弧臨界滑面，掃描示意見圖4。

圖4 三角網格圓弧掃描示意圖Figure 4 Diagram of circular slip surface searching in triangular mesh slope section

圖4中，截面上方的方格交點構成搜索滑面的圓心集，L1和L2線段分別表示掃掠圓弧最下端的上限和下限。將L1至L2之間的距離L進行n等分，對于每一個圓心O，則需要進行n次掃描，每次掃描的半徑為：

(3)

式中:RL1為圓心到L1的垂直距離；i為針對某個圓心的第i次掃描。

提取圖4中圓弧與網格線的交點，每相鄰兩交點之間的土條作為極限平衡法的一個計算單元。圖4中ABCD就是一個計算單元，CD所在三角網格強度參數作為CD線上的強度參數。對每個圓弧滑面，采用朱大勇提出M-P簡明算法進行安全系數求解計算[9]。

4 算例

4.1算例1

圖5(a)為澳大利亞計算機應用協會(ACADS)發布的標準考題截面之一，截面尺寸、各層土體重度γ、粘聚力c及內摩擦角φ均在圖中標示。首先采用本文所提出方法進行邊坡的折線滑面穩定性分析，經試算可以確定當采用8節點7段線滑面時安全系數收斂，此時安全系數FsL=1.345。圖5(b)中虛線為任意隨機生成的8節點7段線初始滑面線。

(a) 邊坡截面尺寸及物理參數

(b) 搜索得到的臨界滑面圖5 算例1參數Figure 5 The Parameters of exapnple 1

同時，對本文圓弧窮盡掃描算法分析本文采用的三角網格邊坡截面，搜索圓弧臨界滑面，計算所得安全系數Fso=1.370。GeoStudio 2007商業軟件計算圓弧滑面安全系數Fsg=1.389。8節點臨界滑面及圓弧滑面見圖5(b)。圖5(b)說明本文提出的截面網格劃分方法具有良好的適應性，能夠針對含不規則形狀土層邊坡進行穩定性分析。雖然本例中多段折線滑面與圓弧滑面很接近，但也很明顯：多段折線滑面假設相對于圓弧滑面假設而言，能搜索得到更為危險的臨界滑面，因而更具有使用優勢。為清楚地顯示，后續示意圖不再畫出三角網格。

4.2算例2

本算例為兩土層中間夾雜一層軟弱斜坡層邊坡[11]，截面尺寸及物理強度參數見圖6(a)。各滑面計算安全系數見表1，各方法得到的臨界滑面見圖6(b)。

(a) 邊坡截面尺寸及物理參數

表1 各臨界滑面計算安全系數Table1 FactorofSafetyforCriticalSlipSurfaces滑面類型安全系數4節點0.422多段折線滑面5節點0.4136節點0.4117節點0.410圓弧滑面本文0.425GeoStudio20070.425

結果表明多段折線滑面在折點數量n=6時安全系數達到收斂。圖6(b)顯示了隨機生成的不同節點數量初始滑面掃描得到的臨界滑面。多次試算顯示，無論初始滑面位置如何，在給定節點數量的情況下，都可以搜尋得到相應的臨界滑面，表明本文所提方法對初始滑面位置不敏感。只要生成初始滑面，就可以完全由程序自身完成臨界滑面的全自動快速搜索。通過對比多段折線滑面和圓弧滑面，同樣發現多段折線滑面展現出較圓弧滑面更為優異的搜索能力。

4.3算例3

圖7(a)為文獻[6]中未能直接進行全局極限滑面搜索的復雜土層邊坡。在改良前，隨機生成的初始滑面在經歷多輪循環后極易落入局部極值，即圖7(b)中的局部臨界滑面S′，對應安全系數Fs′=1.82。因此，為完成此類復雜邊坡的極限滑面搜索，需要通過引入工程師經驗，預判滑面的可能位置，然后比較各可能位置計算之后的安全系數值大小，解算過程比較繁瑣。

(a) 邊坡截面尺寸及物理參數

(b) 搜索得到的臨界滑面圖7 算例3參數Figure 7 The Parameters of exapnple 3

引入模擬退火法之后，再次對這一復雜邊坡進行了多次試算。為方便說明，圖7(b)中的初始隨機滑面選為①②③，分別代表初始滑面位置跨越不同土層區間。計算結果表明，只需要隨機生成初始滑面的任意位置，而無論初始滑面落入哪個區間，本文算法總能得到收斂解，快速且自動搜索得到全局臨界滑面S，對應安全系數Fs=1.55。

5 結論

本文采用Monte Carlo等角/等邊逐段掃描與模脫空時梁段的撓度變化呈縱向分布，有脫空發生時梁段的撓度變化呈橫向分布，無滑道脫空時梁段縱向兩邊的撓度基本為0，有些許下撓的趨勢，中間梁段的撓度有向上翹起的趨勢；滑道脫空后，從該滑道至臨近區域段的梁體呈下撓趨勢，其余梁體出現上翹現象。單滑道脫空時邊滑道脫空引起梁體下撓較大，多滑道脫空中同側邊中滑道脫空引起梁體下撓明顯大于其他脫空工況，達到6 mm以上。有限元分析結果表明：箱梁頂推時，支撐反力及梁體受力對滑道高程差非常敏感，毫米級的誤差都將導致支撐反力發生劇烈變化，梁體局部有應力集中的風險，同時存在箱梁截面發生翹曲變形的可能。在實際頂推過程中，需對滑道高程偏差進行實時監控，通過動態調整保持各滑道橫向始終處于同一高程，以減輕在頂推過程中由于滑道高程差導致的支撐反力差過大造成梁體局部受力集中的問題。

擬退火法相結合進行含軟弱夾層邊坡的多段折線全局臨界滑面搜索。同時，編制了圓弧滑面窮盡掃描算法程序作為對比。以商業軟件Fortran和Matlab為計算手段進行聯合編程實現了Morgenstern-Price極限平衡搜索算法。計算結果表明：

a.截面三角網格劃分方法具有良好的適應性，能較好地對不規則層位幾何形狀進行離散，進而方便后續計算分析。

b.多段折線滑面具有較圓弧滑面更為優秀的極限滑面定位功能，其安全系數隨著滑面段數的增加會逐步達到收斂。

c.不再依賴于工程師經驗進行初始滑面大致范圍選擇，避免了人為因素影響，本文所提算法能快速地、全自動地進行由弱到強不同復雜土質情況的邊坡穩定性計算分析，并能找到全局臨界滑面與安全系數。