张星辉， 康建设， 赵劲松,2， 肖 雷， 曹端超

(1.军械工程学院，石家庄 050003；2.军事交通学院，天津 300161;3.重庆大学，重庆 400030)

有效的退化状态识别是设备健康管理的关键步骤。齿轮及齿轮箱是机械设备中一种必不可少的部件，对其进行状态监测有着重要意义。国内外学者对齿轮箱不同故障模式的识别已经进行了大量的研究，取得了很多成果。而研究单一故障模式从初始发生到失效整个过程的退化规律对维修活动规划有很大意义。齿面磨损是齿轮的常见退化模式，通常会随着运行时间的流逝而逐渐严重，直到齿轮失效。齿面磨损一般不便直接测量，但可通过外部测量设备来确定其磨损状态。目前，隐马尔可夫族模型(Hidden Markov Models，HMMs)[1-4]、支持向量机[5-6]、人工神经网络等模型[7-8]都可用来描述由设备外部测量结果反映内部退化状态的过程，并被广泛应用于各种设备的退化状态识别。贝叶斯网络是表示不确定性专家知识和推理的一种流行方法，非常适合于表达诊断决策问题。回顾贝叶斯网络在设备故障诊断[9-12]中的应用，多是对于整个复杂设备系统进行的简单故障推理，而没有对具体部件故障与外部观测信息之间的关系进行研究。因此，本文构建了混合高斯输出贝叶斯信念网络模型(Mixture of Gaussians Bayesian Belief Network，MoG-BBN)用于磨损状态识别。模型应用变量消元(Variable Elimination，VE)算法求取后验概率，对后验概率需要的参数应用期望最大化(Expectation Maximization，EM)算法进行估计。由于标准的EM算法容易局部收敛，因此本文对算法进行了改进。最后，齿轮磨损实验数据用来对模型进行验证。

1 MoG-BBN基本原理

贝叶斯信念网络(Bayesian Belief Network，BBN)是一种概率推理网络，它以图形节点表示随机变量，节点之间的有向箭头表示随机变量之间的因果关系。BBN是一个在不确定条件下强有力的知识表达和推理方法。BBN中节点所代表随机变量的取值可以是连续的，也可以是离散的。连续变量可以服从任意分布，离散变量的取值可以是两个或多个。本文构建的MoG-BBN模型，其网络结构如图1所示。

该模型可以用如下参数描述：

(1)X表示隐状态，其状态取值分别为1,2,…,K，其代表齿轮的磨损状态；

(2)M表示混合数，其取值范围为1,2,…,b，其代表齿轮箱的某个磨损状态按照第j个高斯分布混合产生观测值；

图1 贝叶斯信念网络

在图1所示的贝叶斯网络中，X是根节点(M和Y的父节点)，M是中间节点(X的子节点，Y的父节点)，Y是叶节点(X和M的子节点)。

2 VE-EM算法

2.1 VE算法

根据贝叶斯信念网链规则得：

P(X,Y,M)=P(X)P(M|X)P(Y|X,M)

(1)

根据条件概率公式得：

(2)

式(2)右侧的分子可表示为：

(3)

式(2)右侧的分母可表示为：

(4)

结合式(1)、(3)、(4)，式(2)可进一步表示为：

(5)

2.2 EM算法

假设训练数据集为{y(1),…,y(n)}，由式(5)可以最终求得数据所对应的故障模式，因此必须求得P(X)，P(M|X)和P(Y|X,M)，这三个量分别对应四个参数：

(1)π：隐状态(设备退化状态)概率分布， πi=P(X=i)，i∈{1,2,…,K}；

(2)C：隐状态混合系数，Cij=P(M=j|X=i)，i∈{1,2,…,a}，j∈{1,2,…,b}；

(3)μ：隐状态产生高斯分布的均值；

(4)∑：隐状态产生高斯分布的方差。

这四个参数可以表示为：λ=(π,C,μ,∑)。

EM算法就是求使得P(Y|λ)最大化的参数λ，然后再通过VE算法求取P(X|Y)。详细的推导如下：

(6)

(7)

其中，q(X)为引入的未知分布，并给出它的约束条件。根据条件概率公式和贝叶斯公式，式(7)可表示为[13]：

(8)

其中的两项分别为：

(9)

(10)

根据Jensen不等式，式(10)可表示为：

(11)

由此可得：lnP(Y|λ)≥L(q,λ)。并且当且仅当q(X)=P(X|Y,λ)时，lnP(Y|λ)=L(q,λ)。其它情况下，L(q,λ)< lnP(Y|λ)。因此，最大化lnP(Y|λ)等同于最大化L(q,λ)。而L(q,λ)只有在q(X)=P(X|Y,λ)时取得最大化，将q(X)=P(X|Y,λt-1)代入式(9)得：

Q(λt,λt-1)+const

(12)

其中：

(13)

(14)

E步：

对所有的(i,j,l)

(15)

根据贝叶斯公式，式(15)可表示为：

P(X(l)=i,M(l)=j|y(l),π,C,μ,∑)=

(16)

M步：更新参数

(17)

(18)

(19)

(20)

重复E步和M步，直到迭代一定的步数或者相邻两次迭代满足|lnP(Y|λt)-lnP(Y|λt-1)|<ξ，ξ取值一般在区间[10-3,10-5]。待获得估计后的参数，就可以通过变量消元算法求取各种状态产生观测值的概率。模型的迭代次数、ξ可以先取一个较大范围的值，如果模型在很小的范围内收敛，那么，在下次训练时，再进一步缩小范围。反之，应加大取值的范围。

2.3 改进后EM算法

由于2.2节的EM算法很容易收敛于局部最优值，而不是全局最优值，因此需要对EM算法进行改进。该改进算法的主要思想是：在每次迭代得到参数 后，利用该参数随机生成一批数据后，对参数再次进行估计。这样，相当于每次迭代都进行了一次参数寻优，避免了参数很快的收敛于局部最优值。其具体过程可以用以下几步来描述：

第1步：(模型参数初始化)以均匀分布随机产生隐状态概率分布和隐状态混合系数，并分别对它们进行归一化。以一组训练数据(观测值)的均值和方差作为隐状态产生高斯分布的均值和方差。

第2步：用式(21)对每一步迭代得到的参数λi进行评价:

(21)

其中，P(Yl|λi)表示模型参数λi产生第l组训练数据的概率。

其中，0<φ(i)<1。

3 曲线距离分析

为了更好地识别齿轮的磨损状态，往往会从多个传感器采集的振动信号中提取多种特征，特征的维数会多达几十甚至上百，这些特征间存在一定的相关性，若把这些特征直接输入识别模型将会导致计算量急剧增大且计算结果不精确。因此，需要对特征向量进行降维，从而提取出少数关键的特征。

主成分分析和独立成分分析是常用的降维方法，它们都假设特征之间的关系是线性的，是线性变换方法。而对于反映齿轮磨损状态的特征向量而言，特征之间往往呈现高度的非线性关系，因此，应用非线性变换方法对其降维更为合理。

CDA[14-15]是用于非线性特征的降维方法，被广泛的应用于图像和语音识别中。该算法通过最小化式(22)来达到降维的目的：

(22)

通常，采用文献[15]提出的递减函数，表示为：

(23)

ECDA可以通过随机梯度下降法求得：

i≠j

(24)

其中：α(t)为t的递减函数。CDA的详细计算过程和步骤可见文献[15]。

4 基于CDA和MoG-BBN的齿轮磨损识别框架

图2 齿轮磨损状态识别框架

5 实验及数据分析

5.1 实验设置

实验齿轮箱组成如图3所示，实验所用齿轮箱型号为ZD10；动力源为电磁调速电机，型号为YCT180-4A；水冷磁粉制动器为齿轮箱提供载荷，型号为FZJ-5。该实验是对输出轴81齿齿轮预置单齿磨损故障，从齿顶量取的磨损量分别为0.2 mm、0.3 mm、0.46 mm和0.6 mm。运行工况分别取200 r/min、400 r/min、600 r/min、800 r/min、1 000 r/min及升速、降速条件下恒载和变载运行。采样频率为20 kHz，采样时间为6 s，每种工况采集20组数据。

图3 齿轮箱结构及传感器位置示意图

图中，①、②、③、④、⑥、⑦号传感器位于轴承座上方，⑤和⑧号传感器位于轴承座下方。定义齿轮的四种磨损状态为：磨损状态1(0.2 mm)、磨损状态2(0.3 mm)、磨损状态3(0.46 mm)和磨损状态4(0.6 mm)。选取转速200 r/min、400 r/min、600 r/min、800 r/min、1 000 r/min，负载0.6 A(磁粉制动器控制电流)工况下的数据作为分析对象。

5.2 磨损状态识别

5.2.1 特征提取

(1)时域特征

该实验提取的时域特征有：均值、标准差、均方根值、峰值、偏斜度、峭度。

(2)齿轮统计特征参数

统计特征参数是NASA技术报告中提出的专门检测齿轮故障的特征[16-18]。这些特征分别定义如下：

(25)

式中：PPx表示时域信号x峰峰值的最大值，Pn表示齿轮啮合频率n阶谐波的幅值，H表示最大的谐波阶数。

(26)

式中：FM4的实质就是取差分信号的峭度值，d表示差分信号(如图4所示)。

(27)

式中：M′表示前M′个时间点测取的是正常齿轮的信号，下同。

(28)

式中：M6A是对差分信号的6阶统计量。

(29)

(30)

式中：M6A是对差分信号的8阶统计量。

(31)

(32)

式中：r表示剩余信号(如图4所示)。

(33)

(34)

式中：s表示包络信号。

(35)

(36)

(37)

图4 特征提取流程

(3)小波包频带能量特征

应用matlab小波工具箱对信号进行小波包3层分解，将第3层小波包分解系数的能量作为齿轮磨损特征，共8个，分别表示0～1.25 kHz，1.25～2.5 kHz，2.5～3.75 kHz，3.75～5 kHz，5～6.25 kHz，6.25～7.5 kHz，7.5～8.75 kHz和8.75～10 kHz。

综上所述，每个通道可提取27个特征，8个通道共216个特征。

5.2.2 特征降维及归一化

经CDA降维后的特征维数为30，每种工况条件下每个磨损状态选取10组数据用来训练，10组数据用来测试。因此，每种工况条件下的训练和测试数据都为40组(磨损状态1、磨损状态2、磨损状态3和磨损状态4各10组)。

对训练集和测试集进行归一化预处理，采用的归一化映射为：

(38)

式中，x,y∈Rn，xmin=min(x),xmax=max(x)。归一化的效果是原始数据被规整到[0,1]范围内，即yi∈[0,1],i=1,2,…,n。该归一化方式称为[0,1]区间归一化。

5.2.3 状态识别结果及分析

应用归一化后的训练和测试数据按照图2所示的识别流程对模型进行验证。模型最大迭代次数为100，ξ取值为1×10-5，模型混合数为2。转速200 r/min条件下，10组磨损状态1的测试数据识别结果如表1所示。

表1 转速200 r/min条件下10组磨损状态1测试数据识别结果

从表中可以看出，测试数据处于磨损状态1的概率最大，识别结果与测试数据的实际状态相一致。验证了该方法的有效性。在训练过程中，模型迭代15步后即找到了最优参数，也就是ξ的值已经小于设定的1×10-5。对于模型混合数的选取，目前还没有跟好的方法，下一步将对模型混合数进行优化研究。

由于篇幅关系，所有工况的识别结果都用图来表示，如图5、6、7、8和9所示。图中，数据的排列为磨损状态1(10组)、磨损状态2(10组)、磨损状态3(10组)和磨损状态4(10组)。概率最大的磨损状态即为识别结果。从图中可以看出，除转速800 r/min和1 000 r/min条件下各有一组数据识别错误外，其它测试数据的识别结果与实际结果都一致。总的识别正确率可以达到99%。在MoG-BBN模型的实际应用中，模型参数对识别结果有很大的影响。如果由训练数据估计得到的模型参数与代表数据的真实参数相差太远，也就是说，通过训练没有得到最优的参数。那么，模型的识别效果会很差。反之，模型识别效果会很好。对实验数据的分析表明，经过改进参数求解算法后的MoG-BBN识别效果非常好。

图5 转速200 r/min条件下MoG-BBN识别结果

图6 转速400 r/min条件下MoG-BBN识别结果

图7 转速600 r/min条件下MoG-BBN识别结果

图8 转速800 r/min条件下MoG-BBN识别结果

图9 转速1 000 r/min条件下MoG-BBN识别结果

为了更进一步验证模型的识别效果，应用HMM对实验数据进行处理并进行对比。对比结果如表2所示。从表2可以看出，不论是运行时间还是识别正确率，本文提出的MoG-BBN模型都要优于HMM模型从而验证了该方法的有效性。HMM结果如图10、11、12、13和14所示。

图10 转速200 r/min条件下HMM识别结果

表2 MoG-BBN和HMM对实验数据分析结果对比

图11 转速400 r/min条件下HMM识别结果

图12 转速600 r/min条件下HMM识别结果

图13 转速800 r/min条件下HMM识别结果

图14 转速1 000 r/min条件下HMM识别结果

6 结 论

本文提出了基于MoG-BBN的齿轮箱齿轮磨损状态识别方法，建立了基于VE-EM的模型推理算法，并对EM算法过程进行了改进，避免了该算法局部收敛问题，针对齿轮磨损特征之间的非线性关系，应用CDA方法对高维特征进行降维。最后，应用模型对齿轮箱预置磨损实验数据进行分析，不同工况条件下测试数据的状态识别结果验证了论文所用方法的有效性，磨损状态识别正确率达到了99%，且识别正确率要明显优于HMM模型，为齿轮箱的健康管理提供了科学依据，也为其它类型设备的退化状态识别提供了借鉴。

参 考 文 献

[1]滕红智,赵建民,贾希胜,等.基于CHMM的齿轮箱状态识别研究[J].振动与冲击,2012,31(5):92-96.

TENG Hong-zhi,ZHAO Jian-min,JIA Xi-sheng,et al.Gearbox state recognition based on continuous hidden Markov model[J].Journal of Vibration and Shock,2012,31(5):92-96.

[2]胡海峰,安茂春,秦国军,等.基于隐半Markov模型的故障诊断和故障预测方法研究[J].兵工学报,2009,30(1):69-75.

HU Hai-feng,AN Mao-chun,QIN Guo-jun,et al.Study on fault diagnosis and prognosis methods based on hidden semi-markov model[J].ACTA ARMAMENTARII,2009,30(1):69-75.

[3]Dong M,He D.A segmental hidden semi-Markov model (HSMM)-based diagnostics and prognostics framework and methodology[J].Mechanical Systems and Signal Processing,2007,21:2248-2266.

[4]Purushotham V,Narayanan S,Prasad S A N.Multi-fault diagnosis of rolling bearing elements using wavelet analysis and hidden Markov model based fault recognition[J].NDT&E International,2005,38:654-664.

[5]Widodo A,Yang B S.Support vector machine in machine condition monitoring and fault diagnosis[J].Mechanical Systems and Signal Processing,2007,21:2560-2574.

[6]Yang B S,Han T,Hwang W W.Fault diagnosis of rotating machinery based on multi-class support vector machines[J].Journal of Mechanical Science and Technology,2005,19(3):846-859.

[7]Huang J C.Remote health monitoring adoption model based on artificial neural networks[J].Expert Systems with Applications,2010,37:307-314.

[8]Santosh T V,Vinod G,Saraf R K,et al.Application of artificial neural networks to nuclear power plant transient diagnosis[J].Reliability Engineering and System Safety,2007,92:1468-1472.

[9]Xu B G.Intelligent fault inference for rotating flexible rotors using Bayesian belief network[J].Expert Systems with Applications,2012,39:816-822.

[10]Dey S,Stori J A.A bayesian network approach to root cause diagnosis of process variations[J].International Journal of Machine Tools and Manufacture,2005,45:75-91.

[11]Verron S,Li J,Tiplica T.Fault detection and isolation of faults in a multivariate process with Bayesian network[J].Journal of Process Control,2010,20:902-911.

[12]Sahin F,Yavuz M C,Arnavut Z,et al.Fault diagnosis for airplane engines using bayesian networks and distributed particle swarm optimization[J].Parallel Computing,2007,33:124-143.

[13]Bishop C M.Pattern recognition and machine learning[M].Springer,2006.

[14]Demartines P,Hérault J.Curvilinear Component Analysis: A self-organizing neural network for nonlinear mapping of data sets[J].IEEE Transactions on Neural Networks,1997,8: 148-154.

[15]Lee A J,Lendasse A,Verleysen M.Curvilinear distance analysis versus Isomap[J].European Symposium on Artificial Neural Networks,2002: 185- 192.

[16]Lei Y G,Zuo M J.Gear crack level identification based on weighted K nearest neighbor classification algorithm[J].Mechanical Systems and Signal Processing,2009,23:1535-1547.

[17]Samuel P D,Pines D J.A review of vibration-based techniques for helicopter transmission diagnostics[J].Journal of Sound and vibration,2005,282:475-508.

[18]Dempsey P J,Zakrajsek J J.Minimizing load effects on NA4 gear vibration diagnostic parameter[R].NASA/TM-2001-210671.