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

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

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

（西南交通大學(xué) 牽引動力國家重點(diǎn)實(shí)驗(yàn)室，四川 成都 610031）

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

故障診斷；自適應(yīng)譜峭度圖；頻譜趨勢；快速譜峭度圖

滾動軸承廣泛應(yīng)用于旋轉(zhuǎn)機(jī)械中，在現(xiàn)代機(jī)械設(shè)備中起著非常重要的作用[1]。在正常情況下，軸承在高速和重載下運(yùn)行，即使輕微損壞對設(shè)備也非常危險(xiǎn)，如果不能立即檢測到損壞，軸承將不穩(wěn)定運(yùn)行，嚴(yán)重時(shí)甚至損壞設(shè)備。因此，有必要實(shí)時(shí)監(jiān)測軸承的運(yùn)行狀態(tài)并進(jìn)行故障診斷。

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

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

針對上述情況，本文提出基于頻譜趨勢的自適應(yīng)劃分頻帶的方法，并將其與峭度指標(biāo)結(jié)合，構(gòu)造出自適應(yīng)譜峭度圖算法，并通過仿真信號及實(shí)測信號證明了該方法的有效性。

1 基本理論

1.1 譜峭度

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

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

1.2 快速譜峭度圖

Antoni和Randall[2]提出了快速譜峭度圖來表示頻譜分割平面上信號的峭度值，它可以被認(rèn)為是中心頻率和帶寬的不同尺度下峭度值的函數(shù)?？焖僮V峭度的“快”體現(xiàn)在其快速頻帶分割技術(shù)中，該技術(shù)使用二元1/3樹結(jié)構(gòu)將頻帶細(xì)分為多個(gè)合理的寬度。其實(shí)現(xiàn)過程如下：

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

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

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

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

圖1 快速譜峭度圖示例結(jié)果

2 自適應(yīng)譜峭度圖

2.1 頻譜趨勢自適應(yīng)劃分頻帶

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

基于上述分析，本文采用大程度減少模態(tài)混疊和不必要迭代篩選且數(shù)據(jù)分解更完整的自適應(yīng)噪聲的集合經(jīng)驗(yàn)?zāi)B(tài)分解（CEEMDAN）來完成頻譜趨勢的尋找。

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

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

（2）應(yīng)用CEEMDAN得到頻譜的IMF分量以及殘差項(xiàng)；

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

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

圖2顯示了原始信號的時(shí)域波形和頻譜，對頻譜進(jìn)行CEEMDAN分解的結(jié)果如圖3所示，圖4展示了部分頻帶劃分結(jié)果，其中藍(lán)色曲線表示殘余項(xiàng)與不同個(gè)數(shù)IMF疊加的結(jié)果、紅色曲線表示為頻譜邊界劃分的位置。

2.2 自適應(yīng)譜峭度圖

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

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

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

（2）計(jì)算相對應(yīng)信號的傅里葉頻譜，并對頻譜進(jìn)行自適應(yīng)噪聲的集合經(jīng)驗(yàn)?zāi)B(tài)分解，得到各IMF分量及殘余項(xiàng)；

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

圖2 振動信號時(shí)域波形及頻譜

圖3 CEEMDAN分解結(jié)果

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

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

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

（7）選取峭度值最大的分量，進(jìn)行Hilbert解調(diào)，進(jìn)行故障分析，當(dāng)包絡(luò)譜中包含有可識別的故障特征頻率及其倍頻時(shí)，則可認(rèn)為該軸承存在相應(yīng)的局部故障。

圖4 頻譜劃分結(jié)果

3 仿真信號分析

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

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

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

圖5 自適應(yīng)譜峭度圖算法流程圖

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

表1 軸承故障仿真信號參數(shù)

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

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

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

4 試驗(yàn)分析

為進(jìn)一步驗(yàn)證該基于頻譜趨勢的自適應(yīng)譜峭度圖的有效性，本節(jié)采用試驗(yàn)數(shù)據(jù)對該方法進(jìn)行驗(yàn)證。

圖6 仿真信號時(shí)域圖、頻譜圖、包絡(luò)譜

圖7 采用基于頻譜趨勢的自適應(yīng)譜峭度圖法處理的結(jié)果（一）

圖8 對仿真信號進(jìn)行快速譜峭度圖計(jì)算的結(jié)果（二）

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

表2 雙列圓錐滾子軸承的基本參數(shù)

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

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

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

圖9 試驗(yàn)臺及故障類型和位置

圖10 外圈故障信號的時(shí)域波形及其傅里葉譜和包絡(luò)譜

5 結(jié)論

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

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

圖11 使用基于頻譜趨勢的自適應(yīng)譜峭度圖法處理的結(jié)果（二）

圖12 對仿真信號進(jìn)行快速譜峭度圖計(jì)算的結(jié)果（二）

[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]孫永鵬，董增壽. 改進(jìn)EMD在軸承故障診斷中的應(yīng)用[J]. 太原科技大學(xué)學(xué)報(bào)，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

國家自然科學(xué)基金（51905453）；中國博士后面上基金（2019M663899XB）；中央高?；究蒲袠I(yè)務(wù)費(fèi)項(xiàng)目（2682020CX50）；牽引動力國家重點(diǎn)實(shí)驗(yàn)室自主課題（2020TPL-T14）

何倩（1996－），女，四川南充人，碩士研究生，主要研究方向?yàn)闄C(jī)械設(shè)備故障診斷與信號分析，E-mail：h243052028@163.com。