偏壓作用下各向異性深埋軟巖隧道的地應力反演

吳幫標，張浩東，于長一，劉愛民

（1.天津大學建筑工程學院，天津 300354；2.中交天津港灣工程研究院有限公司，天津 300222；3.中交第一航務工程局有限公司，天津 300461；4.港口巖土工程技術交通行業重點實驗室，天津300222；5.天津市港口巖土工程技術重點實驗室，天津 300222）

0 引言

初始應力場是隧道變形與失穩的主要影響因素，也是科學指導隧道開挖與支護的前提條件。在工程當中，可用于指導現場施工的初始應力場可視為忽略時間因素（地質年代）的相對穩定應力場，所以直接測量所得的地應力是具有指導意義的。在早期的工程當中，利用某種方法在地下進行實地測量，所得的實測資料，是唯一可供應用的初始地應力值。然而受成本和工期的限制，少數測點又難以滿足工程需要。考慮到在工程中運用少數測點的應力值進行反演，提供一個相對準確的初始應力場，便成為解決工程穩定性問題的一個有效方法[1]。

目前，地應力場的反演方法主要有位移反分析法[2]和應力反分析法[3]。位移反演分析法僅僅適用于地形地貌較簡單，地層分布較為單一的工況。而應力反演分析法又可分為兩類，第一類屬于定性分析的范疇，例如將初始地應力場劃分為自重應力場、構造應力場、溫度應力場和地磁應力場等。產生的方法有邊界荷載調整法、側壓系數法、應力函數法[4-6]等方法。第二類方法則是定量分析，目前采用的分析方法有：多元線性回歸分析法、神經網絡法[7-8]等。上述各種方法都或多或少存在著計算精度和收斂速度方面的問題，難以準確地指導設計和施工。

為此，本文結合云南某單線鐵路各向異性軟巖隧道，在定性分析的基礎上，采用griddata函數插值法對有限元模擬軟件所得的地應力場進行修正。本文主要包括以下幾個部分：首先介紹了工程概況以及研究背景；接著建立了有限元模型對地層的各項參數進行標定；隨之用模擬軟件模擬地層結構，初步得到監測斷面的地應力場；然后利用griddata函數對數值計算結果進行修正，最后將結果與實測值進行對比分析，驗證griddata函數法修正結果的準確性。

1 工程概況

云南某鐵路隧道進口里程DK405+615，出口里程DK415+124，全長9 509 m，施工隧道最大埋深711 m。

隧道整體為構造剝蝕、侵蝕中低山地貌。為了更好地對軟巖隧道的地應力分布進行研究，本文選取的研究對象為隧道里程DK414+647處監測斷面，此斷面處圍巖主要為侏羅紀和平鄉組（J2h）、小紅橋組（J2x）泥巖、砂巖夾頁巖、礫巖等，主要巖性的巖石堅硬程度劃分為軟巖（5 MPa≤Rc≤15 MPa），受地質構造影響較重或嚴重，節理較發育或發育，具有典型的軟巖特征。監測斷面圍巖部分呈塊（石）碎（石）狀鑲嵌結構的部分，接近隧道輪廓線部分巖體表現為較破碎至破碎，呈角礫碎石狀松散結構。

在隧道開挖后，局部發生軟巖大變形，易發生坍塌、掉塊、側向擠入，造成很大的安全隱患。為了預防隧道施工中產生的高風險，在開挖后迅速于監測斷面布置壓力盒，用于監測圍巖壓力的變化情況，充分對該類隧道進行圍巖壓力和支護的應力量測，此時測得的圍巖壓力就等于隧道臨空區最邊緣的應力，其變化曲線可以理解為斷面開挖后應力重分布的變化過程，此應力值將作為地應力反演的一個重要依據。因此，選擇圍巖壓力時程曲線20 d之后壓力值趨于穩定的數值作為圍巖最終壓力值，圖1為DK414+647圍巖壓力值橫斷面分布圖。

圖1 圍巖壓力值橫斷面分布圖（kPa）Fig.1 Cross-section distribution map of surrounding rock pressure value(kPa)

2 有限元模型及參數標定

本文采用有限元軟件進行數值模擬分析，根據工程前期所得的測點圍巖壓力，結合數值模擬實驗驗證，分析地應力分布的影響因素，對局部斷面的地應力場進行反演。

2.1 幾何模型

根據該隧道的地質報告，隧道走向與最大水平主應力的走向之間的夾角始終恒定，即沿著隧道走向，任意斷面的地應力情況均相同，同時二維模型假設荷載沿著隧道軸向均勻分布，故將三維的地應力問題簡化為二維模型進行計算。為消除邊界影響[9]，二維模型選取隧道橫斷面內以隧洞為中心半徑大于3倍洞徑的方形區域，分析范圍為100 m×100 m；開挖隧洞為橢圓形，最高處9.38 m，最寬處7.30 m。計算模型總共劃分單元數為6 878個，結點數為6 997個。

2.2 本構參數

本文主要研究初始地應力的分布，不涉及開挖變形等的討論，因此本構模型選取摩爾庫倫模型；考慮到開挖后的應力松弛，根據極限平衡自穩拱原理[10]，可計算出隧道開挖后自穩拱的最大高度作為拉應力區的半徑最大值，故近似的將在開挖分析步中調整隧洞輪廓線向外延伸2 m的范圍內的圍巖參數。受開挖擾動影響，隧道周圍物理參數開挖前后有差異，將沿隧道徑向取出的土樣分段進行室內實驗，將土樣2 m之內和之外的參數進行分別記錄，最終確定的開挖前后圍巖各項物理力學參數如表1所示。

表1 圍巖物理力學參數Table 1 Physical and mechanical parameters of surrounding rock

2.3 地層參數標定

依據現場地質勘查資料，結合王福和[11]在實際偏壓隧道工程中的模擬結果進行分析，該監測斷面上部受到偏壓的影響，且其圍巖類別多為泥質頁巖，故洞周圍巖的應力值受到上部偏壓和各向異性的影響，因此本文通過數值模擬實驗的邊界荷載調整法來還原地應力場，在模型上部的1/4處、1/2處、3/4處分別施加偏壓荷載，觀察其應力集中區域；各向異性設置為正交各向異性，偏轉角度分為30°、35°、40°、45°、50°、55°、60°7種情況進行實驗（圖2），得出最符合實際情況的邊界荷載及圍巖材料方向。

圖2 實驗方案Fig.2 Experimental program

2.4 數值模擬工況

本次模擬的目的是還原隧道開挖后監測斷面的應力值，因此數值模擬的第一步是還原開挖之前的真實應力場，并且進行地應力平衡，第二步進行隧洞開挖，計算完成后導出模型中對應監測點的應力值，并與實測值進行對比。

2.5 反演邊界條件

根據模擬實驗結果，計劃在進行反演時對模型偏右3/4上邊界施加偏壓荷載，對整個模型上邊界施加上部土體壓強，整個模型施加自身重力，同時設置各向異性材料方向為35°（圖3），模型左右邊界限制X方向上的位移，下部邊界限制X、Y兩個方向上的位移。

圖3 邊界條件Fig.3 Boundary conditions

3 地應力場反演結果

3.1 地應力場初步反演結果

在建立有限元模型的基礎上，通過模擬實驗得到監測斷面的偏壓與各向異性情況，并計算得到監測點處最接近實測值的應力云圖，如圖4所示。

圖4 開挖后監測斷面應力云圖Fig.4 Monitoring section stress cloud map after excavation

從圖中可以看出，通過數值模擬計算出的應力值與實測值較為接近，同樣為左拱腰和右拱腳出現應力集中，應力水平的傾向性較為明顯，符合現場實測資料。由于數值模擬中邊界條件與實際工程中的情況并不相同，且工程中的巖性更為復雜，而使在右拱腰與左拱腳出處應力值有較大的出入，但是應力集中與應力較弱區域的分布情況與實際值是相似的。

盡管模擬出的結果在分布和量值上與實測值相接近，但洞周圍巖的壓力值與實測值之間仍存在著不小的出入，此種誤差的存在，使模擬結果難以作為指導支護施作的主要依據。因此，采用griddata函數作為修正函數結合現場實測值對模擬結果進行修正，使反演值更加接近實際地應力值。

3.2 采用griddata函數法修正地應力場

為了保證地應力的反演結果盡可能的與實測值接近，故在數值模擬的基礎上，結合有限的實測點應力值，運用Matlab中的griddata函數進行修正。

griddata函數可以將位于同一空間坐標系下的散點插值為規則格網，提供了以下4種方法：線性內插法、三次多項式內插法、最鄰近點內插法和Matlab自帶的一種圓滑插值法，可方便地實現結合臨近離散點分布特征的光滑曲面擬合。其程序實現偽代碼如下：

[xi，yi]=meshgrid（min（x）∶dx∶max（x），min（y）∶dy∶max（y）），%，生成規格化格網；

[xi，yi，zi]=griddata（x，y，z，xi，yi，′method′），%，應力值平面擬合；

zi=griddata（x，y，z，xi，yi，′method′），%，求待擬合點的應力值；

surf（xi，yi，zi），%，應力云圖。

從數值模擬的結果分析，地應力場的應力集中區域與實測相對符合，但是隧道斷面邊緣的應力值在部分區域則有些出入，而此處反演有誤差的部分恰好具有實測點的數值，因此，將應力松弛區域以外的模擬應力值與區域內的實測值運用griddata函數再次進行插值計算，可得到與實測值更為接近的地應力場。

3.3 采用griddata函數法修正結果

結合數值模擬計算，采用griddata函數圓滑插值法對斷面地應力場進行修正的結果見表2。

表2 反演應力值結果對比Table 2 Comparison of inversion stress values

從反演結果中可以看到，在隧道左拱腰處與右側墻腳處同樣出現應力集中，與實測和數值模擬計算兩種方法相比，分布區域均類似，偏壓荷載與各向異性影響明顯，數值模擬方法中所得的結論再次得到驗證。同時，反演結果的洞周測點處圍巖應力數值與實測值完全相同，顯然在地應力反演的計算精度上，griddata函數有明顯的優勢。

4 反演應力值對比分析

將用griddata函數法修正后的結果和僅用數值模擬方法計算所得的圍巖壓力與實測值進行比較，均發現在左拱腰與右拱腳處有明顯的應力集中現象，在右拱腰與左拱腳處的應力相對較小，應力集中與應力較弱的區域分布均相同。

從表2的對比結果可以看出，數值模擬軟件的模擬誤差最大為229%，誤差較大，沒有參考價值；右拱腰處的誤差較大，但是應力集中區域的分布類似，可以認為是材料參數與邊界條件引起的誤差；其余點的相對誤差主體在20%以下，超出30%的點為1個。而在地應力測量階段，測量方法本身便會導致測量結果的誤差，據國內外的統計資料所示，初始地應力值的測量結果誤差允許達到25%～30%[12]，故而認為數值模擬方法所得的地應力值與實測值基本符合。

由于插值計算法多為數學模型，難以確定其誤差，故其計算精度可用內外符合精度表示，內符合精度為擬合點與擬合曲線的差值，外符合精度為核檢點與擬合值的插值。griddata函數擬合曲面通過各擬合點，故擬合點的內符合精度為0，據相關學者研究[13]，griddata函數插值點的殘差與外符合精度相較于多面函數有很大的提高，與最小二乘法相比也具有相當的精度。可以認為griddata函數插值法高度的還原了工程實際中的地應力值。本文分析的地形中圍巖情況較為復雜，有明顯的軟巖性質和各向異性，采用griddata函數修正法能夠很好地模擬出隧道中應力集中的區域，表現出其地質構造對地應力分布的影響。

5 結語

針對各向異性的深埋軟巖隧道地應力反演問題，本文首次提出采用griddata函數對模擬軟件得到的地應力進行修正，并結合有限元數值模擬試驗對該問題進行了深入研究，結果如下：

1）本文利用數值模擬軟件對云南某軟巖隧道某監測斷面進行數值模擬，盡管模型邊界條件與巖石物理力學參數并不能與實際完全一致，從而使部分結果存在誤差，但是通過試算法仍然模擬出了與開挖斷面的圍巖壓力值相近的結果，并且能供實際工程參考。

2）griddata函數修正的計算結果具有高精度。本文采用基于griddata函數的函數修正法，得到了其程序實現方法并繪出了計算結果云圖。通過對比分析，griddata函數修正計算結果應力集中區域與數值模擬結果相似，應力值與實測值高度相似，證明griddata函數具有相當高的計算精度，結合現場實測值，能夠為現場施工做出有效的指導。