地球物理学报  2020, Vol. 63 Issue (9): 3442-3451   PDF    
一种新的针对反射波的逆时偏移真振幅成像条件
刘文卿1, 周阳2, 雍学善1, 王小卫1, 王孝1, 胡书华1, 田彦灿1     
1. 中国石油勘探开发研究院西北分院, 兰州 730020;
2. 成都理工大学地球物理学院, 成都 610059
摘要:真振幅成像是一种代表性的定量估计模型参数扰动高波数部分的地震波成像方法.经典的真振幅成像方法在高频近似和理想照明假设条件下求取显式对角Hessian逆矩阵作为偏移振幅加权算子,用以校正波传播过程中的几何扩散效应,得到模型参数扰动的带限估计.真振幅保真成像方法在利用逆时偏移(RTM)框架实现时会产生低波数噪声,影响对高波数参数估计的精度.本文给出了一种新的基于RTM框架的真振幅保真成像条件,该成像条件针对反射波数据,在高频近似下散射模式对应正问题及Bayes反问题框架下导出.与传统基于高频渐进反演的波动方程成像方法类似,利用本文提出RTM成像条件能够保证计算结果与高频近似下反演结果的一致性.同时,利用本文提出RTM真振幅成像条件能够在成像过程中自动保真的消除传统真振幅RTM算法中存在低波数噪声,模型数值实验结果验证了本文方法的正确性和有效性.
关键词: 真振幅      反射波      反演成像      双向波传播算子      低波数噪声      成像条件      逆时偏移     
A new true-amplitude imaging condition of reverse time migration for reflected seismic data
LIU WenQing1, ZHOU Yang2, YONG XueShan1, WANG XiaoWei1, WANG Xiao1, HU ShuHua1, TIAN YanCan1     
1. Research Institute of Petroleum Exploration & Development-Northwest, Petrochina, Lanzhou 730020, China;
2. College of Geophysics, Chengdu University of Technology, Chengdu 610059, China
Abstract: True-amplitude imaging is one of the typical seismic imaging method for quantitatively imaging high-wavenumber model perturbation. Under high-frequency and perfect illumination assumptions, the classical true-amplitude imaging method adopt an explicit inverse diagonal Hessian for removing geometric spreading during wave propagation, aiming at obtaining band-limited model perturbation. When implementing true-amplitude imaging through a RTM (Reverse Time Migration) process, the unwanted low-wavenumber artifacts are produced, which will contaminate the imaging accuracy for true seismic reflectors. In this paper, we propose new inversion imaging condition based on RTM process. The present RTM imaging condition aims at imaging reflected seismic data and derived with linearized surface scattering model, together with Bayes inverse problem framework. As the conventional asymptotic inversion based wave-equation imaging methods, the output of new RTM imaging condition is consistent with the asymptotic inversion in high frequency sense. In the meantime, the proposed imaging condition automatically suppress the RTM artifacts in a way that the phase and amplitude of the target reflector image is maintained. Numerical examples on synthetic data verify effectiveness of the present RTM imaging condition.
Keywords: True-amplitude    Reflected seismic data    Inversion imaging    Two-way wave propagation operator    Low wavenumber noise    Imaging condition    Reverse time migration    
0 引言

地震波保真成像定位于给出模型参数高波数扰动部分的定量估计,是进行精确油气藏描述的重要技术.最小二乘偏移是一种具有代表性的保真成像技术,通过匹配模拟数据与观测数据的振幅来估计参数的高波数扰动.然而,最小二乘偏移的理想实现需将Hessian逆算子作用到常规成像结果上.面对实际资料,由于Hessian矩阵规模巨大,通常采用迭代方式求取Hessian的逆算子,这使得最小二乘偏移成像成为一种计算昂贵的保真成像方法.另外,Hessian逆算子的估计精度很容易受到叠前数据观测不规则、正算子对物理波现象描述不精确以及地震子波未知等多种因素影响,从而出现迭代收敛慢甚至不收敛情况,阻碍了最小二乘偏移成像在大规模实际资料中的实用化进程.

在勘探地震学中,真振幅成像是另一种较为常见的保真成像方法.在高频近似和理想照明假设条件下将Hessian矩阵对角化,得到显式的Hessian逆矩阵作为偏移振幅加权算子.物理上,此时的偏移振幅加权算子校正了波传播过程中的几何扩散影响.与迭代实现的最小二乘偏移相比,真振幅成像由于显式利用了近似的Hessian逆矩阵,能在保证数值计算高效稳健的前提下提高成像精度,获得带限参数扰动的定量估计,在实际生产中获得了较多的应用(Gray, 1997Brandsberg-Dahl et al., 2003Zhang et al., 2014).值得注意的是,真振幅成像公式同样可以作为预条件梯度项融入到最小二乘迭代求解框架下实施,从而进一步提高参数反演成像的精度(Jin et al., 1992; Lambare et al., 1992; Hou and Symes, 2015; Li and Chauris, 2018).根据不同的模型空间定义,真振幅成像可以分为高频近似下的Born近似(ray+Born近似)以及高频近似下的Kirchhoff近似(ray+Kirchhoff近似)两类具体的实现方式(Beydoun and Jin, 1994).其中ray+Born近似下真振幅成像的模型空间假设为短波长的体状散射体,采用体积分描述散射场和模型参数扰动之间的线性关系(Beylkin, 1985).对于ray+Kirchhoff近似下的真振幅成像方法,其假设的模型空间为代表地层分界面的平滑反射界面,采用面积分描述散射场和模型参数扰动之间的线性关系(Bleistein, 1987).除了基于射线类传播算子实现(Miller et al., 1987Gray and Bleistein, 2009),真振幅成像方法也可以基于波动方程类传播算子进行实现.通过将真振幅成像公式中的Green函数项利用波场外推方式进行计算,可以在高频近似意义下使得基于波场成像结果与最小二乘框架下的射线类真振幅成像结果保持一致.Esmersoy和Miller(1989)较早讨论了射线类真振幅成像公式与波动方程类成像方法的联系,指出在高频近似下可以建立两者显式关系.Bleistein等(2005)系统讨论了利用波动方程算子实现ray+Kirchhoff意义下真振幅成像公式的方法.基于类似的思想,Zhang等(2007)Kiyashchenko等(2007)给出了利用单程波算子实现ray+Kirchhoff意义下真振幅成像公式具体实施方案.Zhang和Sun (2009)Xu等(2011)进一步给出了利用双程波传播算子实现ray+Kirchhoff意义下真振幅成像公式的方法,即ray+Kirchhoff近似下的真振幅逆时偏移方法(RTM).Stolk等(2009)Zhang等(2014)讨论了ray+Born近似意义下真振幅成像公式的双程波算子实现, 即ray+Born近似下的真振幅RTM方法.

对于经典意义下的真振幅成像,是在低波数背景参数已知情况下,计算得到高波数的参数扰动,波场成分为震源端和检波器端波矢量夹角较小的部分.由于双向波算子可以模拟空间上全方向传播的波场,因此在利用其实现真振幅成像过程中,会存在源检两端波场对应波矢量夹角较大的情形,对这部分波场进行相关操作会产生低波数能量.这些低波数能量可以用来更新背景速度(Díaz and Sava, 2016),而对于真振幅成像,这些低波数能量会产生成像噪声,需要在成像过程中进行消除(Liu et al., 2011).基于对成像结果的波数域高通滤波或对波场显式角度域分解是当前应用最广泛的压制低波数噪声的方法(Fletcher et al., 2005Guitton et al., 2006Wu and Xie, 2002Zhang and Sun, 2009).压制低波数噪声方法需在成像过程外增加额外计算步骤,另外这些方法在低波数噪声去除过程中容易破坏目标反射层成像结果原始的波形信息,需要进一步设计校正来消除影响,在一定程度上增加了成像方法计算的复杂度.与以上消除低波数噪声具体做法不同的是,Stolk等(2009)提出的基于ray+Born近似的真振幅RTM方法,通过引入波场的隐式方向分解,在波场外推成像过程中自动完成低波数噪声的消除.Whitmore和Crawley(2012)Brandsberg-Dahl等(2013)Stolk等(2009)给出的成像方法应用到实际数据,取得了不错的效果.

本文给出了一种新的真振幅逆时偏移(RTM)成像公式,该成像条件针对反射波数据,能够在保持目标反射层成像结果原始波形信息前提下自动消除RTM成像过程中的低波数噪声.同时,由于该成像条件严格在高频近似下的反问题框架下进行推导,其理论输出结果和射线类高频渐进反演结果具有较好的一致性,是一种定量化的偏移成像方法.我们用基于合成数据的数值实验验证本文提出成像条件的有效性.

1 高频近似下的反演成像及其波动方程实现

在层状反射地层介质假设的前提下,利用基于Kirchhoff近似的线性化散射场比利用基于Born近似的线性化散射场能够更好描述实际接收的反射波数据(Bleistein, 1987).因此,在这种情况下基于ray+Kirchhoff近似下的真振幅成像比ray+Born近似下的真振幅成像在成像精度上更有优势(Beydoun and Jin, 1994).由于本文旨在发展新的针对反射数据的真振幅逆时偏移成像条件,我们首先讨论ray+Kirchhoff近似下真振幅成像公式及其基于波动方程的实现.

1.1 Ray+Kirchhoff近似下的真振幅成像公式

在ray+Kirchhoff近似下,在炮点xs处激发,被检波点xr接收到的二维频率域线性化散射场uS(xr, xs, ω)可以表示为(Xu et al., 2001):

(1)

其中R(x, θr(xr, x, xs))为角度反射系数,θr(xr, x, xs)为反射角,ln分别为震源端和检波点端的射线标号,N(x, xs)表示从震源xs激发到成像点x的射线路径数目,N(xr, x)表示从成像点x到检波点xr接收到的射线路径数目.GnlRK(xr, x, xs, ω)为定义模型空间和数据空间扰动关系的核函数,其具体形式为

(2)

这里,|qnl(xr, x, xs)|为偏移倾角矢量的模(Operto et al., 2000),其定义为

(3)

c0(x)为x处的背景速度.(2)式中的Anl(xr, x, xs)表示从炮点xs出发到地下任一点x然后被xr处接收到射线振幅的改变.相位项Tnl(xr, x, xs)的表达式为

(4)

其中τ(xr, x, xs)为炮点xs出发到地下任一点x然后被xr处接收到射线的总走时.Knl(xr, x, xs)描述射线从炮点xs出射,被检波点xr接收到所经历的一阶焦散点的总个数,在每一个一阶焦散点处,高频近似下的波场发生π/2的相移(Červený,2001).

在Bayes反演框架及高频近似假设下,利用稳相近似及选取合适的数据协方差矩阵,我们可以将Hessian逆矩阵对角化,从而构造散射场计算正问题(1)对应的反问题的解(Jin et al., 1992Ten Kroode et al., 1994),即所谓ray+Kirchhoff近似下的真振幅成像公式,二维情况下其具体形式为(Xu et al., 2001):

(5)

其中θminθmax分别为当前观测系统对地下成像点能达到的最小和最大散射角度值.振幅权重εnlRK的形式为

(6)

J(xs, xr, ω)→(kx, kz, θs)为笛卡尔坐标向相空间坐标转化的行列式,其表达式为(Thierry et al., 1998):

(7)

式中βxsβxr分别为源检两端的起飞角.(5)式中的δ(θrθr)为角度采样函数,使得积分仅对在成像点x处反射角为θr的炮数据进行(Bleistein et al., 2005).

1.2 Ray+Kirchhoff近似下真振幅成像公式的波动方程实现

真振幅成像公式可以通过如下基本步骤利用波动方程算子进行实现.首先将震源两端地震波场在半无限空间中的传播近似为Rayleigh-Sommerfeld外边值问题(Devaney,2012),此时震源两端波场可以利用Rayleigh第二类积分进行表示(Wapenaar and Berkhout, 1989),通过将积分表达的震源两端波场做进一步高频近似,我们可以建立起外推波场与真振幅成像公式(5)的联系,得到利用求解的波动方程外推波场表征的真振幅成像公式,即所谓波动方程真振幅成像条件,上述具体推导可参见Xu等(2011).二维情况下ray+Kirchhoff真振幅成像公式(5)对应的波动方程成像条件为(Zhang et al., 2007):

(8)

其中,PF(x, xs, ω),PB(x, xs, ω)分别为频率域的正传及反传播波场.成像条件式(8)的时间域表征为

(9)

这里H代表Hilbert变换.注意到在一些文献中,给出的反射系数定义为R′(x, θr)=R(x, θr)/|q|(Joncour et al., 2011Xu et al., 2011),则(9)式可以写为

(10)

(10) 式为Zhang等(2007)给出的二维情况下的波动方程真振幅成像条件.(9)式中|q|为水平层状介质假设下的子波拉伸校正算子(Tygel et al., 1994),对保真角度反射系数的获取是有必要的(Kühl and Sacchi, 2003).因此在本文以(9)式作为后续进一步推导新的真振幅RTM成像公式的基础.

2 针对反射波的新的RTM真振幅成像条件

在利用RTM框架实现公式(9)时会出现由于震源两端的波路径重合而产生的低波数能量,对于估计模型高波数参数,这些低波数能量成分是成像噪声应当去除(Liu et al., 2011).为了消除这些低波数噪声,受Stolk等(2009)工作的启发,我们引入角度域滤波算子F,通过将F作用到成像条件上来去除低波数噪声.波数域F的表达式为

(11)

其中θs为震源端和检波点端波场对应波矢量之间的夹角,一般将其称之为散射角(Forgues and Lambaré,1997).从(11)式可以看出随着散射角度的增大,滤波算子权重减小,这是(11)式能压制低波数噪声的基本原理.将(9)式变换到频率波数域并代入(11)式,可以得到带角度滤波的波场相关表达式:

(12)

其中:

(13)

对于(13)式中的,利用频散关系k=ω/c,则有:

(14)

对于(13)式中的,可写成:

(15)

利用向量内积的定义,(15)式可进一步写为:

(16)

将(14)与(16)式利用傅里叶反变换到时空域,并利用时间与空间导数的定义,带角度滤波的波场相关时空域表达式为

(17)

其中

(18)

从公式(18)可看出,角度域滤波算子式(11)的最终实施无需显式进行波场角度分解,而是转化为对波场时间导数和空间梯度的代数运算,且不涉及波场的除法操作,相较于基于波场显式方向分解压制低波数噪声的方法,公式(18)数值计算上更为稳健.另外,由于波场时间导数和空间梯度是逐点求解,相较于波数域滤波方法压制低波数噪声,利用公式(18)能够更好地适应模型的空变(Whitmore and Crawley, 2012).利用(17)式替换(9)式中的波场相关运算,可以得到带角度域滤波的二维ray+Kirchhoff近似下真振幅RTM成像条件,但这并非本文给出成像条件的最终结果.这是因为(11)式中算子|k|2是非单位化算子,直接将其作用到波场上会改变波场的原始波形信息(Zhou and Wang, 2016),为了在消除低波数噪声的同时保持原始成像剖面中有效反射层的波形信息,需要对波场做进一步校正.利用频散关系,针对(11)式我们定义波形校正算子Fc

(19)

将(19)式变换到时间域同(18)式一起代入(9)式最终得到新的RTM真振幅成像条件:

(20)

其中I0I1形式上与(18)式相同,区别在于I0I1中的波场为作用频率域算子1/ω2滤波后的波场.为了提高计算效率,可将对波场的滤波操作转化为对震源子波的频率域滤波或对时间的积分(Xu et al., 2011).公式(20)即为本文给出的新的真振幅RTM成像条件.相较于采用类似(18)式角度域滤波算子的已发表工作(Stolk et al., 2009Whitmore and Crawley, 2012Rocha et al., 2016),本文提出的新成像条件严格在ray+Kirchhoff反问题框架下给出,因此更适用于层状反射地层介质假设下的反射波数据情形,其理论输出结果物理意义明确,即为校正几何扩散后的带限角度反射系数的估计.

另外,从公式(20)可以看出,直接利用此成像条件计算的模型空间为角度域反射系数,为了在不输出角度反射系数的情况下得到有物理意义的叠加结果,可以利用常密度声学近似下的AVA近似关系式(Shuey, 1985)将利用公式(20)得到的不同角度反射系数结果校正到零角度进行叠加成像(Kiyashchenko et al., 2007),最终输出零角度反射系数成像剖面,本文后续的数值实验采用此策略产生最终的成像结果.

3 数值实验

我们利用层状模型验证本文提出的成像条件的正确性和有效性.模型第一层速度为3000 m·s-1,第二层速度为4000 m·s-1.模型的xz方向范围均为0~4000 m, 网格大小两个方向上均为10 m.我们利用真实模型正演得到的数据作为输入炮记录,利用Gussian函数平滑真实速度后的模型作为成像算法采用的背景速度.图 1a为利用不带角度滤波的ray+Kirchhoff近似下真振幅成像条件(9)式成像结果,从成像结果可看出,目标反射层的成像为带限的脉冲函数,这与ray+Kirchhoff真振幅成像的期望输出是一致的(Bleistein et al., 2001).同时,在目标反射层上方存在可以预期的低波数噪声.图 1b为利用本文提出的成像条件(20)式成像结果,可以看出,低波数噪声被消除了,同时成像结果的振幅量级与图 1a保持了一致.图 1c为利用本文提出的成像条件滤出的低波数噪声.图 1d为直接利用不带振幅校正的角度域滤波(18)式的成像结果,虽然低波数噪声被消除了,但是由于角度域滤波算子(11)式的非单位化效应未被消除,成像结果的振幅量级相较于原始成像结果图 1a发生了明显变化.

图 1 单层模型不同成像条件成像结果的对比 (a)利用不带角度滤波的ray+Kirchhoff近似下真振幅成像条件(9)式成像结果;(b)利用本文给出成像条件(20)式的成像结果;(c)利用本文给出成像条件(20)式滤出的低波数噪声;(d)利用不带振幅校正的角度域滤波(18)式的成像结果. Fig. 1 Imaging results for different imaging conditions of a single layer model (a) Results with true-amplitude imaging condition (9) using ray+Kirchhoff approximation without scattering-angle filtering; (b) Results using proposed imaging conditions (20); (c) Imaged low-wavenumber component noise using (20); (d) Imaging results using scattering-angle filtering (18) without amplitude correction.

图 2为不同成像条件的抽道振幅对比,通过对比可看出利用本文提出成像条件(20)式得到的地震道振幅在反射层位置与利用(9)式得到的未进行角度滤波结果吻合很好,而滤出的低波数噪声则和利用(9)式得到结果中低波数部分在振幅上有较好的一致性.这更加有力地证明了利用本文成像条件能够在自动消除低波数噪声的同时保持目标反射层振幅保真和相位的稳定,这对后续进一步的处理解释是十分有利的.

图 2 单层模型成像条件成像结果抽道对比 (a)利用不带角度滤波的真振幅成像条件(9)式成像结果(蓝实线)与利用本文成像条件(20)式成像结果(红虚线)对比;(b)利用不带角度滤波的真振幅成像条件(9)结果(蓝实线)与利用本文成像条件滤出的低波数噪声(红虚线)的对比. Fig. 2 Trace comparison of images with different imaging conditions of a single layer model (a) Using true-amplitude imaging conditions without scattering-angle filtering (9) (blue solid line) versus imaging results using imaging condition (20) (red dashed line) presented in this paper; (b) Comparison of filtered low wavenumber noise using true-amplitude imaging condition without angle filtering (9) results (blue solid line) and using the proposed imaging condition (20) in this paper (red dashed line).

为了进一步证明本文提出的成像条件的优势,我们利用四层模型进行不同成像条件成像对比.图 3a为不带角度滤波的ray+Kirchhoff近似下真振幅成像条件(9)式成像结果.图 3b为利用本文提出的成像条件(20)式成像结果.图 3c为利用本文提出的成像条件(20)式滤出的低波数噪声.图 3d为直接利用不带振幅校正的角度域滤波(18)式的成像结果.与单层模型数值实验结论类似,利用本文提出成像条件(20)式和(18)式均能去除成像条件(9)式成像结果中的低波数噪声.由于本文提出成像条件校正了角度域滤波算子(11)式非单位化效应,因此在去除低波数噪声的同时保持了反射层振幅量级与原始成像结果图 3a的一致性.对比图 3b3d,我们还注意到角度滤波算子(11)式的非单位化效应不仅会改变目标反射层的绝对振幅值,同时也会改变从浅到深不同反射层的相对振幅值,这会直接影响成像结果的保真度.因此在成像条件方面,本文提出成像条件是基于双向波传播算子的反射波波动方程保真成像算法的一个最佳选择.

图 3 四层模型不同成像条件成像结果的对比 (a)利用不带角度滤波的ray+Kirchhoff近似真振幅成像条件(9)式成像结果;(b)利用本文给出成像条件(20)式的成像结果;(c)利用本文给出成像条件(20)式滤出的低波数噪声;(d)利用不带振幅校正的角度域滤波(18)式的成像结果. Fig. 3 Imaging results for different imaging conditions of a four-layer model (a) Results with true-amplitude imaging condition (9) using ray+Kirchhoff approximation without scattering-angle filtering; (b) Results using proposed imaging conditions (20); (c) Imaged low-wavenumber component noise using (20); (d) Imaging results using scattering-angle filtering (18) without amplitude correction.

图 4a为利用不带角度滤波的ray+Kirchhoff近似下真振幅成像条件(9)与利用本文提出的成像条件(20)式成像结果抽道振幅的对比,可以看出图 4a中利用成像条件(9)式得到成像结果(蓝色实线)中的低波数噪声,利用本文提出的成像条件得到成像结果(红色实线)显然不包含低波数噪声,且目标反射层振幅能级和利用(9)式得到结果具有很好的一致性.同时我们注意到图 4a中利用(9)式得到的前两层反射层成像结果振幅要大于利用本文提出成像条件得到的结果,这是因为利用(9)式得到的第一、二层成像结果振幅分别叠加上了来自于第二、第三层反射层产生的低波数噪声的缘故.从图 4a还可以看出,利用(9)式和本文提出成像条件得到的第三层反射层成像结果振幅值由于没有后续反射层产生的低波数噪声干扰完全相互吻合,这也从另一个侧面佐证了前面关于振幅匹配的分析.图 4b为利用不带角度滤波的真振幅成像公式(9)与利用本文提出成像条件滤出的低波数噪声的对比,可以看出利用本文提出成像条件滤出的低波数能量与利用(9)式成像结果中低波数能量部分能够很好地匹配.

图 4 四层模型不同成像条件成像结果抽道对比 (a)利用不带角度滤波的真振幅成像条件(9)式成像结果(蓝色实线)与利用本文提出的成像条件(20)式成像结果(红色实线)的对比;(b)利用不带角度滤波的真振幅成像条件(9)成像结果(蓝色实线)与利用本文提出成像条件滤出的低波数噪声(红色实线)的对比. Fig. 4 Trace comparison of images with different imaging conditions for a four-layer model (a) Using true-amplitude imaging conditions without angle filtering (9) (blue solid line) versus imaging results using imaging condition (20) (red solid line) presented in this paper; (b) Comparison of filtered low wavenumber noise using true-amplitude imaging condition without angle filtering (9) results (blue solid line) and using the proposed imaging condition (20) in this paper (red solid line).

我们利用Marmousi模型进行试算,试算结果如图 5所示.图 5a为利用不带角度滤波的ray+Kirchhoff近似下真振幅成像条件(9)成像结果.图 5b为利用本文提出方法的成像结果.图 5c为直接利用不带振幅校正的角度域滤波(18)式的成像结果.通过对比三种不同成像条件得到的成像结果,我们可以清晰地看出利用本文方法成像结果(图 5b)在消除低波数噪声同时保持了反射层振幅量级与原始成像结果图 5a的一致性,同时不同深度反射层相对振幅关系相较于直接利用不带振幅校正的角度域滤波(18)式成像结果更加合理.图 5d5e、5f分别为利用不带角度滤波的真振幅成像条件式(9)成像结果(蓝色实线)与利用本文提出的成像条件成像结果(红色实线)第40、120、450道的抽道振幅对比结果.可以看出两者的振幅量级是一致的.与三层模型成像结果类似,利用式(9)成像结果中浅层部分反射层成像结果由于和后续反射层产生低波数能量相互叠加,其振幅值相较于利用本文提出成像条件得到的中浅层反射层成像结果振幅值要更大一些,而深层反射层成像结果受后续反射层产生低波数能量影响较小,因此其振幅值与利用本文提出成像条件得到的成像结果振幅值能够很好的匹配.Marmousi模型数值试算结果表明了本文提出成像条件对复杂模型的适应性和正确性.

图 5 Marmousi模型不同成像条件成像及其抽道对比结果 (a)利用不带角度滤波的ray+Kirchhoff近似下真振幅成像条件(9)式成像结果;(b)利用本文给出成像条件(20)式的成像结果;(c)利用不带振幅校正的角度域滤波(18)式的成像结果;(d)、(e)、(f)分别为利用不带角度滤波的ray+Kirchhoff近似下真振幅成像条件(9)式成像结果(蓝色实线)与利用本文提出的成像条件(20)式成像结果(红色实线)第40、120、450道振幅对比. Fig. 5 Comparison of images and channels using different imaging conditions of the Marmousi model (a) True-amplitude imaging conditions (9) imaging results using ray+Kirchhoff approximation without angle filtering; (b) Imaging results using the proposed imaging condition (20); (c) Imaging results using scattering-angle filtering (18) without amplitude correction; (d), (e), and (f) are the amplitude comparison of No. 40, 120, and 450 trace of the images using true-amplitude imaging conditions (9) (blue solid line) with ray+Kirchhoff approximation without angle filtering, and the imaging conditions proposed in this paper (20) (red solid line) correspondingly.

为了进一步验证本文提出成像条件对复杂模型的有效性和适应性,我们设计了一个低速模型.速度模型如图 6a所示.模型在中间部分包含明显的低速异常.逆时偏移采用的平滑速度如图 6d所示.图 6b图 6e为采用不同成像条件的逆时偏移成像结果,很明显,通过应用本文提出的新的反演保真成像条件,可以实现对更深层的反射层振幅的改进.为了定量估计成像条件的有效性,分别选择模型不同空间位置两道成像结果与真实零角度反射系数进行比较,对比结果如图 6c6f所示.可以看出利用本文成像条件(20)式中估计出的带限反射系数峰值振幅和理论值匹配较好,表明本文提出的真振幅RTM成像条件在对储层参数进行定量估计方面具有一定潜力.

图 6 低速模型不同成像条件成像及其抽道对比结果 (a)真实的低速模型; (d) RTM计算用的平滑模型;(b)利用常规ray+Kirchhoff反演成像条件(9)式逆时偏移成像结果;(e)利用本文给出成像条件(20)式的成像结果;(c)利用本文成像条件计算的反射系数第60道结果(红色虚线)与理论零角度反射系数(蓝色实线)对比;(f)利用本文成像条件计算的反射系数第399道结果(红色虚线)与理论零角度反射系数(蓝色实线)对比. Fig. 6 Comparison of RTM images and selected traces with different imaging conditions for low velocity model (a) The true low velocity model; (d) The smoothed model for RTM calculation; (b) RTM stack image using conventional ray+Kirchhoff inversion based imaging condition (9); (e) RTM stack image using RTM imaging condition (20) as proposed in this paper; (c) The reflection coefficient calculated trace No.60 (red dashed line) compared with the theoretical reflection coefficient (blue solid line); (f) The reflection coefficient calculated trace No.399 (red dashed line) compared with the theoretical reflection coefficient (blue solid line).
4 结论与讨论

本文提出了一种新的针对反射波的逆时偏移真振幅成像条件,该反演成像条件严格遵循了ray+Kirchhoff真振幅成像公式,保证了成像结果与射线类渐进反演结果的一致性和计算的稳健性.同时,本文提出的成像条件能够在数据以反射波为主情况下,在RTM成像过程中自动消除低波数噪声.本文成像条件引入了振幅校正算子保证低波数噪声去除过程中目标反射层成像结果波形的保真.本文提出成像条件为后续对成像结果进一步的保真处理奠定了很好的基础.

本文讨论的真振幅RTM成像条件定义在经典的高频渐进反演框架下,该框架仅仅校正波传播过程的几何扩散效应,并不考虑由于观测系统非规则、有限孔径采样以及地下介质非均匀引起的成像点处照明不均问题(Chavent and Plessix, 1999).同时,经典的真振幅成像也忽略了对透射损失的补偿(Frijlink and Wapenaar, 2004).如何在校正几何扩散基础上,校正这些被忽视的波传播效应以进一步提高成像结果保真度是本文下一步研究方向之一.另一方面,本文初步数值实验表明,利用本文给出的真振幅RTM成像条件同样能够保真地输出波场相关中的低波数能量部分,对这部分能量在背景速度更新中的作用,也是值得下一步研究的内容.

References
Beydoun W B, Jin S D. 1994. Born or Kirchoff migration/inversion:what is the Earth's point of view?. Proceedings of SPIE-Mathematical Methods in Geophysical Imaging Ⅱ, 2301: 91-93.
Beylkin G. 1985. Imaging of discontinuities in the inverse scattering problem by inversion of a causal generalized Radon transform. Journal of Mathematical Physics, 26(1): 99-108. DOI:10.1063/1.526755
Bleistein N. 1987. On the imaging of reflectors in the earth. Geophysics, 52(7): 931-942. DOI:10.1190/1.1442363
Bleistein N, Stockwell J W Jr, Cohen J K. 2001. Mathematics of Multidimensional Seismic Inversion. New York: Springer-Verlag Inc.
Bleistein N, Zhang Y, Xu S, et al. 2005. Migration/inversion:think image point coordinates, process in acquisition surface coordinates. Inverse Problems, 21(5): 1715-1744. DOI:10.1088/0266-5611/21/5/013
Brandsberg-Dahl S, De Hoop M V, Ursin B. 2003. Focusing in dip and AVA compensation on scattering-angle/azimuth common image gathers. Geophysics, 68(1): 232-254. DOI:10.1190/1.1543210
Brandsberg-Dahl S, Chemingui N, Whitmore D, et al. 2013. 3D RTM angle gathers using an inverse scattering imaging condition.//83rd Ann. Internat Mtg., Soc. Expi. Geophys.. Expanded Abstracts, 3958-3962.
Červený V. 2001. Seismic Ray Theory. Cambridge: Cambridge University Press.
Chavent G, Plessix R E. 1999. An optimal true-amplitude least-squares prestack depth-migration operator. Geophysics, 64(2): 508-515. DOI:10.1190/1.1444557
Díaz E, Sava P. 2016. Understanding the reverse time migration backscattering:Noise or signal?. Geophysical Prospecting, 64(3): 581-594. DOI:10.1111/1365-2478.12232
Devaney A J. 2012. Mathematical Foundations of Imaging, Tomography and Wavefield Inversion. Cambridge: Cambridge University Press.
Esmersoy C, Miller D. 1989. Backprojection versus backpropagation in multidimensional linearized inversion. Geophysics, 54(7): 921-926. DOI:10.1190/1.1442722
Fletcher R F, Fowler P, Kitchenside P, et al. 2005. Suppressing artifacts in prestack reverse time migration.//85th Ann. Internat Mtg., Soc. Expi. Geophys.. Expanded Abstracts, 2049-2051.
Forgues E, Lambaré G. 1997. Parameterization study for acoustic and elastic ray+Born inversion. Journal of Seismic Exploration, 6(2-3): 253-278.
Frijlink M, Wapenaar K. 2004. Migration with correction for transmission losses in arbitrary inhomogeneous media.//74th Ann. Internat Mtg., Soc. Expi. Geophys.. Expanded Abstracts, 1029-1032. https://www.researchgate.net/publication/240736110_Migration_with_correction_for_transmission_losses_in_arbitrary_inhomogeneous_media
Gray S H. 1997. True-amplitude seismic migration:A comparison of three approaches. Geophysics, 62(3): 929-936. DOI:10.1190/1.1444200
Gray S H, Bleistein N. 2009. True-amplitude Gaussian-beam migration. Geophysics, 74(2): S11-S23. DOI:10.1190/1.3052116
Guitton A, Kaelin B, Biondi B. 2006. Least-squares attenuation of reverse-time-migration artifacts. Geophysics, 72(1): S19-S23.
Hou J, Symes W. 2015. An approximate inverse to the extended Born modeling operator. Geophysics, 80(6): R331-R349. DOI:10.1190/geo2014-0592.1
Jin S D, Madariaga R, Virieux J, et al. 1992. Two-dimensional asymptotic iterative elastic inversion. Geophysical Journal International, 108(2): 575-588. DOI:10.1111/j.1365-246X.1992.tb04637.x
Joncour F, Lambaré G, Svay-Lucas J. 2011. Preserved-amplitude angle domain migration by shot-receiver wavefield continuation. Geophysical Prospecting, 59(2): 256-268. DOI:10.1111/j.1365-2478.2010.00921.x
Kühl H, Sacchi M D. 2003. Least-squares wave-equation migration for AVP/AVA inversion. Geophysics, 68(1): 262-273. DOI:10.1190/1.1543212
Kiyashchenko D, Plessix R E, Kashtan B, et al. 2007. A modified imaging principle for true-amplitude wave-equation migration. Geophysical Journal International, 168(3): 1093-1104. DOI:10.1111/j.1365-246X.2006.03187.x
Lambare G, Virieux J, Madariaga R, et al. 1992. Iterative asymptotic inversion in the acoustic approximation. Geophysics, 57(9): 1138-1154. DOI:10.1190/1.1443328
Li Y B, Chauris H. 2018. Coupling direct inversion to common-shot image-domain velocity analysis. Geophysics, 83(5): R497-R514. DOI:10.1190/geo2017-0825.1
Liu F Q, Zhang G Q, Morton S A, et al. 2011. An effective imaging condition for reverse-time migration using wavefield decomposition. Geophysics, 76(1): S29-S39. DOI:10.1190/1.3533914
Miller D, Oristaglio M, Beylkin G. 1987. A new slant on seismic imaging:Migration and integral geometry. Geophysics, 52(7): 943-964. DOI:10.1190/1.1442364
Operto M, Xu S, Lambaré G. 2000. Can we quantitatively image complex structures with rays?. Geophysics, 65(4): 1223-1238. DOI:10.1190/1.1444814
Rocha D, Tanushev N, Sava P. 2016. Acoustic wavefield imaging using the energy norm. Geophysics, 81(4): S151-S163. DOI:10.1190/geo2015-0486.1
Shuey R T. 1985. A simplification of the Zoeppritz equations. Geophysics, 50(4): 609-614. DOI:10.1190/1.1441936
Stolk C, De Hoop M, Op't Root T J P M. 2009. Inverse scattering of seismic data in the reverse time migration (RTM) approach.//Proceedings of the Project Review, Geo-Mathematical Imaging Group, Purdue University, 1: 91-108. https://www.researchgate.net/publication/228917838_Inverse_scattering_of_seismic_data_in_the_reverse_time_migration_rtm_approach
Ten Kroode A P E, Smit D J, Verdel A. 1994. Linearized inversed scattering in the presence of caustics. Proceedings of SPIE-Mathematical Methods in Geophysical Imaging Ⅱ, 2301: 28-43. DOI:10.1117/12.187493
Thierry P, Operto S, Lambaré G. 1998. Fast 2-D ray+Born migration/inversion in complex media. Geophysics, 64(1): 162-181.
Tygel M, Schleicher J, Hubral P. 1994. Pulse distortion in depth migration. Geophysics, 59(10): 1561-1569. DOI:10.1190/1.1443545
Wapenaar C P A, Berkhout A J. 1989. Elastic Wave Field Extrapolation. Elsevier.
Whitmore N D, Crawley S. 2012. Applications of RTM inverse scattering imaging conditions.//82nd Ann. Internat Mtg., Soc. Expi. Geophys.. Expanded Abstracts.
Wu R S, Xie X B. 2002. Extracting angle domain information from migrated wavefield.//72nd Ann. Internat Mtg., Soc. Expi. Geophys.. Expanded Abstracts. https://www.researchgate.net/publication/238438533_Extracting_angle_domain_information_from_migrated_wavefield
Xu S, Chauris H G, Lambaré G, et al. 2001. Common-angle migration:A strategy for imaging complex media. Geophysics, 66(6): 1877-1894. DOI:10.1190/1.1487131
Xu S, Zhang Y, Tang B. 2011. 3D angle gathers from reverse time migration. Geophysics, 76(2): S77-S92. DOI:10.1190/1.3536527
Zhang Y, Xu S, Bleistein N, et al. 2007. True amplitude angle domain common image gathers from one-way wave equation migrations. Geophysics, 72(1): S49-S58. DOI:10.1190/1.2399371
Zhang Y, Sun J. 2009. Practical issues of reverse time migration:True amplitude gathers, noise removal and harmonic-source encoding. First Break, 26(1): 29-35.
Zhang Y, Ratcliffe A, Roberts G, et al. 2014. Amplitude-preserving reverse time migration:From reflectivity to velocity and impedance inversion. Geophysics, 79(6): S271-S283. DOI:10.1190/geo2013-0460.1
Zhou Y, Wang H Z. 2016. Efficient wave-mode separation in vertical transversely isotropic media. Geophysics, 82(2): C35-C47.