面向人机交互的运动想象脑电信号感知算法研究

王 康，翟弟华，夏元清

（北京理工大学自动化学院，北京100081）

1 引 言

脑机接口（Brain Computer Interface，BCI）可以在人脑与外界的环境设备之间建立直接交流的通道，从而实现大脑对外部设备的交互控制［1］。因可以在控制系统中融合人与机器的智能，脑机接口技术在军事、医疗、娱乐等领域中展现了广阔的前景。21 世纪初美国就在意念控制机器人士兵项目中投入了大量的研究经费。2004年，美军资助了多个研究机构开展“思维控制机器人”项目的研究。2009年，多名美国著名学者撰写了《Opportunities in Neuroscience for Future Army Application》［2］，将脑机接口列为军事领域重点发展技术。2013年，美国国防部提出“阿凡达”项目，拟实现意念远程操控“机器战士”，代替士兵在战场上执行作战任务［3］。近年，我国在此领域也加大研发，并取得了不错的进展。如国防科技大学研制成功脑控汽车，可以根据人脑的思维意识完成启动、转弯和加速等复杂操作［4］。可以相信，脑机接口技术可为机器人、装甲车等武器装备进行智能化赋能，引发武器装备在操控方面的革命。

脑机接口系统主要包括脑电信号的采集与预处理、特征提取、分类识别和输出控制。目前，脑机接口的输入信号主要采用统计特征较为明显的脑电信号，其中运动想象脑电信号通过检测大脑μ 节律和β 波的能量对人脑的意图进行判断［5］，被广泛应用于脑机接口系统。特征提取与分类识别是实现脑机接口系统的核心技术［6］，识别结果可以转换成可被计算机控制指令。运动想象脑电信号常用的特征提取算法包括功率谱估计、自适应自回归模型、小波变换和共空间模式等［7-9］，常用的分类识别算法包括线性判别分析、隐马尔科夫模型、支持向量机和人工神经网络等［10-11］。目前，运动想象脑电信号在二分类任务可以取得良好的结果，基于共空间模式的特征提取算法在右手和舌头的实验中可以取得高达97%的正确率［12］。而在多分类任务中，分类结果表现并不理想，李昕等［13］利用小波变换和支持向量机在多任务的分类识别中最优分类准确率仅为69%；施锦河等［14］利用共空间模式、Hilbert算法和支持向量机在多分类任务中的最优准确率为91%，但实验使用了60 个脑电电极，准确率在实验过程中存在波动，稳定性不足；天津大学的万柏坤等［15］利用了二维的时频域分析、Fisher 判别法和支持向量机进行了多分类任务的研究，最优准确率达到了86%，但使用了60个电极，操作较为复杂，缺乏实用性。

考虑到单一算法的局限性，本文采用了小波包分解算法和一对多共空间模式对脑电信号进行联合特征提取，将原始信号分解为多个子频带，提取了细化的时频域特征，充分考虑了运动想象脑电信号的相关作用频带，并使用人工神经网络对提取的特征进行分类。此外，尽管一些文献提出的算法展现出了良好的效果，但由于受到测试对象和环境等因素的影响，分析的脑电数据存在差异，缺乏统一性，难以得出较为客观的比较结果。因此，本文使用“BCI IV”竞赛中提供的公共标准数据Dataset 2a对算法进行了验证，并与竞赛前几名的数据进行比较，结果表明，本文算法的识别准确率有着比较明显的提升。

2 算法原理

2.1 小波包分解

小波包分解［16］（Wavelet Packet Decomposition，WPD）是在小波变换（Wavelet Transformation，WT）的基础上优化得到的信号分解与重构方法。小波变换可以对信号进行局部化与多尺度分析，在信号分析中具备一定的自适应性。但是，小波变换只是对信号的低频部分进行分解与重构，而小波包变换考虑了信号的高频部分，能够有效降低信号冗余，克服了小波变换的固有缺点。因此，小波包变换对非平稳的神经生物信号、机械振动信号等可以进行更好的时、频域分析。

小波包分解使用低通滤波器hn(n ∈Z)和高通滤波器gn(n ∈Z)对原始的非平稳信号进行滤波，其中gn=(-1)nh1-n，二者为正交关系。正交尺度函数φ(t)和小波函数φ(t)满足：

考虑到公式表述的一般性，本文将正交尺度函数φ(t)和小波函数φ(t)均表示为μ(t)。在尺度0上，尺度函数为μ0，0(t)，尺度1 上的正交尺度函数和小波函数分别为μ1，0(t)和μ1，1(t)，则式（1）在尺度2上可表示为：

对于尺度j，递推可得到表达式为：

其中，函数系μj，n表示小波函数φ(t)的小波包，图1为三层小波包分解的示意图。

2.2 共空间模式

图1 小波包分解图Fig.1 Wavelet packet decomposition diagram

共空间模式［17］（Common Spatial Patterns，CSP）是一种典型的运动想象脑电信号特征提取算法，在二分类任务中取得了良好的效果。CSP算法的原理是通过设计空间滤波器，使一类信号的方差最大化，同时使另一类信号的方差最小化，提取区分度较高的脑电特征。在多分类问题中，可以使用一对多的思想将原始的二分类CSP 算法进行扩展，具体可以将多类信号中的一类作为单独一类，其余类别统一作为另外一类处理，依次计算得到四个二类投影矩阵，构造一个四类的CSP特征提取模型。

假设Xi是大小为N×T 的第i 类原始脑电信号，其中N 表示脑电信号的通道数，T 表示采样点数，Xi的归一化协方差矩阵为：

四分类任务中，混合空间协方差矩阵可以计算为：

其中：R1、R2、R3、R4分别表示对应类别运动想象信号经过多次实验中得到的平均协方差矩阵。R经过特征分解可以得到：

其中：U0表示R 的特征向量，V 表示R 的特征值矩阵。对特征值矩阵V进行降序排列，并对U0进行相应调整，白化矩阵可计算为：

对于四分类任务，一对多CSP 算法在计算脑电信号的投影矩阵时，将其中一类信号作为基准，其余的三类信号则归为另一类，以R1基准的计算为例，另一类R′1可以表示为：

分别计算R1和R′1的平均协方差矩阵，并表示为，然后对二者进行白化变换可得：

分别对S1和S′1进行特征分解：

其中：V1和V′1满足V1+ V′1= I。投影矩阵可以计算为：

同理，Ri对应的投影矩阵可以计算为：

将脑电信号Xi进行投影，可以得到：

对各类别投影后的信号Zi(i = 1，2，3，4)各 取前m 行可以组成一个新的四类CSP 特征提取模型，共有N=4m行，特征向量的计算公式为：

2.3 人工神经网络

人工神经网络［18］由神经元模型组成。神经元模型是一种非线性的信息处理节点［19］，结构如图2所示。

图2 神经元模型结构图Fig.2 Structure diagram of neuron model

其中，x1，…，xn为神经元i 的输入，wi1，..，win为连接权重，yi为神经元的输出，f（·）为激活函数，则模型的输出为：

神经网络有着多种训练方法，其中最为常用的是反向传播算法（back propagation，BP）［20］。

图3 为基于反向传播的单隐层神经网络结构，其中X1，…Xm为神经网络的输入，y1，…，yL为神经网络的输出。

基于BP的神经网络算法属于监督学习算法，学习过程分为正向传播和反向传播两部分。在正向传播过程中，神经网络输入信号经过网络逐层传输到输出层，如果结果符合样本的期望输出，则算法结束；否则，需要计算网络输出与样本的期望输出之间的误差，借助梯度下降法沿网络的反方向对网络各层的权重和阈值进行调整，达到最小化误差的目的［21］。

图3 BP神经网络图Fig.3 BP neural network diagram

3 模型构造方法

针对四分类任务，首先利用一对多CSP 算法构建空间滤波器，使不同类别间的方差最大化，经过投影可得到区分度较高的特征向量。由于运动想象脑电信号的主要信息位于μ节律（8～12Hz）和β波（14～30Hz），原始的脑电信号无法准确提供对应频带的信息，本文首先使用WPD 算法分解原始信号，得到多个细化频带信息，然后使用CSP 算法提取各个频带的脑电特征，最后使用人工神经网络进行分类，算法流程如图4所示。

具体步骤为：

图4 算法流程图Fig.4 Algorithm flow paradigm

步骤1：截取3～6s时间段的脑电信号，并将其滤波至8～30Hz（所选数据集的运动想象任务主要集中于3～6s，μ 节律和β 波频率范围主要在8～30Hz 频带内）；

步骤2：使用WPD 对训练集样本进行分解与重构，得到含有不同频带信息的信号样本；

步骤3：对分解后的样本使用一对多CSP 算法进行特征提取，作为神经网络的输入；

步骤4：将步骤3得到的特征向量输入神经网络进行训练，修正网络的权重和阈值，完成训练；

步骤5：按照步骤2、3 对测试集样本进行处理，然后输入训练好的神经网络得到最终的分类结果。

4 实验结果及分析

4.1 实验数据

本文的实验数据选自“BCI IV”，数据集由Graz University of Technology 提供，共包含9 个志愿者的四分类运动想象脑电信号，四个任务分别对应左手、右手、双脚和舌头。每个志愿者在不同的日期共完成两轮实验，每轮实验采集288组实验数据，由四类运动想象任务等分，每组实验数据的采集过程按照图5的实验范式所示。

图5 实验范式流程图Fig.5 Flow chart of experimental paradigm

用于脑电信号采集的实验设备共有22个电极，摆放位置如图6 所示。设备的采样频率为250Hz，使用0.5～100Hz的带通滤波器对采集的原始信号进行滤波。两轮实验数据分别作为训练样本与测试样本，前者用于模型训练，后者用于模型验证。

4.2 参数选择

本文利用WPD 和CSP 算法联合提取脑电信号特征，需要确定CSP 每类模式的特征数m。在使用人工神经网络进行分类的时候，由于神经网络的参数决定了网络结构的不同，需要确定隐藏层的神经元个数、目标函数的正则化系数等参数。不同参数的选择会影响最终的实验结果。

图6 电极摆放位置图(遵循10～20规则)Fig.6 Electrode position diagram（follow 10～20 rule）

运动想象脑电信号的主要信息位于μ 节律（4～12Hz）和β 波（14～30Hz）中，考虑到两个节律的频带范围和频率分辨率，本文选取“db4”小波作为小波基函数，进行5层小波包分解，选择位于4～32Hz的8个相关子频带进行后续的特征选择。

CSP模式特征数m一般选择在2～12之间，为确定最优的特征数m，将神经网路结构暂定为双隐层，各隐层的神经元个数暂定为100，正则化系数暂定为0.1，调节m，则Kappa值变化如图7所示。

图7 不同m值的Kappa值Fig.7 Kappa rates with various value of m

当CSP模式的特征数m为6时Kappa值最优，固定该参数，然后对神经网络结构进行调整，需要的调整包括隐层数量、隐层的神经元数量、正则化系数以及dropout层的使用。

图8 不同隐层数量的Kappa值Fig.8 Kappa rates with various number of hidden layers

暂定隐层神经元数量为100，正则化系数为0.1，无dropout 层，实验结果如图8所示，可以看出，隐层数目为1时，具有最优的Kappa值0.72。

将隐层数目固定为1，隐层的神经元数量为100，正则化系数为0.1，设置dropout 层的参数调节范围为0～0.5，实验结果如图9所示。

图9 不同dropout参数的分类准确率Fig.9 Kappa rates with various parameter of dropout

当dropout参数为0.4时，最优Kappa值为0.73，固定该参数，暂定正则化系数为0.1，隐层神经元数目的调节范围设置为50～250，实验结果如图10所示。

隐层的神经元数目为100时，有着最优的Kappa值0.73。固定其他参数，调节正则化系数，当正则化系数为0.1时，Kappa 值为0.73；当正则化系数为0.01时，Kappa值为0.70。

4.3 实验结果

图10 不同隐层神经元数目的Kappa值Fig.10 Kappa rates with various number of neurons in hidden layer

选择“db4”小波基函数进行5层小波包分解，设置CSP 模式的特征数为6，选用单隐层神经网络作为分类器，隐层神经元数目、正则化系数和dropout参数分别选择100、0.1、0.4，在基准数据集中的实验结果如表1所示。

表1 实验结果对比Table 1 Comparison of experimental results

可以看出，本文提出的算法与竞赛前3 名的算法准确率相比在整体上有着较为明显的提升，证实了本文算法的有效性。

5 结 论

针对运动想象脑电信号在多分类问题方面特征区分度不高、识别精度低的问题，本文使用了小波包变换对原始脑电信号进行了多层分解，利用一对多CSP 算法对分解后的子频带进行特征提取，利用单隐层神经网络进行分类。最终的实验结果表明，选择合理的小波包分解层数、CSP的m参数和神经网络结构，可以明显提高分类的准确率，与BCI IV竞赛的基准数据集前几名结果的对比中，展现了本文方法的优越性。本文的算法对于实现基于运动想象脑电信号的脑机接口控制具有重要的意义，未来将对在线系统的算法部署进行更为深入的研究。