0 引言

储层孔隙中含有油气水等流体，会引起地震波传播过程中的强振幅衰减和相位畸变^[1]，从而导致传统的声波或者弹性波偏移方法分辨率降低或者失效。关于地震波在黏弹性介质中的传播机理，主要有Maxwell模型、Kelvin-Vogit模型、标准线性固体、广义标准线性体、分数阶常Q模型等^[2-6]。早期的地震资料衰减补偿主要是使用反Q滤波在数据域实现^[7-10]，如Blanch等^[11-12]提出用τ方法模拟常Q，并基于伴随状态法在τ-p域中用线性反演方法恢复黏声介质的短波长分量。考虑到地震波衰减是在传播过程中发生的，融合衰减补偿的偏移方法越来越受到重视。Xin等^[13]和Xie等^[14]利用射线追踪的方法在叠前深度偏移过程中进行衰减补偿；白敏等^[15]和代福材等^[16]在高斯束叠前深度偏移过程中补偿振幅衰减和校正相位畸变；Dai等^[17]、Yu等^[18]、Wang^[19]和Valenciano等^[20]在频率域用单程波动方程偏移进行衰减补偿；陈志德等^[21]引入等效Q值实现频率域波场延拓，基于单程波方程和稳定相点原理实现黏声叠前时间偏移，并提高了补偿算法的稳定性。朱峰等^[22]推导了黏声VTI介质高阶傅里叶有限差分波场延拓算子，实现了适应高角度入射波场的黏声叠前深度偏移方法。

不同学者提出了多种黏声波动方程，在逆时偏移过程中分别对振幅损失和相位畸变进行补偿和校正^[23-27]。李振春等^[28]基于广义标准线性体模型对黏声介质波动方程进行线性化，并借助伴随状态法实现了黏声介质最小二乘逆时偏移(Q-LSRTM)的迭代算法，获得了较高的成像分辨率。Dutta等^[29]、Dai等^[30]基于一阶松弛机制时间域黏声波动方程，用迭代的Q-LSRTM有效补偿和校正了强衰减层导致的振幅损失和相位畸变。徐凯等^[31]在反演框架下，构建了Q-LSRTM多约束正则化目标函数，并采用变分方法求解，消除了震源波场与接收点波场互相关引起的噪声，提高了复杂地质体成像的精度、增强了保幅性。相比于传统声波LSRTM，Q-LSRTM能实现反射体的精确归位且振幅更均衡，收敛速率更快^[32]。但Q-LSRTM也存在一定问题：其一，由于在传统Q-LSRTM中用于反传数据残差的伴随Q传播算子也是衰减的，导致成像结果分辨率降低；其二，基于时间域黏声波动的LSRTM计算成本较高，尤其是需要很多次迭代才能达到预期成像效果。在完全弹性情况下，为了提高传统LSRTM的收敛速度，李闯等^[33]发展了一种预条件偏移算法，提高了盐下构造成像的分辨率和保幅性。

在实际野外地震勘探中，由于稀疏的震源和检波器，记录偏移孔径有限、速度变化快及构造复杂，以及不均匀的地震照明等原因，偏移剖面当中往往存在多种假象^[34]。为了消除偏移假象、提升成像分辨率，偏移反褶积被广泛采用^[35]。偏移成像本质上就是真实反射率的一种模糊，模糊算子即为Hessian算子，而偏移反褶积处理就等效于一种近似的Hessian逆处理。Hu等^[34-35]在空间—波数域设计了一个偏移反褶积算子压制叠后偏移假象；Yu等^[36]将偏移反褶积从叠后偏移拓展到叠前深度偏移。Guitton^[37]直接在空间域构建了一系列匹配滤波器近似Hessian矩阵的逆实现偏移反褶积，后来这些匹配滤波器被Aoki等^[38]作为预条件化因子用于传统的LSRTM，获得了良好的成像效果，这些匹配滤波器也被称为去模糊滤波器。针对多震源LSRTM，Dai等^[39-40]用去模糊滤波器减少串扰噪声并且加速多震源LSRTM的收敛。但目前偏移反褶积的原理和实施都是基于无耗散介质假设，而地下储层往往具有衰减性质。

为了克服Q-LSRTM的缺陷，本文将匹配滤波偏移反褶积与Q-LSRTM相结合。首先基于点扩散函数(PSF)构建黏声去模糊滤波器，然后将其作为Q-LSRTM迭代过程中的预条件化因子，最终实现高分辨的衰减补偿偏移。相比于Q-LSRTM，预条件Q-LSRTM成像结果分辨率更高，反射体振幅更均衡，收敛速度更快。就计算成本来说，尽管预条件Q-LSRTM需要预先计算黏声去模糊滤波器，但其达到预期成像效果需要的迭代次数远比Q-LSRTM少。

1 理论方法 1.1 点扩散函数与去模糊滤波器

在声学介质中，基于Born近似，观测地震记录d的形成过程可以表示为

$ \mathit{\boldsymbol{d}} = \mathit{\boldsymbol{L}}{\mathit{\boldsymbol{m}}_0} $

(1)

式中：L表示线性模拟算子，与观测系统、震源子波和地下介质模型参数有关；m₀是地下反射率模型。对应L，可以获得其伴随偏移算子L^T，则偏移结果m_mig可表示为

$ {\mathit{\boldsymbol{m}}_{{\rm{mig}}}} = {\mathit{\boldsymbol{L}}^{\rm{T}}}\mathit{\boldsymbol{d}} = \overbrace {{\mathit{\boldsymbol{L}}^{\rm{T}}}\mathit{\boldsymbol{L}}}^{模糊算子}{\mathit{\boldsymbol{m}}_0} $

(2)

m_mig其实是真实反射率模型m₀的模糊版本，L^TL就是一个模糊算子，即偏移格林函数^[41]，也被称为点扩散函数(PSF)。由于震源子波、观测系统等多种因素的影响，偏移成像中往往存在较多的假象。为了获得更清晰的成像结果，可以取L^TL的逆并作用于成像结果

$ {\mathit{\boldsymbol{m}}_0} = \overbrace {{{\left( {{\mathit{\boldsymbol{L}}^{\rm{T}}}\mathit{\boldsymbol{L}}} \right)}^{ - 1}}}^{去模糊算子}{\mathit{\boldsymbol{m}}_{{\rm{mig}}}} $

(3)

但在地震成像时，直接计算(L^TL)^-1是不现实的，为此采用迭代LSRTM方法予以求解，但其计算量往往比标准偏移成像大一个数量级。为了加速收敛，设计去模糊算子的近似逆算子(L^TL)^-1，并在LSRTM中作为预条件化算子^[35-39]。

实际应用中，去模糊滤波器基于背景模型即偏移速度场和参考反射率模型获得，以均匀介质速度背景场和均匀分布点散射体反射率模型为例阐述去模糊滤波器的获取过程。图 1a为均匀分布点散射体模型，图 1b为对应的叠前深度偏移成像结果，其中震源和检波器均匀分布于模型上方。单个散射体的PSF响应即为其偏移结果。记参考反射率模型为m_ref，则在Bron近似下利用式(1)可获得正演模拟数据d_ref，继而可以获得标准的参考偏移结果

$ {\mathit{\boldsymbol{m}}_{{\rm{mig - ref}}}} = {\mathit{\boldsymbol{L}}^{\rm{T}}}\mathit{\boldsymbol{L}}{\mathit{\boldsymbol{m}}_{{\rm{ref}}}} = {\mathit{\boldsymbol{L}}^{\rm{T}}}{\mathit{\boldsymbol{d}}_{{\rm{ref}}}} $

(4)

将点散射体模型分解成多个子区域(如图 1a中黑色矩形框)，使点散射体位于矩形窗口的中央，窗口大小为w_x×w_z，足以覆盖其PSF响应的有效能量(图 1b中矩形框)。在这个局部窗口内，假设PSF保持不变。在w_x×w_z局部窗口区域内，去模糊滤波器将参考偏移成像结果和参考反射率模型联系起来^{[38, 40]}，即

$ \begin{array}{l} {\mathit{\boldsymbol{m}}_{{\rm{ref}}}}{\left( {x,z} \right)_i} = \int_{{V_i}} {{F_i}\left( {x - {x_0},z - {z_0}} \right) \times } \\ \;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;{\mathit{\boldsymbol{m}}_{{\rm{mig - ref}}}}{\left( {{x_0},{z_0}} \right)_i}{\rm{d}}V \end{array} $

(5)

式中：i为局部窗口序号；积分限V_i为整个局部窗口区域。式(5)写成矩阵形式为

$ {\left[ {{\mathit{\boldsymbol{m}}_{{\rm{ref}}}}} \right]_i} = {\mathit{\boldsymbol{F}}_i} \otimes {\left[ {{\mathit{\boldsymbol{m}}_{{\rm{mig - ref}}}}} \right]_i} $

(6)

式中：“⊗”表示空间卷积；[m_ref]_i和[m_mig-ref]_i分别为第i个局部窗口内参考反射率模型和参考偏移成像结果。与PSF相同，在每个局部窗口内去模糊滤波器保持不变，定义滤波器的空间尺寸为f_x×f_z。

由卷积的交换律可得

$ {\mathit{\boldsymbol{F}}_i} \otimes {\left[ {{\mathit{\boldsymbol{m}}_{{\rm{mig - ref}}}}} \right]_i} = {\left[ {{\mathit{\boldsymbol{m}}_{{\rm{mig - ref}}}}} \right]_i} \otimes {\mathit{\boldsymbol{F}}_i} $

(7)

将上述卷积运算写成矩阵运算形式，参考偏移结果[m_min-ref]_i变成一个(N+M-1)×N阶卷积矩阵，去模糊滤波器F_i变为N×1阶向量f，其中N=f_x×f_z且M=w_x×w_z，即为

$ {\left[ {{\mathit{\boldsymbol{m}}_{{\rm{mig - ref}}}}} \right]_i}{\mathit{\boldsymbol{f}}_i} = {\left[ {{\mathit{\boldsymbol{m}}_{{\rm{ref}}}}} \right]_i} $

(8)

上式为超定方程，可以采用最小二乘方法求解，在其两端乘以[m_mig-ref]_i^T变成标准方程形式

$ \left[ {{\mathit{\boldsymbol{m}}_{{\rm{mig - ref}}}}} \right]_i^{\rm{T}}{\left[ {{\mathit{\boldsymbol{m}}_{{\rm{mig - ref}}}}} \right]_i}{\mathit{\boldsymbol{f}}_i} = \left[ {{\mathit{\boldsymbol{m}}_{{\rm{mig - ref}}}}} \right]_i^{\rm{T}}{\left[ {{\mathit{\boldsymbol{m}}_{{\rm{ref}}}}} \right]_i} $

(9)

继而采用三角(LU)分解方法求解式(9)即可获得去模糊滤波器，在第i个局部窗口内近似等于(L^TL)^-1。在实际应用中，为了保证计算效果，建议局部窗口长度不超过1.5倍地震波长。

1.2 黏声介质中的去模糊滤波器

将去模糊滤波器应用到偏移成像的过程就是一种偏移反褶积处理。前人对偏移反褶积的研究都集中在无耗散介质当中，但强吸收会对地震波振幅和相位产生明显的衰减和畸变^[1]。为了补偿衰减效应，Dutta等^[29]发展了时间域Q-LSRTM技术，成像结果归位更精确，振幅更均衡。在传统Q-LSRTM中，用于反传数据残差的伴随Q传播算子也是衰减的，因此成像分辨率有一定的损失。

这种分辨率的损失能通过黏声介质中的PSF即偏移格林函数予以解释。给定速度为v₀的均匀介质和位于x_s处角频率为ω的简谐点源，声学格林函数可以表示为

$ G\left( {\mathit{\boldsymbol{x}},{\mathit{\boldsymbol{x}}_s}} \right) = \frac{{\exp \left( {\frac{{{\rm{i}}\omega \left| {{\mathit{\boldsymbol{x}}_s} - \mathit{\boldsymbol{x}}} \right|}}{{{v_0}}}} \right)}}{{\left| {{\mathit{\boldsymbol{x}}_s} - \mathit{\boldsymbol{x}}} \right|}} $

(10)

如果介质是耗散的，则将完全弹性介质中的实速度v₀替换成复速度^[1]

$ v\left( \omega \right) = {v_0}\left( {1 + \frac{1}{{\mathit{\pi }Q}}\ln \frac{\omega }{{{\omega _0}}}} \right)\left( {1 - \frac{{\rm{i}}}{{2Q}}} \right) $

(11)

从而得到黏声介质格林函数

$ \begin{array}{l} G\left( {\mathit{\boldsymbol{x}},{\mathit{\boldsymbol{x}}_s}} \right) = \frac{1}{{\left| {{\mathit{\boldsymbol{x}}_s} - \mathit{\boldsymbol{x}}} \right|}} \times \\ \;\;\;\;\;\;\;\overbrace {\exp \left( {{\rm{i}}\omega \frac{{\left| {{\mathit{\boldsymbol{x}}_s} - \mathit{\boldsymbol{x}}} \right|}}{{{v_0}\xi }}} \right)}^{相位畸变}\overbrace {\exp \left( { - \frac{{\omega \left| {{\mathit{\boldsymbol{x}}_s} - \mathit{\boldsymbol{x}}} \right|}}{{2Q{v_0}\xi }}} \right)}^{振幅/频率衰减} \end{array} $

(12)

式中：Q是品质因子；ω₀是参考角频率；$ \xi = \left[ {1 + \frac{1}{{{\rm{ \mathsf{ π} }}Q}}\ln \left( {\frac{\omega }{{{\omega _0}}}} \right)} \right]\left( {1 + \frac{1}{{4{Q^2}}}} \right)$。式(12)中第一个指数项表示相位畸变，第二个指数项表示振幅/高频损失。据此可得声学和黏声情况下的PSF为

$ \begin{array}{l} {\mathit{\Gamma }_{\rm{A}}} = \int_\omega {{\rm{d}}\omega \sum\limits_s {\sum\limits_r \times } } \\ \;\;\;\;\frac{{\exp \left( {\frac{{{\rm{i}}\omega R}}{{{v_0}\xi }}} \right)\exp \left( { - \frac{{\omega R}}{{2Q{v_0}\xi }}} \right)\exp \left( { - \frac{{{\rm{i}}\omega {R_0}}}{{{v_0}}}} \right)}}{{\left| {{\mathit{\boldsymbol{x}}_s} - \mathit{\boldsymbol{x}}} \right| \times \left| {{\mathit{\boldsymbol{x}}_r} - \mathit{\boldsymbol{x}}} \right| \times \left| {{\mathit{\boldsymbol{x}}_s} - {\mathit{\boldsymbol{x}}_0}} \right| \times \left| {{\mathit{\boldsymbol{x}}_r} - {\mathit{\boldsymbol{x}}_0}} \right|}} \end{array} $

(13)

$ \begin{array}{l} {\mathit{\Gamma }_{\rm{Q}}} = \int_\omega {{\rm{d}}\omega \sum\limits_s {\sum\limits_r \times } } \\ \;\;\;\;\frac{{\exp \left( {\frac{{{\rm{i}}\omega \left( {R - {R_0}} \right)}}{{{v_0}\xi }}} \right)\exp \left( { - \frac{{\omega \left( {R + {R_0}} \right)}}{{2Q{v_0}\xi }}} \right)}}{{\left| {{\mathit{\boldsymbol{x}}_s} - \mathit{\boldsymbol{x}}} \right| \times \left| {{\mathit{\boldsymbol{x}}_r} - \mathit{\boldsymbol{x}}} \right| \times \left| {{\mathit{\boldsymbol{x}}_s} - {\mathit{\boldsymbol{x}}_0}} \right| \times \left| {{\mathit{\boldsymbol{x}}_r} - {\mathit{\boldsymbol{x}}_0}} \right|}} \end{array} $

(14)

式中：R=|x_s-x|+|x_r-x|；R₀=|x_s-x₀|+|x_r-x₀|，x_r为检波点坐标。从式(13)可以看出，当采用声波偏移(声学伴随算子)对黏声数据进行成像时，反射能量将聚焦在R/ξ=R₀的位置，而不是R₀所对应的位置x₀。而当采用黏声偏移的时候，反射体会成像在正确的位置x₀。但是在黏声逆时偏移的过程当中，由于$\exp \left( { - \frac{{\omega {R_0}}}{{2Q{v_0}\xi }}} \right) $项的作用，反向传播的检波器波场会进一步衰减，导致成像分辨率降低。

为了提高Q-LSRTM成像的分辨率并加速最小二乘迭代的收敛，类似于声学去模糊滤波器的构建方法，本文采用局部黏声去模糊滤波器。首先选定一个均匀分布的点散射体模型作为参考反射率模型，然后结合背景速度模型和Q模型合成参考黏声数据d_ref^Q，最后对黏声数据进行Q逆时偏移(Q-RTM)

$ \mathit{\boldsymbol{m}}_{{\rm{mig - ref}}}^{\rm{Q}} = \mathit{\boldsymbol{L}}_{\rm{Q}}^{\rm{T}}{\mathit{\boldsymbol{L}}_{\rm{Q}}}{\mathit{\boldsymbol{m}}_{{\rm{ref}}}} = \mathit{\boldsymbol{L}}_{\rm{Q}}^{\rm{T}}\mathit{\boldsymbol{d}}_{\rm{Q}}^{{\rm{ref}}} $

(15)

如1.1所述，在偏移成像的不同子区域内，通过匹配参考成像数据和参考反射率模型，得到局部去模糊滤波器，即

$ \left[ {{\mathit{\boldsymbol{m}}_{{\rm{ref}}}}} \right] = {\left[ {{\mathit{\boldsymbol{F}}_{\rm{Q}}}} \right]_i} \otimes {\left( {\mathit{\boldsymbol{m}}_{{\rm{mig - ref}}}^{\rm{Q}}} \right)_i} $

(16)

在实际操作当中，背景速度场选定为偏移速度场和偏移Q值模型，并在与实际相同的观测系统和地震子波情况下估计黏声去模糊滤波器，因此，这些黏声滤波器可以近似视为Hessian算子(L_Q^TL_Q)^-1的一种估计。

1.3 黏声最小二乘逆时偏移

基于标准线性固体(SLS)黏弹性模型一阶松弛机制，二维黏声介质波动程可写为^{[11, 42]}

$ \left\{ \begin{array}{l} \frac{{\partial P}}{{\partial t}} + K\left( {\tau + 1} \right)\nabla \cdot \mathit{\boldsymbol{v}} + {\gamma _p} = S\left( {{\mathit{\boldsymbol{x}}_s},t} \right)\\ \frac{{\partial \mathit{\boldsymbol{v}}}}{{\partial t}} + \frac{1}{\rho }\nabla P = 0\\ \frac{{\partial {\mathit{\boldsymbol{\gamma }}_p}}}{{\partial t}} + \frac{1}{{{\tau _\sigma }}}\left( {{\gamma _p} + \tau K\nabla \cdot \mathit{\boldsymbol{v}}} \right) = 0 \end{array} \right. $

(17)

式中：ρ为密度；P为压力；K为体积模量；S为震源子波；γ_p为记忆变量；v为质点速度向量；τ与品质因子Q有关，可表示为

$ \tau = \frac{{{\tau _\varepsilon }}}{{{\tau _\sigma }}} - 1 = \frac{2}{Q}\left( {\frac{1}{Q} + \sqrt {1 + \frac{1}{{{Q^2}}}} } \right) $

(18)

其中τ_σ和τ_ε分别表示应力和应变松弛时间。在Born近似的框架下，对体积模量K做一定扰动δK，则扰动场可以写为

$ \left\{ \begin{array}{l} \frac{{\partial \delta P}}{{\partial t}} + K\left( {\tau + 1} \right)\left( {\nabla \cdot \delta \mathit{\boldsymbol{v}}} \right) + \delta {\gamma _p}\\ \;\;\;\;\;\;\; = - \delta K\left( {\tau + 1} \right)\nabla \cdot \mathit{\boldsymbol{v}}\\ \frac{{\partial \delta \mathit{\boldsymbol{v}}}}{{\partial t}} + \frac{1}{\rho }\nabla \delta P = 0\\ \frac{{\partial \delta {\mathit{\boldsymbol{\gamma }}_p}}}{{\partial t}} + \frac{1}{{{\tau _\sigma }}}\left( {\delta {\gamma _p} + \tau K\nabla \cdot \delta \mathit{\boldsymbol{v}}} \right)\\ \;\;\;\;\;\;\;\; = - \frac{\tau }{{{\tau _\sigma }}}\delta K\nabla \cdot \mathit{\boldsymbol{v}} \end{array} \right. $

(19)

式中忽略了二阶扰动变量。若用压力和记忆变量格林函数G_P(x_r, t; x₀, 0)和G_γp(x_r, t; x₀, 0)，式(19)也可表示为

$ \begin{array}{l} \delta P\left( {{\mathit{\boldsymbol{x}}_r},t;{\mathit{\boldsymbol{x}}_s}} \right)\\ \;\;\;\; = \int_{V0} { - \delta K\left( {{\mathit{\boldsymbol{x}}_0}} \right)\left\{ {\left[ {\tau \left( {{\mathit{\boldsymbol{x}}_0}} \right) + 1} \right]{G_P}\left( {{\mathit{\boldsymbol{x}}_r},t;{\mathit{\boldsymbol{x}}_0},0} \right) * } \right.} \\ \;\;\;\;\;\;\nabla \cdot \mathit{\boldsymbol{v}}\left( {{\mathit{\boldsymbol{x}}_0},t;{\mathit{\boldsymbol{x}}_s}} \right) + \frac{{\tau \left( {{\mathit{\boldsymbol{x}}_0}} \right)}}{{{\tau _\sigma }\left( {{\mathit{\boldsymbol{x}}_0}} \right)}}{G_{{\tau _p}}}\left( {{\mathit{\boldsymbol{x}}_r},t;{\mathit{\boldsymbol{x}}_0},0} \right) * \\ \;\;\;\;\;\;\left. {\nabla \cdot \mathit{\boldsymbol{v}}\left( {{\mathit{\boldsymbol{x}}_0},t;{\mathit{\boldsymbol{x}}_s}} \right)} \right\}{\rm{d}}V \end{array} $

(20)

式中“*”表示时间方向的褶积运算。

在Q-LSRTM的框架下，式(20)的解等效于矩阵—向量运算

$ {\mathit{\boldsymbol{d}}_{\rm{Q}}} = {\mathit{\boldsymbol{L}}_{\rm{Q}}}{\mathit{\boldsymbol{m}}_0} $

(21)

式中：d_Q为Born模拟的黏声地震数据；L_Q为线性黏声模拟算子；m₀为地下介质反射率。用伴随状态法可推导式(19)的伴随方程^{[11, 12, 43]}为

$ \left\{ \begin{array}{l} \frac{{\partial q}}{{\partial t}} + \nabla \cdot \left( {\frac{\mathit{\boldsymbol{u}}}{\rho }} \right) = - \Delta d\left( {{\mathit{\boldsymbol{x}}_g},t;{\mathit{\boldsymbol{x}}_s}} \right)\\ \frac{{\partial \mathit{\boldsymbol{u}}}}{{\partial t}} + \nabla \left[ {K\left( {1 + \tau } \right)q} \right] + \nabla \left( {\frac{1}{{{\tau _\sigma }}}K\tau s} \right) = 0\\ \frac{{\partial s}}{{\partial t}} - \frac{s}{{{\tau _\sigma }}} - q = 0 \end{array} \right. $

(22)

式中：(q，u，s)为(P，v，γ_p)的伴随状态变量；Δd表示每次迭代中实际观测数据和正演数据之间的残差。

成像扰动δm与体积模量扰动δK近似线性相关，而δK可由背景场(通过式(19)计算)和伴随场(通过式(22)计算)之间的零时间延迟互相关得到，即

$ \delta m \approx \delta K = \int_0^T {\left( {1 + \tau } \right)\left( {\nabla \cdot \mathit{\boldsymbol{v}}} \right)q} + \frac{\tau }{{{\tau _\sigma }}}\left( {\nabla \cdot \mathit{\boldsymbol{v}}} \right)s{\rm{d}}t $

(23)

上式可用格林函数表示为

$ \begin{array}{l} \delta m \approx \delta K = \sum\limits_r {\sum\limits_s {\int_0^T {{\rm{d}}t\left\{ {\left[ {\tau \left( \mathit{\boldsymbol{x}} \right) + 1} \right]\nabla \cdot \mathit{\boldsymbol{v}}\left( {\mathit{\boldsymbol{x}},t;{\mathit{\boldsymbol{x}}_s}} \right) \times } \right.} } } \\ \;\;\;\;\;\;\;\left[ {{G_P}\left( {{\mathit{\boldsymbol{x}}_r}, - t;\mathit{\boldsymbol{x}},0} \right) * \Delta P\left( {{\mathit{\boldsymbol{x}}_r},t;{\mathit{\boldsymbol{x}}_s}} \right)} \right] + \\ \;\;\;\;\;\;\;\frac{{\tau \left( \mathit{\boldsymbol{x}} \right)}}{{{\tau _\sigma }\left( \mathit{\boldsymbol{x}} \right)}}\nabla \cdot \mathit{\boldsymbol{v}}\left( {\mathit{\boldsymbol{x}},t;{\mathit{\boldsymbol{x}}_s}} \right) \times \\ \;\;\;\;\;\;\;\left. {\left[ {{G_{{\tau _P}}}\left( {{\mathit{\boldsymbol{x}}_r}, - t;\mathit{\boldsymbol{x}},0} \right) * \Delta P\left( {{\mathit{\boldsymbol{x}}_r},t;{\mathit{\boldsymbol{x}}_s}} \right)} \right]} \right\}\\ \;\;\;\;\;\;\; = \sum\limits_s {\int {{\rm{d}}t\left\{ {\left[ {\tau \left( \mathit{\boldsymbol{x}} \right) + 1} \right]\nabla \cdot \mathit{\boldsymbol{v}}\left( {\mathit{\boldsymbol{x}},t;{\mathit{\boldsymbol{x}}_s}} \right)q\left( {\mathit{\boldsymbol{x}},t;{\mathit{\boldsymbol{x}}_s}} \right) + } \right.} } \\ \;\;\;\;\;\;\;\left. {\left[ {\frac{{\tau \left( \mathit{\boldsymbol{x}} \right)}}{{{\tau _\sigma }\left( \mathit{\boldsymbol{x}} \right)}}\nabla \cdot \mathit{\boldsymbol{v}}\left( {\mathit{\boldsymbol{x}},t;{\mathit{\boldsymbol{x}}_s}} \right)sq\left( {\mathit{\boldsymbol{x}},t;{\mathit{\boldsymbol{x}}_s}} \right)} \right]} \right\} \end{array} $

(24)

用矩阵可以表示为

$ {\mathit{\boldsymbol{m}}_{{\rm{mig}}}} = \mathit{\boldsymbol{L}}_{\rm{Q}}^{\rm{T}}{\mathit{\boldsymbol{d}}_{\rm{Q}}} $

(25)

1.4 预条件Q最小二乘逆时偏移

Q-LSRTM的误差泛函^[29]为

$ \varepsilon = \frac{1}{2}{\left\| {{\mathit{\boldsymbol{L}}_{\rm{Q}}}{\mathit{\boldsymbol{m}}^{\left( k \right)}} - \mathit{\boldsymbol{d}}_{\rm{Q}}^{{\rm{obs}}}} \right\|^2} $

(26)

式中：d_Q^obs表示包含衰减损失的观测数据；m^(k)表示第k次迭代的偏移成像。在高斯牛顿框架下，误差泛函梯度为

$ {\mathit{\boldsymbol{g}}^{\left( k \right)}} = \left( {\mathit{\boldsymbol{L}}_{\rm{Q}}^{\rm{T}}{\mathit{\boldsymbol{L}}_{\rm{Q}}}} \right)\Delta {\mathit{\boldsymbol{m}}^{\left( k \right)}} = \mathit{\boldsymbol{L}}_{\rm{Q}}^{\rm{T}}\left[ {{\mathit{\boldsymbol{L}}_{\rm{Q}}}{\mathit{\boldsymbol{m}}^{\left( k \right)}} - \mathit{\boldsymbol{d}}_{\rm{Q}}^{{\rm{obs}}}} \right] $

(27)

从上式可以看出，每次迭代梯度g^(k)是期望模型更新Δm^(k)的模糊版本，模糊算子为黏声偏移格林函数L_Q^TL_Q。于是，模型梯度可通过黏声去模糊滤波器近似[L_Q^TL_Q]^-1≈F_Q得到，从而预条件化梯度被用于迭代更新方程

$ {\mathit{\boldsymbol{m}}^{\left( {k + 1} \right)}} = {\mathit{\boldsymbol{m}}^{\left( k \right)}} - \alpha {\mathit{\boldsymbol{F}}_{\rm{Q}}} * \left[ {\mathit{\boldsymbol{L}}_{\rm{Q}}^{\rm{T}}\left( {{\mathit{\boldsymbol{L}}_{\rm{Q}}}{\mathit{\boldsymbol{m}}^{\left( k \right)}} - \mathit{\boldsymbol{d}}_{\rm{Q}}^{{\rm{obs}}}} \right)} \right] $

(28)

式中α是更新步长。

2 Marmousi Ⅱ模型测试

为了说明本文提出的预条件Q-LSRTM的有效性，应用Marmousi Ⅱ模型进行测试。图 2a和图 2b分别为正演速度模型和Q值模型，其中Q只取两个值，Q=20表示衰减层，Q=10000可视为弹性介质。整个模型的网格点数为801×351，模型纵横向采样间隔均为10m。在合成黏声地震记录时，炮间隔为50m，共计150炮均匀分布在模型表面，道间距为10m，震源采用主频为15Hz的Ricker子波。采集系统位于地表，且排列长度固定为整个模型的长度。计算中采用时间二阶精度、空间八阶精度的时间域交错网格有限差分法及C-PML边界条件。本文共用到两套数据，一套是基于声波方程的正演结果，另一套是基于一阶松弛机制黏声波动方程的正演结果。对真实速度模型和Q值模型进行二维平滑，得到偏移速度模型(图 2c)和偏移Q值模型(图 2d)。分别执行声波LSRTM、Q-LSRTM、预条件Q-LSRTM三种最小二乘逆时偏移，迭代次数均为20次，其中第1次迭代结果即为声波逆时偏移(RTM)、Q-RTM和预条件Q-RTM的成像结果。

图 3为声波数据声波RTM和LSRTM成像结果，并以此作为基准对比数据。图 4为黏声数据不同方法的偏移结果，可以看出：①利用传统RTM(图 4a)和LSRTM(图 4b)对黏声数据处理获得的成像结果都不能有效恢复衰减层及其下方的反射体振幅，只能粗略体现地下构造形态；②Q-RTM成像结果(图 4c)分辨率比传统RTM(图 4a)还要低，这是由于伴随算子L_Q^T的衰减性质决定的；③与传统LSRTM成像结果相比(图 4b)，Q-LSRTM(图 4d)明显提高了分辨率，但与基准声波LSRTM成像结果(图 3b)相比，分辨率还是相对较低，其原因还是由于伴随算子L_Q^T具有衰减性质；④预条件Q-RTM成像结果(图 4e)增强了衰减层及其下方反射体振幅，但整体噪声有所增强，且最下方背斜同相轴归位不准，这与参考点散射体模型的选取有关；⑤预条件Q-LSRTM成像结果(图 4f)与基准声波LSRTM成像结果(图 3b)较为接近，具有较高的分辨率，与Q-LSRTM成像结果(图 4d)相比，其深层反射体的振幅更加均衡。

图 5、图 6分别为声波数据LSRTM(图 3b)、黏声数据LSRTM(图 4b)、Q-LSRTM(图 4d)、预条件Q-LSRTM(图 4f)成像结果蓝色、红色矩形框局部放大显示，可以看出：预条件Q-LSRTM提高了黏声数据的成像精度和分辨率(蓝色箭头所指)；采用LSRTM对黏声数据进行成像时，不仅存在成像分辨率低的问题，而且由于介质黏弹性导致地震波振幅和相位的畸变，还会导致深层反射体归位出现错乱(图 6b红色箭头所指)。相对而言，尽管Q-LSRTM能有效归位主要反射界面，但分辨率较低，且不能呈现构造细节，预条件Q-LSRTM成像与基准声波LSRTM结果比较接近，分辨率较高。

图 7为声波数据LSRTM(图 5a)、黏声数据LSRTM(图 5b)、Q-LSRTM(图 5c)、预条件Q-LSRTM(图 5d)在x=5.0km处成像结果的中深层归一化波数谱对比，可以看出，传统LSRTM和Q-LSRTM的主波数都相对偏低，只有预条件Q-LSRTM的波数谱与基准声波LSRTM接近，也说明预条件Q-LSRTM提高了黏声数据的成像分辨率。图 8为声波数据LSRTM(图 6a)、黏声数据LSRTM(图 6b)、Q-LSRTM(图 6c)、预条件Q-LSRTM(图 6d)在x=2.6km处成像结果的深层归一化波数谱对比，可以看出，黏声数据的声波LSRTM波数谱严重偏离基准数据谱，Q-LSRTM有所改善，预条件Q-LSRTM能恢复出高波数细节信息，最接近于基准数据。

图 9为黏声数据声波LSRTM、Q-LSRTM、预条件Q-LSRTM(预条件20次)、预条件Q-LSRTM(预条件6次)四种情况下数据残差随着迭代次数的变化曲线，可以看出，四种偏移成像结果的精度是逐渐升高的。就收敛速率来看，预条件Q-LSRTM第4次迭代和第8次迭代时的数据残差接近于Q-LSRTM第8次和第20次迭代的数据残差，因此预条件Q-LSRTM能减少约50%的计算量，尤其在早期的几次迭代当中体现的更为明显。预条件Q-LSRTM收敛速率提升的原因在于黏声去模糊滤波器是Hessian算子逆的一种有效近似估计，尽管需要预先计算去模糊滤波器，但计算成本还是要比Q-LSRTM低得多。兼顾到计算精度和计算效率，参考反射率模型中点散射体的布置应该相对稀疏，由此会引入一些全局噪声，如果预条件在第6次迭代后终止，数据残差会进一步减小，相应地成像精度也有所提高。

为了测试预条件Q-LSRTM对Q值的敏感性，分别应用Q=20、40、80、100四种模型进行偏移成像，结果如图 4f、图 10所示，图 11为对应的参考反射率模型的黏声PSF，可以看出：随着Q值的增大，衰减层及其下部的反射体的成像振幅逐渐减弱，成像分辨率逐渐降低；相反地，当Q值越大时，衰减越弱，参考反射率模型中的点散射体偏移响应即PSF越强，从而用其构建黏声去模糊滤波器时，考虑到的衰减效应越弱、成像精度越低。因此，一个相对精确的偏移Q值模型是预条件Q-LSRTM取得良好成像效果的前提。

3 结论与认识

(1) 在传统Q-LSRTM成像当中，数据残差的反传过程也是衰减的，会导致成像分辨率降低；而基于PSF的预条件Q-LSRTM方法采用黏声去模糊滤波器补偿由于地下强衰减导致的振幅和分辨率降低，能得到较高的成像精度和更均衡的成像振幅。

(2) 由于地下模型是未知的，构建黏声去模糊滤波器时一般选用均匀分布的点散射体模型。基于偏移速度场和Q值模型，采用和实际一致的观测系统和地震子波，通过正演模拟和偏移成像获得黏声PSF，最终采用匹配滤波的思想构建黏声去模糊滤波器。兼顾到计算精度和计算效率，参考反射率模型中点散射体密度和局部窗口长度需根据构造复杂度确定。

(3) 与Q-LSRTM相比，预条件Q-LSRTM的收敛效率得到大幅度提升，在前几次迭代中更为明显。为了获得更高的成像精度和计算效率，预条件Q-LSRTM只在前几次迭代中使用为宜。就计算成本而言，尽管预条件Q-LSRTM需要预先计算黏声去模糊滤波器，但计算耗时仍然比Q-LSRTM小得多。

(4) 与Q-LSRTM类似，为了获得较好的成像效果，预条件Q-LSRTM也需要一个相对精确的Q值模型估计，同时也必须有一个相对精确的偏移速度模型。

(5) 本文的方法可拓展到黏弹各向异性等更复杂介质的偏移成像当中，考虑到计算量，建议采用并行加速计算。