射线追踪计算在地震层析成像、 地震波场正演数值模拟、 地震偏移成像等领域有着非常广泛的应用，是一项基础性的研究. 针对不同的应用情况，前人开展了大量的研究工作，目前已存在多种射线追踪方法，不同的方法各有其优缺点. 如传统射线追踪算法中的打靶法和弯曲(或伪弯曲)法： 当检波器间距较小，且需要追踪大量地震射线时打靶法是有效的； 当震源和检波器连线与震源出射方向的夹角较大时，弯曲法优于打靶法. 但两种方法具有共同的缺点： ① 时常陷入局部解的困扰； ② 当检波器位置变换时需要重新进行射线追踪，因而增加了CPU时间.

目前使用较为广泛的射线追踪方法之一，是基于费马原理和惠更斯原理的最短路径方法(shortest path method，简写为SPM)(Moser，1991). 其优点是算法简单、 数值计算稳健、 全局解. 最短路径方法最早是由Nakanishi和Yamaguchi(1986)提出的，其原理是将速度模型剖分成一系列的网格单元，用网格节点之间的连线近似地震射线路径. 随后不少学者开展了大量的研究工作，主要基于以下三方面： ① 提高计算精度； ② 减少计算用时； ③ 扩展功能(即追踪后续波).

在提高计算精度方面国内外学者做了不少尝试，例如： 增加射线角度覆盖率的Forward Star算法(Moser，1991； Klime，Kvasnika，1994)，减少无用节点的最小走时树算法(Cao，Greenhalgh，1993)，克服之字型射线路径网格节点间引入弯曲的算法(Fischer，Lees，1993； Van Avendonk et al，2001)，动态网格射线追踪算法(张建中等， 2003，2004)，以及减少单元及节点总数的三维改进型最短路径算法(modified shortest-path method，简写为MSPM)(Bai et al，2007).

计算速度的提高主要是在保证计算精度的前提下，尽量减少参与计算节点的总数以及最小走时的挑选方式，例如： 界面网格全局算法(刘洪等，1995)，基于图型结构的射线追踪算法(王辉，常旭，2000)，非均匀介质中射线追踪快速算法(赵爱华等，2000)，双重网格法(赵爱华，丁志峰，2005)，桶排序法(张美根等，2006)，以及三维改进型的最短路径算法(Bai et al，2007).

后续波射线追踪的进展主要有： 二维模型中反射波的追踪计算(Moser，1991)，PS转换波的追踪算法(赵爱华，张中杰，2004)，宽角反射地震波走时模拟的双重网格算法(赵爱华，丁志峰，2005)，PS转换波追踪的界面二次源算法(王童奎等，2007)，以及复杂层状介质中多次反射、 透射及转换波的追踪计算(分区多步改进型最短路径算法)(Bai et al，2009； 唐小平，白超英， 2009a，b).

上述算法的改进在模型划分时大多局限于规则网格，即在模型单元划分时采用了规则(均匀)的概念. 通常模型是由相同的立方体(或长方体)单元组成. 计算精度的提高主要是通过不断缩小单元尺度或增加次级节点的密度来实现的. 然而，实际问题中人们不可能将模型的网格单元划分得过密，否则将导致模型单元和节点数过多，造成计算量过大. 对于复杂地学模型这种规则划分显然是不适用的. 模型单元的常速化，很难引入速度不连续界面. 然而这些不足在模型规则划分的前提下是不可能彻底改变的.

所以必须引进不规则网格的概念，在这方面有关学者做过一些相关的研究，例如： Vesnaver(1996)采用不规则网格单元实现了透射、 反射、 转换波的追踪计算，其中使用了互易原理，即从炮点和检波器开始进行下行波的计算至速度界面，随后由费马原理的稳态时间路径和最小走时原理得到界面处的反射点. Sambridge和Gudmundsson(1998)则采用了更一般的不规则网格单元划分来进行两点间的射线追踪计算.

本文总结了前人的经验，在模型建模时采用另一种不规则网格单元划分方式(以下简称ISPM)，用分区多步技术计算二维、 三维复杂介质(含不规则起伏高、 低速夹层，复杂界面及地表起伏)中的多次透射、 反射及转换波. 采用不规则网格单元的好处有： ① 可以有效地模拟地表起伏或速度间断面; ② 根据不规则网格单元的形状引入次级节点，可有效地减少总节点数，从而提高运算速度; ③ 在不规则单元内可定义速度分布，如莫氏速度不连续界面.

我们曾运用分区多步计算技术与MSPM算法相结合，实现了模型规则网格划分下二维和三维复杂层状介质中多次波的追踪计算(唐小平，白超英， 2009a，b). 本文吸取了该法在计算多次波时的优点，不同的是在模型参数化时采用一种不规则网格单元进行划分，在保证计算精度的前提下提高了运算速度和计算精度.

首先对速度模型参数化，用两种形式进行网格单元剖分，二维模型如图 1所示. 在起伏地表和速度界面处依据各自起伏情况形成非规则网格单元，如梯形网格(图 1c和图 1d)及三角形网格(图 1e)，同时梯形网格单元又可细分成长方形及三角形单元，而远离不规则界面的区域依然用长方形(或正方形)单元. 单元的角点为主节点，主节点之间等间距插入次级节点，次级节点只分布在网格单元的边上，单元内无节点，但炮点和检波器可位于其内.

图 1 含有起伏界面及低速夹层的二维层状速度模型中分区多步ISPM算法计算透射、 反射及转换波示意图(如图 1a所示，模型中次级节点没有给出). 图(b)，(c)，(d)，(e)为模型参数化时采用的4种规则网格及不规则网格. 图(b)—(e)中空心圆表示主节点，实心圆表示所加入的次级节点 Fig. 1 A schematic explanation for computing the transmitted，reflected and converted arrivals using the multistage ISPM in 2-D layered media with irregular interfaces and low velocity layer involved. In diagram(a)note that the secondary nodes are not shown for clarity. Diagrams(b)，(c)，(d) and (e)are 4 different kind of regular and irregular cells used in the model parameterization. White and black circles indicate primary and secondary nodes，respectively

三维模型如图 2所示，可采用规则网格(图 2c)及非规则网格(图 2b). 单元的角点为主节点，在单元的面上插入等间距的次级节点，单元内部无节点. 非规则单元情况下次级节点的间距与规则单元时相同，如图 2b所示.

图 2 模型中主节点和速度界面分布(a)、 模型中不规则网格单元(b)和规则网格单元(c)图(b)和(c)中实心圆表示主节点，空心圆表示所加入的次级节点Fig. 2 Distribution of primary nodes and interfaces in velocity model(a)， and one irregular cell(b) and one regular cell(c). In diagram(b) and (c)black circles are primary nodes and white circles are secondary nodes

根据速度不连续界面的形状将其离散化，为保证算法的精度，离散时界面上离散点的间距应不低于单元面上次级节点的间距，离散点的速度值可由主节点的速度采样值通过线性拉格朗日或B样条插值函数给出，具体计算方法将在下文中给出.

模型分区是指按照模型的地表起伏和速度界面的具体情况将速度模型分成相对独立的层状(起伏)区域，相邻区域由速度界面相连接. 多步计算是指按照所要追踪后续波类型，从炮点所在的区域开始，沿着目标射线所要通过的区域，逐区进行波前扫描. 具体来说，自炮点所在的单元开始进行扫描，等当前单元所有点的走时计算结束后，按照一定的步长进行波前扩展. 由于本文中存在规则和不规则网格单元，所以扩展时两种情况分别进行计算. 待当前区域所有网格单元内节点都计算完毕后，波前停留在该区的速度离散界面上，同时保存速度界面上各离散点的走时和相应的射线路径. 举例来说(如图 1a所示)，如果是追踪透射波P₁P₂(约定： 上标代表上行波，下标代表下行波，数字1和2表示区号)，则调用2区速度模型(如是透射转换波P₁S₂则调用S波速度模型)，从所保存的速度界面中走时最小的点开始(惠更斯原理)，继续2区内波前扩展和射线追踪. 如果是计算反射波P₁P¹(或者是反射转换波P₁S¹)，则调用1区域相应的速度模型，按照上述步骤重复则可实现多次波的追踪计算.

下面介绍计算最小走时和速度插值的方法. 计算节点i，j两点间的最小走时公式为

若i节点不是主节点，且位于不规则单元中，则按以下步骤计算： ① 依据所求点的y轴坐标，在三维不规则网格y轴处做一个切片，转化为二维不规则网格; ② 所求点速度值分两种情况进行计算. 如图 3所示，不规则网格单元大体上分为4种类型，每种情况下网格单元均可分成两部分，如果i节点属于矩形部分，则按公式(2)的二维形式进行计算. 即

图 3 不规则网格单元中4种不同类型(a，b，c，d)下速度插值计算示意图(实心圆代表所求点)Fig. 3 A diagrammatic explanation for computing velocity interpolation in four kinds of the irregular cells. Black dot is the point to be calculated

如果所求点属于不规则部分，先用坐标点①和④(图 3a，c)或者②和③(图 3b，d)按一维线性拉格朗日插值函数得到坐标点⑤的速度值，然后分类采用B样条插值公式计算.

目前基于网格单元的多次波射线追踪主要有分区多步FMM算法(Rawlinson，Sambridge，2004; De Kool et al，2006)和分区多步MSPM算法(Bai et al，2009； 唐小平，白超英， 2009a，b). 对上述两种算法的计算精度与CPU时间进行系统对比，其基本结论是在相同网格间距下，分区多步MSPM算法的计算精度要高于分区多步FMM算法的计算精度，而CPU运行时间则快3—4倍(Bai et al，2009). 这里我们将对分区多步FMM算法与ISPM算法进行对比分析.

首先采用均匀模型，模型大小为100 km×40 km，速度为4.0 km/s. 地表为水平界面，模型30 km处有一水平反射界面. 在模型参数化时，FMM算法中正方形网格单元大小分4种尺度，分别是1 000，500，250 m和25 m. 而ISPM算法中网格单元大小为4 km×4 km，但反射界面上、 下采用三角形单元，也分4种情况在单元边上相应加入3，7，15和31个次级节点，以保证次级节点间距与FMM算法中相应单元节点间距相同. 炮点放在 模型的左上角(0.0，0.0)，100个检波器等间距(1.0 km)置于地面，最小炮检距为1.0 km.

图 4给出了随炮检距变化FMM算法(图 4a)与ISPM算法(图 4b)数值解相对于解析解在4种不同节点间距下反射波走时的相对误差. 表 1给出了两种方法与解析解的最大(或平均)走时误差，最大(或平均)相对误差，以及所用CPU时间(所用PC计算机： DELL Optiplex 755—2.53 GHz). 从图 4及表 1中可以看出，除了节点间距为1 000 m的情况外(此种情形不常用)，ISPM算法的精度是优于FMM算法的. 在ISPM方法中有些点的数值解与解析解误差非常小甚至等于0，所以ISPM算法的误差曲线呈波浪型. 而FMM算法的误差曲线则随炮检距缓慢增加，最后趋近于平稳. 一般而言，ISPM算法的计算速度要比FMM算法提高一个数量级以上.

图 4 随着炮检距变化FMM算法(a)与ISPM算法(b)数值解相对于解析解在4种不同节点间距下反射波走时的相对误差 Fig. 4 Reflected wave travel time errors of the FMM(a) and ISPM(b) solutions against the analytic solutions with varied source-receiver distance in four different grid intervals

表 1(Table 1)

表 1 均匀速度模型中由FMM和ISPM算法所得的反射波走时与解析解的走时误差、 相对误差和CPU时间在4种不同网格节点间距下对比表Table 1 The time errors and the relative time errors of the computed reflecting wave travel times using FMM and ISPM methods against the analytic solutions， and the corresponding consumed CPU times in four different grid spacing

表 1 均匀速度模型中由FMM和ISPM算法所得的反射波走时与解析解的走时误差、 相对误差和CPU时间在4种不同网格节点间距下对比表 Table 1 The time errors and the relative time errors of the computed reflecting wave travel times using FMM and ISPM methods against the analytic solutions， and the corresponding consumed CPU times in four different grid spacing

模型如图 5所示，大小为108 km×40 km，P波速度自地表的4.0 km/s向下线性增加至7.0 km/s(速度梯度为0.075/s). 有两个起伏界面，深度分别是20 km和36 km，界面起伏为±2 km. 模型参数化同2.1节，也分4种不同节点间距． 炮点置于(x=6.0，z=0.0)处，100个检波器等间距(1.0 km)排列于地表，从(7，0)至(106，0)，最小炮检距为1.0 km. 按图 5所示计算多次透射、 反射波的走时. 表 2给出了在4种不同节点间距下两种方法计算图 5所示的多次透射、 反射波走时的最大(或平均)走时差、 最大(或平均)相对值及CPU时间. 从表 2中可以看出，除了1 000 m节点间距情形，两种方法的走时差与相对值基本相同，但CPU用时差逐渐增大. 一般而言，FMM算法的CPU用时是ISPM算法的4—6倍.

表 2(Table 2)

表 2 4种不同节点间距下，FMM与ISPM算法在一维线性增加速度模型中计算图 5所示多次透射、 反射波走时的最大时间差、 平均时间差、 最大相对值、 平均相对值及CPU时间对比表Table 2 The maximum and averaged time differences， and the maximum and averaged relative values，computed using FMM and ISPM approaches，respectively，for multiply transmitted and reflected raypaths，as indicated in Fig. 5， and the corresponding CPU times，in four different grid spacing

表 2 4种不同节点间距下，FMM与ISPM算法在一维线性增加速度模型中计算图 5所示多次透射、 反射波走时的最大时间差、 平均时间差、 最大相对值、 平均相对值及CPU时间对比表 Table 2 The maximum and averaged time differences， and the maximum and averaged relative values， computed using FMM and ISPM approaches，respectively，for multiply transmitted and reflected raypaths，as indicated in Fig. 5， and the corresponding CPU times，in four different grid spacing

图 5 一维线性速度模型及多次反射、 透射波射线路径(图中不同灰色区域代表不同的计算区域，黑色粗线代表反射界面，黑色细线段代表随后计算的射线的路径) Fig. 5 1-D linear velocity model and multiply transmitted and reflected raypaths. The graded grey layers indicate different computational domains，bold black curve denotes reflecting discontinuity， and fold lines show raypaths

为了详细说明FMM算法与ISPM算法的走时差[FMM(t)—ISPM(t)]，图 6给出了两种算法在4种不同节点间距条件下计算上述多次透射、 反射波走时的时间差随炮检距变化的情况. 从图 6中可以看出，除了节点间距为1 000 m的情形外，ISPM算法计算的走时均小于FMM算法，且走时差的大小与反射界面的起伏有很好的相关性. 有关这一现象我们将在以下详细分析.

图 6 不同节点间距条件下FMM与ISPM算法计算多次透射、 反射波(图 5所示)走时差随炮检距变化曲线Fig. 6 Time differences of the computed multiply transmitted and reflected raypaths，as indicated in Fig. 5，with FMM and ISPM methods in four different grid spacing

为了说明上述两种算法在计算多次透射、 反射波走时差的大小与反射界面起伏的相关性，图 7给出了节点间距最小(125 m)时两种算法计算多次透射、 反射波走时-距离曲线. 从图 7中可以看出，ISPM算法的数值解系统地小于FMM算法的数值解. 理论上讲，在均匀或一维线性速度模型中FMM和ISPM算法的数值解要稍大于(或等于)解析解. 换句话说，数值解越小的算法其计算精度也就越高. 造成上述走时差随炮检距的增加而出现的转折(或拐点)现象，主要是由于所给定的不规则界面造成的，即在反射界面凸起的部分反射点汇集，而在反射界面凹陷的部分反射点发散(或者根本没有反射点). 这与反射界面水平时的情形截然不同. 总体来说，ISPM算法计算精度要高于FMM算法，其CPU时间仅是FMM算法的1/5—1/6倍，与使用规则网格单元划分的MSPM算法相比，计算速度也略高(参见Bai等(2009)文中表 3).

图 7 一维线性速度模型中使用最小节点间距时FMM与ISPM算法计算多次透射、 反射波走时(t)-距离(x)曲线图 Fig. 7 The t -x curves for the FMM and ISPM approaches with the finest grid spacing in computing the multiply trans-mitted and reflected raypaths in the 1-D linear velocity model

以Marmousi模型为例(图 8)，模型参数化时规则单元大小为24 km×24 km，而在速度反差较大的区域则采用了不规则单元，单元边上次级节点间距为2 m. 炮点放在了模型顶面的中心点位置，10个检波器均匀分布于地表，最小炮检距为919.2 m. 根据Marmousi模型的速度分布，选取如图 8的不规则界面. 图 8给出了自炮点产生的P波波前传播至不规则界面，经反射后传播至地表时的等时(0.08 s)波前及自炮点经不规则界面反射后到达地表 10个检波器的反射波射线路径. 从图 8中可以清晰地得到3点基本结论： ① 波前传播等振面在高速区被拉伸，而在低速区被压缩； ② 射线的传播方向总是与波前面垂直(即满足Snell定律)； ③ 由于存在高反差速度分布，即使在水平反射界面上也出现了散射点.

图 8 Marmousi模型中自炮点至不规则界面的下行波前与反射波射线路径(上图)，以及自不规则界面反射后上行波前与反射波射线路径(下图). 波前间隔均为0.08 s Fig. 8 The downward wavefront originated from a point source to the irregular interface and related reflected raypaths(upper diagram)， and the upward wavefront reflected from the irregular interface and the entire reflected raypaths(bottom diagram). The wavefront time interval is 0.08 s

所用模型如图 1a所示，模型长为100 km，深度为50 km. 地表为起伏界面，其内有3个起伏反射界面，深度分别为19，26和47 km，界面起伏幅度分别为±2，±2 km和±3 km. P波速度自地表的4.0 km/s向下线性增加至7.0 km/s，梯度为0.06/s，第一层界面与第二层界面之间为低速层，速度为4.8 km/s，与图 1a中的背景速度值相差0.5 km/s. S波速度由公式v_P/v_S=1.732给出，并假定其速度结构和界面起伏与P波相同.

图 9给出了一炮多道接收时的情形，炮点位于0.0 km、 -10.0 km处，10个检波器沿起伏地表等间距(10.0 km)排列，最小炮检距为10 km. 常见的8种地方及近区域地震震相分别为： 直达P波，射线路径如9a所示，受低速区的影响，射线下弯变得平展； 多次透射并分别经第一、 二、 三层界面的反射波，震相分别为： P₁P¹，P₁P₂P²P¹，P₁P₂P₃P³P²P¹(图 9b)； 经炮点激发S波上行至地表，转换成P波后通过多次透射分别经第一、 二、 三层界面的反射波，震相分别为： S¹P₁P¹，S¹P₁P₂P²P¹，S¹P₁P₂P₃P³P²P¹(图 9c)，以及多次透射、 转换、 反射波震相S₁P₂P₃P₃S²P¹. 3.3 三维层状起伏模型(含低速夹层)中多次透射、 转换及反射波射线追踪

三维层状起伏模型，如图 2所示. 模型长、 宽均为100 km，深度为50 km. 模型顶面为起伏地表，地表下有两个起伏反射界面，深度分别是25 km和48 km，起伏幅度分别为±1 km和±2 km. P波速度自地表的4.0 km/s向下线性增加至8.0 km/s，梯度为0.08/s. 第一层起伏界面至30 km处有一低速夹层，速度为5.8 km/s. S波速度由公式v_P/v_S=1.732 给出，并假定其速度结构和界面起伏与P波相同. 模型参数化时，规则立方体单元的边长为5.0 km，不规则网格单元则由界面起伏情形而定，次级节点的间距为0.5 km，地表及两个反射界面也按上述间距离散. 炮点置于： x=50 km，y=50 km，z=-10 km的位置，61 个检波器等间距(5.0 km)分布在模型地表面的3条边上(即地表的左侧、 后侧及右侧的边上).

图 9 一炮多道接收时8种不同类型多次波射线路径.(a)直达波;(b)透射并分别经第一、 二、 三层界面的反射波;(c)地表与第一、 二、 三层界面之间的多次透射、 转换及反射波;(d)经第三层界面的多次透射、 转换及反射波. 图中实心圆代表炮点所在位置，黑色倒三角为地表检波器Fig. 9 The raypaths of 8 different kinds of multiply transmitted，converted reflections on the common source point gathers (a)Direct P wave;(b)transmitted and primary reflections from first，secondary and third interface;(c)multiply transmitted，reflected and converted phases between the top surface and the first interface，between the top surface and the secondary interface， and between the top surface and the third interface;(d)multiply transmitted，converted primary reflection from the third interface. The black circle indicates source location，black triangles denote receivers

首先模拟一次反射的情形，即自炮点激发的P波经第一层起伏界面的反射波(震相P₁P¹)和经第二层起伏界面的反射波(震相P₁P₂P²P¹)，其相应射线路径如图 10a所示.

然后模拟多次透射、 反射及转换波的情形. 震相P¹P₁S¹是炮点激发的P波上行至地表，经反射后下行至第一层反射界面，然后反射且转换成S波上行至检波器. 震相P¹P₁P₂P²S¹是炮点激发的P波上行至地表，经反射后下行至第一层反射界面，透射且转换成S波后继续下行至第二层反射界面，然后反射转换成P波上行至第一层反射界面，最后透射转换成S波上行至检波器. 其相应射线路径见图 10b所示.

图 10(a)透射和一次反射波射线路径．图中黑线为反射震相P1P1，灰线为反射波震相P1P2P2P1 ；(b)多次透射、 转换和反射波的射线路径． 图中黑线为多次反射震相P1P1S1，灰线为多次透射、 转换及反射波震相P1P1S2P2S1 Fig. 10(a)Raypaths of the transmitted and primary reflections. Black lines are the raypaths of P1P1 phase and grey lines are the raypaths of P1P2P2P1 phase；(b)Raypaths of the multiple transmissions，conversions and reflections. Black lines are the raypaths of P1P1S1 phase and grey lines are the raypaths of P1P1S2P2S1 phase

综上所述，用ISPM算法可以实现二维或三维复杂地学模型(含起伏地表、 不规则速度界面、 不规则速度夹层)中任意组合的透射、 反射及转换波的射线追踪，其计算精度要高于有限差分下的FMM算法，且CPU运行时间至少节省5—6倍. 分区多步ISPM算法进行多次波射线追踪的优点如下：

1)采用不规则网格单元划分时可较好地拟合速度模型中的不规则速度界面及起伏地表.

2)算法稳健性高，且具有全局解.

3)计算精度高，运行速度快. 下一步研究的重点则是考虑振幅衰减特征，界面的反射、 透射及转换波系数，进而进行多波偏移成像、 多震相走时和振幅衰减的联合成像.

澳大利亚国立大学地学研究院Nicholas Rawlinson博士为本文提供了分区多步FMM算法程序，在此表示感谢.