基于改進(jìn)EMD方法的系統(tǒng)誤差分離方法研究

姚 剛，趙建軍，姚躍亭，任 喜

(海軍航空工程學(xué)院, 山東 煙臺(tái) 264001)

誤差為測(cè)量值與真值之差。測(cè)量系統(tǒng)中所有測(cè)量結(jié)果均含誤差，其大小直接影響測(cè)量系統(tǒng)精度。測(cè)量誤差分為隨機(jī)誤差與系統(tǒng)誤差。隨機(jī)誤差具有不確定性，其發(fā)生、變化規(guī)律雖無(wú)法確定，但可用統(tǒng)計(jì)方法研究。而系統(tǒng)誤差按確定性規(guī)律變化，其正負(fù)、大小可預(yù)知。為研究、提高測(cè)量系統(tǒng)精度，需對(duì)系統(tǒng)誤差深入分析。系統(tǒng)誤差分離為研究重要內(nèi)容之一。對(duì)復(fù)雜系統(tǒng)，其誤差具有模型難識(shí)別、對(duì)處理結(jié)果影響難估計(jì)的復(fù)雜特性。

系統(tǒng)誤差一般可分為常值誤差、線性漂移誤差、周期性誤差及復(fù)雜規(guī)律變化誤差等。系統(tǒng)誤差分離可視為對(duì)誤差信號(hào)中各頻帶分離過(guò)程。測(cè)量誤差分離方法現(xiàn)使用較多的有最小二乘回歸法、中值濾波法、傅里葉級(jí)數(shù)逼近法、小波及小波包分解法等[1]。各種方法的使用均存在限制，即經(jīng)典最小二乘回歸法的模型階數(shù)較難確定；中值濾波方法精度難以提高；傅里葉級(jí)數(shù)逼近法為對(duì)誤差進(jìn)行諧波分解；小波基函數(shù)選取會(huì)直接影響小波及小波包分解法分離質(zhì)量，增加使用難度。基于上述原因，本文提出用改進(jìn)EMD方法對(duì)系統(tǒng)誤差分析、處理，并進(jìn)行仿真及實(shí)例驗(yàn)證。

1 EMD基本理論

1.1 EMD方法概述

Huang等[2-3]在瞬時(shí)頻率概念基礎(chǔ)上提出的信號(hào)處理方法-基于EMD(Empirical Mode Decomposition)的時(shí)頻分析方法可處理非線性、非平穩(wěn)信號(hào)，分解時(shí)基函數(shù)不確定，只依賴數(shù)據(jù)本身，具有良好的自適應(yīng)特點(diǎn)。EMD方法將數(shù)據(jù)視為由一個(gè)或多個(gè)IMF(Intrinsic Mode Functions，固有模態(tài)函數(shù))組成。由于經(jīng)驗(yàn)?zāi)B(tài)分解方法的優(yōu)點(diǎn)，在生物、海洋、地球科學(xué)、天文學(xué)及工程技術(shù)等領(lǐng)域中廣泛應(yīng)用。實(shí)際的系統(tǒng)誤差往往非平穩(wěn)、非線性，EMD為處理此類復(fù)雜多變的誤差提供了新的方法。

經(jīng)EMD方法分解所得每個(gè)IMF，具有的性質(zhì)為：① 在整個(gè)信號(hào)上，極值點(diǎn)個(gè)數(shù)與過(guò)零點(diǎn)數(shù)目相等或最多相差一個(gè)；② 局部均值為零，即由局部極大值構(gòu)成的上包絡(luò)線與由局部極小值構(gòu)成的下包絡(luò)線的平均值為零。實(shí)際中，上下包絡(luò)均值無(wú)法為零，通常認(rèn)為包絡(luò)均值滿足IMF均值為零的條件為

(1)

1.2 EMD分解步驟

EMD算法假設(shè)對(duì)于何信號(hào)均由若干有限的IMF組成，每個(gè)IMF通過(guò)方法獲得[4-5]。設(shè)原信號(hào)為x(t)，找出整個(gè)信號(hào)x(t)的極值點(diǎn)，并通過(guò)三次樣條擬合出信號(hào)的上下包絡(luò)線，分別為emax(t),emin(t)，計(jì)算上下包絡(luò)線均值為

(2)

將原始信號(hào)減去包絡(luò)線均值，獲得新信號(hào)為

(3)

重復(fù)上述步驟，依次獲得imf2,…,imfn，直至所得余項(xiàng)rn為單調(diào)信號(hào)或小于預(yù)先給定閾值，分解結(jié)束。

原函數(shù)經(jīng)EMD分解得：

(4)

2 基于波形相似度的端點(diǎn)延拓方法

2.1 端點(diǎn)效應(yīng)產(chǎn)生原因

對(duì)誤差信號(hào)進(jìn)行EMD分解時(shí)，需解決端點(diǎn)效應(yīng)問(wèn)題。端點(diǎn)效應(yīng)產(chǎn)生原因即在對(duì)信號(hào)求上下包絡(luò)線時(shí)需對(duì)極值點(diǎn)進(jìn)行擬合。若直接將信號(hào)端點(diǎn)作為信號(hào)的極值點(diǎn)，可使擬合所得包絡(luò)覆蓋整個(gè)信號(hào)。雖滿足經(jīng)驗(yàn)?zāi)B(tài)分解要求條件，但會(huì)導(dǎo)致端點(diǎn)附近包絡(luò)出現(xiàn)較大誤差，且隨篩分過(guò)程的不斷反復(fù)，邊界效應(yīng)產(chǎn)生的誤差會(huì)向信號(hào)內(nèi)部傳播，影響整個(gè)分解過(guò)程。對(duì)較長(zhǎng)數(shù)據(jù)，可采取在分解過(guò)程中通過(guò)拋棄兩個(gè)端點(diǎn)附近數(shù)據(jù)方法削弱端點(diǎn)效應(yīng)。而對(duì)短數(shù)據(jù)而言，只能采取端點(diǎn)延拓方式抑制。

2.2 改進(jìn)的端點(diǎn)延拓方法

目前已有抑制端點(diǎn)效應(yīng)的方法，包括鏡像延拓法、基于神經(jīng)網(wǎng)絡(luò)的延拓法、混沌延拓、基于多項(xiàng)式擬合的延拓方法、基于極值時(shí)間尺度的LS-SVM延拓方法等[6-10]。雖對(duì)端點(diǎn)效應(yīng)抑制有一定效果，但亦存在各自問(wèn)題。鏡像延拓[6]方法處理短數(shù)據(jù)時(shí)效果較差，因其會(huì)截去部分?jǐn)?shù)據(jù)。神經(jīng)網(wǎng)絡(luò)[7]及混沌方法[8]的不足在于速度太慢。而多項(xiàng)式擬合方法[9]適應(yīng)性較差，對(duì)某些數(shù)據(jù)結(jié)果較好，對(duì)另一些數(shù)據(jù)效果會(huì)不理想?；贚S-SVM的端點(diǎn)延拓方法[10]易受內(nèi)積函數(shù)、核函數(shù)、精度參數(shù)及懲罰參數(shù)經(jīng)驗(yàn)選擇影響。為此，本文提出基于波形相似度的端點(diǎn)延拓方法。

端點(diǎn)效應(yīng)抑制原理為尋找信號(hào)之外離端點(diǎn)最近的極值點(diǎn)，與信號(hào)本身的極值點(diǎn)一起進(jìn)行樣條插值，構(gòu)成信號(hào)的上下包絡(luò)線，再截取所需分解數(shù)據(jù)。本質(zhì)上，延拓即為通過(guò)各種數(shù)據(jù)處理方法，對(duì)端點(diǎn)處信號(hào)變化趨勢(shì)進(jìn)行預(yù)測(cè)?；诓ㄐ蜗嗨贫鹊亩它c(diǎn)延拓方法基本思想即能在信號(hào)內(nèi)部找到一段波形與端點(diǎn)處的波形相似，據(jù)相似程度建立數(shù)據(jù)關(guān)系對(duì)原始信號(hào)進(jìn)行延拓?？勺畲笙薅染S護(hù)端點(diǎn)處波形的變化趨勢(shì)，從而能有效抑制端點(diǎn)效應(yīng)的產(chǎn)生。

2.3 波形相似度定義

波形相似度為表明兩波形相似程度的度量。對(duì)信號(hào)x(t)，y(t),(t=1,2…n)，定義兩波形相似度為

(5)

式中：cov(·)為x(t)與y(t)協(xié)方差；σ(x)，σ(y)分別為x(t)與y(t)的方差。

定理1：當(dāng)且僅當(dāng)隨機(jī)變量Y與X之間存在線性關(guān)系Y=a+bX時(shí)，相關(guān)系數(shù)|ρ(X,Y)|=1。

據(jù)定理1知：兩波形相似度較高時(shí)(ρ≥0.8)，可近似將兩者視為線性關(guān)系，用線性回歸方法尋找兩者間線性關(guān)系。

2.4 延拓步驟

圖1 待匹配波示意圖

對(duì)原始信號(hào)延拓時(shí)，左右兩端方法相同，以左端點(diǎn)為例說(shuō)明延拓步驟：

(1) 找出待匹配的波，并確定其長(zhǎng)度。設(shè)待匹配的波為ω，起點(diǎn)為左端點(diǎn)值，向右包括極大值與極小值。當(dāng)?shù)谝粋€(gè)極大值距左端點(diǎn)距離Lmax1較第一個(gè)極小值距左端點(diǎn)距離Lmin1小時(shí)，則待匹配波ω的終點(diǎn)為第一個(gè)極小值點(diǎn)min1，反之則為第一個(gè)極大值點(diǎn)max1。待匹配波ω終點(diǎn)為min1情況見(jiàn)圖1。

(2) 計(jì)算整個(gè)數(shù)據(jù)段中各波段與ω的相似度ρ。為提高運(yùn)算速度，僅計(jì)算極值點(diǎn)附近波的相似度。

(3) 找出相似度最高的子波段ω_match，即相似度|ρ|最大的波。

(4) 用線性回歸方法計(jì)算兩波段數(shù)據(jù)關(guān)系，即

Y=a+bX

(6)

式中：X為ω幅值向量，Y為ω_match幅值向量。

(5) 分別向左延拓一個(gè)極大值(tmax0,vmax0)與一個(gè)極小值(tmin0,vmin0)。設(shè)在ω_match左端最近的極大值、極小值點(diǎn)為(tmaxp,vmaxp)與(tminp,vminp)。延拓的極大值max0與極小值min0的幅值據(jù)式(6)計(jì)算：

vmax0=(maxp-a)/b

(7)

vmin0=(minp-a)/b

(8)

延拓極值對(duì)應(yīng)時(shí)刻由極值距端點(diǎn)距離確定。設(shè)(tmaxp,maxp)，(tminp,minp)距離ω_match左端點(diǎn)距離分別為L(zhǎng)maxp，Lminp，則延拓的極大值、極小值對(duì)應(yīng)時(shí)刻分別為

tmax0=-LmaxP

(9)

tmin0=-LminP

(10)

3 改進(jìn)的EMD分解停止準(zhǔn)則

分解停止準(zhǔn)則為決定EMD分解停止的標(biāo)準(zhǔn)。終止條件太苛刻時(shí)會(huì)產(chǎn)生過(guò)多IMF分量，會(huì)干擾對(duì)信號(hào)的正確判斷；終止條件寬松時(shí)則會(huì)使分解不夠細(xì)致，無(wú)法識(shí)別信號(hào)的變化規(guī)律。對(duì)整個(gè)EMD分解過(guò)程停止需有適當(dāng)?shù)呐袛嘁罁?jù)。通常采用兩個(gè)判斷標(biāo)準(zhǔn)，即當(dāng)最后一個(gè)IMF或殘余分量幅值小于預(yù)設(shè)值時(shí)，整個(gè)分解過(guò)程停止；或當(dāng)殘余分量變成單調(diào)函數(shù)或常數(shù)時(shí)，不能再篩分出IMF時(shí)停止。用主元分析方法(PCA)的累計(jì)方差貢獻(xiàn)率對(duì)分解停止時(shí)機(jī)確定[11]時(shí)需對(duì)IMF及殘余分量進(jìn)行主元分析。即將原始變量投影到新坐標(biāo)系，使新主元互不相關(guān)，從而可用較少維度代表原始變量特征。而用EMD方法所得IMF已具有近似正交性及不相關(guān)性，因此可直接進(jìn)行累計(jì)方差貢獻(xiàn)率計(jì)算。

3.1 IMF近似正交性

由EMD分解過(guò)程知，正交性在實(shí)際意義上是滿足。但理論上尚未給出嚴(yán)格的數(shù)學(xué)證明。通過(guò)EMD分解所得各IMF在局部均應(yīng)相互正交，原因?yàn)槊總€(gè)IMF均由原信號(hào)與其極大、極小值包絡(luò)的局部均值之差獲得，故有

(11)

式(11)嚴(yán)格意義上講不準(zhǔn)確，因均值由信號(hào)極大值、極小值點(diǎn)由曲線擬合求得，與真實(shí)均值存在一定誤差，雖泄露不可避免，但較小。IMF正交性可通過(guò)后驗(yàn)方法給出。原信號(hào)重構(gòu)信號(hào)可表示為

(12)

將殘余分量rn(t)視為第n+1個(gè)IMF分量，對(duì)上式兩邊平方得：

(13)

因此對(duì)原信號(hào)x(t)的正交性指標(biāo)(Index of Orthogonality)可定義[2]為

(14)

式中：T為信號(hào)總長(zhǎng)度。

若分解是正交的，則式(14)中交叉項(xiàng)為0。通過(guò)大量試驗(yàn)驗(yàn)證[2-3]得出，對(duì)一般信號(hào)其正交性指標(biāo)通常不超過(guò)1%，對(duì)短信號(hào)在極端情況下可能達(dá)5%，即EMD分解所得各IMF為近似正交。

3.2 基于累計(jì)方差貢獻(xiàn)率的分解停止準(zhǔn)則

方差貢獻(xiàn)率概念源于主元分析方法。方差貢獻(xiàn)率及累計(jì)方差貢獻(xiàn)率為衡量主成分重要性指標(biāo)。在改進(jìn)的EMD方法中，定義第i個(gè)IMF方差貢獻(xiàn)率為

(15)

該值越大，說(shuō)明第i個(gè)IMF所含信息量越多。其中σ(·)為信號(hào)方差。

累計(jì)方差貢獻(xiàn)率(Cumulative Percent Variance，CPV)為前k個(gè)IMF方差貢獻(xiàn)率之和。實(shí)際應(yīng)用時(shí)需大于一個(gè)期望閾值，從而可更好反映變量波動(dòng)，即

(16)

式中：CL為確定閾值。

據(jù)各IMF分量的近似正交性，每個(gè)IMF分量攜帶信息量的大小可用方差貢獻(xiàn)率衡量。由累計(jì)方差貢獻(xiàn)率確定分解停止時(shí)機(jī)。計(jì)算每次分解IMF分量imf1,imf2,…imfn及殘余分量rn的方差貢獻(xiàn)率，求得每個(gè)分量及殘余分量的貢獻(xiàn)率。當(dāng)殘余分量貢獻(xiàn)率小于給定閾值時(shí)則停止分解。分量方差貢獻(xiàn)率小于4×10-3時(shí)不必再分解[11]。因此本文采用兩個(gè)分解停止判斷標(biāo)準(zhǔn)，即累積方差貢獻(xiàn)率大于設(shè)定的閾值時(shí)，整個(gè)分解過(guò)程停止；或殘余分量變成單調(diào)函數(shù)或常數(shù)時(shí)，不能再篩分出IMF時(shí)停止。

4 仿真實(shí)例

4.1 多種類型系統(tǒng)誤差疊加

系統(tǒng)誤差表現(xiàn)形式可為單一常值誤差、線性漂移誤差或周期性誤差，也可為其組合。本文以兩種形式誤差疊加說(shuō)明改進(jìn)EMD方法的有效性。設(shè)系統(tǒng)誤差有三個(gè)誤差分量構(gòu)成，其中y1(t)為線性誤差(常值誤差可視為線性誤差的特殊形式)，y2(t),y3(t)為周期性誤差。各誤差數(shù)學(xué)表達(dá)式為

y(t)=y1(t)+y2(t)+y3(t)

(17)

式中：y1(t)=ax+b,a=0.5,b=1;y2(t)=2sin(2πt);y3(t)=sin(10πt),(t=0:0.001:3)

為驗(yàn)證改進(jìn)EMD方法的有效性，分別對(duì)上述信號(hào)直接進(jìn)行EMD分解、鏡像延拓、本文方法分解。采用基于累積方差貢獻(xiàn)率的停止準(zhǔn)則，所得結(jié)果見(jiàn)圖2、圖3、圖4(其中實(shí)線為IMF分量，虛線為原始誤差分量)。

由三圖及表1看出，基于累積方差貢獻(xiàn)的停止準(zhǔn)則下，三種方式均能獲得兩個(gè)IMF分量與一個(gè)殘余分量。未對(duì)端點(diǎn)進(jìn)行延拓時(shí)的EMD分解所得分量與原始分量誤差較小，三個(gè)分量中線性誤差最大，此因端點(diǎn)效應(yīng)向里蔓延之結(jié)果；基于鏡像延拓的分解方法，由于鏡像點(diǎn)選擇不合理原因產(chǎn)生一虛假分量，誤差較大，尤其線性誤差淹沒(méi)其中；基于波形相似度的端點(diǎn)延拓方法所得三個(gè)分量在端點(diǎn)處分解效果較好，與各原始信號(hào)誤差最小。

圖2 EMD分解

表1 各分解方法均方差

4.2 隨機(jī)誤差與系統(tǒng)誤差共存

實(shí)際所測(cè)系統(tǒng)中，誤差可由測(cè)量值減去真值獲得，誤差成分必含隨機(jī)誤差與系統(tǒng)誤差。進(jìn)行系統(tǒng)誤差分離時(shí)，需對(duì)隨機(jī)誤差進(jìn)行處理。本節(jié)仿真信號(hào)在式(17)基礎(chǔ)上加入N(0,0.62) 的正態(tài)白噪聲獲得。仿真時(shí)分別對(duì)信號(hào)進(jìn)行直接分解及降噪后分解，實(shí)驗(yàn)結(jié)果見(jiàn)圖5、圖6。圖5為未對(duì)隨機(jī)誤差處理的結(jié)果，分解時(shí)因隨機(jī)誤差頻率較高，產(chǎn)生較多的IMF高頻分量，使分解過(guò)程受到干擾，從而未很好分離出系統(tǒng)誤差。圖6先對(duì)原始信號(hào)進(jìn)行隨機(jī)誤差處理，再進(jìn)行基于波形相似度的延拓分解，不但較完整的分離出各項(xiàng)系統(tǒng)誤差，且端點(diǎn)效應(yīng)亦得到較好抑制；但由于隨機(jī)誤差干擾，端點(diǎn)效應(yīng)未能完全消除。

圖5 直接進(jìn)行EMD分解

圖6 降噪后基于波形相似度延拓方法分解

5 實(shí)例分析

本文用動(dòng)態(tài)測(cè)角儀誤差數(shù)據(jù)[12]驗(yàn)證改進(jìn)EMD方法的分離效果?？紤]采樣點(diǎn)數(shù)較少，對(duì)誤差數(shù)據(jù)先進(jìn)行三次樣條插值，所得誤差曲線見(jiàn)圖7。利用改進(jìn)EMD方法分解所得結(jié)果見(jiàn)圖8。由二圖中看出，原始系統(tǒng)誤差包括線性漂移誤差與周期性誤差。經(jīng)改進(jìn)EMD方法分離得出4個(gè)誤差成分，其中imf1,imf2,imf3為具有周期性誤差諧波，所得殘余分量res為線性漂移誤差。系統(tǒng)誤差得到較好分離，可為測(cè)量?jī)x器的誤差分析提供數(shù)據(jù)支持。

圖7 動(dòng)態(tài)測(cè)角儀誤差圖

圖8 改進(jìn)EMD方法分解

6 結(jié) 論

由于EMD方法為新的數(shù)據(jù)處理方法，不但可處理線性、平穩(wěn)數(shù)據(jù)，且可自適應(yīng)處理非線性及非平穩(wěn)數(shù)據(jù)。本文將改進(jìn)的EMD方法用于系統(tǒng)誤差分離過(guò)程。仿真實(shí)驗(yàn)、實(shí)例分析表明，基于波形相似度的端點(diǎn)延拓方法可有效抑制端點(diǎn)效應(yīng)，累積方差貢獻(xiàn)率可適當(dāng)確定分解停止時(shí)機(jī)，從而較好完成系統(tǒng)誤差的分離，具有一定應(yīng)用前景。

[1] 李世平，付宇，張進(jìn).一種基于EMD的系統(tǒng)誤差分離方法[J].中國(guó)測(cè)試，2011, 37(3): 9-13.

LI Shi-ping，F(xiàn)U Yu，ZHANG Jin.Systematic error separation method based on EMD[J].China Measurement & Test, 2011，37(3): 9-13.

[2] Huang N E, Shen Z, Long S R, et al.The empirical mode decomposition method and the Hilbert spectrum for non-stationary time series analysis[C].Proc.R.Soc.London, U.K., 1998: 903-995.

[3] Huang N E.Review of empirical mode decomposition[C].Proceedings of SPIE.Orlando, USA, 2001: 71-80.

[4] Rilling G, Flandrin P,Goncalves P.On empirical mode decomposition and its algorithms[C].IEEE-EURASIP Workshop on Nonlinear Signal and Image Processing, Grado, Italy, 2003:811-815.

[5] Zhang J,Yan R Q, Gao R X, et al.Performance enhancement of ensemble empirical mode decomposition[J].Mechanical Systems and Signal Processing, 2010(24):2104-2123.

[6] 王婷.EMD算法研究及其在信號(hào)去噪中的應(yīng)用[D].哈爾濱:哈爾濱工程大學(xué),2010.

[7] 邵晨曦，王劍，范金鋒，等.一種自適應(yīng)的EMD端點(diǎn)延拓方法[J].電子學(xué)報(bào),2007,35(10):1944-1948.

SHAO Chen-xi, WANG Jian, FAN Jin-feng, et al.A self- adaptive method dealing with the end issue of EMD[J].Acta Electronica Sinica, 2007, 35(10):1944-1948.

[8] 蔡艷平，李艾華，石林鎖，等.EMD端點(diǎn)效應(yīng)的改進(jìn)型混沌延拓方法及其在機(jī)械故障診斷中的應(yīng)用[J].振動(dòng)與沖擊, 2011, 30(11): 46-52.

CAI Yan-ping, LI Ai-hua, SHI Lin-suo, et al.Processing method for end effects of EMD based on improved chaos forecasting model and its application in machinery fault diagnosis[J].Journal of Vibration and Shock, 2011, 30(11): 46-52.

[9] 許寶杰,張建民,徐小力,等.抑制EMD端點(diǎn)效應(yīng)方法的研究[J].北京理工大學(xué)學(xué)報(bào)，2006, 26(3):196-200.

XU Bao-jie, ZHANG Jian-min, XU Xiao-li, et al.A study on the method of restraining the ending effect of empirical mode decomposition[J].Transactions of Beijing Institute of Technology, 2006, 26(3): 196-200.

[10] 郭明威,倪世宏,朱家海，等.振動(dòng)信號(hào)中的HHT/EMD端點(diǎn)延拓方法研究[J].振動(dòng)與沖擊, 2012, 31(8): 62-66.

GUO Ming-wei, NI Shi-hong, ZHU Jia-hai, et al.HHT/EMD end extension method in vibration signal analysis[J].Journal of Vibration and Shock, 2012, 31(8): 62-66.

[11] 李雪耀，鄒曉杰，張汝波，等.譜熵和主成分分析用于EMD分解研究[J].哈爾濱工程大學(xué)學(xué)報(bào),2009,30(7):797-803.

LI Xue-yao, ZOU Xiao-jie, ZHANG Ru-bo, et al.Research on empirical mode decomposition based on spectrum entropy methods and principal component analysis[J].Journal of Harbin Engineering University, 2009, 30(7):797-803.

[12] 費(fèi)業(yè)泰，盧榮勝.動(dòng)態(tài)測(cè)量誤差修正原理與技術(shù)[M].北京：中國(guó)計(jì)量出版社,2001.