基于灰度預(yù)測的變長度經(jīng)驗?zāi)B(tài)分解方法研究

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

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

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

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

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

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

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

1.1 EMD基本原理

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

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

將原信號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ù)進行上述的篩分工作,直到所得序列滿足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 端點效應(yīng)

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

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

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

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

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

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

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

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

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

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

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

該式的解為

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

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

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

還原值為

此即預(yù)測方程。

2.2 殘差模型

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

及相對殘差序列

進而得到平均相對殘差

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

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

式中,δε為校正權(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)的最大值和最小值。對進行灰度預(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,利用修正模型求得進而求得最后可得到原始數(shù)據(jù)x(0)(t)的全部預(yù)測值。

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

4 模型仿真實驗

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

4.1 正弦疊加信號的分解

實驗中用于分析的信號為

式中,信號的采樣頻率為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)驗?zāi)B(tài)分解方法、基于灰度模型的方法以及本研究提出的改進方法進行信號分解,可以得到如圖4～6所示的結(jié)果。利用式(22)對分解得到的IMF分量與理論結(jié)果進行誤差計算：

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

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

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

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

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

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

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

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

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

5 結(jié) 論

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

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

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