基于VTI介质的隧道超前探测qP波正演模拟研究.pdf

第 49 卷第 3 期 当 代 化 工 Vol.49，No.3 2020 年 3 月 Contemporary Chemical Industry March，2020 基金项目基金项目 国家自然科学基金项目；国家重点研发计划 基金项目，项目号41674107；2017YFB0202904。 收稿日期收稿日期 2019-11-27 作者简介作者简介 董载天（1995-） ，男，山东省济南市人，就读长江大学地球物理与石油资源学院，研究方向工程勘探。E-mail506501452。 通讯作者通讯作者王军民（1960-） ，男，博士，研究方向地球物理仪器开发与研究。E-mailwjm8510。 基于 VTI 介质的隧道超前探测 qP 波正演模拟研究 董载天，王军民 （长江大学， 湖北 武汉 430100） 摘 要 由于隧道超前探测的环境与常规石油地震勘探的差异，不宜直接将过去的模拟方法拿到隧道中 进行应用。 结合层状地层表现出横向各向同性的特性， 推导出在隧道空间适用的二维 VTI 介质一阶 qP 波方程， 分析了其稳定条件。建立断层和软弱夹层的几种数值模型，利用 PML 边界解决模型边界反射，采用交错网格 进行有限差分正演模拟，得到其共炮点道集。结果表明，qP 波有限差分法可以快速地模拟不良的地震响应，规 避伪 S 波和数值不稳定问题，为隧道超前探测资料的解释提供依据，提高预测的准确性。 关 键 词隧道超前探测；VTI 介质；qP 波；PML 中图分类号TV221.2 文献标识码 A 文章编号 1671-0460（2020）03-0683-07 Forward Modeling of Tunnel Advanced Detection qP Wave Based on VTI Medium DONG Zai-tian, WANG Jun-min （Yangtze University, Hubei Wuhan 430100, China） Abstract Due to the difference between the environment of tunnel pre-detection and conventional petroleum seismic exploration, it is not appropriate to directly apply the past simulation to the tunnel. Combined with the transverse isotropy characteristic of layered stratum, the first-order qP wave equation of the two-dimensional VTI medium suitable for tunnel space was derived, and the stability conditions were analyzed. Several numerical models of fault and weak interlayer were established. The boundary reflection of the model was solved by PML boundary, and the finite difference forward modeling was carried out by staggered grid to obtain the common shot gather. The results show that the qP wave finite difference can quickly simulate the bad seismic response, avoid the pseudo S wave and numerical instability problem, and provide a basis for the interpretation of tunnel advance detection data to improve the accuracy of prediction. Key words tunnel advance detection; TI media; qP wave; PML 在我国经济高速发展的当下，大量公路铁路工 程都有隧道施工。在隧道工程施工过程中，由于掌 子面前方地质状况不详，面对破碎带、裂隙和溶洞 的支护不当或支护不及时，经常出现坍塌、突水突 泥、冒顶等重大灾害，造成严重的人员伤亡和财产 损失。显然在隧道工程中，使用超前探测技术了解 掌子面前方地质状况，对工程施工安全极其重要。 超前探测方法从原理上可以划分为地震波法、电 磁法（瞬变电磁、电阻率、激发极化、核磁共振） 、 岩体温度法等，其中地震波法在远距离勘探中有较 好的分辨率，成为目前隧道超前探测中的常用方法 之一。 在隧道地震波法探测中， 常用的方法有 VSP、 HSP ［1］ 、TSP ［2］ 、TVSP ［3］ 、陆地声纳 ［4］ 、HSP ［5］ 、 TRT ［6,7］ 、TST ［8］ 、ISIS ［9］等方法。由于隧道施工现 场的空间狭小且工程工期压力大，所以隧道超前探 测技术必须满足狭小场地条件下进行大距离、高精 度、高效率探测要求 ［10］ ，这就要求在对传统方法技 术（测线与探测对象平行）进行改进（测线与探测 对象垂直）的同时，进一步研究地震波在隧道复杂 条件下的传播机理特性 ［11］ 。因此，本文在层状地层 条件下，构建断层、软弱夹层等不良地质灾害体模 型，从一阶速度−应力 qP 波方程出发，采用交错网 格高阶有限差分进行隧道复杂波场的正演模拟计算， 以充分认识不同条件下地震波响应特征。 1 TI 介质的波动方程 在我国，层状地层分布广泛，约占陆地面积的 77.3 ［12］ ，在隧道施工过程中大多要面对这样的地 质环境，而地下层状在横向各向同性地层往往表现 出横向各向同性的物理响应特征。介于笔者使用波 DOI10.13840/21-1457/tq.2020.03.041 684 当 代 化 工 2020年3月 动方程法进行地震数值模拟，所以推导出横向各向 同性 （TransverseIsotropy） 介质的波动方程尤为重要。 1.1 波动方程导出 波动方程的建立涉及应力、应变和位移这三个 物理量，需要使用弹性体的本构方程、运动微分方 程和柯西方程。由本构方程，可知弹性体应力与应 变关系表达式 , , ,1,2,3 ijijklkl C ei j k l  （1） 式中 应力张量; e 应变张量; C 弹性常数。 弹性常数共有 81 个，由柯西方程可知，应变 张量有对称性 ijji uv ee yx    （2） 其中u、v 位移分量； x, y 笛卡尔坐标。 而应力张量也有此种对称性，即。因此弹性常 数减少为 36 个， 同时弹性常数张量相对于主对角线 对称，有。使用单下标替代双下标 11→1,22→2,33 →3,23（32）→4,13（31）→5,12（21）→6，则可 的表达式 11112131415161 22122232425262 33132333435363 44142434445464 55152535455565 66162636465666 CCCCCCe CCCCCCe CCCCCCe CCCCCCe CCCCCCe CCCCCCe                            （3） 上式可简记为 C e （4） 从弹性体的角度，描述完应变与应力关系、应 变与位移关系后，通过牛顿第二定律可建立运动平 衡微分方程 2 2 U LF t     （5） 式中 弹性体密度，kg/m; U 位移张量; 应力张量; t 时间，s; F 体力张量; L 偏导算子。 偏导算子L为 000 000 000 xzy L yzx zyx              （6） 利用 Cauchy 方程可将（5）式表达为 2 T 2 U L CL UF t     （7） 1.2 TI 介质中 qP 波波动方程 地表岩土物理性质极其复杂不便于数值模拟 研究，地球物理学家为了能借助波动物理研究其波 传播特性，根据各项异性弹性体介质的基本对称性 进行了分类。其中三方、六方各项异性介质常叫做 横向各向同性 （Transverse Isotropy， 简称 TI） 介质， 是目前实际应用研究中使用最广泛的模型。TI 介质 又根据对称轴的水平和垂直细分为， VTI介质和HTI 介质（图 1） 。 图 1 VTI 和 HTI 介质模型 Fig.1 VTI and HTI media models 两介质都需要 5 个独立的弹性参数来描述，其 弹性矩阵分别为 11116613 11662223 313233 VTI 44 44 66 2000 2000 000 00000 00000 00000 CCCC CCCC CCC C C C C             （8） 111113 11223344 31334433 HTI 44 66 66 000 2000 2000 00000 00000 00000 CCC CCCC CCCC C C C C             （9） 分别将两式带入公式（7） ，可得 VTI 介质三维 波动方程 2222 116644 2222 2 2 1166661344 xxxx y z x uuuu CCC txyz u u CCCCCf x yx z            （10a） 2222 116644 2222 22 1166661344 yyyy xz y uuuu CCC txyz uu CCCCCf y xy z           （10b） 第 49 卷第 3 期 董载天，等基于 VTI 介质的隧道超前探测 qP 波正演模拟研究 685 2222 444433 2222 2 2 1344661344 zzzz y x z uuuu CCC txyz u u CCCCCf x yy z            （10c） 由于 VTI 介质与地下实际介质的响应特质较为 一致且研究最为广泛，所以笔者主要使用 VTI 介质 的弹性波波动方程，因此此处不再给出 HTI 介质的 波动方程。 在利用此方程进行TI介质的弹性波正演模拟中， 波动方程的包含三个分量，占用计算机大量内存， 同 时由于纵波速度快于横波， 而网格的空间步长和采样 时间步长主要受最波的小速度影响，为稳定计算， 只 能减小空间和采样的间隔， 使得计算的时间点数与空 间点数大幅提高，进一步减慢了计算速度，不利于计 算机模拟。且目前实际勘探仍以纵波为主，多分量数 据中进行横波分离仍然是个难点问题。 因此笔者借用Alkhalifah提出的声学假设条件， 直接将横波速度设为零，从而在复杂的 VTI 介质 弹性波方程中分裂出 qP 波方程 ［14］ 。该方程不产生 横波，虽然现实中上不可实现，但在运动学上能对 弹性波方程进行有效地近似，便于研究 TI 介质中 qP 波的传播规律。 根据Alkhalifah的假设导出的三维 VTI 介质中 的 qP 波方程为 444 22 42222 444 22 222222 12 2 v v FFF vv ttztx FFF v v tyzxzy              （11） 式中F 为波场函数； t 为时间，s; x、y、z 为三维空间坐标; v 、vv qP 波的动校正速度和垂向速度，m/s; η 各向异性参数。 已知、为 Thomsen 参数 ［15］ 。 由于隧道超前探测观测方向为掌子面正前方， 而常规地震勘探观测方向为地面垂向（如图 2） ，所 以笔者后面推导VTI介质中的波动方程时只保留xy 分量。 因此，由式（11）可得二维 VTI 介质中 qP 波 方程 （12） 图 2 隧道超前探测与地震勘探观测系统差异 Fig.2 Differences between tunnel advance detection and seismic exploration observation 直接求解上式需要消耗大量内存和时间，所以 利用，可将上式转化为 222 222 2 PPP txy     （13） 其中为应力。 借鉴 Hestholm 推导二维 VTI 介质 中 qP 波方程的一阶应力-速度方程形式的方法，可 得 11 , y x v vPP txty       （14a） , y x v v xy        （14b） 2 P t     （14c） 式中 弹性体密度，kg/m; vx、v y x、y 方向速度，m/s; 、 计算过程中的中间变量。 2 交错网格 在实际计算过程中必须对介质模型进行离散， 将介质模型网格化。实际应用过程中，交错网格相 比规则网格有更好的差分算子局部特性，更高的精 度和更好的收敛性。所以笔者采用交错网格（如图 3）将放置在网格节点。 图 3 交错网格划分 Fig.3 Staggered meshing 利用 Taylor 级数展开法，用网格节点上数值的 差商近似代替前式中的导数。 推导出了时间域 2 阶、 空间域 2N阶精度的交错网格离散差分格式。 444 2 42222 44 2222 1 2 2 FFF v ttxty FF txty             686 当 代 化 工 2020 年3 月 11 22 ,,1 , 2 1 , 111 ,,, 222 1 , 1 C [ 1 ]C [] N kk Nk i ji jn in j n i j N kNkk n in ji jni jn n i j t UUP x t PSS y                （15a） 11 22 11111 ,,, 22222 1 11 , 22 111 1,,,1 222 1 11 , 22 1 C [ 1 ]C [] N kk Nk n ijiji n j n ij N kNkk n i njij nij n n ij t VVS x t SQQ y                    （15b） 11 1 22 1111,1, ,, 22 1 121111 ,, 2222 1 C [ ]C [] N kk kkN ni n ji nj ijij n N Nkk n ijnijn n t PPCUU x t CVV y                 （15c） 11 1 22 1111,1, ,, 22 1 121111 ,, 2222 1 C [ ]C [] N kk kkN ni n ji nj ijij n N Nkk n ijnijn n t QQCUU x t CVV y                 （15d） 11 1 22 1144,,1 ,, 22 1 441111 ,, 2222 1 C [ ]C [] N kk kkN ni j ni j n i ji j n N Nkk n in jin j n t SSCUU x t CVV y                 （15e） 式中∆ 空间间隔，m; ∆ 时间间隔，s; i、j、k 空间和时间节点; U、V 速度分量 vx和 vy的离散值，m/s; P、Q、R 应力和的离散值。 3 有限差分模拟中的重要问题 3.1 PML 边界条件 针对隧道超前探测的弹性波数值模拟而设计的 数据模型，仅考虑深埋隧道时，隧道空间为自由边 界，计算区域外部边界则都是人工截断边界。人工 边界会使向外传播的弹性波完全反射回去，从而使 正演模拟计算结果与实际情况大不符。为了尽量降 低人工截断边界所产生的影响，引入 Berenger 针对 电磁波传播提出了完美匹配层（Perfect Matched Layer，简称 PML）吸收边界条件，削弱人工截断 边界的反射，改善正演模拟效果。 根据 PML 的推导思路，可给出式（14）的吸收 边界方程 1 d x x vP x v tx     （16a） 1 y y v P d x v ty      （16b） , y x v v xy        （16c） 2 x x P d x P t      （16d） 2 y y P d y P t      （16e） 式中d（x） 、d（y） 阻尼因子； R 反射系数； 完美匹配层吸收边界的厚度，m; h 完美匹配层吸收边界的深度，m; v 匹配层中的弹性波波速，m/s。 3.2 稳定性 波动方程的差分格式按照时间进行递推，前一 时刻的波函数数值解的舍入误差必然会对下一时刻 的波场产生影响。这就需要对误差传播情况进行估 计，使其不会随着时间推移而快速增长，以致影响 数值模拟的效果进而破坏计算，使模拟无法继续进 行。笔者此处借鉴 Virieux 在 1986 年提出的推导方 法，得到二维 VTI 介质中二阶时间差分精度、2N 阶空间差分精度的稳定条件 2 22 11 max , 1tv x y xy   （17） 当中表示模型中最大速度的平方。当模型横纵 空间步长一致时，亦可将方程改写为 2 1 max , 2 t v x y x    （18） 3.3 震源加载 此部分还有一个重要问题就是加载震源。本文 为保证成像的分辨率，选择使用零相位雷克子波作 为震源。其时间域表达形式为 222 M0 222 M0 [12π ] exp[ π ] s tftt ftt   （19） 式中t0 初始时刻，ms; f 主频, Hz。 4 数值模拟实例 建立隧道空间的模型（如图 4） ，正演模型的尺 寸为 300 m100 m， 模型四周使用 PML 吸收边界。 图 4 隧道数值模型 Fig.4 Tunnel numerical model 第 49 卷第 3 期 董载天，等基于 VTI 介质的隧道超前探测 qP 波正演模拟研究 687 隧道尺寸为 60 m8 m，空腔中充填空气，波 速 340 m/s; 围岩波速 3 600 m/s。 采样间隔为 0.2 ms， 采样时长为 0.1 s。 观测方式如所示，掌子面中心单震源，左右两 侧墙壁各布置 4 个检波器，间隔 2 m，水平最小炮 间距 10 m。震源与检波器皆置于墙内，埋深 2 m。 图 5 观测系统示意图 Fig.5 Schematic diagram of the observation system 通过正演模拟可得到地震响应记录（图 6） ，为 共炮点道集（CSP） 。可以看出直达波能量很强，无 边界反射和伪 S 波干扰。根据费马原理计算的初至 时间在 4 ms 左右，基本与记录吻合。 图 6 无异常体 CSP 道集 Fig.6 No abnormal body CSP gather 4.1 不良地质体模型及其地震响应 隧道施工过程中会遇到很多危害施工安全的不 良地质体，常见的有断层破碎带、软弱夹层等。 由于二维平面不存在倾角，所以后续模型中需要关 注不良地质体的走向。 4.1.1 断层模型及响应 断层破碎带是由于断层上下盘搓动引起断面附 近的岩体发生破碎，形成破碎带。其稳定性较差， 容易引发隧道塌方；在地下水富集区域，破碎带变 成了一条良好水通道，极易引发突水突泥灾害。由 于破碎带与断层的伴生特性且破碎带本身理化性质 过于复杂，笔者此处将断层破碎带的模型设置为断 面两侧阻抗不同， 破碎带并不在模型当中直接体现。 数值模型构建时，首先选择走向与隧道方向垂直， 距掌子面80和120 m的断层， 断层右盘波速为4 000 m/s（图 7） 。 对比两模型的 CSP 道集（如图 8）可以看出， 当断面与掌子面距离减小时，阻抗界面的反射波逐 渐靠近强振幅的直达波，初至有可能被遮蔽，此时 仅靠时域记录是难以判断前方断面距离的；当两者 间距拉大，反射波整体向后时移，初至明显但振幅 减弱， 现场施工噪音较强时， 会难以见到完整波形。 图 7 断层模型（左图掌子面与断面间距为 80 m，右图为 120 m） Fig.7 Fault model （（left hand spacing and section spacing is 80 m, right picture is 120 m）） 图 8 断层模型地震响应 （左图掌子面与断面间距为 80 m， 右图为 120 m） Fig.8 Seismic response of the fault model （（the distance between the face and the section of the left is 80 m, and the picture on the right is 120 m）） 4.1.2 软弱夹层及响应 软弱夹层成因复杂且尚无统一的分类原则， 此处 只针对沉积型软弱夹层进行分析。 此种软弱夹层往往 形成与河流沉积环境，层中粉砂岩、泥岩含量偏高， 强度低质地较软弱 ［17］ 。当构造作用使夹层被剪切、 破碎时， 无论区域是否富水， 都会产生低速低密度的 软弱夹层。 软弱夹层的数值模型参考断层模型，将夹 层设置在掌子面前方 80 m 处且走向与隧道方向垂直， 此时一次反射波不会被直达波遮蔽。分别构建 10、 20 m 厚的夹层，波速为 800 m/s（如图 9） 。 此时获得 CSP 道集（如图 10） 。由于此时夹层 与围岩阻抗差异过大，且检波器内置在墙壁中，导 致阻抗界面反射波振幅较强。对于薄夹层来说，子 波在夹层中旅行时较短，夹层右侧依旧为高阻抗差 界面， 所以 10 m 薄夹层的地震反射记录上出现多次 波。 而层厚为 20 m 的夹层， 子波旅行时被拉长此时 688 当 代 化 工 2020年3月 反射波没有在记录时间内返回到检波器，故道集上 无明显反射波。 图 9 不同厚度软弱夹层模型 （左图夹层厚为 10 m，右图为 20 m） Fig.9 Weak interlayer model with different thicknesses （（the left layer is 10 m thick and the right image is 20 m）） 图 10 不同厚度软弱夹层地震响应 （左图为 10 m 夹层地震 响应，右图为 20 m 夹层地震响应） Fig.10 Seismic response of weak interlayers with different thicknesses （（the left picture shows the seismic response of 10 m interlayer, and the right picture shows the seismic response of 20 m interlayer）） 将软弱夹层与掌子面之间距离拉大至 120 m （模型如图 11） 。 图 11 断面与掌子面距离为 120 m 的软弱夹层模型（左图 夹层厚为 10 m，右图夹层厚为 20 m） Fig.11 Weak interlayer model with a distance of 120 m from the section of the face （（the thickness of the interlayer is 10 m in the left and 20 m in the right））. 从两模型 CSP 道集（如图 12）可以看出，此时 两种模型的地震响应颇为一致，显然夹层右界面的 反射波由于旅行时较长没有在记录时长内返回到检 波器位置。 构建一个右盘为低速岩体的断层（如图 13） ， 此时获得其 CSP 道集（如图 14） ，其与掌子面距离 为 120 m 的软弱夹层的地震响应基本一致。可见固 定采样时长有时无法分辨的不同的不良地质体，此 时需要增加采样时长或者利用其他地质勘探资料来 确定不良地质体形态。 图 12 与掌子面距离为 120 m 的软弱夹层地震响应（左图 夹层厚为 10 m，右图夹层厚为 20 m） Fig.12 Seismic response of the weak interlayer with a distance of 120 m from the face of the hand （（the thickness of the interlayer is 10 m on the left and 20 m on the right））. 图 13 断层模型 （掌子面与断面间距 120 m， 右盘速度 800 m/s） Fig.13 Fault model （（distance between face and section of the face is 120 m, right disk speed is 800 m/s）） 图 14 断层模型地震响应 （断层右盘波速 800 m/s， 距掌子 面 120 m） Fig.14 Seismic response of the fault model （（wave velocity of the right disk of the fault is 800 m/s, 120 m from the face）） （下转第 712 页） 712 当 代 化 工 2020年3月 度 285 ℃， 压力 2.28 MPa 的条件下， 通过脱硫反应 器后，能使混合干气中硫含量从 21 μg/g 左右下降 至 0.03 μg/g 左右，床层温升稳定，满足转化进料 需求。 （4） 本次催化剂活化温度未达到指定温度， 部 分催化剂未深度活化，从初期的运行情况及产品分 析数据来看，产品质量均满足设计要求，后续生产 中需密切关注烯烃、硫含量，但从目前使用情况来 看，催化剂活性良好，产品质量满足要求。 （5） 在本周期运行中， 如果产品各项数据满足 要求，则在今后催化剂开工活化过程中，可将 JT-4 和 JT-1G 载硫型催化剂的活化温度控制在 260 ℃ 左右。 （6） 利用原料中烯烃含量提高催化剂活化时床 层温度有一定效果，为同类生产装置遇到类似问题 提供了参考。 参考文献 ［1］刘海涛，杨桂香，李晓明. JT-1G 加氢催化剂在焦化干气制氢装置 上的应用[J]. 炼油设计，1997，27（6）60-62． ［2］郝树仁，董世达．烃类转化制氢工艺技术[M]. 北京石油工业出 版社，2009. ［3］杜彩霞，周晓奇，霍尚义. JT-1G 型和 JT-4 型炼厂干气加氢精制 催化剂的开发及工业应用[J]. 工业催化，2000，8（2）54-60． ［4］汪加海.预硫化态催化剂在柴油加氢改质装置上的应用[J].当代化 工，2014，43（3）460-463． ［5］翟玉娟，郭英宽，孙志明，等. JT-1G 催化剂在制氢上的工业应用 [J]. 炼油技术与工程，2018，48（9）41-44． ［6］ 龚建友. JT-4 和 JT-1G 型加氢催化剂的预硫化[J]. 工业催化， 2003， 3（11）13-15． ［7］刘凯，郭晓雷，李云鹏.论渣油加氢装置原料油性质对催化剂的影 响[J].化工管理，2014（10）183. ［8］李立权．加氢催化剂硫化技术及影响硫化的因素[J]. 炼油技术与工 程，2007（3）55-62. ［9］ 冯续， 崔芳. 有机硫加氢 （HDS） 催化剂的预硫化[J]. 大氮肥，2003 （1）38-42. （上接第 688 页） 5 结 论 隧道超前探测尚处于发展阶段，由于其探测空 间的复杂以及其与传统地震勘探的差异，使得隧道 超前探测在实际应用过程中存在一些问题。笔者从 弹性体基本理论出发，逐步推导出适用于隧道超前 探测二维正演模拟公式，并模拟了断层、软弱夹层 的地震响应。从中得到一些结论 （1）利用三维 VTI 介质波动方程，根据隧道超 前探测的观测方式，推出了适用于隧道勘探正演模 拟的二维 VTI 介质 qP 方程，后导出其 1 阶应力-速 度方程， 进行时间 2 阶空间 2N阶的有限差分正演模 拟。其计算迅速，规避了伪 S 波干扰和计算不稳定 性，便于快速得到不同地质模型的地震响应特征。 （2）利用 PML 吸收边界，处理了隧道后方人 工边界引起的反射。采用交错网格和，使方程能较 好地模拟断层和软弱夹层的地震响应，减少频散。 （3）二维平面模拟中，采样时长与断层、软弱 夹层的物性参数组合不当，会导致两种不同的不良 地质体产生相同的地震响应。在数值模拟时，可改 变模型参数，使其显现不同的地震响应。而在现场 实际情况来说，需要现场工程师结合其他地质勘探 资料综合分析。 参考文献 ［1］柳杨春. HSP 地质超前预报技术及其应用[J]. 西部探矿工程, 1997 （05） 40-42. ［2］刘志刚, 刘秀峰. TSP（隧道地震勘探）在隧道隧洞超前预报中的 应用与发展[J]. 岩石力学与工程学报, 2003, 22 （8） 1399-1402. ［3］曾昭璜. 隧道地震反射法超前预报[J]. 地球物理学报, 1994 （2）． ［4］ 钟世航. 陆地声纳法的原理及其在铁路地质勘测和隧道施工中的应 用[J]. 中国铁道科学, 1995 （4） 48-55. ［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］Otto R, Button E, Bretterebner H, Schwab P.The application of TRT-trae reflection tomography-at the Unterwald Tunnel[J]. Geophysics, 2002, 20 （2） 51-56. ［7］Yamamoto T, Shirasagi S, Yokota Y, Koizumi Y. Imaging geological conditions ahead of a tunnel face using three-dimensional seismic reflector tracing system[J]. International Journal of the JCRM, 2011, 6 （1） 23-31. ［8］赵永贵. TST 隧道超前预报技术及应用[C].中国地球物理学会第 22 届年会论文集,2006. ［9］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. ［10］宋先海, 顾汉明, 肖柏勋. 我国隧道地质超前预报技术述评[J]. 地 球物理学进展, 2006 （02） 605-613. ［11］张平松,吴健生.中国隧道及井巷地震波法超前探测技术研究分析 [J]. 地球科学进展, 2006 （10） 1033-1038. ［12］张学民. 岩石材料各向异性特征及其对隧道围岩稳定性影响研究 [D]. 中南大学, 2007. ［13］程玖兵, 康玮, 王腾飞. 各向异性介质 qP 波传播描述Ⅰ 伪纯模 式波动方程[J]. 地球物理学报, 2013, 56 （10） 3474-3486. ［14］Alkhalifah, Tariq. An acoustic wave equation for anisotropic media[J]. Geophysics, 2000, 65（4） 44-48. ［15］Thomsen, Leon. Weak Elastic Anisotropy[J]. Geophysics, 1986, 511954-1966. ［16］Clayton R, Engquist B. Absorbing boundary conditions for acoustic and elastic wave equations[J]. Bull. seis