加窗插值快速傅里葉變換在滾動軸承故障診斷中的應用

李心一 謝志江 羅久飛

1.重慶大學機械工程學院，重慶，400044 2.重慶郵電大學先進制造工程學院，重慶，400065

0 引言

滾動軸承是旋轉機構的重要部件之一，其運行狀態的好壞直接影響旋轉機構的工作狀態和使用壽命，因此，滾動軸承故障診斷具有重要的工程現實意義。

在滾動軸承故障診斷中，由損傷點產生的低頻振動稱為軸承的故障特征頻率，采用頻譜分析法進行軸承診斷時，通過查看軸承振動信號頻譜中是否包含故障特征頻率，可以判斷軸承故障類型。研究表明，在滾動軸承故障診斷中，由故障產生的振動信號處于強噪聲背景，且往往出現不同程度的調制現象，因此，準確地提取故障特征頻率是判斷軸承故障類型的關鍵。Hilbert變換是一種常用的信號解調方法，可用于軸承信號的解調，但有一定局限［1］，文獻［2-3］提出了基于迭代Hilbert變換的滾動軸承故障診斷方法。Teager能量算子是在研究非線性語音建模時用于分析和跟蹤窄帶信號的能量所引入的算法，具有比Hilbert變換更好的解調效果，非常適合檢測信號中的沖擊成分，并廣泛應用于軸承故障診斷中。文獻［4］利用Teager能量算子提取軸承故障特征頻率，判斷軸承故障類型。文獻［5-6］將經驗模式分解與Tea?ger能量算子解調方法相結合，用Teager能量算子對單分量的固有模態函數進行解調后分析，從而可有效判斷軸承和齒輪的故障類型。文獻［7］利用Teager能量算子增強倒階次譜提取軸承微弱故障特征。在強噪聲背景下，為了能有效提取滾動軸承故障特征頻率，文獻［8-9］將最小熵解卷積（minimum entropy deconvolution，MED）應用于滾動軸承和齒輪箱的故障診斷，獲得了良好的效果。文獻［10］對MED在滾動軸承故障特征提取中的應用進行了詳細的分析。文獻［11］將Teager能量算子與MED相結合用于滾動軸承的故障診斷中，首先利用MED對振動信號進行降噪預處理，然后通過Teager能量算子求得包絡譜，可有效提取軸承故障特征。

加窗插值快速傅里葉變換（fast Fourier trans?form，FFT）方法最早用于提高FFT的計算精度［12］，后廣泛應用于電力諧波分析中。文獻［13-14］對基于最大旁瓣衰減窗（maximum sidelobe decay window，MSDW）插值FFT做了大量的研究，此類方法較為簡單，易于實現。

為了提高滾動軸承故障頻率的識別精度，本文針對軸承故障信號的特點，并結合MED和Tea?ger能量算子解調的優點，采用D.Belega加窗插值方法，利用常用的MSDW對解調信號進行加權處理并進行頻譜校正。

1 最小熵解卷積

當滾動軸承發生故障時，采集到的信號可表示為

由式（1）可以看出，故障軸承的振動信號 y(n)是故障軸承脈沖沖擊信號u(n)、背景噪聲n(n)與傳遞路徑h(n)的卷積結果。當軸承發生故障時，由于受隨機噪聲、傳遞路徑等的影響，輸出信號y(n)往往不能真實地反映輸入信號u(n)的特性。由文獻［15］可知，熵最小時，信號的峭度最大，所以MED核心思想是利用熵最小原理構造L階逆濾波器 fi(n)，使輸出信號 y(n)經此逆濾波還原為輸入信號u(n)，從而突出沖擊脈沖，即

MED求解方法主要有目標函數法和特征向量法，具體計算步驟見文獻［16］。

由以上分析可知，通過選擇合理的濾波器參數，采用MED對信號進行處理，可降低背景噪聲的干擾，增強信號沖擊成分，非常適合軸承振動信號的濾波預處理。

2 Teager能量算子解調

頻率為 f0的單頻信號以采樣頻率 fs進行采樣，得到離散時間信號

式中，A為幅值；φ為初相位。

則信號 x(n)的Teager能量算子ψ定義為［17］

由式（2）可以看出，Teager能量算子的輸出與信號的頻率和幅值的平方成正比，因此Teager能量算子的輸出對瞬態沖擊成分更敏感。與Hilbert變換相比，Teager能量算子計算簡單，只需計算相鄰3個采樣點，因此具有良好的時間分辨率。

軸承工作時由于局部損傷（剝落、腐蝕、裂紋、擦傷等），振動信號具有周期性沖擊的特點，而Teager能量算子具有增強信號瞬態特征的優點，因此非常適合檢測信號瞬態成分。

在實際工程應用中Hilbert變換和Teager能量算子是2種常用的解調方法，文獻［18］對這2種方法進行了詳細的對比。與Hilbert變換相比，Tea?ger能量算子具有運算量小和解調精度高的優點，因此本文選用Teager能量算子用于軸承振動信號的解調分析。

3 解調譜頻譜校正

最大旁瓣衰減窗是具有良好旁瓣性能的余弦組合窗，其表達式如下：

其中，I為窗函數的項數；ai為窗函數系數并滿足約束條件：

可計算得到2～5項最大旁瓣衰減窗函數系數見表1。由表1可知，Hanning窗是2項最大旁瓣衰減窗，適用于最大旁瓣衰減窗插值FFT方法。圖1所示為I取2、3、4、5，窗長和采樣頻率均為64時，窗函數的幅頻響應曲線，其中橫坐標以頻率點B表示，B等于采樣頻率除以采樣點。由圖1可知，隨著項數的增加旁瓣電平變小，旁瓣衰落速率增大。但是，MSDW的主瓣寬度隨著項數的增加而變寬，從而降低了頻率分辨率。

表1 最大旁瓣衰減窗函數系數Tab.1 Coefficients of MSDW

圖1 2～5項MSDW窗幅頻特性Fig.1 Amplitude-frequency properties of 2～5 term MSDW

文獻［13-14］對基于最大旁瓣衰減窗插值FFT法的頻率估計進行了詳細闡述，與其他經典窗FFT插值法相比，此類算法具有解析解，并且無需計算一元超越方程，插值公式較為簡單也更易實現。不失一般性，考慮式（3）所示的單頻信號，用長度為N的Hanning窗w(n)對該信號進行加權處理，得到加窗信號xw(n)=x(n)w(n)。根據卷積定理，其加窗離散傅里葉變換為

其中，W(·)為Hanning窗的離散時間傅里葉變換，當 N＞＞1時，其表達式為

若只考慮正頻部分，則Xw(k)可表示為

然而，因非同步采樣等原因，在實際情況中k0往往不是整數，因此，信號的精確頻率 f0可表示為

式中，l為峰值頻點k0附近抽樣幅值最大所對應的譜線號，為正整數；q為頻率偏差，其取值范圍為-0.5≤0＜0.5。

可見，得到精確頻率 f0的關鍵是求出頻率偏差值q。有多種方法可以求解q值，文獻［13］介紹了一種十分簡單的三譜線插值方法，該方法利用峰值頻點k0附近的3根離散頻譜的幅值求解q。令l-1為l左邊譜線，l+1為l右邊譜線，記 y1=|X(l-1)|，y2=|X(l)|,y3=|X(l+1)|，另記：

定義比值

根據式（8）和式（9），b可改寫為

把式（7）代入式（13），通過計算得

根據式（14），可求得

求得q后，根據式（9）可求出信號的精確頻率 f0。

4 軸承故障診斷方法

根據前文分析，基于漢寧窗插值FFT法的軸承故障診斷方法包括以下幾個主要步驟：①利用MED方法對振動信號進行降噪預處理；②計算軸承振動信號Teager能量算子輸出；③對Teager能量算子的輸出進行快速傅里葉變換，得到Teager解調譜；④對解調譜加Hanning窗，然后確定l、l+1、l-1；⑤根據三譜線幅值插值計算求解頻率偏值q；⑥得到校正后的軸承故障特征頻率值，判斷故障類型。詳細流程見圖2。

圖2 滾動軸承診斷流程圖Fig.2 Flow chart of rolling bearing diagnosis

5 軸承故障診斷實例

實例分析數據來自美國凱斯西儲大學軸承研究中心軸承故障模擬試驗臺［19］。內圈和外圈試驗采用的軸承為安裝于電機驅動端的6205-2RS JEM SKF深溝球滾動軸承，其軸承節徑dp=39.039 8 mm，滾動體直徑dB=7.940 04 mm，滾動體個數Z=9，接觸角 β=0°。電機負載為0～2.2 kW，轉速約為1720～1797r/min，采樣頻率為12kHz。

5.1 軸承內外圈故障診斷

軸承內圈故障點直徑0.177 8 mm，故障深度0.279 4 mm，電機轉速1 797 r/min。計算軸承內圈故障的通過頻率（the ball passing frequency inner race，BPFI）為162.185 973 Hz。選取數據段 2 049～3 072共1 024個采樣點進行分析，圖3所示為內圈故障信號。從圖3中可以明顯看出，軸承內圈故障信號具有周期性沖擊特性。

圖3 內圈故障信號圖Fig.3 Time domain waveform of inner race fault signal

圖4 所示為圖3信號的Teager解調譜。由圖4可知，Teager能量算子僅能提取出內圈故障特征頻率fBPFI=164.062 5 Hz（并且與理論計算值有較大誤差），其倍頻并不明顯。對圖3信號進行MED降噪預處理，參考文獻［8-10，20］選擇合適的MED參數（濾波器階數L=20，最大循環迭代次數M=30，迭代誤差e=0.01）。圖5所示為經MED濾波后的內圈信號解調譜。與圖3相比，內圈信號經MED預處理后，增強了信號中的沖擊成分，突出了軸承內圈故障特征頻率及其倍頻，但是其精度與理論計算值仍有一定誤差。

圖4 內圈信號解調譜Fig.4 Demodulation spectrum of inner race fault signal

圖5 MED濾波后內圈信號解調譜Fig.5 Demodulation spectrum of inner race fault signal filtered with MED

為了提高軸承故障頻特征頻率識別精度，計算經MED處理后的軸承內圈信號的Teager能量算子的輸出，并對該瞬時Teager能量序列進行傅里葉變換，得到解調譜。然后利用加Hanning窗插值FFT方法對解調頻進行頻譜校正，得到校正后的軸承故障特征頻率見圖6。由圖6可知，經過MED和頻譜校正后，軸承內圈故障特征頻率及其倍頻的精度都得到了有效提高，與理論計算值相吻合。

圖6 內圈故障特征頻率校正值Fig.6 Corrected BPFI of inner race

為了驗證本文所提方法的普遍適用性，排除人為選擇因素和數據樣本等因素的影響，分別選取軸承內圈故障點直徑為0.177 8 mm的隨機數據樣本100個和50個，并分別對應每個樣本數據長度為1 024和2 048。引入均方根誤差（root mean square error，RMSE）指標評估算法有效性，反映校正結果的精密度。RMSE定義如下：

式（16）中，N代表樣本數，fci是每個樣本的校正故障特征頻率值，fj表示計算的理論頻率值。表2、表3所示分別為100個樣本和50個樣本在故障點直徑為0.177 8 mm、電機轉速為1 721～1 797 r/min的RMSE仿真結果。

表2 100個樣本的RMSE（內圈 0.177 8 mm）Tab.2 RMSE of 100 samples（IR 0.177 8 mm）

表3 50個樣本的RMSE（內圈0.177 8 mm）Tab.3 RMSE of 50 samples（IR 0.177 8 mm）

對比表2、表3結果可以看出，未校正的故障頻率RMSE最大值為5.419 789（100個樣本）和2.982 938（50個樣本），最小值也接近2，存在很大誤差。總體上看，在不同電機運行狀態下，經頻譜校正后的故障特征頻率RMSE明顯優于直接由Teager能量算子解調的結果。

軸承外圈選擇故障點直徑為0.177 8 mm，故障深度為0.279 4 mm，故障點位于6點鐘方向。同外圈分析方法一樣，表4、表5所示分別為100個樣本和50個樣本在故障點直徑為0.177 8 mm、電機轉速為1 725～1 796 r/min的RMSE仿真結果。

表4 100個樣本的RMSE（外圈0.177 8mm）Tab.4 RMSE of 100 samples（OR 0.177 8 mm）

表5 50個樣本的RMSE（外圈0.177 8 mm）Tab.5 RMSE of 50 samples（OR 0.177 8 mm）

由表4、表5可知，與其他電機轉速相比，在電機轉速為1 773 r/min時，未校正值和校正值比較接近，這是因為此狀態下未校正特征頻率值為105.468 75 Hz，與理論計算特征頻率值105.930 11 Hz已十分接近。在此狀態下，校正后的RMSE值未有明顯減小。對于大多數情況，經頻譜校正后的RMSE值效果都較好。

5.2 窗函數影響

窗函數的選擇對加窗插值FFT的精度有較大影響。分析采用加窗插值FFT進行軸承故障診斷時窗函數對軸承特征頻率精度的影響，圖7、圖8所示為分別選用2～5項MSDW對軸承內外圈進行加窗插值FFT方法的RMSE值曲線（數據長度均為1 024），可以看出，RMSE值隨著窗函數項數的增加而增大。由圖1可知，2-MSDW具有最窄的主瓣寬度，主瓣寬度決定了加窗序列的頻率分辨率，因此較窄的主瓣寬度能帶來較高的頻率分辨率，從而提高軸承特征頻率的識別精度。

圖7 2～5項MSDW的RMSE（內圈 100樣本）Fig.7 RMSE of 2～5 term MSDW（IR 100 samples）

圖8 2～5項MSDW的RMSE（外圈 100樣本）Fig.8 RMSE of 2～5 term MSDW（OR 100 samples）

6 結論

（1）根據Teager能量算子和故障滾動軸承振動信號的特點，計算振動信號的Teager能量算子輸出，并對其進行傅里葉變換后得到解調譜，可以有效識別軸承的故障特征頻率，但在大多數情況下識別精度不高，與理論計算值存在偏差。

（2）本文介紹的三譜線插值FFT方法計算精度高且易實現，所用的最大旁瓣衰減窗為信號處理最常用的一類窗函數，算法易于編程，非常適合嵌入式系統的應用。

（3）主瓣較窄的窗函數具有較高的頻率分辨率，對故障特征頻率校正有一定影響。

（4）經軸承內外圈診斷實例分析，將MED、Teager能量算子解調和加窗插值FFT法相結合，在大多數情況下選取較少的分析點也能顯著提高軸承特征頻率的識別精度，性能上優于傳統的Teager能量算子解調方法。

［1］ 鄭曉溪，丁康，江利旗.希爾伯特變換解調分析在故障診斷中應用的局限性研究［J］.汕頭大學學報（自然科學版），1999，14（2）：40-46.ZHENG Xiaoxi，DING Kang，JIANG Liqi.Limitaton Studies on Application of Hilbert Transform Demodula?tion Analysis in Fault Diagnosis［J］.Journal of Shan?tou University（Natural Science），1999，14（2）：40-46.

［2］ 秦毅，任兵，李寧，等.改進的迭代希爾伯特變換及滾動軸承故障診斷應用［J］.振動與沖擊，2013，32（11）：83-88.QIN Yi，REN Bing，LI Ning，et al.An Improved Itera?tive Hilbert Transformation and Its Application in Roll?ing Bearing Fault Diagnosis［J］.Journal of Vibration and Shock，2013，32（11）：83-88.

［3］ 鄧四二，王勇，王恒迪.基于IHT的共振解調技術的滾動軸承故障診斷方法［J］.航空動力學報，2012，27（1）：69-74.DENG Sier，WANG Yong，WANG Hengdi.Resonance Demodulation Technique Based on IHT for Rolling Bearings Fault Diagnosis［J］.Journal of Aerospace Power，2012，27（1）：69-74.

［4］ 王天金，馮志鵬，郝如江，等.基于Teager能量算子的滾動軸承故障診斷研究［J］.振動與沖擊，2012，31（2）：1-5.WANG Tianjin，FENG Zhipeng，HAO Rujiang，et al.Fault Diagnosis of Rolling Element Bearings Based on Teager Energy Operator［J］.Journal of Vibration and Shock，2012，31（2）：1-5.

［5］ 李輝，鄭海起，楊紹普.基于EMD和Teager能量算子的軸承故障診斷研究［J］.振動與沖擊，2008，27（10）：15-17.LI Hui，ZHENG Haiqi，YANG Shaopu.Bearing Fault Diagnosis Based on EMD and Teager Kaiser Energy Operator［J］.Journal of Vibration and Shock，2008，27（10）：15-17.

［6］ 趙軍，崔穎，賴欣歡，等.隨機共振降噪下的齒輪微弱故障特征提取［J］.中國機械工程，2014，25（4）：539-546.ZHAO Jun，CUI Ying，LAI Xinhuan，et al.Weak Fea?ture Extraction of Gear Faults Based on Stochastic Res?onance Denoising［J］.China Mechanical Engineering，2014，25（4）：539-546.

［7］ 楊青樂，梅檢民，肖靜，等.Teager能量算子增強倒階次譜提取軸承微弱故障特征［J］.振動與沖擊，2015，34（6）：1-5.YANG Qingle，MEI Jianmin，XIAO Jing，et al.Weak Fault Feature Extraction for Bearings Based on an Or?der Cepstrum Enhanced with Teager Energy Operator［J］.Journal of Vibration and Shock，2015，34（6）：1-5.

［8］ ENDO H，RANDALL R B.Enhancement of Autore?gressive Model Based Gear Tooth Fault Detection Technique by the Use of Minimum Entropy Deconvolu?tion Filter［J］.Mechanical Systems&Signal Process?ing，2007，21：906-919.

［9］ SAWALHI N，RANDALL R B H，ENDO H.The En?hancement of Fault Detection and Diagnosis in Rolling Element Bearings Using Minimum Entropy Deconvolu?tion Combined with Spectral Kurtosis［J］.Mechanical Systems&Signal Processing，2007，21：2616-2633.

［10］ 龔廷愷，袁曉輝，王細洋.最小熵反褶積的數學形態法在滾動軸承故障特征提取中的應用［J］.中國機械工程，2016，27（18）：2467-2471.GONG Tingkai，YUAN Xiaohui，WANG Xiyang.Ap?plications of Mathematical Morphology Method to Fault Feature Extraction of Rolling Bearings Based on Minimum Entropy Deconvolution［J］.China Me?chanical Engineering，2016，27（18）：2467-2471.

［11］ 陳海周，王家序，湯寶平，等.基于最小熵解卷積和Teager能量算子直升機滾動軸承復合故障診斷研究［J］.振動與沖擊，2017，36（9）：45-50.CHEN Haizhou，WANG Jiaxu，TANG Baoping，et al.Helicopter Rolling Bearing Hybrid Faults Diagnosis Using Minimum Entropy Deconvolution and Teager Energy Operator［J］.Journal of Vibration and Shock，2017，36（9）：45-50.

［12］ JAIN V K，COLLINS W L，DAVIS D C.High Accura?cy Analog Measurements via Interpolated FFT［J］.IEEE Transactions on Instrumentation&Measure?ment，2007，28（2）：113-122.

［13］ BELEGA D，DALLET D.Frequency Estimation via Weighted Multipoint Interpolated DFT［J］.IET Sci.Meas.Technol.，2008，2（1）：1-8.

［14］ BELEGA D，DALLET D.Multifrequency Signal Anal?ysis by Interpolated DFT Method with Maximum Side?lobe Decay Windows［J］.Measurement，2009，42（3）：420-426.

［15］ WIGGINS R A.Minimum Entropy Deconvolution［J］.GeophysicalProspectingforPetrole，1980，16（1）：21-35.

［16］ 王宏超，陳進，董廣明.基于最小熵解卷積與稀疏分解的滾動軸承微弱故障特征提取［J］.機械工程學報，2013，49（1）：88-94.WANG Hongchao，CHEN Jin，DONG Guangming.Fault Diagnosis Method for Rolling Bearing’s Weak Fault Based on Minimum Entropy Deconvolution and Sparse Decomposition［J］.Journal of Mechanical En?gineering，2013，49（1）：88-94.

［17］ KAISER J F.On a Simple Algorithm to Calculate the‘Energy’of a Signal［C］//Proceedings of IEEE Inter?national Conference on Acoustics，Speech，Signal Pro?cessing.Albuquerque，1990：381-384.

［18］ POTAMIANOS A，MARAGOS P.A Comparison of the Energy Operator and Hilbert Transform Approach?es for Signal and Speech Demodulation［J］.Signal Processing，1994，37（1）：95-120.

［19］ Case Western Reserve University Bearing Data Cen?ter.Bearing Data Center Fault Test Data［EB/OL］.［2012-11-10］.http：//www.Eecs.Case.edu/laborato?ry/bearing.

［20］ JIANG R L，CHEN J，DONG G M，et al.The Weak Fault Diagnosis and Condition Monitoring of Rolling Element Bearing Using Minimum Entropy Deconvolu?tion and Envelop Spectrum［J］.Engineering Science Engineers，Part C：Journal of Mechanical Engineer?ing Science，2013，227（5）：1116-1129.