将地震反偏移算子记为∇F, 观测地震数据记为d, 那么, 逆时偏移可以表示为∇F′d, 其中“′”表示复转置^[1]。相比于单次逆时偏移, 最小二乘逆时偏移所解决的成像问题可以表示为∇F^†d, 其中“†”表示伪逆运算^[2]。忽略高阶传播项, 地震数据d可以近似线性表示为∇Fδm, 其中δm表示地下的真实介质对比背景介质的模型扰动, 其直观表现为地震像。这样, 单次逆时偏移得到的地震像可以表示为∇F′∇Fδm。可以看出, 相比于真实的地震像, 单次逆时偏移得到的地震像中包含了波传播过程的影响, 体现在∇F′∇F这个算子中, 这个算子在Gauss-Newton方法中常被称为近似Hessian矩阵^[3]。相比于单次偏移, 最小二乘偏移得到的地震像中包含的传播算子, 理论上可以表示为∇F^†∇F≈I, 因此, 精确的最小二乘反演理论上可以给出真振幅真相位的地震像。然而, 相比于单次偏移, 最小二乘偏移涉及求逆运算, 实际运算中往往使用迭代算法求解, 因此, 其计算量取决于迭代次数, 往往是单次偏移计算量的几倍甚至几十倍。为了使最小二乘偏移实用化, 地球物理学家们提出利用炮编码的方法来降低其计算量^[4-5]。从代数理论角度理解，对震源和观测数据同时进行编码，提高了反演成像系统的多解性；从信号的角度分析，编码引入了不同炮集数据之间的交叉噪声, 从而影响反演结果的准确度和信噪比。

压缩感知理论的出现^[6-7]为解决规模庞大的地震反演问题提供了新思路。在地震反演偏移成像问题中, 对于一个压缩的物理观测模型, 压缩感知理论提出了:①精确求解原始观测模型所对应解的充分条件, 从而指导我们最优化设计压缩观测系统; ②如何利用稀疏优化的方法从压缩观测的数据中恢复真实的地震像。本文将介绍作者与HERRMANN等在这个领域共同研究取得的一些重要成果^[8-11]:①最小二乘逆时偏移的快速计算; ②消除逆时偏移中未知的子波参数对反演结果的影响; ③包含表面多次波的地震数据的成像。同时展示了合成和实际数据的实验结果, 以验证本文方法的可行性和有效性。

1 方法和合成数据实验 1.1 快速最小二乘偏移

HERRMANN等^[8]提出了基于压缩感知的最小二乘偏移方法, TU等^[9]在此基础上, 分析了波动方程线性化过程中引入的线性误差对该方法的影响。与传统的炮编码一样, 基于压缩感知的最小二乘偏移也是通过对数据降维来实现波动方程模拟量的减少。与传统的最小二乘偏移不同, 求解压缩感知问题的Basis Pursuit De-Noise(BPDN)模型^[6]引入了稀疏约束, 其在频率域求解的形式如下:

$ \begin{array}{l} \mathop {{\rm{argmin}}}\limits_x {\left\| \mathit{\boldsymbol{x}} \right\|_1}\\ {\rm{subject}}\;{\rm{to}}\sum\limits_{i \in \mathit{\Omega }} {\sum\limits_{j \in \mathit{\boldsymbol{ \boldsymbol{\varSigma} }}} {\left\| {{\mathit{\boldsymbol{d}}_{i,j}} - \nabla \mathit{\boldsymbol{F}}\left[ {{\mathit{\boldsymbol{m}}_0},{w_i}{\mathit{\boldsymbol{s}}_j}} \right] \cdot } \right.} } \\ \quad \quad \left. {{\mathit{\boldsymbol{C}}^*}\mathit{\boldsymbol{x}}} \right\|_2^2 \le {\sigma ^2} \end{array} $

(1)

式中:Ω表示对数据降维时选择的频率子集; Σ表示选择的编码炮集; i和j分别表示频率和炮集的编号; d表示随机降采样的观测数据(在本文中, 根据上下文可判断出是仅包含一次波的观测数据, 还是同时包含一次波和多次波的观测数据); m₀表示光滑的背景模型; w_i表示地震子波w在某个频率上的频谱值; s表示不包含子波效应的随机降采样后的震源波场(即主要包含与观测数据相对应的震源的位置信息, 子波效应包含在w中, 降采样方式和观测数据d一致); C表示曲波变换(C^*表示曲波反变换); x是模型扰动δm在曲波域的系数(δm=C^*x); σ是可调节参数, 用以控制匹配数据的精度, 其选择往往反映数据的信噪比。(1)式的意义是, 在所有可能的曲波系数中, 寻找一个最稀疏的系数x, 经过曲波反变换后得到模型扰动的一个估计, 再经过反偏移, 可以在由σ指定的误差范围内匹配观测数据d。引入曲波变换的意义在于, 由于地质模型在空间域并非稀疏, 曲波变换可将其稀疏化^[12-13]。直观上解释, 上述问题的含义是:在所有可能的曲波系数中, 寻找最稀疏(一范数最小)的解, 这个解经过反曲波变换(变换到空间域中), 得到地下介质的扰动模型, 该模型经过反偏移算子的作用, 能在σ指定的误差范围内匹配观测到的地震数据。HERRMANN等^[8]利用谱梯度投影算法(SPGL1)求解上述问题^[14]。

对比传统的最小二乘偏移算法, 压缩感知引入了稀疏约束(公式(1)第2行)。然而压缩感知对最小二乘偏移的影响还不止于此:压缩感知理论对数据如何进行降采样(同时对应对炮集如何进行编码)也有指导作用。压缩感知理论指出, 在采样过程中引入随机因素, 是最终恢复原始信号的关键条件。信号恢复的充分条件由Restricted Isometry Property(RIP)条件给出^[6]; 在地震信号的处理中, HENNENFENT等^[15]对随机采样的重要性亦作了直观的阐述。因此, 在对观测数据进行降采样时, 应随机选取频率的子集; 在编码炮集时, 也应引入足够的随机因素。详细的降采样方案请参考文献[8]和文献[11]。

另一个问题就是波动方程正演的线性化(即用反偏移代替正演)带来的误差。在数据观测较多, 成像系统是超定系统时, 最小二乘的方法可以自动地将正交于反偏移算子的列空间的噪声排除在解空间之外; 然而, 随着降采样带来的系统欠定性增强, 噪声可能会在解空间投影, 从而给真实解带来干扰。TU等^[9]详细分析了线性误差对成像系统带来的影响, 并提出, 根据特定的规则在不同的迭代之间引入不同的降采样算子(即对数据和炮集重新进行采样)会显著提高成像系统对线性噪声的鲁棒性。

利用SEG/EAGE盐丘模型的一个二维剖面, 对上述算法的有效性作了验证^[10]。实验模型的参数如下:模型3.9km深, 15.7km长, 网格间距约为24m。共设置323个炮点和检波点, 利用互易原理, 将数据设置为每一炮都有共同的检波器来接收。真实速度模型及其平滑之后的背景模型如图 1所示。图 2展示了一系列实验结果。首先, 利用线性数据(即利用反偏移算子作用到真实模型扰动上生成的地震数据)分析无任何噪声的状况下所提算法的有效性。图 2a为对线性数据做单次逆时偏移的结果; 图 2b为对线性数据做快速最小二乘偏移的结果。在快速最小二乘偏移中, 所用合成炮的数量约为原炮数的5%, 所用频率子集的大小约为总频率集合的15%, 我们运行60个SPGL1的迭代, 在这样的反演参数设置下, 总的波动方程的求解数量相当于利用所有数据进行单次逆时偏移, 即图 2a和图 2b所展示的结果包含几乎相同的波动方程求解。接着, 验证存在相干噪声的情况下该算法的效果。用时间域的iWave波场正演引擎来模拟观测波场^[16], 用我们的频率域波场线性反偏移引擎来进行成像, 结果如图 2c所示。这两种波场模拟引擎间存在数值算法上的差别, 同时, iWave正演的观测波场含有高阶散射。在这种情况下, 快速最小二乘成像方法依然取得了较好的效果。和图 2a相比, 图 2b和图 2c均展现出了较高的分辨率。

1.2 震源子波估计

无论是单次逆时偏移还是最小二乘逆时偏移, 都需要预先知道地震子波。如果地震子波估算不准, 那么正传波场就会有误, 从而影响地震成像以及最小二乘偏移中梯度的求取。求取可靠的地震子波在技术上或者计算上面临较大的挑战:依赖统计信息^[17]或者大量的计算^[18]。最小二乘偏移本身是一种反演算法, 因此我们将子波的反演融入其中, 在不给出子波先验信息的情况下, 得到可靠的地震成像结果。

借助在非线性最小二乘问题中变量投影(Variable Projection)的方法^[19-20], 并将其理论推广到包含稀疏约束的非线性最小二乘问题中, 提出了不依赖于子波信息的最小二乘逆时偏移算法^[10]。对于子波w来说, 一旦x确定, 存在下述解析表达:

$ {{\tilde w}_i}\left( \mathit{\boldsymbol{x}} \right) = \frac{{\sum\limits_{j \in \mathit{\boldsymbol{ \boldsymbol{\varSigma} }}} {\left\langle {{\mathit{\boldsymbol{d}}_{i,j}},\nabla \mathit{\boldsymbol{F}}\left[ {{\mathit{\boldsymbol{m}}_0},{\mathit{\boldsymbol{s}}_j}} \right]{\mathit{\boldsymbol{C}}^*}\mathit{\boldsymbol{x}}} \right\rangle } }}{{\sum\limits_{j \in \mathit{\boldsymbol{ \boldsymbol{\varSigma} }}} {\left\langle {\nabla \mathit{\boldsymbol{F}}\left[ {{\mathit{\boldsymbol{m}}_0},{\mathit{\boldsymbol{s}}_j}} \right]{\mathit{\boldsymbol{C}}^*}\mathit{\boldsymbol{x}},\nabla \mathit{\boldsymbol{F}}\left[ {{\mathit{\boldsymbol{m}}_0},{\mathit{\boldsymbol{s}}_j}} \right]{\mathit{\boldsymbol{C}}^*}\mathit{\boldsymbol{x}}} \right\rangle } }} $

(2)

使用交替求解的方式, 每更新一次x, 就更新一次w的估计, 即可在未给定子波的情况下求解最小二乘偏移的梯度^[10]。

继续利用前文所用的模型进行合成数据的实验, 主要结果如图 3所示。图 3a展示了利用线性化数据进行最小二乘逆时偏移的结果, 但正演中使用的地震子波和真实的子波存在一个相位差, 其影响是反射层发生错位(如箭头所示)和整体地震像散焦。图 3b展示了利用相同的线性化数据, 在未给定子波(第一个迭代中的子波起始估计是白谱子波，即时间域的脉冲)的情况下, 利用本文包含子波估计的成像方法得到的地震像。图 3c展示了利用iWave正演的数据进行实验的结果(其它参数与图 3b所用的一致)。从结果来看, 本文方法可以在未知真实子波的情况下, 得到可靠的地震像。对比图 3b, 图 3c和图 2b, 图 2c可以看出, 采用本文方法得到的地震像与知道真实子波情况下得到的地震像具有相似的高精度和分辨率。

1.3 表面多次波成像

海洋表面多次波是海洋地震资料处理时面临的强干扰。逆时偏移本身基于单次散射的假设, 并不能对多次波进行有效成像。我们从描述自由表面反射的自由表面多次波消除公式出发, 研究出对多次波进行成像以及一次波和多次波联合成像的方法^[10-11]。海洋表面多次波的正确成像依赖于:①正确识别表面多次波对应的震源波场。一次波的成像中, 正传波场为点震源波场; 多次波的成像中, 正传波场是海洋表面的整体下行波场, 可以看作一个面震源^{[11, 21-23]}。②新的成像条件。对于含有多阶海洋表面反射的多次波来说, 广泛应用于一次波成像的互相关成像条件会造成各阶次反射之间产生非物理的串扰噪声^[24-25], 形成对正确地震像的干扰。应用反褶积成像条件可以一定程度缓解这个问题^{[24, 26]}, 而真正衰减串扰噪声则需要利用最小二乘反演的方法^{[11, 27]}。

基于我们提出的快速最小二乘逆时偏移的框架, 实现对表面多次波的正确成像。通过对正传震源波场做修正, 用自由表面的下行波场代替点震源波场, 可实现对表面多次波的正确成像^[11]以及一次波和多次波的联合成像^{[10, 28]}。数据的预处理细节参考文献[10]。本文主要介绍一次波和多次波联合成像的方法。这时, 数据匹配的目标变为:

$ \sum\limits_{i \in \mathit{\Omega }} {\sum\limits_{j \in \mathit{\boldsymbol{ \boldsymbol{\varSigma} }}} {\left\| {{\mathit{\boldsymbol{d}}_{i,j}} - \nabla \mathit{\boldsymbol{F}}\left[ {{\mathit{\boldsymbol{m}}_0},{w_i}{\mathit{\boldsymbol{s}}_j} - {\mathit{\boldsymbol{d}}_{i,j}}} \right]{\mathit{\boldsymbol{C}}^*}\mathit{\boldsymbol{x}}} \right\|_2^2} } $

(3)

需要注意的是, 此处的观测数据d同时包含了上行的一次波波场和多次波波场。对于预测数据, 可以看到, 在反偏移算子中, 震源波场w_is_j-d_{i, j}不只是包含点震源s(经过波场传播后预测一次波), 同时包含广义的面震源波场即观测数据d(前面的“-”表示水面的反射系数为-1, 上行的观测数据经过水面反射, 再经过和一次波相同的波场传播后, 预测表面多次波)。除去震源波场的变化, 其它成像过程与前文所述一次波的反演成像过程相同。另外, 这种情况下, 由于震源波场的变化, 子波的估计方法则为:

$ \begin{array}{l} {{\tilde w}_i}\left( \mathit{\boldsymbol{x}} \right) = \\ \frac{{\sum\limits_{j \in \mathit{\boldsymbol{ \boldsymbol{\varSigma} }}} {\left\langle {{\mathit{\boldsymbol{d}}_{i,j}} - \nabla \mathit{\boldsymbol{F}}\left[ {{\mathit{\boldsymbol{m}}_0},{\mathit{\boldsymbol{d}}_{i,j}}} \right]{\mathit{\boldsymbol{C}}^*}\mathit{\boldsymbol{x}},\nabla \mathit{\boldsymbol{F}}\left[ {{\mathit{\boldsymbol{m}}_0},{\mathit{\boldsymbol{s}}_j}} \right]{\mathit{\boldsymbol{C}}^*}\mathit{\boldsymbol{x}}} \right\rangle } }}{{\sum\limits_{j \in \mathit{\boldsymbol{ \boldsymbol{\varSigma} }}} {\left\langle {\nabla \mathit{\boldsymbol{F}}\left[ {{\mathit{\boldsymbol{m}}_0},{\mathit{\boldsymbol{s}}_j}} \right]{\mathit{\boldsymbol{C}}^*}\mathit{\boldsymbol{x}},\nabla \mathit{\boldsymbol{F}}\left[ {{\mathit{\boldsymbol{m}}_0},{\mathit{\boldsymbol{s}}_j}} \right]{\mathit{\boldsymbol{C}}^*}\mathit{\boldsymbol{x}}} \right\rangle } }} \end{array} $

(4)

选用专为验证与表面多次波有关的算法而设计的Sigsbee 2B模型, 对包含海洋表面多次波的数据进行了测试。真实和背景速度模型见图 4。数据正演使用iWave, 使用自由表面边界条件来正演包含多次波和鬼波的观测数据。假设真实地震子波已知, 利用本文算法对一次波和多次波进行联合快速最小二乘偏移成像, 结果如图 5a所示, 其整体的波动方程模拟次数依然控制在单次逆时偏移的水平上。图 5b展示了通过子波估计对一次波和多次波进行联合成像的结果。可见, 图 5a和图 5b结果具有很高的一致性, 与真实模型对比可以看出, 这两个结果中多次波带来的串扰噪声均得到了很好的压制。文献[10]中详细讨论了本文方法对子波初始估计的鲁棒性, 本文方法在子波具有错误的初始振幅和相位的情况下, 均能得到高精度的地震像。

2 实际数据实验

利用英国北海海域某工区实际数据, 对文中所述算法进行验证。图 6a为利用互相关成像条件对多次波进行偏移成像的结果(成像过程中, 震源为面震源, 由表面位置的下行波场构成, 反传数据为多次波), 可以看到不同阶次多次波导致的串扰噪声。图 6b为利用本文快速最小二乘偏移算法对多次波成像的结果。图 7和图 8分别为图 6中区域A和区域B放大后的结果。可以看出, 图 6a中出现的极强的串扰噪声在图 6b中得到了很好的压制。

图 9为不同数据的偏移成像结果。图 9a是一次波成像结果; 图 9b是多次波单独成像结果; 图 9c是一次波与多次波联合成像的结果。可以看出, 多次波单独成像和一次波成像存在较好的可比性, 较强的连续反射层(红色框中的反射层)均能得到较高精度的成像, 但是整体信噪比较低, 有些反射层(黑色箭头所指)没有成像; 对于某些一次波成像没能揭示出的结构(白色箭头所指), 多次波成像能较好的刻画出来。从图 9c可以看出, 一次波与多次波的联合成像, 不仅在数据处理环节避免了多次波和一次波的分离, 而且改进了多次波成像本身信噪比偏低的缺点。

3 结论

本文将压缩感知应用于地震成像中, 通过人工压缩观测系统和观测数据, 大幅降低最小二乘逆时偏移的成本。快速最小二乘逆时偏移是个灵活的框架, 在这个框架下, 介绍了如何利用变量投影技术解决地震偏移中未知子波的问题, 并利用广义面震源实现了海洋多次波的成像以及一次波与多次波的联合成像。模型数据以及海上实际数据应用结果表明, 本文方法在未知地震子波信息的情况下, 利用相当于单次逆时偏移的计算量, 可以得到高分辨率的最小二乘偏移成像结果, 并对多次波进行正确的归位和成像, 为今后更加有效和多样化地利用多次波信息提供了思路。