基于参数优化变分模态分解的滚动轴承故障特征提取方法

郑 义， 岳建海， 焦 静， 郭鑫源

(北京交通大学 机械与电子控制工程学院车辆工程系，北京 100044)

滚动轴承是旋转机械设备中重要部件之一。设备运行时, 磨损、疲劳、腐蚀、过载等原因都可能造成滚动轴承的局部损伤故障[1]。过多的损伤诱因以及恶劣的工作条件导致轴承故障多发，从而降低生产效益。有效的轴承早期故障诊断方法对于提高机械可靠性，降低财产损失有着极大的帮助。

机械故障诊断的关键是提取有效的故障特征信息。然而早期故障产生的微弱冲击特征信息容易被复杂的背景噪声淹没，单一的时域及频域分析难以提取到故障特征，不适合分析具有非平稳性特点的振动信号[2]。Huang等[3]提出的经验模态分解(Empirical Mode Decomposition, EMD)方法作为一种有效的时频分析方法，避免了以往方法基函数的选取问题。它通过信号在时间尺度上的动态特性，自适应地将其分解为有限固有模态函数(Intrinsic Mode Functions, IMFs)的和，在应用中取得较好的效果[4-5]。但这种递归模式分解的方法会导致分解时误差的不断传递，出现模态混叠现象。此外还存在着端点效应、虚假分量、缺乏数学理论基础等问题，影响特征提取效果。

Dragomiretskiy等[6]提出变分模态分解(Variational Mode Decomposition, VMD)方法，这种信号分解方法完全不同于EMD方法，它将模态分量看作具有一定中心频率的有限带宽，分解具有较好的稀疏性。相对于EMD方法来说，VMD方法具有数学理论基础完备，模态混叠现象及端点效应较小等优点，在近年来得到广泛应用。文献[7]将变分模态分解算法成功应用到转子故障诊断中；文献[8]将变分模态分解结合奇异值分解方法应用到轴承故障诊断中，利用奇异值分解的降噪能力对变分模态分解后的重构信号进一步降噪处理，增强故障信息。但VMD方法也存在着一个较为关键且必须解决的问题，那就是它的分解效果受多个参数的影响较为严重，其中较为关键的参数有模态数k和二次惩罚项α。模态数k选取的过大过小会造成过分解、欠分解现象；二次惩罚项α主要对模态分量的带宽造成影响[9-10]。

为实现VMD方法参数的自适应选取，提出利用基于种群启发的蝗虫优化算法(grasshopper optimization algorithm, GOA)[11]对VMD方法进行优化。考虑到旋转机械故障冲击的周期性特征，将对周期冲击较为敏感的相关峭度[12]指标引入优化算法中作为适应度函数。原始故障信号经优化变分模态分解后得到k个模态分量，筛选相关峭度最大的模态分量进行包络解调分析，提取故障特征。仿真信号和实测信号的验证表明，该方法能准确提取出强噪声背景下的轴承故障特征信息。

1 基本原理介绍

1.1 VMD算法

变分模态分解方法将“模态”定义为：围绕着固定中心频率的、有限带宽的调幅调频函数。基于这种“模态”定义，变分模态分解方法将分解过程转化为变分问题的构造和求解过程。

1.1.1 变分问题构造

首先假设每个模态分量是具有一定中心频率wk的有限带宽，变分问题就可以表示为寻求k个模态函数uk(t)，使得每个模态的估计带宽之和最小，约束条件为各模态分量之和为输入信号x(t)。

构造过程中，首先通过Hilbert变换得到信号的单边频谱；随后对解析信号与预估中心频率进行混合，将每个模态分量的频谱调制到相应的预估基频带；计算平移后信号梯度的L2范数平方来估计带宽。整个变分问题的构造过程最终得到的表达式如下

(1)

1.1.2 变分问题的求解

为求解上述变分问题，引入二次惩罚因子α和拉格朗日算子λ(t)，将约束性变分问题变为非约束性变分问题。其中α可以在高斯噪声存在的情况下保证信号的重构精度，拉格朗日算子使得约束条件保持严格性，扩展的拉格朗日表达式如下

L({uk},{wk},λ)∶=

(2)

随后利用交替乘子法(ADMM)迭代求解变分问题，模态分量uk、中心频率wk、和λk的迭代求解公式为

(3)

(4)

(5)

1.2 蝗虫优化算法

VMD算法相较于EMD及其他分解算法，模态分量具有较好的稀疏性，混叠现象较小。但VMD算法的分解结果好坏受多个参数影响，其中模态数k及二次惩罚项α影响较大，是学者们关注的焦点。

本文选取用于优化VMD参数的算法是蝗虫优化算法，属于种群优化算法的一种，由Seyedali等提出。GOA优化算法通过建立用于模拟蝗虫觅食种群行为的数学模型来达到参数寻优的目的。

模拟蝗虫种群行为的数学模型为

Xi=Si+Gi+Ai

(6)

式中:Xi代表第i个蝗虫的位置;Si、Gi和Ai分别代表第i个蝗虫所受的种群影响、重力作用、及风场平流影响。各自的计算公式如下

(7)

s(r)=fe-r/l-e-r

(8)

式中:|xi-xj|代表第i个蝗虫和第j个蝗虫之间的距离，(xi-xj)/dij为单位向量;N代表蝗虫数量，及参数寻优中搜索代理的个数;s(r)用于模拟蝗虫群体中个体之间的影响力(吸引及排斥)。式(8)中，f及l的大小一般分别选取1.5和0.5。

(9)

(10)

(11)

(12)

由式(11)可以看出，蝗虫出现的下一个位置与当前位置、目标位置以及其他蝗虫位置相关，式中的第一部分代表其他蝗虫所处位置对当前位置的作用力，第二部分表示当前位置趋向目标位置的作用力。参数c模拟蝗虫接近食物来源并最终进食时的减速,其变化机制可以保证GOA算法不过快地收敛到目标，避免陷入局部优化，并在优化的最后一步实现快速收敛。

优化算法的另一核心是适应度函数的构造。本文选取对周期故障冲击较为敏感的相关峭度系数作为适应度函数，相关峭度函数表达式如下

(13)

式中，当移位数M和冲击周期T分别为1、0时，相关峭度表达式等价于峭度表达式。因此，相关峭度与峭度的最大差别在于，峭度系数对单个冲击信号更为敏感，而相关峭度系数在考虑到冲击周期信息后，对周期性冲击更为敏感。因此，相关峭度在利用故障特征高峭度指标特点的同时，兼顾了旋转机械故障的周期性特征，更适用于故障特征提取。

2 方法步骤

本文滚动轴承故障特征提取流程图如图1所示。方面步骤如下。

步骤1：参数初始化。在GOA优化算法中，设定模态数k的搜索空间为[2,10]，二次惩罚项α的搜索空间为[2 000,10 000]，参数搜索代理个数为15，迭代次数为10。相关峭度中，移位数M的选取一般不大于8。

步骤2：经GOA优化算法最终确定优化参数[K,α],将优化后得到的参数代入VMD算法中，对目标振动信号进行分解，得到K个模态分量。

步骤3：计算各模态分量的相关峭度值，选取具有最大相关峭度的模态分量作为包含故障信息最多的分量，用于之后的故障特征提取。

图1 轴承故障特征提取流程图

步骤4：对步骤3中选取的最佳模态分量进行包络解调分析，与轴承实际故障特征频率进行比对，进行故障特征提取。

3 单一故障仿真信号验证

为验证本文算法的有效性，使用单一周期冲击序列，并向其中加入强烈白噪声来模拟滚动轴承的单点故障[13]，仿真信号如下

(14)

h(t)=exp(-Ct)sin(2πfnt)

(15)

在故障模型中，采样频率fs=12 000 Hz，采样时间t=0.5 s，故障频率fc=50 Hz，共振频率fn=2 000 Hz，求和部分为周期脉冲成分，A为均值为1.1，方差为0.5的随机数，模拟信号中加入τi来代表滚动体经过故障处因载荷角等因素变化引起的周期波动，正弦指数衰减函数h(t)用于模拟系统的脉冲响应，衰减系数C=200，仿真信号中加入信噪比为-10 dB的高斯白噪声n(t)。仿真信号时域图及包络谱如图2、图3所示。由图2中可以看出，原本周期冲击信号特征在噪声加入后被掩盖，不能观察到周期特征。包络谱中也不能发现故障特征倍频谱线。利用本文提出的特征提取方案，首先进行优化变分模态分解，适应度值优化过程如图4所示。

经GOA优化算法最终确定的k和α的值为[8,7 411]，代入VMD中分解结果如图5所示。经VMD处理后得到8个模态分量，此时问题的关键在于模态分量的选择问题。按照以往EMD及EEMD等分解方法对模态分量的选择标准，一般以峭度(Kurtosis, kt)或相关系数(correlation coefficient, co)作为指示模态分量中含有故障信息的指标。本文提出以含有周期性参数的相关峭度作为筛选模态分量的指标。分别计算各分量的相关峭度指标。为进行对比，同时计算峭度指标及相关系数指标，并进行归一化，各指标如图6所示。

(a) 时域冲击信号

(b) 加噪声后时域信号

图3 仿真信号包络谱

图4 适应度值迭代过程

图6中可以看出，相关峭度和相关系数指标指向第四模态分量，峭度指标指向第七模态分量，分别提取第四和第七模态分量，进行包络解调，包络谱如图7所示。图7(a)中，第四模态分量时域图中出现明显的周期冲击特征，经包络解调后，包络谱中出现明显的前三倍频特征谱线。第七模态分量时域图中并不能明显观察到周期冲击特征，包络谱中也不能发现故障特征谱线。结果证明相关峭度指标正确指示出含有较多故障特征信息的模态分量。

图6 三种指标对比

(a) 第四模态分量时域图

(b) 第四模态分量包络谱

(c) 第七模态分量时域图

(d) 第七模态分量包络谱

进一步对仿真信号进行EMD分解做对比验证，分解结果如图8所示。EMD分解得到11个模态分量和1个残余分量，利用峭度准则选取峭度较大的第2、3、5模态分量进行重构，重构信号及包络谱如图9所示。从重构信号时域图中不能发现明显的周期冲击特征，包络谱中也没有出现与故障相关的频率谱线，提取效果远不如图7(a)中的方法。

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

(b) EMD分解重构信号包络谱

4 单一故障实测信号验证

4.1 外圈故障信号

为进一步验证提出的故障特征提取方法在应用到实际轴承故障信号时的表现情况，采用美国凯斯西储大学电气工程实验室的滚动轴承数据[14-15]。在该实验中，轴承故障是电火花加工技术布置的单点故障，故障损伤直径有0.177 8 mm、0.355 6 mm、和0.533 4 mm三种类型，轴承型号为SKF公司的6203-2RS和6205-2RS深沟球轴承。文中选择其中转速为1 797 r/min，采样频率为12 kHz，试验台驱动端轴承的外圈故障数据来验证该文所提方法的有效性。外圈故障特征频率频率根据轴承参数计算fo=107.4 Hz，转频fr=29.95 Hz。选取故障尺寸最小的0.177 8 mm，外圈故障信号时域图及包络谱如图10所示。由图10(a)可以看出，时域图中不能发现明显的周期冲击特征，包络谱图10(b)中也未能发现与外圈故障特征频率吻合的倍频谱线。利用本文提出的故障特征提取方案，首先进行GOA参数优化，适应度函数优化过程如图11所示。

(a) 外故障信号时域图

(b) 外圈故障信号包络谱

图11 适应度值迭代过程

GOA优化算法最终确定的参数集合为[6,7 273]，代入VMD中进行模态分解，分解结果及筛选指标如图12、13所示。由VMD分解结果及提取的三种筛选指标中可以看出，相关峭度、相关系数、峭度指标分别指向第2、3、6模态分量，提取三个模态分量并进行包络解调，结果如图14所示。由各指标对应的分量包络谱中可以看出，相关峭度指标指向的分量明显含有较多的故障信息，而具有较大峭度指标及相关系数指标的模态分量包络谱中并没有发现明显的故障特征谱线，对比说明了相关峭度指标在正确筛选含有故障成分的模态分量的优势，同时也证明了故障特征提取方法的有效性。

图13 三种指标对比

(a) 相关峭度指标对应的分量包络谱

(b) 相关系数对应的分量包络谱

(c) 峭度系数对应的分量包络谱

对比进行EMD处理外圈故障实测信号，经EMD处理后得到12个模态分量，根据峭度指标筛选出第2、3、5、7模态分量进行重构，重构信号的包络谱如图15所示。从包络谱中我们可以看出，EMD处理后重构信号包络谱中的仅能看出数条与转频相关的倍频谱线，并不能发现故障特征，故障提取效果较差。

图15 重构信号包络谱

4.2 内圈故障信号

同样采用凯斯西储大学数据库中的轴承内圈故障实测数据，选取转速为1 797 r/min,故障尺寸为 0.533 4 mm，计算得内圈故障频率fi=162.2 Hz，故障信号原始时域波形如图16(a)所示。为适应于强噪声背景特征提取的特点，在信号中加入-10 dB的白噪声，内圈故障时域图如图16(b)所示。时域图中，内圈故障冲击特征被强烈的噪声淹没，看不出周期特征，包络谱中也完全不能发现任何内圈故障特征谱线。如图17所示为内圈故障包络谱。

(a) 内圈故障信号

(b) 强噪声加入后信号

图17 内圈故障包络谱

根据GOA优化算法搜索得到的适应度值的优化过程如图18所示，并且最终得到的参数为[7,6 619]，代入变分模态分解算法进行分解得到7个模态分量，各模态分量的时域图如图19所示。计算各模态分量的三种筛选指标如图20所示。从图中看出，三种筛选指标中，峭度指标指向第2模态分量，相关系数和相关峭度指标指向第5模态分量。分别提取各模态分量包络谱进行对比，包络谱如图21所示。

图18 适应度值迭代过程

图20 三种筛选指标对比

由包络谱中可以看出，峭度指标指向的模态分量包络谱中并未发现任何故障信息。与之对比，相关峭度和相关系数指标指向的模态分量中明显发现内圈故障特征，包络谱中除发现转频及前两倍内圈故障特征频率外，也能明显观察到内圈故障频率谱线两侧的转频调制谱线。对比进行EMD方法提取故障特征信息，分解得到11个模态分量，提取峭度较大的第6、7、8模态分量进行重构，重构信号的包络谱如图22所示。由重构包络谱可以看出，并未发现与故障相关的特征谱线，故障特征提取效果不如本文方法。

(a) 峭度指标对应分量包络谱

(b) 相关峭度、相关系数指标对应的包络谱

图22 重构信号包络谱

5 复合故障仿真信号验证

经上述单一故障仿真信号及实测信号的验证，优化变分模态分解算法合理且在故障特征提取中有着较为明显的效果。另外，工程实际生产中也同样存在着复合故障现象且为目前学者研究的热点问题。为拓宽论文方法的应用范围，同时进一步加深验证方案的有效性，采用轴承内圈和外圈的复合故障仿真信号对论文方案进行验证。

在复合故障仿真信号中，外圈故障信号采用与第3节相同的模型，内圈故障信号的数学模型如下

(16)

h1(t)=exp(-C1t)cos(2πfnit)

(17)

A1=Asin(2πfrt)

(18)

在内圈故障仿真信号模型中，采样频率fs=12 000 Hz，采样时间t=0.5 s，故障频率fci=115 Hz，共振频率fni=3 500 Hz，衰减系数C1=1 100，转频fr=50 Hz；外圈故障仿真信号模型中，故障频率fco=65 Hz，共振频率fno=2 500 Hz，衰减系数C=900。最后在复合故障信号中加入信噪比为-10 dB的白噪声，内外圈故障仿真信号以及加入强烈噪声的复合故障信号时域图如图23所示。由图23(c)看出，在加入噪声后，内外圈周期冲击特征被完全淹没。从图24的复合故障信号包络谱中看出，内外圈故障特征频率除受到强烈噪声干扰，除能勉强观察到个别特征谱线外，内外圈故障特征频率及多倍频成分互相混杂在一起，难以对滚动轴承复合故障进行清晰且有效的诊断。

(a) 外圈故障仿真信号

(b) 内圈故障仿真信号

(c) 加入-10 dB噪声的复合故障信号

图24 复合故障信号包络谱

利用文章提出的优化变分模态分解方法进行内外圈复合故障特征提取。对于复合故障信号中的内圈故障，经GOA优化算法确定的VMD参数为[5,3 462]，代入VMD算法中分解得到5个模态分量，计算相关峭度指标后指向第四模态分量，提取第四模态分量进行包络解调分析，结果如图25所示。由图25(a)中已经可以看出明显的周期冲击成分，图25(b)包络谱中更是发现了一阶转频，前三阶内圈故障特征频率，并且可以明显看出受转频调制的边频带，内圈故障提取效果明显。

(a) 第四模态分量时域图

(b) 第四模态分量包络谱

对于复合故障中的外圈故障，经GOA算法确定的VMD参数为[8,6 792]，代入VMD算法中进行分解得到8个模态分量，计算各模态分量相关峭度指标后指向第五模态分量，提取第五模态分量进行包络解调分析，结果如图26所示。

(a) 第五模态分量时域图

(b) 第五模态分量包络谱

由图26可知，时域图中出现明显的周期冲击特征，包络谱中出现前四阶的外圈故障特征频率，算法对外圈故障的提取效果同样明显。

6 结 论

(1) 滚动轴承故障的早期阶段，故障冲击信号微弱，容易淹没在较强的背景噪声中。相关峭度指标优化变分模态分解的方法解决了变分模态分解在参数选取上存在的问题，能够较好地将隐藏在强噪声中的故障信息分离出来。

(2) 提取故障特征的关键之一在于能找到对故障特征敏感有效的指标。在筛选含有较丰富故障特征的模态分量中，由于相关峭度指标对周期冲击特征更有针对性，因此更适合作为模态分量筛选指标。

(3) 利用单一故障仿真信号、单一故障内外圈实测信号以及内外圈复合故障仿真信号对论文方案进行验证。分析结果表明，基于相关峭度优化变分模态分解的滚动轴承故障特征提取方法在强噪声背景下特征提取中具有较为明显的优势。