CN101002675A - 一种基于二阶c0复杂度的癫痫发作预测方法 - Google Patents

Abstract

本发明属临床医学和神经系统技术领域，具体为一种可以自动地对癫痫症发作时间进行预测的方法。该方法通过分析癫痫患者的颅内脑电信号的二阶C₀复杂度，找出癫痫发作前期二阶C₀复杂度的变化特征，再利用统计估计算法识别出这种变化，从而预测癫痫发作时间。本发明具有较高的预测准确率，对于癫痫患者或医生及时地采取有效预防措施具有重要的意义。

Description

一种基于二阶C0复杂度的癫痫发作预测方法

技术领域

本发明属于临床医学和神经系统技术领域，具体涉及一种基于颅内脑电信号的癫痫发作时间预测方法。

技术背景

癫痫发作预测是近年来在临床医学和神经系统科学研究领域中备受关注的问题。如果癫痫发作可以得到有效的预测，则可以研制出能够供临床引用的报警和治疗设备，使患者或医生提前采取有效的预防措施，如电刺激治疗或药物治疗等^[1]。

目前对癫痫发作的研究主要是基于对其脑电信号的分析。在过去的十几年间，许多研究小组已经验证了在癫痫发作前存在一个“发作前期”，这一阶段可能持续几分钟，甚至更久^[2][3]。这些研究成果提示我们如果能及时捕捉到这个发作前的状态，找出其变化特征，就有可能实现癫痫发作的预测。

目前已提出的癫痫发作预测方法大致分为三类：一类是时间域分析法，利用癫痫脑电数据的时域统计特性进行发作预测；一类是频域分析法，将脑电信号划分成不同的频率成分，通过频谱分析得到癫痫发作前兆；第三类为非线性方法，利用非线性动力学理论研究发作前脑电信号波形的变化特征进行预测，主要指标有相关维(Correlation Dimension)^[4]、动态相似性(Dynamic Similarity)^[5]、熵(Entropy)^[6]、相位同步(Phase Synchronization)^[7]、交互相关性(Cross Correlation)^[7]、累积能量(Accumulated Energy)^[8]等。

以往文献中的预测结果有些虽然看起来很理想，但大多是针对自己采集的有限数据得到的，当运用于不同数据库时，其可预测性和预测效果都遭到质疑。2002年4月，Lehnertz和Litt在德国波恩主持举行了第一届癫痫预测国际合作讨论会，有7个国家16个研究中心的51位研究人员参加了会议。Lehnertz和Litt事先向与会人员提供了统一的癫痫脑电数据，以利于比较不同研究者所用的不同方法对癫痫发作预测的有效性。其中8个研究组的报告中，有三个给出了肯定的结果，有四个报告给出了否定的结果，另有一个报告则按照其使用方法的不同而分别给出肯定和否定的结果^[9]。

针对以上问题，本发明提出了基于二阶C₀复杂度的癫痫预测方法，并将该方法运用在网上公开的FDM(Freiburg Center for Data Analysis and Modeling)数据库中^[10]，该方法可以有效地预测癫痫发作时间。

本发明中所涉及的癫痫脑电信号均为颅内脑电信号，它是在手术前，将电极植入癫痫患者的颅内而采集得到的脑电信号。

下面介绍与本发明相关的一些概念：

1、复杂度及C₀复杂度

复杂度是作用于时间序列的一种指标。一段时间序列通过一个特定的复杂度算法处理，可以得到一个复杂度值，该值表现序列的非规则程度，一般地，不规则变化序列的复杂度较高，像常序列、周期序列这些模式比较简单的序列复杂度值较低。在脑电信号分析中，复杂度的分析逐渐受到广泛的关注，如近似熵、Kolmogorov熵和Lempel&Ziv复杂度^[11]都已被用于脑电分析中。这些复杂度算法的缺点是它们都需要对序列做粗粒化，即通过对序列取均值，将序列中每个大于均值的数看作是1，而小于均值的视作0，进而转化为求这个0、1序列的复杂度。这一粗粒化的过程可能导致丢失许多有意义的细节，因此，为克服粗粒化的缺点，顾凡及等提出了C₀复杂度的概念^[12][13]。

C₀复杂度在数据长度较短的情况下也可以得出鲁棒性较好的估计值，并且它对连续信号也有定义，从而无需对原始信号作像二值化这样的过分粗粒化。其主要思想是将信号分解成规则成分和不规则成分两部分，C₀复杂度定义为非规则成分在原信号里所占的比例。

2、二阶复杂度

一阶复杂度描述的是序列的非规则程度，因此它可以反映时间序列的随机性。对于一个平稳或者准平稳时间序列，如果我们用一个滑动窗计算每段原始序列的复杂度值，就可以得到一个由复杂度值组成的新序列。求这个复杂度序列的复杂度定义为二阶复杂度^[14]。

对于一个周期性时间序列而言，其每个时间窗求得的一阶复杂度应近似为零；对于一个严格平稳的时间序列，每个时间窗求得的一阶复杂度近似相同，即接近常量，则从宏观角度来看，其二阶复杂度为一个非常小的值。然而，对于非平稳信号而言，给定长度时间窗计算得到的复杂度是随时间而变化的，则其二阶复杂度为一个较大的值，且该值的大小依赖于序列的非平稳程度。因此，二阶复杂度描述的是时间序列的非平稳性。

参考文献

[1]Witte H，Iasemidis LD，Litt B.Special issue on epileptic seizure prediction.IEEE Trans.Biomed.Eng.，2003，50(5)：537-539

[2]Katz A，Marks DA，McCarthy G，Spencer SS.Does interictal spiking change prior toseizures？Electroencephalogr Clin Neurophysiol，1991，79(2)：153～156

[3]Schiff S，Colella D，Jacyna G.Brain chirps：spectrographic signatures of epileptic seizures.Clin Neurophysiol，2000，111(6)：953～958

[4]Elger CE，Lehnertz K.Seizure prediction by non-linear time series analysis of brainelectrical activity.Eur J Neurosci，1998，10(2)：786～789

[5]Navarro V，Martinerie J，Le Van Quyen M，Clemenceau S，Adam C，Baulac M，Varela F.Seizure anticipation in human neocortical partial epilepsy.Brain，2002，125(3)：640～655

[6]Van Drongelen W，Nayak S，Frim DM，Kohrman MH，Towle VL，Lee HC，McGee AB，Chico MS，Hecox KE.Seizure anticipation in pediatric epilepsy：use of Kolmogorov entropy.Pediatr Neurol，2003，29(3)：207～213

[7]Mormann F，Kreuz T，Andrzejak RG，David P，Lehnertz K，Elger CE.Epileptic seizures arepreceded by a decrease in synchronization.Epilepsy Res，2003，53(3)：173～185

[8]Esteller R，Echauz J，D′Alessandro M，Worrell G，Cranstoun S，Vachtsevanos G，Litt B，Continuous energy variation during the seizure cycle：towards an on-line accumulated energy.Clin Neurophysiol，2005，116(3)：517～526

[9]Ebersole JS，In search of seizure prediction：a critique.Clin Neurophysiol，2005，116(3)：489～492

[10]http://www.fdm.uni-freiburg.de/groups/timeseries/epi/EEGData/[11]Jia WY，Kong N，Li F，Gao XR，Gao SK，Zhang GJ，Wang YP，Yang FS，An epilepticseizure prediction algorithm based on second-order complexity measure.Physiol Meas，2005，26：609～625

[12]陈芳，顾凡及，徐京华，刘增荣，刘仁，一种新的人脑信息传输复杂性的研究。生物物理学报，1998，14：508～512

[13]Shen EH，Cai ZJ，Gu FJ，Mathematical foundation of a new complexity measure.AppliedMathematics and Mechanics，2005，26(9)：1188～1196

[14]Gu FJ，Shen EH，Meng X，Cao Y，Higher order complexity of time series.International J.of Bifurcation and Chaos.2004，14(8)：2979～2990

发明内容

本发明的目的在于提出一种基于二阶C₀复杂度的癫痫发作预测方法。该发明可以有效地自动预测癫痫发作时间，将该算法应用在FDM数据库中，准确率达94.3％以上，平均预测时间为58.4秒。

本发明提出的癫痫发作预测方法，具体步骤如下：

1、脑电数据预处理

由于实验用的癫痫数据通常为术前采集的颅内脑电数据，伪迹干扰较小，因此我们不对脑电数据进行伪迹去除处理。

通过频谱分析和大量实验，我们发现癫痫发作脑电信号的复杂度在20～80Hz频率段变化趋势更为明显。附图1给出了一段采样率为256Hz的癫痫数据在发作前(pre-ictal)、发作过程(ictal)和发作后(post-ictal)的频谱图，从该图中我们可以看出，在发作前和发作结束后，其频谱主要集中在20Hz以下的部分，而在发作过程中，20～80Hz频率段的频谱明显增多。这说明癫痫发作时，可能大量的信息蕴藏在一个频谱相对较高的频率段中。

因此，本发明中，我们首先对病人的癫痫脑电数据x(t)作20～80Hz的带通滤波预处理，得到滤波后的信号

2、C₀复杂度的计算

(1)对预处理后的脑电信号

，计算其C₀复杂度，得到一阶C₀复杂度序列{Cf_i，i＝1，2，…，N}。一般计算C₀复杂度时窗口长度为1-2秒。所以考虑到1-2秒钟的脑电信号可以近似地看成平稳信号，其C₀复杂度趋于一个稳定的值，所以计算C₀复杂度时窗口长度选择1-2秒的数据段。同时，计算复杂度时，每个复杂度值在时间轴上对应数据窗口的终点，这样计算预测点的二阶复杂度所牵涉到的脑电数据就不会落到癫痫发作区。N＝[Y/W]，Y为脑电信号x(t)的长度，W为窗口长度。

(2)计算复杂度序列{Cf_i，i＝1，2，...，N}的C₀复杂度，得到二阶C₀复杂度序列(Cs_i，i＝1，2，...，M}。M＝[N/W]，[]表示取整数。

3、预测算法

本发明采用较为直观的简单统计估算方法预测癫痫发作。对于上述步骤求得的二阶C₀复杂度序列{Cs_i，i＝1，2，…，M)，设滑动窗口长度为L个点，每次移动为D个点，计算平均值m_j为

$m_j=\frac{1}{L}\underset{i=1}{\overset{L}{\mathrm{\Σ}}}C{s}_i,---\left(1\right)$

标准差d_j为

$d_j={\left(\frac{1}{L}\underset{i=1}{\overset{L}{\mathrm{\Σ}}}{({\mathrm{Cs}}_j-m_j)}^2\right)}^{1/2}.---\left(2\right)$

则可以得到两个序列{m_j，j＝1，2，…，K}和{d_j，j＝1，2，…，K}，其中 $K=[\frac{M-L}{D}+1].$

若m_j-1-m_j＞2d_j-1，则发出报警，记录此时相应的二阶C₀复杂度值所对应的时刻。

一般取滑动窗口长度L取为25-35个点，每次移动D取为10-20个点。

本发明的优点：

本发明是基于二阶C₀复杂度的癫痫发作预测算法，其优点在于：一是充分利用癫痫发作前复杂度的变化特征进行预测，具有较高地预测准确率：二是不仅对癫痫病人的颅内脑电信号适用，也适用于癫痫老鼠的颅内脑电信号；三是本发明的全部过程无需人工干预。需指出的是采用本方法预测癫痫发作，尽管预测时间近为几十秒的数量级，但对医生或患者及时采取有效预防措施有着重要的意义。

附图说明

图1癫痫脑电数据(a)发作前、(b)发作过程以及(c)发作后的频谱图。

图2一段癫痫脑电数据的复杂度。其中，垂直的实线表示发作时刻。(a)为一阶C₀复杂度曲线，(b)为二阶C₀复杂度曲线。

图3患者3颅内片状电极的具体位置，其致痫区位于前额左侧的新皮层。

图4癫痫患者1的二阶C₀复杂度变化曲线。其中，垂直的实线表示每次发作的时刻，虚线表示预测报警时刻。(a)为第一次发作，预测时间为64秒；(b)为第二次发作，预测时间为52秒；(c)为第三次发作，预测时间为44秒；(d)为第四次发作，预测时间为48秒，有一次虚警。

图5虚警示例。(a)癫痫患者8一次发作的二阶C₀复杂度曲线；(b)癫痫患者21一次发作的二阶C₀复杂度曲线。

具体实施方式

首先，我们给出C₀复杂度的具体计算公式：

定义一个长度为N的时间序列{x(t)，t＝0，1，2，…，N-1}，则其傅立叶变换序列为

$F_N\left(j\right)=\frac{1}{N}\underset{t=0}{\overset{N-1}{\mathrm{\Σ}}}x\left(t\right)W_N^{-\mathrm{tj}},j=\mathrm{0,1,2},\·\·\·,N-1---\left(3\right)$

其中，W_N＝e^-2πi/N。设{F_N(j)，j＝0，1，2，…，N-1}的均方值为

$G_N=\frac{1}{N}\underset{j=0}{\overset{N-1}{\mathrm{\Σ}}}{\left|F_N\left(j\right)\right|}^2---\left(4\right)$

令

对{

j＝0，1，2，…，N-1}作傅立叶逆变换

$\tilde{x}\left(t\right)=\underset{j=0}{\overset{N-1}{\mathrm{\Σ}}}{\tilde{F}}_N\left(j\right)W_N^{\mathrm{tj}},t=\mathrm{0,1,2},\·\·\·,N-1---\left(6\right)$

则C₀复杂度定义为

$C_0=\frac{\underset{t=0}{\overset{N-1}{\mathrm{\Σ}}}{|x\left(t\right)-\tilde{x}\left(t\right)|}^2}{\underset{t=0}{\overset{N-1}{\mathrm{\Σ}}}{\left|x\left(t\right)\right|}^2}---\left(7\right)$

对于任意一个时间序列{x(t)}，满足0≤C₀≤1。

其次，选取FDM数据库中一段癫痫病人发作过程的脑电信号进行研究。此段信号记录的是颅内脑电信号，采样率为256Hz，信号长度均为46min，选取致痫区的一个导联信号，对其做上文所述的带通滤波预处理后，按上述公式计算其C₀复杂度。附图2(a)为癫痫发作过程的一阶C₀复杂度变化曲线，可以看出，在发作前期，C₀复杂度有下降的趋势，而发作后，复杂度又升高。然而，这一变化趋势并不十分明显。我们对附图2(a)中计算得到的C₀复杂度曲线值再求一次复杂度，即得二阶C₀复杂度，如附图2(b)所示。从附图2(b)中，我们可以清楚地看到，发作前期的二阶C₀复杂度有明显的下降趋势，而在发作时刻出现很陡峭的下降，发作后又呈升高状态。基于这个实验结果，本发明提出了一种可用于临床预测癫痫发作的算法。该预测算法的基本思想是找出从发作间期到发作前期脑电数据的二阶C₀复杂度的变化特点，并识别出这种变化趋势，从而预测癫痫发作。

下面我们分别以癫痫患者和患有癫痫症老鼠的发作数据为例说明具体的实施方式。

1、癫痫病人发作数据

本例所采用的临床数据来自于FDM数据库。该数据是由德国Freiburg大学医学院癫痫病研究中心采集的术前颅内脑电数据。数据库记录了21个患有局部性癫痫症(FocalEpilepsy)病人的颅内脑电数据，采样率为256Hz，其临床特性在表1中列出。对于每个病人，有三个颅内电极位于或接近发作区，如附图3所示为患者3的电极位置及致痫区(Focal Region)和远离区(Remote Region)的选择。本文即选用致痫区导联的脑电信号进行分析。

表1癫痫患者脑电数据的临床特性

Patient	Sex	Age	Seizure type	H/NC	Origin	Electrodes	Seizuresanalyzed
								1	f	15	SP，CP	NC	Frontal	g，s	4
2	m	38	SP，CP，GTC	H	Temporal	d	3
								3	m	14	SP，CP	NC	Frontal	g，s	5
4	f	26	SP，CP，GTC	H	Temporal	d，g，s	5
								5	f	16	SP，CP，GTC	NC	Frontal	g，s	5
6	f	31	CP，GTC	H	Temporo/Occipital	d，g，s	3
								7	f	42	SP，CP，GTC	H	Temporal	d	3
8	f	32	SP，CP	NC	Frontal	g，s	2
								9	m	44	CP，GTC	NC	Temporo/Occipital	g，s	5
10	m	47	SP，CP，GTC	H	Temporal	d	5
								11	f	10	SP，CP，GTC	NC	Parietal	g，s	4
12	f	42	SP，CP，GTC	H	Temporal	d，g，s	4
								13	f	22	SP，CP，GTC	H	Temporo/Occipital	d，s	2

14	f	41	CP，GTC	H，NC	Fronto/Temporal	d，s	4
								15	m	31	SP，CP，GTC	H，NC	Temporal	d，s	4
16	f	50	SP，CP，GTC	H	Temporal	d，s	5
								17	m	28	SP，CP，GTC	NC	Temporal	s	5
18	f	25	SP，CP	NC	Frontal	s	5
								19	f	28	SP，CP，GTC	NC	Frontal	s	4
20	m	33	SP，CP，GTC	NC	Tempo/Parietal	d，g，s	5
								21	m	13	SP，CP	NC	Temporal	g，s	5

注：f-female，m-male；Seizure type：SP＝simple partial，CP＝complex partial，GTC＝generalizedtonic-clonic；Origin：H＝hippocampus，NC＝neocortical；Electrodes：g＝grid，s＝strip，d＝depth.

运用本发明方法对表1列出的21例病人的87次癫痫发作数据进行预测，结果如表2所示。表中TP(True Positive)表示正确预测次数；FP(False Positive)表示虚警次数，又称假阳性次数；MPT(Mean Prediction Time)表示平均预测时间；Sensitivity表示预测的正确率，且Sensitivity＝(TPs/Total Number of Seizures)×100％。

表2癫痫患者临床数据基于二阶C₀复杂度的预测结果

Patient	Number ofseizures	TPs	Sensitivity(％)	FPs	MPT(s)(mean±standarddeviation)
						1	4	4	100	1	52±8.6
2	3	3	100	1	44±10.4
						3	5	5	100	3	36±11.8
4	5	5	100	4	64±12.6
						5	5	5	100	2	41.6±6
6	3	3	100	1	45.3±9.4
						7	3	3	100	1	60±13.1
8	2	1	50	1	40±0
						9	5	5	100	2	79.2±32.2
10	5	5	100	3	52±31.2
						11	4	4	100	2	84±27.9
12	4	4	100	3	79±20.1
						13	2	2	100	2	84±40
14	4	3	75	0	53±5
						15	4	4	100	1	60±21.4
16	5	5	100	3	41.6±7.8

17	5	5	100	1	56±12.1
						18	5	4	80	3	87±33.4
19	4	4	100	2	56±9.4
						20	5	5	100	1	60±19.1
21	5	3	60	0	49.3±1.9
						Total	87	82	94.3	37	58.4±24.5

观察表2中记录的第一个病人的预测结果，该病人4次发作的预测结果如附图4所示，图中虚线表示报警预测的时刻，为了和发作时刻进行比较，用实线给出发作时刻(用SZ表示)。可以看出，该病人的4次发作均能可靠的预测出，故预测准确率为100％，但该组病人的第4次发作出现一次虚警，如附图4(d)所示，预测时间平均为52秒。对于21组数据87次癫痫发作，共有37次虚警，平均预测准确率为94.3％，平均预测时间58.4s，共漏警5次。

附图5给出了其中两次漏警的例子，其余3例与此类似。从图中我们可以看出，实线标出的发作时刻位于复杂度曲线最低点的前端，而如前文所述大多数实验结果显示，当癫痫发作时，复杂度应下降到最低点，这里我们可能不能完全排除数据记录的发作时刻存在某些误差。

2、老鼠癫痫发作数据

本例采用的老鼠癫痫发作数据为清华大学生物医学工程研究所提供的颅内脑电数据，采样频率为200Hz。我们用同样的方法对老鼠脑电数据进行预处理，即对数据做20～80Hz的带通滤波预处理后，计算其二阶C₀复杂度。除第一组数据外，每段数据发作前包含60～70分钟数据段(第一组为31分钟)，预测结果如表3所示。从表中可以看出，对于这四组老鼠癫痫发作的脑电数据，该算法的预测准确率为100％，仅第4组数据出现一次虚警，平均预测时间为51.6秒。

表3癫痫老鼠脑电数据基于二阶C₀复杂度的预测结果

Seizure number	FPs	Prediction time(s)	MPT(s)	Sensitivity(％)
					1	0	75
2	0	31.3
					3	0	56.3
4	1	43.8	51.6	100％

综上所述，本发明提出了基于二阶C₀复杂度的癫痫预测方法，该方法对于患有癫痫症的病人和老鼠的脑电数据都能达到比较理想的预测效果。尽管预测时间平均为几十秒的数量级，但对于癫痫患者及时地采取预防措施仍有着重要的意义，可以为预防和控制癫痫发作提供重要的临床应用价值。

Claims (3)

1、一种基于二阶C₀复杂度的癫痫发作预测方法，其特征在于具体步骤如下：

(一)脑电数据预处理：

病人的癫痫脑电数据x(t)作20～80Hz的带通滤波预处理，得到滤波后的信号

(二)C₀复杂度的计算：

(1)对预处理后的脑电信号

计算其C₀复杂度，得到一阶C₀复杂度序列

{Cf_i，i＝1，2，…，N}；

(2)计算复杂度序列{Cf_i，i＝1，2，…，N}的C₀复杂度，得到二阶C₀复杂度序列{Cs_i，i＝1，2，…，M)；

(三)预测算法

计算平均值m_j：

$m_j=\frac{1}{L}\underset{i=1}{\overset{L}{\mathrm{\Σ}}}{\mathrm{Cs}}_i,---\left(1\right)$

标准差d_j：

$d_j={\left(\frac{1}{L}\underset{i=1}{\overset{L}{\mathrm{\Σ}}}{({\mathrm{Cs}}_i-m_j)}^2\right)}^{1/2}.---\left(2\right)$

则得到两个序列{m_j，j＝1，2，…，K}和{d_j，j＝1，2，…，K}，其中 $K=[\frac{M-L}{D}+1];$ L为滑动窗口长度，D为每次移动长度；

若m_j-1-m_j＞2d_j-l，则发出报警，记录此时相应的二阶C₀复杂度值所对应的时刻。

2、根据权利要求l所述预测方法，其特征在于计算C₀复杂度时窗口长度取为1-2秒，每个复杂度值在时间轴上对应数据窗口的终点。

3、根据权利要求1所述预测方法，其特征在于滑动窗口长度L取为25-35个点，每次移动D取为10-20个点。

Images