基于参数优化VMD和多尺度排列熵的齿轮故障诊断

谢棕，张丽萍

(福州大学 机械工程及自动化学院，福建 福州 350108)

0 引言

齿轮及齿轮箱作为机械设备常用的调节转速和传递转矩的旋转机械设备，不仅能够传递较大的功率和载荷，而且具有较好的可靠性[1]。但是在高精度的切削加工中，当齿轮在变转速、变载荷等复杂工况下工作，极易受到损伤、产生磨损、断齿等情况，使得加工精度大打折扣。因此，齿轮的状态监测和故障诊断变得尤为重要。

齿轮的故障信号多为非线性、非平稳信号。经验模态分解(EMD)和集合经验模态分解(EEMD)在分析非线性、非平稳信号时具有良好的优越性，但是存在严重的端点效应和模态混叠问题[2]。基于此，DRAGOMIRETSKIY K[3]在2014年提出一种自适应、非递归模式的信号分解方法——变分模态分解(variational mode decomposition,VMD)，它很好地解决了EMD算法存在的问题，并且具有更坚实的数学理论基础。但VMD算法在分解过程中会同时受分解参数K和惩罚因子α这两个参数的影响。因此，如何快速、准确地确定这两个参数是研究的重点。刘长良等[4]提出一种观察中心频率的方法确定分解参数K值，当开始出现中心频率相近时，认为出现了过分解，以此来确定K值，但却忽略了惩罚因子α对分解结果的影响。瞿红春等[5]利用手动寻优对VMD参数进行优化，先假定α不变，不断寻优K值使得包络熵最小，再将此K值设为定值，不断寻优α值使得包络熵最小，得出分解参数K和惩罚因子α。但这样得出的K和α并不是最优值，且工作量巨大。基于以上分析，本文提出以IMF分量局部包络熵最小为优化目标，利用改进的蝙蝠算法优化VMD分解参数K和惩罚因子α，进而得到最优的分解参数，通过计算K个IMF分量的多尺度排列熵(multiscale permutation entropy,MPE)构造故障特征向量，输入到极限学习机(extreme learning machine,ELM)进行训练和识别，进而实现齿轮的故障诊断。

1 参数优化VMD

1.1 VMD原理

VMD分解是一种非递归模式的信号分解方法，通过迭代搜寻约束变分模型最优解，将输入信号f分解成k个离散的模态uk(t)，使得每个模态的估计带宽之和最小，约束条件为各模态之和等于输入信号f。

构造相应的约束变分模型如下：

其中：{uk}={u1,…,uk}，表示分解得到的k个模态函数；{wk}={w1,…,wk}，表示各模态分量的中心频率[6]。

为求解式(1)的最优解，引入拉格朗日乘法算子λ(t)以及二次惩罚因子α，将约束性变分问题变为非约束性变分问题。λ(t)让约束条件保持严格性，α能充分减少高斯噪声对信号的干扰：

(2)

利用交替方向乘子算法不断迭代更新各个模态及其中心频率[7]，寻求式(2)的靶点。具体实现算法如下(n为迭代次数)：

2)迭代k=1∶K并更新各模态信号uk和中心频率ωk，i∈[1,K]。

更新uk：

(3)

更新ωk：

(4)

3)根据下式更新λ，Γ为更新因子

(5)

4)若满足下方收敛条件则停止迭代，否则重复1)、2)，ε为>0的正数。

(6)

结束迭代，得到K个IMF分量。

1.2 改进蝙蝠算法

蝙蝠算法(bat algorith,BA)是由英国学者YANG X S[8]在2010年提出的一种群智能优化算法。BA算法把求解问题目标函数的适应度值度量成蝙蝠个体所处位置的优劣。最优解类比蝙蝠的猎物，脉冲响度和频率的变化在一定程度上表明与最优解的靠近程度。

在BA算法中，搜索频率、位置和速度的更新公式如下：

fi=fmin+(fmax-fmin)β

(7)

(8)

(9)

其中：fi、fmax、fmin分别表示第i只蝙蝠当前时刻发出的声波频率以及声波频率的最大值与最小值；β∈[0,1]的随机数；x*表示当前最优位置。

在局部搜索中，如果一个蝙蝠选择了一个最优解，那么将在此最优解附近随机产生一个新解：

Xnew=Xold+εAt

(10)

其中：Xold为从当前最优解中选择的一个解；At为t时刻所有蝙蝠响度的平均值；ε是[-1,1]的一个随机数。

蝙蝠脉冲的响度Ai和速率ri通过下式计算：

(11)

(12)

其中：α和γ是常数，并且0<α<1，γ>0。响度Ai和速率ri的持续迭代更新表明与最优解的靠近程度。

蝙蝠算法能够在算法运行的前期通过将全局优化转换到局部优化来实现算法的快速收敛，这同时也会导致该算法过早地处于停滞阶段。为改善这样的问题，本文在速度更新公式(8)中引入了线性递减惯性权重ω(k)的概念：

(13)

ω(k)=wmax(wmax-wmin)(Tmax-k)/Tmax

(14)

其中：ωmax为初始惯性权重；ωmin为迭代至最大次数时的惯性权重；k为当前迭代次数；Tmax为最大迭代次数。查阅相关文献知，ωmax=0.9，ωmin=0.4时算法性能最好[9]。迭代初期较大的初始惯性权重使算法保持了较强的全局搜索能力，而迭代后期较小的惯性权重有利于算法进行更精确的局部搜索。

1.3 VMD参数优化

VMD分解前，需要预先设定分量个数K和惩罚因子α的值。K和α的取值不同，会有不同的分解结果。K值过小，会造成模态混叠；K值过大，会产生虚假分量。α值过小，各IMF分量的带宽越大；α值过大，各IMF分量的带宽越小。因此，如何确定合适的K和α的值是研究的关键所在。

利用蝙蝠算法对VMD进行参数优化时，需要选定一个适应度函数。本文采用VMD分解后的局部极小包络熵作为蝙蝠算法的适应度函数，搜寻最优的VMD参数组合。信号pj的包络熵Ep可以表示成：

(15)

(16)

(17)

式中：pj是a(j)的概率分布序列；a(j)为信号x(j)经过Hilbert解调后得到的包络信号。

改进BA算法优化VMD参数步骤如下：

步骤1 BA算法基本参数初始化。主要有种群大小sizepop、最大迭代次数maxgen、参数α和K的搜索范围、最大脉冲响度A0、最大脉冲速率r0、脉冲频率范围、音量的衰减系数α以及搜索频率的增强系数γ。

步骤2 初始化种群，以局部极小包络熵作为适应度函数求出一组最优的α和K。

步骤3 根据式(7)、式(13)、式(9)不断迭代更新种群的搜索脉冲频率、速度和位置。

步骤4 生成均匀分布随机数rand，如果rand>ri，则对当前最优解进行随机扰动，产生一个新的解。

步骤5 如果计算出来的包络熵小于所得出的极小包络熵，对最优的α和K以及最小的包络熵进行更新，然后按式(11)、式(12)对响度和频率进行更新。

步骤6 重复步骤3-步骤5，直至达到最大迭代次数以及确定最小的包络熵值，得出最佳参数α和K。

2 多尺度排列熵

排列熵是由BANDT C等[10]提出的一种衡量一维时间序列复杂度的平均熵参数，常用来提取机械故障的特征。但齿轮及齿轮箱运转过程中的故障特征信息分布在多尺度中，仅用单一尺度的排列熵进行分析，会遗漏其余尺度上的故障特征信息。

针对排列熵的不足，COSTA M等[11]提出了多尺度排列熵的概念。多尺度排列熵是在排列熵的基础上将时间序列进行多尺度粗粒化，然后计算不同尺度下粗粒化序列的排列熵，具体计算过程如下：

(18)

式中：s为尺度因子；[N/s]表示取整。

(19)

式中：l表示第l个重构的分量；τ表示延迟时间；m表示嵌入维数。

3)将式(19)得到的重构时间序列按升序排列可得

S(g)=(j1,j2,…,jm)

(20)

(21)

4)当Pg=1/m!时，每种符号序列的概率都相等，此时时间序列的复杂程度最高，排列熵最大，为lnm!。为方便表示，通常会将Hp(m)进行归一化处理，即

Hp=Hp(m)/ln(m!)

(22)

Hp的取值范围为[0,1]，Hp的大小反映了时间序列的复杂和随机程度。Hp越大，说明时间序列越随机，反之，说明时间序列越规则。

3 算法流程

本文提出基于参数优化VMD和多尺度排列熵的齿轮故障诊断方法，具体步骤如下：

1)利用改进的蝙蝠算法优化VMD分解参数K和惩罚因子α，进而得到最优的分解参数，得到K个模态分量IMF1、IMF2、…、IMFK。

2)计算被选模态分量的多尺度排列熵，选取合适尺度的排列熵值作为故障特征向量。

3)将故障特征向量输入到已经训练好的极限学习机中进行训练和识别。

故障诊断流程如图1所示。

图1 齿轮故障诊断流程

4 实验分析

4.1 实验说明

为了验证本文所提方法的有效性，在实验平台下(图2)，采集试验台减速箱齿轮不同状态下的实时振动信号。齿轮减速箱型号为JZQ200，采样频率为8 200 Hz。选择齿轮正常、不平衡、松动、错位、断齿、裂痕这6种不同状态的数据进行实验验证。

图2 齿轮实验平台

4.2 故障特征提取

以减速箱齿轮出现裂痕为分析对象，本文在以下的实验中选取的蝙蝠算法参数均为：sizepop=36，参数α和K的搜索范围low=[1 500 3]，high=[2 500 8]，速度v的范围为[-100 100]，A0=0.5，r0=0.25，脉冲频率范围[0 1]，α=0.8，γ=0.5。蝙蝠算法不同参数组合的优化如表1所示。

表1 蝙蝠算法不同参数组合的优化

图3为局部极小包络熵值随寻优迭代次数变化曲线，优化结果如表2所示。

图3 寻优迭代图

表2 优化结果

作为对比，同时采用手动寻优的方法来求解最优的(K,α)的值。VMD参数设置如下：初始化alpha=2 000，tau=0，DC=0，init=1，tol=le-7。

表3 不同K值下的包络熵值

由于篇幅原因，上述表格每个K值仅陈列5个极小包络熵值，观察到最小的包络熵值仅为5.554 2。对比蝙蝠算法的迭代寻优结果，通过手动寻优求得的最小包络熵值并不理想。

综合分析，改进蝙蝠算法迭代寻优的最优参数K=6，α=247 4。选取裂痕齿轮一组数据中的4 500个样本点进行VMD分解，得到6个IMF分量，分解结果如图4所示。

图4 VMD分解结果

MPE的计算需要设置嵌入维数m，延迟时间τ和尺度因子s。其中维数m取值范围通常是3～7，m取值太小，算法的检测性能降低，若m取值太大，将无法反映时间序列的细微变化。延迟时间τ对时间序列的计算影响较小，尺度因子s最大值一般取>10即可，本文选取m=3，τ=1，s=20。

于是就能求出前面VMD分解后6个IMF分量的排列熵，如图5所示。

图5 裂痕齿轮各IMF分量的MPE值

4.3 齿轮的故障诊断

对测得的正常、不平衡、松动、错位、断齿、裂痕这6种不同状态的数据，每种数据取50组，每组4 500个点。对这50组数据进行VMD分解后并计算模态分量的MPE值，6个MPE分量，每个分量提取20尺度的因子，所以每个样本就得到维度为120的特征。将故障特征向量输入到已经训练好的ELM中，取前40组数据作为已知故障样本，后10组数据作为待识别故障样本，诊断情况如图6所示。

图6 ELM故障识别

图6中，横坐标分别代表训练样本和测试样本，纵坐标代表6种故障类别的标签。测试样本中的空心点代表预测类别。实心点代表实际类别。当测试的数据符合同类故障的标准时，空心点和实心点会重合，代表识别正确。本文测试集总的准确率达到了96.7%，对于绝大多数的故障情况，都能达到很好的识别效果。

5 结语

针对齿轮故障特征在单一尺度难以全面提取的问题，提出一种基于参数优化的VMD和多尺度排列熵的齿轮故障诊断方法。相比较手动寻优，该方法可以搜寻更优的VMD参数组合，更有效地提取出不同故障状态下的信号特征参数。分析结果表明，通过参数优化VMD和多尺度排列熵的齿轮故障诊断方法，极大地提高了齿轮故障诊断的准确性。