基于克里金代理模型與子集模擬的邊坡高效可靠度分析

雷世平，李京澤，劉磊磊*，李云青

（1.湖南省有色地質(zhì)勘查局一總隊，湖南 郴州 423099；2.有色金屬成礦預測與地質(zhì)環(huán)境監(jiān)測教育部重點實驗室（中南大學），湖南 長沙 410083；3.湖南省有色資源與地質(zhì)災害探查湖南省重點實驗室，湖南 長沙 410083；4.中南大學地球科學與信息物理學院，湖南 長沙 410083）

0 引言

近年來，基于概率論與統(tǒng)計學的可靠度分析方法在邊坡穩(wěn)定性分析中得到了廣泛的應用［1-6］?？煽慷确治龇椒▽r土體中客觀存在的不確定性納入邊坡穩(wěn)定性分析框架當中，提高了對邊坡安全性與可靠性評估的準確性與全面性［7］。傳統(tǒng)可靠度分析方法一般通過蒙特卡洛模擬對邊坡失效概率進行計算，但該方法由于需要反復進行大量的邊坡穩(wěn)定性分析，從而存在計算效率低下、計算資源占用過大的問題。針對該問題，代理模型作為提高計算效率、節(jié)省計算資源的方法逐漸受到了關注，該類方法一般采用簡單的顯式函數(shù)或模型代替復雜的邊坡穩(wěn)定性分析，避免了蒙特卡洛模擬中復雜且高度非線性的邊坡穩(wěn)定性分析模型的重復計算，從而極大地提高了可靠度分析的計算效率［8］。

目前常用的代理模型方法有傳統(tǒng)響應面法［9］、多重響應面法［10］和克里金法［11］等。其中，克里金代理模型由于其在計算靈活性和準確性方面的優(yōu)勢而在可靠度分析方法中得到了廣泛的應用［12-14］。例如，GASPAR等［14］提出該方法不僅可以提供未采樣點的估計值，而且可對未采樣點的估計誤差進行量化，進而推動其在邊坡可靠性分析與不確定性量化中的應用。LIU等［15］提出了一種自適應克里金模型，通過在大多數(shù)不確定點逐一選取訓練樣本，達到用最少的樣本擬合邊坡穩(wěn)定性分析的目的，克服了其它模型隨機選取訓練樣本而計算耗時的不足。同時，為了提高克里金模型在邊坡可靠度分析中的適用性與準確性，多種全局優(yōu)化算法如人工蜂群算法［16］、遺傳算法［11］和粒子群優(yōu)化算法［17］近年來亦被用于優(yōu)化克里金模型參數(shù)?？傊?，雖然克里金方法可以有效地用于邊坡可靠度分析，但當邊坡處于小失效概率水平（如失效概率＜10-7）時，傳統(tǒng)蒙特卡洛模擬處理高維矩陣時可能導致計算機內(nèi)存溢出，造成可靠度計算求解困難，從而限制了克里金模型在邊坡可靠度分析中的適用性。

為解決上述問題，本文在建立邊坡穩(wěn)定性克里金代理模型的基礎上，提出結(jié)合子集模擬方法開展高效的邊坡可靠度分析。在該方法中，子集模擬被用于替代傳統(tǒng)克里金可靠度分析中的蒙特卡洛模擬，使得計算矩陣維度大大降低，從而提高計算效率。

1 高效邊坡可靠度分析方法

克里金方法是一種利用鄰近已知位置觀測數(shù)據(jù)進行空間最優(yōu)線性預測的方法，可對區(qū)域化變量求最優(yōu)、無偏內(nèi)插估計值［18-19］。該方法通過假設已知點與未知點之間的距離或方向來反映已知點和未知點之間的空間相關性，并通過高斯過程建模對于插值點進行預測。在邊坡可靠度分析中，克里金法以土體參數(shù)（如抗剪強度）及邊坡響應（如安全系數(shù)）為輸入變量進行參數(shù)建模，隨后將所建立的克里金模型作為代理模型取代耗費計算資源與計算時間的傳統(tǒng)確定性分析模型，進而估計邊坡失效概率。本文中，克里金模型的建立是通過基于MATLAB的開源代碼DACE工具箱進行實現(xiàn)［20］。

為降低輸入變量的矩陣維度以解決小失效概率邊坡的可靠度分析與計算效率問題，本文使用子集模擬代替直接蒙特卡洛模擬進行邊坡失效概率的計算。子集模擬作為一種高效蒙特卡洛模擬方法通過一系列具有較高概率的中間失效事件的概率乘積來表達一個小概率的失效事件F。在邊坡穩(wěn)定性分析中，邊坡的安全系數(shù)為臨界滑面（即所有潛在滑面中安全系數(shù)最小的滑面）的安全系數(shù)。由此，邊坡的功能函數(shù)如下：

式中：x——隨機變量，如土體抗剪強度（粘聚力c、內(nèi)摩擦角φ）；FSmin(x)——邊坡臨界滑面的安全系數(shù)；G(x)——邊坡功能函數(shù)，其反映了邊坡的穩(wěn)定狀態(tài)，當G(x)≥0時，邊坡保持穩(wěn)定狀態(tài)，反之當G(x)＜0時，則認為邊坡失穩(wěn)。

在邊坡可靠度分析當中，邊坡的失效概率Pf表述如下：

式 中：P(·)——事 件 的 概 率；Fi={G(x)＜gi(x)，i=1，2，…，m}——一組中間失效事件，這些事件由閾值FS值的遞減序列g1(x)＞g2(x)＞…＞gm(x)定義；m——子集模擬的層數(shù)。

在計算過程中，gi(x)通過對前一個中間失效事件Fi-1=G(x)＜gi-1(x)的條件樣本的統(tǒng)計分析來確定。其中，臨界閾值gi(x)使得中間失效事件所對應的失效概率均為設定的特定的條件概率值p0。據(jù)LIU等［21］與LI等［22］，子集模擬中的初始蒙特卡洛樣本數(shù)N0=500與條件概率p0=0.1可同時兼顧計算精度與計算效率，因此本文子集模擬中均按照此值進行設定。通過以上方法計算可得邊坡失效概率，隨后可參考《建筑結(jié)構(gòu)可靠性設計統(tǒng)一標準》（GB 50068—2018）［23］（表1）根據(jù)邊坡安全等級進一步確定邊坡的加固與處理方法。

表1 建筑結(jié)構(gòu)安全等級及相應可靠度指標與失效概率[21]Table 1 Building structure safety grade and corresponding reliability index and failure probability

2 計算流程

本文提出的邊坡可靠度分析的計算流程如圖1所示，主要步驟如下：

圖1 基于子集模擬與克里金代理模型的邊坡可靠度分析流程Fig.1 Flow chart of slope reliability analysis based on Subset simulation and the Kriging metamodel

（1）建立邊坡穩(wěn)定性分析模型：基于邊坡幾何參數(shù)和巖土體的物理力學參數(shù)建立邊坡穩(wěn)定性分析模型，可采用有限元法或極限平衡法進行建模。本文采用極限平衡分析軟件SLOPE/W進行建模［24］。

（2）生成訓練樣本：基于土體參數(shù)統(tǒng)計特性，采用拉丁超立方抽樣生成一定數(shù)量的訓練樣本。SILVESTRINI等［25］指 出 訓 練 樣 本 數(shù) 一 般 為10D～15D，其中D為隨機變量個數(shù)［26］。

（3）計算樣本響應值：將訓練樣本輸入至步驟（1）邊坡模型中，計算得到訓練樣本對應的邊坡安全系數(shù)。

（4）構(gòu)建訓練樣本數(shù)據(jù)集：將步驟（2）訓練樣本與步驟（3）計算得到的邊坡安全系數(shù)進行組合，組成訓練樣本數(shù)據(jù)集。

（5）構(gòu)建克里金代理模型：將訓練樣本數(shù)據(jù)集輸入至克里金模型當中，計算模型中的未知參數(shù)，構(gòu)建克里金模型。

（6）可靠度計算：基于步驟（5）建立的克里金代理模型，進行子集模擬，計算不同中間事件的條件概率，并利用公式（2）計算邊坡失效概率。

3 算例

本節(jié)通過2個邊坡算例說明所提出的方法的有效性。將直接蒙特卡洛模擬方法結(jié)果與本文方法結(jié)果進行對比，以驗證本文所提方法對邊坡可靠度計算的有效性。同時，探討模型對訓練樣本數(shù)量和回歸模型選擇的敏感性。

3.1 算例I：單層粘性土質(zhì)邊坡

算例I為一單層粘性土質(zhì)邊坡，該邊坡幾何形狀如圖2所示，坡高為10 m，坡比為1∶1。土體重度γ=20 kN/m3，粘聚力c和摩擦角φ為服從對數(shù)正態(tài)分布的互相關隨機變量，參數(shù)統(tǒng)計特征如表2所示。采用SLOPE/W模塊中簡化Bishop法進行極限平衡條分法邊坡穩(wěn)定性分析，模型共定義潛在滑面5985條?；谕馏w參數(shù)的均值，邊坡的臨界滑面如圖2中所示，其安全系數(shù)為1.207，與文獻中基本一致［11］。

圖2 算例I邊坡幾何結(jié)構(gòu)與穩(wěn)定性分析結(jié)果Fig.2 Slope geometry and slope stability analysis result for Example I

表2 算例I土體參數(shù)統(tǒng)計特性Table 2 Statistics of the soil parameters of Example I

在上述邊坡確定性分析模型基礎上，下面通過基于克里金的子集模擬法對該邊坡開展可靠度分析。首先，根據(jù)10倍隨機變量準則，利用拉丁超立方抽樣技術生成20組標準正態(tài)隨機樣本，并將這些標準正態(tài)隨機樣本按照表2所示土體參數(shù)統(tǒng)計特性轉(zhuǎn)換為對數(shù)正態(tài)隨機樣本，將其輸入圖2所示邊坡穩(wěn)定性分析模型，獲得相應的邊坡穩(wěn)定性響應值，用于訓練克里金模型。構(gòu)建克里金代理模型時，回歸模型選擇二次多項式函數(shù)，相關函數(shù)模型選擇高斯型。隨后，將訓練完畢的克里金代理模型接入子集模擬計算邊坡失效概率，結(jié)果可得該邊坡的失效概率為0.0530。

為驗證本文方法的準確性與計算效率，將本文提出的方法與直接蒙特卡洛模擬、基于克里金的蒙特卡洛模擬的邊坡可靠度計算結(jié)果進行對比，計算結(jié)果見表3。由表3可知，基于克里金模型的蒙特卡洛法與子集模擬方法均可對邊坡失效概率開展高效分析，計算結(jié)果誤差相較于蒙特卡洛模擬法均小于10%。同時，基于克里金模型的蒙特卡洛法與子集模擬方法計算時間相較于傳統(tǒng)蒙特卡洛模擬法大幅度地降低。但需要指出的是，在基于克里金的蒙特卡洛模擬法中計算樣本數(shù)設置為1×106，但若追求更高的計算精度，如當蒙特卡洛樣本數(shù)達到或超過108時，模型內(nèi)存占用將超出計算設備的內(nèi)存容量（16 GB），從而導致內(nèi)存溢出而無法進行計算。間接反應基于克里金的蒙特卡洛模擬在針對小概率失效事件的可靠度分析中存在一定的局限性。與之相比，基于克里金的子集模擬法在相近誤差水平與計算時間的前提下，借助中間失效事件使得計算樣本數(shù)僅1×103，大大降低了計算矩陣的維度，從而有效解決了基于克里金的蒙特卡洛模擬方法內(nèi)存占用過大的問題。

表3 不同分析方法計算的算例I邊坡失效概率結(jié)果Table 3 Probability analysis results obtained by different methods for Example I

進一步地，為探究本文所提方法中訓練樣本數(shù)量對可靠度分析結(jié)果的影響，此處將不同訓練樣本數(shù)量下的計算結(jié)果與蒙特卡洛模擬結(jié)果進行對比。由圖3可知，隨著訓練樣本數(shù)量的增加，基于克里金的子集模擬的失效概率計算結(jié)果基本與蒙特卡洛模擬結(jié)果保持一致，同時計算結(jié)果的波動范圍較小，說明采用10倍隨機變量的樣本數(shù)量已基本可滿足該方法中訓練模型精度需求。

圖3 算例I中失效概率與訓練樣本數(shù)關系Fig.3 Relationship between probability of failure and the number of training samples for Example I

為探究不同相關模型與回歸模型對失效概率計算結(jié)果的影響，此處將不同模型的計算結(jié)果與蒙特卡洛模擬結(jié)果進行對比，計算結(jié)果更接近蒙特卡洛模擬則認為其計算精度更高。其中，訓練樣本數(shù)量均選擇為20組。由表4可知，相較于常數(shù)型與一次函數(shù)回歸模型，選取二次多項式的克里金模型對于失效概率的計算精度更高。在一次函數(shù)與二次多項式回歸模型下的情況，不同相關模型對失效概率的計算結(jié)果影響較小。而在常數(shù)型回歸模型下，高斯型與三次樣條型相關模型的計算精度要顯著高于其他3種相關模型?？偟膩碚f，上述結(jié)果表明二次型多項式回歸模型與高斯型相關模型構(gòu)建的克里金模型更適用于該邊坡穩(wěn)定性的失效概率計算。

表4 算例I不同模型的邊坡失效概率結(jié)果的比較Table 4 Probability analysis results obtained by different models for Example I

3.2 算例II：工程實例邊坡

算例II為一工程實例邊坡，該斜坡位于中國香港西貢區(qū)北部，由距坡頂9 m深的崩積層土體組成［27］。邊坡根據(jù)不同傾角分為5個區(qū)段，邊坡中部存有上覆物所致30 kN點荷載，地下水水位情況及相關工程地質(zhì)情況如圖4所示。土體重度γ=19 kN/m3，粘聚力c和摩擦角φ為服從對數(shù)正態(tài)分布的互相關隨機變量，參數(shù)統(tǒng)計特征如表5所示。采用SLOPE/W模塊中Morgenstern-Price法進行極限平衡條分法邊坡穩(wěn)定性分析，模型共定義潛在滑面5355條?；谕馏w參數(shù)的均值，邊坡的臨界滑面如圖4中所示，其安全系數(shù)為1.330。

圖4 算例II邊坡幾何結(jié)構(gòu)與穩(wěn)定性分析結(jié)果Fig.4 Slope geometry and slope stability analysis result for Example II

表5 算例II土體參數(shù)統(tǒng)計特性Table 5 Statistics of the soil parameters of Example II

通過基于克里金的子集模擬法對該邊坡開展可靠度分析。首先通過拉丁超立方抽樣生成一定數(shù)量（比如20）的訓練樣本，并輸入至圖5所示的邊坡穩(wěn)定性分析模型求解相應的安全系數(shù)值，進而組成一定數(shù)量的訓練樣本集，以此構(gòu)建克里金代理模型。模型構(gòu)建過程中，回歸模型選擇二次多項式函數(shù)，相關函數(shù)模型選擇高斯型。然后，基于構(gòu)建的克里金代理模型，開展子集模擬分析，計算得到該邊坡的失效概率為0.0186。隨后，將本文方法與直接蒙特卡洛模擬法以及基于克里金的蒙特卡洛法進行對比分析，結(jié)果如表6所示。由表6可得，基于克里金的子集模擬法在計算時間上大大低于蒙特卡洛模擬，同時計算精度也高于基于克里金的蒙特卡洛法。此外，在計算精度相近的情況下，基于克里金的子集模擬法（1×103）的計算樣本數(shù)較基于克里金的蒙特卡洛法（1×106）降低了3個數(shù)量級，極大地降低了計算過程中計算資源的消耗，進一步證明了基于克里金的子集模擬法的計算效率更高。

圖5 算例II中失效概率與訓練樣本數(shù)關系Fig.5 Relationship between probability of failure and the number of training samples for Example II

表6 不同分析方法計算的算例II邊坡失效概率結(jié)果Table 6 Probability analysis results obtained by different methods for Example II

接下來，為了探究本文所提方法中訓練樣本數(shù)對可靠度分析結(jié)果的影響，此處將不同訓練樣本數(shù)量下的計算結(jié)果與蒙特卡洛模擬計算結(jié)果進行對比。由圖5可知，隨著訓練樣本數(shù)量的增加，基于克里金的子集模擬對于邊坡失效概率圍繞著蒙特卡洛模擬結(jié)果波動，但波動范圍較小且基本與蒙特卡洛模擬法結(jié)果一致。與算例Ⅰ相同，采用10倍隨機變量數(shù)的樣本數(shù)量基本可滿足該方法中構(gòu)建模型的精度需求，并且隨著訓練樣本的增加，邊坡失效概率變化不顯著。

4 結(jié)論

本文提出了一種基于克里金模型與子集模擬的邊坡可靠度分析方法，并研究了不同回歸模型和相關函數(shù)模型以及訓練樣本數(shù)對該方法精度的影響。該方法利用子集模擬估算邊坡可靠度，可以有效用于高維小失效概率可靠度分析，提高邊坡可靠度分析的計算效率與適用性。采用一個單層粘性土邊坡與一個堅硬地層上的雙層飽和不排水粘性土邊坡驗證了所提方法的有效性。得到的主要結(jié)論如下：

（1）基于克里金模型與子集模擬的邊坡可靠度分析方法能夠有效地替代傳統(tǒng)蒙特卡洛模擬法，準確地估計邊坡的失效概率。基于克里金的子集模擬法計算效率相較于傳統(tǒng)蒙特卡洛方法顯著提高。同時，該方法的計算精度可滿足可靠度分析的需求。

（2）所提方法可有效解決邊坡可靠度分析中的小概率失效事件。尤其在蒙特卡洛模擬法或耦合代理模型的蒙特卡洛模擬法中，過大的樣本數(shù)可能造成計算內(nèi)存占用過大甚至溢出，但子集模擬通過利用子集的概率乘積的方法可有效地降低計算矩陣維度，從而對小失效概率的邊坡可靠度進行有效分析。

（3）構(gòu)建克里金模型時，采用10倍隨機變量數(shù)的訓練樣本即可得到滿足計算精度需求的模型，而額外增加訓練樣本對計算結(jié)果影響較小。而回歸模型與相關模型對克里金模型計算精度有著較大的影響。對比不同情形的計算結(jié)果與文獻經(jīng)驗，本文建議采用二次型回歸模型與高斯型相關模型進行邊坡可靠度分析。