本发明属于雷达技术和声学领域，具体涉及利用多传感器采集数据，提取特征，在特征级别上进行数据融合，进行睡眠分期。
背景技术：
在生物医学工程与大数据等相关技术逐渐应用于临床后，人们发现整夜睡眠信息能在一定程度上映射出人体的健康状况，从而进行疾病的预防和治疗。在睡眠过程中，人的脑电图、心电图等多种生理信号发生变化，而这些信号变化随着入睡的深浅度有一定的相关性变化。医学上根据脑电波、心电、眼电等生理表现人为的将睡眠划分为几个阶段，实际上各个睡眠阶段之间并没有很明确的界限，而是周期之间逐渐转变且有承接性的，整夜的睡眠阶段也是呈现周期性的，有一定的规律。在整夜睡眠中，一般会出现4～6个睡眠周期，前后周期之间有一定的关联性，周而复始、循环往复。
随着人们对健康要求标准的不断提高，睡眠阶段所反映的健康状况也越来越受到重视。睡眠监测产品多采用接触式监测，需要在受试者身体的多个部位贴数个电极，佩戴口鼻气流管和腹部绑带。接触式的监测在很大程度上影响了受试者的生理和心理情况，且测试环境复杂，需要专业医护人员操作，价格相对昂贵。这种方法需要专业人员操作，价格昂贵，测试者需要直接接触仪器，极易产生体感不适。在非接触式的睡眠分期监测方法中患者不需要与任何仪器接触、操作简单，可长期监测，成本较低。但是现有技术中，尚无相关描述。
技术实现要素：
本发明的目的在于提供一种基于多传感器特征优化算法的睡眠分期方法，提高了非接触式睡眠监测产品准确率。
实现本发明目的的技术解决方案为：一种基于多传感器特征优化算法的睡眠分期方法，包括以下步骤：
步骤1、测试者平躺在测试床上，雷达传感器设置于距离人体胸部正上方，用于采集雷达回波信号；音频传感器设置于朝向人脸的方向，用于采集呼吸和鼾声信号；
床上观测区域包括：天花板、床面，天花板和床面相平行，测试者摆平躺在床面上，雷达设置在天花板的中间，雷达天线朝向床面，正对测试者胸腔；雷达采用连续波雷达；音频传感器设置在天花板的中间，麦克风正对测试者头部。
步骤2、同时开启雷达传感器和音频传感器采集信号，测试者入睡，雷达照射到人体胸腔表面形成回波，测试者睡眠产生的鼾声被音频传感器记录；
雷达对物体进行检测具体为：连续波雷达对测试者胸腔进行照射，雷达发射的脉冲信号经过人体胸腔的反射，反射的回波信号中携带着相应的生命体征信息，再由雷达接收机接收到回波信号。
步骤3、利用数字信号处理方法对原始回波进行处理，具体为：对雷达信号进行反正切解调，提取出体动信号、呼吸和心跳信号，使用降噪算法对原始音频进行处理；
步骤4、对步骤3得到的呼吸、体动、心跳提取23种特征并对音频信号提取11种特征，对音频信号的11种特征进行拼接融合，对融合后的数据进行判断，以有无音频特征分两个模型进行判决；
11种特征以分钟为单位进行计算，但是以每t秒为间隔提取特征，所述t的取值范围为20<t<35；
(1)rpm
rpm为每分钟的呼吸次数；
(2)rpm_var
rpm_var为每分钟呼吸次数的方差，设每分钟的呼吸次数为r(n_bre)，整夜信号中每分钟呼吸次数的平均值为
(3)rpm_ada
rpm_ada为每分钟呼吸信号幅度的差值积累，a_bre(n_bre)为呼吸信号幅度，n为每分钟内插值点数，每分钟呼吸信号幅度的差值积累ada的计算公式如下：
(4)rpm_move
rpm_move为呼吸信号中的体动信号特征；
(5)bpm
bpm为每分钟心跳次数；
(6)bpm_var
bpm_var为每分钟心跳次数的方差，b_bre(n_bre_var)为心跳次数，
(7)bpm_ada
bpm_ada为每分钟心跳信号幅值的差值积累，其中a_ada＇(n_bre_ada)为心跳信号幅值，n为每分钟内插值点数，
(8)bpm_move
bpm_move为心跳信号中的体动信号特征，根据心跳信号中体动信号的幅值设定门限，判定超过门限的点数来计算体动特征；每分钟取n个内插点，n取值范围为90<n<110，当前内插点的信号幅值超过门限，且当前内插点前后t秒内信号幅值超过门限的点数超过50个的，此内插点记为一次心跳信号中的体动点，统计t秒内心跳信号的体动点个数记为心跳信号的体动信号特征，t秒取值20<t<35；
(9)rem
rem是取某分钟内前后t秒的呼吸次数的差值之间计算关系的特征，rem具体计算公式如下：
其中，
(10)deep
体动信号deep特征计算公式如下：
其中，
(11)sampen
由雷达回波信号x(t)采样获得雷达回波信号的一维时间序列x(n)，1≤n≤k_r重构为相空间矢量x*(ww)，1≤ww≤k-(m-1)，m为嵌入维数，定义相空间中任意两个矢量x*(w1)和x*(w2)的距离为：
求取样本熵的具体步骤为：
给定相似容限r，r＝hh*sd，hh是取值范围为0.1～0.25的常数，sd是雷达回波信号时间序列x(n)的标准差，对空间中每一个x*(w3)，进行模板匹配：
对w3求平均值，记为δm(r)：
维数m加1，重复上述步骤，得到gm+1(r)；
雷达回波信号时间序列x(n)的样本熵为：
实际中k为有限值，则样本熵估计为
所述23类特征分为两大类，一类是睡眠呼吸相关特征，另一类是鼾声的线性和非线性特征：
a.呼吸相关特征，该类特征共13个，分别是：
1)rpm：每分钟内呼吸次数；
2)bvp：前后t秒内呼吸次数差值，b_bvp(n_bvp)为当前t秒内呼吸次数，b_bvp(n_bvp-1)为前t秒内呼吸次数，b_bvp(n_bvp+1)为后t秒内呼吸次数，前后t秒内呼吸次数差值的计算公式如下：
bvp＝|b_bvp(n_bvp)-b_bvp(n_bvp-1)|
+|b_bvp(n_bvp)-b_bvp(n_bvp+1)|
其中t秒取值20<t<35；
3)bc：呼吸周期，取两次呼吸事件之间的时间差均值；
4)rmse：均方根误差，根据每个音频片段的采样点数，求出该片段内的均方根误差；x_rmse(n_rmse)为当前呼吸信号，
5)rpm_var：每分钟内呼吸次数的方差；xformer为前t秒内呼吸次数，xlatter为后t秒内呼吸次数，t秒取值20<t<35，计算公式如下：
rpm_var＝xformer-xlatter
6)eda：能量差值累积；根据每个音频片段内的采样点数，求出相邻片段的能量累积差值和，x1(n)与x2(n)为相邻片段，n为音频片段内采样点数，计算公式如下：
7)cross_zr：过零率，每分钟内呼吸信号的过零点的次数；
8-10)formant1、formant2、formant3：3个共振峰，共振峰是指在声音信号的频谱图中能量比较集中的区域；取能量较强部分的音频信息有助于获得该片段的特征标记，有利于睡眠阶段的区分；
11-13)formant_var1、formant_var2、formant_var3：分别取(8-10)中3个共振峰的方差；
b.鼾声相关的线性和非线性特征
1)lle：lyapunov李雅普诺夫指数；lyapunov李雅普诺夫指数给出的是动态系统沿其相空间主轴发散或收敛的平均速度；
根据求得的时间延迟
对每个点wj，计算其与最近邻点w′j在第i步前向演化后的距离
dj(i)＝|wj+i-w′j+i|＝dj(0)×eλ×i；
由下式算得最大李雅普诺夫指数,
2)时间延迟
3)嵌入维数m：嵌入维数，是相空间重构中的两个重构参数之一，相空间重构中计算李雅普诺夫指数的参数；利用时间延迟
4)apen：近似熵，是一种计算序列的复杂性和统计量化的参数；
将鼾声事件时间序列的一维时间序列x(n)＝(x1,x2,x3,...,xi,...,xk)，按顺序组成v维矢量vi＝[x(i),x(i+1),...,x(i+v-1)]，i＝1,2,...,k,其中k为鼾声事件时间序列x(n)的长度；对每一个i值计算相矢量vi和其余矢量vj的距离；
dij＝max|x(i+1)-x(j+1)|,1＝0,1,...,v-1
给定阈值r＝a3×sd，其中a3的取值范围为0.1～0.25，sd为鼾声事件时间序列x(n)的标准差；记录每个dij小于阈值r对应的i的个数，并求出与总的v维相矢量个数(k-v+1)的比值，记为
x(n)的近似熵为
apen＝φv(r)-φv+1(r)
5)n：统计嵌入维数m>＝4的个数。
6)d：时间序列的维数；
7-8)去趋势波动分析的阿尔法系数一α1和去趋势波动分析的阿尔法系数二α2：对长度为k的鼾声事件时间序列y2(n)，1≤n≤k，计算其累积差值：
其中，
将y2(n)划分为不重叠的n2个长度为l的区间，l即为时间尺度，n2为区间数量；
对每一段鼾声事件时间序列采用最小二乘法拟合出局部趋势y′2(n)；
剔除y2(n)中每个区间的局部趋势，并计算新序列的均方根：
改变窗口长度l的大小，重复上述步骤，则存在幂律形式的波动：f(n)∝nα画出纵坐标为log[f(n)]，横坐标为log(n)的曲线，则曲线的斜率为时间序列的标度指数α；
取窗口长度为l1时对应的标度指数α为去趋势波动分析的阿尔法系数一α1；
取窗口长度为l2时对应的标度指数α为去趋势波动分析的阿尔法系数二α2；9)样本熵sampen：将鼾声事件时间序列y2(n)，1≤n≤k重构为m维相空间矢量
求取样本熵的具体步骤为：
给定相似容限r1，对空间中每一个
对i2求平均值，记为
维数m加1，重复上述步骤，得到bm+1(r1)；
鼾声事件时间序列的样本熵为：
实际中k为有限值，则样本熵估计为：
10)香农熵h：即信息熵，定义为离散随机事件的出现概率，令每个鼾声事件时间序列y＝{y1,y2…,y3}，对应的概率为p(y＝yi)，则随机变量的信息熵为
时序特征at：若{a1,a2,……,an}为雷达或是音频信号特征组按照相对时间顺序排列的n个序号，对应于所处睡眠时间段，则at＝a/amax为睡眠分期的时序特征。
步骤5、对步骤4得到的特征模型分别加入时序特征，构建特征级融合系统模型，包含雷达剩余片段模型和雷达+音频片段模型，雷达剩余片段模型是单独由雷达传感器的特征数据构建的，雷达+音频片段模型是由雷达传感器和音频传感器共同获取的特征数据构建的；
构建特征级融合系统模型，具体为：
步骤5-1、分别构建两个模型：
(1)雷达剩余片段模型：单独由雷达传感器的特征数据构建的，将雷达特征数据对应有音频特征数据的时间段进行分割，分为雷达剩余片段{t1,t2},{t2,t3}…{tn,tn+1}和雷达片段{t1,t2},{t2,t3}…{tn,tn+1}，用雷达剩余片段{t1,t2},{t2,t3}…{tn,tn+1}训练雷达剩余片段模型；
(2)雷达+音频片段模型：由雷达传感器和音频传感器共同获取的特征数据构建的，将同时有雷达和音频的特征数据融合称为雷达+音频片段，用雷达+音频片段{t1,t2},{t2,t3}…{tn,tn+1}训练雷达+音频片段模型；
步骤5-2、构建特征级融合系统模型——整夜雷达+音频模型：由上述(1)、(2)两个模型共同组成的整体模型，将雷达传感器和音频传感器分别获取的特征进行时间对准、特征融合后送入整夜雷达+音频总模型中，由11个雷达特征训练得到雷达剩余片段模型；由两个传感器同时获取的11个雷达特征和23个音频特征融合，训练得到雷达+音频片段模型；两个模型分别用于得到睡眠分期预测结果，得到雷达剩余片段{t1,t2},{t2,t3}…{tn,tn+1}和雷达+音频片段{t1,t2},{t2,t3}…{tn,tn+1}，将对应时间段的进行拼接，得到特征级融合系统的整夜睡眠的训练数据集。
步骤6、利用特征调整算法计算上述特征所占权重，根据权重筛选特征，之后进入分类器判决；如果是雷达剩余片段模型则用subspaceknn分类器，如果是雷达+音频片段模型则用baggedtrees分类器，得到睡眠分期数据；
利用特征调整算法计算上述特征所占权重，根据权重筛选特征具体为：
步骤6-1、从训练样本集d中随机抽取包含11个雷达特征和23个音频特征的数据片段，设为样本数据r，设m为样本抽样次数，然后从与样本r相同的睡眠分期类别的样本集中找出样本数据r的k个最邻近的样本hj(j＝1,2…,k)，所述训练样本集d包含11个雷达特征和23个音频特征；
步骤6-2、从每个样本数据r的不相同睡眠分期类的样本集中均找出k个距离最邻近的样本mj(c)，c表示除了样本r所属的睡眠分期类别以外的类别；
步骤6-3、更新每个特征的权重，具体计算公式如下：
上式中diff(a,r1,r2)表示样本r1与r2在特征a上的差值，其中mj(c)表示c类中的第j个最邻近的样本，diff(a,r1,r2)的计算公式如下，
所述a表示在样本集d包含11个雷达特征和23个音频特征中的任意特征。
步骤6-4、对特征进行判断，对于relieff算法计算结果为负的特征，在传感器融合中起负作用，将其删除；对于relieff算法计算结果为正的特征，在传感器融合的睡眠分类中保留这些特征。
步骤7、将两个分类器结果按照时间顺序进行拼接，得出最终睡眠分期结果。具体为：
步骤7-1、由步骤6得出的特征筛选结果对步骤5得出的模型进行进一步处理，去除由步骤6计算出权重为负的特征，保留特征权重为正的特征。如果由雷达+音频传感器获得的特征包含音频特征，则执行步骤7-2，否则执行7-3；
步骤7-2、包含雷达和音频35个特征(1个时序特征+11个雷达特征+23个音频特征)送入雷达+音频片段模型，用分类器baggedtrees得出睡眠分期结果；
步骤7-3、只包含雷达特征12个特征(1个时序特征+11个雷达特征)送入雷达片段模型，用分类器subspaceknn得出睡眠分期结果；
步骤7-4、两个分类器结果按照时间顺序进行拼接，最终得到每包特征数据对应预测的标签结果。
本发明与现有技术相比，其显著优点为：1)在非接触式的睡眠分期监测方法中患者不需要与任何仪器接触、操作简单，可长期监测，成本较低；2)本发明使用多个传感采集整夜睡眠数据，从不同的角度进行睡眠分期，分期结果更加科学；3)本发明使用的特征级融合的分类器对各个传感器的判定结果进行融合，利用多个传感器数据，提高准确率；4)特征级融合的容错能力强，也适用于异类传感器。
下面结合说明书附图对本发明做进一步描述。
附图说明
图1为本发明基于多传感器特征优化算法的睡眠分期方法流程图。
图2为本发明方法的实验示意图。
图3为本发明方法的雷达信号处理流程图。
图4为特征级融合系统模型框图。
图5为雷达+音频片段模型流程图。
图6整夜雷达+音频的特征级融合系统模型的简要流程图。
图7分期结果对比图。图(a)为标准睡眠分期图，图(b)为未经relief特征优化算法直接分期的结果图，图(c)为进过relief特征优化算法分期的结果图。
具体实施方式
结合附图，本发明的一种基于多传感器特征优化算法的睡眠分期方法，包括以下步骤:
步骤1、整个系统所需的设备有：(1)雷达传感器：采用连续波雷达，采集人体的生命体征信号(2)音频传感器：采集呼吸和鼾声信。如图2所示，测试者躺在床上，雷达设置于距离人体胸部正上方，音频设备设置于朝向人脸的方向；
步骤2、同时开启雷达传感器和音频传感器采集信号，测试者入睡。雷达照射到人体胸腔表面形成回波，测试者睡眠产生的鼾声被音频传感器记录；待多导睡眠记录仪到达监测停止时间或测试者醒来起身后，关闭雷达及音频设备。
步骤3、利用数字信号处理方法对原始回波进行处理，具体为：对雷达信号进行反正切解调，提取出体动信号、呼吸和心跳信号，使用降噪算法对原始音频进行处理；
步骤3-1本系统采用连续波雷达，是射频为2.475ghz的数字中频雷达。雷达发射电磁波信号，照射到人体后在人体表面反射形成回波信号。由于人在睡眠过程中的心脏跳动和呼吸时的胸腔运动都属于周期性的运动，雷达发射的脉冲信号经过人体胸腔的反射，反射的回波信号中携带着相应的生命体征信息，再由雷达接收机接收到回波信号。
如图3所示，首先对获得的i/q两路回波信号进行反正切解调，根据信号的特点，提取出体动信号，再根据呼吸心跳的特点和区别，选取合适的滤波器，分别滤得呼吸和心跳信号。经过信号处理，得到相应的呼吸、心跳及体动信息，从中提取出能用于判决睡眠分期的特征。
(1)呼吸信号的提取
成年人正常呼吸次数范围为每分钟16-24次，因为需要考虑所有极限情况，所以将呼吸次数的范围定为8-24次(低于8次为呼吸过缓，高于24次为呼吸过速)，呼吸信号的频率范围为0.13-0.4hz。心跳信号频率比呼吸信号频率大很多，因此设计低通频率低于0.13hz的低通滤波器，高频截止频率应设计为小于0.8hz，可以有效去除心跳信号、抑制低频噪声、直流偏移，获得较为纯净的呼吸信号。
(2)心跳信号的提取
成年人正常心跳次数范围为每分钟60-100次，超过100次在医学上被判定为“心动过速”，低于60次被判定为“心动过缓”。呼吸信号的高次谐波是获取心跳信号的主要干扰，考虑到呼吸信号的频率，心跳信号的低频截止频率需要高于呼吸信号的频率，因此，将心跳信号的频率范围定为0.83-3.3hz。设计低通滤波器的截止频率为0.8hz，高通滤波器的截止频率为4hz，采用低通滤波器和高通滤波器的组合可以获得有效的心跳信号。
(3)体动信号的提取
在整夜睡眠过程中，会伴有不同程度的体动发生，轻微的体动类似身体微抖动、手臂摆动等，大体动类似翻身、侧身等。体动信号的检测对于睡眠分期的特征参考也极为重要。根据不同体动的类型特点，可将体动检测方法分为三类：幅值门限检测法、能量谱检测法、均方误差法。
幅值门限检测法：发生体动时，雷达信号的回波幅度会突然增大，且较正常呼吸幅值大很多，体动结束时，雷达信号的回波幅度回到正常值。因此，根据一般情况下的体动信号特性，设置体动的检测门限，高于门限值则认定为发生体动，低于门限值则认定为没有发生体动。
能量谱检测法：能量谱是信号的四种频率特性之一，能量谱可以反应信号的强弱。当有体动事件发生时，能量谱明显增大。离散信号的能量谱公式如下：
均方误差法：是一种时域处理方法，是观测值x(n)与真实值
对于幅值较小的体动信号的能量谱也很小，但是均方误差可以将其检测出来，可以区分较小的体动。当没有体动事件时，均方误差大于0，可能会将幅值较大的呼吸误判为体动，而此时能量谱较小，与有体动时差别较大。所以可以采用能量谱和均方误差相结合的方法，再用幅值门限检测法作为预判的条件，可以较为准确的检测出体动信号。
步骤3-2、音频传感器可以通过检测睡眠中的呼吸和鼾声信号来诊断睡眠中的睡眠呼吸暂停事件以及实现不同睡眠阶段的分类。
呼吸和鼾声的区别在于前者是人在正常过程(清醒或入睡)中的呼吸声音，后者一般为睡眠过程中粗重的鼻息声。睡眠时的音频信号分为有声段和无声段，由于每个测试者的生理特性不同，整夜信号中可能存在大量无声段。在有声段中，存在多种类型的声音，如呼吸声、鼾声、咳嗽声、环境噪声等。因此要将呼吸信号和鼾声信号区分开来，
如图4所示，分为两个步骤：(1)语音片段的截取分割和降噪(2)信号分类，包括：呼吸检测和鼾声检测。
步骤3-2-1、将整夜音频信号以5分钟为一段划分，音频信号包括呼吸、鼾声和环境噪声，通过降噪算法处理来获取较为纯净的呼吸和鼾声信号
步骤3-2-2、利用自适应的门限法来检测呼吸事件和鼾声事件。
呼吸和鼾声信号提取：(1)首先检测出有声段，即声音事件的端点检测方法。常用的端点检测有双门限法、谱熵法、相关法、频带方差法、纵向盒快速变点检测算法等。(2)从有声段中区分鼾声信号和非鼾声信号。
步骤4、对步骤3得到的呼吸、体动、心跳以及音频信号进行特征提取；
步骤4-1、如图3所示，根据各个睡眠阶段的人体生命体征的特点，从雷达信号中提取出用来对睡眠阶段进行划分的特征共有11个。这里的特征以分钟为单位进行计算，但是以每30秒为间隔提取特征。
(1)rpm
rpm为每分钟的呼吸次数
(2)rpm_var
rpm_var为每分钟呼吸次数的方差，设每分钟的呼吸次数为r(n)，整夜信号中每分钟呼吸次数的平均值为
(3)rpm_ada
rpm_ada为每分钟呼吸信号幅度的差值积累，a_bre(n_bre)为呼吸信号幅度，n为每分钟内插值点数，每分钟呼吸信号幅度的差值积累ada的计算公式如下：
(4)rpm_move
rpm_move为呼吸信号中的体动信号特征；
(5)bpm
bpm为每分钟心跳次数；
(6)bpm_var
bpm_var为每分钟心跳次数的方差。b_bre(n_bre_var)为心跳次数，
(7)bpm_ada
bpm_ada为每分钟心跳信号幅值的差值积累。其中a_ada＇(n_bre_ada)为心跳信号幅值，n为每分钟内插值点数。
(8)bpm_move
bpm_move为心跳信号中的体动信号特征。根据心跳信号中体动信号的幅值设定门限，判定超过门限的点数来计算体动特征。每分钟取n个内插点，n取值范围为90<n<110，当前内插点的信号幅值超过门限，且当前内插点前后t秒内信号幅值超过门限的点数超过50个的，此内插点记为一次心跳信号中的体动点，统计t秒内心跳信号的体动点个数记为心跳信号的体动信号特征，t秒取值20<t<35。
(9)rem
rem是取某分钟内前后t秒的呼吸次数的差值之间计算关系的特征。在rem期，睡眠特征rem会明显增大，因此可以用来区分rem期和非rem期。
rem具体计算公式如下：
其中，
(10)deep
根据实际睡眠信号可知，人处于深度睡眠阶段时，体动信号非常小且频率较其他所有阶段都是最小的，因此根据体动信号和呼吸信号的关系，deep特征计算公式如下：
其中，
(11)sampen
由雷达回波信号x(t)采样获得雷达回波信号的一维时间序列x(n)，1≤n≤k_r重构为相空间矢量x*(ww)，1≤ww≤k-(m-1)，m为嵌入维数，定义相空间中任意两个矢量x*(w1)和x*(w2)的距离为：
求取样本熵的具体步骤为：
给定相似容限r，r＝hh*sd，hh是取值范围为0.1～0.25的常数，sd是雷达回波信号时间序列x(n)的标准差，对空间中每一个x*(w3)，进行模板匹配：
对w3求平均值，记为δm(r)：
维数m加1，重复上述步骤，得到gm+1(r)；
雷达回波信号时间序列x(n)的样本熵为：
实际中k为有限值，则样本熵估计为
步骤4-2、根据各个睡眠阶段的人体生命体征的特点，初步从音频信号中提取出用来对睡眠阶段进行划分的特征共有23个。这23个特征分为两大类，一类是睡眠呼吸相关特征，另一类是鼾声的线性和非线性特征。
a.呼吸相关特征，该类特征共13个，分别是：
1)rpm：每分钟内呼吸次数。
2)bvp(breathvarianceparameter)：前后t秒内呼吸次数差值。b_bvp(n_bvp)为当前t秒内呼吸次数，b_bvp(n_bvp-1)为前t秒内呼吸次数，b_bvp(n_bvp+1)为后t秒内呼吸次数，前后t秒内呼吸次数差值的计算公式如下：
bvp＝|b_bvp(n_bvp)-b_bvp(n_bvp-1)|
+|b_bvp(n_bvp)-b_bvp(n_bvp+1)|
其中t秒取值20<t<35。
3)bc(breathcycle)：呼吸周期，取两次呼吸事件之间的时间差均值。
4)rmse(rootmeansquareerror)：均方根误差。根据每个音频片段的采样点数，求出该片段内的均方根误差。x_rmse(n_rmse)为当前呼吸信号，
5)rpm_var：每分钟内呼吸次数的方差。xformer为前t秒内呼吸次数，xlatter为后t秒内呼吸次数，t秒取值20<t<35，计算公式如下：
rpm_var＝xformer-xlatter
6)eda(energydifferenceaccumulation)：能量差值累积。根据每个音频片段内的采样点数，求出相邻片段的能量累积差值和。x1(n)与x2(n)为相邻片段，n为音频片段内采样点数，计算公式如下：
7)cross_zr(crosszero)：过零率。每分钟内呼吸信号的过零点的次数。
8-10)formant1、formant2、formant3：3个共振峰。共振峰是指在声音信号的频谱图中能量比较集中的区域。取能量较强部分的音频信息有助于获得该片段的特征标记，有利于睡眠阶段的区分。
11-13)formant_var1、formant_var2、formant_var3：分别取(8-10)中3个共振峰的方差。
b.鼾声相关的线性和非线性特征
1)lle：lyapunov李雅普诺夫指数。混沌系统中运动对初值条件极为敏感，lyapunov指数是用来描述动力学系统中两个很近的初值产生相似的轨迹，轨迹随时间推移并按指数方式分离的这一现象过程的值，lyapunov指数给出的是动态系统沿其相空间主轴发散或收敛的平均速度。
根据求得的时间延迟
对每个点wj，计算其与最近邻点w′j在第i步前向演化后的距离
dj(i)＝|wj+i-w′j+i|＝dj(0)×eλ×i；
由下式算得最大李雅普诺夫指数,
2)时间延迟
3)嵌入维数m：嵌入维数，是相空间重构中的两个重构参数之一，相空间重构中计算李雅普诺夫指数的参数。利用时间延迟
4)apen：近似熵，是一种计算序列的复杂性和统计量化的参数。
将鼾声事件时间序列的一维时间序列x(n)＝(x1,x2,x3,...,xi,...,xk)，按顺序组成v维矢量vi＝[x(i),x(i+1),...,x(i+v-1)]，i＝1,2,...,k,其中k为鼾声事件时间序列x(n)的长度；对每一个i值计算相矢量vi和其余矢量vj的距离。
dij＝max|x(i+1)-x(j+1)|,1＝0,1,...,v-1
给定阈值r＝a3×sd，其中a3的取值范围为0.1～0.25，sd为鼾声事件时间序列x(n)的标准差。记录每个dij小于阈值r对应的i的个数，并求出与总的v维相矢量个数(k-v+1)的比值，记为
x(n)的近似熵为
apen＝φv(r)-φv+1(r)
5)n：统计嵌入维数m>＝4的个数。
6)d：时间序列的维数。
7-8)去趋势波动分析的阿尔法系数一α1和去趋势波动分析的阿尔法系数二α2：对长度为k的鼾声事件时间序列y2(n)，1≤n≤k，计算其累积差值：
其中，
将y2(n)划分为不重叠的n2个长度为l的区间，l即为时间尺度(窗长)，n2为区间(或窗口)数量；
对每一段鼾声事件时间序列采用最小二乘法拟合出局部趋势y′2(n)；
剔除y2(n)中每个区间的局部趋势，并计算新序列的均方根：
改变窗口长度l的大小，重复上述步骤，则存在幂律形式的波动：f(n)∝nα画出纵坐标为log[f(n)]，横坐标为log(n)的曲线，则曲线的斜率为时间序列的标度指数α。
取窗口长度为l1时对应的标度指数α为去趋势波动分析的阿尔法系数一α1；取窗口长度为l2时对应的标度指数α为去趋势波动分析的阿尔法系数二α2；
9)样本熵sampen：将鼾声事件时间序列y2(n)，1≤n≤k重构为m维相空间矢量
求取样本熵的具体步骤为：
给定相似容限r1，对空间中每一个
对i2求平均值，记为
维数m加1，重复上述步骤，得到bm+1(r1)；
鼾声事件时间序列的样本熵为：
实际中k为有限值，则样本熵估计为
10)香农熵h：即信息熵，定义为离散随机事件的出现概率，令每个鼾声事件时间序列y2(n)，对应的概率为p(y2＝yn)，1≤n≤k，则随机变量的信息熵为
时序特征at：若{a1,a2,……,an}为雷达或是音频信号特征组按照相对时间顺序排列的n个序号，对应于所处睡眠时间段，则at＝a/amax为睡眠分期的时序特征。
步骤5、构建特征级融合系统模型，包含雷达剩余片段模型和雷达+音频片段模型，雷达剩余片段模型是单独由雷达传感器的特征数据构建的，雷达+音频片段模型是由雷达传感器和音频传感器共同获取的特征数据构建的；
步骤5-1、由于音频传感器的不是所有时段都能采集到音频信号的，且只有在音频信号的有声段才能提取到音频特征，所以音频部分的特征是分片段的，即整夜的音频特征是呈片段的。而雷达传感器是所有时段都可以采集到的雷达信号的，即整夜的雷达特征是连续的。特征级融合系统的模型构建基本原理是当没有音频信号的时候单独采用雷达传感器的数据，当两个传感器都能获取信息时同时采用雷达和音频传感器的数据。为了准确的区分数据类型，这里将有音频数据的部分称为“音频片段”，与有音频片段的时间段对应的雷达数据称为“雷达片段”，没有音频数据的时间段对应的雷达数据称为“雷达剩余片段”。
特征级融合系统模型的构建框图如下图5所示，根据传感器的以上特性，首先分别构建两个模型：
(1)雷达剩余片段模型：单独由雷达传感器的特征数据构建的。将雷达特征数据对应有音频特征数据的时间段进行分割，分为雷达剩余片段1、2…n和雷达片段1、2…n。用雷达剩余片段1、2…n训练雷达剩余片段模型。
(2)雷达+音频片段模型：如图5所示，由雷达传感器和音频传感器共同获取的特征数据构建的。将同时有雷达和音频的特征数据融合称为雷达+音频片段，用雷达+音频片段1、2…n训练雷达+音频片段模型。再构建特征级融合系统模型——整夜雷达+音频模型：由(1)、(2)两个模型共同组成的整体模型。其中，雷达传感器提取的特征个数为11个，音频传感器提取的特征个数为23个(不含时序特征)。将雷达传感器和音频传感器分别获取的特征进行时间对准、特征融合后送入整夜雷达+音频总模型中。由11个雷达特征训练得到雷达剩余片段模型。由两个传感器同时获取的11个雷达特征和23个音频特征融合，训练得到雷达+音频片段模型。两个模型分别预测结果得到雷达剩余片段结果1、2…n和雷达+音频片段结果1、2…n，将对应时间段的结果进行拼接，得到特征级融合系统的整夜睡眠分期结果。
步骤6、为了获得更高的分类准确性，以及不同特质对分类器分类的重要程度，需要对多特征进行一个选择优化的处理，选取relieff算法作为特征优化选取的方法。利用特征调整算法计算上述特征所占权重，根据权重筛选特征,之后进入分类器判决。
步骤7、如图6所示，如果是雷达剩余片段模型则用subspaceknn分类器，雷达+音频片段模型则用baggedtrees分类器，得到睡眠分期数据；将两个分类器结果按照时间顺序进行拼接，得出最终睡眠分期结果。
在非接触式的睡眠分期监测方法中患者不需要与任何仪器接触、操作简单，可长期监测，成本较低。
下面结合实施例和对比例对本发明做进一步详细的描述。
现有的睡眠测试分期方法包括以下步骤：
步骤1、测试者胸前佩戴便携式多导睡眠记录仪somnolab2，贴好各路电极，腹部绑好测试绑腹带，左手食指佩戴压力指套。
步骤2、测试者入睡。
步骤3、待测试者清醒，专业医护人员分析数据。
这些睡眠监测产品多采用接触式监测，需要在受试者身体的多个部位贴数个电极，佩戴口鼻气流管和腹部绑带。接触式的监测在很大程度上影响了受试者的生理和心理情况，且测试环境复杂，需要专业医护人员操作，价格相对昂贵。
本发明基于多传感器特征优化算法的睡眠分期方法，包括以下步骤：
步骤1、测试者：一男性实验人员，身高177cm，年龄23岁，体重70kg，bmi指数22.34，实验之前要求实验人员尽量放松心态，不做剧烈运动，饮食规律。实验时该测试者平躺在测试床上，将雷达传感器设置于人体的正上方，用于采集雷达回波信号；将音频传感器设置于朝向人脸的方向，用于采集呼吸和鼾声信号；
步骤2、开启雷达传感器和音频传感器采集整夜数据，整夜数据中，雷达传感器数据时长为435分钟，有效的音频传感器数据时长为341分钟；
步骤3、将多导睡眠记录仪接入电脑，导出睡眠分期的结果，保存计算机端雷达及音频设备采集的信号数据，完成后续的数据处理及睡眠分期等工作。
步骤4、对步骤2采集到的雷达回波信号进行特征提取，其中每分钟呼吸次数的方差s2中的窗长n_var取5；每分钟呼吸信号幅度的差值积累ada_bre中的插值点数n取10；参数rem的t取30；k_rem取2；
步骤5、根据雷达回波信号和音频信号的特征数，构建雷达+音频片段模型和雷达剩余片段模型，利用relief算法筛选特征，将筛选后特征送入两个模型里进行识别分类；雷达+音频片段模型的无时序特征时结果相比单个雷达和单个音频传感器的整夜结果都有明显提高。雷达+音频片段模型的最高准确率为90.8％，平均准确率为80.5％；相比雷达传感器，平均提高了4.07％。相比音频传感器，平均提高了4.01％，说明融合的优势明显。
步骤6、如图7所示，对比psg睡眠分期结果和基于多传感器特征级融合的睡眠分期预测结果，得到最终分期结果的准确率。其中睡眠周期分为2、3、4、5。睡眠周期2代表深度睡眠期(deep)，3代表浅度睡眠期，4代表快速眼动期(rem)，5代表清醒期。特征级融合系统的预测结果与psg的标准分期结果的吻合度达到80.86％。从图(a)中可以看出测试者e3_3的整夜睡眠数据中出现了圈中的3个深度睡眠期(deep)，圈中的3个快速眼动期(rem)，符合睡眠周期规律。图(b)为未优化后的预测图。经过relief特征优化后，准确率进一步提升，从预测图(c)中可以看出这几个周期都能准确预测出来。