張星輝， 康建設(shè)， 趙勁松,2， 肖 雷， 曹端超

(1.軍械工程學(xué)院，石家莊 050003；2.軍事交通學(xué)院，天津 300161;3.重慶大學(xué)，重慶 400030)

有效的退化狀態(tài)識(shí)別是設(shè)備健康管理的關(guān)鍵步驟。齒輪及齒輪箱是機(jī)械設(shè)備中一種必不可少的部件，對(duì)其進(jìn)行狀態(tài)監(jiān)測(cè)有著重要意義。國內(nèi)外學(xué)者對(duì)齒輪箱不同故障模式的識(shí)別已經(jīng)進(jìn)行了大量的研究，取得了很多成果。而研究單一故障模式從初始發(fā)生到失效整個(gè)過程的退化規(guī)律對(duì)維修活動(dòng)規(guī)劃有很大意義。齒面磨損是齒輪的常見退化模式，通常會(huì)隨著運(yùn)行時(shí)間的流逝而逐漸嚴(yán)重，直到齒輪失效。齒面磨損一般不便直接測(cè)量，但可通過外部測(cè)量設(shè)備來確定其磨損狀態(tài)。目前，隱馬爾可夫族模型(Hidden Markov Models，HMMs)[1-4]、支持向量機(jī)[5-6]、人工神經(jīng)網(wǎng)絡(luò)等模型[7-8]都可用來描述由設(shè)備外部測(cè)量結(jié)果反映內(nèi)部退化狀態(tài)的過程，并被廣泛應(yīng)用于各種設(shè)備的退化狀態(tài)識(shí)別。貝葉斯網(wǎng)絡(luò)是表示不確定性專家知識(shí)和推理的一種流行方法，非常適合于表達(dá)診斷決策問題?；仡欂惾~斯網(wǎng)絡(luò)在設(shè)備故障診斷[9-12]中的應(yīng)用，多是對(duì)于整個(gè)復(fù)雜設(shè)備系統(tǒng)進(jìn)行的簡(jiǎn)單故障推理，而沒有對(duì)具體部件故障與外部觀測(cè)信息之間的關(guān)系進(jìn)行研究。因此，本文構(gòu)建了混合高斯輸出貝葉斯信念網(wǎng)絡(luò)模型(Mixture of Gaussians Bayesian Belief Network，MoG-BBN)用于磨損狀態(tài)識(shí)別。模型應(yīng)用變量消元(Variable Elimination，VE)算法求取后驗(yàn)概率，對(duì)后驗(yàn)概率需要的參數(shù)應(yīng)用期望最大化(Expectation Maximization，EM)算法進(jìn)行估計(jì)。由于標(biāo)準(zhǔn)的EM算法容易局部收斂，因此本文對(duì)算法進(jìn)行了改進(jìn)。最后，齒輪磨損實(shí)驗(yàn)數(shù)據(jù)用來對(duì)模型進(jìn)行驗(yàn)證。

1 MoG-BBN基本原理

貝葉斯信念網(wǎng)絡(luò)(Bayesian Belief Network，BBN)是一種概率推理網(wǎng)絡(luò)，它以圖形節(jié)點(diǎn)表示隨機(jī)變量，節(jié)點(diǎn)之間的有向箭頭表示隨機(jī)變量之間的因果關(guān)系。BBN是一個(gè)在不確定條件下強(qiáng)有力的知識(shí)表達(dá)和推理方法。BBN中節(jié)點(diǎn)所代表隨機(jī)變量的取值可以是連續(xù)的，也可以是離散的。連續(xù)變量可以服從任意分布，離散變量的取值可以是兩個(gè)或多個(gè)。本文構(gòu)建的MoG-BBN模型，其網(wǎng)絡(luò)結(jié)構(gòu)如圖1所示。

該模型可以用如下參數(shù)描述：

(1)X表示隱狀態(tài)，其狀態(tài)取值分別為1,2,…,K，其代表齒輪的磨損狀態(tài)；

(2)M表示混合數(shù)，其取值范圍為1,2,…,b，其代表齒輪箱的某個(gè)磨損狀態(tài)按照第j個(gè)高斯分布混合產(chǎn)生觀測(cè)值；

圖1 貝葉斯信念網(wǎng)絡(luò)

在圖1所示的貝葉斯網(wǎng)絡(luò)中，X是根節(jié)點(diǎn)(M和Y的父節(jié)點(diǎn))，M是中間節(jié)點(diǎn)(X的子節(jié)點(diǎn)，Y的父節(jié)點(diǎn))，Y是葉節(jié)點(diǎn)(X和M的子節(jié)點(diǎn))。

2 VE-EM算法

2.1 VE算法

根據(jù)貝葉斯信念網(wǎng)鏈規(guī)則得：

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

(1)

根據(jù)條件概率公式得：

(2)

式(2)右側(cè)的分子可表示為：

(3)

式(2)右側(cè)的分母可表示為：

(4)

結(jié)合式(1)、(3)、(4)，式(2)可進(jìn)一步表示為：

(5)

2.2 EM算法

假設(shè)訓(xùn)練數(shù)據(jù)集為{y(1),…,y(n)}，由式(5)可以最終求得數(shù)據(jù)所對(duì)應(yīng)的故障模式，因此必須求得P(X)，P(M|X)和P(Y|X,M)，這三個(gè)量分別對(duì)應(yīng)四個(gè)參數(shù)：

(1)π：隱狀態(tài)(設(shè)備退化狀態(tài))概率分布， πi=P(X=i)，i∈{1,2,…,K}；

(2)C：隱狀態(tài)混合系數(shù)，Cij=P(M=j|X=i)，i∈{1,2,…,a}，j∈{1,2,…,b}；

(3)μ：隱狀態(tài)產(chǎn)生高斯分布的均值；

(4)∑：隱狀態(tài)產(chǎn)生高斯分布的方差。

這四個(gè)參數(shù)可以表示為：λ=(π,C,μ,∑)。

EM算法就是求使得P(Y|λ)最大化的參數(shù)λ，然后再通過VE算法求取P(X|Y)。詳細(xì)的推導(dǎo)如下：

(6)

(7)

其中，q(X)為引入的未知分布，并給出它的約束條件。根據(jù)條件概率公式和貝葉斯公式，式(7)可表示為[13]：

(8)

其中的兩項(xiàng)分別為：

(9)

(10)

根據(jù)Jensen不等式，式(10)可表示為：

(11)

由此可得：lnP(Y|λ)≥L(q,λ)。并且當(dāng)且僅當(dāng)q(X)=P(X|Y,λ)時(shí)，lnP(Y|λ)=L(q,λ)。其它情況下，L(q,λ)< lnP(Y|λ)。因此，最大化lnP(Y|λ)等同于最大化L(q,λ)。而L(q,λ)只有在q(X)=P(X|Y,λ)時(shí)取得最大化，將q(X)=P(X|Y,λt-1)代入式(9)得：

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

(12)

其中：

(13)

(14)

E步：

對(duì)所有的(i,j,l)

(15)

根據(jù)貝葉斯公式，式(15)可表示為：

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

(16)

M步：更新參數(shù)

(17)

(18)

(19)

(20)

重復(fù)E步和M步，直到迭代一定的步數(shù)或者相鄰兩次迭代滿足|lnP(Y|λt)-lnP(Y|λt-1)|<ξ，ξ取值一般在區(qū)間[10-3,10-5]。待獲得估計(jì)后的參數(shù)，就可以通過變量消元算法求取各種狀態(tài)產(chǎn)生觀測(cè)值的概率。模型的迭代次數(shù)、ξ可以先取一個(gè)較大范圍的值，如果模型在很小的范圍內(nèi)收斂，那么，在下次訓(xùn)練時(shí)，再進(jìn)一步縮小范圍。反之，應(yīng)加大取值的范圍。

2.3 改進(jìn)后EM算法

由于2.2節(jié)的EM算法很容易收斂于局部最優(yōu)值，而不是全局最優(yōu)值，因此需要對(duì)EM算法進(jìn)行改進(jìn)。該改進(jìn)算法的主要思想是：在每次迭代得到參數(shù) 后，利用該參數(shù)隨機(jī)生成一批數(shù)據(jù)后，對(duì)參數(shù)再次進(jìn)行估計(jì)。這樣，相當(dāng)于每次迭代都進(jìn)行了一次參數(shù)尋優(yōu)，避免了參數(shù)很快的收斂于局部最優(yōu)值。其具體過程可以用以下幾步來描述：

第1步：(模型參數(shù)初始化)以均勻分布隨機(jī)產(chǎn)生隱狀態(tài)概率分布和隱狀態(tài)混合系數(shù)，并分別對(duì)它們進(jìn)行歸一化。以一組訓(xùn)練數(shù)據(jù)(觀測(cè)值)的均值和方差作為隱狀態(tài)產(chǎn)生高斯分布的均值和方差。

第2步：用式(21)對(duì)每一步迭代得到的參數(shù)λi進(jìn)行評(píng)價(jià):

(21)

其中，P(Yl|λi)表示模型參數(shù)λi產(chǎn)生第l組訓(xùn)練數(shù)據(jù)的概率。

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

3 曲線距離分析

為了更好地識(shí)別齒輪的磨損狀態(tài)，往往會(huì)從多個(gè)傳感器采集的振動(dòng)信號(hào)中提取多種特征，特征的維數(shù)會(huì)多達(dá)幾十甚至上百，這些特征間存在一定的相關(guān)性，若把這些特征直接輸入識(shí)別模型將會(huì)導(dǎo)致計(jì)算量急劇增大且計(jì)算結(jié)果不精確。因此，需要對(duì)特征向量進(jìn)行降維，從而提取出少數(shù)關(guān)鍵的特征。

主成分分析和獨(dú)立成分分析是常用的降維方法，它們都假設(shè)特征之間的關(guān)系是線性的，是線性變換方法。而對(duì)于反映齒輪磨損狀態(tài)的特征向量而言，特征之間往往呈現(xiàn)高度的非線性關(guān)系，因此，應(yīng)用非線性變換方法對(duì)其降維更為合理。

CDA[14-15]是用于非線性特征的降維方法，被廣泛的應(yīng)用于圖像和語音識(shí)別中。該算法通過最小化式(22)來達(dá)到降維的目的：

(22)

通常，采用文獻(xiàn)[15]提出的遞減函數(shù)，表示為：

(23)

ECDA可以通過隨機(jī)梯度下降法求得：

i≠j

(24)

其中：α(t)為t的遞減函數(shù)。CDA的詳細(xì)計(jì)算過程和步驟可見文獻(xiàn)[15]。

4 基于CDA和MoG-BBN的齒輪磨損識(shí)別框架

圖2 齒輪磨損狀態(tài)識(shí)別框架

5 實(shí)驗(yàn)及數(shù)據(jù)分析

5.1 實(shí)驗(yàn)設(shè)置

實(shí)驗(yàn)齒輪箱組成如圖3所示，實(shí)驗(yàn)所用齒輪箱型號(hào)為ZD10；動(dòng)力源為電磁調(diào)速電機(jī)，型號(hào)為YCT180-4A；水冷磁粉制動(dòng)器為齒輪箱提供載荷，型號(hào)為FZJ-5。該實(shí)驗(yàn)是對(duì)輸出軸81齒齒輪預(yù)置單齒磨損故障，從齒頂量取的磨損量分別為0.2 mm、0.3 mm、0.46 mm和0.6 mm。運(yùn)行工況分別取200 r/min、400 r/min、600 r/min、800 r/min、1 000 r/min及升速、降速條件下恒載和變載運(yùn)行。采樣頻率為20 kHz，采樣時(shí)間為6 s，每種工況采集20組數(shù)據(jù)。

圖3 齒輪箱結(jié)構(gòu)及傳感器位置示意圖

圖中，①、②、③、④、⑥、⑦號(hào)傳感器位于軸承座上方，⑤和⑧號(hào)傳感器位于軸承座下方。定義齒輪的四種磨損狀態(tài)為：磨損狀態(tài)1(0.2 mm)、磨損狀態(tài)2(0.3 mm)、磨損狀態(tài)3(0.46 mm)和磨損狀態(tài)4(0.6 mm)。選取轉(zhuǎn)速200 r/min、400 r/min、600 r/min、800 r/min、1 000 r/min，負(fù)載0.6 A(磁粉制動(dòng)器控制電流)工況下的數(shù)據(jù)作為分析對(duì)象。

5.2 磨損狀態(tài)識(shí)別

5.2.1 特征提取

(1)時(shí)域特征

該實(shí)驗(yàn)提取的時(shí)域特征有：均值、標(biāo)準(zhǔn)差、均方根值、峰值、偏斜度、峭度。

(2)齒輪統(tǒng)計(jì)特征參數(shù)

統(tǒng)計(jì)特征參數(shù)是NASA技術(shù)報(bào)告中提出的專門檢測(cè)齒輪故障的特征[16-18]。這些特征分別定義如下：

(25)

式中：PPx表示時(shí)域信號(hào)x峰峰值的最大值，Pn表示齒輪嚙合頻率n階諧波的幅值，H表示最大的諧波階數(shù)。

(26)

式中：FM4的實(shí)質(zhì)就是取差分信號(hào)的峭度值，d表示差分信號(hào)(如圖4所示)。

(27)

式中：M′表示前M′個(gè)時(shí)間點(diǎn)測(cè)取的是正常齒輪的信號(hào)，下同。

(28)

式中：M6A是對(duì)差分信號(hào)的6階統(tǒng)計(jì)量。

(29)

(30)

式中：M6A是對(duì)差分信號(hào)的8階統(tǒng)計(jì)量。

(31)

(32)

式中：r表示剩余信號(hào)(如圖4所示)。

(33)

(34)

式中：s表示包絡(luò)信號(hào)。

(35)

(36)

(37)

圖4 特征提取流程

(3)小波包頻帶能量特征

應(yīng)用matlab小波工具箱對(duì)信號(hào)進(jìn)行小波包3層分解，將第3層小波包分解系數(shù)的能量作為齒輪磨損特征，共8個(gè)，分別表示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。

綜上所述，每個(gè)通道可提取27個(gè)特征，8個(gè)通道共216個(gè)特征。

5.2.2 特征降維及歸一化

經(jīng)CDA降維后的特征維數(shù)為30，每種工況條件下每個(gè)磨損狀態(tài)選取10組數(shù)據(jù)用來訓(xùn)練，10組數(shù)據(jù)用來測(cè)試。因此，每種工況條件下的訓(xùn)練和測(cè)試數(shù)據(jù)都為40組(磨損狀態(tài)1、磨損狀態(tài)2、磨損狀態(tài)3和磨損狀態(tài)4各10組)。

對(duì)訓(xùn)練集和測(cè)試集進(jìn)行歸一化預(yù)處理，采用的歸一化映射為：

(38)

式中，x,y∈Rn，xmin=min(x),xmax=max(x)。歸一化的效果是原始數(shù)據(jù)被規(guī)整到[0,1]范圍內(nèi)，即yi∈[0,1],i=1,2,…,n。該歸一化方式稱為[0,1]區(qū)間歸一化。

5.2.3 狀態(tài)識(shí)別結(jié)果及分析

應(yīng)用歸一化后的訓(xùn)練和測(cè)試數(shù)據(jù)按照?qǐng)D2所示的識(shí)別流程對(duì)模型進(jìn)行驗(yàn)證。模型最大迭代次數(shù)為100，ξ取值為1×10-5，模型混合數(shù)為2。轉(zhuǎn)速200 r/min條件下，10組磨損狀態(tài)1的測(cè)試數(shù)據(jù)識(shí)別結(jié)果如表1所示。

表1 轉(zhuǎn)速200 r/min條件下10組磨損狀態(tài)1測(cè)試數(shù)據(jù)識(shí)別結(jié)果

從表中可以看出，測(cè)試數(shù)據(jù)處于磨損狀態(tài)1的概率最大，識(shí)別結(jié)果與測(cè)試數(shù)據(jù)的實(shí)際狀態(tài)相一致。驗(yàn)證了該方法的有效性。在訓(xùn)練過程中，模型迭代15步后即找到了最優(yōu)參數(shù)，也就是ξ的值已經(jīng)小于設(shè)定的1×10-5。對(duì)于模型混合數(shù)的選取，目前還沒有跟好的方法，下一步將對(duì)模型混合數(shù)進(jìn)行優(yōu)化研究。

由于篇幅關(guān)系，所有工況的識(shí)別結(jié)果都用圖來表示，如圖5、6、7、8和9所示。圖中，數(shù)據(jù)的排列為磨損狀態(tài)1(10組)、磨損狀態(tài)2(10組)、磨損狀態(tài)3(10組)和磨損狀態(tài)4(10組)。概率最大的磨損狀態(tài)即為識(shí)別結(jié)果。從圖中可以看出，除轉(zhuǎn)速800 r/min和1 000 r/min條件下各有一組數(shù)據(jù)識(shí)別錯(cuò)誤外，其它測(cè)試數(shù)據(jù)的識(shí)別結(jié)果與實(shí)際結(jié)果都一致?？偟淖R(shí)別正確率可以達(dá)到99%。在MoG-BBN模型的實(shí)際應(yīng)用中，模型參數(shù)對(duì)識(shí)別結(jié)果有很大的影響。如果由訓(xùn)練數(shù)據(jù)估計(jì)得到的模型參數(shù)與代表數(shù)據(jù)的真實(shí)參數(shù)相差太遠(yuǎn)，也就是說，通過訓(xùn)練沒有得到最優(yōu)的參數(shù)。那么，模型的識(shí)別效果會(huì)很差。反之，模型識(shí)別效果會(huì)很好。對(duì)實(shí)驗(yàn)數(shù)據(jù)的分析表明，經(jīng)過改進(jìn)參數(shù)求解算法后的MoG-BBN識(shí)別效果非常好。

圖5 轉(zhuǎn)速200 r/min條件下MoG-BBN識(shí)別結(jié)果

圖6 轉(zhuǎn)速400 r/min條件下MoG-BBN識(shí)別結(jié)果

圖7 轉(zhuǎn)速600 r/min條件下MoG-BBN識(shí)別結(jié)果

圖8 轉(zhuǎn)速800 r/min條件下MoG-BBN識(shí)別結(jié)果

圖9 轉(zhuǎn)速1 000 r/min條件下MoG-BBN識(shí)別結(jié)果

為了更進(jìn)一步驗(yàn)證模型的識(shí)別效果，應(yīng)用HMM對(duì)實(shí)驗(yàn)數(shù)據(jù)進(jìn)行處理并進(jìn)行對(duì)比。對(duì)比結(jié)果如表2所示。從表2可以看出，不論是運(yùn)行時(shí)間還是識(shí)別正確率，本文提出的MoG-BBN模型都要優(yōu)于HMM模型從而驗(yàn)證了該方法的有效性。HMM結(jié)果如圖10、11、12、13和14所示。

圖10 轉(zhuǎn)速200 r/min條件下HMM識(shí)別結(jié)果

表2 MoG-BBN和HMM對(duì)實(shí)驗(yàn)數(shù)據(jù)分析結(jié)果對(duì)比

圖11 轉(zhuǎn)速400 r/min條件下HMM識(shí)別結(jié)果

圖12 轉(zhuǎn)速600 r/min條件下HMM識(shí)別結(jié)果

圖13 轉(zhuǎn)速800 r/min條件下HMM識(shí)別結(jié)果

圖14 轉(zhuǎn)速1 000 r/min條件下HMM識(shí)別結(jié)果

6 結(jié) 論

本文提出了基于MoG-BBN的齒輪箱齒輪磨損狀態(tài)識(shí)別方法，建立了基于VE-EM的模型推理算法，并對(duì)EM算法過程進(jìn)行了改進(jìn)，避免了該算法局部收斂問題，針對(duì)齒輪磨損特征之間的非線性關(guān)系，應(yīng)用CDA方法對(duì)高維特征進(jìn)行降維。最后，應(yīng)用模型對(duì)齒輪箱預(yù)置磨損實(shí)驗(yàn)數(shù)據(jù)進(jìn)行分析，不同工況條件下測(cè)試數(shù)據(jù)的狀態(tài)識(shí)別結(jié)果驗(yàn)證了論文所用方法的有效性，磨損狀態(tài)識(shí)別正確率達(dá)到了99%，且識(shí)別正確率要明顯優(yōu)于HMM模型，為齒輪箱的健康管理提供了科學(xué)依據(jù)，也為其它類型設(shè)備的退化狀態(tài)識(shí)別提供了借鑒。

參 考 文 獻(xiàn)

[1]滕紅智,趙建民,賈希勝,等.基于CHMM的齒輪箱狀態(tài)識(shí)別研究[J].振動(dòng)與沖擊,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模型的故障診斷和故障預(yù)測(cè)方法研究[J].兵工學(xué)報(bào),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.