邊坡失效概率估計(jì)的高斯過(guò)程動(dòng)態(tài)響應(yīng)面法

蘇國(guó)韶，趙 偉，彭立鋒，燕柳斌

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

1 引 言

邊坡可靠性研究是當(dāng)前巖土工程界的研究熱點(diǎn)之一[1]。復(fù)雜邊坡工程的可靠性分析中，其功能函數(shù)一般是隱式的，這給常用的一次二階矩法（FORM）、二次二階矩法（SORM）、直接積分法等帶來(lái)了困難，直接蒙特卡羅模擬方法（MCS）適用于求解隱式功能函數(shù)的可靠性問(wèn)題，且計(jì)算精度較高，但為了保證計(jì)算精度，MCS 法要求的抽樣次數(shù)很大，尤其是對(duì)于功能函數(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)面方法仍存在著諸多問(wèn)題需要解決：①二次多項(xiàng)式難以逼近高度非線性的功能函數(shù)，導(dǎo)致計(jì)算誤差較大；ANN、RBF 等存在著小樣本問(wèn)題推廣能力欠佳和過(guò)度擬合等瓶頸問(wèn)題[7]；SVM 存在著合理核函數(shù)及其參數(shù)難以確定的公開(kāi)問(wèn)題[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 方法無(wú)需計(jì)算驗(yàn)算點(diǎn)，對(duì)求解功能函數(shù)呈高度非線性的可靠性問(wèn)題具有較強(qiáng)的適應(yīng)性，為此有必要探求新的基于MCS 法的響應(yīng)面方法。③現(xiàn)有的基于 MCS 的響應(yīng)面方法如ANN-MCS[3]、RBF-MCS[4]、SVM-MCS[10]等存在著過(guò)度依賴預(yù)設(shè)樣本的規(guī)模與分布的瓶頸問(wèn)題。為了降低計(jì)算代價(jià)，此類方法一般采用試驗(yàn)設(shè)計(jì)方法構(gòu)建少量的預(yù)設(shè)樣本，且僅利用一次構(gòu)建的預(yù)設(shè)樣本來(lái)構(gòu)造響應(yīng)面。一旦預(yù)設(shè)樣本過(guò)少或空間分布不合理或外推預(yù)測(cè)樣本遠(yuǎn)離預(yù)設(shè)樣本空間時(shí)，響應(yīng)面的重構(gòu)精度就難以保證，為此有必要在計(jì)算中引入響應(yīng)面動(dòng)態(tài)更新的機(jī)制，以有效降低可靠性分析方法對(duì)預(yù)設(shè)采樣點(diǎn)的依賴性。

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

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

2 高斯過(guò)程回歸

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

假設(shè)有n個(gè)觀察數(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)。

對(duì)于新的輸入X*，GPR 模型的任務(wù)是根據(jù)先驗(yàn)知識(shí)預(yù)測(cè)出與X*相對(duì)應(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階對(duì)稱正定的協(xié)方差矩陣，矩陣中的任一項(xiàng) kij度量了xi和xj的相關(guān)性。

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

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

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

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

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

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

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

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

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

3 高斯過(guò)程動(dòng)態(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個(gè)新的學(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個(gè)計(jì)算過(guò)程，直至收斂。值得指出的是，本文方法中，在MCS 隨機(jī)抽樣中，利用GPR 模型取代真實(shí)功能函數(shù)直接估計(jì)功能函數(shù)值，可以有效減小邊坡重分析的次數(shù)，顯著提高M(jìn)CS 法的計(jì)算效率，從而實(shí)現(xiàn)邊坡可靠性問(wèn)題的快速求解。

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

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

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

GPR 具有堅(jiān)實(shí)的統(tǒng)計(jì)學(xué)習(xí)理論基礎(chǔ)，能很好地處理高維數(shù)、小樣本和非線性的回歸問(wèn)題，具有實(shí)現(xiàn)容易、泛化能力強(qiáng)與參數(shù)自適應(yīng)獲取等優(yōu)點(diǎn)，對(duì)高度非線性函數(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 模型對(duì)加噪后生成的10個(gè)學(xué)習(xí)樣本（x,y）分別進(jìn)行回歸分析。比較各回歸模型在學(xué)習(xí)樣本被噪音腐蝕的情況下擬合真實(shí)函數(shù)的效果。從圖2 可看出，GPR 回歸效果最好，它克服了二次多項(xiàng)式不適用于高度非線性函數(shù)擬合、ANN 在小樣本問(wèn)題中的過(guò)度擬合等明顯問(wèn)題，而且模型超參數(shù)可自適應(yīng)獲取。因此，與多項(xiàng)式和ANN 模型相比較，GPR 模型在小樣本回歸性能上具有一定的優(yōu)越性。

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

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

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

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

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

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

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

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

研究發(fā)現(xiàn)，在學(xué)習(xí)樣本空間之內(nèi)GPR 模型的擬合精度相對(duì)較高，而在學(xué)習(xí)樣本空間之外GPR 模型的擬合精度相對(duì)較低，根據(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)面的擬合精度卻難以保證。為了克服這一問(wèn)題，需要在學(xué)習(xí)樣本空間之外選取最優(yōu)樣本點(diǎn)。

學(xué)習(xí)樣本空間之外的推斷具有較大的不確定性。為了盡量減小結(jié)構(gòu)重分析次數(shù)，本文采用小樣本的回歸策略，但是，以少量有限的信息去推斷完整的信息必將承受一定的風(fēng)險(xiǎn)，了解并控制推斷過(guò)程中的不確定性水平是保證所做推斷有意義的前提。因此，在小樣本推斷的條件下，模型輸出結(jié)果不確定性的合理描述是評(píng)價(jià)推斷合理性的一個(gè)有效手段。在統(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ù)測(cè)變量的均值和方差，方差可作為任意樣本功能函數(shù)預(yù)測(cè)值的不確定性度量指標(biāo)，即方差越大，該樣本的功能函數(shù)預(yù)測(cè)值的不確定性越大。

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

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

由GPR 的基本原理可知，GPR 預(yù)測(cè)方差較大的區(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í)的角度來(lái)看，式（11）等式右邊的分子部分控制了學(xué)習(xí)的精度，分母部分控制了學(xué)習(xí)的廣度。

以圖3為例，假設(shè)A 點(diǎn)與B 點(diǎn)的GPR 預(yù)測(cè)功能函數(shù)值幾乎相等且接近于0，為最優(yōu)樣本點(diǎn)的候選對(duì)象。與A 點(diǎn)相比較，B 點(diǎn)位于學(xué)習(xí)樣本空間之外，其預(yù)測(cè)方差較大，根據(jù)式（11），B 點(diǎn)的適應(yīng)度較小，故B 點(diǎn)為最優(yōu)樣本點(diǎn)。利用B 點(diǎn)的坐標(biāo)及其真實(shí)功能函數(shù)值，作為1個(gè)新的樣本添入原訓(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：通過(guò)邊坡分析程序，計(jì)算功能函數(shù)值g(x1,x2,…,xn)與g(x1,x2,…,xi±fσi,…,xn)，從而構(gòu)造2n+1個(gè)的學(xué)習(xí)樣本，其中x 取均值點(diǎn)，f 取2[2,6]。

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

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

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

步驟5：若前后迭代步之間的失效概率的相對(duì)差值小于1%，則停止迭代，否則，根據(jù)式（11），通過(guò)比對(duì)，從k個(gè)樣本點(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個(gè)高度非線性功能函數(shù)的數(shù)值算例，采用本文方法計(jì)算時(shí)可假設(shè)其功能函數(shù)是隱式的，用來(lái)考查本文方法的可行性。

4.1 數(shù)值算例1

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

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

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

從表1 可看出，預(yù)設(shè)學(xué)習(xí)樣本的GPR 預(yù)測(cè)值與真實(shí)功能函數(shù)值非常接近，說(shuō)明學(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ì)算值見(jiàn)表2。計(jì)算結(jié)果見(jiàn)表3，從中可見(jiàn)，以MCS 法模擬106次的結(jié)果為精確解時(shí)，與RSM、RBF-MCS、ANN-MCS 等響應(yīng)面方法相比較，本文方法的求解精度最高，函數(shù)調(diào)用次數(shù)明顯較少，其主要原因在于：①采用GPR 模型可以很好地?cái)M合高度非線性的真實(shí)功能函數(shù)；②采用GPR 響應(yīng)面動(dòng)態(tài)更新技術(shù)，實(shí)現(xiàn)了極少數(shù)量采樣點(diǎn)條件下GPR 響應(yīng)面擬合精度的迅速提升（參見(jiàn)圖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)面的動(dòng)態(tài)更新過(guò)程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個(gè)預(yù)設(shè)學(xué)習(xí)樣本，計(jì)算迭代4 步后收斂，計(jì)算結(jié)果見(jiàn)表4。從中可見(jiàn)，以MCS 法模擬107次的結(jié)果為精確解時(shí)，F(xiàn)ORM 法與RSM 法的相對(duì)誤差均高達(dá)313%，本文方法的相對(duì)誤差僅為4%，結(jié)果令人滿意。該算例從數(shù)值試驗(yàn)的角度驗(yàn)證了本文方法適用于求解高度非線性隱式功能函數(shù)可靠性問(wèn)題，能克服FORM 法、RSM 法等傳統(tǒng)方法求解此類問(wèn)題時(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ì)特征見(jiàn)表5。

計(jì)算結(jié)果見(jiàn)表6。從中可見(jiàn)，以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個(gè)算例，各算例均采用大型通用邊坡分析商業(yè)軟件SLIDE[16]提供的簡(jiǎn)化Bishop 法計(jì)算邊坡安全系數(shù)K。

5.1 邊坡算例1

有一均質(zhì)邊坡（見(jiàn)圖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個(gè)預(yù)設(shè)學(xué)習(xí)樣本見(jiàn)表8，計(jì)算迭代4 步收斂，各迭代步所獲得的最優(yōu)樣本見(jiàn)表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é)果見(jiàn)表10。從中可見(jiàn)，以MCS 法模擬106次的結(jié)果為相對(duì)精確解時(shí)，本文方法的相對(duì)誤差為2.4%，邊坡分析次數(shù)僅有11 次。若考慮隨機(jī)變量服從相互獨(dú)立的對(duì)數(shù)正態(tài)分布，對(duì)應(yīng)的計(jì)算結(jié)果見(jiàn)表10，從中可見(jiàn)，本文方法的相對(duì)誤差為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ì)邊坡（見(jiàn)圖6），土坡包含二土層，土的重度γ=19 kN/m3，設(shè)黏聚力c，內(nèi)摩擦角φ 是相互獨(dú)立的正態(tài)分布變量，統(tǒng)計(jì)特征如表11 所示，功能函數(shù)設(shè)為g=FS(c,φ)-1。構(gòu)造的7個(gè)預(yù)設(shè)學(xué)習(xí)樣本見(jiàn)表12。程序迭代4 步收斂，各迭代步所獲得的最優(yōu)樣本見(jiàn)表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é)果見(jiàn)表14，從中可知，以MCS 法模擬106次的結(jié)果為相對(duì)精確解時(shí)，本文方法的可靠指標(biāo)相對(duì)誤差為2.9%，邊坡分析次數(shù)僅有12 次，在計(jì)算精度與計(jì)算效率上明顯優(yōu)于RSM、ANN-FORM與二階 HDMR（high dimensional model representation）等不同響應(yīng)面方法。若考慮隨機(jī)變量服從相互獨(dú)立的對(duì)數(shù)正態(tài)分布，對(duì)應(yīng)的計(jì)算結(jié)果見(jiàn)表14，從中可見(jiàn)，與ANN-FORM 法[3]相比較，本文方法的計(jì)算精度較優(yōu)，邊坡分析次數(shù)明顯較少。算例結(jié)果驗(yàn)證本文方法對(duì)于求解邊坡可靠性問(wè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]，邊坡剖面圖及浸潤(rùn)線見(jiàn)圖7。邊坡各土層的力學(xué)參數(shù)見(jiàn)表15，其中黏土層1、黏土層2 與黏土層3 等土層的黏聚力c1、c2、c3為相互獨(dú)立的正態(tài)隨機(jī)變量?？紤]的邊坡失穩(wěn)模式為邊坡滑動(dòng)面與第3 層（黏土層3）的下邊界相切，對(duì)應(yīng)的功能函數(shù)設(shè)為：g=FS(c1,c2,c3)-1。構(gòu)造的7個(gè)預(yù)設(shè)樣本見(jiàn)表16。計(jì)算迭代2 步收斂，各迭代步的最優(yōu)樣本見(jiàn)表17。

計(jì)算結(jié)果見(jiàn)表18。以MCS 法抽樣106次[16]為相對(duì)精確解時(shí)，本文方法的計(jì)算精度明顯優(yōu)于二階HDMR 響應(yīng)面法[18]，邊坡分析次數(shù)約為二階HDMR響應(yīng)面法[18]的1/7。若考慮隨機(jī)變量服從相互獨(dú)立的對(duì)數(shù)正態(tài)分布，對(duì)應(yīng)的計(jì)算結(jié)果見(jiàn)表18，從中可見(jiàn)，本文方法邊坡分析次數(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ì)的高斯過(guò)程動(dòng)態(tài)響應(yīng)面方法，該方法既利用了MCS 法對(duì)于功能函數(shù)的形式、維數(shù)及其分布形式均無(wú)特殊要求，計(jì)算精度高的優(yōu)點(diǎn)，又充分利用了GPR 模型在處理小樣本、非線性回歸問(wèn)題上的優(yōu)勢(shì)，此外，該方法還利用了GPR 模型提供的可對(duì)預(yù)測(cè)結(jié)果不確定度進(jìn)行定量化評(píng)價(jià)的獨(dú)特功能，在此基礎(chǔ)上建立了響應(yīng)面最優(yōu)采樣點(diǎn)的合理配置機(jī)制，并結(jié)合蒙特卡羅隨機(jī)抽樣，實(shí)現(xiàn)了響應(yīng)面的高效率動(dòng)態(tài)更新，進(jìn)而有效地提升了響應(yīng)面法的計(jì)算精度，加快了迭代計(jì)算的收斂速度。

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

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

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

[1]劉曉,唐輝明,熊承仁.邊坡動(dòng)力可靠性分析方法的模式、問(wèn)題與發(fā)展趨勢(shì)[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]蔣水華,李典慶,方國(guó)光.結(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].中國(guó)管理科學(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.