基于克里金模型的邊坡穩(wěn)定可靠度分析方法

羅正東，董 輝，陳 鋮，蘇永華

（1. 湘潭大學(xué) 土木工程與力學(xué)學(xué)院，湖南 湘潭 411105；2. 湖南大學(xué) 巖土工程研究所，湖南 長沙 410082）

1 引 言

由于邊坡結(jié)構(gòu)力學(xué)參數(shù)及破壞模式的復(fù)雜性及不確定性，采用可靠度分析方法對(duì)邊坡的穩(wěn)定性進(jìn)行評(píng)估是目前較為理想的解決途徑，其研究近年來取得了較為理想的成果[1－2]，并發(fā)展了一些近似模型方法來對(duì)隱式功能函數(shù)進(jìn)行顯示化表達(dá)諸如響應(yīng)面法、神經(jīng)網(wǎng)絡(luò)法、支持向量機(jī)（SVM）等，但這些方法在使用過程中均存在一定的限制條件[3－5]。

源于地質(zhì)統(tǒng)計(jì)學(xué)的 Kriging模型是一種只需已知部分信息，而模擬某點(diǎn)待求信息的半?yún)?shù)化插值技術(shù)，與傳統(tǒng)代理模型相比，Kriging模型最主要的特征是其不需要給出模擬函數(shù)的形式，能避免多項(xiàng)式函數(shù)對(duì)結(jié)構(gòu)后續(xù)求解計(jì)算產(chǎn)生的影響，而且其使用非常靈活方便。2005年Kaymaz[6]首次將Kriging模型與一次二階矩分析法結(jié)合求解結(jié)構(gòu)的可靠度指標(biāo)。Won等[7]結(jié)合自適應(yīng)序列樣本方案與 Kriging模型，解決了基于可靠度的優(yōu)化設(shè)計(jì)問題。蘇永華等[8]利用Kriging代理模型，對(duì)邊坡的穩(wěn)定可靠度進(jìn)行了分析研究。

在 Kriging代理模型的擬合過程中必須已知一定數(shù)量的初始樣本點(diǎn)，抽取的訓(xùn)練樣本點(diǎn)代表性強(qiáng)弱將直接影響求解結(jié)果的準(zhǔn)確性。Kriging代理模型典型的試驗(yàn)設(shè)計(jì)方法是拉丁超立方設(shè)計(jì)[9]，但該方法的擬合精度依賴于樣本點(diǎn)的數(shù)量，樣本點(diǎn)越多其精度越高，計(jì)算量也隨之增大。主動(dòng)學(xué)習(xí)方法[10]能主動(dòng)搜索最佳訓(xùn)練樣本點(diǎn)，并將其逐漸加入已知樣本中，使樣本集不斷更新，進(jìn)而使 Kriging代理模型迅速達(dá)到預(yù)期的預(yù)測(cè)精度。

為充分考慮邊坡功能函數(shù)的復(fù)雜性及高度非線性的基本特征，本文在現(xiàn)有的研究基礎(chǔ)上通過引進(jìn)Kriging代理模型，采用具有主動(dòng)偵測(cè)能力的主動(dòng)學(xué)習(xí)和具有布滿整個(gè)空間、不含重復(fù)點(diǎn)的拉丁超立方抽樣進(jìn)行互補(bǔ)和聯(lián)合來抽取樣本點(diǎn)，建立了一種高效的邊坡穩(wěn)定可靠度求解方法。

2 邊坡極限狀態(tài)功能函數(shù)的特征

邊坡穩(wěn)定性分析已發(fā)展了多種方法和理論，但極限平衡分析法仍為當(dāng)前邊坡穩(wěn)定性分析的主要方法[11]。根據(jù)邊坡結(jié)構(gòu)的失穩(wěn)模式及其地質(zhì)特征，常用的極限平衡分析法大體可歸為兩類：一類是適合于圓弧形滑動(dòng)面的瑞典法、Bishop法等，另一類是適合于任意滑動(dòng)面的Janbu法、Spencer法、Sarma法等。在邊坡的穩(wěn)定可靠度研究中，需以極限平衡分析法為基礎(chǔ)展開分析，本文選取適合于任意滑動(dòng)面的Janbu法為例來展示基于Kriging代理模型的可靠度分析方法。

Janbu[12]提出其求解邊坡結(jié)構(gòu)安全系數(shù)的表達(dá)式為

顯見，邊坡穩(wěn)定極限狀態(tài)功能函數(shù)也屬于非線性隱函數(shù)，不便于利用傳統(tǒng)的分析方法來求解其可靠度。

3 隱式功能函數(shù)值的代理方程

3.1 Kriging隨機(jī)過程代理模型

源于地質(zhì)統(tǒng)計(jì)學(xué)的 Kriging代理模型[13]，是一種基于隨機(jī)過程的統(tǒng)計(jì)預(yù)測(cè)方法，以已知樣本信息的動(dòng)態(tài)構(gòu)造為基礎(chǔ)，充分考慮了變量在空間上的相關(guān)特性，建立近似函數(shù)模擬對(duì)象問題在某一點(diǎn)的未知信息，具有平滑效應(yīng)及估計(jì)方差最小的統(tǒng)計(jì)特征[6]。Kriging代理模型的響應(yīng)值與自變量間的對(duì)應(yīng)關(guān)系為式（4）由回歸部分及隨機(jī)過程兩部分組成。

式中：R為由R(η；S)構(gòu)成的對(duì)角元為1，大小為i×i的對(duì)稱矩陣；F為由i個(gè)樣本點(diǎn)處的回歸模型組成的i維向量；f (xnew)為回歸多項(xiàng)式，一般采用不高于二階的多項(xiàng)式，實(shí)際操作中由具體工程情況來確定；r(xnew)為訓(xùn)練樣本與待測(cè)點(diǎn)間的相關(guān)向量，其表達(dá)式為

3.2 Kriging隨機(jī)過程代理模型

4 拉丁超立方試驗(yàn)設(shè)計(jì)

構(gòu)造Kriging隨機(jī)過程代理模型需要功能函數(shù)Z的訓(xùn)練樣本，實(shí)際邊坡工程很難由試驗(yàn)獲得全部參數(shù)樣本值，因此需要通過試驗(yàn)設(shè)計(jì)抽樣法來抽取隨機(jī)變量的樣本數(shù)據(jù)。試驗(yàn)設(shè)計(jì)方法有多種，本文選取拉丁超立方試驗(yàn)設(shè)計(jì)[16]（LHS）抽取需要的樣本。這種方法是一種“充滿空間”的抽樣方法，即使試驗(yàn)樣本點(diǎn)相對(duì)均勻地散布于整個(gè)試驗(yàn)范圍之內(nèi)，并且每個(gè)水平只使用一次。

拉丁超立方試驗(yàn)設(shè)計(jì)作為一種改進(jìn)的蒙特卡洛抽樣方差縮減技術(shù)，其抽樣的均勻性得到了極大體現(xiàn)，該試驗(yàn)設(shè)計(jì)方法的基本原理為

輸出變量：

5 主動(dòng)搜索最可能失效的區(qū)域

5.1 主動(dòng)學(xué)習(xí)搜索規(guī)則

主動(dòng)學(xué)習(xí)法是由Lewis等[17]提出，學(xué)習(xí)過程中依據(jù)一定的標(biāo)識(shí)策略，從未標(biāo)識(shí)的數(shù)據(jù)樣本庫中選擇對(duì)分類器最為有利的樣本進(jìn)行進(jìn)一步訓(xùn)練，從而實(shí)現(xiàn)用最少的樣本點(diǎn)來實(shí)現(xiàn)盡可能高的分類精度。

主、被動(dòng)學(xué)習(xí)間的差異在于具體實(shí)施過程中主動(dòng)學(xué)習(xí)須和外界進(jìn)行信息交互，通過交互過程中反饋的信息來更新初始訓(xùn)練樣本，而被動(dòng)學(xué)習(xí)則依靠學(xué)習(xí)器自身，二者的差異見圖1。

圖1 主動(dòng)學(xué)習(xí)與被動(dòng)學(xué)習(xí)流程Fig.1 Flow of active learning and passive learning

5.2 主動(dòng)學(xué)習(xí)中學(xué)習(xí)函數(shù)的構(gòu)建

擬合好的 Kriging代理模型，需要另取樣本點(diǎn)來檢驗(yàn)擬合模型的精度，以確保其精確性和有效性，經(jīng)驗(yàn)證后的 Kriging代理模型才能用作代理模型進(jìn)行隱式函數(shù)的近似分析，常用檢驗(yàn)方法有多種[18]，本文選取較易操作的平均相對(duì)誤差法來進(jìn)行驗(yàn)證，其表達(dá)式為

5.3 邊坡可靠度的搜索法執(zhí)行步驟

融合Kriging模型、LHS方法及主動(dòng)學(xué)習(xí)搜索法的邊坡可靠度分析方法的實(shí)施步驟：

② 利用 LHS方法抽取少量樣本作為初始訓(xùn)練樣本 S，代入極限平衡法得到樣本的響應(yīng)值（本文采用Geo-Studio軟件中Slope/W模塊進(jìn)行求解，響應(yīng)值采用最小安全系數(shù)）。

③依據(jù)章節(jié)3.1～3.2方法，采用初始訓(xùn)練樣本及其響應(yīng)值建立邊坡功能函數(shù)的初始Kriging模型，根據(jù)式（13）對(duì)樣本點(diǎn)的功能函數(shù)值進(jìn)行預(yù)測(cè)，利用式（20）求解相對(duì)誤差值。

④主動(dòng)學(xué)習(xí)：搜索下一個(gè)最佳訓(xùn)練樣本點(diǎn)*x，不斷更新訓(xùn)練樣本，提高 Kriging模型對(duì)功能函數(shù)的預(yù)測(cè)精度。

⑤將*x代入極限平衡法得其真實(shí)響應(yīng)值，使訓(xùn)練樣本S得到更新，建立新的Kriging模型，不斷重復(fù)步驟④～⑤，直到滿足平均相對(duì)誤差小于0.01便可停止主動(dòng)學(xué)習(xí)。

⑥由 LHS方法及主動(dòng)學(xué)習(xí)搜索出的樣本構(gòu)建最優(yōu)的Kriging模型，最后采用JC法[1]求解邊坡的失效概率Pf。

6 算例分析

已知某土坡[19]剖面幾何形狀如圖2所示。土坡包含2個(gè)土層，其重度均為19 kN/m3。設(shè)土層的黏聚力c及內(nèi)摩擦角φ為相互獨(dú)立的正態(tài)隨機(jī)變量，其分布特征見表1，以Janbu法為基礎(chǔ)進(jìn)行求解。

圖2 邊坡剖面圖(單位：m)Fig.2 Geometry of slope (unit: m)

表1 土層物理力學(xué)性質(zhì)及統(tǒng)計(jì)特征Table 1 Physico-mechanical properties and statistical characteristics of soil layers

表2 算例求解結(jié)果對(duì)比Table 2 Comparison results for example

由表 2可知，當(dāng)以直接 Monte Carlo法計(jì)算2×105次所得結(jié)果為近似精確解時(shí)，本文方法計(jì)算得到的邊坡穩(wěn)定可靠度指標(biāo)相對(duì)誤差為 1.79%，計(jì)算精度高于傳統(tǒng)響應(yīng)面法的求解精度，基本能滿足工程實(shí)際要求，但邊坡穩(wěn)定極限平衡分析計(jì)算工作量與傳統(tǒng)響應(yīng)面法相當(dāng)，且遠(yuǎn)小于Monte Carlo計(jì)算法，尤其在復(fù)雜的隱式功能函數(shù)可靠度分析中將得到較好的體現(xiàn)。

7 結(jié) 語

以Janbu模型為例，利用Kriging隨機(jī)過程的統(tǒng)計(jì)預(yù)測(cè)方法，并結(jié)合拉丁超立方試驗(yàn)設(shè)計(jì)，導(dǎo)出了邊坡功能函數(shù)值的隨機(jī)過程代理模型。

引入主動(dòng)學(xué)習(xí)搜索規(guī)則及收斂準(zhǔn)則，提出了邊坡功能函數(shù)替代方程的更新迭代程序。建立了在最可能失效區(qū)域中依替代隨機(jī)過程方程的驗(yàn)算點(diǎn)法求解失效概率的方法，該方法在失效概率計(jì)算中不需要直接調(diào)用原功能函數(shù)，計(jì)算過程簡(jiǎn)單明了。

在詳細(xì)歸納出上述方法的操作執(zhí)行過程基礎(chǔ)之上，通過工程應(yīng)用分析，表明該方法求解過程簡(jiǎn)明，效率高，更具工程實(shí)用性。

[1] 蘇永華, 趙明華, 蔣德松, 等. 響應(yīng)面方法在邊坡穩(wěn)定可靠度分析中的應(yīng)用[J]. 巖石力學(xué)與工程學(xué)報(bào), 2006,25(7): 1417－1424.SU Yong-hua, ZHAO Ming-hua, JIANG De-song, et al.Application of response surface method to reliability analysis of slope stability[J]. Chinese Journal of Rock Mechanics and Engineering, 2006, 25(7): 1417－1424.

[2] 何婷婷, 尚岳全, 呂慶, 等. 邊坡可靠度分析的支持向量機(jī)法[J]. 巖土力學(xué), 2013, 34(11): 3269－3276.HE Ting-ting, SHANG Yue-quan, Lü Qing, et al. Slope reliability analysis using support vector machine[J]. Rock and Soil Mechanics, 2013, 34(11): 3269－3276.

[3] KRISHNAMURTHY T. Comparison of response surface construction methods for derivative estimation using moving least squares, Kriging and radial basis functions[R]. Austin: American Institute of Aeronautics and Astronautics, 2005.

[4] GOMES H M, AWRUCH A M. Comparison of response surface and neural network with other methods for structural reliability analysis[J]. Structural Safety, 2004,26(1): 49－57.

[5] CHERKASSKY V, MA Y. Practical selection of SVM parameters and noise estimation for SVM regression[J].Neural Networks, 2004, 17(1): 113－126.

[6] KAYMAZ I. Application of Kriging method to structural reliability problems[J]. Structural Safety, 2005, 27(2):133－151.

[7] WON J, CHOI C, CHOI J. Inproved dimesion reduction method in uncertainty analysis using Kriging interpolation[J]. Journal of Mechanical Science and Technology, 2009, 23(5): 1249－1260.

[8] 蘇永華, 羅正東, 張盼鳳, 等. 基于 Kriging的邊坡穩(wěn)定可靠度主動(dòng)搜索法[J]. 巖土工程學(xué)報(bào), 2013, 35(10):1863－1869.SU Yong-hua, LUO Zheng-dong, ZHANG Pan-feng, et al.Active searching algorithm for slope stability reliability based on Kriging model[J]. Chinese Journal of Geotechnical Engineering, 2013, 35(10): 1863－1869.

[9] SIMPSON T W, LIN D K J, CHEN W. Sampling strategies for computer experiments: Design and analysis[J]. International Journal of Reliability and Applications, 2002, 2(3): 209－240.

[10] GIUSEPPE R, DILEK H T. Active Learning: Theory and applications to automatic speech recognition[J]. IEEE Transactions on Speeeh and Audio Proeessing, 2005,13(4): 504－511.

[11] 孔令偉, 陳正漢. 特殊土與邊坡技術(shù)發(fā)展綜述[J]. 土木工程學(xué)報(bào), 2012, 45(5): 141－161.KONG Ling-wei, CHEN Zheng-han. Advancement in the techniques for special soils and slopes[J]. Chinese Civil Engineering Journal, 2012, 45(5): 141－161.

[12] JANBU N. Slope stability computations[C]//Proceedings of Embankment Dam Engineering. New York: John Wiley and Sons, 1973: 47－86.

[13] JONES D R, SCHONLAU M, WELCH W J. Efficient global optimization of expensive black-box functions[J].Journal of Global Optimization, 1998, 13(4): 445－492.

[14] SOREN N L, HANS B N, JACOB S. Aspects of the matlab toolbox DACE[R]. Denmark: Technical University of Denmark, 2002.

[15] ECHARD B, GAYTON N, LEMAIRE M. AK-MCS: An active learning reliability method combining Kriging and Monte Carlo simulation[J]. Structural Safety, 2011,33(2): 145－154.

[16] MCKAY M D, BECKMAN R J, CONOVER W J.Acomparis on of three methods for selecting values of input variables in the analysis of output from a computer code[J]. Technometrics, 1979, 21( 2): 239－245.

[17] LEWIS D, GALE W.A. A sequential algorithm for training text classifiers[C]//Proceedings of the 17th Annual International ACM SIGIR Conference on Research and Development in Information Retrieval.Dublin, Ireland: [s. n.], 1994: 3－12.

[18] 謝延敏. 基于Kriging模型和灰色關(guān)聯(lián)分析的板料成形工藝穩(wěn)健優(yōu)化設(shè)計(jì)研究[D]. 上海: 上海交通大學(xué),2007.

[19] HASSAN A M, WOLFF T F. Search algorithm for minimum reliability index of earth slopes[J]. Journal of Geotechnical and Geoenvironmental Engineering,1999, 125(4): 301－308.