CN101488162B - 一种用于脑电信号自动评估的脑电信号特征提取方法 - Google Patents

Abstract

本发明属于脑电信号处理技术领域，具体是一种用于脑电信号自动评估的脑电信号特征提取方法，在所发明的脑电信号特征提取方法的基础上构造的脑电信号自动评估系统可用于癫痫的自动诊断、监护、疗效评估等场合。脑电信号特征提取方法由以下环节构成：对输入的时间序列进行相空间重构得到高维数据，对高维数据进行坐标变换以达到方向归一化，对高维数据进行平移以达到位置归一化，截取高维数据的一些局部流形，计算各局部流形对应的点集合的统计量作为脑电信号的特征。基于本发明方法的脑电信号评估系统可以在病人不发病时检测脑电异常。

Description

一种用于脑电信号自动评估的脑电信号特征提取方法

技术领域

本发明属于脑电信号处理技术领域，具体涉及一种脑电信号特征提取方法。

背景技术

癫痫是一种常见病，患者发病时常失去知觉，如果患者在驾驶机动车时发病会危急自己和他人生命安全，因此癫痫的诊断非常重要。癫痫有可能是原发性的，也可能是因为脑部其它疾病和损伤引起的并发症，对于后一种情况一般在脑部都存在病灶，而前一种情况脑部没有明显病灶。

现有的用于癫痫诊断的医疗器械的工作原理是：首先由电极帽连接到放大器再连接到计算机组成的硬件系统将脑电信号采集进计算机，然后由经过训练的医生用肉眼观察计算机中记录的脑电图，根据脑电图的波形判断是否有癫痫特征出现。根据国际标准，门诊病人的脑电图需要记录20分钟，医生根据20分钟脑电图的总体涨落趋势判断是否有异常波形出现，一般一名医生一天只能接诊数量非常有限的病人。如果患者发病不是非常频繁，则医生很难在20分钟的脑电图记录中观察到异常，此时患者需要住院观察并接受24小时不间断的脑电图跟踪记录，在这种情况下，医生可能要花数个小时人工分析脑电图，效率非常低。另一个难题是——要训练一个能够正确判断脑电图中有无癫痫特征的医生需要花很多年，目前只有大城市中极少数医生能够较为准确地识别脑电图中的癫痫特征，而中小城市和农村几乎没有在脑电图分析方面训练有素的医生，原因是：脑电信号类似随机信号，从时域波形上很难观察到有规律的特征。

基于以上现状，用机器实现脑电信号模式的自动分析和分类对于癫痫诊断是一个很好的解决方案，带来的好处有：(1)采用机器自动诊断解决了培养一个能够分析脑电图的医生需要花费很多年时间的问题，可以方便地在广大地区部署机器诊断系统。(2)机器诊断比人工诊断效率要提高很多，由于脑电信号的随机性，医生需要记录很长一段信号才能观察出信号中是否存在癫痫特征，这是因为肉眼只能观察到波形的涨落；而机器可以通过计算得到有规律的统计量，对于一段很短的信号也可以判断出其中是否具有癫痫特征。(3)脑电信号不适合人工分析，因为脑电信号随机性强，从时域波形中很难观察到有规律的特征；而机器自动分析可以是不基于波形的，机器通过计算后可以从随机信号中计算出有规律的统计量，因此更适合分析脑电这样具有很强随机性的信号。

目前国内外关于脑电信号自动分析方法和装置的研究很多。早期的研究思路是用分类器(如神经网络)直接对脑电信号进行分类，这方面的工作有陈俊强和刘书朋提出的中国发明专利99124032.4“全自动定量检测脑电图中癫痫样放电的装置”，这种方法的缺点是：由于脑电信号的时域波形几乎是无规律的，不经过特征提取就直接对脑电信号进行分类在实际中很难达到满意的效果。Ina Pichlmayr和Olaf Eckert在[“Method and apparatus for the evaluation ofEEG data”，US patent 5846208]中提出用子带功率谱作为脑电信号的特征、用神经网络作为分类器对脑电信号进行分类。近期的研究大都聚焦在寻找能够区分不同类别脑电信号的特征量，为此，来自物理学、信息科学等领域的工作者对脑电信号特征进行了大量研究，目前达成共识的是——脑电信号的非线性混沌特征对于定量刻画脑电信号和诊断癫痫是非常重要的参量。在物理学领域的杂志上，这方面的工作有：由R.G.Andrzejak等在[PHYSICAL REVIEWE 64(6)：Art.No.061907 Part 1，DEC 2001]提出的“Indications of nonlinear deterministic andfinite-dimensional structures in time series of brain electrical activity：Dependence on recordingregion and brain state”；由J.L.P.Velazquez等在[PHYSICAD-NONLINEAR PHENOMENA，186(3-4)：205-220 DEC 15 2003]中提出的“Dynamical regimes underlying epileptiform events：role ofinstabilities and bifurcations in brain activity”；由T.Gautama等在[PHYSICAL REVIEW E，67(4)：Art.No.046204 Part 2 APR 2003]中提出的“Indications of nonlinear structures in brain electricalactivity”。信息科学领域的研究工作不但包括脑电信号的特征提取、还包括脑电信号的分类，基本的处理方法是：首先利用各种特征提取方法对采集到的脑电信号进行参量计算以得到有规律的模式信息(比如前面提到的一些脑电信号非线性混沌特征)，再将提取出的脑电信号特征参量输入分类器进行分类，分类器的输出即为脑电信号的类别(如：将输入的脑电信号分为正常、癫痫不活跃期、癫痫活跃期中的某一类)，这样可以完成癫痫的自动诊断，这方面的工作有：由K.Lehnertz等在[IEEE ENGINEERING IN MEDICINE AND BIOLOGY MAGAZINE，22(1)：57-63 JAN-FEB 2003]中提出的“Seizure prediction by nonlinear EEG analysis”；由N.F.Guler等在[EXPERT SYSTEMS WITH APPLICATIONS，29(3)：506-514 OCT 2005]中提出的“Recurrent neural networks employing Lyapunov exponents for EEG signals classification”；由I.Guler和E.D.Ubeyli在[IEEE TRANSACTIONS ON INFORMATION TECHNOLOGY INBIOMEDICINE，11(2)：117-126 MAR 2007]中提出的“Multiclass support vector machines forEEG-signals classification”。

发明内容

本发明的目的在于提供一种用于脑电信号自动评估的脑电信号特征提取方法，在本发明提出的脑电信号特征提取方法的基础上可以构造一个用于脑部疾病诊断、监护、疗效评估的脑电信号自动评估系统。

以下先对实现本发明的技术方案所涉及、使用的一些模式识别领域的技术名词、术语作如下定义和解释。

脑电信号的类别：有正常、病症不活跃期、病症活跃期等多个类别。

特征提取：目的是从看似没有规律的脑电波形中抽取出有规律的信息，实际上是通过一种计算方法对脑电信号进行处理以得到一些参量，这些参量在模式识别领域称为特征，经过特征提取得到的这些参量应该具有比较明显的规律性，即：同类别信号提取出的特征参量应该具有相似性，不同类别信号的特征参量具有差异性。这里所称的“特征”是模式识别领域的一个术语，英文名称为“Feature”，它与矩阵论中所指的“特征值”不同，矩阵论中的术语“特征值”的英文名称为“Eigenvalue”。

特征向量：从一个信号中提取出的所有特征参量堆积起来构成一个向量，称为特征向量。这里所称的“特征向量”是模式识别领域的术语，英文称为“Feature Vector”，它不同于矩阵理论中所指的“特征向量”，矩阵论中所指的“特征向量”的英文名称是“Eigenvector”。

分类器：可以看作一种映射，当输入一个特征向量时，分类器输出一个类别号；由于一个特征向量对应一个信号(一个特征向量是从一个信号中提取出来的)，所以分类器输出的类别号也就是输入的特征向量对应的信号的类别号，由此可以实现脑电信号的自动分类，脑电信号的自动分类也称为脑电信号的自动识别或自动诊断。

分类器的训练：分类器一般都有很多参数(如支持向量机)，只有通过训练算法才能把分类器的参数调整为合适的值，当分类器的参数调整到合适的值后，分类器才能较好地工作(以分类正确率衡量)。

模式识别系统：至少包括特征提取和分类器这两个环节，特征提取可以从信号中计算出一些参量并形成特征向量，分类器可以把特征向量映射为某个类别号。

脑电信号分类/脑电信号识别：为脑电信号分类(也称识别)设计的模式识别系统，在此基础上设计的医疗器械可以用于疾病的自动诊断。

主元分析：目的是求流形或点集合的主方向，即所有点投影后坐标值方差最大的方向。

矩：一种统计量，模式识别中的不变矩(英文为Moment Invariants)可以刻画不同几何形状或流形之间的区别。

本发明是一种用于脑电信号自动评估的脑电信号特征提取方法，这里所述的脑电信号特征提取方法至少包括以下计算步骤：

(a)对输入的脑电时间序列S＝[S₁，S₂，...，S_N]进行状态空间重构(也称相空间重构)：先选取参数J和M，一般取5≤M≤30，取J为脑电时间序列自相关函数第一次下降到最大值的1/e时对应的时刻，再任意选取L个满足1≤t(1)＜t(2)＜...＜t(L)≤N-(M-1)J的正整数{t(i)|i＝1，2，...，L}，令X_i＝[S_t(i)，S_t(i)+J，...，S_t(i)+(M-1)J]^T，则称[X₁，X₂，...，X_L]或{X₁，X₂，...，X_L}为时间序列S的一个状态空间重构，{X₁，X₂，...，X_L}可以看作M维空间的一个点集合或流形，[X₁，X₂，...，X_L]可以看作M维空间的一个轨迹；以上状态空间重构也可以以行向量的形式实现，即令X_i＝[S_t(i)，S_t(i)+J，...，S_t(i)+(M-1)J]，i＝1，2，...，L；

(b)对重构的流形进行方向归一化：根据流形在各个方向上分布的特点定出基准方向，参照基准方向对流形进行坐标变换；

(c)对重构的流形进行位置归一化：根据流形的形状和结构特点确定基准点，参照基准点对流形进行平移操作；

(d)截取重构的流形的一些局部区域；

(e)对截取的局部流形或轨迹进行统计量的计算，把计算得到的统计量作为输入的脑电信号的特征。

上面所述的特征提取方法中的计算步骤(b)、(c)、(d)、(e)各有若干种实现方案。上面所述的特征提取方法中的计算步骤(b)的两种实现方案在下面分别用(b1)和(b2)表示，计算步骤(c)的两种实现方案在下面分别用(c1)和(c2)表示，计算步骤(d)的两种实现方案在下面分别用(d1)和(d2)表示，计算步骤(e)的三种实现方案在下面分别用(e1)、(e2)、(e3)表示，具体计算方法如下：

(b1)采用主元分析法定出基准方向，令X＝[X₁，X₂，...，X_L]，X是一个M×L维的矩阵，计算XX^T的M个特征值，并按从大到小的顺序排列得到λ₁≥λ₂≥...≥λ_M，这M个矩阵特征值对应的矩阵XX^T的特征向量{U₁，U₂，...，U_M}作为主轴，令U＝[U₁，U₂，...，U_M]，利用Y＝U^TX对X进行坐标变换，得到新的坐标Y，Y对应经过方向归一化的流形。

(b2)首先计算{X₁，X₂，...，X_L}的算术平均中心 $\stackrel{\‾}{X}=\frac{1}{L}\underset{i=1}{\overset{L}{\mathrm{\Σ}}}{X}_i,$ X_i与

两点可以确定一条直线L_i，设X_j到直线L_i距离为d_ij，求 $\{D_i=\frac{1}{L-1}\underset{j=1}{\overset{L}{\mathrm{\Σ}}}{d}_{\mathrm{ij}}|i=\mathrm{1,2},...,L\},$ D_i是所有点到直线L_i的平均距离，计算 $k=\arg \underset{i}{\min }\left\{D_i\right\},$ D_k是{D_i|i＝1，2，...，L}中的最小值，选 $U_0=\frac{X_k-\stackrel{\‾}{X}}{\left|\right|X_k-\stackrel{\‾}{X}\left|\right|}$ 为第一主轴，这里

表示向量

的模，构造一个经过中心且与第一主轴垂直U₀的超平面，计算{X₁，X₂，...，X_L}在此超平面的投影值{P₁，P₂，...，P_L}，这里{P_i|i＝1，2，...，L}是M-1维向量，对P＝[P₁，P₂，...，P_L]进行主元分析，找到矩阵PP^T的M-1个从大到小排序的特征值λ₁≥λ₂≥...≥λ_M-1和其对应的特征向量U₁，U₂，...，U_M-1，{U₀，U₁，...，U_M-1}构成一个新的坐标系，用U＝[U₀，U₁，...，U_M-1]对X＝[X₁，X₂，...，X_L]进行坐标变换，得到Y＝U^TX，Y对应经过方向归一化的流形。

(c1)用矩的方法计算出流形的中心作为平移操作的基准点，设Y_i＝[Y_i1，Y_i2，...，Y_iM]是M维空间的一个点，计算{Y_i|i＝1，2，...，L}的阶数为[p₁，p₂，...，p_M]的矩 $m_{p_1{p}_2...p_M}=\underset{j=1}{\overset{L}{\mathrm{\Σ}}}\underset{i=1}{\overset{M}{\mathrm{\Π}}}{Y}_{\mathrm{ji}}^{p_i},$ 当{p₁，p₂，...，p_M}中p_k＝1且其它元素为0时记 ${\stackrel{\‾}{m}}_k=m_{p_1{p}_2...p_M},$ 当{p₁，p₂，...，p_M}中所有元素为0时记 ${\stackrel{\‾}{m}}_0=m_{p_1{p}_2...p_M},$ 所求的中心为 $\stackrel{\‾}{Y}={[\frac{{\stackrel{\‾}{m}}_1}{{\stackrel{\‾}{m}}_0},\frac{{\stackrel{\‾}{m}}_2}{{\stackrel{\‾}{m}}_0},\·\·\·,\frac{{\stackrel{\‾}{m}}_M}{{\stackrel{\‾}{m}}_0}]}^T,$ 将所有点的坐标减去

就完成了流形的平移操作，即令 $Y_j\←Y_j-\stackrel{\‾}{Y},$ 这里j＝1，2，...，L。

(c2)计算各坐标轴上各点坐标的算术平均值作为对流形进行平移操作的基准点，基准点的计算方法为 $\stackrel{\‾}{Y}=\frac{1}{L}\underset{i=1}{\overset{L}{\mathrm{\Σ}}}{Y}_i,$ 这里{Y_i|i＝1，2，...，L}是整个流形包含的点，将所有点的坐标减去

就完成了流形的平移操作，即令 $Y_j\←Y_j-\stackrel{\‾}{Y},$ 这里j＝1，2，...，L。

(d1)构造一些与主轴U_i垂直的超平面，称为Poincare截面，这里i∈[1，2，...，M]，Poincare截面与流形构成的轨迹相交，每个Poincare截面与轨迹的交点的计算方法如下：假设一个Poincare截面与主轴U_i垂直相交且交点在U_i轴的坐标值为Y₀，如果轨迹[Y₁，Y₂，...，Y_L]的第j个点Y_j＝[Y_j1，Y_j2，...，Y_jM]的第i个坐标值满足Y_ji＝Y₀，则Y_j为交点，此外，如果轨迹的第j个点和第j+1个点的第i个坐标值满足(Y_ji-Y₀)(Y_(j+1)i-Y₀)＜0，则可以利用线性插值法求得一个交点 $Y_j^{\′}=Y_j+\frac{Y_{\mathrm{ji}}-Y_0}{Y_{\mathrm{ji}}-Y_{(j+1)i}}(Y_{j+1}-Y_j),$ 如此可以求出轨迹穿越每个Poincare截面的交点，轨迹穿越每个Poincare截面的交点构成一个点集合，这个点集合即为所求的局部流形。

(d2)构造一对与主轴U_i垂直的超平面，设这两个Poincare截面与U_i轴的交点在U_i轴上的坐标值分别为Y₀₁和Y₀₂，求落入这两个Poincare截面之间的轨迹上的点，具体计算方法是：如果轨迹[Y₁，Y₂，...，Y_L]的第j个点Y_j＝[Y_j1，Y_j2，...，Y_jM]的第i个坐标值满足min{Y₀₁，Y₀₂}≤Y_ji≤max{Y₀₁，Y₀₂}，则判定Y_j为落入两个Poincare截面之间的点，构造若干对垂直于各个主轴的Poincare截面，再按照上述计算方法求出落入每对Poincare截面之间的点，落入每对Poincare截面之间的点构成一个点集合，这个点集合即为所求的局部流形。

(e1)统计每个局部流形对应的点集合中包含的点的个数，把所述的各个点集合包含的点的个数作为脑电信号的特征。

(e2)分别对每个局部流形对应的点集合求矩，把所述的对各个点集合求出的矩作为脑电信号的特征。

(e3)分别对每个局部流形对应的点集合进行主元分析，把所述的各个点集合经过主元分析得到的矩阵特征值作为脑电信号的特征。

本发明的优点是：

(1)基于人工分析的脑电信号评估的一个缺点是，在病人不发病时无法判断脑电信号有无异常，因为图形的分辨率和人类视觉的辨识能力都是有限的，而基于本发明方法的脑电信号自动评估方法可以在病人不发病时检测出脑电信号的异常。(2)电极与脑部接触的位置、受测者处于睁眼还是闭眼状态等因素都会使采集到的脑电信号具有较大的变化，基于本发明方法的脑电自动评估方法对于电极的位置不敏感，当脑电信号取自与病灶位置相反的位置上的电极时，采用本发明提出的脑电信号自动评估方法依然能够对脑电信号正确分类，另外，无论受测者睁眼还是闭眼，利用本发明提出的脑电信号自动评估方法都能够对脑电信号正确分类，我们采用了三类实际数据进行测试，第一类数据包括两个集合：集合A是正常人睁眼时采集的脑电信号、集合B是正常人闭眼时采集的脑电信号，第二类数据也包括两个集合：集合D是癫痫患者不发病时在病灶位置采集的脑电信号、集合C是癫痫患者不发病时在病灶相反位置采集的脑电信号，第三类即集合E是病人发病时采集的脑电信号，按照本发明提出的方法对以上三类信号进行分类，正确识别率分别为99.9％、99.5％、98.4％，三类平均正确识别率为99.44％。(3)基于本发明方法的脑电信号自动评估方法只需要较短时间的信号就可以分辨出信号的类别，上述所有测试用到的信号的长度均为23.6秒，这表明本发明提出的脑电自动评估方法只需要很短时间的信号记录就可以工作，不需要象人工分析那样长时间记录脑电信号以观察信号的涨落趋势。

附图说明

图1：脑电信号自动评估系统的组成框图

图2：脑电信号特征提取算法的流程图

图中标号：1为电极帽，2为放大器，3为计算机系统，4为脑电信号采集模块，5为特征提取模块，6为分类器，7为显示器，8为USB接口

具体实施方式

本发明是一种用于脑电信号自动评估的脑电信号特征提取方法，在本发明提出的脑电信号特征提取方法的基础上可以构造一个脑电信号自动评估系统，本发明提出的脑电信号特征提取方法是整个脑电信号自动评估系统的一个组成部分，整个脑电信号自动评估系统的组成见附图1，由以下几部分构成：

(1)脑电信号采集：功能是将脑电信号采集进计算机，先由电极帽连接到放大器、再将放大器通过USB接口连接到计算机系统(这里所指的计算机系统除了包括主机外，还包括键盘、鼠标、显示器等输入和输出设备)，调用驱动程序就可以将脑电信号采集进计算机并进行存储和处理，电极帽的电极个数和位置都遵循国际标准。电极帽、放大器、计算机系统在市场上都可以购得，例如：可以购置德国Brain Products公司生产的电极帽、北京中科新拓仪器有限公司生产的脑电放大器、美国惠普公司生产的笔记本电脑。

(2)特征提取：功能是通过一系列的计算步骤从前一步采集到的看似杂乱无章的脑电信号中获取有助于脑电信号分类的有规律的信息，一般由软件模块实现，但也可以用硬件实现，本发明提出的技术方案的计算流程图见附图2，有多种具体的实现方法，后面将列举2个实施例。

(3)分类器：功能是根据前一步计算获得的脑电信号的特征自动判断脑电信号的类别，一般由软件模块实现，但也可以由硬件实现，分类器可以看作一个映射，将前一步计算得到的脑电信号的特征输入分类器，分类器会自动输出一个类别号，将分类器输出的类别号(例如：正常、疾病活跃期、疾病不活跃期)显示在计算机屏幕上就完成了整个脑电信号自动评估过程；分类器输出的关于脑电信号的类别信息可以用于医疗过程中的诊断、监护、疗效评估等；常用的分类器有k近邻分类器、贝叶斯分类器、神经网络、支持向量机等，这里使用支持向量机作为分类器，因为支持向量机有很多开放源代码的软件实现，如LIBSVM(见http://www，csie.ntu.edu.tw/～cjlin/libsvm)；支持向量机的参数需要经过训练调整到较佳值后，在分类时才能较好地工作，LIBSVM软件包里提供有训练工具，支持向量机的参数训练方法可参考Nello Cristianini和John Shawa-Taylor合著的《An introduction to support vector machinesand other kernel-based learning methods》，本书2000年由Cambridge University Press出版。

脑电信号自动评估实施例1：

步骤1：这一步完成脑电信号的采集，将电极帽连接到放大器、将放大器通过USB接口连接到计算机系统、调用驱动程序将脑电信号采集进计算机；

步骤2：这一步完成输入时间序列S＝[S₁，S₂，...，S_N]的状态空间重构，首先选定两个参数J和M，取J为[S₁，S₂，...，S_N]的自相关函数第一次下降到最大值的1/e时对应的时刻，令M＝15，令{t(i)＝i|i＝1，2，...，N-(M-1)J}，令

$X_i=\left[\begin{array}{c}{S}_{t\left(i\right)}\\ S_{t\left(i\right)+J}\\ \·\\ \·\\ \·\\ S_{t\left(i\right)+(M-1)J}\end{array}\right]=\left[\begin{array}{c}{S}_i\\ S_{i+J}\\ \·\\ \·\\ \·\\ S_{i+(M-1)J}\end{array}\right]$

则得到一个矩阵

$X=[X_1,X_2,\·\·\·X_{N-(M-1)J}]=\left[\begin{array}{cccc}{S}_1& S_2& \·\·\·& S_{N-(M-1)J}\\ S_{1+J}& S_{2+J}& \·\·\·& S_{N-(M-2)J}\\ \·& \·& \·& \·\\ \·& \·& \·& \·\\ \·& \·& \·& \·\\ S_{1+(M-1)J}& S_{2+(M-1)J}& \·\·\·& S_N\end{array}\right]$

经过状态空间重构，从原先的时间序列中得到了一个向量序列X＝[X₁，X₂，...，X_N-(M-1)J]，这个向量序列可以看作M维空间的一个轨迹，同时点集合{X₁，X₂，...，X_N-(M-1)J}可看作M维空间的一个流形；

步骤3：这一步用主元分析法对点集合{X₁，X₂，...，X_N-(M-1)J}进行方向归一化并求出主轴，具体如下：首先计算矩阵XX^T的M个特征值(这里的特征值指矩阵的特征值，英文称作“Eigenvalue”)，假设XX^T的M个特征值按照从大到小的顺序排列为λ₁≥λ₂≥...≥λ_M且其对应的特征向量(这里指矩阵的特征向量，英文称作“Eigenvector”)依次分别为U₁，U₂，...，U_M，这里U_i:i∈{1，2，...，M}是M维的列向量，{U₁，U₂，...，U_M}也称为主轴，M个主轴在M维空间张成一个新的坐标系，令矩阵U＝[U₁，U₂，...，U_M]，则原来坐标系中M维列向量X_j在新坐标系中的值为

$Y_j=\left[\begin{array}{c}{Y}_{j1}\\ Y_{j2}\\ \·\\ \·\\ \·\\ Y_{\mathrm{jM}}\end{array}\right]=\left[\begin{array}{c}{U}_1^T{X}_j\\ U_2^T{X}_j\\ \·\\ \·\\ \·\\ U_M^T{X}_j\end{array}\right]$

这里j∈{1，2，...，M}，如此可以求出原先坐标系中所有M维向量X₁，X₂，...，X_N-(M-1)J在新坐标系中的坐标值，X₁，X₂，...，X_N-(M-1)J经过坐标变换后在新坐标系中依次对应Y₁，Y₂，...，Y_N-(M-1)J；

步骤4：这一步对点集合{Y₁，Y₂，...，Y_N-(M-1)J}进行位置归一化，具体如下：首先计算点集合{Y₁，Y₂，...，Y_N-(M-1)J}的阶数为[p₁，p₂，...，P_M]的矩

$m_{p_1{p}_2...p_M}=\underset{j=1}{\overset{N-(M-1)J}{\mathrm{\Σ}}}\underset{i=1}{\overset{M}{\mathrm{\Π}}}{Y}_{\mathrm{ji}}^{p_i}$

当{p_i|i＝1，2，...，M}的取值满足 $p_i=\left\{\begin{array}{cc}1& i=k\\ 0& i\≠k\end{array}\right.$ 时，记 ${\stackrel{\‾}{m}}_k=m_{p_1{p}_2...p_M},$ 当{p_i＝0|i＝1，2，...，M}时，记 ${\stackrel{\‾}{m}}_0=m_{p_1{p}_2...p_M},$ 当k＝1，2，...，M时，计算

${\stackrel{\‾}{Y}}_k=\frac{{\stackrel{\‾}{m}}_k}{{\stackrel{\‾}{m}}_0}$

点集合{Y₁，Y₂，...，Y_N-(M-1)J}的中心由下式定义：

$\stackrel{\‾}{Y}=\left[\begin{array}{c}{\stackrel{\‾}{Y}}_1\\ {\stackrel{\‾}{Y}}_2\\ \·\\ \·\\ \·\\ {\stackrel{\‾}{Y}}_M\end{array}\right]$

对点集合{Y₁，Y₂，...，Y_N-(M-1)J}进行平移操作，将中心

移到坐标原点，具体操作方法是：将{Y₁，Y₂，...，Y_N-(M-1)J}中的每个向量减去即令

$Y_j\←Y_j-\stackrel{\‾}{Y}$

步骤5：这一步截取轨迹的局部流形，具体如下：可以构造一个超平面与第i个主轴U_i垂直相交于某个位置，上述超平面也称为Poincare截面，假设Y₀是这个Poincare截面与U_i轴的交点在U_i轴上的坐标值，计算轨迹[Y₁，Y₂，...，Y_N-(M-1)J]与这个Poincare截面的交点，计算方法如下：如果向量Y_j的第个坐标值满足Y_ji＝Y₀，则Y_j是一个交点；如果(Y_ji-Y₀)(Y_(j+1)i-Y₀)＜0，则利用线性插值求出一个交点

$Y_j^{\′}=Y_j+\frac{Y_{\mathrm{ji}}-Y_0}{Y_{\mathrm{ji}}-Y_{(j+1)i}}(Y_{j+1}-Y_j);$

可以按上述方式构造多个Poincare截面，在每个Poincare截面上可以求出轨迹与这个Poincare截面的交点的集合，轨迹与每个Poincare截面的交点的集合即为一个局部流形；

步骤6：令P(M，i，Y₀)表示与第i个主轴U_i垂直相交于Y₀的一个Poincare截面，P(M，i，Y₀)中的“M”表示M维空间，令C(M，i，Y₀)表示轨迹[Y₁，Y₂，...，Y_N-(M-1)J]穿越Poincare截面P(M，i，Y₀)的次数，C(M，i，Y₀)就是轨迹与Poincare截面P(M，i，Y₀)的交点的个数，令Z(M，i，Y₀)＝C(M，i，Y₀)/[N-(M-1)J]，Z(M，i，Y₀)表示轨迹穿越Poincare截面P(M，i，Y₀)的次数与向量总数[N-(M-1)J]的比值，Z(M，i，Y₀)就是所求的脑电信号的一个特征，由于可以构造不止一个Poincare截面(在M维空间中，当i和Y₀取不同值时可以构造不同的Poincare截面)，所以可以得到一个特征的集合{Z(M，i，Y₀)}，{Z(M，i，Y₀)}就是所求的脑电信号的特征，这些参量可以刻画正常、病症不活跃期、病症活跃期等不同类别脑电信号之间的区别，如果将{Z(M，i，Y₀)}中的元素按一定次序排列就构成脑电信号的一个特征向量；

步骤7：将步骤6中得到的脑电信号的特征向量输入分类器，这里使用支持向量机作为分类器，支持向量机可以用开放源代码软件LIBSVM实现，源代码见http://www.csie.ntu.edu.tw/～cjlin/libsvm，对于输入的每一个脑电信号的特征向量，支持向量机会输出一个类别号，这个类别号指示着步骤1中采集到的脑电信号的类别(即：正常、病症不活跃期、病症活跃期)；支持向量机的参数需要经过训练调整到较佳值后，在分类时才能较好地工作，LIBSVM软件包里提供有训练工具，可以完成支持向量机参数的训练。

脑电信号自动评估实施例2：

步骤1：这一步完成脑电信号的采集，将电极帽连接到放大器、将放大器通过USB接口连接到计算机系统、调用驱动程序将脑电信号采集进计算机；

步骤2：这一步对输入时间序列S＝[S₁，S₂，...，S_N]进行状态空间重构，选定两个参数M和J，取J为[S₁，S₂，...，S_N]的自相关函数第一次下降到最大值的1/e时对应的时刻，令M＝20，令t(i)＝1+2(i-1)，通过状态空间重构可以得到

$X=[X_1,X_2,\·\·\·X_L]=\left[\begin{array}{cccc}{S}_1& S_3& \·\·\·& S_{1+2(L-1)}\\ S_{1+J}& S_{3+J}& \·\·\·& S_{1+2(L-1)+J}\\ \·& \·& \·& \·\\ \·& \·& \·& \·\\ \·& \·& \·& \·\\ S_{1+(M-1)J}& S_{3+(M-1)J}& \·\·\·& S_N\end{array}\right]$

这里L＝1+[N1-(M-1))J/2；

步骤3：这一步对X＝[X₁，X₂，...，X_L]构成的点集合进行方向归一化，具体如下：首先计算{X₁，X₂，...，X_L}的算术平均中心 $\stackrel{\‾}{X}=\frac{1}{L}\underset{i=1}{\overset{L}{\mathrm{\Σ}}}{X}_i,$ X_i与两点可以确定一条直线L_i，设X_j到直线L_i距离为d_ij，求 $\{D_i=\frac{1}{L-1}\underset{j=1}{\overset{L}{\mathrm{\Σ}}}{d}_{\mathrm{ij}}|i=\mathrm{1,2},...,L\},$ 计算 $k=\arg \underset{i}{\min }\left\{D_i\right\},$ 选 $U_0=\frac{X_k-\stackrel{\‾}{X}}{\left|\right|X_k-\stackrel{\‾}{X}\left|\right|}$ 为第一主轴，构造一个经过中心

且与第一主轴U₀垂直的超平面，计算{X₁，X₂，...，X_L}在此超平面的投影{P₁，P₂，...，P_L}，对P＝[P₁，P₂，...，P_L]进行主元分析，即：找到矩阵PP^T的M-1个按从大到小顺序排列的特征值λ₁≥λ₂≥...≥λ_M-1和其对应的特征向量U₁，U₂，...，U_M-1，{U₀，U₁，...，U_M-1}构成一个新的坐标系，用U＝[U₀，U₁，...，U_M-1]对X＝[X₁，X₂，...，X_L]进行坐标变换，得到Y＝U^TX，Y＝[Y₁，Y₂，...，Y_L]对应经过方向归一化的流形。

步骤4：这一步对Y＝[Y₁，Y₂，...，Y_L]构成的流形进行位置归一化，具体如下：首先计算

$\stackrel{\‾}{Y}=\frac{1}{L}\underset{i=1}{\overset{L}{\mathrm{\Σ}}}{Y}_i$

对点集合{Y₁，Y₂，...，Y_L}进行平移操作，将中心

移到坐标原点，具体操作方法是：将{Y₁，Y₂，...，Y_L}中的每个向量减去

即令

$Y_j\←Y_j-\stackrel{\‾}{Y}$

步骤5：构造多对与各个主轴垂直的Poincare截面，假设一对与主轴U_i垂直的Poincare截面与U_i轴的交点在U_i轴上的坐标分别为Y₀₁和Y₀₂，求落入这两个Poincare截面之间的轨迹上的点，具体计算方法是：如果轨迹[Y₁，Y₂，...，Y_L]的第j个点Y_j＝[Y_j1，Y_j2，...，Y_jM]的第i个坐标值满足min{Y₀₁，Y₀₂}≤Y_ji≤max{Y₀₁，Y₀₂}，则判定Y_j为落入这两个Poincare截面之间的点；可以构造多对垂直于各个主轴的Poincare截面，再按照上述计算方法可以求出落入每对Poincare截面之间的点，落入每对Poincare截面之间的点构成一个点集合，这个点集合即为所求的局部流形；

步骤6：令{Z₁，Z₂，...，Z_K}表示通过计算步骤4截取的轨迹[Y₁，Y₂，...，Y_L]落入某个局部的点的集合，这里必然满足{Z₁，Z₂，...，Z_K}∈{Y₁，Y₂，...，Y_L}，点集合{Z₁，Z₂，...，Z_K}的阶数为[P₁，P₂，...，P_M]的矩为

$m_{P_1{P}_2...P_{M-1}}=\underset{j=1}{\overset{K}{\mathrm{\Σ}}}\underset{i=1}{\overset{M}{\mathrm{\Π}}}{Z}_{\mathrm{ji}}^{p_i}$

当{P_i|i＝1，2，...，M}的取值满足 $p_i=\left\{\begin{array}{cc}1& i=k\\ 0& i\≠k\end{array}\right.$ 时，记 ${\stackrel{\‾}{m}}_k=m_{P_1{P}_2...P_{M-1}},$ 当{P_i＝0|i＝1，2，...，M}时，记 ${\stackrel{\‾}{m}}_0=m_{P_1{P}_2...P_{M-1}},$ 计算

$\{{\stackrel{\‾}{Z}}_k=\frac{{\stackrel{\‾}{m}}_k}{{\stackrel{\‾}{m}}_0}|k=\mathrm{1,2},...,M\}$

点集合{Z₁，Z₂，...，Z_K}的阶数为[P₁，P₂，...，P_M]的中心矩定义为

${\mathrm{\μ}}_{P_1{P}_2...P_{M-1}}=\underset{j=1}{\overset{K}{\mathrm{\Σ}}}\underset{i=1}{\overset{M}{\mathrm{\Π}}}{(Z_{\mathrm{ji}}-{\stackrel{\‾}{Z}}_i)}^{P_i}$

由于可以构造多对Poincare截面截取多个局部流形，对于每个局部流形都可以通过上述计算方法求出中心矩，所求得的中心矩的值组成的集合可以作为脑电信号的特征，这个集合中的元素按一定次序排列就构成脑电信号的一个特征向量。

步骤7：将步骤6中得到的脑电信号的特征向量输入分类器，这里使用支持向量机作为分类器，支持向量机可以用开放源代码软件LIBSVM实现，源代码见http://www.csie.ntu.edu.tw/～cjlin/libsvm，对于输入的每一个脑电信号的特征向量，支持向量机会输出一个类别号，这个类别号指示着步骤1中采集到的脑电信号的类别(例如：正常、病症不活跃期、病症活跃期)；支持向量机的参数需要经过训练调整到较佳值后，在分类时才能较好地工作，LIBSVM软件包里提供有训练工具，可以完成支持向量机参数的训练。

基于实施例1的方法进行了如下实验：首先从5个健康人和5个癫痫病人脑部采集脑电信号，共采集了5组信号，每组信号100个样本，信号采样频率173.16Hz，每个信号样本长度为23.6秒；A组信号是5个健康人睁眼时采集的，采集信号时电极位置按照国际标准布置，A组的100个信号样本是从不同的电极采集的；B组100个信号样本是在5个健康人闭眼时采集的；D组100个信号样本是在5个癫痫病人未发病时从病灶处采集的；C组100个信号样本是在5个癫痫病人未发病时从与病灶相反的脑部区域位置采集的；E组100个信号样本是在5个癫痫病人发病时采集的。A和B组的200个信号样本都是从健康人脑部采集的，这里看作一类；C组和D组的200个信号样本都是在癫痫病人不发病时采集的，这里看作一类；D组100个信号样本是癫痫病人发病时采集的，看作一类。实验过程如下：脑电信号自动评估前需要对分类器的参数进行训练，每类各随机取50％的信号样本作为训练样本、其余50％的样本作为测试样本，用训练样本对分类器进行训练，训练完成后，对测试样本进行分类，并计算分类正确率；以上过程重复10次，每次都随机选取各类中50％的样本作为训练样本，将10次的分类正确率求平均就是实验结果，实验结果如下表所示。

A和B	C和D	E	三类平均
				99.9％	99.5％	98.4％	99.44％

实验结果说明本发明提出的技术方案具有如下优点：(1)对信号采集时电极的位置不敏感：C组和D组的信号样本是从病灶和与病灶相反的部位采集的，但是都能被正确识别；(2)对接受测试者睁眼还是闭眼不敏感：A组和B组分别是接受测试者睁眼和闭眼时采集的信号，但是都能被正确识别；(3)病人不发病时的脑电信号能被正确识别：C组和D组的信号样本都是在病人不发病时采集的。

Claims (2)

1.一种脑电信号特征提取方法，其特征在于包含以下几个计算步骤：

(a)对输入的脑电时间序列S＝[S₁，S₂，...，S_N]进行状态空间重构：先选取参数J和M，取5≤M≤30，取J为脑电时间序列自相关函数第一次下降到最大值的1/e时对应的时刻，再任意选取L个满足1≤t(1)＜t(2)＜...＜t(L)≤N-(M-1)J的正整数{t(i)|i＝1，2，...，L}，令X_i＝[S_t(i)，S_t(i)+J，...，S_t(i)+(M-1)J]^T，则称[X₁，X₂，...，X_L]或{X₁，X₂，...，X_L}为时间序列S的一个状态空间重构，把{X₁，X₂，...，X_L}看作M维空间的一个点集合或流形，把[X₁，X₂，...，X_L]看作M维空间的一个轨迹；以上状态空间重构或者用行向量的形式实现，即令X_i＝[S_t(i)，S_t(i)+J，...，S_t(i)+(M-1)J]，i＝1，2，...，L；

(b)对重构的流形进行方向归一化：根据流形在各个方向上分布的特点定出基准方向，参照基准方向对流形进行坐标变换；

(c)对重构的流形进行位置归一化：根据流形的形状和结构特点确定基准点，参照基准点对流形进行平移操作；

(d)截取重构的流形或轨迹的一些局部区域；

(e)对截取的局部流形或轨迹进行统计量的计算，把计算得到的统计量作为输入的脑电时间序列的特征。

2.根据权利要求1所述的方法，其特征在于计算步骤(b)中，采用主元分析法定出基准方向，令X＝[X₁，X₂，...，X_L]，X是一个M×L维的矩阵，计算XX^T的M个特征值，并按从大到小的顺序排列得到λ₁≥λ₂≥...≥λ_M，以这M个矩阵特征值对应的矩阵XX^T的特征向量{U₁，U₂，...，U_M}作为主轴，令U＝[U₁，U₂，...，U_M]，利用Y＝U^TX对X进行坐标变换，得到新的坐标Y，Y对应经过方向归一化的流形；或者：

首先计算{X₁，X₂，...，X_L}的算术平均中心

X_i与两点可以确定一条直线L_i，设X_j到直线L_i距离为d_ij，求

D_i是所有点到直线L_i的平均距离，计算D_k是{D_i|i＝1，2，...，L}中的最小值，选

为第一主轴，这里

表示向量

的模，构造一个经过中心且与第一主轴垂直U₀的超平面，计算{X₁，X₂，...，X_L}在此超平面的投影值{P₁，P₂，...，P_L}，这里{P_i|i＝1，2，...，L}是M-1维向量，对P＝[P₁，P₂，...，P_L]进行主元分析，找到矩阵PP^T的M-1个从大到小排序的特征值λ₁≥λ₂≥...≥λ_M-1和其对应的特征向量U₁，U₂，...，U_M-1，{U₀，U₁，...，U_M-1}构成一个新的坐标系，用U＝[U₀，U₁，...，U_M-1]对X＝[X₁，X₂，...，X_L]进行坐标变换，得到Y＝U^TX，Y对应经过方向归一化的流形。

Images