0 引言

随着油气勘探、开发的深入，勘探难度逐渐加大，在复杂构造、复杂岩性、复杂地表条件等地区进行油气地震勘探面临诸多难题，其中剧烈起伏的地表已成为制约上述地区地震勘探发展的瓶颈之一^[1]。传统的偏移成像方法大多基于水平地表假设，对起伏地形的适应性不强，容易导致复杂地表条件下的地下成像质量不高^[2]。

相比于其他偏移方法，逆时偏移(RTM)应用双程波动方程延拓波场，避免了对波动方程的近似，无倾角限制，是目前较为精确的偏移成像方法，但仍属于常规偏移的范畴，其偏移算子是线性正演算子的共轭转置，而不是它的逆^[3]。当地震数据采集不足或不规则、地下构造复杂、波场带宽有限时，常规偏移方法只能对地下构造模糊成像，无法满足岩性油气藏勘探、开发的需求。最小二乘偏移(Least-Squares Migration, LSM)将成像看作最小二乘意义下的反演问题，通过不断拟合误差泛函，得到振幅保真性更好、分辨率更高、偏移噪声更少的成像结果^[4-6]。早期的LSM主要基于射线理论^[7-8]和单程波理论^[9-13]，近年来基于双程波理论的LSRTM得到广泛关注^[14-19]，但目前的LSRTM算法大多基于水平地表假设，没有考虑起伏地形的影响，在一定程度上限制了其推广应用。

正演模拟是LSRTM的组成单元，有限差分法计算效率高、占用内存小且易简单实现，是当前应用最广泛的求解双程波动方程方法，但在模拟复杂地表及地下复杂界面时，在应用笛卡尔坐标系下的规则网格剖分时，易出现阶梯状边界而形成虚假绕射波。贴体网格能很好地拟合任意起伏地形，宜于处理起伏地表问题^[20]。但贴体网格需要坐标变换，而求解变换后的曲坐标系波动方程时，标准交错网格会引入插值误差而降低模拟精度。为此，本文将全交错网格^[21]引入贴体网格中，实现了起伏地表有限差分正演模拟，进一步推导了贴体网格线性Born正演方程，在此基础上提出了基于贴体全交错网格的起伏地表LSRTM算法，较好地克服了起伏地形的影响。

此外，起伏地形条件的地震数据常含有较多噪声，常用的基于L₂模拟合的LSRTM算法对噪声非常敏感，尤其是当地震数据中含有奇异值时，常规LSRTM结果被噪声严重干扰。L₁模对强噪声的容忍性较L₂模更好，但由于L₁模在零值处不可导，因而常用Huber范数或L₁、L₂混合模替代^[22]。T分布是另一种较L₂模更稳健的方法，已被成功用于噪声数据全波形反演^[23]。Aravkin等^[24]认为，T分布相比Huber范数和混合模，在缺失数据的条件下更稳健，且没有多余参数，因而简单实用，而Huber范数和混合模的结果严重依赖参数选取，需要大量的尝试。本文将T分布推广到起伏地表LSRTM，通过处理含噪数据，证实了方法的有效性。同时考虑到LSRTM的计算量，引入动态相位编码技术，将多炮数据组合成一个超道集，明显节省了计算量、提高了计算效率。

1 贴体网格正演模拟

众所周知，RTM及LSRTM算法的核心为波场延拓，波场模拟方法的优劣直接决定LSRTM的成败，在起伏地形条件下更是如此。因此，在起伏地形条件下开展LSRTM，首先要研究起伏地形条件的正演模拟。传统有限差分法在面对起伏地形时存在困难，而贴体网格则是一种很好的解决方案，在山前带具有广阔的应用前景。

贴体网格是一种适合复杂地表介质的网格离散方法，网格生成的原则是使离散后的网格边界与地表形态吻合，以避免人为产生的阶梯边界引起的虚假散射，贴体网格可以由计算空间到物理空间的坐标变换获得(图 1)。

贴体网格生成之后，笛卡尔坐标系(x, z)与曲线坐标系(ξ, η)的网格点建立了一一对应关系，从而可求得两个坐标系间的映射关系(即雅克比矩阵)ξ_x, ξ_z, η_x, η_z，且有$ {f_x} \equiv \frac{{\partial f}}{{\partial x}}, {f_z} \equiv \frac{{\partial f}}{{\partial z}}$，f代表ξ或η。

利用链式法则和常规笛卡尔坐标系波动方程，可导出曲线坐标系二维声波波动方程^[25]

$ \left\{ {\begin{array}{*{20}{l}} {\rho \frac{{\partial {v_x}}}{{\partial t}} = {\xi _x}\frac{{\partial p}}{{\partial \xi }} + {\eta _x}\frac{{\partial p}}{{\partial \eta }}}\\ {\rho \frac{{\partial {v_z}}}{{\partial t}} = {\xi _z}\frac{{\partial p}}{{\partial \xi }} + {\eta _z}\frac{{\partial p}}{{\partial \eta }}}\\ {\frac{{\partial p}}{{\partial t}} = \rho {v^2}\left( {{\xi _x}\frac{{\partial {v_x}}}{{\partial \xi }} + {\xi _z}\frac{{\partial {v_z}}}{{\partial \xi }} + {\eta _x}\frac{{\partial {v_x}}}{{\partial \eta }} + {\eta _z}\frac{{\partial {v_z}}}{{\partial \eta }}} \right)} \end{array}} \right. $

(1)

式中: v_x和v_z为质点振动速度；p为声压；ρ为介质密度；v为声波速度。

由式(1)可知，曲坐标系下每个变量都要在同一网格点的两个方向计算空间导数，常规的标准交错网格已不满足需求(图 2a)，计算v_x的垂向导数与v_z的横向导数均需复杂的波场插值，降低了模拟精度。由于同位网格的一阶中心差分缺少耗散特性，引起高频振荡现象，因此需要人为滤波，增加了实现的复杂性、降低了模拟精度。全交错网格机制的主要思想是速度和声压的不同分量交错定义在同一网格点(图 2b)。目前全交错网格主要用于水平地表和矩形规则网格，本文将全交错网格引入曲线坐标系^[25]。由全交错网格的网格定义机制可知，其满足式(1)的交错分布特性，避免了插值和滤波运算，提高了模拟精度、降低了实现复杂性。由于同一变量定义在同一网格的两个不同位置，所以需要分别更新、计算，具体差分格式与常规标准交错网格相同，因而简单、易实现。

边界条件是正演模拟的关键，受计算机内存和计算速度等因素限制，数值模拟的范围必须是有限的，除上边界外的其他三个边界都会产生人工边界问题。为此，本文采用完全匹配层吸收边界^[26]，使波能够自由地穿过边界，而不产生反射；上边界为自由边界，本文采用牵引力镜像法实施^[20]。

2 起伏地表LSRTM基本原理 2.1 起伏地表线性化波动方程

为了推导的简便，仿照一阶速度—应力方程与二阶弹性波方程的相互转换方式，可将式(1)重新整理为

$ \begin{array}{l} {s^2}\frac{{{\partial ^2}p}}{{\partial {t^2}}} = \left( {\xi _x^2 + \xi _z^2} \right)\frac{{{\partial ^2}p}}{{\partial {\xi ^2}}} + 2\left( {{\xi _x}{\eta _x} + {\xi _z}{\eta _z}} \right)\frac{{{\partial ^2}p}}{{\partial \xi \partial \eta }} + \\ \;\;\;\;\;\;\left( {\eta _x^2 + \eta _z^2} \right)\frac{{{\partial ^2}p}}{{\partial {\eta ^2}}} + \left( {{\xi _{xx}} + {\xi _{zz}}} \right)\frac{{\partial p}}{{\partial \xi }} + \left( {{\eta _{xx}} + {\eta _{zz}}} \right)\frac{{\partial p}}{{\partial \eta }} \end{array} $

(2)

式中：s为慢度场；$ {f_{xx}} \equiv \frac{{{\partial ^2}f}}{{\partial {x^2}}}, {f_{zz}} \equiv \frac{{{\partial ^2}f}}{{\partial {z^2}}}$，f代表ξ或η。

将慢度场的平方s²分解为背景慢度场平方s₀²与扰动慢度场平方Δs²的叠加，即

$ {s^2} = s_0^2 + \Delta {s^2} $

(3)

由场的叠加原理可知，总波场p可理解为由背景介质产生的背景波场p₀和由扰动介质产生的扰动波场p_s叠加而成，即

$ p = {p_0} + {p_s} $

(4)

p₀与p都满足波动方程，分别将其代入式(2)并相减，然后应用Born近似，由p₀替代p₀+p_s，可得p_s的控制方程

$ \begin{array}{l} s_0^2\frac{{{\partial ^2}{p_s}}}{{\partial {t^2}}} = \left( {\xi _x^2 + \xi _z^2} \right)\frac{{{\partial ^2}{p_s}}}{{\partial {\xi ^2}}} + 2\left( {{\xi _x}{\eta _x} + {\xi _z}{\eta _z}} \right)\frac{{{\partial ^2}{p_s}}}{{\partial \xi \partial \eta }} + \\ \;\;\;\;\;\;\;\;\;\;\;\;\;\;\left( {\eta _x^2 + \eta _z^2} \right)\frac{{{\partial ^2}{p_s}}}{{\partial {\eta ^2}}} + \left( {{\xi _{xx}} + {\xi _{zz}}} \right)\frac{{\partial {p_s}}}{{\partial \xi }} + \\ \;\;\;\;\;\;\;\;\;\;\;\;\;\;\left( {{\eta _{xx}} + {\eta _{zz}}} \right)\frac{{\partial {p_s}}}{{\partial \eta }} - \Delta {s^2}\frac{{{\partial ^2}{p_0}}}{{\partial {t^2}}} \end{array} $

(5)

式(5)即为曲线坐标系线性化Born正演(反偏移)方程。可见，p_s是由p₀与Δs²的相互作用作为二次震源在背景介质中传播的波场，与Δs²呈线性关系，具有明确的物理含义。

定义背景介质中的格林函数G满足

$ \begin{array}{l} s_0^2\frac{{{\partial ^2}G}}{{\partial {t^2}}} = \left( {\xi _x^2 + \xi _z^2} \right)\frac{{{\partial ^2}G}}{{\partial {\xi ^2}}} + 2\left( {{\xi _x}{\eta _x} + {\xi _z}{\eta _z}} \right)\frac{{{\partial ^2}G}}{{\partial \xi \partial \eta }} + \\ \;\;\;\;\;\;\;\;\;\;\;\;\;\left( {\eta _x^2 + \eta _z^2} \right)\frac{{{\partial ^2}G}}{{\partial {\eta ^2}}} + \left( {{\xi _{xx}} + {\xi _{zz}}} \right)\frac{{\partial G}}{{\partial \xi }} + \\ \;\;\;\;\;\;\;\;\;\;\;\;\;\left( {{\eta _{xx}} + {\eta _x}} \right)\frac{{\partial G}}{{\partial \eta }} + \delta \left( {x - {x_s}} \right) \end{array} $

(6)

利用格林函数可将式(5)中的p_s表示为

$ {p_s} = - \int_\Omega {\rm{d}} x\int_t {\rm{d}} tG\Delta {s^2}\frac{{{\partial ^2}{p_0}}}{{\partial {t^2}}} = \int_\Omega {\rm{d}} x\int_t {\rm{d}} tG\frac{{{\partial ^2}{p_0}}}{{\partial {t^2}}}m $

(7)

式中：m=-Δs²，将其定义为模型参数；Ω为积分空间范围。为方便后续的推导，式(7)可写成算子的形式

$ {p_s} = \mathit{\boldsymbol{Lm}} $

(8)

式中$ \mathit{\boldsymbol{L}}{\rm{ = }}\int_\Omega {{\rm{d}}\mathit{x}} \int_t {{\rm{d}}\mathit{tG}\frac{{{\partial ^2}{p_0}}}{{\partial {t^2}}}} $为起伏地表扰动波场的线性正演算子，形式上与水平地表线性正演算子相同，不同之处在于格林函数G是由起伏地表正演方程(式(6))计算获得，考虑了起伏地形的波场传播特征。

2.2 基于T分布的起伏地表LSRTM

基于反演的成像方法寻求最优的地下介质模型，以使正演波场与观测波场残差的模最小，是一个最小范数问题。常规基于L₂模的LSRTM目标泛函为

$ {J_{{{\rm{L}}_2}}}(\mathit{\boldsymbol{m}}) = \frac{1}{2}{\left( {\mathit{\boldsymbol{Lm}} - {p_{{\rm{obs}}}}} \right)^2} $

(9)

式中p_obs为观测记录。

L₁模相比L₂模更稳健，但当数据误差接近于零时，应用L₁范数准则求取目标函数梯度会出现不稳定，因而常用Huber范数或混合模代替，但由于最后的反演结果严重依赖于参数的选取，一般需要多次尝试，大大增加了计算量^[27]。

与混合模不同的是，T分布对参数的依赖性较弱，因而更简单、高效。基于T分布的LSRTM目标泛函为

$ {J_{\rm{T}}}(\mathit{\boldsymbol{m}}) = \frac{{\alpha + 1}}{2}\lg \left[ {1 + \frac{1}{{\alpha {\sigma ^2}}}{{\left( {\mathit{\boldsymbol{Lm}} - {p_{{\rm{obs}}}}} \right)}^2}} \right] $

(10)

式中：α为自由度参数；σ为尺度参数。特别地，当α=1时，式(10)为柯西分布。大量试算表明，α=2、σ=1是一组较好的参数组合，适用于绝大多数情况。

采用梯度导引类算法(最速下降或共轭梯度)求解，需要计算目标泛函关于模型参数的梯度，L₂模和T分布目标泛函的梯度公式可分别写为

$ \frac{{\partial {J_{{{\rm{L}}_2}}}}}{{\partial \mathit{\boldsymbol{m}}}} = {\mathit{\boldsymbol{L}}^*}\left( {\mathit{\boldsymbol{Lm}} - {p_{{\rm{obs}}}}} \right) $

(11)

$ \frac{{\partial {J_{\rm{T}}}}}{{\partial \mathit{\boldsymbol{m}}}} = {\mathit{\boldsymbol{L}}^*}\frac{{\alpha + 1}}{{\alpha {\sigma ^2} + {{\left( {\mathit{\boldsymbol{Lm}} - {p_{{\rm{obs}}}}} \right)}^2}}}\left( {\mathit{\boldsymbol{Lm}} - {p_{{\rm{obs}}}}} \right) $

(12)

式中上角“*”表示共轭转置。由式(11)、式(12)可以看出，两种范数的不同之处在于：L₂模直接利用波场残差计算梯度，T分布则利用加权后的残差记录计算梯度，其他更新流程完全相同。因此，基于常规L₂模的各种加速方法都可直接应用于T分布，如相位编码技术。

LSRTM的计算量过于庞大，从而限制了其推广应用。考虑到LSRTM的计算量与炮数成线性关系，因而通过相位编码技术^[28]将多个炮集组合成一个超道集，可有效减小计算量，在此基础上应用共享存储并行编程(OpenMP)技术，可进一步提高计算效率，本文选用动态震源极性编码^[29]。

3 模型试算

本文给出三个模型算例，以验证本文算法的有效性和适用性。

3.1 倾斜界面模型

以倾斜界面的网格化离散为例，说明贴体网格在刻画不规则界面时的优势。图 3为倾斜界面模型剖分示意图。由图可见：常规矩形网格不能很好地离散倾斜界面，在地表界面处产生很多阶梯状的毛刺(图 3a)；贴体网格则可很好地逼近复杂界面(图 3b)，不仅可以完全拟合真实界面，且在界面处满足很好的正交性，有利于实施边界条件。

图 4为0.3s时刻波场快照。由图可见：常规有限差分算法对倾斜界面的离散作用较差，产生较强的虚假绕射和散射，严重干扰了有效波场(图 4a)；贴体全交错网格算法很好地拟合了倾斜界面(图 4b)。

3.2 起伏地表洼陷模型

图 5为起伏地表洼陷模型，图 6为贴体网格剖分及其山峰处的局部放大图。可见，贴体网格对复杂界面的适应性较强。图 7为本文算法得到的RTM结果及其Laplace滤波结果。由图可见：首先，在本文算法得到的RTM结果中地下构造基本被低频噪声掩盖(图 7a)。其次，RTM结果的Laplace滤波结果虽然可正确成像地下构造，但仍然存在如下问题(图 7b)：①偏移噪声大。Laplace滤波不能彻底去除低频噪声，且还引入了高频噪声；②反射同相轴中间能量强、两侧能量较弱，即振幅均衡性不佳；③由于地下照明强度随深度的增大而减弱，因而RTM结果深部能量较弱，振幅保真性差。

采用本文提出的起伏地表LSRTM算法可以有效地解决常规RTM存在的问题。然而LSRTM的计算量过于庞大，目前的计算机资源无法满足运算需要。多震源技术可有效缓解计算量问题，但会引入串扰噪声，采用动态相位编码技术可很好地压制串扰噪声，在大幅降低计算量的同时，可得到与常规LSRTM算法相当的结果。将101炮地震数据利用震源极性编码方式组合成一个超道集，使计算量相当于单炮情形，从而大大缓解了计算需求。图 8为基于相位编码的起伏地表LSRTM成像结果。由迭代30次的LSRTM成像结果可见，地下构造清晰，在振幅保真性、均衡性、压制低频噪声等方面较常规RTM结果明显改善，但存在由编码引入的较强高频串扰噪声(图 8a)；迭代80次的LSRTM成像结果与理论反射率模型非常接近，有效地压制了串扰噪声(图 8b)。

图 9为残差曲线与起伏地表洼陷模型上界面中点处振幅谱。由图可见：①随着迭代次数增加，数据残差和模型残差都逐渐减小，开始下降趋势明显，随着迭代次数的进一步增大，下降趋势减慢；由于目标泛函是数据空间的拟合，因此数据残差能收敛到更低值(图 9a)。②随着迭代次数增大，高频成分逐渐得到恢复，振幅趋近于真值，表明LSRTM可以反演高频成分(图 9b)。

测试发现，三种范数对随机噪声的容忍度基本相同，噪声较弱时效果较好，噪声较强时效果较差。当数据中含有异常值时，混合模及T分布较L₂模具有较大优势。图 10为含噪单炮记录，图 11为含脉冲噪声的LSRTM结果。由图可见：L₂模LSRTM结果含有较多干扰(图 11a)；混合模(图 11b)和T分布(图 11c)LSRTM结果显著消除了脉冲噪声的影响，但前者需要多次尝试以寻求最优参数，因此T分布LSRTM更简单、实用。

3.3 起伏地表中原模型

图 12为速度模型及其反射率模型。由图可见，起伏地表中原模型不仅地表起伏剧烈，且地下构造复杂，可检验偏移算法的成像效果。采用动态相位编码技术将218炮数据组合成一个超道集，然后利用本文算法迭代计算。图 13为起伏地表中原模型LSRTM结果。由图可见：在1次迭代的LSRTM结果中存在非常强的近地表低频噪声，不能识别地下构造，且由编码引入的高频串扰也很强(图 13a)；20次迭代的LSRTM结果(图 13b)已明显压制了低频噪声，基本可识别地下构造，但仍存在较强的高频串扰，且横向振幅均衡性较差；随迭代次数进一步增大，LSRTM结果质量变好，基本消除了高频串扰和低频噪声，显著提高了振幅均衡性和保真性(图 13c、图 13d)。

4 结束语

针对起伏地表地震偏移成像存在的问题，本文发展了基于贴体全交错网格的起伏地表LSRTM算法。通过理论分析及模型试算，得到以下认识：

(1) 贴体网格能很好地拟合起伏界面，避免了矩形网格中由于阶梯离散产生的虚假绕射波，且可适应任意复杂地表情形。在曲线坐标系下采用全交错网格，避免了标准交错网格的波场插值和同位网格的高频振荡问题，提高了模拟精度，减小了实现复杂度。

(2) 与常规RTM相比，LSRTM可实现真振幅成像，提高了分辨率、减小了偏移噪声，显著提高了振幅均衡性和保真性。本文测试的三种范数在随机噪声下的表现基本相同，但相比常规L₂模，混合模与T分布能较好地压制脉冲噪声的影响，且T分布不依赖于参数选取，因此更简单、实用。

(3) 在起伏地表LSRTM算法中，加入动态相位编码技术不仅能极大地降低计算成本，且可有效压制串扰噪声，在此基础上应用OpenMP技术，可显著提高计算效率，从而将LSRTM的计算成本降低到同常规RTM相同的水平。

尚需指出，本文算法仅进行了理论模型测试，还未用于实际资料处理，且采用的模型近地表变化较为平缓，在横向变速剧烈的近地表条件下算法的适应性还需进一步探讨。此外，如何将其推广到GPU等快速计算设备上进一步提高计算效率，如何通过预条件和正则化算子等增加算法的稳定性并加快收敛速度等是下一步的研究方向。