一种基于振幅调制多尺度熵的睡眠动力学表征方法

1.本发明涉及认知神经科学与信息技术交叉技术领域，具体涉及一种基于振幅调制多尺度熵的睡眠动力学表征方法。


背景技术：

2.近年来，大脑节律的时间复杂性引起了人们对睡眠动力学的极大关注。脑电图是一种以不规则动态振荡为特征的复杂系统，具有非平稳性、间歇性和非线性。而熵的概念被广泛用于复杂性的评估，并且多年来，各种熵度量的方法已与多尺度方式相结合，以解决与系统复杂性相关的问题。越来越多的研究表明，多尺度熵是研究各种生理状态的有效方法，尤其是对心脏电生理活动的研究。但是，与心脏电生理活动不同，脑电信号由大脑尖峰神经元的同步电活动产生，在宽频带范围内波动，是一种快速的神经振荡，在多尺度方式下进行脑电信号的熵分析可能会稀释其复杂性并降低熵测量的性能。
3.睡眠动态过程大约包括五个睡眠阶段：睡眠的五个阶段包括清醒期（awake）、非快速眼动期（n1、n2和n3）和快速眼动期（rem）。直接对睡眠动态过程中的脑电信号进行波动分解只能提供部分估计，例如复杂度。近年来，把振幅调制作为分解的驱动波逐渐引起人们的广泛关注。睡眠阶段的自发转变以脑电信号不同频率的振幅调制为特征，其中高频节律反映了局部区域快速计算和区域计算的协调，相反的，低频节律起着交流的作用，沟通了外部感官输入和内部认知事件。因此，从脑电信号的振幅调制中获得的特征可以更好地描述睡眠动态，减少模糊。
4.传统的多尺度熵方法，直接对原始脑电信号本身进行多尺度熵分析，对于原始脑电信号的高频部分直接进行多尺度熵分析具有一定的模糊性，易受高频噪声的干扰，导致其在对睡眠状态下脑电信号的动态复杂性分析时鲁棒性较低。另外传统的多尺度熵方法，其产生的衍生指标包括多尺度熵曲线下的面积和曲线斜率，这两项衍生指标可以为区分不同的睡眠阶段提供参考，但是也只能进行粗略的睡眠阶段划分，无法实现清晰的睡眠阶段的划分。
5.考虑以上技术，目前尚未有有效的方法可以针对睡眠状态下脑电信号的动态复杂性进行有效的鲁棒估计，并实现清晰的睡眠阶段的划分。


技术实现要素：

6.有鉴于此，本发明提供了一种基于振幅调制多尺度熵的睡眠动力学表征方法，能够利用振幅调制多尺度熵解决快速振荡信号在多尺度方式下可能会削弱信号本身复杂度和熵测度性能的问题，可以更好地描绘睡眠动态，减少模糊，从而实现对睡眠状态下脑电信号的动态复杂性进行有效的鲁棒估计，并实现清晰的睡眠阶段的划分。
7.为达到上述目的，本发明的技术方案包括如下步骤：步骤一，通过睡眠监测仪采集获得睡眠数据，并进行睡眠阶段的划分，获得每个睡眠阶段对应的多通道脑电信号，并进行预处理。
8.步骤二，针对预处理后的多通道脑电信号进行调幅特征的提取，检测得到多通道脑电信号的波峰，对检测到的波峰进行三次样条插值，构建出多通道脑电信号的包络。
9.步骤三，对应每个频带设定不同的尺度因子，根据每一个已设定的尺度因子，提取多通道脑电信号的包络时间序列的多尺度熵。
10.以尺度因子为横轴，振幅调制多尺度熵为纵轴，绘制得到每个睡眠阶段对应的多通道脑电信号的振幅调制多尺度熵曲线，按照频带区分，计算每条曲线的振幅调制多尺度熵曲线下面积和平均曲线斜率。
11.步骤四，对同一频带不同睡眠阶段的振幅调制多尺度熵曲线的曲线下面积进行显著性分析，并对曲线斜率进行显著性分析；利用曲线下面积和曲线斜率的显著性分析结果作为睡眠动力学表征数据。
12.进一步地，步骤一中，通过睡眠监测仪采集获得睡眠数据，并进行睡眠阶段的划分，具体地：针对睡眠数据，以30秒为周期进行人工视觉标记和评分，完成睡眠分期，其中脑电信号采样频率为fs中。
13.进一步地，步骤一中，获得每个睡眠阶段对应的多通道脑电信号，并进行预处理，其中预处理的过程具体为：使用50 hz陷波滤波器去除工频干扰，使用低截止频率0.3 hz和高截止频率为35 hz的带通巴特沃斯滤波器提高信噪比；采用主成分分析方法和快速独立成分分析方法估计脑电信号的独立成分，以排除眼动，眨眼和心律对脑电图的干扰，与眼电信号和心电信号相关性高于0.5的成分被拒绝。
14.进一步地，步骤二中，针对预处理后的多通道脑电信号进行调幅特征的提取，检测得到多通道脑电信号的波峰，具体步骤如下：用顺序统计滤波器order-statistic filter消除所述预处理后的多通道脑电信号的伪峰，检测预处理后的脑电信号真正的峰值；所述检测预处理后的脑电信号真正的峰值，具体为：用滑动tukey窗对预处理后的多通道脑电信号进行加权，滑动窗口每次前移1个样本递归地执行加权步骤，取每个滑动窗口加权结果的最大值作为该窗的输出；所述多通道脑电信号积分功率与所有滑动窗口输出的最大值的交点即为脑电信号的波峰，即为预处理后的脑电信号真正的峰值。
15.进一步地，对应每个频带设定不同的尺度因子，具体为：脑电信号的频带划分为共五个频带，分别为：[2hz、4hz]区间为δ频带、[4hz、8hz]区间为θ频带、[8hz、14hz]区间为α频带、[14hz、28hz]区间为β频带和28hz以上为γ频带。
[0016]
对应每个频带设定不同的尺度因子，具体对应关系为：δ频带对应的尺度因子范围区间为[50,100）、θ频带带对应的尺度因子范围区间为[25,50）、α频带对应的尺度因子范围区间为[14,25）、β频带对应的尺度因子范围区间为[7,14）、γ频带对应的尺度因子范围区间为[1,7）。
[0017]
进一步地，根据每一个已设定的尺度因子，提取多通道脑电信号的包络时间序列的多尺度熵，具体步骤如下：s301）针对已设定的尺度因子τ，对脑电信号包络时间序列进行粗粒度处理，输入被分割为窗口大小为τ的非重叠窗口，n为脑电信号包络时间序列长度，对
每个窗口内的数据进行平均，第j个窗口的输出为，尺度τ取1时就是原始脑电信号包络时间序列；由此获得即为脑电信号包络时间序列的粗粒度序列。
[0018]
s302）递归地计算s301）得到的粗粒度序列的样本熵，构建一个跨时间尺度的样本熵谱；具体地：对于输入长度为l的粗粒度信号，构造l-m个模板向量xm(1),xm(2),
ꢀ…ꢀ
,xm(l-m)，，其中m代表嵌入维数，计算每两个模板向量的距离d[xm(p),xm(q)]。
[0019]
如果两个模板向量的距离小于预定误差r，那么这两个模板向量相似。
[0020]
距离小于预定误差r的模板向量对的数量b
p
，b
p
可能的值有l-m-1个，p满足1≤p≤l-m；则任意两个模板向量在预定误差范围r内相似的频率b
pm
(r)为频率密度为：。
[0021]
在每一个模板向量中增加一个额外的样本，该额外的样本自行设定，即模板向量的维度变成m+1，，同样地，x
m+1
(p)和x
m+1
(q)在预定误差范围r内相似的频率为；其中是增加额外向本后的模板向量中，距离小于预定误差r的模板向量对的数量。
[0022]
频率密度am(r)为：。
[0023]
模板向量的维度为m时，在预定误差r内相似的模板向量对总数b(r)为：。
[0024]
模板向量的维度为m+1时，在预定误差r内相似的模板向量对总数为：。
[0025]
最后，脑电信号包络时间序列的多尺度熵为：。
[0026]
进一步地，d[xm(p),xm(q)]定义为两个向量各元素间的最大差值，即x
p+k
，x
q+k
分别为xm(p),xm(q)中的元素。
[0027]
进一步地，嵌入维数取值为2或3，预设误差r设置为0.25
×
sd，其中sd是脑电信号
包络时间序列的标准差。
[0028]
有益效果：本发明提供的一种基于振幅调制多尺度熵的睡眠动力学表征方法，计算时变振幅调制脑电图的复杂性特征，反映睡眠阶段的转换。使用复杂的头皮脑电信号包络以多尺度方式表达睡眠动态等生理变化，规避了原始脑电信号的宽带特性，可以反映更大尺度的神经元群体。相比于原始脑电信号，本发明提供的一种基于振幅调制多尺度熵的睡眠动力学表征方法，其脑电信号振幅调制的频率更低，消除了高频伪影的干扰，所以与传统的多尺度熵方法相比，振幅调制多尺度熵为睡眠状态下脑电信号的动态复杂性提供了有效的鲁棒估计，脑电信号振幅调制多尺度熵曲线的衍生特征能够在某些睡眠阶段呈现更少的模糊性，并更清晰地分开不同的睡眠阶段。振幅调制多尺度熵解决了快速振荡信号在多尺度方式下可能会削弱信号本身复杂度和熵测度性能的问题，可以更好地描绘睡眠动态，减少模糊。因此本发明适用于提取用于睡眠分期的特征，高效可靠，易于软件化。
附图说明
[0029]
图1为本发明提供的一种基于振幅调制多尺度熵的睡眠动力学表征方法流程图；图2.不同睡眠阶段的多尺度熵曲线图；图3不同睡眠阶段的振幅调制多尺度熵曲线图。
具体实施方式
[0030]
下面结合附图并举实施例，对本发明进行详细描述。
[0031]
为突破在多尺度方式下分析脑电信号这种快速神经振荡会造成复杂性的稀释和熵测量性能降低等技术难点，本发明提出基于振幅调制的多尺度熵分析方法，并用于睡眠分期研究，详细过程如图1，包括如下步骤：步骤一，通过睡眠监测仪采集获得睡眠数据，并进行睡眠阶段的划分，获得每个睡眠阶段对应的多通道脑电信号，并进行预处理。
[0032]
采集过程中，受试者按照实验要求完成实验操作，应用多导睡眠监测仪采集受测者整晚睡眠数据，根据美国医学会的指导方针以30秒为周期进行人工视觉标记和评分，完成睡眠阶段的划分，其中脑电信号采样频率为fs，存储在计算机中。
[0033]
对采集到的多通道脑电信号进行如下预处理：使用50 hz陷波滤波器去除工频干扰，使用低截止频率0.3 hz和高截止频率为35 hz的带通巴特沃斯滤波器提高信噪比。采用主成分分析方法和快速独立成分分析方法估计脑电信号的独立成分，以排除眼动，眨眼和心律对脑电图的污染，与眼电信号和心电信号相关性高于0.5的成分被拒绝。
[0034]
步骤二，针对预处理后的多通道脑电信号进行调幅特征的提取，检测得到多通道脑电信号的波峰，对检测到的波峰进行三次样条插值，构建出多通道脑电信号的包络；其中对预处理后的多维脑电信号进行调幅特征的提取，检测得到多通道脑电信号的波峰具体步骤如下：用顺序统计滤波器（order-statistic filter）消除脑电信号的伪峰，检测预处理后的脑电信号真正的峰值。准确的说，用滑动tukey窗对预处理后脑电信号进行加权，滑动窗口每次前移1个样本递归地执行加权步骤，取每个滑动窗口加权结果的最大值作为该窗的输出。脑电信号积分功率与所有滑动窗口输出的最大值的交点即为脑电信
号的波峰，即为预处理后的脑电信号真正的峰值。
[0035]
对检测到的波峰进行三次样条插值，构建出自适应脑电信号的包络。
[0036]
步骤三，对应每个频带设定不同的尺度因子，根据每一个已设定的尺度因子，提取多通道脑电信号的包络时间序列的多尺度熵；其中根据每一个已设定的尺度因子，提取多通道脑电信号的包络时间序列的多尺度熵，具体步骤如下：s301）用尺度因子τ对脑电信号包络时间序列进行粗粒度处理，输入被分割为窗口大小为τ的非重叠窗口（n为脑电信号包络时间序列长度），对每个窗口内的数据进行平均，例如第j个窗口的输出满足，这一过程类似于移动平均并向下采样。值得注意的是，尺度τ取1时就是原始脑电信号包络时间序列。即为脑电信号包络时间序列的粗粒度序列。
[0037]
s302）递归地计算上一步得到的粗粒度序列的样本熵，构建一个跨时间尺度的样本熵谱。具体来说，对于输入长度为l的粗粒度信号，构造l-m个模板向量xm(1),xm(2),
ꢀ…ꢀ
,xm(l-m)，，其中m代表嵌入维数，本发明实施例中，嵌入维数可取值为2或3。计算每两个模板向量的距离d[xm(p),xm(q)]。
[0038]
其中d[xm(p),xm(q)]定义为两个向量各元素间的最大差值，即其中x
p+k
，x
q+k
分别为xm(p),xm(q)中的元素。
[0039]
如果两个模板向量的距离小于预定误差r，本发明实施例中，预设误差r设置为0.25
×
sd，其中sd是脑电信号包络时间序列的标准差，那么这两个模板向量相似。设b
p
是距离小于预定误差r的模板向量对的数量，所以b
p
可能的值有l-m-1个，p满足1≤p≤l-m。所以任意两个模板向量在预定误差范围r内相似的频率b
pm
(r)为：。
[0040]
频率密度为：。
[0041]
在每一个模板向量中增加一个额外的样本，该额外的样本可自行设定，即模板向量的维度变成m+1，，同样地，x
m+1
(p)和x
m+1
(q)在预定误差范围r内相似的频率a
pm
(r)为。
[0042]
是增加额外向本后的模板向量中，距离小于预定误差r的模板向量对的数量，频
率密度am(r)为。
[0043]
模板向量的维度为m时，在预定误差r内相似的模板向量对总数b(r)可以计算为：。
[0044]
模板向量的维度为m+1时，在预定误差r内相似的模板向量对总数可以计算为：。
[0045]
最后，脑电信号包络时间序列的多尺度熵定义为：。
[0046]
本步骤是针对脑电信号包络时间序列求解多尺度熵，高速振荡频率对应于较短的尺度周期，这表明传统的多尺度熵结果容易受到高频随机噪声的污染，从而降低精度。与此同时，脑电信号的采样率通常在100或200 hz或更高，其中最高的感兴趣波段频率在50 hz左右，例如γ波段。因此，对于脑电信号等具有多源相互作用的高速振荡，每个感兴趣的波段的尺度范围是相当不足的。而脑电信号的包络是低频振荡波，对脑电信号包络进行多尺度熵分析可以反映更大尺度的神经元群体。并且它规避了原始脑电信号的宽带特性，从而为睡眠过程中脑电信号的动态复杂性提供了鲁棒估计。
[0047]
脑电信号的频带划分为共五个频带，分别为：[2hz、4hz]区间为δ频带、[4hz、8hz]区间为θ频带、[8hz、14hz]区间为α频带、[14hz、28hz]区间为β频带和28hz以上为γ频带；对应每个频带设定不同的尺度因子，具体对应关系为：δ频带对应的尺度因子范围区间为[50,100）、θ频带带对应的尺度因子范围区间为[25,50）、α频带对应的尺度因子范围区间为[14,25）、β频带对应的尺度因子范围区间为[7,14）、γ频带对应的尺度因子范围区间为[1,7）。
[0048]
进行上述振幅调制多尺度熵计算之后，以尺度因子为横轴，振幅调制多尺度熵为纵轴，绘制五个睡眠阶段脑电信号的振幅调制多尺度熵曲线，分频带计算五条曲线的振幅调制多尺度熵曲线下面积和平均曲线斜率。
[0049]
步骤四，对同一个脑电频带不同睡眠阶段的振幅调制多尺度熵曲线的曲线下面积进行显著性分析，并对曲线斜率进行显著性分析，其他脑电频带重复此操作。
[0050]
利用曲线下面积和曲线斜率的显著性分析结果作为睡眠动力学表征数据。
[0051]
本发明已在isruc-sleep数据集上进行验证，该数据集中脑电信号采样率为200hz，选取的通道为c4-a1。通过绘制五个睡眠阶段脑电信号的振幅调制多尺度熵曲线，发现与传统的多尺度熵曲线相比，振幅调制多尺度熵曲线轨迹相对稳定，由图2和图3所示。此外，某些脑电频带不同睡眠阶段的振幅调制多尺度熵曲线下面积和平均曲线斜率两个特征均具有显著统计学差异，层次清晰，由表1所示。
[0052]
表1. 使用振幅调制多尺度熵衍生特征预测睡眠阶段的线性混合模型
由此可见，基于振幅调制多尺度熵来表征睡眠动态过程中大脑复杂性和状态的方法可作为睡眠分期的补充工具，有助于提高睡眠分期准确度，在睡眠健康监测中具有一定的潜在价值和应用前景。
[0053]
综上所述，以上仅为本发明的较佳实施例而已，并非用于限定本发明的保护范围。凡在本发明的精神和原则之内，所作的任何修改、等同替换、改进等，均应包含在本发明的保护范围之内。

技术特征：
1.一种基于振幅调制多尺度熵的睡眠动力学表征方法，其特征在于，包括：步骤一，通过睡眠监测仪采集获得睡眠数据，并进行睡眠阶段的初步划分，获得每个睡眠阶段对应的多通道脑电信号，并进行预处理；步骤二，针对预处理后的多通道脑电信号进行调幅特征的提取，检测得到多通道脑电信号的波峰，对检测到的波峰进行三次样条插值，构建出多通道脑电信号的包络时间序列；步骤三，对应每个频带设定不同的尺度因子，根据每一个已设定的尺度因子，提取多通道脑电信号的包络时间序列的多尺度熵；以尺度因子为横轴，振幅调制多尺度熵为纵轴，绘制得到每个睡眠阶段对应的多通道脑电信号的振幅调制多尺度熵曲线，按照频带区分，计算每条曲线的振幅调制多尺度熵曲线下面积和平均曲线斜率；步骤四，对同一频带不同睡眠阶段的振幅调制多尺度熵曲线的曲线下面积进行显著性分析，并对曲线斜率进行显著性分析；利用曲线下面积和曲线斜率的显著性分析结果作为睡眠动力学表征数据。2.如权利要求1所述的一种基于振幅调制多尺度熵的睡眠动力学表征方法，其特征在于，所述步骤一中，通过睡眠监测仪采集获得睡眠数据，并进行睡眠阶段的划分，具体地：针对睡眠数据，以30秒为周期进行人工视觉标记和评分，完成睡眠分期，其中脑电信号采样频率为fs。3.如权利要求1所述的一种基于振幅调制多尺度熵的睡眠动力学表征方法，其特征在于，所述步骤一中，获得每个睡眠阶段对应的多通道脑电信号，并进行预处理，其中预处理的过程具体为：使用50 hz陷波滤波器去除工频干扰，使用低截止频率0.3 hz和高截止频率为35 hz的带通巴特沃斯滤波器提高信噪比；采用主成分分析方法和快速独立成分分析方法估计脑电信号的独立成分，以排除眼动，眨眼和心律对脑电图的干扰，与眼电信号和心电信号相关性高于0.5的成分被拒绝。4.如权利要求1所述的一种基于振幅调制多尺度熵的睡眠动力学表征方法，其特征在于，所述步骤二中，针对预处理后的多通道脑电信号进行调幅特征的提取，检测得到多通道脑电信号的波峰，具体步骤如下：用顺序统计滤波器order-statistic filter消除所述预处理后的多通道脑电信号的伪峰，检测预处理后的脑电信号真正的峰值；所述检测预处理后的脑电信号真正的峰值，具体为：用滑动tukey窗对预处理后的多通道脑电信号进行加权，滑动窗口每次前移1个样本递归地执行加权步骤，取每个滑动窗口加权结果的最大值作为该窗的输出；所述多通道脑电信号积分功率与所有滑动窗口输出的最大值的交点即为脑电信号的波峰，即为预处理后的脑电信号真正的峰值。5.如权利要求1所述的一种基于振幅调制多尺度熵的睡眠动力学表征方法，其特征在于，所述对应每个频带设定不同的尺度因子，具体为：脑电信号的频带划分为共五个频带，分别为：[2hz、4hz]区间为δ频带、[4hz、8hz]区间为θ频带、[8hz、14hz]区间为α频带、[14hz、28hz]区间为β频带和28hz以上为γ频带；对应每个频带设定不同的尺度因子，具体对应关系为：δ频带对应的尺度因子范围区间为[50,100）、θ频带对应的尺度因子范围区间为[25,
50）、α频带对应的尺度因子范围区间为[14,25）、β频带对应的尺度因子范围区间为[7,14）、γ频带对应的尺度因子范围区间为[1,7）。6.如权利要求1所述的一种基于振幅调制多尺度熵的睡眠动力学表征方法，其特征在于，所述根据每一个已设定的尺度因子，提取多通道脑电信号的包络时间序列的多尺度熵，具体步骤如下：s301）针对已设定的尺度因子τ，对脑电信号包络时间序列进行粗粒度处理，输入被分割为窗口大小为τ的非重叠窗口，n为脑电信号包络时间序列长度，对每个窗口内的数据进行平均，第j个窗口的输出为，尺度τ取1时就是原始脑电信号包络时间序列；由此获得即为脑电信号包络时间序列的粗粒度序列；s302）递归地计算s301）得到的粗粒度序列的样本熵，构建一个跨时间尺度的样本熵谱；具体地：对于输入长度为l的粗粒度信号，构造l-m个模板向量x
m
(1), x
m
(2),
ꢀ…ꢀ
, x
m
(l-m)，，其中m代表嵌入维数，计算每两个模板向量的距离d[x
m
(p), x
m
(q)]；如果两个模板向量的距离小于预定误差r，那么这两个模板向量相似；距离小于预定误差r的模板向量对的数量b
p
，b
p
可能的值有l-m-1个，p满足1≤p≤l-m；则任意两个模板向量在预定误差范围r内相似的频率b
pm
(r)为；频率密度为：；在每一个模板向量中增加一个额外的样本，该额外的样本自行设定，即模板向量的维度变成m+1，，同样地，x
m+1
(p)和x
m+1
(q)在预定误差范围r内相似的频率为；其中是增加额外向本后的模板向量中，距离小于预定误差r的模板向量对的数量；频率密度a
m
(r)为：；模板向量的维度为m时，在预定误差r内相似的模板向量对总数b(r)为：；模板向量的维度为m+1时，在预定误差r内相似的模板向量对总数为：；
最后，脑电信号包络时间序列的多尺度熵为：。7.如权利要求6所述的一种基于振幅调制多尺度熵的睡眠动力学表征方法，其特征在于，所述d[x
m
(p), x
m
(q)]定义为两个向量各元素间的最大差值，即x
p+k
，x
q+k
分别为x
m
(p), x
m
(q)中的元素。8.如权利要求6或7所述的一种基于振幅调制多尺度熵的睡眠动力学表征方法，其特征在于，所述嵌入维数取值为2或3，预设误差r设置为0.25
×
sd，其中sd是脑电信号包络时间序列的标准差。

技术总结
本发明公开了一种基于振幅调制多尺度熵的睡眠动力学表征方法，属于认知神经科学和信息技术交叉领域。首先对采集到的已进行睡眠分期的多通道睡眠脑电进行预处理；检测脑电信号的波峰，并包络；将脑电信号包络作为时间序列，计算这一时间序列的多尺度熵，即为脑电信号的振幅调制多尺度熵；绘制脑电信号所有睡眠阶段的振幅调制多尺度熵曲线；对不同睡眠阶段的振幅调制多尺度熵曲线的曲线下面积进行显著性分析，并对不同睡眠阶段的振幅调制多尺度熵曲线的曲线斜率进行显著性分析；利用曲线下面积和曲线斜率的显著性分析结果作为睡眠动力学表征数据，实现对睡眠状态下脑电信号的动态复杂性进行有效的鲁棒估计，并实现清晰的睡眠阶段的划分。段的划分。段的划分。


技术研发人员：史文彬 冯欢 叶建宏
受保护的技术使用者：北京理工大学
技术研发日：2023.06.27
技术公布日：2023/8/1