土質邊坡的單元失效概率與失效模式研究

彭 普，李 澤，張小艷，申林方，許 蕓

(1.昆明理工大學建筑工程學院，云南，昆明 650500；2.昆明理工大學電力工程學院，云南，昆明 650500；3.長江勘測規劃設計研究有限責任公司，湖北，武漢 430010)

邊坡的穩定性是與眾多隨機參數息息相關的，其中土體抗剪強度參數以及地下水的影響尤為重要[1]。天然土體的抗剪強度參數在空間上呈現變異性，研究表明土體參數空間變異性使得邊坡的失效模式和安全性能均具有一定的不確定性，因此需要采用可靠度方法進行研究[2-3]。地下水從兩個方面影響著邊坡的穩定性：一是其使得邊坡內部的滲流場發生變化，滲流場的變化導致邊坡孔隙水壓力發生改變，從而影響著邊坡的安全性能；二是其降低邊坡土體的材料抗剪強度參數，從而影響著邊坡的安全性能[4]。通常情況下地下水位并不是一個確定值，其受眾多因素影響。因此，邊坡的穩定性問題是一個與土體參數空間變異性和地下水位隨機性相關的不確定性問題，在進行邊坡穩定性分析時應綜合考慮這兩個因素帶來的影響。

近年來，邊坡系統可靠性問題受到了廣大學者的關注。基于極限平衡分析方法(LEM)可以直接得到邊坡穩定性系數的數學分布以及邊坡的失效概率，但LEM 在進行邊坡穩定性計算時首先需要假定滑裂面，然后對滑裂面上的土體進行條分并假定條間力的大小及方向將超靜定問題轉化為靜定問題進行求解[5-8]；基于有限元分析方法(FEM)可以得到邊坡的應力-應變分布，還可以獲得比LEM方法更合理的穩定性分析結果，但計算成本較高[9-11]；基于有限元塑性極限分析方法可以在不考慮邊坡加載歷史和材料本構關系的情況下，高效、準確地得到邊坡處于極限狀態下的穩定性系數和失效模式[12-18]。該分析方法主要通過構建靜力許可應力場(下限法(LBM))和機動許可速度場(上限法(UBM))，通過塑性力學的上、下限定理可知，真實解必然位于上限解與下限解之間。張小艷等[19]、李亮等[20]將土體抗剪強度參數設為隨機變量，提出了基于塑性極限分析的邊坡可靠度上、下限法；陳朝暉等[21]提出了考慮參數空間變異性的邊坡可靠度上、下限分析方法。以上研究促進了有限元塑性極限分析方法在邊坡可靠度中的應用。然而，隨機地下水位對邊坡可靠度的影響并未得到系統的研究。

此外，傳統的邊坡失效概率分析采用整體失效概率法，其主要采用邊坡穩定性系數的閾值判斷邊坡是否失效，并據此計算邊坡的失效概率，該方法隱含假定邊坡具有單一的失效模式，這與邊坡存在多種失效模式不相符。為解決這一問題，李典慶等[22]、蔣水華等[23]提出了基于代表性滑動面的邊坡失效概率計算方法，然而在進行邊坡失效模式統計時，隨著計算樣本數的增加計算量會急劇增大，其計算效率和計算精度均較低[24]；張小艷等在文獻[19]中提出了基于上限法的單元失效概率計算方法，利用邊坡穩定性系數和速度場進行邊坡單元失效概率的計算，為邊坡系統失效概率分析提供了一種新思路。

鑒于此，本文將有限元離散技術、上限法理論、相關非高斯隨機場模擬方法以及隨機規劃理論結合起來提出了隨機地下水位作用下考慮參數空間變異性的邊坡可靠度分析方法。

1 隨機地下水位作用下考慮參數空間變異性的邊坡可靠度分析

1.1 隨機地下水位的邊坡穩定滲流場

當前，地下水位通常被設定為某一確定值，在進行邊坡可靠性分析時，邊坡的極限狀態函數不包含地下水位這一隨機變量。然而由于土體顆粒和孔隙在邊坡中隨機分布使得邊坡滲流場具有不確定性，此外降水、地表徑流、地下匯流等對地下水的供給以及蒸發、抽水等對地下水的排出的不確定性使得邊坡地下水位具有不確定性，不確定的地下水位將導致不確定的滲流場。近年來，對于滲流場的研究主要集中在土體材料變化導致的滲透系數變化，進而導致滲流場的變化上[25-26]。對于隨機地下水位對滲流場的影響以及對邊坡可靠性的研究成果依然較少。

本文將地下水位定義為隨機變量，如圖1 所示，假定其服從正態分布[27]，并假定地下水位的上界和下界，地下水位位于上界和下界之間隨機變化。圖1 中：為水位的下界；為水位的上界；為水位的均值。為了簡化計算，本文假設隨機地下水位作用下的滲流場為飽和穩定滲流場，并且不考慮水位突然升高或降低引起的超孔隙水壓力以及對土體抗剪強度參數的影響，且位于浸潤線以上的孔隙水壓力p=0，位于浸潤線以下的孔隙水壓力p=γh，其中 γ為水的容重，h為浸潤線上與計算點處于同一條等勢線上的點到計算點的垂直距離。采用二維穩定滲流公式進行滲流分析[28]，具體公式如下：

圖1 邊坡隨機地下水位作用示意圖Fig.1 Schematic diagram of random groundwater level of soil slope

式中：kx為土體材料x方向的滲透系數；ky為土體材料y方向的滲透系數；Hr為土體內各點的隨機水頭函數。

1.2 描述參數空間變異性的隨機場模型

當前普遍采用自相關函數表征土體的空間變異性，常用的自相關函數類型有指數型、高斯型、對數型、二階自回歸型、指數余弦型以及三角型[29]。考慮邊坡的可靠性對自相關長度較為敏感，對自相關函數的形式不敏感[30-31]，本文選用數學表達式更為簡單的指數型自相關函數進行研究，具體可采用式(2)表示：

式中：xa、xb為空間坐標；xaa、xba、xab、xbb為空間坐標的坐標分量；Lh為水平坐標方向的波動范圍；Lv為豎直坐標方向的波動范圍[31]。

由于土體參數之間存在一定的互相關性和土體參數本身存在一定的自相關性，因而涉及相關非高斯場的離散過程，本文采用類似文獻[32]的方法，采用基于喬列斯基分解的中點法進行隨機場的生成，對于指數型自相關函數，土體黏聚力、內摩擦角的相關對數隨機場可表示為：

式中：m=c,φ，c為土體黏聚力，φ為土體內摩擦角；(xa,ya) ∈?；Hm(xa,ya)為相關非高斯隨機場；為相關標準高斯隨機場；σlnm=為對數正態分布m的均值， σm為對數正態分布m的標準差；ulnm為相應正態變量lnm的均值；σlnm為相應正態變量lnm的標準差。

1.3 邊坡可靠度分析的隨機規劃模型

有限元塑性極限分析兼具FEM 法和極限分析的優點，成功克服了結構體在進行極限分析時靜力許可速度場的困難。SLOAN 等[14-15]將土體采用三角形單元進行離散，通過構建三角形單元塑性流動約束條件、公共邊速度不連續約束條件、三角形單元速度邊界約束條件從而構建機動許可速度場。由上限定理易知，機動許可速度場和外荷載是一一對應的關系，其中最小外荷載無限接近于極限荷載，因此上限法可理解為求解極小值的數學規劃問題。本文在SLOAN 等[14-15]、張小艷等[19]、王均星等[33]研究的基礎上，建立隨機地下水作用下考慮土體參數空間變異性的邊坡可靠性分析隨機規劃模型如下：

式中：Z為邊坡可靠度計算的極限狀態函數；km為邊坡穩定性系數；kγ為土體容重超載系數；We為有限單元的內功功率；Wd為有限單元公共邊的內功功率；WG為自重在有限單元節點速度上所做的外功功率；為孔隙水壓力在有限單元連續體內所作的外功功率；Wdp為孔隙水壓力在有限單元公共邊上所作的外功功率；λe為塑性乘子；λd為決策變量；為有限單元e的塑性流動約束條件矩陣；為相鄰有限單元公共邊d的塑性流動約束條件矩陣；ab為有限單元速度邊界約束條件矩陣；ue為有限單元e的速度矢量；ud為相鄰有限單元公共邊d的速度矢量；ub為邊界上有限單元的速度矢量。

采用三角形單元進行邊坡離散，在隨機地下水位作用下考慮土體參數空間變異性時，不同空間位置的單元將具有不同的抗剪強度參數，式(4)中有限單元的內功功率為：

式中：ce為有限單元e形心處的黏聚力；φe為有限單元e形心處的內摩擦角；Ae為有限單元e的面積，e=(1,2,···,Ne)，Ne為所有有限單元數量。

有限單元公共邊的內功功率為：

式中，ld為公共邊的長度，d=(1,2,···,Nd)，Nd為所有有限單元公共邊的數量。

自重所作的外功功率為：

孔隙水壓力所作的外功功率為：

2 隨機規劃模型的求解策略

式中：tw=(1,2,···,nw)，nw為邊坡地下水位蒙特卡洛隨機數的數量；為邊坡地下水位的第tw個隨機數；μw為邊坡地下水位的均值； σw為邊坡地下水位的標準差；Random為截尾正態分布隨機數生成函數；Normal為地下水位隨機數符合正態分布；Hlb為地下水位的下界，取最低水位；Hub為地下水位的上界，取最高水位。

2)生成邊坡材料抗剪強度參數的隨機場。假設土體材料的黏聚力cr和內摩擦角φr的自相關函數類型為指數型，采用喬列斯基分解的中點法進行隨機場的生成：

本文采用容重超載的方式使邊坡達到極限狀態，隨機地下水位作用下考慮參數空間變異性時，邊坡的容重超載系數為：

式中：kγ(tm,tw)為第tw個地下水位作用下對應著第tm個隨機場的邊坡的容重超載系數；、為土體處于極限狀態時的容重，其與以及相關； γa為土體實際容重；為強度折減以前第tm個非高斯隨機場中有限單元e形心處的黏聚力和內摩擦角。

隨機地下水位作用下考慮參數空間變異性時，邊坡穩定性系數如下：

式中：tw=(1,2,···,nw)；tm=(1,2,···,nm)；、為強度折減以后第tm個非高斯隨機場中有限單元e形心處的黏聚力和內摩擦角。

3)計算邊坡穩定滲流場。將步驟1)生成的地下水位的隨機數tw=(1,2,···,nw)，逐次代入邊坡穩定滲流場的計算公式，進而得到離散單元e三個節點的孔隙水壓力值：，e=(1,2,···,Ne)，tw=(1,2,···,nw)。

圖2 邊坡可靠度上限法流程圖Fig.2 Flow chart of reliability analysis using upper bound method for slope

5)計算邊坡穩定性系數并繪制相關曲線。

3 邊坡系統失效概率和可靠度指標

邊坡工程中普遍采用邊坡穩定性系數來評價邊坡的穩定性。當邊坡穩定性系數km≥1，則認為邊坡穩定；邊坡穩定性系數km<1，則認為邊坡失穩[34]。因此，邊坡的失效功能函數按式(14)進行表示：

式中：tw=(1,2,···,nw)；tm=(1,2,···,nm)；I(tw,tm)為第tw個地下水位作用下對應著第tm個隨機場的邊坡失效功能函數。

根據式(14)邊坡失效功能函數，可以得到邊坡在第tw個地下水位作用下的失效概率：

式中：tw=(1,2,···,nw)；為邊坡在第tw個地下水位作用下的整體失效概率。

式(15)被廣泛應用于邊坡失效概率的計算，然而隱含地假定了邊坡僅具有單一失效模式，這與邊坡存在多種失效模式是不相符的[35-36]，此外該式只考慮了邊坡穩定性系數的影響并未考慮邊坡的失效后果的影響。為克服上述不足，張小艷等[19]提出了單元失效概率法用以分析邊坡的失效概率。該方法充分發揮有限元離散技術在構建速度場上的優越性，認為在每一個失效模式對應的速度場中，當單元速度大于0 時代表此單元已發生塑性流動(失效)、當單元速度等于0 時代表此單元未發生塑性流動(安全)。在此基礎上，本文首次提出根據邊坡中單元的失效信息來估算邊坡的系統失效概率。邊坡系統是由無數個潛在滑動面組成的串聯系統，其系統失效概率近似等于各個代表性失效模式的失效概率之和。當前，采用LEM 法和FEM 法進行邊坡系統可靠度計算時任意潛在失效模式(設為N)均會導致邊坡發生失穩破壞，因此邊坡系統失效的功能函數為：

式中：Zs為邊坡系統失效的功能函數；Zi為邊坡的第i種潛在的失效模式，i=(1,2,···,N)。

根據邊坡系統失效功能函數的定義，可以得到邊坡的系統失效概率估算公式為：

式(17)在進行邊坡系統失效概率計算時需要在大量的潛在失效模式中識別最具有代表性的失效模式，因此，所識別的代表性失效模式對計算結果有較大的影響。塑性極限分析上限法可以直接得到邊坡極限狀態下的失效模式，每種失效模式均是由若干個失效單元組成的區域，鑒于此，本文定義邊坡單元的失效功能函數如下：

式中：Ie(tw,tm)為第tw個地下水位作用下對應著第tm個隨機場中單元e的失效功能函數；uec(tw,tm)為第tw個地下水位作用下對應著第tm個隨機場中單元e的速度。

從分類角度看，采用邊坡穩定性系數的閾值和有限單元的節點速度信息可以識別邊坡所有失效模式。AP(Affinity Propagation)聚類算法非常適合于本文邊坡失效模式的識別，進而進行邊坡系統失效概率的估算，該算法的主要思想是將tw×tm個速度場的失效面積都當作潛在的聚類中心，各失效面積之間連線構成一個網絡，再通過網絡中各條邊的消息傳遞計算出各樣本的聚類中心。第m、j個速度場的失效面積Am、Aj，Aj作為Am聚類中心的能力，記為S(Am,Aj)，S(Am,Aj)越大，兩個點距離越近；定義吸引度R(Am,An)為An作為Am聚類中心的適合程度，如圖3(a)所示為Am選An的過程；定義歸屬度B(Am,An)為Am選擇An作為其聚類中心的適合程度，如圖3(b)所示為An選Am的過程。R(Am,An)與B(Am,An)的和越大，表示An作為聚類中心的可能性越大，Am隸屬于An為聚類中心的可能性也越大。具體的，吸引度迭代公式為：

圖3 AP 聚類算法選點過程Fig.3 AP clustering algorithm selection process

式中：

Rt+1(Am,An)為更新后的R(Am,An)；Rt(Am,An)為更新前的R(Am,An)；λ 為阻尼系數。

歸屬度迭代公式為：

式中：

Bt+1(Am,An)為更新后的B(Am,An)；Bt(Am,An)為更新前的B(Am,An)。

定義單元失效概率為：

根據邊坡可靠度分析上限法數學規劃模型的求解策略中所得nw×nm個邊坡穩定性系數km(tw,tm)，通過統計可以得到邊坡在第tw個地下水位作用下穩定性系數的均值、標準差以及在所有可能發生的地下水位作用下邊坡穩定性系數的均值、標準差。具體為：

式中：tw=(1,2,···,nw)；μk(tw) 、σk(tw)分別為第tw個地下水位作用下邊坡nm個穩定性系數的均值和標準差；μk、 σk分別在所有可能發生的地下水位作用下nw×nm個邊坡穩定性系數的均值和標準差。

4 算例分析

根據本文提出的方法，編制了相應的上限法程序，對一個經典的邊坡算例進行了計算分析并與極限平衡法、有限元法計算結果對比，驗證了本文計算方法的正確性。

4.1 均質邊坡基本信息

本文選取文獻[37]的均質邊坡算例，來驗證本文計算方法的正確性。圖4 為均質邊坡計算模型，邊坡高度為10 m，坡比為1∶1，坡頂寬度為10 m，諸多文獻對此邊坡的失效概率進行過分析[29,32,37]。本文采用三角形有限單元對邊坡進行離散，共得到885 個有限單元、2655 個有限元節點、1279 條公共邊。邊坡左側隨機地下水位為，右側坡腳的地下水位與地表齊平，設置P1、P2、P3 三個孔隙水壓力關鍵點(坐標分別為：P1(5, 5)、P2(10, 5)、P3(15, 5))；設置有4 個有限單元失效概率監測單元，分別為坡頂的E1 單元、邊坡中部的E2、E4 單元、坡腳的E3 單元，其中單元E1、E2、E3 均在邊坡的臨空面上，E1、E4 在同一豎直面上，E2、E4 在同一高程上。設定土體容重為γ=20 kN·m-3、滲透系數為K=5×10-7m/s，兩者均為確定值，其他計算參數見表1。本算例的計算目的是：1)計算隨機地下水位作用下考慮參數空間變異性的邊坡穩定性系數分布類型；2)計算隨機地下水位作用下考慮參數空間變異性的邊坡穩定性系數概率密度曲線和累積概率密度曲線；3)采用AP 聚類算法進行邊坡失效模式的識別并進行邊坡系統失效概率的估算；4)計算邊坡的整體失效概率以及單元失效概率與地下水位之間的關系。

表1 土體參數統計特性Table 1 Statistics of soil parameters

圖4 均質邊坡計算模型 /mFig.4 Homogeneous slope calculation model

4.2 隨機地下水位和滲流場

圖5 地下水位隨機數分布直方圖Fig.5 Histogram of random distribution of groundwater level

圖6 邊坡穩定滲流場Fig.6 Steady seepage field of slope

圖7 關鍵點處孔隙水壓力隨機變化情況Fig.7 The pore water pressure of the monitoring points

4.3 土體抗剪參數空間變異性和隨機場

本文假定土體黏聚力c和內摩擦角φ服從對數正態分布，采用指數型自相關函數模擬土體抗剪強度參數的空間變異性。土體參數隨機場的數量nm取1000；根據式(10)、式(11)，生成邊坡材料抗剪強度參數的1000 個隨機場。圖8 為tm=500時邊坡抗剪強度參數隨機場，在該隨機場下土體黏聚力的最大值為17.7961 kPa，最小值為3.5555 kPa；內摩擦角的最大值為46.8418°，最小值為17.8413°。此外，土體黏聚力和內摩擦角在空間上呈現一定的負相關。

圖8 邊坡抗剪強度參數隨機場Fig.8 Random Field of Slope Shear Strength Parameters

4.4 上限法與LEM 法、FEM 法的對比分析

文獻[37]采用150 項K-L 級數展開的方法離散土體的隨機場，通過5×104次MCS 模擬，得到邊坡的可靠度；文獻[32]采用基于喬列斯基分解的中點法(CD)模擬隨機場，采用拉丁超立方抽樣(LHS)產生土體的黏聚力和內摩擦角樣本，進而進行邊坡可靠度的計算。為驗證本文計算方法在無地下水作用時的有效性，邊坡計算模型以及土體參數均與文獻[32, 37]一致，根據邊坡可靠度分析隨機規劃模型式(4)和隨機規劃模型的求解策略。當無地下水作用時，離散單元三個節點的孔隙水壓力值：均為0，根據第4.3 節土體參數隨機場數量nm=1000，可以得到1000×1=1000 個計算樣本，基于開發的并行計算程序，在一個小型工作站(處理器：AMD 銳龍3970X、32 核心，物理內存：128 GB)上計算1000 個樣本花費了大約0.43 h(單線程計算約17.2 h)，平均每個樣本的計算時間為1.56 s(單線程計算約61.92 s)。此外，使用有限元極限分析軟件OptumG2 的FEM法對該邊坡進行了穩定性分析，計算相同樣本共花費了大約23.4 h，平均每個樣本計算時間為84.24 s。本文并行計算程序相較于單線程計算效率提高約38.69 倍，相較于FEM 法計算效率提高約53 倍。計算結果如表2 所示，本文方法與K-L 法、LHS法相比，所得邊坡穩定性系數的均值、標準差基本相同，所得失效概率略微偏小；本文方法與FEM法相比，所得邊坡穩定性系數的均值、標準差基本相同，所得失效概率略微偏大。

表2 不同方法的邊坡可靠度結果對比Table 2 Comparison of slope reliability results of different methods

為驗證本文計算方法在地下水作用時的有效性，本文選取低、中、高三個地下水位進行對比分析。表3 為本文計算方法與LEM 法、FEM 法計算結果對比情況，其中選取的水位分別為tw=1、tw=25、tw=50。圖9 為三種水位下邊坡穩定性系數分布。

表3 三種水位作用下邊坡的可靠度指標Table 3 Calculation results of reliability index

圖9 邊坡穩定性系數分布特征Fig.9 Distribution characteristics of slope safety factor

計算結果表明：

1)采用UBM 計算所得邊坡穩定性系數均值的上限解大于LEM 法、小于FEM 法所得計算結果，但誤差較小。隨著水位的升高三種計算方法所得邊坡穩定性系數均降低。在tw=1、tw=25、tw=50三種情況下，UMB 法所得邊坡穩定性系數均值均大于相同條件下LEM 法所得的均值，符合上限解的特征。此外，隨著地下水位的升高，邊坡整體失效概率隨之增加，這與實際是相符的。在三種地下水位情況下，采用UBM 法計算所得邊坡整體失效概率均小于LEM 法所得計算結果，但誤差較小，符合上限解的特征，根據上限定理：機動許可速度場對應的荷載是極限荷載的上界，即采用上限法所得邊坡穩定性系數必然大于真實邊坡穩定性系數，因此在統計邊坡的整體失效概率時，會略微低估邊坡的失效概率。然而，在tw=1、tw=25、tw=50三種情況下，UMB 法所得邊坡穩定性系數均值均小于相同條件下FEM 法所得的均值，其原因主要是使用有限元極限分析軟件OptumG2在進行邊坡穩定性分析時，依據計算不收斂、塑性區貫穿、計算點位移發生突變作為邊坡處于極限狀態的判斷標準，因此FEM 法所得邊坡穩定性系數是嚴格的上限解還是下限解無法驗證。

2)圖9(a)為tw=1、tw=25、tw=50三種情況下，UBM 法與LEM 法、FEM 法所得邊坡穩定性系數累積概率密度曲線，不難看出，UBM 法與LEM法所得邊坡穩定性系數累積概率密度曲線都非常接近，誤差較小；UBM 法所得邊坡穩定性系數累積概率密度曲線都位于FEM 法左側，其原因主要是UBM 法和FEM 法在進行邊坡穩定性系數求解時所采用的方法不同，前者通過求解邊坡可靠度上限法數學規劃模型，后者依據邊坡處于極限狀態的判斷標準。此外，隨著地下水位的升高，累積概率密度曲線逐步向左移動。圖9(b)、圖9(c)、圖9(d)分別為tw=1、tw=25、tw=50三種情況下UBM 法所得邊坡穩定性系數頻數分布直方圖，均呈現中間高兩側低的規律，與地下水位頻數分布直方圖相似。在tw=1時邊坡穩定性系數在1.2 附近內出現的頻數最多，最高頻次為172 次；在tw=25時邊坡穩定性系數在1.15 附近內出現的頻數最多，最高頻次為164 次；在tw=50時邊坡穩定性系數在0.9 附近內出現的頻數最多，最高頻次為165 次。

4.5 邊坡穩定性系數分布規律

根據邊坡可靠度分析隨機規劃模型式(4)以及隨機規劃模型的求解策略，共得到50 個地下水位在1000 個隨機場下的邊坡穩定性系數。圖10、圖11分別為邊坡穩定性系數的50 條概率密度曲線以及累積概率密度曲線；圖12 為邊坡穩定性系數的均值與地下水位的關系。通過分析可得到如下規律：

圖10 邊坡穩定性系數概率密度曲線Fig.10 Probability density curve of safety factor of slope

圖11 邊坡穩定性系數累積概率密度曲線Fig.11 Cumulative probability density curve of safety factor of slope

1)邊坡穩定性系數的分布與正態分布一致，且隨著地下水位的升高，邊坡穩定性系數的均值逐漸降低，邊坡穩定性系數的概率密度曲線和累積概率密度曲線均向左移動。此外，隨著地下水位的升高，邊坡穩定性系數的標準差逐漸降低，邊坡穩定性系數的概率密度曲線分布范圍逐漸變窄，累積概率密度曲線逐漸變陡。

2)邊坡穩定性系數的均值與地下水位直接相關，采用多項式擬合得到邊坡穩定性系數的均值與地下水位的量化關系式為：

式(27)表明邊坡穩定性系數的均值與地下水位成負相關。

3)對50 個地下水位在1000 個隨機場時所得的50 000 個邊坡穩定性系數進行統計分析，得到考慮土體參數空間變異性以及地下水位隨機性的邊坡穩定性系數的分布特征，圖13 為隨機地下水位作用下邊坡穩定性系數的頻率分布直方圖。隨機地下水位作用下邊坡穩定性系數的均值為1.127，標準差為0.097。

圖13 邊坡穩定性系數的頻數分布直方圖Fig.13 Frequency distribution histogram of safety factors

4.6 邊坡系統失效概率分析

通過50000 次的模擬，LEM 法默認邊坡僅存在一種失效模式(如圖14 中的LEM 滑動面)，FEM法得到邊坡的4 種失效模式(如圖14 中的FEM 滑動面)。本文UBM 法相較于LEM 法、FEN 法最大的優勢在于可以根據邊坡中單元的失效信息采用AP 聚類分析方法捕捉到邊坡所有的失效模式(如圖14 中的失效區域)，進而進行邊坡系統失效概率的估算，圖14 為邊坡的6 種失效模式AP 聚類結果，表4 列舉了邊坡的6 種失效模式以及對應的失效概率，表5 為tw=1、tw=25、tw=50時邊坡失效模式以及對應的失效概率。由結果可知，失效模式1、失效模式2 失效區域面積較小，屬于淺層滑坡，失效面積為5.45 m2～15.52 m2(失效次數為1193)；失效模式5、失效模式6 失效區域面積較大，屬于深層滑坡，失效面積為89.45 m2～112.65 m2(失效次數為192)；失效模式3、失效模式4 介于淺層滑坡和深層滑坡之間，失效面積為48.72 m2～72.42 m2(失效次數為4430)。

表4 邊坡失效模式以及對應的失效概率Table 4 Slope failure mode and the corresponding failure probability

表5 不同地下水位作用下邊坡失效模式以及對應的失效概率Table 5 Slope failure modes and corresponding failure probabilities under the action of different groundwater levels

圖14 邊坡失效模式AP 聚類結果Fig.14 Schematic diagram of AP clustering results of slope failure mode

地下水位tw=1時，邊坡僅具有單一失效模式，具體為失效模式3(失效次數為8，失效概率為0.8%)；地下水位tw=25時，邊坡的失效模式變為兩種，具體為失效模式3(失效次數為13，失效概率為1.3%)、失效模式4(失效次數為12，失效概率為1.2%)；地下水位tw=50時，邊坡的6 種失效模式均可能發生。其主要原因是：隨著地下水位的升高，邊坡浸潤線相應上移，邊坡土體中飽和區域的面積增加，從而增加了失效模式4 的發生概率。此外，由于邊坡較陡，隨著地下水位的進一步升高，邊坡土體中飽和區域的面積進一步增加，地下水從坡面溢出可能性增加，從而進一步增加了失效模式1、失效模式2、失效模式5、失效模式6 的發生概率。

LEM 法所得的破壞模式僅與本文結果的失效模式3 相吻合，其主要原因在于：采用LEM 法進行邊坡穩定性計算時，事先假定初始滑裂面，并根據抗剪參數的均值搜索得到臨界滑裂面，繼而再基于唯一的臨界滑裂面計算各個抗剪參數樣本對應的邊坡穩定性系數；因此，從理論上講LEM法所得到的臨界滑裂面與各個抗剪參數樣本對應的真實滑裂面存在差距。FEM 法所得的破壞模式與本文結果的失效模式2、失效模式3、失效模式4、失效模式5 相吻合，并未獲得失效模式1 和失效模式6，其主要原因在于：采用FEM 法進行邊坡穩定性計算時，通過構建同時滿足土體靜力平衡條件、應變相容條件以及本構關系的有限元模型，可獲得相較于LEM 法更為合理的結果；由于FEM 法在進行邊坡穩定性計算時需要考慮材料的本構關系，在已知土體本構關系的基礎上可獲得相較于上限法更為精確的結果，但在進行邊坡穩定性系數以及失效模式判斷時包含較多的人為因素。此外，在計算效率上相較與上限法大幅降低。上限法通過有限單元構建邊坡可靠度上限法數學規劃模型，然后采用隨機數學規劃的方法搜索得到邊坡的失穩區域，其結果更加接近于真實情況。計算結果表明：在多種隨機參數作用下，傳統的LEM 法會忽略邊坡的多種失效模式，可能會錯估邊坡的失效風險；FEM 法在已知材料本構關系的情況下可以獲得邊坡所有的失效模式，但土體本構關系復雜，當前并沒有統一的本構模型進行表征。本文上限法可以忽略材料的本構關系，獲得邊坡穩定性系數的上限解和破壞機構，其適用性和計算效率均較高。

需說明的是，本文同時計算了地下水隨機數的數量nw取30、50、70 三種情況。當nw=30時，會損失部分失效模式(失效模式1 和失效模式2)，計算精度較低；當nw=50時，可同時保證計算精度和計算效率；當nw=70時，計算精度較nw=50時雖略微提高，但計算量大幅增加、計算效率大幅降低。限于篇幅，本文僅給出了nw=50的計算結果。

4.7 邊坡單元失效概率分析

根據式(15)計算了50 個隨機地下水位作用下邊坡的整體失效概率。隨著地下水位的逐漸升高，邊坡整體失效概率從0.80%增加到92.70%，邊坡的安全性能大幅度降低。通過分段3 次埃爾米特插值得到邊坡整體失效概率與地下水位的關系曲線(如圖15 所示)。地下水位在5.5 m～11.5 m范圍內，整體失效概率變化幅度較小；地下水位在11.5 m～14 m 范圍內，整體失效概率變化幅度較大；地下水位在14 m～15 m 范圍內，整體失效概率變化幅度又變小。邊坡的整體失效概率隨地下水位的升高整體上呈現“平穩增加急劇上升平穩緩慢上升”的三段式特征。

圖15 邊坡整體失效概率與地下水位的關系Fig.15 The relationship between slope failure probability and groundwater

根據式(22)計算了50 個地下水位作用下邊坡的單元失效概率，圖16 為tw=1、tw=25、tw=50時的單元失效概率等值線；圖17 為特征點的單元失效概率與地下水位的關系曲線。對以上圖形進行分析，可以得到：

圖16 單元失效概率等值線Fig.16 Element failure probability contour

圖17 單元失效概率與地下水位的關系Fig.17 The relationship between the element failure probability and groundwater

1) 單元失效概率等值線可以非常直觀的看出土質邊坡每個部位的失效概率，地下水位tw=1、tw=25、tw=50時邊坡單元失效概率等值線與該水位下邊坡的失效模式輪廓線基本一致，但在失效可能性大小上存在明顯差距，其主要原因在于AP聚類分析根據邊坡中單元的失效信息將具有同一特征的失效面積聚為同一類，而單元失效概率等值線通過對已知單元的失效情況擬合得到。從理論上來說，如果邊坡只有一種失效模式，那么邊坡各單元的失效概率是相同的；如果邊坡具有多種失效模式，則邊坡各單元的失效概率會存在明顯差異，單元失效概率等值線圖也將存在差異；所有失效模式中重疊區域的單元的失效概率最大(如圖16 所示)。

2) 監測單元的單元失效概率與地下水位的關系如圖17 所示。在相同的地下水位作用下，不同的單元失效概率不同；在不同的地下水位作用下，同一單元的失效概率不同，不同位置處的單元失效概率隨著地下水位的變化情況也存在明顯差異。E2 單元和E3 單元隨著地下水位的升高失效概率呈現增加的趨勢，然而E1 單元和E4 單元隨著地下水位的升高失效概率呈現先增大后減少的趨勢，其主要原因在于所有失效模式均包含E2 單元，5 種失效模式(失效模式2、失效模式3、失效模式4、失效模式5、失效模式6)包含E1 單元，4 種失效模式(失效模式3、失效模式4、失效模式5、失效模式6)包含E3 單元，3 種失效模式(失效模式4、失效模式5、失效模式6)包含E4 單元，因此單元失效概率：E2≥E1≥E3≥E4。

5 結論

本文將有限元離散技術、上限法理論、相關非高斯隨機場模擬以及隨機規劃理論結合起來，提出了隨機地下水位作用下考慮參數空間變異性的邊坡可靠度分析方法，并編制了相應的計算程序，獲得了隨機地下水位作用下考慮參數空間變異性的邊坡穩定性系數、失效概率分布規律，為邊坡可靠度分析提供了一種新方法。其主要結論如下：

(1)傳統的LEM 法、FEM 法僅根據邊坡穩定性系數來進行邊坡穩定性的判斷，該方法僅反映了邊坡失效概率的程度。本文采用單元失效概率法對邊坡的失效概率進行研究，不僅可以反映邊坡失效概率的程度，而且可以根據單元的位置信息準確定位單元的失效情況。

(2)采用基于蒙特卡洛的迭代方法求解邊坡可靠度分析的隨機規劃模型時需要進行成千上萬次的模擬，本文基于Matlab 編制了高效的上限法并行程序，大大提高了計算效率。

(3)基于AP 聚類算法得到50 種水位下考慮土體參數變異性的邊坡6 種失效模式以及對應的失效概率，得到邊坡的系統失效概率為11.6180%，然而在多種隨機參數作用下，傳統的LEM 法會忽略邊坡的多種失效模式，可能會錯估邊坡的失效風險，這是因為LEM 法在計算邊坡臨界滑裂面時其計算原理為：首先通過蒙特卡洛方法生成計算所需隨機場，然后對每個隨機場的抗剪強度參數進行平均化取值，再通過計算得到邊坡的臨界滑裂面。因此，LEM 法并沒有得到邊坡所有的失效模式。而UBM 法通過尋求構建機動許可速度場的最小值得到邊坡的失效概率，所得結果更加符合真實情況。

(4)地下水位對考慮參數空間變異性的邊坡整體失效概率和單元失效概率影響均較大，其中整體失效概率隨著水位的升高不斷增大，邊坡單元失效概率等值線可便于巖土工程師更加直觀的了解邊坡各個位置的失效情況。