北京大学学报(自然科学版) 第62 卷 第4 期 2026 年7 月
Acta Scientiarum Naturalium Universitatis Pekinensis, Vol. 62, No. 4 (July 2026)
doi: 10.13209/j.0479-8023.2026.031
中国石油天然气集团有限公司-北京大学基础研究项目“智能物探的超前基础研究与共享框架”资助
收稿日期: 2026-02-01;
修回日期: 2026-04-24
A Reverse Time Migration Imaging Method for Tunnel Ahead Detection Based on Time-Delay Stacking
在煤炭、水利和地铁等地下工程探测场景中,需要提前探测和预知掘进隧道前方的断层、溶洞和破碎带等异常体, 这类工作即为超前探测[1]。对于不同场景的超前探测, 如果能做到高精度成像, 将助力安全生产, 经济意义巨大。
地震偏移方法通过叠加得到一幅能真实地反映地下构造形态的图像[2]。该方法在油气勘探中有大量应用研究, 发展较为成熟[3-12]。同时, 该方法在浅地表结构超前探测领域也取得发展和应用。曾昭璜[13]针对隧道施工掌子面前方简单地质结构探测任务, 提出利用基于地震反射波法得到的负视速度同相轴来识别前方的界面。何振起等[14]将负视速度法和地震横排列法较好地应用到朔黄线长梁山隧道和福州飞鸾岭公路隧道工程中。李永鸿等[15]利用反射波法, 对沪蓉西高速公路的扁担垭隧道工程进行超前地质预测, 施工结果证明预测结果基本上正确。常旭等[16]将逆时偏移(reverse time migration, RTM)方法应用到隧道工程的超前预测中, 并给出一种隧道观测方案来适应逆时偏移成像, 通过数值计算结果验证 RTM 对超前探测的适应性。沈鸿雁等[17]根据地震学几何原理, 研究隧道条件下的时距曲线, 实现一种基于绕射叠加的深度域偏移成像方法, 提高了反射波法的超前探测成像精度。Wang等[18]提出一种基于非结构网格的超前探测RTM 方法, 通过采用隧道双边观测系统, 进一步改善倾斜界面的成像质量。此外, 一些研究者探讨了复杂地质条件下 RTM 的成像问题[19-22]。
尽管超前探测的观测条件和探测目标与地面勘探存在显著的差异, 但二者在偏移成像的核心思路上具有共通性。在地面勘探场景中, 对于不同位置激发的地震记录, 可以通过叠后偏移和叠前偏移方法处理[23]。其中, 叠后偏移通常是先进行共中心点道集(不同的炮点和检波点在地面的中心点位置相同)叠加, 然后对叠加后的地震记录进行偏移, 需要满足水平层状介质和零偏移距自激自收的假设, 具有较高的计算效率; 叠前偏移则直接对原始的未叠加共炮点道集或共检波点道集进行偏移处理, 再将每一次偏移结果叠加起来, 适用于任意复杂地质条件, 但具有更高的计算成本。
相较于基于地表观测的向下探测模式, 超前探测的几何关系发生变化, 相当于将探测目标旋转90°。因此, 常用的基于水平层状介质和零偏移距自激自收假设的叠后偏移方法在超前探测中不再适用。目前, 超前探测领域普遍采用叠前偏移方法。然而, 现有方法仅利用不同激发点位置对应的波场信息, 未能充分利用各震源的激发时刻这一参数。作为波动方程初值问题的重要组成部分, 震源激发时刻对实现高精度波场延拓以及成像归位具有关键作用。如何在偏移框架下有效地引入并利用这一信息, 是进一步提升超前探测成像精度的关键。
本文提出一种基于时延叠加的超前探测逆时偏移成像方法, 通过对不同位置震源的激发时刻设定不同的数学组合, 实现地震波沿特定方向或在目标点相位保持不变, 从而增强振幅, 并结合逆时偏移方法, 提高超前探测成像的信噪比。
RTM 方法的核心是利用地震波信息实现地下构造的成像。当密度连续变化且无明显突变时, 密度梯度项对波场传播的影响较弱, 变密度波动方程可近似为常密度形式。本文以二维常密度声波介质为例进行讨论, 地震波场满足以下标量波动方程:
其中, x 和 z 为空间坐标; t 是地震波传播的时间; v(x,z)是纵波速度, 只与空间坐标相关; p(x, z, t)表示地震波场, 是与时间和空间相关的标量; fs(x-xs,z-zs,t-τs)为震源时间函数, 地震勘探中常用Ricker子波作为震源; xs表示震源的 x 轴坐标, zs表示震源的 z轴坐标, τs表示震源的激发时刻。
地震波场的数值模拟是逆时偏移成像的关键技术。在正向传播过程中, 通常采用有限差分或有限元等数值离散方法求解波动方程, 从而获得每一时刻的正向波场。反向波场的构建则依赖于地震波传播过程的确定性。在弹性波动力学框架下, 若不考虑介质的非弹性衰减, 波场的演化在时间上是对称且可逆的。基于这一物理原理, 可以将预处理后的观测数据(如切除直达波后的共炮点道集)作为边界条件, 通过逆时延拓(即反向求解波动方程)重构波场在地下的传播过程。最终, 将适当的成像条件应用于正向波场与反向波场, 即可实现地下界面的准确归位。
以互相关成像条件为例, 偏移公式为
其中, j 表示时间采样的序号, m 为地震记录时长与时间采样间隔的比值, pS 表示震源激发的正向波场,pR 表示检波器记录切除直达波后逆时延拓产生的反向波场。
当考虑激发时刻信息时, 在二维平面上, 同一条测线上的 n 个震源在不同时刻激发产生的地震波场可以表示为
其中,i表示震源的序号, pi表示第 i 个震源激发产生的地震波场, 是与时间和空间相关的标量; τi表示第 i 个震源的激发时刻, 第 i 个震源与第 1 个震源激发时刻的差值(简称时延), 表示为ΔTi=τi-τ1。
对于不同时刻激发的地震记录, 如何利用激发时刻的信息进一步增强超前探测的成像精度是本文关心的问题。我们通过时延设计, 在特定方向或目标点实现地震波场的相位不变和振幅叠加。
如图1 所示, n 个震源沿隧道(x 轴方向)排列, 第i 个震源与第 1 个震源的距离用
表示。当我们关注某一方向的地震波时, 目标聚焦方向由角度α 定义。
图1 面向目标方向的地震波叠加示意图
Fig.1 Schematic diagram of seismic wave superposition towards the target direction
对于目标方向 α, 第 i 个震源与第 1 个震源在该方向上的波程差由二者之间的距离 Li 在目标方向上的投影决定:
式(4)的几何意义为将不同震源的水平距离投影到目标方向上, 得到其在该方向上的路径差。当 α= 90°时, sinα=1, 波程差变为 Li, 对应叠加最强能量平面波的波前与震源测线垂直的情形(即隧道正前方探测场景)。假设地震波在震源所在介质中以纵波速度 V1 传播, 波程差对应的传播时间差就是不同震源所需设定的时延:
通过设定该时延, 可使所有震源在目标方向 α上的地震波场实现同相位到达, 从而形成能量聚焦, 其物理本质是对震源位置差异导致的波程差进行时延设定, 确保目标方向上各震源波场相位保持一致且振幅叠加。
如图2 所示, 当关注某一目标点 I0 时, 我们也可以通过时延设定, 实现对单一目标点的地震波相位不变、能量叠加到最强, 此时ΔTi应表示为
图2 面向目标点的地震波叠加示意图
Fig.2 Schematic diagram of seismic wave superposition towards the target point
表示第 i 个震源与目标点 I0 的距离, 
表示第1 个震源与目标点 I0 的距离。ΔTi一般为正值, 代表第 i 个震源相较于第 1 个震源为延后激发; 如果结果为负值, 表明第 i 个震源距离成像目标点更远。为简化表达和方便计算, 本文选择距离成像目标点最远的震源作为第 1 个震源, 保证时延都是正值。
基于时延设定公式(式(5)和(6)), 我们可以从物理叠加和数学叠加两个角度来实现这种效果。对于物理叠加, 可以设计对应的可控组合震源, 根据特定方向或目标点, 通过控制不同震源的激发时刻来实现地震波场叠加。相比单个震源, 物理叠加的地震波信号在特定方向或目标点显著增强, 而噪声环境不变。因此, 物理叠加可以同时压制随机噪声和相干噪声。然而, 物理叠加对野外施工布置和震源激发时刻的精度要求较高, 如果精度达不到要求,则效果大打折扣。并且, 每增加一个目标, 就需要多做一次施工, 经济成本非常高。对于数学叠加,可以先在不同的位置独立激发地震波, 然后将每一次独立激发产生的地震记录转换为同一时刻激发的地震记录(本文中转换成 0 时刻), 根据特定方向或目标点来计算不同震源位置的地震记录所需设定的时延, 对时延设定后的地震记录进行叠加, 同样可以压制随机噪声(当不考虑波场衰减效应时, 叠加N 次, 特定方向或目标点信号的振幅增强 N 倍, 随机噪声则增强N倍), 以后每次增加, 均不需要额外施工, 只需再次计算所需设定的时延, 极大地降低施工成本。尽管时延设定既可通过控制各震源的实际激发时刻进行物理实现, 也可通过对独立激发的各炮记录施加相应时移进行数学实现, 但考虑到经济性和工程实施的便利性, 本文主要研究后一种实现方式。
我们利用一个垂直断层模型来验证本文数学叠加方法对特定方向地震波振幅的增强效果。图3 所示模型的横向距离为 800 m, 纵向深度为 400 m。计算时, 设置离散化网格大小为 801×401, 网格间距为 1 m。隧道掌子面与隧道底板的交点设置在模型正中间, 即坐标(400 m, 200 m)处, 其前方 100 m 处设计一个长 400 m, 宽 1 m 的垂直断层。模型的纵波速度为 2000 m/s, 断层的纵波速度为 3000 m/s。第一个震源位置的坐标为(200 m, 200 m), 共有 11 个震源, 每个震源在隧道底板均匀地布设, 间隔为 20 m, 因此最后一个震源位置的坐标为(400 m, 200 m)。在 整个隧道底板布设检波器, 共有 400 个, 间隔为 1 m。若想实现地震波在正前方传播时相位不变且振幅增强, 则式(5)中 α 应取 90°, 则相邻两炮之间的时延为 10 ms。
图3 长 400 m, 宽 1 m 的垂直断层模型及观测系统示意图
Fig.3 Schematic diagram of vertical fault model with the length of 400 m and the width of 1 m and the observation system
在正演过程中, 采用 Ricker 子波作为震源。子波的主频取 100 Hz, 时间采样间隔为 0.2 ms, 地震记录时长为 0.6 s。图4 为利用第一个震源(200 m, 200 m)正演 50, 100, 150, 200, 250 和 300 ms 时刻的波场快照, 图5 为利用式(5)中 α=90°时的数学叠加方法正演 50, 100, 150, 200, 250 和 300 ms 时刻的波场快照。与图4 相比, 图5 中 50 和 100 ms 时刻的地震波在正前方传播时保持相位不变, 振幅不断增强; 当第 11 个震源激发之后(150, 200, 250 和 300 ms 时刻), 图5 中地震波振幅的提升保持稳定。
图4 图3 中速度模型利用第一个震源(200 m, 200 m)正演不同时刻的波场快照
Fig.4 Wavefield snapshots at different times using the first source (200 m, 200 m) for the velocity model in Fig. 3
图5 图3 中速度模型利用式(5)中 α = 90°时的数学叠加方法正演不同时刻的波场快照
Fig.5 Wavefield snapshots at different times using the mathematical implementation of time-delay stacking with α=90° in Eq. (5) for the velocity model in Fig. 3
为验证本文数学叠加方法在实际噪声环境下的有效性, 我们计算 11 个震源在 0 时刻独立激发产生的地震记录, 并切除直达波, 结果如图6 所示。可以发现, 随着震源不断向右移动, 反射波的走时不断向 0 时刻靠近, 证明异常体来自震源右侧。为了更符合实际环境, 在图6 中加入振幅最大值为 10-4的高斯白噪声, 得到图7 所示地震记录。图7 中看不到连续的同相轴信息, 反射波信号完全被随机噪声淹没。为定量地评价图7 中地震资料的质量, 本文依据 GB/T 33684—2017《地震勘探资料解释技术规程》[24], 将信噪比作为核心评价指标, 并采用有效信号与噪声能量比的对数形式[25]来表征:
图6 图3 中速度模型利用 11 个震源在 t = 0 时刻独立激发产生的地震记录切除直达波后的结果
Fig.6 Seismic records generated by independent excitation of 11 sources at t=0 after direct wave removal for the velocity model in Fig. 3
图7 图6 所示结果中加入高斯白噪声后的地震记录
Fig.7 Seismic records after adding Gaussian white noise to the results shown in Fig. 6
其中, ES代表信号窗口内的平均能量, EN代表噪声窗口内的平均能量, SNR 代表信噪比(dB)。
参照地震勘探行业通用经验, 信噪比大于 20 dB 时, 信号品质优良; 信噪比为 9.5～20 dB 时, 主要反射界面可连续追踪; 信噪比为 0～9.5 dB 时, 仅能识别强反射界面; 信噪比小于 0 dB 时, 信号无法识别[25]。经过计算, 图7 中地震记录的信噪比为-25.98 dB, 远小于 0, 说明其信噪比极差, 可以代表实际的强噪声环境。我们基于图7 中地震记录, 利用式(5)中 α=90°时的数学叠加方法, 得到图8 所示地震记录, 其中可以看到明显的连续同相轴。
图8 图7 所示结果利用式(5)中 α = 90°时的数学叠加方法得到的地震记录
Fig.8 Seismic records using the mathematical implementation of time-delay stacking with α=90° in Eq. (5) to the results shown in Fig. 7
我们将背景速度 2000 m/s 作为 RTM 的偏移速度, 利用图7 中每炮的地震记录, 分别执行一次RTM 计算, 再对 11 个偏移结果进行叠加(叠前 RTM方法), 得到图9(a)所示偏移结果; 然后, 利用图8中地震记录, 只执行一次 RTM 计算(本文数学叠加方法), 得到图9(b)所示偏移结果。我们在 z=201 m 处抽取两个偏移结果, 得到波形对比图(图9(c))。选定 x=491～511 m 为信号时窗, 其余为噪声时窗, 利用式(7), 得到叠前 RTM 方法的信噪比为 20.27 dB, 本文数学叠加方法的信噪比为 22.75 dB。上述结果表明, 在关注特定方向时, 相较于常规的叠前偏移方法, 本文数学叠加方法可以进一步提高偏移结果的信噪比, 并在只关注一次特定方向的情况下, 仅计算一次偏移, 节省了计算成本。
图9 图3 中速度模型不同方法的偏移结果和波形对比
Fig.9 Comparison of migration and waveform results by different methods for the velocity model in Fig. 3
我们通过改变断层的宽度来验证本文数学叠加方法在不同模型中的表现。图10(a)中垂直断层的宽度为 5 m, 其余参数与图3 一致。图10(b)和(c)显示, 偏移结果仍然只有一个强能量圆弧, 这是由于垂直断层 5 m 的宽度小于 1/4 个地震波长(30 m)。在z=201 m 处抽取的两个偏移结果的波形(图10(d))进行对比, 选定 x=491～511 m 为信号时窗, 其余为噪声时窗, 利用式(7), 得到叠前 RTM 方法的信噪比为 31.00 dB, 本文数学叠加方法的信噪比为 34.36 dB, 证明本文方法能够增强偏移结果的信噪比。
图10 长 400 m, 宽 5 m 的垂直断层速度模型及利用不同方法的偏移结果和波形对比
Fig.10 A 400 m-long by 5 m-wide vertical fault velocity model, together with its migration results and waveform comparisons obtained by different methods
图11(a)中垂直断层的宽度为 20 m, 图11(b)和(c)显示偏移结果有两个强能量圆弧, 分别对应垂直断层和背景介质的左右界面。在 z=201 m 处抽取两个偏移结果的波形(图11(d))进行对比, 选定 x=491～531 m 为信号时窗, 其余为噪声时窗, 利用式(7), 得到叠前 RTM 方法的信噪比为 27.09 dB, 本文数学叠加方法的信噪比为 30.17 dB。上述结果表明, 对于不同宽度的垂直断层模型, 相较于常规的叠前偏移方法, 本文的数学叠加方法对偏移结果的信噪比均有稳定的提升。
图11 长 400 m, 宽 20 m 的垂直断层速度模型及利用不同方法的偏移结果和波形对比
Fig.11 A 400 m-long by 20 m-wide vertical fault velocity model, together with its migration results and waveform comparisons obtained by different methods
接下来, 我们利用绕射波成像来验证本文数学叠加方法对目标点成像的可行性。如图12 所示, 我们设计一个具有两个绕射点的垂直断层模型, 两个绕射点的坐标分别为(500 m, 120 m)和(500 m, 20 m), 其余参数与图3 一致。断层中点的坐标为(500 m, 70 m)。通过计算可知, 断层中点和第一个震源坐标(200 m, 200 m)连线与隧道底板的夹角为 23.4°, 若想使地震波沿着这个连线的方向传播时相位不变且振幅增强, 则式(5)中 α 应取 66.6°。
图12 具有两个绕射点的垂直断层模型及观测系统
Fig.12 Vertical fault model with two diffraction points and the observation system
分别利用第一个震源(200 m, 200 m)、式(5)中 α =66.6°时的数学叠加方法以及式(6)中目标点为(500 m, 120 m)和(500 m, 20 m)时的数学叠加方法, 正演得到 150 ms 和 200 ms 时刻的波场快照(图13)。可以发现, 在 150 ms 时刻, 3 种数学叠加方法已完成 11 个震源的激发, 并且 4 种波场存在显著的差异; 在200 ms 时刻, 地震波在遇到垂直断层界面时产生反射波(图13 中橙色箭头), 遇到断层端点时产生绕射波(图13 中红色和黄色箭头), 虽然反射波的能量显著强于绕射波, 但是隧道底板的检波器只能接收到绕射波信号, 接收不到能量更大的反射波信号。
图13 利用图12 中速度模型正演得到的波场快照
Fig.13 Wavefield snapshots using the velocity model in Fig. 12
为了验证本文数学叠加方法在实际噪声环境下对绕射波成像的有效性, 我们首先计算 11 个震源在t = 0 时刻独立激发产生的地震记录, 并切除直达波,结果如图14 所示。可以发现, 断层下端点产生的绕射波能量明显大于上端点产生的绕射波能量, 这主要是由距离衰减导致的。
图14 图12 中速度模型利用 11 个震源在 t= 0 时刻独立激发产生的地震记录切除直达波后的结果
Fig.14 Seismic records generated by 11 sources excited independently at t=0 with direct waves removed for the velocity model in Fig. 12
为了更符合实际环境, 在图14 的结果中加入振幅最大值为 10-5 的高斯白噪声, 得到图15 所示地震记录。图15 中看不到连续的同相轴, 绕射波信号完全被随机噪声淹没。利用式(7)得到的计算结果中,图15 中地震记录的信噪比为-29.17 dB, 远小于 0,可代表关注绕射波信号情况下的强噪声环境。基于图15 中地震记录, 分别利用式(5)中 α = 66.6°时的数学叠加方法以及式(6)中目标点为(500 m, 120 m)和(500 m, 20 m)时的数学叠加方法, 得到图16(b)～(d)所示地震记录。相较于第一炮加入随机噪声后的地震记录(图16(a)), 图16(b)和(c)中可看到明显的连续同相轴(红色箭头所示), 对应断层下端点产生的绕射波, 图16(d)中也能看到相对明显的同相轴(黄色箭头所示), 对应断层上端点产生的绕射波。
图15 图14 所示结果加入高斯白噪声后的地震记录
Fig.15 Seismic records with Gaussian white noise added to the results shown in Fig. 14
图16 图15 所示结果利用不同方法得到的地震记录
Fig.16 Seismic records using different methods to the results shown in Fig. 15
同样, 将背景速度 2000 m/s 作为 RTM 的偏移速度, 利用图15 所示每炮地震记录先分别执行一次RTM 计算, 然后对 11 个偏移结果进行叠加(叠前RTM 方法), 得到图17(a)所示偏移结果。然后, 分别利用图16(b)～(d)中地震记录, 各执行一次 RTM计算(本文数学叠加方法), 得到图17(b)～(d)所示偏移结果。可以发现, 4 组 RTM 偏移结果中红色和黄色箭头分别对应模型的两个绕射点。不同之处是, 图17(a)中红色和黄色箭头处圆弧能量相当, 而图17(b)和(c)中红色箭头处圆弧能量显著增强, 图17 (d)中黄色箭头处圆弧能量显著增强。上述结果验证了本文数学叠加方法对目标点成像的有效性。
图17 不同方法的偏移结果对比
Fig.17 Comparison of migration results by different methods
进一步, 我们在 z=120 m 处抽取图17(a)和(b)中两个偏移结果, 得到波形对比图(图18(a))。选定 x= 481～521 m 为信号时窗, 其余为噪声时窗, 利用式(7)得到叠前 RTM 方法的信噪比为 21.10 dB, 式(5)中 α=66.6°时数学叠加方法的信噪比为 21.66 dB。同样, 在 z=120 m 处抽取图17(a)和(c)中两个偏移结 果, 得到波形对比图(图18(b))。选定 x=481～521 m为信号时窗, 其余为噪声时窗, 利用式(7)得到叠前RTM 方法的信噪比为 21.10 dB, 式(6)中目标点为(500 m, 120 m)时数学叠加方法的信噪比为 22.14 dB。可以发现, 相比于式(5)中 α=66.6°时的数学叠加方法, 利用式(6)中目标点为(500 m, 120 m)时的数学叠加方法对断层下端点产生的绕射波能量叠加更强, 偏移结果的信噪比更高。最后, 在 z=20 m 处抽取图17(a)和(d)中两个偏移结果, 得到波形对比图(图18(c))。选定 x=481～521 m 为信号时窗, 其余为噪声时窗, 利用式(7)得到叠前 RTM 方法的信噪比为 14.52 dB, 式(6)中目标点为(500 m, 20 m)时数学叠加方法的信噪比为 16.72 dB。这些结果表明, 相较于常规的叠前偏移方法, 本文提出的数学叠加方法在关注某一目标点时, 也可以进一步提高偏移结果的信噪比。
图18 不同方法的地震波形对比
Fig.18 Comparison of seismic waveforms by different methods
本文提出一种基于时延叠加的超前探测逆时偏移成像方法, 并给出针对特定方向或目标点成像的震源激发时延设定公式。时延叠加既可通过控制各震源的实际激发时刻进行物理实现, 也可通过对独立炮记录施加相应时移进行数学实现, 本文重点推荐后一种方式。数值实验结果表明, 在关注特定方向或目标点时, 相比于常规的叠前偏移方法, 本文方法进一步提高了偏移结果的信噪比。
在多炮地震勘探场景下, 常规的叠前偏移方法与本文提出的“先数学叠加、后偏移”方法均依托多炮数据的叠加机制来提升信噪比。在理想模型假设下, 对于 N炮数据, 由于随机噪声在各炮记录间是非相干的, 其振幅叠加按N的规律增长。然而,由于信号具有相干性, 振幅叠加按 N 的规律增长。因此, 从理论上讲, 两种方法均能够实现信噪比按10 log10 N 的规律提升, 这是两者共同的物理基础。
常规叠前偏移方法采用“单炮独立偏移、照明补偿、成像域叠加”的流程, 由于不同震源位置对目标区域的照明强度存在差异, 信号在多轮独立照明补偿后发生能量损失及叠加失真, 难以达到理想的 N倍增益。同时, 随机噪声在每炮独立偏移过程中均被不同程度地放大, 最终在成像域累积, 导致背景噪声增强。本文数学叠加方法针对特定的方向或目标点, 先在数据域完成精准的时延对齐和相干叠加, 实现有效信号最大化的同相增强, 随机噪声则因非相干特性得到压制。随后, 仅进行单次统一的偏移与照明补偿, 避免多源照明不均带来的信号能量损失及噪声放大效应, 进一步提升成像结果的信噪比。
本文给出的时延设定公式主要由震源所在介质的纵波速度、震源位置、特定方向或目标点构成, 其中震源位置、特定方向或目标点均为已知量, 震源所在介质的纵波速度需要用速度建模方法来反演, 而反演结果存在一定的误差, 给最终的数学叠加效果带来一定的影响, 因此需要配合搜索类方法。在未来的研究中, 我们可以在计算资源允许的范围内, 对实际资料中速度的误差区间进行搜索, 在该区间内计算不同速度下的数学叠加地震记录和偏移结果, 根据偏移结果的信噪比来确定最理想的速度值。
[1] 彭苏萍, 孟召平. 矿井工程地质理论与实践. 北京:地质出版社, 2002
[2] 马在田. 地震成像技术: 有限差分法偏移. 北京:石油工业出版社, 1989
[3] Baysal E, Kosloff D, Sherwood J. Reverse time migration. Geophysics, 1983, 48(11): 1514-1524
[4] Wiggins J. Kirchhoff integral extrapolation and migration of nonplanar data. Geophysics, 1984, 49(8): 1239-1248
[5] Hill N. Gaussian beam migration. Geophysics, 1990,55(11): 1416-1428
[6] Alkhalifah T. Gaussian beam depth migration for anisotropic media. Geophysics, 1995, 60(5): 1474-1484
[7] Gray S. Gaussian beam migration of common-shot records. Geophysics, 2005, 70(4): S71-S77
[8] Dai W, Schuster G. Plane-wave least-squares reversetime migration. Geophysics, 2013, 78(4): S165-S177
[9] Zhou T, Hu W, Ning J. Broadband finite-difference Q-compensated reverse time migration algorithm for tilted transverse isotropic media. Geophysics, 2020, 85(5): S241-S253
[10] Zuo J, Niu F, Liu L, et al. 3D anisotropic P- and S-mode wavefields separation in 3D elastic reverse-time migration. Surveys in Geophysics, 2022, 43(3): 673-701
[11] Zhang D, Zhou T, Yang J, et al. Gaussian beam migration in exploration seismology: methods, advantages and implementation. Frontiers in Earth Science, 2025,13: 1-13
[12] Mao Q, Huang J, Shen Y. Efficient Q-compensated reverse time migration using a new decoupled viscoacoustic wave equation based on the generalized standard linear solid. Geophysics, 2025, 90(5): S145-S160
[13] 曾昭璜. 隧道地震反射法超前预报. 地球物理学报,1994, 37(2): 268-271
[14] 何振起, 李海, 梁彦忠. 利用地震反射法进行隧道施工地质超前预报. 铁道工程学报, 2000, 68(4): 81-85
[15] 李永鸿, 徐光黎, 杨银湖, 等. 地震反射波法技术及其在隧道超前地质预报中的应用研究. 岩土工程学报, 2005, 27(10): 1180-1184
[16] 常旭, 刘伊克, 桂志先. 反射地震零偏移距逆时偏移方法用于隧道超前预报. 地球物理学报, 2006, 49 (5): 1482-1488
[17] 沈鸿雁, 李庆春, 冯宏. 隧道反射地震超前探测偏移成像. 煤炭学报, 2009, 34(3): 298-304
[18] Wang J, Liu J, Cheng F, et al. Reverse time migration imaging of tunnels via the finite element method using an unstructured mesh. Applied Geophysics, 2020, 17 (2): 267-276
[19] 鲁光银, 熊瑛, 朱自强. 隧道反射波超前探测有限差分正演模拟与偏移处理. 中南大学学报(自然科学版), 2011, 42(1): 136-141
[20] Cheng F, Liu J, Qu N, et al. Two-dimensional pre-stack reverse time imaging based on tunnel space. Journal of Applied Geophysics, 2014, 104: 106-113
[21] 张焕钧. 基于弹性波逆时偏移的隧道超前探测成像技术研究[D]. 长春: 吉林大学, 2021
[22] 郭美延. 基于粘声介质的隧道不良地质地震波超前探测逆时偏移成像方法[D]. 济南: 山东大学, 2022
[23] Li Z, Qu Y. Research progress on seismic imaging technology. Petroleum Science, 2022, 19: 128-146
[24] 中华人民共和国国家标准. GB/T 33684—2017地震勘探资料解释技术规程. 北京: 中国标准出版社,2017
[25] 陆基孟, 王永刚. 地震勘探原理. 北京: 中国石油大学出版社, 2009