李建合， 孫偉哲， 蘇國韶,3*

(1.廣西新發(fā)展交通集團(tuán)有限公司， 南寧 530007； 2.廣西大學(xué)土木建筑工程學(xué)院， 南寧 530004； 3.廣西巖溶區(qū)水安全與智慧調(diào)控工程研究中心， 南寧 530004)

交通隧道、引水隧洞或地下發(fā)電廠房等地下巖體工程設(shè)計(jì)與施工中，確定巖體力學(xué)參數(shù)是一項(xiàng)重要的基礎(chǔ)工作。雖然一般巖體力學(xué)參數(shù)可通過室內(nèi)試驗(yàn)或原位測(cè)試得到，但受尺寸效應(yīng)和空間分布不均等方面的影響，試驗(yàn)或測(cè)試結(jié)果與實(shí)際情況往往存在較大出入。如何準(zhǔn)確合理地確定巖體力學(xué)參數(shù)一直以來都是地下工程設(shè)計(jì)與施工中經(jīng)常遇到的問題，其中，利用反演方法推演巖體力學(xué)參數(shù)可以較好地解決這一問題[1]。

通過設(shè)定目標(biāo)函數(shù)、優(yōu)化變量與約束條件，將巖體參數(shù)確定問題轉(zhuǎn)化為優(yōu)化問題，借助優(yōu)化法搜索較優(yōu)的參數(shù)，這類方法稱為巖體力學(xué)參數(shù)優(yōu)化反演法[2]。由于巖體位移監(jiān)測(cè)較為方便，使得基于位移構(gòu)建目標(biāo)函數(shù)與約束條件的優(yōu)化反演法得到了廣泛應(yīng)用[3]。此外，針對(duì)單純形法[4]、共軛梯度法[5]、復(fù)合形法[6]等傳統(tǒng)優(yōu)化算法往往難以獲得全局最優(yōu)解的問題，自20世紀(jì)90年代以來，遺傳算法 (GA)[7-8]、粒子群優(yōu)化算法 (PSO)[9-10]、差分進(jìn)化 (DE)[11]、魚群算法 (IAF)[12]等仿生優(yōu)化算法被應(yīng)用于優(yōu)化反演問題，可以克服傳統(tǒng)優(yōu)化算法的不足。但是，仿生優(yōu)化算法是一種隨機(jī)優(yōu)化算法，對(duì)于非凸全局優(yōu)化問題，往往需要成千上萬次地調(diào)用目標(biāo)函數(shù)以評(píng)價(jià)個(gè)體的優(yōu)劣。復(fù)雜地下工程的巖體力學(xué)參數(shù)優(yōu)化反演中，基于位移變量的目標(biāo)函數(shù)值需要耗時(shí)的數(shù)值計(jì)算才能獲得，成千上萬次地調(diào)用目標(biāo)函數(shù)意味著大量地重復(fù)數(shù)值計(jì)算，導(dǎo)致耗時(shí)激增，影響了仿生優(yōu)化算法的良好應(yīng)用。近年來，有學(xué)者提出采用代理模型替代部分?jǐn)?shù)值計(jì)算以降低數(shù)值計(jì)算調(diào)用次數(shù)，周正軍等[13]提出了土石壩位移反演的響應(yīng)面遺傳算法，顯著提高了計(jì)算效率；關(guān)永平等[14]利用遺傳算法優(yōu)化的BP神經(jīng)網(wǎng)絡(luò)(ANN)反演圍巖力學(xué)參數(shù)，取得了良好成效；徐國文等[15]分別采用粒子群優(yōu)化算法、遺傳算法、交叉驗(yàn)證算法與支持向量機(jī)代理模型結(jié)合進(jìn)行比較研究，證實(shí)了仿生優(yōu)化算法與支持向量機(jī)(SVM)相結(jié)合的反演方法的優(yōu)越性；倪沙沙等[16]在反演土石壩滲流場(chǎng)參數(shù)時(shí)，利用SVM建立滲流參數(shù)與水頭之間的非線性關(guān)系，并通過PSO算法識(shí)別滲流參數(shù)，提高了計(jì)算效率。但是，上述基于代理模型的優(yōu)化反演法尚存在著一些局限性，如，ANN模型易出現(xiàn)過度擬合[17-18]、SVM模型的最優(yōu)超參數(shù)難以合理確定[19-20]、代理模型的預(yù)測(cè)性能過度依賴訓(xùn)練樣本設(shè)置質(zhì)量，訓(xùn)練樣本的數(shù)量與分布設(shè)置不佳易導(dǎo)致優(yōu)化反演陷入局部最優(yōu)等問題。

針對(duì)上述問題，現(xiàn)將高斯過程回歸(Gaussian process regression，GPR)機(jī)器學(xué)習(xí)模型作為數(shù)值計(jì)算的代理模型，將蜜獾仿生優(yōu)化算法(honey badger algorithm，HBA)作為優(yōu)化工具，結(jié)合三維快速拉格朗日數(shù)值計(jì)算(FLAC3D)方法，提出巖體力學(xué)參數(shù)反演的一種代理蜜獾優(yōu)化方法(HBA-GPR-FLAC3D)，以期為復(fù)雜地下工程巖體力學(xué)參數(shù)快速確定提供新的有效手段。

1 HBA和GPR的基本原理

1.1 HBA原理簡(jiǎn)介

HBA是2021年出現(xiàn)的一種新型仿生優(yōu)化算法[21]。與已有仿生優(yōu)化方法相比，HBA具有收斂快、適應(yīng)度函數(shù)調(diào)用次數(shù)的特點(diǎn)。算法思想源于蜜獾覓食行為特點(diǎn)：為了尋找食物來源，蜜獾要么通過自身嗅覺，要么通過跟隨向蜜鳥尋找蜂巢。蜜獾會(huì)在獵物周圍移動(dòng)，選擇合適的區(qū)域搜索獵物的行為稱為“挖掘”； 蜜獾借助向蜜鳥的指引直接定位蜂巢的行為稱為“尋蜜”，如圖1所示。HBA的尋優(yōu)分為挖掘階段和尋蜜階段，在挖掘階段，通常將算子擴(kuò)展到較遠(yuǎn)的區(qū)域來保證搜索空間中包含足夠的潛在目標(biāo)，在尋蜜階段，算子逐步向潛在的目標(biāo)區(qū)域靠攏，進(jìn)行局部搜索。算法原理詳述如下。

圖1 蜜獾的覓食行為特點(diǎn)Fig.1 Characteristics of foraging behavior of honey badger

蜜獾算子的位置表達(dá)式為

xi=lbi+r1(ubi-lbi)

(1)

在蜜獾追蹤獵物時(shí)，蜜獾的追蹤強(qiáng)度與獵物的集中強(qiáng)度以及獵物與第i個(gè)蜜獾之間的距離有關(guān)。HBA中，蜜獾的追蹤強(qiáng)度可由反平方定律表示[22]，當(dāng)目標(biāo)算子集中，或當(dāng)蜜獾與目標(biāo)算子之間的距離縮小，蜜獾的追蹤強(qiáng)度便會(huì)加大，表達(dá)式為

(2)

式(2)中：r2為0～1的隨機(jī)數(shù)；S為目標(biāo)算子的集中強(qiáng)度；di為第i只蜜獾與目標(biāo)間的距離；xprey為獵物的位置。

密度因子α通過控制隨時(shí)間變化的隨機(jī)化，從而確保從全局搜索到局部搜索的平穩(wěn)過渡。使用式(3)更新隨迭代而減小的遞減因子α，以減小隨時(shí)間變化的隨機(jī)化，表達(dá)式為

(3)

式(3)中：tmax為最大迭代次數(shù)；C為大于等于1的常數(shù)(默認(rèn)為2)。

在挖掘階段，蜜獾主要依賴于獵物的氣味強(qiáng)度I、蜜獾與獵物之間的距離di和時(shí)變搜索影響因子α，此外，控制值F也是使蜜獾能夠更好地找到獵物位置的重要因素。蜜獾的動(dòng)作為如圖2所示的心形，其運(yùn)動(dòng)軌跡可由式(4)模擬獲得，即

xnew=xprey+Fβxprey+Fr3αdi×

|cos(2πr4)[1-cos(2πr5)]|

(4)

式(4)中:xprey為當(dāng)前目標(biāo)算子的最佳位置，即全局最佳位置；β一般取6，表征為蜜獾獲得食物的能力；di為獵物與第i只蜜獾的距離；r3、r4和r5為0～1的隨機(jī)數(shù)；F為改變搜索方向的控制值，表達(dá)式為

(5)

在尋蜜階段，蜜獾跟隨向蜜鳥到達(dá)蜂巢的動(dòng)作可通過式(6)模擬得到，即

xnew=xprey+Fr7αdi

(6)

式(6)中：xnew為蜜獾的新位置；F和α分別由式(5)和式(6)確定；r7為0～1的隨機(jī)數(shù)。

藍(lán)色輪廓表示氣味強(qiáng)度；黑色圓形線表示獵物位置圖2 挖掘階段Fig.2 Digging phase

1.2 GPR原理簡(jiǎn)介

GPR是一種基于貝葉斯理論的機(jī)器學(xué)習(xí)方法，具有非參數(shù)和概率性的優(yōu)點(diǎn)[23-24]，可定義為任意有限數(shù)量的具有聯(lián)合高斯分布的隨機(jī)變量的集合，即高斯過程完全由均值函數(shù)m(x)和協(xié)方差函數(shù)kf(x,x′)確定,表達(dá)式為

f(x)～GP[m(x),k(x,x′)]

(7)

對(duì)于回歸問題，設(shè)輸出為隱式函數(shù)f(x)與高斯噪聲ε的和，即

y=f(x)+ε

(8)

(9)

對(duì)于與訓(xùn)練集x具有相同高斯分布的新數(shù)據(jù)集x*，y與預(yù)測(cè)值y*的聯(lián)合先驗(yàn)分布為

(10)

采用貝葉斯公式可導(dǎo)出相應(yīng)的后驗(yàn)分布為

(11)

式(11)中：

(12)

cov(y*)=Kf(x*,x*)-Kf(x,x*)T·

(13)

常用的協(xié)方差函數(shù)有平方指數(shù)協(xié)方差函數(shù)，表達(dá)式為

(14)

式(14)中：σf為函數(shù)變化垂直尺度的信號(hào)方差；σn為噪聲方差；l為長度尺度。超參數(shù)集θ[σf,l,σn]可由最大似然法確定[25]。

2 巖體力學(xué)參數(shù)反演的HBA-GPR-FLAC3D法

2.1 基本思路

優(yōu)化反演就是通過迭代計(jì)算搜索最佳參數(shù)組合，使位移計(jì)算值與位移實(shí)測(cè)值盡可能接近的方法，為此可設(shè)置反演目標(biāo)函數(shù)[24,26]的表達(dá)式為

(15)

研究發(fā)現(xiàn)，采用HBA進(jìn)行優(yōu)化時(shí)，與挖掘階段的全局尋優(yōu)相比，尋蜜階段的局部尋優(yōu)的適應(yīng)度函數(shù)評(píng)價(jià)次數(shù)明顯較多，因此，提升反演效率的關(guān)鍵在于如何顯著降低局部尋優(yōu)過程中的適應(yīng)度函數(shù)評(píng)價(jià)次數(shù)，即減小FLAC3D調(diào)用次數(shù)。為此，在HBA全尋優(yōu)過程中，將有最優(yōu)潛力的局部區(qū)域內(nèi)的算子信息作為訓(xùn)練樣本訓(xùn)練GPR，可獲得目標(biāo)函數(shù)在某個(gè)局部區(qū)域范圍內(nèi)的代理模型，為了加快局部尋優(yōu)，將GPR代理模型作為適應(yīng)度評(píng)價(jià)工具，采用牛頓法進(jìn)行局部尋優(yōu)，獲得當(dāng)前最優(yōu)算子的預(yù)測(cè)值，并調(diào)用FLAC3D計(jì)算獲得當(dāng)前最優(yōu)算子的真實(shí)適應(yīng)度。這樣處理能顯著減少局部尋優(yōu)過程中的FLAC3D調(diào)用次數(shù)，有效提高算法的尋優(yōu)效率。

提出的巖體力學(xué)參數(shù)反演的HBA-GPR-FLAC3D方法基本思路如下：首先，隨機(jī)生成攜帶巖土力學(xué)參數(shù)的初始化算子，并通過有限差分軟件FLAC3D進(jìn)行目標(biāo)函數(shù)的適應(yīng)度評(píng)價(jià)，得到初始算子的適應(yīng)度值；隨后，利用HBA的優(yōu)化策略對(duì)算子進(jìn)行不斷的迭代優(yōu)化，同時(shí)通過FLAD3D對(duì)算子進(jìn)行適應(yīng)度評(píng)價(jià)，在優(yōu)化的過程中，判定當(dāng)前最優(yōu)算子是否滿足收斂準(zhǔn)則，若滿足，則停止迭代并輸出結(jié)果；反之，則繼續(xù)進(jìn)行下一步迭代。

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

2.2.1 采用GPR建立當(dāng)前最優(yōu)算子局部鄰域的近似目標(biāo)函數(shù)

蜜獾優(yōu)化算法首先在模型空間內(nèi)產(chǎn)生初始化算子并進(jìn)行真實(shí)的目標(biāo)函數(shù)適應(yīng)度評(píng)價(jià)，經(jīng)過幾輪全局尋優(yōu)組建尋優(yōu)經(jīng)驗(yàn)數(shù)據(jù)集。之后，選取當(dāng)前最優(yōu)算子較近的多個(gè)算子構(gòu)建成局部尋優(yōu)訓(xùn)練數(shù)據(jù)集，采用GPR對(duì)數(shù)據(jù)集中的樣本進(jìn)行擬合，從而建立表達(dá)決策變量與目標(biāo)函數(shù)值之間的映射關(guān)系的代理模型。

2.2.2 針對(duì)代理模型的局部尋優(yōu)加速

采用牛頓法基于GPR代理模型所求取的一階導(dǎo)數(shù)和二階導(dǎo)數(shù)對(duì)局部最優(yōu)解進(jìn)行擬合，從而求取局部最優(yōu)解。在全局尋優(yōu)過程中，將尋優(yōu)過程中的局部最優(yōu)解代入真實(shí)適應(yīng)度函數(shù)進(jìn)行一次真算，再將其與當(dāng)前最優(yōu)解的真實(shí)適應(yīng)度函數(shù)進(jìn)行比較，如優(yōu)于當(dāng)前最優(yōu)解，則用其替代當(dāng)前最優(yōu)解，并標(biāo)記周邊區(qū)域?yàn)楫?dāng)前最優(yōu)區(qū)域，從而為下一迭代步的局部尋優(yōu)提供指引。

2.2.3 動(dòng)態(tài)更新局部尋優(yōu)數(shù)據(jù)集

在全局尋優(yōu)過程中，構(gòu)建一個(gè)包含所有算子信息的歷史數(shù)據(jù)集，該數(shù)據(jù)集記錄算子在尋優(yōu)過程中所有的歷史位置與相應(yīng)的適應(yīng)度值。在局部尋優(yōu)前，以該歷史數(shù)據(jù)集為基礎(chǔ)建立尋優(yōu)數(shù)據(jù)集。提出了一種動(dòng)態(tài)更新尋優(yōu)數(shù)據(jù)集的方式，以歷史最優(yōu)算子為中心，選取離中心最近的蜜獾算子構(gòu)建局部尋優(yōu)數(shù)據(jù)集，如圖3所示。通過保持?jǐn)?shù)據(jù)集的蜜獾算子總數(shù)不變，適時(shí)調(diào)整中心點(diǎn)位置，從而提高GPR代理模型在尋優(yōu)過程中的預(yù)測(cè)能力。

圖3 局部尋優(yōu)數(shù)據(jù)集的二維模型Fig.3 A 2D model for local optimization datasets

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

具體實(shí)現(xiàn)步驟如圖4所示。

圖4 HBA-GPR方法流程圖Fig.4 Flowchart of HBA-GPR method

步驟1確定待反演的巖體力學(xué)參數(shù)及其取值范圍，通過FLAC3D建立隧洞的數(shù)值計(jì)算模型以及選取相應(yīng)的監(jiān)測(cè)點(diǎn)，將目標(biāo)函數(shù)作為適應(yīng)度評(píng)價(jià)函數(shù)。

步驟2設(shè)置HBA與GPR算法的初始參數(shù)；種群數(shù)量為N，在求解空間中隨機(jī)選取N個(gè)空間位置作為算子的初始坐標(biāo)，代入FLAC3D數(shù)值計(jì)算軟件，計(jì)算監(jiān)測(cè)點(diǎn)的位移值并代入目標(biāo)函數(shù)，獲得相應(yīng)算子的真實(shí)適應(yīng)度。

步驟3根據(jù)HBA算法進(jìn)行全局尋優(yōu)，獲得當(dāng)前最優(yōu)算子；將全局尋優(yōu)過程中所有算子的坐標(biāo)及其真實(shí)適應(yīng)度值保存到歷史數(shù)據(jù)集中；判斷是否滿足收斂準(zhǔn)則，若滿足，則停止迭代并轉(zhuǎn)至步驟8，否則，轉(zhuǎn)至步驟3。

步驟4以當(dāng)前最優(yōu)算子作為中心點(diǎn)，從歷史數(shù)據(jù)集中選取歐氏距離該中心點(diǎn)最近的3N個(gè)算子構(gòu)建成局部尋優(yōu)訓(xùn)練樣本集。

步驟5采用局部尋優(yōu)訓(xùn)練樣本集訓(xùn)練GPR代理模型。

步驟6求解GPR代理模型在當(dāng)前最優(yōu)算子點(diǎn)上一階和二階導(dǎo)數(shù)，采用牛頓法進(jìn)行局部尋優(yōu)，獲得預(yù)測(cè)的局部最優(yōu)算子。

步驟7將預(yù)測(cè)的局部最優(yōu)算子坐標(biāo)代入FLAC3D，計(jì)算監(jiān)測(cè)點(diǎn)位移值，獲得相應(yīng)的真實(shí)適應(yīng)度，并添入歷史數(shù)據(jù)集；若預(yù)測(cè)的局部最優(yōu)算子真實(shí)適應(yīng)度小于當(dāng)前最優(yōu)算子的真實(shí)適應(yīng)度，則將預(yù)測(cè)的局部最優(yōu)算子替代當(dāng)前最優(yōu)算子，否則，當(dāng)前最優(yōu)算子不變。

步驟8判斷當(dāng)前最優(yōu)算子的真實(shí)適應(yīng)度是否滿足收斂準(zhǔn)則，若滿足，則停止迭代并轉(zhuǎn)步驟8，否則，轉(zhuǎn)至步驟3。

步驟9輸出當(dāng)前最優(yōu)算子及其真實(shí)適應(yīng)度，結(jié)束計(jì)算。

3 算例驗(yàn)證

隧洞半徑為3.5 m，隧洞沿坐標(biāo)軸x、z方向圍巖厚度均為16.5 m，沿y方向取單位厚度，均質(zhì)巖性，隧洞的埋深為50 m，圍巖容重為20 kN/m3，初始地應(yīng)力為10 MPa，數(shù)值計(jì)算模型如圖5所示。該模型采用FLAC3D軟件推薦的8節(jié)點(diǎn)單元，共有1 000個(gè)單元體，1 386個(gè)節(jié)點(diǎn)。越靠近開挖面，網(wǎng)格尺度越小。目標(biāo)函數(shù)設(shè)為

(16)

圖5 FLAC3D計(jì)算網(wǎng)格劃分圖 Fig.5 Computational grid of FLAC3D

為驗(yàn)證算法可行性，假設(shè)真實(shí)的巖體力學(xué)參數(shù)為E=2.2 GPa,c=1.1 MPa,φ=30°。沿著隧洞環(huán)向方向分別在0°、22.5°、45°、67.5°、90°處選取5個(gè)監(jiān)測(cè)點(diǎn)(如圖6所示)，采用FLAC3D計(jì)算監(jiān)測(cè)點(diǎn)的水平和垂直位移，將其視為“實(shí)測(cè)位移”，如表1所示。

圖6 隧洞監(jiān)測(cè)點(diǎn)布置圖Fig.6 Layout of monitoring points

表1 各監(jiān)測(cè)點(diǎn)的實(shí)測(cè)位移值Table 1 The measured displacement of points

采用經(jīng)典的 Mohr-Colomb 彈塑性本構(gòu)模型?；贛ohr-Colomb破壞準(zhǔn)則的理想彈塑性本構(gòu)模型是巖石工程領(lǐng)域應(yīng)用最廣泛的本構(gòu)模型之一。該模型輸入?yún)?shù)少，能很好地描述地下工程開挖后的圍巖力學(xué)，可根據(jù)數(shù)值計(jì)算的網(wǎng)格單元是否進(jìn)入塑性狀態(tài)來判斷巖體是否發(fā)生斷裂。

在FLAC3D自帶的FISH語言系統(tǒng)中，位移邊界的約束主要依賴于速度的約束[27]，當(dāng)單元或節(jié)點(diǎn)的速度被約束為0時(shí)，位移即被約束。在該隧洞模型中，右側(cè)和底部節(jié)點(diǎn)的位移自由度在縱向和豎直方向上都被約束住，以保證模型能順利收斂。模型的應(yīng)力邊界模擬重力方向的地應(yīng)力，并通過梯度應(yīng)力施加到模型上。

分別采用HBA-FLAC3D法與HBA-GPR-FLAC3D法反演該隧洞的巖體力學(xué)參數(shù)，兩種方法中的HBA參數(shù)設(shè)置均相同：N=50，β=6，C=2。HBA-GPR-FLAC3D法中的GPR模型初始超參數(shù)均設(shè)置為：lnl=[-1,-1,-1]，lnσf=-1，lnσn=ln(1×10-6)。

基于該隧洞的實(shí)測(cè)巖體力學(xué)參數(shù)，將反演參數(shù)的搜索范圍設(shè)為：2.0 GPa≤E≤2.4 GPa，1.0 MPa≤c≤1.4 MPa，28°≤φ≤32°，均以目標(biāo)函數(shù)f(X)<5×10-2這一條件作為收斂準(zhǔn)則。為了確保上述兩種方法對(duì)比結(jié)果的可信性，降低偶然性的影響，每種算法均進(jìn)行10次獨(dú)立運(yùn)算，取平均適應(yīng)度評(píng)價(jià)次數(shù)作為每種算法的最終結(jié)果。

圖7 地下廠房分期開挖示意圖Fig.7 Schematic diagram of excavation by stages of the power house

計(jì)算結(jié)果如表2所示，與HBA-FLAC3D法相比，HBA-GPR-FLAC3D法的反演精度較高，但FLAC3D調(diào)用次數(shù)與計(jì)算耗時(shí)明顯較少，由此說明本文方法的可行性。

表2 兩種算法計(jì)算結(jié)果對(duì)比Table 2 The results calculated by the two methods

4 工程應(yīng)用

某水電站地下廠房洞室群主要包括主廠房、母線洞、主變室以及尾閘室等大型洞室，圍巖主要為Ⅱ、Ⅲ類花崗巖，布置有四臺(tái)水輪發(fā)電機(jī)組。主廠房尺寸為180 m×25.9 m×53.675 m，主變室尺寸為164 m×17.5 m×18.175 m，母線洞長度為35 m，分期開挖情況如圖7所示。地下廠房共設(shè)置了兩個(gè)位移監(jiān)測(cè)斷面，相應(yīng)監(jiān)測(cè)點(diǎn)布置如圖8所示。選取經(jīng)典彈塑性模型作為其本構(gòu)模型，根據(jù)該本構(gòu)模型可確定待反演的巖體力學(xué)參數(shù)包括彈性模量E、黏聚力c和內(nèi)摩擦角φ。開挖體的數(shù)值計(jì)算網(wǎng)格如圖9所示，數(shù)值計(jì)算區(qū)域?yàn)椋貉豿軸950 m，沿y軸650 m，沿z軸-6～472 m。數(shù)值計(jì)算網(wǎng)格共使用了687 133個(gè)單元和191 323個(gè)節(jié)點(diǎn)，網(wǎng)格通過不同顏色以區(qū)分隧道面積和開挖順序。由于地下洞室圍巖基本為花崗巖，在三維彈塑性數(shù)值計(jì)算中采用彈脆塑性本構(gòu)模型[28]。計(jì)算采用Mohr-Coulomb強(qiáng)度準(zhǔn)則，在建模過程中強(qiáng)度參數(shù)(黏聚力c和摩擦角φ)被動(dòng)態(tài)削弱。

圖8 位移監(jiān)測(cè)斷面Fig.8 The layout map of displacement monitoring points section

圖9 開挖體計(jì)算網(wǎng)格Fig.9 The computing grid of excavated rock

相應(yīng)的監(jiān)測(cè)點(diǎn)布置如圖10所示，選擇M7、M19以及M24監(jiān)測(cè)點(diǎn)的實(shí)測(cè)位移構(gòu)建目標(biāo)函數(shù)，目標(biāo)函數(shù)見式(16)。待反演的巖體力學(xué)參數(shù)有：主廠房上游圍巖彈性模量E1、主廠房1#機(jī)組下游圍巖彈性模量E2、主廠房2# ～ 4#機(jī)組下游圍巖彈性模量E3、主變室圍巖彈性模量E4、主廠房上游圍巖黏聚力c1、主廠房下游圍巖黏聚力c2以及內(nèi)摩擦角φ。參數(shù)搜索區(qū)間設(shè)為：60 GPa≤E1≤80 GPa，10 GPa≤E2≤30 GPa，20 GPa≤E3≤ 40 GPa，100 GPa≤E4≤300 GPa，5 MPa≤c1≤7 MPa，5 MPa≤c2≤7 MPa，40°≤φ≤60°。

圖10 FLAC3D模型中的監(jiān)測(cè)點(diǎn)Fig.10 Monitoring points in the FLAC3D model

經(jīng)過前期調(diào)試工作確定了一種參數(shù)組合作為目標(biāo)函數(shù)的最優(yōu)解：HBA參數(shù)設(shè)置為N=50，β=6，C=2；GPR模型的初始超參數(shù)設(shè)置為lnl=[-1,-1,-1,-1,-1,-1,-1]，lnσf=-1，lnσn=ln(1×10-6)。以FLAC3D調(diào)用次數(shù)大于等于1 000作為收斂準(zhǔn)則。需要指出的是，所提出的模型參數(shù)是適用于彈脆塑性模型的一種參數(shù)組合，如果使用另一種本構(gòu)模型，參數(shù)的取值將會(huì)相應(yīng)發(fā)生變化，以保證結(jié)果的準(zhǔn)確性。

分別采用HBA-FLAC3D法與本文方法反演巖體力學(xué)參數(shù)，結(jié)果如表3、表4所示。對(duì)于M7、M19監(jiān)測(cè)點(diǎn)，HBA-GPR-FLAC3D方法獲得的計(jì)算位移值與實(shí)測(cè)位移值的相對(duì)誤差僅為1.44%、2.65%，而HBA-FLAC3D方法的相對(duì)誤差為7.81%、15.72%；對(duì)于監(jiān)測(cè)點(diǎn)M24，兩種方法的相對(duì)誤差分別為0.06%、0.73%，差別不大。由此說明，在FLAC3D調(diào)用次數(shù)達(dá)1 000次的相同條件下，本文方法的尋優(yōu)精度明顯高于HBA-FLAC3D，說明了本文方法的先進(jìn)性。

表3 兩種方法的巖體力學(xué)參數(shù)反演結(jié)果

表4 計(jì)算位移與實(shí)測(cè)位移對(duì)比Table 4 Comparison of displacements obtained from calculation and measurement

5 結(jié)論

針對(duì)復(fù)雜地下工程巖體力學(xué)參數(shù)優(yōu)化反演中獲取全局最優(yōu)解與計(jì)算耗時(shí)之間的矛盾，提出了一種新的巖體力學(xué)參數(shù)優(yōu)化反演方法，得到以下結(jié)論。

(1)該方法將強(qiáng)全局尋優(yōu)能力的仿生優(yōu)化算法HBA、強(qiáng)擬合能力的GPR機(jī)器學(xué)習(xí)模型以及FLAC3D數(shù)值計(jì)算有機(jī)結(jié)合，具有非侵入性、實(shí)現(xiàn)簡(jiǎn)單、易于工程師掌握的特點(diǎn)。

(2)在圓形隧洞算例和水電站地下廠房實(shí)例中，分別運(yùn)用HBA-FLAC3D算法以及本文提出的HBA-GPR-FLAC3D算法進(jìn)行巖體力學(xué)參數(shù)的反演計(jì)算，結(jié)果表明，本文方法在精度和效率上都有著明顯優(yōu)勢(shì)，適用于單次數(shù)值計(jì)算耗時(shí)大的復(fù)雜地下工程巖體力學(xué)參數(shù)反演問題，為地下工程巖體力學(xué)參數(shù)快速確定提供了一條有效途徑。

由于本文方法的優(yōu)化結(jié)果與蜜獾優(yōu)化算法、高斯過程回歸模型中參數(shù)的設(shè)置有著直接關(guān)系，因此在今后的研究中將探索如何選擇合適的參數(shù)。