脑波信号是伴随着生物电信号而产生的一类信号，承载了人脑活动的很多信息。脑波信号采集方法的发展以及信号处理与分析技术的进步，尤其是近年来机器学习和模式识别技术的发展，使得脑波信号的分析成为研究热点^[1-2]。

P300脑电信号是目前人们研究最广泛的脑电信号之一。当满足在规律的视觉刺激中出现小概率刺激的实验范式时，人脑会产生该响应信号。其中，大概率刺激是一般性刺激，小概率刺激是靶刺激。这种实验范式被称为oddball范式。P300信号在时域波形上有明显的特征，在刺激开始之后约300 ms处有一个正峰^[3]。

利用这种oddball实验范式，可以设计如打字机之类的应用^[4]以实现脑机接口(brain-computer interface, BCI)。图 1所示的是一个6×6的字符矩阵，包含6行和6列。在实验时，这6行和6列会随机闪亮，每次闪亮视为对人眼的一个刺激，共12种刺激。当被试想打的字母所在的行和列闪亮时，为靶刺激，会诱发大脑产生P300信号，而其他行列闪亮时，是一般刺激，大脑不会有相应的反应。脑机接口系统通过接收脑电信号，处理后可以判断出被试所关注的行列，行列交叉即为要打的字。

脑电信号(electroencephalography，EEG)处理的一般流程包含了预处理、特征提取和分类3大步骤：1) 预处理阶段常用方法包括独立分量分析(independent component analysis, ICA)方法和带通滤波法。ICA^[5-6]分解脑电信号得到独立分量，以去除眼动、肌动等带来的生物电伪迹；带通滤波法滤除低频的电极漂移噪声和高频毛刺干扰。2) 特征提取和降维是为了解决EEG信号采样频率高、通道数多、数据维数高的问题，常用方法包括时域降采样法、主成分分析法(principe component analysis，PCA)、一致空间模式(common spatial pattern，CSP)等。时域降采样法保留了波形在时域上的变化趋势，适用于时域上有明显特征的脑电信号^[7]；主成分分析法通过提取脑电信号的主要分量达到降噪和降维的作用；一致空间模式通过寻找一组新的正交基使得原始信号在新的投影空间上具有更好的区分性^[8-9]。

脑电信号分类方法可以分为3大类：

1) 基于统计学习模型的方法。文[10]利用先验知识，事先选择了一组电极通道训练了多个支持向量机(support vector machine, SVM)分类器。文[11]在分类时使用Bootstrap方法，随机采样生成150个训练子集。每个子集用线性判别分析(linear discriminative analysis, LDA)训练分类器，最后投票决定分类结果。文[12]方法采用分类器级联的思路，将训练字符分成多组，对每组训练一个SVM分类器，再将多分类器的得分平均作为最终得分。

2) 基于神经网络的方法^[13-15]。文[13]建立4层的卷积神经网络(convolution neural network, CNN)，前2层分别是空间卷积层和时间卷积层，后2层是全连接层。文[14]利用块归一化技术解决了P300信号检测中的过拟合问题。文[15]利用自编码器自动学习不同个体间鲁棒的特征表达方式。

3) 基于稀疏编码的方法^[16-18]。文[16]通过求解脑电信号的稀疏表示(sparse representation)降低噪声干扰，并在稀疏空间上训练了一个鲁棒的分类器。

脑机接口系统对于脑电信号分类方法的准确性和实时性有着严格的要求。已有工作主要存在两大问题：1) 大多方法都是基于已有数据集进行离线的数据分析和处理。这些方法更加看重分类的准确性，却往往忽视计算复杂度，在实际情况下难以应用。2) 这些方法完全将脑电信号分类问题当作机器学习的分类问题来处理，所使用的分类模型如支持向量机、线性判别分析、卷积神经网络大多是鉴别式模型。这样处理虽然实现简单，可以套用许多经典鉴别式分类模型，但是却忽略了脑电信号本身的特性。脑电信号具有噪声强、维度高的特点，其分类问题不同于一般的分类问题。同时，P300脑电信号有着明显的时域波形特征，因此可以参考雷达信号分类方法，建立一个生成式模型。

1 算法框架与数据集介绍 1.1 算法框架

本文算法框架如图 2所示。首先，将脑电信号用带通滤波器作预处理，并用时域降采样法提取特征；接着，基于信号建模和最大似然估计求解P300正样本信号的主分量，并通过度量学习和梯度下降法迭代求解噪声信号方差，作为匹配滤波器的参数；最后，测试样本利用匹配滤波器模型推导出的线性判定算子得到最终分类结果。

1.2 数据集介绍

BCI国际竞赛Ⅲ中的数据集Ⅱ是当前研究中使用最广泛的一个数据集^[19]。它采集了2个被试的脑电数据，实验范式采用了图 1所示的6行6列的打字机。实验时，这6行和6列随机闪亮。闪亮持续时间为100 ms，之后变暗75 ms。每12次闪亮称为1轮，每个字符重复15轮。因此，对于每个字符，一共有15×12=180个刺激。这些刺激中，30个是正样本刺激，可以诱发P300响应；150个是负样本刺激。每名被试记录了85个字符的脑电波用于训练，100个字符用于测试。因此，训练集包含了2 550个正样本，12 750个负样本；测试集包含了3 000个正样本和15 000个负样本。脑电信号采集设备通道数为60，采样率是240 Hz。

2 P300脑电信号分类 2.1 通用匹配滤波器模型

受雷达信号检测模型启发^[20]，本文建立P300脑电信号分类的生成式模型。将正负样本2类情况建模如下：

$ \begin{array}{*{20}{l}} {{H_0}:x(n) = w(n),}&{n = 0,1, \cdots ,N - 1;}\\ {{H_1}:x(n) = s(n) + w(n),}&{n = 0,1, \cdots ,N - 1.} \end{array} $

(1)

其中：N是采样点数，x(n)是接收到的信号，s(n)是目标信号。w(n)是噪声，可以假设噪声服从Gauss分布，即w~N(0, C)，其中C是噪声协方差矩阵。

式(1)作了H₀和H₁两类假设：在H₀假设下，采集到的脑电信号x(n)是负样本，完全由噪声组成；在H₁假设下，采集到的脑电信号x(n)是正样本，由噪声和标准的P300信号模板两部分组成。

模型分类错误率为

$ {P_{\rm{e}}} = P({H_0}|{H_1})P({H_1}) + P({H_1}|{H_0})P({H_0})P({H_0}). $

(2)

其中：P(H_i)是H_i假设的先验概率，P(H₁|H₀)和P(H₀|H₁)分别是2类假设的错误分类概率。为了使模型分类错误率最低，需要做出使后验概率最大的决策，

$ P({H_1}|\mathit{\boldsymbol{x}}) > P({H_0}|\mathit{\boldsymbol{x}}). $

(3)

根据Bayes公式将后验概率转化为先验概率，

$ P({H_i}|\mathit{\boldsymbol{x}}) = \frac{{P(\mathit{\boldsymbol{x}}|{H_i}) \cdot P({H_i})}}{{P(\mathit{\boldsymbol{x}})}}, $

(4)

可以得到在最小Bayes风险准则下，判定为正样本的分类准则为

$ \frac{{P(\mathit{\boldsymbol{x}}|{H_1})}}{{P(\mathit{\boldsymbol{x}}|{H_0})}} > \frac{{P({H_0})}}{{P({H_1})}}. $

(5)

反之判定为负样本。

式(5)的左边是接收到的脑电信号x在2类假设下先验概率的比值，即似然比。式(5)右边为阈值。在P300打字机的应用中，6行和6列刺激对应的脑电信号中必然只有1个正样本、5个负样本。因此，在实际应用时，只需要求似然比的最大值即可。

根据噪声服从Gauss分布的假设，可以得到在H₀假设下，x~N(0, C)；在H₁假设下，x~N(s, C)，其中s是对P300波形主分量的估计。

定义对数似然比函数为

$ T(\mathit{\boldsymbol{x}}) = \log \frac{{P(\mathit{\boldsymbol{x}}|{H_1})}}{{P(\mathit{\boldsymbol{x}}|{H_0})}}. $

(6)

经过简化，可以得到匹配滤波器统计量，

$ T(\mathit{\boldsymbol{x}}) = {\mathit{\boldsymbol{x}}^{\rm{T}}}{\mathit{\boldsymbol{C}}^{ - 1}}\mathit{\boldsymbol{s}}. $

(7)

在P300打字机的应用中，先在6行中求取T(x)最大值，识别出目标字符所在行；再在6列中求取T(x)最大值，识别出目标字符所在列; 将行列交叉即可判断出目标字符。

式(7)给出了一个非常简洁的线性判别算子。对于测试样本x，只需要进行一次矩阵乘法就可以进行判断并给出分类结果，计算复杂度远低于其他已有的算法。

2.2 P300信号主分量估计

式(7)中的P300正样本主分量s是一个重要参数。由于正样本信号在时域波形上可能会存在混叠，因此不能通过时域平均的方法简单求解。图 3以一个三角波作为P300分量的示例。可以看出，左侧的脑电信号中存在多个P300分量的混叠。图 3左侧红色箭头指示了每个P300分量开始时刻。等号左侧为接收到的P300正样本信号X，等号右侧为指示矩阵D和P300信号主分量s的乘积。指示矩阵D是一个Toeplitz矩阵，由多个单位矩阵组成，单位矩阵和红色箭头对应的P300开始时刻对应，指示着每个P300分量开始时刻。

因此，对正样本脑电信号按式(8)建模，

$ \mathit{\boldsymbol{X}} = \mathit{\boldsymbol{Ds}} + \mathit{\boldsymbol{N}}. $

(8)

其中：$\boldsymbol{X} \in \mathbf{R}^{N_{t} \times N_{s}} $ 为接收到的P300信号，由若干P300分量组成。$\boldsymbol{s} \in \mathbf{R}^{N_{e} \times N_{s}} $ 为信号主分量，$\boldsymbol{D} \in \mathbf{R}^{N_{t} \times N_{e}} $为指示矩阵。N_s是通道数，N_t是时域上的采点数，N_e是P300信号主分量的时域采样点数。N为P300信号中的Gauss白噪声，维数与X相同。可以求解出s的估计值为

$ \mathit{\boldsymbol{\hat s}} = {({\mathit{\boldsymbol{D}}^{\rm{T}}}\mathit{\boldsymbol{D}})^{ - 1}}{\mathit{\boldsymbol{D}}^{\rm{T}}}\mathit{\boldsymbol{X}}. $

(9)

一种特殊情况是，如果正样本信号中多个P300峰不出现混叠，则D中各列正交，$\boldsymbol{D}^{\mathrm{T}} \boldsymbol{D}=\boldsymbol{I}, \hat{\boldsymbol{s}}=\boldsymbol{D}^{\mathrm{T}} \boldsymbol{X} $ 。

2.3 度量学习模型估计最优协方差矩阵

式(7)中的噪声协方差矩阵C是一个关键参数，本文通过度量学习方法估计其最优参数值。由于C是噪声的协方差矩阵，因此必然是一个实对称正定矩阵，且C^-1也是实对称正定矩阵，可以将其作矩阵分解，即存在实矩阵W，使得

$ {{\mathit{\boldsymbol{C}}^{ - 1}} = {\mathit{\boldsymbol{W}}^{\rm{T}}}\mathit{\boldsymbol{W}}.} $

(10)

这样式(7)中线性判定算子可以改写成

$ {T(\mathit{\boldsymbol{x}}) = {\mathit{\boldsymbol{x}}^{\rm{T}}}{\mathit{\boldsymbol{W}}^{\rm{T}}}\mathit{\boldsymbol{Ws}}.} $

(11)

对式(11)增加约束，要求对于正样本，该统计量值较小；对于负样本，统计量值较大。

$ \left\{ {\begin{array}{*{20}{l}} {T(\mathit{\boldsymbol{x}}) < b - 1,\quad {y_i} = 1;}\\ {T(\mathit{\boldsymbol{x}}) > b + 1,\quad {y_i} = - 1.} \end{array}} \right. $

(12)

其中：y_i是样本的标签，b是阈值，i=1, 2, …, N_tr，N_tr是训练样本个数。式(12)可以统一为

$ {y_i}(b - T(\mathit{\boldsymbol{x}})) > 1. $

(13)

定义铰链损失函数为

$ {\rm{loss}} = {\rm{max}}[1 - {y_i}(b - {\mathit{\boldsymbol{x}}^{\rm{T}}}{\mathit{\boldsymbol{W}}^{\rm{T}}}\mathit{\boldsymbol{Ws}}),0]. $

(14)

即如果式(13)成立时，不计入损失项; 如果其不成立时，将误差值计入损失函数。该优化问题的目标函数可以表示为

$ \mathop {{\rm{argmin}}}\limits_\mathit{\boldsymbol{W}} \sum\limits_i {{\rm{max}}} [1 - {y_i}(b - {\mathit{\boldsymbol{x}}^{\rm{T}}}{\mathit{\boldsymbol{W}}^{\rm{T}}}\mathit{\boldsymbol{Ws}}),0]. $

(15)

这个优化过程可以用梯度下降法迭代求解，

$ \boldsymbol{W}_{t+1}=\left\{\begin{array}{ll} \boldsymbol{W}_{t}, & y_{i}\left(b-\boldsymbol{x}^{\mathrm{T}} \boldsymbol{W}^{\mathrm{T}} \boldsymbol{W} \boldsymbol{s}\right)>1 ; \\ \boldsymbol{W}_{t}-\gamma y_{i}, \boldsymbol{W}_{t} \boldsymbol{\varGamma}_{i}, & \text { 其他. } \end{array}\right. $

(16)

其中：Γ_i=xs^T，γ是迭代步长(学习率)。对于满足约束的样本，不更新权重值，否则按照梯度下降的方式来更新。将式(11)变形为

$ T(\mathit{\boldsymbol{x}}) = {\mathit{\boldsymbol{x}}^{\rm{T}}}{\mathit{\boldsymbol{W}}^{\rm{T}}}\mathit{\boldsymbol{Ws}} = {(\mathit{\boldsymbol{Wx}})^{\rm{T}}}(\mathit{\boldsymbol{Ws}}). $

(17)

可以看出，式(17)具有明显的物理意义。T(x)计算的是测试信号x和模板信号s的相关系数，即2个向量间的余弦距离。余弦距离不考虑向量本身的长度，只考虑两者之间的夹角，非常适合脑电信号这类噪声大、时变性强的生物电信号。优化求得的W矩阵相当于对接收到的信号x和模板信号s进行了新的空间投影，在新的投影空间中正负样本的余弦距离进一步拉大，从而进一步增强了分类器的鉴别力。

同时，由于 $ \boldsymbol{x} \in \mathbf{R}^{d \times 1}, \boldsymbol{W} \in \mathbf{R}^{p \times d}$ ，其中d是原始信号的维度，p是新的投影空间的维度，而p远小于d，因此式(17)有效地降低了信号维度。

3 实验结果 3.1 实验设置

实验数据集使用节1.2中介绍的BCI国际竞赛脑电数据集Ⅱ。脑电信号预处理过程包括时域加窗和带通滤波。由于P300信号出现在200~500 ms，因此本文选取刺激开始时刻到之后667 ms之间的时间窗内的信号来作处理和分析。对于原始的脑电数据先进行带通滤波，滤波器选择0.1~10 Hz的Chebyshev I型带通滤波器，滤除低频的直流漂移噪声干扰和高频的工频噪声干扰。对于滤波后的脑电信号再作特征提取和降维。由于带通滤波器的上限截止频率是10 Hz，根据Nyquist采样定理，可以选择20 Hz的频率进行下采样。经过采样，单通道数据的维度减少为20 Hz×0.667 s=14。

本文不同分类器精度的比较实验就是在按照上述方法滤波和降采样处理后的数据集上进行测试的。统计的指标是打字的准确率。由于2名被试都采用100个字符进行测试，因此当其有M个字符打字结果判定正确时，准确率为(M/100)×100%。

3.2 实验结果对比

将本文提出的算法和目前国际领先的算法在BCI国际竞赛上的结果进行比较，如表 1所示。可以看出，随着重复实验轮数的增加，2个被试的识别准确率都在逐渐增高。脑电信号由于背景噪声干扰强，因此需要综合考虑多次重复实验的结果来做出分类决策。同时，由于本文在建立脑电信号分类模型时充分考虑了脑电信号自身噪声大、数据维数高的特点，在估计参数时结合了度量学习方法基于数据驱动求解最优参数，在不同被试的大部分实验轮数上，本文算法取得了最优的准确率。

3.3 算法耗时

本文提出的算法另一个优点在于计算复杂度低。表 2给出了在训练时间和测试时间2个指标上各方法的比较。测试环境为Intel Xeon服务器，32 GB内存，2.40 GHz CPU 16核。对于文[11-12]的算法，由于要训练多组子分类器，因此需要耗费较长的训练和测试时间。对于卷积神经网络算法^[15]，为了拟合复杂的网络，将每个重复实验都作为一个独立的数据而扩充了数据量。模型训练一次需要1 168 s，这在实际情况下难以接受。文[17]算法的求解过程比较复杂，模型训练时间也较长。本文提出的基于匹配滤波器的方法，训练过程根据训练样本估计协方差矩阵C和P300模板信号s，均可以在线性时间复杂度下完成计算，测试时只需要一次线性矩阵乘法，提高了计算效率，在准确率上也高于其他算法。从信息传递速率(information translate rate，ITR)的角度来看，本文方法明显优于其他算法。

4 结论

本文提出了基于通用匹配滤波器的脑电信号分类模型，该模型对两类脑电信号建立了一个生成式模型，并推导出一个简洁的线性判定算子。然后，基于信号建模和最大似然估计求解P300正样本信号的主分量；在估计噪声协方差矩阵时结合了度量学习方法，基于数据驱动迭代求解最优参数。

在公开数据集上不同被试的实验结果表明，本文算法的精度优于其他算法。同时，本文推导出的线性判定算子大大加快了训练和检测速度，同时满足了实际应用系统中高精度和实时性的要求。