基于弹性变分模态分解的癫痫脑电信号分类方法*

景鹏 张学军2)† 孙知信

1) (南京邮电大学电子与光学工程学院,微电子学院,南京 210023)

2) (南京邮电大学,射频集成与微组装技术国家地方联合工程实验室,南京 210023)

3) (南京邮电大学,江苏省邮政大数据技术与应用工程研究中心,南京 210003)

4) (南京邮电大学,国家邮政局邮政行业技术研发中心(物联网技术),南京 210003))

癫痫脑电信号分类对于癫痫诊治具有重要意义.为了实现病灶性与非病灶性癫痫脑电信号的分类,本文利用弹性网回归重构变分模态分解算法,提出弹性变分模态分解算法并将其应用到所提癫痫脑电信号分类方法中.该方法先将原信号分割成多个子信号,并对各子信号进行弹性变分模态分解,然后从分解后的不同变分模态函数中提取精细复合多尺度散布熵作为特征,最后利用支持向量机进行分类.针对癫痫脑电的公共数据集,最终的实验结果表明,准确率、灵敏度和特异度三个性能指标分别达到92.54%,93.22%和91.86%.

1 引 言

癫痫是一种由于人脑神经系统紊乱造成脑功能障碍的疾病,据世界卫生组织2019年统计,全世界有超5000万人患有癫痫[1],对具有抗药性的癫痫患者需要进行侵入性手术治疗,切除大脑中的部分致痫区域,因此手术前准确识别大脑中的致痫区域至关重要.脑电图 (electroencephalograph,EEG)记录了人脑皮层发出的轻微电流信号,癫痫发作时人脑中枢神经系统会出现异常使得同步神经元突然放电,相应的脑电图就会呈现出异常波动,从人脑致痫区域获得的脑电信号称为病灶性脑电信号,从非致痫区域获得的脑电信号称为非病灶性脑电信号[2],因此对病灶性脑电信号和非病灶性脑电信号进行分类,有助致痫区域的识别.

目前国内外对于癫痫脑电信号的分类技术主要由特征提取与特征分类组成,其本质可归结为模式识别问题[3].特征提取先要选择合适的特征,如高阶矩[4]、概率密度函数[5]、自回归模型系数[6]、各类熵[7-10]等均已成功应用到癫痫脑电信号的分类中.文献[11]提出了一种新的复杂性度量特征,即精细复合多尺度散布熵(refined composite multiscale dispersion entropy,RCMDE),其作为特征在处理非线性不平稳生物信号时效果较好.此外,特征提取通常在时域、频域或时频域内进行,其中又以时频域分析法应用较多,如短时傅里叶变换[12]、小波变换[13]、经验模态分解等,经验模态分解作为一种自适应时频处理方法应用较广,但存在端点效应、模态混叠等缺陷[14].由文献[14]提出的变分模态分解 (variational mode decomposition,VMD)是一种新的基于完全非递归思想的时频域信号分析方法,其主要思想是将信号分解成若干个围绕在各个中心频率附近的窄带变分模态分量且中心频率不断变化,与经验模态分解通过不断循环筛分获取本征模态函数不同,VMD通过寻找约束变分模型的最优解自适应获取变分模态函数(variational mode function,VMF).VMD 在构建初始是针对盲信号,使用Tikhonov正则化思想去构建一个约束变分方程[14],但本质是通过对最小二乘回归加残差平方和二次惩罚项,从而达到收缩估计系数的目的,这也可以理解为使用了线性规划中的岭回归(ridge regression)来构建约束方程,虽然岭回归与最小二乘回归相比降低了估计值的方差,提高了估计精度,但是它的回归结果中包含所有的预测变量,没有进行变量选择,因此会影响模型的准确性[15].而套索回归(LASSO regression)使用一次惩罚项来对变量进行收缩,相较于岭回归收缩程度要小,能选出更精确的模型,但是如果预测变量具有群组效应或者预测变量相关性较强时,则会导致估计不稳定[15].Zou和Hastie[16]提出的弹性网回归(elastic net regression)方法综合了岭回归和套索回归的思想,兼有套索回归和岭回归的优点.

本文利用弹性网回归对变分模态分解算法进行改进,构建一种新的时频域信号分析方法,称为弹性变分模态分解 (elastic variational mode decomposition,EVMD),并通过在 EVMD 域中提取精细复合多尺度散布熵作为癫痫脑电信号的特征,然后利用支持向量机 (support vector machine,SVM)实现病灶性脑电信号和非病灶性脑电信号的分类.此外,将本文所提方法应用到实际癫痫脑电数据集中进行仿真分析,以证明该方法的有效性.

2 方 法

应用本文所提方法处理癫痫脑电信号的步骤如图1所示.先利用弹性变分模态分解算法对分割后的癫痫脑电子信号进行分解,从分解后的各VMF中提取RCMDE作为特征,最后利用支持向量机对特征进行分类,进而完成对病灶性脑电信号与非病灶性脑电信号的区分.本节将对所提方法中使用到的主要算法进行叙述.

2.1 弹性变分模态分解

图1 弹性变分模态分解处理癫痫脑电信号的流程图Fig.1.Architecture of processing epileptic EEG by EVMD.

本文提出的弹性变分模态分解是一种完全非递归式的时频域信号分析与处理方法,与原始变分模态分解算法不同的是,EVMD利用弹性网回归代替VMD中构建约束变分方程时使用的岭回归,即EVMD算法构建的约束变分方程中既存在L2正则化项,也存在L1正则化项.本节将对EVMD算法的构造与求解进行详细叙述.

首先,通过以下步骤建立一个约束变分模型：1)通过希尔伯特变换求出每个模态函数的解析信号,进而可以获得信号的单边频谱； 2)通过一个混合指数项将每个模态函数调制到对应中心频率的基频带上； 3)通过弹性网回归方法施加L1正则化与L2正则化来估计带宽,即求信号梯度二范数的平方与梯度一范数的和.建立的约束变分模型为

其中,x为原脑电信号,{uk}={u1,u2,···,uK}为所有模态的集合,{ωk}={ω1,ω2,···,ωK}为对应中心频率的集合,为冲激函数,∂t为对t求偏导.

接着对约束变分模型求最优解.先将拉格朗日乘法算子以及损失项加入到约束变分模型中,从而将约束变分模型转化为一个非约束变分模型,得到的增广拉格朗日函数为

其中λ为拉格朗日乘法算子,α为二次惩罚因子,β为一次惩罚因子,为损失项.

然后利用交替方向乘子法来求解(2)式的最优解,迭代公式为：

迭代停止条件为：

其中,η为拉格朗日乘子更新参数,ε为收敛容限.

下面介绍模态uk与中心频率ωk的详细更新求解步骤,即迭代式(3)式与(4)式的具体推导过程.

A模态uk更新

步骤1将(3)式改写为如下等值公式：

步骤2将(7)式进行整体傅里叶变换：

步骤 3将 (8)式中ω+ωk替换为ω,即将带有惩罚因子的项中的ω用ω-ωk代替,得：

步骤4因为实信号具有埃米尔特对称性,因此将(9)式改为非负频域上的半空间积分：

步骤5将(10)式看为一个二次优化问题进行求解,解得：

B中心频率ωk更新

步骤1因为中心频率ωk不出现在重构保真项中,因此可将(4)式重新写为

步骤2与求解模态uk的操作类似,将(12)式进行傅里叶变换到频域,并最终变换成非负频域上的半空间积分：

步骤3求解(13)式二次优化问题,可得：

至此,完整的EVMD算法已经介绍完毕,其具体流程可以参考算法1所示.

Algorithm 1：EVMD

2.2 精细复合多尺度散布熵

精细复合多尺度散布熵是文献[11]在散布熵的基础上对信号进行多尺度量化得到的,避免了散布熵因单一尺度上处理信号会出现复杂性特征提取不完全等问题.文献[11]详细介绍了RCMDE的原理,这里仅简述其对脑电信号处理时的计算步骤：

步骤1假设脑电信号x经弹性变分模态分解后得到的信号uk是长度为Ψ的时间序列,则uk的第h个粗粒度近似信号为

其中,γ={1,2,···,N},N=Ψ/τ,τ为尺度因子.

步骤2使用正态累积分布函数将粗粒度近似信号a映射到b={b1,b2,···,bN}上：

其中,µ表示均值,σ表示粗粒度近似a的标准差.b的范围从0到1.

步骤3为了实现信号的多尺度化,将b以线性变换的形式映射到{1,2,··,c}中,记为z,即

其中R(*)是取整运算,c代表类别个数.

步骤4若嵌入维数标记为m和时间延迟标记为d,则时间序列定义为

其中i={1,2,···,N-(m-1)d}.

步骤5设每一个时间序列对应一个散布模 式Θv0v1···vm-1,其 中v= 1,2,···,c.若=v0,对应的散布模式为Θv0v1···vm-1.

步骤6计算每个散布模式Θv0v1···vm-1的概率

步骤7根据香农熵的定义,信号uk的精细复合多尺度散布熵为

3 实验结果与分析

3.1 实验数据集

本文使用的癫痫脑电信号数据集是来自瑞士伯尔尼大学神经科提供的一个公共数据集[17].该数据集数据采自5位患有颞叶病灶癫痫的病人,共包括3750对病灶性数据段和3750对非病灶性数据段,每对数据段均包含两列数据,分别采自同一个区域相邻的两个通道.数据段的采样频率为512 Hz,时间为 20 s[18].本文选用该数据集中全部脑电数据段作为实验数据,即选用3750对病灶性脑电数据段与3750对非病灶性脑电数据段,每对数据段均包含两列数据.

3.2 实验

3.2.1 EEG 数据段分割

先将原始脑电数据经过0.5—100 Hz的带通滤波器进行预处理,从而去除部分伪迹与噪声.然后将每对时长为20 s的数据段分割为时长均为1 s的子数据段,且相邻两个子数据段之间的时间重叠率为50%,因此每个数据段共分割为39个子数据段.之所以需要对原始脑电数据进行分割,第一是因为EVMD算法与VMD一样在处理长时间非平稳非线性的脑电信号时模态的谱带会随时间剧烈变化影响分解效果,而将信号分割为短时间的子信号,可以一定程度规避这种局限性.第二是因为可以减小脑电数据中随机波动点的影响,比如假设某数据段中由于噪声等原因出现数个非预期的随机波动点,若对整个数据段进行分析显然这些点会影响最终分析结果,而对若干子数据段进行分析,这些点只会影响部分子数据段的分析结果,而对最终所求结果的影响将降低.

3.2.2 EVMD 分解

因为测试时间和被试的不同,各频带段的脑电信号的中心频率会随着神经活动而发生轻微的变化,简单地固定中心频率的带通滤波无法消除这种变化的影响,而EVMD旨在将信号分解为围绕在中心频率附近的变分模态分量,即变分模态函数VMF,且其中心频率是不断迭代变化的,故可以捕捉到这种变化.

根据2.1节对EVMD的介绍,可以看出EVMD算法在对信号进行分解之前需要分别设定： 分解模态的个数K,二次惩罚因子α,一次惩罚因子β三个参数.对于惩罚因子α,β,其值过大会引起模态重叠,较小会引入噪声,本文根据经验建议将α设定为 2000,β设定为20.对于分解模态的个数K,若K值过大会造成过分解,产生无用分量,同时增加计算复杂度,若K值过小会造成欠分解,使部分带限信号分解不出来造成原信号信息的丢失.为寻求合适的K值,本文先随机选择一对病灶性脑电数据段,并分别进行3重、4重、5重变分模态分解,即将K值分别设为 3,4,5.图2所示为不同K值下,各个变分模态函数分量其中心频率ω随迭代次数的变化曲线.可以看出,当K= 3 时,变分模态函数VMF1—VMF3的中心频率均在40 Hz以下.当K= 4 时,变分模态函数 VMF1—VMF4 的中心频率仍在 40 Hz以下.但当K= 5 时,可以看到VMF5的中心频率超过了50 Hz.一方面由于脑电活动的信息大都包含在低频带(频率 ＜ 40 Hz),另一方面为了防止信号的过分解或欠分解,本文选取K= 4.

图2 不同K值下中心频率随迭代次数的变化曲线Fig.2.The curves of the center frequency with the number of iterations under different K values.

3.2.3 RCMDE 提取特征

根据(23)式,计算RCMDE需要分别设置：嵌入维度m,类别个数c,时间延迟d以及尺度因子τ四个参数.类别个数c若过大会导致具有较大差异的两个量被归为同一类,过小则导致具有较小差异的两个量被归为不同类,根据文献[11]的建议,将c值设为6.对于嵌入维度m,为了保证统计可靠性,文献 [12]中建议cm＜Ψ,其中Ψ为数据长度,本文将采样率 512 Hz,采样时间 20 s的数据段分割为时长为 1 s的子数据段,因此Ψ= 512,而 64= 1296 ＞Ψ,故将m值设为 3.时间延迟d为正整数,但d＞ 1 会导致模态混叠,故取 1.尺度因子τ决定信号粗粒化程度,为了选取合适的τ值,先将τ最大值设为15.

选用数据集中全部3750对病灶性脑电信号和3750对非病灶性脑电信号,对分割后的子数据段进行4重EVMD分解,从分解后得到的4个VMF(VMF1—VMF4)中分别提取15个尺度因子下的RCMDE,其特征熵值的均值加减标准差随尺度因子变化曲线如图3所示.图3中不同尺度因子下的RCMDE均值用曲线相连,每个均值点上下两个点即表示加减标准差.从图3可以看出,从非病灶性脑电信号提取的RCMDE 特征熵值基本比病灶性脑电信号提取的RCMDE熵值大,这说明人脑非致痫区域的神经活动相对致痫区域更活跃,也说明非致痫区域提取的脑电信号相较致痫区域提取的脑电信号更随机、更不平稳,这一结果也与文献[18]的研究结果一致.此外,之所以会出现从不同VMF中提取的两类信号RCMDE特征熵值的均值差异性不同的现象,是因为不同VMF是体现不同中心频率上的脑电信号数据,而人脑非致痫区域与致痫区域在不同频带段上的神经活动会不同,比如经4重EVMD分解下VMF3代表的是中心频率为20 Hz附近频带段的信号,此时两类脑电信号RCMDE熵值的均值的差值与其他VMF相比较小,表明在20 Hz附近人脑非致痫区域与致痫区域的神经活动更相近,产生的脑电信号的差异性较其他频带段较小.但是随着尺度因子的增大,总体非病灶性脑电信号与病灶性脑电信号提取的RCMDE特征熵值的均值的差值会有一个先增大再减小的过程,尤其是VMF2尺度因子等于15时,VMF3尺度因子大于12时,会出现病灶性脑电信号提取的RCMDE特征熵值的均值更大的反常现象,这是因为信号被过度粗粒化从而出现失真,说明尺度因子不能过大.本文从计算量的角度综合考虑将尺度因子τ设为7,即每个数据段提取7个尺度上的RCMDE值.

图3 从各 VMF 中提取的 RCMDE 熵值的均值 (± 标准差)随尺度因子变化曲线Fig.3.The curve of mean value(± SD) of RCMDE computed from VMF.

为了更好评估从非病灶性脑电信号提取的RCMDE特征与从病灶性脑电信号提取的RCMDE特征的区分度,这里使用学生检验的p值来评估其统计学上的差异性,其结果如表1所列.从表2可以看出数据经4重EVMD分解后得到的4个VMF中提取的RCMDE特征p值均小于0.05,因此在统计学上具有显著差异,说明RCMDE可以作为区别病灶性和非病灶性癫痫脑电数据的分类特征.

表1 从各 VMF 中提取的 RCMDE 特征 p 值Table 1.The p values of RCMDE computed from VMF.

表2 EVMD 与 VMD 实验结果对比Table 2.Comparison of experimental result between EVMD and VMD.

3.3 分类结果

选用数据集中全部3750对病灶性癫痫脑电数据段和3750对非病灶性癫痫脑电数据段,按上述方法进行处理,并将得到的每段脑电数据的RCMDE特征送入SVM进行特征分类,SVM选用线性核函数,通过网格搜索法将惩罚参数设定为0.52,并重复进行10次5折交叉验证实验,采用准确度、灵敏度和特异度这三个指标对最终分类结果进行度量,结果如图4所示.由图4可知10次5折交叉验证实验的平均准确度、灵敏度和特异度分别可达 92.54%,93.22%,91.86%.

为了比较本文提出的EVMD算法与原始VMD在处理病灶性癫痫脑电信号与非病灶性信号分类中的性能,将本文实验中的EVMD换成原始VMD,其余步骤不变,重复进行实验,对比实验的结果如表2所列.可以看出,EVMD算法相较原始VMD算法,三个性能指标均有3个百分点以上的提高,因此按本文所提方法进行病灶性癫痫脑电信号与非病灶性癫痫脑电信号的分类实验中,本文所提EVMD算法相较原始VMD算法性能具有一定优势.

图4 10 次 5 折交叉验证实验结果折线图Fig.4.The line chart of the results by 5-fold cross validation for 10 times.

对相同的实际癫痫脑电信号数据集,将本文方法得到的结果与其他文献的结果相比较,结果如表3所列,可以看出本文所提方法性能优于其他方法.

表3 本文方法与其他方法对比Table 3.Comparison between proposed method and previously published methods.

4 结 论

本文利用弹性网回归重构了变分模态分解的约束方程,提出一种新的时频域信号分析与处理算法,称为弹性变分模态分解,并给出了具体推导过程.提出一种基于弹性变分模态分解,从分解后的VMF中提取RCMDE特征,并利用支持向量机对其进行分类的病灶性癫痫脑电信号与非病灶性癫痫脑电信号的分类方法.应用本文提出的方法对公共癫痫脑电数据集进行处理,得到最终准确率为92.54%,灵敏度 93.22%,特异度 91.86%,并与针对同一公共数据集的其他处理方法进行比较,证明了该方法的有效性.然而,弹性变分模态分解的参数(分解模态的个数K、惩罚因子α与β)只能通过经验设定或者根据处理信号不断试错设定,无法自动获取最优参数,因为弹性变分模态分解与原始变分模态分解结构类似,所以现在针对变分模态分解参数的优化方法,如相关性分析优化[21]、粒子群算法优化[22,23]等方法,理论上应该也可以应用到本文提出的弹性变分模态算法上,因此,对弹性变分模态分解的参数优化将是后续研究的一个可能方向.