CN102018503A - 生命探测雷达中的呼吸与心跳信号的提取方法及装置 - Google Patents

Abstract

本发明提供了一种生命探测雷达中的呼吸与心跳信号的提取方法及装置，所述方法包括以下步骤：对雷达接收信号进行预处理；对预处理后的信号进行经验模态分解，得到有限个固有模态函数；对分解出的每个固有模态函数进行快速傅立叶变换；分别计算每个固有模态函数对应的总能量、呼吸频带能量和心跳频带能量；根据所述总能量、呼吸频带能量和心跳频带能量之间的比值关系选取对应的固有模态函数分别重构呼吸和心跳信号。采用本发明提供的方法及装置，能更精确的提取呼吸与心跳信号，且具有自适应的特点，能够适应环境的变化。

Description

生命探测雷达中的呼吸与心跳信号的提取方法及装置

【技术领域】

本发明涉及信号处理技术领域，尤其涉及一种生命探测雷达中的呼吸与心跳信号的提取方法及装置。

【背景技术】

穿墙生命探测雷达融合了雷达和生物医学工程技术，可穿透非金属介质(如砖墙、废墟等)，非接触式、远距离检测到人体的呼吸和心脏跳动等信息，广泛应用于震后救援、特殊病人监护以及反恐斗争等领域。

穿墙生命探测雷达的基本原理是：发射机通过振荡器产生连续波信号cos(2πf₀t)，由发射天线向目标辐射，连续波信号穿透非金属障碍物，遇到目标后返回。回波信号经过接收天线进入低噪声放大器，混频解调后再进行预处理，射频干扰和杂波干扰则被过滤掉，得到基带信号，基带信号经过数模转换后送入信号处理单元，信号处理单元提取出呼吸和心跳信号。如果目标是人体，则由人的呼吸和心跳共同引起的胸腔周期性的运动将导致回波信号产生一个相移，相当于对回波进行了调制。其中回波信号的表达式为：

R(t)＝cos[2πf₀t-4πd₀/λ-4πx(t)/λ+Φ(t-2d₀/c)]

其中，f0为载波频率，λ为载波波长，c为光速，φ(t)为发射机的相位噪声，d0为雷达-目标距离；x(t)＝Δ₁cos(ω₁t)+Δ₂cos(ω₂t)，代表胸腔运动，Δ₁Δ₂ω₁ω₂分别为呼吸和心跳信号的幅度与频率。经过混频解调，得到基带信号为：

B(t)＝cos[θ+4π[Δ₁cos(ω₁t)+Δ₂cos(ω₂t)]/λ+ΔΦ(t)]

式中θ是距离d₀产生的恒定相移，ΔΦ(t)是总的剩余相位噪声。基带信号携带呼吸和心跳信息。对B(t)进行贝塞尔展开，可知B(t)中包含m₁ω₁+m₂ω₂(m₁，m₂为任意整数)组成的谐波分量。经过进一步的滤波和信号处理便可以得到代表呼吸与心跳信息的ω₁、ω₂分量。

雷达接收机采集的人体回波复合信号具有非平稳特征，随机性强，复合信号中心跳信号收到了呼吸运动和体表微动回波信号的压制和干扰，使检测系统无法检测出规则的心跳信号，因此，必须才去有效的信号处理方法将呼吸和心跳信号分离，提取出比呼吸信号能力小得多的心跳信号。

目前的呼吸和心跳信号提取方法通常基于数字滤波技术或快速傅立叶变换(FFT)从频域上入手，因而呼吸信号谐波干扰和目标信号非平稳特性的影像未得到很好解决。通常需采用自适应处理的方法，将提取的呼吸信号作为参考输入，与同时携带呼吸、心跳频率成分的信号进行自适应处理，从而实现消除呼吸谐波干扰，检测出心跳波形的目的。但该方法的前提是必须对呼吸信号基频做精确估计，才能达到理想的结果。也有运用小波变换理论对检测信号进行分析，以求从混合信号中提取规律性较强的心跳信号，但如何选择恰当的小波基函数是难点，并且基函数选定后，不能随输入信号的变化而改变，不具备自适应性。

在实际应用中，呼吸与心跳信号究竟是以何种关系形成是我们选择信号处理关键所在。线性非时变系统比较容易分析，具有精确的数学解析表达式，同时还可以设计出各种线性非时变系统，用以实现各种信号处理功能。线性非时变系统满足叠加原理，在分离加法合成的信号时特别有用。对于占据不同频段的两个信号分量，可用一个线性滤波器在时域上实现分离。然而，雷达接收信号受系统非线性影像，会经常遇到不属于相加性组合的信号，若单处的依靠线性滤波来分离和处理这些信号分量，无法达到预期效果。

【发明内容】

基于此，有必要提供一种更精确的生命探测雷达中的呼吸与心跳信号的提取方法。

一种生命探测雷达中的呼吸与心跳信号的提取方法，包括以下步骤：对雷达接收信号进行预处理；对预处理后的信号进行经验模态分解，得到有限个固有模态函数；对分解出的每个固有模态函数进行快速傅立叶变换；分别计算每个固有模态函数对应的总能量、呼吸频带能量和心跳频带能量；根据所述总能量与呼吸频带能量、心跳频带能量之间的比值关系选取对应的固有模态函数分别重构呼吸和心跳信号。

其中，所述对雷霆接收信号进行预处理的步骤为：对雷达接收信号低通滤波，滤除高频噪声。

其中，按照如下公式对分解出的每个回有模态函数进行快速傅立叶变换：

${\mathrm{IMF}}_j\left(k\right)=\underset{m=0}{\overset{M-1}{\mathrm{\Σ}}}{\mathrm{imf}}_j\left(m\right)e^{-j2\mathrm{\πmk}/M},k=0...M-1,j=1...N$

其中，M表示快速傅立叶变换的点数，N表示分解得到的固有模态函数imf的数量。

其中，在频域上按照如下公式计算每个固有模态函数对应的总能量：

$E_{\mathrm{total}}\left(j\right)=\frac{1}{M}\underset{k=0}{\overset{M-1}{\mathrm{\Σ}}}{\mathrm{IMF}}_j\left(k\right){\mathrm{IMF}}_j^{*}\left(k\right),j=1..N$

按照如下公式在频域上计算每个固有模态函数对应的1Hz～3Hz范围内的心跳频带能量：

$E_{\mathrm{heart}}\left(j\right)=\frac{1}{M}\underset{k_1}{\overset{k_2}{\mathrm{\Σ}}}{\mathrm{IMF}}_j\left(k\right){\mathrm{IMF}}_j^{*}\left(k\right),j=1..N$

其中，设采样频率为f_s，快速傅里叶变换的点数为M，则1Hz对应数字频率为：k₁＝floor(M/f_s)，3Hz对应数字频率为：k₂＝floor(3M/f_s)，floor表示朝负无穷方向舍入取整；

按照如下公式在频域上计算每个固有模态函数对应的0.2Hz～0.8Hz范围内的呼吸频带能量：

$E_{\mathrm{respire}}\left(j\right)=\frac{1}{M}\underset{k_3}{\overset{k_4}{\mathrm{\Σ}}}{\mathrm{IMF}}_j\left(k\right){\mathrm{IMF}}_j^{*}\left(k\right),j=1..N$

其中，0.2Hz对应数字频率为：k₃＝floor(0.2M/f_s)，0.8Hz对应数字频率为：k₄＝floor(0.8M/f_s)，floor表示朝负无穷方向舍入取整。

其中，计算所述总能量与呼吸频带能量、心跳频带能量之间的比值关系为：

$n_1=\underset{1\≤j\≤N}{\arg \min }[\frac{E_{\mathrm{heart}}\left(j\right)}{E_{\mathrm{total}}\left(j\right)}>\mathrm{\δ}]$

$n_2=\underset{1\≤j\≤N}{\arg \max }[\frac{E_{\mathrm{heart}}\left(j\right)}{E_{\mathrm{total}}\left(j\right)}>\mathrm{\δ}]$

$n_3=\underset{1\≤j\≤N}{\arg \min }[\frac{E_{\mathrm{respire}}\left(j\right)}{E_{\mathrm{total}}\left(j\right)}>\mathrm{\δ}]$

$n_4=\underset{1\≤j\≤n}{\arg \max }[\frac{E_{\mathrm{resprie}}\left(j\right)}{E_{\mathrm{total}}\left(j\right)}>\mathrm{\δ}]$

其中，δ为阈值，取值为0.5；

则按照如下公式重构心跳信号： $h\left(t\right)=\underset{j=n_1}{\overset{n_2}{\mathrm{\Σ}}}{\mathrm{imf}}_j\left(t\right);$

按照如下公式重构呼吸信号： $r\left(t\right)=\underset{j=n_3}{\overset{n_4}{\mathrm{\Σ}}}{\mathrm{imf}}_j\left(t\right).$

此外，还有必要提供一种更精确的生命探测雷达中的呼吸与心跳信号的提取装置。

一种生命探测雷达中的呼吸与心跳信号的提取装置，包括：预处理模块，用于对雷达接收信号进行预处理；经验模态分解模块，用于对预处理后的信号进行经验模态分解，得到有限个固有模态函数；傅立叶变换模块，用于对分解出的每个固有模态函数进行快速傅立叶变换；计算模块，用于分别计算每个固有模态函数对应的总能量、呼吸频带能量和心跳频带能量；重构模块，根据所述总能量与呼吸频带能量、心跳频带能量之间的比值关系选取对应的固有模态函数分别重构呼吸和心跳信号。

其中，所述预处理模块为低通滤波器，用于对雷达接收信号进行低通滤波，滤除高频噪声。

其中，所述傅立叶变换模块按照如下公式对分解出的每个有模态函数进行快速傅立叶变换：

${\mathrm{IMF}}_j\left(k\right)=\underset{m=0}{\overset{M-1}{\mathrm{\Σ}}}{\mathrm{imf}}_j\left(m\right)e^{-j2\mathrm{\πmk}/M},k=0...M-1,j=1...N$

其中，M表示快速傅立叶变换的点数，N表示分解得到的固有模态函数imf的数量。

其中，所述计算模块包括：

总能量计算模块，用于按照如下公式在频域上计算每个固有模态函数对应的总能量： $E_{\mathrm{total}}\left(j\right)=\frac{1}{M}\underset{k=0}{\overset{M-1}{\mathrm{\Σ}}}{\mathrm{IMF}}_j\left(k\right){\mathrm{IMF}}_j^{*}\left(k\right),j=1..N;$

心跳频带能量计算模块，按照如下公式在频域上计算每个固有模态函数对应的1Hz～3Hz范围内的心跳频带能量：

$E_{\mathrm{heart}}\left(j\right)=\frac{1}{M}\underset{k_1}{\overset{k_2}{\mathrm{\Σ}}}{\mathrm{IMF}}_j\left(k\right){\mathrm{IMF}}_j^{*}\left(k\right),j=1..N$

其中，设采样频率为f_s，快速傅里叶变换的点数为M，则1Hz对应数字频率为：k₁＝floor(M/f_s)，3Hz对应数字频率为：k₂＝floor(3M/f_s)，floor表示朝负无穷方向舍入取整；

呼吸频带能量计算模块，按照如下公式在频域上计算每个固有模态函数对应的0.2Hz～0.8Hz范围内的呼吸频带能量：

$E_{\mathrm{respire}}\left(j\right)=\frac{1}{M}\underset{k_3}{\overset{k_4}{\mathrm{\Σ}}}{\mathrm{IMF}}_j\left(k\right){\mathrm{IMF}}_j^{*}\left(k\right),j=1..N$

其中，0.2Hz对应数字频率为：k₃＝floor(0.2M/f_s)，0.8Hz对应数字频率为：k₄＝floor(0.8M/f_s)，floor表示朝负无穷方向舍入取整。

其中，所述重构模块包括：

比值计算模块，按照如下公式计算所述总能量、呼吸频带能量和心跳频带能量之间的比值关系：

$n_1=\underset{1\≤j\≤N}{\arg \min }[\frac{E_{\mathrm{heart}}\left(j\right)}{E_{\mathrm{total}}\left(j\right)}>\mathrm{\δ}]$

$n_2=\underset{1\≤j\≤N}{\arg \max }[\frac{E_{\mathrm{heart}}\left(j\right)}{E_{\mathrm{total}}\left(j\right)}>\mathrm{\δ}]$

$n_3=\underset{1\≤j\≤N}{\arg \min }[\frac{E_{\mathrm{resprie}}\left(j\right)}{E_{\mathrm{total}}\left(j\right)}>\mathrm{\δ}]$

$n_4=\underset{1\≤j\≤N}{\arg \max }[\frac{E_{\mathrm{resprie}}\left(j\right)}{E_{\mathrm{total}}\left(j\right)}>\mathrm{\δ}]$

其中，δ为阈值，取值为0.5；

心跳信号重构模块，按照如下公式重构心跳信号： $h\left(t\right)=\underset{j=n_1}{\overset{n_2}{\mathrm{\Σ}}}{\mathrm{imf}}_j\left(t\right);$

呼吸信号重构模块，按照如下公式重构呼吸信号： $r\left(t\right)=\underset{j=n_3}{\overset{n_4}{\mathrm{\Σ}}}{\mathrm{imf}}_j\left(t\right).$

上述生命探测雷达中的呼吸与心跳信号的提取方法及装置，通过对预处理后的信号进行经验模态分解，在时域上重构呼吸和心跳信号，避免了呼吸信号谐波对心跳信号检测的干扰，因而能精确的提取呼吸和心跳信号。

此外，经验模块分解根据信号自身特点分解成有限个固有模态函数，不需要选择基函数，分解出来的固有模态函数分量根据输入信号的变换而变换，具有自适应性，能够适应环境的变换。

【附图说明】

图1为一个实施例中呼吸与心跳信号的提取方法的流程图；

图2为生命探测雷达系统的结构示意图；

图3为一个实施例中呼吸与心跳信号的提取装置的结构示意图；

图4为一个实施例中计算模块的结构示意图；

图5为一个实施例中重构模块的结构示意图。

【具体实施方式】

如图1所示，一种生命探测雷达中的呼吸与心跳信号的提取方法，具体过程如下：

步骤S10，对雷达接收信号进行预处理。具体是对雷达接收信号低通滤波，滤除高频噪声，从而提高信噪比。

步骤S20，对预处理后的信号进行经验模态分解，得到有限个固有模态函数。经验模态分解(Empirical Mode Decomposition，EDM)算法可以处理非线性、非平稳信号，与小波变换相比，其不需要选择基函数，而是根据信号自身特点自适应的分解成有限个固有模态函数(Intrinsic Mode Function，imf)，不同的imf反应了信号在不同时间尺度上的演化特征，因此可以通过处理imf分量，从时域上提取呼吸和心跳信号，避免呼吸信号谐波干扰。

对于接收信号x(t)，EDM通过筛选过程将其分解成多个特征时间尺度的固有模态函数，这些函数成为imf。每个imf分量必须满足两个条件：(1)极值点(包括极大值和极小值)的个数与过零点的个数相等或至多相差一个；(2)由局部极大值点构成的上包络和局部极小值点构成的下包络的均值为零。筛选过程如下：

A.确定ε的取值(0.2-0.3)，j←1；

B.c_j-1(t)←x(t)；

C.提取第j个imf分量：

a.i←1，h_j，i-1(t)←c_j-1(t)；

b.提取h_j，i-1(t)的所有极值点；

c.分别对极大值点和极小值点采用三次样条插值得到h_j，i-1(t)的上包络U_j，i-1(t)和下包络L_j，i-1(t)；

d.计算包络的平均值：m_j，i-1(t)←[U_j，i-1(t)+L_j，i-1(t)]/2；

e.更新：h_j，i(t)←h_j，i-1(t)-m_j，i-1(t)，i←i+1；

f.计算停止迭代准则：

$\mathrm{SD}\left(i\right)=\underset{i=0}{\overset{N}{\mathrm{\Σ}}}\frac{{|h_{j,i-1}\left(t\right)-h_{j,i}\left(t\right)|}^2}{{\left|h_{j,i-1}\left(t\right)\right|}^2}$

g.重复步骤b)～f)，当SD(i)＜ε时，imf_j(t)←h_j，i(t)；

D.更新剩余分量：c_j(t)←c_j-1(t)-imf_j(t)；

E.重复步骤C，j←j+1，当c_j(t)的极值点小于2时，结束整个过程。经过EMD筛选，雷达接收信号x(t)被分解成N个固有模态函数imf_j(t)，j＝1，...，N和一个剩余分量c_N(t)，x(t)可以表示为：

EMD将信号自适应分解成不同特征尺度上的分段固有模态函数imf，这些imf从小的时间尺度到大的时间尺度依序排列，体现了EMD的多分辨率特征。低阶imf代表快的振动模式，高阶imf代表慢的振动模块，在频域上，表现为从高频到低频的层层过滤。

步骤S30，对分解出的每个固有模态函数进行快速傅立叶变换。一个实施例中，对每个固有模态函数进行的快速傅立叶变换为：

${\mathrm{IMF}}_j\left(k\right)=\underset{m=0}{\overset{M-1}{\mathrm{\Σ}}}{\mathrm{imf}}_j\left(m\right)e^{-j2\mathrm{\πmk}/M},k=0...M-1,j=1...N$

其中，M表示快速傅立叶变换的点数，N表示分解得到的固有模态函数imf的数量。

步骤S40，分别计算每个固有模态函数对应的总能量、呼吸频带能量和心跳频带能量。在频域上分别计算每个固有模态函数的总能量为：

$E_{\mathrm{total}}\left(j\right)=\frac{1}{M}\underset{k=0}{\overset{M-1}{\mathrm{\Σ}}}{\mathrm{IMF}}_j\left(k\right){\mathrm{IMF}}_j^{*}\left(k\right),j=1..N.$

由于生命信号的重要结构集中在低频部分，心跳频率分布范围通常在1Hz～3Hz，呼吸频率分布范围通常在0.2Hz～0.8Hz，因此仅需要使用反应生命信号频谱结构特征的固有模态函数来重构呼吸与心跳信号。

在一个实施例中，在频域上计算每个固有模态函数对应的1Hz～3Hz范围内的心跳频带能量E_heart(j)，设采样频率为f_s，快速傅里叶变换的点数为M，则1Hz对应数字频率为：k₁＝floor(M/f_s)，3Hz对应数字频率为：k₂＝floor(3M/f_s)，floor表示朝负无穷方向舍入取整。心跳频带能量的计算公式为：

$E_{\mathrm{heart}}\left(j\right)=\frac{1}{M}\underset{k_1}{\overset{k_2}{\mathrm{\Σ}}}{\mathrm{IMF}}_j\left(k\right){\mathrm{IMF}}_j^{*}\left(k\right),j=1..N.$

在一个实施例中，在频域上计算每个固有模态函数对应的0.2Hz～0.8Hz范围内的呼吸频带能量E_respire(j)，设采样频率为f_s，快速傅里叶变换的点数为M，则0.2Hz对应数字频率为：k₃＝floor(0.2M/f_s)，0.8Hz对应数字频率为：k₄＝floor(0.8M/f_s)，floor表示朝负无穷方向舍入取整。呼吸频带能量的计算公式为：

$E_{\mathrm{respire}}\left(j\right)=\frac{1}{M}\underset{k_3}{\overset{k_4}{\mathrm{\Σ}}}{\mathrm{IMF}}_j\left(k\right){\mathrm{IMF}}_j^{*}\left(k\right),j=1..N.$

步骤S50，根据总能量、呼吸频带能量和心跳频带能量之间的比值关系选择对应的固有模态函数分别重构呼吸和心跳信号。

在一个实施方式中，计算所述总能量、呼吸频带能量和心跳频带能量之间的比值关系为：

$n_1=\underset{1\≤j\≤N}{\arg \min }[\frac{E_{\mathrm{heart}}\left(j\right)}{E_{\mathrm{total}}\left(j\right)}>\mathrm{\δ}]$

$n_2=\underset{1\≤j\≤N}{\arg \max }[\frac{E_{\mathrm{heart}}\left(j\right)}{E_{\mathrm{total}}\left(j\right)}>\mathrm{\δ}]$

$n_3=\underset{1\≤j\≤N}{\arg \min }[\frac{E_{\mathrm{resprie}}\left(j\right)}{E_{\mathrm{total}}\left(j\right)}>\mathrm{\δ}]$

$n_4=\underset{1\≤j\≤N}{\arg \max }[\frac{E_{\mathrm{resprie}}\left(j\right)}{E_{\mathrm{total}}\left(j\right)}>\mathrm{\δ}]$

其中，δ为阈值，取值为0.5。上式表示心跳或者呼吸频段范围内的能量占整个频率能量一半以上的固有模态函数分量体现了生命信号的结构特征，这些固有模态分量则可被用来重构呼吸和心跳信号。

在一个实施例中，存在模式索引n₁，n₂，n₃，n₄，心跳信号h(t)则将由(n₂-n₁+1)个固有模特函数重构：

呼吸信号r(t)则将由(n₄-n₃+1)个固有模态函数重构：

如图2所示，生命探测雷达系统包括发射天线1、振荡器2、接收天线3、低噪声放大器4、混频器5、干扰抑制电路6、数模转换器7和信号处理单元8，其中：振荡器2产生连续波信号，通过发射天线1向目标辐射，连续波信号传统非金属障碍物，遇到目标后返回，得到回波信号，该回波信号经过接收天线3进入到低噪声放大器4中，经过混频器5以及干扰抑制电路6，得到基带信号，基带信号经过数模转换器7后送入到信号处理单元8。

信号处理单元8即呼吸与心跳信号的提取装置，如图3所示，呼吸与心跳信号的提取装置包括预处理模块81、经验模态分解模块82、傅立叶变换模块83、计算模块84和重构模块85，其中：预处理模块81用于对雷达接收信号进行预处理；经验模态分解模块82用于对预处理后的信号进行经验模态分解，得到有限个固有模态函数；傅立叶变换模块83用于对分解出的每个固有模态函数进行快速傅立叶变换；计算模块84用于分别计算每个固有模态函数对应的总能量、呼吸频带能量和心跳频带能量；重构模块85用于根据所述总能量、呼吸频带能量和心跳频带能量之间的比值关系选取对应的固有模态函数分别重构呼吸和心跳信号。

在一个实施例中，预处理模块81是低通滤波器，对雷达接收信号低通滤波，滤除高频噪声，从而提高信噪比。经验模态分解模块82通过经验模态分解(Empirical Mode Decomposition，EDM)算法将雷达接收信号x(t)分解成N个固有模态函数imf_j(t)，j＝1，...，N和一个剩余分量c_N(t)，x(t)可以表示为：

$x\left(t\right)=\underset{j=1}{\overset{N}{\mathrm{\Σ}}}{\mathrm{imf}}_j\left(t\right)+c_N\left(t\right).$

在一个实施方式中，傅立叶变换模块83对每个固有模态函数进行的快速傅立叶变换为：

${\mathrm{IMF}}_j\left(k\right)=\underset{m=0}{\overset{M-1}{\mathrm{\Σ}}}{\mathrm{imf}}_j\left(m\right)e^{-j2\mathrm{\πmk}/M},k=0...M-1,j=1...N$

其中，M表示快速傅立叶变换的点数，N表示分解得到的固有模态函数imf的数量。

如图4所示，一个实施例中，计算模块84包括总能量计算模块841、心跳频带能量计算模块842和呼吸频带能量计算模块843，其中：总能量计算模块841按照如下公式在频域上分别计算每个固有模态函数的总能量：

$E_{\mathrm{total}}\left(j\right)=\frac{1}{M}\underset{k=0}{\overset{M-1}{\mathrm{\Σ}}}{\mathrm{IMF}}_j\left(k\right){\mathrm{IMF}}_j^{*}\left(k\right),j=1..N.$

由于生命信号的重要结构集中在低频部分，心跳频率分布范围通常在1Hz～3Hz，呼吸频率分布范围通常在0.2Hz～0.8Hz，因此仅需要使用反应生命信号频谱结构特征的固有模态函数来重构呼吸与心跳信号。

心跳频带能量计算模块842用于在频域上计算每个固有模态函数对应的1Hz～3Hz范围内的心跳频带能量E_heart(j)，设采样频率为f_s，快速傅里叶变换的点数为M，则1Hz对应数字频率为：k₁＝floor(M/f_s)，3Hz对应数字频率为：k₂＝floor(3M/f_s)，floor表示朝负无穷方向舍入取整。心跳频带能量的计算公式为：

$E_{\mathrm{heart}}\left(j\right)=\frac{1}{M}\underset{k_1}{\overset{k_2}{\mathrm{\Σ}}}{\mathrm{IMF}}_j\left(k\right){\mathrm{IMF}}_j^{*}\left(k\right),j=1..N.$

呼吸频带能量计算模块843用于在频域上计算每个固有模态函数对应的0.2Hz～0.8Hz范围内的呼吸频带能量E_respire(j)，设采样频率为f_s，快速傅里叶变换的点数为M，则0.2Hz对应数字频率为：k₃＝floor(0.2M/f_s)，0.8Hz对应数字频率为：k₄＝floor(0.8M/f_s)，floor表示朝负无穷方向舍入取整。呼吸频带能量的计算公式为：

$E_{\mathrm{respire}}\left(j\right)=\frac{1}{M}\underset{k_3}{\overset{k_4}{\mathrm{\Σ}}}{\mathrm{IMF}}_j\left(k\right){\mathrm{IMF}}_j^{*}\left(k\right),j=1..N.$

如图5所示，在一个实施例中，重构模块85包括比值计算模块851、心跳信号重构模块852和呼吸信号重构模块853，其中：比值计算模块851按照如下公式计算所述总能量、呼吸频带能量和心跳频带能量之间的比值关系：

$n_1=\underset{1\≤j\≤N}{\arg \min }[\frac{E_{\mathrm{heart}}\left(j\right)}{E_{\mathrm{total}}\left(j\right)}>\mathrm{\δ}]$

$n_2=\underset{1\≤j\≤N}{\arg \max }[\frac{E_{\mathrm{heart}}\left(j\right)}{E_{\mathrm{total}}\left(j\right)}>\mathrm{\δ}]$

$n_3=\underset{1\≤j\≤N}{\arg \min }[\frac{E_{\mathrm{resprie}}\left(j\right)}{E_{\mathrm{total}}\left(j\right)}>\mathrm{\δ}]$

$n_4=\underset{1\≤j\≤N}{\arg \max }[\frac{E_{\mathrm{resprie}}\left(j\right)}{E_{\mathrm{total}}\left(j\right)}>\mathrm{\δ}]$

其中，δ为阈值，取值为0.5。上式表示心跳或者呼吸频段范围内的能量占整个频率能量一半以上的固有模态函数分量体现了生命信号的结构特征，这些固有模态分量则可被用来重构呼吸和心跳信号。

该实施例中，存在模式索引n₁，n₂，n₃，n₄，心跳信号重构模块852则通过(n₂-n₁+1)个固有模特函数重构心跳信号：

呼吸信号重构模块853则通过(n₄-n₃+1)个固有模态函数重构呼吸信号：

在仿真中，生命探测雷达系统采用单发单收天线，发射载频为1GHz的连续波，雷达-目标距离为2m。人位于墙后0.5m处，且处于静立，均匀呼吸状态，呼吸频率为0.33Hz(20次/分)，心跳频率为1.2Hz(72次/分)。

仿真结果主要对比参数包括：呼吸检测错误率、心跳检测错误率。其中，错误率＝(实际心跳(呼吸)频率-实测心跳(呼吸)频率)/实际心跳(呼吸)频率。仿真中，在不同信噪比的条件下对实施例中所描述的方案与现有方案进行了对比。

仿真结果如表1、表2所示，我们将现有方案与上述呼吸与心跳信号提取方法在不同信噪比的条件下进行了性能对比。从表1中可以看出，对于呼吸频率的估计本发明的方案获得比较好的性能，相对于现有方案，在相同的信噪比条件下本发明的方案呼吸频率估计错误率要小于现有方案。从表2可以看出，在低信噪比的情况下(SNR＜-6dB)，本发明的方案心跳频率估计错误率要大于现有方案，但是随着信噪比的提高，本发明的方案心跳频率估计错误率小于现有方案。对比表1与表2，本发明的方案对呼吸频率的估计效果好于对心跳频率的估计，这种情况可以解释为：呼吸信号幅度远大于心跳信号幅度，呼吸频率的估计变得更加容易。

表1：呼吸频率(实际值：0.33Hz)在不同信噪比条件下的估计性能对比

表2：心跳频率(实际值：1.2Hz)在不同信噪比条件下的估计性能对比

以上所述实施例仅表达了本发明的几种实施方式，其描述较为具体和详细，但并不能因此而理解为对本发明专利范围的限制。应当指出的是，对于本领域的普通技术人员来说，在不脱离本发明构思的前提下，还可以做出若干变形和改进，这些都属于本发明的保护范围。因此，本发明专利的保护范围应以所附权利要求为准。

Claims (10)

1.一种生命探测雷达中的呼吸与心跳信号的提取方法，包括以下步骤：

对雷达接收信号进行预处理；

对预处理后的信号进行经验模态分解，得到有限个固有模态函数；

对分解出的每个固有模态函数进行快速傅立叶变换；

分别计算每个固有模态函数对应的总能量、呼吸频带能量和心跳频带能量；

根据所述总能量与呼吸频带能量、心跳频带能量之间的比值关系选取对应的固有模态函数分别重构呼吸和心跳信号。

2.如权利要求1所述的生命探测雷达中的呼吸与心跳信号的提取方法，其特征在于，所述对雷达接收信号进行预处理的步骤为：对雷达接收信号低通滤波，滤除高频噪声。

3.如权利要求1所述的生命探测雷达中的呼吸与心跳信号的提取方法，其特征在于，按照如下公式对分解出的每个固有模态函数进行快速傅立叶变换：

${\mathrm{IMF}}_j\left(k\right)=\underset{m=0}{\overset{M-1}{\mathrm{\Σ}}}{\mathrm{imf}}_j\left(m\right)e^{-j2\mathrm{\πmk}/M},k=0...M-1,j=1...N$

其中，M表示快速傅立叶变换的点数，N表示分解得到的固有模态函数imf的数量。

4.如权利要求3所述的生命探测雷达中的呼吸与心跳信号的提取方法，其特征在于，在频域上按照如下公式计算每个固有模态函数对应的总能量：

$E_{\mathrm{total}}\left(j\right)=\frac{1}{M}\underset{k=0}{\overset{M-1}{\mathrm{\Σ}}}{\mathrm{IMF}}_j\left(k\right){\mathrm{IMF}}_j^{*}\left(k\right),j=1..N$

按照如下公式在频域上计算每个固有模态函数对应的1Hz～3Hz范围内的心跳频带能量：

$E_{\mathrm{heart}}\left(j\right)=\frac{1}{M}\underset{k_1}{\overset{k_2}{\mathrm{\Σ}}}{\mathrm{IMF}}_j\left(k\right){\mathrm{IMF}}_j^{*}\left(k\right),j=1..N$

其中，设采样频率为f_s，快速傅里叶变换的点数为M，则1Hz对应数字频率为：k₁＝floor(M/f_s)，3Hz对应数字频率为：k₂＝floor(3M/f_s)，floor表示朝负无穷方向舍入取整；

按照如下公式在频域上计算每个固有模态函数对应的0.2Hz～0.8Hz范围内的呼吸频带能量：

$E_{\mathrm{respire}}\left(j\right)=\frac{1}{M}\underset{k_3}{\overset{k_4}{\mathrm{\Σ}}}{\mathrm{IMF}}_j\left(k\right){\mathrm{IMF}}_j^{*}\left(k\right),j=1..N$

其中，0.2Hz对应数字频率为：k₃＝floor(0.2M/f_s)，0.8Hz对应数字频率为：k₄＝floor(0.8M/f_s)，floor表示朝负无穷方向舍入取整。

5.如权利要求4所述的生命探测雷达中的呼吸与心跳信号的提取方法，其特征在于，计算所述总能量与呼吸频带能量、心跳频带能量之间的比值关系为：

$n_1=\underset{1\≤j\≤N}{\arg \min }[\frac{E_{\mathrm{heart}}\left(j\right)}{E_{\mathrm{total}}\left(j\right)}>\mathrm{\δ}]$

$n_2=\underset{1\≤j\≤N}{\arg \max }[\frac{E_{\mathrm{heart}}\left(j\right)}{E_{\mathrm{total}}\left(j\right)}>\mathrm{\δ}]$

$n_3=\underset{1\≤j\≤N}{\arg \min }[\frac{E_{\mathrm{respire}}\left(j\right)}{E_{\mathrm{total}}\left(j\right)}>\mathrm{\δ}]$

$n_4=\underset{1\≤j\≤n}{\arg \max }[\frac{E_{\mathrm{resprie}}\left(j\right)}{E_{\mathrm{total}}\left(j\right)}>\mathrm{\δ}]$

其中，δ为阈值，取值为0.5；

则按照如下公式重构心跳信号： $h\left(t\right)=\underset{j=n_1}{\overset{n_2}{\mathrm{\Σ}}}{\mathrm{imf}}_j\left(t\right);$

按照如下公式重构呼吸信号： $r\left(t\right)=\underset{j=n_3}{\overset{n_4}{\mathrm{\Σ}}}{\mathrm{imf}}_j\left(t\right).$

6.一种生命探测雷达中的呼吸与心跳信号的提取装置，其特征在于，包括：

预处理模块，用于对雷达接收信号进行预处理；

经验模态分解模块，用于对预处理后的信号进行经验模态分解，得到有限个固有模态函数；

傅立叶变换模块，用于对分解出的每个固有模态函数进行快速傅立叶变换；

计算模块，用于分别计算每个固有模态函数对应的总能量、呼吸频带能量和心跳频带能量；

重构模块，根据所述总能量与呼吸频带能量、心跳频带能量之间的比值关系选取对应的固有模态函数分别重构呼吸和心跳信号。

7.如权利要求6所述的生命探测雷达中的呼吸与心跳信号的提取装置，其特征在于，所述预处理模块为低通滤波器，用于对雷达接收信号进行低通滤波，滤除高频噪声。

8.如权利要求6所述的生命探测雷达中的呼吸与心跳信号的提取装置，其特征在于，所述傅立叶变换模块按照如下公式对分解出的每个有模态函数进行快速傅立叶变换：

${\mathrm{IMF}}_j\left(k\right)=\underset{m=0}{\overset{M-1}{\mathrm{\Σ}}}{\mathrm{imf}}_j\left(m\right)e^{-j2\mathrm{\πmk}/M},k=0...M-1,j=1...N$

其中，M表示快速傅立叶变换的点数，N表示分解得到的固有模态函数imf的数量。

9.如权利要求8所述的生命探测雷达中的呼吸与心跳信号的提取方法，其特征在于，所述计算模块包括：

总能量计算模块，用于按照如下公式在频域上计算每个固有模态函数对应的总能量：

心跳频带能量计算模块，按照如下公式在频域上计算每个固有模态函数对应的1Hz～3Hz范围内的心跳频带能量：

$E_{\mathrm{heart}}\left(j\right)=\frac{1}{M}\underset{k_1}{\overset{k_2}{\mathrm{\Σ}}}{\mathrm{IMF}}_j\left(k\right){\mathrm{IMF}}_j^{*}\left(k\right),j=1..N$

其中，设采样频率为f_s，快速傅里叶变换的点数为M，则1Hz对应数字频率为：k₁＝floor(M/f_s)，3Hz对应数字频率为：k₂＝floor(3M/f_s)，floor表示朝负无穷方向舍入取整；

呼吸频带能量计算模块，按照如下公式在频域上计算每个固有模态函数对应的0.2Hz～0.8Hz范围内的呼吸频带能量：

$E_{\mathrm{respire}}\left(j\right)=\frac{1}{M}\underset{k_3}{\overset{k_4}{\mathrm{\Σ}}}{\mathrm{IMF}}_j\left(k\right){\mathrm{IMF}}_j^{*}\left(k\right),j=1..N$

其中，0.2Hz对应数字频率为：k₃＝floor(0.2M/f_s)，0.8Hz对应数字频率为：k₄＝floor(0.8M/f_s)，floor表示朝负无穷方向舍入取整。

10.如权利要求9所述的生命探测雷达中的呼吸与心跳信号的提取装置，其特征在于，所述重构模块包括：

比值计算模块，按照如下公式计算所述总能量与呼吸频带能量、心跳频带能量之间的比值关系：

$n_1=\underset{1\≤j\≤N}{\arg \min }[\frac{E_{\mathrm{heart}}\left(j\right)}{E_{\mathrm{total}}\left(j\right)}>\mathrm{\δ}]$

$n_2=\underset{1\≤j\≤N}{\arg \max }[\frac{E_{\mathrm{heart}}\left(j\right)}{E_{\mathrm{total}}\left(j\right)}>\mathrm{\δ}]$

$n_3=\underset{1\≤j\≤N}{\arg \min }[\frac{E_{\mathrm{resprie}}\left(j\right)}{E_{\mathrm{total}}\left(j\right)}>\mathrm{\δ}]$

$n_4=\underset{1\≤j\≤N}{\arg \max }[\frac{E_{\mathrm{resprie}}\left(j\right)}{E_{\mathrm{total}}\left(j\right)}>\mathrm{\δ}]$

其中，δ为阈值，取值为0.5；

心跳信号重构模块，按照如下公式重构心跳信号： $h\left(t\right)=\underset{j=n_1}{\overset{n_2}{\mathrm{\Σ}}}{\mathrm{imf}}_j\left(t\right);$

呼吸信号重构模块，按照如下公式重构呼吸信号： $r\left(t\right)=\underset{j=n_3}{\overset{n_4}{\mathrm{\Σ}}}{\mathrm{imf}}_j\left(t\right).$

Images