羅 飛，劉鵬飛，羅 元，朱思蒙

（重慶郵電大學(xué)光電工程學(xué)院，重慶400065）

0 引言

腦-機(jī)接口（Brain-Computer Interface，BCI）系統(tǒng)通過(guò)分析輸入的電生理信號(hào)，將用戶(hù)的意圖解碼成控制指令來(lái)操作輸出設(shè)備［1］。根據(jù)信號(hào)形式的不同，腦電信號(hào)主要分 為 腦 電 圖（ElectroEncephalonGraph，EEG）、腦 磁 圖（MagnetoEncephaloGram，MEG）、功能磁共振成像（Functional Magnetic Resonance Imaging，F(xiàn)MRI）等［2］。其中，EEG 由于其非侵入性和低成本的特點(diǎn)而廣泛應(yīng)用于BCI系統(tǒng)［3］。

運(yùn)動(dòng)想象（Motor Imagery，MI）的原理是想象運(yùn)動(dòng)時(shí)會(huì)對(duì)大腦兩側(cè)的EEG 信號(hào)產(chǎn)生事件相關(guān)去同步/同步現(xiàn)象［4］，已成為EEG 信號(hào)領(lǐng)域的研究熱點(diǎn)。其中，EEG 信號(hào)因其微弱且具有非線(xiàn)性、非平穩(wěn)及時(shí)變敏感等特點(diǎn)，時(shí)-頻域分析和空間濾波被廣泛應(yīng)用于其特征提取過(guò)程［5］。時(shí)-頻域分析主要有短時(shí)傅里葉變換（Short-term Fourier Transform，STFT）［6］、小波變換（Wavelet Transform，WT）［7］和小波包變換（Wavelet Package Transform，WPT）［8］，空間濾波主要為共同空間模式（Common Spatial Pattern，CSP）算法［9］。

然而，基于STFT、WT、WPT的時(shí)-頻域分析方法無(wú)法同時(shí)在時(shí)域和頻域范圍內(nèi)取得較高的分辨率。希爾伯特-黃變換（Hilbert-Huang Transform，HHT）是一種適用于非線(xiàn)性非穩(wěn)定信號(hào)的時(shí)-頻域特征提取方法，該方法通過(guò)經(jīng)驗(yàn)?zāi)B(tài)分解（Empirical Mode Decomposition，EMD）將信號(hào)自適應(yīng)地分解成多個(gè)具有物理意義的固有模態(tài)函數(shù)（Intrinsic Mode Function，IMF），通過(guò)對(duì)各階IMF 作Hilbert 變換，可以獲得具有很高分辨率的時(shí)-頻域特征［10］。但隨著大腦狀態(tài)變化，腦電信號(hào)的時(shí)-頻域特征會(huì)出現(xiàn)波動(dòng)。因此，近年來(lái)相關(guān)研究者綜合考慮多種特征方法進(jìn)行特征提取。Chen 等［11］融合香農(nóng)熵、小波熵和樣本熵進(jìn)行特征提?。粭钅龋?2］基于總體經(jīng)驗(yàn)?zāi)B(tài)分解（Ensemble Empirical Mode Decomposition，EEMD）和近似熵提出一種多特征提取方法。以上方法都表現(xiàn)出很好的自適應(yīng)性和較高的識(shí)別準(zhǔn)確率，但考慮的都是單個(gè)角度特征的融合，不能從多個(gè)角度獲取信號(hào)的更完整描述。因此，本文考慮一種綜合時(shí)-頻-空域特征的方法。

考慮到EMD可以將單個(gè)通道擴(kuò)展成多個(gè)具有物理意義的窄帶信號(hào)IMF，本文在HHT 的基礎(chǔ)上提出一種提取時(shí)-頻-空域特征的方法HCHT（Hilbert-CSP-Huang Transform）。將EMD得到的各階IMF進(jìn)行Hilbert變換，獲得具有很高分辨率的時(shí)-頻域特征，并對(duì)各階IMF 作進(jìn)一步的CSP 分解，即將各階IMF分量合并成新的信號(hào)矩陣，通過(guò)構(gòu)造的IMF信號(hào)矩陣在不同類(lèi)狀態(tài)下的幅值差異構(gòu)造空間濾波器，獲取各階IMF的空間分布特征，從而將HHT 提取的時(shí)-頻域特征擴(kuò)展為時(shí)-頻-空域特征。該方法實(shí)現(xiàn)了對(duì)腦電信號(hào)的多角度描述，且信號(hào)特征都是基于HHT算法，互為補(bǔ)充。在數(shù)據(jù)集上的實(shí)驗(yàn)結(jié)果表明，本文方法平均識(shí)別準(zhǔn)確率更高，且標(biāo)準(zhǔn)差更小，魯棒性更好。最后在智能輪椅平臺(tái)上對(duì)該方法的有效性進(jìn)行了進(jìn)一步的驗(yàn)證。

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

本文采用的數(shù)據(jù)集為BCI Competition II 的Data set III 數(shù)據(jù)集合［13］。實(shí)驗(yàn)樣本記錄來(lái)自一名25 歲的健康女性。EEG信號(hào)從國(guó)際10-20 導(dǎo)聯(lián)系統(tǒng)的C3、C4 和Cz 通道獲得，采樣頻率為128 Hz，濾波范圍為0.5～30 Hz。單次信號(hào)采集過(guò)程如圖1所示。

圖1 單次信號(hào)采集過(guò)程Fig.1 Single signal acquisition process

整個(gè)實(shí)驗(yàn)分為7 組，每組包括40 個(gè)實(shí)驗(yàn)。單次實(shí)驗(yàn)持續(xù)9 s，實(shí)驗(yàn)要求受試者安靜坐在顯示器前，根據(jù)屏幕提示進(jìn)行左手或右手運(yùn)動(dòng)想象（MI）任務(wù)。每次實(shí)驗(yàn)前2 s屏幕顯示空白，提醒受試者放松；t=2 s 時(shí)，顯示器屏幕顯示“+”字形，受試者準(zhǔn)備MI；t=3 s時(shí)，顯示器中央隨機(jī)顯示向左或向右的箭頭，受試者根據(jù)箭頭方向想象左手或右手運(yùn)動(dòng)，持續(xù)6 s；t=9 s時(shí)，單次采集結(jié)束，受試者短暫休息，準(zhǔn)備下一次實(shí)驗(yàn)。

2 基于HCHT的特征提取方法

基于HCHT 的特征提取方法如圖2 所示，具體過(guò)程分為四步：1）將預(yù)處理后的MI 腦電信號(hào)經(jīng)過(guò)EMD 得到多個(gè)IMF；2）提取各階IMF 的Hilbert 瞬時(shí)能量譜（Instanta-neous Energy Spectrum，IES）和邊際能量譜（Marginal Energy Spectrum，MES）分別作為腦電信號(hào)的時(shí)域特征和頻域特征；3）將各階IMF分量合并成新的信號(hào)矩陣，對(duì)各階IMF進(jìn)行進(jìn)一步的CSP分解，獲得表征空間信息的特征向量；4）歸一化第2）步和第3）步提取到的信息，得到表征MI腦電信號(hào)的特征向量。

2.1 經(jīng)驗(yàn)?zāi)B(tài)分解

經(jīng)驗(yàn)?zāi)B(tài)分解其本質(zhì)是通過(guò)信號(hào)的特征時(shí)間尺度判別內(nèi)蘊(yùn)振蕩模式，將信號(hào)自適應(yīng)地分解成多個(gè)具有物理意義的固有模態(tài)函數(shù)IMF［14］。

預(yù)處理后的輸入信號(hào)經(jīng)過(guò)EMD得到如下表達(dá)式：

其中：X(t)為輸入信號(hào)，Ci(t)第i 次篩選得到的IMF 分量，N 為篩選次數(shù)，Rn(t)為最終的剩余分量。

對(duì)單次MI 腦電信號(hào)進(jìn)行經(jīng)驗(yàn)?zāi)B(tài)分解，分解結(jié)果如圖3所示。

圖2 HCHT特征提取過(guò)程Fig.2 Feature extraction process of HCHT

圖3 MI腦電信號(hào)EMD結(jié)果Fig.3 EMD results of EEG indued by MI

分解后的各階IMF 都是蘊(yùn)含有MI腦電信息的窄帶信號(hào)，每一個(gè)IMF分量都可以看作單個(gè)信號(hào)通道。本文對(duì)C3、C4通道的各階IMF 分量進(jìn)行采樣（采樣頻率為原始腦電信號(hào)的采樣頻率），并將兩個(gè)通道的信號(hào)分量合并，構(gòu)造成一個(gè)N×T信號(hào)矩陣，其中N為IMF分量個(gè)數(shù)，T為信號(hào)分量的采樣點(diǎn)數(shù)，信號(hào)矩陣表示為：

其中：n為IMF個(gè)數(shù)，j=L或R。

2.2 Hilbert譜分析

MI 腦電信號(hào)的特征信息主要集中在前三階IMF，因此本文主要對(duì)前三階IMF 進(jìn)行研究。對(duì)每個(gè)IMF 分量進(jìn)行Hilbert變換：

則可求得解析信號(hào)：

其中：Ai(t)為瞬時(shí)幅值，φi(t)為瞬時(shí)相位。

由Ai(t)和φi(t)可進(jìn)一步求取瞬時(shí)頻率ωi(t)，即：

則可描述信號(hào)幅度在時(shí)-頻域的分布情況，即Hilbert譜：

其中，Re為取實(shí)部。

根據(jù)式（6）可進(jìn)一步求取IES和MES：

式中：[ω1，ω2]為信號(hào)的頻率范圍，[t1，t2]為信號(hào)的時(shí)間范圍。

IES 和MES 分別反映了信號(hào)在時(shí)域和頻域上的能量特征。本文分別將IES 和MES 定義為MI 腦電信號(hào)的時(shí)域特征F1∈Rm1×1和頻域特征F2∈Rm2×1，其中m1和m2分別表示時(shí)間點(diǎn)數(shù)和頻率點(diǎn)數(shù)。

2.3 基于CSP的IMF特征提取

共同空間模式來(lái)源于共空域子空間分解（Common Spatial Subspace Decomposition，CSSD），是一種兩分類(lèi)任務(wù)下的空域?yàn)V波算法，能夠從多通道腦電信號(hào)數(shù)據(jù)里面提取出每一類(lèi)的空間分布成分［15］，主要分為兩步：構(gòu)造空間濾波器和特征提取。

根據(jù)2.1 節(jié)信號(hào)矩陣的構(gòu)造方法，將C3、C4 通道前三階IMF分量合并，構(gòu)造左右手MI的信號(hào)矩陣XL和XR。

步驟1 構(gòu)造空間濾波器。

求出左右手MI信號(hào)矩陣的均值空間協(xié)方差矩陣：

其中：Xc，i表示左右手MI 的第i 次實(shí)驗(yàn)，K 表示同一標(biāo)簽的信號(hào)矩陣的實(shí)驗(yàn)數(shù)量，trace()表示矩陣的跡。

求混合空間協(xié)方差矩陣R，并對(duì)其進(jìn)行特征值分解：

進(jìn)一步可求取白化矩陣：

對(duì)白化后的矩陣進(jìn)行特征值分解：

其中：Us是SL和SR的特征向量，λL和λR分別是對(duì)應(yīng)的特征值矩陣。可以證明SL和SR具有相同的特征矩陣Us，且λL與λR之和為單位矩陣。因此，當(dāng)SL的特征值最大時(shí)，SR的特征值最小，可最大限度地區(qū)分兩類(lèi)信號(hào)。則可構(gòu)造出空間濾波器：

步驟2 特征提取。

對(duì)構(gòu)造的信號(hào)矩陣Xj進(jìn)行進(jìn)一步的共同空間模式分解，經(jīng)空間濾波器W投影得到特征矩陣Z，如式（14）所示：

將特征矩陣Z 的前m 行和后m 行作為信號(hào)矩陣的特征，求取特征向量：

vj表征了兩類(lèi)MI 腦電信號(hào)特征，本文由式（10）將其表示為空域特征F3∈R2m×1：

2.4 運(yùn)動(dòng)想象腦電特征向量構(gòu)造

本文的MI腦電信號(hào)由上述三種特征進(jìn)行描述，將求取的IES 和MES 以及CSP 方差向量整合并構(gòu)造輸入特征向量F'={F1，F(xiàn)2，F(xiàn)3}。由于不同的特征具有不同的含義，因此有必要通過(guò)以下方式對(duì)輸入特征向量進(jìn)行歸一化：

其中：μk是第k 個(gè)特征的平均值組成的向量，σk是第k 個(gè)特征的標(biāo)準(zhǔn)差。

最后，將歸一化的輸入特征F={F'1，F(xiàn)'2，F(xiàn)'3}通過(guò)支持向量機(jī)（Support Vector Machine，SVM）分類(lèi)。SVM 的原理是使分離超平面之間的距離最大，并選擇最佳超平面作為決策邊界［16］。由于SVM 具有出色的分類(lèi)性能，因此已廣泛應(yīng)用于腦電信號(hào)處理。

3 線(xiàn)下實(shí)驗(yàn)

3.1 多特征和單一特征的對(duì)比

為了驗(yàn)證本文方法的有效性，將數(shù)據(jù)集隨機(jī)分成10 個(gè)子集，每個(gè)子集包含來(lái)自各個(gè)分類(lèi)的相同比例的樣本。依次選取其中1個(gè)子集用于測(cè)試，其余9個(gè)子集用于訓(xùn)練，重復(fù)10次，并取10次實(shí)驗(yàn)的平均準(zhǔn)確率。為公平比較，使用10折交叉驗(yàn)證法找到每個(gè)過(guò)程中SVM的最佳懲罰系數(shù)C和核函數(shù)參數(shù)σ，并采用10次實(shí)驗(yàn)的平均準(zhǔn)確率和標(biāo)準(zhǔn)差來(lái)測(cè)量性能。單次實(shí)驗(yàn)準(zhǔn)確率如圖4所示，平均分類(lèi)準(zhǔn)確率和標(biāo)準(zhǔn)差如表1所示。

表1 平均準(zhǔn)確率Tab.1 Average accuracy

圖4 單次實(shí)驗(yàn)準(zhǔn)確率對(duì)比Fig.4 Comparison of single test accuracy

從表1 可以看出，當(dāng)使用本文的融合特征F 時(shí)，可以獲得88.0%的平均準(zhǔn)確率，與單一的時(shí)-頻域特征和空域特征相比，分別提高了7.5、10.3 和9.2 個(gè)百分點(diǎn)。這是因?yàn)榕c單一特征相比，融合的特征實(shí)現(xiàn)了信號(hào)不同特征之間的互補(bǔ)，可以獲得更全面的信號(hào)表達(dá)，進(jìn)而可以獲得更好的識(shí)別效果。并且，本文特征提取方法的標(biāo)準(zhǔn)差均小于其他的三種方法，說(shuō)明本文提出的方法具有更強(qiáng)的穩(wěn)健性。

3.2 多種識(shí)別方法結(jié)果比較

表2 還給出了文獻(xiàn)［11-12］兩種多特征提取方法的識(shí)別結(jié)果?？梢钥闯?，當(dāng)通道數(shù)相同時(shí)，本文方法的識(shí)別率分別比文獻(xiàn)［11-12］的方法分別高2.3和4.7個(gè)百分點(diǎn)。這是因?yàn)楸疚姆椒ňC合考慮了腦電信號(hào)的時(shí)-頻-空域信息進(jìn)行特征提取，從而可以獲得信號(hào)特征的更完整表達(dá)，有效提高了MI 任務(wù)的識(shí)別能力。從耗時(shí)來(lái)看，本文方法的平均耗時(shí)比第一名多0.5 s，整體差異不大，但識(shí)別率提高了4.7個(gè)百分點(diǎn)。

同時(shí)，為了體現(xiàn)對(duì)比方法的多樣性，本文還在BCI Competition III 的Data set I數(shù)據(jù)集進(jìn)行了實(shí)驗(yàn)。表2還給出了腦機(jī)接口競(jìng)賽“BCI Competition III Data Set I”前三名識(shí)別方法的結(jié)果（按名次由高到低分別為文獻(xiàn)［17］、文獻(xiàn)［18］和文獻(xiàn)［19］的方法）。本文方法取得的識(shí)別率為87.3%，分別比第二名和第三名高0.3 和1.3 個(gè)百分點(diǎn)，比第一名低3.7 個(gè)百分點(diǎn)。然而，本文提出的特征提取方法僅提取兩個(gè)通道的信號(hào)進(jìn)行分析，在保證分類(lèi)準(zhǔn)確率的情況下，大大減少了通道數(shù)。從測(cè)試集的平均分類(lèi)耗時(shí)對(duì)比可以看出，隨著通道數(shù)的減少，計(jì)算數(shù)據(jù)量大大減少，因此，本文方法的平均耗時(shí)遠(yuǎn)少于以上三種方法，為便攜式BCI 系統(tǒng)的在線(xiàn)采集和識(shí)別提供了可行性。因此，綜合來(lái)看，本文方法均優(yōu)于以上提及的方法。

表2 多種識(shí)別方法結(jié)果比較Tab.2 Result comparison of multiple recognition methods

4 線(xiàn)上實(shí)驗(yàn)

本文在智能輪椅平臺(tái)上采用基于HCHT 的人機(jī)交互系統(tǒng)進(jìn)行了在線(xiàn)實(shí)驗(yàn)。令6 位受試者通過(guò)想象左手或右手運(yùn)動(dòng)來(lái)控制輪椅運(yùn)動(dòng)。采用Emotive傳感器作為腦電信號(hào)采集設(shè)備，16 個(gè)電極的安放位置如圖5 所示。將采樣電極中的FC5 和FC6 電極作為輸入電極，CMS 和DRL 作為參考電極。在輪椅的方向控制中，腦電信號(hào)首先經(jīng)過(guò)0.1～30 Hz 的帶通濾波，然后通過(guò)HCHT 算法對(duì)信號(hào)特征進(jìn)行提取，對(duì)提取到的特征通過(guò)SVM算法進(jìn)行分類(lèi)識(shí)別，并轉(zhuǎn)化為控制指令控制輪椅轉(zhuǎn)向。

圖5 電極安放位置Fig.5 Electrode placement positions

采用基于HHT、CSP 和HCHT 的人機(jī)交互系統(tǒng)，在智能輪椅平臺(tái)上對(duì)6 位受試者進(jìn)行重復(fù)測(cè)試，完成如圖6 所示的“8”字形路線(xiàn)。實(shí)驗(yàn)場(chǎng)地兩側(cè)各設(shè)置一個(gè)障礙物，受試者端坐在輪椅上以0.15 m/s的前進(jìn)速度，從起點(diǎn)出發(fā)，通過(guò)想象左手或右手運(yùn)動(dòng)來(lái)控制輪椅轉(zhuǎn)向，繞過(guò)兩側(cè)的障礙物，得到如圖7 所示的運(yùn)動(dòng)軌跡曲線(xiàn)。從圖7（a）和（b）可以看出，基于HHT 和CSP的BCI系統(tǒng)與圖7（c）的HCHT相比，其運(yùn)動(dòng)軌跡曲線(xiàn)較為混亂，不平滑，且波動(dòng)較大。這是因?yàn)镠CHT 是一種多特征融合的方法，可以獲得MI 腦電信號(hào)特征的更完整表達(dá)，能夠提取出更為準(zhǔn)確的腦電信號(hào)特征，從而獲得更高的分類(lèi)準(zhǔn)確率，實(shí)現(xiàn)對(duì)智能輪椅更精確的控制，受試者能更安全平滑地完成指定路線(xiàn)。除此之外，還可以看出“8”字形的右邊較為混亂，這是由于受試者長(zhǎng)時(shí)間處于MI狀態(tài)會(huì)出現(xiàn)疲勞狀態(tài)，使得信號(hào)的特征值發(fā)生變化，識(shí)別準(zhǔn)確率出現(xiàn)略微下降。

圖6 實(shí)驗(yàn)路徑Fig.6 Experimental path

圖7 基于三種方法的輪椅運(yùn)動(dòng)軌跡曲線(xiàn)Fig.7 Wheelchair trajectory curves based on three methods

圖8 為6 位受試者分別在HHT、CSP 和HCHT 三種方案下控制輪椅完成“8”字形路徑的耗時(shí)對(duì)比?？梢钥闯?，基于HCHT的控制方案耗時(shí)略微多于HHT和CSP，這是因?yàn)镠CHT算法是一種多特征提取算法，同時(shí)提取IES、Hilbert 邊際能量譜以及CSP 方差向量，從多個(gè)角度實(shí)現(xiàn)了對(duì)MI腦電信號(hào)特征的綜合描述，因此耗時(shí)增加。三種方法的耗時(shí)主要集中在EMD 過(guò)程，提取IES、Hilbert 邊際能量譜以及CSP 方差向量的耗時(shí)占比并不大。因此，從圖中可以看出，雖然HCHT 耗時(shí)比HHT和CSP略多，但三者的耗時(shí)差異不大，并不影響整體系統(tǒng)響應(yīng)，且HCHT 算法具有更高的識(shí)別率和穩(wěn)定性，更適用于實(shí)際應(yīng)用環(huán)境。

圖8 采用HHT、CSP以及HCHT控制輪椅運(yùn)動(dòng)的耗時(shí)對(duì)比Fig.8 Time-consuming comparison of using HHT，CSP and HCHT to control wheelchair movement

5 結(jié)語(yǔ)

本文基于HHT 和CSP 提出的MI 腦電信號(hào)特征提取方法HCHT，綜合考慮了腦電的時(shí)間-頻率-空間信息進(jìn)行特征提取。即首先經(jīng)過(guò)EMD 得到各階IMF，接著對(duì)各階IMF 進(jìn)行Hilbert 變換，提取IES 和MES 分別作為時(shí)域特征和頻域特征，并將多個(gè)IMF 合并成新的信號(hào)矩陣，利用CSP 進(jìn)行進(jìn)一步的空間濾波，獲取其空間特征。最后，使用特征向量來(lái)訓(xùn)練SVM分類(lèi)器，并對(duì)MI腦電信號(hào)模式進(jìn)行分類(lèi)。在競(jìng)賽數(shù)據(jù)集上得到的實(shí)驗(yàn)結(jié)果表明，相比單一的時(shí)頻特征和空間特征，本文提出的HCHT 特征提取方法平均準(zhǔn)確率更高，標(biāo)準(zhǔn)差更小，魯棒性更好。最后利用智能輪椅平臺(tái)對(duì)算法進(jìn)行驗(yàn)證，進(jìn)一步驗(yàn)證了算法的有效性。本文方法實(shí)現(xiàn)了對(duì)兩通道左右手MI 腦電信號(hào)的有效識(shí)別，為便攜式BCI技術(shù)提供了一種新的思路，但識(shí)別率和識(shí)別種類(lèi)有待進(jìn)一步的優(yōu)化。