CN108280464A - Dwt和emd融合近似熵的脑电特征提取方法 - Google Patents

Abstract

本发明公开了DWT和EMD融合近似熵的脑电特征提取方法，该方法对预处理后的信号进行小波分解(DWT)，得到各个子带信号，选择频带合适的子带进行经验模式(EMD)分解，得到固有模态函数(Intrinsic Mode Functions,IMFs)，观察并计算每个IMF分量的能量谱，筛选有效的IMF频段(5‑28Hz)，组成新的信号矩阵，对其进行近似熵计算，最后使用支持向量机(Support Vector Machine,SVM)进行分类。分类结果得到9位受试的想象运动平均分类正确率为90％以上，保证了该方法的可行性与有效性。

Description

DWT和EMD融合近似熵的脑电特征提取方法

技术领域

本发明涉及智能信息处理领域，特别涉及一种基于DWT和EMD，同时融合近似熵的脑电信号特征提取方法。

背景技术

传统的运动通道是由大脑神经和肌肉组成，神经传导冲动，肌肉配合完成相应的动作，而脑-机接口(Brain-Computer Interface,BCI)则提供了另一种运动通道，不依赖传统的运动通道，由大脑意识直接与外部设备相连，建立运动通道，用人的大脑意识直接控制外部设备，不需要经过神经传导和肌肉运动，为神经损伤或肌肉损伤的病人提供了一种新的运动方式，不再需要依靠别人的照顾，自己即可完成运动。脑机接口技术领域的发展不仅可以帮助瘫痪病人使用计算机、神经假体、机械臂等电子设备,也实现了包括:运动恢复、通信、环境控制甚至娱乐等其他功能。

脑机接口技术主要包含信号采集、预处理、特征提取、特征分类和接口设备控制等五个步骤。其中，由特征提取得到的特征信号可以标识不同想象运动脑电信号的判别信息，对后续的分类识别有很大的影响，因此特征提取在BCI研究界受到广泛的关注。

有效的特征提取方法是改善识别精度的关键，目前，时间-频率方法被广泛应用于脑信号的研究。传统的时间-频率方法包括：短时傅里叶变换(Short-timeFourierTransform,STFT)，小波变换(WaveletTransform,WT)等，但是这些方法的本质都是基于傅里叶变换，根据海森堡不确定性原理，该方法不可能同时得到时间-频率的良好分辨率。近年，希尔伯特黄变换(Hilbert-Huang Transform,HHT)作为另一种时间-频率分析法已经变得越来越流行，同时它也十分适合分析非线性和非平稳信号。原信号经过经验模式分解(Empirical Mode Decomposition,EMD)被分解为一系列固有模态函数(Intrinsic ModeFunctions,IMFs)，随后对每个固有模态函数进行希尔伯特黄变换，求其相应的能量谱和边际谱，作为特征进行分类。HHT不涉及海森堡不确定性原理可以获得时域和频域的高分辨率。目前被广泛应用于许多信号处理领域，如雷达探测、地震信号和生物医学信号等。

由于生理系统的非线性性质，因此采用非线性分析方法可能有助于更好的揭示其特性和机理。目前多数研究是通过估计信号的某些非线性动力学参数(如关联维数，哥氏熵，李雅普诺夫指数等)来表征信号，但是估计这些参数所需数据点数较多，因此造成这些参数实际应用的困难。近似熵这一概念是Pincus在1991年提出的，很好的解决了该类问题，近似熵只需要较短的数据点就可以估计随机信号和确定信号的特性，具有很好的抗噪和抗干扰能力。对于脑电信号来说，即有很大的随机信号成分，同时又有一些确定的频率成分，因此近似熵对于脑电信号的研究是十分有效和重要。

发明内容

本发明针对上述现有技术的不足，提供一种基于DWT和EMD，同时融合近似熵的脑电特征提取方法，大大提供了模式识别的准确率。该方法是基于如下内容：1)基于小波分解的EEG信号处理；根据频谱筛选子带信号，组成新的信号矩阵；2)基于经验模式分解的EEG信号处理；频谱分析筛选固有模态函数，组成新的信号矩阵；3)对信号矩阵进行近似熵计算；4)支持向量机分类。

本发明解决其技术问题所采取的技术方案是：本发明首先对预处理后信号进行小波分解，根据子带的频谱选择子带信号，再经经验模式分解，根据各阶固有模态函数的频谱分析筛选出集中在想象运动频段的前2阶固有模态函数，对其进行近似熵计算，通过支持向量机进行特征选择得到最终分类结果。

本发明对通道信号的提取进行了改进包括：

对于左手想象运动，提取C3和C4两个通道的数据作为原始信号；对于右手想象运动，只提取C3通道的数据作为原始信号。将右手想象运动的C3通道分别与左手想象运动的C3，C4通道进行对比，提取特征向量，进行分类识别。

本发明对信号的近似熵计算进行了改进包括：

计算8*2000信号矩阵近似熵，其中8可看作通道数(2(C3，C4通道)*2(D4，D3子带)*2(2阶IMFs))，2000为一次试验的采样点个数，对矩阵每一行数据进行滑动窗计算，以起始点开始，以500个数据点为窗口长度，以1个数据点为滑动长度，计算得到该滑动窗的近似熵，得到8*1501的近似熵矩阵，作为特征向量进行分类识别。

方法流程：

步骤1：选取9位受试的脑电信号作为训练集和测试集，分别对单个受试C3、C4两个通道中的信号进行预处理；

步骤2：对预处理后的EEG信号x(t)进行4阶小波分解,得到一系列子带信号A4，D4，D3，D2，D1，信号x(t)的小波变换定义为：

其中Ψ(t)为小波基函数，选取频率在8-30Hz的子带信号D4，D3；

步骤3：对子带信号D4，D3进行经验模式分解；得到一系列固有模态函数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) 公式(3)

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

r(t)＝x(t)-c(t) 公式(4)

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

步骤4：将单次试验的C3，C4通道前2阶IMF分量合并，构成一个8*2000的矩阵X_i，i＝L表示想象左手运动，i＝R表示想象右手运动，其中，8可看作通道数(2(C3，C4通道)*2(D4，D3子带)*2(2阶IMFs))，2000为一次试验的采样点个数，即窗口长度。对矩阵进行近似熵计算；

近似熵算法具体过程如下：

设原始数据为x(1),x(2),…,x(N),共N个点。

步骤1：按序号连续顺序组成一组m维矢量：X(i)＝[x(i),x(i+1),…,x(i+m+1)]，其中i＝1～N-m+1；

步骤2：定义X(i)与X(j)间的距离d[X(i),X(j)]为两者对应元素中差值最大的一个，即：

步骤3：给定阈值r，对每一个i值统计d[X(i),X(j)]小于r的数目及此数目与距离总数N-m+1的比值。记作即：

其中i＝1～N-m+1；

步骤4：先将取对数，再求其对所有i的平均值，记作即：

步骤5：再把维数加1，变成m+1,重复步骤1至步骤4，得和

步骤6：理论上此序列的近似熵为：

一般而言，此极限值以概率1存在，但在实际中N不可能为∞，当N为有限值时，按上述步骤得出的是序列长度为N时近似熵值ApEn的估计值，记作

ApEn(m,r)＝φ^m(r)-φ^m+1(r) 公式(10)

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

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

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

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

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

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

有益效果：

1、本发明利用三个通道小波分解后的子带信号进行EMD分解，对固有模态函数进行近似熵计算，在EMD的基础上加入DWT，很好地解决EMD分解时频带涵盖宽泛的问题。

2、本发明使用滑动窗口对信号矩阵进行分段近似熵计算，再组合各个IMF的近似熵作为特征进行分类识别。

附图说明

图1为本发明方法流程图。

具体实施方式

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

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

步骤1：采集各位受试的脑电信号。选取9位受试的脑电信号作为训练集和测试集，分别对单个受试C3、C4两个通道中的信号进行预处理；

步骤2：对预处理后的EEG信号x(t)进行4阶小波分解,得到一系列子带信号A4，D4，D3，D2，D1，信号x(t)的小波变换定义为：

其中Ψ(t)为小波基函数，选取频率在8-30Hz的子带信号D4，D3；

步骤3：对子带信号D4，D3进行经验模式分解；得到一系列固有模态函数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) 公式(3)

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

r(t)＝x(t)-c(t) 公式(4)

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

步骤4：将单次试验的C3，C4通道前2阶IMF分量合并，构成一个8*2000的矩阵X_i(i＝L表示想象左手运动，i＝R表示想象右手运动)，其中，8可看作通道数(2(C3，C4通道)*2(D4，D3子带)*2(2阶IMFs))，2000为一次试验的采样点个数，即窗口长度。对矩阵进行近似熵计算；

近似熵算法具体过程如下：

设原始数据为x(1),x(2),…,x(N),共N个点。

步骤1：按序号连续顺序组成一组m维矢量：X(i)＝[x(i),x(i+1),…,x(i+m+1)]，其中i＝1～N-m+1；

步骤2：定义X(i)与X(j)间的距离d[X(i),X(j)]为两者对应元素中差值最大的一个，即：

步骤3：给定阈值r，对每一个i值统计d[X(i),X(j)]小于r的数目及此数目与距离总数N-m+1的比值。记作即：

其中i＝1～N-m+1；

步骤4：先将取对数，再求其对所有i的平均值，记作即：

步骤5：再把维数加1，变成m+1,重复步骤1至步骤4，得和

步骤6：理论上此序列的近似熵为：

一般而言，此极限值以概率1存在，但在实际中N不可能为∞，当N为有限值时，按上述步骤得出的是序列长度为N时近似熵值ApEn的估计值，记作

ApEn(m,r)＝φ^m(r)-φm+1(r) 公式(10)

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

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

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

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

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

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

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

Claims (4)

1.DWT和EMD融合近似熵的脑电特征提取方法，其特征在于，所述方法包括如下步骤：

步骤1：选取9位受试的脑电信号作为训练集和测试集，分别对单个受试C3、C4两个通道中的信号进行预处理；

步骤2：对预处理后的EEG信号x(t)进行4阶小波分解,得到一系列子带信号A4，D4，D3，D2，D1，信号x(t)的小波变换定义为:

其中Ψ(t)为小波基函数，选取频率在8-30Hz的子带信号D4，D3；

步骤3：对子带信号D4，D3进行经验模式分解，得到一系列固有模态函数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) 公式(3)

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

r(t)＝x(t)-c(t) 公式(4)

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

步骤4：将单次试验的C3，C4通道前2阶IMF分量合并，构成一个8*2000的矩阵X_i，i＝L表示想象左手运动，i＝R表示想象右手运动，其中，8为通道数(2(C3，C4通道)*2(D4，D3子带)*2(2阶IMFs))，2000为一次试验的采样点个数，即窗口长度，对矩阵信号进行近似熵计算。

2.根据权利要求1所述的DWT和EMD融合近似熵的脑电特征提取方法，其特征在于：所述近似熵算法包括如下步骤：

设原始数据为x(1),x(2),…,x(N),共N个点；

步骤1：按序号连续顺序组成一组m维矢量：X(i)＝[x(i),x(i+1),…,x(i+m+1)]，其中i＝1～N-m+1；

步骤2：定义X(i)与X(j)间的距离d[X(i),X(j)]为两者对应元素中差值最大的一个，即：

步骤3：给定阈值r，对每一个i值统计d[X(i),X(j)]小于r的数目及此数目与距离总数N-m+1的比值，记作即：

其中i＝1～N-m+1；

步骤4：先将取对数，再求其对所有i的平均值，记作即：

步骤5：再把维数加1，变成m+1,重复步骤1至步骤4，得和

步骤6：理论上此序列的近似熵为：

极限值以概率1存在，但在实际中N不可能为∞，当N为有限值时，按上述步骤得出的是序列长度为N时近似熵值ApEn的估计值，记作：

ApEn(m,r)＝φ^m(r)-φ^m+1(r) 公式(10)。

3.根据权利要求1或2所述的DWT和EMD融合近似熵的脑电特征提取方法，其特征在于：所述方法对上述矩阵信号近似熵计算包括：

对8*2000信号矩阵，对矩阵每一行数据进行滑动窗计算，以起始点开始，以500个数据点为窗口长度，以1个数据点为滑动长度，计算得到该滑动窗的近似熵，得到8*1501的近似熵矩阵，作为特征向量进行分类识别。

4.根据权利要求1所述的DWT和EMD融合近似熵的脑电特征提取方法，其特征在于：所述方法使用支持向量机进行特征分类，包括如下步骤：

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

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

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

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

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