CN108416367A - 基于多传感器数据决策级融合的睡眠分期方法 - Google Patents

基于多传感器数据决策级融合的睡眠分期方法 Download PDF

Info

Publication number
CN108416367A
CN108416367A CN201810125718.1A CN201810125718A CN108416367A CN 108416367 A CN108416367 A CN 108416367A CN 201810125718 A CN201810125718 A CN 201810125718A CN 108416367 A CN108416367 A CN 108416367A
Authority
CN
China
Prior art keywords
radar
signal
model
audio
per minute
Prior art date
Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
Granted
Application number
CN201810125718.1A
Other languages
English (en)
Other versions
CN108416367B (zh
Inventor
顾陈
周燕萍
洪弘
蒋洁
李彧晟
孙理
朱晓华
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
Nanjing University of Science and Technology
Original Assignee
Nanjing University of Science and Technology
Priority date (The priority date is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the date listed.)
Filing date
Publication date
Application filed by Nanjing University of Science and Technology filed Critical Nanjing University of Science and Technology
Priority to CN201810125718.1A priority Critical patent/CN108416367B/zh
Publication of CN108416367A publication Critical patent/CN108416367A/zh
Application granted granted Critical
Publication of CN108416367B publication Critical patent/CN108416367B/zh
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Classifications

    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F18/00—Pattern recognition
    • G06F18/20—Analysing
    • G06F18/24—Classification techniques
    • G06F18/241—Classification techniques relating to the classification model, e.g. parametric or non-parametric approaches
    • G06F18/2413—Classification techniques relating to the classification model, e.g. parametric or non-parametric approaches based on distances to training or reference patterns
    • G06F18/24147—Distances to closest patterns, e.g. nearest neighbour classification
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F18/00—Pattern recognition
    • G06F18/20—Analysing
    • G06F18/24—Classification techniques
    • G06F18/241—Classification techniques relating to the classification model, e.g. parametric or non-parametric approaches
    • G06F18/2415—Classification techniques relating to the classification model, e.g. parametric or non-parametric approaches based on parametric or probabilistic models, e.g. based on likelihood ratio or false acceptance rate versus a false rejection rate
    • G06F18/24155—Bayesian classification
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F18/00—Pattern recognition
    • G06F18/20—Analysing
    • G06F18/25—Fusion techniques
    • G06F18/253—Fusion techniques of extracted features

Landscapes

  • Engineering & Computer Science (AREA)
  • Data Mining & Analysis (AREA)
  • Theoretical Computer Science (AREA)
  • Physics & Mathematics (AREA)
  • Bioinformatics & Cheminformatics (AREA)
  • Bioinformatics & Computational Biology (AREA)
  • Computer Vision & Pattern Recognition (AREA)
  • Evolutionary Biology (AREA)
  • Evolutionary Computation (AREA)
  • Artificial Intelligence (AREA)
  • General Engineering & Computer Science (AREA)
  • General Physics & Mathematics (AREA)
  • Life Sciences & Earth Sciences (AREA)
  • Probability & Statistics with Applications (AREA)
  • Radar Systems Or Details Thereof (AREA)
  • Measurement Of The Respiration, Hearing Ability, Form, And Blood Characteristics Of Living Organisms (AREA)

Abstract

本发明公开了一种基于多传感器数据决策级融合的睡眠分期方法，该方法首先使用雷达传感器和音频传感器采集整夜的雷达和音频数据，提取雷达和音频信号特征；根据特征数将其进行数据分类，根据数据分类建立雷达剩余片段模型和雷达+音频片段模型；接着用分类器对雷达剩余片段模型里的雷达特征进行识别分类，得到模型预测结果1，用分类器对雷达+音频片段模型里的雷达和音频特征进行识别分类，得到模型预测结果A和B；然后用朴素贝叶斯模型对模型预测结果A和B进行决策，得到模型预测结果2；最后将模型预测结果1和模型预测结果2进行时序拼接，得到整夜睡眠分期结果。该方法简单易行，准确率高，与实际情况相吻合。

Description

基于多传感器数据决策级融合的睡眠分期方法
技术领域
本发明属于生命体征监测技术领域，是一种非接触式睡眠分期方法。
背景技术
在睡眠过程中，人的脑电图、心电图等多种生理信号发生变化，而这些信号变化随着入睡的深浅度有一定的相关性变化。医学上根据脑电波、心电、眼电等生理表现人为的将睡眠划分为几个阶段，实际上各个睡眠阶段之间并没有很明确的界限，而是周期之间逐渐转变且有承接性的，整夜的睡眠阶段也是呈现周期性的，有一定的规律。
目前睡眠监测产品多采用接触式监测，需要在受试者身体的多个部位贴数个电极，佩戴口鼻气流管和腹部绑带。接触式的监测在很大程度上影响了受试者的生理和心理情况，且测试环境复杂，需要专业医护人员操作，价格相对昂贵。多传感器数据融合的方法进行睡眠分期是一种非接触式分期方法，相比单一传感器数据具有以下优点：多样性、实时性、鲁棒性和高效性。为了最大程度地利用多个传感器的数据获得更为全面多元化的信息，需要对数据进行融合理。对多个用于监测睡眠的传感器进行数据融合，可以提高整个睡眠分期系统的可靠性、实时性和利用率，能获得更为全面的睡眠分期结果，这将使全面睡眠分期研究更具意义。但是现有技术中尚无相关描述。
发明内容
本发明的目的在于提出供一种基于多传感器决策层融合的睡眠分期方法。
实现本发明目的的技术解决方案为：一种基于多传感器数据决策级融合的睡眠分期方法，包括以下步骤：
步骤1、测试者平躺在测试床上，将雷达传感器设置于人体的正上方，用于采集雷达回波信号；将音频传感器设置于朝向人脸的方向，用于采集呼吸和鼾声信号；
步骤2、开启雷达传感器和音频传感器采集信号，获得雷达回波信号X(t)及音频信号Y(t)，之后对雷达回波信号X(t)进行滤波，获得呼吸信号和心跳信号，对雷达回波信号X(t)用幅值门限检测法作为预判的条件，检测出体动信号，对音频信号Y(t)进行降噪和事件检测获得呼吸信号Y1(t)和鼾声信号Y2(t)；
步骤3、对步骤2获得的雷达回波信号X(t)进行特征参数提取，根据呼吸信号得到每分钟呼吸次数N_bre、每分钟呼吸次数的方差S_bre2、每分钟呼吸信号幅度的差值积累Ada_bre、呼吸信号中的体动信号参数move_b，根据心跳信号得到每分钟心跳次数N_hrt、每分钟心跳次数的方差S_hrt2、每分钟心跳信号幅度的差值积累Ada_hrt、心跳信号中的体动信号参数move_h，睡眠分期相关特征REM期参数REM、DEEP期参数DEEP、样本熵SampEn、时序特征at；
所述提取11个雷达回波信号特征以及1个时序特征对睡眠阶段进行划分，12个特征具体如下：
(1)每分钟的呼吸次数N_bre：由相邻峰值点即波峰或波谷的差值计算每次呼吸时长，从而计算每分钟的呼吸次数N_bre；
(2)每分钟呼吸次数的方差S_bre2：
其中，N_bre(g)为第g分钟呼吸的次数，为整夜信号中每分钟呼吸次数的平均值，l_b为计算每分钟呼吸次数的方差的窗长；
(3)每分钟呼吸信号幅度的差值积累Ada_bre：
其中，A_bre(g)为第g分钟呼吸信号幅度，A_bre(g+1)为第g+1分钟呼吸信号幅度，N_i为每分钟内插值点数，取一个常数；
(4)呼吸信号中的体动信号参数：设置呼吸信号中体动信号的平均幅值move_b_thres为门限，判定超过门限的点数来计算体动参数，记为move_b；
(5)每分钟心跳次数N_hrt：计算相邻波谷值之间的时间差来计算每分钟心跳次数N_hrt；
(6)每分钟心跳次数的方差S_hrt2：
其中，N_hrt(g1)为第g1分钟心跳信号的幅值，为整夜心跳信号平均每分钟幅值，N_i为每分钟内插值点数；
(7)每分钟心跳信号幅值的差值积累Ada_hrt(g1)：
其中，A_hrt(g1)为第g1分钟心跳信号幅度，A_hrt(g1+1)为第g1+1分钟心跳信号幅度，N_i为每分钟内插值点数；
(8)心跳信号中的体动信号参数：设置心跳信号中体动信号的平均幅值move_h_thres为门限，判定超过门限的点数来计算体动参数，记为move_h；
(9)取某分钟内前后T秒的呼吸次数的差值之间计算关系的参数REM：
其中，表示当前第g3+g4分钟内前T秒内的呼吸次数，表示当前第g3+g4分钟内后T秒的呼吸次数，T的取值范围是25～30，λ为一个常数；
(10)用来区分深度睡眠阶段的重要特征参数DEEP：
其中，Abody表示睡眠中的体动信号的幅值，Arespir表示睡眠中的呼吸信号的幅值，DEEP表示体动信号在总的呼吸加体动信号总和中所占比例；
(11)样本熵SampEn：由雷达回波信号X(t)采样获得雷达回波信号的一维时间序列X(n)，即雷达回波信号时间序列，1≤n≤k_r，重构为相空间矢量X*(ww)，1≤ww≤k_r-(m*-1)，m*为嵌入维数，k_r为序列长度，定义相空间中任意两个矢量X*(w1)和X*(w2)的距离为：
其中，kk_r为一个常数，其取值范围是1≤kk_r≤m*-1，
求取样本熵的具体步骤为：
给定相似容限r，r＝hh*SD，hh是取值范围为0.1～0.25的常数，SD是
X(n)的标准差，对空间中每一个X*(w3)，计算的数目与总数k_r-m*+1的比值，记为
表示的数目；
对w3求平均值，记为
维数m*加1，重复上述步骤，得到
雷达回波信号时间序列X(n)的样本熵为：
实际中k_r为有限值，则样本熵估计为
(12)时序特征at：若{a1,a2,......,aN_a}为雷达回波信号特征组按照相对时间顺序排列的N_a个序号，对应于所处睡眠时间段，则at＝a/amax为睡眠分期的时序特征。
步骤4、对步骤2获得的音频信号Y(t)进行特征参数提取，包括睡眠呼吸相关参数以及鼾声相关参数，呼吸相关参数为：每分钟呼吸的次数N_bre'、呼吸次数差值BVP、呼吸周期BC、均方根误差RMSE、每分钟内呼吸次数的方差S_bre'2、能量差值累积EDA_bre、过零率即N_zcr、共振峰1即Formant1、共振峰2即Formant2、共振峰3即Formant3、共振峰1的方差即Formant_var1、共振峰2的方差即Formant_var2、共振峰3的方差即Formant_var3；鼾声相关参数为：李雅普诺夫指数LLE、时间延迟嵌入维数m、近似熵ApEn、统计嵌入维数m≥4的个数N_m，鼾声事件时间序列的维数D、去趋势波动分析的阿尔法系数一α1、去趋势波动分析的阿尔法系数二α2、样本熵SampEn*、香农熵H、时序特征bt；
所述对音频信号Y(t)使用能熵比双门限法进行语音端点检测，检测出有效的呼吸事件Y1(t)和有效的鼾声事件Y2(t)，由Y1(t)和Y2(t)进行采样得到呼吸事件时间序列Y1(n)和鼾声事件时间序列Y2(n)，提取出23个音频信号特征以及1个时序特征对睡眠阶段进行划分，24个特征参数分别为：
呼吸相关参数：
(1)每分钟呼吸次数N_bre'：取每分钟的呼吸次数；
(2)每分钟里前后内呼吸次数差值BVP：
BVP＝|N_bre'(s)-N_bre'(s-1)|
其中，N_bre'(s)为第s分钟内呼吸次数，N_bre'(s-1)为第s-1分钟内呼吸次数；
(3)呼吸周期BC，取两次呼吸事件之间的时间间隔；
(4)均方根误差RMSE：
其中，Y1(s1)为第s1个点的呼吸信号的幅值，为呼吸信号幅值的平均值，Ns为采样点数；
(5)每分钟呼吸次数的方差S_bre'2：
其中，N_bre'(s2)为第s2分钟的呼吸次数，为整夜呼吸信号每分钟呼吸次数的平均值，l_b'为计算呼吸次数的方差的窗长，l_b'取一个常数；
(6)能量差值累积EDA_bre：
其中，是第u个呼吸事件的能量，是第u+1个呼吸事件的能量，N_bre'是每分钟呼吸事件总数，Y1u(q)是第u个呼吸事件的第q个采样点对应的呼吸信号幅度，Y1(u+1)(q1)是第u+1个呼吸事件的第q1个采样点对应的呼吸信号幅度，N_p为第u个呼吸事件的总点数，N_p'为第u+1个呼吸事件的总点数；
(7)过零率：呼吸信号的过零点的次数N_zcr；
(8)第一共振峰Formant1；
(9)第二共振峰Formant2；
(10)第三共振峰Formant3；
(11)第一共振峰的方差Formant_var1；
(12)第二共振峰的方差Formant_var2；
(13)第三共振峰的方差Formant_var3；
鼾声相关参数：
(14)时间延迟给定一维的鼾声事件时间序列Y2(n)后，计算序列的自相关函数并作该函数关于时间的函数的图形，当函数值下降到初始值的时，对应的时间就是时间延迟；e取一个极小值；
(15)嵌入维数m：利用时间延迟将一维鼾声事件时间序列嵌入到m维空间中；
(16)LLE：根据时间延迟和嵌入维数m重构相空间W，对每个点Wj寻找其最近邻点Wj'，计算Wj到Wj'的距离dj(0)＝|Wj-Wj'|，其中，j为相空间中点的坐标；
对每个点Wj，计算其与最近邻点Wj'在第i步前向演化后的距离：
dj(i)＝|Wj+i-W′j+i|＝dj(0)×eLLE×i；
由下式算得最大李雅普诺夫指数：
其中，p_n为非零的dj(i)的个数；
(17)ApEn：将鼾声事件时间序列Y2(n)，1≤n≤k，重构为m维相空间矢量Y2'(i1)，1≤i1≤k-(m-1)，k为序列长度，定义相空间中任意两个矢量Y2'(j1)和Y2'(j2)的距离为：
d'j1j2＝max|Y2'(j1+kk)-Y2'(j2+kk)|，1≤j1,j2≤k-m+1,j1≠j2，其中，kk是一个常数，其取值范围是1≤kk≤m-1，
求取近似熵的具体步骤为：
给定相似容限r1，r1＝hh_a*SD_a，hh_a是取值范围为0.1～0.25的常数，SD_a是鼾声事件时间序列Y2(n)的标准差，对空间中每一个Y2'(i1)，计算d'j1j2＜r1的数目与总数k-m+1的比值，记为
N_dj1j2表示d'j1j2＜r1的数目；
将取对数，对i1求平均值，记为φm(r1)：
维数m加1，重复上述步骤，得到φm+1(r1)；
鼾声事件时间序列Y2(n)的近似熵为：
实际中，k为有限值，则近似熵估计为
ApEn(m,r1)＝φm(r1)-φm+1(r1)；
(18)统计嵌入维数m≥4的个数N_m；
(19)鼾声事件时间序列的维数D；
(20-21)去趋势波动分析的阿尔法系数一α1和去趋势波动分析的阿尔法系数二α2：对长度为k的鼾声事件时间序列Y2(n)，1≤n≤k，计算其累积差值：
其中，是时间序列的均值：
将y2(n)划分为不重叠的n2个长度为l的区间，l即为时间尺度(窗长)，n2为区间(或窗口)数量；
对每一段鼾声事件时间序列采用最小二乘法拟合出局部趋势y'2(n)；
剔除y2(n)中每个区间的局部趋势，并计算新序列的均方根：
改变窗口长度l的大小，重复上述步骤，画出纵坐标为log[F(n)]，横坐标为log(n)的曲线，则曲线的斜率为时间序列的标度指数α；
取窗口长度l1时对应的标度指数α为去趋势波动分析的阿尔法系数一α1，
取窗口长度l2时对应的标度指数α为去趋势波动分析的阿尔法系数二α2；
(22)样本熵SampEn*：将鼾声事件时间序列Y2(n)，1≤n≤k重构为m维相空间矢量Y2 *(i2)，1≤i2≤k-(m-1)，定义相空间中任意两个矢量Y2 *(j3)和Y2 *(j4)的距离为：
其中，kk_1是一个常数，其取值范围是1≤kk_1≤m-1，
求取样本熵的具体步骤为：
给定相似容限r1，对空间中每一个计算d* j3j4＜r1的数目与总数k-m+1的比值，记为
N_dj3j4表示d* j3j4＜r1的数目；
对i2求平均值，记为
维数m加1，重复上述步骤，得到Bm+1(r1)；
鼾声事件时间序列的样本熵为：
实际中k为有限值，则样本熵估计为：
(23)香农熵H：即信息熵，定义为离散随机事件的出现概率，令每个鼾声事件时间序列Y2(n)，对应的概率为p(Y2＝yn)，1≤n≤k，则随机变量的信息熵为
(24)时序特征bt：若{b1,b2,......,bN_b}为音频信号特征组按照相对时间顺序排列的N_b个序号，对应于所处睡眠时间段，则bt＝b/bmax为睡眠分期的时序特征。
步骤5、根据步骤3获得的11个雷达回波信号特征以及时序特征at、步骤4获得的23个音频信号特征以及时序特征bt，构建决策级融合系统模型，包含雷达剩余片段模型和雷达+音频片段模型，雷达剩余片段模型是单独由雷达传感器的特征数据构建的，雷达+音频片段模型是由雷达传感器和音频传感器共同获取的特征数据构建的；
构建雷达+音频片段模型和雷达剩余片段模型的基本步骤如下：
步骤5-1、构建雷达+音频片段模型：从整夜的音频数据区间[t1,t2nn]里筛选音频传感器获取的数据，截取有用的音频片段1、音频片段2、……、音频片段NN，每个片段对应的时间区间分别为[t1,t2],[t3,t4],...,[t2nn-1,t2nn]，从整夜雷达数据[t1,t2nn]区间里截取的[t1,t2],[t3,t4],...,[t2nn-1,t2nn]区间里的数据，分别对应雷达数据片段1、雷达片段2、……、雷达片段NN，音频片段包含23个音频特征和1个时序特征共24个特征，雷达片段包含11个雷达回波信号特征和1个时序特征共12个数据特征，把雷达片段特征和音频片段特征拼接成新的融合特征片段：雷达+音频片段；
步骤5-2、构建雷达剩余片段模型：将原始雷达数据中不存在音频数据时间段的雷达数据进行分割，即[t2,t3],[t4,t5],...,[t2nn-2,t2nn-1]区间对应的雷达数据片段，分成雷达剩余片段1、雷达剩余片段2、……、雷达剩余片段NN-1，包含11个雷达回波信号特征、1个时序特征共12个数据特征；
雷达+音频片段模型和雷达剩余片段模型表如下：
由雷达+音频片段1、雷达+音频片段2、……、雷达+音频片段NN构成的模型为雷达+音频片段模型，由雷达剩余片段1、雷达剩余片段2、……、雷达剩余片段NN-1构成的模型即为雷达剩余片段模型。
步骤6、在雷达剩余片段模型里，对雷达数据特征，进行识别分类得到模型预测结果1；在雷达+音频片段模型里，对雷达数据特征进行识别分类，得到模型预测结果A；对音频数据特征进行识别分类，得到模型预测结果B；
所述对雷达剩余片段模型和雷达+音频片段模型里的数据特征进行识别分类，得到模型预测结果，具体为：雷达剩余片段模型里的12个数据特征，用Subspace KNN分类器进行识别分类，得到模型预测结果1，对雷达+音频片段模型里的雷达片段数据的12个数据特征，用Subspace KNN分类器进行识别分类，得到模型预测结果A，对雷达+音频片段模型里的音频片段数据的24个数据特征，用Bagged Trees分类器进行识别分类，得到模型预测结果B。
步骤7、将步骤6获得的模型预测结果A和模型预测结果B通过朴素贝叶斯模型，得到模型预测结果2；具体步骤如下：
步骤7-1、设z＝{z1,z2}为一个待分类项，其中z表示用分类器分类之后的模型预测结果，z1是模型预测结果A的值，即用Subspace KNN分类器对雷达+音频片段模型里的雷达数据片段进行识别分类得到的预测结果，z2是模型预测结果B的值，即用Bagged Trees分类器对雷达+音频片段模型里的音频数据片段进行识别分类得到的预测结果；
步骤7-2、设D1＝{d1,d2,d3,d4}为类别的集合，共有4个类别：深度睡眠期(DEEP)、浅度睡眠期、快速眼动期(REM)、清醒期；
步骤7-3、分别计算条件概率：P(d1|z),P(d2|z),P(d3|z),P(d4|z)，获得睡眠分期的结果，各个睡眠类别下的条件概率计算过程如下：
(1)输入一个已知睡眠类别的待分类项集合，作为训练的样本集，统计各个特征属性在各个睡眠类别下的条件概率估计，即模型预测结果在各个睡眠类别下的条件概率，分别为：
P(z1|d1),P(z2|d1)；
P(z1|d2),P(z2|d2)；…；
P(z1|d4),P(z2|d4)
(2)如果各个特征属性条件独立，即睡眠分期的模型预测结果A和模型预测结果B相互独立，则由贝叶斯定理可以根据模型预测结果，来计算出当前处于某一睡眠阶段的概率：
用以下公式来求得最大值：
步骤7-4、如果P(dpp|z)＝max{P(d1|z),P(d2|z),...,P(dii|z)}，则z∈dpp，pp∈{1,2,3,4}，从而得到睡眠分期的模型预测结果2。
步骤8、将步骤6获得的模型预测结果1和步骤7获得的模型预测结果2按时序拼接，得到最终的睡眠分期结果。
所述将结果1和结果2按照实际的时间顺序t1,t2,...,t2nn-1,t2nn进行拼接，得到最终的分期结果，得到整个决策级融合系统的分期准确率，并且将最终结果与PSG标准睡眠分期结果比对。
本发明与现有技术相比，其显著优点为：1)本发明使用多个传感采集整夜睡眠数据，从不同的角度进行睡眠分期，分期结果更加科学；2)本发明的决策级融合的分类器只需训练单个传感器的数据，数据量小，通过决策直接对各个传感器的判定结果进行融合；3)本发明的决策级融合的容错能力强，也适用于异类传感器。
下面结合说明书附图对本发明做进一步描述。
附图说明
图1为本发明基于多传感器数据决策级融合的睡眠分期方法流程图。
图2为在多传感器数据决策级融合的睡眠分期过程中，用PSG为睡眠分期提供参考的流程图。
图3为通用的决策级融合模型。
图4为雷达回波信号处理流程图。
图5为音频信号处理流程图。
图6为雷达+音频片段模型和雷达剩余片段模型的构建方法图。
图7为决策级融合系统框架图。
图8为基于多传感器数据决策级融合的睡眠分期，与PSG的标准睡眠分期结果对比图，其中图(a)为PSG睡眠分期结果图，图(b)为基于多传感器数据决策级融合的睡眠分期结果图。
具体实施方式
结合附图，本发明的一种基于多传感器决策层融合的睡眠分期方法，包括以下步骤：
步骤1、测试者平躺在测试床上，将雷达传感器设置于人体的正上方，用于采集雷达回波信号；将音频传感器设置于朝向人脸的方向，用于采集呼吸和鼾声信号。
步骤2、获得雷达回波信号X(t)及音频信号Y(t)，由于呼吸和心跳信号的频率不同，设计不同频率段的滤波器对雷达回波信号X(t)进行滤波，从而获得呼吸信号和心跳信号，对雷达回波信号X(t)用幅值门限检测法作为预判的条件，检测出体动信号，由音频信号Y(t)进行降噪和事件检测获得呼吸信号Y1(t)和鼾声信号Y2(t)。
步骤3、对步骤2获得的雷达回波信号X(t)进行特征参数提取，根据呼吸信号得到每分钟呼吸次数N_bre、每分钟呼吸次数的方差S_bre2、每分钟呼吸信号幅度的差值积累Ada_bre、呼吸信号中的体动信号参数move_b，根据心跳信号得到每分钟心跳次数N_hrt、每分钟心跳次数的方差S_hrt2、每分钟心跳信号幅度的差值积累Ada_hrt、心跳信号中的体动信号参数move_h，睡眠分期相关特征REM期参数即REM、DEEP期参数即DEEP、样本熵即SampEn、时序特征at。
考虑到应用于人类的生命体征监测，必须基于不损害人体健康的前提下，结合临床医学对医疗设备的辐射和功耗要求。本系统中的雷达传感器的硬件部分由三个模块组成：一是射频模块，二是中频模块，三是基带模块。射频模块由发射天线、接收天线(一对平板天线)、频率合成器及接收机组成。中频模块包括中频滤波器和放大器以及ADC、FPGA和USB3.0。基带模块包括数据预处理和信号提取。
提取11个雷达回波信号特征和1个时序特征对睡眠阶段进行划分，12个特征具体如下：
(1)每分钟的呼吸次数N_bre：由相邻峰值点即波峰或波谷的差值计算每次呼吸时长，从而计算每分钟的呼吸次数N_bre；
(2)每分钟呼吸次数的方差S_bre2：
其中，N_bre(g)为第g分钟呼吸的次数，为整夜信号中每分钟呼吸次数的平均值，l_b为计算每分钟呼吸次数的方差的窗长；
(3)每分钟呼吸信号幅度的差值积累Ada_bre：
其中，A_bre(g)为第g分钟呼吸信号幅度，A_bre(g+1)为第g+1分钟呼吸信号幅度，N_i为每分钟内插值点数，取一个常数；
(4)呼吸信号中的体动信号参数：设置呼吸信号中体动信号的平均幅值move_b_thres为门限，判定超过门限的点数来计算体动参数，记为move_b；
(5)每分钟心跳次数N_hrt：计算相邻波谷值之间的时间差来计算每分钟心跳次数N_hrt；
(6)每分钟心跳次数的方差S_hrt2：
其中，N_hrt(g1)为第g1分钟心跳信号的幅值，为整夜心跳信号平均每分钟幅值，N_i为每分钟内插值点数；
(7)每分钟心跳信号幅值的差值积累Ada_hrt(g1)：
其中，A_hrt(g1)为第g1分钟心跳信号幅度，A_hrt(g1+1)为第g1+1分钟心跳信号幅度，N_i为每分钟内插值点数；
(8)心跳信号中的体动信号参数：设置心跳信号中体动信号的平均幅值move_h_thres为门限，判定超过门限的点数来计算体动参数，记为move_h；
(9)取某分钟内前后T秒的呼吸次数的差值之间计算关系的参数REM：
其中，表示当前第g3+g4分钟内前T秒内的呼吸次数，表示当前第g3+g4分钟内后T秒的呼吸次数，T的取值范围是25～30，λ为一个常数；
(10)用来区分深度睡眠阶段的重要特征参数DEEP：
其中，Abody表示睡眠中的体动信号的幅值，Arespir表示睡眠中的呼吸信号的幅值，DEEP表示体动信号在总的呼吸加体动信号总和中所占比例；
(11)样本熵SampEn：由雷达回波信号X(t)采样获得雷达回波信号的一维时间序列X(n)，即雷达回波信号时间序列，1≤n≤k_r，重构为相空间矢量X*(ww)，1≤ww≤k_r-(m*-1)，m*为嵌入维数，k_r为序列长度，定义相空间中任意两个矢量X*(w1)和X*(w2)的距离为：
其中，kk_r为一个常数，其取值范围是1≤kk_r≤m*-1，
求取样本熵的具体步骤为：
给定相似容限r，r＝hh*SD，hh是取值范围为0.1～0.25的常数，SD是
X(n)的标准差，对空间中每一个X*(w3)，计算的数目与总数k_r-m*+1的比值，记为
表示的数目；
对w3求平均值，记为
维数m*加1，重复上述步骤，得到
雷达回波信号时间序列X(n)的样本熵为：
实际中k_r为有限值，则样本熵估计为
(12)时序特征at：若{a1,a2,......,aN_a}为雷达回波信号特征组按照相对时间顺序排列的N_a个序号，对应于所处睡眠时间段，则at＝a/amax为睡眠分期的时序特征。
步骤4、音频信号中提取出用来对睡眠阶段进行划分的特征参数共有24个，其中23个音频信号特征分为两大类，一类是呼吸相关参数，另一类是鼾声相关参数，以及时序特征，
呼吸相关参数：
(1)每分钟呼吸次数N_bre'：取每分钟的呼吸次数；
(2)每分钟里前后内呼吸次数差值BVP：
BVP＝|N_bre'(s)-N_bre'(s-1)|
其中，N_bre'(s)为第s分钟内呼吸次数，N_bre'(s-1)为第s-1分钟内呼吸次数；
(3)呼吸周期BC，取两次呼吸事件之间的时间间隔；
(4)均方根误差RMSE：
其中，Y1(s1)为第s1个点的呼吸信号的幅值，为呼吸信号幅值的平均值，Ns为采样点数；
(5)每分钟呼吸次数的方差S_bre'2：
其中，N_bre'(s2)为第s2分钟的呼吸次数，为整夜呼吸信号每分钟呼吸次数的平均值，l_b'为计算呼吸次数的方差的窗长，l_b'取一个常数；
(6)能量差值累积EDA_bre：
其中，是第u个呼吸事件的能量，是第u+1个呼吸事件的能量，N_bre'是每分钟呼吸事件总数，Y1u(q)是第u个呼吸事件的第q个采样点对应的呼吸信号幅度，Y1(u+1)(q1)是第u+1个呼吸事件的第q1个采样点对应的呼吸信号幅度，N_p为第u个呼吸事件的总点数，N_p'为第u+1个呼吸事件的总点数；
(7)过零率：呼吸信号的过零点的次数N_zcr；
(8)第一共振峰Formant1；
(9)第二共振峰Formant2；
(10)第三共振峰Formant3；
(11)第一共振峰的方差Formant_var1；
(12)第二共振峰的方差Formant_var2；
(13)第三共振峰的方差Formant_var3；
鼾声相关参数：
(14)时间延迟给定一维的鼾声事件时间序列Y2(n)后，计算序列的自相关函数并作该函数关于时间的函数的图形，当函数值下降到初始值的时，对应的时间就是时间延迟；e取一个极小值；
(15)嵌入维数m：利用时间延迟将一维鼾声事件时间序列嵌入到m维空间中；
(16)LLE：根据时间延迟和嵌入维数m重构相空间W，对每个点Wj寻找其最近邻点Wj'，计算Wj到Wj'的距离dj(0)＝|Wj-Wj'|，其中，j为相空间中点的坐标；
对每个点Wj，计算其与最近邻点Wj'在第i步前向演化后的距离：
dj(i)＝|Wj+i-W′j+i|＝dj(0)×eLLE×i；
由下式算得最大李雅普诺夫指数：
其中，p_n为非零的dj(i)的个数；
(17)ApEn：将鼾声事件时间序列Y2(n)，1≤n≤k，重构为m维相空间矢量Y2'(i1)，1≤i1≤k-(m-1)，k为序列长度，定义相空间中任意两个矢量Y2'(j1)和Y2'(j2)的距离为：
d'j1j2＝max|Y2'(j1+kk)-Y2'(j2+kk)|，1≤j1,j2≤k-m+1,j1≠j2，其中，kk是一个常数，其取值范围是1≤kk≤m-1，
求取近似熵的具体步骤为：
给定相似容限r1，r1＝hh_a*SD_a，hh_a是取值范围为0.1～0.25的常数，SD_a是鼾声事件时间序列Y2(n)的标准差，对空间中每一个Y2'(i1)，计算d'j1j2＜r1的数目与总数k-m+1的比值，记为
N_dj1j2表示d'j1j2＜r1的数目；
将取对数，对i1求平均值，记为φm(r1)：
维数m加1，重复上述步骤，得到φm+1(r1)；
鼾声事件时间序列Y2(n)的近似熵为：
实际中，k为有限值，则近似熵估计为
ApEn(m,r1)＝φm(r1)-φm+1(r1)；
(18)统计嵌入维数m≥4的个数N_m；
(19)鼾声事件时间序列的维数D；
(20-21)去趋势波动分析的阿尔法系数一α1和去趋势波动分析的阿尔法系数二α2：对长度为k的鼾声事件时间序列Y2(n)，1≤n≤k，计算其累积差值：
其中，是时间序列的均值：
将y2(n)划分为不重叠的n2个长度为l的区间，l即为时间尺度(窗长)，n2为区间(或窗口)数量；
对每一段鼾声事件时间序列采用最小二乘法拟合出局部趋势y'2(n)；
剔除y2(n)中每个区间的局部趋势，并计算新序列的均方根：
改变窗口长度l的大小，重复上述步骤，画出纵坐标为log[F(n)]，横坐标为log(n)的曲线，则曲线的斜率为时间序列的标度指数α；
取窗口长度l1时对应的标度指数α为去趋势波动分析的阿尔法系数一α1，
取窗口长度l2时对应的标度指数α为去趋势波动分析的阿尔法系数二α2；
(22)样本熵SampEn*：将鼾声事件时间序列Y2(n)，1≤n≤k重构为m维相空间矢量Y2 *(i2)，1≤i2≤k-(m-1)，定义相空间中任意两个矢量Y2 *(j3)和Y2 *(j4)的距离为：
其中，kk_1是一个常数，其取值范围是1≤kk_1≤m-1，
求取样本熵的具体步骤为：
给定相似容限r1，对空间中每一个计算d* j3j4＜r1的数目与总数k-m+1的比值，记为
N_dj3j4表示d* j3j4＜r1的数目；
对i2求平均值，记为
维数m加1，重复上述步骤，得到Bm+1(r1)；
鼾声事件时间序列的样本熵为：
实际中k为有限值，则样本熵估计为：
(23)香农熵H：即信息熵，定义为离散随机事件的出现概率，令每个鼾声事件时间序列Y2(n)，对应的概率为p(Y2＝yn)，1≤n≤k，则随机变量的信息熵为
(24)时序特征bt：若{b1,b2,......,bN_b}为音频信号特征组按照相对时间顺序排列的N_b个序号，对应于所处睡眠时间段，则bt＝b/bmax为睡眠分期的时序特征。
步骤5、构建决策级融合系统模型包含雷达剩余片段模型和雷达+音频片段模型，由于音频传感器只有在音频信号的有声段才能提取到音频特征，所以整夜的音频特征是呈片段的，而雷达传感器是所有时段都可以采集到的雷达回波信号的，即整夜的雷达特征是连续的，融合系统的模型构建基本原理是当没有音频信号的时候单独采用雷达传感器的数据，当两个传感器都能获取信息时同时采用雷达和音频传感器的数据，为了准确的区分数据类型，将有音频数据的部分称为“音频片段”，与有音频片段的时间段对应的雷达数据称为“雷达片段”，没有音频数据的时间段对应的雷达数据称为“雷达剩余片段”，构建雷达+音频片段模型和雷达剩余片段模型的基本步骤如下：
(1)构建雷达+音频片段模型：从整夜的音频数据区间[t1,t2nn]里筛选音频传感器获取的数据，截取有用的音频片段1、音频片段2、……、音频片段NN，每个片段对应的时间区间分别为[t1,t2],[t3,t4],...,[t2nn-1,t2nn]，从整夜雷达数据[t1,t2nn]区间里截取的[t1,t2],[t3,t4],...,[t2nn-1,t2nn]区间里的数据，分别对应雷达数据片段1、雷达片段2、……、雷达片段NN，音频片段包含23个音频信号特征、1个时序特征共24个数据特征，雷达片段包含11个雷达回波信号特征、1个时序特征共12个数据特征，把雷达片段特征和音频片段特征拼接成新的融合特征片段：雷达+音频片段；
(2)构建雷达剩余片段模型：将原始雷达数据中不存在音频数据时间段的雷达数据进行分割，即[t2,t3],[t4,t5],...,[t2nn-2,t2nn-1]区间对应的雷达数据片段，分成雷达剩余片段1、雷达剩余片段2、……、雷达剩余片段NN-1，包含11个雷达回波信号特征、1个时序特征共12个数据特征。
步骤6、对雷达剩余片段模型里的12个数据特征用Subspace KNN分类器进行识别分类，得到模型预测结果1，对雷达+音频片段模型里的雷达片段数据的12个数据特征用Subspace KNN分类器进行识别分类，得到模型预测结果A，对雷达+音频片段模型里的音频片段数据的24个数据特征用Bagged Trees分类器进行识别分类，得到模型预测结果B；
步骤7、将模型预测结果A和模型预测结果B通过朴素贝叶斯模型，得到模型预测结果2，使用基于朴素贝叶斯模型的多传感器数据决策级融合框架对A和B进行决策，具体步骤如下：
(1)、设z＝{z1,z2}为一个待分类项，其中z表示用分类器分类之后的模型预测结果，z1是模型预测结果A的值，即用Subspace KNN分类器对雷达+音频片段模型里的雷达数据片段进行识别分类得到的预测结果，z2是模型预测结果B的值，即用Bagged Trees分类器对雷达+音频片段模型里的音频数据片段进行识别分类得到的预测结果；
(2)、设D1＝{d1,d2,d3,d4}为类别的集合，共有4个类别：深度睡眠期(DEEP)、浅度睡眠期、快速眼动期(REM)、清醒期；
(3)、分别计算条件概率：P(d1|z),P(d2|z),P(d3|z),P(d4|z)，获得睡眠分期的结果，各个睡眠类别下的条件概率计算过程如下：
i)输入一个已知睡眠类别的待分类项集合，作为训练的样本集，统计各个特征属性在各个睡眠类别下的条件概率估计，即模型预测结果在各个睡眠类别下的条件概率，分别为：
P(z1|d1),P(z2|d1)；
P(z1|d2),P(z2|d2)；…；
P(z1|d4),P(z2|d4)
ii)如果各个特征属性条件独立，即睡眠分期的模型预测结果A和模型预测结果B相互独立，则由贝叶斯定理可以根据模型预测结果，来计算出当前处于某一睡眠阶段的概率：
用以下公式来求得最大值：
(4)、如果P(dpp|z)＝max{P(d1|z),P(d2|z),...,P(dii|z)}，则z∈dpp，pp∈{1,2,3,4}，从而得到睡眠分期的模型预测结果2。
步骤8、将结果1和结果2按照实际的时间顺序t1,t2,...,t2nn-1,t2nn进行拼接，得到最终的整夜睡眠分期结果，并且将最终结果与PSG标准睡眠分期结果比对。
本发明的方法简单易行，准确率高，与实际情况相吻合。
实施例
本发明基于多传感器数据决策级融合的睡眠分期方法，包括以下步骤：
步骤1、测试者：一男性实验人员，身高177cm，年龄23岁，体重70kg，BMI指数22.34，实验之前要求实验人员尽量放松心态，不做剧烈运动，饮食规律。实验时该测试者平躺在测试床上，将雷达传感器设置于人体的正上方，用于采集雷达回波信号；将音频传感器设置于朝向人脸的方向，用于采集呼吸和鼾声信号；
步骤2、开启雷达传感器和音频传感器采集整夜数据，整夜数据中，雷达传感器数据时长为435分钟，有效的音频传感器数据时长为341分钟；
步骤3、将多导睡眠记录仪接入电脑，导出睡眠分期的结果，保存计算机端雷达及音频设备采集的信号数据，完成后续的数据处理及睡眠分期等工作。
步骤4、对步骤2采集到的雷达回波信号进行特征提取，其中每分钟呼吸次数的方差S_bre2中的窗长l_b取5；每分钟呼吸信号幅度的差值积累Ada_bre中的插值点数N_i取10；参数REM的T取30；λ取2；求取样本熵的hh的值取0.15；
步骤4、对步骤2采集到的音频信号进行特征提取，其中每分钟呼吸次数的方差S_bre'2中的窗长l_b'取5，求取样本熵和近似熵的hh_a的值取0.15，去趋势波动分析的α1系数的l1取13，α2系数的l2取49；
步骤6、根据雷达回波信号和音频信号的特征数，构建雷达+音频片段模型和雷达剩余片段模型，将特征送入两个模型里进行识别分类；
步骤5、用Subspace KNN分类器对雷达剩余片段模型里的雷达回波信号特征进行识别分类，得到模型预测结果1，用Subspace KNN分类器对雷达+音频片段模型里的雷达回波信号特征进行识别分类，得到模型预测结果A，其中Subspace KNN分类器的参数subspaceDimension的值设为6，参数Numbers of learners的值设为50；用Bagged Trees分类器对雷达+音频片段模型里的音频信号特征进行识别分类，得到模型预测结果B，其中Bagged Trees分类器的参数Numbers of learners的值设为5；
步骤6、用朴素贝叶斯模型对模型预测结果A和模型预测结果B进行决策，得到模型预测结果2，将模型预测结果1和模型预测结果2按时序拼接，得到该实验人员整夜睡眠分期的准确率为82.7％；
步骤7、以PSG结果为标准，对比PSG睡眠分期结果和基于多传感器数据决策级融合的睡眠分期预测结果，发现本研究方法与实际情况结果基本吻合。
从图8的结果可以看出，(a)为PSG睡眠分期图，作为睡眠分期标准，(b)预测的整夜睡眠分期图。图(a)中，REM表示快速眼动期、DEEP表示深度睡眠期，各个睡眠周期分别对应：2-深度睡眠期、3-浅度睡眠期，4-快速眼动期，5-清醒期。与图(a)相比，图(b)中，REM期和DEEP期分类的准确性相对较高，且第三个方框中出现了连续的REM周期，与标准睡眠分期图吻合度较高。

Claims (7)

1.一种基于多传感器数据决策级融合的睡眠分期方法，其特征在于，包括以下步骤：
步骤1、测试者平躺在测试床上，将雷达传感器设置于人体的正上方，用于采集雷达回波信号；将音频传感器设置于朝向人脸的方向，用于采集呼吸和鼾声信号；
步骤2、开启雷达传感器和音频传感器采集信号，获得雷达回波信号X(t)及音频信号Y(t)，之后对雷达回波信号X(t)进行滤波，获得呼吸信号和心跳信号，对雷达回波信号X(t)用幅值门限检测法作为预判的条件，检测出体动信号，对音频信号Y(t)进行降噪和事件检测获得呼吸信号Y1(t)和鼾声信号Y2(t)；
步骤3、对步骤2获得的雷达回波信号X(t)进行特征参数提取，根据呼吸信号得到每分钟呼吸次数N_bre、每分钟呼吸次数的方差S_bre2、每分钟呼吸信号幅度的差值积累Ada_bre、呼吸信号中的体动信号参数move_b，根据心跳信号得到每分钟心跳次数N_hrt、每分钟心跳次数的方差S_hrt2、每分钟心跳信号幅度的差值积累Ada_hrt、心跳信号中的体动信号参数move_h，睡眠分期相关特征REM期参数REM、DEEP期参数DEEP、样本熵SampEn、时序特征at；
步骤4、对步骤2获得的音频信号Y(t)进行特征参数提取，包括睡眠呼吸相关参数以及鼾声相关参数，呼吸相关参数为：每分钟呼吸的次数N_bre'、呼吸次数差值BVP、呼吸周期BC、均方根误差RMSE、每分钟内呼吸次数的方差S_bre'2、能量差值累积EDA_bre、过零率即N_zcr、共振峰1即Formant1、共振峰2即Formant2、共振峰3即Formant3、共振峰1的方差即Formant_var1、共振峰2的方差即Formant_var2、共振峰3的方差即Formant_var3；鼾声相关参数为：李雅普诺夫指数LLE、时间延迟嵌入维数m、近似熵ApEn、统计嵌入维数m≥4的个数N_m，鼾声事件时间序列的维数D、去趋势波动分析的阿尔法系数一α1、去趋势波动分析的阿尔法系数二α2、样本熵SampEn*、香农熵H、时序特征bt；
步骤5、根据步骤3获得的11个雷达回波信号特征以及时序特征at、步骤4获得的23个音频信号特征以及时序特征bt，构建决策级融合系统模型，包含雷达剩余片段模型和雷达+音频片段模型，雷达剩余片段模型是单独由雷达传感器的特征数据构建的，雷达+音频片段模型是由雷达传感器和音频传感器共同获取的特征数据构建的；
步骤6、在雷达剩余片段模型里，对雷达数据特征，进行识别分类得到模型预测结果1；在雷达+音频片段模型里，对雷达数据特征进行识别分类，得到模型预测结果A；对音频数据特征进行识别分类，得到模型预测结果B；
步骤7、将步骤6获得的模型预测结果A和模型预测结果B通过朴素贝叶斯模型，得到模型预测结果2；
步骤8、将步骤6获得的模型预测结果1和步骤7获得的模型预测结果2按时序拼接，得到最终的睡眠分期结果。
2.根据权利要求1所述的基于多传感器数据决策级融合的睡眠分期方法，其特征在于，步骤3中提取11个雷达回波信号特征以及1个时序特征对睡眠阶段进行划分，12个特征具体如下：
(1)每分钟的呼吸次数N_bre：由相邻峰值点即波峰或波谷的差值计算每次呼吸时长，从而计算每分钟的呼吸次数N_bre；
(2)每分钟呼吸次数的方差S_bre2：
其中，N_bre(g)为第g分钟呼吸的次数，为整夜信号中每分钟呼吸次数的平均值，l_b为计算每分钟呼吸次数的方差的窗长；
(3)每分钟呼吸信号幅度的差值积累Ada_bre：
其中，A_bre(g)为第g分钟呼吸信号幅度，A_bre(g+1)为第g+1分钟呼吸信号幅度，N_i为每分钟内插值点数，取一个常数；
(4)呼吸信号中的体动信号参数：设置呼吸信号中体动信号的平均幅值move_b_thres为门限，判定超过门限的点数来计算体动参数，记为move_b；
(5)每分钟心跳次数N_hrt：计算相邻波谷值之间的时间差来计算每分钟心跳次数N_hrt；
(6)每分钟心跳次数的方差S_hrt2：
其中，N_hrt(g1)为第g1分钟心跳信号的幅值，为整夜心跳信号平均每分钟幅值，N_i为每分钟内插值点数；
(7)每分钟心跳信号幅值的差值积累Ada_hrt(g1)：
其中，A_hrt(g1)为第g1分钟心跳信号幅度，A_hrt(g1+1)为第g1+1分钟心跳信号幅度，N_i为每分钟内插值点数；
(8)心跳信号中的体动信号参数：设置心跳信号中体动信号的平均幅值move_h_thres为门限，判定超过门限的点数来计算体动参数，记为move_h；
(9)取某分钟内前后T秒的呼吸次数的差值之间计算关系的参数REM：
其中，表示当前第g3+g4分钟内前T秒内的呼吸次数，表示当前第g3+g4分钟内后T秒的呼吸次数，T的取值范围是25～30，λ为一个常数；
(10)用来区分深度睡眠阶段的重要特征参数DEEP：
其中，Abody表示睡眠中的体动信号的幅值，Arespir表示睡眠中的呼吸信号的幅值，DEEP表示体动信号在总的呼吸加体动信号总和中所占比例；
(11)样本熵SampEn：由雷达回波信号X(t)采样获得雷达回波信号的一维时间序列X(n)，即雷达回波信号时间序列，1≤n≤k_r，重构为相空间矢量X*(ww)，1≤ww≤k_r-(m*-1)，m*为嵌入维数，k_r为序列长度，定义相空间中任意两个矢量X*(w1)和X*(w2)的距离为：
其中，kk_r为一个常数，其取值范围是1≤kk_r≤m*-1，
求取样本熵的具体步骤为：
给定相似容限r，r＝hh*SD，hh是取值范围为0.1～0.25的常数，SD是X(n)的标准差，对空间中每一个X*(w3)，计算的数目与总数k_r-m*+1的比值，记为
表示的数目；
对w3求平均值，记为
维数m*加1，重复上述步骤，得到
雷达回波信号时间序列X(n)的样本熵为：
实际中k_r为有限值，则样本熵估计为
(12)时序特征at：若{a1,a2,......,aN_a}为雷达回波信号特征组按照相对时间顺序排列的N_a个序号，对应于所处睡眠时间段，则at＝a/amax为睡眠分期的时序特征。
3.根据权利要求1所述的基于多传感器数据决策级融合的睡眠分期方法，其特征在于，步骤4中对音频信号Y(t)使用能熵比双门限法进行语音端点检测，检测出有效的呼吸事件Y1(t)和有效的鼾声事件Y2(t)，由Y1(t)和Y2(t)进行采样得到呼吸事件时间序列Y1(n)和鼾声事件时间序列Y2(n)，提取出23个音频信号特征以及1个时序特征对睡眠阶段进行划分，24个特征参数分别为：
呼吸相关参数：
(1)每分钟呼吸次数N_bre'：取每分钟的呼吸次数；
(2)每分钟里前后内呼吸次数差值BVP：
BVP＝|N_bre'(s)-N_bre'(s-1)|
其中，N_bre'(s)为第s分钟内呼吸次数，N_bre'(s-1)为第s-1分钟内呼吸次数；
(3)呼吸周期BC，取两次呼吸事件之间的时间间隔；
(4)均方根误差RMSE：
其中，Y1(s1)为第s1个点的呼吸信号的幅值，为呼吸信号幅值的平均值，Ns为采样点数；
(5)每分钟呼吸次数的方差S_bre'2：
其中，N_bre'(s2)为第s2分钟的呼吸次数，为整夜呼吸信号每分钟呼吸次数的平均值，l_b'为计算呼吸次数的方差的窗长，l_b'取一个常数；
(6)能量差值累积EDA_bre：
其中，是第u个呼吸事件的能量，是第u+1个呼吸事件的能量，N_bre'是每分钟呼吸事件总数，Y1u(q)是第u个呼吸事件的第q个采样点对应的呼吸信号幅度，Y1(u+1)(q1)是第u+1个呼吸事件的第q1个采样点对应的呼吸信号幅度，N_p为第u个呼吸事件的总点数，N_p'为第u+1个呼吸事件的总点数；
(7)过零率：呼吸信号的过零点的次数N_zcr；
(8)第一共振峰Formant1；
(9)第二共振峰Formant2；
(10)第三共振峰Formant3；
(11)第一共振峰的方差Formant_var1；
(12)第二共振峰的方差Formant_var2；
(13)第三共振峰的方差Formant_var3；
鼾声相关参数：
(14)时间延迟给定一维的鼾声事件时间序列Y2(n)后，计算序列的自相关函数并作该函数关于时间的函数的图形，当函数值下降到初始值的时，对应的时间就是时间延迟；e取一个极小值；
(15)嵌入维数m：利用时间延迟将一维鼾声事件时间序列嵌入到m维空间中；
(16)LLE：根据时间延迟和嵌入维数m重构相空间W，对每个点Wj寻找其最近邻点Wj'，计算Wj到Wj'的距离dj(0)＝|Wj-Wj'|，其中，j为相空间中点的坐标；
对每个点Wj，计算其与最近邻点Wj'在第i步前向演化后的距离：
dj(i)＝|Wj+i-W′j+i|＝dj(0)×eLLE×i；
由下式算得最大李雅普诺夫指数：
其中，p_n为非零的dj(i)的个数；
(17)ApEn：将鼾声事件时间序列Y2(n)，1≤n≤k，重构为m维相空间矢量Y2'(i1)，1≤i1≤k-(m-1)，k为序列长度，定义相空间中任意两个矢量Y2'(j1)和Y2'(j2)的距离为：
d'j1j2＝max|Y2'(j1+kk)-Y2'(j2+kk)|，1≤j1,j2≤k-m+1,j1≠j2，其中，kk是一个常数，其取值范围是1≤kk≤m-1，
求取近似熵的具体步骤为：
给定相似容限r1，r1＝hh_a*SD_a，hh_a是取值范围为0.1～0.25的常数，SD_a是鼾声事件时间序列Y2(n)的标准差，对空间中每一个Y2'(i1)，计算d'j1j2＜r1的数目与总数k-m+1的比值，记为
N_dj1j2表示d'j1j2＜r1的数目；
将取对数，对i1求平均值，记为φm(r1)：
维数m加1，重复上述步骤，得到φm+1(r1)；
鼾声事件时间序列Y2(n)的近似熵为：
实际中，k为有限值，则近似熵估计为
ApEn(m,r1)＝φm(r1)-φm+1(r1)；
(18)统计嵌入维数m≥4的个数N_m；
(19)鼾声事件时间序列的维数D；
(20-21)去趋势波动分析的阿尔法系数一α1和去趋势波动分析的阿尔法系数二α2：对长度为k的鼾声事件时间序列Y2(n)，1≤n≤k，计算其累积差值：
其中，是时间序列的均值：
将y2(n)划分为不重叠的n2个长度为l的区间，l即为时间尺度(窗长)，n2为区间(或窗口)数量；
对每一段鼾声事件时间序列采用最小二乘法拟合出局部趋势y'2(n)；
剔除y2(n)中每个区间的局部趋势，并计算新序列的均方根：
改变窗口长度l的大小，重复上述步骤，画出纵坐标为log[F(n)]，横坐标为log(n)的曲线，则曲线的斜率为时间序列的标度指数α；
取窗口长度l1时对应的标度指数α为去趋势波动分析的阿尔法系数一α1，
取窗口长度l2时对应的标度指数α为去趋势波动分析的阿尔法系数二α2；
(22)样本熵SampEn*：将鼾声事件时间序列Y2(n)，1≤n≤k重构为m维相空间矢量Y2 *(i2)，1≤i2≤k-(m-1)，定义相空间中任意两个矢量Y2 *(j3)和Y2 *(j4)的距离为：
1≤j3,j4≤k-m+1,j3≠j4，其中，kk_1是一个常数，其取值范围是1≤kk_1≤m-1，
求取样本熵的具体步骤为：
给定相似容限r1，对空间中每一个计算d* j3j4＜r1的数目与总数k-m+1的比值，记为
N_dj3j4表示d* j3j4＜r1的数目；
对i2求平均值，记为
维数m加1，重复上述步骤，得到Bm+1(r1)；
鼾声事件时间序列的样本熵为：
实际中k为有限值，则样本熵估计为：
(23)香农熵H：即信息熵，定义为离散随机事件的出现概率，令每个鼾声事件时间序列Y2(n)，对应的概率为p(Y2＝yn)，1≤n≤k，则随机变量的信息熵为
(24)时序特征bt：若{b1,b2,......,bN_b}为音频信号特征组按照相对时间顺序排列的N_b个序号，对应于所处睡眠时间段，则bt＝b/bmax为睡眠分期的时序特征。
4.根据权利要求1所述的基于多传感器数据决策级融合的睡眠分期方法，其特征在于，步骤5中构建雷达+音频片段模型和雷达剩余片段模型的基本步骤如下：
步骤5-1、构建雷达+音频片段模型：从整夜的音频数据区间[t1,t2nn]里筛选音频传感器获取的数据，截取有用的音频片段1、音频片段2、……、音频片段NN，每个片段对应的时间区间分别为[t1,t2],[t3,t4],...,[t2nn-1,t2nn]，从整夜雷达数据[t1,t2nn]区间里截取的[t1,t2],[t3,t4],...,[t2nn-1,t2nn]区间里的数据，分别对应雷达数据片段1、雷达片段2、……、雷达片段NN，音频片段包含23个音频特征和1个时序特征共24个特征，雷达片段包含11个雷达回波信号特征和1个时序特征共12个数据特征，把雷达片段特征和音频片段特征拼接成新的融合特征片段：雷达+音频片段；
步骤5-2、构建雷达剩余片段模型：将原始雷达数据中不存在音频数据时间段的雷达数据进行分割，即[t2,t3],[t4,t5],...,[t2nn-2,t2nn-1]区间对应的雷达数据片段，分成雷达剩余片段1、雷达剩余片段2、……、雷达剩余片段NN-1，包含11个雷达回波信号特征、1个时序特征共12个数据特征；
雷达+音频片段模型和雷达剩余片段模型表如下：
由雷达+音频片段1、雷达+音频片段2、……、雷达+音频片段NN构成的模型为雷达+音频片段模型，由雷达剩余片段1、雷达剩余片段2、……、雷达剩余片段NN-1构成的模型即为雷达剩余片段模型。
5.根据权利要求1所述的基于多传感器数据决策级融合的睡眠分期方法，其特征在于，步骤6中对雷达剩余片段模型和雷达+音频片段模型里的数据特征进行识别分类，得到模型预测结果，具体为：雷达剩余片段模型里的12个数据特征，用Subspace KNN分类器进行识别分类，得到模型预测结果1，对雷达+音频片段模型里的雷达片段数据的12个数据特征，用Subspace KNN分类器进行识别分类，得到模型预测结果A，对雷达+音频片段模型里的音频片段数据的24个数据特征，用Bagged Trees分类器进行识别分类，得到模型预测结果B。
6.根据权利要求1所述的基于多传感器数据决策级融合的睡眠分期方法，其特征在于，步骤7中将模型预测结果A和模型预测结果B通过朴素贝叶斯模型，得到模型预测结果2，具体步骤如下：
步骤7-1、设z＝{z1,z2}为一个待分类项，其中z表示用分类器分类之后的模型预测结果，z1是模型预测结果A的值，即用Subspace KNN分类器对雷达+音频片段模型里的雷达数据片段进行识别分类得到的预测结果，z2是模型预测结果B的值，即用Bagged Trees分类器对雷达+音频片段模型里的音频数据片段进行识别分类得到的预测结果；
步骤7-2、设D1＝{d1,d2,d3,d4}为类别的集合，共有4个类别：深度睡眠期(DEEP)、浅度睡眠期、快速眼动期(REM)、清醒期；
步骤7-3、分别计算条件概率：P(d1|z),P(d2|z),P(d3|z),P(d4|z)，获得睡眠分期的结果，各个睡眠类别下的条件概率计算过程如下：
(1)输入一个已知睡眠类别的待分类项集合，作为训练的样本集，统计各个特征属性在各个睡眠类别下的条件概率估计，即模型预测结果在各个睡眠类别下的条件概率，分别为：
P(z1|d1),P(z2|d1)；
P(z1|d2),P(z2|d2)；…；
P(z1|d4),P(z2|d4)
(2)如果各个特征属性条件独立，即睡眠分期的模型预测结果A和模型预测结果B相互独立，则由贝叶斯定理可以根据模型预测结果，来计算出当前处于某一睡眠阶段的概率：
用以下公式来求得最大值：
步骤7-4、如果P(dpp|z)＝max{P(d1|z),P(d2|z),...,P(dii|z)}，则z∈dpp，pp∈{1,2,3,4}，从而得到睡眠分期的模型预测结果2。
7.根据权利要求1所述的基于多传感器数据决策级融合的睡眠分期方法，其特征在于，步骤8中将结果1和结果2按照实际的时间顺序t1,t2,...,t2nn-1,t2nn进行拼接，得到最终的分期结果，得到整个决策级融合系统的分期准确率，并且将最终结果与PSG标准睡眠分期结果比对。
CN201810125718.1A 2018-02-08 2018-02-08 基于多传感器数据决策级融合的睡眠分期方法 Active CN108416367B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN201810125718.1A CN108416367B (zh) 2018-02-08 2018-02-08 基于多传感器数据决策级融合的睡眠分期方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN201810125718.1A CN108416367B (zh) 2018-02-08 2018-02-08 基于多传感器数据决策级融合的睡眠分期方法

Publications (2)

Publication Number Publication Date
CN108416367A true CN108416367A (zh) 2018-08-17
CN108416367B CN108416367B (zh) 2021-12-10

Family

ID=63128065

Family Applications (1)

Application Number Title Priority Date Filing Date
CN201810125718.1A Active CN108416367B (zh) 2018-02-08 2018-02-08 基于多传感器数据决策级融合的睡眠分期方法

Country Status (1)

Country Link
CN (1) CN108416367B (zh)

Cited By (14)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN109350075A (zh) * 2018-09-18 2019-02-19 深圳和而泰数据资源与云技术有限公司 一种鼾声检测方法、装置及可读存储介质
CN109523011A (zh) * 2018-11-06 2019-03-26 哈尔滨工业大学（深圳） 一种面向多无人机协同探测的多传感器自适应管理方法
CN109864706A (zh) * 2019-01-08 2019-06-11 廖旺才 用户状态信息的提示方法、装置、计算机设备和存储介质
CN110367986A (zh) * 2019-06-25 2019-10-25 南京理工大学 基于ewt的呼吸信号近似熵特征提取方法
CN110693454A (zh) * 2019-08-23 2020-01-17 深圳大学 基于雷达的睡眠特征事件检测方法、装置及存储介质
CN111603162A (zh) * 2020-05-07 2020-09-01 北京海益同展信息科技有限公司 肌电信号处理方法、装置、智能可穿戴设备及存储介质
CN112066528A (zh) * 2020-08-10 2020-12-11 珠海格力电器股份有限公司 一种空调控制方法、装置、存储介质及空调
CN112168139A (zh) * 2019-07-05 2021-01-05 腾讯科技（深圳）有限公司 一种健康监控方法、装置及存储介质
CN112450881A (zh) * 2020-11-12 2021-03-09 武汉大学 一种基于时序关联性驱动的多模态睡眠分期方法
CN113647909A (zh) * 2021-08-09 2021-11-16 珠海格力电器股份有限公司 睡眠指标监测方法、电子装置、装置及可读存储介质
CN114176521A (zh) * 2021-12-09 2022-03-15 森思泰克河北科技有限公司 基于雷达信息的睡眠分期方法、装置及终端
CN114224320A (zh) * 2021-12-31 2022-03-25 深圳融昕医疗科技有限公司 自适应多通道信号融合的鼾声检测方法、设备及系统
CN114403847A (zh) * 2021-12-17 2022-04-29 中南民族大学 基于腹部和肺部数据相关性的呼吸状态检测方法及系统
CN116369864A (zh) * 2023-06-02 2023-07-04 长春中医药大学 基于数据编码的睡眠监测数据智能管理方法和系统

Citations (6)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US20120088992A1 (en) * 2004-12-23 2012-04-12 Jeffrey Armitstead Method for detecting and discriminating breathing patterns from respiratory signals
KR20140006256A (ko) * 2012-06-29 2014-01-16 전자부품연구원 바이오 레이더 기반 수면 무호흡 모니터링 시스템 및 방법
CN103635134A (zh) * 2011-05-11 2014-03-12 皇家飞利浦有限公司 睡眠阶段注释设备
CN104794505A (zh) * 2015-04-28 2015-07-22 上海大学 一种多通道的脑电信号数据融合降维方法
CN107280674A (zh) * 2017-06-02 2017-10-24 南京理工大学 基于音频录取设备的呼吸模式判决方法
CN107307846A (zh) * 2016-04-27 2017-11-03 南京理工大学 非接触式睡眠分期方法

Patent Citations (6)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US20120088992A1 (en) * 2004-12-23 2012-04-12 Jeffrey Armitstead Method for detecting and discriminating breathing patterns from respiratory signals
CN103635134A (zh) * 2011-05-11 2014-03-12 皇家飞利浦有限公司 睡眠阶段注释设备
KR20140006256A (ko) * 2012-06-29 2014-01-16 전자부품연구원 바이오 레이더 기반 수면 무호흡 모니터링 시스템 및 방법
CN104794505A (zh) * 2015-04-28 2015-07-22 上海大学 一种多通道的脑电信号数据融合降维方法
CN107307846A (zh) * 2016-04-27 2017-11-03 南京理工大学 非接触式睡眠分期方法
CN107280674A (zh) * 2017-06-02 2017-10-24 南京理工大学 基于音频录取设备的呼吸模式判决方法

Non-Patent Citations (1)

* Cited by examiner, † Cited by third party
Title
JUNJUN XIONG ET AL: "Sleep Apnea Detection with Doppler Radar Based", 《ACES》 *

Cited By (24)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN109350075A (zh) * 2018-09-18 2019-02-19 深圳和而泰数据资源与云技术有限公司 一种鼾声检测方法、装置及可读存储介质
CN109350075B (zh) * 2018-09-18 2022-02-01 深圳和而泰数据资源与云技术有限公司 一种鼾声检测方法、装置及可读存储介质
CN109523011B (zh) * 2018-11-06 2021-07-27 哈尔滨工业大学（深圳） 一种面向多无人机协同探测的多传感器自适应管理方法
CN109523011A (zh) * 2018-11-06 2019-03-26 哈尔滨工业大学（深圳） 一种面向多无人机协同探测的多传感器自适应管理方法
CN109864706A (zh) * 2019-01-08 2019-06-11 廖旺才 用户状态信息的提示方法、装置、计算机设备和存储介质
CN110367986A (zh) * 2019-06-25 2019-10-25 南京理工大学 基于ewt的呼吸信号近似熵特征提取方法
CN112168139B (zh) * 2019-07-05 2022-09-30 腾讯科技（深圳）有限公司 一种健康监控方法、装置及存储介质
CN112168139A (zh) * 2019-07-05 2021-01-05 腾讯科技（深圳）有限公司 一种健康监控方法、装置及存储介质
CN110693454A (zh) * 2019-08-23 2020-01-17 深圳大学 基于雷达的睡眠特征事件检测方法、装置及存储介质
CN110693454B (zh) * 2019-08-23 2023-04-25 深圳大学 基于雷达的睡眠特征事件检测方法、装置及存储介质
CN111603162A (zh) * 2020-05-07 2020-09-01 北京海益同展信息科技有限公司 肌电信号处理方法、装置、智能可穿戴设备及存储介质
CN112066528A (zh) * 2020-08-10 2020-12-11 珠海格力电器股份有限公司 一种空调控制方法、装置、存储介质及空调
CN112450881A (zh) * 2020-11-12 2021-03-09 武汉大学 一种基于时序关联性驱动的多模态睡眠分期方法
CN112450881B (zh) * 2020-11-12 2021-11-02 武汉大学 一种基于时序关联性驱动的多模态睡眠分期方法
CN113647909A (zh) * 2021-08-09 2021-11-16 珠海格力电器股份有限公司 睡眠指标监测方法、电子装置、装置及可读存储介质
CN113647909B (zh) * 2021-08-09 2022-05-17 珠海格力电器股份有限公司 睡眠指标监测方法、电子装置、装置及可读存储介质
CN114176521A (zh) * 2021-12-09 2022-03-15 森思泰克河北科技有限公司 基于雷达信息的睡眠分期方法、装置及终端
CN114176521B (zh) * 2021-12-09 2022-07-01 森思泰克河北科技有限公司 基于雷达信息的睡眠分期方法、装置及终端
CN114403847A (zh) * 2021-12-17 2022-04-29 中南民族大学 基于腹部和肺部数据相关性的呼吸状态检测方法及系统
CN114403847B (zh) * 2021-12-17 2022-11-11 中南民族大学 基于腹部和肺部数据相关性的呼吸状态检测方法及系统
CN114224320B (zh) * 2021-12-31 2023-01-17 深圳融昕医疗科技有限公司 自适应多通道信号融合的鼾声检测方法、设备及系统
CN114224320A (zh) * 2021-12-31 2022-03-25 深圳融昕医疗科技有限公司 自适应多通道信号融合的鼾声检测方法、设备及系统
CN116369864A (zh) * 2023-06-02 2023-07-04 长春中医药大学 基于数据编码的睡眠监测数据智能管理方法和系统
CN116369864B (zh) * 2023-06-02 2023-08-08 长春中医药大学 基于数据编码的睡眠监测数据智能管理方法和系统

Also Published As

Publication number Publication date
CN108416367B (zh) 2021-12-10

Similar Documents

Publication Publication Date Title
CN108416367A (zh) 基于多传感器数据决策级融合的睡眠分期方法
CN108388912A (zh) 基于多传感器特征优化算法的睡眠分期方法
CN106725488B (zh) 一种无线场强呼吸检测方法、装置及呼吸检测仪
JP5271718B2 (ja) 複数の時間区分にわたり胎児及び産婦のｅｃｇを特定する方法
Castro et al. Heart sound segmentation of pediatric auscultations using wavelet analysis
CN100418480C (zh) 基于心音分析的心脏病自动分类系统及其心音分段方法
CN107249449B (zh) 用于感测胎儿活动的系统、装置和方法
CN112998701A (zh) 基于毫米波雷达的生命体征检测与身份识别系统和方法
CN112998690B (zh) 一种基于脉搏波多特征融合的呼吸率提取方法
CN107530015B (zh) 一种生命体征分析方法与系统
JP2018512243A5 (zh)
WO2018149147A1 (zh) 一种呼吸率提取方法及装置
CN110960211A (zh) 一种基于嵌入式的心电实时监测系统
Subramaniam et al. A review of significant researches on prediction of preterm birth using uterine electromyogram signal
Casson et al. Toward online data reduction for portable electroencephalography systems in epilepsy
CN114916913B (zh) 一种便携式睡眠呼吸状态实时监测系统及方法
CN115736854A (zh) 一种基于毫米波雷达的呼吸心跳监测系统
CN115474901A (zh) 基于无线射频信号的非接触式起居状态监测方法及系统
CN206979487U (zh) 一种无线场强呼吸检测仪
CN106419884B (zh) 一种基于小波分析的心率计算方法与系统
CN106725487A (zh) 一种呼吸率提取方法及装置
Chatterjee et al. Real–time detection of electrocardiogram wave features using template matching and implementation in FPGA
Bonizzi et al. Sleep apnea detection directly from unprocessed ECG through singular spectrum decomposition
Kaya et al. Abnormality detection in noisy biosignals
Lu et al. Pulse waveform analysis for pregnancy diagnosis based on machine learning

Legal Events

Date Code Title Description
PB01 Publication
PB01 Publication
SE01 Entry into force of request for substantive examination
SE01 Entry into force of request for substantive examination
GR01 Patent grant
GR01 Patent grant