专利名称：一种非线性脑电信号分析方法以及装置的制作方法
技术领域：
本发明涉及一种脑电信号分析方法和装置来识别脑状态，特别涉及一种对脑电序列的预处理和排序熵的分析方法和装置。
背景技术：
脑电信号表达大脑神经细胞的电活动，直接反映脑组织的电生理活动和大脑的功能状态。研究表明，通过对脑电信号的深入分析，可以识别不同的大脑状态。当前，脑电信号分析的主要方法有时域分析、频域分析、时频分析、信息熵和非线性动力学等。概述如下(1)时域分析法直接从时域提取特征，是最早发展起来的方法，它直观性强，物理意义比较明确，因此目前仍有不少使用，主要用来直接提取波形特征，进行峰值检测、直方图分析、方差分析等。实验证明这些参数与脑状态存在着一定的相关性。然而，由于脑电信号测量的个体差异性的影响，EEG的幅值不能直接反映脑电活动，使这些传统的脑电时域参数在实际中很少应用。(2)采用快速傅立叶变换及在此基础上发展起来的频域分析，主要包括δ率、中心频率(Medium power frequency, MPF)、边缘频率(Spectral edge frequency, SEF)等参数。S率是指EEG功率谱中δ波段(0.5-3. 75Hz)功率与α波段(8. 0-13. 5Hz)和β波段(13. 75-30. OHz)功率之和的比值，即慢波功率与快波功率之比，MPF和SEF分别是EEG功率谱中达到总功率的50%和95%时对应的频率。大量研究证实MF和SEF与脑状态有密切关系；然而，这些参数同样在不同个体间存在差异性。(3)脑电信号是时变、非平稳信号，单纯的时域或频域表示都不能准确的表达信号，只有把时间和频率结合起来进行处理，才能更好的从脑电信号中提取有用信息。时频方法在时域和频域同时具有良好的局部化分析以及其它一些重要性质，所以在脑电信号分析中受到重视。常用的时频方法有短时傅里叶变换、小波变换等，它们可以更好的揭示脑电信号的时频分布特性，进而研究大脑神经振荡的相关动力学特性，但如何从中提取有用的脑电特征以描述脑状态的变换，尚需进一步的研究。(4)熵是系统无序性(即混乱度)的一种度量。aiarmon熵根据已知幅值的概率分布预测将来值，计算直接方便，但它未对EEG总功率取归一化，这样也存在个体间差异性问题。傅里叶谱熵通过对信号计算功率谱得到，芬兰Datex-Ohmeda公司基于时频均衡谱熵算法开发了商业化M-entropy模块。然而，谱熵计算基于傅立叶变换，它更适合于分析线性平稳信号，而并不能准确的估计脑电信号的瞬时特征。小波熵将小波分析和信息熵结合，其基本思想是把小波变换系数矩阵处理成一个概率分布序列，用该序列的熵值来反映这个矩阵的稀疏程度，即被分析信号概率分布的有序程度。虽然小波熵算法可以得到信号复杂度在时域上的动态变化信息，但熵的概念是建立在概率的基础上，如果分析的信号较短，则没有统计意义。因此，小波熵在反映信号时变特性方面是比较粗糙的，只能做定性的研究。(5)近年来，随着非线性动力学的发展，越来越多的证据表明大脑是一个非线性系统，而脑电信号即是这个系统的输出。因此，人们把一些非线性动力学的方法应用于脑电信号分析，如关联维数、Lempel-Ziv复杂度、近似熵等。其中近似熵是一种常用的度量序列的复杂性方法，其特点是具有较好的抗干扰和抗噪的能力，对随机信号或是确定性信号都可使用，因此适合用来分析生物信号。但是近似熵的计算则需要相对长的、稳态的、无噪的数据，因此不适合用于实际的脑电状态监控系统的需求。通过分析脑电信号，研究大脑的状态变化仍然是起步阶段。时域和频域方法的应用尽管能看到显著的结果，但准确率不高，另外这些方法应用的前提是假设EEG为线性平稳信号，而实际EEG信号是来自大脑这一非线性动力学系统的非线性、非平稳信号，因此信息熵和非线性动力学方法得到了迅速发展，取得了较好的效果，然而这些分析方法还需要大量的研究来验证其有效性。
发明内容
本发明的目的是提供一种非线性脑电信号分析方法和装置，能够更为准确的对清醒和睡眠状态下的脑电信号分析。本发明的技术方案如下一种非线性脑电信号分析方法，包括以下步骤步骤101，实时采集原始脑电信号；步骤102，对采集的脑电信号进行分段；步骤103，分析各段脑电信号状态并做预处理去除干扰；步骤104，计算各段脑电信号的排序熵，根据排序熵值的大小确定脑的活动状态。步骤102种采用移动窗口技术对脑电信号分段。所述脑电信号分段的分段长度选取10s，相邻数据段重叠75%。步骤103中分析各段脑电信号状态，如果脑电信号为清醒脑电信号，则对清醒脑电信号的预处理采取以下措施去除奇异信号；根据信号的统计特性确定阈值去除异常信号点或段；去除眼动干扰；利用眼动干扰信号的时域、频域特征检测其是否存在，存在的话应用平稳小波变换(SWT)设定阈值去除；去除白噪声、基线漂移以及呼吸伪差；利用小波变换方法分解脑电信号，将最低频带系数置零以去除由于电极阻抗改变引起的基线漂移及呼吸伪差的影响，在其它频带上采用基于经验贝叶斯估计的阈值化方法去除信号中含有的高斯白噪声；去除肌电等干扰；用卡尔曼滤波器训练自适应AR模型去除肌电信号以及其它瞬态高幅信号干扰；如果脑电信号为睡眠脑电信号，则对睡眠脑电信号的预处理采取以下措施去除奇异信号；根据信号的统计特性确定阈值去除异常信号点或段；去除白噪声、基线漂移以及呼吸伪差；利用小波变换方法分解脑电信号，将最低频带系数置零以去除由于电极阻抗改变引起的基线漂移及呼吸伪差的影响，在其它频带上采用基于经验贝叶斯估计的阈值化方法去除信号中含有的高斯白噪声。步骤104中脑电信号的排序熵计算采用以下方法①将脑电信号|x(i)，i = 1，2，···}嵌入到一个m维空间中=Xi = [x(i)， x(i+L),…，x(i+(m-l)L)]，其中m为嵌入维数，L为延迟；②将序列Xi中的元素按照增序排列[x(i+a_l)L) ^ χ(i+(J2-I)L)彡…彡x(i+(jm-l)L)]；当序列中有元素相等，即x(i+(jn_l)L) = χ(i+(Ji2-DL)时，χ按照相应的j的顺序排序，即如果jn < ji2， x(i+(jn-l)U < χ(i+(Ji2-DL)；③将这个重构序列看作一种符号序列，设每种符号出现的概率分别为P1,P2, ···&，其中K彡m!，那么时间序列{x(i)，i = 1，2，···}的排序熵定义为
KHp(Iit) = -^PjInPj
；=1当Pj = 1/m ！时，Hp(m)取得最大值In (m ！)，据此将Hp(m)规一化，表示为0 ^ Hp = Hp/In (m ！)彡 1最终输出排序熵值和绘制排序熵图，根据排序熵值的大小估确定脑活动状态。一种非线性脑电信号分析装置，包括用于实时采集原始脑电信号的采集模块；用于对脑电信号分段的分段模块；用于分析各段脑电信号状态并做预处理去除干扰的预处理模块；用于计算各段脑电信号的排序熵，根据排序熵值的大小确定脑的活动状态的计算确定模块。所述分段模块采用移动窗口技术对脑电信号进行分段。所述脑电信号分段的分段长度选取10s，相邻数据段重叠75%。所述预处理模块分析各段脑电信号状态，如果脑电信号为清醒脑电信号，则对清醒脑电信号进行如下预处理根据信号的统计特性确定阈值去除异常信号点或段；利用眼动干扰信号的时域、频域特征检测其是否存在，存在的话应用平稳小波变换(SWT)设定阈值去除； 利用小波变换方法分解脑电信号，将最低频带系数置零以去除由于电极阻抗改变引起的基线漂移及呼吸伪差的影响，在其它频带上采用基于经验贝叶斯估计的阈值化方法去除信号中含有的高斯白噪声；用卡尔曼滤波器训练自适应AR模型去除肌电信号以及其它瞬态高幅信号干扰；如果脑电信号为睡眠脑电信号，则对睡眠脑电信号进行如下预处理去除奇异信号；根据信号的统计特性确定阈值去除异常信号点或段；去除白噪声、基线漂移以及呼吸伪差；利用小波变换方法分解脑电信号，将最低频带系数置零以去除由于电极阻抗改变引起的基线漂移及呼吸伪差的影响，在其它频带上采用基于经验贝叶斯估计的阈值化方法去除信号中含有的高斯白噪声。所述计算确定模块对脑电信号的排序熵计算为①将脑电信号|x(i)，i = 1，2，···}嵌入到一个m维空间中=Xi = [x(i)， x(i+L),…，x(i+(m-l)L)]，其中m为嵌入维数，L为延迟；②将序列Xi中的元素按照增序排列[x(i+a_l)L) ^ χ(i+(J2-I)L) x(i+(jm-l)L)]；当序列中有元素相等， 即x(i+(jn_l)L) = χ(i+(Ji2-DL)时，χ按照相应的j的顺序排序，即如果jn < ji2， x(i+(jn-l)U < χ(i+(Ji2-DL)；③将这个重构序列看作一种符号序列，设每种符号出现的概率分别为P1,P2, ···&，其中K彡m!，那么时间序列{x(i)，i = 1，2，···}的排序熵定义为
Hp(Tn) = -YjPjXrJyi
；=1当Pj = 1/m ！时，Hp(m)取得最大值In (m ！)，据此将Hp(m)规一化，表示为0 ^ Hp = Hp/In (m ！)彡 1最终输出排序熵值和绘制排序熵图，根据排序熵值的大小估确定脑活动状态。本发明可用于大脑不同认知功能活动、生理状态情况下的脑电信号分析，尤其可以应用于研究睡眠脑活动的变化等。本发明的方法和装置与现有方法相比，有益效果(1)排序熵从相空间的角度出发，对脑电信号的结构变化的模式进行统计分析，揭示脑电信号的不确定性、稳定程度和信息量。(2)本发明是从脑电信号的复杂性特征出发进行分析的，该方法在本质上与脑电信号的特点吻合。(3)排序熵将脑电序列映射为符号序列并恰当的表示和分析，概念简单、计算速度快，抗噪能力强。
图1为本发明的方法工作流程示意图。图2A为从慢波睡眠状态(SWS)到快速眼动睡眠状态(REM)转换期间的原始脑电信号。图2B为图2A的脑电信号预处理后计算的排序熵。图3A为快速眼动睡眠状态(REM)到清醒状态(AWK)转换期间的原始脑电信号。图;3B为图3A的脑电信号预处理后计算的排序熵。图4A是在清醒状态下采集的一段长10秒的含噪脑电信号。图4B是应用步骤3中所述的预处理方法对图4A的脑电信号处理的结果。图 5 为 m = 3 和 L = 1，2，3，4，5(即 PE(3，1)，PE (3，2)，PE (3，3)，PE (3，4)，PE (3， 5))时，SWS和REM状态下的排序熵的比较结果(左列)，以及REM和AWK状态下的排序熵的比较结果(右列)。图6是本发明分析装置的结构框图。
具体实施例方式下面结合附图和一个具体实施例对本发明作进一步的详细说明。图1为本发明的方法工作流程示意图。首先是步骤101，采集脑电信号。本实施例对4只大鼠进行了脑电信号采集，采样频率是128Hz。图2A为一例从慢波睡眠状态(SWS) 到快速眼动睡眠状态(REM)转换期间的脑电信号。图3A为一例快速眼动睡眠状态(REM) 到清醒状态(AWK)转换期间的脑电信号。在步骤102，采用移动窗口技术对采集的脑电信号数据进行分段处理。为了实时跟踪大脑状态的变化，需要对脑电信号进行分段处理。在此实施例中， 对上述脑电数据分段，分段长度选取10s，即N = 1观0(采样频率U8Hz)，相邻数据段重叠 75%进行数据分析。
在步骤103，对分段脑电信号进行预处理。预处理前须判断脑电信号所处的状态。通过脑电中肌电分量的存在性来判断大鼠处于清醒还是意识消失的睡眠状态。肌电分量的出现表明大鼠可对外部激励产生响应，清醒时这种响应会比较显著，而进入睡眠后这种响应会减弱。本实施例借助时频均衡谱熵算法来实现失意时刻的检测。该算法分别测定了反应熵与状态熵，两熵之差为肌电熵即反映了肌电活动。通过监测脑电信号中肌电熵的变化，当熵值低于某一阈值时，即判断此时大鼠由清醒转入失意状态，然后在此前后将采用不同的预处理方案。若大鼠处于清醒期，则对脑电信号采取如下的预处理措施①去除奇异信号；具体可根据信号的统计特性，如均值、方差等确定阈值去除异常信号点或段；②去除眼动干扰；眼动干扰信号表现为高幅慢波，其频率范围一般在2-16HZ。先利用时频域的一些特征检测是否存在眼动干扰，再应用平稳小波变换(SWT)设定阈值去除眼动干扰；③去除白噪声和基线漂移以及呼吸伪差；基于小波变换的方法将信号分解到不同频带，将最低频带 (0-0. 8Hz)系数置零以去除由于电极阻抗改变引起的基线漂移以及呼吸伪差，在其它频带上采用基于经验贝叶斯估计的阈值化方法，去除脑电信号中可能含有的高斯白噪声；④去除肌电等干扰；肌电信号通常表现为高频瞬态大幅度特征，频率范围为32-47HZ，用卡尔曼滤波器训练的自适应的AR模型可有效去除肌电信号以及其它瞬态大幅度信号干扰。若大鼠脑电信号处于睡眠状态，仅采用清醒时预处理措施的第①和③措施，即无需检测和去除眼动和肌电干扰。图4A是在清醒状态下采集的一段长10秒的含噪脑电信号，图4B是应用步骤103 中所述的预处理方法对图4A的脑电信号处理的结果。在步骤104，通过计算排序熵来描述脑状态，其具体细节如下将每一段脑电信号{x(i), i = 1,2, ···}嵌入到一个m维空间中=Xi = [x(i), x(i+L), ...，x(i+(m-l)L)]。其中m为嵌入维数，L为延迟。对序列Xi中的元素按照增序排列[x(i+a-l)L) ^ χ(i+(J2-I)L)彡… ^ χ (i+(Jffl-DL) J0 当有数相等，即 x(i+(jn-l)L) = χ (i+(Ji2-I)L)时，χ 按照相应的 j 的顺序排序，即如果jn < ji2，那么χ (i+ (Jil-I) U <x (i+ (ji2-l) L)。那么任意向量\都能唯一的映射为(J: J2,…丄)。对于m个符号(1，2，…m)可能最多有m ！种排列方式。显然，在m维空间中的每个向量&，都能被映射为m ！种排列方式中的一种。令每种排列方式出现的概率为P1, P2, ···&，其中K彡m!。那么时间序列|x(i)，i = 1，2，···}的排序熵为:
KHp(Iit) = -^PjInPj
；=1这里举一个例子来说明排序熵的计算。如χ = {4，7，9，10，6，11，3}，取111 = 3丄=
1，可重构序列到相空间 X =
47979109106106116113
将其每一行中的数据按增序进行排列并映射为符号，则各行依次可表示为012， 012，120，102，120，那么共有5种排列，且排列号012的次数为2次，排列号120的次数为2
次，排列号102的次数为1次，因此排序熵为
权利要求
1.一种非线性脑电信号分析方法，其特征在于，包括以下步骤 步骤101，实时采集原始脑电信号；步骤102，对采集的脑电信号进行分段；步骤103，分析各段脑电信号状态并做预处理去除干扰；步骤104，计算各段脑电信号的排序熵，根据排序熵值的大小确定脑的活动状态。
2.根据权利要求1所述的一种非线性脑电信号分析方法，其特征在于，步骤102种采用移动窗口技术对脑电信号分段。
3.根据权利要求2所述的一种非线性脑电信号分析方法，其特征在于，所述脑电信号分段的分段长度选取10s，相邻数据段重叠75%。
4.根据权利要求1或2所述的一种非线性脑电信号分析方法，其特征在于，步骤103中分析各段脑电信号状态，如果脑电信号为清醒脑电信号，则对清醒脑电信号的预处理采取以下措施去除奇异信号；根据信号的统计特性确定阈值去除异常信号点或段； 去除眼动干扰；利用眼动干扰信号的时域、频域特征检测其是否存在，存在的话应用平稳小波变换设定阈值去除；去除白噪声、基线漂移以及呼吸伪差；利用小波变换方法分解脑电信号，将最低频带系数置零以去除由于电极阻抗改变引起的基线漂移及呼吸伪差的影响，在其它频带上采用基于经验贝叶斯估计的阈值化方法去除信号中含有的高斯白噪声；去除肌电干扰；用卡尔曼滤波器训练自适应AR模型去除肌电信号以及瞬态高幅信号干扰；如果脑电信号为睡眠脑电信号，则对睡眠脑电信号的预处理采取以下措施 去除奇异信号；根据信号的统计特性确定阈值去除异常信号点或段； 去除白噪声、基线漂移以及呼吸伪差；利用小波变换方法分解脑电信号，将最低频带系数置零以去除由于电极阻抗改变引起的基线漂移及呼吸伪差的影响，在其它频带上采用基于经验贝叶斯估计的阈值化方法去除信号中含有的高斯白噪声。
5.根据权利要求4所述的一种非线性脑电信号分析方法，其特征在于，步骤104中脑电信号的排序熵计算采用以下方法①将脑电信号{x(i)，i = 1,2, ···}嵌入到一个m维空间中=Xi = [x(i), x(i+L),…，x(i+(m-l)L)]，其中m为嵌入维数，L为延迟；②将序列Xi中的元素按照增序排列[x(i+a_l)L) ^ χ(i+(J2-I)L) x(i+(jm-l)L)]；当序列中有元素相等， 即x(i+(jn_l)L) = χ(i+(Ji2-DL)时，χ按照相应的j的顺序排序，即如果jn < ji2， x(i+(jn-l)U < χ(i+(Ji2-DL)；③将这个重构序列看作一种符号序列，设每种符号出现的概率分别为P1,P2, ···&，其中K彡m!，那么时间序列{x(i)，i = 1，2，···}的排序熵定义为
6.一种非线性脑电信号分析装置，其特征在于，包括 用于实时采集原始脑电信号的采集模块；用于对脑电信号分段的分段模块；用于分析各段脑电信号状态并做预处理去除干扰的预处理模块； 用于计算各段脑电信号的排序熵，根据排序熵值的大小确定脑的活动状态的计算确定模块。
7.根据权利要求6所述的一种非线性脑电信号分析装置，其特征在于，所述分段模块采用移动窗口技术对脑电信号进行分段。
8.根据权利要求7所述的一种非线性脑电信号分析装置，其特征在于，所述脑电信号分段的分段长度选取10s，相邻数据段重叠75%。
9.根据权利要求6或7所述的一种非线性脑电信号分析装置，其特征在于，所述预处理模块分析各段脑电信号状态，如果脑电信号为清醒脑电信号，则对清醒脑电信号进行如下预处理根据信号的统计特性确定阈值去除异常信号点或段；利用眼动干扰信号的时域、频域特征检测其是否存在，存在的话应用平稳小波变换 (SWT)设定阈值去除；利用小波变换方法分解脑电信号，将最低频带系数置零以去除由于电极阻抗改变引起的基线漂移及呼吸伪差的影响，在其它频带上采用基于经验贝叶斯估计的阈值化方法去除信号中含有的高斯白噪声；用卡尔曼滤波器训练自适应AR模型去除肌电信号以及瞬态高幅信号干扰； 如果脑电信号为睡眠脑电信号，则对睡眠脑电信号进行如下预处理 去除奇异信号；根据信号的统计特性确定阈值去除异常信号点或段； 去除白噪声、基线漂移以及呼吸伪差；利用小波变换方法分解脑电信号，将最低频带系数置零以去除由于电极阻抗改变引起的基线漂移及呼吸伪差的影响，在其它频带上采用基于经验贝叶斯估计的阈值化方法去除信号中含有的高斯白噪声。
10.根据权利要求9所述的一种非线性脑电信号分析装置，其特征在于，所述计算确定模块对脑电信号的排序熵计算为①将脑电信号{x(i)，i = 1,2, ···}嵌入到一个m维空间中=Xi = [x(i), x(i+L),…，x(i+(m-l)L)]，其中m为嵌入维数，L为延迟；②将序列Xi中的元素按照增序排列[x(i+a_l)L) ^ χ(i+(J2-I)L) x(i+(jm-l)L)]；当序列中有元素相等， 即x(i+(jn_l)L) = χ(i+(Ji2-DL)时，χ按照相应的j的顺序排序，即如果jn < ji2， x(i+(jn-l)U < χ(i+(Ji2-DL)；③将这个重构序列看作一种符号序列，设每种符号出现的概率分别为P1,P2, ···&，其中K彡m!，那么时间序列{x(i)，i = 1，2，···}的排序熵定义为Hp(Tn) = -Y4PjXnPj；=1当P」=1/m ！时，Hp (m)取得最大值ln(m ！)，据此将Hp (m)规一化，表示为 0 ^ Hp = Hp/In (m ！)彡 1最终输出排序熵值和绘制排序熵图，根据排序熵值的大小估确定脑活动状态。
全文摘要
本发明公开一种脑电信号分析方法和装置来识别脑状态，对采集的脑电信号进行分段，分析各段脑电信号状态并做预处理去除干扰；最后计算各段脑电信号的排序熵，根据排序熵值的大小确定脑的活动状态。本发明可用于大脑不同认知功能活动、生理状态情况下的脑电信号分析，尤其可以应用于研究睡眠脑活动的变化等。本发明采用的排序熵从相空间的角度出发，对脑电信号的结构变化的模式进行统计分析，揭示脑电信号的不确定性、稳定程度和信息量。
文档编号A61B5/0476GK102488516SQ201110414388
公开日2012年6月13日 申请日期2011年12月13日 优先权日2011年12月13日
发明者李小俚, 李段, 欧阳高翔 申请人:湖州康普医疗器械科技有限公司