潘隽锴，马玉良，席旭刚，孙明旭，张建海
(1.杭州电子科技大学圣光机联合学院，浙江 杭州 310018；
2.杭州电子科技大学自动化学院，浙江 杭州 310018；
3.济南大学自动化与电气工程学院，山东 济南 250022；
4.杭州电子科技大学计算机学院，浙江 杭州 310018；
5.浙江省脑机协同智能重点实验室，浙江 杭州 310018)
脑科学是研究脑部结构与功能的学科，其核心是以大脑为认知主体的神经科学[1-2]。脑科学的发展是生命科学中的一个重要领域，是当下最为前沿的研究领域之一。脑机接口(brain-computer interface，BCI)技术是指不依赖人体的外围神经和肌肉组织，通过解码大脑的意识活动而实现人脑与外部设备的通信交流[3]。目前脑机接口的发展受限于识别率低、通信速度慢等问题[4]，而基于稳态视觉诱发电位的脑机接口系统(steady-state visual evoked potential，SSVEPBCIs)由于具有较高的识别准确率、信息传输率，以及实验环境配置简单、受试者仅需进行少量的训练等优点，近十年来备受关注[5-6]。
稳态视觉诱发电位是指当人持续注视某一以固定频率闪烁的视觉刺激时，大脑视觉皮层会在刺激频率或谐波频率处产生明显的电位变化。SSVEP的视觉刺激频率范围一般介于4到50 Hz之间，分为低频段(4 Hz～15 Hz)、中频段(15 Hz～30 Hz)和高频段(30 Hz～50 Hz)，SSVEP响应幅值的全局最大值大约在10 Hz出现[7]，高频段的刺激频率能诱发的响应最小，目前大多数系统所采用的视觉刺激主要集中于中低频段。
SSVEP-BCI系统的性能主要取决于以下三个主要因素:视觉刺激呈现、多目标编码方式和目标识别算法[8]。与其他类型的脑机接口类似，信号处理方法、机器学习、深度学习等方法被广泛应用到SSVEP的检测和识别。基于功率谱密度分析(power spectral density analysis，PSDA)的方法如快速傅里叶变换(fast Fourier transform，FFT)被应用于单通道脑电信号(EEG)的频率检测[9]；
基于最小能量组合(minimum energy combination，MEC)的频率识别算法[10]通过寻找空间滤波器将原始EEG进行投影，得到低维组合信号，以降低噪声和其他伪迹信号的影响；
基于多变量同步指数(multivariate synchronization index，MSI)的频率识别算法[11]通过计算两组多通道信号之间的同步指数，筛选最大同步指数对应的频率进行目标识别。Lin等人[12]首先将典型相关分析(canonical correlation analysis，CCA)用于SSVEP-BCI系统的频率识别，由于其具有不需要参数优化，识别率和信息传输率较高的优势，此后基于CCA的改进算法被广泛研究。Poryzala等人[13]提出了CCA的聚类分析(cluster analysis of CCA，CACCA)；
Pan等人[14]提出相位约束的CCA算法(phase-CCA)；
Bin等人[15]提出了基于个体模板的CCA算法(individual template CCA，it-CCA)；
Wang等人[16]提出了扩展典型相关分析(eCCA)，利用个体训练数据增强SSVEP的检测；
Zhang等人[17]提出了多途径CCA算法(Multiway CCA)以寻找合适的参考信号；
Zhang等人[18]又提出了多集CCA(Multiset CCA)以优化来自参考信号的空间特征；
Nakanishi等人[19]提出了任务相关成分分析(task related component analysis，TRCA)，通过在任务期间最大化复现性提取任务相关成分。
当需要考虑受试者校准数据时，大多数相关算法如it-CCA、eCCA、MultisetCCA等选择直接对训练数据平均化以得到包含特定受试者信息的模板信号，这种传统的构造方式可能稍显粗略，从而导致得到的空间滤波器不够十分准确，影响算法的识别结果。因此本文针对电话拨号稳态视觉诱发电位数据集重新选择系数特征组合，在构造个体模板过程中引入反映各训练数据信号质量的权重系数，提出了一种新的方法即coef-eCCA，并通过实验验证了该方法的有效性。
1.1 电话拨号系统数据集
本数据集由加州大学圣迭戈分校(UCSD)提供，所有实验参与者阅读并签署知情同意书[16]。视觉刺激的布局如图1所示，12个目标排列为4×3的矩阵，作为电话拨号系统的虚拟键盘，形状为为6 cm×6 cm的正方形，两两之间水平和垂直间隔分别为5 cm和1.5 cm。实验采用了频率和相位联合编码的方式，相邻的两个目标同时设置不同的频率与相位。刺激频率设置为9.25 Hz到14.75 Hz，等差间隔0.5 Hz；
相位设置为0到1.5π，等差间隔0.5π，各个目标具体对应的频率和相位信息见图1(b)。刺激界面程序在MATLAB下使用PsychToolbox(PTB)工具箱编写。
共有10名视力正常或矫正至正常的健康受试者(9名男性，1名女性)参与了本次实验，实验在昏暗的房间里展开，受试者坐在距离显示器60 cm的舒适椅子上。这项研究进行了模拟在线BCI实验，记录数据进行离线分析，脑电数据由位于覆盖受试者脑部枕叶区域的8个电极采集。对于每位受试者，实验由15个模块组成，每个模块中受试者被要求以随机顺序注视程序指示的一个视觉刺激4 s，完成对应于12个目标的12次试验。在每次试验开始时，一个红色方块会出现在目标刺激位置持续1 s，见图1(a)，受试者需要在1 s内将目光从上一个目标转移到本次试验的目标上。为了减少眼动伪迹信号，受试者被要求在刺激期间避免眨眼。所有数据被降采样到256 Hz以减少数据量。
图1 电话拨号SSVEP-BCI刺激界面
1.2 数据预处理
有研究人员提出在人的视觉系统中存在140 ms的视觉延迟[20]，因此对于Tw s的时间窗，提取的数据段为[0.14 s，0.14＋Tw s]，时间0显示刺激开始。之后使用零相位切比雪夫无限脉冲响应(infinite impulse response，ITR)滤波器对所有数据进行带通滤波，通带频率为6 Hz到80 Hz，刺激频率为10.25 Hz时滤波前后脑电信号时、频域波形如图2所示，滤波后有效地去除了脑电信号中的基线漂移和刺激频率外频段的干扰。由于实验只记录了位于脑皮层视觉区域的8个电极的脑电数据，所以无需进行通道选择操作。
图2 滤波前后脑电信号时、频域波形图
2.1 典型相关分析
典型相关分析算法(CCA)是一种经典的多元统计方法，用于测量两组变量之间的潜在相关性。其目标是寻找一对线性组合，使得经过变换后的两组变量之间相关性最大[21]。考虑两组多维变量X，Y，X∈RC×Ns为一组多通道的脑电信号，Y∈R2Nh×Ns为一组人工正-余弦参考信号，其长度与X相同，设置如下:
式中:fk为刺激频率，fs为采样率，Nh为谐波数量。将脑电信号与正-余弦参考信号进行典型相关分析，通过寻找一对权重向量wx∈RC，wy∈R2Nh，使得映射后的两组典型变量x＝wTx X和y＝wTy Y相关性最大，相关性由二者之间Pearson相关系数确定:
将多通道的时域脑电信号X与包含所有刺激频率的正-余弦参考信号组Y＝[Y1，Y2，…，Yk]做典型相关分析，选择相关系数最大值对应的频率作为SSVEP的目标刺激频率:
2.2 扩展典型相关分析
典型相关分析是一种无监督学习方法，它最初仅被用来检测频率，但随着越来越多的实验范式采用联合频率相位的编码方式，如何在识别过程中有效利用SSVEP的相位信息变得很重要，而CCA无法被用于区分不同的相位[8]。有学者将标准CCA与基于个体模板的CCA相结合，加入训练数据，提出了扩展典型相关分析(eCCA)。除人工构造的正-余弦参考信号外，对典型相关分析算法进行计算结构的调整，通过将受试者的训练数据平均化，构造了个体模板信号:
式中:Xk为训练数据，Ntrain为训练数据个数。由这3种信号可以得到6种空间滤波器的形式——①测试信号X与个体模板之间:wX()，w＾X()；
②测试信号X与正-余弦参考Y之间:wX(XY)，wY(XY)；
③个体模板信号与正-余弦参考信号Y之间:w＾X()，wY()。将这些空间滤波器映射到各种信号上可以得到多种特征。在Masaki等人的研究中，采用了如下5种相关系数的组合方式[16]:
ρ(a，b)即为两组信号a和b的Pearson相关系数，计算方式同式(2)。采用一个集成分类器用于组合这5种特征作为最终特征用于识别:
在计算测试信号与个体模板信号之间的相关系数过程中会产生负数，因此集成分类器添加sign()项用于保留二者间的负相关信息。最后识别刺激目标同式(3)。
2.3 本文提出的方法
2.3.1 特征系数选择
除去扩展典型相关分析中提出的5个相关系数，由6种空间滤波器与测试信号、正-余弦参考信号和模板信号分别映射，总共可以得到10个典型变量，见表1。由这10个典型变量两两之间进行相关系数计算，可以得到45个系数特征。
表1 10个典型变量
为了减少算法在计算过程中消耗的时间，需要降低算法的计算复杂度。各实验数据进行系数特征选择，在保证识别正确率的同时，应尽量减少特征数量。本文选择了如下3个系数特征，为了与标准eCCA对比我们将其命名为pro-eCCA，图3显示了系数特征提取流程。
图3 系数特征提取流程图
2.3.2 重构个体模板
与标准典型相关分析相比，引入个体模板信号后算法的识别性能得到显著提高，这不仅与模板信号中包含重要的相位信息有关，还因为通过将训练数据平均化可以有效地消除噪声影响[22]。本文在此操作的基础上提出了一种更为细致的个体模板构造方式，即在训练数据中引入各次试验权重系数，得到一种新的算法coef-eCCA。
在信号处理领域中，通常采用频谱幅值和信噪比衡量信号质量的好坏，信噪比(signal-noise ratio，SNR)即信号与背景噪声之比。本文分别以训练数据的频谱幅值与信噪比作为评价指标，给予各次试验信号权重系数。定义频率fk处的信噪比为幅频响应曲线中fk处的幅值与附近L个频率的幅值均值之比[23]:
式中:amp(fk)为SSVEP在频率fk处的频谱幅值，L取值为16，相邻频率间隔Δf为0.125Hz。
之后在构造个体模板中分别采用两种计算方式确定权重系数——第一种为计算频谱中所有12个刺激频率对应的评价指标(evaluation index，EI)的平均值，见式(9)；
第二种为针对各刺激频率，计算其基频与各谐波频率所对应的评价指标，见式(10):
式中:Nh为谐波个数，本文设置为5，即针对频率为9.25Hz时计算9.25 Hz、18.5 Hz、27.75 Hz、37 Hz、46.25 Hz对应的评价指标，其他频率以此类推。
最后对各次试验数据的评价指标做归一化处理，得到构造个体模板中各试次所对应的权重系数:
∑EI为各次试验数据评价指标之和，权重系数wtem∈RC×Ntrain，C为通道数，Ntrain为训练数据个数。重构后的个体模板信号为:
交叉验证(Cross Validation)是在建立模型和验证模型参数时常用的方法[24]，将样本数据切分组合成不同的训练集和测试集，用训练集训练模型，用测试集评估模型预测的好坏，交叉验证适用于样本数据不充足的情况下。本文采用留一法交叉验证(Leave-one-out Cross Validation)，即将某一识别目标15次试验中的14项作为训练数据用于构造个体模板，剩下的1项作为测试数据用于验证，循环进行直至每项数据都被当做一次测试数据，以避免实验过程中的随机因素影响识别结果。
一般将识别准确率(accuracy)和信息传输率(information transmission rate，ITR)共同作为评估SSVEP-BCI的性能指标。信息传输率的计算公式为[25]:
式中:N为识别的目标个数，P为识别准确率，T为单次实验目标选择时间，通常为固定值(T＝Tw＋ts，ts＝0.5 s为两次选择之间的目光转移时间)，ITR的单位为bits/min。
3.1 重新选择系数特征
为了显示重新选择系数特征后的扩展典型相关分析(pro-eCCA)的识别性能，将其与典型相关分析(CCA)、标准扩展典型相关分析(eCCA)和任务相关成分分析(TRCA)作对比。Nakanishi等人认为一组多通道脑电信号是由任务相关成分，即诱发信号和任务无关成分，即大脑自发信号，共同线性组合而成[19]。通过TRCA算法对原始EEG进行线性加权，从中得到仅包含任务相关成分的信号，提高信号的信噪比。我们选择了从0.5 s到2.5 s，5个不同的时间窗用于识别，图4显示了这4种方法的识别准确率和信息传输率。
图4 四种方法在不同时间窗下的识别性能对比
传统的标准CCA算法相对来说识别效果最差，在短时间窗下(Tw＝0.5 s，1.0 s)识别准确率仅为19.33%和53.94%；
但随着时间窗长度的增加，算法的识别准确率也会达到80%(Tw＝2.0 s，2.5 s)。针对本数据集重新选择系数特征组合(pro-eCCA)后，识别准确率较标准eCCA有了一定的提高，特别是在短时间窗(Tw＝0.5 s，1.0 s)的情况下，各提高了1.67%。TRCA作为目前SSVEP识别方面最流行的算法之一，识别正确率在短时间窗下有着较大的优势，时间窗为0.5 s时领先pro-eCCA5.78%；
随着时间窗长度的增加其优势不断减小，在2.0 s时TRCA为98.78%，pro-eCCA为98.67%几乎可以与TRCA持平。信息传输率方面eCCA和TRCA在0.5 s时达到最大值，此后随着时间窗长度的增加不断降低；
而CCA在短时间窗下的识别率过低，随着时间窗长度增加，随之提高的识别准确率依然保证了信息传输率的提高，并在Tw＝1.5 s时达到最大值。
为了比较各算法之间的计算复杂度，我们计算了Tw＝2.0 s下每位受试者识别过程中消耗时间的总和，见表2。重新选择系数特征后的eCCA，由于减少了系数特征的个数从而大大降低了计算成本，每位受试者平均消耗时间相比之前减少了41.1%。TRCA算法的计算时间普遍超过eCCA，这可能是由于TRCA在求解空间滤波器时需要计算以使各试次训练数据之间协方差最大化，而本文数据集有15个训练模块，因此消耗了更多的计算时间。pro-eCCA的计算过程相对来说更为简便，因此在保证了识别准确率的同时，消耗的时间更少。
3.2 重构个体模板
为了说明重构个体模板后算法的识别性能，分别用两种计算方式和两个评估指标，在固定时间窗(Tw＝2.0 s)下进行实验，每位受试者的具体结果见表3，表4。
表3 重构个体模板后每位受试者的识别准确率
表4 重构个体模板后每位受试者的计算消耗时间
从四种方式与重新选择系数特征组合后的pro-eCCA的识别准确率对比上看，计算12个刺激频率的频谱幅值用于确定各试次训练数据权重系数的方式取得了最好的结果，共有4位受试者得到提升(S1，S2，S3，S10)，其中S2的提升幅度最大(1.67%)，而S3则提高至100%，但S7的结果有所下降。计算所有刺激频率信噪比的方式使S2获得了提高，S7的结果与用频谱幅值一样也下降了，在其他受试者没有变化的情况下平均识别准确率有略微提高。相比较而言，计算刺激频率的基频及谐波频率的频谱幅值的效果最差，虽然将S3提升至100%，但是S1，S2，S10的结果均有下降，平均准确率也降低至98.33%。最后一种使用刺激基频和谐波频率信噪比的方式提升了S2和S3的准确率，并且没有受试者的结果下降，因此平均识别准确率也有一定的提高，与第一种方式一样均超过了同时间窗下TRCA(98.78%)的结果。为验证实验结果是否具有统计学意义，我们进行了假设检验，比较第一种方式和pro-eCCA的识别结果，除去二者中识别率始终保持100%的受试，剩下结果的p＜0.05，即表明实验结果存在显著性差异。
从计算消耗的时间上看，使用频率幅值作为评估指标与pro-eCCA相比有了些许下降，这是由于计算复杂度的增加不可避免地提高了计算成本，但都依然保持在一个良好的水平，不会影响实际应用。而使用信噪比作为评估指标则大幅度提高了计算成本，严重消耗了计算时间，其原因为计算信噪比是在计算频谱幅值的基础上进行。两种计算方式对比来看，第一种好于第二种，这是由于针对各刺激频率的训练数据，第一种仅需计算所有12个实验刺激频率的评估指标；
而第二种则需根据不同的刺激频率选择相应的谐波频率进行计算，计算结构上更为复杂，所以计算成本更高。
为了可视化重构个体模板带来的效果，图5显示了重构前后两种模板信号的频谱波形图，图5(b)为计算所有刺激频率对应的频谱幅值作为训练数据权重系数的构造方式所得到的模板。对比二者的频谱图，重构后基频和谐波频率外的成分被十分有效地抑制，增强了信号的SSVEP响应，可以更为明显地显现出高频谐波成分。平均化构造模板后，除去6次谐波，基波和其他谐波的信噪比都低于40 dB；
重新构造模板后从基波到谐波的信噪比均高于前者。
图5 重构前后个体模板的频谱波形图
3.3 讨论
为进一步对比算法之间的性能，我们研究了在特定时间窗下(Tw＝2.0 s)，分别改变电极数目和实验试次数量对算法性能的变化情况，实验结果见图6。
图6 不同电极数目和训练试次下算法的识别准确率
在研究电极数目对识别结果影响的实验中，电极取值从3到8，每种条件下电极组合情况依次为:[PO3，POz，PO5](3个电极)，[PO3，POz，PO5，PO4](4个电极)，[PO3，POz，PO5，PO4，PO6](5个电极)，[PO3，POz，PO5，PO4，PO6，O1](6个电极)，[PO3，POz，PO5，PO4，PO6，O1，O2](7个电极)，[PO3，POz，PO5，PO4，PO6，O1，O2，Oz](8个电极)，具体电极分布位置如图7所示。从图6(a)中我们可以看到，提高电极数量，三种算法的性能都得到了提升，并且由于O1，O2，Oz三个电极位置更靠近脑皮层枕叶区域中央，因此增加的这几个电极所采集的数据对于算法性能的提升更显著(电极数由5变为6时，coef-eCCA识别率提升了4.11%)。
图7 电极位置分布
研究试次数对识别影响的实验中，我们分别选取了3、6、9、12和15次(即训练数据试次为2、5、8、11、14)，结果见图6(b)。提高实验试次数量均会提升coef-eCCA和TRCA算法的性能，而在CCA算法中由于各个试次是完全独立的，它们仅用于交叉验证以提供更加可靠的结果，提高试次数没有让CCA算法性能得到提升。相比于eCCA，TRCA算法本身更加依赖试次数量，需要从尽可能多的训练数据中获取有效信息，因此在试次数较少时算法性能明显弱于coef-eCCA。
此外我们还进行了实验，以研究刺激频率范围对算法识别性能的影响。可能由于数据集设置的刺激频率范围较小，从最小频率9.25 Hz到最大频率14.75 Hz中间只间隔了5.5 Hz。将其分为三个频段范围(9.25 Hz～10.75 Hz，11.25 Hz～12.75 Hz，13.25 Hz～14.75 Hz)，从结果上看并无明显差异，因此没有在论文中展示出来，在此作文字说明。
本文在扩展典型相关分析算法传统的个体模板构造方式基础上，提出了一种更为精细的构造方法。首先针对电话拨号系统的SSVEP数据集重新选择特征系数组合，从识别结果上看不同时间窗下的识别准确率均获得提高，并且在计算时间方面较其他对比方法取得了显著的优势。之后提出了一种通过评估指标在构造个体模板过程中引入权重系数的方法coefeCCA，选择两种计算方式和两种评估指标分别进行实验，在没有大幅度影响计算成本的情况下将识别准确率最高提升至99%，验证了改进方法的有效性。未来的工作中针对不同的数据集，还可以尝试不同的计算方式和其他信号评估指标以提高实际识别效果。
猜你喜欢 电极受试者准确率 浅谈新型冠状病毒疫情下药物Ⅰ期临床试验受试者的护理世界最新医学信息文摘(2022年43期)2022-11-19全陶瓷电极可用于锂离子电池汽车工程师(2021年12期)2022-01-18燃料电池电极反应式的书写技巧中学生数理化(高中版.高考理化)(2021年10期)2021-12-06乳腺超声检查诊断乳腺肿瘤的特异度及准确率分析健康之家(2021年19期)2021-05-23多层螺旋CT技术诊断急性阑尾炎的效果及准确率分析健康之家(2021年19期)2021-05-23不同序列磁共振成像诊断脊柱损伤的临床准确率比较探讨医学食疗与健康(2021年27期)2021-05-13颈椎病患者使用X线平片和CT影像诊断的临床准确率比照观察健康体检与管理(2021年10期)2021-01-03巧妙解决燃料电池电极式书写问题高考·中(2019年6期)2019-09-10适度惩罚让人进步现代家长(2016年11期)2016-12-05重要的事为什么要说三遍？党的生活（黑龙江）(2016年8期)2016-08-15