EMD融合WPT和CSP的腦電特征提取方法

張學(xué)軍，劉定宇，霍 延

(1.南京郵電大學(xué) 電子與光學(xué)工程學(xué)院、微電子學(xué)院，江蘇 南京 210023； 2.南京郵電大學(xué) 射頻集成與微組裝技術(shù)國家地方聯(lián)合工程實驗室，江蘇 南京 210023)

0 引 言

腦機接口(brain computer interface，BCI)是涉及眾多學(xué)科和知識領(lǐng)域的控制技術(shù)。腦機接口的概念于二十世紀(jì)末期被首次提出，它建立了人腦與外部設(shè)備的通道，通過采集腦電信號，使用分類算法來識別人腦的指令[1]。腦機接口技術(shù)在近幾十年取得了極大發(fā)展，應(yīng)用領(lǐng)域也越發(fā)廣泛，尤其在運動障礙患者的康復(fù)訓(xùn)練中得到了大量的應(yīng)用，病人通過該技術(shù)提高了與外界的交互能力和自理能力。通過大腦信號，可以直接控制相關(guān)的外部設(shè)備，如計算機、機械小車，機械假肢[2]等，這些應(yīng)用都增強了他們與外界的交流并加快了他們的恢復(fù)。腦機接口技術(shù)還被廣泛地應(yīng)用于其他一些領(lǐng)域，如游戲應(yīng)用和導(dǎo)航[1]。

腦機接口技術(shù)包含信號采集、預(yù)處理、特征提取、分類以及外部設(shè)備控制。特征提取用于識別不同腦電信號的特征信息，因此在腦機接口的研究領(lǐng)域備受關(guān)注。之前研究的比較成熟的方法如短時的傅里葉變換(short-time Fourier transform，STFT)、小波變換(wavelet transform，WT)等等，由于EEG信號的非線性以及非平穩(wěn)性，這些傳統(tǒng)方法不能得到理想的分辨效果。經(jīng)過近年來不斷的研究，希爾伯特黃變換(Hilbert-Huang transform，HHT)得到了發(fā)展，經(jīng)過驗證能夠更好地對EEG信號進行分析。采集的腦電信號，首先經(jīng)過EMD[3]分解后得到一串的IMFs，隨后通過HHT可以得到相應(yīng)的能量譜數(shù)據(jù)，得到的數(shù)據(jù)可以作為一個特征指標(biāo)，用于后續(xù)的分類。HHT被提出后，其優(yōu)秀的性能，集中的特征和較高的分辨率，使其應(yīng)用越發(fā)廣泛[4]。對于EEG信號，目前方法得到的分類效果都無法達(dá)到期望，如需應(yīng)用在BCI系統(tǒng)中，需對其進行更復(fù)雜的處理，可以先進行空間范疇的濾波，獲得高質(zhì)量的信號，提取更優(yōu)質(zhì)的特征信息。常用的算法包括公共空間模式CSP[5]，以及獨立分量分析、濾波器組公共空間模式、判別濾波公共空間模式等改進的CSP算法[6]。傳統(tǒng)的經(jīng)驗?zāi)Ｊ椒纸庠诘皖l區(qū)域首次得到的固有模態(tài)函數(shù)涵蓋過寬的頻率范圍，無法得到單組份的屬性，且不能分開含有低能量成分的信號。為了彌補此類問題，該文提出了一種改進方案，即在傳統(tǒng)的EMD分解之前，使用小波包變換對信號進行處理將原始的信號分為一組窄帶信號，并利用篩選去除不相關(guān)的固有模態(tài)函數(shù)。同時常規(guī)的公共空間模式分解方法在通道數(shù)及頻域方面都有明顯的缺陷[7]。為了彌補這些短板，該文將EMD和CSP算法相結(jié)合，篩選IMF頻段，剩下的計算其能量的特征，使用CSP濾波器，得到其特征的一組數(shù)據(jù)，最后使用SVM進行分類，可以達(dá)到95.9%的平均分類準(zhǔn)確率。

1 相關(guān)算法

1.1 經(jīng)驗?zāi)Ｊ椒纸?/h3>

EMD因為在非線性序列處理方面的優(yōu)越性而被廣泛用于處理非線性數(shù)據(jù)。在時間域中研究其本征波動模式[8]并進行分解是其核心的思想。這種分解過程，可以將其理解為篩選的過程。需要有以下條件：(1)需得到數(shù)據(jù)的極值；(2)需知兩個極值之間的時間間隔定義的特征尺度；(3)如果未滿足上述條件，但是函數(shù)存在拐點并可以通過微分的方法經(jīng)過一次或多次求得極值也可以被認(rèn)為滿足條件，但是需要注意分解結(jié)果也相應(yīng)需要通過積分獲得。

EMD分解具體步驟：

步驟一：讀取信號x(t)的所有極值信息，然后用三次樣條函數(shù)對它們進行處理，得到序列的包絡(luò)線信息，上包絡(luò)vmax(t)，下包絡(luò)vmin(t)，以此進一步得到包絡(luò)線平均值m(t)。

(1)

步驟二：計算原信號與均值的差值c(t)：

c(t)=x(t)-m(t)

(2)

IMF定義的截止條件：在時間域內(nèi)，極值個數(shù)等于過零點的個數(shù)，或最大差值為1，且最大信號幅度形成的包絡(luò)線和最小信號幅度形成的包絡(luò)線算術(shù)平均為零，如果不滿足，就重復(fù)上述過程，剩余量r(t)的公式為：

r(t)=x(t)-c(t)

(3)

步驟三：作為一種新數(shù)據(jù)，剩余量通過同一過程過濾，得到下一個低頻固有模函數(shù)。直到剩余函數(shù)r(t)得到的是單調(diào)函數(shù)或只有一個極值分解的過程才會終止。假設(shè)將原始信號x(t)進行分解，得到了n個本征模函數(shù)和一個剩余函數(shù)r(t)，重構(gòu)信號如下：

(4)

1.2 公共空間模式

目前基于運動圖像的BCI系統(tǒng)中有功率譜密度[9]等方法。CSP就是一種常見的處理方法。其對輸入的腦電數(shù)據(jù)進行映射處理[10]，通過這個過程獲得的特征向量可以得到較好的分辨度。它以兩個不同類別的協(xié)方差矩陣作為基礎(chǔ)，通過對角化處理設(shè)計最優(yōu)化的空間濾波器。對原始EEG數(shù)據(jù)進行映射，獲得最大程度上的分離[11]。具體步驟為：

步驟一：設(shè)定單次想象運動的腦電信號矩陣為X(N*I維)，其中N表示EEG采集的通道數(shù)而I表示各通道的采樣點數(shù)．這里將X進行歸一化處理以此得到相關(guān)的協(xié)方差矩陣：

(5)

步驟二：對合成的協(xié)方差矩陣進行分解。

(6)

其中，Σ為特征值對角矩陣，U0為其對應(yīng)的特征向量矩陣，其隨相應(yīng)的特征值而重新改變。

步驟三：求白化矩陣，白化矩陣P定義如下：

(7)

步驟四：白化協(xié)方差矩陣R1和R2。

(8)

步驟五：主成分分解。

(9)

步驟六：求CSP的空間濾波器W。

W=UTP

(10)

運動圖像腦電信號x，經(jīng)空間濾波z=wx，得到z，特征向量數(shù)為2m。為使兩個樣本得到最大距離，前m個特征向量，取它們與一類運動想象方差的最大值，后m個則取最小值。而后面剩下的m維則使用相反方法[12]。其中第1個和最后1個特征向量，有最大的區(qū)別度，第2個和倒數(shù)第2個相比較小，以此類推。所以應(yīng)選取合適的向量個數(shù)以使其包含最佳的特征信息。

步驟七：特征提取后，再將空間投影后的信號Zp(p=1,…,2m)取對數(shù)，從而使特征值差異更明顯，其中var表示求向量的方差。

(11)

1.3 支持向量機

支持向量機(SVM)[13]，為了得到更好的分離特征向量，將低維信號通過映射的方式放到高維度的特征空間，用超平面的分類方式對數(shù)據(jù)進行下一步的處理。超平面的概念可以理解為，高階的數(shù)據(jù)須由低一維的數(shù)據(jù)進行分割處理，比如，用一條直線劃分一個二維的平面區(qū)域，用一個二維的平面來分割一個三維的立體。所以如何選取一個最好的超平面是其主要的關(guān)注點，做到分類間隔盡量大，這樣準(zhǔn)確率就可以得到提升。同時，核參數(shù)和誤差懲罰因子對向量機的分類也有很大的影響。

2 EMD結(jié)合WPT和CSP算法

2.1 經(jīng)驗?zāi)Ｊ椒纸飧倪M

傳統(tǒng)的EMD在低頻區(qū)域會產(chǎn)生不理想的固有模態(tài)函數(shù)，且首次得到的固有模態(tài)函數(shù)涵蓋太寬的頻率范圍，無法得到單組份的屬性，同時也不能分開含有低能量成分的信號。該文提出了一種改進方法來彌補這些不足，在傳統(tǒng)的EMD之前，先使用小波包變換對信號進行處理，將原始的信號分為一組窄帶的信號，并且利用一個篩選過程從結(jié)果中去除不相關(guān)的固有模態(tài)函數(shù)。

2.1.1 小波包變換過程

為了解決上述問題，這里先使用了小波包變換。

對于一個n層的分解，信號被分為2n個窄帶信號。由于分解是基于頻率而不是基于能量，因此，低能量成分被分解為不同的頻帶保留下來。所以，小波包變換首先用來得到若干窄帶信號，然后對這些窄帶信號進行分解，于是得到的固有模態(tài)函數(shù)同樣有窄頻帶而且它們的瞬時頻率更加接近IMF的真實頻率模式。

由于腦電信號強度較弱，且采集的時候有如眼、肌電等眾多的影響源因素，這些相對的噪聲都會影響腦電信號的質(zhì)量。所以該文使用了小波包變換，以此在頻域?qū)ζ溥M行過濾，有效地提高了腦電信號的信噪比。小波包變換源于小波變換，所以不僅繼承了其優(yōu)點，更是在基礎(chǔ)上有部分改進。它對信號進行正交分解，同時對于頻段的選擇更加準(zhǔn)確，能更好地匹配信號的有用頻段，獲得更有效的頻譜信息，因此其分辨率也得到了提高。小波包變換的這種多分辨率特性可以在所在的頻率范圍內(nèi)鎖定更加有用的信息。在小波多分辨率信號分析中，希爾伯特空間可以變換不一樣的尺度指標(biāo)分解得到所有的小波子空間的正交和，并根據(jù)二進制系統(tǒng)進行細(xì)分，如圖1所示。

圖1 小波包的空間分解

2.1.2 篩選過程

此處將信號和IMF的歸一化相關(guān)系數(shù)μ作為標(biāo)準(zhǔn)來對固有模態(tài)函數(shù)進行篩選。

獲得所有的相關(guān)系數(shù)[14]μi(i=1∶n，n是IMF的數(shù)量)，然后通過下式計算硬閾值：

(12)

其中，η(>1)是一個比例系數(shù)，經(jīng)驗值為10。每一個固有模態(tài)函數(shù)都和λ進行比較，如果一個固有模態(tài)函數(shù)的相關(guān)系數(shù)大于λ，相應(yīng)的固有模態(tài)函數(shù)就保留；如果不符合上述要求，相應(yīng)的固有模態(tài)函數(shù)就被刪除。

該文考慮了兩種篩選過程，第一種用來排除由WPT得到的窄帶信號關(guān)聯(lián)性較差的IMF；第二種用來去除由窄帶信號產(chǎn)生的弱相關(guān)的IMF，這些窄帶信號本身就與原始信號相關(guān)性很差。

2.2 改進EMD結(jié)合WPT和CSP算法流程

該文提出的改進方案，通過采集獲得腦電信號后，經(jīng)由濾波器先進行預(yù)處理。接著在傳統(tǒng)的EMD分解之前，使用小波包變換對信號進行處理，將原始的信號分為一組窄帶信號，篩選去除結(jié)果中相關(guān)較差的固有模態(tài)函數(shù)。同時再將EMD和CSP算法結(jié)合，篩選IMF頻段，剩下的計算其能量的特征，使用CSP濾波器，得到其特征的一組數(shù)據(jù)，最后使用SVM進行分類。算法流程如圖2所示。

圖2 算法流程

3 數(shù)據(jù)集描述

為了便于與經(jīng)典方式進行比較實驗，采用BCI2008年的競賽數(shù)據(jù)。實驗設(shè)備包括電腦及電極帽等腦電采集設(shè)備，受試者頭戴電極帽坐于屏目前約1米處，根據(jù)提示做相關(guān)運動想象同時采集相關(guān)腦電信號。數(shù)據(jù)包含9名右利手的受試者的腦電數(shù)據(jù)。實驗開始時，受試者保持放松并坐于距離屏幕1米處。實驗中，屏幕上會顯示不同方向的箭頭作為想象的提示，受試者根據(jù)提示進行想象運動。在兩天里完成兩個session，每次session6個run，每個run20次實驗，左右手分別10次。所以在一次session中，左右手各60次。

實驗界面和時間序列如圖3所示。

圖3 實驗界面和時間序列

每一次實驗的時間，大概8至9秒，前兩秒放松；第二秒的時候會出現(xiàn)提示音指示實驗開始，屏幕上這個時候會出現(xiàn)一個十字型的光標(biāo)提示預(yù)示著開始進入準(zhǔn)備狀態(tài)，時長1秒；第三秒的時候屏幕上會隨機出現(xiàn)向左向右的箭頭并持續(xù)1.25秒；第四秒開始根據(jù)屏幕的箭頭提示，受試者進行3秒的想象運動并在第7秒時停止。使用Ag/Agcl電極,記錄C3，Cz，C4三個通道的數(shù)據(jù)。C3和C4兩個電極采用的腦電數(shù)據(jù)對應(yīng)運動想象最活躍的腦部部位，同時將Cz作為參考的信號。選用0.5 Hz～100 Hz的帶通濾波器。

4 實驗分析

4.1 小波包分解

根據(jù)事件相關(guān)同步/去同步原理，進行想象運動會在相關(guān)的頻率段產(chǎn)生能量的變動，因此，相關(guān)頻段μ節(jié)律，β節(jié)律的腦電就非常重要。信號頻帶寬100 Hz，這里將腦電信號進行4層分解，分解后得到16個不同的頻率段。其中最小的分辨率如式(13)所示。

Δf=100/16=6.25 Hz

(13)

頻段分解后，子帶信號所對應(yīng)的頻率范圍如表1所示。

表1 各子帶所對應(yīng)的頻率范圍 Hz

當(dāng)進行左右手想象運動的時候，μ節(jié)律變化主要范圍為7 Hz～13 Hz，β節(jié)律在19 Hz～26 Hz。對照表1，μ節(jié)律主要落在子帶(4，1)中，而子帶(4，3)可以看到β節(jié)律的頻段，所以，重構(gòu)(4，1)及(4，3)子帶，提取到μ節(jié)律和β節(jié)律腦電，保留相對有效的腦電信息，可以更好地提取特征。圖4分別展示了兩個子帶重構(gòu)后的波形及頻譜狀態(tài)，可以看出重構(gòu)后的子帶頻譜顯示出更加集中且明顯的特征。

圖4 小波分解子帶信號(上：波形，下：頻譜)

4.2 EMD分解結(jié)果

對腦電信號進行單純的經(jīng)驗?zāi)Ｊ椒纸鈺a(chǎn)生頻段涵蓋寬泛的問題，同時在低頻區(qū)域會產(chǎn)生不理想的IMFs，不能分開含有低能量成分的信號，從而會使結(jié)果產(chǎn)生錯誤。為了彌補這些不足，在EMD之前進行小波變換將信號分為一組窄帶信號，選取合適的子帶信號，分解為頻率更加集中的固有模態(tài)函數(shù)。

利用EMD對WPT處理后的子帶信號進行處理。前兩階IMF的相關(guān)系數(shù)遠(yuǎn)大于剩余的IMF，可見前兩階IMF有EEG信號最主要的信息，將最活躍的頻段包含在內(nèi)，而剩下的IMF則多包含噪聲及無用的信號，因此被舍棄。所以該文選擇僅重構(gòu)前兩階IMF。這里EMD的目的就可以看作是提取μ節(jié)律和β節(jié)律。同時，該方法可以獲取更多有效的信號，提高整個腦電信號的信噪比。圖5和圖6分別展示了(4，1)、(4，3)子帶信號經(jīng)過EMD分解后前兩階的IMF信息。

通過圖5和圖6可以看出，得到的兩個子帶具有比原始信號更大的有效范圍，包含了最活躍的頻段信息，且信噪比也更加優(yōu)越。

(a)時域圖

(b)頻域圖

(a)時域圖

(b)頻域圖

4.3 CSP濾波

在使用EMD之前，該文提取相關(guān)EEG信號，選取C3和C4通道的數(shù)據(jù)，并通過WPT得到一系列窄帶信號。然后對相應(yīng)的子帶信號進行EMD分解得到多個IMFs。選取前兩階作為新的輸入，形成矩陣XI(X=L左手，I=R右手)。選取兩子帶信號(4，1)、(4，3)，經(jīng)EMD分解并選前兩階IMFs。

4.4 實驗結(jié)果分析

將每次試驗得到的數(shù)據(jù)，先經(jīng)WPT分解，選取子帶信號，然后進行EMD分解，選擇IMF，組合并重構(gòu)前兩個IMF分量，接著利用CSP進行特征提取。圖7為單個實驗中的前兩個IMF分量經(jīng)過濾波后的腦地形圖。

圖7 想象運動前兩階IMF的腦地形圖

可以看到，C3和C4電極在兩個不同的想象運動之間的能量差異很大。在想象左手運動中，C4電極的活躍度明顯高于C3電極而當(dāng)想象右手的動作時，兩個電極處于相反的位置，變?yōu)镃3更為活躍，這樣的結(jié)果符合了運動想象中的同步去同步特性。觀察實驗結(jié)果，圖8展示了9位受試者所有實驗的平均分類正確率。

圖8 平均分類準(zhǔn)確率比較

可以發(fā)現(xiàn)，文中的分類方法可以達(dá)到 94%～98%的準(zhǔn)確率。9位受試實驗得到的平均分類準(zhǔn)確度基本上都在94.6%～96.4%，平均準(zhǔn)確率為95.9%。

最后，對比該方法與之前的競賽中不同算法的分辨率。在一樣的數(shù)據(jù)集下，該方法的分類精度相比其他三組在分類準(zhǔn)確度上有了一定提高。同時，該文只對兩個通道信號進行采集，得到了較好的分類精度，也減少了腦電采集的信道數(shù)和數(shù)據(jù)量。 表2展示了之前的算法和文中方法的分類結(jié)果。

表2 BCI競賽成績和文中結(jié)果比較

5 結(jié)束語

為了彌補傳統(tǒng)的EMD算法及CSP算法的缺點，獲得更好的分類效果，提出了EMD結(jié)合WPT、CSP的特征提取方法。先對EEG信號進行預(yù)處理，對原始信號進行濾波去噪，接著用WPT將信號分成一組窄帶信號,處理后的信號經(jīng)EMD得到多個IMF分量。不相關(guān)的固有模函數(shù)被屏蔽和消除。然后將處理后的信號的功率譜密度作為CSP濾波的特征進行計算。最后使用SVM分類。實驗結(jié)果表明，該方法在通道的減少和分類精度的提升上，相對之前的原始方法都有一定的進步。