地球物理学报  2018, Vol. 61 Issue (4): 1413-1420   PDF    
引用本文
陈生昌, 周华敏. 2018. 地震数据的反射波动方程最小二乘偏移. 地球物理学报, 61(4): 1413-1420, doi: 10.6038/cjg2018L0025.  
Chen S C, Zhou H M. 2018. Least squares migration of seismic data with reflection wave equation. Chinese J. Geophys. (in Chinese), 61(4): 1413-1420, doi: 10.6038/cjg2018L0025.  
地震数据的反射波动方程最小二乘偏移
陈生昌1, 周华敏2     
1. 浙江大学地球科学学院, 杭州 310027;
2. 长江水利委员会长江科学院, 武汉 430010
收稿日期 2017-01-09, 2017-05-16 收修定稿
基金项目: 国家自然科学基金项目（41374001，41074133），国家重大专项（2016ZX05033002），中国石化地球物理重点实验室开放研究基金资助
第一作者简介: 陈生昌, 男, 1965年生, 教授, 博士生导师.主要从事勘探地球物理和计算地球物理研究.E-mail:chenshengc@zju.edu.cn
摘要：基于反射波动方程，本文提出了一种估计地下反射率分布的地震数据最小二乘偏移方法.高频近似下，非齐次的一次反射波动方程的源项是由反射率与入射波场的时间一阶导数相互作用产生的.根据反射波动方程，利用线性最小二乘反演方法由地震反射数据重建出地下产生反射波的反射源，再结合波场正演计算出的地下入射波场，得到地下反射率分布的估计.在地下反射源的线性最小二乘反演重建中，我们采用迭代求解方法，并以地震波的检波器单向地下照明强度作为最小二乘优化问题中Hessian矩阵的近似.
关键词: 地震数据      反射波动方程      最小二乘偏移      反射率      迭代重建     
Least squares migration of seismic data with reflection wave equation
CHEN ShengChang1, ZHOU HuaMin2     
1. School of Earth Sciences, Zhejiang University, Hangzhou 310027, China;
2. Changjiang River Scientific Research Institute of Changjiang Water Resources Commission, Wuhan 430010, China
Abstract: We proposed a least squares migration method of seismic data for the estimation of subsurface reflectivity distribution based on reflection wave equation. Under high frequency approximation, the source term of inhomogeneous primary reflection wave equation is generated by the interaction between the reflectivity and the first order derivative of incident wavefield with respect to time. According to the reflection wave equation, we used the linear least squares inversion to recover the subsurface reflection source generating reflection wave from seismic reflection data, then to obtain the estimation of subsurface reflectivity distribution from the recovered reflection source by combining the subsurface incident wavefield calculated by forward modeling. In the reconstruction of the subsurface reflection source with the linear least squares inversion, we adopted the iterative manner for the solution, and used the seismic wave illumination intensity of the receiver direction as the approximation of Hessian matrix of the least squares optimization problem.
Key words: Seismic data    Reflection wave equation    Least squares migration    Reflectivity    Iterative reconstruction    
0 引言

地震数据偏移成像是得到地下构造图像的主要方法技术，特别是地下复杂构造区的构造成像(Gray, et al., 2001；Yilmaz, 2001；Etgen, et al., 2009；Leveille, et al., 2011).随着地震勘探开发目标的日趋复杂化和精细化，地震数据偏移成像结果不仅要满足复杂区域构造成像的要求，还要具有保真性能反映地下反射面岩性变化的特征，使之能在岩性研究方面希望具有类似地震数据全波形反演的功效，以满足地震数据岩性处理解释，特别是复合型复杂油气藏和非常规油气藏地震数据岩性处理解释的需要.

当前的地震数据偏移成像方法起源于20世纪70年代(Claerbout，1971；1976)，其目的是利用波场延拓和成像条件把地表观测到的地震数据转化为反映地下反射面位置的信息.由于当前偏移成像方法中的成像条件(公式)的数学物理意义不明确(也不满足地震波场传播的控制方程)，因此得到的偏移成像结果与反射面岩性变化(反射系数)的关系含糊不清.陈生昌和周华敏(2016)在高频近似条件下，通过引入反射面反射率的定义，推导出基于反射率的反射地震波传播方程，即反射波动方程，然后以此为基础，利用线性反演方法对地震数据偏移进行重新推导，构建了基于反射波动方程的地震数据偏移成像方法，得到的偏移成像结果不仅数学物理意义明确、与反射面岩性变化关系清晰，而且相对于当前偏移成像方法的结果具有准确的相位和更高的分辨率.

当前的地震数据最小二乘偏移是基于线性反演理论对当前地震数据偏移成像方法的发展(Tarantola，1984)，其目的是通过消除地震数据观测系统、地下介质速度变化、地震波几何扩散以及偏移成像公式等对偏移成像结果产生的不利影响，得到一个具有高分辨、高保真、高信噪比的偏移成像结果(Bleistein et al., 2000；Yao and Jakubowicz, 2012a；2012b).当前的最小二乘偏移以一阶Born近似的散射波传播方程为基础，利用迭代反演的方式求取最小二乘偏移成像结果(Tang and Biondi, 2009；Ren et al., 2011；Zhang et al., 2014；黄建平等，2013；周华敏等，2014)，得到的是速度的相对扰动变化，而不是对速度异常体(层)边界(反射面)上物性变化(反射系数或反射率)的成像.我们认为散射波方程仅适合地下介质速度的局部小扰动体产生的散射波，而地下介质中大尺寸(相对于地震波主波长)的速度异常体(层)所产生的地震反射波用散射波方程进行描述是不合适的.此外，在当前的最小二乘偏移求解的迭代反演实现中，把入射波传播、散射波传播综合在一个线性传播算子中，并以该线性传播算子为基础构建具体迭代反演计算公式，没有区分入射波传播算子与散射波传播算子作用的不同.从地表激发、观测的地震波场的传播过程可知，观测波场仅是通过散(反)射波传播算子与地下散射体(反射面)建立联系，因此我们认为在最小二乘偏移的反演中也就仅需要考虑散(反)射波传播算子的作用，而不需要考虑入射波传播算子的作用.

根据陈生昌和周华敏(2016)有关地震数据反射波偏移成像方法的工作，本文在基于反射率反射波动方程的基础上，利用迭代方式求解线性反演问题的方法求取有关地下反射率的最小二乘偏移成像结果，建立基于反射波动方程的地震数据最小二乘偏移方法.在具体的每次迭代中，我们首先利用反源问题求解中的反向传播方法(Devaney, 2012)反演出由反射率与入射波场相互作用而产生的非齐次反射波动方程的源项的近似解，然后再结合地下入射波场得到地下反射率偏移成像结果的近似解.最后把提出的最小二乘偏移方法应用于理论模型的合成地震数据，检验其正确性和有效性.本文所提出的最小二乘偏移方法是对陈生昌和周华敏(2016)所提出的地震数据反射波偏移成像方法的发展.

1 方法原理 1.1 地震波场的传播与反射

在常规地面反射地震勘探中，地震波场的传播过程为：震源在地面激发地震波场，波场向下传播，在地下反射面处产生反射波场，反射波场向上传播到地面被检波器记录.

假设地下速度模型为v(x)，地震波场满足下述的标量波动方程：

(1)

式中，∇2为Laplace算子，有；xs为震源点位置；f(t)为震源时间函数；u(x, t; xs)表示地下地震波场，它包括直达波场、散射波场和反射波场等.方程(1)是一个通用的波动方程，它能描述速度模型v(x)空间变化引起的各种波现象，如直达波、散射波、反射波、透射波等，是进行波场数值模拟和全波形反演的基础方程，但在该方程中地震波场的传播和反射是耦合在一起的，不利于针对地震反射数据的偏移成像.

为了显式地表示反射波的传播和反射，我们把速度模型v(x)表示为下述形式，即

(2)

式中，vm(x)为一个光滑的速度模型，在该速度模型下，波动方程(1)仅存在直达波场，而不产生反射波场，因此该速度模型也称为宏观速度模型；r(x)为基于高频近似定义的地下反射面的反射率分布(陈生昌和周华敏，2016)；⊕表示一种概念相加的运算，因为vm(x)与r(x)有不同的物理量纲.

由速度模型表达式(2)，仅含直达波场而不含反射波场的入射波场所满足的波动方程为

(3)

式中ui(x, t; xs)表示入射波场.

根据陈生昌和周华敏(2016)，对应速度模型表达式(2)的一次反射波场所满足的波动方程为

(4)

式中，ur(x, t; xs)表示一次反射波场；右端项为产生反射波场的虚源，是由反射率r(x)与地下入射波场的时间一阶导数相互作用而产生的.

方程(3)、(4)中的波动算子仅与宏观速度模型vm(x)有关，因此波场(入射波场和反射波场)的走时是由vm(x)决定.由于vm(x)的光滑性，方程(3)、(4)中的波动算子仅是传播波场而不会产生反射波场，因此在反射波动方程(4)中，波场的传播和反射得到了显式的解耦，有利于针对反射波场的反演和成像.

1.2 地震数据的反射波动方程偏移

方程(3)、(4)分别是地震数据偏移中用来描述震源波场正向传播和记录波场反向传播的基础方程.根据陈生昌和周华敏(2016)所提出的偏移成像方法，可得到下述利用地震数据对反射率r(x)进行成像的反射波动方程偏移计算公式：

1) 地下入射波场计算：

(5)

2) 地下反射波场的伴随(逆时)重建(即地下反射源重建)：

(6)

3) 对地下反射率r(x)的近似估计(成像)：

(7)

式(6)中的urobs(xg, t; xs)为观测的共炮道集反射数据.

由于上面的地下反射波场重建是利用伴随(逆时)传播实现的，而不是采用逆(反演)传播实现的，因此计算得到的地下反射波场会受到炮点-检波点的分布、地下速度模型以及波场的几何扩散效应等的影响，不仅使地下反射波场的幅值不能得到准确的重建，还会降低地下反射波场的空间分辨率，进而影响地下反射率成像结果的保幅性和分辨率.另外，在反射率成像计算中，为了计算的稳定性，我们把本应为反射波场ur(x, t; xs)除以∂ui(x, t; xs)/∂t的运算改写为如公式(7)中的形式，也使计算得到的反射率结果不仅缺乏保真性，而且分辨率降低(因为公式(7)不但没有消除数据中的震源子波效应，反而还增加了震源子波效应).

上述的反射波动方程偏移成像计算公式(5)、(6)、(7)与现有的基于Claerbout偏移成像原理所建立的逆时偏移成像计算公式的区别在于式(7).本文的公式(7)是基于反射波动方程和线性反演方法推导出来的，而不是根据Claerbout的“时间一致性成像原理”写出.当然反射率成像公式(7)也是满足Claerbout的“时间一致性成像原理”的，即在公式(7)中ur(x, t; xs)和∂ui(x, t; xs)/∂t的走时是相同的.

1.3 地震数据的反射波动方程最小二乘偏移

为了消除地下反射波场伴随重建式(6)中炮点-检波点的分布、地下速度模型、波场的几何扩散效应以及成像公式(7)对反射率成像的不利影响，以获得保真性好、分辨率高的r(x)成像结果，我们提出基于反射波动方程(4)的最小二乘偏移方法.

反射波动方程(4)所对应的Green函数gr(x, t; xs)，有

(8)

利用方程(8)中的Green函数，把方程(4)写为积分形式，有

(9)

式中，*代表时间褶积；s(x, t; xs)代表方程(4)的右端项，即产生反射波场的虚源，有

(10)

把方程(9)写为频率域形式，有

(11)

式中，为ur(x, t; xs)的Fourier变换；为gr(x, t; x′)的Fourier变换; xs)为s(x′, t; xs)的Fourier变换，有

(12)

其中是地下入射波场ui(x, t; xs)的Fourier变换.根据入射波控制方程(3)，有

(13)

式中，为方程(3)所对应的Green函数的Fourier变换；F(ω)为震源时间函数f(t)的Fourier变换.

方程(11)可视为利用地面反射波场; xs)求解地下反射源的线性积分方程，即反源问题，把它写为矩阵方程形式，有

(14)

式中，为表示频率域的地震一次反射数据的矩阵；Wr代表以Green函数为积分核的传播矩阵；为待求参数矩阵.根据线性反演理论，方程(14)的最小二乘解为

(15)

式中，Wr-g和Wra分别为传播矩阵Wr的最小二乘广义逆矩阵和伴随矩阵.式(15)为地下反射波虚源的最小二乘重建结果.

利用式(15)的反演结果，再结合式(13)表示的地下入射波场，由式(12)就可得到地下反射率r(x)的最小二乘反演估计，即

(16)

式中的为取实部运算.

式(16)的r(x)估计也就是地震反射数据的最小二乘偏移结果，该式虽然在理论上给出了地下反射率r(x)的结果，但存在一些数值计算方面的难题：一是式右端中的逆矩阵(WraWr)-1的大规模计算问题；二是式右端中的波场相除运算会出现计算的不稳定问题.

1.4 最小二乘偏移的迭代实现

为了避免式(15)中逆矩阵的大规模计算问题，我们采用Landweber迭代方法(Kirsch, 2011)求解矩阵方程(14)中的，相应地得到r(x)的迭代解.下述为求解地下反射率r(x)的迭代格式.

初始化，令：r0(x)=0，k=0；

迭代公式：

(17)

(18)

(19)

(20)

(21)

式(18)中的为第k次迭代的数据残差，且有.在式(20)中，α为步长因子，对于α可用不同方法求取，如简单的线性搜索法，也可用((WraWr)-1的近似值，如检波点方向的地震波单向照明强度.

式(21)分母项中的-iω为频率域中的导数运算，它与时间域的导数运算∂/∂t对应.根据波动方程(3)关于右端源项的线性性质，时间域的导数波场∂ui(x, t; xs)/∂t可利用下述波场方程进行计算，即

(22)

式中，有.

为了式(21)计算的稳定性，可把它写为

(23)

式中，为的复共轭；ε为一小的正数.

在最小二乘偏移的迭代过程中，由式(22)得到的地下入射波场的时间一阶导数uid(x, t; xs)是不变的.式(18)中的代表频率域中由反射率得到地震数据的反偏移运算.式(19)中所表示的运算，可通过对数据残差的伴随(逆时)传播来实现，如式(6)的运算.对于式(23)，如果取k=0，α=1，不考虑式中分母的作用，再转换到时间域，则它就退化为成像公式(7).

上述的反射波动方程最小二乘偏移方法与当前常用的基于一阶Born近似的散射波方程的最小二乘偏移方法是不同的，一是采用了不同的波动方程描述反射波的传播，即分别采用反射波动方程与散射波动方程描述反射波的传播，导致了有不同的偏移成像结果、偏移成像公式和反偏移计算公式；二是采用了不同的计算实现方式，本文方法把入射波传播、反射、反射波传播分开考虑，反演实现仅是针对反射波的传播而开展的，而当前的最小二乘偏移方法把入射波传播、散射波传播综合在一个线性传播算子中，因此其反演实现是基于该综合传播算子进行的，没有区分入射波传播与散射波传播作用的不同，如在实现中对于Hessian矩阵的计算.

2 数值试验

为了验证上述所提出的反射波动方程最小二乘偏移成像方法的正确性与有效性，我们应用两个模型的合成地震数据进行数值试验.

2.1 层状模型试验

为了验证本文提出的最小二乘偏移成像方法的正确性，我们用简单的层状模型进行试验，图 1为该模型的速度分布.利用波动方程(1)的有限差分方法进行合成地震数据模拟.在本试验总共模拟了300炮，炮点范围是1.5 km到4.5 km，炮间距10 m，为中间放炮两边接收，每炮301道，道间距10 m，时间采样间隔是1 ms，时间采样点数1700，震源为主频25 Hz的Ricker子波.图 2为用于偏移试验的光滑速度模型，它是通过对图 1的速度进行光滑处理得到的.

图 1 层状模型的速度分布 Fig. 1 The velocity distribution of layers model
图 2 用于偏移的光滑速度模型 Fig. 2 The smooth velocity model for migration

在试验中，我们首先应用1.2节提出的反射波动方程偏移方法对合成的炮道集地震数据进行偏移成像，图 3为利用反射波动方程偏移方法得到的偏移成像结果.对比图 1、3可看出，除了两端由于炮覆盖不足引起的偏移成像效果不好外，反射波动方程偏移方法在层状模型上取得了很好的偏移成像效果.图 4为应用本文提出的反射波动方程最小二乘偏移成像方法对层状模型的偏移成像结果，共进行了16次迭代.对比图 3、4，反射波动方程最小二乘偏移结果的分辨率相对于反射波动方程偏移结果有明显的提高，而且沿层面的振幅均衡性也好.通过该试验验证了本文所提出的反射波动方程最小二乘偏移方法的正确性.

图 3 反射波方程的偏移结果 Fig. 3 The result of reflection wave equation migration
图 4 16次迭代的反射波动方程最小二乘偏移结果 Fig. 4 The result of reflection wave equation least squares migration through 16 iterations
2.2 部分Sigsbee2a模型试验

为了进一步验证本文方法的有效性，我们截取Sigsbee2a速度模型的一部分进行试验，如图 5所示.同样利用波动方程(1)的有限差分方法进行合成地震数据模拟.在本试验总共模拟了300炮，炮点范围是0 km到6 km，炮间距20 m，为中间放炮两边接收，每炮601道，道间距10 m，时间采样间隔是0.5 ms，时间采样点数3000，震源为主频18 Hz的Ricker子波.图 6为用于偏移试验的光滑速度模型，它同样是通过对图 5的速度进行光滑处理得到的.

图 5 部分Sigsbee2a模型的速度分布 Fig. 5 The velocity distribution of part Sigsbee2a model
图 6 用于偏移的光滑速度模型 Fig. 6 The smooth velocity model for migration

在试验中，我们首先也应用1.2节提出的反射波动方程偏移方法对合成的炮道集地震数据进行偏移成像，图 7为利用反射波动方程偏移方法得到的偏移成像结果.由图 7可看出，反射波动方程偏移方法对该模型可获得很好的偏移成像结果.图 8为应用本文提出的反射波动方程最小二乘偏移成像方法对部分Sigsbee2a模型的偏移成像结果，共进行了35次迭代.对比图 7、8，反射波动方程最小二乘偏移结果的分辨率相对于反射波动方程偏移结果的分辨率有明显的提高，而且层位面与断层面的成像更清晰.与图 9所示的部分Sigsbee2a模型的垂向理论反射率对比，可看出图 8的最小二乘偏移结果的振幅保真性相对图 7的偏移结果也有明显的提高.通过该试验进一步验证了本文所提出的反射波动方程最小二乘偏移方法在复杂速度模型上的有效性.

图 7 反射波动方程的偏移结果 Fig. 7 The result of reflection wave equation migration
图 8 35次迭代的反射波动方程最小二乘偏移结果 Fig. 8 The result of reflection wave equation least squares migration through 35 iterations
图 9 图 5速度模型的垂向反射率图 Fig. 9 The vertical reflectivity of the velocity model in Fig. 5

本文所提出的基于反射波动方程的最小二乘偏移方法在实质上是一种获取反射面上反射率估计的线性迭代反演，而当前常见的基于散射波动方程Born近似的最小二乘偏移方法是一种获取地下速度扰动(或速度相对扰动)估计的线性迭代反演，即本文的最小二乘偏移方法得到的是地下反射层界面的物性参数，而当前的最小二乘偏移方法得到的是地下反射层的物性参数.因此，本文的最小二乘偏移方法与当前的最小二乘偏移方法虽然都属于利用迭代方法对地震数据波形线性反演问题的求解，但由于两者所基于的正演方程的不同和反演目标的不同，致使它们之间没有可比性.地震数据的最小二乘偏移是一种利用迭代方法进行求解的线性反演，而当前的地震数据全波形反演是一种利用迭代的线性化反演方法进行求解的非线性反演(Virieux and Operto, 2009；陈生昌和陈国新，2016)，所以本文的最小二乘偏移方法与全波形反演方法也不存在可比性.

3 结论

本文所提出的最小二乘偏移方法是对陈生昌和周华敏(2016)所提出的地震反射波数据偏移成像方法的发展，属于对以反射波动方程为基础的地震数据波形线性反演问题的迭代求解.在最小二乘偏移的具体计算中，考虑了入射波传播算子与反射波传播算子对偏移成像作用的不同，最小二乘偏移中反演的迭代实现仅是针对反射波传播进行的.本文的反射波动方程最小二乘偏移方法与当前基于一阶Born近似的散射波方程的最小二乘偏移方法的主要不同有：1)采用了不同的波动方程描述反射波的传播，导致了有不同的偏移成像结果、偏移成像公式和反偏移计算公式；2)与最小二乘偏移所对应的线性算子不同，本文的最小二乘偏移所对应的线性算子仅由反射波传播算子组成，而当前的最小二乘偏移所对应的线性算子由入射波传播算子和反射波传播算子结合组成，因此构建出了不同的具体实现方式.通过理论模型合成地震数据的数值试验验证了本文的最小二乘偏移方法的正确性和有效性.

参考文献
Bleistein N, Cohen J K, Stockwell J W Jr. 2000. Mathematics of Multidimensional Seismic Inversion, Migration, and Inversion. New York: Springer. http://adsabs.harvard.edu/cgi-bin/nph-data_query?bibcode=2006mmsi.book.....B&db_key=PHY&link_type=ABSTRACT&high=28938
Chen S C, Chen G X. 2016. Full waveform inversion of the second-order time integral wavefield. Chinese Journal of Geophysics (in Chinese), 59(10): 3765-3776. DOI:10.6038/cjg20161021
Chen S C, Zhou H M. 2016. Re-exploration to migration of seismic data. Chinese Journal of Geophysics (in Chinese), 59(2): 643-654. DOI:10.6038/cjg20160221
Claerbout J F. 1971. Toward a unified theory of reflector mapping. Geophysics, 36(3): 467-481. DOI:10.1190/1.1440185
Claerbout J F. 1976. Fundamentals of Geophysical Data Processing. New York: McGraw-Hill Book Co. Inc.
Devaney A J. 2012. Mathematical Foundations of Imaging, Tomography and Wavefield Inversion. Cambridge: Cambridge University Press.
Etgen J, Gray S H, Zhang Y. 2009. An overview of depth imaging in exploration geophysics. Geophysics, 74(6): WCA5-17. DOI:10.1190/1.3223188
Gray S, Etgen J, Dellinger J, et al. 2001. Seismic migration problems and solutions. Geophysics, 66(5): 1622-1640. DOI:10.1190/1.1487107
Huang J P, Li Z C, Kong X, et al. 2013. A study of least squares migration imaging method for fractured-type carbonate reservoir. Chinese Journal of Geophysics (in Chinese), 56(5): 1716-1725.
Kirsch A. 2011. An Introduction to the Mathematical Theory of Inverse Problems. 2nd ed. New York: Springer.
Leveille J P, Jones I F, Zhou Z Z, et al. 2011. Subsalt imaging for exploration, production, and development:A review. Geophysics, 76(5): WB3-WB20. DOI:10.1190/geo2011-0156.1
Ren H R, Wu R S, Wang H Z. 2011. Wave equation least square imaging using the local angular Hessian for amplitude correction. Geophysical Prospecting, 59(4): 651-661. DOI:10.1111/gpr.2011.59.issue-4
Tang Y X, Biondi B. 2009. Least-squares migration/inversion of blended data. //79th Ann. Internat. Mtg., Soc. Expl. Geophys. Expanded Abstracts. http://www.researchgate.net/publication/269067982_Leastsquares_migrationinversion_of_blended_data
Tarantola A. 1984. Linearized inversion of seismic reflection data. Geophysical Prospecting(6): 998-1015.
Virieux J, Operto S. 2009. An overview of full-waveform inversion in exploration geophysics. Geophysics, 74(6): WCC127. DOI:10.1190/1.3237087
Yao G, Jakubowicz H. 2012a. Least-squares reverse time migration. //82nd Ann. Internat. Mtg., Soc. Expl. Geophys. Expanded Abstracts. http://ethos.bl.uk/OrderDetails.do?uin=uk.bl.ethos.616810
Yao G, Jakubowicz H. 2012b. Non-linear least-squares reverse time migration. 82nd//Ann. Internat. Mtg., Soc. Expl. Geophys. Expanded Abstracts
Yilmaz O. 2001. Seismic Data Analysis. Oklahoma: Society Of Exploration Geophysicists.
Zhang Y, Duan L, Xie Y. 2014. A stable and practical implementation of least-squares reverse time migration. Geophysics, 80(1): V23-V31.
Zhou H M, Chen S C, Ren H R, et al. 2014. One-way wave equation least squares migration based on illumination compensation. Chinese Journal of Geophysics (in Chinese), 57(8): 2644-2655. DOI:10.6038/cjg20140823
陈生昌, 陈国新. 2016. 时间二阶积分波场的全波形反演. 地球物理学报, 59(10): 3765–3776. DOI:10.6038/cjg20161021
陈生昌, 周华敏. 2016. 再论地震数据偏移成像. 地球物理学报, 59(2): 643–654. DOI:10.6038/cjg20160221
黄建平, 李振春, 孔雪, 等. 2013. 碳酸盐岩裂缝型储层最小二乘偏移成像方法研究. 地球物理学报, 56(5): 1716–1725. DOI:10.6038/cjg20130529
周华敏, 陈生昌, 任浩然, 等. 2014. 基于照明补偿的单程波最小二乘偏移. 地球物理学报, 57(8): 2644–2655. DOI:10.6038/cjg20140823