CN105232026A - 一种非接触式生命体征检测系统的心跳频率检测算法 - Google Patents

Abstract

本发明提供一种非接触式生命体征检测系统的心跳频率检测算法，步骤包括：对连续波多普勒雷达前端输出的I、Q两路信号分别进行带通滤波，然后采用圆心估计算法进行有用直流分量恢复，采用完全聚类经验模态分解算法将心跳信号与人体抖动信号、呼吸信号和环境干扰噪声分离，求各分解信号的希尔伯特边际谱并进行峰值检测，根据边际谱的谱峰位置和谱峰附近的能量集中度找出与心跳信号对应的希尔伯特边际谱，并且根据谱峰位置得到心跳频率信息。本发明可在人体不稳定、环境干扰噪声大的条件下有效提取心跳信号并且获得准确心跳频率信息，以提高非接触式生命体征检测系统的抗干扰性，满足医务人员对心跳频率检测准确性的需求。

Description

一种非接触式生命体征检测系统的心跳频率检测算法

技术领域

本发明涉及一种基于雷达技术的非接触式生命体征检测系统的信号处理算法，尤其是心跳频率检测算法。

背景技术

非接触式生命特征检测技术将雷达技术和生物医学工程技术融为一体，利用电磁波对人体的呼吸和心跳等生命特征进行非接触式探测。电磁波具有非金属介质穿透特性，因此非接触式生命体征监测系统能够穿透衣服、被褥检测生命信号，不需要与人体皮肤直接接触，这一优势使其在婴儿、有睡眠障碍的病人、大面积烧伤病人、精神病人等特殊群体的生命体征监测中得到应用。此外，经过特殊设计的生物雷达还能穿透砖墙、木板等障碍物探测到人体的微动信号，这使得其在救灾、反恐、军事等领域也有广泛的应用前景。

心跳频率检测的准确度是非接触式生命特征检测系统的最重要的指标之一，由于心跳信号能量很微弱且具有非平稳特性，所以心跳信号很容易淹没在雷达的噪声和杂波中。此外，由于人体在监测过程中很有可能进行随机运动，如人体的抖动、翻身等，这些随机运动产生的扰动干扰信号都具有非平稳特点，其频率很有可能在心跳频率范围内，因此传统滤波方法不能将这些扰动干扰滤除，而且由于传统FFT方法集中度较低、分辨力较差，因此采用FFT测频法不能准确测得心跳频率。

发明内容

为了克服上述现有技术的不足，本发明提供了一种非接触式生命体征检测系统的心跳频率检测算法，该算法能将心跳信号与人体抖动信号、呼吸信号和环境干扰噪声分离，而且能准确提取心跳信号的频率信息。

本发明采用的技术方案如下：

一种用于非接触式生命体征检测系统的心跳频率检测算法，包括如下步骤：

步骤1，对非接触式生命体征检测系统的连续波多普勒雷达前端输出的I、Q两路信号分别进行带通滤波；其中，带通滤波器的低截止频率为0.5Hz、高截止频率为5.2Hz；

步骤2，通过圆心估计算法恢复经过滤波的I、Q两路信号的有用直流分量；

步骤3，对步骤2得到的I、Q两路信号进行反正切解调；

其特征在于，得到解调信号之后的步骤如下：

步骤4，采用完全聚类经验模态分解算法对解调信号进行分解，得到一系列本征模函数；

步骤5，对得到的本征模函数分别做希尔伯特变换，并且求得各本征模函数对应的希尔伯特边际谱；

步骤6，在各本征模函数对应的希尔伯特边际谱上进行峰值检测；

步骤7，根据希尔伯特边际谱的谱峰位置和谱峰附近的能量集中度来判断心跳信号对应的希尔伯特边际谱；

步骤8，在心跳信号对应的希尔伯特边际谱上提取心跳信号频率信息。

进一步的，步骤4中采用完全聚类经验模态分解算法对解调信号进行分解，具体步骤如下：

步骤4.1，构造I个信号x[n]+ε₀wⁱ[n]，其中x[n]是反正切解调得到的信号，wⁱ[n],(i＝1,2,...,I)是分布为N(0,1)的白噪声，利用普通经验模态分解算法分别求解I个信号x[n]+ε₀wⁱ[n]的第一个本征模分量IMF₁ ⁱ[n]，并且计算： ${\mathrm{IMF}}_1\[n\]=\frac{1}{I}\underset{i=1}{\overset{I}{\Σ}}{\mathrm{IMF}}_1^i\[n\]={\stackrel{\‾}{IMF}}_1\[n\];$

步骤4.2，计算第一个剩余分量r₁[n]＝x[n]-IMF₁[n]；

步骤4.3，对r₁[n]+ε₁E₁(wⁱ[n]),i＝1,...,I进行经验模态分解，得到第一个本征模分量；其中，算子E_j(·)的功能是利用普通经验模态分解算法求一个信号的第j个模态，并由此得到： ${\mathrm{IMF}}_2\[n\]=\frac{1}{I}{\mathrm{\Σ}}_{i=1}^I{E}_1(r_1\[n\]+{\mathrm{\ϵ}}_1{E}_1(w^i\[n\]\left)\right);$

步骤4.4，根据k＝2,...,K，计算第k个剩余分量：r_k[n]＝r_(k-1)[n]-IMF_k[n]；

步骤4.5，根据r_k[n]+ε_kE_k(wⁱ[n]),i＝1,...,I，进行经验模态分解，得到第一个本征模分量，并由此得到 ${\mathrm{IMF}}_{(k+1)}\[n\]=\frac{1}{I}{\mathrm{\Σ}}_{i=1}^I{E}_1(r_k\[n\]+{\mathrm{\ϵ}}_k{E}_k(w^i\[n\]\left)\right);$

步骤4.6，重复步骤4.4到步骤4.5步，当剩余本征模分量不再满足可以分解的要求时，即剩余本征模分量的极值不足两个时，停止运算，最后的剩余本征模分量表示为： $R\[n\]=x\[n\]-{\mathrm{\Σ}}_{k=1}^K{\mathrm{IMF}}_k\[n\],$ K是总模态数；

至此，解调信号x[n]被分解为K个本征模函数IMF_k[n]和一个剩余分量R[n]。进一步的，步骤5中求解各本征模函数IMF_k[n]的希尔伯特边际谱，步骤如下：

步骤5.1，令u_k(n)＝IMF_k[n]，首先对u_k(n)做离散希尔伯特变换： $v_k\left(n\right)=\underset{m=0}{\overset{N-1}{\Σ}}h(n-m)u_k\left(n\right),$ 其中N为偶数，且 $h\left(n\right)=\frac{2}{N}{\sin }^2\left(\frac{\mathrm{\π}n}{2}\right)\cot \left(\frac{\mathrm{\π}n}{N}\right);$

步骤5.2，计算 $a_k\left(n\right)=\sqrt{u_k{\left(n\right)}^2+v_k{\left(n\right)}^2};$

步骤5.3，计算 ${\tilde{f}}_k\left(n\right)=\frac{u_k\left(n\right)\left(v_k\right(n+1)-v_k(n-1))-v_k\left(n\right)\left(u_k\right(n+1)-u_k(n-1))}{2t_s(u_k{\left(n\right)}^2+v_k{\left(n\right)}^2)};$ 对上式进行修正 $f_k\left(n\right)=f_{s}acsin({\tilde{f}}_k/f_s);$

步骤5.4，求希尔伯特谱H_k(f,n)＝H_k(f_k(n),n)＝a_k(n)；

步骤5.5，求得各本征模函数IMF_k[n]对应的希尔伯特边际谱： $h_k\left(f\right)=\underset{n=0}{\overset{N-1}{\Σ}}{H}_k(f,n).$

进一步的，步骤7中根据谱峰位置和谱峰附近的能量集中度来判断心跳信号对应的希尔伯特边际谱，具体方法如下：

步骤7.1，判断各边际谱h_k(f)的谱峰对应的频率f_peak,k是否在0.5Hz到3Hz之间；

步骤7.2，对于满足步骤7.1条件的边际谱，计算其在f_peak,k-0.1Hz到f_peak,k+0.1Hz之间各频率点的幅度的总和A_k，并计算其所有频率点的幅度总和B_k，求能量集中度E_k＝A_k/B_k；

步骤7.3，当E_k大于阈值α时，认为h_k(f)是心跳信号的边际谱；若有多个边际谱满足E_k大于阈值α的条件，则认为峰值幅度最大的边际谱是心跳信号的边际谱。

进一步的，步骤8中将心跳信号的希尔伯特边际谱的谱峰位置对应的频率作为心跳信号的频率。

进一步的，所述心跳信号的希尔伯特边际谱峰值对应的是频率为1.2402Hz。

与现有技术相比，本发明可在人体不稳定、环境干扰噪声大的条件下有效提取心跳信号并且获得准确心跳频率信息。

下面结合说明书附图对本发明做进一步描述。

附图说明

图1为本发明的非接触式生命体征检测系统的心跳频率检测算法流程图。图中的输入信号是连续波多普勒雷达前端送来I、Q两路信号。

图2为对I、Q两路信号进行带通滤波、圆心估计以及反正切解调后得到的解调信号。此信号中包含心跳信号以及干扰信号，其中13s之间23s有比较大的人体抖动干扰。

图3为对解调信号进行完全聚类经验模态分解后得到的一系列本征模函数的波形。

图4为第1到6个本征模函数对应的希尔伯特边际谱。(a)、本征模函数1的边际谱，(b)、本征模函数2的边际谱，(c)、本征模函数3的边际谱,(d)、本征模函数4的边际谱，(e)、本征模函数5的边际谱，(f)、本征模函数6的边际谱。

具体实施方式

非接触式生命体征检测系统是基于连续波多普勒雷达探测原理，雷达天线向人体胸腔位置发射电磁波，呼吸和心跳等引起胸腔的起伏会对电磁波的相位进行调制，胸腔反射的电磁波被雷达天线接收后经过雷达前端的处理变为I、Q两路信号输出，对I、Q两路信号进行AD采样后，进行如下处理：

采用本发明提供的非接触式生命体征检测系统的心跳频率检测算法对实验采集到的信号进行处理：

步骤1，对非接触式生命体征检测系统的连续波多普勒雷达前端输出的I、Q两路信号分别进行带通滤波，带通滤波器的低截止频率为0.5Hz、高截止频率为5.2Hz。

步骤2，通过圆心估计算法恢复经过滤波的I、Q两路信号的有用直流分量。

步骤3，对步骤2得到的两路信号进行反正切解调。图2所示为实验实例中得到的解调信号波形，由图可见，此信号中包含心跳信号以及干扰信号，其中13s之间23s有比较大的人体抖动干扰。

步骤4，采用完全聚类经验模态分解算法对解调信号进行分解，得到一系列本征模函数分量，具体计算步骤如下：

步骤4.1，构造I个信号x[n]+ε₀wⁱ[n]，其中x[n]是反正切解调得到的信号，wⁱ[n],(i＝1,2,...,I)是分布为N(0,1)的白噪声。利用普通经验模态分解算法分别求解I个信号x[n]+ε₀wⁱ[n]的第一个本征模分量并且计算 ${\mathrm{IMF}}_1\[n\]=\frac{1}{I}\underset{i=1}{\overset{I}{\Σ}}{\mathrm{IMF}}_1^i\[n\]={\stackrel{\‾}{IMF}}_1\[n\].$

步骤4.2，计算第一个剩余分量r₁[n]＝x[n]-IMF₁[n]。

步骤4.3，对r₁[n]+ε₁E₁(wⁱ[n]),i＝1,...,I进行经验模态分解(算子E_j(·)的功能是求一个信号经验模态分解的第j个模态)，得到第一个本征模分量，并由此得到 ${\mathrm{IMF}}_2\[n\]=\frac{1}{I}{\mathrm{\Σ}}_{i=1}^I{E}_1(r_1\[n\]+{\mathrm{\ϵ}}_1{E}_1(w^i\[n\]\left)\right)$

步骤4.4，对于k＝2,...,K计算第k个剩余分量：r_k[n]＝r_(k-1)[n]-IMF_k[n]

步骤4.5，对r_k[n]+ε_kE_k(wⁱ[n]),i＝1,...,I进行经验模态分解，得到第一个本征模分量，并由此得到 ${\mathrm{IMF}}_{(k+1)}\[n\]=\frac{1}{I}{\mathrm{\Σ}}_{i=1}^I{E}_1(r_k\[n\]+{\mathrm{\ϵ}}_k{E}_k(w^i\[n\]\left)\right).$

步骤4.6，重复步骤4.4到步骤4.5，步，当剩余分量不再满足可以分解的要求时(剩余分量的极值不足两个)，停止运算，最后的剩余分量表示为： $R\[n\]=x\[n\]-{\mathrm{\Σ}}_{k=1}^K{\mathrm{IMF}}_k\[n\],$ K是总模态数。

本实例中，取I＝800，ε₀＝0.2，ε_k＝0.2/std(E_k(wⁱ[n])),k＝1,...,K。图3所示为对实验实例中解调信号进行完全聚类经验模态分解后得到的一系列本征模函数的波形，包括第1—9个本征模函数。

步骤5，对所有本征模函数分别做希尔伯特变换，并且求得各本征模函数对应的希尔伯特边际谱，具体计算步骤如下：

步骤5.1，令u_k(n)＝IMF_k[n]，首先对u_k(n)做离散希尔伯特变换： $v_k\left(n\right)=\underset{m=0}{\overset{N-1}{\Σ}}h(n-m)u_k\left(n\right),$ 其中N为偶数，且 $h\left(n\right)=\frac{2}{N}{\sin }^2\left(\frac{\mathrm{\π}n}{2}\right)\cot \left(\frac{\mathrm{\π}n}{N}\right).$

步骤5.2，计算 $a_k\left(n\right)=\sqrt{u_k{\left(n\right)}^2+v_k{\left(n\right)}^2}.$

步骤5.3，计算 ${\tilde{f}}_k\left(n\right)=\frac{u_k\left(n\right)\left(v_k\right(n+1)-v_k(n-1))-v_k\left(n\right)\left(u_k\right(n+1)-u_k(n-1))}{2t_s(u_k{\left(n\right)}^2+v_k{\left(n\right)}^2)},$ 对上式进行修正 $f_k\left(n\right)=f_{s}acsin({\tilde{f}}_k/f_s).$

步骤5.4，求希尔伯特谱H_k(f,n)＝H_k(f_k(n),n)＝a_k(n)。

步骤5.5，求得各本征模函数IMF_k[n]对应的希尔伯特边际谱 $h_k\left(f\right)=\underset{n=0}{\overset{N-1}{\Σ}}{H}_k(f,n).$

图4为第1—6个本征模函数对应的希尔伯特边际谱，记为边际谱1—6。由于第六个以后的本征模函数能量极小，而且频率远远低于心跳频率，所以不予考虑。

步骤6，在各本征模函数对应的希尔伯特边际谱上进行峰值检测。本实例中，边际谱1—6的谱峰位置分别为：3.5254Hz、1.2598Hz、1.2402Hz、0.8105Hz、0.4395Hz、0.3418Hz。

步骤7，根据各希尔伯特边际谱的谱峰位置和谱峰附近的能量集中度来判断心跳信号对应的希尔伯特边际谱，具体方法如下：

步骤7.1，判断各边际谱h_k(f)的谱峰对应的频率f_peak,k是否在0.5Hz到3Hz之间。本实例中，由步骤6求得的谱峰位置可知，符合条件的为边际谱2、3、4。

步骤7.2，对于满足步骤7.1条件的边际谱(在本实例中是边际谱2、3、4)，计算其在f_peak,k-0.1Hz到f_peak,k+0.1Hz之间各频率点的幅度的总和A_k，并计算其所有频率点的幅度总和B_k，求能量集中度E_k＝A_k/B_k。本实例中，边际谱2、3、4的能量集中度分别为：E₂＝0.4609、E₃＝0.7806、E₄＝0.3242。

步骤7.3，当E_k大于阈值α时，认为h_k(f)是心跳信号的边际谱；若有多个边际谱满足E_k大于阈值α的条件，则认为峰值位置处谱幅度最大的边际谱是心跳信号的边际谱。本实例中，将阈值α设为0.7，由步骤7.2计算得到的能量集中度可知，只有边际谱3满足条件，则此边际谱就是心跳信号的边际谱。

步骤8，在心跳信号对应的希尔伯特边际谱上提取心跳信号频率信息。本实例中心跳信号的希尔伯特边际谱峰值对应的是频率为1.2402Hz，则最终得到的心跳频率为1.2402Hz，此结果与心电监测仪测得的心率74次/分(1.23333Hz)高度一致。

由实施实例可见，本发明提供的非接触式生命体征监测系统的心跳频率检测算法可在人体不稳定、环境干扰噪声大的条件下有效提取心跳信号并且获得准确心跳频率信息。

Claims (6)

1.一种用于非接触式生命体征检测系统的心跳频率检测算法，包括如下步骤：

步骤1，对非接触式生命体征检测系统的连续波多普勒雷达前端输出的I、Q两路信号分别进行带通滤波；其中，带通滤波器的低截止频率为0.5Hz、高截止频率为5.2Hz；

步骤2，通过圆心估计算法恢复经过滤波的I、Q两路信号的有用直流分量；

步骤3，对步骤2得到的I、Q两路信号进行反正切解调；

其特征在于，得到解调信号之后的步骤如下：

步骤4，采用完全聚类经验模态分解算法对解调信号进行分解，得到一系列本征模函数；

步骤5，对得到的本征模函数分别做希尔伯特变换，并且求得各本征模函数对应的希尔伯特边际谱；

步骤6，在各本征模函数对应的希尔伯特边际谱上进行峰值检测；

步骤7，根据希尔伯特边际谱的谱峰位置和谱峰附近的能量集中度来判断心跳信号对应的希尔伯特边际谱；

步骤8，在心跳信号对应的希尔伯特边际谱上提取心跳信号频率信息。

2.根据权利要求1所述的一种非接触式生命体征检测系统的心跳频率检测算法，其特征在于，步骤4中采用完全聚类经验模态分解算法对解调信号进行分解，具体步骤如下：

步骤4.1，构造I个信号x[n]+ε₀wⁱ[n]，其中x[n]是反正切解调得到的信号，wⁱ[n],(i＝1,2,...,I)是分布为N(0,1)的白噪声，利用普通经验模态分解算法分别求解I个信号x[n]+ε₀wⁱ[n]的第一个本征模分量IMF₁ ⁱ[n]，并且计算：

${\mathrm{IMF}}_1\[n\]=\frac{1}{I}\underset{i=1}{\overset{I}{\Σ}}{\mathrm{IMF}}_1^i\[n\]={\stackrel{\‾}{IMF}}_1\[n\];$

步骤4.2，计算第一个剩余分量r₁[n]＝x[n]-IMF₁[n]；

步骤4.3，对r₁[n]+ε₁E₁(wⁱ[n]),i＝1,...,I进行经验模态分解，得到第一个本征模分量；其中，算子E_j(·)的功能是利用普通经验模态分解算法求一个信号的第j个模态，并由此得到：

${\mathrm{IMF}}_2\[n\]=\frac{1}{I}{\mathrm{\Σ}}_{i=1}^I{E}_1(r_1\[n\]+{\mathrm{\ϵ}}_1{E}_1(w^i\[n\]\left)\right);$

步骤4.4，根据k＝2,...,K，计算第k个剩余分量：r_k[n]＝r_(k-1)[n]-IMF_k[n]；

步骤4.5，根据r_k[n]+ε_kE_k(wⁱ[n]),i＝1,...,I，进行经验模态分解，得到第一个本征模分量，并由此得到 ${\mathrm{IMF}}_{(k+1)}\[n\]=\frac{1}{I}{\mathrm{\Σ}}_{i=1}^I{E}_1(r_k\[n\]+{\mathrm{\ϵ}}_k{E}_k(w^i\[n\]\left)\right);$

步骤4.6，重复步骤4.4到步骤4.5步，当剩余本征模分量不再满足可以分解的要求时，即剩余本征模分量的极值不足两个时，停止运算，最后的剩余本征模分量表示为： $R\[n\]=x\[n\]-{\mathrm{\Σ}}_{k=1}^K{\mathrm{IMF}}_k\[n\],$ K是总模态数；

至此，解调信号x[n]被分解为K个本征模函数IMF_k[n]和一个剩余分量R[n]。

3.根据权利要求1所述的一种非接触式生命体征检测系统的心跳频率检测算法，其特征在于，步骤5中求解各本征模函数IMF_k[n]的希尔伯特边际谱，步骤如下：

步骤5.1，令u_k(n)＝IMF_k[n]，首先对u_k(n)做离散希尔伯特变换：

$v_k\left(n\right)=\underset{m=0}{\overset{N-1}{\Σ}}h(n-m)u_k\left(n\right),$ 其中N为偶数，且 $h\left(n\right)=\frac{2}{N}{\sin }^2\left(\frac{\mathrm{\π}n}{2}\right)\cot \left(\frac{\mathrm{\π}n}{N}\right);$

步骤5.2，计算 $a_k\left(n\right)=\sqrt{u_k{\left(n\right)}^2+v_k{\left(n\right)}^2};$

步骤5.3，计算 ${\tilde{f}}_k\left(n\right)=\frac{u_k\left(n\right)\left(v_k\right(n+1)-v_k(n-1))-v_k\left(n\right)\left(u_k\right(n+1)-u_k(n-1))}{2t_s(u_k{\left(n\right)}^2+v_k{\left(n\right)}^2)};$ 对上式进行修正 $f_k\left(n\right)=f_{s}acsin({\tilde{f}}_k/f_s);$

步骤5.4，求希尔伯特谱H_k(f,n)＝H_k(f_k(n),n)＝a_k(n)；

步骤5.5，求得各本征模函数IMF_k[n]对应的希尔伯特边际谱：

$h_k\left(f\right)=\underset{n=0}{\overset{N-1}{\Σ}}{H}_k(f,n).$

4.根据权利要求1所述的一种非接触式生命体征检测系统的心跳频率检测算法，其特征在于，步骤7中根据谱峰位置和谱峰附近的能量集中度来判断心跳信号对应的希尔伯特边际谱，具体方法如下：

步骤7.1，判断各边际谱h_k(f)的谱峰对应的频率f_peak,k是否在0.5Hz到3Hz之间；

步骤7.2，对于满足步骤7.1条件的边际谱，计算其在f_peak,k-0.1Hz到f_peak,k+0.1Hz之间各频率点的幅度的总和A_k，并计算其所有频率点的幅度总和B_k，求能量集中度E_k＝A_k/B_k；

步骤7.3，当E_k大于阈值α时，认为h_k(f)是心跳信号的边际谱；若有多个边际谱满足E_k大于阈值α的条件，则认为峰值幅度最大的边际谱是心跳信号的边际谱。

5.根据权利要求1所述的一种非接触式生命体征检测系统的心跳频率检测算法，其特征在于，步骤8中将心跳信号的希尔伯特边际谱的谱峰位置对应的频率作为心跳信号的频率。

6.根据权利要求5所述的一种非接触式生命体征检测系统的心跳频率检测算法，其特征在于，所述心跳信号的希尔伯特边际谱峰值对应的是频率为1.2402Hz。