多特征融合的运动想象脑电特征提取方法

罗 飞，刘鹏飞，罗 元，朱思蒙

（重庆邮电大学光电工程学院，重庆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技术提供了一种新的思路，但识别率和识别种类有待进一步的优化。