直达波旅行时线性插值算法

直达波旅行时线性插值算法
黄翼坚;朱光明;敦敏
【摘要】文中首先简要介绍了初至波旅行时LTI算法,并分析了该法的计算结果包含非直达波旅行时的原因,由此提出通过对射线方向进行控制实现直达波旅行时计算的方法.直达波旅行时计算步骤与初至波旅行时LTI算法的计算步骤基本一样,只有两个不同之处:其一,旅行时赋初值不同,对震源点以外的节点,初至波旅行时LTI算法赋予大干该节点初至波旅行时的值,而直达波旅行时LTI算法则赋予大干该节点直达波旅行时的值;其二,在每次更新节点旅行时的过程中,一旦发现有波从某个低速地层传向某个高速地层,则在以后的计算中,判断射线是否是从该高速层指向该低速层,如果是则舍去该射线的旅行时.平层速度模型和具有起伏速度界面的多层模型算例说明,本文提出的直达波旅行时LTI迭代算法简单而有效,反射波射线追踪算例说明基于该算法可以进行大角度反射波射线追踪.
【期刊名称】《石油地球物理勘探》
【年(卷),期】2010(045)002
【总页数】5页(P225-229)
【关键词】初至波旅行时;直达波旅行时;线性插值;算法;费马原理;数值模拟
【作者】黄翼坚;朱光明;敦敏
【作者单位】长安大学地质工程与测绘学院,陕西西安,710054;长安大学地质工程与测绘学院,陕西西安,710054;东方地球物理公司物探技术研究中心,河北涿
州,072751
【正文语种】中文
复杂速度模型中地震波旅行时的计算是基尔霍夫积分法叠前深度偏移及层析成像的一个关键步骤。

自20世纪80年代末以来,国内外许多学者研究了初至波旅行时的计算方法。

例如,Vidale[1]提出了初至波旅行时程函方程有限差分解法;Moser[2]提出了最短路径射线追踪及初至波旅行时计算方法;A sawaka和 Kawanaka[3]提出了初至波旅行时线性插值(L TI)算法;刘洪等 [4]提出了界面网全局旅行时计算技术;王华忠等[5]提出了动态规划法旅行时计算方法;赵连锋等[6]提出了有序波前重建法射线追踪技术;秦义龙等[7]利用单频双程波动方程计算初至波旅行时。

但是,在复杂地质构造情况下(例如M armousi模型),由于初至波只携带了较小的地震波能量,因此基于初至波旅行时的偏移方法得不到较精确的地震成像[8]。

为了提高偏移效果, Nichols[9]提出了任意复杂二维介质中地震波主能量旅行时算法,匡斌等 [10]把该算法扩展到三维介质的情况;Shin等[11]提出了利用频率域单程波动方程计算直达波旅行时的方法;Le-Wei Mo和 Harris[12]提出了直达波旅行时程函方程有限差分算法。

本文基于 A sawaka等[3]、Faria等 [13]的工作,提出直达波旅行时线性插值(L TI)迭代算法。

文中首先简要介绍了初至波旅行时L TI算法,并分析了该法的计算结果包含非直达波旅行时的原因,由此提出通过对射线方向进行控制实现直达波旅行时计算的方法。

最后利用两个旅行时算例以及一个反射波射线追踪算例说明了算法的正确性和应用价值。

假设速度模型由多个正方形网格参数化,各网格的速度值为常数。

对其中任一网格节点 C,其四周有8个与其相邻的节点(图1),通过这8个节点中两两相邻的节点对到达节点 C的射线有8种可能模式[13](图1中带箭头虚线所示)。

当节点 C位于模型边界附近时,可能的射线模式数相应减小。

根据文献[3]和文献[13]可以归纳初至波旅行时的L TI算法。

首先,把震源所在节点的旅行时赋予零值,其他网格节点的
旅行时赋予一个大于该节点初至波旅行时的值(如,以该节点到震源距离的若干倍除以平均速度得到的值,或由直射线假设估算的旅行时值)。

然后,对所有网格节点逐个做以下运算：对每个网格节点,依次按照图1所示的8种可能的射线模式计算该节点(对应图中的 C点)的旅行时,并统计这些旅行时的最小者。

若所得的最小值小于该节点旅行时的原有旅行时,则用该最小旅行时更新该节点的旅行时,并记录这个旅行时对应的波在节点C的传入方向。

下面以图1中的射线模式1为例说明节点C 旅行时的算法。

假设A点的旅行时为tA,B点的旅行时为tB,计算点 C的旅行时为tC,穿过线段到达C点的射线与线段交于P点,P点的旅行时记为tP。

假设tP可由tA和tB线性插值得到,即
其中V为A、B和C三点所在网格的速度。

根据费马原理,波沿射线传播使旅行时取得最小值,因此由式(2)取tC对y的偏导数,并令该偏导数为零,得到关于y的一元二次方程
根据y的取值范围[0,H]可知:①当式(3)有一个在[0,H]上的实根 y1时,取y=y1;②当式(3)有两个在[0,H]上的实根 y1和 y2时,y取使得式(2)中的旅行时 tC较小的 y1或 y2值;③当式(3)没有在[0,H]上的实根时,y取使得式(2)中的旅行时tC较小的0或 H值(y=0表示 P点与B点重合,y=H表示P点与A点重合)。

由式(3)可以求出使tC取得最小值的y值,从而可确定 P点的位置,进而可确定模式1的波在C点的入射角度。

以上过程为一次迭代过程,可通过多次迭代取得更精确的结果,终止条件是前、后两次迭代各个网格节点的最大旅行时更新量小于某个指定的较小正数。

直达波是指传播过程中没有经过折射、反射或者绕射的波[12]。

下面首先分析初至波的几种可能的波型,并分析初至波旅行时L TI算法所得结果包含非直达波旅行时的原因,然后对初至波旅行时L TI算法进行修改,从而得到直达波旅行时L TI算法。

如图2所示,根据惠更斯原理,当震源点S位于上面的低速层时,上层中的任一点除了
可以接收到从震源直接传来的直达波(如图2中射线对应的波)以外,还可以接收到以速度界面1上的点为新震源返传回来的广义绕射波,包括折射波(如图2中的折线对应的波)以及从中层传到上层的透射波或其他类型的广义绕射波(如图2中折线所对应的波,注意,该波在N点穿越界面1时不一定满足斯内尔定律),这些波为非直达波。

虽然上述非直达波的传播路径比直达波长,但由于它们经过下面的高速地层,因此其
传播时间有可能比直达波小,从而可能成为初至波。

如果按照前文介绍的初至波L
TI方法计算,由于该算法允许所有波型在所有网格节点向各个方向传播,因此得到的旅行时就会包括这些非直达波的旅行时。

根据以上分析可知,有可能成为初至波的非直达波有两个特点：①这些波从某个高
速地层(记为H层)向某个低速地层(记为L层)传播;②在这些波形成之前,有波从低
速地层(L层)传入高速地层(H层)。

如果在有直达波从某个低速地层进入某个高速
地层以后,禁止所有从高速层返向低速层传播的波传播,而允许其他波传播,就可以在初至波L TI算法结果中去掉这些非直达波的旅行时。

例如,图2中的中间层为高速
地层,该层有从低速地层(上层)传来的直达波(射线对应的波),如果此后禁止波从该
高速层(中层)返向该低速层(上层)传播,就可以避免计算得到两条射线对应的非直达波的旅行时。

而对射线对应的波,虽然它也是从高速层(中层)向低速层(下层)传播,但是由于传入高速层(中层)的直达波不是从低速层(下层)传来,而是从另一低速层(上层)传来,因此这个波的传播不受限制,由图2可见对应的波在传播过程中没有经过折射、反射或者绕射(为直达波),其旅行时是我们想要计算的。

综上所述,可得到直达波旅行时L TI算法,其计算步骤与初至波旅行时L TI算法的计算步骤基本一样,只有两个地方不同。

其中之一是旅行时赋初值的不同,对震源点以
外的节点,初至波旅行时L TI算法赋予大于该节点初至波旅行时的值,而直达波旅行时L TI算法赋予大于该节点直达波旅行时的值;另一个不同点是,在每次更新节点旅行时的过程中,一旦发现有波从某个低速地层传向某个高速地层,则在以后的计算中,
判断射线是否是从该高速层指向该低速层,如果是则舍去该射线的旅行时。

文中用两个旅行时数值模拟算例说明文中所提方法的正确性,然后再给出该法应用于大角度反射波射线追踪的例子。

4.1 双平层模型算例
图3为双平层模型旅行时算例。

由波场快照(t=0.6s)(图3a)可见：直达波和透射波的波前有间断,而折射波的波前把这两者的波前连接起来,形成连续的初至波波前;折射波携带的能量明显比与震源同层的直达波、透射直达波以及反射波弱。

由初至波旅行时等值线图(图3b)和直达波旅行时等值线图(图3c)可见：前者的时间场等值线与波动方程正演结果的初至波波前相一致,包含有折射波部分,后者的时间场等值线与波动方程正演结果的直达波波前相一致,没有折射波部分。

这说明后者是正确的。

4.2 多层模型算例
图4为多层模型旅行时算例,模型由6套地层组成。

由图4可见：①在初至波旅行时等值线图(图4a)中,根据等值线的梯度方向可知,在黄色虚线所圈定的区域内,初至波旅行时时间场包含有从速度界面向上传播的非直达波所对应的旅行时;②在直达波旅行时等值线图(图4b)中,在相同的区域内没有非直达波所对应的旅行时;③沿等值线的负梯度方向,在模型右边界、以20m为等间距的接收点所对应的透射波射线路径中(图4a和图4b中的粉红色曲线),由射线路径可知,按初至波旅行时得到的透射波射线包含有非直达波部分(图4a中白色箭头所指射线),而按直达波旅行时得到的是直达波射线。

图4说明直达波旅行时算法是有效的。

4.3 反射波射线追踪算例
文献[12]用实例说明,在地震波数据偏移方面,利用直达波旅行时较初至波旅行时可获得更好的偏移效果,本文给出了在应用于大角度反射波射线追踪方面,直达波旅行时相对初至波旅行时具有优势的例子。

基于旅行时时间场的反射波射线追踪方法,
不需要对速度模型做光滑处理,可适应于模型速度变化大的情况,计算速度也较快[14,15]。

其做法是：首先计算炮点旅行时时间场,再把震源放于接收点计算接收点旅行时时间场;接着沿速度界面搜索一个使炮点旅行时与接收点旅行时之和最小的点作为反射点,然后根据炮点和接收点时间场,从反射点分别向炮点和接收点做透射波射线追踪即可。

图5为反射波射线追踪算例,由图5可见：在基于直达波旅行时等值线图(图5b)中,红、绿和紫三种颜色曲线分别对应从上到下三个速度界面的反射波射线路径;在基于初至波旅行时等值线图(图5a)中,除了紫色曲线完全对应第三个界面的反射波射线外,红色和绿色曲线中都有一部分对应先在下边的较高速度层传播、再传回上边较低速度层的非直达波(包括折射波)射线(图5a中红色射线以及绿色射线),这是由于在入射角度较大的情况下,图5a的射线追踪所用的旅行时包含有非直达波旅行时所致。

图5说明,当入射角度较大时,基于直达波旅行时可以进行大角度反射波射线追踪,基于初至波旅行时所追踪到的可能是折射波射线。

本文在对非直达波的传播路径进行分析的基础上,通过对传统初至波L TI算法的射线方向进行约束,得到了直达波旅行时L TI迭代算法。

平层速度模型和具有起伏速度界面的多层模型算例说明,本文提出的直达波旅行时L TI迭代算法简单而有效,反射波射线追踪算例说明基于该算法可以进行大角度反射波射线追踪。

文中算法可用于地震波数据偏移、射线追踪以及旅行时层析成像。

[1]Vidale J E.Finite-difference calculation of travel times. Bull Seism Soc Amer,1988,78(6):2062～2076
[2]Moser T J.Shortest path calculation of travel time.
Geophysics,1991,56(1):59～67
[3]Asawaka E,Kawanaka T.Seismic ray tracing using linear traveltime interpolation.Geophysical Prospec-ting,1993,41(1):99～112
[4]刘洪,孟凡林,李幼铭.计算最小走时和射线路径的界面网全局方法.地球物理学
报,1995,38(6)：823～831
[5]王华忠,方正茂,匡斌等.任意介质中的动态规划法地震波三维旅行时计算.地球物
理学报,2001,44(增刊)：179～189
[6]赵连锋,朱介寿,曹俊兴等.有序波前重建法的射线追踪.地球物理学报,2003,46(3)：415～420
[7]秦义龙,张中杰,Shin Changsoo等.利用单频双程波动方程计算初至走时及其振幅.地球物理学报,2005, 48(2)：423～428
[8]Gray SH,May W P.Kirchhoff migration using eikonal equation traveltimes.Geophysics,1994,59:810～817
[9]Nichols D E.Maximum energy traveltimes calculated in the seismic frequency band.Geophysics,1996,61: 253～263
[10]匡斌,王华忠,季玉新等.任意复杂介质中主能量法地震波走时计算.地球物理学报,2005,48(2)：394～398
[11]Shin S,Shin C,Kim W et al.Traveltime calculation using mono-chromatic one-way wave equation.SEG Technical Expanded Abstracts,2000,19:2313～2316
[12]Le-Wei Mo,Harris JM.Finite-difference calculation of direct-arrival traveltimes using the eikonal equation.Geophysics,2002,67(4):1270～1274 [13]Faria E L,Stoffa P L.Traveltime computation in transversely isotropic media.Geophysics,1994,59(2): 272～281
[14]Matsuoka T,Ezaka T.Ray tracing using recip rocity.
Geophysics,1992,57(2):326～333
[15]赵改善,郝守玲,杨尔皓等.基于旅行时线性插值的地震射线追踪算法.石油物
探,1998,37(2)：14～24
A t first in this paper the direct w ave travel time L TI algo rithm was briefly introduced,then the reason for the algorithm’s calculation results to involve non-direct wave travel time was analyzed,so a new algo rithm to calculate the direct wave travel time was realized by controlling ray’s direction.The calculation steps for the direct wave travel time are almost the same as the steps for the first break travel time L TI algorithm excep t the two different aspects as below.The first aspect is the different initial value for the travel time,fo r the nodes o ther than the source point,a value w hich is larger than the value of the first break travel time at the node w as assigned in first break travel time L TIalgorithm,w hile a value w hich is larger than the value of the direct weave travel time at the node was assigned in direct wave travel time L TIalgorithm,the other different aspect exists in the p rocess of refreshing the travel time fo r the nodes,that is once wave ray was found p ropagating from low velocity layer to high velocity layer,then the travel time fo r the w ave ray would be ignored in the next calculations.Calculation examp le for horizontal layer velocity model and multilayer model w ith rugged velocity subsurface show that the direct wave travel time L TIiterative algorithm is samp le but effective,reflection ray tracing calculation examp le demonstrated that large angle reflection w ave ray tracing can be conducted by using the algo rithm.
【相关文献】
[1]Vidale J E.Finite-difference calculation of travel times. Bull Seism Soc
Amer,1988,78(6):2062～2076
[2]Moser T J.Shortest path calculation of travel time. Geophysics,1991,56(1):59～67
[3]Asawaka E,Kawanaka T.Seismic ray tracing using linear traveltime
interpolation.Geophysical Prospec-ting,1993,41(1):99～112
[4]刘洪,孟凡林,李幼铭.计算最小走时和射线路径的界面网全局方法.地球物理学报,1995,38(6)：823～831
[5]王华忠,方正茂,匡斌等.任意介质中的动态规划法地震波三维旅行时计算.地球物理学
报,2001,44(增刊)：179～189
[6]赵连锋,朱介寿,曹俊兴等.有序波前重建法的射线追踪.地球物理学报,2003,46(3)：415～420
[7]秦义龙,张中杰,Shin Changsoo等.利用单频双程波动方程计算初至走时及其振幅.地球物理学报,2005, 48(2)：423～428
[8]Gray SH,May W P.Kirchhoff migration using eikonal equation
traveltimes.Geophysics,1994,59:810～817
[9]Nichols D E.Maximum energy traveltimes calculated in the seismic frequency
band.Geophysics,1996,61: 253～263
[10]匡斌,王华忠,季玉新等.任意复杂介质中主能量法地震波走时计算.地球物理学报,2005,48(2)：394～398
[11]Shin S,Shin C,Kim W et al.Traveltime calculation using mono-chromatic one-way wave equation.SEG Technical Expanded Abstracts,2000,19:2313～2316
[12]Le-Wei Mo,Harris JM.Finite-difference calculation of direct-arrival traveltimes using the eikonal equation.Geophysics,2002,67(4):1270～1274
[13]Faria E L,Stoffa P L.Traveltime computation in transversely isotropic
media.Geophysics,1994,59(2): 272～281
[14]Matsuoka T,Ezaka T.Ray tracing using recip rocity. Geophysics,1992,57(2):326～333
[15]赵改善,郝守玲,杨尔皓等.基于旅行时线性插值的地震射线追踪算法.石油物探,1998,37(2)：14～24
A t first in this paper the direct w ave travel time L TI algo rithm was briefly introduced,then the reason for the algorithm’s calculation results to involve non-direct wave travel time was analyzed,so a new algo rithm to calculate the direct wave travel time was realized by controlling ray’s direction.The calculation steps for the direct wave travel time are almost the same as the steps for the first break travel time L TI algorithm excep t the two different aspects as below.The first aspect is the different initial value for the travel time,fo r the nodes o ther than the source point,a value w hich is larger than the value of the first break travel time at the node w as assigned in first break travel time L TIalgorithm,w hile a value w hich is larger than the value of the direct weave travel time at the node was assigned in direct wave travel time L TIalgorithm,the other different aspect exists in the p rocess of refreshing the travel time fo r the nodes,that is once wave ray was found p ropagating from low velocity layer to high velocity layer,then the travel time fo r
the w ave ray would be ignored in the next calculations.Calculation examp le for horizontal layer velocity model and multilayer model w ith rugged velocity subsurface show that the direct wave travel time L TIiterative algorithm is samp le but
effective,reflection ray tracing calculation examp le demonstrated that large angle reflection w ave ray tracing can be conducted by using the algo rithm.
Direct wave travel time linear interpolation algorithm.
Huang Yi-jian1,Zhu Guang-m ing1 and Guo M in2. 1.College of Geology Engineering and Geomatics, Chang’an University,Xi’an City,Shaanxi Province,710054,Chi na2.Geophysical Technology Research Centre,BGP Inc.of CNPC,Zhuozhou City,Hebei Province, 072751,China。