一種求解信號包絡(luò)曲線端點值的新方法

馬小剛，王永平，杜 歷

（寧夏水利科學(xué)研究院，銀川 750021）

1998年，美國宇航局的華裔科學(xué)家黃鍔（N.E.Huang）博士提出在瞬時頻率概念基礎(chǔ)上用希爾伯特-黃變換（Hilbert-Huang Transform，HHT）處理非線性、非平穩(wěn)信號，包括經(jīng)驗?zāi)Ｊ椒纸猓‥mpirical Mode Decomposition，EMD）和希爾伯特譜分析（Hilbert Spectrum Analysis，HSA）兩個處理過程[1–2]。由于經(jīng)EMD算法分解信號后得到的固有模態(tài)函數(shù)（Intrinsic Mode Function，IMF）分量具有非線性，故該方法非常適合用于分析非線性、非平穩(wěn)信號。然而，該方法采用三次樣條插值來擬合信號數(shù)據(jù)的包絡(luò)曲線，很難實現(xiàn)包絡(luò)曲線準(zhǔn)確到達(dá)兩端邊界處，由此產(chǎn)生分解誤差，這就是所謂的“端點效應(yīng)”。這種誤差一旦產(chǎn)生，會隨著分解層數(shù)的增加而增大，最終使信號的分解結(jié)果無意義[3]。Huang提出用特征波法來解決這一難題，已將該方法申請為美國宇航局專利，但他同時指出該方法僅對一些特定的信號波形有效[4]。隨后，其他學(xué)者對此問題展開研究，提出一些新方法[5–20]，這些新方法對特定信號的處理結(jié)果令人滿意。受鏡像延拓法的啟發(fā)，本文提出一種能有效抑制端點效應(yīng)的新方法，從而實現(xiàn)對非線性、非穩(wěn)定信號的分析和處理。

1 端點效應(yīng)產(chǎn)生機(jī)理

通常端點效應(yīng)的產(chǎn)生是由三次樣條擬合信號數(shù)據(jù)包絡(luò)曲線時引起，這是因為擬合包絡(luò)曲線時需要用到前后兩個相鄰的極值點序列，如此一來就會引發(fā)信號數(shù)據(jù)的極值序列在兩端邊界處的取值問題，即極值序列在邊界兩端處既有可能是極大值點，也有可能是極小值點，還有可能是非極值點。由于是在信號數(shù)據(jù)兩端邊界處取值，因此，也有學(xué)者將其稱之為“邊界效應(yīng)”問題[6–7]，如圖1所示。信號曲線左端邊界處的下包絡(luò)曲線越過上包絡(luò)曲線，出現(xiàn)嚴(yán)重失真的“端點效應(yīng)”問題，而且這種誤差會隨著分解層數(shù)增加不斷增加，將會使得分解結(jié)果誤差過大。由EMD算法的篩分過程可知，第1階IMF分量的準(zhǔn)確與否能嚴(yán)重影響到后續(xù)IMF分量的準(zhǔn)確性。如果第1階IMF分量有較大誤差，那么隨著篩分層數(shù)的增加，誤差會不斷累積傳播，將會導(dǎo)致用EMD算法分解信號后的結(jié)果失去意義。

圖1 端點效應(yīng)

2 最近一次極值法

2.1 最近一次極值法求解原理

最近一次極值鏡像延拓法，可簡稱為“最近一次極值法”。與鏡像延拓法不同，它不是利用原始信號數(shù)據(jù)波形進(jìn)行鏡像延拓，而是先求出原始信號的上下包絡(luò)曲線，然后通過一次極值法直接求解上下包絡(luò)曲線在邊界兩端的極值點，以此達(dá)到包絡(luò)曲線能準(zhǔn)確到達(dá)邊界兩端的目的，從而避免端點效應(yīng)問題的產(chǎn)生。由于該方法是利用原始信號極值序列中離左右兩端最近的一次極值作為拐點和鏡像對稱中心來延拓包絡(luò)曲線，故將其稱之為“最近一次極值鏡像延拓法”，“一次極值”的含義則是指在求包絡(luò)曲線時需對原始信號使用一次極值求解方法才能獲得該信號的極值序列。

設(shè)離散信號

其采樣步長為Δt，X(t)有M個極大值和N個極小值，對應(yīng)的序列下標(biāo)(Im,In)、時間(Tm,Tn)和函數(shù)值(U,V)記為

若要利用最近一次極值法來延拓原始信號的極值序列，一般有如下幾種情況需要討論。

（1）信號序列左端邊界

令左端需要延拓的極值點數(shù)為p，則

Tm(-p+1)、Tm(-p+2)、…、Tm(0)為左端延拓的極大值點所對應(yīng)的時間，xm(-p+1)、xm(-p+2)、…、xm(0)為左端延拓的極大值點所對應(yīng)的函數(shù)值；Tn(-p+1)、Tn(-p+2)、…、Tn(0)為左端延拓的極小值點所對應(yīng)的時間，xn(-p+1)、xn(-p+2)、…、xn(0)為左端延拓的極小值點所對應(yīng)的函數(shù)值。

（2）信號序列右端邊界

令右端需要延拓的極值點數(shù)為q，則為右端延拓的極大值點所對應(yīng)的時間，為右端延拓的極大值點所對應(yīng)的函數(shù)值；為右端延拓的極小值點所對應(yīng)的時間為右端延拓的極小值點所對應(yīng)的函數(shù)值。

綜上所述，最近一次極值法的主要思想和基本原理為：充分利用原始信號局部波形在其邊界附近拐點處的自然走向插補(bǔ)延展極值序列在邊界端點處的取值，旨在確保極值序列能準(zhǔn)確到達(dá)邊界兩端。由于采用的是真實存在的波形拐點來順延極值序列，故延拓部分的波形趨勢能嚴(yán)格遵循信號序列的發(fā)展趨勢，從而避免端點效應(yīng)的產(chǎn)生。

2.2 最近一次極值法算法步驟

原始信號的極值序列到達(dá)兩端邊界處的延拓應(yīng)當(dāng)包括對其上下包絡(luò)曲線左右兩端共計四個極值點的求解，下面以一個非線性、周期性信號x(t)為例，對利用最近一次極值法求解其左端邊界值的具體過程進(jìn)行介紹。

設(shè)原始信號為x(t)：

（1）先求出原始信號x(t)的所有極小值點如圖2中的黑色“+”點所示，由于這些極小值點是第一次求出的，為了以示區(qū)別稱其為一次極小值點

（6）求出水平距離（橫坐標(biāo)值）tm后，以中心對稱軸為鏡像軸找出水平距離（橫坐標(biāo)值）tm在中心對稱軸右側(cè)的鏡像水平距離（橫坐標(biāo)軸）tp，即tp=2tm；

（7）過鏡像水平距離（橫坐標(biāo)值）tp做垂直于橫坐標(biāo)軸并平行于縱坐標(biāo)軸的垂線Lυ，且該垂線與已求出的下包絡(luò)曲線相交于一點P′(tp,xp)，該點即為下包絡(luò)曲線左端點P(t1,xp)關(guān)于中心對稱軸的鏡像點P′(tp,xp)；

（8）再過剛求出的鏡像點P′(tp,xp)做垂直于縱坐標(biāo)軸并平行于橫坐標(biāo)軸的水平線Lh，該水平線與縱坐標(biāo)軸交于一點P(t1,xp)，該點即為所求的信號下包絡(luò)曲線延拓至左端點處的端點值。

需要進(jìn)一步指出的是：由于信號x(t)具有周期性，依據(jù)該信號的周期重復(fù)性，可知圖2中的Pr應(yīng)該就是該信號x(t)下包絡(luò)曲線延拓至左邊界處的精確解，而從圖2上也能直觀地發(fā)現(xiàn)該精確解Pr與用最近一次極值法求出的端點值P恰好是重合的，這也就充分證明了利用最近一次極值法所求信號x(t)的包絡(luò)曲線在邊界處端點值的有效性和準(zhǔn)確性。

圖2 根據(jù)最近一次極值法求端點值

3 仿真算例

對一仿真信號x(t)先用最近一次極值法進(jìn)行延拓處理，然后再進(jìn)行EMD分解，并與其它延拓方法的EMD分解結(jié)果相比較[8–13]，以此來探討最近一次極值法能否有效抑制端點效應(yīng)的產(chǎn)生。

設(shè)仿真信號x(t)由頻率為120 Hz的余弦信號和基頻為20 Hz、調(diào)頻為10 Hz的調(diào)幅調(diào)頻信號疊加而成，具體形式如式（39）所示。

其中，采樣頻率Fs=1 000Hz。原始信號及其組成分量如圖3所示。

下面依次使用偶延拓EMD算法（Even Extension of EMD，EE-EMD）、鏡像延拓EMD算法（Mirror Extension of EMD，ME-EMD）、周期延拓EMD算法（Periodic Extension of EMD，PE-EMD）、最近一次極值EMD算法（The Nearest Simple Extreme Mirror Extension of EMD，NSEME-EMD）分解仿真信號x(t)，分解時間依次為0.172 s、0.219 s、0.344 s、0.547 s，結(jié)果見圖4和圖5。

圖4所示的IMF1對應(yīng)于余弦信號cos(120 πt)分量，IMF2對應(yīng)于調(diào)幅調(diào)頻信號(1+0.2sin(10 πt))·cos(40 πt+0.5sin(20 πt))分量，殘余分量Res可視為EMD算法篩分過程中所產(chǎn)生的分解誤差。

圖5所示的error1表示分解出的IMF1分量與余弦信號cos(120 πt)分量之間的誤差，error2表示分解出的IMF2分量與余弦信號分量之間的誤差。由圖4和圖5可知：EE-EMD算法對端點效應(yīng)的抑制效果較差，ME-EMD算法對端點效應(yīng)的抑制效果稍好，PE-EMD算法對端點效應(yīng)的抑制效果較好，NSEME-EMD算法對端點效應(yīng)的抑制效果最好。

為了準(zhǔn)確評價上述四種EMD算法的分解精度，此處采用基于能量的評價指標(biāo)θi進(jìn)行計算分析[5,9,16,18]

圖3 原始信號及其組成分量

其中：θi為第i個IMF的評價指標(biāo)；RMSi和RMSre為分解得到的第i個IMF和相應(yīng)真實值的有效值。

RMS按下式計算

其中m為信號的采樣點數(shù)；s(j)為信號序列。

通過式（40）可知計算得出的θi≥0 ，如果θi=0，則說明端點效應(yīng)對EMD分解無影響；如果θi＞0，則說明端點效應(yīng)對EMD分解有影響，并且θi值越大影響也越大。

四種EMD算法分解仿真信號的性能評價指標(biāo)計算結(jié)果見表1。

表1 性能評價指標(biāo)計算結(jié)果表/(%)

由表1可知：四種EMD算法的性能指標(biāo)計算結(jié)果從大到小按序排為EE-EMD>ME-EMD>PE-EMD>NSEME-EMD，換句話說，EE-EMD算法在分解復(fù)雜信號x(t)的過程中造成的能量泄露最大，對信號序列x(t)的分解精度最低，ME-EMD算法和PEEMD算法在分解復(fù)雜信號x(t)的過程中造成的能量泄露較小，對信號序列x(t)的分解精度較低，而本文提出的NSEME-EMD算法在分解復(fù)雜信號x(t)的過程中造成的能量泄露最小，對信號序列x(t)的分解精度最高。

圖4 基于四種EMD算法的分解結(jié)果

圖5 根據(jù)四種EMD算法得到的IMF分量的誤差

4 應(yīng)用實例

通過仿真實例驗證了最近一次極值法能有效抑制端點效應(yīng)，故將其應(yīng)用到實際泵站廠房結(jié)構(gòu)振動監(jiān)測信號的分析中去。對寧夏新海水務(wù)有限公司同心分公司三泵站的廠房結(jié)構(gòu)進(jìn)行振動監(jiān)測，從現(xiàn)場采集到的振動信號中選取一條采樣頻率為1 000 Hz、時長為2 s的豎向加速度響應(yīng)信號的時程數(shù)據(jù)序列進(jìn)行分析，如圖6所示。

分別采用Huang提出的傳統(tǒng)EMD方法和本文提出的最近一次極值法對該加速度響應(yīng)信號進(jìn)行分解，具體分解結(jié)果見圖7。

經(jīng)計算，采用傳統(tǒng)EMD方法分解得到前5階各個IMF分量相對應(yīng)的性能評價指標(biāo)θ1～θ5分別為2.576%、3.864%、6.182%、10.510%、18.918%，而采用最近一次極值法分解得到的前5階各個IMF分量相對應(yīng)的性能評價指標(biāo)θ1～θ5分別為1.989%、2.486%、3.232%、4.363%、6.109%，由此可見相比于傳統(tǒng)EMD方法，最近一次極值法對振動信號的分解精度高，能量泄漏少，對端點效應(yīng)的抑制效果好。

圖6 實測豎向加速度響應(yīng)信號時程數(shù)據(jù)序列

圖7 實測豎向加速度響應(yīng)信號分解結(jié)果

5 結(jié)語

（1）通過最近一次極值法得到的邊界端點值能較好匹配原始信號在邊界附近的自然走向，實現(xiàn)原始信號在邊界端點處的光滑過渡，有效降低原始信號極值點在邊界處的波動，從原理上來說不失為一種解決端點效應(yīng)的好方法。

（2）通過仿真算例和應(yīng)用實例對最近一次極值法的分解性能進(jìn)行對比分析，結(jié)果表明該方法能有效解決端點效應(yīng)問題，同時具有簡單易行、適應(yīng)性強(qiáng)、通用性好等特點。