， ，

(廣東海洋大學(xué)，廣東 湛江 524088)

0 引言

柴油機(jī)氣缸內(nèi)氣體壓力周期性變化，運(yùn)動(dòng)部件重力及其往復(fù)慣性力周期性變化，接受功率的部件不能均勻地吸收扭振所產(chǎn)生的激勵(lì)等可能會(huì)導(dǎo)致船用柴油機(jī)曲軸、傳動(dòng)軸以及凸輪軸疲勞折損。如果扭振幅值過大將可能會(huì)激起柴油機(jī)機(jī)架與齒輪箱的橫向振動(dòng)，引起機(jī)艙構(gòu)件的局部振動(dòng)以及上層建筑及船體振動(dòng)。監(jiān)測(cè)船舶軸系的固有頻率和船上有關(guān)的激勵(lì)頻率之間是否出現(xiàn)共振，特別是監(jiān)測(cè)船舶推進(jìn)軸系的低頻振動(dòng)顯得尤為重要。

本文通過數(shù)值仿真實(shí)驗(yàn)，明確固有模態(tài)函數(shù)與系統(tǒng)振動(dòng)頻率的對(duì)應(yīng)關(guān)系，即一個(gè)固有模態(tài)函數(shù)與系統(tǒng)中某一個(gè)頻率不是一一對(duì)應(yīng)的。通過對(duì)扭轉(zhuǎn)振動(dòng)進(jìn)行測(cè)試分析，探討了經(jīng)驗(yàn)?zāi)B(tài)分解在推進(jìn)軸系扭轉(zhuǎn)振動(dòng)分析中的應(yīng)用。

1 經(jīng)驗(yàn)?zāi)B(tài)分解與希爾伯特黃變換

經(jīng)驗(yàn)?zāi)B(tài)分解(Empirical Mode Decomposition, EMD)計(jì)算數(shù)據(jù)的上下包絡(luò)平均值，將包絡(luò)平均值作為瞬時(shí)平均值。用原始信號(hào)減去瞬時(shí)平均值，得到新的數(shù)據(jù)，判斷新數(shù)據(jù)是否滿足固有模態(tài)函數(shù)(Intrinsic Mode Function, IMF)的定義，如果新數(shù)據(jù)是固有模態(tài)函數(shù)，就從原始信號(hào)中減去新數(shù)據(jù)，得到殘留數(shù)據(jù)，直到殘留數(shù)據(jù)少于2個(gè)極值點(diǎn)，終止分解運(yùn)算[1-2]。分解過程見公式1。

如果一個(gè)函數(shù)的極值點(diǎn)數(shù)目和過零點(diǎn)數(shù)目相等或者最多相差1個(gè)；而且這個(gè)函數(shù)在任意點(diǎn)，由局部極大值和局部極小值點(diǎn)構(gòu)成的2條包絡(luò)線平均值為0。那么，這個(gè)函數(shù)稱為固有模態(tài)函數(shù)[2]。

EMD方法是一個(gè)篩選過程,采用信號(hào)局部極值和緊隨其后的樣條函數(shù)擬合來近似分解[3]。EMD是一個(gè)數(shù)據(jù)驅(qū)動(dòng)和自適應(yīng)的過程[4]。每個(gè)線性或非線性模式將具有相同數(shù)量的極值點(diǎn)和零交叉點(diǎn),每個(gè)模式都獨(dú)立于其他模式[5-6]。

x(t)=imf1(t)+r1(t)

=imf1(t)+imf2(t)+r2(t)

=imf1(t)+imf2(t)+imf3(t)+r3(t)

(1)

?

希爾伯特黃變換(Hilbert-Huang Transform, HHT)是將所要分析的數(shù)據(jù)分解為IMF后，再對(duì)每一個(gè)IMF做HHT，從而正確地獲得信號(hào)的瞬時(shí)頻率。瞬時(shí)頻率定義為瞬時(shí)相位的一階導(dǎo)數(shù)[7]。

數(shù)據(jù)x(t)的希爾伯特黃變換定義如下：

(2)

重構(gòu)信號(hào)z(t)為解析信號(hào)，見公式(3)。

z(t)=x(t)+jy(t)

=a(t)eiθ(t)

(3)

那么，瞬時(shí)幅值見公式(4)。

(4)

瞬時(shí)相位見公式(5)。

(5)

瞬時(shí)頻率見公式(6)。

(6)

由于只能對(duì)有限個(gè)樣本進(jìn)行計(jì)算，當(dāng)時(shí)域截?cái)鄷r(shí)，傳統(tǒng)的FFT，不可避免地產(chǎn)生能量泄漏，譜峰值變小，從而產(chǎn)生幅值誤差[8]。

與其他數(shù)學(xué)變換(如快速傅立葉變換(FFT)、小波變換(WAVELET)等)不同，希爾伯特黃變換是一種應(yīng)用在數(shù)據(jù)上的算法，而非理論工具，其處理對(duì)象是非穩(wěn)態(tài)與非線性信號(hào)。這種基本上符合物理意義的分解，是其他變換難以實(shí)現(xiàn)的[1]。HHT變換可以跟蹤瞬時(shí)頻率的變化,在處理變壓器振動(dòng)信號(hào)時(shí)，HHT具有獨(dú)特的優(yōu)勢(shì)[7]。

2 數(shù)值仿真實(shí)驗(yàn)

信號(hào)由不同頻率的正弦波組成，并含有隨機(jī)噪聲。

signal1=sin(2*pi*10*t)+sin(2*pi*50*t)+sin(2*pi*100*t)+sin(2*pi*150*t)+sin(2*pi*200*t)+sin(2*pi*250*t)+sin(2*pi*300*t)+sin(2*pi*350*t)+sin(2*pi*400*t)+sin(2*pi*450*t)+0.1*randn(1, length(t));

信號(hào)signal1的傅立葉變換見圖1。在頻譜圖上，10個(gè)不同頻率明顯被分辨出來。

圖1 signal1與其傅立葉譜

信號(hào)signal1的EMD分解見圖2，從圖中可以看出，第一個(gè)IMF的頻譜只分辨出9個(gè)頻率，第二個(gè)IMF，第三個(gè)IMF雖然有表征，但未能明顯表征10 Hz頻率。EMD對(duì)系統(tǒng)內(nèi)不同頻率振動(dòng)能量差別不大時(shí)，頻率特征識(shí)別能力較差。

圖2 signal1的IMF與其頻譜

另一方面,在采集的數(shù)據(jù)集上，固有模態(tài)函數(shù)(IMF)分量不能自動(dòng)保證有一個(gè)定義明確的物理意義[4]，因此采用EMD分解IMF分量前，需要仔細(xì)選擇所采集數(shù)據(jù)的尺度，以便有明確的解釋。

如果把含有隨機(jī)噪聲的信號(hào)更改如下：

signal2=10*sin(2*pi*10*t)+2*sin(2*pi*50*t)+sin(2*pi*100*t)+20*sin(2*pi*150*t)+8*sin(2*pi*200*t)+30*sin(2*pi*250*t)+20*sin(2*pi*300*t)+sin(2*pi*350*t)+8*sin(2*pi*400*t)+20*sin(2*pi*450*t)+0.1*randn(1, length(t))

不同于信號(hào)signal1，信號(hào)signal2的各頻率成分的振幅(或振動(dòng)能量)相差比較大。在信號(hào)signal2的快速傅立葉變換頻譜圖(圖3)上，8個(gè)不同頻率較明顯地被分辨出來,但100 Hz，350 Hz不能明顯表征。

圖3 signal2與其傅立葉譜

信號(hào)signal2的EMD分解見圖4，從圖4中可以看出，第一個(gè)IMF的頻譜只分辨出9個(gè)頻率，但第二個(gè)IMF對(duì)100 Hz頻率做了補(bǔ)充表征，第三個(gè)IMF對(duì)50 Hz頻率做了補(bǔ)充表征，第四個(gè)IMF對(duì)10 Hz頻率做了補(bǔ)充表征。經(jīng)多次數(shù)值分析實(shí)驗(yàn)，EMD對(duì)各頻率振動(dòng)能量差別大的系統(tǒng)，有較好的頻率特征識(shí)別能力。

圖4 signal2的IMF與其頻譜

如果把不同頻率的振幅再加大，把含有隨機(jī)噪聲的信號(hào)更改如下：

signal3=20*sin(2*pi*10*t)+0.5*sin(2*pi*50*t)+10*sin(2*pi*100*t)+15*sin(2*pi*150*t)+10*sin(2*pi*200*t)+20*sin(2*pi*250*t)+0.1*sin(2*pi*300*t)+20*sin(2*pi*350*t)+10*sin(2*pi*400*t)+20*sin(2*pi*450*t)+0.1*randn(1, length(t))

圖5 signal3與其傅立葉譜

從圖5與圖6可以看出，快速傅立葉變換對(duì)50 Hz,300 Hz沒有明顯表征。而EMD分解對(duì)100 Hz,50 Hz在IMF3做了補(bǔ)充表征，對(duì)10 Hz在IMF4上做了補(bǔ)充表征。但300 Hz頻率補(bǔ)充表征不是很明顯?？傊?，EMD分解寬帶信號(hào)時(shí)，是從高頻到低頻, 而不是從高能量到低能量[3]。這個(gè)特征使得EMD能夠分解出低頻弱能量信號(hào)[3]。

圖6 signal3的IMF與其頻譜

如果信號(hào)更改如下：

a=2*t+3;

signal4=a.*sin(2*pi*200*t+pi*t)

信號(hào)signal4是幅值相位隨時(shí)間變化的信號(hào)，從圖7、圖8中可以看出，快速傅立葉變換與EMD分解都能明顯地表征200 Hz頻率這個(gè)單一頻率。

圖7 signal4與其傅立葉譜

圖8 signal4的IMF與其頻譜

如果信號(hào)更改如下：

a=2*t+3;

signal5=sin(2*pi*200*a.*t+pi*t);

信號(hào)signal5是頻率、相位隨時(shí)間變化的，從圖9、圖10中可以看出，EMD分解能較好地表征低頻弱能量信號(hào)。

圖9 signal5與其傅立葉譜

圖10 signal5的IMF與其頻譜

EMD分解的IMF能量和原信號(hào)的能量(可能)是不完全相等[8]。FFT譜和原信號(hào)的能量是相等的。從EMD的定義與數(shù)值分析實(shí)驗(yàn)表明，一個(gè)固有模態(tài)函數(shù)與系統(tǒng)中某一個(gè)頻率不是一一對(duì)應(yīng)。對(duì)于不復(fù)雜的系統(tǒng)，一般情況下，前4階固有模態(tài)函數(shù)基本上能夠表征系統(tǒng)的模態(tài)。

EMD分解后的固有模態(tài)函數(shù)是不等帶寬的，所以，EMD分解適合非線性信號(hào)分析[6]。

3 實(shí)測(cè)扭振數(shù)據(jù)分析

電機(jī)型號(hào)：Y2-90S-4，電機(jī)功率：1.1 kW，頻率50 Hz,額定轉(zhuǎn)速：1 440 r/min,軸直徑：φ30，中間齒輪齒數(shù)：60。

實(shí)驗(yàn)測(cè)試儀器采用ANZT雙通道扭振分析記錄儀，扭角測(cè)試：量程：0-10°(峰值)；分辨率：1毫度。測(cè)試平臺(tái)見圖11。

圖11 實(shí)驗(yàn)測(cè)試臺(tái)

實(shí)測(cè)數(shù)據(jù)及其傅立葉譜，見圖12。

扭振信號(hào)的頻率成分主要是由低頻的滾振成分和與轉(zhuǎn)速相關(guān)的各諧次成分組成。從圖12、圖13實(shí)測(cè)的數(shù)據(jù)來看，在啟動(dòng)過程中，扭振信號(hào)是一組低頻的變頻信號(hào)。

圖12 實(shí)測(cè)數(shù)據(jù)及其傅立葉譜

圖13 實(shí)測(cè)數(shù)據(jù)的IMF與其頻譜

圖14 實(shí)測(cè)數(shù)據(jù)IMF1的瞬時(shí)頻率

從圖14中看出，所測(cè)扭振信號(hào)頻譜成分復(fù)雜，與頻率變化(非平穩(wěn))的信號(hào)有類似的譜圖，低頻成份較多。與實(shí)際相比，F(xiàn)FT變換僅能夠明晰趨勢(shì)的變化，不能反映啟動(dòng)工況下扭振信號(hào)的實(shí)際物理意義。而EMD，除去其端點(diǎn)效應(yīng)，第一IMF最大扭角達(dá)到350毫度，第二IMF最大扭角達(dá)到210毫度，第三個(gè)IMF最大扭角達(dá)到10毫度，與實(shí)測(cè)數(shù)據(jù)基本相同，基本上表征了扭振實(shí)際過程的特性。從圖13實(shí)測(cè)數(shù)據(jù)IMF1的瞬時(shí)頻率中可以看出，啟動(dòng)過程中，轉(zhuǎn)速達(dá)到200 r/min以上后，扭振頻率沒有太大變化。

4 結(jié)論

EMD分解適合非線性與非平穩(wěn)信號(hào)，所分解的固有模態(tài)函數(shù)與系統(tǒng)中的固有頻率不是一一對(duì)應(yīng)的，通常情況下，前4階固有模態(tài)函數(shù)基本上能夠表征系統(tǒng)的模態(tài)。

在傳統(tǒng)的穩(wěn)態(tài)扭振信號(hào)處理中，首先采集各穩(wěn)態(tài)轉(zhuǎn)速下扭振信號(hào)；然后提取各諧次幅值的譜分析；最后對(duì)各穩(wěn)態(tài)轉(zhuǎn)速下各諧次幅值進(jìn)行整理，并用多項(xiàng)式曲線擬合得到轉(zhuǎn)速振幅圖。

本文由于采用ANZT雙通道扭振分析記錄儀，得到扭振數(shù)據(jù)并將數(shù)據(jù)導(dǎo)入MATLAB，通過編程分析，驗(yàn)證了EMD對(duì)扭振測(cè)試數(shù)據(jù)有較好的分析能力。

[1] Jian Zhang, Ruqiang, et al. Performance enhancement of ensemble empirical mode decomposition[J]. Mechanical System and Signal Processing, 24(2010)：2104-2123

[2] Y. Huang, C. J. Yan, et al. On the difference between empirical mode decomposition and Hilbert vibration decomposition for earthquake motion records[R].15 (World Conference On Earthquake Engineering ) WCEE LISBOA 2012

[3] C. W. Poon, and C. C. Chang. Identification of nonlinear elastic structures using empirical mode decomposition and nonlinear normal modes[J]. Smart Structures and Systems, Vol. 3, No. 2 (2007)

[4] Hui Li, Yuping Zhang, et al. Hilbert-Huang transform and marginal spectrum for detection and diagnosis of localized defects in roller bearings[J]. Journal of Mechanical Science and Technology, 23 (2009)：291-301

[5] Yujun Li, Peter W. Tse, et al. EMD-based fault diagnosis for abnormal clearance between contacting components in a diesel engine. Mechanical Systems and Signal Processing[J].24 (2010)：193-210

[6] Hiralal Khatri, Kenneth Ranney, et al.. Empirical mode decomposition based features for diagnosis and prognostics of systems[J]. Sensors and Electron Devices Directorate. ARL. Army Research Laboratory. Adelphi, MD 20783-1197

[7] Shuyou WU, Weiguo HUANG, et al. extracting power transformer vibration features by a time-scale-frequency analysis method[J]. J. Electromagnetic Analysis & Applications, 2010(2)31-38

[8] 劉慶敏，揚(yáng)午陽，等. 基于經(jīng)驗(yàn)?zāi)B(tài)分解的地震相分析技術(shù)[J].北京：石油地球物理勘探，2010,45(增刊1)：145-149