时频分析在语言想像脑机接口中的应用

王 力 张 雄 仲雪飞 樊兆雯 张 玉 孙 瀚

(东南大学电子科学与工程学院，南京210096)

脑机接口(BCI)技术是指通过检测大脑的活动以实现外部设备控制的一种技术.应用BCI 技术，可绕开外周神经和肌肉组织，直接为大脑创造新的信息交流和控制通路［1］.成本低廉、操作方便的脑电图(EEG)是BCI 系统的主要检测方式［2］.基于EEG 的脑机接口有如下的实验范式:稳态视觉诱发电位、P300 诱发电位、意识任务和运动想像等［3］.其中，前两者主要依赖于额外的设备.意识任务往往需要使用者想像复杂的任务，如想像多位数乘法心算、几何物体旋转或听到一段熟悉的旋律等［4］.运动想像是通过想像肢体运动以诱导感觉运动皮层的EEG 信号发生事件相关(去)同步(ERD/ERS)现象［5］.目前所能分离的运动想像任务最多为四维，分别是左手、右手、脚和舌头［6］.随着分类维数的增加，达到满意的分类精度需要的训练时间也增加.

针对上述情况，研究者们努力探索其他实验范式.DaSalla 等［7］提出用英语元音/a/和/u/作为默读内容，被试默读这些元音实现BCI 系统的操作.文献［8］分析了被试在默读“左”和“壹”这2 个具有不同含义、发音和字形的汉字时的脑电信号.本文将利用时频分析算法来分析汉字默读时的脑电信号，并判断其是否具有类似于运动想像范式的ERD/ERS 现象.共空间模式(CSP)极易受到噪声的干扰，选择合适的滤波范围尤为关键，时频分析可为滤波范围的确定提供理论依据.

1 实验设计

8 位身体和心理均健康的被试均为右利手，年龄为22 ～26 岁，均未参加过任何相关的脑电实验.首先向被试介绍实验目的和相关注意事项，然后被试签署《知情同意书》.实验所需的提示信息采用液晶显示屏完成，被试与显示屏的距离为1 m.

实验为无反馈实验，整个流程的提示反复随机显示“左”和“壹”这2 个汉字.实验时序图见图1.每次实验开始，屏幕中央显示一个“* ”，表示空闲期，持续时间为2 s.然后，显示“+”，时长为1 s，提醒被试做好准备.准备期后显示1 s 的汉字.随后4 s 的想像期需要被试持续默读提示的汉字.想像期过后，单次实验结束，被试可休息1 s.被试在默读时不能发出声音、触动嘴唇或者伸缩舌头.每个汉字在一组实验中各随机出现15 次.每个被试一共完成5 组实验，并在每组实验之间休息5 min.因此，每个汉字共获得75 组实验数据.

图1 实验时序图

2 数据采集

语言作为大脑的高级认知行为，涉及的大脑皮层较为广泛.普遍观点认为语言主要是被大脑左半球处理，包括2 个区域:位于颞区后上方的威尔尼克区和位于额下回后方的布洛卡区;另有研究显示右脑对语言的处理也有贡献［9］.此外，被试默读汉字时可能有想像舌头和咽喉的运动.有别于文献［8］，本文实验采集的数据来自能够覆盖威尔尼克区、布洛卡区和感觉运动皮层的35 个电极.该电极遵照国际10/20 导联系统设置(见图2).采集脑电信号的设备为美国Neuroscan 公司生产的SynAmps2 系统.为了去除眼电的干扰，被试的左右太阳穴和左眼上下位置分别安放2 个双极性电极以记录水平眼电和垂直眼电.接地电极放置在前额以减少50 Hz 工频干扰，采集脑电信号时的参考电极设置在头顶.所有电极的阻抗低于5 kΩ.信号的采样率为250 Hz，采集后的信号经过0.1 ～100 Hz 的带通滤波器滤波后存入计算机.

图2 电极位置

3 分析方法

3.1 时频分析

EEG 信号是一种非平稳时变信号，仅在时间域或者频率域无法全面分析该信号.时频分析能同时在时间和频率范围内对信号进行分析，即把时间域内的信号在时间和频率的平面中展开，将其表示为以时间和频率2 个参数为自变量的函数，进而展现该信号各频率成分沿着时间轴的变化情况，分析结果更加直观和精确.该分析方法既能突出非平稳时变信号的瞬态特征，又可作为BCI 系统的特征提取方法.常用的时频分析算法有短时傅里叶变换和连续小波变换等［10］.

为了研究默读汉字时脑电信号能量相对于空闲状态时的改变情况，采用事件相关谱扰动(ERSP)作为信号的时频分析算法［11］.ERSP 算法是将多次实验数据的能量谱叠加，并计算该能量谱相对于刺激出现前基线能量谱的变化值.如果刺激出现后的能量谱有所提高或降低，则说明其对于时间有锁时关系;反之则没有.为了保留θ 波、α 波和β 波(三者的频率范围分别为4 ～7，8 ～12，13 ～35 Hz)的信息，EEG 经过4 ～45 Hz 带通滤波器滤波后用EEGLAB工具箱绘制ERSP 图［12］，将ERSP 计算值设为p.被试S2 和S6 分别默读2 个汉字时电极P5 的变化最为明显，该电极的分析结果见图3.

图3中时间轴的0 s 位置对应图1中的第3 s，即提示出现的时刻.由图3(a)和(b)可知，被试S2分别默读2 个汉字时，电极P5 的EEG 信号在8 ～16 Hz 范围内的能量相对于基线有所增加.由图3(c)和(d)可知，被试S6 默读2 个汉字时，电极P5的EEG 信号在9 ～24 Hz 范围内的能量相对于基线有所减少.所有被试的EEG 信号能量变化频率范围均可通过这种方法计算得到.默读汉字可诱导EEG 信号的能量增加或减少，可将这一特征作为BCI 系统的控制输出.

在确定能量变化的频率范围后，将图1中8 s 内某个范围的EEG 信号能量设为A，提示发生前的参考时间段(此处选为空闲期中第2 s)内能量设为R.8 s 内的能量相对于提示发生前的能量波动百分数D 为

如果想像期的D 值大于0，则表示默读汉字可引起EEG 信号的能量增加，为事件相关同步(ERS);反之，为事件相关去同步(ERD).ERD/ERS 统称为事件相关(去)同步，计算流程参考文献［5］.

3.2 特征提取

CSP 作为一种监督的特征提取算法，已成功应用于提取运动想像EEG 信号的特征值［13］.CSP 算法将2 个相应的协方差矩阵同时对角化，以得到一个空间滤波器.多通道的EEG 信号能通过该滤波器投影到低维子空间上，投影后的成分包含各通道的权值及其间的相互信息，因此信号的信噪比增加.详细计算过程参考文献［13］.利用CSP 算法计算得到最佳投影矩阵W，W－1的列向量为共空域模式，被认为是脑电信号在时间恒定时的源分布向量.每次实验中的EEG 信号根据电极数和采样数可组成矩阵E，则分解矩阵Z 为

对于每一次的实验数据，仅计算Z 中2m 个向量的方差，以构造特征向量.同时，对2 类协方差矩阵对角化后得到2 组特征值，构建空间滤波器时可选取这2 组特征值中最大的m 个值所对应的特征向量.利用该滤波器对脑电信号滤波后得到矩阵Zj(j=1，2，…，2m)，根据上述计算，矩阵Z 的前m 行和后m 行与Zj相同.特征向量Vj为

图3 被试S2 和S6 分别默读2 个汉字时的ERSP 图

式中，var(·)为方差计算函数;m=2.经过对数变换，将其转换为正态分布.

3.3 特征分类

支持向量机(SVM)是在VC 维和结构风险最小化原理的基础上建立的，并能根据有限样本信息在模型的复杂性和学习能力之间寻求最佳折中［14］.SVM 具有优良的范化能力，特别适用于解决小样本、非线性和高维度模型等问题.对于n 个非线性可分的训练样本xi(i=1，2，…，n)，分类关键是寻找一个满足Mercer 条件的核函数.样本向量可被核函数映射到一个高维空间，然后在该空间中寻求最优线性分类面.常用的核函数包括线性核函数、多项式核函数、径向基核函数和Sigmoid核函数.本实验选择的是径向基核函数K(xi，xc)，其中，xi为训练样本;xc为核函数的中心，其函数表达式为

式中，g 为控制函数作用范围的宽度参数.计算SVM 采用的软件为LIBSVM 工具箱［15］.为了解决线性不可分的状态而引入惩戒因子h，通过10折交叉验证的网格搜索法寻求g 和h 的最优值，两者的搜寻范围均为2－10～210.

4 结果与分析

每个被试默读2 个汉字各75 次，每次的有效数据时长为8 s.首先，将这些数据用共平均参考滤波(每个电极的原始信号分别减去所有电极原始信号的均值);然后，通过时频分析确定数据所适合的频带范围.根据该频带范围计算每个电极事件相关(去)同步的能量波动百分数，并对单次实验数据提取特征值和分类.

被试S2 和S6 的能量波动百分数见图4.

图4中的时间轴对应于图1，根据时频分析得到2 位被试能量波动百分数所用的滤波范围分别为8 ～16 和9 ～24 Hz.被试S2 默读2 个汉字时，能量相对于空闲状态迅速增加，电极P5 具有ERS 特征，能量波动的最大百分数能达到700%.被试S6默读汉字时电极P5 具有ERD 特征，能量波动的最大百分数绝对值小于80%.ERS 在提示出现1 s后才趋于稳定，而ERD 达到稳定状态的速度更快.因此，BCI 系统能更快地检测到ERD，有利于提高信息传输速率.但是ERS 的波动幅度大于ERD，表明其抗噪声干扰的能力更强.2 个不同汉字被默读时的能量变化趋势是一致的，说明通过EEG 信号的α 波和β 波区分默读2 个汉字时的思维活动比较困难.通过时频分析可确定8 位被试的优化频率Fo范围，结果见表1.

图4 被试S2 和S6 的能量波动百分数

表1 8 位被试的优化频率范围 Hz

8 位被试的优化频率范围各不相同，其中4 位被试的ERD/ERS 主要出现在α 波，其他3 位被试则是在α 波和β 波.被试S3 在默读汉字时脑电信号的能量无明显变化，因此优化频率范围设为4 ～35 Hz.这些被试电极所出现的ERD/ERS 也不尽相同，位于被试S1 大脑中央和顶部的电极具有ERS特征，而被试S2 只有顶部电极具有ERS 特征.被试S4 和S7 额部电极具有ERD 特征而顶部电极具有ERS 特征.被试S5 和S8 的额部、中央和顶部电极均具有ERS 特征，但是被试S6 的这些部位电极则具有ERD 特征.由此可知，ERD/ERS 不仅在被试间存在差异，而且同一个被试的不同电极上也出现不同的ERD/ERS.此外，所确定的ERD 或ERS 并不是绝对的，而是总体趋势.例如，被试S4 默读汉字时，电极P3 在出现短暂的ERD 后呈现较长时间的ERS，因此最终判断该电极具有ERS 特征.分析所有ERD/ERS 能量波动最大百分数的结果发现，ERS 的绝对值普遍大于ERD 的绝对值.

虽然文献［8］表明默读2 个汉字时的脑电信号难以被区分开，但是可分辨出被试在默读汉字时和空闲时这2 个状态间的脑电信号.因此，本文在计算分类正确率时，未给出2 个汉字间的结果，仅比较2 个汉字分别对应各自空闲状态的正确率.为了探讨频率范围优化后分类正确率是否得到改进，选择固定频率F 范围(4 ～35 Hz)与表1中的优化频率Fo范围进行比较，两者所用的滤波器均为6 阶巴特沃兹带通滤波器.此外，被试们默读汉字1 s 后EEG信号能量才逐步稳定，故默读汉字的EEG 信号选择为图1中想像期第5 ～第7 s 的数据.选择时间长度同样为2 s 的空闲期数据与其进行对比，将想像期和空闲期的数据分别组成2 类信号，然后用10 ×10 交叉验证方法分别训练和检验空间滤波器和分类器，结果见表2.由表可知，优化频率范围后2 个汉字的分类准确率分别提高2.25%和1.39%.尤其对于被试S1，其默读“壹”的正确率提高7.27%.因此，利用时频分析所优化的频率范围能改善大部分被试的分类结果.被试S3 的ERSP 图没有发现明显的能量变化，故其F 和Fo设为相同值.被试S3 在默读汉字时的EEG 信号与空闲期相比区分度不大.除被试S3 外，其他被试在默读汉字时，EEG 信号的能量明显变化，该特性可用于控制BCI 系统.

表2 语言想像实验范式的分类正确率 %

5 结语

时频分析能直观显示信号各频率内能量随时间的变化情况.本文选择将ERSP 作为时频分析的算法，该算法能突出思维活动时EEG 信号的能量相对空闲时的波动状态，尤为适合于语言想像和运动想像这类能量锁时的实验范式.在确定合适的滤波范围后，能更合理地判断EEG 信号的ERD/ERS 特征，并能减少噪声对CSP 空间滤波器的影响，因此大部分被试的分类正确率在优化滤波频率范围后得到提高.

采集信号的电极是广泛分布于大脑各块皮层的，故信号的空间特性也特别重要.下一步工作是将时间、频率和空间这3 个特征相结合以分析语言想像的EEG 信号.

References)

［1］ Yuan H，He B.Brain-computer interfaces using sensorimotor rhythms:current state and future perspectives［J］.IEEE Transactions on Biomedical Engineering，2014，61(5):1425-1435.

［2］ McFarland D J，Krusienski D J，Sarnacki W A，et al.Emulation of computer mouse control with a noninvasive brain-computer interface［J］.Journal of Neural Engineering，2008，5(2):101-110.

［3］ Hwang H J，Kim S，Choi S，et al.EEG-based braincomputer interfaces:a thorough literature survey ［J］.International Journal of Human-Computer Interaction，2013，29(12):814-826.

［4］ Friedrich E V C，Scherer R，Neuper C.Stability of event-related (de-)synchronization during brain-computer interface-relevant mental tasks［J］.Clinical Neurophysiology，2013，124(1):61-69.

［5］ Pfurtscheller G，Lopes da Silva F H.Event-related EEG/MEG synchronization and desynchronization:basic principles［J］.Clinical Neurophysiology，1999，110(11):1842-1857.

［6］ Gouy-Pailler C，Congedo M，Brunner C，et al.Nonstationary brain source separation for multiclass motor imagery［J］.IEEE Transactions on Biomedical Engineering，2010，57(2):469-478.

［7］ DaSalla C S，Kambara H，Sato M，et al.Single-trial classification of vowel speech imagery using common spatial patterns［J］.Neural Networks，2009，22(9):1334-1339.

［8］ Wang L，Zhang X，Zhong X，et al.Analysis and classification of speech imagery EEG for BCI［J］.Biomedical Signal Processing and Control，2013，8(6):901-908.

［9］ Harpaz Y，Lavidor M，Goldstein A.Right semantic modulation of early MEG components during ambiguity resolution［J］.NeuroImage，2013，82:107-114.

［10］ Wacker M，Witte H.Time-frequency techniques in biomedical signal analysis［J］.Methods of Information in Medicine，2013，52(4):279-296.

［11］ Makeig S.Auditory event-related dynamics of the EEG spectrum and effects of exposure to tones［J］.Electroencephalography and Clinical Neurophysiology，1993，86(4):283-293.

［12］ Delorme A，Makeig S.EEGLAB:an open source toolbox for analysis of single-trial EEG dynamics including independent component analysis［J］.Journal of Neuroscience Methods，2004，134(1):9-21.

［13］ Ramoser H，Muller-Gerking J，Pfurtscheller G.Optimal spatial filtering of single trial EEG during imagined hand movement［J］.IEEE Transactions on Rehabilitation Engineering，2000，8(4):441-446.

［14］ López Chau A，Li X，Yu W.Convex and concave hulls for classification with support vector machine［J］.Neurocomputing，2013，122:198-209.

［15］ Chang C C，Lin C J.LIBSVM—a library for support vector machines［EB/OL］.(2014-04-01)［2014-04-08］.http://www.csie.ntu.edu.tw/ ～cjlin/libsvm.