癫痫是由多种原因导致的慢性脑部疾患，以脑神经元异常放电引起反复、发作性和短暂性的中枢神经系统功能失常为特征.癫痫发作具有突发性，无明显的先兆，是癫痫病人致残和死亡的重要原因^[1].目前对癫痫脑电信号的研究方向主要集中在检测、识别和预测.癫痫发作的预测能够大幅度增强癫痫治疗的效果和提高癫痫患者的生活质量.然而癫痫发作的准确预测是一个有巨大临床价值但仍然没有被解决的问题.

近十多年来,人们从线性分析^[1]、非线性动力学分析^{[2, 3]}、相位同步分析^[4]、时域和频域分析^{[5, 6]}等多种方法分析脑电，希望得到一些特征参量，准确预测癫痫发作.文献[1]采用线性模型方法分析脑电信号，结果表明在癫痫发作前1~6 s有明显变化.文献[3]采用关联维数法分析脑电信号，结果表明在发作前几分钟关联维数显著下降. 文献[4]利用希尔伯特黄变换提取瞬时相位，通过分析相位的变化对癫痫信号进行预测.然而上述方法都是基于单变量或双变量脑电信号分析，没有反映出大脑各区域间的相互作用，所以结果距离临床应用还有一定的距离.文献[5]利用加权频率作为特征值进行癫痫的发作预测，而脑电信号在频域的变化是非常复杂的，仅用单一特征值进行分类缺乏应用的普遍性.文献[6]利用经验模式分解方法提取出幅值和频率的调制带宽作为特征值进行分类.虽然该特征值包含了频域的大部分信息，但是却丢失了重要的时间信息.

本研究针对现有预测方法中信号特征缺乏空间信息和变化方向的问题，考虑到脑电的非线性、非平稳随机特征及癫痫发作时脑电图信号频谱和频率随时间的明显变化,提出将希尔伯特边际谱、希尔伯特加权频率和希尔伯特边际谱的变化方向相结合，组成多变量希尔伯特频域模型,在空间域、时域和频域范围内探索癫痫发作预测的新方法，取得较理想的预测效果.

希尔伯特黄变换是一种新的非线性、非平稳信号分析方法.其核心部分是经验模态分解方法(empirical mode decomposition，EMD).任何复杂信号都可以由EMD分解成有限个本征模态函数(intrinsic mode function,IMF)，再利用希尔伯特变换获得信号的时频分布^{[7, 8]}.设信号为x(t)，Hilbert变换公式为

通过式(2)提取出信号的瞬时幅度

和瞬时相位

根据式(1)求出希尔伯特变换，就可以利用式(5)对时间积分得到希尔伯特边际谱(Hilbert marginal spectrum,HMS)h(ω)：

边际谱表征了每一个频率值所对应的总幅度值，在统计意义上，表征了整个时间跨度内信号在每个频率点上能量累积的分布情况^[9].

对于非平稳信号而言，频率是随时间变化的.因此需要另外一个频率概念——瞬时频率.根据式(4)求出瞬时相位，则瞬时频率公式为式(6).利用式(7)求出希尔伯特加权频率(Hilbert weighted frequency,HWF)作为特征参量之一.f为式(6) 求出的瞬时频率，a为式(3)求出的瞬时幅值^{[7, 8]}, t=1,…,L 为采样点数.

Palus^[10]提出一种基于香农熵的条件互信息论方法，计算两个变量的变化方向指数.联合条件的互信息定义为

SCDI(x,y)的值域为[-1,1].即当SCDI(x,y)>0时,表示导联x驱动导联y；而当SCDI(x,y)<0时 ，表示导联x被导联y驱动；当SCDI(x,y)=0时，导联x和导联y互相影响对称相等.

支持相量机是一种基于监督学习的分类问题工具.在解决小样本、非线性及高维模式识别问题时具有良好的性能.SVM的基本原理是在高维特征空间中寻找一个超平面，使得这个超平面能将输入数据以最大的相互距离隔开.对于非线性问题，可以利用某种非线性变换将非线性问题转化为某个高维空间的线性问题.综合泛函数理论，在最优分类面的求解中加入核函数，就可以实现非线性问题的线性分类.本文采用径向基核函数.

本文采用的癫痫脑电数据采自波士顿儿童医院(CHB-MIT脑电数据库)，18位受试者均为药物难以控制的儿童顽固性癫痫发作患者(9名男孩，9名女孩，年龄在3~18岁).采样频率为256 Hz，电极放置方式采用国际10~20系统法.每个受试者均采集18个导联的脑电信号.每个导联取三组脑电信号，分别为癫痫发作间期组(A)、癫痫发作前期组(发作前5~45 min)(B)、癫痫发作期组(C).脑电数据以段为单位，每段有1 000个数据点，每组取10段脑电数据进行分析.每位受试者的脑电检测时间均超过8 h.

首先利用小波变换对所有导联的三组脑电数据进行处理，提取4个典型的频率成份：α波(8~13 Hz),β波(14~30 Hz),θ波(4~7 Hz),δ波(1~3 Hz).第二步对已得到各子频带做希尔伯特黄变换，分别利用式(5)和式(7)求出各子频带的HMS和HWF，绘制出HMS和HWF的脑电地形图.图 1是2号患者δ波的HMS脑电地形图(随机抽取每组的三段脑电地形图).由图 1可以看出，在癫痫发作间期，HMS在额区和颞区都有出现，但强度较低.在癫痫发作前期，HMS逐渐向枕区进行扩散，强度逐渐升高.在发作状态组，HMS强度明显增强，主要集中在枕区.表 1是根据Kruskal-Wallis检验，对A，B，C三组之间的HMS值作非参数方差，分析P值的平均值.由表 1可以看出，δ波和θ波的HMS值在A，B，C三组之间相比，P值都小于0.05.以α=0.05为检验准则，发作间期、发作前期和发作期δ波和θ波的HMS值有统计意义，即有明显差异.所以选取δ波和θ波的HMS值作为特征值进行癫痫预测.

图 1(Figure 1)

Fig. 1

图 1 δ波的HMS脑电地形图 Fig. 1 Topography of HMS for δ sub-band (a)—A组; (b)—B组; (c)—C组.

表 1

Table 1

表 1 (Table 1)

表 1 HMS和HWF Kruskal-Wallis检验的P值 Table 1 P of HMS and HWF Kruskal-Wallis test 子频率 Kruskal-Wallis检验 A~B B~C A~C δ-HMS 0.003 0.005 0.008 δ-HWF 0.014 0.009 0.016 θ-HMS 0.012 0.009 0.018 θ-HWF 0.005 0.012 0.009 α-HMS 0.063 0.113 0.092 α-HWF 0.095 0.082 0.148 β-HMS 0.073 0.051 0.064 β-HWF 0.184 0.137 0.192

表 1 HMS和HWF Kruskal-Wallis检验的P值 Table 1 P of HMS and HWF Kruskal-Wallis test

图 2是2号患者δ波的HWF脑电地形图(随机抽取每组的三段脑电地形图).由图 2可以看出，在癫痫发作间期，HWF在整个大脑区域的取值范围在0.3~0.5之间，在枕区和颞区取值较高.在癫痫发作前期，HWF的值逐渐下降，取值范围在0.25~0.35之间，越临近癫痫发作期，HWF的值越低.在癫痫发作期，HWF的值明显降低，取值范围在0.1~0.25之间.由表 1可以看出，δ波和θ波的HWF的P值都小于0.05，以α=0.05检验准则，δ波和θ波的HWF值有统计意义，即有明显差异.所以选取δ波和θ波的HWF值作为特征值进行癫痫预测.

图 2(Figure 2)

Fig. 2

图 2 δ波的HWF脑电地形图 Fig. 2 Topography of HWF for δ sub-band (a)—A组; (b)—B组; (c)—C组.

根据前面的分析，提取δ波和θ波的希尔伯特边际谱同步方向指数.

图 3是2号患者第一导联和第五导联的δ波HMS耦合方向的曲线图.该患者在发作间期时，SCDI指数趋近于0，说明两个导联之间HMS的变化处于一种平衡状态.在发作前期时，SCDI指数是在0值附近有比较明显的上下波动，说明两个导联之间的平衡状态被打破.在发作期，SCDI指数大于0，说明第一导联对第五导联的作用增强，所以耦合方向是从第一导联到第五导联.利用上述方法对每位患者各导联之间SCDI指数进行分析，将SCDI指数作为预测特征值.

图 3(Figure 3)

Fig. 3

图 3 第一导联和第五导联δ波的相位耦合方向曲线图 Fig. 3 Phase coupling direction graph for δ between first and fifth route

最后将HMS，HWF和SCDI组成的特征矢量作为多变量希尔伯特频域模型(MHSM)，输入到SVM中，从而实现癫痫发作的预测.

本文采用癫痫发作预测特征方法对预测结果进行评估.癫痫发作预测特征由4个参量组成：癫痫预测范围 (seizure prediction horizon ，SPH)、癫痫发作周期(seizure occurrence period ，SOP)、最大错误预报率(maximum false prediction rate，FPRmax)、癫痫预测敏感度(sensitivity of seizure prediction,SSP)^[4].癫痫预测范围是指从癫痫预警到癫痫发作的时间间隔.SPH应该足够长，目的是使癫痫患者在收到警告之后能够有充足的时间采取一些措施避免癫痫的发生.癫痫发作周期是指癫痫可能发作的时间间隔.SOP应该足够短，目的是缩短患者的等待时间，降低患者的焦虑程度.如果在癫痫发作周期中没有癫痫发作发生，那么这个预警就是错误的预测.在实际的预测中，错误预测是不可避免的.最大错误预报率是指单位时间内可以产生错误预报的最大值.癫痫预测敏感度是指正确预报的癫痫数和癫痫总数的比率.

本研究总共分析了4种方法的预测效果：δ波的MHSM预测、θ波的MHSM预测、原始EEG的HMS值作为特征值的预测^[9]、原始EEG的HWF值作为特征值的预测^[5].图 4是以SPH为参考值的平均敏感度变化趋势图.

图 4(Figure 4)

Fig. 4

图 4 基于SPH的敏感度分布图 Fig. 4 Dependence of averaged sensitivities on SPH

由图 4可以看出，当SPH值超过30 min后，δ波和θ波的MHSM(multivariate HSM)预测的敏感度明显升高，HMS和HWF的敏感度没有明显变化.这个结果表明，δ波和θ波的MHSM预测方法可以得到较大的SPH值，并且预测敏感度也较高.图 5是以SOP为参考值的平均敏感度变化趋势图.当SOP值在5~10 min时，δ波和θ波的MHSM预测的敏感度明显上升；当SOP值大于10 min时，δ波θ波的MHSM预测的敏感度变化比较缓慢，并且敏感值接近90%；当SOP值大于25 min时，HMS和HWF的敏感度变化缓慢，并且敏感值在50%~60%之间. 这个结果表明，δ波和θ波的MHSM预测方法可以得到较短的SOP值，并且预测敏感度也较高.图 6是以FPRmax为参考值的平均敏感度变化趋势图.对于相同的FPRmax值，δ波和θ波的MHSM预测的预测敏感度明显高于HMS和HWF的预测敏感度.对于相同的预测敏感度，δ波和θ波的MHSM预测的错误预报率最低.这个结果表明，δ波和θ波的MHSM预测方法有较低的错误预报率和较高的预测敏感度.

图 5(Figure 5)

Fig. 5

图 5 基于SOP的敏感度分布图 Fig. 5 Dependence of averaged sensitivities on SOP

图 6(Figure 6)

Fig. 6

图 6 基于FPRmax的敏感度分布图 Fig. 6 Dependence of averaged sensitivities on FPRmax

本文提出多变量希尔伯特频域模型方法，对低频癫痫脑电信号进行发作预测，采用癫痫发作预测特征方法对预测结果进行评估.实验结果表明：癫痫预测范围大于30 min时，患者有足够的时间采取措施应对癫痫发作；癫痫发作周期小于10 min，缩短患者等待时间，降低焦虑程度，减轻痛苦；与目前已有的相关方法相比，该方法具有较低的最大错误预报率和较高的预测敏感度.