2021, Vol. 64
2. 同济大学海洋与地球科学学院, 上海 200092;
3. 中国石油勘探开发研究院西北分院, 兰州 730020
2. School of Ocean and Earth Science, Tongji University, Shanghai 200092, China;
3. Research Institute of Petroleum Exploration and Development-Northwest, PetroChina, Lanzhou 730020, China
通过求解双向波方程，逆时偏移(RTM)能够不受传播倾角的限制对陡倾角地质体进行成像.同时，由于在波场反传播过程中校正了几何扩散效应，因此逆时偏移是一种经典意义下的真振幅成像方法(Zhang and Sun, 2009).对于复杂地质体的精确成像，除了波场传播的几何扩散效应，由于观测采集系统不规则、不完备采样以及地下复杂构造等因素，在成像点处通常还会存在照明不均匀的现象(Xie et al., 2006).因此，在逆时偏移中进一步考虑照明补偿会在较大程度上提高对地质体的成像精度.
照明补偿成像可以在偏移成像结果上作用Hessian算子或与其相关的点扩散函数或照明算子的逆来实现.Audebert等(2000)通过在角度域统计射线打击次数近似照明算子实现带照明补偿的克希霍夫偏移成像方法，Gelius等(2002)、Lecomte(2008)通过射线计算点扩散函数，并利用点扩散函数对成像结果进行反褶积来实现对成像结果照明不均匀的补偿.Wu和Chen(2002)以及Xie等(2005)讨论了基于单程波波动方程的照明算子计算及对应的照明补偿方法.Wu等(2004)、陈生昌等(2007)、Cao和Wu(2009)进一步将照明算子分解到入射/散射角度域进行照明补偿成像，取得了更优的照明补偿效果.Yan等(2014)利用逆时偏移在倾角域计算照明算子并提出了一套照明补偿RTM方案.Yan和Xie(2016)利用Poynting矢量对波场进行角度分解，并以此为基础发展了一套角度域照明补偿RTM方法.类似地，Hu等(2001)、Plessix和Mulder(2004)、Ren等(2011)讨论了在最小二乘框架下基于显式计算Hessian及其近似逆算子的照明补偿成像方法.上述波动方程类照明补偿成像方法通常需要利用波动方程显式计算格林函数，对于大规模的三维成像问题，在计算量和存储量上具有不小的挑战.
为了提高利用波动方程进行照明补偿成像的计算效率，Rickett(2003)、Guitton(2004)提出了利用偏移-反偏移的方式估计近似Hessian算子进行照明补偿成像，该方式通过增加一次反偏移和一次偏移的计算代价来避免对格林函数的显式计算和储存.Gherasim等(2010)、Shen等(2011)、Li等(2012)将此种方式推广到了角度域.相较于直接计算格林函数构造照明算子或Hessian算子，利用偏移-反偏移的方式进行照明补偿成像能在一定程度上提高计算效率.对于基于RTM的照明补偿成像，由于一次三维RTM偏移的计算量较为庞大(刘红伟等，2010；许璐等，2019)，因此利用偏移-反偏移的方式进行大规模数据RTM照明补偿成像在计算量上仍然面临较大压力.
另一方面，在最小二乘框架下也可以利用迭代求取Hessian逆算子或其近似算子进行成像的照明补偿，即迭代实现的最小二乘类偏移成像方法技术.该类方法技术的发展经历了射线类迭代最小二乘偏移成像方法(Schuster，1993；Nemeth et al., 1999；Duque et al., 2000；王彦飞等，2009；黄建平等，2013)、单程波算子类迭代最小二乘偏移成像方法(Kühl and Sacchi, 2003; 杨其强和张叔伦，2008; 沈雄君和刘能超，2012；周华敏等，2014)以及迭代实现的最小二乘RTM偏移成像方法(Dai et al., 2012；刘学建和刘依克，2016；李振春等, 2017a, b；方修政等，2018；张攀和毛伟建，2018；周东红等，2020).在实际资料处理中，利用迭代方式估计Hessian逆算子的精度和稳定性易受到噪声不满足先验分布、震源子波空变且未知、初始模型不满足近似线性化反演要求等多种因素影响，从而使得迭代最小二乘类成像方法出现迭代收敛慢甚至不收敛的情况，在一定程度上增加了该方法的实际应用难度.
本文提出一种基于RTM实现的高效照明补偿成像方法.在高频近似意义下，首先从波传播的几何路径出发，定义衡量成像点处照明均匀的雅克比矩阵，并以此为基础进一步给出照明补偿算子的表达式.然后通过波场的边界积分表达，将震源端照明补偿算子利用外推波场进行表达，避免了震源端照明算子中对格林函数显式计算的需求.对于检波点端照明补偿算子，通过引入对参数场的合理近似，将检波点端照明补偿算子中对检波器坐标积分项进行显式表达，在较大程度上提高了计算效率.利用合成数据和实际数据验证了本文提出方法的正确性.
1 方法原理本小节给出基于RTM的高效照明补偿成像方法.首先引入衡量地下照明规则程度的雅克比行列式，然后基于该雅克比行列式给出照明补偿算子的定义.进一步，讨论利用RTM计算引擎高效实现该照明补偿算子的具体算法.
1.1 照明补偿算子为了引入照明雅克比行列式，首先对下文推导中涉及的符号进行定义说明，文中涉及主要符号的几何定义如图 1所示.图 1a为本文涉及的射线传播参数定义.其中ps和pr分别为源端和检波点端在地下成像点x处射线慢度矢量，其定义为震源及检波点两端波场走时场梯度.nqr为照明倾角矢量，其定义为ps与ps的矢量和(Forgues and Lambaré, 1997).ps和pr与竖直方向的夹角βs、βr分别定义为震源及检波点两端的起飞角(Operto et al., 2000).ps和pr之间的夹角θs为成像点处的散射角，代表当前观测系统下该成像点的角度照明范围，照明倾角矢量nqr与竖直方向夹角ψx为照明倾角，代表当前观测系统下角度照明范围的中心位置(Audebert et al., 2005).图 1b为图 1a在成像点x处的局部放大图.
|
图 1 (a) 射线传播参数表征及角度定义.ps、pr分别为源检两端出射射线在地下点x处的慢度矢量，nqr=ps+pr为照明倾角矢量，θs地下点x处的散射角且θs=βr－βs，βs、βr分别为地下点x处的起飞角，ψ为照明倾角且ψ= (βr+βs)/2；(b)图(a)在x处局部放大显示 Fig. 1 (a) Ray parameters for 2D asymptotic inversion imaging and illumination compensation. ps and pr are slowness vectors coming from the source and the receiver, respectively. nqr is the normalized illumination dip vector. βs and βr are the takeoff angles on source and receiver side, respectively. θs=βr－βs is the scattering angle. ψ=(βr+βs)/2 is the illumination dip angle; (b) Enlargement of (a) at x |
对于地下成像点处照明均匀程度，本文利用局部角度域扰动相对于地表坐标扰动量大小来进行度量.图 2展示了在直射线假设情况下利用层状介质和局部速度异常介质对本文给出照明均匀程度定义的几何解释.图 2a为层状介质直射线假设下不同深度局部角度域扰动与地表坐标扰动的关系示意图.从图 2a可以看出当地表坐标具有相同扰动dx时，深部反射地层局部角度域扰动dψ2要小于浅部扰动dψ1.这意味着对于固定观测系统，深部反射地层的角度覆盖范围要小于浅部反射地层覆盖范围，为了使得不同深度反射地层具有一致的照明强度，深部反射地层需要更多的照明补偿.图 2b、c分别为存在及不存在局部高速异常时照明射线路径图.可以看出当地下局部角度域存在相同扰动dψ1时，由于速度异常体对照明射线的偏折作用，存在速度异常体模型中照明射线在地表坐标扰动量dx2要远大于匀速介质中照明射线在地表坐标扰动量dx1.这表明对于固定偏移距的观测系统，速度异常体下方更加难以被地震波场照明，需要更多的照明补偿才能够消除照明不均衡对成像结果的影响.据以上分析，本文定义二维情况下衡量成像点处照明规则程度的雅克比行列式为：
|
(1) |
|
图 2 (a) 在层状介质直射线假设下，不同深度局部角度与扰动域地表扰动的关系示意图；(b) 在常速介质中局部角度域扰动与地表扰动关系示意图；(c) 在包含局部异常情景下，局部角度域扰动与地表扰动关系示意图 Fig. 2 (a) Illustration of relation between local angle and surface coordinate perturbation in layered media with straight ray approximation; (b) Illustration of relation between local angle and surface coordinate perturbation in constant-velocity media; (c) Illustration of relation between local angle and surface coordinate perturbation in constant-velocity media with local anomaly |
其中ψx为成像点x处的照明倾角，|dψx/dxs|及|dψx/dxr|分别代表照明倾角相对地表炮点坐标xs及检波器坐标xr扰动的变化.注意到式(1)右端行列式结果为标量，标量Gs和Gr即为本文给出震源和检波点端衡量照明规则程度的雅克比行列式.
从式(1)可以看出，给定相同照明倾角扰动量dψx，如果地表坐标dxs及dxr扰动越大，则Gs和Gr越小，代表对于固定地震采集系统，地下照明规则程度越差，所需的照明补偿量越大.基于以上分析，由式(1)出发可以进一步定义成像点x处的照明补偿算子Wx的形式为：
|
(2) |
公式(2)中对检波器的积分代表考虑共炮数据中所有检波器接收数据的照明补偿.对照公式(1)和(2)可以看出，当震源和检波点两端照明越不规则，照明补偿算子Wx越大，照明补偿效果越强.
公式(2)即为本文给出的照明补偿算子抽象形式，为了能够在RTM过程中高效地实施照明补偿算子(2)，需要对公式(2)中元素做更具体的表征.利用图 1所示的照明倾角与起飞角的关系ψx=(βr+βs)/2，并注意到震源端射线慢度矢量对应的起飞角βs仅仅与炮点坐标xs有关，检波点端对应的起飞角βr仅仅与检波点坐标xr有关，可以将式(2)改写为：
|
(3) |
据Bleistein等(2005)，在高频近似下，起飞角扰动相对于地表坐标扰动关系为：
|
(4) |
其中A(x, xs)为高频近似下波场从炮点xs出发传播到成像点x处的振幅.c0(xs)为炮点xs处的背景速度.
对于检波点端起飞角相对地表坐标的扰动关系，有类似式(4)的表达形式.将式(4)代入式(3)最终可得照明补偿算子Wx的具体表达形式为：
|
(5) |
其中：
|
(6) |
式(6)中A(xr, x)为高频近似下波场经过成像点x反射被位于xr检波点接收到的振幅，c0(xs)为检波点xr处的背景速度.
1.2 照明补偿算子的RTM实现照明补偿算子(5)和(6)式的计算包含了高频近似下波场的振幅，理论上需要利用动力学射线追踪进行计算(Červený，2001)，这会在逆时偏移过程中增加额外的计算量.本文通过引入波场的边界积分表达和对介质的合理假设，使得能在波场外推过程中高效地利用式(5)、(6)进行照明补偿成像.对于震源端照明算子Gs，有二维情况下震源端波场PF(x, xs, ω)高频近似下的瑞雷第二类积分解形式(Zhang et al., 2003；Devancy，2012)：
|
(7) |
将式(7)代入式(6)并做适当变形可得：
|
(8) |
式(8)即为利用外推波场表达的源端照明补偿校正算子.可以看出，由于源端波场PF(x, xs, ω)在波场外推中已知，式(8)的数值实现仅需额外计算成像点处震源端波场的起飞角βs，起飞角βs可以利用坡印廷矢量在波场外推中快速计算(Yoon and Marfurt, 2006)，因此照明补偿算子式(8)可以在震源端波场外推过程中高效地数值实现.对于检波点端照明补偿算子Gr，受到Plessix和Mulder(2004)工作的启发，采用缓变介质假设，将式(6)中对检波器坐标积分项进行显式表达.缓变介质的格林函数振幅项(Bleistein et al., 2001)为：
|
(9) |
将式(9)代入式(6)并做适当代数运算，最终可以得到检波点端照明补偿算子的形式为：
|
(10) |
其中x、z为成像点的空间坐标，xrmin、xrmax为当前炮数据最小和最大观测孔径.由于在RTM算法的具体实施过程中，成像点坐标以及当前炮数据最小和最大观测孔径均为已知量，因此检波点端照明补偿算子(10)同样可以在RTM成像过程中进行高效地计算.需要注意的是，由于在导出照明补偿算子过程中引入了介质缓变假设，因此本文给出的基于RTM的照明补偿成像方法适用于给定缓变背景模型情况下对模型扰动部分进行成像(Bleistein et al., 2001).
2 数值实验 2.1 模型资料测试首先利用一个四层速度模型来验证本文提出的照明补偿方法的有效性.将利用本文给出的照明补偿算子成像结果与Plessix和Mulder (2004)经典的基于对角Hessian照明补偿算子成像结果进行比较，特别地，对比Plessix和Mulder(2004)文中给出的第三类照明补偿算子成像结果，该类照明补偿算子相较于文中给出的其他几类照明补偿算子有更好的补偿效果(Kiyashchenko et al., 2007).图 3a为四层速度真实模型，四层速度分别为3000 m·s-1、4000 m·s-1、5000 m·s-1和6000 m·s-1.图 3b展示了成像算法采用的平滑模型，该平滑模型通过利用半径为10的高斯平滑滤波器获得.图 3c为在常密度假设情况下利用速度模型计算的真实零角度反射系数(Berteussen and Ursin, 1983).为了更好地对比分析不同照明补偿算子对成像结果的影响，计算了数值实验采用的观测系统的总照明补偿强度，如图 4所示.图 4a、b分别为利用Plessix和Mulder (2004)给出的照明补偿算子及本文给出的照明补偿算子计算得到的总照明补偿强度.可以看出，本文方法计算的总照明补偿强度分布和直接利用近似对角Hessian逆算子计算出的照明补偿强度(Plessix and Mulder, 2004)分布整体上具有较好的一致性，这在某种程度上验证了本文给出的照明补偿算子确为对角Hessian逆算子在高频近似下的一个近似.此外，从图 4a、b的对比可以看出，相较于Plessix和Mulder(2004)给出的对角Hessian逆算子的近似方式，利用本文方法计算得到的总照明补偿强度和速度关系更加密切，不同深度照明倾角信息体现也更加明显.
|
图 3 (a) 四层真实速度模型；(b) 四层平滑速度模型；(c) 由速度场计算得到的真实零角度反射系数模型 Fig. 3 (a) The true four-layer velocity model; (b) The smoothed four-layer velocity for imaging; (c) The true normal reflection coefficient for four-layer model |
|
图 4 四层速度模型的总照明补偿强度对比 (a) 利用Plessix和Mulder (2004)给出的照明补偿算子计算得到的照明补偿强度；(b) 利用本文给出照明补偿算子(8)和(10)式计算得到的照明补偿强度. Fig. 4 Comparison of total illumination compensation intensity for four-layer model (a) The illumination compensation intensity calculated using the operators proposed in Plessix and Mulder (2004); (b) The illumination compensation intensity calculated using the operators (8) and (10). |
图 5a-c分别展示了利用不带照明补偿的RTM算法、利用Plessix和Mulder(2004)照明补偿算子以及本文照明补偿算子(8)和(10)进行照明补偿成像结果.在不同照明补偿策略的RTM成像算法中，本文利用逆散射成像条件来消除RTM中的低波数噪声(Whitmore and Crawley, 2012; Zhou and Wang, 2017).可以看出，相较于不带照明补偿的RTM算法，利用Plessix和Mulder(2004)以及本文的照明补偿算子进行照明补偿成像均能有效提升深层成像能量，使得成像剖面能量更加均衡.为了定量地对比不同照明补偿算子效果，对利用不同照明补偿策略得到的成像剖面进行抽道并与真实零角度反射系数进行对比，结果展示在图 5d-f中.从图 5d可以看出，深层反射层由于照明不足，不做照明补偿直接计算得到反射系数峰值振幅与真实值误差较大.利用Plessix和Mulder (2004)给出照明补偿算子计算得到深层反射系数峰值振幅与理论值的匹配程度有较大程度的提升，如图 5e所示.利用本文给出照明补偿算子计算得到的反射系数抽道对比结果展示在图 5f中，可以看出利用本文给出照明补偿算子能进一步提升深层反射系数计算精度.为了测试本文给出照明补偿算子(8)和(10)在数据包含噪声情况下的计算稳定性，在模拟炮数据中加入满足高斯分布的噪声，原始数据与含噪数据的抽道对比如图 6a所示.图 6c、d分别为利用本文照明补偿算子在原始和含噪数据上单炮成像结果.可以看出相较于不包含噪声的原始数据，含噪声数据单炮照明补偿成像结果包含由数据噪声引起的额外干扰，这些额外成像干扰能量相较于真实反射层成像结果较弱.在所有炮数据成像结果叠加后，这些由数据噪声引起额外成像干扰会被进一步压制.图 6b为利用本文照明补偿算子在原始和含噪数据上所有炮数据叠加成像结果的抽道对比，可以较为清晰看出最终照明补偿成像结果几乎不受数据噪声的影响，这验证了本文给出照明补偿算子对含噪数据的数值计算稳健性.
|
图 5 四层模型不同照明补偿策略成像结果对比 (a) 不带照明补偿补偿的RTM成像结果；(b) 利用Plessix和Mulder(2004)给出的照明补偿算子得到的RTM成像结果；(c) 利用本文给出的照明补偿算子得到的RTM成像结果；(d)-(f) 分别为(a)-(c)抽取中间道的成像结果(红色实线)与真实零角度反射系数(蓝色实线)对比. Fig. 5 Comparisons of imaging results with different illumination compensation operators for four-layer model (a), (b) and (c) Imaging results without illumination, with compensation operator of Plessix and Mulder (2004) and the compensation operator proposed in this paper, respectively; (d), (e) and (f) Comparison for selected central traces between imaging results (red solid line) and theoretical normal reflection coefficient values (blue solid line) corresponding to (a), (b) and (c), respectively. |
|
图 6 四层速度模型包含噪声数据测试 (a) 包含噪数据(蓝色实线)与不含噪数据(红色实线)抽道对比；(b) 包含噪数据叠加成像结果(蓝色实线)与不含音数据叠加成像结果(红色实线)抽道对比；(c) 不含噪数据单炮成像结果；(d) 包含噪数据单炮成像结果. Fig. 6 Noisy data test for four-layer model (a) Comparison of selected traces for noisy data (blue solid line) and clean data (red solid line); (b) Comparison of selected imaged traces with noisy data (blue solid line) and clean data (red solid line); (c) Imaging result with single shot clean data; (d) Imaging result with single shot noisy data. |
为了进一步验证本文照明补偿成像方法对复杂模型的有效性，利用Marmousi模型进行了照明补偿成像测试.图 7a、b分别为Marmousi真实速度以及零角度反射系数模型，利用半径为10的高斯平滑滤波器作用在真实模型上得到平滑模型用于成像.利用不同照明补偿算子计算得到总照明补偿强度展示在图 8中.和四层速度模型测试结论类似，利用本文给出照明补偿算子及Plessix和Mulder(2004)给出照明补偿算子计算得到的总照明补偿强度在整体分布具有较好的一致性，同时利用本文给出照明补偿算子计算得到总照明补偿强度和速度构造有更好的一致性.图 9a-c分别为不带照明补偿、利用Plessix和Mulder (2004)给出补偿算子和利用本文给出照明补偿算子进行基于RTM的照明补偿成像叠加结果.对于不做照明补偿的RTM成像结果，可以清晰看出由于震源附近强波场能量相关导致的成像噪声.通过在RTM过程中施加Plessix和Mulder(2004)给出的照明补偿算子进行成像，震源附近强能量噪声得到了较好的压制，整个成像剖面不同深度能量更加均衡，如图 9b所示.利用本文给出的照明补偿算子成像叠加结果展示在图 9c中，可以较明显看出震源相关噪声得到了进一步压制，深层成像振幅得到了进一步提升.图 9d-f为图 9a-c对应的局部放大结果，可以看出利用本文给出的照明补偿算子对复杂构造区域由于波场未完全干涉形成的成像噪声也具有很好的压制作用，可以得到信噪比更好的复杂区域成像结果. Marmousi模型不同照明补偿策略成像结果抽道对比如图 10所示.图 10a为不做照明补偿成像结果和真实零角度反射系数对比，可以清晰地看出浅层震源处强相关噪声对成像结果的影响.同时，随着成像深度的增加，不做照明补偿成像结果振幅峰值与理论零角度反射系数误差逐渐增大，这与有效照明范围随着成像深度增加逐步减小的趋势是一致的.利用Plessix和Mulder(2004)给出照明补偿算子成像抽道对比结果显示在图 10b中，可以看出利用该照明补偿算子能明显压制浅层震源强相关噪声，有效提升中深层成像结果峰值振幅与真实零角度反射系数匹配程度.图 10c为利用本文照明补偿算子(8)和(10)成像抽道对比结果，相较于图 10b，浅层震源强相关噪声得到了进一步衰减，中深层成像结果峰值振幅更加接近于理论零角度反射系数.
|
图 7 (a) Marmousi真实速度模型；(b) Marmousi模型真实反射系数 Fig. 7 (a) The true Marmousi velocity model; (b) The true reflection coefficient for Marmousi model |
|
图 8 Marmousi速度模型的总照明补偿强度对比 (a) 利用Plessix和Mulder (2004)给出的照明补偿算子计算得到的照明补偿强度；(b) 利用本文给出照明补偿算子(8)和(10)式计算得到的照明补偿强度. Fig. 8 Comparison of total illumination compensation intensity for Marmousi model (a) The illumination compensation intensity calculated using the operators proposed by Plessix and Mulder (2004); (b) The illumination compensation intensity calculated using the operators (8) and (10). |
|
图 9 Marmousi模型不同照明补偿策略成像结果对比 (a) 不带照明补偿补偿的RTM成像结果；(b) 利用Plessix和Mulder(2004)给出的照明补偿算子得到的RTM成像结果；(c) 利用本文给出的照明补偿算子得到的RTM成像结果；(d)-(f) 分别为(a)-(c)的局部放大对比图. Fig. 9 Comparison of imaging results with different illumination compensation operators for Marmousi model (a), (b) and (c)Imaging results without illumination, with compensation operator of Plessix and Mulder (2004) and the compensation operator proposed in this paper, respectively. (d), (e) and (f) Enlarged display of (a), (b) and (c), respectively. |
|
图 10 Marmousi模型不同照明补偿策略成像5250 m处抽道对比结果 (a) 不带照明补偿补偿的RTM成像结果(红色实线)与理论零角度反射系数(蓝色实线)对比结果；(b) 基于Plessix和Mulder(2004)给出照明补偿算子RTM成像结果(红色实线)与理论零角度反射系数(蓝色实线)对比结果；(c) 基于本文给出的照明补偿算子RTM成像结果(红色实线)与理论零角度反射系数(蓝色实线)对比结果. Fig. 10 Comparison of selected traces between imaging results (red solid line) and theoretical normal reflection coefficient values (blue solid line) (a) RTM imaging results calculated without illumination compensation; (b) RTM imaging results calculated with illumination compensation operator of Plessix and Mulder (2004); (c) RTM imaging results calculated with illumination compensation operator proposed in this paper. |
为了检验本文照明补偿成像方法在实际资料上的应用效果，利用西部某二维实际资料进行测试.图 11展示了实际资料的速度场.利用不同照明补偿策略得到的RTM成像叠加剖面展示在图 12a-c中，叠加成像部分放大结果展示在图 12d-f中.可以看出，未进行照明补偿RTM成像叠加结果整个剖面存在由于波场能量未完全干涉造成的成像噪声，且深层反射能量不聚焦.利用Plessix和Mulder(2004)给出照明补偿算子进行照明补偿成像能够较好地提升中深层反射能量，减少非相干成像噪声干扰.利用本文提出的照明补偿算子进行照明补偿成像在成像噪声压制，成像剖面能量均衡以及深层反射层成像结果聚焦性等方面有较明显改进效果，相较于利用Plessix和Mulder(2004)给出照明补偿算子计算结果在成像质量方面有进一步的提升.
|
图 11 中国西部某实际数据速度场 Fig. 11 Velocity field of real data from western China |
|
图 12 (a) 不带照明补偿补偿的RTM成像结果；(b) 利用Plessix和Mulder(2004)给出的照明补偿算子得到的RTM成像结果; (c) 利用本文的照明补偿算子得到的RTM成像结果; (d)-(f) 分别为(a)-(c)对应的局部放大显示结果 Fig. 12 (a), (b) and (c) Imaging results without illumination, with illumination compensation operator given by Plessix and Mulder (2004) and the compensation operator proposed in this paper, respectively; (d), (e) and (f) Zoomed display of partial imaging result of (a), (b) and (c), respectively. |
本文给出的照明补偿算子为对角Hessian逆算子的一个近似，在高频近似意义下矫正了成像点处照明不均匀的影响，提高了成像精度.由于本文给出的照明补偿算子并非完整Hessian算子的逆，因此在利用本文的照明补偿算子进行照明补偿成像后，会有残余的波场传播效应有待进一步消除.对由于照明补偿算子近似Hessian逆算子不精确造成的残留波场传播效应问题，可将本文给出照明补偿算子融入到迭代求解的最小二乘类成像框架中，利用其迭代最小二乘成像方法中构造的预条件梯度项，在后续的迭代过程中消除残余的波场传播效应，进一步提高成像精度.Burgess和Warner(2015)利用Plessix和Mulder(2004)给出的照明补偿算子构造预条件梯度项并应用到全波形反演(FWI)中，成功加速了FWI收敛效率.类似地，从数值实验对比可以看出，利用本文给出照明补偿算子构造预条件梯度项也有望提高最小二乘类成像方法的收敛速度，这是本文后续研究方向之一.
4 结论本文给出了一种新的RTM照明补偿方法.该照明补偿方法利用地下角度域相对于地表坐标的扰动量来刻画局部成像点的照明均匀程度，据此定义照明补偿算子进行成像过程中的照明补偿.通过高频近似下的波场边界积分表达及对介质的合理近似，给出了照明补偿算子基于RTM的高效实现方式.模型和实际资料的测试对比表明，本文提出的照明补偿算子能够在较大程度上提升RTM成像质量.
Audebert F, Froidevaux P, Huard I, et al. 2000. A multi-angle toolbox for restored amplitude images and AVA-Gathers. //70th Ann. Internat Mtg., Soc. Expi. Geophys. . Expanded Abstracts.
|
Audebert F, Nicoletis L, Froidevaux P. 2005. Regularization of illumination in angle domains-A key to true amplitude migration. The Leading Edge, 24(6): 643-654. DOI:10.1190/1.1946222 |
Berteussen K A, Ursin B. 1983. Approximate computation of the acoustic impedance from seismic data. Geophysics, 48(10): 1351-1358. DOI:10.1190/1.1441415 |
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 |
Burgess T, Warner M. 2015. Preconditioning FWI with approximate receiver Green's functions. //85th Ann. Internat Mtg., Soc. Expi. Geophys. . Expanded Abstracts.
|
Cao J, Wu R S. 2009. Fast acquisition aperture correction in prestack depth migration using beamlet decomposition. Geophysics, 74(4): S67-S74. DOI:10.1190/1.3116284 |
Červený V. 2001. Seismic Ray Theory. Cambridge: Cambridge University Press.
|
Chen S C, Ma Z T, Wu R S. 2007. Illumination compensation for wave equation migration shadow. Chinese Journal of Geophysics (in Chinese), 50(3): 844-850. |
Dai W, Fowler P, Schuster G T. 2012. Multi-source least-squares reverse time migration. Geophysical Prospecting, 60(4): 681-695. DOI:10.1111/j.1365-2478.2012.01092.x |
Devaney A J. 2012. Mathematical Foundations of Imaging, Tomography and Wavefield Inversion. Cambridge: Cambridge University Press.
|
Duquet B, Marfurt K J, Dellinger J A. 2000. Kirchhoff modeling, inversion for reflectivity, and subsurface illumination. Geophysics, 65(4): 1195-1209. DOI:10.1190/1.1444812 |
Fang X Z, Niu F L, Wu D. 2018. Least-squares reverse-time migration enhanced with the inverse scattering imaging condition. Chinese Journal of Geophysics (in Chinese), 61(9): 3770-3782. DOI:10.6038/cjg2018L0721 |
Forgues E, Lambaré G. 1997. Parameterization study for acoustic and elastic ray + born inversion. Journal of Seismic Exploration, 6(2-3): 253-277. |
Gelius L J, Lecomte I, Tabti H. 2002. Analysis of the resolution function in seismic prestack depth imaging. Geophysical Prospecting, 50(5): 505-515. DOI:10.1046/j.1365-2478.2002.00331.x |
Gherasim M, Albertin U, Nolte B, et al. 2010. Wave-equation angle-based illumination weighting for optimized subsalt imaging. //80th Ann. Internat Mtg., Soc. Expi. Geophys. . Expanded Abstracts.
|
Guitton A. 2004. Amplitude and kinematic corrections of migrated images for nonunitary imaging operators. Geophysics, 69(4): 1017-1024. DOI:10.1190/1.1778244 |
Hu J X, Schuster G T, Valasek P A. 2001. Poststack migration deconvolution. Geophysics, 66(3): 939-952. DOI:10.1190/1.1444984 |
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. DOI:10.6038/cjg20130529 |
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 |
Lecomte I. 2008. Resolution and illumination analyses in PSDM: A ray-based approach. The Leading Edge, 27(5): 650-663. DOI:10.1190/1.2919584 |
Li Z C, Huang J Q, Huang J P, et al. 2017a. Fast least-squares reverse time migration based on plane-wave encoding for VTI media. Chinese Journal of Geophysics (in Chinese), 60(1): 240-257. DOI:10.6038/cjg20170120 |
Li Z C, Zhou L Y, Huang J P, et al. 2017b. Application of least square reverse time migration in the imaging of fractured-type carbonate reservoirs. Progress in Geophysics (in Chinese), 32(2): 664-671. DOI:10.6038/pg20170229 |
Li Z X, Tang B, Ji S. 2012. Subsalt illumination analysis using RTM 3D dip gathers. //82nd Ann. Internat Mtg., Soc. Expi. Geophys. . Expanded Abstracts.
|
Liu H W, Li B, Liu H, et al. 2010. The algorithm of high order finite difference pre-stack reverse time migration and GPU implementation. Chinese Journal of Geophysics (in Chinese), 53(7): 1725-1733. DOI:10.3969/j.issn.0001-5733.2010.07.024 |
Liu X J, Liu Y K. 2016. Least-squares reverse-time migration of surface-related multiples. Chinese Journal of Geophysics (in Chinese), 59(9): 3354-3365. DOI:10.6038/cjg20160919 |
Nemeth T, Wu C J, Schuster G T. 1999. Least-squares migration of incomplete reflection data. Geophysics, 64(1): 208-221. DOI:10.1190/1.1444517 |
Operto M S, Xu S, Lambaré G. 2000. Can we quantitatively image complex structures with rays?. Geophysics, 65(4): 1223-1238. DOI:10.1190/1.1444814 |
Plessix R E, Mulder W A. 2004. Frequency-domain finite-difference amplitude-preserving migration. Geophysical Journal International, 157(3): 975-987. DOI:10.1111/j.1365-246X.2004.02282.x |
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/j.1365-2478.2011.00947.x |
Rickett J E. 2003. Illumination-based normalization for wave-equation depth migration. Geophysics, 68(4): 1371-1379. DOI:10.1190/1.1598130 |
Schuster G T. 1993. Least-squares cross-well migration. //63rd Ann. Internat Mtg., Soc. Expi. Geophys. . Expanded Abstracts.
|
Shen H, Mothi S, Albertin U. 2011. Improving subsalt imaging with illumination-based weighting of RTM 3D angle gathers. //81st Ann. Internat Mtg., Soc. Expi. Geophys. . Expanded Abstracts.
|
Shen X J, Liu N C. 2012. Split-step least-squares migration. Progress in Geophysics (in Chinese), 27(2): 761-770. DOI:10.6038/j.issn.1004-2903.2012.02.044 |
Wang Y F, Yang C C, Duan Q L. 2009. On iterative regularization methods for migration deconvolution and inversion in seismic imaging. Chinese Journal of Geophysics (in Chinese), 52(6): 1615-1624. |
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, Chen L. 2002. Mapping directional illumination and acquisition-aperture efficacy by beamlet propagators. //72nd Ann. Internat Mtg., Soc. Expi. Geophys. . Expanded Abstracts.
|
Wu R S, Luo M Q, Chen S C, et al. 2004. Acquisition aperture correction in angle-domain and true-amplitude imaging for wave equation migration. //74th Ann. Internat Mtg., Soc. Expi. Geophys. . Expanded Abstracts.
|
Xie X B, Jin S W, Wu R S. 2006. Wave-equation-based seismic illumination analysis. Geophysics, 71(5): S169-S177. DOI:10.1190/1.2227619 |
Xie X B, Wu R S, Huang L J, et al. 2005. Seismic resolution and illumination: A wave-equation-based analysis. //75th Ann. Internat Mtg., Soc. Expi. Geophys. . Expanded Abstracts.
|
Xu L, Zhang Y S, Li B, et al. 2019. Research and application of reverse time migration using wavefield decomposition imaging condition. Progress in Geophysics (in Chinese), 34(4): 1541-1547. DOI:10.6038/pg2019CC0260 |
Yan R, Guan H M, Xie X B, et al. 2014. Acquisition aperture correction in the angle domain toward true-reflection reverse time migration. Geophysics, 79(6): S241-S250. DOI:10.1190/geo2013-0324.1 |
Yan Z, Xie X B. 2016. Full-wave seismic illumination and resolution analyses: a Poynting-vector-based method. Geophysics, 81(6): S447-S458. DOI:10.1190/geo2016-0003.1 |
Yang Q Q, Zhang S L. 2008. Least-squares Fourier finite-difference migration. Progress in Geophysics (in Chinese), 23(2): 433-437. |
Yoon K, Marfurt K J. 2006. Reverse-time migration using the Poynting vector. Exploration Geophysics, 37(1): 102-107. DOI:10.1071/EG06102 |
Zhang P, Mao W J. 2018. Weighted structure-enhancing constrained least-squares reverse time migration of simultaneous-source data. Chinese Journal of Geophysics (in Chinese), 61(10): 4088-4099. DOI:10.6038/cjg2018L0251 |
Zhang Y, Sun J. 2009. Practical issues in reverse time migration: true amplitude gathers, noise removal and harmonic source encoding. First Break, 27(1): 29-35. DOI:10.3997/1365-2397.2009002 |
Zhang Y, Zhang G Q, Bleistein N. 2003. True amplitude wave equation migration arising from true amplitude one-way wave equations. Inverse Problems, 19(5): 1113-1138. DOI:10.1088/0266-5611/19/5/307 |
Zhou D H, Qu Y M, Li Z C, et al. 2020. Multisource least-squares reverse time migration for irregular surface. Progress in Geophysics (in Chinese), 35(4): 1528-1538. DOI:10.6038/pg2020DD0217 |
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 |
Zhou Y, Wang H Z. 2017. Ray+Kirchhoff inversion based RTM imaging condition with implicit directional wavefield decomposition. //87th Ann. Internat Mtg., Soc. Expi. Geophys. . Expanded Abstracts.
|
陈生昌, 马在田, Wu R S. 2007. 波动方程偏移成像阴影的照明补偿. 地球物理学报, 50(3): 844-850. DOI:10.3321/j.issn:0001-5733.2007.03.025 |
方修政, 钮凤林, 吴迪. 2018. 基于逆散射成像条件的最小二乘逆时偏移. 地球物理学报, 61(9): 3770-3782. DOI:10.6038/cjg2018L0721 |
黄建平, 李振春, 孔雪, 等. 2013. 碳酸盐岩裂缝型储层最小二乘偏移成像方法研究. 地球物理学报, 56(5): 1716-1725. DOI:10.6038/cjg20130529 |
李振春, 黄金强, 黄建平, 等. 2017a. 基于平面波加速的VTI介质最小二乘逆时偏移. 地球物理学报, 60(1): 240-257. DOI:10.6038/cjg20170120 |
李振春, 周丽颖, 黄建平, 等. 2017b. 最小二乘逆时偏移在碳酸盐岩缝洞型储层成像中的应用. 地球物理学进展, 32(2): 664-671. DOI:10.6038/pg20170229 |
刘红伟, 李博, 刘洪, 等. 2010. 地震叠前逆时偏移高阶有限差分算法及GPU实现. 地球物理学报, 53(7): 1725-1733. DOI:10.3969/j.issn.0001-5733.2010.07.024 |
刘学建, 刘伊克. 2016. 表面多次波最小二乘逆时偏移成像. 地球物理学报, 59(9): 3354-3365. DOI:10.6038/cjg20160919 |
沈雄君, 刘能超. 2012. 裂步法最小二乘偏移. 地球物理学进展, 27(2): 761-770. DOI:10.6038/j.issn.1004-2903.2012.02.044 |
王彦飞, 杨长春, 段秋梁. 2009. 地震偏移反演成像的迭代正则化方法研究. 地球物理学报, 52(6): 1615-1624. DOI:10.3969/j.issn.0001-5733.2009.06.024 |
许璐, 张永升, 李博, 等. 2019. 逆时偏移波场分解成像条件研究及应用. 地球物理学进展, 34(4): 1541-1547. DOI:10.6038/pg2019CC0260 |
杨其强, 张叔伦. 2008. 最小二乘傅立叶有限差分偏移. 地球物理学进展, 23(2): 433-437. |
张攀, 毛伟建. 2018. 同时震源数据的加权结构增强最小二乘逆时偏移. 地球物理学报, 61(10): 4088-4099. DOI:10.6038/cjg2018L0251 |
周东红, 曲英铭, 李振春, 等. 2020. 起伏地表条件下的混叠数据最小二乘逆时偏移. 地球物理学进展, 35(4): 1528-1538. DOI:10.6038/pg2020DD0217 |
周华敏, 陈生昌, 任浩然, 等. 2014. 基于照明补偿的单程波最小二乘偏移. 地球物理学报, 57(8): 2644-2655. DOI:10.6038/cjg20140823 |