基于自適應(yīng)自相關(guān)譜峭度圖的滾動軸承故障診斷方法

鄭近德 王興龍 潘海洋 童靳于 劉慶運

1.安徽工業(yè)大學(xué)機械工程學(xué)院，馬鞍山，2430322.安徽理工大學(xué)礦山智能裝備與技術(shù)安徽省重點實驗室，淮南，232001

0 引言

滾動軸承是機械設(shè)備中最容易發(fā)生故障的零部件。由于滾動軸承在旋轉(zhuǎn)機械中應(yīng)用廣泛，且是關(guān)鍵零件，因此，開展對滾動軸承故障診斷的研究具有重要意義。

故障軸承在運行過程中通常會產(chǎn)生沖擊，從而得到一系列非線性、非平穩(wěn)的振動信號[1]，該振動信號包含著豐富的故障信息，對該信號進行分析和處理，能夠達到診斷故障的目的。包絡(luò)解調(diào)分析是最常用的故障診斷方法[2]，其難點是尋找最優(yōu)解調(diào)頻帶。針對此問題，DWYER[3]提出了譜峭度方法，尋找最大峭度值所對應(yīng)的頻帶，并對該頻帶內(nèi)信號進行分析得到瞬態(tài)成分[4-5]。ANTONI[6-7]提出了快速譜峭度法，分別通過濾波器組結(jié)構(gòu)[6]和短時傅里葉變換[7]計算譜峭度。LEI等[8]提出改進的譜峭度方法，克服了強噪聲對快速譜峭度方法的影響，改進譜峭度法利用小波包變換從含噪信號中準(zhǔn)確檢測到瞬態(tài)成分，但是此方法在處理具有強非高斯噪聲的振動信號時容易失效。為此，BARSZCZ等[9]提出了Protrugram方法，該方法克服了強非高斯噪聲對峭度值的影響，但在此方法中，帶寬系數(shù)的設(shè)置需依據(jù)人為經(jīng)驗。針對此問題，WANG等[10]提出了增強的Kurtogram方法，其濾波信號包絡(luò)譜的峭度可用來顯示是否存在故障。上述方法在獲取最優(yōu)解調(diào)頻帶方面取得了較好的成果，但是當(dāng)故障頻率在頻譜中不明顯或噪聲干擾較大時，這些方法會受到限制，從而影響識別解調(diào)頻帶的準(zhǔn)確性。

針對上述問題，MOSHREFZADEH等[11]提出了自相關(guān)譜峭度圖(Autogram)方法，該方法以最大重疊離散小波包變換(maximal overlap discrete wavelet packet transform，MODWPT)為基礎(chǔ)，采用二叉樹結(jié)構(gòu)對頻譜進行分割，從而得到一系列解調(diào)頻帶[12]。Autogram方法通過限制原始信號中的非周期脈沖及噪聲來檢測周期性脈沖，能夠有效抑制非周期分量對實際故障頻率的影響，提高檢測最佳頻帶的準(zhǔn)確率。然而，Autogram方法在頻帶分割方面采用二叉樹結(jié)構(gòu)，導(dǎo)致計算最佳解調(diào)頻帶時誤差相對較大。基于此，本文提出了一種自適應(yīng)Autogram方法，即利用順序統(tǒng)計濾波(order statistic filtering，OSF)對信號傅里葉譜進行包絡(luò)[13-14]，克服了某一頻率范圍內(nèi)極大值過于集中的缺點，同時，采用平滑處理，消除了經(jīng)OSF包絡(luò)后產(chǎn)生的“平頂”數(shù)據(jù)。通過尋找包絡(luò)平滑譜的極大值點，選取距相鄰極大值中間位置最近的極小值點作為分割邊界，采用經(jīng)驗小波變換(empirical wavelet transform，EWT)對邊界內(nèi)信號進行重構(gòu)[15-17]，克服原始EWT分割方法的不足[18-19]，同時實現(xiàn)了頻譜的自適應(yīng)分割。

1 自適應(yīng)自相關(guān)譜峭度圖

Autogram方法是一種基于無偏自相關(guān)檢測最佳解調(diào)頻帶的新方法，它能夠有效檢測到淹沒在隨機噪聲中的解調(diào)頻帶及故障特征頻率。該方法以MODWPT為基礎(chǔ)進行頻帶劃分[20]，如圖1所示，各層帶寬固定。通過求取各個頻帶區(qū)域的峭度值，并將值填入到圖1中相對應(yīng)的位置，得到基于Autogram的峭度圖。

圖1 二叉樹劃分頻帶Fig.1 Frequency band is divided by binary tree

由于該方法在分割頻帶時邊界位置已經(jīng)固定，故在劃分過程中不能根據(jù)信號本身特征自適應(yīng)確定分割位置。本文采用EWT自適應(yīng)劃分頻帶，然而EWT中不同的傅里葉譜分割方式對分解結(jié)果具有很大的影響。對于含較多干擾成分的信號，傳統(tǒng)分割方法得到的結(jié)果往往是分割的頻帶過于集中，影響對其他頻率范圍內(nèi)信號的分析，為此，本文以改進的EWT劃分頻譜。

1.1 基于OSF的改進經(jīng)驗小波變換和頻帶分割

在劃分頻帶時，可根據(jù)原始信號傅里葉譜的具體特征尋找劃分邊界，以保證沖擊明顯的分量被包含到所選頻帶中，避免將最大沖擊處直接當(dāng)作邊界的情況出現(xiàn)。具體分割過程如下：①對信號進行傅里葉變換得到其傅里葉譜；②采用OSF對傅里葉譜進行包絡(luò)，得到包絡(luò)譜；③采用平滑處理對包絡(luò)譜進行優(yōu)化，去除一階不可微點，避免其對分割結(jié)果的影響，從而得到平滑包絡(luò)譜；④設(shè)置分割個數(shù)K，尋找平滑包絡(luò)譜的極大值點與極小值點位置，取所有極大值點的前K個極大值，然后將邊界定位到離兩相鄰極大值中間位置最近的極小值處，此邊界即為平滑包絡(luò)譜的分割位置；⑤將邊界進行歸一化；⑥選擇Meyer小波作為基函數(shù)，使用EWT對邊界內(nèi)信號進行重構(gòu)，完成EWT分解。

OSF包絡(luò)法本質(zhì)上是一個具有魯棒性的非線性濾波器，包括最大值濾波器、中值濾波器和最小值濾波器。與耗時的插值包絡(luò)不同，OSF耗時較短且能有效去除脈沖噪聲。該方法包含一個可變的滑動窗口，通過過濾窗口內(nèi)的數(shù)據(jù)，完成對數(shù)據(jù)的包絡(luò)。本文以最大值濾波器來估計上包絡(luò)。

在濾波過程中，首先確定窗口的寬度，然后取窗口內(nèi)數(shù)據(jù)的最大值，即

U(n)=max(Wn)

(1)

式中，U(n)為OSF濾波后結(jié)果；Wn為第n個窗口內(nèi)所有數(shù)據(jù)。

取數(shù)組D={2,3,4,5,6,4,8,3,2,1,7}，默認窗口寬度為3。為防止數(shù)據(jù)長度減小，使用鏡像延拓方法，即以原始數(shù)據(jù)左起第一個數(shù)據(jù)與右起第一個數(shù)據(jù)為基點，擴展長度為1，進行鏡像延拓。此時數(shù)組A變?yōu)閧3,2,3,4,5,6,4,8,3,2,1,7,1}。根據(jù)窗口寬度及式(1)，取{3,2,3}中最大值3，移動窗口，取{2,3,4}中最大值4，{3,4,5}中最大值5，依次類推，得到新的數(shù)組{3,4,5,6,6,8,8,8,3,7,7}。具體濾波過程如圖2所示。

圖2 OSF濾波過程Fig.2 Filtering process through OSF

經(jīng)OSF處理后，包絡(luò)圖中將出現(xiàn)一些“平頂”數(shù)據(jù)，即為一階不可微點。為了對此過程進行詳細說明，以實測軸承數(shù)據(jù)包絡(luò)處理為例，所得結(jié)果如圖3所示，為了方便觀察，選取其中一些“平頂”數(shù)據(jù)進行標(biāo)注。

圖3 OSF包絡(luò)處理Fig.3 Envelope processing through OSF

采用平滑處理進一步優(yōu)化包絡(luò)數(shù)組。本文選取移動平均法進行平滑處理。移動平均法通過將每個數(shù)據(jù)點替換為跨度內(nèi)定義的相鄰數(shù)據(jù)點的平均值來平滑數(shù)據(jù)。此過程等效于低通濾波，即

(2)

式中，ys(i)為第i個點的平滑值；N為ys(i)兩側(cè)的相鄰數(shù)據(jù)點的數(shù)量；2N+1為跨度。

使用移動平均法時需要遵循以下規(guī)則：①跨度必須為奇數(shù)；②需要平滑的數(shù)據(jù)點必須在跨度的中心；③對于無法容納任一側(cè)指定數(shù)量相鄰的數(shù)據(jù)點，跨度需要進行調(diào)整；④當(dāng)無法定義跨度時，端點不平滑。

為了進一步說明移動平均法，取跨度值為5進行平滑處理。根據(jù)以上規(guī)則，前4個元素為

ys(1)=y(1)ys(2)=(y(1)+y(2)+y(3))/3ys(3)=(y(1)+y(2)+y(3)+y(4)+y(5))/5ys(4)=(y(2)+y(3)+y(4)+y(5)+y(6))/5

另外需要注意的是，ys(i)是指排序后的數(shù)據(jù)順序，而不一定是原始數(shù)據(jù)順序。采用該方法對圖3包絡(luò)數(shù)組進行處理，結(jié)果如圖4所示。由圖4可看出，平滑處理有效地去除了一階不可微點，克服了經(jīng)OSF包絡(luò)后無法進一步分割的缺點。

圖4 移動平均平滑處理Fig.4 Moving average smoothing

1.2 計算所產(chǎn)生信號的無偏自相關(guān)

根據(jù)軸承振動信號二階循環(huán)平穩(wěn)性的自協(xié)方差的周期性[21]計算分解信號的無偏自相關(guān)，公式如下：

Rxx(ti,τ)=E(x(ti-τ/2)x(ti+τ/2))

(3)

Rxx(ti,τ)=Rxx(ti+T,τ)

(4)

式中，E(·)表示期望值；τ為延遲因子；x為1.1節(jié)分解得到的信號。

無偏自相關(guān)基于信號平方包絡(luò)進行計算，即

(5)

τ=q/fs

(6)

式中，X為分解信號的平方包絡(luò)；M為信號總長度，q=0,1,…，M-1；fs為采樣頻率。

通過無偏自相關(guān)去除信號中不相關(guān)的分量(如噪聲與無關(guān)脈沖)，可提高解調(diào)頻帶內(nèi)信號的信噪比。經(jīng)過無偏自相關(guān)處理后，噪聲在很大程度上會被消除，輸出結(jié)果更精確。

由式(5)、式(6)可知，隨著τ的增大，用于計算無偏自相關(guān)的數(shù)據(jù)樣本將減少，導(dǎo)致最終結(jié)果的估計方差不足，因此選取無偏自相關(guān)結(jié)果的前部分進行峭度計算。

1.3 檢測最佳解調(diào)頻帶

準(zhǔn)確地檢測到最佳解調(diào)頻帶是軸承故障診斷的關(guān)鍵步驟。首先計算1.2節(jié)中所有包絡(luò)信號無偏自相關(guān)的峭度值，然后尋找最大峭度值對應(yīng)分量所處的頻帶，即為最優(yōu)解調(diào)頻帶。峭度可用來統(tǒng)計數(shù)據(jù)集峰值，能夠檢測相關(guān)信號的異常脈沖。傳統(tǒng)的峭度被定義為

(7)

式中，μx為信號x的均值。

為了量化頻帶內(nèi)脈沖信號的脈沖，對傳統(tǒng)的峭度方程進行修改，即

(8)

選取所得最大峭度對應(yīng)的頻帶，對頻帶內(nèi)信號進行進一步分析，繪制其平方包絡(luò)譜，識別故障特征頻率，從而完成對故障的診斷。

2 仿真信號分析

采用周期性沖擊仿真信號驗證本文方法的有效性，計算公式為

(9)

其中，阻尼系數(shù)g=0.03，共振頻率fn=12 500 Hz，采樣頻率為20 kHz。a(t)為周期性沖擊信號的基本函數(shù)，通過添加信噪比為-12 dB的高斯白噪聲，構(gòu)成特征頻率f1為0.1 kHz的仿真信號，其時域波形如圖5所示(其中紅線為沖擊信號)，圖6為其傅里葉譜。

圖5 仿真信號時域圖Fig.5 Time domain diagram of the simulated signal

圖6 仿真信號傅里葉譜Fig.6 Fourier spectrum of simulated signal

首先，采用所提方法對仿真信號進行分析。所得峭度圖見圖7，圖中顏色代表峭度值大小。當(dāng)分解個數(shù)為5時，最大峭度值所處頻帶的中心頻率為4114.5 Hz，帶寬為737 Hz。對該頻帶內(nèi)信號進行分析，得到其平方包絡(luò)譜，見圖8。由圖8可知，所選頻帶內(nèi)濾波信號的故障特征頻率f1及其倍頻2f1及3f1清晰可見，其周圍無其他分量干擾，實現(xiàn)了故障的準(zhǔn)確診斷，進一步說明了本文方法在檢測最佳解調(diào)頻帶方面的可行性。

圖7 基于自適應(yīng)Autogram所得仿真信號的峭度圖Fig.7 Kurtosis diagram of the simulated signal based on adaptive Autogram

圖8 平方包絡(luò)譜(自適應(yīng)Autogram)Fig.8 Square envelope spectrum(adaptive Autogram)

為了對比，首先采用基于濾波器組結(jié)構(gòu)的快速譜峭度法對仿真信號進行分析，所得峭度圖見圖9。由圖9可知，最大峭度值對應(yīng)解調(diào)頻帶的中心頻率為7265.625 Hz，帶寬為156.25 Hz。計算該頻帶內(nèi)信號的平方包絡(luò)譜，結(jié)果如圖10所示，相比于自適應(yīng)Autogram法所得結(jié)果，其特征頻率被完全淹沒，采用該方法無法得到所需故障特征頻率，說明該方法無法有效檢測到最佳解調(diào)頻帶，導(dǎo)致故障特征頻率不明顯。

圖9 基于快速譜峭度(濾波器組結(jié)構(gòu)) 的仿真信號峭度圖Fig.9 Kurtosis diagram of a simulated signal based on fast kurtogram(filter bank structure)

圖10 平方包絡(luò)譜(濾波器組結(jié)構(gòu))Fig.10 Square envelope spectrum(filter bank structure)

為了對比，采用基于短時傅里葉變換的快速譜峭度計算方法對仿真信號進行分析，所得峭度圖見圖11。由圖11可知，該方法所得最佳解調(diào)頻帶的中心頻率為9375 Hz，帶寬為39 Hz，位于分解等級的第8級。對該頻帶內(nèi)的信號進行分析，得到平方包絡(luò)譜如圖12所示。如同濾波器組結(jié)構(gòu)法處理后得到的結(jié)果，其故障特征頻率仍然被淹沒，說明此方法所得頻帶并非最佳解調(diào)頻帶。

圖11 基于快速譜峭度(短時傅里葉變換) 所得仿真信號的峭度圖 Fig.11 Kurtosis diagram of a simulated signal based on fast kurtogram(short time Fourier transform)

圖12 平方包絡(luò)譜(短時傅里葉變換)Fig.12 Square envelope spectrum (short time Fourier transform)

最后，采用原Autogram方法對仿真信號進行分析，得到峭度圖見圖13。由圖13可知，所得頻帶的中心頻率為4687.5 Hz，帶寬為625 Hz。對該頻帶內(nèi)信號進行分析，得到其平方包絡(luò)譜如圖14所示。由圖14可知，相比于快速譜峭度法，其一倍頻f1與二倍頻2f1較為明顯。相比于本文所提方法，其故障特征頻率不夠明顯，周圍其他干擾分量較多，尤其在其三倍頻3f1處，故障特征頻率完全被淹沒。以上對比結(jié)果說明了所提方法的優(yōu)越性。

圖13 基于原Autogram所得仿真信號的峭度圖Fig.13 Kurtosis diagram of a simulated signal based on original Autogram

圖14 平方包絡(luò)譜(原Autogram)Fig.14 Square envelope spectrum(original Autogram)

3 實驗分析

為了進一步驗證所提診斷方法的可行性與實用性，本節(jié)將對實際故障滾動軸承進行診斷分析。實驗數(shù)據(jù)來自安徽工業(yè)大學(xué)自制滾動軸承模擬故障試驗臺，如圖15所示。實驗軸承型號為6206-2RS1 SKF，滾子個數(shù)為9，內(nèi)徑30 mm，外徑62 mm。使用線切割技術(shù)對軸承進行切割。選取滾動軸承內(nèi)圈切割深度0.4 mm的故障進行實驗，如圖16所示。實驗過程中，設(shè)置載荷為5 kN，實際轉(zhuǎn)速為300 r/min，采樣頻率為8192 Hz，采樣時間為20 s。經(jīng)計算，軸承故障頻率fi=27.15 Hz。

圖15 滾動軸承模擬故障試驗臺Fig.15 Simulated fault test bench for rolling bearings

圖16 滾動軸承內(nèi)圈故障Fig.16 Inner ring failure of rolling bearing

首先，采用本文方法對圖17所示的實測信號進行分析。實測信號傅里葉譜如圖18所示，所得峭度圖見圖19。由圖19可知，當(dāng)分解個數(shù)為13時，檢測到最大峭度值對應(yīng)頻帶的中心頻率為2796.5 Hz，帶寬為2599 Hz。為了驗證所得頻帶的準(zhǔn)確性，對該頻帶內(nèi)信號進一步分析，得到其平方包絡(luò)譜，如圖20所示。由圖20可看出，無關(guān)分量對特征頻率干擾程度低，其故障頻率清晰可見，此結(jié)果有效驗證了自適應(yīng)Autogram方法的可行性。

圖17 實測信號時域波形圖Fig.17 Time domain waveform of measured signal

圖18 實測信號傅里葉譜Fig.18 Fourier spectrum of measured signal

圖19 基于自適應(yīng)Autogram所得實測信號的峭度圖Fig.19 Kurtosis diagram of a measured signal based on adaptive Autogram

圖20 實測信號濾波后平方包絡(luò)譜(自適應(yīng)Autogram)Fig.20 Square envelope spectrum of the filtered measured signal(adaptive Autogram)

為了對比，利用基于濾波器組結(jié)構(gòu)的快速譜峭度法對實測信號進行分析。所得峭度圖見圖21，檢測到最大峭度值對應(yīng)解調(diào)頻帶的中心頻率為2218.6667 Hz，帶寬為341.3333 Hz。對該頻帶內(nèi)信號進一步分析，得到其平方包絡(luò)譜，如圖22所示。由圖22可看出，相比于自適應(yīng)Autogram法，在二倍頻2fi處，周圍干擾成分較多，而三倍頻完全混淆在無關(guān)分量中，故障特征頻率不夠明顯。說明本文方法在檢測最佳解調(diào)頻帶效果上強于基于濾波器組結(jié)構(gòu)的快速譜峭度法。

繼續(xù)采用基于短時傅里葉變換的快速譜峭度法對實測信號進行分析。圖23為所得實測信號的頻帶分布圖，可知，最大峭度對應(yīng)頻帶的中心頻率為3328 Hz，位于分解等級的第4.5級。進一步分析此頻帶內(nèi)信號，得到其平方包絡(luò)譜，如圖24所示。由圖24可知，一倍頻fi處故障特征較為明顯，但其他倍頻被周圍分量淹沒，無法進行準(zhǔn)確診斷。與所提方法進行比較，其二倍頻、三倍頻處特征頻率不明顯，其他分量干擾較大。進一步說明了該方法所選頻帶并非最優(yōu)解調(diào)頻帶。

圖21 基于快速譜峭度(濾波器組結(jié)構(gòu)) 所得實測信號的峭度圖Fig.21 Kurtosis diagram of a measured signal based on fast kurtogram(filter bank structure)

圖22 實測信號濾波后平方包絡(luò)譜(濾波器組結(jié)構(gòu))Fig.22 Square envelope spectrum of the filtered measured signal(filter bank structure)

圖23 基于快速譜峭度(短時傅里葉變換) 所得實測信號的峭度圖Fig.23 Kurtosis diagram of a measured signal based on fast kurtogram(short time Fourier transform)

圖24 實測信號濾波后平方包絡(luò)譜(短時傅里葉變換)Fig.24 Square envelope spectrum of the filtered measured signal(short time Fourier transform)

最后，采用原Autogram對實測數(shù)據(jù)進行分析，圖25為所得峭度圖。由圖25可知，最大峭度對應(yīng)解調(diào)頻帶的中心頻率為2560 Hz，帶寬為1024 Hz。對此頻帶內(nèi)信號進行分析，得到圖26所示的平方包絡(luò)譜。由圖26所知，一倍頻fi處故障特征明顯。與自適應(yīng)Autogram法對比，其二倍頻、三倍頻特征頻率不夠明顯，周圍干擾分量較多，因此，該方法所得頻帶并非最佳解調(diào)頻帶。綜上對比，進一步驗證了所提自適應(yīng)Autogram方法的優(yōu)越性。

圖25 基于原Autogram所得實測信號的峭度圖Fig.25 Kurtosis diagram of a simulated signal based on original Autogram

圖26 實測信號濾波后平方包絡(luò)譜(原Autogram)Fig.26 Square envelope spectrum of the filtered measured signal(origial Autogram)

4 結(jié)論

針對自相關(guān)譜峭度圖(Autogram)無法根據(jù)信號本身特點自適應(yīng)劃分頻帶的缺點，提出了一種自適應(yīng)Autogram方法。新方法在劃分頻帶時，對原始信號的傅里葉包絡(luò)平滑譜進行劃分，通過更改分解個數(shù)以得到不同的分割邊界，采用經(jīng)驗小波變換重構(gòu)邊界內(nèi)分量，從而完成頻帶的自適應(yīng)分割。通過仿真信號與實測數(shù)據(jù)分析，并與現(xiàn)有快速譜峭度及Autogram進行對比，得到如下結(jié)論：

(1)與快速譜峭度法相比，所提方法能夠根據(jù)信號本身特點自適應(yīng)劃分頻帶，最大峭度值對應(yīng)頻帶內(nèi)信號的平方包絡(luò)譜中，故障特征明顯，其頻帶檢測效果強于快速譜峭度。

(2)與原Autogram相比，克服了最大重疊離散小波包變換以二叉樹結(jié)構(gòu)劃分頻帶的缺點，同時頻帶內(nèi)信號故障特征明顯，故障檢測效果明顯提高。

綜上所述，所提方法能夠根據(jù)原始信號傅里葉譜的特點自適應(yīng)劃分頻帶，同時故障特征更加明顯。但本文仍存在不足之處，如在順序統(tǒng)計濾波處理過程中，其窗口寬度需事先設(shè)定，該問題仍需要進一步研究。