油气勘探的最终目标是储层描述。我国油气地震勘探的对象正在由大尺度的构造油气藏向小尺度的(薄互层)岩性油气藏转变, 因此需要采用宽带波阻抗成像结果进行岩性油气藏储层的精确描述^[1]。基于散射和逆散射的全波形反演(Full Waveform Inversion, FWI)成像技术试图利用地表观测的叠前数据估计地下宽波数带的参数场, 基于此反演结果, 直接进行储层描述。然而, 实际情况下, 地震波场与地下介质之间的关系十分复杂, 具有很强的非线性性, 仅仅利用地表观测的有限带宽、有限孔径和无假频数据, FWI的实用化并不成功。当前, 它仅仅能在低频和长偏移距地表观测数据下, 得到较好的背景速度的估计^[2]。另一方面, 勘探地震中地下介质以层状为主, 可以将宽带阻抗分为背景阻抗和反射系数两部分分别进行线性反演。随着油气地震勘探技术逐渐向单点高密度、宽带、宽方位地震勘探方向转变, “两宽一高”的地震数据为精确的背景(偏移)速度和宽带反射系数估计提供了数据基础。结合来自测井数据的密度建模, 采用合理的深度域反演方法, 将背景阻抗和宽带反射系数融合成宽带波阻抗成像结果, 是进行精确油气藏描述的重要技术方向^[1]。

假设已经得到准确的背景速度, 本文聚焦于采用最小二乘叠前深度偏移(LS-PSDM)估计地下保真、高分辨率的反射系数, 为宽带的波阻抗成像奠定基础。LS-PSDM有数据域和成像域两种求解手段, 其本质目标都是估计Hessian矩阵的逆, 作用在常规的成像剖面上, 从而得到真振幅、高分辨率的反演成像结果。原则上, 数据域与成像域LS-PSDM等价, 但是两者的具体实现过程不同, 也会导致不同的求解代价及反演结果。

数据域求解LS-PSDM基于反偏移数据与观测数据的匹配, 利用数据残差计算梯度进行迭代求解。每轮迭代都需要进行反偏移、偏移以及计算步长, 至少是两次常规偏移的计算量。通常需要多达10次的迭代才可以得到满意的结果, 这导致数据域求解LS-PSDM的计算量十分巨大。另外, 由于实际数据中含有噪声、背景速度不准、子波未知等因素会导致迭代收敛缓慢甚至不收敛, 无法得到满意的结果。实际情况下, 数据域LS-PSDM由于其高昂的计算成本以及较慢的收敛速率, 通常不能很好地应用于实际生产中。

相对于数据域的LS-PSDM, 成像域LS-PSDM可以视为图像去模糊问题, 直接对成像剖面进行处理, 避免了数据域迭代中的反偏移、偏移计算, 极大地减少了计算量。成像域LS-PSDM的核心是显式计算Hessian矩阵(逆)。由于全Hessian矩阵十分庞大, 无法直接计算、存储, 所以需要引入近似计算策略使得成像域LS-PSDM得以应用。

本文先对数据域与成像域的LS-PSDM方法原理进行对比, 再介绍一种基于非平稳滤波算子构造近似的Hessian矩阵逆的方法^[3]。采用该方法得到的Hessian矩阵逆的近似可以直接作用在常规成像结果上, 得到成像域LS-PSDM的反演结果, 也可以作为预条件算子, 加速数据域LS-PSDM的迭代收敛。数值试验结果证明, 利用求解非平稳滤波算子近似构造的Hessian矩阵的逆实现成像域的LS-PSDM, 可以有效地均衡常规成像结果的振幅, 并且在一定程度上提高了成像分辨率。

1 数据域与成像域的LS-PSDM

假设波传播算子f能够很好地模拟观测数据, 并且数据中的噪声满足高斯分布, 基于Bayes反演理论的全波形反演(FWI)通过求解如下目标函数估计地下模型参数:

$ \mathop {\min }\limits_\mathit{\boldsymbol{m}} J(\mathit{\boldsymbol{m}}) = \frac{1}{2}\left\| {{\mathit{\boldsymbol{d}}^{{\rm{obs}}}} - f(\mathit{\boldsymbol{m}})} \right\|_2^2 $

(1)

其中, d^obs代表地表观测数据, m为地下模型参数。实际情况下, 波场与地下介质之间的关系十分复杂, 具有很强的非线性性, 因此FWI是一个强非线性问题, 理论上需要采用非线性寻优算法(例如蒙特卡洛方法)来求解。但是, 对于庞大的叠前地震数据, 由于计算成本等因素的限制, 当前并没有可行的非线性寻优算法求解FWI问题。因此, 基于不同的初始模型, FWI通常被近似为一系列线性问题来求解。LS-PSDM为其中一个线性化子问题, 在已知地下模型参数光滑背景m₀的基础上, 估计反射系数。一般地, LS-PSDM利用线性化的正算子(反偏移算子)L建立反射系数Δm与一阶反射数据d^pri间的线性关系:

$ {\mathit{\boldsymbol{d}}^{{\rm{pri}}}} = \mathit{\boldsymbol{L}}\Delta \mathit{\boldsymbol{m}} $

(2)

此时, 可以反演线性目标函数(公式(3)), 估计地下反射系数:

$ \mathop {\min }\limits_{\Delta \mathit{\boldsymbol{m}}} J(\Delta \mathit{\boldsymbol{m}}) = \frac{1}{2}\left\| {{\mathit{\boldsymbol{d}}^{{\rm{pri}}}} - \mathit{\boldsymbol{L}}\Delta \mathit{\boldsymbol{m}}} \right\|_2^2 $

(3)

对于线性目标函数(公式(3)), 数学上采用梯度导引的方法求解。利用正算子L预测当前模型下的数据d^cal, 构造线性正算子的伴随算子L^T, 将观测数据与预测数据的残差反投影到模型空间, 便可以得到线性目标函数(公式(3))的梯度:

$ {\mathop{\rm grad}\nolimits} (\Delta \mathit{\boldsymbol{m}}) = {\mathit{\boldsymbol{L}}^{\rm{T}}}\left( {{\mathit{\boldsymbol{d}}^{{\rm{pri}}}} - {\mathit{\boldsymbol{d}}^{{\rm{cal}}}}} \right) $

(4)

其中, d^cal和L^T需要满足如下关系:

$ \left\langle {\mathit{\boldsymbol{L}}\Delta \mathit{\boldsymbol{m}}, {\mathit{\boldsymbol{d}}^{{\rm{cal}}}}} \right\rangle = \left\langle {{\mathit{\boldsymbol{L}}^{\rm{T}}}{\mathit{\boldsymbol{d}}^{{\rm{cal}}}}, \Delta \mathit{\boldsymbol{m}}} \right\rangle $

(5)

其中, “〈〉”代表内积。梯度是求解线性问题(公式(3))的核心, 不同的迭代算法如最速下降法、共轭梯度法等都是基于梯度构造不同的迭代策略来求解线性问题的方法。然而本质上, 梯度是由正算子L决定的, 只要给定具体的正算子, 梯度便可以用公式(4)计算得到。这说明, 对于线性化的地球物理反演问题, 线性正算子L的选择决定了线性反演的成败, 如果选定的线性正算子具有较好的表达数据的能力, 那么线性反演便具有较好的收敛性。反之, 线性反演会收敛缓慢甚至不收敛, 这对于求解诸如LS-PSDM的大规模地球物理线性反问题是不可接受的。因此, 求解LS-PSDM问题的第一步, 也是最重要的一步, 就是构造一个能够较好地表达一次反射波数据的线性正算子。

在确定了线性正算子的基础上, LS-PSDM剩下的问题便是数值求解线性目标函数(公式(3))。一般有数据域和成像域两种求解手段, 其目标都是估计线性算子L的Hessian矩阵的逆, 并将其作用于常规的成像剖面上, 从而得到真振幅、高分辨率的反演成像结果。

数据域的LS-PSDM基于梯度公式(4), 迭代地求解线性目标函数(公式(3))。当目标函数收敛到给定的阈值范围内(或者达到给定的迭代次数)时, 认为得到了数据域LS-PSDM的解。以最速下降法为例, 对于每一步迭代, 利用线性正算子计算当前模型下的反偏移数据(第一次迭代该步骤省略), 将预测数据与观测的一次反射波数据的残差通过正算子的伴随算子反投影(偏移), 得到本次迭代的梯度。由此可以看出, 数据域LS-PSDM的一次迭代, 至少需要一次全空间的线性正演(用于计算数据残差), 一次全空间的偏移(用于计算梯度), 再加上估计步长的计算量。一般地, 数据域LS-PSDM至少需要10次迭代才能收敛到令人满意的水平^[4], 而对于实际数据而言, 由于线性正算子不能很好地模拟实际波场, 导致数据域的迭代收敛通常很缓慢, 甚至不收敛^[5]。

成像域的LS-PSDM从另一个途径求解目标函数(公式(3))。对于线性目标函数, 当其梯度(公式(4))等于0时, 可得法方程:

$ {\mathit{\boldsymbol{L}}^{\rm{T}}}\mathit{\boldsymbol{L}}\Delta \mathit{\boldsymbol{m}} = {\mathit{\boldsymbol{L}}^{\rm{T}}}{\mathit{\boldsymbol{d}}^{{\rm{pri}}}} $

(6)

令H=L^TL, 代表线性算子L的Hessian矩阵, m₁=L^Td^pri, 代表常规的偏移成像结果。成像域LS-PSDM试图直接求解法方程(6), 估计Hessian矩阵的逆并将其作用到常规成像结果上, 从而得到LS-PSDM的解:

$ \Delta \mathit{\boldsymbol{\tilde m}} = {\left( {{\mathit{\boldsymbol{L}}^{\rm{T}}}\mathit{\boldsymbol{L}}} \right)^{ - 1}}{\mathit{\boldsymbol{L}}^{\rm{T}}}{\mathit{\boldsymbol{d}}^{{\rm{pri}}}} $

(7)

相较于数据域的LS-PSDM需要在全空间进行线性正演、偏移, 成像域LS-PSDM操作灵活, 可以仅对感兴趣的目标区域反演, 进而节约计算成本^[6]。但是, 成像域LS-PSDM需要显式构造Hessian矩阵, 并且对其求逆, 而Hessian矩阵是一个极大规模矩阵, 其规模大小为模型大小的平方, 以目前的计算能力, 无法对Hessian矩阵直接求解及存储, 求逆更是十分困难。所以, 寻求Hessian矩阵逆的快速近似构造方法具有非常重要的实际意义。为了克服这一问题, 使成像域LS-PSDM实用化, 我们提出了一种基于非平稳滤波算子近似Hessian矩阵逆的方法, 大幅降低了计算量, 并且可以得到稳定的成像域LS-PSDM近似解, 成像质量相较常规成像结果有了较大幅度的提升。

2 基于非平稳滤波算子的成像域LS-PSDM

给定观测的一阶反射数据d^pri, 可得常规偏移成像结果m₁:

$ {\mathit{\boldsymbol{m}}_1} = {\mathit{\boldsymbol{L}}^{\rm{T}}}{\mathit{\boldsymbol{d}}^{{\rm{pri}}}} $

(8)

利用线性正算子L, 基于常规偏移成像结果m₁进行反偏移, 获得反偏移数据d₁=Lm₁, 利用偏移算子L^T对反偏移数据d₁再做一次偏移, 获得成像结果m₂:

$ {\mathit{\boldsymbol{m}}_2} = {\mathit{\boldsymbol{L}}^{\rm{T}}}{\mathit{\boldsymbol{d}}_1} $

(9)

此时, m₁与m₂有如下关系:

$ {\mathit{\boldsymbol{m}}_2} = {\mathit{\boldsymbol{L}}^{\rm{T}}}\mathit{\boldsymbol{L}}{\mathit{\boldsymbol{m}}_1} = \mathit{\boldsymbol{H}}{\mathit{\boldsymbol{m}}_1} $

(10)

由(10)式可知, m₂为m₁经过Hessian矩阵模糊作用后的成像结果, 常规成像结果m₁为地下真实反射系数经过Hessian矩阵模糊作用后的结果。由于m₁与m₂已知, 可以构造非平稳匹配滤波算子F, 使得:

$ {\mathit{\boldsymbol{m}}_{\rm{1}}}{\rm{ = }}\mathit{\boldsymbol{F}}{\mathit{\boldsymbol{m}}_{\rm{2}}} $

(11)

其中, F≈(L^TL)^-1, 可认为是Hessian算子逆的低秩近似。

为了求取匹配滤波算子F, 可构造如下最小二乘问题:

$ J(\mathit{\boldsymbol{F}}) = \frac{1}{2}\left\| {\mathit{\boldsymbol{F}}{\mathit{\boldsymbol{m}}_2} - {\mathit{\boldsymbol{m}}_1}} \right\|_2^2 $

(12)

对于非平稳滤波算子, 原则上对于模型空间每一个点i, 存在一组滤波系数。F_{j, i}代表对应模型第i个点的第j个滤波系数。由于非平稳滤波算子F系数的数目大于模型空间的样点数, 所以(12)式为一个欠定问题, 需要额外的先验信息(正则化)对滤波器系数进行约束。由于滤波算子F为Hessian算子逆的近似, 应该包含模型的结构信息, 同时, 模型不同样点处的滤波系数应该光滑变化, 因此, 本文采用具有保边界能力的总变差(total variation, TV)正则化对非平稳滤波算子F的求取进行约束, 以便得到稳定的滤波算子求解结果。由于F≈(L^TL)^-1, 将求得的滤波算子F作用到常规成像结果m₁上, 便可得到成像域LS-PSDM的近似结果:

$ \mathit{\boldsymbol{\hat m}} = \mathit{\boldsymbol{F}}{\mathit{\boldsymbol{m}}_1} \approx {\left( {{\mathit{\boldsymbol{L}}^{\rm{T}}}\mathit{\boldsymbol{L}}} \right)^{ - 1}}{\mathit{\boldsymbol{m}}_1} $

(13)

综上所述, 基于非平稳滤波算子的成像域LS-PSDM方法可总结如下:

1) 给定观测数据d^pri, 进行一次常规偏移, 得到成像结果m₁;

2) 对成像结果m₁进行反偏移, 对反偏移数据再进行一次偏移, 得到成像结果m₂;

3) 利用m₁与m₂估计非平稳滤波算子F;

4) 将估计得到的滤波算子F作用到成像结果m₁上, 得到近似的成像域LS-PSDM反演结果。

上述方法仅需要进行两次偏移, 一次反偏移, 以及求取非平稳滤波算子的计算量。相对于迭代求解数据域LS-PSDM, 计算量大大降低, 并且避免了数据域迭代收敛缓慢甚至不收敛的问题。对于不同的传播、偏移算子, m₁与m₂已经包含了所有波传播的物理规律、算子特性, 所以求取非平稳滤波算子的方式一样。因此, 该方法具备较好的通用性, 适用于各类复杂的波传播、偏移算子(如弹性波算子、各向异性算子等等)。

3 TV正则化约束下非平稳滤波算子的构造

对于匹配滤波问题, 有如公式(11)所示的线性关系。由于滤波算子F为Hessian矩阵逆的近似, 应该包含模型的结构信息, 并且模型相邻样点处的滤波算子系数应该光滑变化, 因此利用TV正则化约束模型不同点对应滤波算子系数的变化。构造如下有约束的滤波算子估计问题:

$ \begin{array}{l} \mathop {\min }\limits_\mathit{\boldsymbol{F}} \sum\limits_j {\int {\int {\sqrt {{{\left( {{\boldsymbol\nabla _x}{\mathit{\boldsymbol{F}}_j}} \right)}^2} + {{\left( {{\boldsymbol\nabla _z}{\mathit{\boldsymbol{F}}_j}} \right)}^2}} } } } {\rm{d}}x{\rm{d}}z\quad {\rm{s}}.{\rm{t}}.\\ \;\;\;\;\;\;\;\left\| {\mathit{\boldsymbol{F}}{\mathit{\boldsymbol{m}}_2} - {\mathit{\boldsymbol{m}}_1}} \right\|\;_2^2 < \varepsilon \end{array} $

(14)

其中, ▽_xF_j=F_j[i+1, k]-F_j[i, k], ▽_zF_j=F_j[i, k+1]-F_j[i, k], 表示模型空间横向、纵向一阶向前差分算子, F_j[i, k]代表模型空间点(i, k)处的第j个滤波器系数, ε代表数据噪声水平。基于Bregman迭代算法^[7], (14)式可采用(15)式和(16)式迭代求解:

$ \begin{array}{l} {\mathit{\boldsymbol{F}}^{k + 1}} = \mathop {\min }\limits_\mathit{\boldsymbol{F}} \sum\limits_j {\int {\int {\sqrt {{{\left( {{\boldsymbol\nabla _x}{\mathit{\boldsymbol{F}}_j}} \right)}^2} + {{\left( {{\boldsymbol\nabla _z}{\mathit{\boldsymbol{F}}_j}} \right)}^2}} } } } {\rm{d}}x{\rm{d}}z + \\ \;\;\;\;\;\;\;\;\;\;\frac{\mu }{2}\left\| {\mathit{\boldsymbol{F}}{\mathit{\boldsymbol{m}}_2} - \mathit{\boldsymbol{m}}_1^k} \right\|_2^2 \end{array} $

(15)

$ \mathit{\boldsymbol{m}}_1^{k + 1} = \mathit{\boldsymbol{m}}_1^k + {\mathit{\boldsymbol{m}}_1} - {\mathit{\boldsymbol{F}}^{\rm{k}}}{\mathit{\boldsymbol{m}}_2} $

(16)

(15) 式为无约束凸优化问题, 可采用Split-Bregman算法^[8], 将基于L1范数的总变差约束项与基于L2范数的数据逼近项解耦。求解(15)式的步骤如下。

首先进行变量替换:

$ \begin{array}{l} \mathop {\min }\limits_{\mathit{\boldsymbol{F}}, {\mathit{\boldsymbol{d}}_x}, {\mathit{\boldsymbol{d}}_z}} \sum\limits_j {\int {\int {\sqrt {{{\left( {{\mathit{\boldsymbol{d}}_\mathit{x}}} \right)}^2} + {{\left( {{\mathit{\boldsymbol{d}}_z}} \right)}^2}} } } } {\rm{d}}x{\rm{d}}z + \frac{\mu }{2}\left\| {\mathit{\boldsymbol{F}}{\mathit{\boldsymbol{m}}_2} - } \right.\\ \;\;\;\left. {\mathit{\boldsymbol{m}}_1^k} \right\|_2^2\quad {\rm{ s}}{\rm{.t}}{\rm{. }}\quad {\mathit{\boldsymbol{d}}_x} = {\boldsymbol\nabla _x}{\mathit{\boldsymbol{F}}_j}, {\mathit{\boldsymbol{d}}_z} = {\boldsymbol\nabla _z}{\mathit{\boldsymbol{F}}_j} \end{array} $

(17)

(17) 式显然与(15)式等价。利用Bregman迭代算法求解(17)式, 有:

$ \begin{array}{l} \mathop {\min }\limits_{\mathit{\boldsymbol{F}}, {\mathit{\boldsymbol{d}}_x}, {\mathit{\boldsymbol{d}}_z}} \sum\limits_j {\int {\int {\sqrt {\mathit{\boldsymbol{d}}_x^2 + \mathit{\boldsymbol{d}}_z^2} } } } {\rm{d}}x{\rm{d}}y + \frac{\mu }{2}\left\| {\mathit{\boldsymbol{F}}{\mathit{\boldsymbol{m}}_2} - \mathit{\boldsymbol{m}}_1^k} \right\|_2^2 + \\ \frac{\lambda }{2}\left\| {{\mathit{\boldsymbol{d}}_x} - {\boldsymbol\nabla _x}{\mathit{\boldsymbol{F}}_j} - \mathit{\boldsymbol{b}}_x^k} \right\|_2^2 + \frac{\lambda }{2}\left\| {{\mathit{\boldsymbol{d}}_z} - {\boldsymbol\nabla _z}{\mathit{\boldsymbol{F}}_j} - \mathit{\boldsymbol{b}}_z^k} \right\|_2^2 \end{array} $

(18)

$ \mathit{\boldsymbol{b}}_x^{k + 1} = \mathit{\boldsymbol{b}}_x^k + \left( {{\boldsymbol\nabla _x}\mathit{\boldsymbol{F}}_j^{k + 1} - \mathit{\boldsymbol{d}}_x^{k + 1}} \right) $

(19)

$ \mathit{\boldsymbol{b}}_z^{k + 1} = \mathit{\boldsymbol{b}}_z^k + \left( {{\boldsymbol\nabla _z}\mathit{\boldsymbol{F}}_j^{k + 1} - \mathit{\boldsymbol{d}}_z^{k + 1}} \right) $

(20)

现在问题集中在求解(18)式, 可分解为如下两步:

$ \begin{array}{*{20}{l}} {\mathit{\boldsymbol{F}} = \mathop {\min }\limits_\mathit{\boldsymbol{F}} \sum\limits_j {\frac{\mu }{2}} \left\| {\mathit{\boldsymbol{F}}{\mathit{\boldsymbol{m}}_2} - \mathit{\boldsymbol{m}}_1^k} \right\|_2^2 + \frac{\lambda }{2}{\mathit{\boldsymbol{d}}_x} - }\\ {\;\;\;\;\;\;{\boldsymbol\nabla _x}{\mathit{\boldsymbol{F}}_j} - \mathit{\boldsymbol{b}}_x^k\left\| {_2^2 + \frac{\lambda }{2}} \right.\left\| {{\mathit{\boldsymbol{d}}_z} - {\boldsymbol\nabla _z}{\mathit{\boldsymbol{F}}_j} - \mathit{\boldsymbol{b}}_z^k} \right\|_2^2} \end{array} $

(21)

$ \begin{array}{*{20}{l}} {\left( {{\mathit{\boldsymbol{d}}_x}, {\mathit{\boldsymbol{d}}_z}} \right) = \mathop {\min }\limits_{{\mathit{\boldsymbol{d}}_x}, {\mathit{\boldsymbol{d}}_z}} \sum\limits_j {\int {\int {\sqrt {\mathit{\boldsymbol{d}}_x^2 + \mathit{\boldsymbol{d}}_z^2} } } } {\rm{d}}x{\rm{d}}y + \frac{\lambda }{2}\left\| {{\mathit{\boldsymbol{d}}_x} - } \right.}\\ {\;\;\;\;\;\;{\boldsymbol\nabla _x}{\mathit{\boldsymbol{F}}_j} - \mathit{\boldsymbol{b}}_x^k\left\| {_2^2 + \frac{\lambda }{2}} \right.\left\| {{\mathit{\boldsymbol{d}}_z} - {\boldsymbol\nabla _z}{\mathit{\boldsymbol{F}}_j} - \mathit{\boldsymbol{b}}_z^k} \right\|_2^2} \end{array} $

(22)

其中, (21)式是一个二次凸优化问题, 可以利用梯度导引法(共轭梯度法)求解。利用广义的shrinkage公式, (22)式有显式的求解算法:

$ \mathit{\boldsymbol{d}}_x^{k + 1} = \max \left( {{\mathit{\boldsymbol{s}}^k} - \frac{1}{\lambda }, 0} \right)\frac{{{\boldsymbol\nabla _x}\mathit{\boldsymbol{F}}_j^{k + 1} + \mathit{\boldsymbol{b}}_x^k}}{{{\mathit{\boldsymbol{s}}^k}}} $

(23)

$ \mathit{\boldsymbol{d}}_z^{k + 1} = \max \left( {{\mathit{\boldsymbol{s}}^k} - \frac{1}{\lambda }, 0} \right)\frac{{{\boldsymbol\nabla _z}\mathit{\boldsymbol{F}}_j^{k + 1} + \mathit{\boldsymbol{b}}_z^k}}{{{\mathit{\boldsymbol{s}}^k}}} $

(24)

其中,

$ {\mathit{\boldsymbol{s}}^k} = \sqrt {{{\left| {{\boldsymbol\nabla _x}\mathit{\boldsymbol{F}}_j^{k + 1} + \mathit{\boldsymbol{b}}_x^k} \right|}^2} + {{\left| {{\boldsymbol\nabla _z}\mathit{\boldsymbol{F}}_j^{k + 1} + \mathit{\boldsymbol{b}}_z^k} \right|}^2}} $

(25)

综上所述, TV正则化约束下, 非平稳匹配滤波器F的计算方法为:

Initialize:F₀=0, b_{i, j}=d_{i, j}=0

While‖Fm₂-m₁‖$_2^2$>ε(求解约束问题公式(14))

for i=1 to N(求解无约束问题公式(15))

$ \begin{array}{l} {\mathit{\boldsymbol{F}}^{k + 1}} = \mathop {\min }\limits_\mathit{\boldsymbol{F}} \sum\limits_j {\frac{\mu }{2}} \left\| {\mathit{\boldsymbol{F}}{\mathit{\boldsymbol{m}}_2} - \mathit{\boldsymbol{m}}_1^k} \right\|_2^2\left. { + \frac{\lambda }{2}} \right\|{\mathit{\boldsymbol{d}}_\mathit{x}} - \\ {\boldsymbol\nabla _x}{\mathit{\boldsymbol{F}}_j} - \mathit{\boldsymbol{b}}_x^k\left\| {_2^2 + \frac{\lambda }{2}} \right\|\left. {{\mathit{\boldsymbol{d}}_z} - {\boldsymbol\nabla _z}{\mathit{\boldsymbol{F}}_j} - \mathit{\boldsymbol{b}}_z^k} \right\|_2^2 \end{array} $

$ {\mathit{\boldsymbol{s}}^k} = \sqrt {{{\left| {{\boldsymbol\nabla _x}\mathit{\boldsymbol{F}}_j^{k + 1} + \mathit{\boldsymbol{b}}_x^k} \right|}^2} + {{\left| {{\boldsymbol\nabla _z}\mathit{\boldsymbol{F}}_j^{k + 1} + \mathit{\boldsymbol{b}}_z^k} \right|}^2}} $

$ \mathit{\boldsymbol{d}}_x^{k + 1} = \max \left( {{\mathit{\boldsymbol{s}}^k} - \frac{1}{\lambda }, 0} \right)\frac{{{\boldsymbol\nabla _x}\mathit{\boldsymbol{F}}_j^{k + 1} + \mathit{\boldsymbol{b}}_x^k}}{{{\mathit{\boldsymbol{s}}^k}}} $

$ \mathit{\boldsymbol{d}}_z^{k + 1} = \max \left( {{\mathit{\boldsymbol{s}}^k} - \frac{1}{\lambda }, 0} \right)\frac{{{\boldsymbol\nabla _z}\mathit{\boldsymbol{F}}_j^{k + 1} + \mathit{\boldsymbol{b}}_z^k}}{{{\mathit{\boldsymbol{s}}^k}}} $

$ \mathit{\boldsymbol{b}}_x^{k + 1} = \mathit{\boldsymbol{b}}_x^k + \left( {{\boldsymbol\nabla _x}\mathit{\boldsymbol{F}}_j^{k + 1} - \mathit{\boldsymbol{d}}_x^{k + 1}} \right) $

$ \mathit{\boldsymbol{b}}_z^{k + 1} = \mathit{\boldsymbol{b}}_z^k + \left( {{\boldsymbol\nabla _z}\mathit{\boldsymbol{F}}_j^{k + 1} - \mathit{\boldsymbol{d}}_z^{k + 1}} \right) $

end

$ \mathit{\boldsymbol{m}}_1^{k + 1} = \mathit{\boldsymbol{m}}_1^k + {\mathit{\boldsymbol{m}}_1} - {\mathit{\boldsymbol{F}}^k}{\mathit{\boldsymbol{m}}_2} $

end

采用以上求解算法, 可以稳健、高效地计算出非平稳滤波器F, 即Hessian算子逆的近似。

4 数值试验

图 1为局部Sigsbee模型。根据该模型背景(光滑)成分和观测数据, 首先进行常规偏移, 得到偏移成像结果m₁, 如图 2所示。基于该成像结果进行反偏移, 对反偏移数据再进行一次偏移成像, 得到第2个成像结果m₂, 如图 3所示。

由图 2和图 3可以看出, 两个成像结果有较大差别, m₂在深层及边界处照明不足的区域幅值更弱, 分辨率更低。根据前文介绍的方法利用m₁和m₂求解非平稳滤波算子F, 可视为Hessian矩阵逆的低秩近似。对于模型每一点, 选取的滤波器相当于一个11×11的二维褶积算子, 计算得到的非平稳滤波器零延迟系数如图 4所示。从图 4可以看出, 计算得到的非平稳滤波算子零延迟系数对模型深层以及两边照明较弱的区域有较大的补偿, 而且携带了模型的构造信息, 这也正体现了Hessian矩阵逆对角线的作用。将求得的非平稳滤波器F作用到常规成像结果m₁上, 可以得到成像域LS-PSDM反演结果, 如图 5所示。对比图 5与图 2可以看出, 成像域LS-PSDM反演结果比常规偏移成像结果振幅更均衡, 弱照明以及深层成像区域的振幅得到了充分的提高, 使得成像结果更保真, 分辨率也得到了提高。

对地表 500m, 3 750m以及7 000m处成像域LS-PSDM反演结果、常规偏移成像结果和真实反射系数进行抽道(图 6到图 8)。由图 6到图 8可以明显地看出, 与常规偏移成像结果相比, 成像域LS-PSDM反演可以得到更保真的反演成像结果, 分辨率也有所提高, 尤其是在弱照明和深层区域。

采用如图 9所示的Marmousi速度模型进一步说明本文方法对成像结果的改善效果。Marmousi速度模型的常规偏移结果如图 10所示。对常规偏移结果再进行反偏移、偏移后的成像结果如图 11所示。利用m₁, m₂构造非平稳滤波器, 对于模型每一点, 选取一个11×11的二维褶积算子, 计算得到的非平稳滤波算子零延迟系数如图 12所示。将求得的非平稳滤波算子F作用到常规偏移成像结果m₁上, 得到成像域LS-PSDM反演结果, 如图 13所示。对比图 13与图 10可以看出, 通过计算近似的Hessian算子逆, 作用到常规偏移成像结果上, 得到的成像域LS-PSDM结果比常规偏移成像结果振幅更均衡, 弱照明区域的振幅得到了充分的提高, 成像结果更保真。

5 讨论

油气地震勘探的目标是识别与描述含油气储层。宽带波阻抗成像更有利于岩性油气藏的(定量)精细描述。“两宽一高”地震数据采集已经成为地震数据采集的核心方法技术, 为宽带波阻抗的估计提供了数据基础。尽管理论上FWI具有直接采用叠前地震数据估计宽带波阻抗的潜力, 但是由于地下波场与介质之间的强非线性关系, 对于实际问题, 很难采用求解非线性的FWI问题得到理想的宽带阻抗成像结果。由于地震勘探中地下介质以层状介质为主, 可以将宽带阻抗分为背景阻抗和反射系数两部分分别进行线性反演, 进而将二者融合成宽带波阻抗成像结果, 这是进行精确油气藏描述的重要技术方向。假设通过层析反演已经得到较为准确的背景速度, 利用测井等信息进行背景密度建模, 可以得到背景阻抗。高保真、高分辨率的反射系数由最小二乘偏移估计得到。利用合理的深度域反演方法, 可以将背景阻抗以及高保真、高分辨率的反射系数的波数成分融合, 并转换为宽带(绝对)波阻抗^[9]。

本文讨论了高保真、高分辨率反射系数估计方法。相对于只能正确地定位反射系数出现空间位置的常规偏移成像, 最小二乘偏移试图将Hessian算子的逆作用在常规偏移成像结果上, 具有补偿由于观测系统不完备以及上覆介质复杂造成的地下照明不均匀, 校正成像振幅, 提高成像结果纵、横向分辨率, 消除成像假象等能力。因此, 最小二乘偏移成像将是今后主流的成像技术。

但是, 实际应用中, 常规数据域迭代的最小二乘偏移成像计算量非常大, 且迭代收敛速度缓慢, 甚至不收敛。究其原因, 可归结为:给定的背景速度不准, 地震子波未知, 线性化的正问题不能很好地预测数据等。鉴于上述原因, LS-PSDM是一个病态性极强的反问题, 正问题的不确定性以及解的非唯一性使得LS-PSDM目标函数的凸性很差, 迭代求解收敛缓慢甚至不收敛, 无法获得理想的成像结果。为了克服数据域迭代求解LS-PSDM计算量过大以及不易收敛的问题, 使得LS-PSDM更实用化, 基于图像去模糊(反褶积)的成像域LS-PSDM试图从另一个途径实现LS-PSDM的功能。成像域LS-PSDM的关键是估计线性算子的Hessian矩阵或者Hessian矩阵的逆, 进而对成像结果进行去模糊或者直接将Hessian矩阵的逆作用在成像剖面上, 提高成像质量。对于线性成像问题, 数据域迭代求解LS-PSDM与成像域的图像去模糊原理上等价, 但相对于数据域迭代的LS-PSDM, 成像域LS-PSDM操作更灵活, 可以仅对感兴趣的目标区域进行反演成像, 而且一旦估计出Hessian算子的逆, 便不需要进行迭代, 也就回避了数据域迭代收敛缓慢和不收敛的问题, 节约大量的计算时间, 更实用。但是, Hessian矩阵是一个超大规模矩阵, 对于实际三维工区, 以目前计算机能力, 直接对其计算和存储基本不可能。因此, 研究成像域LS-PSDM Hessian矩阵(逆)的近似计算具有很现实的意义。

6 结论

本文研究了成像域LS-PSDM Hessian矩阵逆的近似估计方法。利用常规成像剖面以及反偏移、偏移剖面构造非平稳滤波算子, 以近似Hessian矩阵的逆。由于非平稳滤波算子的求解具有较强的欠定性, 考虑到滤波器系数应具备成像结果的构造特征, 本文在求解滤波算子过程中引入具有保留图像边界能力的TV正则化约束, 以得到稳健、合理的滤波器求解结果。求得的非平稳滤波算子作为Hessian矩阵逆的近似直接作用到常规成像剖面上, 即可得到成像域LS-PSDM反演结果。该方法极大地降低了LS-PSDM的计算量, 并且计算稳定, 适用于各类偏移算子。基于局部Sigsbee模型以及Marmousi模型的数值试验结果表明, 求得的非平稳滤波算子能够较好地实现Hessian矩阵逆的功能, 得到的成像域LS-PSDM结果相对于常规偏移成像结果振幅更均衡, 弱照明区域的振幅得到了充分的提高, 使得成像结果更保真, 并且分辨率也有了一定提升。另外, 非平稳滤波算子也可以作为预条件算子, 应用于数据域LS-PSDM的迭代求解过程中, 加速迭代收敛, 以获得更高分辨率、振幅更加保真的反射系数估计结果。求解LS-PSDM得到的高保真、高分辨率反射系数, 为宽带波阻抗成像提供了重要的数据基础。

致谢: 感谢中石油勘探开发研究院及西北分院、中海油研究院和湛江分公司、中国石油化工股份有限公司石油物探技术研究院和胜利油田分公司对波现象与智能反演成像研究组(WPI)研究工作的资助与支持。