巖體力學參數(shù)反演的代理蜜獾優(yōu)化方法

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

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

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

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

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

1 HBA和GPR的基本原理

1.1 HBA原理簡介

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

圖1 蜜獾的覓食行為特點Fig.1 Characteristics of foraging behavior of honey badger

蜜獾算子的位置表達式為

xi=lbi+r1(ubi-lbi)

(1)

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

(2)

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

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

(3)

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

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

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

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

(4)

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

(5)

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

xnew=xprey+Fr7αdi

(6)

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

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

1.2 GPR原理簡介

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

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

(7)

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

y=f(x)+ε

(8)

(9)

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

(10)

采用貝葉斯公式可導出相應的后驗分布為

(11)

式(11)中：

(12)

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

(13)

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

(14)

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

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

2.1 基本思路

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

(15)

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

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

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

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

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

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

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

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

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

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

2.3 實現(xiàn)步驟

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

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

步驟1確定待反演的巖體力學參數(shù)及其取值范圍，通過FLAC3D建立隧洞的數(shù)值計算模型以及選取相應的監(jiān)測點，將目標函數(shù)作為適應度評價函數(shù)。

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

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

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

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

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

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

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

步驟9輸出當前最優(yōu)算子及其真實適應度，結(jié)束計算。

3 算例驗證

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

(16)

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

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

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

表1 各監(jiān)測點的實測位移值Table 1 The measured displacement of points

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

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

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

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

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

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

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

4 工程應用

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

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

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

相應的監(jiān)測點布置如圖10所示，選擇M7、M19以及M24監(jiān)測點的實測位移構(gòu)建目標函數(shù)，目標函數(shù)見式(16)。待反演的巖體力學參數(shù)有：主廠房上游圍巖彈性模量E1、主廠房1#機組下游圍巖彈性模量E2、主廠房2# ～ 4#機組下游圍巖彈性模量E3、主變室圍巖彈性模量E4、主廠房上游圍巖黏聚力c1、主廠房下游圍巖黏聚力c2以及內(nèi)摩擦角φ。參數(shù)搜索區(qū)間設為：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)測點Fig.10 Monitoring points in the FLAC3D model

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

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

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

表4 計算位移與實測位移對比Table 4 Comparison of displacements obtained from calculation and measurement

5 結(jié)論

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

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

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

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