改進經驗小波變換與最小熵解卷積在鐵路軸承故障診斷中的應用

喬志城，劉永強，廖英英

(1. 石家莊鐵道大學 機械工程學院，石家莊 050043；2. 石家莊鐵道大學 省部共建交通工程結構力學行為與系統安全國家重點實驗室(籌)，石家莊 050043；3. 石家莊鐵道大學 土木工程學院，石家莊 050043)

滾動軸承作為影響旋轉機械正常運行的關鍵部件，對其進行狀態監測與故障診斷十分必要。滾動軸承振動信號一般為非線性、非平穩信號，傳統的Fourier變換只適用于頻率不隨時間變化的線性平穩信號，且容易受噪聲干擾影響[1-2]。針對滾動軸承故障信號非線性、非平穩特性以及早期故障特征微弱，易受噪聲干擾的特點，國內外學者提出經驗模態分解(Empirical Mode Decomposition, EMD)的故障特征提取方法。EMD能夠將復雜的軸承振動信號自適應分解為若干個本征模態函數分量，實現微弱故障信號與噪聲的分離。但EMD本身存在嚴重的端點效應、模態混疊等缺陷[3-4]。文獻[5]提出在信號中加入不同幅值的高斯白噪聲的集合經驗模態分解(Ensemble Empirical Mode Decomposition, EEMD)方法，利用高斯白噪聲的均勻分布統計特性，使不同頻率尺度的信號具有連續性，從而有效抑制了模態混疊[6]。文獻[7]提出變分模態分解(Varianational Mode Decomposition, VMD)方法，該方法通過迭代方式搜尋變分模型最優解來確定各分量的中心頻率與帶寬，實現信號頻域剖分及各分量的有效分離。相比于EMD，VMD具有完善的理論基礎與噪聲魯棒性，能夠有效抑制模態混疊[8]。但在使用VMD時需要確定模態分解個數K，K過大或過小都會影響信號分解精度[9-10]。

文獻[11]提出了基于小波變換與窄帶信號分析理論的經驗小波變換(Empirical Wavelet Transform, EWT)。EWT與EMD，VMD等相比[12]，避免了模態混疊與端點效應，且沒有迭代過程，不受參數設置的影響，具有計算量小，分解速度快等優點[13]。但EWT在頻譜劃分過程中，會出現頻譜劃分過多或者劃分不合理等問題。針對這一問題，文獻[14]提出根據頻譜包絡的劃分方法，克服了頻譜劃分集中于高幅值頻段的問題，但需要人工對分量個數進行估計，當面對復雜的頻譜時仍會因為估計不準并產生分解過度或分解不足的現象；文獻[15]提出時頻峭度譜的頻譜劃分方法，但該方法對簡單頻譜有較好的效果，對復雜頻譜的效果有待提升。

本文提出基于互信息的改進EWT頻帶劃分方法，解決傳統EWT頻帶分解過度的問題，根據重構分量的互信息與峭度選擇最優分量進行信號重構，但此時的重構信號中包含較多噪聲，本文采用最小熵解卷積方法(Minimum Entropy Deconvolution, MED)對重構信號進行降噪處理，然后對降噪信號進行包絡分析，能夠有效診斷出鐵路貨車輪對滾動軸承早期微弱故障，具有工程應用價值。

1 理論基礎

1.1 經驗小波變換

經驗小波變換是一種建立在小波理論框架下的自適應小波分析方法。EWT首先對振動信號的頻譜進行分割并構造一組具有適應性的小波濾波器組，然后對不同頻率成分進行分析，提取具有緊支撐特性的調幅-調頻信號[16]。

首先對振動信號進行Fourier變換得到頻譜，將頻譜定義在[0 π]，然后分割為連續的N段，每個頻帶可表示為Λn[ωn-1,ωn]。

根據Little wood-Paley和Meyer理論，經驗小波的尺度函數φn(ω)和小波函數ψn(ω)在頻域的定義為

(1)

(2)

式中：β(x)=x4(35-84x+70x2-20x3)； 0<γ<1；τn=γωn。

經驗小波的細節系數和近似系數分別為

(3)

(4)

結合以上公式，得到振動信號的重構公式

(5)

式中： ＾為傅里葉變換； ˇ為傅里葉逆變換； -為復共軛； *為卷積。

EWT頻帶劃分方法-尺度空間法：將頻譜劃分為N個頻帶，首先在頻譜中確定N+1個邊界點，除端點0和π外，還需確定N-1個點。將除0和π外的極大值點按從大到小順序排列，設找到M個極大值點，此時對應兩種情況：①當M≥N時，此時有足夠數量的邊界點，只保留前N-1個點；②當M≤N時，此時邊界點數量不足，保留所有邊界點的同時，對N值進行重置。

1.2 最小熵解卷積

最小熵解卷積本質：通過迭代方式，尋找一個最優逆濾波器，使包含故障信息的沖擊脈沖信號得到增強，迭代終止條件為峭度值最大[17-18]。假設振動信號出現局部損傷時的振動信號為

y(n)=x(n)*h(n)+e(n)

(6)

式中：y(n)為輸出信號；x(n)為輸入信號；h(n)為傳遞函數；e(n)為隨機噪聲。

根據熵的概念，熵是表征時間序列復雜程度的物理量，當時間序列越復雜，熵值越大；當時間序列規律性越強時熵值越小。因此當包含軸承故障信息的周期沖擊脈沖序列增多時，熵值減小[19]。MED的核心就是尋找一個L階最優濾波器f，使輸出信號y(n)最大程度恢復到輸入信號x(n)，并使其熵值最小[20]，即

(7)

對兩邊求導可得

(8)

式中，L為逆濾波器長度。

以解卷積后得到的x(n)的范數作為目標函數去衡量x(n)熵值的大小。

(9)

MED的目的是求目標函數最大時的最優逆濾波器，對式(9)求偏導并讓其等于零。

(10)

聯立式(8)和式(10)可得

(11)

將式(11)寫成矩陣的形式，即

B=A*F

(12)

進行迭代計算得到的逆濾波器的矩陣g為

g=A-1f

(13)

式中：B為輸入與輸出的互相關矩陣；A為信號y(n)的L×L自相關矩陣。

2 故障診斷方法及步驟

傳統EWT通過尺度空間法對頻譜進行自適應劃分，得到初始分界點。但此時得到的分界點數量較多，將頻譜分成過多的頻帶，對后續分析帶來不便。本文提出基于互信息值的初始分界點進行篩選的方法，將分量互信息值大于均值的相鄰頻帶合并為一條頻帶，分量互信息值小于均值的相鄰頻帶合并為一條頻帶，以此來減少頻帶數量。

互信息發展于信息論中熵的概念，定義為兩個隨機變量間不確定度的差值，能夠很好地衡量兩隨機變量的相關程度，比相關系數更準確[21]。

MI(X,Y)=H(Y)-H(Y|X)

(14)

式中：H(Y)為Y的熵；H(Y|X)為已知X時Y的條件熵。當X與Y相關性越強時，條件熵值H(Y|X)越小，此時互信息值MI(X,Y)越大[22]。圖1為軸承故障診斷流程圖。改進經驗小波變換頻譜劃分流程：

步驟1對原始振動信號進行FFT(Fast Fourier Transform)，求得振動信號頻譜；

步驟2根據尺度空間法確定初始頻帶分界點，獲得初始劃分頻譜；

步驟3計算根據初始分界點劃分頻譜得到的各分量互信息值，若該分量互信息大于均值，則其與相鄰互信息同大于均值的分量合并，若其相鄰分量互信息小于均值，則該分量保持獨立；若該分量互信息值小于均值，則將其與相鄰互信息同小于均值的分量合并，若其相鄰分量互信息大于均值，則該分量保持獨立，根據分量合并情況重新確定分界點；

步驟4根據重新確定的分界點，對頻譜進行重新劃分，從而得到新的分解信號分量；

步驟5最優分量選擇——互信息值最大的分量中包含較多與故障信息無關的振蕩與噪聲成分，因此在選擇最優分量時，首先排除互信息最大的分量，當軸承發生故障時峭度值大于3，因此在剩余峭度值大于3的分量中，選擇峭度大于均值的分量進行信號重構。

圖1 軸承故障診斷流程圖Fig.1 Fault diagnosis flow chart

從圖1可知，本文提出軸承故障診斷方法的思路，即：首先，通過改進的EWT獲得振動信號的分解模態分量，然后根據分量的互信息與峭度優選出最佳分量進行信號重構，為了提高診斷效果，對重構信號進行MED濾波處理，最后對濾波信號進行包絡分析，可以明確診斷出軸承是否存在故障。

3 仿真信號分析

為證明本文方法的有效性，由沖擊、噪聲及諧波三部分構建滾動軸承單點損傷故障模型仿真信號進行驗證，如式(15)～式(18)。

沖擊信號x1

x1=e-αtAsin(2πf1t)

(15)

噪聲信號x2

(16)

諧波信號x3

x3=Csin(2πf2t)

(17)

仿真信號y

y=x1+x2+x3

(18)

式中：A為沖擊幅值，取值1.2 m/s2；B為噪聲幅值，取值3.0 m/s2；C為諧波幅值，取值1.5 m/s2；α為衰減率，取值800；f1為共振頻率，取值3 kHz；f2為諧波頻率，取值25 Hz；fs為采樣頻率，取值10 240 Hz；z為隨機數； 采樣時長為1 s；fo為故障特征頻率，取值64 Hz。圖2為仿真信號各組成成分。圖3為仿真信號時域圖及包絡譜。

圖2 仿真信號各組成成分Fig.2 Components of simulation signal

圖3 仿真信號時域圖及包絡譜Fig.3 Time-domain diagram and envelope spectrum of simulation signal

從圖3(a)能觀察到明顯的沖擊，但圖3(b)所示的包絡譜，故障特征頻率及其倍頻成分很不明顯，難以確診軸承故障。利用EWT尺度空間法對頻譜進行劃分，得到初始頻帶如圖4所示。

圖4 初始頻帶Fig.4 The initial frequency band

從圖4可知，由于分界點較多，頻譜被劃分成53條頻帶。若采用傳統方法，從EWT分解后的53個分量中選擇表1中峭度較大的若干分量進行重構。重構信號如圖5所示。

圖5 重構信號、降噪信號時域圖及包絡譜Fig.5 Reconstructed signal time domain diagram,denoised signal time domain diagram and envelope spectrum

對圖5(a)重構信號進行MED降噪，結果如圖5(b)所示，噪聲雖然得到削弱，沖擊成分變得更加明顯，但從圖5(c)所示的包絡譜中仍然不能確診軸承故障。因此，根據分量互信息值重新對頻帶進行劃分，如圖6所示。

表1 分量峭度值

從圖6(a)可知，只有分量1的互信息值大于均值，因此將其余小于均值的頻帶進行合并，如圖6(b)所示。頻帶數由53減少為2，可知經過頻帶重劃分解決了傳統EWT頻帶劃分過多的問題。重構分量見圖6(c)，根據表2中的分量峭度與互信息，選擇E2分量進行下一步分析。圖7為對E2進行MED降噪得到降噪后時域圖及包絡譜。

圖6 分量互信息、重劃分頻譜及分量時域圖Fig.6 Mutual information of component,repartition of frequency spectrum and time domain of reconstructed component

表2 分量互信息與峭度值

從圖7(a)可知，分量E2經MED降噪后周期沖擊得到明顯增強。從圖7(b)的包絡譜可以明顯觀察到軸承故障特征頻率的倍頻f0～7f0，由此可以確定該軸承存在故障。

圖7 降噪后E2分量時域圖及包絡譜Fig.7 Time-domain diagram and envelope spectrum of E2 after denoise

仿真信號分析證明：軸承故障信號經傳統EWT分解為較多的分量，根據峭度準則進行信號重構與MED降噪，降噪后的包絡譜不能診斷出軸承故障。根據本文方法，在初始邊界的基礎上，根據互信息值對頻帶進行重新劃分，有效減少頻帶數量，將諧波信息與故障沖擊信息分離，從包絡譜中可以診斷出明顯的軸承故障，提高了軸承故障診斷的準確率。

4 實驗信號

為了證明本文方法在實際應用中的有效性，采用貨車輪對軸承實驗臺進行驗證。圖8所示為實驗臺結構圖，實驗選用CA-YD-188型加速度傳感器，采用磁吸座的方式布置在靠近測試軸承的機架上，實驗采用197726型雙列圓錐滾子軸承，軸承參數如表3所示。

圖8 鐵路貨車軸承輪對實驗臺Fig.8 Railway freight car wheelset running test bed

表3 雙列圓錐型滾動軸承參數

4.1 內圈故障

采用含內圈剝離故障的輪對軸承進行實驗，輪對轉速為465 r/min，采樣頻率fs為12.8 kHz，采樣時長為1 s。根據理論計算，軸承內圈故障特征頻率fi為88.25 Hz，轉頻fr為7.75 Hz。內圈故障信號時域圖及包絡譜如圖9所示。

從圖9(a)能夠觀察到時域沖擊成分，但從圖9(b)中僅能觀察到內圈故障頻率的1倍頻，無法判定軸承故障。圖10為尺度空間法得到的初始頻帶。

圖9 內圈信號時域圖及包絡譜Fig.9 Time-domain diagram signal and envelope spectrum of inner race

圖10 初始頻帶Fig.10 The initial frequency band

從圖10可知，EWT將信號頻譜劃分成47個頻帶，表4為選擇的峭度值較大的分量進行信號重構，重構信號如圖11所示。

表4 分量峭度值

圖11 重構信號、降噪信號時域圖及包絡譜Fig.11 Reconstructed signal time domain diagram, denoised reconstructed signal time domain diagram and envelope spectrum

對圖11(a)重構信號進行MED降噪，得到圖11(b)降噪后時域圖，從圖11(c)中僅能勉強找到內圈故障頻率的1倍頻f0及轉頻fr，無法確診軸承存在內圈故障。圖12為根據本文方法對EWT頻帶進行重新劃分。

圖12 分量互信息及重劃分頻譜Fig.12 Band MI and redistributed spectrum

根據圖12(a)分量互信息得到圖12(b)重劃分頻譜。頻帶數量由47減少為8。圖13為重劃分分量時域圖。

圖13 分量時域圖Fig.13 Time-domain diagrams of components

通過對比表5中各分量的互信息與峭度值，發現分量E3互信息值與峭度值均較大，圖14為單獨對E3進行MED降噪及包絡分析結果。

從圖14(a)可知，降噪后故障信號中的噪聲得到削減，沖擊性增強，但從圖14(b)僅能發現內圈故障特征頻率的1倍頻fi與2倍轉頻2fr，不能確診軸承存在內圈故障。

表5 分量互信息與峭度值

圖14 分量E3降噪后時域圖及包絡譜Fig.14 Time domain diagram and envelope spectrum after E3 noise reduction

上述分析證明：互信息值最大的分量與原信號具有最大的相似度，包含較多的諧波與噪聲成分，即使峭度值較大也不利于故障的提取，因此在選擇分量進行重構時，排除互信息值最大的分量，在剩余分量中選擇峭度較大的分量進行重構。

根據最優分量選取原則，排除互信息值最大的E3，結合圖15，在剩余峭度值大于3的分量中選擇峭度值大于均值的E7，E8進行重構。重構信號及MED降噪結果如圖16所示。

圖15 分量信號峭度值Fig.15 Kurtosis of component signals

從圖16(b)可觀察到明顯周期沖擊，從圖16(c)可以觀察到內圈故障的多個倍頻fi～5fi及以轉頻fr=7.8 Hz為間隔的邊頻，據此可判定軸承存在內圈故障。

圖16 重構信號、降噪信號時域圖及包絡譜Fig.16 Reconstructed signal, time-domain diagram and envelope spectrum of denoise signal

內圈故障實驗證明：與傳統方法對比，本文提出改進頻譜重劃分方法，能夠有效減少頻帶數量，方便后續處理；該方法能夠將包含較多諧波與噪聲的互信息值最大分量與其他分量分離，在信號重構時，排除互信息最大的分量，在互信息較小的分量中進行選擇，能夠優選出包含故障信息的最佳分量。

4.2 外圈故障

采用含外圈剝離故障的貨車輪對軸承進行實驗驗證。采樣頻率fs為25.6 kHz，轉速465 r/min，時長為0.5 s。根據理論計算，軸承外圈故障特征頻率為66.75 Hz。圖17為外圈信號時域圖及包絡譜。

圖17(a)可觀察到故障沖擊被淹沒在噪聲中，圖17(b)所示包絡譜中找不到外圈故障特征頻率及其倍頻。證明僅靠時域圖和包絡譜不能診斷軸承故障。

圖17 外圈信號時域圖及包絡譜Fig.17 Time domain diagram and envelope spectrum of outer race

圖18為尺度空間法得到的40個初始頻帶。在初始分界點基礎上，根據互信息對頻帶進行重新劃分，如圖19所示。

圖18 初始頻帶Fig.18 The initial frequency band

圖19 基于分量互信息重新劃分頻帶Fig.19 Band MI and repartited band

根據圖19(a)將互信息值大于或者小于均值的相鄰頻帶進行合并，得到圖19(b)重劃分頻譜，從左至右依次為E1～E7，頻帶數量由40減少為7。圖20為重劃分頻譜時域圖。

圖20 分量時域圖Fig.20 Time domain diagrams of components

從圖20可知，分量E5，E6沖擊最明顯，依據最優分量選取原則，在表6中排除互信息值最大的E2，根據圖21從剩余峭度值大于3的分量中選擇峭度大于均值的E5，E6進行重構。圖22為重構信號及MED濾波后時域圖及包絡譜。

表6 分量互信息與峭度值

圖21 分量信號峭度值Fig.21 Kurtosis of component signals

圖22 重構信號、MED降噪后信號時域圖及包絡譜Fig.22 Reconstruction signal、time-domain diagram and envelope spectrum after MED denoise

從圖22(a)可知，雖然重構信號周期沖擊較明顯，但仍包含較多噪聲。對其進行MED降噪得到圖22(b)，降噪后信號沖擊性得到增強。從圖22(c)可知，外圈故障特征頻率的f0～6f0倍頻，可以判斷軸承存在外圈故障。

為了證明本文方法的優越性，與EEMD方法進行對比。采用EEMD對外圈故障信號進行分解，得到圖23的IMF分量時域圖(考慮篇幅原因，僅展示前6個IMF分量)。

圖23 IMF分量時域圖Fig.23 Time domain diagrams of IMFs

根據表7選擇峭度最大的IMF2，IMF3與IMF5分量進行重構。圖24為重構信號、降噪重構信號時域圖及包絡譜。

表7 IMF分量峭度值

圖24 重構信號、降噪后時域圖及包絡譜Fig.24 Reconstruction signal、time-domain diagram and envelope spectrum after denoise

對圖24(a)重構信號進行MED降噪，圖24(b)為降噪后時域圖，圖24(c)包絡譜中不能找到外圈故障特征頻率及其倍頻，不能確診軸承故障。

外圈故障實驗證明：本文方法能夠有效減少EWT分量數量，并能將故障沖擊成分與諧波和噪聲成分分離，準確診斷出軸承外圈故障。通過與EEMD方法對比，證明本文提出的改進方法具有更高的診斷效率與準確性。

4.3 滾動體故障

采用含滾子剝離故障的貨車輪對軸承進行實驗，輪對轉速為465 r/min，采樣頻率fs為25.6 kHz，采樣時長為0.5 s。根據理論計算，軸承內圈故障特征頻率fb為27.08 Hz。滾動體信號時域圖及包絡譜如圖25所示。

圖25 滾動體信號時域圖及包絡譜Fig.25 Time domain diagram and envelope spectrum of ball signal

從圖25可知，僅通過頻譜分析難以識別出滾動體故障，考慮到篇幅原因，僅對滾動體信號處理結果進行展示。圖26為采用本文方法對滾動體故障軸承處理得到的包絡譜。

圖26 本文方法結果Fig.26 The result of proposed method

從圖26可知，明顯的滾動體特征頻率fb=26.56 Hz及其多個倍頻，可以確診滾動體存在故障。同時證明本文提出方法對滾動體同樣適用。

5 結 論

針對傳統EWT劃分頻譜時得到頻帶過多，為后續最優分量選擇帶來困難的問題，本文在尺度空間法得到原始分界點的基礎上，提出基于互信息的改進頻帶劃分方法，有效減少了頻帶數量。得到的重構信號中往往包含噪聲干擾，使用MED對重構信號進行降噪，提高了信噪比，增強了信號的周期性和沖擊性。本文提出方法與EEMD方法對比，取得了更好的診斷效果，有較大的優勢，并將其應用于貨車輪對軸承的故障診斷。通過改進經驗小波變換的頻譜劃分方法，避免了方法本身的缺點，為EWT的頻譜劃分問題研究提供了一個新的思路，提高了其可推廣應用的范圍。