0 引言

地震波场数值模拟在地震资料采集、处理与解释中具有重要的地位，是地震勘探中的一个有力工具.地震波场数值模拟主要分为两大类：波动方程类数值模拟方法和几何射线类数值模拟方法.波动方程类数值模拟方法实质上是求解地震波动方程，因而模拟的地震波场信息较为丰富，包含了地震波传播的所有信息，但对计算机内存要求较高，计算速度相比于几何射线类方法要慢很多.几何射线类数值模拟方法属于几何地震学方法，由于它将地震波波动理论简化为射线理论，主要考虑的是地震波传播的运动学特征，其主要优点是概念明确，显示直观，运算方便，适应性强；其缺陷是缺少地震波的动力学信息，应用有一定的限制条件，而且存在阴影区和焦散区的问题.随着勘探的地质任务越来越复杂，地下地质体构造类型的多样性、复杂性，在波动方程类数值模拟计算量巨大、对硬件要求较高的情况下，射线追踪方法能发挥其计算快、适应性强的优势.

在射线追踪法中，有两类地震射线追踪问题：一类是一点射线追踪：即试射问题，已知射线初始点位置和初始出射方向求地震波的传播路径问题；另一类是两点射线追踪：即已知射线初始点和另一接收点的位置，不知射线初始出射方向，求两点之间的射线路径问题.

为克服普通射线法和波动方程方法的某些局限，Červeny等人相继发展了一种将波动方程和射线理论相结合的方法，被称做高斯射线束法 (Červeny et al., 1982; Červeny and Pšenčík, 1983, 1984; Červeny, 1985；Klimeš，1984；Weber，1988；Hill，1990).该方法同时考虑波的运动学和动力学特征，适用于复杂的非均匀介质模型，还能考虑介质的吸收作用 (Weber，1988).该法不需要两点射线追踪的计算，具有速度快、精度较高 (可与有限差分的计算结果对比 (George et al., 1987)) 的特点，对焦散区、临界区及暗区等奇异区域都具有较好的效果.

国外方面，有关高斯射线束方法的文章较多，不少研究者对于高斯射线束法合成记录的某些参数选取 (Müller，1984；Weber，1988)、合成记录的精度 (Weber，1988；George et al., 1987) 进行了比较深入的研究，对于高斯射线束合成记录的运动学和动力学射线追踪也有不少文章进行了介绍 (Červeny et al., 1982；Müller，1984).MüHill (1990, 2001) 还提出了高斯射线束偏移的方法.在高斯射线束正演模拟方法的基础上，Cruz等 (2012)提出了投影菲涅尔带约束的高斯射线束正演模拟方法，通过投影菲涅尔带约束，使得高斯射线束在传播过程中保持聚焦和稳定.尽管这些方法还有一些问题 (如参数选取、射线密度等) 值得研究，还需要不断完善和更新，但总的看来，它已经形成了地震模型正演的一个方向.

国内方面，关于高斯射线束方法的研究起步较晚.李瑞忠等 (2006)首先介绍了高斯束方法的基本原理及其在波场的正演模拟和偏移成像方法中的应用；邓飞等 (2009)发展了三维射线快速追踪算法并将其应用于高斯束正演方法，提高了正演的精度和计算效率；孙成禹等 (2011)提出了能量约束下的高斯射线束地震波场正演方法；李辉等 (2015)基于高斯束与高斯波包的Gabor框架实现了散射波模拟方法；黄建平等 (2015)基于有效邻域波场的近似框架提出了适用于三维起伏地表情况下的高斯束正演模拟方法；印兴耀等 (2016)将基于Gabor分解的子波重构技术应用于时—空域高斯束波场计算中，提出了一种可以模拟任意点源函数透射与反射波场的地震正演方法.

本文主要介绍传统射线类正演模拟方法、高斯射线束正演模拟方法和菲涅尔带约束高斯射线束正演模拟方法的原理，并给出相应的模型正演结果进行对比分析，讨论其各自的优缺点及适用范围.

1 基本原理 1.1 传统射线类正演模拟方法

从射线理论可知, 波在介质中传播满足Snell定理、Fermat原理和Huygens原理，所有射线追踪方法基本上是围绕这三个基本原理展开的.

试射法是最早提出且使用最普遍的一种射线追踪方法.其射线追踪过程是：在激发点，给定一系列射线参数初始值，然后根据Snell定理依次进行追踪，在接收点附近选择最接近的两条射线，通过内插，调整初始射线参数值，经过多次的调整修改，可获得满意的结果.这种方法的最大优点是实现了射线的精确追踪，能够避开在盲区中追踪；在相对简单的模型结构中迭代收敛较快，但在复杂结构中收敛较慢，比较耗时，这是该方法的不足.

从数学上看, 试射追踪方法是一个给定初值条件的定解方法, 初值条件是震源激发点的位置和初始射线方向 (初始角).当地质模型建立后, 射线从震源点沿初始射线方向在介质中传播，当遇到中间界面时，射线依据Snell定理透射到新的介质中, 当遇到目标界面时，射线依据Snell定理进行反射，最终到达接收地表；射线在地表的出射点往往并不在接收点上，要获得接收点的射线路径，需要通过不断修改初始射线参数，进行多次重复试射追踪.

如图 1所示为层状介质试射法得到的射线路径，模型共四层，射线路径根据Snell定理进行反射和透射.从图中可以看出，初始条件相同，但每层的出射线分布范围却不同，界面越深，偏离震源点越大，覆盖范围更广；对同一界面，射线入射角度增量相同，出射线分布的点距越来越大.

1.2 高斯射线束正演模拟方法

高斯射线束是一种将波动理论与射线理论结合起来的方法，在运动学射线追踪过程中，同时进行动力学射线追踪，然后根据所求的运动学和动力学参数计算接收点处的走时和振幅，最后将每条高斯束的有效能量叠加起来形成最终的合成记录.该方法即能反映波场的运动学特征，又能体现其动力学特点，且无需费时的两点射线追踪，具有运算速度快、精度高的优点，而且对临界区、阴影区都有较好的效果.

所谓高斯射线束就是指弹性动力学方程集中于射线附近的高频渐近时间调和解，根据Červeny等 (Červeny et al., 1982)，在如图 2所示的射线中心坐标系中，高斯射线束解的波场主分量在频率域具有如下表达形式为

(1)

式中，s代表中心射线Ω上某点到任意参考点的弧长，n代表Ω附近一点到s点的距离；坐标系的基矢量分别为同射线Ω相切的单位切向量t和同Ω垂直并指向Ω同一侧的单位法向量n；ω为圆频率；U_P为P波主分量；i为虚数单位；t为时间参量；τ(s) 为中心射线Ω上某点的走时；τ(s) 和M(s) 可表示为

(2)

而p(s) 和q(s) 为随射线传播而变化的两复值动力学射线参数，其满足动力学射线追踪方程，公式为

(3)

式中，；A(s) 为中心射线上振幅值.

对于层状介质，假定一条射线从炮点s₀出发，经过N个界面后反射，回到地面接收点R，则该振幅的表达式为

(4)

式中，R_i为界面反射或透射系数；ρ为介质密度；Q为某界面的反射点或透射点；v(R) 为传播速度；符号“-”表示在射线通过界面时生成射线 (反射或透射射线) 一侧的量值；α_i和β_i分别为入射射线和生成射线与局部坐标X轴正方向 (界面切线方向) 的夹角.

为方便对高斯束的基本性质与特点进行简单分析，忽略 (1) 式中的exp (-iωt)，并分离的实部和虚部，则可以得到具有更为明显的物理意义的高斯束表达式为

(5)

其中，为高斯束的波前曲率；L(s)=为高斯束同频率有关的有效半宽度.从 (5) 式可以看出，沿射线的动力学射线参数p(s) 和q(s) 要满足两个条件：第一，q(s)≠0，该条件保证高斯束沿射线是处处正则的 (振幅有限)；第二，，该条件保证解集中于中心射线附近.为了满足以上条件，在初始点s₀处选取两组互相独立的初始条件为

(6)

根据龙格-库塔法可以解出射线上任意一点s处的两组线性独立解为

从而可以构造得复数解为

(7)

其中，ε为任意复常数，ε的选取对高斯射线束具有重要意义，它决定了高斯束的波前曲率和有效半宽度，即决定了高斯束的性质.有很多文章对于ε的选取进行了详细的探讨，如Müller (1984)，Weber (1988).特别是Weber (1988)的文章对于比较深入的研究，并用理论模型实验说明了ε的选取对合成地震记录精度的影响，其提出在计算中采用如下选择使得射线束在终点具有最小的半宽度，对计算有利

(8)

通过运动学和动力学射线追踪，可以得到沿中心射线附近的弹性波动方程高频近似P波位移分量在频率域的表达式—高斯射线束，把从震源展开的各高斯射线束在接收点处叠加并通过傅里叶变换即可得到时间域的地震波场

(9)

正演模拟得到的地震记录为式 (9) 的实部.式中，u(R，t) 表示时间域的地震波场，R为接收点；ω为圆频率；Φ(φ) 为与初始入射角有关的权函数，φ是射线的入射角；F(ω) 为所选用震源函数f(t) 的频谱，f(t) 应当为高频函数，以保证高斯射线束法的高频近似性；u_φ(R，ω) 为初始入射角为φ的高斯射线束.Φ(φ) 确定方法是将式 (9) 渐近值与二维线源波动方程精确解的渐近值比较而得出.根据Červeny等 (Červeny et al., 1982) 有：

(10)

v₀为震源处速度.

常用的高斯射线束地震记录合成方法有频谱法、褶积法和波包法.其中，波包法直接在时间域将射线能量进行叠加，效率最高.波包法合成地震记录的离线表达式为

(11)

式中，Δφ为中心射线入射角的间隔.式 (11) 表示了如下的物理意义：波包g(R, φ) 从震源以入射角φ发出，这些波包沿射线传播且与射线紧密相连，它们随着射线的传播连续改变其特性，介质中任意一点的波场是这些波包在该点的叠加.

波包的计算即可以用傅氏变换法，也可以用褶积法求取，下面用一种具有解析形式计算波包的子波函数，该子波函数由下式给出：

(12)

式中参数f_m，γ，v，可根据需要选择，f(t) 对应具有高斯包络的谐波载体，其中，γ控制包络的宽度，f_m为其主频率，该子波一般称作Gabor子波.利用该子波形式，最终可得到计算波包的近似解析表达式为 (Červeny and Pšenčík，1983)：

(13)

式中，f^*=f_m(1-4πf_mG/γ²)，上式为一个近似公式，适用于2πf_mG/㲫1的情况.从式中可以看出高斯波包的主频为f^*，并且在时间和空间都具有高斯包络的特点.

综上所述，用高斯射线束计算波场时，可分为三步：

(1) 做运动学射线追踪，通常的做法是将程函方程化为介质中射线轨迹的微分方程组，然后用龙格-库塔法求足够密的射线轨迹.

(2) 做动力学射线追踪，沿射线求高斯射线束的解，即对多条射线求解微分方程组得到和的值.

(3) 对检波点附近的高斯射线束的贡献进行加权叠加，求得波场值.

1.3 菲涅尔束正演模拟方法

根据Červeny等 (1982)和Müller (1984)可知，高斯射线束的有效半宽度可以表示为

(14)

式中，ω_ref为参考频率；ε为复值初始束参数.

由菲涅尔体射线追踪理论 (Červeny and Soares, 1992) 可知，第一菲涅尔带半径可以近似表示为

(15)

该表达式在震源处会存在不稳定问题，结合均匀介质中的菲涅尔带半径的解析解，可以将其修正为

(16)

式中，λ_avg=2πv_avg/ω_ref，v_avg为模型的平均速度.

根据惠更斯-菲涅尔原理可知，中心射线附近的波场能量主要分布在第一菲涅尔带内，为使高斯束的能量分布符合这一原则，可以使用如下的约束条件：

(17)

将 (14) 式和 (16) 式带入 (17) 式，并考虑到ε=ε₁+iε₂，ε₁=0，和q₁p₂-q₂p₁=1，可以得到：

(18)

与常规高斯射线束方法不同的是，ε不再是一个复常数，而是随射线弧长动态变化的函数，这是菲涅尔带约束高斯射线束正演方法与高斯射线束方法的主要区别.

为了方便进行讨论，将 (18) 式中初始束参数决定的地震波束称为菲涅尔束 (杨继东等，2015)，可以给出该种情况下高斯射线束在频率域的表达式为

(19)

式中，ε(s) 如 (18) 式所示.为了对比高斯射线束与菲涅尔束的传播形态，在均匀介质 (速度为2000 m/s) 情况下，两者传播形态以及有效半宽度对比如图 2所示.通过对比可知，高斯射线束的初始有效半宽度约为，并且在后续的射线路径上以双曲规律增大；而菲涅尔束的初始有效半宽度约为1/4λ，并且在后续的传播路径上以平方根规律增大.因此，相对于高斯束来说，菲涅尔束在整条射线路径上都具有较小的有效束宽度，这使得中心射线附近波场的走时和振幅计算更加准确，从而提高了波动方程解的精度.

利用菲涅尔束代替高斯射线束进行正演模拟，其正演流程与高斯射线束正演完全相同.

2 模型试算

以上介绍了三种射线类正演模拟方法的基本原理，下面通过两个数值算例对比分析三种正演方法的应用效果.

2.1 水平层状模型

首先采用简单的层状介质模型对三种正演方法的正确性进行测试.图 4为水平层状模型及射线路径示意图，水平层状模型大小为8000 m×4000 m，存在三个水平反射界面，分别位于1000、2000、3000 m深度处.观测系统采用中间放炮的激发方式，炮点位置位于4000 m处 (地表黄色星号)，320道接收 (地表绿色星号)，道间距为25 m，记录长度为4000 ms，时间采样间隔为2 ms.图 5a—图 5d分别为利用波动方程有限差分法、试射迭代法、高斯射线束和菲涅尔束方法合成的炮记录，由图可见：波动方程的正演结果存在直达波，而三种射线类正演方法未考虑直达波；四种方法得到的反射波同相轴的位置和形态基本一致.图 6为图 5正演结果的单道提取结果对比.由图可见：在地震波传播过程中，地震子波的振幅和相位均在发生变化，振幅变化主要是由反射、透射以及几何扩散引起，相位变化主要是由介质的非均匀性引起；高斯射线束方法和菲涅尔束方法得到的炮记录与有限差分方法得到的炮记录中反射波的波形及其振幅变化规律基本吻合 (两者采用不同的子波)，而射线正演方法的正演结果未考虑反射、透射和几何扩散的影响，表明高斯射线束正演方法和菲涅尔束正演方法对反射、透射和几何扩散等波场现象的刻画较为准确，而射线正演方法未能刻画反射波的振幅变化；水平层状模型试算证明了三种射线类方法的正确性.

2.2 实际模型——大王庄模型

为了进一步对三种射线类正演方法的适应性进行测试，本文采用实际模型——大王庄模型进行测试.大王庄模型及射线路径如图 7所示，模型大小为20000 m×6000 m，该模型为东部某探区简化得到.观测系统采用中间放炮的激发方式，炮点位置位于10000 m处，480道接收，道间距为25 m，记录长度为6000 ms，时间采样间隔为2 ms.图 8a—图 8e分别为利用波动方程有限差分法、试射迭代法、高斯射线束、菲涅尔束方法合成的炮记录以及实际资料单炮记录，由图可见：四种方法得到的反射波同相轴的位置和形态基本一致；高斯束正演方法和菲涅尔束正演方法对反射波的刻画更为准确，模拟得到的反射波同相轴更加连续且与实际资料单炮记录的同相轴形态较为接近 (图中红色箭头和方框区域).图 9为图 8正演结果的单道提取结果对比.由图可见：在复杂构造情况下，有限差分方法合成的地震记录在有效反射能量之外存在能量较强的毛刺信号 (图 9b红色矩形框所示)，这些毛刺信号是有限差分算法选取的网格尺寸太大、算子精度不够引起的数值噪声以及复杂构造中绕射点产生的绕射波共同引起的结果；由于射线类正演方法不能模拟这类绕射波能量，因此在合成记录中不存在这些能量较强的毛刺信号.通过实际模型试算结果表明，高斯射线束和菲涅尔束正演方法可以较为准确地刻画复杂构造中的地震波传播现象，精度较射线正演方法有明显提高，且通过菲涅尔带约束，高斯束传播更为稳定，更加聚焦，可以一定程度上提高高斯射线束正演精度.

以上两个模型用四种方法进行正演模拟的时间对比见表 1.由表可见，相比于有限差分方法，射线类正演方法的计算效率有极大的提高；高斯射线束和菲涅尔束正演方法的计算效率接近于射线正演方法，但其计算精度有了较大程度的提高.

3 结论与讨论 3.1

解剖了三种射线类正演方法的基本原理，并通过利用第一菲涅尔带半径约束高斯束束宽发展了菲涅尔束正演方法.在此基础上，对三种射线类正演方法进行了对比分析和模型试算，通过测试结果可以得到以下认识：(1) 射线类正演方法相比于波动方程有限差分方法计算效率有极大的提高，在保持一定计算精度的基础上，射线类正演方法更加适用于野外观测系统设计，尤其是三维地震勘探；(2) 高斯射线束正演方法的计算效率基本与传统射线正演方法相当，但其计算精度却有较大程度的提高，地震波场特征更加接近于波动方程方法；(3) 菲涅尔束正演方法通过菲涅尔带约束使高斯束传播更加稳定且能量更加集中，在第一菲涅尔带范围内具有更加明确的物理意义，使得模拟的高斯束波场更接近于波动方程方法，在一定程度上提高了射线类地震波场的模拟精度.

3.2

通过以上分析可以发现，高斯射线束和菲涅尔束正演方法在计算效率和计算精度两方面均满足实际生产的需要，在观测系统设计方面具有较好的应用前景.

致谢 感谢审稿专家提出的宝贵修改意见和编辑部的大力支持！