CN107515424A - 一种基于vmd与小波包的微震信号降噪滤波方法 - Google Patents

Abstract

本发明公开了一种基于VMD与小波包的微震信号降噪滤波方法，属于信号处理技术领域。本发明采取VMD与小波包相结合的方式，借助VMD分解方法的自适应性以及该算法本身具有强大的数学理论基础、抑制高频噪声等的特点，以及小波包具有将频带进行多层次划分，对多分辨分析中没有细分的高频部分进一步分解，能够根据被分析信号的特征，自适应地选择频带的特性，进而提高了时频分辨率；此方法能够在保留微震信号随机性、非平稳、突发瞬态特征的基础上，对微震信号进行滤波，该算法简单易行、效果较为理想，能对矿山含噪微震信号进行有效降噪滤波，具有很好的技术价值和应用前景。

Description

一种基于VMD与小波包的微震信号降噪滤波方法

技术领域

本发明属于信号处理技术领域，具体涉及一种基于VMD与小波包的微震信号降噪滤波方法。

背景技术

岩石发生破裂时诱发微震，形成微震数据，而煤矿井下噪声污染严重，因此微震数据中包含了大量外部噪声，需将微震有效信号从噪声中分离出来。

目前常用的岩石破裂微震信号的降噪滤波方法有经验模态分解(EMD)、集成经验模态分解(EEMD)，小波分析等，这些方法运算速度慢、抗噪性能差、误判率高、拾取精度低、算法实时性不强。如EMD在分解过程中会产生模态混叠现象，即分解得到的一个或多个IMF中包含差异极大的特征时间尺度，信号和噪声混叠在一个或多个IMF中，很难达到有效地降噪滤波效果。

发明内容

针对现有技术中存在的上述技术问题，本发明提出了一种基于VMD与小波包的微震信号降噪滤波方法，设计合理，克服了现有技术的不足，具有良好的效果。

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

一种基于VMD与小波包的微震信号降噪滤波方法，包括以下步骤：

步骤1：读取含噪微震信号的监测数据时序序列X(t)，其中，t＝1,2,...,T；

步骤2：对含噪微震信号X进行VMD分解，得到一系列变分模态分量；

步骤3：对各变分模态分量进行频谱分析，根据模态主频范围以及频谱方差选取含有有用信号的模态；其中频谱方差求解的具体步骤如下：

步骤3.1：分别对各个模态进行频谱分析，计算各模态频谱方差D_k(k＝1,2,3...6)，具体包括如下步骤：

步骤3.1.1：计算各个频谱分量的均值：

其中，N为每个模态的长度，S(ω)(ω＝0,1,2...N-1)为各频谱分量的值；

步骤3.1.2：计算频谱方差值：

步骤3.2：设定硬阈值λ来判定各个模态是否为有用信号；

λ＝max(D_k)/ε (14)；

其中，max(D_k)为模态中最大频谱方差，ε为设置的参数；若模态主频范围小于50HZ且D_k≥λ则判定为有用信号，否则为噪声信号；

步骤4：对选取的含有有用信号的模态分别进行小波包降噪，然后对降噪后的模态进行重构，得到VMD与小波包降噪后的微震信号。

优选地，在步骤2中，对含噪微震信号进行VMD分解，变分约束问题是寻求k个模态函数u_k(t)(k＝1,2,3...6)，要求分解后的各个模态分量的估计带宽之和最小，且各模态之和等于含噪微震信号X，具体的构造为：

通过Hilbert变换，得到每个模态函数u_k(t)的解析信号，目的是得到它的单边频谱:

其中，δ(t)为狄拉克(Dirac)函数，*表示卷积，j²＝-1；

加入将各模态的频谱调制到相应的基频带：

其中，为预估中心频率；

求取解调信号梯度的二范数，估计各模态带宽，则变分约束问题为：

其中，X为含噪微震信号，{u_k}＝{u₁,u₂…u_k}为分解得到的k个变分模态分量，{ω_k}＝{ω₁,ω₂,...ω_k}为k个变分模态分量的中心频率，符号为所有变分模态分量之和；

对变分约束问题求解，引入增广拉格朗日将变分约束问题变为变分非约束问题，其表达式如式(4)所示：

其中，α为二次惩罚因子，λ(t)为拉格朗日乘法算子；

为寻求增广拉格朗日表达式的‘鞍点’来解决式(3)的最小值问题，采用交替方向乘子法优化算法，通过交替更新u_k ⁿ⁺¹、和λⁿ⁺¹来寻求增广拉格朗日表达式的‘鞍点’；其中：

其中，ω_k等同于等同于∑_i≠ku_i(t)ⁿ⁺¹；

利用Parseval/Plancherel傅里叶等距变换，将式(5)转变到频域：

其中，为含噪微震信号X(t)的傅里叶变换，j²＝-1；

将式(6)第一项的ω用ω-ω_k代替并写成非负频率区间积分形式:

将式(7)中的第一项置零得到二次优化问题为:

同理，对于的最小值问题，将中心频率更新问题转换到频域，解得中心频率为:

式中，为当前余项的维纳滤波；为模态功率谱的重心；对进行傅里叶逆变换，那么实部{u_k(t)}即为所求；

求解变分问题的具体步骤如下：

步骤2.1：定义变分模态分量个数K值与惩罚因子α的值；

步骤2.2：初始化n＝0；

步骤2.3：令n＝n+1，执行整个循环；

步骤2.4：执行内层第一个循环，根据式(8)更新u_k；

步骤2.5：令k＝k+1，重复步骤2.4，直到k＝K，结束内层第一个循环；

步骤2.6：执行内层第二个循环，根据式(9)更新ω_k；

步骤2.7：令k＝k+1，重复步骤2.6，直到k＝K，结束内层第二个循环；

步骤2.8：执行外层循环，根据式(10)更新λ；

其中，τ为拉格朗日乘法算子λ(t)的更新步长参数；

步骤2.9：重复步骤2.3至步骤2.8，直到满足迭代停止条件如式(11)所示，结束整个循环，输出结果，得到K个变分模态分量；

其中，ε为求解精度。

优选地，在步骤4中，具体包括如下步骤：

步骤4.1：选择一个合适的小波并确定所需要分解的层次，然后对信号进行小波包分解；

步骤4.2：对于一个给定的嫡标准，计算最佳树，确定最优小波包基；

步骤4.3：选择一个恰当的阈值并对每一个小波包分解系数进行阈值量化；

步骤4.4：根据最低层的小波包分解系数和经过量化处理的系数，进行信号的小波包重构。

本发明原理如下：

为实现微震信号的有效降噪滤波，本发明针对微震信号非平稳、突发瞬态、随机性的特点，读取含噪微震信号进行VMD分解，对分解后的模态分量u_k分别进行频谱分析以及计算频谱方差，将主频范围在50HZ以下且符合频谱方差判定条件的模态视为有用信号进行小波包去噪，之后进行重构即可实现微震信号的降噪滤波。

本发明所带来的有益技术效果：

本方法采取VMD与小波包相结合的方式，借助VMD分解方法的自适应性以及该算法本身具有强大的数学理论基础、抑制高频噪声等的特点，以及小波包具有将频带进行多层次划分，对多分辨分析中没有细分的高频部分进一步分解，能够根据被分析信号的特征，自适应地选择频带的特性，进而提高了时频分辨率；此方法能够在保留微震信号随机性、非平稳、突发瞬态特征的基础上，对微震信号进行滤波，该算法简单易行、效果较为理想，能对矿山含噪微震信号进行有效降噪滤波，具有很好的技术价值和应用前景。

附图说明

图1为本发明一种基于VMD与小波包的微震信号降噪滤波方法的流程图。

图2为含噪微震信号时序序列X(t)的示意图。

图3为含噪微震信号时序序列X(t)经VMD分解后得到的变分模态分量示意图。

图4为各变分模态分量对应的频谱图。

图5为VMD与小波包降噪后的微震信号时序序列示意图。

图6为含噪微震信号与降噪后微震信号的对比图。

图7为原始微震信号频谱图与降噪后微震信号频谱图。

具体实施方式

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

一种基于VMD与小波包的微震信号降噪滤波方法，其流程如图1所示，包括以下步骤：

步骤1：读取含噪微震信号的监测数据时序序列X(t)，其中，t＝1,2,...,T；如图2所示；

步骤2：对含噪微震信号进行VMD分解，得到一系列变分模态分量；如图3所示；

对含噪微震信号进行VMD分解，变分约束问题是寻求k个模态函数u_k(t)(k＝1,2,3...6)，要求分解后的各个模态分量的估计带宽之和最小，且各模态之和等于含噪微震信号X，具体的构造为：

通过Hilbert变换，得到每个模态函数u_k(t)的解析信号，目的是得到它的单边频谱:

其中，δ(t)为狄拉克(Dirac)函数，*表示卷积，j²＝-1；

加入将各模态的频谱调制到相应的基频带：

其中，为预估中心频率；

求取解调信号梯度的二范数，估计各模态带宽，则变分约束问题为：

其中，X为含噪微震信号，{u_k}＝{u₁,u₂…u_k}为分解得到的k个变分模态分量，{ω_k}＝{ω₁,ω₂,...ω_k}为k个变分模态分量的中心频率，符号为所有变分模态分量之和；

对变分约束问题求解，引入增广拉格朗日将变分约束问题变为变分非约束问题，其表达式如式(4)所示：

其中，α为二次惩罚因子，λ(t)为拉格朗日乘法算子；

为寻求增广拉格朗日表达式的‘鞍点’来解决式(3)的最小值问题，采用交替方向乘子法优化算法，通过交替更新u_k ⁿ⁺¹、和λⁿ⁺¹来寻求增广拉格朗日表达式的‘鞍点’；其中：

其中，ω_k等同于等同于∑_i≠ku_i(t)ⁿ⁺¹；

利用Parseval/Plancherel傅里叶等距变换，将式(5)转变到频域：

其中，为含噪微震信号X(t)的傅里叶变换，j²＝-1；

将式(6)第一项的ω用ω-ω_k代替并写成非负频率区间积分形式:

将式(7)中的第一项置零得到二次优化问题为:

同理，对于的最小值问题，将中心频率更新问题转换到频域，解得中心频率为:

式中，为当前余项的维纳滤波；为模态功率谱的重心；对进行傅里叶逆变换，那么实部{u_k(t)}即为所求；

求解变分问题的具体步骤如下：

步骤2.1：定义变分模态分量个数K＝6与惩罚因子α＝5；

步骤2.2：初始化n＝0；

步骤2.3：令n＝n+1，执行整个循环；

步骤2.4：执行内层第一个循环，根据式(8)更新u_k；

步骤2.5：令k＝k+1，重复步骤2.4，直到k＝K，结束内层第一个循环；

步骤2.6：执行内层第二个循环，根据式(9)更新ω_k；

步骤2.7：令k＝k+1，重复步骤2.6，直到k＝K，结束内层第二个循环；

步骤2.8：执行外层循环，根据式(10)更新λ；

其中，τ为拉格朗日乘法算子λ(t)的更新步长参数；

步骤2.9：重复步骤2.3至步骤2.8，直到满足迭代停止条件如式(11)所示，结束整个循环，输出结果，得到K个变分模态分量；

其中，ε为求解精度。

步骤3：对各变分模态分量进行频谱分析，根据模态主频范围以及频谱方差选取含有有用信号的模态；当模态的主频范围在50HZ以下且频谱方差符合判定条件视为有用信号；其中频谱方差求解的具体步骤如下：

步骤3.1：分别对各个模态进行频谱分析，计算各模态频谱方差D_k(k＝1,2,3...6)，具体包括如下步骤：

步骤3.1.1：计算各个频谱分量的均值：

其中，N为每个模态的长度，S(ω)(ω＝0,1,2...N-1)为各频谱分量的值；

步骤3.1.2：计算频谱方差值：

步骤3.2：设定硬阈值λ来判定各个模态是否为有用信号；

λ＝max(D_k)/ε (14)；

其中，max(D_k)为模态中最大频谱方差，ε＝5；若模态主频范围小于50HZ且D_k≥λ则判定为有用信号，否则为噪声信号；

通过分析每个模态的频谱，由图4可以看出u₁、u₂、u₃的主频范围小于50HZ，u₄、u₅、u₆主频范围大于50HZ。各模态的频谱方差如表1所示：

表1

u1	u2	u3	u4	u5	u6
						0.892	0.615	0.596	0.0981	0.0793	0.0637

λ为0.1784，按照步骤3所述u₁、u₂、u₃的主频范围以及频谱方差符合判定条件为有用信号，其余模态为噪声信号。

步骤4：对选取的含有有用信号的模态分别进行小波包降噪，然后对降噪后的模态进行重构，得到VMD与小波包降噪后的微震信号，如图5所示，具体包括如下步骤：

步骤4.1：信号的小波包分解，选择sym8小波并确定所需要分解的层数为5，然后对信号进行小波包分解；

步骤4.2：确定最优小波包基，对于一个给定的shannon嫡标准，计算最佳树；

步骤4.3：小波包分解系数的阈值量化，对于每一个小波包分解系数，选择一个恰当的阈值并对系数进行阈值量化；

步骤4.4：信号的小波包重构，根据最低层的小波包分解系数和经过量化处理的系数，进行信号的小波包重构；

一个给定的信号u_k进行l(l＝1,2,...,5)层小波包分解，在该层分解中可以得到p＝2^l个子频带，若原始信号的最低频率成分为0，最高频率成分为ω_m，每个子带的频率宽度为ω_m/2^l。小波包分解系数重构，可以提取各频带范围内的信号，且总信号可以表示为：

式中：u_k(t)表示分解的某个模态，u_l,m表示第l层分解节点(l,m)上的重构信号， m＝0,1,2,...,p-1。

由于矿震信号频谱特征很明显，其频率分布较为分散，主频主要集中在0-50HZ，然后根据各模态主频范围以及频谱方差判定各模态是否为有用信号，将有用信号进行小波包去噪，然后重构，即可实现含噪微震信号的降噪滤波。

本发明通过VMD进行分解，克服了模态混叠现象，进一步结合小波包进行去噪将噪声分离干净，可在保留原信号尖峰与突变特征的基础上，较好的实现降噪滤波。

当然，上述说明并非是对本发明的限制，本发明也并不仅限于上述举例，本技术领域的技术人员在本发明的实质范围内所做出的变化、改型、添加或替换，也应属于本发明的保护范围。

Claims (3)

1.一种基于VMD与小波包的微震信号降噪滤波方法，其特征在于：包括以下步骤：

步骤1：读取含噪微震信号的监测数据时序序列X(t)，其中，t＝1,2,...,T；

步骤2：对含噪微震信号X进行VMD分解，得到一系列变分模态分量；

步骤3：对各变分模态分量进行频谱分析，根据模态主频范围以及频谱方差选取含有有用信号的模态；其中频谱方差求解的具体步骤如下：

步骤3.1：分别对各个模态进行频谱分析，计算各模态频谱方差D_k(k＝1,2,3...6)，具体包括如下步骤：

步骤3.1.1：计算各个频谱分量的均值：

其中，N为每个模态的长度，S(ω)(ω＝0,1,2...N-1)为各频谱分量的值；

步骤3.1.2：计算频谱方差值：

步骤3.2：设定硬阈值λ来判定各个模态是否为有用信号；

λ＝max(D_k)/ε (14)；

其中，max(D_k)为模态中最大频谱方差，ε为设置的参数；若模态主频范围小于50HZ且D_k≥λ则判定为有用信号，否则为噪声信号；

步骤4：对选取的含有有用信号的模态分别进行小波包降噪，然后对降噪后的模态进行重构，得到VMD与小波包降噪后的微震信号。

2.根据权利说明书1所述的基于VMD与小波包的微震信号降噪滤波方法，其特征在于：在步骤2中，对含噪微震信号进行VMD分解，变分约束问题是寻求k个模态函数u_k(t)(k＝1,2,3...6)，要求分解后的各个模态分量的估计带宽之和最小，且各模态之和等于含噪微震信号X，具体的构造为：

通过Hilbert变换，得到每个模态函数u_k(t)的解析信号，目的是得到它的单边频谱:

其中，δ(t)为狄拉克(Dirac)函数，*表示卷积，j²＝-1；

加入将各模态的频谱调制到相应的基频带：

其中，为预估中心频率；

求取解调信号梯度的二范数，估计各模态带宽，则变分约束问题为：

其中，X为含噪微震信号，{u_k}＝{u₁,u₂…u_k}为分解得到的k个变分模态分量，{ω_k}＝{ω₁,ω₂,...ω_k}为k个变分模态分量的中心频率，符号为所有变分模态分量之和；

对变分约束问题求解，引入增广拉格朗日将变分约束问题变为变分非约束问题，其表达式如式(4)所示：

其中，α为二次惩罚因子，λ(t)为拉格朗日乘法算子；

为寻求增广拉格朗日表达式的‘鞍点’来解决式(3)的最小值问题，采用交替方向乘子法优化算法，通过交替更新u_k ⁿ⁺¹、ω_k ⁿ⁺¹、和λⁿ⁺¹来寻求增广拉格朗日表达式的‘鞍点’；其中：

其中，ω_k等同于 等同于∑_i≠ku_i(t)ⁿ⁺¹；

利用Parseval/Plancherel傅里叶等距变换，将式(5)转变到频域：

其中，为含噪微震信号X(t)的傅里叶变换，j²＝-1；

将式(6)第一项的ω用ω-ω_k代替并写成非负频率区间积分形式：

将式(7)中的第一项置零得到二次优化问题为：

同理，对于的最小值问题，将中心频率更新问题转换到频域，解得中心频率为：

式中，为当前余项的维纳滤波；为模态功率谱的重心；对进行傅里叶逆变换，那么实部{u_k(t)}即为所求；

求解变分问题的具体步骤如下：

步骤2.1：定义变分模态分量个数K值与惩罚因子α的值；

步骤2.2：初始化n＝0；

步骤2.3：令n＝n+1，执行整个循环；

步骤2.4：执行内层第一个循环，根据式(8)更新u_k；

步骤2.5：令k＝k+1，重复步骤2.4，直到k＝K，结束内层第一个循环；

步骤2.6：执行内层第二个循环，根据式(9)更新ω_k；

步骤2.7：令k＝k+1，重复步骤2.6，直到k＝K，结束内层第二个循环；

步骤2.8：执行外层循环，根据式(10)更新λ；

其中，τ为拉格朗日乘法算子λ(t)的更新步长参数；

步骤2.9：重复步骤2.3至步骤2.8，直到满足迭代停止条件如式(11)所示，结束整个循环，输出结果，得到K个变分模态分量；

其中，ε为求解精度。

3.根据权利说明书1所述的基于VMD与小波包的微震信号降噪滤波方法，其特征在于：在步骤4中，具体包括如下步骤：

步骤4.1：选择一个合适的小波并确定所需要分解的层次，然后对信号进行小波包分解；

步骤4.2：对于一个给定的嫡标准，计算最佳树，确定最优小波包基；

步骤4.3：选择一个恰当的阈值并对每一个小波包分解系数进行阈值量化；

步骤4.4：根据最低层的小波包分解系数和经过量化处理的系数，进行信号的小波包重构。