地球物理学报  2017, Vol. 60 Issue (3): 1118-1131   PDF    
引用本文
高正辉, 孙建国, 孙章庆, 等. 2017. 共炮检距高斯波束偏移的一种快速实现算法. 地球物理学报, 60(3): 1118-1131, doi: 10.6038/cjg20170323.  
GAO Z H, SUN J G, SUN Z Q, et al. 2017. A fast algorithm for common-offset Gaussian beam migration. Chinese J. Geophys. (in Chinese), 60(3): 1118-1131, doi: 10.6038/cjg20170323.  
共炮检距高斯波束偏移的一种快速实现算法
高正辉1, 孙建国1 , 孙章庆1, 孙煦1,2, 石秀林1, 刘志强1     
1. 吉林大学地球探测科学与技术学院, 长春 130026;
2. 吉林大学工程训练中心, 长春 130026
收稿日期 2016-04-28, 2017-01-09 收修定稿
基金项目: 国家自然科学基金（41274120、41404085、41504084）及国家科技专项（Sinoprobe09-01）联合资助
作者简介: 高正辉, 男, 1991年生, 吉林大学地球探测与信息技术专业在读博士, 主要从事地震波正演及偏移成像研究.E-mail:gao_jlu_stu@163.com
通讯作者: 孙建国, 男, 1956年生, 博士, 教授, 主要从事地下波动理论与成像技术、地震资料处理方法与解释技术、钻孔电磁波法、地球物理中的天线问题以及岩石物理学和可视化技术在地球物理中的应用等方面的教学和研究工作.E-mail:sun_jg@jlu.edu.cn
摘要: 与共炮高斯波束偏移相比，共炮检距高斯波束偏移具有直接抽取炮检距域共成像点道集的优势.过去，共炮检距高斯波束偏移以损失成像精度的代价采用最速下降法来降低积分的维数，从而提高计算效率.但经过最速下降近似简化的偏移公式仍是频率域的，需要在每个频点进行计算.为此，本文提出一种快速实现算法来避免采用最速下降法.本文通过分析一个水平层状速度模型的偏移过程和Marmousi速度模型的成像结果来检验不同插值方法对快速实现算法的成像精度和计算效率的影响，并建议采用二维三次卷积插值方法.同时本文在Marmousi速度模型下验证了快速实现算法相对于最速下降法在成像精度和计算效率上的优势.此外，本文将采用二维三次卷积插值的快速实现算法应用于Sigsbee2A模型并获得了清晰的盐下图像.
关键词: 共炮检距      高斯波束偏移      快速实现算法      插值方法     
A fast algorithm for common-offset Gaussian beam migration
GAO Zheng-Hui1, SUN Jian-Guo1, SUN Zhang-Qing1, SUN Xu1,2, SHI Xiu-Lin1, LIU Zhi-Qiang1     
1. College of GeoExploration Science and Technology, Jilin University, Changchun 130026, China;
2. Engineering Training Center, Jilin University, Changchun 130026, China
Abstract: Compared with the common-shot Gaussian beam migration, the common-offset Gaussian beam migration has the advantage in direct extraction of offset-domain common-image gathers. In the past, the common-offset Gaussian beam migration adopted the steepest descent method to reduce the dimension of the integrals at the cost of losing precision and to increase computation efficiency. However, the migration formula simplified by the steepest descent method is still in the frequency domain, and it has to be evaluated at each frequency. To solve this problem, we develop a fast algorithm to avoid adopting the steepest descent method. We test the effects of different interpolation methods on the imaging accuracy and the computational efficiency of the fast algorithm through analyzing the migration procedures on a horizontally layered velocity model and the imaging results of the Marmousi velocity model. We suggest to adopt the bicubic convolution interpolation method. The superiority of this fast algorithm on the imaging accuracy and the computational efficiency compared with the steepest descent method is validated with the Marmousi velocity model. Furthermore, we apply this fast algorithm using the bicubic convolution interpolation to the Sigsbee2A model and obtain clear subsalt image.
Key words: Common-offset      Gaussian beam migration      Fast algorithm      Interpolation method     
1 引言

Hill (1990)提出了叠后高斯波束偏移.其通过局部倾斜叠加将地震记录分解成局部平面波，然后通过高斯波束来进行波场延拓.Hale (1992)通过比较叠后的Kirchhoff偏移、倾斜叠加偏移和高斯波束偏移指出高斯波束偏移灵活高效稳定的特点.Hill (2001)将高斯波束应用到叠前数据中，基于互相关成像条件推导出了共炮检距高斯波束偏移的公式，并在只考虑运动学特性的情况下采用最速下降法来进行计算.Gray (2005)基于Hill (2001)采用最速下降法计算共炮检距高斯波束偏移的方法，提出了几种效率和精度不同的计算共炮高斯波束偏移的实现策略.Gray和Bleistein (2009)分别基于互相关成像条件和反褶积成像条件推导出了二维和三维情况下共炮真振幅高斯波束偏移的公式，同时给出了采用最速下降法进行计算的方法.Popov等 (2010)通过高斯波束的叠加直接表示成像点处的格林函数，得到成像点处实际的正向波场和反向波场，再利用互相关成像条件来进行成像.Bleistein和Gray (2010)详细给出了采用最速下降法进行3D叠前高斯波束偏移中Hessian矩阵的计算方法.Popov等 (2011)指出Hill (2001)由于采用最速下降法进行计算，虽然在效率上会有提高，但同时也会引起振幅和相位的精度损失.Casasanta和Gray (2015)分析了在不同射线参数域下采用最速下降法进行计算对于像场能量的影响，指出中心点-炮检距射线参数域能够获得更好的成像结果但也会损失一定的精度.El Yadari (2015)在处理真振幅矢量声波高斯波束偏移时指出采用最速下降法会对成像精度造成损失，并建议避免采用最速下降法进行计算.国内，岳玉波等 (2012)、黄建平等 (2015)及高成等 (2015)也对高斯波束偏移进行了相关研究.

所谓最速下降法，即将关于射线参数的积分从炮点-检波点射线参数域变换到中心点-炮检距射线参数域，进而通过计算积分的鞍点 (走时的虚部最小值所对应的点) 的值来替代整个积分的计算，以达到降低维数的目的.综上，以往的叠前高斯波束偏移大多采用最速下降法进行计算，而这虽然会提高计算效率，但却以损失成像精度为代价 (Popov M M and Popov P M，2011；Casasanta and Gray, 2015；El Yadari，2015).此外，经过最速下降近似简化的偏移公式仍然是频率域的，需要在每个频点进行计算.Popov等 (2010)通过高斯波束的叠加直接表示成像点处的格林函数，得到成像点处实际的正向波场和反向波场，但需要计算每个频率下的格林函数，虽然不会造成精度损失，但效率较低.共炮高斯波束偏移由于不包含明确的炮检距信息，不能直接抽取炮检距域共成像点道集.而共炮检距高斯波束偏移可以直接抽取炮检距域共成像点道集，从而进行偏移速度分析和AVO分析，得到高质量的地下成像结果.为此，对于共炮检距高斯波束偏移，本文提出一种快速实现算法来避免采用最速下降法.为了提出该算法，本文通过改变积分的顺序，推导了共炮检距高斯波束偏移的公式.接着，本文将推导出的公式的两个最内层积分看成是总复走时的实部与虚部的两个二维连续函数.该算法首先确定总复走时的实部与虚部的变化范围和采样间隔；其次，在实部与虚部的所有采样点处构建两个最内层积分的值对应的两个二维数组；最后，根据成像点处实际的总复走时的实部与虚部进行插值得到两个最内层积分的结果.这样，成像精度的损失只受所采用的插值方法精度的影响.本文通过分析一个水平层状速度模型的偏移过程和Marmousi速度模型的成像结果来检验不同插值方法对于快速实现算法的成像精度和计算效率的影响.同时本文在Marmousi速度模型下比较了快速实现算法和最速下降法的成像精度及计算效率.然后，本文将采用二维三次卷积插值的快速实现算法应用于Sigsbee2A模型.

2 理论基础

为了提出快速实现算法，我们通过改变积分的顺序，推导出共炮检距高斯波束偏移公式.

首先根据Gray和Bleistein (2009)给出的未经最速下降法近似的共炮高斯波束偏移公式 (Gray and Bleistein (2009)中的式 (10)) 及，，同时考虑到及，将其变换到共炮检距域，

(1)

式中，x表示成像点位置，ω表示角频率，ωr表示参考频率，L表示高斯窗中心的位置，ΔL表示高斯窗中心的间隔，w0表示高斯波束初始波束半宽，xm表示中心点坐标，h表示半炮检距，L-h和L+h表示单个高斯窗内高斯波束中心射线的出射位置，和分别表示从L-h和L+h处发出中心射线的水平慢度分量，，和分别表示从L-h和L+h处发出的高斯波束，*表示复共轭，i2=-1，

(2)

表示局部倾斜叠加.

Hill (2001)将式 (1) 从域变换到 (pxm, pxh) 域，进而通过最速下降近似将关于射线的二维积分简化成一维积分 (Gray (2005)，Gray和Bleistein (2009)，Bleistein和Gray (2010)，Casasanta和Gray (2015)及El Yadari (2015)也采用同样的方法).但简化后的积分仍然是频率域，需要在每个频点进行计算.这里我们把代入式 (1)，并改变积分的顺序得

(3)

式中，A和T分别表示高斯波束的复振幅和复走时 (为了表示方便，消去GB).

令

(4)

则

(5)

令A=AL-hAL+h，T=TL-h+TL+h，则

(6)

根据Hale (1992)，因为需要得到实值的成像结果，所以在频率域中需要保证每个量共轭对称，即需要

(7)

式中，AR和AI分别表示A的实部和虚部，TR和TI分别表示T的实部和虚部.把式 (7) 代入式 (6) 并整理得

(8)

令

(9)

(10)

则

(11)

将式 (11) 代入式 (4) 便可以得到我们所要推导的共炮检距高斯波束偏移公式

(12)

分析式 (12) 的最内层积分I1a(x; L;pxm; TI; TR) 和I1b(x; L;pxm; TI; TR)(式 (9) 和式 (10)) 可以看出，它们都是关于频率的积分，对于两条成像射线和成像区域内的每一个成像点 (如图 1)，都要对它们进行重复计算，而这会严重降低计算的效率，所以本文的核心是通过一种快速实现算法来计算I1a(x; L;pxm; TI; TR) 和I1b(x; L;pxm; TI; TR).

图 1 单个高斯窗内的成像示意图 (重叠区域表示成像射线x和的成像区域) Fig. 1 Sketch of imaging in one Gaussian window. The overlapping region represents the imaging region of the imaging rays and
3 快速实现算法的具体策略

如图 1，对于两条成像射线和成像区域的某一个成像点x，TR和TI是一对单值，需要对关于频率的积分I1a(x; L;pxm; TI; TR) 和I1b(x; L;pxm; TI; TR) 直接离散求和进行计算.而对于成像区域内的所有成像点，都需要重复同样的步骤，而这会严重降低整个偏移的效率.为了方便起见，我们把这种对I1a(x; L;pxm; TI; TR) 和I1b(x; L;pxm; TI; TR) 直接离散求和的方法叫做离散求和算法.

下面我们通过一种快速实现算法来避免对I1a(x; L;pxm; TI; TR) 和I1b(x; L;pxm; TI; TR) 直接离散求和.分析I1a(x; L;pxm; TI; TR) 和I1b(x; L; pxm; TI; TR) 可以看出，它们的积分核相差相位，所以首先只考虑I1a(x; L;pxm; TI; TR) 的计算问题.再分析I1a(x; L;pxm; TI; TR) 可以看出，经过关于ω的积分之后，该式与ω无关.而对于两条成像射线和成像区域的所有成像点，TR和TI是连续变化的，且有其变化范围.因此对于两条成像射线和成像区域的所有成像点，可以将I1a(x; L;pxm; TI; TR) 看成是TR和TI的二维连续函数，根据TR和TI的变化范围，间隔采样，计算出所有采样点处的积分值，再根据成像点处实际的TR和TI进行插值得到该处的积分值.

首先讨论如何确定TI的变化范围及其采样间隔.考虑到当|ω|TI>5.0时，I1a(x; L;pxm; TI; TR) 的值可以忽略不计，所以我们可以根据参考频率ωr来确定TI的最大值为TImax=5.0/|ωr|.由此，TI的变化范围为[0, 5.0/ωr].TI的采样间隔ΔTI可以由ΔTI=TImax/(NI-1)(NI表示TI的采样点个数) 来确定.

接下来讨论如何确定TR的变化范围及其采样间隔.考虑到本文中的偏移公式是基于互相关成像条件推导而来，所以I1a(x; L;pxm; TI; TR) 可以写成

(13)

式中t表示延迟时间.由式 (9) 可知

(14)

由Fourier变换的位移定理可知

(15)

即对于两条成像射线和成像区域的某一个成像点，根据挑选出Ds(L, pxm, ω)，对每一个频率乘以iω及exp[-|ω|TI]，进行快速Fourier变换得到时间序列I1a(x; L;pxm; TI; t)，再根据TR的大小移动整个时间序列，最后选择t=0时对应的积分值.而在实际计算的过程中，我们不需要对整个时间序列做真正的移动，只需要根据TR的大小，在时间序列中直接挑选出其所对应的值，即为t=0时的积分值.但是对于某一个成像点，TR往往不在时间序列的采样点上，这时可以通过对TR附近时间序列采样点上的积分值进行插值来获取TR处的积分值.我们将时间序列的最大采样时间用tmax来表示，则对于成像区域的所有成像点，TR的变化范围为[0, tmax].因此对于这些成像点，我们都可以根据TR的值在时间序列I1a(x; L;pxm; TI; t) 上进行插值，从而得到该成像点处的积分值.假设时间序列I1a(x; L;pxm; TI; t) 的采样间隔为Δt，则时间序列I1a(x; L;pxm; TI; t) 上的点可以看成是TR=0, Δt, 2Δt，…，(n-1)Δt(n为时间序列I1a(x; L;pxm; TI; t) 的采样点个数) 时的采样值.由此，TR的采样间隔ΔTR即为Δt.Δt可以由计算的最大频率ωmax通过Δt=2π/ωmax来确定，从而ΔTR=Δt=2π/ωmax.

在确定完TR和TI的变化范围及它们的采样间隔之后，根据挑选出Ds(L, pxm, ω)，计算出每个频率对应的iωDs(L, pxm, ω)，对其乘以相应的TI采样点的衰减exp[-|ω|TI，再通过快速Fourier变换得到在每个TR采样点的值.由此我们构建出在TR和TI采样点处的一个二维数组 (如图 2).最后，根据某一成像点处实际的TR和TI进行插值就可以得到该成像点的I1a(x; L;pxm; TI; TR) 积分的结果.

图 2 时构建的关I1a(x; L;pxm; TI; TR) 的二维插值数组 Fig. 2 Constructed two-dimensional interpolation array ofI1a(x; L;pxm; TI; TR) when

在计算I1b(x; L;pxm; TI; TR) 时，只需要在相应的频率上乘以-isgn (ω)，即将上述步骤中的Fourier变换变成Hilbert变换，就可以按照上述方法获得其结果.

4 数值分析 4.1 水平层状速度模型

如图 3，速度模型的横向网格点个数为401，纵向网格点个数为201，纵向和横向网格间距均为5 m.界面深度为500 m，上层速度为2000 m·s-1，下层速度为3000 m·s-1.图 4为我们根据该速度模型得到的炮检距为100 m的地震记录，时间采样点个数为600，时间采样间隔为0.002 s，第一道的位置为50 m，道间距为10 m，共181道.

图 3 水平层状速度模型 Fig. 3 Horizontally layered velocity model
图 4 共炮检距地震记录 (炮检距为100 m) Fig. 4 Common-offset seismic record (offset equals 100 m)

按照Hill (2001)中给出的经验公式，我们将地震记录分成16个相互重叠的高斯窗 (如图 4)，并选取其中第9个高斯窗 (虚线) 内的数据来验证上述快速实现算法的精度.设定参考频率ωr=10×2π rad·s-1，虚走时的采样个数为NI=101，则虚走时的采样间隔ΔTI=0.00796 s.由地震记录的时间采样点个数和时间采样间隔，可以计算出地震记录中的最大频率为500×2π rad·s-1，频率间隔Δω为0.8333×2π rad·s-1，频率采样点个数为600.在实际计算中，为了加密TR的采样密度，我们将最大频率ωmax设为地震记录中最大频率的四倍，即2000×2π rad·s-1，则频率采样点个数为2400.当频率大于500×2π rad·s-1时，该频率对应的值为零.由此，根据ΔTR=Δt=2π/ωmax计算出TR的采样间隔ΔTR=0.0005 s，采样点个数NR=2400.从而，根据快速实现算法，对于每一个方向的局部平面波，我们分别构建出I1a(x; L;pxm; TI; TR) 和I1b(x; L;pxm; TI; TR) 在TR和TI采样点处的二维数组.挑选pxm=-2.859×10-5s·m-1方向的I1a(x; L;pxm; TI; TR) 的数据来检验不同插值方法的精度，插值点在每个网格的中心位置.这里我们分别采用二维线性插值、二维三次卷积插值和二维四次卷积插值 (韩复兴，2009).将对I1a(x; L;pxm; TI; TR) 直接离散求和的结果作为计算插值结果误差的标准.图 5a，5d，5g分别为二维线性插值、二维三次卷积插值和二维四次卷积插值得到的I1a(x; L;pxm; TI; TR) 在插值点的值.图 5b，5e，5h分别为相应插值结果的绝对误差.图 5c，5f，5i为相应插值结果的相对误差.从图 5可以看出二维线性插值的精度最低，二维四次卷积插值的精度最高，而二维三次卷积插值接近二维四次卷积插值的精度水平.

图 5 pxm=-2.859×10-5s·m-1时二维线性插值、二维三次卷积插值和二维四次卷积插值在插值点的插值结果 (a，d，g)，绝对误差 (b，e，h) 及相对误差 (c，f，i) Fig. 5 Interpolation results (a, d, g), absolute errors (b, e, h), relative errors (c, f, i) of two-dimensional linear interpolation, bicubic convolution interpolation, and two-dimensional quartic convolution interpolation at the sampling points when pxm=-2.859×10-5s·m-1

图 6为第9个高斯窗内数据的偏移结果，其中图 6a-6d分别为离散求和算法及采用二维线性插值、二维三次卷积插值和二维四次卷积插值的快速实现算法的偏移结果.图 7a，7c，7e分别为采用三种插值方法的快速实现算法在图 6虚线方框内的成像结果的绝对误差.图 7b，7d，7f为相应的相对误差.通过分析图 7可以看出，采用二维四次卷积插值的快速实现算法的精度最高，采用二维三次卷积插值的快速实现算法次之，而采用二维线性插值的快速实现算法的精度最低.在相同的计算条件下，针对第9个高斯窗内的数据，采用离散求和算法，采用二维线性插值的快速实现算法，采用二维三次卷积插值的快速实现算法及采用二维四次卷积插值的快速实现算法的计算机运行时间分别为11.40 min，2.02 s，2.56 s和3.80 s.图 8a-8d分别为离散求和算法及采用二维线性插值、二维三次卷积插值和二维四次卷积插值的快速实现算法的整个共炮检距地震记录的偏移结果.在相同的计算条件下，相应的计算机运行时间分别为154.28 min，0.46 min，0.59 min和0.87 min.图 9a，9c，9e分别为采用三种插值方法的快速实现算法的偏移结果的每道最大幅值的绝对误差.图 9b，9d，9f分别为相应的相对误差.分析图 9可以得出和图 5及图 7相同的结论，即采用三种插值方法的快速实现算法的精度水平由低到高依次为二维线性插值、二维三次卷积插值、二维四次卷积插值.

图 6 第9个高斯窗内数据的偏移结果 (a) 离散求和算法；(b) 采用二维线性插值的快速实现算法；(c) 采用二维三次卷积插值的快速实现算法；(d) 采用二维四次卷积插值的快速实现算法. Fig. 6 Migration results of the data in the 9th Gaussian window (a) Discrete summation method; (b) Fast algorithm using two-dimensional linear interpolation; (c) Fast algorithm using bicubic convolution interpolation; (d) Fast algorithm using two-dimensional quartic convolution interpolation.
图 7 图 6中虚线方框内偏移结果的绝对误差 (a, c, e) 和相对误差 (b, d, f) (a)(b) 对应图 6b; (c)(d) 对应图 6c; (e)(f) 对应图 6d. Fig. 7 Absolute errors (a, c, e) and relative errors (b, d, f) of the migration results in the dotted box in Fig. 6 (a) and (b) correspond to Fig. 6b; (c) and (d) correspond to Fig. 6c; (e) and (f) correspond to Fig. 6d.
图 8 整个炮检距道集数据的偏移结果 (a) 离散求和算法；(b) 采用二维线性插值的快速实现算法；(c) 采用二维三次卷积插值的快速实现算法；(d) 采用二维四次卷积插值的快速实现算法. Fig. 8 Migration results of whole common-offset gather data (a) Discrete summation method; (b) Fast algorithm using two-dimensional linear interpolation; (c) Fast algorithm using bicubic convolution interpolation; (d) Fast algorithm using two-dimensional quartic convolution interpolation.
图 9 采用快速实现算法的整个炮检距道集数据偏移结果的每道最大幅值的绝对误差 (a, c, e) 和相对误差 (b, d, f)(a)(b) 二维线性插值；(c)(d) 二维三次卷积插值；(e)(f) 二维四次卷积插值. Fig. 9 Absolute errors (a, c, e) and relative errors (b, d, f) of the maximum amplitude of each trace in the migration results of the fast algorithm for the whole common-offset gather data (a) and (b) Two-dimensional linear interpolation; (c) and (d) Bicubic convolution interpolation; (e) and (f) Two-dimensional quartic convolution interpolation.
4.2 Marmousi速度模型

图 10为Marmousi速度模型，模型横向网格点个数为737，横向网格间距为12.5 m，纵向网格点个数为750，纵向网格间距为4 m.基于Marmousi速度模型模拟的地震记录包含240炮记录，放炮方式为中间放炮，第一炮的位置在3000 m处，最后一炮的位置在8975 m处，炮间距为25 m.每一炮的最小炮检距为-2575 m，最大炮检距为2575 m，道间距为25 m，共207道，采样间隔为2 ms，采样点个数为1500.

图 10 Marmousi速度模型 Fig. 10 Marmousi velocity model

根据炮检距信息，将共炮道集抽成共炮检距道集，共抽出207个共炮检距道集，炮检距范围为-2575~2575 m，道间隔为25 m.图 11为对所有共炮检距道集的偏移结果.图 12a，12c，12e分别为x=3750 m处采用三种插值方法的快速实现算法的成像结果的绝对误差.图 12b，12d，12f为相应的相对误差.图 13a，13c，13e为在x=7500 m处采用三种插值方法的快速实现算法的成像结果的绝对误差.图 13b，13d，13f为相应的相对误差.从图 12和图 13可以看出采用二维四次卷积插值的快速实现算法的精度最高，采用二维三次卷积插值的快速实现算法的精度次之，采用二维线性插值的快速实现算法的精度最差.在相同的计算条件下，离散求和算法的计算机运行时间为24.5 h，采用二维线性插值的快速实现算法的计算机运行时间为3.27 min，采用二维三次卷积插值的快速实现算法的计算机运行时间为4.26 min，采用二维四次卷积插值的计算机运行时间为5.97 min.

图 11 Marmousi速度模型的偏移结果 (a) 离散求和算法；(b) 采用二维线性插值的快速实现算法；(c) 采用二维三次卷积插值的快速实现算法；(d) 采用二维四次卷积插值的快速实现算法. Fig. 11 Migration results of the Marmousi velocity model (a) Discrete summation method; (b) Fast algorithm using two-dimensional linear interpolation; (c) Fast algorithm using bicubic convolution interpolation; (d) Fast algorithm using two-dimensional quartic convolution interpolation.
图 12 x=3750 m处分别采用三种插值方法的快速实现算法的成像结果的绝对误差 (a，c，e) 和相对误差 (b，d，f)(a)(b) 二维线性插值；(c)(d) 二维三次卷积插值；(e)(f) 二维四次卷积插值. Fig. 12 Absolute errors (a, c, e) and relative errors (b, d, f) of the imaging results of the fast algorithm using three interpolation methods, respectively, at x=3750 m (a) and (b) Two-dimensional linear interpolation; (c) and (d) Bicubic convolution interpolation; (e) and (f) Two-dimensional quartic convolution interpolation.
图 13 x=7000 m分别采用三种插值方法的快速实现算法的成像结果的绝对误差 (a，c，e) 和相对误差 (b，d，f)(a)(b) 二维线性插值；(c)(d) 二维三次卷积插值；(e)(f) 二维四次卷积插值. Fig. 13 Absolute errors (a, c, e) and relative errors (b, d, f) of the imaging results of the fast algorithm using the three interpolation method, respectively, at x=7000 m (a) and (b) Two-dimensional linear interpolation; (c) and (d) Bicubic convolution interpolation; (e) and (f) Two-dimensional quartic convolution interpolation.

综上，可以看出采用二维三次卷积插值和二维四次卷积插值的快速实现算法精度均较高，而采用二维线性插值的快速实现算法的精度最低.但由于二维线性插值的算子长度最短，所以其具有较高的计算效率.综合考虑成像精度和计算效率，我们建议采用二维三次卷积插值作为快速实现算法的插值方法.

4.3 Marmousi速度模型下同最速下降法的比较

为了验证本文算法相对于最速下降法的优势，我们给出在Marmousi速度模型下最速下降法的偏移结果.图 14为炮检距为-1325 m的共炮检距道集的偏移结果.图 15为炮检距为200 m的共炮检距道集的偏移结果.图 16为所有共炮检距道集的偏移结果.在同上例相同的计算条件下，最速下降法的计算机运行时间为159.7 min，大约是采用二维三次卷积插值的快速实现算法的34.5倍.通过比较单个共炮检距道集的偏移结果可以看出，由于采用最速下降法，只有走时的虚部最小值处的能量才对最终的成像结果有贡献，因此会造成成像精度的损失.通过比较所有共炮检距道集的偏移结果可以看出，由于不同共炮检距道集成像结果的相互弥补，保证了最速下降法的最终结果对于构造成像的准确性，但是在局部会造成成像能量的异常.在提高计算效率方面，最速下降法虽然能够降低积分的维数，但是仍需要在每个频点进行计算.快速实现算法相当于通过插值来代替关于频率的积分的计算，而关于频率的循环次数远比关于射线参数的循环次数多.由此，快速实现算法的成像精度和计算效率均高于最速下降法.

图 14 Marmousi速度模型下炮检距为-1325 m的共炮检距道集的偏移结果 (a) 最速下降法；(b) 采用二维三次卷积插值的快速实现算法. Fig. 14 Migration results of the common-offset gather with offset -1325 m under the Marmousi velocity model (a) Steepest descent method; (b) Fast algorithm using bicubic convolution interpolation.
图 15 Marmousi速度模型下炮检距为200 m的共炮检距道集的偏移结果 (a) 最速下降法；(b) 采用二维三次卷积插值的快速实现算法. Fig. 15 Migration results of the common-offset gather with offset 200 m under the Marmousi velocity model (a) Steepest descent method; (b) Fast algorithm using bicubic convolution interpolation.
图 16 Marmousi速度模型下所有共炮检距道集的偏移结果 (a) 最速下降法；(b) 采用二维三次卷积插值的快速实现算法. Fig. 16 Migration results of all the common-offset gathers under the Marmousi velocity model (a) Steepest descent method; (b) Fast algorithm using bicubic convolution interpolation.
5 Sigsbee2A速度模型下的应用

图 17a为Sigsbee2A速度模型，横向和纵向宽度分别为24.38 km和9.15 km，相应的网格间距分别为7.62 m和11.43 m.基于Sigsbee2A速度模型模拟的地震记录包括500炮记录，放炮方式为左边放炮，第一炮的位置274.32 m，炮间距为45.72 m.单炮最多包含348道地震记录.我们将共炮道集抽成共炮检距道集，共抽出348个，炮检距的范围为0~7932.42 m，道间隔为22.86 m.模型中部高速盐体的存在会对盐下构造的成像造成影响.我们将采用二维三次卷积插值的快速实现算法应用于该模型来测试其对于盐下构造的成像能力.图 17b为我们得到的所有共炮检距道集的偏移结果.从图中可以看出盐下构造成像清晰，在一定程度上减弱了上覆盐体的影响.图 18为在图 17b中所示的四个位置处抽取的炮检距域共成像点道集.从图 18可以看出共成像点道集的同相轴平直，验证了本文共炮检距高斯波束偏移的快速实现算法的准确性.

图 17 Sigsbee2A模型 (a) 速度模型；(b) 采用二维三次卷积插值的快速实现算法的偏移结果. Fig. 17 Sigsbee2A model (a) Velocity model; (b) Migration result of the fast algorithm using bicubic convolution interpolation.
图 18 在图 17b中所示的四个位置抽取的炮检距域共成像点道集 Fig. 18 Offset-domain common-image gathers extracted at four locations shown in Fig. 17b
6 结论

本文提出了共炮检距高斯波束偏移的一种快速实现算法来避免采用最速下降法.为了提出该算法，本文将推导出的共炮检距高斯波束偏移公式的两个最内层积分看成是总复走时的实部与虚部的两个二维连续函数.该算法首先确定总复走时的实部与虚部的变化范围和采样间隔；其次，在相应的实部与虚部的所有采样点处构建两个最内层积分的值对应的两个二维数组；最后，根据成像点处实际的总复走时的实部与虚部进行插值得到两个最内层积分的结果.对于水平层状速度模型的偏移过程和Marmousi速度模型的成像结果的分析表明插值方法的选择对快速实现算法成像精度和计算效率会产生重要影响.综合考量成像精度和计算效率的要求，本文建议采用二维三次卷积插值.此外，本文在Marmousi速度模型下验证了快速实现算法相对于最速下降法在成像精度和计算效率方面的优势.采用二维三次卷积插值的快速实现算法在Sigsbee2A速度模型下的成像结果表明了该算法对于盐下构造良好的成像能力.

参考文献
Bleistein N, Gray S H. 2010. Amplitude calculations for 3D Gaussian beam migration using complex-valued traveltimes. Inverse Problems, 26(8). DOI:10.1088/0266-5611/26/8/085017
Casasanta L, Gray S H. 2015. Converted-wave beam migration with sparse sources or receivers. Geophysical Prospecting, 63(3): 534-551. DOI:10.1111/1365-2478.12226
El Yadari N. 2015. True-amplitude vector-acoustic imaging:Application of Gaussian beams. Geophysics, 80(1): S43-S54. DOI:10.1190/geo2014-0316.1
Gao C, Sun J G, Qi P, et al. 2015. 2-D Gaussian-beam migration of common-shot records in time domain. Chinese J. Geophys. (in Chinese), 58(4): 1333-1340. DOI:10.6038/cjg20150420
Gray S H. 2005. Gaussian beam migration of common-shot records. Geophysics, 70(4): S71-S77. DOI:10.1190/1.1988186
Gray S H, Bleistein N. 2009. True-amplitude Gaussian-beam migration. Geophysics, 74(2): S11-S23. DOI:10.1190/1.3052116
Hale D. 1992. Migration by the Kirchhoff, slant stack, and Gaussian beam methods. CWP-126, Center for Wave Phenomena, Colorado School of Mines, Golden, CO.
Han F X. 2009. On some computational Problems in wavefront construction method[Ph. D. thesis] (in Chinese). Changchun:Jilin University.
Hill N R. 1990. Gaussian beam migration. Geophysics, 55(11): 1416-1428. DOI:10.1190/1.1442788
Hill N R. 2001. Prestack Gaussian-beam depth migratiion. Geophysics, 66(4): 1240-1250. DOI:10.1190/1.1487071
Huang J P, Yuan M L, Li Z C, et al. 2015. The accurate beam migration method without slant stack under dual-complexity conditions and its application (I):Acoustic equation. Chinese J. Geophys. (in Chinese), 58(1): 267-276. DOI:10.6038/cjg20150124
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. DOI:10.1190/1.3361651
Popov M M, Popov P M. 2011. Comparison of the Hill's method with the seismic depth migration by the Gaussian beam summation method. Journal of Mathematical Sciences, 173(3): 291-298. DOI:10.1007/s10958-011-0251-8
Yue Y B, Li Z C, Qian Z P, et al. 2012. Amplitude-preserved Gaussian beam migration under complex topographic conditions. Chinese J. Geophys. (in Chinese), 55(4): 1376-1383. DOI:10.6038/j.issn.0001-5733.2012.04.033
高成, 孙建国, 齐鹏, 等. 2015. 2D共炮时间域高斯波束偏移. 地球物理学报, 58(4): 1333–1340. DOI:10.6038/cjg20150420
韩复兴. 2009.论波前构建法中的几个计算问题[博士论文].长春:吉林大学.
黄建平, 袁茂林, 李振春, 等. 2015. 双复杂条件下非倾斜叠加精确束偏移方法及应用I--声波方程. 地球物理学报, 58(1): 267–276. DOI:10.6038/cjg20150124
岳玉波, 李振春, 钱忠平, 等. 2012. 复杂地表条件下保幅高斯束偏移. 地球物理学报, 55(4): 1376–1383. DOI:10.6038/j.issn.0001-5733.2012.04.033