動態(tài)模式分解及其在軸承早期故障診斷中的應(yīng)用

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

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

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

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

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

1 動態(tài)模式分解與多尺度排列熵

1.1 動態(tài)模式分解

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

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

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

(2)

估計一個最優(yōu)線性算子A實現(xiàn)Xt到Xt+1的映射，使得

Xt+1=AXt

(3)

(4)

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

(5)

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

(6)

(7)

(8)

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

(9)

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

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

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

1.2.1 截斷秩的選取

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

Y=X+ηN

(10)

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

λ=ω(β)λmedian

(11)

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

(12)

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

(13)

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

1.2.2 多尺度排列熵

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

對于一個給定的時間序列x={xk,k=1,2,…,N},首先對其進行粗?；幚?。

(14)

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

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

(15)

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

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

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

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

2 算法在仿真信號中的應(yīng)用

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

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

(16)

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

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

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

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

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

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

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

取熵閾值為0.85，將熵值小于0.85的模式用于信號重構(gòu)，共得259個模式，并對恢復(fù)的信號進行頻域觀察?；謴?fù)信號頻域如圖5所示。

圖5 恢復(fù)信號頻域Fig.5 Frequency spectrum of recovered signal

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

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

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

3 算法在試驗信號中的應(yīng)用

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

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

圖7 Cincinnati 大學(xué)測試試驗臺Fig.7 Cincinnati University test bench

采用第二個數(shù)據(jù)集文件(2004.02.12.19.32.39)第一個信道進行驗證，該信道存在軸承外圈故障，由于噪聲的存在，故障頻率淹沒在背景噪聲中。原信號的采樣頻率fs=20 kHz,數(shù)據(jù)點數(shù)為N=20 480個，為了降低計算時間成本，在原信號基礎(chǔ)上按照等間距對信號進行重新采樣，采樣率降為一半，其頻域幾乎不變，因此設(shè)置新采樣頻率fs=10 kHz,N=10 240。由文獻[28]中的計算公式求解轉(zhuǎn)頻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所示信號看，信號在低頻處出現(xiàn)若干尖峰，為轉(zhuǎn)頻所在頻段，信號特征頻率被背景噪聲淹沒無法識別。對構(gòu)造的漢克矩陣進行初步投影降噪后，通過自適應(yīng)截斷秩方法共獲得752個模式。利用MPE對這些模式進行熵計算，各模式熵值如圖9所示。

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

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

圖10 MPE恢復(fù)頻域Fig.10 Recovered frequency spectrum with MPE

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

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

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

4 結(jié) 論

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

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

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

Vol．41 No．12 2022