基于VMD-SWT联合算法的FMCW雷达生命体征检测

何鹏宇,卓智海

(北京信息科技大学 信息与通信工程学院,北京 100101)

0 引言

生命体征能够反映人体的健康状况,呼吸率和心率是重要的生命体征,可作为诊断某些疾病的重要参数。生命体征检测技术可分为接触式和非接触式两类。传统的生命体征检测方法主要采用接触式监测呼吸和心率[1-2],但接触式设备会限制被测者的行动,使被测者感到不适。而非接触式雷达检测技术可以远距离监测呼吸率和心率,符合人们对生命体征检测技术智能化和舒适度的要求,更适合于长期监测,便于健康管理。

目前用于非接触式生命体征检测的雷达系统主要有三种,分别是连续波(continuous wave,CW)雷达、超宽带(ultra-wideband,UWB)雷达和调频连续波(frequency modulated continuous wave,FMCW)雷达[3-5],文献[6]比较了三种雷达系统的区别。FMCW雷达的发射频率是随时间线性调制的,可以在探测目标距离的同时,利用多普勒效应检测胸腔位移。通过提取中频信号的频率及相位,可以获得被测对象的准确距离和位移信息[3-4]。与微波雷达相比,毫米波雷达的带宽更宽、波长更短,因此具有更高的距离分辨率和灵敏度,适合于非接触式生命体征检测[7]。FMCW毫米波雷达在距离测量和微小位移测量方面都具有优异的性能[5],因此本文采用77 GHz FMCW毫米波雷达进行非接触式生命体征检测。

在处理生命体征信号的过程中,由于呼吸引起的胸部位移远大于心跳引起的胸部位移,较弱的心跳信号容易受到呼吸谐波和杂波的干扰,特别是当心跳信号和呼吸谐波在频域上接近或重叠时,呼吸谐波的频率可能被误认为心率。因此,如何处理呼吸信号引起的谐波干扰是心率检测的主要问题[8-9]。文献[10]提出了一种基于经验模态分解(empirical mode decomposition,EMD)的生命体征检测算法。然而,EMD的模态混合和端点效应问题导致心跳和呼吸谐波在一个模态中共存,影响了心跳信号检测的准确性[11]。文献[12]提出了集合经验模态分解(ensemble empirical mode decomposition,EEMD),通过填充高斯白噪声来提高EMD的性能,但模态分量中的残留噪声会降低心率估计的准确性。文献[13]提出了基于变分模态分解(variational mode decomposition,VMD)的生命体征检测,并通过实验验证了在相同的测试环境下,VMD的效果要优于EMD。但该方法不能提供心跳信号的时频信息,当心率在短时间内出现波动时,检测结果的误差较大。文献[14]提到同步压缩变换在医学领域中用于对脑电图或心电图的特征提取。其中同步压缩小波变换(synchrosqueezed wavelet transform,SWT)是常用的同步压缩变换之一,能够提供高分辨的时频图,文献[15-16]将其用于心电图中的心跳信号时频特征提取。本文将SWT引入了雷达生命体征检测方法中,提出了VMD和SWT的联合算法,经过实验验证,该方法能有效分离呼吸信号和心跳信号,降低呼吸谐波对心率估计的干扰,得到高精度的心跳信号时频曲线。

1 FMCW雷达生命体征检测原理

FMCW雷达发射线性调频信号,正交接收机捕获目标反射的回波信号并与发射信号进行正交混频得到中频信号。提取中频信号的频率和相位可以确定目标与雷达的距离和目标的胸部位移。发射的线性调频信号可以表示为

(1)

式中:AT为信号的幅度;fc为线性调频信号的起始频率;B为带宽;Tc为线性调频信号的持续时间;φ(t)为相位噪声。

接收信号可以认为是发射信号的延迟。若d0为目标到雷达的距离,x(t)为胸部运动引起的位移,则胸部到雷达的距离为R(t)=d0+x(t)。延迟时间为td=2R(t)/c,其中c为光速。则接收到的信号可以表示为

xR(t)=

(2)

式中:AR为接收信号的幅度。

接收信号和发射信号通过两个正交I/Q信道进行混合,然后通过低通滤波器得到中频信号y(t),其可以表示为

y(t)=

(3)

y(t)=ATARexp(j(2πfbt+φb(t)))

(4)

式中:fb=2Bd0/(cTc);φb(t)=4π(d0+x(t))/λc,λc为频率fc对应的波长。

对于身体健康的成年人,心跳的位移幅度为0.1～0.5 mm,呼吸的位移幅度为1～12 mm,心跳频率为0.8～2 Hz,呼吸频率为0.1～0.5 Hz。在相位φb中,胸部位移x(t)和波长λc都处于毫米范围,因此微小的胸部位移也会引起相位φb(t)的明显变化,通过持续发射线性调频信号并提取中频信号相位可以实现对胸部位移x(t)采样,进而估计呼吸率和心率。

2 VMD-SWT联合算法

2.1 概述

基于VMD-SWT联合算法的FMCW雷达生命体征检测流程如图1所示,主要包括目标定位、相位提取、信号分离和频率估计4个步骤。首先,对中频信号的采样数据进行距离快速傅里叶变换(fast Fourier transform,FFT)以获得距离特征(range profile)。然后,根据距离特征的极值确定目标对应的距离单元(range bin)。提取目标所在的距离单元的相位值进行相位展开和相位差分运算。之后对相位差分信号使用VMD算法进行信号分解,得到呼吸信号和心跳信号,通过搜索这两个信号的频谱峰值可以初步估计呼吸率和心率。但心跳信号易受呼吸谐波的干扰,且该方法不能得到心跳信号的时频信息,因此需要使用 SWT进一步处理心跳信号,并采用模极大值法提取小波脊线得到心跳信号的时频曲线,提高对心跳信号的检测精度。

2.2 基于变分模态分解的信号分解方法

变分模态分解算法[17]假设任何信号都是由一系列具有特定中心频率和有限带宽的子信号组成,即本征模态函数(intrinsic mode function,IMF)。IMF定义为调幅调频信号,可以表示为

uk(t)=Ak(t)cos(φk(t))

(5)

式中:Ak(t)和φk(t)分别为第k个模态uk(t)的幅值和相位,幅值Ak(t)和瞬时相位φk′(t)的变化都比较小,且模态uk(t)具有中心频率ωk。模态uk(t)和其中心频率ωk可通过变分问题求解,求解过程需要满足两个条件:①模态带宽之和最小;②模态之和等于原始信号。假设原始信号f(t)被分解成k(正整数)个模态,根据上述条件,约束变分表达式为

(6)

式中:δ(t)为单位冲激函数;*表示卷积。引入拉格朗日乘数λ和惩罚因子α,将约束问题变为非约束问题,使约束最小化问题转为求增广拉格朗日函数L的鞍点,其表达式为

L({uk},{ωk},λ)=

(7)

uk(t)和ωk可以通过交替方向乘子法迭代得到,迭代公式如下:

(8)

(9)

(10)

(11)

(12)

式中:εr和εa分别为相对误差和绝对误差。当两式均满足时,停止迭代。

对相位信号使用VMD方法分解,可获得一组模态分量,根据模态分量的中心频率ωk确定所对应的信号,实现呼吸信号和心跳信号的分离。

2.3 基于同步压缩小波变换的时频分析

同步压缩小波变换以连续小波变换为基础,给定母小波函数ψ(t),对信号s(t)进行连续小波变换,其小波变换系数Ws(a,b)定义为

(13)

根据Plancherel定理,信号s(t)以ψ(t)为母小波函数的连续小波变换Ws(a,b)可以重写为

(14)

(15)

(16)

通过模极大值法可提取同步压缩小波变换系数Ts(ωl,b)的小波脊线。

使用 SWT对VMD分解结果中的心跳信号分量进行处理,并提取小波脊线即可得到心跳信号的时频曲线。

3 实验结果与分析

本文采用德州仪器的毫米波雷达IWR1443,一个帧周期发射一个啁啾(Chirp)信号,将接收信号与发射信号混频得到中频信号,对中频信号采样并进行数据处理。参数配置如表1所示。

表1 参数设置Table 1 Parameter configuration

被测者坐在雷达传感器前,胸部与雷达处于同一水平高度且没有杂物干扰,距离雷达0.3～2 m范围内。实验设备及实验场景如图2所示。同时被测者使用指夹式脉搏血氧仪测量心跳作为参考值。

图2 实验设备及实验场景Fig.2 Experimental equipment and experimental setup

根据实测数据绘制的距离特征如图3所示。图3显示了1 min内测量的1 200帧数据,目标的信息反映在图中较亮的部分。从图中可以看出,目标大约在距离雷达0.7 m的位置。

图3 距离特征Fig.3 Range profile

根据距离特征的极值确定目标所在的距离单元并进行相位提取、相位展开和相位差分,结果如图4所示。图4显示了30 s内目标的胸部位移信息,图中频率低幅度高的分量由呼吸引起,而频率高幅度低的分量主要由心跳引起,同时也存在呼吸信号的谐波以及其他干扰,相位信号的频谱如图5所示。从图5可以看出在0.25 Hz和0.50 Hz处有明显的峰值,分别对应呼吸率和呼吸谐波频率。而心跳频率的峰值不明显,通过求取正常心率范围0.8～2.0 Hz之间的极值,得到峰值为1.17 Hz。

图4 相位信号Fig.4 Phase signal

图5 相位信号频谱Fig.5 Phase signal spectrum

对图4相位信号使用EMD算法进行分解,得到5个分量IMF1～IMF5以及残余分量res,其分解结果如图6左列所示,对每个分量进行傅里叶变换,得到其对应的频谱,结果如图6右列所示。从图中可以看出,IMF3为呼吸信号,IMF2为呼吸谐波信号与心跳信号。EMD的分解模态数不能预先设定,会一直迭代到满足终止条件,模态数不确定导致呼吸信号与心跳信号所处的分量也不确定,对分解后的处理十分不便,且分解结果也不够理想,心跳和呼吸谐波在一个模态中共存。

而VMD算法可以提前设定模态数量,将其预先设定为4,其分解结果如图7左列所示,每个分量对应的频谱如图7右列所示。从图中可以看出,IMF4对应呼吸信号,IMF3对应呼吸谐波信号,IMF2对应心跳信号。与EMD相比,VMD的分解结果更接近正弦波形,其频谱峰值明显,在提取心跳信号方面明显优于EMD。在IMF2中可以看到2个较明显的峰值分别为1.00 Hz和1.17 Hz,出现两个峰值的原因可能是呼吸信号的谐波干扰,或者是心跳信号的频率在测量时间内有波动,其在两个峰值频率所保持的时间较长。由于VMD的分解结果不能提供时频信息,无法得知心跳信号在测试时间内的变化趋势,因此需要对心跳信号使用SWT进一步处理。

图6 EMD模态分量及频谱Fig.6 EMD mode components and spectrum

图7 VMD模态分量及频谱Fig.7 VMD mode components and spectrum

对VMD的模态IMF2使用SWT进行时频分析,其结果如图8所示。通过模极大值法提取小波脊线,可得到心跳信号的时频曲线,其结果如图9所示。

图8 心跳信号的SWT时频面表示Fig.8 SWT time-frequency representation of heartbeat signal

图9 心跳信号的时频曲线Fig.9 Time-frequency curve of heartbeat signal

将图9的心跳频率转换为每分钟心跳次数,指夹式脉搏血氧仪测量的心率作为参考值ref,两者对比的结果如图10所示。从图中可以看出,两条曲线虽然不能完全重合,但其变化趋势基本一致,且SWT的结果与参考值偏差也较小。因此VMD-SWT联合算法可以获得准确的心跳信号时频曲线。

图10 SWT时频曲线与指夹式脉搏血氧仪检测结果对比Fig.10 Comparison of SWT time-frequency curves with finger pulse oximeter detection results

为了进一步验证所提算法的准确性和有效性,对10名测试者进行了实验测试。每位测试者采集3 min的数据,并从中提取出1 min的数据长度进行呼吸率和心率估计,并使用均方根误差(root mean square error,RMSE)对心率估计进行性能比较。实验结果如表2所示,其中EMD和VMD的结果为心跳信号的频谱峰值所对应的心率,指压式脉搏血氧仪和VMD-SWT的联合算法为1 min内的心率平均值。由表2可以看出,VMD-SWT联合算法的准确性高于其他方法。

表2 VMD-SWT联合算法与其他方法的心率检测结果对比Table 2 Comparison of heart rate detection results between VMD-SWT joint algorithm and other methods 次/min

4 结束语

本文基于FMCW毫米波雷达实现生命体征监测,针对雷达生命体征检测中心跳信号易受到呼吸谐波干扰的现象,提出了基于VMD-SWT的联合算法,对VMD分解结果中的心跳信号使用SWT做进一步处理,获取心跳信号的时频信息,得到更准确的心率估计,并通过实验将该方法与EMD和VMD方法进行对比。实验结果表明,基于VMD-SWT的联合算法能够准确有效地分离呼吸信号和心跳信号,与EMD和VMD方法相比,该方法可以提供准确的时频信息,降低呼吸谐波对检测结果的干扰,提高雷达生命体征检测精度和稳定性。