工业界常用的基于伴随算子的地震偏移成像算法假设地震数据中只存在一次反射波(一次波), 而多次反射波(多次波)的存在不符合传统成像算法的假设条件, 造成有效波振幅、频率、相位畸变, 影响地震解释^[1]。多次波广泛存在于陆地和海洋地震数据中, 特别是在海洋地震数据中, 由于存在海面和海底两个强反射界面, 因而导致海洋地震数据的表面多次波干扰问题更为严重。对于陆上地震数据而言, 当地下存在强反射界面, 诸如基岩面、火成岩侵入体、不整合面等, 地震波在强反射界面间或其内部发生多次反射, 形成层间多次波^[2]。由于多次波传播路径更加复杂, 因而压制层间多次波更有挑战性。多次波压制方法主要分为两类: 滤波法和预测相减法^[3]。滤波类方法利用一次波和多次波在变换域内的差异压制多次波, 包括F-K滤波^[4]、Radon变换滤波方法^[5-8]等。通常当一次波和多次波速度差异较大时, 滤波类方法效果较理想; 但当地下介质复杂时, 滤波类方法很难在保护有效信号的同时压制多次波。预测相减法首先基于波动理论得到预测多次波模型, 然后通过自适应相减压制多次波。这类算法更加适用于复杂介质情况。BERKHOUT等^[9]与VERSCHUUR等^[10]基于反馈迭代模型, 提出了表面多次波去除方法(SRME), 该方法已经成为工业上成熟应用的表面多次波去除方法。针对层间多次波的复杂波场特征, 诸多国内外地球物理研究人员发展了相应的层间多次波压制方法, 包括逆散射级数法^{[7, 11-13]}、共聚焦点技术^[14-15]和Marchenko层间多次波压制技术^[16-17]等。

但是, 在实际地震数据采集过程中, 观测系统往往存在缺陷, 震源子波常常不一致等。这些因素使得预测多次波模型往往存在相位、频带、振幅的差异^[9]。因此, 自适应相减方法作为预测相减法不可或缺的关键步骤, 能够校正预测多次波模型中的误差。工业界最常用的基于最小二乘匹配滤波器的多次波自适应相减方法^[10]利用预测多次波直接匹配观测数据, 以去除多次波后的数据的L₂范数最小为优化目标。该方法效率很高, 可以快速估计出观测数据中的多次波。但是其仅在一次波和多次波正交假设下, 才能取得较好的效果。随后, 基于L₁范数的匹配滤波器方法被提出^[18]。该方法可以更好地处理一次波和多次波不正交的数据。独立成分分析方法也被应用于多次波自适应相减问题^[19-20]。

基于模式学习的多次波自适应相减方法相继被提出^[21-24]。模式学习方法首先从预测多次波中学习多次波特征, 然后利用学习到的多次波特征去重构原始数据, 这样就可以得到原始数据中的多次波, 进而将其从原始数据中减去即可得到一次波数据。模式学习方法将多次波衰减分为多次波模式学习和自适应相减两个独立过程, 能够更好地保护有效信号。本文基于混合L₂-L₁范数模式学习方法^[24], 衰减地震数据中的多次波。同时, 针对预测多次波模型中存在的相位和频带误差, 引入伪地震数据^[25-28], 在不损伤有效信号的同时, 进一步减少多次波残留。最后利用模拟数据和实际数据验证本文方法的有效性。

采用SRME预测多次波模型m∈ $ \mathbb{R}$^h×w, 其中, w表示地震道数, h表示地震道采样点数。我们利用主成分分析(Principle Component Analysis, PCA)方法来提取预测多次波的关键模式。首先将预测多次波模型m划分为n个互相重叠的大小为d=k×l的图像块, 再将这些图像块向量化之后合成一个矩阵M∈$ \mathbb{R}$^d×n, 其中, n=(h－k+1)×(w－l+1)。对该矩阵进行PCA特征编码, 有:

式中: B∈$ \mathbb{R}$^d×p表示学习到的p个d维特征向量构成的特征矩阵; X∈$ \mathbb{R}$^p×n表示对应的特征系数。特征矩阵B的构建过程如下。

算法1:SVD求解公式。

1) 初始化: 通常λ=0.1。

2) 对数据矩阵M作SVD分解, 有[U, Σ, V]=SVD(M)。

3) 设Σ=diag({σ_i}_1≤i≤N), 选择前p个最大的特征值, 使得σ_p+1 < λσ₁, 且σ_p≥λσ₁。

4) 设U=[u₁, u₂, …, u_n], 选择前p个最大的特征值对应的特征向量, 构成特征矩阵B=[u₁, u₂, …, u_p]。

考虑分频编码策略, 用带通滤波器将预测多次波模型分为低、中、高3个不同的频段, 即m_LF, m_MF, m_HF。相应的低、中、高频数据矩阵记为M_LF, M_MF, M_HF。利用算法1学习分频段特征矩阵B_LF, B_MF, B_HF。

最终的分频特征矩阵为:

由此, 我们可以得到预测多次波m的特征矩阵B^[24]。在完成模式学习并学习到多次波的特征之后, 我们便可以利用学习到的特征去重构观测数据中包含的多次波, 进而分离出观测数据的一次波, 即自适应相减过程。本质上是估计一组系数, 使得多次波关键模式线性组合后, 与观测数据的差异最小。

与模式学习过程相对应, 我们将观测数据s(包含多次波和一次波)划分为若干个大小为k×l的重叠图像块, 再将这些图像块向量化之后组合成一个矩阵S∈$ \mathbb{R}$^d×n。在数据矩阵S中, 每个列向量对应一个原数据s中的图像块。我们用L₁范数来重构多次波^[24]:

令P=S－BY即一次波数据矩阵, 其中P∈$ \mathbb{R}$^d×n, 公式重写为:

由拉格朗日乘数法, 引入一个罚函数, 我们将公式(4)转化为一个无约束的优化问题:

其中, μ为惩罚系数, 且μ>0。当μ较小时, 公式(5)求解不准确; 当μ较大时, 公式(5)求解数值不稳定。为避免μ调参复杂, 我们利用裂步Bregman算法^[29]求解公式(5)。引入中间变量A^k, 且A¹=0, 则:

在裂步Bregman算法^[29]中, 公式(7)用来添加约束, 避免选择一个较大的μ带来的数值不稳定问题。通常, 我们选择μ=1。公式(6)两个变量P和Y是解耦的, 我们可以利用交替下降法近似求解公式(6):

公式(8)和公式(9)的解析解为:

其中, prox_1/μ¹(y)=sign(y)max(|y|－1/μ, 0), 即为收缩阈值算子。

裂步Bregman算法^[29]求解自适应相减问题(公式(3))的步骤总结如下。

算法2:裂步Bregman算法^[29]求解公式。

1) 初始化: k=1, A¹=0, P¹=0, P⁰=S, D=(B^TB)^－1。

2) While‖P^k－P^k－1‖₂>γ do。

3) Y^k+1=DB^T(S+A^k－P^k)。

4) P^k+1=prox_1/μ¹(S+A^k－BY^k+1)。

5) A^k+1=A^k+(S－BY^k+1－P^k+1)。

6) k=k+1。

7) End while。

8) 输出: $ \mathit{\boldsymbol{\tilde P}} = {\mathit{\boldsymbol{P}}^k}, \mathit{\boldsymbol{\tilde M}} = \mathit{\boldsymbol{B}}{\mathit{\boldsymbol{Y}}^k}$。

算法2主要包括3步, 即公式(7)、公式(10)和公式(11)。其中, $ \mathit{\boldsymbol{\tilde P}}$和$ \mathit{\boldsymbol{\tilde M}}$分别表示估计的一次波和多次波。综上所述, 模式学习方法将多次波衰减分为多次波模式学习和多次波自适应相减两个独立过程, 能够更好地保护有效信号。但是, 由于预测多次波模型中存在相位和频带误差^[9], 因此, 估计的一次波中会存在多次波残留。我们进一步引入伪地震数据, 提供多次波模型额外的相位和频带信息, 以实现更有效的多次波衰减。

MONK^[25]考虑到预测多次波的相位和频带误差, 对预测多次波模型道m进行数学变换, 生成伪地震数据, 拓展了传统的单道匹配相减算法。数学变换具体包括预测多次波模型道m的梯度道$ {\mathit{\boldsymbol{\dot m}}}$, Hilbert变换道m^H及其梯度道$ {{\mathit{\boldsymbol{\dot m}}}^{\rm{H}}}$^[25-28]。伪地震数据引入了额外的多次波模型相位和振幅信息, 能够减少多次波残留。在模式学习阶段, 我们先得到预测多次波模型道m及其梯度道$ {\mathit{\boldsymbol{\dot m}}}$, Hilbert变换道m^H及其梯度道${{\mathit{\boldsymbol{\dot m}}}^{\rm{H}}} $构成的数据矩阵M, ${\mathit{\boldsymbol{\dot M}}} $, M^H, $ {{\mathit{\boldsymbol{\dot M}}}^{\rm{H}}}$。利用算法1得到相应的特征矩阵B_m, $ {{\mathit{\boldsymbol{\dot B}}}_m}$, B_m^H, $ \mathit{\boldsymbol{\dot B}}_m^{\rm{H}}$。最终的特征矩阵为:

在得到特征矩阵B后, 利用算法2重构观测数据中的多次波, 进而分离一次波和多次波。图 1展示了基于伪地震数据的多次波模式学习和自适应相减方法流程。

图 1(Fig. 1)

图 1 基于伪地震数据的多次波模式学习和自适应相减方法流程

采用Sigsbee模拟数据验证算法的有效性^[23]。图 2a, 图 2b和图 2c分别展示了真实一次波、预测多次波模型和观测数据。我们对真实多次波进行如下处理生成预测多次波模型: 纵向错动3个网格, 横向错动2个网格, 然后与主频20Hz, 长30个单位网格的雷克子波进行褶积。图 3a, 图 3b和图 3c分别展示了采用分频模式学习方法估计的一次波、多次波和误差。图 4a, 图 4b和图 4c分别展示了采用伪地震数据模式学习方法估计的一次波、多次波和误差。对比图 3c和图 4c可以发现, 伪地震数据模式学习方法误差更小, 在未损伤一次波的基础上, 多次波残留更少。定量来看, 分频模式学习方法和伪地震数据模式学习方法估计的一次波信噪比分别为19.03dB和20.06dB, 伪地震数据模式学习方法的信噪比提升了1dB。图 2中箭头所指位置强能量多次波造成分频模式学习方法存在多次波残留, 而伪地震数据模式学习方法明显残留更少。为了更清晰地对比, 我们将图 2至图 4中红色框部分放大显示(图 5至图 7)。由于绕射多次波预测不准, 分频模式学习方法残留绕射多次波严重。伪地震数据模式学习方法几乎完全压制了绕射多次波。因此, 相比于分频模式学习方法, 伪地震数据模式学习方法在不损伤一次波的基础上, 能够有效避免多次波残留, 特别是对绕射多次波压制效果更佳。

图 2(Fig. 2)

图 2 Sigsbee模拟数据 a 真实一次波; b 预测多次波; c 观测数据

图 3(Fig. 3)

图 3 采用分频模式学习方法估计的一次波(a)、多次波(b)和误差(c)

图 4(Fig. 4)

图 4 采用伪地震数据模式学习方法估计的一次波(a)、多次波(b)和误差(c)

图 5(Fig. 5)

图 5 图 2红色框位置的放大显示结果 a 真实一次波; b 预测多次波; c 观测数据

图 6(Fig. 6)

图 6 图 3红色框位置的放大显示结果 a 分频模式学习估计的一次波; b 分频模式学习估计的多次波; c 分频模式学习的误差

图 7(Fig. 7)

图 7 图 4红色框位置的放大显示结果 a 伪地震数据模式学习估计的一次波; b 伪地震数据模式学习估计的多次波; c 伪地震数据模式学习的误差

我们对海上实际地震资料进行了多次波压制。图 8a和图 8b分别显示了观测数据和预测多次波模型; 图 9a和图 9b分别显示了采用分频模式学习方法估计的一次波和多次波; 图 10a和图 10b分别显示了采用伪地震数据模式学习方法估计的一次波和多次波。对比图 8至图 10中箭头所指处可以看出, 强能量多次波导致分频模式学习方法存在一个低频残留; 而伪地震数据模式学习方法不存在这一问题。放大显示图 8至图 10中的矩形框1和矩形框2, 如图 11至图 16所示。图 11至图 13中箭头所指位置尖灭构造生成复杂多次波。由于多次波预测过程中多次褶积的影响, 存在相位误差。分频模式学习方法存在一个明显和多次波走势相同的残留。伪地震数据模式学习方法减少了多次波残留。图 14至图 16中箭头所指位置绕射多次波发育。但由于大角度绕射多次波预测误差较大, 导致分频模式学习方法并没有很好地压制掉大角度的绕射多次波, 而伪地震数据模式学习方法能够缓解这一问题。但是对于更大角度处, 预测多次波时移误差更大, 仍然存在部分残留。

图 8(Fig. 8)

图 8 海上实际地震资料 a 观测数据; b 预测多次波

图 9(Fig. 9)

图 9 采用分频模式学习方法估计的一次波(a)和多次波(b)

图 10(Fig. 10)

图 10 采用伪地震数据模式学习方法估计的一次波(a)和多次波(b)

图 11(Fig. 11)

图 11 图 8中矩形框1的放大显示 a 观测数据; b 预测多次波

图 12(Fig. 12)

图 12 图 9中矩形框1的放大显示 a 分频模式学习估计的一次波; b 分频模式学习估计的多次波

图 13(Fig. 13)

图 13 图 10中矩形框1的放大显示 a 伪地震数据模式学习估计的一次波; b 伪地震数据模式学习估计的多次波

图 14(Fig. 14)

图 14 图 8中矩形框2的放大显示结果 a 观测数据; b 预测多次波

图 15(Fig. 15)

图 15 图 9中矩形框2的放大显示结果 a 分频模式学习估计的一次波; b 分频模式学习估计的多次波

图 16(Fig. 16)

图 16 图 10中矩形框2的放大显示结果 a 伪地震数据模式学习估计的一次波; b 伪地震数据模式学习估计的多次波

本文提出了一种基于伪地震数据的多次波模式学习和自适应相减方法。该方法一方面能够继承模式学习方法的优势, 更好地保护有效信号; 另一方面, 利用了伪地震数据中额外的相位和高频信息, 避免了多次波残留。Sigsbee模拟数据实验结果表明, 相比于分频模式学习方法, 本文提出的伪地震数据模式学习方法一次波信噪比能够提升1dB。海上实际数据应用结果表明, 分频模式学习方法在多次波复杂情况下效果不佳, 本文方法能够更好地消除大角度的绕射多次波。但是当绕射多次波预测时移误差更大时, 仍然存在一定残留。事实上, 可以通过进一步调参来消除更大时移误差, 但会损伤部分有效信号。因此, 对于实际地震数据, 多次波预测误差变化较大时, 必须要在保护一次波和压制多次波之间进行折衷处理。