一种运动状态下心率计算方法、系统及计算机存储介质

专利名称：一种运动状态下心率计算方法、系统及计算机存储介质

专利类型：发明专利

专利申请号：CN202210069711.9

专利申请（专利权）人：中国人民解放军陆军军医大学
权利人地址：重庆市沙坪坝区高滩岩正街30号

专利发明（设计）人：龚渝顺, 魏良, 李永勤, 谢佳玲, 何密, 陈碧华, 王建杰

专利摘要：本发明提出了一种运动状态下心率计算方法、系统及计算机存储介质。该方法为：采集心电信号，实时同步截取时长为t的数据，对其进行预处理后得到心电信号数据；对预处理后的心电信号数据进行膨胀运算、粗粒化，得到类正弦心率信号；对类正弦心率信号进行傅里叶变换，得到频谱；于频谱中寻找到峰值频率，根据峰值频率计算心率和峰值频率能量占比；若该占比大于设定值，则输出心率；若该占比小于设定值，计算当前心率与上一次心率计算结果之间的差值，若差值小于α，则输出心率；若差值大于α，则对频谱以峰值频率为中心进行抑制，将抑制后的频谱信号作为新的频谱重新执行该步骤，直至达到心率输出标准或迭代次数达到预设次数后，输出当前心率。本方法提高了心率检测的准确性。

主权利要求：一种运动状态下心率计算方法、系统及计算机存储介质
技术领域
[0001] 本发明涉及心率检测领域，具体涉及一种运动状态下心率计算方法、系统及计算机存储介质。
背景技术
[0002] 目前穿戴式心率检测已经逐步融入了人们的日常生活，基于光电溶剂脉搏波(PPG)的心率检测方法，在智能手表、手腕中得到了广泛的应用，穿戴式心率检测技术在睡眠监护、日常活动监测以及运动评估等领域发挥着越来越重要的作用。
[0003] 但由于受到PPG自身信号特点的制约，基于PPG的心率检测极易受到运动干扰的影响而失准，尤其在身体动作幅度较大、运动较为剧烈的情况下，欲从PPG信号中准确分析出心率信息是一项非常困难的任务。
[0004] 与PPG信号不同，心电信号中R波的频谱主要分布在较高的频率范围上，与日常活动引起的干扰频带之间的重叠较少，因此在受到较强运动干扰时，也能够保留足够的有效信息。随着基于导电织物的穿戴式心电监护技术的崛起，使基于心电信号的心率检测成为可能。虽然从心电信号在抗运动干扰方面与PPG相比具有天然的优势，但是人体剧烈的活动，仍然会引入大量的运动干扰，尤其对于非粘贴类的导电织物型心电电极就更加明显，从而影响心率检测的准确性和可靠性。
发明内容
[0005] 为了克服上述现有技术中存在的缺陷，本发明的目的是提供一种运动状态下心率计算方法、系统及计算机存储介质。
[0006] 为了实现本发明的上述目的，本发明提供了一种运动状态下心率计算方法，包括以下步骤：
[0007] 采集心电信号，并实时同步截取时长为t的数据，对其进行预处理后得到心电信号数据x(n)；
[0008] 对心电信号数据x(n)进行膨胀运算并粗粒化提取，得到类正弦心率信号z(n)；
[0009] 对类正弦心率信号z(n)进行傅里叶变换，得到频谱s(f)；
[0010] 于频谱s(f)中寻找到峰值频率，根据峰值频率计算心率和峰值频率能量占比；若峰值频率能量占比大于设定值，则输出心率；若峰值频率能量占比小于设定值，则将当前心率与上一次心率计算结果进行比较，若两者的差值的绝对值小于α，α为正实数，则输出心率；若两者的差值的绝对值大于α，则对频谱s(f)以峰值频率为中心进行抑制得到抑制后的频谱信号s1(f)，将抑制后的频谱信号s1(f)作为新的频谱s(f)重新执行该步骤，直至达到心率输出标准或迭代次数达到预设次数后，输出当前计算得到的心率。
[0011] 本方法不去定位R波位置，而是通过对整段信号的运算和变换，将R波出现的心率信息变换为拟正弦信号，通过信号的频谱分布中的峰值频率来反推心率值，并在严重干扰的情况下，利用峰值频率能量占比以及连续时间的心率比对，判断心率计算的正确性，通过干扰频谱能量的抑制和心率频率的重定位，对计算有误的心率进行了修正，提高了心率检测的准确性。
[0012] 该运动状态下心率计算方法的优选方案：获取心电信号数据x(n)后，识别心电信号数据x(n)中是否存在高幅度尖刺干扰，如存在高幅度尖刺干扰，则对高幅度尖刺干扰进行抑制，然后对抑制后的心电信号数据y(n)进行膨胀运算并粗粒化提取，得到类正弦心率信号z(n)；如不存在高幅度尖刺干扰，则对心电信号数据x(n)进行膨胀运算并粗粒化提取，得到类正弦心率信号z(n)。
[0013] 该运动状态下心率计算方法的优选方案：所述预处理包括对截取的时长为t的数据进行带通滤波，得到心电信号数据x(n)。经过滤波后的信号突出了R波特征，抑制了其他心电成分以及运动干扰成分。
[0014] 该运动状态下心率计算方法的优选方案：对高幅度尖刺干扰进行抑制的步骤为：
[0015] 构建加权函数w(n)对心电信号数据x(n)中的高幅度尖刺干扰进行抑制，得到信号y(n)：y(n)＝w(n)·x(n)；
[0016] 加权函数w(n)的计算步骤为：
[0017] 首先定位心电信号数据x(n)中所有大于幅度阈值Th的数据样本位置posi，并计算其宽度wide，以长度L为窗宽，1个采样点为移动步长对心电信号数据x(n)进行峰峰值ppv的计算，得到连续的峰峰值信号xp(n)，并使xp(n)数据长度与心电信号数据x(n)保持一致，其中L＝δ×fs，fs为采样率，δ为系数；
[0018] 加权函数 其中，hanning(·)为窗函数，Th为识别高幅度尖刺干扰的幅度阈值。

[0019] 该运动状态下心率计算方法的优选方案：膨胀运算的具体步骤为：
[0020] 以L为长度的移动窗，1个采样点为移动步长，对x(n)或y(n)进行峰峰值ppv的计算，得到连续的峰峰值信号xp1(n)，L＝δ×fs，fs为采样率，δ为系数；以均值 为阈值计算信号xp1(n)的占空比 其中 为满足条件 的数据样点数；

[0021] 修正移动窗口长度 η为系数；
[0022] 以Lnew为长度的移动窗，1个采样点为移动步长，对x(n)或y(n)进行峰峰值ppv的计算，得到连续的峰峰值信号xp2(n)。
[0023] 粗粒化提取的具体步骤为：对xp2(n)进行长度为Lnew移动平均滤波，得到类正弦心率信号z(n)。
[0024] 这里进行两次移动窗操作使得信号的占空比接近于0.5，这样可以使后面计算傅里叶变换时心率对应的峰值频率单频成分更加突出。
[0025] 该运动状态下心率计算方法的优选方案：根据上一次心率计算结果对α进行分段设置，上一次心率计算结果越高，α越小。这提高了心率修正的准确性。
[0026] 本发明提出了一种运动状态下心率计算系统，包括采集心电信号的心电监护设备、存储器和控制器，所述心电监护设备与控制器通信连接，向控制器发送实时采集的心电信号，所述存储器与控制器通信连接，所述存储器用于存放至少一可执行指令，所述可执行指令使所述控制器对所述心电信号执行如上述的运动状态下心率计算方法对应的操作。
[0027] 本发明提出了一种计算机存储介质，所述存储介质中存储有至少一可执行指令，所述可执行指令使处理器执行如上述的运动状态下心率计算方法对应的操作。
[0028] 本发明的有益效果是：本发明在计算心率时，不去定位R波位置，而是通过对整段信号的运算和变换，将R波出现的心率信息变换为拟正弦信号，通过信号的频谱分布中的峰值频率来反推心率值，如果因干扰导致个别心跳位置出现偏差、遗漏，仅仅会对频率分布的旁瓣能量造成影响，使得峰值频率稍许降低，但是并不会改变主要的心率值。此外，本发明对高幅度尖刺干扰抑制，并通过连续心率比对，对心电频谱分布进行分析和修正，抑制强干扰引起的干扰频带能量，完成对心率的修正，进一步提高了强干扰状态下心率检测的准确性。
[0029] 本发明的附加方面和优点将在下面的描述中部分给出，部分将从下面的描述中变得明显，或通过本发明的实践了解到。
附图说明
[0030] 本发明的上述和/或附加的方面和优点从结合下面附图对实施例的描述中将变得明显和容易理解，其中：
[0031] 图1是该运动状态下心率计算方法的流程示意图；
[0032] 图2是轻微干扰心电信号心率计算波形图；
[0033] 图3是中度干扰心电信号心率计算波形图。
[0034] 图4是严重干扰心电信号心率计算波形图；
[0035] 图5是高幅度干扰心电信号心率计算波形图。
具体实施方式
[0036] 下面详细描述本发明的实施例，所述实施例的示例在附图中示出，其中自始至终相同或类似的标号表示相同或类似的元件或具有相同或类似功能的元件。下面通过参考附图描述的实施例是示例性的，仅用于解释本发明，而不能理解为对本发明的限制。
[0037] 如图1所示，本发明提供了一种运动状态下心率计算方法的实施例，该实施的具体步骤如下：
[0038] 首先，利用穿戴式心电监护设备采集心电信号，并实时同步截取时长为t的数据，对其进行预处理后得到心电信号数据x(n)。
[0039] 本实施例中，t为正数，通常取3‑10秒，同时为了突出R波特征，抑制其他心电成分以及运动干扰成分，预处理包括对同步截取时长为t的数据进行带通滤波，带通滤波的带宽优选但不限于为5.5～75Hz，带通滤波得到的数据即为心电信号数据x(n)。
[0040] 获取心电信号数据x(n)后，识别心电信号数据x(n)中是否存在高幅度尖刺干扰，如存在高幅度尖刺干扰，则对高幅度尖刺干扰进行抑制，然后对抑制后的心电信号数据y(n)进行膨胀运算并粗粒化提取，得到类正弦心率信号z(n)；如不存在高幅度尖刺干扰，则对心电信号数据x(n)进行膨胀运算并粗粒化提取，得到类正弦心率信号z(n)。
[0041] 本实施例中通过设置识别高幅度尖刺干扰的幅度阈值Th以判断心电信号数据x(n)中是否存在高幅度尖刺干扰，Th优选但不限于为4mv，若心电信号数据x(n)中存在大于Th的样点，则认为存在高幅度尖刺干扰，若心电信号数据x(n)中不存在大于Th的样点，则认为不存在高幅度尖刺干扰。
[0042] 这里对高幅度尖刺干扰进行抑制的步骤为：
[0043] 构建加权函数w(n)对心电信号数据x(n)中的高幅度尖刺干扰进行抑制，得到信号y(n)：y(n)＝w(n)x(n)。
[0044] 加权函数w(n)的计算步骤为：
[0045] 首先定位x(n)中所有大于幅度阈值Th的数据样本位置posi，并计算其宽度wide，以长度L为窗宽，1个采样点为移动步长对x(n)进行峰峰值ppv的计算，得到连续的峰峰值信号xp(n)，并优选但不限于采用padding的方法使xp(n)数据长度与x(n)保持一致，其中L＝δ×fs，fs为采样率，δ为系数，优选但不限于为0.08。padding的意思就是在信号边缘上进行补点操作，使滑窗运算后的输出点数与输入点数相同。
[0046] 加权函数 其中，hanning(·)为窗函数。

[0047] 膨胀运算的具体步骤为：
[0048] 以L为长度的移动窗，1个采样点为移动步长，对x(n)或y(n)进行峰峰值ppv的计算(如没有高幅度尖刺干扰，则对x(n)进行峰峰值ppv的计算，如有高幅度尖刺干扰，则对高幅度尖刺干扰抑制后的y(n)进行峰峰值ppv的计算)，得到连续的峰峰值信号xp1(n)，L＝δ×fs，fs为采样率，δ为系数；以均值 为阈值计算信号xp1(n)的占空比其中 为满足条件 的数据样点数。

[0049] 修正移动窗口长度 η为系数，优选但不限于为0.9。
[0050] 以Lnew为长度的移动窗，1个采样点为移动步长，对x(n)或y(n)进行峰峰值ppv的计算，得到连续的峰峰值信号xp2(n)，再对xp2(n)进行长度为Lnew移动平均滤波，得到类正弦心率信号z(n)。
[0051] 然后对类正弦心率信号z(n)进行傅里叶变换，得到频谱s(f)，为提高频谱信号的分辨率，本实施例优选但不限于进行1024*16点数的傅里叶变换。
[0052] 于频谱s(f)中寻找到峰值频率fpeak，根据峰值频率fpeak计算心率HR和峰值频率能量占比ratio，其中，心率HR＝60fpeak，峰值频率能量占比 fpeak为峰值频率，μ为常数，通常取值范围为(0.2，0.7)，本实施例中优选但不限于为0.3，s(fi)为频谱s(f)中对应频率fi时的频谱能量。若峰值频率能量占比ratio大于设定值，本实施例中设定值优选但不限于为0.23，则输出心率；若峰值频率能量占比ratio小于设定值，则将当前心率HR与上一次心率计算结果HRlast进行比较，若两者的差值的绝对值小于α，α为正实数，则输出心率；若两者的差值的绝对值大于α，则对频谱s(f)以峰值频率fpeak为中心进行抑制，本实施例中抑制方法优选但不限于采用s(fj)＝0,j∈(fpeak‑0.5,fpeak+0).5，得到抑制后的频谱信号s1(f)，将抑制后的频谱信号s1(f)作为新的频谱s(f)重新执行该步骤计算峰值频率fpeak，估算心率HR，直至达到心率输出标准或迭代次数达到预设次数后，输出当前计算得到的心率。本实施例中，迭代次数优选但不限于为3次。

[0053] 关于α的取值，可根据上一次心率计算结果对α进行分段设置，上一次心率计算结果越高，α越小。
[0054] 本实施例中，
[0055] 本实施例中的术语解释：抑制就是指减少其能量幅度，在数值上就是减小数据的大小，膨胀理解为将心电信号R波扩展变宽，粗粒化理解为以固定的尺度(采样点数)连续取平均值，以创建不同尺度下(或分辨率)的信号，峰峰值是信号在一定尺度上波峰与波谷的幅值差值。
[0056] 下面以具体实验为例进行介绍：
[0057] 以轻微干扰为例，如图2所示，图2中(a)为截取的心电信号(4s)，图2中(b)是经过滤波后的信号，突出R波特征，抑制其他心电成分以及运动干扰成分，图2中(c)对(b)做膨胀运算(利用固定窗口进行峰峰值移动滤波；在运算时需要循环做两次，第一次得到一个类似(c)的图，再调整新的窗口长度，重新对(b)做一次运算即可得到(c)，其目的是为了使得(c)中的占空比接近于0.5，这样可以使后面计算傅里叶变换时心率对应的峰值频率单频成分更加突出)；图2中(d)对(c)进行加了滑窗的移动平均滤波，即粗粒化提取；图2中(e)傅里叶变换，找到峰值频率，计算心率。
[0058] 图2中的ratio＝0.42，大于了0.23，所以认为干扰影响不大，可直接输出结果。
[0059] 下面对中度干扰、严重干扰以及高幅度干扰进行介绍，基本也是这个流程，但是会有一些细微的变化。
[0060] 中度干扰，如图3所示，由于受到了较大运动干扰影响(从原始波形可以看出)，在处理后的波形中带有部分干扰，频谱分布中，出现了第二条干扰条带，但是主要能量仍然集中在峰值频率(心率对应)处。计算的能量占比ratio小于了0.23，认为受到了干扰，因此与前一次心率进行比较，其差值为5，小于步骤6中的阈值，认为计算可信。
[0061] 严重干扰，如图4所示，严重干扰使得信号中的噪声很大，频谱中有大部分的干扰能量分布，出现的两个峰值中，最大的一个对应的是心率224，其能量占比0.11远小于0.23，因此比对前一次心率121，差值123，大于了步骤6中给出的阈值，因此对频谱分布进行峰值频率消除(虚线部分)；之后再次计算峰值频率对应的心率122，ratio为0.17，再次比对前次计算结果，其差值为1，在步骤6中的阈值范围内，因此认可结果122。
[0062] 高幅度干扰，如图5所示，信号中包含有一个短时的高幅度干扰，对其进行消除，得到(c)波形图；可以看到高幅度尖刺干扰得到有效抑制，虽然其左边的R波也被消减，但是并不影响整体的计算结果。
[0063] 本发明还提出了一种运动状态下心率计算系统的实施例，其包括采集心电信号的心电监护设备、存储器和控制器，所述心电监护设备与控制器通信连接，向控制器发送实时采集的心电信号，所述存储器与控制器通信连接，所述存储器用于存放至少一可执行指令，所述可执行指令使所述控制器对所述心电信号执行如上述的运动状态下心率计算方法对应的操作。
[0064] 本发明还提出了一种计算机存储介质的实施例，所述存储介质中存储有至少一可执行指令，所述可执行指令使处理器执行如上述的运动状态下心率计算方法对应的操作。
[0065] 在本说明书的描述中，参考术语“一个实施例”、“一些实施例”、“示例”、“具体示例”、或“一些示例”等的描述意指结合该实施例或示例描述的具体特征、结构、材料或者特点包含于本发明的至少一个实施例或示例中。在本说明书中，对上述术语的示意性表述不一定指的是相同的实施例或示例。而且，描述的具体特征、结构、材料或者特点可以在任何的一个或多个实施例或示例中以合适的方式结合。
[0066] 尽管已经示出和描述了本发明的实施例，本领域的普通技术人员可以理解：在不脱离本发明的原理和宗旨的情况下可以对这些实施例进行多种变化、修改、替换和变型，本发明的范围由权利要求及其等同物限定。

专利地区：重庆

专利申请日期：2022-01-21

专利公开日期：2023-09-01

专利公告号：CN114521880B