王江萍，常鵬飛

(西安石油大學(xué) 機(jī)械工程學(xué)院，陜西 西安 710065)

引 言

機(jī)械系統(tǒng)的振動(dòng)信號(hào)蘊(yùn)藏著豐富的狀態(tài)信息，對(duì)機(jī)械系統(tǒng)振動(dòng)模態(tài)進(jìn)行分析，獲得系統(tǒng)的模態(tài)參數(shù)，是認(rèn)識(shí)機(jī)械系統(tǒng)動(dòng)態(tài)特性設(shè)計(jì)合理與否，運(yùn)行狀態(tài)正常與否的關(guān)鍵。在機(jī)械系統(tǒng)的振動(dòng)分析中，利用測(cè)得振動(dòng)信號(hào)對(duì)機(jī)器零部件的固有頻率進(jìn)行識(shí)別，進(jìn)而對(duì)系統(tǒng)固有頻率進(jìn)行分析，以確定故障的部位的故障診斷方法得到了廣泛的應(yīng)用[1]。而如何從振動(dòng)信號(hào)中得到機(jī)械系統(tǒng)的固有頻率一直以來都是國(guó)內(nèi)外的一個(gè)研究熱點(diǎn)。

實(shí)際工程中得到的機(jī)械振動(dòng)響應(yīng)信號(hào)一般是非平穩(wěn)信號(hào)，時(shí)域識(shí)別法和頻域識(shí)別法一般只能處理平穩(wěn)振動(dòng)信號(hào)，要識(shí)別非平穩(wěn)信號(hào)，需要采用時(shí)頻分析方法。典型的時(shí)頻分析方法有短時(shí)傅里葉變換[2]、小波變換[3]、經(jīng)驗(yàn)?zāi)B(tài)分解(Empirical Mode Decomposition，EMD)[4]、局部均值分解 (Local Mean Decomposition，LMD)[5]、集合經(jīng)驗(yàn)?zāi)B(tài)分解(Ensemble Empirical Mode Decomposition，EEMD)[6]等。短時(shí)傅里葉變換會(huì)產(chǎn)生譜泄漏，使所估計(jì)的譜值發(fā)生偏差；小波變換需要事先確定小波基，缺乏自適應(yīng)性；EMD和LMD方法存在模態(tài)混疊、過包絡(luò)、欠包絡(luò)等一系列問題[7]。

變分模態(tài)分解(Variational Mode Decomposition,VMD)[8]是Dragomiretskiy等人于2014年提出的一種新的自適應(yīng)信號(hào)分解方法，該方法在分解過程中通過循環(huán)迭代求取約束變分問題的最優(yōu)解來確定分解得到的固有模態(tài)分量的頻率中心及帶寬，實(shí)現(xiàn)信號(hào)各頻率成分的有效分離。與EMD和LMD相比，VMD方法具有牢固的數(shù)學(xué)理論基礎(chǔ)，并且具有較高的運(yùn)算效率及良好的噪聲魯棒性[9]。本文以懸臂梁振動(dòng)加速度響應(yīng)仿真信號(hào)為例，通過比較仿真結(jié)果和VMD、EMD、EEMD方法的固有頻率識(shí)別結(jié)果，驗(yàn)證了基于VMD的固有頻率識(shí)別方法的有效性。

1 變分模態(tài)分解原理

VMD將固有模態(tài)函數(shù)(IMF)定義為一個(gè)調(diào)幅-調(diào)頻信號(hào)，其表達(dá)式為：

uk(t)=Ak(t)cos(φk(t))。

(1)

式中：t為時(shí)間，Ak(t)為uk(t)的瞬時(shí)幅值，φk(t)為瞬時(shí)相位函數(shù)，ωk(t)=dφk(t)，ωk(t)為uk(t)的瞬時(shí)頻率。

假設(shè)多分量信號(hào)f由K個(gè)有限帶寬的IMF分量uk組成，且各IMF的中心頻率為ωk，VMD方法建立的約束變分模型為

(2)

uk滿足

(3)

為求取上述約束變分問題的最優(yōu)解，引入增廣Lagrange函數(shù)，即

(4)

式中：α為懲罰因子；λ為L(zhǎng)agrange乘子。

利用交替方向乘子算法求取上述增廣Lagrange函數(shù)的鞍點(diǎn)，即為式(3)約束變分模型的最優(yōu)解，從而將f分解為K個(gè)窄帶IMF分量。分解過程如下：

(2)使n=n+1，開始執(zhí)行迭代循環(huán)。

(5)

(6)

(7)

(5)對(duì)于給定的判別精度ε0，若

(8)

則停止迭代，完成分解。否則重復(fù)步驟(2)至(5)。

2 VMD方法的固有頻率識(shí)別

2.1 有限元建模與瞬態(tài)動(dòng)力學(xué)分析

利用ANSYS軟件建立如圖1所示的懸臂梁模型，設(shè)置為長(zhǎng)1 m，寬0.1 m，高0.1 m，彈性模量為2×1011Pa，密度為7 800 kg/m3，泊松比為0.3。網(wǎng)格邊長(zhǎng)為0.02 m，Z=0處為約束端。

瞬態(tài)動(dòng)力學(xué)分析，也叫時(shí)間歷程分析，是用來確定結(jié)構(gòu)在隨時(shí)間變化的載荷作用下的結(jié)構(gòu)動(dòng)力響應(yīng)的方法。瞬態(tài)動(dòng)力學(xué)分析可以真實(shí)地模擬結(jié)構(gòu)所受載荷的情況，對(duì)研究系統(tǒng)的動(dòng)力學(xué)特性有著重要的意義。利用ANSYS對(duì)圖1模型進(jìn)行瞬態(tài)動(dòng)力學(xué)分析，施加階躍載荷，集中力為-10 N，方向沿Y，分析時(shí)間為0.409 6 s，得到4 096個(gè)振動(dòng)加速度與時(shí)間相對(duì)應(yīng)的數(shù)據(jù)，如圖2所示。

圖1 懸臂梁有限元模型Fig.1 Finite element model of cantilever beam

圖2 原始信號(hào)Fig.2 Original signal

2.2 VMD分解與希爾伯特變換確定固有頻率

利用式(5)～(8)對(duì)原始信號(hào)進(jìn)行VMD分解，得到固有模態(tài)函數(shù)，固有模態(tài)函數(shù)中包含固有頻率的信息，對(duì)各階固有模態(tài)函數(shù)進(jìn)行希爾伯特變換，得到對(duì)應(yīng)的瞬時(shí)頻率。對(duì)瞬時(shí)頻率曲線進(jìn)行線性擬合，可根據(jù)擬合后的直線來識(shí)別固有頻率[13]，即可近似認(rèn)為擬合直線的變動(dòng)中心是系統(tǒng)的固有頻率。

對(duì)IMF即u(t)作希爾伯特變換：

(9)

其中，PV代表柯西主值。

構(gòu)造解析信號(hào)z(t)：

z(t)=u(t)+jH[u(t)]=a(t)ejΦ(T)。

(10)

其中，幅值函數(shù)為

(11)

相位函數(shù)為

(12)

瞬時(shí)頻率函數(shù)為

(13)

VMD算法需預(yù)先設(shè)定分解模態(tài)個(gè)數(shù)，本文采用觀察各模態(tài)中心頻率的方法[10]，利用式(10)～(13)得到不同分解階數(shù)時(shí)的前4階固有頻率，見表1，可以看出當(dāng)分解階數(shù)為11階時(shí)前4階固有頻率不再變化，本文取模態(tài)個(gè)數(shù)K=11，根據(jù)文獻(xiàn)[10-12]當(dāng)α取默認(rèn)值2 000時(shí)可滿足大多數(shù)工況實(shí)際需求。分解得到11階固有模態(tài)函數(shù)，圖3所示為前4階。并且得到對(duì)應(yīng)的瞬時(shí)頻率，如圖4所示。

表1 VMD方法識(shí)別固有頻率Tab.1 Identification of natural frequency based on VMD

圖3 VMD 固有模態(tài)函數(shù)Fig.3 VMD intrinsic mode function

圖4 VMD瞬時(shí)頻率Fig.4 VMD instantaneous frequency

3 基于ANSYS、EMD、EEMD方法的固有頻率識(shí)別

3.1 ANSYS模態(tài)分析

模態(tài)分析用于確定設(shè)計(jì)結(jié)構(gòu)或機(jī)器部件的振動(dòng)特性，即結(jié)構(gòu)的固有頻率和振型，它們是承受動(dòng)態(tài)載荷結(jié)構(gòu)設(shè)計(jì)中的重要參數(shù)，常用有限元分析法仿真計(jì)算結(jié)構(gòu)的固有頻率。對(duì)圖1所示模型利用有限元分析軟件ANSYS進(jìn)行模態(tài)分析，采用的分析方法是分塊蘭索斯法，Z=0處為固定端，可得到其前四階固有頻率，結(jié)果見表2。

表2 ANSYS模態(tài)分析計(jì)算固有頻率Tab.2 Calculation of natural frequency based on ANSYS modal analysis

在求取一些模型較為簡(jiǎn)單的機(jī)械結(jié)構(gòu)時(shí)，利用ANSYS軟件既快捷又準(zhǔn)確，但在實(shí)際工程中機(jī)械系統(tǒng)的模型往往比較復(fù)雜，難以建立，其分析結(jié)果受網(wǎng)格劃分密集程度的影響，模型的工況約束與實(shí)際的工況約束可能不同，所以其分析結(jié)果往往是實(shí)際結(jié)構(gòu)振動(dòng)模態(tài)的近似。

3.2 EMD方法識(shí)別固有頻率

經(jīng)驗(yàn)?zāi)B(tài)分解(EMD)固有頻率識(shí)別方法的基本思想與VMD方法有相似之處，首先對(duì)原始信號(hào)進(jìn)行經(jīng)驗(yàn)?zāi)B(tài)分解，獲得IMF，對(duì)其進(jìn)行希爾伯特變換，得到對(duì)應(yīng)的瞬時(shí)頻率曲線，最后根據(jù)曲線來識(shí)別固有頻率。

EMD方法假設(shè)任何信號(hào)都是由不同的固有模態(tài)函數(shù)組成，任意兩個(gè)模態(tài)函數(shù)之間是相互獨(dú)立的，這樣任何一個(gè)信號(hào)就可以分解為有限個(gè)固有模態(tài)函數(shù)之和，其中IMF都應(yīng)滿足以下條件:

(1)整個(gè)數(shù)據(jù)段內(nèi)，極值點(diǎn)的個(gè)數(shù)和過零點(diǎn)的個(gè)數(shù)必須相等或最多相差一個(gè);

(2)任意時(shí)刻，由局部極大值點(diǎn)形成的包絡(luò)線和由局部極小值點(diǎn)形成的包絡(luò)線的均值為零。

對(duì)相同的原始信號(hào)(圖2)利用EMD進(jìn)行分解，其分解過程如下:

(1)確定信號(hào)f(t)中所有的局部極值點(diǎn)，利用三次樣條插值函數(shù)擬合形成信號(hào)f(t)的上、下包絡(luò)線；

(2)將上包絡(luò)線和下包絡(luò)線的均值記作m1，求得：

h1=f(t)-m1。

(14)

(3)判斷h1是否為IMF;若h1不滿足IMF的條件，把h1作為原始數(shù)據(jù)重復(fù)步驟(1)直至滿足IMF的條件。記c1為信號(hào)f(t)第1個(gè)滿足IMF條件的分量，

c1=h1。

(15)

(4)將c1從f(t)中分離出來，得到

r1=f(t)-c1。

(16)

(5)將r1作為原始數(shù)據(jù)重復(fù)步驟(1)～(4)，得到信號(hào)f(t)的n個(gè)滿足IMF條件的分量。所以有

(17)

當(dāng)rn成為一個(gè)單調(diào)函數(shù)，不能再從中提取滿足IMF條件的分量時(shí)，循環(huán)結(jié)束。則信號(hào)f(t)可表示為

(18)

其中rn稱為殘余函數(shù)。

對(duì)原始信號(hào)進(jìn)行EMD分解可得到10階固有模態(tài)。圖5所示為前4階。

圖5 EMD固有模態(tài)函數(shù)Fig.5 EMD intrinsic mode function

重復(fù)公式(9)～(13)得到前4階固有模態(tài)瞬時(shí)頻率，如圖6所示。

機(jī)械系統(tǒng)的固有頻率不隨時(shí)間而變化，且頻率不會(huì)為負(fù)。由于計(jì)算方法的缺陷，從圖6中可見，只有0.15～0.409 6 s時(shí)間段各階固有模態(tài)瞬時(shí)頻率相對(duì)穩(wěn)定且都為正值。此時(shí)間段內(nèi)各階瞬時(shí)頻率對(duì)應(yīng)的擬合直線變動(dòng)中心即固有頻率，見表3。

其中IMF1的瞬時(shí)頻率計(jì)算值對(duì)應(yīng)于懸臂梁的第二階固有頻率，IMF2的瞬時(shí)頻率計(jì)算值對(duì)應(yīng)于懸臂梁的第一階固有頻率，但I(xiàn)MF3及IMF4的瞬時(shí)頻率得到的計(jì)算結(jié)果明顯有誤。

表3 EMD方法識(shí)別固有頻率Tab.3 Identification of natural frequency based on EMD

圖6 EMD瞬時(shí)頻率Fig.6 EMD instantaneous frequency

上述計(jì)算結(jié)果，符合Huang關(guān)于EMD正交性的說明[14]：分解得到的成分并不滿足全局正交性，只是滿足局部正交性。其實(shí)這與分解得到的IMF的性質(zhì)有關(guān)，分解得到的相鄰IMF成分雖然是嚴(yán)格按照頻率從高到低產(chǎn)生的，但是其意思并不是說低階IMF的頻率一定比高階IMF的高，正確的理解是低階IMF中的某個(gè)局部的頻率比高階IMF中相同局部的頻率要高，這也正好反映了EMD算法局部性強(qiáng)的本質(zhì)所在。所以在本研究中，IMF1的瞬時(shí)頻率的計(jì)算值>IMF2的瞬時(shí)頻率的計(jì)算值>IMF3的瞬時(shí)頻率的計(jì)算值，IMF3的瞬時(shí)頻率的計(jì)算值

3.3 EEMD方法識(shí)別固有頻率

EEMD是對(duì)EMD算法的改進(jìn)[15]。通過在原始信號(hào)當(dāng)中添加時(shí)域均值為零、頻域功率譜均勻分布的白噪聲，使信號(hào)具有相對(duì)均勻的極值點(diǎn)，保持了信號(hào)在時(shí)間尺度上的連續(xù)性，從而避免了EMD 分解中會(huì)出現(xiàn)的模態(tài)混疊現(xiàn)象[6]。其固有頻率識(shí)別的基本思想與EMD方法相同。

EEMD分解過程如下:

(1)為待分析信號(hào)f(t)中加入正態(tài)分布的白噪 聲，得到混合信號(hào)

fi(t)=f(t)+ni(t) (i=1,2,3，…，M) 。

(19)

(2)采用EMD分解方法對(duì)混合信號(hào)進(jìn)行分析，可以得到如下相應(yīng)的K個(gè)IMF分量，記為cij(t)(j=1,2,3…，K)，殘余項(xiàng)記為ri(t)，其中cij(t)表示第i次加入高斯白噪聲后分解得到的第j個(gè)IMF分量。

(3)將每次分解得到的IMF分量進(jìn)行集合平均，由于加入的白噪聲不相關(guān)，其統(tǒng)計(jì)均值為0，所以最后得到的IMF分量為

(20)

對(duì)原始信號(hào)進(jìn)行EEMD分解，利用文獻(xiàn)[16]的方法確定白噪聲的幅值系數(shù)取0.1，集合平均次數(shù)M=100?？傻玫?2階固有模態(tài)，圖7所示為前4階。利用相同的方法對(duì)IMF進(jìn)行分析，得到前4階固有模態(tài)瞬時(shí)頻率(如圖8所示)，以及固有頻率(見表4)，其中IMF4的瞬時(shí)頻率計(jì)算值對(duì)應(yīng)于懸臂梁的第一階固有頻率，IMF3的瞬時(shí)頻率計(jì)算值對(duì)應(yīng)于懸臂梁的第二階固有頻率。

圖7 EEMD固有模態(tài)函數(shù)Fig.7 EEMD intrinsic mode function

圖8 EEMD瞬時(shí)頻率Fig.8 EEMD instantaneous frequency

表4 EEMD方法識(shí)別固有頻率Tab.4 Identification of natural frequency based on EEMD

由于EEMD算法本身沒有突破EMD算法的實(shí)質(zhì)，EEMD算法中加入的白噪聲是無限帶寬信號(hào)，會(huì)影響原信號(hào)的頻譜分布。所以EEMD方法同樣有EMD方法的一些缺點(diǎn)，包括會(huì)出現(xiàn)負(fù)頻率、瞬時(shí)頻率波動(dòng)較大、可以識(shí)別的階數(shù)少以及IMF瞬時(shí)頻率計(jì)算值的階數(shù)與固有頻率的階數(shù)無對(duì)應(yīng)關(guān)系等。

表5是各種識(shí)別方法的結(jié)果以及相較于ANSYS模態(tài)分析結(jié)果的誤差。

表5 固有頻率識(shí)別結(jié)果及誤差Tab.5 Results and errors of identification of natural frequency using diffirent methods

由表5可以看出，基于VMD的識(shí)別方法得到的前四階固有頻率與ANSYS軟件得到的結(jié)果接近，有較高的實(shí)際利用價(jià)值。而EMD與EEMD得到的結(jié)果雖然誤差也較小但只能識(shí)別前兩階，實(shí)際利用價(jià)值較小。

4 結(jié) 論

(1)利用有限元分析求取結(jié)構(gòu)固有頻率，受模型建立、有限單元網(wǎng)格劃分密集程度、工況約束條件設(shè)定等因素的影響，其分析結(jié)果往往是實(shí)際結(jié)構(gòu)振動(dòng)模態(tài)的近似。

(2)VMD方法充分利用結(jié)構(gòu)工作中的振動(dòng)信號(hào)，反映工程實(shí)際，計(jì)算結(jié)果較準(zhǔn)確。VMD方法可以識(shí)別前4階固有頻率，與ANSYS的仿真結(jié)果基本一致，而且IMF瞬時(shí)頻率計(jì)算值的階數(shù)與固有頻率的階數(shù)一一對(duì)應(yīng)，有較大的實(shí)際利用價(jià)值。

(3)EMD與EEMD方法雖然同樣利用結(jié)構(gòu)工作中的振動(dòng)信號(hào)，但在識(shí)別固有頻率尤其是高階固有頻率時(shí)，存在瞬時(shí)頻率波動(dòng)較大、會(huì)出現(xiàn)負(fù)值、產(chǎn)生模態(tài)混疊現(xiàn)象等問題，從而產(chǎn)生虛假固有頻率，計(jì)算可靠性受到限制。由于其分解得到的瞬時(shí)頻率由高到低，在實(shí)際工程中無法直接確定機(jī)械系統(tǒng)的固有頻率所對(duì)應(yīng)的IMF的瞬時(shí)頻率計(jì)算值，所以其實(shí)際利用價(jià)值較小。

參 考 文 獻(xiàn)：

[1] 陳雪峰，李兵，胡橋，等.基于小波有限元的裂紋故障診斷[J].西安交通大學(xué)學(xué)報(bào)，2004,38(3)：295-298.

CHEN Xuefeng,LI Bing,HU Qiao,et al.Crack fault diagnosis based on wavelet finite elements[J].Journal of Xi'an Jiaotong University,2004,38(3):295-298.

[2] 徐明遠(yuǎn)，劉增力.MATLAB仿真在信號(hào)處理中的應(yīng)用[M].西安:西安電子科技大學(xué)出版社，2007.

[3] 楊福生.小波變換的工程分析與應(yīng)用[M].北京：科學(xué)出版社，1999.

[4] 莫平杰，楊世錫，曹沖鋒.振動(dòng)模態(tài)固有頻率和阻尼比的EMD識(shí)別方法[J].機(jī)電工程，2011,28(4)：392-396,428.

MO Pingjie,YANG Shixi,CAO Chongfeng.EMD method on identifying natural frequency and damping ratio of vibration mode[J].Mechanical & Electrical Engineering Magazine,2011,28(4):392-396,428.

[5] 程軍圣，朱文峰，李寶慶.基于LMD的模態(tài)參數(shù)識(shí)別方法[J].計(jì)算機(jī)工程與應(yīng)用，2014,51(14)：214-218.

CHENG Junsheng,ZHU Wenfeng,LI Baoqing.Method of identifying analog parameters based on LMD[J].Computer Engineering and Applications,2014,51(14):214-218.

[6] WU Z H,HUANG N E.Ensemble empirical mode decomposition:a noise assisted data analysis method[J].Advances in Adaptive Data Analysis,2009,1(1):1-41.

[7] 程軍圣，鄭近德，楊宇.一種新的非平穩(wěn)信號(hào)分析方法——局部特征尺度分解法[J].振動(dòng)工程學(xué)報(bào)，2012，25(2):215-220.

CHENG Junsheng,ZHENG Jinde,YANG Yu.A nonstationary signal analysis approach:the local characteristic-scale decomposition method[J].Journal of Vibration Engineering,2012,25(2):215-220.

[8] DRAUOMIRETSKIY K,ZOSSO D.Variational mode decomposition[J].IEEE Transactions on Signal Proccssing,2014，62(3):531-544.

[9] 唐貴基，王曉龍.變分模態(tài)分解方法及其在滾動(dòng)軸承早期故障診斷中的應(yīng)用[J]振動(dòng)工程學(xué)報(bào),2016，29(4)：638-647.

TANG Guiji,WANG Xiaolong.Variational mode decomposition method and its application on incipient fault diagnosis of rolling bearing[J].Journal of Vibration Engineering,2016，29(4)：638-647.

[10] 劉長(zhǎng)良，武英杰，甄成剛.基于變分模態(tài)分解和模糊C均值聚類的滾動(dòng)軸承故障診斷[J].中國(guó)電機(jī)工程學(xué)報(bào)，2015,35(13)：3358-3365.

LIU Changliang,WU Yingjie,ZHEN Chenggang.Rolling bearing fault diagnosis based on variational mode decomposition and fuzzy C means clustering[J].Proceedings of the CSEE,2015,35(13):3358-3365.

[11] 岳應(yīng)娟，孫鋼，蔡艷平，等.變分模態(tài)分解在軸承故障診斷中的應(yīng)用[J].軸承，2016(8):50-54,65.

YUE Yingjuan，SUN Gang，CAI Yanping，et al.Application of variational mode decomposition in fault diagnosis for bearings[J].Bearing,2016 (8):50-54,65.

[12] 范磊，衛(wèi)志農(nóng)，李慧杰,等.基于變分模態(tài)分解和蝙蝠算法——相關(guān)向量機(jī)的短期風(fēng)速區(qū)間預(yù)測(cè)[J].電力自動(dòng)化設(shè)備，2017,37(1):93-100.

FAN Lei，WEI Zhinong，LI Huijie，et al.Short-term wind speed interval prediction based on VMD and BA-RVM algorithm[J].Electric Power Automation Equipment,2017,37(1):93-100.

[13] TSUKUNE H,TSUKAMOTO M,MATSUSHITA T,et al.Modular manufacturing[J].Journal of Intelligent Manufacturing,1993,4(2):163-181.

[14] 張健.關(guān)于EMD算法中模態(tài)混疊問題的分析及改進(jìn)方法的研究[D].合肥：中國(guó)科學(xué)技術(shù)大學(xué)，2008.

[15] 苗晟，王威廉，姚紹文.Hilbert-Huang變換發(fā)展歷程及其應(yīng)用[J].電子測(cè)量與儀器學(xué)報(bào)，2014,28(8):812-818.

MIAO Sheng,WANG Weilian,YAO Shaowen.Historic development of HTT and its applications[J].Journal of Electronic Measurement and Instrument，2014,28(8):812-818.

[16] 王鳳利，邢輝，邱赤東,等.基于改進(jìn)自適應(yīng) EEMD 的柴油機(jī)氣缸磨損診斷[J].內(nèi)燃機(jī)學(xué)報(bào),2017,35(1):89-95.

Wang Fengli，Xing Hui，Qiu Chidong，et al.Wear diagnosis of diesel cylinder liner via improved adaptive EEMD method[J].Transactions of CSICE，2017,35(1):89-95.