CN105259574B - 基于一次波稀疏约束的多道预测反褶积方法 - Google Patents

Abstract

本发明属于地震勘探技术中地震信号处理领域，具体公开了一种基于一次波稀疏约束的多道预测反褶积方法。该方法首先确定多道预测反褶积中2D预测滤波器的有限支撑域和相应的数学模型，降低所求解的2D预测滤波器的系数个数，然后构建对一次波施加稀疏约束的优化问题，并采用快速迭代收缩阈值算法求解2D预测滤波器，实现对多次波的压制。相比于传统的多道预测反褶积方法需要估计滤波器系数空间的所有滤波器系数，并对一次波施加能量最小化约束来求解2D预测滤波器而言，本发明方法能减小所求解的滤波器系数个数，有效地均衡一次波的保护和多次波的压制，同时降低求解优化问题的计算复杂度。

Description

基于一次波稀疏约束的多道预测反褶积方法

技术领域

本发明属于地震勘探技术中地震信号处理领域，具体涉及一种基于一次波稀疏约束的多道预测反褶积方法。

背景技术

在海洋地震勘探中，预测反褶积用来消除水层多次波。多道预测反褶积方法比单道预测反褶积方法能更好地消除多次波(M.T.Taner,“Long period sea-floor multiplesand their suppression,”Geophysical Prospecting,vol.28,no.1,pp.30-48,Feb.1980.)。多道预测反褶积方法采用2D预测滤波器将原始数据的多道进行结合来预测多次波。为避免可能存在的一次波损伤，多道预测反褶积采用相同的2D预测滤波器来同时预测多道中的多次波。因此，可以在相互重叠的时间和空间数据窗口内利用多道预测反褶积进行多次波压制，并能有效地适应地震数据的时间和空间变化特性。多道预测反褶积方法的数学模型为：

$\begin{array}{c}{v}_{q+i,j}=u_{q+i,j}-\underset{k=-K}{\overset{K}{\Σ}}\underset{l=-L}{\overset{L}{\Σ}}{x}_{k,l}{u}_{k+i,l+j}\\ \begin{array}{cc}i=1,2,...,T_0& j=1,2,...,X_0\end{array}\end{array},---\left(1\right)$

或

v＝u-Ux， (2)

其中，表示估计一次波，表示原始数据，x＝[x_-K,-L…x_{-K, L}x_-K+1,-L…x_-K+1,L…x_K,-L…x_K,L]^T表示2D预测滤波器，U表示数据的褶积矩阵。

在方程(1)中，q为时间采样点个数，表示时间方向的预测步长；2K+1和T₀为时间采样点个数，分别表示2D预测滤波器的时间长度和数据窗口的时间长度；2L+1和X₀为道数，分别表示2D预测滤波器的空间长度和数据窗口的空间长度。

在多道预测反褶积中，2D预测滤波器的系数绝对值中的大值集中在沿同相轴倾角方向的有限范围内(J.Claerbout.(2006,Feb 28).Image estimation by example:Geophysical soundings image construction.[Online].Available:http://www.reproducibility.org/RSF/book/gee/)，称之为滤波器的有限支撑域。在滤波器的系数空间中，有限支撑域内的滤波器系数对多道预测反褶积的多次波预测起主要作用。通常，与单道预测反褶积中的1D预测滤波器相比，2D预测滤波器能够更好地消除多次波，但需要更多的计算时间。

传统多道预测反褶积对估计一次波施加能量最小化约束。另外，为避免滤波器估计的不稳定性，滤波器系数也假设满足能量最小化约束。相应的优化问题为

$\arg \underset{x}{min}\left|\right|u-Ux|{|}_2^2+\mathrm{\λ}|\left|x\right|{|}_2^2,---\left(3\right)$

其中，λ为正则化参数。方程(3)中的2D预测滤波器可以采用最小二乘算法进行求解：

x＝(U^TU+λI)^-1U^Tu。 (4)

最小二乘算法需要一次波和多次波的正交性假设。当一次波和多次波相互重叠或有强一次波同相轴存在时，最小二乘算法会产生残余多次波或造成一次波的损伤。在单道预测反褶积中，对一次波施加稀疏约束，并采用迭代重加权最小二乘算法或快速迭代收缩阈值算法求解1D预测滤波器，能有效地均衡一次波的保护和多次波的压制。其中，快速迭代收缩阈值算法采用1D收缩阈值算子促进一次波的非高斯性，相比于迭代重加权最小二乘算法，在保持计算精度的同时，能进一步提高计算效率(L.Liu and W.Lu,“A fast L1linearestimator and its application on predictive deconvolution,”IEEE Geosciencesand Remote Sensing Letters,vol.12,no.5,pp.1056-1060,May.2015.)。然而，传统的多道预测反褶积方法需要求解滤波器系数空间的所有滤波器系数，并对一次波施加能量最小化约束来求解2D预测滤波器，不能有效地均衡一次波的保护和多次波的压制。

发明内容

本发明的目的在于提出一种基于一次波稀疏约束的多道预测反褶积方法，该方法首先确定2D预测滤波器的有限支撑域，降低所求解的滤波器系数的个数，然后构建对一次波施加稀疏约束的优化问题，并采用快速迭代收缩阈值算法降低优化问题求解的计算复杂度，同时能有效地均衡一次波的保护和多次波的压制。

为了实现上述目的，本发明采用如下技术方案：

基于一次波稀疏约束的多道预测反褶积方法，包括如下步骤：

a设置变量初始值，需要设置初始值的变量包括时间方向的预测步长q，2D数据窗口的时间长度T₀，2D数据窗口的空间长度X₀，2D预测滤波器的时间长度2D预测滤波器的空间长度滤波器阈值ρ，一次波阈值s_α，阻尼因子β和迭代次数

b输入一个待处理的2D数据窗口中的数据u，根据比u的旅行时小q个时间采样点的2D数据窗口内的数据以及数据窗口的长度参数、2D滤波器的长度参数构造褶积矩阵U，并确定滤波器的有限支撑域为：

$\hat{x}=U^{T}u,\hat{x}=\left\{{\hat{x}}_{i,j}\right\},\mathrm{\Ω}=\left\{(i,j)\right|\left|{\hat{x}}_{i,j}\right|>\mathrm{\ρ}\×max\left(\right|\hat{x}\left|\right)\},$

其中，表示2D预测滤波器的粗略估计，表示向量中下标为(i,j)的元素，Ω表示滤波器的有限支撑域；

c根据2D预测滤波器的有限支撑域，确定多道预测反褶积的数学模型为：

v＝u-U_Ωx_Ω，

其中，v表示一次波，x_Ω只包含有限支撑域内的滤波器系数，U_Ω为相应的褶积矩阵；计算逆矩阵I表示单位矩阵；

d利用步骤c得到的逆矩阵对地震道集的2D数据窗口的数据u进行处理；

e判断该地震道集中所有数据窗口内的数据u是否全部处理完毕；如果否，返回步骤b；如果全部处理完毕，则首先采用2D汉宁窗将每一个2D数据窗口中估计的一次波进行加权，并融合为一个道集然后采用同样的方式将2D汉宁窗融合为一个道集最终的一次波估计结果表示为：其中，/表示逐个元素的相除操作。

在步骤d中，利用逆矩阵对地震道集的2D数据窗口的数据u进行处理的具体过程为：

d1设置迭代次数m＝0，利用步骤c得到的逆矩阵求取一次波的初始估计值：

$v^{\left(0\right)}=u-U_{\mathrm{\Ω}}\left({\stackrel{\‾}{U}}_{\mathrm{\Ω}}\right({U_{\mathrm{\Ω}}}^{T}u\left)\right);$

d2令m＝m+1,对上一步估计的一次波进行收缩阈值操作：

${\tilde{v}}^{\left(m\right)}=T_{\mathrm{\α}}\left(v^{\left(m\right)}\right),$

其中，T_α为收缩阈值算子，定义为：

T_α(v)_i,j＝(|v_i,j|-s_αC)₊sgn(v_i,j)，

i＝1,2,…,T₀；j＝1,2,…,X₀；

其中，表示估计一次波的收缩阈值结果，v^(m)表示估计的一次波，v_i,j表示向量v中下标为(i,j)的元素，v＝{v_i,j}，C＝max(|v_i,j|)，

d3对收缩阈值结果进行更新：

${\tilde{y}}^{\left(m\right)}=\left\{\begin{array}{cc}{\tilde{v}}^{\left(1\right)},& m=1\\ {\tilde{v}}^{\left(m\right)}+\left(\right(t^{\left(m\right)}-1)/t^{(m+1)})({\tilde{v}}^{\left(m\right)}-{\tilde{v}}^{(m-1)}),& m>1\end{array}\right.,$

其中，表示更新后的收缩阈值结果，序列{t^(m)}表达为t⁽¹⁾＝1，

d4求取更新后的一次波估计结果：

$v^{\left(m\right)}=u-U_{\mathrm{\Ω}}\left({\stackrel{\‾}{U}}_{\mathrm{\Ω}}\right({U_{\mathrm{\Ω}}}^T\left(u-{\tilde{y}}^{\left(m\right)}\right)\left)\right);$

d5判断迭代次数m是否达到迭代次数如没达到，返回步骤d2；如果达到，输出当前数据窗口的一次波估计结果。

本发明具有如下优点：

本发明方法首先确定多道预测反褶积中2D预测滤波器的有限支撑域和相应的数学模型，降低所求解的2D预测滤波器的系数个数，然后构建对一次波施加稀疏约束的优化问题，并采用快速迭代收缩阈值算法求解2D预测滤波器，实现对多次波的压制。相比于传统的多道预测反褶积方法需要估计滤波器系数空间的所有滤波器系数，并对一次波施加能量最小化约束来求解2D预测滤波器而言，本发明方法能减小所求解的滤波器系数个数，有效地均衡一次波的保护和多次波的压制，同时降低求解优化问题的计算复杂度。

附图说明

图1为本发明中基于一次波稀疏约束的多道预测反褶积方法的流程示意图；

图2a为偏移距等于40米的共偏移距道集图；

图2b为将图2a中每一道延迟400毫秒的共偏移距道集图；

图2c为在图2a中人为地加入一个倾斜的一次波同相轴的共偏移距道集图；

图2d为真实的一次波的共偏移距道集图；

图2e为滤波器有限支撑域的示意图；

图3a为本发明中方法估计的一次波图；

图3b为本发明中方法去除的多次波图；

图3c为不利用滤波器有限支撑域的基于快速迭代收缩阈值算法的多道预测反褶积方法估计的一次波图；

图3d为不利用滤波器有限支撑域的基于快速迭代收缩阈值算法的多道预测反褶积方法去除的多次波图；

图4a为利用滤波器有限支撑域的基于迭代重加权最小二乘算法的多道预测反褶积方法估计的一次波图；

图4b为利用滤波器有限支撑域的基于迭代重加权最小二乘算法的多道预测反褶积方法去除的多次波图；

图4c为不利用滤波器有限支撑域的基于迭代重加权最小二乘算法的多道预测反褶积方法估计的一次波图；

图4d为不利用滤波器有限支撑域的基于迭代重加权最小二乘算法的多道预测反褶积方法去除的多次波图；

图5a为传统的基于最小二乘算法的多道预测反褶积方法估计的一次波图；

图5b为传统的基于最小二乘算法的多道预测反褶积方法去除的多次波图；

图5c为传统的基于快速迭代收缩阈值算法的单道预测反褶积方法估计的一次波图；

图5d为传统的基于快速迭代收缩阈值算法的单道预测反褶积方法去除的多次波图。

具体实施方式

本发明的基本思想是：

逐个2D数据窗口地进行多次波压制，首先确定2D预测滤波器的有限支撑域，然后构建相应的褶积矩阵和数学模型，并构建对一次波施加稀疏约束的优化问题：

$\arg \underset{x_{\mathrm{\Ω}}}{min}\left|\right|u-U_{\mathrm{\Ω}}{x}_{\mathrm{\Ω}}|{|}_1,$

其中，u为原始数据，x_Ω只包含有限支撑域内的滤波器系数，U_Ω为相应的褶积矩阵。求解上式中的优化问题对2D预测滤波器进行估计，实现对该2D数据窗口中一次波的估计，最后将所有2D数据窗口内的一次波估计结果进行合并，得到该道集的一次波估计结果。

下面结合附图以及具体实施方式对本发明作进一步详细说明：

结合图1所示，基于一次波稀疏约束的多道预测反褶积方法，包括如下步骤：

a设置变量初始值，需要设置初始值的变量包括时间方向的预测步长q，2D数据窗口的时间长度T₀，2D数据窗口的空间长度X₀，2D预测滤波器的时间长度2D预测滤波器的空间长度滤波器阈值ρ，一次波阈值s_α，阻尼因子β和迭代次数

b输入一个待处理的2D数据窗口中的数据u，根据比u的旅行时小q个时间采样点的2D数据窗口内的数据以及数据窗口的长度参数、2D滤波器的长度参数构造褶积矩阵U，并确定滤波器的有限支撑域为：

$\hat{x}=U^{T}u,\hat{x}=\left\{{\hat{x}}_{i,j}\right\},\mathrm{\Ω}=\left\{(i,j)\right|\left|{\hat{x}}_{i,j}\right|>\mathrm{\ρ}\×max\left(\right|\hat{x}\left|\right)\},$

其中，表示2D预测滤波器的粗略估计，表示向量中下标为(i,j)的元素，Ω表示滤波器的有限支撑域；

c根据滤波器的有限支撑域来确定多道预测反褶积的数学模型为：

v＝u-U_Ωx_Ω，

其中，v表示一次波，x_Ω只包含有限支撑域内的滤波器系数，U_Ω为相应的褶积矩阵；计算逆矩阵I表示单位矩阵；

本发明中方法构建对一次波施加稀疏约束的优化问题：

$\arg \underset{x_{\mathrm{\Ω}}}{min}\left|\right|u-U_{\mathrm{\Ω}}{x}_{\mathrm{\Ω}}|{|}_1;$

d利用步骤c得到的逆矩阵对地震道集的2D数据窗口的数据u进行处理；

本发明方法采用快速迭代收缩阈值算法求解上述公式中的2D预测滤波器，实现多次波压制，即：在每一步迭代，快速迭代收缩阈值算法采用收缩阈值算子促进一次波的稀疏性，得到一次波的收缩阈值结果，并将当前步和前一步的收缩阈值结果进行线性组合来对当前步的收缩阈值结果进行更新，然后采用更新后的收缩阈值结果来估计一次波；

其具体处理过程为：

d1设置迭代次数m＝0，利用步骤c得到的逆矩阵求取一次波的初始估计值：

$v^{\left(0\right)}=u-U_{\mathrm{\Ω}}\left({\stackrel{\‾}{U}}_{\mathrm{\Ω}}\right({U_{\mathrm{\Ω}}}^{T}u\left)\right);$

d2令m＝m+1,对上一步估计的一次波进行收缩阈值操作：

${\tilde{v}}^{\left(m\right)}=T_{\mathrm{\α}}\left(v^{\left(m\right)}\right),$

其中，T_α为收缩阈值算子，定义为：

T_α(v)_i,j＝(|v_i,j|-s_αC)₊sgn(v_i,j)，

i＝1,2,…,T₀；j＝1,2,…,X₀；

其中，表示估计一次波的收缩阈值结果，v^(m)表示估计的一次波，v_i,j表示向量v中下标为(i,j)的元素，v＝{v_i,j}，C＝max(|v_i,j|)，

d3对收缩阈值结果进行更新：

${\tilde{y}}^{\left(m\right)}=\{\begin{array}{cc}{\tilde{v}}^{\left(1\right)},& m=1\\ {\tilde{v}}^{\left(m\right)}+\left(\right(t^{\left(m\right)}-1)/t^{(m+1)})({\tilde{v}}^{\left(m\right)}-{\tilde{v}}^{(m-1)}),& m>1\end{array},$

其中，表示更新后的收缩阈值结果，序列{t^(m)}表达为t⁽¹⁾＝1，

d4求取更新后的一次波估计结果：

$v^{\left(m\right)}=u-U_{\mathrm{\Ω}}\left({\stackrel{\‾}{U}}_{\mathrm{\Ω}}\right({U_{\mathrm{\Ω}}}^T\left(u-{\tilde{y}}^{\left(m\right)}\right)\left)\right);$

d5判断迭代次数m是否达到迭代次数如没达到，返回步骤d2；如果达到，输出当前数据窗口的一次波估计结果。

e判断该地震道集中所有数据窗口内的数据是否全部处理完毕；如果否，返回步骤b；如果全部处理完毕，则首先采用2D汉宁窗将每一个2D数据窗口中估计的一次波进行加权，并融合为一个道集然后采用同样的方式将2D汉宁窗融合为一个道集最终的一次波估计结果表示为：其中，/表示逐个元素的相除操作。

在仿真实验中，利用2D模型数据验证本发明的有效性：

图2a为偏移距等于40米的共偏移距道集图。根据图2a中每一道的互相关值，时间预测长度被选为400毫秒。将图2a中每一道的旅行时延迟400毫秒可以得到图2b中的道集。为进一步验证本发明能有效地去除残余多次波，一个倾斜的一次波同相轴人为地加入到原始数据中，相应的道集显示在图2c中。真实的一次波道集显示在图2d中。预测反褶积通过校正图2b中道集和图2c中道集的时间、空间和振幅差异来求解预测滤波器，实现多次波压制。图3a和图3b为本发明方法得到的一次波估计结果和所去除的多次波图。图2e中的白色区域表示滤波器的有限支撑域，该滤波器对应于图2c中白色方框标明的2D数据窗口。在本发明中，图2e中有限支撑域内的滤波器系数得到估计。图3c和图3d为不利用滤波器有限支撑域的基于快速迭代收缩阈值算法的多道预测反褶积方法得到的一次波估计结果和所去除的多次波图。图4a和图4b为利用滤波器有限支撑域的基于迭代重加权最小二乘算法的多道预测反褶积方法得到的一次波估计结果和所去除的多次波图。图4c和图4d为不利用滤波器有限支撑域的基于迭代重加权最小二乘算法的多道预测反褶积方法得到的一次波估计结果和所去除的多次波图。对于图3a、图3c、图4a和图4c中的一次波估计结果，信噪比分别为45.77、45.19、45.32和44.65，计算时间分别为234毫秒、292毫秒、487毫秒和659毫秒。对于基于快速迭代收缩阈值算法的多道预测反褶积，本发明中确定预测滤波器有限支撑域的策略在保持计算精度的同时，降低了19.9％的计算时间。对于基于迭代重加权最小二乘算法的多道预测反褶积，本发明中确定预测滤波器有限支撑域的策略在保持计算精度的同时，降低了26.1％的计算时间。另外，所引入的快速迭代收缩阈值算法比传统的迭代重加权最小二乘算法，在保持计算精度的同时，能有效地降低计算量。传统的基于最小二乘算法的多道预测反褶积方法估计的一次波显示在图5a中，信噪比为36.60，去除的多次波显示在图5b中。图5a中的白色箭头表明最小二乘算法在有强一次波存在的地方会造成残余多次波。图5c和图5d为传统的基于快速迭代收缩阈值算法的单道预测反褶积方法估计的一次波和去除的多次波图。对于图5c中估计的一次波，信噪比为17.21。在这个例子中，单道预测反褶积会造成一次波损伤。相对于传统的基于快速迭代收缩阈值算法的单道预测反褶积方法和基于最小二乘算法的多道预测反褶积方法，本发明中方法能更好地均衡一次波的保护和多次波的压制。

其中，图2a至2d、图3a至3d、图4a至4d及图5a至5d中横坐标Trace Number表示道号，纵坐标Time表示时间，单位为毫秒(ms)。图2e中横坐标Trace Number表示道号，纵坐标Time Sample Number表示时间采样点个数。

当然，以上说明仅仅为本发明的较佳实施例，本发明并不限于列举上述实施例，应当说明的是，任何熟悉本领域的技术人员在本说明书的教导下，所做出的所有等同替代、明显变形形式，均落在本说明书的实质范围之内，理应受到本发明的保护。

Claims (2)

1.基于一次波稀疏约束的多道预测反褶积方法，其特征在于，包括如下步骤：

a设置变量初始值，需要设置初始值的变量包括时间方向的预测步长q，2D数据窗口的时间长度T₀，2D数据窗口的空间长度X₀，2D预测滤波器的时间长度2D预测滤波器的空间长度滤波器阈值ρ，一次波阈值s_α，阻尼因子β和迭代次数

b输入一个待处理的2D数据窗口中的数据u，根据比u的旅行时小q个时间采样点的2D数据窗口内的数据以及数据窗口的长度参数、2D滤波器的长度参数构造褶积矩阵U，并确定滤波器的有限支撑域为：

$\hat{x}=U^{T}u,\hat{x}=\left\{{\hat{x}}_{i,j}\right\},\mathrm{\Ω}=\left\{(i,j)\right|\left|{\hat{x}}_{i,j}\right|>\mathrm{\ρ}\×max\left(\right|\hat{x}\left|\right)\},$

其中，表示2D预测滤波器的粗略估计，表示向量中下标为(i,j)的元素，Ω表示滤波器的有限支撑域；

c根据2D预测滤波器的有限支撑域，确定多道预测反褶积的数学模型为：

v＝u-U_Ωx_Ω，

其中，v表示一次波，x_Ω只包含有限支撑域内的滤波器系数，U_Ω为相应的褶积矩阵；计算逆矩阵I表示单位矩阵；

d利用步骤c得到的逆矩阵对地震道集的2D数据窗口的数据u进行处理；

e判断该地震道集中所有数据窗口内的数据u是否全部处理完毕；如果否，返回步骤b；如果全部处理完毕，则首先采用2D汉宁窗将每一个2D数据窗口中估计的一次波进行加权，并融合为一个道集然后采用同样的方式将2D汉宁窗融合为一个道集最终的一次波估计结果表示为：其中，/表示逐个元素的相除操作。

2.根据权利要求1所述的基于一次波稀疏约束的多道预测反褶积方法，其特征在于，步骤d中，利用逆矩阵对地震道集的2D数据窗口的数据u进行处理的具体过程为：

d1设置迭代次数m＝0，利用步骤c得到的逆矩阵求取一次波的初始估计值：

$v^{\left(0\right)}=u-U_{\mathrm{\Ω}}\left({\stackrel{\‾}{U}}_{\mathrm{\Ω}}\right({U_{\mathrm{\Ω}}}^{T}u\left)\right);$

d2令m＝m+1,对上一步估计的一次波进行收缩阈值操作：

${\tilde{v}}^{\left(m\right)}=T_{\mathrm{\α}}\left(v^{\left(m\right)}\right),$

其中，T_α为收缩阈值算子，定义为：

T_α(v)_i,j＝(|v_i,j|-s_αC)₊sgn(v_i,j)，

i＝1,2,…,T₀；j＝1,2,…,X₀；

其中，表示估计一次波的收缩阈值结果，v^(m)表示估计的一次波，v_i,j表示向量v中下标为(i,j)的元素，v＝{v_i,j}，C＝max(|v_i,j|)，

d3对收缩阈值结果进行更新：

${\tilde{y}}^{\left(m\right)}=\left\{\begin{array}{cc}{\tilde{v}}^{\left(1\right)},& m=1\\ {\tilde{v}}^{\left(m\right)}+\left(\right(t^{\left(m\right)}-1)/t^{(m+1)})({\tilde{v}}^{\left(m\right)}-{\tilde{v}}^{(m-1)}),& m>1\end{array}\right.,$

其中，表示更新后的收缩阈值结果，序列{t^(m)}表达为t⁽¹⁾＝1，d4求取更新后的一次波估计结果：

$v^{\left(m\right)}=u-U_{\mathrm{\Ω}}\left({\stackrel{\‾}{U}}_{\mathrm{\Ω}}\right({U_{\mathrm{\Ω}}}^T\left(u-{\tilde{y}}^{\left(m\right)}\right)\left)\right);$

d5判断迭代次数m是否达到迭代次数如没达到，返回步骤d2；如果达到，输出当前数据窗口的一次波估计结果。