基于WOA-VMD和快速谱峭度的轴承故障诊断

杨栋,黄民,马超

(北京信息科技大学 机电工程学院,北京 100192)

0 引言

滚动轴承是旋转机械的重要支撑部分之一,长期处于复杂工况运行环境下,不可避免地会发生损伤,进而影响到整个机械设备的可靠运行[1]。据统计,在旋转机械故障造成的事故中,有30%是轴承故障导致,而在我国现投入使用的轴承有约40%需要检验,约有33%需要更换[2]。滚动轴承在使用过程中,所处的工作环境恶劣,工况复杂多变,多种因素会造成其损伤,导致故障产生。因此,为了保证机械设备能够安全、稳定、高效地运行,延长设备寿命,实现预知维修,需要对滚动轴承进行快速而有效的故障诊断。

在滚动轴承发生故障时,振动信号为非平稳、非线性信号[3],并且包含大量的噪声干扰。因此,如何有效地提取轴承故障特征频率,是滚动轴承故障诊断的关键。变分模态分解(variational mode decomposition,VMD)以其独特的自适应、非递归的特点被广泛应用到信号分解中[4]。与经验模态分解(empirical mode decomposition,EMD)相比,VMD能够很好地抑制模态混叠和端点效应。武英杰等[5]针对EMD和局部均值分解(local mean decomposition,LMD)的模态混叠问题,提出将VMD应用到风电机组故障诊断中,发现对不平衡故障具有较好的诊断效果。刘长良等[6]根据VMD的中心频率和宽带确定分解模态数,利用模糊C均值聚类完成滚动轴承的故障诊断。Antoni等[7-8]提出快速谱峭度算法,利用谱峭度原理,将短时傅里叶变换应用到谱峭度分析中,通过任意频率中心与宽带组合的窄带定位滚动轴承故障信息的最优频带,并且证明了该法在非平稳信号的特征检测中具有较好的效果,在故障诊断领域得到成功应用。孙萧等[9]将谱峭度与互补集合经验模态分解(complementary ensemble EMD,CEEMD)结合,通过对比信号不同性能指标对故障的敏感度,选择合适的特征完成轴承故障识别。针对VMD如何选择本征模态分量(intrinsic mode functions,IMF)的问题,很多研究人员做了相关研究。徐甜甜等[10]使用经验值设置VMD的参数,并结合谱峭度完成IMF选择,有效判断出齿轮箱轴承故障。张颖等[11]采用VMD和谱峭度结合的方法分析风电机组轴承故障,成功提取轴承故障频率。何凯等[12]等提出一种基于改进VMD和快速谱峭度的方法用于故障诊断。鲸鱼优化算法(whale optimization algorithm,WOA)受启发于座头鲸的捕猎方式,相比于蚁群算法、遗传算法等其他智能算法具有寻优能力强、结构简单、调节参数少等优点[13]。刘国辉[14]将WOA结合概率神经网络,与小波变换进行综合,对轴承故障损伤程度进行诊断与检测。李慧等[15]提出一种改进的WOA并与支持向量机相结合,对轴承故障诊断具有较好效果。

综上,在选择VMD模态分解数k和惩罚因子α时,大多数学者采用观察法和经验预测法,对于不同的故障数据无法进行适应性修改。针对上述问题,本文利用WOA优化VMD参数,并结合快速谱峭度进行先期信号滤波处理,以期达到最优的分解效果。

1 预备知识

1.1 变分模态分解

VMD算法是一种自适应、完全非递归的信号处理方法。该算法利用迭代搜索变分模型最优解来确定每个IMF的中心频率和宽带[16],实现固有模态的有效分离、信号的频域划分,进而完成信号的有效分解,具有较好的鲁棒性。其核心思想是构造变分问题和求解变分问题。

VMD构造变分问题时,通过约束变分模型使得每个IMF中心频率的带宽之和最小。相应的约束变分模型如下:

(1)

VMD求解变分问题时,为了求解式(1)的约束最优问题,引入拉格朗日(Lagrange)算子将约束变分问题转化为非约束问题。Lagrange增广函数为:

(2)

式中:α为二次惩罚因子;λ为Lagrange算子。

最后求解Lagrange增广函数极值,更新每个模态分量:

(3)

更新每个模态分量的中心频率:

(4)

根据上述原理,可以看出,VMD分解效果的优劣由模态分解次数k和惩罚因子α决定。参数k和α的选择将直接影响到算法的分解程度、噪声含量、模态重复率等。因此,如何合理地确定参数k和α对VMD分解效果至关重要[17]。

1.2 鲸鱼优化算法

WOA主要包括收缩包围机制、螺旋收缩机制、随机搜索机制三部分。由于目标猎物在搜索空间的位置是未知的,算法假定当前搜索代理是目标猎物或最优解,完成最佳搜索代理定义后,其他搜索代理会尝试对最佳搜索代理更新它们的位置[18]。本文以包络熵极小值作为适应度函数,包络熵值越低,样本序列复杂度就越小,样本的自我相似性就越高[19],即分解得到的IMF的噪声少,包含的特征信息多。

信号的包络熵计算如下:

(5)

式中:a(i)为经过模态分解后各分量的包络信号;ε(i)是a(i)归一化后得到的概率分布序列;N为采样点数。

算法步骤如下:

1)鲸鱼种群个体数N,个体位置X=[k,α],最大迭代次数tmax等参数初始化。

2)计算每个鲸鱼适应度,确定当前搜索代理(猎物)位置。

3)产生一个[0,1]间的随机数p,当p<0.5且|A|<1时,群体中的鲸鱼个体向猎物发起包围,即算法的收缩包围机制,如式(6)所示:

(6)

式中:Xj*(t)为当前最佳搜索代理位置向量,j=1,2,…,d;D为在第j维空间位置的鲸鱼个体向当前最佳搜索代理靠近的长度向量;|·|表示对向量中的每个元素取绝对值;Xj(t)为当前搜索代理位置向量;A和C为控制系数;r是[0,1]内的随机数;t为当前迭代次数;tmax为 最大迭代次数;α在迭代过程中从2线性减小到0。

4)当p<0.5且|A|≥1时,鲸鱼随机选择其他鲸鱼位置来搜索猎物,即随机搜索机制,如式(7)所示:

(7)

式中:Xrandj(t)为随机个体位置向量。

当p≥0.5时,选择螺旋收缩机制,如式(8)所示:

(8)

式中:b为限定对数螺旋形状的常数;l为区间[-1,1]内的随机数。

5)判断算法是否结束(即t=tmax),若是,输出最优搜索代理并结束;否则,根据对应公式更新α、A、C的值,继续进行迭代更新。

WOA优化VMD参数具体流程如图1所示。优化过程中选取包络熵作为适应度函数,计算每个鲸鱼个体适应度,通过收敛因子大小选择迭代公式进行迭代更新,直到满足终止条件,输出最优VMD参数。

图1 WOA优化VMD流程

1.3 快速谱峭度

Antoni通过Wold-Crammer提出非平稳信号谱峭度的定义[20],通过计算频域上每条谱线的峭度值,以便发现信号中非平稳成分的存在,并指出这些非平稳成分在哪个频带。

Wold-Crammer的表达式可以分解任何零均值非平稳信号x(n):

(9)

式中:dZx(f)是标准正交谱增量;H(n,f)是x(n)在频率f上的复包络。

谱峭度被定义为四阶归一化累积量:

(10)

从式中可以看出平稳过程的峭度值是频率的常数函数并且平稳信号的高斯过程的峭度值为零。瞬态冲击信号占主导的频率段的谱峭度较大,而平稳高斯噪声占主导的频率段谱峭度较小,因此,可以通过峭度值的大小,找到瞬态冲击信号所在频率段。

在存在平稳加性噪声的情况下,式(10)变为

(11)

式中:ρ(f)为信号的信噪比,其大小决定了Kx+b的取值。ρ(f)较大时,Kx+b(f)近似等于零;ρ(f)较小时,Kx+b(f)近似等于Kx(f)。因此,在整个频域上通过计算信号Kx+b(f)的谱峭度最大时的频带,可以得到非平稳信号的峭度值最大时的频带。

1.4 时域特征指标

轴承故障信号包含多种信号成分,无任何规律可循。选择合适的时域特征指标作为轴承诊断的特征参数,能够准确地表征信号的状态,进而判断轴承运行状态,因此可利用时域特征指标挑选所含故障信息丰富的IMF分量。信号分解后,IMF分量是否能够充分反应原信号的特征取决于它们之间的相关程度,因此选取相关系数法来判断IMF分量的有效性。为了避免单值参考系数存在误判或不充分的问题,同时选取峭度值作为评判指标。在轴承发生故障时,随着冲击能量的增加,峭度值会明显增加,可以准确地表征IMF分量是否包含较多的冲击成分。

相关系数计算:

(12)

式中:x(n)为原始信号,u(n)为IMF;IMF与原始信号间的相关系数Rxu∈[-1,1],︳Rxu︳值越大,分量与原始信号间的相关程度越大;M为信号序列中的点数。

峭度计算:

(13)

2 滚动轴承故障诊断方法

本文将WOA-VMD和快速谱峭度相结合,在对信号进行滤波的前提下,根据不同的信号特点选取VMD分解参数,用于处理信噪比低、故障特征不突出的滚动轴承振动信号,可以准确地提取出滚动轴承故障特征并完成故障模式诊断。

WOA-VMD结合快速谱峭度轴承诊断具体步骤如下(如图2):

1)利用快速谱峭度对原始振动信号进行计算分析,选取最优中心频率及带宽。

2)利用快速谱峭度对故障信号分析的结果构造带通滤波器,完成轴承故障振动信号的滤波降噪。

3)通过WOA优化VMD的参数,寻找基于该故障信号下的最优分解模态数k与惩罚因子α。

4)通过计算各IMF与原始滤波信号的相关系数与峭度值,挑选最优IMF。

5)对最优分量进行包络谱解调分析,实现轴承故障诊断。

图2 滚动轴承故障特征提取流程

3 实验验证与结果分析

3.1 实验说明

使用美国凯斯西储大学轴承数据中心实验数据来验证本文所提方法的可行性。实验平台由电动机、扭矩传感器、测功器、电子控制器组成。风扇端轴承型号为SKF6203-2RS,使用加速度传感器采集轴承振动信号,轴承损伤采用电火花单点损伤。实验选取转速为1 750 r/min、采样频率为12 kHz、损伤程度为0.533 4 mm的风扇端轴承故障数据。轴承参数如表1所示,实验装置如图3所示。

表1 风扇端轴承参数

图3 轴承实验台

将表1轴承参数代入到轴承故障特征频率计算公式,即:

(14)

(15)

(16)

式中:fi为轴承内圈故障频率;fo为轴承外圈故障频率;fr为轴承转频;Fv为轴承转速。求得轴承转频为29.2 Hz,轴承外圈故障频率为89.14 Hz,轴承内圈故障频率为144.46 Hz。

3.2 滚动轴承内圈故障分析

实验选取轴承内圈故障数据为例来验证所提方法的可靠性与有效性。图4所示为轴承内圈故障数据的时频域图。由图中可以看出,当轴承内圈在发生故障时,噪声对故障冲击信号有很大抑制作用,故障特征不明显,难以提取。

图4 轴承内圈故障信号时频域图

为了保证振动信号带通滤波效果,对轴承原始振动信号进行快速谱峭度分析。图5所示为轴承内圈故障信号快速谱峭度图,其中,原始信号最优中心频率fc为2 343.75 Hz,最优带宽Bw为187.5 Hz。

图5 轴承内圈故障信号快速谱峭度

根据快速谱峭度分析结果,以2 156.25～2 537.25 Hz为通带构造带通滤波器对原信号进行滤波。经滤波后的信号时域图与频谱图如图6所示。可以看到经过滤波后,故障信号的冲击成分明显,并有一定的规律性。WOA迭代次数设置为30,种群大小为10,利用WOA优化VMD参数,并分解滤波后的故障信号,得到一组IMF分量(IMF1～IMF3),如图7所示。

图6 轴承内圈故障信号滤波后时频域图

图7 IMF分量时域图

在完成WOA-VMD分解后,计算每个IMF分量与滤波信号的相关系数和峭度值,结果如表2所示。可以看出,IMF2分量的相关系数与峭度值最大,即包含故障冲击成分最丰富。对IMF2进行包络谱分析,如图8所示,轴承的转频为29.296 9 Hz,轴承内圈故障频率为145.02 Hz;与理论计算转频29.2 Hz、故障频率144.46 Hz非常吻合,其他3条主频分别为转频的二倍频58.593 8 Hz、三倍频86.425 8 Hz、四倍频115.723 Hz,以转频调制边带规律明显,由此可以判断出轴承处于内圈故障状态。

表2 风扇端轴承参数

图8 IMF2包络谱

3.3 对比实验分析

为进一步证明本文方法的有效性,直接将原始轴承故障信号进行VMD分解,模态分解数与惩罚因子根据经验值设置为6和2 000。分解得到一组IMF分量(IMF1～IMF6)。计算所有分量相关系数与峭度值,根据计算结果,选取最优分量进行包络谱分析,如图9所示。从图中可以看出,与内圈故障频率144.46 Hz相吻合的145.02 Hz处的谱线不明显,该处谱线仍然有噪声干扰,不突出。以转频(29.2 Hz)为间隔的调制边带规律也不明显。因此针对轴承内圈故障信号处理,本文方法具有明显优势。

图9 IMF分量频域图

4 结束语

滚动轴承在发生故障时,其振动信号不仅包含故障冲击成分,同时也包含大量背景噪声,具有非线性、非平稳等特点,常用信号处理方法无法有效地提取故障成分。本文利用WOA对VMD参数进行优化,结合快速谱峭度分析,分解轴承故障振动信号,判断轴承故障模式,并利用公开数据集验证本文方法的有效性,具体结论如下:

1)快速谱峭度分析结合带通滤波对初始故障信号进行预处理,能够有效剔除噪声,突出故障特征信号,进而使WOA-VMD分解得到的IMF效果更好。

2)通过相关系数法结合峭度值分析能够选取包含故障冲击成分最多的IMF,进一步提高故障诊断的可靠性。

3)经过公开试验台故障轴承数据实验验证与对比分析,本文方法得出的故障转频为29.296 9 Hz,与理论转频完全一致,故障频率为145.02 Hz,与理论故障频率144.46 Hz,相差0.56 Hz,与理论数据非常吻合,可以准确判断故障类型,并且倍频谱线清晰准确,转频调制边带规律明显。因此,本文方法能够有效地诊断出轴承故障。