CN106691425A - 一种运动手环的腕部心率监测方法 - Google Patents

Abstract

本发明公开了一种运动手环的腕部心率监测方法，首先对原始光学脉搏波信号进行阈值处理，消除与心率不相关信号。然后通过奇异频谱分析的方法，对同步加速度信号进行处理，得到运动伪影参考信号。将运动伪影信号与预处理光学信号进行自适应滤波处理，将光学信号中的运动伪影成分消除，得到纯净脉搏波信号。接着使用自适应阈值迭代的方法，对纯净脉搏波信号进行处理，提高其频谱精度。最后针对在脉搏波频谱的峰值选择中可能出现的情况进行充分分析，以此确定脉搏波峰值。通过准确判定峰值与其经历的时间，进而可以计算当前时间下的心率值。本方法能有效的减少运动伪影对光学脉搏波信号产生的干扰，也进一步提高了心率估算的精度。

Description

一种运动手环的腕部心率监测方法

技术领域

本发明涉及可穿戴设备技术领域，尤其涉及一种腕部心率监测设备的运动伪影消除与心率值计算方法。

背景技术

随着生活水平的提高和移动互联网的发展，人们对健康越发重视，全民运动成为了一种潮流。心率作为一种常见的人体体征，在临床医学以及运动科学中均有重要的研究价值。与此同时，可穿戴设备的发展使人们在运动中实时测量人体心率成为可能。目前，越来越多的可穿戴设备拥有心率测定功能。然而，多数产品在运动中测算的心率数据与实际值之间较大误差，如何在剧烈运动中测得高精度心率值也成为目前热门研究领域。

现有的心率监测设备分为三种：胸带式、指夹式和腕式。出于运动的考虑，前两者不利于人们运动时的便携性，人们更多的选择腕式测量方法。腕部的测量方式又分为压电式、光电透射式、光电反射式。出于设备的简化，以光学体积脉搏波描记法技术(PPG)为代表的光电反射原理，越来越成为主流技术。由于人体血液与人体组织对绿光吸收引起光强度衰减变化，与心房搏动变化相一致，PPG技术在越来越多的心率监测设备中得到运用实践。出于步数及与智能手机相配合的多种功能的需要，加速度计芯片也越来越多的出现在各类运动手环、腕表中。

目前，多数腕部监测设备的心率估算误差较大，其原因主要在于运动伪影的引入与心率峰值判断策略的不完善。实现运动手环的高精度心率监测的主要任务在于，运动伪影的消除与心率峰值的判断选择。多项研究表明，三轴加速度信号中包含运动伪影信号成分。利用同步加速度信号去除运动伪影也成为主流研究方向。

发明内容

本发明要解决的技术问题是针对人体运动中，腕部心率监测设备采集信号中对运动伪影的消除，同时全面分析多种峰值判断场景，有效提高当前心率估算精度。

本发明的基本原理为：与脉搏波监测同步记录的加速度信号中，可以含有类似于运动伪影的波形成分。对加速度信号进行奇异值分解并重建后，得到运动伪影参考信号；通过对两路PPG信号与运动伪影参考信号进行自适应运动伪影消除处理，得到不受运动伪影信号影响的纯净PPG信号；多次自适应阈值迭代可以提高纯净PPG信号的精度，为峰值选择提供更为精确的数据；在心率峰值选择策略中，本方法对脉搏波功率谱可能出现的多种情况进行分析，从而进一步提高心率值估算精度。

本发明提供针对在运动情况下的腕部PPG信号进行运动伪影消除和心率峰值选择的方法，具体实现步骤如下：

步骤(1)，人体脉搏正常工作范围是40～210BPM，对应频率为0.6～3.5Hz，而运动伪影及其它噪声相应频率远超范围。首先需要对采集的双通道原始PPG信号进行阈值处理，以此去除部分与脉搏波不相关的噪声和部分运动伪影成分。进行阈值处理后的两个通道PPG信号分别记作S_PPG1和S_PPG2；

步骤(2)，三轴加速度数据同时包含运动伪影成分，对三轴加速度数据进行处理，产生运动伪影参考信号。对三轴加速度数据的具体操作包括嵌入、奇异值分解、分组、归一化处理、心率相关频谱移除五个步骤，流程如图2所示，具体如下：

步骤(2.1)，将每个时间窗口下加速度时间序列，映射为二维矩阵中。

步骤(2.2)，对步骤(2.1)中的二维矩阵进行奇异值分解，使之成为若干线性独立的分量。

步骤(2.3)，对所有线性独立的分量进行分组，使各分组满足同频或谐频相关特性。

步骤(2.4)，对分组进行对角线平均化处理，重建长度为M的时间序列。加速度信号能再次表达为g个时间序列的总和。

步骤(2.5)，使用周期图方法，估算加速度数据的频谱，用F_acc表示。在该频谱中同时包括与心率相关的频率分量，如果该部分被当作运动伪影从PPG信号中移除，将造成心率信号的失真。将前一时间窗中估算的心率基频和谐频的位置索引从F_acc中移除，记作F_accq，该信号即运动伪影参考信号。

步骤(3)，将F_accq分别与S_PPG1、S_PPG2进行多阶自适应滤波处理，将后者中与F_accq相同的频率成分去除，得到净化过的PPG信号。多阶自适应滤波处理原理图如图3。

步骤(4)，对纯净PPG信号进行自适应阈值迭代处理，使其恢复原始特性成为高分辨率去噪频谱，为后续峰值选择决策机制提供高精度信号。具体原理如图4。

步骤(5)，在非运动状态下，预先测量一个不受运动伪影干扰的心率值，以此作为心率初值。

步骤(6)，在峰值选择策略中，对如下三种情况分别判断：频谱存在明显峰值峰值选择；频谱中未出现无明显峰值，但存在谐波相关的频率选择；其他情况。

步骤(7)，判断当前时间窗峰值点个数与所经历的时间长度，以此计算当前时间窗的心率值。

有益效果

本发明采用的针对人体腕部剧烈运动时的PPG信号中运动伪影消除以及心率值估算的方法，通过结合同步加速度计信号运算处理，能明显减少运动伪影对PPG信号造成的干扰；同时，对多种峰值选择场景进行分析，有效的提高对心率估算结果的准确性。

附图说明

图1为本发明的方法原理图。

图2为奇异频谱分析流程图。

图3为基于多阶自适应滤波器的运动伪影消除方法。

图4为自适应阈值迭代方法的原理图。

具体实施方式

步骤(1):首先将采集的PPG原始信号输入带通滤波器(BPF)，进行阈值处理。本方法采用0.4～5Hz的带通滤波器，一方面尽可能的去除非心率频带信号，另一方面为了考虑人体心率异常情况的发生。该操作可以在一定程度上去除运动伪影和其他不相关噪声。得到的双路PPG信号分别记为S_PPG1、S_PPG2。

步骤(2)，通过奇异频谱分析来分解三轴加速度的方式，生成参考MA信号。具体包括嵌入、奇异值分解、分组、归一化处理、心率相关频谱移除五个步骤，具体如下：

步骤(2.1)，对于每个维度的加速度数据，均可看成与时间相关的序列其中各分量y_i,i∈(1,M)为连续时间上的加速度信号在该维度上的分量，M为时间序列的长度，本方法中M取值为8，T在此处代表矩阵的转置；将时间序列y映射为一个L*K矩阵中，此处将该矩阵称为L轨迹矩阵，其中L、K分别为所映射的矩阵的行数和列数，L、K与M的关系满足K＝M+1-L,L＜M/，2本方法中M取值为1000，L取值为100，K取值为901。该L轨迹矩阵用符号Y表示，如下式所示：

步骤(2.2)，在该步骤中，将上述步骤生成的三组加速度分量的L轨迹矩阵分别进行奇异值分解，具体原理如下公式，

其中，σ_i是第i个奇异值，u_i是相应的左奇异向量，v_i是相应的右奇异向量，各Y_i之间均线性独立。本方法中d的取值为100，因此i的取值范围为1到100。

步骤(2.3)，将个d一维行矩阵Y_i重新分配为g组，其中g≤d，如下式所示。本方法中g的取值为100，p的取值范围为1到100。每组包含的矩阵均满足同频或者谐频的特性。

步骤(2.4)，在该步骤中，每组通过对角线平均化来重建一个时间序列y_p，它的长度为M。原始信号y能再次表达为g个时间序列的总和。

步骤(2.5)，使用周期图方法，从确定的主频中估算加速度数据的每个信道的频谱。其中，主频率是在给定的频谱中其幅值高于最大幅值一半的谱峰。在频谱中所选的主频位置索引的集合，用F_acc表示。将先前时间窗中估算的心率基频和谐频点从F_acc中移除，记作F_accq，该信号当作运动伪影参考信号。此举避免将与心跳相关频率分量被当作运动伪影成分，在后续步骤中从PPG原始信号中被移除，造成PPG信号失真。

步骤(3)，将通过步骤(1)的带通滤波器的PPG信道1、PPG信道2分别与一阶自适应滤波器连接，各个阶的自适应滤波器连接顺次串联，各个参数运动伪信号作为待滤除分量分别与各自适应滤波器连接，最终由第N个自适应滤波器输出净化过的双路PPG信号。

步骤(4)，对纯净PPG信号进行自适应阈值迭代处理。具体原理如图4。

步骤(4.1)，初次阈值处理时先进行初始化操作，输入全零块。

步骤(4.2)，对输入进行离散转化，使其从稀疏域转化为信息域。

步骤(4.3)，将噪声信号的信息域样本替换为精确地信息域样本。

步骤(4.4)，对替换后的信号进行离散反变化，使其从信息域变换为稀疏域。

步骤(4.5)，对稀疏域信号进行阈值处理，具体如公式

其中，X_in表示离散反变化产生的稀疏域信号，阈值系数α＝0.1，阈值系数β＝2，i＝5表示迭代次数。

步骤(4.6)，将X_out作为步骤(4.2)的输入，进行下一轮的自适应阈值迭代。

步骤(5)，初始化系统，保持腕部静止，记录前8秒时间窗的频谱的最高值作为心率初值。

步骤(6)，在某一时刻时间窗口中，用N_cur表示当前时间窗内估算的心率频率点的个数，用N_prev来表示先前时间窗估算的心率频率位置。在当前时间窗N_cur中设置三个频率变化范围，找出三个频率范围即第一频率范围R₁、第二频率范围R₂、第三频率范围R₃，如下式所示：

R₁＝[N_prev-ε₁,N_prev+ε₁]

R₂＝[2×N_prev-ε₂,2×N_prev+ε₂]

R₃＝[3×N_prev-ε₃,3×N_prev+ε₃]

其中N_prev为500，ε₁，ε₂，ε₃为各频率范围的领域半径，本方法中分别为60，60，100。

用P₁、P₂、P₃表示R₁、R₂、R₃的最高峰值频率位置；将S₁、S₂、S₃记为P₁、P₂、P₃点处相应的信号频谱值。因此，信道一的频率位置值表示为P₁₁、P₂₁、P₃₁，S₁₁、S₂₁、S₃₁则分别表示为频率位置值对应的信号频谱值；信道二的频率峰值位置值表示为P₁₂、P₂₂、P₃₂，S₁₂、S₂₂、S₃₂则分别表示为频率位置值对应的信号频谱值。综合考虑三种可能出现的情况，需完成以下步骤：

步骤(6.1)，若S₁₁、S₂₁、S₃₁、S₁₂、S₂₂、S₃₂中某一值明显大于其他显著峰值点(定义如公式)，将当前时间窗的心率频率值当做相应的基频。N_cur定义如以下公式所示，

其中T为算法最优而预定义的阈值，本方法中T＝0.6。

步骤(6.2)，若两个PPG信道均没有明显的波峰，在P₁₁、P₂₁、P₃₁、P₁₂、P₂₂、P₃₂中找出一个谐波相关的峰值对，并将它们的基频作为当前时间窗心率频率指数，具体如下公式所示，

其中，δ为一任意小整数，本方法中δ＝9。

步骤(6.3)，若以上两种情况均未出现，则建立一个10s的时间窗，用以连接当前和先前窗口，同样在该区域定义上述三个谐波变化范围。用Q₁₁、Q₂₁、Q₃₁、Q₁₂、Q₂₂、Q₃₂表示这些频率变换范围的最高峰值频率位置指数。将所有可利用的基础频率的平均值作为N_cur。如以下公式所示。

步骤(7)，对当前时间窗内的峰值计数，并以B表示；计算峰值间使用时间，以时间t来表示。因此，当前时刻的估算心率值H，如公式表示，其单位为BPM。

H＝60B/t

结束。

Claims (2)

1.一种运动手环的腕部心率监测方法，与脉搏波监测同步记录的加速度信号中，含有类似于运动伪影的波形成分；对加速度信号进行奇异值分解并重建后，得到运动伪影参考信号；通过对两路PPG信号与运动伪影参考信号进行自适应运动伪影消除处理，得到不受运动伪影信号影响的纯净PPG信号；多次自适应阈值迭代可以提高纯净PPG信号的精度，为峰值选择提供更为精确的数据；在心率峰值选择策略中，本方法对脉搏波功率谱可能出现的多种情况进行分析，从而进一步提高心率值估算精度；

其特征在于：具体实现步骤如下，

步骤(1)，人体脉搏正常工作范围是40～210BPM，对应频率为0.6～3.5Hz，而运动伪影及其它噪声相应频率远超范围；首先需要对采集的双通道原始PPG信号进行阈值处理，以此去除部分与脉搏波不相关的噪声和部分运动伪影成分；进行阈值处理后的两个通道PPG信号分别记作S_PPG1和S_PPG2；

步骤(2)，三轴加速度数据同时包含运动伪影成分，对三轴加速度数据进行处理，产生运动伪影参考信号；对三轴加速度数据的具体操作包括嵌入、奇异值分解、分组、归一化处理、心率相关频谱移除五个步骤，具体如下：

步骤(2.1)，将每个时间窗口下加速度时间序列，映射为二维矩阵中；

步骤(2.2)，对步骤(2.1)中的二维矩阵进行奇异值分解，使之成为若干线性独立的分量；

步骤(2.3)，对所有线性独立的分量进行分组，使各分组满足同频或谐频相关特性；

步骤(2.4)，对分组进行对角线平均化处理，重建长度为M的时间序列；加速度信号能再次表达为g个时间序列的总和；

步骤(2.5)，使用周期图方法，估算加速度数据的频谱，用F_acc表示；在该频谱中同时包括与心率相关的频率分量，如果该部分被当作运动伪影从PPG信号中移除，将造成心率信号的失真；将前一时间窗中估算的心率基频和谐频的位置索引从F_acc中移除，记作F_accq，该信号即运动伪影参考信号；

步骤(3)，将F_accq分别与S_PPG1、S_PPG2进行多阶自适应滤波处理，将后者中与F_accq相同的频率成分去除，得到净化过的PPG信号；

步骤(4)，对纯净PPG信号进行自适应阈值迭代处理，使其恢复原始特性成为高分辨率去噪频谱，为后续峰值选择决策机制提供高精度信号；；

步骤(5)，在非运动状态下，预先测量一个不受运动伪影干扰的心率值，以此作为心率初值；

步骤(6)，在峰值选择策略中，对如下三种情况分别判断：频谱存在明显峰值峰值选择；频谱中未出现无明显峰值，但存在谐波相关的频率选择；

步骤(7)，判断当前时间窗峰值点个数与所经历的时间长度，以此计算当前时间窗的心率值。

2.根据权利要求1所述的一种运动手环的腕部心率监测方法，其特征在于：

步骤(1):首先将采集的PPG原始信号输入带通滤波器，进行阈值处理；本方法采用0.4～5Hz的带通滤波器，一方面尽可能的去除非心率频带信号，另一方面为了考虑人体心率异常情况的发生；该操作可以在一定程度上去除运动伪影和其他不相关噪声；得到的双路PPG信号分别记为S_PPG1、S_PPG2；

步骤(2)，通过奇异频谱分析来分解三轴加速度的方式，生成参考MA信号；具体包括嵌入、奇异值分解、分组、归一化处理、心率相关频谱移除五个步骤，具体如下：

步骤(2.1)，对于每个维度的加速度数据，均可看成与时间相关的序列其中各分量y_i,i∈(1,M)为连续时间上的加速度信号在该维度上的分量，M为时间序列的长度，本方法中M取值为8，T在此处代表矩阵的转置；将时间序列y映射为一个L*K矩阵中，此处将该矩阵称为L轨迹矩阵，其中L、K分别为所映射的矩阵的行数和列数，L、K与M的关系满足K＝M+1-L,L＜M/2，本方法中M取值为1000，L取值为100，K取值为901；该L轨迹矩阵用符号Y表示，如下式所示：

步骤(2.2)，在该步骤中，将上述步骤生成的三组加速度分量的L轨迹矩阵分别进行奇异值分解，具体原理如下公式，

$Y=\underset{i=1}{\overset{d}{\mathrm{\Σ}}}{Y}_i,d\stackrel{\mathrm{\Δ}}{=}\min \{L,K\}$

其中，σ_i是第i个奇异值，u_i是相应的左奇异向量，v_i是相应的右奇异向量，各Y_i之间均线性独立；本方法中d的取值为100，因此i的取值范围为1到100；

步骤(2.3)，将个d一维行矩阵Y_i重新分配为g组，其中g≤d，如下式所示；本方法中g的取值为100，p的取值范围为1到100；每组包含的矩阵均满足同频或者谐频的特性；

$Y=\underset{p=1}{\overset{g}{\mathrm{\Σ}}}{Y}_{I_p}$

步骤(2.4)，在该步骤中，每组通过对角线平均化来重建一个时间序列y_p，它的长度为M；原始信号y能再次表达为g个时间序列的总和；

步骤(2.5)，使用周期图方法，从确定的主频中估算加速度数据的每个信道的频谱；其中，主频率是在给定的频谱中其幅值高于最大幅值一半的谱峰；在频谱中所选的主频位置索引的集合，用F_acc表示；将先前时间窗中估算的心率基频和谐频点从F_acc中移除，记作F_accq，该信号当作运动伪影参考信号；此举避免将与心跳相关频率分量被当作运动伪影成分，在后续步骤中从PPG原始信号中被移除，造成PPG信号失真；

步骤(3)，将通过步骤(1)的带通滤波器的PPG信道1、PPG信道2分别与一阶自适应滤波器连接，各个阶的自适应滤波器连接顺次串联，各个参数运动伪信号作为待滤除分量分别与各自适应滤波器连接，最终由第N个自适应滤波器输出净化过的双路PPG信号；

步骤(4)，对纯净PPG信号进行自适应阈值迭代处理；

步骤(4.1)，初次阈值处理时先进行初始化操作，输入全零块；

步骤(4.2)，对输入进行离散转化，使其从稀疏域转化为信息域；

步骤(4.3)，将噪声信号的信息域样本替换为精确地信息域样本；

步骤(4.4)，对替换后的信号进行离散反变化，使其从信息域变换为稀疏域；

步骤(4.5)，对稀疏域信号进行阈值处理，具体如公式

$X_{out}=\left\{\begin{array}{c}{X}_{in},\left|X_{in}\right|\≥{\mathrm{\βe}}^{-\mathrm{\α}i}\\ 0,\left|X_{in}\right|\≤{\mathrm{\βe}}^{-\mathrm{\α}i}\end{array}\right.$

其中，X_in表示离散反变化产生的稀疏域信号，阈值系数α＝0.1，阈值系数β＝2，i＝5表示迭代次数；

步骤(4.6)，将X_out作为步骤(4.2)的输入，进行下一轮的自适应阈值迭代；

步骤(5)，初始化系统，保持腕部静止，记录前8秒时间窗的频谱的最高值作为心率初值；

步骤(6)，在某一时刻时间窗口中，用N_cur表示当前时间窗内估算的心率频率点的个数，用N_prev来表示先前时间窗估算的心率频率位置；在当前时间窗N_cur中设置三个频率变化范围，找出三个频率范围即第一频率范围R₁、第二频率范围R₂、第三频率范围R₃，如下式所示：

R₁＝[N_prev-ε₁,N_prev+ε₁]

R₂＝[2×N_prev-ε₂,2×N_prev+ε₂]

R₃＝[3×N_prev-ε₃,3×N_prev+ε₃]

其中N_prev为500，ε₁，ε₂，ε₃为各频率范围的领域半径，本方法中分别为60，60，100；

用P₁、P₂、P₃表示R₁、R₂、R₃的最高峰值频率位置；将S₁、S₂、S₃记为P₁、P₂、P₃点处相应的信号频谱值；因此，信道一的频率位置值表示为P₁₁、P₂₁、P₃₁，S₁₁、S₂₁、S₃₁则分别表示为频率位置值对应的信号频谱值；信道二的频率峰值位置值表示为P₁₂、P₂₂、P₃₂，S₁₂、S₂₂、S₃₂则分别表示为频率位置值对应的信号频谱值；综合考虑三种可能出现的情况，需完成以下步骤：

步骤(6.1)，若S₁₁、S₂₁、S₃₁、S₁₂、S₂₂、S₃₂中某一值明显大于其他显著峰值点，将当前时间窗的心率频率值当做相应的基频；N_cur定义如以下公式所示，

$N_{cur}=\frac{p_{ij}}{i}\⇔\∀(k,l)\≠(i,j):S_{ij}\×T>S_{kl}$

其中T为算法最优而预定义的阈值，本方法中T＝0.6；

步骤(6.2)，若两个PPG信道均没有明显的波峰，在P₁₁、P₂₁、P₃₁、P₁₂、P₂₂、P₃₂中找出一个谐波相关的峰值对，并将它们的基频作为当前时间窗心率频率指数，具体如下公式所示，

$N_{cur}=\frac{\frac{p_{ij}}{i}+\frac{p_{kt}}{k}}{2}\&DoubleLeftRightArrow;\frac{p_{ij}}{i}-\frac{p_{kt}}{k}<\mathrm{\&delta;}$

其中，δ为一任意小整数，本方法中δ＝9；

步骤(6.3)，若以上两种情况均未出现，则建立一个10s的时间窗，用以连接当前和先前窗口，同样在该区域定义上述三个谐波变化范围；用Q₁₁、Q₂₁、Q₃₁、Q₁₂、Q₂₂、Q₃₂表示这些频率变换范围的最高峰值频率位置指数；将所有可利用的基础频率的平均值作为N_cur；如以下公式所示；

$N_{cur}=\frac{1}{12}(P_{11},\frac{P_{21}}{2},\frac{P_{31}}{3},P_{12},\frac{P_{22}}{2},\frac{P_{32}}{3},Q_{11},\frac{Q_{21}}{2},\frac{Q_{31}}{3},Q_{12},\frac{Q_{22}}{2},\frac{Q_{32}}{3})$

步骤(7)，对当前时间窗内的峰值计数，并以B表示；计算峰值间使用时间，以时间t来表示；因此，当前时刻的估算心率值H，如公式表示，其单位为BPM；

H＝60B/t

结束。