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

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

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

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

1 地應(yīng)力場多元回歸反演分析方法

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

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

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

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

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

同理，也可以得出該觀測點(diǎn)的水平大主應(yīng)力和小主應(yīng)力的回歸方程

式中：下標(biāo)l、s分別表示水平大主應(yīng)力、水平小主應(yīng)力，其他符號含義同式（1）。

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

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

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

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

1.2 工況1有限元模擬計(jì)算

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

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

水平地應(yīng)力隨埋深的變化

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

σv、σhmax和σhmin都是主應(yīng)力并且水平方向的兩個主應(yīng)力大小相等。有限元計(jì)算不同觀測點(diǎn)3個方向的主應(yīng)力見表1。

表1 有限元計(jì)算工況1不同觀測點(diǎn)3個方向的主應(yīng)力

1.3 工況2有限元模擬計(jì)算

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

表2 有限元計(jì)算工況2不同觀測點(diǎn)3個方向的主應(yīng)力

1.4 工況3有限元模擬計(jì)算

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

表3 有限元計(jì)算工況3不同觀測點(diǎn)3個方向的主應(yīng)力

1.5 工況4有限元模擬

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

表4 有限元計(jì)算在工況4不同觀測點(diǎn)3個方向的主應(yīng)力

1.6 超定方程組求解

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

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

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

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

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

2 工程應(yīng)用

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

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

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

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

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

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

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

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

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

表6 地應(yīng)力反演預(yù)測值

3 結(jié)論

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

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

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

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