龔楊凱12，盧翠芳12，黃杰12，蘇國(guó)韶*12

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

0 引言

在地下工程穩(wěn)定性分析中，由于巖體介質(zhì)的高度復(fù)雜性和顯著的尺度效應(yīng)，室內(nèi)及現(xiàn)場(chǎng)的巖石力學(xué)試驗(yàn)往往不能夠合理地獲得巖體力學(xué)參數(shù)，如何合理地確定巖體力學(xué)參數(shù)一直是一個(gè)比較棘手的現(xiàn)實(shí)問(wèn)題[1]。利用巖體開(kāi)挖過(guò)程監(jiān)測(cè)到的位移或破壞區(qū)等實(shí)測(cè)信息進(jìn)行反分析，進(jìn)而推求定巖體參數(shù)的巖體參數(shù)分析方法是解決上述問(wèn)題的有效途徑。但是，對(duì)于復(fù)雜巖體工程，反分析的目標(biāo)優(yōu)化函數(shù)常具有表達(dá)式未知、高度非線性、多極值等特征[2]，傳統(tǒng)優(yōu)化方法難以獲取全局最優(yōu)解。近年來(lái)，學(xué)者們采用的遺傳算法(GA)、粒子群算法(PSO)、人工蜂群算法(ABC)等隨機(jī)全局優(yōu)化算法進(jìn)行反分析，取得了良好成效[3-5]。但對(duì)于洞室群等大型巖體工程的參數(shù)反分析，為保證數(shù)值計(jì)算的精度，計(jì)算單元致密且數(shù)量龐大，導(dǎo)致單次數(shù)值計(jì)算的耗時(shí)較大，若采用隨機(jī)全局優(yōu)化算法進(jìn)行反分析，常需要成千上萬(wàn)次地進(jìn)行數(shù)值計(jì)算，因計(jì)算耗時(shí)巨大導(dǎo)致所謂的高計(jì)算代價(jià)問(wèn)題。將機(jī)器學(xué)習(xí)模型與隨機(jī)全局優(yōu)化算法相結(jié)合是解決高計(jì)算代價(jià)問(wèn)題的有效途徑，利用機(jī)器學(xué)習(xí)模型替代數(shù)值計(jì)算模型，并建立巖體參數(shù)與數(shù)值計(jì)算結(jié)果的非線性映射關(guān)系，可顯著提高計(jì)算效率，其中，基于神經(jīng)網(wǎng)絡(luò)—遺傳算法(ANN-GA)以及支持向量機(jī)—遺傳算法(SVM-GA)的反分析方法應(yīng)用較為廣泛[1、6-11]，但這些方法尚存在著神經(jīng)網(wǎng)絡(luò)不適用于小樣本、合理的網(wǎng)絡(luò)結(jié)構(gòu)與超參數(shù)難以確定、易限于局部最優(yōu)解等局限性問(wèn)題。

高斯過(guò)程優(yōu)化(Gaussian Process Optimization, GPO)是近年出現(xiàn)的一種用于解決高計(jì)算代價(jià)問(wèn)題的優(yōu)化算法[12]。本文提出一種基于GPO的巖體參數(shù)反分析方法，為大型復(fù)雜地下工程巖體參數(shù)的合理確定提供一條新的途徑。

1 GPO算法的基本原理

Efficient Global Optimization(EGO)算法源于1998年，是GPO算法的前身[12、13]。GPO的原理源于廣泛應(yīng)用于求解函數(shù)表達(dá)式未知的優(yōu)化問(wèn)題的貝葉斯優(yōu)化方法[14-16]。GPO的基本思路是：利用少量訓(xùn)練樣本建立回歸模型，通過(guò)貝葉斯后驗(yàn)分布對(duì)待測(cè)區(qū)間進(jìn)行概率預(yù)測(cè)，依據(jù)某種準(zhǔn)則選取目標(biāo)函數(shù)極值的預(yù)測(cè)最優(yōu)解，將該解及其相應(yīng)真實(shí)函數(shù)值作為新樣本添入到訓(xùn)練樣本中，在利用新訓(xùn)練樣本集調(diào)整回歸模型，不斷重復(fù)此過(guò)程以獲得全局最優(yōu)解。它的基本實(shí)現(xiàn)步驟簡(jiǎn)述如下：

①根據(jù)少量訓(xùn)練樣本建立高斯回歸模型。

已知訓(xùn)練樣本集(x,y)，高斯過(guò)程回歸模型為[17-18]：

y=f(x)+ε=xTw+ε，

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

式中：K(X，X*)是測(cè)試X*點(diǎn)與訓(xùn)練樣本集合的所有輸入點(diǎn)X的n×1階協(xié)方差矩陣，k(X*，X*)是測(cè)試點(diǎn)X*本身的協(xié)方差矩陣，可分別簡(jiǎn)寫(xiě)為K(X*)、k(X*)。根據(jù)貝葉斯原理，給定新的輸入X*、訓(xùn)練集的輸入值X和觀察目標(biāo)值y的條件下，推斷出f*預(yù)測(cè)值的均值和方差為：

②基于LCB準(zhǔn)則獲取預(yù)測(cè)最優(yōu)點(diǎn)。

通過(guò)產(chǎn)生一個(gè)新的訓(xùn)練樣本(預(yù)測(cè)最優(yōu)解)提高回歸效果。如何選擇合理的新訓(xùn)練樣本是GPO算法的核心問(wèn)題。新樣本的選擇可通過(guò)對(duì)獲取函數(shù)(Acquisition Function)的極值來(lái)獲得，常用的獲取函數(shù)有PI(Probability of Improvement)、EI(Expected of Improvement)、LCB(Lower Confidence Bounds)等三種[11]，本文采用LCB準(zhǔn)則：

k=3-2.9×i/M，

其中：i為當(dāng)前迭代步iter，M為最大迭代步Maxiter。αLCB(x)可采用優(yōu)化算法求解，本文采用PSO算法。

③將預(yù)測(cè)最優(yōu)解代入真實(shí)函數(shù)或者數(shù)值計(jì)算模型，獲得一個(gè)新的訓(xùn)練樣本，將此樣本添入訓(xùn)練樣本集。

④重復(fù)步驟①～③，直至滿足收斂條件。

下面給出一個(gè)利用GPO算法求解一維多峰值函數(shù)最小值的示例，GPO采用SE協(xié)方差函數(shù)[10]。函數(shù)f(x)=xsin10πx+2-0.0516在[0, 0.9]區(qū)間的函數(shù)曲線見(jiàn)圖1(a)。假設(shè)有5個(gè)初始訓(xùn)練樣本(圖1(a))，采用GPO進(jìn)行尋優(yōu)，僅經(jīng)過(guò)6步迭代，就能迅速地獲得全局最優(yōu)解：x=0.7513,fmin=1.1977。

(a) 第0步

(b) 第2步

(c) 第4步

(d) 第6步

圖1 一維函數(shù)的GPO尋優(yōu)過(guò)程
Fig.1 An optimization processof one-dimensional function using GPO algorithm

2 GPO算法性能測(cè)試

為驗(yàn)證GPO算法的尋優(yōu)性能，將它分別與隨機(jī)全局優(yōu)化法PSO、機(jī)器學(xué)習(xí)與隨機(jī)全局優(yōu)化算法相結(jié)合的優(yōu)化算法GP-PSO[1]進(jìn)行對(duì)比。函數(shù)極小化采用 Sphere(單峰非連續(xù)函數(shù))、Schwefel(單峰連續(xù)函數(shù))和Griewank(多峰連續(xù)函數(shù))等三種常用的測(cè)試函數(shù)：

設(shè)優(yōu)化問(wèn)題維度n=5。為便于比較，PSO隨機(jī)產(chǎn)生的初始樣本作為GP-PSO、GPO的初始訓(xùn)練樣本。f1(x)、f3(x)的收斂準(zhǔn)則均為1×10-4，f2(x)的收斂準(zhǔn)則為1×10-3。

三種算法中的PSO算法參數(shù)設(shè)置：學(xué)習(xí)因子c1=2.1、c2=2.0，最大速度Vmax=[1, 1, 1, 1, 1]。GP采用SE協(xié)方差函數(shù)。計(jì)算結(jié)果見(jiàn)表1，表中三種算法的計(jì)算結(jié)果為10次計(jì)算的均值?？梢?jiàn)，與PSO相比，GP-PSO與GPO算法的函數(shù)調(diào)用次數(shù)降低了一個(gè)數(shù)量級(jí)；與GP-PSO相比，GPO的函數(shù)調(diào)用次數(shù)更少，例如，對(duì)于Schwefel函數(shù)，GPO通過(guò)22次實(shí)現(xiàn)收斂，比GP-PSO的函數(shù)調(diào)用次數(shù)少20 %。由此可見(jiàn)，GPO的優(yōu)化效率更佳。

表1 三種算法的函數(shù)調(diào)用次數(shù)對(duì)比Tab.1 Comparison of the number of function calls among three algorithms

3 巖體參數(shù)反分析的GPO方法

對(duì)GPO算法編制MATLAB程序，利用MATLAB命令調(diào)用并打開(kāi)FLAC3D數(shù)值計(jì)算軟件，通過(guò)FLAC3D的fish內(nèi)嵌語(yǔ)言編制數(shù)據(jù)接口文件讀取巖體參數(shù)變量并進(jìn)行數(shù)值計(jì)算獲取圍巖位移場(chǎng)，由此構(gòu)建適應(yīng)度函數(shù)表達(dá)式：

圖2 高斯過(guò)程優(yōu)化參數(shù)反分析方法的偽代碼圖
Fig.2 Pseudo code of parameter back analysis using GPO

4 工程算例

有一半徑為3 m的圓形斷面隧洞，水平和豎向初始地應(yīng)力分別為15 MPa、10 MPa。假定巖體的本構(gòu)模型為彈塑性模型，待反演的巖體力學(xué)參數(shù)分別為彈性模量E、粘聚力c和內(nèi)摩擦角φ。為驗(yàn)證本文方法的可行性，假設(shè)巖體的“真實(shí)力學(xué)參數(shù)”為：E=2.2 GPa、c=4.0 MPa、φ=4×10°。將此參數(shù)代入FlAC3D數(shù)值計(jì)算模型(見(jiàn)圖3，位移邊界采用法向約束)，計(jì)算其位移場(chǎng)，將開(kāi)挖面上的1、2、3、4、5測(cè)點(diǎn)的位移計(jì)算值視為“實(shí)測(cè)位移”。為驗(yàn)證本文方法的先進(jìn)性，分別采用PSO、GP-PSO與GPO三種方法進(jìn)行巖體力學(xué)參數(shù)反分析。

圖3 隧洞的FLAC3D數(shù)值計(jì)算模型Fig.3 FLAC3D numerical model of the tunnel

PSO初始參數(shù)設(shè)置：c1=2.1、c2=2.0，最大速度Vmax=[1, 1, 1]，尋優(yōu)區(qū)間為：E=[1.0～4.0]GPa、c=[3.0～5.0]MPa、φ= [3～5]×10°。GP-PSO與GPO均采用SE協(xié)方差函數(shù)。為便于比較，三種方法的初始訓(xùn)練樣本相同，見(jiàn)表2。

計(jì)算結(jié)果見(jiàn)表3。從中可知，當(dāng)收斂準(zhǔn)則設(shè)為f<0.3或數(shù)值計(jì)算最大次數(shù)300時(shí)，從數(shù)值計(jì)算次數(shù)來(lái)看，GPO分別為PSO、GP-PSO的14 %、60 %；從計(jì)算耗時(shí)來(lái)看，GPO分別為PSO、GP-PSO的16 %、58 %。當(dāng)收斂準(zhǔn)則為f<0.2時(shí)或數(shù)值計(jì)算最大次數(shù)200時(shí)，PSO算法無(wú)法收斂，GP-PSO的數(shù)值計(jì)算次數(shù)和計(jì)算耗時(shí)約為GPO的2倍。由此可見(jiàn)本文方法的先進(jìn)性。

表2 初始訓(xùn)練樣本Tab.2 Initial training samples

表3 三種方法的計(jì)算結(jié)果對(duì)比Tab.3 Comparison of the results of three methods

5 結(jié)論

本文提出一種基于GPO算法的巖體參數(shù)反分析方法。研究表明，本文方法是可行的，與基于PSO算法的反分析方法相比，本文方法的數(shù)值計(jì)算重分析次數(shù)顯著減小，同時(shí)克服了PSO的計(jì)算代價(jià)較高、尋優(yōu)結(jié)果波動(dòng)較大的不足；與基于GP-PSO算法的反分析方法相比，本文方法的數(shù)值計(jì)算重分析次數(shù)有一定程度降低，更不容易陷入局部最優(yōu)，且具有參數(shù)少、原理簡(jiǎn)單、實(shí)現(xiàn)容易等優(yōu)點(diǎn)。因此，本文方法在解決大規(guī)模、高計(jì)算代價(jià)的巖體參數(shù)反分析問(wèn)題上具有良好的發(fā)展?jié)摿蛻?yīng)用前景。需指出的是，初始訓(xùn)練樣本的建立對(duì)于尋優(yōu)結(jié)果具有重要影響，關(guān)于如何生成高質(zhì)量的初始訓(xùn)練樣本是下一步需開(kāi)展研究的課題。