基于Gibbs抽样的轴承故障诊断方法

王 岩，罗 倩，邓 辉

(北京信息科技大学 信息与通信工程学院，北京 100101)(*通信作者电子邮箱wyleo7@qq.com)

0 引言

在旋转机械中，滚动轴承是易于损坏的零部件之一。据有关资料统计，旋转机械的故障中振动故障占70%，而30%的振动故障是由滚动轴承故障引起的[1]，因此滚动轴承故障诊断的理论和应用研究一直是一个重点。

随着大数据和人工智能时代的到来，越来越多的机器学习算法被用于轴承故障诊断中的模式识别，如Ali等[2]采用神经网络，焦卫东等[3]使用支持向量机(Support Vector Machine, SVM)进行轴承故障诊断，但它们或多或少地存在一些缺陷，比如支持向量机中针对每一种轴承故障都设置一个判别函数进行诊断，造成输入空间中存在无法分类的区域，即诊断片面性的问题[4]。为解决此问题，本文提出了一种新的基于Gibbs抽样的轴承故障诊断算法，该方法将概率统计学应用到轴承故障诊断中。Gibbs抽样是马尔可夫链蒙特卡罗(Markov Chain Monte Carlo， MCMC)理论的一个分支，用来获取一系列近似等于指定多维概率分布(比如2个或者多个随机变量的联合概率分布)观察样本的算法。使用Gibbs抽样建立每一种轴承状态的分布模型，针对一个待诊断的样本，分别计算出属于每一种轴承状态的概率，选取最大的判定为该类别，从而达到故障诊断的目的，解决了轴承故障诊断中存在的诊断片面性问题。

基于Gibbs抽样的滚动轴承故障诊断思路如下：首先对轴承故障振动信号进行局部特征尺度分解(Local Characteristic scale Decomposition， LCD)得到若干个内禀尺度分量(Intrinsic Scale Component， ISC)，再对轴承故障振动信号和ISC分别提取时域特征，对所有时域特征使用特征的类间标准差和类内标准差的比值筛选出敏感特征并组成特征集。针对不同种类的轴承故障，使用基于Gibbs抽样算法训练特征集产生不同的多维高斯分布模型。对于待诊断的故障数据，分别计算在每个模型中的概率密度，通过后验分析得到概率，选取概率最大的结果将待诊断数据诊断为该故障类别。

1 数据预处理和特征提取

1.1 数据LCD预处理

准确提取轴承振动信号特征是轴承故障诊断中重要的一步，如果直接从这些非平稳或非线性振动信号中提取特征势必影响故障诊断的效果，因此特征提取前必须进行信号预处理。本文采用局部特征尺度分解进行信号预处理，它能够同时在时域和频域提供非平稳信号的局部化信息，且避免了经验模态分解(Empirical Mode Decomposition, EMD)的频率混淆问题[5]和局部均值分解(Local Mean Decomposition， LMD)的信号突变问题[6]，所以它对非平稳信号处理效果比较好。

LCD是一种基于极值点的局部尺度参数自适应的非平稳信号分解方法。该方法以任意相邻的两个极值点为跨度，通过分段的形式对信号进行线性变换来构造基线函数，从而将原始信号自适应地分解为若干个ISC分量。其具体步骤为:将原始信号通过两个循环分解为具有物理意义的内禀尺度分量ISC和一个剩余项。第一个循环是寻找基线信号，并经循环迭代求出ISC分量的过程，ISC分量应满足ISC分量条件；第二个循环是将原始信号减掉已获得的ISC分量后的剩余信号作为原始信号再次经循环迭代得到ISC分量。重复上述步骤直至所得剩余项为一个单调函数为止[7]。对于任一时间序列，通过LCD分解可以得到若干个ISC和一个剩余项。LCD分解本质是对非平稳信号进行平稳化处理，将信号中不同尺度的波动或趋势逐级分解。因为信号的能量主要集中在前几个ISC分量中，因此本文只取前4个ISC分量。

1.2 特征提取

在LCD预处理之后，对原始振动信号和4个ISC分量分别进行时域特征提取。故障特征可以分为无量纲特征和有量纲特征两大类[8-9]。无量纲指标因其不受工况影响的特点，对剥落、压痕等故障有很好的诊断效果；有量纲指标对磨损类故障诊断效果比较好。因此，为了能够更加全面、准确地诊断故障，本文对原始信号和4个ISC信号分别提取了6个无量纲特征(偏斜度、峭度、峰值、波形指标、脉冲、裕度)和5个有量纲特征(有效值、均值、标准差、方根幅值、最大值)，共得到55个特征。

1.3 特征筛选方法

如果使用上述55个特征作为一个特征集，计算量过大且诊断效果不一定最佳，所以需要从中挑选出敏感特征用于模型训练。本文使用了基于故障类间、类内标准差比值的特征筛选方法。

类内标准差计算如下：

(1)

类间标准差计算如下：

(2)

对于某一特征，其类间标准差与类内标准差的比值，称为该特征的敏感度指标。计算公式如下：

χ=Si/So

(3)

其中：χ表示敏感度指标，特征对不同状态的区别度与敏感度成正相关，某一特征的χ值越大，该特征对不同状态的区别度越高。

2 基于Gibbs抽样的轴承故障诊断

本文只涉及到轴承运行的4种状态：内圈故障状态、外圈故障状态、滚动体故障状态和正常状态。基于Gibbs抽样的轴承故障诊断基本思路是：从上述55个特征中选取排名靠前的若干敏感特征组成特征集，针对轴承的每一个状态，对这些特征建立多维概率分布，使用Gibbs算法对多维概率分布进行参数估计，最终得到4个多维概率分布模型。对于一个待诊断轴承数据，分别计算在4种模型中的概率密度，通过后验分析得到概率，选取最大概率判定待诊断数据为该轴承状态类别。

2.1 特征值分布

建立特征的多维概率分布，首先需要分析每个特征在不同状态下的一维分布，然后构建出合适的多维概率分布。生成轴承不同状态下各特征的分布直方图，通过直观观察发现绝大多数直方图服从高斯分布。进一步，使用Quantile-Quantile正态检验方法[10]和Michael拟合优度检验方法[11]检验各特征分布，发现除信号峭度特征在任何轴承状态下不满足高斯分布外，其他特征都服从高斯分布，这些特征都落在90%接受区间内，可以认为这些特征都服从高斯分布，因此，针对轴承的每一个状态，可以建立特征集的多维高斯分布。

2.2 多维高斯分布参数估计

使用若干个排名靠前的特征组成特征集构建多维高斯分布，对于轴承故障的第j个状态的多维高斯分布模型Vj的参数有均值μj、协方差矩阵Σj，对这两个参数进行参数估计。高斯分布均值μj的先验概率服从高斯分布，协方差矩阵Σj的先验概率服从逆Wishart分布，通过对参数先验分布进行Gibbs抽样，可以得到这两个参数具体的值[12]。故首先对每个多维高斯分布分别建立参数的先验分布。

1)从逆Wishart先验分布中随机抽样协方差矩阵Σj，下标j表示轴承第j个状态：

(4)

(5)

2)从高斯先验分布中抽样协方差矩阵μj，如式(6)：

(6)

(7)

第j个轴承状态的多维高斯分布参数估计的Gibbs算法通过将上述两个步骤运算N次，可以分别得到每个参数的N个抽样结果，最终μj和Σj参数值由式(8)～(9)得到，其中舍弃了前N/2的抽样结果[13]。

(8)

(9)

轴承故障第j个状态由均值μj和协方差矩阵Σj构成的多维高斯分布模型为Vj：

(10)

其中：X是待诊断数据根据训练特征集提取的时域特征；D是特征集中包含特征的个数，即多维高斯分布的维数；uj表示X在第j个轴承状态模型中计算得到概率密度。

得到概率密度后，需要进行后验概率分析。在轴承故障诊断之前，假设该轴承被诊断为4种状态的先验概率Zj是相等的，即Zj=0.25。由贝叶斯定理[14]可得，X被诊断为第j个状态的后验概率pj为：

(11)

2.3 轴承故障诊断

使用轴承不同状态的振动信号样本数据经由LCD并提取时域特征后，选择D个敏感特征组成特征集，训练产生4个多维高斯分布模型。对于待诊断的数据，先根据训练特征集提取的时域特征分别计算式(10)得到概率密度u1、u2、u3和u4，再使用式(11)分别计算概率p1、p2、p3和p4，选取最大概率诊断为该待诊断数据的类别。

3 仿真及实测数据诊断结果分析

使用了凯斯西储大学的滚动轴承振动数据进行仿真(数据来源：http://csegroups.case.edu/Bearingdatacenter/pages/download-datafile)。振动信号由驱动端加速度传感器采集，轴承型号为SKF6205，电机转速为1 772 r/s，使用电火花加工技术在驱动端轴承上布置了单点故障，内圈、滚动体及外圈故障直径都为0.007 inch(1 inch=2.54 cm)，外圈故障位于6点钟位置，采样频率为12 kHz，实际故障频率：内圈为159.93 Hz，滚动体为139.20 Hz，外圈为105.88 Hz。训练模型时，将轴承振动信号若干个采样点作为一组，从凯斯西储大学原始数据中提取一共240组组成训练样本，其中轴承4种状态各占60组。测试模型时，从原始数据中另外提取一共240组组成测试数据，其中轴承4种状态也各占60组。将训练数据的时域特征按照敏感度从高到低排列，其中排名前10的特征如表1所示，逐个将特征加入特征集进行模型训练，可以得到55组不同特征集的模型。测试数据按照训练特征集提取特征，分别计算其在55组不同模型中的诊断正确率，可以得到特征集取不同数量的特征时该算法的诊断正确率。

表1 特征敏感度排名

从图1中可以看出，当特征集的特征个数增加到2时，诊断正确率就达到100%。图1中的诊断正确率随着特征数的增加呈下降趋势，因为特征按照本文提到的特征筛选排名方法，逐个顺序加入特征集并计算正确率。按此排名顺序，特征对于区分不同轴承故障状态的敏感度随着特征数增加而下降，第20个特征之后加入的特征敏感度较差，会对整体诊断正确率产生影响，因此图1中诊断正确率随着特征数的增加而下降，也验证了特征筛选排名方法的有效性。

图1 诊断正确率与特征数关系

将240组测试集数据分别代入使用4个特征作为特征集的训练模型中，按照诊断结果的标签对4维数据使用了多维标度法(Multi-Dimensional Scaling， MDS)[15]映射在二维空间中，如图2所示。

图2 高维测试集数据的二维映射(240组测试集数据)

从图2中可以看出，4种状态数据特征分布较为鲜明，说明本文提出的算法可以很好地判别4种不同的轴承状态。

对前4个特征组成的特征集进一步分析，使用上述240组训练集产生轴承4种不同状态的模型：内圈故障V1、外圈故障V2、滚动体故障V3和正常状态V4。每个模型的参数μj和Σj计算结果列于表2，在图3中使用多维标度法将不同模型的参数μj映射在二维空间中。

表2 多维高斯分布模型参数

图3 高维高斯分布均值的二维映射

为了验证本文提出的算法优于基于支持向量机的轴承故障诊断算法，使用凯斯西储大学轴承数据中采样频率为48 kHz的数据，振动信号由风扇端加速度传感器采集，轴承型号为SKF6205，电机转速为1 797 r/s，使用电火花加工技术在驱动端轴承上布置了单点故障，内圈、滚动体及外圈故障直径都为0.021英寸，外圈故障位于3点钟位置，实际故障频率：内圈为162.20 Hz，滚动体为141.09 Hz，外圈为107.30 Hz。使用本文中的特征提取、筛选及排名的方法，在轴承故障诊断中分别使用了基于Gibbs抽样的诊断方法和基于支持向量机的诊断方法。原始数据中缺少采样频率为48 kHz的轴承正常状态数据，因此只建立了有关内圈故障、外圈故障及滚动体故障的模型，从凯斯西储大学原始数据中提取一共180组数据组成训练样本，其中轴承3种状态各占60组。测试模型时，从原始数据中另外提取一共180组组成测试数据，轴承3种状态也各占60组。特征按照本文中的方法进行排名并依次顺序加入特征集，计算每个特征集有关模型的诊断正确率。

如图4所示：在特征数比较少时，两种诊断方法都能达到最高100%的诊断正确率；随着特征数的增加，特征对轴承不同状态的敏感度降低，二者的诊断正确率都下降，但基于Gibbs抽样诊断算法下降的幅度明显比支持向量机小；特别在特征数达到43时，基于Gibbs抽样方法诊断正确率为82.8%，而基于支持向量机方法诊断正确率为71.7%。与基于支持向量机诊断算法相比，基于Gibbs抽样诊断算法正确率提升了11.1个百分点；当特征数大于43后，两种算法的诊断正确率差值进一步增大。说明Gibbs抽样诊断算法对于敏感度低的特征诊断效果好于支持向量机，特别对于高维复杂数据的离群点及支持向量机输入空间中存在无法分类的区域[4]，Gibbs抽样诊断算法通过概率计算依然可以实现对轴承故障的诊断。

图4 两种算法的诊断正确率比较

凯斯西储大学的数据是实验室条件下得到的数据，与真实的数据还有一定差别，使用该方法对中国铁路某局实际的滚动轴承数据进行了诊断。数据来自中国铁路某局甲2014年某月。该数据只有轴承故障和正常两个状态，并没有对故障进行细分，数据是经过共振解调[15]预处理过的。图5是其振动信号200个抽样点的时间序列波形。

图5 轴承振动信号时间序列

使用60组故障数据和60组正常数据作为训练样本，另外60组故障数据和60组正常数据作为测试数据。图6是该局使用本文提出的诊断方法在不同数量特征集下的诊断正确率，由图可知取14个特征组成特征集时诊断正确率最大，为95.8%。

图6 某局甲的故障诊断正确率

使用14个特征组成特征集对测试集数据进行诊断，诊断结果的标签使用了多维标度法(Multi-Dimensional Scaling， MDS)[15]显示二维空间中，如图7所示。从实际数据仿真结果可以看出，本文方法依然有较高的故障诊断正确率。

图7 高维测试集数据的二维映射(14个特征组成特征集)

4 结语

本文提出了一种基于Gibbs抽样的轴承故障诊断算法，通过理论分析及数据仿真结果分析表明:

1)本文算法将概率统计学运用到轴承故障诊断中，采用概率模型代替支持向量机中的判别函数进行故障诊断，参考每种模型概率密度的计算结果，通过后验分析得到概率，以判别概率最大为准则，对于支持向量机输入空间中无法分类的区域[4]也能有效判别，避免了判别偶然性和片面性的出现，使得识别结果更加客观。与基于SVM的轴承诊断方法相比，在特征数为43时诊断正确率提升了11.1个百分点。

2)提出了将Gibbs抽样算法应用于滚动轴承故障诊断，可以实现从信号处理到分类的智能故障诊断。

3)本文分别使用凯斯西储大学实验室条件下的数据和中国铁路机车轴承的实际采集数据进行了算法验证，结果表明本文方法能够有效地诊断轴承故障状态。