CN115886742B - 一种基于变分心肺耦合的睡眠动力学分析方法 - Google Patents

一种基于变分心肺耦合的睡眠动力学分析方法 Download PDF

Info

Publication number
CN115886742B
CN115886742B CN202310227600.0A CN202310227600A CN115886742B CN 115886742 B CN115886742 B CN 115886742B CN 202310227600 A CN202310227600 A CN 202310227600A CN 115886742 B CN115886742 B CN 115886742B
Authority
CN
China
Prior art keywords
function
eigenmode
time sequence
frequency
coupling
Prior art date
Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
Active
Application number
CN202310227600.0A
Other languages
English (en)
Other versions
CN115886742A (zh
Inventor
史文彬
冯欢
叶建宏
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
Beijing Institute of Technology BIT
Original Assignee
Beijing Institute of Technology BIT
Priority date (The priority date is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the date listed.)
Filing date
Publication date
Application filed by Beijing Institute of Technology BIT filed Critical Beijing Institute of Technology BIT
Priority to CN202310227600.0A priority Critical patent/CN115886742B/zh
Publication of CN115886742A publication Critical patent/CN115886742A/zh
Application granted granted Critical
Publication of CN115886742B publication Critical patent/CN115886742B/zh
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Images

Classifications

    • YGENERAL TAGGING OF NEW TECHNOLOGICAL DEVELOPMENTS; GENERAL TAGGING OF CROSS-SECTIONAL TECHNOLOGIES SPANNING OVER SEVERAL SECTIONS OF THE IPC; TECHNICAL SUBJECTS COVERED BY FORMER USPC CROSS-REFERENCE ART COLLECTIONS [XRACs] AND DIGESTS
    • Y02TECHNOLOGIES OR APPLICATIONS FOR MITIGATION OR ADAPTATION AGAINST CLIMATE CHANGE
    • Y02ATECHNOLOGIES FOR ADAPTATION TO CLIMATE CHANGE
    • Y02A90/00Technologies having an indirect contribution to adaptation to climate change
    • Y02A90/10Information and communication technologies [ICT] supporting adaptation to climate change, e.g. for weather forecasting or climate simulation

Landscapes

  • Measurement And Recording Of Electrical Phenomena And Electrical Characteristics Of The Living Body (AREA)

Abstract

本发明涉及一种基于变分心肺耦合的睡眠动力学分析方法,属于睡眠医学和信息技术领域。本发明对心跳间隔时间序列和心电图衍生呼吸信号进行变分模态分解,具有对噪声和采样率的鲁棒性,并且频带分离效果好,解决了模态混叠问题,能将信号分解成几乎没有频带重叠的成分,能避免后续画出的心肺耦合图谱出现虚假耦合的现象。本发明结合希尔伯特变换能极大地提高心肺耦合图谱的时频分辨率,得到细节化的心肺耦合图谱和准确的心肺耦合指数,能提高诊断阻塞性睡眠呼吸暂停低通气综合征的准确性。

Description

一种基于变分心肺耦合的睡眠动力学分析方法
技术领域
本发明涉及一种基于变分心肺耦合的睡眠动力学分析方法,属于睡眠医学和信息技术领域。
背景技术
阻塞性睡眠呼吸暂停低通气综合征(OSAHS)是老年人常见病,并逐渐呈现年轻化趋势。OSAHS与高血压、心律失常、心力衰竭等心血管疾病有潜在相关性,严重情况下患者可能在睡眠过程中窒息死亡,因此,对OSAHS进行及时的诊断和治疗具有重要意义。多导睡眠监测是诊断OSAHS的权威方法。但是,多导睡眠监测需要专业的大型设备,价格昂贵。此外,在使用多导睡眠监测获得脑电图后,需要专业的睡眠技师花费大量时间手动识别脑电波的特征波形,识别睡眠中发生的呼吸暂停或低通气事件。因此,有必要找到一种更加准确便捷的方法来研究OSAHS。
目前根据OSAHS患者睡眠期间心跳间隔时间序列的变化表现出与正常人不同的特征这一特点,出现了大量利用能代表心跳间隔时间序列变化规律的心率变异性来诊断OSAHS的方法。然而,某些疾病(例如充血性心力衰竭或自主神经病变)和药物(例如副交感神经溶解剂)与心率变异性的显著降低有关。因此仅根据心率变异性难以将OSAHS区分于以上这些疾病和用药情况。除此之外,从体表心电图中提取的心电图衍生呼吸(EDR)信号,能作为OSAHS研究的补充信息,因此研究人员提出结合心率变异性和EDR信号的心肺耦合算法(CPC)来检测睡眠中的呼吸暂停-低通气事件。但传统的心肺耦合算法基于短时傅里叶变换或小波变换,这一算法的时频分辨率较低,无法提供可靠的瞬时频率,不可避免地模糊了睡眠图谱的细节。后来提出的基于希尔伯特-黄变换的心肺耦合算法虽然在一定程度上提高了时频分辨率,但存在一些明显的缺点,如对噪声和采样率的敏感性以及模态混叠问题,在睡眠图谱中会出现虚假的耦合结果,降低了对OSAHS诊断的准确性。
发明内容
本发明的技术解决问题是:克服现有技术的不足,提出一种基于变分心肺耦合的睡眠动力学分析方法,该方法将变分模态分解和希尔伯特变换结合到心肺耦合技术中,提出的一种新的心肺动力学测量方法——变分心肺耦合(VMD-CPC)。
本发明的技术解决方案是:
一种基于变分心肺耦合的睡眠动力学分析方法,该方法的步骤包括:
第一步,采集受试者的心电信号,并对采集到的心电信号进行预处理,去除杂讯和干扰;
第二步,提取第一步预处理后心电信号的心跳间隔时间序列(R-R间期信号)和心电图衍生呼吸信号(EDR信号);
第三步,对第二步提取的心跳间隔时间序列(R-R间期信号)和心电图衍生呼吸信号(EDR信号)分别进行变分模态分解,得到心跳间隔时间序列(R-R间期信号)的本征模态函数集合和心电图衍生呼吸信号(EDR信号)的本征模态函数集合;
第四步,对第三步得到的心跳间隔时间序列(R-R间期信号)的本征模态函数集合进行希尔伯特变换,然后根据希尔伯特变换结果计算心跳间隔时间序列(R-R间期信号)的本征模态函数集合的瞬时包络和瞬时相位,最后再根据计算得到的瞬时相位得到心跳间隔时间序列(R-R间期信号)的本征模态函数集合的瞬时频率;
第五步,对第三步得到的心电图衍生呼吸信号(EDR信号)的本征模态函数集合进行希尔伯特变换,然后根据希尔伯特变换结果计算心电图衍生呼吸信号(EDR信号)的本征模态函数集合的瞬时包络和瞬时相位,最后再根据计算得到的瞬时相位得到心电图衍生呼吸信号(EDR信号)的本征模态函数集合的瞬时频率;
第六步,根据第四步得到的心跳间隔时间序列(R-R间期信号)的本征模态函数集合的瞬时频率将心跳间隔时间序列(R-R间期信号)的本征模态函数的时间序列转换为以瞬时频率和时间为因变量的二维解析函数;
第七步,根据第五步得到的心电图衍生呼吸信号(EDR信号)的本征模态函数集合的瞬时频率将心电图衍生呼吸信号(EDR信号)的本征模态函数的时间序列转换为以瞬时频率和时间为因变量的二维解析函数;
第八步,根据第六步得到的二维解析函数和第七步得到的二维解析函数计算得到两信号二维解析函数的互功率谱和相干性,再根据两信号二维解析函数的互功率谱和相干性得到心肺耦合指数;
第九步,以时间为x轴,频率为y轴,第八步得到的心肺耦合指数为z轴,以颜色不同代表心肺耦合指数的大小,得到可视化的睡眠期间心肺耦合图谱,并根据得到的心肺耦合图谱得到心肺耦合特征,完成基于变分心肺耦合的睡眠动力学分析。
所述第一步中,对采集的心电信号进行预处理时,使用50Hz的陷波滤波器去除心电信号中的工频干扰,使用中值滤波算法去除基线漂移,使用通带截止频率为80 Hz、阻带截止频率为100 Hz的巴特沃斯低通滤波器去除肌电干扰。
所述第二步中,提取的心电信号的心跳间隔时间序列为心电图上所有R峰间的时间间隔;
所述R峰的定位方法为心电信号的积分功率谱与心电信号的积分功率阈值的交点;
所述积分功率谱P(t)为:
Figure SMS_1
其中,x(t)代表心电信号,g(t-τ)代表高斯函数;
[δ 1,δ 2]为心电信号的频率范围;
f代表频率,t代表时间,τ代表高斯函数g(t)每次平移的步长;
心电信号的积分功率阈值的确定方法为:对积分功率谱进行滑动Tukey窗口加权,取每次加权数据中的最大值作为每个窗口的积分功率阈值。
所述第二步中,提取的心电信号衍生的呼吸信号为心电图上所有QRS波的面积。
所述第三步中,设心跳间隔时间序列的本征模态函数集合为{u R,k } ={u R,1, …,u R,KR },K R 为心跳间隔时间序列的本征模态函数的个数;本征模态函数集合的确定方法为求出变分模态分解过程中的约束变分问题的解;
所述约束变分问题为:
Figure SMS_2
其中,y R 为心跳间隔时间序列;
{ω R,k } ={ω R,1 , …,ω R,KR },表示这K R 个本征模态函数的中心频率的集合;
t代表时间,δ(t)代表狄拉克函数;j为虚数单位;
所述变分问题的解的确定方法为:
使用二次惩罚项α R 和拉格朗日乘子λ R (t),将约束变分问题转化为增广拉格朗日函数,所述的增广拉格朗日函数为:
Figure SMS_3
其中,α R 是一个预设值,λ R (t)将在求解过程中变化,初始值为0;
在傅里叶域内计算这个增广拉格朗日函数关于u R,k ω R,k 的鞍点,并在傅里叶域内迭代更新{u R,k }、{ω R,k }以及λ R (t)直到满足迭代停止条件,得到本征模态函数集合{u R,k }={u R,1, …,u R,KR };
n+1步的更新公式为:
Figure SMS_4
停止条件为
Figure SMS_5
其中, R,k (ω)代表u R,k (t)的傅里叶变换;
R (ω)代表y R 的傅里叶变换;
λ̂ R (ω)代表λ R (t)的傅里叶变换;
ε R 是一个预设值;
达到停止条件后得到的{u R,k } ={u R,1, …,u R,KR }即为心跳间隔时间序列的本征模态函数集合。
所述第三步中,设心电信号衍生的呼吸信号的本征模态函数集合为{u E,k } ={u E,1,…,u E,KE },本征模态函数集合的确定方法为求出变分模态分解过程中的约束变分问题的解;K E 为心跳间隔时间序列的本征模态函数的个数;
所述约束变分问题为:
Figure SMS_6
其中,y E 为心电图衍生呼吸信号;
{ω E,k } ={ω E,1 , …,ω E,KE },表示这K E 个本征模态函数的中心频率的集合;
t代表时间,δ(t)代表狄拉克函数;j为虚数单位;
所述变分问题的解的确定方法为:
使用二次惩罚项α E 和拉格朗日乘子λ E (t),将约束变分问题转化为增广拉格朗日函数,所述的增广拉格朗日函数为:
Figure SMS_7
其中,α E 是一个预设值,λ E (t)将在求解过程中变化,初始值为0;
在傅里叶域内计算这个增广拉格朗日函数关于u E,k ω E,k 的鞍点,并在傅里叶域内迭代更新{u E,k }、{ω E,k }以及λ E (t)直到满足迭代停止条件,得到本征模态函数集合{u E,k }={u E,1, …,u E,KE };
n+1步的更新公式为:
Figure SMS_8
停止条件为
Figure SMS_9
其中, E,k (ω)代表u E,k (t)的傅里叶变换;
E (ω)代表y E 的傅里叶变换;
λ̂ E (ω)代表λ E (t)的傅里叶变换;
ε E 是一个预设值;
达到停止条件后得到的{u E,k } ={u E,1, …,u E,KE }即为心电信号衍生的呼吸信号的本征模态函数集合。
所述第四步中,心跳间隔时间序列的本征模态函数集合的瞬时包络为:
Figure SMS_10
R i (t)表示心跳间隔时间序列,H R,i (t)表示心跳间隔时间序列的希尔伯特变换,K R 表示分解出的心跳间隔时间序列的本征模态函数集合中本征模态函数的个数;
心跳间隔时间序列的本征模态函数集合的瞬时相位为:
Figure SMS_11
心跳间隔时间序列的本征模态函数集合的瞬时频率为:
Figure SMS_12
所述第五步中,呼吸信号的本征模态函数集合的瞬时包络为:
Figure SMS_13
其中,E j (t)表示呼吸信号,H E,j (t)表示呼吸信号的希尔伯特变换,K E 表示分解出的心跳间隔时间序列的本征模态函数集合中本征模态函数的个数;
呼吸信号的本征模态函数集合的瞬时相位为:
Figure SMS_14
呼吸信号的本征模态函数集合的瞬时频率为:
Figure SMS_15
所述第六步中,二维解析函数A为:
Figure SMS_16
所述第七步中,二维解析函数B为:
Figure SMS_17
所述第八步中,互功率谱为:
Figure SMS_18
相干性为:
心肺耦合指数为:
Figure SMS_19
其中,<>表示在给定频率下多次测量的平均值。
有益效果
(1)本发明同时考虑OSAHS患者心率变异性和呼吸的关系,能将OSAHS与充血性心力衰竭、自主神经病变等疾病以及副交感神经溶解剂等用药情况区分开来,加快诊断速度、缩小诊断范围。且只需要便携式或可穿戴式心电监测工具,易于软件化。
(2)本发明对心跳间隔时间序列(R-R间期信号)和心电图衍生呼吸信号(EDR信号)进行变分模态分解,具有对噪声和采样率的鲁棒性,并且频带分离效果好,解决了模态混叠问题,能将信号分解成几乎没有频带重叠的成分,能避免后续画出的心肺耦合图谱出现虚假耦合的现象。
(3)本发明结合希尔伯特变换能极大地提高心肺耦合图谱的时频分辨率,得到细节化的心肺耦合图谱和准确的心肺耦合指数,能提高OSAHS诊断的准确性。
附图说明
图1为一名受试者的睡眠期间心肺耦合图谱;
图2为一名健康受试者及三名患者的睡眠期间心肺耦合图谱,以及四位受试者发生呼吸暂停或低通气事件的时间记录图。
具体实施方式
下面结合附图和实施例对本发明作进一步说明。
实施例
为了量化睡眠时心电和呼吸的耦合关系,提出变分心肺耦合方法,详细过程如下:
步骤1)受试者按照实验要求完成实验操作,采集受测者整晚睡眠数据,心电(ECG)信号采样频率为fs,存储在计算机中。
步骤2)对上述心电信号进行预处理。包括使用50 Hz的陷波滤波器去除原始心电信号中的工频干扰,使用中值滤波算法去除基线漂移,和使用通带截止频率为80 Hz,阻带截止频率为100 Hz的巴特沃斯低通滤波器去除肌电干扰。
步骤3)从预处理后的心电信号提取出代表心率变异性的R-R间期信号和代表呼吸的EDR信号。心电信号的成分是固定的,包括P波、QRS波和T波。R-R间期信号是心电信号中所有相邻R峰之间的时间间隔。具体过程为:对心电信号进行Gabor变换,得到心电信号的时频谱。在心电信号的频率范围[δ 1,δ 2]内计算它的积分功率得到其功率谱,功率谱为
Figure SMS_20
(1)
其中f代表频率,t代表时间,τ代表高斯函数g(t)每次平移的步长,x(t)代表心电信号,g(t-τ)代表高斯函数。
计算[δ 1,δ 2]频带内积分功率的阈值。阈值利用阶统计滤波器计算,具体做法为对积分功率进行滑动Tukey窗口加权,然后取加权数据中的最大值作为阈值。在用阶统计滤波方法计算阈值时,窗口每次向前移动1点,直到整个信号被分析。积分功率谱与其阈值的交点为R峰,R峰间的时间间隔即为R-R间期信号。提取出R峰后,将QRS波面积作为特征值来估计EDR信号。通过三次样条插值将R-R间期信号和EDR信号重新采样至4 Hz。最后使用截止频率在0.04 Hz ~ 0.045 Hz 之间的高通滤波器滤除4 Hz R-R间期信号中的低频非周期趋势项,该趋势项可能导致时频分析的失真。
步骤4)对提取出的R-R间期信号和EDR信号逐一进行变分模态分解(VMD)。具体地:VMD非递归地将信号分解为有限个本征模态函数的集合,其中本征模态函数的定义为窄带调幅调频(AM-FM)信号。VMD类比维纳滤波去噪来重构原始信号,计算每个模态的中心频率,将每个模态调制到其中心频率,并通过解决变分问题来最小化各模态的实际带宽之和,使频带分离效果最佳。变分模态分解处理过程中的约束变分问题为
Figure SMS_21
(2)
其中,y为待分解的信号(R-R间期信号或EDR信号);{u k } = {u 1, …,u K },表示分解得到的K个模态的集合;{ω k } = {ω 1 ,…,ω K },表示这K个模态的中心频率,K为参数,由用户自行设置。t代表时间,δ(t)代表狄拉克函数;j为虚数单位。为了使问题无约束,引入二次惩罚项和拉格朗日乘子λ(t),其中二次惩罚项的权重α与噪声等级成反比。最终得到的无约束变分问题是一个增广拉格朗日函数
Figure SMS_22
(3)
利用交替乘子法(ADMM)在傅里叶域内计算这个增广拉格朗日函数关于u k ω k 的鞍点,并在傅里叶域内迭代更新{u k }、{ω k }以及λ得到带宽和最小的解。其中第n+1步的更新公式为
Figure SMS_23
其停止条件为
Figure SMS_24
(7)
其中 k (ω)代表u k (t)的傅里叶变换,(ω)代表y的傅里叶变换,λ̂(ω)代表λ(t)的傅里叶变换,ε是一个预设值,一般设定为10-7
达到停止条件后得到的{u k }即为需要的本征模态函数的集合。
步骤5)经过上一步的运算,对于R-R间期信号和EDR信号,将各自得到K R K E 个模态。为了量化两信号的耦合程度,需要先对R-R间期信号和EDR信号进行希尔伯特变换得到每一个模态的解析函数形式。进一步提取每一个模态对应的包络和相位,其中瞬时包络为
Figure SMS_25
R i (t)代表R-R间期信号,H R,i (t)代表R-R间期信号的希尔伯特变换,E j (t)代表EDR信号,H E,j (t)代表EDR信号的希尔伯特变换。
瞬时相位为
Figure SMS_26
每一个模态对应的瞬时频率为
Figure SMS_27
这里得到的瞬时包络和瞬时相位都是关于时间的一维函数,将瞬时包络、瞬时相位与瞬时频率匹配,以获得瞬时包络和瞬时相位关于时间和频率的二维函数。R-R间期信号和EDR信号的二维包络表示为A R (t, f)和A E (t, f),二维相位表示为θ R (t, f)θ E (t, f)。最终二维解析信号的复指数形式为
Figure SMS_28
步骤6)在评估R-R间期信号与EDR信号之间的耦合强度时,需要考虑两个关键因素。如果在给定频率下,两信号都有比较大的振荡幅度,那么两信号很可能相互耦合,这种效应可以通过计算互功率谱来测量;如果在给定频率上的两个振荡保持恒定的相位关系,那么两信号也很可能相互耦合,这种效应可以通过计算相干性来测量。因此使用互功率谱和相干性的乘积来量化心肺耦合的程度。
互功率谱为
Figure SMS_29
相干性为
Figure SMS_30
心肺耦合指数为
Figure SMS_31
其中<>表示在给定频率下多次测量的平均值,因为相干性是统计结果。32点窗口用于平均相干性。
步骤7)画出心肺耦合图谱,如图1所示,在心肺耦合图谱中能直观地观察不同时段或不同频带的心肺耦合指数大小,颜色越浅代表心肺耦合指数越大,颜色越深代表心肺耦合指数越小。如图2所示,画出受试者的心肺耦合图谱以及他们在实际受试过程中正常呼吸和发生呼吸暂停或低通气事件的时间记录图,其中白色部分代表正常呼吸,灰色部分代表发生呼吸暂停或低通气事件,比较健康受试者和OSAHS患者心肺耦合图谱的不同特征,以及心肺耦合图谱的不同特征与受试者测试过程中真实发生的睡眠事件的对应情况。对已知的受试者的呼吸紊乱指数(AHI)和本发明提出的心肺耦合衍生指标进行Spearman相关性分析。其中,AHI是指平均每小时睡眠中的呼吸暂停加低通气次数,AHI数值越大,代表患病程度越重。心肺耦合衍生指标选择的是归一化高频耦合功率(HFnorm),即高频(0.1Hz ~ 0.4Hz)心肺耦合指数之和与总(0.001 Hz ~ 0.4 Hz)心肺耦合指数之和的比值的时间平均值。
以著名的睡眠数据库Apnea-ECG的研究结果为例,该数据库中心电信号采样率为200Hz,并记录了每位受试者夜间睡眠的呼吸紊乱指数(AHI)和每位受试者发生呼吸暂停或低通气事件的时间。如图2所示,(a)图为健康受试者,(b)、(c)和(d)图为OSAHS患者,在对比健康受试者和OSAHS患者的心肺耦合图谱后发现,健康人的心肺耦合图谱呈现占优势耦合频段自发地在高频耦合(0.1 Hz ~ 0.4 Hz)和低频耦合(0.01 Hz ~ 0.1 Hz)间切换,而对于OSAHS患者,占优势耦合频段大部分时间以低频(0.01 Hz ~ 0.1 Hz)为主,并且患病程度越重,耦合频带越往低频集中。因此,根据心肺耦合图谱呈现出的不同模式可以初步判断受试者是否是OSAHS患者。对比患者的心肺耦合图谱以及他们在实际受试过程中正常呼吸和发生呼吸暂停或低通气事件的时间记录图发现,低频(0.01 Hz ~ 0.1 Hz)心肺耦合指数较小的时刻对应着患者正常呼吸的时刻,低频(0.01 Hz ~ 0.1 Hz)心肺耦合指数较大的时刻对应着患者睡眠过程中发生呼吸暂停或低通气事件的时刻。因此,从心肺耦合图谱能直接估计出患者出现呼吸暂停或低通气事件的时刻。对所有受试者的AHI和HFnorm的Spearman相关性分析结果显示,HFnorm与AHI呈强负相关关系(r= -0.8378,p<0.0001),因此HFnorm可以作为判断患病严重程度的指标。建议基于变分模态分解的心肺耦合分析方法可以作为研究睡眠情况下心肺交互作用,也就是心血管系统与呼吸系统之间内在的协调机制和作用的可靠工具。
综上所述,以上仅为本发明的较佳实施例而已,并非用于限定本发明的保护范围。凡在本发明的精神和原则之内,所作的任何修改、等同替换、改进等,均应包含在本发明的保护范围之内。

Claims (6)

1.一种基于变分心肺耦合的睡眠动力学分析方法,其特征在于该方法的步骤包括:
第一步,采集受试者的心电信号,对采集的心电信号进行预处理,并绘制预处理后心电信号的心电图;
第二步,提取第一步采集到的心电信号的心跳间隔时间序列和心电图衍生的呼吸信号;
第三步,对第二步提取的心跳间隔时间序列和呼吸信号分别进行变分模态分解,得到心跳间隔时间序列的本征模态函数集合和呼吸信号的本征模态函数集合;
第四步,对第三步得到的心跳间隔时间序列的本征模态函数集合进行希尔伯特变换,然后计算心跳间隔时间序列的本征模态函数集合的瞬时包络、瞬时相位和瞬时频率;
第五步,对第三步得到的呼吸信号的本征模态函数集合进行希尔伯特变换,然后计算呼吸信号的本征模态函数集合的瞬时包络、瞬时相位和瞬时频率;
第六步,根据第四步得到的心跳间隔时间序列的本征模态函数集合的瞬时频率将心跳间隔的本征模态函数的时间序列转换为以瞬时频率和时间为因变量的二维解析函数A;
第七步,根据第五步得到的呼吸信号的本征模态函数集合的瞬时频率将呼吸信号的本征模态函数的时间序列转换为以瞬时频率和时间为因变量的二维解析函数B;
第八步,计算第六步得到的二维解析函数A和第七步得到的二维解析函数B的互功率谱和相干性,并进一步得到心肺耦合指数;
第九步,以时间为x轴、频率为y轴、第八步得到的心肺耦合指数为z轴,得到可视化的睡眠期间心肺耦合图谱,并根据得到的心肺耦合图谱得到心肺耦合特征,完成基于变分心肺耦合的睡眠动力学分析;
所述第二步中,提取的心电信号的心跳间隔时间序列为心电图上所有R峰间的时间间隔;
所述R峰的定位方法为心电信号的积分功率谱与心电信号的积分功率阈值的交点;
所述积分功率谱P(t)为:
Figure QLYQS_1
其中,x(t)代表心电信号,g(t-τ)代表高斯函数;
12]为心电信号的频率范围;
f代表频率,t代表时间,τ代表高斯函数g(t)每次平移的步长;
心电信号的积分功率阈值的确定方法为:对积分功率谱进行滑动Tukey窗口加权,取每次加权数据中的最大值作为每个窗口的积分功率阈值;
所述第二步中,提取的心电信号衍生的呼吸信号为心电图上所有QRS波的面积;
所述第三步中,设心跳间隔时间序列的本征模态函数集合为
Figure QLYQS_2
Figure QLYQS_3
KR为心跳间隔时间序列的本征模态函数的个数;本征模态函数集合的确定方法为求出变分模态分解过程中的约束变分问题的解;
所述约束变分问题为:
Figure QLYQS_4
Figure QLYQS_5
其中,yR为心跳间隔时间序列;
Figure QLYQS_6
表示这KR个本征模态函数的中心频率的集合;
t代表时间,δ(t)代表狄拉克函数;j为虚数单位;
所述变分问题的解的确定方法为:
使用二次惩罚项αR和拉格朗日乘子λR(t),将约束变分问题转化为增广拉格朗日函数,所述的增广拉格朗日函数为:
Figure QLYQS_7
其中,αR是一个预设值,λR(t)将在求解过程中变化,初始值为0;
在傅里叶域内计算这个增广拉格朗日函数关于uR,k和ωR,k的鞍点,并在傅里叶域内迭代更新{uR,k}、{ωR,k}以及λR(t)直到满足迭代停止条件,得到本征模态函数集合
Figure QLYQS_8
第n+1步的更新公式为:
Figure QLYQS_9
Figure QLYQS_10
Figure QLYQS_11
停止条件为
Figure QLYQS_12
其中,
Figure QLYQS_13
代表uR,k(t)的傅里叶变换;
Figure QLYQS_14
代表yR的傅里叶变换;
Figure QLYQS_15
代表λR(t)的傅里叶变换;
εR是一个预设值;
达到停止条件后得到的
Figure QLYQS_16
即为心跳间隔时间序列的本征模态函数集合;
所述第三步中,设心电信号衍生的呼吸信号的本征模态函数集合为
Figure QLYQS_17
本征模态函数集合的确定方法为求出变分模态分解过程中的约束变分问题的解;KE为心跳间隔时间序列的本征模态函数的个数;
所述约束变分问题为:
Figure QLYQS_18
Figure QLYQS_19
其中,yE为心电图衍生呼吸信号;
Figure QLYQS_20
表示这KE个本征模态函数的中心频率的集合;
t代表时间,δ(t)代表狄拉克函数;j为虚数单位;
所述变分问题的解的确定方法为:
使用二次惩罚项αE和拉格朗日乘子λE(t),将约束变分问题转化为增广拉格朗日函数,所述的增广拉格朗日函数为:
Figure QLYQS_21
其中,αE是一个预设值,λE(t)将在求解过程中变化,初始值为0;
在傅里叶域内计算这个增广拉格朗日函数关于uE,k和ωE,k的鞍点,并在傅里叶域内迭代更新{uE,k}、{ωE,k}以及λE(t)直到满足迭代停止条件,得到本征模态函数集合
Figure QLYQS_22
第n+1步的更新公式为:
Figure QLYQS_23
Figure QLYQS_24
Figure QLYQS_25
停止条件为
Figure QLYQS_26
其中,
Figure QLYQS_27
代表uE,k(t)的傅里叶变换;
Figure QLYQS_28
代表yE的傅里叶变换;
Figure QLYQS_29
代表λE(t)的傅里叶变换;
εE是一个预设值;
达到停止条件后得到的
Figure QLYQS_30
即为心电信号衍生的呼吸信号的本征模态函数集合。
2.根据权利要求1所述的一种基于变分心肺耦合的睡眠动力学分析方法,其特征在于:
所述第一步中,对采集的心电信号进行预处理时,使用50Hz的陷波滤波器去除心电信号中的工频干扰,使用中值滤波算法去除基线漂移,使用通带截止频率为80Hz、阻带截止频率为100Hz的巴特沃斯低通滤波器去除肌电干扰。
3.根据权利要求1所述的一种基于变分心肺耦合的睡眠动力学分析方法,其特征在于:
所述第四步中,心跳间隔时间序列的本征模态函数集合的瞬时包络为:
Figure QLYQS_31
Ri(t)表示心跳间隔时间序列,HR,i(t)表示心跳间隔时间序列的希尔伯特变换,KR表示分解出的心跳间隔时间序列的本征模态函数集合中本征模态函数的个数;
心跳间隔时间序列的本征模态函数集合的瞬时相位为:
Figure QLYQS_32
心跳间隔时间序列的本征模态函数集合的瞬时频率为:
Figure QLYQS_33
4.根据权利要求3所述的一种基于变分心肺耦合的睡眠动力学分析方法,其特征在于:
所述第五步中,呼吸信号的本征模态函数集合的瞬时包络为:
Figure QLYQS_34
其中,Ej(t)表示呼吸信号,HE,j(t)表示呼吸信号的希尔伯特变换,KE表示分解出的心跳间隔时间序列的本征模态函数集合中本征模态函数的个数;
呼吸信号的本征模态函数集合的瞬时相位为:
Figure QLYQS_35
呼吸信号的本征模态函数集合的瞬时频率为:
Figure QLYQS_36
5.根据权利要求4所述的一种基于变分心肺耦合的睡眠动力学分析方法,其特征在于:
所述第六步中,二维解析函数A为:
Figure QLYQS_37
6.根据权利要求5所述的一种基于变分心肺耦合的睡眠动力学分析方法,其特征在于:
所述第七步中,二维解析函数B为:
Figure QLYQS_38
所述第八步中,互功率谱为:
Figure QLYQS_39
Figure QLYQS_40
相干性为:
Figure QLYQS_41
心肺耦合指数为:
CPC(t,f)=<ΓR,E(t,f)>2ΛR,E(t,f);
其中,<>表示在给定频率下多次测量的平均值。
CN202310227600.0A 2023-03-10 2023-03-10 一种基于变分心肺耦合的睡眠动力学分析方法 Active CN115886742B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202310227600.0A CN115886742B (zh) 2023-03-10 2023-03-10 一种基于变分心肺耦合的睡眠动力学分析方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202310227600.0A CN115886742B (zh) 2023-03-10 2023-03-10 一种基于变分心肺耦合的睡眠动力学分析方法

Publications (2)

Publication Number Publication Date
CN115886742A CN115886742A (zh) 2023-04-04
CN115886742B true CN115886742B (zh) 2023-06-09

Family

ID=86496464

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202310227600.0A Active CN115886742B (zh) 2023-03-10 2023-03-10 一种基于变分心肺耦合的睡眠动力学分析方法

Country Status (1)

Country Link
CN (1) CN115886742B (zh)

Citations (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN105982664A (zh) * 2015-02-26 2016-10-05 张政波 基于单导联ecg的心肺耦合分析方法
CN115630290A (zh) * 2022-11-17 2023-01-20 北京理工大学 基于同步挤压变换的心肺耦合特征提取方法与系统

Family Cites Families (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US7734334B2 (en) * 2004-05-17 2010-06-08 Beth Israel Deaconess Medical Center, Inc. Assessment of sleep quality and sleep disordered breathing based on cardiopulmonary coupling

Patent Citations (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN105982664A (zh) * 2015-02-26 2016-10-05 张政波 基于单导联ecg的心肺耦合分析方法
CN115630290A (zh) * 2022-11-17 2023-01-20 北京理工大学 基于同步挤压变换的心肺耦合特征提取方法与系统

Also Published As

Publication number Publication date
CN115886742A (zh) 2023-04-04

Similar Documents

Publication Publication Date Title
Chua et al. Application of higher order statistics/spectra in biomedical signals—A review
Adnane et al. Development of QRS detection algorithm designed for wearable cardiorespiratory system
O’Brien et al. A comparison of algorithms for estimation of a respiratory signal from the surface electrocardiogram
US7862515B2 (en) Apparatus for detecting sleep apnea using electrocardiogram signals
Kesper et al. ECG signal analysis for the assessment of sleep-disordered breathing and sleep pattern
US6415174B1 (en) ECG derived respiratory rhythms for improved diagnosis of sleep apnea
CN100413461C (zh) 一种获取呼吸暂停事件和睡眠结构图信息的方法
EP1562472A1 (en) Method, apparatus and system for characterizing sleep
JP2009545345A (ja) 心電信号の処理方法、およびその装置
Bricout et al. Accelerometry-derived respiratory index estimating apnea-hypopnea index for sleep apnea screening
CN106901694A (zh) 一种呼吸率提取方法及装置
CN113974576B (zh) 一种基于心磁图的睡眠质量监测系统及监测方法
Kocak et al. Automated detection and classification of sleep apnea types using electrocardiogram (ECG) and electroencephalogram (EEG) features
Sarkar et al. A novel approach towards non-obstructive detection and classification of COPD using ECG derived respiration
Gu et al. An improved method with high anti-interference ability for R peak detection in wearable devices
CN115630290B (zh) 基于同步挤压变换的心肺耦合特征提取方法与系统
Dong et al. An integrated framework for evaluation on typical ECG-derived respiration waveform extraction and respiration
CN115886742B (zh) 一种基于变分心肺耦合的睡眠动力学分析方法
Khavas et al. Robust heartbeat detection using multimodal recordings and ECG quality assessment with signal amplitudes dispersion
Chatlapalli et al. Accurate derivation of heart rate variability signal for detection of sleep disordered breathing in children
Negi et al. Development of a real-time breathing-rate monitor using difference operation method and adaptive windowing on dry-electrode ECG signal
Srinivasa et al. Application of entropy techniques in analyzing heart rate variability using ECG signals
Maier et al. Extraction of respiratory myogram interference from the ECG and its application to characterize sleep-related breathing disorders in atrial fibrillation
Zhao et al. Derivation of respiratory signals from single-lead ECG
Bassiouni et al. Combination of ECG and PPG Signals for Healthcare Applications: A Survey

Legal Events

Date Code Title Description
PB01 Publication
PB01 Publication
SE01 Entry into force of request for substantive examination
SE01 Entry into force of request for substantive examination
GR01 Patent grant
GR01 Patent grant