2012, Vol. 55
2. 国土资源部应用地球物理综合解释理论开放实验室-波动理论与成像技术实验室, 长春 130026;
1. College for Geoexploration Science and Technology, Jilin University, Changchun 130026, China;
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
2. ;
1. ;
2.
地震数据处理技术效果不好是起伏地表地区地震勘探经常会面临的一个挑战，其主要原因是现有的地震数据处理技术很多都基于水平地表假设，所以研究专门针对起伏地表的三维地震数据处理技术是很有必要的[1].地震波走时计算技术是很多地震数据处理技术[2-7]中必不可少的正演模拟工具，因此对起伏地表条件下的三维地震波走时计算技术进行研究，可以为研发能很好适应起伏地表条件的三维地震数据处理技术打下基础，进而有助于提高起伏地表地区的三维地震勘探的效果.
计算地震波走时的方法主要有两大类:射线追踪和数值求解Eikonal方程.其中，射线追踪方法主要包括:两点射线追踪法[8]､波前构建法[9-10]､最短路径射线追踪法[11-12]､插值法[13-14]､辛几何算法[15-16]等;数值求解Eikonal方程的方法主要是指有限差分法[17-22].关于三维起伏地表条件下的地震波走时计算，国内外学者主要采用射线追踪类方法(包括:两点射线追踪法[23]､插值法[24]和最短路径射线追踪法[25]等)直接或间接地研究过该问题，并取得了一些研究成果.实际上，由于各种差分格式和各种波前扩展方式的提出和发展，数值求解Eikonal方程的方法已具有计算效率高､精度易调控(可以通过差分格式阶数的改变来调控算法的精度)､适用于复杂介质模型､稳定性好､易于编程实现等优点[18-22].
综合利用有限差分法所具备的优点，笔者曾基于对该方法的改进求解过二维起伏地表条件下的地震波走时计算问题[26-27].与二维空间中的算法相比，有限差分法在处理三维空间中的起伏地表问题时，必然会遇到一些新的问题，但算法的基本原理决定了其具有处理三维空间问题的能力.因此，本文提出了三维空间中算法的具体内容:① 三维空间中起伏地表模型的不等距网格建立方法;② 综合运用不等距差分格式､迎风差分格式､Huygens 原理和Fermat原理推导建立三维空间中适应于不等距网格的局部走时计算公式;③ 基于对常规窄带技术的改进来获得三维起伏地表条件下算法的整体实现步骤;④对方法进行精度分析，并通过引入网格加密技术在不过多牺牲计算效率的情况下实现计算精度的大幅度提高;⑤通过一个计算实例检验算法在处理三维起伏地表､地下复杂地质构造时的有效性和适应性.
2 三维起伏地表模型的网格剖分为了准确地描述三维起伏地表的形态，并为后续数值计算打下基础，采用不等距网格剖分三维起伏地表模型.如图 1所示，不等距网格在进行三维起伏地表模型剖分时，在地表附近的区域采用纵向上网格间距不等的网格，在除地表附近以外的空间采用规则的立方体网格.不等距网格剖分三维起伏地表模型具有如下特点:① 它能准确定位地表上采样点的高程位置，而不用作“利用就近的规则网格节点代替地表采样点"的近似处理;②与地表附近的四面体网格和其他网格形态相比，不等距网格相对简洁，这也会让地表附近的局部走时算法变得相对简洁(这也是有限差分法的一个优点).
|
图 1 剖分三维起伏地表模型的不等距网格 Fig. 1 The model discretization for the 3D undulating surface model by the grids with nonuniform grid spacing |
为了让有限差分法能够适应如图 1所示的不等距网格，首先需要建立该网格中的局部走时计算公式，当然该局部走时计算公式首先必须满足地震波的传播规律，同时还需要在面对不等距网格时具有无条件稳定性.综合分析有限差分法中各种差分格式后发现迎风差分格式具有无条件稳定性[18]，但是该差分格式不能直接应用于不等距网格中.为了解决该问题，本文通过在迎风差分格式中融入不等距差分格式､Huygens原理和Fermat原理来建立三维起伏地表条件下的局部走时计算公式.其中迎风差分格式的思想用于控制算法的稳定性，不等距差分格式用于适应不等距网格，而Huygens原理和Fermat原理则用于保证算法在地震波传播规律的控制下进行.
下面首先简要地回顾一下常规迎风差分格式的基本原理.为了便于阐述，首先引入三维空间中的Eikonal方程:
|
(1) |
这里t(x，y，z)､s(x，y，z)均是空间坐标(x，y，z)的函数，它们分别表示走时和慢度.
如图 2所示，在立方体网格中，采用迎风差分格式[18]对方程(1)进行离散可得
|
图 2 常规迎风差分格式的局部走时计算公式 Fig. 2 The local traveltime computation formula of the conventional upwind finite-difference scheme |
|
(2) |
其中，max(A，B，C)为求取A，B，C中的最大值的函数，
例如1阶x方向上(另外两个方向与x方向上的形式相同)的差分算子可表示为
|
(3) |
其中h为立方体网格的网格间距.
公式(2)中的max(l)项是根据被算点周围点的走时值的已知情况及大小来求取的:
|
(4) |
同理可以求取max(y)和max(z).因为在计算过程中，当前被计算点的所有邻近点中必然包括窄带技术(关于窄带技术的相关内容将在下一节中详细阐述)选取的扩展点(该点的走时值为已知)，所以max(x)､max(y)和max(z)必然不会同时为零(若它们同时为零，为造成(2)式的不成立，进而无法保持计算的稳定性).将3个方向上选取的差分格式代入(2)中并解该一元二次方程即可获得计算点(i，j，k)的走时值的公式:
|
(5) |
其中
公式(5)即为三维空间的规则立方体网格中的常规迎风差分格式的局部走时计算公式.该公式显然不能处理由三维起伏地表引起的地表附近的不规则网格问题.因此需要针对不规则网格推导建立新的局部计算公式，在此采用的方法是在常规迎风差分格式中，综合融入Fermat原理､Huygens原理和不等距差分格式.
|
图 3 三维起伏地表条件下的局部走时计算公式 Fig. 3 The local traveltime computation formula in 3D undulating surface condition |
如图 3所示，根据地表附近的被算点在不等距网格中所处位置的不同，可将处于不规则网格中的节点类型分为:①位于地表表面上的网格节点(如图 3中的点H);②与地表表面上的网格节点直接相邻的网格节点(该节点到地表的距离小于一个网格间距，即如图 3中的M点).下面将分别推导这两类网格节点的局部走时计算公式:
①当M点为当前被算点时;
设HM= d，点X1､X2､Y1､Y2､M､Z1､H的走时值分别为
|
(6) |
其中:tX= min(tX1，tX2)，tY= min(tY1，tY2);然后，在M点Z方向向下的方向直接采用常规迎风差分格式:
|
(7) |
最后，分别求解(6)､(7)对应的一元二次方程，并根据Fermat原理求取二者中的较小者作为最终的计算结果，有
|
(8) |
式中:θ =h2/d2.公式(8)即为计算M点的走时值时的局部走时计算公式，该公式是通过在迎风差分格式中融入不等距差分格式，并综合应用Fermat原理而获得的.
②当H点为当前被算点时.
设点H1､H2､H3､H4 的走时值分别为tH1､tH2､tH3､tH4.此时H点处的局部走时公式的建立问题将变得更加复杂，因为在该点处无法建立迎风差分格式，同时H点在X､Y两个方向上的网格节点均为地表上的点(如图 3 所示的点H1､H2､H3､H4)，并且绝大多数情况下与H点不在同一高程上.为了处理这些问题，首先，同样在M点处将不等距差分格式融入迎风差分格式中得公式(6)，不过此时的未知量变为了tH，而解此以tH为未知量的一元二次方程有:
|
(9) |
然后，在H点的X､Y两个方向上的地表点H1､H2､H3､H4 处分别采用Huygens原理，即将这四个点分别作为子震源点，则有:
|
(9) |
其中HHi表示点Hi与点H之间的距离;最后，利用Fermat原理，综合考虑(9)､(10)有:
|
(11) |
公式(11)即为计算H点的走时值时的局部走时计算公式，该公式是通过在迎风差分格式中融入不等距差分格式，并综合应用Fermat原理和Huygens原理而获得的.
综合上述推导可知，三维起伏地表条件下走时计算的局部计算公式，因为被算点所处的位置不同而不同.当计算规则网格节点的走时值时，采用公式(5);当计算地表附近的不规则网格中的节点的走时值时则采用公式(8)和(11)，其中公式(8)用于计算地表表面上的网格节点的走时值，而公式(11)则用于计算与地表表面上节点直接相邻的网格节点的走时值.
4 算法的整体实现步骤以上分别建立了三维起伏地表模型和推导了局部走时计算公式，除此之外还需要提出算法的整个实现步骤.不同于水平地表问题，这里算法的整体实现步骤需要在保证满足地震波传播规律(保证算法的合理性，即算法的实现策略不破坏地震传播过程的因果关系)的同时，还能处理起伏地表条件造成的不规则计算边界的问题.综合分析该问题，并充分结合笔者在处理二维空间中的类似问题时的研究经验[26-27]可知，作为有限差分法中的一种稳定､高效且灵活的波前扩展方式，窄带技术完全具有处理三维起伏地表问题的潜力[18-21].因此下面将以窄带技术的基本思想为基础，以设定新的网格节点类型的方式来建立三维起伏地表条件下的算法的整体实现步骤.
如图 4a所示，常规窄带技术在实现时把整个计算区域的所有网格节点分为三种类型，包括:①完成点，表示走时计算已经完成的节点，原文为“Accepted"[18]在此意译为“完成点"，即图 4a中黑色填充的圆点;②窄带点，表示近似波前上的节点，即图 4a中灰色填充的圆点;③远离点，表示走时计算还未进行的节点，即图 4a中白色填充的圆点.窄带技术的实现分为初始化和循环计算两个部分.初始化时，仅有震源点被设定为完成点;与震源点直接邻近的几个网格节点被设定为窄带点，它们构成初始窄带(即初始波前);除此之外的那些网格节点均被设定为远离点.循环计算包括如下步骤:
|
图 4 窄带技术 (a)常规窄带技术；（b)三维起伏地表条件下的窄带技术 Fig. 4 The narrow band technique (a)The conventional narrow band technique； (b)The narrow band technique in 3D undulating surface condition. |
步骤1:通过比较找出窄带内走时值最小的窄带点作为子震源点，并将其属性重新设定为完成点;
步骤2:对子震源直接邻近的6 个网格节点的类型进行判断，如果其类型为远离点则采用公式(5)计算它的走时值，并将其类型重新设定为窄带点;如果其类型为窄带点则再次计算它的走时值，并与其原来的走时值比较，较小者被保留;如果其类型为完成点则不对其作任何处理.
步骤3:判断窄带点是否还存在，是则跳回步骤1重复循环计算，否则计算结束.
基于以上阐述可知，常规窄带技术通过近似描述波前的窄带的不断推进，来模拟波前在介质中的传播过程.其中窄带技术实现过程中各个网格节点的类型的不断交替变化是窄带技术的灵魂，因为正是这样的交替变化才实现了窄带的不断推进.很显然这样的网格类型的不断交替变化在起伏地表条件下也是可以实现的，当然必须对其进行一些必要的改进.改进的目的是保证窄带的不断推进只在地表以下进行(因为地震波不能穿越地空界面进入空气中)，同时还能根据节点的位置对应地选取计算公式(因为常规的计算公式在地表附近会失效).
为了达到上述目的，在此采用的策略是在窄带技术中设定两种新类型的网格节点，即如图 4b的地表点(图 4b中起伏地表表面的点)和地表上点(地表以上的网格节点，图 4b没有画出，但计算时由于采用规则的数组，所以在计算时存在地表以上的点).当然有别于其它原有的三种类型的网格节点，这里的地表点的属性在计算过程中可能分别为完成点､窄带点或远离点(如图 4b所示地表上网格节点的填充颜色有三种可能情况，即黑色､灰色､白色).与常规窄带技术相比，起伏地表条件下的窄带技术只需要对上述步骤2进行改进:
步骤2:对子震源直接邻近的6 个网格节点的类型进行判断，如果其类型为远离点则采用公式(5)､(8)或(11)计算它的走时值(公式的选择是通过判断该点是否为地表点或为与地表点直接相邻的节点来实现的)，并将其类型重新设定为窄带点;如果其类型为窄带点则再次计算它的走时值，并与其原来的走时值比较，较小者被保留;如果其类型为完成点或地表上点则不对其作任何处理.
通过改进步骤2而获得的起伏地表条件下的窄带技术是否能适应不规则计算边界并满足波前传播的规律呢?当然能，因为:①起伏地表条件下的窄带技术保留步骤1不变，即找出波前上的最小走时点作为子震源点的方法是与Huygens原理和Fermat原理相符的;②修改后的步骤2中加入“如果其类型为完成点或地表上点则不对其作任何处理"的操作是为了保持窄带不能扩展进入地表以上区域，这与地震波不能穿越地空界面进入空气中的物理事实是相符的.
5 计算精度分析及实例为了分析上述提出的算法的计算精度，在此采用如图 5a所示的模型.模型的地表为一山峰，整个计算区域的大小为8.0km×8.0km×6.0km, 地下介质的地震波速度为1.0km/s, 计算时分别采用50.0m､20.0m 的网格间距，震源点置于(4.0km, 4.0km, 6.0km)点处，并且为了使得地表上各个点的走时值可以通过解析公式计算而获得，特设计地表上的各点可以通过直线直接与震源点相连接.
|
图 5 计算精度分析 (a)三维起伏地表模型;(b)网格加密示意图;(cl,c2,c3)网格间距分别为50. 0 m､20. 0 m及加密E 相对误差分布;(dl,d2,d3)网格间距分别为50. 0 m､20. 0 m及加密时Y=4. 0 km剖面内相对误差分布 Fig. 5 The computational accuracy analysis (a) The 3D undulating surface model; (b) An example for the process of reducing the grid spacing with 狀=5; (cl, c2, c3) The relative error distribute on the earth;s surface when the grid spacing is 50. 0 m, 20. 0 m and reduced respectively; (dl, d2, d3) The relative error distribute in the section Y=4. 0 km when the grid spacing is 50. 0 m, 20. 0 m, and reduced respectively. |
图 5c1､图 5c2分别给出了采用50.0m､20.0m为网格间距时，地表上各个网格节点的走时值的数值计算结果相对于解析值的相对误差.图 5d1､图 5d2分别给出了采用50.0 m､20.0 m 为网格间距时，剖面Y=0.0km 内的走时值的数值计算结果相对于解析值的相对误差.这些相对误差的分布均体现出一个共同的特点:震源附近区域的相对误差最大，远离震源的区域相对误差小，并且离震源越远的区域误差越小.
基于上述误差分布规律，为了在不过多牺牲计算效率的情况下大幅度提高算法的计算精度，在此引入一种震源附近的网格加密技术[28].如图 5b 所示，该网格加密技术对震源附近一定区域内的网格单元依次进行2n×2n，2n-1×2n-1，…21×21 的细分，这样的细分使得在实际计算时网格在震源附近区域采用很小的网格间距，而随着离震源距离的增加网格间距逐渐增加，直到网格间距变为常规网格间距.本文选择n=5对常规网格间距为20.0 m 的情况进行网格加密处理(关于网格加密的基本原理和实现策略请读者详见参考文献[28])，图 5c3､图 5d3分别给出了地表上和Y=0.0km 剖面上各网格节点的相对误差，与图 5c2､图 5d2 对比分析发现:算法的整体计算精度有大幅度提高，尤其是在震源附近区域.此外，网格加密后算法的整体计算量仅增加了3.71%.
以上得出的精度分析结果表明:本文提出的算法能够很好地解决三维起伏地表条件下的地震波走时计算问题，并且算法能够获得很好的计算精度.与此同时网格加密技术能在不过多牺牲计算效率的前提下大幅度提高计算精度.
为了检验本文算法能否稳定､灵活及有效地处理三维起伏地表问题，下面将给出一个三维起伏地表模型算例，该模型由Marmousi模型改造而获得，模型具有地表起伏剧烈､地下(包括近地表)地质构造复杂等特点.如图 6a所示，模型的大小为9.0km×9.0km×6.0km, 计算时将震源置于(0.0 km, 0.0km, 5.0km)处，采用20.0m 的网格间距.图 6(b—d)分别给出了X=0.0km､X-Y=0.0km､Y=0.0km剖面的等时线分布.分析图 6(b—d)的计算结果可以得出:本文算法的计算结果能与地震波在起伏地表条件下的传播规律相符，同时该计算实例也充分表明本文提出的算法在处理三维起伏地表问题时，能够保证很好的稳定性､灵活性和有效性，并且算法能同时处理三维剧烈地表起伏､地下(包括近地表)复杂地质构造等问题.
|
图 6 三维起伏地表模型中的地震波走时分布 (a)三维起伏地表模型;(b—d)X=0.0 km､X—Y=0.0 km､Y=0.0 km剖面内的地震波走时分布. Fig. 6 The traveltime contours are superimposed on the 3D undulating surface model (a)The 3D undulating surface model； (b—d)The traveltime contours are superimposed on the section X=0. 0 km, Y=0. 0 km, and Y=0. 0 km respectively. |
三维起伏地表条件下的地震波走时计算技术是研究三维起伏地表地区很多地震数据处理技术的基础性工具.为了获得适应于三维起伏地表且精度高的走时算法，本文提出了一种不等距迎风差分法.算法通过在常规迎风差分格式中综合引入不等距差分格式､Huygens原理和Fermat原理来获取地表附近的局部走时计算公式，同时通过改进常规窄带技术来获取算法的整体实现步骤.针对算法的精度分析及算例分析表明本文提出的算法具有如下特点:①由于采用不等距网格剖分起伏地表模型，并采用有限差分法作为局部走时算法，所以方法简洁且易于编程实现;②方法的实现步骤能够满足地震波传播的基本规律，即Huygens原理､Fermat原理以及地震波在地空边界处的传播特点;③ 算法能够保证很高的计算精度，同时网格加密技术能在不过多牺牲计算效率的情况下大幅度提高计算精度;④算法能稳定､灵活及有效地处理剧烈起伏的三维地表､地下(包括近地表)复杂地质构造等问题.
| [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) (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(5): 1045-1076. |
| [5] | Sun J G. On the limited aperture migration in two dimensions. Geophysics , 1998, 63(5): 984-994. |
| [6] | 井西利, 杨长春, 王世清. 一种改进的地震反射层析成像方法. 地球物理学报 , 2007, 50(6): 1831–1836. Jing X L, Yang C C, Wang S S. A improved seismic reflection tomographic method. Geophys. (in Chinese) (in Chinese) , 2007, 50(6): 1831-1836. |
| [7] | 瞿辰, 周蕙兰, 赵大鹏. 使用纵波和横波走时层析成像研究菲律宾海板块西边缘带和南海地区的深部结构. 地球物理学报 , 2007, 50(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. Geophys. (in Chinese) (in Chinese) , 2007, 50(6): 1757-1768. |
| [8] | Julian B R, Gubbins D. Three-dimensional seismic ray tracing. .Geophys. , 1977, 43: 95-113. |
| [9] | Vinje V, Iversen E, Gj?ystdal H. Traveltime and amplitude estimation using wavefront construction. Geophysics , 1993, 58(8): 1157-1166. DOI:10.1190/1.1443499 |
| [10] | Han F X, Sun J G, Sun Z Q, et al. Positioning of grid points in wavefront construction. Applied Geophysics , 2009, 6(3): 248-258. DOI:10.1007/s11770-009-0024-z |
| [11] | Moser T J. Shortest path calculation of seismic rays. Geophysics , 1991, 56(1): 59-67. DOI:10.1190/1.1442958 |
| [12] | 张美根, 贾豫葛, 王妙月, 等. 界面二次源波前扩展法全局最小走时射线追踪技术. 地球物理学报 , 2006, 49(4): 1169–1175. Zhang M G, Jia Y G, Wang M Y, et al. A global minimum traveltime ray tracing algorithm of wavefront expanding with interface points as secondary sources. Geophys. (in Chinese) (in Chinese) , 2006, 49(4): 1169-1175. |
| [13] | Asakawa E, Wanaka T. Seismic ray tracing using linear traveltime interpolation. Geophysical Prospecting , 1993, 41(1): 99-111. DOI:10.1111/gpr.1993.41.issue-1 |
| [14] | Wang H Z, Ma Z T. Traveltime calculation in 3D media with arbitray velocity distribution. 69th Ann. Internat. Mtg., Soc. Expl. Geophys., Expanded Abstracts, 1999: 1778-1781. Traveltime calculation in 3D media with arbitrary velocity distribution |
| [15] | 秦孟兆, 陈景波. Maslov渐近理论与辛几何算法. 地球物理学报 , 2000, 43(4): 523–533. Qin M Z, Chen J B. Maslov asymptotic theory and symplectic algorithm. Geophys. (in Chinese) (in Chinese) , 2000, 43(4): 523-533. |
| [16] | 高亮, 李幼铭, 陈旭荣, 等. 地震射线辛几何算法初探. 地球物理学报 , 2000, 43(3): 402–409. Gao L, Li Y M, Chen X R, et al. An attempt to seismic ray tracing with symplectic algorithm. Geophys. (in Chinese) (in Chinese) , 2000, 43(3): 402-409. |
| [17] | Vidale J E. Finite-difference calculation of travel times. Soc. Am. , 1988, 78(6): 2062-2076. |
| [18] | Sethian J A, Popovici A M. 3-D traveltime computation using the fast marching method. Geophysics , 1999, 64(2): 516-523. DOI:10.1190/1.1444558 |
| [19] | Sethian J A. Fast marching methods. SIAM Review , 1999, 41(2): 199-235. DOI:10.1137/S0036144598347059 |
| [20] | Rawlinson N, Sambridge M. Wave front evolution in strongly heterogeneous layered mediausing the fast marching method. Geophys , 2004, 156(3): 631-647. DOI:10.1111/gji.2004.156.issue-3 |
| [21] | Rawlinson N, Sambridge M. Multiple reflection and transmission phases in complex layered media using a multistage fast marching method. Geophysics , 2004, 69(5): 1338-1350. DOI:10.1190/1.1801950 |
| [22] | 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. http://cn.bing.com/academic/profile?id=2047541605&encoded=0&v=paper_preview&mkt=zh-cn |
| [23] | 蒋先艺. 基于二维与三维复杂结构模型正演的地震数据采集设计方法研究. 成都: 成都理工大学地球物理学院, 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. http://www.oalib.com/references/19481174 |
| [24] | 刘少勇, 王华忠, 张兵. 起伏地表Kirchhoff积分法叠前深度偏移方法研究与应用. 北京: 科学出版社, 2010 : 49 -54. Liu S Y, Wang H Z, Zhang B. Kirchhoff integral PSDM for rugged topography: technology and application (in Chinese). Beijing: Science Press, 2010 : 49 -54. |
| [25] | 白超英, 黄国娇, 李忠生. 三维复杂层状介质中多震相走时联合反演成像. 地球物理学报 , 2011, 54(1): 182–192. Bai C Y, Huang G J, Li Z S. Simultaneous inversion combining multiple-phase traveltimes within 3D complex layered media. Geophys. (in Chinese), (in Chinese) , 2011, 54(1): 182-192. |
| [26] | 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., Expl. Geophys., Expanded Abstracts, 2009: 2667-2671. http://cn.bing.com/academic/profile?id=2320201466&encoded=0&v=paper_preview&mkt=zh-cn |
| [27] | 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 |
| [28] | Sethian J A, Popovici A M. 3-D traveltime computation using the fast marching method. Geophysics , 1999, 64(2): 516-523. DOI:10.1190/1.1444558 |