王江萍,常鵬飛
(西安石油大學(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)的固有頻率一直以來(lái)都是國(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)混疊、過(guò)包絡(luò)、欠包絡(luò)等一系列問(wèn)題[7]。
變分模態(tài)分解(Variational Mode Decomposition,VMD)[8]是Dragomiretskiy等人于2014年提出的一種新的自適應(yīng)信號(hào)分解方法,該方法在分解過(guò)程中通過(guò)循環(huán)迭代求取約束變分問(wèn)題的最優(yōu)解來(lái)確定分解得到的固有模態(tài)分量的頻率中心及帶寬,實(shí)現(xiàn)信號(hào)各頻率成分的有效分離。與EMD和LMD相比,VMD方法具有牢固的數(shù)學(xué)理論基礎(chǔ),并且具有較高的運(yùn)算效率及良好的噪聲魯棒性[9]。本文以懸臂梁振動(dòng)加速度響應(yīng)仿真信號(hào)為例,通過(guò)比較仿真結(jié)果和VMD、EMD、EEMD方法的固有頻率識(shí)別結(jié)果,驗(yàn)證了基于VMD的固有頻率識(shí)別方法的有效性。
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)
為求取上述約束變分問(wèn)題的最優(yōu)解,引入增廣Lagrange函數(shù),即
(4)
式中:α為懲罰因子;λ為L(zhǎng)agrange乘子。
利用交替方向乘子算法求取上述增廣Lagrange函數(shù)的鞍點(diǎn),即為式(3)約束變分模型的最優(yōu)解,從而將f分解為K個(gè)窄帶IMF分量。分解過(guò)程如下:
(2)使n=n+1,開始執(zhí)行迭代循環(huán)。
(5)
(6)
(7)
(5)對(duì)于給定的判別精度ε0,若
(8)
則停止迭代,完成分解。否則重復(fù)步驟(2)至(5)。
利用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í)間歷程分析,是用來(lái)確定結(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
利用式(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ù)擬合后的直線來(lái)識(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階固有頻率,見(jiàn)表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
模態(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é)果見(jiàn)表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)的近似。
經(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ù)曲線來(lái)識(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ù)和過(guò)零點(diǎn)的個(gè)數(shù)必須相等或最多相差一個(gè);
(2)任意時(shí)刻,由局部極大值點(diǎn)形成的包絡(luò)線和由局部極小值點(diǎn)形成的包絡(luò)線的均值為零。
對(duì)相同的原始信號(hào)(圖2)利用EMD進(jìn)行分解,其分解過(guò)程如下:
(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)中分離出來(lái),得到
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ù),不能再?gòu)闹刑崛M足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中可見(jiàn),只有0.15~0.409 6 s時(shí)間段各階固有模態(tài)瞬時(shí)頻率相對(duì)穩(wěn)定且都為正值。此時(shí)間段內(nèi)各階瞬時(shí)頻率對(duì)應(yīng)的擬合直線變動(dòng)中心即固有頻率,見(jiàn)表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正交性的說(shuō)明[14]:分解得到的成分并不滿足全局正交性,只是滿足局部正交性。其實(shí)這與分解得到的IMF的性質(zhì)有關(guān),分解得到的相鄰IMF成分雖然是嚴(yán)格按照頻率從高到低產(chǎn)生的,但是其意思并不是說(shuō)低階IMF的頻率一定比高階IMF的高,正確的理解是低階IMF中的某個(gè)局部的頻率比高階IMF中相同局部的頻率要高,這也正好反映了EMD算法局部性強(qiáng)的本質(zhì)所在。所以在本研究中,IMF1的瞬時(shí)頻率的計(jì)算值>IMF2的瞬時(shí)頻率的計(jì)算值>IMF3的瞬時(shí)頻率的計(jì)算值,IMF3的瞬時(shí)頻率的計(jì)算值 EEMD是對(duì)EMD算法的改進(jìn)[15]。通過(guò)在原始信號(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分解過(guò)程如下: (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所示),以及固有頻率(見(jiàn)表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算法本身沒(méi)有突破EMD算法的實(shí)質(zhì),EEMD算法中加入的白噪聲是無(wú)限帶寬信號(hào),會(huì)影響原信號(hào)的頻譜分布。所以EEMD方法同樣有EMD方法的一些缺點(diǎn),包括會(huì)出現(xiàn)負(fù)頻率、瞬時(shí)頻率波動(dòng)較大、可以識(shí)別的階數(shù)少以及IMF瞬時(shí)頻率計(jì)算值的階數(shù)與固有頻率的階數(shù)無(wú)對(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à)值較小。 (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)象等問(wèn)題,從而產(chǎn)生虛假固有頻率,計(jì)算可靠性受到限制。由于其分解得到的瞬時(shí)頻率由高到低,在實(shí)際工程中無(wú)法直接確定機(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)混疊問(wèn)題的分析及改進(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.3.3 EEMD方法識(shí)別固有頻率
4 結(jié) 論