CN110659621A - 一种基于变分模态分解和排列熵的联合降噪方法 - Google Patents

Abstract

本发明公开了一种基于变分模态分解和排列熵的联合降噪方法，属于信号处理技术领域。该联合降噪方法包括如下步骤：首先，读取含噪信号x(t)；计算x(t)的尺度指数a值，根据a值确定变分模态分解(VMD)分解个数K；其次，对含噪信号x(t)进行K层VMD分解，得到一系列变分模态分量u_k；然后，计算各变分模态分量u_k的排列熵值，进一步确定K的值；最后，剔除噪声u_k分量，对剩余的u_k分量进行重构，得到降噪滤波后的信号。本发明方法具有自适应确定分解个数，有效去除噪声分量，鲁棒性和实时性强的优点，能对信号进行有效降噪滤波处理，对非平稳信号降噪方法具有很好的技术价值和应用前景。

Description

一种基于变分模态分解和排列熵的联合降噪方法

技术领域

本发明属于信号处理技术领域，具体涉及一种基于变分模态分解和排列熵的联合降噪方法。

背景技术

由众多传感器设备获取到的监测信号，具有强噪声、非平稳性等特征，对监测信号进行降噪处理是对信号进行数据分析的必要前序步骤，需将真实、有效的信号从噪声中分离出来。

目前常用的信号分解方法有经验模态分解(EMD)、集成经验模态分解(EEMD)、局部均值分解(LMD)和小波变换等。EMD、EEMD或者LMD属于递归“筛选”模态，存在端点效应和模态混叠现象，且受采样频率影响，分解误差较大。小波变换需选取合适的基波函数和阈值才能达到较好的降噪效果，且不同的小波基分析同一个信号会产生不同的分解结果。

发明内容

本发明的目的在于提供一种基于变分模态分解和排列熵的联合降噪方法，其有效地解决了模态混叠问题，具有坚实的理论基础，有很好的噪声鲁棒性，具有良好的降噪效果。

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

一种基于变分模态分解和排列熵的联合降噪方法，依次包括如下步骤：

步骤1：读取含噪信号x(t)，其中，t＝1,2，…，N，N为信号的采样点个数；

步骤2：采用DFA去趋势波动分析方法，计算所述的含噪信号x(t)的尺度指数a值，根据所述的尺度指数a值来确定参数J的值以及VMD变分模态分解的分解个数K的范围，具体包括以下子步骤：

步骤2.1：计算x(t)，t＝1,2，…，N的累积时间序列y(j)，j＝1,2，…，N

式(1)中：

步骤2.2：将累积时间序列y(j)分成N_n＝[N/n]个非重叠时间序列段y_n(j)，每个时间序列段长度为n，将n称为尺度指标；

步骤2.3：利用2阶最小二乘法对每个时间序列段y_n(j)进行拟合，时间序列段y_n(j)描述为式(2)：

y_n(j)＝a_n·j²+b_n·j+c_n (2)；

步骤2.4：波动函数F(n)描述为式(3)：

步骤2.5：取尺度指标n为4≤n≤16，按照固定步长1递增，重复步骤2.2～步骤2.4，得到波动函数F(n)随尺度指标n变化的曲线，如果F(n)与n之间存在式(4)的幂律关系：

F(n)∝n^a (4)；

式(4)中，a即为要求解的尺度指数；

步骤2.6：若a≤0.75，则J＝1；若0.75＜a≤1.0，则J＝2；若1.0＜a≤1.5，则J＝3；若a＞1.5时，则J＝4；K的取值范围为s_J～e_J；其中s₁＝2，e₁＝4；s₂＝5，e₂＝8；s₃＝9，e₃＝14；s₄＝e₄＝15；

步骤3：对含噪信号x(t)进行K层VMD分解，得到一系列变分模态分量u_k；

步骤4：计算各变分模态分量u_k的排列熵值H_P(k)，k＝1，2，…，K，若H_P(1：K)＜θ的个数大于J，则继续执行步骤5，否则K＝K+1，则跳转执行步骤3；

步骤5：剔除噪声u_k分量，对剩余的u_k分量进行重构，得到降噪滤波后的信号

作为本发明的一个优选方案，在步骤3中对含噪信号x(t)采用VMD方法分解为一系列变分模态分量u_k＝{x_k(t)，t＝1,2，…，N}，k＝1,2，…，K，其中二次惩罚因子取值为2000。

作为本发明的另一个优选方案，在步骤4中，计算各变分模态分量u_k的排列熵值，具体包括如下子步骤：

步骤4.1：对变分模态分量u_k的时间序列x_k(t)，t＝1,2，…，N进行相空间重构，得到矩阵X_k为：

式(5)中，m为嵌入维数，τ为延迟时间，M为重构空间中重构分量的个数，且M＝N-(m-1)τ；

步骤4.2：矩阵X_k的每一行都是一个重构分量，由每一个重构分量可以得到一组符号序列：

P(l)＝[j₁，j₂，…，j_m] (6)；

式中l＝1，2，...，q，且q≤m！，m个不同的符号序列[j₁，j₂，...，j_m]共有m！种排列；

步骤4.3：计算每一种符号序列出现的概率P₁，P₂，...，P_m，按照Shannon信息熵的形式，时间序列x_k(t)的q种不同符号序列的排列熵可定义为：

对H_P(k)(m)进行归一化处理，即：

步骤4.4：获得各变分模态分量u_k的排列熵值为H_P(k)，k＝1,2，…，K；

步骤4.5：若H_P(1：K)＜θ的个数大于J，θ＝0.75，则继续执行步骤5，否则K＝K+1，则跳转执行步骤3。

进一步的，在步骤5中，具体包括如下子步骤：

步骤5.1：应用式(9)计算排列熵值最大的变分模态分量下标值：

s＝{s|H_P(s)＝max(H_P(k))}，k＝1，2，…，K (9)；

步骤5.2：降噪后的信号可以表示为式(10)：

本发明原理如下：

为实现非平稳监测信号的有效降噪滤波，本发明针对非平稳信号的随机性、高噪声等特点。首先，读取含噪信号x(t)；计算x(t)的尺度指数a值，根据a值确定变分模态分解(VMD)分解个数K；其次，对含噪信号x(t)进行K层VMD分解，得到一系列变分模态分量u_k；然后，计算各变分模态分量u_k的排列熵值，进一步确定K值；最后，剔除噪声u_k分量，对剩余的u_k分量进行重构，得到降噪滤波后的信号。

VMD算法在分解信号时需要预先设定分解所得变分模态分量的个数，分解个数K设置的不同，最终降噪效果也将不同，本发明采用DFA方法分析含噪信号的长程相关性并确定分解个数K值；同时，排列熵具有较高的鲁棒性和快速简便的算法，可用于区分有效监测信号与噪声。

与现有技术相比，本发明带来了以下有益技术效果：

采用DFA方法确定VMD分解个数，根据含噪信号的长程相关特性自适应确定K值；以各变分模态分量的排列熵值为特征筛选噪声分量。本方法不但可以滤除高频噪声，对有效信号间的中频噪声也可以很好的滤除，因此，本发明所提出的滤波降噪方法，能对实际信号进行自适应分解，最大程度上对信号进行滤波，充分保留了信号随机性和非平稳的特征，该算法具有算法简单、自适应性和实时性强的特点，在对非平稳信号降噪方法具有很好的技术价值和应用前景。

附图说明

下面结合附图对本发明做进一步说明：

图1为本发明一种基于变分模态分解和排列熵的联合降噪方法的流程图；

图2为典型非平稳颠簸信号时序序列的示意图；

图3为典型非平稳含噪颠簸信号时序序列x(t)的示意图；

图4中(4a)、(4b)、(4c)、(4d)、(4e)、(4f)、(4g)为含噪信号x(t)经VMD分解后得到的7个变分模态分量示意图；

图5为降噪后的颠簸信号时序序列示意图；

图6中(6a)、(6b)、(6c)分别为原始颠簸信号、含噪颠簸信号、降噪后颠簸信号的对比图。

具体实施方式

本发明提出了一种基于变分模态分解和排列熵的联合降噪方法，为了使本发明的优点、技术方案更加清楚、明确，下面结合具体实施例对本发明做详细说明。

如图1所示，一种基于变分模态分解和排列熵的联合降噪方法，具体包括如下步骤：

步骤1：读取含噪信号x(t)，其中，t＝1，2，…，N，N为信号的采样点个数；

步骤2：采用去趋势波动分析(DFA)方法，计算x(t)的尺度指数a值，根据a值确定参数J的值以及变分模态分解(VMD)分解个数K的范围，具体步骤如下：

步骤2.1：计算x(t)，t＝1，2，…，N的累积时间序列y(j)，j＝1，2，…，N

式(1)中：

步骤2.2：将累积时间序列y(j)分成N_n＝[N/n]个非重叠时间序列段y_n(j)，每个时间序列段长度为n，将n称为尺度指标；

步骤2.3：利用2阶最小二乘法对每个时间序列段y_n(j)进行拟合，时间序列段y_n(j)描述为式(2)：

y_n(j)＝a_n·j²+b_n·j+c_n (2)；

步骤2.4：波动函数F(n)描述为式(3)：

步骤2.5：取尺度指标n为4≤n≤16，按照固定步长1递增，重复步骤2.2～步骤2.4，得到波动函数F(n)随尺度指标n变化的曲线，如果F(n)与n之间存在式(4)的幂律关系：

F(n)∝n^a (4)；

式(4)中，a即为要求解的尺度指数；

步骤2.6：若a≤0.75，则J＝1；若0.75＜a≤1.0，则J＝2；若1.0＜a≤1.5，则J＝3；若a＞1.5时，则J＝4；K的取值范围为s_J～e_J；其中s₁＝2，e₁＝4；s₂＝5，e₂＝8；s₃＝9，e₃＝14；s₄＝e₄＝15；

步骤3：对含噪信号x(t)进行K层VMD分解，得到一系列变分模态分量u_k；

步骤4：计算各变分模态分量u_k的排列熵值H_P(k)，k＝1，2，…，K，若H_P(1：_K)≥θ的个数大于J，则继续执行步骤5，否则K＝K+1，则跳转执行步骤3，具体步骤如下：

步骤4.1：对变分模态分量u_k的时间序列x_k(t)，t＝1，2，…，N进行相空间重构，得到矩阵X_k为：

式中，m为嵌入维数，τ为延迟时间，M为重构空间中重构分量的个数，且M＝N-(m-1)τ；

步骤4.2：矩阵X_k的每一行都是一个重构分量，由每一个重构分量可以得到一组符号序列：

P(l)＝[j₁，j₂，…，j_m] (6)；

式中l＝1，2，...，q，且q≤m！，m个不同的符号序列[j₁，j₂，...，j_m]共有m！种排列；

步骤4.3：计算每一种符号序列出现的概率P₁，P₂，...，P_m，按照Shannon信息熵的形式，时间序列x_k(t)的q种不同符号序列的排列熵可定义为：

对H_P(k)(m)进行归一化处理，即：

步骤4.4：获得各变分模态分量u_k的排列熵值为H_P(k)，k＝1，2，…，K；

步骤4.5：若H_P(1：K)＜θ的个数大于J，θ＝0.75，则继续执行步骤5，否则K＝K+1，则跳转执行步骤3；

步骤5：剔除噪声u_k分量，对剩余的u_k分量进行重构，得到降噪滤波后的信号

具体步骤如下：

步骤5.1：应用式(9)计算排列熵值最大的变分模态分量下标值：

s＝{s|H_P(s)＝max(H_P(k))}，k＝1，2，…，K (9)；

步骤5.2：降噪后的信号可以表示为：

图2所示为典型非平稳颠簸信号的时序序列，加噪后的颠簸信号(SNR＝10dB)如图3所示。

步骤1：获取以时间(ms)为横轴，振幅为纵轴，加噪后的颠簸信号(SNR＝10dB)的时间序列x(t)，t＝1，2，...，2048；

步骤2：采用去趋势波动分析(DFA)方法，计算x(t)的尺度指数a＝0.91，则J＝2，确定变分模态分解(VMD)分解个数K＝5～8；

步骤3：对含噪信号x(t)进行5层VMD分解，得到5个变分模态分量u₁，u₂，u₃，u₄，u₅；

步骤4：计算各变分模态分量的排列熵值，如表1所示：

表1各变分模态分量的排列熵值

u<sub>1</sub>	u<sub>2</sub>	u<sub>3</sub>	u<sub>4</sub>	u<sub>5</sub>
					0.4001	0.6596	0.7582	0.8238	0.7834

因HP_(1：5)＜0.75的个数为2，不大于J，则K加1，K＝6，跳转执行步骤3，经验证，K＝6依然不满足该判断条件；

直到K＝7时，如图4所示，(4a)、(4b)、(4c)、(4d)、(4e)、(4f)、(4g)分别是对含噪信号x(t)经VMD分解后得到的7个变分模态分量，计算各变分模态分量的排列熵值，如表2所示：

表2各变分模态分量的排列熵值

u<sub>1</sub>	u<sub>2</sub>	u<sub>3</sub>	u<sub>4</sub>	u<sub>5</sub>	u<sub>6</sub>	u<sub>7</sub>
							0.4258	0.4098	0.6184	0.7501	0.7703	0.8117	0.6787

由表2可知，当K＝7时，片_P(1：5)＜0.75的个数为3，大于J，则继续执行步骤5；

步骤5：根据式(9)计算得到排列熵值最大的变分模态分量下标值为s＝6；

降噪后的颠簸信号可以表示为：

降噪滤波后的微震信号如图5所示，降噪后的微震信号的信噪比提升为12.64dB。

非平稳信号具有随机性、强噪声的特点，其频率分布较为分散，有效信号主要集中于中、低频部分，高频部分往往为随机噪声，因此可通过对由VMD分解后的各变分模态分量进行排列熵计算，将排列熵值较高的高频部分的u_k分量进行剔除，将排列熵值较低的表征有效信号的u_k分量重构，即可实现含噪非平稳信号的消噪滤波。

本发明的优点：采用DFA方法分析含噪信号的长程相关性，由尺度指数a值自适应确定分解个数K值；根据各分量的排列熵值，在充分保留了微震信号本身的非平稳特征的基础上，有效减少和消除信号中混杂的噪声成份。

如图6所示，通过原始颠簸信号(6a)、含噪颠簸信号(6b)以及降噪后颠簸信号(6c)的对比可以看出，应用VMD和排列熵进行非平稳信号的联合降噪滤波完整地保留了原信号的波形特征，同时也较好地保留了原信号的尖峰与突变部分，得到了较好的降噪滤波效果。

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

Claims (4)

1.一种基于变分模态分解和排列熵的联合降噪方法，其特征在于，依次包括如下步骤：

步骤1：读取含噪信号x(t)，其中，t＝1，2，…，N，N为信号的采样点个数；

步骤2：采用DFA去趋势波动分析方法，计算所述的含噪信号x(t)的尺度指数a值，根据所述的尺度指数a值来确定参数J的值以及VMD变分模态分解的分解个数K的范围，具体包括以下子步骤：

步骤2.1：计算x(t)，t＝1，2，…，N的累积时间序列y(j)，j＝1，2，…，N

式(1)中：

步骤2.2：将累积时间序列y(j)分成N_n＝[N/n]个非重叠时间序列段y_n(j)，每个时间序列段长度为n，将n称为尺度指标；

步骤2.3：利用2阶最小二乘法对每个时间序列段y_n(j)进行拟合，时间序列段y_n(j)描述为式(2)：

y_n(j)＝a_n·j²+b_n·j+c_n (2)；

步骤2.4：波动函数F(n)描述为式(3)：

步骤2.5：取尺度指标n为4≤n≤16，按照固定步长1递增，重复步骤2.2～步骤2.4，得到波动函数F(n)随尺度指标n变化的曲线，如果F(n)与n之间存在式(4)的幂律关系：

F(n)∝n^a (4)；

式(4)中，a即为要求解的尺度指数；

步骤2.6：若a≤0.75，则J＝1；若0.75＜a≤1.0，则J＝2；若1.0＜a≤1.5，则J＝3；若a＞1.5时，则J＝4；K的取值范围为s_J～e_J；其中s₁＝2，e₁＝4；s₂＝5，e₂＝8；s₃＝9，e₃＝14；s₄＝e₄＝15；

步骤3：对含噪信号x(t)进行K层VMD分解，得到一系列变分模态分量u_k；

步骤4：计算各变分模态分量u_k的排列熵值H_P(k)，k＝1,2，…，K，若H_P(1：K)＜θ的个数大于J，则继续执行步骤5，否则K＝K+1，则跳转执行步骤3；

步骤5：剔除噪声u_k分量，对剩余的u_k分量进行重构，得到降噪滤波后的信号

2.根据权利要求1所述的一种基于变分模态分解和排列熵的联合降噪方法，其特征在于，在步骤3中对含噪信号x(t)采用VMD方法分解为一系列变分模态分量u_k＝{x_k(t)，t＝1，2，…，N}，k＝1，2，…，K，其中二次惩罚因子取值为2000。

3.根据权利要求1所述的一种基于变分模态分解和排列熵的联合降噪方法，其特征在于，在步骤4中，计算各变分模态分量u_k的排列熵值，具体包括如下子步骤：

步骤4.1：对变分模态分量u_k的时间序列x_k(t)，t＝1,2，…，N进行相空间重构，得到矩阵X_k为：

式(5)中，m为嵌入维数，τ为延迟时间，M为重构空间中重构分量的个数，且M＝N-(m-1)τ；

步骤4.2：矩阵X_k的每一行都是一个重构分量，由每一个重构分量可以得到一组符号序列：

P(l)＝[j₁，j₂，…，j_m] (6)；

式中l＝1，2，...，q，且q≤m！，m个不同的符号序列[j₁，j₂，...，j_m]共有m！种排列；

步骤4.3：计算每一种符号序列出现的概率P₁，P₂，...，P_m，按照Shannon信息熵的形式，时间序列x_k(t)的q种不同符号序列的排列熵可定义为：

对H_P(k)(m)进行归一化处理，即：

步骤4.4：获得各变分模态分量u_k的排列熵值为H_P(k)，k＝1,2，…，K；

步骤4.5：若H_P(1：K)＜θ的个数大于J，θ＝0.75，则继续执行步骤5，否则K＝K+1，则跳转执行步骤3。

4.根据权利要求1所述的一种基于变分模态分解和排列熵的联合降噪方法，其特征在于，在步骤5中，具体包括如下子步骤：

步骤5.1：应用式(9)计算排列熵值最大的变分模态分量下标值：

s＝{s|H_P(s)＝max(HP_(k))}，k＝1,2，…，K (9)；

步骤5.2：降噪后的信号可以表示为式(10)：

Images