陸振宇，陸旭峰，楊瑞洪，常珊

（1.江蘇理工學(xué)院電氣信息工程學(xué)院生物信息與醫(yī)藥工程研究所，常州213001；2.揚(yáng)州工業(yè)職業(yè)技術(shù)學(xué)院，揚(yáng)州225127）

0 引言

腦機(jī)接口（Brain Computer Interface，BCI）是綜合神經(jīng)科學(xué)、信號(hào)處理、機(jī)器學(xué)習(xí)和信息技術(shù)等領(lǐng)域的跨學(xué)科技術(shù)。BCI 的目的是把大腦活動(dòng)轉(zhuǎn)變成對設(shè)備的控制指令，具有恢復(fù)失去的感官和運(yùn)動(dòng)功能的潛力，如治療帕金森病的腦深度電刺激，截肢者和脊髓損傷病人的假肢、輪椅，還有中風(fēng)和肌肉萎縮患者使用的鼠標(biāo)和單詞拼寫設(shè)備[1]。BCI 已經(jīng)在游戲娛樂、機(jī)器人、生物識(shí)別和教育等應(yīng)用領(lǐng)域表現(xiàn)出強(qiáng)大的發(fā)展前景。

目前，運(yùn)用較為廣泛的信號(hào)是當(dāng)受試者自主想象某種特定動(dòng)作所產(chǎn)生的神經(jīng)活動(dòng)，稱之為基于運(yùn)動(dòng)想象（Motor Imagery，MI）的腦電信號(hào)。研究表明，當(dāng)受試者進(jìn)行手、腳、舌頭動(dòng)作的運(yùn)動(dòng)想象，大腦皮層相對應(yīng)區(qū)域的μ頻帶（8-13Hz）表現(xiàn)顯著，而對側(cè)區(qū)域其節(jié)律能量下降，該現(xiàn)象被稱為事件相關(guān)同步/去同步（ERD/ERS）現(xiàn)象，這也為信號(hào)的特征提取提供了依據(jù)[2-4]。

BCI 的實(shí)現(xiàn)主要涉及信號(hào)處理和模式識(shí)別。信號(hào)處理包括數(shù)據(jù)的預(yù)處理，如濾波、偽跡去除，還有關(guān)鍵的特征提取技術(shù)。文獻(xiàn)中的特征提取方法一般有，頻域分析有自回歸（AR）模型、功率譜估計(jì)；時(shí)頻分析有小波變換、小波包變換、局部均值分解；空間域有主成分分析（PCA）、獨(dú)立分量分析（ICA）和共空間模式（CSP），用于多導(dǎo)聯(lián)的腦電信號(hào)[5-11]。模式識(shí)別是利用機(jī)器學(xué)習(xí)技術(shù)，對輸入特征進(jìn)行分類。常用的二分類技術(shù)有線性判別分析（LDA）、神經(jīng)網(wǎng)絡(luò)（NN）、支持向量機(jī)（SVM）、基于隨機(jī)森林和Boosting 的集成技術(shù)。

通過對相關(guān)文獻(xiàn)的研究，信號(hào)變換和分解的方法多樣，一般在信號(hào)處理后選取能量、熵等作為特征，但變換方法及提取的特征單一化，無法多方面地表示腦電信號(hào)。本文選取提取算法易于理解、計(jì)算量小、效果良好的AR 模型、小波包變換、CSP 方法對信號(hào)進(jìn)行處理，提取信號(hào)的多方面特征，采用SVM、LDA 分類器進(jìn)行模式識(shí)別，并對分類結(jié)果評價(jià)。

1 實(shí)驗(yàn)數(shù)據(jù)

本文實(shí)驗(yàn)的數(shù)據(jù)是來源于2003 年BCI 競賽的Graz 數(shù)據(jù)集，由奧地利Graz 科技大學(xué)醫(yī)學(xué)系、生物醫(yī)學(xué)工程研究所提供[12]。實(shí)驗(yàn)范式是受試者放松地坐在手扶椅上，根據(jù)屏幕顯示做出相應(yīng)的反饋。前2s 為準(zhǔn)備時(shí)間，在第2s 時(shí)，蜂鳴器提示實(shí)驗(yàn)開始，同時(shí)屏幕出現(xiàn)1s 的“+”符號(hào)，3-9s 是受試者按給出的箭頭方向進(jìn)行對應(yīng)的運(yùn)動(dòng)想象。實(shí)驗(yàn)采集了C3、CZ、C4 通道的腦電信號(hào)，采樣頻率為128Hz，并通過0.5-30Hz 的帶通濾波。整個(gè)數(shù)據(jù)集共有280 個(gè)樣本，其中左右手運(yùn)動(dòng)想象各140，訓(xùn)練和測試樣本各占一半，數(shù)據(jù)集分布均勻。

2 特征提取算法

2.1 AR模型

AR 模型可用如下公式表示：

式中，ε（n）是均值為零、方差為σ2的白噪聲，p 為模型階次，ap（i）為自回歸模型參數(shù)[4-5]。本文選取競賽數(shù)據(jù)C3、C4 通道，時(shí)間段在4-7s 的信號(hào)，基于最終預(yù)測誤差準(zhǔn)則，確定階次為5，使用Burg 算法提取其AR 模型系數(shù)和系數(shù)差作為特征。

2.2 小波包變換

小波包變換是為了克服小波分解在高頻段的頻率分辨率較差，而在低頻段表現(xiàn)較好的基礎(chǔ)上提出的。它是一種更精細(xì)的信號(hào)分析的方法，提高了信號(hào)的時(shí)域分辨率[13]。

競賽數(shù)據(jù)采樣頻率為128Hz，依據(jù)Nyquist 采樣定理，腦電信號(hào)頻率范圍為1～64Hz，使用小波基函數(shù)對C3、C4 信號(hào)進(jìn)行4 層小波包分解，分解頻寬為4Hz。分解后的第4 層小波包節(jié)點(diǎn)頻率范圍如表1 所示，具有代表性的μ節(jié)律信號(hào)為節(jié)點(diǎn)（4,3）。圖2 為采用bior3.3 小波基，對單個(gè)樣本的C3、C4 通道信號(hào)做小波包變換得到的時(shí)域和頻域圖。

表1 第4 層小波包分解的節(jié)點(diǎn)及頻帶范圍

圖2 單樣本節(jié)點(diǎn)（4,3）的bior3.3小波包變換

腦電信號(hào)經(jīng)小波包分解后，節(jié)點(diǎn)（i,j）信號(hào)的小波包能量由計(jì)算該節(jié)點(diǎn)的小波包系數(shù)平方和得到。

式中，i 為分解層數(shù)，j 為該層下的節(jié)點(diǎn)，di,j（n）為節(jié)點(diǎn)（i,j）的小波包系數(shù)。

根據(jù)能量守恒定律，分解后的各節(jié)點(diǎn)能量總和等于原信號(hào)的能量，則相對小波包能量為節(jié)點(diǎn)信號(hào)能量與該層節(jié)點(diǎn)總能量的比值：

式中，Eraw為原信號(hào)能量，Etotal為原信號(hào)分解i 層后各節(jié)點(diǎn)能量總和，pi,j則為節(jié)點(diǎn)（i,j）能量的歸一化，稱為相對能量。

信息熵可以表征信號(hào)的不確定度，結(jié)合小波包變換，可得出分解后的節(jié)點(diǎn)信號(hào)小波包熵為：

本文選取C3、C4 通道信號(hào)做小波包變換，以一次實(shí)驗(yàn)所記錄的C3 信號(hào)為例，具體步驟如下：

（1）采用固定時(shí)間窗的方法，以128 點(diǎn)的時(shí)間窗，將長度為1152 點(diǎn)的信號(hào)分割成9 段；

（2）在一段128 點(diǎn)的信號(hào)上做小波包分解，小波基選擇bior3.3[4]，分解層數(shù)為4，選擇節(jié)點(diǎn)（4,3）的重構(gòu)信號(hào)，由式2 計(jì)算出小波包能量，式4 得出小波包相對能量，式5 得出小波包熵。

（3）C4 通道重復(fù)1、2 步驟，一個(gè)樣本的小波包變換結(jié)果為2×9 的矩陣，按樣本類別進(jìn)行疊加平均，得到C3、C4 通道腦電信號(hào)的小波包能量、相對能量和熵的變化，圖3 為右手運(yùn)動(dòng)想象下的平均小波包能量和熵走勢圖。

圖3 右手運(yùn)動(dòng)想象的bior3.3平均小波包變換特征

由圖3 可知，4-7s 的C3、C4 通道小波包相對能量和小波包熵區(qū)分度明顯，走勢基本一致，數(shù)值大小不一樣。本文將其二維向量特征轉(zhuǎn)化成一維向量，一個(gè)樣本C3、C4 通道的小波包能量差、相對能量差、熵差△WPE 作為分類的特征，圖4 為小波包熵差的走勢圖。

圖4 左右手運(yùn)動(dòng)想象的bior3.3平均小波包熵差

2.3 共空間模式

CSP 是尋找空間濾波器，使得濾波處理后的數(shù)據(jù)與其中一類方差達(dá)到最大，而與另一類方差達(dá)到最小的算法?；谶\(yùn)動(dòng)想象的BCI，依賴特定的頻帶能量，對于兩類的左右手運(yùn)動(dòng)想象信號(hào)，經(jīng)過CSP 處理得到的特征向量增強(qiáng)了兩類的區(qū)別，可區(qū)分程度最大化[14]。

單次腦電信號(hào)樣本可由Ns×Nc的二維矩陣E 表示，其中，Ns是采樣點(diǎn)數(shù)，Nc是通道數(shù)，則歸一化協(xié)方差矩陣為：

式中ET表示E 的轉(zhuǎn)置，表示矩陣的秩。按樣本類別求疊加平均的協(xié)方差矩陣，然后總樣本的協(xié)方差矩陣為：

式中Cov1、Cov2 分別表示類別1、2 的協(xié)方差矩陣。對總樣本協(xié)方差矩陣作特征值分解，得到特征向量矩陣U、特征值對角陣λ：

將U、λ進(jìn)行降序排列，求得白化矩陣：

對類別1 的協(xié)方差矩陣使用白化矩陣進(jìn)行轉(zhuǎn)變：

對其作特征值分解，得到特征向量矩陣U1、特征值對角陣λ1：

同樣對U1、λ1進(jìn)行降序排列，得到類別1 的投影矩陣W：

取W 的前m 行、后m 行，構(gòu)成空間濾波器W’，原信號(hào)E 經(jīng)過空間濾波器得到：

對E’按行求方差，并取對數(shù)，得到特征值f：

一般地還有，在方差的歸一化后取對數(shù)，得到特征值f：

本文選取4-7s 的運(yùn)動(dòng)想象腦電信號(hào)，對其做CSP變換，按式（15）、（16）分別計(jì)算，各自得到2 個(gè)特征值，用于區(qū)分左右手運(yùn)動(dòng)想象。

圖5 左右運(yùn)動(dòng)想象各70個(gè)樣本的特征值散點(diǎn)圖

3 模式判別

3.1 線性判別分析

線性判別分析（LDA）是線性二分類器，將輸入向量映射到一個(gè)超平面，該超平面將輸入空間劃分為兩個(gè)半空間：每個(gè)空間代表一種類別。決策邊界由超平面公式?jīng)Q定：

式中w 為法向量，w0為閾值，均由訓(xùn)練數(shù)據(jù)得出。對于新的輸入特征向量，代入式16，如果g（x）≥0，為類別1，反之為類別2。

3.2 支持向量機(jī)

支持向量機(jī)（SVM）是尋找兩類樣本間隔最大的分離超平面，以獲得最好的泛化能力的一種分類器。本文采用libsvm 工具，選擇C-SVM 類型，利用徑向基核函數(shù)實(shí)現(xiàn)腦電信號(hào)特征的非線性映射，使用網(wǎng)格搜索獲得最優(yōu)參數(shù)C、gama。

4 實(shí)驗(yàn)結(jié)果分析

為了驗(yàn)證本文提出方法的有效性，將單特征和多特征融合提取方法運(yùn)用在整個(gè)數(shù)據(jù)集上，采用SVM、LDA 分類器對其分類，結(jié)果如表2、3，其中識(shí)別率是10折交叉驗(yàn)證的平均值。

通過表2 可以直觀地發(fā)現(xiàn)：單特征選取方面，AR系數(shù)和系數(shù)差都能較好地被區(qū)分，識(shí)別率最高為87.5%；小波包熵差的效果高于能量差、相對能量差，原因在于小波包熵的數(shù)量級(jí)比后兩者大，從圖3 可知，小波包能量的范圍在0.25～0.75，而小波包熵是1.2～3.0；而共空間模式下的兩種特征在分類器上表現(xiàn)了較大差異，未歸一化的特征在LDA 上的表現(xiàn)最好，原因在于數(shù)據(jù)的數(shù)值較大，而對于SVM，無論特征是否歸一化，準(zhǔn)確率都是73.93%，表現(xiàn)一般。

通過表2 的分析，選取分類效果好的AR 模型系數(shù)、小波包熵差、未歸一化的CSP 特征，進(jìn)行排列組合，用于多特征融合，得到表3 的結(jié)果。其中，分類器方面，LDA 的識(shí)別率均比SVM 略高；特征融合方面，三者組合的準(zhǔn)確率最高，達(dá)到了91.43%，而兩兩組合的準(zhǔn)確率比各自的單特征高，均在90%左右。

表2 單特征的平均識(shí)別率（單位：%）

表3 多特征融合的分類準(zhǔn)確率（單位：%）

5 結(jié)語

本文提出了多方面對運(yùn)動(dòng)想象腦電信號(hào)進(jìn)行特征提取及融合，并使用分類器對其評價(jià)。AR 模型代表了腦電信號(hào)的功率譜估計(jì)，僅使用C3、C4 兩通道信號(hào)的系數(shù)，也可以得到較好的分類結(jié)果；小波包變換是定位到給定頻帶，計(jì)算熵值可以反映其能量的差異性，選擇C3、C4 通道信號(hào)的熵差作為特征，也得到了較好的效果；CSP 算法是提取腦電信號(hào)的空間特征，適用于較多導(dǎo)聯(lián)，且需要訓(xùn)練數(shù)據(jù)構(gòu)造濾波器，本文的競賽數(shù)據(jù)是3 導(dǎo)聯(lián)，使用CSP 僅得到兩個(gè)特征，分類效果不穩(wěn)定。而多特征融合的方式，則彌補(bǔ)了單特征的缺陷，在分類效果上也更好，這為運(yùn)動(dòng)想象腦電信號(hào)的特征工程提供了思路。在分類器的選擇上，SVM 雖然泛化能力高，但計(jì)算量大，對數(shù)據(jù)敏感，而LDA 實(shí)現(xiàn)簡單，運(yùn)算速度快，效果也較好，適用于在線BCI。