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

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

1 研究背景

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

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

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

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

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

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

另外，Kriging模型采用已知樣本x1，x2，…，xm的響應值 y1，y2，...，ym的線性組合來預測待測樣本點x的響應值，即

圖1 鏈表篩分法程序結構框圖

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

圖2 二分法程序結構框圖

2.3子集模擬法結構的失效概率可以表示為

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

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

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

在子集模擬中，一個很小的結構失效概率pf被認為是一系列較大的條件失效概率的乘積［7］

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

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

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

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

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

3 土石壩邊坡可靠度計算實例

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

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

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

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

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

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

圖5 不同條件樣本點響應值

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

圖6 不同條件樣本點的響應值預測效果

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

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

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

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

另外，如果采用直接蒙特卡羅方法計算邊坡失效概率，則只需判斷FS(x)＜1是否為真，而不用去計算FS(x)，所以需要的總迭代數(shù)只相當于驗算強度折減系數(shù)為1時的迭代數(shù)。如圖9（a）所示判斷第一層最小邊坡安全系數(shù)樣本所需要的迭代數(shù)是9，而利用順序搜索法計算其安全系數(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，而利用順序搜索法計算其安全系數(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ù)遠遠少于計算FS(x)所需要的迭代總數(shù)。但是考慮到本算例的失效概率在10-4這個數(shù)量級上，所以必須進行104數(shù)量級以上的直接蒙特卡羅模擬抽樣才能得到滿意結果。而子集模擬法在本算例的初始樣本抽取即第一層抽樣中只需要計算300個樣本點的FS值，在其余各分層抽樣需要各自計算270個樣本點的FS值。所以不采用Krig?ing代理模型的子集模擬法總共需要計算1110個樣本點的FS值。而采用Kriging代理模型后，只需要計算348個樣本點的FS值（T1其中48個樣本點為第二到四層中為提高Kriging預測精度所增設的樣本點）。所以圖10給出了方法1（T1直接蒙特卡羅法），方法2（T2子集模擬法+鏈表篩分法），方法3（T3子集模擬法+二分法+Krig?ing）和方法4（T4子集模擬法+鏈表篩分法+Kriging）的時間對比。從圖中顯示，方法4即本文所提方法所用時間最少，約為方法1的一半。如果遇到失效概率更小的算例時，則本文方法的優(yōu)點會更加明顯。

圖10 各方法的時間對比圖

4 結論

本文提出了一種土石壩邊坡小概率失效計算方法。該方法利用鏈表篩分算法自動計算出按邊坡安全系數(shù)大小升序排列的樣本集。這一樣本集合可以作為子集模擬法的第一層有序樣本集。然后第一層樣本集可以用來訓練Kriging代理模型。接著在子集模擬法的中間事件的樣本模擬中采用Kriging代理模型來預測條件樣本的邊坡安全系數(shù)，這樣可以避免對條件樣本點的彈塑性有限元計算。另外，該方法所采用的子集模擬部分引入了合理的中間失效事件將失效概率表達為一系列較大的條件失效概率乘積，因此只需要在子集模擬的每一層計算中抽取少量樣本點，就能使抽樣中心快速向失效域移動。對于本文所提供的10-4數(shù)量級的邊坡穩(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］.機械工程學報，2009，45（4）：33-38.

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

［4］宋述芳，呂震宙.高維小失效概率可靠性分析的序列重要抽樣法［J］.西北工業(yè)大學學報，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.