脑电图（EEG）分析已被广泛应用于疾病的诊断，针对癫痫患者的脑电检测可及时对患者的发病情况作出判断，具有很强的实用价值，因此急需癫痫脑电自动检测、诊断分类技术。为实现患者正常期、癫痫发作间期和发作期各时段脑电的快速、高精度自动检测分类，本文提出一种基于样本熵（SampEn）与小波包能量特征提取结合纠错编码（ECOC）Real AdaBoost算法的脑电自动分类识别方法。将输入信号的样本熵值和4层小波包分解后的部分频段能量作为特征，并用纠错编码和Real AdaBoost算法相结合的方式对其进行分类。本文采用德国波恩大学癫痫数据库实验数据（含正常人清醒、睁眼与清醒、闭眼，癫痫患者间歇期致痫灶外与致痫灶内及癫痫发作期5组脑电信号）进行了方法有效性检验。研究结果表明，该方法有较强的脑电特征分类识别能力，尤其对癫痫间歇期脑电信号识别率提升显著，上述5组3个时期不同特征脑电信号的平均识别率可达96.78%，优于文献已报道的多种算法且有较好稳定性与运算速度及实时应用潜力，可在临床上对癫痫疾病的预报及检测起到良好的辅助决策作用。
引用本文： 张健钊, 姜威, 贲晛烨. 基于样本熵与小波包能量特征提取和Real AdaBoost算法的正常期、癫痫间歇与发作期的脑电自动检测. 生物医学工程学杂志, 2016, 33(6): 1031-1038. doi: 10.7507/1001-5515.20160166 复制
引言
癫痫(epilepsy)是由于大脑神经元异常放电所导致的短暂大脑功能障碍。癫痫患者在间歇期和发作期的脑电图(electroencephalogram,EEG)会有较为明显的异常表现，出现棘波、尖波及棘慢复合波等特征波形，因而脑电图作为一项重要手段已广泛应用于癫痫检测中。医学专家可以对脑电仪采集的脑电信号采取人工判读的方法，从中找出癫痫的特征波形，从而实现对癫痫的诊断。但这种方法耗时耗力，而且不同专家对同一段波形的判断也不尽一致，因此针对癫痫信号的自动检测和识别分类显得尤为必要。
近年来，众多学者对癫痫脑电的自动诊断与分类展开了多项研究。小波变换[1]、经验模态分解(empirical mode decomposition,EMD)[2]、熵[3-4]等常用信号分析方法均被应用于癫痫脑电的特征提取上；而目前流行的分类算法如人工神经网络(artificial neural network,ANN)[4]、支持向量机(support vector machine,SVM)[5]、Boosting[6]等被广泛用于脑电信号分类问题中，并达到了较高的分类精度。近年来广受关注的极限学习机(extreme learning machine,ELM)因其训练快速、分类性能出色等优点也受到了研究者们的青睐[7-8]。虽然这一领域的研究已取得了较为显著的成果，但是仍面临很多问题。在特征提取方面，仅局限于时、频域或非线性方法用于描述脑电信号会显得过于单一，使得分类器训练不够充分导致最终分类结果识别率不高，泛化能力不强，因此部分研究者尝试将多种特征应用于脑电信号的判别上[8-9]。在模式识别方面，基于支持向量机的分类方法虽然对于小样本问题有着较强的泛化能力，但对于规模较大的样本训练速度缓慢。而极限学习机虽然针对传统神经网络方法作了改进，但其初始输入参数随机产生，无法保证参数最优，使得最终分类效果受到一定影响。另外，目前的多数研究集中于脑电信号的二分类问题上，对信号的多时期分类的研究仍然较少。本文在此基础上，结合集成学习算法不依赖参数、运算速度快的优点提出一种对癫痫各时期脑电信号检测的方案。首先对截取的信号片段进行滤波处理去除相应高频成分，并将改进的快速样本熵(sample entropy,SampEn)计算方法引入其中，将样本熵值与部分节点的小波包能量组成特征向量送入基于纠错编码(error-correcting output codes,ECOC)的Real AdaBoost多类分类器中，并采用累计置信度作为判决依据，最终实现了对正常期、发作间期和发作期脑电信号的分类。本文所述方案的基本流程如图 1所示。
图1
脑电信号识别流程
Figure1.
Process of EEG signal recognition
1 实验数据
本文使用的实验数据来源于德国波恩大学癫痫研究中心数据库[10]，共分为5组数据集，编号分别为A和B(正常人清醒、睁眼时脑电数据和清醒、闭眼时脑电数据)、C和D(癫痫发作间期致痫灶外和致痫灶内脑电数据)及E(癫痫发作期脑电数据)。每组共有100段脑电信号，采样频率为173.6 Hz，每段共计4 097点数据，本文采用A~E组数据进行三个时期分类。将A、B组标记为正常时期信号数据，C、D组标记为发作间期信号数据，E组标记为发作期信号数据。由于原始数据段较长，实验前先将每段数据最后一个点舍弃，并等分为4段，每段1 024点数据，每个数据集得到400段脑电信号样本，共计2 000个样本。
2 特征提取方法
2.1 样本熵及其优化
样本熵是由Richman等[11]提出的度量时间序列信号自我相似程度及复杂度的方法，是对于近似熵(approximate entropy,ApEn)的改进。较近似熵而言，样本熵对数据长度及丢失数据不敏感，只需要较短的时间序列即可得到较为准确的估计。其基本思想是通过改变嵌入维度以计算序列产生新信息概率的大小。样本熵可以用参数m、 r、 N表示，m为嵌入维度，r为相似容限，N为时间序列长度。对于同一个序列，构造包含m个和m+1个矢量的新序列，分别计算矢量之间距离小于r的平均数量Bm(r)和Bm+1(r)，最终结果表示为$-\ln \left( \frac{{{B}^{m+1}}\left( r \right)}{{{B}^{m}}\left( r \right)} \right)$。通常情况下，m=1或2，r=0.1×std~0.25×std(std为序列标准差)。时间序列越复杂，自我相似程度越低，熵值越大。脑电信号作为一种有代表性的非平稳信号，其非线性特征的提取对于识别往往具有良好的效果，目前样本熵理论也被广泛应用于脑电信号分析中。实验表明，癫痫发作间期和发作期信号成分复杂度减小，样本熵值较正常时期明显降低，如图 2所示，为其中一组参数下(m=2,r=0.2×std,N=1 024)的两个时期样本熵分布规律比较。可见，发作期、发作间期与正常时期样本熵区分效果较好，适合作为分类特征。
图2
样本熵在两个时期的分布规律
Figure2.
Distribution of the sample entropy during the two periods
为满足实时计算的要求，洪波等[12]提出了一种快速近似熵算法，引入二值矩阵，通过矩阵元素间简单的逻辑运算简化了运算过程,本文在此基础上进一步优化了样本熵的运算过程。由于通常情况下取m=2，因此以下方法在m=2的情况下进行计算，具体步骤如下：
(1) 对给定的N点序列{u(i)}、阈值r及嵌入维度m，初始化(N-1)×(N-1)的二值矩阵M，定义第i行第j列的元素为aij，则：
| ${{a}_{ij}}=\left\{ \begin{array}{*{35}{l}} 1 & u\left( i \right)-u\left( j+1 \right)r且j\ge i \\ 0 & 其他 \\ \end{array} \right.$ |
其中0＜i,j＜N；
(2) 初始化N-m+1维的行向量V，定义其中每个元素为vk且vk =0(1≤k≤N-m+1)。对矩阵M的每一行i(1≤i≤N-m)，初始化si=0，则：
| ${{v}_{j+1}}\text{=}\left\{ \begin{array}{*{35}{l}} {{v}_{j+1}}+1 & {{a}_{ij}}\cap {{a}_{\left( i+1 \right)(j+1)}}=1 \\ {{v}_{j+1}} & {{a}_{ij}}\cap {{a}_{\left( i+1 \right)(j+1)}}=0 \\ \end{array} \right.$ |
| ${{s}_{i}}\text{=}\left\{ \begin{array}{*{35}{l}} {{s}_{i}}+1 & {{a}_{ij}}\cap {{a}_{\left( i+1 \right)(j+1)}}=1 \\ {{s}_{i}} & {{a}_{ij}}\cap {{a}_{\left( i+1 \right)(j+1)}}=0 \\ \end{array} \right.$ |
其中i≤j≤N-m。则两矢量间最大距离的平均值为：
| ${{B}^{m}}_{i}\left( r \right)=\frac{{{s}_{i}}+{{v}_{i}}}{N-m}$ |
(3)对所有的Bmi(r)求平均，结果记为Bm(r)，则：
| ${{B}^{m}}\left( r \right)=\frac{\sum\limits_{i=1}^{N-m}{{{B}^{m}}_{i}\left( r \right)+{{v}_{N-m+1}}/\left( N-m \right)}}{N-m+1}$ |
(4) 令m=m+1，重复执行第(2)(3)步，其中式(2)和式(3)修改为：
| ${{v}_{j+1}}\text{=}\left\{ \begin{array}{*{35}{l}} {{v}_{j+1}}+1 & {{a}_{ij}}\cap {{a}_{\left( i+1 \right)(j+1)}}\cap {{a}_{\left( i+2 \right)(j+2)}}=1 \\ {{v}_{j+1}} & ~{{a}_{ij}}\cap {{a}_{\left( i+1 \right)(j+1)}}\cap {{a}_{\left( i+2 \right)(j+2)}}=0 \\ \end{array} \right.$ |
| ${{s}_{i}}\text{=}\left\{ \begin{array}{*{35}{l}} {{s}_{i}}+1 & {{a}_{ij}}\cap {{a}_{\left( i+1 \right)(j+1)}}\cap {{a}_{\left( i+2 \right)(j+2)}}=1 \\ {{s}_{i}} & {{a}_{ij}}\cap {{a}_{\left( i+1 \right)(j+1)}}\cap {{a}_{\left( i+2 \right)(j+2)}}=0 \\ \end{array} \right.$ |
其余步骤不变，得到Bm+1(r)；
(5) 由于序列长度N为有限值，所以根据以下公式计算样本熵：
| $\text{SampEn}\left( m,r \right)=-\ln \frac{{{B}^{m+1}}\left( r \right)}{{{B}^{m}}\left( r \right)}$ |
由以上步骤可知，该方法省去了矩阵中重复部分的构造，并且通过引入向量V，避免了重复计算。原始算法和改进后快速算法的耗时如表 1所示。可见对于较长的序列，运算也可在1 s左右的时间完成，运算速度提高9倍以上。
表1
原方法与改进方法耗时比较
Table1.
Comparison of time-consuming between the original method and the improved method
2.2 小波包分解
小波分析作为一种时频分析的主要工具已被广泛应用于各个领域。相比于傅里叶变换和短时傅里叶变换，小波分析具有多分辨率分析的优点，可以在多个尺度上反映信号的局部细节。但传统的小波变换只对每一次分解后信号的低频成分做进一步分解而忽略了高频成分，因而不能充分地反映信号细节。小波包分解是一种能够继续将高频细节进一步分解的方法，因此对信号的描述更为精细。本文为提取脑电特征，采用Mallet快速算法对原始信号进行分解，并对相应节点的小波包系数进行重构得到最终系数。如图 3所示为三层小波包分解所得到的小波包分解树。脑电图出现异常时，出现大量特征波(棘波、尖波、慢波)，相应频带能量显著增强，两种状态下脑电节律有很大差异。因此，小波包频带能量可以有效表征正常与异常两种状态。频带j的小波包能量为：
| ${{E}_{j}}=\sum\limits_{i=1}^{N}{\|{{n}_{i}}{{\|}^{2}}}$ |
图3
三层小波包分解的小波包树
Figure3.
Wavelet packet tree of the three layer wavelet packet decomposition
其中N为相应频带系数个数，ni为小波包系数。
3 分类方法
AdaBoost算法是由Freund等[13]在20世纪90年代提出的集成学习算法。它的核心思想是对于同一训练集训练若干个弱分类器，并最终将弱分类器组合成强分类器。该算法通过每次迭代更新样本权重，用将错分的样本权重增大的方式让分类器更聚焦于难以被正确分类的样本上。但AdaBoost算法弱分类器的输出区间被限制在{+1,-1}之内，略显粗糙。Schapire等[14]提出的Real AdaBoost算法对此进行了改进，引入置信度的概念，弱分类器输出从原来的+1和-1拓展到了实数的范围，以反映某样本属于某个类别的可能性大小。二分类问题的Real AdaBoost算法步骤如下：
(1) 已知训练集集合S={(x1,y1),(x2,y2),…,(xn,yn)}，弱分类器集合为H。其中xi为特征向量，yi∈{+1,-1}。初始化样本权重D1(i)=1/n，i=1,2,3,…,n；
(2) 循环执行① ~④ 步，t为循环轮数，且t=1,2,3,…,T。
① 对样本划分为若干个互不相交的子空间，划分后区间个数为m。对划分后的每一区间计算+1和-1类的累计样本权重，记为W+1jt和WW-1jt，1≤j≤m；
② 计算归一化因子${{Z}_{t}}=\sum\limits_{j=1}^{m}{\sqrt{W_{+1}^{jt}W_{-1}^{jt}}}$，定义弱分类器函数ht(x)，对任意属于某区间的x有${{h}_{t}}\left( x \right)=\frac{1}{2}\text{ln}(W_{+1}^{jt}+\varepsilon )/(W_{-1}^{jt}+\varepsilon )$，其中1≤j≤m，ε为一很小的正数，选取${{h}_{t}}\left( x \right)=\arg \underset{h\in H}{\mathop{\min }}\,{{Z}_{t}}$；
③ 更新样本权重：${{D}_{t+1}}\left( i \right)=\frac{{{D}_{t}}\left( i \right)}{{{Z}_{t}}}{{e}^{-{{y}_{i}}{{h}_{t}}({{x}_{i}})}}$；
④ 当达到初始设定条件时停止迭代，执行第(3)步，否则继续执行循环体。
(3) 输出强分类器$H\left( x \right)=sign(\sum\limits_{t=1}^{T}{{{h}_{t}}\left( x \right)+b})$。b为平滑因子，默认情况下为0。
基于纠错编码的分类方法是一种能将多类问题转化成多个二类问题的方法[15]。其基本思想是对于一个K类问题，对每个类别均进行长度为L的二进制编码(log2K＜L≤2K-1-1)并生成一个K行L列的编码表，在码字的每一位上均进行二类分类并输出一个类别标签。经过L次分类之后生成一个码字并与表中的K个码字进行比较，最终的输出类别被判为与表中汉明距离(Hamming distance)最小的码字所对应的类别，此过程称为解码(decoding)。如表 2所示，是一个三类问题的编码方案，其中ci代表类别，fi为每个二分类问题的分类器函数。通常情况下，二分类器可选用决策树(decision tree,DT)、支持向量机等。但是当类别数目较少或码长较短时，最后的分类结果往往会出现多种可能性，也就是生成的编码表对该情况不具备纠错能力。如表 2所示的编码方案，当生成的码字为“100”时，可以对应c1或c3两个类别(汉明距离相等)，这将对最终的分类带来干扰，导致样本被错分的可能性增大。
表2
一个三类问题的编码表
Table2.
Code table of a three-class problem
由于Real AdaBoost算法的输出值衡量了样本属于某类的可能性大小,即置信度的累计量，所以可采用累计置信度作为解码规则。具体步骤为:
(1) 生成一个K行N列编码表，对每个码位均采用Real AdaBoost算法训练一个二类分类器，分类器个数为N，类别标签分别用+1、-1表示；
(2) 输入一个样本x并通过分类器输出累计置信度yi(x)，1≤i≤N；
(3) 给定一个码字{cj1,cj2,…,cjN} (1≤j≤K，cji∈{-1,+1})，对于所有的j，计算某个码字的累计置信度，即：
| ${{Y}_{j}}\left( x \right)=\sum\limits_{i=1}^{N}{{{c}_{ji}}{{y}_{i}}\left( x \right)}$ |
(4)取Yj(x)最大值所对应的类别为样本最终所属类别，即：
| $K\left( x \right)=\arg \underset{j}{\mathop{\max }}\,{{Y}_{j}}\left( x \right)$ |
利用置信度作为解码依据可改善传统基于汉明距离解码所存在的缺陷，对结果的描述更为精细。分类效果较好的分类器往往在最后的决策中占据主要地位，弥补了某些分类器在分类效果上的不足，因而具备一定的纠错能力，有效提高了预测精度。
4 数据处理及结果分析
首先对每段脑电信号进行低通滤波处理。根据采样定理，信号最大频率约为86.8 Hz。为减弱噪声等因素干扰，滤除高频成分，保留60 Hz以下的频段范围，如图 4所示，显示为处理后的各数据集波形。其中A和B分别为正常人清醒、睁眼时脑电数据和清醒、闭眼时脑电数据；C和D分别为癫痫发作间期致痫灶外和致痫灶内脑电数据；E为癫痫发作期脑电数据。
图4
处理后的各数据集脑电信号波形
Figure4.
EEG waveforms of each data set after being processed
样本熵在宏观上对非线性序列有较强的表征能力而无法对细节进行描述，而小波包分解则在细节上具有出色的描述能力，因此本文将两者结合作为特征向量，对输入脑电信号进行样本熵计算及4层小波包分解。文献[7]曾通过实验指出，样本熵参数选择m=2，r=0.2×std时分类效果较佳，因此本文采取上述参数进行计算。小波基函数选择Daubechies4小波。人类脑电活动范围在0.5~30 Hz之间，为消除脑电记录过程中出现的伪差干扰，同时需要尽可能包含脑电的各节律成分，选取(4,1)~(4,7)节点对应频带(层次结构如图 3所示)的小波包系数进行重构并计算每个节点的能量的十万分之一，与样本熵值共同组成8维特征向量送入分类器进行检测。如表 3所示，列出了各数据集随机抽取300个样本的特征平均值和小波包分解后各节点频率范围。其中E1~E7依次代表(4,1)~(4,7)节点对应频带能量的十万分之一。
表3
各数据集部分样本特征平均值和相应频带频率范围
Table3.
Mean value of each feature and frequency range in corresponding band of partial samples
该分类问题采用的编码表如表 4所示。为尽量保证训练集各类别数量的平衡，将正负样本按照适当比例进行随机抽取，并从每个数据集剩余样本中各随机抽取100个样本作为测试集，最终识别率取20次实验平均值。如表 5所示，显示了训练样本的具体分配情况。其中，Real AdaBoost算法采用单层决策树作为每次样本划分的弱分类器，并规定当训练误差收敛或达到最大迭代次数时即停止训练。
表4
该脑电分类问题所采用的编码表
Table4.
Code table used in this EEG classification problem
表5
各数据集训练样本个数分配情况
Table5.
Distribution of the training samples in each dataset
采用传统基于汉明距离的解码方法和本文方法所得到的各时期识别率和平均识别率的比较如图 5所示，结果用条形图和多次实验的结果标准差表示。可以看出，对于正常时期脑电信号，两者差距不大，但对于发作间期和发作期的分类，分类准确度分别提高2%和3%以上，且识别率的上下浮动较小，三个时期的综合识别率最高可达到97.50%。说明从整体分类效果而言，利用置信度作为解码规则可充分利用集成学习算法包含的信息，将输出值而不是简单的类别标签输出作为判别依据，有效提高了多分类问题中的识别率并具有更佳的稳定性。
图5
不同解码方法识别率比较
Figure5.
Comparison of the accuracy between different decoding methods
为了进一步显示采用多特征的优势，我们列出了采用不同特征提取方法所得出的三个时期的识别率及标准差，如图 6所示。从两者的比较中可以看出，两种特征提取方法均可以有效表征三个时期的脑电信号，但仅采用小波包分解后各频带能量作为特征的识别率明显低于小波包分解与样本熵相结合作为特征的平均识别率，三个时期的识别率分别为93.55%、89.60%、95.50%和97.73%、95.10%、97.50%。对三个时期的分类精度分别提升了4.2%、5.5%和2.0%，尤其对发作间期信号识别率有明显提升。这是因为样本熵对正常期及异常期(发作间期与发作期)的区分度较好，作为补充特征对最终效果影响较大。由此可见，信号的非线性特征和时频信息相结合的特征提取方案更为有效准确。
图6
不同特征提取方法识别率比较
Figure6.
Comparison of the accuracy between different feature extraction methods
其它同类型的研究在各数据集上的分类效果如表 6所示。可见在类似的问题上，所采用的特征提取方法和分类方法均对最终的分类效果有影响，而多特征结合则具有更大优势，结合时域和频域等信息共同表征脑电信号将是类似问题研究的重点。另外，对于较复杂的分类问题，分类算法的选取及优化往往对最终结果产生不可忽视的影响。综合来看，本文所采取的方法即使是在多种数据集的情况下也可达到可观的分类准确度，对癫痫各时期诊断有较强的参考价值。
表6
部分同类研究在各数据集上的分类识别率
Table6.
Recognition rate of partial same kind of research on datasets used above
5 结论
为解决癫痫多时期脑电分类问题，本文将Real AdaBoost算法与纠错编码相结合，并将样本熵和小波包能量作为特征，实现了对正常、发作间期和发作期脑电信号的有效识别。实验表明，采用多种特征可有效解决单一特征识别率不高的问题，且利用快速样本熵算法可满足实时性检测的要求。另外，用累计置信度代替汉明距离作为解码规则可充分利用每个二类分类器输出所包含的信息，并成功地把二类问题推广到多类问题上，因此本文提出的方法对癫痫脑电检测较为有效，在实际中有较强的应用价值，也为其它脑电分类问题，例如运动想象脑电分类、情绪分类等多分类问题提供了有效的解决方案。但是同样可以看到在多类问题中，发作间期脑电信号的识别率仍然略低于正常时期和发作期。在今后的工作中，我们将尝试其它的特征提取方法，提高发作间期脑电信号的识别率，以达到与其他两个时期识别率相当的水平。
引言
癫痫(epilepsy)是由于大脑神经元异常放电所导致的短暂大脑功能障碍。癫痫患者在间歇期和发作期的脑电图(electroencephalogram,EEG)会有较为明显的异常表现，出现棘波、尖波及棘慢复合波等特征波形，因而脑电图作为一项重要手段已广泛应用于癫痫检测中。医学专家可以对脑电仪采集的脑电信号采取人工判读的方法，从中找出癫痫的特征波形，从而实现对癫痫的诊断。但这种方法耗时耗力，而且不同专家对同一段波形的判断也不尽一致，因此针对癫痫信号的自动检测和识别分类显得尤为必要。
近年来，众多学者对癫痫脑电的自动诊断与分类展开了多项研究。小波变换[1]、经验模态分解(empirical mode decomposition,EMD)[2]、熵[3-4]等常用信号分析方法均被应用于癫痫脑电的特征提取上；而目前流行的分类算法如人工神经网络(artificial neural network,ANN)[4]、支持向量机(support vector machine,SVM)[5]、Boosting[6]等被广泛用于脑电信号分类问题中，并达到了较高的分类精度。近年来广受关注的极限学习机(extreme learning machine,ELM)因其训练快速、分类性能出色等优点也受到了研究者们的青睐[7-8]。虽然这一领域的研究已取得了较为显著的成果，但是仍面临很多问题。在特征提取方面，仅局限于时、频域或非线性方法用于描述脑电信号会显得过于单一，使得分类器训练不够充分导致最终分类结果识别率不高，泛化能力不强，因此部分研究者尝试将多种特征应用于脑电信号的判别上[8-9]。在模式识别方面，基于支持向量机的分类方法虽然对于小样本问题有着较强的泛化能力，但对于规模较大的样本训练速度缓慢。而极限学习机虽然针对传统神经网络方法作了改进，但其初始输入参数随机产生，无法保证参数最优，使得最终分类效果受到一定影响。另外，目前的多数研究集中于脑电信号的二分类问题上，对信号的多时期分类的研究仍然较少。本文在此基础上，结合集成学习算法不依赖参数、运算速度快的优点提出一种对癫痫各时期脑电信号检测的方案。首先对截取的信号片段进行滤波处理去除相应高频成分，并将改进的快速样本熵(sample entropy,SampEn)计算方法引入其中，将样本熵值与部分节点的小波包能量组成特征向量送入基于纠错编码(error-correcting output codes,ECOC)的Real AdaBoost多类分类器中，并采用累计置信度作为判决依据，最终实现了对正常期、发作间期和发作期脑电信号的分类。本文所述方案的基本流程如图 1所示。
图1
脑电信号识别流程
Figure1.
Process of EEG signal recognition
1 实验数据
本文使用的实验数据来源于德国波恩大学癫痫研究中心数据库[10]，共分为5组数据集，编号分别为A和B(正常人清醒、睁眼时脑电数据和清醒、闭眼时脑电数据)、C和D(癫痫发作间期致痫灶外和致痫灶内脑电数据)及E(癫痫发作期脑电数据)。每组共有100段脑电信号，采样频率为173.6 Hz，每段共计4 097点数据，本文采用A~E组数据进行三个时期分类。将A、B组标记为正常时期信号数据，C、D组标记为发作间期信号数据，E组标记为发作期信号数据。由于原始数据段较长，实验前先将每段数据最后一个点舍弃，并等分为4段，每段1 024点数据，每个数据集得到400段脑电信号样本，共计2 000个样本。
2 特征提取方法
2.1 样本熵及其优化
样本熵是由Richman等[11]提出的度量时间序列信号自我相似程度及复杂度的方法，是对于近似熵(approximate entropy,ApEn)的改进。较近似熵而言，样本熵对数据长度及丢失数据不敏感，只需要较短的时间序列即可得到较为准确的估计。其基本思想是通过改变嵌入维度以计算序列产生新信息概率的大小。样本熵可以用参数m、 r、 N表示，m为嵌入维度，r为相似容限，N为时间序列长度。对于同一个序列，构造包含m个和m+1个矢量的新序列，分别计算矢量之间距离小于r的平均数量Bm(r)和Bm+1(r)，最终结果表示为$-\ln \left( \frac{{{B}^{m+1}}\left( r \right)}{{{B}^{m}}\left( r \right)} \right)$。通常情况下，m=1或2，r=0.1×std~0.25×std(std为序列标准差)。时间序列越复杂，自我相似程度越低，熵值越大。脑电信号作为一种有代表性的非平稳信号，其非线性特征的提取对于识别往往具有良好的效果，目前样本熵理论也被广泛应用于脑电信号分析中。实验表明，癫痫发作间期和发作期信号成分复杂度减小，样本熵值较正常时期明显降低，如图 2所示，为其中一组参数下(m=2,r=0.2×std,N=1 024)的两个时期样本熵分布规律比较。可见，发作期、发作间期与正常时期样本熵区分效果较好，适合作为分类特征。
图2
样本熵在两个时期的分布规律
Figure2.
Distribution of the sample entropy during the two periods
为满足实时计算的要求，洪波等[12]提出了一种快速近似熵算法，引入二值矩阵，通过矩阵元素间简单的逻辑运算简化了运算过程,本文在此基础上进一步优化了样本熵的运算过程。由于通常情况下取m=2，因此以下方法在m=2的情况下进行计算，具体步骤如下：
(1) 对给定的N点序列{u(i)}、阈值r及嵌入维度m，初始化(N-1)×(N-1)的二值矩阵M，定义第i行第j列的元素为aij，则：
| ${{a}_{ij}}=\left\{ \begin{array}{*{35}{l}} 1 & u\left( i \right)-u\left( j+1 \right)r且j\ge i \\ 0 & 其他 \\ \end{array} \right.$ |
其中0＜i,j＜N；
(2) 初始化N-m+1维的行向量V，定义其中每个元素为vk且vk =0(1≤k≤N-m+1)。对矩阵M的每一行i(1≤i≤N-m)，初始化si=0，则：
| ${{v}_{j+1}}\text{=}\left\{ \begin{array}{*{35}{l}} {{v}_{j+1}}+1 & {{a}_{ij}}\cap {{a}_{\left( i+1 \right)(j+1)}}=1 \\ {{v}_{j+1}} & {{a}_{ij}}\cap {{a}_{\left( i+1 \right)(j+1)}}=0 \\ \end{array} \right.$ |
| ${{s}_{i}}\text{=}\left\{ \begin{array}{*{35}{l}} {{s}_{i}}+1 & {{a}_{ij}}\cap {{a}_{\left( i+1 \right)(j+1)}}=1 \\ {{s}_{i}} & {{a}_{ij}}\cap {{a}_{\left( i+1 \right)(j+1)}}=0 \\ \end{array} \right.$ |
其中i≤j≤N-m。则两矢量间最大距离的平均值为：
| ${{B}^{m}}_{i}\left( r \right)=\frac{{{s}_{i}}+{{v}_{i}}}{N-m}$ |
(3)对所有的Bmi(r)求平均，结果记为Bm(r)，则：
| ${{B}^{m}}\left( r \right)=\frac{\sum\limits_{i=1}^{N-m}{{{B}^{m}}_{i}\left( r \right)+{{v}_{N-m+1}}/\left( N-m \right)}}{N-m+1}$ |
(4) 令m=m+1，重复执行第(2)(3)步，其中式(2)和式(3)修改为：
| ${{v}_{j+1}}\text{=}\left\{ \begin{array}{*{35}{l}} {{v}_{j+1}}+1 & {{a}_{ij}}\cap {{a}_{\left( i+1 \right)(j+1)}}\cap {{a}_{\left( i+2 \right)(j+2)}}=1 \\ {{v}_{j+1}} & ~{{a}_{ij}}\cap {{a}_{\left( i+1 \right)(j+1)}}\cap {{a}_{\left( i+2 \right)(j+2)}}=0 \\ \end{array} \right.$ |
| ${{s}_{i}}\text{=}\left\{ \begin{array}{*{35}{l}} {{s}_{i}}+1 & {{a}_{ij}}\cap {{a}_{\left( i+1 \right)(j+1)}}\cap {{a}_{\left( i+2 \right)(j+2)}}=1 \\ {{s}_{i}} & {{a}_{ij}}\cap {{a}_{\left( i+1 \right)(j+1)}}\cap {{a}_{\left( i+2 \right)(j+2)}}=0 \\ \end{array} \right.$ |
其余步骤不变，得到Bm+1(r)；
(5) 由于序列长度N为有限值，所以根据以下公式计算样本熵：
| $\text{SampEn}\left( m,r \right)=-\ln \frac{{{B}^{m+1}}\left( r \right)}{{{B}^{m}}\left( r \right)}$ |
由以上步骤可知，该方法省去了矩阵中重复部分的构造，并且通过引入向量V，避免了重复计算。原始算法和改进后快速算法的耗时如表 1所示。可见对于较长的序列，运算也可在1 s左右的时间完成，运算速度提高9倍以上。
表1
原方法与改进方法耗时比较
Table1.
Comparison of time-consuming between the original method and the improved method
2.2 小波包分解
小波分析作为一种时频分析的主要工具已被广泛应用于各个领域。相比于傅里叶变换和短时傅里叶变换，小波分析具有多分辨率分析的优点，可以在多个尺度上反映信号的局部细节。但传统的小波变换只对每一次分解后信号的低频成分做进一步分解而忽略了高频成分，因而不能充分地反映信号细节。小波包分解是一种能够继续将高频细节进一步分解的方法，因此对信号的描述更为精细。本文为提取脑电特征，采用Mallet快速算法对原始信号进行分解，并对相应节点的小波包系数进行重构得到最终系数。如图 3所示为三层小波包分解所得到的小波包分解树。脑电图出现异常时，出现大量特征波(棘波、尖波、慢波)，相应频带能量显著增强，两种状态下脑电节律有很大差异。因此，小波包频带能量可以有效表征正常与异常两种状态。频带j的小波包能量为：
| ${{E}_{j}}=\sum\limits_{i=1}^{N}{\|{{n}_{i}}{{\|}^{2}}}$ |
图3
三层小波包分解的小波包树
Figure3.
Wavelet packet tree of the three layer wavelet packet decomposition
其中N为相应频带系数个数，ni为小波包系数。
3 分类方法
AdaBoost算法是由Freund等[13]在20世纪90年代提出的集成学习算法。它的核心思想是对于同一训练集训练若干个弱分类器，并最终将弱分类器组合成强分类器。该算法通过每次迭代更新样本权重，用将错分的样本权重增大的方式让分类器更聚焦于难以被正确分类的样本上。但AdaBoost算法弱分类器的输出区间被限制在{+1,-1}之内，略显粗糙。Schapire等[14]提出的Real AdaBoost算法对此进行了改进，引入置信度的概念，弱分类器输出从原来的+1和-1拓展到了实数的范围，以反映某样本属于某个类别的可能性大小。二分类问题的Real AdaBoost算法步骤如下：
(1) 已知训练集集合S={(x1,y1),(x2,y2),…,(xn,yn)}，弱分类器集合为H。其中xi为特征向量，yi∈{+1,-1}。初始化样本权重D1(i)=1/n，i=1,2,3,…,n；
(2) 循环执行① ~④ 步，t为循环轮数，且t=1,2,3,…,T。
① 对样本划分为若干个互不相交的子空间，划分后区间个数为m。对划分后的每一区间计算+1和-1类的累计样本权重，记为W+1jt和WW-1jt，1≤j≤m；
② 计算归一化因子${{Z}_{t}}=\sum\limits_{j=1}^{m}{\sqrt{W_{+1}^{jt}W_{-1}^{jt}}}$，定义弱分类器函数ht(x)，对任意属于某区间的x有${{h}_{t}}\left( x \right)=\frac{1}{2}\text{ln}(W_{+1}^{jt}+\varepsilon )/(W_{-1}^{jt}+\varepsilon )$，其中1≤j≤m，ε为一很小的正数，选取${{h}_{t}}\left( x \right)=\arg \underset{h\in H}{\mathop{\min }}\,{{Z}_{t}}$；
③ 更新样本权重：${{D}_{t+1}}\left( i \right)=\frac{{{D}_{t}}\left( i \right)}{{{Z}_{t}}}{{e}^{-{{y}_{i}}{{h}_{t}}({{x}_{i}})}}$；
④ 当达到初始设定条件时停止迭代，执行第(3)步，否则继续执行循环体。
(3) 输出强分类器$H\left( x \right)=sign(\sum\limits_{t=1}^{T}{{{h}_{t}}\left( x \right)+b})$。b为平滑因子，默认情况下为0。
基于纠错编码的分类方法是一种能将多类问题转化成多个二类问题的方法[15]。其基本思想是对于一个K类问题，对每个类别均进行长度为L的二进制编码(log2K＜L≤2K-1-1)并生成一个K行L列的编码表，在码字的每一位上均进行二类分类并输出一个类别标签。经过L次分类之后生成一个码字并与表中的K个码字进行比较，最终的输出类别被判为与表中汉明距离(Hamming distance)最小的码字所对应的类别，此过程称为解码(decoding)。如表 2所示，是一个三类问题的编码方案，其中ci代表类别，fi为每个二分类问题的分类器函数。通常情况下，二分类器可选用决策树(decision tree,DT)、支持向量机等。但是当类别数目较少或码长较短时，最后的分类结果往往会出现多种可能性，也就是生成的编码表对该情况不具备纠错能力。如表 2所示的编码方案，当生成的码字为“100”时，可以对应c1或c3两个类别(汉明距离相等)，这将对最终的分类带来干扰，导致样本被错分的可能性增大。
表2
一个三类问题的编码表
Table2.
Code table of a three-class problem
由于Real AdaBoost算法的输出值衡量了样本属于某类的可能性大小,即置信度的累计量，所以可采用累计置信度作为解码规则。具体步骤为:
(1) 生成一个K行N列编码表，对每个码位均采用Real AdaBoost算法训练一个二类分类器，分类器个数为N，类别标签分别用+1、-1表示；
(2) 输入一个样本x并通过分类器输出累计置信度yi(x)，1≤i≤N；
(3) 给定一个码字{cj1,cj2,…,cjN} (1≤j≤K，cji∈{-1,+1})，对于所有的j，计算某个码字的累计置信度，即：
| ${{Y}_{j}}\left( x \right)=\sum\limits_{i=1}^{N}{{{c}_{ji}}{{y}_{i}}\left( x \right)}$ |
(4)取Yj(x)最大值所对应的类别为样本最终所属类别，即：
| $K\left( x \right)=\arg \underset{j}{\mathop{\max }}\,{{Y}_{j}}\left( x \right)$ |
利用置信度作为解码依据可改善传统基于汉明距离解码所存在的缺陷，对结果的描述更为精细。分类效果较好的分类器往往在最后的决策中占据主要地位，弥补了某些分类器在分类效果上的不足，因而具备一定的纠错能力，有效提高了预测精度。
4 数据处理及结果分析
首先对每段脑电信号进行低通滤波处理。根据采样定理，信号最大频率约为86.8 Hz。为减弱噪声等因素干扰，滤除高频成分，保留60 Hz以下的频段范围，如图 4所示，显示为处理后的各数据集波形。其中A和B分别为正常人清醒、睁眼时脑电数据和清醒、闭眼时脑电数据；C和D分别为癫痫发作间期致痫灶外和致痫灶内脑电数据；E为癫痫发作期脑电数据。
图4
处理后的各数据集脑电信号波形
Figure4.
EEG waveforms of each data set after being processed
样本熵在宏观上对非线性序列有较强的表征能力而无法对细节进行描述，而小波包分解则在细节上具有出色的描述能力，因此本文将两者结合作为特征向量，对输入脑电信号进行样本熵计算及4层小波包分解。文献[7]曾通过实验指出，样本熵参数选择m=2，r=0.2×std时分类效果较佳，因此本文采取上述参数进行计算。小波基函数选择Daubechies4小波。人类脑电活动范围在0.5~30 Hz之间，为消除脑电记录过程中出现的伪差干扰，同时需要尽可能包含脑电的各节律成分，选取(4,1)~(4,7)节点对应频带(层次结构如图 3所示)的小波包系数进行重构并计算每个节点的能量的十万分之一，与样本熵值共同组成8维特征向量送入分类器进行检测。如表 3所示，列出了各数据集随机抽取300个样本的特征平均值和小波包分解后各节点频率范围。其中E1~E7依次代表(4,1)~(4,7)节点对应频带能量的十万分之一。
表3
各数据集部分样本特征平均值和相应频带频率范围
Table3.
Mean value of each feature and frequency range in corresponding band of partial samples
该分类问题采用的编码表如表 4所示。为尽量保证训练集各类别数量的平衡，将正负样本按照适当比例进行随机抽取，并从每个数据集剩余样本中各随机抽取100个样本作为测试集，最终识别率取20次实验平均值。如表 5所示，显示了训练样本的具体分配情况。其中，Real AdaBoost算法采用单层决策树作为每次样本划分的弱分类器，并规定当训练误差收敛或达到最大迭代次数时即停止训练。
表4
该脑电分类问题所采用的编码表
Table4.
Code table used in this EEG classification problem
表5
各数据集训练样本个数分配情况
Table5.
Distribution of the training samples in each dataset
采用传统基于汉明距离的解码方法和本文方法所得到的各时期识别率和平均识别率的比较如图 5所示，结果用条形图和多次实验的结果标准差表示。可以看出，对于正常时期脑电信号，两者差距不大，但对于发作间期和发作期的分类，分类准确度分别提高2%和3%以上，且识别率的上下浮动较小，三个时期的综合识别率最高可达到97.50%。说明从整体分类效果而言，利用置信度作为解码规则可充分利用集成学习算法包含的信息，将输出值而不是简单的类别标签输出作为判别依据，有效提高了多分类问题中的识别率并具有更佳的稳定性。
图5
不同解码方法识别率比较
Figure5.
Comparison of the accuracy between different decoding methods
为了进一步显示采用多特征的优势，我们列出了采用不同特征提取方法所得出的三个时期的识别率及标准差，如图 6所示。从两者的比较中可以看出，两种特征提取方法均可以有效表征三个时期的脑电信号，但仅采用小波包分解后各频带能量作为特征的识别率明显低于小波包分解与样本熵相结合作为特征的平均识别率，三个时期的识别率分别为93.55%、89.60%、95.50%和97.73%、95.10%、97.50%。对三个时期的分类精度分别提升了4.2%、5.5%和2.0%，尤其对发作间期信号识别率有明显提升。这是因为样本熵对正常期及异常期(发作间期与发作期)的区分度较好，作为补充特征对最终效果影响较大。由此可见，信号的非线性特征和时频信息相结合的特征提取方案更为有效准确。
图6
不同特征提取方法识别率比较
Figure6.
Comparison of the accuracy between different feature extraction methods
其它同类型的研究在各数据集上的分类效果如表 6所示。可见在类似的问题上，所采用的特征提取方法和分类方法均对最终的分类效果有影响，而多特征结合则具有更大优势，结合时域和频域等信息共同表征脑电信号将是类似问题研究的重点。另外，对于较复杂的分类问题，分类算法的选取及优化往往对最终结果产生不可忽视的影响。综合来看，本文所采取的方法即使是在多种数据集的情况下也可达到可观的分类准确度，对癫痫各时期诊断有较强的参考价值。
表6
部分同类研究在各数据集上的分类识别率
Table6.
Recognition rate of partial same kind of research on datasets used above
5 结论
为解决癫痫多时期脑电分类问题，本文将Real AdaBoost算法与纠错编码相结合，并将样本熵和小波包能量作为特征，实现了对正常、发作间期和发作期脑电信号的有效识别。实验表明，采用多种特征可有效解决单一特征识别率不高的问题，且利用快速样本熵算法可满足实时性检测的要求。另外，用累计置信度代替汉明距离作为解码规则可充分利用每个二类分类器输出所包含的信息，并成功地把二类问题推广到多类问题上，因此本文提出的方法对癫痫脑电检测较为有效，在实际中有较强的应用价值，也为其它脑电分类问题，例如运动想象脑电分类、情绪分类等多分类问题提供了有效的解决方案。但是同样可以看到在多类问题中，发作间期脑电信号的识别率仍然略低于正常时期和发作期。在今后的工作中，我们将尝试其它的特征提取方法，提高发作间期脑电信号的识别率，以达到与其他两个时期识别率相当的水平。