利用最短路径算法进行地震波射线追踪

杨春雨　王小磊　李栋

射线追踪方法作为一种快速有效的波场近似计算方法，不仅对于地震波理论研究具有重要意义，而且也直接应用于地震波反演及偏移成像等过程。该文在收集、整理国内外有关研究资料的基础上，着重介绍了最短路径算法，并通过对Marmousi模型的计算，验证了最短路径方法的高精度，及对复杂模型的适应性。

射线追踪的方法种类较多。传统方法有基于初值问题的打靶法（Shooting method， i.e.， Julian and Gubbins，1977）和基于边值问题的弯曲法（Bending method， i.e.， Um and Thurber， 1987），但是他们不能处理介质中较强的速度变化，有时无法求出全局最小走时，计算效率较低，阴影区内无射线。随着射线追踪方法的发展，出现了直接从Huygens原理或Fermat原理出发，采用等价波前描述地震波场特征的方法。Vidale（1988，1990）和Podvin（1990）等人则从程函方程出发，首先求出走时场分布，再计算走时场的最速下降方向的办法，得到每一条接收点到震源的射线路径。随后，Qin（1992）等人对Vidale的方法作了改进，提出了波前扩展方法。黃联捷、李幼铭、吴如山（1992）基于Huygens 原理提出WFRT方法，根据计算精度的要求，逐次细化网格，而求得最小走时。Sava 和Fomel（1998）提出了HWT（huygencs wavefront tracing）法。Moser（1991）提出了根据费马原理的最短路径法。

自上世纪八十年代后期发展起来的基于网格单元的射线追踪算法由于其诸多优点，因而倍受人们的青睐。与传统的射线追踪算法相比而言，基于网格单元的算法具有四大优点：①可利用波振面向外扩展传播的原理一次性计算出速度模型中所有网格节点的射线走时及相应的路径，并能正确的追踪检波器位于射线阴影区的衍射波射线；②算法数值计算稳定，所得到的解总是全局最佳射线路径及相应最小走时；③多炮多道接收时表现得更为高效，精度也比传统射线追踪方法要高；④在连续介质中网格单元算法始终能找到初至波走时，而传统射线方法则只能找出唯一走时，且很难判别这个走时是属于初至波还是属于后续波（De Kool et al，2006）。目前，基于网格单元的射线追踪算法已成功用于地震定位、地震偏移成像和地震层析成像中。

在众多基于网格单元的射线追踪算法中，最短路径射线追踪算法，简称SPM以其精度高、计算稳健以及较强的模型适应性而广受推崇。

2 最短路径算法

最短路径算法是在网格化节点速度场中计算所有节点走时的一种十分有效的射线追踪方法。这种方法是由网格节点的连线作为具有走时的射线段，采用类似Dijkstra算法来求取给定点到网格上所有节点的最短路径。根据费马原理两点间最短路径（最小走时路径）相应于实际射线路径。

最短路径网格通常是由单元或中心节点来参数化。Nakanishi和Yamaguchi （1986）把速度场参数化为由常速度单元组成，而节点定义在单元边界上。该算法的优势在于每对节点间的走时可以很容易的计算 ，这里d是两节点间的距离，s是含节点单元的波慢度。计算精度可通过缩小单元尺度或增加单元边界上节点的个数来提高。

另一种形成网格节点的方法是利用规则的速度节点，并进行线性连接（Moser，1991）。两连接节点A和B间的走时可简单地由： 得到，这里 和 分别为节点A和B处的波慢度值，而d则为两节点间的距离。

一旦网格节点结构及节点间走时计算的方式确定之后，接下来就是计算整个走时场和相应射线路径。Dijkstra （1959）最早提出原始的网格理论算法，该文原载于中国社会科学院文献信息中心主办的《环球市场信息导报》杂志http：//www.ems86.com总第547期2014年第15期-----转载须注名来源其中计算时间正比于 ，而 为所有节点总数。此算法的概念十分简单，即：有总数为Q的未知走时节点，起初 Q 含M个元素，而 P 是空集，将Q集节点的走时设置为任意大的数。算法将炮点所在节点加入P集开始，然后计算临近节点的走时。这些组成了可能的走时，然后算法从中挑选最小走时，将其加入 P 集直至所有 Q 集内节点计算完毕为止。如果 Q 集内的节点在上次循环中已有计算走时，则选择具有最小走时使之更新。完整的走时场可通过M次迭代得到，射线路径可通过记录节点更新的顺序获得。

3.模型实例

为了检验最短路径算法的有效性、精确性及对复杂模型的适应性，我们选择了较为复杂的marmousi模型。将模型划分为三个计算区域。炮点位于模型左上角（0m，0m）处，910个检波器等间距（100m）的分布在地表，最大炮检距为9100m。我们分别使用波场模拟方法和最短路径算法计算来自第一、二个不规则界面的一次反射波，并计算两种算法计算走时的相对误差Ere。根据计算结果，两种算法的相对误差在0.2%之内，而且，最短路径计算走时均小于波场模拟计算走时，更接近真实解，并且走时差大小与反射界面的起伏呈很好的相关性。说明最短路径算法能够更好的模拟界面的起伏特征，保证计算的精度。

4 结论

通过使用最短路径算法对marmousi模型的试算，我们可以得出以下结论：

最短路径算法是一种十分有效的射线追踪算法，具有精度高，计算稳健的特点；最短路径算法能够适应十分复杂的地质模型，具有较强的实用性；计算过程中最短路径算法需要较多内存空间，计算耗时较长，需要改进。

（作者单位：1.长安大学地质工程与测绘学院；2.长安大学地球科学与国土资源学院）