基于HRV频域分析的阻塞性睡眠呼吸暂停综合征的检测系统的制作方法

本发明属于医学数据监测技术领域，具体的说，涉及一种基于HRV频域分析的阻塞性睡眠呼吸暂停综合征的检测系统。

背景技术：

睡眠呼吸障碍疾病是一种常见的睡眠疾病，临床表现为夜间睡眠打鼾呼吸暂停和白天嗜睡。如不及时治疗会导致高血压、冠心病、糖尿病、中风、神经衰弱等多种疾病。我国目前至少4000万患有睡眠呼吸障碍疾病的患者。该病会导致患者白天无法集中注意力，记忆力理解力下降，易怒、烦躁、焦虑。

临床上把睡眠呼吸暂停分为三种：①阻塞性：呼吸暂停时，胸腹式呼吸仍存在；②中枢性：呼吸暂停时胸腹式呼吸也消失；③混合性：一次呼吸暂停时，开始出现中枢型呼吸暂停，然后出现阻塞型呼吸暂停。

多导睡眠图（PSG）是目前最常用的睡眠检测手段，通过夜间连续的呼吸，动脉血氧饱和度，脑电图，心电图等指标的检测，是国际公认的诊断睡眠呼吸暂停的金标准。多导睡眠图监测虽然是目前最常用的睡眠检测手段，但是其缺点也是很明显，它需要患者在价格昂贵的睡眠实验室中呆1~2晚，需要专业技术人员帮其佩戴多种信号检测设备，费事费力且价格高昂。睡眠实验室价格高昂且数量极少，患者需要等待很长达数月才能入住睡眠实验室。因此大部分患有睡眠呼吸障碍疾病的患者无法及时确诊。另外病人入住到陌生的环境会影响正常睡眠，会影响检测结果。对一般阻塞性睡眠呼吸暂停（OSA）病人，最好能在自然的家庭睡眠环境下进行低成本的检测。

研究人员力图通过心电信号，口鼻呼吸气流，鼾声等单一信号来取代PSG，使得信号的采集更为简单易行，检测设备价格更为便宜。之所以能通过心电信号检测阻塞型睡眠呼吸障碍，是因为呼吸阻塞时，心率过缓；呼吸恢复时，心跳加快。如此导致了心跳间隔序列(RR)成锯齿波，称之为：心率的周期性变化(cyclic variation of heart rate, CVHR)。CVHR反映到频率上是在超低频部分的0.02Hz左右功率谱的变化。通常研究人员采用模式识别的方法通过检测CVHR来识别OSA病人。

心率变异性HRV分析通常分为时域分析、频域分析和非线性分析。HRV的时域分析主要是计算如下统计指标：

①SDNN：全部正常窦性心博间期（NN）的标准差；

②SDANN：全程按五分钟分成连续的时间段，先计算每5分钟的NN间期平均值，再计算所有平均值的标准差；

③RMSSD：全程相邻NN间期之差的均方根值；

④: 全部NN间期中，相邻的NN间期之差大于50ms的心搏数。

非线性分析方法有：散点图分析法、分形维数分析法、复杂度分析法、近似熵分析法。

频域分析法根据频段可划分为超低频功率（ULF）频段（≤0.003Hz），极低频频段（VLF，0.003Hz-0.04Hz），低频功率（LF，0.04-0.15Hz），高频频段（HF，0.15-0.4Hz）。其中，副交感神经调节主要影响高频（HF）组成，低频部分受交感和副交感神经神经系统的共同影响。LF/HF是交感-迷走平衡和交感调节的度量标准。

HRV的谱分析方法可以分为参数方法与非参数方法。其中，参数方法有：自回归谱分析法(AR)：非参数方法有周期图法（periodogram）、Welch法、Lomb-Scargle 功率谱法。我们知道HRV信号是非线性信号，且采样频率非等间隔，因此在这些算法中，Lomb-Scargle算法的优势在刻画HRV信号的功率谱的优势最为明显：无需重采样，更为准确。

HRV频域分析的优点是：敏感、精确。与时域分析相比，对HRV序列异常值更加敏感。频域分析不但可以用来评价个体交感神经和迷走神经的生理变异和对心脏的影响，也可以分析病人植物神经某一成分的损伤情况。

技术实现要素：

为了克服现有技术的不足，本发明的目的在于提供一种基于HRV频域分析的阻塞性睡眠呼吸暂停综合征的检测系统。

本发明的技术方案具体介绍如下。

本发明提供一种基于HRV频域分析的阻塞性睡眠呼吸暂停综合征的检测系统，其由心电信号采集模块、心电信号预处理模块、心跳间隔信号提取模块、心跳间隔预处理模块、 网络转换模块和拓扑特征计算模块组成；其中：

心电信号采集模块，获取患者睡眠期间的ECG信号；

心电信号预处理模块：通过带通滤波器滤除噪声以及基线漂移；

心跳间隔提取模块：检测R波，从而获取HRV信号；

心跳间隔预处理模块：通过滑动平均滤波器对HRV信号滤波，滤掉异常值；

网络生成模块：将HRV信号等间隔换分成N个片段，每个片段作为网络中一个节点；计算每个HRV片段的功率谱密度，通过动态时间扭曲距离计算片段之间的0-0.04Hz超低频之间的差异性；设定阈值，如果扭曲距离矩阵D中元素，则节点i与节点j相连，否则不连，通过设定阈值将DTW矩阵变成一个无权无向网络；

拓扑特征计算模块：计算无权无向网络的拓扑特征值，并确定拓扑特征值工作的最佳阈值；其中：拓扑特征值包括平均度k、全局效率E、聚类系数C和模块性Q，其计算方式如下：网络平均度k为网络节点的平均值，即

网络节点的局部聚类系数定义为：

这里为节点i的边数；

网络的局部聚类系数C定义为：

C =

全局聚类系数定义T为：

全局效率E定义为：

这里为节点与节点之间最短路径。

本发明中，采用Lomb-Scargle算法计算HRV片段的功率谱密度。

本发明中，通过受试者工作特征曲线ROC确定无权无向网络的拓扑特征值工作的最佳阈值。

本发明的有益效果在于： 本发明的检测系统操作方法简单，检测准确度高。

附图说明

图1是健康的长度为5分钟的HRV信号图。

图2是呼吸暂停的5分钟的HRV信号图。

图3是正常HRV片段对应的PSD以及低频部分的放大图。

图4是呼吸暂停时HRV片段的PSD以及低频部的放大图。

图5 本发明检测方法流程图。

图6 欧式距离与DTW距离图示。

图7 扭曲路径图示 。

图8 (a)两个正常HRV片段的超低频波形的DTW距离；(b)两个睡眠呼吸障碍时HRV片段的超低频波形的DTW距离；(c)正常HRV片段与睡眠呼吸障碍HRV片段的超低频波形的DTW距离。

图9 原始ECG信号与去除噪声以及基线漂移后的ECG信号图示。

图10 R波检测图示。

图11 原始的RR序列与滤波后的RR序列图示。

图12 度值的ROC曲线图。

具体实施方式

下面结合附图和实施例对本发明技术方案进行详细说明。

数据库介绍：我们所使用的数据来自MIT-Physionet的Apnea-ECG Database，分为训练集和测试集。训练集和测试集中各有35个个体，记录了每个个体的长度为8小时的夜间ECG信号。ECG信号抽样频率为100Hz。这些个体年龄为43.8±10.8，体重53到135公斤。AHI的范围为0到93.5。这里我们根据AHI=5为标准来判断个体的健康与否，即AHI≥5认为是OSA病人，否则为健康个体。训练集中有22个OSA病人，13个健康个体；测试集中23个OSA病人，12个健康个体。

图5为检测方法的流程图。本发明采用的基于HRV频域分析的阻塞性睡眠呼吸暂停病人进行检测时，运用的检测系统和方法具体如下：

（1）心电信号采集模块，获取患者睡眠期间的ECG信号；

（2）心电信号预处理模块：通过0.5-40Hz带通滤波器，滤除肌电噪声、工频噪音以及基线漂移。图9展示了原始ECG信号以及处理后的ECG信号。ECG信号抽样频率从100Hz上抽样到500Hz。

（3）心跳间隔提取模块：检测ECG信号的R波，提取HRV信号。图10为一段长度为5

分钟的ECG信号上的R波检测。图1是健康的长度为5分钟的HRV信号图。图2是呼吸暂停的5分钟的HRV信号图。

（4）心跳间隔预处理模块：通过滑动滤波器，窗长为5次心跳，窗内的RR信号中位数记为mRR，如果RR>1.3*mRR或者RR<1.3*mRR,则认为该RR为异常值，用mRR替换。图11展示了原始的RR序列与滤波后的RR序列。

（5）网络生成模块分为两部分：

计算HRV片段之间的相似性和无权无向网络生成。

一、计算HRV片段之间相似性：

在这里我们采用分析HRV频域特点来判断两个HRV片段的相似性。

HRV信号切割成不重叠HRV片段，每个HRV片段长度为5分钟，总共获得70个HRV片段。通过Lomb-Scargle算法计算HRV片段的功率谱密度(PSD)。

Lomb-Scargle频谱(LS)计算公式为：

这里

通过LS算法我们获取了每个HRV片段的功率谱密度。获得的功率谱密度根据频带可以划分为超低频功率谱（0-0.04Hz），低频功率谱（0.04-0.15Hz），高频功率谱（0.15Hz-0.4Hz）。

图3为图1中的健康HRV片段的功率谱密度以及在超低频部分的放大图。

图4为图2中的睡眠呼吸障碍（OSA）发病时HRV片段的功率谱密度以及在超低频部分的放大图。

这里我们只关注超低频信号部分，这是因为OSA片段的低频波形(图4)与正常片段的低频功率谱的波形（图3）形有明显的差异。因此我们用动态时间扭曲距离（dynamic time warping, DTW）计算了每个个体的70个HRV片段之间的低频功率谱波形之间差异。

动态时间扭曲距离用于衡量两个时间序列的相似性，常用于语言识别。与传统的欧式距离相比，该方法对于刻画某些时间序列相似性上更为灵活准确（见图6）。

动态时间扭曲(DTW) 的计算方法如下：时间序列和，计算两个时间序列与之间的距离，扭曲路径W为距离矩阵d构成的连续路径。

扭曲路径需要满足三个条件：

①边界条件:，,即要求扭曲路径必须从矩阵的起始位置开始和结束位置处结束。

②连续性：要求扭曲路径每一步都是连续的。

③单调性：即时间轴上是单调的。

虽然满足以上条件的路径很多，但是存在着一条最优的路径，DTW方法就是找一条动态最优的连续路径，即

基于动态最优原理可知，从起始点(1,1)到(m,n)之间的局部最优解通过递归搜索获得，

构造累计距离矩阵D，

采用DTW距离计算HRV片段之间超低频相似性，得到一个70*70的距离矩阵D。图8展示了正常HRV片段之间的超低频部分的DTW距离(D=0.22)；两个睡眠呼吸障碍时HRV片段之间的超低频部分的DTW距离(D=0.06) ；一个正常HRV片段与一个睡眠呼吸障碍HRV片段之间的超低频部分的DTW距离(D=0.54)。可以非常直观的看到，正常的HRV片段的超低频部分极为相似，睡眠呼吸障碍HRV的超低频部分也极为相似，正常HRV片段与睡眠呼吸障碍HRV片段的超低频部分波形差别很大。

二、HRV信号的无权无向网络构建：

网络的节点为70个长度为5分钟的HRV片段。

网络的边：设定阈值，如果距离矩阵D中元素，则节点i与节点j相连，否则不连。通过设定阈值可以将DTW矩阵变成一个无权无向网络。

加权网络转换到无权无向网络的阈值的设定非常重要，影响OSA病人检测的准确度。实验表明，阈值在[0.44 0.55]范围内均能取得不错的结果。

(6) 拓扑特征计算模块：计算网络拓扑特征值，并确定拓扑特征值工作的最佳阈值。

为了确定哪一个网络指标能够有效的区分OSA病人，我们测试了如下的网络拓扑

特征值：平均度k、全局效率E、聚类系数C、模块性Q。首先我们简单介绍一下这几个指标的计算方法。

网络平均度k为网络节点的平均值，即

网络节点的局部聚类系数定义为

这里为节点i的边数。

网络的局部聚类系数C为每个节点局部聚类系数Ci的平均值，定义为：

C =

全局聚类系数定义T为

全局效率E定义为

这里为节点与节点之间最短路径。

首先测试单个指标对OSA区分效果。为了判断单个拓扑特征值在何时能取得最佳的区分效果，我们采用受试者工作特征曲线（ROC）来判断阈值。通过ROC方法（受试者工作特征曲线）为每一个拓扑特征值寻求最好的区分阈值使得假阳性率与假阴性率值最低。图12为度值的ROC曲线，当网络阈值时，网络平均度阈值时，网络平均度能达到最佳的区分效果，即度值大于34认定为OSA病人，而度值小于34认定为健康个体。在训练集、测试集中都能达到完美的区分效果(敏感性Se为100%，特异性Sp为100%，准确性Ac为100%)。具体结果见表1的第一行。另外，全局效率E，全局聚类系数C，局部聚类系数T，模块性Q的区分阈值，以及在训练集测试集中区分的敏感性、特异性、准确性也在表1中一一列出。

表1是 单个拓扑特征值在训练集与测试集中的区分效果。

表1

通过ROC方法获得全局效率E最佳区分阈值，在训练集中的敏感性Se为100%，特异性Sp为100%，准确性Ac为100%，在测试集中的敏感性Se为100%，特异性Sp为100%，准确性Ac为100%。局部聚类系数=0.82，在训练集中敏感性Se为100%，特异性Sp为100%，准确性Ac为97.14%，在测试集中敏感性Se为100%，特异性Sp为100%，准确性Ac为100%。全局聚类系数，在训练集中敏感性Se为95.45%，特异性Sp为92.31%，准确性Ac为94.29%，在测试中敏感性Se为95.65%，特异性Sp为100%，准确性Ac为97.14%。模块性Q在时，取得最佳的区分效果，在训练集中，敏感性Se为90.91%，特异性Sp为92.31%，准确性Ac为91.43%；在测试集中，敏感性Se为100%，特异性Sp为100%，准确性Ac为100%。