多特征融合的運動想象腦電特征提取方法

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

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

0 引言

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

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

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

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

1 實驗數(shù)據(jù)

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

圖1 單次信號采集過程Fig.1 Single signal acquisition process

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

2 基于HCHT的特征提取方法

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

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

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

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

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

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

圖2 HCHT特征提取過程Fig.2 Feature extraction process of HCHT

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

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

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

2.2 Hilbert譜分析

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

則可求得解析信號：

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

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

則可描述信號幅度在時-頻域的分布情況，即Hilbert譜：

其中，Re為取實部。

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

式中：[ω1，ω2]為信號的頻率范圍，[t1，t2]為信號的時間范圍。

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

2.3 基于CSP的IMF特征提取

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

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

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

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

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

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

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

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

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

步驟2 特征提取。

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

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

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

2.4 運動想象腦電特征向量構(gòu)造

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

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

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

3 線下實驗

3.1 多特征和單一特征的對比

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

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

圖4 單次實驗準(zhǔn)確率對比Fig.4 Comparison of single test accuracy

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

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

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

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

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

4 線上實驗

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

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

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

圖6 實驗路徑Fig.6 Experimental path

圖7 基于三種方法的輪椅運動軌跡曲線Fig.7 Wheelchair trajectory curves based on three methods

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

圖8 采用HHT、CSP以及HCHT控制輪椅運動的耗時對比Fig.8 Time-consuming comparison of using HHT，CSP and HCHT to control wheelchair movement

5 結(jié)語

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