長大深埋高鐵隧道三維地應力場反演方法及應用
——以銀河山隧道為例

李 鵬 袁 維 張光明 王欽剛 李國圣

(1. 中國鐵路設計集團有限公司,天津 300308; 2. 城市軌道交通數字化建設與測評技術國家工程研究中心,天津 300308;3. 石家莊鐵道大學土木工程學院,石家莊 050043)

引言

在長大深埋隧道的設計和施工過程中,準確評估隧道區域巖體的初始地應力場對隧道設計和施工至關重要。 根據初始地應力場分布特征設計合理的開挖方案以及支護方式,可有效防止巖爆、大變形、塌方等事故的發生[1-3]。 在工程實踐中,難以直接測量隧道開挖影響區域內各點全部的地應力,只能通過水壓致裂等方法測出若干關鍵測點的地應力,然后采用相應的數學方法按照“點→線→體”的思路對隧道影響區域內的地應力進行預測,以獲得初始地應力場的大致分布規律,為工程設計和施工提供指導。 這種獲取地應力分布規律的方法被稱為地應力反演[4-6]。

地應力反演方法有很多種,目前應用較廣泛的方法有應力回歸分析法[7]、邊界荷載調整法[8]、應力試算法[9]等。 已有學者進行相關研究,張敏等采用應力試算法確定模型加載的邊界條件,并結合FLAC3D 分析得到桑珠嶺隧址區范圍的初始地應力場[10];張晨等對BP 神經網絡進行優化,結合數值模擬和現場數據,反演計算實際工程的地應力狀態[11];王志云等采用線性方程組和線性回歸相結合的數學方法,通過數值模擬和現場數據驗證該方法的可靠性[12];劉春康等建立某隧道三維數值模型,并基于多元線性回歸法反演獲得該隧道軸線處的初始地應力場[13];徐安等基于實測數據反推得到浩吉鐵路九嶺山隧道和桐木隧道初始應力場應力邊界條件,并通過數值模擬方法獲得初始應力場分布特征[14]。

以某高鐵銀河山隧道為研究對象,建立其三維數值分析模型,設定多組非線性荷載邊界條件模式樣本數據庫,通過數值模擬方法獲得地應力實際觀測點的地應力值,將邊界條件數據和地應力計算值分別作為輸入和輸出樣本,通過BP 神經網絡進行訓練,并結合遺傳算法尋優功能得到地應力實測數據對應的最優應力邊界條件,將該邊界條件輸入到銀河山隧道的三維數值分析模型中進行正算,獲得該隧道模型的三維初始地應力場。 此外,采用多元線性回歸分析方法對該隧道地應力進行反演,并與上述反演結果進行對比驗證。

1 基本原理

1.1 神經網絡基本簡介

神經網絡為模擬人類大腦神經系統的一種信息處理系統,該網絡由許多神經元組成,用于大量數據信息的學習和處理,通過對輸入、輸出樣本的自學習,并進行網絡結構和連接權值訓練修正,從而建立復雜的非線性預測模型,達到對新輸入數據進行預測的目的[15-16]。

BP 神經網絡是一種廣泛應用于工程領域的模型,其算法是將一組樣本的輸入、輸出變成非線性優化問題,通過調整BP 網絡中的連接權值、網絡規模使實際輸出和預期輸出之間誤差的平方最小,從而實現任意精度逼近任何非線性函數。 按照層次構造,BP 網絡包含一個輸入層、一個輸出層和一個或多個隱含層,神經元只和與該層緊鄰的上、下層的各神經元連接,BP 網絡結構結構見圖1。

1.2 遺傳算法基本簡介

遺傳算法效仿生物的進化、遺傳過程而提出的隨機搜索方法,即從某一初始值群體出發,經過生存競爭→優勝劣汰→適者生存的過程,基于復制、雜交、變異等操作手段,不斷迭代演化逼近最優解,直到實現全局最優。遺傳算法包含編碼、初始群體的生成、適應度評估檢測、復制、雜交、變異共6 個基本步驟,流程見圖2[17-18]。

圖2 遺傳算法計算流程

2 工程概況

某高鐵銀河山隧道進口位于河北省保定市阜平縣,出口位于山西省忻州市,穿越太行山脈剝蝕中山區,全長12 549.74 m。 為了探明該區域地應力狀態,在隧道主軸線上鉆孔,鉆孔深度為290 m,靜水位9.6 m,鉆孔直徑78 mm。 在孔內116~220 m 深度內選定了4 個測試段,深度分別為116.50 ~ 117.10 m、130.20~ 130.80 m、148.00 ~ 148.60 m 和161.70 ~162.30 m。 根據水力壓裂的基本原理,通過計算得到最大水平主應力值(SH)和最小水平主應力值(Sh);此外,垂直應力(Sv)采用上覆巖體的容重進行估算。 因此按巖石密度(2.65 g/cm3)估算出各測段處的垂直應力值,水壓致裂法地應力測試結果見表1。

表1 水壓致裂法地應力測試結果

3 三維地應力場非線性反演

3.1 基本流程

基于“反演+正算”原理的地應力場非線性反演方法計算步驟如下。

(1)根據地形平面、鉆孔剖面等地質資料建立該隧道三維地質力學數值模擬模型。

(2)將應力邊界條件作為輸入參數,設置不同的應力邊界條件參數,采用有限差分法對不同應力邊界條件下的數值模型開展巖體地應力計算分析。

(3)提取上述數值模擬結果中地應力測點處的計算值,將邊界條件參數作為輸入樣本、各個測點的地應力計算值作為輸出樣本,組建訓練樣本數據庫。

(4)利用BP 神經網絡的自學習功能,對上述樣本進行訓練、檢驗,構建非線性預測模型。

(5)確定隧道模型實際應力邊界參數的搜索范圍,將輸入參數在搜索范圍內進行離散化,將每組參數輸入到訓練好的非線性預測模型中,從而得到對應的地應力預測值。

(6)將地應力實測值作為期望,通過最佳逼近目標函數將地應力計算值與期望值比較,采用遺傳算法尋找使目標函數取得最小值的參數組,即參變量全局尋優。

(7)將第(6)步得到的最優參數組輸入到隧道數值分析模型中進行“正”算,從而獲得地應力場。

3.2 數值模型

采用有限差分軟件Flac3D 建立三維數值分析模型。 模型以地應力測孔為中心,沿隧道縱向朝兩邊延伸5 000 m、垂直隧道兩邊各延伸3 000 m,豎直向下模擬至高程89 m。 銀河山隧道地表地質模型見圖3(單元數為52 109,節點數為10 026)。

圖3 銀河山隧道地表地質模型

銀河山隧道區域地質條件復雜,巖體分區、分層較多,地質歷史上構造運動較強烈,地表風化、剝蝕等外在地質營力非常明顯,巖體應力條件極其復雜。 因此,地應力反演采用彈塑性本構模型進行計算,模擬時相關參數見表2。

表2 巖體力學參數

3.3 地應力反演

數值模型的邊界水平投影是規則的長方形,2 組平面分別沿東西向和南北向,恰好與該區域的最大、最小主應力方向平行,模型邊界上的剪應力可忽略不計。為簡便起見,只在模型的4 個側面邊界上施加法向應力荷載。

設定24 組非線性荷載邊界條件模式樣本數據(見表3),并利用FLAC3D 軟件對每組參數樣本進行數值分析。 將其中20 組計算結果用于神經網絡訓練,剩下的4 組則用于訓練效果檢驗。

表3 邊界條件參數樣本數據

樣本計算結束后,從上至下選擇4 個測點的最大和最小的水平主應力計算值作為學習訓練的輸出項,而將應力邊界參數作為對應輸入項。 輸入值的格式為aG、ax、bx、cx、dx、ay、by、cy、dy,輸入值見表3;輸出值的格式為測點1 的SH、測點1 的Sh、測點2 的SH、測點2 的Sh、測點3 的SH、測點3 的Sh、測點4 的SH、測點4 的Sh,輸出值見表4。

表4 地應力計算值樣本數據MPa

把樣本數據分為訓練數據和檢驗數據,其中,訓練數據為1~ 20 號,檢驗數據為21~ 24 號。 在正式學習訓練前,選定合適的允許誤差、學習速率、沖量系數和最大迭代次數等必須參數。 選定后,開展神經網絡自學習訓練,并對訓練效果進行同步檢驗。 如果效果不理想,則需要對模型進行再次補充訓練,直至檢驗效果達到要求。 4 個檢驗樣本預測值和計算值對比見表5。

表5 神經網絡學習訓練檢驗

根據地應力實測結果,估算出應力邊界條件的取值范圍,將其作為遺傳算法的參數值搜索范圍。 應用遺傳算法優異的全局尋優能力,在參數值域范圍內搜索最逼近地應力實測值的應力邊界條件參數組,其結果為aG= 0.813 731,ax= 0.005 402,bx= 0.113 426,cx= 0.056 345,dx= 0.362 836,ay= 0.008 859,by=0.073 725,cy=0.171 401,dy=0.382 576。 將以上參數搜索結果代入表3,并將應力邊界計算值施加到數值模型中再次進行數值模擬“正算”,從而獲取銀河山隧道地應力場的分布結果。 隧道縱剖面的水平最小主應力和最大主應力見圖4、圖5。

圖4 水平最小主應力(單位:Pa)

圖5 水平最大主應力(單位:Pa)

最后,將地應力實測值與反演結果進行對比,結果見表6。

表6 地應力實測值與反演結果的對比

由表6 可知,非線性反演方法得到地應力值與實測值比較接近,最大水平主應力、最小水平主應力的最大相對誤差分別是29%和26%,平均偏差分別為20.01%和19.80%。 根據相關經驗,地應力現場測試的誤差范圍一般為20%~ 30%[19-20],由此可知,本次反演結果可靠。

4 對比分析

4.1 荷載工況簡化模型

地應力場由自重應力場和構造應力場疊加而成,多元回歸分析方法是分別構建自重和構造應力場,然后采用最小二乘法進行回歸分析并完成回歸顯著性檢驗,之后建立合適的線性應力回歸函數,同時將二維應力場擴展至三維應力場。 該方法從局部到整體,能夠反映大范圍內的地應力場特征,方法簡單且應用廣泛。

假設巖體中任意觀測點的最大和最小水平主應力分別與巖體自重(情況1)、小主應力方向的構造應力(情況2)和大主應力方向的構造應力(情況3)有關,本研究中水平初始邊界最大、最小構造應力荷載值分別為8 MPa、6 MPa,見圖6[12]。 因此,各個觀測點的地應力多元線性回歸方程可以表示為

圖6 典型荷載工況簡化加載模型

式中,^σ為觀測點的水平最大主應力或水平最小主應力;σv、σu、σw分別為觀測點在工況1、工況2、工況3 下的最大或最小水平主應力計算值;d為誤差;a、b、c為待回歸系數。

4.2 數值模擬計算

(1)工況1

利用實測的巖體密度、彈性模量和泊松比,計算重力產生的自重應力場。 網格模型中,固定住該模型z向的位移,數值模擬時在模型的底部平面施加法向位移約束(即z方向位移恒為0)。 同時,將垂直于x軸的2 個平面的x軸方向位移設置為0;將垂直于y軸的2 個平面的y方向位移設置為0。

(2)工況2

巖體的力學參數和網格模型同情況1,在工況2 時取重力加速度為0。 固定住該模型z向的位移,數值模擬時把約束加在模型底面;在垂直于x軸的2 個平面施加沿x軸方向位的6 MPa 均勻分布的力;固定垂直于y軸的2 個平面的y方向位移,即令其在平行y軸方向不能移動。

(3)工況3

巖體的力學參數和網格模型同情況1,重力加速度取為0。 邊界約束為在垂直于y方向的2 個邊界面上施加8 MPa 均勻分布的構造應力;但是,垂直于x方向的2 個邊界曲面施加0 位移約束,dx=0,垂直位移約束施加在模型底部,dz=0。

采用彈性本構模型分別對3 種工況進行數值模擬計算,由此得到各種工況4 個觀測點的最大、最小計算水平主應力見表7。

表7 不同工況最大、最小水平主應力 MPa

4.3 回歸模型參數求解

根據上述計算結果,確定了地應力場多元線性回歸超定方程組的系數矩陣。

式中,A為超靜定線性方程組的系數矩陣(8×4);L為回歸方程的系數矢量(4×1);Sm為現場實測數據矢量(8×1)。 矩陣A及矢量Sm的參數見表8。

表8 矩陣A 及矢量Sm 相關參數

通過多元回歸分析,獲得的回歸系數分別為a=10.925 0,b=10.372 0,c=-6.611 3,d=-1.688 5,復相關系數為R=0.867 3;表明回歸效果顯著,回歸計算模式合理。 將3 種工況的觀測點計算值代入式(1),得到Sm的擬合計算值為[4.124 8,4.119 1,4.361 8,4.858 3,5.486 0,5.427 1,10.727 4,6.395 4]。SH、Sh反演值和實際值對比見圖7、圖8。

圖7 最大水平主應力實測值與非線性反演值、回歸分析值對照

圖8 最小水平主應力實測值與非線性反演值、回歸分析值對照

由圖7 可知,針對最大水平主應力,回歸分析的結果與實測值的平均偏差為28.3%,非線性反演分析的結果與實測值的平均偏差為20.01%。

由圖8 可知,針對最小水平主應力而言,回歸分析的結果與實測的平均偏差為21.03%,非線性反演分析的結果與實測值的平均偏差為19.8%。 由此可見,非線性反演方法計算結果優于多元回歸分析方法。

5 結論

基于“反演+正算”原理,建立地應力場反演三維數值分析模型,設定多組非線性荷載邊界條件模式,通過數值模擬方法獲得測點地應力計算值,將邊界條件數據和地應力計算值分別作為輸入和輸出樣本,再通過BP 神經網絡將樣本進行訓練,并結合遺傳算法尋優功能得到地應力實測數據對應的最優應力邊界條件,之后,將該邊界條件輸入到地應力場反演三維數值分析模型中進行計算,獲得該模型三維地應力場。 采用上述方法對某高鐵銀河山隧道地應力場進行反演,同時與多元回歸分析方法的反演結果進行對比驗證。結果表明,地應力場具有高度的非線性特征,與多元線性回歸分析方法相比,非線性反演方法得到的結果與實測值的平均偏差分別為20.01%(最大水平主應力)、19.80%(最小水平主應力),明顯小于多元線性回歸分析方法得到的平均偏差。