CN102274022B - 一种基于脑电信号的睡眠状态监测方法 - Google Patents

Abstract

本发明公开了一种基于脑电信号的睡眠状态监测方法，可以实现连续的睡眠深度监测和睡眠阶段分级。此方法通过对睡眠脑电信号采用传统的时频域分析方法和新引入的非线性动力学分析方法提取到一批能表征睡眠深度的特征参数。再通过套索回归方法选择特征参数的最佳子集，通过统计学验证方法，确定模型系数。最终建立睡眠指数自动监测模型，可以根据相应的睡眠阶段分级标准实现睡眠阶段的自动分级。本发明更加准确地反映脑电活动和睡眠状态之间的关系，便于后续采用相关生理信号的监测来改善睡眠质量和睡眠监护。

Description

一种基于脑电信号的睡眠状态监测方法

技术领域

本发明涉及信号处理技术领域，尤其涉及一种在睡眠状态下的脑电信号(EEG)的特征提取和信息融合方法。

背景技术

睡眠作为一种复杂的生理过程，是机体复原、整合和巩固的重要环节。目前国际上普遍使用R&K睡眠EEG分期规则，根据睡眠时脑电信号的表现，将睡眠分为:觉醒期、非快速眼动睡眠期(non-rapid eye movement，NREM)和快速眼动睡眠期(rapid eye movement，REM)。其中NREM又可分为1、2、3、4期。睡眠分期研究，对于睡眠状态分析以及睡眠质量的科学评估，具有重要的应用价值。  

脑电数字化分析已广泛用于大脑皮层的功能监测，应用单一的数字化EEG 参数分析大脑在睡眠状态下的镇静程度更为直观方便。脑电信号（EEG）反映脑细胞群自发而有节律的电活动，一般用波幅、频率和相位等特征来描述。当进入睡眠状态时，EEG的基本特征随着睡眠程度加深或变浅呈顺序性变化，呈现一定的函数关系，因而可用来反映睡眠深度。

在已有的临床睡眠分析中，通常采用时频域分析的方法，提取时域、频域和高阶谱变量等特征参数作为脑电图定量分析指标，其通过特定的算法建立睡眠深度量化模型。

根据现有脑神经生理学EEG产生机制的研究，EEG信号起源于一个高度的非线性系统，不仅在中枢神经系统每个分层发现许多的反馈环路，而且单个神经元自身也表现出高度非线性因素。在神经细胞膜上可以观察到混沌行为，神经放电转化遵循分叉规律，而混沌和分叉行为属于非线性科学的范畴。因此EEG信号是大量神经细胞的非线性耦合，是一个高度非线性多单元连接的复合体，EEG活动具有确定性混沌特性，大脑是复杂、自组织的非线性动力学系统。

在分析脑电信号时，非线性方法较于传统的时频域分析方法，有着独特的优势。因此，本发明中不仅采用传统的时频域分析的方法提取到麻醉特征参数，包括EEG四个节律的比重、总功率谱重心 、β比率 ，还从非线性动力学的角度刻画麻醉状态的EEG变化规律，得到基本的非线性参数，包括近似熵、复杂度、信息熵；另外还有非线性混沌参数，包括最小嵌入维、关联维数、最大李雅普洛夫指数。更加准确地反映脑电活动和睡眠状态之间的关系。

发明内容

本发明的目的是针对现有技术的不足，提供一种基于脑电信号的睡眠状态监测方法。

本发明的目的是通过以下技术方案来实现的：一种基于脑电信号的睡眠状态监测方法，该方法包括以下步骤：

（1）结合多导睡眠监护仪，获取原始脑电信号：用多导睡眠监护仪记录监护者从入睡到醒来的连续脑电数据，采样频率为256kz；采集得到睡眠脑电样本，记录被试者由清醒到深度睡眠整个过程中的脑电信号；采集不同性别、不同年龄层的实验者的数据，建立睡眠脑电数据仓库；

（2）采用时频域方法和非线性动力学分析方法对睡眠脑电数据仓库中的睡眠脑电信号进行睡眠深度特征参数提取，待提取的特征参数包括时频域参数、非线性动力学参数和非线性混沌参数；

（3）对提取得到的特征参数做信息融合处理，以期得到能反映麻醉程度指标的量化方法；采用统计学的方法，在大量样本集合的验证下，确定最优模型来量化睡眠深度指数计算模型；该模型的输入为样本的睡眠深度特征参数，输出为0到100的连续无标度睡眠深度指数；

（4）根据上步计算得到的连续的睡眠深度指数，按照R&K睡眠EEG分期规则，根据睡眠深度指数与睡眠状态的关系进行睡眠阶段的分级。

本发明的有益效果是：现有的多导睡眠分析系统在基于脑电信号（EEG）分析睡眠结构时，通常借助眼电图（EOG）和肌电图（EMG）和一些其他的生理参数。这为睡眠状态监测的灵活性带来了一定的影响，本发明只通过对脑电信号分析提取更有效果的特征参数，来进行睡眠状态的判别，输出睡眠深度指数（1-100的无标度数值）。睡眠深度指数可以实时反映监护者的睡眠状态，预测睡眠深度趋势。这一功能可以加强现有睡眠分析软件系统。更加准确地反映脑电活动和睡眠状态之间的关系，便于后续采用相关生理信号的监测来改善睡眠质量，和睡眠监护。

附图说明

图1为本发明的自动睡眠状态监测系统结构图；

图2为本发明的近似熵算法流程图；

图3为本发明的信息熵算法流程图；

图4为本发明复杂度算法流程；

图5 为本发明的关联维的算法流程图；

图6 为本发明的睡眠脑电信号样本；

图7为本发明的连续睡眠深度指数计算结果图；

图8为本发明的标准脑电镇静程度曲线图；

图9为本发明的连续睡眠深度指数计算结果图；

图10为本发明的六等级睡眠阶段分级图。

具体实施方式

非线性动力学理论在分析脑电信号的关联维、Lyapunov指数、近似熵和复杂度等非线性参数上已有相关的理论研究基础，我们在此基础上，提取脑电信号的非线性特征。通过Cao方法等途径选取正确的延迟时间和嵌入维数，获得脑电信号高维空间中的混沌吸引子拓扑结构，分别进行其关联维、李亚普诺夫指数、近似熵、信息熵及复杂度等非线性参数的计算和分析。

在大量样本的统计学分析部分，采用 LASSO 套索估计方法揭示睡眠深度特征参数对睡眠深度指数影响的联合效应，并确定睡眠深度特征参数的最小子集来建立麻醉指数的评价模型。

LASSO 方法用于逼近线形模型。LASSO 最小化受约束的残留平方和误差。由于这一约束的本质特性，LASSO 可以使得一些系数严格为零，从而能够给出可解释的模型。

     

LASSO回归中的罚项

用

代替，使得解为y的非线性组合，计算时用二次规划算法。如果t选择的足够小，

相应地增大。会使得一些系数等于零。

本发明提出一种基于睡眠脑电信号的自动睡眠状态监测系统的设计方法：

该方法包括四个步骤，他们的顺序关系见图1。

在第一步中，结合现有的多导睡眠监护仪，获取原始脑电信号。采集不同性别、不同年龄层的实验者的数据，建立睡眠脑电数据仓库。

在第二步中，采用时频域方法和非线性动力学分析方法对睡眠脑电数据仓库中的睡眠脑电信号进行睡眠深度特征参数提取。待提取的特征参数包括时频域参数、非线性动力学参数、非线性混沌参数。

在第三步中，对提取得到的特征参数做信息融合处理，以期得到能反映麻醉程度指标的量化方法。采用统计学的方法，在大量样本集合的验证下，确定最优模型来量化睡眠深度指数计算模型。该模型的输入为样本的睡眠深度特征参数，输出为0到100的连续无标度睡眠深度指数。

在第四步中，根据上步计算得到的连续的睡眠深度指数，按照临床上采用的R&K睡眠EEG分期规则，根据睡眠深度指数与睡眠状态的关系进行睡眠阶段的分级。

第一步：脑电信号采集

用通用的多导睡眠监护仪记录监护者从入睡到醒来的连续脑电数据，采样频率为256kz。

采集得到睡眠脑电样本，记录被试者由清醒到深度睡眠整个过程中的脑电信号。

第二步：睡眠状态特性参数提取

任意截取连续的1000个点的睡眠脑电数据，按照下文中的计算方法提取三类脑电信号参数：频域参数、时域参数和非线性混沌参数。其中，频域参数为总功率和频率重心；时域参数包括近似熵、信息熵和复杂度；非线性混沌参数包括关联维和最大李亚普诺夫指数。

选定确定长度为1000个点作为构成滑动窗（滑动窗的长度可以根据实际测量数据的信噪比来调节）。顺次计算在整个睡眠脑电区间上的多个特征参数。得到的特征参数进入第三步的计算步骤。

总功率和频率重心的计算方法如下：

在脑电信号的临床和生理学研究中的关注的频率带为 0.5~35Hz，通常划分为四个频段：α 频段(7.5~12.5Hz)、β 频段(12.5Hz~35Hz)、δ 频段(0.5Hz~3.5Hz)、θ 频段(5~7.5Hz) 。

总功率：

每个节律所占总功率的比重计算如下：

;

是四个频段

所占总功率的比值。

、

分别为四个频段的上下范围。

是EEG信号的功率谱，是脑电信号的最大有效频率范围，这里取35Hz（可以根据实际测量装置的敏感性来调整）。

近似熵的计算方法如下：

近似熵的算法流程见图2。

1、对

点序列

，按序列顺序重构成一组

维矢量（为模式维数）：

2、定义

与

之间的距离

为两者对应元素中差值最大的一个，即：

3、给定阈值

，对每一个

值统计

小于的数目及此数目与距离总数的比值，记作

，即：

4、先将

取对数，再求其对所有

的平均值，记作

，即：

5、再把维数加1，变成

维，重复上述步骤，计算得

和

6、此序列的近似熵为：

一般来说，此极限以概率1存在。但实际计算中，往往以一定长度的序列作为

的估计值。记作

取

，

是原始序列的标准差。

信息熵的计算方法如下：

我们使用三种符号代表脑电信号的变化方向：

            

其中

表示脑电信号减小；

表示脑电信号基本持平；而

表示脑电信号增大。这三个值仅仅表示三种状态，其数值本身没有任何意义。该过程为符号化过程。

对于符号化的方向信号，为揭示其时序和结构的规律，使用滑动窗的方法构造字长宽度为

的向量序列：

譬如，当字长宽度为

时，连续变化共有种可能的模式。统计每种模式出现的概率：

，其中

是第

种模式出现的次数。

计算结果就是信息熵的值。

复杂度的计算方法如下：

第一步：对序列进行粗粒化处理，得到一个长度与原序列相同的符号序列。设已知的时间序列为

，记该时间序列的最大值和最小值分别为

、

。用表示将原序列中的数据粗粒化的段数，

和字符集，定义序列中的元素为

由此建立起一个具有

个不同字符的符号序列

。

第二步：扫描粗粒化处理后的重构序列，按照特定的算法计数以前没有出现过的新子串，所有子串的数目即为复杂度的绝对值；将此值用随机信号的复杂度值(

为序列长度)进行归一化，消除数据长度的影响，得到最后的复杂度。

关联维的计算方法如下：

关联维可以作为非线性系统复杂程度的量化说明，反映了原始信号的不同幅度的变化。关联维越高，系统越复杂。数学上严格的关联维求取过程需要关联积分曲线以及该曲线的微分曲线，同时还是个序列长度趋向无穷大距离向量趋向无限小的极限过程。为了解决关联维的计算机算法实现问题,Grassberg和Procaccia 提出了用于计算关联维数的GP算法。

关联积分

表示的意义为重构后的相空间中距离小于标度

的点对占所有点对的比例，在已知延时

和嵌入维数

时，标度的取值就直接影响到关联维

的值。当标度

取得过大，所有点对的距离都不会超过

，关联积分

=1，这样的

反应不了系统的动态特性；

取得过小，则几乎所有点对的距离都大于

，也反映不了系统的特性，因此标度

有一定的范围限制。标度

与关联积分之间有以下关系式成立：，变换得：

，我们做出关联积分

与标度的双对数曲线，即关联积分曲线图。通过最小二乘法拟合该直线，直线的斜率就是所求的关联维数。

最大李亚普诺夫指数LLE计算方法如下：

李亚普诺夫指数（Lyapunov指数）是衡量系统动力学特性的一个重要定量指标，它表征了系统在相空间中相邻轨道间收敛或发散的平均指数率。我们使用Wolf法计算。具体算法如下：

 对给定的混沌时间序列首先进行延时相空间重构，并搜索给定轨道上每个点（即嵌入矢量）的最近邻域点，即

，其中为时间序列平均周期，可以通过对时间序列的FFT估计得到。最大李亚普诺夫指数几何意义是量化初始闭轨道的指数发散和估计系统的总体混沌水平，据此有

，

两边取对数有：

，

由此可见，大致上等于上面这组直线的斜率，因此LLE可以通过最小二乘拟合该直线的“平均斜率”而得到，即

，

这里

表示对所有求平均。

第三步：多变量生理信号融合处理

选取第二步中对一段脑电数据样本计算得到的多个时频域特征参数和非线性混沌参数作为输入变量，临床评价的睡眠深度指数作为响应变量，对睡眠深度指数计算模型建模。

建模过程如下：

选取套索回归（lasso）来计算每个特征参数的权重系数。在该算法的实现中，采用一类解决凸优化的坐标下降（coordinate descent）方法来实现Lasso回归算法，该算法比传统算法更快速，更能适应融合套索回归分析（可参见相关文献）。

本发明采用glmnet软件包（开源，可以网上搜索得到），在MATLAB中导入软件包（ PACKAGE），来计算特征参数的权重系数。Lasso 和elastic-net（弹性网）回归产生线性模型。作为本软件的一个核心算法包之一，glmnet package 能够实现套索回归、弹性网回归、Logistic 回归和多元回归。该算法采用坐标下降方法，速度很快。而且在能通过对样本的交叉 验证智能选取最优模型，实现线性模型条件下的最精准预测。

首先要得到睡眠脑电和其他辅助监测设备得到的相关信息，医生可根据临床经验来确定受试者的睡眠状态，确定对应时间下的睡眠指数

（定义为1到100连续变化的离散值），作为本发明中睡眠指数计算模型的确立依据。

定义为上步中在第i个滑动窗中计算得到的特征参数的向量，

为在此滑动窗下医生根据临床确定的睡眠指数。

这里，

且。可以等价为在限制条件

下，

为可调参数，最小化

。使用上述软件包，计算得到最优解。确定各个参数的权重系数，即为睡眠深度指数模型。

通过上述步骤，得到的睡眠指数计算模型如下：

表1 睡眠深度指数计算模型系数表

   睡眠指数的计算模型可以写成以下线性模型：

上式中，

为表中所有参数的值构成的向量， 

为线性模型的截距，

是线性模型的系数（见上表）。

为由模型计算得到的睡眠深度指数。

对应原始脑电数据的第i个滑动窗。

在全部的脑电数据上采用上述方法中的滑动窗，计算全部区间上的特征参数

，带入上式，即得到整个监测时间区间上的睡眠深度指数值。为连续变化的0到100的数值曲线。

第四步：睡眠深度分级

国际上普遍使用R&K睡眠EEG分期规则，根据睡眠时脑电信号的表现，将睡眠分为:觉醒期、非快速眼动睡眠期(non-rapid eye movement，NREM)和快速眼动睡眠期(rapid eye movement，REM)。其中非快速眼动睡眠期又可分为NREM1、NREM2、NREM3、NREM4期。NREM4期对应最深度的睡眠状态。

表2：睡眠阶段特征表

应用R&K睡眠EEG分期规则，得到睡眠深度指数和睡眠阶段之间的关系（见上表），对应睡眠阶段按照从0到100的数字划分，0对应最深度睡眠，100代表完全清醒。可以反映出睡眠程度细微的变化过程，按照临床统计得到的睡眠深度指数和睡眠阶段的对应关系，可以灵活地调整相应睡眠阶段划分。也为预测睡眠程度变化提供了依据。

实施例1：

采集到原始睡眠脑电信号（图6），图中大致反映了睡眠脑电的幅度变化趋势。采样频率为256Hz，持续睡眠一个小时，睡眠程度由浅入深。经历清醒、浅睡眠阶段、中睡眠阶段和深睡眠阶段。每五秒计算一次睡眠深度指数。

用本发明中介绍的方法进行计算，得到连续的睡眠深度指数（图7）。

   根据睡眠深度指数，按照R&K睡眠EEG分期规则将睡眠状态划分为四个等级：1、清醒状态；2、浅睡眠状态（快速眼动）；3、中睡眠阶段；4、深睡眠阶段，见图8。

在对本实施例中用到的原始脑电信号采集的同时，我们还记录了大脑皮层活动程度曲线，用于表示被监护者的镇静程度（图9）。对比图7中我们计算得到的连续的睡眠深度指数，有很高的一致性。这说明我们的方法可以监测睡眠状态，同时还能反映大脑活动状态。

如果需要更为精细的睡眠阶段划分，根据睡眠深度指数，按照R&K睡眠EEG分期规则将睡眠状态划分为六个等级：觉醒期；快速眼动睡眠期（REM）； NREM1； NREM2； NREM3； NREM4。对应六等级睡眠阶段分级图中就是：2、觉醒期；3、快速眼动睡眠期（REM）；4、NREM1；5、NREM2；6、NREM3；由于NREM4数据在临床中较难得到，本例中没有NREM4阶段（图10）。

Claims (1)

1.一种基于脑电信号的睡眠状态监测方法，其特征在于，该方法包括以下步骤：

（1）结合多导睡眠监护仪，获取原始脑电信号：用多导睡眠监护仪记录监护者从入睡到醒来的连续脑电数据，采样频率为256kz；采集得到睡眠脑电样本，记录被试者由清醒到深度睡眠整个过程中的脑电信号；采集不同性别、不同年龄层的实验者的数据，建立睡眠脑电数据仓库；

（2）采用时频域方法和非线性动力学分析方法对睡眠脑电数据仓库中的睡眠脑电信号进行睡眠深度特征参数提取，待提取的特征参数包括时频域参数、非线性动力学参数和非线性混沌参数；

（3）对提取得到的特征参数做信息融合处理，以期得到能反映麻醉程度指标的量化方法；采用统计学的方法，在大量样本集合的验证下，确定最优模型来量化睡眠深度指数计算模型；该模型的输入为样本的睡眠深度特征参数，输出为0到100的连续无标度睡眠深度指数；

（4）根据上步计算得到的连续的睡眠深度指数，按照R&K睡眠EEG分期规则，根据睡眠深度指数与睡眠状态的关系进行睡眠阶段的分级；

所述步骤（2）中，所述频域参数为总功率和频率重心；时域参数包括近似熵、信息熵和复杂度；非线性混沌参数包括关联维和最大李亚普诺夫指数；

    所述步骤（3）中，所述模型通过以下子步骤来建立：

（A）首先确定对应时间下的睡眠指数                                                ，定义为1到100连续变化的离散值，作为睡眠指数计算模型的确立依据；

（B）定义

为步骤（2）中在第i个滑动窗中计算得到的特征参数的向量，

为在此滑动窗下医生根据临床确定的睡眠指数：

；

这里，

且

；等价为在限制条件

下，

为可调参数，最小化

；计算得到最优解，确定各个参数的权重系数，即为睡眠深度指数模型。

Images