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

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

（重慶郵電大學光電工程學院，重慶400065）

0 引言

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

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

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

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

1 實驗數據

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

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

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

2 基于HCHT的特征提取方法

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

2.1 經驗模態分解

經驗模態分解其本質是通過信號的特征時間尺度判別內蘊振蕩模式，將信號自適應地分解成多個具有物理意義的固有模態函數IMF［14］。

預處理后的輸入信號經過EMD得到如下表達式：

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

對單次MI 腦電信號進行經驗模態分解，分解結果如圖3所示。

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

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

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

其中：n為IMF個數，j=L或R。

2.2 Hilbert譜分析

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

則可求得解析信號：

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

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

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

其中，Re為取實部。

根據式（6）可進一步求取IES和MES：

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

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

2.3 基于CSP的IMF特征提取

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

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

步驟1 構造空間濾波器。

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

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

求混合空間協方差矩陣R，并對其進行特征值分解：

進一步可求取白化矩陣：

對白化后的矩陣進行特征值分解：

其中：Us是SL和SR的特征向量，λL和λR分別是對應的特征值矩陣。可以證明SL和SR具有相同的特征矩陣Us，且λL與λR之和為單位矩陣。因此，當SL的特征值最大時，SR的特征值最小，可最大限度地區分兩類信號。則可構造出空間濾波器：

步驟2 特征提取。

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

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

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

2.4 運動想象腦電特征向量構造

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

其中：μk是第k 個特征的平均值組成的向量，σk是第k 個特征的標準差。

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

3 線下實驗

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

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

表1 平均準確率Tab.1 Average accuracy

圖4 單次實驗準確率對比Fig.4 Comparison of single test accuracy

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

3.2 多種識別方法結果比較

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

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

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

4 線上實驗

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

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

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

圖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 方差向量，從多個角度實現了對MI腦電信號特征的綜合描述，因此耗時增加。三種方法的耗時主要集中在EMD 過程，提取IES、Hilbert 邊際能量譜以及CSP 方差向量的耗時占比并不大。因此，從圖中可以看出，雖然HCHT 耗時比HHT和CSP略多，但三者的耗時差異不大，并不影響整體系統響應，且HCHT 算法具有更高的識別率和穩定性，更適用于實際應用環境。

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

5 結語

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