韩 锦,董博文,刘 邈,,许敏鹏,,明 东,
(1.天津大学精密仪器与光电子工程学院,天津 300072;2.天津大学医学工程与转化医学研究院,天津 300072)
脑‑机接口(Brain‑computer interface,BCI)是一门涉及多领域的交叉技术,通过采集脑电、脑磁和脑血氧等信号获取大脑的活动信息,利用信号处理技术解码大脑意图,从而将中枢神经系统活动转化为人工输出,以达到替代、修复、增强、补充或者改善中枢神经系统正常输出的目的[1‑4]。基于BCI的脑‑控技术正是通过测量和分析大脑信号并将其转化为外部设备的控制信号,从而完成相应任务,为严重运动障碍患者与外界交流提供了可能,也为智能化控制技术提供了新的技术路径。
研究人员利用大脑不同维度的信息特征开展脑‑控技术多项工作,并取得了较大进展。2010年,Royer等[5]首次采用头皮脑电(Electroencephalography,EEG)信号和智能控制策略,实现对虚拟直升机的三维控制;2013年,Lafleur等[6]开发了脑控无人机系统,操作员通过运动想象(Motor imagery,MI)控制无人机的飞行模式,并可通过实时的视觉反馈信息来调整其运动方式;2015年,Karavas等[7]开发了一种新型脑‑控系统,使飞行员能够同时操控多架无人机进行编队飞行;2017年,Meriño等[8]设计了一种异步脑控无人机系统,实现了对无人机6个自由度的控制,相对基于MI的无人机系统,该系统具有更多控制指令的潜在优势;2018年,Chen等[9]提出了15指令的脑‑控机械臂系统,验证了便携式脑‑控技术的可行性;2019年,Edelman等[10]利用无创BCI技术,首次实现了对外部机械臂的连续控制,该系统可以帮助患者实现生活自理,证明了脑‑控技术可以改善患者生活的生活质量。
虽然脑‑控技术发展迅速,但多采用对单人的脑信息解码,即单人脑控方式,仍处于较初级的形态。面对复杂条件下的工作任务,如在航天领域,航天员在着服情况下,操控舱外机械臂等场景,单人脑控存在执行效率差、可控自由度低的问题,难以高效完成操控任务[11‑13]。针对此问题,本研究利用时分多址和频分多址的混合编码方法,设计了“双人协同”操控策略,通过任务相关成分分析算法(Task‑related component analysis,TRCA)解码P300与稳态视觉诱发电位(Steady‑state visual evoked potential,SS‑VEP)混合脑电特征,开发了108指令高速率“双人协同式”脑‑控机械臂系统,实现对汉字一笔一划的书写,有效扩展了BCI交互方式,为协同式BCI研究提供了参考。
共有10名健康受试者(4名女性和6名男性,年龄在21~28岁,都是右利手)参加离线实验,其中8人参加了在线实验。天津大学机构审查委员会批准了本次实验。根据赫尔辛基宣言,所有受试者均已充分了解所有实验程序,签署了知情同意书,并声明他们已清楚研究所有可能的后果与收益。
脑‑控机械臂系统由BCI子系统和机械臂控制子系统组成。图1展示了整个系统的控制框图,其中BCI子系统包括信号采集、预处理和视图解码算法3部分。系统流程如下:受试者通过注视视图编码界面,诱发相应的特征脑电,经过预处理和视图解码算法解码受试者信号,然后将识别结果以视‑听觉反馈的形式先呈现于视图编码界面,受试者可在3 s内判断结果是否正确。若错误,可通过检测咬牙信号取消该次指令,并开始下一轮次的编码刺激;若正确,该结果会通过TCP/IP通信协议,经坐标映射算法转换为机械臂的6‑D空间坐标,后将其传入指令队列进行缓存,实现对机械臂的实时控制。
1.2.1 视图编码界面
视图编码策略采用频分多址和时分多址相结合的混合编码方法,可同时诱发P300和稳态视觉诱发电位脑电特征。如图2所示,首先将108个指令平均分为12个3×3模块,每个模块都是一个独立的P300拼写器,其9个字符按照伪随机顺序闪烁,闪烁形式为不同灰度的矩形,矩形的垂直方向视角为1.49°,水平方向视角为1.78°。每个字符闪烁时间为200 ms,刺激间隔(Inter‑stimulus inter‑val,ISI)为100 ms。所有模块同时开始闪烁,所以在1 s之内108个字符可以完成一次完整的遍历闪烁,称为1个轮次(round)。图3展示了模块1内9个字符的刺激流程图,其中红色箭头代表特定时刻刺激界面的灰度分布情况。与传统P300编码模式不同,字符刺激矩形的灰度以正弦曲线的方式变化,且每个模块的频率和初相位都不相同。为了便于后续对脑电特征进行模式识别,刺激频率的设计均大于12 Hz,且频率间隔以0.2 Hz增加,频率范围为12.4~14.6 Hz,初相位的步长设置可通过公开数据库仿真得到[14],本课题采样步长为0.35π。该视觉刺激是在27寸LCD显示器上呈现,分辨率为1 920像素×1 080像素,刷新率为120 Hz。
图2 BCI子系统——视图编码界面108指令分布Fig.2 BCI subsystem—Distribution of 108 targets on the stimulation interface
图3 混合编码刺激流程图Fig.3 Flow chart of the hybrid encoding stimulus para‑digm
为实现对汉字的书写,本文设计了汉字书写策略,将目标汉字会按笔画进行拆分,每个笔画由首尾两个端点及其端点间的连线指定。视图编码界面的刺激块中心称为端点,则可以通过对目标端点的识别顺序和端点之间连接关系的规定,实现对汉字笔画的书写。因此,在本文中,为简化书写复杂度,每两个端点称为一对,端点对内会自动实现连线操作,而端点对间不进行连线。例如汉字“福”的第一笔对应图4中编号为(10,14)的端点对,第二笔对应(9,18)端点对;对于汉字“书”的第一笔对应图5中编号为(14,24)的端点对,第二笔对应(24,57)端点对。依据此原理,可以实现对汉字一笔一划的书写。图中红色编号和红色背景分别代表需选择一次或两次作为目标端点。
图4 汉字“福”书写示例Fig.4 Chinese character writing example—Fu
图5 汉字“书”书写示例Fig.5 Chinese character writing example—Shu
1.2.2 视图解码算法
本研究特征主要包括P300、SSVEP,以及混合特征,相应的解码算法依次主要采用逐步线性判别分析(Step‑wise linear discriminant analysis,SWLDA)、扩展型典型相关分析(Extended canonical corre‑lation analysis,ECCA)和TRCA算法。其中SWLDA是机器学习领域经典算法之一,实现过程也较为简单,因此下面主要对ECCA和TRCA进行介绍。
(1)任务相关成分分析算法
对混合特征和SSVEP特征的解码,采用TRCA,其核心思想是对于同类样本,寻找一个投影矩阵w=[w j1,w j2,w j3,…,w Nch],使得同类别试次间的协方差和最大[15‑17],其中j1和j2代表导联编号,Nch是采用的导联数量,具体计算过程如下。
首先假设记录到的脑电数据为x(t)∈RNc×Np,其中Np是数据点个数,则同类别内所有试次在最优方向投影后协方差和为
式中:Cov(a,b)代表变量a与b的协方差,S代表不同试次投影前的协方差和,Nt代表同类别下所有的试次个数,h1和h2代表试次编号。为对式(1)进行求解,加入以下限定条件
式中Q代表同类别样本不同导联间的协方差和。综合式(1,2),求解投影矩阵可转化为以下优化问题
通过拉格朗日乘数法可求得,最优解是矩阵Q-1S对应的特征向量。最后,对于Nf个刺激频率,即对于Nf个类别样本,采用列向拼接将投影方向进行集成,得到集成投影矩阵W如下
(2)扩展型典型相关分析算法
典型相关分析(Canonical correlation analysis,CCA)是通过寻找两组高维度向量的投影矩阵,使向量在矩阵投影后的相关性最大,具体过程如下[14,16]。
假设脑电信号为X∈RNch×Np,Yf∈R2Nh×Np为正余弦参考模板,有
式中Nh和Np分别代表刺激频率的谐波个数(本文设置为5)以及数据点个数。假设两类信号X和Y f对应的投影矩阵分别为W x和W y,两类CCA过程为
式中:E[·]代表数学期望,ρ代表皮尔逊相关系数。ECCA基于CCA方法,引入了受试者个体信息,即脑电模板开发的解码算法[15]。ECCA的空间滤波器由以下3部分组成:测试样本和个体平均模板经CCA得到的空间滤波器测试数据和正余弦模板经CCA得到的空间滤波器个体平均模板和正余弦模板经CCA获得的空间滤波器其中i代表刺激频率的索引,则通过空间滤波和融合获得的相关系数向量r i为
依据文献[16],对系数进行加权融合,即
式中k代表式(7)中4个系数的索引,最终通过选取最大对应的频率为目标识别频率。
1.2.3 机械臂控制子系统
该子系统采用Universal Robots公司的UR10系列六轴机械臂,最大移动距离为1 300 mm,最大载质量为10 kg,最大移动速度为1 300 mm/s。对于机械臂的移动可通过TCP/IP协议进行控制,位置由6‑D坐标系(X,Y,Z,RX,RY,RZ)指定,其中X代表机械臂末端相对于基点的高度(基点标定为机械臂的工作台高度),Y和Z代表真实三维空间中水平面的笛卡尔坐标,RX、RY、RZ分别代表了X、Y、Z三轴的旋转角度。在该系统中,BCI子系统可以对机械臂进行实时控制。
1.2.4 坐标映射算法
视图编码界面为2‑D平面(长×宽=m×n),而机械臂为6‑D空间,因此在BCI子系统解码脑电信息得到结果后,需要经过坐标映射算法,将2‑D平面坐标转换为机械臂的6‑D空间坐标。在机械臂标定后,X、RX、RY、RZ是固定常数值。因此,假设机械臂的操作平台的长×宽为M×N,则对于视图编码界面(二维平面,长×宽=m×n)任意一点(x,y)到机械臂书写平面(即YZ平面,长×宽=M×N)的映射实际上是一个线性映射,坐标应为(X,Y-y/(n-1)×N,Z+x/(m-1)×M,RX,RY,R Z)[18]。
脑电数据通过Neuroscan Synamps2系统进行采集,脑电采集系统包含64个脑电通道,4个双极通道,2个高电平通道,此通道与脑电通道完全独立隔离,支持外部设备信号的采集。每个通道采用24位的模数转换芯片,可以准确记录脑电的变化;最高支持20 KHz采样率,实验时,导联阻抗低于10 kΩ,采样率设置为1 000 Hz,经带通滤波0.1~200 Hz,50 Hz陷波滤波器处理后保存。基于国际10/20系统,共采集了Fz、Cz、Pz、PO7、PO8、O1、O2和Oz 8个导联,参考电极为左乳突M 1,接地电极为FPz和Fz连线的中点。
在分类识别时,主要包括两步:(1)模块的识别;(2)模块内字符的识别。当两步骤均识别正确时,则可以输出正确的靶目标,否则目标识别错误。本系统靶字符的识别,主要基于SSVEP和P300特征特性,采用基于滤波器组的TRCA方法和线性判别分析(Linear discriminant analysis,LDA)解码算法,除此之外,还会采用文献[15]和文献[17]提出的混合特征处理方法,该方法的理论基础请见文献[19]。因此,对于模块的识别,主要采用SSVEP特征和混合特征两种处理办法;对于模块内字符识别,主要采用P300特征和混合特征两种处理办法;对于靶字符的识别,主要采用SSVEP‑P300和混合特征两种处理办法。脑电特征的处理参数如下。
对于P300特征,预处理阶段将采样率降为25 Hz,上述8导联数据经过切比雪夫I型滤波器组处理,滤波器设置参数为[1 Hz,10 Hz],特征提取时间窗为[50 ms,800 ms],0 ms代表刺激开始时刻。
对于SSVEP特征,预处理阶段将采样率降为250 Hz,上述8导联数据经过切比雪夫I型滤波器组处理,滤波器设置参数为[XHz,92 Hz](X=11,22,34,46,58,70,82),特征提取时间窗为[140 ms,
340 ms]。
对于混合特征,预处理阶段将采样率降为250 Hz,上述8导联数据经过切比雪夫I型滤波器组处理,滤 波 器 设 置 参 数 为[XHz,92 Hz](X=1,11,22,34,46,58,70,82),特 征 提 取 时 间 窗 为[50 ms,450 ms]。
图6 信号处理流程图Fig.6 Flow chart of signal processing
实验时,受试者坐在距离显示器60 cm的椅子上,眼睛注视着预先设定的目标字符,并默数目标字符闪烁的次数,共5次。实验分为离线和在线实验。在离线实验中,指定的目标字符首先由红色三角提示,三角形视角为0.79°,持续0.7 s。然后所有模块内的字符会同时开始连续5轮的闪烁,持续5 s,即受试者注视相同的目标5次。每轮刺激包括1个目标刺激和8个非目标刺激。所有受试者需要按照伪随机的顺序依次注视108个字符,这些字符分成3组完成(每组36个字符),每组之间可以进行休息,总实验时间约为13 min。离线实验共采集了5组被试,每组2人,共10名。在线实验共采集了4组共8名被试。每组被试需要合作完成对汉字“福”的书写。具体如下,依据汉字书写策略,汉字“福”共16笔画,将其等分,每组的2名被试各完成8个笔画,即依据提示顺序,每名被试完成16个像素端点的注视和识别即可,当识别输出正确后,机械臂将移动到相应位点进行书写。
BCI领域常用的评价指标为信息传输速率(Information‑transfer rate,ITR),其表达式为
式中:N为指令个数,P为系统的分类正确率,T为输出单个指令所用的时间。在本系统中N为108,由于包含5个轮次的闪烁,所以T为1.7、2.7、3.7、4.7和5.7 s(转移视线时间:0.7 s;闪烁刺激时间:1、2、3、4和5 s)。需要注意的是,在线实验中,T需要增加额外的4 s,用于机械臂的移动控制或取消指令输出。
本研究通过时间‑频率‑相位联合编码设计,可以诱发P300和SSVEP特征,并通过相应的机器学习算法对诱发的脑电特征进行识别与分类。文献[20‑21]的研究结果表明,瞬态ERP成分在频域主要集中在1~4 Hz内,即为低频成分。由于本研究设计时将频率刺激设置为12 Hz以上,所以可以通过不同的滤波器处理,分别观察瞬态ERP和SSVEP成分。图7展示了在Cz、Pz、Oz三导联处的瞬态ERP波形。为排除SSVEP的干扰,数据经过带宽为1~11 Hz的带通滤波器处理,并降采样为250 Hz。时间轴0 ms代表刺激开始时刻。从图7可以观察到明显的N100成分,除此之外,在Cz、Pz导联有明显的P300成分,在Oz导联有明显的P100成分。图8展示了12.4 Hz和14.6 Hz频率刺激在Oz导联的波形图,为排除瞬态ERP的干扰,数据经过带宽为11~17 Hz的带通滤波器处理,并降采样为250 Hz。时间轴0 ms代表刺激开始时刻。从图8可以看到在刺激开始(0 ms)后,脑电波形的振荡幅度开始增大,图8(b)波形响应频率大于图8(a),说明振荡节律与刺激频率相关。如图7、8灰色波形所示,由于实验设计中的刺激间隔为100 ms,因此背景脑电会包含10 Hz的固定频率振荡。上述结果表明,已经成功诱发瞬态ERP成分和SSVEP成分,可用于下一步的分类识别。
图7 多导联的瞬态ERP平均波形Fig.7 Average transient ERP at Cz,Pz and Oz
图8 Oz导联处不同频率的SSVEP平均波形Fig.8 Average SSVEP of different stimulation frequencies at Oz channel
离线结果采用留一法(Leave‑one‑out,LOO),包括离线正确率和ITR。图9展示了采用单一特征和混合特征,不同闪烁轮次的模块识别正确率、模块内字符正确率和靶字符正确率。其中误差棒为标准误差,“*”代表p<0.05,“**”代表p<0.01,“***”代表p<0.001,统计检验方法为单因素重复测量方差分析。
从图9可以看出,正确率随着闪烁轮次的增加呈上升趋势。基于混合特征的TRCA方法在模块识别、模块内字符识别方面显著高于其他方法。在靶字符识别方面(图9(c)),基于混合特征的TRCA算法在第1轮次到第5轮次,靶字符识别正确率依次为73.80%、88.98%、93.70%、94.91%、96.29%,显著优于基于传统特征的TRCA‑SWLDA算法,证明了本文提出的混合特征解码算法的先进性。在后续的在线实验中,被试依据各自离线实验结果,设定闪烁轮次进行实验。
图9 不同轮次下模块、模块内字符和靶字符的分类识别正确率Fig.9 Classification accuracy of module recognition,character recognition within modules and target character recognition against the number of rounds
表1列举了8名被试的在线实验结果。在线实验中,两名被试为一组,需要合作完成对汉字“福”的书写,其中被试S1与S2是一组,S3与S4是一组,以此类推。依据书写策略,汉字“福”对应刺激界面的32个位点,被试需要按照位点的提示顺序,控制机械臂移动到相应位点,每名被试共完成16个位点。依据各自被试离线结果,确定在线单指令输出时间,一般为4.7或5.7 s(转移视线时间:0.7 s;闪烁时间:1或2 s;机械臂移动时间:3 s)。从结果中可以看出,被试最高在线正确率达到94.12%,且对应的ITR最高,为77.05 b/min。8人的平均在线正确率为87.92%,对应的平均ITR为66.00 b/min。在线结果表明了双人协同式脑‑控机械臂汉字书写系统的可行性和有效性。
表1 在线实验结果Table 1 Results of online experiment
本文采用时分多址和频分多址联合编码方法,设计双人协同操控策略,开发了基于108指令的双人协同式脑‑控机械臂系统,通过任务相关成分分析算法解码诱发的P300和SSVEP混合特征,将被试的意图控制指令映射到6‑D机械臂,实现了对机械臂的在线实时控制,完成对汉字一笔一划的书写。相比于单人脑控方式,双人协同方式具有更高的操控效率,其相关技术也为协同BCI提供了技术支撑。