2012, Vol. 55
2. 国土资源部应用地球物理综合解释理论开放实验室—波动理论与成像技术实验室, 长春 130026
2. Laboratory for Integrated Geophysical Interpretation Theory of the Ministry for Land and Resources of China-Laboratory for Wave Theory and Imaging Technology, Changchun 130026,China
我国大量的地震勘探工作在具有复杂地形的区域开展，成像效果差或基本不成像是复杂地形地区地震数据处理常会遇到的一个难题.其中一个很重要的因素是现有的大多数地震波理论和地震成像技术均是基于水平地表的，并且静校正等技术在面对很多复杂地形时难以取得好的效果.所以对直接基于复杂地形的地震成像技术进行研究是很有必要的[1].地震波走时在偏移成像[2-5]､层析成像[6-7]等地震数据处理技术中是一个非常重要的物理量，所以对复杂地形条件下的地震走时计算方法进行研究具有一定的理论意义和实际价值[1].
基于水平地表的走时计算方法主要有两点射线追踪法(包括试射法和弯曲法)[8-9]､波前构建法[10-11]､有限差分法[12-16]､插值法[17-18]､最短路径射线追踪法[19-21]等.在这些算法中，试射法､弯曲法和最短路径射线追踪法有着不错的计算精度但是其计算效率相对较低.有限差分法由于各种高效的波前扩展方式[13-14]的引入而具有相对较高的计算效率，同时各种差分格式[15]的引入在很大程度上也能够保证算法的稳定性和计算精度.线性插值法的局部走时计算公式的原理简单，计算精度相对较高且适应于任意的复杂网格，所以常被引入到有限差分方法[16]中作为局部走时计算公式，但是纯粹的插值法[17]的扩展方式是不满足地震波传播规律的，并且计算效率也相对不理想.综上所述，这些算法各有千秋，但它们却有一个共同点:均在水平地表条件下被提出，且均需要做出一些改进才能用于解决复杂地形问题.
与水平地表问题不同，复杂地形问题需要算法在保证地震波传播规律和计算精度的前提下，还要对地形带来的不规则计算边界问题做灵活､稳定的处理.因此算法需要考虑以下两方面问题:①算法在不规则边界处必须采用稳定的局部走时计算公式；②算法的整体实现策略能够灵活地适应不规则边界.关于复杂地形条件下的地震波走时计算问题，国内外学者已经做过一些相关研究[22-34]，例如:Sun等基于对水平地形条件下的快速推进法和线性插值法的改进，分别提出了不等距网格有限差分法[31-32]和混合网格线性插值法[33-34].在这两种方法中，分别采用不等距网格和混合网格对复杂地形附近的介质进行剖分，这两种网格剖分方式虽然能细致地刻画地表的起伏形态，但是也会使得算法的计算公式和实现过程比水平地表的情况要复杂很多.针对这种情况，本文首先将常规的阶梯网格和迎风差分格式组合在一起，提出了一种基于阶梯网格的算法，即阶梯网格迎风差分法.然后将这种算法与不等距网格有限差分法和混合网格线性插值法作对比，并分析了三种算法各自的特点.最后利用一个计算实例来检验算法对复杂地形､近地表及地下复杂介质等复杂地质条件的适应性和稳定性.
2 复杂地形条件下的网格剖分及局部走时算法复杂地形条件下走时计算的核心是对地表起伏的不规则边界问题的处理.该问题实际上就是对含不规则边界计算区域采用怎样的网格剖分方式，以及采用怎样与之相对应的局部走时计算公式.针对复杂地形的网格剖分问题，最简洁的方法是采用阶梯状的正方形网格，因此，将首先提出一种基于阶梯网格的算法，该算法采用改进后的迎风差分格式作为局部走时计算公式.然后再简要地回顾不等距网格有限差分法和混合网格线性插值法的网格剖分方式和局部走时计算公式.
2.1 阶梯网格迎风差分法常规的地震波走时算法一般采用正方形网格剖分速度模型，然而在面对复杂地形时，起伏地表上真实的高程点(如图 1a中用“￭"表示的节点)一般不能刚好落在规则网格节点处.此时采用地表以下且离地表真实的高程点最近的规则网格节点(如图 1a中用“●"表示的节点)来近似代替地表高程点，这样地表的形态就被刻画为用规则正方形网格组成的阶梯状的网格.下面将建立阶梯网格中的局部走时计算公式.
|
图 1 阶梯网格迎风差分法的模型剖分及局部走时算法 (a)模型的网格剖分；(b)局部走时算法. Fig. 1 The model discretization and local traveltimecomputation scheme of upwind finite difference method with stepwise grid (a) The grids for model discretization； (b) The scheme for local traveltime computation. |
常规的用于近似Eikonal方程中走时函数梯度项的迎风差分格式[13]可表述为
|
(1) |
其中，▽为nabla算子，vi，j､Ti，j分别为网格节点(i，j)处的速度值､走时值，Dn+xTi，j､Dn-xTi，j､Dn+zTi，j､Dn-zTi，j分别为网格节点(i，j)处x和z方向上的走时函数T的n阶向前､向后差分算子，在正方形网格中以一阶差分格式为例有:
|
(2) |
其中，h为网格间距.
常规迎风差分格式(1)在计算点(i，j)的局部走时时，需要通过判断与其最邻近的四个网格节点的走时值的已知情况来选取差分格式的方向[13]，如图 1b中的情况1.
同理可知max(Dn-zTi，j，-Dn-zTi，j，0).当然x､z方向上的差分格式的选择不可能同时都为零(若都为零，(1)式不成立，无法计算)，因为被算点的其中一个邻居点一定是窄带技术(后续内容将详细阐述)选取的局部最小走时点，而该点的走时值必为已知.最后，将两个方向上选取的差分格式代入(1)式中即可计算点(i，j)的走时值.
然而当面对复杂地形问题时，即在阶梯网格中采用迎风差分时，被算点周围的部分网格节点将出现与图 1b中的情况1不一致的网格分布情况，即可能出现如图 1b中的情况2~4，被算点周围的部分网格节点(图 1b中“×"标注的节点)根本不在计算区域内.为了解决该问题，在此采用如下处理方式:
|
(3) |
直接把那些不在计算区域内的网格节点的走时值视为公式(3)中的“未知"的情况.该处理方式的结果是在走时计算时永远不会用到那些不在计算区域内的点.而采用这种处理方式的理由:迎风差分法在计算某点的走时值时，实际上用到的是与该点相邻且走时值已知的网格节点(若同一个方向的两个相邻点的走时值均已知，则选择走时值较小的网格节点)的走时信息，这样做是遵循波前的传播规律的(波前的走时信息总是从走时值小的位置传向走时值大的位置).而那些如图 1b的2~4情况中“×"标注的节点均位于空气中，因为地震波不进入空气中，所以这些点均不含已被算过的走时信息，即说明这种永远“未知"的处理方式是符合物理事实的.
2.2 不等距网格有限差分法如图 2a所示，不等距网格有限差分法在对模型进行网格剖分时直接保留地表上真实高程点(如图 2a中用“￭"标注的点)的准确位置.这样的网格剖分方式会在地表的不同位置产生不同的网格间距，同时也会给局部走时计算带来一定的困难，因为此时地表上的节点的走时计算不能采用像规则网格节点一样的迎风差分格式.为了解决该问题，在迎风差分格式中引入不等距有限差分格式[31-32].
|
图 2 不等距有限差分法的网格剖分及局部走时算法 (a)模型的网格剖分；(b)局部走时算法. Fig. 2 The model discretization and local traveltime computation scheme of finite difference methodwith nonuniform grid spacing (a) The grids for model discretization； (b) The scheme for local traveltime computation. |
如图 2a所示，根据地表高程点所处位置地形情况的不同，可将局部不等距网格的类型划分为三类，即图 2b中的情况1~3(情况1~3分别表示与地表高程点X最邻近的地表以下的点A的左右两个网格节点均在地表以下､一个在地表以下､均不在地表以下的情况).下面将分别讨论这三种情况下X点和A点的局部走时算法.
2.2.1 计算地表上点X的走时值当网格为图 2b中的情况1 时，此时由于X点位于地表上，无法构建如图 1b情况1所示的差分格式(2)，此时在A点的x方向上采用(3)式的差分方法，而在z方向采用不等距差分格式:
|
(4) |
其中，TA､TB､TC､TX分别为A､B､C､X点的走时值，hAX为线段AX的长，SA为点A处介质的平均慢度值.根据因果关系TX> TA(已完成走时计算的网格节点的走时值不大于当前正进行走时计算的网格节点的走时值)有:
|
(5) |
其中，θAX=hAX/h.
当网格为图 2b中的情况2和3时，A点处x方向上用“×"标注的节点的处理方式可由图 1b中的情况2~4同理而得.公式(5)实际上近似考虑了地震波由A点传向X点时的情况，若考虑地震波从横向方向传入时，即当相邻地表上的节点(图 2b中的点L和点R)的走时值为已知时，此时基于Huygens原理和Fermat原理有:
|
(6) |
TL､TR分别为点L､R的走时值，dLX､dRX分别为点L､R到X点的距离，SX为点X处介质的平均慢度值.
2.2.2 计算与地表点最邻近的A点的走时值当网格为情况1时，由于在z方向上max((TA-TX)/hAX，- (TD-TA)/h，0)因为hAX≠h而无法确定，此时基于Fermat原理两次采用公式:
|
(7) |
其中，一次Ti=Tx､hi=hAX；另外一次Ti= TD､hi=h，然后求解二次方程(7)并选取两次结果中较小者作为TA最终的解.
综上所述，不等距网格有限差分法在进行整个复杂地形模型的走时计算时，当被计算网格节点为规则网格节点则采用(3)式的差分格式，当被计算点在地表附近时则分别对应地采用(5)~(7)式的差分格式.
2.3 混合网格线性插值法如图 3a所示，在对模型进行网格剖分时，混合网格线性插值法采用三角网格和正方形网格相结合的混合网格描述地表形状，而在局部走时计算时则采用线性插值法[33-34].
|
图 3 混合网格线性插值法的网格剖分及局部走时算法 (a)模型的网格剖分；(b-d)局部走时算法. Fig. 3 The model discretization and local traveltime computation scheme of linear interpolation method with hybrid grid (a)The grids for model discretization； (b-d)The scheme for local traveltime computation. |
如图 3b所示，常规线性插值法是在正方形网格中推导局部走时计算公式的[17]，然而在面对复杂地形问题时需要推导三角网格中的局部走时计算公式，即如图 3c所示，设点A､B的走时为TA､TB 且为已知，现求点C的走时值TC.根据线性插值法的假设，线段v上的走时值满足线性分布且到达点C的射线经过线段AB上的一点E.根据费马原理E点是线段AB上使得TC 取最小值的那一点，从而若要求取TC，则首先需要确定E点在线段AB上的位置.为此首先设点E为线段AB上的任意一点，经C点向线段AB引垂线交线段AB于点D(垂足)，并以点D为坐标原点建立局部直角坐标系，该坐标系的x轴与线段AB重合，y轴与线段DC重合(如图 3.d 所示).设线段AB､BD､DC､EC的长分别为d1､d2､d3､d4，E在局部坐标系中的坐标为(x，0).若以B点为参考点，则通过线性插值得:
|
(8) |
其中，ΔT=TA-TB，-(d1-d2)≤x≤d2，因此无论x>0(E点位于线段BD之间)还是x＜0(E点位于线段DA之间)，d2-x均表示线段BE的长.利用公式(8)可得到TC的表达式为
|
(9) |
其中:S为三角形ABC所在区域的平均慢度.
如上所述，公式(9)对线段AB上的任意一点均成立.为了找到使TC取最小值(即C点的最小走时)的x值(E点的坐标)，根据费马原理求式(9)关于x的最小值，即求解方程$\partial {{T}_{C}}/\partial x=0$.将求出的x值带回(9)式中便可得三角网格中计算C点走时值TC的局部走时计算公式[33-34]:
|
(10) |
如图 3b所示，同理可得在规则的正方形网格中的局部走时计算公式[33-34]:
|
(11) |
综上所述，线性插值法在进行局部走时计算时，当被计算网格节点为规则网格节点时，采用(11)式的计算公式，当被计算点为地表上的点时则采用(10)式的计算公式.
3 复杂地形条件下的波前扩展方式为了灵活稳定地适应不规则计算边界，选取窄带技术作为复杂地形条件下的波前扩展方式[31-34].常规的窄带技术是一种满足波前传播规律的波前扩展方式，它采用窄带来近似模拟波前，通过窄带的扩展演化来模拟地震波在介质中的传播过程[13].但是常规的窄带技术也是基于水平地表提出的，为了将其应用到复杂地形条件下，需要在窄带技术中引入新的网格节点类型.
如图 4所示，若地表这一不规则计算边界不存在，常规窄带技术将计算区域内的网格节点分为三种类型:接受点､窄带点､远离点.图中黑色点为接受点，表示走时值已经计算完成的网格节点；灰色点为窄带点，表示窄带内的网格节点，即近似波前上的网格节点；白色点为远离点，表示走时计算还未完成的点.窄带技术的基本实现策略是:通过窄带内网格节点的纳入和去除来实现窄带的扩展，进而达到近似模拟波前传播的目的(常规窄带技术的实现步骤详见参考文献[13]).
为了让窄带技术适应不规则计算边界，在此引入两种新的网格节点类型:地表点(图 4中用正方形表示的网格节点)和地表上点(图 4中用“×"表示的网格节点).对常规的窄带技术的实现步骤进行改进可得出复杂地形条件下窄带技术的实现步骤:
|
图 4 复杂地形条件下的窄带技术 Fig. 4 The narrow band technique in complex topographical conditions |
① 初始化:震源点为惟一的接受点，走时值为零；与震源点直接邻近的网格节点设为窄带点，它们构成初始窄带，其走时值通过解析公式求得；其余网格节点设为远离点，其走时值设为足够大(远大于可能会算出的最大的走时值).
②迭代计算:步骤一从窄带内选取最小走时点，并将其属性由窄带点改为接受点(将该点从窄带中移除)；步骤二对选出的最小走时点周围的网格节点做出判断，若其为远离点，则将其属性改为窄带点(纳入窄带中)，并利用(1)~(11)式中对应的公式计算其走时值；若其为窄带点，则利用(1)~(11)式中对应的公式更新其走时值(再计算一遍其走时值，并与原来的走时值做比较，若新计算出的走时值较小，则用其替换原来的走时值，否则保持原来的走时值不变)；若其为接受点和地表上点，则保持网格节点原有属性不变；步骤三判断窄带是否为空，若是则计算结束，否则跳回步骤一继续重复计算.
4 算法分析及计算实例 4.1 算法的计算精度对比为了研究上述三种算法在复杂地形条件下的计算精度，下面给出一个山峰模型.如图 5 所示，模型大小为8.0km×4.0km, 地震波在介质中的传播速度为1.0km/s, 计算时采用的网格间距为10.0m, 震源点的位置为(4.0km, 0.0km).图 5a显示了分别采用三种算法计算出的地表上的走时值的相对误差分布，从这个误差图可以得出:三种算法的误差均是在震源处较大，而离震源较远的区域较小.各算法之间，混合网格线性插值法的计算精度最高，阶梯网格迎风差分法的计算精度最低并且其相对误差曲线呈锯齿振荡状，而其他两种算法除在震源附近误差分布有振荡外，其余部分变化相对平稳.图 5(b~d)，分别显示了混合网格线性插值法､不等距网格有限差分法､阶梯网格迎风差分法在地表以下区域走时计算的相对误差分布.从图中可以发现整个计算区域三种算法均能保证良好的计算精度(绝大部分区域的相对误差均低于1.0%)，而混合网格线性插值法的计算精度最高.
|
图 5 复杂地形条件下三种走时算法的精度对比 (a)地表上三种算法走时计算的相对误差分布对比；(b)〜(d)混合网格线性插值法､不等距网格有限差分法和阶梯网格迎风差分法地表以下区域相对误差分布. Fig. 5 The computational accuracy comparison of three traveltime computation methods in complex topographical conditions (a) The comparison on the earth's surface of threemethods； (b)〜(d) The relative error distributing in the areas under the earth's surface of linearinterpolation method with hybrid grid, finite difference method with nonuniform grid spacing, and upwind finite difference method with stepwise grid. |
为了对比三种算法的计算效率，采用图 5所示的山峰模型，同时为了产生不同的计算量分别采用网格间距20､10､5､2m 剖分速度模型.表 1显示了三种算法在面对不同计算量时的CPU耗时，通过分析可以发现:阶梯网格的计算效率最高，其次是不等距网格有限差分法､计算效率最差的是混合网格线性插值法.
|
表 1 三种算法的计算效率(CPU耗时/s)对比 Table 1 The efficiency (CPU time/s) comparison of the three methods |
上述对比研究了三种可行的计算复杂地形条件下地震波走时的方法，这三种算法均都包含2个核心内容:模型的网格剖分和局部走时算法.
在模型的网格剖分方面，三种方法分别采用了阶梯网格､不等距网格和混合网格.从地表形态描述的准确性方面讲，阶梯网格最为粗糙；不等距网格准确的定位了地表的高程点；而混合网格法在准确定位地表高程点的同时，还通过三角网格建立了地表点之间以及地表点与规则网格节点间的位置关系，因此混合网格的方法更好地刻画了地表形态.不过混合网格却比其他两种方法更为复杂，同时也会给计算带来一定程度的复杂性.所以选取怎样的网格剖分模型，应该综合考虑算法的精度要求和地形的特点，比如说当地表高程变化较小时可以考虑采用较简单的阶梯网格剖分模型，而地表起伏很剧烈时应该考虑采用更为精细的混合网格剖分模型.
在局部走时计算方法方面，三种算法分别采用了迎风差分法､不等距差分法､线性插值法.当地表为水平时，前面两种方法实际上是一致的.在这几种局部走时算法中，迎风差分法是最简洁的，其次是不等距网格差分法，线性插值法由于需要同时考虑三角网格和正方形网格中的插值公式而更加复杂一些.三种算法均达到良好的计算精度，其中线性插值法的计算精度相对是最高的，其次是不等距网格有限差分法，最差的是阶梯网格迎风差分法.
从总体上讲，混合网格线性插值法较好地刻画了地表形态和采用了更为复杂的局部走时计算公式，所以获得了相对较高的计算精度，但是算法相应的复杂程度也提高了，所以其计算效率最低；阶梯网格迎风差分法的计算精度虽然最低，但其计算效率最高；而不等距网格有限差分法的计算精度和计算效率均居中.
4.4 计算实例为了验证三种算法在处理复杂地形问题的有效性，下面将给出一个经过复杂地形修正后的Marmousi速度模型(如图 6 所示)，该模型大小为8.0km×12.0km, 震源位于(6.0km, 1.0km)，计算时采用的网格间距为10.0m.通过观察可以发现该模型在地表处，除了地形剧烈起伏外地表附近的介质也很复杂，同时地下介质结构也很复杂.同样三种算法的计算结果基本吻合，算法只是在一些地表处和介质速度变化较大的区域有微小的差别.算法的计算结果说明三种算法在处理地表的剧烈起伏和近地表及地下复杂介质时均有较好的适应性和稳定性.
|
图 6 经过夏杂地形修正后的Marmous丨速度模型的等时线分布 (a)不等距网格有限差分法(白色)和混合网格线性插值法(黑色)；(b)阶梯网格迎风差分法(白色)和混合网格线性插值法(黑色). Fig. 6 The traveltime contours are superimposed on the Marmousi velocity model including strong surface topography (a) Finite difference method with nonuniform grid spacing (white) and linear interpolation scheme with hybrid grid (black)；(b) Upwind finite diference method with stepwise grid (white) and linear rnterpolation scheme with hybrid grid (black). |
复杂地形条件下的地震波走时算法对于研究直接基于复杂地形的地震偏移成像技术意义重大.为了获得精度高且适应于复杂地形的走时算法，提出了一种阶梯网格迎风差分法，并将该方法与不等距网格有限差分法和混合网格线性插值法进行对比研究.相关数值分析表明:(1)混合网格线性插值法计算精度最高，其次是不等距网格有限差分法和阶梯网格迎风差分法；(2)阶梯网格迎风差分法的计算效率最高，其次是不等距网格有限差分法和混合网格线性插值法；(3)究竟选取哪种算法作为给定复杂地形模型的地震波走时算法，应该综合考虑地形的特点､研究对象的计算精度及计算效率的要求等因素；(4)三种算法在面对复杂地形､近地表及地下复杂介质等复杂地质条件时均能保证很好的适应性和稳定性.
| [1] | 孙建国. 复杂地表条件下地球物理场数值模拟方法评述. 世界地质 , 2007, 26(3): 345–362. Sun J G. Methods for numerical modeling of geophysical fields under complex topographical conditions: a critical review. Global Geology (in Chinese) , 2007, 26(3): 345-362. |
| [2] | Sun J G. True-amplitude weight functions in 3D limited-aperture migration revisited. Geophysics , 2004, 69(4): 1025-1036. DOI:10.1190/1.1778245 |
| [3] | Sun J G. Limited-aperture migration. Geophysics , 2000, 65(2): 584-595. DOI:10.1190/1.1444754 |
| [4] | Sun J G. On the aperture effect in 3D Kirchhoff-type migration. Geophysical Prospecting , 1999, 47(6): 1045-1076. DOI:10.1046/j.1365-2478.1999.00152.x |
| [5] | Sun J G. On the limited aperture migration in two dimensions. Geophysics , 1998, 63(3): 984-994. DOI:10.1190/1.1444409 |
| [6] | 井西利, 杨长春, 王世清. 一种改进的地震反射层析成像方法. 地球物理学报 , 2007, 50(6): 1831–1836. Jing X L, Yang C C, Wang S Q. A improved seismic reflection tomographic method. Chinese J. Geophys. (in Chinese) , 2007, 50(6): 1831-1836. |
| [7] | 瞿辰, 周蕙兰, 赵大鹏. 使用纵波和横波走时层析成像研究菲律宾海板块西边缘带和南海地区的深部结构. 地球物理学报 , 2007, 60(6): 1757–1768. Qu C, Zhou H L, Zhao D P. Deep structure beneath the west margin of Philippine Sea Plate and South China Sea from P and S wave travel time tomography. Chinese J. Geophys. (in Chinese) , 2007, 60(6): 1757-1768. |
| [8] | Julian B R, Gubbins D. Three-dimensional seismic ray tracing. J. Geophys. , 1977, 43(1-2): 95-113. |
| [9] | 徐涛, 徐果明, 高尔根, 等. 三维复杂介质的块状建模和试射射线追踪. 地球物理学报 , 2004, 47(6): 1118–1126. Xu T, Xu G M, Gao E G, et al. Block modeling and shooting ray tracing in complex 3-D media. Chinese J. Geophys. (in Chinese) , 2004, 47(6): 1118-1126. |
| [10] | Vinje V, Iversen E, Gjystdal H. Traveltime and amplitude estimation using wavefront construction. Geophysics , 1993, 58(8): 1157-1166. DOI:10.1190/1.1443499 |
| [11] | Han F X, Sun J G, Sun Z Q, et al. Positioning of grid points in wave front construction. Applied Geophysics , 2009, 6(3): 248-258. DOI:10.1007/s11770-009-0024-z |
| [12] | Vidale J E. Finite-difference calculation of travel times. Bull. Seism. Soc. Am. , 1988, 78(6): 2062-2076. |
| [13] | Sethian J A, Popovici A M. 3-D travel time computation using the fast marching method. Geophysics , 1999, 64(2): 516-523. DOI:10.1190/1.1444558 |
| [14] | Rawlinson N, Sambridge M. Wave front evolution in strongly heterogeneous layered media using the fast marching method. Geophys. J. Int. , 2004, 156(3): 631-647. DOI:10.1111/gji.2004.156.issue-3 |
| [15] | Sun J G, Yang H, Han F X. A finite-difference scheme for solving the eikonal equation with varying grid spacing. 77th Ann. Internat. Mtg., Soc. Expl. Geophys., Expanded Abstracts, 2007: 2120-2124. |
| [16] | Schneider W A Jr, Ranzinger K A, Balch A H, et al. A dynamic programming approach to first arrival traveltime computation in media with arbitrarily distributed velocities. Geophysics , 1992, 57(1): 39-50. DOI:10.1190/1.1443187 |
| [17] | Asakawa E, Kawanaka T. Seismic ray tracing using linear traveltime interpolation. Geophysical Prospecting , 1993, 41(1): 99-111. DOI:10.1111/gpr.1993.41.issue-1 |
| [18] | Wang H Z, Ma Z T, Hubral P. Traveltime calculation in 3D media with arbitrary velocity distribution. 69th Ann. Internat. Mtg., Soc. Expl. Geophys., Expanded Abstracts 1999: 1778-1781. |
| [19] | Moser T J. Shortest path calculation of seismic rays. Geophysics , 1991, 56(1): 59-67. DOI:10.1190/1.1442958 |
| [20] | Cao S H, Greenhalgh S. Calculation of the seismic first-break time field and its ray path distribution using a minimum traveltime tree algorithm. Geophys. J. Int. , 1993, 114(3): 593-600. DOI:10.1111/gji.1993.114.issue-3 |
| [21] | 张美根, 贾豫葛, 王妙月, 等. 界面二次源波前扩展法全局最小走时射线追踪技术. 地球物理学报 , 2006, 49(4): 1169–1175. Zhang M G, Jia Y G, Wang M Y, et al. A global minimum traveltime raytracing algorithm of wavefront expanding with interface points as secondary sources. Chinese J. Geophys. (in Chinese) , 2006, 49(4): 1169-1175. |
| [22] | 李满树, 方伍宝, 周腾, 等. 初至波走时信息在复杂地区近地表速度反演中的应用. 石油物探 , 2004, 43(1): 72–75. Li M S, Fang W B, Zhou T, et al. Inversion of near-surface velocity through first-arrival traveltime in complex region. Oil Geophysical Prospecting (in Chinese) , 2004, 43(1): 72-75. |
| [23] | 刘少勇, 王华忠, 张兵. 起伏地表Kirchhoff积分法叠前深度偏移方法研究与应用. 2010年国际石油地球物理技术交流会会议专刊, 科学出版社 , 2010: 49–54. Liu S Y, Wang H Z, Zhang B. Kirchhoff integral PSDM for rugged topography: technology and application (in Chinese). Conference on Exploration Geophysics in the West of China, Science Press (in Chinese) , 2010: 49-54. |
| [24] | 刘国峰, 刘洪, 李博, 等. 起伏地表直接叠前时间偏移. 石油地球物理勘探 , 2010, 45(2): 196–200. Liu G F, Liu H, Li B, et al. Direct pre-stack time migration on rugged topography. Oil Geophysical Prospecting (in Chinese) , 2010, 45(2): 196-200. |
| [25] | 张建中, 陈世军, 徐初伟. 动态网络最短路径射线追踪. 地球物理学报 , 2004, 47(5): 899–904. Zhang J Z, Chen S J, Xu C W. A method of shortest path raytracing with dynamic networks. Chinese J. Geophys. (in Chinese) , 2004, 47(5): 899-904. |
| [26] | 张东, 谢宝莲, 杨艳, 等. 一种改进的线性走时插值射线追踪算法. 地球物理学报 , 2009, 52(1): 200–205. Zhang D, Xie B L, Yang Y, et al. A ray tracing method based on improved linear traveltime interpolation. Chinese J. Geophys. (in Chinese) , 2009, 52(1): 200-205. |
| [27] | 蒋先艺. 基于二维与三维复杂结构模型正演的地震数据采集设计方法研究. 成都: 成都理工大学地球物理学院, 2004. Jiang X Y. The Study of Seismic Acquisition Method Based on Modeling of 2D and 3D Complex Geologic Structure (in Chinese). Chengdu: College of Geophysics of Chengdu University of Technology, 2004. |
| [28] | 赵瑞, 白超英. 复杂层状模型中多次波快速追踪——一种基于非规则网格的最短路径算法. 地震学报 , 2010, 32(4): 433–444. Zhao R, Bai C Y. Fast multiple ray tracing within complex layered media: the shortest path method based on irregular grid cells. Acta Seismologica Sinica (in Chinese) , 2010, 32(4): 433-444. |
| [29] | 黄国娇, 白超英. 二维复杂层状介质中地震多波走时联合反演成像. 地球物理学报 , 2010, 53(12): 2972–2981. Huang G J, Bai C Y. Simultaneous inversion with multiple traveltimes within 2-D complex layered media. Chinese J. Geophys. (in Chinese) , 2010, 53(12): 2972-2981. |
| [30] | 白超英, 黄国娇, 李忠生. 三维复杂层状介质中多震相走时联合反演成像. 地球物理学报 , 2011, 54(1): 182–192. Bai C Y, Huang G J, Li Z S. Simultaneous inversion combining multiple-phase travel times within 3D complex layered media. Chinese J. Geophys. (in Chinese) , 2011, 54(1): 182-192. |
| [31] | Sun J G, Sun Z Q, Han F X. A finite difference scheme for solving the eikonal equation including surface topography. Geophysics , 2011, 76(4): T53-T63. DOI:10.1190/1.3580634 |
| [32] | Sun J G, Sun Z Q, Han F X. A finite difference scheme for solving the eikonal equation including surface topography. 79th Ann. Internat. Mtg., Soc. Expl. Geophys., Expanded Abstracts, 2009: 2667-2671. |
| [33] | Sun J G, Sun Z Q, Han F X. A method for computing the 2D seismic traveltime including the surface topography by combining the linear interpolation and the narrow band technique. CPS/SEG 2009 International Geophysical Conference & Exposition, 2009: ID-1187. |
| [34] | 孙章庆, 孙建国, 韩复兴. 复杂地表条件下基于线性插值和窄带技术的地震波走时计算. 地球物理学报 , 2009, 52(11): 2846–2853. Sun Z Q, Sun J G, Han F X. Traveltimes computation using linear interpolation and narrow band technique under complex topographical conditions. Chinese J. Geophys. (in Chinese), (in Chinese) , 2009, 52(11): 2846-2853. |