王 恒， 季 云， 朱龍彪， 劉 肖

(南通大學(xué)機械工程學(xué)院 南通，226019)

引 言

為了保證機械設(shè)備長期安全運行，工業(yè)現(xiàn)場需要對機械設(shè)備進行狀態(tài)監(jiān)測和性能評估，實現(xiàn)預(yù)知維護和管理。隱馬爾可夫模型為一雙重隨機過程：一個是Markov鏈，用來描述隱藏狀態(tài)間的轉(zhuǎn)移概率，此為基本隨機過程；另一個是描述每個隱藏狀態(tài)下產(chǎn)生觀測值的一般隨機過程。HMM這一雙重隨機特性，能夠很好地描述機械設(shè)備運行過程中觀測到的退化征兆信號(如振動、轉(zhuǎn)速和位移等)與隱藏的衰退狀態(tài)之間的隨機關(guān)系，在機械設(shè)備性能退化評估與預(yù)測中得到廣泛應(yīng)用[1-2]。

傳統(tǒng)HMM的定義和估計過程存在嚴(yán)重不足，例如，HMM參數(shù)訓(xùn)練學(xué)習(xí)中需要預(yù)先確定模型的數(shù)學(xué)結(jié)構(gòu)，不能解決在模型訓(xùn)練過程中出現(xiàn)的過適應(yīng)或欠適應(yīng)問題。劉新民等[3]建立了基于支持向量機(support vector machine, 簡稱SVM)與HMM串聯(lián)結(jié)構(gòu)的故障診斷模型，實現(xiàn)了對直升機減速箱的狀態(tài)識別和診斷。張繼軍等[4]利用HMM自身的狀態(tài)識別和轉(zhuǎn)移回溯能力，結(jié)合多智能體遺傳算法(multi-agent genetic algorithm, 簡稱MAGA)實現(xiàn)了溫控放大器的狀態(tài)分類。上述文獻中HMM的初始狀態(tài)數(shù)都是依據(jù)傳統(tǒng)經(jīng)驗設(shè)定的，缺乏科學(xué)性和通用性。曾慶虎[5]提出了基于最小描述長度準(zhǔn)則(minimum description length, 簡稱MDL)學(xué)習(xí)算法，通過MDL優(yōu)化調(diào)整狀態(tài)數(shù)，但模型狀態(tài)數(shù)還是需要預(yù)先確定，計算過程較為復(fù)雜。騰紅智等[6]基于CHMM對齒輪箱全壽命過程的退化狀態(tài)識別進行了研究，提出了基于K均值算法和交叉驗證相結(jié)合的狀態(tài)數(shù)優(yōu)化方法，但K均值聚類算法仍需要預(yù)先確定狀態(tài)數(shù)，且交叉驗證對不同的退化狀態(tài)數(shù)都要進行訓(xùn)練并檢驗分類錯誤率，計算時間長，效率低。張星輝等[7]建立了一組聚類方法評價指標(biāo)，利用K均值聚類算法對狀態(tài)特征進行聚類，通過指標(biāo)評定結(jié)果從中選取模型的最優(yōu)狀態(tài)數(shù)。但是設(shè)備運行過程中，隨著監(jiān)測數(shù)據(jù)的更新，退化狀態(tài)數(shù)也需要隨之不斷更新，關(guān)于如何有效地確定最優(yōu)退化狀態(tài)數(shù)還需要進一步深入研究。

針對HMM參數(shù)的優(yōu)化問題，提出了一種基于分層狄利克雷過程(hierarchical Dirichlet process, 簡稱HDP)和HMM相結(jié)合的機械零件性能退化評估方法，將HDP引入HMM中，結(jié)合HMM良好的分析和建模能力，實現(xiàn)運行過程中的退化狀態(tài)識別。軸承退化性能評估結(jié)果表明了該方法的有效性和適用性。

1 HDP模型

1.1 DP定義

狄利克雷過程定義為關(guān)于一組分布或者隨機測度的分布，假設(shè)參數(shù)服從一類樣本空間上的寬先驗分布，參數(shù)的后驗分布通過采樣推斷，該狄利克雷模型及其擴展模型則具有良好的聚類特性。近幾年，狄利克雷模型已經(jīng)在機器學(xué)習(xí)、生物信息學(xué)、文本聚類和圖像分割等方面有較好的應(yīng)用。Ferguson首次提出狄利克雷過程的定義，G0為測度空間Θ上的隨機概率測度，參數(shù)α為正實數(shù)。對于測度空間Θ的任意有限劃分A1,,Ar，如果存在如下關(guān)系

(G(A1),，G(Ar))～Dir(αG0(A1)，,αG0(Ar))

(1)

則G服從由基分布G0和參數(shù)α組成的DP，即

G～DP(α,G0)

(2)

1.2 HDP過程

狄利克雷過程可以實現(xiàn)一組數(shù)據(jù)的聚類，但在研究多組數(shù)據(jù)的聚類問題時，單純利用狄利克雷混合模型無法進行建模分析，這大大限制了其應(yīng)用。針對這個問題，引入分層狄利克雷過程[8-9]。分層狄利克雷過程將基礎(chǔ)分布G0擴展為一個服從狄利克雷過程的隨機概率測度G，即G～DP(γ,H)，多個數(shù)據(jù)源Gj～DP(α,G)將共享基礎(chǔ)分布G的離散原子，實現(xiàn)多組關(guān)聯(lián)數(shù)據(jù)的聚類。

HDP無法實現(xiàn)對其過程的采樣，在實際應(yīng)用中往往采用不同形式的構(gòu)造實現(xiàn)HDP過程的應(yīng)用。HDP的Stick-breaking構(gòu)造可以分為兩層進行[10]。第1層Dirichlet過程為

(3)

其中：H為基礎(chǔ)分布；β為第1層分布的權(quán)重；k為任意初始聚類數(shù)；φk為抽樣的原子。

第2層每組Gj的Stick-breaking構(gòu)造方式為

(4)

其中：Gj為擴展分布；aj為第2層分布中每組的權(quán)重。

在HDP模型中，狄利克雷過程作為參數(shù)的先驗分布存在，假設(shè)模型中的觀測數(shù)據(jù)xji表示第j組第i個觀測數(shù)據(jù)，其分布形式定義如式(5)，兩層狄利克雷過程分別用上述的Stick-breaking構(gòu)造[10]。

(5)

HDP模型參數(shù)的含義如下：xji為觀測數(shù)據(jù)；k為類別標(biāo)簽的初始聚類數(shù)；γ，α為HDP模型的超參數(shù)；β，aj為HDP模型每層的權(quán)重系數(shù)；λ為基礎(chǔ)分布的采樣參數(shù)；φk為隨機抽樣的離散原子。

2 基于HDP-CHMM的設(shè)備性能退化評估算法

2.1 基于HDP的退化隱狀態(tài)數(shù)確定

隱狀態(tài)數(shù)K的確定是HMM模型訓(xùn)練和測試的關(guān)鍵，但目前HMM隱狀態(tài)數(shù)大多根據(jù)經(jīng)驗人為設(shè)定，或者通過訓(xùn)練一個最可能的隱狀態(tài)數(shù)來達到近似的目的，這些方法很難將各種類別都考慮到，且新的數(shù)據(jù)中也可能有未知類型出現(xiàn)，依靠訓(xùn)練樣本得到的固定模型結(jié)構(gòu)對新的觀測數(shù)據(jù)的適用性和涵蓋性不強。HDP算法不依賴于訓(xùn)練樣本，而且隨著數(shù)據(jù)的變化，模型結(jié)構(gòu)能夠?qū)崿F(xiàn)自適應(yīng)調(diào)整，實現(xiàn)動態(tài)聚類?；贖DP的隱狀態(tài)數(shù)確定算法流程圖如圖1所示，主要步驟如下。

圖1 基于HDP模型的隱狀態(tài)數(shù)確定Fig.1 The number of hidden states of HDP

1) 初始化HDP模型的參數(shù)，任意給定隱狀態(tài)數(shù)K，迭代次數(shù)M，超參數(shù)α，γ。

2) 假設(shè)輸入的多組觀測數(shù)據(jù)xji服從多項式分布，分布密度函數(shù)為f(·|θ)；觀測數(shù)據(jù)分布的參數(shù)θji服從其共軛分布DP分布，分布密度函數(shù)為h(·)；HDP的權(quán)重系數(shù)β和aj分別通過式(3)和式(4)中的Stick-breaking構(gòu)造先驗分布，觀測數(shù)據(jù)xji的條件分布為

(6)

其中：fkxji(xji)為觀測數(shù)據(jù)x中除了xji外的條件概率分布。

3) 基于馬爾科夫鏈蒙特卡羅(Markov chain Monte Carlo, 簡稱MCMC)統(tǒng)計模擬方法簡化條件概率的計算，并實現(xiàn)模型后驗參數(shù)的更新，利用Augmented representation后驗采樣算法[11]實現(xiàn)對參數(shù)β的更新(ajk類似)，根據(jù)超參數(shù)γ和α對β進行采樣

(7)

其中：m.1,,m.k利用直接分配后驗采樣算法[12]進行采樣。

(8)

其中：mjk為第j組的聚類數(shù)目；Mjk為除了第j組的第k類之外的聚類數(shù)；nj.k為第j組觀測數(shù)據(jù)中屬于k類的數(shù)據(jù)總數(shù)；s(n,m)為Stirling數(shù)；γ為權(quán)重的第K+1個值。

通過對β和m的交叉采樣實現(xiàn)參數(shù)更新，并在Stick-breaking構(gòu)造中運用直接分配后驗采樣算法對觀測數(shù)據(jù)的指示因子Zji進行采樣，每次只更新一個數(shù)據(jù)的聚類屬性，當(dāng)某個類簇中元素個數(shù)為0時，K減1，繼續(xù)迭代，待聚類數(shù)目穩(wěn)定時停止迭代，獲得最終的狀態(tài)數(shù)目K。

4) 超參數(shù)γ和α決定DP過程的離散程度，初始值的選取對聚類結(jié)果影響很大，大大影響了HDP模型的通用性[13]。因此，在進行后驗參數(shù)的更新時，同時對超參數(shù)進行更新。設(shè)α～Γ(a,b)，對任意k=1,2，，n，P(k丨α)服從Beta分布，則

(9)

其中：β′為Beta分布。

從β′(α+1,n)中抽樣參數(shù)η，則參數(shù)α和η的聯(lián)合分布為

p(α,η|k)∝p(α)αk-1(α+n)ηα(1-η)n-1

(10)

參數(shù)α的后驗分布為

p(α|η,k)∝αa+k-2(α+n)e-α(b-log(η))

(11)

則

(α|η,k)～aηG(a+k,b-log(η))+

(1-aη)G(a+k-1,b-log(η))

(12)

其中：aη為權(quán)重。

(α|η，k)服從兩個Gamma分布之和，且分布系數(shù)分別為aη和1-aη，超參數(shù)α和γ分別從后驗分布aη/(1-aη)中不斷迭代采樣實現(xiàn)更新。

2.2 基于CHMM的機械設(shè)備退化狀態(tài)評估

基于CHMM的機械設(shè)備退化狀態(tài)評估算法的流程圖如圖2所示，主要步驟如下。

圖2 基于CHMM的設(shè)備退化狀態(tài)識別Fig.2 Identification of equipment degradation status based on CHMM

1) 根據(jù)步驟1～4確定的狀態(tài)數(shù)K作為CHMM模型的輸入?yún)?shù)，在一定約束條件下，隨機初始化模型的其他參數(shù)(初始狀態(tài)Q，轉(zhuǎn)移矩陣A，觀測值矩陣B，混合高斯數(shù)N)。

2) 運用CHMM模型進行建模，將觀測數(shù)據(jù)作為模型的輸入，通過Baum-Welch算法對觀察數(shù)據(jù)進行訓(xùn)練，利用EM算法求概率參數(shù)模型的最大似然估計，重估權(quán)值w、均值μ和方差ξ，通過不斷迭代估計模型參數(shù)λ=(Q,A,B)。

3) 根據(jù)訓(xùn)練得到的CHMM模型，利用Viterbi算法計算已知模型參數(shù)λ最可能的隱藏狀態(tài)序列，即P(O|λ)的最大值，通過路徑回溯求得每個觀察序列最可能的狀態(tài)，即得到退化狀態(tài)轉(zhuǎn)移曲線。

3 應(yīng) 用

應(yīng)用研究采用美國USFI/UCR的智能維護中心提供的軸承全壽命數(shù)據(jù)[14]，圖3為實驗裝置。4個ZA-2115雙列軸承并列安裝在同一軸上，由恒定轉(zhuǎn)速2 kr/min的直流電機驅(qū)動，在軸承和軸上加載約26 671 N的彈性徑向載荷，采用加速度傳感器采集振動信號，采樣頻率為20 kHz。軸承1在連續(xù)運轉(zhuǎn)約163 h外圈出現(xiàn)嚴(yán)重損傷，共采集982組數(shù)據(jù)，采用軸承1的全壽命數(shù)據(jù)進行建模與分析。

圖3 軸承加速壽命實驗裝置示意圖Fig.3 Schematic diagram of the bearing life device

DP模型的擴展模型HDP利用分層構(gòu)造DP原理，能夠?qū)崿F(xiàn)多組特征值的聚類。筆者基于HDP模型對軸承偏度指標(biāo)(skewness)和峭度指標(biāo)(kurtosis)進行訓(xùn)練，以獲取CHMM模型中的最佳退化狀態(tài)數(shù)。以軸承峭度指標(biāo)為例，分析HDP模型中超參數(shù)α的不同初始值對聚類結(jié)果的影響，由圖4可知，當(dāng)α分別取2，20，2 000時，該模型聚類結(jié)果均能收斂到相同的值，可見，HDP模型的超參數(shù)初始值選取對聚類結(jié)果不敏感。

圖4 不同α值影響分析Fig.4 Effect analysis of different valuesα

選取峭度指標(biāo)和偏度指標(biāo)作為HDP模型的觀測數(shù)據(jù)。假設(shè)觀測數(shù)據(jù)服從Multinomial分布，觀測數(shù)據(jù)分布的參數(shù)服從DP分布。設(shè)定初始聚類數(shù)目K=50，聚集參數(shù)α=20，γ=1，迭代次數(shù)M=300，通過Stick-breaking構(gòu)造HDP，利用不斷采樣更新參數(shù)，實現(xiàn)基于HDP的自動聚類，結(jié)果如圖5所示，聚類結(jié)果趨近于5，且聚類結(jié)果比圖4中單一特征量輸入更穩(wěn)定。因此，將K=5作為CHMM模型隱狀態(tài)數(shù)進行模型的訓(xùn)練和測試。

利用混合高斯模型來擬合各狀態(tài)下的觀測值概率密度函數(shù)，在構(gòu)造連續(xù)隱馬爾科夫模型時，高斯分量數(shù)目N取3，基于Baum-Welch算法重估參數(shù)(Q,A,B)，獲得經(jīng)過多次迭代的重估模型。利用Viterbi算法計算P(O丨λ)的最大值，即得到退化狀態(tài)轉(zhuǎn)移曲線。

圖5 基于HDP模型的多組觀測數(shù)據(jù)混合聚類Fig.5 Multi group observation data mixed clustering of HDP

基于HDP-CHMM獲得的軸承退化狀態(tài)轉(zhuǎn)移曲線如圖6所示。軸承從正常狀態(tài)到失效狀態(tài)的全壽命歷程中一共出現(xiàn)了5次不同狀態(tài)，分別為正常狀態(tài)、早期退化狀態(tài)1、中度退化狀態(tài)2、嚴(yán)重退化狀態(tài)3和失效狀態(tài)。通過該模型可以找出軸承運行時的早期故障點在576號文件處，嚴(yán)重故障點在930號文件處，并能識別軸承在運行過程中的一系列退化狀態(tài)，為軸承的早期維護和維修提供理論指導(dǎo)。

圖6 基于HDP-CHMM的軸承退化狀態(tài)識別Fig.6 Identification of bearing degradation status based on HDP-CHMM

為了驗證HDP-CHMM算法的有效性，與文獻[15]中基于K-S檢驗的軸承退化狀態(tài)識別方法進行了對比。由圖6，7可知，軸承1全壽命數(shù)據(jù)在基于K-S檢驗的退化評估中，早期故障點在第533序列處，嚴(yán)重故障點在962處，與基于HDP-CHMM算法識別的早期故障和嚴(yán)重故障點大體一致。在軸承性能退化階段，由于振動波動異常劇烈，基于距離的K-S檢驗法將軸承振幅跳變劇烈處的狀態(tài)劃分過多，在工作過程中不利于軸承實際狀態(tài)的判定。

圖7 基于K-S檢驗的軸承退化狀態(tài)識別Fig.7 Identification of bearing degradation status based on K-S test

4 結(jié) 論

1) 提出了一種基于HDP-CHMM的機械設(shè)備性能退化評估方法，利用HDP的分層聚類特性，將軸承多個特征指標(biāo)作為模型的輸入，使聚類結(jié)果更加準(zhǔn)確，解決了傳統(tǒng)的CHMM狀態(tài)數(shù)必須預(yù)先設(shè)定的不足，實現(xiàn)了模型結(jié)構(gòu)的自適應(yīng)調(diào)整。

2) 滾動軸承全壽命數(shù)據(jù)建模結(jié)果表明，HDP-CHMM可以有效地找出軸承運行時的早期故障點和嚴(yán)重故障點，便于設(shè)置預(yù)警機制，并能夠識別軸承在運行過程中的一系列退化狀態(tài)，為基于狀態(tài)的設(shè)備退化評估提供了一種新的方法。

3) HDP-CHMM算法和K-S檢驗算法對比可知，HDP-CHMM算法可以對軸承實際運行過程中的狀態(tài)轉(zhuǎn)移進行建模，有效識別出軸承運行中的不同退化狀態(tài)，比K-S檢驗算法基于距離的診斷方法更符合設(shè)備實際退化過程。