基于偏最小二乘法的事件相关电位单次提取研究

严瀚莹 吴 帆 姜忠义 邹 凌

(常州大学信息科学与工程学院 江苏 常州 213164) (常州市生物医学信息技术重点实验室 江苏 常州 213164)

0 引 言

脑电(electroencephalogram,EEG)因其非侵入式、分辨率高、成本低等优点，在临床以及科学研究领域发挥着重要作用。事件相关电位(Event-related potentials,ERP)为EEG的一种衍生物，是由外部施加的特定刺激引起的电位变化[1-2]。目前研究最为广泛的成分是P300,它通常是由受试者在oddball范式中辨认靶刺激时所诱发，潜伏期为300 ms左右的正向电位，主要分布在头顶皮层处[3]。自P300成分被发现以来，它已经被广泛应用于认知心理学、神经科学等领域。

EEG信号在记录时，人脑头皮上任意一点所收集的都是该大脑区域各种活动的非线性混合信号。而诱发电位信噪比很低，很易被强背景噪声所淹没，因此它的提取是ERP研究中的主要挑战之一。目前诱发电位提取技术中，独立成分分析(independent component analysis,ICA)应用最为广泛[4]。但对于成分的选取，ICA需依靠主观经验，而由于利用多导联信号提取诱发电位增加了计算复杂度，先前利用ICA提取诱发电位的研究中只研究了少量导联。例如，Lee等[5]只选取了单个导联来提取P300成分。除ICA之外，其他使用较广泛的方法还有：小波分解、经验模式分解(empirical mode decomposition,EMD)等[6-7]。虽然小波分解在阿尔兹海默症患者等脑电特征提取中效果显著，EMD也在ERP的提取中得到了应用，但这两种方法也存在着一定的局限性。小波分解可能会丢失部分有效成分，选择合理的阈值较为困难，且实时性差[8]。而EMD方法操作性强，依赖于研究者的经验。除上述传统方法，约束ICA、稀疏贝叶斯时空模型、残差迭代分解等方法也被应用于诱发电位的成分提取[9-11]。约束ICA,相比于ICA具有更好的鲁棒性，但也无法解决多导计算量大的问题。贝叶斯时空模型解决了多导信号特征提取困难的问题，但是这种基于多尺度建模的问题，需要先验且算法复杂度高。残差迭代分解法解决了ERP延时性问题，但ERP提取的幅值误差等同样需要考虑。近年来，极限学习机(extreme learning machine,ELM)作为神经网络新算法也应用于脑机接口的特征提取及分类[12]。偏最小二乘(partial least squares,PLS)算法可以对复杂混合体系进行快速分析，最初应用于经济学和化学计量学。PLS算法模型相对简单，对处理变量多样本少的数据具有明显的优势[13]，适用于单次多导脑电信号的ERP提取。2016年Cheung等[14]已经将该方法引入神经影像学，用于脑磁图数据的统计分析，但目前将PLS应用于诱发电位成分提取的研究还比较匮乏。

本文首先构造不同信噪比的模拟脑电数据，采用PLS算法提取P300、N100成分，验证了PLS方法的鲁棒性。基于oddball实验范式采集真实脑电数据，使用小波分解、EMD、ELM及PLS算法提取P300成分，并通过相关系数、峰值误差和峰值潜伏期偏移量指标进一步验证PLS算法的有效性。

1 实验数据及方法

1.1 实验数据

1.1.1模拟脑电数据

由刺激诱发的EEG信号中包含了诱发电位、自发脑电信号以及采集信号时所带的噪声信号。构造模拟脑电数据的过程主要分为三个步骤：

首先利用式(1)对EP信号进行仿真，其中采样点2 016,采样频率为1 008 Hz,采样时间为2 s,波形如图1(a)所示，横坐标为时间(单位毫秒)，纵坐标为幅值，范围为-8～6 μV。

f(t)=-10e(-2t)sin(4πt)

(1)

图1 模拟脑电信号

其次，采用AR模型仿真出自发脑电信号v(n),如式(2)所示，其中u(n)为零均值高斯噪声。之后对真实脑电中的噪声信号进行仿真，这里使用白噪声，采样点为2 016，采样频率为1 008 Hz,采样时间为2 s。将两个信号相加作为诱发电位的背景噪声，如图1(b)所示。

v(n)=u(n)+a1×u(n-1)+a2×u(n-2)+

a3×u(n-3)+a4×u(n-4)

(2)

最后，将仿真自发脑电信号、白噪声和仿真诱发电位，以-9 dB的低信噪比混合作为实验数据，如图1(c)所示，幅值约为-50～50 μV,图中基本看不出N100和P300的波形特征。共仿真20道模拟脑电数据，每道时长为2 s。

1.1.2真实脑电数据

(1) 数据采集。为了验证偏最小二乘方法的可行性，共采集了8名健康被试者(18～22岁)的头皮脑电数据。实验所使用的信号采集设备为美国EGI公司64导脑电帽及放大器，参考电极取Cz点，采样频率为1 000 Hz。本次实验之前所有被试均签署了知情同意书。

(2) 实验范式。实验采用oddball范式，该范式是ERP实验常用范式之一。通常oddball范式在试验中会随机呈现作用于同一感官的两种刺激，两种刺激出现的概率不同。将概率较大的刺激作为标准刺激，概率小的刺激作为靶刺激。本次实验以被试自身证件照为靶刺激(出现概率为20%)、以陌生人证件照为标准刺激(出现概率为80%)，要求被试忽略陌生人证件照，对自身证件照出现时作出相应的按键反应来进行实验。其中靶刺激120次，标准刺激为480次。在每个被试实验前做一次预实验，剔除预选照片中被试较为熟悉的照片。

(3) 数据预处理。为了衡量诱发电位提取的效果，将相应导联上处理后叠加平均所得的诱发电位信号作为参考信号。脑电信号的预处理步骤为：① 降采样，将采样率降为500 Hz;② 滤波，范围为0.1～40Hz;③ 去除眼电伪迹，使用EEGLAB工具箱中的ICA工具;④ 对信号进行分段，刺激前200 ms至刺激后1 000 ms;⑤ 对每种刺激下的信号进行叠加平均，最终得到相对纯净的诱发电位波形。

1.2 方法设计

提取诱发电位的过程实际上就是从记录的脑电信号中去除自发脑电和其他噪声信号。小波分解与经验模态分解都是信号去噪的常用方法，在诱发电位的提取中也应用广泛。极限学习机是神经网络新方法。将偏最小二乘方法与以上三种方法的效果进行对比可验证其有效性。总体流程如图2所示。

图2 P300提取流程图

1.2.1小波变换

离散小波变换(Discrete Wavelet Transformation,DWT)方法不论在频域还是时域上都能较好地反映信号特征。小波分解所采用的窗函数是可变的，增加了时间平移和尺度伸缩两个参数。因此，使用小波分解的方法对信号的高频信息及低频信息都能进行分析。对于分析类似于脑电信号的非平稳信号，小波分解是一种较为合适的方法。

利用小波分解进行成分提取共包含三个方面：基本小波的选择、阈值函数的选取、阈值的选择。在提取成分时小波基选取CoifN小波函数，该小波基相比dbN更能保证信号在分解、重构时不会产生失真和相位延迟。分解层数为7层，根据不同被试脑电信号特性选取分解层进行阈值函数处理，其他分解层置零，最后对信号进行重构。

1.2.2经验模态分解

EMD是根据信号自身的特点，将该信号分解成若干个本征模态函数IMF之和，是一种自适应性的信号处理方法。每一个本征函数必须满足以下两个条件：在整段的数据上，极值点的个数和过零点的个数必须相等或最多相差一个；在数据的任意时间段，上、下包络线的平均值必须是零。

每一个IMF计算步骤如下：

首先，找到原信号x(t)的所有极值点，通过三次样条插值将所有极值点连接起来，拟合出上、下包络线e+(t)和e-(t)。将上、下包络线的均值作为x(t)的均值包络m1(t):

(3)

(4)

将c1(t)从原信号x(t)中分离出来，得到第二个IMF分量c2(t),重复之前的步骤n次后，得到第n个IMF分量cn(t)。当其余量rn(t)小于预设值，或当rn(t)为单调函数或常量时，循环终止。最后得到：

(5)

式中：rn(t)为残余量，代表信号的平均趋势或均值。

1.2.3极限学习机

极限学习机(ELM)是一种针对单隐层前馈神经网络(Single-hidden Layer Feedforward Neural Network,SLFN)的新算法。相较于传统SLFN算法，ELM方法具有学习速度快、泛化性能好等优点。对于一单隐层神经网络，输入样本(Xi,ti)，Xi=[xi1,xi2,…,xin]T∈Rn，ti=[ti1,ti2,…,tim]T∈Rm。通过一个有L个隐层节点的单隐层神经网络，表示为：

(6)

式中：g(x)为激活函数，Wi为输入权重，βi为输出权重，bi是第i个隐层单元的偏置。ELM学习的目标是使得模型输出oj尽可能逼近样本实际输出tj，矩阵表示为：

Hβ=T

(7)

1.2.4偏最小二乘法

实现偏最小二乘回归的基本思想，其实是求如下目标函数最大化的优化问题：

max{Cov(t1，u1)}=max(X0w1,Y0c1)

(8)

该问题求解的算法思想可以理解为解释隐变量与反映隐变量之间的回归关系方程。假定m个自变量{x1,x2,…,xm}和n个因变量{y1,y2,…,yn},构成自变量与因变量的数据表X={x1,x2,…,xm}和Y={y1,y2,…,yn}。在X与Y中分别提取出代表各自特征的成分t1和u1,在提取t1和u1成分时，要求t1和u1应尽可能大地携带它们各自数据表中的变异信息，而且满足t1与u1的相关程度能够达到最大。第一次迭代后，成分t1和u1被提取，X进行对t1的回归，而Y进行对u1进行回归。若回归方程已经达到满意的精度，则成分确定；否则将利用X被t1以及Y被u1解释后的残余信息进行第二轮成分t2和u2提取，继续实施X和Y对t2和u2的回归，该过程一直迭代直到精度满足要求为止。若最终对X共提取了k个成分t1,t2,…,tk,再通过Y对t1,t2,…,tk的回归，最后都可转化为Y对自变量x1,x2,…,xm的回归方程，完成了偏最小二乘的回归建模。

运用到脑电信号中，自变量即为含噪的脑电信号，这里我们用的是64导脑电信号，则X={x1,x2,…,xn}中，n为64，而经过预处理及多试次叠加平均所得的脑电信号便被近似为纯净ERP信号，表示为因变量Y={y1,y2,…,yn},偏最小二乘方法集成了主成分分析、多元线性回归和相关性分析的优点，它主要是通过最小化误差的平方和找到一组数据的最佳函数匹配。用最简单的方法求得一些绝对不可知的真值，而令误差平方之和为最小。脑电信号的单次提取，以减少实验获得质量较好的ERP成分为目的。算法执行按如下操作：

第一次建模，自变量输入为单次trail-1的原始信号，设为X1,其经过预处理后则为Y1,经过PLS建模后得到重构系数矩阵sol1。

第二次建模，自变量输入为单次trail-2的原始信号，设为X2,预处理后的X1与X2迭加平均作为Y2,建模得到sol2,用sol1模型与X2得到的预测信号为y2,ER2为y2与Y2之间的误差。

第n次建模，自变量输入为单次trail-n的原始信号Xn,前n次预处理后的脑电信号的均值作为Yn,建模得到soln,若soln-1模型与Xn得到的预测信号yn与Yn之间的误差小于限定误差ER,则soln将作为最终实验模型，此时迭代次数为n。

1.3 评价指标

为了检测PLS方法应用于诱发电位提取的有效性，使用了如下几个评价指标：相关性系数、峰值误差、峰值潜伏期偏移量。

1.3.1相关系数

相关系数主要考虑两个信号在一段时间内的同步程度。用在信号重构中，可用于检测信号重构的质量。

(9)

1.3.2峰值误差

对于ERP的提取中，我们所关注的成分的直接表现就是波形的幅值，峰值误差表示了模型提取的成分波与理想成分波在峰值上的误差的绝对值。表示为：

Errpeak=|xm-yn|m∈(k,l),n∈(k,l)

(10)

在特定区间(k,l)内，信号X在m时刻取得最高值(或最低值)，而信号Y在n时刻取得最高值(或最低值)，则其峰值差的绝对值为峰值误差。

1.3.3峰值潜伏期偏移量

除了信号的幅值以外，峰值相位的准确度也是检验成分提取质量的标准之一。峰值潜伏期偏移量则是在峰值确定下时刻m与时刻n点所对应的实际时间：

Err=|m-n|×t

(11)

式中：t为信号样本点间的时间最小间隔。

2 结果分析

2.1 模拟数据

对于模拟的含噪脑电信号，首先使用PLS方法迭代共12次得到最佳提取模型，之后利用该模型提取模拟含噪脑电信号中的N100、P300成分。观察提取信号的时域波形可以直观反映提取效果，图3为信噪比为-7 dB模拟信号的理想ERP与PLS方法提取波形的对比图。实线波形为理想ERP信号，虚线波形为PLS方法提取波形。从图中可以直观看出，PLS方法提取的诱发电位与理想ERP波形相似性很高，N100与P300成分也都很明显。

图3 PLS提取波形与理想ERP对比图

为了衡量PLS方法的稳定性，对不同信噪比的模拟含噪脑电信号用PLS方法进行成分提取。将提取后的成分潜伏期与理想ERP成分潜伏期作差，该差值即为PLS方法成分提取的潜伏期偏移量。对信噪比为0 dB至-9 dB的模拟脑电信号进行成分提取，N100、P300的潜伏期偏移量如图4所示。可以看出N100、P300潜伏期偏移量都在7 ms之内，误差较小。其中N100潜伏期偏移量较P300偏大，P300潜伏期最大偏移量为4 ms而N100与真实位置偏差最大达到6.2 ms。但最后噪声再增加而N100、P300偏移量将保持稳定。

图4 PLS模型在不同信噪比下模拟数据的潜伏期偏移量

2.2 真实脑电数据

将PLS应用于真实的脑电数据，将提取结果与小波变换、经验模式分解以及极限学习机方法的结果进行对比。图5为同一个被试Pz电极处利用四种不同方法提取的时域波形，(a)为该被试单试次原始脑电波形，(b)为四种方法提取波形与叠加平均(Ensemble Averaging,EA)后的ERP波形对比。(b)中PLS方法在P300处与EA波形最为接近，ELM结果次之，EMD方法峰值与EA相差过大而DWT方法峰值较EA有明显偏移。PLS方法较其他三种方法时域波形结果较好。

图5 不同方法诱发电位提取结果

为了定量衡量不同方法时域波形的提取效果，对P300成分较为明显的P3、P4、C3、C4导联计算提取信号与EA的相关系数。从图6中可明显看出EMD方法相关系数最低，PLS以及ELM方法在四个导联上的相关系数都在0.7之上，小波分解方法除C3导联外相关系数都比PLS方法低。PLS和ELM较DWT和EMD方法与EA有更强相关性，且它们在多个电极的ERP提取上都表现出良好性能。

图6 不同方法提取信号相关系数

对四种方法在8个被试的P300提取中的偏移误差进行计算，如表1所示。PLS方法的潜伏期偏移量与幅值偏移量较其他三种方法平均值都更小，方差统计也显示出其更稳定的性能。

表1 四种方法在P300提取中的偏移误差比较

3 结 语

在仿真数据的实验中，偏最小二乘模型在特征提取上具有很好的抗干扰性能，能成功提取N100和P300成分。而在真实脑电ERP的提取中，虽然小波、经验模态分解以及极限学习机在使用大量试次(120次)叠加平均建模的基础上提取到了准确的P300成分，但偏最小二乘模型在少试次(最少9次)的建模中得到了更加稳定的模型。在后续的验证试验中，PLS与ELM提取的P300成分，在相关性指标上均具有较好表现。而PLS在信号峰值误差、潜伏期误差等指标上都优于其他三个模型。PLS模型实现了少试次精确建模，解决了多导ERP成分的单试次提取困难，为脑电ERP研究以及BCI脑机接口的发展提供了技术支持。