陸振宇,陸旭峰,楊瑞洪,常珊
(1.江蘇理工學院電氣信息工程學院生物信息與醫藥工程研究所,常州213001;2.揚州工業職業技術學院,揚州225127)
腦機接口(Brain Computer Interface,BCI)是綜合神經科學、信號處理、機器學習和信息技術等領域的跨學科技術。BCI 的目的是把大腦活動轉變成對設備的控制指令,具有恢復失去的感官和運動功能的潛力,如治療帕金森病的腦深度電刺激,截肢者和脊髓損傷病人的假肢、輪椅,還有中風和肌肉萎縮患者使用的鼠標和單詞拼寫設備[1]。BCI 已經在游戲娛樂、機器人、生物識別和教育等應用領域表現出強大的發展前景。
目前,運用較為廣泛的信號是當受試者自主想象某種特定動作所產生的神經活動,稱之為基于運動想象(Motor Imagery,MI)的腦電信號。研究表明,當受試者進行手、腳、舌頭動作的運動想象,大腦皮層相對應區域的μ頻帶(8-13Hz)表現顯著,而對側區域其節律能量下降,該現象被稱為事件相關同步/去同步(ERD/ERS)現象,這也為信號的特征提取提供了依據[2-4]。
BCI 的實現主要涉及信號處理和模式識別。信號處理包括數據的預處理,如濾波、偽跡去除,還有關鍵的特征提取技術。文獻中的特征提取方法一般有,頻域分析有自回歸(AR)模型、功率譜估計;時頻分析有小波變換、小波包變換、局部均值分解;空間域有主成分分析(PCA)、獨立分量分析(ICA)和共空間模式(CSP),用于多導聯的腦電信號[5-11]。模式識別是利用機器學習技術,對輸入特征進行分類。常用的二分類技術有線性判別分析(LDA)、神經網絡(NN)、支持向量機(SVM)、基于隨機森林和Boosting 的集成技術。
通過對相關文獻的研究,信號變換和分解的方法多樣,一般在信號處理后選取能量、熵等作為特征,但變換方法及提取的特征單一化,無法多方面地表示腦電信號。本文選取提取算法易于理解、計算量小、效果良好的AR 模型、小波包變換、CSP 方法對信號進行處理,提取信號的多方面特征,采用SVM、LDA 分類器進行模式識別,并對分類結果評價。
本文實驗的數據是來源于2003 年BCI 競賽的Graz 數據集,由奧地利Graz 科技大學醫學系、生物醫學工程研究所提供[12]。實驗范式是受試者放松地坐在手扶椅上,根據屏幕顯示做出相應的反饋。前2s 為準備時間,在第2s 時,蜂鳴器提示實驗開始,同時屏幕出現1s 的“+”符號,3-9s 是受試者按給出的箭頭方向進行對應的運動想象。實驗采集了C3、CZ、C4 通道的腦電信號,采樣頻率為128Hz,并通過0.5-30Hz 的帶通濾波。整個數據集共有280 個樣本,其中左右手運動想象各140,訓練和測試樣本各占一半,數據集分布均勻。
AR 模型可用如下公式表示:

式中,ε(n)是均值為零、方差為σ2的白噪聲,p 為模型階次,ap(i)為自回歸模型參數[4-5]。本文選取競賽數據C3、C4 通道,時間段在4-7s 的信號,基于最終預測誤差準則,確定階次為5,使用Burg 算法提取其AR 模型系數和系數差作為特征。
小波包變換是為了克服小波分解在高頻段的頻率分辨率較差,而在低頻段表現較好的基礎上提出的。它是一種更精細的信號分析的方法,提高了信號的時域分辨率[13]。
競賽數據采樣頻率為128Hz,依據Nyquist 采樣定理,腦電信號頻率范圍為1~64Hz,使用小波基函數對C3、C4 信號進行4 層小波包分解,分解頻寬為4Hz。分解后的第4 層小波包節點頻率范圍如表1 所示,具有代表性的μ節律信號為節點(4,3)。圖2 為采用bior3.3 小波基,對單個樣本的C3、C4 通道信號做小波包變換得到的時域和頻域圖。

表1 第4 層小波包分解的節點及頻帶范圍

圖2 單樣本節點(4,3)的bior3.3小波包變換
腦電信號經小波包分解后,節點(i,j)信號的小波包能量由計算該節點的小波包系數平方和得到。

式中,i 為分解層數,j 為該層下的節點,di,j(n)為節點(i,j)的小波包系數。
根據能量守恒定律,分解后的各節點能量總和等于原信號的能量,則相對小波包能量為節點信號能量與該層節點總能量的比值:

式中,Eraw為原信號能量,Etotal為原信號分解i 層后各節點能量總和,pi,j則為節點(i,j)能量的歸一化,稱為相對能量。
信息熵可以表征信號的不確定度,結合小波包變換,可得出分解后的節點信號小波包熵為:

本文選取C3、C4 通道信號做小波包變換,以一次實驗所記錄的C3 信號為例,具體步驟如下:
(1)采用固定時間窗的方法,以128 點的時間窗,將長度為1152 點的信號分割成9 段;
(2)在一段128 點的信號上做小波包分解,小波基選擇bior3.3[4],分解層數為4,選擇節點(4,3)的重構信號,由式2 計算出小波包能量,式4 得出小波包相對能量,式5 得出小波包熵。
(3)C4 通道重復1、2 步驟,一個樣本的小波包變換結果為2×9 的矩陣,按樣本類別進行疊加平均,得到C3、C4 通道腦電信號的小波包能量、相對能量和熵的變化,圖3 為右手運動想象下的平均小波包能量和熵走勢圖。

圖3 右手運動想象的bior3.3平均小波包變換特征
由圖3 可知,4-7s 的C3、C4 通道小波包相對能量和小波包熵區分度明顯,走勢基本一致,數值大小不一樣。本文將其二維向量特征轉化成一維向量,一個樣本C3、C4 通道的小波包能量差、相對能量差、熵差△WPE 作為分類的特征,圖4 為小波包熵差的走勢圖。


圖4 左右手運動想象的bior3.3平均小波包熵差
CSP 是尋找空間濾波器,使得濾波處理后的數據與其中一類方差達到最大,而與另一類方差達到最小的算法。基于運動想象的BCI,依賴特定的頻帶能量,對于兩類的左右手運動想象信號,經過CSP 處理得到的特征向量增強了兩類的區別,可區分程度最大化[14]。
單次腦電信號樣本可由Ns×Nc的二維矩陣E 表示,其中,Ns是采樣點數,Nc是通道數,則歸一化協方差矩陣為:

式中ET表示E 的轉置,表示矩陣的秩。按樣本類別求疊加平均的協方差矩陣,然后總樣本的協方差矩陣為:

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

將U、λ進行降序排列,求得白化矩陣:

對類別1 的協方差矩陣使用白化矩陣進行轉變:

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

同樣對U1、λ1進行降序排列,得到類別1 的投影矩陣W:

取W 的前m 行、后m 行,構成空間濾波器W’,原信號E 經過空間濾波器得到:

對E’按行求方差,并取對數,得到特征值f:

一般地還有,在方差的歸一化后取對數,得到特征值f:

本文選取4-7s 的運動想象腦電信號,對其做CSP變換,按式(15)、(16)分別計算,各自得到2 個特征值,用于區分左右手運動想象。

圖5 左右運動想象各70個樣本的特征值散點圖
線性判別分析(LDA)是線性二分類器,將輸入向量映射到一個超平面,該超平面將輸入空間劃分為兩個半空間:每個空間代表一種類別。決策邊界由超平面公式決定:

式中w 為法向量,w0為閾值,均由訓練數據得出。對于新的輸入特征向量,代入式16,如果g(x)≥0,為類別1,反之為類別2。
支持向量機(SVM)是尋找兩類樣本間隔最大的分離超平面,以獲得最好的泛化能力的一種分類器。本文采用libsvm 工具,選擇C-SVM 類型,利用徑向基核函數實現腦電信號特征的非線性映射,使用網格搜索獲得最優參數C、gama。
為了驗證本文提出方法的有效性,將單特征和多特征融合提取方法運用在整個數據集上,采用SVM、LDA 分類器對其分類,結果如表2、3,其中識別率是10折交叉驗證的平均值。
通過表2 可以直觀地發現:單特征選取方面,AR系數和系數差都能較好地被區分,識別率最高為87.5%;小波包熵差的效果高于能量差、相對能量差,原因在于小波包熵的數量級比后兩者大,從圖3 可知,小波包能量的范圍在0.25~0.75,而小波包熵是1.2~3.0;而共空間模式下的兩種特征在分類器上表現了較大差異,未歸一化的特征在LDA 上的表現最好,原因在于數據的數值較大,而對于SVM,無論特征是否歸一化,準確率都是73.93%,表現一般。
通過表2 的分析,選取分類效果好的AR 模型系數、小波包熵差、未歸一化的CSP 特征,進行排列組合,用于多特征融合,得到表3 的結果。其中,分類器方面,LDA 的識別率均比SVM 略高;特征融合方面,三者組合的準確率最高,達到了91.43%,而兩兩組合的準確率比各自的單特征高,均在90%左右。

表2 單特征的平均識別率(單位:%)

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