腦機(jī)接口的廣義核線性判別分析方法研究

王金甲 胡 備

（燕山大學(xué)信息科學(xué)與工程學(xué)院，秦皇島 0 66004）

引言

腦機(jī)接口（brain-computer interface，BCI）主要目的就是希望建立人腦和計算機(jī)（或儀器設(shè)備）的直接通信，進(jìn)而實現(xiàn)人腦對外部設(shè)備的控制［1］。從復(fù)雜腦電EEG信號中提取可以有效分類的特征就成為了腦機(jī)接口中一個嚴(yán)峻挑戰(zhàn)。

對于腦機(jī)接口左右手想象運(yùn)動的腦電EEG信號特征選擇和提取，很多文獻(xiàn)中都有研究報道。研究者采取很多腦電特征對想象運(yùn)動進(jìn)行分類，如共空間模式（common spatial pattern，CSP）、核共空間模式（kernel CSP，KCSP）、核線性判別分析（kernel discriminant analysis，KDA）、自回歸模型參數(shù)等［2－5］。CSP方法針對腦機(jī)接口中想象運(yùn)動的分類效果已被廣泛認(rèn)可［2］，但其缺點在于對兩類樣本嚴(yán)格的線性假設(shè)，針對這一問題文獻(xiàn)［3］提出了KCSP，引入核函數(shù)很好解決了線性要求限定，新問題是計算復(fù)雜性；文獻(xiàn)［6］采用聚類方法解決計算復(fù)雜性問題，這在一定程度上丟失了部分特征。對此本研究提出了一種基于核方法和廣義奇異值分解［7］（generalized singular value decomposition，GSVD）的廣義核線性判別分析方法（GKLDA）。首先通過一個非線性映射，將訓(xùn)練數(shù)據(jù)映射到一個高維特征空間，解決了CSP中對原始數(shù)據(jù)線性的要求，使得不同類別EEG信號在該空間線性可分（或近似線性可分），采用線性判別分析，尋找一個最佳投影方向，即為最優(yōu)空域濾波器。在解該濾波器時，為了最大限度的釋放不同類間的非線性模式，同時解決小樣本采樣問題，采用GSVD方法解除了奇異性限定。最后將數(shù)據(jù)通過該最優(yōu)空域濾波器，提取最佳非線性特征進(jìn)行分類。分類器采用fisher線性判別分析分類器。文獻(xiàn)［8］提出了一種廣義判別分析（generalized discriminant analysis，GDA），本研究提出將該方法應(yīng)用在腦電分類中。在實驗部分，采用了線性、多項式及高斯等三種不同的核函數(shù)，將GKLDA與KCSP、KDA、GDA及普通 LDA分別用于相同的腦機(jī)接口公開數(shù)據(jù)進(jìn)行特征提取對比分類結(jié)果，結(jié)果表明 GKLDA能取得很好的效果，值得推廣。

1 CSP方法及其核推廣KCSP

CSP方法所基于的事實是，假定大腦中每一類運(yùn)動都是由一種或幾種源所觸發(fā)的。假設(shè)有兩類運(yùn)動記錄下來的EEG信號為X1和X2，對其各自空間協(xié)方差矩陣R1，R2取均值，得到平均協(xié)方差矩陣和。對復(fù)合空間協(xié)方差矩陣 R =+做特征值分解，U和λ分別為其特征向量和特征值。定義白化矩陣 P，P=λ－1/2UT。分別對 PPT和PPT做特征值分解，可以證明兩者有共同的特征向量，且特征值對角陣之和為單位陣。兩者最大m個特征值所對應(yīng)的特征向量記為V1m和V2m，即可得空域濾波器 F，F(xiàn)=［V1mV2m］P。

KCSP方法中，首先通過一個非線性映射將EEG數(shù)據(jù)映射到高維空間，即X→Φ。類比于CSP，在這個新空間中有R1=Φ1W1Φ1，從而R=R1+R2=ΦtWΦtT。W是所有樣本數(shù)n的倒數(shù)組成的對角陣，其目的在于求協(xié)方差矩陣的平均值。為了在高維空間完成CSP方法，首先要解決R的特征值分解問題。但是直接分解是不可能的，因為通常情況下R維數(shù)特別高。文獻(xiàn)［6］中給出一個思路，即通過構(gòu)造特征值分解因式，用已知或可求變量代替未知或不易求變量，通過構(gòu)造核矩陣的形式代替高維空間的運(yùn)算。最終濾波器表達(dá)式為：

式中，V是KW特征向量，K=ΦtTΦt是核矩陣；P=Λ－1Ψ－1/2VTWΦtT為白化矩陣，Λ是KW特征值；Ψ是一個對角矩陣，它的第j列元素為vjTWvj，vj是V的第j列。細(xì)節(jié)內(nèi)容詳見文獻(xiàn)［3，6］。

2 LDA方法及其核推廣KDA

以兩類數(shù)據(jù)為例。設(shè)有一集合X，包含n個d維樣本x1，x2，…xn，其中n1個屬于w1類，n2個屬于w2類，n1+n2=n，LDA目標(biāo)是找到線性投影方向，使得投影后類內(nèi)散度最小，類間散度最大。Sw為總類內(nèi)離散度矩陣，Sb為類間離散度矩陣，最佳投影方向w通過下式目標(biāo)函數(shù)最大化：

式（2）為廣義Rayleigh商形式，可以用Lagrange乘子法求解，轉(zhuǎn)化為方程Sbw*=λSww*求解問題。解該方程即為式（2）的極值。細(xì)節(jié)內(nèi)容見文獻(xiàn)［9］。

3 基于核方法和GSVD的廣義核線性判別分析

矩陣與其轉(zhuǎn)置的乘積一樣，離散度矩陣也分解成其平方根因子乘積的形式。

此時，在映射空間中，目標(biāo)函數(shù)可寫為如下Rayleigh商形式

式（6）可轉(zhuǎn)化為求解下式廣義特征值問題

據(jù)再生核理論，可知式（7）的解空間可限定在由{φ（x）}張成的空間中。令

式中 Y =［φ（x1）…φ（xn）］，α=［α1…αn］T。將φ分別代入式（7），左乘 YT，得

此時，GTY，HTY可以通過下式計算：

本研究提出一個廣義核線性判別分析方法。首先，對于EEG數(shù)據(jù)X，大小為M×T×n，M為通道數(shù)，T為每次試驗的采樣點數(shù)，n為試驗次數(shù)。兩類試驗次數(shù)分別為n1、n2，即n1+n2=n，n為EEG信號所有試驗的次數(shù)。將數(shù)據(jù)重新排列為二維數(shù)據(jù)，每一列為一次實驗的所有采樣點，即（M×T）×n。假定通過一個非線性的映射函數(shù)將輸入空間的數(shù)據(jù)映射到一個希爾伯特空間φ：R→F即X→φ（X），據(jù)式（3），在映射后的空間中投影方向向量目標(biāo)函數(shù)變?yōu)?/p>

散度矩陣與協(xié)方差矩陣相似，都是半正定的。正像在許多信號處理問題中協(xié)方差矩陣分解為一個

式（9）等價于：

為了釋放KG，KH的非奇異限定，對，)采取廣義奇異值分解 G SVD，即UTα=［Γ0］，VTα = ［Γ0］，由此可得，

式（11）和式（12）中的α即為式（5）的解，也是所求最優(yōu)空域濾波器。

對于新輸入EEG信號Z，可獲得該數(shù)據(jù)特征，

式中，j=1，2，…m，m為α的列數(shù)，一般情況下等于n，但適當(dāng)?shù)倪x取m有利于正確率的提高，筆者對m的選擇做了進(jìn)一步實驗。式（13）中濾波后的EEG信號其能量可作為特征向量用于分類。

4 實驗方法

實驗分為兩個部分，實驗一采用BCI競賽一、二和三的部分?jǐn)?shù)據(jù)分別進(jìn)行了實驗，實驗二采用自己的數(shù)據(jù)做了實驗。分別采用了GKLDA、KCSP、KDA、GDA等4種方法進(jìn)行對比；并且給出相應(yīng)數(shù)據(jù)已知的最優(yōu)正確率作為對比；此外采用了線性核函數(shù)，多項式核函數(shù)和高斯核函數(shù)等作了對比。其中多項式核函數(shù)k（x，y）=（xTy）d，d=1.5，d的選擇是經(jīng)過大量實驗后的最優(yōu)結(jié)果。高斯核函數(shù)k（x，y）=exp（－x－y2/2σ2），σ由十倍交叉驗證獲得。實驗中的分類投影圖為GKLDA方法提取特征后的一維投影，其中橫軸為樣本數(shù)，縱軸為提取的特征降到一維后的值。

實驗步驟如下：首先是數(shù)據(jù)預(yù)處理，即將三維數(shù)據(jù)（通道×采樣點×試驗次數(shù)）變?yōu)閮删S（［通道×采樣點］×試驗次數(shù)），符合GKLDA的處理格式；其次對數(shù)據(jù)實施前文所述的GKLDA以及需要比較的幾種方法，提取各類別數(shù)據(jù)的特征；最后將各類的特征送入分類器，得出該數(shù)據(jù)集的正確率。

5 實驗結(jié)果

5.1 實驗一

1.數(shù)據(jù)采用BCI競賽一的數(shù)據(jù)集，本數(shù)據(jù)集來自一名受試者，其任務(wù)是根據(jù)提示用手指敲擊指定鍵盤，共有413個訓(xùn)練樣本，其中兩類樣本數(shù)分別為194和219，100個測試樣本。27導(dǎo)信號，采樣率為100Hz。GKLDA降維后的投影見圖1，其中橫軸為試驗次數(shù)n，縱軸為提取特征后將特征維數(shù)降為一維時的投影幅值，不同類別分別用星號和圓點表示，而最后分類并沒有降維。

4種方法下的訓(xùn)練和測試結(jié)果見表1。高斯核中KCSP的參數(shù)2σ2=1×1010，KDA的參數(shù)2σ2=4×106，GDA的參數(shù)2σ2=1×106，GKLDA的參數(shù)2σ2=1×106。GKLDA正確率較高。本數(shù)據(jù)集用普通LDA特征和CSP提取特征的分類錯誤率結(jié)果分別為為11.0％和23.0％。競賽第一名的錯誤率為4.0％。

圖1 數(shù)據(jù)一的投影結(jié)果。（a）Linear核訓(xùn)練數(shù)據(jù)；（b）Linear核測試數(shù)據(jù)；（c）Polynomial核訓(xùn)練數(shù)據(jù)；（d）Polynomial核測試數(shù)據(jù)；（e）Gaussian核訓(xùn)練數(shù)據(jù)；（f）Gaussian核測試數(shù)據(jù)Fig.1 The projection of data I.（a）Training data with linear kernel；（b）Testing data with linear kernel；（c）Training data with polynomial kernel；（d）Testing data with polynomial kernel；（e）Training data with Gaussian kernel；（f）Testing data with Gaussian kernel

表1 數(shù)據(jù)一在4種方法下的分類錯誤率（％）Tab.1 The error rate（％）on datasetⅡ with four different methods

2.數(shù)據(jù)采用BCI競賽二的數(shù)據(jù)集Ⅳ，該數(shù)據(jù)集為左右手運(yùn)動，受試者坐在椅子上，手臂放松放在桌子上，其任務(wù)是敲擊指定按鍵。共有樣本416個，其中316個為訓(xùn)練集，左右手分別為159和157，其他100個作為測試集。電極數(shù)為28，采樣率100Hz，時間長度為0.5 s。圖2所示GKLDA降維后的投影結(jié)果，其中橫軸為試驗次數(shù)n，縱軸為提取特征后將特征維數(shù)降為一維時的投影幅值，不同類別分別用星號和圓點表示。而最后分類并沒有降維。

4種方法下的訓(xùn)練和測試結(jié)果見表2。高斯核中 K CSP的參數(shù) 2 σ2=6.4×107，KDA的參數(shù) 2 σ2=2.5×107，GDA的參數(shù) 2 σ2=6.4×107，GKLDA的參數(shù)2σ2=6.4×107。GKLDA正確率較高。本數(shù)據(jù)集用普通LDA特征提取直接分類，結(jié)果錯誤率為28.0％；文獻(xiàn)［11］中采用 C SSD和 F DA結(jié)合的方法，錯誤率為16.0％。

3.數(shù)據(jù)采用 BCI競賽三的數(shù)據(jù)集ⅢB中的S4b，該數(shù)據(jù)集為在線反饋的想象左右手運(yùn)動，反饋時間為第4 s到第7 s。共有樣本1 080個，其中540個為訓(xùn)練集，左右手各為270，其他540個作為測試集。電極數(shù)為2，采樣率125Hz，時間長度為3 s。圖3所示為GKLDA降維后的投影結(jié)果，其中橫軸為試驗次數(shù)n，縱軸為提取特征后將特征維數(shù)降為一維時的投影幅值，其中（c）和（d）為多項式核提取特征后的二維投影，不同類別分別用星號和圓點表示。而最后分類并沒有降維。

表2 數(shù)據(jù)二在4種方法下的分類錯誤率（％）Tab.2 The error rate（％）on dataset II with four different methods

4種方法下的訓(xùn)練和測試結(jié)果見表3，高斯核中KCSP的參數(shù)2σ2=2 500，KDA的參數(shù)2σ2=2 500，GDA的參數(shù)2σ2=10 000，GKLDA的參數(shù)2σ2=10 000。GKLDA正確率較高。實驗中試著改變（13）式中α的列數(shù)，在線性核下有明顯效果，減少α的列數(shù)能明顯提高正確率，當(dāng)α全選時錯誤率為40.0％，而當(dāng)α為210列時，錯誤率降為21.3％，所以在此α的列數(shù)取為210；多項式核情況下沒有明顯效果，所以在此α全選；高斯核函數(shù)α列數(shù)取為450，比全部的540正確率提高近15.0％，所以在此α保留450列。本數(shù)據(jù)集用普通LDA特征提取的分類錯誤率為39.1％；文獻(xiàn)［12］中采用Chiu＇s FIS（模糊推理系統(tǒng)）方法，錯誤率為25.3％；文獻(xiàn)［13］中采取自回歸的方法，錯誤率為17.4％。

圖3 數(shù)據(jù)三的投影結(jié)果。（a）Linear核訓(xùn)練數(shù)據(jù)；（b）Linear核測試數(shù)據(jù)；（c）Polynomial核訓(xùn)練數(shù)據(jù)；（d）Polynomial核測試數(shù)據(jù)；（e）Gaussian核訓(xùn)練數(shù)據(jù)；（f）Gaussian核測試數(shù)據(jù)Fig.3 The Projection of data III.（a）Training data with linear kernel；（b）Testing data with linear kernel；（c）Training data with polynomial kernel；（d）Testing data with polynomial kernel；（e）Training data with Gaussian kernel；（f）Testing data with Gaussian kernel

5.2 實驗二

數(shù)據(jù)格式說明：受試者為一健康男性，采用gtec信號放大器，C3和C4導(dǎo)聯(lián)，時間長為3 s，采樣率為125Hz，任務(wù)為根據(jù)提示想象左右手運(yùn)動。共采集樣本300，其中訓(xùn)練樣本200，左右手各100，測試樣本100。4種方法下的訓(xùn)練和測試結(jié)果見表4，多項式核中設(shè)置d=1.5，分類效果較為明顯；高斯核中KCSP的參數(shù)2σ2=10 000，KDA的參數(shù)2σ2=2 500，GDA的參數(shù)2σ2=10 000，GKLDA的參數(shù)2σ2=10 000。GKLDA正確率較高。本數(shù)據(jù)集用普通LDA特征的分類錯誤率為13.0％。

表3 數(shù)據(jù)三在4種方法下的分類錯誤率（％）Tab.3 The error rate（％）on dataset Ⅲ with four different methods

表4 數(shù)據(jù)四在3種方法下的分類錯誤率（％）Tab.4 The error rate（％）on dataset Ⅳ with four different methods

6 討論和結(jié)論

上述3組公開數(shù)據(jù)集中，所提出的GKLDA方法均能在適當(dāng)選擇核函數(shù)和其適當(dāng)參數(shù)后取得較好結(jié)果。相較于其他核方法（KCSP、KDA、GDA）和普通線性判別分析，競賽二中本方法略優(yōu)，而競賽一競賽三中本方法正確率明顯超過其他核方法。每組實驗結(jié)果中，本方法并不是最優(yōu)的，最優(yōu)的方法往往只適合一組數(shù)據(jù)，不具有通用性。而GKLDA往往能在適當(dāng)調(diào)節(jié)參數(shù)后適應(yīng)各類數(shù)據(jù)（也即個體差異），這是其他無參數(shù)方法所不具備的。所以從通用性角度來講，GKLDA更為突出。

核函數(shù)的參數(shù)選擇方面，多項式核d=1.5；高斯核函數(shù)的σ參數(shù)選擇過程如下，在交叉驗證時首先以10的冪次為階數(shù)，大致找到一個范圍。不失一般性以數(shù)據(jù)三為例，首先將訓(xùn)練數(shù)據(jù)按9∶1的比例分為訓(xùn)練數(shù)據(jù)和測試數(shù)據(jù)，然后參數(shù)σ的選擇范圍在10－4～105，據(jù)此，可以判斷該數(shù)據(jù)的參數(shù)范圍大致應(yīng)該在101～104之間，且最有可能在102附近，故下一步可進(jìn)行交叉驗證選取最優(yōu)參數(shù)。α的列數(shù)選擇方面，根據(jù)多組實驗結(jié)果，多數(shù)情況下是選取全部列數(shù)為優(yōu)，數(shù)據(jù)三中α的列數(shù)選取也是經(jīng)過大量測試選取的，因為訓(xùn)練數(shù)據(jù)中某些樣本的識別率較低（這從其降維后的投影中可以看出來），故在選取濾波器時稍有差異。從3種核函數(shù)的實驗結(jié)果可以看出，多項式核則較為不穩(wěn)定，結(jié)果稍差；線性核沒有可調(diào)參數(shù)，故而有較強(qiáng)的普遍性，實驗結(jié)果較好；而高斯核有參數(shù)可調(diào)（σ），對不同的樣本（即不同的人）有較強(qiáng)的針對性，故實驗結(jié)果最好。

腦電信號分類中，由于信號與信號源之間的復(fù)雜性，其大多數(shù)有用特征跟源之間都是非線性的關(guān)系，所以這就限制了CSP和線性判別分析的應(yīng)用。本研究提出的廣義核線性判別分析方法，核心在于結(jié)合了核方法和廣義奇異值分解，核方法的加入解決了一般方法對線性要求的局限性，在核函數(shù)的選擇上，由于高斯核函數(shù)的參數(shù)可變性，針對個體差異性找出相應(yīng)參數(shù)值，取得了較好的結(jié)果，又由于GSVD的應(yīng)用而減弱了對非奇異的嚴(yán)格要求，從而解決了小樣本采樣問題，使分類結(jié)果更加有效。實驗表明該方法在腦電信號處理中值得推廣。未來進(jìn)一步的工作可以放在將該方法推廣至多類別分類中。

［1］Wolpaw.JR，Birbaumer N，McFarlandDJ，etal.Braincomputer interfaces for communication and control［J］.Clinical Neurophysiology，2002，113（6）：767－791.

［2］Wang Yijun，Gao Shangkai，Gao Xiaorong.Common spatial pattern method for channel selection in motor imagery based brain-computer interface［C］//IEEE EMBS 27th Annual Conference 2005.Shanghai：IEEE，2005：5392－5395.

［3］Sun Shiliang，Zhang Changshui.An optimal kernel feature extractor and its application to EEG signal classification［J］.Neurocomputing 2005，69（13－15）：1743－1748.

［4］高湘萍，許丹，吳小培.基于核Fisher判別分析的意識任務(wù)識別新方法［J］.計算機(jī)技術(shù)與發(fā)展，2006，16（9）：82－84.

［5］王金甲，周麗娜，趙玉超.基于MEG的腦機(jī)接口特征提取方法研究［J］.儀器儀表學(xué)報，2010，31（7）：1460－1465.

［6］Nasihatkon B，Boostani R，Jahromi M.An efficient hybrid linear and kernel CSP approach for EEG feature extraction［J］.Neurocomputing.2009，73（1－3）：432－437.

［7］趙海龍，穆志純，張霞.基于LDA/GSVD和支持向量機(jī)的人耳識別［J］.上海理工大學(xué)學(xué)報，2009，31（6）：601－608.

［8］Baudat G，Anouar F.Generalized discriminant analysis using a kernel approach［J］.Neural Comput，2000，12（10）：2385－2404.

［9］邊肇祺，張學(xué)工.模式識別［M］.北京：清華大學(xué)出版社，2000：84－90.

［10］Wang Yanxia，Ruan Qiuqi.Kernel fisher discriminant analysis forpalmprintrecognition［C］// IEEE 18th International Conference on Pattern Recognition.Washington DC：IEEE，2006：457－460.

［11］Wang Yijun；Zhang Zhiguang，Li Yong，et al.BCI competition 2003—data set IV：an algorithm based on CSSD and FDA for classifying single-trial EEG［J］.IEEE transactions on biomedical engineering 2004，51（6）：1081－1086.

［12］Lotte F，Lécuyer A，Lamarche F，et al.Studying the use of fuzzy inference systems for motor imagery classification［J］.IEEE Transactions on Neural System and Rehabilitation Engineering，2007，15（2）：322－324.

［13］Brodu N，Lotte F，Lécuyer A.Comparative study of bandpower extraction techniques for motor imagery［C］，IEEE Symposium Serieson Computational Intelligence，Cognitive Algorithms，Mind，and Brain，2011：95 －100.

［14］劉廣權(quán)，黃淦，朱向陽.共空域模式方法在多類別分類中的應(yīng)用［J］.中國生物醫(yī)學(xué)工程學(xué)報，2009，28（6）：935－938.