闫秘,白雪梅,张晨洁,郭家赫
(长春理工大学 电子信息工程学院,长春 130022)
脑-机接口(brain-computer interface,BCI)是一种将大脑活动转换为计算机或其他设备命令的通信系统,可以代替人的肢体或言语器官实现人与外界的交流以及对外部环境的控制[1]。有很多种生理电信号可以应用在BCI系统中,其中P300信号以其成分特征鲜明,易于识别且不需要受试者进行训练便可得到的优点受到了国内外研究人员的青睐。P300是事件相关电位(Event Related Potentials,ERP)中重要的内源性成分,可以清楚地反映出大脑认知过程中神经电位的生理变化。当人受到小概率事件刺激并进行认知加工后约300~500 ms时会产生一个可以在头皮上记录到的区别于未受刺激脑电信号(electroencephalography,EEG)的波峰。通过识别EEG中是否含有P300成分,可以判断大脑的真实意图,以便于后续工作的进行,如外周神经和肌肉损伤的患者应用脑电信号驱动轮椅运行和丧失语言能力的患者通过脑电信号进行字符拼写等。
近年来,随着科技的发展及社会需求的增加,国内外研究人员提出了很多算法来提高对P300的分类准确率。特征提取是脑电信号预处理、特征提取与分类这三个过程中最重要的环节,很大程度上决定了最终的分类准确率,因此,找到合适的方法对P300特征进行提取是十分重要的。常见的特征提取方法有共空间模型(Common Spatial Pattern,CSP)法[2]、自适应回归模型(Adaptive Autoregressive Model,AAR)[3]、谱分析法、计算波形参数法、小波变换法[4]、时域能量熵和独立成分分析法(Independent Component Analysis,ICA)[5]等。其中CSP方法仅考虑两类任务在空间上的可分性,极大的限制了BCI的实用性;AAR的参数虽然能够清晰地反映大脑的状态,但该方法更倾向于对平稳信号分析;频谱分析法,如快速傅里叶变换(fast Fourier transforms,FFT)和功率谱估计等,这些算法具有简单易实现的优点,但FFT无法有尺度变化且功率谱分析分辨率有限;计算峰值、波形面积等指标作为时域特征能够直观的展现出目标波形的特点,但通常情况下这些指标在不同的波形中差异较大,无法统成为一个标准。目前小波变换、计算能量熵和ICA等特征提取方法受到了广泛关注。小波变换解决了FFT无法进行尺度变换的缺点且具有低频时频率分辨率高的特点;时域能量熵值可以清楚地区分出信号是否存在P300波峰;ICA能够在未知混合信号中分离出原信号且可同时处理多导信号。这些方法均单独利用了P300时、频、空域特征,无法充分体现出P300的本质特征。
针对上述问题,本文结合信息融合中相关概念提出了一种多域特征提取方法,分别提取P300的时域特征、频域特征和空域特征并将这些特征融合成特征集合,全面表征P300脑电信号的特征,再将这个特征集合送入支持向量机(support vector machine,SVM)中进行分类,有效提高分类准确率。
本文采用的数据来自于2004年BCI Competition III data set II中的P300字符拼写实验数据[6]。该实验是在1988年由Farwell和Donchin提出的,数据中包含两个受试者A和B的多次实验数据。
实验过程中,受试者头部戴有64通道的脑电帽,电极按照国际10-20系统标准放置,如图1所示,采样频率设置为240 Hz。受试者面对如图2所示的6×6字符矩阵,该矩阵中的行和列随机地连续闪烁,每行或每列闪烁一次,被称为一个trail,十二个行和列均闪烁一次称为一个block,其中包含2次目标刺激和10次非目标刺激。对于每个字符,这样的过程重复15次,称为一个run,为了方便理解,给出如图3所示时序图。在实验开始前已经为受试者设定了待拼写的目标字符(矩阵左上方),当目标字符所在的行或列闪烁后受试者的EEG中就会产生一个P300波形,受试者默记这些行和列闪烁的次数,而对其他字符的行和列闪烁不作计数。经过数据处理后,SVM分类算法能够根据上述特点判断出目标字符所在的行或列,进而推断出目标字符。
图1 EEG采集电极位置图
图2 诱发界面及其行列编号
图3 P300拼写实验时序图
脑电信号是强随机性、微幅值的非平稳信号,极易受到外界因素干扰,因此采集到的脑电信号中往往惨杂着各种各样的噪声干扰,如工频干扰、心电伪迹及肢体运动产生的伪迹等,其信噪比低[7]。预处理过程是脑电信号处理过程中的第一步,目的是降低夹杂在原始EEG中的噪声和伪迹干扰,使特征提取算法能够提取到反映大脑真实意图的特征。
有关研究表明,并不是所有电极采集到的EEG都存在明显的P300波形,而顶枕处电极组合对P300分类具有有效性[8],所以本文采用Cz,C1,C2,Pz(对应序号:10,11,12,51)进行后续处理,这样可以极大减少待处理的数据量。然后使用经验模态分解(Empirical Mode Decomposition,EMD)算法和FastICA算法对采集到的EEG进行去噪声和伪迹的操作,图4为去除噪声和伪迹前后的对比图,原始EEG中尖峰处为眼电伪迹,经过预处理后较好地实现了降低噪声和去除伪迹的目标。
图4 去噪声和伪迹前后对比图
特征提取为脑电信号处理过程的中间步骤,对分类准确率的高低起着至关重要的作用。在特征提取之前对信号进行叠加平均,使P300信号更加明显。在本实验刺激范式下,闪烁时间持续100 ms,距离下次闪烁的时间间隔为75 ms,即两次闪烁刺激的时间间隔是175 ms。根据P300的时间特性及刺激的时间特性,截取每次闪烁刺激后650 ms的数据段,并将相同刺激对应的数据进行15次叠加平均。
假设一次刺激产生的脑电信号[9]为:
式中,xi(t)为记录到的原始脑电信号;si(t)为P300信号,单次刺激的功率为P;ni(t)为噪声,方差为σ2,均值为0;则单次刺激的信噪比为:P/σ2。经过N次叠加平均后:
式中,X(t)为经过N次叠加平均后的脑电信号;S(t)为经过N次叠加平均的P300信号。
经过N次叠加平均后P300信号的功率依然为P,而噪声的方差(功率)为σ2/N,则叠加平均后的信噪比为N⋅P/σ2,相比于未进行叠加平均的信号的信噪比提高N倍。
以Cz导联为例,对数据分段并进行15次叠加平均后得到图5所示的效果图。从图中可以看出,经过叠加平均的信号可以将P300更加清晰的显示出来。该操作为特征提取提供了有力帮助,使得P300信号的特征提取更为准确,为分类准确率的提高提供了保障。
图5 包含P300的信号和不包含P300的信号分段叠加平均前后对比
P300脑电信号是与认知活动过程密切相关的诱发脑电信号,在时域上有明显的特点,即,小概率目标事件刺激受试者约300 ms后会出现一段在整段信号中波幅最高的正电位,如图5所示。而时域能量熵反映的是信号能量在时域上的分布情况,可以较好地展示出含有P300成分的脑电信号与不含有P300成分的脑电信号在300 ms左右幅值的差异。根据信息论中对香浓熵的定义,时域能量熵可由式(3)计算:
从式(3)可以看出,某段信号的能量在时域上所占信号总能量的比值越大,那么该段信号的能量熵值Ei越小。受试者受到刺激后产生P300信号,此时间段信号的波幅明显提升,在时域上信号的波幅与能量的大小正相关,因此该段能量占整个信号段能量的比重大,时域能量熵小。而没有诱发出P300的脑电信号在整个信号段上波幅没有明显变化,能量分布较为均匀。
以Cz导联为例,将650 ms的数据分为10段,每段65 ms,即1~65 ms的数据为第1段数据,586~650 ms的数据为第10段数据,分别计算分段后含有P300成分和不含有P300成分的脑电信号的时域能量,并根据公式(3)计算每段信号的时域能量熵,以N1-N10命名。含P300成分与不含P300成分的时域能量熵如图6所示。
图6 有无P300成分的时域能量熵
图6中带有三角形标记的曲线的第5段、第6段、第7段的时域能量熵达到了最小值,分别为2.501 3、1.013 4和3.442 2,即该段能量占整个信号段能量的比重大,对应着含有P300波幅的脑电信号。因此,选择N5、N6、N7这三个时域能量熵值作为时域特征向量。
小波变换(Wavelet Transform,WT)是一种时频分析方法,能够对信号的频率逐层分解,直至达到目标信号所在的频率范围。EEG信号经过小波分解得到不同尺度下的高频细节系数与低频近似系数,这些系数反映着EEG在特定频段的特征,由于P300信号频率范围低,所以用低频近似系数作为频域特征,而小波变换具有在低频部分频率分辨率高的特点,因此非常适合处理EEG这类非平稳的随机信号。
信号f(x)的离散小波变换和逆变换有如下定义:
式中,f(x)为待处理的信号;dj,k为小波系数;为母小波;j为分解尺度;k为时间平移量。
考虑到算法运行速度等现实问题,本文选择Mallat算法对信号进行有限层分解,母小波选择与P300波形相近的Daubechies4小波。有关研究表明,P300信号的频率范围为2~12 Hz,甚至更低频段[10]。而实验数据的采样频率是240 Hz,因此需要对原始信号进行4层小波分解,得到0~15 Hz的低频段信号,即包含P300的频段,从这个频段提取的低频近似系数特征对应着P300的频域特征。以Cz导联为例,最终得到16个低频近似系数,命名为X1-X16。
图7 有无P300成分的小波近似系数
如图7所示,带有圆形标记的曲线表示的是包含P300成分的小波近似系数,星形标记的曲线表示的是不包含P300成分的小波近似系数。因为小波近似系数能够清晰地区分出信号是否含有P300成分,所以选择差别明显的X1-X5和X9-X12这9个小波近似系数作为频域特征向量。
EEG能够全方位辐射到头皮上,但其在大脑不同区域的分布情况不同,所以可以对多导EEG空间分布特性进行分析。以空间信息作为特征,能够达到对EEG有效成分进行空间定位的目的。独立成分分析是目前最受关注的EEG空间分析方法,属于盲源分离算法[11](Blind Source Separation,BSS)的一种。它假设源信号之间相互独立,且源信号和信号混合矩阵都是未知的,观测信号是经过线性混合后的源信号,这样可以用一组在最大程度上逼近源信号的相互独立的分量来表示观测到的信号。
假设X为n维观测信号矢量,S为独立的m维未知源信号矢量,A为混合矩阵,W为混合矩阵A的逆,即A=W-1。ICA的目的就是寻找解混矩阵W,使观测信号X通过解混矩阵W的投影变换得到输出独立分量U,使U为未知源信号S的最优估计。为了将ICA具体化,需要用数学表达式表达:ICA的线性模型为X=AS,U=WX=WAS。ICA的原理如图8所示。
图8 ICA原理图
由于事先选定了C1、C2、Cz、Pz四个通道的数据进行处理,因此得出式(6)用于求解混合矩阵W-1:
混合矩阵W-1中的矩阵元素Wij反映了第j个独立分量Uj到观测信号X的投影,而X是在大脑上不同电极测得,即混合矩阵中的元素可以表征EEG的空间分布特性,所以使用混合矩阵中的4×4=16个元素作为空域特征向量,以W1-W16命名。
通过时域、频域及空域的特征提取进一步降低了数据处理量,得到了Cz导联上的3个时域特征向量和9个频域特征向量;Pz导联上的2个时域特征向量和10个频域特征向量;C1导联上的3个时域特征向量和9个频域特征向量;C2导联上的3个时域特征向量和9个频域特征向量以及所选4导上的16个空域特征向量组成的多域特征集合。
特征分类为脑电信号处理过程中的最后一步。在BCI系统中,对EEG进行处理的最终目的是得到能够反应大脑真实意图的控制信号,使大脑准确控制外部设备运行成为现实。评价EEG处理效果的最终指标是特征分类准确率,而特征提取的有效性直接影响着分类结果的准确性,即分类结果是特征提取效果的一个反映。
支持向量机(Support Vector Machine,SVM)是一种有监督学习模型,在解决高维及非线性模式识别中表现出许多特有的优势,能够有效地避免过学习、维数灾难、局部极小等传统分类方法存在的问题,在小样本条件下仍然具有良好的泛化能力,因此被广泛应用到EEG分类当中[12]。SVM的基本原理是:通过核函数对输入数据进行非线性变换,将其映射到一个高维特征空间中使这些数据线性可分,然后在这个高维空间中求取最优线性分类面,并使分类间隔最大,将不同的特征尽可能正确地分开。
核函数的选择是SVM分类的关键。核函数可以将难以划分的低维空间向量集映射到高维空间中并确保不会增加这个映射过程的计算复杂度。由于径向基RBF核函数与其他核函数相比具有参数较少、复杂度低,更容易进行数值运算的优点[13],所以本文使用径向基RBF核函数,其表达式为:
在实验过程中,A和B两个受试者均做了5组实验,每个受试者识别185个字符,其中85个字符带有对应于各次闪烁的身份信息,即可以直接看出该段数据对应的字符,而另外100个字符不带有上述信息,即无法判断出应用算法识别出的字符是否和给予受试者刺激的字符是否一致。因此,本文选择受试者A的85个字符数据作为训练集,受试者B的85个字符数据作为测试集。SVM对不同特征进行分类的过程如图9所示。
图9 分类流程图
3.2.1 多域特征提取的分类结果
表1为在不同叠加平均次数下本文的分类准确率和其它文献分类准确率的对比结果,可以看出对P300进行多域特征提取,能够明显提高SVM的分类准确率。在本文中进行5次叠加平均后的分类准确率高于对照组[14-15]进行10次叠加平均的分类准确率。
表1 不同叠加平均次数下的P300分类准确率/%
3.2.2 单域特征提取的分类结果
表2为仅使用时域能量熵作为特征的P300分类准确率,表3为仅使用小波近似系数作为特征的P300分类准确率,表4为仅使用混合矩阵元素作为特征的的P300分类准确率。
表2 以时域能量熵为特征的P300分类准确率/%
表3 以小波近似系数为特征的P300分类准确率/%
表4 以混合矩阵中元素为特征的P300分类准确率/%
从上述3个表中可以得出以下结论:
(1)P300的分类准确率随着叠加次数的增加而增加,说明叠加平均次数越高越能够将P300的特征显示出来。
(2)对于单一特征提取方法,小波分析提取近似系数作为特征的方法得到的分类准确率高于以时域能量熵为特征的分类准确率,而使用ICA的空域特征提取方法相比于其他两种方法得到的分类准确率最低。
针对P300脑电信号微弱及单一的特征提取方法无法全面表征P300特性的问题,本文提出一种多域特征提取方法。该方法以信号在时域上的能量信息作为时域特征,以小波分解后得到的近似系数作为频域特征,以ICA分解得到的混合矩阵中的元素作为空域特征。从多角度提取特征,使得P300的本质特征得到近乎全面的表达,不但降低了数据处理量,而且为SVM分类的结果提供了保证。
实验结果表明,采用本文提出的多域特征提取方法得到的分类准确率明显高于单一特征提取的分类准确率。
本文在数据处理过程中,首先对有用电极进行选择,剔除无关电极,使待处理的数据量大幅减少,根据相关研究选择P300相对明显的电极C1,C2,Cz,Pz进行处理,使数据处理量为原始数据的1/16,大大节省了数据处理的时间。然后在特征提取之前进行数据截取及相同刺激下的数据叠加平均操作,提高了信噪比,使P300的特征更加明显。