隧道反射波超前探测有限差分正演模拟与偏移处理

鲁光银，熊瑛，朱自强

(中南大学 地球科学与信息物理学院，湖南 长沙，410083)

摘  要：针对地下工程领域隧道不良地质灾害超前探测问题，构建断层、倾斜软弱夹层等地质灾害体的隧道介质模型，研究隧道反射波超前探测的波场传播机理与偏移成像问题。利用一阶速度－应力弹性波方程和高阶交错网格差分计算方法，导出隧道工程反射波探测数值模拟的差分计算格式、数值稳定性条件以及边界条件，对上述模型进行正演模拟，得到相应模型不同时刻的波场快照与地震合成记录；通过对这些模拟结果进行分析，验证高阶交错网格差分法在隧道波场正演中的有效性。在导出弹性波方程逆时延拓交错网格差分计算格式基础上，通过求解程函数方程，获取激发时间成像条件，进而得到隧道全波场逆时偏移算法及偏移处理程序，并应用于隧道波场偏移成像中。研究结果表明：采用高阶交错网格有限差分逆时偏移算法能使反射波归位, 绕射波收敛, 从而大大提高了隧道工程反射波探测的分辨率，为隧道超前探测资料的解释提供理论依据。

关键词：反射波法；交错网格；正演计算；逆时偏移

中图分类号：P631.4⁺25          文献标志码：A         文章编号：1672-7207(2011)01-0136-06

Detection simulation ahead of tunnel face and reverse-time migration with reflection wave method

LU Guang-yin, XIONG Yin, ZHU Zi-qiang

(School of Earth Sciences & Info-physics, Central South University, Changsha 410083, China)

Abstract: The propagation and imaging problem of detection ahead of mining tunnel with the reflection wave method in underground engineering were discussed. Based on the fact that the detection and forecasting for badness geological structures ahead of the tunneling are difficult, two tunneling conceptual geological models with such geological hazards in tunnels as fault, tilted weak interlayer in front of tunnel face were constructed. After high-order finite-difference format of wave equation, the free boundary conditions and artificial boundary conditions were deduced based on one-order velocity-stress elastic wave equation, the detection simulation ahead of tunnel face for these two models was implemented. Then, the wave field snapshots and forward records were acquired and discussed, and the validity of high-order finite-difference forward simulation method in tunnel wave field was illuminated. Staggered-grid finite-difference scheme for inverse time calculation was derived from Eikonal equation, the seismic travel time at each underground point in these two models was computed by upwind-finite-difference format, which was the excitation-time imaging condition in reverse-time migration. Then the reverse-time migration algorithm of tunnel wave field and migration processing program were gotten, and synthetic data of those models were applied. The wave-field characters and migrated sections were analyzed with different models using reflection wave methods. The results show that this migration algorithm can make the reflection wave return to the original position and diffraction wave converge, which offers the logical proof for reflection wave method for detection ahead of tunnel.

Key words: reflection wave method; staggered-grid; forward simulation; reverse-time migration

在隧道地铁等地下工程修筑过程中，由于掌子面前方地质资料不详而造成支护不当或支护不及时，经常遭遇坍塌、突水突泥、冒顶等不良地质灾害，极易造成重大人员伤亡和财产损失。地下工程不良地质灾害超前预报问题已经成为土木工程领域面临的重大技术难题^[1]。超前预报方法一般分为破坏性方法与非破坏性探测方法。破坏性探测方法如超前钻孔、超前导坑法等具有费时、费力等不足，严重影响工程进度，一般只是在地质情况特别复杂时偶尔使用^[2]。在非破坏性探测方法中，常用的方法有陆地声纳法^[3]、TVSP法^[4]、HSP方法^[5]、TRT方法^[6]、SSP方法^[7]、ISIS方法^[8]、TST方法^[9]，常用的设备有TSP设备^[10]和DTC-150 防爆超前探测仪^[11]等。由于隧道超前地质预报的特点决定了现场工作必须在狭小的隧道施工现场实时完成，因此，超前探测技术必须满足狭小场地条件下进行大距离、高精度、高效率探测要求^[12]，这就要求在对传统方法技术(测线与探测对象平行)进行改进(测线与探测对象垂直)的同时，进一步研究地震波在隧道复杂条件下的传播机理特性与偏移成像技   术^[13]。目前，国内外学者开展的这方面研究^{[4, 14-16]}要么针对简单的声波方程，无法模拟隧道复杂条件下地震波场传播机理，要么采用传统的标量波波场处理方法，由于隧道复杂波场难以完全分离，极易影响处理精度。为此，本文作者在构建断层、倾斜软弱夹层等不良地质灾害体的隧道介质模型后，从一阶速度-应力弹性波方程出发，采用交错网格高阶有限差分进行隧道复杂波场的正演模拟计算，并采用矢量波波场处理方法实现隧道复杂波场的逆时偏移成像，以提高隧道反射波超前探测的分辨率。

1  交错网格高阶有限差分格式

在二维非均匀各向同性介质中，一阶速度-应力弹性波方程可以表示为如下矢量形式：

              (1)

式中：；v_x和v_z分别为x和z方向的速度分量；为面(体)应力分量；A₁和A₂分别为仅与介质密度和弹性系数有关的系数矩阵，

           (2)

        (3)

和为拉梅常数，，；v_p和v_s分别为介质纵波与横波波速；为介质密度。

对于方程(1)，采用高阶交错网格进行有限差分数值求解。交错网格法采用的网格剖分方案如图1所示。将v_x定义在网格结点(i, j)上，和定义在网格结点(i+1/2, j)上，定义在网格结点(i, j+1/2)上，定义在网格结点(i+1/2, j+1/2)上。质点震动速度分量和定义在时间离散采样点t=(k+1/2)上；，和定义在时间离散采样点t=k上。设和分别为X和Z轴方向的空间采样间隔；为时间采样间隔；U和V分别代表速度分量v_x和v_z的离散值，P，Q和R分别代表应力，和的离散值，则可推导出方程(1)的交错网格任意偶数阶差分格式：

(4)

在求解方程(4)时，当空间满足4阶精度，时间满足2阶精度时，系数，。

图1  交错网格剖分示意图

Fig.1  Diagram of staggered-grid

2  稳定性条件与吸收边界条件

设置稳定性条件是高阶有限差分数值求解过程中必须考虑的问题。对于二维各向同性介质，要使方程(3)稳定，在空间满足4阶精度、时间满足2阶精度时，稳定性条件为：

≤

 

≤

              (5)

对于隧道超前探测弹性波数值模拟计算模型，若仅考虑深埋隧道，则计算区域外部边界都是人工截断边界，隧道空间为自由边界。人工截断边界将向外传播的弹性波完全反射回去，从而使正演模拟计算与逆时偏移结果与实际情况不符。为此，本文采用董良国等^[17-18]提出的吸收边界条件和加阻尼边界的复合边界条件。人工边界示意图如图2所示。在模型的上、下、左、右边界及左顶角、右顶角、左下角、右下角加入人工吸收边界条件；同时，在吸收边界内部增加衰减带，G为衰减因子，其表达式如下：

  0≤i≤20      (6)

模型中，设隧道边界范围为：47＜x≤53，0＜     z＜36，隧道宽7 m；掌子面坐标z＝36 m。则对于隧道区域自由边界上，隧道侧壁L和R边界满足：

                  (7)

N_x为x轴方向网格节点数，N_z为z轴方向网格节点数

图2  人工边界示意图

Fig.2  Diagram of artificial boundary

隧道掌子面F边界满足：

                  (8)

由方程(8)，有

           (9)

对于交错网格，v_x和v_z不可能同时处在自由边界，只需处理式(9)中1个方程即可。则可推导出隧道掌子面自由边界F满足差分格式：

      (10)

同理，可以求出隧道左、右侧壁L和R自由边界方程的差分格式。

3  逆时偏移处理

波动方程的逆时传播问题实际上是弹性波正演的反问题。由于隧道侧壁各检波点在各个时刻的波场值是已知的，假设在t=t_L(t_L为隧道侧壁各接收点的最大记录时间)时，隧道模型中除检波点外的各点位移和质点振动速度均为0，在t＞t_L 时，检波点和隧道模型中各点的波场矢量均为0，由此可以将波动方程的逆时偏移问题转化为如下边值问题：

≤

           (11)

式中：P为隧道超前探测的波场矢量；为检波点位置()接收到的矢量波场值；t为时间。采用与方程(4)相同的网格对弹性波方程进行高阶差分离散，即可得到逆时传播的交错网格任意偶数阶差分格式。

论文采用激发时间成像条件，地震波走时采用熊金良等^[19]提出的逆风差分算法求解二维各向同性介质中的程函数方程：

           (12)

其中：t为时间(s)；s(x, z)为二维模型的慢度分布(m^-1∙s)；x和z为空间坐标。求解方程(11)时，需考虑模型中隧道自由边界与吸收边界影响。

4  计算实例

为验证上述高阶有限差分数值模拟方法的有效性与正确性，设计垂直断层与倾斜软弱夹层模型，采用主频为90 Hz的正弦波作为震源，模拟地震波在隧道介质模型中的传播。

模型1(图3(a))为1个垂直断层模型，断层位置位于z=120 m处。震源设置于掌子面正前方，坐标为(150，62)，模拟断层上盘参数为：纵波波速v_p=3.688 km/s，横波波速v_s=2.098 km/s，密度ρ=2.50 g/cm³；断层下盘参数为：v_p=4.513 km/s，v_s=2.434 km/s，ρ=2.70 g/cm³。每道时间采样点数N_t=2 000，采样时间间隔∆t=0.1 ms，记录长度为0.20 s(下同)。

图3  断层模型正演记录及偏移结果

Fig.3  Simulation and migration results for fault model

图4  倾斜软弱夹层模型正演记录及偏移结果

Fig.4  Simulation and migration results for gradient feeble interlayer model

图3(b)和3(c)所示分别为断层正演计算过程中弹性波在0.04 s时的水平分量与垂直分量的波场。波场清晰反映出断层界面上的反射波与透射纵波、透射横波、瑞雷波等复杂波场。根据正演计算结果，采用激发时间成像条件进行逆时偏移，水平分量与垂直分量偏移结果如图3(d)和3(e)所示。从图3(d)和3(e)可见：偏移成像结果清晰，且准确反映断层界面位置(z=120 m)，但垂直面的分量偏移效果比水平分量的好。

模型2为倾斜软弱夹层模型(图4(a))，软弱夹层位于模型对角线附近，与隧道轴线呈45º夹角，宽40 m。软弱夹层介质参数为：v_p=2.384 km/s，v_s=1.330 km/s，ρ=2.17 g/cm³；软弱夹层前面的岩层参数为：v_p=3.688 m/s，v_s=2.098 m/s，ρ=2.50 g/cm³；软弱夹层后面的岩层参数为：v_p=4.114 m/s，v_s=2.320 km/s，ρ=2.60 g/cm³。

图4(b)和4(c)所示分别为弹性波在0.045 s时的水平分量与垂直分量的波场。可见：波场清晰反映出该倾斜软弱夹层两层界面上的多次反射波与透射波及隧道边界上传播的瑞雷波。利用正演模拟结果作为逆时偏移的初始条件，进行逆时偏移，图4(d)和4(e)所示分别为水平分量和垂直分量的逆时偏移结果。从图4(d)和4(e)可见：无论水平分量还是垂直分量，皆较准确地反映出倾斜软弱夹层界面；上界面能量较强，下界面能量较弱；水平分量较垂直分量强，界面更加清晰。同时，在软弱夹层Ⅰ的镜像部位出现假象，但其能量较弱，且越靠近左侧，能量越弱，据此可以判断该软弱夹层与隧道轴向的倾向关系(断层判断与此类似)。

5  结论

(1) 利用一阶速度-应力弹性波方程，进行高阶有限差分波动方程正演计算，能够较好地描述变化剧烈的速度间断面；在数值计算过程中，计算程序简单，易于实现，与二阶波动方程相比能够减少计算量，提高计算效率。

(2) 采用交错网格，外部边界采用吸收边界条件和加阻尼边界的复合边界条件，隧道内部采用自由边界，弹性波方程能够完好地模拟隧道超前探测中的断层、软弱夹层等不良地质灾害体，反映出隧道波场传播机理和传播特征，对于满足一定精度的正演模拟计算与逆时偏移问题，计算速度较常规网格快，网格频散明显减少。

(3) 采用激发时间成像条件进行弹性波方程的逆时偏移成像，可将隧道反射波能量归位到其空间的真实位置，获取隧道掌子面前方断层、倾斜软弱夹层等简单不良地质灾害体的真实构造图像。

(4) 当隧道掌子面前方存在倾斜软弱夹层(或断层)时，虽然在其隧道轴线镜像部位存在假象，但根据逆时偏移成像剖面能量分布情况，可以判断软弱夹层(或断层)倾向。

参考文献：

[1] 张平松, 刘盛东, 吴健生. 隧道及井巷工程超前探测模拟及其偏移技术研究[J]. 岩石力学与工程学报, 2007, 26(7): 2847-2851.
ZHANG Ping-song, LIU Sheng-dong, WU Jian-sheng. Study on detecting simulation ahead of tunnel and laneway and its migration technique[J]. Chinese Journal of Rock Mechanics and Engineering, 2007, 26(7): 2847-2851.

[2] Madani H. Tunneling[M]. 2nd ed. Persian: Tehran Polytechnic University Press, 1998: 15-70.

[3] 钟世航. 陆地声纳法的原理及其在铁路地质勘测和隧道施工中的应用[J]. 中国铁道科学, 1995, 16(4): 48-55.
ZHONG Shi-hang. The principle and application of land sonar method in railway geological survey and tunnel construction[J]. China Railway Science, 1995, 16(4): 48-55.

[4] 曾昭璜. 隧道地震反射法超前预报[J]. 地球物理学报, 1994, 37(2): 268-271.
ZENG Zhao-huang. Prediction ahead of the tunnel face by the seismic reflection methods[J]. Chinese Journal of Geophysics, 1994, 37(2): 268-271.

[5] Inazaki T, Isahai H, Kawamura S, et al. Stepwise application of horizontal seismic profiling for tunnel prediction ahead of the face[J]. The Leading Edge, 1999, 18(12): 1429-1431.

[6] Neil D, Haramy K, Hanson D, et al. Tomography to evaluate site conditions during tunneling[C]//3rd National Conference of the Geo-Institute, American Society of Civil Engineers. Geotechnical Special Publication, 1999: 13-17.

[7] Kneib G, Kassel A, Lorenz K. Automated seismic prediction ahead of the tunnel boring machine[J]. First Break, 2000, 18(7): 295-302.

[8] Borm G, Giese R, Klose C, et al. ISIS integrated seismic imaging system for the geologic prediction ahead in underground construction[C]//65th Annual Conference and Exhibition, EAGE, Extended Abstracts. Norway, 2003: 887-899.

[9] Zhao Y G, Jiang H, Zhao X P. Tunnel seismic tomography method for geological prediction and its application[J]. Applied Geophysics, 2006, 3(2): 69-74.

[10] Dickmann T, Sander B. Drivage concurrent tunnel seismic prediction[J]. Felsbau-rock and Soil Engineering, 1996, 14(1): 406-411.

[11] 梁庆华, 宋劲. 矿井多波多分量地震勘探超前探测原理与实验研究[J]. 中南大学学报: 自然科学版, 2009, 40(5): 1392-1398.
LIANG Qing-hua, SONG Jin. Advanced detection theory and experimental research of multi-wave and multi-component seismic exploration in mine[J]. Journal of Central South University: Science and Technology, 2009, 40(5): 1392-1398.

[12] 宋先海, 顾汉明, 肖伯勋. 我国隧道地质超前预报技术述评[J]. 地球物理学进展, 2006, 21(2): 605-613.
SONG Xian-hai, GU Han-ming, XIAO Bo-xun. Overview of tunnel geological advanced prediction in China[J]. Progress in Geophysics, 2006, 21(2): 605-613.

[13] 张平松, 吴健生. 中国隧道及井巷地震波法超前探测技术研究分析[J]. 地球科学进展, 2006, 21(10): 1033-1038.
ZHANG Ping-song, WU Jian-sheng. Research and analysis of forward prediction technology using seismic reflection wave in tunnel and laneway in china[J]. Advance in Earth Science, 2006, 21(10): 1033-1038.

[14] 何振起, 李海, 梁颜忠. 利用地震反射法进行隧道施工地质预报[J]. 铁道工程学报, 2000(4): 81-85.
HE Zhen-qi, LI Hai, LIANG Yan-zhong. Geological super leading forecast during tunnel construction by utilizing seismic response analysis methods[J]. Journal of Railway Engineering Society, 2000(4): 81-85.

[15] 常旭, 刘伊克, 桂志先. 反射地震零偏移距逆时偏移方法用于隧道超前预报[J]. 地球物理学报, 2006, 49(5): 1482-1488.
CHANG Xu, LIU Yi-ke, GUI Zhi-xian. Zero-offset reverse time migration for prediction ahead of tunnel face[J]. Chinese Journal of Geophysics, 2006, 49(5): 1482-1488.

[16] Bohlen T, Lorang U. Rayleigh-to-shear wave conversion at the tunnel face From 3D-FD modeling to ahead-of-drill exploration[J]. Geophysics, 2007, 72(6): T67-T79.

[17] 董良国, 马在田, 曹景忠. 一阶弹性波方程交错网格高阶差分法稳定性研究[J]. 地球物理学报, 2000, 43(6): 856-864.
DONG Liang-guo, MA Zai-tian, CAO Jing-zhong. A study on stability of the staggered-grid high-order difference method of first-order elastic wave equation[J]. Chinese Journal of Geophysics, 2000, 43(6): 856-864.

[18] Cerjian C, Kosloff D, Kosloff R, et al. A nonreflecting boundary condition for discrete acoustic and elastic wave equations[J]. Geophysics, 1985, 50(3): 705-708.

[19] 熊金良, 李国发, 王濮, 等. 地震波走时计算的逆风差分算法[J]. 物探化探计算技术, 2004, 26(4): 295-298.
XIONG Jin-liang, LI Guo-fa, WANG Pu, et al. Upwind finite difference solution of eikonal equation[J]. Computing Techniques for Geophysical and Geochemical Exploration, 2004, 26(4): 295-298.

(编辑 陈灿华)

                                 

收稿日期：2009-10-12；修回日期：2010-01-15

基金项目：国家自然科学基金资助项目(10871208)；湖南省科技计划项目(05SK3077)

通信作者：鲁光银(1976-)，男，湖北宜昌人，博士，副教授，从事隧道井巷地质灾害探测与评价技术研究；电话：13975894898；E-mail: luguangyin2006@126.com