一种基于心跳二次谐波重构的心率估计方法与流程

1.本发明涉及一种心率估计方法，具体是一种基于心跳二次谐波重构的心率估计方法，属于智能健康监测技术领域。


背景技术：

2.心率监测医学监护领域具有重要意义，心率变化和心脏疾病密切相关，心率是最直接反映我们心脏健康的标志。据世界卫生组织统计，心血管病死亡率居高不下，并呈现年轻化趋势，所以，不管是对于上班族来说还是中老年人来说，实时监测心率都是非常必要的，通过连续准确的心率监测和分析，才能及时发现心脏活动的异常数据，为患者争取宝贵的治疗时间。
3.目前，心率监测主要依靠心电图仪器等专业监测仪器，需要通过传感器、电极等与人体接触才能监测出心跳频率，这种方法虽然测量精确度高，但是其需要专业医生的操作，且操作步骤繁琐，患者在测量过程中行动不便，不利于进行实时的心率监测。因此，非接触式体征监测仪器应运而生，非接触体征监测仪器具有监测方便、应用场景更加广泛等优点，特别是对婴幼儿、传染病患者、严重烧伤患者以及睡眠呼吸暂停患者等非接触式体征监测是不可或缺的。
4.近年来，基于毫米波雷达的非接触式生命体征监测越来越受到人们的关注，毫米波雷达测量心率具有精度高、抗干扰能力强、非接触式等显著优点。目前，常用的模态分解方法是心跳信号分离方法，但是经验模态分解在胸壁唯一信号分为为几个本征模态函数的过程中会产生模态混合问题，进而导致一个本征模态函数里同时包含呼吸谐波和心跳，影响心跳估计准确度。另外，还有人提出通过其他方法来削弱呼吸及其谐波对心跳的干扰的影响，如：1)基于hough变换的呼吸谐波抑制算法，再应用基于线性预测编码(lpc)的高分辨率心跳推导算法来实现对心跳的高分辨率估计；2)通过心跳信号的增强方法，在通过自相关函数获取峰值信息，进而获得心率估计，但是上述方法均难以有效地去除呼吸谐波对心跳分量的干扰。因此，如何解决呼吸谐波严重干扰心率估计，以及由于其频率与心跳信号重叠难以分离的问题，提高心率的估计的精度与抗干扰能力，成为目前亟需解决的技术问题。


技术实现要素：

5.针对上述现有技术存在的问题，本发明提供一种基于心跳二次谐波信号重构来估计心率的方法，能够解决呼吸谐波及其频率对心率的干扰，提高心率估计的精度及抗干扰能力。
6.为实现上述目的，本发明提供一种基于心跳二次谐波重构的心率估计方法，包括以下步骤：
7.s1利用毫米波雷达发射连续调频波，通过回波信号与发射信号混频后得到含有人体胸腔的位移信息的中频信号，中频信号以距离维*慢时间维的格式被存储于一个数据矩阵中；
8.s2对步骤s1中获得的中频信号进行预处理，得到目标距离单元中包含胸腔位移信息的瞬时相位差信号，即胸壁位移信号；
9.s3对胸壁位移信号进行分解与二次谐波频段提取：首先对胸壁位移信号进行小波包分解，得到不同频段的子信号，然后选择对1.6hz-4hz的心跳二次谐波频段的子信号进行提取；
10.s4对步骤s3中提取的二次谐波信号加权重构：分别计算出步骤s3提取出的心跳二次谐波信号频段各个子信号的能量，并算出各个子信号的能量与二次谐波频段所有子信号能量总和的比值，该比值作为心跳二次谐波重构时的子信号的加权系数，对各子信号进行加权并相加，完成心跳二次谐波信号重构；
11.s5功率谱估计：通过加权交叠平均法将步骤s4中重构的二次谐波信号数据分段处理，相邻的每段进行重叠，然后对数据加窗进行功率谱估计，最终得到心率估计值。
12.与现有技术相比，本发明通过避开在0.82hz-2hz获取由心跳引起的部分胸壁位移信号hr(t)，在二次谐波频段(1.6-4hz)中提取由心跳引起的部分胸壁位移信号hr(t)的二次谐波信号，有效地去除了呼吸谐波的干扰，使得心率估计更为准确；同时通过分段重叠和加窗的方法进行功率谱估计，最终得到心率估计值，减少了心跳二次谐波方差，解决了呼吸频率对心率的干扰，进一步提高了心跳估计的精度及抗干扰能力。
附图说明
13.图1是本发明一个实施例的流程图；
14.图2是本发明一个实施例所得的心率实时曲线与真实值的拟合结果。
具体实施方式
15.下面详细描述本发明的实施例，所述实施例的示例在附图中示出，其中自始至终相同或类似的标号表示相同或类似的元件或具有相同或类似功能的元件。下面通过参考附图描述的实施例是示例性的，旨在用于解释本发明，而不能理解为对本发明的限制。
16.下文的公开提供了许多不同的实施例或例子用来实现本发明的不同结构。为了简化本发明的公开，下文中对特定例子的部件和设置进行描述。当然，它们仅仅为示例，并且目的不在于限制本发明。此外，本发明可以在不同例子中重复参考数字和/或字母。这种重复是为了简化和清楚的目的，其本身不指示所讨论各种实施例和/或设置之间的关系。此外，本发明提供了的各种特定的工艺和材料的例子，但是本领域普通技术人员可以意识到其他工艺的可应用性和/或其他材料的使用。
17.下面参考图1描述根据本发明实施例的一种基于心跳二次谐波重构的心率估计方法，包括以下步骤：
18.s1利用毫米波雷达发射连续调频波，通过回波信号与发射信号混频后得到含有人体胸腔的位移信息的中频信号，中频信号以距离维*慢时间维的格式被存储于一个数据矩阵中，具体步骤如下：
19.s1.1毫米波雷达发射连续调频信号，调频脉冲信号发射过程中被其发射路径上的物体阻挡后进行反射，其中，线性调频连续波雷达的发射信号表示为：
[0020][0021]
式中：a
t
为发射功率，φ(t)是相位噪声，b为线性调频信号的带宽，fc为初始频率，t为发射时间，tc为线性调频脉冲的宽度；
[0022]
反射信号表示为：
[0023]
xr(t)＝αa
t
cos(2πfc(t-td)+πk(t-td)2+φ(t-td))
ꢀꢀꢀ
(2)
[0024]
式中：td是发射信号与反射信号的时间延迟，r(t)是雷达与目标的径向距离，这里的目标是毫米波雷达发射路径上的物体，α是信号传播过程中幅度的衰减，t为发射时间，
[0025]
s1.2将发射信号与接收信号进行差频得到中频信号w(t)：
[0026][0027]
其中，
[0028][0029]
式中：fb为中频信号频率，φb(t)为中频信号相位，δφ(t)为相位噪声，ar为接收信号的幅度；c为光速；
[0030]
由于在实际应用中，通常很小，可以忽略不计，所以，令结合公式(4)可得公式：
[0031][0032]
式中：中频信号相位与雷达与目标的径向距离r(t)成正相关变化，当目标运动位移为δr时，前后相位相差则可推出公式(7)：
[0033][0034]
由此，得到了包含有胸腔位移信息的中频信号w(t)，从而得到中频信号频率fb和中频信号相位φb(t)的信息，进而获得被测目标的距离信息和胸腔振动频率信息。然而，在实际被测环境中，存在着静态物体等其他因素干扰目标，使得中频信号w(t)中的中频信号频率fb分量不同，为得到测量目标的胸腔位移信号s(t)，需要进一步对中频信号w(t)进行预处理。
[0035]
s2对步骤s1中获得的中频信号进行预处理，得到目标距离单元中包含胸腔位移信息的瞬时相位差信号，即胸壁位移信号；
[0036]
s2.1距离傅里叶变换：对每个啁啾采用快速傅里叶变换(fft)来获得胸壁的目标距离仓，在慢时间维对中频信号进行fft处理之前，采用汉宁窗对慢时信号进行处理，该窗可以对慢时信号进行旁瓣抑制，同时保留振动相位信息；
[0037]
所述目标距离仓解释如下：已知毫米波雷达与受试者的真正距离，因此，可以通过
距离傅立叶变换分离不同距离的信号，根据毫米波雷达与受试者的真正距离，在毫米波雷达采集的所有数据中筛选出受试者所在距离的、含有人体胸腔的位移信息的中频信号，称为目标距离仓。为进一步理解目标距离仓，进一步举例说明：如毫米波雷达与受试者的距离为2m，则毫米波雷达在采集的所有数据中筛选出距离毫米波雷达2m处含有人体胸腔的位移信息的中频信号，该数据称为目标距离仓。
[0038]
s2.2直流偏置补偿：在雷达进行体征监测过程中，杂波反射、器件电路不稳定等会导致直流输出，在目标距离仓中，通过减去慢时间窗的平均值来消除直流偏置对心跳估计的影响，获得优化的目标距离仓；；
[0039]
s2.3相位提取与解缠绕：将优化的目标距离仓采用复信号表示，该复信号的实部和虚部包含了目标距离仓的相位信息，利用实部和虚进行反正切运算得到目标距离仓的相位值。由于在进行反正切运算获取目前的相位值的过程中会产生相位缠绕问题，因此，通过解缠绕算法来对目标距离仓的相位数据进行校正以获取准确的相位信息；
[0040]
s2.4相位差计算：在获取准确的目标距离仓的相位值之后，通过对相邻相位进行连续展开计算相位差，消除相位漂移；
[0041]
s2.5噪声去除：采用移动平均滤波器对相位差信号进行处理，以消除部分噪声的影响，得到去除噪声杂波的目标距离单元中包含胸腔位移信息的瞬时相位差信号s(t)，即胸壁位移信号s(t)；
[0042]
s3胸壁位移信号s(t)分解与二次谐波频段提取：首先对胸壁位s(t)移信号进行小波包分解，得到不同频段的子信号，考虑到心跳频段(0.8hz-2hz)存在呼吸二次、三次谐波的强干扰，本实施例选择对1.6hz-4hz的心跳二次谐波频段的子信号进行提取，具体步骤如下：
[0043]
s3.1步骤s2得到的胸壁位移信号s(t)包含由心跳引起的部分胸壁位移信号hr(t)(以下简称“hr(t)”)，、呼吸引起的部分胸壁位移信号br(t)(以以下简称“br(t)”)以及其余噪声干扰n(t)，经过步骤s2.1-s2.5处理后，n(t)已经可忽略不计。
[0044]
为了分离不同频段的信号，首先对胸壁位移信号s(t)进行小波包分解，设小波包分解层数为n，胸壁位移信号s(t)的采样率为fs，通过小波包分解，胸壁位移信号s(t)将会被分解成2n个子信号，且根据公式δf＝fs/2n各个子信号带宽都为δf，按照频段大小从低到高排列，则各子信号的频段分别为[0，δf]、[δf，2δf]
……
、[(2
n-1)δf，2nδf]，同时从低到高的频段排列顺序，各子信号分别记为s1(t)、s2(t)、
……
、
[0045]
hr(t)的频段是0.8hz-2hz，br(t)的频段是0.1-0.5hz，虽然在0.8hz-2hz频段提取hr(t)不会受到br(t)的干扰，但是在br(t)的二次谐波频段(0.2-1hz)以及三次谐波频段(0.3hz-1.5hz)会和hr(t)局部重叠进而降低提取hr(t)的准确度，因此，需要避开在0.8hz-2hz获取hr(t)，转向在hr(t)二次谐波频段(1.6hz-4hz)中提取hr(t)的二次谐波信号，将hr(t)的二次谐波信号记为hr2(t)。由于三次谐波往后呼吸谐波强度干扰很小，所以不考虑呼吸三次谐波往后的频段。因此，本实施例通过避开在0.82hz-2hz获取hr(t)，而是在二次谐波频段(1.6-4hz)中提取hr(t)的二次谐波信号，有效地去除了呼吸谐波的干扰，使得心率估计更为准确。
[0046]
s3.2在2n个子信号中找出在心跳二次谐波信号hr2(t)频段的子信号，设子信号分别为sm(t)、
……
sn(t)(0＜m＜n＜
2n
)，则可得到表达式：
[0047]
hr2(t)＝sm(t)+s
m+1
+
……
+sn(t)
ꢀꢀꢀ
(8)
[0048]
s4二次谐波信号能量加权重构：分别计算出步骤s3提取出的心跳二次谐波信号hr2(t)频段各个子信号的能量，并算出各个子信号的能量与二次谐波频段hr2(t)所有子信号能量总和的比值，该比值作为心跳二次谐波重构时的子信号的加权系数，对各子信号进行加权并相加，完成心跳二次谐波信号重构，具体步骤如下：
[0049]
s4.1通过信号能量公式得到由心跳引起的部分胸壁位移信号hr(t)的二次谐波频段各子信号的能量，分别记为am(t)、
……
an(t)(0＜m＜n＜
2n
)；
[0050]
式中，x(t)为信号，a(t)是x(t)的能量；
[0051]
s4.2：根据各个子信号的能量与二次谐波信号hr2(t)频段所有子信号能量和的比值，在公式(8)的基础上对各子信号进行加权并相加，加权相加处理之后的心跳二次谐波信号记为d(t)，d(t)表示为：
[0052][0053]
经过能量加权重构后，二次谐波信号频段中二次谐波信号hr2(t)得到进一步加强，其余位于二次谐波频段的干扰信号被大大削弱，在接下来的功率谱估计中，计算的hr2(t)频率更加准确。
[0054]
s5功率谱估计：通过加权交叠平均法(welch)将重构的二次谐波信号数据分段处理，相邻的每段进行重叠，然后对数据加窗进行功率谱估计，最终得到心率估计值，具体步骤如下：
[0055]
s5.1将加权相加处理之后的心跳二次谐波信号d(t)分成l小段，其中，每小段x个点，相邻小段交叠y/2个点，段数l为：
[0056][0057]
s5.2对各小段加同样的平滑窗w(n)后求傅里叶变换得到第i小段的功率谱ti(e
jω
)：
[0058][0059]
s5.3根据步骤s5.2所得数据求所有段的平均功率谱p
x
(e
jw
)：
[0060][0061]
其中，为归一化因子。
[0062]
s5.4根据步骤s5.3得到的平均功率谱计算出心跳信号的频率，具体方法如下：
[0063]
设功率谱谱峰对应的频率为f2h，则将心跳信号频率的估计值记为f2h/2，
[0064]
本发明步骤s5通过分段重叠和加窗的方法，本步骤的功率谱估计结果是逐步一致估计，并且减少了心跳二次谐波方差，使得心跳估计更为准确。
[0065]
本实施例方法主要应用在室内场景，如办公室中，本实施例毫米波雷达信号为iwr1843，毫米波雷达被放置在一个固定的支架上，受试者正对毫米波雷达，离雷达波雷达的距离在0.6m～2.0m的范围内。在测试过程中，受试者保持身体静止，进行正常的呼吸，毫米波雷达配备3个tx天线和4个rx天线，运行频率范围为76ghz-81ghz，毫米波雷达设置为
tdm-mimo模式，将iwr1843与dca1000evm连接起来。开启毫米波雷达，毫米波雷达开始发射连续调频波并接收反射信号，将混频得到的中频信号通过以太网发送至上位机，上位机采用本实施例的方法高速进行中频信号数据处理，最终实现准确的心率估计。为了验证本实施例方法的准确性，实验人员在进行方法的测试实验时利用心电芯片ads1292同步心电数据采集，作为地面真实值，以下为实验人员进行的一分钟数据采集，如图2所示为本实施例方法所得的心率实时曲线与真实值的拟合结果示意图，其中，实线为真实值，虚线为心率实时曲线，通过图2可以明显的看出通过本发明方法测得的心率实时曲线与真实值基本相同，因此，本发明方法心率估计具有较高的准确率。
[0066]
根据本发明实施例的基于心跳二次谐波重构的心率估计方法的其他构成例如毫米波雷达等以及操作对于本领域普通技术人员而言都是已知的，这里不再详细描述。
[0067]
在本发明的描述中，需要理解的是，术语“中心”、“纵向”、“横向”、“长度”、“宽度”、“厚度”、“上”、“下”、“前”、“后”、“左”、“右”、“竖直”、“水平”、“顶”、“底”、“内”、“外”、“顺时针”、“逆时针”、“轴向”、“径向”、“周向”等指示的方位或位置关系为基于附图所示的方位或位置关系，仅是为了便于描述本发明和简化描述，而不是指示或暗示所指的装置或元件必须具有特定的方位、以特定的方位构造和操作，因此不能理解为对本发明的限制。
[0068]
此外，术语“第一”、“第二”仅用于描述目的，而不能理解为指示或暗示相对重要性或者隐含指明所指示的技术特征的数量。由此，限定有“第一”、“第二”的特征可以明示或者隐含地包括一个或者更多个该特征。在本发明的描述中，“多个”的含义是两个或两个以上，除非另有明确具体的限定。
[0069]
在本发明中，除非另有明确的规定和限定，术语“安装”、“相连”、“连接”、“固定”等术语应做广义理解，例如，可以是固定连接，也可以是可拆卸连接，或成一体；可以是机械连接，也可以是电连接，还可以是通信；可以是直接相连，也可以通过中间媒介间接相连，可以是两个元件内部的连通或两个元件的相互作用关系。对于本领域的普通技术人员而言，可以根据具体情况理解上述术语在本发明中的具体含义。
[0070]
在本发明中，除非另有明确的规定和限定，第一特征在第二特征“上”或“下”可以是第一和第二特征直接接触，或第一和第二特征通过中间媒介间接接触。而且，第一特征在第二特征“之上”、“上方”和“上面”可是第一特征在第二特征正上方或斜上方，或仅仅表示第一特征水平高度高于第二特征。第一特征在第二特征“之下”、“下方”和“下面”可以是第一特征在第二特征正下方或斜下方，或仅仅表示第一特征水平高度小于第二特征。
[0071]
在本说明书的描述中，参考术语“一个实施例”、“一些实施例”、“示例”、“具体示例”、或“一些示例”等的描述意指结合该实施例或示例描述的具体特征、结构、材料或者特点包含于本发明的至少一个实施例或示例中。在本说明书中，对上述术语的示意性表述不必须针对的是相同的实施例或示例。而且，描述的具体特征、结构、材料或者特点可以在任一个或多个实施例或示例中以合适的方式结合。此外，在不相互矛盾的情况下，本领域的技术人员可以将本说明书中描述的不同实施例或示例以及不同实施例或示例的特征进行结合和组合。
[0072]
尽管已经示出和描述了本发明的实施例，本领域的普通技术人员可以理解：在不脱离本发明的原理和宗旨的情况下可以对这些实施例进行多种变化、修改、替换和变型，本发明的范围由权利要求及其等同物限定。