孟 明,楊國(guó)雨,高云園,甘海濤,羅志增

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

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

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

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

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

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

1 原理和算法

1.1 EEMD算法

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

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

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

(1)

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

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

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

(2)

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

1.2 DSS算法原理

1.2.1 盲源分離基本模型

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

(3)

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

X=AS+N

(4)

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

Y=WX

(5)

1.2.2DSS理論框架

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

步驟1 求源信號(hào)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ù)前一次得到的源信號(hào)S與混合矩陣A對(duì)本次觀測(cè)信號(hào)的概率進(jìn)行估計(jì),p(S|A,X)為源信號(hào)的后驗(yàn)概率,p(S)為前一次源信號(hào)S的先驗(yàn)概率。

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

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

(7)

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

(8)

(9)

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

s=wTX

(10)

s+=f(s)

(11)

w+=Xs+T

(12)

(13)

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

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

1.3 近似熵

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

近似熵的具體計(jì)算步驟為:

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

步驟2 計(jì)算時(shí)間序列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)]為兩者對(duì)應(yīng)元素中差值最大的一個(gè),即:

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,對(duì)每一個(gè)i值統(tǒng)計(jì)d[H(i),H(j)]小于r的數(shù)目及此數(shù)目與距離總數(shù)N-m+1的比值,即:

(15)

步驟6 定義φm(r):

(16)

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

步驟8 計(jì)算近似熵ApEn:

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

(17)

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

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

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

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

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

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

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

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

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

2 實(shí)驗(yàn)結(jié)構(gòu)與分析

2.1 仿真實(shí)驗(yàn)

(18)

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

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

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

圖2 EEMD分解加噪信號(hào)后的IMF分量

圖3 DSS各獨(dú)立源分量

圖4 DSS各獨(dú)立源分量頻譜圖

表1 DSS各獨(dú)立分量頻譜近似熵和時(shí)域近似熵

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

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

(19)

(20)

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

圖8 原始C3通道EEG分解后IMF分量時(shí)域圖

2.2 EEG消噪實(shí)驗(yàn)

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

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

實(shí)驗(yàn)選擇C3通道運(yùn)動(dòng)想象部分的腦電信號(hào)作為研究對(duì)象,其波形如圖7所示。

圖6 實(shí)驗(yàn)設(shè)計(jì)流程時(shí)序圖

圖7 C3通道運(yùn)動(dòng)想象部分腦電信號(hào)

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

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

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

圖10 3種消噪方法頻譜對(duì)比圖

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

3 結(jié)論

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