基于多特征融合的運(yùn)動想象腦電信號識別研究

姜月，鄒任玲

上海理工大學(xué)醫(yī)療器械與食品學(xué)院，上海200093

前言

腦機(jī)接口（Brain Computer Interface,BCI）是不依賴于大腦的正常輸出通路（外周神經(jīng)系統(tǒng)和肌肉組織），即可實現(xiàn)大腦與外界環(huán)境交流的一種人機(jī)交互方式，它通過對外部設(shè)備的控制使人們的意愿得到表達(dá)。BCI 系統(tǒng)利用采集的特定模式的腦電信號（Electroencephalogram,EEG）分析識別，轉(zhuǎn)換為相應(yīng)的控制命令，實現(xiàn)預(yù)期功能。BCI能否準(zhǔn)確實現(xiàn)用戶意愿，關(guān)鍵在于腦電識別率，直接影響著最終的控制準(zhǔn)確度，而算法優(yōu)化就是提高識別率的關(guān)鍵。

Gaur 等［1］采用經(jīng)驗?zāi)B(tài)分解（Empirical Mode Decomposition,EMD）提取頻段能量等作為特征，采用線性判別分析分類，獲得了70.2%的平均正確率。針對單一特征通常分類效果一般，目前很多研究開始組合使用不同特征提取方法。Peng 等［2］采用希爾伯特-黃變換（Hilbert-Huang Transform,HHT）提取自回歸（Auto Regressive,AR）模型參數(shù)等作為特征，使用支持向量機(jī)（Support Vector Machine,SVM）分類，識別率達(dá)到91.65%。金海龍等［3］利用HHT提取邊際譜能量差，并結(jié)合EEG時域復(fù)雜度作為特征量，使用BP神經(jīng)網(wǎng)絡(luò)分類，達(dá)到87.14%的識別率。Chatterjee等［4］提取統(tǒng)計特征，小波的能量、熵及均方根，功率譜密度的平均功率、頻帶功率共26維特征向量，分別使用SVM 和多層感知器得到85.00%和85.71%的識別準(zhǔn)確率。將不同特征結(jié)合更好地體現(xiàn)了EEG的有用信息，為分類器提供更優(yōu)的輸入。

本研究提出利用EMD 分解EEG 得到內(nèi)蘊(yùn)模態(tài)函數(shù)（Intrinsic Mode Function, IMF），提取其前三階的能量值，用共空間模式（Common Spatial Pattern,CSP）得到空域特征，求取6 階AR 模型系數(shù)，將組合后的20 維時-頻-空特征向量采用主成分分析（Principle Component Analysis,PCA）降維，新的低維特征用SVM分類，并與其它算法比較，驗證本研究方法的有效性。

1 數(shù)據(jù)來源

為了客觀評價算法的有效性，本研究采用BCI2003中由Graz科技大學(xué)提供的腦電數(shù)據(jù)（Data set III）。受試者為1名25歲健康女性，任務(wù)為根據(jù)屏幕的隨機(jī)提示以想象左右手運(yùn)動的方式控制目標(biāo)左右移動。如圖1所示，t=0～2 s時，受試者處于安靜放松狀態(tài)；t=2 s時，出現(xiàn)聲音提示實驗開始，并有“+”顯示在屏幕上且保持1 s；t=3 s時，向左或向右的箭頭出現(xiàn)，受試者想象左手或右手運(yùn)動以使方塊左右移動，直至9 s結(jié)束。信號采集使用差分放大器和Ag/AgCl電極，采集通道為國際標(biāo)準(zhǔn)10-20導(dǎo)聯(lián)系統(tǒng)中的C3、Cz和C4這3個通道，如圖2所示。其中，C3和C4位于大腦初級感覺皮層運(yùn)動區(qū)，能有效反映肢體運(yùn)動或想象運(yùn)動時的相關(guān)節(jié)律變化，Cz位于大腦中央部位，作為參考電極。信號采樣率為128 Hz，進(jìn)行0.5～30.0Hz帶通濾波。

圖1 實驗時序示意圖Fig.1 Experimental sequence diagram

圖2 電極放置位置Fig.2 Position of electrodes

2 特征提取方法

本研究采用的特征提取方法是提取3 種特征值后進(jìn)行多特征融合，可以更好突出EEG 特征信息。下面將特征提取各部分算法的實現(xiàn)分別介紹。

2.1 EMD

EMD 是一種針對非線性、非平穩(wěn)信號的時頻分析方法，EEG 正是這種信號。與同為時頻分析的小波變換相比，它不需要選取基函數(shù)與分解層數(shù)，具有自適應(yīng)性［5］。它從信號時間尺度特征出發(fā)，進(jìn)行平穩(wěn)化處理，將原始信號分解成一系列不同尺度的單一分量，即IMF［6-7］。IMF 應(yīng)滿足兩個條件：（1）在整個序列中極值點(diǎn)數(shù)量與過零點(diǎn)數(shù)量須相等或至多相差一個；（2）在任意時間點(diǎn)上、下包絡(luò)線均值為零，即信號關(guān)于時間軸局部對稱［8］。

EMD對腦電信號時間序列x(t) 的處理步驟如下：

（1）確定x(t) 的所有局部極值點(diǎn)，用三次樣條插值函數(shù)進(jìn)行擬合，由局部極大值和極小值點(diǎn)分別形成x(t)的上包絡(luò)線m+(t)和下包絡(luò)線m-(t)，并計算兩包絡(luò)線的均值為：

（3）計算x(t)減去c1(t)后的殘差r1(t)，將r1(t)作為新序列，重復(fù)步驟（1）和步驟（2），可依次得到第二階、第三階IMF分量等，即c2(t),c3(t),…,cn(t)，直至rn(t)為單調(diào)函數(shù)且無法繼續(xù)從中提取分量，分解則結(jié)束。此時，可將x(t)表示為所有IMF 分量和殘差之和，即：

其中，n為得到的IMF個數(shù)。

EEG信號被分解成一系列IMF后，突出了不同時間尺度的局部特征，且隨階數(shù)增大所代表的頻率逐漸降低。一般分解EEG后的IMF數(shù)目為5～8個［9］，但并非所有都包含重要信息，只需選擇部分IMF即可。

2.2 CSP

CSP主要針對二分類問題，非常適合左右手運(yùn)動想象的特征構(gòu)造，其基本原理是對兩個協(xié)方差矩陣同時對角化，并提取兩種任務(wù)的空間成分，構(gòu)建兩類空間濾波器，使兩類差異最大化［10］。實現(xiàn)步驟如下：

（1）分別計算兩類信號的協(xié)方差矩陣：

其中，Ei表示想象左手和右手的數(shù)據(jù)矩陣，為跡，表示矩陣對角線元素之和。

（4）進(jìn)行白化處理，構(gòu)建白化矩陣P：

利用P 使Ci變成如下形式：

（5）由于S1和S2具有相同特征向量，再次特征值分解，得到：

其中，B 為兩者共同的特征向量，且λ1+λ2=I（I 為單位矩陣），也就是當(dāng)其中一類特征值最大時，另一類特征值必為最小。

（6）按照特征值大小將對應(yīng)特征向量從大到小排序，取前m 個特征向量組成矩陣B1，其余為B2，則構(gòu)成一對空間濾波器分別為：

利用它對EEG 矩陣Ei濾波，則Ei被變換為Zi=Wi×Ei。

（7）將Zi作如下運(yùn)算作為特征值：

其中，p=1,2,…,2m(2 m ＜n)。將所有fp構(gòu)成最終特征值矩陣F={f1,f2,…,f2m} ，即為一組腦電特征［11-12］。

2.3 AR模型

AR 模型屬于時域分析，基本思想是：假設(shè)可用自回歸過程近似于真實的EEG，從而可用較少的參數(shù)反映較多譜信息［13］?；诖思僭O(shè)，選取合適的階次和參數(shù)，可使AR 過程最大程度逼近EEG。AR 模型的參數(shù)辨識簡單，實時性較好，因此AR 系數(shù)作為腦電特征量非常合適［14］。

對于一般隨機(jī)信號，AR模型可以表示為如下形式：

其中，y(n)是信號的第n 個采樣值，p 是模型階數(shù)，{ai} 是模型系數(shù)，r(n)為零均值白噪聲的殘差。本研究運(yùn)用Burg算法進(jìn)行系數(shù)估計，主要原理如下：

令ep(n)和bp(n)分別為AR模型的前向和后向預(yù)測誤差，Pe和Pb分別為它們的功率，具體表達(dá)式這里不進(jìn)行贅述。令兩功率之和為：

當(dāng)階數(shù)i=1～p 時，ei(n)與bi(n)有如下遞推關(guān)系：

其中，e0(n)=b0(n)=y(n)，ki為反射系數(shù)，ki=aii。

Burg 算法的目的是使Peb相對于ki最小。系數(shù)ai則可由Levinson-Durbin遞推算法求出：

選擇階數(shù)是影響AR系數(shù)特征的關(guān)鍵，本研究選擇的階數(shù)為6。

2.4 基于PCA的特征降維

上述3種算法得到的特征向量組合后維數(shù)較高，不利于分類器識別，使算法不夠準(zhǔn)確快速，因此選擇PCA 算法實現(xiàn)降維。PCA 根據(jù)方差最大化原理，用一組線性無關(guān)且相互正交的新向量來對原向量進(jìn)行表征，這組新向量即主成分，是原來向量的線性組合［9］。實現(xiàn)步驟如下［15-16］：

（1）設(shè)原始特征集為{x1,x2,…,xn}，先進(jìn)行去均值處理，得到中心化后的矩陣：

（2）計算X 的協(xié)方差矩陣C，C=XXT。對C 特征值分解，C=UλUT，得到特征值λ 和特征向量U ，λ 為對角陣，表示主成分方差的大小。

（3）將λ 排序：λ1≥λ2≥…≥λd，取前d′個特征值對應(yīng)的特征向量構(gòu)成投影矩陣W=(u1,u2,…,ud′)。其中，d′為降維后的維數(shù)，由各成分的累積貢獻(xiàn)率決定，即：

一般Crate≥85%即可。

（4）根據(jù)W 求出低維空間的新特征量F=X*W 。

經(jīng)PCA處理后的新特征最大程度代表了有用信息，剔除對應(yīng)于d-d′個特征值的那些特征向量，使樣本采樣密度增大，一定程度起到去噪作用［16］，作為下一步分類器的輸入更佳。

3 數(shù)據(jù)處理及結(jié)果

3.1 預(yù)處理

為了使數(shù)據(jù)盡可能去除其他因素的影響，具有代表性，首先對原始EEG信號預(yù)處理。Cz電極一般作為參考電極，故本研究只選取C3和C4通道，圖3為某次樣本C3和C4通道原始信號。研究表明，人們在實施或僅是想象單側(cè)手運(yùn)動時，其對側(cè)運(yùn)動感覺區(qū)的μ節(jié)律（8～13 Hz）和β節(jié)律（18～30 Hz）會出現(xiàn)幅值降低的現(xiàn)象，稱為事件相關(guān)去同步（Event-Related Desynchronization,ERD）；其同側(cè)運(yùn)動感覺區(qū)的μ節(jié)律和β節(jié)律出現(xiàn)幅值升高的現(xiàn)象，稱為事件相關(guān)同步（Event-Related Synchronization,ERS）［17］。因此濾除噪聲及非相關(guān)的大腦活動頻段，使信號更好體現(xiàn)運(yùn)動想象特征。本研究選擇FIR等波紋濾波器對包含140次信號的樣本進(jìn)行8～30 Hz帶通濾波，濾波器參數(shù)設(shè)置為：通帶截止頻率分別為8和30 Hz，阻帶截止頻率分別為6和32 Hz，通帶衰減為1 dB，阻帶衰減為50 dB，密度因子為20。圖4為濾波后的C3和C4通道信號，濾波后的EEG較好地消除噪聲信號及無關(guān)思維信號的影響。

圖3 原始EEG信號Fig.3 Original electroencephalogram(EEG)signals

由于數(shù)據(jù)采集過程中，第3秒才給予運(yùn)動想象提示，且受試者有一定反應(yīng)時間。且在圖4中，該組信號對應(yīng)類別號為1，即左手運(yùn)動想象，可看出在第500個采樣點(diǎn)后，與左手同側(cè)的C3通道幅值增大，對側(cè)的C4 通道幅值則較小，體現(xiàn)了一定的ERS/ERD 現(xiàn)象。因此本研究經(jīng)過對比后，選取9 s數(shù)據(jù)中的4～8 s作為使用。數(shù)據(jù)采樣率為128 Hz，單次實驗時長9 s，則每次實驗中每通道有1 152 個點(diǎn)，故應(yīng)截取其中的第512～1 023 個點(diǎn)。此時數(shù)據(jù)為512×2×140，其中，512為采樣點(diǎn)數(shù)，2為通道數(shù)，140為實驗次數(shù)。

圖4 濾波后的EEG信號Fig.4 Filtered EEG signals

3.2 特征提取與融合

（1）采用EMD 對預(yù)處理后的EEG 信號分解，得到一系列IMF。由于EMD 是一種自適應(yīng)分解方法，因此無須提前確定分解階數(shù)。圖5和圖6分別是C3和C4通道信號的IMF波形。

如圖5和圖6所示，信號被自適應(yīng)地分解為imf1～imf6 共6 階IMF 和殘差（res），殘差表示最終的信號趨勢情況。從圖中可以看出，C3和C4的趨勢相反。在6階IMF中，隨著信號分解，頻率不斷降低，有用信息逐漸減少，因此僅選取能量最為集中的前三階IMF進(jìn)行特征值提取。由于ERS/ERD現(xiàn)象，C3和C4 通道的能量有所差別，故使用能量值作為特征。能量的計算公式為［18］：

其中，Ej,k為第j 個通道第k 個IMF 分量的能量，c(i)是IMF的第i個值，n為長度。每次樣本可得到6維能量特征向量：

（2）使用CSP 提取空域特征值，它實質(zhì)上是構(gòu)造一組空間濾波器。對于二分類，可得到一對空間濾波器，則特征值為2 維，即每次樣本得到的特征向量為：F2={f1,f2} 。

圖5 C3通道的各階IMFFig.5 Each order intrinsic mode functions(IMF)of channel C3

圖6 C4通道的各階IMFFig.6 Each order IMF of channel C4

（3）用AR 系數(shù)作為時域特征，本研究基于Burg算法進(jìn)行模型系數(shù)估計，經(jīng)過對比后確定階數(shù)為6，則可相應(yīng)得到6 個系數(shù)。對每次樣本的C3 和C4 通道分別提取6個AR系數(shù)后，共可得到12維特征向量F3={a1,a2,a3,…,a12} ，其中，a1～a6為C3 通道的AR系數(shù)，a7～a12為C4通道的AR系數(shù)。

經(jīng)過3種算法對EEG進(jìn)行特征提取后，將每種特征向量進(jìn)行合并，共可得到20 維時-頻-空特征向量由于高維特征向量不利于分類器識別，且每種特征的重要程度不同，即對運(yùn)動想象EEG分類的貢獻(xiàn)率不同，用PCA降維以保留累積貢獻(xiàn)率大于85%的主成分特征，得到6維的主成分特征向量F'。此時的特征向量F'最大程度體現(xiàn)了腦電特征，不僅減少了分類器的輸入量，同時去除了貢獻(xiàn)程度小的冗余信息，更加利于識別。

3.3 基于SVM的分類識別

基于統(tǒng)計學(xué)習(xí)理論的SVM是一種監(jiān)督學(xué)習(xí)方法，對非線性、小樣本及高維數(shù)據(jù)可取得理想的分類效果［19］，基本思想是找到一個最優(yōu)分類超平面，使它盡可能將兩類樣本點(diǎn)分開，且使間隔最大化［20］。本研究的SVM分類模型為C-SVC，經(jīng)過對比后，選擇RBF核函數(shù)。此外要確定SVM中的兩個參數(shù)，懲罰因子C和核參數(shù)γ。本研究采用網(wǎng)格搜索法進(jìn)行參數(shù)尋優(yōu)，找到最佳的C和γ組合，分別以C和γ的取值范圍作為網(wǎng)格邊長，用合適的步進(jìn)值交叉成一系列網(wǎng)格，則每個交叉點(diǎn)就是一組C和γ值，通過對這些點(diǎn)進(jìn)行搜索，計算每一處的分類準(zhǔn)確率，以找到使準(zhǔn)確率最高時的C和γ值。在［0,1］內(nèi)對特征向量歸一化處理，首先粗略優(yōu)化，在較大范圍的搜索后將范圍設(shè)置小一些，再精細(xì)優(yōu)化，設(shè)置懲罰因子Cmin=-4，Cmax=4，核參數(shù)gmin=-4，gmax=4，步進(jìn)值Cstep和gstep均為0.5。得到如圖7所示的參數(shù)選擇結(jié)果等高線圖，可看出在不同的C和g組合時的準(zhǔn)確率。得到Cbest=1.414 2，gbest=0.707 1，由于程序在實現(xiàn)時搜索的值為2的指數(shù)，即log2c 和log2g，即實際C的最優(yōu)值為21.4142，γ 的最優(yōu)值為20.7071，此時的最大分類準(zhǔn)確率為91.9%。

圖7 SVM參數(shù)選擇結(jié)果圖（等高線圖）Fig.7 Support vector machine(SVM)parameter selection result(contour map)

3.4 結(jié)果對比分析

為驗證本研究算法的有效性，做了兩組對比實驗。第一組是與單一特征及兩兩特征結(jié)合時的識別率對比，均使用SVM 識別，得到如表1所示的結(jié)果。如表1所示，使用單一特征值時，CSP的識別率最高，EMD的識別率最低；在特征兩兩結(jié)合時，識別率均大于80%，較單一特征有所提高，其中，CSP和AR模型結(jié)合時的識別率最高，同時由于單一特征時CSP 效果最好，在兩兩結(jié)合時，與CSP結(jié)合的兩種算法較另一種都更高；而本研究采用的多特征融合方法，達(dá)到91.9%的準(zhǔn)確率，明顯高于其他6種，證明本研究特征提取算法的有效性。

表1 不同特征下的SVM識別率對比Tab.1 Comparison of SVM recognition rate based on different features

第二組對比是將本研究的特征用另外兩種分類器識別。選擇BP神經(jīng)網(wǎng)絡(luò)和概率神經(jīng)網(wǎng)絡(luò)概率神經(jīng)網(wǎng)絡(luò)進(jìn)行對比，均使用本研究特征作為輸入，BP神經(jīng)網(wǎng)絡(luò)、概率神經(jīng)網(wǎng)絡(luò)和SVM 識別率分別為85.4%、83.6%、91.9%，BP神經(jīng)網(wǎng)絡(luò)和PNN的識別率均低于本研究通過參數(shù)尋優(yōu)后的SVM，證明本研究SVM分類的有效性和可行性，說明將其與本研究的融合特征結(jié)合是相對匹配的。

4 結(jié)論

本研究對BCI2003 數(shù)據(jù)預(yù)處理后分別使用EMD、CSP、AR 模型3 種算法，得到時-頻-空特征向量，用PCA 降維后，用SVM 分類，得到91.9%的識別率。同時做了兩組對比實驗，證明本研究多特征融合及SVM分類的有效性。運(yùn)動想象腦電信號的識別作為BCI的關(guān)鍵環(huán)節(jié)，對BCI的功能準(zhǔn)確度起到重要作用，而特征值的表征能力與腦電識別率密切相關(guān)。本研究的特征提取方法較以往單一特征得到較好的效果，可用于BCI對高識別率的要求。此外該算法可繼續(xù)研究運(yùn)動想象多分類，增加算法普適性。