2019, Vol. 62
2. 中国科学院大学, 北京 100049;
3. 中国石油集团公司东方地球物理公司物探技术研究中心, 河北涿州 072751
2. University of Chinese Academy of Sciences, Beijing 100049, China;
3. Research & Development Center, BGP, CNPC, Zhuozhou Hebei 072751, China
地震射线理论对于描述地震波的传播起着重要的作用.射线追踪技术在勘探地震, 天然地震等领域都有着广泛的应用, 比如地震勘探资料的偏移成像, 地球物理层析成像、震源定位等.射线追踪旨在解决给定激发点与接收点的两点射线追踪问题, 对于射线追踪方法的研究已有很长的历史.传统射线追踪方法主要包括以斯奈尔定理为基础的试射法初值射线追踪和以费马原理为基础的弯曲法(扰动法)边值射线追踪.众多学者(Pereyra et al., 1980；Um and Thurber, 1987；Sambridge and Kennett, 1990；Farra, 1993；Mao and Stuart, 1997；徐涛等, 2004；李飞等, 2013；Xu and Mao, 2018)针对两点射线追踪问题提出了许多解决方法, 其出发点都是基于试射法和扰动法.基于这两种方法的射线追踪技术对于速度分界面能够用曲面函数进行描述的块状构造数值模拟具有优势, 能够较好地适用于反射波和多次波的快速正演, 较少应用于叠前深度偏移成像.
20世纪80年代以来, Kirchhoff积分法叠前深度偏移在复杂构造成像上发挥了重要的作用, 这使得网格类的射线追踪方法得到了快速发展.网格射线追踪主要包括基于图论和费马原理的最短路径算法, 基于求解程函方程的有限差分法和基于求解射线追踪方程组的波前类方法.Nakanishi和Yamaguchi(1986)提出了最短路径法网格射线追踪, Moser(1991)对最短路径算法做了全面的研究引起了人们的关注.随后, Fischer和Lees(1993)、Klime和Kvasni ka(1994)、刘洪等(1995)、王辉和常旭(2000)、Van Avendonk等(2001)、张美根等(2006), 从计算精度和计算效率出发, 充分发展和完善了最短路径算法.最短路径算法是一种无条件稳定的地震波射线追踪方法, 可以无盲区地求得每一个网格点上的初至路径和走时.但该方法通过网格节点之间的连线近似射线路径, 射线的出射角不能随入射角连续变化, 求得的路径难以满足斯奈尔定理, 走时计算的误差会随着传播长度的增加不断累积, 射线路径和走时计算的精度甚至低于传统算法, 且难以处理变速介质中存在的多重射线路径.
Vidale(1988, 1990)最早提出了程函方程的有限差分走时计算方法, 开拓了一条网格法射线追踪的新途径.该方法用矩形网格剖分速度场, 从震源开始通过矩形盒式扩展的方式计算走时, 模拟波前面的传播.在速度变化平缓的介质中有限差分法可以获得较精确的走时和射线路径, 当速度的非均匀性较强时会导致算法不稳定, 模拟的波前不准确, 复杂介质中射线的多路径问题也难以解决. Van Trier和Symes(1991)、Podvin和Lecomte(1991)、Qin等(1992)对该方法进行了不同程度的改进, 依然存在计算效率低的缺点.为了提高射线追踪的效率, 增强对复杂速度模型的适应性, Vinje等(1993)在二维光滑介质中采用波前构建法实现了射线路径、地震波初至走时和振幅的计算.波前构建法是一种基于求解射线追踪方程组(Červený et al., 1977)的波前扩展射线追踪方法, 该方法通过解射线追踪方程组获得初始射线路径和波前面, 然后在初始波前面上插入新的射线, 逐步递推, 模拟波前的扩展；在波前扩展的过程中根据多条射线上的走时, 插值计算规则网格点上的走时.Vinje等(1996)将波前构建法射线追踪由二维推广到三维, 为了提高运动学和动力学射线追踪方程组的求解精度, Vinje采用四阶泰勒展开代替原先的二阶泰勒展开, 提高了解的精度, 但是降低了计算效率.Ettrich和Gajewski(1996)采用四阶Runge-Kutta法求解射线追踪方程组, 这种方法同时具有很高的精度和计算效率，被一直沿用至今.波前构建法能够更好地描述非均匀介质中的地震波, 它能够获得射线路径, 可以处理复杂介质中的多路径, 计算多走时；但是在波前面上插入新的射线严重影响了计算效率, 在计算机实现上也存在一定困难, 尤其是在三维介质中.在求解射线追踪方程组的波前类射线追踪方法中, 基于单条射线的走时外插方法克服了波前构建法需要在波前面上插入射线而导致的计算效率低下和实现困难的缺陷, 能够更好地适应三维介质的叠前偏移成像.
为更好地进行复杂构造的高质量叠前偏移成像, 本文从射线追踪方程组的计算效率和走时外插的计算精度出发, 提出了一套快速的射线追踪算法和一种高精度的走时外插计算方法.本文提出采用基于线性多步法的预测-校正公式求解方程组, 给出了具体的实现过程, 与经典四阶Runge-Kutta公式相比, 有效地提高了计算效率.在走时的外插处理中, 通过射线路径上相邻的已知点确定圆台, 利用圆台的两个底面圆心上的走时计算规则网格节点上的走时, 与传统矩形体走时外插方法相比, 提高了计算精度和稳定性.根据圆台走时外插方法, 采用稀疏分布的射线也能够获得高精度的走时表, 节省了射线追踪的计算量, 对高精度的三维成像非常有利.最后通过对二维SEG/EAGE盐丘模型和三维复杂构造的逆散射偏移成像(Beylkin and Burridge, 1990；Mao et al., 2013；Ouyang et al., 2015；李武群等, 2017), 验证了该方法的有效性和适用性.
1 射线追踪方程组的高效求解在三维各向同性介质中, 地震波的传播路径满足以下运动学射线追踪方程：
|
(1) |
上式中, p =(px, py, pz)T, r =(x, y, z)T, t为地震波传播时间, v为传播速度；x, y, z是射线上的点r在直角坐标系下的三个分量, px, py, pz是慢度矢量p在直角坐标系下的三个分量.(1)式是一个一阶常微分方程组, 可写成以下分量形式：
|
(2) |
包含m个分量方程的一阶常微分方程组具有以下向量表达的一般形式
|
(3) |
式中, Y =(y1, y2, …, ym)T, Y是由m个分量函数构成的向量, F是其各分量随自变量τ的导数构成的向量, Y0为初始τ=τ0各分量函数构成的初值向量.在(3)式的数值求解方法当中, 四阶Runge-Kutta公式(以下简称RK)是高精度的单步法, 广泛应用于微分方程组的数值解.以线性多步法为基础的预测-校正公式(Hamming, 1959；孙志忠等, 2011)充分利用前面若干步的信息, 也能够获得与四阶Runge-Kutta同阶精度的解.此外, 线性多步法减少了所需计算偏导数的次数, 有效地提高了计算效率, 常用的预测-校正方法有Adams预测-校正公式, Milne-Hamming预测-校正公式.
求解(3)式的四阶RK公式为
|
(4) |
式中, n=0, 1, 2, …, h为离散自变量的步长, K1, K2, K3, K4满足下式
|
(5) |
求解(3)式的Adams预测-校正公式(以下简称Adams)为
|
(6) |
求解(3)式的Milne-Hamming预测-校正公式(以下简称MH)为
|
(7) |
(6) 式和(7)式中, Yn+1为Yn+1的预测值.预测-校正公式分为两个步骤, 本文以(6)式(Adams预测-校正公式)为例, 说明求解运动学射线追踪方程组(2)式的具体实现过程.首先是预测,
|
(8) |
上式中, ht是射线追踪时间步长, xk, yk, zk是第k步(k≥3)射线端点坐标的三个分量, (px)k, (py)k, (pz)k是第k步射线端点对应的慢度矢量的三个分量, 
|
(9) |
其中, 
|
(10) |
按照同样的过程, (6)式也可以求解动力学射线追踪方程组.(6)式充分运用了前面若干步的计算结果, 减少了逐步推进求解微分方程组的过程中所需计算偏导数的次数, 是一种高精度数值求解微分方程组的快速方法(孙志忠等, 2011).(7)式(MH公式)也是基于多步法的高精度快速求解公式, 实现过程与(6)式(Adams公式)类似.为更好地比较三种方法的计算结果和效率, 我们分别采用(4)式(RK公式), (6)式和(7)式求解(2)式.图 1给出了二维盐丘速度模型与三维盐丘速度模型(切片显示)及这两个模型下三种方法计算得出的射线路径；图 1a为二维模型, 图 1b为其对应的射线路径；图 1c为三维模型, 图 1d为其对应的射线路径.在图 1b与图 1d中, 黑色实线为RK公式的计算结果, 绿色的虚线为MH公式的计算结果, 红色的虚线为Adams公式的计算结果, 可以看出, (4)式, (6)式和(7)式的计算结果几乎完全一致.进一步, 在三维盐丘速度模型中, 指定初始入射角与射线追踪时间步长, 分别采用
|
图 1 盐丘速度模型与三种方法求得的射线路径 (a)二维盐丘速度模型；(b)二维射线路径；(c)三维盐丘速度模型；(d)三维射线路径. Fig. 1 Salt velocity models and the ray paths computed by three methods (a) 2D salt velocity model; (b) 2D ray paths; (c) 3D salt velocity model; (d) 3D ray paths. |
Adams预测-校正公式与四阶Runge-Kutta公式求解微分方程组获得两条射线路径, 表 1给出了6个不同的传播时间(长度), 两条射线的偏差百分比(射线的偏差与传播长度的比值)；从表 1中可以看到, 即使传播7000 m, 两条射线的偏差不超过0.005%, 在数值误差允许的范围内, 两种方法计算的结果精度相当.在图 1的二维盐丘速度模型中, 给出了计算20万个入射角的耗时；在图 2的三维盐丘速度模型中, 给出了计算501×501个入射角(入射倾角和方位角均为501个)的耗时.在相同的运行环境下, 三种方法的计算耗时见表 2.可以看出, 在计算结果一致的情况下, 基于线性多步法的预测-校正公式具有计算效率上的优势.
|
表 1 两条射线的偏差百分比 Table 1 The percentage deviation of two rays |
|
图 2 网格节点上走时外插示意图 (a)盒式外插；(b)圆台外插. Fig. 2 Sketch map of traveltime extrapolation at grids (a) Box extrapolation; (b) TC extrapolation. |
|
|
表 2 三种方法计算时间对比 Table 2 Comparison of computation time for three methods |
在采用多步公式完成了射线追踪方程组的逐步推进求解后, 就获得了射线路径上的各个时间采样点对应的空间位置.为更好地适应地震偏移成像处理, 需要将射线路径上的初至走时外插到规则网格节点上.传统的走时外插方法基于矩形体, 该方法以射线路径上的已知点为立方体(二维介质中则是正方形)中心, 根据炮点到已知点的射线路径长度l和比例因子ε确定立方体的半边长a/2(a=2εl), 然后利用已知点的走时来近似立方体内网格节点上的走时, 实现网格节点上走时的外插计算.图 2a是传统矩形体盒式走时外插方法(以下简称盒式外插)示意图, 盒式外插方法的精度和稳定性与比例因子ε有十分密切的关系.ε的选定依赖于射线的密度、走时表的网格间距大小和射线追踪的时间步长.在地震偏移成像中, 盒式外插方法需要确定十分合适的比例因子ε.一方面, 应当将射线路径上的走时信息外插到更多的网格节点上, 这就要求ε不能太小；另一方面, 要求网格节点上的走时尽量准确, 这就要求ε不能太大.盒式外插方法计算的走时受ε影响较大, 在三维偏移成像的应用中, 考虑到三维射线追踪的计算量, 盒式外插方法难以获得十分准确的走时表, 不利于高精度的三维偏移成像.
为提高网格节点上走时计算的准确度, 本文充分运用地震波波前与射线垂直的传播特性, 提出了一种基于圆台的走时外插计算方法.该方法按照射线路径上的已知点把路径切割成多个线段, 以线段为轴确定圆台, 通过圆台(二维介质中则是等腰梯形)上下底的两个圆心上的走时计算圆台内网格节点上的走时.图 2b是本文圆台走时外插方法(以下简称圆台外插)的示意图.下面将具体说明圆台外插的实现过程.
图 3中, 设A(x1, y1, z1)和B(x2, y2, z2)是三维空间中某条射线上相邻的两个已知点.地震波从震源传到A点的走时为iht, 传播长度为li, 从震源传到B点的旅行时为(i+1)ht, 传播长度为li+1.则线段AB能够确定一个上底半径为ri=εli、下底半径为ri+1=εli+1、高为
|
图 3 圆台走时外插计算示意图 Fig. 3 Sketch map of TC traveltime extrapolation |
设M(xp, yp, zp)是模型中的任意一个网格节点, 则M落在圆台内需要满足两个条件：一是M点必须夹在圆台上底面所在的平面和圆台下底面所在的平面之间；二是M点到直线AB的距离小于经过这个点的圆台横截面的半径(将其记为r).满足这两个条件的网格节点落在圆台内, 不满足这两个条件的网格节点则落在圆台外或圆台表面.
容易看出, 圆台上表面和下表面所在的平面方程分别为
|
(11) |
其中
|
(12) |
|
(13) |
则M(xp, yp, zp)夹在两个平行平面之间需要满足的条件为
|
(14) |
假设M满足条件(14), 并且将它在直线AB上的投影取为点N(xt, yt, zt).由空间几何的知识可知N满足以下方程
|
(15) |
其中
另一方面, 很容易得出N点的初至走时t为
|
(16) |
根据波前与射线垂直, 用N点的初至走时来近似圆台内网格节点M的初至走时.至此, 实现了射线上的初至走时外插到规则网格节点的计算.
3 偏移成像数值试验为更好地比较两种走时外插方法(盒式外插方法与圆台外插方法)的计算准确度, 本节将两种方法应用到二维和三维逆散射偏移成像数值试验.我们选用了相同的射线分布(射线的数量、各条射线的初始入射角、射线追踪时间步长和射线追踪方法均相同), 只改变由射线到网格节点、走时的外插计算方法.对两种外插方法获得的走时数据以相同的成像方法分别成像, 根据成像结果比较两种外插方法的计算准确度, 说明圆台外插方法的准确性.
3.1 二维模型偏移测试图 4是二维盐丘速度模型, 模型的横向(x方向)有649个采样点, 垂向(z方向)有300个采样点, 横向采样间隔为80 ft, 垂向采样间隔为40 ft.炮点数量为325, 炮间距为160 ft, 检波点数量为649, 检波点间距为80 ft, 地震数据通过波动方程正演获得.表 3是射线追踪采用的两组参数, θmin和θmax分别为初始入射倾角的最小值和最大值, nθ是均匀离散初始入射倾角的数量, dt是射线追踪的时间步长, ε是走时外插选定的比例因子.图 5a与图 5b是表 3的第一组参数下, 两种走时外插方法获得的偏移成像结果, 图 5a是盒式外插的偏移结果, 图 5b是圆台外插的偏移结果.可以看出, 在盒式外插的偏移结果中, 同相轴不清晰, 位置也不准确, 尤其是深部构造, 整个剖面噪声较大, 分辨率低, 这是网格节点上的走时计算不准确造成的, 说明盒式外插方法在表 3的第一组参数下不适用；相反地, 圆台外插方法获得了较准确的成像结果, 同相轴清晰, 归位准确, 整个剖面信噪比和分辨率都较高.提高盒式外插的走时表准确度需要增加射线追踪的计算成本(增加射线密度, 减小射线追踪时间步长, 减小比例因子ε的值).为此, 表 3给出了入射角加密后的第二组参数, 以及第二组参数下两种外插方法的偏移成像结果图 5c和图 5d.图 5c为盒式外插的偏移结果(黑色虚线标识了一条同相轴的准确位置), 图 5d是圆台外插的偏移结果(黑色虚线标识了一条同相轴的准确位置).比较图 5a与图 5c可以看出, 对于盒式外插方法, 相比表 3第一组参数下获得的成像结果, 第二组参数下获得的成像结果得到明显改善, 但是依然存在同相轴位置不够准确的现象.比较图 5c与图 5d可以看出, 在表 3的第二组参数下, 圆台外插方法获得的成像结果, 同相轴更清晰, 位置更准确.对比圆台外插方法在表 3的第一组参数下获得的成像结果(图 5b)与在第二组参数下获得的成像结果(图 5d), 可以看出, 两个结果非常一致.这说明圆台外插方法计算的走时受射线追踪参数的影响很小, 具有非常好的稳定性, 采用稀疏分布的射线即可获得高精度的走时表, 这对复杂模型的高精度成像是非常有利的.盒式外插方法计算的走时受射线追踪参数影响较大, 难以确定合适的参数, 即使增加计算成本也难以保证精度, 不利于高精度的地震成像.
|
图 4 二维SEG/EAGE盐丘速度模型 Fig. 4 2D SEG/EAGE salt velocity model |
|
图 5 表 3的两组射线追踪参数下, 两种走时外插方法的逆散射偏移成像结果 (a)第一组参数下盒式外插的成像结果；(b)第一组参数下圆台外插的成像结果；(c)第二组参数下盒式外插的成像结果；(d)第二组参数下圆台外插的成像结果；(c)与(d)中黑色的虚线标记了一条准确的同相轴. Fig. 5 Inverse scatting migration results of 2 traveltime extrapolation methods with 2 sets of ray tracing parameters in Table 3 (a) Migration result of box extrapolation with the first set of parameters; (b) Migration result of TC extrapolation with the first set of parameters; (c) Migration result of box extrapolation with the second set of parameters; (d) Migration result of TC extrapolation with the second set of parameters; An accurate interface location is marked with a black dashed line both in (c) and in (d). |
|
|
表 3 二维射线追踪参数 Table 3 2D ray tracing parameters |
图 6是三维速度模型与地表炮点分布, 模型的横向(x方向)与纵向(y方向)各有200个采样点, 垂向(z方向)有150个采样点, 三个方向的采样间隔均为20 m.一共9条炮线, 线间距为400 m, 每条炮线上17个炮点, 炮间距为200 m.每炮81×81道接收, 横向、纵向接收点间距均为40 m, 起始炮点与接收点的坐标均为(400 m, 400 m, 0 m), 地震数据通过波动方程正演获得.表 4是两组射线追踪参数, 起始倾角的最小值和最大值分别为θmin和θmax, 采样数量为nθ, 起始方位角的最小值和最大值分别为φmin和φmax, 采样数量为nφ, dt是射线追踪的时间步长, ε是走时外插选定的比例因子.图 7是目标偏移速度切片, 图 7a与图 7b分别是y=2000 m的速度切片与x=2000 m的速度切片.图 8为表 4的第一组参数下盒式外插方法的偏移结果(图 8中的白色虚线标识了同相轴的准确位置), 图 8a是y=2000 m切片的偏移结果, 图 8b是x=2000 m切片的偏移结果.图 9为表 4的第一组参数下圆台外插方法的偏移结果(图 9中黑色的虚线标识了同相轴的准确位置), 图 9a是y=2000 m切片的偏移结果, 图 9b是x=2000 m切片的偏移结果.可以看到, 两种外插方法获得的成像结果整体上比较接近, 细节上有所不同.通过图 8a与图 8b可以看出, 盒式外插方法获得的成像结果, 同相轴的位置与准确位置(图中白色虚线)有较小的偏离.通过图 9a与图 9b可以看出, 圆台外插方法获得的成像结果, 同相轴的位置与准确位置匹配得很好.相比盒式外插方法, 圆台外插方法更加符合地震波传播的运动学特征, 该方法在外插的过程中考虑了地震波的传播方向引起的走时变化, 使网格节点上求得的走时更加准确.对于圆台外插, 即使给定较大的外插比例因子, 也可以获得高精度的走时表.图 10是表 4的两组参数下, 圆台外插方法获得的两个成像结果, 其中图 10a是第一组参数下获得的偏移成像结果, 图 10b是第二组参数下获得的偏移成像结果, 成像切片位置为y=2000 m.可以看出, 在射线追踪参数变化较大的情况下, 圆台外插方法获得的偏移成像结果都具有很高的精度, 且一致性非常好.圆台外插方法容许射线追踪参数在较大的范围内变化, 是一种精度高、稳定性强的走时外插方法, 对于高精度的三维成像(尤其是深部构造)非常有利.
|
图 6 三维速度模型(a)与炮点分布(b) Fig. 6 3D velocity model (a) and shots distribution on the surface (b) |
|
|
表 4 三维射线追踪参数 Table 4 3D ray tracing parameters |
|
图 7 目标速度切片 (a) y=2000 m的速度切片；(b) x=2000 m的速度切片. Fig. 7 Targeted velocity slices (a) Velocity slice at y=2000 m; (b) Velocity slice at x=2000 m. |
|
图 8 表 4的第一组参数下盒式外插方法的逆散射偏移结果 (a) y=2000 m的偏移结果；(b) x=2000 m的偏移结果；白色的虚线标记了一条准确的同相轴. Fig. 8 Inverse scattering migration results of box extrapolation method with the first set of parameters in Table 4 (a) Migration results at y=2000 m; (b) Migration results at x=2000 m; An accurate interface location is marked with a white dashed line. |
|
图 9 表 4的第一组参数下圆台外插方法的逆散射偏移结果 (a) y=2000 m的偏移结果；(b) x=2000 m的偏移结果；黑色的虚线标记了一条准确的同相轴. Fig. 9 Inverse scattering migration results of TC extrapolation method with the first set of parameters in Table 4 (a) Migration results at y=2000 m; (b) Migration results at x=2000 m.An accurate interface location is marked with a black dashed line. |
|
图 10 表 4的两组参数下, 圆台外插方法在y=2000 m处的逆散射偏移结果 (a)第一组参数的偏移结果；(b)第二组参数的偏移结果. Fig. 10 Inverse scattering migration results of TC extrapolation method with two sets of parameters in Table 4 at y=2000 m (a) Migration result with the first set of parameters; (b) Migration result with the second set of parameters. |
为提高射线追踪方程组的计算效率和成像网格节点上初至走时的计算精度, 使其更好地适应复杂构造的偏移成像, 本文提出了基于线性多步法的快速射线追踪方法和基于圆台的高精度走时外插计算方法.鉴于线性多步法预测-校正公式能充分利用多步信息, 减少了速度偏导数的计算次数, 本文将其应用到射线追踪方程组的求解中, 提高了射线追踪的计算速度.在网格节点的走时表计算中, 本文充分考虑地震波波前与射线垂直的传播特性, 将射线切割成多个线段, 并以每个线段为轴确定圆台, 然后将圆台轴上已知的走时外插到圆台内的网格节点上；由此, 可明显提高走时外插计算精度.这种圆台外插算法受射线追踪参数的影响较小, 在较大的参数变化范围内都能够获得高精度的走时表, 具有很好的稳定性.圆台走时外插方法在射线分布稀疏的情况下依然可获得高精度的走时表, 克服了传统矩形体走时外插方法对射线追踪参数敏感以及难以选择合理射线追踪参数的缺陷, 对高精度的三维地震偏移成像非常有利.最后通过对二维和三维模型的逆散射偏移成像数值试验, 验证了算法的有效性和适用性.
Beylkin G, Burridge R. 1990. Linearized inverse scattering problems in acoustics and elasticity. Wave Motion, 12(1): 15-52. DOI:10.1016/0165-2125(90)90017-X |
Červený V, Molotkov I A, Pšenčík I. 1977. Ray Method in Seismology. Praha: Univerzita Karlova.
|
Ettrich N, Gajewski D. 1996. Wave front construction in smooth media for prestack depth migration. Pure and Applied Geophysics, 148(3-4): 481-502. DOI:10.1007/BF00874576 |
Farra V. 1993. Ray tracing in complex media. J. Appl. Geophys., 30(1-2): 55-73. DOI:10.1016/0926-9851(93)90018-T |
Fischer R, Lees J M. 1993. Shortest path ray tracing with sparse graphs. Geophysics, 58(7): 987-996. DOI:10.1190/1.1443489 |
Hamming R W. 1959. Stable predictor-corrector methods for ordinary differential equations. Journal of the ACM, 6(1): 37-47. |
Klimeš L, Kvasnička M. 1994. 3-D network ray tracing. Geophys. J. Int., 116(3): 726-738. |
Li F, Xu T, Wu Z B, et al. 2013. Segmentally iterative ray tracing in 3-D heterogeneous geological models. Chinese Journal of Geophysics (in Chinese), 56(10): 3514-3522. DOI:10.6038/cjg20131026 |
Li W Q, Mao W J, Ouyang W, et al. 2017. Second-order Born approximation and its application to inverse scattering GRT inversion. Progress in Geophysics (in Chinese), 32(2): 645-656. DOI:10.6038/pg20170227 |
Liu H, Meng F L, Li Y M. 1995. The interface grid method for seeking global minimum travel-time and the correspondent raypath. Chinese Journal of Geophysics (in Chinese), 38(6): 823-832. |
Mao W J, Stuart G W. 1997. Rapid multi-wave-type ray tracing in complex 2-D and 3-D isotropic media. Geophysics, 62(1): 298-308. DOI:10.1190/1.1444131 |
Mao W J, Ouyang W, Li X L. 2013. Nonlinear migration inversion with second-order Born approximation.//75th EAGE Annual International Conference and Exhibition. Extended Abstracts, Tu-P06-13.
|
Moser T J. 1991. Shortest path calculation of seismic rays. Geophysics, 56(1): 59-67. |
Nakanishi I, Yamaguchi K. 1986. A numerical experiment on nonlinear image reconstruction from first-arrival times for two-dimensional island arc structure. J. Phys. Earth., 34(2): 195-201. DOI:10.4294/jpe1952.34.195 |
Ouyang W, Mao W J, Li X L. 2015. Approximate solution of nonlinear double-scattering inversion for true amplitude imaging. Geophysics, 80(1): R43-R55. |
Pereyra V, Lee W H K, Keller H B. 1980. Solving two-point seismic-ray tracing problems in a heterogeneous medium:Part 1. Bull. Seismol. Soc. Am., 70(1): 79-99. |
Podvin P, Lecomte I. 1991. Finite difference computation of traveltimes in very contrasted velocity models:A massively parallel approach and its associated tools. Geophys. J. Int., 105(1): 271-284. |
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 |
Sambridge M S, Kennett B L N. 1990. Boundary value ray tracing in a heterogeneous medium:a simple and versatile algorithm. Geophys. J. Int., 101(1): 157-168. |
Sun Z Z, Wu H W, Yuan W P, et al. 2011. An Elementary Numerical Analysis (in Chinese). 5th ed. Nanjing: Southeast University Press.
|
Um J, Thurber C. 1987. A fast algorithm for two-point seismic ray tracing. Bull. Seismol. Soc. Am., 77(3): 972-986. |
Van Trier J, Symes W W. 1991. Upwind finite-difference calculation of traveltimes. Geophysics, 56(6): 812-821. DOI:10.1190/1.1443099 |
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. Geophysics, 66(2): 648-653. |
Vidale J E. 1988. Finite-difference calculation of travel times. Bull. Seismol. Soc. Am., 78(6): 2062-2076. |
Vidale J E. 1990. Finite-difference calculation of traveltimes in three dimensions. Geophysics, 55(5): 521-526. DOI:10.1190/1.1442863 |
Vinje V, Iversen E, Gjøystdal H. 1993. Traveltime and amplitude estimation using wavefront construction. Geophysics, 58(8): 1157-1166. DOI:10.1190/1.1443499 |
Vinje V, Iversen E, Åstebøl K, et al. 1996a. Estimation of multivalued arrivals in 3D models using wavefront construction-Part I. Geophysical Prospecting, 44(5): 819-842. DOI:10.1111/j.1365-2478.1996.tb00175.x |
Wang H, Chang X. 2000. 3-D ray tracing method based on graphic structure. Chinese Journal of Geophysics (in Chinese), 43(4): 534-541. |
Xu Q R, Mao W J. 2018. An efficient ray-tracing method and its application to Gaussian beam migration in complex multilayered anisotropic media. Geophysics, 83(5): T281-T289. DOI:10.1190/geo2017-0402.1 |
Xu T, Xu G M, Gao E G, et al. 2004. Block modeling and shooting ray tracing in complex 3-D media. Chinese Journal of Geophysics (in Chinese), 47(6): 1118-1126. |
Zhang M G, Cheng B J, Li X F, et al. 2006. A fast algorithm of shortest path ray tracing. Chinese Journal of Geophysics (in Chinese), 49(5): 1467-1474. |
李飞, 徐涛, 武振波, 等. 2013. 三维非均匀地质模型中的逐段迭代射线追踪. 地球物理学报, 56(10): 3514-3522. DOI:10.6038/cjg20131026 |
李武群, 毛伟建, 欧阳威, 等. 2017. 二阶Born近似及其在逆散射GRT反演中的应用. 地球物理学进展, 32(2): 645-656. DOI:10.6038/pg20170227 |
刘洪, 孟凡林, 李幼铭. 1995. 计算最小走时和射线路径的界面网全局方法. 地球物理学报, 38(6): 823-832. DOI:10.3321/j.issn:0001-5733.1995.06.015 |
孙志忠, 吴宏伟, 袁慰平, 等. 2011. 计算方法与实习. 5版. 南京: 东南大学出版社.
|
王辉, 常旭. 2000. 基于图形结构的三维射线追踪方法. 地球物理学报, 43(4): 534-541. DOI:10.3321/j.issn:0001-5733.2000.04.014 |
徐涛, 徐果明, 高尔根, 等. 2004. 三维复杂介质的块状建模和试射射线追踪. 地球物理学报, 47(6): 1118-1126. DOI:10.3321/j.issn:0001-5733.2004.06.027 |
张美根, 程冰洁, 李小凡, 等. 2006. 一种最短路径射线追踪的快速算法. 地球物理学报, 49(5): 1467-1474. DOI:10.3321/j.issn:0001-5733.2006.05.026 |