一种基于运动想象的脑-机接口时空滤波器迭代算法

段 放 高小榕

(清华大学医学院生物医学工程系，北京 100084)

引言

脑-机接口(brain-computer interface，BCI)是一种不经过外围肌肉运动，通过直接提取大脑活动特征来指挥外部设备的通信方式［1］。在脑-机接口系统中所使用的监测大脑活动的脑功能成像方法多为头皮脑电图(electroencephalogram，EEG)信号，其具有对使用者无损伤的优势。

运动想象是在肢体不实际运动的情况下对肢体某个部分的潜在想象［2］。基于运动想象的脑-机接口系统的生理依据是使用者在想象不同肢体运动时，自主控制其大脑运动感觉皮层与所想象肢体对应区域μ节律或β节律的变化。如何提取μ节律或β节律的能量差异这种特征是运动想象脑-机接口系统的关键问题。

共空间模式(common spatial patterns，CSP)已被证实是运动想象脑-机接口系统特征提取的有效方法［3］，其生成的空域滤波器可有效提高信号的信噪比。但CSP只能在信号进行了比较合理的时域滤波条件下，才能有效地发挥其性能［3］。因此通常在使用CSP前对信号的时域带通滤波变得十分关键。一般情况下，常采用性能不佳的较宽通带滤波(如7～30 Hz)，或依据受试脑电信号特点对滤波器人为离线设计，这必然限制了系统的性能提升。

为了避免过多离线人工干预，很多机器学习方法被引入到CSP的改进中来。例如共空频模式(common spatial spectral pattern，CSSP)、稀疏空频模式(common sparse spectral spatial pattern，CSSSP)、滤波器库空域模式(filter bank common spatial pattern，FBCSP)等［4－8］。这些方法都通过一定程度的遍历搜索确定出脑电信号的有效频率带宽，但由于涉及遍历搜索计算量均较大或优化求解不确定性高等问题，很难在线实现。

本研究提出一种分别对空域滤波器和时域滤波器进行迭代优化的特征提取方法，避免了遍历搜索的时间不确定性，并给出了收敛的必然性证明。减少了特征提取中的人工干预，有利于在线应用。

1 计算方法

1.1 共空间模式

对于两类多导信号xi(t)∈RC，C为导联数，i=1，2。信号的协方差矩阵为:

l为两类信号的长度。CSP生成一组空间加权系数w∈RC，使经过空域加权滤波后的信号yi(t)=wTxi(t)更具有辨识力。其信号能量为:

所谓辨识力体现在分别最大化两者能量比E1/E2或E2/E1。不失一般性，设为:

这个优化问题可通过求解广义特征值问题Σ1w=λΣ2w获得结果，其中λ即为两信号能量比。限于篇幅，具体的推导过程可参考文献［3］。

由于通常CSP只在经过特定频率滤波的情况下才能发挥其有效性，需分别引入一个给定长度的有限冲激响应滤波器(Finite Impulse Response filter，FIR filter)H。确定一个有效的H是很关键的，Dornhege等人也提出了使用 FIR滤波器与 CSP结合最大化能量比的思路(CSSSP)。但他没有给出FIR滤波器的求解方法，而采用线性梯度搜索的方法通过穷举的方式确定FIR滤波器。即使对长度较短的 FIR 滤波器也无疑计算量是巨大的［5，7］，因此不适用于在线系统。而较短的滤波器也无法提供相对窄的滤波频带，CSP在较宽的滤波频带下性能是有限的，因此CSSSP并不能得到足够好的效果。

1.2 共迭代时空模式

给出一种分别整合空域滤波器、FIR滤波器的机器学习方法来提高CSP的性能。原CSP问题为两个形式相同的最大化能量比问题。

1.2.1 问题提出

图8为母线电流响应法的实验波形。首先,对开关组合S1和S2,S5和S6分别进行开通和关断,测得的电流响应分别为i3和i2,如图8(a)所示;然后,对开关组合S3和S4,S5和S6分别进行开通和关断,测得的电流响应分别是i1和i2,如图8(b)所示。

对于给定长度的 FIR滤波器h(n)，n=0，…，L－ 1，向前补0l－ L项得到 H=h(n)，n=0，…，l－1，H ∈ Rl。则其滤波器循环行列式为 A，A ∈ Rl×l，其中:

经过 FIR滤波器滤波后的信号 zi=xiA。原CSP问题为最大化能量比问题表述为:

将滤波前的信号带入式(1)得到:

式(2)为ICSTP所想解决的优化问题。

1.2.2 优化求解

在机器学习中最优解往往会引起过拟合的不稳定问题，所以可以不追求全局最优解。下面给出利用基于迭代的局部最优解求解方法。

首先尝试在某一个初始FIR滤波器H的基础上优化空域滤波器，式(2)转化为

这是一个原始的CSP问题，可通过解广义特征值问题求得空域滤波器 w。对于原始信号 xi(t)，经过空域滤波后的信号为yi(t)=wxi(t)，代入式(2)得到:

得 yiA=HRi，其中 Ri∈ Rl×l为 yi的逆向循环矩阵。

对于长度为 L的滤波器，H中前 l－L项为0。现将 H长度还原为，只需考虑 Ri的后行，即 Qi∈RL×l，式(2)优化转化为

可以发现经时空滤波器互换以上优化问题与传统CSP优化问题形式相同，因此次优化可使用优化传统CSP方法的求解广义特征值问题的方法求得新的滤波器H。将新的H变换为A代入式(3)再求解空域滤波器w。

通过式(3)与式(4)交替迭代求得最终需要得到的FIR滤波器H与空域滤波器w。出于脑电信号的特点，我们可以使用任意7～30 Hz FIR滤波器作为初始滤波器。迭代次数的选择可以通过以下两点由使用者直接确定:即到达预先设置的迭代次数或相邻两次迭代的空域滤波器w的变化很小。对于运动想象问题有两个能量比最大化问题E1/E2或E2/E1，分别对应于增强左手或右手想象信号能量。因此我们最终根据两个问题分别返回1组对应的空域滤波器与 FIR 滤波器解 {wi，Hi}，i=1，2。整体算法流程见图1。

1.2.3 算法的收敛性

作为迭代算法，其收敛性是非常重要的。下面给出这种迭代求解方法收敛性的证明。

每次对式(3)的优化是在式(4)的优化所得结果H基础上进行的，而对式(4)的优化则是在式(3)的优化所得结果w基础上进行的。每次优化均极大化两信号能量比，易知对于第i次迭代有

2 实际脑电信号评估

图1 算法流程Fig.1 Diagram of algorithm

为检验ICSTP的运动想象特征提取能力，使用人工选择滤波频段的CSP算法进行特征提取作为比较对象，评价ICSTP算法的性能。

2.1 人工进行频段选择

CSP算法的目标是使数据更有辨识力，因此要选择有足够辨识力的频段作为CSP方法的预先滤波频段。利用r2值来评价两类数据在频域上是否有辨识力［1］，使用r2较高的频段作为 CSP的预先滤波频段。

2.2 特征分类器

模式识别中特征提取与特征分类是两个连贯的过程。本算法作为一个特征提取的方法，仍需要在特征提取后应用相对应的特征分类方法。对运动想象的脑电数据分类中线性分类器具有足够好的性能［10］，因此使用线性软间隔支持向量机(Support Vector Machine，SVM)。SVM 分类器源于统计学习原理，由于其基于结构风险最小化的特点使其有相对强的抗过拟合能力［11］，这样较有利于小样本学习问题。之所以使用线性核的SVM，是为了防止样本学习在高维非线性空间产生的过拟合现象。

2.3 脑电实验数据

图2 数据处理流程Fig.2 Diagram of data processing

用以检验算法的数据是本实验室采集的运动想象脑电数据。7组脑电数据分别从年龄在21～24岁之间的7位志愿者处采集，数据是在有视觉反馈提示的情况下获得的。每次试验时间长度为8 s。在2 s的无提示准备时间后屏幕上出现方向提示，提示时间1 s。之后受试通过左手或右手的运动想象控制屏幕上光标根据提示进行垂直运动。运动想象时间持续5 s。5 s想象后根据光标，最终位置给出正确与否的提示。每个受试分两天共记录240次试验数据，其中左手想象右手想象各120次用于离线分析。整个实验采用 BioSemi ActiveTwo脑电放大器系统记录数据。根据国际10/20系统，在第一运动区与辅助运动区设置了32个电极。参考电极为左耳垂，采样频率256 Hz。用以离线分析的数据的时间窗选定为试验的3.5～7.5 s。为CSP算法与ICSTP算法的公平性使用全部32导数据用以检验算法而不再各自单独挑选导联。更具体实验时序图和脑电电极安放的实际位置图可参考文献［7］。

2.4 ICSTP设置与迭代条件

进行一次时空优化作为一次迭代。在整个数据实验过程设定10次迭代为最终迭代次数。同时为检验收敛速度，对于实验所用32导数据，设定基于空域滤波器的迭代停止条件运动想象人为设计滤波器带宽通常小于5 Hz［12］。为获得小于5 Hz的通带宽度，设定FIR滤波器约为1/5采样频率，长度为50点。

2.5 数据实验结果

对数据进行了10×3随机分组划分实验。即将两类数据平均分为相等的3份，每份中两类数据数量相等。每次测试将其中两份作为训练集，最后一份作为测试集。同时给出宽带滤波(7～35 Hz)和人工挑选频段的后CSP结果，说明频段选择的重要性。作为一种相似的迭代算法，给出文献［7］中对相同数据使用ISSPL算法的实验结果作比较。分类正确率与迭代次数见表1，数据处理流程见图2。

表1中黑体字部分为正确率最高的情况。可以看到仅有两名受试在人工选择滤波频段后进行CSP的结果好于ICSTP。而CSP优于ICSTP的情况是在受试运动想象能力较强的情况下，且两者结果差距很小。而对于两名想象运动能力相对较弱的受试，ICSTP的结果明显优于 CSP较多。同时本算法与ISSPL相比仍保持向对高的正确率。

表1 分类正确率±方差与迭代次数±方差Tab.1 Classification accuracies±standard deviation and iteration times±standard deviation

3 讨论

可以发现平均迭代次数均不超过4次。所有的数据实验过程均未达到预先设定的迭代10次的迭代次数最大值。在很少的迭代次数之后空域滤波器已无明显变化，算法收敛速度快。

3.1 与ISSPL相比较的优势

同样作为迭代算法，本算法在正确率上也普遍高于ISSPL。这可能是由于本算法在收敛性上能够得到保证。ISSPL中的将分类器设计与频率窗的学习相结合，采用支持向量机计算频率窗函数同时确定决策函数分类面，其空域滤波器的设计采用CSP最大能量比的思想，频域滤波器期望取得尽量大的间隔。迭代算法的两个迭代步骤没有统一的优化目标，这样算法的收敛性难以证明。本算法在两个步骤均最大化能量比，已证明了算法收敛。本算法与ISSPL还有一个不同点是本算法首先进行特征提取，将特征维度压缩到2维后再进行分类器学习。而ISSPL在学习频率窗函数的同时进行分类函数学习。分类决策函数所需要确定的系数与频率窗长度相同，维数高，相比本算法可能容易出现过拟合。

除正确率外，算法的计算时间也非常重要。ISSPL算法每次进行频域滤波器计算需求解一个二次规划问题。此二次规划问题的维度为2×频率窗长 +样本个数。而本算法只需解决一个广义特征值问题，矩阵的阶次为FIR滤波器长度。消耗时间较少。同时ISSPL算法中的正则化项的数值也需要通过交叉验证来确认。对于迭代算法每次迭代过程的正则化项是否应该使用相同的值也仍然是一个开放的问题。使用一台配置为Inter Core2 2.8 GHz CPU、4 GB内存的普通个人电脑考察算法计算时间。文献［5］中所描述CSSSP算法中单次优化求解需要约1 min，由于CSSSP需要交叉验证，最终约需要20 min。对于ISSPL算法，忽略算法的正则化系数的选取时间，仅在确定一个系数的条件下进行ISSPL算法运算。5次迭代使用时间约为80 s。本算法经过5次迭代选取特征后在设计分类器，总耗时约15 s。可见本算法更有利于在线系统。

3.2 特征提取效果

图3 频域特征。(a)各导联频域 分布图;(b)左手想象特征滤波器的幅频响应;(c)右手想象特征滤波器幅频响应Fig.3 Feature in frequency domain.(a)r2-value of all channels in frequency domain;(b)magnitude response of left-hand movement imagery filter;(c)magnitude response of right-hand movement imagery filter

观察图3(a)为受试WW r2分布图。可以发现其脑电在13 Hz与14.5 Hz的r2值出现两个峰值。在经过ICSTP迭代学习后知道了出现这种现象的原因。即这个受试的在进行左侧想象和右侧想象时对应的μ节律是不同的。比较图3中(b)与(c)，可以发现ICSTP方法计算出的FIR滤波器很好地吻合了r2分布图的结果。强化左侧想象的FIR滤波器与强化右侧想象的FIR滤波器分别与r2分布图中两个辨识度最高的频段相吻合。

4 结论

本研究提出了一种基于迭代优化的运动想象特征提取方法。并且数据实验证实了算法的有效性。这种算法能在无人工干预的情况下，快速计算出适用于不同受试的运动想象时域滤波器与空域滤波器。作为一种迭代算法，并给出了迭代收敛性的证明。并在数据实验中证实了其收敛的快速性。与传统的运动想象空频滤波器特征提取方法相比有着计算速度快的优势。这无疑有利于在线系统的实现，特别是应用于在线系统的在线更新问题上。

(致谢 感谢本实验室吴畏与闫铮提供的重要建议。感谢本实验室王毅军提供了研究数据)

［1］Wolpaw JR， BirbaumerN， McFarlandDJ， etal.Braincomputer interfaces for communication and control［J］.Clinical Neurophysiology，2002，113(6):767－791.

［2］De Vries S，Mulder T.Motor imagery and stroke rehabilitation:a critical discussion ［J］.Journal of Rehabilitation Medicine，2007，39(1):5 －13.

［3］Blankertz B，Tomioka R，Lemm R，et al.Optimizing spatial filters for robust EEG single-trial analysis［J］.IEEE Signal Proc Magzine，2008，25(1):41－56.

［4］Lemm S，Blankertz B，Curio G，et al.Spatio-spectral filters for improved classification of single trial EEG ［J］.IEEE Trans Biomed Eng，2005，52(9):1541－1548.

［5］Dornhege G，Blankertz B，Krauledat M，et al.Combined optimization of spatial and temporal filters for improving brain–computer interfacing［J］.IEEE Trans Biomed Eng，2006，53(11):2274－2281.

［6］Novi Q，Guan C，Dat TH，et al.Sub-band Common Spatial Pattern (SBCSP)forBrain-ComputerInterface［C］//Proceedings of 3rd International IEEE/EMBS Conference on Neural Engineering.Hawaii:IEEE，2007:296－300.

［7］Wu Wei，Gao Xiaorong，Hong Bo，et al.Classifying Single-Trial EEG During Motor Imagery by Iterative Spatio-Spectral Patterns Learning(ISSPL)［J］.IEEE Trans Biomed Eng，2008，55(6):1733－1743.

［8］Ang KK，Chin ZY，and Guan C.Filter Bank Common Spatial Pattern(FBCSP)in Brain-Computer Interface［C］//Proceedings of IJCNN 2008［C］.HongKong:IEEE，2008:2391－2398.

［9］Fukunaga F and Koontz W.Applications of the Karhunen-Loève expansion to feature selection and ordering［J］.IEEE Trans Computers，1970，19(5):311－318.

［10］Dornhege G，Blankertz B，Curio G，et al.Boosting bit rates in noninvasive EEG single-trial classifications by feature combination and multiclass paradigms［J］.IEEE Trans Biomed Eng，2004，51(6):993－1002.

［11］Vapnik VN.Statistical Learning Theory［M］.New York:Wiley，1998.

［12］Pfurtscheller G，Lopes Da Silva FH.Event-related EEG/MEG synchronization and desynchronization:basic principles［J］.Clin Neurophysiol，1999，110(10):1842－1857.