郭苗苗 齐志光 王 磊 徐桂芝
(省部共建电工装备可靠性与智能化国家重点实验室(河北工业大学),天津 300130)
脑机接口(Brain-Computer Interface,BCI)技术可以为患有运动功能障碍的病人提供一种新的信息交流方式和控制路径,从而能够绕过周围的神经及肌肉组织,直接让大脑能够与外部辅助设备进行交互,进而提高病人的生活品质[1]。目前基于运动想象的脑机接口系统已经具有了令人满意的分类精度和比较广泛的应用,比如控制四旋翼虚拟直升机等[2- 4]。不过,操作层面少是自发脑电脑机接口最大的不足,迄今为止只有左右手、舌头以及脚四个有限的运动想象任务操作层面[5]。除此之外,还可以通过设想回旋的几何物体以及思考多位数乘法等意识任务实现自主操控脑电信号[6]。
不仅运动想象和以上所述意识任务在平时生活中很常见,默读以及口语交际也是很频繁的生活行为。假如可以提取出由默读或者说话意图此类较为高级的思想行为而诱发形成的脑电图(Electroencephalogram,EEG)的有效特性,将其作为BCI系统的操控输入,则可以丰富脑机接口的范式,并有助于失语症病人的康复。大多数失语症病人口齿不够清晰,然而头脑思路却不受影响。尤其是随着我国疾病谱的改变,约有70%~75%的脑血管病患者出现言语障碍[7],失语症患者正在快速递增。语言的交流在人类交际活动中占了很大的比重,假如脑机接口技术可以解码病人的内心默读,则可以在很大程度上提升其生活质量。不过大脑的构造极其繁杂,迄今为止只凭脑电信号尚做不到对内心语言思想的解码[8]。所以有学者最开始通过默读字母的方式,例如Leuthardt等人经过试验记录受试者分别在默读或朗读OO,EE以及AH和EH等字母时的皮层脑电信号来实现对BCI的操控[9]。Dasalla等人则是提取默读英语元音/a/和/u/时的脑电信号来达到对脑机接口的控制[10]。言语属于高端的认识活动,是人类通过后天的学习获得的本领,因此不一样的学习过程以及由于本身成分所造成的差别都可以引起个人之间的不同[11]。汉字属于表意和平面型的文字,其组成成分和构造形式与拼音形语言有着非常显著的差异[12]。本文选择“喝”、“右”、“吃”与“冷”四个常用汉字进行默读,提取其脑电信号作为脑机接口系统的输入信号,将其脑电信号的特征向量通过共空间模式(Common Spatial Pattern,CSP)和Fisher分类器进行提取和配对,并对算法中的电极数目、分析时间及频率选择进行优化,对比分析优化前后的匹配正确率。
9位被试均身体健康且为右利手,年龄为19到25岁,而且都没有参与过相关的脑机接口实验,实验之前被试均签署了《知情同意书》。本实验通过显示屏提示信息的输出,受试者与液晶显示屏之间的间隔为70 cm。
本实验过程中,显示屏反复随机出现“喝”、“右”、“吃”、“冷”四个汉字。其中一组实验时序图如图1所示,首先显示屏的中央出现一个持续时间为1000 ms的“*”来表示空闲时期的开始;紧跟其后3000 ms的空闲持续期要求受试者继续保持空闲状态;之后显示1000 ms的汉字。显示之后3000 ms的想象持续期要求受试者不断的默读显示屏显示的汉字。想象期过后,单次实验结束。受试者在进行默读时要注意避免发声以及嘴唇和舌头的活动。每一组实验中,四个汉字各随机的呈现15次。9个受试者分别需要完成6组实验,每组实验完成之后会有6 min的休息时间。最终可以获得“喝”、“右”、“吃”、“冷”四个汉字以及空闲状态的实验数据,一共5类,每类分别为90组。
图1 实验中单个汉字显示时序图Fig.1 Schematic illustration of a single trial
本实验采用64导电极帽,电极排列方式按照10-20国际标准。选用美国NeuroScan公司的SynAmps2放大器记录脑电数据,带通滤波的频率范围设置为0.1 Hz到100 Hz,采样率选取1 kHz。实验时,在左眼上下两侧分别放置两个分散的电极用于测量垂直眼电,在双眼两侧放置另外两个分散的电极用于测量水平眼电,并选择头顶作为参考电极。所有电极的电阻均被调至5 kΩ 以下。EEG数据由运行于Windows7平台上的Curry8软件进行记录。
眼电会对采集到的脑电信号产生很强的干扰,首先用Curry8软件的预处理功能来去除脑电信号中眼电信号的干扰成分。正常成人清醒时所测到的EEG主要包含了θ波(4~7 Hz)、α波(8~15 Hz)和β波(16~31 Hz),因此选择4~45 Hz的带通滤波器进行滤波。经过滤波后的EEG信号可通过事件相关谱扰动的方法进行时频分析,按照EEG信号能量的变化选择出与汉字默读相关的时间与频率范围。与运动想象不同的是,除运动皮层外大脑皮层的其他区域也会因为默读汉字而被激活,所以时频分析之后,将获得的EEG信号的特征向量通过共空间模式进行提取,再通过空域分析选择特征通道并对电极数目做出优化,最后使用Fisher分类器进行配对,得到匹配正确率。具体的信号处理流程如图2所示。
图2 信号处理流程图Fig.2 Flow chart of signal processing
脑电信号的能量会因为外界刺激的作用而做出相应的改变,经过事件相关谱扰动(Event Related Spectral Perturbation,ERSP)[13]分析后的EEG信号可以确定其平均功率谱在频域内的改变状况,因此可以利用ERSP对脑电信号进行时频分析。
ERSP具体的计算公式为:
(1)
XK(c,f,t,k)表示时频分布值,t是时间,k表示脑电信号被分成较短的并且相互重合的数据段,c代表导联,f代表频率值。上式计算的值为数据段总数为n的ERSP。
本文选择事件相关谱扰动对EEG信号进行时频分析,用于研究汉字默读所引起的EEG信号变化状况。汉字提示后,脑电信号能量与空闲时期相比,若有显著增强或减弱,则证实脑电信号能量的变化在时间域上有锁时关系,进而将这一特性用于共空间模式进行特征提取,通过Fisher分类器进行配对所得的结果用于脑机接口系统的控制条件。使用MATLAB将所有受试者经预处理后的EEG信号绘制事件相关谱扰动图[14],从而在时间域和频率域上获得各个受试者脑电信号能量变化状况。
共空间模式(CSP)常被用来获取两类模式的特征值,属于监督式的空间滤波算法[15]。CSP能够从多个通道的脑电信号中获取到每类空间分散的要素,其中涵盖了每个通道的权值和彼此间的相互信息。该方法获取的空间滤波器是把两个协方差矩阵同时联合对角化,即第一类方差最大化的同时使第二类方差最小化。本文为多分类任务,针对四个汉字及空闲状态共5类任务的匹配,通过“一对多”模式,将一个汉字作为一类,另外四个任务作为另一类,构建空间滤波器。
步骤1 将两类EEG信号进行归一化,求得的协方差矩阵为:
(2)
步骤2 求出混合协方差矩阵:
C=Cu+Cd
(3)
步骤3 主分量分解矩阵C:
C=UAUT
(4)
矩阵U由矩阵C的特征向量组成,矩阵A是矩阵C特征值的对角矩阵。
步骤4 白化变换:
P=A-1/2UT
(5)
Su=PCuPT,Sd=PCdPT
(6)
步骤5 由Cu和Cd经过白化变换后得到具有相同特征向量的Su和Sd,再分别对其进行主分量分解:
Su=By1BT,Sd=By2BT
(7)
矩阵y1和矩阵y2的和为单位矩阵,Su和Sd具有相同特征向量B。由此可以得出结论:如果Sd的特征值越小,Su的特征值越大,反过来依然正确。将脑电信号进行白化变换,再将其投射于矩阵B的前m以及后m列,则能够获取最优的匹配特性。
步骤6 求投影矩阵W,即空间滤波器:
W=(BTP)T
(8)
则可以将每次实验的EEG数据E分解为Z=WE。
步骤7 求脑电信号经过滤波后的特征值,将分解后的EEG数据ZP(P=1,2,…,2m (9) 其中,矩阵ZP的方差为var(ZP)。 Fisher准则函数是为探求一个投射方位,让d维空间的两类样本依据此方位投射于一条直线。投射后,这两类样本各自均值的差值达到最大,而样本自身的离散程度为最小[17]。本文通过ERSP算法优化滤波范围与时间范围后的脑电EEG数据,经过CSP空间滤波提取到的特征矩阵fp,利用Fisher算法进行特征匹配。Fisher算法的运算步骤如下: 步骤1 求均值向量mi: (10) 式中,Ni是类wi的样本个数。 (11) 步骤3 求样本类间离散度矩阵Sb: Sb=(m1-m2)(m1-m2)T (12) (13) 步骤5 将训练集内所有的样本进行投影: y=(w*)TX (14) 步骤6 计算在投影空间上的分割阈值y0: (15) 步骤7 对于给定的X,计算它在w*上的投影点y: y=(w*)TX (16) 步骤8 根据决策规则匹配,有 y>y0⟹X∈W1 y (17) 对预处理之后的信号,采用事件相关谱扰动计算默读汉字时其能量的时频分布。大脑左半球额叶的布洛卡区与颞上回的韦尼克区等与语言加工处理密切相关,研究这些脑区的变化,在一定程度上能够说明语言活动的大脑机制,电极F5与P5分别靠近威尔尼克区与布洛卡区,所以选择电极F5和F6,以及P5和P6分析,并且经过对比分析,这四个导联得出的ERSP图较为相似,且相对于其他导联效果更为明显和清晰。以被试1位于大脑左半球的F5和右半球的F6电极为例,在默读实验中,两个电极的时频结果如图3所示。 图3 被试1空闲状态和默读4个汉字时的ERSP图Fig.3 The ERSP diagram of subject 1 从图3中可以明显看出被试1的脑电信号能量随默读时间增强或减弱,区分较为明显,可以将时频分析结果用于特征选择。ERSP图中的横坐标表示时间,纵坐标为脑电频率。空闲时期ERSP图的横坐标的0 s与图1中的第0 s相对应,即“*”号开始呈现的时刻。默读时期ERSP图横坐标的0 s与图1中的第4 s相对应,也就是汉字开始呈现的时刻。从图3中(a)可以看出,被试1在空闲状态时,电极F5、F6的脑电信号能量在9~14 Hz的频率区域内能量较基线有明显的减弱。由图3中(b)可知,被试1默读汉字“喝”时,F5与F6电极的脑电信号能量在9~16 Hz频率区域,500 ms到2100 ms时间范围内较基线有明显的增强。由图3中(c)可知,被试1默读汉字“右”时,EEG信号在0到1500 ms时间范围内能量并未发生明显波动,而在1500 ms之后,8~14 Hz频率区域能量较基线有明显增强。由图3中(d)可知,被试1默读汉字“吃”时,EEG信号在300 ms到1800 ms时间内,8~14 Hz频率区域内能量较基线均有明显减弱,而在其他时间无显著变化。由图3中(e)可知,被试1默读汉字“冷”时,EEG信号在300 ms到2000 ms时间内,6~13 Hz和16~22 Hz频率区域内能量较基线均有明显减弱,而在其他时间无特别明显变化。 被试默读汉字所引起的脑电信号能量变动的差异主要体现在α波和β波,而且随时间动态变化,反映了被试默读不同汉字时的思维活动。当ERSP的绝对值超过3 dB时,能量相对基线的变化显著,通过时频能量特征的显著变化更为精确的确定9位被试的时频特征,并与固定的4~45 Hz滤波范围和0~3500 ms时间范围内的脑电特征做对比,改进滤波范围与改进时间范围如表1所示。9位被试的改进滤波范围和改进时间范围各不相同,其中3位受试者能量变化主要选取在α波频段,剩余6位受试者选取范围均包括α波频段以及β波频段。能量变化趋势在个体之间存在差异,而且同一被试不同时间也会呈现出不一样的能量变化趋势。 表1 9位被试的改进频率与改进时间范围 从图3能够看出,大脑皮层的某些区域由于默读汉字而被激活致使脑电信号能量发生变化。为了分析四个汉字默读时EEG信号的空域特征差异,选取改进时间和改进频率后的EEG数据进行对比,以受试者1为例,图4为其分别默读四个汉字的CSP脑电地形图。该脑电地形图来自CSP计算空闲状态与默读状态的W-1列向量的第一列和最后一列的差值[12]。按照共空间模式的理论,此两组数据最具差异性。从图4中可以看出分别默读不同汉字时,其地形图中信号变化有明显的不同。异于左右手运动想象的是汉字默读所诱发的脑电信号在各个脑区均有不同程度的变化。此外与改进频率和改进时间范围相同的是,被试CSP脑电地形图也具有明显的个体差异性。因此不同的人汉字默读时很难选取出同样的最佳电极,不利于基于语言BCI的推行,故可以根据信号间的差异选择出作用效果明显的电极,当CSP的系数大于0.1时,信号间的差异明显,可用于特征通道的选取。被试1选择电极AF3、AF4、F6、C5、C6、P3、PO3、P6为特征通道。电极AF3、AF4、F6位于大脑的额叶区反映了默读时大脑的精神活动以及语言想象激活的布洛卡区。电极P3、PO3、PO6位于枕叶区域则反应了默读时激活的视觉功能。其余电极则是由视觉刺激所诱发顶叶与枕叶区的活动[18]。 图4 被试1的CSP脑电地形Fig.4 EEG topography of CSP for subject 1 为了尽可能在使用电极数目最少的情况下,寻求最高的匹配正确率,找到最佳的导联组合。将64个电极以匹配准确率进行排列,按照CSP脑电地形图选择出的一组特征电极,依次在该组特征电极的每个电极周围同时增加一个导联,直至增加至全导联,查看每组匹配准确率。以被试1为例,按照CSP脑电地形图选择出的特征电极1- 8导作为第一组,每次在8个特征电极的每一个电极附近增加一个导联,即每次递加8导,第二组导联数则为1-16导,第三组导联数为1-24导,第四组导联数为1-32导,第五组导联数为1- 40导,第六组导联数为1- 48导,第七组导联数为1-56导,第八组则为1- 64全导联,各组相应的匹配。被试1各导联组对应的匹配准确率如图5所示。 从图5中可以看出,随着导联组内电极数量的增多,匹配准确率呈降低趋势。导联数的增多有可能会增加CSP空间滤波器的噪声误差,所以EEG信号经过CSP脑电地形图的空域分析选择出最佳导联可以降低噪声误差,进而提高匹配准确率。将9位被试的脑电信号通过空域分析确定其最佳导联个数,所得结果如表2所示。 图5 被试1各导联组匹配准确率Fig.5 The matching accuracy of each channels group for subject 1 被试编号最佳导联个数18个27个38个49个511个68个76个810个97个 本文选择比较刺激前后以及默读不同汉字时EEG信号改变状况,针对空闲状态以及默读4个汉字共5类EEG信号进行区分,匹配结果如图6所示。为提高计算匹配准确率,选择K-折交叉验证(K-fold Cross Validation),这里K选择为10,将90组数据集随机分为10组数量相等的子数据集,其中9组子数据集用于构建共空间滤波器,剩余1组则用于验证匹配。然后将各子数据集轮流替换获取10个不一样的匹配结果,再取平均以及得出标准差。 从图6(a)能够看出,汉字默读会诱发脑电信号产生显著变化,空闲状态匹配正确率基本都超过80%。其中在频率和时间范围未进行改进时,即选择空闲期的前3.5 s与想象期的前3.5 s进行对比的平均准确率为81.3%,只改进时间范围不改进滤波范围的平均准确率为82.7%,只改进滤波范围不改进时间范围的平均准确率为83.6%,同时改进时间范围和滤波范围后的结果为85.1%,尤其是被试5匹配准确率都高于90%。分别改进时间范围与改进滤波范围后所有被试空闲状态与默读状态的平均匹配正确率都有所升高,并且只改进滤波范围相对只改进时间范围正确率要高。当滤波范围与时间范围共同改进之后,匹配准确率更高,相对于未改进时提高了3.8%。特别是被试1、被试9默读与空闲时匹配正确率均提高超过5%。所以通过时频分析而改进的时间范围和滤波范围可以改进多数被试匹配效果。 从图6(b)、(c)、(d)、(e)四组结果中可以看出汉字的匹配准确率并没有空闲状态区分度高。在4~45 Hz固定滤波范围内,改进时间范围较未改进时间范围时,四个汉字的匹配准确率总平均提升为1.46%。当只改进滤波范围时,时间范围不做改进,匹配准确率总平均提升为2.12%,由此可知改进滤波范围同样提升了匹配正确率,并且比改进时间范围的匹配正确率更高。同时改进时间范围和滤波范围时汉字间的匹配准确率总平均提升为3.37%,这解释了经过改进时间与滤波范围后能有效利用脑电信号的α波与β波频段来辨别默读两个不同的汉字,从而可以利用默读某个汉字的范式当作脑机接口操控条件。 图6 匹配准确率Fig.6 Matching accuracy 本文设计了汉字默读实验,选择事件相关谱扰动算法对EEG信号进行时频分析,精准提取时频特征,利用共空间模式进行空域分析,并通过电极优化算法选择出最佳导联组合,最终采用Fisher分类器进行默读特征的分类。ERSP可以凸显认知活动时脑电信号频域与时域的相对变化状况,适合用于研究汉字默读所引起的EEG信号变化状况。通过改进时间范围与滤波范围,并且根据最具差异性特征的CSP脑电地形图选择合适的特征通道,降低CSP空间滤波器的噪声误差。结果表明:四个汉字的平均分类正确率均达到70%。改进时间范围与改进频率范围均能提高汉字的平均匹配准确率,并且只改进频率范围比只改进时间范围分类效果更好。同时,两者均进行改进则匹配准确率则会更高,准确率平均提升达到3.37%。 本项研究有助于语言脑机接口系统的发展,对于言语障碍的医疗康复具有一定的推动作用。目前西方发达国家又将康复过程分为3个阶段:急性期、恢复期和维持期。基于默读发音想象的BCI系统尤其适合于康复阶段的恢复期和维持期,不仅可以改善恢复效果和安全度、提高病人生活质量、减少后遗症和医疗事故,还可以显著减短住院时间进而削减医疗费用。下一步将完成语言在线脑机接口康复系统的设计以辅助语言障碍患者的康复。3.3 Fisher分类
4 实验结果及其分析
4.1 时频域分析结果
4.2 空域分析结果
4.3 Fisher分类结果
5 结论