克里格法对隧道围岩稳定性分析

薛洪涛

(中铁十九局集团第七工程有限公司，广西 百色 533000)



克里格法对隧道围岩稳定性分析

薛洪涛

(中铁十九局集团第七工程有限公司，广西 百色 533000)

隧道围岩稳定性一直是隧道施工控制的重点，采用克里格算法对隧道衬砌结构稳定性计算进行了研究。根据不同的变异函数对随机变量特征进行表达，结合二次正交组合法，建立了隧道衬砌结构功能函数的插值算法，避免了隐式函数无法直接求解的问题；对隧道力学参数的取值,提出了基于超前地质预报中的地震波法获取方法，克服了隧道力学参数取值的模糊性。采用克里格法与采用传统蒙特卡洛方法的对比结果表明：迭代次数减少,失效概率的绝对误差仅为0.0152%，相对误差为4.26%。证明了该方法的合理性。

克里格法；隧道；二次正交组合法；超前地质预报

引 言

隧道作为很多公路或铁路的控制性工程，其围岩稳定性是重点，一旦出现工程事将损失巨大。研究者采用理论方法、现场监测和数值模拟等多种手段对隧道围岩稳定性进行了大量研究，得到的成果颇多。对于理论算法，1951年，南非地质学者Danie Krige提出了克里格模型[1-2]。该模型采用随机过程算法，能够得到方差估计较小的无偏估计，且具有半参数化的特点[3-4]，本文拟采用该方法应用于隧道围岩稳定性分析。

对于隧道力学参数的取值，计算时往往采用室外取试样后在室内做试验确定，得到的参数往往存在误差。而且隧道赋存的岩土环境本来就复杂多样，隧道的力学参数在同一断面的不同位置及不同断面的相同位置，力学参数都是变化的。随着隧道超前地质预报的发展，出现了一系列新技术，比如地质雷达、钻孔摄像和地震波法等手段。特别是地震波法，即采用瑞士的TSP203 plus，可获取测试掌子面前方岩体的泊松比μ、内摩擦角φ、弹性模量E和粘聚力c等一系列力学参数，其他参数也可经过计算得到[5]。这为隧道围岩稳定计算提供了新思路，即根据获得的实测数据分析围岩稳定，有更高的精准性与实时性[6]。

根据不同的变异函数对随机变量特征进行表达，结合二次正交组合法，建立了隧道衬砌结构功能函数的插值算法，避免了隐式函数无法直接求解的问题，超前地质预报中的地震波法获取隧道力学参数对隧道衬砌可靠度进行计算。

1 衬砌结构稳定功能函数

隧道结构主动压力经历弹性阶段、弹塑性阶段、塑性-分离阶段和松动阶段四个变形阶段，如图1所示。实践表明，D点是支护的最佳时刻，此时围岩压力PD达到最小，用Pamin表示。

图1 围岩形变与压力发展过程

Pamin= γr0[Rmax/r0-1]

(1)

式中，Rmax为最大的围岩允许松动区半径，γ为岩体重度，r0为隧道等代圆半径。

根据卡斯特纳方程，围岩最大松动圈半径Rmax：

(2)

式中，c1、φ分别为已开挖破裂区围岩的粘聚力和内摩擦角，σz为侧压系数λ=1的原岩应力。

c1=c+τaAs/SaSb

(3)

式中：c为围岩本身粘聚力，τa为锚杆抗剪强度，As为锚杆的横截面积，Sa,Sb分别为锚杆的横向和纵向间距。

根据可靠度有关定义，则极限状态方程为：

Z=Pi-Pamin=0

(4)

式中，Pi为支护结构即衬砌提供的支护阻力。

由式(1)～式(3)可知，Pamin为γ、r0、φ、c、σz、Sa、Sb、Pamin的函数，即

Pamin=f(γ,r0,φ,c,σz,Sa,Sb,Pamin)

(5)

Pamin为隐式函数，在确定性分析中，只能采用迭代方法求解。另外，由于式(5)不能表达为基本参数的明晰解析式，也就不可获取Pamin的概率分布，从而无法采用直接概率积分及已有的可靠度计算方法，必须寻找等效计算方法。本文引入克里格插值模型解决这一问题。

2 克里格插值算法

2.1 变异函数

设Z(x)为一随机函数，x∈Rn,x1,x2, ... ,xm是n维空间中样本点的位置，Z(x1)，Z(x2)…，Z(xm)为相应的联合分布的随机变量。克里格根据随机变量变异的空间位置特征，提出采用变异函数γ(h)进行模拟。在有限样本的情况下，γ(h)值可以计算：

(6)

Nh指在(xi+h,xi)之间用来计算n维空间中样本的变异函数值的样本对数，下标h表示Nh是分离距离h的函数；Z(xi+h)指点xi偏离h处随机变量的实现值；Z(xi) 指点xi处随机变量的实现值；h为步长，为样本点的空间间隔距离。

变异函数可用于表征随机变量的空间变异结构或空间连续性，目前比较常用的有球形模型、指数模型、高斯模型和幂函数模型等。采用式(6)计算得出h及γ(h)值，做出h-γ(h)试验变异图。在试验变异图的基础上，再对其拟合一个最优的理论变异函数模型。

2.2 算法

普通克里格法估计的一般公式：

(7)

式中，Z(xi)为随机函数在xi位置的已知值；Z*(xi)为在x0位置的随机函数估计值；n为估计过程中已知值的个数；λi为分配给Z(xi)残差的权重。

使用式(6)进行插值估计时，克里格法需满足2个要求：(1)选取Ki，使Z*(x0)的估计无偏；(2)估计误差的方差σz最小。

根据要求(1)有：

E[Z*(x0)-Z(X0)]=0

(8)

将式(8)代入式(7)，经整理可得：

(9)

根据要求(2)有：

σ2=Pamin=g(γ,r0,φ,c,σz,Sa,Sb)

(10)

将式(7)代入(8)，简化得：

(11)

式中，γ(h)为变异函数，由式(7)可得λi，根据拉格朗日原理和驻点的一阶偏导数原理，对式(11)的极值运算可得：

(12)

式中，μ为拉格朗日乘数，λi满足：

(13)

由式(12)、式(13)得：

(14)

式中，γij=γ(xi-xj)即距离为xi和xj之间的变异函数值，式(14)可简化为

[F][L]=[B]

(15)

式中，[F]为式(14)左边的系数矩阵，[L]为式(14)左边的待求解的向量，[B]为式(14)右边的已知向量。解这些线性方程组，即可得到所有的权重1,..., λn和拉格朗日乘数μ。

由计算所得的权重和拉格朗日乘数，克里格法的估计方差S可通过式(7)求得估计值Z*计算：

(16)

3 结合二次正交法的克里格插值可靠度计算

3.1 试验设计抽样方法

二次回归方程的一般形式为：

k=1，2，…，m-1(j≠k)

(17)

其中，a,{bj},{bkj},{bjj}为回归系数，可以看出该方程共有 (m+1)(m+2)/2项，因此试验次数n≥(m+1)(m+2)/2。

二次回归正交组合设计的总试验次数为：

n=mc+mγ+m0

(18)

根据正交性的要求，可以推导出星号臂长度γ必须满足关系式：

(19)

可见，星号臂长度γ表征试验区间，与因素数m、零水平试验次数m0及二水平试验数mc有关。

3.2 确定合适的二次回归正交组合设计

根据因素m选择合适的正交表进行交换，明确二水平试验方案，确定二水平试验次数mc和星号试验次数mγ。

根据二次回归正交组合设计表确定的试验方案，共进行n次试验，得到n个试验指标。有多少个未知数，便至少需进行多少次试验计算回归系数，建立含规范变量的回归方程。回归系数的计算：

常数项：

一次项偏回归系数：

交互项偏回归系数：

二次项偏回归系数：

总平方和：

自由度：dfT=n-1

一次项偏回归平方和：

交互项偏回归平方和：

二次项偏回归平方和：

各种偏回归平方和的自由度都为1。

残差平方和：

SSe=SST-SSR

自由度：

dfe=dfT-dfR

4 可靠性计算

4.1 TSP软件获取岩石力学参数

TSP203探测是根据人工制造一系列地震波的回波原理，在以开挖的隧道边墙布置若干个炮眼并置炸药半截与孔内，通过瞬发电雷管引爆每个药孔形成一系列的轻微震源，地震波经过掌子面前方返回被接收元件接受并传输于电脑。经处理分析软件可分析出掌子面后未开挖不良地质体，如：软弱破碎带、溶洞和沟槽等不良地质构造，同样可测掌子面为参考，根据波速、动态模量等一系列岩体力学参数判断将开挖里程段的围岩级别[5-6]。

提取的反射层可以显示小范围内岩石构造的显著变化。岩石特性的大幅度变化(与地震的主波长相当，主波长即波速除以主频)将通过速度函数，以较简单、可靠的方式进行评估。TSPwin通过纵横波的速度计算岩石力学参数，对于某些参数，则必须借助于经验关系式[7-8]。可从速度剖面直接计算出纵横波速度比和泊松比，TSP win中的岩石密度ρ通过纵横波速度的经验公式计算，求得密度后，可以计算出体积模量k、剪切模量μ、拉梅参数λ以及动态杨氏模量Edyn，纵横波速度比和泊松比是无量刚的，其他参数可由公式推导[6-8]。对于变质岩、火山岩、火成岩、沉积岩四大岩类，不同类型的岩石使用不同的经验公式。对于所有的由经验公式计算出的弹性参数，取由两种公式所求得结果的平均值。

计算过程中，不能由TSP获得的参数，即粘聚力、内摩擦角、锚杆的横向和纵向间距，需要根据地质资料、围岩分类和反分析来确定力学参数。等代圆半径r0变异性不大，取隧道等代圆半径即可，而σz原岩应力由隧道埋深h计算，参数变化见表1。

表1 围岩基本参数统计特征

4.2 克里格法的实现过程

计算方法1：

计算方法2：

选取普通克里格方法建立功能函数Pamin的插值模型；随机从随机变量(x1,x2, ... ,xm)中抽取一组变量代入Pamin的插值模型中计算对应的Pamin值，再将Pamin的值代入极限状态方程中求出功能函数Z值，重复n次，即可获得ζ1,ζ2...ζn的n个值计算出其均值和方差。

(20)

(21)

由可靠度定义式求出可靠度指标。

因为Pamin是关于γ、r0、φ、c、σz、Sa、Sb、Pamin的函数，一共含有7个随机变量。若在计算时同时考虑每一个变量，则计算工作量较大，过程繁琐。在满足实际工程需要的前提下，可以从中选取若干因素，将其视为随机变量进行分析。在利用TSP系统进行地质预报的过程中，可以获得岩体密度、泊松比以及静和动态弹性模量等参数，部分其他参数取自地质资料、公式计算和反分析，其中可以用作可靠度计算的参数为岩体密度γ(kN/m3)。

运用二次回归正交试验设计方法和正交试验设计方法，对将岩体密度视为随机变量γ取60个水平进行抽样做可靠度分析[11]，并将结果进行比选。

二次回归正交方法抽样的样本方差为：

Var(f(x))=0.4743

正交试验设计方法抽样的样本方差为：

Var(f(x))=1.3743

二次回归正交方法能够取得较小的样本方差，因此更适合本文分析。

用二次回归正交方法分别获取60、80、100、120个样本点计算结果。

由表2可以看出，当样本点分别为20、30、40、50与60时，其失效概率及可靠度的计算结果大体相近。为了验证克里格法的正确性，将60个样本点的计算结果与蒙特卡洛算法的结果进行比较，其中克里格法与蒙特卡洛法的误差为：

(22)

相对误差为：

(23)

其中，Pfk为克里格法的误差，Pfm为蒙特卡洛法的误差。由计算结果可知选用克里格法进行可靠度分析与已有方法结果较为接近。

表2 不同样本点时的可靠度比较

5 结束语

在隧道结构稳定功能函数的基础上，引入了克里格插值算法，并结合二次回归正交组合统计，建立了基于克里格算法的隧道围岩稳定性评价方法。

利用TSP获得的围岩力学参数，对实例围岩稳定性进行了可靠度分析。结果表明，当样本点个数为20、30、40、50、60时，其围岩失效概率与可靠度结果较为接近。同蒙特卡洛法进行比较时，其误差在5%之内，表明采用克里格法进行围岩可靠度计算的方法是可靠的。

[1] 苏永华,罗正东,张盼凤,等.基于Kriging的边坡稳定可靠度主动搜索法[J].岩土工程学报,2013(10):1863-1869.

[2] 靳国栋,刘衍聪,牛文杰.距离加权反比插值法和克里格插值法的比较[J].长春工业大学学报,2003,24(3):53-57.

[3] 张团峰,王家华.试论克里格估计与随机模拟的本质区别[J].西安石油学院学报,1997,12(2):52-55.

[4] 张崎,李兴斯.基于Kriging模型的结构可靠性分析[J].计算力学学报,2006,23(2):175-179.

[5] 张平松,吴健生.中国隧道及井巷地震波法超前探测技术研究分析[J].地球科学进展,2006,21(10):1033-1038.

[6] 杨正华,张文波,王卫东.浅析地震波法用于隧道病害的诊断与预测[J].灾害学,2008(1):27-30.

[7] 李坚.TSP 法在铁路客运专线隧道超前地质预报工作中的应用前景[J].铁道勘察,2006,31(6):45-49.

[8] 李华,李富,鲁光银.TSP 法与探地雷达相结合在隧道超前地质预报中的应用研究[J].工程勘察,2009,37(7):86-90.

[9] BOUCNEAU G,MEIRVENNE V M,Thas O,et al.Integrating properties of soil map delineations into ordinary kriging[J].European Journal of Soil Science,1998,49:213-229.

[10] 张仁铎.空间变异理论及应用[M].北京:科学出版社,2005.

[11] REZAEE H,ASGHARI O,YAMAMOTO J K.On the reduction of the ordinary kriging smoothing effect[J].Journal of Mining and Environment,2011,2(2):102-117.

Analysis of the Stability of Surrounding Rock in Tunnel Based on Kriging Method

XUEHongtao

(No.7 Engineering China Rail Way 19th Bureau Group, Baise 533000, China)

The stability of surrounding rock in tunnel is always the key point during construction, which was studied through the application of Kriging algorithm. The characteristic of random variable was expressed according to different variation function. In combination with quadratic orthogonal combination method, the interpolation algorithm of the function of the tunnel lining structure was established, which avoid the problem that the implicit function cannot be solved directly. For the value of tunnel mechanical parameters, earthquake wave method based on advance geology forecast was put forward, which overcome the fuzziness of tunnel mechanical parameter values. At last, compared Monte Carlo method with kriging method, the results showed that the number of iterations was reduced, and the absolute error of the failure probability was only 0.0152%, the relative error was 4.26%. This method was proved to be rational.

Kriging method; tunnel; quadratic orthogonal combination method; advance geological prediction

2015-07-24

国家自然科学基金(41130742)

薛洪涛(1975-)，男，辽宁沈阳人，工程师，主要从事隧道施工技术方面的研究，(E-mail)1412887537@qq.com

1673-1549(2016)01-0071-05

10.11863/j.suse.2016.01.15

TU94

A