CN112869733B - 一种心冲击图实时心搏间期测算方法 - Google Patents
一种心冲击图实时心搏间期测算方法 Download PDFInfo
- Publication number
- CN112869733B CN112869733B CN202110025943.XA CN202110025943A CN112869733B CN 112869733 B CN112869733 B CN 112869733B CN 202110025943 A CN202110025943 A CN 202110025943A CN 112869733 B CN112869733 B CN 112869733B
- Authority
- CN
- China
- Prior art keywords
- time
- signal
- interval
- bcg
- peak
- 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
Links
Images
Classifications
-
- A—HUMAN NECESSITIES
- A61—MEDICAL OR VETERINARY SCIENCE; HYGIENE
- A61B—DIAGNOSIS; SURGERY; IDENTIFICATION
- A61B5/00—Measuring for diagnostic purposes; Identification of persons
- A61B5/103—Detecting, measuring or recording devices for testing the shape, pattern, colour, size or movement of the body or parts thereof, for diagnostic purposes
- A61B5/11—Measuring movement of the entire body or parts thereof, e.g. head or hand tremor, mobility of a limb
- A61B5/1102—Ballistocardiography
-
- A—HUMAN NECESSITIES
- A61—MEDICAL OR VETERINARY SCIENCE; HYGIENE
- A61B—DIAGNOSIS; SURGERY; IDENTIFICATION
- A61B5/00—Measuring for diagnostic purposes; Identification of persons
- A61B5/02—Detecting, measuring or recording pulse, heart rate, blood pressure or blood flow; Combined pulse/heart-rate/blood pressure determination; Evaluating a cardiovascular condition not otherwise provided for, e.g. using combinations of techniques provided for in this group with electrocardiography or electroauscultation; Heart catheters for measuring blood pressure
- A61B5/024—Detecting, measuring or recording pulse rate or heart rate
- A61B5/02405—Determining heart rate variability
-
- A—HUMAN NECESSITIES
- A61—MEDICAL OR VETERINARY SCIENCE; HYGIENE
- A61B—DIAGNOSIS; SURGERY; IDENTIFICATION
- A61B5/00—Measuring for diagnostic purposes; Identification of persons
- A61B5/72—Signal processing specially adapted for physiological signals or for diagnostic purposes
- A61B5/7203—Signal processing specially adapted for physiological signals or for diagnostic purposes for noise prevention, reduction or removal
- A61B5/7207—Signal processing specially adapted for physiological signals or for diagnostic purposes for noise prevention, reduction or removal of noise induced by motion artifacts
-
- A—HUMAN NECESSITIES
- A61—MEDICAL OR VETERINARY SCIENCE; HYGIENE
- A61B—DIAGNOSIS; SURGERY; IDENTIFICATION
- A61B5/00—Measuring for diagnostic purposes; Identification of persons
- A61B5/72—Signal processing specially adapted for physiological signals or for diagnostic purposes
- A61B5/7225—Details of analog processing, e.g. isolation amplifier, gain or sensitivity adjustment, filtering, baseline or drift compensation
-
- A—HUMAN NECESSITIES
- A61—MEDICAL OR VETERINARY SCIENCE; HYGIENE
- A61B—DIAGNOSIS; SURGERY; IDENTIFICATION
- A61B2576/00—Medical imaging apparatus involving image processing or analysis
- A61B2576/02—Medical imaging apparatus involving image processing or analysis specially adapted for a particular organ or body part
- A61B2576/023—Medical imaging apparatus involving image processing or analysis specially adapted for a particular organ or body part for the heart
Landscapes
- Health & Medical Sciences (AREA)
- Life Sciences & Earth Sciences (AREA)
- Engineering & Computer Science (AREA)
- Biomedical Technology (AREA)
- Medical Informatics (AREA)
- Cardiology (AREA)
- Physiology (AREA)
- Veterinary Medicine (AREA)
- Physics & Mathematics (AREA)
- Signal Processing (AREA)
- Biophysics (AREA)
- Pathology (AREA)
- Public Health (AREA)
- Heart & Thoracic Surgery (AREA)
- General Health & Medical Sciences (AREA)
- Molecular Biology (AREA)
- Surgery (AREA)
- Animal Behavior & Ethology (AREA)
- Artificial Intelligence (AREA)
- Psychiatry (AREA)
- Computer Vision & Pattern Recognition (AREA)
- Power Engineering (AREA)
- Dentistry (AREA)
- Oral & Maxillofacial Surgery (AREA)
- Measurement And Recording Of Electrical Phenomena And Electrical Characteristics Of The Living Body (AREA)
Abstract
本发明提供一种心冲击图实时心搏间期测算方法,包括如下步骤:S1、进行BCG信号采集;S2、对采集的BCG信号进行预处理,包括滤波去热噪声及基线处理、形态滤波去呼吸干扰处理;S3、实时定位每一个BCG信号J峰,测量实时JJ间期;S4、计算瞬时心率、多时间尺度平均心率;S5、显示BCG波形、瞬时心率、多时间尺度平均心率;解决了传统BCG实时心搏检测上限范围过低(低于120次/min)的问题,还解决短时体动干扰(<4s)连续心搏间期检测误差过大的问题。
Description
技术领域
本发明涉及医疗测算技术领域,具体涉及一种心冲击图实时心搏间期测算方法。
背景技术
心冲击图描记了心脏活动及血液运行时,心脏和大血管冲击和反弹而引起躯体微弱运动。心冲击图曲线形态丰富,包含了H、I、J、K、L共5个波段。得益于传感器技术迅猛发展,心冲击图的获取可以在用户处于静息情况下被相对准确的采集。鉴于心冲击图与心脏收缩力、心排血量、血流速度、大血管弹性等有关,对心冲击图信号的采集与定位,对进一步分析心率变异性以及心功能检查等具有一定意义。
参阅附图1BCG信号时域波形图,由图1可知,一个完整的BCG信号包括除主峰J之外的H、I、K、L等次峰,其中,H、L次峰方向于主峰J峰一致,在不同心率的情况下,H-J、J-L间期时间大致范围在150ms-350ms,这意味着,在心率>120次/min的情况下,当前BCG信号的H峰、L峰可能分别与上下相邻BCG信号的L峰、H峰重叠,重叠后信号可能产生高于J峰的形态,给实时J峰定位带来极大挑战。
现有技术中,存在两点问题,问题一:传统瞬时心率测算依据为识别每一个BCG信号的J峰,并借助相邻BCG信号的J-J间期计算当前心搏的瞬时心率。问题二:传统BCG信号定位技术方案的大致思路为:以当前BCG信号(第i个)J峰时刻为基准,前向固定时间尺度遍历,搜索下一个BCG信号(第i+1)J峰。以文献常用的搜索时间尺度0.5s—1.5s为例,最高适用于心搏120次/min的场景,当心搏速率高于120次/min时,上述前向遍历区间内不仅忽略了相邻下一BCG信号(第i+1)J峰,同时却覆盖了最多高达3个非相邻BCG-J峰(第i+2、i+3甚至i+4个BCG信号),导致无法准确测量心率>120次/min的相邻J-J间期。
发明内容
为解决上述问题,本发明提供了一种心冲击图实时心搏间期测算方法,包括如下步骤:
S1、进行BCG信号采集;
S2、对采集的BCG信号进行预处理,包括滤波去热噪声及基线处理、形态滤波去呼吸干扰处理;
S3、实时定位每一个BCG信号J峰,测量实时JJ间期,实现两个目标:
目标一:心搏测量上限范围有传统方法上限的120次/min提升至150次/min;
目标二:对<4s短时体动区间BCG信号JJ间期预测补偿;
S4、计算瞬时心率、多时间尺度平均心率;
S5、显示BCG波形、瞬时心率、多时间尺度平均心率。
进一步地,所述步骤S2中,滤波去热噪声及基线处理具体为:采用50Hz陷波器,去除信号50Hz工频干扰;设计截止频率为0.01Hz至15Hz的带通滤波器,去除基线及15Hz以上的热噪声干扰。
进一步地,所述步骤S2中,形态滤波去呼吸干扰处理包括如下步骤:S21、设计截止频率为20Hz的LPF滤波器;
S22、对滤波后信号进行开运算;
S23、对S22中输出信号进行闭运算,输出为原始采集信号中的基线漂移信号;
S24、从原始信号中移除第三步输出的基线漂移信号,分离呼吸干扰,获得去干扰后的BCG信号;
S25、原始信号中标识出体动间期,具体为:计算缓存信号幅值分布,借助信号分布函数中最大、最小与主要分布均值,排除因临床导致信号过小、体动导致信号过大两种状态。
进一步地,所述步骤S3中,所述目标一的实现步骤是:
首先定义信号处理堆栈时间窗长度为T=10秒,信号实时进入时间窗,更新时间频率为1s,直至信号填满堆栈开始运算;
S311、检测堆栈内信号是否满足计算条件;
S312、堆栈内BCG信号模板初始化;
S313、确定基于初始J峰定位的动态搜索区间;
S314、搜索动态区间内所有BCG信号。
进一步地,所述步骤S311具体为:
设定BCG幅值上限为A_max,当堆栈T内信号超出A_max对应时长>4s,则定义超出A_max对应连续时间尺度内的信号区间为体动区间,并定义该区间为不应期,不参与BCG信号运算;
如超出A_max对应信号时间不足4s,对相应区间BCG间期测算采用自回归模型填充;
所述步骤S312具体为:
信号采集过程,当连续时间尺度T’内采集信号均满足峰值均小于A_max的条件,视为信号初始化运算的基本条件;
对T’区间信号进行固定子时间窗切割,相邻子窗交叠50%,子时间窗初始化间隔为1s;
定位每个1s子窗定没所有上峰和下峰,当相邻上峰及下峰满足时间间隔大于100ms且小于350ms的约束条件;
计算每一个相邻上峰至下峰构成直角三角形区域面积,选取每个1s子窗内最大三角形面积对应上峰时间;
构造出T’时间间期内形态距与对应时间的2*N维矩阵A,N为依托子窗定位T’时间间期内相邻上下峰构成最大面积直角三角形的个数;其中,矩阵第一行A[1,:]为面积,第二行A[2,:]为对应上峰时间;
计算A[1,:]构成的1*N维矢量值:D_m=mean{A[1,:]},计算A[2,:]构成1*N维矢量中的间期均值:T_m=mean{A[2,:]},选定矩阵A中同时不满足以下两个条件的所有样本并剔除:
条件(1):A[1,:]中样本应满足[a1*D_m,a2*D_m],其中a1为小于1的非负数,a2为大于1的非负数;
条件(2):A[2,:]中样本应满足[b1*D_m,b2*D_m],其中b1为小于1的非负数,b2为大于1的非负数;
选取矩阵A中同时满足满足非上述条件的第一个三角形上峰对应时间t_0作为初始BCG信号J峰的时间,t_0时刻的幅值为初始BCG信号J峰的幅值,定义初始J峰为BCG_J0。
当矩阵A同时不满足条件(1)和条件(2)的样本多于30%,则认为当前T’时间尺度的所有样本可信度不满足检测初始化要求,丢弃所有数据,并对下一T’时间尺度信号重复上述检测,直至满足要求;所述步骤S313具体为:
以初始化BCG信号J峰t_0时间起点,测算T’时间尺度内初始化心搏间期平均值Tb_mean,前向遍历一个动态时间区间内所有信号,其中,动态时间取值范围取决于Beat范围,并做如下操作定义:
当Beat范围>150次/min,重复步骤S312计算过程;
当Beat范围为(100,150]次/min,动态时间区间范围:[0.375,1.2];
当Beat范围为(50,99]次/min,动态时间区间范围:[0.5,1.5];
当Beat范围为(0,50]次/min,重复步骤S312计算过程;
所述步骤S314具体为:
当平均心率达[100,150]次/min范围时,最小相邻心搏JJ间期为0.4s,最大相邻心搏JJ间期为0.6s,此时基于初始化BCG_J0的前向时间尺度范围为[0.375,1.2]秒,判定为:在无心律不齐的正常窦性心律情况下,该时间尺度内存在1-3个BCG信号;
对于遍历区间真实存在BCG信号情况下的检测步骤如下:
遍历对应时间尺度范围内所有上下峰,计算每一个相邻上下峰构成的形态距,其中计算纳入标准为:相邻上下峰对应时间间隔大于100ms且小于250ms;
定义所有满足纳入标准的上下峰形态距集合为矢量D=[d1,…,dK],dk为第k个满足纳入标准的形态距,K>=1为时间尺度内所有纳入形态距的个数;
将D中所有样本与已检测到最近的M个BCG信号J-K峰形态距的均值d_mean进行对比,选出矢量D中所有[0.707*d_mean,1.414*d_mean]范围内形态距对应的峰值时间;
计算满足上述形态距条件峰值时间与相邻已检出BCG信号J峰,构成若干相邻JJ间期T_j=[tj_1,…tj_N],N为所有JJ间期的个数;将T_j矢量中所有样本与已检测到最近的M个BCG信号的JJ间期均值tj_mean进行对比,选出以最临近减除BCG_J0为基准时间,区间[0.85*tj_mean,1.15*tJ_mean]内最大形态距样本对应上峰,定义为前向遍历区间内的第一个检出BCG信号J峰,定义为BCG_J1;
以第一个检出BCG信号J峰对应时间替上一个BCG信号J峰作为初始化时间起点,重复步骤S314,直至检出新一个BCG信号J峰。
进一步地,当遍历区间不满足步骤S314中确定BCG信号J峰的检测标准时,对前向遍历区间信号开展有效性分析,包括:
遍历区间内信号的平均功率与已检出BCG信号的同样时间尺度的8个区间进行对比;当遍历区间信号平均功率低于之前8个同样尺度区间平均功率的下限时,原因为信号输出不满足检测标准,显示单元输出“信号无效”,将当前遍历时间尺度定义为信号不应期,并重复步骤S312;
当遍历区间满足检测标准,可能由于残留的身体微震干扰或残留呼吸干扰导致BCG信号无法检出,此时,根据最近已检出的连续M个BCG信号的JJ间期平均值,在遍历区间中赋值,定义赋值时间为“虚拟”BCG信号J峰时间,计算对应“虚拟”JJ间期,并重复步骤S314进行下一个遍历区间的搜索。
进一步地,所述步骤S3中,所述目标二的实现步骤是:
S321:基于双向检测的体动尾迹干扰时长的确定;
S322:补偿体动尾迹期间JJ间期;
S323:<4s不应期的JJ间期补偿及相邻T’时间尺度的信息应予以指针标注并缓存,为系统非实时二次重定位分析提供时间索引。
本发明还提供了一种心冲击图实时采集系统,使用权利要求1至7任一项所述的心搏间期测算方法,包括BCG信号预处理单元、BCG信号实时心搏间期测量单元和心率波形显示单元,所述BCG信号预处理单元与所述BCG信号实时心搏间期测量单元相连,所述BCG信号实时心搏间期测量单元与所述心率波形显示单元相连,所述BCG信号预处理单元包括陷波器、带通滤波器、LPF滤波器,所述陷波器设为最先接受信号,所述带通滤波器设为第二接受信号,所述LPF滤波器设为第三接受信号,所述心率波形显示单元与处理器连接,所述心率波形显示单元包括独立显示模组、终端电子设备。
进一步地,所述陷波器为50Hz陷波器,所述带通滤波器为0.01Hz至15Hz的带通滤波器,所述LPF滤波器为截止频率为20Hz的LPF滤波器。
进一步地,所述心率波形显示单元可以通过数据传输线与处理器连接,也可以通过无线方式与处理器连接,所述心率波形显示单元的终端电子设备包括手机、PAD、个人电脑。
本发明的心冲击图实时心搏间期测算方法,采集BCG信号后先对进行预处理,去除噪声、基线处理及去呼吸干扰,随后实时定位每一个BCG信号J峰,测量实施JJ间期,以实现心搏测量上限范围有传统方法上限的120次/min提升至150次/min、实现对<3s的小体动区间进行心搏间期预测补偿,提升BCG检测覆盖率。解决了传统BCG实时心搏检测上限范围过低(低于120次/min)的问题,还解决短时体动干扰(<4s)连续心搏间期检测误差过大的问题。
附图说明
图1为BCG信号时域波形图;
图2为本发明方法的BCG信号采集系统功能框图;
图3为本发明BCG信号预处理中形态学滤波器算法流程图;
图4为本发明目标2中混叠<4s体动残留干扰的BCG信号图。
具体实施方式
下面将结合本发明实施例中的附图,对本发明实施例中的技术方案进行清楚、完整的描述,显然,所描述的实施例仅仅是本发明一部分实施例,而不是全部的实施例。本领域普通人员在没有做出创造性劳动前提下所获得的所有其他实施例,均属于本发明的保护范围。
本发明提供了一种心冲击图实时心搏间期测算方法,旨在解决传统BCG实时心搏检测上限范围过低的问题,同时还为解决短时体动干扰(<4s)连续心搏间期检测误差过大的问题作出了具体设计,本发明的步骤如下:
S1、进行BCG信号采集;
S2、对采集的BCG信号进行预处理,包括滤波去热噪声及基线处理、形态滤波去呼吸干扰处理;
S3、实时定位每一个BCG信号J峰,测量实时JJ间期,实现两个目标:
目标一:心搏测量上限范围有传统方法上限的120次/min提升至150次/min;
目标二:对<4s短时体动区间BCG信号JJ间期预测补偿;
S4、计算瞬时心率、多时间尺度平均心率;
S5、显示BCG波形、瞬时心率、多时间尺度平均心率。
步骤S2中,滤波去热噪声及基线处理具体为:采用50Hz陷波器,去除信号50Hz工频干扰;设计截止频率为0.01Hz至15Hz的带通滤波器,去除基线及15Hz以上的热噪声干扰。
步骤S2中,形态滤波去呼吸干扰处理:BCG信号源自心冲击力,因而采集过程与呼吸引起的胸腹运动力在时间维度呈叠加态,形态滤波去呼吸干扰处理包括如下步骤:
S21、设计截止频率为20Hz的LPF滤波器;可以为多阶FIR滤波器、巴特沃斯滤波器等。
S22、对滤波后信号进行开运算;具体使用先腐蚀再膨胀的方式;
S23、对S22中输出信号进行闭运算,使用先膨胀再腐蚀的方式,输出为原始采集信号中的基线漂移信号;
S24、从原始信号中移除第三步输出的基线漂移信号,分离呼吸干扰,获得去干扰后的BCG信号。
前述形态滤波整体设计流程参阅附图3,本发明BCG信号预处理中形态学滤波器算法流程图。
S25、原始信号中标识出体动间期,具体为:计算缓存信号幅值分布,借助信号分布函数中最大、最小与主要分布均值,排除因临床导致信号过小、体动导致信号过大两种状态。
步骤S3中,所述目标一的实现步骤是:
首先定义信号处理堆栈时间窗长度为T=10秒,信号实时进入时间窗,更新时间频率为1s,直至信号填满堆栈开始运算。
S311、检测堆栈内信号是否满足计算条件;
S312、堆栈内BCG信号模板初始化;
S313、确定基于初始J峰定位的动态搜索区间;
S314、搜索动态区间内所有BCG信号。
步骤S3中,所述目标一的实现步骤是:
首先定义信号处理堆栈时间窗长度为T=10秒,信号实时进入时间窗,更新时间频率为1s,直至信号填满堆栈开始运算。
S311、检测堆栈内信号是否满足计算条件;
S312、堆栈内BCG信号模板初始化;
S313、确定基于初始J峰定位的动态搜索区间;
S314、搜索动态区间内所有BCG信号。
步骤S311具体为:
设定BCG幅值上限为A_max(单位:mV),当堆栈T内信号超出A_max对应时长>4s(占时间窗40%及以上),则定义超出A_max对应连续时间尺度内的信号区间为体动区间,并定义该区间为不应期,不参与BCG信号运算。
如超出A_max对应信号时间不足4s,对相应区间BCG间期测算采用自回归模型填充。
步骤S312具体为:
信号采集过程,当连续时间尺度T’(例如T’<10s)内采集信号均满足峰值均小于A_max的条件,视为信号初始化运算的基本条件。对T’区间信号进行固定子时间窗切割,相邻子窗交叠50%,子时间窗初始化间隔为1s。
定位每个1s子窗定没所有上峰和下峰,当相邻上峰及下峰满足时间间隔大于100ms且小于350ms的约束条件。
计算每一个相邻上峰至下峰构成直角三角形区域面积,选取每个1s子窗内最大三角形面积对应上峰时间。
构造出T’时间间期内形态距与对应时间的2*N维矩阵A,N为依托子窗定位T’时间间期内相邻上下峰构成最大面积直角三角形的个数。注意:矩阵第一行A[1,:]为面积,第二行A[2,:]为对应上峰时间。
计算A[1,:]构成的1*N维矢量值:D_m=mean{A[1,:]},计算A[2,:]构成1*N维矢量中的间期均值:T_m=mean{A[2,:]},选定矩阵A中同时不满足以下两个条件的所有样本并剔除:
条件(1):A[1,:]中样本应满足[a1*D_m,a2*D_m],其中a1为小于1的非负数,a2为大于1的非负数(例如a1=0.707,a2=1.414)。
条件(2):A[2,:]中样本应满足[b1*D_m,b2*D_m],其中b1为小于1的非负数,b2为大于1的非负数(例如b1=0.707,b2=1.414)。选取矩阵A中同时满足满足非上述条件的第一个三角形上峰对应时间t_0作为初始BCG信号J峰的时间,t_0时刻的幅值为初始BCG信号J峰的幅值,定义初始J峰为BCG_J0。
当矩阵A同时不满足条件(1)和条件(2)的样本多于30%,则认为当前T’时间尺度的所有样本可信度不满足检测初始化要求,丢弃所有数据,并对下一T’时间尺度信号重复上述检测,直至满足要求。
步骤S313具体为:
以初始化BCG信号J峰t_0时间起点,测算T’时间尺度内初始化心搏间期平均值Tb_mean(T’时间间隔的平均心率Beat=60/Tb_mean后四舍五入取整数),前向遍历一个动态时间区间内所有信号,其中,动态时间取值范围取决于Beat范围,并做如下操作定义:
当Beat范围>150次/min,重复步骤2计算过程;
当Beat范围为(100,150]次/min,动态时间区间范围:[0.375,1.2](单位:秒);
当Beat范围为(50,99]次/min,动态时间区间范围:[0.5,1.5](单位:秒);
当Beat范围为(0,50]次/min,重复步骤2计算过程。
步骤S314具体为:
当平均心率达[100,150]次/min范围时,最小相邻心搏JJ间期为0.4s,最大相邻心搏JJ间期为0.6s,此时基于初始化BCG_J0的前向时间尺度范围为[0.375,1.2]秒,可知,在无心律不齐的正常窦性心律情况下,该时间尺度内存在1-3个BCG信号。
对于遍历区间真实存在BCG信号情况下的检测方案如下:
遍历对应时间尺度范围内所有上下峰,计算每一个相邻上下峰构成的形态距,计算纳入标准为:相邻上下峰对应时间间隔大于100ms且小于250ms。
定义所有满足纳入标准的上下峰形态距集合为矢量D=[d1,…,dK],dk为第k个满足纳入标准的形态距,K>=1为时间尺度内所有纳入形态距的个数。
将D中所有样本与已检测到最近的M(M为大于1的正整数,例如M=8)个BCG信号J-K峰形态距的均值d_mean进行对比,选出矢量D中所有[0.707*d_mean,1.414*d_mean]范围内形态距对应的峰值时间。
计算满足上述形态距条件峰值时间与相邻已检出BCG信号J峰,构成若干相邻JJ间期T_j=[tj_1,…tj_N],N为所有JJ间期的个数。
将T_j矢量中所有样本与已检测到最近的M个BCG信号的JJ间期均值tj_mean进行对比,选出以最临近减除BCG_J0为基准时间,区间[0.85*tj_mean,1.15*tJ_mean]内最大形态距样本对应上峰,定义为前向遍历区间内的第一个检出BCG信号J峰,定义为BCG_J1。
以第一个检出BCG信号J峰对应时间替上一个BCG信号J峰作为初始化时间起点,重复步骤S314,直至检出新一个BCG信号J峰。
当遍历区间不满足步骤4中确定BCG信号J峰的检测标准时,对前向遍历区间信号开展有效性分析,包括:
遍历区间内信号的平均功率与已检出BCG信号的同样时间尺度的8个区间(相邻区间50%交叠)进行对比;当遍历区间信号平均功率低于之前8个同样尺度区间平均功率的下限(经验阈值,根据已检区间平均功率设定)时,原因为信号输出不满足检测标准,显示单元输出“信号无效”,将当前遍历时间尺度定义为信号不应期,并重复步骤S312。
当遍历区间满足检测标准,可能由于残留的身体微震干扰或残留呼吸干扰导致BCG信号无法检出,此时,根据最近已检出的连续M个BCG信号的JJ间期平均值,在遍历区间中赋值,定义赋值时间为“虚拟”BCG信号J峰时间,计算对应“虚拟”JJ间期,并重复步骤4进行下一个遍历区间的搜索。
附图4为本发明目标2中混叠<4s体动残留干扰的BCG信号图,即为存在<4s体动片段的案例,其中,体动时间间隔<4s(身体微小振颤、吞咽、深呼吸等动作引起),发生于固定时间尺度(图4案例中尺度为10s)的中间时刻,持续时长约为2.1-2.3s,其中涵盖2个BCG信号。在无体动期,BCG信号J峰均已四角形标注。针对上述体动尾迹干扰的补偿,步骤S3中,所述目标二的实现步骤是:
S321:基于双向检测的体动尾迹干扰时长的确定;由于在前向体动尾迹区间搜索BCG信号J峰,不满足上下峰形态距约束条件,且同时不满足JJ间期时间约束条件。因此,上述体动尾迹期间被视为不应期。
定义不应期相邻前后BCG信号J峰间期与不应期相邻前后两个JJ间期均值之差为实际体动尾迹持续时长。
对不应期时长>4s的情况,不进行关于BCG信号的任何相关运算。在不应期结束后,初始化心率为基于不应期前T’时间的平均心率,并执行S312-S314初始化BCG信号J峰及相关检测。
S322:补偿体动尾迹期间JJ间期;有两种方法补偿体动尾迹期间JJ间期。
方法一:体动前Q个BCG信号JJ间期与体动后Q个JJ间期(Q为大于等于1的正整数),对体动区间进行三次样条差值。
方法二:从前T’时间尺度中,选出满足[0.92*t_mean,1.16*t_mean]范围内的M(M>2)个BCG信号JJ间期,采用M阶自回归模型进行迭代补偿,其中自回归模型遗忘因子为小于1的非负数。
借助方法1和方法2,构造补偿体动伪迹对应区间的JJ间期序列,分别为T_J=[TJ1,TJ2,…,TJN]和T_J’=[TJ1,TJ2,…,TJN]其中,TJi为堆栈第i和i+1个BCG信号J-J峰间期。计算T_J和T_J’两个序列时间尺度的绝对离差,如绝对离差大于阈值c(c为>0的实数),表明单一方法补偿数据不稳健,伪迹期间间期补偿无效,取均值填充;当对离差小于阈值,由两种方法补偿间期的均值,即0.5*(T_J+T_J’)替代。
S323:<4s不应期的JJ间期补偿及相邻T’时间尺度的信息应予以指针标注并缓存,为系统非实时二次重定位分析提供时间索引。
结合参阅附图2本发明的方法的BCG信号采集系统功能框图,本发明还提供了一种心冲击图实时采集系统,包括BCG信号预处理单元、BCG信号实时心搏间期测量单元和心率波形显示单元,BCG信号预处理单元与BCG信号实时心搏间期测量单元相连,BCG信号实时心搏间期测量单元与心率波形显示单元相连,BCG信号预处理单元包括陷波器、带通滤波器、LPF滤波器,BCG信号处理单元具备缓存设置,对信号检测过程中所有不应期及其相邻时间尺度信号进行标注及缓存,具备二次离线(非实时)分析能力。
陷波器设为最先接受信号,所述带通滤波器设为第二接受信号,所述LPF滤波器设为第三接受信号,所述心率波形显示单元与处理器连接,心率波形显示单元包括独立显示模组、终端电子设备。
其中,陷波器为50Hz陷波器,带通滤波器为0.01Hz至15Hz的带通滤波器,LPF滤波器为截止频率为20Hz的LPF滤波器。
心率波形显示单元可以通过数据传输线与处理器连接,也可以通过无线方式与处理器连接,所述心率波形显示单元的终端电子设备包括手机、PAD、个人电脑等。
BCG信号实时心搏间期测量单元实时定位每一个BCG信号J峰,测量实时JJ间期,实现前述步骤3中的两个目标。
心率波形显示单元包括以下三个功能:
功能1:可以实时显示BCG波形(滞后时间T’);
功能2:可以显示瞬时心率(根据上一检出BCG与当前检出BCG信号的JJ间期测算)。
功能3:可以显示多时间尺度的平均心率(根据c*T’时间尺度的平均JJ间期测算,c为大于等于1的正整数)。
此外,显示单元同时具备以下内容显示功能:
1.当平均心率>100次/min,显示单元同步输出“心动过速”或相近内容文字信息,同时显示Beat平均心率数值。
2.当平均心率<50次/min,显示单元同步输出“心动过缓”或相近内容文字信息,同时显示Beat平均心率数值。
3.当信号检测不满足条件或处于>3s不应期,显示单元同步输出“信号无效”或相近内容文字信息,不显示心率数值。
以上所述为本发明的较佳实施例而已,但本发明不应局限于该实施例和附图所公开的内容,所以凡是不脱离本发明所公开的精神下完成的等效或修改,都落入本发明保护的范围。
Claims (6)
1.一种心冲击图实时心搏间期测算方法,其特征在于,具体包括如下步骤:
S1、进行BCG信号采集;
S2、对采集的BCG信号进行预处理,包括滤波去热噪声及基线处理、形态滤波去呼吸干扰处理;
S3、实时定位每一个BCG信号J峰,测量实时JJ间期,实现两个目标:
目标一:心搏测量上限范围由传统方法上限的120次/min提升至150次/min;
所述目标一的实现步骤是:
首先定义信号处理堆栈时间窗长度为T=10秒,信号实时进入时间窗,更新时间频率为1s,直至信号填满堆栈开始运算;
S311、检测堆栈内信号是否满足计算条件;
S312、堆栈内BCG信号模板初始化;
S313、确定基于初始J峰定位的动态搜索区间;
S314、搜索动态区间内所有BCG信号;
所述步骤S311具体为:
设定BCG幅值上限为A_max,当堆栈T内信号超出A_max对应时长>4s,则定义超出A_max对应连续时间尺度内的信号区间为体动区间,并定义该区间为不应期,不参与BCG信号运算;
如超出A_max对应信号时间不足4s,对相应区间BCG间期测算采用自回归模型填充;
所述步骤S312具体为:
信号采集过程,当连续时间尺度T’内采集信号均满足峰值均小于A_max的条件,视为信号初始化运算的基本条件;
对T’区间信号进行固定子时间窗切割,相邻子窗交叠50%,子时间窗初始化间隔为1s;
定位每个1s子窗定位所有上峰和下峰,当相邻上峰及下峰满足时间间隔大于100ms且小于350ms的约束条件;
计算每一个相邻上峰至下峰构成直角三角形区域面积,选取每个1s子窗内最大三角形面积对应上峰时间;
构造出T’时间间期内形态距与对应时间的2*N维矩阵A,N为依托子窗定位T’时间间期内相邻上下峰构成最大面积直角三角形的个数;其中,矩阵第一行A[1, :]为面积,第二行A[2, :]为对应上峰时间;
计算A[1, :]构成的1*N维矢量值:D_m = mean{A[1, :]},计算A[2, :]构成1*N维矢量中的间期均值:T_m = mean{A[2, :]},选定矩阵A中同时不满足以下两个条件的所有样本并剔除:
条件(1):A[1, :]中样本应满足[a1*D_m, a2*D_m],其中a1为小于1的非负数,a2为大于1的非负数;
条件(2):A[2, :]中样本应满足[b1*D_m, b2*D_m],其中b1为小于1的非负数,b2为大于1的非负数;
选取矩阵A中同时满足非上述条件的第一个三角形上峰对应时间t_0作为初始BCG信号J峰的时间,t_0时刻的幅值为初始BCG信号J峰的幅值,定义初始J峰为BCG_J0;
当矩阵A同时不满足条件(1)和条件(2)的样本多于30%,则认为当前T’时间尺度的所有样本可信度不满足检测初始化要求,丢弃所有数据,并对下一T’时间尺度信号重复上述检测,直至满足要求;
所述步骤S313具体为:
以初始化BCG信号J峰t_0时间起点,测算T’时间尺度内初始化心搏间期平均值Tb_mean,前向遍历一个动态时间区间内所有信号,其中,动态时间取值范围取决于Beat范围,并做如下操作定义:
当Beat范围>150次/min,重复步骤S312计算过程;
当Beat范围为(100, 150]次/min,动态时间区间范围:[0.375, 1.2];
当Beat 范围为(50, 99]次/min,动态时间区间范围:[0.5, 1.5];
当Beat范围为(0, 50]次/min,重复步骤S312计算过程;
所述步骤S314具体为:
当平均心率达[100, 150]次/min范围时,最小相邻心搏JJ间期为0.4s,最大相邻心搏JJ间期为0.6s,此时基于初始化BCG_J0的前向时间尺度范围为[0.375, 1.2]秒,判定为:在无心律不齐的正常窦性心律情况下,该时间尺度内存在1-3个BCG信号;
对于遍历区间真实存在BCG信号情况下的检测步骤如下:
遍历对应时间尺度范围内所有上下峰,计算每一个相邻上下峰构成的形态距,其中计算纳入标准为:相邻上下峰对应时间间隔大于100ms且小于250ms;
定义所有满足纳入标准的上下峰形态距集合为矢量D = [d1,…, dK],dk为第k个满足纳入标准的形态距,K >= 1为时间尺度内所有纳入形态距的个数;
将D中所有样本与已检测到最近的M个BCG信号J-K峰形态距的均值d_mean进行对比,选出矢量D中所有[0.707*d_mean, 1.414*d_mean]范围内形态距对应的峰值时间;
计算满足上述形态距条件峰值时间与相邻已检出BCG信号J峰,构成若干相邻JJ间期T_j = [tj_1, …tj_N],N为所有JJ间期的个数;
将T_j矢量中所有样本与已检测到最近的M个BCG信号的JJ间期均值tj_mean进行对比,选出以最临近减除BCG_J0为基准时间,区间[0.85*tj_mean, 1.15*tj_mean]内最大形态距样本对应上峰,定义为前向遍历区间内的第一个检出BCG信号J峰,定义为BCG_J1;
以第一个检出BCG信号J峰对应时间替上一个BCG信号J峰作为初始化时间起点,重复步骤S314,直至检出新一个BCG信号J峰;
当遍历区间不满足步骤S314中确定BCG信号J峰的检测标准时,对前向遍历区间信号开展有效性分析,包括:
遍历区间内信号的平均功率与已检出BCG信号的同样时间尺度的8个区间进行对比;当遍历区间信号平均功率低于之前8个同样尺度区间平均功率的下限时,原因为信号输出不满足检测标准,显示单元输出“信号无效”,将当前遍历时间尺度定义为信号不应期,并重复步骤S312;
当遍历区间满足检测标准,可能由于残留的身体微震干扰或残留呼吸干扰导致BCG信号无法检出,此时,根据最近已检出的连续M个BCG信号的JJ间期平均值,在遍历区间中赋值,定义赋值时间为“虚拟”BCG信号J峰时间,计算对应“虚拟”JJ间期,并重复步骤S314进行下一个遍历区间的搜索;
目标二:对<4s短时体动区间BCG信号JJ间期预测补偿;
所述目标二的实现步骤是:
S321:基于双向检测的体动尾迹干扰时长的确定;
S322:补偿体动尾迹期间JJ间期;
S323:<4s不应期的JJ间期补偿及相邻T’时间尺度的信息应予以指针标注并缓存,为系统非实时二次重定位分析提供时间索引;
S4、计算瞬时心率、多时间尺度平均心率;
S5、显示BCG波形、瞬时心率、多时间尺度平均心率。
2.根据权利要求1所述的心冲击图实时心搏间期测算方法,其特征在于:所述步骤S2中,滤波去热噪声及基线处理具体为:采用50Hz陷波器,去除信号50Hz工频干扰;设计截止频率为0.01Hz至15Hz的带通滤波器,去除基线及15Hz以上的热噪声干扰。
3.根据权利要求1所述的心冲击图实时心搏间期测算方法,其特征在于:所述步骤S2中,形态滤波去呼吸干扰处理包括如下步骤:
S21、设计截止频率为20Hz的LPF滤波器;
S22、对滤波后信号进行开运算;
S23、对S22中输出信号进行闭运算,输出为原始采集信号中的基线漂移信号;
S24、从原始信号中移除第三步输出的基线漂移信号,分离呼吸干扰,获得去干扰后的BCG信号;
S25、原始信号中标识出体动间期,具体为:计算缓存信号幅值分布,借助信号分布函数中最大、最小与主要分布均值,排除因临床导致信号过小、体动导致信号过大两种状态。
4.一种心冲击图实时采集系统,使用权利要求1至3任一项所述的心冲击图实时心搏间期测算方法,其特征在于,包括BCG信号预处理单元、BCG信号实时心搏间期测量单元和心率波形显示单元,所述BCG信号预处理单元与所述BCG信号实时心搏间期测量单元相连,所述BCG信号实时心搏间期测量单元与所述心率波形显示单元相连,所述BCG信号预处理单元包括陷波器、带通滤波器、LPF滤波器,所述陷波器设为最先接受信号,所述带通滤波器设为第二接受信号,所述LPF滤波器设为第三接受信号,所述心率波形显示单元与处理器连接,所述心率波形显示单元包括独立显示模组、终端电子设备。
5.根据权利要求4所述的心冲击图实时采集系统,其特征在于,所述陷波器为50Hz陷波器,所述带通滤波器为0.01Hz至15Hz的带通滤波器,所述LPF滤波器为截止频率为20Hz的LPF滤波器。
6.根据权利要求4所述的心冲击图实时采集系统,其特征在于,所述心率波形显示单元可以通过数据传输线与处理器连接,也可以通过无线方式与处理器连接,所述心率波形显示单元的终端电子设备包括手机、PAD、个人电脑。
Priority Applications (1)
Application Number | Priority Date | Filing Date | Title |
---|---|---|---|
CN202110025943.XA CN112869733B (zh) | 2021-01-08 | 2021-01-08 | 一种心冲击图实时心搏间期测算方法 |
Applications Claiming Priority (1)
Application Number | Priority Date | Filing Date | Title |
---|---|---|---|
CN202110025943.XA CN112869733B (zh) | 2021-01-08 | 2021-01-08 | 一种心冲击图实时心搏间期测算方法 |
Publications (2)
Publication Number | Publication Date |
---|---|
CN112869733A CN112869733A (zh) | 2021-06-01 |
CN112869733B true CN112869733B (zh) | 2021-12-24 |
Family
ID=76047424
Family Applications (1)
Application Number | Title | Priority Date | Filing Date |
---|---|---|---|
CN202110025943.XA Active CN112869733B (zh) | 2021-01-08 | 2021-01-08 | 一种心冲击图实时心搏间期测算方法 |
Country Status (1)
Country | Link |
---|---|
CN (1) | CN112869733B (zh) |
Families Citing this family (4)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
CN113476039B (zh) * | 2021-07-30 | 2024-03-19 | 深圳数联天下智能科技有限公司 | 心冲击信号的采集方法、装置、存储介质及计算机设备 |
CN114010186B (zh) * | 2022-01-11 | 2022-03-18 | 华南师范大学 | 一种心冲击图信号定位方法和计算机设备 |
CN114767127B (zh) * | 2022-06-21 | 2022-09-30 | 广州中科新知科技有限公司 | 一种心冲击图信号的处理方法及系统 |
CN117752326B (zh) * | 2023-11-14 | 2024-07-19 | 武汉光谷航天三江激光产业技术研究院有限公司 | 一种基于bcg信号的递归图和波形分布的心脏状态判断方法 |
Citations (8)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
CN102274020A (zh) * | 2011-06-30 | 2011-12-14 | 东北大学 | 一种低功耗的动态心电监护仪及其控制方法 |
CN107569226A (zh) * | 2017-09-27 | 2018-01-12 | 广州中科新知科技有限公司 | 基于压电传感获取hrv的方法及应用 |
ES2656765A1 (es) * | 2016-07-27 | 2018-02-28 | Universitat Politécnica de Catalunya | Método y aparato para detectar eventos sistólicos mecánicos a partir del balistocardiograma |
CN107951472A (zh) * | 2017-12-11 | 2018-04-24 | 青岛海信医疗设备股份有限公司 | 心脏冲击信号测量装置、方法及坐具 |
CN108294745A (zh) * | 2018-03-07 | 2018-07-20 | 武汉大学 | 多导联心电图信号中p波、t波起止点检测方法及系统 |
CN108922626A (zh) * | 2018-08-22 | 2018-11-30 | 华南师范大学 | 一种体征参数模型构建方法和体征参数评价方法 |
CN110916636A (zh) * | 2019-11-22 | 2020-03-27 | 新绎健康科技有限公司 | 一种基于动态二阶差分阈值的bcg信号心率计算方法及系统 |
CN111091116A (zh) * | 2019-12-31 | 2020-05-01 | 华南师范大学 | 一种用于判断心律失常的信号处理方法及系统 |
-
2021
- 2021-01-08 CN CN202110025943.XA patent/CN112869733B/zh active Active
Patent Citations (8)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
CN102274020A (zh) * | 2011-06-30 | 2011-12-14 | 东北大学 | 一种低功耗的动态心电监护仪及其控制方法 |
ES2656765A1 (es) * | 2016-07-27 | 2018-02-28 | Universitat Politécnica de Catalunya | Método y aparato para detectar eventos sistólicos mecánicos a partir del balistocardiograma |
CN107569226A (zh) * | 2017-09-27 | 2018-01-12 | 广州中科新知科技有限公司 | 基于压电传感获取hrv的方法及应用 |
CN107951472A (zh) * | 2017-12-11 | 2018-04-24 | 青岛海信医疗设备股份有限公司 | 心脏冲击信号测量装置、方法及坐具 |
CN108294745A (zh) * | 2018-03-07 | 2018-07-20 | 武汉大学 | 多导联心电图信号中p波、t波起止点检测方法及系统 |
CN108922626A (zh) * | 2018-08-22 | 2018-11-30 | 华南师范大学 | 一种体征参数模型构建方法和体征参数评价方法 |
CN110916636A (zh) * | 2019-11-22 | 2020-03-27 | 新绎健康科技有限公司 | 一种基于动态二阶差分阈值的bcg信号心率计算方法及系统 |
CN111091116A (zh) * | 2019-12-31 | 2020-05-01 | 华南师范大学 | 一种用于判断心律失常的信号处理方法及系统 |
Non-Patent Citations (2)
Title |
---|
Adaptive Beat-to-Beat Heart Rate Estimation;Christoph Bruser;《IEEE TRANSACTIONS ON INFORMATION TECHNOLOGY IN BIOMEDICINE》;20110930;全文 * |
Ballistocardiogram Signal Processing: A Literature;Ibrahim Sadek;《Review Paper》;20180703;全文 * |
Also Published As
Publication number | Publication date |
---|---|
CN112869733A (zh) | 2021-06-01 |
Similar Documents
Publication | Publication Date | Title |
---|---|---|
CN112869733B (zh) | 一种心冲击图实时心搏间期测算方法 | |
CN107714023B (zh) | 基于人工智能自学习的静态心电图分析方法和装置 | |
CN103860162B (zh) | 一种起搏信号检测方法、系统和心电检测设备 | |
CN107837082A (zh) | 基于人工智能自学习的心电图自动分析方法和装置 | |
CN108601543B (zh) | 一种ecg信号处理方法及装置 | |
CN101828918B (zh) | 基于波形特征匹配的心电信号r波峰检测方法 | |
US20090143693A1 (en) | Method and apparatus for generating determination indexes for identifying ecg interfering signals | |
JP2008520384A (ja) | 限られたプロセッサ資源で呼吸数を実時間判定するための方法及びシステム | |
KR101910982B1 (ko) | 개인화된 생체 신호 패턴을 이용한 생체 신호의 동잡음 제거 방법 및 장치 | |
CN103961089B (zh) | 基于分段直线拟合的窦性心率震荡趋势检测方法 | |
EP0458896A1 (en) | Wrist worn heart rate monitor | |
CN108903929B (zh) | 心率检测修正的方法、装置、存储介质和系统 | |
CN105286853A (zh) | 基于可穿戴设备的疾病检测方法及装置、可穿戴设备 | |
CN104408718B (zh) | 一种基于双目视觉测量的步态数据处理方法 | |
CN109700450A (zh) | 一种心率检测方法及电子设备 | |
US20180014756A1 (en) | Method for Determining Types of Human Motor Activity and Device for Implementing Thereof | |
CN106815570A (zh) | 一种基于动态模式识别的心电信号st‑t段识别方法 | |
EP3003140A1 (en) | Electrocardiogram analysis | |
CN106361325B (zh) | 一种便携式心电仪所测单导联心电图的筛查识别系统 | |
CN109009086A (zh) | 一种心电信号r波检测方法及系统 | |
CN105444763A (zh) | 一种imu室内定位方法 | |
KR20210001575A (ko) | 잡음 환경 내 실시간 심장 박동수 탐지 장치 및 그 방법 | |
CN111345814B (zh) | 一种心电信号中心拍的分析方法、装置、设备和存储介质 | |
CN108992063A (zh) | 一种心电信号的质量判断方法及系统 | |
US20170311899A1 (en) | Apparatus and method for identifying movement in a patient |
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 |