脑机接口(brain computer interface，BCI)作为新一代人机交互技术，可实现人脑与计算机或其他设备之间的直接通信[1]．脑电信号可以通过视觉、听觉等不同诱发形式产生．当受到固定频率的视觉刺激时，大脑皮层产生连续的响应称为稳态视觉诱发电位稳态视觉诱发电位(steady state visual evoked potential，SSVEP)[2]．基于SSVEP的BCI系统(SSVEP-BCI)是目前BCI领域研究的热点之一[3]．基于SSVEP的BCI系统须设计刺激界面来诱发脑电信号．已有研究表明：刺激界面的众多变量，如形状、频率、颜色和分散程度[4-5]，都能够显著影响SSVEP的响应强度，进而决定脑机接口的性能．目前，主要以实验结果判断刺激界面的优劣评估给定刺激界面的刺激效果．传统的刺激界面优化设计建立在使用DOE方法安排实验采样的基础上，存在实验次数多、成本高、周期长的问题[6]，主要原因是：a. 已有的DOE方法中实验方案制定和刺激界面方案优化独立进行，导致出现较多的无效实验点[7]；b. 诱发SSVEP的刺激变量类型众多、取值范围各异，若要完整测试所有刺激界面变量组合，则需大量实验且使被试者产生认知和视觉负荷[8]；c. 在SSVEP-BCI技术尚未成熟的背景下，可以在真实SSVEP-BCI系统上执行的实验数量有限．总之，以更少的实验获取最优的刺激界面方案是进行SSVEP刺激界面的实验设计的关键．实验设计的方法可划分为空间填充设计和序贯式实验设计．空间填充设计先对样本空间进行全局性建模，再开展实验方案，主要有正交设计[9]、拉丁超立方实验设计[10]等．序贯式实验设计通过序贯地增加实验样本，对实验方案逐步优化达到最优解，其中有黄金分割法、贝叶斯优化[11](Bayesian optimization，BO)等．现有的序贯式实验设计能够有效结合实验方案设计和目标优化过程，充分利用实验过程信息．然而，在序贯式实验的初始阶段仍存在采样次数较多的问题．针对这一问题，最新的研究思路是在实验设计中融入先验知识．文献[12]证明了在计算机仿真实验上集成先验知识的有效性，提出了一种结合仿真数据和形状知识的导数信息进行贝叶斯优化的方法．基于以上研究成果，SSVEP实验设计中一般存在刺激界面变量对SSVEP响应影响的共性知识，借助此类知识有目的地安排实验点，有利于提高实验效率，也能够显著降低实验成本．本研究提出一种知识驱动的SSVEP序贯实验方案优化技术．首先利用概率模型表征SSVEP刺激特性的先验知识．利用概率积分变换对样本空间进行翘曲，将较好刺激界面方案存在可能性高的区域进行扩充，同时压缩可能性低的区域．在此基础上，建立序贯式实验设计的优化模型，对该模型进行基于贝叶斯优化方法的实验设计方案优化．为了验证所提出方法的有效性，采用文献[5]中的SSVEP实验设计问题及数据，分别利用正交实验设计、拉丁超立方实验设计和所提方法开展实验设计，对比分析所需的实验次数．1 问题描述在SSVEP刺激界面的实验方案设计中，将刺激界面的变量作为自变量，SSVEP信号响应作为判别刺激界面优劣的标准．SSVEP刺激界面的实验中包含两类输入，分为可控输入xi和不可控输入δ．可控输入是在d个刺激变量中选取并组合成n组刺激界面，不可控输入因素δ包括环境亮度、噪音等．该过程的最终输出结果为可观测的输出结果f(xi)和不可观测的被试者生理负荷ωi．SSVEP刺激界面实验方案的黑箱问题表示如图1所示．10.13245/j.hust.210720.F001图1SSVEP刺激界面实验方案的黑箱问题表示只考虑可控输入因素对可观测的输出结果的影响．设SSVEP的刺激变量是一组d维的可观测设计变量x，即x={x1,x2,…,xd}，实验变量组成的样本空间为x∈Rd．通过实验设计方法寻找最优的SSVEP刺激界面x*这一实际问题，可以表述为一个输出为SSVEP响应最大化的优化问题x*=argmax             x∈Rdf(x), (1)式中f(x)为目标函数，它通常由单一的SSVEP响应性能指标表征，如SSVEP分类精度、响应幅值、能量效率、识别准确率或信息传输速率等．由于刺激界面变量和SSVEP响应指标之间的映射关系未知，而且很难获取解析表达式，因此SSVEP的刺激界面优化为黑箱问题，其特性如下．a. 黑箱特性：在实验寻优的过程中，只能观测到实验的输入与输出，不能获取过程中的梯度信息．b. 随机特性：由于存在不可控的输入因素和随机性的实验采样，因此相同条件下每次实验的SSVEP响应值有所差异．c. 全局优化：整个实验样本空间上存在多极值情况，不能假设目标函数f的局部极大值个数或凸性，理想情况下须寻求目标函数的全局最大值．要求解表征该实验设计黑箱函数的优化问题，常用方法是利用有限的样本数据构建代理模型[13]．本研究基于贝叶斯优化构建高斯过程代理模型，进而进行序贯式实验设计．2 知识驱动的序贯式实验设计方法传统DOE实验设计仅从样本空间中指定部分实验点，无法有效获取最优变量组合，存在一定盲目性．这里使用边设计边实验的序贯式实验设计方法，以知识和数据为指导逐步地开展实验设计，使得获取的刺激变量组合对SSVEP响应具有最大的贡献．知识驱动的SSVEP刺激界面实验方案优化方法如图2所示．10.13245/j.hust.210720.F002图2知识驱动的SSVEP刺激界面实验方案优化方法2.1　嵌入先验知识的高斯过程代理模型采用贝叶斯优化进行SSVEP的实验设计．贝叶斯优化[11]基于贝叶斯定理，利用已知有限的SSVEP样本数据对实验设计的黑箱函数进行全局优化p(f |Di)∝p(Di|f)p(f)，(2)式中：p(f)与p(Di|f)分别为f(x)的先验和后验概率分布；p(f |Di)为似然函数；Di为已观测样本点集合，满足递归关系：D1={(x1,f(x1))};Di={Di-1,(xi,f(xi))}．(3)贝叶斯优化的代理模型采用高斯过程．高斯过程由一个均值函数m：X→R和一个半正定的协方差函数k(x,x')组成．该模型是表征SSVEP刺激界面的变量与响应之间映射关系的关键：f(x)∼GP(m(x),k(x,x')); (4)m(x)=E[f(x)],k(x,x')=E[(f(x)-m(x))(f(x')-m(x'))]. (5)采用基准的平方指数协方差函数k(x,x')=σ02exp[-||x-x'||2/(2l2)]，(6)式中：l为与函数平滑度相关的长度尺度变量；σ02是函数方差．为了利用概率模型表征的先验知识来指导实验采样过程，采用概率积分变化定义一个翘曲空间，并通过修改协方差函数将翘曲融入到高斯过程中．设v为服从连续分布的刺激变量随机变量，p(v)为其概率密度函数(probability density function，PDF)，Φ为其累积分布函数(cumulative distribution function，CDF)，则任意随机变量u=Φ(v)遵循均匀分布U(0,1)．概率积分变换为Φu(u)=p(u≤u)=p(Φv(v)≤u)=  p(v≤Φv-1(u))=Φv(Φv-1(u))=u．(7)高斯分布的概率密度函数与累计分布函数曲线如图3所示．图中v∈[3,7]的区间内占样本空间的36.4%，目标函数值不佳的样本区域称为低似然区域，v∈[0,3]⋃[7,11]的区域内占样本空间的63.6%．经过概率积分变换后的累计分布函数，当v∈[3,7]时，Φ∈[0.05,0.84]，此时该高似然区域占新样本空间的79%；而v∈[0,3]⋃[7,11]时，Φ∈[0.00,0.05]⋃[0.84,1.00]，此时低似然区域仅占新样本空间的21%．可见，对样本空间的翘曲使原始样本空间的最优值区域被扩大，其他区域被缩减．10.13245/j.hust.210720.F003图3高斯分布的概率密度函数与累计分布函数曲线将包含最优刺激变量的知识整合到协方差函数中，引入高斯过程模型．将Φ(x)替换协方差函数中的x，使协方差函数中变量之间的距离根据知识强度放缩，新的协方差函数形式[14]为k˜a(xa,xa')=σ02exp-||Φa(xa)-Φa(xa')||22l2;k˜(x,x')=∏i=1dk˜a(xa,xa'), (8)式中：d为输入空间维数；Φa(xa)和Φa(xa')分别为点xa和xa'沿维数a的CDF变换．式(2)~(8)给出了建立并更新贝叶斯优化中基于翘曲协方差函数的高斯过程代理模型，使得后续采样倾向于从高似然区域产生数量更多的实验样本．通过融入知识，可在一定程度上避免实验在低似然区域内的低价值探索，有效缩减了采集无效样本点的次数．2.2　基于EI采集函数的采样策略采集函数的作用是在高斯过程建模SSVEP的变量与目标函数值后，通过权衡选择下一组实验样本点，进行下一次实验评估．使用EI采集函数E(x)=(max(f(x)-f(x+),0))，式中：f(x+)为当前最优样本的函数值；x+为当前最优样本点，即x+=argmax             xi∈x1  : tf(xi)．分析评估EI效果的函数定义为Z=(μ(x)-f(x+)-ξ)/σ(x)    ( σ(x) 0), 0    (   σ(x) =0);E(x)=(μ(x)-f(x+)-ξ)Φ(Z)+σ(x)P(Z)，式中：μ(x)和σ(x)分别为高斯过程后验分布的均值和标准差；Φ(x)和P(x)分别为代理模型中样本对应的累积分布函数与概率密度函数；ξ用于确定采集函数的探索量．3 SSVEP刺激界面实验设计验证3.1　研究对象与参数设置以式(1)中的模型为待优化模型，贝叶斯优化基于python的BayesianOptimization模块．实验中选取文献[5]中的SSVEP实验方案的数据源进行实验设计．该数据中包含20名被试者数据，包含枕顶区的O1，O2，Pz，P3，P4，P7和P8七个通道，实际计算时经统计学检验并取测量均值，意味着该算例的最优界面是针对全体被试者的群体最优界面．主要刺激界面变量包括刺激频率、颜色、尺寸、间距和形状，表1设定了此次实验每个刺激变量的取值范围(颜色取值已映射到0~255的连续值空间)，用于定义实验样本空间．10.13245/j.hust.210720.T001表1SSVEP刺激界面实验刺激变量设定刺激变量变量范围先验分布先验分布参数刺激频率x1∈[0，25]高斯分布μ=15.0，σ=4.0刺激颜色x2∈[0，255]高斯分布μ=255.0，σ=90刺激块边长x3∈[0°，10°]高斯分布μ=5.5，σ=3.0刺激间距x4∈[0°，10°]高斯分布μ=5.0，σ=2.5刺激形状x5∈{方，圆}根据已有的文献或专家知识，采用概率模型对先验知识进行量化表征，建立先验分布．由于刺激形状为无法以数值度量的离散标称值，因此并不指定其先验分布，根据类别独立建立样本空间．采用高斯分布建立刺激变量的先验知识分布，其参数设定见表1．3.2　知识驱动方法的有效性验证将传统的贝叶斯优化方法与知识驱动的方法进行对比．知识驱动的方法将每一维度的先验分布融入翘曲协方差函数中，建立高斯过程代理模型．传统的贝叶斯优化方法使用式(6)的基准平方指数协方差函数．SSVEP响应在很大程度上取决于在一定刺激频率下的SSVEP信号功率与环境信号功率之间的信噪比．使用刺激频率(ω)为14和25 Hz时Oz单通道的信噪比，随机采样5个实验数据作为初始样本建立代理模型，共采样迭代35次趋于平稳，测定两种方法每次迭代的SSVEP信噪比(r)随迭代次数(m)的变化曲线如图4所示．10.13245/j.hust.210720.F004图4两种方法每次迭代的SSVEP信噪比随迭代次数的变化曲线由图4可知：当迭代次数小于15时，两种方法的信噪比波动都较为明显，但在15~35范围内都平稳上升．这是由于初始样本点较少时，EI采集函数须权衡未探测区域和已知较优区域，因此在迭代初期尤其是频率为25 Hz时，两种方法的差别不够明显；迭代中后期掌握了一定的样本空间信息，会使采样偏向已知较优区域．总体上看，知识驱动的优化方法达到最优信噪比的迭代次数更少．在SSVEP-BCI系统方面，SSVEP的输出性能度量常使用分类精度和信息传输速率(ITR)[15]．以分类精度为响应目标，终止条件为迭代19次，SSVEP分类精度(q)随迭代次数变化曲线见5(a)；以信息传输速率为响应目标，终止条件为迭代24次，信息传输速率(V)随迭代次数变化曲线见5(b)．从图5可以看出：采用知识驱动的优化方法的SSVEP精度和ITR均值曲线随迭代次数增加明显上升，不确定性逐渐减少．SSVEP的分类精度达90%时，传统贝叶斯优化方法需要18次，而知识驱动的优化方法仅需13次；ITR达28 bit/min时，传统贝叶斯优化方法需要20次，而知识驱动的优化方法仅需15次．这表明达到相同的SSVEP响应时，知识驱动的优化方法通过先验知识的翘曲融入，有助于采集函数更快找到最佳实验方案，减少了采样次数．10.13245/j.hust.210720.F005图5两种方法的迭代过程对比3.3　实验设计方法对比分析根据文献[5]中的刺激变量选择过程，设定两个阶段的SSVEP刺激变量的分析实验用于确定最优变量，实验1中的刺激变量的因子水平如表2所示，在实验1的基础上选取较优的变量值并应用于其他20个被试者设计实验2，实验2中的刺激变量的因子水平如表3所示．10.13245/j.hust.210720.T002表2实验1因子水平表刺激变量(因子)刺激水平一二三四五刺激频率/Hz14172530刺激颜色白(255)红(125)绿(106)黄(246)蓝(66)刺激块边长/(°)0.571.492.603.72刺激间隔/(°)02.302.605.27刺激形状方形圆形10.13245/j.hust.210720.T003表3实验2因子水平表刺激变量(因子)刺激水平一二三刺激颜色白(255)红(125)黄(246)刺激块边长/(°)2.603.72刺激间隔/(°)2.302.60SSVEP的输出性能度量采用事件相关扰动指标，即通过刺激频率处或其谐波处的功率相对增加测量SSVEP强度，ERSP(event related spectral perturbation，事件相关谱扰动)反映特定频带内由刺激引起的频谱功率变化，定义如下[5]Pω=(P¯ω+-P¯ω-)/P¯ω-,式中Pω+和Pω-分别为给定刺激频率ω下在视觉刺激期间和非视觉刺激期间记录的平均信号功率．将正交实验设计方法、拉丁超立方实验设计方法和知识驱动的优化设计方法，应用于该文献的数据样本空间对SSVEP刺激界面开展优化．三种实验设计方法的条件设定如下：正交实验设计中考虑刺激尺寸与间隔的交互效应；拉丁超立方实验设计中由专家根据经验设定样本子空间个数；知识驱动的优化设计只利用该实验的变量范围，并不指定每个变量的离散值．根据文献[5]中达到最优刺激界面的SSVEP响应标准，设定贝叶斯优化算法的终止条件为P=10．得到如表4所示的三种方法的实验次数对比．10.13245/j.hust.210720.T004表4三种方法的实验次数对比方法第一阶段第二阶段总计正交实验设计26632拉丁超立方实验设计20525知识驱动的优化设计51015正交实验设计和拉丁超立方实验设计方法都须将样本空间离散化，因此对于连续刺激变量的描述能力较差，很难达到全局最优解．正交实验设计方法所需的实验次数过多，拉丁超立方方法的实验次数相比于正交实验设计次数有所降低，但针对周期较长的SSVE-BCI交互实验仍会导致被试者疲劳．此外，该方法在设计效率和样本性能方面很难均衡，须凭借经验或算法设定样本子空间．对于知识驱动的优化设计方法，其实验次数对比正交实验设计和拉丁超立方实验设计分别降低了53%和44%，优化效果显著．同时，该方法能够保证每个连续的实验点都有概率被采集函数采样评估，从而实现寻找最优刺激变量的目的．4 结论针对SSVEP刺激界面变量优化复杂和交互实验中的被试负荷问题，提出以知识驱动的贝叶斯优化方法进行序贯实验设计，有效减少实验次数，显著提升实验效率．a. 分析了SSVEP实验设计的相关刺激变量并建立实验设计问题的黑箱函数模型；建立了实验变量的先验知识分布模型，采用概率积分变换定义了翘曲的样本空间，按照变量对函数的决策度量分配空间．b. 在贝叶斯优化框架下修改高斯过程的协方差函数，使得先验知识在翘曲样本空间中融入实验设计的代理模型，并通过采集函数生成新的实验评价点．c. 通过传统贝叶斯优化和知识驱动的优化方法的迭代对比，验证所提方法有效性；对比了知识驱动的优化实验设计、正交实验设计和拉丁超立方实验设计，评估其应用于SSVEP刺激界面求解最优值所用的实验次数，验证了该方法的优化效果．