基于SSA-VMD-MCKD的强背景噪声环境下滚动轴承故障诊断

2023-02-22甄龙信董前程张云鹏

振动与冲击 2023年3期

任良，甄龙信，赵云，董前程，张云鹏

(燕山大学河北省特种运载装备重点实验室，河北秦皇岛 066004)

滚动轴承是机械设备的关键部件，应用非常广泛，实际生产中其产生的故障往往会造成非常严重的后果[1]，其工作环境通常比较恶劣，故障信号特征往往被背景噪声和机械转动信号覆盖，不易提取[2]，因此，需要发展在强背景噪声干扰下的故障诊断方法。

滚动轴承故障信号通常呈现非线性、非平稳的特点且表现为周期性冲击特征。随着信号处理技术的发展，滚动轴承故障诊断的研究取得一定成效。文献[3]利用经验小波分解(empirical wavelet transfor, EWT)提取故障特征，再通过独立分量分析(independent component analysis, ICA)算法降噪并成功进行轴承故障诊断。文献[4]利用小波分解与FastICA算法对轴承信号联合降噪，然而这种经典降噪算法存在小波基难以选择的问题，且实际操作中FastICA算法用于轴承信号降噪的效果一般。经验模态分解(empirical mode decomposition, EMD)是由Huang等[5]提出的一种自适应分解的时频分析方法，但其存在模态混叠的缺陷，为减弱这种缺陷，文献[6]提出集成经验模态分解(ensemble empirical mode decomposition, EEMD)。文献[7]先采用最大相关峭度解卷积(maximum correlated kurtosis deconvolution，MCKD)算法降噪，再用EEMD分解，成功提取柔性薄壁轴承故障特征信号。文献[8]提出变分模态分解(variational mode decomposition，VMD)算法，相比EMD缺乏严格的理论推导，VMD具有严格的数学理论支持且可避免模态混叠的产生。文献[9]将VMD成功应用于提取轴承故障信号。然而VMD的二次惩罚因子alpha和分解层数K等参数取值对VMD分解结果影响较大，故文献[10]采用相对熵最小化作为目标函数对VMD进行参数优化，优化思路为先设定一个alpha寻最优K值，再利用最优的K值寻最优alpha值，该算法存在难以找到全局最优参数的弊端。

文献[11]提出MCKD算法，该算法通过构建和寻找最优滤波器以突出信号中的连续周期脉冲信号，但滤波器长度L、移位数M等参数对MCKD算法的效果影响较大，文献[12]采用粒子群算法对MCKD和VMD进行参数优化，并成功提取了滚动轴承早期微弱故障特征。文献[13]提出一种处理非线性非平稳信号的新方法——奇异谱分析(singular spectrum analysis，SSA)，通过对信号轨迹矩阵的奇异值分解，将信号重构成不同趋势的分量，奇异值分解可有效将信号和噪声分量分开，故该算法可用于信号降噪[14-15]。

为充分利用SSA、VMD、MCKD在降噪和特征提取方面的优势，本文提出将三者结合用于强背景噪声环境下的滚动轴承早期故障诊断，为提高诊断准确性，以包络熵极小值为目标函数，利用鲸鱼优化算法(whale optimization algorithm，WOA)优化VMD和MCKD参数。

1 相关理论

1.1 SSA算法

SSA是一种处理非线性非平稳时间序列信号的方法，可实现信号的去噪、动态重构和特征提取。本文利用SSA进行信号去噪，其核心在于利用奇异值分解原理将特征信号与噪声分离，主要步骤为如下。

步骤1输入一维时域信号[s1,s2,s3,…,sn]。

步骤2确定合适窗口长度l,得到如下轨迹矩阵

(1)

步骤3进行奇异值分解，将S分解为如下形式

S=αΣVT

(2)

式中：α为l×l的矩阵;Σ为l×(n-l+1)的矩阵;V为(n-l+1)×(n-l+1)的矩阵。

计算轨迹矩阵的协方差矩阵

C=SST

(3)

(4)

步骤4分解所得的l个分量代表着不同的趋势成分，本文通过求l个分量与源信号的互相关系数进行筛选重构。

1.2 VMD算法

VMD算法将信号x(t)分解成K个互不相关的稀疏子信uk(t)，具体分解步骤如下。

步骤1首先定义本征模态函数(intrinsic mode function，IMF)

uk(t)=Ak(t)cos[φk(t)],k=1,2,…

(5)

步骤2将uk(t)进行Hilbert变换得到单边谱

(6)

式中:δ(t)为狄利雷克分布函数;*为卷积。

步骤3将信号频谱移频至基带上

(7)

步骤4基于梯度二范数的平方估计各个uk(t)的带宽，保证所有IMF带宽之和最小

(8)

L({uk},{wk},λ)=alpha·

(9)

式中：alpha为二次惩罚因子;η(t)为拉格朗日乘子。

步骤5引进交替方向乘子算法对式(8)、式(9)进行迭代寻优，可得K个IMF，如式(10)所示

(10)

由式(5)可知K影响瞬时频率的预估，不恰当的K取值将导致模态混叠的产生，由式(9)和式(10)可知alpha值会影响IMF的带宽，故参数alpha和K的值对VMD分解效果影响较大[16]。

1.3 MCKD算法

MCKD算法能有效提取周期性冲击脉冲分量，算法的本质是在构建最合适的滤波器用于突出被噪声掩盖的特征信号。假设产生周期性冲击脉冲信号s=[s1,s2,…,sn],信号在传递过程中会混有各类噪声。观测信号的表达式为

x=h·s+n

(11)

式中：x为观测信号;s为原始周期脉冲信号;h为系统传递函数;n(t)为噪声。

构建滤波器f=[f1,f2,…,fL]，以相关峭度最大为目标函数，寻找最优滤波器，用于分离出原始周期脉冲信号，先忽略噪声，则表达式为

(12)

式中,L为滤波器长度。

1.4 基于WOA算法的VMD和MCKD参数优化

本文利用WOA[17]寻找VMD的参数alpha和K以及MCKD算法的参数L和M的最优值，均以包络熵极小值为目标函数，包络熵越小则代表信号中含特征信号越多。

信号x(i)的包络熵Ep的计算如式(17)所示

(17)

式中：a(i)为原信号x(i)经Hilbert解调后的包络信号;ε(i)为a(i)的归一化形式;N为信号x(i)的长度。

鲸鱼群体狩猎时，鲸鱼个体有两种捕食策略:一种是直接包围猎物;另一种是环形游动并产生汽泡形成气泡网驱赶猎物至气泡网中心，算法流程如下。

步骤1鲸鱼种群个体数量、位置、进化代数等参数初始化，第i个鲸鱼个体位置为

Xi=r·(ub-lb)+lb

(18)

式中：r∈[0 1];Xi的取值范围为[lb,ub],lb为待寻优参数的下边界，ub为上边界。

步骤2p<0.5且|A|<1时，WOA按照式(19)搜索

(19)

式中：r1,r2,p为[0 1]的随机数;t为迭代次数;tmax为最大迭代次数;X(t)为当前解位置;X*(t)为当前最优解位置;A,C为系数。

步骤3当p<0.5且|A|≥1时，WOA根据式(20)进行随机搜索迭代更新

(20)

式中，Xrand为随机选取鲸鱼个体位置。

步骤4当p≥0.5时，WOA根据式(21)进行螺旋收缩方式迭代更新

(21)

式中：D为个体与猎物之间的距离;b为螺线常数;l∈[-1,1]且为随机数。

步骤5若t>tmax则迭代停止，输出最优寻优结果，若不满足则返回步骤2。

WOA优化VMD和MCKD算法时，目标函数O为包络熵Ep极小值，即作为WOA算法的适应度函数。表达式如下

o=min{Ep}

(22)

参数优化的表达式如下

(23)

式中,lb和ub为待优化的参数的取值边界。

具体流程如图1所示。

图1 WOA优化VMD，MCKD流程图

2 故障诊断流程

本文提出采用SSA-VMD-MCKD方法提取强背景噪声干扰下的滚动轴承微弱故障信号特征并进行故障诊断，步骤如下。

步骤1先将信号经SSA分解降噪，采用时域互相关准则将分解产生的不同趋势的子信号筛选重构。时域互相关准则是衡量两个时域信号相关程度的有效指标，计算公式如下

(24)

步骤2重构信号经WOA优化的VMD分解提取故障特征，WOA优化时设定的参数大小为：种群规模为100，最大迭代次数为10，alpha的寻优范围为[0,10 000]，K的寻优范围为[2,10]。采用峭度准则，将分解产生的IMF筛选重构，通常轴承正常运转时信号的峭度值接近3，大于3时则认为含有较多的冲击成分，故本文选择值大于3的IMF进行重构，峭度值计算公式如下

(25)

式中：d为信号长度;u为信号均值;σ为标准差。

步骤3对重构信号包络解调，确定最大幅值所对的特征频率视为fi，根据式(16)计算MCKD的参数T，重构信号经WOA优化的MCKD分析，以加强故障特征从而确定故障。WOA优化时设定的参数为：种群规模为100，最大迭代次数为10，L的寻优范围为[20,500]，M的寻优范围为[1,7]。

步骤4对经MCKD处理的信号进行包络解调，诊断故障。

总结以上诊断步骤，绘制故障诊断流程图如图2所示。

图2 基于参数优化的SSA-VMD-MCKD的故障诊断步骤

3 仿真信号分析

3.1 单一轴承仿真信号故障诊断

为验证本文所提方法的有效性，构建单一轴承周期性冲击故障仿真信号，表达式如下

(26)

式中：s(t)为周期性冲击故障信号分量;fn为固有频率;A0为位移常数;g为阻尼系数;n(t)为高斯白噪声信号。

取fn=1 600，A0=0.8，g=0.12，冲击故障的时间间隔T为0.006 25 s,故特征频率fi=1/T=160 Hz,采样频率为12 kHz,采样点数为7 800。未添加高斯白噪声的故障信号的时域图，如图3所示。由图3可看出明显的周期冲击特征，采用短时傅里叶分析，画出信号的时频图，如图4所示，可知图4不存在零星的噪点，表明不含噪声。冲击故障的包络谱图，如图5所示。由图5可清晰看出故障特征频率fi及其倍频。

为模拟强背景噪声环境，添加-14 dB的高斯白噪声干扰。添加-14 dB噪声后的故障信号的时域图，如图6(a)所示。与图3相比，可知信号特征已被噪声严重污染。强背景噪声轴承信号的时频图，如图6(b)所示。可以看到图中存在大量噪点，故障信号的频谱图和包络谱图，如图6(c)、图6(d)所示。显然，信号特征频率及其倍频已被强背景噪声掩盖无法识别。

(a) 添加-14 dB噪声的冲击信号时域图

采用本文所提SSA-VMD-MCKD方法对仿真故障信号进行处理，先将信号进行SSA分解，窗口长度l值越大，算法运行的时间越长且不合适的取值会破坏故障特征，经多次验证窗口长度l=13时算法的效率较好，即将信号分解成13个不同趋势的子信号，计算每个子信号与原信号的时域互相关系数，结果如图7所示。通常相关系数大于0.5表示与原信号相似性较高，为避免故障特征信号丢失，本文设置阈值为0.5，故选取子信号2、信号4、信号5、信号7、信号8、信号11进行信号重构。再利用WOA优化算法寻找VMD的alpha和K的最优值，WOA优化算法的种群规模设为100，最大迭代次数为10，参数优化的目标函数即适应度函数为包络熵极小值，alpha的寻优范围为[0,10 000],K的寻优范围为[2，10],计算机参数为Intel(R) Core(TM) i5-10300H CPU@2.50 GHz,运行内存为16 GB,利用MATLAB2014a处理，寻得最优参数组合为[12,4],故设定VMD分解的alpha=12，K=4，对重构信号利用VMD分解以提取故障特征，分解所得IMF的峭度值如图8所示。由图8可知IMF1和IMF4峭度值大于3，故选取这两个信号进行重构。

图7 SSA分解各子信号的互相关系数统计图

图8 IMF的峭度统计图

重构信号的时域图和频谱图如图9(a)、图9(b)所示，与图6(a)、图6(c)对比可知经SSA-VMD分解后信号中的毛刺和噪声频带减少很多，表明降噪效果明显，为了更直观看出降噪效果，作出时频图，如图9(c)所示，对比图6(b)可知时频图上的中高频噪点消除很多，即大部分中高频噪声被滤除，但仍存在少量噪声。经参数优化的VMD分解后的重构信号的包络谱图，如图9(d)所示。选择图中最大峰值处频率作为特征频率fi，故fi=160，根据式(16)计算出MCKD算法的参数T=75；采用WOA优化算法优化MCKD的参数L和M，WOA算法的种群数设为100，最大迭代次数为10，参数优化的目标函数即适应度函数为包络熵值极小值，参数L的寻优范围为[20,500],M的范围为[1,7],寻得最优参数为[231,4],故设定MCKD算法的参数L=231，M=4，T=75。再将VMD分解后的重构信号经MCKD算法分析以加强故障特征，最后对信号进行包络解调，诊断故障。最终所得信号的包络谱图，如图10所示。显然，图中特征谱线十分突出，可轻易识别故障的特征频率fi及其2倍～6倍频，从而证实了所提诊断方法的有效性。

(a) VMD分解后重构信号时域图

图10 本文方法处理后信号包络谱图

3.2 验证所提方法的合理性与必要性

为了验证采用WOA算法优化MCKD的参数L和M的合理性，现改变L的参数，保持M不变观察结果，即随机令L=142，M值不变，经MCKD算法分析的信号包络谱图，如图11所示。从图11可以看出,特征频率及其倍频周围存在较多谱线干扰，且5倍频未成功提取。再保持L值不变，改变M值，即随机取M=1，L值不变，经MCKD算法分析的信号包络谱图，如图12所示。图12已无法发现故障特征频率及其倍频，表明本文采用WOA优化MCKD的参数L和M结果的合理性。

Sugita等已证明VMD与MCKD结合的合理性，为证明强背景噪声下SSA与VMD-MCKD结合的必要性，现将仿真信号直接通过WOA优化的VMD和MCKD算法进行分析，为了使对比试验更具对比性及说服力，WOA优化的种群规模、最大迭代次数以及VMD和MCKD的参数寻优范围均与之前设置相同。经参数优化VMD-MCKD算法分析的信号包络谱图，如图13所示。显然，图13突出的谱线频率为146 Hz，并非故障特征频率，表明强背景噪声下，以包络熵为目标函数，直接通过参数优化的VMD-MCKD算法分析，故障特征提取效果不佳，故间接证明了本文采用SSA结合参数优化的VMD-MCKD的必要性，强背景噪声干扰下，本文采用SSA分解的目的是降低噪声以便于特征提取。

4 实测信号分析

4.1 实测轴承信号故障诊断

采用凯斯西储大学公开的滚动轴承故障模拟平台试验数据进行分析。试验平台如图14所示，电机驱动端安装SKF6205-2RS型深沟球轴承，其参数如表1所示。采用电火花加工技术设置故障，电机转速为1 772 r/min,采样频率为12 kHz，本文选取7 800点进行分析。电机的转频fr及内圈故障特征频率fi的计算公式如下

图14 滚动轴承故障模拟试验平台

(27)

式中,n为电动机转速。

依据表1参数，由式(27)可计算故障特征频率的理论值为fi=159.93 Hz。

表1 滚动轴承结构参数

为了模拟强背景噪声环境，将试验信号中再添加-5 dB的高斯白噪声。信号的时域图、频谱图和包络谱图，如图15(a)～图15(c)所示。由图可知故障特征无法辨别且特征频率及其倍频已完全被转频及其倍频以及噪声干扰频带掩盖无法进行故障诊断。故障信号的时频图，如图15(d)所示。可观察到图中噪点密布,表明含噪较多。

(a) 实测信号时域图

采用本文所提方法，先将信号经SSA分解，窗口长度l设为13，分解所得子信号与原故障信号的时域互相关系数如图16所示。取相关系数大于0.5的子信号重构，即取子信号1、信号2、信号3、信号4、信号6、信号7进行信号重构。采用WOA优化VMD参数alpha和K，计算机参数与3.1节相同，WOA的种群规模设为100，最大迭代次数为10，参数优化的目标函数即适应度函数为包络熵极小值，alpha的寻优范围[0,10 000],K的寻优范围为[2,10],最优参数组合为[2 017,5]，故设置alpha=2 017,K=5,将重构信号通过VMD分解，分解所得的各个IMF的峭度值如图17所示，可知选取IMF1,IMF2,IMF3,IMF5进行信号重构。

图16 SSA分解各子信号的互相关系数统计图

图17 IMF的峭度统计图

重构信号的时域图、时频图，如图18(a)、图18(b)所示。图中可观察到冲击信号特征且时频图中噪点大幅减少，表明所提方法降噪效果明显。重构信号的包络谱图，如图18(c)所示。取图中最大幅值处所对的频率为特征频率fi,其值为160，根据式(16)计算出MCKD算法的参数T的值，可得T=75，利用WOA优化MCKD，WOA算法的种群规模设为100，最大迭代次数为10，L的寻优范围为[20,500],M的寻优范围[1,7],寻优耗时为1 726.306 3 s，最终所得最优参数为[399，4],将重构信号经MCKD算法，按照最优参数设置，L=399，M=4。经参数优化的MCKD算法处理所得信号的包络图，如图18(d)所示。图18(d)故障的特征频率及其倍频显而易见，故可确定滚动轴承故障位置在轴承内圈。

(a) VMD分解重构信号时域图

4.2 与EWT-ICA方法结果比较

吕跃刚等[3]对故障信号先进行EWT分解，计算分解所得IMF的峭度值和时域互相关系数，选取峭度值及相关系数值都较大的IMF作为快速独立分量分析(FastICA)算法的观测信号，其余信号则作为该算法的虚拟噪声通道信号，经FastICA算法解混后，再将特征信号包络解调，此方法在弱背景噪声环境下，故障诊断效果较好。为验证本文方法的优越性，采用吕跃刚等的方法对添加-5 dB高斯白噪声的实测故障信号进行分析，并将结果与本文所提方法对比。首先，将信号通过EWT分解采用“locmax”准则分解为5个IMF,每个信号的相关系数及峭度值如表2所示，选取IMF1,IMF2,IMF5作为观测信号，然后将剩余IMF作为虚拟噪声通道信号输入FastICA算法，最后将EWT-ICA算法分解所得源信号包络解调，包络谱如图19所示。可知图19只出现特征频率160 Hz，其倍频并未出现，表明EWT-ICA算法在强背景噪声环境下无法完全提取故障特征相关信息，间接表明本文方法的优越性。