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

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

0 引言

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

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

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

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

1.1 克里金插值法

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

在隧道工程選線、勘察及設(shè)計階段,為了避免線路反復調(diào)整及隧道地質(zhì)信息和地應力實測數(shù)據(jù)缺乏等不良影響,將已有的地應力數(shù)據(jù)(不需要隧道地質(zhì)內(nèi)容)作為已知點,在隧道邊界范圍內(nèi),利用克里金基本公式能較快捷高效地求出隧道全域內(nèi)各點的地應力值,即可得到隧道軸線3個主應力結(jié)果?？死锝鸩逯捣ǖ幕竟饺缡?1)所示。

(1)

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

1.2 實驗變異函數(shù)

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

(2)

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

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

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

(3)

(4)

(5)

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

1.3 克里金方程組

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

(6)

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

1.4 隧道地應力插值方法的實現(xiàn)流程

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

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

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

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

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

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

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

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

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

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

2 工程應用

2.1 工程概況

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

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

2.2 交叉驗證

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

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

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

(7)

(8)

(9)

(10)

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

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

2.3 主應力插值及分析

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

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

(11)

(12)

(13)

(14)

表1 變異函數(shù)模型誤差分析結(jié)果

表2 高斯模型參數(shù)

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

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

圖3 最大水平主應力擬合結(jié)果

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

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

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

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

σθ=k·σH。

(15)

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

(16)

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

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

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

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

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

(17)

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

(18)

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

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

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

3 結(jié)論與建議

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

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

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

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