一種EEMD-HHT與時頻重排結合的轉速曲線估計新方法

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

(1. 天津工業大學 電子與信息工程學院，天津 300387； 2. 天津市光電檢測技術與系統重點實驗室，天津 300387)

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

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

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

1 EEMD-HHT和時頻重排算法原理

1.1 EEMD-HHT算法

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

1.2 時頻重排算法

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

(1)

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

(2)

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

(3)

(4)

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

(5)

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

2 EEMD-HHT與時頻重排結合的估計方法

針對以上問題，本文提出了將EEMD-HHT與時頻重排相結合的方法來估計參考軸的轉速曲線，轉頻曲線估計算法實現的流程圖如圖1所示。其詳細步驟如下：

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

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

(6)

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

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

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

(7)

式中：X表示某個IMF分量，M為IMF的個數。ek越大表示第K個IMF分量所占比重越大，越小則表示所占比重越小。由于振動信號是由參考軸轉動提供動力產生的，所以其對應的IMF相關系數與能量占比兩個數據都可能更大。綜合考慮兩種計算結果，選取相關系數值和能量占比都較大的IMF分量作為參考軸的本征模態函數。

(3) 對選擇的IMF分量作Hilbert變換，得到其各個點的瞬時頻率值。

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

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

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

(8)

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

(9)

(6) 用第5步得到IMF前后兩端截取各點的瞬時頻率值代替第3步Hilbert變換得到的對應點處的值，組成一個新的頻率數組。

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

3 應用實例

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

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

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

對振動信號直接進行EMD分解，得到的IMF如圖4所示。從圖中可以看到，相鄰IMF間產生了模態混疊現象，使IMF失去了真實意義。

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

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

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

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

表1 各IMF與原信號間的相關系數值

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

圖6 能量譜圖Fig.6 Energy spectrum

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

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

篩選后的IMF與EEMD分解得到的IMF對應關系，如表2所示。

表2 篩選后的IMF與原IMF對應關系

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

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

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

圖9 兩端數據直接提取得到的時頻瘠點Fig.9 Ridge points of both ends extracted directly

圖10 兩端數據滿足條件的瘠點Fig.10 Qualified ridge points of both ends

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

圖11 兩端數據的擬合瘠線Fig.11 Fitting ridges of both ends

圖12 聯合頻率曲線和其擬合曲線Fig.12 Joint frequency curves and its fitting curves

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

圖13 兩種方法估計的轉頻曲線Fig.13 Speed curve estimated by the two methods

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

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

圖15 階比譜Fig.15 Order spectrum

4 結 論

(1) 針對HHT變換過程中的模態混疊、端點效應，以及傳統時頻分析時頻聚焦性差，計算效率低等問題，提出EEMD-HHT與時頻重排相結合的轉速曲線估計方法。

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