hh中的包絡(luò)線和分段冪函數(shù)法

hh中的包絡(luò)線和分段冪函數(shù)法

本征模函數(shù)方法信號分析方法廣泛應(yīng)用于振動工程中。1998年美國科學家Norden.E.Huang等人提出了一種主要用于分析非平穩(wěn)信號的新信號分析方法——希爾伯特-黃變換(Hilbert-Huangtransform,HHT),其主要創(chuàng)新是本征模函數(shù)(Intrinsicmodefunction,IMF)概念的提出和經(jīng)驗?zāi)B(tài)分解(Empiricalmodefunction,EMD)方法的引入:通過EMD,將信號分解成IMF(一般為有限數(shù)目)的和,對每個IMF分別進行Hilbert變換譜分析就可以獲得有意義的瞬時頻率,從而給出非平穩(wěn)信號中頻率隨時間變化的精確表達。目前這一方法已得到許多應(yīng)用,正吸引越來越多的人對其進行研究。HHT的關(guān)鍵問題是通過信號極值點擬合信號的包絡(luò)線,但目前采用的包絡(luò)線(擬合)算法還存在一些不足,需要進一步研究與完善。1信號極值的估計HHT中的IMF是滿足以下兩個條件的信號:(1)整個信號中,零點數(shù)與極點數(shù)相等或至多相差1;(2)信號上任意一點,由局部極大值點確定的包絡(luò)線和由局部極小值點確定的包絡(luò)線的均值為零。用EMD將信號分解成IMF的和的步驟是:先識別出信號(設(shè)為s(t))的所有極值點,然后用所有極大值點和所有極小值點分別擬合出s(t)的上包絡(luò)線esup(t)和下包絡(luò)線elow(t),滿足elow(t)≤s(t)≤esup(t)(1)elow(t)≤s(t)≤esup(t)(1)記上、下包絡(luò)線的平均值為m(t)m(t)=(esup(t)+elow(t))/2(2)m(t)=(esup(t)+elow(t))/2(2)將s(t)減去m(t)得c(t)=s(t)-m(t)(3)c(t)=s(t)?m(t)(3)將c(t)視為新的s(t),反復重復以上過程就可以篩選出原信號中的第一階IMFc1(t)(由判止準則判斷)。m1(t)=s(t)-c1(t)(4)m1(t)=s(t)?c1(t)(4)將m1(t)視為新的s(t),用同樣的方法可以依次篩選出原信號中的其他IMF。s(t)最終可表示為s(t)=n∑k=1ck(t)+mn(t)(5)s(t)=∑k=1nck(t)+mn(t)(5)其中mn(t)稱為s(t)的余項,代表信號的平均趨勢。對式(5)中的每個IMF應(yīng)用式(6)fik(t)=12πddt[argzk(t)](6)fik(t)=12πddt[argzk(t)](6)就可以分別獲得它們的有物理意義的瞬時頻率,其中上標“i”表示瞬時,zk(t)為分量ck(t)的解析信號,它通過Hilbert變換計算,即zk(t)=ck(t)+jΗ[ck(t)](7)zk(t)=ck(t)+jH[ck(t)](7)其中H[ck(t)]為ck(t)的Hilbert變換,即Η[ck(t)]=ck(t)*1πt=Ρ.V.∫∞-∞ck(t-τ)πτdτ(8)這里P.V.表示柯西主值積分?？梢钥闯?通過信號極值點擬合信號包絡(luò)線的算法影響著EMD的全過程,最終也影響著瞬時頻率的計算,因而是HHT中的關(guān)鍵問題。但如何從理論上嚴格確定HHT中的包絡(luò)線算法是迄今還未解決的問題,目前采用的是經(jīng)驗算法,其中最常用的經(jīng)驗算法是三次樣條插值法。實踐證明,三次樣條插值法的確是一種相對較理想的包絡(luò)線算法,在許多情況下能獲得比較滿意的分析效果。但三次樣條插值法容易出現(xiàn)過沖問題,如圖1(a)中A,B和C處所示?！斑^沖”不僅容易使上、下包絡(luò)線的平均值曲線發(fā)生偏移,甚至使EMD過程所要求的條件式(1)也得不到滿足,如圖1(a)中的B處所示,因此需對HHT中的包絡(luò)線算法進行深入探索和改進。作者在研究中發(fā)現(xiàn),三次樣條插值法容易造成過沖的根本原因是其“柔性”不夠。事實上,對于一個有經(jīng)驗的人(如木工師傅)來說,對一條給定的曲線,可以憑經(jīng)驗畫出其上、下兩條“很好”的包絡(luò)線。例如,對于圖1(a)中A,B和C處所出現(xiàn)的過沖,木工師傅們在該處會設(shè)法使包絡(luò)線更具“柔性”一些。但在算法中,“柔性”和“光滑性”是一對矛盾,因此算法中使包絡(luò)線“更具柔性”意味著要犧牲“光滑性”,而三次樣條插值曲線是二階光滑的(即二階連續(xù)可導),因此作者想到用只有一階光滑性(即插值曲線一階連續(xù)可導)的典型插值法——Akima插值法進行實驗。Akima插值法依據(jù)如下原理:對于一個有經(jīng)驗的人來說,當它在考慮由一系列數(shù)據(jù)點所確定的曲線時,往往只要注意所需的一段曲線附近的幾個數(shù)據(jù)點就行,因為這幾個數(shù)據(jù)點的位置完全可以描繪出這段曲線的形狀。所以用Akima插值法來作為HHT中的包絡(luò)線算法具有一定的可行性。但作者經(jīng)過實驗發(fā)現(xiàn),雖然用Akima插值法是可行的,但有時所擬合的曲線又顯得“過于柔軟”,甚至出現(xiàn)明顯的“折點”,如圖1(b)中的A,B和C處所示。2分段物理法的原理在HHT中的包絡(luò)線算法問題上,三次樣條插值法柔性不夠,而Akima插值法“過于柔軟”(即光滑性不夠)。下面先介紹曲線擬合中的拋物線參數(shù)樣條插值法。對插值點P1,P2,…,Pn進行拋物線參數(shù)樣條插值的步驟是:將所有插值點中任意3個相鄰點Pi-1,Pi,Pi+1用拋物線參數(shù)樣條插值,將Pi,Pi+1,Pi+2也用拋物線參數(shù)樣條插值。設(shè)所得拋物線分別為Si(t),Si+1(t),則Pi與Pi+1之間的曲線用Si(t),Si+1(t)重新插值即得拋物線參數(shù)樣條曲線。Pi與Pi+1之間的最終插值曲線的參數(shù)形式為S(t)=(1-2t)Si(t+0.5)+2tSi+1(t)(9)或S(t)=(-t+4t2-4t3)Ρi-1+(t-10t2+12t3)Ρi+(t+8t2-12t3)Ρi+1+(-2t2+4t3)Ρi+2(10)其中0≤t≤0.5。已經(jīng)證明這樣的曲線是一階光滑曲線。拋物線參數(shù)樣條插值法是一種常用的曲線擬合算法,如MicrosoftWord中的曲線繪圖就是使用的該算法。用拋物線參數(shù)樣條法所擬合出的曲線既具有很好的柔性,又具有很好的光滑性,因此具有HHT中包絡(luò)線算法所希望的性質(zhì),但拋物線參數(shù)樣條插值法不能保證所插值曲線是單值的,如圖2中A,B,C,D為插值點,插值曲線在x0處顯然是多值的。作者利用拋物線參數(shù)樣條插值法的原理提出一種新的單值插值法——分段冪函數(shù)法。分段光滑定理:設(shè)P1(x1,y1),P2(x2,y2),P3(x3,y3),P4(x4,y4)為平面內(nèi)4點,且x1≤x2≤x3≤x4,S2(x),S3(x)分別為經(jīng)P1,P2,P3和P2,P3,P4所插值出的單值光滑曲線(即其一階連續(xù)可導),則曲線S(x)=x3-xx3-x2S2(x)+x-x2x3-x2S3(x)(11)在區(qū)間(x1,x4)內(nèi)一階光滑。證明:顯然只需證明S(x)在點P2,P3一階連續(xù)可導即可。由式(11)得S′(x)=1x3-x2(-S2(x)+(x3-x)S′2(x)+S3(x)+(x-x2)S′3(x))所以S′(x2)=1x3-x2(-S2(x2)+(x3-x2)S′2(x2)+S3(x2))=1x3-x2(-y2+(x3-x2)S′2(x2)+y2)=S′2(x2)故S(x)在點P2一階連續(xù)可導,同理可證S(x)在點P3一階連續(xù)可導。分段冪函數(shù)法的步驟與拋物線參數(shù)樣條插值法類似,只是分段冪函數(shù)法采用一種冪函數(shù)的方法對插值點P1,P2,…,Pn中的相鄰點Pi-1,Pi,Pi+1和Pi,Pi+1,Pi+2分別插值,再依據(jù)式(11)對兩段曲線進行“拼接”。在冪函數(shù)法中對任意3個插值點P1,P2,P3插值的方法如下:如圖3(a)所示,設(shè)過點P2且與P1P3平行的直線分別交直線x=x1和x=x3于點A和B,先構(gòu)造冪函數(shù)y=xβ(β>1且β的取值使y=xβ的定義域為(-∞,+∞),值域為[0,+∞)),再構(gòu)造函數(shù)y={→A′Ρ′1→A′Cxβx≤0→B′Ρ′3→B′Dxβx≥0={→AΡ1(x1-x2)βxβx≤0→BΡ3→(x3-x2)βxβx≥0={(xx1-x2)β→AΡ1x≤0(xx3-x2)β→AΡ1x≥0(12)其中A′P′1=AP1,B′P′3=BP3,C,D分別為曲線y=xβ與直線x=x1-x2和x=x3-x2的交點,如圖3(b)所示。設(shè)過點P2且與P1P3平行的直線的方程為y=g(x),式(12)最右邊的函數(shù)為f(x),則對插值點P1,P2,P3采用式(13)進行插值y=f(x-x2)+g(x)(13)即S2(x)={(x-x2x1-x2)β→AΡ1+y3-y1x3-x1(x-x2)+y2x≤x2(x-x2x3-x2)β→AΡ1+y3-y1x3-x1(x-x2)+y2x≥x2(14)其中→AΡ1=(x3-x2)y1-(x2-x1)y3x3-x1。容易證明,由式(14)確定的S2(x)為單值光滑曲線。分段冪函數(shù)法插值中參數(shù)β的取值很關(guān)鍵,它可以是整數(shù)和小數(shù),可以有無窮多種取值,但作者經(jīng)過大量實驗發(fā)現(xiàn)當β取2.5時效果較為理想,這大概是因為β略低于三次樣條插值中各分段函數(shù)的次數(shù)(3次),因而更具柔性。3包絡(luò)線算法hhd分析圖4給出了用分段冪函數(shù)法對圖1(a,b)中同一信號s所擬合出的包絡(luò)線,從圖中A,B,C,D處可以看出,新包絡(luò)線合算法使三次樣條插值法的過沖問題和Akima插值法的“折點”問題均得到了較明顯的改善。圖5(a～c)分別給出了用三次樣條插值法、分段冪函數(shù)法和Akima插值法作為包絡(luò)線算法時對圖1中信號s進行EMD分解的結(jié)果,圖6則給出了它們各自前兩階IMF分量的瞬時頻率(由于第三階IMF分量不是主要分量,故這里不對其瞬時頻率進行分析)。從圖5,6可以看出,用3種包絡(luò)線算法所分析出的結(jié)果在總體上很一致,表明3種算法總體上都比較合理;但在細節(jié)部分則存在明顯差異。從圖5中還難以說明這種差異代表了孰優(yōu)孰壞,但從圖6(a,c)可以看出,采用三次樣條插值法時所得到的第一階IMF分量的瞬時頻率在圖中A和B處出現(xiàn)了明顯跳躍,采用Akima插值法時所得到的第二階IMF分量的瞬時頻率在圖中C和D處出現(xiàn)了較大波動,而從圖1中卻看不出原信號具有這種急劇變化的頻率(實際上,圖1中的信號是用三次樣條曲線所擬合出的曲線,理論上也不存在急劇變化的頻率),這說明,用三次樣條插值法或Akima插值法作為包絡(luò)線算法進行HHT分析時會出現(xiàn)明顯的虛假頻率,其原因正是由于三次樣條插值法存在過沖,而Akima插值法存在“折點