深埋長隧道縱剖面初始地應力插值計算方法及應用

任 洋, 吳岳華, *, 李天斌

(1. 成都理工大學環境與土木工程學院, 四川 成都 610059; 2. 成都理工大學 地質災害防治與地質環境保護國家重點實驗室, 四川 成都 610059)

0 引言

地應力是隧道及地下工程勘察、設計、施工過程中最為關鍵的因素之一。獲得隧道地應力最直接有效的辦法是現場地應力測試,但存在測試難度大、工期長、造價高等問題[1-2],并且僅靠少量實測點不能滿足深埋長隧道地應力實際需求。后來工程技術人員通過地質建模,采用數值反演獲得隧道地應力[3-4]; 另外,隨著人工智能和計算機技術的興起,地應力智能反演研究也越來越多[5-6],各種地應力反演和分析方法得到了很好的發展和應用[7-9]。但地應力數值反演需要比較詳細、準確的地質資料和巖土體參數,在隧道工程勘察和設計階段相關資料缺乏或不太準確,加之隧道線路和設計方案等修改頻繁,往往需要重新建立隧道地質模型并計算,耗費大量時間。因此,探尋一種能滿足隧道工程選線、勘察和設計等階段地應力使用需求和精度要求,且能快捷高效地獲得隧道軸線初始地應力的方法是十分必要的。

地質統計分析及空間插值方法不受地質資料等限制,可以快捷高效地獲得隧道軸線地應力結果。克里金插值方法是諸多地質統計及插值方法中常用的方法,其能直接給出插值誤差并評價插值的準確性,效率較高,實用性較強[10-11]。劉亞靜等[12]采用克里金插值方法估計礦體儲量; 劉軒等[13]運用普通克里金法對同位素測氡探火數據進行優化處理,提高對火源位置判斷的準確性; 楊雪峰等[14]采用克里金法獲得海域三維溫度場的空間分布特征; 鄧朗妮等[15]采用優化克里金插值方法計算土方量; 徐池等[16]選取球狀模型作為變異函數,采用優化泛克里金法對短波頻率進行重構,豐富了軍事通訊的技術; Yanto等[17]運用克里金插值方法判斷滑坡的易發區域,為滑坡治理提供參考; Hao等[18]對克里金插值方法進行了改進,并準確預測了儲層孔隙度; Mahdi等[19]運用已有的84個數據構建克里金模型,得到估算巖石節理抗剪強度的最佳模型; Liu 等[20]利用滑坡位移監測數據實現了時空最優組合,建立了時空變形場的預測模型,為滑坡預警提供了科學依據。克里金插值方法在地質統計中的應用較廣泛,但目前對于地形起伏深埋長隧道的地應力插值計算研究極少。

本文構建了一種基于地質統計克里金法的深埋長隧道縱剖面初始地應力插值計算方法。該方法選取最優變異函數模型進行地應力插值計算,獲得隧道縱剖面主應力量值及最大水平主應力方向; 再考慮隧道剖面內不同地層及斷層的影響,對插值計算結果進行修正,最終獲得隧道洞軸線的主應力量值及最大水平主應力方向。

1 基于克里金法的深埋長隧道地應力插值計算方法

1.1 克里金插值法

克里金插值法是由南非工程師D. G. Krige提出的一種空間插值方法。其基本假設是一點的屬性值與周圍點的屬性值有關,并且可以由其周圍點的屬性值推導得出,它是以變異函數為計算工具,并結合結構分析的一種空間相關性很強的最優、無偏估算方法[21]。

在隧道工程選線、勘察及設計階段,為了避免線路反復調整及隧道地質信息和地應力實測數據缺乏等不良影響,將已有的地應力數據(不需要隧道地質內容)作為已知點,在隧道邊界范圍內,利用克里金基本公式能較快捷高效地求出隧道全域內各點的地應力值,即可得到隧道軸線3個主應力結果。克里金插值法的基本公式如式(1)所示。

(1)

式中:Z(X0)為估測值;λi為第i個權重系數;Z(Xi)為已知點的屬性值。

1.2 實驗變異函數

實驗變異函數是克里金插值方法的重要研究工具。2個點實測地應力值之差的平方的一半被稱為這2個點的半方差,即

(2)

式中:γij為半方差;Zi為i點的地應力測量值;Zj為j點的地應力測量值。

在n個地應力測點中不重復地選取2個測點的地應力值計算半方差,任意2個點的半方差值與距離相關,通過求得所有測點兩兩之間的半方差與距離,得到C(n,2)個與距離相關的半方差。其中,C(n,2)為組合的線性寫法。利用C(n,2)個半方差值與C(n,2)個對應的距離,擬合實驗變異函數。

實驗變異函數常用的計算模型有[22]: 球狀模型、指數模型和高斯模型,其表達式分別為式(3)、式(4)和式(5)。

(3)

(4)

(5)

式(3)—(5)中:γ(h)為實驗變異函數;C0為基臺值;C為拱高,表示區域化變量在空間上變化的極大值;a為變程,表示區域化變量具有關聯性的范圍;h為區域化變量到待估測點的距離。

1.3 克里金方程組

克里金方程組包含n+1個方程,用以求解n個權重系數,該方法無偏估計的條件是所有權重系數λi的和為1。運用拉格朗日數乘法構造求解函數并得到克里金方程組,如式(6)所示。

(6)

式中:γ(Xi,Xj)為第i個與第j個已知測點的半方差值;γ(Xi,X0)為估算點與已知測點i的實驗變異函數值;u為拉格朗日乘數;λi為搜索領域內各測點對估算點的權重系數。

1.4 隧道地應力插值方法的實現流程

基于克里金法的深埋長隧道地應力插值計算流程如圖1所示,其具體實現流程如下。

圖1 基于克里金法的深埋長隧道地應力插值計算流程

1)數據篩選及處理。收集深埋長隧道縱剖面及地應力測點資料(不需要地質內容)。對實測地應力數據進行統計分析和LOG或BOX-COX變換,使其服從正態分布。將地應力數據代入式(2)計算兩兩測點的半方差,得到C(n,2)個與距離相關的半方差,用于擬合實驗變異函數。

2)擬合實驗變異函數。確定C(n,2)個距離中最長的距離Hmax,將C(n,2)個數據按距離Hmax/H1分組,計算各組的平均距離及平均半方差,得到Hmax/H1組平均距離與平均半方差的關系,用最小二乘法擬合實驗變異函數。

3)求解權重系數。得到實驗變異函數后即可得到待估算點與已知點的半方差,在以待估算點為中心的搜索范圍內,計算估算點與所有已知點兩兩的半方差,再結合第2)步計算得到所有已知點之間的半方差值,代入式(6)中建立m+1個方程,并計算m個權重系數(m為以待估算點為中心,搜索范圍內所有已知點的個數)。

4)求估算值。將求解得到的m個權重系數代入式(1)進行插值計算,得到該點地應力插值計算結果。

5)優化各個變異函數模型。滯后距H1的選擇將影響變異函數的質量,不斷改變滯后距H1,直至得到最優的實驗變異函數。

6)誤差分析。模型誤差分析的目的是在常用的3種實驗變異函數模型中確定最優變異函數模型作為最終計算模型。在各測點處重復第3)、4)步驟得到各測點處的估算值,利用測點處估算值和測量值計算出評估變異函數模型質量的相關指標(平均誤差θ、均方根誤差RMS、標準化均方根誤差RMSS和平均標準誤差ASE),選取地應力擬合精度最高的變異函數模型。

7)結果分析。選取最優變異函數模型的輸出結果進行分析,并將插值計算結果與實測地應力值及地應力方向進行擬合對比。目前常用的地應力數值反演和智能反演等方法的誤差在20%左右[23-24],本文以相對誤差在20%以內作為地應力插值計算的精度。

8)提取隧道軸線的初始地應力量值及最大水平主應力方向,考慮不同地層和斷層的影響,對隧道軸線地應力插值結果進行修正,獲得隧道洞軸線地應力分布特征。

2 工程應用

2.1 工程概況

郭達山某深埋長隧道軸線長約10 km,最大埋深約1 760 m,隧道軸線范圍內共布置12個地應力測點,如圖2所示。

圖2 隧道縱剖面及地應力測點布置圖

2.2 交叉驗證

用球狀模型、指數模型和高斯模型分別計算最大水平主應力SH、最小水平主應力Sh、垂直主應力Sv及最大水平主應力方向的實驗變異函數,通過不斷改變滯后距H1計算出最大水平主應力SH、最小水平主應力Sh、垂直主應力Sv及最大水平主應力方向的最優球狀模型、指數模型和高斯模型。

采用交叉驗證來檢驗各模型的質量。交叉驗證會移除樣本中的1個測點,利用剩下的測點計算出移除點的估算值,重復該操作直至獲得所有測點的估算值。這樣可以用估算值和測量值計算出評估變異函數模型質量的相關指標。

評估變異函數模型質量的指標包括: 平均誤差θ、均方根誤差RMS、標準化均方根誤差RMSS和平均標準誤差ASE,計算公式分別見式(7)—(10)。平均誤差越接近于0、標準化均方根誤差越接近于1、平均標準誤差與均方根誤差越接近,說明變異函數模型的質量越好[25]。

(7)

(8)

(9)

(10)

式(7)—(10)中:Z(S1)為地應力實測值;Z(S2)為插值估算值;n為測量樣本數;δ為樣本標準差。

在球狀模型、指數模型和高斯模型中分別計算出最大水平主應力SH、最小水平主應力Sh、垂直主應力Sv及最大水平主應力方向D的誤差結果,如表1所示。在3個模型中,3個主應力及最大水平主應力方向D的平均誤差均接近于0。高斯模型中的標準化均方根誤差最接近于1,且均方根誤差與平均標準誤差最接近,所以高斯模型是本案例中最優變異函數模型。

2.3 主應力插值及分析

采用交叉驗證迭代不斷更新滯后距H1,并用最小二乘法擬合計算最大水平主應力SH、最小水平主應力Sh、垂直主應力Sv及最大水平主應力方向D的高斯模型的各參數,計算得到主應力變異函數模型的參數,如表2所示。

將參數代入模型中得到最大水平主應力SH、最小水平主應力Sh、垂直主應力Sv及最大水平主應力方向D的變異函數模型(高斯模型)公式中,如式(11)—(14)所示。

(11)

(12)

(13)

(14)

表1 變異函數模型誤差分析結果

表2 高斯模型參數

由于數據點過多,限于篇幅僅列出埋深較大的③、④、⑧號3個鉆孔的地應力測量值、插值數據及對應的擬合偏差度,如表3所示。最大水平主應力擬合結果如圖3所示。其中,擬合偏差度=(插值結果-測量值)/測量值。

將以上3個鉆孔共29個測點的最大水平主應力SH、最小水平主應力Sh和垂直主應力Sv的插值數據與實測數據進行對比可知: 1)大部分測點的最大水平主應力值擬合偏差度為[-0.15,0.15],誤差控制在±15%以內,擬合度約為85%,僅鉆孔淺部的個別測點擬合度為76%。初步分析認為鉆孔淺部地應力受斜坡地形及測試結果本身誤差等影響。2)大部分測點的最小水平主應力值擬合偏差度為[-0.20,0.15],鉆孔淺部個別測點擬合度較差。3)大部分測點的垂直主應力值擬合偏差度為[-0.20,0.10]。

圖3 最大水平主應力擬合結果

2.4 隧道軸線地應力特征分析

最優實驗變異函數插值結果與實測地應力值的擬合度較好(12個測孔的平均擬合度超過80%,其中65個測點的擬合度大于85%)。

獲得隧道縱剖面主應力量值結果,最大水平主應力SH等值線見圖4。

從圖4中提取出隧道軸線最大水平主應力值,考慮不同地層和斷層的密度及巖體儲能能力的差異,構建應力修正關系式,如式(15)所示。修正系數k與不同地層及斷層的彈性模量和密度有關,構建相應的關系式,如式(16)所示。根據相關試驗及科研成果,隧道軸線不同地層及斷層的物理力學參數如表4所示。

σθ=k·σH。

(15)

k=a′·E+b′·ρ。

(16)

式(15)—(16)中:σθ為修正后的主應力值;σH為最大水平主應力;a′為巖體彈性模量對k的貢獻系數,為常數;b′為巖體密度對k的貢獻系數,為常數;E為巖體彈性模量;ρ為巖體密度。

圖4 最大水平主應力等值線圖

表4 隧道軸線不同地層及斷層的密度和彈性模量

將郭達山某隧道地應力數值計算結果與本文的插值結果進行對比,得到隧道軸線各地層初始地應力修正公式(如式(17)所示)及斷層處初始地應力修正公式(如式(18)所示)。

σ1=(0.004 2E+0.000 3ρ)σH。

(17)

σ2=(5.798 6E-0.004 7ρ)σH。

(18)

式(17)—(18)中:σ1為各地層的主應力修正值;σ2為斷層處的主應力修正值。

根據以上公式對隧道軸線地應力進行修正,獲得隧道軸線最大水平主應力結果。繪制隧道軸線最大水平主應力值及方向,如圖5所示。隧道軸線最大水平主應力值達52.50 MPa,斷層處地應力明顯較低; 最大水平主應力方向以NW向為主,局部段為NE向,其中,NW向的方位角為N49.70W～N83.78W,NE向的方位角為N58.31E～N79.60E。

圖5 隧道軸線最大水平主應力值及方向

3 結論與建議

針對隧道工程地應力的使用需求,構建了一種深埋長隧道縱剖面初始地應力插值計算方法,并進行了工程應用,獲得了以下主要結論:

1)依據少量的實測地應力數據,通過最優變異函數模型進行插值計算,考慮不同地層及斷層的影響,對地應力插值結果進行修正,獲得隧道軸線主應力量值及最大水平主應力方向。深埋長隧道縱剖面初始地應力插值計算方法不受勘察資料及線路和設計方案往返調整等影響,可快速獲得隧道軸線的地應力結果。

2)對多個計算模型的質量進行評價,采用交叉驗證迭代計算出各指標的誤差值,在多個變異函數模型中選出最優函數模型用于地應力插值計算。通過工程實例應用可知,插值計算得到的地應力值與實測地應力值的擬合度較高,其中最大水平主應力的擬合度約為85%。

本文主要研究隧道縱剖面的地應力插值計算方法,僅考慮了隧道軸線內不同地層和斷層對地應力的影響,獲得的是隧道軸線地應力值及其方向。隧址區地應力分布呈三維狀態,且各地層和斷層的空間分布對地應力有一定影響,建議后續研究中考慮地質條件及空間分布等因素,以獲得隧道地應力的空間分布規律。