蘇國韶，趙 偉，彭立鋒，燕柳斌

（廣西大學(xué) 土木建筑工程學(xué)院 工程防災(zāi)與結(jié)構(gòu)安全教育部重點(diǎn)實(shí)驗(yàn)室，南寧 530004）

1 引 言

邊坡可靠性研究是當(dāng)前巖土工程界的研究熱點(diǎn)之一[1]。復(fù)雜邊坡工程的可靠性分析中，其功能函數(shù)一般是隱式的，這給常用的一次二階矩法（FORM）、二次二階矩法（SORM）、直接積分法等帶來了困難，直接蒙特卡羅模擬方法（MCS）適用于求解隱式功能函數(shù)的可靠性問題，且計(jì)算精度較高，但為了保證計(jì)算精度，MCS 法要求的抽樣次數(shù)很大，尤其是對于功能函數(shù)值需要借助于計(jì)算耗時(shí)的有限元獲取時(shí)，大量的重分析導(dǎo)致計(jì)算量巨大，從而導(dǎo)致MCS 法在工程應(yīng)用中受到了極大限制。

利用少量的采樣點(diǎn)，采用二次多項(xiàng)式的經(jīng)典響應(yīng)面（RSM）[2]、人工神經(jīng)網(wǎng)絡(luò)（ANN）[3]、徑向基函數(shù)（RBF）[4]、克里金（Kriging）[5]以及支持向量機(jī)（SVM）[6]等回歸工具構(gòu)建隱式功能函數(shù)的響應(yīng)面，然后結(jié)合FORM、SORM 及MCS 等常規(guī)方法進(jìn)行可靠性分析，可有效降低結(jié)構(gòu)重分析次數(shù)，顯著提升計(jì)算效率，目前已成為復(fù)雜邊坡可靠性分析的重要途徑。但這類響應(yīng)面方法仍存在著諸多問題需要解決：①二次多項(xiàng)式難以逼近高度非線性的功能函數(shù)，導(dǎo)致計(jì)算誤差較大；ANN、RBF 等存在著小樣本問題推廣能力欠佳和過度擬合等瓶頸問題[7]；SVM 存在著合理核函數(shù)及其參數(shù)難以確定的公開問題[6]，為此有必要探求采用回歸性能更好的回歸工具，實(shí)現(xiàn)高度非線性隱式功能函數(shù)的高精度重構(gòu)。②當(dāng)功能函數(shù)非線性較高時(shí)，采用FORM、SORM 法往往難以收斂或獲得合理的驗(yàn)算點(diǎn)，這必將導(dǎo)致基于FORM 或SORM 的響應(yīng)面方法計(jì)算誤差較大[8－9]，MCS 方法無需計(jì)算驗(yàn)算點(diǎn)，對求解功能函數(shù)呈高度非線性的可靠性問題具有較強(qiáng)的適應(yīng)性，為此有必要探求新的基于MCS 法的響應(yīng)面方法。③現(xiàn)有的基于 MCS 的響應(yīng)面方法如ANN-MCS[3]、RBF-MCS[4]、SVM-MCS[10]等存在著過度依賴預(yù)設(shè)樣本的規(guī)模與分布的瓶頸問題。為了降低計(jì)算代價(jià)，此類方法一般采用試驗(yàn)設(shè)計(jì)方法構(gòu)建少量的預(yù)設(shè)樣本，且僅利用一次構(gòu)建的預(yù)設(shè)樣本來構(gòu)造響應(yīng)面。一旦預(yù)設(shè)樣本過少或空間分布不合理或外推預(yù)測樣本遠(yuǎn)離預(yù)設(shè)樣本空間時(shí)，響應(yīng)面的重構(gòu)精度就難以保證，為此有必要在計(jì)算中引入響應(yīng)面動態(tài)更新的機(jī)制，以有效降低可靠性分析方法對預(yù)設(shè)采樣點(diǎn)的依賴性。

高斯過程是新近發(fā)展起來的一種機(jī)器學(xué)習(xí)模型[11]，可分為高斯過程分類與高斯過程回歸（Gaussian process regression，GPR）兩種模型，它有嚴(yán)格的統(tǒng)計(jì)學(xué)習(xí)理論基礎(chǔ)，對于處理高維數(shù)、小樣本、非線性等復(fù)雜回歸問題具有良好的適應(yīng)性，在諸多領(lǐng)域中得到了廣泛應(yīng)用[12－14]。值得指出的是，與ANN、SVR 等常用的機(jī)器學(xué)習(xí)模型相比較，GPR 的優(yōu)勢主要體現(xiàn)在可給出具有概率意義的回歸結(jié)果以及模型最優(yōu)超參數(shù)可自適應(yīng)獲取。

本文提出一種邊坡失效概率快速估計(jì)的高斯過程動態(tài)響應(yīng)面法，該方法利用GPR 模型實(shí)現(xiàn)小樣本條件下非線性隱式功能函數(shù)的高精度逼近，通過迭代計(jì)算構(gòu)建一種GPR 響應(yīng)面擬合誤差的快速縮減機(jī)制，從而實(shí)現(xiàn)響應(yīng)面擬合精度的自適應(yīng)提高，進(jìn)而結(jié)合MCS 法計(jì)算邊坡的失效概率，為高效快速地求解邊坡的失效概率提供新的途徑。

2 高斯過程回歸

從統(tǒng)計(jì)學(xué)的角度來看，GPR是這樣的一個隨機(jī)過程：其任意有限變量集合的分布都是高斯分布，即對任意整數(shù)n≥1 及任意的一族隨機(jī)變量x，與其對應(yīng)的t時(shí)刻的過程狀態(tài)f(x)的聯(lián)合概率分布服從n維高斯分布。GPR 的全部統(tǒng)計(jì)特征完全由它的均值和協(xié)方差函數(shù)來確定。

假設(shè)有n個觀察數(shù)據(jù)的訓(xùn)練集D={(xi,yi)|i=1,…}n，xi為d 維輸入矢量，觀察目標(biāo)值yi∈R 。如果X 表示d×n 維輸入矩陣，y 表示輸出矢量，那么訓(xùn)練集 D=(X,y)。

對于新的輸入X*，GPR 模型的任務(wù)是根據(jù)先驗(yàn)知識預(yù)測出與X*相對應(yīng)的輸出值 y*。

帶高斯噪聲ε 的標(biāo)準(zhǔn)線性回歸模型為

式中：f(X)為回歸函數(shù)值；ε的各元素均為符合高斯分布的獨(dú)立隨機(jī)變量，均值為0，方差為。

觀察目標(biāo)值y 的先驗(yàn)分布為

式中：I為單位矩陣；K=K(X,X)為n×n階對稱正定的協(xié)方差矩陣，矩陣中的任一項(xiàng) kij度量了xi和xj的相關(guān)性。

n個訓(xùn)練樣本的觀察目標(biāo)值y 和m個測試樣本的回歸函數(shù)輸出f*所形成的聯(lián)合高斯先驗(yàn)分布為

式中：K(X,X*)為測試點(diǎn) X*與訓(xùn)練集的所有輸入點(diǎn)X 的n×m階協(xié)方差矩陣，可簡寫為K(X*)；k(X*,X*)為測試點(diǎn)X*自身的m×m協(xié)方差矩陣，可簡寫為 k(X*)。

GPR 模型常用的協(xié)方差函數(shù)為

式中的最優(yōu)超參數(shù)θ=(l,σf,σn)可通過極大似然法自適應(yīng)獲得，對數(shù)似然函數(shù)的形式為

對訓(xùn)練集進(jìn)行學(xué)習(xí)，根據(jù)貝葉斯原理，利用觀察到的真實(shí)數(shù)據(jù)不斷更新概率預(yù)測分布，即給定測試點(diǎn)X*、訓(xùn)練集的輸入值X 和觀察目標(biāo)值y 的條件下，推斷出 f*的最大可能的預(yù)測分布為

其中：預(yù)測分布的均值和方差為

GPR 模型是一個具有概率預(yù)測意義的模型，它能將認(rèn)知的先驗(yàn)信息和樣本信息有效結(jié)合并以概率分布的形式描述不確定性。不僅可以給出預(yù)測變量的均值，還能給出其方差和指定概率的置信區(qū)間，描述變量發(fā)生的不確定性程度。GPR 模型的最大特色之一在于以概率分布的形式描述預(yù)測結(jié)果，利用GRP 模型提供的預(yù)測方差，能十分方便地對回歸與外推預(yù)測結(jié)果進(jìn)行不確定性評價(jià)。

下面用一個函數(shù)回歸的例子來說明GPR 模型的不確定性評價(jià)功能。假設(shè)回歸的目標(biāo)函數(shù)為y=0.5sin(2x)，根據(jù)該式隨機(jī)構(gòu)造考慮噪音腐蝕條件下x∈[-2.5，2.5]的10個學(xué)習(xí)樣本，利用學(xué)習(xí)樣本建立GPR 模型，然后利用GPR 模型對x∈[-3，3]區(qū)間函數(shù)曲線進(jìn)行預(yù)測，預(yù)測效果見圖1，圖中的灰色條帶區(qū)域表示為置信水平為95%的置信區(qū)間（即2 倍標(biāo)準(zhǔn)差）。從圖1 可見，在x∈[-2.5，2.5]區(qū)間內(nèi)，灰色帶寬較小，表明了在訓(xùn)練樣本區(qū)間內(nèi)，預(yù)測結(jié)果的預(yù)測方差較小，即不確定性較?。辉趚∈[-2.5，2.5]區(qū)間外，距離訓(xùn)練樣本所在區(qū)間越遠(yuǎn)，灰色帶寬越大，說明距離訓(xùn)練樣本區(qū)間越遠(yuǎn)，預(yù)測結(jié)果的不確定性越大，進(jìn)而導(dǎo)致出現(xiàn)較大誤差的概率急劇增大（x=-3 或3時(shí)，曲線的偏差最大）。由此可見，GPR 模型的預(yù)測方差為預(yù)測結(jié)果的不確定性定理化評價(jià)提供了可能。

圖1 GPR 回歸預(yù)測示例Fig.1 Example for regression using GPR

3 高斯過程動態(tài)響應(yīng)面法

3.1 基本思路

本文方法的基本思路是：首先，利用邊坡分析程序生成少量學(xué)習(xí)樣本，用以訓(xùn)練GPR 模型，利用訓(xùn)練后的GPR 模型重構(gòu)功能函數(shù)，使隱式功能函數(shù)顯式化；然后，將GPR 模型作為真實(shí)功能函數(shù)的替代工具，結(jié)合MCS 法計(jì)算失效概率，在MCS 大量抽樣中選取失效域中的最接近于極限狀態(tài)方程且不確定性較大的點(diǎn)作為最優(yōu)樣本點(diǎn)，利用邊坡分析程序計(jì)算其真實(shí)功能函數(shù)值，由此獲得1個新的學(xué)習(xí)樣本，將該樣本添入原學(xué)習(xí)樣本集，GPR 模型重新學(xué)習(xí)，由此實(shí)現(xiàn)GPR 響應(yīng)面在極限狀態(tài)曲面附近擬合精度的自適應(yīng)提高，進(jìn)而再次結(jié)合MCS 法計(jì)算失效概率。不斷重復(fù)上述的GPR 模型學(xué)習(xí)、應(yīng)用MCS法計(jì)算失效概率以及搜索最優(yōu)樣本點(diǎn)等3個計(jì)算過程，直至收斂。值得指出的是，本文方法中，在MCS 隨機(jī)抽樣中，利用GPR 模型取代真實(shí)功能函數(shù)直接估計(jì)功能函數(shù)值，可以有效減小邊坡重分析的次數(shù)，顯著提高M(jìn)CS 法的計(jì)算效率，從而實(shí)現(xiàn)邊坡可靠性問題的快速求解。

3.2 關(guān)鍵技術(shù)

（1）構(gòu)建基于小樣本回歸策略的GPR 響應(yīng)面

在有限采樣點(diǎn)條件下，建立高性能的響應(yīng)面，實(shí)現(xiàn)真實(shí)隱式功能函數(shù)的高精度擬合與泛化，是降低可靠性分析方法對采樣點(diǎn)規(guī)模依賴的關(guān)鍵所在。

GPR 具有堅(jiān)實(shí)的統(tǒng)計(jì)學(xué)習(xí)理論基礎(chǔ)，能很好地處理高維數(shù)、小樣本和非線性的回歸問題，具有實(shí)現(xiàn)容易、泛化能力強(qiáng)與參數(shù)自適應(yīng)獲取等優(yōu)點(diǎn)，對高度非線性函數(shù)擬合具有良好的適用性。因此，僅以少量采樣點(diǎn)為代價(jià)，采用模型最優(yōu)超參數(shù)可自適應(yīng)獲取的GPR 模型重構(gòu)功能函數(shù)，可以方便快捷地建立滿足高效率精度要求的響應(yīng)面代理模型。

為了展示GPR 的回歸性能，下面以sinc 函數(shù)y=sin|x|/|x|（x=0時(shí)，y=1）為例，分別采用二次多項(xiàng)式、ANN 模型（BP 神經(jīng)網(wǎng)絡(luò)）和GPR 模型對加噪后生成的10個學(xué)習(xí)樣本（x,y）分別進(jìn)行回歸分析。比較各回歸模型在學(xué)習(xí)樣本被噪音腐蝕的情況下擬合真實(shí)函數(shù)的效果。從圖2 可看出，GPR 回歸效果最好，它克服了二次多項(xiàng)式不適用于高度非線性函數(shù)擬合、ANN 在小樣本問題中的過度擬合等明顯問題，而且模型超參數(shù)可自適應(yīng)獲取。因此，與多項(xiàng)式和ANN 模型相比較，GPR 模型在小樣本回歸性能上具有一定的優(yōu)越性。

（2）利用GPR 的動態(tài)學(xué)習(xí)更新響應(yīng)面

響應(yīng)面方法的計(jì)算精度主要取決于響應(yīng)面對極限狀態(tài)曲面的逼近程度。傳統(tǒng)的ANN-MCS[3]、RBF-MCS[4]、SVM-MCS[8]等響應(yīng)面法中，僅利用一次構(gòu)建的預(yù)設(shè)學(xué)習(xí)樣本來構(gòu)造響應(yīng)面。研究發(fā)現(xiàn)，當(dāng)預(yù)設(shè)樣本數(shù)過少或分布不合理時(shí)，對于高度非線性功能函數(shù)，高精度擬合是很難實(shí)現(xiàn)的。因此，為了提高計(jì)算效率，需要建立響應(yīng)面擬合誤差的高效修正機(jī)制。

圖2 各種回歸模型的擬合結(jié)果比較Fig.2 Regression result comparisons among different models

響應(yīng)面法的關(guān)鍵是響應(yīng)面對真實(shí)極限狀態(tài)曲面（g(x)=0 對應(yīng)的曲面）能很好地?cái)M合。筆者構(gòu)造一種迭代方式，以實(shí)現(xiàn)響應(yīng)面的動態(tài)更新：將GPR 模型作為真實(shí)功能函數(shù)的替代工具，結(jié)合MCS 法進(jìn)行大量隨機(jī)抽樣，然后從大量抽樣中選取對真實(shí)極限狀態(tài)曲面擬合效果影響最大點(diǎn)的作為最優(yōu)樣本點(diǎn)，計(jì)算該點(diǎn)的真實(shí)功能函數(shù)值，由此建立新的學(xué)習(xí)樣本，添入到原學(xué)習(xí)樣本集，應(yīng)用GPR 模型重新進(jìn)行學(xué)習(xí)，由此更新GPR 響應(yīng)面。重復(fù)上述過程，不斷利用上一迭代步的GPR 響應(yīng)面，結(jié)合MCS 大量抽樣，獲得新的最優(yōu)樣本點(diǎn)，通過GPR 模型學(xué)習(xí)不斷更新響應(yīng)面。

通過響應(yīng)面的動態(tài)更新，不僅能夠極大地降低響應(yīng)面對預(yù)設(shè)學(xué)習(xí)樣本的依賴性，而且能夠快速地改善響應(yīng)面的擬合效果，加快迭代計(jì)算的收斂速度。值得指出的是，在上述迭代中，每一迭代步僅需進(jìn)行1 次結(jié)構(gòu)分析，確保了極低的計(jì)算代價(jià)。因此，這個迭代方式可以實(shí)現(xiàn)響應(yīng)面對真實(shí)極限狀態(tài)曲面的擬合誤差的高效修正。

（3）利用GPR 不確定性評價(jià)功能獲取最優(yōu)樣本點(diǎn)

本文方法中，實(shí)現(xiàn)響應(yīng)面動態(tài)更新的關(guān)鍵在于對真實(shí)極限狀態(tài)曲面擬合效果影響最大點(diǎn)——最優(yōu)樣本點(diǎn)的確定。對于最優(yōu)樣本點(diǎn)的評判標(biāo)準(zhǔn)，一種直觀地想法是，鑒于極限狀態(tài)曲面附近的樣本點(diǎn)對響應(yīng)面的擬合精度影響最大，可將失效域內(nèi)最接近于真實(shí)極限狀態(tài)曲面的樣本點(diǎn)x*作為最優(yōu)樣本點(diǎn)，由此建立最優(yōu)樣本點(diǎn)的適應(yīng)度函數(shù)：

對于真實(shí)極限狀態(tài)函數(shù)，g(xi)=0，式（9）可改寫為

研究發(fā)現(xiàn)，在學(xué)習(xí)樣本空間之內(nèi)GPR 模型的擬合精度相對較高，而在學(xué)習(xí)樣本空間之外GPR 模型的擬合精度相對較低，根據(jù)式（10）估計(jì)的最優(yōu)樣本點(diǎn)往往落在學(xué)習(xí)樣本空間內(nèi)，GPR 重新學(xué)習(xí)后，學(xué)習(xí)樣本空間內(nèi)的GPR 響應(yīng)面將更逼近于真實(shí)極限狀態(tài)曲面，但在學(xué)習(xí)樣本空間之外，GPR 響應(yīng)面的擬合精度卻難以保證。為了克服這一問題，需要在學(xué)習(xí)樣本空間之外選取最優(yōu)樣本點(diǎn)。

學(xué)習(xí)樣本空間之外的推斷具有較大的不確定性。為了盡量減小結(jié)構(gòu)重分析次數(shù)，本文采用小樣本的回歸策略，但是，以少量有限的信息去推斷完整的信息必將承受一定的風(fēng)險(xiǎn)，了解并控制推斷過程中的不確定性水平是保證所做推斷有意義的前提。因此，在小樣本推斷的條件下，模型輸出結(jié)果不確定性的合理描述是評價(jià)推斷合理性的一個有效手段。在統(tǒng)計(jì)理論中，方差可以當(dāng)作不確定性的一種度量。例如在經(jīng)濟(jì)學(xué)中，人們常依據(jù)方差進(jìn)行投資項(xiàng)目的技術(shù)經(jīng)濟(jì)分析，在同等期望值條件下，方差值越小，則意味著投資項(xiàng)目的風(fēng)險(xiǎn)越小、回報(bào)較為穩(wěn)定[15－16]。GRP 模型的最大特色在于能給出預(yù)測變量的均值和方差，方差可作為任意樣本功能函數(shù)預(yù)測值的不確定性度量指標(biāo)，即方差越大，該樣本的功能函數(shù)預(yù)測值的不確定性越大。

為此，將預(yù)測方差 σ2引入式（10），可得到基于GPR 不確定性評價(jià)的最優(yōu)樣本點(diǎn)的適應(yīng)度函數(shù)：

上式的物理意義是：失效域內(nèi)，既接近于極限狀態(tài)曲面且預(yù)測不確定性較大的點(diǎn)為最優(yōu)樣本點(diǎn)。

由GPR 的基本原理可知，GPR 預(yù)測方差較大的區(qū)域，一般出現(xiàn)在學(xué)習(xí)樣本空間以外。根據(jù)式（11）進(jìn)行樣本點(diǎn)遴選，就能在學(xué)習(xí)樣本空間之外選取接近于極限狀態(tài)曲面且方差較大的最優(yōu)樣本點(diǎn)，由此有效地提升學(xué)習(xí)樣本空間以外的GPR 響應(yīng)面的擬合精度。從機(jī)器學(xué)習(xí)的角度來看，式（11）等式右邊的分子部分控制了學(xué)習(xí)的精度，分母部分控制了學(xué)習(xí)的廣度。

以圖3為例，假設(shè)A 點(diǎn)與B 點(diǎn)的GPR 預(yù)測功能函數(shù)值幾乎相等且接近于0，為最優(yōu)樣本點(diǎn)的候選對象。與A 點(diǎn)相比較，B 點(diǎn)位于學(xué)習(xí)樣本空間之外，其預(yù)測方差較大，根據(jù)式（11），B 點(diǎn)的適應(yīng)度較小，故B 點(diǎn)為最優(yōu)樣本點(diǎn)。利用B 點(diǎn)的坐標(biāo)及其真實(shí)功能函數(shù)值，作為1個新的樣本添入原訓(xùn)練樣本集，重新訓(xùn)練GRP 模型，就可以實(shí)現(xiàn)學(xué)習(xí)樣本空間之外的極限狀態(tài)響應(yīng)面的有效修正。

圖3 最優(yōu)樣本點(diǎn)選取示意圖Fig.3 Sketch of selecting the optimum sample points

3.3 實(shí)現(xiàn)步驟

假設(shè)結(jié)構(gòu)的功能函數(shù)為Z=g(x)，基本隨機(jī)變量x 的聯(lián)合概率密度函數(shù)為f(x)，則結(jié)構(gòu)失效概率為

式中：I[g(x)]為g(x)的示性函數(shù)，當(dāng)g(x)＜0時(shí)取1，否則取0。

采用MCS 法進(jìn)行大量隨機(jī)抽樣，設(shè)樣本容量為N，第k 次抽樣的樣本為xk，則pf的估計(jì)值為

根據(jù)式（7），利用GPR模型逼近真實(shí)功能函數(shù)：

將式（14）代入式（13），可得到基于GPR響應(yīng)面的失效概率估計(jì)式為

本文方法的實(shí)現(xiàn)步驟如下：

步驟1：通過邊坡分析程序，計(jì)算功能函數(shù)值g(x1,x2,…,xn)與g(x1,x2,…,xi±fσi,…,xn)，從而構(gòu)造2n+1個的學(xué)習(xí)樣本，其中x 取均值點(diǎn)，f 取2[2,6]。

步驟2：采用GPR 模型對學(xué)習(xí)樣本進(jìn)行學(xué)習(xí)（自回歸），建立真實(shí)功能函數(shù)的GPR 響應(yīng)面，由此建立隨機(jī)變量與功能函數(shù)值之間的非線性映射關(guān)系。

步驟3：根據(jù)各隨機(jī)變量的統(tǒng)計(jì)特征，隨機(jī)抽取N個（本文N=106）樣本點(diǎn)，根據(jù)式（14），利用學(xué)習(xí)后的GPR 響應(yīng)面推求相應(yīng)地近似功能函數(shù)值。

步驟4：根據(jù)式（15）估計(jì)失效概率，若失效概率為0，返回步驟1，f 取上一次取值的1/2，增加2n個學(xué)習(xí)樣本，重復(fù)步驟1～4，直至失效概率不為0。

步驟5：若前后迭代步之間的失效概率的相對差值小于1%，則停止迭代，否則，根據(jù)式（11），通過比對，從k個樣本點(diǎn)中確定最優(yōu)樣本點(diǎn)。

步驟6：將最優(yōu)樣本點(diǎn)及其真實(shí)功能函數(shù)值作為新的最優(yōu)樣本，添入原學(xué)習(xí)樣本集，返回步驟（2）。

根據(jù)上述步驟，編制了MATLAB 計(jì)算程序。

4 方法驗(yàn)證

下面給出了2個高度非線性功能函數(shù)的數(shù)值算例，采用本文方法計(jì)算時(shí)可假設(shè)其功能函數(shù)是隱式的，用來考查本文方法的可行性。

4.1 數(shù)值算例1

有一非線性功能函數(shù)為

隨機(jī)變量x1和x2均服從標(biāo)準(zhǔn)正態(tài)分布。

表1 列出了5個預(yù)設(shè)學(xué)習(xí)樣本，采用GPR 模型進(jìn)行學(xué)習(xí)，自適應(yīng)獲得的最優(yōu)超參數(shù)對數(shù)值為lnθ=(1.50,2.91,7.86,-6.75)，預(yù)設(shè)學(xué)習(xí)樣本的GPR 預(yù)測值可按下式計(jì)算：

從表1 可看出，預(yù)設(shè)學(xué)習(xí)樣本的GPR 預(yù)測值與真實(shí)功能函數(shù)值非常接近，說明學(xué)習(xí)后的GPR 模型可以在學(xué)習(xí)樣本空間內(nèi)很好地逼近真實(shí)功能函數(shù)，從而可以替代真實(shí)功能函數(shù)進(jìn)行蒙特卡羅抽樣。

表1 預(yù)設(shè)學(xué)習(xí)樣本及其功能函數(shù)值Table 1 The learning samples and function values

計(jì)算迭代7 步后收斂，各迭代步的最優(yōu)樣本點(diǎn)及其失效概率計(jì)算值見表2。計(jì)算結(jié)果見表3，從中可見，以MCS 法模擬106次的結(jié)果為精確解時(shí)，與RSM、RBF-MCS、ANN-MCS 等響應(yīng)面方法相比較，本文方法的求解精度最高，函數(shù)調(diào)用次數(shù)明顯較少，其主要原因在于：①采用GPR 模型可以很好地?cái)M合高度非線性的真實(shí)功能函數(shù)；②采用GPR 響應(yīng)面動態(tài)更新技術(shù)，實(shí)現(xiàn)了極少數(shù)量采樣點(diǎn)條件下GPR 響應(yīng)面擬合精度的迅速提升（參見圖4）。

表2 各迭代步的計(jì)算結(jié)果Table 2 Results at different iterative steps

表3 數(shù)值算例1 的計(jì)算結(jié)果Table 3 Results of numerical example 1

圖4 GPR 響應(yīng)面的動態(tài)更新過程Fig.4 Dynamic updating process of GPR response surface

4.2 數(shù)值算例2

有一高次非線性功能函數(shù)為

隨機(jī)變量x1、x2和x3均服從標(biāo)準(zhǔn)正態(tài)分布。

生成7個預(yù)設(shè)學(xué)習(xí)樣本，計(jì)算迭代4 步后收斂，計(jì)算結(jié)果見表4。從中可見，以MCS 法模擬107次的結(jié)果為精確解時(shí)，F(xiàn)ORM 法與RSM 法的相對誤差均高達(dá)313%，本文方法的相對誤差僅為4%，結(jié)果令人滿意。該算例從數(shù)值試驗(yàn)的角度驗(yàn)證了本文方法適用于求解高度非線性隱式功能函數(shù)可靠性問題，能克服FORM 法、RSM 法等傳統(tǒng)方法求解此類問題時(shí)計(jì)算誤差偏大或計(jì)算不收斂的局限性。

表4 數(shù)值算例2 的計(jì)算結(jié)果Table 4 Results of numerical example 2

4.3 數(shù)值算例3

有一非線性功能函數(shù)為

隨機(jī)變量x1、x2和x3的統(tǒng)計(jì)特征見表5。

計(jì)算結(jié)果見表6。從中可見，以MCS 法模擬107次的結(jié)果為精確解時(shí)，與RSM 法相比較，本文方法計(jì)算精度較高，且函數(shù)調(diào)用次數(shù)明顯較少。

表5 隨機(jī)變量統(tǒng)計(jì)特征Table 5 Properties of stochastic variables

表6 數(shù)值算例3 的計(jì)算結(jié)果Table 6 Results of numerical example 3

5 邊坡算例

下面給出邊坡可靠性分析的3個算例，各算例均采用大型通用邊坡分析商業(yè)軟件SLIDE[16]提供的簡化Bishop 法計(jì)算邊坡安全系數(shù)K。

5.1 邊坡算例1

有一均質(zhì)邊坡（見圖5），設(shè)土的重度γ、黏聚力c 和內(nèi)摩擦角φ 是相互獨(dú)立的正態(tài)變量，統(tǒng)計(jì)特征如表7 所示。功能函數(shù)設(shè)為g=FS(c，φ，γ)-1，F(xiàn)S(·)為隱式函數(shù)，其函數(shù)值等于邊坡安全系數(shù)K。構(gòu)造的7個預(yù)設(shè)學(xué)習(xí)樣本見表8，計(jì)算迭代4 步收斂，各迭代步所獲得的最優(yōu)樣本見表9。

圖5 單層邊坡剖面圖（單位：m）Fig.5 Geometry of one-layer slope(unit:m)

表7 單層邊坡隨機(jī)變量統(tǒng)計(jì)特征Table 7 Properties of stochastic variables of one-layer slope

表8 單層邊坡預(yù)設(shè)學(xué)習(xí)樣本Table 8 Initial learning samples of one-layer slope

表9 單層邊坡各迭代步的最優(yōu)樣本Table 9 Optimum sample at different iterative steps of one-layer slope

計(jì)算結(jié)果見表10。從中可見，以MCS 法模擬106次的結(jié)果為相對精確解時(shí)，本文方法的相對誤差為2.4%，邊坡分析次數(shù)僅有11 次。若考慮隨機(jī)變量服從相互獨(dú)立的對數(shù)正態(tài)分布，對應(yīng)的計(jì)算結(jié)果見表10，從中可見，本文方法的相對誤差為1.1%，邊坡分析次數(shù)僅有14 次。

表10 邊坡算例1 的計(jì)算結(jié)果Table 10 Results of slope example 1 of one-layer slope

上述分析表明，本文方法計(jì)算精度令人滿意，與MCS 法相比較，獲取較高計(jì)算精度條件下所需要的邊坡分析次數(shù)非常少。

5.2 二層邊坡

有一非均質(zhì)邊坡（見圖6），土坡包含二土層，土的重度γ=19 kN/m3，設(shè)黏聚力c，內(nèi)摩擦角φ 是相互獨(dú)立的正態(tài)分布變量，統(tǒng)計(jì)特征如表11 所示，功能函數(shù)設(shè)為g=FS(c,φ)-1。構(gòu)造的7個預(yù)設(shè)學(xué)習(xí)樣本見表12。程序迭代4 步收斂，各迭代步所獲得的最優(yōu)樣本見表13。

圖6 二層邊坡剖面圖（單位：m）Fig.6 Geometry of two-layer slope(unit:m)

表11 二層邊坡隨機(jī)變量統(tǒng)計(jì)特征Table 11 Properties of stochastic variables of two-layer slope

表12 二層邊坡預(yù)設(shè)學(xué)習(xí)樣本Table 12 Initial learning samples of two-layer slope

計(jì)算結(jié)果見表14，從中可知，以MCS 法模擬106次的結(jié)果為相對精確解時(shí)，本文方法的可靠指標(biāo)相對誤差為2.9%，邊坡分析次數(shù)僅有12 次，在計(jì)算精度與計(jì)算效率上明顯優(yōu)于RSM、ANN-FORM與二階 HDMR（high dimensional model representation）等不同響應(yīng)面方法。若考慮隨機(jī)變量服從相互獨(dú)立的對數(shù)正態(tài)分布，對應(yīng)的計(jì)算結(jié)果見表14，從中可見，與ANN-FORM 法[3]相比較，本文方法的計(jì)算精度較優(yōu)，邊坡分析次數(shù)明顯較少。算例結(jié)果驗(yàn)證本文方法對于求解邊坡可靠性問題的可行性與優(yōu)越性。

表13 二層邊坡各迭代步的最優(yōu)樣本Table 13 Optimum sample at different iterative steps of two layer slope

表14 邊坡算例2 的計(jì)算結(jié)果Table 14 Results of slope example 2

5.3 多層邊坡

Congress St.Cut model 是驗(yàn)證邊坡可靠性分析方法的經(jīng)典算例之一[18,20]，邊坡剖面圖及浸潤線見圖7。邊坡各土層的力學(xué)參數(shù)見表15，其中黏土層1、黏土層2 與黏土層3 等土層的黏聚力c1、c2、c3為相互獨(dú)立的正態(tài)隨機(jī)變量?？紤]的邊坡失穩(wěn)模式為邊坡滑動面與第3 層（黏土層3）的下邊界相切，對應(yīng)的功能函數(shù)設(shè)為：g=FS(c1,c2,c3)-1。構(gòu)造的7個預(yù)設(shè)樣本見表16。計(jì)算迭代2 步收斂，各迭代步的最優(yōu)樣本見表17。

計(jì)算結(jié)果見表18。以MCS 法抽樣106次[16]為相對精確解時(shí)，本文方法的計(jì)算精度明顯優(yōu)于二階HDMR 響應(yīng)面法[18]，邊坡分析次數(shù)約為二階HDMR響應(yīng)面法[18]的1/7。若考慮隨機(jī)變量服從相互獨(dú)立的對數(shù)正態(tài)分布，對應(yīng)的計(jì)算結(jié)果見表18，從中可見，本文方法邊坡分析次數(shù)僅有10 次，但計(jì)算精度令人滿意。

圖7 多層邊坡剖面圖（單位：m）Fig.7 Geometry of multilayer slope(unit:m)

表15 多層邊坡土體計(jì)算參數(shù)Table 15 Material properties of soil of multilayer slope

表16 多層邊坡預(yù)設(shè)學(xué)習(xí)樣本Table 16 The initial learning samples of multilayer slope

表17 多層邊坡各迭代步的最優(yōu)樣本Table 17 The optimum sample at different iterative steps of multilayer slope

表18 邊坡算例3 的計(jì)算結(jié)果Table 18 Results of slope example 3

6 結(jié) 論

（1）本文將GPR 響應(yīng)面法與MCS 法相結(jié)合，提出了一種邊坡失效概率估計(jì)的高斯過程動態(tài)響應(yīng)面方法，該方法既利用了MCS 法對于功能函數(shù)的形式、維數(shù)及其分布形式均無特殊要求，計(jì)算精度高的優(yōu)點(diǎn)，又充分利用了GPR 模型在處理小樣本、非線性回歸問題上的優(yōu)勢，此外，該方法還利用了GPR 模型提供的可對預(yù)測結(jié)果不確定度進(jìn)行定量化評價(jià)的獨(dú)特功能，在此基礎(chǔ)上建立了響應(yīng)面最優(yōu)采樣點(diǎn)的合理配置機(jī)制，并結(jié)合蒙特卡羅隨機(jī)抽樣，實(shí)現(xiàn)了響應(yīng)面的高效率動態(tài)更新，進(jìn)而有效地提升了響應(yīng)面法的計(jì)算精度，加快了迭代計(jì)算的收斂速度。

（2）數(shù)值算例與邊坡算例的研究結(jié)果表明，本文方法是可行的，它克服了經(jīng)典響應(yīng)面法求解高度非線性功能函數(shù)可靠性問題時(shí)計(jì)算誤差大的不足，解決了MCS 法計(jì)算效率低、現(xiàn)有的“響應(yīng)面+MCS”混合方法的計(jì)算精度過度依賴預(yù)設(shè)樣本的規(guī)模與分布等瓶頸問題。

（3）與傳統(tǒng)邊坡可靠性分析方法相比較，本文方法具有計(jì)算精度高、估計(jì)速度快、易于與已有的邊坡分析軟件相結(jié)合的優(yōu)點(diǎn)，計(jì)算過程中無需對功能函數(shù)求導(dǎo)數(shù)，故實(shí)現(xiàn)較為容易，便于工程應(yīng)用，尤其適用于結(jié)構(gòu)分析代價(jià)較高且具有高度非線性隱式功能函數(shù)的可靠性問題。

（4）受篇幅限制，本文算例僅考慮了隨機(jī)變量相互獨(dú)立、邊坡滑面假設(shè)為圓弧的情況，對于隨機(jī)變量相關(guān)性與非規(guī)則邊坡滑面的情況，將另文闡述。

[1]劉曉,唐輝明,熊承仁.邊坡動力可靠性分析方法的模式、問題與發(fā)展趨勢[J].巖土力學(xué),2013,34(5):1217－1234.LIU Xiao,TANG Hui-ming,XIONG Chen-ren.Patterns,problems,and development trends of analysis methods for slope dynamic reliability[J].Rock and Soil Mechanics,2013,34(5):1217－1234.

[2]KIM S H,NA S W.Response surface method using vector projected sampling points[J].Structural Safety,1997,19:3－19.

[3]CHO S E.Probabilistic stability analyses of slopes using the ANN-based response surface[J].Computers and Geotechnics,2009,36(5):787－797.

[4]DENG J.Structural reliability analysis for implicit performance function using radial basis function network[J].International Journal of Solids and Structures,2006,43(11－12):3255－3291.

[5]LUO X F,LI X,ZHOU J,et al.A Kriging-based hybrid optimization algorithm for slope reliability analysis[J].Structure Safety,2012,34(1):401－406.

[6]ZHAO H B.Slope reliability analysis using a support vector machine[J].Computers and Geotechnics,2008,35(3):459－467.

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

[8]LI D Q,CHEN Y F,LU W B,et al.Stochastic response surface method for reliability analysis of rock slopes involving correlated non-normal variables[J].Computers and Geotechnics,2011,38(1):58－68.

[9]蔣水華,李典慶,方國光.結(jié)構(gòu)可靠度分析的響應(yīng)面法和隨機(jī)響應(yīng)面法的比較[J].武漢大學(xué)學(xué)報(bào)(工學(xué)版),2012,45(1):46－53.JING Shui-hua,LI Dian-qing,PHOON Kokkwang.A comparative study of response surface method and stochastic response surface method for structure reliability analysis[J].Engineering Journal of Wuhan University,2012,45(1):46－53.

[10]HURTADO J E.An examination of methods for approximating implicit limit state functions from the viewpoint of statistical learning theory[J].Structural Safety,2004,26(3):271－293.

[11]RASMUSSEN C E，WILLIAMS C K I.GaussianProcesses for Machine Learning[M].Massachusetts:MIT Press,2006.

[12]SEEGER M.Gaussian processes for machine learning[J].International Journal of Neural System,2004,14(2):69－106.

[13]GIROLAMI M,ROGERS S.Variational Bayesian multinomial probit regression with Gaussian process priors[J].Neural Computation,2006,18(8):1790－1817.

[14]SU G S，YAN L B，SONG Y C.Gaussian process for non-linear displacement time series prediction of landslide[J].Journal of China University of Geosciences,2007,18：219－212.

[15]MARKOWITZ H.Portfolio selection[J].Journal of Finance,1952,7:77－91.

[16]朱奉云,邱菀華,劉善存.投資組合均值-方差模型和極小極大模型的實(shí)證比較[J].中國管理科學(xué),2002,10(6):13－17.ZHU Feng-yun,QIU Wan-hua,LIU Shan-cun.An empirical comparison between mean-variance model and minimax model in Chinese stock market[J].Chinese Journal of Management Science,2002,10(6):13－17

[17]呂震宙,楊子政.基于神經(jīng)網(wǎng)絡(luò)的可靠性分析新方法[J].機(jī)械強(qiáng)度,2006,28(5):699－702.LU˙˙Zhen-zhou,YANG Zi-zheng.New reliability analysis method based on artificial neural network[J].Journal of Mechanical Strength,2006,28(5):699－702.

[18]ROCSCIENCE Inc.Slide 2D limit equilibrium slope stability for soil and rock slopes[R].[S.l.]:[s.n.],2006.

[19]XU B,LOW B K.Probabilistic stability analyses of embankments based on finite element method[J].Journal of Geotechnical and Geoenvironmental Engineering,2006,132(11):1444－1454.

[20]CHOWDHURY R N,RAO B N.Probabilistic stability assessment of slopes using high dimensional model representation[J].Computers and Geotechnics,2010,37(7－8):876－884.