基于同步挤压小波变换的振动信号自适应降噪方法

2018-08-02陶新民王若彤

振动与冲击 2018年14期

沈微，陶新民，高珊，常瑞，王若彤

(东北林业大学工程技术学院，哈尔滨 150040)

在旋转机械振动信号的获取过程中，由于受到测试环境、设备以及人为因素等影响，获取的振动信号往往受到随机噪声的干扰。当设备发生故障时，这些干扰会使反应故障特征的有效信息变得非常微弱。因此，对原始信号进行降噪，提高信噪比是有效提取故障特征的前提条件和必要环节。然而，由于旋转机械内部结构的复杂性，测得的振动信号具有非平稳性、非线性等特征，进而加大了降噪难度。

针对非平稳信号的降噪问题，国内外学者已展开了大量深入研究。其大致可分为基于傅里叶变换[1]、基于小波变换[2-6]以及基于经验模式分解(Empirical Mode Decomposition, EMD)[7-9]等方法。这些方法都有各自的优缺点，其中傅里叶变换方法理论成熟，但难以处理噪声抑制和信号边缘保护两者之间的矛盾，对频率重叠的信号也难以区分；小波变换方法虽然降噪效果优于传统滤波方法且具有多分辨性能，但其时频谱在频率方向会出现能量发散现象，严重影响降噪性能。基于EMD时频分析降噪方法因其具有很强的频率带选层性能，且具有自适应分解的特点，因此通常情况下其降噪性能相对于传统小波变换方法而言有较大幅度的提高。但该方法极易受信号突变的影响，导致算法分解时出现混叠现象和端点效应，进而影响降噪效果。针对以上不足，Fladndrin等[10]将噪声辅助方法运用到经验模式分解中实现抗混分解，有效抑制了混叠现象，进而提出集成经验模式分解方法。由于信号经EEMD分解得到的固有模式分量(Intrinsic Mode Functions,IMF)其物理意义更加明显，因此学者们相继提出了诸多基于EEMD分解的降噪方法。其中文献[11]考虑到噪声多出现在高频段，提出通过剔除前几阶IMF分量只保留剩余分量的方式进行信号重构实现降噪目的。然而该方法无法满足噪声出现在低频段的情况，再者若原始信号中含有有效的高频信息，这种方法也会导致有用信息的丢失。文献[12-13]提出利用每一个IMF分量与原有信号互相关值的大小来选取IMF分量，但该方法因需要原始信号参与计算，而原始信号本身又含有噪声成分，因此该方法可能会引入虚假分量影响降噪性能。文献[14]利用振动信号周期性的特点，提出依据IMF自相关函数和原始信号自相关函数的相关系数大小进行IMF分量的选择，该方法只适合于周期性明显的信号，不具有普适性。需要说明的是EEMD虽然较好地抑制了混叠现象，但该方法由于缺少严格的数学物理基础，其稳定性会受模态分量分解过程的影响，从而造成分解结果不准确。

针对该问题，Daubechies等[15]以小波变换为基础，通过对小波系数进行压缩重组的方法实现高精度、高分辨率的时频分析[16-17]。作为新的时频分析方法最初被应用于油气检测[18]、地震弱信号提取[19]、地震信号面波去除[20]、谐波信号提取[21]等领域，而在振动信号降噪领域中的研究和应用相对较少。鉴于此，本文利用SST分解得到的分量具有高精度、高分辨率，窄带以及抗混叠等特点，根据噪声分量与真实分量频率复杂度的不同以及噪声分量自相关系数的分布特征，有选择性地进行信号重构实现降噪目的。实验最后将本文方法同其他流行的方法进行比较，结果表明本文方法的降噪性能更强。

1 同步挤压小波变换及时频性能对比分析

1.1 同步挤压小波变换简介

通常情况下传统时变信号都可以表示成多个谐波信号的叠加形式，即信号f(t)可表示为：

(1)

式中：Ak(t)为第k个谐波分量的瞬时振幅，θk(t)为第k个谐波分量的瞬时相位，e(t)为噪声或误差，K为可分解分量的个数。

同步挤压小波变换是在小波变换的基础上，根据时间-尺度平面各元素绝对值的大小，对平面中的能量进行重新分配，最后通过映射公式将时间-尺度平面转换为时间-频率平面，达到提高时间和空间分辨率的目的。本文将利用同步挤压小波变换实现振动信号降噪。同步挤压小波变换算法步骤如下：

算法首先对待分析的信号f(t)进行连续小波变换(CWT)，得到的小波系数为Wf(a,b)，表达式为：

(2)

式中：a为尺度因子，b为平移因子；ψ*为共轭小波函数，根据Plancherel定理，在频率域等价变换为：

(3)

(ai,tm),ai=2i/nvΔt,i=1,2,…，Lnv

(4)

式中：nv为自定义变量，决定了尺度系数的个数。tm为采样时间间隔点，采样时间间隔为Δt。L为最大尺度，其中f(t)的长度为n=2L+1。根据式(1)，假设f(t)为通常意义上的单一谐波函数f(t)=Acos(wt)，其傅里叶变换为：

(5)

代入式(3)中，其连续小波变换表示为：

(6)

(7)

(8)

其中median为中值函数，则在中心频率wl上同步挤压小波变换值Tf(wf,b)为：

(9)

其中(Δa)i=ai-ai-1。同步挤压小波反变换为：

(10)

1.2 小波变换与同步挤压小波变换的时频性能对比

由上述分析可知传统小波变换会在瞬时频率附近出现发散现象，导致时频谱变得十分模糊。本文以信号f(t)=1.0 cos(2π×240t)+1.0 cos(2π×310t)为例进行传统小波变换，并与同步挤压小波变换进行对比分析。其中算法参数设置如下： Morlet为基本小波，nv=64，采样频率为1 000 Hz，采样点数1 000。实验结果如图1，2所示。

图1 传统小波变换时频图Fig.1 Time-frequency spectrum after wavelet transform

图2 同步挤压小波变换时频图Fig.2 Time-frequency spectrum after SST

从图1中我们发现，信号经传统小波变换后得到的时频谱能量在瞬时频率两侧出现发散现象，小波脊线较模糊；图2是同步挤压小波变换的时频分析结果，可以看出同步挤压小波变换的时频分析效果更加清晰，在瞬时频率方向更加聚焦，得到的小波脊线也更细更明显。通过对公式(6)的分析可以看出，单一频率的谐波信号经传统小波变换后得到的结果在b点位置，沿着瞬时频率w0上下震荡(ejwb)导致频谱信息发散。同步挤压小波变换正是基于这一点，对传统小波时频谱中的能量进行压缩后重排，实现时间-尺度平面到时间-频率平面的转换。最终的结果正如图1和图2所示，同步压缩小波变换时频谱上的能量更集中，频率分辨率更高。

2 基于同步挤压小波变换的自适应降噪算法

2.1 传统相关降噪方法分析

由于传统的EEMD分解方法把信号分解为频率由高到低的一系列IMF分量和一个残余项，基于EEMD的强降噪法假设混有噪声的原始信号经EEMD分解后得到的前几个高频IMF分量为噪声，将其剔除并利用剩余的IMF分量重构信号实现降噪目的。在应用此方法进行降噪时，最关键的问题就是IMF分量的选取。如果选择不当就会造成降噪效果欠佳，甚至由于原有真实IMF分量被剔除，导致信号中有用信息的丢失。为了解决上述问题，学者们又尝试利用IMF分量与原始信号相似性作为评价指标判断IMF分量的真实性，相继提出了基于互相关信息和相关系数的EEMD降噪算法。这些方法都需要与原信号进行比较分析，而原信号已被噪声污染，用它与IMF分量进行相似性的比较显然容易导致误判，因此可以说上述算法都没有从根本上实现噪声分量与真实分量的判别。

2.2 基于IMF分量信号瞬时频率稳定度降噪方法

由于噪声信号属于非线性、非平稳信号，其瞬时频率不稳定，随时间变动会出现较大幅度的波动；而谐波信号属于线性平稳信号，瞬时频率相对平稳。因此可以通过度量瞬时频率的稳定性来判定IMF分量是否为噪声信号。其中IMF分量瞬时频率的求解方法如下：

原始信号经过同步挤压小波变换后，得到同步挤压变量Tf(wl,b)，然后计算各个小波脊线所含的主频率wq*(k)，设区间：

则可求原始信号的各个谐波分量(IMF)信号fk(t)：

(11)

当各个谐波分量提取完成后，对fk(t)进行Hilbert变换：

(12)

令Zk(t)=Xk(t)+jYk(t)=a(t)ejθk(t)

其中

(13)

则求出fk(t)的瞬时频率：

(14)

为了能定量地度量各IMF分量瞬时频率的稳定性，这里采用四分位距作为衡量瞬时频率离散度(波动)的指标。首先将经过Hilbert变换后得到的各IMF分量的瞬时频率值按升序排列，然后利用处于数据集3/4位置的上四分位数(75%，Q3)减去1/4位置的下四分位数(25%，Q1)即为四分位距。四分位距度量法能够规避数据集中异常值影响极差对数据离散度的判断问题。具体公式为：

IRQ=Q3-Q1

(15)

2.3 基于IMF分量信号自相关峰度阈值降噪方法

考虑到现实中出现的噪声普遍具有白噪声特性，其自相关函数的特点十分明显，即只与自身相关，其他情况几乎为零。这与一般信号自相关函数特性不同。为了进一步说明这一特点，我们以20 Hz和30 Hz正弦叠加谐波信号为例，进行自相关函数的计算，其结果如图3所示。从图中可以看出一般信号的自相关函数虽在零点处最大，但其余点并不立即衰减为零，而是会有一个缓慢下降的过程；相反噪声信号的自相关函数只在零点处最大，其他点处几乎为零。原因是由于随机噪声各时间点取值具有一定的随机性，既表现为弱相关性；而一般信号各时间点取值具有一定关联，表现为较强的相关性。

信号的自相关函数在一定程度上反应了信号与其自身在不同时间点的相关程度，是一种时间域的统计度量方法，其具体定义如下：

Rx(t,t+τ)=E[x(t)x(t+τ)]

(16)

(17)

式中：Rx(0)表示信号与自身在同一时刻的相关程度，显然，对于任何信号而言该值都最大。

图3 随机噪声和一般信号的自相关函数分布特征对比Fig.3 Auto-correlative function distribution characteristics of random noise and normal signal

为了能定量地描述相关函数的分布特点，我们引入峰度系数作为描述相关函数分布特征的指标。峰度系数(Kurtosis)用来度量数据在中心聚集程度，即分布曲线顶端尖峭或扁平程度，值越大则越陡峭。在正态分布情况下，峰度系数值是0。正的峰度系数说明数据更集中，具有比正态分布更长的尾部；负的峰度系数说明数据不那么集中，具有比正态分布更短的尾部。

(18)

式中：N代表样本个数；σ表示数据集合的标准方差。上例中的一般谐波信号自相关函数的峰度值为4.062 2，而随机信号自相关函数的峰度值为114.038 1。可见，噪声信号的自相关函数分布极其陡峭，如图3所示。因此以自相关函数峰度值作为区分噪声IMF分量和真实IMF分量的指标在理论上是可行的。

2.4 基于SST的振动信号自适应降噪方法

本文将上述两种方法结合后提出一种基于同步挤压小波变换的振动信号自适应降噪方法，具体流程如图4所示。

图4 基于SST自相关峰度阈值和瞬时频率稳定度的降噪方法流程图Fig.4 Flow of de-noising method based on SST self autocorrelation’s kurtosis and instantaneous frequency

算法首先对原始信号进行SST分解，分解后得到不同的本征模态分量; 然后利用Hilbert 变换求解每一个本征模态分量的瞬时频率曲线并计算其四分位距值，如果该值大于给定阈值则剔除对应的本征模态分量。最后针对筛选后得到的每一个本征模态分量，求解自相关函数并计算其峰度，依据给定阈值判断是否满足随机噪声自相关函数的分布特征，将大于阈值的本征模态分量剔除。为了防止出现分解不完全的情况，本文对于残余量的判断只采用基于自相关函数峰度阈值的判定方法。由于本文方法不要求噪声只出现在高频部分，因此可以保留原始信号有用的高频信号分量，同时由于本文方法也无需同受噪声污染的原信号进行相似性判断，避免了真假IMF分量的误判问题。经多次试验本文设定自相关峰度阈值为10，瞬时频率稳定度阈值为100。

3 实验分析及对比

3.1 仿真实验降噪性能对比分析

为了验证本文降噪方法的可行性和有效性，实验中首先利用仿真信号进行降噪性能的对比分析。实验环境为i7-3.4G，内存4G，利用Matlab2010b进行实验。其中仿真信号我们采用的是含有调幅和调频的多谐波信号，并同基于EEMD的降噪算法进行性能对比。其中仿真信号的表达如下：

s=0.3sin(2π×40t)+

1.2(1+0.5sin(2π×15t))cos(2π×300t)+

sin(2π×120t+cos(2π×20t))

(19)

其中混入的噪声是方差为0.5的高斯白噪声。

为了定量地描述降噪性能本文采用信噪比作为衡量指标，其定义如下：

SNR=10 lgPs/Pn

(20)

实验首先利用EEMD算法对含噪声的信号进行分解，其中加入辅助噪声的标准偏差为0.05，集成数目为5。分解得到的前4个IMF结果如图5所示。不难发现从第二分量开始出现了严重的混频现象，并含有虚假频率成分。这一现象说明EEMD方法分解得到的IMF分量频带较宽，不具备窄带特征，进而验证EEMD算法也没能完全消除模式混叠效应，故对IMF分量进行Hilbert变换不具备严格的物理意义，导致时频图出现虚假频率。接下来对含噪声的信号进行同步挤压小波变换，参数设置同上，分解层数为7，经SST分解后得到的前4个分解结果如图6所示。从图6中的频率信息可以看出，原始信号的三个有效谐波成分已经在前三个IMF分量中呈现出来。该实验结果表明同步挤压小波变换在分解精度、准确度上都优于EEMD算法，因此在此基础上进行信号的降噪处理其性能无疑比基于EEMD的降噪方法性能更优。

图5 经EEMD分解后得到的IMF分量以及频率信息Fig.5 IMF components and frequency information after EEMD

图6 经SST分解后得到的IMF分量以及频率信息Fig.6 IMF components and frequency information after SST

接下来，对每一个经SST分解得到的IMF分量进行Hilbert变换并计算瞬时频率的四分位距值以及每一个IMF分量的自相关系数峰度值。结果如图7所示。

从图7中的曲线可以看出，经SST分解得到的各个IMF分量瞬时频率四分位距值和自相关系数峰值随着分解层数的增加基本呈现增大趋势，且两者基本保持一致，即都同时满足或不满足阈值条件，其中前5个IMF分量同时满足阈值条件。本实验中为了便于与基于EEMD降噪方法进行性能比较，我们选择前4个满足阈值条件的IMF分量进行信号合成。计算两种不同方法的SNR指标分别为17.078 4和5.043 7，由此可见本文提出的基于SST瞬时频率稳定度和自相关函数峰度阈值降噪方法的性能要优于基于EEMD的降噪方法。

图7 IMF分量的瞬时频率四分位距和自相关系数峰度值Fig.7 IMFs’ frequency IRQ and autocorrelation kurtosis

为了进一步观察选取不同数量IMF分量进行信号合成对降噪效果的影响，我们采用逐步累加IMF分量的形式来观察SNR指标变化。实验结果如图8所示，不难看出该实验结果与图7中IMF分量的判定结果相吻合。

图8 不同IMF分量累加后的SNR值Fig.8 SNR of different IMF components accumulation

这也说明本文提出的基于瞬时频率稳定度和相关函数峰度阈值进行合成IMF分量选择的可行性和准确性。图9显示了经本文算法降噪后的信号时频分析图，可以看出原始信号中的三个谐波成分已清晰地包含在降噪后的信号中，这也进一步验证了本文算法的优越性。

图9 经本文方法降噪后信号时频分析图Fig.9 time-frequency after the proposed de-noising method

3.2 工程实践降噪性能对比分析

为考察本文算法在工程应用中的有效性，利用航空物流传送设备故障测试平台收集振动信号。其中主要组成部分为：驱动电机，变速齿轮，传送轴承，皮带，模拟负载以及压电加速度振动传感器和采集终端等。其中传感器安装在传动轴承的外环部位。实验中轴承型号为N205EM (外径52 mm，内径25 mm，滚动体直径7.5 mm，个数为12个)。采样点为1 024个，采样频率2 kHz。采样信号的时序图如图10所示。由于实际应用中采集系统受外界因素，如电磁场及传感器牢固程度等，以及系统测量误差的影响，使得收集到的信号难免存在误差，因此为了说明本文方法处理工程问题的有效性，实验中加入方差为0.02的高斯白噪声，采用本文降噪方法对含噪的振动信号进行降噪处理，其中分解层数为10，其他参数设置同上。

图10 原始信号的时域波形图Fig.10 Time domain wave of original signal

计算每一个分解后IMF分量的瞬时频率四分位距值和自相关函数峰度值，其中分量号按照每个IMF分量瞬时频率的四分位距值升序排列，结果如图11所示。图12是对IMF分量逐级累加后得到的SNR曲线，结合图11的结果不难发现，只有第4个IMF分量满足阈值条件，因此依据本文的降噪方法，参与信号合成的IMF向量只需第4个分量。其得到的SNR结果为4.605 4。图中随着IMF分量的依次累加，SNR值会随之下降，这也表明了其他IMF分量含有噪声成分，参与信号合成后会严重影响降噪性能。另一方面，随着SNR逐渐降低也可以侧面地说明IMF分量的瞬时频率四分位距值及自相关函数峰度值可以作为判断一个IMF分量真实性的依据。图13为经本文降噪方法处理前后的时频分析对比结果，不难发现经本文降噪后的信号时频图更精细更清晰。

图11 IMF分量的瞬时频率四分位距和自相关系数峰度值Fig.11 IMFs’ frequency IRQ and autocorrelation kurtosis

图12 不同IMF分量累加后的SNR值Fig.12 SNR of different IMF components accumulation

图13 去噪前后信号的时频图对比Fig.13 Comparison of time-frequency before and after de-noising

为了验证不同分解层数对本文算法降噪性能的影响，试验中采用10～100，间隔为10的分解策略进行性能对比分析。为了消除随机影响，我们采用10次循环结果的平均值作为性能指标，从图14的结果可以看出，初期随着分解层数的增加，SST分解越来越完备，信噪比也会随之提高；但达到一定程度后，随着分解层数进一步加大信噪比将趋于平稳甚至略有下降的趋势，这也说明随着SST分解越来越完备，本文算法的降噪性能也就趋于稳定。另外，随着分解层数的增加耗时也会呈现线性增长，因此实际应用中只需要选取适中的分解层数即可达到性能最优。