石秀林, 孙建国, 孙辉, 刘明忱, 刘志强

吉林大学地球探测科学与技术学院, 长春　130026



基于delta波包叠加的时间域深度偏移

石秀林, 孙建国*, 孙辉, 刘明忱, 刘志强

吉林大学地球探测科学与技术学院, 长春130026

摘要delta波包可由高斯波束经傅里叶逆变换得到，是高斯波束在时空域的表达.它最早出现在合成理论地震图的研究中，本文将其应用于偏移领域.通过delta波包叠加表达时间域格林函数，可将高斯波束偏移由频率域转换到时间域，再结合Rayleigh积分和激励时间成像条件，本文给出了基于delta波包叠加的深度偏移算法.该偏移算法可在时间域直接计算，但因包含褶积运算，成像时将耗费大量的计算时间.针对这一问题，本文提出了将褶积简化为乘积的近似公式.近似后的偏移算法，不仅保留了高斯波束偏移的优点，而且计算效率得到显著提升.文中通过两个数值算例验证了上述结论.

关键词delta波包； 高斯波束； 褶积运算； 数值近似； 时间域深度偏移

1引言

20世纪90年代，Hill(1990)将高斯波束应用于偏移成像，提出了基于倾斜叠加的高斯波束偏移方法.该偏移方法通过τ -p变换将波场分解为沿不同方向传播的平面波，经高斯波束传播后，在目标区域叠加成像.Hill的高斯波束偏移算法不仅计算效率高，并且具有可处理多走时，计算焦散区振幅等优点.近些年来，许多勘探地震学者对其进行了研究和改进(Hill,2001; Nowack et al.,2003; Albertin et al.,2004; Gray,2005; Gray and Bleistein,2009; 岳玉波,2011).但我们同时注意到，若利用高斯波束表达点源波场的格林函数，则无需τ -p变换，就可进行偏移计算.基于此种思路，Popov等(2010)提出了另一种严格遵守Kirchhoff型积分理论的高斯波束偏移算法.它不用进行Hill(2001)和Gray(2005)等人提出的最陡下降近似，理论上要优于Hill的算法，但缺点是耗费更多的计算时间(孙建国，2012).

基于delta波包的偏移算法可在时间域直接计算，并保留了高斯波束偏移的诸多优点，如可计算多走时、计算焦散区振幅等，但因包含褶积运算，成像时将耗费大量的计算时间.针对这一问题，本文根据褶积算子性质，通过泰勒函数展开，将褶积运算简化为乘积运算.理论上，虚部越小，误差越小，近似公式越准确.相反，虚部越大，误差越大.但当虚部变大时，由于delta波包的振幅衰减明显，误差对于叠加结果的影响变小.所以在叠加充分区域，可近似认为成像结果是准确的.在叠加不充分区域，近似误差可能影响成像结果.近似偏移算法，几乎不需褶积运算(仅需对地震记录的时间导数进行一次Hilbert变换)，计算时间仅为近似前的0.2～0.25倍，计算效率得到了显著提升.本文最后通过Marmousi和Sigsbee2a两个数值算例验证了新偏移算法的成像效果.

2基本原理

2.1延拓波场

假设S代表地表，xs=(xs,ys,zs)为震源，xg=(xg,yg,zg)为地表检波器的位置，U(xg,t;xs)为检波器xg接受到由震源xs引起的地震反射波场.x′=(x′,y′,z′)为地下成像点.根据第一种Rayleigh积分，成像点x处的延拓波场U(x,x′;t)满足表达式(BornandWolf，1999；孙建国，2012)

(1)

其中，G(x,x′;t)为检波器xg到成像点x′的时间域格林函数.∂/∂n代表沿地表外法线方向求导.角标*表示复数共轭.

在高斯波束求和方法中，需要先计算频率域格林函数，再经傅里叶逆变换得时间域格林函数.其中频率域格林函数GGB(x,x′;ω)满足表达式

(2)

其中，γ=(γ1,γ2)为射线坐标，通常可以取球坐标系或笛卡尔坐标系，区域D为射线坐标的积分范围.UGB(ω)为单条高斯波束的波场表达式.Φ(γ)为权函数，满足Φ(γ)=ΦR(γ)+isgn(ω)ΦI(γ).

当地表S水平，利用高频渐近方法，忽略掉高阶项，导数∂GGB(x,x′;ω)/∂ n可近似为(Popovetal.，2010)

(3)

其中，vg为检波器xg处的地表速度，θg为检波器xg处射线出射方向与z坐标轴正方向的夹角.

将公式(3)代入(1)，可得

(4)

上式为高斯波束方法计算延拓波场U(x,x′;t)的基本表达式.

2.2delta波包

与高斯波束求和方法不同，本文在计算时间域格林函数GGB(x,x′;t)时，采用先对单条高斯波束进行傅里叶逆变换，再在时间域叠加求和的顺序.

(5a)

(5b)

其中，gDP(t)为delta波包.

用delta波包叠加表示的时间域格林函数GDP(x,x′;t)为

(6)

其中，g(γ,t)满足表达式

g(γ,t)=Φ(γ)*gDP(γ,t)

(7)

公式(7)表明利用高斯波束理论中的权函数Φ(γ)、振幅AGB和走时τGB，可计算时间域g(γ,t).为了以后推导方便，去掉角标，使用振幅A和走时τ表示.这里振幅A和走时τ为复数，满足运动学和动力学射线追踪方程.

(8)

(8)式为delta波包方法计算延拓波场U(x,x′;t)的公式，其中g*(γ,t)满足关系式

(9)

2.3成像公式

使用激励时间成像条件，当t=0时，成像条件为I(x′,xs)=U(x,x′;t=0)

(10)

其中，g*(γ,t)满足关系式(9)，并且

(11)

这里，Φs、As和τs分别代表由震源xs到成像点x′的权函数、振幅和走时；Φg、Ag和τg分别代表由检波点xg到成像点x′的权函数、振幅和走时.

3波包解析式

通过之前的推导可知，利用公式(9)—(11)可进行时间域深度偏移.但从数值计算角度考虑，仍存在两个问题：一是delta波包表达式(9)仍需积分计算；二是成像公式(10)中褶积过程会耗费大量计算时间.

(12)其中

(13)

这里，函数b(t)为高斯波束中包含复走时τ=τR+iτI项的傅里叶逆变换.需要注意的是，当权函数Φ或振幅A与频率ω有关时，若权函数Φ或振幅A的时间域解析式存在，仍可以写成类似公式(13)的形式，但这并不属于本文的研究范围，本文将只讨论权函数Φ或振幅A与频率ω无关的情况.

本文将b(t)称为褶积算子，褶积算子的性质取决于复走时τ.由于复走时实部τR和虚部τI对算子的影响并不相同，为了便于研究，将b(t)进一步分解为两个子算子，满足b(t)=b1(t)*b2(t)，其中

(14a)

(14b)

(15)

图1给出了不同τI下褶积算子b2(t)随时间t的变化曲线.结合公式(15)可知，b2(t)为偶函数，峰值为1/πτI.τI越小，峰值越大，能量也越集中于零点附近；τI越大，峰值越小，能量越发散.

图1　不同τI时褶积算子b2(t)Fig.1　The amplitude of operator b2(t) with different τI

将g*(γ,t)的解析式(12)代入延拓波场U(x,x′;t)表达式(8)，可得

(16)

其中

f(t)=Re(Φ*A*)x(t)-Im(Φ*A*)h(t),

(17)

这里

(18a)

(18b)

h(t)为x(t)的Hilbert变换.公式(16)—(18b)可在时间域直接计算延拓波场U(x,x′;t)，无需傅里叶变换.

4数值近似

虽然利用公式(16)—(18b)可在时间域直接进行深度偏移，但公式中存在的褶积运算f(t)*b(t)会耗费大量的计算时间.从提升计算效率的角度考虑，本节将根据算子b2(t)的性质，对褶积运算进行数值近似.

4.1数值积分近似

将子算子(14a)和(14b)代入公式(16)得

(19)

设K(t)=f(t+τR)*b2(t)，将其改写成关于变量ξ的积分形式为

(20)

f(ξ+τR)=f(t+τR)+f′(t+τR)(ξ-t)

(21)

其中，η介于ξ与t之间.所以K(t)可近似为

(22)仔细观察上式，K(t)被写成三项相加求和的形式.因为(ξ-t)在ξ=t处为奇函数，b2(t-ξ)在ξ=t处为偶函数，当τI≠0时，在积分范围[t-Δt,t+Δt]内，第二项的积分值为零.则K(t)可近似为第一项，当τI≠0时，得

(23)

截断误差由第三项决定，有

(24)

4.2参数取值

通过上节的推导，K(t)被近似为

(25)

其中

(26)

这里，振幅衰减程度由衰减项A衰决定.A衰为Δt和τI的比值，并且A衰∈[0,1].设β=Δt/τI，β越小，A衰越小，振幅衰减越明显.相反，β越大，A衰趋近于1.

理论上，积分区间[t-Δt,t+Δt]越小，f(t+τR)的一阶泰勒展开越准确，所以本文建议Δt只取几个时间采样间隔.此时，当τI较小时，因为b2(t)能量聚焦于零点附近，所以近似公式(25)仍保有较高精度；当τI较大时，b2(t)的能量发散，近似公式误差变大，但由于delta波包振幅快速衰减，其对叠加结果的影响却变小.综上所述，当Δt只取几个时间采样间隔时，既满足数值近似的精度要求，也满足偏移的需要.

关于Δt的取值，本文还建议最好能保证衰减项A衰与实际振幅衰减程度相近.这里以Ricker子波为例(如图2)，褶积计算的K(t)为理论值(见图2a)，近似公式(25)计算的K(t)为近似值(见图2b)，图2c给出了K(t)峰值振幅随τI变化曲线.由图2c可以看出，随着τI的增大，峰值振幅逐渐变小.Δt越大，振幅衰减越慢；Δt越小，振幅衰减越迅速.当Δt=0.005 s时，K(t)幅值衰减程度与理论值相近.此时，K(t)近似值(图2b)与褶积值(图2a)的形态也是一致的.

图2　(a) K(t)的理论值(不同τI); (b) K(t)的近似值(不同τI)，Δt取0.005 s; (c) K(t)峰值随τI变化曲线(不同Δt)Fig.2　(a) The theoretical values of K(t) with different τI; (b) The approximate values of K(t) with different τI, Δt=0.005 s; (c) The peak amplitudes of K(t) varies with τI under different Δt

4.3近似成像公式

将公式(25)、(26)和(17)代入公式(19)，延拓波场U(x,x′;t)满足表达式

(27)

所以，成像条件I(x′,xs)可重新整理为

(28)其中

(29)这里，x(t)和h(t)分别满足公式(18a)和(18b).利用成像条件(28)和(29)可直接在时间域进行深度偏移.近似后的成像条件几乎不需褶积运算(仅需对x(t)进行Hilbert变换得到h(t)，单道地震记录只需计算一次)，计算效率得到很大的提升.

5数值算例

本节将展示两个数值算例.第一个为Marmousi模型，图3为速度模型，横向网格节点数为737，网格间距12.5m；纵向网格节点数为750，网格间隔4m.使用的地震记录共240炮，第一炮位置在x=3000m处，每炮207道，道间距和炮间距均为25m，中间放炮，两边接收，炮点位置与第104道接收位置相同.每道时间采样点为750个，采样间隔4ms.

图3　Marmousi速度模型Fig.3　Stratigraphic interval velocity model for Marmousi dataset

图4为单炮的时间域偏移结果(Δt=0.005s)，图4a、4c为包含褶积运算的偏移结果，图4b、4d为使用近似公式的偏移结果.通过对比可以看出，两种算法的偏移结果基本一致，但在细节处，前者的信噪比要高于后者，图像更为清晰些.

图4　(a) 第100炮褶积运算的偏移结果(x=5475 m)； (b) 第100炮近似公式的偏移结果(x=5475 m)； (c) 第150炮褶积运算的偏移结果(x=6725 m)； (d) 第150炮近似公式的偏移结果(x=6725 m)Fig.4　(a) The No.100 migrated shot record obtained by using the the exact method (x=5475 m)； (b) The No.100 migrated shot record obtained by using the approximated method (x=5475 m)； (c) The No.150 migrated shot record obtained by using the the exact method (x=6725 m)； (d) The No.150 migrated shot record obtained by using the approximated method (x=6725 m)

图5为叠加后的偏移剖面(Δt=0.005s).可以看出，两种方法偏移结果一致，三个断层与两个背斜构造清晰可见，较浅背斜右侧翼，因为倾角过大，边界模糊.较深背斜中的小储层构造清晰可见.之前单炮记录(近似公式，图4b、4d)中的噪声，经叠加后被有效压制.表1给出了两种方法的计算效率(CPU耗时)对比，在使用MPI进行多线程计算条件下，使用近似公式的计算时间约为近似前的0.2～0.25倍，计算效率得到了明显提升.

图5　(a) 褶积运算的叠加剖面；(b)近似公式的叠加剖面Fig.5　(a) Stacked image obtained by using the exact method； (b) Stacked image obtained by using the approximated method

第二个算例为Sigsbee2A模型.图6为速度模型.横向网格节点数为2133，网格间隔为11.43m，纵向网格节点数为1201，网格间隔为7.62m.地震记录共500炮，最大道数为348道，炮间距和道间距均为22.86m.

图6　Sigsbee2A速度模型Fig.6　Stratigraphic interval velocity model for Sigsbee2A dataset

第100炮第150炮叠加剖面近似前21.619.63312.5近似后4.23.9705.6

图7为给出了两种算法的叠加剖面(Δt=0.005 s).通过对比可以看出，两种方法对盐丘的陡倾边界以及盐体下断层构造与绕射体进行了较好的成像，但盐体左翼及右翼最深处的下方构造，近似后的成像质量相比近似前，信噪比要低，干扰更大一点，这是由于盐下区域照明不足，delta波包叠加不充分，噪声没有得到有效压制.最后给出两种方法的计算效率(CPU耗时)对比，在使用MPI进行多线程计算的条件下，近似前的偏移时间为119小时15分，近似后的偏移时间为25小时13分.

图7　(a) 褶积运算的叠加剖面；(b)近似公式的叠加剖面Fig.7　7(a) Stacked image obtained by using the exact method； (b) Stacked image obtained by using the approximated method

6结论

本文将delta波包应用到偏移领域，给出了基于delta波包叠加的深度偏移算法，并根据褶积算子的性质，进一步将褶积运算简化为乘积形式.在近似公式中，衰减项A衰由Δt和τI的比值β所决定，本文建议Δt的取值为几个时间采样间隔，并最好保证A衰与实际振幅衰减程度相近.最后，通过Marmousi和Sigsbee2A两个数值实例验证了本文提出的基于delta波包叠加的深度偏移算法，保留了高斯波束偏移的优点，对复杂介质具有良好的成像能力，并且近似算法的偏移结果与近似前基本一致，但计算时间约为近似前的0.2～0.25倍，计算效率得到显著提升.该偏移算法适用于多种反射地震采集方式(如共炮，共炮检距等)，易于进行面向目标成像，可拓展到3维地震数据偏移，在勘探地震学领域有着广泛的应用前景.同时，该算法仍存在诸多有待解决的问题，如Δt的影响，盐下区域成像改善等，需要进一步研究与讨论.

References

Albertin U, Yingst D, Kitchenside P, et al. 2004. True-amplitude beam migration. ∥74th Annual International Meeting, SEG. Expanded Abstracts, 398-401.

Born M, Wolf E. 1999. Principles of Optics: Electromagnetic Theory of Propagation, Interference and Diffraction of Light. Cambirdge: Cambridge University Press.

Compile Group of Mathmatics Handbook. 1979. Mathmatics handbook(in Chinese). Beijing: People′s Education Press.

Gray S H. 2005. Gaussian beam migration of common-shot records. Geophysics, 70(4): S71-S77. Gray S H, Bleistein N. 2009. True-amplitude Gaussian-beam migration. Geophysics, 74(2): S11-S23. Hill N R. 1990. Gaussian beam migration. Geophysics, 55(11): 1416-1428. Hill N R. 2001. Prestack Gaussian-beam depth migration. Geophysics, 66(4): 1240-1250. Kogelnik H. 1965. On the propagation of Gaussian beams of light through lenslike media including those with a loss or gain variation. Applied Optics, 4(12): 1562-1569.

Nowack R L, Sen M K, Stoffa P L. 2003. Gaussian beam migration for sparse common-shot and common-receiver data. ∥73th Annual International Meeting, SEG. Expanded Abstracts, 1114-1117.

Popov M M. 1982. A new method of computation of wave fields using Gaussian beams. Wave Motion, 4(1): 85-97.

Popov M M, Semtchenok N M, Popov P M, et al. 2010. Depth migration by the Gaussian beam summation method. Geophysics, 75(2): S81-S93. Sun J G. 2012. The history, the state of the art and the future trend of the research of Kirchhoff-type migration theory: a comparison with optical diffraction theory, some new results and new understanding, and some problems to be solved. Journal of Jilin University (Earth Science Edition) (in Chinese), 42(5): 1521-1552.

Yue Y B. 2011. Study on Gaussian beam migration methods in complex medium [Ph. D. thesis] (in Chinese). Qingdao: China University of Petroleum(East China).Zhang Y L. 2003. Engineering Mathmatics: Integral Transformation (in Chinese). Beijing: Higher Education Press.

附中文参考文献

《数学手册》编写组. 1979. 数学手册. 北京: 人民教育出版社.

孙建国. 2012. Kirchhoff型偏移理论的研究历史、研究现状与发展趋势展望——与光学绕射理论的类比、若干新结果、新认识以及若干有待于解决的问题. 吉林大学学报: 地球科学版, 42(5): 1521-1552.

岳玉波. 2011. 复杂介质高斯束偏移成像方法研究[博士论文]. 青岛: 中国石油大学(华东).

张元林. 2003. 工程数学·积分变换. 北京: 高等教育出版社.

(本文编辑胡素芳)

基金项目国家自然科学基金(41274120)，国家科技专项(Sinoprobe 09-01)，国家自然科学基金(41404085，41504084)资助.

作者简介石秀林, 男, 1985年生, 吉林大学地球探测与信息技术专业在读博士, 主要从事地震偏移成像研究.E-mail: xiulinshi@126.com *通讯作者孙建国, 男, 1956年生, 博士, 教授, 主要从事波动理论与成像技术、地震资料处理方法与解释技术等方面的教学和研究工作．E-mail: sun_jg@jlu.edu.cn

doi:10.6038/cjg20160727 中图分类号P631

收稿日期2015-10-10，2016-03-09收修定稿

The time-domain depth migration by the summation of delta packets

SHI Xiu-Lin, SUN Jian-Guo*, SUN Hui, LIU Ming-Chen, LIU Zhi-Qiang

GeoExplorationofScienceandTechnology,JilinUniversity,Changchun130026,China

AbstractDelta packets may serve as building blocks of a wavefield in space-time domain, just as Gaussian beams do in space-frequency domain. This paper will present a depth migration using the summation of delta packets and develop an approximate method that can speed up the convolution between the data and the delta packet. By definition, a delta packet is the inverse Fourier transformation of Gaussian beam under consideration. Thus, the Green′s function in the space-time domain can be asymptotically expressed as an integral superposition of delta packets. In the Rayleigh integral, the choice of Green′s function determines the migration method. Choosing a Green′s function of the delta packet summation form will lead to a new depth migration in space-time domain. The space-time formulation contains a convolution with respect to time between the filtered seismogram and the delta packet, which makes the new depth migration very time-consuming. Here the convolution integration is rewritten in terms of Taylor expansion, and is further reduced to product formula.

In the approximate expression of convolution, the time-shifted seismogram only uses multiplication by the term Aattenuation, so that the computational cost is much lower than before approximation. The term Aattenuationindicates attenuation degree of amplitudes, and depends on ratio of parameter Δt (the half width of the integration interval) and imaginary part of the complex travel time. When the value of Δt takes no more than several time sampling intervals, such as Δt=0.005 s, there will be a satisfactory approximation. Numerical tests using the Marmousi and the Sigsbee2A models show that the migration quality for the approximated method is comparable to the one with the exact method. However, the computation time can reduce by about 75%. The depth migration using the summation of delta packets can be recognized as a flexible and effective migration technique. And it is also easily extended to 3-D depth migration and suitable for different kinds of acquisition geometries.

KeywordsDelta packet; Gaussian beam; Convolution; Numerical approximation; Time-domain depth migration

石秀林, 孙建国, 孙辉等. 2016. 基于delta波包叠加的时间域深度偏移. 地球物理学报，59(7):2641-2649,doi:10.6038/cjg20160727.

Shi X L, Sun J G, Sun H, et al. 2016. The time-domain depth migration by the summation of delta packets. Chinese J. Geophys. (in Chinese),59(7):2641-2649,doi:10.6038/cjg20160727.