基于VMD与CSP的脑电特征提取方法

刘 帅，乌日￣开西·艾依提

(新疆大学机械工程学院，新疆 乌鲁木齐 830047)

1 引言

脑机接口(Brain Computer Interface， BCI)为人类与周围环境进行交互提供了一种不依赖于外周神经和肌肉的新的输出通道[1，2]。脑电图(Electroencephalographic，EEG)具有无创、高时间分辨率、低成本等优点，被广泛应用于脑机接口中。研究表明，当人在想象左手或右手运动时，大脑的感觉运动皮层Mu节律和Beta节律能量会产生变化，这一现象被称为事件相关去同步(Event-related Desynchronization， ERD)和事件相关同步(Event-related Synchronization， ERS)[3]。近年来，运动想象脑机接口技术被用来训练神经肌肉损伤的患者，帮助其恢复肢体的运动感知功能[4，5]。由于运动想象脑电信号是由大脑活动产生的，信噪比低，容易受到外界噪声的影响，难以获得较高的识别精度。如何对脑电信号进行有效的特征提取和分类成为一个难点问题。

目前，国内外研究者提出了若干运动想象脑电特征提取和分类算法，常用的特征提取算法有自回归模型(Auto Regressive，AR)、小波包变换(Wavelet Package Transform，WPD)、功率谱密度(Power Spectral Density，PSD)、共空间模式(Common Spatial Pattern，CSP)等[6-8]。常用的分类算法有支持向量机(Support Vector Machine，SVM)、线性判别分析(Linear Discriminant Analysis，LDA)等[9]。杨默涵等[10]提出了一种基于总体经验模态分解(Ensemble Empirical Mode Decomposition，EEMD)的运动想象脑电处理方法，将分解后的分量利用希尔伯特变换提取边际谱及瞬时能谱，选取较优频带波幅特征作为特征，取得了较好的分类效果。Shalu Chaudhary等[11]采用两种时频分析方法将脑电信号转化为图像，输入到基于AlexNet模型的深度卷积神经网络(Deep Convolutional Neural Networks，DCNN)进行特征提取与分类，发现连续小波变换(Continuous Wavelet Transform， CWT)转化的图像具有更高的识别率。Yu Z等[12]提出了一种时间约束稀疏组空间模式(Temporally Constrained Sparse Group Spatial Patterns，TSGSP)特征提取算法，通过同时优化CSP中的滤波器频带和时间窗口提高运动想象脑电信号的分类精度。

为了进一步提高在少数通道下的运动想象脑电信号的分类准确率，本文提出一种变分模态分解(Variational Mode Decomposition，VMD)结合CSP的特征提取方法。将脑电信号进行VMD分解，得到多个具有物理意义的固有模态函数(Intrinsic Mode Function，IMF)，通过巴氏距离计算各阶固有模态函数与原始信号之间的相关程度，选取合适的固有模态函数组合成为新的信号矩阵作为输入，再通过CSP构造空间滤波器，得到各IMF的空间分布特征，最后使用支持向量机进行分类。

2 数据集描述

本文采用两个数据集，第一个是BCI Competition II的Data set III的数据集[13]，数据集记录了一位25岁的健康女性，实验要求被试坐在屏幕前，根据屏幕上随机出现的左右箭头提示进行相应的左右手运动想象。数据采用双导联方式记录，采集了C3，Cz和C4三个通道的EEG信号。采样频率为128Hz，对信号进行0.5-30Hz滤波。实验共记录7组，每组包含40次实验，每次实验持续9s，共280次实验，选其中的140次实验作为训练样本，剩余的140次实验作为测试样本。单次实验的采集范式如图1所示。

图1 Data set III单次实验采集范式

第二个数据集来自三位身体健康的志愿者，均为右利手，且均为第一次参加脑电实验。使用Neuracle公司的NeuSenW16通道无线脑电采集系统进行实验数据的采集，脑电帽的电极按照国际标准10-20系统放置，仅采集了C3、C4两个通道的数据。参考电极位于头顶部，接地电极位于前额。采样频率设置为1000Hz，初始带通滤波为0.5-100Hz，陷波滤波50Hz。

在实验过程中，受试者佩戴脑电帽被要求坐在距离电脑屏幕大约1米的椅子上，双手自然放置在桌子上，保持放松状态。在想象过程中，要求被试尽量避免眨眼、眼动、身体以及头部移动。单次实验数据采集范式如图2所示。试验开始前2s为准备阶段，屏幕显示空白；2s后有一段短暂的提示音并且屏幕显示“+”字形，提示被试实验马上开始；4-7s屏幕会随机出现向左或向右的箭头，被试需要根据提示进行左右手的运动想象；7s后有2-4s的随机休息时间；为了防止被试疲劳，每隔20次实验休息2min。实验共分为3组，每组100次实验，共300次实验，其中200次实验作为训练样本，100次实验作为测试样本。

图2 实验室单次实验数据采集范式

3 改进的特征提取方法

因个体差异较大，每位受试者在运动想象时，ERD/ERS发生的频率段不相同。传统的共空间模式是根据固定频率段的脑电信号进行特征提取，且需要大量通道的数据。针对这一问题，本文通过结合VMD和CSP用于改善特征提取算法。改进后的特征提取流程如图3所示。首先，将滤波后的信号进行VMD分解成K个IMF分量，然后通过计算各IMF分量与原始信号的巴氏距离，选择巴氏距离最小的3阶IMF分量进行重构，采用CSP滤波提取特征，得到特征向量F。

3.1 变分模态分解

变分模态分解通过迭代搜索变分模型的最优解来确定每个模态分量的中心频率和有限带宽，可以有效抑制经验模态分解(Empirical Mode Decomposition，EMD)方法中存在的端点效应和模态混叠现象，具有坚实的理论基础，且对采样和噪声的鲁棒性更强[14]。

图3 VMD与CSP相结合的特征提取流程

在VMD算法中，定义固有模态函数(IMF)为一个调幅-调频信号，通过混合一个预估中心频率，将每个模态的频谱调制到相应的基频带

(1)

计算梯度平方L2范数对信号解调，得到各个模态函数的带宽，则相约束变分问题如下：

(2)

式中，{uk}={u1，u2，…uk}为分解后的k个模态分量，{ωk}={ω1，ω2，…ωk}为每个模态对应的中心频率，f为输入的脑电信号。

使用二次惩罚因子α和Lagrange乘法算子λ，将约束变分问题转变成非约束变分问题，增广型的拉格朗日表达式如下

L({uk}，{ωk}，λ)

(3)

利用交替方向乘子法(Alternate Direction Method of Multipliers，ADMM)，并结合Parseval/Plancherel 傅里叶等距变换，优化得到各模态分量和中心频率，求解增广拉格朗日函数的最优解，交替寻优迭代un+1，ωn+1和λn+1的表达式如下：

(6)

迭代终止条件

(7)

最终得到K个有限带宽的IMF分量。

3.2 模态函数选择

利用巴氏距离用于测量IMF分量与原始信号的相似性。对于相同定义域X上的离散概率分布p和q，巴氏距离为

DB(p，q)=-ln(BC(p，q))

(8)

(9)

其中，p表示第K个IMF分量，q表示原始信号，BC(p，q)为离散概率分布 Bhattacharyya系数。

3.3 传统共空间模式算法

共空间模式是一种对两分类任务下的空域滤波特征提取算法，能够从多通道的脑电信号中提取每一类的空间分布成分。共空间模式的原理是利用矩阵对角化，找到一组最优空间滤波器对信号进行投影，使得一类方差最大化，另一类方差最小化，从而得到具有较高区分度的特征向量[15]。

单次运动想象脑电实验数据由一个N×M的矩阵X组成，其中N表示通道数，M表示采样点数，其归一化后的协方差矩阵计算如下

(10)

(11)

对混合协方差矩阵Vc进行特征分解

(12)

式中，Uc为特征向量矩阵，λc为特征值矩阵。构造白化矩阵P

(13)

(14)

式中，Sl和Sr具有共同的特征矢量E，且对应的特征值之和为1，则Sl最大特征值对应的特征向量是Sr最小特征值对应的特征向量。对Sl和Sr进行特征分解可得

Sl=EλlETSr=EλrET

(15)

式中，对角矩阵λl和λr之和为单位矩阵I。由此构造投影矩阵W

W=ETP

(16)

则对输入矩阵X进行滤波可得Z=WX，取Z的前m列和后m列构成特征矩阵Zp，p=1，2，…，2m。

(17)

3.4 支持向量机

支持向量机是一种有效的二分类模型，具有较强的泛化能力。其基本原理是将输入的向量映射到高维的特征空间，并在此空间寻找一个最优分类超平面，使两类样本的分类间隔最大[16]。本文将改进的CSP特征提取出来的特征向量作为输入，使用LIBSVM工具包进行分类。训练集用来训练模型，测试集用来评估最优超参数训练出来的模型，同时对模型的参数进行寻优。采用网格搜索算法寻找最佳惩罚系数c和核参数g组合。

4 实验结果分析

4.1 数据预处理

根据ERD/ERS现象，在做肢体运动想象时，感觉运动区域的节律变化主要发生在Mu节律(8-13Hz)和Beta节律(14-30Hz)，本文使用6阶Butterworth滤波器对原始EEG信号进行8-30Hz的滤波。

利用连续小波变换对信号进行时频分析，通过时频图观察运动想象期间ERD发生的时间段。想象左右手运动时的时频图与脑地形图如图4和图5所示。想象左手运动时，在8-12Hz频段对侧的C4通道3.5-7s时间段的能量显著降低，发生ERD现象；而同侧的C3通道3.5-7s时间段的能量显著增加，发生ERS现象。与之相反，在想象右手运动时，对侧的C3通道3.5-7s发生ERD现象，同侧的C4通道3.5-7s发生ERS现象。对于BCI Competition II的Data set III的数据集，本文选择单次实验的3.5-7s时间段的脑电信号进行后续的特征提取。数据集中每次实验共采集9s数据，采样频率128Hz，每次实验每个通道共有1152个采样点。选取其中3.5-7s的数据，此时数据为448×3×140。

图4 想象左手运动时频图与脑地形图

图5 想象右手运动时频图与脑地形图

4.2 特征提取与分类

将预处理后的脑电信号进行VMD自适应分解，得到K个IMF分量，每阶IMF分量都有各自的中心频率和有限带宽。然后计算各IMF分量与原始信号的巴氏距离来选择合适的IMF分量。本文K取值为7。以想象左手为例，C3通道的脑电信号经VMD分解后的IMF波形及各分量对应的频谱图如图6所示。从频谱图中可以看出，自适应分解后的IMF分量主要集中在8-30Hz范围，不同于EMD分解时出现的相同成分的信息出现在不同IMF分量中，VMD分解后的各阶分量均为有限带宽，有效克服了EMD中的模态混迭现象。

图6 VMD分解后的各阶IMF分量及各分量对应的频谱图

根据式(8)计算各IMF分量与原始信号的巴氏距离。以单次实验数据为例，表1为C3和C4通道经VMD分解后各IMF分量与原始信号的巴氏距离。巴氏距离越小，对应的模态分量与原始信号越相似。观察表1可知，C3通道巴氏距离最小的三阶IMF分量分别为IMF2、IMF3和IMF4；Cz通道巴氏距离最小的三阶IMF分量分别为IMF2、IMF1和IMF3；C4通道巴氏距离最小的三阶IMF分量分别为IMF2、IMF4和IMF7。

表1 IMF分量与原始信号的巴氏距离

本文选取每个通道数据分解后与原始信号巴氏距离最小的三阶IMF分量重新构造新的输入矩阵作为特征提取的输入，原输入矩阵重新构造成一个448×9的矩阵，再通过CSP滤波提取特征得到特征向量，最终将特征向量输入到支持向量机中进行分类。

4.3 实验对比

为了证明本文方法的有效性，在相同脑电数据集情况下，与其它文献的方法进行了比较，对比结果如表2所示。可以看出，在使用同一数据集的情况下，相比于其它三种算法，本文所提方法的分类结果取得了91.43%的分类准确率。说明本文所提出的特征提取算法能够有效提高运动想象脑电的分类准确率。

表2 本文方法与其它方法分类准确率对比

为了进一步验证该方法的有效性，在实测数据中进行了验证，与传统的特征提取算法CSP进行了比较。如表3所示，对实验室实测的三名被试的C3和C4通道的脑电信号分别进行CSP特征提取和VMD融合CSP的特征提取，并用支持向量机进行分类。从表3可以看出，相比于传统的CSP特征提取算法，经VMD分解后的信号再进行CSP滤波提取特征能够获得更高的分类识别率，三名被试的分类准确率分别提高了6%、8%、6%。

表3 实测数据的分类识别率比较

5 结论

针对运动想象脑电信号个体差异较大，传统的CSP特征提取算法需要大量通道的数据的问题，本文提出了一种VMD结合CSP的特征提取算法。通过巴氏距离可以选择与脑电信号最相关的固有模态函数构造信号矩阵进行特征提取，最后使用SVM分类。在公共数据集上取得了91.43%的分类准确率，在实测数据上相比于传统的共空间模式特征提取方法，取得了更好的分类效果。

实验结果表明该方法在少数通道的情况下能够有效提取运动想象脑电信号的特征，取得了较好的分类准确率，证明了该方法的有效性，为运动想象脑机接口应用奠定了基础。在下一阶段，本文提出的方法将用于实现在线脑机接口。