基于Dirichlet過程混合模型的滾動軸承運行狀態識別

瞿家明，周易文，王恒，黃希

(南通大學 機械工程學院，江蘇 南通 226019)

如果能夠在軸承運行過程中及時檢測其性能退化的程度、跟蹤早期故障，可以有針對性地組織生產和設備維護，有效防止異常失效的發生。

針對滾動軸承的狀態監測和退化評估的研究集中在神經網絡模型[1]、支持向量機[2]、支持向量數據描述[3]、流形學習[4]等基于數據驅動的方法，但這些方法均無法實時顯示數據對應的退化狀態階段。隱馬爾科夫模型(Hidden Markov Model,HMM)具有狀態隱含、觀測序列可見的雙重隨機屬性，很好地描述了設備運行過程中的衰退狀態與觀測到的征兆信號(如振動、轉速和位移等)之間的隨機關系，在滾動軸承性能退化評估與預測中得到了廣泛應用[5-7]。但HMM在定義和學習過程中需要預設設備所經歷的狀態數，而實際應用中由于缺乏相應的先驗知識，并不能準確地給出模型的退化狀態數，限制了HMM的應用場合[8]。

在滾動軸承運行過程中，狀態需要隨著監測數據的更新而不斷更新，如何有效地利用監測數據實現運行狀態數的自動識別還需要深入研究。因此，提出了一種基于Dirichlet過程混合模型(Dirichlet Process Mixture Model，DPMM)的狀態識別算法，并通過實例驗證其有效性和可行性。

1 DP混合模型

Dirichlet過程(Dirichlet Process，DP)是一種典型的非參數Bayes模型[9]，主要用于非參數問題中的先驗分布，DP可以擬合任意類型的概率分布，且與多項式分布互為共軛分布，因此在觀測值的基礎上，DP后驗分布便于分析與計算。近年來，DP模型已成為機器學習、文本處理和自然語言處理領域中的研究熱點，廣泛應用于各種聚類問題的研究中[10-12]。

1.1 DP的定義

假設G0是測度空間Θ上的隨機概率測度，參數α為正實數[12]。對于測度空間Θ的任意有限劃分A1,A2,…,Ar，若存在如下關系

(G(A1),G(A2),…,G(Ar))～

DDir(αG0(A1),…,αG0(Ar))，

(1)

則G服從由基分布G0和參數α組成的Dirichlet過程，即

G：GDP(α,G0)，

(2)

式中：DDir表示Dirichlet分布。

1.2 DP的構造

在實際應用中往往采用不同形式的構造實現DP的應用。截棍構造(Stick-Breaking Construction，SBC)可以用于獨立構造服從DP的隨機樣本[13]，其設有2個參數：聚集參數α、基礎分布G0。隨機概率質量πk可以通過如下方式構造：對長度為1的棒在比例β1處切割，并將切掉的這部分長度賦值給π1，而后對剩余長度為(1-β1)的棒在其比例β2處切割, 并將切掉的棒的長度賦值給π2，然后按照相同的方式對剩余的棒在比例βk處切割,并將切掉的棒的長度賦值給πk，記為πk：GEM(α)，如圖1所示。

圖1 SBC示意圖Fig.1 Diagram of SBC

隨機概率分布G構造為

(3)

其中，隨機概率質量πk通過以α為參數的β分布產生，隨機原子序列θk從基礎分布G0抽樣。

與SBC類似，中國餐館過程(Chinese Restaurant Process,CRP)構造如下[13]：設中國餐館內可以容納無限多張桌子，所有桌子上都貼上標簽1,2,…,K，進入餐館的顧客可以挑一張桌子坐下。θi被比作第i個進入餐廳的顧客，而不同的φk值表示餐桌。第1個顧客以概率1就座于一張新桌子，第i個顧客以概率mn/(n-1+α)就座于一張已坐下mn個顧客的舊桌子，n為進入餐館的顧客總數；或以概率α/(n-1+α)就座于一張新桌子，即K增加1，而φk:G0,θi=φk，如圖2所示，圓圈代表桌子，方框代表進入餐館的顧客。

圖2 CRP構造示意圖Fig.2 Construction diagram of CRP

從Dirichlet過程的構造可以看出, 無論哪種方式, 均體現了其良好的聚類性質。

1.3 DPMM

DP表現了良好的聚類性質，但其只能將具有相同值的數據聚為一類，如果2組數據不相等，不管兩者間的相似性多強，利用DP均無法實現聚類，因此就需要引入DPMM[13]。在DPMM中，DP作為參數的先驗分布存在，假設觀測數據為xi，其分布服從

(4)

式中：F(θi) 為觀測數據xi服從以θi為參數的分布；當G服從Dirichlet分布時，該模型稱為Dirichlet過程混合模型。

2 基于DPMM的滾動軸承狀態數識別算法研究

基于DPMM的滾動軸承狀態數識別算法流程如圖3所示，其具體步驟如下：

圖3 基于DPMM的滾動軸承狀態識別算法流程Fig.3 Flow chart of algorithm for state recognition of rolling bearing based on DPMM

1) 對軸承振動數據進行特征提取；

2) 隨機初始化DPMM的參數，即初始聚類數目N、迭代次數M和聚集參數α；

3) 構造觀測序列，使其服從Gaussian分布或多項式分布；

4) 觀測序列分布參數服從Gaussian-Wishart分布或DP分布；

5) 通過截棍過程、中國餐館過程構造獲得Dirichlet先驗分布；

6) 通過Gibbs抽樣實現參數后驗更新，當某個類簇中元素個數為0時，N減1，繼續迭代步驟5～6，待聚類數目穩定時停止迭代，獲得軸承最終運行狀態數N。

將軸承觀測數據作為DPMM的輸入，記為o1,o2,…,on，則

(5)

假設觀測序列服從Gaussian分布O～Nd(μ,S)，其參數服從共軛分布Gaussian-Wishart分布μ,R～Wd(m,S)，每個分布都具有相同的超參數和各自的參數。其中，軸承數據任意分組，維數為d，均值為μ，均方差為S，相關系數為R；軸承數據服從分布參數的維數為v，均值為m，均方差為S，相關系數為r。

假設觀測數據服從Gaussian分布，即

(6)

多組數據似然函數為

p(O丨μ,R)=(2π)-nd/2|R|n/2·

(7)

參數服從Gaussian-Wishart分布，即

p(R)=2-vd/2π-d(d-1)/4|S|v/2|R|(v-d-1)/2·

(8)

(9)

(10)

根據Bayes公式，觀測數據及其分布參數的聯合概率可以表示為似然函數與參數先驗的乘積，也可表示為后驗分布與邊緣函數的乘積，即

p(μ,R,O)=p(O|μ,R)·p(μ,R)=

p(μ,R|O)·p(O)。

(11)

觀測數據和分布參數的聯合概率為

(12)

利用Gibbs采樣實現后驗參數的更新

(13)

軸承觀測數據的邊緣概率為

(14)

由(12)，(13)和(14)式推導出后驗分布公式為

(15)

(9)式與(15)式具有相同的分布形式，對比可知后驗分布的更新意味著超參數(r,v,m,S)的更新。在(14)式中，需要計算Γ項的比率以及|S|和|S″|的值。為簡化計算，假設更新前后v不變，將(14)式進行Γ項擴展得

(16)

通過Gibbs采樣實現參數后驗更新，在每個類簇中每次增加(或減少)1個觀測序列，利用Cholesky分解更新S的值，然后根據新的參數值計算每個類簇的似然概率。在迭代過程中，不斷更新每個類簇的均值和方差，若類簇中的值為空，則類簇總數減1，否則將繼續迭代更新。DPMM算法不依賴于訓練樣本，而且隨著觀測數據的變化，模型結構能夠自適應調整，從而實現動態聚類，自動識別軸承的運行狀態。

3 實例研究

為驗證DPMM算法的有效性與可行性，對滾動軸承狀態監測數據[14]進行驗證和分析，選擇負載0 kW，轉速1 797 r/min下SKF6205-2RS型軸承的振動加速度數據。采用電火花加工對軸承進行破壞，故障直徑分別為0.177 8，0.355 6，0.533 4和0.711 2 mm，軸承基本參數及其運行狀態分別見表1和表2。將軸承正常狀態數據和4種人為制造的狀態數據分別取容量相同的樣本相互連接作為觀測數據。

表1 試驗軸承技術參數Tab.1 Technical parameters of test bearing

表2 滾動軸承運行狀態劃分Tab.2 Partition of operation states of rolling bearing

在DPMM算法中，設初始聚類數N=100，聚集參數α=20，迭代次數M=300，取軸承峭度指標作為觀測數據構造觀測序列服從Gaussian分布，參數服從Gaussian-Wishart分布，通過SBC和CRP分別構造DP過程。基于DPMM的滾動軸承運行狀態識別結果如圖4所示，經過300次的迭代，聚類結果趨于穩定且收斂到5，與已知的滾動軸承狀態數(表2)相一致，說明DPMM模型能夠有效識別滾動軸承運行狀態，為軸承退化評估與壽命預測提供了一種新方法。

圖4 基于峭度指標的軸承運行狀態數識別結果Fig.4 Recognition result of number of operation states of bearing based on kurtosis index

聚集參數α不同取值時的聚類結果如圖5所示，從圖中可以看出，無論參數α的初始值如何選取，該模型聚類結果均能收斂到相同的值。進一步分析特征值的選擇對識別結果的影響，取軸承均方根作為觀測數據，其他參數不變，基于DPMM的運行狀態數識別結果如圖6所示。

圖5 不同α值的聚類結果Fig.5 Clustering results of differentα

根據圖4及圖6的結果，對峭度指標和均方根分別獲得滾動軸承運行狀態數進行統計分析，獲得其概率直方圖(圖7)，從圖中可以看出，狀態數聚類為5的概率最大，可以認為最終狀態數為5。

圖6 基于均方根的軸承狀態數識別結果Fig.6 Recognition result of number of states of bearing based on RMS

圖7 狀態數聚類結果概率直方圖Fig.7 Probability histogram of clustering results of number of states

綜合分析可得：構造方法不影響結果的穩定性，參數α不同取值時均能快速實現聚類，結果穩定。DPMM算法對于參數的預設沒有要求，也不依賴滾動軸承特征值、聚集參數等初始參數的選擇，具有很強的適應性和穩定性。

4 結束語

故障診斷算法中退化狀態數的確定一直以來依靠的是設備維護的經驗。這對于算法實現設備退化過程描述和觀測數據的處理帶來了一定的困難。DPMM拋開傳統經驗，利用數據自身的特性進行聚類，為進一步分析設備退化過程以及每個階段的退化研究提供了一種新的方式。

初始參數的選擇是模型訓練的難點，其影響著模型收斂的速度和準確性。DPMM引入了Bayes聚類算法，從另一個角度實現了參數需要人為預先設定的問題，從數據分析的角度也有利于模型的收斂速度和準確性，解決了模型參數預設的問題。另外，DPMM能夠隨著數據變化實現自適應調整，避免了故障診斷過程中數據變化而模型不變進而可能出現的預測錯誤，具有一定的研究價值。