閆秘，白雪梅，張晨潔，郭家赫

（長春理工大學 電子信息工程學院，長春 130022）

腦-機接口（brain-computer interface，BCI）是一種將大腦活動轉換為計算機或其他設備命令的通信系統(tǒng)，可以代替人的肢體或言語器官實現(xiàn)人與外界的交流以及對外部環(huán)境的控制［1］。有很多種生理電信號可以應用在BCI系統(tǒng)中，其中P300信號以其成分特征鮮明，易于識別且不需要受試者進行訓練便可得到的優(yōu)點受到了國內(nèi)外研究人員的青睞。P300是事件相關電位（Event Related Potentials，ERP）中重要的內(nèi)源性成分，可以清楚地反映出大腦認知過程中神經(jīng)電位的生理變化。當人受到小概率事件刺激并進行認知加工后約300～500 ms時會產(chǎn)生一個可以在頭皮上記錄到的區(qū)別于未受刺激腦電信號（electroencephalography，EEG）的波峰。通過識別EEG中是否含有P300成分，可以判斷大腦的真實意圖，以便于后續(xù)工作的進行，如外周神經(jīng)和肌肉損傷的患者應用腦電信號驅(qū)動輪椅運行和喪失語言能力的患者通過腦電信號進行字符拼寫等。

近年來，隨著科技的發(fā)展及社會需求的增加，國內(nèi)外研究人員提出了很多算法來提高對P300的分類準確率。特征提取是腦電信號預處理、特征提取與分類這三個過程中最重要的環(huán)節(jié)，很大程度上決定了最終的分類準確率，因此，找到合適的方法對P300特征進行提取是十分重要的。常見的特征提取方法有共空間模型（Common Spatial Pattern，CSP）法［2］、自適應回歸模型（Adaptive Autoregressive Model，AAR）［3］、譜分析法、計算波形參數(shù)法、小波變換法［4］、時域能量熵和獨立成分分析法（Independent Component Analysis，ICA）［5］等。其中CSP方法僅考慮兩類任務在空間上的可分性，極大的限制了BCI的實用性；AAR的參數(shù)雖然能夠清晰地反映大腦的狀態(tài)，但該方法更傾向于對平穩(wěn)信號分析；頻譜分析法，如快速傅里葉變換（fast Fourier transforms，F(xiàn)FT）和功率譜估計等，這些算法具有簡單易實現(xiàn)的優(yōu)點，但FFT無法有尺度變化且功率譜分析分辨率有限；計算峰值、波形面積等指標作為時域特征能夠直觀的展現(xiàn)出目標波形的特點，但通常情況下這些指標在不同的波形中差異較大，無法統(tǒng)成為一個標準。目前小波變換、計算能量熵和ICA等特征提取方法受到了廣泛關注。小波變換解決了FFT無法進行尺度變換的缺點且具有低頻時頻率分辨率高的特點；時域能量熵值可以清楚地區(qū)分出信號是否存在P300波峰；ICA能夠在未知混合信號中分離出原信號且可同時處理多導信號。這些方法均單獨利用了P300時、頻、空域特征，無法充分體現(xiàn)出P300的本質(zhì)特征。

針對上述問題，本文結合信息融合中相關概念提出了一種多域特征提取方法，分別提取P300的時域特征、頻域特征和空域特征并將這些特征融合成特征集合，全面表征P300腦電信號的特征，再將這個特征集合送入支持向量機（support vector machine，SVM）中進行分類，有效提高分類準確率。

1 實驗數(shù)據(jù)描述及預處理

1.1 實驗數(shù)據(jù)描述

本文采用的數(shù)據(jù)來自于2004年BCI Competition III data set II中的P300字符拼寫實驗數(shù)據(jù)［6］。該實驗是在1988年由Farwell和Donchin提出的，數(shù)據(jù)中包含兩個受試者A和B的多次實驗數(shù)據(jù)。

實驗過程中，受試者頭部戴有64通道的腦電帽，電極按照國際10-20系統(tǒng)標準放置，如圖1所示，采樣頻率設置為240 Hz。受試者面對如圖2所示的6×6字符矩陣，該矩陣中的行和列隨機地連續(xù)閃爍，每行或每列閃爍一次，被稱為一個trail，十二個行和列均閃爍一次稱為一個block，其中包含2次目標刺激和10次非目標刺激。對于每個字符，這樣的過程重復15次，稱為一個run，為了方便理解，給出如圖3所示時序圖。在實驗開始前已經(jīng)為受試者設定了待拼寫的目標字符（矩陣左上方），當目標字符所在的行或列閃爍后受試者的EEG中就會產(chǎn)生一個P300波形，受試者默記這些行和列閃爍的次數(shù)，而對其他字符的行和列閃爍不作計數(shù)。經(jīng)過數(shù)據(jù)處理后，SVM分類算法能夠根據(jù)上述特點判斷出目標字符所在的行或列，進而推斷出目標字符。

圖1 EEG采集電極位置圖

圖2 誘發(fā)界面及其行列編號

圖3 P300拼寫實驗時序圖

1.2 腦電信號預處理

腦電信號是強隨機性、微幅值的非平穩(wěn)信號，極易受到外界因素干擾，因此采集到的腦電信號中往往慘雜著各種各樣的噪聲干擾，如工頻干擾、心電偽跡及肢體運動產(chǎn)生的偽跡等，其信噪比低［7］。預處理過程是腦電信號處理過程中的第一步，目的是降低夾雜在原始EEG中的噪聲和偽跡干擾，使特征提取算法能夠提取到反映大腦真實意圖的特征。

有關研究表明，并不是所有電極采集到的EEG都存在明顯的P300波形，而頂枕處電極組合對P300分類具有有效性［8］，所以本文采用Cz，C1，C2，Pz（對應序號：10，11，12，51）進行后續(xù)處理，這樣可以極大減少待處理的數(shù)據(jù)量。然后使用經(jīng)驗模態(tài)分解（Empirical Mode Decomposition，EMD）算法和FastICA算法對采集到的EEG進行去噪聲和偽跡的操作，圖4為去除噪聲和偽跡前后的對比圖，原始EEG中尖峰處為眼電偽跡，經(jīng)過預處理后較好地實現(xiàn)了降低噪聲和去除偽跡的目標。

圖4 去噪聲和偽跡前后對比圖

2 腦電信號特征提取

特征提取為腦電信號處理過程的中間步驟，對分類準確率的高低起著至關重要的作用。在特征提取之前對信號進行疊加平均，使P300信號更加明顯。在本實驗刺激范式下，閃爍時間持續(xù)100 ms，距離下次閃爍的時間間隔為75 ms，即兩次閃爍刺激的時間間隔是175 ms。根據(jù)P300的時間特性及刺激的時間特性，截取每次閃爍刺激后650 ms的數(shù)據(jù)段，并將相同刺激對應的數(shù)據(jù)進行15次疊加平均。

假設一次刺激產(chǎn)生的腦電信號［9］為：

式中，xi(t)為記錄到的原始腦電信號；si(t)為P300信號，單次刺激的功率為P；ni(t)為噪聲，方差為σ2，均值為0；則單次刺激的信噪比為：P/σ2。經(jīng)過N次疊加平均后：

式中，X(t)為經(jīng)過N次疊加平均后的腦電信號；S(t)為經(jīng)過N次疊加平均的P300信號。

經(jīng)過N次疊加平均后P300信號的功率依然為P，而噪聲的方差（功率）為σ2/N，則疊加平均后的信噪比為N?P/σ2，相比于未進行疊加平均的信號的信噪比提高N倍。

以Cz導聯(lián)為例，對數(shù)據(jù)分段并進行15次疊加平均后得到圖5所示的效果圖。從圖中可以看出，經(jīng)過疊加平均的信號可以將P300更加清晰的顯示出來。該操作為特征提取提供了有力幫助，使得P300信號的特征提取更為準確，為分類準確率的提高提供了保障。

圖5 包含P300的信號和不包含P300的信號分段疊加平均前后對比

2.1 時域特征提取

P300腦電信號是與認知活動過程密切相關的誘發(fā)腦電信號，在時域上有明顯的特點，即，小概率目標事件刺激受試者約300 ms后會出現(xiàn)一段在整段信號中波幅最高的正電位，如圖5所示。而時域能量熵反映的是信號能量在時域上的分布情況，可以較好地展示出含有P300成分的腦電信號與不含有P300成分的腦電信號在300 ms左右幅值的差異。根據(jù)信息論中對香濃熵的定義，時域能量熵可由式（3）計算：

從式（3）可以看出，某段信號的能量在時域上所占信號總能量的比值越大，那么該段信號的能量熵值Ei越小。受試者受到刺激后產(chǎn)生P300信號，此時間段信號的波幅明顯提升，在時域上信號的波幅與能量的大小正相關，因此該段能量占整個信號段能量的比重大，時域能量熵小。而沒有誘發(fā)出P300的腦電信號在整個信號段上波幅沒有明顯變化，能量分布較為均勻。

以Cz導聯(lián)為例，將650 ms的數(shù)據(jù)分為10段，每段65 ms，即1～65 ms的數(shù)據(jù)為第1段數(shù)據(jù)，586～650 ms的數(shù)據(jù)為第10段數(shù)據(jù)，分別計算分段后含有P300成分和不含有P300成分的腦電信號的時域能量，并根據(jù)公式（3）計算每段信號的時域能量熵，以N1-N10命名。含P300成分與不含P300成分的時域能量熵如圖6所示。

圖6 有無P300成分的時域能量熵

圖6中帶有三角形標記的曲線的第5段、第6段、第7段的時域能量熵達到了最小值，分別為2.501 3、1.013 4和3.442 2，即該段能量占整個信號段能量的比重大，對應著含有P300波幅的腦電信號。因此，選擇N5、N6、N7這三個時域能量熵值作為時域特征向量。

2.2 頻域特征提取

小波變換（Wavelet Transform，WT）是一種時頻分析方法，能夠?qū)π盘柕念l率逐層分解，直至達到目標信號所在的頻率范圍。EEG信號經(jīng)過小波分解得到不同尺度下的高頻細節(jié)系數(shù)與低頻近似系數(shù)，這些系數(shù)反映著EEG在特定頻段的特征，由于P300信號頻率范圍低，所以用低頻近似系數(shù)作為頻域特征，而小波變換具有在低頻部分頻率分辨率高的特點，因此非常適合處理EEG這類非平穩(wěn)的隨機信號。

信號f(x)的離散小波變換和逆變換有如下定義：

式中，f(x)為待處理的信號；dj,k為小波系數(shù)；為母小波；j為分解尺度；k為時間平移量。

考慮到算法運行速度等現(xiàn)實問題，本文選擇Mallat算法對信號進行有限層分解，母小波選擇與P300波形相近的Daubechies4小波。有關研究表明，P300信號的頻率范圍為2～12 Hz，甚至更低頻段［10］。而實驗數(shù)據(jù)的采樣頻率是240 Hz，因此需要對原始信號進行4層小波分解，得到0～15 Hz的低頻段信號，即包含P300的頻段，從這個頻段提取的低頻近似系數(shù)特征對應著P300的頻域特征。以Cz導聯(lián)為例，最終得到16個低頻近似系數(shù)，命名為X1-X16。

圖7 有無P300成分的小波近似系數(shù)

如圖7所示，帶有圓形標記的曲線表示的是包含P300成分的小波近似系數(shù)，星形標記的曲線表示的是不包含P300成分的小波近似系數(shù)。因為小波近似系數(shù)能夠清晰地區(qū)分出信號是否含有P300成分，所以選擇差別明顯的X1-X5和X9-X12這9個小波近似系數(shù)作為頻域特征向量。

2.3 空域特征提取

EEG能夠全方位輻射到頭皮上，但其在大腦不同區(qū)域的分布情況不同，所以可以對多導EEG空間分布特性進行分析。以空間信息作為特征，能夠達到對EEG有效成分進行空間定位的目的。獨立成分分析是目前最受關注的EEG空間分析方法，屬于盲源分離算法［11］（Blind Source Separation，BSS）的一種。它假設源信號之間相互獨立，且源信號和信號混合矩陣都是未知的，觀測信號是經(jīng)過線性混合后的源信號，這樣可以用一組在最大程度上逼近源信號的相互獨立的分量來表示觀測到的信號。

假設X為n維觀測信號矢量，S為獨立的m維未知源信號矢量，A為混合矩陣，W為混合矩陣A的逆，即A=W-1。ICA的目的就是尋找解混矩陣W，使觀測信號X通過解混矩陣W的投影變換得到輸出獨立分量U，使U為未知源信號S的最優(yōu)估計。為了將ICA具體化，需要用數(shù)學表達式表達：ICA的線性模型為X=AS，U=WX=WAS。ICA的原理如圖8所示。

圖8 ICA原理圖

由于事先選定了C1、C2、Cz、Pz四個通道的數(shù)據(jù)進行處理，因此得出式（6）用于求解混合矩陣W-1：

混合矩陣W-1中的矩陣元素Wij反映了第j個獨立分量Uj到觀測信號X的投影，而X是在大腦上不同電極測得，即混合矩陣中的元素可以表征EEG的空間分布特性，所以使用混合矩陣中的4×4=16個元素作為空域特征向量，以W1-W16命名。

通過時域、頻域及空域的特征提取進一步降低了數(shù)據(jù)處理量，得到了Cz導聯(lián)上的3個時域特征向量和9個頻域特征向量；Pz導聯(lián)上的2個時域特征向量和10個頻域特征向量；C1導聯(lián)上的3個時域特征向量和9個頻域特征向量；C2導聯(lián)上的3個時域特征向量和9個頻域特征向量以及所選4導上的16個空域特征向量組成的多域特征集合。

3 分類結果與討論

特征分類為腦電信號處理過程中的最后一步。在BCI系統(tǒng)中，對EEG進行處理的最終目的是得到能夠反應大腦真實意圖的控制信號，使大腦準確控制外部設備運行成為現(xiàn)實。評價EEG處理效果的最終指標是特征分類準確率，而特征提取的有效性直接影響著分類結果的準確性，即分類結果是特征提取效果的一個反映。

3.1 支持向量機

支持向量機（Support Vector Machine，SVM）是一種有監(jiān)督學習模型，在解決高維及非線性模式識別中表現(xiàn)出許多特有的優(yōu)勢，能夠有效地避免過學習、維數(shù)災難、局部極小等傳統(tǒng)分類方法存在的問題，在小樣本條件下仍然具有良好的泛化能力，因此被廣泛應用到EEG分類當中［12］。SVM的基本原理是：通過核函數(shù)對輸入數(shù)據(jù)進行非線性變換，將其映射到一個高維特征空間中使這些數(shù)據(jù)線性可分，然后在這個高維空間中求取最優(yōu)線性分類面，并使分類間隔最大，將不同的特征盡可能正確地分開。

核函數(shù)的選擇是SVM分類的關鍵。核函數(shù)可以將難以劃分的低維空間向量集映射到高維空間中并確保不會增加這個映射過程的計算復雜度。由于徑向基RBF核函數(shù)與其他核函數(shù)相比具有參數(shù)較少、復雜度低，更容易進行數(shù)值運算的優(yōu)點［13］，所以本文使用徑向基RBF核函數(shù)，其表達式為：

3.2 基于SVM的P300分類

在實驗過程中，A和B兩個受試者均做了5組實驗，每個受試者識別185個字符，其中85個字符帶有對應于各次閃爍的身份信息，即可以直接看出該段數(shù)據(jù)對應的字符，而另外100個字符不帶有上述信息，即無法判斷出應用算法識別出的字符是否和給予受試者刺激的字符是否一致。因此，本文選擇受試者A的85個字符數(shù)據(jù)作為訓練集，受試者B的85個字符數(shù)據(jù)作為測試集。SVM對不同特征進行分類的過程如圖9所示。

圖9 分類流程圖

3.2.1 多域特征提取的分類結果

表1為在不同疊加平均次數(shù)下本文的分類準確率和其它文獻分類準確率的對比結果，可以看出對P300進行多域特征提取，能夠明顯提高SVM的分類準確率。在本文中進行5次疊加平均后的分類準確率高于對照組［14-15］進行10次疊加平均的分類準確率。

表1 不同疊加平均次數(shù)下的P300分類準確率/%

3.2.2 單域特征提取的分類結果

表2為僅使用時域能量熵作為特征的P300分類準確率，表3為僅使用小波近似系數(shù)作為特征的P300分類準確率，表4為僅使用混合矩陣元素作為特征的的P300分類準確率。

表2 以時域能量熵為特征的P300分類準確率/%

表3 以小波近似系數(shù)為特征的P300分類準確率/%

表4 以混合矩陣中元素為特征的P300分類準確率/%

從上述3個表中可以得出以下結論：

（1）P300的分類準確率隨著疊加次數(shù)的增加而增加，說明疊加平均次數(shù)越高越能夠?qū)300的特征顯示出來。

（2）對于單一特征提取方法，小波分析提取近似系數(shù)作為特征的方法得到的分類準確率高于以時域能量熵為特征的分類準確率，而使用ICA的空域特征提取方法相比于其他兩種方法得到的分類準確率最低。

4 結論

針對P300腦電信號微弱及單一的特征提取方法無法全面表征P300特性的問題，本文提出一種多域特征提取方法。該方法以信號在時域上的能量信息作為時域特征，以小波分解后得到的近似系數(shù)作為頻域特征，以ICA分解得到的混合矩陣中的元素作為空域特征。從多角度提取特征，使得P300的本質(zhì)特征得到近乎全面的表達，不但降低了數(shù)據(jù)處理量，而且為SVM分類的結果提供了保證。

實驗結果表明，采用本文提出的多域特征提取方法得到的分類準確率明顯高于單一特征提取的分類準確率。

本文在數(shù)據(jù)處理過程中，首先對有用電極進行選擇，剔除無關電極，使待處理的數(shù)據(jù)量大幅減少，根據(jù)相關研究選擇P300相對明顯的電極C1，C2，Cz，Pz進行處理，使數(shù)據(jù)處理量為原始數(shù)據(jù)的1/16，大大節(jié)省了數(shù)據(jù)處理的時間。然后在特征提取之前進行數(shù)據(jù)截取及相同刺激下的數(shù)據(jù)疊加平均操作，提高了信噪比，使P300的特征更加明顯。