基于XGBoost的初始地應力場修正反演法

劉敬文，白金朋，王 建，徐 韌，徐偉男

(1.河海大學水利水電學院，江蘇 南京 210098； 2.中國電建集團北京勘測設計研究院有限公司，北京 100024； 3. 無錫市惠山區(qū)水利局，江蘇 無錫 214174)

初始地應力場是地下廠房等工程設計與施工所需考慮的重要因素之一[1]，然而，受場地、經費等現(xiàn)實條件的限制，不可能大規(guī)模開展原位地應力測量，只能利用有限的實測地應力樣本來反演廠房區(qū)的地應力分布，這是目前確定初始地應力場的主要途徑。

常用的反演方法主要有多元線性回歸分析法[2-3]和智能優(yōu)化算法[4-6]等，其中，多元線性回歸分析法采用最小二乘法求解回歸系數且解唯一，故被廣泛采用[7]。但是，多元線性回歸分析法是基于彈性力學的疊加原理，更適用線彈性材料的情況。對于復雜地質條件下的巖體應力場，部分巖體結構不滿足彈性假設，力學非線性較為顯著，需要開展修正以解決線性回歸可能存在的擬合精度問題。此前已有部分學者研究過修正、優(yōu)化反演的方法，郭運華等[8]在最小二乘回歸的基礎上，根據屈服條件修正偏應力大小，對不滿足強度條件的單元地應力進行修正，以滿足彈性假設條件。郭明偉等[9]采用優(yōu)化位移邊界條件擬合工程區(qū)域初始地應力場，通過調整位移邊界模式得到其最優(yōu)組合并依此擬合工程區(qū)域的初始地應力場。閆相禎等[10]利用阻尼最小二乘法建立了應力場反演的優(yōu)化約束模型，提出了一種用于反演地應力場邊界力的優(yōu)化分析方法。然而，上述方法大都從數值模擬的角度出發(fā)，通過改變模型邊界的約束參量實現(xiàn)修正反演的目的。為了進一步提高反演計算精度，本文在有限元多元線性回歸分析的基礎上，采用極端梯度提升(extreme gradient boosting,XGBoost)算法，以不同構造應力下有限元的計算結果為特征函數，以應力實測值與線性回歸結果的差值為目標函數，對多元線性回歸的反演結果進行修正，其核心思想是基于實測應力場與線性回歸應力場的差值構造一個殘差修正場，這樣做一方面有效減小了使用線彈性模型回歸計算帶來的誤差，同時避免了直接使用XGBoost純數學算法擬合應力場，更加準確擬合實際應力場。

XGBoost算法是一種基于Boosting的集成學習算法，該算法具有計算復雜度低、運行速度快、精度高等優(yōu)點。由于迭代初始值采用單位載荷法，XGBoost有效地避免了地層幾何構造和邊界條件的限制，同時，對不滿足彈性假設的一些單元的地應力進行了修正，使其更符合實際情況。因此，使用XGBoost修正反演結果，一方面可以有效修正使用線彈性模型回歸計算帶來的應力偏差，提高反演結果的精度，另一方面利用該集成算法所得的黑箱模型可以有效降低誤差、增強數據可靠性，從而為提高地應力場反演計算的精度和數據的可靠性提供了一種新的解決方案。

1 修正反演分析原理和方法

1.1 多元線性回歸法

(1)

根據最小二乘法原理，只需求得使殘差平方和為最小值時的n個待定回歸系數L=(L1,L2,…,Ln)T，計算域內任一點P的回歸初始地應力就可由該點各工況有限元法計算值迭加而得

(2)

式中：j=1,2,…,n對應n個初始地應力分量。

1.2 XGBoost算法修正反演原理與模型構建

XGBoost算法是梯度提升機器(gradient boosting machine，GBM)算法的擴展，是一種同時具備非線性模型和樹模型特點的優(yōu)化模型，能夠同時完成回歸和分類任務，且已在疾病診斷、交通流量預測、商品銷售量預測等領域的應用取得成功[11]，目前在水利及巖土領域應用較少。XGBoost算法由多棵決策樹(cart)組成，通過決策樹集成實現(xiàn)機器學習，所有決策樹的回歸值相累加即為模型回歸值。與一般算法不同的是，XGBoost算法通過對損失函數二階泰勒展開以快速逼近目標函數，并加入了正則項，調節(jié)參數避免過擬合，訓練速度快，擬合及預測精度高[12]。

在構建地應力回歸的修正反演模型時，XGBoost算法通過在回歸函數中不斷加入關于各因子的新函數來逼近實測值，以解決地應力回歸時非線性擬合不足的問題，即：

(3)

XGBoost算法的目標函數定義如下：

(4)

式中：C為損失函數，用于評估回歸值與真實值之間的損失；Ω(fk)為正則化函數，表示第k棵決策樹的復雜度，用于控制模型計算效率，避免過擬合。正則化函數定義如下：

(5)

式中：γ、λ為正則項懲罰系數；T為該樹的葉子節(jié)點個數；ω為葉子節(jié)點對應的權重。

可以對式(4)改寫并按二階泰勒級數展開，進一步簡化后得到：

式中：gi為損失函數的一階梯度統(tǒng)計；hi為損失函數的二階梯度統(tǒng)計。通過迭代找到f(x)使損失函數最小則模型訓練完成。

綜上所述，本次XGBoost算法修正反演的建模流程如圖1所示，具體步驟如下：

圖1 基于XGBoost算法修正反演的計算流程

步驟1根據地質條件，對廠房及周邊區(qū)域進行有限元建模。開展基于實測地應力的反演分析，對模型進行各構造模式對應的加載和約束，獲取各構造模式下的應力場，并利用多元線性回歸分析反演得到整個廠房及周邊區(qū)域的初始地應力分布。

步驟2提取有限元計算的應力結果和多元線性回歸值的殘差，對這些數據進行處理，整理成合適的影響因子，共同組成XGBoost輸入樣本的數據集，并劃分出訓練集和測試集。對模型進行參數調優(yōu)，選取最優(yōu)參數結果構建基于XGBoost的修正反演模型。

步驟3將測試集中的數據代入步驟2中構建的模型，并進行交叉驗證，分別比較多元線性回歸值、模型修正反演值與實測值的均方根誤差RMSE和相對誤差均值RSE，對XGBoost修正反演模型做出評價。最終利用該模型對廠房及周邊區(qū)域內的初始地應力回歸值進行修正，得到修正后的整體初始地應力場。

2 工程實例

2.1 工程概況及地質條件

某抽水蓄能電站位于遼寧省大連市境內，其上水庫為東、北、西三面環(huán)山的天然庫盆地形，主壩修筑于庫盆南側的主溝內，副壩位于庫盆東北側低矮分水嶺處，均為瀝青混凝土面板堆石壩；輸水發(fā)電系統(tǒng)深埋于上、下水庫間山體內，圍巖以微風化～新鮮變質長石砂巖為主。

在進行地應力回歸分析時，計算模型的范圍為：上、下游邊界以地下廠房為中心，向南北方向各延伸400 m；左、右邊界以地下廠房為中心，向東西方向各延伸500 m；上、下邊界上至地表，下至高程-170 m。

初始地應力場反演區(qū)域和三維數值模型如圖2所示。模型主要考慮fp25、fp28、fp31、fp32等較大的斷層，有限元模型共有節(jié)點169 349個，單元967 079個。本次計算將斷層簡化為勻質連續(xù)體，采用彈性模型模擬巖體及斷層的力學行為，涉及的力學參數包括彈性模量、泊松比和密度等。本次計算采用的巖體和斷層力學參數見表1。

圖2 地應力反演三維有限元模型

某抽水蓄能電站在地下廠房斜井段(鉆孔ZK213)、高壓岔管段(鉆孔ZK201)、下平段(鉆孔ZK204)、廠房區(qū)(鉆孔ZK202、ZK206)、地表孔ZK149共選定6個鉆孔(每個鉆孔有5個測點)內進行了水壓致裂法地應力測試。鉆孔位置見圖3，地應力

表1 計算巖體和斷層力學參數

試驗實測結果匯總表見表2，表中SH為水平方向最大壓應力；Sh為水平方向最小壓應力；SV為豎向應力；負號表示為壓力。

圖3 地應力測孔位置示意圖

表2 實測地應力

2.2 基于多元線性回歸的一次反演

根據實測地應力的結果，采用多元線性回歸進行初始地應力場的反演，通過模擬若干組構造模式，并對其進行線性組合，從而獲得最終地應力。由于水壓致裂法無法測得鉛錘面內的剪切應力，一般采取以下5種構造模式：①自重；②x方向均勻分布、三角形分布的擠壓構造力；③y方向均勻分布、三角形分布的擠壓構造力；④xy向的剪切應力；⑤yx向的剪切應力。由這5種構造模式進行多元線性回歸計算，以確定自重系數L1，x方向水平構造系數L2、L3，y方向水平構造系數L4、L5，xy向剪切應力系數L6、yx向剪應力系數L7。將32個測點的192組應力數據代入公式進行計算，求得回歸系數結果如下：L1=1.189、L2=1.498、L3=0.553、L4=1.322、L5=-0.526、L6=1.652、L7=-1.470。復相關系數r=0.974，回歸效果較好。偏相關系數V1=0.965，V2=0.915，V3=0.364，V4=0.940，V5=0.440，V6=0.701，V7=0.752。其中V1、V2、V4較大，說明所施加的自重應力場、x方向均布荷載、y方向均布荷載作用較顯著。

求得回歸系數后，代入式(2)求得在實測位置的回歸計算值，部分回歸結果與實測值的比較見表3(表中數據壓應力為負，拉應力為正)。雖然本次線性回歸反演的復相關系數較高，但標準差卻高達0.97 MPa，說明反演擬合精度存在一定的不足，需要進行修正。

表3 部分測點的應力分量實測值與回歸值對比 單位：MPa

2.3 基于XGBoost的修正反演

2.3.1數據預處理

將2.2中的每一種構造模式代入有限元，作為邊界條件進行正演計算，得到有限元計算結果，并提取L1～L7七個系數分別對應工況下的應力因子，作為輸入變量；將2.2中多元線性回歸結果中應力分量的殘差作為輸出變量。

輸入樣本總數為180組，將輸入樣本集的前150組劃分為訓練集，后30組劃分為測試集。

2.3.2模型的參數調優(yōu)

對各關鍵參數進行初值定義，由于模型參數max_depth(最大樹深)對突變數據較為敏感，對模型收斂的穩(wěn)定性影響較大，因此針對max_depth進行敏感性分析，以權衡參數的最優(yōu)取值。

將max_depth的優(yōu)化范圍設為(2，8)，比較和權衡不同樹深取值下模型的收斂速度與精度，調參結果顯示，當最大樹深max_depth=4時，模型的收斂速度與精度都較高，收斂穩(wěn)定性較好，故本次最大樹深參數選定為4，其余參數為默認值。

模型設置最大迭代次數為250次，最多經過190次迭代即可收斂。圖4為模型訓練迭代過程。

圖4 模型訓練迭代過程

2.3.3模型的結果與評價

求得測試集的修正反演值后，與實測值、回歸值的對比如表4所示。

由表4可知：多元線性回歸結果與實測值之間的均方根誤差為1.13 MPa，XGBoost修正反演結果與實測值之間的均方根誤差為0.48 MPa，擬合精度提高約60%；從相對誤差的平均值看，XGBoost修正反演精度也得到不同程度提高。可見XGBoost修正反演的結果更精確且更具有代表性，對多元線性回歸的反演結果進行了二次優(yōu)化，同時，也在一定程度上減少了非協(xié)調回歸因子的影響，在多元線性回歸的基礎上提升了反演精度以及數據的可靠性。

為探究本文方法在工程實際中樣本量相對較少時的反演精度，以及評估模型的質量、避免過擬合與欠擬合現(xiàn)象的發(fā)生，采用兩種方法加以驗證：

a.保持原先30組測試集不變，將訓練集數量由150組減少到60組，其余參數不變，得到修正反演結果與實測值之間的均方根誤差為0.54 MPa，相比于多元線性回歸結果的均方根誤差1.13 MPa，擬合精度提高約50%。由此可見，當樣本數量較少時，該修正反演方法的精度依舊較高，效果良好。

表4 鉆孔測試集主應力的反演分析值對比

b.保持訓練集與測試集總數不變，重新劃分一組數據集，進行交叉驗證。取樣本集的后150組為訓練集，前30組為測試集，其余參數不變，得到修正反演結果與實測值之間的均方根誤差為0.65 MPa，相比于多元線性回歸結果的均方根誤差1.56 MPa，擬合精度提高約58%，說明模型穩(wěn)定性高、泛化能力好，能夠很好地滿足初始地應力修正反演的工作需求。

本文模型基于Python3.7進行，計算機內存12.0 GB、4核CPU，計算967 079個單元的初始應力共計用時948.5 s，反演效率高，滿足工程需要。

圖5和圖6分別為主廠房軸線、垂直主廠房軸線剖面的第一、第三主應力等值線云圖。由圖5和圖6可知，主廠房開挖區(qū)域的主應力從上到下基本遞增，等值線在淺層受地形變化影響較大，說明修正反演計算的應力場充分反映了地形、地貌的影響。

圖5 主廠房軸線剖面第一主應力云圖

圖6 主廠房軸線剖面第三主應力云圖

3 結 語

在地應力反演的傳統(tǒng)多元線性回歸方法基礎上，引入XGBoost算法并基于實測應力場與線性回歸應力場的差值構造殘差修正場，對多元線性回歸的反演結果進行了二次優(yōu)化，既提高了反演精度，同時也避免了直接使用純數學的XGBoost算法來擬合應力場物理概念不明確、容易產生過擬合的問題，更加符合力學原理。

工程實例計算表明，采用本文方法，擬合精度提高約60%，且計算穩(wěn)定性高、泛化能力好，能夠很好地滿足初始地應力修正反演的需求。