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

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

0 引言

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

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

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

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

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

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

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

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

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

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

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

2 計(jì)算流程

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

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

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

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

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

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

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

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

3 算例

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

3.1 算例I：?jiǎn)螌诱承酝临|(zhì)邊坡

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

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

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

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

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

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

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

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

3.2 算例II：工程實(shí)例邊坡

算例II為一工程實(shí)例邊坡，該斜坡位于中國(guó)香港西貢區(qū)北部，由距坡頂9 m深的崩積層土體組成［27］。邊坡根據(jù)不同傾角分為5個(gè)區(qū)段，邊坡中部存有上覆物所致30 kN點(diǎn)荷載，地下水水位情況及相關(guān)工程地質(zhì)情況如圖4所示。土體重度γ=19 kN/m3，粘聚力c和摩擦角φ為服從對(duì)數(shù)正態(tài)分布的互相關(guān)隨機(jī)變量，參數(shù)統(tǒng)計(jì)特征如表5所示。采用SLOPE/W模塊中Morgenstern-Price法進(jìn)行極限平衡條分法邊坡穩(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)計(jì)特性Table 5 Statistics of the soil parameters of Example II

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

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

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

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

4 結(jié)論

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

（1）基于克里金模型與子集模擬的邊坡可靠度分析方法能夠有效地替代傳統(tǒng)蒙特卡洛模擬法，準(zhǔn)確地估計(jì)邊坡的失效概率?；诳死锝鸬淖蛹M法計(jì)算效率相較于傳統(tǒng)蒙特卡洛方法顯著提高。同時(shí)，該方法的計(jì)算精度可滿足可靠度分析的需求。

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

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