1 引言

在油气勘探中，地震偏移是获取地下复杂构造成像的重要处理手段.然而，传统的偏移算子为正演算子的共轭算子 (Claerbout，1992)，而非逆算子.这种近似，在受有限孔径以及采集系统的不规则性等的影响下，往往会产生偏移假象 (Nemeth et al., 1999).这种影响可以通过计算Hessian矩阵的逆来进行补偿.然而在实际应用中，由于Hessian矩阵维度巨大，直接计算所需内存与成本较高，通常只近似Hessian矩阵的主对角元素 (Docherty，1991；Gray，1997；Shin et al., 2001)，其中震源照明补偿为其中一种近似 (Plessix and Mulder, 2004).然而，这种方式只能部分补偿振幅失真，而不能校正由稀疏采样、有限孔径所产生的假象.为了解决这个问题，可以通过反演的方式来考虑Hessian矩阵的影响，这种方法称为最小二乘偏移 (Schuster，1993；Nemeth et al., 1999).

最小二乘偏移最初应用于Kirchhoff偏移 (Nemeth et al., 1999；Duquet et al., 2000；王彦飞等，2009)，然后发展到单程波偏移算法中 (Wang et al., 2005；Tang，2009；周华敏等，2014)，近几年又发展到了基于双程波动方程的最小二乘逆时偏移 (LSRTM)(Dai et al., 2012；Dai and Schuster, 2013；李振春等，2014；黄建平等，2015).相对于传统的偏移算法，LSRTM具有以下优点 (Zeng et al., 2014)：振幅更加均衡、可以减弱由采集脚印或非均匀照明造成的偏移假象、可以压制旁瓣来提高空间分辨率.这种改进是通过不断迭代来逐渐纠正偏移假象，使得模拟数据与观测数据之间的误差达到最小.然而，在LSRTM梯度计算时，虽然假定相对光滑的背景速度场，但是当背景速度场变化剧烈时，仍会产生背向反射，从而产生低频噪声.虽然LSRTM可以通过迭代的方式逐渐消除低频噪声的影响，但是更好的方法是在梯度计算时直接去除.传统RTM中低频噪声的压制方法 (杜启振等，2013) 有：坡印廷矢量法 (Yoon et al., 2004)、波场分解法 (Liu et al., 2011)、角道集叠加法、带补偿的拉普拉斯滤波等.其中带补偿的拉普拉斯滤波相对于其他方法计算量小，又能较好的去除低频噪声.因此，本文将这种方法应用到LSRTM去除梯度计算中的低频噪声中，避免通过迭代的方式消除这种低频噪声的影响，从而加速LSRTM的收敛速度.

LSRTM计算过程中需要反复求解波动方程，因而算法的精度与效率主要依赖于求解波动方程的精度与效率.常用的算法是基于有限差分法，这种方法通常是在结构化的规则网格上进行差分离散，因而对地下复杂构造的描述存在阶梯近似，并且在实施过程中，由于空间采样往往由低速介质所限定，因而对于地下高速介质，则存在过采样 (Zhang，2004).有限元方法，则可以很好的描述地下复杂构造，但是这种方法存在计算量大的缺点.而格子法 (Zhang and Liu, 1999)，是一种基于非结构、非规则网格的算法，其兼具有限元精细刻画地下界面的优点，又具有较高的计算效率，并且可以结合随速度的自适应网格剖分使得整体的计算规模降低，现已成功应用到强速度对比介质 (Zhang，2004)、裂缝介质 (Zhang and Gao, 2009) 等的地震波模拟以及全波形反演中 (魏哲枫等，2014).因此，本文选用格子法作为LSRTM求解波动方程的方法，以利用其精细刻画地下界面以及随速度自适应剖分的优点.

本文首先推导了LSRTM的基本原理，并给出了求解梯度时去除低频噪声的方法，然后通过非结构化网格离散，推导出格子法具体的离散格式，最后通过模型试算验证了方法的正确性和有效性.

2 方法原理 2.1 最小二乘逆时偏移的基本原理

常密度声波方程可以表示为

(1)

其中：t为时间，x为空间坐标，x_s为震源位置，c₀(x)为背景速度场，p₀(t, x, x_s) 为震源函数f(t; x_s) 所产生的地震波场.根据扰动理论，扰动速度场δc (x)产生扰动场δp(x, t; x_s)，并满足下述方程：

(2)

其中，扰动模型m(x)定义 (Dai and Schuster, 2013) 如下：

(3)

为了计算扰动场δp(x, t; x_s)，首先需要根据方程 (1) 求解背景场p₀(t, x, x_s)，然后通过扰动模型m(x)结合方程 (2) 即可求出.

利用伴随状态法 (Plessix，2006)，对于单炮记录d(x_g, t; x_s) 进行梯度求解，可以表述为

(4)

其中，q(x, t; x_s) 为伴随场，并满足伴随方程：

(5)

为了简化公式，方程 (2) 可以使用矩阵的形式描述，公式为

(6)

其中：d为反射波场，m为扰动模型，L为反偏移算子.类似的，地震成像可以表述为

(7)

其中：m_mig为成像剖面，L^T为偏移算子.LSRTM旨在求解扰动模型m(x)，使得目标函数J (m)最小，公式为

(8)

本文采用预条件共轭梯度算法，算法描述如下：

(9)

其中：P为震源照明补偿.由上表述可知，传统的逆时偏移算法相当于最小二乘逆时偏移算法的初次迭代结果.

在LSRTM迭代求解中，虽然假定相对光滑的背景速度场，但是使用 (4) 式求解梯度时，仍存在低频噪声，尤其在浅层或强反射界面的上部，类似传统RTM，这种低频噪声可以通过带补偿的拉普拉斯滤波消除 (Zhang and Sun, 2008).在波数域，拉普拉斯滤波的形式为

(10)

其中，θ为反射角.因此，拉普拉斯算子相当于对求得的梯度应用角度域衰减因子cos²θ，但是改变了频谱与振幅.这种改变可以通过对震源子波在时间域积分两次 (刘红伟等，2010)，然后对求得的梯度进行c₀²/4标定进行补偿.

2.2 非规则网格离散

格子法的核心是建立基于积分近似的微分方程弱形式，下面以方程 (1) 为例进行非规则网格离散.

如图 1所示为格子法局部网格示意图，其中实线所包围的三角形为格子法的基本单元，虚线连接三角形形心与边的中心所构成的小闭合区域 (如虚线1-2-3…11-12-1所围区域) 为格子法的基本积分区域，类似于有限体积法的控制体积.

在k点的邻域对方程 (1) 进行面积分，则：

(11)

对上式左端应用集中质量原理，公式为

(12)

其中：Ω为基本积分区域，M_k是节点k周围每个网格单元中之值总和的三分之一 (被虚线分割出来的部分)，而为节点k处背景场p₀的二阶时间导数.

对方程 (11) 右端第一项运用高斯散度定理，可得：

(13)

其中：n为节点k周围三角形单元的个数，S_kl为节点k周围第l个三角形单元中的虚线段.

对于格子法的典型单元，如三角形ijk，通过对各个节点运用一阶插值函数，并计算一阶空间偏导数，这里引入空间偏导数算子D_x、D_z，分别代表x方向、z方向上的导数，则：

(14)

其中：A代表三角形单元ijk的面积，(p₀)_i下标处的i表示与节点i相关的物理量，(p₀)_j、(p₀)_k类似，定义b_i=(z_k-z_j)/2，a_i=(x_k-x_j)/2为几何参数，并且按照i→j→k→i方向替换即可得到b_j、a_j、b_k、a_k.

从式 (14) 可以得出，D_xp₀、D_zp₀在每个三角形内为常数，因此，方程 (13) 可以简化为

(15)

其中：(a_k)_l、(b_k)_l为节点k周围第l个单元中对应节点k的几何参数，(D_xp₀)_l、(D_zp₀)_l是p₀在节点k周围第l个单元中的空间偏导数.把 (12)、(15) 式带入 (11) 式得：

(16)

上式中对于节点k处背景场的二阶时间偏导数(∂²p₀/∂t²)_k可以采用二阶中心差分求解公式为

(17)

上式即为方程 (1) 的非规则网格离散形式，其余方程的离散格式类似.

3 模型试算

为了验证方法的正确性，分别以倾斜层模型与Marmousi模型进行测试.其中，倾斜模型主要测试格子法能够避免阶梯近似引入的虚假绕射、拉普拉斯滤波的应用效果，以及验证格子法LSRTM的正确性；通过对Marmousi模型测试，主要测试格子法能够结合随速度相关的自适应剖分来降低整体计算规模；格子法LSRTM提高空间分辨率、补偿深部照明不均匀性等的能力，以及拉普拉斯滤波对收敛速度的影响.

3.1 倾斜层模型数值测试

所用倾斜层速度模型如图 2所示，模型大小为：水平方向为3 km，深度方向为2 km.第一层速度为2000 m·s^-1，第二层速度为2800 m·s^-1，速度对比强烈.总共激发3炮，炮间距为500 m，第一炮的水平位置为1500 m.采用全覆盖接收，共301个检波器，道间距为10 m，均匀分布在地表面.选用25 Hz的雷克子波作为震源子波进行模拟，模拟时长1.5 s.

图 3所示为第二炮750 ms波场快照，由于采用非规则网格离散，避免了规则网格离散中阶梯近似所引入的虚假绕射.

偏移所用速度模型为原始速度进行平滑以后的模型，为了显示带补偿的拉普拉斯滤波对低频噪声的消除，首先对单炮结果进行测试，其中第二炮第一次迭代后的结果如图 4所示，虽然所用速度为平滑以后的速度，在图 4a中仍能见到一些低频噪声，采用带补偿的拉普拉斯滤波以后，低频噪声得到明显消除，并且滤波以后的相位与原始成像结果的相位一致.

叠加结果如图 5所示，相对于传统的RTM，LSRTM通过压制旁瓣，使得成像的分辨率明显提高.

在梯度计算中，往往会在强反射界面的上部产生很强的低频噪声，如果低频噪声在计算过程中不加去除，则会使得迭代不断更新由低频噪声对炮记录的影响，因而即使经过多次迭代以后，误差函数仍处于一个相对高值，正如图 6所示，即使经过10次迭代，绿色曲线的误差函数只与红色曲线的第二次迭代结果的误差函数相当.

3.2 Marmousi模型数值测试

为了测试本文方法对复杂模型的应用效果，对Marmousi模型进行测试.速度模型如图 7所示，模型大小为水平方向9.2 km，深度3 km.震源选用20 Hz的雷克子波，炮点埋深4 m，第一炮的水平位置为1.2 km，炮间距为80 m，总共90炮.每炮至多540道，放置于地表面，其中道间距为8 m.

在求解波动时需要满足频散关系，对于已知的震源主频，网格尺度往往由局部速度决定.然而对于有限差分，由于采用结构化、规则网格离散，网格大小一致，因此网格尺度受限于最低速度，而非局部速度，从而使得计算规模增大，尤其是速度模型中存在高速度对比的情况.而非结构化、非规则网格可以利用与速度相关的自适应剖分技术，如图 8所示，对于浅层低速介质，剖分出的网格小，而对于深层高速介质，剖分出的网格大，从而可以使得整体的计算规模降低.

偏移所用速度如图 9所示，虽然使用的是平滑以后的速度模型，在成像剖面上，仍会存在大量的低频噪声，这种低频噪声掩盖了真实的成像结果 (图 10所示).如果这种低频噪声不在成像的过程中去除，则会影响目标函数的收敛速度.正如图 14绿色曲线所示，在第一次迭代以后，收敛速度变得非常慢.

图 11为LSRTM求解梯度时采用拉普拉斯滤波并进行补偿后的结果，其中低频噪声得到较好的去除，浅层以及深层构造已经显现出来，相位也与图 10相一致，从而证明带补偿的拉普拉斯滤波的有效性.

对比图 11与图 12可知，由于受不均匀照明的影响，LSRTM第一次迭代结果浅层振幅较强，深部振幅较弱，尤其是图中黑色箭头所指部位.经过十二次迭代以后 (图 12所示)，深部能量得到更好的补偿，同相轴的连续性更好，并且深部振幅与浅层振幅相比也更加均衡，剖面的空间分辨更高.

为了展示LSRTM具有提高空间分辨率的作用，对图 11、12中白色虚线所围区域进行放大显示，如图 13a和13b所示，箭头所指黑色旁瓣部分，经过迭代，得了很好的压制，使得LSRTM结果具有更高的空间分辨率.

对比图 14所示收敛曲线，在计算梯度时，经过滤波及补偿以后，由于消除了低频噪声的影响，使得迭代更多的校正由于有限孔径、振幅不准确所带来的影响，从而使得收敛曲线相对较快，残差更小.

4 结论

本文发展了基于非结构化网格的最小二乘逆时偏移方法，并把拉普拉斯滤波以及补偿应用于最小二乘逆时偏移中，在算法实现的基础上，对简单倾斜层和复杂的marmousi模型进行了测试，得到了如下结论：通过非结构化网格离散，更好的利用与速度自适应剖分的优点；与传统RTM类似，在求解LSRTM梯度中，虽然使用相对光滑的偏移速度模型，在浅层以及强反射界面的上部，仍会产生低频噪声，这种低频噪声，可以通过拉普拉斯滤波及补偿进行消除，否者会影响目标函数的收敛速度；相比传统的RTM，LSRTM能够补偿由有限孔径、不规则采集等所引起的照明不均匀性，能够使得浅层与深层的振幅更加均衡，能够通过压制旁瓣来提高空间分辨率.