武杰， 盧振連，， 馬洪儒， 朱艷芳， 吳耀春， 薛曉峰， 姜闊勝

（1. 安陽工學(xué)院 機械工程學(xué)院，河南 安陽 455000；2. 安徽理工大學(xué) 機械工程學(xué)院，安徽 淮南 232000；3. 運城學(xué)院 機電工程系，山西 運城 044000）

0 引言

軸承作為煤礦機械設(shè)備關(guān)鍵運動部件，直接影響設(shè)備的健康狀態(tài)[1-2]。因此，對煤礦機械設(shè)備軸承進行故障診斷具有重要意義。傳統(tǒng)的軸承故障診斷方法主要分為經(jīng)典的時域統(tǒng)計分析方法和以傅里葉理論為核心的全局域變換方法。由于煤礦機械設(shè)備運行過程中振動大、沖擊強，軸承振動信號呈現(xiàn)出瞬態(tài)非平穩(wěn)的特性[3-4]，傳統(tǒng)的軸承故障診斷方法難以識別故障特征。

傳統(tǒng)階次跟蹤方法是通過有鍵相硬件設(shè)備獲取非平穩(wěn)信號的瞬時頻率，再計算出瞬時相位[5-7]，進一步通過重采樣獲取角度域的循環(huán)平穩(wěn)信號，從而實現(xiàn)非平穩(wěn)信號到循環(huán)平穩(wěn)信號的轉(zhuǎn)變。但煤礦機械設(shè)備工作環(huán)境惡劣，不利于有鍵相硬件設(shè)備的安裝，導(dǎo)致瞬時頻率獲取困難，階次跟蹤難以進行，嚴重影響了對軸承的故障診斷。無鍵相階次跟蹤方法可通過振動信號直接獲得瞬時頻率[8-11]，克服了對有鍵相硬件設(shè)備的依賴，能夠很好地實現(xiàn)對軸承故障特征的提取。但該方法運用的前提是通過振動信號獲取準確的瞬時頻率，進而計算出瞬時相位。然而，煤礦機械設(shè)備振動劇烈，轉(zhuǎn)速呈現(xiàn)出強烈的非線性特性，傳統(tǒng)的基于短時傅里葉變換（Short Time Fourier Transform，STFT）的時頻分析方法在估計該工況下的瞬時頻率時魯棒性較差。為解決該問題，文獻[12]提出了一種基于Chirplet 變換的瞬時頻率估計方法，該方法利用Chirplet 變換來逼近信號，從而實現(xiàn)對瞬時頻率的細化，在一定程度上提高了瞬時頻率的估計精度，但由于振動信號的非線性調(diào)制，在速度波動劇烈的情況下，時頻分辨率較低，導(dǎo)致計算的瞬時相位精度不高[13-14]。文獻[15]針對劇烈的非線性調(diào)制，提出了瞬時故障特征頻率的概念，從軸承包絡(luò)信號的時頻譜中提取出瞬時頻率，但在轉(zhuǎn)速急劇變化的情況下，包絡(luò)信號的轉(zhuǎn)頻由于能量較弱，無法準確估計。

本文提出了一種基于諧波匹配補償（Harmonic Matching Compensation，HMC）的時頻分析方法，可實現(xiàn)對軸承劇烈非平穩(wěn)運行工況下瞬時頻率的準確估計；進一步采用無鍵相階次跟蹤方法，實現(xiàn)對軸承故障的準確識別。

1 基于HMC 和無鍵相階次跟蹤的軸承故障診斷原理

基于HMC 和無鍵相階次跟蹤的軸承故障診斷流程如圖1 所示。首先，采用基于HMC 的時頻分析方法對振動信號進行處理，準確估計出瞬時頻率。其次，由于在瞬態(tài)非平穩(wěn)工況下，信號中混疊的各階諧波分量隨轉(zhuǎn)速一起波動，通過Vold-Kalman 濾波方法提取諧波分量信號。再次，計算出瞬時相位，進而實現(xiàn)信號的角度域重采樣。最后，對重采樣的信號進行快速傅里葉變換，從而進行階次譜分析，識別出軸承故障特征。

圖1 基于HMC 和無鍵相階次跟蹤的軸承故障診斷流程Fig. 1 Bearing fault diagnosis process based on harmonic matching compensation and keyless phase order tracking

1.1 基于HMC 的瞬時頻率估計

在瞬態(tài)非平穩(wěn)工況下，對于一個有限的能量信號x(t)，其HMC 被定義為

式中：t，f分別為時間和頻率；S(t)為瞬時頻率的解調(diào)函數(shù)；hT(t)為 窗函數(shù)，其中心為τ、寬度為T，為了便于區(qū)分信號中的諧波分量，選擇具有低混淆特性的漢寧窗。

對式（1）進一步變形，可得

在HMC 計算過程中，對于有限能量信號x(t)，先通過STFT 對瞬時頻率進行粗略估計，之后采用泰勒公式將粗略估計的瞬時頻率f(t)展開：

式中：n為階數(shù)；Rn(t)為泰勒余項。

在式（3）中，余項的大小取決于余項中t的階次。當t處于高階次時，余項值不斷減小，同時無限逼近瞬時頻率f(t)。因此，可以用關(guān)于t的負指數(shù)冪的形式來近似余項，同時通過調(diào)節(jié)參數(shù)c實現(xiàn)對瞬時頻率的匹配：

式中：f0為冪指數(shù)函數(shù)的常數(shù)項；d為指數(shù)因子。

由式（4）可知，粗略估計出的瞬時頻率曲線為指數(shù)冪形式。因此，瞬時頻率的解調(diào)函數(shù)S(t)可以被定義為

則式（2）中解調(diào)操作因子 e xp(?j2πS(t))可表示為

為了方便計算，d可設(shè)定為0.5。對于和c，可通過將瞬時頻率曲線上任意2 個時頻點 (t1，f1)和(t2，f2)代入式（5）進行求解。因此，可在每個漢寧窗中，分別獲得相應(yīng)的解調(diào)函數(shù)S(t)和解調(diào)操作因子exp(?j2πS(t))。 將 解 調(diào) 函 數(shù)S(t)和 解 調(diào) 操 作 因 子exp(?j2πS(t))分別轉(zhuǎn)換成向量形式和。因此，經(jīng)過HMC 變換解調(diào)后的信號為

式中 real[·]為實數(shù)運算操作。

1.2 基于Vold-Kalman 濾波的諧波分量提取

Vold-Kalman 濾波器是一種可有效跟蹤信號波動變化的時變?yōu)V波器，能跟隨信號特性自適應(yīng)分離和提取諧波[16-17]。相對于傳統(tǒng)的帶通濾波器，Vold-Kalman 濾波器的中心帶寬可以跟隨轉(zhuǎn)速的波動自適應(yīng)調(diào)整，這種性質(zhì)使其很適用于對煤礦機械設(shè)備軸承信號的跟蹤濾波。在提取第k次諧波分量的同時，為有效抑制其他相鄰諧波分量，Vold-Kalman 濾波器的瞬時帶寬Ck需滿足如下條件：

式中Ok為第k次諧波和其相鄰諧波之間的距離。

1.3 瞬時相位計算

當?shù)趉次諧波分量e(n)的時域信號從原始振動信號中被提取出來后，采用Hilbert 變換計算第k次諧波的瞬時相位：

式中： unwrap[·]為 相位解卷積操作，可將e(n)的瞬時相 位 從 區(qū) 間 [?π，π]映 射 到 相 位 累 積 空 間；(n)為e(n)經(jīng)Hilbert 變換后的結(jié)果。

理論上，第k次諧波信號的瞬時相位 φk(n)與原始信號瞬時相位 φ(n)之間嚴格遵循線性比例關(guān)系。但信號在實際的傳遞過程中，受振動傳遞路徑的影響，瞬時相位會發(fā)生偏移。因此，式（9）可修正為

式中：a為 φk(n) 和 φ(n)之 間的線性比例系數(shù)； φs(n)為由振動傳遞路徑引入的相位偏移。

由于振動傳遞路徑引入的相位偏移較大，需要對相位進行修正。修正后的實際瞬時相位為

1.4 信號角度域重采樣

在角度域重采樣過程中，為了避免產(chǎn)生頻率混淆， ?θ需要滿足如下約束：

式中fs為信號采樣頻率。

通過時間域和角度域的映射關(guān)系，以 ?θ為采樣間隔，可實現(xiàn)原始時間域信號在角度域的等角度重采樣，從而獲得角度域重采樣信號。

2 仿真驗證

根據(jù)煤礦機械設(shè)備軸承的實際瞬態(tài)非平穩(wěn)運行工況，設(shè)置軸承振動仿真信號[18-19]：

假設(shè)軸承外圈固定在軸承座中，內(nèi)圈隨轉(zhuǎn)軸一起旋轉(zhuǎn)，且在軸承外圈上存在1 處局部剝落缺陷。因此，軸承的故障脈沖信號s(t)可設(shè)置為[19]

假設(shè)軸承故障信號中包含2 個關(guān)于轉(zhuǎn)頻的諧波分量，其幅值分別為B1=0.005g（g為重力加速度），B2=0.01g， 初始相位分別為 φ1=π/6， φ2=?π/3。在測試過程中軸承經(jīng)歷了劇烈的轉(zhuǎn)速上升和轉(zhuǎn)速下降過程，因此瞬時頻率[18]為

將軸承故障特征頻率對外圈轉(zhuǎn)頻歸一化處理[20]，設(shè)軸承的外圈故障特征頻率為轉(zhuǎn)軸轉(zhuǎn)頻的4.35 倍，則故障脈沖的平均角度間隔為360°/4.35≈82.76°。設(shè)置信號的采樣頻率為20 kHz，采樣持續(xù)時間為2 s，噪聲干擾為?3 dB 的白噪聲。

軸承瞬時轉(zhuǎn)速模擬曲線如圖2 所示，可模擬煤礦機械設(shè)備作業(yè)過程中轉(zhuǎn)速劇烈波動工況。

圖2 軸承瞬時轉(zhuǎn)速模擬曲線Fig. 2 Simulation curve of bearing instantaneous speed

仿真信號的時域波形如圖3 所示。由于轉(zhuǎn)速劇烈波動，各諧波分量及噪聲干擾混淆，難以發(fā)現(xiàn)由軸承缺陷所激發(fā)的故障沖擊特征。此外，受軸承轉(zhuǎn)速劇烈波動影響，傳統(tǒng)的時域和頻域統(tǒng)計指標分析不再適用。

圖3 仿真信號時域波形Fig. 3 Time domain waveform of simulated signal

采用本文方法對仿真信號進行處理。首先，通過快速譜峭度分析法[21]獲得仿真信號的譜峭度，如圖4 所示，進而識別出軸承的最佳共振頻帶。然后，根據(jù)最佳共振頻帶設(shè)計中心頻率為1 875 Hz、帶寬為416.67 Hz 的帶通濾波器，提取出仿真信號的共振解調(diào)序列，如圖5（a）所示。最后，采用Hilbert 變換得到仿真信號的包絡(luò)，如圖5（b）所示。

圖4 仿真信號譜峭度Fig. 4 Spectral kurtosis of simulated signal

圖5 仿真信號的共振解調(diào)序列和包絡(luò)Fig. 5 Resonance demodulation sequence and envelope of simulated signal

為驗證基于HMC 的時頻分析方法在轉(zhuǎn)速劇烈波動工況下對瞬時頻率估計的有效性，分別采用HMC 方法和STFT 方法對瞬時頻率進行估計，得到的時頻譜如圖6 所示。可看出采用HMC 方法得到的時頻譜中時頻聚集性顯著且分辨率高，而采用STFT 方法得到的時頻譜在劇烈轉(zhuǎn)速波動和噪聲雙重干擾下較模糊，且時頻聚集性和分辨率較差。

圖6（a）中有2 個能量突出的諧波分量，對諧波分量Ⅰ進行瞬時頻率估計，結(jié)果如圖7 所示?？煽闯龉烙嫷乃矔r頻率曲線與實際的瞬時頻率曲線（通過式（16）計算所得）高度重合，瞬時頻率估計值與實際值之間的最大相對誤差不超過1%。

圖6 不同方法下仿真信號時頻譜Fig. 6 Time-frequency spectrum of simulated signal under different methods

圖7 仿真信號瞬時頻率估計結(jié)果Fig. 7 Instantaneous frequency estimation result of simulated signal

以HMC 方法精確估計出的瞬時頻率為中心頻率，利用Vold-Kalman 濾波器對諧波分量Ⅱ進行自適應(yīng)提取，結(jié)果如圖8 所示。

圖8 Vold-Kalman 濾波器提取的諧波分量Fig. 8 Harmonic components extracted by Vold-Kalman filter

從圖8 可看，出提取出的諧波分量波形平滑，幅值為0.01g，與設(shè)置的仿真信號幅值B2=0.01g相同。根據(jù)式（8），取Vold-Kalman 濾波器的瞬時帶寬為 0.5階次。

依據(jù)式（9）和式（11），估計出瞬時相位，并與實際瞬時相位對比，結(jié)果如圖9 所示?？煽闯龉烙嫷乃矔r相位曲線與實際瞬時相位曲線高度重合，瞬時相位估計值與實際值的最大相對誤差不超過2%。

圖9 仿真信號瞬時相位估計結(jié)果Fig. 9 Instantaneous phase estimation result of simulated signal

根據(jù)估計出的瞬時相位，通過式（12）獲得時間域和角度域的映射關(guān)系，進而完成角度域的重采樣。再將重采樣的角度域信號進行快速傅里葉變換，得到包絡(luò)階次譜，如圖10 所示?？煽闯霰疚姆椒ê陀墟I相階次跟蹤方法得到的包絡(luò)階次譜高度吻合，通過本文方法可清晰找到軸承外圈故障特征階次（4.35 階次）及相關(guān)倍頻特征階次，證明軸承外圈存在故障。

圖10 仿真信號包絡(luò)階次譜Fig. 10 Envelope order spectrum of simulated signal

3 試驗驗證

為驗證本文方法對煤礦機械設(shè)備軸承故障診斷的有效性，搭建軸承模擬試驗系統(tǒng)，如圖11 所示。

圖11 軸承模擬試驗系統(tǒng)Fig. 11 Bearing simulation test system

將軸承外圈通過過盈配合固定于軸承座內(nèi)，軸承內(nèi)圈與軸一起旋轉(zhuǎn)。電動機與伺服驅(qū)動系統(tǒng)提供動力輸入，驅(qū)動軸承轉(zhuǎn)動。采用3 個靈敏度為50 mV/g的單軸振動加速度傳感器采集軸承振動信號。轉(zhuǎn)速計為光電脈沖編碼器，用來采集軸承轉(zhuǎn)頻信息。信號采集系統(tǒng)主要設(shè)備為Coco?80 數(shù)字采集儀。為了模擬煤礦機械設(shè)備運行時軸承運行工況，設(shè)置軸承經(jīng)歷瞬間啟動、載荷及轉(zhuǎn)速迅速上升和瞬間停止3 個過程。設(shè)置采樣頻率為25.6 kHz。

軸承幾何參數(shù)見表1。試驗設(shè)置軸承外圈存在人為設(shè)置的磨損故障缺陷，如圖12 所示。

表1 軸承幾何參數(shù)Table 1 Bearing geometry parameters

圖12 軸承外圈故障Fig. 12 Outer ring fault of bearing

通過歸一化操作[20]可得軸承外圈故障特征階次：

式中：z為滾動體個數(shù)；fr為外圈轉(zhuǎn)頻；r為節(jié)圓半徑；D為滾動體直徑； ?為接觸角。

將表1 數(shù)據(jù)代入式（17），計算出軸承外圈故障特征階次為3.072 7。

采集的原始振動信號、轉(zhuǎn)速和截取的試驗信號如圖13 所示。由于轉(zhuǎn)速的劇烈瞬態(tài)波動，從圖13（a）和圖13（c）中無法識別出任何關(guān)于軸承故障的沖擊特征。

圖13 試驗采集信號Fig. 13 Test acquisition signal

通過快速譜峭度分析法獲得截取的試驗信號的譜峭度，如圖14 所示，進而可識別出軸承的最佳共振頻帶。根據(jù)獲得的最佳共振頻帶設(shè)計中心頻率為1 718.8 Hz、帶寬為312.5 Hz 的帶通濾波器，進而提取出截取的試驗信號的共振解調(diào)序列，并采用Hilbert 變換得到截取的試驗信號的包絡(luò)，如圖15所示。

圖14 試驗信號譜峭度Fig. 14 Spectral kurtosis of test signal

圖15 試驗信號的共振解調(diào)序列和包絡(luò)Fig. 15 Resonance demodulation sequence and envelope of test signal

分別采用HMC 方法和STFT 方法處理試驗信號，得到時頻譜，如圖16 所示。可看出采用STFT 方法處理后的時頻譜中識別不出任何頻率波動的趨勢，而由于HMC 方法在瞬態(tài)條件下的優(yōu)異性能，可清楚識別出瞬時頻率的變化趨勢。

圖16（a）中有2 個諧波能量較為顯著。為了確定這2 個諧波分量的位置，在0.058 38 s 處對其進行檢測，其中心頻率分別為 975.3 Hz 和953.6 Hz，表明2 個諧波分量之間的頻率間隔為21.7 Hz。因此，確定2 個諧波分量分別為試驗信號的第45（975.3 Hz/21.7 Hz≈45）個 和 第44（953.6 Hz/21.7 Hz≈44）個 分量。本文通過第45 個諧波分量估計試驗信號的瞬時頻率，結(jié)果如圖17 所示。可看出瞬時頻率估計曲線與實際的瞬時頻率曲線高度重合，瞬時頻率估計值與實際值之間的最大相對誤差不超過1%。

圖16 不同方法下試驗信號時頻譜Fig. 16 Time-frequency spectrum of test signal under different methods

圖17 試驗信號瞬時頻率估計結(jié)果Fig. 17 Instantaneous frequency estimation result of test signal

采用HMC 方法準確估計出瞬時頻率后，通過Vold-Kalman 濾波器對第45 個諧波分量進行自適應(yīng)提取。根據(jù)式（9）和式（10），估計出瞬時相位，并與實際瞬時相位對比，結(jié)果如圖18 所示。可看出估計的瞬時相位曲線與實際的瞬時相位曲線高度重合，瞬時相位估計值與實際值之間的最大相對誤差不超過2%。

圖18 試驗信號瞬時相位估計結(jié)果Fig. 18 Instantaneous phase estimation result of test signal

根據(jù)估計出的瞬時相位，通過式（12）獲得時間域和角度域的映射關(guān)系，進而完成角度域重采樣。對重采樣的包絡(luò)信號進行快速傅里葉變換，獲得包絡(luò)階次譜，如圖19 所示?？煽闯霰疚姆椒ê陀墟I相階次跟蹤方法得到的包絡(luò)階次譜高度吻合，通過無鍵相階次跟蹤方法可清晰找到軸承外圈故障特征階次（3.072 7 階次）及相關(guān)倍頻特征階次，由此可明確推斷出軸承外圈存在故障。

圖19 試驗信號包絡(luò)階次譜Fig. 19 Envelope order spectrum of test signal

為進一步驗證本文方法的有效性，對軸承進行內(nèi)圈故障加工，如圖20 所示。應(yīng)用本文方法得到軸承內(nèi)圈故障振動信號的包絡(luò)階次譜，如圖21 所示。可看出本文方法和有鍵相階次跟蹤方法所獲得的故障特征階次相同，但本文方法獲得的故障特征階次幅值明顯高于有鍵相階次跟蹤方法獲得的故障特征階次幅值，更有利于識別故障特征。

圖20 軸承內(nèi)圈故障Fig. 20 Inner ring fault of bearing

圖21 軸承內(nèi)圈故障振動信號的包絡(luò)階次譜Fig. 21 Envelope order spectrum of vibration signal of bearing inner ring fault

4 結(jié)論

1） 基于HMC 的時頻分析方法不依賴傳統(tǒng)的局部平穩(wěn)假設(shè)要求，可自適應(yīng)地匹配調(diào)整瞬時頻率信號，并對瞬時頻率信號進行自適應(yīng)補償，進而降低了信號的非線性度，避免了傳統(tǒng)基于STFT 的時頻分析方法在瞬態(tài)非平穩(wěn)運行工況下時頻譜模糊的問題，從而精確估計出瞬時頻率。

2） 將基于HMC 的時頻分析方法與無鍵相階次 跟蹤相結(jié)合，應(yīng)用于軸承故障診斷，能夠準確估計瞬時頻率，獲得較為精確的瞬時相位，通過包絡(luò)階次譜可準確且明顯地表征出軸承故障特征階次，從而有效識別故障特征。