利用高频脑电的局灶性癫痫患者癫痫发作检测

王栋,李扩,刘晓芳,闫相国,王刚

(1.西安交通大学生物医学信息工程教育部重点实验室,710049,西安;2.西安交通大学第一附属医院神经外科,710061,西安)

癫痫是一种常见的慢性脑部疾病,其发作表现为大脑神经元过度同步化放电,并导致短暂性中枢神经系统功能失常。脑电图(EEG)是一种非侵入性的测量脑电活动的工具,其中蕴含大量的脑功能信息,对大脑疾病的诊断具有很高的价值,也可以用作癫痫发作的检测[1]。在多数情况下,癫痫发作无法在短期内预测,因此需要在较长的时间内连续记录脑电。长程脑电图监测可以有效地提供关于大脑活动、癫痫发作频率等的信息,有助于临床医生对癫痫的诊断和治疗。长程脑电经过专家的目测检查,可以确定癫痫发作的情况,但是这是一项相当乏味而耗时的任务,繁重的工作量容易引起医护人员疲劳而导致人工检测结果的准确性下降[2],并且这很依赖于临床医生的自身经验和主观判断,不同的医护人员所作的诊断有可能不一致,结果的可重复性较低。除此之外,如果在癫痫发作后较短的时间内检测到,也利于医护人员采取相应的措施来减少对患者的伤害。因此,计算机自动检测技术一直是癫痫脑电图监测工作中的迫切需求。

近几十年来,出现了各种各样的癫痫自动检测的方法,其中特征提取方法主要包括时域、频域、时频域和非线性的方法[3]。多方研究表明,脑电波信号表现为一种非线性、非平稳的随机过程,因此关联维[4]、Lyapunov指数[5]、熵[6]等度量特征已经被引入到信号分析和特征提取中,小波变换也被广泛地应用到了该领域[7-8]。除了单一通道的脑电,越来越多的人开始关注多通道脑电之间的相互关系,关注脑区各部分的信息交流[9-10]。有向传递函数(DTF)是通过脑电来计算脑区信息交流的方法之一[11]。现有的大多数的癫痫发作检测针对的都是低频段,忽略头皮脑电的高频成分。

本文首先利用小波分解将脑电信号中的高频γ波段提取出来,利用该波段脑电信号的流出信息强度作为特征来进行癫痫发作检测。结果表明,癫痫发作前后,高频段的脑电流出信息强度差异很大,利用长程脑电的高频成分对局灶性癫痫患者癫痫发作检测有着很高的灵敏性和特异性。

1 数据和方法

1.1 癫痫患者和脑电数据

本文长程脑电数据是从西安交通大学第一附属医院10例临床癫痫患者身上采集的。数据采集仪器是NIHON KOHDEN公司的EEG-1100脑电图机,使用19导电极,电极位置按照国际10—20标准,以CZ电极为参考电极,数据的通带截止频率为0.5～60 Hz,采样频率为200 Hz。癫痫发作的开始点和终止点由西安交通大学第一附属医院的两位具有临床经验的医生通过视觉检测和评估患者的实际临床记录和表现进行了标注。我们又人工剔除了原始脑电数据中的肌电等伪迹脑电。本文获得的已采集脑电数据的总时长为121.8 h,包含了44次发作,其中癫痫发作期的平均时间为98.7 s。表1给出了每位患者的相关脑电数据信息。

1.2 癫痫检测算法

癫痫脑电数据分为发作期脑电和非发作期脑电,发作检测的目的是在长程脑电中识别发作期脑电。已有研究表明,癫痫发作与脑电的高频成分有关系,癫痫发作时会有高频脑电成分出现。在癫痫发作期和非发作期,头皮脑电的高频成分在不同电极间的信息流有可能是不同的,所以本文利用DTF算法来计算高频脑电在不同脑区的信息流通情况。信息流是有方向的,即流出信息和流入信息,本文将各个电极的流出信息强度作为特征,用于癫痫检测。本文利用高频脑电成分进行癫痫发作检测的流程如图1所示。首先将每位病人19通道的脑电数据进行分割,分为发作脑电数据和非发作期脑电数据。对每个通道的脑电在滑动的时间窗内用小波分解的方式提取高频成分(γ波段),然后对19通道的高频脑电建立多变量自回归模型(MVAR)。利用有向传递函数来提取该频段上各通道的信息交流特征,求出各通道的流出信息强度用于特征降维,最后经SVM分类,通过五重交叉验证得到检测结果。

1.2.1 高频脑电信号的提取 为了使用小波变换的方式提取出高频波段的信号,根据小波分解的原理,当信号的频率为256 Hz时,对小波进行两层分解,即可得到大于32 Hz的脑电信号。小波分解中小波的选取非常重要,Db小波群对非平稳信号(比如脑电信号)具有更好的分解效果[7-8],而db4小波更适合分析癫痫脑电信号[12]。因为原始脑电数据的采样频率为200 Hz,所以先进行了一个256 Hz的重采样,然后选择db4小波对各个通道的脑电信号进行两层的小波分解,就能得到高频脑电信号(>32 Hz)。

1.2.2 基于有向传递函数的特征提取 由于采集的脑电信号最小频率为0.5 Hz,为了保证能达到最低分辨率,对脑电信号的分割需要不少于2 s[10]。本文对各个频段的脑电信号进行2 s无重叠的分割,利用有向传递函数进行特征提取的第一步需要建立多变量自回归模型(MVAR),然后在求得MAVR模型系数的基础之上,利用有向传递函数(DTF)进行特征值提取,该特征能够反映不同导联脑电信号之间的信息流动情况。

定义t时刻的N导联(癫痫脑电的导联数N=19)脑电信号为

X(t)=[X1(t),X2(t),…,XN(t)]T

(1)

式中:Xn(t)(n=1,2,…,N)是第n个通道的脑电信号。通过建立多通道自回归模型,脑电序列可以表示如下

(2)

式中:p为MVAR模型的阶数;Ar为N×N的系数矩阵,r=1,2,…,p;E(t)为估计误差,理想情况下是均值为0的非相关白噪声。模型的阶数p可以通过Akaike信息准则(AIC)来确定[13]

F(p)=ln|∑(p)|+2pN2/M

(3)

式中:∑(p)为残差的协方差矩阵;N为通道数;M是全部用于计算的数据点数。选择最小的F对应的阶数,系数矩阵Ar的估计可以用Levison Wiggins Robinson(LWR)算法[14]求得。

对式(2)两边作傅里叶变换,得到

(4)

(5)

传递矩阵定义为

(6)

进而可以获得在频率f上从导联j到导联i的信息流强度为

Dij(f)=

(7)

式中:Hij(f)为矩阵H(f)的第i行第j列元素;hi(f)为矩阵的第i列;Dij(f)表示在频率f时,从导联j到导联i的信息流的强度和方向。

1.2.4 支持向量机分类 支持向量机(SVM)是一种新的机器学习方法,基本思想是:对于线性不可分样本,经非线性变换,将其映射到另一个高维空间中,在变换后的空间中寻找一个最优的分界面(超平面),使之线性可分[15]。支持向量机对于二分类问题有着明显的优势,本文选择径向基核函数用于发作期脑电和非发作期脑电的分类,采用五重交叉验证的方法得到最终结果:将每一个病人的脑电数据平均分为5组,每次用4组数据做训练,一组数据用来测试,重复5次,保证每组都只测试1次,将5次的结果进行平均作为最终测试结果。

2 结 果

算法的性能是根据每个患者临床信息的分析统计结果的平均值来评估的,以避免由于单个患者有多次癫痫发作而引起的可能的偏差。通过比较本文算法对EEG片段做出的癫痫发作、非发作标记以及医生对此做出的标记的差别来计算算法的性能。其中,真阳性TP是算法判断正确的癫痫发作期片段,假阳性FP是算法判断错误的癫痫发作期片段,真阴性TN是算法判断正确的癫痫非发作期片段,假阴性FN是算法判断错误的癫痫非发作期片段。算法的性能可以由以下指标进行评估[14,16],这些指标的定义如下。

(1)正确率

%

(8)

(2)选择性

%

(9)

(3)灵敏性

%

(10)

(4)特异性

%

(11)

(5)平均检出率

ADR=0.5(SPE+SEN)×100%

(12)

本文获得的已采集脑电数据的总时长为121.8 h,包含了44次发作,其中癫痫发作期的平均时间长度为98.7 s。以2 s时间窗作为分割,训练样本的平均非发作期片段数为11 959,平均发作片段数为157,测试样本的平均非发作期片段数为3 313,平均发作期片段数为39。表2为利用γ波段进行癫痫发作检测的结果,10位病人的平均正确率为98.4%,平均选择性为60.7%,平均敏感性为93.4%,平均特异性为98.4%,平均检出率为95.9%。

表2 基于高频脑电的局灶性癫痫患者 癫痫发作检测结果

3 讨 论

3.1 基于DTF不同波段癫痫检测结果对比

图2 基于有向传递函数(DTF)算法的不同波段癫痫检测结果

为了讨论不同子频带的特征对癫痫检测结果的影响,本文用不同子频带的特征进行SVM分类。首先利用5层小波分解得到4种典型的脑电信号波段,即δ波(0.5～4 Hz)、θ波(4～8 Hz)、α波(8～16 Hz)、β波(16～32 Hz),然后采用DTF算法求出这4个子波上的流出信息特征,进行癫痫发作检测。将所有频带的检测结果进行比较,如图2所示。通过单因素方差分析,全频带分类结果的5项指标与γ波段分类结果的5项指标无显著性差异(P>0.05)。低频段的4个子波段分类结果的5项指标与γ波段分类结果的5项指标有着显著性差异(P<0.05)。这说明了不同子波段对于癫痫脑电的检测都提供了不同量的特征;另一方面,说明了癫痫发作与非发作相比,γ波段对于癫痫脑电的检测提供更多的特征,即γ波段的信息流强度变化显著。

3.2 基于部分有向相干(PDC)不同波段癫痫检测结果对比

为探究癫痫发作前后γ波段的信息流动变化是否显著,本文采用了另一种脑功能连接的PDC算法。相较于DTF算法,PDC算法只会发现一些直接的因果关系,比如A、B、C 3者之间的关系,用PDC可能只能发现A到C的有向连接,而用DTF可能会发现A到B、B到C和A到C都有因果性联系。因此,本文采用PDC算法进行脑功能连接,求出脑电各个频段上的流出信息强度特征,用以癫痫发作的检测,检测结果如图3所示。单因素方差分析的结果与DTF算法检测结果基本一致,其中全频带分类结果的5项指标与γ波段分类结果的5项指标无显著性差异(P>0.05),其他4个子波段分类结果的5项指标与γ波段分类结果的5项指标有着显著性差异(P<0.05)。

图3 基于部分有向相干(PDC)算法的不同波段癫痫检测结果

3.3 癫痫发作前后流出信息强度变化

为了探究癫痫发作前后流出信息强度的变化情况,寻找DTF和PDC算法在γ波段癫痫发作检测效果优于其他子波段的原因,将癫痫发作前后各个波段上的各通道的流出信息强度进行了对比,并选取了一位病人(患者2),给出了在一次癫痫发作前后10 s的不同波段上的脑区流出信息强度的脑地形图,如图4所示。癫痫发作的前后10 s,δ波、θ波、α波、β波、γ波的流出信息强度变化不同。癫痫发作前10 s,各频段的流出信息强度也不相同,δ波和θ波流出信息强度主要散布在右侧前额叶、右侧颞叶和枕叶;α波流出信息强度主要散布后侧颞叶和枕叶;β波流出信息强度集中在额叶、顶叶和枕叶;γ波流出信息强度都很小,全频带流出信息强度也都很小,少量集中在枕叶。癫痫发作第10 s,δ波流出信息量主要在后额叶和右侧后颞叶;θ波流出信息强度主要集中在后额叶、顶叶和左侧前额叶;α波、β波、γ波和全频带流出信息强度集中的脑区位置很相似,主要在后额叶和顶叶。全频带的流出信息强度变化和γ波段变化很相似,而且它们在癫痫发作前后流出信息强度变化差异很大。特别在γ波,在癫痫发作之前,各个脑区的流出信息强度都比较弱,发作开始后,出现明显的信息流在某些脑区汇聚和加强。这也可能是图2中用γ和全频段检测结果的5项指标无显著性差异且比其他子频段的检测结果好的原因。

Ayala等用颅内脑电(iEEG)的γ波段的功率谱进行癫痫检测取得了不错的效果,并指出癫痫发作与γ波的功率谱有一定关系[17]。Lu等用头皮脑电的高频部分对癫痫病灶进行定位,也说明了癫痫的发作与高频脑电活动有关[18]。本文实验结果显示:γ波在癫痫发作前后流出信息强度变化显著,说明了癫痫发作与γ波甚至更高的高频振荡有着密切的关系。对于局灶性癫痫患者,在癫痫发作时高频成分的变化,导致了高频波段的信息流在某一脑区聚集和增强,而且相对规律,这种特征可以被捕捉到并用于癫痫发作检测。相对于其他子波段,在癫痫发作前后,信息流的聚集和变化情况表现为不规律状态,故没有用高频成分检测效果好。

3.4 滑动时间窗对检测结果的影响

本文分别在滑动时间窗为2、3和5 s时进行了癫痫发作检测。当窗口为3 s时,检测结果的平均正确率为98.1%,平均选择性为55.2%,平均敏感性为93.4%,平均特异性为98.2%,平均检出率为95.7%。当窗口为5 s时,检测结果的平均正确率为98.4%,平均选择性为61.2%,平均敏感性为96.6%,平均特异性为98.5%,平均检出率为97.6%。结果显示,检测的各项指标无显著性差异(P>0.05),所以我们认为该方法与时间窗无关。为了减少运算复杂度和在癫痫发作后最短的时间检测到,本文采用2 s的时间窗。

图4 病人2的癫痫发作前后各个波段的脑电流出信息强度脑地形图

3.5 与其他检测方法比较

图5 不同检测方法的结果对比

选择AR模型[19]和近似熵[20]两种现存的癫痫发作检测方法在相同数据集上进行实验。计算近似熵时时间窗为2 s,然后将每一导的信号求近似熵,将19导得到的特征进行组合。对于AR模型方法,时间窗口也是2 s,设置的阶数为3,并将19导得到的特征进行组合。对比结果如图5所示,AR模型方法检测结果的平均正确率为72.6%,平均选择性为4.9%,平均敏感性为51.2%,平均特异性为72.9%,平均检出率为62.1%。近似熵方法检测结果的平均正确率为93.8%,平均选择性为23.9%,平均敏感性为85.5%,平均特异性为93.8%,平均检出率为89.6%。本文方法检测结果的5项指标与AR模型和近似熵方法的检测结果相比有显著的提升(P<0.05)。

4 结 论

本文利用长时程头皮脑电的高频成分对局灶性癫痫患者进行癫痫发作检测,将多通道的脑电信号通过小波变换提取γ波段,然后在移动的时间窗内建立MVAR模型,采用有向传递函数算法求出γ波段上不同导联的流出信息强度,最后经过SVM分类器进行分类,识别出发作期脑电和非发作期脑电。通过五重交叉验证得到的平均正确率为98.4%,平均选择性为60.7%,平均敏感性为93.4%,平均特异性为98.4%,平均检出率为95.9%。通过和其他波段脑电信号癫痫发作检测结果比较,发现γ波段的检测效果最好。另外,通过对癫痫发作前后不同脑区在不同子频段上的流出信息强度进行对比与分析,发现癫痫发作时在γ波段上信息流在某一脑区集中和增强,而且相对规律,这种特征能被本文方法捕捉并可以实现癫痫发作的检测。对大脑高频成分活动的研究可能对癫痫发作起因、发作检测、发作预测和病灶定位有着重要意义。

[1] SUBASI A. Epileptic seizure detection using dynamic wavelet network [J]. Expert Systems with Applications, 2005, 29(2): 343-355.

[2] SHARMA R, PACHORI R B. Classification of epileptic seizures in EEG signals based on phase space representation of intrinsic mode functions [J]. Expert Systems with Applications, 2015, 42(3): 1106-1117.

[3] ACHARYA U R, SREE S V, SWAPNA G, et al. Automated EEG analysis of epilepsy: a review [J]. Knowledge-Based Systems, 2013, 45(3): 147-165.

[4] GENG S, ZHOU W. Nonlinear feature comparison of EEG using correlation dimension and approximate entropy [C]∥2010 3rd International Conference on Biomedical Engineering and Informatics. Piscataway, NJ, USA: IEEE, 2010: 978-981.

[5] GULER N, UBEYLI E, GULER I. Recurrent neural networks employing Lyapunov exponents for EEG signals classification [J]. Expert Systems with Applications, 2005, 29(3): 506-514.

[7] KHAN Y U, GOTMAN J. Wavelet based automatic seizure detection in intracerebral electroencephalogram [J]. Clinical Neurophysiology, 2003, 114(5): 898-908.

[8] 马东华, 郑旭媛, 王真. 基于形态成分分析的癫痫脑电棘波检 [J]. 生物医学工程杂志, 2013, 30(4): 710-713. MA Donghua, ZHENG Xuyuan, WANG Zhen. Detection of epileptic spike wave in EEG signal Based on morphological component analysis [J]. Journal of Biomedical Engineering, 2013, 30(4): 710-713.

[9] VAN MIERLO P, PAPADOPOULOU M, CARRETTE E, et al. Functional brain connectivity from EEG in epilepsy: seizure prediction and epileptogenic focus localization [J]. Progress in Neurobiology, 2014, 121: 19-35.

[10]WANG G, SUN Z, TAO R, et al. Epileptic seizure detection based on partial directed coherence analysis [J]. Journal of Biomedical and Health Informatics, 2016, 20(3): 873-879.

[11]KAMINSKI M J, BLINOWSKA K J. A new method of the description of the information flow in the brain structures [J]. Biological Cybernetics, 1991, 65(3): 203-210.

[12]ADELI H, ZHOU Z, DADMEHR N. Analysis of EEG records in an epileptic patient using wavelet transform [J]. Journal of Neuroscience Methods, 2003, 123(1): 69-87.

[13]AKAIKE H. A new look at the statistical model identification [J]. IEEE Transactions on Automatic Control, 1974, 19(6): 716-723.

[15]吴烜, 李京. 基于支持向量机的测厚仪CS值电压漂移故障判定及处理 [J]. 计算技术与自动化, 2014, 33(1): 42-45. WU Xuan, LI Jing. Thickness gauge CS value voltage drift down approach based on support vector machines algorithm determines [J]. Computing Technology and Automation, 2014, 33(1): 42-45.

[16]AARABI A, GREBE R, WALLOIS F. A multistage knowledge-based system for EEG seizure detection in newborn infants [J]. Clinical Neurophysiology, 2007, 118(12): 2781-2797.

[17]AYALA M, CABRERIZO M, JAYAKAR P, et al. Subdural EEG classification into seizure and nonseizure files using neural networks in the gamma frequency band [J]. Journal of Clinical Neurophysiology, 2011, 28(1): 20-29.

[18]LU Y, WORRELL G A, ZHANG H C, et al. Noninvasive imaging of the high frequency brain activity in focal epilepsy patients [J]. IEEE Transactions on Biomedical Engineering, 2014, 61(6): 1660-1667.

[19]CHISCI L, MAVINO A, PERFERI G, et al. Real-time epileptic seizure prediction using AR models and support vector machines [J]. IEEE Transactions on Biomedical Engineering, 2010, 57(5): 1124-1132.

[20]YENTES J M, HUNT N, SCHMID K K, et al. The appropriate use of approximate entropy and sample entropy with short data sets [J]. Annals of Biomedical Engineering, 2013, 41(2): 349-65.