基于KICA和MAMA-EMD的眼電偽跡去除方法

張圓圓 孫炎珺 李明愛,2,3

0 引言

腦電信號(electroencephalogram，EEG)是由電極檢測到并能反映神經(jīng)元活動的生物電信號，因其具有低成本、可移植和高時間分辨率等優(yōu)點被廣泛應(yīng)用于腦疾病診斷和腦機接口(brain-computer interface，BCI)系統(tǒng)中[1-2]。然而，在實際采集過程中，EEG常受到各種非大腦活動的偽跡干擾，其中眼電(electrooculogram，EOG)是EEG中最主要的干擾成分，它往往以大脈沖或尖峰的形式出現(xiàn)在EEG信號中，這為后續(xù)EEG的分析和應(yīng)用帶來了極大的困難。因此，如何在不丟失有用的神經(jīng)信息的前提下，準確分離并去除EEG中的EOG偽跡具有重要的研究意義[3]。

盲源分離(blind source separation，BSS)是一種最為常用的EOG偽跡去除方法[4]，它首先將多通道EEG分離為若干個源信號分量，然后通過剔除EOG相關(guān)分量重構(gòu)無偽跡的EEG信號。常用的BSS方法包括獨立分量分析(independent component analysis，ICA)[5]、典型相關(guān)分析(canonical correlation analysis，CCA)[6]、核獨立分量分析(kernel independent component analysis，KICA)[7]等。不足的是，單一的BSS技術(shù)直接將檢測到的EOG相關(guān)分量置零，這難免會丟失一部分有用的腦電數(shù)據(jù)。因而，一些研究學(xué)者將盲源分離技術(shù)與信號分解方法相結(jié)合以進一步提高偽跡去除的準確性。在文獻[8]中，作者將ICA與離散小波變換(discrete wavelet transform，DWT)相結(jié)合，提出一種將ICA用于小波域的DWICA方法，用于去除腦電信號中的EOG成分。結(jié)果顯示，DWICA不僅收斂速度快，而且具有較好的魯棒性。Li等[9]提出一種基于FastKICA和DWT的EOG偽跡自動去除方法(FKD)，該方法首先對受污染的EEG信號進行FastKICA分離，識別與EOG相關(guān)的獨立分量(independent components，ICs)；然后利用DWT對其進行分解以提取其中的腦電數(shù)據(jù)。結(jié)果表明，F(xiàn)KD在線性和非線性混合模型中都具有較好的偽跡去除效果。然而，小波變換需要根據(jù)信號自身特性提前選擇最優(yōu)的母小波函數(shù)和分解層數(shù)，這在實際應(yīng)用中較為困難且缺乏自適應(yīng)性[10]。經(jīng)驗?zāi)B(tài)分解(empirical mode decomposition，EMD)是一種自適應(yīng)的信號分解方法，它可將時間序列分解為一組固有模態(tài)函數(shù)(intrinsic mode functions，IMFs)。在文獻[11]中提出一種ICA與EMD相結(jié)合的眼電偽跡去除方法，該方法不僅能夠自動去除EOG成分，同時有效地保留了腦電數(shù)據(jù)。Yang等[12]提出一種基于CCA和集成經(jīng)驗?zāi)B(tài)分解(ensemble EMD，EEMD)的方法(CCA-EEMD)，用于去除EEG信號中的EOG偽跡，結(jié)果顯示，CCA-EEMD相比于ICA和CCA具有明顯的優(yōu)勢。然而，尖峰噪聲是影響EMD分解性能的一個重要因素，它會導(dǎo)致模態(tài)分裂效應(yīng)，即一個有物理意義的成分被分解到多個IMF中，不準確地提取所需的IMF最終會影響EOG的去除效果。掩膜最小弧長經(jīng)驗?zāi)B(tài)分解(masking-aided minimum arclength EMD，MAMA-EMD)是最近提出的一種改進的EMD算法，它可將尖峰分離到第1個IMF中，從而減輕其對篩選過程的影響，提高后續(xù)IMF推導(dǎo)的準確性[13]。

針對眼電偽跡的尖峰狀干擾特性，本文提出一種基于KICA和MAMA-EMD的EOG偽跡自動去除方法，記為KICMME。利用KICA分離EOG相關(guān)的IC，并使用MAMA-EMD對其進行分解，剔除眼電相關(guān)的IMFs。在半模擬和真實腦電數(shù)據(jù)上進行實驗研究，驗證本文方法的有效性。

1 基于KICMME的眼電偽跡去除方法

基于KICMME的眼電偽跡去除方法，其主要步驟及涉及的基本原理如下。

1.1 基于KICA分離原始腦電信號

設(shè)X=[x1(t),x2(t),…,xL(t)]T∈RL×M為原始EEG信號,S=[s1(t),s2(t),…,sL(t)]T∈RL×M為未知源信號,M為采樣點數(shù)，L為導(dǎo)聯(lián)數(shù)。KICA的數(shù)學(xué)模型為:

X=AS

(1)

Y=WX

(2)

式中:A∈RL×L為未知的混合矩陣。與經(jīng)典的ICA算法不同，KICA使用核函數(shù)代替向量之間的內(nèi)積，并采用低階近似法(如Gram矩陣)最小化對比函數(shù)，求取分離矩陣W，因而具有更好的分離穩(wěn)定性。

基于KICA，將受污染的EEG信號分離為多個ICs，記為：Y=[y1(t),y2(t),…,yL(t)]T∈RL×M。

1.2 基于峰度確定EOG相關(guān)的ICs

峰度是用于描述信號分布陡緩程度的統(tǒng)計量，通常，與EOG相關(guān)的獨立分量的峰度值遠高于EEG分量。對于第i個IC，其峰度：

(3)

式中:mn(n=2,4)為數(shù)據(jù)的第n個中心矩。

mn=E[(x-m1)n]

(4)

式中：m1為平均值；E表示期望函數(shù)。

計算每個IC的峰度值，并設(shè)定閾值st，將峰度值高于st的ICs識別為與EOG相關(guān)的ICs，記作：EY=[ey1(t),ey2(t),…,eyu(t)]T∈Ru×M，u為EOG相關(guān)IC的個數(shù)。

1.3 基于MAMA-EMD分解EOG相關(guān)的IC

MAMA-EMD算法通過向原始信號中添加一對掩膜信號，并根據(jù)最小弧長準則調(diào)整尖峰極值的高度，從而在第一個IMF中分離出尖峰。

圖1為最小弧長準則的示意圖，對于一個尖峰極大值點/極小值點ri(tas,xas)∈G,G為尖峰極值的集合，其最優(yōu)替換極值點r′i(tas,x′as)可通過以下步驟計算。

圖1 最小弧長準則示意圖Figure 1 The diagram of minimum arclength criterion

(1) 尋找尖峰極值ri的鄰近極大值/極小值點zj(taj,xaj)(j=1,2,…,Nz)。

(2) 利用3次樣條插值在{zj}和點(tas,y)之間進行插值，構(gòu)建包絡(luò)線U(t|·)。

(3) 根據(jù)以下公式計算最優(yōu)替換極值x′as:

x′as=arg miny{F(U(t|yaj=xaj;yas=y))}

(5)

式中:

(6)

代表包絡(luò)線U(t|·)的長度。通常，基于尖峰極值點附近的10個極值點(即Nz=10)構(gòu)建上/下包絡(luò)線就足夠了。

依據(jù)文獻[13]，對于給定信號g(t)，MAMA-EMD算法的主要步驟可總結(jié)為以下4個步驟。

(1) 選擇合適的振幅aM和頻率fM，生成掩膜信號ω(t)=aMsin(2πfMt)。

(4) 從g(t)中減去c1，推導(dǎo)剩余IMF的步驟與EMD算法相同。

基于MAMA-EMD算法,對1.2節(jié)中EOG相關(guān)的IC:eyi(t)(1≤i≤u)進行分解，得到一組IMFs,記為：eyi(t)=[ci1(t),ci2(t),…,civ(t)],v表示IMF的個數(shù)。采用MAMA-EMD對信號進行分解消除了尖峰對于IMF推導(dǎo)過程的干擾，從而能夠提高EOG相關(guān)IMFs的識別準確性。

1.4 基于低頻功率占比識別EOG相關(guān)的IMFs

由于EOG的能量主要集中在低頻范圍內(nèi)(f<8 Hz)[14],本文通過計算各個IMF的低頻功率占比自動識別與EOG相關(guān)的IMF。低頻功率占比定義為：k=∑P(f<8)/∑P(f)，P為功率，設(shè)定適合的閾值kt，將k值高于kt的IMFs視為EOG偽跡并去除，保留剩余EEG相關(guān)的IMFs。

1.5 重構(gòu)EOG偽跡去除后的EEG信號

計算1.4節(jié)中EEG相關(guān)的IMFs之和，并用于替換對應(yīng)EOG相關(guān)的IC。最后，基于KICA逆變換將所有無眼電成分的ICs重構(gòu)為“純凈”的EEG信號。

2 實驗結(jié)果與分析

2.1 基于半模擬腦電數(shù)據(jù)的實驗研究

2.1.1 數(shù)據(jù)的構(gòu)建

實驗數(shù)據(jù)來源于BCI Competition IV 2b數(shù)據(jù)集中的受試者S1，包含了3個導(dǎo)聯(lián)(C3、Cz和C4)的EEG信號，采樣頻率為250 Hz，總次數(shù)為120。從該數(shù)據(jù)庫中隨機選取1導(dǎo)記錄的EOG信號，濾波至0～8 Hz，作為純凈EOG信號。根據(jù)EEG和EOG信號之間的雙向污染特性[15]，建立如下模型:

EEGcon(i)=EEGcle(i)+piEOGcle

(7)

EOGcon=EOGcle+∑qiEEGcle(i)

(8)

式中：EEGcle(i)和EOGcle分別表示第i導(dǎo)純凈的EEG信號和EOG信號；EEGcon(i)和EOGcon則代表被污染的第i導(dǎo)EEG和EOG信號，i={1,2,3}為導(dǎo)聯(lián)序號。pi為EOG對第i導(dǎo)的EEG信號的干擾系數(shù)，qi為第i導(dǎo)EEG對EOG信號的干擾系數(shù)。為了確保實驗的合理性和可靠性，干擾系數(shù)pi和qi均隨機生成，且pi,qi∈[0.01,0.3]。圖2展示了一組受污染的EEG和EOG數(shù)據(jù)?？梢钥闯觯? s和2 s附近，3個導(dǎo)聯(lián)的EEG信號都被EOG嚴重污染，產(chǎn)生2個尖峰干擾。

圖2 受污染的EEG和EOG信號Figure 2 The contaminated EEG and EOG signals

2.1.2 眼電偽跡的自動識別與去除

(1) EOG相關(guān)IC的確定。首先采用KICA算法將受污染的EEG和EOG信號分離為4個獨立分量IC1～IC4，如圖3所示。為了選擇合適的閾值，本文根據(jù)120次試驗的數(shù)據(jù)分別計算了4個IC的平均峰度，結(jié)果如表1所示，IC4的平均峰度最大，為12.010，且明顯高于其他的ICs。將臨界閾值st設(shè)為10；進而IC4被自動檢測為與EOG相關(guān)的獨立分量。

圖3 基于KICA分離得到的獨立分量IC1～IC4Figure 3 The IC1-IC4 separated by KICA

表1 基于120次試驗數(shù)據(jù)的4個IC的平均峰度Table 1 The average kurtosis of ICs based on 120 trials

(2) 基于MAMA-EMD分解眼電IC。本文使用了EMD和MAMA-EMD兩種算法對IC4進行分解，分解結(jié)果見圖4。從圖4(a)中可以看出，在1s和2s附近，由于兩個尖峰脈沖的影響，EMD分解中產(chǎn)生了模態(tài)分裂效應(yīng)，導(dǎo)致眼電產(chǎn)生的尖峰成分被混合在IMF2和IMF3中；而在圖4(b)中，MAMA-EMD算法通過最小弧長法則替換尖峰極值，在篩選的過程中將尖峰分離至IMF1中，從而緩解了模態(tài)分裂效應(yīng)，提高了后續(xù)的分解性能。

圖4 EMD和MAMA-EMD對IC4的分解結(jié)果Figure 4 The decomposition results of IC4 based on EMD and MAMA-EMD

(3) IMF的功率譜分析。將IC4分解為一組IMFs后，計算了各IMF的功率譜，結(jié)果見圖5。由于EMD推導(dǎo)出的IMF的頻段覆蓋范圍是隨著IMF階次的增加，由高到低排列的，因此，圖5中僅展示了前6個IMF的功率譜。從圖5(a)中可以看出，EMD分解得到的IMF2和IMF3在尖峰的影響下存在嚴重的高低頻成分混合現(xiàn)象，因此很難準確地識別和剔除與EOG相關(guān)的IMF。而在圖5(b)中，除了IMF2、IMF3和IMF4的功率分布在高頻段以外，其余IMF的功率都主要集中在低頻范圍內(nèi)，符合EOG的頻域特性。這說明MAMA-EMD算法能夠很好地消除尖峰對信號分解過程的影響，克服了模態(tài)分裂的問題。

圖5 基于EMD和MAMA-EMD得到的IMF1～IMF6的功率譜對比Figure 5 The comparison of power spectrum for IMF1-IMF6 derived by EMD and MAMA-EMD

(4) 眼電去除與EEG信號恢復(fù)。為了準確地識別出與EOG高度相關(guān)的IMF，本文計算了各IMF的低頻功率占比k，結(jié)果如圖6所示，圖中藍色虛線表示閾值kt=0.9。由圖6可見，除IMF2、IMF3和IMF4外，其余IMFs的k值均大于閾值kt，因此，將IMF2、IMF3和IMF4視為腦電數(shù)據(jù)予以保留，而其余IMFs被識別為EOG相關(guān)的IMFs，從信號中剔除。最后，基于MAMA-EMD和KICA算法的逆變換，重構(gòu)去除EOG偽跡后的EEG信號。圖7展示了1.5 ～2.5 s內(nèi)的純凈腦電數(shù)據(jù)以及使用KICMME方法重構(gòu)的EEG，二者具有很好的一致性，這表明3個導(dǎo)聯(lián)EEG中的眼電偽跡均得到較好的去除，并且基于KICMME恢復(fù)的EEG信號能夠有效保留有用的神經(jīng)信息。

圖6 各IMF的低頻功率占比Figure 6 The low frequency power proportion of each IMF

圖7 純凈的EEG和重構(gòu)的EEG的對比Figure 7 Comparison of pure and reconstructed EEG

2.1.3 基于均方誤差和信噪比的性能評估與方法比較

本節(jié)使用均方誤差(mean square error, MSE)和信噪比(signal-to-noise ratio, SNR)這兩項指標進一步評估KICMME在EOG偽跡去除上的性能。 對于每次試驗的腦電數(shù)據(jù)，均隨機生成10組干擾系數(shù)，共構(gòu)建1 200組受污染的EEG和EOG信號，計算1 200次眼電去除實驗的平均MSE和SNR值。圖8展示了使用KICMME以及多種基于BSS或EMD方法(包括KICA-EEMD、CCA-EEMD[13]、FKD[9]、ICA-EMD[11]、KICA[9]和經(jīng)典ICA[6])去除眼電后3個導(dǎo)聯(lián)EEG信號的平均MSE和SNR值。

由圖8可見，經(jīng)典ICA算法的MSE值最高，SNR最低，這表明ICA難以準確識別并去除EEG信號中EOG偽跡。與ICA相比，KICA由于引入了核函數(shù)使得性能得到改善。不足的是，這兩種方法直接將與EOG相關(guān)的IC置零，這難免會丟失部分有用的腦電數(shù)據(jù)。使用ICA-EMD和KICA-EEMD時，通過EMD和EEMD算法進一步將EOG相關(guān)的IC分解為一組IMF，提取其中的EEG數(shù)據(jù)。然而，EMD和EEMD難以抑制由尖峰噪聲引起的模態(tài)分裂效應(yīng)，導(dǎo)致提取的腦電數(shù)據(jù)中含有部分眼電成分，這使得ICA-EMD和KICA-EEMD這兩種方法的偽跡去除效果也不是最優(yōu)的。相比之下，MAMA-EMD算法能夠?qū)⒓夥宸蛛x至第一個IMF中，緩解了模態(tài)分裂問題，有利于更精確地提取并保留腦電IMF，因此， KICMME的MSE值更低，SNR值更高。此外，圖8還給出了FKD和CCA-EEMD這兩種方法的實驗結(jié)果，對比可以發(fā)現(xiàn)KICMME方法具有最小的均方誤差0.82μV2和最高的信噪比12.51 dB，這表明KICMME相比于其他方法能夠在去除EOG偽跡的同時，保留更多有用的腦電數(shù)據(jù)。

圖8 多種眼電偽跡去除方法的平均MSE和SNR對比Figure 8 Comparison of average MSE and SNR for various EOG artifact removal methods

2.2 基于真實EEG信號的眼電去除實驗

2.2.1 數(shù)據(jù)集介紹

真實腦電數(shù)據(jù)來自于BCI Competition IV I公開數(shù)據(jù)集[16]，該數(shù)據(jù)集共記錄了7名健康受試者的運動想象腦電信號(motor imagery EEG, MI-EEG)。對于每名受試者，想象任務(wù)為左手運動、右手運動和腳部運動中的兩類，每類任務(wù)進行100次。采集的EEG信號進行0.05～49 Hz帶通濾波，并降采樣至100 Hz。選取與運動想象最相關(guān)的15個導(dǎo)聯(lián)(F3、F1、Fz、F2、F4、C5、C3、C1、Cz、C2、C4、C6、P1、Pz、P2)，截取0～4 s的運動想象期數(shù)據(jù)，進行EOG偽跡去除實驗。

2.2.2 EOG偽跡去除與結(jié)果展示

對于真實腦電信號的EOG偽跡去除過程類似于2.1.2節(jié),其中閾值st和kt分別設(shè)為15和0.8。圖9展示了原始受污染的多通道EEG信號，圖10對應(yīng)為使用KICMME方法恢復(fù)的無EOG偽跡的腦電數(shù)據(jù)，通過對比可以看出，KICMME方法對于真實記錄的EEG信號也具有較好的偽跡去除效果。

圖9 真實受污染的多通道腦電信號Figure 9 Real multichannel contaminated EEG signals

圖10 基于KICMME方法重構(gòu)的腦電信號Figure 10 The reconstructed EEG by using KICMME

2.2.3 基于分類性能的效果評估

由于真實記錄的EEG信號中隱含的純凈腦電數(shù)據(jù)是未知的，因而無法使用MSE和SNR等常用指標來衡量EOG偽跡去除方法的性能。就BCI系統(tǒng)而言，不同類別的EEG信號能否被正確識別直接決定了系統(tǒng)的穩(wěn)定性，且已有研究表明，EOG偽跡會降低BCI設(shè)備的分類性能[17]。因此，在去除EOG偽跡后，利用共空間模式(common spatial pattern, CSP)算法提取MI-EEG信號的空間特征，并采用支持向量機(support vector machine, SVM)分類器進行模式分類，基于分類性能評估KICMME方法應(yīng)用于真實EEG信號時的偽跡去除效果。

基于20次5折交叉驗證計算多種偽跡去除方法恢復(fù)的MI-EEG信號的平均分類準確率和平均Kappa值，結(jié)果如表2所示，表中數(shù)據(jù)為7名受試者的平均值。從表中可以看出，基于KICMME得到的平均分類準確率達到了91%，平均Kappa值為0.82，相比于其他基于BSS或EMD的偽跡去除方法具有明顯的提升。

表2 基于不同偽跡去除方法的平均分類準確率和平均Kappa值(n=7)Table 2 The average classification accuracies (Acc) and average Kappa values of different artifact removal methods(n=7)

3 討論與結(jié)論

EOG對EEG的尖峰狀污染是神經(jīng)信息分析所面臨的一個重要問題，它可能會干擾腦疾病的診斷結(jié)果，甚至誤導(dǎo)BCI系統(tǒng)的實際應(yīng)用。本文提出一種能夠自動識別和去除多通道EEG信號中EOG偽跡的新方法，即KICMME。該方法首先利用KICA分離出EOG相關(guān)的ICs，然后基于MAMA-EMD對其進行進一步分解，提取有用的腦電數(shù)據(jù)。在半模擬和真實腦電信號上的實驗結(jié)果表明，MAMA-EMD能夠有效消除EMD分解過程中由尖峰引起的模態(tài)分裂效應(yīng)，使得與EOG相關(guān)的IMF的識別與去除更加準確；同時，本文提出的KICMME方法在MSE、SNR和分類性能上均優(yōu)于其他基于BSS的EOG偽跡去除方法，如CCA-EEMD、FKD、ICA-EMD、KICA和經(jīng)典ICA等。這一改進為EEG信號的預(yù)處理提供了一個新的思路，并將擴大其在腦科學(xué)研究和BCI系統(tǒng)中的應(yīng)用。

本文的局限性在于KICMME首先對EEG信號進行KICA分離，因此無法處理單通道受污染的EEG信號。在未來的研究中，課題組將考慮將MAMA-EMD算法應(yīng)用于原始信號，然后結(jié)合BSS算法去除單通道EEG信號中的EOG偽跡。