劉寶，唐雨琦,2，蔡夢(mèng)迪，薄迎春，張?jiān)?/p>

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

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

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

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

對(duì)不同刺激所產(chǎn)生的EEG 信號(hào)的分析處理是目前的研究熱點(diǎn)，不同的刺激也會(huì)產(chǎn)生不同的EEG 信號(hào)特征，其中以對(duì)運(yùn)動(dòng)想象刺激所產(chǎn)生的EEG 信號(hào)特征研究居多。由于通過(guò)人本身的運(yùn)動(dòng)想象刺激，誘發(fā)大腦感覺(jué)運(yùn)動(dòng)皮層的腦電發(fā)生節(jié)律性變化的特征較明顯，可以通過(guò)對(duì)不同運(yùn)動(dòng)想象刺激下的EEG 特征進(jìn)行判別，以實(shí)現(xiàn)不依賴外部刺激的自發(fā)的外部設(shè)備控制，這具有一定的實(shí)際意義和研究?jī)r(jià)值，因此，運(yùn)動(dòng)想象EEG 信號(hào)的分析研究及應(yīng)用也成為了研究熱點(diǎn)。本文對(duì)EEG信號(hào)特征提取與分類算法的研究也建立在對(duì)運(yùn)動(dòng)想象EEG信號(hào)特征的研究基礎(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)人主動(dòng)進(jìn)行肢體動(dòng)作以及運(yùn)動(dòng)想象時(shí)，會(huì)激活大腦皮層EEG 對(duì)應(yīng)區(qū)域并改變EEG頻帶能量。與運(yùn)動(dòng)想象相關(guān)的頻段為α波、中間快波及β波，即8～30 Hz，其中以8～12 Hz的μ節(jié)律以及18～26 Hz 的β節(jié)律為主。在大腦做單側(cè)的運(yùn)動(dòng)或運(yùn)動(dòng)想象時(shí)，大腦同側(cè)運(yùn)動(dòng)感知區(qū)μ節(jié)律和β節(jié)律幅度會(huì)升高，即產(chǎn)生事件相關(guān)同步(eventrelated synchronization，ERS)現(xiàn)象，而對(duì)側(cè)μ/β節(jié)律幅度會(huì)降低，即出現(xiàn)事件相關(guān)去同步(eventrelated synchronization，ERD)現(xiàn)象[14]。因此，ERD/ERS 現(xiàn)象可以作為對(duì)運(yùn)動(dòng)想象信號(hào)進(jìn)行識(shí)別的有效依據(jù)。本文主要針對(duì)想象左手、想象右手、想象舌頭和想象腳這4類EEG信號(hào)進(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)動(dòng)想象EEG信號(hào)特征提取與分類算法

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

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

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

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

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

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

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

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

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

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

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

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

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

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

2.2.1 一級(jí)特征提取單元的算法原理

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

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

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

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

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

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

圖3 EEG信號(hào)特征提取算法Fig.3 EEG signal feature extraction algorithm

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

2.2.2 二級(jí)特征提取單元的算法原理

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

可見(jiàn)，在由PSO 優(yōu)化的SVM 參數(shù)選擇中，是通過(guò)粒子跟蹤發(fā)現(xiàn)的最佳參數(shù)(γbest,m，Cbest,m)和現(xiàn)在的位置Xm，計(jì)算出交叉驗(yàn)證識(shí)別率，然后更新找尋方向，如式(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)把整個(gè)樣本平均分割為k個(gè)互補(bǔ)包含的子集s1,s2,…,sk；

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

① 初始化i=1；

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

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

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

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

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

3 實(shí)驗(yàn)結(jié)果對(duì)比分析

3.1 實(shí)驗(yàn)條件與數(shù)據(jù)來(lái)源

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

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

圖9 試驗(yàn)范例定時(shí)方案Fig.9 Timing scheme of paradigm

3.2 4類運(yùn)動(dòng)想象EEG信號(hào)特征

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

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

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

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

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

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

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

圖11 原始信號(hào)時(shí)域波形與單邊頻譜圖Fig.11 Original signal time domain waveform and unilateral spectrum

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

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

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

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

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

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

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

圖13 節(jié)點(diǎn)(6,4)和節(jié)點(diǎn)(6,7)合成信號(hào)的時(shí)域波形圖與單邊頻譜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)行分類識(shí)別。首先將4×750維的特征向量中的每一行進(jìn)行算術(shù)求和，即得到一組4×1維的特征向量作為最終的輸入特征。然后將訓(xùn)練集輸入分類器中進(jìn)行訓(xùn)練，得到每名被試的(γbest，Cbest)。PSO 參數(shù)設(shè)置如表4 所示。最終優(yōu)化得到的每名被試的(γbest，Cbest)以及交叉驗(yàn)證識(shí)別率，如表5所示。圖15所示為被試T03和T09訓(xùn)練時(shí)PSO的平均適應(yīng)度。

圖14 被試T03的4類運(yùn)動(dòng)想象EEG信號(hào)特征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)練之后，將測(cè)試數(shù)據(jù)的特征向量輸入SVM 分類器進(jìn)行分類。分類結(jié)果采用目前通用的Kappa系數(shù)K進(jìn)行衡量：

式中：p0為準(zhǔn)確率。Kappa系數(shù)越高，則對(duì)應(yīng)的識(shí)別率越高。將本文提出的算法分類結(jié)果與BCI競(jìng)賽第1 名算法分類結(jié)果[19]及文獻(xiàn)[22]中算法分類結(jié)果進(jìn)行對(duì)比，結(jié)果如表6 所示。由表6 可知：在BCI競(jìng)賽中，第1名使用帶通濾波器去除偽記進(jìn)行預(yù)處理，再通過(guò)構(gòu)建PW-CSP 濾波提取特征，最后通過(guò)樸素貝葉斯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。由此可見(jiàn)，基于CSP-PSO-SVM的算法具有一定的優(yōu)勢(shì)。

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

表6 本文算法與其他算法分類結(jié)果對(duì)比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算法的分類識(shí)別2部分。

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