地震勘探中多次波问题普遍存在, 多次波的存在使得一次反射波的相位、振幅和频率等动力学特征失真, 对构造解释和地震反演都造成了不利影响。当一次波和多次波相互重叠时, 会更加难以识别和处理多次波。因此, 在地震资料处理中, 有效压制多次波一直是地震处理研究的重点和难点之一。

BERKHOUT^[1]提出了描述复杂多次波系统的反馈理论框架, 奠定了反馈迭代多次波压制方法的数学物理基础, 提出的多维反演算法通过数据矩阵相乘预测多次波, 能够适应任意复杂的地下结构。但消除多次波时需要已知震源子波。VERSCHUUR^[2]引入基于一次波能量最小假设的自适应相减方法, 成功地从实际地震数据中估计出震源子波, 极大地推动了表面相关多次波压制方法(surface-related multiple elimination, SRME)的发展。为了避免因求解震源子波而产生的非线性优化问题, VERSCHUUR等^[3]提出了迭代SRME方法, 在每次迭代过程中采用最小二乘匹配多次波模型与原始数据, 确定震源子波, 将原来的非线性问题转化为线性问题, 增强了方法的实用性。为了避免自适应相减可能损伤有效信号, SRME被重新定义为反演一次波脉冲响应和震源子波的全波形反演问题, 形成了基于稀疏反演的一次波估计方法(estimation of primaries by sparse inversion, EPSI)^[4-5]。

传统的EPSI方法需要精确地包含有一次波的时窗和反演参数设定, 存在一定的不稳定性。因此, 需要对EPSI方法进行改进, 将其转化成为双凸优化问题^[6], 即反演一次波脉冲响应过程中采用的目标函数是一个凸函数, 应用的约束集合同样也是一个凸函数, 通过交替优化迭代过程, 对一次波脉冲响应采用谱梯度投影(SPGL1)^[7]算法求解, 同时引入帕累托优化(pareto)曲线。改进的EPSI方法^[8-10]在求解过程中, 需要在输入记录的多维互相关数据中拾取与最强反射一次波相一致的同相轴来估计初始的震源子波。该工作不仅费时费事, 还可能由于识别错误和拾取参数设置误差, 得不到合适的震源子波。冯飞^[11]使用常规最小平方正交三角(QR)分解算法求解震源子波。宋家文等^[12]和白兰淑等^[13]通过匹配滤波方法求取地震子波。本文提出的自适应稀疏反演多次波压制方法(adaptive estimation of primaries by sparse inversion, AEPSI), 不需要人工拾取和时窗设置, 完全由数据驱动解决初始震源子波及震源子波迭代更新问题。

1 方法原理 1.1 自适应稀疏反演多次波压制原理

为了便于理解和推导, EPSI频率域方程表示为:

$ \mathit{\boldsymbol{P}} = \mathit{\boldsymbol{GQ}} + \mathit{\boldsymbol{GRP}} $

(1)

式中:P为频率域数据矩阵; G为频率域一次波脉冲响应; R为频率域自由表面反射算子; Q为频率域震源矩阵。如果用q(t)表示对应于时间域数据矩阵p的震源特性, 并假定q(t)对所有的炮和道是不变的, 则Q=w(ω)I。其中, w(ω)是频率的标量函数, 是q(t)的频率域表示; I是单位矩阵。假定自由表面反射算子R=－I。方程通过逆傅里叶变换返回到时间域, G的时间域表示为g, P的时间域表示为p, 时间域子波q(t)简记为q, 则时间域数据矩阵p可表示为:

$ \mathit{\boldsymbol{p}} = \mathit{\boldsymbol{M}}(\mathit{\boldsymbol{g}},\mathit{\boldsymbol{q}};\mathit{\boldsymbol{p}}) $

(2)

由于稀疏约束是施加在时间域, 所以用作用在时间域的M(g, q; p)函数讨论问题会更方便。如果基于数据矩阵p的方程(2)用隐式表示, 则有:

$ \mathit{\boldsymbol{p}} = \mathit{\boldsymbol{M}}(\mathit{\boldsymbol{g}},\mathit{\boldsymbol{q}}) $

(3)

对于方程(3)来说, 如果给定变量g、q中的一个, 则M(g, q)变成了相对另一个变量的线性算子。因此, M(g, q)也叫做双线性算子, 可以写成:

$ \begin{array}{*{20}{l}} {\mathit{\boldsymbol{M}}(\mathit{\boldsymbol{g}},\mathit{\boldsymbol{q}}): = \mathit{\boldsymbol{F}}_t^{\rm{H}}{\rm{BlockDia}}{{\rm{g}}_\omega }[(w(\omega )\mathit{\boldsymbol{I}} - }\\ {{\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} \mathit{\boldsymbol{ P}}{)^{\rm{T}}} \otimes \mathit{\boldsymbol{I}}]{\mathit{\boldsymbol{F}}_t}\mathit{\boldsymbol{g}} = \mathit{\boldsymbol{p}}} \end{array} $

(4)

式中:$ \otimes $表示Kronecker积; F_t是傅里叶变换算子, 表示在时间域进行傅里叶变换, 也就是变换为频率切片, 可对任何时间域信号进行处理; F_t^H为F_t共轭算子, 作用是将频率域变换到时间域。其中, F_tg表示对时间域一次波脉冲响应g进行傅里叶变换。可以看出, M(g, q)是一个包含有震源子波w(ω)和总的上行波场P的线性算子。

假定变量g的估计值为$ \mathit{\boldsymbol{\tilde g}}$, 则得到关于q的线性算子:

$ {\mathit{\boldsymbol{M}}_{\mathit{\boldsymbol{\tilde g}}}} = \mathit{\boldsymbol{M}}(\mathit{\boldsymbol{\tilde g}},\mathit{\boldsymbol{q}}) $

(5)

假定变量q的估计值为$ \mathit{\boldsymbol{\tilde q}}$, 则得到关于g的线性算子:

$ {\mathit{\boldsymbol{M}}_{\mathit{\boldsymbol{\tilde q}}}} = \mathit{\boldsymbol{M}}(\mathit{\boldsymbol{g}},\mathit{\boldsymbol{\tilde q}}) $

(6)

根据(5)式和(6)式, 定义:

$ {\mathit{\boldsymbol{M}}_{\mathit{\boldsymbol{\tilde g}}}}\mathit{\boldsymbol{q}}: = {\left( {\frac{{\partial \mathit{\boldsymbol{M}}}}{{\partial \mathit{\boldsymbol{q}}}}} \right)_{\mathit{\boldsymbol{\tilde g}}}}\mathit{\boldsymbol{q}} = \mathit{\boldsymbol{M}}(\mathit{\boldsymbol{\tilde g}},\mathit{\boldsymbol{q}}) $

(7a)

$ {\mathit{\boldsymbol{M}}_{\mathit{\boldsymbol{\tilde q}}}}\mathit{\boldsymbol{g}}: = {\left( {\frac{{\partial \mathit{\boldsymbol{M}}}}{{\partial \mathit{\boldsymbol{g}}}}} \right)_{\mathit{\boldsymbol{\tilde q}}}}\mathit{\boldsymbol{g}} = \mathit{\boldsymbol{M}}(\mathit{\boldsymbol{g}},\mathit{\boldsymbol{\tilde q}}) $

(7b)

(7a)式中左边表示已知变量g的估计值$ \mathit{\boldsymbol{\tilde g}}$, 则p=M(g, q)是关于q的显式线性关系。(7b)左边表示已知变量q的估计值$ \mathit{\boldsymbol{\tilde q}}$, 则p=M(g, q)是关于g的显式线性关系。因此, 可使用新的符号来表示EPSI过程。它是求解下列标准形式约束优化问题的方程:

$ \mathop {{\rm{minimize}}}\limits_{\mathit{\boldsymbol{g}},\mathit{\boldsymbol{q}} \in \varLambda } {\left\| {\mathit{\boldsymbol{p}} - \mathit{\boldsymbol{M}}(\mathit{\boldsymbol{g}},\mathit{\boldsymbol{q}})} \right\|_2}\quad {\rm{ s}}{\rm{.t}}{\rm{. }}{\left\| \mathit{\boldsymbol{g}} \right\|_0} \le \rho $

(8)

式中:q∈Λ表示满足任何所需的或先验约束施加在震源上。这里假定q在时间域是短时窗的震源子波, 并允许非因果的描述记录数据中可能出现的时移。当然, 也可以施加其它需要的震源特性。

根据迭代求解约束优化问题的方法, 则可以通过下列两个子问题的迭代求解EPSI约束优化问题:

$ \begin{array}{*{20}{c}} {{{\mathit{\boldsymbol{\tilde g}}}_{k + 1}} \leftarrow \mathop {{\rm{ argmin }}}\limits_\mathit{\boldsymbol{g}} {{\left\| \mathit{\boldsymbol{g}} \right\|}_1}\quad {\rm{ s}}{\rm{.t}}{\rm{. }}\quad \left\| {\mathit{\boldsymbol{p}} - } \right.}\\ {{{\left. {{\mathit{\boldsymbol{M}}_{{{\mathit{\boldsymbol{\tilde q}}}_k}}}\mathit{\boldsymbol{g}}} \right\|}_2} \le \sigma } \end{array} $

(9)

$ {\mathit{\boldsymbol{\tilde q}}_{k + 1}} \leftarrow \mathop {{\rm{argmin}}}\limits_q \left\| {\mathit{\boldsymbol{p}} - {\mathit{\boldsymbol{M}}_{{{\mathit{\boldsymbol{\tilde g}}}_{k + 1}}}}\mathit{\boldsymbol{q}}} \right\|_2^2\quad {\rm{ s}}{\rm{.t}}{\rm{. }}\quad \mathit{\boldsymbol{q}} \in \varLambda $

(10)

假定已知震源子波$ \mathit{\boldsymbol{\tilde q}}$_k, 则可计算M_{$ \mathit{\boldsymbol{\tilde q}}$k}。(9)式是一个基追踪降噪(basis pursuit de-noising, BPDN)问题, 表示成L₁范数约束问题(通常称为Lasso问题)为:

$ \mathop {{\rm{argmin}}}\limits_\mathit{\boldsymbol{g}} {\left\| {\mathit{\boldsymbol{p}} - {\mathit{\boldsymbol{M}}_{{{\mathit{\boldsymbol{\tilde q}}}_k}}}\mathit{\boldsymbol{g}}} \right\|_2}\quad {\rm{ s}}{\rm{.t}}{\rm{. }}{\left\| \mathit{\boldsymbol{g}} \right\|_1} \le {\tau _k} $

(11)

这里τ_k是L₁范数约束, 使用SPGL1方法求解。Lasso问题的解σ(τ_k)会快速下降到Pareto曲线, 如果有一系列明确的约束τ₀ < τ₁ < τ₂ < … < τ_k, 那么就能收敛到最小约束并使目标函数达到选定的误差水平, 最终收敛到BPDN问题的解。应用牛顿求根法, 在BPDN问题的Pareto曲线上, 可以找到约束系列τ_k。每一个Lasso问题σ(τ_k), 对应得到一个一次波脉冲响应$ \mathit{\boldsymbol{\tilde g}}$_k, 并且$ \mathit{\boldsymbol{\tilde g}}$_k是稀疏的。

得到$ \mathit{\boldsymbol{\tilde g}}$_k后, 进行震源子波的估计更新, 以便计算M_{$ \mathit{\boldsymbol{\tilde q}}$k}, 用于下一次迭代求解Lasso问题。Lasso问题施加了L₁范数约束在一次波脉冲响应$ \mathit{\boldsymbol{\tilde g}}$_k上, 但没有直接限制$ \mathit{\boldsymbol{\tilde g}}$_k中一次波反射同相轴的振幅。由于问题的稀疏性, 反射同相轴的振幅可能被低估了。因此, 为了使估计的震源子波具有正确的振幅, 需要重新调整$ \mathit{\boldsymbol{\tilde g}}$_k的振幅。假定记录中的同相轴具有同样的震源特性, 用单一因子s_k就能调整$ \mathit{\boldsymbol{\tilde g}}$_k满足:

$ \mathop {{\rm{argmin}}}\limits_\mathit{\boldsymbol{q}} {\left\| {\mathit{\boldsymbol{p}} - {\mathit{\boldsymbol{M}}_{{{\mathit{\boldsymbol{\tilde q}}}_{k - 1}}}}({s_k}{{\mathit{\boldsymbol{\tilde g}}}_k})} \right\|_2} $

(12)

调整因子s_k为:

$ {s_k} = \frac{{\mathit{\boldsymbol{\tilde g}}_k^{\rm{H}}\mathit{\boldsymbol{M}}_{{{\mathit{\boldsymbol{\tilde q}}}_{k - 1}}}^{\rm{H}}\mathit{\boldsymbol{p}}}}{{\left\| {{\mathit{\boldsymbol{M}}_{{{\mathit{\boldsymbol{\tilde q}}}_{k - 1}}}}{{\mathit{\boldsymbol{\tilde g}}}_k}} \right\|_2^2}} $

(13)

式中:$ \mathit{\boldsymbol{\tilde g}}$_k^H为对应$ \mathit{\boldsymbol{\tilde g}}$_k的共轭算子; M_{$ \mathit{\boldsymbol{\tilde q}}$k－1}^H为对应M_{$ \mathit{\boldsymbol{\tilde q}}$k－1}的共轭算子。得到$ \mathit{\boldsymbol{\tilde q}}$_k后, 在下一次求解$ \mathit{\boldsymbol{\tilde g}}$_k+1前, s_k$ \mathit{\boldsymbol{\tilde g}}$_k要恢复到$ \mathit{\boldsymbol{\tilde g}}$_k, 以保证返回到了之前的Pareto曲线点。

上述求解过程从假定已知震源子波$ \mathit{\boldsymbol{\tilde q}}$_k开始, 因此, 可用一个零相位(也可以是最小相位或混合相位)理论子波作为初始震源子波。这就引出两个问题:①使用何种理论子波及如何确定子波初始振幅？②怎样使理论子波向实际地震子波逼近？

1.2 初始震源子波的确定

实际中的地震子波是一个很复杂的问题, 因为地震子波与地层岩石性质有关, 地层岩石性质本身就是一个复杂体。为了研究方便, 对地震子波进行模拟, 这就是理论子波。目前普遍认为Ricker提出的地震子波数学模型具有广泛的代表性, 因此, Ricker子波可作为初始震源子波。Ricker子波公式如下:

$ q(t) = A(1 - 2{\pi ^2}f_{\rm{m}}^2{t^2}){{\rm{e}}^{ - {\pi ^2}f_{\rm{m}}^2{t^2}}} $

(14)

式中:t为时间; A为振幅; f_m为主频。图 1展示了主频从15Hz到45Hz, 振幅A为1.0的Ricker子波。

假定地震子波已知为q(t), 且振幅A=1.0, 由(4)式有:

$ \begin{array}{l} \mathit{\boldsymbol{M}}(\mathit{\boldsymbol{g}},\mathit{\boldsymbol{q}}): = \\ {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} \mathit{\boldsymbol{F}}_t^{\rm{H}}{\rm{BlockDia}}{{\rm{g}}_\omega }[{(w(\omega )\mathit{\boldsymbol{I}} - \mathit{\boldsymbol{P}})^{\rm{T}}} \otimes \mathit{\boldsymbol{I}}]{{\rm{F}}_t}\mathit{\boldsymbol{g}}\\ {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} \begin{array}{*{20}{l}} {\mathit{\boldsymbol{ = F}}_t^{\rm{H}}{\rm{BlockDia}}{{\rm{g}}_\omega }[{{(w(\omega )\mathit{\boldsymbol{I}})}^{\rm{T}}} \otimes \mathit{\boldsymbol{I}}]{\mathit{\boldsymbol{F}}_t}\mathit{\boldsymbol{g + }}}\\ {{\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} \mathit{\boldsymbol{F}}_t^{\rm{H}}{\rm{BlockDia}}{{\rm{g}}_\omega }[{{( - \mathit{\boldsymbol{P}})}^{\rm{T}}} \otimes \mathit{\boldsymbol{I}}]{\mathit{\boldsymbol{F}}_t}\mathit{\boldsymbol{g}}} \end{array}\\ {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} = {\mathit{\boldsymbol{M}}_q}\mathit{\boldsymbol{g}} = {\mathit{\boldsymbol{Q}}_q}\mathit{\boldsymbol{g}} + {\mathit{\boldsymbol{P}}_q}\mathit{\boldsymbol{g}} = \mathit{\boldsymbol{p}} \end{array} $

(15)

这里定义算子:

$ \left\{ {\begin{array}{*{20}{l}} {{\mathit{\boldsymbol{M}}_q} = \mathit{\boldsymbol{F}}_t^{\rm{H}}{\rm{BlockDia}}{{\rm{g}}_\omega }[{{(w(\omega )\mathit{\boldsymbol{I}} - \mathit{\boldsymbol{P}})}^{\rm{T}}} \otimes \mathit{\boldsymbol{I}}]{\mathit{\boldsymbol{F}}_t}}\\ {{\mathit{\boldsymbol{Q}}_q} = \mathit{\boldsymbol{F}}_t^{\rm{H}}{\rm{BlockDia}}{{\rm{g}}_\omega }[{{(w(\omega )\mathit{\boldsymbol{I}})}^{\rm{T}}} \otimes \mathit{\boldsymbol{I}}]{\mathit{\boldsymbol{F}}_t}}\\ {{\mathit{\boldsymbol{P}}_q} = \mathit{\boldsymbol{F}}_t^{\rm{H}}{\rm{BlockDia}}{{\rm{g}}_\omega }[{{( - \mathit{\boldsymbol{P}})}^{\rm{T}}} \otimes \mathit{\boldsymbol{I}}]{\mathit{\boldsymbol{F}}_t}} \end{array}} \right. $

(16)

假定已知一次波脉冲响应g的估计值为$ \mathit{\boldsymbol{\tilde g}}$。则(15)式可以写成:

$ \mathit{\boldsymbol{\tilde p}} = {\mathit{\boldsymbol{M}}_\mathit{\boldsymbol{q}}}\mathit{\boldsymbol{\tilde g}} = {\mathit{\boldsymbol{Q}}_q}\mathit{\boldsymbol{\tilde g}} + {\mathit{\boldsymbol{P}}_q}\mathit{\boldsymbol{\tilde g}} = {\mathit{\boldsymbol{p}}_{\rm{p}}} + {\mathit{\boldsymbol{p}}_{\rm{m}}} $

(17)

式中:$ \mathit{\boldsymbol{\tilde p}}$=M_q$ \mathit{\boldsymbol{\tilde g}}$表示对整体地震数据p的估计; p_p=Q_q$ \mathit{\boldsymbol{\tilde g}}$表示对一次波的估计; p_m=P_q$ {\tilde g}$表示对多次波的估计。

为了使子波初始振幅A与实际数据相匹配, 则要求整体地震数据的估计值$ \mathit{\boldsymbol{\tilde p}}$与实际地震记录p在能量水平上接近, 用均方根数值来衡量。设rms(·)表示均方根函数, s_{$ \tilde p $}表示比例系数, 则有:

$ {s_{\tilde p}} = \frac{{ {\rm{rms}} (\mathit{\boldsymbol{p}})}}{{ {\rm{rms}} (\mathit{\boldsymbol{\tilde p}})}} $

(18)

s_{$ {\tilde p}$}表示对$ \mathit{\boldsymbol{\tilde g}}$的调整, 即$ \mathit{\boldsymbol{\tilde g}}$_{$ \mathit{\boldsymbol{\tilde p}}$}=s_{$ {\tilde p}$}$ \mathit{\boldsymbol{\tilde g}}$。

同理, 要求多次波的估计值p_m和一次波估计值p_p, 与实际地震记录p也在能量水平上接近, 分别用s_m和s_p表示比例系数, 则有:

$ \begin{array}{*{20}{l}} {{s_{\rm{m}}} = \frac{{{\rm{rms}}(\mathit{\boldsymbol{p}})}}{{{\rm{rms}}({\mathit{\boldsymbol{p}}_{\rm{m}}})}}}\\ {{s_{\rm{p}}} = \frac{{{\rm{rms}}(\mathit{\boldsymbol{p}})}}{{{\rm{rms}}({\mathit{\boldsymbol{p}}_{\rm{p}}})}}} \end{array} $

(19)

这里, s_m表示对$ \mathit{\boldsymbol{\tilde g}}$_{$ {\tilde p}$}的调整, 即$ \mathit{\boldsymbol{\tilde g}}$_m=s_m$ \mathit{\boldsymbol{\tilde g}}$_{$ {\tilde p}$}, 则一次波的估计值p_p=Q_q$ \mathit{\boldsymbol{\tilde g}}$_m, 其中包含震源子波q(t)和做过二次调整的一次波脉冲响应$ \mathit{\boldsymbol{\tilde g}}$_m。因此, 可以用s_p调整q(t), 有:

$ \mathit{\boldsymbol{\tilde q}}(t) = {s_{\rm{p}}}\mathit{\boldsymbol{q}}(t) $

(20)

由上述过程得到A=s_p, 则初始震源子波$ \mathit{\boldsymbol{\tilde q}}$(t)在能量级别上与输入的地震记录相匹配。

1.3 自适应逼近实际地震子波

通常在EPSI算法中, 使用SPGL1方法求解, 其中子波求解采用常规的最小平方QR分解算法或匹配滤波算法。这里采用无约束非线性优化算法。定义目标函数为:

$ f(\mathit{\boldsymbol{q}}) = \left\| {\mathit{\boldsymbol{p}} - {\mathit{\boldsymbol{M}}_{{{\mathit{\boldsymbol{\tilde g}}}_{k + 1}}}}\mathit{\boldsymbol{q}}} \right\|_2^2 $

(21)

则无约束非线性优化问题为:

$ \mathop {{\rm{ argmin }}}\limits_\mathit{\boldsymbol{q}} f(\mathit{\boldsymbol{q}}) $

(22)

在该问题((22)式)中, 子波q是要求解的变量。求解该问题有多种算法, 我们采用置信域优化算法, 用能量调整后的初始震源子波$ \mathit{\boldsymbol{\tilde q}} $(t)作为优化问题的初始点, 对$ \mathit{\boldsymbol{\tilde q}} $(t)不加任何约束, 但$ \mathit{\boldsymbol{\tilde q}} $(t)的长度由理论子波确定。在迭代过程中, 震源子波的主频、振幅、相位等特性自适应地逐步逼近实际地震子波, 因此, 本文提出的方法叫做自适应稀疏反演多次波压制方法。

实现流程如下。

1) 输入地震记录p, 设定残差值σ(‖p‖₂的5%~10%)。

2) 根据地震记录的主频f_m计算振幅A=1的Ricker子波q(t)。

3) 多维互相关得到一次波脉冲响应的估计值$ \mathit{\boldsymbol{\tilde g}}$。根据(18)式和(19)式计算s_{$ {\rm{\tilde p}} $}、s_m和s_p。

4) 得到初始$ \mathit{\boldsymbol{\tilde q}}$(t)=s_pq(t)。

5) 迭代计数k=0, L₁范数约束τ₀=0。

6) 循环开始。

7) 由σ和τ_k, 通过SPGL1确定τ_k+1。

8) ${\mathit{\boldsymbol{\tilde g}}_{k + 1}} \leftarrow \mathop {{\rm{argmin}}}\limits_\mathit{\boldsymbol{g}} {\left\| {\mathit{\boldsymbol{p}} - {\mathit{\boldsymbol{M}}_{\tilde qk}}\mathit{\boldsymbol{g}}} \right\|_2}{\rm{s}}{\rm{.t}}.\;\;{\left\| \mathit{\boldsymbol{g}} \right\|_1} \le {\tau _{k + 1}} $。

9) 按照(13)式, 调整$ \mathit{\boldsymbol{\tilde g}}$_k+1:$ \mathit{\boldsymbol{\tilde g}}$_k+1=s_k+1$ \mathit{\boldsymbol{\tilde g}}$_k+1。

10) 采用置信域优化算法更新地震子波: ${\mathit{\boldsymbol{\tilde q}}_{k + 1}} \leftarrow \mathop {{\rm{argmin}}}\limits_\mathit{\boldsymbol{q}} \left\| {\mathit{\boldsymbol{p}} - {\mathit{\boldsymbol{M}}_{_{\tilde gk + 1}}}\mathit{\boldsymbol{q}}} \right\|_2^2{\rm{s}}{\rm{.t}}{\rm{.}}\;\mathit{\boldsymbol{q}} \in \mathit{\Lambda } $。

11) 重新调整为$ \mathit{\boldsymbol{\tilde g}}$_k+1:$ \mathit{\boldsymbol{\tilde g}}$_k+1=$ \mathit{\boldsymbol{\tilde g}}$_k+1/s_k+1。

12) k=k+1。

13) 直到:‖p－M($ \mathit{\boldsymbol{\tilde g}}$, $ \mathit{\boldsymbol{\tilde q}}$)‖₂≤σ或达到预定的迭代次数, 停止循环。

14) 输出:一次波估计p_p=Q_q$ \mathit{\boldsymbol{\tilde g}}$_k、多次波估计p_m=P_q$ \mathit{\boldsymbol{\tilde g}}$_k和震源子波估计$ \mathit{\boldsymbol{\tilde q}}$_k。

2 应用算例 2.1 简单岩丘模型

设计一个如图 2所示的顶部有水层的岩丘地质模型。该模型的水层深度在200m左右, 横向变化的盐丘层大致位于深度400~800m。通过时间二阶、空间四阶有限差分声波程序进行正演模拟, 使用主频为25Hz的零相位Ricker子波。震源线与接收线精确重合, 道间距为15m。

2.1.1 应用效果测试分析

抽取排列在岩丘高点附近的150道数据进行处理。图 3a是测线中部相邻两炮原始数据, 水层底部和盐丘的顶部产生的一阶及高阶多次波(图中白色箭头指示处)严重影响了一次波同相轴。图 3b是采用本文方法压制多次波后的对应单炮记录, 图中白色箭头指示处的多次波已得到有效压制, 一次波得到很好的恢复, 1.6s以下空白, 没有信息, 这与模型1700m以下没有反射层相对应。图 3c是采用EPSI方法压制多次波后的对应两炮记录, 白色箭头处的多次波也得到压制。对比图 3b和图 3c可以看出, 图 3b中的绕射波更完整, 同时最右边的箭头和最下边两个箭头处, 原来被多次波掩盖的有效波得到更好的恢复, 说明采用本文方法压制多次波时, 有效波得到更好的保护和恢复。

图 4a为原始数据的零炮检距剖面, 多次波严重影响资料品质。图 4a中存在两个产生较强多次波的界面, 一个是水层顶面, 一个是水层底面, 图中白色箭头指示了表面相关一阶和高阶多次波。水层顶面产生表面相关多次波, 水层底面产生层间多次波。图 4b为采用本文方法进行多次波压制后的零炮检距剖面, 可以看出, 白色箭头所指水层表面产生的各阶多次波得到大幅压制, 被掩盖的一次波同相轴得到了很好的恢复, 多次波压制效果良好。在0.8s左右存在由水底产生的层间多次波。在1.2s左右处, 存在由水底产生的高阶层间多次波。对于这些层间多次波, 可采用其它层间多次波压制方法进一步压制。图 4c为采用EPSI进行多次波压制后的零炮检距剖面。可以看出, 图 4c整体效果与图 4b相当, 但图 4b中的绕射波信息更完整, 同时最下面的两个白色箭头处的弱反射也完整地显现出来了。

2.1.2 初始子波误差容忍度分析

为了分析本文方法对初始子波误差的容忍度, 基于同样的数据集, 分别采用主频为20Hz的零相位、最小相位和混合相位Ricker子波作为初始输入子波, 应用本文方法进行处理和分析。

不同相位子波输入的多次波压制结果如图 5所示。图 5a为原始单炮记录(局部), 含有严重的表面相关多次波, 图中白色箭头指示了主要多次波。图 5b至图 5d分别为零相位子波、最小相位子波和混合相位Ricker子波输入时的多次波压制结果。从图 5可以看出, 表面相关多次波得到有效压制(如图中白色箭头所示), 一次波得到较好恢复, 不同相位子波输入对多次波压制效果影响较小, 表明本文方法不受输入子波类型影响, 对初始子波误差的容忍度较高。

不同相位子波输入时的子波迭代更新结果如图 6所示。图 6a为零相位子波迭代更新结果, 初始子波(k=1)经过4次迭代, 子波的主频和振幅趋于稳定, 并且向实际理论子波逼近。求取的子波(k=4)和模拟正演使用的理论子波(主频为25Hz的零相位Ricker子波)一致。图 6b为最小相位子波迭代更新结果, 初始子波(k=1)经过5次迭代, 求取的子波(k=5)和模拟正演使用的理论子波一致。图 6c为混合相位子波迭代更新结果, 初始子波(k=1)经过5次迭代, 求取的子波(k=5)和模拟正演使用的理论子波一致。

采用主频为20Hz的零相位、最小相位和混合相位Ricker子波作为初始输入子波, 应用本文方法进行处理, 最终子波都收敛到正演模拟数据的实际理论子波, 而多次波压制效果没有明显变化, 说明子波迭代更新结果受产生数据的实际子波控制, 本文方法对子波误差具有较高的容忍度。

2.2 实际数据算例

实际地震数据来自苏伊士湾海上地区(Gulf of Suez region), 浅海底产生了强烈的多次波, 每炮178道, 采样间隔为4ms, 记录长度为4s。初始子波使用20Hz零相位Ricker子波。

图 7a为包含多次波的原始单炮记录, 整个数据体被强烈的多次波覆盖。图 7b为采用本文方法进行多次波压制后的单炮记录, 可以看出, 不同阶次的多次波都得到很好的压制, 一次波得到恢复。

图 8a为原始数据叠加剖面, 多次波在剖面中普遍存在。图 8b为采用本文方法进行多次波压制后的叠加剖面, 可以看出, 海洋表面相关各阶多次波得到很好的压制。

采用本文方法迭代求解时的子波变化情况如图 9所示。可以看出, 初始子波(k=1)经过5次迭代, 子波趋于稳定收敛, 第4次迭代(k=4)和第5次迭代(k=5)的子波几乎重合。

3 结论

本文以EPSI多次波压制方法频域理论公式为基础, 将EPSI过程表示为双线性算子一次波脉冲响应和震源子波的约束优化问题, 提出了一种自适应稀疏反演多次波压制方法(AEPSI), 在迭代求解的过程中, 震源子波的主频、振幅、相位等特性, 自适应逐步逼近实际地震子波, 避免了通过输入记录的多维互相关数据来估计初始震源子波的不确定性以及求取一次脉冲响应时的人工拾取和时窗设置。模型数据和实际资料处理结果表明, 本文方法正确有效, 稳定收敛, 对子波误差具有较高的容忍度, 整个处理过程完全由数据驱动, 无需人工干预, 具有较好的应用前景。