基于CDA和MoG-BBN的齒輪磨損狀態識別研究

張星輝， 康建設， 趙勁松,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.