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

(1.昆明理工大學(xué)建筑工程學(xué)院，云南，昆明 650500；2.昆明理工大學(xué)電力工程學(xué)院，云南，昆明 650500；3.長(zhǎng)江勘測(cè)規(guī)劃設(shè)計(jì)研究有限責(zé)任公司，湖北，武漢 430010)

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

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

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

鑒于此，本文將有限元離散技術(shù)、上限法理論、相關(guān)非高斯隨機(jī)場(chǎng)模擬方法以及隨機(jī)規(guī)劃理論結(jié)合起來(lái)提出了隨機(jī)地下水位作用下考慮參數(shù)空間變異性的邊坡可靠度分析方法。

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

1.1 隨機(jī)地下水位的邊坡穩(wěn)定滲流場(chǎng)

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

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

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

式中：kx為土體材料x方向的滲透系數(shù)；ky為土體材料y方向的滲透系數(shù)；Hr為土體內(nèi)各點(diǎn)的隨機(jī)水頭函數(shù)。

1.2 描述參數(shù)空間變異性的隨機(jī)場(chǎng)模型

當(dāng)前普遍采用自相關(guān)函數(shù)表征土體的空間變異性，常用的自相關(guān)函數(shù)類型有指數(shù)型、高斯型、對(duì)數(shù)型、二階自回歸型、指數(shù)余弦型以及三角型[29]?？紤]邊坡的可靠性對(duì)自相關(guān)長(zhǎng)度較為敏感，對(duì)自相關(guān)函數(shù)的形式不敏感[30-31]，本文選用數(shù)學(xué)表達(dá)式更為簡(jiǎn)單的指數(shù)型自相關(guān)函數(shù)進(jìn)行研究，具體可采用式(2)表示：

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

由于土體參數(shù)之間存在一定的互相關(guān)性和土體參數(shù)本身存在一定的自相關(guān)性，因而涉及相關(guān)非高斯場(chǎng)的離散過(guò)程，本文采用類似文獻(xiàn)[32]的方法，采用基于喬列斯基分解的中點(diǎn)法進(jìn)行隨機(jī)場(chǎng)的生成，對(duì)于指數(shù)型自相關(guān)函數(shù)，土體黏聚力、內(nèi)摩擦角的相關(guān)對(duì)數(shù)隨機(jī)場(chǎng)可表示為：

式中：m=c,φ，c為土體黏聚力，φ為土體內(nèi)摩擦角；(xa,ya) ∈?；Hm(xa,ya)為相關(guān)非高斯隨機(jī)場(chǎng)；為相關(guān)標(biāo)準(zhǔn)高斯隨機(jī)場(chǎng)；σlnm=為對(duì)數(shù)正態(tài)分布m的均值， σm為對(duì)數(shù)正態(tài)分布m的標(biāo)準(zhǔn)差；ulnm為相應(yīng)正態(tài)變量lnm的均值；σlnm為相應(yīng)正態(tài)變量lnm的標(biāo)準(zhǔn)差。

1.3 邊坡可靠度分析的隨機(jī)規(guī)劃模型

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

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

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

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

有限單元公共邊的內(nèi)功功率為：

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

自重所作的外功功率為：

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

2 隨機(jī)規(guī)劃模型的求解策略

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

2)生成邊坡材料抗剪強(qiáng)度參數(shù)的隨機(jī)場(chǎng)。假設(shè)土體材料的黏聚力cr和內(nèi)摩擦角φr的自相關(guān)函數(shù)類型為指數(shù)型，采用喬列斯基分解的中點(diǎn)法進(jìn)行隨機(jī)場(chǎng)的生成：

本文采用容重超載的方式使邊坡達(dá)到極限狀態(tài)，隨機(jī)地下水位作用下考慮參數(shù)空間變異性時(shí)，邊坡的容重超載系數(shù)為：

式中：kγ(tm,tw)為第tw個(gè)地下水位作用下對(duì)應(yīng)著第tm個(gè)隨機(jī)場(chǎng)的邊坡的容重超載系數(shù)；、為土體處于極限狀態(tài)時(shí)的容重，其與以及相關(guān)； γa為土體實(shí)際容重；為強(qiáng)度折減以前第tm個(gè)非高斯隨機(jī)場(chǎng)中有限單元e形心處的黏聚力和內(nèi)摩擦角。

隨機(jī)地下水位作用下考慮參數(shù)空間變異性時(shí)，邊坡穩(wěn)定性系數(shù)如下：

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

3)計(jì)算邊坡穩(wěn)定滲流場(chǎng)。將步驟1)生成的地下水位的隨機(jī)數(shù)tw=(1,2,···,nw)，逐次代入邊坡穩(wěn)定滲流場(chǎng)的計(jì)算公式，進(jìn)而得到離散單元e三個(gè)節(jié)點(diǎn)的孔隙水壓力值：，e=(1,2,···,Ne)，tw=(1,2,···,nw)。

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

5)計(jì)算邊坡穩(wěn)定性系數(shù)并繪制相關(guān)曲線。

3 邊坡系統(tǒng)失效概率和可靠度指標(biāo)

邊坡工程中普遍采用邊坡穩(wěn)定性系數(shù)來(lái)評(píng)價(jià)邊坡的穩(wěn)定性。當(dāng)邊坡穩(wěn)定性系數(shù)km≥1，則認(rèn)為邊坡穩(wěn)定；邊坡穩(wěn)定性系數(shù)km<1，則認(rèn)為邊坡失穩(wěn)[34]。因此，邊坡的失效功能函數(shù)按式(14)進(jìn)行表示：

式中：tw=(1,2,···,nw)；tm=(1,2,···,nm)；I(tw,tm)為第tw個(gè)地下水位作用下對(duì)應(yīng)著第tm個(gè)隨機(jī)場(chǎng)的邊坡失效功能函數(shù)。

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

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

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

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

根據(jù)邊坡系統(tǒng)失效功能函數(shù)的定義，可以得到邊坡的系統(tǒng)失效概率估算公式為：

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

式中：Ie(tw,tm)為第tw個(gè)地下水位作用下對(duì)應(yīng)著第tm個(gè)隨機(jī)場(chǎng)中單元e的失效功能函數(shù)；uec(tw,tm)為第tw個(gè)地下水位作用下對(duì)應(yīng)著第tm個(gè)隨機(jī)場(chǎng)中單元e的速度。

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

圖3 AP 聚類算法選點(diǎn)過(guò)程Fig.3 AP clustering algorithm selection process

式中：

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

歸屬度迭代公式為：

式中：

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

定義單元失效概率為：

根據(jù)邊坡可靠度分析上限法數(shù)學(xué)規(guī)劃模型的求解策略中所得nw×nm個(gè)邊坡穩(wěn)定性系數(shù)km(tw,tm)，通過(guò)統(tǒng)計(jì)可以得到邊坡在第tw個(gè)地下水位作用下穩(wěn)定性系數(shù)的均值、標(biāo)準(zhǔn)差以及在所有可能發(fā)生的地下水位作用下邊坡穩(wěn)定性系數(shù)的均值、標(biāo)準(zhǔn)差。具體為：

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

4 算例分析

根據(jù)本文提出的方法，編制了相應(yīng)的上限法程序，對(duì)一個(gè)經(jīng)典的邊坡算例進(jìn)行了計(jì)算分析并與極限平衡法、有限元法計(jì)算結(jié)果對(duì)比，驗(yàn)證了本文計(jì)算方法的正確性。

4.1 均質(zhì)邊坡基本信息

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

表1 土體參數(shù)統(tǒng)計(jì)特性Table 1 Statistics of soil parameters

圖4 均質(zhì)邊坡計(jì)算模型 /mFig.4 Homogeneous slope calculation model

4.2 隨機(jī)地下水位和滲流場(chǎng)

圖5 地下水位隨機(jī)數(shù)分布直方圖Fig.5 Histogram of random distribution of groundwater level

圖6 邊坡穩(wěn)定滲流場(chǎng)Fig.6 Steady seepage field of slope

圖7 關(guān)鍵點(diǎn)處孔隙水壓力隨機(jī)變化情況Fig.7 The pore water pressure of the monitoring points

4.3 土體抗剪參數(shù)空間變異性和隨機(jī)場(chǎng)

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

圖8 邊坡抗剪強(qiáng)度參數(shù)隨機(jī)場(chǎng)Fig.8 Random Field of Slope Shear Strength Parameters

4.4 上限法與LEM 法、FEM 法的對(duì)比分析

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

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

為驗(yàn)證本文計(jì)算方法在地下水作用時(shí)的有效性，本文選取低、中、高三個(gè)地下水位進(jìn)行對(duì)比分析。表3 為本文計(jì)算方法與LEM 法、FEM 法計(jì)算結(jié)果對(duì)比情況，其中選取的水位分別為tw=1、tw=25、tw=50。圖9 為三種水位下邊坡穩(wěn)定性系數(shù)分布。

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

圖9 邊坡穩(wěn)定性系數(shù)分布特征Fig.9 Distribution characteristics of slope safety factor

計(jì)算結(jié)果表明：

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

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

4.5 邊坡穩(wěn)定性系數(shù)分布規(guī)律

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

圖10 邊坡穩(wěn)定性系數(shù)概率密度曲線Fig.10 Probability density curve of safety factor of slope

圖11 邊坡穩(wěn)定性系數(shù)累積概率密度曲線Fig.11 Cumulative probability density curve of safety factor of slope

1)邊坡穩(wěn)定性系數(shù)的分布與正態(tài)分布一致，且隨著地下水位的升高，邊坡穩(wěn)定性系數(shù)的均值逐漸降低，邊坡穩(wěn)定性系數(shù)的概率密度曲線和累積概率密度曲線均向左移動(dòng)。此外，隨著地下水位的升高，邊坡穩(wěn)定性系數(shù)的標(biāo)準(zhǔn)差逐漸降低，邊坡穩(wěn)定性系數(shù)的概率密度曲線分布范圍逐漸變窄，累積概率密度曲線逐漸變陡。

2)邊坡穩(wěn)定性系數(shù)的均值與地下水位直接相關(guān)，采用多項(xiàng)式擬合得到邊坡穩(wěn)定性系數(shù)的均值與地下水位的量化關(guān)系式為：

式(27)表明邊坡穩(wěn)定性系數(shù)的均值與地下水位成負(fù)相關(guān)。

3)對(duì)50 個(gè)地下水位在1000 個(gè)隨機(jī)場(chǎng)時(shí)所得的50 000 個(gè)邊坡穩(wěn)定性系數(shù)進(jìn)行統(tǒng)計(jì)分析，得到考慮土體參數(shù)空間變異性以及地下水位隨機(jī)性的邊坡穩(wěn)定性系數(shù)的分布特征，圖13 為隨機(jī)地下水位作用下邊坡穩(wěn)定性系數(shù)的頻率分布直方圖。隨機(jī)地下水位作用下邊坡穩(wěn)定性系數(shù)的均值為1.127，標(biāo)準(zhǔn)差為0.097。

圖13 邊坡穩(wěn)定性系數(shù)的頻數(shù)分布直方圖Fig.13 Frequency distribution histogram of safety factors

4.6 邊坡系統(tǒng)失效概率分析

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

表4 邊坡失效模式以及對(duì)應(yīng)的失效概率Table 4 Slope failure mode and the corresponding failure probability

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

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

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

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

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

4.7 邊坡單元失效概率分析

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

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

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

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

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

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

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

5 結(jié)論

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

(1)傳統(tǒng)的LEM 法、FEM 法僅根據(jù)邊坡穩(wěn)定性系數(shù)來(lái)進(jìn)行邊坡穩(wěn)定性的判斷，該方法僅反映了邊坡失效概率的程度。本文采用單元失效概率法對(duì)邊坡的失效概率進(jìn)行研究，不僅可以反映邊坡失效概率的程度，而且可以根據(jù)單元的位置信息準(zhǔn)確定位單元的失效情況。

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

(3)基于AP 聚類算法得到50 種水位下考慮土體參數(shù)變異性的邊坡6 種失效模式以及對(duì)應(yīng)的失效概率，得到邊坡的系統(tǒng)失效概率為11.6180%，然而在多種隨機(jī)參數(shù)作用下，傳統(tǒng)的LEM 法會(huì)忽略邊坡的多種失效模式，可能會(huì)錯(cuò)估邊坡的失效風(fēng)險(xiǎn)，這是因?yàn)長(zhǎng)EM 法在計(jì)算邊坡臨界滑裂面時(shí)其計(jì)算原理為：首先通過(guò)蒙特卡洛方法生成計(jì)算所需隨機(jī)場(chǎng)，然后對(duì)每個(gè)隨機(jī)場(chǎng)的抗剪強(qiáng)度參數(shù)進(jìn)行平均化取值，再通過(guò)計(jì)算得到邊坡的臨界滑裂面。因此，LEM 法并沒有得到邊坡所有的失效模式。而UBM 法通過(guò)尋求構(gòu)建機(jī)動(dòng)許可速度場(chǎng)的最小值得到邊坡的失效概率，所得結(jié)果更加符合真實(shí)情況。

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