基于CSP-PSO-SVM的運(yùn)動想象EEG信號特征提取與分類算法

劉寶，唐雨琦,2，蔡夢迪，薄迎春，張月

(1.中國石油大學(xué)(華東)信息與控制工程學(xué)院，山東青島，266580；2.中國電波傳播研究所，山東青島，266109)

腦機(jī)接口(brain-computer interface，BCI)是一種直接通過外部設(shè)備與大腦建立連接通路的新型人機(jī)接口方式，包括信號采集，信號分析與外部設(shè)備控制3個部分[1]。BCI分為植入式和非植入式，其中非植入式BCI 采集的是大腦頭皮的腦電信號(electro-encephalogram，EEG)，并且非植入式的采集方式有無損傷、使用便捷和成本低等優(yōu)點[2]，因此，基于EEG 信號采集的BCI 技術(shù)具有很好的研究價值，在醫(yī)療領(lǐng)域、人工智能領(lǐng)域和娛樂領(lǐng)域都具有廣闊的應(yīng)用前景。BCI技術(shù)中的核心環(huán)節(jié)就是對EEG 信號進(jìn)行分析處理，主要包括對于EEG信號的特征提取和識別分類[3]。常用的EEG信號特征提取方法主要有時域分析法[4]、頻域分析法、時頻混合分析法[5]以及共空間模式(common spatial pattern，CSP)，其中以CSP的應(yīng)用最廣，效果也最好。CSP主要是利用矩陣對角化原理構(gòu)建空間濾波器，然后通過濾波器提取信號中方差最優(yōu)的特征。傳統(tǒng)的CSP 主要用于解決二分類問題，難以滿足對EEG 信號多分類的需求，因此，衍生了“一對一”CSP(pair-wise CSP，PW-CSP)[6]以及“一對多”CSP(over versus rest CSP，OVR-CSP)[7]。WU等[7]在對OVR-CSP 應(yīng)用時得到了80%的平均準(zhǔn)確率。除此之外，CHENG 等[8]提出了一種正則化CSP算法，并取得了較好效果。不過，雖然CSP具有良好的特征提取性能，但其在實際應(yīng)用過程中仍存在對偽跡敏感度較高等問題，因此，依然無法滿足現(xiàn)實中對EEG 信號分類快速性與準(zhǔn)確率的雙重要求。常用的EEG 信號分類算法主要有線性判別式(linear discriminant analysis，LDA)、神經(jīng)網(wǎng)絡(luò)算法(neural networks，NN) 和支持向量機(jī)(support vector machine，SVM)。其中，LDA 因為其判別函數(shù)是線性的，因此，對于不可分的特征向量線達(dá)不到理想的分類效果[9]；神經(jīng)網(wǎng)絡(luò)算法可用于解決非線性問題，但其存在學(xué)習(xí)速率低、網(wǎng)絡(luò)結(jié)構(gòu)調(diào)整困難等問題；而SVM 既能解決線性問題又可以解決非線性問題，對噪聲具有魯棒性，能很好解決非線性分類及小樣本問題，并且在EEG 分類問題中有較好的效果[10]。SCHLOGL 等[11]基于EEG分類問題，將SVM與其他分類算法進(jìn)行對比，發(fā)現(xiàn)其有較好的分類效果。但是，SVM 參數(shù)調(diào)整的自適應(yīng)性與快速性還有待進(jìn)一步提高，以滿足EEG 信號分類識別中對實時性的要求。目前，對多類別EEG信號的特征提取依然十分困難，并且對EEG 信號的分類識別也存在準(zhǔn)確性較低等問題。為此，本文作者對EEG 信號特征進(jìn)行研究分析，提出一種基于CSP-PSO-SVM 算法的EEG信號特征提取與分類算法，以提高EEG 信號的特征提取準(zhǔn)確性和分類效果；最后，采用2008 年BCI競賽的數(shù)據(jù)集2A對該算法進(jìn)行驗證。

1 EEG信號生理學(xué)基礎(chǔ)

EEG 信號圖是指以從人的頭皮表層采集到的電信號以電位作為縱軸，以時間作為橫軸，所得到的電位與時間的相互關(guān)系平面圖。由于人的大腦活動十分復(fù)雜，并且大腦各神經(jīng)細(xì)胞的分布也不是均勻的，所以，由單個電極采集的EEG 信號其實是大腦不同分區(qū)和各種類型的腦電波所疊加而成的信號。EEG信號具有頻率、波幅和位相3個基本特征[12]。在大腦進(jìn)行不同的思維活動時，EEG信號的基本特征會產(chǎn)生相應(yīng)變化。EEG 信號是通過非植入式的BCI系統(tǒng)采集得到的，由于需要采集電極與頭皮進(jìn)行有效接觸以采集EEG 信號，因此電極的導(dǎo)聯(lián)組合十分重要[13]。目前，在EEG 研究中電極安放通常采用國際腦電圖學(xué)會標(biāo)定的10-20電極導(dǎo)聯(lián)定位標(biāo)準(zhǔn)，其中，10-20表示每個電極與鄰近電極的距離為大腦垂直中心線(鼻根至枕外粗隆的前后連線)或水平中心線(左右耳廓UI 高點之間的連線)的10%或20%。電極名稱用電極位置所在的頭部分區(qū)的英文單詞首字母加上國際通用的阿拉伯?dāng)?shù)字來表示，如表1所示。表1中，以中心線為基準(zhǔn)，左側(cè)部位電極名稱中的數(shù)字為奇數(shù)，右側(cè)為偶數(shù)，正中位為零點(zero)，如A1 和A2 分別為左、右耳垂電極名稱。接近中線用較小數(shù)字，較外側(cè)用較大數(shù)字。10-20電極導(dǎo)聯(lián)定位方法示意圖如圖1所示。

對不同刺激所產(chǎn)生的EEG 信號的分析處理是目前的研究熱點，不同的刺激也會產(chǎn)生不同的EEG 信號特征，其中以對運(yùn)動想象刺激所產(chǎn)生的EEG 信號特征研究居多。由于通過人本身的運(yùn)動想象刺激，誘發(fā)大腦感覺運(yùn)動皮層的腦電發(fā)生節(jié)律性變化的特征較明顯，可以通過對不同運(yùn)動想象刺激下的EEG 特征進(jìn)行判別，以實現(xiàn)不依賴外部刺激的自發(fā)的外部設(shè)備控制，這具有一定的實際意義和研究價值，因此，運(yùn)動想象EEG 信號的分析研究及應(yīng)用也成為了研究熱點。本文對EEG信號特征提取與分類算法的研究也建立在對運(yùn)動想象EEG信號特征的研究基礎(chǔ)之上。

表1 10-20電極導(dǎo)聯(lián)標(biāo)準(zhǔn)中不同部位的電極名稱Table 1 Electrode name for different parts in the 10-20 electrode lead standard

研究表明，當(dāng)人主動進(jìn)行肢體動作以及運(yùn)動想象時，會激活大腦皮層EEG 對應(yīng)區(qū)域并改變EEG頻帶能量。與運(yùn)動想象相關(guān)的頻段為α波、中間快波及β波，即8～30 Hz，其中以8～12 Hz的μ節(jié)律以及18～26 Hz 的β節(jié)律為主。在大腦做單側(cè)的運(yùn)動或運(yùn)動想象時，大腦同側(cè)運(yùn)動感知區(qū)μ節(jié)律和β節(jié)律幅度會升高，即產(chǎn)生事件相關(guān)同步(eventrelated synchronization，ERS)現(xiàn)象，而對側(cè)μ/β節(jié)律幅度會降低，即出現(xiàn)事件相關(guān)去同步(eventrelated synchronization，ERD)現(xiàn)象[14]。因此，ERD/ERS 現(xiàn)象可以作為對運(yùn)動想象信號進(jìn)行識別的有效依據(jù)。本文主要針對想象左手、想象右手、想象舌頭和想象腳這4類EEG信號進(jìn)行分析和研究。

圖1 10-20電極導(dǎo)聯(lián)定位方法示意圖Fig.1 Schematic diagram of 10-20 electrode lead positioning method

2 基于CSP-PSO-SVM 的運(yùn)動想象EEG信號特征提取與分類算法

為提高多類EEG 信號識別準(zhǔn)確率，通過對EEG 信號特征的研究，提出一種基于CSP-PSOSVM 的特征提取與分類算法，該算法主要包括預(yù)處理、特征提取與分類3 個部分，整體結(jié)構(gòu)如圖2所示。

圖2 EEG信號特征提取及分類算法整體結(jié)構(gòu)圖Fig.2 Overall structure of EEG signal feature extraction and classification algorithm

該算法各部分的具體組成及處理流程如下。

1)預(yù)處理部分。該部分首先提取EEG 信號中的刺激反應(yīng)區(qū)間信號，剔除數(shù)據(jù)中明顯的干擾，然后利用小波包變換進(jìn)行偽跡濾波以及EEG 特征頻段信號提取，最后按照一定的比例將實驗數(shù)據(jù)分成測試數(shù)據(jù)與訓(xùn)練數(shù)據(jù)。

2)特征提取部分。該部分首先利用訓(xùn)練數(shù)據(jù)構(gòu)造多類“一對一”的CSP 空間濾波器即“PWCSP”濾波器。然后，將測試數(shù)據(jù)與訓(xùn)練數(shù)據(jù)通過“PW-CSP”濾波器進(jìn)行濾波處理；最后，經(jīng)過Hilbert 變換和歸一化處理得到最終的EEG 數(shù)據(jù)特征。

3)特征分類部分。該部分首先用訓(xùn)練數(shù)據(jù)對應(yīng)的特征值輸入SVM進(jìn)行訓(xùn)練，且SVM的參數(shù)值由粒子群算法(particle swarm optimization，PSO)進(jìn)行全局優(yōu)化得到，然后用測試集的特征值對訓(xùn)練好的分類器進(jìn)行精度測試。

2.1 基于小波包變換的EEG信號預(yù)處理

因為EEG 信號很微弱且易受眼電、心電的偽跡干擾，需要對其進(jìn)行預(yù)處理以去除偽跡與噪聲。由于小波包變換具有無冗余、無疏漏優(yōu)勢，能實現(xiàn)對信號的局部時頻分析，可以滿足對包含大量細(xì)節(jié)成分的EEG 信號的分析需求，并且無需考慮額外的眼電信號采集通道，因此，在對比各類預(yù)處理算法后，本文最終選擇用小波包變換來實現(xiàn)特征頻段的提取。

對EEG 信號進(jìn)行小波包變換其實可以視為通過選擇一個形狀可變且能做平移的信號分析窗即小波基函數(shù)對EEG 信號進(jìn)行掃描分析。小波包函數(shù)的表達(dá)式如下[15]：

式中：un(t)為t時刻的信號在分解后的節(jié)點n處的表達(dá)式；l為位移；H0和H1為1對共軛正交的小波包濾波器，滿足H1(k)=(-1)kH0(1-k)。且當(dāng)n=0時，可得到尺度函數(shù)u0(t)=φ(t)，小波基函數(shù)u1(t)=Φ(t)。

小波包變換是一個類似二叉樹的結(jié)構(gòu)。若將小波包分解中的節(jié)點標(biāo)注為(j,p)(其中j為小波包分解的深度或者層數(shù)，p為節(jié)點的頻帶順序)，則在進(jìn)行j層小波包分解之后，第j個節(jié)點所對應(yīng)的頻帶范圍為其中，fs為采樣頻率。

對EEG 信號進(jìn)行小波包變換時，首先需要根據(jù)EEG 信號的特征選擇合適的小波基函數(shù)，然后需要確定EEG 信號的特征頻段，以確定小波包分解的層數(shù)與重構(gòu)節(jié)點。

2.2 基于CSP-Hilbert 變換的4 類運(yùn)動想象EEG信號特征提取

在對EEG 信號進(jìn)行預(yù)處理之后，需要對其進(jìn)行特征提取?；趯EG 信號本身特征的研究與分析，本文提出基于CSP-Hilbert 變換的特征提取算法，以提取出數(shù)據(jù)中歸屬于某一類別的特征分量。該算法主要包括一級特征提取單元、二級特征提取單元和終極特征提取單元3 個部分。其中，一級特征提取單元由共空間模式算法構(gòu)成，二級特征提取單元由Hilbert 變換構(gòu)成，終極特征提取單元由歸一化和算數(shù)求和處理構(gòu)成，具體流程如圖3所示。

2.2.1 一級特征提取單元的算法原理

此部分主要是利用CSP 算法得到不同類別EEG 信號的空間特征。CSP 算法的主要思想是利用矩陣的對角化特征找到1個最優(yōu)的空間因子，使得2類數(shù)據(jù)的方差值差異最大化，且類內(nèi)的方差最小，該算法多用于EEG信號特征處理[16-17]。假設(shè)X1和X2分別為2類N×T維的EEG數(shù)據(jù)(N為通道數(shù)，T為每通道的采樣點數(shù)，且有N＜T)，則有

式中：S1和S2分別為2類相互線性獨立的源信號幅值；SR為2 種類型任務(wù)下所共同擁有的源信號幅值，CR為與SR相應(yīng)的共有的空間模式；假設(shè)S1和S2分別由m1和m2個源所構(gòu)成，則C1和C2便由S1和S2相關(guān)的m1和m2個共同空間模式組成。

由于每個空間模式都是1個N×1維的向量，下面用這個向量來表示單個的源信號所引起的信號在N個源即N個通道上的分布權(quán)重。CSP算法的目標(biāo)是設(shè)計空間濾波器F1和F2得到空間因子W。經(jīng)過CSP空間濾波的EEG數(shù)據(jù)可以表示為

式中：XCSP為CSP空間濾波后所得到的EEG特征數(shù)據(jù)；Xinitial為EEG原始信號幅值。

由原理可知CSP 只能解決2 類數(shù)據(jù)的分類問題，為了將多分類問題轉(zhuǎn)化為二分類問題，多采用構(gòu)建“一對一”或者“一對多”共空間濾波器即PW-CSP 或者OVR-CSP 濾波器，其中，PWCSP濾波器效果相對較好，因此，本文中采用構(gòu)建PW-CSP 濾波器的方式進(jìn)行特征提取。對于d類EEG 數(shù)據(jù)，則需要建立[d(d-1)/2]組PW-CSP 濾波器，每類對應(yīng)其中(d－1)組濾波器。以4類EEG數(shù)據(jù)為例，需要分別構(gòu)建6組PW-CSP濾波器，每一類對應(yīng)其中的3 組濾波器，如圖4 所示，其中，W1～W6為6個空間濾波器。

先以訓(xùn)練數(shù)據(jù)為數(shù)據(jù)輸入，對4 類EEG 數(shù)據(jù)構(gòu)建PW-CSP 濾波器后，將訓(xùn)練數(shù)據(jù)和測試數(shù)據(jù)統(tǒng)一進(jìn)行濾波處理即可分別得到訓(xùn)練數(shù)據(jù)和測試數(shù)據(jù)的一級空間特征。

圖3 EEG信號特征提取算法Fig.3 EEG signal feature extraction algorithm

圖4 一級特征提取單元示意圖Fig.4 Schematic diagram of the 1st-level feature extraction unit

2.2.2 二級特征提取單元的算法原理

一級特征提取單元提取到的是不同類別EEG信號的空間特征，但這些特征不能完全滿足EEG信號的分類要求，因此，需要再進(jìn)行下一步處理。根據(jù)帕斯維爾定理，信號時域總能量與頻域總能量存在相等的關(guān)系，因此，可以通過Hilbert 變換求取信號包絡(luò)，即通過Hilbert 變換將信號頻域能量的變化轉(zhuǎn)化為時域能量變化，將ERDERS 現(xiàn)象映射到時域特征上。Hilbert 變換可理解為將一個原始信號和另一個信號h(t)=1/(2π)做卷積的結(jié)果，其本質(zhì)上是1 個轉(zhuǎn)向器，即將信號相位推遲90°。Hilbert變換也是研究ERD/ERS現(xiàn)象的重要手段，其表達(dá)式為

式中：x(t)為原始信號幅值；τ為時移變量；?(t)為原始信號的Hilbert變換結(jié)果。將經(jīng)過第1部分處理之后的信號進(jìn)行Hilbert 變換，得到如式(4)所示的Hilbert 變換輸出，進(jìn)而得到EEG 信號在時域的二級能量特征。同樣以4 類EEG 數(shù)據(jù)為例，其Hilbert變換過程如圖5所示。

2.2.3 終極特征提取單元的算法原理

在經(jīng)過前2個特征提取單元的處理之后，還需要進(jìn)行歸一化處理，以降低基線不穩(wěn)、危機(jī)干擾對信號特征的影響。歸一化處理是將信號以時間序列為基準(zhǔn)，將信號幅值轉(zhuǎn)化為0～1之間。同樣以4類EEG數(shù)據(jù)為例，若某一時刻一組訓(xùn)練數(shù)據(jù)經(jīng)過二級特征提取單元的處理之后得到的二級能量特征為6組向量則其歸一化處理為

圖5 二級特征提取單元示意圖Fig.5 Schematic diagram of the 2nd-level feature extraction unit

2.3 基于PSO-SVM 的4 類運(yùn)動想象EEG 信號分類算法

在獲得EEG 信號的特征之后，需要利用分類器進(jìn)行分類識別。由于SVM 能夠解決機(jī)器學(xué)習(xí)中普遍存在的局部極值問題和高維度問題，并具有很強(qiáng)的泛化能力，同時，在解決小樣本問題、非線性問題和模式識別問題等方面也有很好的效果，因此，本文采用SVM 作為非線性EEG 信號分類器。實踐表明，SVM 的性能與核函數(shù)類型及其參數(shù)選擇以及支持向量機(jī)算法中的容忍度因子密切相關(guān)，因此，核函數(shù)的選擇以及容忍度因子的參數(shù)調(diào)整十分重要。本文采用粒子群優(yōu)化(particle swarm optimization，PSO)的方式對支持向量機(jī)中的參數(shù)選擇進(jìn)行優(yōu)化調(diào)整并以交叉驗證(cross validation，CV)的方式衡量算法的泛化能力，以提高算法的自適應(yīng)性與快速性。

圖6 終級特征提取單元示意圖Fig.6 Schematic diagram of the final-level feature extraction unit

以SVM實現(xiàn)EEG信號特征分類識別的基本思想如下：首先通過1 個核函數(shù)構(gòu)成的非線性映射Φ(x)將y類EEG信號特征x映射到高維空間，然后找到1個超平面使得所有點到超平面的距離大于一定值[18]，對最優(yōu)超平面的構(gòu)造歸結(jié)為凸二次規(guī)劃問題：

式中：C為容忍度因子；xi為EEG 信號的特征向量；yi為xi所對應(yīng)的類別；ω為構(gòu)造的超平面的法向量；b為超平面的偏置向量；ξ為松弛向量。

該凸二次規(guī)劃問題可轉(zhuǎn)換為對其對偶規(guī)劃問題進(jìn)行求解：

式中：K(x,xi)為核函數(shù)；αi(i=1,2,…,m)為拉格朗日向量，若αi＞0，則稱對應(yīng)EEG信號的特征向量xi為支持向量。對式(7)求解可得到最終的超平面和決策函數(shù)分別如式(8)和式(9)所示：

由SVM 原理可知，最后得到的決策函數(shù)為1個符號函數(shù)，因此，只能解決二分類的問題。而對于EEG 信號識別這類多分類問題，通常采用直接法和間接法2種方法來求解。本文采用間接法中的“一對一”法(one-versus-one，OVO-SVMs)。OVO-SVMs 通過在任意2 類樣本之間設(shè)計1 個SVM 分類器，對于d類EEG 數(shù)據(jù)，則需要建立個SVM 分類器，每類對應(yīng)其中(d-1)個SVM 分類器，然后以投票法來確定樣本最終類別，例如對于第1類和第2類的SVM分類器，若將信號輸入判定為第1 類，則第1 類的票數(shù)加1，反之則第2 類的票數(shù)加1。以44 類EEG 信號為例，SVM的設(shè)計與投票過程如圖7所示，圖7中SVMi-j代表第i類和第j類之間的分類器。

由于EEG 信號具有非線性、高維度等特征，因此，本文選擇映射維度不受限制的徑向基函數(shù)(radial basis function，RBF)作為核函數(shù)，其表達(dá)式為

由式(7)和式(10)可知：RBF核函數(shù)中的函數(shù)寬度γ以及SVM 算法中的容錯率因子C對于最終的模型建立極為重要，其直接決定SVM 的性能。因此，需要找到(γ，C)的最優(yōu)參數(shù)以得到最佳SVM分類模型。對于參數(shù)(γ，C)的優(yōu)化確定，一般采用網(wǎng)格搜索或者經(jīng)驗法選取，但在算法訓(xùn)練速率上還有待提高，因此，本文采用PSO 實現(xiàn)全局搜索優(yōu)化。

PSO主要是來源于鳥類覓食的最優(yōu)策略[18]，對于SVM 參數(shù)優(yōu)化來說，鳥類覓食可以類比為對參數(shù)(γ，C)最優(yōu)參數(shù)的尋找。首先將鳥抽象為沒有質(zhì)量和體積的粒子，并延伸到多維空間，即初始化一群隨機(jī)粒子作為隨機(jī)解。假設(shè)延伸至D維，某個粒子I在D維空間的位置表示為矢量Xi=(x1,x2,…，xD)，飛行速度表示為矢量Vi=(v1，v2，…，vD)。每個粒子都有一個目標(biāo)函數(shù)決定的適應(yīng)值(fitness value)，在每次迭代的飛行過程中積累自己的飛行經(jīng)驗來確定下一步運(yùn)動。本文中的目標(biāo)函數(shù)為以k-CV(k-folder cross-validation )方法進(jìn)行交叉驗證下的平均識別率。k-CV 是將數(shù)據(jù)集分為k個子集，每個子集輪流作為一次驗證集，其余作為訓(xùn)練集，進(jìn)行k次訓(xùn)練驗證，將k次驗證的平均識別率作為最終的交叉驗證識別率。

圖7 SVM分類過程示意圖Fig.7 Schematic diagram of SVM classification process

式中：A為最終的交叉驗證識別率；k為子集的數(shù)量；ai為第i次的驗證識別率。

可見，在由PSO 優(yōu)化的SVM 參數(shù)選擇中，是通過粒子跟蹤發(fā)現(xiàn)的最佳參數(shù)(γbest,m，Cbest,m)和現(xiàn)在的位置Xm，計算出交叉驗證識別率，然后更新找尋方向，如式(12)和式(13)所示。

式中：m=1,2,…,M；M為算法迭代的次數(shù)；ωp為慣性因子且非負(fù)；Vm+1為粒子更新后的速度；Vm為粒子當(dāng)前的速度；rand(·)為介于(0,1)之間的隨機(jī)數(shù)；Xm+1為粒子更新后的位置；Xm為粒子的當(dāng)前位置；c1和c2為學(xué)習(xí)因子。

整體的算法流程圖如圖8 所示。由圖8 可知，算法可分為5步：

1)讀取樣本數(shù)據(jù)，隨機(jī)產(chǎn)生1組(γ，C)作為初始位置；

2)把整個樣本平均分割為k個互補(bǔ)包含的子集s1,s2,…,sk；

3)根據(jù)當(dāng)前的(γ，C)訓(xùn)練SVM 計算識別交叉驗證下的平均識別率，即計算k個子集輪流作為訓(xùn)練集，其余子集作為驗證集時，依據(jù)所得到的k次驗證識別率的平均值以及式(12)和式(13)更新位置。

① 初始化i=1；

② 子集留作檢驗集，其余子集合并起來作為訓(xùn)練集，訓(xùn)練SVM；

③ 計算i子集的準(zhǔn)確率，直到i=k+1；

圖8 PSO優(yōu)化SVM參數(shù)算法流程圖Fig.8 PSO optimized SVM parameter algorithm flow chart

4)以交叉驗證的平均識別率作為適應(yīng)值，并記憶個體與群體所對應(yīng)的最佳位置為(γbest，Cbest)，根據(jù)式(12)和式(13)搜尋更好的(γ，C)；

5) 重復(fù)步驟2)～4)，直到滿足最大迭代次數(shù)為止。

3 實驗結(jié)果對比分析

3.1 實驗條件與數(shù)據(jù)來源

本文采用離線數(shù)據(jù)集對所提出的算法進(jìn)行實驗驗證，包括3個步驟：EEG數(shù)據(jù)預(yù)處理、特征提取與識別分類。仿真計算機(jī)配置條件為：CPU，Intel(R) Core(TM) i5 - 5200U(2.20GHz)； RAM，8.0 GB；操作系統(tǒng)，Windows 10；仿真環(huán)境為MATLABR2018b。

采用2008 年BCI 競賽中的標(biāo)準(zhǔn)4 類運(yùn)動想象EEG 數(shù)據(jù)集2A[19]來驗證本文所設(shè)計的算法效果。該數(shù)據(jù)集EEG 信號由22 路Ag/AgCl 電極測量，電極分布參照國際10-20導(dǎo)聯(lián)標(biāo)準(zhǔn)。采集時以左乳頭作為參考電位，右乳頭作為地，采樣頻率為250 Hz。共對9名被試(T01，T02，…，T09)進(jìn)行4類運(yùn)動想象的試驗，每個人進(jìn)行2個獨立周期的試驗，每個周期記錄288 組數(shù)據(jù)，各類別72 組。在仿真實驗時，用其中一個周期的采集數(shù)據(jù)作為訓(xùn)練集，另一個周期的數(shù)據(jù)作為測試集。實驗范例如圖9 所示，被試在每次試驗的3～6 s 的時間內(nèi)進(jìn)行運(yùn)動想象。

圖9 試驗范例定時方案Fig.9 Timing scheme of paradigm

3.2 4類運(yùn)動想象EEG信號特征

當(dāng)人進(jìn)行運(yùn)動想象時會出現(xiàn)ERD/ERS現(xiàn)象[20]，因此，ERD/ERS 可作為運(yùn)動想象EEG 信號的識別依據(jù)。但是，ERD/ERS 現(xiàn)象并非出現(xiàn)在所有的頻段。為了便于后期實驗的特征提取，需要對ERD/ERS 現(xiàn)象出現(xiàn)的頻段進(jìn)行研究。為此，本文繪制出被試T01的一組4類運(yùn)動想象的對數(shù)功率譜密度以及其在頻率為6.6，13.2，19.7 和26.3 Hz 的腦地形圖，如圖10所示。

由圖10 可見：ERD/ERS 現(xiàn)象在頻率為13.2，19.7 和 26.3 Hz 時較明顯。所以，ERD/ERS 現(xiàn)象主要出現(xiàn)在α波和β波，并且由文獻(xiàn)[21]可知 EEG 信號的頻譜在9～13 Hz的變化最大。因此，在預(yù)處理時選擇提取EEG中9～13 Hz頻段的成分。

3.3 預(yù)處理結(jié)果與分析

由于該實驗數(shù)據(jù)為連續(xù)采集，因此，要先根據(jù)實驗范例提取出刺激之后3 s反應(yīng)時間內(nèi)的EEG信號，在提取出之后發(fā)現(xiàn)實驗數(shù)據(jù)中存在部分缺失值。9 名被試的缺失值數(shù)量見表2，其中被試T04和被試T09的缺失值數(shù)量最多。由于缺失值會影響最后算法的實現(xiàn)效果，因此，本文采取均值插補(bǔ)的方式對缺失值進(jìn)行處理。

在處理EEG 信號中的明顯偽跡之后，可以依據(jù)提出的基于小波包變換預(yù)處理算法對EEG 信號進(jìn)行預(yù)處理。此處以被試T09的1組左手的運(yùn)動想象EEG數(shù)據(jù)中C3通道的數(shù)據(jù)為例，首先采用快速傅里葉變換求取信號頻率分布的方式對原始數(shù)據(jù)進(jìn)行頻譜分析。圖11 所示為該通道的原始信號時域波形與單邊頻譜圖。由圖11 可見：原始信號時域波形雜亂無章，無法直觀地看出信號特征，在頻域中頻率范圍主要集中在30 Hz以內(nèi)，而本文所需求的信號頻段為9～13 Hz。因此，需要對其進(jìn)行下一步處理。

圖10 4類運(yùn)動想象對數(shù)功率譜密度Fig.10 Four types of motion imaginary logarithmic power spectral density

表2 每名被試的4類運(yùn)動想象EEG數(shù)據(jù)中的缺失值(NaN)數(shù)量Table 2 Number of NaN values in EEG data for each subjects with four types of exercise imaginary 個

圖11 原始信號時域波形與單邊頻譜圖Fig.11 Original signal time domain waveform and unilateral spectrum

由于實驗數(shù)據(jù)的采樣頻率均為250 Hz，且目標(biāo)頻率段為9～13 Hz。因此，本文采用較為常用的Daubechies小波基中的db6小波基函數(shù)對EEG信號進(jìn)行6層分解，分解過程如圖12所示。

圖12 6層小波包分解示意圖Fig.12 Diagram of 6-layer wavelet packet decomposition

式中：m0為轉(zhuǎn)換函數(shù)；ωm為轉(zhuǎn)換值。可以發(fā)現(xiàn)，在經(jīng)過6 層分解之后所得到的第(6,4)和(6,7)這2 個 節(jié) 點 的 頻 段 分 別 為 11.718 8～13.671 9 Hz 和9.578 2～11.718 8 Hz，已滿足對目標(biāo)頻段范圍的要求。因此，選取第(6,4)和(6,7)這2個節(jié)點進(jìn)行信號重構(gòu)之后，再合成最終的信號。

重構(gòu)信號時域波形與單邊頻譜圖如圖13所示。由圖13 可見：經(jīng)過小波包變換后，可得到原始EEG 信號中所需頻段范圍信號特征，并且所得信號在時域上的波形更光滑。

3.4 特征提取結(jié)果與分析

對預(yù)處理之后的數(shù)據(jù)，利用CSP-Hilbert 變換進(jìn)行特征提取實驗，得到每一組EEG 信號的終極EEG 信號特征。并且由算法原理可知，終極特征向量是由該組EEG 信號中歸屬于每一類的特征成分所構(gòu)成。以被試T03的特征提取過程為例，在經(jīng)過3 層特征提取之后，得到了4×750 維的特征向量，每一行即代表每一類的特征成分。表3所示為均等分布的8個采樣點對應(yīng)的特征值以及每一類特征分量之和。

將被試T03的4類運(yùn)動想象EEG信號特征以時間為橫軸，相對幅值為縱軸，得到了被試T03特征分布圖，如圖14所示。由圖14可以更直觀地看出特征成分的分布關(guān)系。由表3 和圖14 可知：信號特征中某類的特征分量中最小的分量與信號類別相對應(yīng)。但是在具體應(yīng)用時，會存在特征分量不足以進(jìn)行區(qū)分的情況，比如有的信號存在最小的特征分量并不對應(yīng)所屬類別，因為其余特征分量在某一點存在遠(yuǎn)小于其他特征分量的極值點。因此，還需要引入分類器加以判別區(qū)分，進(jìn)一步提高EEG信號分類識別的準(zhǔn)確率。

圖13 節(jié)點(6,4)和節(jié)點(6,7)合成信號的時域波形圖與單邊頻譜Fig.13 Time domain waveform and unilateral spectrum of composite signal of node(6,4)and node(6,7)

表3 被試T03的部分特征值及其特征分量之和Table 3 Some eigenvalues and the sum of characteristic components of subject T03

3.5 分類結(jié)果與分析

本文采用基于PSO 優(yōu)化的SVM 分類器進(jìn)行分類識別。首先將4×750維的特征向量中的每一行進(jìn)行算術(shù)求和，即得到一組4×1維的特征向量作為最終的輸入特征。然后將訓(xùn)練集輸入分類器中進(jìn)行訓(xùn)練，得到每名被試的(γbest，Cbest)。PSO 參數(shù)設(shè)置如表4 所示。最終優(yōu)化得到的每名被試的(γbest，Cbest)以及交叉驗證識別率，如表5所示。圖15所示為被試T03和T09訓(xùn)練時PSO的平均適應(yīng)度。

圖14 被試T03的4類運(yùn)動想象EEG信號特征Fig.14 Characteristics of four types of motion imagnation EEG signals of subject T03

表4 PSO參數(shù)設(shè)置Table 4 PSO parameter setting

表5 每名被試經(jīng)PSO優(yōu)化得到的γbest和CbestTable 5 γbestand Cbest optimized by PSO for each subject

在完成訓(xùn)練之后，將測試數(shù)據(jù)的特征向量輸入SVM 分類器進(jìn)行分類。分類結(jié)果采用目前通用的Kappa系數(shù)K進(jìn)行衡量：

式中：p0為準(zhǔn)確率。Kappa系數(shù)越高，則對應(yīng)的識別率越高。將本文提出的算法分類結(jié)果與BCI競賽第1 名算法分類結(jié)果[19]及文獻(xiàn)[22]中算法分類結(jié)果進(jìn)行對比，結(jié)果如表6 所示。由表6 可知：在BCI競賽中，第1名使用帶通濾波器去除偽記進(jìn)行預(yù)處理，再通過構(gòu)建PW-CSP 濾波提取特征，最后通過樸素貝葉斯Parzen 窗口分類器進(jìn)行分類，其分類結(jié)果的Kappa 系數(shù)平均值為0.57，最大值為0.77，最小值為0.27；文獻(xiàn)[22]采用KNN 進(jìn)行分類，其Kappa 系數(shù)平均值為0.59，最大值為0.76，最小值為0.21；而本文采用的CSP-PSO-SVM 方法所得Kappa 系數(shù)平均值為0.60，最大值為0.91，最小值為0.34。由此可見，基于CSP-PSO-SVM的算法具有一定的優(yōu)勢。

圖15 PSO迭代優(yōu)化適應(yīng)度曲線Fig.15 PSO iterative optimization fitness curves

表6 本文算法與其他算法分類結(jié)果對比Table 6 Comparison of classification results from the proposed algorithm and other algorithms

4 結(jié)論

1)提出了一種基于CSP-PSO-SVM 的分類算法，該算法包括基于CSP-Hilbert 變換的特征提取以及基于PSO優(yōu)化的SVM算法的分類識別2部分。

2)與其他算法相比，本文方法提高了EEG 信號分類識別的準(zhǔn)確率，最高可達(dá)到93.07%。