0 引言

弹性波动理论的发展以及长期的勘探实践均已表明，地震波是以矢量波的形式在地球介质中传播的。建立在声波介质假设上的问题只能主要围绕纵波速度和密度进行，忽略了横波，不符合精细勘探的要求^[1-3]。为此针对精细岩性油气藏勘探，发展弹性介质假设下的成像方法迫在眉睫。逆时偏移是当前偏移方法中较为精确的深度偏移方法且发展迅速，已经从叠后发展到叠前，从二维到三维，从声波到弹性波，从各向同性到各向异性^[4-7]。弹性波逆时偏移主要分为2个过程，首先构建震源波场和检波波场，利用成像条件得到成像结果。1986年Chang等^[8]首次完成了各向同性介质弹性波叠前逆时偏移；随后Chang等^[9-10]又深入研究，完成了2D和3D弹性波逆时偏移；董良国等^[11]分析不同地形起伏情况下自由边界的情况，模拟出地表起伏情况下弹性波复杂的传播现象，实现了弹性波起伏地表自由边界条件的数值计算；为了深入研究弹性波波场传播规律，陈可洋^[12]将波印廷矢量应用于多波多分量各向异性介质弹性波波动方程中，实现了上行波、下行波、左行波和右行波的波场分离。针对弹性波成像的串扰问题，王维红等^[13]依据各向同性介质一阶速度-应力方程组，利用Helmholtz分解提取纯纵波和纯横波波场，使用震源归一化的互相关成像条件获得纯波成像，避免了纵横波串扰问题。张伟等^[14]利用高阶交错网格有限差分数值方法构建了矢量的纵横波波场，将震源归一化的内积成像条件应用于分离后的纯纵波和横波矢量场，由此得到的转换波成像避免了传统弹性波成像方法中出现的极性反转。成像条件是影响弹性逆时偏移的另一个重要因素。杜启振等^[15]基于激发时间成像条件实现了横向各向同性介质中的多波多分量叠前逆时偏移；张晓语等^[16]基于能量守恒定理及能量密度，实现了弹性波能量互相关成像，能够压制背向散射；张智等^[17]提出稳定的激发振幅成像条件，在震源波场和检波器波场的传播过程中，保存最大能量密度的时刻和相应的波场值，归一化后获得角度依赖的反射系数成像剖面，与归一化成像条件相比，稳定激发振幅成像条件具有更小的计算量。

常规的逆时偏移使用的是反偏移算子的共轭，而不是它的逆，因此造成成像结果模糊化，成像精度难以满足精细油气勘探的精度要求^[18-19]。为了弥补常规逆时偏移的缺陷，从最小二乘思想的角度重新审视逆时偏移的过程^[20-25]。近年来，对最小二乘逆时偏移(LSRTM)的研究逐渐引起了国内外学者的关注。其中，为了提升效率，减少耗时，Dai等^[26-27]引入了相位编码；针对高精度成像这一目标，黄建平等^[28]优化了LSRTM算法；LSRTM还可应用于黏声介质，李振春等^[29]完成了黏声介质LSRTM算法，改善成像质量。考虑到LSRTM的优势，国内外学者将弹性波逆时偏移(ERTM)与最小二乘偏移(LSM)结合发展出弹性波最小二乘逆时偏移(EL‐ SRTM)来进行复杂构造的保幅成像处理^[30-31]。

从弹性波波动方程出发，推导弹性波线性正演以及逆时偏移算子表达式；然后从误差函数出发，推导Hessian矩阵的表达式，以此作为梯度预条件算子，实现一种预条件弹性波最小二乘逆时偏移(PELSRTM)算法；最后，通过SEG/EAGE二维盐丘模型测试P-ELSRTM对复杂模型成像时的有效性及其优势，并对比和分析常规ELSRTM及P-ELSRTM对低信噪比炮数据的适应性，以期证明预处理对反演过程中的改进效果。

1 弹性波最小二乘逆时偏移 1.1 线性化反偏移和偏移算子

在非均匀各向同性弹性介质中，二维弹性波一阶速度-应力方程^[30-31]可表示为

$ \left\{ {\begin{array}{*{20}{l}} {\rho \frac{{\partial {\kern 1pt} {u_{\rm{x}}}}}{{\partial {\kern 1pt} t}} - \left( {\frac{{\partial {\kern 1pt} {\sigma _{{\rm{xx}}}}}}{{\partial {\kern 1pt} x}} + \frac{{\partial {\kern 1pt} {\sigma _{{\rm{xz}}}}}}{{\partial {\kern 1pt} z}}} \right) = 0}\\ {\rho \frac{{\partial {\kern 1pt} {u_{\rm{z}}}}}{{\partial {\kern 1pt} t}} - \left( {\frac{{\partial {\kern 1pt} {\sigma _{{\rm{zz}}}}}}{{\partial {\kern 1pt} z}} + \frac{{\partial {\kern 1pt} {\sigma _{{\rm{xz}}}}}}{{\partial {\kern 1pt} x}}} \right) = 0}\\ {\frac{{\partial {\kern 1pt} {\sigma _{{\rm{xx}}}}}}{{\partial {\kern 1pt} t}} - \lambda \left( {\frac{{\partial {\kern 1pt} {u_{\rm{x}}}}}{{\partial {\kern 1pt} x}} + \frac{{\partial {\kern 1pt} {u_{\rm{z}}}}}{{\partial {\kern 1pt} z}}} \right) - 2\mu \frac{{\partial {u_{\rm{x}}}}}{{\partial x}} = {S_{{\rm{xx}}}}}\\ {\frac{{\partial {\kern 1pt} {\sigma _{{\rm{zz}}}}}}{{\partial {\kern 1pt} t}} - \lambda \left( {\frac{{\partial {\kern 1pt} {u_{\rm{x}}}}}{{\partial {\kern 1pt} x}} + \frac{{\partial {\kern 1pt} {u_{\rm{z}}}}}{{\partial {\kern 1pt} z}}} \right) - 2\mu \frac{{\partial {\kern 1pt} {u_{\rm{z}}}}}{{\partial {\kern 1pt} z}} = {S_{{\rm{zz}}}}}\\ {\frac{{\partial {\kern 1pt} {\sigma _{{\rm{xz}}}}}}{{\partial {\kern 1pt} t}} - \mu \left( {\frac{{\partial {\kern 1pt} {u_{\rm{x}}}}}{{\partial {\kern 1pt} z}} + \frac{{\partial {\kern 1pt} {u_{\rm{z}}}}}{{\partial {\kern 1pt} x}}} \right) = 0} \end{array}} \right. $

(1)

式中：S_xx和S_zz分别为水平方向和垂直方向上的震源；u_x和u_z分别为介质中在水平方向和垂直方向上质点振动速度分量；σ_xx，σ_xz，σ_zz分别代表应力分量；ρ为介质质量密度，kg/cm³，常数；λ和μ均为拉梅常数。拉梅常数λ和μ、密度ρ、P波速度v_P、S波速度v_S之间的关系为：λ = ρ(v_P² - 2 v_S²)，μ = ρv_S²。

与背景介质和扰动介质相对应的是背景波场P =[u_x，u_z，σ_xx，σ_zz，σ_xz]和扰动波场δu =[Δu_x，Δu_z，Δσ^xx，Δ σ^zz，Δσ_xz]。总波场同样满足波动方程，即

$ \left\{ \begin{array}{l} \begin{array}{*{20}{c}} {\rho \frac{{\partial ({u_{\rm{x}}} + \Delta {u_{\rm{x}}})}}{{\partial {\kern 1pt} t}} - \left[ {\frac{{\partial ({\sigma _{{\rm{xx}}}} + \Delta {\sigma _{{\rm{xx}}}})}}{{\partial {\kern 1pt} x}} + \frac{{\partial ({\sigma _{{\rm{xz}}}} + \Delta {\sigma _{{\rm{xx}}}})}}{{\partial {\kern 1pt} z}}} \right] = 0}\\ {\rho \frac{{\partial ({u_{\rm{z}}} + \Delta {u_{\rm{z}}})}}{{\partial {\kern 1pt} t}} - \left[ {\frac{{\partial ({\sigma _{{\rm{zz}}}} + \Delta {\sigma _{{\rm{zz}}}})}}{{\partial {\kern 1pt} z}} + \frac{{\partial ({\sigma _{{\rm{xz}}}} + \Delta {\sigma _{{\rm{xz}}}})}}{{\partial {\kern 1pt} x}}} \right] = 0} \end{array}\\ \begin{array}{*{20}{l}} {\frac{{\partial ({\sigma _{{\rm{xx}}}} + \Delta {\sigma _{{\rm{xx}}}})}}{{\partial {\kern 1pt} t}} - (\lambda + \delta \lambda )\left[ {\frac{{\partial ({u_{\rm{x}}} + \Delta {u_{\rm{x}}})}}{{\partial {\kern 1pt} x}} + \frac{{\partial ({u_{\rm{z}}} + \Delta {u_{\rm{z}}})}}{{\partial {\kern 1pt} z}}} \right] - 2(\mu + \delta \mu )\frac{{\partial ({u_{\rm{x}}} + \Delta {u_{\rm{x}}})}}{{\partial {\kern 1pt} x}} = {S_{{\rm{xx}}}}}\\ {\frac{{\partial ({\sigma _{{\rm{zz}}}} + \Delta {\sigma _{{\rm{zz}}}})}}{{\partial {\kern 1pt} t}} - (\lambda + \delta \lambda )\left[ {\frac{{\partial ({u_{\rm{x}}} + \Delta {u_{\rm{x}}})}}{{\partial {\kern 1pt} x}} + \frac{{\partial ({u_{\rm{z}}} + \Delta {u_{\rm{z}}})}}{{\partial {\kern 1pt} z}}} \right] - 2(\mu + \delta \mu )\frac{{\partial ({u_{\rm{z}}} + \Delta {u_{\rm{z}}})}}{{\partial {\kern 1pt} z}} = {S_{{\rm{zz}}}}} \end{array}\\ \frac{{\partial ({\sigma _{{\rm{xz}}}} + \delta {\sigma _{{\rm{xz}}}})}}{{\partial t}} - (\mu + \delta \mu )\left[ {\frac{{\partial ({u_{\rm{x}}} + \delta {\kern 1pt} {u_{\rm{x}}})}}{{\partial {\kern 1pt} z}} + \frac{{\partial ({u_{\rm{z}}} + \delta {\kern 1pt} {u_{\rm{z}}})}}{{\partial {\kern 1pt} x}}} \right] = 0 \end{array} \right. $

(2)

将式(2)与式(1)相减，并且在模型扰动足够小的假设下，用背景波场P代替总波场P_total，可得到基于Born近似的关于扰动波场δu的波动方程，可以表示为

$ \left\{ \begin{array}{l} \begin{array}{*{20}{l}} {\rho \frac{{\partial {\kern 1pt} \Delta {\kern 1pt} {u_{\rm{x}}}}}{{\partial {\kern 1pt} t}} - \left( {\frac{{\partial {\kern 1pt} \Delta {\kern 1pt} {\kern 1pt} {\sigma _{{\rm{xx}}}}}}{{\partial {\kern 1pt} x}} + \frac{{\partial {\kern 1pt} \Delta {\kern 1pt} {\kern 1pt} {\sigma _{{\rm{xz}}}}}}{{\partial {\kern 1pt} z}}} \right) = 0}\\ {\rho \frac{{\partial {\kern 1pt} \Delta {\kern 1pt} {u_{\rm{z}}}}}{{\partial {\kern 1pt} t}} - \left( {\frac{{\partial {\kern 1pt} \Delta {\kern 1pt} {\kern 1pt} {\sigma _{{\rm{zz}}}}}}{{\partial {\kern 1pt} z}} + \frac{{\partial {\kern 1pt} \Delta {\kern 1pt} {\kern 1pt} {\sigma _{{\rm{xz}}}}}}{{\partial {\kern 1pt} x}}} \right) = 0} \end{array}\\ \frac{{\partial {\kern 1pt} \Delta {\kern 1pt} {\kern 1pt} {\sigma _{{\rm{xx}}}}}}{{\partial {\kern 1pt} t}} - \lambda \left( {\frac{{\partial {\kern 1pt} \Delta {\kern 1pt} {u_{\rm{x}}}}}{{\partial {\kern 1pt} x}} + \frac{{\partial {\kern 1pt} \Delta {\kern 1pt} {u_{\rm{z}}}}}{{\partial {\kern 1pt} z}}} \right) - 2\mu \frac{{\partial {\kern 1pt} \Delta {\kern 1pt} {u_{\rm{x}}}}}{{\partial {\kern 1pt} x}} = (\lambda + 2\mu ){m_{\rm{P}}}\frac{{\partial {\kern 1pt} {u_{\rm{x}}}}}{{\partial {\kern 1pt} x}} + [(\lambda + 2\mu ){m_{\rm{P}}} - 2\mu {m_{\rm{S}}}]\frac{{\partial {\kern 1pt} {u_{\rm{z}}}}}{{\partial {\kern 1pt} z}}\\ \frac{{\partial {\kern 1pt} \Delta {\kern 1pt} {\kern 1pt} {\sigma _{{\rm{zz}}}}}}{{\partial {\kern 1pt} t}} - \lambda \left( {\frac{{\partial {\kern 1pt} \Delta {\kern 1pt} {u_{\rm{x}}}}}{{\partial {\kern 1pt} x}} + \frac{{\partial {\kern 1pt} \Delta {\kern 1pt} {u_{\rm{z}}}}}{{\partial {\kern 1pt} z}}} \right) - 2\mu \frac{{\partial {\kern 1pt} \Delta {\kern 1pt} {u_{\rm{z}}}}}{{\partial {\kern 1pt} z}} = (\lambda + 2\mu ){m_{\rm{P}}}\frac{{\partial {\kern 1pt} {u_{\rm{z}}}}}{{\partial {\kern 1pt} z}} + [(\lambda + 2\mu ){m_{\rm{P}}} - 2\mu {m_{\rm{S}}}]\frac{{\partial {\kern 1pt} {u_{\rm{x}}}}}{{\partial {\kern 1pt} x}}\\ \frac{{\partial {\kern 1pt} \Delta {\kern 1pt} {\kern 1pt} {\sigma _{{\rm{xz}}}}}}{{\partial {\kern 1pt} t}} - \mu \left( {\frac{{\partial {\kern 1pt} \Delta {\kern 1pt} {u_{\rm{x}}}}}{{\partial {\kern 1pt} z}} + \frac{{\partial {\kern 1pt} \Delta {\kern 1pt} {u_{\rm{z}}}}}{{\partial {\kern 1pt} x}}} \right) = \mu {m_{\rm{S}}}\frac{{\partial {\kern 1pt} {u_{\rm{x}}}}}{{\partial {\kern 1pt} z}} + \mu {m_{\rm{S}}}\frac{{\partial {\kern 1pt} {u_{\rm{z}}}}}{{\partial {\kern 1pt} x}} \end{array} \right. $

(3)

式中：Δu_x和Δu_z均为质点振动速度扰动项；Δσ_xx，Δσ_xz，Δσ_zz均为应力扰动项；P，S波反射系数模型可用速度扰动δv_P，δv_S和背景速度v_P，v_S表示为m_P (x) = 2 δv_P/v_P，m_S (x) = 2 δv_S/v_S，反映了介质速度扰动的相对大小。

从式(3)可以看出，要得到质点振动速度扰动项和应力扰动项需要经历2个过程：第1个过程是先通过式(1)得到背景波场P₀；第2个过程是利用式(3)进行计算，通过背景慢度模型和新的虚拟震源项再次激发，再传播到地下x位置处形成散射波场。另外，从式(3)的右端表达式可以看出，质点振动速度扰动项和应力扰动项的产生是由背景介质参数与介质速度扰动作为新的震源，而且两者与介质速度扰动项以及应力扰动项呈线性关系。

为了便于表示，用矩阵形式将式(3)简写成^[32-35]

$ L{\kern 1pt} m = d $

(4)

式中：m为反射系数模型的矩阵形式；d表示矩阵形式的观测数据，d =(d_x，d_z，0，0，0)T；L为基于Born近似下的线性反偏移算子，它建立了介质参数与正演数据之间的联系。

采用时间2阶、空间10阶精度的交错网格有限差分法对式(1)进行差分离散。在进行波场模拟时，由于人工截断边界会导致边界反射，采用完全匹配层(PML)吸收边界，来压制边界反射提高模拟区域的精度与信噪比。

本文借助伴随状态法得到伴随算子及相应的伴随波场控制方程，由伴随状态法可得

$ \langle {L^*}w,\delta u\rangle = \langle w,L{\kern 1pt} \delta {\kern 1pt} u\rangle $

(5)

式中：w为伴随波场向量，$w = {\left[{{w_{\rm{x}}}, {w_{\rm{z}}}, {{\hat \sigma }_{{\rm{xx}}}}, {{\hat \sigma }_{{\rm{zz}}}}, {{\hat \sigma }_{{\rm{xz}}}}} \right]^T}$；L^*为L的伴随算子。

则伴随方程可表示为

$ \left\{ \begin{array}{l} \begin{array}{*{20}{c}} {\rho \frac{{\partial {\kern 1pt} {w_{\rm{x}}}}}{{\partial {\kern 1pt} t}} - \frac{{\partial [(\lambda + 2\mu ){{\hat \sigma }_{{\rm{xx}}}}]}}{{\partial {\kern 1pt} x}} - \frac{{\partial (\lambda {\kern 1pt} {{\hat \sigma }_{{\rm{zz}}}})}}{{\partial {\kern 1pt} x}} - \frac{{\partial (\mu {\kern 1pt} {{\hat \sigma }_{{\rm{xz}}}})}}{{\partial {\kern 1pt} z}} = - \delta {d_{\rm{x}}}({x_{\rm{g}}},t;{x_{\rm{s}}})}\\ {\rho \frac{{\partial {\kern 1pt} {w_{\rm{z}}}}}{{\partial {\kern 1pt} t}} - \frac{{\partial [(\lambda + 2\mu ){{\hat \sigma }_{{\rm{zz}}}}]}}{{\partial {\kern 1pt} z}} - \frac{{\partial (\lambda {\kern 1pt} {{\hat \sigma }_{{\rm{xx}}}})}}{{\partial {\kern 1pt} z}} - \frac{{\partial (\mu {\kern 1pt} {{\hat \sigma }_{{\rm{xz}}}})}}{{\partial {\kern 1pt} x}} = - \delta {d_{\rm{z}}}({x_{\rm{g}}},t;{x_{\rm{s}}})} \end{array}\\ \begin{array}{*{20}{l}} {\frac{{\partial {\kern 1pt} {\kern 1pt} {{\hat \sigma }_{{\rm{xx}}}}}}{{\partial {\kern 1pt} t}} - \frac{{\partial {\kern 1pt} {w_{\rm{x}}}}}{{\partial {\kern 1pt} x}} = 0}\\ {\frac{{\partial {\kern 1pt} {\kern 1pt} {{\hat \sigma }_{{\rm{zz}}}}}}{{\partial {\kern 1pt} t}} - \frac{{\partial {\kern 1pt} {w_{\rm{z}}}}}{{\partial {\kern 1pt} z}} = 0}\\ {\frac{{\partial {\kern 1pt} {\kern 1pt} {{\hat \sigma }_{{\rm{xz}}}}}}{{\partial {\kern 1pt} t}} - \frac{{\partial {\kern 1pt} {w_{\rm{x}}}}}{{\partial {\kern 1pt} z}} - \frac{{\partial {\kern 1pt} {w_{\rm{z}}}}}{{\partial {\kern 1pt} x}} = 0} \end{array} \end{array} \right. $

(6)

式中：δd为残差向量，$\delta d = {\left[{\begin{array}{*{20}{l}} {\delta {d_x}}&{\delta {d_z}}&0&0&0 \end{array}} \right]^T}$。

为了得到与记录数据最佳匹配的偏移结果，把最小二乘逆时偏移引入反演中，其误差函数定义为

$ \chi (m) = \frac{1}{2}{\left\| {Lm - {d_{{\rm{obs}}}}} \right\|^2} = \frac{1}{2}{\left\| {\delta d} \right\|^2} $

(7)

根据伴随状态方法，求解式(7)可得λ，μ的梯度表达式为

$ \left\{ {\begin{array}{*{20}{l}} {\delta {\kern 1pt} \lambda = \frac{{\partial {\kern 1pt} \chi }}{{\partial {\kern 1pt} \lambda }} = - \int {\left[ {\left( {\frac{{\partial {u_{\rm{x}}}}}{{\partial x}} + \frac{{\partial {u_{\rm{z}}}}}{{\partial z}}} \right)({{\hat \sigma }_{{\rm{xx}}}} + {{\hat \sigma }_{{\rm{zz}}}})} \right]{\rm{d}}{\kern 1pt} t} }\\ {\delta {\kern 1pt} \mu = \frac{{\partial {\kern 1pt} \chi }}{{\partial {\kern 1pt} \mu }} = - \int {\left[ {2\frac{{\partial {u_{\rm{x}}}}}{{\partial x}}{{\hat \sigma }_{{\rm{xx}}}} + 2\frac{{\partial {u_{\rm{z}}}}}{{\partial z}}{{\hat \sigma }_{{\rm{zz}}}} + \left( {\frac{{\partial {u_{\rm{x}}}}}{{\partial z}} + \frac{{\partial {u_{\rm{z}}}}}{{\partial x}}} \right){{\hat \sigma }_{{\rm{xz}}}}} \right]{\rm{d}}{\kern 1pt} t} } \end{array}} \right. $

(8)

v_P和v_S的梯度可由链式公式求得，即

$ \left\{ {\begin{array}{*{20}{l}} {\delta {v_{\rm{P}}} = 2{v_{\rm{P}}}\rho \frac{{\partial {\kern 1pt} \chi }}{{\partial {\kern 1pt} \lambda }}}\\ {\delta {v_{\rm{S}}} = - 4{v_{\rm{S}}}\rho \frac{{\partial {\kern 1pt} \chi }}{{\partial {\kern 1pt} \lambda }} + 2{v_{\rm{S}}}\rho \frac{{\partial {\kern 1pt} \chi }}{{\partial {\kern 1pt} \mu }}} \end{array}} \right. $

(9)

梯度计算是最小二乘成像的核心问题之一，对梯度做预条件，加速问题的收敛。本文中用于后续迭代更新的梯度为g = (δv_P，δv_S)。将式(9)代入式(7)可得δv_P和δv_S的梯度表达式，即

$ \left\{ {\begin{array}{*{20}{l}} {\delta {v_{\rm{P}}} = - 2{v_{\rm{P}}}\rho \int {\left[ {\left( {\frac{{\partial {\kern 1pt} {u_{\rm{x}}}}}{{\partial {\kern 1pt} x}} + \frac{{\partial {\kern 1pt} {u_{\rm{z}}}}}{{\partial {\kern 1pt} z}}} \right)({{\hat \sigma }_{{\rm{xx}}}} + {{\hat \sigma }_{{\rm{zz}}}})} \right]{\rm{d}}{\kern 1pt} t} }\\ {\delta {v_{\rm{S}}} = - 2{v_{\rm{S}}}\rho \int {\left[ { - 2\frac{{\partial {\kern 1pt} {u_{\rm{x}}}}}{{\partial {\kern 1pt} x}}{{\hat \sigma }_{{\rm{xx}}}} - 2\frac{{\partial {\kern 1pt} {u_{\rm{z}}}}}{{\partial {\kern 1pt} z}}{{\hat \sigma }_{{\rm{zz}}}} + \left( {\frac{{\partial {\kern 1pt} {u_{\rm{x}}}}}{{\partial {\kern 1pt} z}} + \frac{{\partial {\kern 1pt} {u_{\rm{z}}}}}{{\partial {\kern 1pt} x}}} \right){{\hat \sigma }_{{\rm{xz}}}}} \right]{\rm{d}}{\kern 1pt} t} } \end{array}} \right. $

(10)

1.2 梯度预条件算子

式(7)所对应的解的估计称最小二乘偏移成像的结果，真实反射模型的最小二乘解：

$ m = {({L^*}L)^{ - 1}}{L^*}{d_{{\rm{obs}}}} $

(11)

传统的偏移方法是通过反偏移算子的共轭转置实现的，会造成地震信息缺失。从式(9)可看出，最小二乘偏移反演通过Hessian算子L^* L对偏移成像进行修正，难点就在于求解一个精确的Hessian并求逆。Hessian算子H是误差函数在m处对模型的二阶导数，即

$ \begin{array}{*{20}{l}} {H = \frac{{{\partial ^2}\chi }}{{\partial {\kern 1pt} {m^2}}} = \frac{{{\partial ^2}\chi }}{{\partial {\kern 1pt} {m^2}}}(d - {d_{{\rm{obs}}}}) + \frac{{\partial {\kern 1pt} \chi }}{{\partial {\kern 1pt} m}}\frac{{\partial {\kern 1pt} \chi }}{{\partial {\kern 1pt} m}}}\\ {{\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} = \frac{{\partial {\kern 1pt} {L^*}}}{{\partial {\kern 1pt} m}}(d - {d_{{\rm{obs}}}}) + {L^*}L} \end{array} $

(12)

从式(12)可以看出，H由2部分构成，其中，$\frac{\partial L^{*}}{\partial m}\left(d-d_{\mathrm{obs}}\right)$为其非线性项且依赖于模型，由于非线性项的计算过于复杂，计算成本过高，通常使用Hessian算子的线性项L^* L来近似代替全汉森矩阵。通过Hessian算子的逆矩阵作为梯度预处理算子，可以得到分辨率更高的成像结果，并且加快收敛。

式(1)用矩阵形式可表示为

$ S{\kern 1pt} w = f $

(13)

其中：

$ \begin{array}{l} S = \left[ {\begin{array}{*{20}{c}} {\rho \frac{\partial }{{\partial {\kern 1pt} t}}}&0&{ - \frac{\partial }{{\partial {\kern 1pt} x}}}&0&{ - \frac{\partial }{{\partial {\kern 1pt} z}}}\\ 0&{\rho \frac{\partial }{{\partial {\kern 1pt} t}}}&0&{ - \frac{\partial }{{\partial {\kern 1pt} z}}}&{ - \frac{\partial }{{\partial {\kern 1pt} x}}}\\ { - (\lambda + 2\mu )\frac{\partial }{{\partial {\kern 1pt} x}}}&{ - \lambda \frac{\partial }{{\partial {\kern 1pt} z}}}&{\frac{\partial }{{\partial {\kern 1pt} t}}}&0&0\\ { - \lambda \frac{\partial }{{\partial {\kern 1pt} x}}}&{ - (\lambda + 2\mu )\frac{\partial }{{\partial {\kern 1pt} z}}}&0&{\frac{\partial }{{\partial t}}}&0\\ { - \mu \frac{\partial }{{\partial {\kern 1pt} z}}}&{ - \mu \frac{\partial }{{\partial {\kern 1pt} x}}}&0&0&{\frac{\partial }{{\partial {\kern 1pt} t}}} \end{array}} \right],\\ w = \left[ {\begin{array}{*{20}{l}} {{u_{\rm{x}}}}\\ {{u_{\rm{z}}}}\\ {{\sigma _{{\rm{xx}}}}}\\ {{\sigma _{{\rm{zz}}}}}\\ {{\sigma _{{\rm{xz}}}}} \end{array}} \right],f = \left[ {\begin{array}{*{20}{l}} 0\\ 0\\ {{S_{{\rm{xx}}}}}\\ {{S_{{\rm{zz}}}}}\\ 0 \end{array}} \right]。\end{array} $

对式(13)求导可得

$ {\frac{{\partial {\kern 1pt} S}}{{\partial {\kern 1pt} m}}w + S\frac{{\partial {\kern 1pt} w}}{{\partial {\kern 1pt} m}} = 0} $

(14)

移项：

$ {\frac{{\partial {\kern 1pt} w}}{{\partial {\kern 1pt} m}} = - {S^{ - 1}}\frac{{\partial {\kern 1pt} S}}{{\partial {\kern 1pt} m}}w} $

(15)

通过汉森矩阵的线性项来近似代替全汉森矩阵，即

$ H = {L^*}L = {\left( {\frac{{\partial {\kern 1pt} w}}{{\partial {\kern 1pt} m}}} \right)^*}\left( {\frac{{\partial {\kern 1pt} w}}{{\partial {\kern 1pt} m}}} \right) = {\left( {{S^{ - 1}}\frac{{\partial {\kern 1pt} S}}{{\partial {\kern 1pt} m}}{\kern 1pt} {\kern 1pt} w} \right)^*}\left( {{S^{ - 1}}\frac{{\partial {\kern 1pt} S}}{{\partial {\kern 1pt} m}}w} \right) $

(16)

用拉梅常数表示的预条件算子表达式为

$ \left\{ {\begin{array}{*{20}{l}} {{H^\lambda } = \int {\frac{{{{({{\dot \sigma }_{{\rm{xx}}}} + {{\dot \sigma }_{{\rm{zz}}}})}^2}}}{{2{{(\lambda + \mu )}^2}}}{\rm{d}}{\kern 1pt} t} }\\ {{H^\mu } = \int {\left[ {\frac{{\dot \sigma _{{\rm{xz}}}^2}}{{{\mu ^2}}} + \frac{{{{({{\dot \sigma }_{{\rm{xx}}}} + {{\dot \sigma }_{{\rm{zz}}}})}^2}}}{{2{{(\lambda + \mu )}^2}}} + \frac{{{{({{\dot \sigma }_{{\rm{xx}}}} - {{\dot \sigma }_{{\rm{zz}}}})}^2}}}{{2{\mu ^2}}}} \right]{\rm{d}}{\kern 1pt} t} } \end{array}} \right. $

(17)

1.3 预条件弹性波最小二乘逆时偏移

利用共轭梯度法迭代求解式(7)得到模拟数据与观测数据的差别最小的最佳反射系数模型，以P表示预条件算子，梯度g = (δv_P, δv_S)，具体实现流程如下：

(1) 输入初始模型(初始模型赋0值，即第1次迭代的结果为逆时偏移的偏移剖面)，观测数据d_obs，梯度误差阈值δ。

(2) 计算第i次迭代共轭梯度dkⁱ和共轭方向修正因子β

$ \left\{ {\begin{array}{*{20}{l}} {d{\kern 1pt} {k^i} = P{\kern 1pt} {g^i} + \beta {\kern 1pt} d{\kern 1pt} {k^{i - 1}}}\\ {\beta = \frac{{{g^i}P{\kern 1pt} {g^i}}}{{{g^{i - 1}}P{\kern 1pt} {g^{i - 1}}}}} \end{array}} \right. $

(18)

(3) 计算步长α并对成像结果进行更新，得到m^{i + 1}

$ \left\{ {\begin{array}{*{20}{l}} {\alpha = \frac{{{{(d{\kern 1pt} {k^i})}^T}{g^i}}}{{{{(L{\kern 1pt} d{\kern 1pt} {k^i})}^T}L{\kern 1pt} d{\kern 1pt} {k^i}}}}\\ {{m^{i + 1}} = {m^i} - \alpha {\kern 1pt} d{\kern 1pt} {k^i}} \end{array}} \right. $

(19)

(4) 判断是否满足终止迭代条件，若满足，则退出迭代并保存成像结果；不满足，则回到第(2)，(3)步继续迭代直到满足迭代终止条件，最后一次迭代即为最终的成像结果。

2 模型试算 2.1 盐丘模型

采用SEG/EAGE提供的盐丘模型来测试常规ELSRTM和P-ELSRTM对模型的成像效果，纵横波真实速度场和反射系数模型如图 1所示。设定纵横波速度比为1.73，密度为常数。模型横向网格点数为645，纵向网格点数为150，采样间隔为10 m。利用该模型主要考察P-ELSRTM成像算法是否能够使盐下构造准确成像以及对振幅的补偿效果：由于盐丘的屏蔽作用，导致盐下构造的振幅较小，另外，盐丘侧翼发育的高陡小构造会引起严重的构造假象。观测系统采用全接收方式，每炮的网格间隔数为10 m，共60炮。震源类型为爆炸震源，震源函数采用主频为25 Hz的雷克子波。

ERTM成像结果如图 2所示，逆时偏移成像中使用全波动方程模拟波场，可以实现对陡倾角构造的成像，但是在ERTM剖面中含有较强的低波数噪声和震源效应[图 2(a)和图 2(b)]。应用拉普拉斯算子去噪之后，相比未作滤波的ERTM结果得到了改善[图 2(c)和图 2(d)]，模型浅部的低频噪音得到了压制，但成像能量不均衡。使用预条件算子和滤波之后，相比滤波后的ERTM成像剖面质量得到进一步提高[图 2(e)和图 2(f)]，盐下构造清晰，且模型浅部的低频噪音和震源效应也得到了更好地压制。

为了讨论预条件ELSRTM的有效性及优点，分别对常规ELSRTM和预条件P-ELSRTM成像结果进行对比(图 3)：①未作滤波的常规ELSRTM在迭代次数为5和30次的成像结果中都残留较严重的低波数噪音，影响对盐丘之上沉积层的识别[图 3 (a)-(d)]，而应用高通滤波预条件算子后，成像剖面中的低频噪音得到了压制，剖面振幅更加均衡[图 3(e)-(h)]，证明了滤波算子的有效性；②滤波后的常规ELSRTM成像结果中盐体边界清晰，盐下成像相对于ERTM结果更为准确，但仍存在深部成像能量不足的现象，并且对顶部低频噪声的压制不彻底[图 3(e)-(h)]；③应用照明补偿预条件算子和高通滤波后的结果[图 3(i)-(l)]表明，采用预条件成像方法实现了对复杂模型深部小断层较为清晰的成像，盐下断层及基底均能正确成像，并且成像剖面信噪比更高、能量更加均衡。

分别从常规ELSRTM和P-ELSRTM的成像结果中抽取偏移距为4 900 m的地震单道以对比各方法的保幅性(迭代30次)，结果如图 4所示。从图 4可看出：①在模型的浅部区域由于低频噪音和震源效应的存在，常规ELSRTM的振幅曲线会有明显跳动，而P-ELSRTM的振幅曲线抖动幅度较小；②在模型的500 m，830 m，1 430 m，4 500 m(反射界面)处，常规ELSRTM的振幅与真实反射系数之间有较大的差异，而P-ELSRTM的振幅与真实反射系数更加接近。因此，在浅-中深层P-ELSRTM都比常规ELSRTM具有更好的保幅性。

常规ELSRTM和P-ELSRTM的成像结果与真实结果的归一化残差随迭代次数变化的结果如图 5所示。从图 5可以看出：P-ELSRTM的成像反演结果在前60次迭代中都能稳定收敛，最终收敛到22.6%；常规ELSRTM随着迭代次数的增加只能收敛到52%便停止收敛。因此，使用梯度预条件算子能够使复杂模型成像结果的残差收敛更快且达到更低的值。

在这次测试中，CPU型号为Intel(R)Xeon(R) CPU E5-2650 v2@2.60 GHz，P-ELSRTM和常规EL‐ SRTM均串行运行，迭代一次所花费的时间分别为4 667 s和4 642 s。基于不同最优化算法的LSRTM结果的对比如表 1所列。由于计算汉森算子所花费的计算量相对于偏移运算较小，因此表 1中忽略了汉森算子的计算量，仅对比了偏移运算的计算量。P-ELSRTM算法与ELSRTM算法的计算量相近，但前者具有更高的成像质量、更快的收敛速度及更好的稳定性。

野外地震资料采集的数据中不可避免地含有环境以及人为因素造成的随机噪声。为了测试PELSRTM对低信噪比数据的适应性，对炮记录添加低信噪比的随机噪声，使其X分量和Z分量信噪比为0.5 dB(图 6)。从ELSRTM和P-ELSRTM的成像结果(图 7)可看出，对于低信噪比数据，ELSRTM和P-ELSRTM结果中都含有类似于随机噪声的偏移假象，但相比常规ELSRTM，P-ELSRTM能够更好地压制低频噪音。从残差收敛曲线可以看出，P-ELSRTM的收敛速度较常规ELSRTM快，并且迭代60次后其残差较小(图 8)。因此，P-ELSRTM在处理含噪数据方面有较大优势，抗噪性更好。

3 结论与讨论

(1) P-ELSRTM和ELSRTM能够有效压制浅层低频噪音，消除震源效应，有效改善常规ERTM的成像质量，具有更高的成像分辨率。

(2) P-ELSRTM在中深部具有更好的保幅性和分辨率，有利于深部及盐下构造的识别；P-ELSRTM比常规ELSRTM的残差曲线更快地收敛。在对含噪数据的成像测试中，P-ELSRTM比常规ELSRTM可以收敛到一个更低的残差水平。

(3) 本文使用一阶位移波动方程得到的纵横波是耦合在一起的，因此通过迭代无法完全消除纵横波串扰噪音。针对该问题，进行纵波和横波波场矢量波场分离可进一步提高成像精度。