張曉波劉保華于凱本劉 苗楊志國于盛齊宗 樂

(1.自然資源部 國家深?；毓芾碇行?山東 青島266237；2.青島海洋科學(xué)與技術(shù)試點(diǎn)國家實(shí)驗(yàn)室 海洋地質(zhì)過程與環(huán)境功能實(shí)驗(yàn)室,山東 青島266061；3.中國海洋大學(xué) 環(huán)境科學(xué)與工程學(xué)院,山東 青島266100；4.青島海洋科學(xué)與技術(shù)試點(diǎn)國家實(shí)驗(yàn)室 公共關(guān)系部,山東 青島266061)

對于一般的非線性非平穩(wěn)信號,信號特征多表現(xiàn)為時變和頻帶混疊,導(dǎo)致有效信號往往存在于某一個或幾個頻帶中,因此從原始信號中分離和提取有效信號是信號處理的重要研究內(nèi)容之一。經(jīng)驗(yàn)?zāi)B(tài)分解(Empirical Mode Decomposition,EMD)是由美籍華人Huang等[1]于1998年提出的一種信號譜分解方法。該方法的核心思想是任何復(fù)雜的信號都可以被分解成有限個固有模態(tài)函數(shù)(Intrinsic Mode Function,IMF)和1個剩余分量之和。因?yàn)樵摲椒ㄊ窃诨谛盘柋旧硖卣鞯幕A(chǔ)實(shí)現(xiàn)的分解,它與基于傅里葉變換的方法相比,具有更強(qiáng)的自適應(yīng)性和局部時頻分析能力,因此在分析非線性和非平穩(wěn)信號中具有很高的實(shí)用價值。

雖然經(jīng)驗(yàn)?zāi)B(tài)分解被認(rèn)為適合分解非線性和非平穩(wěn)信號,但由于其提出時間較短,缺乏理論支撐,需要解決的問題很多,例如在應(yīng)用經(jīng)驗(yàn)?zāi)B(tài)分解時信號邊界很容易產(chǎn)生信號分解失真問題,即端點(diǎn)效應(yīng)。其產(chǎn)生原因是在基于EMD方法進(jìn)行信號分解時,由于端點(diǎn)處為非極值點(diǎn),在計算信號包絡(luò)時會出現(xiàn)端點(diǎn)處不能完全包絡(luò)以及由三次樣條插值引起的過沖和欠沖等問題,這會使得分解后的IMF函數(shù)誤差過大,從而對信號分解結(jié)果的準(zhǔn)確度產(chǎn)生嚴(yán)重的影響。

目前國內(nèi)外學(xué)者已提出了多種壓制端點(diǎn)效應(yīng)的方法,如特征波法[1],波形匹配法[2],鏡像圓周對稱延拓法[3-4],極值點(diǎn)鏡像延拓法[5],基于時間序列預(yù)測模型(Auto Regression Moving Average,ARMA)的延拓方法[6],基于神經(jīng)網(wǎng)絡(luò)的時間序列預(yù)測方法[7-8]以及基于灰度預(yù)測模型的延拓方法[9]等。但這些方法在應(yīng)用中均存在一定程度的不足。例如,特征波法的問題在于非平穩(wěn)非線性信號的兩端具有的波形特征不一定相同,所以選擇合適的特征波比較困難；波形匹配法在實(shí)際應(yīng)用中計算效率相對低下；鏡像圓周對稱延拓法則將端點(diǎn)直接做為極值點(diǎn)進(jìn)行處理,反復(fù)迭代計算過程中會引起端點(diǎn)值的失真；極值點(diǎn)鏡像延拓法通常把距離端點(diǎn)最近的極值點(diǎn)作為鏡面點(diǎn)對信號進(jìn)行延拓,當(dāng)端點(diǎn)不是極值點(diǎn)時還需對信號進(jìn)行截斷處理,如果數(shù)據(jù)較短,則端點(diǎn)效應(yīng)的壓制會失效；基于時間序列預(yù)測模型的延拓方法和基于神經(jīng)網(wǎng)絡(luò)的延拓方法兩種方法預(yù)測效果較好,但存在計算效率低下和對大數(shù)據(jù)量和信號周期化的依賴等問題；基于灰度預(yù)測模型的延拓方法根據(jù)已知的極值點(diǎn)位置和數(shù)據(jù)利用灰度預(yù)測模型預(yù)測未知的極值點(diǎn)位置和數(shù)據(jù),與前述其他方法相比,該方法在計算效率和抑制效果等方面均有顯著改善,但由于灰度預(yù)測方法本身的限制,其效果仍有待提升。

本文圍繞經(jīng)驗(yàn)?zāi)B(tài)分解中的端點(diǎn)效應(yīng)問題,對常規(guī)基于灰度預(yù)測模型的端點(diǎn)延拓方法[9]加以改進(jìn),實(shí)現(xiàn)了基于灰度模型的變長度經(jīng)驗(yàn)?zāi)B(tài)分解方法,并利用仿真數(shù)據(jù)進(jìn)行結(jié)果對比。

1 EMD基本原理及端點(diǎn)效應(yīng)

1.1 EMD基本原理

EMD算法基于任何信號都是由有限個固有模態(tài)函數(shù)組成,每一個固有模態(tài)函數(shù)都可以通過以下方法得到：

假設(shè)原始信號為x(t),找出信號的全部極值點(diǎn),將所有局部極大值點(diǎn)和所有局部極小值點(diǎn)分別進(jìn)行擬合得到envmax(t)和envmin(t),這樣使所有的信號數(shù)據(jù)包含在2條包絡(luò)線之間。2條包絡(luò)線的平均值記為m1(t),其表達(dá)式為

將原信號x(t)減去m1(t)得到一個新信號h1(t),即h1(t)=x(t)-m1(t),此過程稱為篩分。若h1(t)滿足IMF條件[5],則h1(t)為第1個固有模態(tài)函數(shù)；若h1(t)不滿足IMF條件,則將所得的h1(t)序列作為待處理數(shù)據(jù),重復(fù)進(jìn)行上述的篩分工作,直到所得序列滿足IMF條件為止。這樣便得到信號x(t)的第1個固有模態(tài)函數(shù),記為c1。將c1從x(t)中分離出來,可以得到去除掉高頻成分后的剩余部分r1=x(t)-c1,然后把r1作為原始數(shù)據(jù)重復(fù)以上過程,可以得到x(t)的第2個固有模態(tài)函數(shù)c2。以此類推,如果最終得到剩余分量(rn)比預(yù)定的誤差要小,或當(dāng)rn變成1個單調(diào)函數(shù)不能再繼續(xù)篩分時,則停止上述分解過程。最終原始信號可以表示為

式中,ci為IMF分量。

1.2 端點(diǎn)效應(yīng)

在EMD分解過程中,對極小值和極大值包絡(luò)的求取通常采用插值的方法。當(dāng)信號端點(diǎn)處不為極值點(diǎn)時,包絡(luò)計算將產(chǎn)生明顯誤差,由此分解得到的第i個IMF分量比第i-1個IMF分量在端點(diǎn)處會由于誤差累計產(chǎn)生更大的誤差。與此同時,誤差還會由數(shù)據(jù)兩端向內(nèi)傳播,從而導(dǎo)致整個分解過程失去意義。

這里對信號端點(diǎn)作舉例說明。設(shè)原始信號為式中,t∈[0,0.3],采樣頻率為10 000 Hz,對該信號進(jìn)行經(jīng)驗(yàn)?zāi)B(tài)分解時需要求取其極大值包絡(luò)、極小值包絡(luò)以及包絡(luò)的均值(圖1)。

由圖1可見,如果端點(diǎn)處不是極值點(diǎn),那么求取的包絡(luò)在端點(diǎn)處就會出現(xiàn)“欠沖”或“過沖”現(xiàn)象,從而導(dǎo)致計算過程中在端點(diǎn)處出現(xiàn)誤差。

圖1 經(jīng)驗(yàn)?zāi)B(tài)分解的端點(diǎn)效應(yīng)Fig.1 Edge effect of empirical mode decomposition

2 基于灰度模型的經(jīng)驗(yàn)?zāi)B(tài)分解方法

根據(jù)前文所述,為了壓制EMD的端點(diǎn)效應(yīng),需對端點(diǎn)處的極值點(diǎn)進(jìn)行有效預(yù)測[6]。因此這里利用一階單變量灰度預(yù)測模型算法對端點(diǎn)處極值點(diǎn)實(shí)現(xiàn)預(yù)測,并基于廣泛使用的三次樣條插值方法進(jìn)行包絡(luò)計算。

2.1 一階單變量灰度預(yù)測模型

定義X(0)為預(yù)測模型的數(shù)據(jù)集合

X(1)為X(0)的一次累加生成數(shù)序列

Z(1)為X(1)的緊鄰均值生成序列

假定建模序列同號,則其累加生成數(shù)序列為單調(diào)序列。如果將這個單調(diào)序列用指數(shù)曲線進(jìn)行擬合,則該序列可以寫成的微分方程為

該式的解為

式中,a,b是待定系數(shù)。其中a,b是待定系數(shù)。

式(7)通常被稱為上式的白化方程或影子方程。設(shè)為未知向量則的最小二乘估計滿足

根據(jù)以上論述,式(9)所示的灰色微分方程的時間響應(yīng)序列為

還原值為

此即預(yù)測方程。

2.2 殘差模型

對于上文中得到的預(yù)測模型,其精確程度可通過殘差值進(jìn)行檢驗(yàn)。根據(jù)x(0)(i)與建立絕對殘差序列

及相對殘差序列

進(jìn)而得到平均相對殘差

給定α,當(dāng),且φn＜α成立時,稱所求出的預(yù)測模型為殘差合格模型。

此外,利用殘差值還可對預(yù)測模型值進(jìn)行修正,令殘差值序列為

式中,δε為校正權(quán)重系數(shù),為殘差模型的未知向量利用可得修正后的預(yù)測模型值。

則修正模型寫作

2.3 灰度預(yù)測流程

對于符號相同的數(shù)據(jù),利用灰度預(yù)測方法往往能夠得到較好的預(yù)測結(jié)果,因此,為使原始數(shù)據(jù)滿足這個條件,可對原始數(shù)據(jù)x(0)(t)(t=1,2,3,…,n)做如下變換：

式中,xmax和xmin是x(0)(t)的最大值和最小值。對進(jìn)行灰度預(yù)測得到預(yù)測數(shù)據(jù)然后利用反變換式

得到原始數(shù)據(jù)x(0)(t)的預(yù)測值。

灰度預(yù)測模型的算法流程如下：預(yù)測的次數(shù)為m,預(yù)測的序列為

3)根據(jù)式(10)構(gòu)建矩陣B和Y,并計算,求得[a1,b1]T；

4)對于t=1,2,3,…,n,將[a1,b1]T代入式(11)求得

6)計算原始序列x(0)(t)與的相對殘差序列φ,并計算平均相對殘差。給定α,當(dāng),且φn＜α成立時,稱模型為殘差合格模型；當(dāng)φ＞α或φn＞α?xí)r,計算殘差模型輸入數(shù)據(jù)ε(0)(t),并計算殘差模型和修正模型,利用修正模型求得,然后重復(fù)步驟5)和6),直至得到殘差合格模型。

7)對于t=n+1,…,n+m,利用修正模型求得進(jìn)而求得最后可得到原始數(shù)據(jù)x(0)(t)的全部預(yù)測值。

2.4 基于灰度預(yù)測的端點(diǎn)延拓

基于灰度預(yù)測的端點(diǎn)延拓方法主要用于極值點(diǎn)序列的延拓,利用已知的n個極值點(diǎn)位置和數(shù)據(jù),預(yù)測m個未知的極值點(diǎn)位置和數(shù)據(jù)。以極大值的延拓為例,在EMD每次迭代的過程中,當(dāng)?shù)玫綐O大值點(diǎn)序列之后,分別利用兩端的4個已知極大值構(gòu)建灰度預(yù)測模型,然后向兩端分別延拓2個極大值點(diǎn),從而得到新的極大值點(diǎn)序列,極小值點(diǎn)延拓同理可得。

端點(diǎn)延拓的具體步驟如下：假設(shè)用于預(yù)測的n個極值點(diǎn)的位置序列為lk(k=1,2,…,n),極值序列為vk(k=1,2,…,n)；預(yù)測出的m個極值點(diǎn)位置序列為lk(k=n+1,n+2,…,n+m),極值序列為vk(k=n+1,n+2,…,n+m)。對于i=1,2,3,…,m,第i個預(yù)測的極值點(diǎn)位置li+n可以用{li,li+1,…,li+n-1}進(jìn)行預(yù)測,第i個預(yù)測的極值vi+n可以用{vi,vi+1,…,vi+n-1}進(jìn)行預(yù)測。

3 基于灰度預(yù)測的經(jīng)驗(yàn)?zāi)B(tài)分解的改進(jìn)

3.1 常規(guī)基于灰度預(yù)測的經(jīng)驗(yàn)?zāi)B(tài)分解方法存在的問題

常規(guī)灰度預(yù)測模型在端點(diǎn)延拓方面相對于鏡像延拓,對稱延拓,基于AR模型延拓和基于神經(jīng)網(wǎng)絡(luò)延拓方法有更好的效果,但灰度預(yù)測方法在應(yīng)用中仍存在一些問題：

1)存在矩陣B為0的情況：在灰度預(yù)測的過程中,需要求解矩陣B,然而由于對數(shù)據(jù)進(jìn)行了標(biāo)準(zhǔn)化處理,使得矩陣B中的第1列數(shù)據(jù)存在全為0的情況,所以導(dǎo)致計算出現(xiàn)問題。

2)極值點(diǎn)位置預(yù)測不準(zhǔn)確的情況：極值點(diǎn)位置序列存在預(yù)測過程中仍大于0的情況,沒有對端點(diǎn)處的值做有效預(yù)測。例如,存在對極值點(diǎn)位置的預(yù)測始終大于0的情況,這樣對端點(diǎn)處的包絡(luò)求取仍存在問題。

3)未考慮端點(diǎn)極值情況導(dǎo)致結(jié)果不準(zhǔn)確：對復(fù)雜多變的信號,特別是在端點(diǎn)處有劇烈變化的信號,若不考慮端點(diǎn)極值情況,預(yù)測的結(jié)果將導(dǎo)致端點(diǎn)處欠包絡(luò)。與之相較,極值點(diǎn)鏡像延拓法,極值點(diǎn)對稱延拓法則能有效解決該問題。

4)信號端點(diǎn)處的污染往往無法避免,并會影響下一階固有模態(tài)函數(shù)的精度。

3.2 變長度經(jīng)驗(yàn)?zāi)B(tài)分解方法

針對以上問題,本研究提出基于類對稱延拓和灰度模型結(jié)合的方法進(jìn)行數(shù)據(jù)的經(jīng)驗(yàn)?zāi)B(tài)分解處理,在EMD分解過程中,對信號做端點(diǎn)延拓實(shí)際上都是依靠原有數(shù)據(jù)的特征做出的預(yù)測,對于非線性非穩(wěn)定信號,如果端點(diǎn)處的變化毫無規(guī)律又變化劇烈,通過各種預(yù)測方法得到的結(jié)果可能都是徒勞。由于預(yù)測不可能完全準(zhǔn)確,可能多次迭代之后端點(diǎn)處的值出現(xiàn)發(fā)散,并將向內(nèi)污染整個序列。鑒于此,研究通過在端點(diǎn)處向外對稱延拓一部分原始信號,使其變成新的信號,在每得到一個IMF分量之后,在兩端切掉一定長度的序列,這段序列為預(yù)測不準(zhǔn)確的序列,切除之后,這種預(yù)測不準(zhǔn)確序列不會被代入下一個IMF的求解,從而最大程度的保證中間信號不被污染。

基于灰度預(yù)測的變長度經(jīng)驗(yàn)?zāi)B(tài)分解流程：

1)設(shè)原始數(shù)據(jù)的長度為N,原始數(shù)據(jù)序列為xi,其中i=0,1,2,…,N-1。設(shè)各端延拓的長度為Length,本實(shí)驗(yàn)中令Length=N/10,IMF的最大階數(shù)為M,則每次切掉的長度為step=Length/M。

2)建立新的序列,其中i=0,1,2,…,N+2Length-1,對新序列進(jìn)行賦值；對于k=0,1,2,…,N-1,令x⌒k+Length=xk。

3)根據(jù)端點(diǎn)處數(shù)據(jù)的大小判斷延拓的方式,以左端點(diǎn)為例：

式中,當(dāng)x0的幅值較大時,以該點(diǎn)為對稱點(diǎn)做對稱延拓；當(dāng)其幅值介于原始信號最小值ximin與原始信號最大值ximax之間時,以該點(diǎn)為中心點(diǎn)做類中心對稱延拓,δ∈[0,1]為權(quán)系數(shù)。右端點(diǎn)處延拓方式同理可得。

4)然后對新序列進(jìn)行經(jīng)驗(yàn)?zāi)B(tài)分解,其端點(diǎn)延拓方式選擇基于灰度模型的延拓。

5)當(dāng)?shù)玫絢階IMF分量之后,將中下標(biāo)為i=Length,Length+1,…,Length+N-1的N個點(diǎn)的值,賦值給,其中i=0,1,…,N-1。

6)將步驟5)中得到的k階段IMF分量從步驟4)所述的新序列中減掉,得到剩余分量,其中i=0,1,2,…,N+2Length-1。然后,在剩余分量的兩端分別減去step個點(diǎn)的值后得到,且令Length=Length-step,則序列長度為i=0,1,2,…,N+2Length-1,以為新序列再次執(zhí)行第4),5)和6)步驟,直到滿足EMD分解結(jié)束的條件。

4 模型仿真實(shí)驗(yàn)

模型仿真實(shí)驗(yàn)的內(nèi)容包括基于本研究的方法和灰度預(yù)測的EMD方法分別對正弦疊加信號的分解以及非線性非平穩(wěn)信號分解效果的比較。

4.1 正弦疊加信號的分解

實(shí)驗(yàn)中用于分析的信號為

式中,信號的采樣頻率為500 Hz,合成后的信號如圖2所示。在理想條件下,此信號應(yīng)該被分解成2個正弦IMF分量x11=sin(10πt)和x12=6sin(50πt)以及一個剩余分量x13=0,各分量信號如圖3所示。

圖2 合成信號x 1(t),t∈[0,1]Fig.2 Synthesized signals x 1(t),t∈[0,1]

分別用常規(guī)經(jīng)驗(yàn)?zāi)B(tài)分解方法、基于灰度模型的方法以及本研究提出的改進(jìn)方法進(jìn)行信號分解,可以得到如圖4～6所示的結(jié)果。利用式(22)對分解得到的IMF分量與理論結(jié)果進(jìn)行誤差計算：

式中,IMFi(t)代表在t采樣點(diǎn)處第i階固有模態(tài)函數(shù)的值。

經(jīng)計算得到,基于灰度預(yù)測的經(jīng)驗(yàn)?zāi)B(tài)分解所得結(jié)果的誤差為(±3.337 8),本研究提出的方法的誤差為(±0.025 4)。由此可見,本研究提出的改進(jìn)方法的效果要優(yōu)于改進(jìn)之前的結(jié)果。

圖3 信號子成分x 11,x 12以及x 13Fig.3 Signal subspace x 11,x 12 and x 13

圖4 常規(guī)經(jīng)驗(yàn)?zāi)B(tài)分解的結(jié)果Fig.4 Results of conventional empirical mode decomposition

圖5 基于灰度模型的經(jīng)驗(yàn)?zāi)B(tài)分解的結(jié)果Fig.5 Results of empirical mode decomposition based on grayscale model

圖6 基于本文提出的改進(jìn)方法分解的結(jié)果Fig.6 Results of decomposition based on the improved method proposed in this study

4.2 非線性非平穩(wěn)信號的分解

本實(shí)驗(yàn)中用于分析的信號為

信號的采樣頻率為1 Hz,其波形如圖7所示。在理想條件下,此信號可以被分解成2個正弦IMF分量x21(t)=15sin[(t/80)1.5π]和x22(t)=3sin[(t/100)0.8π]及1個線性剩余分量x23(t)=t/500+5,各分量波形如圖8所示。

這里分別用常規(guī)經(jīng)驗(yàn)?zāi)B(tài)分解方法、基于灰度模型的方法以及本文提出的方法進(jìn)行信號分解,得到如圖9～11所示的結(jié)果。對比3種方法的處理結(jié)果可以看出：在波形吻合度方面,利用本文提出的方法所得到的固有模態(tài)函數(shù)與原信號吻合程度最高,而且對信號兩端的相位信息處理更準(zhǔn)確,基于灰度模型方法所得結(jié)果的吻合程度次之,常規(guī)經(jīng)驗(yàn)?zāi)B(tài)分解方法得到的固有模態(tài)函數(shù)的個數(shù)與原信號不一致,其吻合程度最差；在計算誤差方面,利用式(23)分別計算基于灰度模型的方法以及本文所提出方法的結(jié)果誤差,其中基于灰度預(yù)測的經(jīng)驗(yàn)?zāi)B(tài)分解所得結(jié)果的誤差為(±527.8),本研究所提出方法的誤差為(±311.619 5),因此本文所提出方法的結(jié)果誤差更小。綜上,本研究提出方法的分解效果要優(yōu)于未改進(jìn)之前方法的結(jié)果。

圖7 原始信號x 2(t),t∈[0,3 000]Fig.7 Original signal x 2(t),t∈[0,3 000]

圖8 原始信號的分量x 21(t),x 22(t),x 23(t),t∈[0,3 000]Fig.8 Orignal signal subspace x 21(t),x 22(t),x 23(t),t∈[0,3 000]

圖9 常規(guī)經(jīng)驗(yàn)?zāi)B(tài)分解的結(jié)果Fig.9 Results of conventional empirical mode decomposition

圖10 基于灰度模型經(jīng)驗(yàn)?zāi)B(tài)分解的結(jié)果Fig.10 Results of empirical mode decomposition based on gray model

圖11 基于灰度模型的變長度經(jīng)驗(yàn)?zāi)B(tài)分解的結(jié)果Fig.11 Results of length-varying empirical mode decomposition based on gray model

5 結(jié) 論

針對經(jīng)驗(yàn)?zāi)B(tài)分解方法中的端點(diǎn)效應(yīng)問題展開研究,對常規(guī)基于灰度預(yù)測模型的端點(diǎn)延拓方法存在的問題加以改進(jìn),并提出了基于灰度模型的變長度經(jīng)驗(yàn)?zāi)B(tài)分解方法,然后利用仿真數(shù)據(jù)進(jìn)行對比。模型實(shí)驗(yàn)結(jié)果表明：

1)對于正弦疊加信號,常規(guī)基于灰度預(yù)測的經(jīng)驗(yàn)?zāi)B(tài)分解所得結(jié)果的誤差為(±3.337 8),本文提出的方法的誤差為(±0.025 4),應(yīng)用改進(jìn)方法之后的分解精度得到顯著提升；

2)對于非線性非平穩(wěn)信號,常規(guī)經(jīng)驗(yàn)?zāi)B(tài)分解方法基于灰度預(yù)測的經(jīng)驗(yàn)?zāi)B(tài)分解所得結(jié)果的誤差為(±527.8),用本文提出的方法的誤差為(±311.619 5),應(yīng)用改進(jìn)方法之后的分解精度仍然得到提升,同時,本文提出的方法對信號端點(diǎn)處相位信息的處理更準(zhǔn)確。