癫痫是第二大神经系统疾病，表现为由大脑神经活动超同步引起的神经系统紊乱. 全世界约1%的人口患有癫痫^[1]，影响约5000万人^[2]. 现有的治疗方法主要是抗癫痫药物治疗和手术切除病灶. 没有靶向的药物会对神经系统产生较大副作用，其中约30%的患者无法通过药物进行控制. 癫痫发作时的突然抽搐和意识丧失使患者生活在无助和焦虑中，这些症状有可能导致更严重的伤害，如跌倒、溺水、烧伤甚至死亡. 如果能在癫痫发作前向护理人员及时发出警报，可以确保患者的安全，提高患者的生活质量. 脑电图具有较高的时间分辨率，可以反映各个脑区的状态变化，具有实时性好、便捷的特征，是监测神经活动的重要工具，已广泛应用于癫痫相关的医疗研究中.

癫痫发作预测算法一般包括特征提取和分类，所提取的特征涵盖线性、非线性和空间信息. 许多线性特征，如峰度^[3]、能量^[4]、统计特征^[5]和功率谱分析^[6–9]等相关研究证实脑电信号在发作前期和发作间期均有变化^[10–14]. Aarabi等^[3]研究表明，患者的脑电信号在发作前会出现峰度增加和方差减小.Teixeira等^[5]使用功率谱预测和检测癫痫，准确率达到75.8%. 此外，小波变换^[15]、自回归系数^[16]的分析方法在癫痫发作预测中得到广泛应用.非线性特征如李雅普诺夫指数^[17]和相关维数可以描述大脑的动态变化，被广泛应用于脑电信号分析^[18]. 研究证实，在癫痫发作前存在脑电活动由复杂向简单转化的状态^[19]. 在癫痫发作前数分钟至数十分钟脑电活动的李雅普诺夫指数、相关维数动力学特征指标显著减小，预示癫痫即将发作^[20].

现有研究表明，大脑各区域的超同步放电是癫痫的重要特征^[21]，基于空间信息的脑疾病研究取得了丰富的成果. Florian等^[22]发现患者发作前脑电信号相位同步有明显变化，可以用于预测癫痫发作. Kaspar等^[23]通过相关矩阵法研究60例癫痫患者的颅内多通道脑电图信号，分析癫痫发作过程中大脑不同区域之间的相互作用. 研究发现，在癫痫发作前各区域间相关性较弱，在发作时相关性显著增强. Williamson等^[24]采用时空相关结构分析预发作期的脑电信号，发现癫痫发作前大脑出现周期性放电.

随着机器学习的发展，分类器的使用使得癫痫发作的预测结果准确性有了显著提高. 针对患者个体训练预测模型是重大突破，可以获得更好的预测结果. 为了挑选出合适的特征子集，研究人员提出不同的优化方案，如最小冗余最大相关特征选择^[25]、全局冗余最小化、遗传算法^[26]等. Zhang等^[27]提取功率谱和功率谱比值，再通过主成分分析和分支定界法实现特征挑选. Birjandtalab等^[28]采用随机森林算法进行特征选择，结合t-分布邻域嵌入算法和功率谱密度特征实现癫痫检测.

综上所述，在癫痫发作预测方面已有较多成果. 然而，基于脑电图的癫痫发作预测还难以实现临床应用，还有一些挑战须解决，特别是：1）癫痫发作往往起源于某一区域并不断向周围传播，脑区之间的同步活动和空间连通性在传播过程中不断变化，现有基于多变量的方法没有考虑大脑内部的空间同步和连通信息. 2）癫痫发作预测是二分类问题，通常以发作前期和发作间期的特征为正、负样本进行训练. 所选择的发作前期正样本是否合理将直接影响分类器的性能，目前临床上没有明确方法判定最佳预发作时段.

本研究主要解决以下关键问题：1）通过计算脑功能连接网络参数来描述癫痫发作预测和检测中经常被忽略的空间信息和导联间的同步活动，再结合时频特征与非线性特征从多方面描述癫痫发作前的异常状态. 2）通过预发作数据段选取方法从发作前期自动选择出与发作间期相比具有最大线可分性的数据段作为正样本，并验证选取的合理性和有效性.

数据来源为MIT头皮脑电（sEEG）公共数据库^[29]和首都医科大学宣武医院患者脑电数据集. MIT公共数据库收集于波士顿儿童医院，由顽固性癫痫患者的脑电数据组成. 该数据采用国际10-20脑电图电极位置和命名系统进行记录，共收集22例（5名男性，3.0~22.0岁；17名女性，1.5~19.0岁），所有信号均以256 Hz进行采样. 该数据库主要针对多次发作的患者，其中一次发作用于训练模型，其余至少一次发作用于测试，共分析14名患者的68次癫痫发作. 首都医科大学宣武医院经志愿者允许后使用其脑电数据进行研究，并得到伦理委员会的批准. 该数据采用国际10-20脑电图系统进行记录，包含从19名患者中收集的94次癫痫发作，总时长为278 h. 脑电图在256 Hz下采集和记录，电极数量为26. 所有患者均包含一次以上的癫痫发作.

所提取的特征包括功率谱密度、功率谱熵和脑功能网络参数. 数据处理流程如图1所示.

特征提取自半重叠的2 s滑动窗. 综合考虑大脑状态的变化周期和频率分析的长度，将窗长设置为2 s，并采用0.5 Hz的高通滤波和47~53 Hz的带通滤波去除噪声干扰.

1）功率谱密度. 功率谱密度（power spectral density, PSD）表示信号在频域中的功率分布，用于描述癫痫发作引起的大脑状态变化. 脑电节律性活动一般分为5个频带，即 $ \omega $={δ, θ, α, β, γ}. 为了更细致地分析γ频段的信息，将其进一步分为4个子带：γ₁(30~50 Hz)、γ₂(50~70 Hz)、γ₃(70~90 Hz)、γ₄(90~128 Hz). 考虑了功率谱密度及其归一化值. 对于每个时间窗，第i个子带的功率谱表达式为

(1) $ G\left(f\right)=\frac{1}{{f}_{\mathrm{s}}N}{\left|\sum\nolimits _{n=0}^{N-1}{x}_{n}^{(c)}{\rm{{exp}}}\;(-{\rm j}2\text{π} fn)\right|}^{2};\;-\frac{{f}_{\rm s}}{2}<f<\frac{{f}_{\rm s}}{2}. $

(2) $ {p}_{i}=\sum\nolimits _{f={f}_{1}}^{{f}_{2}}G\left(f\right);\; i=\mathrm{1,2},\cdots ,8. $

式中： $ {x}_{n}^{(c)} $为第c导联的第n个采样点，N为信号的长度， $ {f}_{\mathrm{s}} $为采样频率， $ {f}_{1} $、 $ {f}_{2} $为第i个子带对应的频率范围.

归一化功率谱为各子带的功率与总功率之比，计算方法为该子带的功率除以总功率P_tot：

(3) $ {p}_{{\mathrm{n}\mathrm{o}\mathrm{r}\mathrm{m}}_{i}}=\frac{{p}_{i}}{{p}_{\mathrm{t}\mathrm{o}\mathrm{t}}}=\frac{\displaystyle\sum\nolimits _{f={f}_{1}}^{{f}_{2}}G\left(f\right)}{\displaystyle\sum\nolimits _{f=1}^{{{f}_{\mathrm{s}}}/{2}}G\left(f\right)};\;i=\mathrm{1,2},\cdots ,8. $

2）功率谱熵. 功率谱熵（power spectral entropy，PSE）能较好地描述癫痫发作的频谱分布，可以理解为测量频域中的时间不确定性和复杂性. 第i子带的功率谱概率密度分布函数和功率谱熵的表达式分别为

(4) $ P\left(f\right)={G\left(f\right)}\Bigg/{\sum\nolimits _{f=1}^{{{f}_{\mathrm{s}}}/{2}}G\left(f\right)}. $

(5) $ {p}_{{\mathrm{e}\mathrm{n}}_{i}}=-\sum\nolimits _{f={f}_{1}}^{{f}_{2}}P\left(f\right){\rm{ln}}\;P\left(f\right);\;i=\mathrm{1,2},\cdots ,8. $

3）脑功能网络参数. 神经元同步活动是癫痫发作的特征之一^[21]. 通过对不同导联记录的脑电图信号进行时间相关性分析来分析脑功能的连通性，可以反映不同区域之间的统计连通性关系^[30].

脑功能连接网络参数（functional connectivity network，FCN）分析方法包括2个步骤：基于皮尔逊相关系数建立相关矩阵并对其二值化；采用平均路径长度和聚类系数来度量网络特性. 相关矩阵表达式如下:

(6) $ {\rho }_{i,j}=\frac{{\rm{cov}}\;\left({x}_{i},{x}_{j}\right)}{{\sigma }_{{x}_{i}}{\sigma }_{{x}_{j}}};\;i,j=1,\cdots ,c\;. $

式中：c为导联的个数； $x_i $为第 $i $导联的信号； ${\sigma}_{x_i} $为第 $i $导联信号的标准差；cov为协方差，表达式如下：

(7) $ {\rm{cov}}\;\left( {{x_i},{x_j}} \right) = E\left[ {\left( {{x_i} - {\mu _{{x_i}}}} \right)\left( {{x_j} - {\mu _{{x_j}}}} \right)} \right]. $

式中： $\mu_{x_i} $为第i导联信号的均值.

通过确定的阈值将相关矩阵转换为二值脑网络矩阵，小于阈值的矩阵值设为0，表示2个导联之间没有连接；大于阈值的矩阵值设为1，表示2个导联之间有连接. 阈值计算方法为单样本t检验，具体步骤如下.

1）设相关矩阵二值化的阈值为T，样本均值为μ，在显著性水平α下，假设： $ {H}_{0}:\mu <T,{H}_{1}:\mu \geqslant T $.

$ {H}_{0} $拒绝域为

(8) $ t = \left( {\bar X - T} \right)/\sqrt {{S^2}/n} \geqslant{t_\alpha }\left( {n - 1} \right). $

式中： $ \bar X $为样本均值，即对相关矩阵的所有值求均值； $ {S}^{2} $为样本方差，即对相关矩阵所有值求方差；n为样本数，即导联数的平方；置信度α=0.95. 以每个窗所计算的相关矩阵为一组样本，计算T，得到每个窗的阈值.

2）为了使阈值划定更具有统计意义，对于一位病患，挑选当前时刻前120个窗所求得的T进行平均，作为当前时刻的脑功能网络的划分阈值. 每移动一个窗，划分阈值更新一次.

同一位患者在不同时期所构建的二值脑网络连接矩阵如图2所示. 图中，黑色代表0，表示2个导联之间没有连接；白色代表1，表示2个导联之间有连接. 每位患者的大脑功能连接性都存在明显差异.

平均路径长度为第i个导联与其他导联间的平均距离，用于衡量网络传递信息的能力：

(9) $ {L}_{i}=\frac{1}{N}\sum\limits _{j\in V}{d}_{ij}. $

式中：V为所有导联的序号集合， $ {d}_{ij} $为第i导联到第j导联的路径长度.

聚类系数表示一个节点的所有邻居之间互为邻居的可能性，用于衡量节点趋向于聚集在一起的程度. 第i导联的聚类系数表达式为

(10) $ {C}_{i}=\frac{2{e}_{i}}{{k}_{i}\left({k}_{i}-1\right)}. $

式中： $ {e}_{i} $为第i导联的邻居们相互连接的边数， $ {k}_{i}({k}_{i}-1) $为连接边数的上限.

综上所述，每2 s从一个导联中提取26个特征，包括8个绝对功率谱值、8个归一化功率谱值、8个功率谱熵值以及脑功能网络参数的平均路径长度和聚类系数. 与仅计算功率谱密度相比，使用这3类特征的关键优势在于：这些特征描述了癫痫发作的3个不同维度，包括功率在频率域的分布、非线性动力学信息和大脑不同区域间的同步活动. 不同特征之间的相关性较小，可以得到更多的信息.

癫痫脑电数据可以分为几个时段：发作间期、发作前期、发作期和发作后期. 癫痫发作前脑电信号发生显著变化的时段为发作前期，可以持续几分钟到几个小时，医学领域对这一时段的时间长度并无统一划分. Teixeira等^[5]分别研究发作前期的时间长度选择10、20 、30、40 min时的预测效果，发现发作前期时间越长，预测的准确性越高. 按照癫痫脑电分析研究的常用时段划分，定义癫痫发作时刻前60 min为发作前期（pre-ictal period）；临床发作时间标定的癫痫起始、终止时间之内为发作期（ictal period）；癫痫终止时刻后30 min为发作后期（post-ictal period），介于上一次发作后期与此次发作前期之间的时段为发作间期（inter-ictal period）.

本研究提出预发作数据段选取（pre-ictal period selection，PPS）方法，以最小类内距离和最大类间距离为标准，从发作前期的数据中选取连续的10 min作为预发作期. 该数据与从发作间期随机挑选的10 min数据相比，其特征具有最大线性可分性，通过窗长为10 min，包含1 s重叠的滑动窗进行挑选. 该方法使用类内距离 $ {S}_{\rm{w}} $和类间距离 $ {S}_{\rm{{b}}} $来衡量，通过保持类内距离最小和类间距离最大来寻找最具判别性的预发作数据段.

(11) $ {S}_{\rm{{w}}}=\sum\limits _{{ x}\in {X}_{i}}{\left({{{x}}-{{\mu}} }_{i}\right)}^{\rm T}\left({{{x}}-{{\mu}} }_{i}\right), $

(12) ${S}_{\rm{{b}}}=\sum\limits _{i=\mathrm{0,1}}{\left({{\mu }}_{i}-{{\mu}} \right)}^{\rm T}\left({{\mu }}_{i}-{{\mu}} \right). $

式中： $ {x } $为m维特征向量， $ {x }=[{x}_{1},{x }_{2},\cdots ,{x}_{m}{]}^{\rm{{T} }}$；X₀、X₁分别表示负样本和正样本集合， $ {{\mu }}_{1} $、 $ {{\mu }}_{0} $、 $ {{\mu}} $均为m维向量，分别为正样本均值、负样本均值和所有样本均值. 预发作数据段可以认为是与发作间期相比类间距离最大、类内距离最小的发作前期数据段. 因此，挑选准则表达式为

(13) $ J={\left|{S}_{\rm{{w}}}+{S}_{\rm{{b}}}\right|}/{\left|{S}_{\rm{{w}}}\right|}. $

最终选取J最大的10 min作为预发作期进行特征选择和分类. 理论上所选预发作数据段为发作前期内与发作间期相比线性可分性最好的时段，但有时会受到噪声的影响，因此设置J平均值的5倍作为上限以避免噪声.

不同患者的癫痫类型和病灶位置不同，其信号特征和有效电极也不同. 算法^[4,6]仅考虑癫痫发作的不同时期而没有考虑个体差异，因而效果不佳. 为了实现针对病患个人进行预测，特征挑选和导联挑选得到广泛应用.

首先使用弹性网^[31]对每一导联分别进行特征挑选. 弹性网表达式如下：

(14) $ {\rm{min}}\;\mathop \sum \nolimits_{i = 1}^m \left( {{y_i} - {{{w}}^{\rm{T}}}{{{x}}_i}} \right) + {\lambda _2}\left\| {{w}} \right\|_2^2 + {\lambda _1}{\left\| {{w}} \right\|_1}, $

(15) $ \left\| {{w}} \right\|_2^2 = \mathop \sum \nolimits_{j = 1}^m w_j^2, $

(16) $ {\left\| {{w}} \right\|_1} = \mathop \sum \nolimits_{j = 1}^m \left| {{w_j}} \right|. $

式中： ${({{ x}_i},{{ y}_i})}$为一个样本； ${{ x}_i} $为特征矩阵， $ {{{x}}_i} = \left[ {{x_{i1}},} \right.$ $\left. {{x_{i2}}, \cdots ,{x_{id}}} \right]$； $y_i $为响应；w为投影矩阵； $ {\lambda}_1 $、 $ {\lambda}_2$为正则化参数.

弹性网是LASSO和岭回归的一个凸组合，当 $ {\lambda }_{2} $=0时，式(14)等价于LASSO回归；当 $ {\lambda }_{1} $=0时，式(14）等价于岭回归. 为了降低过拟合风险，同时获得较好的特征挑选结果，设置 $ {\lambda }_{1}={\lambda }_{2}=0.5 $. 在进行特征挑选时，使用10折交叉验证使得特征挑选更加合理.

根据所选特征进行导联挑选，以避免无效和冗余特征的影响. 从所有导联中选择K个有效的导联，包括子集搜索和子集评价过程. 定义导联集 $ {S}_{1} $存放剩余导联，初始化为包含所有导联的集合；定义导联集 $ {S}_{0} $存放被挑选导联，初始化为空集. 导联挑选的具体步骤如下：1）依次对 $ {S}_{1} $中导联的有效特征计算目标函数J，如式(13)所示. 选择J最大的导联放入 $ {S}_{0} $中，同时将该导联在 $ {S}_{1} $中剔除；2）将 $ {S}_{1} $中每一导联的有效特征依次与 $ {S}_{0} $中所有导联的有效特征组合后计算目标函数J，选择J最大的导联放入 $ {S}_{0} $中，并在 $ {S}_{1} $中剔除；3）循环 $ K $次，最终 $ {S}_{0} $中包含 $ K $个导联，为最终挑选导联. 本研究定义K=4.

为了实现针对病患个人进行预测，患者的发作次数应多于一次，其中1次发作用于训练，其余至少1次发作用于测试. 以挑选出的预发作期的特征为正样本，发作间期的特征为负样本，输入至XGBoost分类器中. 认为癫痫发作前60 min内的警报是真阳性，发作后30 min内的数据因可能存在发作后脑电图抑制现象而不对其进行预测.

采用击中率（recall，R）和假阳性率（false positive rate，FPR）衡量算法性能. 击中率是准确预测的癫痫发作次数与癫痫发作总次数之比，反映算法的查全程度；假阳性率是每小时内发生的虚警个数与所有未发作样本数之比，反映算法是否容易产生虚警. 对每位患者分别训练模型，采用交叉验证方法，每次实验依次以其中一次发作的数据训练模型，再对所有次发作进行预测，计算该次试验的击中次数. 对多次试验的结果求平均值，作为平均击中次数，并计算击中率. 算法在MIT公共数据库和宣武医院数据集上进行测试，在MIT数据库上的击中率为95.76%，假阳性率为0.1073 h⁻¹；在宣武医院数据集上的击中率为97.80%，假阳性率为0.0453 h⁻¹. MIT数据库与宣武医院数据集的预测结果如表1、2所示. 表中，T为数据时长，N_e为有效发作次数，n_h为平均击中次数.

提出新的预发作数据段选取方法，对发作前期的数据进行挑选，使挑选出的样本具有更强的判别性，从而提高预测的准确性. 为了更好地说明预发作数据段选取的挑选效果，对所选预发作数据段和发作前时刻的特征进行对比，如图3所示，依次展示了MIT数据库中第10、14、18、22号患者的正、负样本对比结果. 图中，第1列为J曲线，其中第1个圆圈表示J的最大值，在时间轴上对应于挑选的10 min数据；第2个圆圈表示J的最后一个值，在时间轴上对应于最接近发作的10 min数据. 第2列为以挑选的10 min数据作为正样本与发作间期负样本特征进行对比，第3列为以最接近发作的10 min数据作为正样本与相同的发作间期负样本特征进行对比. 图中，“黑色*”表示发作前期的特征，“灰色○”表示发作间期的特征，所有特征均进行了归一化.

比较图3的第2、3列可知，相比于以发作前时刻作为正样本，若以挑选后的预发作数据段作为正样本，正样本与发作间期负样本的特征重叠更少，可分离程度更大，特征具有更好的线性可分离性. 观察J的曲线可知，从发作前1 h到发作时刻，J不断变化，最大值并不完全出现于发作前时刻，说明发作间期数据与发作前期数据段的可分离程度不一定会随着癫痫临近而增加.

将不同组的数据输入XGBoost分类器中进行预测，分析预发作数据段选取对分类器性能的影响. 实验组为所有患者均进行挑选，输入所挑选的预发作期（有PPS）和发作间期作为正、负样本；对照组为所有患者均不进行挑选，选择发作前时刻（无PPS）和相同的发作间期作为正、负样本. 算法在MIT数据库上进行验证，验证结果如图4、5所示. 经预发作数据段选取后，所有患者的平均击中率由90.28%提升至95.76%，假阳性率由0.1486 h⁻¹降低至0.1073 h⁻¹. 击中率和假阳性率均得到改善，分类器性能有所提升.

提取脑功能网络参数以表征大脑内部的空间信息和导联间的同步活动. 为了验证脑功能网络参数的有效性，将不同组合的数据输入XGBoost分类器中进行预测. 在相同的预发作数据段选取、特征挑选和导联挑选的情况下，实验组使用含有脑功能网络参数的26个特征（有FCN），对照组使用不含有脑功能网络参数的24个特征（无FCN）. 算法在MIT数据库上进行验证，验证结果如图6、7所示.

使用包含脑功能网络参数的特征进行预测，所有患者的平均击中率由90.22%提升至95.76%，假阳性率由0.1261 h⁻¹降低至0.1073 h⁻¹. 引入脑功能网络参数可以分析癫痫发作前的空间信息，表征空间状态变化，更加符合癫痫发作时的大脑生理模型，从而提高预测的准确性.

目前已有许多癫痫发作预测的方法. 如表3所示为2013—2017年的6篇癫痫发作预测相关文献，对比了不同方法在MIT数据库上的结果，这些研究均针对患者个体训练模型. 与现有其他方法相比，本研究的方法具有更高的击中率和更低的假阳性率. 可以得出，本研究算法的预测效果更好.

提出基于患者特异性的癫痫发作预测算法，并提出新的预发作数据段选取方法来选择最具判别性的预发作期时段. 提取脑电信号的时频特征和空间特征，经过预发作数据段选取、特征挑选和导联挑选后对模型进行训练，最终在MIT数据库和宣武医院数据集上的击中率分别为95.76%、97.80%，假阳性率分别为0.1073、0.0453 h⁻¹. 与现有研究相比，本研究的特征提取更全面，预发作数据段选取方法更合理，预测模型的击中率高，假阳性率低，具有较好的预测效果.

使用头皮脑电图进行癫痫发作预测已有较多研究成果，与颅内脑电图（intracranial EEG, iEEG）相比，头皮脑电图包含更多的噪声. iEEG的采样率更高，包含更多的大脑模式信息. 随着神经电刺激技术的发展，基于iEEG实现癫痫发作预测可以与药物释放或电刺激治疗相结合，形成闭环控制，成为癫痫治疗的新方法. 未来将针对iEEG验证该方法的有效性.