CN114676734A - 一种基于wpd-emd-wpd的微震信号降噪方法 - Google Patents

Abstract

本发明公开了一种基于WPD‑EMD‑WPD的微震信号降噪方法，包括以下步骤：步骤S1：获取含噪微震信号X＝{s(1),s(2),...,s(N)}，采用采用小波包变换的降噪方法对含噪微震信号进行预降噪,得到预降噪后的微震信号X′＝{s(1)′,s(2)′,...,s(N)′}，其中N为信号长度；步骤S2：将预降噪后的微震信号X′＝{s(1)′,s(2)′,...,s(N)′}进行经验模态分解，得到一组频率由高到低的IMF分量：imf₁,imf₂,...,imf_M；通过各个IMF分量与预降噪后微震信号X′的相关系数找出噪音与有效信号的分界IMF分量；步骤S3：对分界之前的IMF分量进行小波包变换降噪,降噪后的IMF分量与剩余IMF分量重构得到降噪后的信号X″＝{s(1)″,s(2)″,...,s(N)″}。本发明克服了小波包变换缺乏自适应性，又解决了低信噪比下EMD分解质量低的问题，为微震信号的后续分析提供了准确条件。

Description

一种基于WPD-EMD-WPD的微震信号降噪方法

技术领域

本发明涉及水电站安全与信号处理交叉领域，特别涉及低信噪比、随机非平稳水电站洞室微地震信号的噪声压制方法，是一种基于WPD-EMD-WPD的微震信号降噪方法.

背景技术

微震监测技术起源于矿山工程，微震作为一种三维“体”监测方法，利用传感器接收岩石微破裂信号，进而捕捉岩体失稳破坏的前兆信息。微震监测技术能够实时捕捉岩体内部破裂信息，直接揭示裂缝萌生、发育、扩展、贯通直至宏观变形破坏的演化过程，是研究岩体损伤演化过程的重要工具。现在主要用煤矿开采，大型水电站，地下洞室开挖等大型工程。由于微震监测系统的使用，我国每年由于煤矿开采而死亡的人数剧减。然而实际工程中工况复杂，各通道拾取的信号会受到噪音干扰，影响初至时刻拾取、震源定位及参数反演。因此，微震信号降噪是开展微震监测及岩体稳定性评价的基础，具有非常重要的工程意义。

微震信号具有非平稳性、随机性、低信噪比、复杂性等特点，傅里叶变换虽能在一定程度抑制噪声，但其主要适用于平稳随机信号，对于突变的微震信号降噪效果不佳。目前用于非平稳信号的降噪方法有小波变换、经验模态分解(EMD分解)；但是，小波包变换存在缺乏自适应性的问题，低信噪比下EMD分解存在质量低的问题，对微震信号的降噪带来了诸多不便。

发明内容

本发明的目的在于克服现有技术的不足，提供一种基于WPD-EMD-WPD的微震信号降噪方法，克服了小波包变换缺乏自适应性，又解决了低信噪比下EMD分解质量低的问题，为微震信号的后续分析提供了准确条件。

本发明的目的是通过以下技术方案来实现的：一种基于WPD-EMD-WPD的微震信号降噪方法，包括以下步骤：

步骤S1：获取含噪微震信号X＝{s(1),s(2),...,s(N)}，采用小波包变换的降噪方法对含噪微震信号进行预降噪,得到预降噪后的微震信号X′＝{s(1)′,s(2)′,...,s(N)′}，其中N为采样点个数；

步骤S2：将预降噪后的微震信号X′＝{s(1)′,s(2)′,...,s(N)′}进行经验模态分解，得到一组频率由高到低的IMF分量：imf₁,imf₂,...,imf_M；通过各个IMF分量与预降噪后微震信号X′的相关系数找出噪音与有效信号的分界IMF分量；

步骤S3：对分界IMF分量之前的IMF分量进行小波包变换降噪,降噪后的IMF分量与剩余IMF分量重构得到降噪后的信号X″，设得到的信号X″＝{s(1)″,s(2)″,...,s(N)″}。

其中，WPD是指是小波包降噪方法，首先采用小波包降噪(WPD)方法进行预降噪,然后再用经验模态分解法(EMD)分解出IMF,通过各IMF与预降噪后微震信号X′的相关系数找出噪音与有效信号的分界IMF分量,最后对分界之前的IMF分量进行小波包降噪(WPD),降噪后的IMF分量与剩余IMF分量重构得到降噪后的信号；故本申请中整个流程本称为WPD-EMD-WPD降噪。

进一步地，所述步骤S1中，采用小波包变换的降噪方法对含噪微震信号进行预降噪时，选择能反映原信号特点的db9小波和3层的分解层数；

所述步骤S3中，对分界IMF分量之前的IMF分量进行小波包变换降噪时，选择db9小波和3层的分解层数。

进一步地，所述步骤S2包括以下子步骤：

S201.将预降噪后的微震信号X′＝{s(1)′,s(2)′,...,s(N)′}进行经验模态分解，得到一组频率由高到低的IMF分量：imf₁,imf₂,...,imf_M；

S202.计算第i个IMF分量与预降噪后的微震信号X′之间的相关系数：

其中，Cov为协方差函数，Var为方差函数；

第i个IMF分量相当于是由很多个点组成,把每个时间点对应的值带入公式计算得到该IMF分量与信号X′的相关系数；

S203.在i＝1,2,…,M时，重复执行步骤S202，得到每一个IMF分量与信号X′，其中与信号X′相关系数最大的IMF分量即为分界IMF分量。

进一步地，所述步骤S3包括：

S301.设imf₁,imf₂,...,imf_M中，第j个分量为分界IMF分量；j为1～M之间的整数；对于imf₁,imf₂...,imf_j中的每一个IMF分量进行小波包变换降噪，得到imf₁′,imf₂′...,imf_j′；

S302.将降噪后的IMF分量imf₁′,imf₂′...,imf_j′与剩余IMF分量imf_j+1,imf_j+2,...,imf_M相加，实现对降噪后的信号的重构，得到重构后的信号X″：

X″＝imf₁′+imf₂′+...+imf_j′+imf_j+1+imf_j+2,...,imf_M；

设得到的信号表示为X″＝{s(1)″,s(2)″,...,s(N)″}。

进一步地，所述微震信号降噪方法还包括：计算含噪微震信号X与降噪后的信号X″之间的信偏比，作为降噪效果的判断依据，信偏比SDR计算方式如下：

其中，s(n)为含噪微震信号X中第n个采样点的值，s(n)″为信号X″中第n个采样点的值，n＝1,2,…N。

本发明的有益效果是：本发明采用EMD分解结合小波包变换降噪方法具有自适应性,无需通过反复设置小波基函数,小波分解层数来达到降噪效果最优。解决了小波分解存在基函数类型选取,分解层数等不确定性问题；与EMD分解结合小波包的降噪效果较传统单一降噪算法效果较好,主要表现在信号的高频噪声被有效抑制,信号的细节信息丰富,保留了较多的有用信号。提出的信偏比能更好的反映降噪效果，且与信噪比和均方根误差有良好的对应性，信偏比越小这意味着噪声被抑制越多，从而与未含噪声的信号偏离程度越小，降噪效果越明显。

附图说明

图1为本发明的方法流程图；

图2为实施例中原始信号和加噪信号的波形图；

图3为实施例中原始信号和加噪信号的频谱图；

图4为实施例中小波包降噪波形图和频谱图；

图5为实施例中EMD降噪波形图和频谱图；

图6为实施例中本发明降噪波形图和频谱图；

图7为实施例中原始微震信号与本发明降噪后波形、频谱的对比示意图；

图8为实施例中微震信号EMD降噪波形图和频谱图；

图9为实施例中微震信号小波包降噪波形图和频谱图；

图10为实施例中原始微震信号和本发明降噪的微震信号时频谱图；

图11为小波分解与重构的原理示意图；

图12为实施例中个IMF分量的相关系数示意图。

具体实施方式

下面结合附图进一步详细描述本发明的技术方案，但本发明的保护范围不局限于以下所述。

本发明首先采用小波包对含噪微震信号预降噪该步骤抑制了噪声对EMD分解的影响，减轻了模态混叠现象；然后将预降噪后的微震信号EMD分解为一组频率由高到低的IMF，通过各IMF与预降噪后微震信号的相关系数找出噪音与有效信号的分界IMF分量，将分界之前的IMF分量进行小波包降噪；最后将降噪后的IMF分量与剩余IMF分量重构得到降噪后的信号。具体的：

如图1所示，一种基于WPD-EMD-WPD的微震信号降噪方法，包括以下步骤：

步骤S1：获取含噪微震信号X＝{s(1),s(2),...,s(N)}，采用小波包降噪方法对含噪微震信号进行预降噪,得到预降噪后的微震信号X′＝{s(1)′,s(2)′,...,s(N)′}，其中N为采样点个数；

步骤S2：将预降噪后的微震信号X′＝{s(1)′,s(2)′,...,s(N)′}进行经验模态分解，得到一组频率由高到低的IMF分量：imf₁,imf₂,...,imf_M；通过各个IMF分量与原信号的相关系数找出噪音与有效信号的分界IMF分量；

步骤S3：对分界之前的IMF分量进行小波包降噪,降噪后的IMF分量与剩余IMF分量重构得到降噪后的信号X″，设得到的信号X″＝{s(1)″,s(2)″,...,s(N)″}。

其中，WPD是指是小波包降噪方法，首先采用小波包降噪(WPD)方法进行预降噪,然后再用经验模态分解法(EMD)分解出IMF,通过各IMF与预降噪后微震信号X′的相关系数找出噪音与有效信号的分界IMF分量,最后对分界之前的IMF分量进行小波包降噪(WPD),降噪后的IMF分量与剩余IMF分量重构得到降噪后的信号；故本申请中整个流程本称为WPD-EMD-WPD降噪。

进一步地，所述步骤S1中，采用小波包降噪方法对对含噪微震信号进行预降噪时，选择能反映原信号特点的db9小波和3层的分解层数。

如图11所示，S₀相当于原始信号(X)，选择合适的小波基(db9小波)和分解层数(3层)就可以对S₀进行分解，图中就是进行了三层分解；

S₀被分解成了L₃,H₃,L₄,H₄,L₅,H₅,L₆,H₆.(S₀＝L₃+H₃+L₄+H₄+L₅+H₅+L₆+H₆),对分解系数阈值量化就是对L₃,H₃,L₄,H₄,L₅,H₅,L₆,H₆进行处理得到L’₃,H’₃,L’₄,H’₄,L’₅,H’₅,L’₆,H’₆，噪声的小波系数较小，并且噪声的小波系数要小于信号的小波系数，通过预先选取一个合适的阀值，大于阀值的小波系数被认为是有信号产生的，应予以保留，小于阀值的则认为是噪声产生的，置为零从而达到去噪的目的。重构的过程就是把L’₃,H’₃,L’₄,H’₄,L’₅,H’₅,L’₆,H’₆相加得到降噪后的信号S’₀。

所述步骤S3中，对分界IMF分量之前的IMF分量进行小波包变换降噪时，选择db9小波和3层的分解层数。

进一步地，所述步骤S2包括以下子步骤：

S201.将预降噪后的微震信号X′＝{s(1)′,s(2)′,...,s(N)′}进行经验模态分解，得到一组频率由高到低的IMF分量：imf₁,imf₂,...,imf_M；

在本申请的实施例中EMD分解的具体过程包括：

步骤A1：找出原始信号a(t)所有极值点，利用3次样条差值函数分别连接极大值点与极小值点，形成上下两条包络线；

步骤A2：计算上下包络线均值b(t)；

步骤A3：求差值c₁(t)＝a(t)-b₁(t)；

步骤A4：检查c₁(t)是否符合IMF的限定条件：若满足IMF的限定条件，r₁(t)＝a(t)-c₁(t)；若不满足条件，则将c₁(t)作为新的信号，重复步骤A1～A3，直到c₁(t)满足IMF条件为止；

步骤A4中，IMF的限定条件为：

条件1：imfi(t)必须满足信号局部极值数目之和与过零点数目相等或最多相差1个；

条件2：imfi(t)必须满足信号的局部极值所定义下包络线的平均值为零；

步骤A5：若r₁是单调函数，分解结束，若r₁不是单调函数，r₁作为原函数重复步骤A1～A5；分解结束后，信号被分解为n个IMF分量和1个趋势项r_n,表示为：

S202.计算第i个IMF分量与预降噪后的微震信号X′之间的相关系数：

其中，Cov为协方差函数，Var为方差函数；

第i个IMF分量相当于是由很多个点组成,把每个时间点对应的值带入公式计算得到该IMF分量与信号X′的相关系数；

S203.在i＝1,2,…,M时，重复执行步骤S202，得到每一个IMF分量与信号X′，其中与信号X′相关系数最大的IMF分量即为分界IMF分量。

进一步地，所述步骤S3包括：

S301.设imf₁,imf₂,...,imf_M中，第j个分量为分界IMF分量；j为1～M之间的整数；对于imf₁,imf₂...,imf_j中的每一个IMF分量进行小波包变换降噪，得到imf₁′,imf₂′...,imf_j′；

小波包降噪过程中，分解层数为3；小波基类型为db9小波

S302.将降噪后的IMF分量imf₁′,imf₂′...,imf_j′与剩余IMF分量imf_j+1,imf_j+2,...,imf_M相加，实现对降噪后的信号的重构，得到重构后的信号X″：

X″＝imf₁′+imf₂′+...+imf_j′+imf_j+1+imf_j+2,...,imf_M；

设得到的信号表示为X″＝{s(1)″,s(2)″,...,s(N)″}。

如图12所示，通过EMD分解可以得到多个imf，本实施例中是8个，如图是分解得到的其中一个imf₁，imf₁相当于是由很多个点组成,把每个时间点对应的幅值带入下面的公式可计算得到imf₁与原始信号的相关系数。用相同的道理得到其他imf与预降噪后微震信号X′的相关系数，最后可以得到各个imf与预降噪后微震信号X′的相关系数。可以看出IMF3就是分解分量(相关系数达到了最大值)，然后本文就把前两个imf继续进行小波包降噪后，再与后六个imf分量重构；

进一步地，所述微震信号降噪方法还包括：计算含噪微震信号X与降噪后的信号X″之间的信偏比，作为降噪效果的判断依据，信偏比SDR计算方式如下：

其中，s(n)为含噪微震信号X中第n个采样点的值，s(n)″为信号X″中第n个采样点的值，n＝1,2,…N。

在本申请的第一个实施例中，采用仿真信号进行去噪验证

采用模拟地震的Ricker子波进行仿真分析，表达式为：

f(t)＝(1-2π²f_p ²t²)exp-(π²f_p ²t²)

式中f(t)为振幅，fp为谱峰频率。取fp＝30Hz，采样频率为1kHz，采样点数为1000。在Ricker子波中加入信噪比为4dB的高斯白噪声。实施例中原始信号和加噪信号的波形图如图2所示，其中图2(a)为原始波形图，图2(b)为加噪波形图；

图3为实施例中原始信号和加噪信号的频谱图，其中图3(a)为原始信号频谱图，可以看出主频为30HZ，图3(b)为加噪信号频谱图，可以看出存在大量高频噪声。如图4所示，小波包变换降噪有效抑制了高频噪声，其中，图4(a)是小波包变换降噪后的波形图，如图4(a)所示，该法虽能使波形更加平滑，但仍存在毛刺现象；图4(b)是小波包降噪后的频谱图，如图4(b)所示，该法虽一定程度上抑制了高频噪声，但部分高频分量仍存在。图5(a)是EMD降噪波形图，图5(b)是EMD降噪频谱图；

本发明去噪后的波形图与频谱图和三种方法降噪后的评价指标值分别如图6所示、表1所示，其中图6(a)是本发明降噪波形图，图6(b)是本发明降噪的频谱图；

表1降噪效果对比分析

从图4、6的去噪效果图反映出，小波包变换去噪后的信号存在毛刺现象和残留部分高频噪声，而本发明去噪后的信号更接近原始信号，残留的高频噪声更少。

表1展示了更为详细的结果：与其他两种方法相比，本发明降噪后的信号，信偏比最低，信噪比最高，均方根误差最小，这反映出了本发明的实用性和更加优异的去噪性能，也说明信偏比与其他的降噪指标有良好的对应性，给评价降噪效果提供一种新思路。

在本申请的第二个实施例中，采用微震信号进行去噪验证

实际数据来源于一大型水电工程地下洞室微震监测，从监测数据中随机抽取含噪微震信号进行分析，应用本文发明对微震信号波形进行降噪处理，其降噪前后的波形和频谱图如图7(a)～(d)所示，并计算3种不同方法下微震信号降噪前后的信偏比SDR，计算结果见表2。从图7(a)可以看出，原始微震信号含有大量的随机噪声，这对初至时刻拾取、频率特征分析、震源机制解释等后续研究工作产生极大干扰。图8为实施例中微震信号EMD降噪波形图和频谱图；其中，图8(a)是微震信号EMD降噪波形图，其中，图8(b)是微震信号EMD降噪频谱图；图9为实施例中微震信号小波包降噪波形图和频谱图，其中，图9(a)是微震信号小波包降噪波形图，图9(b)是微震信号小波包降噪频谱图；

如图10展示了降噪前后微震信号的时频谱图，图10(a)是原始微震信号时频谱图，图10(b)是本发明降噪的微震信号时频谱图，通过本发明方法降噪处理图9(b)后，微震信号波形和频谱更加清晰，优势频谱更加明显。从表2可以看出，WPD-EMD-WPD方法降噪后的信号有较低的信偏比，说明降噪后的信号与不含噪声的微震信号偏离程度较小。

表2微震信号降噪效果分析

综上，本发明提高了对微震信号的去噪能力，克服了小波包变换缺乏自适应性，避免了传统EMD分解产生的模态混叠现象。通过含噪Ricker子波去噪以及实际的微震信号去噪两种实验表明，本发明方法能够进一步提高去噪信号的信噪比，降低均方根误差，适合低信噪比随机非平稳的微震信号降噪，并且提出的信噪比与其他降噪效果指标有良好的对应性，为评价降噪效果提供一种新思路。

上述说明示出并描述了本发明的一个优选实施例，但如前所述，应当理解本发明并非局限于本文所披露的形式，不应看作是对其他实施例的排除，而可用于各种其他组合、修改和环境，并能够在本文所述发明构想范围内，通过上述教导或相关领域的技术或知识进行改动。而本领域人员所进行的改动和变化不脱离本发明的精神和范围，则都应在本发明所附权利要求的保护范围内。

Claims (5)

1.一种基于WPD-EMD-WPD的微震信号降噪方法，其特征在于：包括以下步骤：

步骤S1：获取含噪微震信号X＝{s(1),s(2),...,s(N)}，采用小波包变换的降噪方法对含噪微震信号进行预降噪,得到预降噪后的微震信号X′＝{s(1)′,s(2)′,...,s(N)′}，其中N为采样点个数；

步骤S2：将预降噪后的微震信号X′＝{s(1)′,s(2)′,...,s(N)′}进行经验模态分解，得到一组频率由高到低的IMF分量：imf₁,imf₂,...,imf_M；通过各个IMF分量与预降噪后微震信号X′的相关系数找出噪音与有效信号的分界IMF分量；

步骤S3：对分界IMF分量之前的IMF分量进行小波包变换降噪,降噪后的IMF分量与剩余IMF分量重构得到降噪后的信号X″，设得到的信号X″＝{s(1)″,s(2)″,...,s(N)″}。

2.根据权利要求1所述的一种基于WPD-EMD-WPD的微震信号降噪方法，其特征在于：所述步骤S1中，采用小波包变换的降噪方法对含噪微震信号进行预降噪时，选择能反映原信号特点的db9小波和3层的分解层数；

所述步骤S3中，对分界IMF分量之前的IMF分量进行小波包变换降噪时，选择db9小波和3层的分解层数。

3.根据权利要求1所述的一种基于WPD-EMD-WPD的微震信号降噪方法，其特征在于：所述步骤S2包括以下子步骤：

S201.将预降噪后的微震信号X′＝{s(1)′,s(2)′,...,s(N)′}进行经验模态分解，得到一组频率由高到低的IMF分量：imf₁,imf₂,...,imf_M；

S202.计算第i个IMF分量与预降噪后的微震信号X′之间的相关系数：

其中，Cov为协方差函数，Var为方差函数；

第i个IMF分量相当于是由很多个点组成,把每个时间点对应的值带入公式计算得到该IMF分量与信号X′的相关系数；

S203.在i＝1,2,…,M时，重复执行步骤S202，得到每一个IMF分量与信号X′，其中与信号X′相关系数最大的IMF分量即为分界IMF分量。

4.根据权利要求1所述的一种基于WPD-EMD-WPD的微震信号降噪方法，其特征在于：所述步骤S3包括：

S301.设imf₁,imf₂,...,imf_M中，第j个分量为分界IMF分量；j为1～M之间的整数；对于imf₁,imf₂...,imf_j中的每一个IMF分量进行小波包变换降噪，得到imf₁′,imf₂′...,imf_j′；

S302.将降噪后的IMF分量imf₁′,imf₂′...,imf_j′与剩余IMF分量imf_j+1,imf_j+2,...,imf_M相加，实现对降噪后的信号的重构，得到重构后的信号X″：

X″＝imf₁′+imf₂′+...+imf_j′+imf_j+1+imf_j+2,...,imf_M；

设得到的信号表示为X″＝{s(1)″,s(2)″,...,s(N)″}。

5.根据权利要求1所述的一种基于WPD-EMD-WPD的微震信号降噪方法，其特征在于：所述微震信号降噪方法还包括：计算含噪微震信号X与降噪后的信号X″之间的信偏比，作为降噪效果的判断依据，信偏比SDR计算方式如下：

其中，s(n)为含噪微震信号X中第n个采样点的值，s(n)″为信号X″中第n个采样点的值，n＝1,2,…N。

Images