基于EEMD與DSS-ApEn的腦電信號消噪方法*

孟 明,楊國雨,高云園,甘海濤,羅志增

(杭州電子科技大學(xué)智能控制與機器人研究所,杭州 310018)

頭皮腦電EEG(Electroencephalogram)由于其易采集、無創(chuàng)性以及較好的時間分辨率,在最近被廣泛關(guān)注的腦機接口應(yīng)用中,有著不可替代的作用。但由于頭皮腦電是一種非平穩(wěn)非線性極其微弱的隨機信號[1],極易被心電,眼電,電磁干擾,工頻干擾等大量外界干擾信號淹沒[2],所以消噪就成了腦電信號預(yù)處理階段中最重要的步驟之一。

文獻[3]提出的經(jīng)驗?zāi)B(tài)分解EMD(Empirical Mode Decomposition)是一種自適應(yīng)的數(shù)據(jù)分解方法,能夠較好地處理隨機非平穩(wěn)信號,具有良好的自適應(yīng)性、信號的完備性等[4]。因此,EMD成為現(xiàn)代信號處理分析領(lǐng)域的一個研究熱點,然而EMD方法存在模態(tài)混疊現(xiàn)象:當(dāng)信號的極值點分布不均勻時,會導(dǎo)致一個相似尺度信號分布在不同的內(nèi)蘊模態(tài)函數(shù)IMF(Intrinsic Mode Function)分量中或者一個IMF分量中出現(xiàn)多個尺度信號,從而使IMF分量失去了原本的物理意義。

Wu和Huang在EMD基礎(chǔ)上改進提出了集合經(jīng)驗?zāi)B(tài)分解[5]EEMD(Ensemble Empirical Mode Decomposition),通過加入白噪聲改變信號的極值點分布,進而有效地改善EMD的模態(tài)混疊現(xiàn)象,極大提高了信號的信噪比。但是基于EEMD的信號消噪方法通過直接去除EEMD分解結(jié)果中的前幾個高頻IMF分量來消噪,導(dǎo)致在去掉高頻IMF分量中噪聲的同時去除了其中的有效信息成分[6]。

降噪源分離DSS(De-noising Source Separation)是盲信號分離廣泛應(yīng)用的一種技術(shù)。Valpola等人在研究基于貝葉斯準(zhǔn)則的快速獨立分量分析算法的過程中,提出了DSS的基本理論概念[7]。DSS本質(zhì)是在源信號復(fù)雜結(jié)構(gòu)未知的的條件下,根據(jù)信號的統(tǒng)計特性,選擇合適的降噪函數(shù)將復(fù)雜信號分解為若干組成分量,從而實現(xiàn)源信息的挖掘。文獻[8]利用基于空間濾波的DSS方法分析處理腦電信號,有效地降低了刺激誘發(fā)反應(yīng)范式中的噪聲。王元生等[9]把DSS應(yīng)用于機械工程領(lǐng)域,有效地提取出轉(zhuǎn)子振動故障信號特征。但DSS得到的各獨立源分量的順序不確定,文獻[8]通過人工觀察區(qū)分獨立源分量,文獻[9]通過計算源信號和消噪信號的相關(guān)系數(shù)來區(qū)分。以上區(qū)分方法存在消噪過程中不能直接確定獨立源分量是噪聲信號還是有效信號的缺陷。

本文在使用近似熵ApEn(Approximate Entropy)區(qū)分DSS結(jié)果中各獨立源信號是有效信號還是噪聲信號的基礎(chǔ)上,提出一種基于EEMD與DSS-ApEn相結(jié)合的腦電信號消噪方法,并與獨立EEMD消噪方法以及基于EEMD與改進提升小波的消噪方法進行比較。

1 原理和算法

1.1 EEMD算法

EEMD主要利用了白噪聲在各個頻段具有能量一致的統(tǒng)計特性來解決EMD的模態(tài)混疊現(xiàn)象[10],其主要原理是:在不連續(xù)的信號中,低頻分量的極值點間隔分布稀疏、而高頻分量的則分布密集。當(dāng)在待處理信號中加入白噪聲,信號的低頻分量的極值點分布發(fā)生改變,整個頻帶中的極值點間隔分布均勻,保證每個固有模態(tài)函數(shù)時域的連續(xù)性,從而避免了模態(tài)混淆。EEMD分解的具體步驟[11]為:

步驟1 往待處理信號s(t)中添加均值為0、標(biāo)準(zhǔn)差為常數(shù)的白噪聲n1(t),得到待分解的信號s1(t)。

步驟2 EMD算法分解信號s1(t),得到n個IMF分量imfi1(t)和余項rn1(t),即:

(1)

式中,n為分解層數(shù)。

步驟3 重復(fù)以上兩個步驟T-1次,但每次都加入同分布的不同白噪聲。

步驟4 對相應(yīng)的IMF分量進行總體平均運算處理,以消除白噪聲對IMF分量的影響,最終得到各IMF分量為:

(2)

式中,imfi(t)為總體平均后的第i個IMF分量;T為總體平均次數(shù);imfij(t)為第j次加入白噪聲后獲得的第i個IMF分量。

1.2 DSS算法原理

1.2.1 盲源分離基本模型

設(shè)n個信號源S={s1(t),s2(t),…,sn(t)}發(fā)出的信號在m個不同位置觀測到的混合信號分別為X={x1(t),x2(t),…,xm(t)},即:

(3)

式中,aij為混合系數(shù),ni為第i個源信號所對應(yīng)的噪聲,則可以將上式轉(zhuǎn)換為矩陣形式:

X=AS+N

(4)

因此,源分離處理就是要在源信號S和混合系數(shù)矩陣A均未知的情況下,求解分離矩陣W,使得在觀測信號中分離出對源信號的估計Y=[y1(t),y2(t),…,yn(t)],即:

Y=WX

(5)

1.2.2DSS理論框架

期望值最大化EM(Expectation-Maximization)算法是經(jīng)典的盲源分離算法之一,其具體的算法可以分為兩步:

步驟1 求源信號S期望值

q(S)=p(S|A,X)=p(X|A,S)p(S)/p(X|A)

(6)

式中,p(X|A)則是在已知混合矩陣情況下的條件概率,p(X|A,S)為根據(jù)前一次得到的源信號S與混合矩陣A對本次觀測信號的概率進行估計,p(S|A,X)為源信號的后驗概率,p(S)為前一次源信號S的先驗概率。

步驟2 尋求期望值中的最大值

利用式(6)中求得的期望值,就可以根據(jù)EM算法的求得新的混合矩陣Anew,即:

(7)

式中,CXS和CSS是根據(jù)期望函數(shù)q(S)計算得到的協(xié)方差矩陣,即:

(8)

(9)

在EM算法的迭代求解過程中,所有的獨立成分分量按照統(tǒng)計特性被一次性求解得到。相對EM算法中分量被同時求出,Hyv?rinen等[12]在2001年提出利用預(yù)白化的方法逐次提取各獨立成分源信號進行迭代求解,充分體現(xiàn)了組成成分的獨立和重要性,這是提出DSS理論的重要理論基礎(chǔ)。基本求解過程框架如下:

s=wTX

(10)

s+=f(s)

(11)

w+=Xs+T

(12)

(13)

式(10)通過似然估計函數(shù)計算源信號的噪聲估計,式(11)中s+對應(yīng)EM算法中對期望q(S)的估計,即為降噪過程。所有含有式(11)降噪環(huán)節(jié)的源分離方法被定義為DSS。式(12)為利用降噪處理后的源信號s+對分離矩陣w的重新估計。式(13)完成歸一化。DSS對每個源信號的分離過程,就是在給定式(10)中的wT之后,通過不斷迭代式(11)～式(13)求得每個源信號的分離信號。

DSS理論的關(guān)鍵在于根據(jù)待處理的信號特點,選擇或者設(shè)計適合的降噪函數(shù)。文獻[13]模擬實驗分析的結(jié)果表明,基于正切降噪函數(shù)的源分離方法從非線性耦合信號中提取的分量與源信號時域相關(guān)系數(shù)最高,所以本文選取正切降噪函數(shù)f(s)=s-tanh(s)作為DSS的降噪函數(shù)。

1.3 近似熵

DSS結(jié)果中的各獨立源分量可以看作各信號源,需要找出含有噪聲的獨立信號源。在信息論中,熵指數(shù)能夠表征系統(tǒng)的復(fù)雜程度。近似熵[14]可以衡量時間序列復(fù)雜度,用于度量信號產(chǎn)生新模式的概率,越復(fù)雜的時間序列對應(yīng)的近似熵越大,故本文采用近似熵作為腦電分量與噪聲分量的判據(jù)。

近似熵的具體計算步驟為:

步驟1 對系統(tǒng)進行等間隔采樣,得到時間序列u(1),u(2),…,u(N);

步驟2 計算時間序列u(i)的標(biāo)準(zhǔn)差SD;

步驟3 構(gòu)造一組m維矢量:H(1),H(2),…,H(N-m+1),其中H(i)=[u(i),u(i+1),…,u(i+m-1)];

步驟4 定義H(i)與H(j)間的距離d[H(i),H(j)]為兩者對應(yīng)元素中差值最大的一個,即:

d[H(i),H(j)]=max[|u(i+k)-u(j+k)|],1≤i,j≤N-m+1,k=0,…,m-1

(14)

步驟5 給定閾值r,對每一個i值統(tǒng)計d[H(i),H(j)]小于r的數(shù)目及此數(shù)目與距離總數(shù)N-m+1的比值,即:

(15)

步驟6 定義φm(r):

(16)

步驟7 將比較序列長度參數(shù)m加1,即m=m+1,重復(fù)步驟3～7,得到φm+1(r);

步驟8 計算近似熵ApEn:

ApEn=φm(r)-φm+1(r)

(17)

ApEn的值顯然與m,r的取值有關(guān)。Pincus[15]根據(jù)實踐,建議取值m=2,r=0.1×SD。

1.4 基于EEMD與DSS-ApEn相結(jié)合的消噪方法

采用EEMD對信號進行分解,得到多個IMF分量,前幾個分量主要包含信號中的高頻成分,后面的分量主要包含信號中的低頻成分,噪聲強度隨著IMF分量層次的增加也會越來越弱,其中,腦電信號的低頻成分以有效信息為主,而高頻成分含大量的噪聲[16]。EEMD消噪方法通過直接去除EEMD分解中的前幾個高頻IMF分量來消噪,但是這樣會導(dǎo)致消除高頻中的噪聲成分的同時亦剔除了其中的有效信息成分。

DSS能有效的提取各獨立源分量,但是DSS結(jié)果中各獨立源分量的順序不確定,不能直接確定獨立源分量對應(yīng)的信號源是噪聲信號還是腦電信號。由于噪聲信號頻率的復(fù)雜程度比腦電信號的要大很多,故采用近似熵來判斷獨立成分中的噪聲信號。

針對以上情況,本文將EEMD與DSS-ApEn相結(jié)合,其具體消噪步驟如下。

步驟1 應(yīng)用EEMD將待消噪信號分解為IMF集,去掉最高頻IMF分量得到新的IMF集Q;

步驟2 對步驟1中新的IMF集Q應(yīng)用DSS算法,得出一組獨立源分量S,并求出各獨立源分量的頻譜;

步驟3 計算步驟2中頻譜的近似熵,選擇近似熵最大的獨立源分量作為噪聲信號濾除,得到一組新的獨立源分量S′;

步驟4 將混合矩陣A乘以獨立源分量S′變換出新的IMF集,并將該IMF集累加重構(gòu),得到消噪后的腦電信號。

2 實驗結(jié)構(gòu)與分析

2.1 仿真實驗

(18)

標(biāo)準(zhǔn)信號和加入高斯白噪聲后得到的信號如圖1 所示,加噪信號的信噪比為1 dB,圖2是EEMD分解加噪信號后的IMF分量。

按照1.4小節(jié)中的步驟2進行處理,去掉最高頻IMF分量,對剩余的IMF應(yīng)用DSS得到各獨立源分量如圖3所示,并求出各獨立源分量頻譜如圖4所示,各獨立源分量頻譜近似熵和時域近似熵如表1 所示。

圖1 標(biāo)準(zhǔn)信號和信噪比為1 dB時加噪標(biāo)準(zhǔn)信號

圖2 EEMD分解加噪信號后的IMF分量

圖3 DSS各獨立源分量

圖4 DSS各獨立源分量頻譜圖

表1 DSS各獨立分量頻譜近似熵和時域近似熵

從圖4中可以看出IC6是噪聲信號,但是觀察表1中各獨立源分量的頻譜近似熵和時域近似熵,頻譜近似熵可以很好的區(qū)分出噪聲信號,而時域近似熵的區(qū)分度不是很大,故本文采用頻譜近似熵區(qū)別各獨立源信號。表1中IC6的頻譜近似熵最大,故第6個獨立源分量是噪聲信號,濾除噪聲信號后得到一組新的獨立源分量。按照1.4小節(jié)中的步驟4進行重構(gòu),即可得出消噪后的信號。本文消噪方法與獨立EEMD分解后直接舍棄高頻IMF分量的方法、基于EEMD與改進提升小波的消噪方法[18]的消噪效果如圖5所示。

使用信噪比SNR(Signal to Noise Ratio)和均方根誤差RMSE(Mean Squared Error)作為評價加噪信號去噪效果指標(biāo),分別定義為:

(19)

(20)

表2 3種消噪方法加了不同信噪比噪聲的仿真信號的實驗數(shù)據(jù)

圖8 原始C3通道EEG分解后IMF分量時域圖

2.2 EEG消噪實驗

本文進一步采用實際采集的腦電信號數(shù)據(jù)進行消噪效果驗證。應(yīng)用美國NeuroScan公司的scan4.3腦電采集系統(tǒng)作為腦電信號的采集設(shè)備,采集C3、C4和Cz 3個電極處的腦電信號,采樣頻率為250 Hz,精度為32 bit。受試者為身體健康的在校研究生。腦電電極按照國際標(biāo)準(zhǔn)導(dǎo)聯(lián)10-20系統(tǒng)放置,以左側(cè)乳突為參考電極,右側(cè)乳突為接地電極。

本實驗通過采集想象左手、右手所產(chǎn)生的腦電信號進行消噪處理。每次實驗持續(xù)9 s。0～2 s時,受試者保持身體放松狀態(tài);2 s時,系統(tǒng)發(fā)出“滴”提示音并且在顯示屏中間顯示一個十字光標(biāo),提示受試者集中注意力;3 s～6 s時,受試者根據(jù)屏幕上隨機出現(xiàn)的左(←)、右(→)兩個方向的誘導(dǎo)圖片分別想象左手揮動、右手揮動這兩種運動,直至第6 s時系統(tǒng)發(fā)出“嘟”的提示音,受試者結(jié)束想象。一次完整實驗的采樣時序如圖6所示。

實驗選擇C3通道運動想象部分的腦電信號作為研究對象,其波形如圖7所示。

圖6 實驗設(shè)計流程時序圖

圖7 C3通道運動想象部分腦電信號

圖7中,橫坐標(biāo)為信號的采樣點個數(shù),縱坐標(biāo)為腦電信號的幅值。在原始C3通道EEG經(jīng)過EEMD分解后,得到8個IMF分量和一個余量res,所得IMF分量的時域圖如圖8所示。從圖8可以看出,各階IMF分量包含不同的時間特征尺度,且隨著階次的增加,頻率越小。

按照本文1.4小節(jié)中的消噪步驟進行消噪。本文消噪方法與獨立EEMD分解后直接舍棄高頻IMF分量的方法、基于EEMD與改進提升小波的消噪方法的消噪時域圖如圖9所示、頻域圖如圖10所示。

圖9 3種消噪方法時域?qū)Ρ葓D

圖10 3種消噪方法頻譜對比圖

結(jié)合圖9和圖10所示,雖然經(jīng)前兩者算法處理后的信號噪聲得到了明顯的消除,但是EEMD分解消噪方法得到的信號毛刺較多,且30 Hz以上的高頻部分還有保留,而基于EEMD與改進提升小波消噪方法消噪后所得的信號過于平滑,頻率在20 Hz～30 Hz細節(jié)信號濾除過多,發(fā)生了信號失真,部分細節(jié)信號被當(dāng)做噪聲消除了。而采用基于EEMD與DSS-ApEn相結(jié)合的消噪方法消噪后的腦電信號波形相對清晰,更重要的是原始信號的細節(jié)特征也被很好地保留下來。這是由于EEMD分解后僅對噪聲主導(dǎo)的高頻IMF分量進行消噪,從而有效保留了有用的腦電信號細節(jié)信息。DSS在迭代求解獨立源分量的過程中,將利用前一次的求解結(jié)果作為后一次求解的先驗知識,具有很強的自適應(yīng)性,所以DSS可以更好地從高頻IMF分量中分離出細節(jié)信息,達到更好保留細節(jié)信息的目的。

3 結(jié)論

本文提出一種基于EEMD與DSS-ApEn相結(jié)合的消噪方法。對EEMD分解得到的IMF分量處理方法進行改進,將濾除最高頻分量后的IMF作為DSS的輸入。對DSS得到的各獨立源分量,選擇頻譜近似熵最大的獨立源分量作為噪聲信號濾除。采用不同消噪方法的仿真和實際腦電信號消噪實驗結(jié)果表明,本文的消噪方法具有更強的噪聲抑制能力。良好的消噪效果和保留信號細節(jié)的能力,為后續(xù)腦電信號的特征提取與模式識別打下堅實的基礎(chǔ)。