2017, Vol. 60
2. 长安大学计算地球物理研究所, 西安 710054
2. Institute of Computing Geophysics, Chang'an University, Xi'an 710054, China
地震射线追踪技术是指在高频近似的条件下，准确地计算和描述地震波主能量传播轨迹(即射线)的方法，该技术广泛应用于建立全球走时表(Jeffreys and Bullen, 1940; Kennett and Engdahl, 1991; Kennett, 2005)、地震定位(Husen and Kissling, 2001; Bai and Greenhalgh, 2006; 赵爱华等, 2015)以及走时层析成像(Aki and Lee, 1976; 华标龙和刘福田, 1995; 黄国娇和白超英, 2013a)等工作，是一种重要的正演方法技术.
地震射线追踪方法有很多(张钋等, 2000; Cerveny, 2001; 李强和白超英, 2012)，传统的方法有基于初值问题的打靶法和基于边值问题的弯曲法以及相应的改进方法.打靶法首先给定射线出射角度，根据由震源发出的一条射线到达接收点的情况对射线出射角进行调整，直至射线正好落在检波点上(马争鸣和李衍达, 1991; Cerveny, 2001)；弯曲法则是从震源与接收点之间的一条假想初始路径开始，根据最小走时原理对初始路径进行扰动，从而得到该炮检排列的最小走时及路径(Julian and Gubbins, 1977; 田玥和陈晓非, 2005; 李飞等, 2013).
上世纪80年代后期基于全局解的波前扩展算法开始流行，较为实用的主要有三种：其一是有限差分求解程函方程法(简称FD，Vidale, 1988, 1990)，Qin等(1992)在FD法的基础上实现了扩展波前的递推算法.该算法随后被扩展到快速行进法(简称FMM，Sethian，1999)，Rawlinson和Sambridge(2004)在FMM算法中引入分区多步计算技术，实现了多震相地震射线的追踪计算.其二是线性旅行时插值法(简称LTI，Asakawa and Kawanaka, 1993)，张东等(2009)提出基于LTI和网格界面剖分的三维地震射线追踪算法；马德堂等(2014)将LTI算法应用于三维各向异性介质中初至波射线追踪.其三是最短路径算法(Nakanishi and Yamaguchi, 1986; Moser, 1991；张建中等, 2004; 张美根等, 2006; Bai et al., 2007).王辉和常旭(2000)引入Bresenham画线算法计算节点走时，结合使用快速排序法和插入排序算法选取子波源点，提高了走时和射线路径的计算精度以及计算效率；Zhao等(2004)采用动态调整波前传播区域的方法，可有效减少计算量，同时保证较高的计算精度；Bai等(2009)将改进型最短路径法与分区多步计算技术相结合，实现了2D/3D复杂层状地学模型中多震相地震射线的追踪计算.
上述算法大多建立在直角坐标系下，当计算区域范围较大时(震中距大于200 km)，地球曲率的影响不可忽略(Snoke and Lahr, 2001).Bijwaard和Spakman(1999)在三维球坐标系下分别采用最短路径法和弯曲法计算全球震相；粱春涛等(2001)采用试射法追踪球坐标系下射线路径；de Kool等(2006)将FMM法推广到球坐标系下，可快速准确地计算多种震相；Huang等(2013b)采用改进型最短路径算法，结合分区多步计算技术，实现了全球49种主要震相的精确计算，并将之应用于多参数同时反演研究(黄国娇等, 2015).事实上地球并不是规则的球体，而是两极略扁的椭球体，Dziewonski和Gilbert(1976)研究表明，椭球体和球体上射线路径和走时存在一定差异，主要受径向速度梯度影响.进行高精度地震波走时计算的过程中，有必要进行相应的校正，通常是先在球坐标系下计算走时，然后对走时进行椭球校正(Dziewonski and Gilbert, 1976).该方法所用走时校正公式较为复杂，校正系数求取困难，适用性受到限制.
本文将分区多步最短路径算法(唐小平和白超英, 2009; Bai et al., 2009)推广至椭球坐标系中，实现了椭球坐标系下多震相射线追踪计算.为检验算法的计算精度和有效性，采用AK135地球模型计算P波和PcP波走时，计算结果与AK135理论走时表的最大误差小于0.1 s，证明该算法精确有效.
2 椭球坐标系下分区多步最短路径算法 2.1 网格剖分对于一般的椭球体，需要明确其三个相互正交轴的长度，才能准确描述该椭球体的形状.作为特殊情况，三个轴中某两个轴长度相等时，为旋转椭球体；三个轴长度均相等时，则为球体.因此，椭球体网格剖分与球体的剖分有相似的方面，可按照实际需求，选用合适的网格对椭球体进行剖分.在经度或纬度方向上，网格间距不变；但由于椭球面上各点到椭球中心的距离并不相等，因此深度方向的间距并不相同，这一点有别于球体.为了便于剖分，一般不直接给定深度方向的间距，而是按照径向网格数将椭球面上一点到中心的距离平均分配.在图 1所示坐标系统中，椭球面上任意点的位置存在如下转换关系:
|
图 1 坐标系统示意图 Fig. 1 Diagrammatic illustration to show the coordinates system |
|
(1) |
其中x、y、z为直角坐标，R为该点到椭球中心的距离，θ∈[－90°, 90°]为纬度，负值表示南纬，正值表示北纬；φ∈[－180°, 180°]为经度，负值表示西经，正值表示东经.在直角坐标系下，椭球面方程为
|
(2) |
其中a、b、c分别为椭球三个半轴的长度.将公式(1) 代入公式(2)，经过整理即可得到椭球面上任意点到椭球中心的距离的表达式：
|
(3) |
图 2a为剖分结果示意图，所采用椭球体的三个半轴长度分别为a=500 km, b=600 km, c=400 km.图 2b-2d分别展示了过中心网格单元、过极点网格单元、以及其它网格单元的形状.网格单元的角点为主节点(图 2中圆点所示)，速度采样在主节点上进行.
|
图 2 模型参数化示意图 (a)椭球剖分示意图; (b)过中心网格单元; (c)过极点网格单元; (d)其他网格单元.图中圆点为主节点. Fig. 2 Diagrammatic illustration to show model parameterization (a) Cells in ellipsoid; (b) A cell near center; (c)A cell near pole; (d) Other cell. In figure the dots represent primary nodes. |
为了提高计算精度，需在主节点之间插入次级节点(图 3中灰点所示)，次级节点仅分布在单元的表面，单元内部无次级节点，但震源和台站可位于其中(Bai et al., 2007).图 3a-3c分别展示了0°经圈、90°经圈以及0°纬圈剖面上网格单元和主、次节点的分布情况，可以看出，越靠近中心的区域，两主节点间的弧长越短，插入的次级节点数越少；反之，主节点间的弧长越大，插入的次级节点数越多.
|
图 3 模型剖面情况 (a)—(c)分别是0°经圈、90°经圈、以及0°纬圈剖面上网格单元和主、次节点的分布情况; (d) 0°纬圈剖面网格和界面分布情况.图中黑点为主节点, 灰点为次级节点, 虚线代表界面. Fig. 3 Slice of the model parameterization (a)—(c) Show cross section view at longitude=0°, longitude=90°, latitude=0°, respectively; (d) Cells and interfaces at latitude=0°. In figure black dots represent primary nodes, gray dots represent secondary nodes and dashed lines represent interfaces. |
速度界面主要由速度模型中的不连续面确定，并将之离散化.为保证计算精度，界面上离散点的间距应不低于次级节点的间距.如果介质中速度连续或者分布复杂，界面的精确位置很难用速度差异来判断，此时可人为给定界面.图 3d展示了0°纬圈剖面上网格和界面分布情况，其中虚线所示为两个界面，它们将模型空间划分为三个相互独立的区域.为使成图简洁、清晰，图 3d中主、次级节点均未给出.
需要注意的是，模型中心附近网格相互连接，使得中心成为一个奇异点，可能导致较大计算误差.为避免这个问题，实际处理中可将奇异点附近无限小半径范围内的区域移除，形成一个球面，从而使中心附近网格分离开来(Huang et al., 2014).
2.2 速度插值由于速度采样只在主节点进行，次级节点(包括震源和台站)的速度值由其所在单元的主节点速度值通过插值得到.根据网格单元形状的不同，分别采用不同的插值方法：
(1) 若节点位于图 2b或图 2c所示的不规则单元上，则需先将该单元剖分为多个四面体(四面体的四个角点必须是主节点)，然后再由节点所在四面体的四个主节点上定义的四面体体积插值函数给出(Kasson et al., 1993):
|
(4) |
其中，x为待求点坐标，vj表示节点所在四面体单元各主节点的速度值，uj为体积坐标(参见图 4所示).图 4中P点在四面体T1T2T3T4内的体积坐标u1可以定义为：u1=[PT2T3T4]/[T1T2T3T4]，其中[PT2T3T4]、[T1T2T3T4]分别表示四面体PT2T3T4和T1T2T3T4的体积，对于u2, u3, u4以此类推.
|
图 4 四面体单元内速度插值示意图 Fig. 4 Velocity interpolation schematic map in a tetrahedral cell |
(2) 若节点位于图 2d所示的规则单元上，其速度可由单元的八个主节点上的速度值通过三次线性拉格朗日插值得到(Bai et al, 2007)：
|
(5) |
其中xi(i=1, 2, …, 8) 和x分别为该单元内8个主节点和待求次级节点的矢量坐标，v(xi)(i=1, 2, …, 8) 为该单元主节点的速度值.
2.3 节点走时计算方法由于单元网格间距较小，对于次级震源i点附近任意节点j处的走时(i和j位于相同单元内)，可用以下公式(两点间距离除以平均速度)计算：
|
(6) |
其中，ti是地震波到达i节点的最小走时，Cj为待求点j的所有邻点的集合，V(xi)和V(xj)分别是i和j节点上的速度值，D(xi, xj)是节点i与节点j之间的距离.椭球体内两点距离可由余弦定理求出：
|
(7) |
式中R由公式(3) 求取，zi和zj分别是i、j两点的深度，γ为点i、j与球心组成的圆心角，且cosγ由下式计算(熊介, 1988):
|
(8) |
追踪地震射线时，为减少计算量，避免对射线不经过的区域进行计算，可以先对模型进行分区处理：依据界面起伏情况，将速度模型分成相互独立的壳状(或层状)区域，相邻区域由界面连接.然后根据所要追踪地震波的震相类型，从震源所在的区域开始，分步、依次计算射线将要经过的某一个或多个区域，直至波前到达检波器.以图 3d所示模型为例，假设震源和台站均位于地表，现追踪来自第一个界面的反射转换波P1S1(上标代表上行波，下标代表下行波，数字表示分区号)，计算步骤为：1) 选择第一区作为计算区域(下行波经过的区域)，并调用相应的P波速度模型，由震源开始计算，直至波前下行至界面一；2) 仅保留界面一上离散点的走时和路径，同样选择第一区作为计算区域(上行波经过的区域)，调用相应的S波速度模型，并从界面一上选取走时最小节点作为次级震源，计算该区域内走时和路径.按照上述原理可实现多震相地震射线的追踪计算，详细介绍可参考文献(唐小平和白超英, 2009).
3 射线追踪实例 3.1 计算精度和效率分析为分析算法的计算精度，首先构建一个均匀椭球模型进行试算，地震波速度4.0 km·s-1.考虑到计算量和内存开销，所选模型不大，三个轴的半径分别为：a=90 km，b=100 km，c=80 km.模型在经度、纬度和深度方向各剖分40个网格，共64000个网格单元，68921个主节点.以点(0°, 90°, 0 km)为震源(坐标依次为经度、纬度和深度，下同)，分别计算不同次级节点间距下的直达波走时场，并与解析解作对比.计算所用处理器型号为：Intel(R) Core(TM) i5-4570@3.2 GHz，结果见表 1.从表 1可明显看出，随着次级节点间距减小，总节点数明显增多，最大误差和平均误差均显著减小，但计算耗时更长.表 1仅讨论了网格间距一定时，减小次级节点间距时的情形，而当次级节点间距一定时，在一定范围内减小网格间距并不能显著提高计算精度.在均匀模型中，各离散点的速度采样值不会受到网格间距的影响，计算精度主要受射线出射方向控制.单元各边上插入的次级节点越多，射线可能的出射方向也就越多，射线路径将更加符合真实情况，故而计算误差越小.因此，可根据速度模型的复杂情况选用合适的网格单元尺寸(速度变化平缓时可选用稍大的网格单元)和合适的次级节点间距，以兼顾计算的效率和精度.
|
表 1 不同次级节点间距下的误差情况 Table 1 Traveltime errors with different interval of secondary nodes |
图 5直观反映了次级节点间距为0.250 km时计算所得的走时场和误差分布情况.图 5a为0°纬圈剖面上的走时场，图 5b为对应的误差分布，误差呈环状分布，中心附近值最大，由中心向外围逐渐减小；图 5c为0°经圈剖面上的走时场，图 5d为对应的误差分布，误差呈扫帚状分布，中心点下方区域误差较大.误差分布与节点和网格的分布规律相关，出现上述现象是由于靠近中心附近次级节点较少，使得射线出射方向受限，导致计算误差累积明显.
|
图 5 走时场及误差分布剖面图 (a)—(b) 0°纬圈剖面上走时场和误差分布; (c)—(d) 0°经圈剖面上走时场和误差分布. Fig. 5 Traveltime contour and error distribution at cross section view (a)—(b) Traveltime contour and error at latitude=0°; (c)—(d) Traveltime contour and error at longitude=0°. |
为了更加符合真实情况，我们构建一均匀区域模型(三个半轴长度a=b=c=6371 km)进行计算，模型速度4.0 km·s-1，模型大小20°×20°×1000 km，网格大小1°×1°×100 km，次级节点间距2.5 km，震源置于模型地表正中位置(图 6中星号所示).计算结果最大走时误差0.052 s，平均走时误差17.52 ms.图 6为Latitude=0°处的结果剖面，其中图 6a为走时场，图 6b为误差分布.可以发现，与震源呈45°方向上的误差稍大，这与节点分布规律相关.本例中虽然采用了较大的网格，但由于主节点间插入的次级节点较多(39个次级节点)，计算误差依然较小，这与前面的讨论结果一致.
|
图 6 纬度0°剖面结果 (a)走时场; (b)误差分布. Fig. 6 Traveltime contour map (a) and error (b) at cross section of latitude=0° |
选取一个区域模型进行试算：经度范围(-10°, 10°)、纬度范围(-50°, 50°)、深度范围(0, 2900 km)，即模型大小20°×100°×2900 km，网格大小1°×1°×25 km，次级节点间距2.78 km.为便于与AK135走时表对比，采用球对称模型(a=b=c=6371 km)及AK135速度模型.震源位于(0°, -50°, 0 km)，20个台站沿0°经线由南往北每间隔5°均匀布置.图 7给出了P波(虚线)和PcP波(实线)射线路径(为使成图清晰，对P波路径做了镜像处理).图 8为计算走时与AK135理论走时表的差异(Tcal-Tak135)，可以看出，不管是P波还是PcP波，最大绝对走时差均小于0.1 s，表明该算法精确有效.
|
图 7 P波(虚线)和PcP波(实线)射线路径 Fig. 7 Ray path of P (dashed lines) and PcP (solid lines) wave |
|
图 8 P波和PcP波走时与AK135走时表的差异 Fig. 8 Traveltime error between calculated results and AK135 traveltime tables for P and PcP wave |
将上述模型由球体改为旋转椭球体模型(a=b=6378.140 km, c=6356.755 km)，依然采用AK135速度模型(速度随深度变化)其余参数不变，再次计算P波和PcP波走时.两者结果中射线路径差异较小，这里不再给出，但走时差(T椭球-T球)显著(见图 9).详细分析图 9，可以得出这样的结论：对于P波而言，震中距小于20°时，两者走时差很小，震中距在20°~70°范围内，(T椭球-T球)之差大约稳定在0.15 s的水平，而当震中距大于70°时，则上述走时差由正值变为负值，且随着震中距的进一步增大，这一负值也进一步增大.而对于PcP波而言，震中距小于30°时，(T椭球-T球)之差大约稳定在0.15 s的水平，而当震中距大于30°时，这一走时差由正值向负值线性变化.震中距95°时，两者走时差达到最大值(接近0.6 s).上述结论说明，对于不同的震相而言，应分不同的震中距进行椭球校正.
|
图 9 椭球坐标系下计算所得P、PcP走时与AK135走时表的误差 Fig. 9 Error between the traveltimes calculated in ellipsoidal coordinates and AK135 traveltime Table for P and PcP waves |
球坐标系下计算所得走时与实际情况之间存在一定误差，需要进行椭球走时校正.AK135全球理论走时表中的椭球校正采用Dziewonski和Gilbert(1976)给出的走时校正公式，并按照5°间距给出主要震相的校正系数，其中P波校正范围(0°, 95°)，PcP波校正范围(0°, 90°).据此，我们给出了椭球坐标系下计算的P、PcP走时与Ak135理论走时表中P波和PcP波经过校正后的走时之差(T椭球-TAK135校)，见图 10所示.理论上来说，若椭球校正合理，则(T椭球-TAK135校)的值应在0值附近波动.但从图 10中可以看出，其椭球校正在震中距25°~70°范围内基本合理(两者之差接近0值)，而对于震中距小于25°时，对P波的校正过大，而对PcP波的校正过小；而当震中距大于70°时，无论是对P波还是PcP波的校正都不够.因此，椭球坐标下走时计算可以直接用来对AK135走时表进行相关震相的椭球走时校正，这也是进行椭球坐标系下射线追踪的另一个应用.
|
图 10 椭球坐标系下计算所得P、PcP走时与经椭球校正后AK135走时表的误差 Fig. 10 Error between the traveltimes calculated in ellipsoidal coordinates and corrected AK135 traveltime Table for P and PcP waves |
除了常见的地震震相，理论上来说，该算法可以追踪任意震相的地震射线.下例在参考椭球体(a=b=6378.140 km, c=6356.755 km)上进行，所选模型(见图 11)：经度范围(-10°, 10°)、纬度范围(-5°, 5°)、深度范围(0, 700 km)，即模型大小20°×10°×700 km，网格大小0.5°×0.5°×50 km，次级节点间距5.5 km.模型中410 km和660 km深度处各存在一反射界面，界面呈正弦状起伏，最大扰动幅度分别是10 km和15 km.从模型地表到模型底部，P波背景速度由5.8 km·s-1线性递增到10.8 km·s-1，且
|
图 11 直达P波(虚线)和反射P660P波(实线)射线路径 Fig. 11 Ray path of P (dashed lines) and reflected P660P (solid lines) wave |
本例中对4种震相的地震射线进行了追踪计算.首先追踪直达P波射线路径(图 11中虚线所示)以及由震源出发经界面二反射到达台站的反射P660P波射线路径(图 11中实线所示).然后计算多次反射波射线路径.图 12中灰线所示为由震源出发的P波，依次经地表和界面一反射后到达台站的pP410P波；图 12中黑线所示为由震源出发的P波，经地表反射后到达界面二，然后反射转换为S波并到达台站的pP660S波.
|
图 12 多次波pP410P(灰线)和pP660S(黑线)射线路径 Fig. 12 Ray path of pP410P(gray lines) and pP660S(black lines) wave |
考虑到地球本身非球对称的特性，将分区多步最短路径算法应用于椭球坐标系中，实现了椭球坐标系下多震相射线追踪.为保证计算精度，主节点之间须插入一定数量的次级节点，但由于每个网格的实际大小不同(距离球心越近，网格越小)，需在每个网格单元面上插入不同数量的次级节点，可在保证计算精度的前提下减少节点总数.均匀模型中计算结果表明，在网格单元剖分不变的情况下，该算法计算精度较高，且允许通过减小次级节点间距来提高计算精度，但同时应注意到计算量将会增大；采用AK135地球模型时，计算所得P波和PcP波走时与AK135走时表的误差小于0.1 s，体现了该算法较高的计算精度和实用性；使用球对称模型和椭球模型计算所得的结果差异显著(特别是震中距较大时)，这种情形下有必要采用椭球坐标系.另外，该算法适用于任意形状的椭球体，能够追踪计算多种震相的地震波射线路径和走时，且不需要额外进行走时校正，算法简单，并能够减少计算量.另一方面，椭球坐标下的走时计算可直接用于对AK135全球理论走时表进行椭球校正.
| Aki K, Lee K W. 1976. Determination of three-dimensional velocity anomalies under a seismic array using first P arrival times from local earthquakes:1.A homogeneous initial model. Journal of Geophysical Research, 81: 4381-4399. DOI:10.1029/JB081i023p04381 | |
| Asakawa E, Kawanaka T. 1993. Seismic ray tracing using linear traveltime interpolation. Geophysical Prospecting, 41(1): 99-111. DOI:10.1111/gpr.1993.41.issue-1 | |
| Bai C Y, Greenhalgh S. 2006. 3D Local earthquake hypocenter determination with an irregular shortest-path method. Bulletin of the Seismological Society of America, 96(6): 2257-2268. DOI:10.1785/0120040178 | |
| Bai C Y, Greenhalgh S, Zhou B. 2007. 3D ray tracing using a modified shortest-path method. Geophysics, 72(4): T27-T36. DOI:10.1190/1.2732549 | |
| Bai C Y, Tang X P, Zhao R. 2009. 2-D/3-D multiply transmitted, converted and reflected arrivals in complex layered media with the modified shortest path method. Geophysical Journal International, 179(1): 201-214. DOI:10.1111/gji.2009.179.issue-1 | |
| Bijwaard H, Spakman W. 1999. Fast kinematic ray tracing of first-and later-arriving global seismic phases. Geophysical Journal International, 139(2): 359-369. DOI:10.1046/j.1365-246x.1999.00950.x | |
| Cerveny V. 2001. Seismic Ray Theory. Cambridge, UK: Cambridge University Press. | |
| De Kool M, Rawlinson N, Sambridge M. 2006. A practical grid-based method for tracking multiple refraction and reflection phases in three-dimensional heterogeneous media. Geophysical Journal International, 167(1): 253-270. DOI:10.1111/gji.2006.167.issue-1 | |
| Dziewonski A M, Gilbert F. 1976. The effect of small, aspherical perturbations on travel times and a re-examination of the corrections for ellipticity. Geophysical Journal International, 44(1): 7-17. DOI:10.1111/j.1365-246X.1976.tb00271.x | |
| Hua B L, Liu F T. 1995. Joint reflection imaging of velocity and interface-theory and method. Chinese Journal of Geophysics, 38(6): 750-756. | |
| Huang G J, Bai C Y. 2013a. Simultaneous inversion of three model parameters with multiple classes of arrival times. Chinese Journal of Geophysics, 56(12): 4215-4225. DOI:10.6038/cjg20131224 | |
| Huang G J, Bai C Y, Greenhalgh S. 2013b. Fast and accurate global multiphase arrival tracking:the irregular shortest-path method in a 3-D spherical earth model. Geophysical Journal International, 194(3): 1878-1892. DOI:10.1093/gji/ggt204 | |
| Huang G J, Bai C Y, Li X W, et al. 2014. 3-D simultaneous inversion for velocity and reflector geometry using multiple classes of arrivals in a spherical coordinate frame. Journal of Seismology, 18(1): 123-135. DOI:10.1007/s10950-013-9406-z | |
| Huang G J, Bai C Y, Qian W. 2015. Simultaneous inversion of three model parameters with multiple phases of arrival times in spherical coordinates. Chinese Journal of Geophysics, 58(10): 3627-3638. DOI:10.6038/cjg20151016 | |
| Husen S, Kissling E. 2001. Local earthquake tomography between rays and waves:fat ray tomography. Physics of the Earth and Planetary Interiors, 123(2-4): 127-147. DOI:10.1016/S0031-9201(00)00206-5 | |
| Jeffreys H, Bullen K E. 1940. Seismological tables. London: British Association for the Advancement of Science. | |
| Julian B R, Gubbins D. 1977. Three-dimensional seismic ray tracing. J. Geophys. Gubbins, 43(1): 95-113. | |
| Kasson J M, Plouffe W, Nin S I. 1993. Tetrahedral interpolation technique for color space conversion.//Proceedings of SPIE 1909-Device-Independent Color Imaging and Imaging Systems Integration. San Jose, CA:SPIE. | |
| Kennett B L N, Engdahl E R. 1991. Traveltimes for global earthquake location and phase identification. Geophysical Journal International, 105(2): 429-465. DOI:10.1111/gji.1991.105.issue-2 | |
| Kennett B N. 2005. Seismological Tables:ak135. Canberra, Australia: Research School of Earth Sciences, Australian National University. | |
| Li F, Xu T, Wu Z B, et al. 2013. Segmentally iterative ray tracing in 3-D heterogeneous geological models. Chinese Journal of Geophysics, 56(10): 3514-3522. DOI:10.6038/cjg20131026 | |
| Li Q, Bai C Y. 2012. Review on seismic wavefront and raytracing in complex media. Progress in Geophysics, 27(1): 92-104. DOI:10.6038/j.issn.1004-2903.2012.01.011 | |
| Liang C T, Zhu J S, Cao J M, et al. 2001. The methods and applications of 3D ray tracing on spherical coordinate. Computing Techniques for Geophysical and Geochemical Exploration, 23(3): 199-205, 231. | |
| Ma D T, Yang C Y, Li X X. 2014. Preliminary waves ray tracing by linear travel-time interpolation method in 3-D transversely isotropic medium. Progress in Geophysics, 29(3): 1201-1205. DOI:10.6038/pg20140327 | |
| Ma Z M, Li Y D. 1991. Two-step ray tracing method. Acta Geophysica Sinica, 34(4): 501-508, 533. | |
| Moser J T. 1991. Shortest path calculation of seismic rays. Geophysics, 56(1): 59-67. DOI:10.1190/1.1442958 | |
| Nakanishi I, Yamaguchi K. 1986. A numerical experiment on nonlinear image reconstruction from first-arrival time for two-dimensional island arc structure. Journal of Physics of the Earth, 34(2): 195-201. DOI:10.4294/jpe1952.34.195 | |
| Qin F H, Luo Y, Olsen K B, et al. 1992. Finite-difference solution of the Eikonal equation along expanding wavefronts. Geophysics, 57(3): 478-487. DOI:10.1190/1.1443263 | |
| Rawlinsonand N, Sambridge M. 2004. Multiple reflection and transmission phases in complex layered media using a multistage fast marching method. Geophysics, 69(5): 1338-1350. DOI:10.1190/1.1801950 | |
| Sethian J A. 1999. Fast Marching Methods. SIAM Review, 41(2): 199-235. DOI:10.1137/S0036144598347059 | |
| Snoke J A, Lahr J C. 2001. Locating earthquakes:at what distance can the earth no longer be treated as flat. Seismological Research Letters, 72(5): 538-541. DOI:10.1785/gssrl.72.5.538 | |
| Tang X P, Bai C Y. 2009. Multiple ray tracing within 3-D layered media with the shortest path method. Chinese Journal of Geophysics, 52(10): 2635-2643. DOI:10.3969/j.issn.0001-5733.2009.10.024 | |
| Tian Y, Chen X F. 2005. A rapid and accurate two-point ray tracing method in horizontally layered velocity model. Acta Seismologica Sinica, 27(2): 147-154. | |
| Vidale J. 1988. Finite-difference calculation of travel times. Bulletin of the Seismological Society America, 78(6): 2062-2076. | |
| Vidale J. 1990. Finite-difference calculation of travel time in three dimensions. Geophysics, 55(5): 521-526. DOI:10.3321/j.issn:0001-5733.2000.04.014 | |
| Wang H, Chang X. 2000. 3-D ray tracing method based on graphic structure. Chinese Journal of Geophysics, 43(4): 535-541. DOI:10.3321/j.issn:0001-5733.2000.04.014 | |
| Xiong J. 1988. Ellipsoidal Geodesy. Beijing: Publishing House of Liberation Army. | |
| Zhang D, Fu X R, Yang Y, et al. 2009. 3-D Seismic ray tracing algorithm based on LTI and partition of grid interface. Chinese Journal of Geophysics, 52(9): 2370-2376. DOI:10.3969/j.issn.0001-5733.2009.09.023 | |
| Zhang J Z, Chen S J, Xu C W. 2004. A method of shortest path raytracing with dynamic networks. Chinese Journal of Geophysics, 47(5): 899-904. DOI:10.3321/j.issn:0001-5733.2004.05.023 | |
| Zhang M G, Cheng B J, Li X F, et al. 2006. A fast algorithm of shortest path ray tracing. Chinese Journal of Geophysics, 49(5): 1467-1474. DOI:10.3321/j.issn:0001-5733.2006.05.026 | |
| Zhang P, Liu H, Li Y M. 2000. The situation and progress of ray tracing method research. Progress in Geophysics, 15(1): 36-45. DOI:10.3969/j.issn.1004-2903.2000.01.002 | |
| Zhao A H, Zhang Z J, Teng J W. 2004. Minimum travel time tree algorithm for seismic ray tracing:improvement in efficiency. Journal of Geophysics and Engineering, 1(4): 245-251. DOI:10.1088/1742-2132/1/4/001 | |
| Zhao A H, Ding Z F, Bai Z M. 2015. Improvement of the ray-tracing based method calculating hypocentral loci for earthquake location. Chinese Journal of Geophysics, 58(9): 3272-3285. DOI:10.6038/cjg20150922 | |
| 华标龙, 刘福田. 1995. 界面和速度反射联合成像-理论与方法. 地球物理学报, 38(6): 750–756. | |
| 黄国娇, 白超英. 2013a. 多震相走时联合三参数同时反演成像. 地球物理学报, 56(12): 4215–4225. DOI:10.6038/cjg20131224 | |
| 黄国娇, 白超英, 钱卫. 2015. 球坐标系下多震相走时三参数同时反演成像. 地球物理学报, 58(10): 3627–3638. DOI:10.6038/cjg20151016 | |
| 李飞, 徐涛, 武振波, 等. 2013. 三维非均匀地质模型中的逐段迭代射线追踪. 地球物理学报, 56(10): 3514–3522. DOI:10.6038/cjg20131026 | |
| 李强, 白超英. 2012. 复杂介质中地震波前及射线追踪综述. 地球物理学进展, 27(1): 92–104. DOI:10.6038/j.issn.1004-2903.2012.01.011 | |
| 粱春涛, 朱介寿, 曹家敏, 等. 2001. 球坐标系下三维射线追踪的方法和应用. 物探化探计算技术, 23(3): 199–205, 231. | |
| 马德堂, 杨春雨, 李绪宣. 2014. 基于双线性插值的三维横向各向同性介质初至波射线追踪. 地球物理学进展, 29(3): 1201–1205. DOI:10.6038/pg20140327 | |
| 马争鸣, 李衍达. 1991. 二步法射线追踪. 地球物理学报, 34(4): 501–508, 533. | |
| 唐小平, 白超英. 2009. 最短路径算法下三维层状介质中多次波追踪. 地球物理学报, 52(10): 2635–2643. DOI:10.3969/j.issn.0001-5733.2009.10.024 | |
| 田玥, 陈晓非. 2005. 水平层状介质中的快速两点间射线追踪方法. 地震学报, 27(2): 147–154. | |
| 王辉, 常旭. 2000. 基于图形结构的三维射线追踪方法. 地球物理学报, 43(4): 535–541. | |
| 熊介. 1988. 椭球大地测量学. 北京: 解放军出版社. | |
| 张东, 傅相如, 杨艳, 等. 2009. 基于LTI和网格界面剖分的三维地震射线追踪算法. 地球物理学报, 52(9): 2370–2376. DOI:10.3969/j.issn.0001-5733.2009.09.023 | |
| 张建中, 陈世军, 徐初伟. 2004. 动态网络最短路径射线追踪. 地球物理学报, 47(5): 899–904. | |
| 张美根, 程冰洁, 李小凡, 等. 2006. 一种最短路径射线追踪的快速算法. 地球物理学报, 49(5): 1467–1474. DOI:10.3321/j.issn:0001-5733.2006.05.026 | |
| 张钋, 刘洪, 李幼铭. 2000. 射线追踪方法的发展现状. 地球物理学进展, 15(1): 36–45. DOI:10.3969/j.issn.1004-2903.2000.01.002 | |
| 赵爱华, 丁志峰, 白志明. 2015. 基于射线追踪技术计算地震定位中震源轨迹的改进方法. 地球物理学报, 58(9): 3272–3285. DOI:10.6038/cjg20150922 | |