局部特征尺度分解改進(jìn)算法

盛 沛，崔偉成，楊永彬，許愛強(qiáng)

（1.海軍航空大學(xué)，山東煙臺264001；2.91291部隊(duì)，海南三亞572000）

傳統(tǒng)的以傅里葉變換為基礎(chǔ)的信號分析與處理方法只適用于平穩(wěn)信號，無法表述信號的時(shí)頻局部特性。為了準(zhǔn)確將信號進(jìn)行分解，從而得到分量的局部特征，眾多學(xué)者開展了時(shí)頻分析領(lǐng)域的研究。近年來，以經(jīng)驗(yàn)?zāi)B(tài)分解（Empirical Mode Decomposition，EMD）[1-3]為代表的分析方法逐漸成為熱點(diǎn)，EMD與希爾伯特變換（Hilbert Transform，HT）技術(shù)相結(jié)合的希爾伯特-黃變換（HHT），是一種有效的自適應(yīng)時(shí)頻分析方法，被廣泛地應(yīng)用于旋轉(zhuǎn)機(jī)械部件故障、心/腦電圖等的特征提取與診斷領(lǐng)域[4]。但EMD 在使用過程中仍然存在一些問題，比如，過/欠包絡(luò)、端點(diǎn)效應(yīng)、模態(tài)混疊等。因此，大量學(xué)者開始致力于對EMD 進(jìn)行改進(jìn)和替代的算法研究。

其中，比較有代表性的有Smith 于2005 年提出了局部均值分解（Local Mean Decomposition，LMD）方法[5]，以及程軍圣團(tuán)隊(duì)提出的局部特征尺度分解方法（Local Characteristic-Scale Decomposition，LCD）[6-7]。

2013年，程軍圣等通過仿真信號對比了LCD及EMD，證明前者在端點(diǎn)效應(yīng)、計(jì)算時(shí)間等方面都優(yōu)于后者[8]。

但是，LCD 仍然存在均值曲線定義不盡合理、分解過程模態(tài)混疊等問題。從算法過程來看，LCD方法的均值點(diǎn)構(gòu)造本質(zhì)上是反復(fù)利用3個(gè)相鄰極值點(diǎn)進(jìn)行線性插值、反復(fù)迭代，包含了3/4波長的信息。從這一角度來看，該方法的分解質(zhì)量必然優(yōu)于包含半波長信息的LMD 方法。但是，標(biāo)準(zhǔn)LCD 方法采用線性的方法估計(jì)包絡(luò)曲線時(shí)，沒有考慮曲線段的“上凸”或“下凹”特性，這必然造成了信息的“丟失”。

基于以上思路，本文嘗試對LCD算法中均值點(diǎn)求取部分做出適當(dāng)改進(jìn)，目的在于完善LCD的基本理論和分解過程。

1 LCD算法中均值曲線的構(gòu)造

LCD 方法的原理[9]是基于單分量信號“局部關(guān)于零均值對稱”的結(jié)論，其基礎(chǔ)在于構(gòu)造均值曲線，并將均值曲線不斷從原始信號中分離，直至信號為內(nèi)稟尺度分量，其構(gòu)造如圖1所示。

圖1 標(biāo)準(zhǔn)LCD算法中局部均值曲線的構(gòu)造Fig.1 Construction of local mean curve in standard LCD

局部均值曲線的構(gòu)造方法如下。

1）將任意2 個(gè)相同類型的極值點(diǎn)( τk,Xk)、( τk+2,Xk+2)連接成線段，求出τk+1時(shí)刻對應(yīng)的Ak+1。

2）由( τk+1,Ak+1)與( τk+1,Xk+1) 得到局部均值點(diǎn)

式中，a ∈( 0,1) 為一常量，典型地，a=0.5。

3）依據(jù)時(shí)刻τk將x( t )分割成若干區(qū)間，在每個(gè)區(qū)間（任意2個(gè)相鄰極值點(diǎn)之間）對x( t )進(jìn)行線性變換，式（2）中，Hk( t )表示信號的第k 個(gè)區(qū)間的局部均值曲線。

4）將Hk( t )依次連接得到均值曲線。

對任意實(shí)信號x( t )，將均值曲線不斷從信號中扣除，就可得到ISC分量，這是LCD的思路。

2 LCD方法的改進(jìn)

如前所述，估計(jì)一個(gè)極值點(diǎn)處的均值點(diǎn)應(yīng)考慮相鄰極值點(diǎn)的信息；合理地添加臨近的極值點(diǎn)作為約束會(huì)優(yōu)化LCD 均值點(diǎn)的估計(jì)效果。但添加的約束極值點(diǎn)不能過多，主要原因有：①LCD方法定義的基線為“局部”均值曲線，計(jì)算一個(gè)均值點(diǎn)考慮過多的極值點(diǎn)有悖于LCD 的初衷；②增加過多的約束極值點(diǎn)會(huì)增加估計(jì)的限制條件，從而降低LCD 對低頻、短數(shù)據(jù)集的分解能力。因此，在計(jì)算某一個(gè)極值點(diǎn)對應(yīng)的均值點(diǎn)時(shí)，可將該極值點(diǎn)和相鄰3 個(gè)不同類型的極值點(diǎn)作為計(jì)算條件，提高均值點(diǎn)的估計(jì)精度[10-11]。

以估計(jì)極小值點(diǎn)( τmin1,Xmin1) 處的均值點(diǎn)為例，如圖2所示，可估計(jì)上包絡(luò)曲線，進(jìn)而求出上包絡(luò)曲線在τmin1時(shí)刻的值X′max1，對Xmin1、X′max1求均值可得到均值點(diǎn)的估計(jì)值。為了使上包絡(luò)曲線足夠光滑（在極值點(diǎn)( τmax2,Xmax2)處二階微分連續(xù)），選用3 次樣條插值（Cubic Spline Interpolation，CSI）方法插值上包絡(luò)曲線。顯然，采用類似的方法也可求出下包絡(luò)曲線，進(jìn)而得到極大值( τmax2,Xmax2)處的局部均值。

圖2 改進(jìn)的均值曲線構(gòu)造方法Fig.2 Construction of local mean curve in improved LCD

上述方法中，插值的方向以及向前還是向后延拓并不是影響分解質(zhì)量的關(guān)鍵因素，但的確會(huì)對分解結(jié)果的端點(diǎn)產(chǎn)生較大影響。在實(shí)際應(yīng)用中，可根據(jù)需要靈活選取[12]。更進(jìn)一步的，2 種方法可以綜合起來得到如下雙向形式：

式（3）中：Lb及Lf分別為后向和前向插值方法獲得的均值點(diǎn)。

實(shí)際上，EMD、LCD 及其改進(jìn)算法的基本思想是一致的，都是根據(jù)信號的極值點(diǎn)分布情況，利用循環(huán)迭代的方式，得到所需分量，其過程如圖3所示。

圖3 LCD、EMD和LMD方法的統(tǒng)一框架Fig.3 General framework of LCD，EMD and LMD

參照標(biāo)準(zhǔn)LCD算法，上述改進(jìn)算法的具體過程如下。

1）將信號x( t )的所有極值點(diǎn)Xk及對應(yīng)的時(shí)刻τk（k=1,2,…,M）分別構(gòu)造時(shí)間序列，為了減少分解的端點(diǎn)效應(yīng)，對2 個(gè)時(shí)間序列進(jìn)行延拓。具體的延拓方式可在兩端各增加一個(gè)極值點(diǎn)[4]：

延拓后的時(shí)間序列分別為Xk、τk，其中：k=0,1,…,M+1。

2）取極值點(diǎn)X( τk-1) 、X( τk)、X( τk+1)及X( τk+3)，其中，k=1,2,…,M-2。

3）用分段3 次樣條方法計(jì)算局部包絡(luò)曲線u( t ),t ∈[τk-1,τk+3] ，u( t )的約束條件：u( τk-1) =X( τk-1) ；u( τk+1) =X( τk+1) ；u( τk+3)=X( τk+3)；u′( τk-1) =s′( τk-1) ；u′( τk+1) =s′( τk+1)；u′( τk+3)=s′( τk+3)；u″( τk-1) =s″( τk-1)；u″( τk+1)=s″( τk+1) ；u″( τk+3)=s″( τk+3)。

式（6）中，s( t )為原始信號。

局部上包絡(luò)滿足u( t )≥s( t )，局部下包絡(luò)滿足u( t )≤s( t )。

4）按式（7）分別計(jì)算τk時(shí)刻的局部均值點(diǎn)Lb( τk)及Lf( τk)，（k=1,2,…,M-2）：

5）按式（2）給出的分段線性方法計(jì)算局部均值曲線Hk( t )。

6）由Hk( t )依次連接成均值曲線H1( t )，并將均值曲線H1( t )從原始信號中分離出來，即

若h1( t )是一個(gè)ISC 分量，輸出的ISC1( t )=h1( t )。否則，將h1( t )作為原始信號，將步驟1）～6）重復(fù)循環(huán)k-1次，得到內(nèi)稟尺度分量ISC1( t )=h1k( t )。

7）將ISC1( t )從信號x( t )中分離出來，可得到1 個(gè)新的剩余信號r1( t )，即

8）將r1( t )視為原始數(shù)據(jù)，將步驟1）～6）重復(fù)循環(huán)n-1 次，直至rn( t )單調(diào)（無極值點(diǎn)）。就將x( t )分解為n 個(gè)ISC 分量和剩余信號之和[8]，即

3 改進(jìn)算法性能驗(yàn)證

3.1 改進(jìn)算法與原算法分解效果比較

為了驗(yàn)證均值點(diǎn)改進(jìn)算法的效果，考察如下示例信號：

式（11）中，x( t )為仿真信號，由信號x1( t )、x2( t )合成，時(shí)域波形如圖4所示。

圖4 仿真信號Fig.4 Simulation signal

對該仿真信號進(jìn)行標(biāo)準(zhǔn)LCD 分解以及改進(jìn)的雙向插值LCD分解，結(jié)果如圖5、6所示。

圖5 標(biāo)準(zhǔn)LCD方法Fig.5 Standard LCD

圖6 改進(jìn)LCD方法Fig.6 Improved LCD

從分解結(jié)果來看：

1）標(biāo)準(zhǔn)LCD 算法與改進(jìn)算法均將原始信號分解為3個(gè)分量ISC分量及剩余分量r，由于算法原理導(dǎo)致頻率較高的x3( t )被先分解出，分解結(jié)果與原始信號總體趨勢一致；

2）改進(jìn)算法中，2 個(gè)分量在端點(diǎn)處失真較小，且ISC2分解效果較原算法略高，是因?yàn)楦倪M(jìn)方法減小了局部均值點(diǎn)的估計(jì)誤差，得到了更準(zhǔn)確的均值曲線。

3.2 噪聲對分解結(jié)果的影響

在實(shí)際工程應(yīng)用中，由于受環(huán)境、采集設(shè)備等因素的影響，算法所面對的信號一定是混有噪聲的。因此，分析噪聲對算法的影響顯得尤為必要。為了不失一般性，將原始信號頻率一并納入考察范圍，考察指標(biāo)為各算法在加噪情況下所分解得到的各分量與原始信號中相應(yīng)分量的互相關(guān)系數(shù)?？疾烊缦滦盘枺?/p>

式（12）中：noise(t,SNR)表示信噪比為SNR 的高斯白噪聲[13]，SNR按每個(gè)采樣點(diǎn)計(jì)算，單位dB；f1、f2分別為可變頻率。

具體實(shí)現(xiàn)步驟為：

1）設(shè)定初始頻率f1、f2，初始信噪比SNR，利用式（12）得到初始化的含噪信號；

2）用指定算法對含噪信號進(jìn)行分解；

3）將信號x1( t )、x2( t )分別于所有分量求取互相關(guān)系數(shù)，并分別取最大值作為最終的互相關(guān)系數(shù)[14]；

4）分別增大頻率f1、f2，信噪比SNR重復(fù)步驟2）～3）。

按照上述步驟，得到3種算法[15-17]分解效果隨SNR及頻率變化圖如圖7～9所示，這里只給出x1分量相關(guān)系數(shù)評價(jià)指標(biāo)，x2分量與之類似?？梢钥闯觯?/p>

a）當(dāng)信噪比大于約20 dB 以后，各算法分解效果十分理想，先后達(dá)到1。其中，改進(jìn)的LCD算法相關(guān)性指標(biāo)更早達(dá)到1。

b）原始信號頻率對各算法分解效果有一定影響，使指標(biāo)達(dá)到1的速度變緩。

c）步驟1）中每次都會(huì)對純凈信號重新加入高斯白噪聲，因而計(jì)算結(jié)果會(huì)有起伏波動(dòng)且每次仿真結(jié)果都略有初入。這也是為了更好地檢查各算法的魯棒性。顯然，EMD魯棒性最差，改進(jìn)算法較優(yōu)。

圖7 EMD方法分解效果x1 相關(guān)性指標(biāo)Fig.7 Decomposition effect for x1 correlation index in EMD

圖8 LCD方法分解效果x1 相關(guān)性指標(biāo)Fig.8 Decomposition effect for x1 correlation index in LCD

圖9 改進(jìn)的LCD方法分解效果x1 相關(guān)性指標(biāo)Fig.9 Decomposition effect for x1 correlation index in improved LCD

需要說明的是，上述仿真過程中兩分量信號是一并遞增的。若將兩信號頻率間拉開一定距離，分解效果將有較大變化，在這里不多做討論。

4 結(jié)束語

本文所提出的改進(jìn)算法在給出具體實(shí)現(xiàn)步驟的同時(shí)，通過仿真實(shí)例驗(yàn)證了其分解質(zhì)量的確有一定提升。盡管如此，該方法仍有其不足之處，如改進(jìn)算法也存在模態(tài)混疊問題等。下一步，將結(jié)合抗混疊等方法進(jìn)一步完善該算法。另外，在實(shí)際應(yīng)用中，還可以考慮先對信號進(jìn)行去噪處理。比如，結(jié)合小波閾值去噪和卡爾曼濾波。