脑-机接口(brain-computer interface, BCI)系统可以在不依赖外周神经肌肉系统的前提下，建立起人脑和外界环境之间直接交流通路，即通过分析大脑活动中的意图，然后将其翻译成控制外部设备的指令。BCI在许多领域有着广阔的应用前景，比如为那些有运动疾病的人提供与外界环境交流的能力(脑控轮椅、机械臂、字符拼写器等^[1-4])。用于非侵入式BCI系统的控制信号包括P300、感觉运动节律(sensorimotor rhythm, SMR)、稳态视觉诱发电位(steady-state visual evoked potential，SSVEP)以及它们之间的组合^[5-6]，其中，基于SSVEP的BCI具有高信噪比和高信息传输率(information transfer rate, ITR)等优点。在典型的基于SSVEP的BCI系统中，多个刺激以不同的频率闪烁(如图形或者闪光刺激，一般大于6 Hz)，当用户注视其中一个靶刺激时，大脑视觉皮层会在靶刺激频率的基频和倍频处产生SSVEP特征成分，通过检测SSVEP中的显性频率就可以识别出靶刺激。

评价BCI系统的性能有很多指标，比如ITR、准确率、输出指令的速度以及指令集的大小等，其中ITR是目前最常用的评价指标^[7]。用于提高SSVEP-BCI信息传输率的方法有很多，比如优化刺激的编码方式和解码SSVEP的算法等。通过优化刺激范式可以增加指令集的数目，并保证SSVEP信号识别的准确率。Zhang等提出多重频率序列编码方法，实现了用少量频率设计出较多控制指令的目的^[8]。Wang等提出了以更灵活的频率来产生稳健闪烁刺激的频率近似法^[9]。Chen等提出联合频率-相位调制(joint frequency-phase modulation, JFPM)方法，设计出有40个字符的拼写器^[10]；并结合视觉刺激的颜色和亮度特征开发出相互调制的频率编码方式，实现了更好的分类表现^[11]。另外，Yin等将行/列P300和两步SSVEP范式整合在一起，提出了混合刺激呈现范式，得到了64个可选指令^[12]。用于目标识别的SSVEP检测算法的进步同样可以提高ITR，比如早期的以快速Fourier变换为代表的能量谱密度分析方法，实现了27.15 b/min的平均ITR^[13]；后来的研究证明空间滤波技术更加有效，且可以提高SSVEP信号的信噪比，其中应用最为广泛的是典型相关分析(canonical correlation analysis，CCA)方法，实现了58 b/min的平均ITR^[14]。许多基于CCA的算法大幅提升了SSVEP-BCI的性能。比如，Wang等^[15]和Chen等^[10]分别将个体模板和滤波器分析整合到CCA方法中，将平均ITR分别提升到104.6 b/min和270 b/min。最近出现的任务相关成分分析(task-relevant component analysis, TRCA)算法更是将平均ITR提升到325.33 b/min ^[16]。

在保证准确率和指令集的前提下，缩短单个指令的输出时间(即缩短数据长度)可以得到更高的ITR。传统的优化时间方法是静态自适应(fixed stopping, FS)策略，即根据离线分析结果，将最高ITR所对应的单试次呈现刺激的时间选定为最优目标选择时间。但是，由于脑电图(electroencephalography, EEG)数据信号的非线性和不稳定性，每个试次的最优时间是不同的，因此更优的解决方案是动态地为每个试次确定时间。为此，可以将以前只有单次检测SSVEP并输出分类结果的过程，分解为多时间点的顺序检测，在每个时间点都作如下判断：若分类结果可信，则输出指令，顺序检测停止；若不可信，则不输出指令，顺序检测继续，并等待下一个时间点继续判断。相比于单次检测并输出的过程，这种顺序检测的过程可以充分利用传统SSVEP检测算法的累积效应，并有效减小噪声的干扰。Schreuder等提出的基于P300的BCI系统的动态自适应(dynamic stopping, DS)策略，就是一个顺序检测过程^[17]。分类结果的置信度通过基于Bayes的概率模型进行评估，这种自适应DS策略的有效性已经在很多研究中得到了验证^[17-22]。然而，目前还没有基于Bayes的DS策略对于提高SSVEP-BCI性能的可行性和有效性的相关研究。本文通过公开数据库对所提出的方法和用固定数据长度的FS策略进行了模拟在线测试，并将分类准确率、所用的数据长度以及ITR作为评价指标，以验证所提出方法的有效性。

1 方法与实验 1.1 受试者

本文采用Wang等公开的数据库来测试所提出的方法^[23]。被试共有35名(其中：17名女性，18名男性；年龄在17~34岁，平均年龄为22岁)，健康且视力正常或者校正之后正常。在这些被试中，有8名被试之前参加过SSVEP-BCI拼写器的相关实验，其余27名被试是初次参加此类实验。全部被试在实验之前均被要求阅读和签署清华大学研究伦理委员会认可的知情同意书。

1.2 实验设计和流程

用户界面是一个包含40个字符(26个英文字母、10个数字和4个字符)的5×8刺激矩阵。每个刺激是一个140×140像素的方块，呈现在液晶显示屏上，显示屏的屏幕尺寸为23.6 cm，屏幕长宽比为16: 9，分辨率为1 920×1 080像素，刷新频率为60 Hz。字符呈现在刺激中间一个32×32像素的方块中。每两个相邻刺激的水平和垂直距离均为50像素。此外，刺激频率与相位均为线性递增；刺激按照从上到下、从左到右的方式依次排列，对应的频率和相位依次按0.2 Hz和0.5 π的步长增加，第一个刺激的频率和相位为8 Hz和0 π，最后一个刺激的频率和相位为15.8 Hz和1.5 π。在本研究中，采用正弦采样的刺激方法，并用JFPM方法编码视觉闪烁^[10]。刺激程序是使用Psychophysics Toolbox Version 3在MATLAB(MathWorks, Inc.)环境下开发的。

实验包含6轮，每一轮包含40个试次，分别对应于40个字符。每个试次开始于持续时间为0.5 s的用于提示目标字符的红色方块，被试被要求在这0.5 s的提示时间内尽可能快地将视线转移到屏幕上的目标字符，提示时间结束的同时，所有字符对应的方块开始以相应的频率闪烁5 s，然后屏幕出现0.5 s的空白。对于每一轮，被试被要求在刺激闪烁的时间内尽量避免眨眼。为了减缓视觉疲劳，相邻轮次之间会有几分钟的休息。实验被安排在一个微暗、有照明的隔音实验室中进行。

1.3 数据收集和信号处理

EEG数据采集的为Pz、PO5、PO3、POz、PO4、PO6、O1、Oz和O2导联的数据，采样率是1 000 Hz，脑电放大器是Synamp2系统(Neuroscan, Inc.)。参考电极放在头顶，接地电极放在前额叶。所有的数据首先被降采样到250 Hz，然后根据时间触发器进行提取。

一个50 Hz的陷波滤波器被用来去除电源工频噪声干扰。除此之外，还用到了滤波器组分析，所用到的滤波器是Chebyshev Ⅰ型无限脉冲响应滤波器，对于第n个频带，其上截止频率和下截止频率分别是90 Hz和n×8 Hz (n=1，2，…，5)，并且使用MATLAB中的filtfilt()函数进行正向和反向滤波。

1.4 模拟在线实验

模拟在线实验模拟的是真实在线实验的过程。在真实在线实验中，数据是被实时收集的，并随着时间的推移对一个个数据片段进行处理。相比于真实在线实验，模拟在线实验的优点在于可以在同一测试数据上进行不同方法的对比。考虑到视觉通路上的延迟，本文将刺激呈现之后0.14~5.14 s之间的数据段定义为一个试次，并用于后面的分析。所有试次的前5轮被分别用来建立基于TRCA的分类模型，剩下的1轮被用来进行模拟在线测试。在测试过程中，每个试次被一个接一个片段地送到建立好的分类模型中。在动态自适应策略中，每处理完一个数据片段，预测的分类结果和相应的置信度就会被更新。数据片段的原始长度是0.5 s，接下来每次增加的数据长度Δt是0.1 s。一旦预测结果的置信度达到了输出阈值或者输入的数据长度超过了预设的最大值t_max(比如1 s)，那么这个数据片段输入过程将会停止，并且输出当前预测结果。

1.5 分类模型

最近，整合TRCA方法在识别SSVEP上展现了杰出的性能^[16]。在TRCA中，空间滤波器u可以通过最大化试次之间的协方差得到，

$ \begin{align} & \sum\limits_{\begin{smallmatrix} {{h}_{1}}, {{h}_{2}}=1 \\ \ \ {{h}_{1}}\ne {{h}_{2}} \end{smallmatrix}}^{{{N}_{\text{t}}}}{{{C}_{{{h}_{1}}, {{h}_{2}}}}}=\sum\limits_{\begin{smallmatrix} {{h}_{1}}, {{h}_{2}}=1 \\ \ \ {{h}_{1}}\ne {{h}_{2}} \end{smallmatrix}}^{{{N}_{\text{t}}}}{\sum\limits_{{{j}_{1}}, {{j}_{2}}=1}^{{{N}_{\text{c}}}}{{{u}_{{{j}_{1}}}}{{u}_{{{j}_{2}}}}}}\text{Cov}\left( \mathit{\boldsymbol{x}}_{{{j}_{1}}}^{\left( {{h}_{1}} \right)}, \mathit{\boldsymbol{x}}_{{{j}_{2}}}^{\left( {{h}_{2}} \right)} \right)= \\ & \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ {{\mathit{\boldsymbol{u}}}^{\text{T}}}\mathit{\boldsymbol{Su}}. \\ \end{align} $

(1)

其中：x∈R^N_s, N_s是每个试次中采样点的数目，j表示通道的索引，N_c是通道的总数目，h表示训练试次的索引，N_t是训练试次的数目，

$ \mathit{\boldsymbol{S}}={{\left( {{S}_{{{j}_{1}}, {{j}_{2}}}} \right)}_{1\le {{j}_{1}}, {{j}_{2}}\le {{N}_{\text{c}}}}}, {{S}_{{{j}_{1}}{{j}_{2}}}}=\sum\limits_{\begin{smallmatrix} {{h}_{1}}, {{h}_{2}}=1 \\ \ \ {{h}_{1}}\ne {{h}_{2}} \end{smallmatrix}}^{{{N}_{\text{t}}}}{\text{Cov}\left( \mathit{\boldsymbol{x}}_{{{j}_{1}}}^{\left( {{h}_{1}} \right)}, \mathit{\boldsymbol{x}}_{{{j}_{2}}}^{\left( {{h}_{2}} \right)} \right)}. $

为了得到一个有限的解，将$ \sum\limits_{j=1}^{{{N}_{\text{c}}}}{{{u}_{j}}{{\mathit{\boldsymbol{x}}}_{j}}} $之间的协方差限制为1,

$ \sum\limits_{{{j}_{1}}, {{j}_{2}}=1}^{{{N}_{\text{c}}}}{{{u}_{{{j}_{1}}}}{{u}_{{{j}_{2}}}}}\text{Cov}\left( \mathit{\boldsymbol{x}}_{{{j}_{1}}}^{{}}, \mathit{\boldsymbol{x}}_{{{j}_{2}}}^{{}} \right)={{\mathit{\boldsymbol{u}}}^{\text{T}}}\mathit{\boldsymbol{Qu}}=1. $

(2)

其中：Q=(Q_j1j2)_{1≤j1, j2≤Nc}, Q_j1j2=Cov(x_j1, x_j2)。最终的目标函数是

$ \mathit{\boldsymbol{\hat u}} = \arg \mathop {\max }\limits_\mathit{\boldsymbol{u}} {\kern 1pt} \frac{{{\mathit{\boldsymbol{u}}^{\rm{T}}}\mathit{\boldsymbol{Su}}}}{{{\mathit{\boldsymbol{u}}^{\rm{T}}}\mathit{\boldsymbol{Qu}}}}. $

(3)

最优解$ \mathit{\boldsymbol{\hat{u}}} $是矩阵Q^-1S的特征向量。整合的空间滤波器U定义为

$ \mathit{\boldsymbol{U}}=\left[ {{{\mathit{\boldsymbol{\hat{u}}}}}_{1}}, {{{\mathit{\boldsymbol{\hat{u}}}}}_{2}}, \cdots , {{{\mathit{\boldsymbol{\hat{u}}}}}_{{{N}_{\text{f}}}}} \right]. $

(4)

其中：${{\mathit{\boldsymbol{\hat{u}}}}}$_k(k=1, 2, …, N_f)代表第k个刺激的空间滤波器，N_f表示刺激的数目。被用来识别目标字符的特征系数为

$ \tilde{r}=\rho \left( {{\mathit{\boldsymbol{X}}}^{\text{T}}}\mathit{\boldsymbol{U}}, {{{\mathit{\boldsymbol{\bar{X}}}}}^{\text{T}}}\mathit{\boldsymbol{U}} \right). $

(5)

其中：ρ()表示Pearson相关系数，X表示测试数据的试次，X表示训练数据试次的均值。

通过使用滤波器组分析的方法^[10]，可以更好地利用嵌入到SSVEP谐波成分中的信息。原始的SSVEP信号X会首先被分解为不同的子频带成分(X⁽ⁿ⁾, n=1, 2, …, N)，然后根据式(6)对所有子带成分的特征系数(如${\tilde{r}}$_k⁽¹⁾, ${\tilde{r}}$_k⁽²⁾, …, ${\tilde{r}}$_k^(N))的平方进行加权求和：

$ {{{\tilde{r}}}_{k}}=\sum\limits_{n=1}^{N}{w\left( n \right)}\cdot {{\left( \tilde{r}_{k}^{\left( n \right)} \right)}^{2}}, k=1, 2, \cdots , {{N}_{\text{f}}}. $

(6)

其中：n表示子带的索引，N表示子带的总数目，k表示刺激的索引，加权系数w(n)被定义为w(n)=n^-1.25+0.25。在所有的指令类别中，预测的目标类别为最大的累加特征系数所对应的类别。

1.6 动态自适应策略

静止自适应策略对所有的试次使用固定的数据长度，而动态自适应策略动态地为每一个试次优化数据长度。由于不同试次之间的SSVEP特征存在着差异，因此每个试次都存在着各自的最优数据长度，动态自适应策略旨在为每个试次找到最优的数据长度，从而得到更高的信息传输率。本研究探讨了可能应用于SSVEP-BCI的基于Bayes估计的动态自适应策略。

就像1.5节所介绍的那样，预测的分类结果通过累加特征系数决定，因此假设特征系数值和分类准确率之间存在着相关性。基于Bayes的动态自适应策略试图通过计算给定特征系数值下正确预测的概率，来确定分类结果的置信度。特别地，就像图 1所显示的那样，利用所有的训练试次，分别为正确预测和错误预测构建关于特征系数的似然概率密度函数(probability density function, PDF)。这里的特征系数是经过标准化的最大累加相关系数，

$ {{d}_{\text{m}}}=\frac{{{{\tilde{r}}}_{\text{m}}}}{\sum\limits_{i=1}^{{{N}_{\text{f}}}}{{{{\tilde{r}}}_{i}}}}. $

(7)

其中${\tilde{r}}$_m是所有类别中最大的累加特征系数。假设正确预测写作H₁，错误预测写作H₀。似然概率密度函数p(d_m|H₁, t)和p(d_m|H₀, t)用Gauss核密度估计来生成，其中t表示用于分类的数据长度。对不同的数据长度，正确预测和错误预测所对应的先验概率密度分布p(H₁|t)和p(H₀|t)用训练数据在不同数据长度下的分类准确率来估计。

在模拟在线实验过程中，对于每一个新的数据片段，先计算出预测的目标字符所对应的标准化累加特征系数，然后通过Bayes推论更新对应的被正确预测的后验概率，

$ \begin{align} &\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ P\left( {{H}_{1}}|{{d}_{\text{m}}}, t \right)= \\ &\frac{P\left( {{d}_{\text{m}}}|{{H}_{1}}, t \right)P\left( {{H}_{1}}|t \right)}{P\left( {{d}_{\text{m}}}|{{H}_{1}}, t \right)P\left( {{H}_{1}}|t \right)+P\left( {{d}_{\text{m}}}|{{H}_{0}}, t \right)P\left( {{H}_{0}}|t \right)}. \\ \end{align} $

(8)

一旦后验概率达到了阈值P_thre，预测的结果就会被输出，而阈值是通过网格法搜索离线确定的。如果当前的试次一直达不到阈值，那么预测结果也会在数据片段的长度达到预设最大值t_max之后被输出。

1.7 性能评估

除了分类准确率之A_c外，ITR也是目前最常用的评估BCI性能的参数之一，它同时考虑了准确率、刺激的数目和时间3个指标^[7]：

$ \begin{align} &\text{ITR}=\left( {{\log }_{2}}{{N}_{\text{f}}}+{{A}_{\text{c}}} \right.{{\log }_{2}}{{A}_{\text{c}}}+\left( 1-{{A}_{\text{c}}} \right)\cdot \\ &\left. \ \ \ \ \ \ \ \ \ \ {{\log }_{2}}\left( \frac{1-{{A}_{\text{c}}}}{{{N}_{\text{f}}}-1} \right) \right)\times 60/T. \\ \end{align} $

(9)

其中T是完成一次输出所需要的时间，包括目标注视时间(即所用到的数据长度)，以及本文中额外的用于实际应用的视线转移时间0.55 s^[23]。除了实际的数据长度，还存在一个理论的数据长度，它可以在最短时间内提供正确的分类结果。如果动态自适应策略中基于Bayes的概率估计方法对分类结果的正误判断不出错的话，那么错误就只会发生在预设的最大数据长度所对应的时刻。对所有的模拟在线测试的试次，可以计算出最优的平均数据长度、准确率，进而计算出理论信息传输率。

2 结果 2.1 模拟在线测试结果

图 2为动态自适应策略与整合TRCA分类模型的分类准确率和ITR的对比。本文并没有使用JFPM方法将刺激的相位间隔从0.5π改成0.35π。从图 2a中可以看出，由于DS策略在相同的数据长度上得到了比整合TRCA更高的准确率，因此图 2b中DS策略得到了更高的信息传输率。整合TRCA得到了一个峰值为(169.14±14.95)b/min的平均ITR，而DS策略得到的平均ITR为(230.2±65.8) b/min。配对t检验显示DS策略所得到的最高ITR显著优于整合TRCA(t₃₄=11.55, p＜0.001)。

表 1显示的为FS策略、DS策略以及理论值的对比。如表 1所示，DS策略比FS策略有更高的平均ITR。配对t检验显示，对于ITR和准确率，DS策略要显著优于FS策略(ITR: t₃₄=8.3, p＜0.001；准确率: t₃₄=4.7, p＜0.001)。但是，对于所用到的数据长度，两种方法在统计上没有显著差异(t₃₄=-1.7, p=0.1)。另外，DS策略的结果，包括数据长度、准确率和ITR，都显著差于理论值(所有的情况p＜0.001)。

2.2 ITR提升分析

图 3展示的是DS策略比对应的FS策略的ITR提升情况。可以发现，不同的被试有着不同的ITR提升比例，所有被试提升比例的均值是13.5%，其中提升比例最高的是第12名被试，达到40.8%。并且，使用DS策略，所有被试的ITR比使用FS策略的都没有出现负增长。这些结果表明, 使用TRCA分类模型的DS策略能够进一步提升被试的BCI表现。

2.3 输出过程分析

图 4显示的是在模拟在线测试期间使用DS策略在不同时刻的累计输出比例。为了对比，理论值的输出情况也显示了出来。对于正确的输出比例，在各个时刻，理论情况的输出都要显著高于DS策略(所有的情况p＜0.001)。但是，DS策略对应的输出比例有更快的增长趋势，而且在1 s处实现了和理论值相当的输出比例。图 4b所示为不同方法的错误输出比例。在理论情况下，1 s以前没有错误输出；在1 s处，使用DS策略的错误输出显著多于理论情况(t₃₄=5.2, p＜0.001)。

3 讨论

关于DS策略的研究主要集中在基于P300的BCI系统中，将其引入到SSVEP-BCI的相关研究相对较少。自适应的方法在早期被用来根据被试当前的在线实时表现，动态地限制子试次呈现的数目，从而减少刺激呈现的时间，进而达到优化ITR和准确率的目的^[18-20]。随后，Mainsah等提出了基于Bayes推断的用于P300拼写器的数据驱动算法，它基于用户EEG数据的信噪比，自适应地选择子试次的数量，从而提高拼写速度^[21]。最近，Aydin等提出了动态自适应算法来确定所需要的刺激序列的数目，并用受试者工作特征(receiver operating characteristics, ROC)曲线的最优工作点来确定阈值，结果显示动态自适应策略能够显著提升系统性能，并降低刺激序列的数目和出现错误的比例^[22]。至于DS策略在SSVEP-BCI上的应用，Yin等提出了一个动态优化刺激时间以提升系统性能的方法^[6]。在本研究中，基于Bayes估计的DS策略被整合到SSVEP-BCI中，除了以累加的相关系数作为特征，还使用了依赖于时间的信息来计算分类结果的置信度。另外，有3名被试的数据在使用DS策略处理之后所得到的ITR，与使用FS策略的情况相比并没有提高(见图 3中的第2、8和23名被试)，这可能是由于相对较低的分类准确率并没有能提供足够的信息用来评估分类的置信度。

4 结论

本文提出的Bayes动态自适应策略能够显著提升SSVEP-BCI的性能。该方法通过建立关于EEG特征值的正确预测和错误预测概率分布模型，结合相应的先验概率分布，从而利用Bayes推论估计未知EEG片段的正确预测概率，并与预设的阈值进行比较，以决定是否在当前时刻输出预测的分类结果，最终实现了自适应地为每个个体的每个试次优化数据长度的目的。结果显示：相比于传统的FS策略，本文方法可以在几乎不损失分类准确率的基础上，缩短了所用到的平均数据长度，从而得到更高的ITR。未来的研究将聚焦于将动态自适应策略整合到实际的在线SSVEP-BCI系统中。