CN115630290B

CN115630290B - 基于同步挤压变换的心肺耦合特征提取方法与系统

Info

Publication number: CN115630290B
Application number: CN202211465021.1A
Authority: CN
Inventors: 叶建宏; 王依宁; 史文彬
Original assignee: Beijing Institute of Technology BIT
Current assignee: Beijing Institute of Technology BIT
Priority date: 2022-11-17
Filing date: 2022-11-17
Publication date: 2023-04-07
Anticipated expiration: 2042-11-17
Also published as: CN115630290A

Abstract

本发明公开了基于同步挤压变换的心肺耦合特征提取方法与系统，能够突破传统心肺耦合图谱关键节律模糊的技术难点，实现对睡眠呼吸暂停事件的检测。具体方案为：从单导联心电信号中提取心跳间期R‑R时间序列和替代呼吸EDR时间序列。对心跳间期R‑R时间序列和替代呼吸EDR时间序列分别进行同步挤压变换，分别得到R‑R时间序列的同步挤压系数和EDR时间序列的同步挤压系数。利用R‑R时间序列的同步挤压系数和EDR时间序列的同步挤压系数计算心肺耦合CPC值，绘制睡眠期间心肺耦合图谱。基于睡眠期间心肺耦合图谱提取多类心肺耦合特征，利用机器学习分类的方法对所提取的多类心肺耦合特征进行学习分类，以识别正常呼吸事件和睡眠呼吸暂停事件。

Description

基于同步挤压变换的心肺耦合特征提取方法与系统

技术领域

本发明涉及睡眠健康监测和信息技术交叉技术领域，具体涉及一种基于同步挤压变换的心肺耦合特征提取方法与系统。

背景技术

睡眠是人体重要生理过程，在个体认知记忆、机体生长、心理健康等方面具有重要作用。长期的睡眠缺乏或睡眠障碍，可能引起精力不足和嗜睡等现象，甚至可能导致心脑血管疾病的发生。睡眠呼吸暂停疾病（Sleep Apnea）是最常见的睡眠障碍之一，严重影响健康，据统计，全球约10亿人患有该疾病。睡眠健康问题日益成为突出的医疗及公共卫生问题。然而，多导睡眠图作为临床诊断睡眠障碍的黄金标准，却面临着成本高昂，过多的外加硬件设备对人体造成负荷等问题。当前，针对便携式睡眠健康监测的重大发展需求，越来越多的研究致力于开发更精准、更高效的算法或参数评估睡眠障碍性疾病。

目前，学者们提出了多种基于体表心电信号中心脏跳动间隔（R-R）检测睡眠呼吸事件的方法。然而，这些方法在应用于心率变异性较低的受试者时会受到限制，例如受药物治疗影响或长期心率变异性低的个体。此外，研究表明仅对单一生理变量如心率变异性进行分析，无法准确和全面地描述自主神经调节功能状态。因此，近年来，基于心电信号的心肺耦合研究受到了广泛关注。哈佛医学院研究团队于2005年首次提出了心肺耦合（CPC）概念，通过融合心率和呼吸时间序列双模态信息，准确衡量心肺耦合强度，实现睡眠质量定量测量并检测睡眠呼吸事件。与基于传统多导睡眠图的睡眠诊断技术相比，心肺耦合数据的采集方法简单便携，在可穿戴健康管理以及临床辅助诊断领域具有巨大发展潜力。

然而，原始的心肺耦合算法基于傅里叶分析，傅里叶分析本质假设输入信号为平稳信号，而心电信号通常为非平稳、非线性信号，无法满足这一要求。此外，使用傅里叶变换无法同时获得高频率分辨率和高时间分辨率，不可避免地导致图谱结果中状态变化的模糊。因此，该方法通常只适合提供睡眠质量的一般信息（如给出检测为睡眠呼吸紊乱的记录时间百分比），而不是精确的睡眠呼吸紊乱检测。另一种常见的信号处理方法为小波变换法，与傅里叶变换相比，小波变换具有自适应时频窗口特性，然而由于变换过程中母小波的引入，导致生成的时频分布通常无法收拢于真实瞬时频率，因此，大部分情况下基于小波变换的时频谱较为模糊。此外，受限于海森堡测不准原理，基于小波变换输出的时频图谱的时间分辨率和频率分辨率总是相互制约，即无法在同一时刻保证高时间分辨率和高频率分辨率。为了解决上述问题，一些研究提出将新的时频分析方法应用于心肺耦合分析中，以提高心肺耦合有效性，这些时频分析方法包括时频重分配法、基于经验模态分解的希尔伯特-黄变换等。然而，这些方法依然存在各自弊端。时频重分配法虽然可以改善时频分析结果，然而在重建、分解信号上具有明显劣势。经验模态分解，虽具有较好自适应性，但若原始信号存在较大程度噪声，则使用该算法获取的本征模态函数易存在模态混叠缺点，即无法通过常规分解提取理想情况下的单模态分量，从而导致后续分析错误。此外，尽管许多科学家致力于研究经验模态分解的数学理论基础，该方法目前尚缺乏完善数学理论依据。上述存在问题导致当前心肺耦合图谱常存在模糊现象，可视化效果较差，关键节律动态变化无法获取。目前该问题尚未得到有效解决。

发明内容

有鉴于此，本发明提供了一种基于同步挤压变换的心肺耦合特征提取方法与系统，能够突破传统心肺耦合图谱关键节律模糊的技术难点，实现对睡眠呼吸暂停事件的检测。

为达到上述目的，本发明的技术方案包括如下步骤。

步骤1：从单导联心电信号中提取心跳间期R-R时间序列和替代呼吸EDR时间序列。

步骤2：对心跳间期R-R时间序列和替代呼吸EDR时间序列分别进行同步挤压变换，分别得到R-R时间序列的同步挤压系数和EDR时间序列的同步挤压系数。

步骤3：利用R-R时间序列的同步挤压系数和EDR时间序列的同步挤压系数计算心肺耦合CPC值，绘制睡眠期间心肺耦合图谱。

步骤4：基于睡眠期间心肺耦合图谱提取多类心肺耦合特征，利用机器学习分类的方法对所提取的多类心肺耦合特征进行学习分类，以识别正常呼吸事件和睡眠呼吸暂停事件。

进一步地，步骤1中，单导联心电信号为针对待测的睡眠期间心电信号进行预处理获得，针对待测的睡眠期间心电信号进行预处理的具体过程包括：对整条心电信号按照每分钟进行分段，对分段后的心电信号使用带通滤波器，保留0.5-30Hz频段，移除其它频段；由此获得单导联心电信号。

进一步地，步骤1中，从单导联心电信号中提取心跳间期R-R时间序列和替代呼吸EDR时间序列，具体为：

在单导联心电信号中标记心跳，计算相邻两次心跳时间间隔，得到初始R-R时间序列信号。

在固定时间间隔内估计心电图QRS波群与校正后心电图基线所包围区域的面积，衡量该波群面积变化，获得初始EDR时间序列信号。

针对初始R-R时间序列信号和初始EDR时间序列信号，使用41个数据点窗口的滑动平均滤波器去除由于错误心跳检测产生的异常值，当窗口中的中心点位于平均值20%之外时，将被移除；最后，使用三次样条插值以4Hz频率均匀地对R-R和EDR这两个信号进行重新采样，重采样之后二者采样间隔均为∆ t，获得最终提取的心跳间期R-R时间序列和替代呼吸EDR时间序列。

进一步地，步骤1中，使用三次样条插值以4Hz频率均匀地对R-R和EDR这两个信号进行重新采样，重采样之后二者采样间隔均为∆ t，之后还包括：在重采样之后的R-R时间序列和重采样之后的EDR时间序列的两侧均进行填充，使填充后信号长度为2的幂次倍。

进一步地，步骤2，对心跳间期R-R时间序列和替代呼吸EDR时间序列分别进行同步挤压变换，分别得到R-R时间序列的同步挤压系数和EDR时间序列的同步挤压系数，包括如下具体步骤：

S201：心跳间期R-R时间序列和替代呼吸EDR时间序列进行连续小波变换，分别得到R-R时间序列的小波系数 W _R( a, t)和EDR时间序列的小波系数 W _E( a, t)，表达式为：

(2)

(3)

其中 a为尺度， t为时间，Ψ(·)为母小波函数， *表示取共轭；在这一步中，一个关键参数为母小波函数的种类Ψ，其可以设置为Morlet、Bump或Morse小波函数等多种形式； R( u)、 E( u)分别为心跳间期R-R时间序列和替代呼吸EDR时间序列， u为积分变量；

S202：由小波系数求出瞬时频率，其表达式为：

(4)

(5)

其中 ω _R( a, t)、 ω _E( a, t)分别为心跳间期R-R时间序列和替代呼吸EDR时间序列于每一尺度 a和时间 t对应的瞬时频率； γ为可自定义噪声门限，Im表示取虚部； Inf表示无穷大；

S203：针对心跳间期R-R时间序列和替代呼吸EDR时间序列的对频率轴进行重排挤压后得到同步挤压系数矩阵；

其中同步挤压系数表达式为：

(6)

(7)

其中， T _R( ω _l, t)和 T _E( ω _l, t)为最终获得的同步挤压系数， W _R( a, t)、 W _E( a, t)分别为心跳间期R-R时间序列和替代呼吸EDR时间序列的小波系数， ω _R( a, t)、 ω _E( a, t)则分别为心跳间期R-R时间序列和替代呼吸EDR时间序列的瞬时频率， t为时间， a为尺度， B _l为根据瞬时频率划分的统计区间， ω _l为划分后的瞬时频率；{ a: ω _R( a, t)∈ B _l}、{ a: ω _E( a, t)∈ B _l}分别表示针对心跳间期R-R时间序列和替代EDR时间序列将其频率统计区间 B _l内的所有瞬时频率挤压到指定频率点 ω _l；

进一步地，步骤3中，计算心肺耦合CPC值，绘制睡眠期间心肺耦合图谱，具体为：

首先计算心跳间期R-R时间序列和替代呼吸EDR时间序列的互功率谱，记作Γ( R, E)，表达式如下：

(8)

其中，和分别为R-R和EDR时间序列的同步挤压系数；*表示复共轭；

计算这两个信号之间的相干性，记为 C；相干性的定义为平均交叉谱的平方除以单个信号的平均谱功率的乘积，即：

(9)

其中<>表示对原始谱图中的频率进行平均或对给定频率下的多个测量值进行平均；

最后，计算给定频率下的心肺耦合指标 CPC：

CPC=<Γ( R, E)>²· C (10)

基于上述心肺耦合CPC值，绘制睡眠期间心肺耦合图谱。

进一步地，步骤4，基于睡眠期间心肺耦合图谱提取多类心肺耦合特征，多类心肺耦合特征包含；总功率TP、高频功率HF、低频功率LF、超低频功率VLF、低高频功率比LF/HF、标准化低频功率LFnorm和标准化高频功率Hfnorm；其中，总功率TP定义为频段范围0-0.40Hz内所有心肺耦合值之和；HF，LF和VLF分别为高频带0.15-0.40Hz，低频带0.04-0.15Hz和超低频带0.003-0.04Hz内所有心肺耦合值之和；LF/HF为低频功率LF与高频功率HF的比值；LFnorm为低频功率LF与总功率TP的比值；HFnorm为高频功率HF与总功率TP的比值。

进一步地，利用机器学习分类的方法对所提取的多类心肺耦合特征进行学习分类之前还包括使用最小冗余-最大相关算法对全部提取的心肺耦合特征进行排序。

本发明另外一个实施例还提供了基于同步挤压变换的心肺耦合特征提取系统，包括时间序列提取模块、同步挤压变换模块、心肺耦合图谱绘制模块、心肺耦合特征提取模块以及机器学习模块；

时间序列提取模块，用于接收单导联心电信号，并从中提取心跳间期R-R时间序列和替代呼吸EDR时间序列；

同步挤压变换模块，用于对心跳间期R-R时间序列和替代呼吸EDR时间序列分别进行同步挤压变换，分别得到R-R时间序列的同步挤压系数和EDR时间序列的同步挤压系数。

心肺耦合图谱绘制模块，用于利用R-R时间序列的同步挤压系数和EDR时间序列的同步挤压系数计算心肺耦合CPC值，绘制睡眠期间心肺耦合图谱。

心肺耦合特征提取模块，用于基于睡眠期间心肺耦合图谱提取多类心肺耦合特征。

机器学习模块，用于利用机器学习分类的方法对所提取的多类心肺耦合特征进行学习分类，以识别正常呼吸事件和睡眠呼吸暂停事件。

有益效果：

1.本发明提供的一种基于同步挤压变换的心肺耦合特征提取方法，基于同步挤压变换于构建可视化图谱的高分辨率优势，首次提出了一套基于同步挤压变换的心肺耦合特征提取新算法，包括以下步骤：对采集到的原始单导联心电信号进行预处理，避免工频噪声以及无关频率成分对结果造成影响；从预处理后的心电信号提取心脏跳动间隔（R-R）间期时间序列和替代呼吸（EDR）时间序列；利用同步挤压变换算法分别对EDR时间序列和R-R间期时间序列进行处理，获得上述两时间序列的同步挤压系数；根据R-R和EDR时间序列的同步挤压系数计算心肺耦合（CPC）值，构建睡眠期间心肺耦合图谱；基于睡眠心肺耦合图谱提取多类心肺耦合特征，进行特征选择，完成特征重要性排序；最终，构建有效特征子集作为机器学习分类器的输入，输出当前时段睡眠呼吸暂停事件的判断结果。本发明基于所提出新方法实现了对睡眠呼吸暂停事件的精准检测，突破了传统心肺耦合图谱关键节律模糊的技术难点。

2.本发明提供的一种基于同步挤压变换的心肺耦合特征提取系统，基于以上提出一套基于同步挤压变换的心肺耦合特征提取新方法，适用于睡眠呼吸暂停事件的精准检测，高效可靠，易于软件化，通过对方法的功能分析，设计软件功能模块，使得该方法易于实现，可望为便携式睡眠健康监测提供可行新途径。

附图说明

图1为本发明提供的一种基于同步挤压变换的心肺耦合特征提取方法技术流程图。

图2为基于同步挤压变换的心肺耦合算法的部分处理结果，包括单导联心电信号及生成的心脏跳动间隔R-R间期和替代呼吸EDR信号。

具体实施方式

下面结合附图并举实施例，对本发明进行详细描述。

本发明提出基于同步挤压变换（SST）的心肺耦合分析技术。SST作为较新颖的时频处理算法之一，适用于处理非线性和非平稳时间序列。该技术具有完善的数学理论依据，其核心思想可概括为找寻某一时刻的瞬时频率，进而将此前变换所得结果“挤压”到该频率附近，从而避免结果的弥散现象，提供更优的时频聚焦特性。与经验模态分解相比，SST克服了其缺乏数学理论依据的缺点。此外，SST算法支持由多模态频域分量重构信号，亦体现了时频重分配法不具有的独特优势。事实上，SST算法在生物医学信号处理中有重要且广泛的应用，可用于从信号中提取有用特征或量化其动力学行为，以利后续统计分析。目前，SST已被用于单导联心电图自动QRS波群检测，从脑电图预测情绪状态，心房颤动患者的EDR信号提取等研究。然而，SST算法尚未被应用于CPC分析中。

本发明基于同步挤压变换算法优化传统心肺耦合分析技术，创新性地提出一套基于同步挤压变换的心肺耦合算法，并将所提出算法应用于睡眠呼吸暂停事件的检测。

本发明的技术流程图如图1所示，详细过程如下：

步骤1）采集受测者睡眠过程中的单导联心电信号，采样频率设为 f _s，存储在计算机中。本发明中可以利用开源数据集，例如Physionet数据库中的Sleep Apnea数据集。

对上述原始心电信号进行预处理，详细流程包括：对整条心电信号按照每分钟进行分段，避免数据过长导致计算效率下降；对分段后的心电信号使用带通滤波器，保留0.5-30Hz频段（心电信号处理中的常规频段），移除其它频段，避免工频噪声以及无关频率分量对结果造成影响。

从预处理后的心电信号中提取心跳间期（R-R）时间序列和替代呼吸（EDR）时间序列。具体来说，在单导联心电信号中标记心跳，计算相邻两次心跳时间间隔，得到R-R时间序列；在固定时间间隔内估计心电图QRS波群与校正后心电图基线所包围区域的面积，衡量该波群面积变化，获得EDR时间序列信号。在从心电信号中提取正常R-R间期时间序列和EDR信号后，使用41个数据点窗口的滑动平均滤波器去除由于错误心跳检测产生的异常值。当窗口中的中心点位于平均值20%之外时，将被移除。最后，使用三次样条插值以4Hz频率均匀地对R-R和EDR这两个信号进行重新采样。重采样之后二者采样间隔均为∆ t，获得最终提取的心跳间期R-R时间序列和替代呼吸EDR时间序列。

步骤2）对心跳间期R-R时间序列和替代呼吸EDR时间序列分别进行同步挤压变换，具体来说：

通过对心电信号进行预处理，可由心电信号获得：（1）心跳间期R-R信号，记作 R( n)，和（2）替代呼吸EDR信号，记作 E( n)。假设 R( n)和 E( n)的采样间隔为∆ t，信号长度为 N（ n为采样点，取值1到 N）。为了避免边界效应，在 R( n)和 E( n)的两侧进行填充，常用的边界填充方法有零填充、常数填充、对称填充和周期填充等。同时，为提高计算效率，通常使填充后信号长度为2的幂次倍，这里假设填充后信号长度为 T=2^L+1，其中 L为非负整数。

本步骤2包括如下具体步骤：

S201：对R-R和EDR信号进行连续小波变换得到小波系数。连续小波变换主要思想可概括为在一系列离散的尺度 a _j上，计算原始信号与母小波函数的卷积。尺度 a _j通常以对数划分，这也是小波变换的一大优点，即在低频拥有更高的频率分辨率。小波变换尺度划分如式（1）所示，其中 n _v为自定义常数，影响频率离散化程度，即小波变换时下一频率等于前一频率乘以2^{1/
nv}。在实际应用中， n _v通常取32或64。

(1)

其中 j取值1到 L×n _v，表示共 L×n _v个尺度。

随后，分别对R-R和EDR时间序列计算小波系数，表达式为：

(2)

(3)

其中 W _R( a, t)、 W _E( a, t)分别为心跳间期R-R时间序列和替代呼吸EDR时间序列的小波系数， a为尺度， a取一系列值 a _j，记作 a={ a _j, j=1,2,…, L×n _v}， t为时间，Ψ(·)为母小波函数， *表示取共轭， u为积分变量。在这一步中，一个关键参数为母小波函数的种类Ψ，其可以设置为Morlet、Bump或Morse小波函数等多种形式。公式（2）（3）为小波系数常规计算方式。

S202：由小波系数求出瞬时频率，其表达式为：

(4)

(5)

其中 ω _R( a, t)、 ω _E( a, t)分别为心跳间期R-R时间序列和替代呼吸EDR时间序列于每一尺度 a和时间 t对应的瞬时频率。 γ为可自定义噪声门限，Im表示取虚部， Inf表示无穷大。噪声门限 γ的作用为减少噪声波动对结果造成影响，当某个时频点的小波系数极其微弱时，则认为是由噪声所引起。因此，通过将该点的瞬时频率设置为无穷大来忽略该点数据。

S203：对频率轴进行“重排挤压”后得到同步挤压系数矩阵。首先，将频率轴划分为一系列指定频率点，其中 ω ₀>0且 ω _l+1> ω _l。这些指定频率点表示最终希望将时频能量“挤压”到的位置，其划分方式支持用户自定义，可为线性或对数划分。在设置指定频点时应注意，受到奈奎斯特采样定律的限制，可观测到的最大频率为 ω=1/2∆ t。同时，由于采样时间有限，可设置的最小频率为 ω =1/ n∆ t。因此，以对数划分方式为例， ω _l可设置为2^{l∆
ω} ω ,l=1,2,…, Ln _v-1，其中，间隔∆ ω=1/( Ln _v-1)log₂( n/2)。进一步，由这些频率点可确定若干统计区间，记作 B _l={ ω'∈ R:| ω'-ω _l| <| ω'-ω _l'|, l≠ l'}，该式表明对于落入相同统计区间 B _l内的任何瞬时频率 ω'，都将被划分至同一最近指定频点 ω _l，而不是其它距离 ω'更远的频点 ω _l'。随后，通过将某一频率统计区间 B _l内的所有瞬时频率“挤压”到指定频率点 ω _l，提高时频分辨率，获得更优可视化时频图谱，表达式为：

(6)

(7)

其中， T _R( ω _l, t)和 T _E( ω _l, t)为最终获得的同步挤压系数， W _R( a, t)、 W _E( a, t)分别为心跳间期R-R时间序列和替代呼吸EDR时间序列的小波系数（式（2）和式（3））， ω _R( a, t)、 ω _E( a, t)则分别为心跳间期R-R时间序列和替代呼吸EDR时间序列的瞬时频率（式（4）和式（5））， t为时间， a为尺度， B _l为根据瞬时频率划分的统计区间， ω _l为划分后的瞬时频率。{ a: ω _R( a, t)∈ B _l}、{ a: ω _E( a, t)∈ B _l}分别表示针对心跳间期R-R时间序列和替代EDR时间序列将其频率统计区间 B _l内的所有瞬时频率挤压到指定频率点 ω _l。

步骤3）利用R-R时间序列的同步挤压系数和EDR时间序列的同步挤压系数计算心肺耦合（CPC）值，绘制睡眠期间心肺耦合图谱。

首先计算心跳间期（R-R）时间序列和替代呼吸（EDR）时间序列的互功率谱，记作Γ( R, E)，表达式如下：

(8)

其中，和分别为R-R和EDR时间序列的同步挤压系数。*表示复共轭。

另一方面，计算这两个信号之间的相干性，记为 C。相干性的定义为平均交叉谱的平方除以单个信号的平均谱功率的乘积，即：

(9)

其中<>表示对原始谱图中的频率进行平均或对给定频率下的多个测量值进行平均。这是因为相干性是一种统计性度量，所以进行统计平均是必须的。在原始文献中，通过将每个观察窗口划分为子窗口，再对子窗口求平均获得平均谱功率。

最后，计算给定频率下的心肺耦合指标，该指标通过整合交叉谱功率和相干性来量化R-R和EDR时间序列之间的耦合强度，心肺耦合指标定义如下：

CPC=<Γ( R, E)>²· C (10)

基于上述心肺耦合CPC值，绘制睡眠期间心肺耦合图谱。通过观察心肺耦合图谱中的耦合模式，识别受测者是否患有睡眠呼吸暂停疾病。

步骤4）基于睡眠期间心肺耦合图谱提取多类心肺耦合特征，利用机器学习分类的方法对所提取的多类心肺耦合特征进行学习分类，以识别正常呼吸事件和睡眠呼吸暂停事件。

基于睡眠期间心肺耦合图谱提取多类心肺耦合特征。在心肺耦合分析中，通常在获得完整的心肺耦合图谱后，将高频带（0.15-0.40Hz），低频带（0.04-0.15Hz）和超低频带（0.003-0.04Hz）内的所有心肺耦合值分别求和作为各个频段的心肺耦合功率，记作HF，LF和VLF，并计算低频带心肺耦合值与高频带心肺耦合值的比值LF/HF，从而为后续睡眠呼吸暂停事件的检测提供参考。此外，对上述指标进行标准化以后获得的新指标也被广泛应用于与睡眠相关的研究中。表1归纳了常见的心肺耦合指标及其对应的频段范围。由表1可见，TP定义为指定频段范围内所有心肺耦合值之和，即总功率；HF、LF和VLF则依次代表高频、低频、和超低频范围内的功率；LFnorm和HFnorm分别表示标准化的低频和高频功率；LF/HF则定义为低频与高频功率比。

表1心肺耦合特征提取

进一步地，使用最小冗余-最大相关算法对全部提取特征进行排序，将所选特征输入多种机器学习分类模型，进行每分钟睡眠呼吸事件检测。

以上提供的基于同步挤压变换的心肺耦合特征提取方法，适用于睡眠呼吸暂停事件的精准检测，高效可靠，易于软件化，使用该方法时，可以根据每个步骤的功能，进行模块化软件设计，具体如下：

基于同步挤压变换的心肺耦合特征提取系统，包括时间序列提取模块、同步挤压变换模块、心肺耦合图谱绘制模块、心肺耦合特征提取模块以及机器学习模块。

时间序列提取模块，用于接收单导联心电信号，并从中提取心跳间期R-R时间序列和替代呼吸EDR时间序列。

本发明已在Physionet数据库中的Sleep Apnea数据集上进行验证。通过绘制不同受试者于睡眠时期的心肺耦合图谱，发现与传统心肺耦合算法相比，基于同步挤压变换的心肺耦合图谱具有更优可视化效果。此外，新型心肺耦合特征与睡眠呼吸暂停严重程度（即睡眠呼吸暂停低通气指数—AHI值）存在显著相关性，以及在不同疾病严重程度的组间显著统计学差异，为后续睡眠呼吸暂停检测提供支撑。综上，基于同步挤压变换的心肺耦合算法可作为临床常规诊断的补充工具，有助于提高睡眠呼吸事件诊断能力，在便携可穿戴式睡眠健康监测中具有一定的潜在价值和应用前景。

综上所述，以上仅为本发明的较佳实施例而已，并非用于限定本发明的保护范围。凡在本发明的精神和原则之内，所作的任何修改、等同替换、改进等，均应包含在本发明的保护范围之内。

Claims

1.基于同步挤压变换的心肺耦合特征提取方法，其特征在于，包括如下步骤：

步骤1：从单导联心电信号中提取心跳间期R-R时间序列和替代呼吸EDR时间序列；

步骤2：对所述心跳间期R-R时间序列和所述替代呼吸EDR时间序列分别进行同步挤压变换，分别得到R-R时间序列的同步挤压系数和EDR时间序列的同步挤压系数；

步骤3：利用所述R-R时间序列的同步挤压系数和EDR时间序列的同步挤压系数计算心肺耦合CPC值，绘制睡眠期间心肺耦合图谱，具体为：

首先计算心跳间期R-R时间序列和替代呼吸EDR时间序列的互功率谱，记作Γ(R,E)，表达式如下：

(8)

计算这两个信号之间的相干性，记为C；相干性的定义为平均交叉谱的平方除以单个信号的平均谱功率的乘积，即：

(9)

最后，计算给定频率下的心肺耦合指标CPC：

CPC=<Γ(R,E)>²·C (10)

基于上述心肺耦合CPC值，绘制睡眠期间心肺耦合图谱；

2.如权利要求1所述的基于同步挤压变换的心肺耦合特征提取方法，其特征在于，所述步骤1中，所述单导联心电信号为针对待测的睡眠期间心电信号进行预处理获得，针对待测的睡眠期间心电信号进行预处理的具体过程包括：对整条心电信号按照每分钟进行分段，对分段后的心电信号使用带通滤波器，保留0.5-30Hz频段，移除其它频段；由此获得单导联心电信号。

3.如权利要求1所述的基于同步挤压变换的心肺耦合特征提取方法，其特征在于，所述步骤1中，所述从单导联心电信号中提取心跳间期R-R时间序列和替代呼吸EDR时间序列，具体为：

在单导联心电信号中标记心跳，计算相邻两次心跳时间间隔，得到初始R-R时间序列信号；

在固定时间间隔内估计心电图QRS波群与校正后心电图基线所包围区域的面积，衡量该波群面积变化，获得初始EDR时间序列信号；

针对初始R-R时间序列信号和初始EDR时间序列信号，使用41个数据点窗口的滑动平均滤波器去除由于错误心跳检测产生的异常值，当窗口中的中心点位于平均值20%之外时，将被移除；最后，使用三次样条插值以4Hz频率均匀地对R-R和EDR这两个信号进行重新采样，重采样之后二者采样间隔均为∆t，获得最终提取的心跳间期R-R时间序列和替代呼吸EDR时间序列。

4.如权利要求3所述的基于同步挤压变换的心肺耦合特征提取方法，其特征在于，所述步骤1中，使用三次样条插值以4Hz频率均匀地对R-R和EDR这两个信号进行重新采样，重采样之后二者采样间隔均为∆t，之后还包括：在重采样之后的R-R时间序列和重采样之后的EDR时间序列的两侧均进行填充，使填充后信号长度为2的幂次倍。

5.如权利要求1所述的基于同步挤压变换的心肺耦合特征提取方法，其特征在于，所述步骤2，对所述心跳间期R-R时间序列和所述替代呼吸EDR时间序列分别进行同步挤压变换，分别得到R-R时间序列的同步挤压系数和EDR时间序列的同步挤压系数，包括如下具体步骤：

S201：所述心跳间期R-R时间序列和所述替代呼吸EDR时间序列进行连续小波变换，分别得到R-R时间序列的小波系数W _R(a,t)和EDR时间序列的小波系数W _E(a,t)，表达式为：

(2)

(3)

其中a为尺度，t为时间，Ψ(·)为母小波函数，*表示取共轭；在这一步中，一个关键参数为母小波函数的种类Ψ，设置为Morlet、Bump或Morse小波函数；R(u)、E(u)分别为心跳间期R-R时间序列和替代呼吸EDR时间序列，u为积分变量；

S202：由小波系数求出瞬时频率，其表达式为：

(4)

(5)

其中ω _R(a,t)、ω _E(a,t)分别为心跳间期R-R时间序列和替代呼吸EDR时间序列于每一尺度a和时间t对应的瞬时频率；γ为可自定义噪声门限，Im表示取虚部；Inf表示无穷大；

其中同步挤压系数表达式为：

(6)

(7)

其中，T _R(ω _l,t)和T _E(ω _l,t)为最终获得的同步挤压系数，W _R(a,t)、W _E(a,t)分别为心跳间期R-R时间序列和替代呼吸EDR时间序列的小波系数，ω _R(a,t)、ω _E(a,t)则分别为心跳间期R-R时间序列和替代呼吸EDR时间序列的瞬时频率，t为时间，a为尺度，B _l为根据瞬时频率划分的统计区间，ω _l为划分后的瞬时频率；{a:ω _R(a,t)∈B _l}、{a:ω _E(a,t)∈B _l}分别表示针对心跳间期R-R时间序列和替代EDR时间序列将其频率统计区间B _l内的所有瞬时频率挤压到指定频率点ω _l。

6.如权利要求1所述的基于同步挤压变换的心肺耦合特征提取方法，其特征在于，所述步骤4，基于睡眠期间心肺耦合图谱提取多类心肺耦合特征，所述多类心肺耦合特征包含：总功率TP、高频功率HF、低频功率LF、超低频功率VLF、低高频功率比LF/HF、标准化低频功率LFnorm和标准化高频功率Hfnorm；其中，总功率TP定义为频段范围0-0.40Hz内所有心肺耦合值之和；HF，LF和VLF分别为高频带0.15-0.40Hz，低频带0.04-0.15Hz和超低频带0.003-0.04Hz内所有心肺耦合值之和；LF/HF为低频功率LF与高频功率HF的比值；LFnorm为低频功率LF与总功率TP的比值；HFnorm为高频功率HF与总功率TP的比值。

7.如权利要求6所述的基于同步挤压变换的心肺耦合特征提取方法，其特征在于，所述利用机器学习分类的方法对所提取的多类心肺耦合特征进行学习分类之前还包括使用最小冗余-最大相关算法对全部提取的心肺耦合特征进行排序。

8.基于同步挤压变换的心肺耦合特征提取系统，其特征在于，包括时间序列提取模块、同步挤压变换模块、心肺耦合图谱绘制模块、心肺耦合特征提取模块以及机器学习模块；

所述时间序列提取模块，用于接收单导联心电信号，并从中提取心跳间期R-R时间序列和替代呼吸EDR时间序列；

所述同步挤压变换模块，用于对所述心跳间期R-R时间序列和所述替代呼吸EDR时间序列分别进行同步挤压变换，分别得到R-R时间序列的同步挤压系数和EDR时间序列的同步挤压系数；

所述心肺耦合图谱绘制模块，用于利用所述R-R时间序列的同步挤压系数和EDR时间序列的同步挤压系数计算心肺耦合CPC值，绘制睡眠期间心肺耦合图谱；具体为：

(8)

(9)

最后，计算给定频率下的心肺耦合指标CPC：

CPC=<Γ(R,E)>²·C (10)

基于上述心肺耦合CPC值，绘制睡眠期间心肺耦合图谱；

所述心肺耦合特征提取模块，用于基于睡眠期间心肺耦合图谱提取多类心肺耦合特征；

所述机器学习模块，用于利用机器学习分类的方法对所提取的多类心肺耦合特征进行学习分类，以识别正常呼吸事件和睡眠呼吸暂停事件。