基于小波包頻帶稀疏編碼的非完備信息條件下軸承狀態識別

馬云飛，賈希勝，胡起偉，郭馳名，邢 鵬

(1.武警士官學校 軍械系，杭州 310023；2.陸軍工程大學石家莊校區 裝備指揮與管理系，石家莊 050003)

作為裝備傳動系統的常見組件，軸承在各類型裝備中發揮重要作用，及早發現軸承故障并進行相應維修能夠有效提高裝備的可用性和戰備完好性[1]。因此，利用故障預測與健康管理[2]技術對軸承進行退化狀態識別具有十分重要的意義。

目前，傳統的退化狀態識別模型建立在完備先驗知識的基礎上，即需要知道各退化狀態的分界點，得到各個退化狀態數據。但在實際中，很難準確地對軸承全壽命數據進行狀態階段劃分，導致完備先驗知識下的狀態識別模型適用性較差[3]。而新軸承在經歷短暫磨合期后，處于正常潤滑階段，以及軸承發生了較為嚴重的故障，處于嚴重摩擦階段的情況下，監測者很容易區分出上述階段，從而獲得大量這兩個階段的狀態信息。基于此，本文提出了基于正常潤滑和嚴重摩擦兩種狀態信息的非完備先驗知識狀態評估模型。

近年來，稀疏編碼理論廣泛應用于降噪[4]，圖像超分辨率重建[5]，圖像融合[6]。Wright等[7]根據稀疏字典的判別性，提出了稀疏表示分類方法(sparse representation-based classification,SRC)。此后，學者將SRC理論應用于狀態監測領域，開展了大量研究。Zhang等[8]提出基于壓縮感知的軸承故障診斷模型，只需要采樣和存儲少量壓縮觀測數據，直接使用這些數據即可進行故障診斷。Tang等[9]基于稀疏表示方法和隨機降維投影，對原始信號進行降維采樣和SRC故障識別。Yu等[10]將時域振動信號轉換到頻域以減少噪音的破壞，然后利用L1范數和L2范數的組合在字典上對樣本進行稀疏編碼。苗中華等[11]對人工軸承數據集的振動信號進行稀疏編碼特征提取，試驗結果表明稀疏特征具有良好的可分性。余建波等[12]等提出基于固有時間尺度分解和稀疏編碼收縮相結合的特征提取方法。王維剛等[13]將分類器訓練反饋結果融入字典學習過程，將重構誤差項、稀疏編碼判別項和分類誤差項進行整合優化。

盡管已有稀疏編碼方法在故障分類領域取得了很大進展，但是：① 直接對時頻域信號進行稀疏編碼不夠精細，不能完全反映信號本質特征；② 已有研究需要所有故障狀態的完備信息進行訓練，但在實際中完備信息較難獲得；③ 少有文獻將稀疏編碼應用于裝備全壽命數據評估。

針對上述問題，本文提出一種基于小波包頻帶的稀疏編碼方法，在此基礎上建立非完備信息條件下軸承退化狀態評估模型。選取美國辛辛那提大學智能維修中心(intelligent maintenance systems,IMS)的軸承數據對新的稀疏編碼效果進行驗證，仿真結果表明新稀疏編碼比傳統稀疏編碼能夠更早地發現軸承故障。

1 稀疏編碼原理

稀疏編碼(sparse coding)的基本原理是將原始信號表示為一系列預先選定的字典原子線性組合。假設原始信號y∈Rm，預先訓練得到冗余字典H={h1,h2,…,hn}，包含n個字典原子(冗余字典每一列稱為原子)。若將y看作字典原子的線性組合，則可將其分解為

(1)

式中，x=[x1,x2,…,xn]T為稀疏系數。當原始信號y和字典H已知，求解稀疏系數x的過程稱為稀疏編碼。為求解原始信號的最優稀疏表示，稀疏編碼通常需要解決以下兩方面問題：① 字典學習，稀疏編碼的前提是確定能表征任意信號的原子集合，即完備字典；② 系數求解，在給定完備字典H基礎上，求解稀疏系數，常用算法有基追蹤[14],正交匹配追蹤[15]，基于迭代回歸的凸松弛算法[16]等。

常見的字典可分為固定字典和自適應字典。其中自適應字典通過對獲得的原始信號樣本進行學習獲得，因此稀疏表示能力更強。本文采用K奇異值分解(K-singular value decomposition，K-SVD)算法進行自適應字典學習。關于K-SVD算法的實現細節可以參考文獻[17]。

2 基于小波包頻帶的稀疏編碼方法

考慮到原始時頻域信號并不能完全反映出信號本質特征，而通過小波包變換處理軸承數據，能夠盡可能多的獲取高低頻信息，更有利于在軸承退化過程中捕捉到細微的變化。如圖1所示，本文將小波包變換與稀疏編碼方法相結合，改進方法主要有3項內容。

圖1 基于小波包頻帶稀疏表示與灰色B型絕對關聯度的稀疏編碼方法Fig.1 The sparse coding approach based on wavelet packet sparse representation and absolute grey relational degree of B-mode

(1)小波包變換和最優頻帶選取

首先，對原始信號進行4層小波包變換，并根據每層結點的能量選取最優頻帶，最優頻帶能夠反映最主要的信號特征。假設得到的16個小波包頻帶為{w1,w2,…,w16}，則可以計算每個頻帶的能量比重為

(2)

式中，lwp為每層頻帶的長度，能量比越大表明該層頻帶包含的信號信息越多，因此根據τ(i)選出最優頻帶。

(2)小波包頻帶稀疏分解與壓縮重構

選擇軸承正常狀態數據作為訓練樣本，對每一個選定的最優頻帶，利用K-SVD字典學習算法分別訓練一個過完備字典，得到字典集合{H1,H2,…,H}(假設ε為選擇的最優頻帶數量)。計算新樣本特征時，選用一種基于Laplace先驗的多任務貝葉斯壓縮感知(Laplace prior-based correlated-sparse-block Byesian compressive sensing,Lap-CBCS)算法[18]，可根據多組數據之間的相關性聯合重構，重構精度較高。

具體的，對每個測試樣本進行小波包變換得到頻帶wpi，根據設定好的觀測矩陣Φ，得到觀測向量

μ=Φ·wi

(3)

之后對wpi利用訓練好的稀疏字典Hi進行稀疏分解。如果訓練樣本與測試樣本結構相似，那么稀疏分解會十分有效，重構誤差較小。稀疏分解過程如下

μi=Φi·wi=Φi·Hi·θi,i=1,…,

(4)

式中，θi為稀疏系數。由于μi和Φi·Hi已知，本文利用Lap-CBCS算法對多個最優頻帶進行聯合重構。對每個測試樣本，最優頻帶的重構誤差δi可以表示如下

(5)

(3)灰色B型絕對關聯度計算重構誤差向量距離

由于重構誤差向量是多維數據，為了對其進行降維，利用灰色B型絕對關聯度計算樣本與正常軸承重構誤差平均值之間的距離作為稀疏編碼值。灰色B型絕對關聯度(absolute grey relational degree of B-mode,AGRDB)[19]方法基于灰色關聯理論計算兩個客觀事物之間距離，目前廣泛應用于工程領域。假設計算兩個向量S1和S2之間的RAB，計算方法如下

(6)

其中

(7)

(8)

(9)

3 非完備信息條件下軸承退化狀態評估模型

如圖2，本文建立一種基于正常狀態和嚴重摩擦狀態的軸承退化評估模型，具體如下。

圖2 非完備信息條件下軸承退化狀態評估模型Fig.2 Flow chart of bearing degradation assessment under incomplete information conditions

(1)退化特征提取。本文采用小波包頻帶稀疏編碼方法進行退化特征提取，提取方法見第2章。在特征提取前，應先對正常軸承信號進行小波包分解，最優頻帶選擇，以及K-SVD字典學習。

(2)模型訓練。對正常信號和嚴重摩擦信號樣本進行小波包頻帶稀疏編碼，并計算RAB距離，得到軸承狀態邊界分別為RAB，normal,RAB,failure。

(3)狀態監測。對待測試信號采用相同的特征提取方法計算出RAB,sample，則當前退化狀態評價指標Eunknown可以表示為

(10)

由式(10)可知，Eunknown可以用來描述測試信號與正常軸承信號的差異程度：當Eunknown→0時，表示測試信號的狀態接近于正常；當Eunknown→1時，表示測試信號的狀態接近于嚴重摩擦。此外，可以預先設定報警的門限值Q：當Eunknown>Q時，則表明軸承退化嚴重，發布預警信息并制定維修計劃；若Eunknown≤Q，則繼續監測下一個樣本。

4 公開數據驗證

采用辛辛那提大學智能維修中心全壽命軸承數據進行方法驗證。該試驗裝置包含4個(Rexnord)ZA-2115雙排滾子軸承，在每個軸承的外殼上安裝了PCB 353 V33型號的高靈敏度石英ICP加速度計。總共進行3次試驗，數據描述如表1所示。軸的轉速為2 000 r/min，采樣頻率為20 kHz，數據采集設備采用NI公司DAQ 6062E采集卡。

表1 IMS軸承數據集參數Tab.1 The parameters of IMS bearing data sets

(1)小波包分解與最優頻帶選擇

選取3組試驗中的故障軸承數據進行分析，首先對其進行4層小波包分解，得到對應的頻譜和能量分布如圖3所示。按照能量從大到小排序，選擇其中最靠前的8個小波包頻帶作為最優頻帶得到表2所示結果。

表2 IMS軸承數據集小波包分解最優頻帶Tab.2 The best frequency bands of IMS bearing data sets after wavelet packet decomposition

(2)退化特征提取

對選取的最優頻帶分別訓練一個稀疏字典，由于全壽命數據前一部分數據認為是正常狀態，選取每組數據前200～前500個樣本作為訓練樣本。為獲取過完備字典，要求字典的規格為m×n(m>n)，本文選取的稀疏字典規格如表3所示，其中線性組合數量為稀疏表示選取的原子數量。為保證收斂，K-SVD學習的迭代次數設置為50。

表3 稀疏字典規格Tab.3 Sparse dictionary specification

訓練完成后，對全壽命數據每個數據樣本分別按照第2章所示方法進行特征提取，得到小波包分解8個頻帶的重構誤差。如圖4所示，以第二組數據集全壽命數據為例，采用模糊C算法對獲得的重構誤差進行聚類。

圖4 小波包頻帶重構誤差特征聚類效果Fig.4 The clustering effect of reconstruction errors for wavelet packet frequency bands

(a)第一組數據

首先將特征歸一化，并設置類別數為3，將所有特征值聚成3堆，其中類別1為正常狀態，類別2為中間狀態，而類別3為嚴重摩擦狀態。為保證聚類既考慮特征值又考慮時間順序因素，除了獲得的8通道小波包特征，還增加了時間因子序列{0/984,1/984,…,983/984}作為第9個特征。由聚類結果可知，重構誤差特征在沒有先驗知識的情況下，可以將全部樣本按照退化時間順序聚為3類。除了中間狀態和嚴重摩擦狀態臨界處有部分樣本混淆外，大致反映了軸承退化規律，由此可以驗證小波包頻帶重構誤差具有判別性。

(3)狀態監測與預警

得到稀疏編碼特征后，按照式(6)計算出軸承狀態邊界分別為RAB，normal,RAB，failure。其中，計算RAB，normal時選用正常狀態數據樣本，計算RAB，failure時選用嚴重故障數據樣本。樣本區間可根據全壽命退化數據的變化趨勢確定(例如軸承退化最后階段發生突變)。如圖5～圖7，數據集1和數據集2很容易根據曲線變化確定出正常運行和嚴重摩擦狀態，從而得到訓練樣本，樣本區間如表4所示。但觀察數據集3發現其退化趨勢不明顯，無法區分出正常運行和嚴重摩擦狀態邊界，該數據集不適合用本文方法進行狀態識別。

表4 訓練樣本區間Tab.4 The training sample interval

圖5 數據集1訓練樣本選取Fig.5 The training samples selection of data set No.1

圖6 數據集2訓練樣本選取Fig.6 The training samples selection of data set No.2

圖7 數據集3退化特征序列Fig.7 The degradation feature sequence of data set No.3

獲得待測試樣本數據后，可計算出該樣本稀疏編碼特征，并得到RAB,sample。根據式(10)可求出最終的退化狀態指標，并與閾值Q進行比較。本文設置兩級預警，其中黃色預警線對應閾值Q=0.07，表明軸承只出現了輕微異常，處于早期故障狀態；紅色預警線對應閾值Q=0.7，表明軸承已經發生了嚴重故障。如圖8和圖9所示，采用兩種方法計算數據集2的稀疏編碼特征：① 直接進行K-SVD字典學習與壓縮重構，即傳統稀疏編碼;② 采用小波包頻帶稀疏編碼，獲得對應重構誤差與AGRDB稀疏編碼。對兩種方法得到的退化特征，分別畫出對應的黃色預警線和紅色預警線。根據局部放大圖可以發現，利用傳統稀疏編碼方法得到的兩級預警點分別為580，959，而利用本文小波包頻帶稀疏編碼得到的兩級預警點分別為549，935，本文方法兩級預警臨界點均早于傳統方法，有利于軸承嚴重摩擦故障與早期故障的預警。此外，在圖8(a)和圖9(a)中700～800出現了一個較大波動，與其他特征對比，發現是數據集本身工況變化等原因造成的。。

(a)退化狀態監測

如果根據實際需要，兩級預警線閾值Q發生改變，根據兩種方法以及小波包能量退化曲線分別得到的3種預警點如表5和表6所示(當Q≥0.5時，設置紅色預警線；當Q<0.5時，設置黃色預警線)。發現小波包頻帶稀疏編碼的絕大部分預警臨界點小于傳統稀疏編碼特征和小波包能量特征。同理，對數據集1計算不同閾值Q下的預警臨界點得到表7。觀察兩組數據集結果，發現當閾值Q較小時，小波包頻帶稀疏編碼的預警臨界點可能會大于傳統稀疏編碼，這主要是因為此時軸承處于正常狀態，退化趨勢不明顯，小波包稀疏編碼無法反應出退化特征微小變化。此外，本文在原始數據基礎上增加高斯隨機白噪聲，噪聲等級按照信噪比計算

表5 數據集2紅色預警線臨界點對比Tab.5 The comparison of red warning curve critical points for data set No.2

表6 數據集2黃色預警線臨界點對比Tab.6 The comparison of yellow warning curve critical points for data set No.2

表7 數據集1預警線臨界點對比Tab.7 The comparison of warning curve critical points for data set No.1

(13)

式中，u和u′分別為原始信號和帶噪音信號。如圖10，得到6種噪音等級下的稀疏編碼特征趨勢，重點觀察兩級預警線區域的特征變化，發現噪音等級對稀疏編碼特征的影響不大，且不同特征曲線的整體趨勢基本一致，表明本文所提特征具有一定的抗噪性。

(a)退化狀態監測

圖10 不同噪音下小波包頻帶稀疏編碼特征趨勢Fig.10 The trend of sparse coding on wavelet packet subbands under different noises

此外，對比傳統稀疏編碼和小波包頻帶稀疏編碼特征計算需要的時間復雜度。編碼特征計算包括兩部分：稀疏字典訓練(K-SVD算法)和重構誤差求解(Lap-CBCS算法)。由表8可知，盡管小波包頻帶稀疏編碼進行字典學習的時間較短，但其重構時間較長。這主要是因為本文所提特征需要對小波包分解后的每個頻帶分別進行重構，盡管每個頻帶長度小于原始信號長度，但把若干個頻帶重構時間加和后大于對原始信號進行重構的時間。

表8 傳統稀疏編碼與小波包頻帶稀疏編碼時間復雜度對比Tab.8 The time complexity comparison of traditional sparse coding and sparse coding on wavelet packet frequency bands

5 結 論

本文提出一種基于小波包頻帶稀疏分解和灰色B型絕對關聯度的編碼方法，并構建了非完備信息條件下軸承退化狀態評估模型，根據稀疏編碼計算退化特征評價指標，通過與故障閾值比較進行預警發布。利用IMS中心的軸承公開數據進行仿真實驗，結果表明新稀疏編碼特征相較于傳統稀疏編碼特征能夠更早的發現故障，有利于及時對軸承進行預防性維修，提高軸承壽命。盡管新稀疏編碼特征的計算復雜度要高于傳統方法，但考慮到數據傳輸到狀態監測后臺可以利用的計算資源較多，且在某些特定場景需要盡早發現故障征兆，因此新方法在實際中存在一定的應用價值。