基于殘差平方和取舍常數項的地應力場反演分析

蒙 偉，何 川，嚴 健，吳枋胤，周子寒，寇 昊

(西南交通大學 交通隧道工程教育部重點實驗室，四川 成都 610031)

近年來，隨著復雜艱險山區高速公路、鐵路等重大工程和西部大開發相繼實施，在這些區域逐漸涌現出大量長大、深埋隧道等地下工程[1]。該區域大量地應力資料表明，豎直地應力小于水平地應力的現象非常普遍，即中國西部地區普遍存在擠壓構造運動[2-3]。故對此區域隧道等地下工程巖體的初始地應力場進行反演分析具有重要的實際意義[4]。

在反演過程中，郭懷志等[5]給出了能夠獲得唯一解的反演分析方法，故此方法被廣泛地沿用至今。文獻[5-9]、[10-14]對巖體的初始地應力場進行反演分析采用的回歸模型分別為

(1)

(2)

由式(1)、式(2)可知，依據初始地應力場回歸模型是否含有常數項，回歸模型可分為兩類：①不含常數項的回歸模型;②含常數項的回歸模型。然而，截至目前，在巖體初始地應力場回歸反演過程中，含與不含常數項回歸模型被單獨地用于反演巖體的初始地應力場。由于這兩類回歸模型的反演結果不相同，因此，僅采用其中一類回歸模型進行反演會存在一個問題，即在這兩類回歸模型中，反演結果較差的回歸模型很可能被用于反演。為解決此問題，提高巖體初始地應力場的反演結果，在反演過程中，有必要對如何選擇最佳反演結果的回歸模型進行分析研究。

本文為客觀地選擇最佳反演結果的回歸模型，結合概率統計理論提出了衡量巖體初始地應力場反演結果的量化指標。然后在郭懷志教授提出的反演分析方法的基礎上，在含與不含常數項回歸模型中，基于量化指標提出了客觀獲得最佳反演結果的反演流程。最后，以某復雜艱險山區鐵路隧道為依托對反演流程進行工程應用分析。

1 取舍回歸模型常數項的反演機理

1.1 衡量反演結果的量化指標

由概率統計相關理論可知，殘差平方和可描述隨機因素影響引起的偏差，即殘差平方和可反映實測值與計算值之間的偏差程度。

(3)

假設y由一個自變量組成，樣本容量為5，展示含與不含常數項回歸模型實測值與計算值的偏差程度見圖1。

圖1 含與不含常數項回歸模型的偏差程度

由圖1可知，由于實測值是固定不變的，當用含與不含常數項回歸模型去匹配實測值時，實測值與計算值之間的偏差程度是不相同的。同理，在巖體初始地應力場回歸反演過程中，采用含與不含常數項回歸模型的反演結果也不相同，故有必要對如何選擇最佳反演結果的回歸模型進行分析研究。

通過對比式(3)與式(1)、式(2)可知，不含與含常數項回歸模型的應力實測值與計算值的偏差程度分別為

(4)

(5)

應力是矢量，除了大小還具有方向，即合理的巖體初始地應力場反演結果應滿足應力實測值與計算值在大小和方向上基本一致[9]。因此，需采用應力大小殘差平方和與應力方向殘差平方和兩個指標表示應力實測值與計算值之間的偏差程度。應力方向殘差平方和為

(6)

為綜合評價應力實測值與計算值的大小和方向偏差程度，需把應力大小殘差平方和與應力方向殘差平方和兩個指標放到一個計算式中，而應力大小殘差平方和的單位(Pa2或MPa2)與應力方向殘差平方和的單位[(°)2]不相同，如果采用應力大小殘差平方和加應力方向殘差平方和的形式則還需對其進行標準化處理。為方便工程人員以后應用此方法，故采取應力大小殘差平方和乘以應力方向殘差平方和的形式，計算式為

(7)

1.2 基于殘差平方和取舍常數項的反演流程

在反演過程中，首先同時采用含與不含常數項回歸模型進行反演，然后通過式(7)客觀地選擇其中最佳反演結果的回歸模型(較目前僅采用其中一類回歸模型的反演流程，可避免采用反演結果較差的回歸模型進行反演)，最后采用該回歸模型即可得到與實測原位地應力較為吻合的巖體初始地應力場。如果采用最小二乘法求解回歸系數，則具體的實施方案見流程圖2。

圖2 基于殘差平方和取舍常數項的地應力場反演流程

2 工程概況

某隧道全長約16 km，最大埋深1 350 m，為某新建復雜艱險山區鐵路隧道的關鍵控制性工程，屬于單洞單線隧道，處于中國地殼板塊構造運動最強烈的區域之一(青藏高原)。隧址區主要分布的是閃長巖等硬巖，彈性模量為40 GPa，隧址區屬構造剝蝕區，山勢陡峻，溝谷深切，最大切割深度為2 100 m，典型高山峽谷地貌，加上地殼板塊構造運動強烈，在隧道開挖過程中頻繁發生中等及強烈巖爆。

在隧道DK186+183—192附近伴有強烈巖爆發生，為研究該區域附近復雜的初始地應力場，本文以DK183+000—188+000深埋峽谷段為例進行反演分析，隧道穿越的巖體物理參數見表1。該區域通過水壓致裂法測得的原位地應力見表2[17]。此外，隧道還存在高地溫，記錄的最高地溫為89.5 ℃，且地溫高居我國之首[15-16]。

表1 巖體物理參數

表2 實測原位地應力

由表2可知，最大水平主應力大于垂直應力，故在該水壓致裂鉆孔附近存在強烈的水平構造運動，與隧道在拱頂附近發生巖爆現象基本一致，表明由于板塊構造運動，該鉆孔附近區域巖體的初始地應力場較自重應力場更加復雜，且DK-SZLSD-2鉆孔5個測點的原位地應力較難反映該局部區域巖體初始地應力場的宏觀分布規律，故有必要對該局部區域巖體的初始地應力場進行反演分析。

3 工程應用

3.1 建立三維數值模型

數值模型的網格及地形見圖3，共計2 853 904個單元，3 863 303個節點，模型底面為隧道軸線標高以下1 000 m，頂面為巖體實際地形。為方便獲取垂直于隧道軸線方向的應力，數值模型以隧道軸線方向為x軸，垂直隧道軸線方向為y軸，沿隧道埋深鉛直方向為z軸。

圖3 劃分網格后的數值模型(單位：m)

3.2 定義邊界條件

因原位地應力采用水壓致裂法進行測量，見表2。受其原理限制，所測得的原位地應力無鉛垂面內的剪切應力，故剪切構造荷載僅考慮水平xoy平面內的構造荷載[5, 11, 14]，施加的自重及構造荷載大小見圖4。施加方法：除數值模型的頂部邊界，其他非加載的邊界均施加法向位移零的約束；對圖4(a)數值模型施加豎向9.8 m/s2重力加速度；對圖4(b) 數值模型施加沿x軸的1 MPa水平均布荷載；對圖4(c) 數值模型施加沿y軸的1 MPa水平均布荷載；對圖4(d) 數值模型施加x方向的1 kPa/m水平三角形荷載；對圖4(e)數值模型施加y方向的1 kPa/m水平三角形荷載；對圖4(f)施加水平剪切荷載，依據文獻[18]提出的公式，可對該數值模型垂直于x軸的邊界施加5 000 m/100 000=0.05 m的剪切位移，對該數值模型垂直于y軸的邊界施加3 000 m/100 000=0.03 m的剪切位移，即可模擬得到該模型的水平剪切構造應力場。

圖4 荷載分布

3.3 轉換原位地應力坐標系

如圖3所示，為方便獲取垂直于隧道軸線方向的應力，數值模型以隧道軸線方向為x軸，因此數值模型采用的xoy坐標系與實測原位地應力采用的NoW坐標系不一致，故需要對實測原位地應力進行坐標變換。應力由舊坐標系轉換到新坐標系時，彈性力學應力轉軸式為

σu′v′=αu′uαv′vσuv

(8)

式中：σuv、σu′v′分別為轉換前后的應力張量；αu′u、αv′v為轉換系數。

含與不含常數項回歸模型的反演結果見表3。

表3 反演結果

3.4 回歸反演巖體初始地應力場

依據式(1)和式(2)的回歸模型，結合圖4的自重和構造荷載，以及溫度應力[16]，隧道采用的不含與含常數項回歸模型分別為

(9)

(10)

使應力實測值與計算值的殘差平方和最小即可求解得到回歸系數。在式(2)的含常數項回歸模型中，即需使式(5)最小即可，為使殘差平方和最小，則需求殘差平方和分別關于回歸系數的偏導數，并使偏導數為零，即有

s=1,2,…,n

(11)

化簡式(11)有

(12)

式(12)為通過最小二乘法求解回歸系數的正規方程組，為方便求解回歸模型的回歸系數，可將式(12)改寫成直觀的矩陣形式為

(13)

同理，對于式(1)的不含常數項回歸模型，通過最小二乘法求解回歸系數的矩陣形式正規方程組為

(14)

由式(13)、式(14)可知，在通過最小二乘法求解回歸系數時，所采用的正規方程組是不相同的。因此，在求解回歸系數的過程中需要注意，若采用含常數項回歸模型，則應采用式(13)的正規方程組求解回歸系數；若采用不含常數項回歸模型，則應采用式(14)的正規方程組求解回歸系數。而在目前的文獻中，大都選擇式(14)求解回歸系數，卻得到了含常數項的回歸模型，這在今后巖體初始地應力場的反演中應避免。

通過回歸模型及回歸系數的顯著性檢驗后得到的結果分別為

(15)

(16)

式中：σT≈0.016ZMPa。其中，Z=H-h，H為埋深，h為巖體恒溫帶離陸地表面的距離[16]。

回歸得到式(15)的相關系數為0.997，式(16)的相關系數為0.992，由此可見，兩類回歸模型子構造應力場與巖體初始地應力場之間線性相關的密切程度都較高。

3.5 基于殘差平方和取舍回歸模型常數項

依據式(4)、式(5)可得不含、含常數項的回歸模型應力實測值與計算值的大小殘差平方和分別為

(17)

(18)

(19)

tanα=(σ1-σx)/τxy

(20)

(21)

由式(7)可得含與不含常數項回歸模型應力的大小和方向偏差程度，分別為1 185.75、2 442.96。

由此可見，隧道采用含常數項回歸模型的應力大小和方向偏差程度較不含常數項回歸模型更小，故采用含常數項回歸模型的反演結果更好。因此，從反演結果出發，本文建議采用含常數項回歸模型對此局部區域的巖體初始地應力場進行反演分析。如果不采用圖2的反演流程，則有可能采用不含常數項回歸模型進行反演，而采用此類模型反演得到的反演結果較含常數項回歸模型更差。含與不含常數項回歸模型的反演結果對比見圖5～圖10。圖中，εi=實測值-計算值。

圖5 σx殘差平方和

圖6 σy殘差平方和

圖7 σz殘差平方和

圖8 τxy殘差平方和

圖9 含常數項回歸模型方位角的殘差平方和

圖10 不含常數項回歸模型方位角殘差平方和

由圖5～圖10可知，兩類回歸模型應力實測值與計算值在量值上存在一定的差異。含常數項回歸模型應力大小最大殘差為1.5 MPa，應力方向最大殘差為-7.8°；不含常數項回歸模型應力大小最大殘差為1.7 MPa，應力方向最大殘差為10.3°。即采用含常數項回歸模型反演得到應力實測值與計算值的大小和方向偏差程度更小，反演結果較不含常數項回歸模型更好，與式(22)、式(23)的結果一致。

綜上所述，在巖體初始地應力場反演過程中，采用含與不含常數項回歸模型得到的反演結果是不相同的。因此，僅采用其中一類回歸模型進行反演，很可能會導致反演結果較差的回歸模型被用于反演，從而導致反演得到的巖體初始地應力場不是最優的。故為提高巖體初始地應力場的反演結果，在反演過程中，以應力的大小和方向偏差程度為衡量標準，采用圖2的反演流程可客觀地獲得與實測原位地應力較為吻合的巖體初始地應力場。

4 結論

(3)在巖體初始地應力場反演中，若采用含常數項回歸模型，則應采用式(13)的正規方程組求解回歸系數；若采用不含常數項回歸模型，則應采用式(14)的正規方程組求解回歸系數。