土石壩邊坡小概率失效計(jì)算方法

杜建剛，徐佳成，羅顯楓，李 昕

（1.中國(guó)國(guó)際工程咨詢公司 能源業(yè)務(wù)部，北京 100048；2.大連理工大學(xué) 建設(shè)工程學(xué)部，遼寧 大連 116023 3.湖北理工大學(xué) 土木建筑工程學(xué)院，湖北 黃石 435003）

1 研究背景

在工程實(shí)際中經(jīng)常需要對(duì)一些小概率失效的土石壩進(jìn)行安全性評(píng)價(jià)分類，例如水庫(kù)土石壩的壩坡穩(wěn)定性安全評(píng)價(jià)（規(guī)范要求其失效概率數(shù)量級(jí)小于10-5）。面對(duì)此類情況時(shí)，如果采用直接蒙特卡羅隨機(jī)有限元模擬方法則必須抽取105數(shù)量級(jí)以上的樣本才能估算其失效概率，導(dǎo)致計(jì)算時(shí)間冗長(zhǎng)［1］。而子集模擬法是一種針對(duì)小概率失效問(wèn)題的分層抽樣算法，在每一層中引入合理的中間失效事件，將小概率問(wèn)題表示成一系列較大的條件概率的乘積［2-4］。而計(jì)算每一層的條件概率所需要抽取的樣本數(shù)量非常少，因此它計(jì)算失效概率所要的樣本總數(shù)量要遠(yuǎn)小于直接蒙特卡羅模擬方法需要的樣本數(shù)量。另一方面，很多商業(yè)巖土仿真軟件常采用基于土的彈塑性本構(gòu)的強(qiáng)度折減法計(jì)算邊坡安全系數(shù)，該方法本質(zhì)上是一種彈塑性增量有限元方法，而且在邊坡安全系數(shù)計(jì)算中還經(jīng)常使用基于二分搜索法的強(qiáng)度折減法來(lái)反復(fù)搜索邊坡安全系數(shù)的上下界限以至消耗機(jī)時(shí)過(guò)長(zhǎng)［5］，而鏈表篩分算法無(wú)需反復(fù)迭代搜索邊坡安全系數(shù)的上下界限，從而能快速計(jì)算多個(gè)樣本點(diǎn)的安全系數(shù)。因此可以在子集模擬法的第一層初始樣本點(diǎn)的響應(yīng)值計(jì)算中采用鏈表篩分算法來(lái)節(jié)省計(jì)算機(jī)時(shí)。而在子集模擬法的其余分層模擬中，可通過(guò)經(jīng)過(guò)初始樣本點(diǎn)訓(xùn)練后的Kriging代理模型來(lái)預(yù)測(cè)條件樣本點(diǎn)的響應(yīng)值，從而能進(jìn)一步節(jié)省機(jī)時(shí)。

本文根據(jù)以上分析，提出了基于鏈表篩分法和kringing代理模型的混合子集模擬法，并通過(guò)算例說(shuō)明本方法的有效性。

2 基于鏈表篩分法和Kriging代理模型的混合子集模擬法

2.1鏈表篩分法鏈表篩分法能自動(dòng)計(jì)算出按邊坡安全系數(shù)FS大小升序排列的樣本集，其算法原理框圖見(jiàn)圖 1 所示［5］。從圖中可看到土性參數(shù)黏聚力c=(c1，c2，...，cN)和內(nèi)摩擦角φ =(φ1，φ2，...，φN)是輸入宗量，是對(duì)邊坡土性參數(shù)c和φ分別進(jìn)行直接蒙特卡羅模擬時(shí)產(chǎn)生的樣本。當(dāng)邊坡土性參數(shù)的賦值為樣本對(duì)(cj，φj)時(shí)可計(jì)算得到相應(yīng)的邊坡安全系數(shù)FS(cj，φj)，(j=1，2，…，N)。圖1所示的邊坡安全系數(shù)初始順序搜索區(qū)間為FOS(l)=(M/LMAX)*l，(l=1，2，...)，LMAX（l代表加載的增量步序列號(hào)，LMAX代表最大增量步數(shù)）。然后利用彈塑性有限元增量加載的計(jì)算特性，在第l次增量步計(jì)算中篩選出達(dá)到最大迭代次數(shù)i=IMAX的樣本點(diǎn)（i代表迭代序列號(hào)），該樣本點(diǎn)對(duì)應(yīng)的邊坡安全系數(shù)就是FOS(l)。因此，鏈表篩分算法實(shí)際具有按邊坡安全系數(shù)從小到大增量順序搜索這一特性，所以輸入土性參數(shù)c和φ，可以自動(dòng)計(jì)算出按邊坡安全系數(shù)從小到大順序排列的安全系數(shù)FS。而如圖2所示商業(yè)軟件中采用的二分搜索法要進(jìn)行直接蒙特卡羅模擬時(shí)，只能按樣本產(chǎn)生的順序(cj，φj)來(lái)計(jì)算對(duì)應(yīng)的FS(cj，φj)(j=1，2，…，N)。另外，二分法常常需要把搜索上下界限設(shè)置為固定常數(shù)（圖2中t1=0，t2=5），這樣會(huì)導(dǎo)致當(dāng)邊坡的實(shí)際安全系數(shù)遠(yuǎn)小于二分法初始上界時(shí)，二分法搜索目標(biāo)安全系數(shù)耗時(shí)過(guò)長(zhǎng)。

2.2Kriging模型Kriging模型假定系統(tǒng)樣本x的響應(yīng)值為隨機(jī)過(guò)程函數(shù)y(x)，由回歸模型和隨機(jī)誤差組成［6］，即

式中；βj為回歸系數(shù)；fj(x)為基函數(shù)；z(x)為二階平穩(wěn)的隨機(jī)過(guò)程。

另外，Kriging模型采用已知樣本x1，x2，…，xm的響應(yīng)值 y1，y2，...，ym的線性組合來(lái)預(yù)測(cè)待測(cè)樣本點(diǎn)x的響應(yīng)值，即

圖1 鏈表篩分法程序結(jié)構(gòu)框圖

由于Kriging對(duì)y^(x)采用的是無(wú)偏估計(jì)，則需使其預(yù)測(cè)值的均方誤差E[(y^(x)-y(x))2]達(dá)到最小值。這樣便可求出系數(shù)c1，c2，...，cm及y^(x)。詳細(xì)推導(dǎo)步驟可參見(jiàn)相關(guān)文獻(xiàn)，在此不做贅述。

圖2 二分法程序結(jié)構(gòu)框圖

2.3子集模擬法結(jié)構(gòu)的失效概率可以表示為

式中： fX(x)=fX(x1，x2，…，xn)是隨機(jī)變量X=(X1，X2，…，Xn)T的聯(lián)合概率密度函數(shù)；y(x)是極限狀態(tài)函數(shù)；y(x)＜0代表結(jié)構(gòu)失效。

對(duì)于邊坡穩(wěn)定失效時(shí)，極限狀態(tài)函數(shù)y(x)可以表示為

式中：x一般是邊坡土層的土性參數(shù)；FS(x)是把x輸入邊坡穩(wěn)定分析程序后計(jì)算得到的邊坡安全系數(shù)。

在子集模擬中，一個(gè)很小的結(jié)構(gòu)失效概率pf被認(rèn)為是一系列較大的條件失效概率的乘積［7］

式中：Fi={ }FS＜fsi，(i=1，2，…，m)， 表示一系列中間失效事件；P(Fi)=P(FS＜fsi)表示中間失效事件Fi發(fā)生的概率；fsi表示中間失效事件Fi發(fā)生的閾值；這些閾值的大小順序依次為fs1＞fs2＞…fsm=1；這些中間失效事件滿足包含關(guān)系F1?F2?…?Fm=F；P(Fi|Fi-1)=P(FS＜fsi|FS＜fsi-1)表示已知中間失效事件Fi-1發(fā)生時(shí)，中間失效事件Fi發(fā)生的概率。由于P(Fi|Fi-1)的概率大于P(F)，所以需要的樣本數(shù)量也要少，另外在子集模擬中，用馬爾可夫蒙特卡羅法（Markov Chain Monte Carlo，MCMC）［8］計(jì)算P(Fi|Fi-1)所需要的樣本數(shù)量也遠(yuǎn)小于直接蒙特卡羅方法。

2.4基于混合子集模擬法的邊坡失效概率計(jì)算如研究背景中所提到的，可以利用鏈表篩分算法自動(dòng)計(jì)算出按邊坡安全系數(shù)大小升序排列的樣本集。這一樣本集合可以作為子集模擬法的第一層有序樣本集，并用該樣本集來(lái)訓(xùn)練Kriging模型。然后在子集模擬法的中間事件的樣本模擬中采用Kriging模型預(yù)測(cè)各個(gè)條件樣本點(diǎn)的邊坡安全系數(shù)。具體做法如下

（1）步驟1。用直接蒙特卡羅模擬方法產(chǎn)生N個(gè)獨(dú)立同分布的樣本

（2）步驟2。利用鏈表篩分法計(jì)算按升序排列的邊坡安全系數(shù)即樣本的響應(yīng)值，(j=1，2，…，N)。取其升序排列后的第Np0+1個(gè)值作為第一層中間失效事件發(fā)生的閾值fs1。

（6）步驟 6。重復(fù)步驟 4和 步驟 5直到進(jìn)入第m層失效域此時(shí)條件概率的估計(jì)值為是落在失效域Fm的樣本數(shù)量）。

3 土石壩邊坡可靠度計(jì)算實(shí)例

3.1計(jì)算實(shí)例本節(jié)用一個(gè)關(guān)于土石壩邊坡穩(wěn)定失效概率計(jì)算的例子來(lái)說(shuō)明所提方法的有效性。這一算例來(lái)自于文獻(xiàn)［9］。邊坡的幾何形狀及土層劃分如圖3所示。邊坡各土層的土性參數(shù)如表1所示。表中屬于隨機(jī)變量的土質(zhì)指標(biāo)均認(rèn)為服從對(duì)數(shù)正態(tài)分布。有限元計(jì)算采用八節(jié)點(diǎn)四邊形單元的平面應(yīng)變模型。為了和文獻(xiàn)［9］對(duì)比，這里的有限元模型也劃分為1 407個(gè)單元。約束條件是邊坡的底邊界固定約束，左右邊界水平約束。土體本構(gòu)采用理想彈塑性模型和摩爾-庫(kù)倫破壞準(zhǔn)則。

3.2本文方法計(jì)算結(jié)果圖4給出第一層全部樣本的頻率直方圖（本例中N=300），從頻率直方圖中能看出邊坡安全系數(shù)最小值約為1.2，邊坡安全系數(shù)最大值約為3.9。這說(shuō)明直接蒙特卡羅模擬法的抽樣中心遠(yuǎn)離失效域FS≤1。圖5（a）顯示出由子集模擬法第二層的第一個(gè)樣本產(chǎn)生的條件樣本點(diǎn)的響應(yīng)值因?yàn)樵诒纠屑俣╬0=0.1，所以馬爾可夫鏈的長(zhǎng)度L=9。從圖5（a）中可以看出由于滿足Metropolis-Hasting算法的接受準(zhǔn)則，初始狀態(tài)需要進(jìn)行更新，因此移動(dòng)到。而與相等是因?yàn)椴粷M足Metropolis-Hasting算法的接受準(zhǔn)則，所以即時(shí)狀態(tài)保持不變。圖5（b）顯示出由子集模擬法第四層的第一個(gè)樣本產(chǎn)生的條件樣本點(diǎn)的響應(yīng)值從圖5

（b）中可以看出這些樣本點(diǎn)全部落在失效域FS≤1。這也說(shuō)明子集模擬法通過(guò)分層模擬條件樣本點(diǎn)，容易抽取到落入FS≤1失效域內(nèi)的樣本點(diǎn)。

圖3 邊坡的幾何形狀及土層分區(qū) （單位：m）

圖4 首層全部樣本響應(yīng)值FS(x(1))的頻率直方圖

表1 各土層的土性參數(shù)

圖5 不同條件樣本點(diǎn)響應(yīng)值

為了說(shuō)明Kriging代理模型在分層模擬中的預(yù)測(cè)效果，圖6（a）給出了Kriging代理模型對(duì)條件樣本點(diǎn)的響應(yīng)值的預(yù)測(cè)值與有限元法算出的邊坡安全系數(shù)（即目標(biāo)值）的對(duì)比。圖中9個(gè)圓圈的中心位置非?？拷c坐標(biāo)軸成45°夾角的直線（其中有3個(gè)圓圈的中心是重合的），說(shuō)明Kriging代理模型對(duì)分層模擬的第二層條件樣本點(diǎn)的響應(yīng)值預(yù)測(cè)效果顯著。圖6（b）給出了Kriging代理模型對(duì)條件樣本點(diǎn)的響應(yīng)值的預(yù)測(cè)值與有限元法算出的邊坡安全系數(shù)（即目標(biāo)值）的對(duì)比。圖中9個(gè)圓圈的中心位置非?？拷c坐標(biāo)軸成45°夾角的直線（其中有4個(gè)圓圈的中心是重合的），說(shuō)明Kriging代理模型對(duì)分層模擬的第四層條件樣本點(diǎn)的響應(yīng)值預(yù)測(cè)效果同樣顯著。

圖6 不同條件樣本點(diǎn)的響應(yīng)值預(yù)測(cè)效果

圖7 順序搜索FS和二分搜索FS

圖8 FS順序搜索和二分搜索FS

3.3 本文方法與其它方法比較圖7（a）顯示了鏈表篩分法計(jì)算FS()時(shí)在彈塑性有限元分析的每一增量步中所需要的迭代數(shù)。從圖中可以看出當(dāng)強(qiáng)度折減系數(shù)接近目標(biāo)邊坡安全系數(shù)時(shí)，則對(duì)應(yīng)該強(qiáng)度折減系數(shù)的增量步中所需要的迭代數(shù)也迅速增加。如圖所示，順序搜索的總迭代數(shù)為2+2+2+2+2+2+2+4+9+9+25+35+67+135+216+500=1014。圖7（b）顯示了二分法計(jì)算FS)時(shí)對(duì)每一次的強(qiáng)度折減系數(shù)計(jì)算時(shí)所花的迭代數(shù)。從圖中可以看到二分搜索的總迭代數(shù)為500+500+2+16+235+500+400+500=2653。圖8（a）顯示了鏈表篩分法計(jì)算FS()時(shí)在彈塑性有限元分析的每一增量步中所需要的迭代數(shù)。從圖中能看出當(dāng)強(qiáng)度折減系數(shù)接近目標(biāo)邊坡安全系數(shù)時(shí)，則對(duì)應(yīng)該強(qiáng)度折減系數(shù)的增量步中所需要的迭代數(shù)也迅速增加。如圖所示，總迭代數(shù)為2+2+2+2+2+9+15+18+26+43+117+295+500=1033。圖8（b）顯示了二分法計(jì)算FS()時(shí)對(duì)每一次的強(qiáng)度折減系數(shù)計(jì)算時(shí)所花的迭代數(shù)。從圖中可以看到總迭代數(shù)為500+500+17+500+500+35+300+500=2852。

圖9 順序搜索min()和順序搜索max()

另外，如果采用直接蒙特卡羅方法計(jì)算邊坡失效概率，則只需判斷FS(x)＜1是否為真，而不用去計(jì)算FS(x)，所以需要的總迭代數(shù)只相當(dāng)于驗(yàn)算強(qiáng)度折減系數(shù)為1時(shí)的迭代數(shù)。如圖9（a）所示判斷第一層最小邊坡安全系數(shù)樣本所需要的迭代數(shù)是9，而利用順序搜索法計(jì)算其安全系數(shù)所經(jīng)歷的迭代總數(shù)為2+2+2+2+2+2+2+2+9+9+14+17+21+28+43+89+154+252+500=1152。如圖9（b）所示判斷第一層最大邊坡安全系數(shù)樣本所需要的迭代數(shù)是2，而利用順序搜索法計(jì)算其安全系數(shù)所經(jīng)歷的迭代總數(shù)為2*30+6+13+18+21+28+37+41+56+60+68+145+226+304+367+433+500=2383。因此，判斷FS(x)＜1是否為真所需要的迭代數(shù)遠(yuǎn)遠(yuǎn)少于計(jì)算FS(x)所需要的迭代總數(shù)。但是考慮到本算例的失效概率在10-4這個(gè)數(shù)量級(jí)上，所以必須進(jìn)行104數(shù)量級(jí)以上的直接蒙特卡羅模擬抽樣才能得到滿意結(jié)果。而子集模擬法在本算例的初始樣本抽取即第一層抽樣中只需要計(jì)算300個(gè)樣本點(diǎn)的FS值，在其余各分層抽樣需要各自計(jì)算270個(gè)樣本點(diǎn)的FS值。所以不采用Krig?ing代理模型的子集模擬法總共需要計(jì)算1110個(gè)樣本點(diǎn)的FS值。而采用Kriging代理模型后，只需要計(jì)算348個(gè)樣本點(diǎn)的FS值（T1其中48個(gè)樣本點(diǎn)為第二到四層中為提高Kriging預(yù)測(cè)精度所增設(shè)的樣本點(diǎn)）。所以圖10給出了方法1（T1直接蒙特卡羅法），方法2（T2子集模擬法+鏈表篩分法），方法3（T3子集模擬法+二分法+Krig?ing）和方法4（T4子集模擬法+鏈表篩分法+Kriging）的時(shí)間對(duì)比。從圖中顯示，方法4即本文所提方法所用時(shí)間最少，約為方法1的一半。如果遇到失效概率更小的算例時(shí)，則本文方法的優(yōu)點(diǎn)會(huì)更加明顯。

圖10 各方法的時(shí)間對(duì)比圖

4 結(jié)論

本文提出了一種土石壩邊坡小概率失效計(jì)算方法。該方法利用鏈表篩分算法自動(dòng)計(jì)算出按邊坡安全系數(shù)大小升序排列的樣本集。這一樣本集合可以作為子集模擬法的第一層有序樣本集。然后第一層樣本集可以用來(lái)訓(xùn)練Kriging代理模型。接著在子集模擬法的中間事件的樣本模擬中采用Kriging代理模型來(lái)預(yù)測(cè)條件樣本的邊坡安全系數(shù)，這樣可以避免對(duì)條件樣本點(diǎn)的彈塑性有限元計(jì)算。另外，該方法所采用的子集模擬部分引入了合理的中間失效事件將失效概率表達(dá)為一系列較大的條件失效概率乘積，因此只需要在子集模擬的每一層計(jì)算中抽取少量樣本點(diǎn)，就能使抽樣中心快速向失效域移動(dòng)。對(duì)于本文所提供的10-4數(shù)量級(jí)的邊坡穩(wěn)定的失效概率求解范例中，只需要進(jìn)行子集模擬的四層抽樣計(jì)算，就能和直接蒙特卡羅模擬方法的結(jié)果非常接近。這足以說(shuō)明本文所提方法的有效性。此外，通過(guò)本文方法與直接蒙特卡羅方法的計(jì)算時(shí)間分析比較中，可以看出本文方法更適合計(jì)算邊坡小概率失效問(wèn)題。

［1］LI D，XIAO T，Cao Z，et al.Enhancement of random finite element method in reliability analysis and risk assess ment of soil slopes using Subset Simulation［J］.Landslides，2016，13（2）：293-303.

［2］宋述芳，呂震宙.基于馬爾可夫蒙特卡羅子集模擬的可靠性靈敏度分析方法［J］.機(jī)械工程學(xué)報(bào)，2009，45（4）：33-38.

［3］宋述芳，呂震宙.基于子集模擬和重要抽樣的可靠性靈敏度分析方法［J］.力學(xué)學(xué)報(bào)，2008，40（5）：654-662.

［4］宋述芳，呂震宙.高維小失效概率可靠性分析的序列重要抽樣法［J］.西北工業(yè)大學(xué)學(xué)報(bào)，2006，24（6）：782-784.

［5］LUO X，CHENG T，LI X，et al.Slope safety factor search strategy for multiple sample points for reliability analy?sis［J］.Engineering Geology，2012，129：27-37.

［6］LUO X，LI X，ZHOU J，et al.A Kriging-based hybrid optimization algorithm for slope reliability analysis［J］.Structural Safety，2012，34（1）：401-406.

［7］AU S，BECK J L.Estimation of small failure probabilities in high dimensions by subset simulation［J］.Probabi?listic Engineering Mechanics，2001，16（4）：263-277.

［8］LI H，CAO Z.Matlab codes of Subset Simulation for reliability analysis and structural optimization［J］.Structur?al and Multidisciplinary Optimization，2016，54（2）：391-410.

［9］CHO S E.Probabilistic stability analyses of slopes using the ANN-based response surface［J］.Computers and Geotechnics，2009，36（5）：787-797.