随着地震勘探开发的深入发展, 复杂地质条件下的地震勘探已经成为目前学术界与工业界的研究重点与热点。地球介质广泛存在粘滞性, 地震波在地下传播时会发生衰减, 深层能量呈指数形式下降, 如若忽略掉粘弹介质的吸收衰减作用, 会造成成像结果偏差大, 频散严重, 深层能量不足且保幅性差, 进而严重影响后续的反演解释结果。

针对粘介质的吸收衰减作用, 人们首先采用反Q滤波方法来补偿高频能量, 但是该方法只适用于简单介质。随着偏移技术的发展, 利用Q偏移方法来恢复因吸收衰减作用而损失的能量逐渐成为研究的热点。如RIBODETTI等^[1]提出了基于射线理论的Q偏移方法, 主要是利用射线理论做高频近似。随着逆时偏移技术的发展, 基于波动方程理论的Q偏移也取得了长足的进步。CAUSSE等^[2]将粘声介质引入逆时偏移中, 利用时间域波动方程实现了简单的粘声逆时偏移。ZHU等^[3]在成像之前先进行相位校正, 再做基于粘声介质的逆时偏移, 较好地解决了相位异常的问题。ZHANG等^[4]推导了基于粘声介质常Q模型的微分方程, 有效地补偿了深层能量衰减与相位频散影响。同时, 针对高频不稳定的情况引入了规则化算子, 有效解决了反向延拓过程中的不稳定问题。与单程波偏移方法相比, 逆时偏移有保幅方面的优势。

传统的逆时偏移成像只是简单的数据共轭转置, 并不是正演的逆, 而最小二乘偏移是在非线性问题比较弱的情况下, 利用Born近似对非线性问题局部寻优迭代求解的一种反演方法。但是最小二乘偏移的计算量较大, 需要大量的存储空间, 因此其发展较为缓慢。很多学者在反演思想的框架下, 将逆时偏移与最小二乘偏移思想相结合。TARANTOLA^[5]最早提出了最小二乘偏移思想, LEBRAS等^[6]、LAMBARE等^[7]在其基础上进行了完善, 使用最速下降法反演相对于背景速度的速度扰动。NEMETH等^[8]首先将最小二乘算法引入地震偏移以求取地下反射系数, 然后将Kirchhoff最小二乘偏移方法应用于不完整的反射地震数据偏移。目前, 许多学者将单程波方程和逆时偏移与最小二乘理论相结合, 实现了基于波动方程的最小二乘偏移^[9-16]。黄建平等^[17]将最小二乘偏移方法应用于碳酸盐岩裂缝型储层预测, 详细分析了最小二乘偏移理论在碳酸盐岩成像中的优势, 并将其应用于地下裂缝模型测试。刘玉金等^[18]通过引入预条件算子实现了扩展最小二乘逆时偏移, 进一步提高了最小二乘逆时偏移的精度和确定性。QU等^[19]实现了弹性波最小二乘逆时偏移, 更加真实地模拟了地下实际情况, 取得了较好的成像效果。

本文在反演思想的框架下, 引入最小二乘理论, 给出了基于广义标准线性固体模型(GSLS)的粘声介质波动方程, 并推导了对应的反偏移算子与伴随算子, 将粘声介质逆时偏移与最小二乘理论结合, 实现了基于粘声衰减补偿的最小二乘逆时偏移(LSRTM)。通过层状模型与Marmousi模型试算验证了该方法的正确性及有效性。

1 方法原理 1.1 粘声介质波动方程

基于GSLS模型的波动方程可以写成如下形式:

$ \begin{array}{l} \frac{{{\partial ^2}p}}{{\partial {t^2}}} = \frac{\partial }{{\partial t}}\left\{ {{M_R}\left[ {1 - \sum\limits_{i = 1}^I {\left( {1 - \frac{{{\tau _{\varepsilon i}}}}{{{\tau _{\sigma i}}}}} \right){{\rm{e}}^{ - t/{t_{\sigma i}}}}} } \right]H\left( t \right)} \right\} * \\ \;\;\;\;\;\;\;\;\;\nabla \cdot \left( {\frac{1}{\rho }\nabla \rho } \right) + f\left( t \right) \end{array} $

(1)

式中:I是标准线性固体的个数, ▽是梯度算子, ▽·是散度算子, *是卷积, M_R是松弛模量, τ_εi和τ_σi是松弛时间, H(t)为单位阶跃函数, f(t)为震源项。

采用(1)式的交错网格有限差分形式可以进行波动方程模拟。基于一阶速度-应力交错网格的有限差分模拟是目前使用较多的交错网格模拟方式, 交错网格可以有效代替普通矩形网格, 更好地模拟地下真实情况, 提高正演模拟的精度。

1.2 粘声介质逆时偏移所对应的反偏移算子

最小二乘逆时偏移的关键步骤是反偏移, 主要是通过背景波场与慢度扰动的互相关产生所需要的反射波。反偏移并不像一般的正演过程那样需要界面信息, 可以在光滑的背景速度场中获取反射波。

根据波场叠加原理, 有:

$ p\left( {x,t} \right) = {p_0}\left( {x,t} \right) + {p_{\rm{s}}}\left( {x,t} \right) $

(2)

这里, p₀为背景波场, p_s为扰动波场。背景波场的传播方程与普通波场传播方程一致, 但是传播速度存在较大区别。背景波场的速度场是平滑速度场, 方程表示为:

$ \begin{array}{l} \frac{1}{{v_0^2}}\frac{{{\partial ^2}{p_0}}}{{\partial {t^2}}} = \left( {1 + \tau } \right)\nabla \cdot \nabla {p_0} - \frac{\tau }{{{\tau _\sigma }}}\left[ {{{\rm{e}}^{ - \tau /{\tau _\sigma }}}H\left( t \right)} \right] * \\ \;\;\;\;\;\;\;\;\;\;\;\;\;\;\left( {\nabla \cdot \nabla {p_0}} \right) + f\left( t \right) \end{array} $

(3)

定义慢度扰动为:

$ m\left( x \right) = \frac{1}{{{v^2}}} - \frac{1}{{v_0^2}} $

(4)

式中:v为全波场速度, v₀为背景波场速度。扰动波场的传播方程也与普通方程类似, 而且速度场也是平滑后的速度场。扰动波场传播与背景波场传播最大的不同是震源不同, 扰动波场的震源不再是传统的地震子波震源, 而是由背景波场与慢度扰动的乘积代替:

$ \begin{array}{*{20}{c}} {\frac{1}{{{v^2}}}\frac{{{\partial ^2}{p_{\rm{s}}}}}{{\partial {t^2}}} = \left( {1 + \tau } \right)\nabla \cdot \nabla {p_{\rm{s}}} - \frac{\tau }{{{\tau _\sigma }}}\left[ {{{\rm{e}}^{ - \tau /{\tau _\sigma }}}H\left( t \right)} \right] * }\\ {\left( {\nabla \cdot \nabla {p_{\rm{s}}}} \right) - m\left( x \right)\frac{{{\partial ^2}{p_0}}}{{\partial {t^2}}}} \end{array} $

(5)

为了简化上述方程, 对该方程应用Born近似, 得到基于Born近似的粘声介质最小二乘逆时偏移反偏移算子:

$ \left\{ \begin{array}{l} \left\{ {\frac{1}{{v_0^2}}\frac{{{\partial ^2}}}{{\partial {t^2}}} - \left( {1 + \tau } \right)\nabla \cdot \nabla + \frac{\tau }{{{\tau _\sigma }}}\left[ {{{\rm{e}}^{ - \tau /{\tau _\sigma }}}H\left( t \right)} \right] * } \right.\\ \;\;\;\;\;\;\;\left. {\left( {\nabla \cdot \nabla } \right)} \right\}{p_0}\left( {x,t} \right) = f\left( t \right)\\ \left\{ {\frac{1}{{v_0^2}}\frac{{{\partial ^2}}}{{\partial {t^2}}} - \left( {1 + \tau } \right)\nabla \cdot \nabla + \frac{\tau }{{{\tau _\sigma }}}\left[ {{{\rm{e}}^{ - \tau /{\tau _\sigma }}}H\left( t \right)} \right] * } \right.\\ \;\;\;\;\;\;\;\left. {\left( {\nabla \cdot \nabla } \right)} \right\}{p_{\rm{s}}}\left( {x,t} \right) = m\left( x \right)\frac{{{\partial ^2}{p_0}\left( {x,t} \right)}}{{{\partial ^2}t}} \end{array} \right. $

(6)

扰动波场p_s(x, t)由两个过程计算得到:一是震源f(t)和背景速度v₀产生p₀(x, t); 二是利用m(x)[∂²p₀(x, t)]/∂²t和背景速度v₀产生p_s(x, t)。从公式(6)可以看出, 背景场中没有扰动场; 扰动场中只存在一次扰动场。但Born近似依赖于地下尺度大小, 当存在小尺度体时, Born近似效果较好。这是当前最小二乘叠前深度偏移成像的基础。

1.3 粘声介质逆时偏移对应的伴随算子

粘声介质逆时偏移包括两个过程:波的正向传播与波的反向传播。随着时间的增加, 粘声介质地震波在地下的正向传播是稳定的。在波的反推过程中, 随着时间的减少往往会出现不稳定问题。解决该问题的方法主要有两种:一种是在反传过程中引入规则化算子, 压制反传过程中的高频震荡。另一种是通过正演波动方程的伴随方程实现波场的逆时延拓, 有效地避免反推过程中的不稳定问题。本文采用第二种方法解决不稳定问题。

去掉震源项, 正演算子T波动方程可以写成:

$ \begin{array}{l} \frac{1}{{{v^2}}}\frac{{{\partial ^2}p}}{{\partial {t^2}}} - \left( {1 + \tau } \right)\nabla \cdot \nabla p + \frac{\tau }{{{\tau _\sigma }}}\left[ {{{\rm{e}}^{ - \tau /{\tau _\sigma }}}H\left( t \right)} \right] * \\ \;\;\;\;\;\;\;\left( {\nabla \cdot \nabla p} \right) = 0 \end{array} $

(7)

令正演算子T的伴随算子为T^*, 两者的关系满足 <Tx, y>= <x, T^*y>, 则:

$ \begin{array}{l} < \frac{1}{{{v^2}}}\frac{{{\partial ^2}p}}{{\partial {t^2}}} - \left( {1 + \tau } \right)\nabla \cdot \nabla p + \frac{\tau }{{{\tau _\sigma }}}\left[ {{{\rm{e}}^{ - \tau /{\tau _\sigma }}}H\left( t \right)} \right] * \\ \;\;\;\left( {\nabla \cdot \nabla p} \right),u > = < p,\frac{1}{{{v^2}}}\frac{{{\partial ^2}u}}{{\partial {t^2}}} - \nabla \cdot \nabla \left( {1 + \tau } \right)u + \\ \;\;\;\frac{\tau }{{{\tau _\sigma }}}\left[ {{{\rm{e}}^{ - t/{\tau _\sigma }}}H\left( { - t} \right)} \right] * u > \end{array} $

(8)

因此, 正演的伴随算子T^*也是一个二阶微分方程, 表示为:

$ \frac{1}{{{v^2}}}\frac{{{\partial ^2}p}}{{\partial {t^2}}} = \nabla \cdot \nabla \left( {1 + \tau } \right)p - \left( {\nabla \cdot \nabla \bar r} \right) $

(9)

这里,

$ \bar r = \frac{\tau }{{{\tau _\sigma }}}\left[ {{{\rm{e}}^{ - t/{\tau _\sigma }}}H\left( { - t} \right)} \right] * p $

(10)

对公式(10)两边同时对时间求导, 则:

$ \frac{{\partial \bar r}}{{\partial t}} = \frac{\tau }{{{\tau _\sigma }}}\left[ { - p + \frac{1}{{{\tau _\sigma }}}{e^{ - t/{\tau _\sigma }}}H\left( { - t} \right)} \right] * p $

(11)

为了消除公式(11)中褶积的影响, 运用递归卷积方法得到的一阶线性微分方程:

$ \frac{{\partial \hat r}}{{\partial t}} = - \frac{\tau }{{{\tau _\sigma }}}p + \frac{1}{{{\tau _\sigma }}}\hat r $

(12)

公式(9)和公式(12)构成了粘声介质反向传播方程。为了便于粘声介质微分方程进行交错网格数值模拟, 引入u_x和u_z, 将其空间二阶微分方程降为空间一阶微分方程, 对空间微分进行交错网格求解, 再建立复频移完美匹配层(PML)边界方程, 进行数值模拟。

1.4 粘声介质最小二乘理论

最小二乘逆时偏移的主要思想是利用波动方程来模拟反射系数的响应, 得到正演记录, 然后将其与输入的实际地震数据对比, 如果存在误差, 就采用最小二乘方法校正, 直至误差减小到一个合理的范围内。

定义最小二乘优化的目标函数:

$ f\left( \mathit{\boldsymbol{m}} \right) = \frac{1}{2}\left\| {\mathit{\boldsymbol{Lm}} - \mathit{\boldsymbol{d}}} \right\|_{{{\rm{H}}_2}}^2 + \frac{a}{2}\left\| \mathit{\boldsymbol{m}} \right\|_{{{\rm{H}}_1}}^2 $

(13)

式中:m代表反射系数模型, 即成像结果; L为Born近似下的反偏移算子; d为实际输入数据(野外炮记录); a为正则化参数; H₁和H₂代表 1范数和2范数。当目标函数f(m)最小化时, 我们认为此时求得的m就是最优解。求解(13)式的首要问题是计算目标函数关于参数变化的导数, 即目标函数的梯度, 本文采用以下变分方法。

对任意的介质扰动δm和扰动步长t, 总存在如下目标函数变分:

$ \begin{array}{l} f\left( {\mathit{\boldsymbol{m}} + t\delta \mathit{\boldsymbol{m}}} \right) = \frac{1}{2}\left\| {\mathit{\boldsymbol{L}}\left( {\mathit{\boldsymbol{m}} + t\delta \mathit{\boldsymbol{m}}} \right) - \mathit{\boldsymbol{d}}} \right\|_{{{\rm{H}}_2}}^2 + \frac{a}{2} \cdot \\ \left\| {\mathit{\boldsymbol{m}} + t\delta \mathit{\boldsymbol{m}}} \right\|_{{{\rm{H}}_1}}^2 = \frac{1}{2} < \mathit{\boldsymbol{Lm}} - \mathit{\boldsymbol{d}} + t\mathit{\boldsymbol{L}}\delta \mathit{\boldsymbol{m}},\mathit{\boldsymbol{Lm}} - \mathit{\boldsymbol{d}} + \\ t\mathit{\boldsymbol{L}}\delta \mathit{\boldsymbol{m}}{ > _{{{\rm{H}}_2}}} + \frac{a}{2} < \mathit{\boldsymbol{m}} + t\delta \mathit{\boldsymbol{m}},\mathit{\boldsymbol{m}} + t\delta \mathit{\boldsymbol{m}}{ > _{{{\rm{H}}_1}}} = \frac{1}{2} \cdot \\ < \mathit{\boldsymbol{Lm}} - \mathit{\boldsymbol{d}},\mathit{\boldsymbol{Lm}} - \mathit{\boldsymbol{d}}{ > _{{{\rm{H}}_2}}} + \frac{t}{2}\left( { < \mathit{\boldsymbol{Lm}} - \mathit{\boldsymbol{d}},\mathit{\boldsymbol{L}}\delta \mathit{\boldsymbol{m}}{ > _{{{\rm{H}}_2}}} + } \right.\\ \left. { < \mathit{\boldsymbol{L}}\delta \mathit{\boldsymbol{m}},\mathit{\boldsymbol{Lm}} - \mathit{\boldsymbol{d}}{ > _{{{\rm{H}}_2}}}} \right) + \frac{{{t^2}}}{2} < \mathit{\boldsymbol{L}}\delta \mathit{\boldsymbol{m}},\mathit{\boldsymbol{L}}\delta \mathit{\boldsymbol{m}}{ > _{{{\rm{H}}_2}}} + \\ \frac{a}{2} < \mathit{\boldsymbol{m}},\mathit{\boldsymbol{m}}{ > _{{{\rm{H}}_1}}} + \frac{a}{2}t\left( { < \mathit{\boldsymbol{m}},\delta \mathit{\boldsymbol{m}}{ > _{{{\rm{H}}_1}}} + < \delta \mathit{\boldsymbol{m}},} \right.\\ \left. {\mathit{\boldsymbol{m}}{ > _{{{\rm{H}}_1}}}} \right) + \frac{{a{t^2}}}{2} < \delta \mathit{\boldsymbol{m}},\delta \mathit{\boldsymbol{m}}{ > _{{{\rm{H}}_1}}} = \frac{1}{2}\left( {\mathit{\boldsymbol{Lm}} - \mathit{\boldsymbol{d}},\mathit{\boldsymbol{Lm}} - } \right.\\ \left. {\mathit{\boldsymbol{d}}{ > _{{{\rm{H}}_2}}} + a < \mathit{\boldsymbol{m}},\mathit{\boldsymbol{m}}{ > _{{{\rm{H}}_1}}}} \right) + \frac{t}{2}\left( { < \mathit{\boldsymbol{Lm}} - \mathit{\boldsymbol{d}},\mathit{\boldsymbol{L}}\delta \mathit{\boldsymbol{m}}} \right.\\ \left. {{ > _{{{\rm{H}}_2}}} + < \mathit{\boldsymbol{L}}\delta \mathit{\boldsymbol{m}},\mathit{\boldsymbol{Lm}} - \mathit{\boldsymbol{d}}{ > _{{{\rm{H}}_2}}}} \right) + \frac{{at}}{2}\left( { < \mathit{\boldsymbol{m}},\delta \mathit{\boldsymbol{m}}{ > _{{{\rm{H}}_1}}} + } \right.\\ \left. { < \delta \mathit{\boldsymbol{m}},\mathit{\boldsymbol{m}}{ > _{{{\rm{H}}_1}}}} \right) + \frac{{{t^2}}}{2}\left( {\left\| {\mathit{\boldsymbol{L}}\delta \mathit{\boldsymbol{m}}} \right\|_{{{\rm{H}}_2}}^2 + a\left\| {\delta \mathit{\boldsymbol{m}}} \right\|_{{{\rm{H}}_{\rm{1}}}}^2} \right) \end{array} $

(14)

对上式定义的变分求偏导数:

$ \begin{array}{l} \frac{{\partial f\left( \mathit{\boldsymbol{m}} \right)}}{{\partial t}}\left| {_{t = 0}} \right. = \left\{ {\left( { < \mathit{\boldsymbol{Lm}} - \mathit{\boldsymbol{d}},\mathit{\boldsymbol{L}}\delta \mathit{\boldsymbol{m}}{ > _{{{\rm{H}}_2}}} + } \right.} \right.\\ a < \mathit{\boldsymbol{m}},\delta \mathit{\boldsymbol{m}}{ > _{{{\rm{H}}_1}}} + t\left( {\left\| {\mathit{\boldsymbol{L}}\delta \mathit{\boldsymbol{m}}} \right\|_{{{\rm{H}}_2}}^2 + } \right.\\ \left. {\left. {a\left\| {\delta \mathit{\boldsymbol{m}}} \right\|_{{{\rm{H}}_{\rm{1}}}}^2} \right)} \right\}\left| {_{t = 0}} \right. = \left( { < \mathit{\boldsymbol{Lm}} - \mathit{\boldsymbol{d}},\mathit{\boldsymbol{L}}\delta \mathit{\boldsymbol{m}}{ > _{{{\rm{H}}_2}}} + } \right.\\ \left. {a < \mathit{\boldsymbol{m}},\delta \mathit{\boldsymbol{m}}{ > _{{{\rm{H}}_1}}}} \right) = \left[ { < {\mathit{\boldsymbol{L}}^ * }\left( {\mathit{\boldsymbol{Lm}} - \mathit{\boldsymbol{d}}} \right),\delta \mathit{\boldsymbol{m}}{ > _{{{\rm{H}}_1}}} + } \right.\\ \left. {a < \mathit{\boldsymbol{m}},\delta \mathit{\boldsymbol{m}}{ > _{{{\rm{H}}_1}}}} \right] = < {\mathit{\boldsymbol{L}}^ * }\left( {\mathit{\boldsymbol{Lm}} - \mathit{\boldsymbol{d}}} \right)a\mathit{\boldsymbol{m}},\delta \mathit{\boldsymbol{m}}{ > _{{{\rm{H}}_1}}} = \\ < \left( {{\mathit{\boldsymbol{L}}^ * }\mathit{\boldsymbol{L}} + af} \right)\mathit{\boldsymbol{m}} - {\mathit{\boldsymbol{L}}^ * }\mathit{\boldsymbol{\hat d}},\delta \mathit{\boldsymbol{m}}{ > _{{{\rm{H}}_1}}} = \\ < \mathit{\boldsymbol{g}}f\left( \mathit{\boldsymbol{m}} \right),\delta \mathit{\boldsymbol{m}}{ > _{{{\rm{H}}_1}}} \end{array} $

(15)

推导出g(m)=(L^*L+af)m－L^*$\mathit{\boldsymbol{\hat d}}$, 即误差函数相对于模型变量的梯度。式中, L^*是波场正向传播的伴随算子, 代表波场的反向传播。一般利用该梯度构造最速下降迭代方法, 以求解Born近似后的线性问题的解。迭代公式为:

$ \begin{array}{l} {\mathit{\boldsymbol{m}}^{\left( {k + 1} \right)}} = {\mathit{\boldsymbol{m}}^{\left( k \right)}} - \mathit{\boldsymbol{mg}}\left( \mathit{\boldsymbol{m}} \right) = {\mathit{\boldsymbol{m}}^{\left( k \right)}} - \mathit{\boldsymbol{mL}} * \\ \;\;\;\;\;\;\;\left( {\mathit{\boldsymbol{Lm}} - {\mathit{\boldsymbol{d}}^{{\rm{obs}}}}} \right) = {\mathit{\boldsymbol{m}}^{\left( k \right)}} - \mathit{\boldsymbol{mL}} * \left( {\mathit{\boldsymbol{L}}{\mathit{\boldsymbol{m}}^{\left( k \right)}} - {\mathit{\boldsymbol{d}}^{{\rm{obs}}}}} \right) = \\ \;\;\;\;\;\;\;{\mathit{\boldsymbol{m}}^{\left( k \right)}} - \mathit{\boldsymbol{m}}{\mathit{\boldsymbol{L}}^ * }\left( {{\mathit{\boldsymbol{d}}^{{\rm{cal}}}} - {\mathit{\boldsymbol{d}}^{{\rm{obs}}}}} \right) \end{array} $

(16)

该最优化问题也可以采用共轭梯度法来求解, 实现思想与最速下降法类似, 就是以模型m₀为初始模型, 建立迭代公式, 求取梯度与步长, 最后检测误差是否在允许范围内。在目标函数对模型参数m求导时, 误差的梯度为:

$ \mathit{\boldsymbol{g}}\left( \mathit{\boldsymbol{m}} \right) = \frac{{\partial f\left( \mathit{\boldsymbol{m}} \right)}}{{\partial \mathit{\boldsymbol{m}}}} = {\mathit{\boldsymbol{L}}^ * }\left( {\mathit{\boldsymbol{Lm}} - \mathit{\boldsymbol{d}}} \right) $

(17)

其迭代公式为:

$ {\mathit{\boldsymbol{m}}^{\left( {k + 1} \right)}} = {\mathit{\boldsymbol{m}}^{\left( k \right)}} - \alpha \mathit{\boldsymbol{g}}\left( \mathit{\boldsymbol{m}} \right) $

(18)

$ \alpha = \frac{{{{\left( {{\mathit{\boldsymbol{g}}^k}} \right)}^{\rm{T}}}{\mathit{\boldsymbol{g}}^k}}}{{{{\left( {\mathit{\boldsymbol{L}}{\mathit{\boldsymbol{g}}^k}} \right)}^{\rm{T}}}\mathit{\boldsymbol{L}}{\mathit{\boldsymbol{g}}^k}}} $

(19)

式中:α为迭代步长。

在最小二乘叠前深度偏移中, 一般令梯度等于零(实际梯度不等于零)时, 可以得到线性方程:

$ \left( {{\mathit{\boldsymbol{L}}^ * }\mathit{\boldsymbol{L}} + af} \right)\mathit{\boldsymbol{m}} = {\mathit{\boldsymbol{L}}^ * }\mathit{\boldsymbol{d}} $

(20)

该方程求解的核心问题是求解Hessian矩阵的逆, 但是Hessian矩阵的逆无法显式求解, 因此难以得到其精确解。由于迭代求解的算法不是直接显式求解Hessian矩阵的逆, 而是通过隐式求解的方法得到, 因此最小二乘的方法被广泛应用于偏移与全波形反演中。

具体实现流程如图 1所示。

2 模型试算 2.1 层状模型试算

采用如图 2所示水平层状介质模型(Q均值为101)进行了测试。模型纵、横向采样点数为150×150, 纵、横向网格间距为10m。震源为主频20Hz的雷克子波, 共50炮, 炮间距30m, 炮点分布在网格点上。每炮150个检波器, 采用全接收排列, 检波器间隔10m, 采用粘声介质交错网格正演模拟得到炮记录。图 3给出了声波最小二乘逆时偏移第1次迭代结果、粘声波最小二乘逆时偏移第1次迭代结果、声波最小二乘逆时偏移第40次迭代结果以及粘声波最小二乘逆时偏移第40次迭代结果。

从图 3可以看出, 无论是逆时偏移还是最小二乘逆时偏移, 采用粘声波方程的效果明显好于采用普通声波方程的效果。由于采样射线不足, 常规声波逆时偏移对边界的刻画存在一定不足(图 3a, 图 3c), 但粘声波逆时偏移的结果分辨率较高, 对浅层噪声压制效果较好(图 3b, 图 3d)。由图 3c和图 3d可以看出, 经过40次迭代, 声波最小二乘逆时偏移与粘声波最小二乘逆时偏移效果都有了较大提高, 同相轴能量加强, 能量分布更加均衡。但是声波最小二乘逆时偏移存在较强的浅层噪声, 边界处成像效果较差, 深层能量不足, 同相轴不够清晰。而粘声波最小二乘逆时偏移能明显压制浅层噪声, 成像分辨率更高。特别是在第二层(深层)界面, 同相轴更加清晰, 分辨率较高, 能量更加均衡, 有效补偿了吸收衰减损失的能量。

图 4给出了最小二乘逆时偏移数据残差随迭代次数变化的曲线。可以看出, 随着迭代次数的增加, 残差稳定下降, 直至最后收敛。粘声最小二乘逆时偏移由于对模型粘性吸收衰减进行了补偿, 收敛得更快。图 5, 图 6, 图 7分别对比了x=250m, 750m, 1250m处粘声波最小二乘逆时偏移第1次和第40次迭代结果、声波最小二乘逆时偏移第1次和第40次迭代结果与真实反射系数。由于观测系统分布及覆盖次数等问题, 不同位置能量可能存在一定差异, 图 5和图 7能量较弱, 图 6能量较强。对比图 5, 图 6和图 7可以看出, 第1次迭代采用声波和粘声波方程偏移的结果均与真实反射系数差距较大, 但是粘声波逆时偏移明显能量更强, 效果更好, 更加接近真实反射系数。经过40次迭代以后, 声波和粘声波偏移效果都有提高, 但是很明显, 粘声波最小二乘逆时偏移已经接近真实反射系数, 能量更加均衡, 尤其是第二层(深层)界面, 粘声波最小二乘逆时偏移能与真实反射系数接近, 保幅性较好。说明该方法能很好地消除地层吸收等因素对地震波振幅及相位的影响, 有效提高成像分辨率, 实现保幅成像。

2.2 Marmousi模型试算

采用国际通用的Marmousi模型进行了试算, 图 8是Marmousi模型的速度场和Q值模型, 该模型含有较多的断层、背斜等, 在速度为2000~4000m/s范围内存在较强的衰减。模型纵、横向采样点数为737×380, 纵、横向网格间距为10m。震源为主频20Hz的雷克子波, 共184炮, 炮间距为40m, 炮点分布在网格点上。每炮737个接收点, 接收点间距10m, 采用粘声介质交错网格正演模拟得到炮记录, 时间长度为6.0s, 采样间隔为1.0ms。

对声波最小二乘逆时偏移方法和粘声波最小二乘逆时偏移方法进行了对比测试。图 9a为声波最小二乘逆时偏移第1次迭代结果, 图 9b为粘声波最小二乘逆时偏移第1次迭代结果, 可以看出, 无论是声波最小二乘逆时偏移, 还是粘声波最小二乘逆时偏移, 第1次迭代结果都存在较强的浅层噪声, 成像结果模糊不清, 分辨率较低, 但是在模型底部区域, 粘声波最小二乘逆时偏移的能量比声波最小二乘逆时偏移强, 同相轴更加明显, 说明该方法具有一定的能量补偿作用。图 9c与图 9d分别是声波最小二乘逆时偏移与粘声波最小二乘逆时偏移第50次迭代后的结果。可以看出, 随着迭代次数的增加, 两种方法都能有效地压制低频噪声, 成像结果分辨率更高, 同相轴更加清晰。但在高衰减值区域, 粘声波最小二乘逆时偏移结果的能量明显强于声波最小二乘偏移结果, 层位更加清晰, 同相轴更加明显, 成像位置准确, 构造刻画较好, 能量更加均衡, 保幅性更好, 说明该方法具有较好的振幅衰减补偿作用, 能够有效补偿因为地下介质吸收衰减而造成的能量损失, 避免成像位置错乱等问题。图 10对比了第301道数据偏移结果, 可以看出, 本文方法偏移的结果更加接近于真实反射系数, 中深层能量对比更加均衡。

3 结束语

本文基于GSLS波动方程, 在反演理论框架下, 引入粘声介质最小二乘理论, 推导了粘声介质最小二乘逆时偏移中的反偏移算子, 针对偏移中的不稳定现象, 构建了对应的粘声逆时偏移伴随算子, 建立了粘声介质最小二乘逆时偏移的目标函数, 给出了最小二乘逆时偏移的实现流程与步骤；通过层状模型和Marmousi模型试算, 验证了粘声介质最小二乘逆时偏移方法的正确性与有效性。研究结果表明：该方法可以合理消除地层吸收、透射及几何扩散等因素对地震振幅、频率及相位的影响, 提高成像分辨率；消除常规逆时偏移方法由于观测系统对地下一些地区照明弱导致的能量不均衡对成像效果的影响, 以及震源波场和接收点波场互相关引起的噪声, 改善复杂介质条件下地质体成像的振幅属性, 实现保幅、保真、高精度成像。

但是, 实际生产中地震资料质量较差, 噪声较多, 如何在低信噪比情况下取得良好的成像效果将是我们下一步研究的课题。而且实际资料中通常会存在倾角约束、井约束等先验信息, 如何利用先验信息降低对背景速度的依赖也是我们下一步努力的方向。