基于EMD模态相关和形态学降噪的齿轮故障诊断研究

，

(石家庄铁道大学 机械工程学院，河北 石家庄 050043)

齿轮是旋转机械设备的重要组成部分，其运行状态是否正常直接影响机械设备的整体性能。据统计，传动机械中80%的故障是由齿轮引起的，旋转机械中齿轮故障占其故障的10%左右[1]。当齿轮发生故障时，故障信号受到外部环境影响而被噪声信号淹没，信噪比低，早期难以发现，加之故障冲击引起机械系统的非线性振动，导致测得的齿轮故障信号也是非线性、非平稳的，因此应用传统的信号处理方法无法获得满意的效果。

为了提取齿轮故障信息，消噪、降噪是信号处理的关键，许多专家学者提出一些有效的方法：傅里叶滤波、小波降噪，Hilbert-Huang变换等，但都有各自的局限性[2-4]。EMD作为一种自适应的信号处理算法，它可以分解得到频率由高到低的IMF分量，通过不同分量的重构可以构成不同种类的滤波器，从而达到降噪的目的。但是IMF分量的重构一直没有一个统一的标准，且重构方法因人而异，得到的效果不尽相同[5]。

数学形态学[6]是Matheron和Serra共同提出的一种新的算法，该方法可以很好地抑制噪声，其算法简单、实用有效[7]，被广泛应用于信号处理领域。如郝如江等[8-9]提出了多尺度的形态滤波器的降噪算法，并在轴承故障特征提取方法中取得良好效果。但形态学滤波器存在如何选择最优结构元素的问题。

针对上述问题，本文提出了基于模态相关法的分选准则和基于峭度准则优化结构元素尺度的自适应形态滤波算法，二者结合更好地达到了降噪、消噪的目的。以仿真数据和齿轮故障实测信号进行了验证，结果表明该方法可以更好地提取出齿轮故障特征信息。

1 自适应形态学滤波

1.1 形态学滤波器

形态滤波的原理是以信号的时域特征为基础，通过选择一定形状的结构元素与被分析的信号进行适当的形态运算，从而达到降噪、突出信号特征的目的。它有4个基本算子：腐蚀、膨胀、形态开和形态闭[10]。由于论文的研究对象均是一维离散信号，因此仅给出一维的形态学算子定义。

设采样得到的一维多值信号f(n)和一维结构元素序列g(m)分别为定义域F=(0,1,2，…,N-1)和G=(0,1,2,…,M-1)上的离散函数，且N≥M，则f(n)关于g(n)的腐蚀和膨胀分别定义为

(fΘg)(n)=min[f(n+m)-g(m)],m∈G

(1)

(f⊕g)(n)=max[f(n-m)+g(m)],m∈G

(2)

式中，Θ为腐蚀算子；⊕为膨胀算子。

数学形态学的这4种基本运算都可以提取信号的轮廓信息，但其作用不同。腐蚀运算可以抑制正冲击，平滑负冲击；膨胀运算则相反；开运算用于滤除信号上方的峰值噪声，去除信号边缘的毛刺；闭运算则用于平滑或抑制信号下方的波谷噪声，填补信号的漏洞和裂纹。也就是说，腐蚀和开运算的效果类似，膨胀和闭运算的效果则与之相反。实际应用中，常用的形态滤波器是这几种算子的组合。

(1)平均(AVG)滤波器

(AVG)(f)=(f·g+f∘g)/2

(3)

(2)差值(DIF)滤波器

(DIF)(f)=f·g-f∘g

(4)

开运算抑制正冲击、闭运算抑制负冲击，两种滤波结果进行平均可以起到降噪的目的。差值滤波器可以提取信号中的正、负冲击。

1.2 基于峭度准则优化形态学结构元素

结构元素是形态学的基本算子，形态学滤波效果的好坏很大程度上取决于结构元素，而结构元素的选择在于其形状和尺寸(应根据待分析信号的特点而定)的确定，结构元素与待提取的特征形状越相似，提取效果越好。结构元素可以选择为圆盘形、三角形、余弦、正方形、扁平形等。文献[11]的研究结果表明，结构元素形状的选取对滤波效果影响不大，但整体来说余弦、半圆、三角优于其它形状的结构元素。因此为了计算简单，选取幅值为1的余弦结构元素。当其形状和高度确定之后，结构元素长度的最优解就成了形态滤波中的重要工作。

对于一个离散变量x，归一化的4阶中心矩称为峭度(kurtosis)，定义为

(5)

式中，E(x-μ)4是4阶数学期望，μ是均值；σ是标准差。在一定范围内，故障信号的峭度值随着故障的严重程度而增加，峭度的这一特性适用于齿轮故障特征提取。

本文利用峭度指标(峭度越大，滤波误差越小，则提取效果越好)作为衡量标准来自适应地寻优最佳尺度值。结合实验工况和理论计算的齿轮故障频率，将结构元素尺度设定在[2,100]，即从2循环到100，找到峭度值最大时对应的尺度值，即为最优解。

2 EMD滤波理论

EMD分解方法是将一个待分析信号分解为若干个频率由高到低的内禀模态函数，分解的整个过程就是一个“筛分”的过程，表现出EMD本身的一个自适应滤波特性。根据它的这个特点，我们可以将所需要的IMF分量重构，从而得到某一个频段的信息，这便相当于一个滤波器。

在实际应用中，EMD算法的滤波效果具有一定的局限性，具体表现为：

(1)EMD将信号分解为N个频率由高到低的IMF分量，其中阶数低的对应高频成分，一般认为这些高频成分是由噪声主导的，重构信号时直接去除，但是具体去除哪几个，没有相应的准则来判断。

(2)即使人为地判断出哪些是噪声主导的模态分量，但对实际测得的信号而言，这些分量中也包含一些故障信号的高频成分，直接去除这些分量会导致故障信息的缺失，影响后续齿轮故障的判断。

针对上述问题，提出基于模态相关分选准则和形态学滤波相结合的降噪方法，该方法的具体步骤以及用于齿轮故障信号的分析如下。

3 基于模态相关分选准则和形态学滤波相结合的降噪方法

相关系数主要反映两个函数的相关程度，相关系数越大，表明两个函数之间的相关程度越高，相反，相关程度越小。经过大量实验发现，各阶模态与原始信号的相关系数具有较好的数据稳定性。本文将选用各阶模态与原始信号的相关系数作为噪声模态与信号模态的判断准则，称之为模态相关分选准则。其相关系数定义如下[12]

(6)

式中， cov()代表协方差。

文献[12]用具体的数学推导过程证明了该准则的合理性。另外，本文采用各个IMF分量的自相关函数对其进行了进一步地验证。现就自相关函数的定义和相关特性介绍如下：

自相关函数作为信号时域中的一种统计特性，主要反映信号在两个不同时刻时的相关程度，其定义为

Rx(t1,t2)=E(x(t1)x(t2))

(7)

实际应用中一般对其归一化处理，即

(8)

式中，x(t)为随机信号；τ为时间差，即τ=t1-t2。

按式(7)和式(8)分别得到高斯白噪声和一般信号的自相关函数，如图1所示。其中，一般信号是由频率均为10 Hz的正弦和余弦成分叠加而成，并将其归一化处理。

图1 自相关函数图

由图1(a)可知，高斯白噪声的自相关函数在零点处为1，其它点处迅速衰减到零；由图1(b)可知，一般信号的自相关函数在零点处取得最大值，其它点处不是迅速衰减到零，而是缓慢变化的，这与噪声信号存在显著的区别，因此可以根据各个IMF分量的自相关函数来判断模态相关法分选准则的正确性。

对于判断出的噪声起主导作用的IMF分量，传统的EMD滤波算法是将其直接去除，事实上这些分量中还包含少量有用信号的高频成分，为了尽可能多地保留有用信号的特征信息，本文利用自适应形态学滤波器对其进行降噪处理，然后将处理后的IMF分量与未处理的IMF分量进行重构得到滤波重构信号。

综上，所述方法流程如图2所示：

(1)给定被噪声污染的信号x(t)，将其进行EMD分解，得到N个IMF分量；

(2)分别计算各个分量与原始信号的相关系数，做出相关系数分布图，以及各个分量的自相关函数图；

(3)从相关系数图中找到噪声主导分量与信号主导分量的分界点k，并根据自相关函数的特点，验证所得分界点的正确性；

(4)对判断出的噪声主导的分量IMF1～IMFk，进行自适应形态学滤波处理，得到处理后的分量IMF1′～IMFk′；

图2 基于EMD模态相关和形态学降噪方法流程图

4 实验分析研究

4.1 仿真数据分析

为了验证所述方法的有效性，下面以仿真信号为例来说明，其中采样频率为1 024 Hz，采样点数为1 000，假设原始信号为

图3 仿真信号时域图

x1(t)=sin(2π×10t)+cos(2π×10t)

(9)

该信号由频率均为10 Hz的正弦信号和余弦信号叠加而成，如图3所示。其加入噪声之后的信号为

x(t)=x1(t)+x2(t)

(10)

图4 仿真信号加噪时域图

式中，x2(t)为高斯白噪声；x(t)为加噪之后的信号(信噪比为2)，如图4所示。其降噪的目的就是要得到一个最接近信号x(t)的估计信号x′(t)，使得x(t)-x′(t)最小。其具体仿真步骤及结果分析如下：

(1)EMD将信号分解为9个IMF分量，计算每个IMF分量与原始信号的相关系数，如图5所示。最后，做出各个IMF分量的自相关函数，如图6所示。

(2)由图5可知，噪声模态和信号模态之间的分界点为4，即前4个IMF分量是噪声占主导的分量，后5个分量是信号占主导的分量，并且由图6可知，前4个分量的自相关函数图满足噪声自相关函数的衰减特性，从第5个分量开始满足一般信号自相关函数的特性，进一步证明所提方法的准确性。

图5 各IMF分量与原始信号的相关系数

图6 前5个IMF分量的自相关函数

(3)由于前4个分量中也包含一些有用信号的高频成分，所以对其进行自适应形态学滤波处理，通过峭度准则优化形态学结构元素尺度，选择峭度最大时对应的尺度值，如图7所示，选取尺度为13的余弦结构元素。最后，将滤波之后的分量与后面的5个分量进行重构，得到滤波重构信号，如图8所示，由图可以看出，重构信号与原始信号之间相差并不大，几乎重合，充分证明了所提方法在降噪方面的可行性。

图7 不同尺度形态滤波的峭度值

图8 滤波重构信号与原始信号对比图

4.2 齿轮故障信号的特征提取分析

为了进一步证明所述方法的可行性，采用实测齿轮故障信号进行分析研究，实验采用江苏千鹏诊断系统-QPZZ-II，其中输入轴转速为870 r/min，采样频率为5 kHz，采样点数为65 536，故障齿轮(小齿轮)为3.2 mm的裂纹，齿轮参数和故障频率如表1所示，采集到的齿轮故障信号的时域波形和频谱图如图9所示。

表1 齿轮参数及故障频率

图9 实测信号时域波形和频谱图

由图9可以看出，信号的频率成分被淹没在噪声之中，无法找到信号的故障频率，对其直接进行自适应形态学滤波处理，得到的时域波形和经形态学差值滤波器解调之后的结果如图10所示。从图10(b)中可以找到信号的故障特征频率(14.5 Hz)和倍频及其边频信息，但幅值较小，二倍频和三倍频不是特别明显。

采用本文所提方法进行处理，首先将采集得到的信号进行EMD分解，得到19个IMF分量，并利用公式(6)计算每一个分量与原信号的相关系数，做出相关系数图，如图11所示。由图可知，噪声模态和信号模态的分界点是第9个分量，并对每一个IMF分量做了自相关分析验证了本文所提方法的准确性。由于前9个分量中也可能包含齿轮故障信号的高频成分，直接去除将导致故障频谱信息不完整，所以引入自适应形态滤波算法，将其进行滤波处理，通过峭度准则优化形态学结构元素尺度，如图12所示，选取尺度为37的余弦结构元素。最后，将滤波之后的分量与未处理的分量进行重构，得到滤波重构信号，然后对其进行Hilbert包络解调，结果如图13所示。由图可知，本文所述方法有效地提取出了齿轮的故障频率、倍频以及边频信息。边频信息数量多、幅值不高、均匀分布，倍频成分比较明显，幅值较高、数量多，与所设的齿轮裂纹故障相吻合。并与图10(b)对比发现，本文方法明显比直接进行形态学滤波的效果更佳，得到的频谱图中也更能凸出齿轮的故障特征频率，其前三阶故障频率幅值分别提高107%、243%和190%(见表2)。充分证明，本文所提基于模态相关分选准则和自适应形态学滤波的降噪新方法可以更好地对信号进行降噪处理，可以更有效地应用于齿轮故障诊断研究。

图10 形态滤波处理后的结果

图11 各阶分量与原始齿轮信号的相关系数

图12 不同尺度形态滤波的峭度值

图13 滤波重构解调后的频谱图

表2 故障特征频率幅值比较

5 结论

将EMD模态相关法和形态学滤波相结合，首先对信号进行EMD分解，计算分解得到的IMF分量与原始信号的相关系数，找到噪声主导模态与信号主导模态的分界点，并利用各个IMF分量的自相关函数验证该分界点的准确性。其次将选取的噪声模态利用自适应形态学滤波进行降噪处理。最后，将降噪后的模态与剩余模态重构，并对其进行Hilbert包络解调得到故障特征频率。仿真数据和齿轮裂纹故障实验分析均表明该方法可以更加有效地提取齿轮故障特征频率。主要结论有：

(1)针对IMF重构中噪声模态和信号模态分界点的选取问题，提出了基于相关系数的模态分选准则，并利用各个IMF分量的自相关函数进行验证，得出该准则的准确性。

(2)针对形态学滤波器中结构元素的选取问题，提出了基于峭度准则优化结构元素的自适应算法，实验结果表明该方法可以提高形态学的滤波能力。

(3)将EMD模态相关分选准则和自适应形态学滤波相结合能够明显地从齿轮振动信号中提取出故障特征频率，比直接进行形态学滤波效果更佳。