CN104077749A - 一种基于轮廓波变换的地震数据去噪方法 - Google Patents

Abstract

本发明涉及一种基于轮廓波变换的地震数据去噪方法，其特征在于包括如下步骤：1、读取地震数据并进行拉普拉斯金字塔分解；2、将拉普拉斯金字塔分解得到的各个尺度的高频子带输入方向滤波器组，得到各个尺度的高频子带在各个方向上的分量；3、根据地震信号中噪声的特性选择噪声模型，然后对其进行塔形方向滤波器组分解并求取阈值，用阈值对各个尺度的高频子带在各个方向上的分量进行滤波；4、对滤波后的各个尺度的高频子带在各个方向上的分量进行塔形方向滤波器组逆变换，得到去噪后的地震信号。本发明将轮廓波变换应用于地震数据处理，克服了小波变换、曲波变换等方法在去除噪声的同时对有用信号造成损伤的缺点，提高了地震数据的保真度，进而压制噪音，提取有效信息。

Description

一种基于轮廓波变换的地震数据去噪方法

技术领域：

本发明涉及一种基于轮廓波变换的地震数据去噪方法，属石油地震勘探地震资料处理技术领域。

背景技术：

随着科技的不断发展，石油的战略地位日趋重要。高密度地震技术是国内外发展比较迅速的物探技术之一，广泛应用于石油资源的勘探。然而，随着地下地质构造复杂性的增加，油气勘探的难度越来越大，对地震数据处理的质量要求也越来越高。

在早期，传统小波变换一直是图像去噪领域乃至更广阔的图像处理领域的重要手段。但是，图像信号是一种二维信号，而小波的最优表示则局限于一维信号，由可分离小波构成的二维函数不能有效地描述数字图像信号的形状信息，即几何特征，因此在小波域不能有效地将噪声信号和图像物体的几何特征有效地区别开来。

轮廓波变换的技术日益成熟，其应用逐渐推广到各个领域。如申请号为CN201210008436的专利申请公开了一种基于轮廓波变换模型的视频指纹方法，特征是先将一段视频进行标准化处理，对每帧图像进行分块；然后对每个图像块采用多尺度和多方向轮廓波变换处理，并采用隐马尔科夫树技术提取模型参数；再将每块的标准方差矩阵进行奇异值分解，提取最大的奇异值，将最大奇异值级联并进行归一化后构成视频指纹矢量。其提高了视频指纹的稳健性、识别率和搜索效率，可广泛用于视频文件的复制检测、检索、识别和版权保护。申请号为CN201010288552的专利申请公开了一种多源图像融合方法，在前期特征提取阶段将基于小波核的支持度变换和抗混叠轮廓波结合在一起，一方面增加了支持度变换无法提取的方向信息，另一方面消除了轮廓波变换中产生的混叠现象；在图像融合阶段，采用脉冲耦合神经网络构建高频判断器对各融合图像低频信号进行融合判断，采用绝对值最大选取规则对各融合图像高频信号进行融合判断，实现了对高低频信号的有效融合。

将轮廓波变换应用于地震数据处理，可以克服小波变换、曲波变换等变换方法在去除噪声的同时对有用信号造成损伤的缺点，可提高地震数据的保真度，对提高地震勘探的分辨率具有十分重要的意义。

发明内容：

为了克服现有技术的不足，本发明的目的在于提供一种基于轮廓波变换的地震数据去噪方法，将轮廓波变换应用于地震数据处理，克服了小波变换、曲波变换等方法在去除噪声的同时对有用信号造成损伤的缺点，提高了地震数据的保真度。

本发明是通过如下技术方案来实现上述目的的。

本发明所提供的一种基于轮廓波变换的地震数据去噪方法，包括如下步骤：

(1)、读取地震数据，对地震数据进行拉普拉斯金字塔分解得到地震信号在各个尺度的低频子带和表示地层细节部分的高频子带；

(2)、将拉普拉斯金字塔分解得到的各个尺度的高频子带输入方向滤波器组，得到各个尺度的高频子带在各个方向上的分量；拉普拉斯金字塔和方向滤波器组合称为塔形方向滤波器组；对信号进行轮廓波变换即对信号进行塔形方向滤波器组分解；

(3)、根据地震信号中噪声的特性选择噪声模型，然后对噪声模型进行塔形方向滤波器组分解；根据分解得到的分解系数求取阈值，并用阈值对由步骤2中得到的各个尺度的高频子带在各个方向上的分量进行滤波；

(4)、对滤波后的各个尺度的高频子带在各个方向上的分量进行塔形方向滤波器组逆变换，得到去噪后的地震信号。

所述的对地震数据进行拉普拉斯金字塔分解是对原始地震数据采用拉普拉斯金字塔进行尺度分解，在每一尺度都产生一个原信号下采样的低频子带和高频子带；

所述的将拉普拉斯金字塔分解得到的各个尺度的高频子带输入方向滤波器组，得到各个尺度的高频子带在各个方向上的分量；是将得到的高频子带通过方向滤波器组，可以产生2ⁱ个方向子带，其中i表示尺度；

利用方向滤波器组可以获得轮廓波各个尺度的方向分解；

所述的用阈值对由步骤2中得到的各个尺度的高频子带在各个方向上的分量进行滤波；是根据对噪声模型进行塔形方向滤波器组分解得到轮廓波系数求取的阈值，对步骤2中得到的各个尺度的高频子带在各个方向上的分量进行滤波，去除噪音成分，保留有用信息，达到滤波的目的；

所述的对滤波后的各个尺度的高频子带在各个方向上的分量进行塔形方向滤波器组逆变换，得到去噪后的地震信号，是其轮廓波系数经过滤波处理后，其数据复杂度大大降低，噪音信号有效的被去除，故而得到去噪后的地震数据。

本发明与现有的技术相比，将轮廓波变换应用于地震数据处理，对地震数据和噪声模型利用拉普拉斯金字塔和方向滤波器组进行分解，得到轮廓波系数，然后设置阈值去噪，克服了小波变换、曲波变换等方法在去除噪声的同时对有用信号造成损伤的缺点，提高了地震数据的保真度，进而压制噪音，提取有效信息。

附图说明:

图1是本发明基于轮廓波变换的地震数据去噪方法的流程示意图。

图2是拉普拉斯金字塔分解的流程图。

图3是梅花形滤波器组的结构示意图。

图4是原始地层剖面图1。

图5是原始地层剖面图2。

图6是原始地层剖面图3。

图7是拉普拉斯金字塔逆变换的流程图。

图8是原始地层剖面图1经轮廓波处理后的剖面图。

图9是原始地层剖面图2经轮廓波处理后的剖面图。

图10是原始地层剖面图3经轮廓波处理后的剖面图。

具体实施方式：

下面结合附图和实施例对本发明作进一步说明。

实施例：

本发明在该实例中的具体步骤如下：

步骤101：读取文件中的地震数据，对地震数据进行拉普拉斯金字塔分解，得到地震信号在各个尺度的低频子带和表示地层细节部分的高频子带。

先调用自编的读取数据文件的函数程序，读取地震数据的文件头，可以得到地震数据的道数为512道和每道的采样点数512个采样点；接着对地震数据进行拉普拉斯金字塔分解。

进行拉普拉斯金字塔分解时，首先，对原始地震信号采用低通滤波器进行低通滤波，然后进行下采样得到低频信号，即得到原始地震信号的近似分量a；再对该低频信号进行上采样，使该近似分量与原始地震信号等长；然后用综合滤波器对上采样的低频信号进行低通滤波，以平滑上采样所插进的0。将滤波后的地震信号与原始地震信号进行差分，得到拉普拉斯金字塔分解后的高频信号b。整个分解流程如图2所示。这个过程可以在下采样等得到的低通信号上迭代进行，即得到信号的多尺度分量并捕捉奇异点。

步骤102：将拉普拉斯金字塔分解得到的各个尺度的高频子带输入方向滤波器组，得到各个尺度的高频子带在各个方向上的分量；拉普拉斯金字塔和方向滤波器组合称为塔形方向滤波器组；对信号进行轮廓波变换即对信号进行塔形方向滤波器组分解；

对得到的高频子带采用二维方向滤波器组来捕捉方向信息，将其输入方向滤波器组即可得到2j个方向分量，其中j表示尺度；

本发明采用的是由梅花形滤波器组和旋转的信号重采样模块组成的方向滤波器组来进行捕捉方向信息。梅花形滤波器组的结构如图3所示，

其中：梅花形滤波器组为采样矩阵且使信号旋转，有Q₀和Q₁两种形式，其中Q₀使信号旋转45度，Q₁使信号旋转-45度，Q₀和Q₁的形式如下：

$Q_0=\begin{array}{}\left(\begin{array}{cc}1& -1\\ 1& 1\end{array}\right)\end{array}$ $Q_1=\left(\begin{array}{cc}1& 1\\ -1& 1\end{array}\right)$

本发明采用Q₀的形式。

另外信号重采样模块是在梅花形滤波器组对信号分解之前进行的，并且在合成阶段后进行一个相反的操作，其主要作用就是对信号进行重采样。信号被采样后旋转，同时宽度增大为原来的两倍。在该模块中可采用以下四种采样矩阵

$R_0=\left(\begin{array}{cc}1& 1\\ 0& 1\end{array}\right)$ $R_1=\left(\begin{array}{cc}1& -1\\ 0& 1\end{array}\right)$ $R_2=\left(\begin{array}{cc}1& 0\\ 1& 1\end{array}\right)$ $R_3=\left(\begin{array}{cc}1& 0\\ -1& 1\end{array}\right)$

从中可看出，有两组互为转置矩阵，故本发明采用R₀和R₁两个矩阵，再加上转置变化即可实现对信号进行重采样的目的。

步骤103：根据地震信号中噪声的特性选择合适的噪声模型，然后对噪声模型进行塔形方向滤波器组分解；根据噪声模型的塔形方向滤波器组分解得到的分解系数求取阈值，并用该阈值对由步骤2中得到的各个尺度的高频子带在各个方向上的分量进行滤波。

由图4可以看出，原始地震数据受随机噪声影响严重，所以我们选取了包含系统噪声、次生噪声、环境噪声等的随机噪声作为噪声模型，根据随机噪声在一定的区域范围内遵循统计规律的这一特性对其进行滤除。

首先，将噪声模型进行塔形方向滤波器组分解，并将噪声模型的分解系数转化为Matlab中的胞元数组；由于轮廓波变换是线性变换，在随机噪声的优势方向子带上反映的随机噪声的分解系数值较大；且由于轮廓波的稀疏性，在该子带上反映随机噪声的稀疏的数目有限。

采用常规的阈值求取方法，其表达式为：

th＝3σ

$\mathrm{\σ}=\frac{\mathrm{Median}\left(\right|c_{i,j}\left|\right)}{k},i=l_k...l_q;j=m,...n,$

其中：k表示尺度，c_i,j为随机噪声的第i个需要处理尺度上的第j个随机噪声方向子带，Median()为取中值函数，th表示阈值，σ表示随机噪声的标准差，对于不同的方向子带，σ值不同。

接着，对由步骤2中得到的各个尺度的高频子带在各个方向上的分量进行滤波；对于地震信号，轮廓波变换的最低尺度层对应的是最低频带，其得到的分解系数主要反映的是地震信号的整体轮廓，一般不对最低尺度层进行处理；当地震信号的高频子带在各个方向上的分量的分解系数大于阈值th时，认为该稀疏反映的是随机噪声成分应去除掉；当地震信号的高频子带在各个方向上的分量的分解系数小于th时，则应保留。

在最佳尺度层中，该层有效信号的能量比R达到某一较高的值(在实验中,R取为0.9)时,认为该层上有效信号成分非常多，噪声成分少,则不对该尺度层做处理，以避免有效信号的损失。

(4)、对滤波后的各个尺度的高频子带在各个方向上的分量进行塔形方向滤波器组逆变换，得到去噪后的地震信号。

步骤104：对滤波后的各个尺度的高频子带在各个方向上的分量进行塔形方向滤波器组逆变换，得到去噪后的地震信号。

对滤波后的胞元数组转化为塔形方向滤波器组分解系数，然后再对这些分解系数进行方向滤波器组逆变换，最后再对方向滤波器组逆变换的结果进行拉普拉斯金字塔逆变换，如图7，即得到去噪后的地震信号。

本发明在如表1所示的测试环境下进行实验。

表 1：实验测试环境

图4、图5、图6是原始地层剖面图，从图中可以看出，原始数据随着地层深度的增加受随机噪声的影响越发严重，以至影响有用地层信息的识别。因此，我们设计了随机噪声的噪声模型，将地震数据和噪声模型都进行轮廓波变换，得到了地震数据和噪声模型的轮廓波系数，再根据噪声模型的轮廓波系数确定噪声在轮廓波域集中分布的尺度和方向，通过噪声模型的轮廓波系数设置阈值，并用此阈值对地震数据的轮廓波系数进行滤波，图8、图9、图10是原始地层剖面图经轮廓波处理后的剖面图。与图4、图5、图6的原始地层剖面图进行对比可知，通过本发明的方法处理后，地震数据中的地层信息变得更加清晰，部分在处理前无法识别的薄地层在处理后显现出来，表2为地震数据处理前后的信噪比对比结果。

表 2：地震数据处理前后的信噪比对比结果

从表2的测试结果可以看出，本发明的方法能够有效地压制随机噪声干扰，保护有用的反射地震信号，通过本发明的方法处理后，去噪后的地震信号的细节信息变得更加明显，信噪比明显提高。

通过以上实验和分析可知，将本发明的方法运用于地震信号的去噪过程中，能够充分发掘噪声在轮廓波域中的分布特点，压制噪声干扰。

Claims (5)

1.一种基于轮廓波变换的地震数据去噪方法，其特征在于包括如下步骤：

(1)、读取地震数据，对地震数据进行拉普拉斯金字塔分解得到地震信号在各个尺度的低频子带和表示地层细节部分的高频子带；

(2)、将拉普拉斯金字塔分解得到的各个尺度的高频子带输入方向滤波器组，得到各个尺度的高频子带在各个方向上的分量；拉普拉斯金字塔和方向滤波器组合称为塔形方向滤波器组；对信号进行轮廓波变换即对信号进行塔形方向滤波器组分解；

(3)、根据地震信号中噪声的特性选择噪声模型，然后对噪声模型进行塔形方向滤波器组分解；根据分解得到的分解系数求取阈值，并用阈值对各个尺度的高频子带在各个方向上的分量进行滤波；

(4)、对滤波后的各个尺度的高频子带在各个方向上的分量进行塔形方向滤波器组逆变换，得到去噪后的地震信号。

2.根据权利要求1所述的一种基于轮廓波变换的地震数据去噪方法，其特征在于所述的对地震数据进行拉普拉斯金字塔分解是对原始地震数据采用拉普拉斯金字塔进行尺度分解，在每一尺度都产生一个原信号下采样的低频子带和高频子带。

3.根据权利要求2所述的一种基于轮廓波变换的地震数据去噪方法，其特征在于将拉普拉斯金字塔分解得到的各个尺度的高频子带输入方向滤波器组，可以产生2ⁱ个方向子带，其中i表示尺度。

4.根据权利要求3所述的一种基于轮廓波变换的地震数据去噪方法，其特征在于所述的方向滤波器组由梅花形滤波器组和旋转的信号重采样模块组成。

5.根据权利要求4所述的一种基于轮廓波变换的地震数据去噪方法，其特征在于采用常规的阈值求取方法，其表达式为：

th＝3σ

$\mathrm{\σ}=\frac{\mathrm{Median}\left(\right|c_{i,j}\left|\right)}{k},i=l_k...l_q;j=m,...n,$

其中：k表示尺度，c_i,j为随机噪声的第i个需要处理尺度上的第j个随机噪声方向子带，Median()为取中值函数，th表示阈值，σ表示随机噪声的标准差，对于不同的方向子带，σ值不同。