基于心率变异性与机器学习的睡眠呼吸事件分类*

梁九兴，张湘民，黄少雄，曾令紫，罗语溪

（1.中山大学工学院，广东广州510006；2.中山大学附属第六医院，广东广州510655）

睡眠呼吸暂停低通气综合征（SAHS）是由于在睡眠期间反复出现的呼吸暂停或低通气状态，其中每次呼吸暂停或低通气持续时间至少达到10 s，且暂停或低通气指数（AHI）达到 5次/h以上[1-2]。SAHS是睡眠呼吸障碍中最为常见的疾病之一。根据大量国内外相关临床病学调查研究，在美国约有4%的男性及2%的女性患有SAHS[3]。而在我国相关的调查研究中，亦普遍认为人群中约有4%患有该疾病[4-5]，由此可知SAHS普遍存在于广大人群中。除此以外，SAHS具有巨大的潜在危害性。这是由于反复的呼吸暂停或者低通气状态，使机体长期处于缺氧及二氧化碳的潴留可能引致诸如低氧血症、高碳酸症[6]、肺动脉高压[7]、脑卒中[8]、心血管疾病[9-10]、心衰[11]等临床症状。

虽然SAHS疾病存在巨大的潜在危害性且普遍存在于大众中，但中重度SAHS患者中有93%的女性及82%的男性却没有得到及时的诊断[12-13]。SAHS的发病率之高与有条件被诊断患者的比例之低是临床和流行病学调查工作的主要矛盾。随着国民对SAHS的认识程度和国民收入水平的不断提高，越来越多的疑似患者将会去医院就诊并前往睡眠实验室进行监测，上述结果必然导致通过多导睡眠（polysomnography，PSG）检查疑似 SAHS患者人数的上升。尽管PSG是进行各种睡眠障碍疾病诊断及治疗评价的金标准[1，14]，但是，PSG监测具有生理负荷高、设备检查环境要求高、检查和分析技术复杂、人力资源消耗大和费用相对昂贵的特点。开发低生理负荷、便携式的监测设备及相应的诊断算法已经成为大势所趋。在此过程中，针对异常睡眠呼吸事件的识别算法是其技术的关键因素。

通过提取心率变异性的时域、频域及非线性域特征参数，并采用机器学习对训练样本进行模型训练，可以对测试样本进行预测分类。一些研究人员已把心率变异性与各机器学习算法应用于多种预测分类场合，如睡眠分期[15]和与心血管疾病相关的病征[16-17]的预测问题，但尚未涉及对睡眠呼吸事件的预测识别研究。心率变异性与呼吸节律存在着紧密的联系，因此有可能通过心率变异性各特征参数实现对异常睡眠呼吸事件的预测分类，并结合短时耗的分类算法，实现对SAHS病患低生理负荷的筛查。

1 睡眠呼吸事件分类模型的建立

本研究的方法包括数据提取、心电数据预处理、心率变异性特征值的提取、机器学习模型的训练及分类输出等步骤，具体流程详见图1。

图1 基于心率变异性的睡眠呼吸事件分类流程图Fig.1 The flow chart of the classification method of SAHSbased on HRV

1.1 实验数据的采集

本研究共有19例SAHS患者的整晚多导睡眠监测（PSG）实验数据。在进行PSG监测之前医生需提前告知在参与本监测实验前6 h内不能服用任何药物和饮用含有酒精的饮料以及含有可卡因等会导致精神振奋的食品。并且，参与本实验的对象均主诉无其他心血管类病史。本实验均在中山大学附属第六医院睡眠呼吸障碍诊疗中心进行。19例SAHS患者的人口统计学和睡眠监测基本信息，如表1所示。

表1 人口统计学和睡眠监测基本信息Table 1 Demographic and sleep monitoring basic data

为保障本实验具有足够的正常呼吸事件、异常呼吸事件样本参与机器学习，本实验的入组和退出准则如下：1）AHI控制在10～60次/h；2）整晚有效监测的总睡眠时间长度不应小于300 min；3）剔除心电信号和口鼻呼吸气流信号伪差明显的患者；4）剔除患有失眠、癫痫、心血管等病史的患者。

1.2 特征向量的提取

从正常睡眠呼吸状态及异常睡眠呼吸事件状态下的多导睡眠监测数据中提取对应的心电信号。提取心电信号时以60 s为一个数据截取窗长度。为了提取心率变异性各特征参数作为分类学习的特征向量，首先需对截取的心电信号进行一维墨西哥帽小波变换用以去除基线漂移[18]。其次，采用差分阈值算法检测心电信号的QRS复合波后通过区域最大值检测R波峰，再通过最小心搏时间间隔去除异搏[19]。然后，提取RR间期序列。在心率变异性的频域特征提取过程中，还需对RR间期序列进行3次样条插值及重采样处理[20-21]，并采用自回归模型的Burg算法对插值处理后的RR间期序列进行现代功率谱密度估计 。最后，从功率谱密度估计中提取各频域的特征参数。

本实验中，输入分类器进行学习训练和分类的特征向量共有12个特征值，如表2所示。由于个体之间各心率变异性参数之间的差异较大，为了能够消除每位实验参与者之间的个体性差异以及方便分类器的特征输入，需对每一位实验参与者的特征参数按照式（1）进行归一化处理：

x指代输入心率变异性的特征序列。归一化处理之后，对每一组的12个特征参数根据PSG标准进行特征向量的标定。

表2 心率变异性的特征向量Table 2 Eigenvector of heart rate variability

1.3 概率神经网络分类算法

概率神经网络（PNN）是一种有导师监督学习的机器学习方法，最早于1988年由Specht[23-24]提出。它是基于贝叶斯概率决策准则而构造的一种概率密度分类估计及其并行处理的前馈型神经网络。因此，概率神经网络既具有神经网络所拥有的一般属性，亦具有较高的分类准确精度，很好的泛化能力以及由并行处理带来的快速学习能力[25]。正因为概率神经网络具有上述优势，该机器学习模式被广泛地应用于各类分类预测问题。概率神经网络的算法架构如图2所示[26]。

PNN的架构分别由输入层，中间模式层，中间求和层以及输出层组成。PNN分类算法包含以下3个基本步骤：内核计算，类条件概率计算及最大类条件概率的选择。设某一训练样本集T＝｛xCi，j｝，对某一分类问题训练样本集中下标j指代的是特征样本归属Cj的类别。如果样本集T中共有NC组类别，则该样本可分为NC个子集，此时Ci＝｛xCi，1，xCi，2，…，xCi，NC｝。其具体计算流程可由图3进行归纳[26]。

图2 PNN算法架构Fig.2 The architecture of PNN algorithm

图3 PNN算法流程Fig.3 Flow diagram of PNN algorithm

1.4 预测结果评价

采用准确率、灵敏度、特异性、受试者工作特性曲线及曲线下面积（AUC）等评价指标进行PNN分类算识别性能的评估。

TP指代的是正确地分类出异常睡眠呼吸事件的例数，TN指代的是正确地分类出正常睡眠呼吸状态的例数，FP指代的是错误地把正常睡眠呼吸状态分类成异常睡眠呼吸事件的例数，FN指代的是错误地把异常睡眠呼吸事件分类成正常睡眠呼吸的例数[27]。

受试者工作特性曲线（receiver operating characteristic curve，ROC）[28-29]以测试灵敏度为纵坐标、1-特异性或假阳性率（false positive rate，FPR）为横坐标曲线图，用以表示测试的真阳性率（true positive rate，TPR）与其假阳性率（FPR）之间的关系。它是评估分类器整体性能的一种方法[30]。实际应用中，通过计算 ROC曲线下面积（AUC）进行量化评估。AUC以0～1之间的数值表征从测试样本中随机选择某一正例比起选择某一负例的概率要高[31]。因此，AUC数值越大，则对应分类器的整体性能越好[32]。

2 实验结果与分析

共有5 910组心率变异性特征参数参与机器学习分类实验。其中，标记为正常睡眠呼吸状态有2 816组，而标记为异常睡眠呼吸事件有3 094组。异常睡眠呼吸事件与正常睡眠呼吸状态数据总体上相对平衡。在机器学习中，为了体现本算法的泛化性能，实验随机抽取一半的数据集作为训练样本，剩余一半作为测试样本。

2.1 概率神经网络参数训练

在概率神经网络（PNN）对有无异常睡眠呼吸事件的预测分类研究中，需要对径向基函数的扩展速度进行选取。本实验选取了在0.02～0.4之间，步进扩展速度为0.02，对每一种扩展速度随机进行10次计算分类，其运行结果如图4所示。结果表明，径向基函数的扩展速度为0.16时可达到最优准确率。因此，PNN算法中径向基函数的扩展速度取为0.16。

2.2 概率神经网络分类结果

以0.16作为径向基函数扩展速度的最优值进行预测分类。在该参数下随机计算10次，得到的预测分类性能，如表3所示。

图4 径向基函数扩展速度预测准确率的影响Fig.4 The effect of the spread of the radial basis function on prediction accuracy

表3 PNN算法的分类预测结果Table 3 Classification of PNN algorithm

在表3中，PNN算法的十次分类预测的平均准确率为75.97%，其标准方差为0.63%；平均灵敏度为82.51%，其标准方差为1.33%；平均特异度为76.22%，其标准方差为1.70%；平均AUC为0.793 6，其标准方差为0.005 0。

将PNN分类算法与广义回归神经网络（GRNN）、极限学习机（ELM）算法进行对比。3种分类算法采取同样的数据、同样的训练和测试数据长度，运行10次后所得的性能评估结果，如表4所示。

表4 PNN、GRNN及ELM算法的预测性能Table 4 Prediction performance of PNN，GRNN and ELM algorithms

用于睡眠呼吸事件识别的算法应有较好的准确性与较少的计算时耗。本实验采用的计算机性能参数为：中央处理器Interl Core i5-3470，处理器工作频率3.2 GHz，内存8.0 GB。我们对比了PNN、GRNN及ELM3种机器学习算法所消耗的计算时间。3种分类器进行10次分类预测消耗的计算时间，如图5所示。

图5 3种分类算法的测试耗时Fig.5 The time-consuming test of three classification algorithms

PNN、GRNN及ELM分类器进行十次计算所耗费的平均时间分别为（0.63±0.06）、（0.97±0.16）及（3.49±0.21）s。由此可知，3种分类器的耗时长短顺序依次为：PNN＞GRNN ＞ELM。综上所述，3种机器学习算法的准确率从高到低的序位为ELM＞GRNN＞PNN、灵敏度从高到低为PNN＞GRNN＞ELM、特异度从高到低为PNN＞GRNN＞ELM、AUC从高到低为PNN＞GRNN＞ELM、消耗的计算时间从短到长为PNN＞GRNN＞ELM。因此，若仅准确率为目标则应选择ELM作为最优算法。若综合灵敏度、特异性、AUC及预测时间耗时则PNN分类算法最优。

3 结 论

本研究应用PNN分类算法判别了有无异常睡眠呼吸事件，并对该分类算法的性能进行了多维度的评估，并将之与GRNN、ELM分类算法进行了比较。结果表明：1）PNN算法具有较高的准确性；2）与GRNN和ELM算法相比，PNN算法的灵敏度、特异度及AUC面积最优，并且其分类耗时最短。因此，PNN算法在对异常睡眠呼吸事件的分类识别上具有较好的性能，能为低生理负荷的SAHS筛查提供了一种新的方法，具有一定的临床实际应用意义。

[1]BERRY R B，BROOKSR，GAMALDO C E，et al.The AASM manual for the scoring of sleep and associated events：rules， terminology and technical specifications[M].2nd ed.American：Academy of Sleep Medicine，2012：1-59.

[2]RANKLIN K A，LINDBERG E.Obstructive sleep apnea is a common disorder in the population-a review on the epidemiology of sleep apnea[J].JThorac Dis，2015，7（8）：1311-1322.

[3]Al Lawat N M，PATEL SR，AYASN T.Epidemiology，risk factors，and consequences of obstructive sleep apnea and short sleep duration[J].Prog Cardiovasc Dis，2009，51（4）：285-293.

[4]上海市医学会呼吸病学分会睡眠呼吸疾病学组.上海市30岁以上人群阻塞性睡眠呼吸暂停低通气综合征流行病学调查[J].中华结核和呼吸杂志，2003，26（5）：268-265.Sleep Respiratory Disorder Study Group.Respiratory Disease Branch，Shanghai Medical Association.Prevalence of obstructive sleep apnea-hypopnea syndrome in Chinese adults aged over 30 yr in Shanghai[J].Chin J Tuberc Respir Dis，2003，26（5）：268-265.

[5]胡庆磊，杜翠萍，杨扬，等.上海市普陀区20岁以上人群阻塞性睡眠呼吸暂停低通气综合征流行病学调查[J].中国眼耳鼻喉科杂志，2017，17（1）：49-54.HU Qinglei，DU Cuiping，YANG Yang，et al.Prevalence of obstructive sleep apnea hypopnea syndrome in Chinese adults aged over 20 years in Putuo district of Shanghai[J].China JOphthalmol and Otorhinolaryngol，2017，17（1）：49-54.

[6]STANSBURY R C，STROLLO P J.Clinical manifestations of sleep apnea[J].J Thorac Dis，2015，7（9）：298-310.

[7]FLORASJS.Hypertension and sleep apnea[J].Can J Cardiol，2015，31（7）：889-897.

[8]LYONSO D，RYAN CM.Sleep apnea and stroke[J].Can JCardiol，2015，31（7）：918-927.

[9]BADRAN M，YASSIN BA，FOX N，et al.Epidemiology of sleep disturbances and cardiovascular consequences[J].Can J Cardiol，2015，31（7）：873-879.

[10]PUNGINATHN D.Obstructive sleep apnea and cardiovascular risk[J].Therapeutics＆Clinical Risk Management，2015，3（6）：1105-1111.

[11]LYONS O D，BRADLEY T D.Heart failure and sleep apnea[J].Can J Cardiol，2015，31（7）：898-908.

[12]李利，麦慧娟，张素.阻塞性睡眠呼吸暂停低通气综合征患者就医状况调查[J].中华护理杂志，2014，49（2）：197-201.LI Li，MAI Huijuan，ZHANG Su.Station of hospital presentation in patients with obstructive sleep apnea hypopnea syndrome[J].Chin J Nurs，2014，49（2）：197-201.

[13]YOUNGT，EVANSL，FINN L，et al.Estimation of the clinically diagnosed proportion of Sleep Apnea Syndrome in middle-aged men and women[J].Sleep，1997，20（9）：705-706.

[14]BERRY R B，BUDHIRAJA R，GOTTLIEB D J，et al.Rules for scoring respiratory events in sleep：update of the 2007 AASM Manual for the Scoring of Sleep and Associated Events.Deliberations of the Sleep Apnea Definitions Task Force of the American Academy of Sleep Medicine[J].JClin Sleep Med，2012，8（5）：597-619.

[15]WILLEMEN T，VAN DEUN D，VERHAERT V，et al.An evaluation of cardiorespiratory and movement features with respect to sleep-stage classification[J].IEEE J Biomed Health，2014，18（2）：661-669.

[16]HOUSHYARIFAR V，AMIRANI M C.An approach to predict Sudden Cardiac Death（SCD）using time domain and bispectrum features from HRV signal[J].Bio-Med Mater Eng，2016，27（2/3）：275-285.

[17]PATIDAR S，PACHORI R B，ACHARYA U R.Automated diagnosis of coronary artery disease using tunable-Q wavelet transform applied on heart rate signals[J].Knowl-Based Syst，2015，82：1-10.

[18]HUANG G，NEWCHURCH M J，KUANG S，et al.Definition and determination of ozone laminae using Continuous Wavelet Transform（CWT）analysis[J].Atmospheric Environment，2015，104：125-131.

[19]YEH Y C，WANG W J.QRS complexes detection for ECG signal：the difference operation method[J].Comput Methods Programs Biomed，2008，91（3）：245-254.

[20]JANG D G，HAHN M，JANG JK，et al.A comparison of interpolation techniques for RR interval fitting in AR spectrum estimation[C]//IEEE Biomedical Circuits and Systems Conference（BioCAS），2012：13264332.

[21]SINGH D，VINOD K，SAXENA S C.Sampling frequency of the RR interval time series for spectral analysis of heart rate variability[J].J Med Eng Technol，2004，28（6）：263-272.

[22]PICHON A，ROULAUD M，ANTOINE-JONVILLE S，et al.Spectral analysis of heart rate variability：interchangeability between autoregressive analysis and fast Fouriertransform[J].J Electrocardiol，2006，39（1）：31-37.

[23]SPECHT D F.Probabilistic neural networks and the polynomial adaline as complementary techniques for classification[C]//IEEE International Conference on Neural Networks.1990：111-121.

[24]SPECHT D F.Probabilistic neural networks for classification，mapping，or associative memory[C]//IEEE International Conference on Neural Networks，1988：I525-I532.

[25]郁磊，史峰，王辉，等.MATLAB智能算法30个案例分析[M].北京：北京航空航天大学出版社，2011.

[26]ZEINALI Y，STORY B A.Competitive probabilistic neural network[J].Integrated Computer-Aided Engineering，2017，24（2）：105-118.

[27]FAWCETT T.An introduction to ROC analysis[J].Pattern Recognition Letters，2006，27（8）：861-874.

[28]BRADLEY A P.The use of the area under the ROC curve in the evaluation of machine learning algorithms[J].Pattern Recognition，1997，30（7）：1145-1159.

[29]GREINER M，SMITH R D，PFEIFFER D.Principles and practical application of the receiver-operating characteristic analysis for diagnostic tests[J].Preventive Veterinary Medicine，2000，45（1/2）：23-41.

[30]PARK S H，GOO J M，JO C H.Receiver Operating Characteristic（ROC）curve practical review for radiologists[J].Korean J Radiol，2004，5（1）：11-18.

[31]WALTER S D.The partial area under the summary ROC curve[J].Stat Med，2005，24（13）：2025-2040.

[32]HUANG J，LING C X.Using AUC and accuracy in evaluating learning algorithms[J].IEEE Transactions on Knowledge and Data Engineering，2005，17（3）：299-310.