林海波，龔 璐，張 毅，羅 元

（重慶郵電大學 信息無障礙工程研發(fā)中心，重慶400065）

0 引 言

腦電信號主要是皮層內(nèi)大量神經(jīng)元突觸后電位同步總和所形成的，是許多神經(jīng)元共同活動的結(jié)果［1］。目前已出現(xiàn)使用腦電控制輪椅、機器人等［2－4］，但由于腦電信號的非線性和非平穩(wěn)性，致使信號的正確識別率較低，嚴重影響了特征提取算法在人機交互中的應(yīng)用。

近幾年一些特征提取方法包括AR 模型［5］、共同空間模式［6，7］、小波變換［8］、樣本熵［9］、希爾伯特－黃變換 （Hilbert－Huang transform，HHT）［10，11］等被用來分析腦電信號，雖然在一定程度上提高了信號的正確識別率，但還不能滿足實時控制的要求。如果直接采用HHT 對腦電信號進行特征提取，那么在人為選取固有模態(tài)函數(shù) （intrinsic mode function，IMF）時，會導致在重構(gòu)信號中加入噪聲，丟失有用信息。同時大腦是典型的非線性系統(tǒng)，腦電信號中還體現(xiàn)出明顯的非線性特征，樣本熵則能較好地反映腦電信號的非線性特征。基于此，本文對HHT 算法進行改進，并結(jié)合樣本熵提取想象左右手運動腦電信號特征，使用支持向量機 （support vector machine，SVM）對兩類腦電信號進行分類。

1 腦電信號特征提取方法

1.1 HHT

經(jīng)驗?zāi)B(tài)分解 （empirical mode decomposition，EMD）是HHT 算法中最核心的一步，EMD 是將復雜的非平穩(wěn)信號分解為有限個IMF的和［12，13］。假設(shè)任何復雜的非平穩(wěn)信號都是通過一組不同的IMF組成，每個IMF的極值點數(shù)和過零點數(shù)相等或者最多相差一個，同時在任意時刻，其上下包絡(luò)線對稱于時間軸，且任何兩個IMF之間是相互獨立的。那么IMF的疊加可以重構(gòu)原信號

式中：X（t）——原始信號；ci（t）——進行EMD 時第i次篩選出的IMF，它代表了信號從高到低各個頻率段的成分，rn（t）為最后的殘余分量。

對IMF進行希爾伯特變換

式中：P——柯西主值，然后構(gòu)造出解析信號，即

得到隨時間變化的瞬時幅值ai（t）瞬時相位θi（t）和瞬時頻率fi（t）

這樣原始數(shù)據(jù)可表示成

1.2 樣本熵

HHT 能提取出腦電信號的時頻特征，但不能反映出腦電信號的非線性特征。而樣本熵是一種非線性動力學中的分析方法，通過衡量時間序列中產(chǎn)生新模式概率的大小，得出腦電信號的復雜度，來描述腦電信號的非線性特征，但它卻不能反映腦電信號的時頻特征。其算法步驟如下：

（1）設(shè)原始序列為u（1），u（2），u（3），…，u（N－1），u（N），共N 個采樣點，按序號連續(xù)順序組合成一組m 維的矢量，即Xm（i）＝［u（i），u（i＋1），…u（i＋m－1）］，i＝1，2，…，N－m。

（2）定義矢量Xm（i）與Xm（j）之間的距離為兩個矢量對應(yīng)元素相減后差值最大的一個

式中：i，j＝1～N－m，i≠j。

（3）給定相似容限r(nóng)（r ＞0），統(tǒng)計距離d［Xm（i），Xm（j）］＜r，（i，j＝1～N－m，i≠j）的個數(shù)記作（r），以及此數(shù)目與總的距離數(shù)N－m－1 的比值，記作Bmi（r），即

（4）當i從1取到N－m 時，求Bmi（r）的平均值，記作Bm（r）

（5）將矢量維數(shù)增加1，構(gòu)造m＋1維矢量，按照以上所述步驟 （1）～ （4）計算 出（r），然后再求出Bm＋1（r），則序列的樣本熵值為 SampEn（m，r） ＝在實際計算時，N 取有限值，計算得到樣本熵估計值為

2 腦電數(shù)據(jù)的采集

實驗中所使用的傳感器是Emotiv，其佩戴方式如圖1所示，其采樣頻率為128Hz，各個傳感器位置按照國際10－20標準電極安放法放置，如圖2 所示，其中 “DRL”和“CMS”為參考電極。實驗選擇采集FC5和FC6 兩個通道的腦電信號，采集了想象左手運動和想象右手運動的腦電信號樣本各120組，每組樣本持續(xù)時間2s。

圖1 Emotiv信號采集儀及佩戴方式

3 實驗數(shù)據(jù)分析與處理

人在進行一側(cè)的肢體運動想象時，在大腦同側(cè)的初級感覺運動皮層區(qū)的μ/β節(jié)律 （8－26Hz）幅值出現(xiàn)增高，這種現(xiàn)象稱為時間相關(guān)同步 （event－related synchronization，ERS），而在大腦對側(cè)的初級感覺運動皮層區(qū)的μ/β節(jié)律幅值出現(xiàn)降低，這種現(xiàn)象稱為事件相關(guān)去同步 （event－related desynchronization，ERD）。

圖2 Emotiv電極安放位置

特征提取過程主要采用HHT 和樣本熵兩個算法。HHT 算法實現(xiàn)對腦電信號8－26Hz的濾波處理，提取代表μ/β節(jié)律部分的IMF，并計算由有效IMF 重構(gòu)的腦電信號的能量均值，同時還計算出重構(gòu)腦電信號的樣本熵值組合成特征向量，整個改進后的特征提取算法流程如圖3所示。

圖3 特征提取算法流程

3.1 改進HHT的特征提取

首先對原始腦電數(shù)據(jù)進行零均值化，去除原始腦電信號中低頻直流分量的干擾，然后對腦電信號樣本進行EMD，得到一組IMF，如圖4所示是受試者的一次想象右手運動的腦電數(shù)據(jù)樣本在FC6通道經(jīng)過EMD 后得到的各個IMF。

圖4中可以看到，腦電信號中最先被分離出來的是高頻部分，隨著IMF階數(shù)的增加，被分離出來的信號頻率逐漸降低。不同的腦電信號樣本所得到的IMF階數(shù)可能會不一樣，這充分體現(xiàn)出了EMD 對信號的自適應(yīng)性。

在傳統(tǒng)HHT 的特征提取過程中，通常是選取前幾階IMF來進行分析，為了克服該問題，本文提出了選取有效IMF的兩個條件。

（1）引入相關(guān)系數(shù)，以各個IMF 與進行EMD 前的原始腦電信號的相關(guān)系數(shù)為判定依據(jù)，來選取有效的IMF，即

圖4 EMD 后的各個IMF

式中：r——一個IMF 與原始腦電信號的相關(guān)系數(shù)，c——IMF，x——原始腦電信號，n＝256。相關(guān)系數(shù)越大，表明該IMF中含原始腦電信號的有效成分就越高。滿足選取IMF的條件1為

（2）確定IMF對應(yīng)的頻率范圍。對滿足式 （13）的每一個IMF，計算IMF中各個時刻的瞬時頻率，由于一組樣本有256個數(shù)據(jù)，那么一個IMF 可以計算出256個瞬時頻率，分別統(tǒng)計每一個IMF中的瞬時頻率fi（t）在μ/β節(jié)律頻帶內(nèi)的個數(shù)，即滿足（8≤fi（t）≤26 Hz），并記為mj，j是滿足式 （13）的IMF的序號。計算mj與該IMF中瞬時頻率個數(shù)的比值，設(shè)定閾值δ。即選取的IMF需要滿足的條件2為

當δ＝0.7時，得到的特征效果較好。綜上所述選取的IMF必須滿足式 （13）和式 （14）。

對選取的有效IMF進行信號重構(gòu)，即

式中：y——重構(gòu)后的信號，i——有效IMF的序號。

同時對有效IMF進行希爾伯特變換得到希爾伯特邊際譜，如圖5所示。從圖5 （a）、（b）兩圖可以看到想象左手運動腦電信號的幅值在FC5通道要高于FC6通道，而想象右手運動腦電信號幅值在FC6通道要高于FC5通道，其頻率基本在8－26Hz內(nèi)，符合μ/β節(jié)律頻帶范圍。

圖5 改進HHT 算法的希爾伯特邊際譜

本文也對傳統(tǒng)HHT 算法中選取的前3階IMF 進行了希爾伯特變換，得到希爾伯特邊際譜如圖6 所示。在圖6（a）、（b）兩圖可以看到希爾伯特邊際譜混入了較明顯的高頻和低頻噪聲，想象左右手運動在FC5和FC6兩個通道區(qū)別并不是很明顯。

圖6 傳統(tǒng)HHT 算法的希爾伯特邊際譜

根據(jù)上述結(jié)果，將FC5和FC6通道重構(gòu)的腦電信號的能量均值作為提取的時頻特征，即

3.2 結(jié)合樣本熵的特征提取

樣本熵算法中取m ＝2，相似容限r(nóng)＝0.2SD （SD 表示原始序列的標準差），分別計算在FC5 和FC6 通道由式（15）重構(gòu)的腦電信號的樣本熵，提取特征腦電信號的非線性特征，并表示為SampEnFC5、SampEnFC6。

如圖7 （a）、（b）兩圖是想象左右手運動的腦電信號各120組樣本在FC5和FC6通道計算得到的樣本熵。

圖7 有效IMF重構(gòu)信號的樣本熵值

在圖7中可以觀察到，只有較少的一部分樣本在兩個通道還不能區(qū)分，大部分的數(shù)據(jù)樣本都表現(xiàn)出在想象左手運動時FC5通道的樣本熵值低于FC6通道，而想象右手運動時FC6通道的樣本熵值低于FC5通道，這恰好與腦電信號中μ/β節(jié)律幅值變化相反，說明使用樣本熵能夠提取出想象左右手運動腦電信號的非線性特征，可以有效提高信號的正確識別率。

將樣本熵算法對腦電信號提取的2維特征向量和改進HHT 算法提取的2維特征向量組成最終的4維特征向量，形成腦電信號的多特征提取。即 ［EFC5，EFC6，SampEn－FC5，SampEnFC6］。

3.3 算法性能測試

本文采用支持向量機作為模式分類器對提取的特征向量進行分類，其中核函數(shù)采用徑向基核函數(shù)，為了能夠合理地測試分類器的識別率，使用交叉驗證來獲得最優(yōu)的核參數(shù)γ和懲罰因子C，將樣本集隨機地分成5份，輪流將其中4份作為訓練樣本，另外1份作為測試樣本，獲得了5次測試結(jié)果的平均離線識別率。表1給出了使用改進HHT結(jié)合樣本熵、傳統(tǒng)HHT 結(jié)合樣本熵以及在單一HHT、樣本熵算法下得到的正確識別率。

表1 不同特征提取方法下的正確識別率

在文獻 ［5］中采用AR 模型算法獲得最高正確識別率為87.857%，而在文獻 ［7］中采用CSP 算法獲得最高正確識別率僅為82.86%，這表明本文提出的改進HHT 結(jié)合樣本熵的特征提取方法是有效的。

將上述的分類結(jié)果作為智能輪椅的控制指令，在智能輪椅平臺上進行實驗。當輪椅開始動作后保持自行前進，受試者通過想象左手運動控制輪椅左轉(zhuǎn)、想象右手運動控制輪椅右轉(zhuǎn)。

本文讓5名受試者操作分別用HHT 結(jié)合樣本熵和改進HHT 結(jié)合樣本熵的人機交互系統(tǒng)控制智能輪椅完成指定路徑，如圖8所示，并記錄輪椅走過的路徑軌跡，進而觀察在實時控制過程中EEG 人機交互系統(tǒng)的穩(wěn)定性。

圖8 實驗路徑

實驗中每位受試者都按照圖8中箭頭所指方向控制輪椅行走。圖9給出了5位受試者操作智能輪椅時記錄下來的路徑。

由圖9可以看出采用HHT 和樣本熵的人機交互系統(tǒng)控制輪椅時，軌跡曲線有較大的波動且不光滑，是因為受試者控制輪椅時存在動作誤識別所導致的輪椅原地打轉(zhuǎn)、

圖9 5名受試者操作智能輪椅時的運動軌跡

向非意識方向動作等現(xiàn)象；而采用本文提出算法的人機交互系統(tǒng)控制輪椅時，輪椅的行走曲線光滑并且波動較小，可以順利第避開障礙物。這說明該人機交互系統(tǒng)對運動想象的意圖識別率更高、穩(wěn)定性更好。

4 結(jié)束語

本文提出了一種改進HHT 和樣本熵結(jié)合的特征提取方法，分別將HHT 結(jié)合樣本熵和改進HHT 結(jié)合樣本熵的兩種方法對想象左右手運動腦電信號進行了分析，并提取了兩類腦電信號的特征組成特征向量，最后在智能輪椅平臺上進行了實驗。研究結(jié)果表明，改進HHT 結(jié)合樣本熵的特征提取算法對想象左右手運動腦電信號的正確識別率明顯高于HHT 結(jié)合樣本熵的特征提取算法；采用基于改進算法的智能輪椅系統(tǒng)擁有更好的實時性與穩(wěn)定性。

［1］Long Jinyi，Li Yuanqi，Wang Hongtao，et al.A hybrid brain computer interface to control the direction and speed of a simulated or real wheelchair［J］.Neural Systems and Rehabilitation Engineering，2012，20 （5）：720－729.

［2］Postelnicu C，Talaba D.P300－based brain－neuronal computer interaction for spelling applications ［J］.Biomedical Engineering，2013，60 （2）：534－543.

［3］Zhang Haihong，Yang Huijuan，Guan Cuntai.Bayesian learning for spatial filtering in an EEG－based brain－computer interface［J］.Neural Networks and Learning Systems，2013，24（7）：1049－1060.

［4］Chae Y，Jeong J，Jo S.Toward brain－actuated humanoid robots：Asynchronous direct control using an EEG－based BCI［J］.Robotics，2012，28 （5）：1131－1144.

［5］HUANG Sijuan，WU Xiaoming.Feature extraction and classification of EEG based on energy characteristics ［J］.Chinese Journal of Sensors and Actuators，2010，23 （6）：782－785 （in Chinese）.［黃思娟，吳效明.基于能量特征的腦電信號特征提取與分類 ［J］.傳感技術(shù)學報，2010，23 （6）：782－785.］

［6］Robinson N，Vinod A，Ang K，et al.EEG－based classification n of fast and slow hand movements using wavelet－CSP algorithm ［J］.Biomedical Engineering，2013，60 （8）：2123－2132.

［7］LIU Chong，ZHAO Haibin，LI Chunsheng，et al.CSP/SVM－based EEG classification of imagined hand movements［J］.Journal of Northeastern University（Natural Science），2010，31（8）：1098－1101（in Chinese）.［劉沖，趙海濱，李春勝，等.基于CSP與SVM 算法的運動想象腦電信號分類［J］.東北大學學報：自然科學版，2010，31 （8）：1098－1101.］

［8］XU Baoguo，SONG Aiguo，F(xiàn)EI Shumin.Feature extraction and classification of EEG in online brain－computer interface［J］.ACTA Electronica Sinica，2011，39 （5）：1025－1030 （in Chinese）. ［徐寶國，宋愛國，費樹岷.在線腦機接口中腦電信號的特征提取與分類方法［J］.電子學報，2011，39 （5）：1025－1030.］

［9］ZHOU Peng，GE Jiayi，CAO Hongbao，et al.Classification of motor imagery based on sample entropy ［J］.Information and Control，2008，37 （2）：191－196 （in Chinese）. ［周鵬，葛家怡，曹紅寶，等.基于樣本熵的運動想象分類研究 ［J］.信息與控制，2008，37 （2）：191－196.］

［10］Chang H，Lee P，Lo M，et al.Inter－trial analysis of postmovement beta activities in EEG signals using multivariate empirical mode decomposition ［J］.Neural Systems and Rehabilitation Engineering，2013，21 （4）：607－615.

［11］Park C，Looney D，Rehman N，et al.Classification of motor imagery BCI using multivariate empirical mode decomposition［J］.Neural Systems and Rehabilitation Engineering，2013，21 （1）：10－22.

［12］Bajaj V，Pachori R.Classification of seizure and nonseizure EEG signals using empirical mode decomposition ［J］.Information Technology in Biomedicine，2012，16 （6）：1135－1142.

［13］Park C，Looney D，Kidmose P，et al.Time－frequency analysis of EEG asymmetry using bivariate empirical mode decomposition ［J］.Neural Systems and Rehabilitation Engineering，2011，19 （4）：366－373.