基于特征筛选的经验模式分解脑电信号分类方法

张学军， 汪 敏， 胡晓雯

(1.南京邮电大学电子与光学工程学院， 南京 210023； 2.南京邮电大学射频集成与微组装技术国家地方联合工程实验室， 南京 210023；3．南京医科大学生物医学工程与信息学院， 南京 211166)

脑电(electroencephalogram，EEG)信号是一组包含人类大脑活动信息的电位差，展示了有关脑电流的数据。EEG的测量可以通过放置在头皮上的传感器或使用颅内电极来获得。EEG应用领域广泛，如情绪识别、脑-机接口(brain-computer interface，BCI)等[1]。

基于运动想象的BCI以其无创、适用性好、可移植性强等优点吸引了众多研究者的兴趣。BCI为人们提供了一种使用脑电信号与外部辅助设备交互的方式，在生物医学工程和神经修复中有着广泛的应用[2]。首先，受试者需要在大脑中进行基于运动想象(motion imagination，MI)的特定运动，其次，基于MI的BCI会对特定运动的脑电信号进行采集、分类，将其转换为手、脚等不同想象任务的控制信号。

BCI研究的目标是开发帮助残疾用户与他人交流系统。BCI系统是一个连续闭环系统，通常由脑信号获取、预处理、特征提取、分类、输出命令和反馈五个部分组成。预处理方法如滤波或盲源分离算法可以减少不同类型的脑电噪声和眼电尾迹。

由于脑电信号具有噪声强、强度低、对周围环境敏感、数据维数大和分布复杂等特点，在实际应用中还不能很好地分类。为了减轻噪声脑电信号干扰的影响，可以利用各种信号处理方法从脑电信号中提取的参数对脑电信号的分类。基于傅里叶变换的谱参数对脑电信号的分析在脑电信号的分类上取得了良好的效果。然而，傅里叶变换没有时频的局部化性能，从而傅里叶域在信号中不表现出任何时域特征[3]。Gabor提出使用短时傅里叶变换(short time fourier transform，STFT)解决表示时域特征的问题，但STFT并不能对信号进行多分辨率分析，这是因为STFT使用相同带宽的滤波器对所有频率的信号进行分解[4]。通常利用小波分析来解决。在小波分析中，通过形成具有不同带宽的带通滤波器来促进多分辨率时频分析。研究人员发现小波分析对于各种信号处理应用是非常有用的工具。在脑电信号处理方面，Zhang等[5]提出了用于脑电信号处理的小波CSP(common space pattern)算法。姚悦等[6]提出小波变换结合二阶盲辨识的眼电伪迹自动去除方法，其分离效果不受 源信号高斯性影响，在抑制脑电信号眼电尾迹方面取得了良好的效果。尽管上述文献中使用小波进行信号处理，但都是基于被分析信号是平稳的。对脑电信号分析的研究表明，脑电信号的频率成分在一段时间内会发生变化。经验模态分解(empirical mode decomposition，EMD)是一种基于时频的经验方法，对于非平稳信号的时频分析是有效的[7]，它将信号分解为若干固有模函数(intrinsic mode function，IMF)，IMF是振荡分量。

为了控制BCI系统来识别用户的活动，并将其转换为命令，在大多数现有的BCI中依赖于分类算法，不同的分类算法用于不同的BCI应用(例如线性分类器、神经网络、非线性贝叶斯分类器、最近邻分类器)。线性分类器是使用线性函数来分类的算法。BCI系统的设计主要采用两种线性分类器，即线性判别分析(linear discriminant analysis，LDA)和支持向量机(support vector machine，SVM)。但是使用LDA的前提是假设数据是正态分布的，对于EEG数据分类敏感度不高。选用SVM分类器，SVM的主要思想是构造一个最优的超平面，使分割平面和数据之间的差值最大化。使用核函数可以向高维空间进行映射，解决非线性的分类[8]。 SVM最佳参数可以通过网格搜索算法进行调整，但是精度较低。

现提出一种基于小波变换和EMD结合排列熵的特征筛选的分类方法。利用小波变换提取脑电运动想象窄带信号，再进一步EMD提取窄带信号的IMF分量。为了减少冗余特征对算法精度的影响，计算每个IMF的排列熵，进行Perason系数比较筛选出合适的特征，最后运用遗传算法(genetic algorithm, GA)对SVM分类方法进行优化。

1 特征提取与分类方法

1.1 小波分解

小波变换(wavelet transform，WT)是描述非平稳信号的有力数学工具，是一种变换分析方法，能够在时间、空间频率的局部化分析，通过伸缩平移运算对信号逐步进行多尺度细化，最终达到高频处时间细分，低频处频率细分，能自动适应时频信号分析的要求[9]。它具有自适应性和多分辨率能力，因此适合于分解时频分辨率不同的脑电信号。

小波变换是将母小波函数作位移τ后，再在不同的尺度a下，与待分析信号X(t)作内积，即

(1)

式(1)中:a> 0，称为尺度因子;τ反映位移,可正可负；ψ()为母小波函数；WTX为小波变换后的函数；t为时间。离散小波变换DWT对尺度参数按幂级数进行离散化处理，对时间进行均匀离散化取值如二进制离散化尺度时间为2，4，6，…,2n(要求采样率满足尼奎斯特采样定理)。

利用Mallat算法进行小波分解过程如图1所示。

CA为低频信息、近似分量；CD为高频、细节分量

1.2 经验模态分解

EMD具有直观性和自适应性，对EEG等非平稳信号表现良好。EMD的目的是将信号分解为一组固有模函数IMF。IMF定义为极值数和过零点数相等(或最多相差一个)的函数，其包络由所有局部极大值和极小值定义，相对于零对称[10]。IMF表示一个简单的振荡模式，作为傅里叶分析中使用的简单谐波函数的对应项[11]。

给定一个信号x(n)，EMD的起点是识别所有的局部极大值和极小值。所有局部极大值以三次样条曲线作为上包络eu(n)连接，同样，所有局部极小值以样条曲线作为下包络el(n)连接。两个包络的平均值表示为ml(n)=[eu(n)+el(n)]/2，则获得第一个原始IMF分量hl(n)为

hl(n)=x(n)-ml(n)

(2)

上述提取IMF的过程称为筛选过程。由于hl(n)仍然在零交叉点之间包含多个极值，因此对hl(n)再次执行筛选过程。将此过程重复应用于原始IMFhk(n)，直到得到满足固有模态函数条件[12]。

rl(n)=x(n)-hl(n)

(3)

式(3)中：rl(n)为残余分量，包含一些有用的信息。因此，可以将残余当作一个新的信号，并应用上述过程来获得新的IMF。

ri-1(n)-ci(n)=ri(n),i=2，3,…，N

(4)

当残余rN(n)是常数、单调斜率或只有一个极值的函数时，整个过程终止。结合式(3)和式(4)中的方程得到原始信号：

(5)

1.3 排列熵

考虑标量时间序列X(t)(t∈1,2,…,N)，非线性数据分析的第一步是相空间重构，最常用的方法是利用延迟时间嵌入定理。排列熵的原理是引入一个时间延迟和嵌入维度[13]。

在这种方法中，时间序列的值被转换成一个延迟向量：

Xi→[Xi-(d-1)τ，Xi-(d-2)τ，…，

Xi-τ，Xi]

(6)

式(6)中:d为嵌入维度;τ为时间延迟。这会将N个标量转换为具有重叠项的N-τ(d-1)向量。

可以按递增顺序排列d维延迟向量中的值，以实现有序模式:

[xi-rd-1τ≤xi-rd-2τ≤…≤xi-r1τ≤xi-r0τ]

(7)

它们的相等发生在设置rl

根据香农熵原理可得排列熵定义为

(8)

显然，当P(j)的递增或递减序列达到下限时,0

(9)

排列熵是时间序列的一种复杂度量，可以快速、简便地计算有序模式。当处理值之间的顺序关系时，排列熵对于噪声是鲁棒的。

1.4 基于遗传算法的支持向量机参数优化

基于径向基函数(radial basis function，RBF)核函数的支持向量机分类器在分析高维数据方面具有优势，SVM核函数的参数对分类器的性能影响很大，需定义两个参数惩罚因子c和核参数g，目前还缺乏有效确定参数选取的结构方法。因此，将遗传算法应用于所提出的支持向量机模型中，优化参数选择[14]。遗传算法是通过世代搜索而不是单点搜索来获得最优解或准最优解；它具有全局寻优能力；它是种群变化的并行过程，具有内在的并行性。遗传算法的处理对象是参数集被编码的个体，而不是参数本身，这一特点使得遗传算法得到了广泛的应用。

基于遗传算法的改进支持向量机利用训练样本集的输入，遗传算法搜索核函数及其训练参数。基于遗传算法的整个优化过程如图2所示。

图2 遗传算法优化过程

遗传算法基于种群中适者生存原则，通过世代传递遗传信息来保留遗传信息。图3给出了所提出的支持向量机模型的框架。利用遗传算法在支持向量机中寻找两个参数的最佳组合，得到较小的分类误差。

图3 GA-SVM的流程图

2 实验数据

研究使用的4个数据集来自BCI Competition 2008 data sets 2b数据，由奥地利格拉茨理工大学提供。数据集由9名受试者的左右手运动图像数据组成。每个运动图像数据由三个电极(C3、Cz和C4)提取，采样频率为250 Hz的双极记录法进行记录[15]。滤波器采用0.5～100 Hz的带通滤波器和50 Hz的陷波滤波器。在前两个实验中，每个实验都从固定十字架和听觉刺激开始。3 s后，一个可视箭头指示左或右运动想象任务，此命令持续1.25 s，然后进行MI任务。MI持续4 s后，受试者有至少1.5 s的短暂休息时间。图4所示为数据集实验时间线。

图4 实验进程

采用Ag/AgCl电极记录脑电图如图5所示。

图5 采样电极位置

3 实验结果

3.1 小波分解结果

离散小波变换(DWT)通过将信号分解为一个粗糙的近似值，并从时域信号的连续高通滤波和低通滤波中获得的详细信息。对于许多信号，低频成分相当重要，它常常蕴含着信号的特征，而高频成分则给出信号的细节或差别[16]。根据事件相关同步/去同步原理，进行想象运动会在相关的频率段产生能量的变动，因此，α节律、β节律的脑电就会被定位。对于EEG信号在30 Hz以上可以采用低通滤波器滤除。在本研究中，使用小波变换来分解C3、C4和Cz通道滑动窗内的每个脑电片段。矩形窗口的大小是相邻窗口重叠的1/2。利用Daubechies 4-tap小波将每段脑电信号分解为4个层次，分解后得到16个不同的频率段。

频段分解后，子带小波系数所对应的频率范围如表 1 所示。

表1 每个子带频带分布

获得低频子带CAL，高频子带CDL，CDL-1，…,CD1。结果存储在低通信道{A4}和高通信道{D1、D2、D3和D4}中。

子带D4、D3在脑电信号频段，则选择这两个子带信号进行信号重构。如图6所示小波分解后的D3子带信号的频段集中分部在8～16 Hz，如图7所示D4则分布在16～30 Hz。

图6 D3小波分解子带信号

图7 D4小波分解子带信号

3.2 EMD分解结果

EMD应用筛选过程将信号分解为IMFs和一个残基的余数，可以看作是一个单调的斜率，一个只有一个极值的函数，或者最后一个IMF。考虑到信号的主观性质、数据长度和极值数目等因素，IMF的数目将不是恒定的。IMF表示简单的振荡模式；而高阶IMF表示慢振荡模式，低阶IMF表示快振荡模式。所有脑电信号可以覆盖的IMF数量为8，如图8所示为D3子带分解的前8个IMF分量的波形图[图8(a)]和频谱图[图8(b)]，图9所示为D4子带分解的前8个IMF分量的波形图[图9(a)]和频谱图[图9(b)]。

f为频率

Pearson系数用来衡量两个独立的服从正态分布的连续变量之间的相关性。表2为每个IMF的脑电信号的排列熵特征的Pearson系数。两个特征之间的Pearson系数绝对值越大，表示特征间冗余越大，两个特征越不应该被同时选中；反之，说明两个特征间的相关性越小[16]。表2中绿色底纹所标注的Pearson系数绝对值都小于0.04，说明具有特征极弱相关或者无相关。因此，选择IMF2、IMF3、IMF7、IMF5的特征进行分类。

表2 每个排列熵特征的Pearson系数

3.3 GA-SVM分类

使用GA算法对支持向量机参数的优化，可以使SVM更好对不同的特征进行分类[18]。

针对支持向量机存在两个模型参数c和g，运用遗传算法选择最优的c和g对支持向量机进行优化。调用MATLAB遗传算法工具箱，实现逐步启发式优化，经反复试验确定最大的进化代数maxgen=100，种群最大数为sizepop=40(种群规模)，交叉概率为pc=0.9，变异概率为pm=0.01，并且设置支持向量机进行10折交叉验证，得到最优c=0.000 913 62，g=0.634 15。

如图10为使用GA算法优化前后分类准确率对比变化图，可得未用GA算法的最高准确率可到99.22%，最低78.56%，平均分类准确率97.30%。使用GA算法之后分类准确类最高100.00%，最低96.70%，平均达97.64%。推出GA算法优化了g和c使准确率更加平稳，减少了因个体差异性误判的可能性[19]。

图10 GA对排列熵参数优化前后准确率对比

3.4 实验总结

为了突出本文方法的优越性，将其与其他最新的研究进行了比较。表3给出了本文方法和其他方法的分类结果的比较。本文方法在使用较少通道数得到比其他研究方法更优的分类结果。Pearson系数筛选特征分类结果进行对比，特征筛选后平均准确率到达97.30%。GA优化SVM和特征筛选都起作用，有效地提高了MI-EEG的分类精度[9]。

表3 与2008 BCI竞赛结果比较

4 结论

本文的总体目标是提出一种应用于脑电信号特征分类的方法。为了去除不太重要的数据信息，降低计算复杂度，特征提取方法只使用优势IMF。利用小波变换提取脑电运动想象窄带信号，再经EMD提取窄带信号的IMF分量，计算IMF的排列熵作为分类特征。对IMF各个排列熵特征的Pearson系数进行比较，选取较小的IMF的排列熵作为特征值。然后用支持向量机对其进行分类，运用遗传算法对支持向量机参数进行优化，实现了对脑电信号较高精度的分类。在进一步的研究中，将研究不同的特征选择和优化方法，并与GA-SVM方法进行比较，并将GA-SVM方法应用于BCI系统中。