基于LMD和CSP的多域融合腦電信號分類方法

陳 舒 周 青

1(廣州市公用事業(yè)技師學(xué)院信息傳媒產(chǎn)業(yè)系 廣東 廣州 510100) 2(中山大學(xué)信息科技學(xué)院 廣東 廣州 510006)

0 引 言

腦-機接口(Brain Computer Interface, BCI)是一種將人腦與外部計算機等輔助設(shè)備直接聯(lián)系起來的實時通信系統(tǒng),BCI技術(shù)直接對腦電信號進(jìn)行采集和分析,提取其中包含的運動意圖信息并進(jìn)一步將其解譯為驅(qū)動外部設(shè)備的命令,能夠幫助運動障礙患者恢復(fù)與外界的信息交互,提升其生活便利性[1-3]。腦電圖(Electro-Encephalon-Graph, EEG)由于其低成本和非侵入等優(yōu)點成為當(dāng)前BCI系統(tǒng)中應(yīng)用最為廣泛的一種腦電信號。

EEG信號的產(chǎn)生機理復(fù)雜,是一種典型的非平穩(wěn)、非線性和微弱的時變信號。如何從中提取反映患者不同思維活動的差異信息并將其與特定的控制命令聯(lián)系起來是BCI技術(shù)的難點,也是當(dāng)前研究的熱點[4]。目前常用的方法可以分為時-頻變換類方法和空間濾波類方法,其中典型的時-頻變換類方法有短時傅里葉變換(Short Time Fourier Transformation, STFT)、小波變換ST和經(jīng)驗?zāi)B(tài)分解(Empirical Mode Decomposition, EMD)等。其中,FT將時域EEG信號轉(zhuǎn)換到時-頻域,能夠較好地刻畫信號的局部特征,但是受制于海森堡測不準(zhǔn)原理,該方法無法同時獲得較高的時間分辨率和頻率分辨率[5];小波變換作為一種經(jīng)典的非平穩(wěn)、非線性信號分析方法,能夠?qū)π盘栠M(jìn)行不同“尺度”的細(xì)化分析,被廣泛應(yīng)用于EEG信號分類中,然而“小波基”函數(shù)的選擇和分解層數(shù)的確定對小波方法的性能影響較大[6-8];EMD不需要預(yù)先設(shè)定“基函數(shù)”,直接從原始信號中提取本征模函數(shù),能夠同時獲得高時間分辨率和頻率分辨率,但是邊緣效應(yīng)的存在制約了該類方法的應(yīng)用[9]。與時-頻變換類方法不同的是,空間濾波類方法從提升EEG信號空間分辨率的角度出發(fā),通過多通道的方式對不同腦域的特征信息進(jìn)行提取和分析,典型的空間濾波類方法有獨立主分量分析法,共空間模式法(Common Spatial Pattern, CSP)、空間稀疏共空間模式(Spatially Sparse Common Spatial Pattern, SSCSP)、濾波器組共空間模式(Filter Bank Common Spatial Pattern, FBCSP)等多種CSP變化方法[10-13]。該類方法相對于時-頻變換類方法由于采用了更多觀測通道,通?？梢垣@得較高的分類結(jié)果,但是多通道的數(shù)據(jù)獲取是以增加系統(tǒng)復(fù)雜度為代價的,而且該類方法丟失了頻率信息。

上述分析表明,現(xiàn)有方法雖然能夠提取出反映不同運動腦電信號中包含的差異信息,但是由于這些方法僅僅是從某個特定的變換域?qū)EG進(jìn)行分析,提取的特征較為單一,沒能包含不同維度的信息,具有一定的局限性,同時在低信噪比條件下會出現(xiàn)分類識別性能下降的問題[14-15]。

針對上述問題,本文提出一種基于局部均值分解(LMD)和CSP的多域融合腦電信號分類方法,首先利用LMD對EEG信號進(jìn)行分析,得到各階乘積分量(PF)的同時提取時-頻域特征,然后對各階PF分量進(jìn)一步作CSP分解,獲取EEG信號的空間分布特征,從而將時-頻域特征擴展為時-頻-空域特征。最后采用相關(guān)向量機(RVM)[16]對特征選擇和分類函數(shù)進(jìn)行聯(lián)合優(yōu)化,剔除冗余特征并得到最優(yōu)的特征組合?；贐CI競賽數(shù)據(jù)集的實驗結(jié)果表明,所提方法可以實現(xiàn)對EEG信號的多角度描述,能夠獲得更優(yōu)的分類性能。

1 實驗數(shù)據(jù)

本文實驗選用與文獻(xiàn)[10]一致BCI2003 Competition Ⅱ中DatasetⅢ標(biāo)準(zhǔn)數(shù)據(jù)集。實驗中對一名25歲健康女性左右手運動想象(MI)EEG信號進(jìn)行采集和記錄。EEG信號采用國際標(biāo)準(zhǔn)的10～20導(dǎo)聯(lián)系統(tǒng)的C3和C4通道,采用差分電極。采樣頻率為128 Hz,并對頻率為0.5～30 Hz范圍內(nèi)的信號進(jìn)行濾波。實驗中,受試者按要求坐在屏幕前,根據(jù)屏幕上出現(xiàn)的光標(biāo)提示進(jìn)行想象左右手運動,單次實驗流程如圖1所示。具體說明如下:

1) 0 s≤T<2 s,受試者前方屏幕顯示空白,提示受試者放松。

2)T=2 s,受試者前方屏幕出現(xiàn)“+”光標(biāo),提示受試者做好準(zhǔn)備。

3) 3 s≤T<9 s,受試者前方屏幕隨機出現(xiàn)“←”或“→”,受試者根據(jù)提示想象左手或者右手運動。

4)T=9 s,受試者前方屏幕再次顯示空白,提示受試者單次實驗結(jié)束。

圖1 單次實驗數(shù)據(jù)采集時序

整個數(shù)據(jù)集中包含7組,每組40次實驗,總共280組實驗數(shù)據(jù),包含想象左右手運動數(shù)據(jù)各140組。由于每次實驗從第3 s正式開始,因此本文取每組實驗數(shù)據(jù)中第3～9 s對應(yīng)的768個采樣點作為有效數(shù)據(jù)。

2 時-頻-空域特征提取

2.1 基于LMD的時-頻域特征提取

EEG信號的幅度和頻率均隨著時間的變化而變化,是一種典型的調(diào)頻調(diào)幅信號,其瞬時頻率為信號中調(diào)頻成分的瞬時頻率,幅度為信號中調(diào)幅成分的瞬時幅度,因此如果能將EEG信號中的調(diào)頻成分和調(diào)幅成分進(jìn)行分離,就能進(jìn)一步對信號瞬時頻率和瞬時幅度的細(xì)節(jié)信息進(jìn)行分析,從而可以提取更多對EEG信號進(jìn)行分類的有用信息。

LMD是在EMD基礎(chǔ)上發(fā)展起來的一種新的非平穩(wěn)復(fù)雜信號自適應(yīng)分解算法,LMD可以將信號分解為若干個具有物理意義PF分量之和的形式,而獲得PF分量的本質(zhì)是從原始信號中分離出純調(diào)頻信號和包絡(luò)信號,因此LMD非常適合于分析EEG信號。

對任意給定EEG信號,對其進(jìn)行LMD分解的步驟如下:

(1) 計算均值函數(shù)m11(t)。

Step1對輸入信號s(t)進(jìn)行分析,計算其中包含的全部極值點ni,i=1,2,…,W。

Step2根據(jù)式(1)計算s(t)的均值序列mi:

Step3利用移動平均法對mi構(gòu)成的曲線進(jìn)行平滑,得到輸入信號s(t)的均值函數(shù)m11(t)。

(2) 計算包絡(luò)估計函數(shù)a11(t)。

Step1根據(jù)式(2)計算包絡(luò)估計值ai:

Step2利用移動平均法對ai構(gòu)成的曲線進(jìn)行平滑,得到包絡(luò)估計函數(shù)a11(t)。

(3) 計算差值信號h11(t)。

從原始信號中減去均值函數(shù)m11(t),得到差值信號h11(t)=s(t)-m11(t)。

(4) 計算純調(diào)頻信號s1j(t)。

Step1對差值信號h11(t)按式(3)進(jìn)行幅度解調(diào)處理,得到s11(t):

Step2將s11(t)作為步驟(1) 的Step 1中的輸入信號s(t),重復(fù)進(jìn)行步驟(1)-步驟(3),直到第j次迭代后,步驟(2)得到的包絡(luò)函數(shù)a1j(t)=1,此時得到的s1j(t)即為純調(diào)頻信號。

(5) 計算幅值函數(shù)a1(t)。

(6) 計算第一個PF分量f1(t)。

將a1(t)與s1j(t)相乘,得到PF分量f1(t)=a1(t)s1j(t)。

(7) 計算剩余信號u1(t)。

從原始信號s(t)中減去f1(t)得到殘留信號u1(t)=s(t)-f1(t)。

(8) 計算所有L個PF分量。

令s(t)=u1(t),重復(fù)進(jìn)行步驟(1)-步驟(7),直到獲得L個PF分量,且剩余信號U(t)為單調(diào)函數(shù)時,迭代終止,此時原始信號可以表示為:

(4)

圖2(a)和圖2(b)分別給出了利用上述LMD方法對想象左手運動和想象右手運動腦電信號進(jìn)行分析得到的分解結(jié)果。其中第1行子圖為原始信號,第2-第5行子圖為分解得到的PF1-PF4,最后一行子圖為剩余信號。從圖2可以看出,隨著分解層數(shù)的增加,PF信號的頻域特性逐漸減少,時域特性逐漸增加,因此PF分量中包含了信號的時域和頻域特性。同時對比圖2(a)和圖2(b)可以看出,兩種腦電信號分解結(jié)果呈現(xiàn)出了一定的差異性,因此我們可以提取以下時-頻域特征對LMD分解結(jié)果的差異性進(jìn)行描述。

(a) 想象左手運動

(b) 想象右手運動圖2 LMD分解結(jié)果

特征1:波形寬度特征:

式中:find(·)為計算滿足括號內(nèi)要求的樣本點個數(shù),i=1,2,…,4。波形寬度特征能夠體現(xiàn)信號頻率分布差異。

特征2:波形熵特征:

(6)

特征3:波形四階距特征:

(7)

對PF1-PF4分量分別提取上述三維特征,構(gòu)成12維特征向量如表1所示。

表1 LMD分解得到的時-頻域特征向量

2.2 基于CSP的空域特征提取

CSP是一種針對兩類任務(wù)的經(jīng)典空域濾波方法,能夠從多通道EEG數(shù)據(jù)中提取每一類信號的空間分布模式,其基本原理是通過矩陣對角算法構(gòu)建空間濾波器,將多通道EEG數(shù)據(jù)投影到低維子空間并使兩類樣本的協(xié)方差矩陣差異最大化。

對于實驗中獲得的N通道想象左、右手運動腦電信號數(shù)據(jù)矩陣分別為EL∈RN×T和ER∈RN×T,T為每個通道的采樣點數(shù)。對其進(jìn)行CSP分解的步驟如下:

(1) 計算歸一化協(xié)方差矩陣。根據(jù)式(8)可以計算得到想象左、右手運動腦電信號數(shù)據(jù)矩陣對應(yīng)的歸一化協(xié)方差矩陣RL和RR。

式中:trace(·)表示矩陣求跡運算。

(2) 構(gòu)建復(fù)合協(xié)方差矩陣并進(jìn)行特征值分解。

Step1對RL和RR求和得到復(fù)合協(xié)方差矩陣R=RL+RR。

Step2對R進(jìn)行特征值分解得到特征值矩陣Σ。

R=UΣU

(9)

式中:U為特征向量矩陣;Σ為特征值構(gòu)成的對角矩陣。

(3) 對RL和RR進(jìn)行白化處理。

Step1根據(jù)式(10)計算白化矩陣P。

(10)

Step2根據(jù)式(11)利用白化矩陣P對RL和RR進(jìn)行白化處理。

(4) 對SL和SR進(jìn)行特征分解。

根據(jù)式(12)對SL和SR進(jìn)行特征分解:

(12)

根據(jù)式(13)利用空間濾波器W對多通道信號矩陣E進(jìn)行空域濾波得到特征矩陣:

Z=WE

(13)

計算特征矩陣前m行和后m行對應(yīng)的方差構(gòu)成空域特征向量:

特征4:空域特征。

(14)

式中:var(·)為求方差運算。本文取m=9,即總共有18維特征構(gòu)成空域特征向量。

3 基于RVM的特征選擇和分類

得到特征向量后,需要設(shè)計合適的分類算法才能獲得期望的分類性能,目前常用的分類算法有線性分類器、支撐向量機分類器和貝葉斯分類器等。本文第2節(jié)提取的30維時-頻-空域特征雖然實現(xiàn)了特征域的擴展,但是高維特征向量中難免會存在一些冗余特征,這些特征的存在會增加分類算法的復(fù)雜度,降低實時性,因此本文采用RVM進(jìn)行特征選擇和分類識別的聯(lián)合優(yōu)化,降低算法復(fù)雜度的同時,提升分類性能。

RVM理論是在SVM基礎(chǔ)上發(fā)展起來的一種基于貝葉斯框架的監(jiān)督類學(xué)習(xí)方法,其基本原理也是利用核函數(shù)將低維空間中的線性不可分問題映射到高維空間中的線性可分問題,但不同于SVM的是,RVM具有更高的稀疏性,核函數(shù)的選擇不受摩西準(zhǔn)則[16]限制,并且能夠提供概率式的預(yù)測。

對于想象左手運動和想象右手運動腦電信號分類問題,完整的RVM模型可以表示為:

式中:xn(n=1,2,…,N,N為樣本總數(shù))為第n個訓(xùn)練樣本的特征向量;K(x,xn)為核函數(shù);w=[w1,w2,…,wn]T為服從0均值、協(xié)方差矩陣為α-1I高斯分布的權(quán)重向量;ε為服從0均值、協(xié)方差矩陣為γ-1I高斯分布的噪聲分量,與信號分量相互獨立。為了構(gòu)建完整的貝葉斯模型,RVM進(jìn)一步假設(shè)α和γ服從超參數(shù)為a0、b0、c0和d0的伽馬分布。

圖3利用RVM算法對上述時-頻-空域33維特征向量進(jìn)行特征選擇迭代終止時的歸一化權(quán)值,可以看出經(jīng)過特征選擇后,權(quán)值較大的有6維特征,對應(yīng)4維時-頻域特征Feature11、Feature12、Feature24和Feature33,以及2維空域特征Feature47和Feature416,體現(xiàn)了本文方法信息提取維度的多元化。

圖3 RVM特征選擇結(jié)果

4 實驗及結(jié)果分析

圖4給出了本文所提基于LMD和CSP的多域融合腦電信號分類方法流程,可以看出算法包含訓(xùn)練和測試兩個階段。

圖4 本文分類算法流程

在訓(xùn)練階段,首先對訓(xùn)練樣本進(jìn)行LMD分解,得到PF分量并提取12維時-頻域特征,然后將PF分量作為CSP的多通道數(shù)據(jù)進(jìn)行空域濾波并提取18維空域特征,然后利用RVM分類器對30維時-頻-空域特征進(jìn)行特征選擇和分類。

在測試階段,對測試樣本進(jìn)行LMD分解,然后提取訓(xùn)練階段RVM分類器選出的時-頻域特征,同樣將PF分量作為CSP的多通道數(shù)據(jù)進(jìn)行空域濾波,然后提取訓(xùn)練階段RVM分類器選出的空域特征,并將其與時-頻域特征一起構(gòu)成特征向量送入RVM分類器進(jìn)行分類。

4.1 單一維度特征和時-頻-空三維特征對比

為了驗證時-頻-空域融合特征的分類性能,將數(shù)據(jù)集隨機分成5個子集,采用5組交叉驗證法進(jìn)行實驗,即依次選取1組子集作為訓(xùn)練樣本,剩余4組子集作為測試樣本,記錄每次實驗的分類結(jié)果如圖5所示,表2給出了5次實驗的統(tǒng)計結(jié)果。為了進(jìn)行對比,圖5和表2給出了在相同條件下分別單獨采用時-頻域特征、空域特征進(jìn)行分類的結(jié)果。從圖5和表2的結(jié)果可以看出,所提時-頻-空域特征相對于單一的時-頻域特征和空域特征分類性能提升了5%,并且對不同數(shù)據(jù)集的分類穩(wěn)定性更好。同時只采用RVM分類器選擇的6維時-頻-空域特征與直接采用30維特征的分類性能接近、方差更小,表明RVM選擇的6維特征對分類性能的貢獻(xiàn)率要遠(yuǎn)超其他24維特征。

圖5 不同特征單次分類結(jié)果

表2 5次實驗平均分類結(jié)果

4.2 不同算法分類性能對比

表3給出了本文方法與近年來其他文獻(xiàn)中介紹方法在相同條件下的分類結(jié)果和F1得分指標(biāo),其中文獻(xiàn)[17]采用的是小波包能量特征以及線性分類器,文獻(xiàn)[18]采用的是小波熵特征以及SVM分類器;文獻(xiàn)[19]采用的是5階AR模型系數(shù)以及線性分類器。從表3可以看出,本文方法可以獲得最優(yōu)的分類性能。

表3 不同方法的分類性能

除了分類性能外,算法執(zhí)行時間也是評價算法優(yōu)劣的一項重要指標(biāo)。表3給出了不同方法在MATLABr2015平臺的算法執(zhí)行時間,可以看出本文算法的實時性明顯優(yōu)于文獻(xiàn)[18]和文獻(xiàn)[19]所述方法,僅比文獻(xiàn)[17]采樣線性分類器慢0.4 s,表明本文方法兼具高分類性能和高算法實時性的優(yōu)點。

4.3 不同算法低信噪比條件下魯棒性對比

微弱性是腦電信號的一個顯著特點,其電壓通常為μV量級,容易受到噪聲污染。因此噪聲穩(wěn)健性是評估分類算法優(yōu)劣的另一項關(guān)鍵指標(biāo)。為了驗證本文方法在不同信噪比條件下的分類性能,采取向?qū)嶒灁?shù)據(jù)中加入高斯白噪聲的方式進(jìn)行驗證。

圖6給出了上述方法在不同信噪比條件下的分類性能?？梢钥闯?本文方法在低信噪比條件下的表現(xiàn)要明顯優(yōu)于其他對比方法,當(dāng)信噪比高于15 dB時,本文方法可以獲得優(yōu)于90%的分類正確率。這是由于本文方法從時-頻-空多個維度進(jìn)行特征提取,降低了噪聲分量的影響,同時RVM分類器具有的噪聲魯棒性也提升了低信噪比條件下的魯棒性。

圖6 分類性能隨信噪比的變化曲線

5 結(jié) 語

針對運動腦電信號非平穩(wěn)、時變和微弱性的特點,本文提出一種基于LMD和CSP的多域特征融合分類方法,綜合考慮了EEG信號的時-頻-空域信息。算法首先利用LMD對EEG信號進(jìn)行時-頻域特征提取,并將LMD分解得到的PF分量作為CSP的多通道數(shù)據(jù)進(jìn)行空域特征提取,從而得到30維時-頻-空域特征向量。然后利用RVM對高維特征向量進(jìn)行自適應(yīng)特征選擇與分類,從30維特征中自動選取對分類性能貢獻(xiàn)最大的6維特征,在降低算法復(fù)雜度的同時保證了分類性能。基于實際數(shù)據(jù)的實驗結(jié)果表明,相對于單一維度時-頻特征和空域特征,本文方法識別率更高、標(biāo)準(zhǔn)差更小、魯棒性更好,與現(xiàn)有方法的對比結(jié)果表明,本文方法可以獲得更高的分類性能和噪聲穩(wěn)健性,具有較好的應(yīng)用前景。