基于LFMCW雷达的心率估计算法的制作方法

本发明涉及雷达检测领域，具体涉及一种基于lfmcw雷达的心率估计算法。

背景技术：

随着生活水平的提高，人们对自身的健康状况更加关心，生命体征的监测技术的要求也越来越高，非接触生命体征监测的技术也受到了众多学者的关注。超声波、wifi、相机和雷达等都在非接触生命体征监测中得到应用，其中超声波、wifi和雷达都是基于多普勒效应。相比超声波设备功率高，噪声大，wifi的信号处理不方便等缺点，雷达在非接触生命体征监测中受到广大学者的青睐。生物雷达的研究也分为三种，分别是单音连续波雷达(continuouswave，cw)，脉冲超宽带雷达(impulseradioultra-wideband，ir-uwb)和调频连续波雷达(frequency-modulatedcontinuouswave，fmcw)。

cw雷达的工作原理是发射天线发射一个单音信号，通过接收物体的反射信号来获取运动目标的相位信息，从而获得多普勒信息和距离信息，这种雷达没有距离分辨率，所以距离隔离性差，检测微小运动的能力也差。多普勒雷达还存在直流偏移和i/q不平衡的问题，这种系统上的缺陷是难以弥补的。

ir-uwb雷达发射的带宽较大，因此可以实现非常高的距离分辨率。其特点是发射天线发射占空比很窄的脉冲。由于占空比低，信号的能量和信噪比(snr)也低，距离精度因此也降低。

基于fft的心率估计和基于levenberg-marquardt(lm)拟合的心率估计是当前lfmcw雷达测量心率的主要方法，这两种方法存在稳定性差，精准度低等现状。

请看图1a，为原始的微动信号的幅度-时间图示。基于该微动信号，分别进行基于fft的心率估计，以及进行基于levenberg-marquardt(lm)拟合的心率估计。

请看图1b，为基于fft进行心率估计时的频谱图。在图1b所示的频谱图中，当0.8hz～2hz对应的通带内存在干扰时，这种方法无法解决噪声干扰问题，准确率低。请看图1c，为图1b中0～2hz的频段的放大，可以看出，心跳为88的志愿者的估计结果为79，与实际使用接触式心律仪检测后获得的实时结果相对比，准确率仅为89.77％。

除此之外，这种方法还需要较高的采样率或较长采样时间。当采样率为20hz时，需要采样点12000点，即需要采集10分钟数据才能实现心跳分辨率为1跳。

请看图1d，为基于levenberg-marquardt(lm)拟合的心率估计法进行降噪之后获取的频谱图，其中图1e为对图1d中的0.8～2hz的频段进行放大后获取到的图谱。可以看出，心跳为88的志愿者的估计结果也为79，与实际使用接触式心律仪检测后获得的实时结果相对比，准确率也仅为89.77％。该方法能够在一定程度上解决噪声问题，但准确率仍然不高，微动干扰下性能提升不大。

技术实现要素：

本发明的目的在于提供一种基于lfmcw雷达的心率估计算法，具有更高的稳定性，以及更高的精准度。

为了解决上述技术问题，以下提供了一种基于lfmcw雷达的心率估计算法，包括以下步骤：提供lfmcw雷达，所述lfmcw雷达能够进行检测信号的收发；使用所述lfmcw雷达向目标物发送n个检测信号，并对应接收n个反射信号；对各个反射信号和与各个反射信号对应的检测信号进行混频，获取n个中频信号，形成原始数据矩阵；对所述原始数据矩阵进行傅里叶变换，获取距离矩阵；获取所述目标物在所述距离矩阵中的下标；根据所述目标物在所述距离矩阵中的下标，获取所述目标物的原始相位信号；根据所述原始相位信号获取到所述目标物的微动信号；从所述微动信号中筛选心跳信号，并进行心率估计。

可选的，构成n×m的原始数据矩阵m[nslow,mfast]，其中nslow＝1，2，……n，mfast＝1，2，……m，其中m是对每个检测信号进行采样时的采样点数。

可选的，对所述原始数据矩阵进行傅里叶变换，获取距离矩阵时，对所述原始数据矩阵进行快时间维度的傅里叶变换，得到所述距离矩阵rp[nslow,mfast]。

可选的，计算所述目标物的位置所在的检测信号在所述距离矩阵中的下标mpeople，所述目标物与所述lfmcw雷达之间的距离为其中为mpeople的频率，

可选的，根据所述目标物在所述距离矩阵中的下标，获取所述目标物的原始相位信号时，所述原始相位信号为swrap(nslow)＝rp[nslow,mpeople]，为慢时间维度的相位。

可选的，在根据所述原始相位信号获取到所述目标物的微动信号前，还包括以下步骤：对所述目标物的原始相位信号的相位跳变进行修正。

可选的，根据修正了相位跳变后的所述目标物的原始相位信号，获取所述微动信号。

可选的，所述修正了相位跳变后的所述目标物的原始相位信号为：

sunwrap(nslow)＝unwrap(swrap(nslow))；

所述微动信号为：

可选的，使用基于pe的meemd滤波器从所述微动信号中筛选获得心跳信号。

可选的，通过峰值检测算法获得心率估计值hr。

可选的，所述lfmcw雷达发送的检测信号为：

其中fc是中心频率，是线性调频斜率，b是lfmcw雷达的带宽，τchirp是所述检测信号斜率上升的时间，是所述检测信号的初始相位，r(τ)为目标物与lfmcw雷达的距离。

可选的，经过距离为r(τ)的目标物反射后，所述lfmcw雷达接收到的反射信号如下所示：

其中srx(t)为所述反射信号，σ表示所述反射信号的幅度，由反射物体的雷达散射截面和传播损耗共同决定；将所述反射信号和检测信号进行混频，获取的所述中频信号为：

其中，sif(t)为所述中频信号，fif＝4πγr(τ)/c与距离r(τ)成正比，且4πfcr(τ)/c为慢时间相位，为残余相位；

所述距离r(τ)满足：

可选的，对所述中频信号进行快时间维度的傅里叶变换，以获取距离频谱图，从而获取距离为r(τ)的目标物的原始相位信号。

可选的，使用基于pe的meemd滤波器从所述微动信号中筛选获得心跳信号时，包括以下步骤：步骤1：在所述微动信号中添加零均值白噪声信号，获得第一信号和第二信号，所述第一信号为所述微动信号与所述零均值白噪声信号的和，所述第二信号为所述微动信号与所述零均值白噪声信号的差；步骤2：分别对所述第一信号和第二信号进行经验模态分解，得到第m阶imf分量，分别为和其中i＝1，2，……ne，，且进行ne次经验模态分解后得到的平均第m阶imf分量im(t)为：步骤3：设置pe阈值thmeemd，计算平均第m阶imf分量im(t)的归一化排列熵其中k代表嵌入维数为k的pe；步骤4：若则判定所述平均第m阶imf分量im(t)为平稳信号，所述平均第m阶imf分量im(t)包含心跳信号，并进入步骤5，否则，判定所述平均第m阶imf分量im(t)为噪声信号或者干扰，并进行所述微动信号与所述平均第m阶imf分量im(t)的相减运算，并重新开始所述步骤1；步骤5：对所述微动信号进行经验模态分解，设置pe区间，且若所述的归一化排列熵位于所述pe区间内，则判断所述平均第m阶imf分量im(t)为心跳信号，否则判定为杂波。

可选的，所述pe阈值thmeemd为0.6，所述pe区间为[0.31，0.44]。

本发明的基于lfmcw雷达的心率估计算法从所述微动信号中筛选心跳信号，并进行心率估计，能够去除微动信号中的噪声和干扰，使得最终获取到的心跳信号更加准确，具有更高的稳定性，以及更高的精准。

附图说明

图1a至图1f为采用各种方法进行心率估计时的示意图。

图2为本发明的一种具体实施方式中的基于lfmcw雷达的心率估计算法步骤流程示意图。

图3a为采样率为20sps时50至120beat/min的pe区间的示意图。

图3b为采样率为4000sps时50至120beat/min的pe区间的示意图。

图4为本发明的一种具体实施方式中的使用基于pe的meemd滤波器从所述微动信号中筛选获得心跳信号时的步骤流程示意图。

具体实施方式

以下结合附图和具体实施方式对本发明提出的一种基于lfmcw雷达的心率估计算法作进一步详细说明。

请参阅图2，为本发明的一种具体实施方式中的基于lfmcw雷达的心率估计算法步骤流程示意图。

在该具体实施方式中，提供了一种基于lfmcw(linearfrequencymodulatedcontinuouswave)雷达的心率估计算法，包括以下步骤：提供lfmcw雷达，所述lfmcw雷达能够进行检测信号的收发；使用所述lfmcw雷达向目标物发送n个检测信号，并对应接收n个反射信号；对各个反射信号和与各个反射信号对应的检测信号进行混频，获取n个中频信号，形成原始数据矩阵；对所述原始数据矩阵进行傅里叶变换，获取距离矩阵；获取所述目标物在所述距离矩阵中的下标；根据所述目标物在所述距离矩阵中的下标，获取所述目标物的原始相位信号；根据所述原始相位信号获取到所述目标物的微动信号；从所述微动信号中筛选心跳信号，并进行心率估计。

在该具体实施方式中，使用lfmcw雷达结合心跳筛选的方法来实现人体识别和心跳检测。所述lfmcw雷达综合了cw雷达和ir-uwb雷达的优点，既能够基于相位测量距离，以获得更精确的距离测量结果，有更好的距离探测精度，又能够随着带宽增加，获取更精细的距离分辨率。并且，由于lfmcw雷达发射的是连续波，因此可以传输更多的能量，信噪比也得到提高，因此距离精度进一步提高。

在该具体实施方式中，所述lfmcw雷达为大带宽的lfmcw雷达，工作在77ghz频段，发送的检测信号是毫米波，可以通过设计合适的参数，来获取高距离分辨率，从而可以隔离被探测的目标物周围的杂波，并检测到目标物的微小运动，达到更加精确的心率估计效果。

在一种具体实施方式中，所述lfmcw雷达的带宽可高达4ghz。实际上也可根据需要设置所述lfmcw雷达的带宽范围。

在一种具体实施方式中，构成n×m的原始数据矩阵m[nslow,mfast]，其中nslow＝1，2，……n，mfast＝1，2，……m，其中m是对每个检测信号进行采样时的采样点数。

在一种具体实施方式中，对所述原始数据矩阵进行傅里叶变换，获取距离矩阵时，对所述原始数据矩阵进行快时间维度的傅里叶变换，得到所述距离矩阵rp[nslow,mfast]，距离矩阵就是多个检测信号的距离谱图的堆叠。

在一种具体实施方式中，计算所述目标物的位置所在的检测信号在所述距离矩阵中的下标mpeople，所述目标物与所述lfmcw雷达之间的距离为其中为mpeople的频率，

在一种具体实施方式中，根据所述目标物在所述距离矩阵中的下标，获取所述目标物的原始相位信号时，所述原始相位信号为swrap(nslow)＝rp[nslow,mpeople]，为慢时间维度的相位。

在该具体实施方式中，获取下标的方法是通过比较距离矩阵中每一行的峰值的下标，求取n个检测信号的下标之后，出现最多的下标就是最终的下标，一般情况下n个检测信号的下标是相同的。

在一种具体实施方式中，所述检测信号为线性调频信号(chirp),线性调频信号可以对载波频率进行调制以增加信号的发射带宽并在接收时实现脉冲压缩。由于线性调频信号具有较高的距离分辨力，当在速度上无法区分多目标时，可以通过增加目标距离测试解决多目标的分辨问题。同时在抗干扰方面，线性调频信号可以在距离上区分干扰和目标，因而可以有效地对抗拖曳式干扰，这使得线性调频信号在雷达波形设计中得到了广泛的应用。

在一种具体实施方式中，在根据所述原始相位信号获取到所述目标物的微动信号前，还包括以下步骤：对所述目标物的原始相位信号的相位跳变进行修正。

在一种具体实施方式中，根据修正了相位跳变后的所述目标物的原始相位信号，获取所述微动信号。

在一种具体实施方式中，所述修正了相位跳变后的所述目标物的原始相位信号为：

sunwrap(nslow)＝unwrap(swrap(nslow))；

所述微动信号为：

需要注意的是，本领域的技术人员能够毫无疑义的知道该如何进行unwrap操作来纠正相位跳变，在此不多做赘述。

在一种具体实施方式中，使用基于pe(permutationentropy，排列熵)的meemd(改进的经验模态分解)滤波器从所述微动信号中筛选获得心跳信号。排列熵是一种一维时间序列的复杂度的平均熵度量方式。较小的噪声本质上不会改变含噪信号的复杂度，所以可以认为任意现实世界时间序列计算排列熵。由于排列熵具有快速和健壮的特点，因此在存在大量数据集且没有时间进行预处理和参数微调时是可取的。

在该具体实施方式中，在筛选之前需要通过meemd算法对微动信号进行分解，使得分解出来的信号中包括心跳信号，才能通过排列熵pe去将心跳信号筛选出来。

pe的计算方法如下：

对序列{x(i)，i＝1，2，……n}进行相空间重构得到下式所示序列x：

其中m是嵌入维数，λ是时间延迟，x(k)＝{x(k)，x(k+λ)，……x(k+(m－1)λ)}，k＝1，2，……n－(m－1)λ，将x(k)按升序排列，如下式所示：

xsort(k)＝{x(i1)，……x(iq)，……x(im)}；

其中iq∈{k，k+λ，……，k+(m-1)λ}为下标，x(i1)<……<x(iq)<……<x(im)，表示从x(i1)到x(iq)再到x(im)是一个递增数列，xsort(k)为升序序列，遇到等大的值，则根据下标升序排序。序列(i1，……iq，……im)记为indexsort(k)＝j，j∈{1，2，……m！}，代表(i1，……iq，……im)这个下标序列在1～m的全排列矩阵中的第j行。下式(1)定义了排列熵。下式(2)定义了归一化排列熵。

其中pj为全排列矩阵中的第j行出现的频率，ln^m！为排列熵为hp(m)的最大值。hp_norm(m)越大说明序列越随机。

在采样率确定的情况下，范围为50至120beat/min的心跳信号的归一化排列熵hp_norm(m)与心率呈正比，这与pe的定义相符。与此同时，采样率的降低，正比关系整体趋势不变，虽然在局部会有波动，但是排列熵的区间用于筛选信号仍是有效的。如图3a所示，当采样率为20sps(samplepersecond，每秒采样次数)时，pe值线性分布在0.31～0.44之间。如图3b所示，当采样率升高到4000sps，周期现象变强，pe值降低，线性分布在0.107～0.111之间。可以根据这个性质筛选meemd分解之后的分量，获得心跳信号。

在一种具体实施方式中，通过峰值检测算法获得心率估计值hr。峰值检测算法可由本领域的技术人员轻易获知，在此不多做赘述。

在一种具体实施方式中，所述lfmcw雷达发送的检测信号为：

其中fc是中心频率，是线性调频斜率，b是lfmcw雷达的带宽，τchirp是所述检测信号斜率上升的时间，是所述检测信号的初始相位，r(τ)为目标物与lfmcw雷达的距离。

在一种具体实施方式中，经过距离为r(τ)的目标物反射后，所述lfmcw雷达接收到的反射信号如下所示：

其中srx(t)为所述反射信号，σ表示所述反射信号的幅度，由反射物体的雷达散射截面和传播损耗共同决定；将所述反射信号和检测信号进行混频，获取的所述中频信号为：

其中，sif(t)为所述中频信号，fif＝4πγr(τ)/c与距离r(τ)成正比，且4πfcr(τ)/c为慢时间相位，为残余相位；

所述距离r(τ)满足：

在一种具体实施方式中，对所述中频信号进行快时间维度的傅里叶变换，以获取距离频谱图，从而获取距离为r(τ)的目标物的原始相位信号。

请参阅图4，为本发明的一种具体实施方式中的使用基于pe的meemd滤波器从所述微动信号中筛选获得心跳信号时的步骤流程示意图。在该具体实施方式中，使用基于pe的meemd滤波器从所述微动信号中筛选获得心跳信号时，包括以下步骤：步骤1：在所述微动信号中添加零均值白噪声信号，获得第一信号和第二信号，所述第一信号为所述微动信号与所述零均值白噪声信号的和，所述第二信号为所述微动信号与所述零均值白噪声信号的差，需要注意的是，图4中使用s(t)来描述该微动信号；步骤2：分别对所述第一信号和第二信号进行经验模态分解，得到第m阶imf分量，分别为和其中i＝1，2，……ne，且进行ne次经验模态分解后得到的平均第m阶imf分量im(t)为：步骤3：设置pe阈值thmeemd，计算平均第m阶imf分量im(t)的归一化排列熵其中k代表嵌入维数为k的pe；步骤4：若则判定所述平均第m阶imf分量im(t)为平稳信号，所述平均第m阶imf分量im(t)包含心跳信号，并进入步骤5，否则，判定所述平均第m阶imf分量im(t)为噪声信号或者干扰，并进行所述微动信号与所述平均第m阶imf分量im(t)的相减运算，并重新开始所述步骤1；步骤5：对所述微动信号进行经验模态分解，设置pe区间，且若所述的归一化排列熵位于所述pe区间内，则判断所述平均第m阶imf分量im(t)为心跳信号，否则判定为杂波。

在一种具体实施方式中，所述pe阈值thmeemd为0.6，所述pe区间为[0.31，0.44]。

在该具体实施方式中，获取距离矩阵中的下标mpeople，通过速度-距离谱图识别被测人体，再通过仿真实验获得心跳区间50至120次/分钟时对应的pe区间为[0.31，0.44]。通过该pe区间，可以有效消除消除微动信号中的干扰和噪声。

如图1f所示，采用该具体实施方式中的基于pe的meemd心跳筛选算法进行降噪后获取的微动信号的频谱图。可以看出，图1f中的微动信号含有87个周期，这与实际使用接触式心律仪检测后获得的实时结果接近，准确率达到98.86％。

此外，采用该具体实施方式中的基于pe的meemd心跳筛选算法，采样率为20hz时，仅需采样1200点便可估计心跳，极大的减少计算量，提高了心跳测量的实时性和准确性。

以上所述仅是本发明的优选实施方式，应当指出，对于本技术领域的普通技术人员，在不脱离本发明原理的前提下，还可以做出若干改进和润饰，这些改进和润饰也应视为本发明的保护范围。