基于VMD和MRVM变负荷工况下的滚动轴承故障诊断*

2019-12-31周凤星黎会鹏严保康

振动、测试与诊断 2019年6期

关键词：果蝇边际轴承

徐波, 周凤星, 黎会鹏, 严保康, 刘毅, 严丹

(1. 武汉科技大学信息科学与工程学院武汉，430081)

(2. 黄冈师范学院电子信息学院黄冈，438000)

(3. 华中科技大学机械科学与工程学院武汉，430074)

引言

滚动轴承故障的智能化诊断本质上是一个模式识别过程，主要包括故障特征提取和故障识别。有效的特征提取是提高故障诊断模型性能的关键，同时对于变负荷下不同的故障特征，故障分类器的鲁棒性也至关重要。

滚动轴承故障振动信号多为非平稳信号，经验模态分解 (empirical mode decomposition, 简称EMD)[1]作为有力的非平稳信号处理工具一经提出，受到机械故障诊断相关学者的广泛关注[2-8]。张玲玲等[9]提出了基于EEMD和模糊c均值聚类算法，并对曲轴轴承的轻微磨损、严重磨损进行了有效诊断。受EMD方法的启示，Smith[10]提出另一种自适应信号分解方法——局部均值分解 (localean decompositon，简称 LMD)，也引起了广泛关注。武哲等[11]提出了基于LMD自适应多尺度形态学和Teager能量算子的轴承故障诊断方法，能够很好地提取出故障频率特征。张淑清等[12]提出了基于LMD近似熵和FCM聚类的机械轴承故障诊断方法。虽然EMD和LMD信号处理方法在非线性、非平稳信号处理中取得了很好地效果，但两者本质上均属于递归模式分解，很难从根本上解决一些固有问题。为此，Dragomiretskiy等[13]提出了一种新的自适应信号处理方法——VMD，该方法有效缓解或避免了EMD和LMD方法存在的一系列不足。刘长良等[14]引入变分模态分解，并结合模糊c均值聚类方法对滚动轴承在变负荷工况下的故障进行诊断，但外圈故障诊断精度不高。唐贵基等[15]引入粒子群优化算法对变分模态分解两个重要参数进行最优选择，取得了很好的诊断效果。但粒子群优化算法容易早熟和陷入局部最优的缺陷，限制了全局寻优能力。针对上述问题，笔者提出改进的混沌果蝇优化算法[16]对VMD两个重要参数的最优组合值进行搜索，用优化后的VMD算法对变负荷下的轴承故障信号进行处理，获得相应的IMF分量，再采用互信息法剔除虚假分量，计算有效IMF分量的二维边际谱熵值并作为MRVM的输入特征向量，用于变负荷下轴承的多种故障类型及故障程度的识别。

RVM是依据贝叶斯理论框架提出的新型分类器[17]。相比SVM方法，其解具有更强的稀疏性[18]。然而，RVM在本质上仍属于二分类器，对于多类轴承故障及故障程度的智能识别需要对多类型的数据进行精确地分类。目前，实现多分类的方法主要有一对一、一对多、决策树等方法，但它们都没有很好解决不可分区域问题，影响了分类精度。笔者提出“嵌套一对一”多分类算法构造一个精度高的多分类相关向量机，将任意两种载荷状态下的混合实验数据的二维边际谱熵值作为多分类相关向量机的训练样本，用其他载荷状态下的实验数据的二维边际谱熵值作为测试样本，进行故障诊断模型的验证，证明了所提方法在变负荷工况下的轴承故障诊断的有效性。

1 变分模态分解原理

VMD算法中，本征模态函数IMF被重新定义为一个调幅-调频信号，通过构造并求解约束变分问题将原始信号分解为指定个数的IMF分量，则对应的约束变分模型表达式为

(1)

其中：{uk}={u1,u2,…,uk}为分解得到的K个IMF分量；{ωk}={ω1,ω2,…,ωk}为各个分量的频率中心；f为输入信号；e-jωkt为预估中心频率。

为求取上述约束变分问题，引入增广Lagrange函数

L({uk},{ωk}):=

(2)

其中：α为惩罚参数；λ为Lagrange乘子。

利用交替方向乘子算法求取上述增广Lagrange函数的鞍点，为式(2)约束变分模型的最优解，而将原始信号分解为K个窄带IMF分量。

2 改进混沌果蝇优化算法

实际工程中实测的轴承故障的振动信号组成十分复杂，选定合适的参数K和α是VMD算法有效提取轴承故障特征的关键所在，同时参数K和α对VMD的性能是交互影响的，因此需要一个能够对K和α进行并行优化的群体智能算法。

2.1 果蝇优化算法

针对参数优化问题，粒子群优化算法[19]、遗传算法[20]、人工免疫算法[21]及蚁群算法[22]等群体优化算法得到广泛应用，但是它们的计算过程复杂且难以理解。Pan[23]提出了一种仿效果蝇觅食行为的果蝇优化算法(fruit fly optimization algorithm，简称 FOA)，该算法的突出特点是搜索过程简单,易于理解和实现。其主要步骤如下。

1) 果蝇优化算法的参数初始化：果蝇数量N、最大迭代次数T及果蝇的位置(Xaxis,Yaxis)。其中果蝇位置(Xaxis,Yaxis)可以通过求解如下方程来初始化

(3)

其中：rand()用于产生一个[0,1]区间内的随机值。

2) 更新每只果蝇的位置(Xi,Yi)

(4)

3) 计算对应的气味浓度判断值

(5)

4) 根据适应度函数计算气味浓度值

S(i)=Fitness(Si)

(6)

5) 最大气味浓度的选择

Index=argmax[S(i)]

(7)

6) 更新最大气味浓度

(8)

7) 检查终止条件：将当前最大气味浓度值与先前最大的气味浓度值进行比较，如果当前的最大浓度值不再优于前一个浓度，则满足终止条件，停止迭代过程。

2.2 改进Logistic混沌映射算法

FOA虽然有突出的优点，但算法的收敛精度对初始值非常敏感。为此，利用混沌思想初始化种群策略来改善果蝇优化算法，增强算法的收敛精度和收敛速度，提高全局搜索能力[24]。经典的Logistic混沌映射虽然形式简单，但有非常复杂的动力学特性[25]，也存在无穷不动点、稳定窗和空白窗、满映射状态的参数区域过小等问题。针对Logistic映射的上述问题，笔者提出基于控制参数μ与初始值z(0)的映射自变量区间分段方法，将满射范围扩大到整个参数区间。

经典Logistic混沌映射的数学模型如下

z(n+1)=μz(n)(1-z(n))

(9)

其中：μ为控制参数；z∈[0,1]。

z(n+1)=

(10)

笔者提出的改进映射的分布区间比经典的Logistic映射分布区间大，其遍历性更好。

2.3 改进混沌果蝇优化算法

使用改进混沌果蝇优化算法对VMD组合参数[α,K]进行最优组合值选择，但是优化前需确定一个适应度函数。笔者在信息熵的基础上定义信号的VMD边际谱熵Hp，其定义如下

(11)

其中：Pi为第i个频率对应幅值出现的概率；h(i)为第i个IMF分量的边际谱。

为了便于分析，将边际谱熵值进行归一化处理，HE=HP/lnL，其中L为h(i)序列的长度。根据边际谱熵的定义可知，如果IMF分量包含较多的噪声成分，其对应的边际谱熵值较大；反之，如果IMF分量主要包含轴承故障的周期性冲击成分，则VMD边际谱熵很小。将边际谱熵的最小值称为局部极小值，对应的IMF分量作为最佳分量，则VMD边际谱熵局部极小值定义为

(12)

笔者将VMD边际谱熵的局部极小值作为寻优过程中的适应度值，以局部极小值最小化作为最终的寻优目标。其优化步骤如下：

1) 初始化果蝇群体个数N、迭代次数T及上限值Tmax；

2) 初始化混沌的相关参数，随机生成2个群体个数为N的果蝇种群Uα和UK(VMD方法有2个关键参数α和K需要优化)、适应度值Fitbest≈0；

4) 根据FOA的算法并结合适应度值计算最优浓度值及对应的位置；

5) 根据预先设定的最大迭代次数Tmax，判断当前迭代次数是否达到要求，未达到则返回步骤2继续执行；若达到则停止迭代，返回α和K的最优参数值。

3 仿真实验与对比分析

3.1 仿真信号建模

滚动轴承在出现点蚀或裂纹故障时，会产生周期性的冲击信号。利用模拟强噪声背景下的轴承早期故障的仿真信号进行分析[26]，仿真信号的表达式为

x(t)=x1(t)+x2(t)+xn(t)

(13)

其中，x1(t)为周期性脉冲衰减信号，频率为12Hz；x2(t)为余弦组合信号，频率分别为25和15Hz；xn(t)为添加的高斯白噪声信号。

图1所示为仿真信号的时域波形。

图1 轴承早期故障的仿真信号

3.2 参数寻优分析

为了验证笔者提出的改进混沌FOA参数寻优的有效性与优越性，与标准的FOA和混沌FOA进行对比。首先设定相关参数：果蝇数量N=12，混沌搜索最大次数Kmax= 24，最大的迭代次数Tmax=200。详细过程如图2和表1所示。

寻优搜索过程中边际谱熵局部极小值随种群进化代数变化情况，以及搜索得到的对应的VMD最优参数组合[αopt，Kopt]见表1。从图2可以观察到，ICFOA的搜索速度最快，信号重构误差最小，收敛精度也最高。

图2 3种优化算法的搜索过程图

表1 3种优化算法的对比分析

Tab.1 Comparative analysis of three optimization algorithms

方法局部极小值迭代次数[α,K]重构误差FOA0.093130[1 200,9]1.15×10-4CFOA0.07689[920,7]6.17×10-5ICFOA0.05659[860,6]1.02×10-5

3.3 EEMD，LMD和优化VMD方法的对比分析

为了进一步验证基于ICFOA算法优化的VMD方法的有效性和优越性，分别使用EEMD，LMD和VMD方法对图1中的仿真信号进行处理，结果如图3所示。

图3 基于EEMD,LMD和VMD方法的仿真信号的分量及相应的边际谱

由图3可以观察到，3种方法都能够有效地提取仿真信号中包含25和15Hz的低频余弦分量的特征频率，以及12Hz的微弱冲击信号的特征频率及其倍频频率。对比3种方法获得特征频率的幅值可知，笔者提出的基于改进混沌果蝇优化算法的VMD方法提取出的25和15Hz的低频余弦分量的特征频率幅值，以及12Hz的微弱冲击信号的特征频率及其相应的倍频频率的幅值明显大于EEMD和LMD的特征频率幅值，表明所提方法能够为故障诊断提供更加有效的故障特征信息。

4 多分类相关向量机

4.1 相关向量机基本原理

(14)

其中：K(x,xn)为核函数；x=[x1,x2,···,xN]。

预测样本便签值为tn=yn+εn，附加噪声ε～N(0,σ2)。权值向量w=[w0,w1,···,wN]Τ，并添加附加条件wj(0≤j≤N)服从wj～N(0,σ2)分布，可得分布函数

(15)

其中：α为权值向量w先验分布的超参数向量。

(16)

若每次观测多为独立事件，得到观测结果为t的概率为

(17)

4.2 “嵌套一对一”多分类算法

轴承故障诊断需要对多类型的数据进行分类，因此，需要对RVM进行多分类扩展,其中最常用的有一对多和一对一方法。但无论是一对一算法还是一对多算法，都没有很好解决不可分区域问题。因此，笔者给出一种基于一对一多分类算法的嵌套策略来处理不可分区域问题，提高分类精度。实现步骤如下：

3) 如果落入不可分区域的样本类别大于等于3，把这些样本看成多分类问题，并利用这些样本根据一对一构造超平面；

4) 重复步骤2～3，直到该算法收敛;

5) 如果最后形成的不可分区域只包含一类样本，将该区域指派给该类，如果包含两类样本，利用二分类RVM来划分该区域并指派给相应的类别。

对于一个给定的C分类问题，假设样本数为l，当基于一对一策略构造完超平面后，对获得的超平面进行编号，如图4所示。最后，综合上述方法构建一个高精度的多分类RVM模型。

图4 嵌套一对一判别函数

为了验证嵌套一对一算法的有效性，把该算法与一对一算法、模糊判别算法、有向无环图算法在5个UCI数据集上进行比较，结果如表2、表3所示。由表2可知，嵌套一对一算法在解决不可分区域问题时表现出比一对一算法、模糊判别法、有向无环图等方法更好的性能。由表3可知，对于整个数据集，嵌套一对一算法的分类精度比其他方法更高。因此，嵌套一对一算法的总体分类精度高于其他方法。

表2 不可分区域的分类精度

表3 总体分类精度

4.3 特征向量的选择

边际谱能够准确反映信号的实际频率成分的分布情况以及信号频谱的不确定性程度。文献[22]验证了二维特征向量指标与一维特征向量指标相比具有很大优势，能够更有效地区分正常和故障信号。因此，笔者计算了VMD的IMF分量的二维边际谱熵作为多分类RVM的输入特征向量。

5 诊断流程

在实际工程中，在某一未知负载下采集到完整的故障样本是非常困难的。因此，使用已有工况下的故障数据，识别其他工况下的故障及类型才具有实际应用价值。另外，不同的载荷工况会影响振动信号的幅度，导致对损伤程度产生误诊。针对上述问题，笔者提出了用于变负荷工况下的轴承故障的智能诊断方法，诊断步骤如图5所示。

6 实例分析

6.1 实验平台

鉴于采集实际工程设备在变负荷工况下运行时的轴承故障数据十分困难，笔者使用滚动轴承故障模拟实验平台采集的变负荷工况下故障数据作为分析数据。该平台的驱动端轴承型号为SKF6205，采样频率为12kHz，采样数据序列长度为4 096，采用电火花加工不同程度的点蚀故障(向采集的故障中添加强随机噪声以模拟实际工程设备的强背景噪声环境)。其中，电机负荷包括：负荷0(0W，负荷为0N, 转速为1 797 r/min)、负荷1(768W，负荷为400N, 转速为1 772 r/min)、负荷2(1 536W，负荷为800N, 转速为1 750 r/min)和负荷3(2 304W，负荷为1 200N, 转速为1 730 r/ min)。

6.2 变负荷工况下的实验数据

为了验证所提出的方法在变负荷工况下故障诊断的有效性与鲁棒性，用上述4种载荷中的任意2种载荷工况下的实验数据作为训练样本，用剩下的2种载荷的实验数据作为测试样本，每个实验数据样本长度取4 096，详细的实验数据描述如表4所示。

6.3 故障特征提取

为了简化分析过程，笔者选择滚动体微弱故障的实验数据作为分析数据，如图6所示。根据滚动轴承结构参数及故障特征频率的理论计算公式，零载荷下滚动体的故障特征频率理论计算值fb0=

图5 变负荷工况下的轴承故障诊断流程图

表4 变负荷工况下的实验数据

Tab.4 Experimental data of variable load conditions

混合训练样本(0+1 536W)混合测试样本(768W+2 304W)点蚀半径/mm(损伤程度)轴承健康状态故障标识120+12080+800正常0120+12080+800.18外圈故障1120+12080+800.36外圈故障2120+12080+800.54外圈故障3120+12080+800.18内圈故障4120+12080+800.36内圈故障5120+12080+800.54内圈故障6120+12080+800.18滚动体故障7120+12080+800.36滚动体故障8120+12080+800.54滚动体故障9

图6 零载荷下的滚动体故障信号

141.31Hz，1 200N载荷下滚动体故障的特征频率理论计算值fb2= 135.3Hz。结果如图7所示。

图7 不同点蚀半径的滚动体故障在0和800N两种载荷下的边际谱

从图7可以观察到，使用EEMD和LMD方法处理获得的边际谱中故障特征频率几乎被其他无关频率成分淹没，而使用ICFOA优化的VMD方法处理可以得到故障特征频率及其二倍频。因此，验证了笔者提出的ICFOA优化的VMD方法能够更有效地提取轴承滚动体的微弱故障的特征频率。

6.4 故障诊断

为了验证所提方法的有效性，使用0和800N两种载荷条件下轴承的正常、外圈故障、内圈故障和滚动体故障信号的二维边际谱熵进行混合，作为训练样本；用400N与1 200N两种载荷为未知载荷，作为测试样本进行测试。实验结果如表5、图8所示。将400N和1 200N载荷下的实验数据的二维边际谱熵作为测试集合分别进行测试，由图8决策输出结果可知，能够对不同载荷下轴承的外圈、内圈和滚动体等多类轴承故障进行很好的识别，但对不同点蚀半径(模拟不同受损程度)的故障识别存在一定的误差。