0 引言

癫痫是一种突然、反复发作的大脑功能障碍慢性病, 由于大脑异常电活动的起始部位及传播方式不同, 导致癫痫临床表现复杂多样, 包括短暂的感觉障碍、肢体抽搐、意识丧失、行为障碍等, 对患者身体、精神带来严重损伤^[1]. 脑电信号是脑神经元间信息交互产生的微弱的电信号, 能够反映大脑活动情况, 基于脑电信号的癫痫发作的预测、诊断和治疗, 其有效性和可行性已被证实^[2].

癫痫发作的预测是一个模式识别的过程, 需要数据前期的预处理, 特征提取, 以及训练分类模型识别发作前期、发作期、发作间期等不同状态. 深度学习已普遍地用于医学图像和生物电信号的处理, 在拥有大量数据的情况下深度学习的表现在很多方面超越了传统的特征提取加机器学习. 文献[3]采用CHB-MIT数据集, 将EEG信号分为发作间期、发作前期和发作期, 在时域和频域下分别绘制像素图, 再利用CNN学习特征后进行分类. 文献[4]将原始信号直接输入13层CNN预测癫痫发作. 文献[5]采用Freiburg医院数据集以及CHB-MIT数据集, 将EEG信号通过短时傅里叶转换为时频能量图, 通过CNN学习特征对癫痫发作状态进行分类. 文献[6]将多通道的癫痫EEG信号振幅谱分为19个频率子带特征, 使用CNN算法提取特征, 结合SVM进行分类. 文献[7]将长短期记忆网络(LSTM)引入使用EEG信号的癫痫发作预测中, 并利用一个两层LSTM网络来评价预测癫痫性能.

癫痫发作的准确预测, 对患者而言可以达到提前预警的作用, 有利于医生提前制定治疗计划. 但现有研究对多通道脑电信号分析较少, 癫痫发作前预测的准确率有待提高, 并且预测时间太短. 针对这些问题, 本文提出多尺度符号化排列传递熵, 分析多通道间的耦合关系, 生成同步矩阵图(SMD), 并筛选癫痫发作重要通道, 减少冗余特征对模型学习的干扰; 此外, 还利用功率谱密度能量图(PSDED)表征脑电信号特征, 与表示脑通道交互特征的同步矩阵图(SMD)进行融合, 结合脑电通道特征以及通道间交互的特征, 采用深度卷积神经网络模型对特征图进行训练, 对癫痫状态进行分类和预测, 分类及预测性能都得到了提升.

1 材料与方法

本文提出了多尺度符号化排列传递熵, 提取多通道癫痫EEG信号同步特征, 进行双向耦合性分析. 分析患者癫痫发作期以及发作间期的脑电信号, 生成同步矩阵图, 并通过显著性分析的方法分析两种状态的脑电信号, 筛选具有显著变化的通道. 对筛选通道后的脑电信号通过功率谱密度分析生成PSDED图表征脑电特征, 与描述脑通道交互特征的同步矩阵图(SMD)融合后用深度卷积神经网络进行训练, 再将训练后的模型应用于癫痫发作预测模型, 验证预测性能, 算法流程如图 1所示.

1.1 数据来源

数据取自美国麻省理工学院与波士顿儿童医院共同采集的癫痫数据集(CHB-MIT)^[8]. 患者均为停药后采集的数据, 对EEG进行长期监测, 记录EEG以及确定发作时刻. 其中24位患者的脑电数据, 总计时间约为980个小时, 且每个患者癫痫次数均大于3次. 数据集中的脑电信号为皮层脑电(EEG), 遵照“10-20标准”进行电极排布, 准确地涵盖了全头皮多通道, 更具全面性、可解释性. 数据说明中包含了癫痫发作的开始时间以及结束时间, 由专业的神经学医生标注. EEG的采样频率为256 Hz, 每例脑电数据中包含23个导联, 采用双极导联方法采集.

1.2 数据处理

原始EEG信号在大脑皮层上采集, 会引入人体内噪声干扰, 如眼电伪迹、肌电干扰等. 为了更好地提取EEG的特征, 需要对其进行消噪处理^[9]. 本文采用小波阈值滤波对EEG进行消噪. 使用‘db6’小波函数对EEG进行6层分解. 采用滑窗为4 s, 步长为2 s的方式选取患者的癫痫EEG进行分析. 本文将患者的EEG分为4类: 发作间期(interictal), 发作前30 min (preictal I), 发作前10 min(preictal II)以及发作期(seizure).

1.3 多尺度符号化排列传递熵

人脑是一个跨越多尺度的复杂系统, 脑电信号是一种跨越多尺度时空相互作用的非线性复杂信号, 具有较高的复杂度和噪声. 对癫痫脑电信号时空层次的同步特征进行分析, 有助于理解脑区间的同步关系, 本文采用多尺度排列传递熵对EEG通道间的信息交互进行分析.

时间序列在符号化时会捕捉大尺度的特征, 虽然损失了EEG部分细节特征, 但符号序列保存了更重要的动力学特征^[3]. 采用排列模式^[10-11]降低分析的复杂度, 减少计算量, 同时增强信号的抗噪能力和鲁棒性. 本文采用静态划分法对原始的时间序列进行粗粒化, 根据原始序列数值范围划分对应符号, 长度为$ n $的原始序列$ X $的符号化公式如下:

$ \begin{align} X^{s}=\begin{cases} S-1, x_{i}\geqslant X_{s-1};\\ \ \ \ \ \ \ \ \ \vdots \\ 1, X_{2} > x_{i}\geqslant X_{1};\\ 0, x_{i} <X_{1}. \end{cases} \end{align} $

(1)

其中: $ S $是符号的数量, $ S $的值越大对应的符号越多, 转换后的符号序列与原始序列越相似. 对$ X^s $进行相空间重构, 有

$ \begin{align} X^{s}_{t}=\{x^{S}_{t}, x^{S}_{t+\tau}, \ldots , x^{S}_{t+(m-1)\tau}\}. \end{align} $

(2)

其中: $ m $为嵌入维数, $ \tau $为延迟因子. 对时间重构序列按升序排列, 即

$ \begin{align} &x^{S}_{t+(j_{1}-1)\tau}\leqslant x^{S}_{t+(j_{2}-1)\tau} \leqslant\\ & x^{S}_{t+(j_{3}-1)\tau}\leqslant\ldots \leqslant x^{S}_{t+(j_{m}-1)\tau}, \end{align} $

(3)

会得到一组符号序列$ J=\{j_{1}, j_{2}, \ldots, j_{m}\} $为重构序列后对应元素原始位置的索引, 共有$ m $!种排列可能. 对$ m $!种排列对应$ \{0, 1, 2, \ldots, m!-1\} $进行符号化, 可以得到新的时间序列.

传递熵是反映两个序列双向信息流强度的方法^[12]. 通过信息熵反映序列间的因果关系, 表现为: 在时间序列$ Y $的历史时间尺度下, 时间序列$ X $的历史时间尺度包含的$ Y $的未来信息量^[13], 对于任意两个时间序列$ X $和$ Y $, 其概率分布函数分别是$ p(x) $和$ p(y) $，联合分布为$ p(x, y) $. 可定义$ X $对$ Y $的传递熵, 设时间序列$ X_{t} $和$ Y_{t} $的长度为$ n $, 各自从$ t-1 $到$ t-p $的历史记为$ X^{-}=X^{t-p}_{t-1}, Y^{-}=Y^{t-p}_{t-1} $. $ X_{t} $对$ Y_{t} $的传递熵$ T_{X\rightarrow Y} $定义为$ X $与$ Y^{-} $在$ X^{-} $下的互信息.

$ \begin{align} &T_{X\rightarrow Y} =I(X;Y^{-}|X^{-})=\\ &H(X|X^{-})-H(X|X^{-}, Y^{-})=\\ &\sum\limits_{x, y}^{} p(x, x^{t-p}_{t-1}, y^{t-p}_{t-1}){\rm log}\dfrac{p(x_{t}|x^{t-p}_{t-1}, y^{t-p}_{t-1})}{p(x_{t}|x^{t-p}_{t-1})}. \end{align} $

(4)

1.4 特征提取 1.4.1 同步矩阵图

虽然不同个体的癫痫发作EEG不同, 但他们在空间、频谱和时间特征上是相似的, 在Shoeb^[14]的研究中也有提及. 通过多尺度符号化排列传递熵生成的同步矩阵图, 可以很好地反映出脑电信号中多通道间的信息传递量以及交互特征, 体现脑区之间的耦合性变化. 本文通过多尺度排列传递熵得到的值生成$ n\times n $的同步矩阵, 表示$ n $个脑电通道间的耦合性强度, 同步矩阵为

$ \begin{align} X=\begin{bmatrix} T_{11} &T_{12} & T_{13} &\ldots & T_{1n}\\ T_{21} &T_{22} & T_{23} &\ldots & T_{2n}\\ T_{31} &T_{32} & T_{33} &\ldots & T_{3n}\\ \vdots &\vdots & \vdots &\ddots & \vdots\\ T_{n1} &T_{n2} & T_{n3} &\ldots & T_{nn} \end{bmatrix}, \end{align} $

(5)

其中$ T_{ab} $为通道$ a $到通道$ b $的有向耦合强度. 对角线为信号自身的耦合强度, 所以为0. 图 2为将同步矩阵归一化后的生成图. 可以发现, 患者在癫痫发作期和发作间期时部分通道间的耦合强度会有所变化.

1.4.2 功率谱密度能量图

针对脑电信号的频率能量特征, 对去噪后的多通道EEG采用功率谱密度(power spectral density, PSD)分析后提取能量变化的特征. PSD分析可以反映EEG的频率特性, 是常用的特征提取方法. 对于一个长度有限的非连续信号$ u(n)=0, 1, \ldots, N-1 $, 先进行离散傅里叶变换, 即

$ \begin{align} U(k)=\sum\limits_{n=0}^{N-1}u(n)W^{kn}_{N}. \end{align} $

(6)

其中: $ W_{N}={\rm e}^{-j2\pi / N} $, $ k=0, 1, \ldots, N-1 $. $ W_{N} $具有周期性, 所以$ U(k) $同样包含周期性, 二者周期大小为$ N $. 常用的PSD求取方法采用周期图法, 即

$ \begin{align} {\rm PSD}=U^{2}(k)/N. \end{align} $

(7)

其中: $ U(k) $为非连续信号的离散傅里叶变换, $ N $为快速傅里叶变换的采样点数, PSD为功率谱估计. 通过PSD分析可以体现出脑电信号的能量变化, 利用这一特性提出了癫痫脑电信号提取特征图方法^[15], 称为功率谱密度能量图法(power spectrum density energy diagrams, PSDED). 通过在EEG频域中的能量提取相应的特征, PSD可以表示信号功率在频域的分布. 共分4类脑电信号: 发作间期、发作前I (发作前30 min)、发作前II (发作前10 min)和发作期. 将$ n $个通道的EEG信号等分为32个频段, 并对不同频段的PSD函数进行积分, 最后组成一个$ n\times 32 $的二维矩阵, 将得到的二维矩阵进行归一化处理后生成能量图, 如图 3所示.

1.5 深度卷积神经网络

采用癫痫脑电信号分类方法(epileptic EEG signal classification, EESC)^[15], 利用深度卷积神经网络(DCNN)和从原始多通道EEG迁移学习的PSDED对4种不同的癫痫状态进行分类. 所提出的方法如图 4所示. 网络深度的增加伴随着大量的学习参数, 容易出现过拟合现象, 而且计算量会呈指数级增长. Inception网络降低了网络的复杂度, 核心内容是找到图片最佳的局部稀疏结构, 在同一层网络采用多个不同大小的卷积核, 增加网络的宽度, 改善网络的性能. Inception v3将$ n\times n $的卷积核分解为非对称的$ 1\times n $和$ n\times 1 $的卷积核, 可以有效减少网络参数个数. 当CNN网络层数达到一定数量后继续增加网络深度不仅不能提升分类效果, 而且会出现退化问题, 得到更差的分类性能. 残差网络(ResNet)可以有效解决网络过深导致的退化问题, 其中ResNet152是网络较深、分类性能出色的ResNet网络. 而Inception-ResNet-v2模型包含了ResNet和Inception网络的优点, 由3个模块组成, 分别是为了捕获更多图像特征的Stem模块, 可以增加网络层同时减少参数计算量的Inception-ResNet模块, 减少卷积核尺寸同时增加卷积核数量的Reduction模块. 该网络在解决退化问题的同时加快了网络训练的速度. 选取这3个非常具有代表性的DCNN: Inception-ResNet-v2、Inception-v3和ResNet152^[16], 在保持网络深度的同时有效地降低了训练复杂度, 并且解决了网络过深导致的退化问题. 在迁移学习框架中, 加载对应的ImageNet预训练权重, 在PSDED数据集上进行微调, 所得的3个模型视为3个特征提取器, 分别在PSDED数据集上提取瓶颈特征, 即输入全连接层之前的特征. 将3个DCNN模型提取得到的特征进行拼接后再输入后续的全连接层进行训练和降维, 最后再经过softmax进行多分类.

2 实验过程及结果 2.1 EEG耦合性分析及通道筛选

癫痫病人EEG相较于正常人EEG更为复杂, 分析时选取嵌入维数$ m=3 $, 一般选择延迟时间$ \tau=1 $进行分析^{[11, 17]}. 选择患者在癫痫发作间期和发作期的脑电进行多尺度分析. 选择合适的尺度研究EEG的耦合性极其重要, 患者在癫痫发作期和发作间期的传递熵随尺度的变化如图 5所示, 可以看出在尺度65左右时的传递熵值较大, 耦合性较强. 因此本文选取尺度为65进行实验.

为了优化通道选择, 利用EEG同步矩阵确定不同状态下变化明显的脑电通道. 本文采用单因素方差分析(ANOVA)检验^[18-19], 判断通道在对发作期与发作间期两种状态下是否具有显著性差异. 为保证通道筛选的一致性, 通道选择的策略全部为采用单样本方差进行分析, 并筛选通道. 当该通道显著性结果$ p<0.05 $时, 表明该通道在发作间期和发作期时变化明显, 保留该通道作为重要通道. 若结果$ p\geqslant0.05 $时, 则表明该通道不具有显著性差异, 将不再具体讨论该通道. 最后可以得到通道筛选后的同步矩阵, 如图 6所示.

使用符号化排列传递熵算法分析多通道癫痫EEG中各通道间两两耦合关系, 构造多通道EEG信号的同步矩阵, 通过显著性分析进行通道的筛选. 由于不同患者癫痫类型会有所区别, 癫痫发作源以及影响的通道也会不同, 所以需要对每个患者的EEG做具体的分析. 针对不同的患者采取不同的通道选取策略, 保证选取通道后不同状态提取的特征会有更加显著的差别. 保留显著性分析后剩余的通道, 根据不同患者共同的重要通道进行筛选, 最后选择12个通道. 由于T8-P8、P7-T7重复采集, 实际为10个通道. 表 1为筛选后的脑电通道.

2.2 癫痫发作预测结果

将通道筛选后的数据生成优化功率谱密度能量(optimization PSDED, OPSDED), 且特征图为$ 12\times 32 $, OPSDED表示EEG的能量特征. 结合多尺度符号化排列传递熵得出$ 12\times 12 $同步矩阵图(synchronization matrix diagram, SMD), SMD作为EEG通道的耦合关系表示的信息交互特征. 将两种特征图通过矩阵拼接融合成为$ 12\times 44 $的新特征, 再使用ResNet152、Inception-v3和Inception-ResNet-v2迁移特征融合后的分类方法, 进行模型的训练. 最后使用预测评价系统对预测方法进行评价.

2.2.1 评价系统

本文遵循文献[20]中提到的评价系统, 采用癫痫发作预测范围(seizure prediction horizon, SPH)和发作发生期(seizure occurrence period, SOP)两个指标, 其中SPH是指癫痫发作前的预测时间段, SOP指的是癫痫一定会发作的时间段. 为了避免算法误报, 当连续预测到5次癫痫发生后, 警报才会响起. 为了做出正确的预测, 癫痫发作必须在SPH之后且在SOP范围内. 同样, 当预测系统检测到癫痫发作, 但在SOP期间没有发生癫痫时, 会出现假警报. 当警报响起时, 将持续到SOP结束. 在临床医疗应用中, SPH必须足够长, 以允许足够的干预或预防措施. SPH也称为干预时间. 相反, SOP不宜太长, 以减轻患者的焦虑, 所以选择SPH较长, SOP较短. 本文采用的预测评价系统采用两种发作预测时间: 1) $ {\rm SPH}=10 $ min, $ {\rm SOP}=10 $ min; 2) $ {\rm SPH}=30 $ min, $ {\rm SOP}=10 $ min.

本文采用的评价指标为预测敏感性(prediction sensitivity, PSEN)和错误预测率(false prediction rate, FPR). 在医学统计中, 敏感性是指在确诊病症的不误诊(假阳性)概率, FPR是每小时错误预测的次数, 单位为h$ ^{-1} $或者(1/h).

$ \begin{align} &{\rm PSEN}=\dfrac{\text{正确预测癫痫次数}}{\text{癫痫总共发生次数}}\times 100 \%, \end{align} $

(8)

$ \begin{align} &{\rm FPR}=\dfrac{\text{错误报警次数}}{\text{预测时间}}. \end{align} $

(9)

2.2.2 分类结果

为了验证通道筛选的有效性, 本文对比了无筛选方法的分类性能, 由表 2的实验结果发现PSDED作为特征图的分类算法准确率会高于SMD方法. PSDED相较于SMD在发作前两种状态的准确率有4.3 %和3.5 %的增益, 较为明显. 但在发作间期和发作期的准确率没有显著增益. 结合通道筛选后的PSDED方法, 在总体准确率上有了显著的提高, 尤其在发作前两种状态的准确率会显著提升. 将通道筛选后的PSDED与SMD融合, 结合DCNN训练后的准确率得到了进一步的提升, 发作间期达到98.2 %, 发作前I达到95.4 %, 发作前II达到95.1 %, 以及发作期为98.6 %. 筛选融合后的方法在癫痫发作期有大幅提升, 并且在难以区分的两种状态(发作前30 min、发作前10 min)得到较大提升.

2.2.3 预测结果

本文所提出方法的分类性能得到了显著的提升, 但本文重点为算法的预测能力. 这里采用的训练及测试策略与分类方法有区别, 不再对全部样本进行随机划分7:3的训练集和测试集. 一个癫痫发作周期为: 发作间期、发作前期、发作期, 每位患者至少3次以上癫痫发作周期, 选择患者前两次癫痫发作周期为训练样本, 包含60 min的癫痫发作前期数据和周期内所包含的发作间期数据, 没有训练的数据用作测试样本, 预测训练模型如图 8所示.

对比通道筛选前和通道筛选后算法对EEG预测的性能, 表 3的结果显示: 当$ {\rm SPH}=10 $ min, $ {\rm SOP}=10 $ min时, 采用PSDED方法结合深度卷积神经网络的预测56次癫痫中有47次预测成功, 所以敏感性达到了83.7 %, 误检率达到了0.17/ h. 而在通道筛选后的PSDED方法的预测敏感性提高到93.63 %, 误检率下降到0.09/ h.

将PSDED与SMD进行特征图融合, 再进行模型的训练, 分析融合后的结果与通道筛选后得到的结果, 如表 4所示. 结果表明: 特征图融合后的PSEN可以达到92.03 %, 误检率达到了0.07/ h; 而经过通道筛选后的PSEN为96.66 %, 误检率为0.03/ h, 反映出经过通道筛选后会有更好的预测性能.

当$ {\rm SPH}=30 $ min, $ {\rm SOP}=10 $ min时, 结果如表 5所示. PSDED与SMD融合后结合DCNN预测的方法, 在通道筛选前的算法预测56次癫痫中有48次预测成功, 所以敏感性达到了85.71 %, 误检率达到了0.13/ h; 而在通道筛选后的PSDED方法在预测敏感性提高到93.17 %, 误检率下降到0.05/ h.

从EEG中提取的特征被广泛用于预测癫痫发作. 本文对比了筛选通道前后算法以及特征图融合方法的分类和预测性能: 1) 筛选通道后的分类和预测性能会好于原始算法, 说明进行通道选择后, 可以舍去部分会影响分类决策的干扰特征; 2) 通过特征图融合, 结合了EEG能量变化的特征以及通道间信息交互的特征, 从两个角度体现了EEG的特征, 更好地体现了不同状态下的EEG特征.

2.2.4 预测性能比较

对比同样适用CHB-MIT数据集的算法, 如表 6所示. 其中文献[21]设定SPH为0, 而SOP设定为50 min, 并解释为报警后50 min内癫痫发作则算成功预测一次. 文献[7]也设定SPH为0, 只考虑SOP, 与本文采用的预测系统有所差别. 对比表明: 在考虑SPH的情况下, 本文的预测算法SOP时间会少于其他研究, 这样更有利于及时提醒患者并减少等待时间, 且得到了较好的预测敏感度和较低的误检率, 模型有预测癫痫发作的能力.

3 讨论分析

通过脑电通道优化筛选后得出了更好的癫痫发作预测结果, 为了验证通道筛选有利于癫痫发作的预测, 对患者的耦合性进行了详细的分析. 使用多尺度符号化排列传递熵对24个病历记录进行分析, 并提取每个病历的重要通道. 共计19例病历中C3-P3通道在发作期及发作间期的变化较大, 被筛选为重要通道, 通道呈现图 9(a)所示变化的趋势, 癫痫发作时熵值会下降, 通道间的耦合强度会下降. 其余的病历变化趋势如图 9(b)所示, 发作期的熵值无明显变化.

18例病历在T8-P8通道被筛选为重要通道, 熵值变化趋势如图 10(a)所示, 癫痫发作期的熵值会上升, 通道间的耦合强度会增强. 3例病历中T8-P8通道变化趋势如图 10(b)所示, 发作期与发作间期的熵值变化不明显, 其余病历T8-P8通道变化如图 10(c)所示, 发作期的熵值会下降, 耦合强度会减弱. 由于CHB-MIT没有具体说明数据集中有哪些癫痫类型的患者, 这里通过实验结果推测出大部分患者是颞叶区癫痫发病的类型, 所以癫痫时熵值会上升; 部分患者癫痫发作与颞叶区关系较小.

分析24例病历, 大部分患者在癫痫发作期时T8-P8通道的排列传递熵值大于发作间期, 有更强耦合性. T8-P8所属颞叶区, 最常见的是颞叶癫痫, 多表现为精神游离、失忆等^[23]. 而大部分患者的C3-P3通道在发作期的熵值会小于发作间期, 患者基本会失去对身体的控制. 有研究表明: 癫痫患者的左前额叶子的同步强度比正常人弱, 患者发作期的左侧额叶和右侧额叶的耦合性降低^[24]. 文献[25]通过部分定向相干法建立脑网络对CHB-MIT癫痫数据加以分析, 实验结果表明发作期的颞叶区相关性变强, 而额叶区、顶叶区及枕叶区耦合性降低, 且重点分析T8-P8相关性在癫痫发作时会增强. 通过对部分重要脑电通道的耦合强度和病理角度的详细分析, 基本可以解释筛选脑电通道可以带来更好的预测性能, 原因可能是筛选通道后, 提取的特征可以更好地反映不同状态下脑电信号的差异性, 以便深度卷积神经网络可以训练出更好的预测模型.

4 结论

本文利用多尺度符号化排列传递熵对患者发作间期及发作期各通道EEG的耦合强度的变化, 生成同步矩阵并通过显著性分析进一步筛选重要EEG通道, 减少了冗余的特征对预测算法的影响. 针对优化筛选通道后的脑电信号再生成PSDED特征图和同步矩阵图SMD. 结合通道能量特征和不同通道间信息交互的特征, 产生融合特征图, 使用拥有更好学习效率的ResNet152、Inception-v3和Inception-ResNet-v2三个DCNN模型, 将提取得到的特征进行拼接后再输入后续的全连接层进行训练和降维, 最后再经过softmax进行分类. 对训练好的模型, 在预测系统中进行测试, 采用了SPH为10 min, SOP为10 min的设定, 最后结果预测敏感性达到了96.66 %, 误检率达到0.03/ h; SPH为30 min, SOP为10 min的设定, 最后结果预测敏感性达到了93.17 %, 误检率达到0.05/ h. 预测的性能超越了没有进行通道优化筛选且没有融合特征的算法. 通过符号化排列传递熵算法通道优化选择后, 再进行有针对性的训练为癫痫的检测与预测带来了新的优化方法. 此研究提高了癫痫发作的预测性能, 有益于患者的预防和治疗.