一種基于本征波匹配的EMD邊界處理方法

鐘佑明，趙 強(qiáng)，周建庭

（1.重慶交通大學(xué) 信息科學(xué)與工程學(xué)院，重慶 400074；2.重慶交通大學(xué) 土木工程學(xué)院，重慶 400074）

希爾伯特－黃變換（Hilbert－Huang Transform，簡稱HHT）是上世紀(jì)末由Huang等人［1－2］為更精準(zhǔn)地分析非線性非平穩(wěn)信號(hào)而提出的一種信號(hào)分析方法。這一方法已被許多學(xué)者應(yīng)用于各種工程和生活領(lǐng)域，取得了較為理想的效果。HHT包含若干關(guān)鍵技術(shù)，EMD（Empirical Mode Decomposition，即經(jīng)驗(yàn)?zāi)B(tài)分解）邊界處理就是其中之一。目前已經(jīng)有不少EMD邊界處理的方法，如特征波法［1］、波形匹配預(yù)測法［3－4］、神經(jīng)網(wǎng)絡(luò)分析方法［5］、鏡像閉合延拓法［6］、多項(xiàng)式擬合法［7］、包絡(luò)極值延拓法［8］、相似極值延拓法等［9］，雖然這些方法都有一定效果，但同時(shí)又有明顯的理論缺陷，如特征波法、鏡像閉合延拓法、包絡(luò)極值延拓法等都僅利用了與邊界相鄰的少數(shù)一些內(nèi)部波來預(yù)測信號(hào)邊界的包絡(luò)或均值，使得它們的應(yīng)用效果不太理想。本文第一作者在文獻(xiàn)［10］中曾提出了一種改進(jìn)的波形匹配預(yù)測法——本征波匹配預(yù)測法。但由于當(dāng)時(shí)對(duì)該方法的有效性還有些懷疑，尤其是還不知道如何合理構(gòu)造仿真信號(hào)來驗(yàn)證該方法，因此文獻(xiàn)［10］只是提出了初步算法而沒有進(jìn)行深入的算例驗(yàn)證。本文最終彌補(bǔ)了這一重要不足，使文獻(xiàn)［10］的研究得到較大程度地深化和完善。

1 EMD及其邊界處理問題

EMD的主要內(nèi)容是通過篩選將信號(hào)分解為有限個(gè)固有模態(tài)函數(shù)（IMF）和趨勢項(xiàng)之和，其中IMF被認(rèn)為是構(gòu)成非線性非平穩(wěn)信號(hào)的基本信號(hào)，它具有唯一的瞬時(shí)頻率，而趨勢項(xiàng)反映信號(hào)變化的趨勢［1］。所謂IMF，是指滿足如下兩個(gè)條件的信號(hào)：

① 信號(hào)中，零點(diǎn)數(shù)與極點(diǎn)數(shù)相等或至多相差1；

② 信號(hào)上任意一點(diǎn)，由局部極大值點(diǎn)確定的包絡(luò)線和由局部極小值點(diǎn)確定的包絡(luò)線的均值都為零。

EMD的過程簡介如下：對(duì)于給定的實(shí)信號(hào)s（t），例如圖1（a）所示信號(hào)，首先確定出其上的所有極大值點(diǎn)和極小值點(diǎn)，然后將所有極大值點(diǎn)和所有極小值點(diǎn)分別用一條光滑曲線（如三次樣條曲線）聯(lián)接起來，使兩條曲線間包含所有的原信號(hào)，如圖1（b）虛線所示。將這兩條曲線分別作為s（t）的上、下包絡(luò)線，計(jì)算出它們的平均值曲線m1（t），如圖1（b）點(diǎn)劃線所示。用s（t）減去m1（t）得到：

h1（t）的波形如圖1（c）所示。以上過程稱為“篩選”。理想情況下，所篩選出的h1（t）應(yīng)該是一個(gè)IMF，因?yàn)閔1（t）的構(gòu)造過程就是使它滿足IMF的條件。但實(shí)際上往往還需對(duì)h1（t）按以上方法進(jìn)行多次篩選［1］，直到其滿足某個(gè)終止準(zhǔn)則（如SD＜0.2［1］），就認(rèn)為成為了一個(gè)IMF，記為c1（t），如圖1（d）所示。將s（t）減去已經(jīng)篩選出的IMFc1（t）后再用同樣的方法可以依次篩選出s（t）中的其余IMF及趨勢項(xiàng)。最終s（t）可以表示為：

圖1 EMD過程示意圖Fig.1 The illustration of EM

從以上介紹可以看出，EMD是利用極值點(diǎn)擬合信號(hào)包絡(luò)線或均值曲線對(duì)信號(hào)進(jìn)行分解的。實(shí)際中獲得的信號(hào)往往是被截?cái)嗪蟮男盘?hào)，這類信號(hào)內(nèi)部的極值點(diǎn)總是能獲得，但邊界外的極值點(diǎn)或邊界包絡(luò)與均值只有經(jīng)過預(yù)測才能得出，這就是HHT中的EMD邊界處理問題，它是HHT的一個(gè)重要問題［1］。如引言所述，目前已經(jīng)有若干EMD邊界處理方法，但這些方法還有待改善。筆者通過分析認(rèn)為，在已有的EMD邊界處理方法中波形匹配預(yù)測法是一種相對(duì)更實(shí)用的方法，因?yàn)樗门c邊界波最相匹配的內(nèi)部波來預(yù)測信號(hào)的邊界均值，與許多實(shí)際信號(hào)具有波形重復(fù)性的規(guī)律相一致。但目前的波形匹配預(yù)測法是采用邊界復(fù)合波與內(nèi)部復(fù)合波的匹配性來預(yù)測邊界均值的，使其不適用于有調(diào)幅或疊加波的復(fù)合信號(hào)［10］。

2 基于本征波匹配的EMD邊界處理方法

針對(duì)波形匹配預(yù)測法的不足，本文第一作者在文獻(xiàn)［10］中曾提出了一種改進(jìn)的波形匹配預(yù)測法——本征波匹配預(yù)測法。它是基于對(duì)IMF本身結(jié)構(gòu)的深入認(rèn)識(shí)和所建立的本征波概念而提出的。

一個(gè)IMFc（t）可以看成是由cosφ（t）經(jīng)包絡(luò)線a（t）調(diào)制而成的，如圖2（a）、圖 2（b）所示，cosφ（t）中包含了c（t）的全部頻率信息，因此稱之為本征信號(hào)［10］。而cosφ（t）中的頻率往往是時(shí)變的，不同時(shí)刻可以有不同的頻率，每處的頻率有一個(gè)對(duì)應(yīng)尺度的子波，如圖2（a）中所示的從一個(gè)極大值點(diǎn)到下一個(gè)相鄰的極大值點(diǎn)形成的信號(hào)段ABC（后面本文理論的應(yīng)用例子表明，對(duì)這種“子波”無需清晰定義，只需大致定義即可），且同一cosφ（t）中不同時(shí)刻的子波可以相互疊加。因此可以認(rèn)為本征信號(hào)cosφ（t）又是由一些不同尺度的子波疊加構(gòu)成，這些子波在文獻(xiàn)［10］中稱為本征波。這樣，在HHT中任何一個(gè)信號(hào)f（t）都可以看成是由一些本征波相互疊加形成本征信號(hào)，本征信號(hào)經(jīng)過調(diào)幅后形成IMF，如圖2（b）所示，IMF再相互疊加就形成復(fù)合信號(hào)f（t），而本征波在復(fù)合信號(hào)中必有一個(gè)一一對(duì)應(yīng)的復(fù)合波，如圖2（c）所示A″B″C″。

圖2 HHT中的復(fù)合信號(hào)形成過程Fig.2 The forming process of composite signal in HHT

下面以起始點(diǎn)的EMD邊界處理為例對(duì)本征波匹配預(yù)測法的算法作一介紹。

設(shè)信號(hào)起始點(diǎn)及連續(xù)相鄰的兩個(gè)極值點(diǎn)分別為s（0）、e（t1）、e（t2），則這三點(diǎn)構(gòu)成一個(gè)邊界復(fù)合波s（0）e（t1）e（t2），其隱含一個(gè)邊界本征波，讓s（0）e（t1）·e（t2）在極值點(diǎn)序列中以步長2逐漸向右移動(dòng)，所移之處都有一個(gè)由頂點(diǎn)e（ti－1），e（ti），e（ti+1）構(gòu)成的內(nèi)部復(fù)合波e（ti－1）e（ti）e（ti+1）［隱含一個(gè)內(nèi)部本征波］與之對(duì)應(yīng)，如圖3（a）所示。本征波匹配預(yù)測法就是為了尋找一個(gè)其本征波與邊界波的本征波匹配誤差最小的內(nèi)部波來預(yù)測邊界均值。為此，首先將波形s（0）e（t1）e（t2）減去直線y（t）=e（t1）（相當(dāng)于去掉疊加波），再除以e（t1）－e（t2）（相當(dāng)于消除調(diào)幅波），最終得波形s0e1e2，它可以被認(rèn)為就是波形s（0）e（t1）e（t2）所隱含的近似邊界本征波，如圖3（b）所示；同理，將波形e（ti－1）e（ti）e（ti+1）減去直線y（t）=e（ti）并除以e（ti）－e（ti+1）得波形e（ti－1）e（ti）e（ti+1）所隱含的近似內(nèi)部本征波ei－1eiei+1。這樣邊界波s（0）e（t1）e（t2）與內(nèi)部波e（ti－1）e（ti）e（ti+1）之間的匹配誤差就可以采用如下本征波匹配誤差估計(jì)公式計(jì)算：

圖3 本征波匹配預(yù)測法示意圖Fig.3 The illustration of character wave matching prediction

當(dāng)計(jì)算完所有誤差E（i）后，確定出最小誤差E（k），此處的內(nèi)部波即被認(rèn)為與起始邊界波最匹配，因此可以利用此內(nèi)部波的均值預(yù)測起始點(diǎn)處的均值，預(yù)測公式如下：

這里不僅給出了估計(jì)起始點(diǎn)處均值的公式，還給出了估計(jì)第一個(gè)極值點(diǎn)對(duì)應(yīng)均值的公式，以提高邊界均值曲線擬合的穩(wěn)定。用同樣的方法可以估計(jì)出信號(hào)右端點(diǎn)處的均值。

根據(jù)以上介紹可以看出，本征波匹配預(yù)測法是利用信號(hào)邊界波和內(nèi)部波之間的本征波匹配性來預(yù)測信號(hào)邊界處的包絡(luò)值和均值，既繼承了波形匹配預(yù)測法的優(yōu)點(diǎn)，同時(shí)又克服了波形匹配預(yù)測法在計(jì)算匹配誤差時(shí)忽略了信號(hào)調(diào)制與疊加對(duì)波形匹配誤差的影響問題，使得本征波匹配預(yù)測法的適用范圍得以大幅擴(kuò)大。

3 本征波匹配預(yù)測法的算例驗(yàn)證

為了驗(yàn)證本征波匹配預(yù)測法的可行性和進(jìn)步性，這里分別用本征波匹配預(yù)測法和兩種具有代表性的EMD邊界處理方法——包絡(luò)極值延拓法與波形匹配法對(duì)圖4（a）、圖5（a）所示的兩個(gè)含調(diào)幅與疊加波的仿真信號(hào)和圖6（a）所示的一實(shí)測空調(diào)振動(dòng)信號(hào)進(jìn)行EMD，所得結(jié)果依次如圖4（b）、圖4（c）、圖4（d）、圖5（b）、圖5（c）、圖5（d）和圖6（b）、圖6（c）、圖 6（d）所示，其中篩選終止準(zhǔn)則為SD＜0.2。由圖4（d）可以看出，以本征波匹配預(yù)測法進(jìn)行EMD邊界處理時(shí)，所得趨勢項(xiàng)與原始信號(hào)的總體趨勢非常一致，且所得IMF1，IMF2的左右兩邊界波形也比較合理，而由圖4（b）、圖4（c）可以看出，以包絡(luò)極值延拓法或波形匹配法進(jìn)行EMD邊界處理時(shí)，所得IMF與趨勢項(xiàng)在右端均容易出現(xiàn)明顯變異。觀察圖5也可以得到同樣的結(jié)論。再觀察圖6可以看出，以包絡(luò)極值延拓法進(jìn)行EMD邊界處理時(shí)，所得IMF2的左端、趨勢項(xiàng)的左右端均出現(xiàn)了明顯變異；以波形匹配法進(jìn)行EMD邊界處理時(shí)所得IMF2左端也變異明顯，雖然所得趨勢項(xiàng)的左右端變異不明顯，但同用本征波匹配預(yù)測法所得圖6（d）所示的趨勢項(xiàng)相比，仍不如后者更接近肉眼所觀察到的原信號(hào)兩端的變化趨勢。另外與圖6（b）、圖6（c）相比，圖6（d）中IMF2的左端變異也有明顯改善。因此，無論對(duì)于仿真信號(hào)還是實(shí)測信號(hào)，用本征波匹配預(yù)測法進(jìn)行EMD邊界處理均能取得比其他EMD邊界處理方法更滿意的效果。

圖4 三種EMD邊界處理方法在對(duì)一仿真信號(hào)進(jìn)行EMD時(shí)的結(jié)果比較Fig.4 The EMD result comparison of three sorts of EMD end process methods to a simulated signal

圖5 三種EMD邊界處理方法在對(duì)另一仿真信號(hào)進(jìn)行EMD時(shí)的結(jié)果比較Fig.5 The EMD result comparison of three sorts of EMD end process methods to another simulated signal

圖6 三種EMD邊界處理方法在對(duì)實(shí)測信號(hào)進(jìn)行EMD時(shí)的結(jié)果比較Fig.6 The EMD result Comparison of three sorts of EMD end process methods to a test signal

4 結(jié)論

本文著重對(duì)文獻(xiàn)［10］中提出的一種基于本征波匹配的EMD邊界處理方法——本證波匹配預(yù)測法進(jìn)行了介紹，特別是介紹了將該方法和其他兩種有代表性的EMD邊界處理方法——包絡(luò)延拓法和波形匹配預(yù)測法應(yīng)用于仿真信號(hào)和和實(shí)測信號(hào)的算例對(duì)比，從而驗(yàn)證了本證波匹配預(yù)測法的先進(jìn)性。本文工作使文獻(xiàn)［10］中的相關(guān)研究得到較大程度地深化和完善。

［1］Huang N E，Zheng S，Long S R，et al.The empirical mode decomposition and the Hilbert spectrum for nonlinear and non-stationary time series analysis［J］.P roc.R.Soc.Lond.A ，1998，454：903－995.

［2］Pen Z K，Tse P W，Chu F L.An improved Hilbert Huang transform and its application in vibration signal analysis［J］.Journal ofSoundandVibration， 2005， 286 （122）：187－205.

［3］蓋 強(qiáng)，馬孝江，張海勇，等.一種消除局域波法邊界效應(yīng)的新方法［J］.大連理工大學(xué)學(xué)報(bào)，2002，42（1）：115－117.

［4］王 婷，楊莘元，李冰冰.一種改善EMD端點(diǎn)效應(yīng)的新方法［J］.哈爾濱理工大學(xué)學(xué)報(bào)，2009，14（5）：23－26.

［5］鄧擁軍，王 偉，錢成春，等.EMD方法及Hilbert變換中邊界問題的處理［J］.科學(xué)通報(bào)，2001，46（3）：257－263.

［6］黃大吉，趙進(jìn)平，蘇紀(jì)蘭.希爾伯特－黃變換的端點(diǎn)延拓［J］.海洋學(xué)報(bào)，2003，25（1）：1－11.

［7］劉慧婷，張 旻，程家興.基于多項(xiàng)式擬合算法的EMD端點(diǎn)問題的處理［J］.計(jì)算機(jī)工程與應(yīng)用，2004，40（16）：84－86.

［8］徐世艷.基于延拓技術(shù)的經(jīng)驗(yàn)?zāi)B(tài)分解改進(jìn)算法［J］.吉林大學(xué)學(xué)報(bào)，2009，27（4）：359 －365.

［9］沈 路，周曉軍，張志剛，等.Hilbert-Huang變換中的一種端點(diǎn)延拓方法［J］.振動(dòng)與沖擊，2009，28（8）：168－175.

［10］鐘佑明.希爾伯特－黃變換局瞬信號(hào)分析理論的研究［D］.重慶：重慶大學(xué)，2002.