0 引言

双程波动方程可以模拟出任意复杂介质的地震波传播形态，因此，相对于其他偏移方法，基于双程波动方程的叠前深度逆时偏移方法——逆时偏移(RTM)，它对于陡倾角、速度变化剧烈的复杂地质构造具有更强的振幅保真性和更高的成像质量.双程波成像方法理论早在20世纪80年代就被一些地球物理学家提出(Baysal et al., 1983; Whitmore, 1983)，但受制于当时有限的计算能力，逆时偏移没有在工业生产中得到广泛应用.随着计算机硬件技术的发展以及波场数字模拟算法的显著提高(Zhang and Yao, 2013; Wang et al., 2014, 2017; Yao et al., 2016)，逆时偏移逐步突破了制约其发展的瓶颈，成为工业生产中最重要的成像技术之一.

虽然逆时偏移方法具有高精度、高保真的成像能力，但是基于目前广泛使用的互相关成像条件(Claerbout, 1971)会产生干扰成像质量的低频噪音，特别是当速度模型中存在较强的速度梯度异常值或者较强的波阻抗反射界面时，这是因为正演模拟的震源波场和逆时反传的检波点波场都会在这些波阻抗反射界面上产生较强的逆散射能量，这些逆散射能量满足互相关成像条件，因此，成像过程会产生较强的低频噪音，从而干扰成像质量(Yoon et al., 2004; Liu et al., 2011; 杜启振等, 2013).为了消除这种低频噪音，Baysal等(1984)使用恒定波阻抗代入波动方程，得到了垂直入射方向无反射的双程波动方程，以此消除波阻抗界面达到压制逆散射能量的目的.Loewenthal等(1987)提出了直接平滑模型的慢度来压制界面强反射能量方法，其给出的数值模拟结果显示平滑慢度可以有效压制逆散射能量.Fletcher等(2006)使用双程波动方程，在无反射双程波动方程的基础上(Baysal et al., 1984)，引入一个方向阻尼项作为边界条件，克服了无反射双程波动方程在非垂直入射反射界面仍然会背向反射波场的缺点，这种方法的缺点是处理复杂模型时需要人工选择，处理难度较大.Yoon和Marfurt(2006)应用能流密度矢量构造大角度衰减因子来压制大角度成像，从而压制大角度成像的低频噪音.Zhang和Sun(2009)分析了Laplace滤波的物理意义，并建议偏移前先进行Laplace振幅校正，再对偏移后剖面进行Laplace滤波，以压制低频噪音.Liu等(2007, 2011)提出了对双程波动方程数值模拟的波场进行相应的单程波分量分解，选择合适的上下行波分量进行成像可以有效地压制低频噪音. Op′t Root等(2012)提出了逆散射成像条件，可以显著地压制低频噪音.Whitmore和Crawley(2012)、Brandsberg-Dahl等(2013)成功地应用逆散射成像条件去除了逆时偏移算法中的低频噪音，提高了成像的质量.

随着勘探目标越来越复杂，对于成像结果的精确性、分辨率和保幅性也提出了更高的要求.常规的一些偏移方法，例如逆时偏移、Kirchhoff偏移等，虽然在地震数据处理过程中是主要的成像方法，但是这些偏移方法都是应用地震波正演算子的伴随算子而非真正意义上逆算子来成像.采用伴随算子近似正演算子的逆算子，其精确性通常会受到地震数据噪音、不规则采样、有限偏移距等影响(Claerbout, 1992).最小二乘偏移(Least-Squares Migration，LSM)通过不断地迭代修正成像模型来拟合实际观测的地震记录，此过程隐含考虑了Hessian矩阵对振幅的影响，可以减少常规偏移过程中产生的假象，提高成像分辨率以及振幅精度等(Schuster, 1993; Nemeth et al., 1999).Tarantola(1984a, 1984b)在理论上推导了以波场残差的二范数平方为目标函数的最小二乘偏移算法，通过局部梯度算法来最小化目标函数，不断地修正反演成像剖面，并证明了常规的偏移结果只是最小二乘偏移的第一步迭代.Nemeth等(1999)对不完整数据进行了最小二乘Kirchhoff偏移并取得了很好的效果，进一步展示了最小二乘偏移的优势.目前，最小二乘逆时偏移(Least-Squares Reverse-Time Migration，LSRTM)方法的研究也取得了很大进展(Wong et al., 2011; Dai et al., 2012; Dai and Schuster, 2013; Wu et al., 2016; Yao and Jakubowicz, 2016; Yang et al., 2016; Yao and Wu, 2017; Yao et al., 2017).但是，基于伴随状态法求取的最小二乘逆时偏移的梯度是通过互相关成像条件计算的，因此，当地质体存在强的反射界面或者地震记录中含有大量折射波时，最小二乘逆时偏移的梯度中会产生很强的低频噪音，在迭代过程中会干扰成像精度、降低收敛速度.

为了压制这种低频噪音对最小二乘逆时偏移的影响，我们研究了基于线性逆散射成像条件的最小二乘逆时偏移方法.该方法用线性Born逆散射成像条件代替了常规的互相关成像条件，从而有效地压制了低频噪音对最小二乘逆时偏移成像质量的干扰，加快反演的收敛速度，提高反演的效率，使逆散射成像条件从逆时偏移扩展到最小二乘逆时偏移.实际算例验证了这种基于逆散射成像条件的最小二乘逆时偏移方法的有效性.

1 逆散射成像条件的最小二乘逆时偏移

频域中常密度声波方程可以写为

(1)

式中，v(x)为速度模型，ω是圆周频率，P(x; x_s, ω)为频率压力波场，S(ω)为震源函数.根据声波方程可以知道，速度模型的扰动会导致波场的扰动，令

(2)

相应波场为

(3)

式中，v₀(x)为背景参考速度模型，δv(x)为速度模型的扰动，P₀(x; x_s, ω)为背景速度模型产生的波场，δP(x; x_s, ω)为扰动的速度产生的扰动波场.

把公式(2)和公式(3)带入公式(1)得到

(4)

对上式的速度项进行泰勒级数展开并取一阶近似：

(5)

基于线性Born近似，结合公式(5)对公式(4)进行化简可以得到

(6)

式中，ΔS(x, ω)为散射震源，满足如下表达式：

(7)

式中，m(x)为本文偏移成像模型(速度相对扰动)，本文均使用振幅表述其数值大小，表达式如下：

(8)

因为格林函数可以表示为以下方程的解，

(9)

所以公式(6)中扰动波场可以用积分形式表示如下：

(10)

对上式进行展开，可以得到基于声波方程线性Born近似的反射地震数据的正问题(也叫反偏移方程)，可以表示如下：

(11)

式中，x_r为检波点位置、x_s为震源位置，d(x_r, ω)为去除直达波的观测地震记录. S(ω)为震源子波，G₀(x_r; x, ω)为从散射点x传播到检波点x_r的格林函数，G₀(x; x_s, ω)为从震源x_s传播到散射点x的格林函数.Born近似对微小扰动是精确的，因此背景速度必须要求是相对准确的.为了使反偏移方程式(11)以及后续表达式尽量简洁，可以把(11)式简写为

(12)

式中，为线性Born近似正演算子，通常也称为反偏移算子.地震数据偏移算子可以看成是地震数据正演算子的伴随算子(Claerbout, 1992)，逆时偏移算子也可以看成是线性Born近似正演算子的伴随算子.因此逆时偏移成像可以表示成

(13)

式中，L^*为偏移算子，是线性Born正演算子L的伴随算子.上式是常规的偏移，具体形式如下：

(14)

式中，p(x)=G₀(x; x_s, ω)S(ω)为正传震源波场，q(x)=G₀^*(x; x_r, ω)d(x_r, ω)为逆时反传地震记录得到的波场，表示对各地震道和每次单炮正演进行累加.

上式是基于互相关成像条件的成像表达式，会受到低频噪音的干扰.Op′t Root等(2012)提出的线性逆散射成像条件可以显著压制低频噪音.为了结合线性Born近似假设条件和线性逆散射条件的优势，我们研究了基于线性逆散射成像条件的最小二乘逆时偏移方法.这里定义线性Born散射算子的伴随算子为A^*，因此基于逆散射成像条件的逆时偏移为m=A^*d.在频域其具体表达式如下：

(15)

式中：∇为关于空间的梯度算子，d(x_r, ω)的定义与公式(11)相同，p(x)和q(x)的定义与公式(14)相同. B(x)可以理解为衰减逆散射能量的权系数矩阵，当B(x)为每一炮的权系数矩阵时，其和散射角度密切相关，每一炮、不同散射角的权系数矩阵都不一样；若B(x)为多炮叠加后的权系数矩阵时，其可以使用同一个权系数矩阵表示；若对各向同性介质且多炮叠加后加权，B(x)可以使用数值1近似.

常规偏移成像方法一般使用正演算子的伴随算子(计算量小，稳定性好等优点)代替其逆算子，对构造进行偏移成像.使用常规偏移成像方法对不规则采样或含噪音的地震数据进行成像会产生偏移假象，降低成像精度.最小二乘逆时偏移相对于常规逆时偏移具有更好的成像精度和保幅性.基于最小二乘反演的成像方法通过扰动模型参数，使反偏移模拟的地震数据和观测的地震数据在二范数意义下残差最小.其目标函数可以定义为

(16)

式中，d_syn为反偏移预测数据，d_obs为实际记录的地震记录.基于线性Born近似的最小二乘逆时偏移方法，反偏移数据和模型扰动之间是线性关系，可以采用共轭梯度法来实现最小化目标函数.其反演的更新方向为梯度的共轭方向，d_k，其表达式为

(17)

式中，下标k为迭代次数，β_k为共轭梯度系数，不同共轭梯度算法β_k选取也不同.目标函数关于模型参数的梯度可以通过伴随方法求取(Plessix，2006)，其表达式如下：

(18)

式中，q_k(x)=G₀^*(x; x_r, ω)(d_syn(x_r; x_s, ω)-d_obs(x_r; x_s, ω))_k为第k次迭代反传的残差波场，p(x)=G₀(x; x_s, ω)S(ω)为正传震源波场.对比公式(14)和(18)可以发现：常规逆时偏移和最小二乘逆时偏移都应用了互相关成像条件，但是前者反传的是地震记录，而后者反传的是地震记录的残差.当地质模型中存在强反射界面或者偏移数据中包含折射波时，通过公式(18)计算目标函数梯度含有很强的低频噪音，从而影响反演的收敛速度和成像质量.为了压制这种低频噪音，我们引入逆散射成像条件来代替上式中的互相关成像条件.因此，基于逆散射成像条件的目标函数梯度可以表示如下：

(19)

式中，p(x)与q_k(x)分别为正传震源波场与第k次迭代反传残差波场，均和公式(18)定义相同.对比公式(15)和(19)可以发现：基于逆散射成像条件的常规逆时偏移方法反传的是地震记录，而基于逆散射成像条件求取目标函数梯度反传的是地震记录残差.

根据公式(19)计算目标函数的梯度，根据公式(17)得到目标函数的更新方向.通过以下迭代公式对偏移模型进行更新：

(20)

式中，a_k为迭代步长.为了求取迭代步长，把迭代公式(20)带入目标函数公式(16)，利用线性Born正演算子的性质，可以得到

(21)

上式对a_k求导, 并令导数为零, 可以得到迭代步长计算公式：

(22)

2 GPU切片存储策略及计算效率

逆散射成像条件类似于互相关成像条件，为了降低存储量，我们采用边界波场重建震源波场存储策略(王保利等，2012；方修政等，2015；Ke et al., 2018).CPU/GPU协同计算方法，可以提高有限差分数值模拟计算效率(Shi and Wang, 2016)，考虑到GPU显存以及共享存储器空间内存有限的特点，我们采用分时间切片重建震源波场策略.在整个采样时间轴划分出一定数目的时间切片，例如：NC个Checkpoint点，保存这些时间切片点处(NC个Checkpoint点处)的模拟计算区域波场到CPU端或者GPU端.利用保存的时间切片点处模拟计算区域波场，可以正演模拟出相邻两个时间切片之间任意时刻的震源波场(NB个时间离散采样点震源波场)，只保存这些震源波场的边界波场到GPU端，利用保存到GPU端的边界震源波场和时间切片点处的模拟计算区域波场，可以逆时重建出相邻两个时间切片之间任意时刻的震源波场，以此类推可以重建出整个历史时刻震源波场.由于采用的是震源波场重建策略，因此可以避免存储所有历史时刻的震源波场，降低存储压力.GPU端只保存边界震源波场，可以显著缓解显存压力，避免了GPU端和CPU端数据频繁交换.分时间切片存储策略示意图如图 1所示.这种分时间切片存储策略优点是可以根据GPU显存大小，灵活划分每两个时间切片之间的时间点数NB(需要存储的边界震源波场离散点数目)，可以利用有限的GPU显存正演比较大的模型.由于可以根据GPU显存大小，人为地设定每两个时间切片之间的边界震源波场存储数量，因此分时间切片重建震源策略可以推广到三维.

和互相关成像条件对比：逆散射成像条件需要额外计算震源波场和检波点波场的时间导数和空间梯度(理论上：对程序优化可以避免时间导数求取和空间梯度直接计算求取，减少计算耗时；本文算法直接计算求取时间导数和空间梯度，对比了算法耗时).本文逆散射成像条件最小二乘逆时偏移算法，每次迭代需要四次正演，分别为边界波场重建震源波场、逆时反传残差地震记录波场、反偏移模拟共轭梯度波场(两次正演)，因此，算法最耗时部分是有限差分正演数值模拟过程.相对于正演数值模拟耗时，波场偏导数和时间导数的计算耗时是可以接受的.分时间切片重建震源波场策略，可以根据GPU的显存大小灵活地划分时间切片点数，最优化地利用CPU/GPU协同并行加速技术实现波场模拟，避免了CPU与GPU频繁的数据交换耗时.表 1给出了单炮最小二乘逆时偏移，不同模型尺寸和不同时间采样点的平均耗时对比.本文通过分时间切片存储策略，逆时重建出震源波场，运用逆散射成像条件得到目标泛函的下降方向，反偏移模拟共轭梯度波场，进一步得到迭代步长，不断迭代更新偏移成像模型剖面，直到满足终止条件，输出成像结果.基于线性逆散射成像条件的最小二乘逆时偏移流程如图 2所示.

3 模型试算

为了验证逆散射成像条件最小二乘逆时偏移算法的有效性，本文对不同的速度模型进行测试，整个算法测试过程中，没用采用任何滤波算法处理.为了提高计算效率，算法均采用CPU/GPU协同并行加速计算方法.

3.1 算法有效性及精度测试

为了验证算法的有效性和测试算法的精度，我们分别对层状介质和Sigsbee截取的部分速度模型进行了测试.

层状介质模型参数如下：横向网格点数200，纵向网格点数160，网格间距10 m.计算参数：采用主频25 Hz的Ricker子波作为震源，时间采样步长1 ms，总记录时长1 s.地表采用200检波器全接收方式，地表 1000 m处放置1炮，采用空间12阶精度、时间2阶精度有限差分进行波场数值模拟，完全匹配层(PML)吸收边界条件(Komatitsch and Tromp, 2003).图 3a为层状速度模型，图 3b为对真实层状速度平滑处理的偏移速度模型.图 3c是单炮逆散射成像条件RTM成像结果，图 3d是50次迭代单炮逆散射成像条件LSRTM成像结果，对比可以发现LSRTM可以提高深层的振幅信息，振幅均衡性更好.分别取图 3c和图 3d地表 1000 m处垂向单道振幅曲线，归一化后进行对比，如图 4a所示.图 4a为逆散射成像条件RTM、LSRTM和理论偏移成像模型(速度相对扰动模型)地表 1000 m处垂向单道曲线对比图.对比可以发现，逆散射成像条件LSRTM相对于RTM振幅均衡性更好，更加接近理论偏移成像模型.

Sigsbee截取的部分速度模型参数如下：横向网格点数600，纵向网格点数400，网格间距10 m.计算参数：采用主频25 Hz的Ricker子波作为震源，时间采样步长1 ms，总记录时长5 s.地表采用600检波器全接收方式，地表均匀放置30炮，炮间距为200 m，采用空间12阶精度、时间2阶精度有限差分进行波场数值模拟，完全匹配层(PML)吸收边界条件.图 5a为Sigsbee部分速度模型，图 5b为对真实速度模型平滑处理的偏移速度模型.图 5c是理论偏移成像模型.图 5d为50次迭代逆散射成像条件LSRTM成像结果.分别取图 5c和图 5d地表 3000 m处和4000 m处垂向单道振幅曲线，归一化后进行对比，如图 4b和图 4c所示.图 4b为地表 3000 m处垂向单道逆散射成像条件LSRTM和理论偏移成像模型(速度相对扰动模型)对比曲线，图 4c为地表 4000 m处垂向单道逆散射成像条件LSRTM和理论偏移成像模型对比曲线.对比可以发现，逆散射成像条件LSRTM可以很好地匹配理论偏移成像模型.图 6a为图 5d的前20次迭代归一化目标函数收敛曲线.

3.2 算法抗低频噪音测试

常规LSRTM和逆散射条件LSRTM相对于RTM具有提高成像分辨率的能力.当速度模型存在强反射界面时，或者大偏移距地震记包含大量的折射波时，常规的互相关成像条件LSRTM会因为梯度中低频噪音的存在，影响收敛速度和成像质量；逆散射条件LSRTM可以有效地压制梯度中的低频噪音，提高收敛速度和成像质量.本文分别对地堑速度模型和二维盐丘(Salt2d)速度模型测试了逆散射条件LSRTM抗低频噪音的能力.

地堑状速度模型参数如下：横向网格点数340，纵向网格点数170，网格间距10 m.计算参数：采用主频25 Hz的Ricker子波作为震源，时间采样步长1 ms，总记录时长2 s.地表采用340检波器全接收方式，地表均匀放置68炮，炮间距为50 m，采用空间12阶精度、时间2阶精度有限差分进行波场数值模拟，完全匹配层(PML)吸收边界条件.图 7为地堑状速度模型测试结果.图 7a为地堑速度模型，图 7b为对真实速度平滑处理的偏移速度模型.图 7c是常规互相关成像条件RTM成像结果，有很强的低频噪音干扰.图 7d为12次迭代常规互相关成像条件LSRTM成像结果，同图 7c对比可以发现LSRTM可以压制低频噪音，提高成像质量，但是初始有限次迭代并不能很好地去除低频噪音，尤其是地堑构造倾斜反射层处的噪音，这些噪音干扰层位的有效更新，破坏偏移成像模型的振幅和相位.图 7e为逆散射成像条件RTM成像结果，可以看出逆散射成像条件能有效地压制低频噪音干扰.图 7f为12次迭代逆散射成像条件LSRTM成像结果，同图 7e对比可以看出，逆散射成像条件LSRTM可以有效地提高成像质量，压制逆时偏移假象，显著提高成像质量.对比图 7d和图 7f：逆散射成像条件LSRTM相对于互相关成像条件LSRTM，浅层区域成像结果具有相对较少的低频噪音，反射界面成像清楚，反射层位振幅均衡，验证了逆散射成像条件LSRTM的抗低频噪音能力.

Salt2d模型参数如下：模型横向网格点数649，纵向网格点数170.计算参数：采用主频25 Hz的Ricker子波作为震源，时间采样步长1 ms，总记录时长2 s.地表采用649个检波器全接收方式，地表均匀放置129炮，炮间距为50 m.在计算中，我们采用和上面例子相同的有限差分算法和边界吸收条件.本文选取了Salt2d模型进行了LSRTM测试，主要测试了互相关成像条件LSRTM和逆散射成像条件LSRTM对强反射界面复杂速度体的成像质量.图 8为Salt2d模型测试结果.图 8a为Salt2d速度模型，可以看出Salt2d速度模型中存在高速盐丘异常体，在盐丘的边界存在强反射界面.基于互相关成像条件的RTM会产生较强的低频噪音.图 8b为对真实速度模型平滑处理的偏移速度模型.图 8c是常规的逆时偏移剖面，盐丘顶部较大速度梯度区域成像结果被低频噪音严重干扰.图 8d为11次迭代常规互相关成像条件LSRTM结果，可以看出低频噪音的压制和成像质量均得到一定程度提升，但是深部振幅没得到有效更新，例如：盐下成像质量没有得到保证.图 8e为逆散射条件逆时偏移剖面，与图 8c常规互相关条件逆时偏移剖面相比较，可以看出低频噪音被有效地压制.图 8f为11次迭代逆散射成像条件LSRTM结果，可以看出振幅更加均衡.对比图 8d和图 8f：前11次迭代反演，逆散射条件LSRTM相对于互相关成像条件LSRTM，偏移成像模型得到了有效的更新，尤其盐丘边界及盐丘下的界面成像质量得到了有效提高.图 6b为图 8d和图 8f归一化目标函数收敛曲线.可以看出：逆散射成像条件LSRTM目标函数前几次迭代下降量较大，两者整体趋势基本一致.综合以上：初始相同有限次迭代，在成像质量和收敛速度方面，逆散射成像条件LSRTM比常规互相关成像条件LSRTM具有优势.

3.3 成像质量对比分析

为了分析讨论逆散射成像条件LSRTM和互相关成像条件LSRTM成像质量差别，对图 8d和图 8f局部区域进行了放大对比，如图 9所示.图 9a和9b分别为互相关成像条件LSRTM结果和逆散射成像条件LSRTM结果局部放大图，截取于岩丘左上边界(图 8d和图 8f实线白色矩形框内部分).对比岩丘上边界可以看出：逆散射成像条件LSRTM(图 9b)，岩丘上边界振幅更加均衡，连续性更好；互相关成像条件LSRTM仍有少量低频噪音干扰岩丘上边界振幅.理论上，二范数意义最小二乘优化方法经过不断迭代更新，最终低频噪音可以被消除.但是，低频噪音不会明显改变预测地震数据的振幅，最小化地震数据振幅差异的LSRTM反演对低频噪音不是很敏感，因此，初始有限次迭代主要用于去除低频噪音，而偏移成像模型没有得到有效的更新.同时剩余的低频噪音会影响岩丘上边界振幅均衡性和连续性.图 9c和9d分别为互相关成像条件LSRTM和逆散射成像条件LSRTM偏移成像局部放大图，截取于岩丘下部(图 8d和图 8f虚线白色矩形框内区域).对比图 9c和图 9d中白色实线椭圆和黑色虚线椭圆区域可以发现：初始相同有限次迭代，逆散射成像条件LSRTM成像精度较高，可以更新出椭圆区域处的构造.

为了直观对比逆散射成像条件和互相关成像条件LSRTM梯度，以及讨论低频噪音对于梯度更新和偏移成像模型更新的影响，我们对比了单炮不同成像条件第一次迭代梯度.图 10a和10b分别为Salt2d模型互相关成像条件和逆散射成像条件LSRTM第一次迭代梯度，炮点位于地表 3 km处.可以看出，互相关成像条件LSRTM，强逆散射能量成像路径含有很强的低频噪音；逆散射成像条件LSRTM，岩丘上边界区域强逆散射能量成像路径上，低频噪音被很好地压制.互相关成像条件LSRTM，初始有限次梯度迭代更新过程中，反偏移很难辨别出岩丘边界构造散射波信息，岩丘上边界构造信息被低频噪音干扰，很难得到有效更新.逆散射成像条件LSRTM，由于岩丘上边界范围的低频噪音被很好地压制，反偏移可以辨别岩丘边界构造散射波信息，因此岩丘上边界构造可以得到相对有效的更新.同样岩丘下部构造，互相关成像条件LSRTM目标泛函梯度中低频能量很强，导致其散射波能量相对较弱，初始有限次迭代，反偏移很难得到其散射波能量，因此得不到有效更新，导致岩丘下部偏移成像模型质量得不到有效提高.逆散射成像条件LSRTM，目标泛函梯度中低频噪音被很大程度地压制，低频噪音能量相对较小，反偏移得到的散射波场中，相对均衡包含了岩丘下部梯度散射波能量，因此迭代更新过程中，岩丘下部可以得到有效的更新.

4 结论

针对常规互相关成像条件的最小二乘逆时偏移方法，容易受到低频噪音影响，减慢收敛速度和降低成像质量，因此本文发展了基于逆散射成像条件的最小二乘逆时偏移.利用可以有效压制低频噪音的逆散射成像条件，代替常规互相关成像条件计算目标函数梯度.通过理论分析和数值模拟计算，我们发现，当速度体中存在强反射界面时，基于常规互相关成像条件逆时偏移方法会产生低频噪音，基于常规互相关成像条件的最下二乘逆时偏移方法收敛速度会变慢，一部分迭代主要用于去除低频，无法保证成像剖面得到有效迭代更新.与此相反，基于逆散射成像条件的最小二乘逆时偏移，可以有效压制低频噪音导致的成像剖面无效率更新，加快成像剖面有效收敛速度，提高成像质量.