运动准备电位单次检测技术研究

逯 鹏， 牛 新， 刘素杰， 胡玉霞， 胡航航(郑州大学 电气工程学院，河南 郑州 450001)

0 引言

目前，提高BCI系统实时性的主要途径是提高算法的效率[1].脑科学研究发现，大脑在自主运动前的1 500～2 000 ms已经被激活[2],在大脑额叶皮层区和初级运动皮层区诱发电位，称为运动准备电位(readiness potentials,RP)[3].RP由两个部分组成：早期的RP于自主运动前1 500 ms诱发一个缓慢上升的负电位；晚期的RP于自主运动前400 ms诱发一个陡峭的斜坡电位，如图1所示[4-5].RP波幅范围是2～10 μF，信噪比低、幅值小、易受干扰.此外，提取RP的过程易受肌肉收缩力度、收缩速度等影响[6].

图1 运动相关电位Fig.1 Motor relative potential

RP检测技术主要有Lew等[7]采用运动前的脑电信号对自主右手腕扩展动作进行了预测，结果表明在运动起始前的0.62 s对动作的预测正确率平均达到(62.75±10)%；Xu等[8]采用MRCP以位置保留映射方法提取EEG，采用LDA分类器进行分类，获得真阳性为(79±11)%的检测率.

笔者针对RP单次特征提取困难的问题，基于EMG信号确定运动准备的有效时间段，基于有效时间段内的脑地形图序列，定位激活脑区并选取关联通道，采用小波包变换与功率谱分析相融合的方法分析低频(0.1～1.0 Hz)信号，以信号的能量、均值和方差为特征，并采用支持向量机进行运动单次RP检测.

1 试验与数据采集

提示信息采用声音、图片等会诱发相关的脑电信号，会对RP信号造成污染[9].本试验设计一种受试者手腕自主运动的试验范式，避免提示信息产生的噪声污染.试验范式流程如下：受试者坐于计算机屏幕前，两手自然地垂放在椅子上，受试者注视计算机屏幕，开始提示音后，受试者手腕上下摆动5次，手腕自主摆动的时间间隔约为2 000 ms.2 min后，会有结束提示音，结束这一组试验，如图2所示.

图2 实施范式流程图Fig.2 Experimental flowchart

试验采用Neuroscan SynAmps 64导联脑电采集系统，10～20国际标准电极分布，采样频率设置为1 000 Hz，放大器的带通滤波设置为0.01～100 Hz；试验中，脑电采用单极记录方式，将右侧乳突(M2) 设置为参考电极，同时记录水平眼电和垂直眼电,电极阻抗小于5 kΩ.

为了单次检测RP，将信号分为两部分：无运动意图信号和有运动意图信号，如图3所示.把提示音前的1 000 ms作为0类；把运动开始前的1 000 ms作为1类，用于RP检测.

图3 用于训练分类器的数据段Fig.3 Selected EEG samples to build the training set of the movement intention intention classifier

2 方法

2.1 预处理

脑电信号中包含了大量的噪声干扰和伪迹信号，其中RP易受眼电伪迹、肌电伪迹干扰.采用二阶统计量和时序结构特征法去除眼电伪迹、肌电伪迹，以降低其影响[10].

2.2 确定有效时间段

运动准备的有效时间段指的是运动起始前的1 000 ms,如何确定运动的起始时刻是技术难点.

对EMG信号采用高通滤波器(高通5 Hz，低通50 Hz，24 dB/octave)滤波，EMG信号在运动时刻的能量有大幅的增加，选取这一点作为运动的起始时刻.基于EMG信号的能量在运动起始时刻突增这一现象，对受试者分别做9组试验，确定该受试者的最优的阈值，达到阈值的这一时刻，作为运动的起始时刻.其中一个受试者的EMG信号图形如图4所示，红色为最优阈值.

图4 阈值确定运动起始时刻Fig.4 Threshold determines the movement starting time

2.3 定位激活脑区并选取关键通道

大脑运动准备过程的大脑机理的研究尚不明确，导致与运动准备电位相关的通道也不能确定[11].

针对这一问题，我们采集了被试的全通道EEG，平均叠加EEG中的RP成分后，绘制出整个过程的脑地形图，分析时序关系和各个时刻激活的相关区域.其中一个受试者的脑地形图序列如图5所示.

图5 脑地形图序列(低频信号)Fig.5 The sequence of brain topographic maps(low frequency)

在图5中，运动的起始时刻前1 000 ms为零时刻，从0 ms时刻到运动起始时刻950 ms绘制脑地形图，每隔50 ms画一幅.在0～100 ms期间，大脑处于静息的状态；在550～800 ms期间，内侧前额及额中央部观察到ERP的差异，这些活动反映运动意图的预测过程.

从脑地形图序列分析发现C1、C2、C3、C4、C5、C6、CP1、CP2、CP3、CP4、CP5、CP6、FC1、FC2、FC3、FC4、FC5、FC6、P5、P3这20个通道在运动准备时，有较大的振幅.20个通道所在脑区是后顶叶和额叶皮质区，该脑区与运动规划和运动意图关联较大[12-13]，因此我们选择此20个通道提取EEG信号.

2.4 确定敏感频率点和频率段

在运动准备过程的有效时间段内，大脑的感觉运动皮层会出现EEG的节律性活动，包括低频直流成分[14-15].基于小波包变换的性质，采用多次迭代的小波转换进一步分析输入EEG的低频直流成分(0.1～1 Hz).

首先对采样频率为1 000 Hz，运动准备相关的20个通道，运动开始时刻前的EEG数据采用小波包进行12层分解，将0.1～1 Hz频率段EEG划分为7个频率段，如表1所示.

表1 频带范围Tab.1 The frequency range

根据7个频率段分量，对每一个通道均进行叠加平均，绘制各频段的功率谱，其中主要通道的功率谱图如图6所示.横坐标表示时间，纵坐标表示表中的序号，运动的起始时刻为1 000 ms.在序号2所对应的0.24～0.36 Hz频率段，运动起始前的1 000 ms的能量比其他时间段的能量高，认为此为运动准备过程的有效特征.

图6 0.1～1 Hz功率谱图Fig.6 0.1～1 Hz power spectrum

图7 不同受试者的相关通道功率谱图Fig.7 The power spectrum of Correlation channel of different subjects

2.5 通道优化

采用2.3节中选取的20个通道采取EEG数据，提取0.24～0.36 Hz频率分量作进一步分析，绘制了其中6名受试者能量随时间变化的规律.如图7所示，横坐标表示时间，纵坐标表示通道，运动的起始时刻为1 000 ms.不同的被试者之间与运动准备相关的通道有差异，但是都包含在选择的20个通道中.因此，选择20个通道的数据进行特征提取.

2.6 特征提取

基于2.2、2.3、2.4、2.5节中对于RP时域、频域特性分析结果，提出针对RP单次特征提取算法.

(1)对EMG信号进行处理，确定运动的起始时刻.

(2)输入已选择的20个通道，运动的起始时刻前1 000 ms的 EEG数据X(i).

(3)去除EEG的基线，采用二阶统计量和时序结构特征法去除眼电伪迹、肌电伪迹.

(4)选取‘db4′小波包基函数将EEG数据X(i)进行12层分解.

(5)重构低频0.24～0.36 Hz频率分量S12,3，其中S12,3是小波重构第12层第3节频率分量.

(6)计算S12,3对应的能量Ew，均值Mw和方差Vw.

(7)对每个通道重复4～6步骤，获得60维特征向量.

3 结果分析

3.1 训练模型

SVM 的核参数直接关系到分类器识别正确率的高低.经验值和试算方法确定核参数的计算效率低，且难以找到最优解[16].笔者采用网格搜索法获得最佳参数c和g，得到最佳模型.基于RP单次特征提取算法，选取每名受试者80%数据训练SVM模型，如图8所示.

图8 参数c和g的搜索结果Fig.8 The search results of parameter c and g

3.2 分类识别

分别对15名受试者进行9组试验，并对试验数据进行处理，剔除有明显错误、有明显波动等无效数据，使无运动意图类和有运动意图类两类样本比例为1∶1，取每名受试者80%数据训练，剩余20%数据进行测试，基于3.1节中所述的RP单次特征提取算法，采用SVM分类器分别进行单次RP检测，其结果如表2所示.不同的受试者单次RP检测准确率存在差异.15名受试者9组试验中最高检测率是77.5%～91.3%；每名受试者的9组平均检测率是68.2%～ 91.2%.

表2 单次RP检测正确率Tab.2 The single test results of RP %

4 结论

通过实时的BCI系统采集运动执行前的EEG和EMG信号，根据EMG信号确定运动准备的有效时间段，由脑地形图定位激活脑区，采用小波包变换的方法提取低频0.24～0.36 Hz频率分量，分析其幅值和能量的特征，用支持向量(SVM)分类识别.通过单次检测RP，以确定受试者的运动意图，为异步BCI系统的实现提供了新的技术途径.笔者通过离线分析只对运动前1 000 ms的数据研究RP，不同受试者检测准确率存在差异.今后将进一步研究克服不同受试者的差异，提高预测RP的准确率，将其应用于在线系统中.