基于EEMD与DSS-ApEn的脑电信号消噪方法*

孟 明,杨国雨,高云园,甘海涛,罗志增

(杭州电子科技大学智能控制与机器人研究所,杭州 310018)

头皮脑电EEG(Electroencephalogram)由于其易采集、无创性以及较好的时间分辨率,在最近被广泛关注的脑机接口应用中,有着不可替代的作用。但由于头皮脑电是一种非平稳非线性极其微弱的随机信号[1],极易被心电,眼电,电磁干扰,工频干扰等大量外界干扰信号淹没[2],所以消噪就成了脑电信号预处理阶段中最重要的步骤之一。

文献[3]提出的经验模态分解EMD(Empirical Mode Decomposition)是一种自适应的数据分解方法,能够较好地处理随机非平稳信号,具有良好的自适应性、信号的完备性等[4]。因此,EMD成为现代信号处理分析领域的一个研究热点,然而EMD方法存在模态混叠现象:当信号的极值点分布不均匀时,会导致一个相似尺度信号分布在不同的内蕴模态函数IMF(Intrinsic Mode Function)分量中或者一个IMF分量中出现多个尺度信号,从而使IMF分量失去了原本的物理意义。

Wu和Huang在EMD基础上改进提出了集合经验模态分解[5]EEMD(Ensemble Empirical Mode Decomposition),通过加入白噪声改变信号的极值点分布,进而有效地改善EMD的模态混叠现象,极大提高了信号的信噪比。但是基于EEMD的信号消噪方法通过直接去除EEMD分解结果中的前几个高频IMF分量来消噪,导致在去掉高频IMF分量中噪声的同时去除了其中的有效信息成分[6]。

降噪源分离DSS(De-noising Source Separation)是盲信号分离广泛应用的一种技术。Valpola等人在研究基于贝叶斯准则的快速独立分量分析算法的过程中,提出了DSS的基本理论概念[7]。DSS本质是在源信号复杂结构未知的的条件下,根据信号的统计特性,选择合适的降噪函数将复杂信号分解为若干组成分量,从而实现源信息的挖掘。文献[8]利用基于空间滤波的DSS方法分析处理脑电信号,有效地降低了刺激诱发反应范式中的噪声。王元生等[9]把DSS应用于机械工程领域,有效地提取出转子振动故障信号特征。但DSS得到的各独立源分量的顺序不确定,文献[8]通过人工观察区分独立源分量,文献[9]通过计算源信号和消噪信号的相关系数来区分。以上区分方法存在消噪过程中不能直接确定独立源分量是噪声信号还是有效信号的缺陷。

本文在使用近似熵ApEn(Approximate Entropy)区分DSS结果中各独立源信号是有效信号还是噪声信号的基础上,提出一种基于EEMD与DSS-ApEn相结合的脑电信号消噪方法,并与独立EEMD消噪方法以及基于EEMD与改进提升小波的消噪方法进行比较。

1 原理和算法

1.1 EEMD算法

EEMD主要利用了白噪声在各个频段具有能量一致的统计特性来解决EMD的模态混叠现象[10],其主要原理是:在不连续的信号中,低频分量的极值点间隔分布稀疏、而高频分量的则分布密集。当在待处理信号中加入白噪声,信号的低频分量的极值点分布发生改变,整个频带中的极值点间隔分布均匀,保证每个固有模态函数时域的连续性,从而避免了模态混淆。EEMD分解的具体步骤[11]为:

步骤1 往待处理信号s(t)中添加均值为0、标准差为常数的白噪声n1(t),得到待分解的信号s1(t)。

步骤2 EMD算法分解信号s1(t),得到n个IMF分量imfi1(t)和余项rn1(t),即:

(1)

式中,n为分解层数。

步骤3 重复以上两个步骤T-1次,但每次都加入同分布的不同白噪声。

步骤4 对相应的IMF分量进行总体平均运算处理,以消除白噪声对IMF分量的影响,最终得到各IMF分量为:

(2)

式中,imfi(t)为总体平均后的第i个IMF分量;T为总体平均次数;imfij(t)为第j次加入白噪声后获得的第i个IMF分量。

1.2 DSS算法原理

1.2.1 盲源分离基本模型

设n个信号源S={s1(t),s2(t),…,sn(t)}发出的信号在m个不同位置观测到的混合信号分别为X={x1(t),x2(t),…,xm(t)},即:

(3)

式中,aij为混合系数,ni为第i个源信号所对应的噪声,则可以将上式转换为矩阵形式:

X=AS+N

(4)

因此,源分离处理就是要在源信号S和混合系数矩阵A均未知的情况下,求解分离矩阵W,使得在观测信号中分离出对源信号的估计Y=[y1(t),y2(t),…,yn(t)],即:

Y=WX

(5)

1.2.2DSS理论框架

期望值最大化EM(Expectation-Maximization)算法是经典的盲源分离算法之一,其具体的算法可以分为两步:

步骤1 求源信号S期望值

q(S)=p(S|A,X)=p(X|A,S)p(S)/p(X|A)

(6)

式中,p(X|A)则是在已知混合矩阵情况下的条件概率,p(X|A,S)为根据前一次得到的源信号S与混合矩阵A对本次观测信号的概率进行估计,p(S|A,X)为源信号的后验概率,p(S)为前一次源信号S的先验概率。

步骤2 寻求期望值中的最大值

利用式(6)中求得的期望值,就可以根据EM算法的求得新的混合矩阵Anew,即:

(7)

式中,CXS和CSS是根据期望函数q(S)计算得到的协方差矩阵,即:

(8)

(9)

在EM算法的迭代求解过程中,所有的独立成分分量按照统计特性被一次性求解得到。相对EM算法中分量被同时求出,Hyvärinen等[12]在2001年提出利用预白化的方法逐次提取各独立成分源信号进行迭代求解,充分体现了组成成分的独立和重要性,这是提出DSS理论的重要理论基础。基本求解过程框架如下:

s=wTX

(10)

s+=f(s)

(11)

w+=Xs+T

(12)

(13)

式(10)通过似然估计函数计算源信号的噪声估计,式(11)中s+对应EM算法中对期望q(S)的估计,即为降噪过程。所有含有式(11)降噪环节的源分离方法被定义为DSS。式(12)为利用降噪处理后的源信号s+对分离矩阵w的重新估计。式(13)完成归一化。DSS对每个源信号的分离过程,就是在给定式(10)中的wT之后,通过不断迭代式(11)～式(13)求得每个源信号的分离信号。

DSS理论的关键在于根据待处理的信号特点,选择或者设计适合的降噪函数。文献[13]模拟实验分析的结果表明,基于正切降噪函数的源分离方法从非线性耦合信号中提取的分量与源信号时域相关系数最高,所以本文选取正切降噪函数f(s)=s-tanh(s)作为DSS的降噪函数。

1.3 近似熵

DSS结果中的各独立源分量可以看作各信号源,需要找出含有噪声的独立信号源。在信息论中,熵指数能够表征系统的复杂程度。近似熵[14]可以衡量时间序列复杂度,用于度量信号产生新模式的概率,越复杂的时间序列对应的近似熵越大,故本文采用近似熵作为脑电分量与噪声分量的判据。

近似熵的具体计算步骤为:

步骤1 对系统进行等间隔采样,得到时间序列u(1),u(2),…,u(N);

步骤2 计算时间序列u(i)的标准差SD;

步骤3 构造一组m维矢量:H(1),H(2),…,H(N-m+1),其中H(i)=[u(i),u(i+1),…,u(i+m-1)];

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

d[H(i),H(j)]=max[|u(i+k)-u(j+k)|],1≤i,j≤N-m+1,k=0,…,m-1

(14)

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

(15)

步骤6 定义φm(r):

(16)

步骤7 将比较序列长度参数m加1,即m=m+1,重复步骤3～7,得到φm+1(r);

步骤8 计算近似熵ApEn:

ApEn=φm(r)-φm+1(r)

(17)

ApEn的值显然与m,r的取值有关。Pincus[15]根据实践,建议取值m=2,r=0.1×SD。

1.4 基于EEMD与DSS-ApEn相结合的消噪方法

采用EEMD对信号进行分解,得到多个IMF分量,前几个分量主要包含信号中的高频成分,后面的分量主要包含信号中的低频成分,噪声强度随着IMF分量层次的增加也会越来越弱,其中,脑电信号的低频成分以有效信息为主,而高频成分含大量的噪声[16]。EEMD消噪方法通过直接去除EEMD分解中的前几个高频IMF分量来消噪,但是这样会导致消除高频中的噪声成分的同时亦剔除了其中的有效信息成分。

DSS能有效的提取各独立源分量,但是DSS结果中各独立源分量的顺序不确定,不能直接确定独立源分量对应的信号源是噪声信号还是脑电信号。由于噪声信号频率的复杂程度比脑电信号的要大很多,故采用近似熵来判断独立成分中的噪声信号。

针对以上情况,本文将EEMD与DSS-ApEn相结合,其具体消噪步骤如下。

步骤1 应用EEMD将待消噪信号分解为IMF集,去掉最高频IMF分量得到新的IMF集Q;

步骤2 对步骤1中新的IMF集Q应用DSS算法,得出一组独立源分量S,并求出各独立源分量的频谱;

步骤3 计算步骤2中频谱的近似熵,选择近似熵最大的独立源分量作为噪声信号滤除,得到一组新的独立源分量S′;

步骤4 将混合矩阵A乘以独立源分量S′变换出新的IMF集,并将该IMF集累加重构,得到消噪后的脑电信号。

2 实验结构与分析

2.1 仿真实验

(18)

标准信号和加入高斯白噪声后得到的信号如图1 所示,加噪信号的信噪比为1 dB,图2是EEMD分解加噪信号后的IMF分量。

按照1.4小节中的步骤2进行处理,去掉最高频IMF分量,对剩余的IMF应用DSS得到各独立源分量如图3所示,并求出各独立源分量频谱如图4所示,各独立源分量频谱近似熵和时域近似熵如表1 所示。

图1 标准信号和信噪比为1 dB时加噪标准信号

图2 EEMD分解加噪信号后的IMF分量

图3 DSS各独立源分量

图4 DSS各独立源分量频谱图

表1 DSS各独立分量频谱近似熵和时域近似熵

从图4中可以看出IC6是噪声信号,但是观察表1中各独立源分量的频谱近似熵和时域近似熵,频谱近似熵可以很好的区分出噪声信号,而时域近似熵的区分度不是很大,故本文采用频谱近似熵区别各独立源信号。表1中IC6的频谱近似熵最大,故第6个独立源分量是噪声信号,滤除噪声信号后得到一组新的独立源分量。按照1.4小节中的步骤4进行重构,即可得出消噪后的信号。本文消噪方法与独立EEMD分解后直接舍弃高频IMF分量的方法、基于EEMD与改进提升小波的消噪方法[18]的消噪效果如图5所示。

使用信噪比SNR(Signal to Noise Ratio)和均方根误差RMSE(Mean Squared Error)作为评价加噪信号去噪效果指标,分别定义为:

(19)

(20)

表2 3种消噪方法加了不同信噪比噪声的仿真信号的实验数据

图8 原始C3通道EEG分解后IMF分量时域图

2.2 EEG消噪实验

本文进一步采用实际采集的脑电信号数据进行消噪效果验证。应用美国NeuroScan公司的scan4.3脑电采集系统作为脑电信号的采集设备,采集C3、C4和Cz 3个电极处的脑电信号,采样频率为250 Hz,精度为32 bit。受试者为身体健康的在校研究生。脑电电极按照国际标准导联10-20系统放置,以左侧乳突为参考电极,右侧乳突为接地电极。

本实验通过采集想象左手、右手所产生的脑电信号进行消噪处理。每次实验持续9 s。0～2 s时,受试者保持身体放松状态;2 s时,系统发出“滴”提示音并且在显示屏中间显示一个十字光标,提示受试者集中注意力;3 s～6 s时,受试者根据屏幕上随机出现的左(←)、右(→)两个方向的诱导图片分别想象左手挥动、右手挥动这两种运动,直至第6 s时系统发出“嘟”的提示音,受试者结束想象。一次完整实验的采样时序如图6所示。

实验选择C3通道运动想象部分的脑电信号作为研究对象,其波形如图7所示。

图6 实验设计流程时序图

图7 C3通道运动想象部分脑电信号

图7中,横坐标为信号的采样点个数,纵坐标为脑电信号的幅值。在原始C3通道EEG经过EEMD分解后,得到8个IMF分量和一个余量res,所得IMF分量的时域图如图8所示。从图8可以看出,各阶IMF分量包含不同的时间特征尺度,且随着阶次的增加,频率越小。

按照本文1.4小节中的消噪步骤进行消噪。本文消噪方法与独立EEMD分解后直接舍弃高频IMF分量的方法、基于EEMD与改进提升小波的消噪方法的消噪时域图如图9所示、频域图如图10所示。

图9 3种消噪方法时域对比图

图10 3种消噪方法频谱对比图

结合图9和图10所示,虽然经前两者算法处理后的信号噪声得到了明显的消除,但是EEMD分解消噪方法得到的信号毛刺较多,且30 Hz以上的高频部分还有保留,而基于EEMD与改进提升小波消噪方法消噪后所得的信号过于平滑,频率在20 Hz～30 Hz细节信号滤除过多,发生了信号失真,部分细节信号被当做噪声消除了。而采用基于EEMD与DSS-ApEn相结合的消噪方法消噪后的脑电信号波形相对清晰,更重要的是原始信号的细节特征也被很好地保留下来。这是由于EEMD分解后仅对噪声主导的高频IMF分量进行消噪,从而有效保留了有用的脑电信号细节信息。DSS在迭代求解独立源分量的过程中,将利用前一次的求解结果作为后一次求解的先验知识,具有很强的自适应性,所以DSS可以更好地从高频IMF分量中分离出细节信息,达到更好保留细节信息的目的。

3 结论

本文提出一种基于EEMD与DSS-ApEn相结合的消噪方法。对EEMD分解得到的IMF分量处理方法进行改进,将滤除最高频分量后的IMF作为DSS的输入。对DSS得到的各独立源分量,选择频谱近似熵最大的独立源分量作为噪声信号滤除。采用不同消噪方法的仿真和实际脑电信号消噪实验结果表明,本文的消噪方法具有更强的噪声抑制能力。良好的消噪效果和保留信号细节的能力,为后续脑电信号的特征提取与模式识别打下坚实的基础。