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

(1. 成都理工大學(xué)環(huán)境與土木工程學(xué)院, 四川 成都 610059; 2. 成都理工大學(xué) 地質(zhì)災(zāi)害防治與地質(zhì)環(huán)境保護(hù)國家重點(diǎn)實(shí)驗(yàn)室, 四川 成都 610059)

0 引言

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

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

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

1 基于克里金法的深埋長隧道地應(yīng)力插值計(jì)算方法

1.1 克里金插值法

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

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

(1)

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

1.2 實(shí)驗(yàn)變異函數(shù)

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

(2)

式中:γij為半方差;Zi為i點(diǎn)的地應(yīng)力測量值;Zj為j點(diǎn)的地應(yīng)力測量值。

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

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

(3)

(4)

(5)

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

1.3 克里金方程組

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

(6)

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

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

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

圖1 基于克里金法的深埋長隧道地應(yīng)力插值計(jì)算流程

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

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

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

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

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

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

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

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

2 工程應(yīng)用

2.1 工程概況

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

圖2 隧道縱剖面及地應(yīng)力測點(diǎn)布置圖

2.2 交叉驗(yàn)證

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

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

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

(7)

(8)

(9)

(10)

式(7)—(10)中:Z(S1)為地應(yīng)力實(shí)測值;Z(S2)為插值估算值;n為測量樣本數(shù);δ為樣本標(biāo)準(zhǔn)差。

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

2.3 主應(yīng)力插值及分析

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

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

(11)

(12)

(13)

(14)

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

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

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

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

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

2.4 隧道軸線地應(yīng)力特征分析

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

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

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

σθ=k·σH。

(15)

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

(16)

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

圖4 最大水平主應(yīng)力等值線圖

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

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

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

(17)

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

(18)

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

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

圖5 隧道軸線最大水平主應(yīng)力值及方向

3 結(jié)論與建議

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

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

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

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