一種新的圖譜域滾動軸承早期故障檢測與識別方法

陳子旭， 朱振杰， 盧國梁

(1. 山東大學 機械工程學院，濟南 250061； 2. 山東大學 高效潔凈機械制造教育部重點實驗室，濟南 250061；3. 山東大學 機械工程國家級實驗教學示范中心，濟南 250061)

滾動軸承是旋轉機械的關鍵部件，其損壞會直接影響設備的運行過程，導致設備停機、甚至引發人員傷亡。對滾動軸承的健康狀態進行實時、在線診斷具有重要意義[1]。

滾動軸承頻發的故障大多為局部缺陷形式的潛在故障，類型較復雜，其機械信號特點有：①強烈的非線性、非平穩特性；②常伴有強噪聲干擾；②調制源弱導致早期故障特征通常很微弱。這些都給軸承故障在早期階段的檢測與識別帶來很大困難[2]。軸承早期故障檢測與識別的關鍵是從振動信號中提取故障特征信息。傳統的時域方法通過提取軸承振動信號的峰度、峭度、均方根值等統計指標，對信號進行動態描述，但其容易受到噪聲干擾，魯棒性不足[3]。以傅里葉變換為基礎的傳統頻域分析方法在處理非線性、非平穩信振動號時，常常無法準確描述信號頻率隨時間的動態特征，制約了其在早期故障檢測中的應用[4]。為了彌補上述方法的不足，時頻分析方法被相繼提出，代表性的時頻分析方法包括短時傅里葉變換、小波變換等。這些方法被認為可提供信號完整的時頻分布，因此已被廣泛應用于軸承故障診斷領域[5]。隨著數字信號處理技術的不斷進步，近些年來，在傳統時頻分析方法基礎上，自適應時頻分析方法得到了較大發展。Huang等[6]提出的經驗模態分解(empirical modal decomposition，EMD)解決了傳統時頻分析存在的非自適應問題，更加適合分析多分量信號，但其也存在端點效應、模態混淆的缺點。Smith[7]提出的局部均值分解(local mean decomposition, LMD)是在EMD的基礎上，對端點效應、模態混淆問題進行了改善，能獲得原始信號更加精確的時頻分布。盡管LMD方法在軸承故障診斷領域獲得了較多應用，但其仍存在兩個問題：

(1) 由LMD分解得到多尺度分量信號數據量龐大、復雜，因此仍需進一步計算以提取更低維度的數據特征。Tian等[8]提出基于乘積函數(product functions，PF)分量構建漢克爾矩陣，利用奇異值分解壓縮數據尺度以獲得更加穩定的特征向量值；王洪明等[9]通過計算全矢能量熵從PF分量中提取有用信息，獲取得到了低維度的數據特征。

(2) 由LMD分解可以獲取信號的強數據特征，而軸承早期故障特征微弱，因此需要對LMD進行增強處理以獲取軸承早期微弱故障特征。最大相關峰度反褶積處理能夠有效消除噪聲對PF分量的影響，用于增強故障特征信息[10]；基于PF分量構建Teager能量譜能夠有效識別滾動軸承故障信號中的沖擊特性，從而能夠提取微弱故障特征[11]；王志堅等[12]提出了一種將掩膜信號法與LMD相結合的LMD-MS(local mean decomposition-mask signal)方法，成功提取了軸承的微弱故障特征。

隨著圖論的不斷發展，圖模型作為一種高效的數學模型，已被應用于機械設備狀態監測領域[13]。Lu等[14]采用圖模型的方式進行數據建模，通過對模型進行距離度量實現對旋轉機械的動態描述。Wang等[15]首先采用無向加權圖的建模方式基于信號的周期譜進行建模，通過計算模型中位圖和距離指標來描述機械設備的動態狀態。Wang等[16]在之前工作的基礎上進一步擴展了圖建模策略，構建了時空圖模型并通過比較試驗驗證了該方法的有效性。受以上研究的啟發，本文提出了一種新的圖譜域滾動軸承早期故障檢測與識別方法。本方法實現了兩方面改進：①建立的異常度動態指標實現了低維度軸承健康狀態動態描述；②基于PF分量構建圖模型進行圖譜分析，實現了對LMD的增強處理。通過上述改進，本文所提方法可有效提取軸承早期故障的微弱特征。

1 方法簡介

方法主要包含以下3個步驟，如圖1所示。

圖1 方法流程圖

步驟1滾動軸承振動信號動態建模。采用LMD對滾動軸承振動信號進行分解，將分解得到的多尺度分量信號作為輸入構建圖模型，計算相鄰圖模型的相似性指標以得到能描述滾動軸承健康狀態的動態指標。

步驟2早期故障檢測。基于上述構建得到的動態指標，基于假設檢驗的方法，采用拉依達準則對軸承早期故障進行決策。

步驟3故障識別。對于檢測到的故障異常，進一步采用模式識別技術進行故障類型的識別。

2 滾動軸承健康狀態動態建模

滾動軸承由于復雜、惡劣的工況導致其振動信號經常受到噪聲等冗余信息的干擾，并且不利于微弱故障特征的提取。LMD可有效將噪聲等冗余信息同有效信息隔離，進行多尺度分析。同時，隨著圖論的發展，圖建模方式已經被證實在提取微弱故障信息方面的重要作用。因此，針對以上方法的優勢，本文基于LMD和圖論方法，對滾動軸承健康狀態進行動態建模。

2.1 局部均值分解

作為一種自適應分解方法，LMD能將信號分解為一系列具有物理意義的分量信號，反映信號能量在多尺度空間下的分布情況。針對監測獲取的滾動軸承振動信號，記為x(t)，其經過LMD分解的步驟如下。

步驟1根據所有的局部極值點ni，求出所有局部極值平均值mi和包絡估計值ai，用滑動平均法處理后得到局部均值函數m11(t)和a11(t)。

步驟3將上述迭代過程的包絡函數相乘得到包絡信號，即a1(t)=a11(t)·a12(t)…a1n(t)。

步驟4將純調頻信號s1n(t)與包絡信號a1(t)相乘，得到第一個PF分量，即FP1(t)=s1n(t)·a1(t)。

步驟5將FP1(t)從x(t)中分離出來，得到u1(t)，將其作為振動信號重復以上計算過程得到u2(t)，直到uq(t)為單調函數時停止迭代。

步驟6滾動軸承振動信號經過上述迭代過程，被分解為一系列PF分量和一個殘余分量。

(1)

本文首先對滾動軸承振動信號采用無重疊滑動窗的方式進行LMD分解，原始信號被分解為一系列PF分量，即x(t)→{FP1(t),FP2(t),…,uq(t)}。從LMD定義來看，殘余分量uq(t)為一單調函數，其包含的振動信號極其微弱，具有的能量接近于零；并且振動信號分解后的前5個分量幾乎包含了滾動軸承振動信號的所有頻率成分，其能量和與信號的總能量近似相等[17]。因此，本文選取前5個PF分量用于圖建模。

2.2 圖建模

圖模型通常由一系列的節點V和邊L組成，即G={V,L}。通過對每對節點vi和vj之間計算距離，比如歐幾里得距離、曼哈頓距離、赫凌格距離等，以表示兩節點所構成的邊lij的權重。圖模型可以有效的表示相鄰節點以及互連節點之間的連接信息，并且應用到信號的多尺度分析中。本文將分解得到的前5個PF分量作為圖模型的節點序列，并采用最常用的歐氏距離度量方式，賦予圖模型每條邊權重進行圖建模，如圖2所示。本文以t1時刻為例，給出了圖建模過程示意圖，其建模過程分為以下幾個步驟。

圖2 圖建模示意圖

步驟1將分量FP1～FP5作為圖模型的節點序列，5個分量在t1時刻的取值分別為{FP1,1,FP2,1,FP3,1,FP4,1,FP5,1}，作為圖模型的節點{v1,v2,v3,v4,v5}，每對節點vi和vj之間的連接作為圖模型的邊lij，其中i,j∈[1,5],Z。

步驟2計算每對節點vi和vj之間的歐式距離dij，作為邊lij的權重對其進行賦值。

di,j=|vi-vj|

(2)

步驟3將計算得到的歐式距離dij如式(3)所示，以鄰接矩陣的形式表示，得到圖模型的矩陣表示。

(3)

使用來自每個PF分量的多尺度信息在每個數據點上構建圖模型。因此，每個滑動窗內部的振動信號則可以表示為一系列的圖模型，即{G1,G2,…,Gl}，l為滑動窗的長度，本文將l設為5 000。構建的圖模型除了可以描述時序特性，還包含信號不同尺度之間的相關信息。

在圖論中，均值圖是用于度量一組圖模型平均水平的常見方式。對于在每個滑動窗內構建得到的圖模型序列，計算其均值圖用以對該段數據狀態進行描述，計算公式如下

(4)

2.3 相似性度量

構建得到的圖模型具有較高的維度，同時，為了更好地描述滾動軸承動態特性、便于故障檢測的決策，引入了相似性度量指標。本文認為LMD分解得到的各分量對于提取故障信息具有相同的貢獻度，因此，采用歐式距離度量方式計算相鄰圖模型之間相似性，得到異常度分數序列，對滾動軸承運行過程中狀態變化的可能性進行量化。計算公式如下所示

(5)

3 早期故障檢測

當滾動軸承正常運轉時，其振動信號分布接近正態分布，此狀態下圖模型結構穩定，異常度分數服從高斯分布；而當滾動軸承在某一時刻發生故障時，隨之而來的是故障特征頻率的出現和沖擊成分的增多，其振動信號將偏離正態分布。基于此，以高斯分布為假設，采用拉依達準則對異常度分數進行假設檢驗(本文采用6σ準則)，實現早期故障時刻的定位。原假設和備擇假設如下所示

(6)

式中，μh-1和σh-1分別為歷史數據{s1,s2,…,sh-1}的平均值和標準差，計算公式如下

(7)

(8)

假設H0設定：當前時刻數據與歷史數據服從同一個正態分布，即滾動軸承處于正常運行狀態。備擇假設H1設定：當前時刻數據出現了異于正常的狀態變化，該時刻即為檢測到的由正常轉向故障的狀態變化點。實線為異常度分數曲線，虛線為上下6σ控制限，圈表示檢測到的狀態變化點，如圖3所示。當滾動軸承運行出現了異于正常的狀態變化時，圖模型結構的變化會導致異常度分數發生較大變化，異常度分數不再服從先前分布，通過判斷當前時刻異常度分數值是否處于此前分布的置信區間，便可以確定當前時刻滾動軸承是否發生狀態變化。

圖3 早期故障檢測示意圖

4 故障識別

基于第3章早期故障檢測，一旦檢測到滾動軸承運行狀態發生異常，可利用模式識別方法(如支持向量機(support vector machines, SVM)、人工神經網絡、隨機森林、深度學習等)實現故障類型的識別。在常見的識別方法中，K最近鄰(k-nearest neighbors, KNN) 算法被廣泛應用。KNN算法通過測量不同特征值之間的距離，依照多數表決的決策規則進行分類，即由待測樣本的K個最近鄰的訓練實例中的多數類別決定待測樣本的類別。本文給出了KNN算法用于分類的示意圖，如圖4所示。首先計算測試數據與各訓練數據之間的距離，按照距離的遞增關系進行排序，如果一個樣本在特征空間中的K個最相似(即特征空間中最近鄰的前K個點)的樣本中的大多數屬于某一個類別，則該樣本也屬于這個類別。KNN是一種非參的、惰性算法模型，其簡單易用，對異常值不敏感，可用于非線性分類，同時具有較低的運算復雜度，無需大量數據進行網絡訓練[18]。

圖4 KNN算法原理圖

本文采用KNN開展故障類型識別的試驗。需要特別指出的是，本文主要關注點在于提出了一種基于圖的振動信號建模與故障表征方法，因此未對其他模式識別方法進行驗證。實際上，本文方法也可與其他模式識別方法進行結合，以進一步提高其準確率。一旦檢測到滾動軸承運行狀態發生異常，則狀態變化點之后的數據為待識別數據，將待測數據的均值圖結合KNN算法對故障類型進行識別。值得提出的是，均值圖表示為鄰接矩陣的形式，在輸入KNN之前需將其轉換為行向量。

5 試驗及討論

5.1 試驗數據

分別在XJTU-SY[19]和凱斯西儲大學(Case Western Reserve University, CWRU)[20]數據集上進行了驗證。其中，XJTU-SY數據集用于驗證早期故障檢測，該數據集是由西安交通大學設計科學研究所和基礎構件研究所以及長興晟陽科技有限公司提供的滾動軸承加速退化試驗數據。該試驗數據在3種工況下取得：①轉速2 100 r/min，徑向力12 kN；②轉速2 250 r/min，徑向力11 kN；③轉速2 400 r/min，徑向力10 kN，每個工況下有5個軸承。數據采樣頻率為25.6 kHz，采樣間隔為1 min，每個采樣持續時間為1.28 s，失效原因包括外圈磨損、內圈磨損、滾動體損壞、保持架斷裂等。每一個測試軸承的詳細信息，如表1所示。有關該試驗系統以及試驗詳情見文獻[19]。

表1 XJTU-SY軸承數據集信息

CWRU數據集用于驗證故障類型識別，CWRU數據集來自凱斯西儲大學軸承數據中心。數據集包含正常數據、采樣率12 kHz的驅動端故障數據、采樣率48 kHz的驅動端故障數據、風扇端故障數據。本文采用48 kHz驅動端故障數據進行試驗，試驗數據在4種工況下獲得：①負載0，轉速1 797 r/min；②負載745.69 W，轉速1 772 r/min；③負載1 491.39 W，轉速1 750 r/min；④負載2 237.07 W，轉速1 730 r/min。每種工況下包含滾動軸承4種工作狀態：①正常(Normal)；②外圈磨損(OF)；③內圈磨損(IF)；④滾動體磨損(RF)，其中，3種故障類型均是通過放電加工產生，每種故障類型包括3種不同的破壞尺寸，如表2所示。數據集共包括10種不同類型的數據。有關該數據集更為詳細的介紹見文獻[20]。

表2 CWRU 驅動端故障數據類型

5.2 故障檢測試驗

XJTU-SY數據集既包含緩慢失效的故障信號，又包含突然失效的故障信號。本文方法在故障檢測階段的兩個例子，如圖5、圖6所示。圖5(a)、圖6(a)為XJTU-SY數據集中數據Bearing1_3和Bearing3_3，分別對應這兩種不同退化形式的信號。圖5(b)、圖6(b)為本文檢測結果，其中圈為檢測得到的狀態變化點。若故障發生后，方法發出預警，認為方法成功檢測到狀態變化，并且預警時刻與真正開始發生故障時刻的延遲越小則檢測效果越好；若方法發出預警，但軸承此時實際并未發生故障，則認為方法此時為誤報。對于突然失效的故障信號，其時域波形發生較大突變，通過肉眼便能確定故障開始時刻[21]；而對于緩慢失效的故障信號，通常對應軸承退化的多個階段，很難確定故障開始時刻[22]。因此針對緩慢失效故障信號，對檢測到的狀態變化點所在滑動窗數據時域波形和頻譜分析，并與正常數據進行比較，以確定方法檢測的準確性。本文以Bearing1_3數據為例進行說明，該軸承在試驗時發生了軸承外圈故障，正常數據和檢測到的故障數據的時域波形和頻譜圖，如圖7、圖8所示。從時域來看，檢測到的故障數據與正常數據相比波形上有較大差距，出現了較多故障發生時典型的調幅和沖擊成分。該數據集軸承一倍轉頻理論值為34.38 Hz，從頻譜圖來看，正常數據中主要包含35.85 Hz，102.40 Hz，291.90 Hz，327.70 Hz頻率成分，分別對應軸承轉動的基頻和倍頻；檢測到的故障數據頻譜中，除了一倍轉頻成分35.85 Hz外，還出現了107.50 Hz及其二倍頻215.10 Hz、三倍頻322.60 Hz等倍頻成分，該頻率成分接近外圈故障特征頻率理論值107.91 Hz，由此可以獲知本文方法成功做出狀態變化預警。

圖5 Bearing1_3檢測結果

圖6 Bearing3_3檢測結果

圖7 正常數據波形及頻譜

圖8 故障數據波形及頻譜

為了進一步驗證方法的優勢，本文將所提方法與廣泛采用的均方根值(root mean square, RMS)、峭度(Kurtosis)、偏度(Skewness)、均方誤差(mean square error, MSE)、標準差(standard deviation, SD)、方差(variation, Var)、特征能量(signature energy, SE)、特征能量比(signature energy ratio, SER)共8種指標進行了比較。將每種方法首次檢測到的異常警報作為軸承早期故障，利用查準率、查全率和綜合測度(F1_score)3個指標對每個方法進行綜合評價，計算公式如下

(9)

(10)

(11)

所有方法的評價結果，如表3所示。其中，本文方法同SE、SER、MSE、Var指標均達到了查全率100%，查準率100%，綜合測度1.0，優于RMS、Kurtosis、SD、Skewness指標。進一步地，本文在表4給出了本文方法同SE、SER、MSE、Var 4種方法預警時刻的對比。以表4中Bearing1_3為例，本文方法在第386個滑動窗處發出預警，早于其他方法預警時刻。從表4可知，除Bearing1_2和Bearing3_2，本文方法發出故障預警的時刻均早于其他方法，說明本文方法能夠更有效提取滾動軸承振動信號中的早期微弱故障信息。

表3 試驗對比結果

表4 不同方法的預警時刻

在實際應用中，計算時間同樣對早期故障檢測有著重要影響[23]。本文以“Bearing1_3”為測試數據，在MATLAB R2017b，i5-5200U CPU 2.2 GHz，RAM 4.00 GB操作環境下，對不同方法計算時間進行測試。測試數據前386個滑動窗內數據每個滑動窗內的平均計算時間，如表5所示。幾種方法前386個滑動窗數據總的計算時間，如圖9所示。相比4種對比方法，本文方法計算時間較高，并且隨著數據量增加升高明顯。但是，值得提出的是，本文方法在第386個滑動窗處已成功發出預警，而其他方法在較晚時刻，甚至是完全失效階段才發出預警。相比其他方法，本文方法對早期故障階段的信息更為敏感，能在故障發生后的較早階段發出故障預警，并且對于緩慢失效過程尤為明顯。因此，本文通過LMD和圖論方法提取信號多尺度信息，相比其他方法能夠更有效提取早期微弱故障特征。在實際應用中，在設備故障的早期階段發出預警可幫助操作人員及時排查設備故障原因、避免損失進一步擴大，為停機操作和制定故障維修計劃爭取更多的時間。作為未來工作方向之一，我們將繼續改進所提方法，例如采用并行計算等方式提高方法計算效率，使其更加適用于實際應用。

表5 平均計算時間

圖9 測試數據計算時間

5.3 故障識別試驗

首先利用t-SNE方法將故障數據映射到低維空間進行三維可視化。不同形狀的散點表示不同類型的軸承數據，如圖10所示。如IF7為軸承內圈故障，故障尺寸為0.177 8 mm，其他具體圖例含義同表2。由圖10可知，相同類型的數據具有較好的聚集性，并且不同類型數據之間具有較好的分離性。這說明所建圖模型可以有效的對不同類數據進行描述，便于分類器進行分類實現故障診斷。

圖10 故障數據在特征空間的三維可視化

本文將構建的圖模型的50%用于訓練，50%用于測試。在數據集給出的4種負載條件下，對這10種類型數據利用KNN算法進行分類識別。本文方法在4種負載條件下，分別能達到90%、100%、100%、100%的準確率，為了進一步說明方法的優越性，本文同現有的方法進行了對比，本文方法和對比方法的試驗結果，如表6所示。從表6可知，對于在負載745.69 W、1 491.39 W、2 237.07 W條件下，本文方法的表現為最優。在0條件下，楊云博等[24-26]準確率優于本文方法。楊云博等利用一種改進的人工魚群算法對SVM進行優化，通過提取到的振動信號時域特征進行故障診斷；王太勇等提出了一種基于DenseNet的卷積核dropout(KD)智能故障診斷模型KD-DenseNet。其中，楊云博等只考慮了正常、外圈磨損、內圈磨損、滾動體磨損4種故障類型，并且并未將故障破壞尺寸考慮在內；王太勇等試驗中的故障類型只有7種，沒有考慮故障破壞尺寸為0.533 4 mm時的故障類型。本文認為故障類型數目的不同是導致這兩種方法準確率均高于所提方法的原因。并且，需要指出，在實際應用中，軸承的故障類型多種多樣，不同的工況會導致不同的故障破壞尺寸，本文方法將CWRU數據集中的10種故障類型全部考慮在內，分類更加精確，具有更好的現實意義。Zeng等提出了一種基于柔性凸殼的最大余量分類方法，并將其應用于滾動軸承的故障診斷，該方法將83%的數據用于訓練，17%的數據用于測試，在負載0時，其準確率達到99.60%，并且故障種類數為10。相比該方法，本文方法使用了更少的訓練數據，具有更高的實際應用價值。在實際應用過程中，軸承故障類型多種多樣，軸承故障類數據又通常難以大量獲取，因此，使用較少的故障數據完成故障診斷任務更加符合實際應用需求。因此，相比其他方法，本文所提方法具有較大優勢。

表6 故障識別結果

6 結 論

針對滾動軸承健康狀態監測問題，本文提出了一種新的圖譜域滾動軸承早期故障檢測與識別方法。通過基于PF分量構建圖模型進行圖譜分析，能夠有效獲取滾動軸承早期故障特征，通過計算相鄰圖模型間的相似性指標對軸承健康狀態進行低維度動態描述，實現了滾動軸承早期故障檢測。利用LMD對滾動軸承振動信號進行多尺度分析，通過圖建模的方式對得到的多尺度數據增強，能夠很好地用于故障類型描述，達到滾動軸承故障識別的目標。在后續工作中，將研究提高方法的運算效率，使其更好的應用于實際場合中。