CN110688964A - 一种基于稀疏分解的小波阈值与emd联合降噪方法 - Google Patents

Abstract

本发明是一种基于稀疏分解的小波阈值与EMD联合降噪方法。本发明首先利用EMD分解获取信号的固有模态函数IMF分量，然后计算各IMF分量与原始信号的相关性，根据相关性大小确定IMF噪声频段，采用小波阈值去噪方法对噪声频段进行处理，最后对处理后的信号进行重构得到去噪信号。通过模拟实验分析，证明该方法具有智能选择噪声频段的能力，是一种更适于鱼雷信号的去噪方法。本发明可以在降噪的同时还能成功地保留信号特征。在保留了EMD运算较快、降噪性能好的优点的同时，采用稀疏化小波阈值降噪改善了各保留IMF的噪声性能，从而提高噪声抑制水平。

Description

一种基于稀疏分解的小波阈值与EMD联合降噪方法

技术领域

本发明涉及联合降噪技术领域，是一种基于稀疏分解的小波阈值与EMD联合降噪方法。

背景技术

在信号处理中，噪声是要影响信号传输的主要因素，于是要从采集的信号中，尽可能分离噪声和有效信号，提高信噪比和分辨率。噪声是任何信号中最常见的干扰波，形成原因复杂多样，而且噪声的表现形式是没有规律的振动。由于它的视速度无法预测，传播方向不固定，频率分布范围更是无法确定，这给分离有效信号和随机噪声带来极大的挑战。

经验模态分解(Empirical Mode Decomposition,EMD)以自身数据特性为基础，实现自适应时频处理。信号经过EMD分解，可以被分解成数量有限的本征模态函数(IntrinsicMode Function,IMF)分量，联合Hilbet变换法对信号进行分析，该过程一方面可以确保信号分解的唯一性，另一方面又可以在时频域内得到优化良好的局部信号。

J.Yu and Z.Zhang使用LMD(局部均值分解)和EMD方法对信号进行去噪并比较两种方法的去噪性能，相关数据表明，使用EMD可获得更高的相关度，滤波后的信号具有更高的SNR。通常LMD分解信号后会出现模态混叠现象，造成分解效果不佳，而且，对于大部分信号而言，选取合适的小波基，然后采用小波变换后会获得稀疏的小波系数，就是含有较少的大幅度系数。但是，当信号带有噪声干扰时，将大大削弱小波变换之后小波系数具有的稀疏性。

发明内容

本发明为实现抑制噪声干扰和提高计算效率，本发明提供了一种基于稀疏分解的小波阈值与EMD联合降噪方法，本发明提供了以下技术方案：

一种基于稀疏分解的小波阈值与EMD联合降噪方法，包括以下步骤：

步骤1：采用EMD法对带躁信号s(t)进行分解，得到固有模态函数，根据固有模态函数确定原信号；

步骤2：通过计算每一个IMF分量与带噪信号的相关性，确定每一个IMF分量中是否包含信息数据，

步骤3：找到信号模态和噪声模态之间的分界点k；

步骤4：使用小波阈值降噪算法进行降噪处理，对原信号进行N层小波变换，得到对应的小波轮廓系数和小波系数，对小波轮廓系数和小波系数进行阈值处理，得到硬阈值和软阈值；

步骤5：对小波系数进行稀疏化处理，利用小波轮廓系数和小波系数进行信号重构，得到抑噪后的信号。

优选地，所述步骤1具体为：

步骤1.1：采用EMD法对带噪信号s(t)的全部极小值与极大值点进行分解拟合，得到带噪信号的下包络线e_l(t)及上包络线e_h(t)，计算上下包络线的平均值，通过下式表示上下包络线的平均值：

其中，m(t)为上下包络线的平均值；

步骤1.2：采用带噪信号减去上下包络线的平均值得到新一轮的带噪信号，通过下式表示新一轮的带噪信号：

imf(t)＝s(t)-m(t) (2)

其中，imf(t)为新一轮的带噪信号；

步骤1.3：由IMF判别准则筛选新一轮的带噪信号，得到的带噪信号的第一阶固有模态函数，直至得到所有的IMF分量为止，根据固有模态函数确定原信号，通过下式确定原信号：

其中，y(t)为原信号，imf_j(t)为固有模态函数，j＝1,2,…,N，N为时间尺度，r_N(t)为残余函数。

优选地，所述步骤2具体为：确定每一个IMF分量与原始信号的相关系数，根据相关系数确定确定每一个IMF分量中是否包含信息数据，通过下式确定每一个IMF分量与原始信号的相关系数：

其中，R(y(t),imf_j(t))为每一个IMF分量与原始信号的相关系数。

优选地，所述步骤3具体为：

步骤3.1将IMF按顺序对应于由高到低处在不一致频率段上的分量，即阶数高的IMF对应于信号的低频部分，阶数低的IMF对应于信号的高频部分；

步骤3.2：根据固有模态函数确定适信号模态和噪声模态之间的分界点k，通过下式确定所述分解点k：

优选地，根据所述分界点k将分解点前的IMF分量进行组合，第1层到第k－1层均为噪声，k-1层后的IMF分量含有有用信号。

优选地，所述步骤4具体为：

步骤4.1：使用小波阈值降噪算法进行降噪处理，选取小波基对信号进行分解，确定小波分解层数，对原信号进行N层小波变换，得到对应的小波轮廓系数和小波系数，通过下式表示小波轮廓系数和小波系数：

swa_j＝argmax||W_f(a,b)||₀ (6)

swd_j＝ΦW_f(a,b) (7)

其中，swa_j为小波轮廓系数，swd_j为小波系数，Φ为测量矩阵，a为伸缩因子，b为平移因子，W_f(a,b)为小波微积分；

步骤4.2：对小波轮廓系数swa_j和小波系数swd_j进行阈值处理，得到硬阈值和软阈值，通过下式表示硬阈值和软阈值：

其中，H_s为硬阈值，H_h为软阈值，λ为阈值标准值，B为信号长度，δ为噪声标准差。

优选地，所述步骤5具体为:

步骤5.1：采用OMP算法进行小波系数的稀疏化处理，选择随机矩阵Φ，所述随机矩阵满足RIP特性，N维信号在正交基作用下实现稀疏化，通过下式表示稀疏化后的N维信号：

x＝ψswa_j (11)

其中，x为稀疏化后的N维信号，ψ为正交基；

步骤5.2：初始化小波细节系数使得

β_j为稀疏系数；当

ε为残差阈值，则令

并使用迭代式进一步计算sωd_j ^(k)；否则，令

步骤5.3：利用获得的第k个IMF分量的小波细节系数sωd_j ^(k)和硬阈值进行小波逆变换，通过下式表示小波逆变换过程：

sωd_j′＝H_s(sωd_j ^(k)+μΦ^T(y-Φsωd_j ^(k)))μ＞0 (12)

其中，sωd_j′为小波逆变换的小波细节系数，y为M维测量向量，μ为迭代次数；

步骤5.3：通过线型测量矩阵φ对小波轮廓系数进行测量，得到远小于信号长度的M维测量向量，通过下式表示M维测量向量：

y＝φψswa_j＝Aswa_j (13)

其中，A为传感矩阵；

步骤5.4：采用信号重构算法对信号进行重构，通过下式表示重构过程：

其中，y′(t)为重构信号，r′_N(t)为去噪后的IMF分量。

本发明具有以下有益效果：

本发明采用基于稀疏分解的小波阈值函数与经验模态分解方法联合噪声抑制技术，是建立在EMD和小波阈值的基础上的，并通过稀疏表示，因此，可以在降噪的同时还能成功地保留信号特征。在保留了EMD运算较快、降噪性能好的优点的同时，采用稀疏化小波阈值降噪改善了各保留IMF的噪声性能，从而提高噪声抑制水平。

附图说明

图1是一种基于稀疏分解的小波阈值与EMD联合降噪方法流程图；

图2是硬阈值与软阈值函数图；

图3是小波阈值去噪流程图。

具体实施方式

以下结合具体实施例，对本发明进行了详细说明。

具体实施例一：

按照图1所示，本发明提供一种基于稀疏分解的小波阈值与EMD联合降噪方法，包括以下步骤：

步骤1：采用EMD法对带躁信号s(t)进行分解，得到固有模态函数，根据固有模态函数确定原信号；

在EMD是在假设任何信号都可用有限个本征模态函数(Intrinsic ModeFunction，IMF)分量组成的前提下，使信号分解后产生一组有着不一致时间尺度的固有模态分量，使得各个IMF都是窄带信号。同时对于每个IMF都要符合以下两个判别准则：

a)在所有数据区间中，函数的极值点和过零点的数目一定是至多相差一个或相同的关系；

b)任何一点处，全部极小值点所生成的下包络线与全部极大值点所生成的上包络线之间的平均值保持为0，即为该信号的上下包络线关于时间轴处于对称的位置。

EMD分解实质上为一种筛选的过程，而该过程有两个作用：让数据波形完美对称与除掉叠加波。对信号s(t)进行EMD分解的步骤如下：

1)得到信号s(t)的全部极小值点与极大值点，对应将它们利用拟合函数拟合为信号的下包络线e_l(t)及上包络线e_h(t)，然后求得上下包络线的平均值，即

然后用原信号s(t)减去m(t)，记为：

imf(t)＝s(t)-m(t)

2)使得imf(t)成为新一轮的原信号s(t)，重复上述的过程，由IMF判别准则筛选得到原信号的第一阶本征模态函数，用s(t)减去imf₁(t)，记为：

m₁(t)＝s(t)-imf₁(t)

若m₁(t)视为新一轮的s(t)，采用同一步骤直到得到所有的IMF分量为止。因此，原信号y(t)最终可表示为：

其中，r_N(t)为残余函数，表示信号的趋势；imf₁(t),imf₂(t),…,imf_N(t)其意义为信号中不一致时间尺度的分量，尺度按升序排列。

步骤2：通过计算每一个IMF分量与带噪信号的相关性，确定每一个IMF分量中是否包含信息数据；

由于重要的信号会在采集的含噪信号中占据较大的比重，因此可以通过计算每一个IMF分量与含噪信号的相关性来确定每一个IMF分量是否包含重要的信息。每一个IMF分量与原始信号的相关系数可以表示为：

步骤3：找到信号模态和噪声模态之间的分界点k；

根据固有模态函数确定适信号模态和噪声模态之间的分界点k，通过下式确定所述分解点k：

各个IMF按顺序对应于由高到低处在不一致频率段上的分量，即阶数高的IMF对应于信号的低频部分，阶数低的IMF对应于信号的高频部分。而针对多数带噪信号而言，信号的主要能量聚集于低频处，而高频处通常很少包含有用信息。EMD去噪时通过判定条件分解到第k层，并认为第1层到第k－1层均为噪声，k-1层后的IMF分量含有有用信号。

步骤4：使用小波阈值降噪算法进行降噪处理，对原信号进行N层小波变换，得到对应的小波轮廓系数和小波系数，对小波轮廓系数和小波系数进行阈值处理，得到硬阈值和软阈值；

小波变换原理：

设ψ(t)∈L²(R)(L²(R)表示平方可积的实数空间，即能量有限的信号空间)。当ψ(t)满足条件：

一种基于稀疏分解的小波阈值与EMD联合降噪方法称ψ(t)为一个基本小波或母小波，将母小波函数ψ(t)经伸缩和平移后，就可以得到一个小波序列：

其中a为伸缩因子，b为平移因子。

对于任意的函数f(t)∈L²(R)的连续小波变换为：

其中W_f(a,b)为f(t)的小波积分或频谱。

使用小波阈值降噪算法]进行降噪处理时，首先选取适当的小波函数对信号进行分解，其次对分解出的参数进行阈值处理，选取合适的阈值进行分析，最后利用降噪处理后的参数进行逆小波变换。其逆变换为：

小波变换的时频窗是可以由伸缩因子a和平移因子b来调节的，平移因子b，可以改变窗口在相平面时间轴上的位置，而伸缩因子b的大小不仅能影响窗口在频率轴上的位置，还能改变窗口的形状。小波变换对不同的频率在时域上的取样步长是可调节的，在低频时，小波变换的时间分辨率较低，频率分辨率较高：在高频时，小波变换的时间分辨率较高，而频率分辨率较低。

一般来说，小波降噪过程可以分为3个步骤如图3进行：

a)选择合适的小波基并确定小波分解的层数N，然后对信号进行N层小波变换，即可得到相应的小波变换系数swa_j和swd_j；

b)通过选择合适的阈值和阈值函数，对小波分解后的系数swa_j和swd_j进行阈值处理，获得软阈值和硬阈值分别为H_s和H_h；

c)利用经过小波稀疏化后的小波细节系数swd_j′进行信号重构，再通过小波逆变换，得到去噪后的信号。

小波阈值去噪的整个过程中，选择什么样的阈值或阈值函数也是一非常重要的问题，也将直接影响信号去噪后的效果。当函数在±λ处连续时选用软阈值，这样对小波系数重构有良好的平滑效果；当函数高阶可导时选用硬阈值，这样在一定程度上保存了初始信号的噪声分量。

小波变换得到的小波轮廓系数和小波系数估计值分别为：

swa_j＝argmax||W_f(a,b)||₀ swd_j＝ΦW_f(a,b)

其中，Φ为测量矩阵。

常用的阈值处理函数有：硬阈值和软阈值如图2，硬阈值函数为：

软阈值函数为：

其中，λ为阈值标准值，B为信号长度，δ为噪声标准差。

在实际工程应用中，有用的信号通常表现为低频信号或是一些比较平稳的信号，而噪声信号表现为高频信号。对于小波变换的多层分解，噪声部分通常包含在高频成分中。因此，从信号学的角度看，小波去噪是一个信号滤波的问题。尽管在很大程度上小波去噪可以看成是低通滤波，但由于在去噪后，还能成功地保留信号特征，所以在这一点上又优于传统的低通滤波器。由此可见，小波去噪实际上是特征提取和低通滤波的综合。

步骤5：对小波系数进行稀疏化处理，利用小波轮廓系数和小波系数进行信号重构，得到抑噪后的信号。

这里一种基于稀疏分解的小波阈值与EMD联合降噪方法首先引入压缩感知(Compressive Sensing，CS)算法，其意在充分利用小波系数的整体结构特性估计重要小波系数但不丢失重要信息，压缩感知包括三个核心步骤：采用OMP算法完成小波系数的稀疏化、观测矩阵设计、信号重构算法。压缩感知方法的前提条件是信号具有稀疏性或近似稀疏性，即N维信号x在正交基ψ下实现稀疏化表示，θ为K稀疏性向量，即存在K个非0元素

x＝ψswa_j

采用OMP算法完成小波系数的稀疏化，算法流程如下：

1)选取合适的随机矩阵Φ_M×N，使得其满足RIP特性；

2)初始化

β_j为稀疏系数；

3)假如

ε为残差阈值，则令

并使用迭代式进一步计算sωd_j ^(k)，重复步骤3；否则，令

转至步骤4；

4)如果n＝1，转至步骤5；否则，令n＝n-1，转至步骤2；

5)利用获得的第k个IMF分量的小波系数估计值sωd_j ^(k)进行小波逆变换

迭代式为：

sωd_j′＝H_s(sωd_j ^(k)+μΦ^T(y-Φsωd_j ^(k)))μ＞0

其中，swd_j ^(k)n＝0,1,2,3....表示第n个IMF分量的小波系数估计值，μ为迭代次数。

通过线性测量矩阵φ对信号进行测量，可以得到远小于信号长度的M(M＜＜N)维测量向量y：

y＝φψswa_j＝Aswa_j

其中，A称为传感矩阵。x在正交基ψ下的稀疏化可以描述为：

swd_j′＝argmin||y||₀

最后，可以通过信号重构算法对信号进行重构：

r′_N(t)为去噪后的IMF分量。为了由M维测量向量y准确恢复出N维原始信号，压缩感知对测量信号有比较严苛的要求，即著名的受限等容特性(Restricted IsometryProperty，RIP)。

以上所述仅是一种基于稀疏分解的小波阈值与EMD联合降噪方法的优选实施方式，一种基于稀疏分解的小波阈值与EMD联合降噪方法的保护范围并不仅局限于上述实施例，凡属于该思路下的技术方案均属于本发明的保护范围。应当指出，对于本领域的技术人员来说，在不脱离本发明原理前提下的若干改进和变化，这些改进和变化也应视为本发明的保护范围。

Claims (7)

1.一种基于稀疏分解的小波阈值与EMD联合降噪方法，其特征是：包括以下步骤：

步骤1：采用EMD法对带躁信号s(t)进行分解，得到固有模态函数，根据固有模态函数确定原信号；

步骤2：通过计算每一个IMF分量与带噪信号的相关性，确定每一个IMF分量中是否包含信息数据；

步骤3：找到信号模态和噪声模态之间的分界点k；

步骤4：使用小波阈值降噪算法进行降噪处理，对原信号进行N层小波变换，得到对应的小波轮廓系数和小波系数，对小波轮廓系数和小波系数进行阈值处理，得到硬阈值和软阈值；

步骤5：对小波系数进行稀疏化处理，利用小波轮廓系数和小波系数进行信号重构，得到抑噪后的信号。

2.根据权利要求1所述的一种基于稀疏分解的小波阈值与EMD联合降噪方法，其特征是：所述步骤1具体为：

步骤1.1：采用EMD法对带噪信号s(t)的全部极小值与极大值点进行分解拟合，得到带噪信号的下包络线e_l(t)及上包络线e_h(t)，计算上下包络线的平均值，通过下式表示上下包络线的平均值：

其中，m(t)为上下包络线的平均值；

步骤1.2：采用带噪信号减去上下包络线的平均值得到新一轮的带噪信号，通过下式表示新一轮的带噪信号：

imf(t)＝s(t)-m(t) (2)

其中，imf(t)为新一轮的带噪信号；

步骤1.3：由IMF判别准则筛选新一轮的带噪信号，得到的带噪信号的第一阶固有模态函数，直至得到所有的IMF分量为止，根据固有模态函数确定原信号，通过下式确定原信号：

其中，y(t)为原信号，imf_j(t)为固有模态函数，j＝1,2,…,N，N为时间尺度，r_N(t)为残余函数。

3.根据权利要求1所述的一种基于稀疏分解的小波阈值与EMD联合降噪方法，其特征是：所述步骤2具体为：确定每一个IMF分量与原始信号的相关系数，根据相关系数确定确定每一个IMF分量中是否包含信息数据，通过下式确定每一个IMF分量与原始信号的相关系数：

其中，R(y(t),imf_j(t))为每一个IMF分量与原始信号的相关系数。

4.根据权利要求1所述的一种基于稀疏分解的小波阈值与EMD联合降噪方法，其特征是：所述步骤3具体为：

步骤3.1 将IMF按顺序对应于由高到低处在不一致频率段上的分量，即阶数高的IMF对应于信号的低频部分，阶数低的IMF对应于信号的高频部分；

步骤3.2：根据固有模态函数确定适信号模态和噪声模态之间的分界点k，通过下式确定所述分解点k：

5.根据权利要求4所述的一种基于稀疏分解的小波阈值与EMD联合降噪方法，其特征是：根据所述分界点k将分解点前的IMF分量进行组合，第1层到第k－1层均为噪声，k-1层后的IMF分量含有有用信号。

6.根据权利要求1所述的一种基于稀疏分解的小波阈值与EMD联合降噪方法，其特征是：所述步骤4具体为：

步骤4.1：使用小波阈值降噪算法进行降噪处理，选取小波基对信号进行分解，确定小波分解层数，对原信号进行N层小波变换，得到对应的小波轮廓系数和小波系数，通过下式表示小波轮廓系数和小波系数：

swa_j＝argmax||W_f(a,b)||₀ (6)

swd_j＝ΦW_f(a,b) (7)

其中，swa_j为小波轮廓系数，swd_j为小波系数，Φ为测量矩阵，a为伸缩因子，b为平移因子，W_f(a,b)为小波微积分；

步骤4.2：对小波轮廓系数swa_j和小波系数swd_j进行阈值处理，得到硬阈值和软阈值，通过下式表示硬阈值和软阈值：

其中，H_s为硬阈值，H_h为软阈值，λ为阈值标准值，B为信号长度，δ为噪声标准差。

7.根据权利要求1所述的一种基于稀疏分解的小波阈值与EMD联合降噪方法，其特征是：所述步骤5具体为:

步骤5.1：采用OMP算法进行小波系数的稀疏化处理，选择随机矩阵Φ，所述随机矩阵满足RIP特性，N维信号在正交基作用下实现稀疏化，通过下式表示稀疏化后的N维信号：

x＝ψswa_j (11)

其中，x为稀疏化后的N维信号，ψ为正交基；

步骤5.2：初始化小波细节系数使得

β_j为稀疏系数；当

ε为残差阈值，则令

并使用迭代式进一步计算sωd_j ^(k)；否则，令

步骤5.3：利用获得的第k个IMF分量的小波细节系数sωd_j ^(k)和硬阈值进行小波逆变换，通过下式表示小波逆变换过程：

sωd_j′＝H_s(sωd_j ^(k)+μΦ^T(y-Φsωd_j ^(k)))μ＞0 (12)

其中，sωd_j′为小波逆变换的小波细节系数，y为M维测量向量，μ为迭代次数；

步骤5.3：通过线型测量矩阵φ对小波轮廓系数进行测量，得到远小于信号长度的M维测量向量，通过下式表示M维测量向量：

y＝φψswa_j＝Aswa_j (13)

其中，A为传感矩阵；

步骤5.4：采用信号重构算法对信号进行重构，通过下式表示重构过程：

其中，y′(t)为重构信号，r′_N(t)为去噪后的IMF分量。

Images