叢 華，崔 超，劉遠(yuǎn)宏，馮輔周

(裝甲兵工程學(xué)院，北京 100072)

0 引言

齒輪箱作為重要的傳動部件，被廣泛應(yīng)用于各類機(jī)械設(shè)備中，對其進(jìn)行故障診斷與分類研究具有重要的實際意義［1］。振動信號是齒輪箱故障特征信息的有效載體，對振動信號的分析可實現(xiàn)不停機(jī)操作下的齒輪箱故障診斷。但振動信號特別容易受到噪聲干擾，特別是對復(fù)雜設(shè)備，眾多運動部件同時產(chǎn)生振動激勵，使得實際獲取信號的信噪比不高［2］，而目前信號主要的時頻分析方法有短時傅里葉變換、Wigner-Ville 分布、小波尺度分析等，這些方法的時頻分辨率都受到不確定性原理的限制，因此無法同時在時域和頻域得到很好的分辨率。故障識別模型以神經(jīng)網(wǎng)絡(luò)模型應(yīng)用最廣，但運用中需要大量的訓(xùn)練樣本，在訓(xùn)練樣本有限的條件下，泛化能力較差;隱馬爾可夫模型(HMM)具有很強(qiáng)的特征分類能力，最早應(yīng)用于語音識別領(lǐng)域［3］，它所具有的兩層隨機(jī)結(jié)構(gòu)與基于振動的故障診斷機(jī)理較為接近，但目前的研究多用于離散觀察序列，對于連續(xù)信號，需要首先進(jìn)行離散化處理。

排列熵反映了一維時間序列復(fù)雜度，用于信號特征提取可以很好地放大系統(tǒng)的弱變信號，同時檢測出復(fù)雜系統(tǒng)的動力學(xué)突變;CHMM 的輸出序列不存在量化處理，能夠比較精確地表示原始信號，有利于提高識別精度。結(jié)合兩種模型的優(yōu)點，提出一種基于排列熵和連續(xù)隱馬爾可夫(CHMM)的齒輪故障特診斷方法，并通過試驗驗證該方法的有效性。

1 CHMM 建模

首先介紹隱馬爾可夫(HMM)模型。HMM 是一種統(tǒng)計分析模型，用來描述一個含有隱含未知參數(shù)的馬爾可夫過程，其研究對象是一個數(shù)據(jù)序列［4］。模型可記為:λ=(N，M，π，A，B)，或簡寫為:λ=(π，A，B)。其中:N 為模型中Markov 鏈的狀態(tài)數(shù);M 為每個狀態(tài)對應(yīng)的可能觀測值數(shù);π為初始狀態(tài)概率分布矢量;A 為狀態(tài)轉(zhuǎn)移概率矩陣;B 為觀測值概率矩陣。

連續(xù)隱馬爾可夫模型(CHMM)是HMM 的一種改進(jìn)算法。所謂CHMM，是指觀測值為一個連續(xù)隨機(jī)變量，任一狀態(tài)對應(yīng)的觀測值的觀測概率由一個觀測概率密度函數(shù)表示，該密度函數(shù)可以通過高斯概率密度函數(shù)模擬產(chǎn)生。在實際應(yīng)用中，常常用幾個高斯概率密度函數(shù)的線性組合模擬觀測序列的產(chǎn)生。如果高斯概率密度函數(shù)足夠多，則混合高斯密度可以逼近任意的概率分布函數(shù)［5］。每個高斯概率密度函數(shù)都有各自的均值和協(xié)方差矩陣，這些參數(shù)可以通過大量的觀測樣本特征統(tǒng)計得到。記觀測概率密度函數(shù)為b j(o)。

式中:o 表示觀測向量D ×T(D 為維數(shù)，T 為觀測序列長度);M 為每個狀態(tài)包含的混合高斯元個數(shù);cjl表示第j 個狀態(tài)第l 個混合高斯元的權(quán)值(即混合系數(shù));G 表示正態(tài)高斯概率密度函數(shù);μjl表示第j 個狀態(tài)第l 個混合高斯元的均值矢量;Ujl表示第j 個狀態(tài)第l 個混合高斯元的協(xié)方差矩陣。一個具有混合密度形式的CHMM 可以用五元組來表示:

機(jī)械系統(tǒng)的振動觀測序列通常都是連續(xù)變化的信號，盡管可以通過標(biāo)量量化將一個特征矢量信號進(jìn)行離散并符號化，但有可能降低信號所包含的信息特征。而CHMM 直接以特征矢量作為觀測序列，可最大限度地保留信號的特征信息，因此采用CHMM 進(jìn)行故障預(yù)測可以獲得更高的精度。

2 排列熵算法

設(shè)一離散時間序列{x(i)，i=1，2，3，…，n}，對其中任意元素進(jìn)行相空間重構(gòu)［6］，得到

式中，m、τ 分別為嵌入維數(shù)和延遲時間。將x(i)中的m 個重構(gòu)分量{x(i)，x(i +τ)，…，x(i +(m-τ)}按照升序重新排列，即

如果重構(gòu)分量中存在相等的值，即

按照j1，j2值的大小來排序，也就是j1＜j2當(dāng)時，有

因此，對于任意一個時間序列重構(gòu)所得的矩陣Y 中每一行Y(j)都可以得到一組符號序列:

其中，l=1，2，…，k，且k≤m!，m 維相空間映射不同的符號序列(j1，j2…，jm)總共有m!，符號序列S(l)是其中的一種排列［7］。計算每一種符號序列出現(xiàn)的概率p1，p2…，pk，此時，時間序列X(i)的k種不同符號序列的排列熵可以按照Shannon 熵的形式定義為:

當(dāng)Pj=1/m!時，Hp(m)就達(dá)到最大值ln(m!)。為了方便，通常用ln(m!)將Hp(m)進(jìn)行歸一化處理，即

3 故障診斷試驗

3.1 試驗裝置

試驗在自主搭建的變速箱系統(tǒng)試驗臺上進(jìn)行。試驗臺主要包括變速箱、轉(zhuǎn)速控制臺、三相異步電動機(jī)、電磁測功儀、扭矩測試儀、功率控制儀、傳感器和采集設(shè)備等。變速箱為某型坦克變速箱，功率輸入由轉(zhuǎn)速為1 200 r/min、額定功率為5 kW的交流異步電動機(jī)驅(qū)動，轉(zhuǎn)速由轉(zhuǎn)速控制臺調(diào)節(jié)，在終端連接一吸收功率為45 kW 的電磁測功儀，起到負(fù)載的作用(圖1)。信號采集通過基于虛擬儀器平臺搭建的數(shù)據(jù)采集系統(tǒng)完成。

試驗中變速箱掛Ⅲ檔，主動齒輪齒數(shù)為14，從動齒輪齒數(shù)為27，預(yù)加載的扭力為100 N·m。對正常、輕微磨損、嚴(yán)重磨損、斷齒4 種狀態(tài)進(jìn)行試驗。在被動齒輪某齒齒面上，利用磨床在齒面上磨掉少許深度，模擬輕微磨損故障，磨損深度1～2 mm，模擬嚴(yán)重磨損故障，在某齒的齒端處，通過銑床銑掉一個齒，模擬斷齒故障(圖2)。

圖1 試驗裝置Fig.1 Experimental facilities

圖2 齒輪故障Fig.2 Failed gear

3.2 數(shù)據(jù)預(yù)處理和特征提取

試驗中傳感器選用DYTRAN 公司生產(chǎn)的3215M1 型單自由度加速度傳感器，測點布置在靠近輸出軸端的箱體表面外側(cè)，該測點距離振源較近，能夠較好地獲取檔位齒輪故障位置附近的振動信號。采用編碼器控制采集，使每次采集的相位都一致，采樣頻率為20 kHz，數(shù)據(jù)采集寬度為10 s，每種狀態(tài)采集5 個數(shù)據(jù)樣本，每個樣本文件含有20 萬個采樣點。圖3 所示為每組第一個樣本的振動情況?？梢钥闯?，齒輪存在磨損和斷齒時信號的振幅明顯高于正常狀態(tài)，并且存在周期性沖擊的信號。根據(jù)參數(shù)計算得到故障齒輪的嚙合頻率為233.5 Hz。

原始信號含有大量的噪聲，致使有用信息淹沒在噪聲里，所以首要進(jìn)行降噪處理。采用小波相關(guān)濾波法進(jìn)行降噪，該算法是利用各層信號之間的相關(guān)特征將信號中的重要特征信息與噪聲區(qū)分開，從噪聲中檢出重要的信號邊緣，并移除噪聲，使得信噪比大大提高［4，8］。

采用前面介紹的排列熵算法對降噪信號進(jìn)行分析，嵌入維數(shù)m=4，延遲時間τ=1(試驗發(fā)現(xiàn)m 取4～6，τ 取1～3，結(jié)果差異無顯著意義)。將樣本數(shù)據(jù)10 萬個點每200 個點為一組，共分為400 個組，計算的排列熵變化情況如圖4 所示。

排列熵反映了一維時間序列復(fù)雜度，在這里運用排列熵可以反映齒輪振動信號的復(fù)雜度。在轉(zhuǎn)速、負(fù)載相同的情況下，齒輪振動信號的排列熵越大，振動信號越隨機(jī);反之，排列熵越小，振動越平穩(wěn)。當(dāng)齒輪發(fā)生磨損或斷齒時，機(jī)械設(shè)備發(fā)生異常必然產(chǎn)生異常頻率成分，由于有缺陷的齒輪在嚙合過程中存在的低頻、低振幅所激發(fā)的高頻、高振幅共振，排列熵值能夠有效反映這些突變信號特征。計算得到各樣本的排列熵參數(shù)如表1所示。

3.3 CHMM 建模

由于CHMM 的初始概率π 和狀態(tài)轉(zhuǎn)移矩陣A 的初始值選取對模型的訓(xùn)練結(jié)果影響不大，這里采用等概率方式產(chǎn)生［9］，即各模型初始狀態(tài)轉(zhuǎn)移矩陣和初始觀測值概率矩陣分別為

圖3 原始振動信號Fig.3 Vibration signals of the gear

圖4 降噪信號排列熵Fig.4 Permutation entropy of the noise reduction signals

表1 樣本排列熵參數(shù)Table 1 Permutation entropy parameter of the examples

續(xù)表1

而CHMM 中觀測值概率矩陣的初值對模型的訓(xùn)練結(jié)果影響較大，并且需要考慮不同樣本對輸出結(jié)果的影響，這就涉及到CHMM訓(xùn)練的問題，即混合高斯概率密度函數(shù)中的均值、方差和權(quán)系數(shù)應(yīng)該如何初始化的問題［10］。采用的方法是:將幾個屬于同一觀察樣本的特征矢量組成一個大的矩陣(試驗中將每個狀態(tài)的前3 個樣本組成一個大的觀察序列)，然后對這個矩陣進(jìn)行分段，對每一段的特征矢量進(jìn)行K-調(diào)和均值聚類，得到連續(xù)混合高斯概率密度函數(shù)。從理論上講，使用混合高斯元個數(shù)越多，擬合模型的精度越高，然而計算速度會變得越慢［11］?？紤]到這一問題，在CHMM 的建模過程中，對每段特征矢量，選取4個高斯元來構(gòu)建連續(xù)混合高斯概率密度函數(shù)，每個狀態(tài)模型的狀態(tài)數(shù)取為4。按照上述方法可對正常狀態(tài)下觀測值概率矩陣的混合高斯概率密度函數(shù)中的均值矩陣M、協(xié)方差矩陣C 和權(quán)系數(shù)矩陣W 進(jìn)行如下形式的初始化(其他狀態(tài)類似):

圖5 給出了4 種狀態(tài)下模型的訓(xùn)練曲線，最大迭代次數(shù)為50，4 種模型的對數(shù)似然概率隨著迭代次數(shù)的增大而增大，在20 步內(nèi)均達(dá)到收斂。

3.4 試驗結(jié)果分析

各種狀態(tài)的模型建立以后，用各組數(shù)據(jù)的后2 個樣本診斷故障，驗證模型的準(zhǔn)確性，模擬試驗的測試統(tǒng)計結(jié)果見表2。表中數(shù)值為最大對數(shù)似然概率值，數(shù)值越大，表示越接近模擬的狀態(tài)，從表2 中可以看出，算法對每種故障的診斷結(jié)果較為理想。在實際診斷中，可以通過多次采集數(shù)據(jù)進(jìn)行訓(xùn)練，進(jìn)一步提高故障診斷的精度。

表2 故障診斷試驗結(jié)果Table 2 Results of fault diagnosis experiment

圖5 CHMM 訓(xùn)練圖Fig.5 Training diagram of CHMM

4 結(jié)論

1)以齒輪振動信號為研究對象，運用小波相關(guān)法進(jìn)行降噪，并利用排列熵算法較好地提取了目標(biāo)齒輪的信號特征。

2)應(yīng)用CHMM 進(jìn)行齒輪的故障診斷，減少了對連續(xù)觀察序列進(jìn)行離散化處理這一步驟，提高了診斷精度，并通過齒輪箱試驗系統(tǒng)驗證該方法的有效性。

［1］丁康，李巍華，朱小勇.齒輪及齒輪箱故障診斷使用技術(shù)［M］.北京:機(jī)械工業(yè)出版社，2005:10-26.

［2］李斌，郭瑜，劉亭偉，等.基于獨立分量分析與包絡(luò)階比分析的齒輪箱多振源特征提?。跩］.振動與沖擊，2012，31(19):68-72.

［3］國立新，莫福源，李昌立.基于連續(xù)高斯混合密度HMM 的漢語全音節(jié)語音識別研究［J］.聲學(xué)學(xué)報，1995，20(5):321-329.

［4］陸汝華，王魯達(dá).基于狀態(tài)加權(quán)合成的HMM 滾動軸承故障診斷［J］.軸承，2011(10):53-56.

［5］王國峰，李玉波，秦旭達(dá)，等.基于TVAR-HMM 的滾動軸承故障診斷［J］.天津大學(xué)學(xué)報，2010，43(2):168-173.

［6］印欣運，何永勇，彭志科，等.小波熵及其在狀態(tài)趨勢分析中的應(yīng)用［J］.振動工程學(xué)報，2004，17(2):165-169.

［7］馮輔周，司愛輝，饒國強(qiáng)，等.基于小波相關(guān)排列熵的軸承早期故障診斷技術(shù)［J］.機(jī)械工程學(xué)報，2012，29(3):73-79.

［8］王文歡，王細(xì)洋，萬在紅.基于細(xì)化譜和隱馬爾可夫模型的齒輪故障分類方法［J］.失效分析與預(yù)防，2014，9(1):30-34.

［9］Atlas L，Ostendorf M，Bernard G D.Hidden Markov models for monitoring machining tool-wear ［C］.IEEE International Conference on Acoustics，Speech and Signal Processing，2000:3887-3890.

［10］Liu B，Riemenschneider S，Xu Y.Gearbox fault diagnosis using empirical mode decomposition and Hilbert spectrum［J］.Mechanical Systems and Signal Processing，2006，20(3):718-734.

［11］Hilton M L，Ogden R T.Data analytic wavelet threshold selection in 2-D signal denoising［J］.IEEE Transaction on Signal Processing，1997，45(2):496-500.