李同錄，王劉華，張常亮，李 萍

（長安大學地質工程與測繪學院，陜西西安710054）

0 引言

1915—1916年，Peterson在對瑞典Goetebborg采石場邊坡破壞形態(tài)研究的基礎上，首先提出了將均質土坡滑動面近似看作是一個圓弧的邊坡穩(wěn)定性分析方法，此后瑞典巖土工程學會對大量邊坡特別是鐵路工程邊坡的調查也進一步證實了該觀點［1］。目前已經有多種基于圓弧滑動面的邊坡穩(wěn)定性分析方法被提出，并在工程中廣泛應用［2］。

然而，對于一個特定邊坡，計算模型是在給定滑動面的情況下建立的。穩(wěn)定性評價還需要從無數個給定的圓弧滑動面中找出穩(wěn)定系數最小的一個作為臨界滑動面，這是一項繁瑣而復雜的工作。早期由于運算能力低，邊坡穩(wěn)定性分析主要采用圖解法。Fellenius首先提出了圓弧條分法，并利用圖解法繪制成圖表，以便實際應用［3］；Taylor提出的摩擦圓法也是基于圖解法制成圖表，這些圖表在中國的土力學教材中都有介紹［4］。

隨著計算機的應用，處理大量的滑弧滑動面已經不再是難題。當圓心和半徑大概的范圍給定，就可以對圓心區(qū)域網格化，對半徑按步長增加，計算每個圓心網格點和給定半徑的對應圓弧滑動面的穩(wěn)定系數。根據所有圓心網格點的穩(wěn)定系數便可繪制穩(wěn)定系數等值線，由此可得出穩(wěn)定系數的最小值及與之相對應的滑弧參數。如果斜坡是均質土坡，那么只有一個極值，但如果由不止一種土層組成，且各土層性質存在明顯區(qū)別，那么可能存在幾個極值，其中最小的一個就是所期望的值，目前普遍使用的商業(yè)軟件Geoslope就是采用這種算法［5－6］。一個圓有3個自由度，一般是圓心的x、y坐標和半徑R。在某些情況下，如粘聚力c趨于0時，半徑R可能趨于無窮大，x和y坐標遠離斜坡，因此不容易得到一個包括臨界圓中心的限制區(qū)域（圖1）。在這種情況下，以上方法是不嚴密的。針對該方法的不足，李同錄等提出了一種更為嚴格的二維搜索方法［7］，筆者將該方法進一步推廣到三維滑動面進行搜索。

1 臨界滑動面搜索的優(yōu)化算法

圖1 滑弧中心無窮遠與半徑無窮大的特殊情況Fig.1 Special Case That the Center and Radius of Slip Circle Is Far from the Slope

三維臨界滑動面的搜索以二維為基礎，因此先介紹二維方法。建立如圖2所示的直角坐標系，假定x軸的正方向指向坡內，z軸垂直，設點O（xO，zO）是給定滑弧的圓心（其中xO和zO分別代表點O的橫坐標和縱坐標，其他點類推），點A（xA，zA）和B（xB，zB）分別是剪入點和剪出點，通過A點做圓弧的切線并與x軸相交于點T（xT，0），首先，將圓心坐標xO，zO和半徑RO等3個參數用A、B和T等3個點的坐標表示。

圖2 由A、B兩點和切線AT確定的圓弧滑動面Fig.2 Illustration of Slip Surface Determined by Two Points A，Band Tangent Line AT

則圓弧滑面高程z的表達式為

滑坡坡面函數zs為

假定xA、xB給定，則zA、zB可以通過式（5）求得，且zT＝0。因此，滑弧實際是由xA、xB與xT3個參數控制的，且它們都可以控制在合理的有限域內。從數學和幾何方面而言，xA和xB的定義域分別為xA∈（xE，＋∞）、xB∈（－∞，xF），然而此處不難給出一個合適的xA的右邊界值和xB的左邊界值。假定這兩個邊界值分別為xD、xG，則xA和xB的定義域分別為xA∈（xE，xD）、xB∈（xG，xF）。對于xT，由于滑弧不可能越過最右側的AI線，所以它的右邊界值應該為xI＝xA，而它的左邊界值為xH。由于滑弧上的任何一點都不可能出現(xiàn)在邊坡坡面的上方，基于該限制條件，很容易確定xH的橫坐標，并得到xT的定義域xT∈（xH，xA）。對于任何圓弧條分法，穩(wěn)定系數F都可以表示為通用函數式中：γ為土的重度；c為粘聚力；φ為內摩擦角。 這3個參數反映土的屬性，可通過試驗結果或者經驗確定；坐標xA、xB與xT為變量。接下來找到xA、xB與xT的一個確定值，使穩(wěn)定系數極小，該問題用黃金分割法來求解。

如果邊坡土體是均勻的，F(xiàn)是凹函數，那么在整個循環(huán)過程中對應每一個xA、xB與xT分別應用黃金分割法，在循環(huán)運行結束后，就可以得到使F達到極小值所對應的xA、xB與xT的值，黃金分割法的整個過程如圖3。最后，由式（1）～（3）計算臨界滑弧的圓心和半徑。

圖3 用黃金分割法搜索土坡臨界滑動面的流程Fig.3 Process for Searching the Critical Slip Surface by Golden Block Optimization Method

如果邊坡由不同土層組成，那么F可能會存在局部極值。在這種情況下，黃金分割法可能會找到某一個極值點，但是并不一定是F的最小值。因此，首先按一定步長改變xA、xB與xT，找出所有局部極值域，然后在每一個極值域中運用黃金分割法求極小值，所有極小值中的最小值就是最終的結果。

基于上述討論，對于給定的土質邊坡，搜索其臨界滑動面的具體步驟。

（1）確定xA、xB與xT的合理區(qū)間，即xA∈（xE，xD）、xB∈（xG，xF）、xT∈（xH，xA）。

（2）按xA、xB、xT的順序，采用黃金分割法由外向內劃分黃金分割點，當到最里圈時就可以確定xA、xB與xT的一個設定。

（3）將xA、xB與xT代入式（1）～（4）計算xO、zO以及滑弧的半徑RO，并檢驗滑弧方程式z。

（4）將式（1）～（3）及相關土性參數代入式（6），計算F。

（5）重復以上（2）到（4），每一次循環(huán)都可以得到3個F值，通過比較這些值縮減相應變量的區(qū)間。

（6）最后，當xA、xB與xT的所有區(qū)間長度減小到給定誤差以內時循環(huán)結束。此時得到F的最小值及相應的xA、xB與xT值。

（7）將xA、xB與xT代入式（1）～（3）計算xO、zO和RO，并畫出相應的圓。

該算法與具體方法無關，對所有基于圓弧滑動面的方法都適用（如簡化Bishop法）。

2 三維臨界滑動面搜索

二維極限平衡法適合側向延伸遠且坡面形態(tài)簡單的邊坡，該類邊坡屬于平面應變問題，取其中一個單位寬度的窄條來分析。二維方法沒有考慮邊坡的側向約束力，代表了最保守的情況。如果邊坡地形條件復雜，如發(fā)育在溝谷或凸出的山梁上，邊界效應十分突出，采用二維方法可能低估了其穩(wěn)定性，此時應該采用三維方法。

對于均質土坡，三維臨界滑動面一般假定為旋轉橢球體。這一方面是大量破壞面觀察近似反映出該種形態(tài)特征，另一方面該假定保持了與二維圓弧滑動面的假定一致性，因為旋轉橢球的任一縱斷面都是圓（圖4），旋轉橢球體的方程可用下式表示

令式（7）的y＝0，則簡化為中軸面的圓弧方程。式（7）比圓弧方程多了一個橢球體的橫向半徑Rb，由此橢球滑動面的方程為

圖4 三維邊坡破壞面形態(tài)Fig.4 Fracture Sufrace Pattern of Three－dimensional Earth Slope

對于任一三維穩(wěn)定性計算方法，其穩(wěn)定系數可表示為

由式（9）可見，三維滑動面的搜索比二維多了一個自由度Rb，即旋轉橢球的橫向半徑。式（9）適合具有旋轉橢球滑動面的任何三維極限平衡法。

以一種具體方法為例，說明三維臨界滑動面的搜索。目前的三維極限平衡法，基本上是基于二維極限平衡法的擴展。張常亮等建立了三維通用極限平衡法［8］，給定特定的假定條件，可以簡化為具體方法。為了和二維簡化Bishop法比較，將通用法簡化為三維Bishop法［9］其計算公式為

式中：dW為微條柱的體力；dAz為微條柱底滑面的面積；dNz為微條柱底部法向力；u為孔隙水壓力；FM為基于整體力矩平衡的穩(wěn)定系數；αx為xOz平面上滑動面與x軸的夾角；αy為yOz平面上滑動面與y軸的夾角；αx、αy分別用滑動面方程求導來確定，即

Ry為橢球在垂直y軸、對應y取值的截面圓的半徑，由式（7）求得，即

當y＝0時

將以上參數代入式（10）就簡化為二維簡化Bishop法的公式，即

式中：mα＝cosαx＋tanφsinαx／FM。

式（10）的分子和分母項都可以用數值積分，如Gauss積分求解，xL、xH、yL、yH分別代表x和y方向的積分邊界。數值積分時，可不考慮潛在滑動區(qū)的實際邊界，積分邊界取x和y方向邊界（即橢球和地面的交線）的最大、最小值，相當于在一個矩形域求解。計算時對每一積分點做判斷，若地面高程小于橢球面，則不參與計算，其結果等效于按實際邊界積分。

二維圓弧滑動面搜索時，F(xiàn)M隨著xA、xB、xT的變化總能找到極小值。通過文獻［10］的一個算例來分析FM隨Rb的變化規(guī)律。該坡體為一均質土坡，土性和幾何參數如圖5，并假定橫向無限延伸。首先利用二維簡化Bishop法找出二維臨界滑動面，其圓心為（36.6，27.4），半徑RO為24.4m，穩(wěn)定系數為2.08。

圖5 用于三維計算的一個算例Fig.5 Worked Case for Three－dimentional Calculation

分別取Rb為RO的0.5、1.0、1.5、2.0、2.5、3.0、4.0、6.0倍，計算三維穩(wěn)定系數。給定Rb后，式（9）的自由變量與式（6）相同，除計算公式采用式（10）外，與二維算法相同（圖6）。

圖6 算例和案例中橢球長短軸之比與穩(wěn)定系數的關系Fig.6 Correlation Between the Ratios of Major and Minor Axises and the Factors of Safety for Worked and Real Cases

由圖6可見，邊坡的三維穩(wěn)定系數隨著滑動面橫向寬度的增大而減小，當長軸與短軸之比Rb／RO小于3時，三維效應十分明顯，三維穩(wěn)定系數隨著Rb／RO的增大而急劇減小；Rb／RO大于3時，穩(wěn)定系數的變化趨于平穩(wěn)；當Rb→∞時，三維FM趨近于二維FM。李同錄等利用三維Sarma法也得出了類似的結果［11］。

以上結果表明，三維穩(wěn)定系數隨著滑動面寬度增大而單調減小，不存在極小值，但存在極限值，當寬度無限大時，極限值逼近二維穩(wěn)定系數，二維穩(wěn)定系數是三維的下限。但三維穩(wěn)定系數隨滑動面寬度變化具有明顯的轉折點，若用橢球滑動面的橫縱半徑之比大于3時，邊坡穩(wěn)定性分析可簡化為二維問題；小于3時，應當按三維分析。前者如道路工程中開挖的長邊坡、河流侵蝕形成的直線岸坡等；后者如發(fā)育在溝槽或梁峁處的自然邊坡或者基坑的拐角處等。

橫向無限延伸的邊坡實際屬于二維問題，所以三維結果逼近二維。具有三維屬性的邊坡，邊坡破壞范圍可結合邊坡實際邊界條件加以界定，再利用橢球滑動面擬合，以便計算穩(wěn)定系數。給定約束條件后，也可能存在極值。

3 工程實例

現(xiàn)以西安咸陽國際機場污水排放口黃土高邊坡為例說明空間滑動面搜索方法的具體應用。該邊坡位于涇陽縣高莊鎮(zhèn)傅家村黃土塬的北緣，咸陽機場的污水經過凈化處理后，用暗管引到此處塬頂，然后從塬頂引到坡腳，再經明渠流入涇河。該邊坡高85m，坡度40°～50°，最陡處可達65°。斜坡北側為黃土塬，南側坡腳為涇河一級階地；東側為人工開挖的弧形高邊坡，坡腳建有泵站，用于汲取涇河水，坡頂為廢棄的磚瓦廠；西側為天然沖溝，寬30余米，深22～50m，溝頭周圍為墓地（圖7）。

圖7 西安咸陽國際機場排污口邊坡的衛(wèi)星影像Fig.7 Satellite Image of the Slope near the Outlet of the Sewage of Xi’an Xianyang International Airport

咸陽國際機場污水渠于2003年開始運營，此前未做詳細勘查。2006年3月春灌時，有流水加泥從排放口西約60m處坡腳的裂縫中流出，淹沒塬下部分麥田，在斜坡上也形成沖溝。隨后到塬頂查看，發(fā)現(xiàn)大量灌溉水流入一張開的裂縫，裂縫總長近100m，局部寬度達120cm，平行邊坡近東西向展布，距坡頂邊緣43m。調查當地村民得知，該裂縫在2001年3月春灌時就已出現(xiàn)，呈鋸齒狀延伸，出露長度約100m，寬10～60cm，停止?jié)驳睾罅芽p未繼續(xù)發(fā)展，并因農耕被掩埋。其后，該裂縫在每年灌溉時，仍時斷時續(xù)出現(xiàn)。為了防止灌溉水的進一步灌入，2006年春灌后，當地村民沿裂縫開挖溝槽，用隔水土工布和灰土夯填進行處理。但2008年春灌后，發(fā)現(xiàn)該裂縫再次張開，說明該斜坡處于不斷蠕變中。該斜坡滑動可能直接造成污水排放設施的損毀，因此引起了機場管理部門的高度重視，要求對其做穩(wěn)定性評價。

構成該邊坡的地層主要為第四系風積黃土夾古土壤層。黃土呈棕黃色，結構均勻，具大孔隙，垂直節(jié)理發(fā)育；古土壤呈褐紅色，具有粒狀結構。自坡頂到坡底可見L1—S9層黃土－古土壤序列，下伏河流沖積相的致密黏土。坡面一般覆蓋10cm至數米的坡積黃土狀土，邊坡地層剖面如圖8。

圖8 西安咸陽國際機場排污口邊坡地質斷面Fig.8 Geological Profile of the Slope near the Outlet of the Sewage of Xi’an Xianyang International Airport

該邊坡的變形主要由灌溉水下滲引起。黃土下部的致密黏土是相對隔水帶，灌溉水沿坡頂裂縫集中下滲，到該層頂面被阻滯，使地下水位大幅度上升，圖8斜坡上鉆孔揭露的地下水位高程位于坡腳高程以上，灌溉時水位可能會更高。因此該邊坡剪出口應位于潛水面附近的軟弱帶。

根據鉆孔試樣的物理力學性質測試結果，土的重度和有效強度參數取值為：γ＝18.7kN／m3；γsat＝19.8kN／m3；c′＝45kPa；φ′＝30°。

為了分析灌溉水對邊坡穩(wěn)定性的影響，分地下水位低和高兩種情況計算。前者是指地下水位在致密黏土中，即剪出口以下；后者指地下水位在坡腳以上的黃土中，按實測水位計算，地下水的影響范圍為灌溉區(qū)域，兩側的磚瓦廠和墓地沒有地表水補給，可認為地下水位深。

首先根據地形和地表裂縫展布確定中軸面，再給定不同的橢球長軸和短軸的比例計算穩(wěn)定系數（圖5）。根據滑坡兩側地形條件的限制和后緣裂縫的展布范圍，當RO為138.8m、Rb為267.5m時，滑動面橢球和地面的交線與后緣裂縫和側緣邊界最為吻合，因此應取Rb／RO為2.5對應的穩(wěn)定系數為實際穩(wěn)定系數值。由圖6可見，不考慮灌溉的情況下，邊坡穩(wěn)定系數為1.12，處于基本穩(wěn)定狀態(tài)；考慮灌溉導致地下水位升高的情況下，邊坡的穩(wěn)定系數為1.04，由于鉆孔所測水位不在灌溉期，灌溉期水位可能更高，邊坡穩(wěn)定系數在灌溉期更低，此時邊坡接近臨界狀態(tài)。由此可見，該邊坡的變形主要由農田灌溉引起。

為了分析潛在破裂面寬度對穩(wěn)定性的影響，圖6給出了不同的Rb／RO對應的穩(wěn)定系數，可以看出，不考慮灌溉時，穩(wěn)定系數的變化和圖5的算例規(guī)律相同，隨Rb／RO的增大，穩(wěn)定系數減小，Rb／RO大于2.5時，趨于一個穩(wěn)定值1.12。考慮灌溉時，總體規(guī)律也一樣，但在Rb／RO為2.5時，達到極小值1.04，隨后緩慢增大，此時的滑動面與實際邊界條件吻合最好，這說明在考慮實際地下水和地形等邊界條件的情況下，三維穩(wěn)定性計算也存在極小值。與不考慮地下水的情況比較表明，該極小值和地下水位局部上升有關，灌溉引起的高地下水位分布在張裂縫帶所在的有限范圍，潛在滑動面底部與地下水位面達到最大吻合時，穩(wěn)定系數達到極小，潛在滑動面再擴大，兩側水位降低，滑動面在水位線以上，則整體穩(wěn)定性提高，穩(wěn)定系數增大。由此可見，三維邊坡穩(wěn)定性評價時，現(xiàn)場確定邊界條件尤為重要。

4 結語

（1）均質土坡三維潛在滑動面可按旋轉橢球體來近似擬合，這和實際觀測基本吻合，也和二維圓弧滑動面的假定一致，由此可將二維滑動面的搜索方法推廣到三維滑動面。

（2）對三維旋轉橢球滑動面邊坡穩(wěn)定系數的計算結果表明，側向無限長邊坡，三維穩(wěn)定系數沒有極小值，但有極限值。因此，當邊坡側向延伸遠，使旋轉橢球長短軸之比大于3時，邊坡穩(wěn)定性分析可簡化為二維問題；當側向受地形條件或巖土條件限制，邊坡局限于一定范圍時，應按三維分析。

（3）按三維分析時，應結合邊坡的實際邊界條件，包括地形條件、巖土性質和地下水動態(tài)等，以便使邊坡潛在破裂面與邊界條件吻合，此時可能存在穩(wěn)定系數的極小值。

［1］ Taylor D W.Stability of Earth Dams［J］.Journal of Boston Society of Civil Engineering，1937，24（3）：197－246.

［2］ Bishop A W.The Use of Slip Circle in the Stability Analysis of Slopes［J］.Geotechnique，1954，5（1）：7－17.

［3］ Fellenius W.Calculation of the Stability of Earth Dams［C］∥United States Government Printing Office.Transaction of the 2ndCongress for Large Dams.Washington DC：United States Government Printing Office，1936：445－462.

［4］ 陳希哲.土力學地基基礎［M］.第2版.北京：清華大學出版社，1989.

［5］ 鄭 濤，張玉燈，毛新生.基于Geoslope軟件的土質邊坡穩(wěn)定性分析［J］.水利與建筑工程學報，2008，6（1）：6－8.

［6］ 王 勝，李云華，蘇 謙.基于Geoslope的高陡路基穩(wěn)定性分析［J］.路基工程，2008，20（5）：138－140.

［7］ 李同錄，鄧宏科，李 萍，等.搜索簡單土坡潛在滑動面的一種新方法［J］.長安大學學報：地球科學版，2003，25（3）：56－59.

［8］ 張常亮，李同錄，李 萍.三維極限平衡法通用形式的建立及應用［J］.地球科學與環(huán)境學報，2010，32（1）：98－105.

［9］ Hungr O.An Extension of Bishop's Simplified Method of Slope Stability Analysis to Three Dimensions［J］.Geotechnique，1987，37（1）：113－117.

［10］ Zhang X.Three－dimensional Stability Analysis of Concave Slopes in Plane View［J］.Journal of Geotechnical Engineering，1988，114（6）：658－671.

［11］ 李同錄，王艷霞，鄧宏科.一種改進的三維邊坡穩(wěn)定性分析方法［J］.巖土工程學報，2003，25（5）：611－614.