一种基于PE和TVF-EMD的超宽带雷达生命体征信号去噪算法

一种基于pe和tvf-emd的超宽带雷达生命体征信号去噪算法
技术领域
1.本发明涉及超宽带雷达生命探测技术领域，特别涉及一种基于pe和tvf-emd的超宽带雷达生命体征信号去噪算法。


背景技术：

2.超宽带雷达探测技术作为一种新兴的生命体征探测手段，本质上是发射电磁波，因而具备极强的穿透能力，可以穿透非金属障碍物，并在人体表面发生反射，同时还具有距离分辨率高、穿透能力强、抗干扰能力强等优势，且不容易受到天气、温度、光照等外部环境因素的影响，是一种非常理想的生命体征探测方法；但是在实际的探测环境中并不只有被观察对象，还存在其他干扰，需要对雷达原始回波进行去噪优化处理以提高信噪比，因此基于超宽带雷达生命体征信号去噪算法的研究，对提高雷达探测生命体征的救援效率及维护社会安定等方面具有十分重要的意义。


技术实现要素：

3.针对现有技术的不足，本发明提出一种基于pe和tvf-emd的超宽带雷达生命体征信号去噪算法。
4.本发明所采取的技术方案是一种基于pe和tvf-emd的超宽带雷达生命体征信号去噪算法，总体流程如图1所示，包括以下8个步骤。
5.步骤1：建立生命体征信号数学模型，具体步骤如下。
6.步骤1.1：假设d0为人体胸部表面到天线的距离，δd为人体呼吸和心跳产生的胸腔周期性变化，dr是由于呼吸产生的胸腔周期性变化，dh是由于心跳产生的胸腔周期性变化，则雷达天线到人体胸部的瞬时距离表达式为。
7.d(t)＝d0+δd＝d0+dr+dh＝d0+arcos(2πfrt)+ahcos(2πfht)
ꢀꢀꢀꢀ
(1)
8.其中，t为慢时间，ar和ah分别为呼吸运动和心跳运动产生的振幅，fr和fh分别为呼吸运动的频率和心跳运动的频率。
9.步骤1.2：假设在检测场景中，除了人体的呼吸运动和心跳运动以外，其他物体均保持静止状态，则雷达信号的冲激响应为。
[0010][0011]
其中，τ为快时间，avδ(τ-τv(t))是人体目标响应，av表示人体微动回波信号幅度，τv(t)是人体微动回波在快时间方向上的时延，是周围静止目标响应之和，ai表示周围静止物体回波信号幅度，τi是周围静止物体回波在快时间方向上的时延，则τv(t)可以表示为。
[0012]
[0013]
其中，v为电磁波的传播速度，τ0为雷达天线与人体之间的固定时延，τr为人体呼吸运动时延，τh为心跳运动时延。
[0014]
步骤1.3：假设超宽带脉冲雷达发射天线所发出的信号为p(τ)，则接收天线收到的回波信号为。
[0015][0016]
步骤1.4：对雷达回波信号进行离散处理得到生命体征信号的数学模型。
[0017][0018]
其中，δ
t
为快时间采样间隔，τ＝mδ
t m＝0,1,2
…
m-1为对快时间进行离散化，m为快时间采样点数；ts为慢时间采样间隔，t＝nt
s n＝0,1,2
…
n-1为对慢时间进行离散化，n为慢时间采样点数。
[0019]
步骤2：对原始雷达回波进行去噪预处理以提高信噪比，具体步骤如下。
[0020]
步骤2.1：道信号相减法，将每一行雷达回波信号与前一行雷达回波信号相减，则该过程中保持恒定不变的量就会被消除；对于雷达回波矩阵s(m,n)，m为雷达回波延迟时间的采样点数，n为扫描时间的采样点数；假设快时间序列为yj，j＝1,2,3
…
,n，其中n为采样点数，则道信号相减法的结果yj'为。
[0021]
yj'＝y
j+1-yjꢀꢀꢀꢀ
(6)
[0022]
步骤2.2：减平均法，对雷达回波数据沿扫描时间方向取平均值，即直流分量；直流分量由静止物体所产生，由于静止物体与探测装置间的距离固定不变，而人体的呼吸信号与心跳信号是呈周期性变化的，因此可以用原始数据减去这一平均值就可以达到去除直流分量的效果，则消除直流分量的结果为。
[0023][0024]
其中，可近似看为背景杂波中的直流分量。
[0025]
步骤2.3：线性趋势抑制，根据扫描时间方向上的线性最小二乘拟合来估计所得到的背景杂波中具有潜在性的线性趋势，随后从雷达回波数据中减掉该线性趋势，就可以达到去除背景杂波的目地，去除结果的计算公式如下。
[0026]
w＝ω
t-x(x
t
x)-1
x
t
ω
t
ꢀꢀꢀ
(8)
[0027]
其中，x＝[x1,x2]，x1＝[1,2,
…
,n]
t
，x2＝[1,1,
…
,1]
t
。
[0028]
步骤2.4：自动控制增益，首先选择2λ+1的时间窗，并根据其中的能量计算相应的增益系数，来完成自适应控制的目的；设第n1帧回波信号为r(τ,n1)，则。
[0029][0030]
re(t,τ)＝gm(t,τ)
×
r(t,τ)
ꢀꢀꢀꢀꢀ
(10)
[0031]
其中gm(t,τi)为增益系数，re(t,τ)为处理后的人体生命体征信号。
[0032]
步骤2.5：快时间滤波，根据小波变换的小波阈值去噪是一种常见的信号去噪方式，如何选择合适的阈值是雷达回波信号去噪的关键；采用广义交叉验证函数来确定阈值，
该方法与传统的d-j阈值相比无需依赖输入输出数据，具有更大的优势，gcv公式表达式如下。
[0033][0034]
其中n表示小波系数的数，n0表示小波系数为0的数量，w表示含噪声的回波信号，w
t
表示降噪后的回波信号。
[0035]
步骤2.6：慢时间滤波，采用巴特沃斯低通滤波器滤去高频的噪声信号，它的平方幅度函数定义为。
[0036][0037]
其中n为滤波器的阶数，ωc为3db截至频率，ε是控制带通波纹幅度的参数，由于人体呼吸频率范围是0.1-0.5hz，心跳频率范围是0.8-3.0hz，所以其频带选择为0.1-3.0hz。
[0038]
步骤3：雷达回波矩阵慢时间序列中的目标区域范围和非目标区域范围的信号幅值不同，采用排列熵(permutation entropy，pe)算法来检测人体目标的位置，具体步骤如下。
[0039]
假设快时间为m1，对长度为n的慢时间序列x[m1,i]＝{y(1),y(2),
…
,y(i)}，i＝1,2,
…
,n，y∈wm×n，进行相空间重构可得。
[0040][0041]
其中τd为延迟时间，q为嵌入维度，k＝n-(q-1)τd；该矩阵的每一行都是重构分量，共有k个重构分量，将第j个重构分量中的元素按升序排列，可以得到新的重构分量。
[0042]
x(i+(j
1-1)τd)≤x(i+(j
2-1)τd)≤
…
≤x(i+(j
q-1)τd)
ꢀꢀꢀꢀ
(14)
[0043]
若其中存在相等的分量，则按ji中下标i的大小进行排列即j
1 j2的排列顺序，对于任意时间序列经过空间重构后得到的重构矩阵中的任意分量可以得到一组序列。
[0044]
s(l)＝(j1,j2,
…
,jq)
ꢀꢀꢀ
(15)
[0045]
其中l＝1,2,
…
,k，k≤q！，m维相空间映射可以产生m！种s(l)排序。
[0046]
计算每种序列s(l)不同排列顺序的出现概率pk，概率和为1，则按照shannon熵的形式，时间序列x[m1,i]的k种不同符号序列的排列熵可以定义为。
[0047][0048]
当即每种排列顺序的概率存在且相等时排列熵的值最大为ln(m！)，因此对排列熵h
p
进行归一化处理可得。
[0049][0050]hp
的值表示慢时间序列x[m1,i]的复杂程度，h
p
的值越小，时间序列变化越规则，反之越随机；在雷达回波矩阵中，目标位置范围在慢时间方向上的数据变化比较规则对应的人体目标的pe值较小，可以通过寻找pe值结果的最低点来判断人体目标的位置p
pos
，公式表示如下。
[0051]
range＝(v
×
p
pos
×
tf)/2
ꢀꢀꢀ
(18)
[0052]
其中v＝3
×
108m|s，tf为快时间方向上的采样间隔。
[0053]
假设人体胸腔的横向距离为d，计算得到胸腔距离在信号中所占用的点数为。
[0054][0055]
进而可以得到生命体征信号矩阵δ的表达式为。
[0056][0057]
步骤4：范围估计，对人体位置处的信号进行选取以及信号的再结合。
[0058]
步骤5：对结合过的信号进行带通滤波的处理。
[0059]
步骤6：对经过预处理后的信号进行tvf-emd分解，具体操作步骤如下。
[0060]
步骤6.1：局部截止频率
[0061]
b-样条拟合滤波器的截止频率是随着时间变化而发生改变的，它主要通过构造多项式样条的方式来逼近输入信号，可表示如下形式。
[0062][0063]
其中βn(t)为b-样条函数，c(k)为b-样条相关系数，为逼近误差，用b-样条阶次n和节点m来确定b-样条相关系数时可使逼近误差最小，公式表示为。
[0064][0065]
其中[]
↑m为上采样操作，即在每两个采样点之间插入m个采样点。
[0066]
根据公式(21)，c(k)可表示为。
[0067][0068]
其中[]
↓m为下采样操作，即每隔m个点进行采样；为预滤波器，公式表示如下。
[0069][0070]
综上，b-样条拟合器可表示为如下形式。
[0071][0072]
根据上述公式可将b-样条逼近看作一种特殊形式的低通滤波器，通常情况下，用
节点m来确定b-样条滤波器的局部截止频率，信号执行时变滤波的时间则由截止频率来确定，进行时变滤波经验模态分解的具体实现步骤如下。
[0073]
步骤6.1.1：利用hilbert变换计算输入信号x(t)的瞬时幅值a(t)与瞬时频率
[0074][0075][0076]
其中x(t)为输入信号x(t)的hilbert变换。
[0077]
步骤6.1.2：确定瞬时幅值a(t)的局部极大值序列a({t
max
})与局部极小值序列a({t
min
})。
[0078]
z(t)作为多分量信号的解析信号可表示为。
[0079][0080]
由此可得。
[0081][0082][0083]
其中表示第i个分量的瞬时相位，ai(t)表示第i个分量的幅值。
[0084]
假设在t
min
处取得局部最小值可得。
[0085][0086]
将公式(31)代入(29)和公式(30)可得。
[0087]
a(t
min
)＝|a1(t
min
)-a2(t
min
)
ꢀꢀꢀ
(32)
[0088][0089]
根据导数运算可得。
[0090]a′1(t
min
)-a'2(t
min
)＝0
ꢀꢀꢀ
(34)
[0091]
同理，可得到如下公式。
[0092][0093]
a(t
max
)＝|a1(t
max
)+a2(t
max
)|
ꢀꢀꢀ
(36)
[0094][0095]a′1(t
max
)+a'2(t
max
)＝0
ꢀꢀꢀ
(38)
[0096]
步骤6.1.3：分别对a({t
max
})与a({t
min
})进行插值拟合可得局部极值拟合函数β1(t)与β2(t)并计算瞬时均值函数a1(t)与瞬时包络函数a2(t)。
[0097][0098][0099]
步骤6.1.4：计算瞬时频率分量和
[0100][0101]
其中η1(t)与η2(t)分别为对进行插值后的结果。
[0102]
步骤6.1.5：计算局部截止频率
[0103][0104]
步骤6.1.6：重新对局部截止频率进行调整来解决间歇问题。
[0105]
局部截至频率可能会受到噪声的影响，因此对进行重构可得到一个新的信号。
[0106][0107]
步骤6.2：局部均值函数
[0108]
利用时变滤波器对输入信号进行滤波来获取局部均值，并将f(t)的极值点作为构造时变滤波器的节点，可以使滤波器的截止频率与局部截止频率保持一致，对输入信号x(t)进行b-样条插值逼近可以得到逼近结果m(t)。
[0109]
步骤6.3：停止准则
[0110]
判断准则如下。
[0111][0112][0113][0114]
其中b
l
(t)为瞬时带宽，为加权均值瞬时频率，ξ为给定的带宽阈值，若满足公式(45)则该信号为局部窄带信号；根据人体的呼吸频率为0.1hz-0.5hz，心跳频率为
0.8hz-3.0hz可将回波信号分为n个imf分量并在频域上对每个分量进行呼吸和心跳能量百分比的计算，计算公式如下。
[0115][0116][0117]
其中q(j)为第j个imf的频域能量，qr(j)为第j个imf呼吸频域能量，qh(j)为第j个imf心跳频域能量，δr与δh分别为呼吸和心跳的判断阈值。
[0118]
步骤7：将符合要求的imf呼吸分量与心跳分量进行重构，公式如下。
[0119]
sr(t)＝∑imf(t)
ꢀꢀꢀꢀ
(50)
[0120]
sh(t)＝σimf(t)
ꢀꢀꢀꢀ
(51)
[0121]
步骤8：利用快速傅里叶变换(fft)对重构后的信号进行频谱分析。
[0122]
采用上述技术方案所产生的有益效果在于：本发明提供一种基于pe和tvf-emd的超宽带雷达生命体征信号去噪算法。该方法通过计算沿慢时间方向的雷达接受脉冲的pe值估算人类目标和雷达之间的距离。采用tvf-emd算法将组合后的雷达回波信号自适应地分解为imfs。根据每个imf在呼吸和心跳频带内的能量百分比，重构呼吸和心跳信号并对其进行fft，从而获得呼吸和心跳的频率。
附图说明
[0123]
图1为本发明一种基于pe和tvf-emd的超宽带雷达生命体征信号去噪算法流程图。
[0124]
图2为本发明生命体征信号模型的图示。
[0125]
图3为本发明预处理算法性能图。
[0126]
其中，(a)原始雷达回波矩阵图；(b)rps处理结果图；(c)ms处理结果图；(d)agc处理结果图；(e)线性趋势抑制结果对比图；(f)快时间滤波结果对比图；(g)慢时间滤波结果对比图。
[0127]
图4为本发明不同位置pe值显示人体位置示意图。
[0128]
其中，(a)0.5m处pe值距离曲线；(b)1.0m处pe值距离曲线；(c)1.5m处pe值距离曲线。
[0129]
图5为本发明其他对比算法的分解图。
[0130]
其中，(a)emd算法分解图；(b)eemd算法分解图；(c)ceemd算法分解图。
[0131]
图6为本发明生命体征信号的tvf-emd分解图。
[0132]
图7为本发明各imf分量的呼吸心跳能量占比图。
[0133]
图8为本发明重构的呼吸信号与心跳信号。
[0134]
图9为本发明呼吸信号与心跳信号的频谱图。
[0135]
其中，(a)呼吸信号频谱图；(b)心跳信号频谱图。
具体实施方式
[0136]
下面结合附图对本发明的具体实施方式作进一步详细描述。
[0137]
在雷达实际的探测过程中，雷达回波信号主要包括目标人体的生命体征信号及周围环境中的噪声干扰和一些其他的杂波信号，生命探测的主要目的是根据人体呼吸、心跳等生命信号的特点在复杂的回波信号中提取出人体的生命体征信息，并对此过程中产生的噪声及杂波干扰采用合适的去噪算法进行滤波降噪的处理，以提高信噪比，整体的算法流程如图1所示；首先在预处理阶段通过道信号相减法对背景杂波进行简单消除；利用直接减平均法消除背景杂波中的直流分量；利用线性趋势抑制的方法消除雷达回波数据中的线性趋势；再通过自动控制增益降低多径效应，从慢时间方向上增强微弱的呼吸信号与心跳信号来提升人体生命体征的信噪比；利用广义交叉验证来确定阈值同真实数据及噪声能量之间的无关性；最后用巴特沃斯低通滤波器滤除回波信号中的高频噪声；然后通过对慢时间方向上数据的峰度、标准差、方差等统计特征进行分析，采用pe算法来检测人体目标的位置；并采用tvf-emd算法对提取后的生命体征信号进行分解，计算分解后各imf分量在呼吸频带与心跳频带内各自的能量占比，选择合适的信号进行重构，最后利用fft变换对重构后的信号进行频谱分析。
[0138]
本实施方式的生命体征信号模型如图2所示，具体步骤如下。
[0139]
假设d0为人体胸部表面到天线的距离，δd为人体呼吸和心跳产生的胸腔周期性变化(dr是由于呼吸产生的胸腔周期性变化，dh是由于心跳产生的胸腔周期性变化)，则雷达天线到人体胸部的瞬时距离表达式为。
[0140]
d(t)＝d0+δd＝d0+dr+dh＝d0+arcos(2πfrt)+ahcos(2πfht)
ꢀꢀꢀꢀ
(1)
[0141]
其中，t为慢时间，ar和ah分别为呼吸运动和心跳运动产生的振幅，fr和fh分别为呼吸运动的频率和心跳运动的频率。
[0142]
假设在检测场景中，除了人体的呼吸运动和心跳运动以外，其他物体均保持静止状态，则雷达信号的冲激响应为。
[0143][0144]
其中，τ为快时间，avδ(τ-τv(t))是人体目标响应，av表示人体微动回波信号幅度，τv(t)是人体微动回波在快时间方向上的时延，是周围静止目标响应之和，ai表示周围静止物体回波信号幅度，τi是周围静止物体回波在快时间方向上的时延，则τv(t)可以表示为。
[0145][0146]
其中，v为电磁波的传播速度，τ0为雷达天线与人体之间的固定时延，τr为人体呼吸运动时延，τh为心跳运动时延。
[0147]
假设超宽带脉冲雷达发射天线所发出的信号为p(τ)，则接收天线收到的回波信号为。
[0148][0149]
对雷达回波信号进行离散处理得到生命体征信号的数学模型。
[0150]
[0151]
其中，δ
t
为快时间采样间隔，τ＝mδ
t m＝0,1,2
…
m-1为对快时间进行离散化，m为快时间采样点数；ts为慢时间采样间隔，t＝nt
s n＝0,1,2
…
n-1为对慢时间进行离散化，n为慢时间采样点数。
[0152]
本实施方式的预处理算法性能图如图3所示，具体步骤如下。
[0153]
道信号相减法，对雷达回波信号中的背景杂波进行消除的最简单又直接的方法就是用每一行雷达回波信号与前一行雷达回波信号相减，在这一过程中保持恒定不变的量就会被消除，对于雷达回波矩阵s(m,n)，m为雷达回波延迟时间的采样点数，n为扫描时间的采样点数，假设快时间序列为yj，j＝1,2,3
…
,n，其中n为采样点数，则道信号相减法的结果yj'为。
[0154]
yj'＝y
j+1-yjꢀꢀꢀꢀ
(6)
[0155]
减平均法，对雷达回波数据沿扫描时间方向取平均值，即直流分量；直流分量由静止物体所产生，由于静止物体与探测装置间的距离固定不变，而人体的呼吸信号与心跳信号是呈周期性变化的，因此可以用原始数据减去这一平均值就可以达到去除直流分量的效果，则消除直流分量的结果为。
[0156][0157]
其中，可近似看为背景杂波中的直流分量。
[0158]
线性趋势抑制，对于由平稳信号和线性趋势混合在一起的背景杂波可采用线性趋势抑制的方法，具体步骤如下：可以根据扫描时间方向上的线性最小二乘拟合来估计所得到的背景杂波中具有潜在性的线性趋势，随后从雷达回波数据中减掉该线性趋势，就可以达到去除背景杂波的目地，去除结果的计算公式如下。
[0159]
w＝ω
t-x(x
t
x)-1
x
t
ω
t
ꢀꢀꢀ
(8)
[0160]
其中，x＝[x1,x2]，x1＝[1,2,
…
,n]
t
，x2＝[1,1,
…
,1]
t
。
[0161]
自动控制增益，该方法主要起到提高生命体征信号能量值的作用，它会根据回波信号矩阵中的能量值的反馈结果来进行自动调节，从而加强生命体征信号，忽略其他无用的信号来达成提高信噪比的目的，具体步骤为：首先选择2λ+1的时间窗，并根据其中的能量计算相应的增益系数，来完成自适应控制的目的；设第n1帧回波信号为r(τ,n1)，则。
[0162][0163]
re(t,τ)＝gm(t,τ)
×
r(t,τ)
ꢀꢀꢀꢀꢀ
(10)
[0164]
其中gm(t,τi)为增益系数，re(t,τ)为处理后的人体生命体征信号。
[0165]
快时间滤波，根据小波变换的小波阈值去噪是一种常见的信号去噪方式，如何选择合适的阈值是雷达回波信号去噪的关键；采用广义交叉验证函数来确定阈值，该方法与传统的d-j阈值相比无需依赖输入输出数据，具有更大的优势，gcv公式表达式如下。
[0166][0167]
其中n表示小波系数的数量，n0表示小波系数为0的数量，w表示含噪声的回波信
号，w
t
表示降噪后的回波信号。
[0168]
慢时间滤波，采用巴特沃斯低通滤波器滤去高频的噪声信号，它的平方幅度函数定义为。
[0169][0170]
其中n为滤波器的阶数，ωc为3db截至频率，ε是控制带通波纹幅度的参数，由于人体呼吸频率范围是0.1-0.5hz，心跳频率范围是0.8-3.0hz，所以其频带选择为0.1-3.0hz。
[0171]
本实施方式的不同位置pe值显示人体位置示意图如图4所示，雷达回波矩阵慢时间序列中的目标区域范围和非目标区域范围的信号幅值不同，采用排列熵算法来检测人体目标的位置，排列熵是一种可以用来表示信号的复杂性以及衡量信息的不确定性的指标，通过比较相邻时间序列的值来检测时间序列的动态变化，时间序列越规律，其对应的排列熵会越小；反之时间序列越复杂则其对应的排列熵越大，由于人体的生命体征信号较探测环境中的噪声杂波相对简单，因此人体目标对应的pe值会较小，这一特性可以作为检测人体目标位置的主要依据，具体步骤如下。
[0172]
假设快时间为m1，对长度为n的慢时间序列x[m1,i]＝{y(1),y(2),
…
,y(i)}，i＝1,2,
…
,n，y∈wm×n，进行相空间重构可得。
[0173][0174]
其中τd为延迟时间，q为嵌入维度，k＝n-(q-1)τd，该矩阵的每一行都是重构分量，共有k个重构分量，将第j个重构分量中的元素按升序排列，可以得到新的重构分量。
[0175]
x(i+(j
1-1)τd)≤x(i+(j
2-1)τd)≤
…
≤x(i+(j
q-1)τd)
ꢀꢀꢀ
(14)
[0176]
若其中存在相等的分量，则按ji中下标i的大小进行排列即j
1 j2的排列顺序，对于任意时间序列经过空间重构后得到的重构矩阵中的任意分量可以得到一组序列。
[0177]
s(l)＝(j1,j2,
…
,jq)
ꢀꢀꢀ
(15)
[0178]
其中l＝1,2,
…
,k，k≤q！，m维相空间映射可以产生m！种s(l)排序。
[0179]
计算每种序列s(l)不同排列顺序的出现概率pk，概率和为1，则按照shannon熵的形式，时间序列x[m1,i]的k种不同符号序列的排列熵可以定义为。
[0180][0181]
当即每种排列顺序的概率存在且相等时排列熵的值最大为ln(m！)，因此对排列熵h
p
进行归一化处理可得。
[0182][0183]hp
的值表示慢时间序列x[m1,i]的复杂程度，h
p
的值越小，时间序列变化越规则，反之越随机；在雷达回波矩阵中，目标位置范围在慢时间方向上的数据变化比较规则对应的人体目标的pe值较小，可以通过寻找pe值结果的最低点来判断人体目标的位置p
pos
，公式表示如下。
[0184]
range＝(v
×
p
pos
×
tf)/2
ꢀꢀꢀ
(18)
[0185]
其中v＝3
×
108m/s，tf为快时间方向上的采样间隔。
[0186]
假设人体胸腔的横向距离为d，计算得到胸腔距离在信号中所占用的点数为。
[0187][0188]
进而可以得到生命体征信号矩阵δ的表达式为。
[0189][0190]
本实施方式的其他对比算法的分解图如图5所示，具体步骤如下。
[0191]
采用emd、eemd、ceemd算法对公式h(t)＝sin(20πf1t)+0.5sin(2πf2t)的模拟混合信号进行分解，其中f1＝10，f2＝100，采样频率为500hz，进而比较说明各算法的优缺点。
[0192]
经验模态分解(empirical mode decomposition，emd)是一种对频率随时间变化的非平稳信号进行局部频谱分析的常用方法，主要利用信号的极值点信息将其分解为若干个本征模态函数(intrinsic mode function，imf)，由图5(a)的emd算法分解图可知，该算法存在模态混叠的问题，主要有一个单独的imf信号中含有不同的时间尺度或者相同的时间尺度出现在不同的imf中两种现象道信号相减法。
[0193]
为避免模态混叠问题，提出了一种基于集总经验模式分解(ensemble empirical mode decomposition，eemd)的概念，即在原始信号中添加正态分布的白噪声，并将其看作一个整体，然后再对其进行emd分解，由图5(b)的eemd算法分解图可知，它的不足之处在于加入的白噪声虽在集总平均之后基本抵消但仍存在残留，该残留噪声不可忽略，并且集总平均次数一般在几百次以上十分耗时。
[0194]
针对以上问题又提出了互补集合经验模态分解(complete ensemble empirical mode decomposition，ceemd)，该算法加入一组白噪声后又减去了一组白噪声，这样以来噪声会随着集总平均次数的增加而减少，但同时也存在迭代次数过多的问题，ceemd算法的分解图如图5(c)所示。
[0195]
本实施方式的生命体征信号的tvf-emd分解图如图6所示，该信号被自适应地分解成一系列的本征模态函数，按照低阶到高阶、高频到低频的顺序依次排列开来，每一个本征模态函数都表示出信号在不同的频率范围内的振荡特性，并且阶数会随着imf频率的升高而降低，imf频率越高阶数越低，反之频率越低阶数越高，具体步骤如下。
[0196]
局部截止频率：b-样条拟合滤波器的截止频率是随着时间变化而发生改变的，它主要通过构造多项式样条的方式来逼近输入信号，可表示如下形式。
[0197][0198]
其中βn(t)为b-样条函数，c(k)为b-样条相关系数，为逼近误差，用b-样条阶次n和节点m来确定b-样条相关系数时可使逼近误差最小，公式表示为。
[0199][0200]
其中[]
↑m为上采样操作，即在每两个采样点之间插入m个采样点。
[0201]
根据公式(21)，c(k)可表示为。
[0202][0203]
其中[]
↓m为下采样操作，即每隔m个点进行采样；为预滤波器，公式表示如下。
[0204][0205]
综上，b-样条拟合器可表示为如下形式。
[0206][0207]
根据上述公式可将b-样条逼近看作一种特殊形式的低通滤波器，通常情况下，用节点m来确定b-样条滤波器的局部截止频率，信号执行时变滤波的时间则由截止频率来确定，进行时变滤波经验模态分解的具体实现步骤如下。
[0208]
1.利用hilbert变换计算输入信号x(t)的瞬时幅值a(t)与瞬时频率1.利用hilbert变换计算输入信号x(t)的瞬时幅值a(t)与瞬时频率
[0209][0210]
其中x(t)为输入信号x(t)的hilbert变换。
[0211]
2.确定瞬时幅值a(t)的局部极大值序列a({t
max
})与局部极小值序列a({t
min
})。
[0212]
z(t)作为多分量信号的解析信号可表示为。
[0213][0214]
由此可得。
[0215][0216][0217]
其中表示第i个分量的瞬时相位，ai(t)表示第i个分量的幅值。
[0218]
假设在t
min
处取得局部最小值可得。
[0219][0220]
将公式(31)代入(29)和公式(30)可得。
[0221]
a(t
min
)＝|a1(t
min
)-a2(t
min
)|
ꢀꢀ
(32)
[0222][0223]
根据导数运算可得。
[0224]a′1(t
min
)-a'2(t
min
)＝0
ꢀꢀꢀ
(34)
[0225]
同理，可得到如下公式。
[0226][0227]
a(t
max
)＝|a1(t
max
)+a2(t
max
)|
ꢀꢀꢀ
(36)
[0228][0229]
3.分别对a({t
max
})与a({t
min
})进行插值拟合可得局部极值拟合函数β1(t)与β2(t)并计算瞬时均值函数a1(t)与瞬时包络函数a2(t)。
[0230][0231][0232]
4.计算瞬时频率分量和
[0233][0234][0235]
其中η1(t)与η2(t)分别为对进行插值后的结果。
[0236]
5.计算局部截止频率
[0237][0238]
6.重新对局部截止频率进行调整来解决间歇问题。
[0239]
局部截至频率可能会受到噪声的影响，因此对进行重构可得到一个新的信号。
[0240]
[0241]
局部均值函数：利用时变滤波器对输入信号进行滤波来获取局部均值，并将f(t)的极值点作为构造时变滤波器的节点，可以使滤波器的截止频率与局部截止频率保持一致，对输入信号x(t)进行b-样条插值逼近可以得到逼近结果m(t)。
[0242]
停止准则如下，其中b
l
(t)为瞬时带宽，为加权均值瞬时频率，ξ为给定的带宽阈值。
[0243][0244][0245][0246]
本实施方式的各imf分量的呼吸心跳能量占比如图7所示，横轴表示imf的序列号、纵轴表示各imf分量的能量百分比，对于心跳信号序列imf5、imf6的能量占比均达到90％以上，其余imf在心跳频带的能量占比较小，因此可选用以上两个本征模态函数进行心跳信号的重构；对于呼吸信号序列imf8、imf9、imf10、imf11、imf12的能量占比均达到了90％，其余imf在呼吸频带的能量占比较小，因此可选用以上五个本征模态函数对呼吸信号进行重构操作，具体步骤如下。
[0247]
将回波信号分为n个imf分量并在频域上对每个分量进行呼吸和心跳能量百分比的计算，计算公式如下。
[0248][0249][0250]
其中q(j)为第j个imf的频域能量，qr(j)为第j个imf呼吸频域能量，qh(j)为第j个imf心跳频域能量，δr与δh分别为呼吸和心跳的判断阈值。
[0251]
本实施方式的重构的呼吸信号与心跳信号如图8所示，具体步骤如下。
[0252]
将符合要求的imf呼吸分量与心跳分量进行重构，公式如下。
[0253]
sr(t)＝∑imf(t)
ꢀꢀꢀꢀ
(50)
[0254]
sh(t)＝∑imf(t)
ꢀꢀꢀꢀ
(51)
[0255]
重构呼吸和心跳信号之后，分别使用fft方法进行频域分析，本实施方式的呼吸信号与心跳信号的频谱如图9所示。
[0256]
虽然以上描述了本发明的具体实施方式，但是本领域内的熟练的技术人员应当理解，这些仅是举例说明，可以对这些实施方式做出多种变更或修改，而不背离本发明的原理和实质，本发明的范围仅由所附权利要求书限定。

技术特征：
1.一种基于pe和tvf-emd的超宽带雷达生命体征信号去噪算法，其特征在于，该方法包括以下8个步骤：步骤1：建立生命体征信号数学模型；步骤2：对原始雷达回波进行去噪预处理以提高信噪比；步骤3：雷达回波矩阵慢时间序列中的目标区域范围和非目标区域范围的信号幅值不同，采用pe算法来检测人体目标的位置；步骤4：范围估计，对人体位置处的信号进行选取以及信号的再结合；步骤5：对结合过的信号进行带通滤波的处理；步骤6：对经过预处理后的信号进行tvf-emd分解；步骤7：将符合要求的imf呼吸分量与心跳分量进行重构；步骤8：利用fft变换对重构后的信号进行频谱分析。2.根据权利要求1所述的一种基于pe和tvf-emd的超宽带雷达生命体征信号去噪算法，其特征在于，所述步骤1的过程如下：步骤1.1：假设d0为人体胸部表面到天线的距离，δd为人体呼吸和心跳产生的胸腔周期性变化(d
r
是由于呼吸产生的胸腔周期性变化，d
h
是由于心跳产生的胸腔周期性变化)，则雷达天线到人体胸部的瞬时距离表达式为：d(t)＝d0+δd＝d0+d
r
+d
h
＝d0+a
r
cos(2πf
r
t)+a
h
cos(2πf
h
t)
ꢀꢀꢀꢀ
(1)其中，t为慢时间，a
r
和a
h
分别为呼吸运动和心跳运动产生的振幅，f
r
和f
h
分别为呼吸运动的频率和心跳运动的频率；步骤1.2：假设在检测场景中，除了人体的呼吸运动和心跳运动以外，其他物体均保持静止状态，则雷达信号的冲激响应为：其中，τ为快时间，a
v
δ(τ-τ
v
(t))是人体目标响应，a
v
表示人体微动回波信号幅度，τ
v
(t)是人体微动回波在快时间方向上的时延，是周围静止目标响应之和，a
i
表示周围静止物体回波信号幅度，τ
i
是周围静止物体回波在快时间方向上的时延，则τ
v
(t)可以表示为：其中，v为电磁波的传播速度，τ0为雷达天线与人体之间的固定时延，τ
r
为人体呼吸运动时延，τ
h
为心跳运动时延；步骤1.3：假设超宽带脉冲雷达发射天线所发出的信号为p(τ)，则接收天线收到的回波信号为：步骤1.4：对雷达回波信号进行离散处理得到生命体征信号的数学模型：
其中，δ
t
为快时间采样间隔，τ＝mδ
t
m＝0,1,2
…
m-1为对快时间进行离散化，m为快时间采样点数；t
s
为慢时间采样间隔，t＝nt
s
n＝0,1,2
…
n-1为对慢时间进行离散化，n为慢时间采样点数。3.根据权利要求1所述的一种基于pe和tvf-emd的超宽带雷达生命体征信号去噪算法，其特征在于，所述步骤2的过程如下：步骤2.1：道信号相减法，对雷达回波信号中的背景杂波进行消除的最简单又直接的方法就是用每一行雷达回波信号与前一行雷达回波信号相减，在这一过程中保持恒定不变的量就会被消除，对于雷达回波矩阵s(m,n)，m为雷达回波延迟时间的采样点数，也就是快时间点数；n为扫描时间的采样点数，也叫做慢时间点数；假设快时间序列为y
j
，j＝1,2,3
…
,n，其中n为采样点数，则道信号相减法的结果y
′
j
为：y
′
j
＝y
j+1-y
j
ꢀꢀꢀꢀ
(6)步骤2.2：减平均法，对雷达回波数据沿扫描时间方向取平均值，即直流分量；直流分量由静止物体所产生，由于静止物体与探测装置间的距离固定不变，而人体的呼吸信号与心跳信号是呈周期性变化的，因此可以用原始数据减去这一平均值就可以达到去除直流分量的效果，则消除直流分量的结果为：其中，可近似看为背景杂波中的直流分量；步骤2.3：线性趋势抑制，根据扫描时间方向上的线性最小二乘拟合来估计所得到的背景杂波中具有潜在性的线性趋势，随后从雷达回波数据中减掉该线性趋势，就可以达到去除背景杂波的目地，去除结果的计算公式如下：w＝ω
t-x(x
t
x)-1
x
t
ω
t
ꢀꢀꢀꢀ
(8)其中，x＝[x1,x2]，x1＝[1,2,
…
,n]
t
，x2＝[1,1,
…
,1]
t
；步骤2.4：自动控制增益，选择2λ+1的时间窗，并根据其中的能量计算相应的增益系数，来完成自适应控制的目的；设第n1帧回波信号为r(τ,n1)，则：r
e
(t,τ)＝g
m
(t,τ)
×
r(t,τ)
ꢀꢀꢀꢀ
(10)其中g
m
(t,τ
i
)为增益系数，r
e
(t,τ)为处理后的人体生命体征信号；步骤2.5：快时间滤波，根据小波变换的小波阈值去噪是一种常见的信号去噪方式，如何选择合适的阈值是雷达回波信号去噪的关键；采用广义交叉验证函数来确定阈值，该方法与传统的d-j阈值相比无需依赖输入输出数据，具有更大的优势，gcv公式表达式如下：其中n表示小波系数的数量，n0表示小波系数为0的数量，w表示含噪声的回波信号，w
t
表示降噪后的回波信号；步骤2.6：慢时间滤波，采用巴特沃斯低通滤波器滤去高频的噪声信号，它的平方幅度
函数定义为：其中n为滤波器的阶数，ω
c
为3db截至频率，ε是控制带通波纹幅度的参数，由于人体呼吸频率范围是0.1-0.5hz，心跳频率范围是0.8-3.0hz，所以其频带选择为0.1-3.0hz。4.根据权利要求1所述的一种基于pe和tvf-emd的超宽带雷达生命体征信号去噪算法，其特征在于，所述步骤3的过程如下：假设快时间为m1，对长度为n的慢时间序列x[m1,i]＝{y(1),y(2),
…
,y(i)}，i＝1,2,
…
,n，y∈w
m
×
n
，进行相空间重构可得：其中τ
d
为延迟时间，q为嵌入维度，k＝n-(q-1)τ
d
。该矩阵的每一行都是重构分量，共有k个重构分量，将第j个重构分量中的元素按升序排列，可以得到新的重构分量：x(i+(j
1-1)τ
d
)≤x(i+(j
2-1)τ
d
)≤
…
≤x(i+(j
q-1)τ
d
)
ꢀꢀꢀꢀꢀ
(14)若其中存在相等的分量，则按j
i
中下标i的大小进行排列即j1j2的排列顺序，对于任意时间序列经过空间重构后得到的重构矩阵中的任意分量可以得到一组序列：s(l)＝(j1,j2,
…
,j
q
)
ꢀꢀꢀꢀꢀ
(15)其中l＝1,2,
…
,k，k≤q！，m维相空间映射可以产生m！种s(l)排序；计算每种序列s(l)不同排列顺序的出现概率p
k
，概率和为1，则按照shannon熵的形式，时间序列x[m1,i]的k种不同符号序列的排列熵可以定义为：当即每种排列顺序的概率存在且相等时排列熵的值最大为ln(m！)，因此对排列熵h
p
进行归一化处理可得：h
p
的值表示慢时间序列x[m1,i]的复杂程度，h
p
的值越小，时间序列变化越规则，反之越随机；在雷达回波矩阵中，目标位置范围在慢时间方向上的数据变化比较规则对应的人体目标的pe值较小，可以通过寻找pe值结果的最低点来判断人体目标的位置p
pos
，公式表示如下：range＝(v
×
p
pos
×
t
f
)/2
ꢀꢀꢀꢀ
(18)其中v＝3
×
108m/s，t
f
为快时间方向上的采样间隔；
假设人体胸腔的横向距离为d，计算得到胸腔距离在信号中所占用的点数为：进而可以得到生命体征信号矩阵δ的表达式为：5.根据权利要求1所述的一种基于pe和tvf-emd的超宽带雷达生命体征信号去噪算法，其特征在于，所述步骤6的过程如下：步骤6.1：局部截止频率b-样条拟合滤波器的截止频率是随着时间变化而发生改变的，它主要通过构造多项式样条的方式来逼近输入信号，可表示如下形式：其中β
n
(t)为b-样条函数，c(k)为b-样条相关系数，为逼近误差，用b-样条阶次n和节点m来确定b-样条相关系数时可使逼近误差最小，公式表示为：其中[]
↑
m
为上采样操作，即在每两个采样点之间插入m个采样点；根据公式(21)，c(k)可表示为：其中[]
↓
m
为下采样操作，即每隔m个点进行采样；为预滤波器，公式表示如下：综上，b-样条拟合器可表示为如下形式：根据上述公式可将b-样条逼近看作一种特殊形式的低通滤波器，通常情况下，用节点m来确定b-样条滤波器的局部截止频率，信号执行时变滤波的时间则由截止频率来确定，进行时变滤波经验模态分解的具体实现步骤如下：步骤6.1.1：利用hilbert变换计算输入信号x(t)的瞬时幅值a(t)与瞬时频率利用hilbert变换计算输入信号x(t)的瞬时幅值a(t)与瞬时频率利用hilbert变换计算输入信号x(t)的瞬时幅值a(t)与瞬时频率其中x(t)为输入信号x(t)的hilbert变换；步骤6.1.2：确定瞬时幅值a(t)的局部极大值序列a({t
max
})与局部极小值序列a({t
min
})z(t)作为多分量信号的解析信号可表示为：
由此可得：由此可得：其中表示第i个分量的瞬时相位，a
i
(t)表示第i个分量的幅值；假设在t
min
处取得局部最小值可得：将公式(31)代入(29)和公式(30)可得：a(t
min
)＝|a1(t
min
)-a2(t
min
)|
ꢀꢀꢀꢀ
(32)根据导数运算可得：a
′1(t
min
)-a'2(t
min
)＝0
ꢀꢀꢀꢀ
(34)同理，可得到如下公式：a(t
max
)＝|a1(t
max
)+a2(t
max
)|
ꢀꢀꢀꢀ
(36)a
′1(t
max
)+a'2(t
max
)＝0
ꢀꢀꢀꢀ
(38)步骤6.1.3：分别对a({t
max
})与a({t
min
})进行插值拟合可得局部极值拟合函数β1(t)与β2(t)并计算瞬时均值函数a1(t)与瞬时包络函数a2(t)(t)步骤6.1.4：计算瞬时频率分量和和和其中η1(t)与η2(t)分别为对进行插值后的结果；
步骤6.1.5：计算局部截止频率计算局部截止频率步骤6.1.6：重新对局部截止频率进行调整来解决间歇问题局部截至频率可能会受到噪声的影响，因此对进行重构可得到一个新的信号：步骤6.2：局部均值函数利用时变滤波器对输入信号进行滤波来获取局部均值，并将f(t)的极值点作为构造时变滤波器的节点，可以使滤波器的截止频率与局部截止频率保持一致，对输入信号x(t)进行b-样条插值逼近可以得到逼近结果m(t)；步骤6.3：停止准则判断准则如下：判断准则如下：判断准则如下：其中b
l
(t)为瞬时带宽，为加权均值瞬时频率，ξ为给定的带宽阈值，若满足公式(45)则该信号为局部窄带信号；根据人体的呼吸频率为0.1hz-0.5hz，心跳频率为0.8hz-3.0hz可将回波信号分为n个imf分量并在频域上对每个分量进行呼吸和心跳能量百分比的计算，计算公式如下：计算公式如下：其中q(j)为第j个imf的频域能量，q
r
(j)为第j个imf呼吸频域能量，q
h
(j)为第j个imf心跳频域能量，δ
r
与δ
h
分别为呼吸和心跳的判断阈值。6.根据权利要求1所述的一种基于pe和tvf-emd的超宽带雷达生命体征信号去噪算法，其特征在于，所述步骤7的过程如下：将符合要求的imf呼吸分量与心跳分量进行重构，公式如下：s
r
(t)＝∑imf(t)
ꢀꢀꢀꢀ
(50)s
h
(t)＝∑imf(t)
ꢀꢀꢀꢀ
(51)。

技术总结
本发明公开了一种基于PE和TVF-EMD的超宽带雷达生命体征信号去噪算法，属于雷达生命探测技术领域。该方法首先采用道信号相减法、减平均法、线性趋势抑制、自动控制增益、快时间滤波、慢时间滤波等方法对原始雷达回波信号中的杂波与噪声信号进行预处理；然后通过对慢时间方向上数据的峰度、标准差、方差等统计特征进行分析，采用PE算法来检测人体目标的位置。该优化算法可以便捷地、精准地定位信号发生突变的时刻，并能放大信号的细微变化。采用TVF-EMD算法对提取后的生命体征信号进行分解，计算分解后各IMF分量在呼吸频带与心跳频带内各自的能量占比，选择合适的信号进行重构，最后利用FFT变换对重构后的信号进行频谱分析。利用本发明所提出的雷达生命信号去噪方法，能准确定位目标并提取生命信息，去噪结果十分可观。去噪结果十分可观。去噪结果十分可观。


技术研发人员：李鑫 周婧 杨桢 吴方泽 段雨昕 李昊
受保护的技术使用者：辽宁工程技术大学
技术研发日：2023.04.06
技术公布日：2023/7/25