陆振宇，陆旭峰，杨瑞洪，常珊

（1.江苏理工学院电气信息工程学院生物信息与医药工程研究所，常州213001；2.扬州工业职业技术学院，扬州225127）

0 引言

脑机接口（Brain Computer Interface，BCI）是综合神经科学、信号处理、机器学习和信息技术等领域的跨学科技术。BCI 的目的是把大脑活动转变成对设备的控制指令，具有恢复失去的感官和运动功能的潜力，如治疗帕金森病的脑深度电刺激，截肢者和脊髓损伤病人的假肢、轮椅，还有中风和肌肉萎缩患者使用的鼠标和单词拼写设备[1]。BCI 已经在游戏娱乐、机器人、生物识别和教育等应用领域表现出强大的发展前景。

目前，运用较为广泛的信号是当受试者自主想象某种特定动作所产生的神经活动，称之为基于运动想象（Motor Imagery，MI）的脑电信号。研究表明，当受试者进行手、脚、舌头动作的运动想象，大脑皮层相对应区域的μ频带（8-13Hz）表现显著，而对侧区域其节律能量下降，该现象被称为事件相关同步/去同步（ERD/ERS）现象，这也为信号的特征提取提供了依据[2-4]。

BCI 的实现主要涉及信号处理和模式识别。信号处理包括数据的预处理，如滤波、伪迹去除，还有关键的特征提取技术。文献中的特征提取方法一般有，频域分析有自回归（AR）模型、功率谱估计；时频分析有小波变换、小波包变换、局部均值分解；空间域有主成分分析（PCA）、独立分量分析（ICA）和共空间模式（CSP），用于多导联的脑电信号[5-11]。模式识别是利用机器学习技术，对输入特征进行分类。常用的二分类技术有线性判别分析（LDA）、神经网络（NN）、支持向量机（SVM）、基于随机森林和Boosting 的集成技术。

通过对相关文献的研究，信号变换和分解的方法多样，一般在信号处理后选取能量、熵等作为特征，但变换方法及提取的特征单一化，无法多方面地表示脑电信号。本文选取提取算法易于理解、计算量小、效果良好的AR 模型、小波包变换、CSP 方法对信号进行处理，提取信号的多方面特征，采用SVM、LDA 分类器进行模式识别，并对分类结果评价。

1 实验数据

本文实验的数据是来源于2003 年BCI 竞赛的Graz 数据集，由奥地利Graz 科技大学医学系、生物医学工程研究所提供[12]。实验范式是受试者放松地坐在手扶椅上，根据屏幕显示做出相应的反馈。前2s 为准备时间，在第2s 时，蜂鸣器提示实验开始，同时屏幕出现1s 的“+”符号，3-9s 是受试者按给出的箭头方向进行对应的运动想象。实验采集了C3、CZ、C4 通道的脑电信号，采样频率为128Hz，并通过0.5-30Hz 的带通滤波。整个数据集共有280 个样本，其中左右手运动想象各140，训练和测试样本各占一半，数据集分布均匀。

2 特征提取算法

2.1 AR模型

AR 模型可用如下公式表示：

式中，ε（n）是均值为零、方差为σ2的白噪声，p 为模型阶次，ap（i）为自回归模型参数[4-5]。本文选取竞赛数据C3、C4 通道，时间段在4-7s 的信号，基于最终预测误差准则，确定阶次为5，使用Burg 算法提取其AR 模型系数和系数差作为特征。

2.2 小波包变换

小波包变换是为了克服小波分解在高频段的频率分辨率较差，而在低频段表现较好的基础上提出的。它是一种更精细的信号分析的方法，提高了信号的时域分辨率[13]。

竞赛数据采样频率为128Hz，依据Nyquist 采样定理，脑电信号频率范围为1～64Hz，使用小波基函数对C3、C4 信号进行4 层小波包分解，分解频宽为4Hz。分解后的第4 层小波包节点频率范围如表1 所示，具有代表性的μ节律信号为节点（4,3）。图2 为采用bior3.3 小波基，对单个样本的C3、C4 通道信号做小波包变换得到的时域和频域图。

表1 第4 层小波包分解的节点及频带范围

图2 单样本节点（4,3）的bior3.3小波包变换

脑电信号经小波包分解后，节点（i,j）信号的小波包能量由计算该节点的小波包系数平方和得到。

式中，i 为分解层数，j 为该层下的节点，di,j（n）为节点（i,j）的小波包系数。

根据能量守恒定律，分解后的各节点能量总和等于原信号的能量，则相对小波包能量为节点信号能量与该层节点总能量的比值：

式中，Eraw为原信号能量，Etotal为原信号分解i 层后各节点能量总和，pi,j则为节点（i,j）能量的归一化，称为相对能量。

信息熵可以表征信号的不确定度，结合小波包变换，可得出分解后的节点信号小波包熵为：

本文选取C3、C4 通道信号做小波包变换，以一次实验所记录的C3 信号为例，具体步骤如下：

（1）采用固定时间窗的方法，以128 点的时间窗，将长度为1152 点的信号分割成9 段；

（2）在一段128 点的信号上做小波包分解，小波基选择bior3.3[4]，分解层数为4，选择节点（4,3）的重构信号，由式2 计算出小波包能量，式4 得出小波包相对能量，式5 得出小波包熵。

（3）C4 通道重复1、2 步骤，一个样本的小波包变换结果为2×9 的矩阵，按样本类别进行叠加平均，得到C3、C4 通道脑电信号的小波包能量、相对能量和熵的变化，图3 为右手运动想象下的平均小波包能量和熵走势图。

图3 右手运动想象的bior3.3平均小波包变换特征

由图3 可知，4-7s 的C3、C4 通道小波包相对能量和小波包熵区分度明显，走势基本一致，数值大小不一样。本文将其二维向量特征转化成一维向量，一个样本C3、C4 通道的小波包能量差、相对能量差、熵差△WPE 作为分类的特征，图4 为小波包熵差的走势图。

图4 左右手运动想象的bior3.3平均小波包熵差

2.3 共空间模式

CSP 是寻找空间滤波器，使得滤波处理后的数据与其中一类方差达到最大，而与另一类方差达到最小的算法。基于运动想象的BCI，依赖特定的频带能量，对于两类的左右手运动想象信号，经过CSP 处理得到的特征向量增强了两类的区别，可区分程度最大化[14]。

单次脑电信号样本可由Ns×Nc的二维矩阵E 表示，其中，Ns是采样点数，Nc是通道数，则归一化协方差矩阵为：

式中ET表示E 的转置，表示矩阵的秩。按样本类别求叠加平均的协方差矩阵，然后总样本的协方差矩阵为：

式中Cov1、Cov2 分别表示类别1、2 的协方差矩阵。对总样本协方差矩阵作特征值分解，得到特征向量矩阵U、特征值对角阵λ：

将U、λ进行降序排列，求得白化矩阵：

对类别1 的协方差矩阵使用白化矩阵进行转变：

对其作特征值分解，得到特征向量矩阵U1、特征值对角阵λ1：

同样对U1、λ1进行降序排列，得到类别1 的投影矩阵W：

取W 的前m 行、后m 行，构成空间滤波器W’，原信号E 经过空间滤波器得到：

对E’按行求方差，并取对数，得到特征值f：

一般地还有，在方差的归一化后取对数，得到特征值f：

本文选取4-7s 的运动想象脑电信号，对其做CSP变换，按式（15）、（16）分别计算，各自得到2 个特征值，用于区分左右手运动想象。

图5 左右运动想象各70个样本的特征值散点图

3 模式判别

3.1 线性判别分析

线性判别分析（LDA）是线性二分类器，将输入向量映射到一个超平面，该超平面将输入空间划分为两个半空间：每个空间代表一种类别。决策边界由超平面公式决定：

式中w 为法向量，w0为阈值，均由训练数据得出。对于新的输入特征向量，代入式16，如果g（x）≥0，为类别1，反之为类别2。

3.2 支持向量机

支持向量机（SVM）是寻找两类样本间隔最大的分离超平面，以获得最好的泛化能力的一种分类器。本文采用libsvm 工具，选择C-SVM 类型，利用径向基核函数实现脑电信号特征的非线性映射，使用网格搜索获得最优参数C、gama。

4 实验结果分析

为了验证本文提出方法的有效性，将单特征和多特征融合提取方法运用在整个数据集上，采用SVM、LDA 分类器对其分类，结果如表2、3，其中识别率是10折交叉验证的平均值。

通过表2 可以直观地发现：单特征选取方面，AR系数和系数差都能较好地被区分，识别率最高为87.5%；小波包熵差的效果高于能量差、相对能量差，原因在于小波包熵的数量级比后两者大，从图3 可知，小波包能量的范围在0.25～0.75，而小波包熵是1.2～3.0；而共空间模式下的两种特征在分类器上表现了较大差异，未归一化的特征在LDA 上的表现最好，原因在于数据的数值较大，而对于SVM，无论特征是否归一化，准确率都是73.93%，表现一般。

通过表2 的分析，选取分类效果好的AR 模型系数、小波包熵差、未归一化的CSP 特征，进行排列组合，用于多特征融合，得到表3 的结果。其中，分类器方面，LDA 的识别率均比SVM 略高；特征融合方面，三者组合的准确率最高，达到了91.43%，而两两组合的准确率比各自的单特征高，均在90%左右。

表2 单特征的平均识别率（单位：%）

表3 多特征融合的分类准确率（单位：%）

5 结语

本文提出了多方面对运动想象脑电信号进行特征提取及融合，并使用分类器对其评价。AR 模型代表了脑电信号的功率谱估计，仅使用C3、C4 两通道信号的系数，也可以得到较好的分类结果；小波包变换是定位到给定频带，计算熵值可以反映其能量的差异性，选择C3、C4 通道信号的熵差作为特征，也得到了较好的效果；CSP 算法是提取脑电信号的空间特征，适用于较多导联，且需要训练数据构造滤波器，本文的竞赛数据是3 导联，使用CSP 仅得到两个特征，分类效果不稳定。而多特征融合的方式，则弥补了单特征的缺陷，在分类效果上也更好，这为运动想象脑电信号的特征工程提供了思路。在分类器的选择上，SVM 虽然泛化能力高，但计算量大，对数据敏感，而LDA 实现简单，运算速度快，效果也较好，适用于在线BCI。