地震成像的目标是估计出地下介质的几何形态和弹性参数值, 为精确描述油气藏提供可靠的资料。当前生产中常用的偏移成像技术, 包括Kirchhoff射线类方法^[1-4]、单程波波动方程类方法^[5-8]和双程波波动方程类方法^[9-11], 其偏移算子都是波场正演传播算子的共轭转置, 而非逆算子。由于地震数据观测孔径有限、地下模型复杂以及地震波带宽有限, 采用常规偏移技术仅仅能够得到比较模糊的构造信息, 无法对复杂油气藏进行精细成像, 并且成像振幅不均衡, 影响后续岩性参数反演的可靠性和精度。最小二乘偏移是一种线性化的地震反演成像技术, 在最小二乘反演框架下估计一个最优的偏移成像结果, 相当于对常规偏移结果作用一个Hessian逆矩阵算子, 可以降低或者去除采集照明不均或者采集脚印所导致的偏移假象, 校正偏移结果中存在的振幅误差, 提高地震成像的分辨率。因此, 近年来最小二乘偏移成为地震成像领域的研究热点, 相继发展了最小二乘Kirchhoff偏移^[12]、最小二乘单程波偏移^[13]和最小二乘逆时偏移^[14-16]等技术。但是, 最小二乘偏移基于数据拟合方式逐次迭代寻求最优成像值, 波场传播算子不能准确预测出实际数据, 迭代反演计算量巨大, 收敛比较困难, 三维最小二乘逆时偏移尤为严重。国内外学者为此开展了大量的研究工作, 如TANG^[17]采用随机相位编码方式求取Hessian矩阵, 并提出了面向目标的最小二乘偏移方法; GUITTON^[18]利用初始偏移结果和反偏移-再偏移结果估计非稳相滤波器来近似Hessian逆矩阵, 提高了最小二乘偏移反演的效果; 王彦飞^[19]讨论了最小二乘偏移的预条件与正则化方法, 提高了反演的稳定性; ZENG等^[20]针对预测数据与观测数据的不一致性, 引入了一个可信度参数作为反演的预条件; ZHANG等^[21]、DUTTA^[22]、DUAN等^[23]将互相关误差泛函引入到最小二乘偏移成像框架中, 降低了最小二乘偏移方法对震源子波和振幅因素的敏感性。为了提高最小二乘偏移的计算效率, DAI等^[24]提出了平面波编码的最小二乘逆时偏移方法; 黄金强等^[25]将平面波最小二乘偏移思路应用于各向异性介质; 田文彬等^[26]采用导向滤波算子压制平面波最小二乘偏移假象, 同时提高了收敛效率和成像精度; 柯璇等^[27]则将GPU加速器应用于二维最小二乘逆时偏移算法中, 加速效果显著。

本文针对三维最小二乘逆时偏移的数据匹配和计算效率问题, 采用基于双程波傅里叶有限差分的反偏移方法、共轭梯度反演步长三点抛物搜索寻优方法和GPU加速偏移/反偏移核函数计算方法, 形成了一套具有实际处理能力的三维数据域最小二乘逆时偏移技术方案, 利用实际资料验证了方法的有效性。

1 方法原理

最小二乘偏移成像是指由已知的地震观测数据、假设的波传播方程和背景速度模型对地下介质的高波数扰动进行估计。由于观测数据信息不完整、波传播方程无法精确描述真实波场传播、震源子波未知以及观测数据存在噪声等因素影响, 地震反演无法获得精确的参数模型。现代地震反演理论采用贝叶斯估计思想对地震反演过程及反演解进行评价。贝叶斯公式为:

$ P(\mathit{\boldsymbol{m}}/\mathit{\boldsymbol{d}}) = \frac{{P(\mathit{\boldsymbol{d}}/\mathit{\boldsymbol{m}})P(\mathit{\boldsymbol{m}})}}{{P(\mathit{\boldsymbol{d}})}} $

(1)

式中:m为地下参数, d为地震数据, 两者关系由波场传播算子L建立, d=Lm; P(m/d)是后验概率密度函数, P(d)是数据空间的概率密度函数, P(m)是模型空间的先验概率密度函数, P(d/m)是条件概率密度函数。假设P(m)和P(d/m)符合高斯分布, P(d)符合标准的正态分布(观测系统是规则的), 那么后验概率密度可表示为:

$ P(\mathit{\boldsymbol{m}}/\mathit{\boldsymbol{d}}) \propto \exp \left\{ { - \frac{1}{2}{{\left( {L\mathit{\boldsymbol{m}} - {\mathit{\boldsymbol{d}}_{{\rm{obs}}}}} \right)}^{\rm{T}}}\mathit{\boldsymbol{C}}_{\rm{D}}^{ - 1}\left( {L\mathit{\boldsymbol{m}} - {\mathit{\boldsymbol{d}}_{{\rm{obs}}}}} \right) \cdot } \right.\left. \\ {{{\left( {\mathit{\boldsymbol{m}} - {\mathit{\boldsymbol{m}}_{{\rm{prior}}}}} \right)}^{\rm{T}}}\mathit{\boldsymbol{C}}_{\rm{M}}^{ - 1}\left( {\mathit{\boldsymbol{m}} - {\mathit{\boldsymbol{m}}_{{\rm{prior}}}}} \right)} \right\}{\rm{ }} $

(2)

式中, C_M表示模型空间的协方差矩阵, C_D表示数据空间的协方差矩阵, m_prior表示先验模型信息, d_obs表示观测数据。反演问题的求解可以认为是在P(m)、P(d)、P(d/m)的约束下, 使P(m/d)最大化。记代价函数为:

$ E(\mathit{\boldsymbol{m}}){\rm{ }} = \frac{1}{2}{\left( {L\mathit{\boldsymbol{m}} - {\mathit{\boldsymbol{d}}_{{\rm{obs}}}}} \right)^{\rm{T}}}\mathit{\boldsymbol{C}}_{\rm{D}}^{ - 1}\left( {L\mathit{\boldsymbol{m}} - {\mathit{\boldsymbol{d}}_{{\rm{obs}}}}} \right) \cdot {\left( {\mathit{\boldsymbol{m}} - {\mathit{\boldsymbol{m}}_{{\rm{prior}}}}} \right)^{\rm{T}}}\mathit{\boldsymbol{C}}_{\rm{M}}^{ - 1}\left( {\mathit{\boldsymbol{m}} - {\mathit{\boldsymbol{m}}_{{\rm{prior}}}}} \right){\rm{ }} $

(3)

那么, 可将后验概率密度最大化问题转化为代价函数最小化问题, 于是反演问题变成了一个非线性的、无约束的最优化问题。假设模型空间、数据空间的元素相互独立, 每个元素的方差相等, 且不考虑先验模型信息, 则得到经典的最小二乘偏移目标函数:

$ E(\mathit{\boldsymbol{m}}){\rm{ }} = \frac{1}{2}{\left( {L\mathit{\boldsymbol{m}} - {\mathit{\boldsymbol{d}}_{{\rm{obs}}}}} \right)^{\rm{T}}}\left( {L\mathit{\boldsymbol{m}} - {\mathit{\boldsymbol{d}}_{{\rm{obs}}}}} \right) $

(4)

其物理意义是使得观测数据与预测数据的L₂范数误差最小。

最小二乘反演问题通常采用数据域迭代算法(如最速下降法、共轭梯度法等)求解, 在一定假设条件下由初始值逐步逼近真实解。共轭梯度法最小二乘逆时偏移成像迭代公式为:

$ {\mathit{\boldsymbol{m}}_{k + 1}} = {\mathit{\boldsymbol{m}}_k} - {\mu _{k + 1}}{\mathit{\boldsymbol{d}}_k} $

(5)

式中, m_k、m_k+1分别是第k、k+1次迭代后的成像值。步长μ_k+1、共轭梯度方向d_k可分别表示为:

$ {{\mu _{k + 1}} = \frac{{{{\left( {{\mathit{\boldsymbol{d}}_k}} \right)}^{\rm{T}}}{\mathit{\boldsymbol{g}}_k}}}{{{{\left( {L{\mathit{\boldsymbol{d}}_k}} \right)}^{\rm{T}}}L{\mathit{\boldsymbol{d}}_k}}}} $

(6)

$ {{\mathit{\boldsymbol{d}}_k} = \left\{ {\begin{array}{*{20}{l}} { - {\mathit{\boldsymbol{g}}_k}}&{k = 1}\\ { - {\mathit{\boldsymbol{g}}_k} + {\alpha _k}{\mathit{\boldsymbol{d}}_{k - 1}}}&{k \ge 2} \end{array}} \right.} $

(7)

此处, g_k是泛函梯度, 由数据残差偏移成像得到。α_k是共轭梯度方向修正因子, 可采用下式构造:

$ {\alpha _k} = \frac{{{{({\mathit{\boldsymbol{g}}_k} - {\mathit{\boldsymbol{g}}_{k - 1}})}^{\rm{T}}}{\mathit{\boldsymbol{g}}_k}}}{{\mathit{\boldsymbol{g}}_{k - 1}^{\rm{T}}{\mathit{\boldsymbol{g}}_{k - 1}}}} $

(8)

数据域共轭梯度法最小二乘逆时偏移具体实现流程如图 1所示。

2 傅里叶有限差分反偏移方法

基于双程波方程的逆时偏移和反偏移是最小二乘逆时偏移算法的核心, 通常采用高阶差分波场传播算子实现。然而, 由于实际资料成像计算网格比较大, 采用有限差分方法存在较为严重的数值频散, 计算得到的反偏移数据与实际地震采集数据存在比较明显的频率差异。最小二乘偏移期望正演模拟算子能够完全模拟观测数据, 避免由于算子精度不够导致的数据不一致问题。如果无法避免这种问题出现, 应该尽量减少频散效应。为了降低反偏移计算中的数值频散, 我们可以通过提高有限差分的阶数或者减小计算网格的大小, 也可以采用更高精度的波场传播算法, 如伪谱法^[28-33]、拟解析法^[34]和傅里叶有限差分法^[35]等来实现。本文利用傅里叶有限差分法, 在时间-空间域与时间-波数域交替进行波场传播, 从而实现双程波反偏移正演模拟。该方法的主要思想是利用速度分解, 求解如下波场传播方程^[36]:

$ \begin{array}{*{20}{l}} {u(\mathit{\boldsymbol{x}}, t + \Delta t) + u(\mathit{\boldsymbol{x}}, t - \Delta t) - 2u(\mathit{\boldsymbol{x}}, t) = }\\ {2\int_{ - \infty }^{ + \infty } {\widehat u} (\mathit{\boldsymbol{k}}, t)(\cos (|\mathit{\boldsymbol{k}}|v\Delta t) - 1){{\rm{e}}^{ - {\rm{i}}\mathit{\boldsymbol{kx}}}}{\rm{d}}\mathit{\boldsymbol{k}}} \end{array} $

(9)

式中, x是空间位置; k是波数矢量; u(x, t)、u(x, t-Δt)、u(x, t+Δt)表示3个相邻时刻的时间-空间域波场; ${\widehat u}$ (k, t)表示时间-波数域波场。速度可分解为平均速度和扰动速度, 基于平均速度利用傅里叶正、反变换进行波场传播, 再通过有限差分算法校正速度扰动的影响, 实现双向波场在变速介质中的传播。

图 2是某工区实际资料反偏移记录的局部放大结果, 3个方向的计算网格大小分别是25、25、10m。由图 2可见, 采用常规高阶有限差分法反偏移的结果存在较为严重的数值频散, 在记录上表现为频散噪声严重(图 2a), 而基于傅里叶有限差分法的反偏移结果则较好地压制了频散噪声(图 2b)。

3 共轭梯度反演步长寻优方法

三维数据域最小二乘逆时偏移运算成本非常高, 除了因为双程波方程波场传播算子计算量大以外, 另一个更重要的原因是迭代收敛速度慢。共轭梯度法相对于最速下降法而言, 是一种收敛速度更快的反演迭代算法。然而, 经典的共轭梯度法反演步长计算公式是基于目标泛函梯度为0推导得到的, 与实际情况不符, 迭代更新效果往往不是很好。本文采用三点抛物曲线拟合步长寻优法, 假设误差泛函在解的邻域为抛物线型, 利用3个试探步长所对应的3个点拟合一条抛物线, 将抛物线极值点对应的步长作为迭代更新的最优步长。

假设3个步长α₀、α₁、α₂位于一个单峰区间内, 对应的目标函数值依次是E_α0、E_α1、E_α2, 于是可以构造一个抛物函数。计算其抛物线最小值对应的步长:

$ \alpha = \frac{1}{2}\frac{{\left( {\alpha _1^2 - \alpha _2^2} \right){E_{a0}} + \left( {\alpha _2^2 - \alpha _0^2} \right){E_{a1}} + \left( {\alpha _0^2 - \alpha _1^2} \right){E_{a2}}}}{{\left( {{\alpha _1} - {\alpha _2}} \right){E_{a0}} + \left( {{\alpha _2} - {\alpha _0}} \right){E_{a1}} + \left( {{\alpha _0} - {\alpha _1}} \right){E_{a2}}}} $

(10)

构造抛物曲线时, 3个试探步长点选择方法如下:α₀=0, E_α0为上一次成像更新后的误差, α₁是经典共轭梯度法公式计算出的步长, α₂=2α₁。如果不满足E_α2>E_α1, 则α₂=4α₁, 逐步增大α₂, 直到E_α2>E_α1为止。

从上述过程可以看出, 步长计算中最少需要两次试探步长后的反偏移模拟来拟合抛物线, 如果所有炮都进行反偏移, 则计算量巨大, 对于三维最小二乘偏移而言不可接受。为了减少计算量, 可以随机选择部分炮进行正演模拟和曲线拟合。图 3对比了二维模型数据最小二乘偏移迭代中全部炮的误差曲线和部分炮的误差曲线, 可见两者总体趋势一致, 说明通过部分炮的误差来拟合抛物曲线是可行的。

4 GPU加速计算方法

在最小二乘逆时偏移中, 无论是高阶有限差分逆时偏移算子还是傅里叶有限差分反偏移算子, 都存在计算量巨大的问题。GPU加速是应对这一问题的有效途径。

高阶有限差分算法的GPU并行方案比较简单, 在Y方向做循环, 在X、Z方向做多线程并行计算。逆时偏移中采取震源波场存储方案^[37], 仅存储边界波场, 在记录波场反传的同时利用边界数据重构震源波场, 这大大降低了震源波场的直接存储量。然而, 由于GPU内存较小, 在每一步波场延拓中需要将边界波场数据传送到CPU进行存储。

对于傅里叶有限差分算法的GPU加速, 由于每个时刻波场延拓包含正、反两次三维傅里叶变换, 计算效率很低, 因此采用GPU计算平台提供的三维傅里叶变换——cuFFT函数实现FFT, 加速20倍以上。图 4是基于GPU的三维傅里叶有限差分法反偏移流程, 整体反偏移计算效率提高12倍以上。

5 实际资料应用 5.1 西部探区碎屑岩成像

将本文三维数据域最小二乘逆时偏移技术应用于西部某工区实际地震资料处理。该区中浅层存在一套碎屑岩储层, 对成像分辨率要求较高。实际数据共4.8×10⁴炮, 计算网格为30m×30m×10m, 成像深度为1×10⁴m, Inline、Xline偏移孔径分别为6000m、5000m, 共完成了3次最小二乘逆时偏移成像迭代。图 5对比了逆时偏移和最小二乘逆时偏移的效果, 图 6对比了逆时偏移和最小二乘逆时偏移成像剖面的频谱。

从图 5、图 6可以看出, 最小二乘逆时偏移的分辨率更高, 储层内部刻画更清楚。图 7是基于逆时偏移和最小二乘逆时偏移成像剖面的沿层相干属性分析结果, 可以看出最小二乘逆时偏移结果中雁形断裂特征更加清晰。

5.2 南方山地地震成像

南方某山地地震资料整体信噪比偏低, 常规成像效果不佳。选择面积为60km²的数据进行最小二乘偏移成像试处理, 数据量为7190炮, 计算网格为25m×25m×10m, 成像深度为1.2×10⁴m, Inline、Xline偏移孔径分别为5000m、4000m, 共完成5次迭代。图 8对比了逆时偏移和最小二乘逆时偏移效果, 可见经过最小二乘逆时偏移迭代反演, 成像剖面信噪比有一定提升, 同相轴连续性变好, 断层成像更加清晰。

6 结束语

本文从反演方法理论、反偏移方法、共轭梯度步长寻优算法、GPU加速方法等几个方面阐述了一套具有实际处理能力的三维数据域最小二乘逆时偏移技术实施方案, 并通过两个探区实际地震资料处理验证了方法的效果和优势, 获得以下结论和认识:

1) 基于双程波傅里叶有限差分的反偏移方法能够有效压制数值频散, 提高反偏移模拟数据与实际观测数据的匹配程度, 有利于改善最小二乘逆时偏移成像效果。

2) 导致数据域最小二乘偏移效率低下的一个重要原因是反演迭代次数多, 必须提高共轭梯度步长寻优的效果和效率, 部分随机炮误差曲线拟合和三点抛物搜索是达到这一目的的有效技术手段。

3) 双程波波场传播算子计算量大是制约三维数据域最小二乘逆时偏移效率的另一个重要因素, 借助GPU并行计算能够显著提升三维数据域最小二乘逆时偏移的计算效率, 提升该技术实际应用的可行性。

4) 最小二乘逆时偏移能够提高地震成像分辨率, 改善同相轴的连续性, 整体成像效果优于逆时偏移。

由于三维数据域最小二乘逆时偏移计算效率低、实际数据匹配困难, 目前其规模化推广应用仍然较为困难。下一步的研究方向是先验信息的利用、Hessian矩阵的近似等, 以逐步打造一个规模化、实用化的最小二乘偏移反演成像技术。