深埋长隧道纵剖面初始地应力插值计算方法及应用

任 洋, 吴岳华, *, 李天斌

(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%。

本文主要研究隧道纵剖面的地应力插值计算方法,仅考虑了隧道轴线内不同地层和断层对地应力的影响,获得的是隧道轴线地应力值及其方向。隧址区地应力分布呈三维状态,且各地层和断层的空间分布对地应力有一定影响,建议后续研究中考虑地质条件及空间分布等因素,以获得隧道地应力的空间分布规律。