榮 鋒， 陳 寧， 郭翠娟， 楊明珠

(1. 天津工業(yè)大學(xué) 電子與信息工程學(xué)院，天津 300387； 2. 天津市光電檢測技術(shù)與系統(tǒng)重點(diǎn)實(shí)驗(yàn)室，天津 300387)

旋轉(zhuǎn)機(jī)械廣泛應(yīng)用于社會生產(chǎn)的各個(gè)行業(yè)，保障其安全穩(wěn)定的運(yùn)行有著重大意義。目前，對旋轉(zhuǎn)機(jī)械的狀態(tài)監(jiān)測與故障診斷的重要手段是對振動的監(jiān)測與分析[1]。旋轉(zhuǎn)機(jī)械在運(yùn)行狀態(tài)下的轉(zhuǎn)速嚴(yán)格來說是波動的，尤其是升降速過程，此時(shí)的振動信號屬于非平穩(wěn)信號，不適合于常規(guī)的傅里葉頻譜分析[2]。階比分析是工程實(shí)際中常用的旋轉(zhuǎn)機(jī)械變轉(zhuǎn)速狀態(tài)下振動信號的分析方法，其關(guān)鍵在于將時(shí)域非平穩(wěn)信號通過等角度重采樣轉(zhuǎn)化為角度域的平穩(wěn)信號，這個(gè)過程稱為階比跟蹤[3-4]。階比跟蹤具體來說是指準(zhǔn)確獲得參考軸的基準(zhǔn)頻率以及得到等角度重采樣時(shí)刻的過程。傳統(tǒng)的階比跟蹤方法有硬件階比跟蹤和計(jì)算階比跟蹤，兩者都需要有特定的測量轉(zhuǎn)速的硬件裝置，在不便安裝的場合，兩種方法都無能為力。近年來，隨著學(xué)者不斷的研究，出現(xiàn)了無需硬件測速裝置的軟件階比跟蹤技術(shù)，其關(guān)鍵在于從振動信號中提取出參考軸的轉(zhuǎn)速曲線，從而計(jì)算得到重采樣時(shí)刻。該方法極大削弱了階比分析對硬件的依賴，得到廣泛的應(yīng)用。

由振動信號獲取轉(zhuǎn)速曲線的方法有很多種，郭瑜等首先提出利用短時(shí)傅里葉變換(STFT)或Wigner-Ville分布(WVD)的時(shí)頻分析階比跟蹤方法，通過對振動信號的時(shí)頻分布進(jìn)行峰值搜索，從而獲得參考軸的轉(zhuǎn)速曲線。曹書峰等[5]提出了基于時(shí)頻融合分布的轉(zhuǎn)速估計(jì)方法，通過將兩種時(shí)頻分析的時(shí)頻分布進(jìn)行融合，之后再進(jìn)行峰值搜索和最小二乘擬合得到參考軸的轉(zhuǎn)速曲線。李輝等[6]利用Hilbert-Huang變換(HHT)獲得轉(zhuǎn)速曲線。以上文獻(xiàn)雖然都有一定的成果，但是也各自存在一些問題，比如時(shí)頻分析存在時(shí)頻聚焦性差或交叉項(xiàng)干擾等問題[7]，而且在數(shù)據(jù)量較大時(shí)，產(chǎn)生的時(shí)頻矩陣會占用極大的計(jì)算機(jī)硬件內(nèi)存，大大降低計(jì)算效率；HHT變換過程中的經(jīng)驗(yàn)?zāi)Ｊ椒纸?Empirical Mode Decomposition，EMD)會使相鄰的本征模態(tài)函數(shù)(Intrinsic Mode Function，IMF)之間產(chǎn)生模態(tài)混疊，且Hilbert變換會在信號兩端產(chǎn)生嚴(yán)重震蕩，甚至可能出現(xiàn)無意義的負(fù)頻率，影響轉(zhuǎn)速曲線的估計(jì)精度[8-9]。

為此，本文提出一種基于EEMD的HHT與時(shí)頻重排算法相結(jié)合的轉(zhuǎn)速曲線估計(jì)方法。首先由EEMD得到參考轉(zhuǎn)軸所對應(yīng)的IMF，之后分兩步得到該IMF的瞬時(shí)頻率：先對IMF整體做HHT變換，得到其各個(gè)點(diǎn)的瞬時(shí)頻率值；再截取IMF兩端各一定長度的數(shù)據(jù)，分別進(jìn)行時(shí)頻重排變換，在得到的兩個(gè)時(shí)頻重排矩陣中進(jìn)行限制搜索條件的瘠提取，并對瘠線擬合，得到IMF兩端截取點(diǎn)的瞬時(shí)頻率。之后用第二步得到的瞬時(shí)頻率值代替第一步得到的對應(yīng)點(diǎn)處的瞬時(shí)頻率值，并由最小二乘法將“拼接”頻率擬合成一條光滑的曲線，即為參考軸轉(zhuǎn)頻曲線，利用轉(zhuǎn)頻與轉(zhuǎn)速間的關(guān)系即可得到轉(zhuǎn)速曲線。通過對實(shí)測信號分析來詳細(xì)說明算法實(shí)施步驟，分析結(jié)果表明本文方法行之有效。

1 EEMD-HHT和時(shí)頻重排算法原理

1.1 EEMD-HHT算法

EEMD-HHT是對HHT的一種改進(jìn)，其利用EEMD分解方法，有效的抑制了HHT變換中產(chǎn)生的模態(tài)混疊現(xiàn)象[10-11]，使得分解得到的IMF更加純凈。EEMD利用白噪聲在整個(gè)頻域內(nèi)分布均勻的特性，在原始信號中加入白噪聲，使得信號在不同尺度上連續(xù)。其在分解時(shí)要設(shè)定執(zhí)行EMD的次數(shù)M，及疊加的隨機(jī)白噪聲序列的幅值標(biāo)準(zhǔn)差σ。EEMD-HHT的主要思想是把一個(gè)時(shí)間序列信號由EEMD分解成一系列不同尺度的單分量IMF，之后對每個(gè)IMF進(jìn)行Hilbert變換，得到各IMF所對應(yīng)的瞬時(shí)頻率值和幅值，從而得到振幅-頻率-時(shí)間的三維譜分布。EEMD-HHT有效的抑制了HHT的模態(tài)混疊現(xiàn)象，但是分解過程仍會產(chǎn)生虛假分量成分，使得有效本征模函數(shù)不好識別，而且Hilbert變換在信號端點(diǎn)處造成的頻率震蕩依舊存在。

1.2 時(shí)頻重排算法

STFT變換和WVD變換是兩種應(yīng)用最為廣泛的時(shí)頻分析方法。但是前者無法兼顧時(shí)間分辨率和頻率分辨率，時(shí)頻聚焦性相對較差，后者分析多成分信號分量時(shí)，會產(chǎn)生交叉項(xiàng)干擾。時(shí)頻重排算法最初是為了改進(jìn)譜圖的可讀性[12]。信號x(t)的譜圖定義為其STFT變化的幅值平方[13]，表達(dá)式為

(1)

式中：h(t)為STFT所用的窗函數(shù)。譜圖表達(dá)式可以寫成信號的WVD和分析窗的WVD的二維卷積形式

(2)

此分布一定程度上減弱了信號WVD產(chǎn)生的相關(guān)項(xiàng)，但是降低了時(shí)頻分辨率，是以邊緣性質(zhì)和一階矩有偏為代價(jià)的。上式表明，WVDh(t-s,f-a)是在點(diǎn)(t,f)附近設(shè)定了一個(gè)鄰域來分配信號WVD的加權(quán)平均值。而重排算法就是將時(shí)頻分布面上任意一點(diǎn)(t,f)處的譜值從原來位置移動到該點(diǎn)附近信號能量的重心

(3)

(4)

重排后譜圖上任一點(diǎn)(t′,f′)的值是所有重排到這一點(diǎn)的原譜圖值的和，即

(5)

式中：δ(t)為沖擊函數(shù)。重排矩陣不僅利用了短時(shí)傅里葉變換的幅值信息，還利用了其相位信息，提高了時(shí)頻聚焦性，改進(jìn)了短時(shí)傅里葉變換所得到譜圖的可讀性，并且可以直接由重排譜提取瘠線，求得信號的瞬時(shí)參數(shù)。但是當(dāng)數(shù)據(jù)量很大時(shí)，時(shí)頻矩陣會占據(jù)大量內(nèi)存空間，使得計(jì)算效率嚴(yán)重降低，計(jì)算時(shí)間大大增加。

2 EEMD-HHT與時(shí)頻重排結(jié)合的估計(jì)方法

針對以上問題，本文提出了將EEMD-HHT與時(shí)頻重排相結(jié)合的方法來估計(jì)參考軸的轉(zhuǎn)速曲線，轉(zhuǎn)頻曲線估計(jì)算法實(shí)現(xiàn)的流程圖如圖1所示。其詳細(xì)步驟如下：

(1) 將振動信號通過EEMD分解為一系列的IMF。

(2) 提取參考軸所對應(yīng)的IMF分量。由于EEMD分解會產(chǎn)生虛假分量，因此本文提出采用相關(guān)系數(shù)與能量譜結(jié)合的方法，準(zhǔn)確提取所需IMF。具體算法如下：首先，將所有IMF和原信號進(jìn)行歸一化處理，以防止對一些幅值較小的真實(shí)IMF做出誤判。再分別計(jì)算各IMF分量與原振動信號間的相關(guān)系數(shù)值，相關(guān)系數(shù)公式為

(6)

式中：X和Y分別表示某個(gè)IMF分量和原始信號。相關(guān)系數(shù)越大代表該分量與原信號相關(guān)性越好，越能表示真實(shí)的信號成分，反之代表與原信號相關(guān)性不大，可能為虛假分量[14]。

圖1 算法流程圖Fig.1 Flow chart of algorithm

其次，定義某IMF分量中各點(diǎn)的平方和作為其能量的標(biāo)志。分別計(jì)算各IMF的能量值，并做歸一化處理，以各IMF所占的百分比作為判斷指標(biāo)，公式如下

(7)

式中：X表示某個(gè)IMF分量，M為IMF的個(gè)數(shù)。ek越大表示第K個(gè)IMF分量所占比重越大，越小則表示所占比重越小。由于振動信號是由參考軸轉(zhuǎn)動提供動力產(chǎn)生的，所以其對應(yīng)的IMF相關(guān)系數(shù)與能量占比兩個(gè)數(shù)據(jù)都可能更大。綜合考慮兩種計(jì)算結(jié)果，選取相關(guān)系數(shù)值和能量占比都較大的IMF分量作為參考軸的本征模態(tài)函數(shù)。

(3) 對選擇的IMF分量作Hilbert變換，得到其各個(gè)點(diǎn)的瞬時(shí)頻率值。

(4) 截取IMF前后兩端各一定長度數(shù)據(jù)做時(shí)頻重排變換，得到兩個(gè)時(shí)頻重排矩陣。截取的數(shù)據(jù)長度可根據(jù)實(shí)際情況具體確定，本文提供如下一種方法：觀察第3步得到的IMF時(shí)頻曲線，將兩端突變嚴(yán)重的數(shù)據(jù)段作為截取數(shù)據(jù)段，截取長度最好為2的N次方，以方便計(jì)算機(jī)計(jì)算，但不要過大，以免占用過多內(nèi)存，影響計(jì)算效率。重排矩陣的列代表時(shí)間，行代表頻率，其交叉點(diǎn)為該點(diǎn)處的能量值[15-16]。

(5) 分別對兩個(gè)時(shí)頻矩陣進(jìn)行瘠線提取。雖然時(shí)頻重排算法提高了短時(shí)傅里葉變換的聚焦性，但是由于IMF兩端點(diǎn)處的震蕩，重排矩陣中會產(chǎn)生明顯偏離轉(zhuǎn)軸頻率的值。所以本文提出設(shè)置頻率搜索范圍進(jìn)行瘠線提取，具體步驟為：

首先，設(shè)定能量閾值Ea，式中Esum為重排矩陣的總能量，N為IMF的長度。

(8)

然后，以增速過程為例，對于前端數(shù)據(jù)，從(ti,fi)中找出滿足fd

(9)

(6) 用第5步得到IMF前后兩端截取各點(diǎn)的瞬時(shí)頻率值代替第3步Hilbert變換得到的對應(yīng)點(diǎn)處的值，組成一個(gè)新的頻率數(shù)組。

(7) 對新的頻率數(shù)組重新進(jìn)行最小二乘擬合，獲得光滑的轉(zhuǎn)速曲線。

3 應(yīng)用實(shí)例

將本文所提方法應(yīng)用于某雙跨轉(zhuǎn)子平臺的轉(zhuǎn)子不平衡故障分析，以實(shí)際信號分析來說明本文算法的分析過程，并對其有效性和優(yōu)越性進(jìn)行驗(yàn)證。試驗(yàn)中先為質(zhì)量均勻轉(zhuǎn)子添加適當(dāng)配重，然后接通電機(jī)電源，將轉(zhuǎn)速調(diào)至1 500 r/min，待轉(zhuǎn)軸速度穩(wěn)定時(shí)，設(shè)定增速至2 000 r/min，并開始采集加速度傳感器獲得的振動信號。采樣速率為2 000 Hz，采樣時(shí)間為5 s，數(shù)據(jù)長10 000點(diǎn)。圖2為升速階段轉(zhuǎn)子不平衡狀態(tài)下振動信號的時(shí)序波形圖，此信號為典型的非平穩(wěn)信號，對其直接進(jìn)行傅里葉變換會出現(xiàn)“頻率模糊”的現(xiàn)象，掩蓋了部分故障信息，如圖3所示。

圖2 振動信號時(shí)域波形Fig.2 Time domain waveform of vibration signal

圖3 振動信號頻譜Fig.3 Spectrum of vibration signal

對振動信號直接進(jìn)行EMD分解，得到的IMF如圖4所示。從圖中可以看到，相鄰IMF間產(chǎn)生了模態(tài)混疊現(xiàn)象，使IMF失去了真實(shí)意義。

因此，用本文所提方法對參考軸的轉(zhuǎn)速曲線進(jìn)行估計(jì)，之后用于階比分析，獲得階比譜，檢驗(yàn)方法的有效性。首先，設(shè)置EEMD中隨機(jī)白噪聲的標(biāo)準(zhǔn)差σ=0.5，執(zhí)行次數(shù)M=200，得到不包含殘差在內(nèi)的共13個(gè)IMF，其中IMF1為原始信號，IMF2～13為分解得到的IMF分量，RES為殘差，如圖5所示?？梢?，EEMD分解結(jié)果中出現(xiàn)了很多虛假分量，需要對其進(jìn)行篩選。

圖4 EMD分解得到的IMFFig.4 IMF obtained by EMD decomposition

圖5 EEMD分解得到的IMFFig.5 IMF obtained by EEMD decomposition

分別求各IMF與原始信號的相關(guān)系數(shù)μi(i=1,2,3,…,13)，如表1所示。由表中數(shù)據(jù)得，除作為原始信號的IMF1外，IMF4、IMF5和IMF6與原始信號的相關(guān)性相對較大。

表1 各IMF與原信號間的相關(guān)系數(shù)值

計(jì)算各IMF的能量值，并做歸一化處理，得到如圖6所示的能量譜柱狀圖。由于IMF1為原始信號分量，所以忽略其影響。從圖中可以看到，除IMF1外，IMF4、IMF5和IMF6所占比重較大。

圖6 能量譜圖Fig.6 Energy spectrum

參考文獻(xiàn)[14]，設(shè)置閾值ε=max(μi)/10，將大于ε的μi所對應(yīng)的IMF保留下來，認(rèn)為該IMF與原始信號相關(guān)性較好，是真實(shí)分量；將小于該閾值的IMF合并到殘差中，認(rèn)為其與原始信號相關(guān)性較差，是虛假分量。經(jīng)過篩選得到的新的IMF，如圖7所示。

圖7 篩選得到的IMFFig.7 The filtered IMF

篩選后的IMF與EEMD分解得到的IMF對應(yīng)關(guān)系，如表2所示。

表2 篩選后的IMF與原IMF對應(yīng)關(guān)系

由表2可得，除原始信號IMF1外，篩選后的IMF3′與原始信號間的相關(guān)系數(shù)最大，相關(guān)性最好，并且其能量占比也同時(shí)達(dá)到最大，所以綜合相關(guān)系數(shù)與能量比值，選取IMF3′作為參考軸的振動分量。對IMF3′做Hilbert變換，得到的時(shí)頻曲線及其擬合曲線如圖8所示?？梢钥吹接捎贖ilbert變換的“端點(diǎn)效應(yīng)”，IMF3′時(shí)頻曲線兩端頻率明顯突變，尤其是末端一些點(diǎn)出現(xiàn)了如圖中“豎線”所示的極大偏差，其擬合曲線兩端的值也因此大幅失真。

圖8 Hilbert變換得到的時(shí)頻曲線及其擬合曲線Fig.8 Time-frequency curve and fitting curve obtained by Hilbert transform

觀察IMF3′的時(shí)頻曲線，分別截取其兩端各1 024個(gè)點(diǎn)進(jìn)行時(shí)頻重排變換。對兩個(gè)重排矩陣直接進(jìn)行瘠線提取，得到如圖9所示的前后兩端截取數(shù)據(jù)的時(shí)頻瘠點(diǎn)。從圖中可以看出，瘠點(diǎn)有一定的波動，并不都處于轉(zhuǎn)頻范圍內(nèi)，應(yīng)該進(jìn)行限制條件的瘠線提取。由于轉(zhuǎn)速范圍為1 500～2 000 r/min，即轉(zhuǎn)軸頻率在25 Hz到33.3 Hz，所以在時(shí)頻矩陣中設(shè)置頻率搜索下限fd=24 Hz，上限fu=34 Hz。由本轉(zhuǎn)子平臺電機(jī)的增速特性，選取m和n的值為2，即對于前端數(shù)據(jù)，在24～26 Hz內(nèi)進(jìn)行搜索；對于后端數(shù)據(jù)，在32～34 Hz內(nèi)進(jìn)行搜索。前后兩端數(shù)據(jù)滿足條件的時(shí)頻瘠點(diǎn)，如圖10所示。

圖9 兩端數(shù)據(jù)直接提取得到的時(shí)頻瘠點(diǎn)Fig.9 Ridge points of both ends extracted directly

圖10 兩端數(shù)據(jù)滿足條件的瘠點(diǎn)Fig.10 Qualified ridge points of both ends

將滿足條件的瘠點(diǎn)進(jìn)行最小二乘擬合，得到前后各1 024點(diǎn)的擬合瘠線，如圖11所示。把截取各點(diǎn)擬合得到的瞬時(shí)頻率值與Hilbert變換得到的瞬時(shí)頻率值進(jìn)行“拼接”。由于Hilbert變換得到的數(shù)據(jù)中間部分的瞬時(shí)頻率仍有一定波動，所以需要將聯(lián)合的頻率進(jìn)一步擬合，以得到光滑的轉(zhuǎn)頻曲線。聯(lián)合頻率曲線及其二階擬合得到的光滑轉(zhuǎn)頻曲線，如圖12所示。

圖11 兩端數(shù)據(jù)的擬合瘠線Fig.11 Fitting ridges of both ends

圖12 聯(lián)合頻率曲線和其擬合曲線Fig.12 Joint frequency curves and its fitting curves

將HHT方法和本文方法估計(jì)的轉(zhuǎn)頻曲線與真實(shí)曲線作比較，如圖13所示。從圖中可以看到，前者最大誤差為7.94 Hz，最大相對誤差率為23.82%；后者最大誤差為0.71 Hz，最大相對誤差率為2.84%，本文方法的兩種參數(shù)指標(biāo)較HHT擬合方法更小，在兩端點(diǎn)處明顯改善了頻率的突變現(xiàn)象，提高了擬合精度，驗(yàn)證了本文方法提取轉(zhuǎn)速曲線相對于HHT方法的優(yōu)越性。

圖13 兩種方法估計(jì)的轉(zhuǎn)頻曲線Fig.13 Speed curve estimated by the two methods

由轉(zhuǎn)頻曲線計(jì)算重采樣時(shí)刻，并對原始信號進(jìn)行等角度重采樣，重采樣結(jié)果如圖14所示。之后，對重采樣數(shù)據(jù)作FFT分析，得到如圖15所示的階比譜。從階比譜中可以看到，信號能量主要集中于基頻，但是有明顯的二階和三階等高次頻率存在，整個(gè)頻譜呈“樅樹形”，符合轉(zhuǎn)子不平衡故障的特征，驗(yàn)證了本方法的有效性。

圖14 重采樣波形Fig.14 Resampling waveform

圖15 階比譜Fig.15 Order spectrum

4 結(jié) 論

(1) 針對HHT變換過程中的模態(tài)混疊、端點(diǎn)效應(yīng)，以及傳統(tǒng)時(shí)頻分析時(shí)頻聚焦性差，計(jì)算效率低等問題，提出EEMD-HHT與時(shí)頻重排相結(jié)合的轉(zhuǎn)速曲線估計(jì)方法。

(2) 利用該方法對某雙跨轉(zhuǎn)子平臺的轉(zhuǎn)子不平衡故障進(jìn)行分析，提取參考軸對應(yīng)的轉(zhuǎn)頻曲線，并與HHT方法直接提取的轉(zhuǎn)頻曲線對比，驗(yàn)證了本方法的優(yōu)越性。將時(shí)域非平穩(wěn)振動信號進(jìn)行等角度間隔重采樣轉(zhuǎn)化為角域平穩(wěn)信號，獲得階比譜，驗(yàn)證了本文所提方法在階比分析中的有效性。