CN109858537A - 改进的eemd与csp结合的脑电信号特征提取方法 - Google Patents

Abstract

本发明提供一种改进的EEMD与CSP结合的脑电信号特征提取方法，包括：对脑电信号集进行预处理，将其划分为训练集和测试集；对预处理后的脑电信号进行集合经验模式分解，得到若干个固有模态函数；对脑电信号进行FIR滤波，每一个固有模态函数被分解为q个子带信号，各个固有模态函数中频段相同的子带信号构成一个子带信号集，得到q个子带信号集；将q个子带信号集分别经过公共空间模式进行特征提取，得到q个特征向量，将q个特征向量组合为一个特征矩阵；训练支持向量机；将待分类脑电信号的特征矩阵输入支持向量机进行分类。本发明在CSP的基础上结合了EEMD的频域信息，并且经过滤波器组过滤优化，能够更好的用于运动想象特征提取，获得更高的分类精度。

Description

改进的EEMD与CSP结合的脑电信号特征提取方法

技术领域

本涉及智能信息处理领域，特别涉及一种集合经验模式分解算法融合FIR滤波器与公共空间模式算法结合的脑电信号特征提取方法。

背景技术

脑机接口技术为是一种基于脑电信号(EEG)为人类与外部设备之间提供通信与控制的接口，通过脑机接口技术可以达到通过人脑控制外界环境的目的。脑机接口的研究主要是运用一种非侵入方法来反映出脑活动和精神状态，脑机接口为神经肌肉系统功能受损严重的病人提供了一种新的与外界交流的手段，它代替了传统的神经肌肉通道直接通过意念来与电脑程序交流并能够控制外部设备。近年来，由于神经学科，模式识别，信号处理和电子测量技术的推动，脑机接口技术的引起了越来越多的关注。

脑机接口技术主要包含五个步骤：信号采集；信号预处理，特征提取，模式识别与分类以及控制外部装置。如何有效的提取出脑电图的特征并改善分类精度是BCI(BrainComputer Interface，脑机接口)技术研究的关键。从时域和频域单独分析脑电信号是最直接的简单的特征提取方法，常见的时域特征分析方法主要包括：波形特征描述和自回归AR模型；常见的频域分析方法主要分为经典谱估计和现代谱估计。由于时频域算法具有一定的缺陷和局限性，该方法不能同时得到时间-频率良好的分辨率。脑电图信号包含大量的信息，复杂程度较高，仅仅在单一频域中分析并不能较好的分析出信号的特征。因此在几年的研究中，往往将两到三种不同分析域的方法相结合。CSP作为一种特征提取算法被广泛的应用。Z.J.Koles于1991年首次用CSP(common spatial patterns，公共空间模式)方法来提取EEG(Electroencephalogram，脑电波)信号中的非正常成分；文献首次将CSP运用于运动想象脑电图的分类中；但是CSP因为单独使用而缺乏频域信息而具有局限性；2007年Quadrianto Novi等人提出了滤波器组公共空间模式，将信号经过选择好的滤波器组后再经过CSP，其方法融入了频域信息，并获得了较高的分类精度；Yang Banghua等人提出了小波包分解与CSP相结合的方法提取运动想象脑电信号特征；而在国内，李明爱等人采用的改进后的CSP也获得了较好的分类性能。

CSP特征提取算法已经取得了广泛的应用，但是其有着输入导联数多并且缺乏频域等缺陷。本文在CSP的基础上结合了EEMD(Ensemble Empirical Mode Decomposition，集合经验模态分解)的频域信息，并且经过滤波器组过滤优化，能够更好的用于运动想象特征提取，获得更高的分类精度。

发明内容

本发明目的在于针对传统公共空间模式算法的不足，提出了一种公共空间模式算法结合集合经验模式分解的脑电信号特征提取方法。该方法是基于如下内容：1、基于集合经验模式分解的脑电信号处理；2、频域能量分析；3、根据频谱分析筛选固有模态函，并经过FIR滤波器重构新的信号矩阵；4、经过公共空间模式分解，解决公共空间模式多输入、缺频域信息的问题；5、利用支持向量机分类。

本发明解决其技术问题所采取的技术方案是：

一种改进的EEMD与CSP结合的脑电信号特征提取方法，包括如下步骤：

步骤1：对所采集脑电信号进行预处理，所述脑电信号集由左手运动想象信号和右手运动想象信号组成；

步骤2：对预处理后的脑电信号进行集合经验模式分解，得到若干个固有模态函数；

步骤3：对步骤2所得脑电信号进行FIR滤波，其中每一个固有模态函数被分解为 q个子带信号，各个固有模态函数中频段相同的子带信号构成一个子带信号集，得到q 个子带信号集；

步骤4：将q个子带信号集分别经过公共空间模式进行特征提取，得到q个特征向量，将q个特征向量组合为一个特征矩阵；

步骤5：将若干个脑电信号的特征矩阵划分为训练集和测试集，利用所述训练集训练支持向量机，并用所述测试集进行验证；

步骤6：利用步骤1-4得到待分类脑电信号的特征矩阵，输入步骤5所得支持向量机进行特征分类。

优选的，步骤1中预处理包括ICA(Independent Component Analysis，独立成分分析) 去伪迹、8-32Hz带通滤波。

优选的，步骤2中集合经验模式分解的步骤包括：

步骤2.1：原始脑电信号记为：x(t)＝[x₁(t),x₂(t),…,x_n(t)]∈R^N×n，其中n代表通道数目，N代表样本点，x_n(t)代表第n个通道的脑电信号；

步骤2.2：产生白噪声将白噪声增加到信号段x_i(t)，i＝(1,2,…n)，x_i(t)代表第i个通道的脑电信号，得到即i＝1代表集合经验模式分解过程的开始；

步骤2.3：对进行经验模式分解得到k个固有模态函数IMFs，即： x_i,j,1(t),x_i,j,2(t),x_i,j,3(t),…x_i,j,k(t)，其中，x_i,j,l(t)，l∈(1.2.…,k)代表第i个通道脑电信号的第j次经验模式分解得到的第l个IMF，j∈(1,2,…M),M代表经验模式分解进行的总次数；

步骤2.4：将得到的IMFs组成矩阵形式得到：x_i,j(t)＝[x_i,j,1(t),x_i,j2(t),,…x_i,j,k(t)]；

步骤2.5：重复步骤2.2-2.4，直到j＝M，计算所得到的M个矩阵的平均值，即：

为第i个通道经过集合经验模式分解后的IMFs矩阵；x_i,j(t)表示第i个通道第j次经验模式分解后得到的IMFs矩阵；

经集合经验模式分解后，每一个通道的脑电信号x_i(t)，i＝(1,2,…,n)被分解为k个 IMFs，即：

x_i,r(t)∈R^N×1(r＝1,2,…k)

因此有：

x(t)＝[x_1,1(t),x_1,2(t),…,x_1,k(t),…x_n,1(t),x_n,2(t),…x_n,k(t)]。

优选的，步骤3中对脑电信号进行FIR滤波包括如下步骤：

步骤3.1：脑电信号经过集合经验模式分解之后，选取第i个通道的前g(0＜g≤k)个 IMFs，其中每一个IMF被FIR滤波器分解为q个子带信号；

因此有：

x_i,j(t)＝[x_i,j,1(t),x_i,j2(t),…,x_i,j,q(t)]∈R^N×q

故第i个通道的信号x_i(t)表示为：

x_i(t)＝[x_i,1(t),x_i,2(t),…x_i,qg(t)]∈R^N×(qg)；

步骤3.2：将选取的g个IMF对应的相同频段的信号重构成一个新的矩阵，因此信号最终重构为q个矩阵。

优选的，步骤4中将q个子带信号集分别经过公共空间模式进行特征提取包括如下步骤：

步骤4.1：脑电信号的标准化空间协方差矩阵为：

其中X_H和X_F分别代表预处理后的左手运动想象信号和右手运动想象信号，维度为n*N，其中n代表通道数目，N代表每个通道的样本点，X^T代表X的转置，trace(*) 为矩阵*对角元素的总和；

步骤4.2：平均协方差矩阵计算为：

其中Y^k＝{H,F}代表第k次试验的类标签，H和F分别表示左手运动想象和右手运动想象；

步骤4.3：复合协方差矩阵对角化分解为：

其中U₀为矩阵R的特征向量矩阵，∑为对应的特征值的对角化矩阵；

步骤4.4：白化矩阵：

P＝∑^1/2U₀ ^T

将平均协方差矩阵白化进行白化：

步骤4.5：S_H与S_F有着共同的特征向量U，两者所对应特征值构成的矩阵和为单位矩阵I；

S_F＝U∑_FU^T S_H＝U∑_HU^T

∑_H+∑_F＝I

S_H最大特征值的特征向量对应S_F的最小特征值，反之亦然；白化的脑电信号在对应于∑_H与∑_F最大特征值的特征向量上进行变换，得到两类信号矩阵的最优分离方差；

步骤4.6：投影矩阵W表示为：

W＝U^TP

步骤4.7：利用投影矩阵W将原始脑电信号进行转化：

Z₀＝WX

矩阵Z＝[z₁,z₂,……,z_2m]∈R^N×2m由Z₀的前和后m行组成，定义特征向量为：

f＝[f₁,f₂.…,f_2m]^T∈R^2m×1

其中：

优选的，步骤1中所述脑电信号为所述脑电信号为BCICompetition 2008 datasets 2b 数据集中的C3、C4通道数据。

有益效果：1、本发明将脑电信号进行EEMD分解后得到的固有模态函数进行FIR 滤波器滤波，再经过CSP进行特征提取，在CSP的基础上加入EEMD的频域信息，很好地解决CSP缺乏频域信息的问题。

2、本发明将集合经验模式分解后的多阶固有模态函数重构为新的信号，每一个IMF 经过FIR滤波器滤波得到多个子带信号，IMFs相同频段的信号组合成一个信号集，提取每一个子带信号集的CSP特征，并将所有子带信号集的CSP特征向量进行组合形成新的特征矩阵，得到更具判别性的特征，获得较好的特征分类结果。

3、进行公共空间模式分解，在仅使用C3，C4两个通道的情况下，获得较好的特征分类结果，解决一般CSP算法需要大量输入通道问题。

附图说明

图1本发明方法流程图；

图2 EEG信号采集流程图；

图3为左手运动想象IMF分量频谱图；其中(1)为C3通道前4阶IMFs及频谱图；(2)为C3通道后4阶IMFs及频谱图；(3)C4通道前4阶IMFs及频谱图；(4)为C4通道后4阶IMFs及频谱图；

图4为右手运动想象IMF分量频谱图；其中(1)为C3通道前4阶IMFs及频谱图；(2)为C3通道后4阶IMFs及频谱图；(3)C4通道前4阶IMFs及频谱图；(4)为C4通道后4阶IMFs及频谱图；

图5为脑地形图；其中(1)为左手运动想象第一阶IMF能量的脑地形图；(2)为右手运动想象第一阶IMF能量的脑地形图；(3)为左手运动想象第二阶IMF能量的脑地形图；(4)为右手运动想象第二阶IMF能量的脑地形图；(5)为左手运动想象第三阶IMF 能量的脑地形图；(6)为右手运动想象第三阶IMF能量的脑地形图。

具体实施方式

下面结合说明书附图对本发明作进一步的详细说明。

如图1所示，本发明所述方法包括如下步骤：

步骤1：选择BCICompetition 2008 data sets 2b的数据集作为实验数据，9个受试者进行实验，记录脑电信号的三个电极分别为C3，Cz和C4，选择250Hz的采样频率。其中Cz为参考电极，而C3、C4电极位置则包含了想象手运动最丰富的信息，本发明使用C3，C4通道的数据进行分析。

选择C3、C4电极不仅包含运动想象最丰富的信息，并且减少了通道的数目，降低了数据处理的复杂性，解决了一般CSP算法需要大量输入通道问题。所以分别对单个受试C3、C4两个通道中的信号进行预处理；包括ICA去伪迹、8-32Hz带通滤波；

步骤2：对预处理后的EEG信号x(t)进行EEMD；得到一系列固有模态函数IMF_i(i 为固有模态函数的阶数)并绘制所有固有模态函数能谱图；图3和图4为受试01左右手想象运动时，单次trial的脑电信号EEMD分解后的IMFs与频谱图。

EEMD算法包括如下步骤：

(1)：设置EEMD的一个集合值M。原信号记录为：

x(t)＝[x₁(t),x₂(t),…,x_n(t)]∈R^N×n，其中n代表通道数目，N代表样本点。

(2)：产生白噪声将白噪声增加到信号段x_i(t),i＝(1,2,…n)，得到即i＝1代表EEMD过程的开始。

(3)：对得到的不同的进行经验模式分解(EMD)得到一系列的固有模态函数IMFs，即：x_i,j,1(t),x_i,j2(t),x_i,j,3(t),…x_i,j,k(t)。其中，x_i,j,l(t)l∈(1.2.…,k)代表第i个通道 EEG信号的第j次EMD得到的第l个IMF，j∈(1,2,…M)。

(4)：将得到的IMFs组成矩阵形式得到：

x_i,j(t)＝[x_i,j,1(t),x_i,j2(t),,…x_i,j,k(t)] (1)

(5)：重复步骤(2)-(4)，直到j＝M。计算所得到的M个矩阵的平均值，即：

通过使用EEMD，每一个通道的EEG信号x_i(t)i＝(1,2,…,n)被分解为k个IMFs，即：

x_i,j(t)∈R^N×1(j＝1,2,…k) (3)

因此有：

x(t)＝[x_1,1(t),x_1,2(t),…,x_1,k(t),…x_n,1(t),x_n,2(t),…x_n,k(t)] (4)

步骤3：在一组实验中，重构后的左手运动想象C3与C4通道信号分别为3*2000 的矩阵，组合为6*2000的矩阵，其中6为IMF的个数，可以看作通道数，2000为采样点数。经过6HZ的四个带通滤波器组成的FIR滤波器后得到四个子带信号集为四个 6*2000的矩阵，为新的重构信号，右手运动想象数据处理过程相同。

FIR滤波器的步骤如下：

(1)信号经过EEMD之后，选取第i个通道的前g(g≤k)个IMFs，其中每一个IMF 被FIR滤波器分解为q个子带信号。

因此有：x_i,j(t)＝[x_i,j,1(t),x_i,j2(t),…,x_i,j,q(t)]∈R^N×p故第i个通道的信号x_i(t)可以近似表示为：x_i(t)＝[x_i,1(t),x_i,2(t),…x_i,qg(t)]∈R^N×(qg)。

(2)将选取的g个IMF对应的相同频段的信号重构成一个新的矩阵。因此信号最终重构为q个矩阵。

(3)将q个矩阵分别经过CSP进行特征提取，将得到的q个不同的特征向量组合成一个新的特征矩阵。

步骤4：每次运动想象处理后的信号相同频段对应的子带信号集作为CSP的输入进行特征提取得到一个特征向量，故一组实验可以得到四个特征向量组成的特征矩阵。一个session包含了120次实验，运动想象后得到60组向量矩阵，分为40组训练向量矩阵和20组测试向量矩阵，对每一组经过上述处理后的信号进行CSP分解。

CSP分解具体算法步骤如下：

(1)其中X_H和X_F代表两类运动想象的预处理的脑电图信号，维度为N*T,其中N 代表通道数目，T代表每个通道的样本点。EEG信号的标准化空间协方差矩阵为：

其中X^T代表X的转置，trace(A)为矩阵A对角元素的总和。

(2)Y^k＝{H,F}代表第K次试验的类标签。

两类的平均协方差矩阵计算为：

(3)复合协方差矩阵对角化分解为：

其中U₀为矩阵R的特征向量矩阵，∑为对应的特征值的对角化矩阵。

(4)白化矩阵：

P＝∑^-1/2U₀ ^T (8)

将平均协方差矩阵白化进行白化：

(5)S_H与S_F有着共同的特征向量，两者所对应特征值构成的矩阵和为单位矩阵。

S_F＝U∑_FU^T S_H＝U∑_HU^T (10)

∑_H+∑_F＝I (11)

S_H最大特征值的特征向量对应S_F的最小特征值，反之亦然。白化的EEG信号在对应于∑_H与∑_F最大特征值的特征向量上进行变换，可以得到两类信号矩阵的最优分离方差。

(6)投影矩阵W表示为:

W＝U^TP (12)

(7)利用投影矩阵W将原始EEG信号进行转化:

Z₀＝WX (13)

矩阵Z＝[z₁,z₂,……,z_2m]∈R^N×2m由Z₀的前和后m行组成，定义特征向量为：

f＝[f₁,f₂.…,f_2m]^T∈R^2m×1 (14)

其中：

步骤4：使用支持向量机进行特征分类，其具体步骤如下：

(1)设置已知的训练集：T＝{(x_i,y_i)},i＝1,2,……l，其中，x_i∈Rⁿ，y_i∈{1,-1}。

(2)构建并解决最优化问题：

得到最优解a^*＝(a_i ^*,……,a_l ^*)^T。

(3)计算选择a^*的正分量a_j ^*来计算

(4)构建最优分离超平面为w^*x+b^*＝0，并得到判别函数：

f(x)＝sgn(w^*x+b^*) (17)

SVM算法是最大化两类间隔的线性判别方法，可以提高分类器的分类精度。选择最优分类平面时，应该使离超平面最近的数据点到超平面的距离尽可能的大，即要求将数据正确分类的同时得到的分类间隔最大。核参数γ和误差惩罚因子C是影响SVM性能的主要参数。γ的取值影响空间变换后的数据分布，惩罚因子C则决定了支持向量机的收敛速度及推广能力；因此，对γ和C的选择很大程度上影响了脑电信号的识别率。

(5)采用网格化交叉验证方法进行γ和C最优参数的选择，将训练集作为原始的数据集，在一定范围内改变核函数和惩罚因子的值，运用交叉验证方法进行分类，选择分类准确率最高的γ和C作为最佳参数。

(6)确定γ和C后，将测试集输入训练好的支持向量机进行特征分类得到分类精度。

(7)图5为受试01左右手想象运动时，单次trial的脑电信号EEMD分解后前三阶IMF分量经过CSP空域滤波后特征的脑地形图，由图可见：在进行两类想象运动时， C3、C4电极能量相差较大，同时，想象左手运动和时C4电极明显比C3电极活跃，想象右手运动时则刚好相反，表现为C3电极比C4活跃，这与运动想象中时间按相关同步/去同步特性刚好相符。由此说明，该特征可以有效表征出“源”的空间位置信息，可以作为分类特征。

综上所述，仅为本发明较佳的具体实施方式，但本发明的保护范围并不局限于此。在发明所披露的技术范围内，可轻易想到的变化或替换，都应涵盖在本发明所揭露的技术范围之内。因此，本发明的保护范围应以权利要求书的保护范围为准。

Claims (6)

1.一种改进的EEMD与CSP结合的脑电信号特征提取方法，其特征在于，包括如下步骤：

步骤1：对所采集脑电信号进行预处理，所述脑电信号集由左手运动想象信号和右手运动想象信号组成；

步骤2：对预处理后的脑电信号进行集合经验模式分解，得到若干个固有模态函数；

步骤3：对步骤2所得脑电信号进行FIR滤波，其中每一个固有模态函数被分解为q个子带信号，各个固有模态函数中频段相同的子带信号构成一个子带信号集，得到q个子带信号集；

步骤4：将q个子带信号集分别经过公共空间模式进行特征提取，得到q个特征向量，将q个特征向量组合为一个特征矩阵；

步骤5：将若干个脑电信号的特征矩阵划分为训练集和测试集，利用所述训练集训练支持向量机，并用所述测试集进行验证；

步骤6：利用步骤1-4得到待分类脑电信号的特征矩阵，输入步骤5所得支持向量机进行特征分类。

2.根据权利要求1所述一种改进的EEMD与CSP结合的脑电信号特征提取方法，其特征在于，步骤1中预处理包括ICA去伪迹、8-32Hz带通滤波。

3.根据权利要求1所述一种改进的EEMD与CSP结合的脑电信号特征提取方法，其特征在于，步骤2中集合经验模式分解的步骤包括：

步骤2.1：原始脑电信号记为：x(t)＝[x₁(t),x₂(t),…,x_n(t)]∈R^N×n，其中n代表通道数目，N代表样本点，x_n(t)代表第n个通道的脑电信号；

步骤2.2：产生白噪声将白噪声增加到信号段x_i(t)，i＝(1,2,…n)，x_i(t)代表第i个通道的脑电信号，得到即i＝1代表集合经验模式分解过程的开始；

步骤2.3：对进行经验模式分解得到k个固有模态函数IMFs，即：x_i,j,1(t),x_i,j,2(t),x_i,j,3(t),…x_i,j,k(t)，其中，x_i,j,l(t)，l∈(1.2.…,k)代表第i个通道脑电信号的第j次经验模式分解得到的第l个IMF，j∈(1,2,…M)，M代表经验模式分解进行的总次数；

步骤2.4：将得到的IMFs组成矩阵形式得到：x_i,j(t)＝[x_i,j,1(t),x_i,j2(t),,…x_i,j,k(t)]；

步骤2.5：重复步骤2.2-2.4，直到j＝M，计算所得到的M个矩阵的平均值，即：

为第i个通道经过集合经验模式分解后的IMFs矩阵；x_i,j(t)表示第i个通道第j次经验模式分解后得到的IMFs矩阵；

经集合经验模式分解后，每一个通道的脑电信号x_i(t)，i＝(1,2,…,n)被分解为k个IMFs，即：

x_i,r(t)∈R^N×1(r＝1,2,…k)

因此有：

x(t)＝[x_1,1(t),x_1,2(t),…,x_1,k(t),…x_n,1(t),x_n,2(t),…x_n,k(t)]。

4.根据权利要求1所述一种改进的EEMD与CSP结合的脑电信号特征提取方法，其特征在于，步骤3中对脑电信号进行FIR滤波包括如下步骤：

步骤3.1：脑电信号经过集合经验模式分解之后，选取第i个通道的前g(0＜g≤k)个IMFs，其中每一个IMF被FIR滤波器分解为q个子带信号；

因此有：

x_i,j(t)＝[x_i,j,1(t),x_i,j2(t),…,x_i,j,q(t)]∈R^N×q

故第i个通道的信号x_i(t)表示为：

x_i(t)＝[x_i,1(t),x_i,2(t),…x_i,qg(t)]∈R^N×(qg)；

步骤3.2：将选取的g个IMF对应的相同频段的信号重构成一个新的矩阵，因此信号最终重构为q个矩阵。

5.根据权利要求1所述一种改进的EEMD与CSP结合的脑电信号特征提取方法，其特征在于，步骤4中将q个子带信号集分别经过公共空间模式进行特征提取包括如下步骤：

步骤4.1：子带信号集中脑电信号的标准化空间协方差矩阵为：

其中X_H和X_F分别代表预处理后的左手运动想象信号和右手运动想象信号，维度为n*N，其中n代表通道数目，N代表每个通道的样本点，X^T代表X的转置，trace(*)为矩阵*对角元素的总和；

步骤4.2：平均协方差矩阵计算为：

其中Y^k＝{H,F}代表第k次试验的类标签，H和F分别表示左手运动想象和右手运动想象；

步骤4.3：复合协方差矩阵对角化分解为：

其中U₀为矩阵R的特征向量矩阵，∑为对应的特征值的对角化矩阵；

步骤4.4：白化矩阵：

P＝∑^-1/2U₀ ^T

将平均协方差矩阵白化进行白化：

步骤4.5：S_H与S_F有着共同的特征向量U，两者所对应特征值构成的矩阵和为单位矩阵I；

S_F＝U∑_FU^T S_H＝U∑_HU^T

∑_H+∑_F＝I

S_H最大特征值的特征向量对应S_F的最小特征值，反之亦然；白化的脑电信号在对应于∑_H与∑_F最大特征值的特征向量上进行变换，得到两类信号矩阵的最优分离方差；

步骤4.6：投影矩阵W表示为：

W＝U^TP

步骤4.7：利用投影矩阵W将原始脑电信号进行转化：

Z₀＝WX

矩阵Z＝[z₁,z₂,……,z_2m]∈R^N×2m由Z₀的前和后m行组成，定义子带信号集的特征向量为：

f＝[f₁,f₂.…,f_2m]^T∈R^2m×1

其中：

6.根据权利要求1所述的一种改进的EEMD与CSP结合的脑电信号特征提取方法，其特征在于，步骤1中所述脑电信号为BCI Competition 2008data sets 2b数据集中的C3、C4通道数据。