基于EMD和能量算子的模態參數識別在行星齒輪箱中的應用

李康強 ， 馮志鵬

(北京科技大學 機械工程學院，北京 100083)

行星齒輪箱由于其結構緊湊，傳動比大，承載能力強，在車輛、直升機和風力發電等裝備中應用廣泛。作為動力傳動系統中的關鍵環節，行星齒輪箱的安全性對于整個系統的穩定運行具有重要意義。因此，研究行星齒輪箱動力學特性也變得非常重要[1-2]。

通常模態分析通過試驗采集輸入輸出信號經過參數識別來獲得模態參數(模態頻率、模態阻尼比和模態振型)[3-4]。行星齒輪箱由于其結構復雜、多自由度、多模態，直接進行模態參數識別易造成模態混疊產生較大誤差。而EMD能夠自適應分解出具有局部特征信號的本質模函數(IMF)，可以初步解決模態混疊。因此首先應用EMD進行模式分解出多個本質模函數便于下一步分析和定階。針對模態頻率的識別，傳統方法為根據輸入輸出信號求出頻響函數，再由傅里葉變換得出模態頻率[5-6]。能量算子能追蹤信號能量計算瞬時頻率，而高階能量算子在此基礎提高了計算精度，因此我們提出EMD與高階能量算子(HEO)相結合的方法，利用高階能量算子[7-9]的優勢估算EMD分解出的模函數的頻率，同時根據估算出的頻率來確定模態階次。由于EMD分解過程可能會產生虛假分量，因此計算各個IMF與原信號的相關性，同時根據各個模函數的頻率來識別及去除虛假分量。

對于阻尼比的識別的研究，目前的計算方法主要通過參數識別理論，通過實際測量的數據識別阻尼比。傳統方法有對數衰減率法[10]、ITD(Ibrahim Time Domain Technique，ITD)法[11]、STD(Spare Time Domain Technique，STD)法[12]、半功率帶寬法[13]等，這些方法都存在精度與魯棒性不能并存的缺點。由于阻尼現象就是反映了振動結構耗散的能量，所以我們可以從能量損耗的角度來計算阻尼比，Marceau[15]提到了信號能量，并應用信號能量方法來估計電力系統的暫態穩定極限。曾儲惠等[16]提出應用信號能量識別阻尼比方法，創新性和魯棒性較強。受益于前人的研究，發揮能量算子追蹤系統能量的優勢和周期衰減的恒定性，我們提出一種新的基于能量算子計算結構阻尼比的方法：半周期能量算子法(Half Cycle Energy Operator, HCEO)。用能量比來估算阻尼比，其抗噪性和魯棒性將大為提高。

通過分析仿真信號驗證了本方法在識別模態頻率和阻尼比的有效性。應用于實驗實測行星齒輪箱模態信號，成功估算了各階次模態頻率及阻尼比。

1 基本理論

1.1 經驗模態分解(EMD)

經驗模式分解[17-18]可以根據信號自身的特點自適應地獲得基函數，通過對非線性、非平穩信號的分解獲得一系列表征信號特征時間尺度的固有模式函數。經驗模式分解實際上是對信號進行了平穩化和線性化處理。

本質模式函數滿足以下兩個條件：

1)整個IMF中，極值點數與過零點數相等或至多相差一個；

2)在任意時刻，由局部極大值點和局部極小值點分別確定的上、下包絡線的均值均為零。即信號關于時間軸局部對稱。

EMD分解具體實施步驟如下：

(1) 確定信號所有的局部極大值點和極小值點，用三次樣條曲線分別擬合所有的局部極大值點和極小值點作為上包絡線和下包絡線。計算上下包絡線的平均值m1，用原始信號x(t)減去上下包絡線平均值得到

h1=x(t)-m1

(1)

如果h1滿足IMF的兩個條件，那么h1就是原始信號的第一個IMF分量。

(2)如果h1不滿足IMF的條件，則將其視為原始信號，重復步驟(1)的篩選，直到得到的hk滿足IMF的條件。記c1=hk，則c1為原始信號x(t)的第一個IMF分量。

(3)原始信號減去c1得到第一階剩余信號r1。將其作為第二階原始信號，重復步驟(1)、(2)，得到第二個滿足IMF條件的分量c2，重復循環n次，直到當rn成為一個單調函數不能再分解時，循環結束。得到原始信號的n個滿足IMF條件的分量。分解后原信號可以表示為：

(2)

其中本質模式函數c1,c2…cn代表著從高到低的頻率組分。rn為殘余分量。

由以上經驗模式分解的分解過程可以看出，經驗模式分解的基函數是一系列幅值和頻率可變的正弦或余弦函數，它依賴信號本身自適應地獲得，因此分解效果更好。

1.2 能量算子(EO)

對于連續信號x(t)，Teager能量算子Ψ定義為：

(3)

對離散信號可將用差分替代微分，可得離散能量算子：

Ψ[x(n)]=[x(n)]2-x(n+1)x(n-1)

(4)

對于任意連續信號，高階微分能量算子定義為：

(5)

式中：k為高階能量算子的階次；x(k)為相對于時間的第k階導數。

同理離散信號的離散高階能量算子為：

Υk(x[n])=x[n]x[n+k-2]-x[n-1]x[n+k-1]

(6)

對于簡諧信號x=Acos(ωt+φ)，由四階能量算子和二階能量算子可以估算出信號的頻率和幅值：

(7)

需要注意的是，根據高階能量算子估算出的頻率為瞬時頻率，還需根據結果截取平穩階段再求平均得到模態頻率。

1.3 半周期能量算子法(HCEO)

由于阻尼就是反映了振動結構耗散的能量，而能量算子具有追蹤信號能量的性質，基于以上理論提出根據半周期性的能量描述來提取阻尼比。

半周期能量定義為：

(8)

式中：ti為某一點過零點時刻；td為振蕩周期。

對于單自由度自由衰減振動響應模型

x(t)=Ae-ξωntsin(ωdt+φ)

(9)

(10)

繼而延伸通過n個半周期來計算阻尼比。

(11)

而能量算子具有跟蹤系統能量的性質。

(12)

則式(10)可化為：

(13)

令

(14)

又由

(15)

則聯立式(12)～(14)可得：

(16)

式中：Ψx(1)、Ψx(2)、…Ψx(2N)代表的物理含義為：信號x(t)每半個周期的能量算子跟蹤能量值。其中周期可由式(17)可以計算得出。

(17)

2 分析步驟

(1) 為了減少計算量和避免弱衰減和零衰減部分對結果的影響，對原始信號截取明顯衰減的部分進行EMD分解得到多個IMF分量。

(2) 離散信號的EMD分解可能造成取到異常極值得到虛假分量，因此我們計算各IMF與原信號的相關系數，根據相關性結合估算出的頻率值來確定虛假分量。其中相關性閾值設定為λ=max(covf(i))/β,β=2～10，其中β的取值可以根據實際頻率值確定，一般取5或者10。如果閾值小于λ，考慮可能為虛假分量。

(3) 由高階能量算子估算每個IMF分量的頻率，根據分量波形特征及對應的頻率可以對分量虛假性進行進一步判斷。

(4) 對各個IMF分量用HCEO法估算阻尼比。同時與傳統方法對比分析后結合步驟(2)和步驟(3)確定模態階次及各階次模態頻率與模態阻尼比。

3 仿真信號分析

不失一般性，仿真信號由幾個單自由度自由衰減振動信號疊加組成：

(18)

式中：Am為第m階幅值；ξm為第m階阻尼比；ωnm為第m階無阻尼固有頻率；ωdm為第m階有阻尼固有頻率；φm為初相位，仿真信號中分別加入信噪比SNR=10 dB、30 dB和50 dB來進行對比分析。具體取值見表1。

表1 仿真信號參數Tab.1 Parameters of simulated signal

仿真信號衰減時域如圖1所示。首先對信號進行EMD分解，分解結果如圖2，分別用傳統方法(頻響函數傅里葉變換(FRFT))、希爾伯特變換法(HT)和高階能量算子(HEO)估算分解后的各階次模態頻率，結果見表3。對應各階次IMF分別用半周期能量算子法(HCEO)、傳統半周期法(HCE)、對數衰減率法(LD)和半功率帶寬法(HF)計算結構阻尼比，結果見表4。其中計算結果的階次順序是按照EMD分解出的IMF順序排列，所以模態階次1～6階也是從高階到低階排列。

圖1 仿真信號(SNR=30 dB) Fig.1 simulated signal (SNR=30 dB)

圖2 EMD分解結果(SNR=30 dB) Fig.2 IMFs of EMD (SNR=30 dB)

SNR/dBIMF1IMF2IMF3IMF4IMF5IMF6100.21080.47900.39750.19750.82660.1209300.56790.26670.10770.73740.1386-0.1601500.35480.48370.60290.14250.11850.1075

表3 仿真信號模態頻率Tab.3 Modal frequencies of simulated signal

表4 仿真信號阻尼比Tab.4 Damping ratios of simulated signal

圖3 HEO誤差趨勢 Fig.3 Errors of HED

圖4 HCEO誤差趨勢 Fig.4 Errors of HCED

取其中信噪比 SNR=30 dB 時的信號的EMD分解結果如圖2所示，通過IMF時域圖可以看到c3、c5和c6的波形特征不明顯和沒有呈現很好的衰減性特征，初步判斷可能是虛假分量。各個IMF分量與原信號的相關系數如表2所示，由于仿真信號由三個分量組成，結合計算出的頻率可以設置β=5，得出閾值λ=0.147，小于閾值的IMF為c3、c5和c6，結合表3的頻率估計統計可以斷定真實分解的分量只有c1、c2和c4。結合頻率估計和相關性分析可以判斷 SNR=10 dB時的真實階次為c2、c3和c5。當SNR=50 dB 時的真實階次為c1、c2和c3。提取真實階次重新統計模態頻率與阻尼比見表5和表6。通過傳統方法(FRFT)、希爾伯特變換法(HT)和高階能量算子(HEO)三種方法估算頻率能力的對比可以看出，結果相差無幾，說明實際運算中估算頻率的這幾種方法可以綜合借鑒。針對我們使用的HEO方法來說，圖3中可以看出，HEO方法總體誤差率較低，而且精度隨著信噪比增加而升高。對于比較半周期能量算子法(HCEO)、傳統半周期法(HCE)、對數衰減率法(LD)和半功率帶寬法(HF)計算結構阻尼比如表6所示，可以看出相比HCE法和LD法，HCEO方法總體上是比較穩定，魯棒性較好。其中表4中出現負值的結果是由于上半周期的總能量接近或者略小于下半周期的總能量，所以其比值的對數為負，導致最終計算結果會出現負數。而HF方法是對幅頻特性的阻尼比識別，所以噪聲對其影響較小。但是這種方法對于小阻尼比的時候精度較高，隨著阻尼比的增大，計算精度也大幅下降，而且選用的采樣點數和采樣頻率會對其結果影響很大，甚至可能出現無解的情況。從圖4中的HCEO方法估算阻尼比誤差趨勢可以看出，其精度與信噪比大小成正比，而且較其他兩種方法穩定，計算效果良好。

表5 仿真信號實際階次模態頻率Tab.5 Modal frequencies of simulated signal

表6 仿真信號實際階次阻尼比Tab.6 Damping ratios of simulated signal

4 實測模態試驗信號分析

4.1 實驗系統

某行星齒輪箱模態試驗采集系統如圖5所示，加速度傳感器分別布置在齒輪箱箱體的水平徑向、水平軸向和垂直徑向見圖中所標。所用輸入、輸出傳感器參數如表7所示。

圖5 實驗系統 Fig.5 Experimental system

激勵(輸入)使用力錘進行諧振頻率測試和模態分析實驗。利用分布在箱體水平徑向、水平軸向和垂直方向的加速度傳感器采集(輸出)振動信號。使用力錘分別從所布傳感器的對立面的水平徑向、水平軸向和垂直方向對行星齒輪箱進行錘擊實驗，并采集敲擊過程中的振動信號。為保證信號的質量，每個方向重復三次實驗。需要注意的是靜態測試中使用力錘敲擊時應避免二次敲擊。

表7 傳感器參數Tab.7 Parameters of sensors

4.2 X方向信號分析

從水平徑向(X方向)對行星齒輪箱進行敲擊并采集振動信號，時域信號如圖6所示，利用EMD對其進行分解，結果如圖7所示。

圖6 時域波形 Fig.6 Time-domain wave

圖7 EMD分解結果 Fig.7 IMFs of EMD

取前6個IMF分別計算與原信號的相關系數如表8。用傳統方法(FRFT)、希爾伯特變換法(HT)和高階能量算子(HEO)估算模態頻率，結果見表9。對應各階次IMF分別用半周期能量算子法(HCEO)、傳統半周期法(HCE)、對數衰減率法(LD)和半功率帶寬法(HF)計算結構阻尼比，結果見表10。

通過圖7我們無法判斷虛假分量。首先估算頻率如表9所示，根據表9的頻率特征和我們已經掌握的行星齒輪箱的其他數據，在此可以取閾值參數β=10，得出閾值λ=0.094 5。通過表8的相關系數可以看出，只有前三個分量大于這個閾值。而通過表9的頻率估計統計結果可看出，第一階次和第二階次的不同方法所求的頻率相差無幾。對所有IMF計算阻尼比的結果如表10所示，可以看出第一階次中各個方法的數據遠小于其他階次數據，其中HCE方法的第一階次出現負值，這是由于上半周期的總能量接近或者略小于下半周期的總能量，所以其比值的對數為負導致最終計算結果會出現負數，再次結合表9中頻率結果可以判斷第一階次為虛假分量。最終判斷分解出的分量只有c2、c3結果符合真實階次。

表8 相關系數Tab.8 Correlation

表9 X方向信號模態頻率Tab.9 Modal frequencies of X direction signal

表10 X方向信號阻尼比Tab.7 Damping ratios of X direction signal

對于比較結構阻尼結果，可以看出對于HCEO、HCE和LD方法的c2、c3、c4結果比較穩定，但是第4階次之后阻尼比就大幅增大。而其中HCEO的阻尼比結果比較穩定。相比其他三種方法，頻域HF方法結果總體上偏大，這是由于在低頻率的計算過程中尋找半功率點的誤差較大所導致。最終得到真實階次模態及其參數如表11所示。

表11 X方向模態參數Tab.11 Modal parameters of X direction signal

4.3 Y方向信號分析

從水平軸向(Y方向)對行星齒輪箱進行敲擊并采集振動信號，時域振動信號如圖8所示，利用EMD對其進行分解，分解結果如圖9所示。

圖8 時域波形 Fig.8 Time-domain wave

圖9 EMD分解結果 Fig.9 IMFs of EMD

取前6個IMF分別用傳統方法(FRFT)、希爾伯特變換法(HT)和高階能量算子(HEO)估算模態頻率，結果見表13。對應各階次IMF分別用半周期能量算子法(HCEO)、傳統半周期法(HCE)、對數衰減率法(LD)和半功率帶寬法(HF)計算結構阻尼比，結果見表14。

表12 相關系數Tab.12 Correlation

表13 Y方向信號模態頻率Tab.13 Modal frequencies of Y direction signal

表14 Y方向信號阻尼比Tab.14 Damping ratios of Y direction signal

在水平軸向測試過程中，由于振動傳遞路徑一部分是遍布凹槽的外箱體，另一部分是經過太陽輪和行星輪系的傳遞路徑，所以Y方向振動信號遠比X方向信號更為復雜。

首先計算相關性閾值為0.054，通過表12的相關系數可以看出，IMF6與原信號為負相關，而且取絕對值之后也是小于相關性閾值，所以可直接判斷為虛假分量。通過圖9可以看出c2沒有明顯的衰減波形，初步判斷可能是分解出的虛假分量，而通過表14可以發現c2的幾種方法的阻尼比均較高，由于通常我們實驗用的行星齒輪箱的阻尼比不會超過0.1，再結合c2時域波形所以可以判斷c2為虛假分量。通過表12的相關系數可以看出c3與原信號相關性最高，同時通過表14可以發現c3的阻尼比與X方向的結果接近，可以判斷c3得出的結果為正確的。而分量c2、c4、c5、c6均可判斷為虛假分量。最終判斷分解出的分量只有c1、c3結果符合真實階次。得到真實階次模態及其參數如表15所示。

表15 Y方向模態參數Tab.15 Modal parameters of Y direction signal

4.4 Z方向信號分析

由于垂直方向無法建立敲擊點與采集點為對立面關系，所以敲擊點和傳感器都在齒輪箱體的正上方。從垂直徑向(Z方向)對行星齒輪箱進行敲擊并采集振動信號，時域信號如圖10所示，利用EMD對其進行分解，結果如圖11所示。

圖10 時域波形 Fig.10 Time-domain wave

圖11 EMD分解結果 Fig.11 IMFs of EMD

取前6個IMF分別計算相關系數見表16。用FRFT法、HT法和HEO法估算模態頻率結果見表17。對應各階次IMF分別用HCEO法、HCE法、LD法和HF法計算結構阻尼比，結果見表18。

表16 相關系數Tab.16 Correlation

表17 Z方向信號模態頻率Tab.17 Modal frequencies of Z direction signal

表18 Z方向信號阻尼比Tab.18 Damping ratios of Z direction signal

由于敲擊點離傳感器非常近，所以時域波形衰減信號非常清晰。通過圖11的EMD分解結果可以看出第一階次分解效果最好，其他階次暫時無法估計是否虛假分量。通過表16的相關系數可以很明顯看出c1和c2與原函數相關性很高，其他分量相關性較小(相關性閾值為0.072)。通過表17的模態頻率估計統計可以看出第一階次和第二階次的結果相近，結合X方向與Y方向的分析結果和第二階IMF不完整的時域波形可以判斷第二階次為虛假分量。其他分量通過時域波形(自由衰減不明顯或者波形不完整)、頻率和阻尼比綜合考慮，同樣都可以判斷為虛假模態分量。最終得到Z方向信號的真實階次模態及其參數如表19。

表19 Z方向模態參數Tab.19 Modal parameters of Z direction signal

5 結 論

本文針對多自由度行星齒輪箱的模態參數識別問題，提出一種將經驗模式分解與能量算子相結合的模態識別方法：用EMD結合高階能量算子估算模態頻率；用EMD結合半周期能量算子估算結構阻尼比。通過本方法可以識別模態頻率和結構阻尼比，同時可以確定模態階次和判斷模態混迭。使用不同信噪比的仿真信號驗證了本方法，通過對比傳統方法，體現本算法的簡單易行，抗噪性和魯棒性較強。然后將方法應用于行星齒輪箱模態實驗信號中，包括水平徑向、水平軸向和垂直方向對行星齒輪箱的錘擊實驗信號，綜合發現行星齒輪箱在不同方向計算的模態頻率略有差別，但是阻尼比相差不大，均在0.03～0.05之間。本方法能有效提取各階次模態頻率和模態阻尼比并判斷出 虛假分量，驗證了方法的有效性與實用性。

[ 1 ] PARKER R G, AGASHE V, VIJAYAKAR S M. Dynamic response of a planetary gear system using a finite element/contact mechanics model[J]. Journal of Mechanical Design, 2000, 122(3): 304-310.

[ 2 ] FENG Z, ZUO M J. Vibration signal models for fault diagnosis of planetary gearboxes[J]. Journal of Sound and Vibration, 2012, 331(22): 4919-4939.

[ 3 ] 許本文, 焦群英. 機械振動與模態分析基礎[M]. 北京: 機械工業出版社, 1998.

[ 4 ] 曹樹謙, 張文德, 蕭龍翔. 振動結構模態分析: 理論、實驗與應用[M]. 天津： 天津大學出版社, 2014.

[ 5 ] LIN J, PARKER R G. Sensitivity of planetary gear natural frequencies and vibration modes to model parameters[J]. Journal of Sound and Vibration, 1999, 228(1): 109-128.

[ 6 ] LIN J, PARKER R G. Natural frequency veering in planetary gears[J]. Mechanics of Structures and Machines, 2001, 29(4): 411-429.

[ 7 ] MARAGOS P, POTAMIANOS A．Higher order differential energy operators[J]．IEEE Signal Processing, 1995, 2(8)： 152-154．

[ 8 ] FENG Z, QIN S, LIANG M. Time-frequency analysis based on Vold-Kalman filter and higher order energy separation for fault diagnosis of wind turbine planetary gearbox under nonstationary conditions[J]. Renewable Energy, 2016, 85: 45-56.

[ 9 ] FAGHIDI H, LIANG M. Bearing fault identification by higher order energy operator fusion: A non-resonance based approach[J]. Journal of Sound & Vibration, 2016, 381: 83-100.

[10] COOPER J E. Extending the logarithmic decrement method to analyse two degree of freedom transient responses[J]. Mechanical Systems and Signal Processing, 1996, 10(4): 497-500.

[11] IBRAHIM S R, MIKULCIK E C. A method for the direct identification of vibration parameters from the free response[J]. The Shock and Vibration Bulletin, 1977, 47: 183-

198.

[12] IBRAHIM S R. An approach for reducing computational requirements in modal identification[J]. Aiaa Journal, 1986, 24(10): 1725-1727.

[13] SEYBERT A F. Estimation of damping from response spectra[J]. Journal of Sound and Vibration, 1981, 75(2): 199-

206.

[14] LE D H, CHENG J, YANG Y, et al. Gears fault diagnosis method using ensemble empirical mode decomposition energy entropy assisted ACROA-RBF neural network[J]. Journal of Computational & Theoretical Nanoscience, 2016, 13(5): 3222-3232.

[15] MARCEAU R J, SIRANDI M, SOUMARE S, et al. A review of signal energy analysis for the rapid determination of dynamic security limits[J]. Canadian Journal of Electrical and Computer Engineering, 1996, 21(4): 125-132.

[16] 曾儲惠, 黃方林, 柳成蔭, 等. 基于信號能量分析的結構阻尼比識別方法[J]. 振動與沖擊, 2003, 22(2): 66-68.

ZENG Chuhui, HUANG Fanglin, LIU Chengyin, et al. Damping ratio identification method based on signal energy analysis[J]. Journal of Vibration and Shock, 2003, 22(2): 66-68.

[17] HUANG N E, SHEN Z, LONG S R, et al. The empirical mode decomposition and the Hilbert spectrum for nonlinear and non-stationary time analysis[J]. Proc. R. Soc. Lond. A, 1998, 454(1971): 903-955.

[18] LOH C H. Application of the empirical mode decomposition-hilbert spectrum method to identify near-fault ground-motion characteristics and structural responses[J]. Journal of Loss Prevention in the Process Industries, 2016, 91(5): 1339-1357.