一种基于逆时偏移的高效照明补偿成像方法

周阳， 曹俊兴， 王兴建， 胡江涛， 王华忠， 刘文卿

1 油气藏地质及开发工程国家重点实验室， 成都理工大学, 成都 610059 2 同济大学海洋与地球科学学院， 上海 200092 3 中国石油勘探开发研究院西北分院， 兰州 730020

0 引言

通过求解双向波方程，逆时偏移(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.这表明对于固定偏移距的观测系统，速度异常体下方更加难以被地震波场照明，需要更多的照明补偿才能够消除照明不均衡对成像结果的影响.据以上分析，本文定义二维情况下衡量成像点处照明规则程度的雅克比行列式为：

图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

(1)

其中ψ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的具体表达形式为：

Wx=4IsIr,

(5)

其中：

(6)

式(6)中A(xr,x)为高频近似下波场经过成像点x反射被位于xr检波点接收到的振幅，c0(xs)为检波点xr处的背景速度.

1.2 照明补偿算子的RTM实现

(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)

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.

2.2 实际资料测试

为了检验本文照明补偿成像方法在实际资料上的应用效果，利用西部某二维实际资料进行测试.图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.

3 讨论

本文给出的照明补偿算子为对角Hessian逆算子的一个近似，在高频近似意义下矫正了成像点处照明不均匀的影响，提高了成像精度.由于本文给出的照明补偿算子并非完整Hessian算子的逆，因此在利用本文的照明补偿算子进行照明补偿成像后，会有残余的波场传播效应有待进一步消除.对由于照明补偿算子近似Hessian逆算子不精确造成的残留波场传播效应问题，可将本文给出照明补偿算子融入到迭代求解的最小二乘类成像框架中，利用其迭代最小二乘成像方法中构造的预条件梯度项，在后续的迭代过程中消除残余的波场传播效应，进一步提高成像精度.Burgess和Warner(2015)利用Plessix和Mulder(2004)给出的照明补偿算子构造预条件梯度项并应用到全波形反演(FWI)中，成功加速了FWI收敛效率.类似地，从数值实验对比可以看出，利用本文给出照明补偿算子构造预条件梯度项也有望提高最小二乘类成像方法的收敛速度，这是本文后续研究方向之一.

4 结论

本文给出了一种新的RTM照明补偿方法.该照明补偿方法利用地下角度域相对于地表坐标的扰动量来刻画局部成像点的照明均匀程度，据此定义照明补偿算子进行成像过程中的照明补偿.通过高频近似下的波场边界积分表达及对介质的合理近似，给出了照明补偿算子基于RTM的高效实现方式.模型和实际资料的测试对比表明，本文提出的照明补偿算子能够在较大程度上提升RTM成像质量.