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

喬志城，劉永強，廖英英

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

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

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

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

1 理論基礎(chǔ)

1.1 經(jīng)驗小波變換

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

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

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

(1)

(2)

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

經(jīng)驗小波的細節(jié)系數(shù)和近似系數(shù)分別為

(3)

(4)

結(jié)合以上公式，得到振動信號的重構(gòu)公式

(5)

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

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

1.2 最小熵解卷積

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

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

(6)

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

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

(7)

對兩邊求導(dǎo)可得

(8)

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

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

(9)

MED的目的是求目標函數(shù)最大時的最優(yōu)逆濾波器，對式(9)求偏導(dǎo)并讓其等于零。

(10)

聯(lián)立式(8)和式(10)可得

(11)

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

B=A*F

(12)

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

g=A-1f

(13)

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

2 故障診斷方法及步驟

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

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

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

(14)

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

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

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

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

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

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

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

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

3 仿真信號分析

為證明本文方法的有效性，由沖擊、噪聲及諧波三部分構(gòu)建滾動軸承單點損傷故障模型仿真信號進行驗證，如式(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為隨機數(shù)； 采樣時長為1 s；fo為故障特征頻率，取值64 Hz。圖2為仿真信號各組成成分。圖3為仿真信號時域圖及包絡(luò)譜。

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

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

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

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

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

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

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

表1 分量峭度值

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

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

表2 分量互信息與峭度值

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

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

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

4 實驗信號

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

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

表3 雙列圓錐型滾動軸承參數(shù)

4.1 內(nèi)圈故障

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

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

圖9 內(nèi)圈信號時域圖及包絡(luò)譜Fig.9 Time-domain diagram signal and envelope spectrum of inner race

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

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

表4 分量峭度值

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

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

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

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

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

通過對比表5中各分量的互信息與峭度值，發(fā)現(xiàn)分量E3互信息值與峭度值均較大，圖14為單獨對E3進行MED降噪及包絡(luò)分析結(jié)果。

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

表5 分量互信息與峭度值

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

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

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

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

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

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

內(nèi)圈故障實驗證明：與傳統(tǒng)方法對比，本文提出改進頻譜重劃分方法，能夠有效減少頻帶數(shù)量，方便后續(xù)處理；該方法能夠?qū)^多諧波與噪聲的互信息值最大分量與其他分量分離，在信號重構(gòu)時，排除互信息最大的分量，在互信息較小的分量中進行選擇，能夠優(yōu)選出包含故障信息的最佳分量。

4.2 外圈故障

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

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

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

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

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

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

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

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

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

表6 分量互信息與峭度值

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

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

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

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

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

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

表7 IMF分量峭度值

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

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

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

4.3 滾動體故障

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

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

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

圖26 本文方法結(jié)果Fig.26 The result of proposed method

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

5 結(jié) 論

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