基于经验小波变换的地震资料噪声压制方法

覃发兵 徐振旺 啜晓宇 张小明 郭乃川 董玉文 陈 伟

（ 1长江大学管理学院；2中国石油辽河油田公司勘探开发研究院；3河北煤炭科学研究院；4 中国石油集团东方地球物理公司大港物探处；5中海石油（中国）有限公司天津分公司渤海石油研究院；6中国石油集团东方地球物理公司研究院资料处理中心；7非常规油气湖北省协同创新中心 ）

目前实际生产中常用的地震资料去噪方法有频率域滤波方法、频率波数域滤波方法、频率空间域滤波方法、基于Radon变换的去噪方法、聚束滤波方法、基于小波分解和重建的去噪方法，以及局部径向道中值滤波方法和傅里叶相关系数滤波方法等[1-3]；其他去噪方法还有多项式拟合、K—L变换和矢量分解等方法。其中频率域滤波方法主要是在地震资料采集时进行去假频处理和对含工业干扰的地震资料进行限波处理[4-5]。频率波数域滤波方法即f—k滤波主要用于去除地震资料中的面波，它优于一维的滤波方法在于采用了扇形滤波器[6-7]。频率空间域滤波方法包括由国九英提出的ω—x域算子外推方法[8]和蔡加铭等提出的f—x域算子外推方法[9]，这两种外推方法能有效去除面波和线性干扰波。基于Radon变换的去噪方法主要有线性拉东变换去噪方法、抛物线拉东变换去噪方法和双曲线拉东变换去噪方法，线性法方法适合于具有线性时差的同相轴，双曲线方法适合于具有正常时差的同相轴，抛物线方法居中。聚束滤波方法从叠加基础上发展而来，其能对约束条件进行直接控制，能根据数据调整模型并避免畸变，所以这种方法实际上是一种基于模型的波场分解方法。基于小波分解和重建的去噪方法适用于深层地震资料的面波去除。基于Radon变换的去噪方法、聚束滤波方法、基于小波分解和重建的去噪方法，以及局部径向道中值滤波方法由Pangs等[10]提出，不会轻易丢失数据，不需要复杂的插值算法，能有效去除线性噪声。傅里叶相关系数滤波方法由Douglas[11]提出，它利用了相关系数谱，能够完整地保存有效信号并减少地震数据中的线性倾斜相关能量。

上述的常规地震资料去噪方法均在一定程度上存在着噪声去除不干净、有效信号丢失等问题，归根到底是因为这些算法不是自适应算法，即算法本身不能识别出信号中的噪声成分，尤其在面对非线性非平稳的地震信号时，常规的去噪方法很难提高地震资料的信噪比，且容易产生虚假信号和假频。

为了解决非线性非平稳信号处理问题，需要寻找自适应算法。最有可能实现自适应算法的是基于连续尺度展开的基追踪小波理论体系[12]，该体系提供了自适应的时频分布，且在很多方面都有很好的应用，但其在对信号分段时需要固定周期，从而限制了自适应性。小波体系中的另外一种自适应算法是Malvar—Wilson小波，这是一种试着将时域信号进行分段处理，使得每一段信号包含有不同频谱信息的方法[13]。尽管这个想法非常好，但是在时间域对信号进行分割是很难有效做到的。Meyer等[14]提出了一种叫作梳状波的方法，它是通过在傅里叶域建立一个自适应的滤波器组来实现的。但是它依然建立在前述的Malvar—Wilson小波的基础上，只是在频率域将信号分割，而不是在时间域分割信号，而且其构造方法较复杂且依赖于事先设好的分割方法。最后要提的是Daubechies所做的名为“同步子波”的工作[15]，该方法包含经典的小波分析和时频域信息再分配思想。这种算法可以得到更加精确的时频分布，可以提取特殊的“模态”。上面提到的常规自适应方法，要么是预先设计好分割策略，要么就是将经典的小波分析进行巧妙的输出。

美籍华人Huang在1998年首次提出经验模态分解（EMD）方法用于对非线性非平稳信号进行自适应分解[16]，这种算法可以将信号分解成不同的固有模态分量（Intrinic Mode Function，简称IMF），Huang认为每个固有模态分量均代表原始信号中具有独立性质的分量。经验模态分解能够解决大部分的非线性非平稳问题，并且成功应用在各领域。但是Huang将分解后的IMF分别求取振幅谱后发现，对于强非线性非平稳信号而言，各IMF在频率域是很难完全分开的，即EMD算法具有模态混叠现象，由此也导致EMD分解得到的IMF数量太多，即存在分解出不包含在原始信号中的非本征信号问题。为了解决模态混叠现象，Huang将具有一定总体的高斯随机白噪加入到原始信号中并进行分解，且重复多次取平均，最后得到的IMF在频域不再重叠，克服了模态混叠现象，该方法就是总体经验模态分解（Ensemble Empirical Mode Decomposition，EEMD）算法[17-19]。EEMD算法解决了EMD算法存在的模态混叠问题，但是计算量大，且分解出的IMF的本质特性依旧依赖于EMD算法，因此其分解后的IMF也包含一些原始信号中不具有的本征信号，破坏了自适应性。近年在信号领域流行的经验小波变换（Empirical Wavelet Transform， 简称EWT）是一种全新的自适应分解算法[20-21]，其相较于EMD和EEMD算法能更好地分解出原始信号中固有的本征信号，具有更高的自适应性。另外，EWT算法建立在成熟的小波理论基础上，具有充分的数学理论基础，并且其借助于小波分解快速算法使得自身具有较高的计算效率。考虑到EWT是基于EMD的最新自适应算法，本文首先介绍了EMD和EEMD的基本原理，然后给出了EWT的基本框架，最后提出基于EWT的噪声压制算法流程并给出例证。本文首次将EWT算法应用到地震资料噪声压制领域，从数值模拟和实际资料两方面分别进行了研究。

1 方法原理

1.1 经验模态分解（EMD）和总体经验模态分解（EEMD）

对于给定的信号s（t)，经验模态分解的算法流程如下[8]：

（1）确定原信号s（t)所有的极大值点和极小值点，将所有极大值点和极小值点分别用光滑曲线联接起来，将这两条曲线分别作为s（t)的上下包络线，并计算上下包络线的平均值曲线m（t)，用原始信号s（t)减去m（t)得到h1（t)，即：

往往h1（t)还不是一个严格意义上的IMF分量。

（2）将h1（t)作为原信号重复步骤⑴的过程，得到：

反复筛选k次，直到h1，k（t)变为一个固有模态分量，最后一次迭代如下：

这样就从原信号中分离出了第一个固有模态分量c1（t)，即：

本文采用Huang等人提出的一种仿柯西收敛准则来终止循环，首先令：

当0.2≤SD≤0.3时循环结束。

（3）从原始信号中减去c1（t)得到第一阶剩余分量r1（t)，即：

一般r1（t)中还含有很多固有模态分量，因此需要将r1（t)作为原始信号重复上面的步骤，这样依次得到第二阶固有模态分量至第N阶固有模态分量，以及第二阶剩余分量至第N阶剩余分量，即：

当第N阶固有模态分量cN（t)或者当第N阶剩余分量rN（t)（能量）足够小时终止整个分解过程。

综合上述分解过程，可知原信号分解如下：

总体经验模态分解的理论基础是经验模态分解。总体经验模态分解的算法如下：

（1）给原始信号s（t)加上白噪声w（t)得到总体S（t)：

（2）对S（t)进行经验模态分解，得到各个固有模态分量：

（3）给目标信号加入不同的白噪声wj（t)：

重复以上步骤，得到各总体的固有模态分量组：

（4）取相应固有模态的均值作为最终的固有模态组：

（5）取相应rn的均值作为最终的剩余分量：

1.2 经验小波变换

小波变换在目前信号处理中用得最多，其具有完善的理论，并被广泛应用到了各个领域。现将信号s（t)的傅里叶变换及其逆变换分别记为s和s。在时间域，一个小波集｛Ψu，ν｝是由具有零均值的母小波函数Ψ加上尺度因子ν（ν＞0）和平移因子u （u ∈R)来定义的，即有：

将信号s（t)与上述的小波函数作内积运算即可得到信号s（t)的小波变换，即有：

如果尺度因子ν是一个连续的变量，则公式（14）即为连续小波变换；如果ν是离散的，令ν=aj，则公式（14）变为：

公式（15）即为离散小波变换的表达式。小波变换的一个有用的性质是其可以被当作一个滤波器组（每一个滤波器对应一个尺度），一般认为尺度因子ν=2j。如果从傅里叶角度出发，构造小波变换等于构造一系列的带通滤波器。为了自适应地处理信号，可以重点研究信号的局部成分对应的频域信息。1.1节所述单分量信号IMF的性质表明IMF的频谱具有紧支性。为清楚起见，仅仅只考虑实信号（频谱关于圆频率ω=0对称），但下文的讨论依然能简单地通过在正频率域和负频率域分别建立不同的滤波器来推广到复信号领域。考虑一个周期为2π的归一化的傅里叶轴，为了满足香农定理（Shannon theory），将讨论区间限制在ω∈[0，π]。

首先将傅里叶区间[0，π]分割成N个连续的小区间，共有N+1个间断点，第n个间断点用ωn来表示，第一个间断点为ω0=0，最后一个间断点为ωN=π（图 1），第n个小区间为 Λn=[ωn-1，ωn]，从而显而易见有Λn=[0，π]。以ωn为中心，可以定义宽度为2τn的过渡带Tn（图1中灰色阴影部分）。

图1 傅里叶频率轴分割示意图

经验小波是定义在区间Λn上的带通滤波器。借用Littlewood—Paley小波和Meyer小波中的构造方法，对于∨n＞0，通过公式（16）和公式（17）来定义经验尺度函数和经验小波：

公式（16）和公式（17）中β（x）的函数值要么为0，要么为1，且有：

许多函数均满足公式（18）中的性质，一般选用下面的表达式：

对于τn参数的选择则有几种可能性，一般认为τn是ωn的一部分，即有：τn=γ ωn，且0＜γ＜1。这样，对于∨n＞0，公式（16）和公式（17）可以简化为：

怎样将傅里叶谱进行分段在经验小波变换中至关重要，其直接关系到对原始信号分解后的自适应程度。经验小波变换将原始信号进行不同的分割，比如对某个中心频率的紧支撑部分进行分割。假设断点数目为N，这意味着需要N+1个边界。除了起点0和终点π以外，还需要N-1个边界。为了找到这些边界，首先对信号频谱的局部极大值点进行降序排列（起点0和终点π包括在内）。假设找到了M个极大值点，将会出现下面两种情况：

（1)M≥N：算法发现了足够的极值点以便于分割原始信号，但本文只取前N-1个极大值点；

（2)M＜N：信号没有预期的那么多模态，保留这M个极值点，并添加一些近似值直到极值点达到N个。

这样，加上起点0和终点π，就定义了每一个间断的边界ωn作为相连两个极大值点的中点。从前面的论述可以知道如何才能建立一个紧支撑的经验小波。仿照经典的小波变换理论来定义经验小波变换，其细节系数由信号与经验小波的内积得到：

近似系数（采用传统记法Wsε（0，t））通过信号与尺度函数的内积得到：

这里Ψn（ω）和φ1（ω)为公式（20）和公式（21）所定义。这样原始信号可以通过下面的公式来重构：

为简单起见，经验小波变换所蕴含的经验模态函数可以定义为：

本文利用经验小波变换进行地震噪声压制处理是按单道进行，其具体流程如下：

（1）选取合适的小波函数并利用EWT算法对目标单道地震信号s（t)进行自适应分解得到其各个频率尺度的固有模态分量s0（t)，s1（t)，…，sk（t)，…；

（2）对原始单道地震信号以及分解出的固有模态分量信号作频谱分析，根据本文给定的主频计算公式求取各自的主频，并根据原始地震信号的主频设定频率阈值范围，主频的计算公式如下：

式中 df——信号主频；

f——频率分量；

s（f)——振幅谱。

（3）选取主频值在阈值范围内的固有模态分量进行重构，最终获取噪声得到压制的地震信号；

（4）对二维或三维地震数据体中的每一道均执行步骤（1)至（3)，得到最终的噪声压制数据体。

2 数值试验

2.1 单道地震信号去噪模拟

本文提出的噪声压制方法是基于单道进行的，本节选取文献[1]中的单道合成地震记录来测试基于EWT算法的地震资料噪声压制方法，以小波变换、EMD、EEMD3种去噪方法作为对比。所选取的单道地震信号如图2a所示，该信号由主频为20Hz的余弦信号在局部叠加不同主频的子波信号形成的：在0.3s处添加频率为100Hz的Morlet子波信号，在 1.07s和 1.1s添加 30Hz的雷克子波信号，在1.3～1.7s之间添加 7Hz、30Hz和 40Hz 3种不同频率叠加的复合信号。需要说明的是，1.3～1.7s之间添加的7Hz频率成分是不连续的，且均持续不到一个周期的时间，分别在1.37s、1.51s和1.65s出现。图2a给出的单道信号具有典型的非线性非平稳特征，可以用来很好地验证自适应性算法的优劣。为了更好地说明问题，本文定义了如下噪信比公式：

式中 NSR——噪信比；

Snoise——噪声；

Ssignal——信号；

对加噪后的地震记录分别进行EMD、EEMD、小波变换和EWT去噪测试，其中EEMD所添加高斯噪声的相对标准差为0.6，总体个数为300，结果如图3所示。从图中很清晰地看到，EMD和EEMD对叠加的高频子波（相对背景信号的主频而言）很难识别，而小波变换和EWT能较好地识别，如图3中红框所示。但是小波变换对高频子波的形态勾画还是不如EWT精确，图中能看到小波变换的结果使得高频子波的旁瓣振幅较大，而EWT压制了旁瓣，后者与实际子波形态吻合较好。从图3中蓝框还能明显看出，当添加的非平稳子波主频与背景信号主频相差不大时，小波变换很难恢复添加的子波的形态，而EMD、EEMD和EWT此时具有较好的子波恢复能力。结果表明，EWT能在较好地恢复信号简谐成分的同时更好地刻画信号的局部非平稳细节，另外，EWT吸收了EMD/EEMD和小波变换共同的优点，即在非线性非平稳的情况下既能分辨信号中的高频成分，又能识别低频成分，这一点从EWT的定义中可以看出，即EWT是建立在小波框架基础上的自适应分解算法。比较EEMD和EWT对信号的分解结果（图4、图5）可以看出，对于同样的信号，EWT分解出的固有模态分量要少于EEMD，且与原始信号中所含有的频率成分基本对应，第二个分量和第三个分量的结果主要包含在1.65s时的40Hz子波和其他更高频率的组分，第四个分量反映了1.1s附近主频为30Hz雷克子波、1.4s时主频为30Hz子波及40Hz子波的残余部分，20Hz余弦背景信号主要反映在第五个分量中。分解出的固有模态分量与原始信号组分几乎对应，主要源于EWT分解结果不包含模态混叠的成分，这一点还可以从分解出的固有模态分量的主频分布看出。

图2 人工合成的地震记录(a)及加入噪信比为0.3的高斯白噪声后的测试信号(b)

图3 4种去噪结果

图4 EEMD对信号的分解

图5 EWT对信号的分解结果

通过去噪前后的时频谱可以看出去噪效果。借助S变换时频谱制作方法求取了原始不含噪声的单道信号、加噪后单道信号的S变换时频谱（图6)。从图6a可以看出，原始不含噪声的信号时频谱比较“干净”，时频谱能完整地体现各掺杂的子波成分，其中20Hz背景信号的时频轴分布非常稳定，其他主频的子波在图中也能清晰分辨出来；图6b显示了加噪后的S变换时频谱灰度图，由于噪声的存在，图中能很明显看出噪声的时频能量团，如图中黑色椭圆所示，另外噪声还影响了20Hz背景信号的时频轴光滑度，如图6b中蓝框所示，背景信号的时频轴出现了毛刺，且首尾出现了明显的能量泄露，其他子波的能量团也出现了能量泄漏的现象。EMD、EEMD、小波变换以及EWT 4种方法去噪后信号的S变换时频谱在统一色标下的时频图如图7所示。从图7中很清晰地看到，EEMD方法得到的时频谱（图7b)显示的子波信号和背景信号要比EMD结果（图7a)更加清晰，但是对0.3s处100Hz的高频子波不能识别，如图中蓝框所示。小波变换结果（图7c）同样也不能识别高频子波，但其显示的背景信号同相轴不存在首尾能量泄漏现象，如图中蓝框所示。EWT结果（图7d)能更加清晰地分辨高频子波，且背景信号同相轴不存在能量泄漏现象。对比4种方法可以看出，在噪声压制方面EWT方法具有小波变换和EMD、EEMD两种时频体系的优点，且克服了缺点，是一种优秀的噪声压制方法。

图6 原始无噪声信号(a)和原始含噪信号的S变换时频谱(b)

图7 4种去噪结果的S变换时频谱

2.2 二维合成地震剖面去噪模拟

为了考察各种方法在空间上的连续性，对合成的二维地震剖面（图8a)加上随机噪声后（图8b)分别用EMD、EEMD、小波变换以及EWT进行噪声压制数值实验，结果如图9所示，其中EEMD方法所添加高斯噪声的相对标准差为0.6，总体个数为300。二维合成地震剖面（图8a)包含线性连续同相轴、线性间断同相轴、弯曲连续同相轴、断层共4种地质特征；当加入噪声后，线性同相轴中的间断点（图8b中红色椭圆和红色箭头所示）以及断层构造（图8b中蓝色椭圆所示）很难被识别出来。EMD、EEMD去噪结果（图9a、b）相比小波变换（图9c）和EWT（图9d）分辨率较低，信噪比低，但是各同相轴的间断点更清晰，这也体现了EMD和EEMD是一类能解决非线性问题的关键方法，而EEMD相对EMD同相轴更加稳定，损失能量较少。小波变换去噪结果很难识别间断点（图9c中红色、蓝色椭圆以及箭头所示），但地震剖面分辨率较高。EWT去噪结果无论是从分辨率还是识别间断点能力上都是最强的，其在二维地震剖面噪声压制时吸收了EMD/EEMD和小波变换的优点，在空间连续性方面表现优越。对4种方法去噪结果与原始不含噪声地震剖面的残差剖面进行比较可以发现，EWT损失有效信号最少，且EWT残差剖面很难见到有效信号的“影子”，是一种保幅高效的噪声压制方法。

图8 合成的二维地震剖面(a)及加上随机噪声后的剖面(b)

图9 合成二维4种不同方法的去噪结果

3 实际资料处理

为了考察各种方法在实际资料噪声压制的能力，选用了一个含明显随机噪声的地震剖面，如图10所示。该剖面含有明显背斜和断层构造（图10中红色箭头所示），但是受随机噪声的干扰同相轴连续性较差，间断点不清晰，很难在该剖面上进行人工解释。与数值模拟类似，分别采用EMD、EEMD、小波变换及EWT对该实际资料进行噪声压制实验，结果如图11所示，其中EEMD方法所添加高斯噪声的相对标准差为0.6，总体个数为300。EMD、EEMD去噪结果（图11a、b）相比小波变换（图11c）和EWT（图11d）分辨率较低，信噪比低，但是断层和背斜构造更清晰，这里也能看出EEMD相对EMD同相轴更加稳定，损失能量较少。小波变换去噪结果很难识别背斜的顶（图11c中红色箭头所示)，但地震剖面分辨率相对EMD/EEMD较高。EWT去噪结果无论是从分辨率还是识别特殊构造的能力上都是最强的，尤其是同相轴的连续性大幅提高，这为人工地震解释做了很好的铺垫。分别作出4种方法去噪结果与原始不含噪声地震剖面的残差剖面（图12）可以发现，EMD、EEMD、小波变换的去噪结果明显损失了有效信息，残差剖面上有明显的有效信号的痕迹，而EWT损失有效信号最少，其残差剖面很难见到有效信号的“影子”，体现了其在处理实际资料时具有明显的保幅性。EMD、EEMD、小波变换以及EWT4种方法在处理地震资料所耗时间分别为460.32s、1606.41s、363.42s、521.58s，可见 EWT 算法在保证去噪效果的同时计算效率也较高，能够适用于大型实际资料噪声压制处理任务。

图10 含噪声的实际地震资料

图11 实际资料的4种去噪结果

图12 4种去噪结果的残差剖面

4 结论及建议

EMD和EEMD算法的地震资料噪声压制效果要比常规的信号处理算法优越，但是EMD具有模态混叠现象，EEMD在信号重构后的结果具有不完备性，且二者没有充实的数学理论基础。近年在信号领域流行的EWT算法是一种全新的自适应分解算法，其相较于EMD/EEMD算法能更好地分解出原始信号中的固有模态函数，具有更高的自适应性。另外，EWT算法建立在成熟的小波理论和经验模态分解的基础上，具有充分的数学理论基础，并且其借助于小波分解快速算法使得自身具有较高的计算效率。通过与EMD、EEMD及小波变换3种方法的对比可以看出，EWT方法在噪声压制后的地震剖面信噪比更高，尤其是在间断点和断层的识别方面具有明显的优势，且计算效率较高，是一种优越的噪声压制新方法。