引言

射线追踪技术在地震层析成像以及混凝土超声波射线层析成像等领域具有重要作用. 目前常用的射线追踪方法主要有两点射线追踪算法(包括试射法)(Julian，Gubbins，1977； 徐涛等，2004； 田玥，陈晓非，2005)、 弯曲法(Thurber，Ellsworth，1980; Xu et al， 2006，2010)、 有限差分解程函方程法(Vidale，1988; Qin et al，1992; 李振春等，2004)、 最短路径法(Moser，1991; 赵爱华等，2000； Zhao et al，2004; Bai et al，2010; 赵爱华，徐涛，2012)和LTI(linear travel-time interpolation)射线追踪算法(Asakawa，Kawanaka，1993; 赵改善等，1998; Cardarelli，Cerreto，2002; 黄靓，黄政宇，2002; 聂建新，杨慧珠，2003; 张建中等， 2003，2004; 黄靓，2008； 张东等，2009; 卢江波，方志，2014)等． 其中，LTI射线追踪算法因其计算精度较高、 计算速度较快，且适用于任意复杂的速度介质模型，在地震层析成像等领域得到广泛应用． 但原始LTI算法(Asakawa，Kawanaka，1993)所采用的由震源向模型边界外推的计算方式，不能正确追踪逆向传播的射线，影响地震层析成像的精度．

针对这一问题，不少研究人员提出了多种改进算法. 例如： 黄靓和黄政宇(2002)及黄靓(2008)提出了扩张收缩LTI算法，卢江波和方志(2014)提出了扩张收缩LTI改进算法，张东等(2009)提出了循环计算LTI算法； 此外，张建中等(2003，2004)将波前扩展方式与LTI基本方程相结合，提出了动态网络最短路径射线追踪算法. 这些算法中，扩张收缩LTI算法、 扩张收缩LTI改进算法以及循环计算LTI算法都需要迭代计算，对于复杂的速度模型以及网格划分较为精细的模型，其迭代次数较多，计算效率偏低. 动态网络最短路径射线追踪算法能够一次计算出所有节点的走时，具有相对较高的计算效率，且也能有效解决上述问题，但该算法在进行波前扩展时，所选择的插值线段不够合理，无效重复计算较多，计算量大; 此外，该算法采用快速排序法及插入排序法管理波前阵列节点，效率较低．

本文针对动态网络最短路径射线追踪算法存在的不足对其进行改进. 首先依据波的传播规律以及LTI算法的基本方程，选择更为合理的插值线段，并排除该算法中的无效重复计算； 然后采用传统的二叉树堆排序算法对波前阵列节点进行管理，以提高该算法的计算效率．

1.   LTI算法的基本方程

如图 1所示，射线经单元下边界的AB节段到达C节点，射线与AB的交点为D． A点及B点走时分别为t_A和_B，节段AB长为L，单元慢度为s，C点距A点的水平及竖向距离均为已知，分别为x_C，y_C． 建立以A点为原点的局部坐标系，确定A，B，C，D点的局部坐标． 现推导C点走时以及交点D距A点长度r的计算公式(Asakawa，Kawanaka，1993;赵改善等，1998)．

图  1  经过边界AB到达C点的射线路径图

Figure  1.  The diagram of ray path from segment AB to C

下载: 全尺寸图片 幻灯片

由线性追踪算法的基本假设可得D点的走时为

根据D点的走时，结合单元慢度以及各点的局部坐标等条件，可得C点的走时为

将式(1)代入式(2)得

根据费马原理，t_C对r的一阶偏导数应满足等于零的条件，即(设Δt=_B-t_A)

当L²s²-Δt²>0时，解方程(4)可得

若r≥0且r≤L，则

若r < 0或r>L，则计算r=0和r=L时的t_C值，并取两者较小值作为最终t_C值．当L²s²-Δt²≤0时，t_C的计算方法与r < 0或r>L时的情况相同．

2.   LTI算法及动态网络最短路径射线追踪算法存在的不足及改进

LTI算法的逐列和逐行扫描均采用从震源单元所在列或行向模型边界进行递推扫描的方式(Asakawa，Kawanaka，1993)，不能正确追踪如图 2所示存在逆向射线(GR)的路径： 逐列扫描时，射线只能从震源S所在列往右边界传播，而要正确追踪路径SEFGR，射线必须要向左传播； 同样，逐列扫描后再逐行扫描，也不能正确追踪路径SEFGR．

图  2  LTI算法存在的问题

Figure  2.  Problems of LTI algorithm

下载: 全尺寸图片 幻灯片

对于这一问题，张建中等(2003，2004)提出的动态网络最短路径射线追踪算法能予以有效解决. 但是该算法在进行波前扩展时，选择的插值线段不够合理，无效重复计算较多． 现以图 3所示均质模型为例，对该算法的基本步骤及其计算策略进行详细说明，同时依据波的传播规律以及LTI算法的基本方程对该算法进行改进，改进算法的基本步骤及其计算策略如图 4所示. 具体步骤及改进分析如下:

图  3  动态网络最短路径射线追踪算法的基本步骤及其计算策略(a)计算震源所在单元各节点的走时；(n)当前扩展点(例如A₁)的邻点(例如A₁₂)为未计算过情况；(c)当前扩展点(例如A₁)的邻点(例如A₂)为计算过但未做过扩展点的情况；(d)当前扩展点(例如A₁)的邻点(例如犛)为已做过扩展点的情况；(e)所有节点均已做过扩展点的情况

Figure  3.  Basic steps and computing strategy of the path ray tracing algorithm with dynamic networks (a)Compute node's traceltimes of cell at which source is located;(b)Current extension point's(for example the A₁)neighborhood node(for example the point A₁₂)did not compute;(c)Current wxtendion point's(for example the point A₁)neighbrohood node(for example of the A₂)has been computed but was not the extendion point;(d)Current extendion point's(for example the point A₁)neighbrohood node(for example of the S)was extendion point;(e)All nodes were extendion points

下载: 全尺寸图片 幻灯片

1)首先计算震源S所在单元中所有节点的走时，如图 3a及图 4a所示，然后将这些节点加入波前阵列中用于下一步计算．

图  4  改进算法的基本步骤及其计算策略．各子图题同图 3

Figure  4.  Bsaic steps and computing strategy of the improved algorithm presented in thid study. The captions of the subfigures are the same as those in Fig. 3

下载: 全尺寸图片 幻灯片

2)从波前阵列中找出走时最小的节点作为当前扩展点，假定节点为A₁，然后对A₁的所有邻点(S，A₂，A₁₂)一一进行分析，如图 3b--d及图 4b--d所示，并根据邻点的不同状态采用不同的计算策略，具体如下：

① 邻点未计算过的策略. 当前扩展点不与该邻点形成插值段，而是以当前扩展点为震源，对单元内还未做过扩展点的节点进行计算． 对邻点A₁₂进行分析时，即采用这一策略，如图 3b所示． 事实上，对于邻点未计算过的情况，只需利用当前扩展点A₁的走时信息对该邻点进行计算(图 4b). 考虑波沿单元边界的速度与通过单元内部的速度不同，将Ⅱ号单元中的节点分为A₁A₄、 A₁A₁₄边界上的节点以及其它节点； 对于其它节点，显然在A₁₂做扩展点，并与A₁形成插值段时进行计算更为合理有效； 对于A₁A₄、 A₁A₁₄边界上的节点，与最短路径法的计算方法相同，只需计算节点A₂以及A₁₂，而节点A₂的计算则可放在对邻点A₂的分析中，因此，此时只需对节点A₁₂进行计算．

② 邻点已计算过但未做过扩展点(不确定该邻点是否已得到理论最小走时)的策略. 当前扩展点与该邻点形成走时插值段，并应用LTI基本方程计算插值段所在单元中还未做过扩展点的节点． 对邻点A₂进行分析时，即采用这一策略，具体如图 3c所示． 事实上，对于邻点已计算过但未做过扩展点的情况，只需利用当前扩展点A₁的走时信息对该邻点进行计算(图 4c). 由LTI基本方程式(3)可知，在邻点A₂计算过但未做过扩展点的情况下，节点A₁与A₂形成插值段的意义不明确． 而且节点A₂做扩展点时，节点A₂与A₁也会形成插值段，此时，节点A₁、 A₂均确定得到最小走时． 显然，与未做过扩展点的邻点A₂形成插值段既不合理也无必要． 但考虑节点A₂的理论射线路径可能经过节点A₁(如Ⅱ号单元的速度远大于Ⅰ号单元)，因此，在不形成插值段的情况下，还需利用节点A₁的走时信息对节点A₂进行计算．

③ 邻点已做过扩展点(可以确定该邻点已得到理论最小走时)的策略. 当前扩展点与该邻点形成走时插值段，并应用LTI基本方程计算插值段所在单元中还未做过扩展点的节点． 对邻点S进行分析时，即采用这一策略，具体如图 3d所示． 实际上，在当前扩展点与已做过扩展点的邻点形成插值段时，并不需要对插值段所在单元中所有未做过扩展点的节点进行计算，此时仅需考虑位于当前扩展点一侧且不与插值段在同一边界的节点(图 4d). 而且在对这些节点进行计算时，还可以通过规定节点的计算顺序，并增加一些判断条件，进一步排除无效重复计算. 具体为： 先计算位于插值段A₁S对边A₄A₇，且与当前扩展点A₁处于同一位置的节点A₄，然后再依次计算其它节点A₃、 A₂； 在节点计算前对节点进行判断，以节点A₄为例，若节点已获得的走时小于当前扩展点邻点的走时与节点至当前扩展点的走时之和，即t_A4 < t_S+l_A1A4/v_I(t_A4和t_S分别为节点A₄及节点S已获得的节点走时，l_A1A4为节点A₁与A₄的距离，v_I为Ⅰ号单元的速度)，则该节点可不计算； 在节点计算后，以节点A₄为例，若节点满足费马原理的射线过当前扩展点A₁，则该节点后面的节点A₃、 A₂可不再计算． 下面对节点计算范围以及计算前、 后的判断条件进行说明：

a)计算范围的确定. 首先，对于与插值段处于同一边界的节点，显然只需考虑与当前扩展点相邻的节点，当邻点已做过扩展点时(如图 3d中的节点S)，不需要计算，若邻点未做过扩展点，则可按照上面的计算策略①或②进行计算，此时无需考虑； 然后，对位于邻点一侧且不与插值段处于同一边界的节点(A₅—A₉)，可知在节点S已做过扩展点的情况下，若A₅—A₉节点的理论射线路径过A₁S插值段，则交点必然为节点S． 现对节点S的左邻点A₁₁进行分析，若节点A₁₁的理论最小走时小于或等于节点S的理论最小走时，则节点A₆—A₉的理论射线路径不过A₁S插值段，因为这些节点通过A₁₁计算得到的走时较通过A₁S插值段更小； 若A₁₁的理论最小走时大于节点S的理论最小走时(图 3d)，则当A₁₁做扩展点时，节点S必已做过扩展点，A₁₁将与节点S形成插值段，通过对A₁₁一侧的节点进行计算，节点A₆—A₉能得到最小走时. 无论哪种情况均不用利用A₁S插值段对节点A₆—A₉进行计算． 对于节点A₅，可放在节点S做扩展点时进行计算，考虑节点S做扩展点时，节点A₁₁和A₁可能均未做过扩展点，而节点A₁₁和A₁做扩展点时，又只考虑位于节点A₁₁或A₁一侧的节点，A₅将不能得到最小走时，因此在节点S为当前扩展点、 而其邻点A₁₁和A₁均未做过扩展点的情况下，除对节点A₁₁和A₁进行计算外，还需增加节点A₅的计算(具体见计算策略④)．

b)节点计算前的判断条件. 以图 4d中的节点A₄为例进行说明. 显然，邻点S的走时为插值段A₁S上的最小走时，节点A₄到插值段A₁S的最短距离为A₁A₄. 根据LTI基本方程式(1)和(2)可知，若t_A4 < t_S+l_A1A4/v_I，则节点A₄的理论射线路径必然不过A₁S插值段，可不计算． 需要说明的是，在A₁为当前扩展点、 S为其邻点的情况下，节点A₂—A₄不满足这一判断条件，都需要计算．

c)节点计算后的判断条件. 假定A₄满足费马原理的射线过当前扩展点A₁. 显然，节点A₃和A₂满足费马原理的射线也过当前扩展点A₁. 如果A₂—A₄的理论射线路径过A₁S插值段，那么这些节点的理论射线路径必然沿A₁A₄边界传播，与最短路径法的计算方法相同，只需对A₁的邻点A₂进行计算即可. 而节点A₂的计算则可放在对邻点A₂的分析中. 因此，当A₄满足费马原理的射线过当前扩展点A₁时，A₄后面的节点A₃和A₂均不用计算．

以上分析针对的是当前扩展点为单元端节点的情况，对于当前扩展点为中间节点的情况亦有同样的结论. 以节点A₁₁为当前扩展点、 S为其已做过扩展点的邻点为例，节点的计算顺序从A₆到A₉，若某节点满足费马原理的射线过当前扩展点A₁₁，假定为A₆，则根据LTI基本方程可以证明，A₇—A₉满足费马原理的射线也过当前扩展点A₁₁． 如果A₇—A₉的理论射线路径过A₁₁S插值段，那么这些节点的理论射线路径必然过节点A₁₁． 现对节点A₁₁的左邻点A₁₀进行分析. 如果节点A₁₀的理论最小走时小于或等于节点A₁₁的理论最小走时，则节点A₇—A₉的理论射线路径必然不过A₁₁S插值段，因为这些节点通过A₁₀计算得到的走时较通过A₁₁计算的走时更小； 如果节点A₁₀的理论最小走时大于节点A₁₁的理论最小走时，那么当节点A₁₀做扩展点时，节点A₁₁必然已做过扩展点，A₁₀与A₁₁将形成插值段，通过对A₁₀一侧的节点进行计算，A₇—A₉可得到最小走时； 无论哪种情况，A₆后面的节点A₇—A₉均不用计算．

④ 如果当前扩展点为单元的中间节点，且其两个邻点均未做过扩展点，则需利用当前扩展点的走时信息，对位于当前扩展点对边，且与当前扩展点处于同一位置的节点进行计算． 以图 4d为例，如果I号单元右边界上的中间节点A₂为当前扩展点，且其两个邻点A₁及A₃均未做过扩展点，则需利用A₂的走时信息对A₉及A₁₅进行计算． 需要特别说明的是，这一条计算规则是改进算法特有的. 动态网络最短路径射线追踪算法没有这一规则．

对于以上4种计算策略，所有节点在计算后，都需进行如下处理： 若节点之前未被计算过，则更新其走时并将其加入波前阵列中； 若节点之前已计算过，而新计算出的走时小于原来的走时，则更新其走时． 在分析完当前扩展点的所有邻点后，将当前扩展点从波前阵列中删除．

3)重复步骤2)，直到模型中所有节点均做过扩展点，此时所有节点均得到最小走时，图 3e、 图 4e为最终状态．

对比图 3与图 4可以看出，动态网络最短路径射线追踪算法存在较多的无效重复计算，经改进后，该算法的节点计算量明显减少． 此外，动态网络最短路径射线追踪算法采用快速排序法和插入排序法管理波前阵列节点，并认为这种排序方法要优于二叉树堆排序算法，且两者的运行时间分别为O(n)和O(nlgn)，式中n表示当前波前阵列的节点数(王辉，常旭，2000； 张建中等张建中等(2003，2004)在引用王辉和常旭(2000)一文时，认为n为总节点数. 实际上是误解了原文的意思.，2003，2004). 但是，这一结论的前提是波前阵列节点始终保持从小到大的排列顺序． 事实上，采用堆排序算法对波前阵列进行管理并不要求波前阵列节点保持这种顺序，而仅要求堆结构在插入、 更新以及删除节点后仍保持最小堆的性质即可，整个过程的运行时间为O(nlgn)(Cormen et al，2009)． 因此二叉树堆排序算法实际上要优于插入排序算法． 当然，除二叉树堆排序算法外，还存在性能更为优越的算法，如Fibonacci堆、 quake堆以及Brodal队列等(Brodal，2013)，但这些算法较为复杂，实用性还存在不足． 综合考虑，本文仍采用传统的二叉树堆排序算法对波前阵列节点进行管理．

3.   数值算例

为说明原始LTI算法存在的问题，比较各算法的计算效率，以及验证动态网络最短路径射线追踪算法和本文改进算法的有效性，建立了如图 5所示的存在高速区和低速区的连续介质模型，其整体尺寸为2500 m×600 m. 其中深度在400—600 m之间的区域为高速区，速度为6000 m/s； 深度在200—400 m之间且水平方向在500—2200 m的区域为低速区，速度为500 m/s； 其余区域的速度随深度增加而线性增加，深度z处的速度v(z)=v₀(1+βz)，式中v₀取2 000 m/s，β取0.001． 震源S位于模型的左上角，接收点R均匀布设于模型上边界，间距为100 m； 单元大小为5 m×5 m，共有6×10⁴个单元，单元边界划分段数为1．

图  5  模型速度

Figure  5.  Velocity model

下载: 全尺寸图片 幻灯片

单元内波传播速度按以下方式选取： 高速区及低速区单元波传播速度分别取6000 m/s和500 m/s，其它区域单元波传播速度按波沿单元竖向边界传播耗时相同的原则确定v(i)= 10/ln[(200+i)/(199+i)]，式中v(i)表示第i行单元的速度． 数值模型接收点的理论射线路径、 各算法的射线追踪结果以及各算法初至波走时的误差分别如图 6—9所示. 其中，图 6的理论射线路径可通过速度随深度线性增加情况下地震波的射线方程式(陆基孟，王永刚，2009)分析计算得到； 图 8中相对误差的计算公式为(t₁－t₂)/t₁ ×100%，式中，t₁为单元边界划分30段并采用最短路径法计算的节点走时，t₂为各算法的计算结果； 图 9中绝对误差的计算公式为t₂－t₁．

图  6  各接收点的理论射线路径

Figure  6.  Theoretical ray paths of all receiving points

下载: 全尺寸图片 幻灯片

图  7  原始LTI算法(a)、 动态网络最短路径射线追踪算法(b)和本文改进算法(c)的射线追踪结果

Figure  7.  Ray path tracing results by the three algorithms(a)Original LTI algorithm;(b)The shortest path ray tracing algorithm with dynamic networks;(c)Improved algorithm presented in this study

下载: 全尺寸图片 幻灯片

图  8  原始LTI算法(a)、 动态网络最短路径射线追踪算法(b)和本文改进算法(c)的走时相对误差a图中x < 500 m的区域没有出现如b，c图所示的等高线，是因为该区域内的相对误差所对应的颜色为浅色，未能有效显示所致. 图 9也存在类似的情况

Figure  8.  Travel time relative errors by three algorithms(a)Original LTI algorithm;(b)The shortest path ray tracing algorithm with dynamic networks;(c)Improved algorithm presented in this study. The reason why the contour in the x < 500 m do not present in Fig. 8a like Figs.8b or 8c is that the color which corresponds to the value of relative error in the x < 500 m was too light to present in Fig. 8a. The similar situation is presented in Fig. 9

下载: 全尺寸图片 幻灯片

图  9  原始LTI算法(a)、 动态网络最短路径射线追踪算法(b)和改进算法(c)的走时绝对误差

Figure  9.  Travel time absolute error by three algorithms (a)Original LTI algorithm;(b)The shortest path ray tracing algorithm with dynamic networks;(c)Improved algorithm presented in this study

下载: 全尺寸图片 幻灯片

对比图 7a与图 6可以看出，原始LTI算法对于大多数接收点均能正确追踪其射线路径，但当接收点的理论射线路径存在逆向传播的射线时(如图 6中的c，d，e，f点)，原始LTI算法不能对其进行正确的射线追踪． 图 7b，c中各接收点的射线路径与图 6中的理论射线路径一致，表明动态网络射线追踪算法和本文改进算法均能考虑逆向传播的射线，能够正确处理复杂介质模型中接收点的射线追踪问题．

从图 8和图 9可以看出，原始LTI算法在存在逆向传播射线的区域，其走时相对误差和绝对误差较大，最大分别达到43.36%和297.8 ms； 而动态网络最短路径射线追踪算法和本文改进算法则不存在这一问题，均能够正确计算各节点的最小走时，其最大相对误差和绝对误差分别为3.96%和1.5 ms．

为对比3种算法的计算效率，以图 5所示模型为计算对象，单元边界划分段数为1，2，4，10，20和30等6种情况，分别记录各算法的计算时间(不包括算法的前处理过程). 其中，动态网络最短路径射线追踪算法采用插入排序法及二叉树堆排序法对波前节点进行管理，本文改进算法只采用二叉树堆排序法． 计算机CPU主频为4.3 GHz，对比结果如表 1所示.

表  1  3种算法计算效率的对比

Table  1.  Comparison of computational efficiency for three algorithms

下载: 导出CSV 

| 显示表格

从表 1可以看出，本文改进算法具有较高的计算效率，是动态网络最短路径射线追踪算法的4.5—30倍，是原始LTI算法的2—6.5倍； 当动态网络最短路径射线追踪算法采用堆排序算法时，本文改进算法的计算效率是其3.5—15倍． 为验证本文改进算法在更为复杂模型上的有效性，下面采用本文改进算法对Marmousi速度模型进行计算，该模型尺寸为9192 m×2904 m，震源位于模型的上表面(4800 m，0)，计算时采用的单元尺寸为24 m×24 m，单元边界划分段数为4，计算结果如图 10所示. 其中图 10b为本文改进算法的相对误差，其计算公式为(t₁-t₂)/t₁×100%，式中，t₁为单元边界划分30段并采用最短路径法计算的节点走时，t₂为改进算法计算的走时． 计算结果验证了本文改进算法对复杂速度模型的有效性．

图  10  本文改进算法对Marmousi速度模型的初至波计算结果(a)初至波等时线;(b)相对误差

Figure  10.  The first arrival calculations result based on Marmousi velocity model by using the improved algorithm presented in this study (a)Isochrons of first arrivals;(b)The relative error

下载: 全尺寸图片 幻灯片

4.   讨论与结论

动态网络最短路径射线追踪算法能够有效解决原始LTI算法不能正确追踪逆向传播射线的问题，且只需一次扩展计算就能得到所有节点的最小走时，具有较高的计算效率. 但该算法存在较多的无效重复计算，且其波前阵列节点的管理方法效率较低. 为此本文根据波的走时信息以及LTI基本方程提出了改进动态网络最短路径射线追踪算法，排除了动态网络最短路径射线追踪算法中大量的节点计算，同时采用传统二叉树堆排序法管理波前阵列节点，较大幅度提高了该算法的计算效率. 数值算例表明，本文提出的改进算法具有较高的计算效率，是动态网络最短路径射线追踪算法的4.5—30倍，是原始LTI算法的2—6.5倍； 当动态网络最短路径射线追踪算法采用堆排序算法时，本文改进算法的计算效率是其3.5—15倍． 此外，若将本文的改进算法用于三维模型，则效果将更为显著.