地震波初至走时的计算方法综述 赵烽帆 ¹ , 马婷 ^2,3 , 徐涛 ²

1. 中国地震台网中心, 北京 100045;
2. 中国科学院地质与地球物理研究所, 岩石圈演化国家重点实验室, 北京 100029;
3. 中国科学院大学, 北京 100039

收稿日期: 2013-9-9 修回日期: 2014-1-11 .

基金项目：中国地震局地震科技星火计划青年项目（XH14055Y）、中国地震局监测预报司震情跟踪合同制定向工作任务（2014020412）和国家自然科学基金（41174075）联合资助.

作者简介：赵烽帆，女，1983年生，工程师，主要从事震源物理及地震活动性研究.（E-mail：[email protected]）

通讯作者：徐涛，男，1978年生，副研究员，主要从事地震波传播和壳幔结构成像研究.（E-mail：[email protected]）

摘要 ：在地震波场中，初至波到时信息由于初至震相可追踪、易识别性，在地震学领域占有重要的位置，广泛地应用于叠前偏移、叠前速度分析、地震走时层析成像及地震定位等.本文主要介绍了四类具有代表性的计算初至波走时的方法：（1）基于高频近似射线理论方法，如最短路径方法（SPM），及修正后的最短路径方法（MSPM）；（2）基于程函方程的数值解方法，如有限差分方法（FD）、快速推进法（FMM）和快速扫描法（FSM）；（3）基于惠更斯原理的波前构建法（WFC）；（4）基于频率域波动方程数值解法（FWQ）.最短路径方法计算精度较高，稳定性较好，但其需要采用更多的网格节点，因此计算效率低；程函方程数值解法无需计算射线路径，具有计算效率高、稳定性较好、易于实现等优势，但其计算精度较低，可以通过引入高阶差分格式得到提高；波前构建法计算精度高，稳定性好，但其需要在射线网格和规则网格之间做网格转换，因此计算效率较低；频率域波动方程方法能适应任意复杂介质，但其计算精度和计算效率较低. 关键词 ： 最短路径法 快速推进法 快速扫描法 波前构建法 频率域波动方程法 A review of the travel-time calculation methods of seismic first break ZHAO Feng-fan ¹ , MA Ting ^2,3 , XU Tao ²

1. China Earthquake Networks Center, Beijing 100045, China;
2. State Key Laboratory of Lithospheric Evolution, Institute of Geology and Geophysics, Chinese Academy of Sciences, Beijing 100029, China;
3. University of Chinese Academy of Sciences, Beijing 100039, China

Abstract : In the seismic wave field, the first arrival travel-time plays an important role in the field of seismology, which is due to the first arrival phases can be easily traced and identified. The seismic first arrival travel-time field is widely used in seismic prestack migration, velocity analysis, seismic travel-time tomography and earthquake location. This paper mainly introduces four typical seismic first arrival travel-time calculation methods, that is, (1) seismic ray method based on the high-frequency approximation, the commonly used methods are Shortest Path Method (SPM) and Modified Shortest Path Method (MSPM); (2) the method based on the numerical solution of eikonal equation, which mainly include Finite Difference Method (FD), Fast Marching Method (FMM) and Fast Sweeping Method (FSM); (3) Wavefront Construction (WFC) method based on the Huygens principle; and (4) the method based on the frequency domain wave equation (FWQ). The first one has higher computational accuracy and better stability, while it needs more grid points, which may result in low computational efficiency. The calculation of ray paths is not required in the second one, so it has advantages in computational efficiency, stability and realization, but the computational accuracy is lower, which can be improved by introducing the high-order finite difference scheme. The third one can calculate traveltime with high accuracy and stability, while it requires the grid transition between ray grid and regular grid, which may result in low computational efficiency. The last one can adapt any complex medium, but the computational accuracy and efficiency are lower. Key words : first break shortest path eikonal equation finite difference fast marching fast sweeping wavefront construction frequency domain wave equation  1.2 修正后的最短路径方法

Bai等(2007) 提出了三维改进后的最短路径算法(Modified SPM，MSPM)，在保证计算精度的前提下大大提高了计算速度.MSPM算法与传统的SPM在确保射线角度的覆盖率方面有着本质的区别.为了增加射线角度的覆盖，SPM在波前点的计算时需要向外扩大大量的网格点(Forward Star技术)，计算效率不高，同时离波前点很远的网格点的走时只考虑了两点的速度；为了更好的平衡射线计算精度与计算时间的矛盾， Bai等(2007) 提出的MSPM，通过在单元插入次级节点来有效的解决射线角度覆盖问题，同时这种方法有效减少单元个数，适应大速度模型.次级节点、炮点和检波点(如不与网格重合)的速度值通过拉格朗日三次线性插值得到，可更好的适应复杂介质.传统的SPM和MSPM的网格的角度覆盖情况(以二维为例，见 图 2 )有很大区别，前者是通过向外扩大网格点(Forward Star技术)，来增加射线的角度覆盖；后者是通过插入很多个次级节点，大大增加了射线的角度覆盖以获取最小走时.相比于SPM方法，MSPM方法能适应大网格的速度模型，有更快的计算速度和较高的计算精度.具体计算过程如下：在整个计算过程中，把所有网格节点的集合 N分为四个子集.P：已获得最小走时的节点，即已作过次级震源的节点集合；Q：以计算出从P中至少一个节点传来的走时，还没有作为次级震源节点的集合；R：在N中除去P、Q后剩余的节点集合；F(i)：与震源或次震源直接相连的节点集合.具体计算算法如下：(1)初始化.P=Φ，t(s)=0，s为震源，Φ为空集；Q=N，t(i)=∞，i∈N .(2)选择.在Q中选择走时最小的节点i，i∈Q .(3)更替.计算从i点传到j点的旅行时，若该值比原值小，则用该值取代原值，否则保持不变.即，对j∈F(i)∩Q，t(j)= min [t(j)，t(i)+d _ij ]；对j∈F(i)∩R，t(j)=t(i)+ d _ij .将节点i从Q转至P.迭代判断.如果P=N或Q=Φ停止迭代，否则转向(2).确定从震源点到某一点的波传播射线路径，是从该点开始向震源逆向搜寻处理的过程.在从震源点开始扩展波前的过程中，存储下各个最小走时节点的前一次级震源点.确定射线路径时，从接收节点开始，依次找出前一级次级震源点的位置，直到震源点结束. 图 3 显示了该方法对 Marmousi模型的一个应用实例(Bai， et al.， 2007).该模型包含了许多复杂的特征，例如盐丘构造，断层，薄层，以及强弯曲界面( 图 3 a ).炮点位于模型底部(x=6000 m ，z=-2800 m )处，384个检波器的排列(间隔24 m )位于模型顶部. 图 2 传统的SPM网格与改进后的MSPM网格的射线角度覆盖情况 (a)传统Forward Star技术；(b)MSPM网格，主节点为波前点；(c)MSPM网格，次级节点为波前点；其中实心圆为主节点，空心圆为次级节点. Fig. 2 Ray angle coverage for traditional SPM grid and MSPM grid (a)Traditional Forward Star technique;(b)MSPM grid，primary nodes are wavefront points;(c)MSPM grid，secondary nodes are wavefront points. Solid circles denote primary nodes，hollow circles denote secondary nodes. 2 基于程函方程的数值解方法

2.1 经典有限差分法

有限差分法( Finite Difference Method，FD) 最早由Vidale在1988年提出，并于1990年将该方法推广到三维.该算法采用有限差分离散程函方程，并结合一定的波前扩展方式来计算模型的初至走时.因此，FD最核心的两个问题是：(1)采用怎样的差分格式离散Eikonal方程；(2)采用怎样的波前扩展方式.关于这两个问题，前者主要影响算法的计算精度和稳定性，后者主要影响算法的计算效率和合理性.

二维程函方程为

如 图 4 a 所示，设地震波到达A、B ₁ 、B ₂ 点的走时分别为t ₀ 、t ₁ 、t ₂ ，方程中的两个微分项可用有限差分近似表示为

式中h表示网格单元边长。将以上两式带入方程(1)，可以得到 C ₁ 点的走时( 图 4 a )为

该式为平面波外推公式.

图 3 2D Marmousi模型应用实例( Bai et al ., 2007 ) 上图：Marmousi模型；下图：沿顶部72 m 间隔检波器排列的射线路径. Fig. 3 Example for 2D Marmousi model ( Bai et al ., 2007 ) Upper diagram: Marmousi model; bottom diagram: corresponding selected raypaths for receivers at 72 m spacing along top surface.

为保证波前面曲率较大时的走时计算精度， Vidale 又提出了一种球面波外推公式.取A点为坐标原点，设波阵面曲率中心的坐标及走时分别为(－x _s ，－z _s )和t _s ，则点A，B ₁ ，B ₂ 和C ₁ 的走时分别为

上式为球面波外推公式.计算中可单独使用(3)式或(4)式进行走时外推，或者将二者结合使用. 图 4 有限差分法算法示意图( Vidale，1988 ) (a)网格节点A及周围的8个点；(b)外推计算二维节点走时示意图；(c)节点走时计算的顺序. Fig. 4 Schematic diagram of finite difference algorithm( Vidale，1988 ) (a)The grid point A and the eight points in the ring surrounding point A;(b)Schematic diagram of extrapolation of 2D node travel-time;(c)Sequence of calculating the node travel-time. 图 5 高速半空间上存在低速层模型的有限差分法走时结果( Vidale，1988 ) 波前位置间隔3-sec，从震源位置A到检波器B1，B2，B3的最小走时射线路径也显示在图中 Fig. 5 FD results for a high-velocity half-space below a low-velocity layer( Vidale，1988 ) Location of the wavefront is at 3 s intervals. Minimum travel-time paths from a source at A to receivers at B1，B2， and B3 are also show. 在算法的整体实现方面，Vidale采用的是波前扩展方式，如 图 4b 所示.在环内的点的走时为已知，环外的为未知.按四角的点顺序地求出四面的解，这里从右边开始来求解 图 4b 中矩形框中的点的走时. 图 4c 所示为一个环的一边的求解例子.在该行上的点按从左至右的顺序检查并确定相对最小点.相对最小即是对相邻点的时间的相对最小.为确定边缘上的第一个点的时间，必须使用公式(1)的非中心的有限差分.应用平面波的公式为

其中，t ₃ 为未知时间，t ₀ 为在内侧行的相对最小时间，t ₁ 和t ₂ 为时间等于t ₀ 点两侧点的时间.在每个相对最小的点开始沿每行求解旅行时间直到遇到相对最大的值.完成整行的从左到右的扫描后，再按相反的方向扫描，然后剩余的点从相对最小到相对最大的顺序求解.相对最小是指从两边求解出两个旅行时取其中最小.一旦四条边的时间以此方式求解出，则四角的时间也可得出.随后就可以计算下一环.重复上述方法，整个二维节点上的旅行时就都可计算出来了.

有限差分法计算走时可用于地震定位， Kirchhoff 偏移，走时层析成像等. 图 5 显示了高速半空间(8.0 km/s)上存在一个30 km厚的低速层(6.0 km/s)模型的走时结果.其可以追踪到折射波的射线路径.

其中, 为空间步长，且

在计算区域的边界处采用一阶差分近似.例如，在左侧边界点t _1，j ，在x方向上用一阶差分，程函方程可离散为

图 6 描述了用快速推进法( Sethian，1999 ; Sethian and Popovici， 1999 ; 兰海强等，2012b )求解离散的程函方程的详细过程.首先初始化震源点，其走时在后续计算中保持不变.我们把震源点周围的点的走时按从小到大的顺序存在一个数组里，那么它就构成了这一时刻的波前.取出其中的最小值，固定它的值(在后续迭代中保持不变)，所有周围未被固定的点的值可按下式进行迭代：

将每一个新计算出来的网格点都放进波前数组里，利用堆选排技术保存旅行时，然后再从堆栈里取出最小值，更新周围未被固定的点的走时，依次重复，直至遍布整个计算空间.

图 6 快速推进法更新步骤( 兰海强等，2012b ) 网格点可以被分为：活动点、试验点和远点.活动点(黑色)是走时已经知道的点.试验点(灰色)是活动点周围的点，这些点上的走时将被计算.试验点的集合被称作窄带.随着窄带的推进，窄带内的点逐渐被更新为活动点.远点表示走时还未计算的点.当走时传播时，远点逐渐变为试验点.图(a~f)显示了走时传播的顺序：(a)黑色的点表示初始的震源点；(b)计算黑点周围的点的走时，这些点由远点变为试验点；(c)取出试验点中走时最小的点(假设为‘A’)；(d)计算A点周围的点的走时，这些点由远点变为试验点；(e)取出试验点中走时最小的点(假设为‘D’)；(f)计算D点周围的点的走时，这些点由远点变为试验点.如此循环直至遍布整个计算空间. Fig. 6 Updating procedure for the fast marching method( Lan et al ., 2012b ) Grid points are classified as ‘alive’，‘trial’ and ‘far’. ‘Alive’ points(in black)are those where the values of T are known. ‘Trial’ points(in grey)are those around the curve(of ‘alive’ points)，at which the propagation is to be computed. The set of trial points is called the narrow b and . To compute the propagation，points in the narrow b and are updated to ‘alive’，while the narrow b and advances. ‘Far’ points(in white)are those at which the propagation has not yet been computed. During the propagation，the far points are converted into trial points. Figures(a~f)provide an explanation of this sequence: in(a)the black point(alive)represents the initial point; in(b)the value of T is computed in the neighborhood of the black point; the points in this neighborhood are converted from far(white)to trial(gray)points; in(c)the trial point with the smallest value of T is chosen(for example ‘A’); in(d)values of T are computed in the neighborhood of point A，converting them from far to trial; In(e)the trial point is chosen that has the smallest value of T(for example，‘D’); in(f)points in the neighborhood of D are converted from far to trial， and the procedure continues in this manner. 2.3 快速扫描法

快速扫描法( Zhao， 2005 ， 2007 ; 兰海强等，2012b )的主要思想是基于因果关系将走时场传播的方向分成有限个组，对于每一组分别利用Gauss-Seidel迭代方法求解非线性逆风差分格式离散化后的方程组.每一次Gauss-Seidel迭代也称为一次扫描，每次扫描(Gauss-Seidel迭代)按一定的方向求解沿该方向传播的走时场，这里以二维为例.

图 7 二维程函方程中的快速扫描方法( 兰海强等，2012b ) (a)二维程函方程中的快速扫描方法.两条虚线分别表示 x、y 轴，震源位于原点，每个象限中的数字分别表示适合计算该象限走时场的扫描顺序;(b)第一象限走时场的快速扫描法计算.虚线2上的点的走时只依赖与虚线1上的点的走时；坐标轴上的点的走时与坐标轴外点的走时无关(如 x 轴正半轴上的点的走时只依赖于它左侧点的走时，y轴正半轴上的点的走时只依赖于它下方点的走时)；而4个象限被x、y轴隔开，因此在扫描过程中在某一象限传播的走时不会穿过x、y轴而在其它象限传播，即在扫描过程中某个象限内点的走时只依赖于本象限内点的走时，与其它象限点的走时无关. Fig. 7 The fast sweeping algorithm for the eikonal equation in two dimensions.( Lan et al .，2012b ) (a)The fast sweeping algorithm for the eikonal equation in two dimensions. The two dotted lines denote x and y axis， respectively; the source is located at the origin，the numbers in each quadrant mean the appropriate sweep order for this quadrant;(b)The traveltime calculation with the fast sweeping method in the first quadrant. Traveltimes at grid points on the dashed line two are determined by values at grid points on dashed line one，etc. Moreover the four quadrants are separated by two grid lines，i.e.，the x and y axes. The traveltimes of grid points on these two lines do not depend on any values of grid points off these two lines. So the propagation of traveltimes in one quadrant during the corresponding sweep cannot cross these two lines into other quadrants.

二维情形走时场传播的所有方向可分为四组：右上、左上、左下和右下.初值在这四个方向上依次传播.我们可以通过四组不同顺序的Gauss-Seidel迭代，即(1)i=1:I，j=1:J，(2)i=I:1，j=1:J，(3)i=I:1，j=J:1，(4)i=1:I，j=J:1，并结合逆风差分格式分别求解沿这四组方向传播的走时. 图 7 说明了对于一个独立初始点的走时函数为什么可以经过四次不同顺序的扫描而收敛.事实上，第一象限内(i，j)点的走时t ^h _i，j 只依赖于它的左侧点(i-1，j)的走时t ^h _i－1，j 及下方点(i，j-1)的走时t ^h _i，j－1 ，而二者已经计算得到了正确值，这是因为在第一次扫描中二者可以递归地追溯到初始点.这样，在第一次扫描过后，我们得到了第一象限内和x、y正半轴上所有节点的正确的走时(x、y正半轴上所有节点的走时分别只依赖于它的左侧点或下方点的走时).同样，在第二次扫描过后，第二象限内和x负半轴上所有节点的正确的走时也已得到，且第一象限内和x、y正半轴上的所有节点的走时因为已经达到最小且满足离散方程组，它们的值在第二次扫描中不会发生改变.第三次扫描过后，第三象限内和y轴负半轴的所有节点也得到正确的走时，且前两次扫描得到的正确值保持不变.第四次扫描过后，我们就得到了满足离散方程组的所有节点的正确的走时. 图 8 Marmousi模型快速推进法和快速扫描法计算结果( 兰海强等，2012b ) (a)Marmousi模型及用快速推进法和快速扫描法分别计算的走时等值线．红线表示快速推进法计算的结果，黑线表示快速扫描法计算的结果，蓝色的框表示被放大的区域，右侧的色标表示模型速度;(b)为把(a)的局部区域(如(a)中蓝色的框所示的区域)进行放大的结果;(c)两种方法(快速推进法和快速扫描法)计算结果的差异图．走时等值线与差异图单位均为秒. Fig. 8 Marmousi model results computed by the fast marching method and fast sweeping method( Lan et al ., 2012b ) (a)Marmousi model and traveltimes(contours)computed by the fast marching method and fast sweeping method，respectively. Red and black lines are the traveltimes computed by the fast marching method and fast sweeping method，respectively; colourbar shows the velocity;(b)is the enlarged result of a local area(indicating by the blue box in Fig.(a))of Fig.(a);(c)Differences of the traveltimes computed by the two methods，respectively. The units in the traveltime contours and error map both are second(s). 该方法计算步骤主要包含初始化、交替扫描和收敛判断三个步骤.具体如下：

(1)初始化：在震源点上赋初值为0，震源点周围的四个点上赋值为精确值以避免一阶数值误差，这些点上的走时在迭代过程中保持不变.其它网格点上的初值设为一个较大的值(远大于所有计算点最终能算出的走时值)，这些点上的走时在后续计算时会被更新.

(2)交替扫描：在交替的4个方向上，使用 Gauss-Seidel 型迭代法求解离散的非线性方程组.在每一个走时未被固定的网格点，根据(9)式计算出t，与原值比较，取较小值作为当前值t ^new _i，j = min (t ^old _i，j ，t).我们分别按四个方向彻底扫描一遍(1)i=1:I，j=1:J，(2)i=I:1，j=1:J，(3)i=I:1，j=J:1，(4)i=1:I，j=J:1.

(3)收敛判断：如果 t ^new －t ^old ≤δ，那么迭代收敛，结束；如果不收敛，那么继续迭代直至收敛为止.δ是给定的控制迭代停止的极小值.我们可以得到方程(1)的解为 接下来，我们展示一个对比快速推进法和快速扫描法的算例. 图 8a 为Marmousi模型的试算结果.由图可以看出，两种方法计算的结果依然吻合的很好，几乎看不出差别(选取 图 8a 的局部区域放大，见 图 8b )， 图 8c 为两种方法计算的走时的差异图，二者差异很小，这说明对速度纵横向变化剧烈的Marmousi模型，快速推进法和快速扫描法仍具有很好的稳定性和精度.在相同精度的情况下，快速扫描法所耗费的CPU时间约为快速推进法的1/3.

3 基于惠更斯原理的波前构建法

Vinje等(1993) 首次提出了波前构建法(Wavefront Construction，WFC).如 图 9 所示，与其他射线追踪方法相比，该方法将时间参数作为算法在实现过程中的自变量.算法是通过对运动学和动力学射线追踪方程组进行数值分析来模拟波前的传播过程的，并通过各种插值方法来求取地震波的射线路径、走时及振幅等参数.

图 9 波前构建法示意图 图中连接在一起的实线梯形网格为射线网格，虚线矩形网格为规则网格，其走时信息是利用射线网格上的已知走时信息通过插值来求取的. Fig. 9 Schematic diagram of wavefront construction method The trapezoidal grids connected together by solid lines are ray grids. Dashed rectangular grids are regular grids. The travel-time on regular grids is interpolated using known travel-time on the ray grids.

波前构建法主要由两套方程组组成：一是运动学射线追踪方程组，该方程组主要用来计算射线路径及走时；另一个是动力学射线追踪方程组，该方程组主要用来计算射线振幅信息及相位.在计算的过程中，两套方程组可以单独使用，也可以同时使用，依次获得地震波的走时、射线路径及振幅信息.

图 10 带有强速度差介质模型的波前构建法计算结果( Vinje，1993 ) 上图：下层为高速介质的均匀介质中有一圆柱体的模型.圆柱体速度为2 km/s，上层为4 km/s，下层速度为8 km/s.分界面用0.25 km的算子进行光滑.下图：波前构建法计算结果.波前间隔为0.025 s.检波器排列位于模型(0.1，1.1)到(4.8，1.1)km之间. Fig. 10 WFC result for a model with strong variation in velocity( Vinje，1993 ) Upper diagram: A cylindrical body is buried in a homogeneous background with a high-velocity lower layer. The velocities are 2 km/s in the cylinder，4 km/s in the upper layer and 8 km/s in the lower layer. The interfaces are smoothed by a 0.25 km operator. Bottom diagram: The result is calculated by wavefront construction. The wavefront interval is 0.025 s. A receiver line is placed in the model from(0.1，1.1) to (4.8，1.1)km. 该算法的基本思想(以二维情况为例，如 图 9 所示)为：若定义时间τ _i 上的由所有射线组成的波前面为i，那么下一个波前面i+1就可以由波向前传播一个固定的时间长度Δ τ 后的所有射线组成.波在传播过程中的射线参数是通过二元四阶的 Runge—Kutta 法对运动学射线追踪系统常微分方程进行数值分析来确定的.波前构建法的实现过程如下：(1)速度模型的读取及光滑处理.在应用波前构建法射线追踪之前，首先必须将速度模型读入到程序中，在实际地震资料处理中，速度信息一般是以规则网格数据形式给出的，为了能够在波前构建法中应用，必须在读入后对模型进行光滑处理，而在这个过程中，必须考虑不同的光滑算子、次数、因子等因素对模型带来的影响.(2)震源的初始化.在经过光滑处理后，必须通过一定的初始化程序来产生初始的波前，具体的过程是根据震源的初始位置、射线的初始密度、射线的初始角度来计算出初始波前的位置和初始射线的几何属性.(3)数值法计算波前传播.由于波前构建法以时间为参数，随着时间的推移，波前不断往前传播，在此过程中必须选择好插值方法.在上述步骤完成之后，得到的是以规则网格数据形式给出的射线走时和慢度矢量. 图 10 显示了该方法对带有强速度差介质模型的计算结果.震源位于(4.5，1.0)km处(靠近圆柱体).射线路径和波前(0.025 s间隔)显示在 图 10 的下图中.这里存在来自下层高速介质的折射波以及圆柱体后的屏蔽区.

4 基于频率域波动方程法

均匀各向同性介质中频率域声波方程可以表示为

其中，U表示波场，ω表示角频率，v表示速度场； Δ ² 为拉普拉斯算符.其有限差分求解公式可表示为( Shin et al.， 2001)

其中，M，C和K分别是n×n的(n是模型网格节点数)质量矩阵，衰减矩阵和刚度矩阵( Shin et al.， 2001)，f和U分别是n×1的震源和波场数据矢量，m是由模型参数(速度和密度)组成的矢量.方程(12)可以进一步简化为

式中，S=K+iωC－ω ² M.在实际中，把矩阵S分解成一个上三角矩阵和一个下三角矩阵的乘积，然后通过追赶法求解频率域波场U.

这里，为了压制频率域的折叠效应，使用复数频率，即ω ^* =ω+iα.其中，ω ^* 是复数角频率，α是衰减因子.根据傅里叶变换理论可知，复数频率的虚部相当于衰减因子，即对时间域的信号乘以因子e ^－αt .因此衰减因子足够大的时候，时间域的波场可以近似地表示为一个脉冲函数

其中，A ~(x，z，t)=A(x，z，t)e ^－αt ，A(x，z，t)是地下介质某一深度点上的振幅，τ(x，z)是震源到地下某一深度点的地震波走时，δ是脉冲函数.方程(14)在频率域可表示为

公式(15)又可分解为

其中，X=A ~(x，z)cos(－ωτ(x，z))，Y=A ~(x，z)sin(－ωτ(x，z)).由(7)式可以得到初至走时的计算公式，即

同时，初至走时的振幅可表示为

这样，先计算得到波场，然后利用公式(17)和(18)即可计算初至走时及其振幅.

正确的走时和振幅的计算需要选择合适的频率和衰减因子.大的衰减因子可以更好地压制初至走时以后的能量，提高初至走时拾取的精度，但是一个大的衰减因子需要更小的网格间距来减小频散效应，致使计算量增大；小的衰减因子可以提高计算效率，但是初至走时拾取的误差增大.同时，由于频率域波场相位的较大的折叠效应，较大的频率会导致误差的急剧增大.频率和衰减因子的选取详见 秦义龙等(2005) . 图 11 显示了该方法应用于Marmousi模型的计算结果.

图 11 频率域波动方程法Marmousi模型走时计算结果( 秦义龙等，2005 )震源位于距地 表 6 .8 km处，图中数字单位为s Fig. 11 Marmousi model result computed by frequency-domain wave equation method( Qin et al ., 2005 ) The source is located at 6.8 km of the surface， unit in figure is second. 5 结 论

本文主要介绍了目前计算初至波走时的四类主要方法：基于高频近似射线理论的最短路径方法、基于程函方程的数值解方法、基于惠更斯原理的波前构建法、基于频率域波动方程方法.对于上述四类方法，最短路径法将图形理论和地震波射线理论有机结合，并采用网格剖分模型，算法有相对较高的计算精度和较好的稳定性，并且在复杂介质中有较好的适应性，但其计算量较大且计算效率低，另外，SPM是通过减小网格或者扩大网格点(Forward Star技术)来增加角度覆盖次数以获取最小走时，MSPM法是通过插入多个次级节点来增加角度覆盖次数以获取最小走时，与SPM相比有更快的计算速度和较高的计算精度.有限差分法以小的矩形网格单元划分表征模型空间，近似认为地震波在网格单元内以平面波形式传播，在此基础上通过程函方程的有限差分算子计算出地震波从震源出发到所有网格点的初至走时，算法主要包括差分格式和波前扩展方式两个方面，在所有差分格式中，迎风差分格式计算简洁且无条件稳定，在波前扩展方式方面，快速推进法和快速扫描法对不均匀介质依然有很好的稳定性和适用性，且二者的计算精度相当，而快速扫描法的效率较高，因此，程函方程数值解法已具有稳定性好，计算效率高，易于实现等优势，且其还可以通过引入高阶差分格式来提高算法的计算精度.波前构建法可以获取介质模型的走时信息、射线路径信息以及振幅信息，并且能灵活处理地震波在复杂介质中传播出现的多值走时问题，其算法稳定性较好，计算精度较高，但其计算效率较低.频率域波动方程法通过在模拟算法中加入一个衰减因子来压制初至走时之后的能量，从而利用频率域波场的相位项和振幅项计算走时和振幅，其能适应于任意复杂介质，但是计算效率较低. 致 谢 感谢中国科学院地质与地球物理研究所滕吉文院士的悉心指导；感谢白志明副研究员、兰海强博士、刘有山博士生有益的讨论和建议.

Alkhalifah T, Fomel S. 2001. Implementing the fast marching eikonal solver: spherical versus Cartesian coordinates [J]. Geophysical Prospecting, 49(2): 165-178 . Alton K, Mitchell I M. 2008. Fast marching methods for stationary Hamilton-Jacobi equations with axis-aligned anisotropy [J]. SIAM Journal on Numerical Analysis, 47(1): 363-385. Bai C, Greenhalgh S, Zhou B. 2007. 3D ray tracing using a modified shortest-path method [J]. Geophysics, 72(4): T27-T36. Bulant P, Klime L. 1999. Interpolation of ray theory traveltimes within ray cells [J]. Geophysical Journal International, 139(2): 273-282 . Cao S, Greenhalgh S. 1993. Calculation of seismic first-break time field and its ray path distribution using a minimum traveltime tree algorithm [J]. Geophysical Journal International, 114(3): 593-600 . Cerveny V. 2001. Seismic ray theory [M]. Cambridge University Press. Coman R, Gajewski D. 2000. 3D Wavefront Construction Method with Spherical Interpolation [C]. 62nd EAGE Conference & Exhibition. Coman R, Gajewski D. 2001. 3D Traveltime and migration weight computation using wavefront oriented ray tracing [C]. 63rd EAGE Conference & Exhibition. Coman R A. 2003. Computation of multivalued traveltimes in three-dimensional heterogeneous media [J]. Hamburg: University of Hamburg. Dellinger J, Symes W W. 1997. Anisotropic finite-difference traveltimes using a Hamilton-Jacobi solver [C]. Expanded Abstracts, Society of Exploration Geophysicists. 1786-1789 . Ettrich N, Gajewski D. 1996. Wave front construction in smooth media for prestack depth migration [J]. Pure and Applied Geophysics, 148(3/4): 481-502. Fischer R, Lees J M. 1993. Shortest path ray tracing with sparse graphs [J]. Geophysics, 58(7): 987-996. Fomel S. 1997. A variational formulation of the fast marching eikonal solver [J]. SEP-95: Stanford Exploration Project, 127-147 . Gibson Jr R L, Durussel V, Lee K J. 2005. Modeling and velocity analysis with a wavefront-construction algorithm for anisotropic media [J]. Geophysics, 70(4): T63-T74. Han F X, Sun J G, Yang H. 2008. Ray-tracing of wavefront construction by bicubic convolution interpolation [J]. Journal of Jilin University (Earth Science Edition) (in Chinese), 38(2): 336-340. Han F X, Sun J G, Sun Z Q. 2009. Study on grid point positioning and attribute evaluation with the method of wavefront construction [J]. Progress in Geophysics (in Chinese), 24(5): 1748-1756. Hole J. 1992. Nonlinear high-resolution three-dimensional seismic travel time tomography [J]. Journal of Geophysical Research, 97(B5): 6553-6562 . Huang L J, Li Y M, Wu R S. 1992. The wave-front ray tracing method for image reconstruction [J]. Chinese Journal of Geophysics (in Chinese), 35(2): 223-233. Huang Z Y, Zhao J Z. 2004. Three points Fermat ray tracing technique [J]. Progress in Geophysics (in Chinese), 19(1): 201-204. Jeong W K, Whitaker R T. 2008. A fast iterative method for eikonal equations [J]. SIAM Journal on Scientific Computing, 30(5): 2512-2534 . Klime L, Kvasni č ka M. 1994. 3-D network ray tracing [J]. Geophysical Journal International, 116(3): 726-738. Kim S, Cook R. 1999. 3-D traveltime computation using second-order ENO scheme [J]. Geophysics, 64(6): 1867-1876 . Lambare G, Lucio P S, Hanyga A. 1996. Two-dimensional multivalued traveltime and amplitude maps by uniform sampling of a ray field [J]. Geophysical Journal International, 125(2): 584-598 . Lan H, Zhang Z. 2011a. Three-Dimensional Wave-Field Simulation in Heterogeneous Transversely Isotropic Medium with Irregular Free Surface [J]. Bulletin of Seismological Society of America, 101(3): 1354-1370 . Lan H, Zhang Z. 2011b. Comparative study of free-surface boundary condition in two-dimensional finite-difference elastic wave-field simulation [J]. Journal of Geophysics and Engineering, 8: 275-286 . Lan H Q, Zhang Z, Xu T, et al, 2012a. Effects due to the anisotropic stretching of the surface-fitting grid on the traveltime computation for irregular surface by the coordinate transforming method [J]. Chinese Journal of Geophysics (in Chinese), 55(10): 3355-3369. Lan H Q, Zhang Z, Xu T, et al. 2012b. A comparative study on the fast marching and fast sweeping methods in the calculation of first-arrival traveltime field [J]. Progress in Geophysics (in Chinese), 27(5): 1863-1870. Li F, Xu T, Wu Z B, et al. 2013. Segmentally iterative ray tracing in 3D heterogeneous geological models [J]. Chinese Journal of Geophysics (in Chinese), 56(10): 3514-3522, dio:10.6038/cjg20131026 Liu H, Meng F L, Li Y M. 1995. The interface grid method for seeking global minimum travel-time and the corespondent raypath [J]. Chinese Journal of Geophysics (in Chinese), 38(6): 823-832. Liu Y Z, Dong L G, Li P M, et al. 2009. Fresnel volume tomography based on the first arrival of the seismic wave [J]. Chinese Journal of Geophysics (in Chinese), 52(9): 2310-2320 . Liu Y F, Lan H Q. 2012. Study on the numerical solutions of the eikonal equation in curvilinear coordinate system [J]. Chinese Journal of Geophysics (in Chinese), 55(6): 2014-2026 . Lu H Y, Liu Y K, Chang X. 2013. MSFM-based travel-times calculation in complex near-surface model [J]. Chinese Journal of Geophysics (in Chinese), 56(9): 3100-3108. Ma T, Zhang Z. 2014. Calculating ray paths for first-arrival travel times using a topography-dependent eikonal equation solver [J]. Bulletin of Seismological Society of America, 2014, 104(3), doi: 10.1785/0120130172. Moser T J. 1991. Shortest path calculation of seismic rays [J]. Geophysics, 56(1): 59-67. Nakanishi L, Yamaguchi K. 1986. A numerical experiment on nonlinear image reconstruction from first-arrival time for two-dimension island arc structure [J]. Journal of Physics of the Earth, 34(2): 195-201. Nichols D E. 1996. Maximum energy traveltimes calculated in the seismic frequency band [J]. Geophysics, 61(1): 253-263 . Osher S, Shu C W. 1991. High-order essentially nonoscillatory schemes for Hamilton-Jacobi equations [J]. SIAM Journal on Numerical Analysis, 28(4): 907-922 . Pan J S, Zhang X K, Zhao C B, et al. 2006. Computation method of seismic wave travel times and rays in 2D complicated velocity models and its application [J]. Journal of Seismological Research (in Chinese), 29(3): 245-250. Podvin P, Lecomte I. 1991. Finite difference computation of traveltimes in very contrasted velocity models: a massively parallel approach and its associated tools [J]. Geophysical Journal International, 105(1): 271-284 . Qian J, Zhang Y T, Zhao H K. 2007a. A fast sweeping method for static convex Hamilton-Jacobi equations [J]. Journal of Scientific Computing, 31(1): 237-271 . Qian J, Zhang Y T, Zhao H K. 2007b. Fast sweeping methods for Eikonal equations on triangular meshes [J]. SIAM Journal on Numerical Analysis, 45(1): 83-107 . Qin F, Luo Y, Olsen K B, et al. 1992. Finite-difference solution of the eikonal equation along expanding wavefronts [J]. Geophysics, 57(3): 478-487 . Qin Y L, Zhang Z J, Shin Changsoo, et al. 2005. First-arrival traveltime and amplitude calculation from monochromatic two-way wave equation in frequency domain [J]. Chinese Journal of Geophysics (in Chinese), 48(2): 423-428. Rawlinson N, Sambridge M. 2004a. Wave front evolution in strongly heterogeneous layered media using the fast marching method [J]. Geophysical Journal International, 156(3): 631-647 . Rawlinson N, Sambridge M. 2004b. Multiple reflection and transmission phases in complex layered media using a multistage fast marching method [J]. Geophysics, 69(5): 1338-1350. Sava P, Fomel S. 2001. 3-D traveltime computation using Huygens wavefront tracing [J]. Geophysics, 66(3): 883-889 . Schneider Jr W A, Ranzinger K A, Balch A H, et al. 1992. A dynamic programming approach to first arrival traveltime computation in media with arbitrarily distributed velocities [J]. Geophysics, 57(1): 39-50 . Schneider Jr W A. 1995. Robust and efficient upwind finite-difference traveltime calculations in three dimensions [J]. Geophysics, 60(4): 1108-1117 . Sethian J A. 1999. Fast marching methods [J]. SIAM review, 41(2):199-235. Sethian J A, Popovici A M. 1999. 3-D traveltime computation using the fast marching method [J]. Geophysics, 64(2): 516-523 . Shin C, Yoon K, Marfurt K J, et al. 2001. Efficient calculation of a partial-derivative wavefield using reciprocity for seismic imaging and inversion [J]. Geophysics, 66(6): 1856-1863. Shin C, Ko S, Marfurt K J, et al. 2003a. Wave equation calculation of most energetic traveltimes and amplitudes for Kirchhoff prestack migration [J]. Geophysics, 68(6): 2040-2042. Shin C, Ko S, Kim W, et al. 2003b. Traveltime calculations from frequency-domain downward-continuation algorithms [J]. Geophysics, 68(4): 1380-1388 . Sun Z Q, Sun J G, Han F X, et al. 2007. Traveltime computation using fast marching method from rugged topography [J]. Progress in Exploration Geophysics (in Chinese), 30(5): 392-395 . Tang X P, Bai C Y. 2009. Multiple ray tracing within 3-D layered media with the shortest path method [J]. Chinese Journal of Geophysics (in Chinese), 52(10): 2635-2643. van Trier J, Symes W W. 1991. Upwind finite-difference calculation of traveltimes [J]. Geophysics, 56(6): 812-821 . Van Avendonk H J A, Harding A J, Orcutt J A, et al. 2001. Hybrid shortest path and ray bending method for traveltime and raypath calculations [J]. Geophysics, 66(2): 648-653 . Vidale J. 1988. Finite-difference calculation of travel times [J]. Bulletin of the Seismological Society of America, 78(6): 2062-2076. Vidale J. 1990. Finite-difference calculation of traveltimes in three dimensions [J]. Geophysics, 55(5): 521-526 . Vinje V, Iversen E, Gjoystdal H. 1993. Traveltime and amplitude estimation using wavefront construction [J]. Geophysics, 58(8): 1157-1166 . Vinje V, Iversen E, stebl K, et al. 1996a. Estimation of multivalued arrivals in 3D models using wavefront construction-Part I [J]. Geophysical Prospecting, 44(5): 819-842 . Vinje V, Iversen E, stebl K, et al. 1996b. Estimation of multivalued arrivals in 3D models using wavefront construction-Part II [J]. Geophysical Prospecting, 44(5): 843-858 . Vinje V, stebl K, Iversen E, et al. 1999. 3-D ray modeling by wavefront construction in open models [J]. Geophysics, 64(6): 1912-1919. Wang H, Chang X. 2000. 3-D ray tracing method based on graphic structure [J]. Chinese Journal of Geophysics (in Chinese), 43(4): 534-541. Xu K, Wu L, Wang M Y. 1998. Improved Moser ray tracing [J]. Progress in Geophysics (in Chinese), 13(4): 60-66. Xu T, Xu G M, Gao E G, et al. 2004. Block modeling and shooting ray tracing in complex 3D media [J]. Chinese Journal of Geophysics (in Chinese), 47(6): 1118-1126 . Xu T, Ning J R, Liu C C, et al. 2007. Influence of the self-organization of the Earth interior upon the traveltime and amplitude of seismic wave [J]. Chinese Journal of Geophysics (in Chinese), 50(4): 1174-1181. Xu T, Xu G, Gao E, et al. 2006. Block modeling and segmentally iterative ray tracing in complex 3D media [J]. Geophysics, 71(3): T41-T51. Xu T, Zhang Z, Gao E, et al. 2010. Segmentally Iterative Ray Tracing in Complex 2D and 3D Heterogeneous Block Models [J]. Bulletin of the Seismological Society of America, 100(2): 841-850 . Xu T, Li F, Wu Z B, et al, 2014. A successive three-point perturbation method for fast ray tracing in complex 2D and 3D geological models [J]. Tectonophysics, doi: 10.1016/j.tecto.2014.02.012. Yang H, Sun J G, Han F X, et al. 2010. Fast algorithm of the expanding wavefronts finite-difference traveltime calculation based on the three branch tree structure heap sorts [J]. Journal of Jilin University (Earth Science Edition) (in Chinese), 40(1): 188-194. Zhang J Z, Chen S J, Yu D X. 2003. Improvement of shortest path ray tracing method [J]. Progress in Geophysics (in Chinese), 18(1): 146-150. Zhang J Z, Chen S J, Xu C W. 2004. A method of shortest path raytracing with dynamic networks [J]. Chinese Journal of Geophysics (in Chinese), 47(5): 899-904. Zhang M G, Cheng B J, Li X F, et al. 2006a. A fast algorithm of shortest path ray tracing [J]. Chinese Journal of Geophysics (in Chinese), 49(5): 1467-1474. Zhang M G, Jia Y G, Wang M Y, et al. 2006b. A global minimum traveltime raytracing algorithm of wavefront expanding with interface points as secondary sources [J]. Chinese Journal of Geophysics (in Chinese), 49(4): 1169-1175. Zhang Z, Wang G, Teng J, et al. 2000. CDP mapping to obtain the fine structure of the crust and upper mantle from seismic sounding data: an example for the southeastern China [J]. Physics of the Earth and Planetary Interiors, 122(1-2): 133-146 . Zhang Z J, Qin Y L, Chen Y, et al. 2004. Reconstruction of the semblance section for the crust and mantle reflection structure by wide-angle reflection seismic data [J]. Chinese Journal of Geophysics (in Chinese), 47(3): 469-474. Zhao A H, Zhang Z J. 2004. Fast calculation of converted wave travel time in 3-D complex media [J]. Chinese Journal of Geophysics (in Chinese), 47(4): 702-707. Zhao A, Zhang Z, Teng J. 2004. Minimum travel time tree algorithm for seismic ray tracing: improvement in efficiency [J]. Journal of Geophysics and Engineering, 1(4): 245-251 . Zhao A H, Ding Z F. 2005. A double-grid algrithm for calculating traveltimes of wide-angle reflection waves [J]. Chinese Journal of Geophysics (in Chinese), 48(5): 1141-1147. Zhao H. 2005. A fast sweeping method for eikonal equations [J]. Mathematics of computation, 74(250): 603-628. Zhao H. 2007. Parallel implementations of the fast sweeping method [J]. Mathematica Numerica Sinica, 25(4): 421-429 . Zhao L F, Zhu J S, Cao J X, et al. 2003. Ray tracing using ordinal wavefront reconstruction method [J]. Chinese Journal of Geophysics (in Chinese), 46(3): 415-420 . Zhao R, Bai C Y. 2010. Fast multiple ray tracing within complex layered media: The shortest path method based on irregular grid cells [J]. Acta Seismological Sinica (in Chinese), 32(4): 433-444. Zhou H B, Zhang G Q. 1994. Accurate calculation of first-arrival traveltime in the complex areas [J]. Chinese Journal of Geophysics (in Chinese), 37(4): 515-520. Zhu J M, Wang L Y. 1992. Finite difference calculation of seismic travel times [J]. Acta Geophysical Sinica (in Chinese), 35(1): 86-92 . 韩复兴, 孙建国, 杨昊. 2008. 基于二维三次卷积插值算法的波前构建射线追踪[J]. 吉林大学学报(地球科学版), 38(2): 336-340 . 韩复兴, 孙建国, 孙章庆. 2009. 波前构建法中网格点相对定位及属性计算研究[J]. 地球物理学进展, 24(5): 1748-1756 . 黄联捷, 李幼铭, 吴如山. 1992. 用于图像重建的波前法射线追踪[J]. 地球物理学报, 35(2): 223-233 . 黄中玉, 赵金州. 2004. 矩形网格三点 Fermat 射线追踪技术[J]. 地球物理学进展, 19(1): 201-204 . 兰海强, 张智, 徐涛, 等. 2012a. 贴体网格各向异性对坐标变换法求解起伏地表下地震初至波走时的影响[J]. 地球物理学报, 55(10): 3355-3369 . 兰海强, 张智, 徐涛, 等. 2012b. 地震波走时场模拟的快速推进法和快速扫描法比较研究[J]. 地球物理学进展, 27(5): 1863-1870 . 李飞, 徐涛, 武振波, 等. 2013. 三维非均匀地质模型中的逐段迭代射线追踪[J]. 地球物理学报, 56(10): 3514-3522 . 刘洪, 孟凡林, 李幼铭. 1995. 计算最小走时和射线路径的界面网全局方法[J]. 地球物理学报, 38(6): 823-832 . 刘玉柱, 董良国, 李培明, 等. 2009. 初至波菲涅尔体地震层析成像[J]. 地球物理学报, 52(9): 2310-2320 . 刘一峰, 兰海强. 2012. 曲线坐标系程函方程的求解方法研究[J]. 地球物理学报, 55(6): 2014- 2026 . 卢回忆, 刘伊克, 常旭. 2013. 基于 MSFM 的复杂近地表模型走时计算[J]. 地球物理学报, 56(9): 3100-3108 . 潘纪顺, 张先康, 赵成斌, 等. 2006. 二维复杂介质中地震波走时和射线的计算方法及其应用[J]. 地震研究, 29(3): 245-250 . 秦义龙, 张中杰, Shin Changsoo, 等. 2005. 利用单频双程波动方程计算初至走时及其振幅[J]. 地球物理学报, 48(2): 423-428 . [100] 孙章庆, 孙建国, 韩复兴, 等. 2007. 波前快速推进法起伏地表地震波走时计算[J]. 勘探地球物理进展, 30(5): 392-395 . [101] 唐小平, 白超英. 2009. 最短路径算法下三维层状介质中多次波追踪[J]. 地球物理学报, 52(10): 2635-2643 . [102] 王辉, 常旭. 2000. 基于图形结构的三维射线追踪方法[J]. 地球物理学报, 43(4): 534-541 . [103] 许琨, 吴律, 王妙月. 1998. 改进Moser法射线追踪[J]. 地球物理学进展, 13(4): 60-66 . [104] 徐涛, 徐果明, 高尔根, 等. 2004. 三维复杂介质的块状建模和试射射线追踪[J]. 地球物理学报, 47(6): 1118-1126 . [105] 徐涛, 宁俊瑞, 刘春成, 等. 2007. 地球介质自组织性对地震波走时和振幅的影响[J]. 地球物理学报, 50(4): 1174-1181 . [106] 杨昊, 孙建国, 韩复兴, 等. 2010. 基于完全三叉树堆排序的波前扩展有限差分地震波走时快速算法[J]. 吉林大学学报 (地球科学版), 40(1): 188-194 . [107] 张建中, 陈世军, 余大祥. 2003. 最短路径射线追踪方法及其改进[J]. 地球物理学进展, 18(1): 146-150 . [108] 张建中, 陈世军, 徐初伟. 2004. 动态网络最短路径射线追踪[J]. 地球物理学报, 47(5): 899-904 . [109] 张美根, 程冰洁, 李小凡, 等. 2006a. 一种最短路径射线追踪的快速算法[J]. 地球物理学报, 49(5): 1467-1474 . [110] 张美根, 贾豫葛, 王妙月, 等. 2006b. 界面二次源波前扩展法全局最小走时射线追踪技术[J]. 地球物理学报, 49(4): 1169-1175 . [111] 张中杰, 秦义龙, 陈赟, 等. 2004. 由宽角反射地震资料重建壳幔反射结构的相似性剖面[J]. 地球物理学报, 47(3): 469-474 . [112] 赵爱华, 张中杰. 2004. 三维复杂介质中转换波走时快速计算[J]. 地球物理学报, 47(4): 702-707 . [113] 赵爱华, 丁志峰. 2005. 宽角反射地震波走时模拟的双重网格法[J]. 地球物理学报, 48(5): 1141-1147 . [114] 赵连锋, 朱介寿, 曹俊兴, 等. 2003. 有序波前重建法的射线追踪[J]. 地球物理学报, 46(3): 415-420 . [115] 赵瑞, 白超英. 2010. 复杂层状模型中多次波快速追踪--一种基于非规则网格的最短路径算法[J]. 地震学报, 32(4): 433-444 . [116] 周洪波, 张关泉. 1994. 复杂构造区域的初至波走时计算[J]. 地球物理学报, 37(4): 515-520 . [117]