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

1 概述

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

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

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

2 邊坡安全系數(shù)計算方法

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

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

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

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

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

3.1局部搜索掃描算法

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

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

3.2全局搜索掃描算法

3.2.1多段折線模擬退火法

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

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

a.初始溫度選擇。

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

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

(1)

b.溫度下降方法。

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

Ti+1=αTi

(2)

c.算法終止條件。

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

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

3.2.2圓弧窮盡掃描算法

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

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

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

(3)

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

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

4 算例

4.1算例1

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

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

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

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

4.2算例2

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

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

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

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

4.3算例3

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

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

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

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

5 結(jié)論

本文采用Monte Carlo等角/等邊逐段掃描與模脫空時梁段的撓度變化呈縱向分布，有脫空發(fā)生時梁段的撓度變化呈橫向分布，無滑道脫空時梁段縱向兩邊的撓度基本為0，有些許下?lián)系内厔?，中間梁段的撓度有向上翹起的趨勢；滑道脫空后，從該滑道至臨近區(qū)域段的梁體呈下?lián)馅厔荩溆嗔后w出現(xiàn)上翹現(xiàn)象。單滑道脫空時邊滑道脫空引起梁體下?lián)陷^大，多滑道脫空中同側(cè)邊中滑道脫空引起梁體下?lián)厦黠@大于其他脫空工況，達到6 mm以上。有限元分析結(jié)果表明：箱梁頂推時，支撐反力及梁體受力對滑道高程差非常敏感，毫米級的誤差都將導(dǎo)致支撐反力發(fā)生劇烈變化，梁體局部有應(yīng)力集中的風(fēng)險，同時存在箱梁截面發(fā)生翹曲變形的可能。在實際頂推過程中，需對滑道高程偏差進行實時監(jiān)控，通過動態(tài)調(diào)整保持各滑道橫向始終處于同一高程，以減輕在頂推過程中由于滑道高程差導(dǎo)致的支撐反力差過大造成梁體局部受力集中的問題。

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

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

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

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