0 引言

射线走时和射线路径是构成零阶几何射线场(亦称为几何光学场)的两个基本要素。其中，射线走时对应高频波场的相位，射线路径对应高频波场的能量传播轨迹。因此，自从射线理论被引入到地震学中以来，对射线走时和射线路径计算方法的研究就一直没有间断过，尤其是在20世纪80年代以后，在反射地震偏移和反演研究的推动下^[1-2]，对射线走时和射线路径的研究一度成为研究热点，有很多研究者和工程师参与到其中，提出了很多方法和技术^[3-9]。形成这种局面的主要原因是基于射线理论的偏移和反演方法具有对输入数据体观测装置适应性强和计算速度快等优点^[1-2]。再有，在偏移成像中的射线追踪与常规的两点射线追踪有所不同，只须设定射线起点的位置，而不须设定射线终点的位置。这是因为基于射线理论的偏移在具体实施上是个一点对多点过程，任意一个成像点上的射线走时和射线方向均可利用相邻点上的相关信息通过插值得到。然而，由于这种只考虑初值的射线追踪方法不能在计算过程中随时插入由非均匀介质中的聚焦和反聚焦现象所产生的新射线(散射射线)，从而会在成像区域形成射线覆盖盲区。因此，在20世纪90年代以后，无论是工业界还是学术界均倾向于采用波前构建法来替代常规的初值射线追踪法，并逐步建立了基于波前构建的Kirchhoff型偏移和波束偏移^[1-2]。但是，波前构建法在实现时需要利用复杂的数据结构，而且对于多路径走时的存储也比较复杂^[10-14]。所以，在20世纪90年代以后，基于程函方程数值分析的射线走时计算方法逐步受到关注和重视^[1-5]。

在历史上，利用程函方程数值解计算地震波走时的基本思想早已有之，但是由于没有相应的需求而一直不为人所关注。1988年，Vidale首先将前人的基本思想付诸实践^[3-4]，提出了基于程函方程有限差分(FD)数值解的地震波走时计算方法。与上文提到的初值方法相比，FD法具有计算速度快和没有计算盲区等特点，从而可以作为一种理想的地震波走时计算方法替代传统的射线追踪法。然而，在无人为干预时，利用有限差分等基于网格的程函方程数值分析方法只能沿着波的前进方向计算走时，而不能沿着波的回折方向计算走时。换句话说，在无人为干预时，利用任何一种基于程函方程数值分析的方法均只能计算沿着波的前进方向的透射(波)走时，而不能同时计算由该透射波所激发的反射(波)走时，更不能计算由该透射波所激发的、其传播方向经过多次改变的多次透射(波)和多次反射(波)走时。针对这个问题，Rawlinson等^[5]在2004年提出了一种在人为干预下利用程函方程数值分析方法计算多次透射和多次反射走时的分区多级方案，在这种方案中，走时计算采用了基于窄带技术的快速推进法(fast marching method，FMM)^[5]。为计算多次反射走时，Rawlinson等^[5]将每一个带有光滑弯曲界面的非均匀层视为一个独立的计算单元，而每一个层面均被视为窄带(层面窄带)。在层面窄带中，快速推进过程被反复重置为初始状态(重新初始化)，并根据预定的传播方向重新进行推进计算。为了计算射线路径，Rawlinson等^[5]采用了逆向梯度法，即首先利用有限差分网格节点上的走时值计算走时梯度，然后再根据梯度方向计算射线路径的局部坐标，最后将所有的局部坐标点连成一体，得到射线路径。

2006年，de Kool等^[6]将上述计算方法推广到三维。2017年，孙章庆等^[7]将这种多级计算方法与群推进法组合在一起，提出了一种针对复杂山地多波型走时计算的多级次群推进迎风混合法。类似地，唐小平等^[8]将多级计算方法与最短路径算法组合在一起，提出了一种基于最短路径法的、针对三维层状介质中多次波的走时与射线路径计算方法。在理论上，分区多级计算的基本思想可以用于任何一种基于程函方程数值分析的地震波走时计算方法，因此存在有无限多的组合方式。

无论是在Rawlinson等的开创性工作中，还是在后续的跟踪完善性工作中，走时梯度都是直接利用网格点走时值来计算的。由此而产生的问题是：①梯度计算的精度和稳定性受有限差分网格常数的制约，并且随计算区域和网格密度的变化而变化，尤其是在复杂(曲)边界附近或是在变网格中更是如此；②梯度计算的精度和稳定性受单个网格点走时精度的制约。为解决这个问题，张东等^[15]采用B样条插值多项式，将对离散数据的数值求导问题转化为对分片光滑函数的求导问题。2014年，张婷婷等^[16]将B样条插值用于三维层状介质中的多次波走时和多次波射线路径计算，为在三维条件下研究多次反射提供了一个新的计算途径。然而，根据张婷婷等^[17]的计算结果，基于B样条插值多项式的计算方法会在速度突变区域内产生一定的误差，因而难以在速度剧烈变化区域内得到满意的计算结果。本文针对这个问题：①利用Chebyshev正交多项式逼近走时场，使分区多级计算的误差在最小二乘的意义下达到最小，从而避开速度突变所带来的走时突变问题；②对Chebyshev逼近函数求导，得到射线路径在所考虑区域内的局部坐标。

1 基于程函方程有限差分解的分区多级走时计算

根据Rawlinson等^[5]的工作，分区多级计算的主要步骤是：①在模型中划分出反射区和透射区。②将反射(透射)区内的网格点划分为接收点、活动点和远离点(图 1)，其中：接收点位于窄带以外的上风区，是已经完成走时计算的点；活动点位于窄带之内，是当前正在进行走时计算的点；远离点位于窄带以外的下风区，是还没有进行走时计算的点。③利用公式

计算所考虑活动点上走时梯度的模。式中：S_{i, j, k}为网格慢度；|∇t|_{i, j, k}为走时梯度的模；D_n^+lt_{i, j, k}和D_n^－lt_{i, j, k}(l=(x, y, z))分别表示点(i, j, k)处走时在对应方向上的n阶向前和向后差分算子。当n=1时，

式中，h为网格间距。④将最小走时值赋予所考虑的活动点，并将其标记为接收点，移出窄带。⑤将距所考虑活动点最近的远离点移入窄带，重复③—④所描述的计算过程，直到全部远离点均变为接收点。

图 2给出了上述分区多级计算步骤的进一步说明。①首先由震源开始，计算反射(透射)面及以上区域内所有网格点处的走时(图 2a)。②将界面(反射(透射)面)重置为窄带(图 2b)。③为计算反射波，将界面(反射(透射)面)以上的点设为远离点，按照上文中的计算步骤重启FMM(图 2c)。为计算透射波，将界面(反射(透射)面)以下层内的点设为远离点，按照上文中的计算步骤重启FMM(图 2d)。重复图 2所描述的计算过程，即可得到任意级次的多次反射和多次透射。

2 基于Chebyshev多项式的离散走时逼近

Chebyshev逼近是基于Chebshev多项式的一种逼近方法，是为解决零偏差问题而提出的一种基于正交多项式的最佳逼近方法^[18]。在历史上，Chebyshev逼近和Chebshev多项式由俄罗斯数学家P L Chebyshev在1854年提出，其最初目的是解决连杆设计问题^[19]。根据其数学表现形式可将Chebyshev提出的多项式分为两种，即第一类和第二类Chebyshev多项式。其中，第一类Chebyshev多项式起源于多倍角余弦和正弦函数展开，用T_n(x)表示。具体地讲，T_n(x)=cos[ncos^－1(x)]。第二类Chebyshev多项式用U_n(x)表示，与T_n(x)的关系为U_n(x)=[dT_n+1(x)/dx]/(n+1)。在实际应用中，利用下列递推公式计算T_n(x)和U_n(x)，即T₀(x)=1、T₁(x)=x、T_n+1(x)=2xT_n(x)－T_n－1(x)、U₀(x)=1、U₁(x)=2x、U_n+1(x)=2xU_n(x)－U_n－1(x)。

作为一种正交多项式，Chebyshev多项式与其他正交多项式相比，例如与Legendre、Laguerre和Hermite等正交多项式相比具有许多独特的性质^[9]。事实上，Chebyshev多项式特别适合于进行最佳平方逼近。此外，在一定条件下可将连续函数展开为Chebyshev级数等等。因此，为在最小平方意义下逼近来自于快速推进法的走时场，在下文中选取Chebyshev多项式为逼近函数。如图 3所示，在诸多正交多项式中，Chebyshev多项式的方差最小，即逼近误差亦为最小。

现在利用第一类Chebyshev多项式T_n(x)来构造3维离散走时τ(x, y, z)的近似解析表达式t(x, y, z)。为此目的，引入系数序列a′m，m=0, 1, …, 8，使得

展开后得到

式中，a_p(p∈0, 1, 2, …, 26)是利用最小二乘法得到的插值系数，即

式中：τ(x_i, y_j, z_k)为网格节点上的走时值；φ_p(x_i, y_j, z_k)(p∈0, 1, 2, …, 26)为逼近函数中与系数对应的正交项；M, N, L分别为x, y, z方向的网格节点数。

一旦完成了公式(4)所描述的逼近过程，即可根据走时函数t(x, y, z)的梯度函数∇t(x, y, z)求取当前追踪点的射线方向。根据定义，

3 数值计算实例 3.1 常梯度速度模型

首先考虑图 4中所示的常梯度速度模型。模型尺度为x×y×z=500 m×500 m×250 m，网格间距为10.0 m。震源点坐标是(1.0 m，1.0 m，1.0 m)；接收点均匀地分布在模型表面上，其间距为20.0 m；速度分布函数为v=500+60z。式中：z为深度；v为速度。

对比图 4a、b可见：Chebyshev逼近法所得到的射线路径在源点附近优于B样条插值法；在其他区域，两者所得射线路径在整体上相差不大。为便于进一步对比，取y=25.0 m处的xz剖面进行对比，其结果示于图 5，可以看出，Chebyshev逼近法所给出的走时具有较小的相对误差。另外，在一台CPU为酷睿i5的台式机(CPU频率为2 450 MHz；内存为6 G RAM)上计算时，B样条插值法所用的CPU时间为5.191 s，而Chebyshev逼近法所用的CPU时间为4.205 s，其计算效率相差约20%。因此，综合来看，在这个模型中，Chebyshev逼近法要优于B样条插值法。

3.2 Marmousi模型数值计算

为检验在上文中所提出的方法在复杂模型中的计算精度，首先利用三维Marmousi模型进行对比研究。图 6给出了三维Marmousi模型中一个剖面上的对比结果，不难看出：运动学射线追踪虽然能准确地计算出射线路径，但是存在射线路径盲区，从而不能求出模型内的全部走时；与此相反，在Chebyshev逼近法中，射线路径是通过走时来计算的，而且是程函方程的有限差分解，存在于整个速度模型之中，因此，在Chebyshev逼近法中没有射线覆盖盲区问题。值得指出的是，虽然利用基于程函方程有限差分法的Chebyshev逼近法能算出整个模型中的射线路径，但是考虑到过多的射线条数会影响对比效果，所以在图 6中只在射线追踪法的覆盖区内进行了对比。

3.3 多次透射一次反射

图 7a、b分别给出4层水平层状模型中的透射走时和在第4层底边界的反射走时。图 7c给出的是在第4层底边界利用Chebyshev逼近法得到的一次反射射线路径，其中透射部分的射线均垂直穿过在图 7a中所给出的等时线，反射部分的射线均垂直穿过图 7b中所给出的等时线。图 7d是利用Chebyshev逼近法求得的所有4个层面上的入射与反射射线路径。图 7e给出的是Chebyshev逼近法与样条插值法之间关于射线路径计算的相对误差对比，可以看出，Chebyshev逼近法的精度要高于B样条插值法的精度。

3.4 层间多次反射

为考察Chebyshev逼近法计算多次反射的效果，首先选取Marmousi模型中的一个小部分，构成局部层状介质(图 8)。由图 8可见：高速区射线路径偏向层面方向，而低速区射线偏向层面的法线方向。这与Snell定律所要求的完全一致，说明基于Chebyshev逼近的射线路径计算是正确的。

鉴于在图 8所示模型中，Chebyshev逼近和B样条插值之间的误差对比分析与图 7e中的结果类似，故在此不再进行专门的讨论。

其次，考察高速度反差条件下Chebyshev逼近法的计算效果。为此选用二维Sigsbee模型，其计算结果展示在图 9中。为了保证能在盐丘底部形成有效的多次反射，将震源置于远离高速体的位置，并以模型的底界面为一次反射面。首先，利用多级次FMM计算出整个模型中的多次反射走时，再从模型顶界面开始进行逆向梯度计算，进而得到如图 9所示的射线路径。为方便起见，只计算了模型底界面和盐丘底界面之间的2次反射射线路径。尽管如此，也可以清晰地看出在高速体存在时的层间多次反射始终满足Snell定律，且由高速体向低速区域呈“扫把”形出射。这进一步地说明了Chebyshev逼近法计算一次和多次反射射线路径的有效性。

4 结论

针对利用B样条插值进行走时梯度计算所存在的问题(在速度突变时会出现较大的计算误差)，本文提出利用Chebyshev多项式替代B样条对走时场进行逼近，并在此基础上结合多级快速推进法(FMM)建立了一种计算三维多次反射射线的方法。根据简单的水平层状介质和复杂的Marmousi和Sigsbee模型中的数值计算结果可以得到下列结论：

1) 利用Chebyshev多项式逼近走时可以有效地解决离散数据求导过程中出现的不稳定问题，以及B样条插值法在速度突变时所引起的大误差问题，并能较好适应分区多级次算法得到的多次反射、透射走时场。

2) Chebyshev逼近法在简单和复杂模型中均能保证所计算出的射线路径具有较高的精度，能够满足偏移成像的要求。