张大可， 马隽， 王立英， 王钢*， 蔡靖， 孙玉冰

(1.北华大学电气与信息工程学院， 吉林 132021； 2. 北华大学基础医学院， 吉林 132021； 3. 清华大学附属垂杨柳医院， 北京 100020； 4.吉林大学仪器科学与电气工程学院， 长春 130061)

睡眠呼吸暂停(sleep apnea，SA)是指在连续7 h睡眠中发生30次以上的呼吸暂停，每次气流中止10 s以上(含10 s)，或平均每小时低通气次数(呼吸紊乱指数)超过5次，而引起慢性低氧血症及高碳酸血症的临床综合征，是一种睡眠时呼吸停止的睡眠障碍[1-2]，主要表现为睡眠过程中反复发生上气道塌陷、阻塞，睡眠打鼾和白天嗜睡，常见于中年以上肥胖男性，对患者的心血管系统和内分泌系统等有着严重的损害，增大肿瘤患病的可能性[3]。

目前关于睡眠呼吸暂停的诊断主要是利用多导睡眠仪(polysomnography，PSG)[4]对患者的脑电、心电(electrocardiogram, ECG)、血氧饱和度或眼电等生理信号进行监测，然而PSG价格昂贵，且佩戴电极繁杂，佩戴过程烦琐，影响患者正常睡眠[5]，监测过程费时费力，往往需要数位医生共同分析监测结果，无法给予患者及时有效的诊断反馈。

Guilleminault等[6]最早利用单导ECG信号对睡眠呼吸暂停实现了诊断。近年来，研究者们发现，睡眠呼吸暂停患者在呼吸暂停发作时，会产生明显的心率周期性变化现象。梁九兴等[7]利用心率变异性(heart rate variability, HRV)的时、频域特征构建概率神经网络(probabilistic neural network，PNN)模型，用于检测SA，使用10折交叉验证实现准确率达75.97%。Cafer等[8]通过心电信号估算心电伴生呼吸(derived respiration，EDR)信号，进而利用EDR信号估算受试者瞬时心率序列，对3 min信号对应的瞬时心率序列进行小波分解，再利用非线性自回归型人工神经网络分类器构建分类模型，准确率可达78.3%。范文兵等[9]通过多尺度熵证明ECG信号多尺度熵指数能够区分SA。葛靖等[10]采用长短时记忆模型，通过心电的间接信号提取特征，识别睡眠呼吸暂停，准确率可达87.4%。然而，利用上述方法进行睡眠呼吸暂停检验或分类的准确率不高，导致误检的可能性增加。

现采用粒子群优化-支持向量机(particle swarm optimization-support vector machine，PSO-SVM)方法[11]，分析处理受试者的心电信号，实现对睡眠呼吸暂停的准确检测，以期不仅减轻医务工作者的负担，也及时分析患者的检测结果并得出结论，提高工作效率。

1 方法

1.1 总体框架

提出的睡眠呼吸暂停检测方法的主要计算流程如图1所示，首先读取数据文件，根据注释文件对睡眠呼吸暂停状态进行分类，再根据期望数据要求对心电信号进行预处理，去除噪声和因其他原因产生的杂乱信号，进而利用峰值检波提取心电信号的R波，并计算心率和心率变异性(heart rate variablity，HRV)，对上述信号进行特征提取，再利用方差阈值法对特征进行筛选，最后利用粒子群算法优化支持向量机，从而实现对睡眠呼吸暂停的检测。

图1 总体框图Fig.1 Overall block diagram

1.2 数据预处理

数据预处理主要包含两方面：对心电信号进行切割分段；对心电信号进行心率变异性的计算。

首先对ECG按照60 s时间窗分段，对应睡眠呼吸暂停或正常睡眠标签；在分段后，利用峰值检测寻找R波位置，并记录R波产生时间。将R波产生时间相减，得到受试者RR间期信号，进而得到受试者的心率变化，再将相减后的R波时间做差分，即可得到心电信号的心率变异性HRV序列[12]。

1.3 特征提取

根据心电信号和心率变异性的定义和特点，以及睡眠呼吸暂停产生过程中患者心率波动，对心率变异性的时域、频域、时频域和非线性特征[13]进行了提取，时域特征主要包含均值(MEAN)、总体标准差(SDNN)、均值标准差(SDANN)和差值均方的平方(r-MSSSD)等[14]，频域特征主要包含HRV频谱特征、功率谱密度(PSD)、信号总功率(TP，f≤0.4 Hz)、甚低频段功率(VLF，0.003 3 Hz≤f≤0.04 Hz)、低频段功率(LF，0.04 Hz≤f≤0.15 Hz)、高频段功率(HF，0.15 Hz≤f≤0.4 Hz)和低频段与高频段功率比(LF/HF)[15]，时频域特征主要包含瞬时能量和瞬时中位频率[16]，非线性特征主要包含LZ复杂度、近似熵(ApEn)、样本熵(SampEn)、信息熵(InEn)、边际谱熵(MSEn)和能量谱熵(ESEn)[17]。部分特征参数的计算公式如下所示。

(1)

(2)

(3)

(4)

(5)

(6)

(7)

式(7)中：pm为第m个边际谱在整个边际谱中所占的百分比，即

(8)

(9)

(10)

式中：hx为边际谱；H(m,t)为Hilbert谱；Pk为第k个能量谱在整个能量谱中所占的百分比，即

(11)

(12)

式中：ex(k)为能量谱；H2(k,t)为能量谱密度。

1.4 特征选择

在进行特征提取后，由于上述特征间存在冗余，且存在与睡眠呼吸暂停相关性小的特征，因此采用方差阈值法，计算各个特征的方差，删除方差低于预设阈值的特征。

方差阈值法(variance thresholding, VT)是一种利用方差实现特征提取的方法，通过程序设置允许通过的方差最小值，即可对方差小于该值的特征进行去除，从而筛选与检测对象相关性大的特征，提高检测的准确率。方差表示数据分布的可变性，公式为

(13)

若一组数据的方差为0，则表示该组数据在预测过程中无意义，只能增加模型的复杂性，而不会增加模型的预测能力，方差阈值法就是将无法提高模型预测能力的特征筛除的过程。

然而，由于不同特征的分布不同，导致计算所得方差的尺度也不尽相同，无法直接利用方差判别其有效性，因此首先对特征进行归一化处理，即每一个特征feati除以特征的均值，得到归一化特征z_feat，再求取其方差，实现特征选择。

(14)

(15)

1.5 SA检测模型建立

由于ECG信号的非线性特征，以及SA检测标签仅分为“呼吸正常”和“呼吸紊乱”两类，在特征选择后，利用支持向量机(SVM)对睡眠呼吸暂停进行训练和检测。SVM是一种基于间隔最大化原则实现二分类的机器学习方法[18-19]，其目标函数和约束条件可表示为

(16)

式(16)中：L(a)为拉格朗日函数；ai为拉格朗日系数；xi为数据向量；yi为xi对应的数据类别；K(xixj)为核函数。基于径向基核函数的非线性SVM具有较好的性能，因此拟选择径向基核函数；C为惩罚因子，C越大则惩罚越重，且泛化能力会随之降低；高斯核超函数g表示对径向范围的控制。

由于惩罚因子C和高斯核超函数g会严重影响SVM的分类能力，采用粒子群算法(PSO)实现参数寻优，以提高SVM分类器的泛化性能。粒子群算法是一种基于迭代的优化算法，结构简单，具有全局寻优能力，能够以高概率和高收敛速度找到最优解。PSO算法中根据式(17)和式(18)更新粒子的速度和位置[19-20]。

(17)

(18)

利用粒子群算法优化支持向量机的程序流程图如图2所示。

图2 PSO-SVM程序流程图Fig.2 PSO-SVM program flow chart

2 结果与分析

2.1 数据源

所采用的数据来源于2000年心脏病学计算会议和PhysioNet联合进行的一项竞赛中首次公开的Physionet Apnea-ECG数据库[21]。数据库中共包含70组采样频率为100 Hz的心电图，总持续时间为401～578 min，正常呼吸时间在11～535 min变化，呼吸紊乱时间在0～534 min变化。A类为呼吸暂停组，即患者组，定义为呼吸紊乱超过100 min的受试者，共包含16名29～63岁受试者的40条记录，其中15名男性，1名女性；C类为对照组，即健康组，定义为呼吸紊乱少于5 min的受试者，共包含11名年龄为27～42岁的受试者的20条记录，其中6名男性，5名女性。

在数据库提供的测试集文件中，列举了20条患者组和10条健康组的ECG信号和标签，ECG信号存储于.dat文件，ECG信号的标签存储于.apn文件，即记录每分钟是否存在睡眠呼吸暂停事件，若存在，则标记为A；否则标记为N。

2.2 数据预处理

部分原始心电信号如图3所示，由于数据库定义每分钟心电信号存在一个标签，对心电信号按照60 s的时间窗分段，每段数据包含60×100个数据点，对应睡眠呼吸暂停(A)或正常睡眠(N)标签。取每位受试者测试时间前423 min即423×60×100=253 800个采样点的测试数据，按60 s分为423段，分别对应每分钟的SA标签，取14位受试者共计5 922组样本数据。

图3 原始心电信号Fig.3 Original ECG signal

图4为利用峰值检波检测R波位置结果；图5为利用RR间期计算受试者1 min心率变化；图6为利用RR间期计算受试者1 min心率变异性波动；图7为心率变异性频谱图。

图4 心电信号R波定位Fig.4 R wave location of ECG

图5 受试者1 min内心率波动Fig.5 1 min heart rate fluctuation of subjects

图6 受试者1 min内心率变异性波动Fig.6 1 min HRV fluctuation of subjects

图7 HRV频谱图Fig.7 HRV spectrum

2.3 特征提取

图8为心电信号RR间期的均值、标准差、均值标准差和标准差均值的变化情况；图9为某受试者HRV功率谱密度；图10为HRV信号总功率、低频段功率、高频段功率和功率比的变化情况；图11为HRV边际谱熵和能量谱熵的变化情况。可以看出，受试者处于不同睡眠呼吸阶段，其心电信号各项特征存在明显的波动变化，可用于进行睡眠呼吸暂停的标记。

图8 时域特征Fig.8 Time domain features

图9 HRV功率谱密度对数Fig.9 Logarithm of HRV power spectral density

图10 频域特征Fig.10 Frequency domain characteristics

图11 非线性特征Fig.11 Nonlinear characteristics

2.4 特征选择

通过方差阈值法，按照各特征的归一化方差排序，共筛选了10个特征放入检测算法中，其中包含4个时域特征(均值、标准差、均值标准差、差值均方的平方)、3个频域特征(信号总功率、低频段功率、高频段功率)、1个时频域特征(瞬时中位频率)和2个非线性特征(边际谱熵、能量谱熵)。

2.5 模型评价

利用粒子群算法对支持向量机的参数进行选取，进而实现睡眠呼吸暂停的检测，采用准确率(Acc)衡量SVM检测表现[22-23]，公式为

(18)

式(18)中：TP表示真正结果数；TN表示真负结果数；FP表示假正结果数；FN表示假负结果数。

图12为粒子群算法适应度曲线，经过300次迭代，针对SVM的参数优化的最优解的准确率可达97.89%，优化后的惩罚因子C=100，高斯核超参数g=0.01。

图12 PSO适应度曲线Fig.12 PSO fitness curve

图13为利用PSO-SVM实现睡眠呼吸暂停的部分结果与其真实值的对比情况，测试结果良好。

图13 部分训练结果对比Fig.13 Comparison of some training results

PSO-SVM与SVM相比，实现了SVM最优参数的选取，减小了算法的时间复杂度[24]，假设支持向量个数为N，输入向量维度为dl，训练样本点个数为1，则基于径向基核函数SVM的计算复杂度为O(N3+N2l+Ndll+O(dl)N)；假设存在M个N维粒子，进行T次迭代，则PSO-SVM的时间复杂度为O(TMN+dll2+O(dl)N)。图14绘制了利用PSO-SVM和SVM检测SA的接受者操作特征(receiver operating characteristics，ROC)曲线，展示了PSO-SVM的良好准确率，曲线下的面积(area under curve，AUC)为96.94%。表1中列举了不同方式SVM算法的执行时间，以本文训练数据规模为基础，PSO-SVM可以满足对准确率和执行效率的要求。

图14 PSO-SVM与SVM的ROC曲线对比图Fig.14 ROC curve comparison between PSO-SVM and SVM

表1为单独利用SVM不同核函数检测SA和利用PSO-SVM检测SA的结果参数，可以看出，利用粒子群算法对支持向量机进行优化，提高了支持向量机检测的准确率，可达94.0%。

表1 训练结果对比Table 1 Comparison of training results

3 结论

利用粒子群算法优化支持向量机对睡眠呼吸暂停进行检测，通过单一的心电信号，降低了多种信号采集处理的复杂性和困难程度；利用方差阈值法实现ECG信号的特征选择，减少检测模型的输入，提高检测效率；通过PSO-SVM对睡眠呼吸暂停建立预测模型，使检测准确率提升至94.0%，优于单一的SVM检测模型。实验结果表明，单一的心电信号通过求取心率变异性及其他特征，可以用于睡眠呼吸暂停的检测，且分类性能良好。

不足之处在于数据来源单一且测量时间较早，与当前国内患者可能存在一定的信号差异；此外，由于数据标签限制，仅实现了睡眠呼吸暂停的二分类检测，无法根据模型结果判断其严重程度。在未来的工作中，可以进一步扩展数据来源，深度优化计算模型，提高泛化能力和检测范围。