多元線性回歸的地應力場反演分析及應用

王志云，于 申，李守巨，楊朔明，樊可為

（1.大連海洋大學海洋與土木工程學院，遼寧 大連 116023；2.大連理工大學 工業(yè)裝備結(jié)構(gòu)分析國家重點實驗室，遼寧 大連 116024）

巖體內(nèi)的初始應力場直接影響隧道和硐室圍巖的穩(wěn)定性；因此，合理確定地下巖體的應力場對隧道圍巖穩(wěn)定性評估和支護襯砌結(jié)構(gòu)參數(shù)選取至關(guān)重要。在漫長的地質(zhì)年代地殼構(gòu)造運動作用下，巖層會產(chǎn)生褶皺、彎曲、斷裂和錯動，使得巖層種除了自重產(chǎn)生的應力場之外，還會有2個水平方向的構(gòu)造應力和1個垂直方向的構(gòu)造應力。構(gòu)造應力的大小和方向也可能會隨著埋深的增加而發(fā)生變化；作為一種最簡單的假設，認為構(gòu)造應力的大小不隨深度的增加而變化，是一個常數(shù)。該假設可能使地應力的擬合精度下降，但也使得地應力場反演算法的復雜程度大大降低，得到的地應力場分布結(jié)果也基本上滿足工程建設精度要求。郭喜峰等[1]以清遠抽水性能電站為例，研究了地應力場的多元回歸分析反演方法，建立了基于最小二乘法原理的地應力場參數(shù)計算的法方程的系數(shù)矩陣；周春華等[2]以鄂西某隧道工程為例，對地應力場反演方法及圍巖穩(wěn)定性進行了分析；謝洪強等[3]對錦屏水電站壩區(qū)初始地應力場進行了回歸反演分析；李永松等[4]研究了深圳抽水性能電站的應力場反演方法，提出了4種荷載工況下的加載方法；何少云等[5]基于地應力實測結(jié)果，采用多元線性回歸方法和三維有限差分法，研究了某水電站地下廠房區(qū)域的初始地應力場反演問題；徐中衛(wèi)等[6]采用多元回歸方法、神經(jīng)網(wǎng)絡和遺傳算法3種方法，研究了蟠龍抽水蓄能電站地下工程區(qū)域的初始地應力場反演問題；李守巨等[7～8]采用優(yōu)化算法研究了地下廠房巖體應力場的反演問題。

本文提出一種基于求解超定線性方程組的地應力場反演方法，建立地應力場多元線性回歸的數(shù)學模型，根據(jù)現(xiàn)場觀測地應力場數(shù)據(jù)，反演確定線性回歸方程的系數(shù)，采用有限元模型模擬得到某地下廠房區(qū)域的整體應力場分布。

1 地應力場多元回歸反演分析方法

為確定地下巖體內(nèi)的初始應力場和進行開挖后的圍巖穩(wěn)定性，采用水壓致裂法或應力解除法進行地應力現(xiàn)場測試，一般需要幾個鉆孔，以便對測試結(jié)果相互驗證。某地下廠房處于微風化花崗巖巖體內(nèi)，在DH008鉆孔內(nèi)布置5個觀測點，每個觀測點觀測3個主應力數(shù)值，共15個觀測數(shù)據(jù)。假設巖體內(nèi)任意觀測點的3個主應力分別與巖體的自重（工況1）、大主應力方向（有限元模型x坐標）的構(gòu)造應力（工況2）、小主應力方向（有限元模型y坐標）的構(gòu)造應力（工況3）和垂直方向（有限元模型z坐標）的構(gòu)造應力（工況4）有關(guān)。根據(jù)現(xiàn)場觀測，大主應力方向為N30E；因此，有限元計算模型取x坐標方向為N30E，y坐標方向為N60W，z坐標方向為垂直向上。這樣，在有限元計算模型邊界上，x、y、z方向均為主應力方向，與坐標軸平行的有限元模型邊界上沒有剪應力。見圖1。

圖1 4種典型荷載工況簡化加載模型

1.1 地應力場多元線性回歸的數(shù)學模型

基于多元回歸分析的基本原理，以第一個觀測點的垂直主應力為例，將垂直主應力的回歸計算值SRv作為函數(shù)，把有限元在4種荷載工況下分別計算求得的該觀測點垂直主應力計算值Sv1、Sv2、Sv3、Sv4作為自變量，則該觀測點地應力的垂直主應力線性回歸方程

式中：Lz為巖石自重的多元回歸系數(shù)；Lx為水平大主應力方向構(gòu)造應力的回歸系數(shù)；Ly為水平小主應力方向構(gòu)造應力的回歸系數(shù)；Lt為垂直構(gòu)造應力的回歸系數(shù)。

同理，也可以得出該觀測點的水平大主應力和小主應力的回歸方程

式中：下標l、s分別表示水平大主應力、水平小主應力，其他符號含義同式（1）。

巖體地應力場反演的目的在于使經(jīng)過回歸分析得到的所有觀測點3個主應力與現(xiàn)場觀測值的殘差平方和最小，定義目標函數(shù)J

式中：i表示第i個觀測點，i=1，2，3，4，5；Sm為主應力現(xiàn)場觀測值。

為求解該極小值，未知數(shù)為回歸方程的回歸系數(shù)，令目標函數(shù)J對于回歸系數(shù)的導數(shù)等于0，根據(jù)復合函數(shù)求導規(guī)則，推導出

一個觀測孔有5個觀測點，每個觀測點3個主應力觀測值，一共15個方程，4個未知數(shù)，方程的個數(shù)大于未知數(shù)的個數(shù)，屬于超定方程組。MATLAB軟件為解決超定方程組問題提供了數(shù)值計算精度的保證。根據(jù)地勘報告和地應力現(xiàn)場觀測結(jié)果，確定巖體基本力學參數(shù)：密度2 600 kg/m3、彈性模量9.2 GPa、泊松比0.27、抗壓強度65.3 MPa。初步假定3個方向上的構(gòu)造應力作為有限元數(shù)值模擬所需假定的邊界條件：垂直方向構(gòu)造應力5.00 MPa、水平大主應力方向構(gòu)造應力6.0 MPa、水平小主應力方向構(gòu)造應力3.0 MPa。

1.2 工況1有限元模擬計算

采用巖體實測密度、彈性模量和泊松比，計算在自重作用下所產(chǎn)生的自重應力場。其中，有限元模型底部施加垂直位移為零約束，dz=0；在與x軸垂直的兩個平面施加x方向位移為零約束，dx=0；在y軸垂直的兩個平面施加y方向位移為零約束，dy=0。在自重應力作用下，巖體內(nèi)的垂直地應力隨埋深的變化

式中：ρ為巖石的密度；g為重力加速度；h為埋深。

水平地應力隨埋深的變化

式中：μ為巖石的泊松比。

σv、σhmax和σhmin都是主應力并且水平方向的兩個主應力大小相等。有限元計算不同觀測點3個方向的主應力見表1。

表1 有限元計算工況1不同觀測點3個方向的主應力

1.3 工況2有限元模擬計算

取巖體的密度為0，在與x方向垂直的兩個邊界面施加均勻分布6 MPa的構(gòu)造應力；而與y方向垂直的兩個邊界面施加位移為零約束，dy=0；在有限元模型的底部施加垂直位移為零約束，dz=0。不同觀測點3個方向的主應力見表2。

表2 有限元計算工況2不同觀測點3個方向的主應力

1.4 工況3有限元模擬計算

取巖體密度為0，在與y方向垂直的兩個邊界面施加均勻分布3 MPa的構(gòu)造應力；而與x方向垂直的兩個邊界面施加位移為零約束，dx=0，在有限元模型的底部施加垂直位移為零約束，dz=0。不同觀測點3個方向的主應力見表3。

表3 有限元計算工況3不同觀測點3個方向的主應力

1.5 工況4有限元模擬

取巖體密度為0，在模型頂部施加5 MPa的局部壓力；在與y方向垂直的兩個邊界面施加位移為零約束，dy=0；而與x方向垂直的兩個邊界面施加位移為零約束，dx=0；在有限元模型的底部施加垂直位移為零約束，dz=0。不同觀測點3個方向的主應力見表4。

表4 有限元計算在工況4不同觀測點3個方向的主應力

1.6 超定方程組求解

根據(jù)地應力場多元線性回歸的數(shù)學模型和4種工況下有限元數(shù)值模擬結(jié)果，確定地應力場多元線性回歸超定方程組系數(shù)矩陣。

式中：A為超定線性方程組的系數(shù)矩陣（15×4）；L為回歸方程的系數(shù)矢量（4×1）；Sm為現(xiàn)場觀測數(shù)據(jù)矢量（15×1）。

超定方程組的系數(shù)矩陣是根據(jù)式（5）～（7）及表1-表4，得到MATLAB格式的系數(shù)矩陣，排列順序依次為觀測點1垂直主應力在4個工況有限元計算值、大主應力在4個工況有限元計算值、小主應力在4個工況有限元計算值；觀測點2、觀測點3、觀測點4、觀測點5的3個主應力4個工況有限元計算值。

另外一種計算回歸方程系數(shù)的方法是將式（10）兩端同時左乘系數(shù)矩陣A的轉(zhuǎn)置矩陣AT。

經(jīng)過變換后，式（11）的系數(shù)矩陣為4×4方陣，直接求解該線性方程組就可以得到回歸方程的系數(shù)。

2 工程應用

在鉆孔編號DH008不同高程布置5個地應力測試點，每個觀測點觀測到的3個主應力，該組觀測數(shù)據(jù)是地應力場反演的依據(jù)。見表5。

表5 DH008地應力場現(xiàn)場觀測數(shù)據(jù)

地應力場多元線性回歸超定方程組的右端項（列向量）是根據(jù)表5按照順序與系數(shù)矩陣排序相對應：觀測點1的垂直主應力、大主應力和小主應力；觀測點2、觀測點3、觀測點4、觀測點5的垂直主應力、大主應力和小主應力，即Sm列向量。于是，超定方程組的解的MATLAB表達式為

采用MATLAB軟件，求解地應力場線性回歸分析超定方程組，得到4個線性回歸系數(shù)。

由于地應力現(xiàn)場觀測誤差的存在，直接求解超定方程組得到的解有時偏離工程實際。考慮到在地應力現(xiàn)場測試中，垂直主應力的精度最高，采取如下策略求解：根據(jù)垂直主應力的觀測數(shù)據(jù)，先確定與垂直地應力相關(guān)的兩個回歸系數(shù)（Lz和Lt）；再根據(jù)水平大主應力和小主應力的觀測數(shù)據(jù)，確定另外兩個回歸系數(shù)（Lx和Ly）。從表2和表3可以看出，工況2和工況3對垂直主應力沒有任何影響；因此，將式（5）簡化為

根據(jù)5個觀測點的垂直主應力數(shù)據(jù)，可以列出5個線性方程，求解該超定方程組，求出兩個回歸系數(shù)Lz=0.982 7和Lt=0.918 7后，式（6）和式（7）改寫為

同樣根據(jù)5個觀測點的水平大主應力和小主應力數(shù)據(jù)，求解10個超定方程組，得到另外兩個回歸系數(shù)Lx=0.859 4和Ly=0.604 4。

為驗證所得應力場參數(shù)的精度，根據(jù)表5中垂直地應力現(xiàn)場觀測數(shù)據(jù)，線性擬合得到垂直地應力隨埋深的變化

對比式（16）和表5中反演所得到的垂直構(gòu)造應力可以發(fā)現(xiàn)，反演分析和擬合分析所得到的垂直方向的構(gòu)造應力分別為4.59、4.57 MPa，說明反演方法所得到的地應力場參數(shù)是十分準確的。回歸反演確定的構(gòu)造應力：垂直方向應力為4.593 5 MPa、水平大主應力為5.156 4 MPa、水平小主應力為1.813 2 MPa，將回歸分析反演得到的回歸系數(shù)和有限元計算得到的表1-表4帶入到式（1）～式（3），得到5個觀測點地應力反演預測值。見表6。

表6 地應力反演預測值

3 結(jié)論

1）采用多元線性回歸分析方法，建立了地應力場回歸分析的數(shù)學模型，基于有限元數(shù)值模擬方法和現(xiàn)場觀測的地應力數(shù)據(jù)反演確定了地應力場參數(shù)。該方法為準確估計該區(qū)域的地應力場三維分布特性提供了依據(jù)，為后期地下廠房的圍巖穩(wěn)定性評估奠定了基礎。

2）計算結(jié)果表明，反演分析和擬合分析所得到的垂直方向的構(gòu)造應力分別為4.59、4.57 MPa，說明反演方法所得到的地應力場參數(shù)是十分準確的。相對比而言，垂直地應力場構(gòu)造應力擬合精度最高，最大相對誤差＜2%；水平小主應力方向的構(gòu)造應力最大誤差＜7%，水平大主應力方向的構(gòu)造應力擬合精度有一定的誤差。

3）本文僅給出了基于單個鉆孔5個觀測點地應力觀測數(shù)據(jù)反演得到的地應力場數(shù)據(jù)；多個觀測孔現(xiàn)場觀測數(shù)據(jù)時，反演方法與單孔的反演地應力場方法類似，只不過隨著數(shù)據(jù)增多，超定方程組的行數(shù)更多，線性回歸方程的系數(shù)個數(shù)不變，觀測數(shù)據(jù)有時可能會出現(xiàn)冗余，甚至有些觀測數(shù)據(jù)相互不協(xié)調(diào)或者相互矛盾的情況。

4）由于巖體地質(zhì)構(gòu)造運動的復雜性和不確定性，在某個方向上的構(gòu)造應力可能不是常數(shù)，即有限元計算模型邊界上的側(cè)壓力（構(gòu)造應力）分布可能不是矩形，而是一個梯形分布，這時就會增加線性回歸方程系數(shù)的個數(shù)，但會使得反演預測精度進一步提高。