動態模式分解及其在軸承早期故障診斷中的應用

文 明， 黨 章,3， 余 震， 呂 勇 魏國前

(1. 武漢科技大學 冶金裝備與控制技術教育部重點實驗室，武漢 430081；2. 武漢科技大學 機械傳動與制造工程湖北省重點實驗室，武漢 430081；3. 武漢科技大學 國家實驗機械教學示范中心，武漢 430081)

作為機器中常見的零部件，軸承被廣泛應用于各種機械設備中，大部分造成設備停機或者重大事故的滾動軸承在其服役早期通常含有不易發現的微弱故障[1]。對滾動軸承早期故障進行有效診斷不僅能夠排除設備故障隱患，而且對于提高設備運行的經濟性和可靠性意義重大[2]。傳感器采集得到的原始系統動力學響應信號通常包含表征機器健康狀況的特征信息，這些特征信息是滾動軸承健康狀態判斷和故障診斷的依據。早期故障軸承由于故障特征不明顯，設備運行過程中激起的能量低，且信號在傳遞過程中與設備本征信號及外部環境噪聲動態耦合，使得早期故障信號具有非線性、非平穩、故障特征微弱等特點,用傳統的故障診斷方法難以提取[3]。

近年來，針對復雜非線性、非平穩信號，已經提出了各種有效的時頻分析方法和模式分解方法，并應用于滾動軸承早期故障信號的故障診斷中[4]。短時傅里葉變換和小波分析方法基底具有構造簡單、對稱等特點，但其結構固定，缺乏變化，不適于分析復雜信號[5]。近十幾年出現的經驗模式分解(empirical mode decomposition, EMD)[6]、局部均值分解(local mean decomposition, LMD)[7]、變模式分解(variational mode decomposition, VMD)[8]及其變種模式分解方法[9-11]為信號分析帶來了方便，但EMD和LMD缺乏嚴格的理論公式推導，且分解得到的本征模式函數中存在低秩成分和稀疏成分混疊現象。當軸承存在早期故障時，由于其故障成分能量不高，噪聲相對于特征成分較強，其表征故障的低秩成分往往淹沒在強背景噪聲中，使得這些模式分解方法的有效性和效果受限，有時甚至不能提取出故障特征成分。VMD相對于EMD，LMD，其理論性更強，對噪聲魯棒性更好，但其分解效果依賴于參數的有效選取。這些信號分析方法在某些場景下具有良好的應用效果，但應用于強背景噪聲的早期故障信號時，使用效果將不可避免受到影響。

相比上述模式分解方法，動態模式分解(dynamic mode decomposition, DMD)[12]作為一種基于奇異值分解(singular value decomposition, SVD)的無方程數據驅動頻率分析方法，能夠精確提取復雜動態系統的時空特征。其最早起源于流體動力學領域，近年來被廣泛應用于流體場變量的特征提取、趨勢預測以及復雜系統控制等方面[13]。DMD算法結合了SVD和多模式表征的優點，通過一個動態數據分析模型，能夠將蘊含機械動力學特性的信號(數據)分解成一系列具有內在時空特征的單頻非正交模式[14]，每個模式對應于原始信號的低秩成分(有用成分)或稀疏成分(噪聲成分)，即使是具有強背景噪聲的早期故障信號，也同樣能將其分解為系列單頻模式，避免了EMD和LMD方法中的模式混疊問題。DMD方法將提取出的單頻模式與時間矩陣進行內積，得到的重構信號能夠高概率地恢復出原始信號的動力學特性。基于數據驅動的DMD方法為機械故障診斷中的復雜信號分析提供了一種新的研究思路和手段。

1 動態模式分解與多尺度排列熵

1.1 動態模式分解

假設信號采集系統以等時間間隔(Δt=ti+1-ti)對原始連續系統進行采樣,得到一維時間序列x=(x1,x2,…,xN)∈R1×N,以步長為1的定長窗口滑動選取x中元素并按行排列為Hankel矩陣，表示為

式中：m+n-1=N；X∈Rm×n。將X矩陣分解成兩個連續漢克子矩陣Xt，Xt+1，將其表示為時間上前后的兩個快照，表示為

Xt=[X1X2...Xn-1]
Xt+1=[X2X3...Xn]

(2)

估計一個最優線性算子A實現Xt到Xt+1的映射，使得

Xt+1=AXt

(3)

(4)

式中：U∈Rm×r,VT∈Rr×m分別為左右特征向量矩陣，且U×UT=I,V×VT=I;Δ為主對角線上含有不為零元素的矩陣，即特征值矩陣。通過最小化Xt+1與AXt之間的Frobenium范數誤差求解最佳近似解,表示為

(5)

則最佳近似解可由式(5)解得

(6)

(7)

(8)

式中，φ為標準DMD模式。令ωi=ln(λi)/Δt，則原信號X的近似重構解為

(9)

(1)將原始信號X按照式(1),式(2)選取維度m進行相空間重構，得到式(2)中兩個連續快照。

(4)按照式(3)～式(9)計算各模式和恢復信號。

1.2 截斷秩選取與多尺度排列熵

1.2.1 截斷秩的選取

DMD能夠精確地將信號分解成多個代表系統低秩成分或稀疏成分的模式，但選擇哪些模式用于恢復重構，即截斷秩r的選取是個重點問題。截斷秩r的大小決定了用于恢復的模式數量，如果r取值太大，可能將一些不必要的模式用于重構導致恢復信號包含大量噪聲，不利于特征提取；如果r太小，那么恢復的信號將丟失代表早期故障特征的有用成分。許多學者曾提出并使用若干基于奇異值的硬閾值選取方法，如奇異值差分譜[16]、奇異值能量譜[17]及奇異值標準譜[18]等，但相比而言，根據信號自身特點的軟閾值選取方法更加合理， Gavish等[19]提出了自適應截斷秩r的選取方法。對于包含噪聲和系統動力學信息的觀測矩陣

Y=X+ηN

(10)

式中：N為獨立分布的高斯白噪聲;X為低秩矩陣，η為噪聲水平系數。對于實際工程信號而言,信號噪聲大小往往是未知的，Gavish等指出，其幅值可以根據Y的奇異值中值進行估計。Gavish等應用奇異值中值λmedian估計噪聲矩陣ηN,以噪聲矩陣Y的相關參數λ為閾值，小于閾值的特征值所對應的模式去除，大于閾值的模式保留。閾值可以定義為

λ=ω(β)λmedian

(11)

ω(β)≈0.56β3-0.95β2+1.82β+1.43

(12)

式中：β為Y的行與列之比;ω(β)為閾值系數計算公式;λmedian為奇異值中值;λ為最佳截斷秩閾值。該文獻提供了一種根據信號自身特點選取截斷秩的方法，通過該方法確定截斷秩并恢復信號后，信號頻域可讀性很高，但低能量特征頻率(如內圈故障頻率邊頻)容易被當作背景噪聲去除。傳統以奇異值能量作為硬閾值選取截斷秩的方法往往不能準確的找到合適的閾值，有時為了避免部分特征頻率被“丟棄”選擇了過大的比值，使得故障識別過程中存在過多噪聲干擾。為了克服這個不足，根據式(11)確定截斷秩后，重新尋找一個自適應的奇異值能量閾值，把式(11)確定的奇異值能量與全部奇異值能量比值的均方根作為新的能量閾值，取得自適應能量截斷秩，該方法能夠更為準確的找到合理的比值，彌補了傳統奇異值能量閾值選取不夠精確的問題。假設通過式(11)得到的截斷秩為r，則修正后的奇異值能量比值表示為

(13)

式中，p=min(2m,n-1)。

1.2.2 多尺度排列熵

排列熵[20]是一種衡量時間序列復雜度的算法,代表了系統的復雜程度，但其無法在多個時間尺度上進行數據復雜度的測量。多尺度排列熵(multi-scale permutation entropy, MPE)[21]是在排列熵基礎上進行改進的一種反映信號不確定性和不規則性的非線性分析方法，具有更好的穩定性和更強的抗噪聲能力[22]，當存在動態觀測噪聲時，這種方法很有用。MPE是一種在多個時間尺度上測量系統復雜度的方法，它具有算法簡單，運算速度快等特點，其基本思想是將時間序列首先進行多尺度粗粒化，然后再進行排列熵值計算。熵值越小，表明系統時間序列趨于有序，呈現出周期性，如軸承、齒輪系統正常運轉產生的周期振動信號，反之趨于無序，如機械系統中廣泛存在的高斯白噪聲。

對于一個給定的時間序列x={xk,k=1,2,…,N},首先對其進行粗粒化處理。

(14)

式中：s為尺度因子，s=1,2,…；當s=1時，粗粒化序列即為原始序列,MPE與排列熵值等價;[N/s]為對N/s往下取整；當s≥2時，原始序列被粗粒化長度為[N/s]的時間序列。

對向量中的每兩個數據之間進行升序或降序考慮,每個重構向量有d!個可能的順序模式。多尺度排列熵值H最終被定義為

(15)

式中，pi為某一個特定的順序模式在所有重構向量中的頻率，圖1描述了針對一個時間序列pi的計算過程。正則參數ln(d!)保證了H(N,d,τ)的最大取值為1。

圖1 MPE計算流程Fig.1 Calculation process for MPE

本文在改進的全局DMD基礎上采用MPE算法計算DMD所得系列單頻模式的復雜度，通過設定閾值，將大于該閾值的模式作為噪聲過濾，小于閾值的模式認為包含有原始系統的動力學特征保留，并將其用于信號重構，以進行信號特征提取。本文技術路線如圖2所示。

圖2 基于MPE的DMD技術流程圖Fig.2 Flow chat of DMD based on MPE

2 算法在仿真信號中的應用

將本文提出的基于自適應截斷秩和MPE的DMD應用于Randall等[23]提出的軸承仿真信號模型，該模型由方程式(16)描述，根據軸承工作特點，當軸承內圈相對旋轉時，不同位置的轉子將承受不同的交變載荷，當軸承出現點蝕，膠合等缺陷時，軸承將產生瞬時沖擊信號，且振動信號以自身高頻振蕩衰減。

假設以等時間間距Δt進行信號采樣，S(t)為軸承的自然衰減頻率函數，Ak為第k個沖擊響應的幅值，N(t)為均值為零的背景噪聲，仿真信號模型表示為

(16)

式中：ak為第k個沖擊能量幅值;γ和φA為初始相位；fm為調制頻率；B為軸承系統相關的衰減系數；τk為波動引起的滯后時間；cA為一個常數，軸承內圈故障為fi，轉頻為fr，式(16)中的參數選取如表1所示。在軸承早期故障信號中，故障特征信息能量較低，容易被背景噪聲淹沒，因此仿真信號增加了-1 dB的高斯白噪聲，采樣頻率設置為fs=3 000 Hz。

表1 軸承仿真信號參數選擇Tab.1 Parameter selection for bearing simulation signal

試驗中設置m=2 000,文獻[24]建議采樣點數N>5d!=3 600,因此取N=4 000。軸承內圈故障仿真信號時域頻域如圖3所示。

圖3 仿真信號時域與頻域圖Fig.3 Simulation signal in time domain and frequency domain

從圖3可以看出，信號頻域存在尖峰干擾，故障頻率需要經過初步識別才能找到，且其轉頻和特征頻率的邊頻尋找困難，能量較小的邊頻淹沒在背景噪聲中。根據1.1節式(1)、式(2)構造漢克時間序列矩陣，按照步驟1～步驟3將Xt,Xt+1矩陣進行投影降噪，截斷秩根據提出的自適應能量軟閾值確定，使得信號在初步降噪時保留住特征頻率而減少噪聲的干擾，然后根據截斷秩對系統特征算子A進行特征估計。為了進一步濾除其中的一些噪聲，獲得更加可辨識的故障信號，采用1.2.2節中的MPE對模式進行閾值篩選，熵大于閾值的模式認為包含較多噪聲濾除，熵小于閾值的模式進行保留，通過最佳截斷秩篩選，共有r=309個模式，然后對篩選出的模式進行MPE值計算。

MPE是系統復雜性的度量，其值大小與尺度因子s和嵌入維度d有著明顯的關系。若d太小，重構向量中包含狀態太少，算法失去意義和有效性。若d太大相空間重構會均勻化時間序列，不僅導致計算費時，也不能反映出時間序列的微小變化。根據文獻[25]的建議，MPE中嵌入維數d取值為3～7，時延τ對時間序列的影響較小[26]，一般取τ=1。并據此經驗選取d=6,s=5,將選取的模式進行熵計算，得到各模式的熵分布如圖4所示。

圖4 各模式MPE值Fig.4 MPE value of each mode

取熵閾值為0.85，將熵值小于0.85的模式用于信號重構，共得259個模式，并對恢復的信號進行頻域觀察。恢復信號頻域如圖5所示。

圖5 恢復信號頻域Fig.5 Frequency spectrum of recovered signal

通過MPE篩選后，恢復信號頻域噪聲相比于原信號少很多，雖然頻域依舊存在若干不明干擾成分，但能量較強的特征頻率和能量較弱的邊頻都能夠找到，說明了基于MPE的DMD算法在早期故障信號特征提取中的可行性。為了表明該方法的故障特征提取效果，將該方法與EMD和SVD去噪效果進行對比。對信號采用EMD后，將含有特征頻率的本征模式函數進行恢復得到經驗模式恢復信號，SVD方法的截斷秩通過式(11)確定，通過各方法恢復的信號頻域如圖6所示。

圖6 不同方法處理后頻譜圖Fig.6 Frequency spectrum of signal processed by different methods

從各恢復信號頻域看，EMD恢復頻域依舊存在很多噪聲，且特征頻率邊頻辨識困難，沒有很好的從背景噪聲中分離。SVD雖然去噪效果很好，但是也去除了淹沒在噪聲中的邊頻，對信號進行過度去噪，導致丟失了有用成分。相比于這兩種傳統的去噪方法，MPE算法對模式進行閾值降噪，其恢復效果是這兩種傳統方法的折衷，通過閾值篩選過濾掉噪聲模式，進一步去除了背景噪聲的干擾，從恢復的頻域看，特征頻率及其邊頻能夠全部找到，且噪聲干擾更小，識別難度更低，準確度更高。通過正交投影降噪和MPE去噪，信號得到了良好的去噪效果，恢復信號特征頻率辨識更為簡單方便，綜合識別的準確性和難易程度看，本文方法在強背景噪聲下的軸承故障成分提取中具有一定的優越性。

3 算法在試驗信號中的應用

實際工程信號比仿真信號更為復雜，信號成分更加多樣，將本文方法應用于實際信號中，以驗證其工程應用效果。該節采用了Cincinnati大學公布的軸承試驗數據[27]。

Cincinnati大學試驗臺結構如圖7所示，4個Rexnord ZA-2115型雙列軸承安裝在一個由交流電機帶動的2 000 r/min的軸上，軸承每列有16個滾動體，軸承節距和滾動體直徑分別為28.15 mm，3.31 mm，軸承接觸角為15.17°，信號由裝在軸承箱的高靈敏度石英ICP加速度計(PCB 353B33)收集，采用NI DAQ Card 6062E收集信號。

圖7 Cincinnati 大學測試試驗臺Fig.7 Cincinnati University test bench

采用第二個數據集文件(2004.02.12.19.32.39)第一個信道進行驗證，該信道存在軸承外圈故障，由于噪聲的存在，故障頻率淹沒在背景噪聲中。原信號的采樣頻率fs=20 kHz,數據點數為N=20 480個，為了降低計算時間成本，在原信號基礎上按照等間距對信號進行重新采樣，采樣率降為一半，其頻域幾乎不變，因此設置新采樣頻率fs=10 kHz,N=10 240。由文獻[28]中的計算公式求解轉頻fr=33.3 Hz,外圈故障為fo=236.17 Hz。m取5 000，圖8為其信號時域和頻域圖。

圖8 Cincinnati試驗信號時頻域Fig.8 Cincinnati experimental signal in time domain and frequency domain

從圖8所示信號看，信號在低頻處出現若干尖峰，為轉頻所在頻段，信號特征頻率被背景噪聲淹沒無法識別。對構造的漢克矩陣進行初步投影降噪后，通過自適應截斷秩方法共獲得752個模式。利用MPE對這些模式進行熵計算，各模式熵值如圖9所示。

圖9 各模式MPE值Fig.9 MPE value of each mode

設置閾值為0.85，將高于閾值的模式進行濾除，剩余得到289個模式。將這些模式進行重構恢復，得到恢復信號，其頻域如圖10所示。

圖10 MPE恢復頻域Fig.10 Recovered frequency spectrum with MPE

雖然恢復頻域中仍然存在噪聲，但是通過r階投影去噪和MPE過濾后，相比于原信號頻域，噪聲已經大為減少，頻域中僅存在少量無法辨認的高能頻率，且故障頻率的二倍諧頻也被找到，故障頻率及其諧頻已經易于辨識。因此結合自適應截斷秩和MPE的動態模式分解方法在軸承故障提取中具有一定的可行性及有效性。同樣將EMD，SVD應用于該信號進行提取效果對比，各頻域如圖11所示。

圖11 各方法處理后頻域Fig.11 Frequency spectrum of signal processed by different methods

從圖11可知，相比于原始頻域，EMD雖然能夠提取特征頻率，但降噪效果不明顯，背景噪聲圍繞在特征頻率附近,沒有得到有效去除。在早期故障信號的強噪聲背景下，通過SVD方法恢復的頻域圖中，特征頻率及其諧波被當作噪聲去除，導致故障特征無法提取。MPE能夠根據熵值篩選出不同的單頻模式分量，熵值大于閾值的模式被視為噪聲去除，低于閾值的模式被視為有用成分保留，在去除噪聲同時，能夠較為方便的找到特征頻率。雖然MPE不能去除全部噪聲，從特征識別結果中能夠看出該方法具有優于傳統去噪方法的優勢，具備一定的應用前景。總體來說，本文提出的基于自適應截斷秩和MPE的DMD方法能夠更為有效的提取出含有早期強背景噪聲故障信號中的故障特征頻率，具有較強的應用前景。

4 結 論

在軸承早期故障信號中，噪聲往往充斥在整個頻域，由于早期故障頻率能量較弱，噪聲相對較強，軸承故障成分常常淹沒于本底噪聲之中而難以直接識別。為有效剔除早期軸承故障信號中的強背景噪聲，準確提取出故障頻率，本文提出選取改進的全局算子DMD方法將信號分解成系列單頻模式，克服了傳統模式分解的一些局限，為后續噪聲濾除帶來方便。進一步提出通過自適應最佳截斷秩方法和MPE閾值法對原始系統模式成分進行過濾篩選，對含噪信號中的有序模式進行了高概率恢復。

本文通過對漢克時間序列矩陣進行r階正交投影，構造了去噪漢克矩陣，在初步去噪的同時解決自適應截斷秩選取的問題，使得信號在保留特征頻率的同時包含更少的噪聲。對系列表征動力學特征信息的模式分量使用MPE進行閾值篩選，進一步提取出含噪更少的模式。通過對篩選出的模式進行恢復，能夠更加方便簡單地提取出原信號的特征頻率。仿真信號和試驗信號的試驗結果表明，本文提的基于自適應截斷秩和MPE的DMD方法在強背景噪聲下的故障信號特征提取中具有良好的應用潛力。

雖然通過本文方法能夠較好地對早期軸承故障信號進行降噪并提取特征頻率，但重構信號依舊存在若干不明成分，同時，在MPE閾值的選取上，也需要針對主觀恢復效果進行人工選取。因此，基于典型工況下不同軸承信號模式分量的熵值信息值得進一步研究，以提出自適應熵閾值選取方法。

Vol．41 No．12 2022