CN110163128A - 改进的emd算法结合小波包变换及csp算法的脑电信号分类方法 - Google Patents

Abstract

本发明提供一种改进的EMD算法结合小波包变换及CSP算法的脑电信号分类方法，该方法首先对预处理后信号进行小波包变换，然后对重构后的窄带信号进行经验模式分解，根据各阶固有模态函数的相关系数筛选出集中在想象运动频段的固有模态函数进一步进行CSP滤波，通过支持向量机进行特征选择得到最终分类结果。本发明在经验模式分解之前利用小波包变换对EEG信号进行频域滤波，可有效提高EEG的信噪比。本发明利用C3,C4通道EMD分解后的固有模态函数进行CSP滤波，在CSP的基础上加入EMD的频域信息，很好地解决CSP缺乏频域信息的问题。

Description

改进的EMD算法结合小波包变换及CSP算法的脑电信号分类 方法

技术领域

本发明涉及智能信息处理领域，特别涉及一种改进的EMD算法结合小波包变换及CSP算法的脑电信号分类方法。

背景技术

脑机接口(Brain Computer Interface，BCI)，是涉及众多学科和知识领域的控制技术。脑机接口是在人脑与外部设备间进行直接的信息传递，通过采集与提取大脑产生的脑电信号来判别人的意图。脑机接口技术的发展发展不仅帮助瘫痪病人使用计算机、神经假体、机器臂等电子设备，也实现了包括：运动恢复、通信、环境控制甚至娱乐等其他功能。

脑机接口技术主要包含五个步骤，为信号采集、预处理、特征提取、特征分类和接口设备控制。特征提取可以识别不同想象运动的脑电信号的判别信息，因此在脑机接口的研究领域备受关注。传统的时间-频率方法有：短时傅里叶变换(Short-time FourierTransform,STFT)，小波变换(Wavelet Transform,WT)等，但是由于EEG信号非线性和非平稳性，这些方法不可能同时得到时间-频率的良好分辨率。近年，希尔伯特黄变换(Hilbert-Huang Transform,HHT)作为另一种时间-频率分析法十分适合分析非线性和非平稳信号。原信号经过经验模式分解(Empirical Mode Decomposition,EMD)被分解为一系列固有模态函数(Intrinsic Mode Functions,IMFs)，随后对每个固有模态函数进行希尔伯特黄变换，求其相应的能量谱和边际谱，作为特征进行分类。HHT不涉及海森堡不确定性原理可以获得时域和频域的高分辨率，目前它被应用于多个领域，如雷达探测，地震信号和生物医学信号等。

同时由于EEG信号的空间分辨率低，其构成的BCI系统需要进行有效的空间滤波，从而确保能从受试者的相关脑域中提取特征信息。常用的算法有：共空间模式 (CommonSpatial Pattern,CSP)、独立主成分分析(Independent Component Analysis, ICA)、滤波器CSP(Filter Bank Common Spatial Pattern,FBCSP)、判别滤波CSP (DiscriminantFiltering Common Spatial Pattern,DFBCSP)等多种CSP的改进算法。

但是传统的经验模式分解在低频区域首次得到的固有模态函数涵盖太宽的频率范围，无法得到单组份的属性，而且不能分开含有低能量成分的信号。为了弥补这些不足，本发明提出了一种改进方法，即在经验模式分解之前，使用小波包变换将信号分为一组窄带信号，并利用筛选过程来从结果中去除不相关的固有模态函数。同时，常规的公共空间模式分解方法需要大量的输入通道、缺乏频域信息，发展受到限制。为了克服这些缺点，将经验模式分解和公共空间模式算法结合，改变CSP滤波器成分选择方式，筛选有效的IMF频段，使用改进的CSP滤波器进行滤波获取特征。本发明提出的小波包变换结合经验模式分解和共空间模式的方法，很好地解决了上述问题，为更好地进行想象运动判别提供了一个可选的途径。

发明内容

本发明目的在于针对传统经验模式分解及公共空间模式算法的不足，提出了一种经验模式分解结合小波包变换和公共空间模式算法的脑电信号特征提取方法。该方法是基于如下内容：1、基于小波包变换的EEG信号处理；2、重构后的信号进行经验模式分解得到固有模态函数；3筛选固有模态函数，并进行公共空间模式分解，解决CSP多输入、缺频域信息的问题；4、支持向量机分类。

本发明解决其技术问题所采取的技术方案是：本发明的系统首先对预处理后信号进行小波包变换，然后对重构后的窄带信号进行经验模式分解，根据各阶固有模态函数的相关系数筛选出集中在想象运动频段的固有模态函数进一步进行CSP滤波，通过支持向量机进行特征选择得到最终分类结果。

本发明在传统经验模式分解EMD算法及常规的公共空间模式分解CSP上做了如下改进：

在使用经验模态分解EMD之前，将C3,C4通道的数据进行小波分解，得到一系列的窄带信号，再将合适子带的信号进行EMD分解，得到多个本征模态函数IMFs，取其包含μ和β节律的前两阶作为新的输入信号，则相当于共有8个通道的数据，构成一个 8*2000的矩阵Xi(i＝L表示想象左手运动，i＝R表示想象右手运动)，其中，原始脑电信号经小波分解，选取两子带信号(4，1)，(4，3)，经EMD分解，选取前两阶IMFs，共8个，2000为一次试验的采样点个数，进行公共空间模式分解。

方法流程：

步骤1：选取9位受试的脑电信号作为训练集和测试集，分别对单个受试者C3、 C4通道的数据作为分析，采样频率为250Hz，采样数据通过0.5-100Hz带通滤波器滤波进行预处理，50Hz陷波滤波器。

步骤2：对预处理后的EEG信号进行小波包变换过程。由于数据集EEG信号的频带宽为100Hz。采用db4小波基对EEG信号进行小波包4层分解，分解后得到16个频率成分，将信号分解为若干窄带信号。对选取的小波包子带进行重构，从而提取到相关性较强的脑电，很好的滤除噪声的干扰，保留了运动想象脑电信息。

步骤3：对处理后的EEG信号进行经验模式分解；得到一系列固有模态函数IMF_i(i为固有模态函数的阶数)；

EEG信号进行经验模式分解的具体步骤如下：

(1)判断每个x(t)的局部极值，用三次样条曲线进行曲线拟合，局部极大值形成上包络e_max(t),局部极小值形成下包络e_min(t)。

(2)求e_max(t)和e_min(t)的均值：

(3)计算输入信号x(t)和m(t)的差值：

c(t)＝x(t)-m(t) (2)

如果c(t)不能满足IMF的定义截至条件，重复上述过程(1)-(3)，否则，提取c(t)作为固有模态函数，剩余量r(t)计算如下：

r(t)＝x(t)-c(t) (3)

(5)剩余量作为一个新的数据经过相同的筛选过程以获得下一个更低频率的固有模态函数。直到剩余函数r(t)为一个单调函数或者仅有一个极致时，分解过程停止。假设原始信号x(t)被分解为n个固有模态函数和一个剩余函数量r(t)，可以得重构信号：

步骤4：信号进行EMD分解后，得到多个本征模态函数IMFs，取其包含μ和β节律的频段，并筛选前两阶作为新的输入信号，则相当于共有8个通道的数据，构成一个8*2000的矩阵Xi(i＝L表示想象左手运动，i＝R表示想象右手运动)，即原始脑电信号经小波分解，选取两子带信号经EMD分解，选取前两阶IMFs，进行公共空间模式分解。

公共空间模式算法具体过程如下：

(1)X_L和X_R表示预处理过的EEG矩阵(下标L和R分别代表左右想象运动和右手想象运动)，N*T维，N为通道数，T为采样点数。标准化的空间协方差表示为：

(2)对复合协方差矩阵进行对角化分解：

其中∑为特征值对角矩阵，U₀为其对应的特征向量矩阵。将∑的特征值按降序排列，其对应的U₀也重新排列。

白化矩阵：

对平均协方差矩阵进行变形：

S_L和S_R有共同的特征向量，满足下式：

S_L＝B∑_LB^T (11)

S_R＝B∑_RB^T (12)

∑_L+∑_R＝I (13)

(3)S_L最大特征值对应的特征向量对应S_R最小特征值，反之亦然。对应白化EEG 的最大特征值的特征向量进行变换，可以获得两个信号矩阵的最优分离方差。投影矩阵 W表示为：

W＝U^TP (14)

此处W即空间滤波器，原始的EEG信号x(t)经过W可以被投影得到新的数据集Z₀，如式(13)所示：

Z₀＝x(t)·W (15)

在传统的CSP算法中，构建矩阵Z＝[z₁,z₂,…,z_2m]∈R^N×2m由Z₀的前行和后m行构成，特征向量f＝[f₁,f₂,…,f_2m]^T∈R^2m×1即定义为：

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

(1)对于非线性EEG问题，将特征集通过非线性变换转化为另一个空间中的线性问题，构造最优分类面。相应的最优决策函数为：

其中N为支持向量个数，αi为Lagrangue乘子。从而目标函数变为使下式最小化：

核参数γ和误差惩罚因子C是影响SVM性能的主要参数。γ的取值影响空间变换后的数据分布，惩罚因子C则决定了支持向量机的收敛速度及推广能力；因此，对γ和 C的选择很大程度上影响了脑电信号的识别率。

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

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

步骤6：将待分类脑电信号输入已训练好的支持向量机中，得到其分类结果。

有益效果：1、本发明在经验模式分解之前利用小波包变换对EEG信号进行频域滤波，可有效提高EEG的信噪比。小波包是在小波变换的基础上发展而来，在信号分析中，小波包比小波变换更为精细。它不仅能够将信号在全频段内进行正交分解，而且可以根据信号特征自适应地选择相应频段，使之与信号频谱相匹配，具有更高的时频分辨率。利用小波包变换的多分辨率特性，可以选择EEG原始信号中的最佳分量组合关系，从而提取出包含有用信息频率范围的信号。

2、本发明利用C3,C4通道EMD分解后的固有模态函数进行CSP滤波，在CSP的基础上加入EMD的频域信息，很好地解决CSP缺乏频域信息的问题。

3、本发明将经验模式分解后的多阶固有模态函数看作多输入信号进行公共空间模式分解，在仅使用C3，C4通道的情况下，获得较好的特征分类结果，解决一般CSP算法需要大量输入通道问题。

附图说明

图1本发明方法流程图。

具体实施方式

下面结合实施例对本发明做进一步说明。

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

步骤1：选取9位受试的脑电信号作为训练集和测试集，分别对单个受试者C3、 C4通道的数据作为分析，采样频率为250Hz，采样数据通过0.5-100Hz带通滤波器滤波进行预处理，50Hz陷波滤波器。

步骤2：对预处理后的EEG信号进行小波包变换过程。由于数据集EEG信号的频带宽为100Hz。采用db4小波基对EEG信号进行小波包4层分解，分解后得到16个频率成分。基于频率，小波包变换后将信号分解为若干窄带信号。最小分辨率Δf＝100 Hz/16＝6.25Hz。因此，对小波包分解系数进行重构，得到的小波包各个子带信号所对应的频率范围有，子带(4，0)对应0～6.25Hz，子带(4，1)对应6.25～12.5Hz，子带 (4，2)对应12.5～18.75Hz，子带(4，3)对应18.75～25Hz等等，其中，μ节律脑电包含在小波包子带(4，1)中，β节律脑电包含在小波包子带(4，3)中，因此，对小波包子带(4，1)、(4，3)进行重构，从而提取到μ节律和β节律脑电，为后期的特征提取和分类识别提供有力保证。

步骤3：对处理后的EEG信号进行经验模式分解；得到一系列固有模态函数IMF_i(i为固有模态函数的阶数)；

EEG信号进行经验模式分解的具体步骤如下：

(1)判断每个x(t)的局部极值，用三次样条曲线进行曲线拟合，局部极大值形成上包络e_max(t)，局部极小值形成下包络e_min(t)。

(2)求e_max(t)和e_min(t)的均值：

(3)计算输入信号x(t)和m(t)的差值:

c(t)＝x(t)-m(t) (2)

如果c(t)不能满足IMF的定义截至条件，重复上述过程(1)-(3)，否则，提取c(t)作为固有模态函数，剩余量r(t)计算如下：

r(t)＝x(t)-c(t) (3)

(6)剩余量作为一个新的数据经过相同的筛选过程以获得下一个更低频率的固有模态函数。直到剩余函数r(t)为一个单调函数或者仅有一个极致时，分解过程停止。假设原始信号x(t)被分解为n个固有模态函数和一个剩余函数量r(t)，可以得重构信号：

步骤4：将合适子带的信号进行EMD分解，得到多个本征模态函数IMFs，取其包含μ和β节律的前两阶作为新的输入信号，选取两子带信号(4，1)，(4，3)，经EMD 分解，选取前两阶IMFs，共8个，2000为一次试验的采样点个数，进行公共空间模式分解。

公共空间模式算法包括如下步骤：

(1)X_L和X_R表示预处理过的EEG矩阵(下标L和R分别代表左右想象运动和右手想象运动)，N*T维，N为通道数，T为采样点数。标准化的空间协方差表示为：

(2)对复合协方差矩阵进行对角化分解：

其中∑为特征值对角矩阵，U₀为其对应的特征向量矩阵。将∑的特征值按降序排列，其对应的U₀也重新排列。

白化矩阵：

对平均协方差矩阵进行变形：

S_L和S_R有共同的特征向量，满足下式：

S_L＝B∑_LB^T (11)

S_R＝B∑_RB^T (12)

∑_L+∑_R＝I (13)

(3)S_L最大特征值对应的特征向量对应S_R最小特征值，反之亦然。对应白化EEG 的最大特征值的特征向量进行变换，可以获得两个信号矩阵的最优分离方差。投影矩阵 W表示为：

W＝U^TP (14)

此处W即空间滤波器，原始的EEG信号x(t)经过W可以被投影得到新的数据集Z₀，如式(13)所示：

Z₀＝x(t)·W (15)

在传统的CSP算法中，构建矩阵Z＝[z₁,z₂,…,z_2m]∈R^N×2m由Z₀的前行和后m行构成，特征向量f＝[f₁,f₂,…,f_2m]^T∈R^2m×1即定义为：

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

(1)对于非线性EEG问题，将特征集通过非线性变换转化为另一个空间中的线性问题，构造最优分类面。相应的最优决策函数为：

其中N为支持向量个数，αi为Lagrangue乘子。从而目标函数变为使下式最小化：

核参数γ和误差惩罚因子C是影响SVM性能的主要参数。γ的取值影响空间变换后的数据分布，惩罚因子C则决定了支持向量机的收敛速度及推广能力；因此，对γ和 C的选择很大程度上影响了脑电信号的识别率。

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

(3)确定γ、C，将测试集输入训练好的支持向量机进行特征分类。

步骤6：将待分类脑电信号输入已训练好的支持向量机中，得到其分类结果。

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

Claims (4)

1.改进的EMD算法结合小波包变换及CSP算法的脑电信号分类方法，其特征在于，包括如下步骤：

步骤1：采集左右手运动的C3、C4通道的脑电信号，并对脑电信号进行预处理；

步骤2：对预处理后的脑电信号进行小波包变换，得到一系列的窄带信号；

步骤3：对窄带信号进行经验模式分解，得到一系列固有模态函数IMF_i，i为固有模态函数的阶数；

步骤4：根据各阶固有模态函数的相关系数筛选出集中在想象运动频段的固有模态函数进行CSP滤波，得到脑电信号的特征向量；

步骤5：将特征向量分为训练集和测试集训练支持向量机；

步骤6：将待分类脑电信号输入已训练好的支持向量机中，得到其分类结果。

2.根据权利要求1所述的改进的EMD算法结合小波包变换及CSP算法的脑电信号分类方法，其特征在于，

步骤3中，对脑电信号进行经验模式分解的具体步骤如下：

(1)判断每个脑电信号x(t)的局部极值，用三次样条曲线进行曲线拟合，局部极大值形成上包络e_max(t)，局部极小值形成下包络e_min(t)；

(2)求e_max(t)和e_min(t)的均值：

(3)计算输入信号x(t)和m(t)的差值：

c(t)＝x(t)-m(t) (2)

如果c(t)不能满足IMF的定义截至条件，重复上述过程(1)-(3)，否则，提取c(t)作为固有模态函数，剩余量r(t)计算如下：

r(t)＝x(t)-c(t) (3)

(4)剩余量作为一个新的数据经过相同的筛选过程以获得下一个更低频率的固有模态函数；直到剩余函数r(t)为一个单调函数或者仅有一个极致时，分解过程停止；假设原始信号x(t)被分解为n个固有模态函数和一个剩余函数量r(t)，可以得重构信号：

3.根据权利要求1所述的改进的EMD算法结合小波包变换及CSP算法的脑电信号分类方法，其特征在于，步骤4中，选取前两阶IMFs进行公共空间模式分解。

4.根据权利要求1所述的改进的EMD算法结合小波包变换及CSP算法的脑电信号分类方法，其特征在于，步骤4中，公共空间模式算法具体过程如下：

(1)X_L和X_R表示预处理过的EEG矩阵，下标L和R分别代表左右想象运动和右手想象运动，N*T维，N为通道数，T为采样点数；标准化的空间协方差表示为：

(2)对复合协方差矩阵进行对角化分解：

其中∑为特征值对角矩阵，U₀为其对应的特征向量矩阵，将∑的特征值按降序排列，其对应的U₀也重新排列；

白化矩阵：

对平均协方差矩阵进行变形：

S_L和S_R有共同的特征向量，满足下式：

S_L＝B∑_LB^T (11)

S_R＝B∑_RB^T (12)

∑_L+∑_R＝I (13)

(3)S_L最大特征值对应的特征向量对应S_R最小特征值，反之亦然；对应白化EEG的最大特征值的特征向量进行变换，可以获得两个信号矩阵的最优分离方差；投影矩阵W表示为：

W＝U^TP (14)

此处W即空间滤波器，原始的EEG信号x(t)经过W可以被投影得到新的数据集Z₀，如式(13)所示：

Z₀＝x(t)·W (15)

在传统的CSP算法中，构建矩阵Z＝[z₁,z₂,…,z_2m]∈R^N×2m由Z₀的前行和后m行构成，特征向量f＝[f₁,f₂,…,f_2m]^T∈R^2m×1即定义为：