麻雀算法参数优化VMD联合K-SVD滚动轴承故障诊断

褚 惟，王贵勇，刘 韬，王振亚

（1.昆明理工大学 机电工程学院，昆明 650500；2.内蒙古第一机械集团有限公司，内蒙古 包头 014000）

滚动轴承作为机械装备中的重要零部件被广泛应用于工程领域[1]。轴承的健康状态与机械装备的正常运转紧密联系，一旦发生故障便会影响到生产安全和经济效益。由于滚动轴承的振动信号蕴含了丰富的设备状态信息，因此对轴承进行振动分析已成为状态监测分析的基本手段[2]。

目前基于正交基的信号分解方法，如傅里叶变换[3]，小波变换[4]等，往往需要大量的基函数才能对故障信号进行完整表达，一定程度上难以简洁和自适应地进行信号分解。而稀疏分解以过完备字典为基础，能根据分析信号自身特点尽可能地优选字典基函数来表示信号，使所得信号更为稀疏简洁，避免了基函数的冗余，更好地刻画信号的内在特征。常用的稀疏字典有基于原始原子伸缩、平移等变换所创建的字典，如小波基字典，傅里叶基字典，离散余弦基(Discrete Cosine Transform，DCT)字典等[5]。这类方法的主要思想是设计具有特征相似度的字典，因此在与信号良好匹配的情况下是有效的，但处理某些复杂多变信号时，其性能会降低。另一类方法则是根据信号本身进行自适应学习，即基于信号本身的学习字典算法，如最佳方向法(Method of Option Directions，MOD)[6]，K-SVD[7]等。K-SVD算法通过奇异值分解进行字典更新，避免了MOD算法的大规模矩阵求逆过程，大大提高了运算效率[8]。

在工程环境中，由于工况条件的复杂多变，滚动轴承故障发生时信号往往具有非线性、非平稳、弱周期性、低信噪比等特性，使得经典的K-SVD 算法在自身学习过程中易受到噪声干扰产生虚假原子从而影响表达效果[9]。VMD[10]通过限制信号各分量的带宽，对不同中心频率的信号成分进行分离，可有效地去除噪声干扰。因此，以经VMD分解的本征模函数(Intrinsic Mode Function，IMF)分量进行K-SVD字典学习来稀疏表达能避免K-SVD 学习过程产生的虚假原子，更有效进行故障识别。首先，通过麻雀搜索算法(Sparrow Search Algorithm，SSA)优化VMD 的分解层数k和平衡因子α，并将所得优化值代入VMD中进行分解；其次利用平方包络谱峭度(Kurtosis of Squared Envelope Spectral，KSES)[11]遴选分解出最优IMF分量；最后对最优IMF分量进行K-SVD字典学习，提取有效故障信息。

1 理论基础

1.1 变分模态分解

VMD主要分为变分问题的构建和求解，即：

对式(1)引入平衡参数α和拉格朗日乘数λ使该变分问题无约束化。因此增广拉格朗日量可表示为：

用交替乘子法可得式(1)解对应式(2)的鞍点。再预先确定k，初始化参数并通过式(3)、式(4)对和ωk进行迭代更新。

更新和ωk后，再以式(5)对λ进行更新。

当满足迭代精度ε时停止，即：

式(6)中，参数k保证了分解模式数的适当性和准确性；参数α与信号的重构精度相关，二者的选择对VMD 算法的分解效果尤为重要。当k和α过大时，容易造成模态混叠，反之会造成有用信息的缺失。因此本文引入麻雀搜索算法对两者进行优化。

1.2 麻雀搜索算法

SSA算法是Xue等[12]于2020年提出的一种新的优化算法，可以概括为寻找-跟随-预警的抽象模型，它模拟了麻雀的觅食过程以获取待优化问题的解。

设M是麻雀觅食的优化搜索空间集，存在N只麻雀个体，第x只麻雀处在M集的位置表示为Sx=[sx1,…,sxd,sxM]，x=1,2,3,…,N，sxd表示麻雀x在M集中居d维的位置，寻找者位置更新表达式为：

其中：t为当前迭代数；K为最大迭代数；β是属于区间(0,1]的均匀随机数；P是服从N(0,1)分布的随机数；J为1×d的单位矩阵；R2∈[0,1]为预警值，KZ∈[0.5,1]为安全值。而跟随者位置更新可表达为：

式中：swtd为麻雀群进行t次迭代时处在d维最差的位置；反之sct+1d为t+1次迭代时最好的位置；当x＞n/2，表示适应度较低，需扩大搜索范围；当x≤n/2时，表示适应度较高，可在sc位置周围随机觅食。预警麻雀位置迭代为：

式中：δ为服从N(0,1)的随机数；V为[-1,1]的随机数，表示麻雀移动的方向，同时控制步长；e是避免分母为零而设的极小值；hx为处于位置x时麻雀的适应度值，hw为当前麻雀的最差适应度值，hg则为最优值。通常，预警麻雀个数占总种群的15%，为兼顾优化准确性和计算效率，本文设置种群数和最大迭代数[N,M]=[30,20][13]。

在对VMD 的k、α进行优化时，需考虑SSA 算法中一关键点，即适应度函数值的构建。本文选取包络熵为麻雀优化算法的适应度函数值，包络熵[14]可以很好评价信号的稀疏性，反映所研究信号分解情况的概率分布特性。

1.3 平方包络谱峭度

在故障发生时，信号中的瞬时能量变化主要受故障冲击和噪声脉冲的影响。对于任意的模态分解分量，假设SE()为平方包络信号，其方差的变化可很好表现故障信号的瞬时能量波动，表达式为：

式(10)中E(·)为数学期望，那么平方包络信号的峭度可表示为：

由式(11)可知，当d()增大时，峭度也会随之增大，而故障冲击和噪声脉冲的峭度值都偏大，当故障信号瞬态冲击循环频率过高时，其有效值会明显增大，但信号的瞬时能量变化范围反而会减小，d()降低，信号的峭度会减小，导致采用传统峭度指标容易误判最优模态。由于低信噪比条件下故障特征所具备的包络谱结构易被干扰或淹没，使其值与其他分量值相差不大，不利于筛选。因此，为了凸显故障循环冲击成分的有效性，选择了平方包络谱峭度作为筛选指标。

综上分析，当分量信号中的瞬态冲击循环频率比较高时，信号的有效值会增加，d([n])较小，分量信号具有较小的峭度值，但平方包络谱峭度却能呈现较大值，此时信号的瞬时能量变化最大，克服了传统峭度对单周期瞬态冲击敏感但对多周期冲击响应存在不足的缺陷，有利于检验故障信号中的循环冲击特性[15]。

1.4 K-SVD算法

K-SVD算法分为稀疏编码和字典更新两部分。给定D∈Rm×K，当K大于m时，称过完备字典。首先初始化并固定字典D，其中字典D常取DCT 字典。对信号矩阵Y=[y1,y2,y3,…,yn]∈Rm×n，稀疏编码过程可表示为：

式中：δ为逼近误差阈值(‖·‖q为lq范数)，X=[x1,x2,…,xn]∈RK×n为待计算的稀疏系数矩阵。信号矩阵Y通常采用Hankel矩阵，其结构表示为：

式中：H为转换算子，m=N-n+1。

当得到系数矩阵后，特征矩阵可重构为：

式中：dk、xk分别表示D的第k列和X的第k行。可以看作是分解后的第k个分量。与原Hankel矩阵具有相同的维数。可以被视为退化的Hankel矩阵，因此H记为H-1的逆运算为：

式(18)中：P等于在给定条件下i和j的总组合数。当字典更新时，K-SVD利用残差矩阵中的主成分按顺序更新原子。残差矩阵表示为：

将SVD应用于残差矩阵，分离不同奇异值对应的分量。式中：Ekj=表示为第j个奇异分量；Δ=diag(σ1,σ2,…，σn)是奇异值的降序对角矩阵，σj表示第j个奇异值。uj、vj分别是U、V的第j列。在KSVD 中，dk被U的第一列更新，xk被σ1vT1取代，通过往复迭代，字典D中的所有原子都可以依次更新。

2 滚动轴承故障诊断流程

针对经典K-SVD 算法在学习过程中易引入虚假原子导致信号稀疏不彻底以及VMD 中参数难以确定的问题，提出了基于麻雀算法优化VMD参数与K-SVD的联合诊断方法，其流程如图1所示。

图1 基于SSA-VMD联合K-SVD的故障诊断流程图

具体步骤如下：

步骤1：初始化麻雀算法种群数、最大迭代参数[N,M]、VMD 参数[k,α]优化范围。如满足迭代条件，转入下一步；否则更新各麻雀位置，继续寻优；

步骤2：将优化参数组合[k,α]代入VMD 进行分解，得到k个本征模态分量IMFs；

步骤3：计算各本征模态的平方包络谱峭度，选择平方包络峭度值最大的IMF分量；

步骤4：选择步骤3优选的IMF分量相空间构建Hankel矩阵。初始化DCT字典，设置最大迭代次数L和编码阈值δ，对Hankel 矩阵信号进行字典学习，迭代过完备字典和稀疏编码系数，重构信号并包络解调。

3 仿真信号分析

根据轴承内圈故障特点构建仿真信号[16]：

设置内圈故障仿真信号xi(t)幅值A0=1，采样频率fs=12 000 Hz，共振频率fn=3000 Hz，仿真模型的衰减系数B=1000，故障特征频率fi=120 Hz，转频fr=20 Hz，τi为服从μ=0、δ2=0.5%×fr的正态分布随机滑动系数，同时加入信噪比为-15 dB 的高斯白噪声n(t)，分析的信号长度为1×4 096。

内圈故障仿真信号时域波形如图2所示，可以看到冲击成分被高斯白噪声完全淹没，难以获取到冲击规律和有效信息。其对应的内圈故障仿真信号包络谱如图3所示，从中虽然能看到特征频率，但其周围存在严重的噪声干扰谱线，且转频和其他特征倍频也难以发现，无法有效进行故障特征识别。

图2 仿真信号原始时域波形

图3 仿真信号原始包络谱

基于本文方法，首先初始化麻雀算法参数，其中k取[2,10]，平衡因子α取[200,3000]进行麻雀寻优[17]。SSA迭代的适应度函数包络熵变化趋势如图4所示，18次迭代后出现了最小包络熵值3.622 7，此时通过优选得到[k,α]=[6,2214]。将该参数组代入VMD 对仿真信号进行分解，如图5所示，各IMF分量时域波形差异性不明显。

图4 仿真信号适应度函数迭代曲线

图5 仿真信号VMD分解结果

为筛选出包含冲击信息最多的模态分量，选择了平方包络谱峭度并与常见指标对比，归一化结果如表1所示，其中，样本熵指向模态1，欧式距离和峭度指向模态6。但模态1和模态6为虚假模态，无法给出故障特征频率。基于平方包络谱峭度最大指标确定了最优模态4，其频谱中心频率与所仿真信号的共振频率3 000 Hz 相同，可以看出平方包络谱峭度相较于传统峭度指标具有更好的性能。

表1 仿真信号IMF分量遴选指标值

选取模态4分量构建Hankel矩阵，初始化字典，设置字典维数m和原子个数n分别为140 和200，最大迭代次数L=10，编码阈值δ=1.9[18]。利用KSVD 算法进行字典学习，最优稀疏波形及其包络谱如图6和图7所示，从包络谱中可以明显观察到仿真信号的转频20 Hz、内圈故障频率120 Hz及其2倍特征频率240 Hz。

图6 仿真最优模态信号K-SVD稀疏时域波形

图7 本文所提方法仿真信号包络谱

4 实验信号分析

实验数据来源于美国辛辛那提大学IMS中心[19]。实验轴承为ZA2115 双列滚子轴承，每列滚子数为16，滚子直径为8.4 mm，节径为71.5 mm，接触角为15.17°，实验时转速为2 000 r/min，采样频率为20 kHz，采样间隔为10 min，共984组，每组信号长度为1×20 480。实验结束后发现转轴上的轴承1外圈出现损坏，依据故障轴承参数，计算得该轴承外圈故障频率fo=236 Hz，本文选取实验数据中的第400组振动信号（分析的信号长度为1×5 120）进行验证分析。原始信号的时域波形如图8所示，可以看到时域波形故障冲击幅值小、周期性不明显，且伴随大量背景噪声。其包络谱如图9所示，特征频率230.6 Hz几乎被淹没，难以进行故障识别。

图8 实验信号原始时域波形

图9 实验信号原始包络谱

采用SSA 算法对原始信号进行迭代搜索，迭代优化过程如图10所示，可看到迭代至第9 次及其之后，最小包络熵值稳定在3.628 5，即得到最终优化参数[k,α]=[5,316]，以此参数进行VMD信号分解，最终得到5个分解模态，其时域波形如图11所示。

图10 实验信号适应度函数迭代曲线

图11 实验信号VMD分解结果

为定量选择最优模态，计算模态的筛选指标平方包络谱峭度，进行归一化处理并与常见指标对比，得到的模态筛选指标如表2所示。可以看到样本熵指向模态1，欧式距离指向模态3，经过处理发现模态1、3为虚假模态。而峭度和本文所提指标均指向模态5，以此易筛选得出模态5为最优模态分量。

表2 实验信号IMF分量遴选指标值

为了更好说明平方包络谱峭度较峭度的优越性，给出峭度与本文所提指标对比图如图12所示，虽然两指标均指向模态5，但本文指标所判定的最优分量在整体均值线的上方，而峭度指标波动较大，除最佳模态分量5外同样存在其他分量在整体均线上方的情况，以此说明了本文所提指标具有更好的稳定性且相较于其他常见指标具有更好的筛选能力。

图12 峭度与本文所提指标对比图

基于上述VMD预处理结果，以模态5分量相空间构建Hankel 矩阵。初始化字典，取字典维数、原子个数分别为m=100,n=130，并设置最大迭代次数、编码阈值分别为L=10，δ=0.1。利用K-SVD学习字典对所构造的Hankel 矩阵信号进行稀疏编码和字典学习，获得已恢复至时间序列的稀疏重构信号时域波形如图13所示，可以看出经字典学习的稀疏信号冲击特征和周期性明显，对其进行包络分析如图14所示，易观察到明显的故障频率峰值230.6 Hz，证明所提方法提取轴承故障特征的有效性和泛化性。

图13 实验最优模态信号基于K-SVD稀疏结果

图14 经所提方法处理的实验信号包络谱

为了说明结合VMD算法的必要性，现与设置相同参数的经典K-SVD 算法进行对比。原始信号经经典K-SVD 处理后的包络谱如图15所示，由图15可明显看出，故障特征频率的幅值小、包络谱线峰值不突出，且在低频段分布有大量的混淆频率，难以识别故障。由此看出采用经典K-SVD 方法提取微弱故障特征存在困难，本文所提方法可有效避免此局限性。

图15 经经典K-SVD处理后实验信号包络谱

5 结语

本文针对VMD中模态分解层数k和平衡因子α难以选择的问题，提出了以SSA 算法进行迭代寻优方法。同时，针对经VMD分解后的最优模态难以选择问题，引入平方包络谱峭度遴选最优模态分量，通过对最优模态的K-SVD 字典学习和包络检波捕捉低信噪比条件下的轴承故障特征信息。仿真和试验结果表明平方包络谱峭度指标具有良好的适用性，所提采用SSA优化VMD联合K-SVD的诊断方法能够很好地在低信噪比环境中提取滚动轴承故障，具有良好的泛化性和一定的工程实际意义。