基于可視圖圖譜幅值熵的滾動軸承故障診斷方法

陳 芒， 于德介， 高藝源

(湖南大學 汽車車身先進設計制造國家重點實驗室，長沙 410082)

對機械振動信號進行特征提取是機械故障診斷的關鍵步驟[1]。目前，機械振動信號特征提取技術主要是基于傳統的時域分析、頻域分析、時-頻域分析方法[2]。然而，由于滾動軸承故障振動信號的非平穩性和非線性，傳統的信號處理方法往往難以有效地提取復雜的故障特征[3-4]。因此，有必要進一步研究有效的軸承故障特征提取方法。

近年來，圖信號處理(graph signal processing，GSP)技術隨著圖譜理論的發展而迅速興起，該項技術主要是將傳統的信號處理方法拓展到圖信號的分析處理。目前國內外學者先后提出了許多圖信號處理方法，如圖離散信號處理(discrete signal processing on graph，DSPG)、圖傅里葉變換(graph Fourier transform，GFT)、圖短時傅里葉變換(graph short-time Fourier transform，GSTFT)、圖小波變換(graph wavelet transform，GWT)和圖經驗模態分解(graph empirical mode decomposition，GEMD)等[5-9]。目前，GSP技術已廣泛應用于計算機網絡、氣候變化和圖像處理等領域[10-11]，然而該技術在機械故障診斷領域中的應用并不多見。

GSP方法與傳統的振動信號處理方法不同，前者研究的是振動數據集關系的圖結構，后者研究的是振動數據集本身[12]。復雜網絡是一種特殊的圖結構，對于非線性時間序列，復雜網絡能更好地捕捉其動力學特征。可視圖算法(visibility graph，VG)[13-16]是一種特殊的能將離散時間序列轉變成復雜網絡的算法。對于振動信號，其本身就是典型的離散時間序列，因此可以通過可視圖建網方法將振動信號轉變成可視圖復雜網絡，從而在復雜網絡的頂點域、圖譜域和頂點-圖譜域對故障振動信號分析處理。近年已有一些學者初步研究了基于復雜網絡頂點域的機械故障特征提取方法。孫斌等[17]針對離心泵振動信號的非線性及其非平穩特性，提出了一種基于可視圖網絡節點重要性度量的離心泵振動故障診斷方法。陳安華等[18]提出了基于復雜網絡社團聚類的故障模式識別方法。

本文將復雜網絡與GSP技術引入故障診斷領域，針對滾動軸承振動信號轉變的可視圖信號，在圖譜域對故障振動信號分析處理。熵是一種描述系統不確定性程度的量，能夠反映出系統中量的分布情況[19]。本文根據信息熵理論，定義了圖譜幅值熵(graph spectrum amplitude entropy，GSAE)，并將其作為滾動軸承的單一故障特征，提出了基于可視圖圖譜幅值熵(graph spectrum amplitude entropy of visibility graph，GSAEVG)和馬氏距離(Mahalanobis distance，MD)[20]的滾動軸承故障診斷方法。應用實例表明，本文方法不僅能有效準確地識別滾動軸承不同故障，還能有效區分不同程度內圈故障，且區分效果明顯優于基于傳統的熵指標的故障診斷方法。

1 可視圖

1.1 圖的基本概念

一個無向、連通、加權圖可用G=(V,E,W)來表示，其中：V為頂點的集合；E為邊的集合；W為加權的鄰接矩陣。 若vi和vj這2個頂點之間有邊eij=(vi,vj)連接，那么這條邊的權值就用wij來表示，若沒有邊連接，則wij=0。在實際應用中，wij通常由人為來設定，一般用以下3種方式來定義

W1∶wij=1

(1)

(2)

W3∶wij=‖xi-xj‖2

(3)

式中，xi和xj分別為頂點vi和vj的函數值。式(1)定義的權值都為1，完全忽視了頂點之間的差異；式(2)定義的權值均小于等于1，忽略了部分頂點之間的差異；式(3)定義的權值為兩頂點之間的平方歐氏距離，可以如實反映頂點間的差異。

在鄰接矩陣基礎上，建立能更好反映圖中蘊含在頂點之間關系的圖拉普拉斯矩陣，定義為

L=D-W

(4)

式中，D為度對角矩陣，表示各頂點連接的邊數，其對角元素為di=∑wij。

1.2 可視圖建網算法

可視圖建網算法中，網絡中的每個節點與離散時間序列數據中的每個時間點對應 ，其基本思想如圖1所示。在圖1(a)中，用5個直方條表示一個離散時間序列的前5個數據點，直方條的高度表示每個時間點數值大小；圖1(b)中的每個實點與時間序列中的數據頂點相對應，假如2個直方條的頂端相互可視，則對應的兩點相連。可視性準則如下：

圖1 時間序列可視圖建網Fig.1 Visibility graph of time series

若離散時間序列中任意2個點(va,xa)與(vb,xb)相互可視，那么對任意點(vc,xc)，其中va

(5)

可視圖算法建立的復雜網絡有如下性質：①網絡為無向網絡；②每個點至少與其左右相鄰的節點相連；③橫軸及縱軸經過仿射變換或者坐標比例尺度改變，可視性仍然保持不變。可視圖建網算法能夠保持原有時間序列的一些固有特征，即可以將離散時間序列映射成規則網絡，進而可以捕捉時間序列的幾何結構。

在圖1中，頂點的集合可表示為V={v1,v2,v3,v4,v5}，邊的集合可表示成E={(v1,v2),(v1,v3),(v1,v4),(v1,v5),(v2,v3),(v2,v4),(v3,v4),(v4,v5)}。若采用式(1)定義邊的權值，那么圖的鄰接矩陣為

圖的度對角矩陣為

2 圖譜幅值熵

2.1 圖傅里葉變換

(6)

則GFT的逆變換的定義為

(7)

式中，r為特征值和特征向量的階次，特征值和階次一一對應，而階次與階次之間是線性關系，特征值與特征值之間是非線性關系，在圖譜中通常使用階次。

GFT可以把一個復雜的圖信號分解成一系列不同階次的特征向量的疊加，它建立了圖信號與階次之間的對應關系。因此，圖信號可以通過GFT將其從頂點域變換到圖譜域，進而可以實現類似傳統頻域中“頻率”的概念分析[21]。

2.2 圖譜幅值熵

Shannon[22]首次把熵的概念引入到信息理論中，提出利用信息熵來對事物包含的信息量及反映其狀態變化的信息量進行描述。

(8)

根據信息熵理論，本文定義圖信號f(n)的圖譜幅值熵為

(9)

圖譜幅值熵度量了圖信號的階次分布均勻程度，體現了圖信號能量分布的圖譜域復雜度。當圖信號能量集中在少數幾個階次區域時，圖譜幅值熵Sf取值較小；當圖信號能量在整個階次區域分布比較均勻時，圖譜幅值熵值Sf取值較大。當滾動軸承出現故障，其振動信號必然發生變化，由此轉換得到的可視圖信號必然不同。不同的故障狀態會得到不同的可視圖信號，其包含的故障信息也不同，階次分布的均勻程度亦不同。因此可通過提取可視圖信號的圖譜幅值熵這一特征指標來指示故障是否發生。

3 滾動軸承故障特征提取和分類

3.1 特征提取

首先將采集的軸承振動信號通過可視圖建網方法得到可視圖信號；然后通過式(3)計算可視圖的鄰接矩陣W； 再通過式(4)計算其拉普拉斯矩陣L； 最后通過式(6)和式(9)計算圖譜幅值熵Sf，并把該指標作為滾動軸承的特征參數。

3.2 馬氏距離分類

馬氏距離是由印度統計學家P.C.Mahalanobis提出的，表示數據的協方差距離，它可以有效的計算2個未知樣本集的相似度，進而實現對未知樣本的預測。馬氏距離不受量綱影響，即不受變量之間的相關性干擾。依據實際應用背景，將馬氏距離定義為

(10)

對于滾動軸承故障信號，不同類型故障信號的特征參數之間的馬氏距離較大，同種類型故障信號的特征參數之間的馬氏距離較小。因此，可以將圖譜幅值熵作為滾動軸承振動信號的單一特征參數，選擇馬氏距離作為分類器，從而對滾動軸承故障進行分類。

4 故障診斷流程

本文先將采集的滾動軸承振動信號依據可視圖建網算法轉換成可視圖信號，采用式(3)構造鄰接矩陣。然后將圖譜幅值熵作為滾動軸承特征參數，再通過馬氏距離實現故障分類。詳細診斷流程如下：

步驟1分別采集滾動軸承正常，滾動體故障，內圈故障，外圈故障振動信號，每種狀態進行m次采樣，總共有4m個樣本。將其分成測試樣本和訓練樣本兩組，其中，每種狀態訓練樣本各k個，剩下的4m-4k個樣本作為測試樣本。

步驟2將每個樣本數據轉換成可視圖，求出每個樣本數據圖譜幅值熵Sf，并將該指標作為特征參數。

步驟4根據式(10)，計算每個測試樣本與4種狀態訓練樣本的馬氏距離，測試樣本與4種訓練樣本的MD值分別記作d1，d2，d3，d4。

步驟5比較d1，d2，d3，d4，取其中最小的前k個值所對應的訓練樣本類型為對應的測試樣本故障類型，進而判斷滾動軸承的狀態，達到故障類型識別的目的。

5 應用實例

5.1 試驗一

為了驗證本文的方法的有效性，采用美國凱斯西儲大學電氣工程實驗室的滾動軸承實測數據[21-22]，對滾動軸承不同類型故障進行識別。測試軸承型號為6205-2RS JEM SKF深溝球軸承，采樣頻率fs=12 000 Hz，軸承轉速為1 772 r/min，采樣長度N=2 048。人為設置軸承點蝕故障，故障直徑為0.177 8 mm，深度為0.279 4 mm，通過加速度傳感器分別采集軸承正常，滾動體故障，內圈故障，外圈故障狀態的振動信號，4種狀態下某一樣本振動信號時域波形如圖2所示。

圖2 不同狀態軸承振動信號時域波形Fig.2 Time-domain waveforms of bearings’ vibration signals with different states

每種狀態軸承振動信號各截取20個樣本，對每種狀態取5個樣本作為訓練樣本，其余15個樣本作為測試樣本，總共有20個訓練樣本，60個測試樣本。根據診斷流程步驟計算各個狀態下4個訓練樣本特征參數的均值和方差。測試樣本序列的排列順序依次為正常狀態樣本(1～15)，滾動體故障樣本(16～30)，內圈故障樣本(31～45)，外圈故障樣本(46～60)。按照式(10)計算出每一個測試樣本與4種狀態訓練樣本的馬氏距離，基于可視圖圖譜幅值熵的馬氏距離判別結果如圖3所示。

圖3中：d1為60個測試樣本與正常狀態訓練樣本的馬氏距離；d2為60個測試樣本與滾動體故障訓練樣本的馬氏距離；d3為60個測試樣本與內圈故障訓練樣本的馬氏距離；d4為60個測試樣本與外圈故障訓練樣本的馬氏距離。

從圖3可以看出，4種測試樣本判別結果明顯，取馬氏距離最小值所對應的訓練樣本的狀態為測試樣本的狀態類別。由圖3(a)可以看出，序列1～15的測試樣本d值最小，可以判定序列1～15的測試樣本為正常狀態樣本；由圖3(b)可以看出，序列16～30的測試樣本d值明顯小于其他樣本序列，可以判定序列16～30的測試樣本為滾動體故障樣本；同理，由圖3(c)可以判定序列31～45的測試樣本為內圈故障樣本；由圖3(d)可以判定序列46～60的測試樣本為外圈故障樣本。這判別結果與預期是完全一致的，由此可見，利用可視圖圖譜幅值熵進行特征提取并用馬氏距離判別法判別故障診斷的方法能有效識別滾動軸承故障類型。

圖3 4種不同故障類型的馬氏距離判別結果Fig.3 Classification results of rolling bearings with four fault types using Mahalanobis distance

上述方法產生的結果都是基于式(3)構造可視圖鄰接矩陣得到的，為了證明該權值構造鄰接矩陣的有效性，分別用式(1)、式(2)構造可視圖鄰接矩陣，采用同樣的分析處理方法對試驗一中的試驗數據進行分析，得到的結果與基于式(3)構造可視圖鄰接矩陣得到的結果進行對比，其結果如表1所示。

從表1可以看出，基于W1構造的鄰接矩陣方法無法準確地識別各類故障，因為其忽視了各頂點間的差異；基于W2構造的鄰接矩陣方法識別率也不高，因為該權值忽視了部分頂點間的差異；而基于平方歐氏距離構造鄰接矩陣的方法能夠對滾動軸承故障進行準確有效地識別，且對各類故障的正確識別率均為100%。這就表明基于平方歐氏距離構建鄰接矩陣的方法優于其他權值構建鄰接矩陣的方法，這是因為該權值能如實反映頂點間的差異，進而能捕捉不同故障類型可視圖結構的差異，從而對軸承的不同故障類型進行有效準確地識別。

表1 不同鄰接矩陣對滾動軸承故障數據識別正確率Tab.1 Fault identification accuracies of rolling bearings using different adjacency matrices %

為了表明本文對滾動軸承故障特征提取方法的優越性，將本文的特征提取方法與基于排列熵(permutation entropy, PE)[23]及多尺度排列熵(multiscale permutation entropy, MPE)[24]的滾動軸承故障特征提取方法對比。在同樣的試驗條件下，取上述相同的正常狀態、滾動體故障、內圈故障和外圈故障4種試驗數據進行分析，將其分別編號為1，2，3，4，其中，每種狀態訓練樣本5組，測試樣本15組。首先分別利用排列熵、多尺度排列熵和可視圖圖譜幅值熵提取滾動軸承信號的故障特征；然后根據診斷流程用馬氏距離判別函數作為分類器進行故障識別。其中，基于PE的特征提取方法，嵌入維數為m=6，時延為t=2；基于MPE的特征提取方法，嵌入維數為m=6，時延為t=2，尺度因子為s=4。分類結果如圖4～圖6所示。

對比圖4、圖5與圖6結果可以看出，基于PE以及基于MPE的滾動軸承故障特征提取方法，均出現一定的錯誤。然而，無論針對上述哪種狀態軸承，基于本文可視圖圖譜幅值熵的特征提取方法均有很好的分類效果，故障識別正確率為100%，表明用GSAEVG對軸承振動信號進行特征提取能夠準確的捕捉故障信息，該方法優于基于排列熵以及多尺度排列熵的特征提取方法。

圖4 基于PE特征提取的分類結果Fig.4 Classification results based on the feature extraction of PE

圖5 基于MPE特征提取的分類結果Fig.5 Classification results based on the feature extraction of MPE

圖6 基于GSAEVG特征提取的分類結果Fig.6 Classification results based on the feature extraction of graph spectrum amplitude entropy

5.2 試驗二

為了進一步驗證該方法的有效性，將其應用于滾動軸承模擬故障試驗臺的實測信號。試驗臺如圖7所示，試驗軸承型號均為SKF 6206-2RS1/C3，采樣頻率為8 192 Hz，采樣長度為2 048，電動機轉速為900 r/min。為模擬軸承故障，使用線切割加工技術分別在深溝球軸承內圈、外圈上切割單點故障，其中內圈故障深度為0.4 mm，外圈設置了2種不同程度故障，深度分別為0.2 mm，0.3 mm。通過加速度傳感器分別采集軸承正常，內圈，不同程度外圈故障狀態的振動信號，4種狀態下某一樣本振動信號時域波形如圖8所示。

圖7 滾動軸承模擬故障試驗臺Fig.7 Rolling bearing simulation failure test bench

圖8 不同故障狀態振動信號時域波形Fig.8 Time-domain waveforms of bearings’ vibration signals with different states

采用上述軸承相同的故障診斷方法。先將訓練樣本進行訓練，計算出訓練樣本特征參數及方差；然后根據式(10)計算出每一個測試樣本與每一類訓練樣本的d值；最后根據d值對測試樣本進行分類。在圖9中，測試樣本序列的排列順序依次為正常狀態(1～15)，內圈故障(16～30)，外圈故障1(31～45)，外圈故障2(46～60)。圖9中：d1為60個樣本與正常狀態訓練樣本的d值；d2為60個測試樣本與內圈故障訓練樣本的d值；d3為60個測試樣本與外圈故障1訓練樣本的d值；d4為60個測試樣本與外圈故障2訓練樣本的d值。同理可知，馬氏距離最小值所對應的訓練樣本的狀態為測試樣本的狀態類別，從圖9可以看出，每類故障都能有效的進行識別。

圖9 不同故障類型的馬氏距離判別結果Fig.9 Classification results of rolling bearings with different fault types using Mahalanobis distance

為了進一步體現本文方法的優越性，在此次試驗基礎上，同樣地將本文的特征提取方法與基于PE及MPE的滾動軸承故障特征提取方法對比。其中，基于PE與MPE方法的參數與試驗一中的參數相同。取上述相同的4種試驗數據(正常、內圈故障、外圈故障1、外圈故障2)進行分析，將其分別編號為1，2，3，4，其中，每種狀態訓練樣本5組，測試樣本15組。首先分別利用PE，MPE和GSAEVG提取滾動軸承信號的故障特征；然后根據診斷流程用馬氏距離判別函數作為分類器進行故障識別。分類結果如表2所示。

表2 基于不同方法特征提取的分類結果Tab.2 Classification results based on feature extraction from different methods

由表2可以看出，基于PE以及基于MPE的滾動軸承故障特征提取方法，均出現一定的錯誤。然而，無論針對上述哪種狀態軸承，基于本文GSAEVG的特征提取方法均有很好的分類效果，故障識別錯誤數為0，進一步表明用圖譜幅值熵對軸承振動信號進行特征提取能夠準確的捕捉故障信息，其特征提取效果優于傳統的基于排列熵以及多尺度排列熵的特征提取方法。

6 結 論

針對軸承故障特征難以提取的問題，本文將復雜網絡與圖信號處理引入機械故障診斷領域，提出了基于可視圖圖譜幅值熵和馬氏距離的滾動軸承故障診斷方法。主要結論如下：

(1) 將軸承振動信號轉換為可視圖信號，通過GFT將可視圖信號從頂點域變換到圖譜域，進而在圖譜域對軸承振動信號進行特征提取，從而為故障診斷的特征提取方法提供了一條不同于傳統的時域、頻域或時頻域分析的新途徑。

(2) 可視圖圖譜幅值熵度量了可視圖信號的階次分布的均勻程度，體現了可視圖信號能量分布的圖域復雜度，不同故障類型可視圖信號能量分布圖域復雜度不同。通過可視圖圖譜幅值熵能有效地提取出軸承的不同故障特征。應用實例分析結果表明，與傳統的排列熵及多尺度排列熵特征提取方法比較，本文方法能夠更有效準確地提取滾動軸承故障特征。

(3) 需要說明的是，本文是在有監督分類基礎上進行滾動軸承故障診斷，接下來筆者將研究低標簽及無標簽下的機械故障診斷，即在已知部分故障樣本標簽或無故障樣本標簽基礎上進行機械故障診斷。