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

全文HTML

--> --> -->

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.弹性变分模态分解

本文提出的弹性变分模态分解是一种完全非递归式的时频域信号分析与处理方法, 与原始变分模态分解算法不同的是, EVMD利用弹性网回归代替VMD中构建约束变分方程时使用的岭回归, 即EVMD算法构建的约束变分方程中既存在L2正则化项, 也存在L1正则化项. 本节将对EVMD算法的构造与求解进行详细叙述.
首先, 通过以下步骤建立一个约束变分模型: 1)通过希尔伯特变换求出每个模态函数的解析信号, 进而可以获得信号的单边频谱; 2)通过一个混合指数项将每个模态函数调制到对应中心频率的基频带上; 3)通过弹性网回归方法施加L1正则化与L2正则化来估计带宽, 即求信号梯度二范数的平方与梯度一范数的和. 建立的约束变分模型为

$\begin{split}&\mathop {\min }\limits_{\{ {u_k}\} ,\{ {\omega _k}\} } \left\{ \sum\limits_k {\left[ {\left\| {{\partial _t}\left[ {\left({\delta (t) + \frac{{\rm{j}}}{{{\rm{\pi }}t}}} \right) \times {u_k}(t)} \right]{{\rm{e}}^{ - {\rm{j}}{\omega _k}t}}} \right\|_2^2} \right.} \right.\\&\left.{{ + }}\left. {{{\left\| {{\partial _t}\left[ {\left({\delta (t) + \frac{{\rm{j}}}{{{\rm{\pi }}t}}} \right) \times {u_k}(t)} \right]{{\rm{e}}^{ - {\rm{j}}{\omega _k}t}}} \right\|}_1}} \right] \right\}\\& \text{s.t.}\sum\limits_k {{u_k}(t) = x},\\[-20pt]\end{split}$

其中, x为原脑电信号, $\{ {u_k}\} = \{ {u_1}, {u_2}, \cdots, {u_K}\}$

为所有模态的集合, $\{ {\omega _k}\} = \{ {\omega _1}, {\omega _2}, \cdots, {\omega _K}\}$

为对应中心频率的集合, $\displaystyle \sum\nolimits_k = \sum\nolimits_{k = 1}^K {}$

, $\delta (t)$

为冲激函数, ${\partial _t}$

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

$\begin{split}&\varGamma (\{ {u_k}\} ,\{ {\omega _k}\} ,\lambda) \\=\;& \alpha \sum\limits_k {\left\| {{\partial _t}\left[ {\left({\delta (t) + \frac{{\rm{j}}}{{{\rm{\pi }}t}} \times {u_k}(t)} \right.} \right]{{\rm{e}}^{ - {\rm{j}}{\omega _k}t}}} \right\|_2^2} \\&+ \beta \sum\limits_k {\left\| {{\partial _t}\left[ {\left({\delta (t) + \frac{{\rm{j}}}{{{\rm{\pi }}t}}} \right) \times {u_k}(t)} \right]{{\rm{e}}^{ - {\rm{j}}{\omega _k}t}}} \right\|_1^{}} \\&||x(t) - \sum\limits_k {{u_k}(t)} ||_2^2 + [\lambda ,x(t) - \sum\limits_k {{u_k}(t)} ],\end{split}$

其中$\lambda $

为拉格朗日乘法算子, $\alpha $

为二次惩罚因子, $\beta $

为一次惩罚因子, $||x(t) - \sum\limits_k {{u_k}(t)} ||_2^2$

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

$u_k^{n + 1} \leftarrow \mathop {\arg \min }\limits_{{u_k}} \varGamma (\{ u_{i < k}^{n + 1}\},\{ u_{i \geqslant k}^n\},\{ \omega _i^n\},{\lambda ^n}),$

$\omega _k^{n + 1} \leftarrow \mathop {\arg \min }\limits_{{\omega _k}} \varGamma (\{ u_i^{n + 1}\},\{ \omega _{i < k}^{n + 1}\},\{ \omega _{i \geqslant k}^n\},{\lambda ^n}),$

${\lambda ^{n + 1}} \leftarrow {\lambda ^n} + \eta \times \left(x - \sum\limits_k {u_k^{n + 1}}\right).$

迭代停止条件为:

$\sum\limits_k {||u_k^{n + 1} - u_k^n||_2^2/} ||u_k^n||_2^2 < \varepsilon ,$

其中, $\eta $

为拉格朗日乘子更新参数, $\varepsilon $

为收敛容限.
下面介绍模态${u_k}$

与中心频率${\omega _k}$

的详细更新求解步骤, 即迭代式(3)式与(4)式的具体推导过程.
A模态${u_k}$

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

$\begin{split}u_k^{n + 1} =\;& \mathop {\arg \min }\limits_{{u_k}}\! \Bigg\{ \!\alpha \left\| {{\partial _t}\!\left[\! {\left({\delta (t) \!+\! \frac{{\rm{j}}}{{{\rm{\pi }}t}}}\! \right)\! \times\! {u_k}(t)} \right]\!{{\rm{e}}^{ - {\rm{j}}{\omega _k}t}}} \right\|_2^2 \\&{{ + }}\beta \left\| {{\partial _t}\left[ {\left({\delta (t) + \frac{{\rm{j}}}{{{\rm{\pi }}t}}} \right) \times {u_k}(t)} \right]{{\rm{e}}^{ - {\rm{j}}{\omega _k}t}}} \right\|_1^{}\\& {{{ + }}\left\| {x(t) - \sum\limits_i {{u_i}(t) + \frac{{\lambda (t)}}{2}} } \right\|_2^2} \Bigg\};\\[-20pt]\end{split}$

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

$\begin{split}\hat u_k^{n + 1} =\;& \mathop {\arg \min }\limits_{{{\hat u}_k}} {{\Bigg\{ }}\alpha ||{\rm{j}}\omega [(1 \!+\! {\mathop{\rm{sgn}}} (\omega \!+\! {\omega _k})){{\hat u}_k}(\omega \!+\! {\omega _k})]{{||}}_2^2\\&{{ + }}\beta ||{\rm{j}}\omega [(1 + {\mathop{\rm{sgn}}} (\omega + {\omega _k})){{\hat u}_k}(\omega + {\omega _k})]{{||}}_1^{}\\& {{{ + }}\left\| {\hat x(\omega) - \sum\limits_i {{{\hat u}_i}(\omega) + \frac{{\hat \lambda (\omega)}}{2}} } \right\|_2^2} \Bigg\};\\[-20pt]\end{split}$

步骤3　将(8)式中$\omega + {\omega _k}$

替换为$\omega $

, 即将带有惩罚因子的项中的$\omega $

用$\omega - {\omega _k}$

代替, 得:

$\begin{split}\hat u_k^{n + 1} =\;& \mathop {\arg \min }\limits_{{{\hat u}_k}} {{\Bigg\{ }}\alpha ||{\rm{j}}(\omega - {\omega _k})[(1 + {\mathop{\rm{sgn}}} (\omega)){{\hat u}_k}(\omega)]{{||}}_2^2\\&{{ + }}\beta ||{\rm{j}}(\omega - {\omega _k})[(1 + {\mathop{\rm{sgn}}} (\omega)){{\hat u}_k}(\omega)]{{||}}_1^{}\\&{{ + }} {\left\| {\hat x(\omega) - \sum\limits_i {{{\hat u}_i}(\omega) + \frac{{\hat \lambda (\omega)}}{2}} } \right\|_2^2} \Bigg\};\\[-20pt]\end{split}$

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

$\begin{split}\hat u_k^{n + 1} =\;& \mathop {\arg \min }\limits_{{{\hat u}_k}} \int_0^{ + \infty } 2 \alpha {(\omega - {\omega _k})^2}|{\hat u_k}(\omega){{|}}_{}^2\\&{{ + }}\beta {\rm{j}}(\omega - {\omega _k}){\hat u_k}(\omega)\\&{{ + }}\left| {\hat x(\omega) - \sum\limits_i {{{\hat u}_i}(\omega) + \frac{{\hat \lambda (\omega)}}{2}} } \right|_{}^2{\rm{d}}\omega ;\end{split}$

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

$\hat u_k^{n + 1}(\omega) = \frac{{\hat x(\omega) - \sum\limits_{i \ne k} {{{\hat u}_i}(\omega) + \frac{{\hat \lambda (\omega)}}{2}} }}{{1 + 2\alpha {{(\omega - \omega {}_k)}^2} + \beta {\rm{j}}(\omega - {\omega _k})}}. $

B中心频率${\omega _k}$

更新
步骤1　因为中心频率${\omega _k}$

不出现在重构保真项中, 因此可将(4)式重新写为

$\begin{split}\omega _k^{n + 1} =\;& \mathop {\arg \min }\limits_{{\omega _k}}\! \left\{\! \alpha \left\| {{\partial _t}\!\left[\! {\left({\delta (t) \!+\! \frac{{\rm{j}}}{{{\rm{\pi }}t}}} \!\right) \!\times \!{u_k}(t)} \right]{{\rm{e}}^{ - {\rm{j}}{\omega _k}t}}} \right\|_2^2\right.\\&\left.{{ + }}\beta \left\| {{\partial _t}\left[ {\left({\delta (t) + \frac{{\rm{j}}}{{{\rm{\pi }}t}}} \right) \times {u_k}(t)} \right]{{\rm{e}}^{ - {\rm{j}}{\omega _k}t}}} \right\|_1^{} \right\};\end{split}$

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

$\begin{split}\omega _k^{n + 1} = \;&\mathop {\arg \min }\limits_{{\omega _k}} {{\bigg\{ }}\int_0^{ + \infty } 2 \alpha {(\omega - {\omega _k})^2}|{\hat u_k}(\omega){{|}}_{}^2\\&{{ + }}\beta {\rm{j}}(\omega - {\omega _k}){\hat u_k}(\omega){\rm{d}}\omega\bigg\} ;\end{split}$

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

$\omega _k^{n + 1} = \dfrac{{\displaystyle\int_0^{ + \infty } {\omega |{{\hat u}_k}(\omega){|^2} + \frac{{\beta {\rm{j}}}}{{2\alpha }}|{{\hat u}_k}(\omega)|{\rm{d}}\omega } }}{{\displaystyle\int_0^{ + \infty } | {{\hat u}_k}(\omega){|^2}{\rm{d}}\omega }}.$

至此, 完整的EVMD算法已经介绍完毕, 其具体流程可以参考算法1所示.
Algorithm 1: EVMD
　Initialize $\hat u_k^1$

, $\omega _k^1$

, ${\hat \lambda ^1}$

, n ← 0
repeat
　　n ← n + 1
　for k = 1: K do

$\begin{split}&{\rm{Update}}\; {\hat u_k} : \hat u_k^{n + 1}(\omega) \leftarrow \\&\frac{{\hat x(\omega) - \displaystyle\sum\limits_{i < k} {\hat u_i^{n + 1}(\omega) - \displaystyle\sum\limits_{i > k} {\hat u_i^{n + 1}(\omega)} + \frac{{{{\hat \lambda }^n}(\omega)}}{2}} }}{{1 + 2\alpha {{(\omega - \omega _k^n)}^2} + \beta {\rm{j}}(\omega - \omega _k^n)}}\end{split}$

$\begin{split}&{\rm{Update}}\;{\omega _k} : \;\omega _k^{n + 1} \leftarrow \\&\frac{{\displaystyle\int_0^{ + \infty } {\omega |\hat u_k^{n + 1}(\omega){|^2} + \frac{{\beta {\rm{j}}}}{{2\alpha }}|\hat u_k^{n + 1}(\omega)|{\rm{d}}\omega } }}{{\displaystyle\int_0^{ + \infty } | \hat u_k^{n + 1}(\omega){|^2}{\rm{d}}\omega }}　\qquad\end{split}$

　end for
　for all $\omega \geqslant 0$

${\hat \lambda ^{n + 1}}(\omega) \leftarrow {\hat \lambda ^n}(\omega) + \eta \times \left({\hat x(\omega) - \displaystyle\sum\limits_k {\hat u_k^{n + 1}} (\omega)} \right)$

　end for
until convergence: $\displaystyle\sum\limits_k \dfrac{||u_k^{n + 1} \!-\! u_k^n||_2^2}{ ||u_k^n||_2^2} \! < \! \varepsilon$

2.2.精细复合多尺度散布熵

-->

2.2.精细复合多尺度散布熵

精细复合多尺度散布熵是文献[11]在散布熵的基础上对信号进行多尺度量化得到的, 避免了散布熵因单一尺度上处理信号会出现复杂性特征提取不完全等问题. 文献[11]详细介绍了RCMDE的原理, 这里仅简述其对脑电信号处理时的计算步骤:
步骤1　假设脑电信号x经弹性变分模态分解后得到的信号$u_k^{}$

是长度为$\varPsi $

的时间序列, 则$u_k^{}$

的第h个粗粒度近似信号为

$a_{h,\gamma }^\tau = \frac{1}{\tau }\sum\limits_{\chi = h + (\gamma - 1)\tau }^{h + \gamma \tau - 1} {{u_{k,\chi }}} ,$

其中, $\gamma = {{\{ 1, 2, }}\cdots{{, }}N{{\}, }}\;N{{ = }} {\varPsi }/{\tau }$

, $\tau $

为尺度因子.
步骤2　使用正态累积分布函数将粗粒度近似信号a映射到$b = \{ {b_1}, {b_2}, \cdots, {b_N}\}$

上:

$b_\gamma ^\tau = \frac{1}{{\sigma \sqrt {2{\rm{\pi }}} }}\int_{ - \infty }^{a_\gamma ^\tau } {{{\rm{e}}^{\tfrac{{ - {{(t - \mu)}^2}}}{{2{\sigma ^2}}}}}} {\rm{d}}t,$

其中, μ表示均值, σ表示粗粒度近似a的标准差. b的范围从0到1.
步骤3　为了实现信号的多尺度化, 将b以线性变换的形式映射到{1, 2, ···, c}中, 记为z, 即

$z_\gamma ^c = R(c \cdot {b_\gamma } + 0.5),$

其中R(*)是取整运算, c代表类别个数.
步骤4　若嵌入维数标记为m和时间延迟标记为d, 则时间序列$z_i^{m, c}$

定义为

$z_i^{m,c} = \{ z_i^c,z_{i + d}^c,\cdots,z_{i + (m - 1)d}^c\} ,$

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

步骤5　设每一个时间序列$z_i^{m, c}$

对应一个散布模式${\varTheta _{{v_0}{v_1}\cdots{v_{{\rm{m}} - 1}}}}$

, 其中$v{{ = 1, 2, }}\cdots, c$

. 若$z_i^c = {v_0}$

, $z_{i + d}^c = {v_1}$

, $z_{i + (m - 1)d}^c = {v_{m - 1}}$

, 则$z_i^{m, c}$

对应的散布模式为${\varTheta _{{v_0}{v_1}\cdots{v_{{\rm{m}} - 1}}}}$

.
步骤6　计算每个散布模式${\varTheta _{{v_0}{v_1}\cdots{v_{{\rm{m}} - 1}}}}$

的概率

$p({\varTheta _{{v_0}{v_1}\cdots{v_{{\rm{m}} - 1}}}}) = \frac{{\# ({\varTheta _{{v_0}{v_1}\cdots{v_{{\rm{m}} - 1}}}})}}{{N - (m - 1)d}},$

其中$\# ({\varTheta _{{v_0}{v_1}\cdots{v_{{\rm{m}} - 1}}}})$

为$z_i^{m, c}$

对应散布模式的个数.
步骤7　根据香农熵的定义, 信号$u_k^{}$

的精细复合多尺度散布熵为

$ \begin{split} &{\rm{RCMDE}}(u_k^{},m,c,d,\tau) \\ =\;& - \displaystyle \sum\limits_{\varTheta = 1}^{{c^m}} {\bar p({\varTheta _{{v_0}{v_1}\cdots{v_{m - 1}}}})} \ln \bar p({\varTheta _{{v_0}{v_1}\cdots{v_{m - 1}}}}), \end{split} $

其中, $\bar p({\varTheta _{{v_0}{v_1}\cdots{v_{m - 1}}}}) = \dfrac{1}{\tau }\displaystyle \sum\limits_1^\tau {p_h^\tau }$

表示粗粒度近似信号a的散布模式$\varTheta $

的平均值.

3.实验结果与分析

3.1.实验数据集

-->

3.1.实验数据集

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

3.2.实验

-->

3.2.实验

3.2.1.EEG数据段分割

-->

3.2.1.EEG数据段分割

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

3.2.2.EVMD分解

-->

3.2.2.EVMD分解

因为测试时间和被试的不同, 各频带段的脑电信号的中心频率会随着神经活动而发生轻微的变化, 简单地固定中心频率的带通滤波无法消除这种变化的影响, 而EVMD旨在将信号分解为围绕在中心频率附近的变分模态分量, 即变分模态函数VMF, 且其中心频率是不断迭代变化的, 故可以捕捉到这种变化.
根据2.1节对EVMD的介绍, 可以看出EVMD算法在对信号进行分解之前需要分别设定: 分解模态的个数K, 二次惩罚因子$\alpha $

, 一次惩罚因子$\beta $

三个参数. 对于惩罚因子$\alpha $

, $\beta $

, 其值过大会引起模态重叠, 较小会引入噪声, 本文根据经验建议将$\alpha $

设定为2000, $\beta $

设定为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值下中心频率随迭代次数的变化曲线
Figure2. The curves of the center frequency with the number of iterations under different K values.

3

3.2.3.RCMDE提取特征

-->

3.2.3.RCMDE提取特征

根据(23)式, 计算RCMDE需要分别设置: 嵌入维度m, 类别个数c, 时间延迟d以及尺度因子τ四个参数. 类别个数c若过大会导致具有较大差异的两个量被归为同一类, 过小则导致具有较小差异的两个量被归为不同类, 根据文献[11]的建议, 将c值设为6. 对于嵌入维度m, 为了保证统计可靠性, 文献[12]中建议c^m < $\varPsi $

, 其中$\varPsi $

为数据长度, 本文将采样率512 Hz, 采样时间20 s的数据段分割为时长为1 s的子数据段, 因此$\varPsi $

= 512, 而6⁴ = 1296 > $\varPsi $

, 故将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熵值的均值(± 标准差)随尺度因子变化曲线
Figure3. 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可以作为区别病灶性和非病灶性癫痫脑电数据的分类特征.

	VMF1	VMF2	VMF3	VMF4
p值	1.17 × 10^–4	2.38 × 10^–2	4.42 × 10^–2	7.50 × 10^–3

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

指标	准确度/%	灵敏度/%	特异度/%
EMVD	92.54	93.22	91.86
VMD	89.49	87.56	88.12

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

2

3.3.分类结果

-->

3.3.分类结果

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

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

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

	提取方法	特征	分类器	灵敏度、特异度、准确度/%
文献[7]	EMD	平均香农熵、样本熵等	LSSVM	90, 84, 87
文献[8]	TQWT	模糊熵	LSSVM	83.86, 85.46, 84.67
文献[19]	MFDFA	多重分形	SVM	92.5, 92.7, 92.2
文献[20]	EWT	RPS, CTM	LSSVM	88, 92, 90
本文	EVMD	RCMDE	SVM	93.22, 91.86, 92.54

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

本站小编 Free考研考试/2021-12-29

English Abstract

eEpileptic electroencephalogram signal classification method based on elastic variational mode decomposition

Corresponding author:Zhang Xue-Jun, xjzhang@njupt.edu.cn

全文HTML

2.1.弹性变分模态分解

2.2.精细复合多尺度散布熵

3.1.实验数据集

3.2.实验

3.2.1.EEG数据段分割

3.2.2.EVMD分解

3.2.3.RCMDE提取特征

3.3.分类结果

相关话题/数据 信号 实验 文献 优化

领限时大额优惠券,享本站正版考研考试资料!

基于Zernike模型系数优化的椭球型窗口光学系统像差校正

铷原子系综自旋噪声谱实验研究

桌面飞秒极紫外光原子超快动力学实验装置

椭球胶体在圆球胶体体系中扩散行为的实验研究

电池材料数据库的发展与应用

石榴石型固态电解质表界面问题及优化的研究进展

基于PE型压机中子衍射高温高压组装的优化设计与实验验证

基于时谱信号分析的在轨空间目标姿态感知

退火效应增强铁磁异质结太赫兹发射实验及机理

X波段高重频长脉冲高功率多注相对论速调管放大器的设计与实验研究