CN108829257A - 一种运动想象脑电信号的特征提取方法 - Google Patents

Abstract

本发明公开了一种运动想象脑电信号的特征提取方法。根据平均功率谱确定对运动想象脑电信号进行特征提取的最优时间段，然后对这一时间段内的运动想象脑电信号进行四层双树复小波分解，利用每个子带的复小波系数进行信号重构，通过计算重构后信号最有时间段内数据的平均能量特征作为运动想象脑电信号的时频特征；提出IL‑MVU算法对重构后的脑电信号的最有时间段内数据进行降维，将降维后得到的低维向量作为运动想象脑电信号的非线性特征，最后通过对最优时间段内运动想象脑电信号的时频特征与非线性特征进行标准化处理与特征融合，得到运动想象脑电信号的特征向量。本发明极大地降低了算法的时间消耗，有助于提高MI‑EEG信号的分类准确率。

Description

一种运动想象脑电信号的特征提取方法

技术领域

本发明属于脑电信号处理方法，具体应用于脑-机接口(Brain-ComputerInterface,BCI)系统中对运动想象脑电信号的特征提取，对地标点最大方差展开(Landmark Maximum Variance Unfolding,L-MVU)进行改进，提出了一种增量式地标点最大展开方法(Incremental,L-MVU，IL-MVU)，并与双树复小波变换(Dual Tree ComplexWavelet Transform,DTCWT)相结合，实现对运动想象脑电信号的特征提取与融合。

背景技术

运动想象脑电信号(Motor Imagery Electroencephalography,MI-EEG)中包含大量的生理信息，并且与意识状态密切相关。因此，对MI-EEG的识别成为脑-机接口系统中的关键，而从运动想象脑电信号中获得特征的优劣将直接影响对其的识别精度。

运动想象脑电信号是一种非线性、非平稳及时频特性比较明显的信号，如何充分利用这些特征是精确提取特征的关键。双树复小波变换由两个并行的离散小波变换(Discrete Wavelet Transform,DWT)构成，这两个小波变换分解产生的小波系数作为复小波系数的实部与虚部，构成实部树与虚部树。双树复小波变换实现了实部树与虚部树信息的互补，既能够实现从多个尺度上对信号进行分解，又具有完全重构性，并且可以有效消除频率混叠。但是仅仅提取运动想象脑电信号的时频特征并不能得到利于分类的完备的特征，流形学习(Manifold Learning)的出现能够从非线性角度提取特征，同时解决了高维度特征不利于分类的缺点。流形学习能够从高维数据中恢复其内在的低维流形结构，在维度约简的同时最大化保留原始数据信息量完整，并且低维数据更利于可视化。L-MVU作为一种典型的流形学习方法，相比于其他流形学习方法比如等距映射(Isometric Mapping,Isomap)和局部线性嵌入算法(Locally Linear Embedding,LLE)，能够在降维的同时保留距离与角度的信息，从而提取到更有利于分类的特征，已经在MI-EEG信号提取上得到了一些应用。

然而，涉及到实际应用，而对于上述提到的L-MVU，也存在一些缺点。第一，L-MVU仅仅能对给定数据集进行维数约减，无法实时处理新来的数据，使得L-MVU对于样本外数据的泛化能力不足，不利于BCI系统的在线实现；第二，L-MVU提取的非线性特征无法体现MI-EEG的时频特性，对于时频特征提取的缺失会对最终的模式分类产生一定的影响。

发明内容

针对现有流形学习方法在运动想象脑电信号特征提取领域中应用存在的不足，本发明提出了一种基于双树复小波变换与增量式地标点最大方差展开的运动想象脑电信号特征提取的方法。利用此方法首先在重构特定频段脑电信号上避免了引入不期望的频段；其次，保证了对运动想象脑电信号特征提取的紧致性与完备性，提高了分类准确率；最后，提出L-MVU的增量式版本IL-MVU应用于MI-EEG的特征提取，使得对于样本外数据的泛化能力大大增强，并且在测试时间上也相应的缩短，更有利于本方法在实际应用中的在线实现。

本方法的技术路线如下：

首先，根据平均能量谱确定对运动想象脑电信号进行特征提取的最优时间段，然后对原始运动想象脑电信号进行四层双树复小波分解，利用每个子带的复小波系数进行信号重构，通过计算重构后信号的平均能量，将其作为运动想象脑电信号的时频特征；同时，利用IL-MVU算法对重构后的脑电信号进行降维，将降维后得到的低维向量作为运动想象脑电信号的非线性特征，接着通过对时频特征与非线性特征进行标准化处理，得到运动想象脑电信号的特征向量，最后将特征向量输入LDA分类器进行分类，并通过遍历的方式得到IL-MVU算法的最优参数。

脑电信号特征提取的生理学依据是：当大脑皮质某区域被激活时，该区域的代谢和血流增加，导致脑电波α节律(8～13Hz成分)和β节律(14～30Hz成分)幅度的降低，称为事件相关去同步(Event-Related Desynchronization,ERD)；在大脑静息或惰性状态下，α节律和β节律表现出幅度明显增高，被称为事件相关同步化(Event-RelatedSynchronization,ERS)。因此，人在想象单侧手运动时，其对侧相应初级感觉运动皮层的脑电α节律和β节律幅值会降低；而同侧对应节律的EEG幅值会升高。这为后续EEG信号提取时频特征特定子频带的确定，提供了理论基础。脑电信号的这种事件相关去同步与事件相关同步化现象成为分析和判断左右手运动想象脑电信号的最根本依据。

基于上述分析，本方法具体按照如下步骤实现：

步骤1，信号预处理。

根据MI-EEG信号产生的特点，MI-EEG的α节律和β节律中包含了特征较为明显的信息，并且C3和C4两通道的MI-EEG信号的事件相关同步/事件相关去同步现象最为明显，因此首先根据所有实验的C3和C4两导信号的平均能量谱分析确定脑电特征的有效时间范围，ch导联平均能量由下式计算：

其中，N代表运动想象任务次数，da^ch(u，v)代表MI-EEG信号第u次运动想像任务ch导联的第v个采样点的信号幅值，ch表示导联名称，ch＝{c3，c4}。依此得到想象左右手运动的平均能量谱，根据平均能量谱选取ERS/ERD现象最为明显的时间段OT[min，max]；

步骤2，对原始脑电信号进行L层双树复小波分解并对各子带的小波系数进行重构。

双树复小波变换由两次离散小波变换构成，两次离散小波变换分别构成实部树与虚部树，信号S(t)经双树复小波分解后得到两组小波系数与尺度系数，分别记为和j＝1,2,…,L，由此得到双树复小波分解后各层复小波系数与复尺度系数表示如下：

其中，i为虚数单位。

假设由重构得到j尺度下的细节分量Dt_j，由重构得到近似分量A_L，则对于信号S^ch(t)的L层双树复小波分解可以表示为：

从而，将信号S^ch(t)划分为多个子带信号。设信号S^ch(t)的采样频率为f_s，则 各分量分别所对应的频带范围依次为：[0，f_s/2^L+1]，[f_s/2^L+1，f_s/2^L]，[f_s/2^L，f_s/2^L-1]，…，[f_s/2²，f_s/2]。

步骤3，时频特征提取。

对步骤2中重构好的子带信号A_L，Dt_L，Dt_L-1…Dt₁进行截取，将其最优时间段OT内的数据分别记为设置长为2s的滑动时间窗，计算的能量如下：

针对C3和C4导联信号，选取频带范围与α节律和β节律相近的两个分量，分别记为并对其进行能量归一化：

每次滑动一个采样点，计算滑动时间窗内两个子带信号的归一化能量获得相应能量序列和其长度为(max-min-2fs+2)。进而，计算每个能量序列的均值

最后，由下式计算C3与C4导联在α节律和β节律的平均能量差AE_α,AE_β：

因此，对于每次运动想象任务，可以得到一个二维时频能量特征，记为F₁，表示如下：

F₁＝[AE_α AE_β]^T∈R^2×1 (9)

步骤4，非线性特征提取。

计算的平均功率谱，ry＝{α，β}，选取ERS/ERD现象最为明显的节律，得到(ry＝α或ry＝β)，以作为IL-MVU的输入。

步骤4.1，假设表示第u次运动想像任务的第v个采样点，为了得到更明显的初始特征，定义dS(u，v)如下式：

从而，获得第u次运动想像任务的初始高维特征D_u∈R^{(max-min+1)×1}如下：

D_u＝[dS(u，min)，dS(u，min+1)，...，dS(u，max)]^T，u＝1，2，...，N (11)

进而，构成初始高维特征集D∈R^{(max-min+1)×N}，并将其分为高维训练特征集为Tr∈和高维测试特征集为两部分。

步骤4.2，针对高维训练特征集Tr设定L-MVU算法初始参数：地标点个数为m，低维特征维数为d，计算权值矩阵的近邻数为r，求解SDP问题的近邻数为k。将高维训练特征集Tr输入L-MVU算法，得到Tr的低维表示如下式：

步骤4.3，提出IL-MVU算法，并用于求取高维测试集Te的低维特征集Y^te。

针对高维测试特征集Te的一个新数据Te_u(u＝1，2，...，N-n₁)，首先在高维训练数据集Tr中求得其p个近邻点，记为p_uv(v＝1，2，...，p)，得到Te_u的近邻点集合Ne_u为：

Ne_u＝{p_u1，p_u2，...，p_up} (13)

通过最小化式(14)中的ε，得到Te_u的重建权值向量W_u＝[W_u1，W₂，...，W_ip]。

式(14)满足约束条件：

在Y^tr中找到与近邻点p_uv，相对应的低维表示,记为q_uv,按照高维空间中距离很近的点经过流行学习方法降维后在低维空间中依旧距离很近的原则，设重建权值向量W_i不变，则高维测试集特征Te_u的低维非线性特征计算如下：

从而，获得低维测试特征集为：

步骤4.4，对于任意一次运动想象，可以获得一个d维非线性特征F₂∈R^d×1。因此，可得到相应的低维特征训练集及特征测试集

步骤5，特征融合。

将步骤3得到的特征F₁与步骤4得到的特征F₂进行特征融合，得到本方法提取的特征向量F∈R^(2+d)×1，记作：

其中，将F₁放大100倍是为了与F₂达到同一数量级，以平衡两类特征的作用，有利于提高分类精度。

与现有技术相比，本发明具有如下优点：

(1)本发明针对传统的流形学习方法对于训练集外数据泛化能力不强的缺点，在L-MVU方法的基础上提出了其增量式版本，经过实验验证，IL-MVU在分类精度上较普通版本有一定提升，并且在测试集的特征提取速度上取得了量级的加快。

(2)本发明在时频特征提取方面采用双树复小波变换方法分解与重构各子带信号，与传统的离散小波变换相比，双树复小波变换在能够更精确的提取相关频带信息并重构，不易引入不期望的频率信息，使特征提取更加精确，并且在提取能量特征方面采用了归一化平均能量序列再次求平均的方法，避免同一想象运动在不同的信号采集周期内所获取的脑电信号强度不一致造成的特征不稳定，使得到的能量特征更易于区分；

(3)本发明相比传统方法提取单一特征的做法，从时频特征与非线性特征两方面提取脑电信号的特征，取得的特征包含了更多表征脑电信号特点的信息，在使用LDA进行模式分类的实验中，取得了十折交叉验证后90.71％的分类正确率，相较其他其他传统的流形学习方法与时频方法具有优越性。

附图说明

图1.1为想象左手运动原始信号C3、C4导联平均能量谱。

图1.2为想象右手运动原始信号C3、C4导联平均能量谱。

图2.1为想象左手运动S_α重构信号C3、C4导联平均能量谱。

图2.2为想象右手运动S_α重构信号C3、C4导联平均能量谱。

图2.3为想象左手运动S_β重构信号C3、C4导联平均能量谱。

图2.4为想象右手运动S_β重构信号C3、C4导联平均能量谱。

图3为电极位置分布。

图4为实验时序图。

图5为本方法的实施流程图。

具体实施方式

本发明中具体实验是在Windows 10(64位)系统使用Matlab2017a的仿真环境下进行的。

本发明采用的MI-EEG数据集来源于第三届脑机接口竞赛数据集一，由奥地利格拉茨理工大学BCI研究中心提供。整个实验由280次实验组成，其中140次用于训练，140次用于测试，采用AgCl作为电极，采样频率是128Hz，电极放置如图2所示，电极放置遵循国际标准的10-20导联系统的C3、CZ和C4三个导联通道。每次实验持续9s，具体时序如图3所示。在t＝0～2s时，受试者保持休息状态；在t＝2s时，显示器上持续显示一个十字光标，并给出短暂的提示音，实验开始；在t＝3s时十字光标被一个随机产生的向左或向右的箭头代替，同时要求受试者接下来根据箭头引导想象左右手的运动；在t＝4～7s时，受试者想象左右手运动；在t＝8～9s时，受试者休息，准备下一次实验。

结合具体脑电数据集，如图5所示流程图，本发明算法具体实施步骤如下所述：

步骤1，信号预处理。根据式(1)计算计算MI-EEG数据集的平均能量，结果如图1.1与图1.2所示，由图中分析，C3和C4导的左右手MI-EEG信号的平均功率在相应大脑的运动感觉区会表现出明显的ERD和ERS现象，并且在第450-第1000采样点时段内表现最为突出，为此选择最优时间段的采样点为OT[450,1000]。

步骤2，对原始脑电信号进行四层双树复小波分解，并且对各子带进行重构，得到重构后的原始脑电的信号分量A₄,Dt₄,Dt₃,Dt₂,Dt₁,分别对应的频带范围是0～4Hz，4～8Hz,8～16Hz,16～32Hz,32～64Hz。

步骤3，时频特征提取。

将步骤2中重构好的原始MI-EEG信号分量A₄,Dt₄,Dt₃,Dt₂,Dt₁截取最优时间段OT内的数据分别记为然后选取频带范围与α节律和β节律相近的两个分量的C3和C4导联的信号分别记为

接着对做如下处理：首先选取长为2s的滑动时间窗(256采样点)，窗口以一个采样点作为滑动间隔，每滑动一次，利用式(5)，式(6)计算时间窗内两个子带信号的归一化能量直到最后一个窗口，由此得到四个能量序列 然后根据式(7)计算得到其的平均值最后由式(8)计算得到平均能量差AE_α,AE_β，并由式(5)重组得到时频特征F₁。

步骤4，非线性特征提取。

对于步骤3中定义的分量，分别如步骤1中进行平均能量谱绘制，如图2.1～图2.4所示。然后选取最优时间段内ERS/ERD现象最为明显的分量数据作为非线性特征提取的输入。

步骤4.1，如式(10)，式(11)计算得到的初始高维特征集D，并划分高维训练特征集Tr∈R^551×252高维测试特征集Te∈R^551×28

步骤4.2，对于高维训练特征集Tr,设置L-MVU算法参数d＝3,k＝4,r＝40,m＝5,得到高维训练特征集的低维表示Y^train。

步骤4.3，对于高维测试特征集中的数据Te_u(i＝1，2，...，28)，从高维训练特征集Tr中求得如式(13)所示的集合Ne_u，并由式(14)计算出重建权值W_u。按照高维空间中距离很近的点在经过流行学习方法降维后在低维空间中依旧距离很近的原则，设重建权值向量W_u不变，由式(15)求得测试集数据的低维非线性特征

步骤4.4，由步骤4.2所得的低维训练非线性特征与步骤4.3所得低维测试非线性特征记为F₂

步骤5，将步骤3所得F₁与步骤4所得F₂如式(18)进行特征融合，得到最终的脑电特征F。

为对本发明提出的基于双树复小波变换与增量式地标点最大方差展开的运动想象特征提取方法的有效性进行验证，利用LDA分类器对本发明提取的特征进行了模式分类与增量式L-MVU参数寻优，参数寻优采用遍历的方法，最终得到的最优参数为：d＝5,k＝26,r＝48,m＝14,最高十折分类正确率为90.71％，具体结果如表1所示。

表1与相关研究结果的对比

同时，在测试集低维非线性特征的产生时间上，本方法较传统流形学习方法也表现出了显著的优越性，生成低维非线性特征的时间如表二所示

表2与传统流行学习方法的时间消耗对比

Claims (2)

1.一种运动想象脑电信号的特征提取方法，其特征在于：

首先，根据平均能量谱确定对运动想象脑电信号进行特征提取的最优时间段，然后对原始运动想象脑电信号进行四层双树复小波分解，利用每个子带的复小波系数进行信号重构，通过计算重构后信号的平均能量，将其作为运动想象脑电信号的时频特征；同时，利用IL-MVU算法对重构后的脑电信号进行降维，将降维后得到的低维向量作为运动想象脑电信号的非线性特征，接着通过对时频特征与非线性特征进行标准化处理，得到运动想象脑电信号的特征向量，最后将特征向量输入LDA分类器进行分类，并通过遍历的方式得到IL-MVU算法的最优参数；

脑电信号特征提取的生理学依据是：当大脑皮质某区域被激活时，该区域的代谢和血流增加，导致脑电波α节律即8～13Hz成分和β节律即14～30Hz成分幅度的降低，称为事件相关去同步即ERD；在大脑静息或惰性状态下，α节律和β节律表现出幅度明显增高，被称为事件相关同步化即ERS；因此，人在想象单侧手运动时，其对侧相应初级感觉运动皮层的脑电α节律和β节律幅值会降低；而同侧对应节律的EEG幅值会升高；这为后续EEG信号提取时频特征特定子频带的确定，提供了理论基础；脑电信号的这种事件相关去同步与事件相关同步化现象成为分析和判断左右手运动想象脑电信号的最根本依据。

2.根据权利要求1所述的基于DTCWT和Incremental L-MVU的运动想象脑电信号的特征提取方法，其特征在于：

步骤1，信号预处理；

根据MI-EEG信号产生的特点，MI-EEG的α节律和β节律中包含了特征较为明显的信息，并且C3和C4两通道的MI-EEG信号的事件相关同步/事件相关去同步现象最为明显，因此首先根据所有实验的C3和C4两导信号的平均能量谱分析确定脑电特征的有效时间范围，ch导联平均能量由下式计算：

其中，N代表运动想象任务次数，da^ch(u，v)代表MI-EEG信号第u次运动想像任务ch导联的第v个采样点的信号幅值，ch表示导联名称，ch＝{c3c4}；依此得到想象左右手运动的平均能量谱，根据平均能量谱选取ERS/ERD现象最为明显的时间段OT[min，max]；

步骤2，对原始脑电信号进行L层双树复小波分解并对各子带的小波系数进行重构；

双树复小波变换由两次离散小波变换构成，两次离散小波变换分别构成实部树与虚部树，信号S(t)经双树复小波分解后得到两组小波系数与尺度系数，分别记为和j＝1,2,…,L，由此得到双树复小波分解后各层复小波系数与复尺度系数表示如下：

其中，i为虚数单位；

假设由重构得到j尺度下的细节分量Dt_j，由重构得到近似分量A_L，则对于信号S^ch(t)的L层双树复小波分解可以表示为：

从而，将信号S^ch(t)划分为多个子带信号；设信号S^ch(t)的采样频率为f_s，则 各分量分别所对应的频带范围依次为：[0，f_s/2^L+1]，[f_s/2^L+1，f_s/2^L]，[f_s/2^L，f_S/2^L-1]，…，[f_s/2²，L_s/2]；

步骤3，时频特征提取；

对步骤2中重构好的子带信号A_L，Dt_L，Dt_L-1…Dt₁进行截取，将其最优时间段OT内的数据分别记为设置长为2s的滑动时间窗，计算的能量如下：

针对C3和C4导联信号，选取频带范围与α节律和β节律相近的两个分量，分别记为 并对其进行能量归一化：

每次滑动一个采样点，计算滑动时间窗内两个子带信号的归一化能量获得相应能量序列和其长度为(max-min-2fs+2)；进而，计算每个能量序列的均值

最后，由下式计算C3与C4导联在α节律和β节律的平均能量差AE_α,AE_β：

因此，对于每次运动想象任务，可以得到一个二维时频能量特征，记为F₁，表示如下：

F₁＝[AE_α AE_β]^T∈R^2×1 (9)

步骤4,非线性特征提取；

计算的平均功率谱，ry＝{α，β}，选取ERS/ERD现象最为明显的节律，得到以作为IL-MVU的输入；

步骤4.1，假设表示第u次运动想像任务的第v个采样点，为了得到更明显的初始特征，定义dS(u，v)如下式：

从而，获得第u次运动想像任务的初始高维特征D_u∈R^{(max-min+1)×1}如下：

D_u＝[dS(u，min)，dS(u，min+1)，...，dS(u，max)]^T,u＝1,2,…,N (11)

进而，构成初始高维特征集D∈R^{(max-min+1)×N}，并将其分为高维训练特征集为 和高维测试特征集为；

步骤4.2，针对高维训练特征集Tr设定L-MVU算法初始参数：地标点个数为m,低维特征维数为d，计算权值矩阵的近邻数为r，求解SDP问题的近邻数为k；将高维训练特征集Tr输入L-MVU算法，得到Tr的低维表示如下式：

步骤4.3，提出IL-MVU算法，并用于求取高维测试集Te的低维特征集Y^te；

针对高维测试特征集Te的一个新数据Te_u(u＝1，2，...，N-n₁)，首先在高维训练数据集Tr中求得其p个近邻点,记为p_uv(v＝1，2，...，p),得到Te_u的近邻点集合Ne_u为：

Ne_u＝{p_u1，p_u2，...，p_up} (13)

通过最小化式(14)中的ε,得到Te_u的重建权值向量W_u＝[W_u1，W₂，...，W_ip]；

式(14)满足约束条件：

在Y^tr中找到与近邻点相对应的低维表示,记为q_uv,按照高维空间中距离很近的点经过流行学习方法降维后在低维空间中依旧距离很近的原则，设重建权值向量W_i不变，则高维测试集特征Te_u的低维非线性特征计算如下：

从而，获得低维测试特征集为：

步骤4.4，对于任意一次运动想象，可以获得一个d维非线性特征F₂∈R^d×1；因此，可得到相应的低维特征训练集及特征测试集

步骤5，特征融合；

将步骤3得到的特征F₁与步骤4得到的特征F₂进行特征融合，得到本方法提取的特征向量F∈R^(2+d)×1，记作：

其中，将F₁放大100倍是为了与F₂达到同一数量级，以平衡两类特征的作用，有利于提高分类精度。