基于CEEMDAN与PSO-MOMEDA的滚动轴承故障诊断方法*

2022-11-09郭金泉吴欣然钟建华

组合机床与自动化加工技术 2022年10期

郭金泉，吴欣然，钟建华

(1.福州大学机械工程及自动化学院，福州 350108；2.福建省力值计量测试重点实验室，福州 350100)

0 引言

随着旋转机械在生产制造、交通建设、能源设备等工业场景中应用越来越广泛，对其性能和工作效率提出了更高的要求[1]。而滚动轴承作为旋转机械最关键同时也是最容易受到损伤的部件，其工作环境复杂多变，有时甚至伴随强烈噪声，故障信息中的冲击成分常淹没在强噪声中，十分微弱，难以提取。因此，如何有效提取故障特征，诊断轴承故障对提高工作效率及消除安全隐患具有重要意义[2]。

此前常使用的时频域方法如短时傅里叶变换，小波变换等存在一些固有的不足，包括窗函数不能同时满足既有高频成分又有低频成分的信号分析，小波基函数和小波分解层数的选择不具有自适应性等[3-5]。HUANG等[6]提出经验模态分解(empirical mode decomposition，EMD)实现了非线性信号的自适应分解，将信号分解成若干个频率不同的本征模态函数(intrinsic mode function，IMF)分量和一个残余信号，实现了重大突破。之后众多学者基于EMD原理对其进行不断地改进，WU等[7]在EMD分解前加入白噪声，提出集成经验模态分解(EEMD)。TORRES等[8]提出了自适应噪声完备集合经验模态分解(complete ensemble empirical mode decomposition with adaptive noise，CEEMDAN)算法，进一步抑制了EMD的模态混叠现象，解决了EEMD计算量大、分解完备性差等问题，具有明显的优越性，近几年被广泛使用。

MCDONALD等[9]提出了多点最优最小熵解卷积(multipoint optimal minimum entropy deconvolution adjusted，MOMEDA)方法，以最大D-范数为目标，并引入目标向量t来促进理想周期脉冲序列的反褶积，采用非迭代方式找寻最佳逆滤波器的信号降噪方法，克服了最小熵解卷积(minimum entropy deconvolution，MED)方法在解卷积过程只能突出局部几个故障脉冲成分，并不能反映故障的真实情况[10]。而最大相关峭度解卷积(maximum correlated kurtosis deconvolution，MCKD)也存在一些固有缺点，限制了它在许多实际应用中的使用，例如移位数的定义和与故障相关的循环频率的估计必须非常准确，相关峭度的提出是经验的，没有对其统计性质进行研究等[11]。但是MOMEDA方法中的故障周期搜索范围T和滤波器长度L是影响滤波效果好坏的重要参数，且一般依靠人为经验选取，因此采用粒子群(PSO)算法[12-13]对这两个参数进行优化选用。

基于以上理论背景，为了充分体现MOMEDA在增强故障冲击成分方面的优点，本文尝试将CEEMDAN和MOMEDA相结合，根据峭度值和相关系数筛选最优分量，求包络谱以提取故障特征频率，对滚动轴承的故障进行诊断。为实现MOMEDA的两个重要参数的自适应寻优，并考虑到两个参数的相互影响，采用粒子群优化算法，以包络谱峰值因子作为适应度函数对两个参数寻优。通过仿真信号和风机滚动轴承实验台实测信号分析，与CEEMDAN+PSO-MCKD方法和其他方法对比，证明CEEMDAN和MOMEDA相结合的必要性，能够对轴承故障特征准确提取和诊断。

1 方法原理及诊断步骤

1.1 CEEMDAN算法原理

设定CEEMDAN第k个IMF分量为IMFk，算子Ek(·)表示EMD分解得到的第k个IMF分量，则CEEMDAN的具体算法步骤如下：

步骤1：设滚动轴承时间序列信号为x(t)，Ni(t)为添加的单位方差均值为0的高斯白噪音，m为添加白噪音的次数，εk为所添加高斯白噪声的标准差，对x(t)添加m组高斯白噪音后进行EMD分解；

(1)

然后对分解得的一系列IMF进行算术平均计算得到IMF1：

(2)

步骤2：计算第一个残差；

r1(t)=x(t)-IMF1(t)

(3)

接着在r1(t)中添加高斯白噪音进行EMD分解；

(4)

然后对分解得的一系列IMF进行算术平均计算得到IMF2：

(5)

计算第二个残差：

r2(t)=r1(t)-IMF2

(6)

步骤3：根据步骤1和步骤2依次类推IMFk和rk(t)；

(7)

(8)

rk(t)=rk-1(t)-IMFk

(9)

步骤4：重复执行以上步骤，直至无法满足EMD的条件，分解终止，最后的残差为R(t)，即信号x(t)可写做：

(10)

1.2 峭度-相关系数筛选分量方法

CEEMDAN将振动信号分解为一系列IMF分量，通常只有几个IMF分量包含大量故障信息，而剩余分量则是与故障无关或由于分解迭代误差和环境噪声等影响产生的虚假分量和噪声分量。因此，采取峭度-相关系数准则，剔除一些与故障无关的IMF，提高故障诊断的准确性。

峭度[14]是一个无量纲参数，可以反映随机变量分布特性，对信号的冲击成分敏感。当轴承处于正常状态时，其振动信号近似正态分布，峭度系数约为3；而当轴承出现故障时，冲击成分增多，峭度随之增大。其计算公式为：

(11)

式中，E(x)为信号的期望值；μ和σ分别为信号的均值和标准差。

相关系数[15]反映各IMF与原信号之间的相关性，可以判别出虚假分量并剔除与原信号相关性较小IMF。其计算表达式为：

(12)

式中，c1(t),c2(t),…,cn(t)代表各IMF分量。

1.3 MOMEDA原理

设采集到的信号组成为：

x(t)=h*y(t)+e(t)

(13)

式中，x(t)为采集到的信号；y(t)为故障冲击信号；h为传递函数；e(t)为噪声。

由上式可见，由于噪声和传递路径，原故障信号y(t)成为了传感器采集到的信号x(t)。因此MOMEDA的目的就是寻找一个逆滤波器f，使y(t)能被还原出来，即：

(14)

MOMEDA在D-范数的基础上提出了多点D-范数，并求最大多点D-范数，定义为：

(15)

(16)

式中，t为目标向量，是由0和1组成的常量向量，定义了解卷积目标脉冲的位置和权重。

上述优化问题等价于求解方程：

(17)

令X0=[M1M2…Mk]，则式(17)可以简化为：

(18)

(19)

将y=X0f代入式(19)，得：

(20)

取上式的特解，记为：

(21)

式(21)即为一组最优滤波器，代入式(14)中可恢复原始故障信号。

1.4 方法步骤

本文采用CEEMDAN+PSO-MOMEDA的方法对噪声下的滚动轴承振动信号进行故障诊断，具体步骤如下：

步骤1：使用CEEMDAN对振动信号进行分解，并根据峭度-相关系数准则筛选最优分量；

步骤2：使用PSO对MOMEDA的参数[L,T]进行寻优。PSO设定参数：粒子个数为30，最大迭代次数为20，最大权重为0.9，最小权重为0.4，学习因子为1.5。滤波器长度取值范围100≤L≤1000；窗函数选取长度为3的矩形窗；故障周期搜索范围[Ts，Tf]，且Ts和Tf分别为故障周期搜索的初始值和最终值，一般T的计算式为T=fs/fi，其中fs为采样频率，fi为故障特征频率(To=fs/fo为外圈故障周期、Ti=fs/fi为内圈故障周期、Tb=fs/fb为滚动体故障周期)，取Ti≤Ts≤Tb，Tf取值略大于整个传动系统中理论轴承最大故障特征频率对应的故障周期并取整；

步骤3：PSO的适应度函数采用包络谱峰值因子(crest factor of envelope spectrum，Ec)[16]，其计算公式为：

(22)

式中，X(n)(n=1,2,…,N)为信号包络谱幅值序列。Ec越大，周期冲击特性越强，故障特征越明显，因此选取Ec最大值确定MOMEDA最佳参数。

步骤4：对最优分量使用最佳参数进行MOMEDA解卷积，再对解卷积后的信号包络谱分析，再将计算出的轴承理论故障特征频率值与频谱中的突出峰值。

方法流程图如图1所示。

图1 方法流程图

2 仿真信号分析

为验证本文所提方法的可行性，构造了一个模拟内圈故障仿真信号：

(23)

式中，s(t)为周期性冲击成分；n(t)为高斯白噪声；幅值为A0为0.5；τi为第i次冲击相对于周期T的微小波动；衰减系数C为800；共振频率fn为4000 Hz；转频fr为30 Hz，内圈故障特征频率fi=1/T=120 Hz，随机波动服从零均值正态分布，标准差为转频的0.5%，仿真信号信噪比为-14 dB，采样频率fs为12 800 Hz，采样点数N为8192。不添加噪声及添加噪声的仿真信号时域波形和频域图如图2所示。图2a可以观察出明显的冲击成分，而在图2b中冲击成分淹没在噪声中，图3c也无法观察出周期性的突出谱线，无法准确辨认故障特征。

(a) 不加噪时域波形 (b) 加噪时域波形

采用本文所提方法对仿真信号进行分析。首先使用CEEMDAN算法分解仿真信号，得到13个分量，并计算各分量峭度及相关系数如图3和图4所示，虽然分量8的峭度最大，但是相关系数十分小，可以推断分量8含添加噪声较多，因此选择相关系数最大且峭度同样也很大的分量1作为最优分量。

图3 各分量峭度值图4 各分量相关系数

其次根据第3节确定Ts的取值范围为[100,400]，再对[L,Ts]寻优，结果如图5所示。图5a为指标Ec随种群进化迭代数的变化结果，可知PSO寻优过程中种群进化到第17代得到最佳参数组合[586.5,107.1]。图5b为最优分量解卷积后信号波形图，图5c为最优分量解卷积后信号包络频谱图。可以看出解卷积后的信号中故障特征频率及其倍频对应的谱线清晰可见，故障特征被准确提取，证明了方法对仿真信号的有效性。

(a) 迭代曲线 (b) 解卷积后信号时域波形

为证明峭度-相关系数准则筛选分量的准确性，对峭度最大的分量8进行PSO-MOMEDA解卷积包络谱分析，如图6所示，可以观察到分量8包络频谱虽然也能提取出故障特征频率及其倍频，但是杂频较多且幅值较小，与上文推测相符，也说明了通过峭度-相关系数准则选择分量1为最优分量具有优越性。

为说明MOMEDA算法在此方法中的必要性，直接对最优分量进行包络解调，如图7所示，故障特征频率被提取出来了，同时还存在一些噪声产生的杂频，因此需要经过MOMEDA的进一步滤波。

图6 分量8解卷积后信号包络谱图7 CEEMDAN最优分量包络谱

3 实测信号验证

在Spectra Quest公司设计的风力发电模拟系统(wind turbine simulator，WTS)进行轴承内圈故障试验验证本文所提方法。实验装置如图8所示。轴承型号为ER-12K，电机的输入转频为0.798 Hz，采样频率fs为12 800 Hz，采样点数为12 801点。相关参数如表1所示。