去除癫痫脑电信号运动伪迹的变分模态分解-自适应熵阈值方法

张力行,张四聪,徐光华,2,李焕发,吴永程

(1.西安交通大学机械工程学院,712000,西安;2.西安交通大学机械制造系统工程国家重点实验室,710049,西安;3.西安交通大学第一附属医院,710061,西安)

癫痫是大脑神经元突发性异常放电导致短暂大脑功能障碍的慢性疾病。目前国内患病人数已达900万左右,且每年有约40万的新增患者[1],癫痫已成为神经科的第二大疾病。虽然癫痫可以使用药物和手术进行治疗,但是仍有30%的顽固性癫痫病人无法被治愈,存在着癫痫随时可能发作,造成患者受伤、残疾甚至死亡的风险[2]。由此,迫切需要对癫痫进行提前预测,有效避免癫痫突然发作带来的伤害。

癫痫作为大脑中枢功能障碍疾病,其诊断与预测主要依托脑电信号进行。早期由于癫痫病人脑电信号采集设备体积大,主要采用病人卧床状态下脑电信号对癫痫进行诊断[3],由于脑电信号作为微伏级微弱信号,滤除采集信号中混有的眼电、心电等伪迹成为脑电信号去噪的主要研究方向[4-6]。近年来,便携式脑电设备的出现,设备更小巧、支持无线信号传输、功率较低且使用的电极数较少[7],解决了正常生活人群的脑电采集问题。国内外多家公司推出的便携式脑电设备已被有效使用,例如Emotiv公司开发的EPOC+无线脑电仪[8],采用盐水电极、支持14个通道的脑电信号采集,并提供了手机端和电脑端的应用软件开发平台,为基于脑电分析的情绪识别、娱乐游戏和医疗康复等领域新技术发展奠定了基础,也使得顽固性癫痫病人脑电预测技术发展成为可能。便携设备采集的癫痫脑电信号中,由于患者活动时电缆摇摆或电极与头皮相对移动,又增加了脑电信号的运动伪迹,其功率谱覆盖低频到高频相对较宽的范围[9],对准确进行癫痫的脑电信号预测产生了较大的阻碍。因此,研究有效去除癫痫脑电信号运动伪迹的消噪方法十分必要。

目前,对脑电图(EEG)中运动伪迹去除主要采用盲信号分离技术,常用的方法有独立成分分析(ICA)[10-11]、典型成分分析(CCA)[12]、经验模式分解(EMD)[13]和模板匹配法[14-16]。文献[11]中使用Informax ICA方法进行盲源分离,通过观察被试者在4种行走速度下脑电地形图、自相关函数和谱图来手动识别和去除运动伪迹;文献[12]认为运动伪迹具有更低的自相关性,使用CCA方法得到多组自相关性逐渐降低的向量基,然后将识别出的运动伪迹分量置0,逆变换恢复去伪迹后的信号。

由于ICA、CCA等方法只能用于多通道脑电信号的分析且要求信号非高斯分布,不适合单通道信号预处理。文献[17]采用EMD结合ICA的方法,有效地将ICA方法去伪迹的范围从多通道扩展到单通道,但是由于EMD方法存在较严重的频谱混叠的现象(模态混叠),会导致有用脑电信号较严重地流失。文献[13]结合集成经验模式分解(EEMD)、CCA和独立向量分析(IVA)方法,能够较好地分离出肌肉伪影,但是该方法基于固定阈值,无法确定是否足以在大型受试者群体中自动去除肌肉伪影。

此外,文献[14]中提出一种基于步态伪迹模板的自适应滤波器,实现对步态伪迹的去除。虽然基于步态伪迹模板的方法可以较好地识别出步态伪迹成分,但是对于一些简单运动,例如扭头等运动伪迹无法有效地去除,具有一定的局限性。

本文提出一种基于变分模态分解-自适应熵阈值法(VMD-AET)的运动伪迹识别方法。利用VMD分解有效克服EMD方法存在频率混叠导致有用EEG信号流失的问题,可以较好地实现频谱分离;基于自适应阈值的选择实现更多类型的运动伪迹去除;并使用交大一附院的癫痫数据集进行有效性的验证。将本文方法与ICA-阈值法、EMD-阈值法的去伪迹效果进行比较,验证了本文方法的有效性。

1 VMD-AET方法去除运动伪迹

1.1 VMD方法原理

VMD方法是一种自适应、完全非递归的模态变分和信号处理方法,具有可以确定模态分解个数的优点[18]。VMD分解在搜索和求解过程中可以自适应地匹配每种模态的最佳中心频率和有限带宽,实现信号的频域划分、分量的有效分离,最终获得变分问题的最优解。VMD方法不仅实现了提升复杂度高和非线性强的时间序列的平稳性,从而获得包含多个不同频率尺度且相对平稳的子序列的目的,而且克服了EMD方法存在端点效应和模态分量混叠的问题[19]。VMD方法的核心思想是构建和求解变分问题[20],其求解过程如下。

首先构造变分问题。假设原始信号f被分解为k个分量,分解序列为具有中心频率的有限带宽的模态分量。求各模态的估计带宽之和,使其最小。约束条件为所有模态之和等于原始信号,相应表达式为

(1)

求解式(1),引入Lagrange乘法算子λ,将约束变分问题转变为非约束变分问题,得到增广Lagrange表达式为

L({uk},{ωk},λ)=

(2)

式中:α为二次惩罚因子,作用是降低高斯噪声的干扰。利用交替方向乘子(ADMM)迭代算法结合帕塞瓦尔定理、傅里叶等距变换,优化各模态分量和中心频率,并搜寻增广Lagrange函数的鞍点。交替寻优迭代后的uk、ωk和λ的表达式如下

(3)

(4)

(5)

VMD方法在求解过程中因为最小化分量频带宽度而克服了频带混叠的现象,基于单通道求解的原理优于ICA方法,只适用于正定问题。综上所述,VMD方法在频带分离效果和使用条件上满足去除脑电信号运动伪迹的需求。

1.2 VMD-AET方法

非线性动力学方法[21]已经被证实在反映脑电信号的复杂性特征方面具有显著优势而备受关注,其中样本熵[22]、近似熵、模糊熵[23]和能量熵[24]常作为脑电信号去伪迹特征。由于运动伪迹与静态分量在能量分布上的差异性,因此本文使用能量熵进行求解。能量熵的计算如下

对脑电信号f(t)进行VMD分解

(6)

式中:uk(t)为第k个VMF分量;K为VMD分量的个数。

VMF分量的能量为

(7)

VMF能量谱占总能量谱的比重为

(8)

每个VMF分量的能量熵为

Hk=-PklgPk

(9)

由能量熵定义可知,各频带能量占总能量越多时能量熵越小。3种状态下的时域信号如图1所示,观察静坐状态和运动状态下的EEG可以发现,运动伪迹分量的幅值远大于正常状态下的脑电幅值,因此推知运动伪迹具有更高的能量占比和更小的能量熵,正常状态下的脑电信号则相反。

图1 3种状态下的时域信号Fig.1 Time-domain signal diagram under three states

进一步观察图1可以发现,静坐状态的EEG信号变化速度缓慢,上下点头时变化更快,跑步时变化最快。对运动状态下脑电信号进行傅里叶变换,并对扭头频率和步频进行计算,可以发现,频谱图的特征频率与计算得到的运动本身的频率相近。由于实际运动时癫痫患者可能按照任意自发的速度进行扭头、走路或跑步,因此导致了脑电信号中运动伪迹分量的特征频率分布较广。由于VMD方法具有较好的频带分离功能,可以很好地将脑电信号中运动伪迹分量分离出来,因此本文结合VMD和能量熵提出一种VMD-AET运动伪迹去除方法。

VMD-AET方法的详细流程如图2所示。首先,对含有运动伪迹的EEG信号进行50 Hz的递归滤波器陷波滤波和0～64 Hz的巴特沃斯带通滤波(一般认为与癫痫有关的脑电信号频率在此范围),去除工频干扰、提取有用的脑电频带信号;在综合考虑计算复杂度和算法性能的基础上,设置VMD的参数k=6以进行6层的变分模式分解来得到VMF分量;然后,分别求VMF分量的能量熵并进行降序排列;将最小到最大的能量熵作为阈值,将能量熵小于阈值的VMF分量作为运动伪迹分量进行去除;比较不同阈值情况下的伪迹去除率和信噪比提升指标,可以发现,当阈值设置为最大能量熵时,去运动伪迹效果最好;最后,将剩下的正常脑电分量进行VMD逆变换,得到不含运动伪迹的EEG。

图2 VMD-AET方法流程图Fig.2 Flow chart of VMD-AET algorithm

1.3 信号评价指标

为进一步精确地评价本文方法和其他方法去除运动伪迹的效果,引入伪迹去除率和信噪比提升两个评价指标[25]对信号的去噪效果进行评价。前者可以较好地反映去噪后的信号恢复到去噪前形态的能力,后者反映了信噪比的提升能力。下面详细介绍这两项指标。

伪迹去除率为

(10)

式中:Cref为参考信号的自相关;Crec为去伪迹后的信号与参考信号的互相关;Cart为参考信号与伪迹信号的互相关。伪迹去除率反映了去噪后的信号与参考信号的相似程度,它越大说明去噪后与去噪前信号形态越相似。

信噪比提升为

(11)

式中:Verbf和Veraf分别为去噪前后信号的方差。信噪比提升越大说明去噪效果越好。

2 实验验证与分析

为了验证本文方法的有效性,首先设计在实验室采集5名健康被试者不同运动状态下的脑电数据,以验证该方法去除运动伪迹的有效性,进而对西安交通大学第一附属医院的5名癫痫患者发作前活动癫痫脑电数据进行分析,以验证该方法的实际应用能力。

2.1 实验室数据采集与分析

2.1.1 数据采集 实验采用奥地利Guger Technologies公司研发的12导有线脑电采集设备,以1 200 Hz的频率进行连续采样,滤波范围为0.01～100 Hz。脑电电极的排布方式依据国际标准10/20系统,以左耳为参考电极。为了解不同运动方式和速度的脑电信号的特征,本文选择5名被试者(4男1女,年龄在21～24岁之间),让其进行静坐、上下点头、左右扭头、0.5 m/s行走(走路1)、1.2 m/s行走(走路2)、1.4 m/s行走(走路3)、1.6 m/s行走(走路4)、1.8 m/s行走(走路5)、2 m/s跑步,并在这9种状态下采集每名被试者在静坐和不同运动状态下的脑电信号。为了消除不同运动状态之间的影响,要求被试者在每一种运动状态测试结束后进行一定时间的休息,监测其心率恢复到正常后再进行下一运动状态的测试。各种运动状态说明如下。

状态1 静坐。让被试者在椅子上静坐。

状态2 上下点头。让被试者以正常速度上下点头。

状态3 左右点头。让被试者以正常速度左右扭头。

状态4 慢走(走路1)。让被试者以0.5 m/s的速度在跑步机上慢走。

状态5 正常走(走路2)。让被试者以1.2 m/s的速度在跑步机上行走。

状态6 正常走(走路3)。让被试者以1.4 m/s的速度在跑步机上行走。

状态7 正常走(走路4)。让被试者以1.6 m/s的速度在跑步机上行走。

状态8 快走(走路5)。让被试者以1.8 m/s的速度在跑步机上快走。

状态9 跑步。让被试者以2 m/s的速度在跑步机上跑步。

图3为实验采集到的一名健康被试者的脑电数据。由图3可以看出,在静坐状态下脑电信号的幅值最小且波动较小;在扭头阶段,幅值和波动的变化速度增加;在走路和跑步阶段,脑电信号的幅值随着运动速度的增加变化越来越大,在跑步状态下信号的幅值最大且变化速度最快。

图3 一被试者在静坐和8种运动状态下的脑电信号Fig.3 EEG signals of a subject in sedentary and eight motion states

对EEG信号进行傅里叶变换,图4为左右扭头状态下信号的频谱。经过步频和速率的计算公式可知被试者的实际扭头频率大约为0.5 Hz,与特征频率基本一致,计算其他运动状态下也有类似的结论。

图4 一名被试者在左右扭头状态下的信号频谱Fig.4 Spectrum of a subject in the state of left and right head-turning

2.1.2 数据处理与分析 运动伪迹去除的具体步骤如下。

步骤1对信号进行滤波。本文选择了二阶递归滤波器陷波滤波去除50 Hz工频干扰,以及巴特沃斯滤波器进行低通滤波,频率设置范围为0～64 Hz。

经滤波处理可得到5名被试者的EEG信号。图5为一名被试者滤波前后的时域,可以看出,滤波前EEG信号中存在的50 Hz工频干扰和100 Hz的整流信号在滤波后均被有效地消除。由此说明陷波滤波和巴特沃斯滤波器均达到了有效的滤波效果。

(a)原信号时域

(b)原信号频域

(c)滤波后信号时域

(b)滤波后信号频域图5 信号滤波前后的时频谱Fig.5 Time and frequency spectrums before and after filtering

步骤2选择VMD参数并进行VMD分解得到VMF分量。图6a为原信号和VMD分解后各个分量的时域信号。图6b为VMD分量的频谱,图6c为EMD分量的频谱,经过对比可以看出,VMD比EMD方法具有更好的频带分离能力。

(a)被试者的脑电原信号和VMD分解后的时域

(b)VMD分解后各分量的频谱

(c)EMD分解后信号的频谱图6 VMD分解前后时域图与VMD、EMD频谱Fig.6 Time domain diagram before and after VMD decomposition, and VMD and EMD spectrums

步骤3求VMF分量的能量熵。图7为由每种运动状态下各VMF分量的能量熵绘制的折线图,阈值以下被认为是脑电信号运动伪迹分量。

(a)静坐和前4种运动状态

(b)静坐和后4种运动状态图7 9种状态下各个VMF分量的能量熵折线图以及阈值设置Fig.7 Energy entropy broken line diagram and threshold setting of each VMF component in nine states

步骤4对能量熵进行降序排列。分别将最小到最大的能量熵作为阈值,并将能量熵小于阈值的VMF分量作为运动伪迹进行去除。

步骤5比较不同阈值情况下所有被试者的伪迹去除率和信噪比提升。如图8所示,为设置的5种阈值情况下得到的去伪迹信号的α、β的箱型图。可以看出,当阈值选择为最大的能量熵时,伪迹去除率和信噪比提升指标明显高于其他情况,说明此时去除脑电信号中运动伪迹的效果最好。

(a)不同阈值设置下α

(b)不同阈值设置下β图8 不同阈值设置下的脑电信号去运动伪迹评价指标α、β的箱型图Fig.8 Two evaluation indicators of motion artifact removal versus threshold

步骤6将去除运动伪迹后的脑电分量进行VMD逆变换,得到不含运动伪迹的EEG信号。

运用VMD-AET方法对实验室采集到的5名健康被试者不同运动状态下的脑电数据进行运动伪迹去除,并使用伪迹去除性能评价指标进行评价。本文将静坐时的脑电信号作为参考信号,对除静坐外的8种状态的运动伪迹去除效果进行分析。如图9所示,实验室采集到的健康被试者在每种运动状态下均能达到去除脑电信号中运动伪迹的效果,在被试者跑步时伪迹去除率和信噪比提升是最高的,达到5.345%和10.926 dB,两项指标的平均值分别为4.445%和7.556 dB,由此验证了本文方法对于不同运动状态下的脑电信号中的运动伪迹去除的有效性。

图9 健康被试者8种运动状态下的去伪迹结果及平均值Fig.9 Average values of artifact removal evaluation for healthy subjects in eight motion states

2.2 临床数据验证与分析

2.2.1 临床数据 为测试本文方法的临床应用效果,对采集到的西安交通大学第一附属医院5名癫痫患者的6次发作前活动癫痫脑电信号数据进行去运动伪迹效果测试和验证。临床头皮脑电数据的详细信息见表1。电极安放采用国际10/20系统连接方式,以双耳为参考电极,以每秒200 Hz的采样频率记录12导联脑电数据。以移动脑电设备(导线长度为1.5 m左右)进行数据采集,病人可以在该导线范围内进行小范围的活动,脑电信号中包含了扭头、慢走在内的运动伪迹信号。选取病人在卧床时的脑电信号作为参考信号。视频信息和脑电数据同步记录,以便由神经外科专家对发作开始和结束的时间进行标注,作为脑电信号分段的依据。

表1 临床头皮脑电数据的详细信息Table 1 Detailed information of clinical scalp EEG

2.2.2 临床验证与分析 将本文方法应用于临床采集的5名癫痫患者6次发作前的脑电数据上,得到伪迹去除率为5.54%,信噪比提升约为10.35 dB,说明本文方法适用于临床癫痫病人脑电信号的运动伪迹去除。

将本文方法与文献[10]的ICA-阈值法、文献[13]的EMD-阈值法效果进行对比的结果如图10所示。EMD-阈值法的伪迹去除率和信噪比提升最小,可以解释为EMD方法频带混叠较为严重导致有用EEG信号被去除较多、有效信息保留较少。ICA-阈值法的效果接近本文方法,说明ICA-阈值法去除脑电信号中运动伪迹效果与本文方法相近,但是因为ICA方法去伪迹时间较长,整体效果不如本文方法。

图10 本文方法与ICA-阈值法、EMD-阈值法去除癫痫病人脑电信号的运动伪迹的效果对比Fig.10 Comparison of the proposed method with ICA-threshold and EMD-threshold methods in motion artifacts removal from EEG signals of epileptic patients

综上所述,在临床癫痫病人发作前期脑电数据上,VMD-AET方法相比于其他方法具有更好的脑电信号中运动伪迹去除效果。

3 结 论

针对盲源分离法约束条件为正定问题和EMD方法频带划分能力差、步态模板提取法的适用范围小的问题,本文提出了一种基于VMD-AET的单通道运动伪迹去除方法,采用VMD方法对原始脑电信号进行分解,得到EEG的VMF分量信号;分别求取每个VMF的能量熵并将其设为阈值;通过对比不同阈值下的伪迹去除效果,自适应确定阈值大小取代以往的固定阈值方法。实验结果表明,当阈值取最大的情况下伪迹去除效果最好。将本文方法与常用的ICA-阈值法、EMD-阈值法对癫痫病人的脑电数据进行脑电信号运动伪迹去除效果进行比较,由信噪比提升指标最高可以得到本文方法的去伪迹能力较强,由伪迹去除率指标最高可以得出去除运动伪迹后的信号与静坐状态下的信号具有更高的相似性,由此说明信号具有较高的还原能力。