駱金晨，姜 月，胡秀枋，鄒任玲,2*

(1.上海理工大學(xué) 醫(yī)療器械與食品學(xué)院,上海 200093;2.民政部神經(jīng)功能信息與康復(fù)工程重點(diǎn)實(shí)驗(yàn)室 ，上海 200093)

腦-機(jī)接口(Brain-Computer Interface,BCI)將人類的意圖轉(zhuǎn)換為控制信號，實(shí)現(xiàn)人腦與外部設(shè)備之間建立信息交流與控制[1]。該系統(tǒng)在輔助控制[2]、疾病康復(fù)[3-4]、生活娛樂[5]等領(lǐng)域均有著廣泛的應(yīng)用前景。而BCI技術(shù)實(shí)現(xiàn)的關(guān)鍵依賴于腦電信號的高識別率，因此腦電信號(Electroencephalogram,EEG)提取方法的優(yōu)化是目前BCI研究的重點(diǎn),尤其基于運(yùn)動想象腦電信號特征提取和分類器的優(yōu)化識別。

近幾年來，多特征融合因能表征更多的有效腦電信息而成為多分類運(yùn)動想象的研究熱點(diǎn)之一。其中代表性的有: Sreeja等[6]提出一種基于稀疏表示的分類方法來對運(yùn)動想象腦電信號進(jìn)行多分類識別，該方法直接從分段的運(yùn)動想象數(shù)據(jù)中計(jì)算小波能量并構(gòu)造字典，使用來自字典的稀疏表示進(jìn)行分類，獲得了90%的平均正確率。陸振宇[7]等人提出一種基于多特征融合的運(yùn)動想象腦電分類方法，提取AR模型系數(shù)、小波包熵差和方差作為特征，經(jīng)LDA獲得91.43%的識別準(zhǔn)確率。KWON-WOO H[8]等人提出了一種通過膠囊網(wǎng)絡(luò)(CapsNet)學(xué)習(xí)EEG的各種屬性來實(shí)現(xiàn)運(yùn)動想象腦電信號分類方法，實(shí)驗(yàn)驗(yàn)證了CapsNet優(yōu)于卷積神經(jīng)網(wǎng)絡(luò)(CNN)和其他傳統(tǒng)機(jī)器學(xué)習(xí)方式，并獲得了78.44%的平均正確率。鄒曉紅[9]等人提出一種基于局部均值分解(LMD)和多尺度熵(MSE)相結(jié)合的腦電信號分類識別方法，將6維的融合特征向量作為SVM的輸入，相比于采用傳統(tǒng)的共空間模式(Common Spatial Pattern,CSP)特征提取算法和AR特征提取算法，最優(yōu)識別率能達(dá)到85.21%。

在閱讀相關(guān)文獻(xiàn)的基礎(chǔ)上,本研究提出在EMD基礎(chǔ)上與希爾伯特(Hilbert Transform，HT)相結(jié)合的希爾伯特-黃變換(Hilbert-Huang Transform，HHT)，獲得重構(gòu)前三階內(nèi)蘊(yùn)模態(tài)函數(shù)(Intrinsic Mode Function,IMF)能量、邊際譜熵，用“一對一”共空間模式(OVO-CSP)獲得空間域特征，選取近似熵(Approximate Entropy)、模糊熵(Fuzzy Entropy)和樣本熵(Sample Entropy)作為非線性動力學(xué)特征，將組合的時-頻-空-非線性融合特征向量采用主成分分析(Principle Component Analysis,PCA)降維，后作為優(yōu)化后SVM分類器的輸入，獲取算法仿真識別率,來驗(yàn)證本研究方法在四分類腦電信號識別研究的可靠性。整個算法流程圖(見圖1)。

圖1 算法流程圖Fig.1 Flow chart of algorithm

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

實(shí)驗(yàn)數(shù)據(jù)采用2005年第三屆國際BCI競賽中由Graz科技大學(xué)提供的腦電數(shù)據(jù)(data set IIIa)，共采集了三名受試者的腦電信號，k3b，k6b和l1b。以第一位受試者k3b作為本次研究對象。在實(shí)驗(yàn)采集過程中，受試者以放松狀態(tài)坐在椅子上，分別執(zhí)行想象左手、右手、腳或舌頭動作，提示是隨機(jī)出現(xiàn)的，實(shí)驗(yàn)由至少6輪組成，每輪試驗(yàn)40次，四種提示在每輪中分別被顯示10次。如圖2所示，試驗(yàn)開始后，前2 s為安靜狀態(tài)，在t=2 s時，有一個聲音刺激表明實(shí)驗(yàn)開始，十字形“+”出現(xiàn)，從t=3 s，一個向左、向右、向上或向下的箭頭被顯示1 s，同時，要求受試者分別想象左手、右手、舌頭或腳部運(yùn)動，直到箭頭在t=7 s時消失。信號采用Neuroscan公司的64通道腦電放大器進(jìn)行記錄，采集通道為其中的60個通道，信號采樣率為250 Hz,采用陷波濾波器進(jìn)行1-50 Hz信號濾波。

圖2 實(shí)驗(yàn)范式的時序圖Fig.2 Experimental paradigm sequence diagram

2 特征提取方法

特征提取對于多分類運(yùn)動想象識別有著至關(guān)重要的影響。本研究進(jìn)行左手、右手、舌頭、腳的運(yùn)動想象四分類數(shù)據(jù)特征提取，其采用的特征提取方法是分別提取時-頻特征、空間域特征以及非線性動力學(xué)特征后進(jìn)行多特征融合處理。

2.1 希爾伯特-黃變換

HHT因具有完全自適應(yīng)性優(yōu)勢，非常適合分析非線性、非平穩(wěn)的EEG信號。本研究采用EMD基礎(chǔ)上與希爾伯特(Hilbert Transform，HT)相結(jié)合的希爾伯特-黃變換融合算法，主要由經(jīng)驗(yàn)?zāi)B(tài)分解 (Empirical mode decomposition, EMD)和 HT算法組成[10]。將預(yù)處理后的EEG信號經(jīng)EMD分解生成IMF，再進(jìn)行Hilbert譜分析得到每個IMF的能量譜和邊際譜熵[11]。

HHT對預(yù)處理后的腦電信號x(t)的處理步驟如下[12-13]：

(1)找出x(t)的所有極值點(diǎn)并進(jìn)行插值，計(jì)算包絡(luò)線平均值：

(1)

(2)用x(t)減去包絡(luò)線平均值m1(t)得到疑似IMF分量h1(t)，判斷h1(t)是否符合IMF兩個條件，若條件符合，則將h1(t)視為一個IMF分量c1(t)，若不符合，將h1(t)代替x(t)重復(fù)上述步驟k次找出c1(t)，用x(t)減去c1(t)得到一個去掉高頻成分的新信號r1(t)。

r1(t)=x(t)-c1(t)

(2)

(3)將r1(t)看做是x(t)，重復(fù)上述步驟，由此x(t)可表示為一組IMF和一個殘余項(xiàng)和

(3)

(4)將每一個IMF進(jìn)行Hilbert變換得到相應(yīng)的Hilbert譜疊加得到EEG信號的Hilbert譜，其Hilbert變換為：

(4)

2.2 CSP優(yōu)化擴(kuò)展

采用“一對一”(One Versus One，OVO)方法對CSP算法做多分類擴(kuò)展[14]。采用OVO-CSP方法對四分類運(yùn)動想象腦電信號進(jìn)行空域特征提取。該方法將多分類拆分為若干個二分類問題，故對用于二分類的CSP傳統(tǒng)算法具體實(shí)現(xiàn)過程做說明。

(1)求兩類數(shù)據(jù)的空間協(xié)方差矩陣

(5)

(6)

(3)對混合空間協(xié)方差矩陣按式進(jìn)行特征值分解、白化處理得具有相同特征向量的S1和S2，再對特征向量S1和S2分別進(jìn)行特征值分解處理。

S1=Bλ1BTS2=Bλ2BT

(7)

B是S1和S2共同的特征向量，特征值之和為1。

(4)構(gòu)建空間濾波器后對腦電信號矩陣Ei濾波得Zi，將Zi進(jìn)行如下運(yùn)算后作為特征值：

(8)

其中p=1,2,....,2m(2m

2.3 近似熵

采用近似熵作為非線性動力學(xué)特征之一，是基于近似熵具有量化EEG信號的規(guī)律性和不可預(yù)測性優(yōu)勢,可以表示EEG信號的復(fù)雜度，反映信號中新信息發(fā)生的可能性大小[16]。近似熵對預(yù)處理后EEG信號x(t)的處理步驟如下：

(1)對N維原始信號時間序列為等時間間隔采樣，重構(gòu)m維向量X(1),X(2),...,X(N-m+1)，其中X(i)=[u(i),u(i+1),...,u(i+m-1)]。

(2)對于i=1,2,...,N-m+1，統(tǒng)計(jì)滿足以下

滿足條件的向量個數(shù)

(9)

(3)定義

(10)

(4)則可將近似熵定義為

ApEn=Φm(r)-Φm+1(r)

(11)

式中通常設(shè)置參數(shù)m=2或m=3，m=3能更細(xì)致地重構(gòu)系統(tǒng)的動態(tài)演化過程；r值主要取決于應(yīng)用場合，通常選擇r=0.2*std(std為時間序的標(biāo)準(zhǔn)差)。

2.4 模糊熵

采用模糊熵作為非線性動力學(xué)特征之一，模糊熵是于沐涵等[17]人針對提取熵值不連續(xù)問題而提出的一種改進(jìn)算法，在繼承樣本熵優(yōu)點(diǎn)的同時，減少對時間序列長度的依賴性。其因具有較好的連續(xù)性和魯棒性，可有效用于腦電時間序列的分析中。模糊熵對預(yù)處理后腦電信號x(t)的處理步驟如下：

(1)給定的N維信號時間序列與近似熵相同，定義相空間維數(shù)為m(m

X(i)=[u(i),u(i+1),...,u(i+m-1)]-u0(i)

(12)

(2)引入模糊關(guān)系函數(shù)A(x)，并計(jì)算

(13)

(14)

(4)則可將模糊熵定義為

FuzzyEn=lnΦm(r)-lnΦm+1(r)

(15)

2.5 樣本熵

采用樣本熵作為非線性動力學(xué)特征之一，是Richman等[18]人基于樣本熵是在近似熵的基礎(chǔ)上改進(jìn)而來的一種新的統(tǒng)計(jì)族，它是一種衡量信號時間序列復(fù)雜度的算法[19]。樣本熵對預(yù)處理后腦電信號x(t)的處理步驟如下：

(16)

(3)令k=m+1，重復(fù)樣本熵的第一步和第二步，可得

(17)

(4)則可將樣本熵定義為

(18)

3 運(yùn)動想象腦電信號的處理

3.1 預(yù)處理

腳部運(yùn)動想象與中央的Cz區(qū)域相關(guān)，其頻率段為7～8 Hz和20～24 Hz，舌頭運(yùn)動想象與C3、C4較接近，其頻率段為10～11 Hz。本研究通過相關(guān)的計(jì)算值貢獻(xiàn)，確定使用24、25、28(C3)、31(Cz)、34(C4)、41共六個通道數(shù)據(jù)。采用FIR等波紋濾波器分別對包含180次的試驗(yàn)樣本進(jìn)行7～30 Hz帶通濾波。本次試驗(yàn)選取每次8 s試驗(yàn)中的4～7 s數(shù)據(jù)進(jìn)行研究是因?yàn)? s以前的EEG信號與運(yùn)動想象無關(guān)，并且與受試者的反應(yīng)時間發(fā)生重疊。左右手運(yùn)動想象時預(yù)處理后的六個通道信號(見圖3)。

3.2 特征提取與融合

(1)將每次試驗(yàn)的6個通道的原始數(shù)據(jù)經(jīng)EMD分解，以某次試驗(yàn)第一個通道的EMD分解情況為本次特征提取對象，該通道的IMF的波形(見圖4a)，可明顯看出能量主要集中在前三階IMF。因EMD具有完備性，為避免出現(xiàn)特征維數(shù)過高的情況，可將EMD分解后信號進(jìn)行重構(gòu)處理。該信號前三階IMF正向重構(gòu)后波形，重構(gòu)后的信號包含絕大部分有用腦電信號(見圖4b)。將得到的6維能量特征向量F1=[E1,E2,E3,E4,E5,E6]進(jìn)行Hilbert變換，得到時頻譜和邊際譜，并計(jì)算邊際譜熵，得到每次樣本的6維邊際譜熵特征向量F2=[S1,S2,S3,S4,S5,S6]。Hilbert譜對時間積分的Hilbert邊際譜圖(見圖5)，可明顯看出每個頻率點(diǎn)的累積幅值分布，較準(zhǔn)確的反應(yīng)了信號的實(shí)際頻率成分。

圖3 左手運(yùn)動想象的預(yù)處理后信號Fig.3 Preprocessed signal of left hand movement imagination

圖4 EMD分解得到IMF和 IMF1至IMF3正向重構(gòu)后信號Fig.4 EMD decomposition to IMF and IMF1 to IMF3 are forward reconstructed

(2)使用OVO-CSP算法提取多空域特征值，分別對 “左手與右手”、“左手與舌頭”、“左手與腳”、“右手與舌頭”、“右手與腳”和“舌頭與腳”6個組合進(jìn)行CSP處理，獲得 48維特征向量，即F3=[f1,f2,f3,...,f48]。

(3)使用近似熵、模糊熵跟樣本熵作為本研究的非線性動力學(xué)特征。三種算法中的嵌入維數(shù)m和相似容限值r的參數(shù)設(shè)分別設(shè)為2和0.50，獲得3維非線性動力學(xué)特征向量，即F4=[ApEn,FuEn,SamEn]。180次樣本近似熵、模糊熵跟樣本熵的對比圖(見圖6)，可以看出藍(lán)色的近似熵整體值偏低，綠色的模糊熵相對偏高，且同一熵值在不同樣本之間的差異也較明顯，證明將其作為特征向量是可行的。

時-頻、空間域和非線性動力學(xué)特征向量融合成63維特征向量F，高維特征向量由于容易出現(xiàn)“維數(shù)災(zāi)難”，且不利于后期分類器分類，故需要通過降維來剔除對分類貢獻(xiàn)值較小的特征向量，本次試驗(yàn)把特征向量貢獻(xiàn)率設(shè)置為85.00%，將特征向量F經(jīng)PCA降維得到16維新的特征向量。

圖5 HHT邊際譜圖Fig.5 HHT marginal spectrum

圖6 三種非線性動力學(xué)特征對比圖Fig.6 Comparison of three nonlinear dynamics characteristics

4 基于SVM的腦電信號分類識別

目前，運(yùn)動想象識別分類主要采用(Linear discriminant analysis,LDA)線性判別分析和SVM兩種方法[20]，有文獻(xiàn)表明，基于參數(shù)尋優(yōu)的SVM在單次運(yùn)動想象識別分類中的效果較好[21]。本研究選取了網(wǎng)格搜索算法(Grid Search，GS)、遺傳算法(Genetic Algorithm，GA)、粒子群優(yōu)化算法(Particle Swarm Optimization，PSO)三種SVM尋優(yōu)法進(jìn)行對比研究。

基于GS尋優(yōu)法得到的GS-SVM多分類參數(shù)選擇結(jié)果等高線(見圖7)，不同的懲戒因子C與核參數(shù)g組合得到不同的準(zhǔn)確率，當(dāng)懲戒因子C為2.00，核參數(shù)g為0.35時可取得最大分類準(zhǔn)確率為81.70%。基于GA尋優(yōu)法得到的GA參數(shù)尋優(yōu)適應(yīng)度曲線(見圖8a)，從圖中可以看出隨著迭代次數(shù)的不斷增大，適應(yīng)度即分類準(zhǔn)確率也在不斷增加，最佳適應(yīng)度在進(jìn)化至45代后趨于穩(wěn)定，當(dāng)懲戒因子C約為2.56、核參數(shù)g約為0.29時最優(yōu)，可得最大分類準(zhǔn)確率為91.70%?；赑SO尋優(yōu)法得到的PSO參數(shù)尋優(yōu)適應(yīng)度曲線(見圖8b)，從圖中可以看出最佳適應(yīng)度隨著進(jìn)化代數(shù)的不斷增加而增加，在36代后趨于穩(wěn)定，在進(jìn)化至200代時迭代結(jié)束，當(dāng)懲戒因子C約為8.50、核參數(shù)g約為24.45時最優(yōu)，可得到最大分類準(zhǔn)確率為93.00%。

圖7 多分類GS-SVM參數(shù)選擇結(jié)果(等高線圖)Fig.7 Multi-class GS-SVM parameter selection results (contour map)

圖8 GA參數(shù)尋優(yōu)適應(yīng)度曲線圖和PSO參數(shù)尋優(yōu)適應(yīng)度曲線圖Fig.8 GA parameter optimization fitness curves and PSO parameter optimization fitness curve

將數(shù)據(jù)集中的k6b、l1b和k3b三個受試者數(shù)據(jù)經(jīng)多特征提取融合后，分別送入GS-SVM、GA-SVM、PSO-SVM分類算法中識別分類。實(shí)驗(yàn)顯示PSO-SVM的識別準(zhǔn)確率最高，其最大分類準(zhǔn)確率為93.30%，故選取PSO-SVM為本研究的最終識別分類算法。三種參數(shù)尋優(yōu)算法下的識別率對比情況(見表1)。

表1 三種參數(shù)尋優(yōu)算法下的識別率Table 1 Recognition rates of three parameter optimization algorithms %

5 結(jié)論與分析

為了驗(yàn)證本研究提出算法的有效性，一是以BCI2005 Data set IIIa的k3b數(shù)據(jù)集為例進(jìn)行驗(yàn)證，二是通過運(yùn)動想象腦電采集實(shí)驗(yàn)獲得四名受試者數(shù)據(jù)進(jìn)行驗(yàn)證。

對k3b數(shù)據(jù)使用PSO-SVM分別對單一特征、兩兩組合特征、三組合特征和四組合特征進(jìn)行識別分類結(jié)果(見表2)，當(dāng)四種特征結(jié)合時，識別率達(dá)到了89.90%，但該特征向量維數(shù)過高，易造成“維數(shù)災(zāi)難”，而本研究采用的多特征融合方法，達(dá)到了93.30%的準(zhǔn)確率，明顯高于其它15種方法，驗(yàn)證了將多特征融合與PSO-SVM結(jié)合作為多分類腦電信號識別的有效性。

表2 不同特征的識別率對比Table 2 Comparison of recognition rates of different features %

采用美國Neuro Scan公司的腦電信號采集系統(tǒng)進(jìn)行運(yùn)動想象腦電采集實(shí)驗(yàn)，共獲得四名身體健康的受試者腦電數(shù)據(jù)，包括2名男生和2名女生，分別記作A、B、C、D。實(shí)驗(yàn)一天內(nèi)完成，將實(shí)驗(yàn)分兩輪進(jìn)行，每人每輪進(jìn)行100次運(yùn)動想象，其中想象左手拍籃球、右手拍籃球、舌頭動、右腳踢足球各25次，構(gòu)成四組各200次樣本的數(shù)據(jù)集，選擇140個樣本為訓(xùn)練集，60個樣本為測試集。將四名受試者數(shù)據(jù)經(jīng)多域特征融合處理后分別送入OVO-SVM、GS-SVM、GA-SVM、PSO-SVM、ELM不同分類器識別，實(shí)驗(yàn)結(jié)果可以看出基于多特征融合的PSO-SVM取得四名受試者平均識別率為72.96%，在五種分類器中識別率最高(見圖9)。

圖9 多分類數(shù)據(jù)使用不同分類器識別率對比圖Fig.9 Comparison of multi-class data using different classifier recognition rates

利用EMD、HHT、OVO-CSP、近似熵、模糊熵和樣本熵進(jìn)行多域特征提取，獲取EEG信號的時-頻域、空間域和非線性動力學(xué)特征信息，后經(jīng)PCA降維后，將多特征融合向量送入PSO-SVM算法中分類。利用國際標(biāo)準(zhǔn)數(shù)據(jù)集BCI2005 Data set IIIa提供的四分類腦電數(shù)據(jù)使用本文提出的算法，最終得到93.30%的識別率。通過兩種具體方式驗(yàn)證了多特征融合及PSO-SVM分類算法的有效性，為便攜式腦機(jī)接口提供了一種新思路。本研究不僅對非穩(wěn)定、非線性類信號的分析研究具有一定的借鑒意義，而且對于腦電信號與其它多生理電信號融合研究也具有一定的參考價值，未來可增加誘發(fā)腦電的刺激種類，引入其他類型的信號，以獲得多種待識別的模式。