CN103558633A - 一种基于稀疏反演的层间多次波压制方法 - Google Patents

Abstract

本发明涉及一种基于稀疏反演的层间多次波压制方法，其包括以下步骤：1）在叠加剖面上选取介于强反射界面之间的基准面，在炮集上分离来自基准面上地层的波场数据δP和下地层的波场数据

2）拾取δP初至一次波振幅峰值对应的脉冲估计基准面以上地层脉冲响应的初始值，并通过匹配滤波估计震源子波初始值，并设定反演脉冲响应与震源子波的迭代次数；3）迭代更新基准面以上地层的脉冲响应、基准面以下地层的脉冲响应以及震源子波；4）匹配层间多次波预测结果与原始地震数据得到振幅校正因子，利用振幅校正因子更新脉冲响应和震源子波；5）以步骤4）的参数更新结果作为新的初始值，继续重复步骤3）～4）直至达到给定迭代次数，停止迭代；6）从地震数据P中减去预测得到的层间多次波M得到一次波P₀。

Description

一种基于稀疏反演的层间多次波压制方法

技术领域

本发明涉及地震勘探领域，特别是关于一种适用于地震勘探数据中层间多次波的有效衰减与反射信号的保真恢复的基于稀疏反演的层间多次波压制方法。

背景技术

层间多次波严重干扰目标层有效信号，增加储层解释和反演的难度。由于形态复杂，与有效信号速度差异小等特点，层间多次波衰减一直是勘探地球物理研究的热点与难点。

目前，工业界采用的波动方程预测相减法通过先预测出层间多次波模型，然后从地震数据中自适应减去的步骤压制层间多次波。由于层间多次波往往与有效信号交叉甚至重合，违背了常规自适应相减方法的前提假设，因此在实际应用中，波动方程预测相减法会损伤有效信号，影响地震资料保真度，从而进一步影响后续的处理、反演与解释。

发明内容

针对上述问题，本发明的目的是提供一种能够检验一次波与层间多次波的物理一致性，最大程度地拟合观测地震数据的基于稀疏反演的层间多次波压制方法。

为实现上述目的，本发明采取以下技术方案：一种基于稀疏反演的层间多次波压制方法，其包括以下步骤：1）在叠加剖面上选取介于强反射界面之间的基准面，在炮集上分离来自基准面上地层的波场数据δP和下地层的波场数据

2）拾取δP初至一次波振幅峰值对应的脉冲估计基准面以上地层脉冲响应的初始值，并通过匹配滤波估计震源子波初始值，并设定反演脉冲响应与震源子波的迭代次数；3）迭代更新基准面以上地层的脉冲响应、基准面以下地层的脉冲响应以及震源子波；4）匹配层间多次波预测结果与原始地震数据得到振幅校正因子，利用振幅校正因子更新脉冲响应和震源子波；5）以步骤4）的参数更新结果作为新的初始值，继续重复步骤3）～4）直至达到给定迭代次数，停止迭代；6）从地震数据P中减去预测得到的层间多次波M得到一次波P₀。

所述步骤3）迭代更新基准面以上地层的脉冲响应、基准面以下地层的脉冲响应以及震源子波的具体计算过程为：3.1）更新基准面以上地层的脉冲响应δX：利用最速下降法得到δX在i+1次迭代时的最速下降方向：

$\mathrm{\Δ\δX}={{\mathrm{VS}}_i}^H-{\stackrel{\‾}{P}V}^H{\stackrel{\‾}{X}}_i$

变换ΔδX到时间域，拾取每道振幅峰值对应的脉冲后变换至频率域，得到脉冲响应的更新方向ΔδX，乘以不依赖于频率的步长α，更新脉冲响应与拟合残差：

δX_i+1＝δX_i+αΔδX

V＝V+αΔV₁

其中，

$\mathrm{\α}=-\frac{\underset{\mathrm{\ω}}{\mathrm{\Σ}}\mathrm{vec}{\left(\mathrm{\Δ}{V}_1\right)}^H\mathrm{vec}\left(V\right)}{\underset{\mathrm{\ω}}{\mathrm{\Σ}}\mathrm{vec}{\left(\mathrm{\Δ}{V}_1\right)}^H\mathrm{vec}\left(\mathrm{\Δ}{V}_1\right)}$

$\mathrm{\Δ}{V}_1=-\mathrm{\Δ\δ}\underset{\‾}{X}{S}_i+{\stackrel{\‾}{X}}_i\mathrm{\Δ\δ}{\underset{\‾}{X}}^H\stackrel{\‾}{P};$

3.2）更新基准面以下地层的脉冲响应利用最速下降法计算

在i+1次迭代时的最速下降方向：

$\mathrm{\Δ}\stackrel{\‾}{X}={{\mathrm{VS}}_i}^H-V{\stackrel{\‾}{P}}^H\mathrm{\δ}{X}_{i+1}$

同样的，变换

到时间域，拾取每道最大振幅值对应的脉冲后变回频率域得到

的更新方向乘以不依赖于频率的更新步长β，更新脉冲响应与拟合残差：

${\stackrel{\‾}{X}}_{i+1}={\stackrel{\‾}{X}}_i+\mathrm{\β\Δ}\underset{\‾}{\overset{\‾}{X}}$

V＝V+βΔV₂

其中，

$\mathrm{\β}=-\frac{\underset{\mathrm{\ω}}{\mathrm{\Σ}}\mathrm{vec}{\left(\mathrm{\Δ}{V}_2\right)}^H\mathrm{vec}\left(V\right)}{\underset{\mathrm{\ω}}{\mathrm{\Σ}}\mathrm{vec}{\left(\mathrm{\Δ}{V}_2\right)}^H\mathrm{vec}\left(\mathrm{\Δ}{V}_2\right)}$

$\mathrm{\Δ}{V}_2=-\mathrm{\Δ}\underset{\‾}{\overset{\‾}{X}}{S}_i+\mathrm{\Δ}\underset{\‾}{\overset{\‾}{X}}\mathrm{\δ}{X_{i+1}}^H\stackrel{\‾}{P};$

3.3）匹配脉冲响应

与

更新震源子波。

所述步骤4）振幅校正因子的计算过程为：

$\mathrm{\γ}=1-\frac{\underset{\mathrm{\ω}}{\mathrm{\Σ}}\mathrm{vec}{\left(\mathrm{\Δ}{V}_3\right)}^H\mathrm{vec}\left(V\right)}{\underset{\mathrm{\ω}}{\mathrm{\Σ}}\mathrm{vec}{\left(\mathrm{\Δ}{V}_3\right)}^H\mathrm{vec}\left(\mathrm{\Δ}{V}_3\right)}$

式中， $\mathrm{\Δ}{V}_3={\stackrel{\‾}{X}}_{i+1}\mathrm{\δ}{X_{i+1}}^H\stackrel{\‾}{P}.$

所述步骤4）利用振幅校正因子更新脉冲响应和震源子波，具体为：

$\begin{array}{c}{S}_{i+1}=\frac{1}{\sqrt{\mathrm{\γ}}}{S}_{i+1}\\ \mathrm{\δ}{X}_{i+1}=\sqrt{\mathrm{\γ}}\mathrm{\δ}{X}_{i+1}.\\ {\stackrel{\‾}{X}}_{i+1}=\sqrt{\mathrm{\γ}}{\stackrel{\‾}{X}}_{i+1}\end{array}$

所述步骤6）中预测得到的层间多次波M的计算公式为：

本发明由于采取以上技术方案，其具有以下优点：1、本发明是在波动方程预测相减法基础上提出的一种全波形反演方法，通过从地震数据中直接减去预测得到的层间多次波估计一次波，避免了自适应相减，更好地保护有效信号。2、本发明在脉冲响应的更新过程中利用层间多次波估计一次波，能够检验一次波与层间多次波的物理一致性，最大程度地拟合观测地震数据。本发明可以广泛应用于地震勘探中。

附图说明

图1是本发明基于稀疏反演的层间多次波压制方法流程示意图；

图2是本发明的叠前炮集的同相轴分离示意图，其中，图2（a）是原始炮集示意图，图2（b）是基准面以上地层的反射波场示意图，图2（c）是基准面以下地层的反射波场示意图；

图3是本发明炮集上层间多次波压制效果示意图，其中，图3（a）是原始炮集示意图，图3（b）是压制层间多次波后的示意图，图3（c）是压制掉的层间多次波；

图4是本发明叠加剖面上层间多次波压制效果示意图，其中，图4（a）是原始叠加剖面示意图，图4（b）是压制层间多次波后的示意图，图4（c）是压制掉的层间多次波示意图。

具体实施方式

下面结合附图和实施例对本发明进行详细的描述。

地震数据可以排列成以炮点位置，检波点位置和旅行时为坐标的三维数据体，沿时间方向作傅里叶变换，可以将三维数据体变换到频率域，其中单频切片对应的矩阵称为数据矩阵，本发明中采用粗体符号P表示。

如图1所示，本发明的基于稀疏反演的层间多次波压制方法，包括以下步骤：

1、在叠加剖面上选取介于强反射界面之间的基准面，在炮集上分离来自基准面上地层的波场数据δP和下地层的波场数据

实际地下介质中，能够产生明显层间多次波的强反射界面数量有限，因此可以在叠加剖面上选取基准面将强反射界面划分到不同的反射层，从而以反射层为基本物理单元预测层间多次波。利用基准面的自激自收时间和上覆地层的动校正速度，计算叠前炮集中基准面对应的时间分割线，分离出预测层间多次波所需要的分别来自基准面上地层的波场数据δP和下地层的波场数据

如图2所示，沿着计算得到的白色虚线可以将如图2（a）所示的原始炮集分割成上下两部分，得到如图2（b）与如图2（c）所示的分别来自基准面上、下地层界面的反射。

2、拾取δP初至一次波振幅峰值对应的脉冲估计基准面以上地层脉冲响应的初始值，并通过匹配滤波估计震源子波初始值，并设定反演脉冲响应与震源子波的迭代次数。

地震数据可以通过脉冲响应与震源子波表征为一次波与多次波之和（公式1），因此脉冲响应与震源子波可以作为未知参数通过全波形反演方法进行迭代估计。为了驱动该算法，必须给定至少一个未知参数的初值。由于反射波的初至能量较强，通过拾取δP初至同相轴振幅峰值对应的脉冲，可以得到基准面以上地层脉冲响应的初始估计值δX_i＝0，匹配该脉冲响应初值δX_i＝0与地震数据δP可以得到震源子波的初值S_i＝0，而基准面以下地层波场

对应的脉冲响应初值

可以直接设为0。

反演脉冲相应与震源子波的迭代次数由地震记录长度确定，一般可以选择为地震记录长度的10倍，即地震记录长度为2s则迭代次数设定为20次。

3、迭代更新基准面以上地层的脉冲响应、基准面以下地层的脉冲响应以及震源子波。

脉冲响应与震源子波的更新是整个发明的核心部分。压制自由表面多次波后的地震记录P可以利用数据矩阵乘积表示为一次波与层间多次波之和：

$P=\mathrm{\δXS}+\stackrel{\‾}{X}S-\stackrel{\‾}{X}\mathrm{\δ}{X}^H\stackrel{\‾}{P}---\left(1\right)$

式中，δXS表示基准面以上地层的反射波场，

表示基准面以下地层的反射波场，表示与基准面相关的层间多次波，其中δX^H表示δX的共轭转置。

第i次迭代后的拟合残差可以表示为：

$V=P-\mathrm{\δ}{X}_i{S}_i-{\stackrel{\‾}{X}}_i{S}_i+{\stackrel{\‾}{X}}_i\mathrm{\δ}{X_i}^H\stackrel{\‾}{P}---\left(2\right)$

3.1）更新基准面以上地层的脉冲响应δX。

利用最速下降法可以得到δX在i+1次迭代时的最速下降方向：

$\mathrm{\Δ\δX}={{\mathrm{VS}}_i}^H-{\stackrel{\‾}{P}V}^H{\stackrel{\‾}{X}}_i---\left(3\right)$

变换ΔδX到时间域，拾取每道振幅峰值对应的脉冲后变换至频率域，得到脉冲响应的更新方向ΔδX，乘以不依赖于频率的步长α，更新脉冲响应与拟合残差：

δX_i+1＝δX_i+αΔδX   (4)

V＝V+αΔV₁   (5)

其中，

$\mathrm{\α}=-\frac{\underset{\mathrm{\ω}}{\mathrm{\Σ}}\mathrm{vec}{\left(\mathrm{\Δ}{V}_1\right)}^H\mathrm{vec}\left(V\right)}{\underset{\mathrm{\ω}}{\mathrm{\Σ}}\mathrm{vec}{\left(\mathrm{\Δ}{V}_1\right)}^H\mathrm{vec}\left(\mathrm{\Δ}{V}_1\right)}---\left(6\right)$

$\mathrm{\Δ}{V}_1=-\mathrm{\Δ\δ}\underset{\‾}{X}{S}_i+{\stackrel{\‾}{X}}_i\mathrm{\Δ\δ}{\underset{\‾}{X}}^H\stackrel{\‾}{P};---\left(7\right)$

3.2）更新基准面以下地层的脉冲响应

利用最速下降法可以计算

在i+1次迭代时的最速下降方向：

$\mathrm{\Δ}\stackrel{\‾}{X}={{\mathrm{VS}}_i}^H-V{\stackrel{\‾}{P}}^H\mathrm{\δ}{X}_{i+1}---\left(8\right)$

同样的，变换

到时间域，拾取每道最大振幅值对应的脉冲后变回频率域得到

的更新方向

乘以不依赖于频率的更新步长β，更新脉冲响应与拟合残差：

${\stackrel{\‾}{X}}_{i+1}={\stackrel{\‾}{X}}_i+\mathrm{\β\Δ}\underset{\‾}{\overset{\‾}{X}}---\left(9\right)$

V＝V+βΔV₂   (10)

其中，

$\mathrm{\β}=-\frac{\underset{\mathrm{\ω}}{\mathrm{\Σ}}\mathrm{vec}{\left(\mathrm{\Δ}{V}_2\right)}^H\mathrm{vec}\left(V\right)}{\underset{\mathrm{\ω}}{\mathrm{\Σ}}\mathrm{vec}{\left(\mathrm{\Δ}{V}_2\right)}^H\mathrm{vec}\left(\mathrm{\Δ}{V}_2\right)}---\left(11\right)$

$\mathrm{\Δ}{V}_2=-\mathrm{\Δ}\underset{\‾}{\overset{\‾}{X}}{S}_i+\mathrm{\Δ}\underset{\‾}{\overset{\‾}{X}}\mathrm{\δ}{X_{i+1}}^H\stackrel{\‾}{P};---\left(12\right)$

3.3）匹配脉冲响应

与更新震源子波

震源子波S_i+1可以通过在时间域计算具有一定长度（100～200毫秒）的滤波器匹配脉冲响应

获得，其求取过程与维纳滤波完全一致，即实现脉冲响应与地震数据在最小二乘意义下的最佳匹配，其为现有技术，在此不再赘述。

4、匹配层间多次波预测结果与原始地震数据得到振幅校正因子，利用振幅校正因子更新脉冲响应和震源子波。

根据步骤3更新的三个参数（基准面以上地层的脉冲响应、基准面以下地层的脉冲响应以及震源子波），求取振幅校正因子γ，使得i+1次迭代得到的层间多次波与原始地震数据中的层间多次波振幅一致，并根据求取的振幅校正因子，进一步更新上述三个未知参数：

$\mathrm{\Δ}{V}_3={\stackrel{\‾}{X}}_{i+1}\mathrm{\δ}{X_{i+1}}^H\stackrel{\‾}{P}---\left(13\right)$

$\mathrm{\γ}=1-\frac{\underset{\mathrm{\ω}}{\mathrm{\Σ}}\mathrm{vec}{\left(\mathrm{\Δ}{V}_3\right)}^H\mathrm{vec}\left(V\right)}{\underset{\mathrm{\ω}}{\mathrm{\Σ}}\mathrm{vec}{\left(\mathrm{\Δ}{V}_3\right)}^H\mathrm{vec}\left(\mathrm{\Δ}{V}_3\right)}---\left(14\right)$

$\begin{array}{c}{S}_{i+1}=\frac{1}{\sqrt{\mathrm{\γ}}}{S}_{i+1}\\ \mathrm{\δ}{X}_{i+1}=\sqrt{\mathrm{\γ}}\mathrm{\δ}{X}_{i+1}\\ {\stackrel{\‾}{X}}_{i+1}=\sqrt{\mathrm{\γ}}{\stackrel{\‾}{X}}_{i+1}\end{array}---\left(15\right)$

5、以步骤4的参数更新结果作为新的初始值，继续重复步骤3～4直至达到给定迭代次数，停止迭代。

6、从地震数据P中减去预测得到的层间多次波M得到一次波P₀。

利用地震数据与脉冲响应的数据矩阵乘法（公式16）估计层间多次波M，并从地震记录P中直接减掉M得到压制层间多次波后的地震记录P₀。

$M=-\stackrel{\‾}{X}\mathrm{\δ}{X}^H\stackrel{\‾}{P}---\left(16\right)$

P₀＝P-M   (16)

如图3所示的本发明在叠前炮集上的层间多次波压制效果，与图3（a）所示的原始炮集相比，图3（b）所示的一次波剖面中，层间多次波得到了非常理想的压制，深层一次反射得到了有效地恢复，而且衰减的层间多次波剖面中没有发现有效信号的损伤，展现了广义EPSI优秀的保幅能力，其中，图3（c）是压制掉的层间多次波。

如图4所示的本发明在叠加剖面上的层间多次波压制效果示意图，其中，图4（a）是原始叠加剖面示意图，图4（b）是压制层间多次波后示意图，图4（c）是压制掉的层间多次波。综上可以看出，广义EPSI在压制多次波的同时，能够保真地恢复有效信号。

上述各实施例仅用于说明本发明，其中方法各实施步骤等都是可以有所变化的，凡是在本发明技术方案的基础上进行的等同变换和改进，均不应排除在本发明的保护范围之外。

Claims (6)

1.一种基于稀疏反演的层间多次波压制方法，其包括以下步骤：

1）在叠加剖面上选取介于强反射界面之间的基准面，在炮集上分离来自基准面上地层的波场数据δP和下地层的波场数据

2）拾取δP初至一次波振幅峰值对应的脉冲估计基准面以上地层脉冲响应的初始值，并通过匹配滤波估计震源子波初始值，并设定反演脉冲响应与震源子波的迭代次数；

3）迭代更新基准面以上地层的脉冲响应、基准面以下地层的脉冲响应以及震源子波；

4）匹配层间多次波预测结果与原始地震数据得到振幅校正因子，利用振幅校正因子更新脉冲响应和震源子波；

5）以步骤4）的参数更新结果作为新的初始值，继续重复步骤3）～4）直至达到给定迭代次数，停止迭代；

6）从地震数据P中减去预测得到的层间多次波M得到一次波P₀。

2.如权利要求1所述的一种基于稀疏反演的层间多次波压制方法，其特征在于：所述步骤3）迭代更新基准面以上地层的脉冲响应、基准面以下地层的脉冲响应以及震源子波的具体计算过程为：

3.1）更新基准面以上地层的脉冲响应δX：

利用最速下降法得到δX在i+1次迭代时的最速下降方向：

$\mathrm{\Δ\δX}={{\mathrm{VS}}_i}^H-{\stackrel{\‾}{P}V}^H{\stackrel{\‾}{X}}_i$

变换ΔδX到时间域，拾取每道振幅峰值对应的脉冲后变换至频率域，得到脉冲响应的更新方向ΔδX，乘以不依赖于频率的步长α，更新脉冲响应与拟合残差：

δX_i+1＝δX_i+αΔδX

V＝V+αΔV₁

其中，

$\mathrm{\α}=-\frac{\underset{\mathrm{\ω}}{\mathrm{\Σ}}\mathrm{vec}{\left(\mathrm{\Δ}{V}_1\right)}^H\mathrm{vec}\left(V\right)}{\underset{\mathrm{\ω}}{\mathrm{\Σ}}\mathrm{vec}{\left(\mathrm{\Δ}{V}_1\right)}^H\mathrm{vec}\left(\mathrm{\Δ}{V}_1\right)}$

$\mathrm{\Δ}{V}_1=-\mathrm{\Δ\δ}\underset{\‾}{X}{S}_i+{\stackrel{\‾}{X}}_i\mathrm{\Δ\δ}{\underset{\‾}{X}}^H\stackrel{\‾}{P};$

3.2）更新基准面以下地层的脉冲响应

利用最速下降法计算

在i+1次迭代时的最速下降方向：

$\mathrm{\Δ}\stackrel{\‾}{X}={{\mathrm{VS}}_i}^H-V{\stackrel{\‾}{P}}^H\mathrm{\δ}{X}_{i+1}$

同样的，变换

到时间域，拾取每道最大振幅值对应的脉冲后变回频率域得到

的更新方向

乘以不依赖于频率的更新步长β，更新脉冲响应与拟合残差：

${\stackrel{\‾}{X}}_{i+1}={\stackrel{\‾}{X}}_i+\mathrm{\β\Δ}\underset{\‾}{\overset{\‾}{X}}$

V＝V+βΔV₂

其中，

$\mathrm{\β}=-\frac{\underset{\mathrm{\ω}}{\mathrm{\Σ}}\mathrm{vec}{\left(\mathrm{\Δ}{V}_2\right)}^H\mathrm{vec}\left(V\right)}{\underset{\mathrm{\ω}}{\mathrm{\Σ}}\mathrm{vec}{\left(\mathrm{\Δ}{V}_2\right)}^H\mathrm{vec}\left(\mathrm{\Δ}{V}_2\right)}$

$\mathrm{\Δ}{V}_2=-\mathrm{\Δ}\underset{\‾}{\overset{\‾}{X}}{S}_i+\mathrm{\Δ}\underset{\‾}{\overset{\‾}{X}}\mathrm{\δ}{X_{i+1}}^H\stackrel{\‾}{P};$

3.3）匹配脉冲响应

与

更新震源子波。

3.如权利要求2所述的一种基于稀疏反演的层间多次波压制方法，其特征在于：所述步骤4）振幅校正因子的计算过程为：

$\mathrm{\γ}=1-\frac{\underset{\mathrm{\ω}}{\mathrm{\Σ}}\mathrm{vec}{\left(\mathrm{\Δ}{V}_3\right)}^H\mathrm{vec}\left(V\right)}{\underset{\mathrm{\ω}}{\mathrm{\Σ}}\mathrm{vec}{\left(\mathrm{\Δ}{V}_3\right)}^H\mathrm{vec}\left(\mathrm{\Δ}{V}_3\right)}$

式中， $\mathrm{\Δ}{V}_3={\stackrel{\‾}{X}}_{i+1}\mathrm{\δ}{X_{i+1}}^H\stackrel{\‾}{P}.$

4.如权利要求2所述的一种基于稀疏反演的层间多次波压制方法，其特征在于：所述步骤4）利用振幅校正因子更新脉冲响应和震源子波，具体为：

$\begin{array}{c}{S}_{i+1}=\frac{1}{\sqrt{\mathrm{\γ}}}{S}_{i+1}\\ \mathrm{\δ}{X}_{i+1}=\sqrt{\mathrm{\γ}}\mathrm{\δ}{X}_{i+1}.\\ {\stackrel{\‾}{X}}_{i+1}=\sqrt{\mathrm{\γ}}{\stackrel{\‾}{X}}_{i+1}\end{array}$

5.如权利要求3所述的一种基于稀疏反演的层间多次波压制方法，其特征在于：所述步骤4）利用振幅校正因子更新脉冲响应和震源子波，具体为：

$\begin{array}{c}{S}_{i+1}=\frac{1}{\sqrt{\mathrm{\γ}}}{S}_{i+1}\\ \mathrm{\δ}{X}_{i+1}=\sqrt{\mathrm{\γ}}\mathrm{\δ}{X}_{i+1}.\\ {\stackrel{\‾}{X}}_{i+1}=\sqrt{\mathrm{\γ}}{\stackrel{\‾}{X}}_{i+1}\end{array}$

6.如权利要求2或3或4或5所述的一种基于稀疏反演的层间多次波压制方法，其特征在于：所述步骤6）中预测得到的层间多次波M的计算公式为：

$M=-\stackrel{\‾}{X}\mathrm{\δ}{X}^H\stackrel{\‾}{P}.$

Images