自適應譜峭度圖及其在滾動軸承故障診斷中的應用

何倩，易彩，吳文逸，林建輝

自適應譜峭度圖及其在滾動軸承故障診斷中的應用

何倩，易彩，吳文逸，林建輝

（西南交通大學 牽引動力國家重點實驗室，四川 成都 610031）

提出了一種新的自適應譜峭度圖方法。首先采用CEEMDAN方法分解信號的傅里葉頻譜，然后對得到的低頻IMF分量與殘差項進行分層累加計算，對每一層的頻譜趨勢尋找極小值點，使其作為頻譜劃分邊界，將劃分后的各頻段通過正交濾波器并利用峭度指標對劃分后信號中的故障信息進行評估，構造出一種新的自適應譜峭度圖。該方法能夠克服快速譜峭度圖因遵循固定劃分規則易使檢測出的頻帶不能包含完整的故障信息或是提取到無效分量的情況，并且分別利用軸承外圈故障的仿真信號和試驗信號，驗證了該方法的有效性。

故障診斷；自適應譜峭度圖；頻譜趨勢；快速譜峭度圖

滾動軸承廣泛應用于旋轉機械中，在現代機械設備中起著非常重要的作用[1]。在正常情況下，軸承在高速和重載下運行，即使輕微損壞對設備也非常危險，如果不能立即檢測到損壞，軸承將不穩定運行，嚴重時甚至損壞設備。因此，有必要實時監測軸承的運行狀態并進行故障診斷。

Dwyer[2]提出了一種非常強大的統計工具——譜峭度（Spectral Kurtosis，SK），可以代替功率譜密度來表示信號的高斯成分及其在頻域中的位置。雖然譜峭度的應用領域較廣，但其正式的數學定義尚未完善。Antoni等[3]對SK進行了深入而系統的研究，詳細闡述了相關理論，并給出數學定義，即能量歸一化四階譜累積量。此外，Antoni[4]證明了上述定義對噪聲中非平穩信號的特征檢測有很好的效果，提出快速譜峭度圖（Fast Kurtogram），并廣泛應用于故障診斷領域。

濾波器選擇的核心是保證所設計濾波器的邊界能夠完全覆蓋故障共振頻帶。由于故障信號的諧振頻率和帶寬可能是任意的，因此能夠靈活改變濾波器帶寬和中心頻率的濾波器是濾波器選擇的主要方向。由于快速譜峭度方法中由二叉樹或二叉樹結構的濾波器組構成的濾波器在每一級的帶寬和中心頻率是固定的，因此在工程實踐中可能無法獲得最佳共振帶。Barszc等[5]提出了根據固定窄帶內信號的平方包絡譜（Squared Envelope Spectrum，SES）的峭度選擇最佳解調頻帶，雖然這個過程可修改中心頻率，但由于不能調整帶寬的缺陷，結果仍然容易出現局部最優解。經驗小波變換是從采集到的信號中提取有意義模態的有力方法之一[6]，可將其視為一系列濾波器組，將一個信號分解為若干子信號。尺度空間[7]表示作為一種識別EWT（Empirical Wavelet Transform，經驗小波變換）邊界的無參數方法，由于其自適應性而被廣泛關注，但是相對來說，基于尺度空間的自適應經驗小波變換速度慢且不準確，且無效分量較多。

針對上述情況，本文提出基于頻譜趨勢的自適應劃分頻帶的方法，并將其與峭度指標結合，構造出自適應譜峭度圖算法，并通過仿真信號及實測信號證明了該方法的有效性。

1 基本理論

1.1 譜峭度

Dwyer首先將峭度應用于STFT的實部和虛部，并因此引入了頻域峭度（Frequency Domain Kurtosis，FDK）的概念[8]。SK最初定義為其頻率分量的峭度[9]，并與不同光譜頻率的振幅變化進行比較，因此該統計參數表示信號的脈沖性如何隨頻率變化。相反，Antoni[10]基于Wold-Cramer分解定義了SK，該分解將任何隨機非平穩過程描述為因果、線性和時變系統的輸出，即：

Antoni已經證明非平穩過程具有非高斯性是其基本特性[10]，然后將SK清楚地表示為CNS（Conditionally Non-Stationary，條件非平穩）過程的能量歸一化四階譜累積量，對應的表達式為：

1.2 快速譜峭度圖

Antoni和Randall[2]提出了快速譜峭度圖來表示頻譜分割平面上信號的峭度值，它可以被認為是中心頻率和帶寬的不同尺度下峭度值的函數。快速譜峭度的“快”體現在其快速頻帶分割技術中，該技術使用二元1/3樹結構將頻帶細分為多個合理的寬度。其實現過程如下：

步驟1：在不同層次上分割頻帶，然后形成一個具有不同帶寬和中心頻率的窄帶二維平面。

步驟2：過濾掉每個窄帶信號，并計算過濾信號的峭度值。峭度圖是使用這些峭度值填充二維平面生成。

步驟3：峭度值最大的頻帶位于峭度圖中作為最佳共振帶，如圖1所示，同時可以獲得所選共振帶的中心頻率和帶寬。

步驟4：對提取的濾波信號進行包絡分析，用于軸承故障診斷。

圖1 快速譜峭度圖示例結果

2 自適應譜峭度圖

2.1 頻譜趨勢自適應劃分頻帶

許多方法都可以估計給定信號的趨勢。傳統的方法有線性或多項式回歸分析、移動平均、基于趨勢估計的傅里葉濾波等，它們的共同局限性是需要提前給出基函數或跨度，且只能在給定情況下有效地估計它們。Wu等[11]提出的基于EMD（Empirical Mode Decomposition，經驗模態分解）的趨勢是一種數據驅動的、完全自適應的趨勢提取方法，首先對信號進行EMD分解，然后將最后幾個IMF和殘差進行合并得到EMD趨勢。EMD可以將復雜信號自適應地分解多個固有模態函數（Intrinsic Mode Function，IMF），但存在嚴重的模態混疊問題[12]。Torres等[13]提出了一種自適應高斯噪聲的集合經驗模態分解方法（Complete Ensemble Empirical Mode Decomposition with Adaptive Noise，CEEMDAN），進一步提高了分解信號的準確性和完整性，有效地克服了模態混疊問題。

基于上述分析，本文采用大程度減少模態混疊和不必要迭代篩選且數據分解更完整的自適應噪聲的集合經驗模態分解（CEEMDAN）來完成頻譜趨勢的尋找。

基于頻譜趨勢的自適應頻帶劃分的具體步驟如下：

（1）利用FFT（Fast Fourier Transform，快速傅里葉變換）獲得輸入信號的頻譜；

（2）應用CEEMDAN得到頻譜的IMF分量以及殘差項；

（3）將最后一項的殘差項作為第一項，剩余的IMF按照頻率成分從低到高排列，依次從第一項殘余項進行累加計算，每累加一個IMF作為一次頻譜趨勢，即每層頻譜趨勢；

（4）尋求每次頻譜趨勢的極小值位置，使其作為邊界點，分割輸入信號的頻譜。

圖2顯示了原始信號的時域波形和頻譜，對頻譜進行CEEMDAN分解的結果如圖3所示，圖4展示了部分頻帶劃分結果，其中藍色曲線表示殘余項與不同個數IMF疊加的結果、紅色曲線表示為頻譜邊界劃分的位置。

2.2 自適應譜峭度圖

針對快速譜峭度圖固定劃分所帶來的問題，本文提出基于頻譜趨勢的頻帶劃分方法和峭度結合求解諧振頻率和故障頻率的方法——自適應譜峭度圖，算法流程如圖5所示。

自適應譜峭度圖的軸承故障診斷的具體步驟如下：

（1）獲取軸箱上固定的振動加速度傳感器采集得到的振動加速度數據；

（2）計算相對應信號的傅里葉頻譜，并對頻譜進行自適應噪聲的集合經驗模態分解，得到各IMF分量及殘余項；

（3）將殘余項作為第一項，并按照頻率成分升序排列IMF分量，依次疊加組合殘余項與低頻IMF分量，得到每一層相對應的頻譜趨勢；

圖2 振動信號時域波形及頻譜

圖3 CEEMDAN分解結果

（4）計算每一層頻譜趨勢的局部極小值作為頻譜分割邊界，完成頻譜劃分；

（5）將劃分后的各頻段通過正交濾波器；

（6）計算每一層各頻段對應的峭度值，并將其鋪成二維平面，得到自適應譜峭度圖；

（7）選取峭度值最大的分量，進行Hilbert解調，進行故障分析，當包絡譜中包含有可識別的故障特征頻率及其倍頻時，則可認為該軸承存在相應的局部故障。

圖4 頻譜劃分結果

3 仿真信號分析

為驗證本文所提方法的有效性，構造一組單故障軸承仿真信號[14-15]如式（4）所示。

軸承故障仿真信號的參數如表1所示。信號采樣頻率為10 kHz，采樣時間為1 s。同時將信噪比為-10 dB的高斯白噪聲添加在仿真信號中。

仿真信號時域圖、頻譜圖、包絡譜如圖6所示。從圖6（c）可以看出，信號的中心頻率位于2500 Hz左右，共振帶長大衣位于2000～3000 Hz之間，與設置的仿真信號結果一致。從圖6（d）可以看到包含故障噪聲的故障特征頻率及其倍頻，也從側面印證了通過包絡譜可以更為清晰地了解故障情況。

圖5 自適應譜峭度圖算法流程圖

式中：A為第個故障脈沖的振幅；為阻尼比；為真實的故障脈沖周期，＝F/，F為采樣頻率，為軸承故障特征頻率；為隨機滑移效應；()為單位階躍函數；f為共振頻率。

表1 軸承故障仿真信號參數

采用基于頻譜趨勢的自適應譜峭度圖法處理該信號，結果如圖7所示。從圖7（a）可以看到，隨著IMF低頻分量的增加，每一層極小值點的個數也隨之增加，頻譜被劃分的數量也隨之增加。其中峭度值最大的頻帶位于層數4，中心頻率為2422 Hz，帶寬為834 Hz，提取此分量的波形如圖7（b）所示，對此分量做包絡分析，可以看到圖7（c）包絡譜中包含特征頻率，且可以找到2倍、3倍、4倍、5倍等高倍頻，故障信息明顯，且中心頻率f為2422 Hz，與實際中心頻率非常接近。考慮到噪聲的影響，該誤差可以接受。因此，該方法適用于含噪聲及周期性沖擊信號。

對上述仿真信號進行快速譜峭度圖計算，結果如圖8所示。從圖8（a）可看到峭度值最大的頻段位于層數4，中心頻率為2343.75 Hz，帶寬為312.5 Hz，與實際中心頻率相差接近150 Hz。提取此分量，波形如圖8（b）所示，對此分量做包絡分析，從圖8（c）包絡譜中可以看到特征頻率及2倍頻、3倍頻，說明快速譜峭度對于軸承故障診斷具有一定的效果。

通過兩種方法的對比分析，可以看到雖然兩者都能夠診斷出故障特征，但從整體上來說，通過使用基于頻譜趨勢的自適應譜峭度圖可以更好地反映出仿真信號信息，即計算所得的中心頻率及帶寬與實際仿真信號中心頻率更為接近，且幅值整體更為突顯，也印證了快速譜峭度提因遵循頻帶固定劃分，從而導致得到的最優頻帶包含不夠豐富的故障信息。

4 試驗分析

為進一步驗證該基于頻譜趨勢的自適應譜峭度圖的有效性，本節采用試驗數據對該方法進行驗證。

圖6 仿真信號時域圖、頻譜圖、包絡譜

圖7 采用基于頻譜趨勢的自適應譜峭度圖法處理的結果（一）

圖8 對仿真信號進行快速譜峭度圖計算的結果（二）

試驗臺運行的軸箱來自高速列車車輛，試驗數據來自臺架試驗軸箱軸承，該軸承為雙列圓錐滾子軸承，軸承型號為CEI-2962，其幾何尺寸如表2所示，試驗臺故障類型及位置如圖9所示，故障為人為設置，其中外圈損壞的深度為1 mm，寬度分別為1 mm、3 mm、5 mm，采樣頻率為10 kHz，軸旋轉頻率f為5.1 Hz，采樣周期為10 s。軸承外圈對應的特征頻率分別為83.33 Hz。

表2 雙列圓錐滾子軸承的基本參數

圖10為采集外圈故障信號的時域波形及其傅里葉譜和包絡譜。從頻譜中可以發現信號有較強的諧波干擾，干擾頻率為50 Hz，同時無法從振動信號的包絡譜中看到軸承外圈故障特征頻率，因此，需對信號進行處理，提取到故障特征信息，從而判斷故障類型。

使用基于頻譜趨勢的自適應譜峭度圖處理該信號，結果如圖11所示。從圖11（a）中可以看到峭度值最大的頻帶位于層數2，提取此分量的波形圖如圖11（b）所示，對該頻帶內的信號進行分析，可以看到圖11（c）中包絡譜包含特征頻率，且可以找到2倍、3倍、4倍等倍頻，故障信息明顯且豐富，精確地檢測到實際的故障特征。

當使用快速譜峭度圖處理上述信號時，所得結果如圖12所示，圖12（a）顯示了最大峭度值出現在中心頻率625 Hz、帶寬416.6667 Hz的頻帶中，濾波信號包絡和包絡譜分別如圖12（b）和（c）所示，與原始信號的包絡譜類似，從該包絡譜中無法檢測到任何軸承故障信息。因此，快速譜峭度圖無法檢測到所需的頻帶，很難確定軸承是否存在外圈故障，效果較差。

圖9 試驗臺及故障類型和位置

圖10 外圈故障信號的時域波形及其傅里葉譜和包絡譜

5 結論

本文針對快速譜峭度圖中采用固定劃分容易使檢測出的頻帶不能包含完整的故障信息或是提取到無效分量的問題，提出使用頻譜趨勢引導故障共振頻帶探測的方法，構造出基于頻譜趨勢的自適應譜峭度圖。

從頻譜趨勢理論給出了具體分析過程，可以自適應地得到譜峭度圖中每個節點的頻帶劃分，并分別使用本文所提方法與快速譜峭度圖方法對軸承故障仿真信號與試驗信號進行了驗證，表明采用頻譜趨勢劃分頻帶的方法可以精確地對共振頻帶進行劃分。

圖11 使用基于頻譜趨勢的自適應譜峭度圖法處理的結果（二）

圖12 對仿真信號進行快速譜峭度圖計算的結果（二）

[1]谷玉海，朱騰騰，饒文軍，等. 基于EMD二值化圖像和CNN的滾動軸承故障診斷[J]. 振動，測試與診斷，2021，41（1）：105-113.

[2]DWYER R F. Fourth-order spectra of Gaussian amplitude-modulated sinusoids[J]. Journal of the Acoustical Society of America，1991，90（2）：918-926.

[3]Antoni J. Randall R B. The spectral kurtosis：Application to the vibratory surveillance and diagnostics of rotating machines[J]. Mechanical Systems & Signal Processing，2006，20（2）：308-331.

[4]Antoni J. Fast computation of the kurtogram for the detection of transient faults[J]. Mechanical Systems & Signal Processing，2007，21（1）：108-124.

[5]T. Barszcz，A. Jab?oński. A novel method for the optimal band selection for vibration signal demodulation and comparison with the Kurtogram[J]. Mechanical Systems and Signal Processing，2011，25（1）：431-451.

[6]J. Gilles. Empirical Wavelet Transform[J]. IEEE Trans. Signal Process，2013，61（16）：3999-4010.

[7]J. Gilles，K. Heal. A parameterless scale-space approach to find meaningful modes in histograms - Application to image and spectrum segmentation[J]. International Journal of Wavelets，Multiresolution and Information Processing，2014，12（6）：1 450 044.

[8]Feng Z P，Liang M，Chu F L. Recent advances in time-frequency analysis methods for machinery fault diagnosis: A review with application examples[J]. Mechanical Systems and Signal Processing，2013，38（1）：165-205.

[9]Dwyer R. Detection of non-Gaussian signals by frequency domain kurtosis estimation[C]. ICASSP'83，IEEE International Conference on Acoustics, Speech, and Signal Processing，IEEE，1983（8）：607-610.

[10]Antoni N J. The spectral kurtosi: A useful tool for characterizing non-stationary signals[J]. Mechanical Systems & Signal Processing，2006，20（2）：282-307.

[11]Wu Z，Huang NE，Long SR，et al. On the trend, detrending, and variability of nonlinear and nonstationary timeseries[J]. Proceedings of the National Academy of Sciences of the United States of America，2007（104）：14889-14894.

[12]孫永鵬，董增壽. 改進EMD在軸承故障診斷中的應用[J]. 太原科技大學學報，2018，39（1）：25-30.

[13]Torres M E，Colominas M A，Schlotthauer G，et al. A complete ensemble empirical mode decomposition with adaptive noise[C]. 2011 IEEE international conference on acoustics, speech and signal processing (ICASSP)，IEEE，2011：4144-4147.

[14]Bard J. A Review of: “Engineering Optimization: Theory and Practice” [J]. Iie Transactions，1997，29（9）：802-803.

[15]McDonald G L，Zhao Q，Zuo M J. Maximum correlated Kurtosis deconvolution and application on gear tooth chip fault detection[J]. Mechanical Systems and Signal Processing，2012（33）：237-255.

Adaptive Spectral Kurtosis Map and Its Application in Fault Diagnosis of Rolling Bearings

HE Qian，YI Cai，WU Wenyi，LIN Jianhui

(State Key Laboratory of Traction Power, Southwest Jiaotong University, Chengdu 610031, China )

A new adaptive spectral kurtosis map method is proposed in this paper. First, the CEEMDAN methods are applied to decompose signal Fourier spectrum. Then the low-frequency IMF component and residuals are obtained and performed layered incremental calculation. The minimum points of frequency spectrum trend of each layer are found and made as spectrum division boundaries. The divided frequency bands pass through orthogonal spectrum filter and the fault information in the signal is evaluated by the kurtosis index to construct a new adaptive spectrum kurtosis map. This method can overcome the problem that the detected frequency band may not contain complete fault information or may extract invalid components by following fixed division rules. And the effectiveness of the method was verified by using the simulation signal and the test signal of bearing outer ring failure respectively.

fault diagnosis；adaptive spectral kurtosis map；spectral trend；fast spectral kurtosis map

TH133.33；TN911.7

A

10.3969/j.issn.1006-0316.2022.11.002

1006-0316 (2022) 11-0009-10

2022-02-25

國家自然科學基金（51905453）；中國博士后面上基金（2019M663899XB）；中央高校基本科研業務費項目（2682020CX50）；牽引動力國家重點實驗室自主課題（2020TPL-T14）

何倩（1996－），女，四川南充人，碩士研究生，主要研究方向為機械設備故障診斷與信號分析，E-mail：h243052028@163.com。