癫痫脑电的自动分类对于癫痫的诊断和治疗具有重要意义。本文提出了一种基于小波多尺度分析和极限学习机的癫痫脑电分类方法。首先，利用小波多尺度分析对原始脑电信号进行多尺度分解，提取出不同频段的脑电信号。然后采用Hurst指数和样本熵两种非线性方法对原始脑电信号和小波多尺度分解得到的不同频段脑电信号进行特征提取。最后，将得到的特征向量输入到极限学习机中，实现癫痫脑电分类的目的。本文采用的方法在区分癫痫发作期和发作间期时取得了99.5%的分类准确率。结果表明，本方法在癫痫的诊断和治疗中具有很好的应用前景。
引用本文： 崔刚强, 夏良斌, 梁建峰, 涂敏. 基于小波多尺度分析和极限学习机的癫痫脑电分类算法. 生物医学工程学杂志, 2016, 33(6): 1025-1030,1038. doi: 10.7507/1001-5515.20160165 复制
引言
癫痫是一种由多种病因引起的慢性脑功能障碍疾病，其典型特征是大脑局部神经元反复地、突然地过度放电，导致中枢神经系统功能失常[1]。全世界有1%~2%的人口患有癫痫[2]。癫痫患者临床表现为肌肉抽搐、意识丧失等，频繁的癫痫发作往往会给患者的身体带来极大的伤害，甚至危及生命，而且给患者家庭和社会带来沉重的负担。因此，加强癫痫疾病的前期诊断和后期治疗工作，具有极其重要的意义。
脑电在临床癫痫疾病的诊断、病灶定位、疗效评估中有着不可替代的作用。癫痫发作表现为高幅同步节律波，包括棘波、尖波、棘慢复合波和尖慢复合波等[3]。目前，癫痫疾病的诊断工作主要由经验丰富的医生肉眼观察长时间的脑电信号来完成，工作量巨大，而且受主观因素影响较大，不同医生给出的诊断结果有可能并不相同。因此，有必要开发一种癫痫脑电自动分类的方法。
癫痫脑电分类任务可以分为特征提取和分类两部分工作。非线性动力学相关方法已经在癫痫脑电信号的非线性特征提取工作中有了很好的应用。已有研究者将Lyapunov指数[4]、Hurst指数[5]、近似熵[6]、分形维数[7]、样本熵[8]、功率谱熵[9]、Lempel-Ziv复杂度[8]等非线性方法用于癫痫脑电信号的特征提取中。目前，癫痫脑电分类中使用较多的是K-最邻近分类器[10]、贝叶斯分类器[11]、决策树[12]、人工神经网络(artificial neural network,ANN)[11]、支持向量机(support vector machine,SVM)[13-15]等，并且已经取得了一定的成果。文献[14]比较了SVM、贝叶斯分类器、线性判别分析、K-最邻近分类器、ANN等几种分类方法在同一数据集上的分类效果，最终ANN取得了最高的97.7%的分类准确率。
Huang等[16]提出了一种新颖的单隐层前馈神经网络的机器学习算法，即极限学习机(extreme learning machine,ELM)，在学习机训练开始时，就随机给定它的隐藏层偏置和输入权值，且无需在训练过程中调整。与反向传播(back propagation,BP)神经网络相比，极限学习机既防止了学习过程中陷入局部最优解，又具有泛化能力强、训练速度快的优点，有较好的应用前景。
本文提出了一种基于小波多尺度分析和极限学习机的癫痫脑电分类算法。首先，对原始脑电信号进行小波多尺度分解，提取出不同频带的脑电信号；然后，利用Hurst指数和样本熵两种方法进行非线性特征提取；最后，应用极限学习机分类器对癫痫脑电进行分类，并在波恩大学的癫痫脑电数据库上测试分类器的性能。
1 基于小波多尺度分析的特征提取方法
1.1 小波多尺度分析
小波多尺度分析是一种将信号分解到多个尺度上的子空间的方法，分解后的信号在各个子空间同时具有时域和频域的分辨率，可以方便地分析信号的时频特性，在生物医学信号处理中有很好的应用。
对于信号S来说，就是将其分解为尺度空间的A1近似部分和细节空间的D1细节部分，然后对A1部分再次进行尺度和细节的分解，得到A2和D2两部分，再对A2进行如此分解。这里的近似信号对应于信号的低频部分，而细节信号对应于信号的高频部分。与傅里叶分析不同的是，高频部分是分层次的，是在不同分辨率下逐步产生的。图 1所示为信号S经过小波3层分解的示意图。
图1
小波3层分解示意图
Figure1.
Scheme of 3 layer wavelet decomposition
1.2 Hurst指数
Hurst指数是英国科学家H. E. Hurst在20世纪中叶提出的一种判断时间序列是否具有时间依赖性的参数，现在也用于研究混沌时间序列的预测问题[17]。本文采用R/S法，即重标极差方法(rescaled range analysis)。步骤如下：
对于给定的时间序列xi(i=1,2,…,N)，假定计算长度n，将时间序列划分为长度为n的M个子序列，第m个子序列的第t个数据记为xt,m(t=1,…,n,m=1,…,M) ，第m个子序列的重标极差比为
| ${{\left( R/S \right)}_{n,m}}=\frac{{{R}_{n,m}}}{{{S}_{n,m}}}$ |
式中，Rn,m为第m个子序列的极差，Sn,m为第m个子序列的标准差：
| ${{R}_{n,m}}=\underset{1\le k\le n}{\mathop{\max }}\,\sum\limits_{t=1}^{k}{k({{x}_{t,m}}-\overline{{{x}_{n,m}}})}-\underset{1\le k\le n}{\mathop{\min }}\,\sum\limits_{t=1}^{k}{k({{x}_{t,m}}-\overline{{{x}_{n,m}}})}$ |
| ${{S}_{n,m}}=\sqrt{\frac{1}{n}\sum\limits_{t=1}^{k}{{{({{x}_{t,m}}-\overline{{{x}_{n,m}}})}^{2}}}}$ |
式中，xn,m为第m个子序列。时间序列xi对应于计算长度n的重标极差为
| ${{\left( R/S \right)}_{n}}=\frac{1}{M}\sum\limits_{m=1}^{M}{{{\left( R/S \right)}_{n,m}}}$ |
改变计算长度n，重新计算重标极差比，可以得到一组计算长度与重标极差的序列。Hurst认为，重标极差序列(R/S)n与计算长度n之间存在指数关系如下：
| ${{\left( R/S \right)}_{n}}=C\times {{n}^{H}}$ |
式中，C为常数，H即为Hurst指数。通过取双对数：
| $\lg \left( R/S \right)=\lg \left( C \right)+H\lg \left( n \right)$ |
采用最小二乘法拟合公式(6)，可以估算出Hurst指数。
1.3 样本熵
Richman等[18]提出了一种评价给定时间序列复杂度的新方法，即样本熵(sample entropy,SamEn)，它是近似熵的一种改进方法。它是条件概率(conditional probability,CP)严格的自然对数，可用SamEn(m,r,N)来表示。其中,N为长度，r为相似容限，维数为m及m+1。样本熵的具体算法如下：
将序列x(1),x(2),…,x(N)按顺序组成m维矢量，其中1≤i≤N-m+1，即
| ${{X}_{m\left( i \right)}}=\left[ x\left( i \right),x\left( i+1 \right),\ldots ,x\left( i+m-1 \right) \right]$ |
定义矢量Xm(i)和矢量Xm(j)之间的距离为d[Xm(i),Xm(j)]，即
| $d[{{X}_{m}}\left( i \right),{{X}_{m}}\left( j \right)]=\max x\left( i+k \right)-x\left( j+k \right)$ |
式中，1≤k≤m-1,1≤i,j≤N-m+1,i≠j。
给定相似容限r(r≥0)，对于每一个1≤i≤N-m，统计出d[Xm(i),Xm(j)]≤r的数目与矢量总数N-m-1的比值，记作
| ${{B}^{m}}_{i}\left( r \right)=\frac{1}{N-m-1}num\{d[{{X}_{m}}\left( i \right),{{X}_{m}}\left( j \right)]\le r\}$ |
式中，1≤j≤N-m,i≠j。求其对所有i的平均值
| ${{B}^{m}}\left( r \right)=\frac{1}{N-m}\sum\limits_{i-1}^{N-m}{{{B}^{m}}_{i}\left( r \right)}$ |
将维数加1，即对于m+1点矢量，得到Bm+1(r)。
理论上此序列的样本熵为
| $\text{SamEn}\left( m,r \right)=\underset{N\to \infty }{\mathop{\lim }}\,\{-\ln ({{B}^{m+1}}\left( r \right)/{{B}^{m}}\left( r \right))\}$ |
当N取有限值时，上式表示为
| $\text{SamEn}\left( m,r \right)=-\ln \{{{B}^{m+1}}\left( r \right)/{{B}^{m}}\left( r \right)\}$ |
SamEn的值随嵌入维数m和相似容限r的变化而变化，在特定的研究环境下，参数是一定的。一般情况下，m=1或2，r=0.1SD~0.25SD(SD是给定时间序列的标准差)。
2 基于极限学习机的癫痫脑电分类算法
2.1 极限学习机
假设有N个训练样本为(xi,yi)∈RN×Rm，其中，xi∈RN为输入，yi∈Rm为输出。具有M个隐藏层节点的标准单隐层前馈神经网络的数学模型为
| $\sum\limits_{j-1}^{M}{{{\beta }_{j}}g({{\omega }_{j}}{{x}_{i}}+{{b}_{j}})={{o}_{i}},i=1,2,\ldots ,N}$ |
其中，ωj为连接输入神经元与第j个隐藏层神经元的输入权值向量，βj为连接第j个隐藏层神经元与输出神经元的输出权值向量，oi为实际输出向量，bj为隐藏层神经元的偏置，g(·)是隐藏层神经元的激活函数。
如果该模型能够零误差地逼近训练样本的输出yi，即$\sum\limits_{i-1}^{M}{{{o}_{i}}-{{y}_{i}}=0}$，则存在βj，ωj，bj使下式成立：
| $\sum\limits_{j=1}^{M}{{{\beta }_{j}}g({{\omega }_{j}}{{x}_{i}}+{{b}_{j}})={{y}_{i}},i=1,2,\ldots ,N}$ |
将上式简化成紧凑的矩阵形式：
| $H\beta =Y$ |
其中H被称为神经网络隐藏层的输出矩阵。
Huang等[16]的研究表明，当神经元的激活函数任意可微时，单隐层前馈神经网络的训练误差可以接近无限小的正数ε，此时极限学习机的输入权值向量ωj和隐藏层偏置bj在训练过程中可保持不变，且可以随机分配。因此，训练过程等价为求线性系统的最小二乘解${\hat{\beta }}$，即
| $\left\| H\hat{\beta }-T \right\|=\underset{\beta }{\mathop{\min }}\,\left\| H\beta -T \right\|$ |
其解$\hat{\beta }\text{=}{{H}^{+}}Y$，H+是隐藏层输出矩阵H的Moore-Penrose广义逆。
2.2 癫痫脑电分类方法
本研究提出了一种基于小波多尺度分析和极限学习机的癫痫脑电分类算法。首先，利用小波多尺度分析对原始脑电信号进行5层多尺度分解，从中提取出4个不同频段的脑电信号。其次，利用两种非线性方法Hurst指数和样本熵分别对原始脑电信号和提取出的4个不同频段的脑电信号进行特征提取，得到了10维特征向量。最后，将非线性特征提取得到的10维特征向量作为极限学习机的输入，实现癫痫脑电分类的目的。详细流程图如图 2所示。
图2
癫痫脑电分类方法流程图
Figure2.
Flowchart of classification algorithm for epileptic EEG
3 实验结果与分析
本文实验数据来自德国波恩大学癫痫研究中心的脑电数据库。脑电信号采集系统采样频率为173.61 Hz，滤波带宽为0.53~40 Hz，脑电数据已通过人工去除伪迹。该数据集包含5个子集(记为Z、O、N、F和S)，其中，数据子集Z和O分别是健康人睁眼、闭眼时的正常脑电信号；数据子集N和F分别是癫痫患者发作间期病灶外和病灶区的脑电信号；数据子集S为癫痫患者发作期的脑电信号。数据集N、F和S均是由带状电极植入颅内，贴在海马结构表面记录得到的，其中癫痫发作期脑电信号(S)和癫痫发作间期的病灶区脑电信号(F)均是由带状电极贴在海马结构病灶区记录得到的，而癫痫发作间期的病灶外脑电信号(N)是带状电极贴在病灶区之外的另一个大脑半球海马结构处得到的。每个子集包含100段脑电信号，每段4 096点。将每段数据平均分为4小段，分别求每一小段的特征值，再对得到的4个特征值求均值，作为该段数据的特征值。本文中小波多尺度分析采用的是db3小波基函数。样本熵分析中，m=2，r=0.20SD。表 1为本文中涉及到的缩写及其含义。
表1
文中的缩写及含义说明
Table1.
Abbreviations and its meanings in this paper
利用Hurst和样本熵分别对癫痫发作期和发作间期的原始脑电信号和提取出的各频段脑电信号进行特征提取。为了方便比较癫痫发作期和发作间期脑电信号S两种特征值的差异，绘制了原始脑电信号发作期和发作间期特征值的分布序列和箱图，如图 3和图 4所示。由于篇幅所限，本文中仅展示了原始脑电信号发作期和发作间期特征值的分布序列和箱图,并未展示其他子频段信号特征值的对比情况。图 3是癫痫发作期和发作间期脑电信号特征值Hurst_S的分布序列和箱图。从图 3(a)可以看出，对大多数样本，发作间期的Hurst特征值明显大于发作期的Hurst特征值。从图 3(b)可以看出，发作间期的Hurst特征值的均值明显大于发作期的Hurst特征值均值。Hurst特征值对比结果表明，癫痫脑电发作期和发作间期均为持续稳定序列，发作间期的持续稳定性更强。图 4是癫痫发作期和发作间期脑电信号特征值SamEn_S的分布序列和箱图，对大多数样本，发作间期的样本熵特征值明显大于发作期的样本熵特征值，发作间期的样本熵特征值的均值也明显大于发作期的样本熵特征值的均值。样本熵特征值对比结果表明，发作间期的癫痫脑电比发作期的癫痫脑电更加不规则、更加复杂，产生新模式的概率也更高。
图3
癫痫发作期和发作间期脑电信号特征值Hurst_S的分布序列和箱图
(a)Hurst_S的分布序列；(b)Hurst _S的箱图
Figure3. Distribution series and box graphs of Hurst_S between epileptic ictal and interictal EEG(a) distribution series of Hurst_S； (b) box graphs of Hurst_S
图4
癫痫发作期和发作间期脑电信号特征值SamEn_S的分布序列和箱图
(a) SamEn_S的分布序列；(b)SamEn_S的箱图
Figure4. Distribution series and box graphs of SamEn_S between epileptic ictal and interictal EEG(a) distribution series of SamEn_S； (b) box graphs of SamEn_S
表 2是对癫痫脑电10个特征值的单一特征值分类能力的评价结果。从表中可以得到以下结论:① 任意一个特征值对癫痫脑电的分类准确率均大于50%；② Hurst_S获得了单一特征值的最高分类准确率89.5%；③ 利用原始脑电信号特征值得到的分类准确率大于从原始脑电信号提取出的各频段脑电信号特征值的分类准确率；④ 总体来看，脑电信号的Hurst特征值的分类性能优于样本熵特征值的分类性能。结果表明，本文提取的非线性特征能够较好地反映大脑在发作期和发作间期时的非线性动力学特性。
表2
癫痫脑电单一特征值分类能力评价
Table2.
Classification ability assessments of each single feature of the epileptic EEG
按照上述特征提取方法，可以得到100个癫痫患者发作期的样本，以及100个癫痫患者发作间期的样本，其中每个样本均包含一个10维特征向量。将10维特征向量作为极限学习机的输入，采用交叉留一法(假设有N个样本，将其中一个样本作为测试集，其它N－1个样本作为训练集，这个步骤一直持续到每个样本都被当作一次测试集为止，这样得到N个分类器和N个测试结果，用这N个测试结果的平均值来评价分类器的性能)，先训练极限学习机分类模型，然后测试极限学习机的分类性能。本文中极限学习机的隐藏层节点均为50个。表 3是基于不同核函数的癫痫脑电发作期和发作间期分类结果，结果表明基于sig和radbas核函数的极限学习机均能达到很好的分类结果，其中训练集上分类准确率均为100%，测试集上分类准确率为99.5%。
表3
基于不同核函数的癫痫脑电发作期和发作间期分类结果
Table3.
Classification results of epileptic ictal and interictal EEG based on different kernel functions
表 4给出了本文采用的方法与其他几种方法关于癫痫脑电发作期和发作间期的分类结果，表明本文采用的小波多尺度分析联合非线性特征提取和极限学习机得到的分类结果，优于其他非线性特征和ANN、SVM等的分类结果。
表4
本方法与其他方法关于癫痫脑电发作期和发作间期分类结果的比较
Table4.
Comparison of classification results of epileptic ictal and interictal EEG between using the method proposed in this paper and using other methods
4 结论
非线性动力学理论在脑电信号分析中有着很高的研究价值。Hurst和样本熵两种非线性方法均可使用较短的数据进行分析，有很好的抗噪能力，可用于生理信号的分析。本文利用Hurst和样本熵两种非线性方法对发作期和发作间期的脑电信号及由小波多尺度分解得到的各频段脑电信号进行了特征提取，结果表明所有的非线性特征值都能达到50%以上的分类准确率，说明本文使用的非线性特征提取方法能够反映发作期和发作间期的脑电信号的非线性动力学差异。将特征提取的10维特征向量作为输入，比较了基于不同核函数的极限学习机的分类性能，发现基于sig和radbas核函数的极限学习机均能达到很好的分类结果，其中训练集上分类准确率均为100%，测试集上分类准确率为99.5%。本文取得的分类结果优于其他研究中的结果，表明由小波多尺度分析联合非线性特征构建的极限学习机分类器能够很好地完成癫痫脑电分类任务，而且极限学习机分类器具有很好的泛化能力。
引言
癫痫是一种由多种病因引起的慢性脑功能障碍疾病，其典型特征是大脑局部神经元反复地、突然地过度放电，导致中枢神经系统功能失常[1]。全世界有1%~2%的人口患有癫痫[2]。癫痫患者临床表现为肌肉抽搐、意识丧失等，频繁的癫痫发作往往会给患者的身体带来极大的伤害，甚至危及生命，而且给患者家庭和社会带来沉重的负担。因此，加强癫痫疾病的前期诊断和后期治疗工作，具有极其重要的意义。
脑电在临床癫痫疾病的诊断、病灶定位、疗效评估中有着不可替代的作用。癫痫发作表现为高幅同步节律波，包括棘波、尖波、棘慢复合波和尖慢复合波等[3]。目前，癫痫疾病的诊断工作主要由经验丰富的医生肉眼观察长时间的脑电信号来完成，工作量巨大，而且受主观因素影响较大，不同医生给出的诊断结果有可能并不相同。因此，有必要开发一种癫痫脑电自动分类的方法。
癫痫脑电分类任务可以分为特征提取和分类两部分工作。非线性动力学相关方法已经在癫痫脑电信号的非线性特征提取工作中有了很好的应用。已有研究者将Lyapunov指数[4]、Hurst指数[5]、近似熵[6]、分形维数[7]、样本熵[8]、功率谱熵[9]、Lempel-Ziv复杂度[8]等非线性方法用于癫痫脑电信号的特征提取中。目前，癫痫脑电分类中使用较多的是K-最邻近分类器[10]、贝叶斯分类器[11]、决策树[12]、人工神经网络(artificial neural network,ANN)[11]、支持向量机(support vector machine,SVM)[13-15]等，并且已经取得了一定的成果。文献[14]比较了SVM、贝叶斯分类器、线性判别分析、K-最邻近分类器、ANN等几种分类方法在同一数据集上的分类效果，最终ANN取得了最高的97.7%的分类准确率。
Huang等[16]提出了一种新颖的单隐层前馈神经网络的机器学习算法，即极限学习机(extreme learning machine,ELM)，在学习机训练开始时，就随机给定它的隐藏层偏置和输入权值，且无需在训练过程中调整。与反向传播(back propagation,BP)神经网络相比，极限学习机既防止了学习过程中陷入局部最优解，又具有泛化能力强、训练速度快的优点，有较好的应用前景。
本文提出了一种基于小波多尺度分析和极限学习机的癫痫脑电分类算法。首先，对原始脑电信号进行小波多尺度分解，提取出不同频带的脑电信号；然后，利用Hurst指数和样本熵两种方法进行非线性特征提取；最后，应用极限学习机分类器对癫痫脑电进行分类，并在波恩大学的癫痫脑电数据库上测试分类器的性能。
1 基于小波多尺度分析的特征提取方法
1.1 小波多尺度分析
小波多尺度分析是一种将信号分解到多个尺度上的子空间的方法，分解后的信号在各个子空间同时具有时域和频域的分辨率，可以方便地分析信号的时频特性，在生物医学信号处理中有很好的应用。
对于信号S来说，就是将其分解为尺度空间的A1近似部分和细节空间的D1细节部分，然后对A1部分再次进行尺度和细节的分解，得到A2和D2两部分，再对A2进行如此分解。这里的近似信号对应于信号的低频部分，而细节信号对应于信号的高频部分。与傅里叶分析不同的是，高频部分是分层次的，是在不同分辨率下逐步产生的。图 1所示为信号S经过小波3层分解的示意图。
图1
小波3层分解示意图
Figure1.
Scheme of 3 layer wavelet decomposition
1.2 Hurst指数
Hurst指数是英国科学家H. E. Hurst在20世纪中叶提出的一种判断时间序列是否具有时间依赖性的参数，现在也用于研究混沌时间序列的预测问题[17]。本文采用R/S法，即重标极差方法(rescaled range analysis)。步骤如下：
对于给定的时间序列xi(i=1,2,…,N)，假定计算长度n，将时间序列划分为长度为n的M个子序列，第m个子序列的第t个数据记为xt,m(t=1,…,n,m=1,…,M) ，第m个子序列的重标极差比为
| ${{\left( R/S \right)}_{n,m}}=\frac{{{R}_{n,m}}}{{{S}_{n,m}}}$ |
式中，Rn,m为第m个子序列的极差，Sn,m为第m个子序列的标准差：
| ${{R}_{n,m}}=\underset{1\le k\le n}{\mathop{\max }}\,\sum\limits_{t=1}^{k}{k({{x}_{t,m}}-\overline{{{x}_{n,m}}})}-\underset{1\le k\le n}{\mathop{\min }}\,\sum\limits_{t=1}^{k}{k({{x}_{t,m}}-\overline{{{x}_{n,m}}})}$ |
| ${{S}_{n,m}}=\sqrt{\frac{1}{n}\sum\limits_{t=1}^{k}{{{({{x}_{t,m}}-\overline{{{x}_{n,m}}})}^{2}}}}$ |
式中，xn,m为第m个子序列。时间序列xi对应于计算长度n的重标极差为
| ${{\left( R/S \right)}_{n}}=\frac{1}{M}\sum\limits_{m=1}^{M}{{{\left( R/S \right)}_{n,m}}}$ |
改变计算长度n，重新计算重标极差比，可以得到一组计算长度与重标极差的序列。Hurst认为，重标极差序列(R/S)n与计算长度n之间存在指数关系如下：
| ${{\left( R/S \right)}_{n}}=C\times {{n}^{H}}$ |
式中，C为常数，H即为Hurst指数。通过取双对数：
| $\lg \left( R/S \right)=\lg \left( C \right)+H\lg \left( n \right)$ |
采用最小二乘法拟合公式(6)，可以估算出Hurst指数。
1.3 样本熵
Richman等[18]提出了一种评价给定时间序列复杂度的新方法，即样本熵(sample entropy,SamEn)，它是近似熵的一种改进方法。它是条件概率(conditional probability,CP)严格的自然对数，可用SamEn(m,r,N)来表示。其中,N为长度，r为相似容限，维数为m及m+1。样本熵的具体算法如下：
将序列x(1),x(2),…,x(N)按顺序组成m维矢量，其中1≤i≤N-m+1，即
| ${{X}_{m\left( i \right)}}=\left[ x\left( i \right),x\left( i+1 \right),\ldots ,x\left( i+m-1 \right) \right]$ |
定义矢量Xm(i)和矢量Xm(j)之间的距离为d[Xm(i),Xm(j)]，即
| $d[{{X}_{m}}\left( i \right),{{X}_{m}}\left( j \right)]=\max x\left( i+k \right)-x\left( j+k \right)$ |
式中，1≤k≤m-1,1≤i,j≤N-m+1,i≠j。
给定相似容限r(r≥0)，对于每一个1≤i≤N-m，统计出d[Xm(i),Xm(j)]≤r的数目与矢量总数N-m-1的比值，记作
| ${{B}^{m}}_{i}\left( r \right)=\frac{1}{N-m-1}num\{d[{{X}_{m}}\left( i \right),{{X}_{m}}\left( j \right)]\le r\}$ |
式中，1≤j≤N-m,i≠j。求其对所有i的平均值
| ${{B}^{m}}\left( r \right)=\frac{1}{N-m}\sum\limits_{i-1}^{N-m}{{{B}^{m}}_{i}\left( r \right)}$ |
将维数加1，即对于m+1点矢量，得到Bm+1(r)。
理论上此序列的样本熵为
| $\text{SamEn}\left( m,r \right)=\underset{N\to \infty }{\mathop{\lim }}\,\{-\ln ({{B}^{m+1}}\left( r \right)/{{B}^{m}}\left( r \right))\}$ |
当N取有限值时，上式表示为
| $\text{SamEn}\left( m,r \right)=-\ln \{{{B}^{m+1}}\left( r \right)/{{B}^{m}}\left( r \right)\}$ |
SamEn的值随嵌入维数m和相似容限r的变化而变化，在特定的研究环境下，参数是一定的。一般情况下，m=1或2，r=0.1SD~0.25SD(SD是给定时间序列的标准差)。
2 基于极限学习机的癫痫脑电分类算法
2.1 极限学习机
假设有N个训练样本为(xi,yi)∈RN×Rm，其中，xi∈RN为输入，yi∈Rm为输出。具有M个隐藏层节点的标准单隐层前馈神经网络的数学模型为
| $\sum\limits_{j-1}^{M}{{{\beta }_{j}}g({{\omega }_{j}}{{x}_{i}}+{{b}_{j}})={{o}_{i}},i=1,2,\ldots ,N}$ |
其中，ωj为连接输入神经元与第j个隐藏层神经元的输入权值向量，βj为连接第j个隐藏层神经元与输出神经元的输出权值向量，oi为实际输出向量，bj为隐藏层神经元的偏置，g(·)是隐藏层神经元的激活函数。
如果该模型能够零误差地逼近训练样本的输出yi，即$\sum\limits_{i-1}^{M}{{{o}_{i}}-{{y}_{i}}=0}$，则存在βj，ωj，bj使下式成立：
| $\sum\limits_{j=1}^{M}{{{\beta }_{j}}g({{\omega }_{j}}{{x}_{i}}+{{b}_{j}})={{y}_{i}},i=1,2,\ldots ,N}$ |
将上式简化成紧凑的矩阵形式：
| $H\beta =Y$ |
其中H被称为神经网络隐藏层的输出矩阵。
Huang等[16]的研究表明，当神经元的激活函数任意可微时，单隐层前馈神经网络的训练误差可以接近无限小的正数ε，此时极限学习机的输入权值向量ωj和隐藏层偏置bj在训练过程中可保持不变，且可以随机分配。因此，训练过程等价为求线性系统的最小二乘解${\hat{\beta }}$，即
| $\left\| H\hat{\beta }-T \right\|=\underset{\beta }{\mathop{\min }}\,\left\| H\beta -T \right\|$ |
其解$\hat{\beta }\text{=}{{H}^{+}}Y$，H+是隐藏层输出矩阵H的Moore-Penrose广义逆。
2.2 癫痫脑电分类方法
本研究提出了一种基于小波多尺度分析和极限学习机的癫痫脑电分类算法。首先，利用小波多尺度分析对原始脑电信号进行5层多尺度分解，从中提取出4个不同频段的脑电信号。其次，利用两种非线性方法Hurst指数和样本熵分别对原始脑电信号和提取出的4个不同频段的脑电信号进行特征提取，得到了10维特征向量。最后，将非线性特征提取得到的10维特征向量作为极限学习机的输入，实现癫痫脑电分类的目的。详细流程图如图 2所示。
图2
癫痫脑电分类方法流程图
Figure2.
Flowchart of classification algorithm for epileptic EEG
3 实验结果与分析
本文实验数据来自德国波恩大学癫痫研究中心的脑电数据库。脑电信号采集系统采样频率为173.61 Hz，滤波带宽为0.53~40 Hz，脑电数据已通过人工去除伪迹。该数据集包含5个子集(记为Z、O、N、F和S)，其中，数据子集Z和O分别是健康人睁眼、闭眼时的正常脑电信号；数据子集N和F分别是癫痫患者发作间期病灶外和病灶区的脑电信号；数据子集S为癫痫患者发作期的脑电信号。数据集N、F和S均是由带状电极植入颅内，贴在海马结构表面记录得到的，其中癫痫发作期脑电信号(S)和癫痫发作间期的病灶区脑电信号(F)均是由带状电极贴在海马结构病灶区记录得到的，而癫痫发作间期的病灶外脑电信号(N)是带状电极贴在病灶区之外的另一个大脑半球海马结构处得到的。每个子集包含100段脑电信号，每段4 096点。将每段数据平均分为4小段，分别求每一小段的特征值，再对得到的4个特征值求均值，作为该段数据的特征值。本文中小波多尺度分析采用的是db3小波基函数。样本熵分析中，m=2，r=0.20SD。表 1为本文中涉及到的缩写及其含义。
表1
文中的缩写及含义说明
Table1.
Abbreviations and its meanings in this paper
利用Hurst和样本熵分别对癫痫发作期和发作间期的原始脑电信号和提取出的各频段脑电信号进行特征提取。为了方便比较癫痫发作期和发作间期脑电信号S两种特征值的差异，绘制了原始脑电信号发作期和发作间期特征值的分布序列和箱图，如图 3和图 4所示。由于篇幅所限，本文中仅展示了原始脑电信号发作期和发作间期特征值的分布序列和箱图,并未展示其他子频段信号特征值的对比情况。图 3是癫痫发作期和发作间期脑电信号特征值Hurst_S的分布序列和箱图。从图 3(a)可以看出，对大多数样本，发作间期的Hurst特征值明显大于发作期的Hurst特征值。从图 3(b)可以看出，发作间期的Hurst特征值的均值明显大于发作期的Hurst特征值均值。Hurst特征值对比结果表明，癫痫脑电发作期和发作间期均为持续稳定序列，发作间期的持续稳定性更强。图 4是癫痫发作期和发作间期脑电信号特征值SamEn_S的分布序列和箱图，对大多数样本，发作间期的样本熵特征值明显大于发作期的样本熵特征值，发作间期的样本熵特征值的均值也明显大于发作期的样本熵特征值的均值。样本熵特征值对比结果表明，发作间期的癫痫脑电比发作期的癫痫脑电更加不规则、更加复杂，产生新模式的概率也更高。
图3
癫痫发作期和发作间期脑电信号特征值Hurst_S的分布序列和箱图
(a)Hurst_S的分布序列；(b)Hurst _S的箱图
Figure3. Distribution series and box graphs of Hurst_S between epileptic ictal and interictal EEG(a) distribution series of Hurst_S； (b) box graphs of Hurst_S
图4
癫痫发作期和发作间期脑电信号特征值SamEn_S的分布序列和箱图
(a) SamEn_S的分布序列；(b)SamEn_S的箱图
Figure4. Distribution series and box graphs of SamEn_S between epileptic ictal and interictal EEG(a) distribution series of SamEn_S； (b) box graphs of SamEn_S
表 2是对癫痫脑电10个特征值的单一特征值分类能力的评价结果。从表中可以得到以下结论:① 任意一个特征值对癫痫脑电的分类准确率均大于50%；② Hurst_S获得了单一特征值的最高分类准确率89.5%；③ 利用原始脑电信号特征值得到的分类准确率大于从原始脑电信号提取出的各频段脑电信号特征值的分类准确率；④ 总体来看，脑电信号的Hurst特征值的分类性能优于样本熵特征值的分类性能。结果表明，本文提取的非线性特征能够较好地反映大脑在发作期和发作间期时的非线性动力学特性。
表2
癫痫脑电单一特征值分类能力评价
Table2.
Classification ability assessments of each single feature of the epileptic EEG
按照上述特征提取方法，可以得到100个癫痫患者发作期的样本，以及100个癫痫患者发作间期的样本，其中每个样本均包含一个10维特征向量。将10维特征向量作为极限学习机的输入，采用交叉留一法(假设有N个样本，将其中一个样本作为测试集，其它N－1个样本作为训练集，这个步骤一直持续到每个样本都被当作一次测试集为止，这样得到N个分类器和N个测试结果，用这N个测试结果的平均值来评价分类器的性能)，先训练极限学习机分类模型，然后测试极限学习机的分类性能。本文中极限学习机的隐藏层节点均为50个。表 3是基于不同核函数的癫痫脑电发作期和发作间期分类结果，结果表明基于sig和radbas核函数的极限学习机均能达到很好的分类结果，其中训练集上分类准确率均为100%，测试集上分类准确率为99.5%。
表3
基于不同核函数的癫痫脑电发作期和发作间期分类结果
Table3.
Classification results of epileptic ictal and interictal EEG based on different kernel functions
表 4给出了本文采用的方法与其他几种方法关于癫痫脑电发作期和发作间期的分类结果，表明本文采用的小波多尺度分析联合非线性特征提取和极限学习机得到的分类结果，优于其他非线性特征和ANN、SVM等的分类结果。
表4
本方法与其他方法关于癫痫脑电发作期和发作间期分类结果的比较
Table4.
Comparison of classification results of epileptic ictal and interictal EEG between using the method proposed in this paper and using other methods
4 结论
非线性动力学理论在脑电信号分析中有着很高的研究价值。Hurst和样本熵两种非线性方法均可使用较短的数据进行分析，有很好的抗噪能力，可用于生理信号的分析。本文利用Hurst和样本熵两种非线性方法对发作期和发作间期的脑电信号及由小波多尺度分解得到的各频段脑电信号进行了特征提取，结果表明所有的非线性特征值都能达到50%以上的分类准确率，说明本文使用的非线性特征提取方法能够反映发作期和发作间期的脑电信号的非线性动力学差异。将特征提取的10维特征向量作为输入，比较了基于不同核函数的极限学习机的分类性能，发现基于sig和radbas核函数的极限学习机均能达到很好的分类结果，其中训练集上分类准确率均为100%，测试集上分类准确率为99.5%。本文取得的分类结果优于其他研究中的结果，表明由小波多尺度分析联合非线性特征构建的极限学习机分类器能够很好地完成癫痫脑电分类任务，而且极限学习机分类器具有很好的泛化能力。