基于改进变分模态分解的生命体征检测

doi:10.13232/j.cnki.jnju.2022.04.012

基于改进变分模态分解的生命体征检测

韩宇, 张兴敢^,

南京大学电子科学与工程学院，南京，210023

Vital signs detection based on improved variational mode decomposition

Han Yu, Zhang Xinggan^,

School of Electronic Science and Engineering, Nanjing University, Nanjing, 210023, China

通讯作者: E⁃mail：zhxg@nju.edu.cn

收稿日期: 2022-04-28

基金资助:

国家自然科学基金. 61976113

Received: 2022-04-28

摘要

变分模态分解方法常被用于生命体征信号检测，针对该方法模态个数确定存在主观性的问题，提出一种改进变分模态分解方法进行非接触式测量生命体征信号.利用毫米波雷达获得生命体征回波信号，经过降噪、目标距离单元定位等处理获得生命体征信号.生命体征信号成分主要包括心跳与呼吸信号.其中心跳信号能量小，受呼吸信号影响大，本方法利用余量与模态分量重构信号的自相关系数确定最优分解模态数量，避免了模态数量的过分解和欠分解，达到自适应分离生命体征信号的目的.实验结果表明，与变分模态分解方法相比，该方法可自适应确定模态分解个数.

关键词： 生命体征检测 ; 毫米波雷达 ; 变分模态 ; 模态分解

Abstract

Variational modal decomposition method is often used for vital sign signal detection. Aiming at the problem of subjectivity in determining the number of modes in this method,this paper proposes an improved variational modal decomposition method for non⁃contact measurement of vital sign signals. Here,the millimeter wave radar is used to obtain the vital sign echo signal,and the vital sign signal is obtained through noise reduction,target distance unit positioning and other processing. The vital sign signal components mainly include heartbeat and respiration signals. Among them,the energy of the heartbeat signal is small and is greatly affected by the breathing signal. The method proposed uses the autocorrelation coefficient of the residual and modal components to reconstruct the signal to determine the optimal number of decomposed modes,which avoids over⁃decomposition and under⁃decomposition of the number of modes. The purpose of adaptively separating vital sign signals is achieved. The experimental results show that,compared with the variational modal decomposition method,the proposed method can adaptively determine the number of modal decompositions.

Keywords： vital signs detection ; millimeter wave radar ; variational modes ; mode decomposition

PDF (1591KB) 元数据多维度评价相关文章导出 EndNote| Ris| Bibtex 收藏本文

本文引用格式

韩宇, 张兴敢. 基于改进变分模态分解的生命体征检测. 南京大学学报（自然科学）[J], 2022, 58(4): 680-688 doi:10.13232/j.cnki.jnju.2022.04.012

Han Yu, Zhang Xinggan. Vital signs detection based on improved variational mode decomposition. Journal of nanjing University[J], 2022, 58(4): 680-688 doi:10.13232/j.cnki.jnju.2022.04.012

非接触式、远距离场景下的人体生命体征检测是目前研究的热点，在地震救灾、健康检测等领域有广泛应用.毫米波雷达因其频带宽、波束窄、分辨率高、便于携带等特点，常用于生命体征探测.同时毫米波信号在空气中传播衰减小，可穿透衣物、废墟等非金属材料，在面对身体严重烧伤或烫伤、有传染性疾病以及重度感染时进行非接触生命探测和健康监视更具有优势^［1］.

经验模态分解（Empirical Mode Decomposition，EMD）和变分模态分解（Variational Mode Decomposition，VMD）是常用的生命体征信号分离方法.EMD方法无严格的数学理论基础，存在模式混叠问题.2009年EEMD（Ensemble Empirical Mode Decomposition）方法被提出^［2-3］，改善了EMD分解中的模式混叠问题.针对EEMD方法噪声去除不完全的问题，CEEMD（Complete Ensemble Empirical Mode Decomposition）方法在2011年被提出^［4］.Dragomiretskiy and Zosso^［5］提出的VMD方法是一种完全非递归的信号分解算法，相较于EMD算法分解过程中模态分量个数的不确定性，VMD方法可人为规定模态分解个数^［6］.VMD方法的模态个数确定具有主观性，影响了信号分解的可靠性，针对VMD方法可靠性低的问题，本文提出了改进的VMD方法，该方法可自适应确定所分解的模态分量个数.

本文首先分析了人体生命信号特征及调频连续波雷达原理，利用AWR1642EVM毫米波雷达开发平台获得被测目标回波信号并进行信号预处理，最后对仿真信号和实测信号进行算法验证.针对毫米波雷达获得的目标静态生命体征信号，本文主要研究呼吸和心跳频率的提取.实验结果表明，本方法解决了VMD的不可靠性，避免了人工设置模态个数产生的误差，比EMD分离出的呼吸和心跳信号更纯净.

1 毫米波雷达检测信号理论基础

1.1　人体体征信号特点

本文将生命体征活动近似看作由呼吸和心跳共同引起的胸腔表面的位移振动^［7］.在目标平静时，健康成年人心跳次数一般为每分钟60~100次，在运动过程中，心率比平时速度快.对于呼吸频率，普通成年人的呼吸频率为每分钟12~20次.表1为成人体征信号相关指标，从表中可看出人体生命信号是微弱和低速的运动信号.

表1 体征信号相关指标

Table 1 Signs and signals related indicators

体征信号	频率（Hz）	幅度（mm）
呼吸速率（成人）	0.1~0.5	1~12
心跳速率（成人）	0.8~3.0	0.1~0.5

新窗口打开| 下载CSV

生命信号一般包括心跳信号、呼吸信号以及背景噪声，心跳和呼吸信号可近似看作频率稳定的正弦波信号^［8］，生命信号可由式（1）表示：

\begin{array}{l} x (t) = \\ A_{h} s i n (2 π f_{h} t + φ_{h}) + A_{b} s i n (2 π f_{b} t + φ_{b}) + z \end{array}

(1)

其中，A_h为心跳信号振幅，f_h为心跳信号频率，A_b为呼吸信号振幅，f_b为呼吸信号频率，z为噪声.

1.2　FMCW雷达工作原理

毫米波FMCW雷达连续发射的线性调频脉冲信号为 $s (t)$ ：

s (t) = e^{j (2 π f_{c} t + π \frac{B}{T} t^{2})}

(2)

其中， $f_{c}$ 为起始频率，B为雷达调制带宽，T为线性调频脉冲信号持续时间.

与雷达间隔R处某单一目标接收到电磁波并反射回波，雷达与目标间存在传播距离，传播过程中产生 $t_{d}$ 的时延，则雷达接收回波信号 $r (t)$ 为：

r (t) = e^{j [2 π f_{c} (t - t_{d}) + π \frac{B}{T} {(t - t_{d})}^{2}]}

(3)

其中 $t_{d} = \frac{2 R}{c}$ ，c为光速.

将发射信号 $s (t)$ 与回波信号 $r (t)$ 进行混频获得差频信号 $b (t)$ ，如式（4）所示：

b (t) = s (t) r (t) = e^{j (\frac{2 π B t_{d}}{T} t + 2 π f_{c} t_{d} - π \frac{B}{T} {t_{d}}^{2})}

(4)

$r (t)$ 仍为正弦信号，代入 $t_{d} = \frac{2 R}{c}$ ，时延二次项 $t_{d}^{2}$ 极小，因此可忽略，推导出差频信号近似表达，如式（5）所示：

b (t) \approx e^{j (4 π \frac{B R}{c T} t + \frac{4 π}{λ} R)} = e^{j (f_{b} t + ϕ_{b})}

(5)

其中，频率 $f_{b} = \frac{4 π B R}{c T}$ ，相位 $ϕ_{b} = \frac{4 π R}{λ}$ ，λ为波长.

由式（5）可看出，差频信号相位部分 $ϕ_{b}$ 与目标径向距离R为线性关系，目标的极小位移会引起相位的较大变化，为了测量小尺度振动，在目标距离单元测量FMCW信号随时间的相位变化^［9］.连续测量值的相位随物体移动距离ΔR的变化如式（6）所示.可以看出，相位变化很好地描述了目标微小位移，如当 $λ = 4 m m$ 时，目标位移 $Δ R = 1 m m$ ，其相位改变为 $Δ ϕ_{b} = π$ .

Δ ϕ_{b} = \frac{4 π}{λ} Δ R

(6)

2 信号预处理

2.1　动目标检测

在实际测量中，有许多静止分量的干扰信号，当这些干扰信号能量高于待检测目标能量时，真实目标被淹没，检测结果会发生严重偏差，甚至无法正确提取人体生命信号.本节采用动目标检测方法，主要利用脉冲去除静态杂波分量.在雷达回波中，每个脉冲的干扰成分相同，运动目标相对雷达的距离变化会引起相位改变.前后脉冲对消，可滤除信号中的静止杂波分量，留下运动目标分量.将原始矩阵每行的每个数减去该行的平均值，得到原始信号滤除静止杂波分量后的数据，该矩阵去除了静止杂波分量，保留了运动信息.实验结果见图1和图2.

图1

新窗口打开| 下载原图ZIP| 生成PPT

图1 原始雷达回波信号

(a)回波脉冲三维距离像；(b)回波脉冲二维距离像

Fig.1 Raw radar echo

图2

新窗口打开| 下载原图ZIP| 生成PPT

图2 去除静态杂波后的雷达回波信号

(a)回波脉冲三维距离像；(b)回波脉冲二维距离像

Fig.2 Radar echo signal after removing static clutter

图1为原始雷达回波的距离像，可以看出静态杂波分量比较强.图2为雷达回波信号经过动目标检测后的距离像，可看出静态杂波得到明显抑制.

2.2　距离门切换

从图2可知，提取目标所在距离单元位置是通过距离单元能量确定的.实验中，目标胸壁起伏、轻微人体位移等均会导致距离单位的跳变，导致结果存在误差甚至错误.本节通过距离门切换解决上述问题.在实验中由于人体轻微晃动有可能使距离门改变，从而引起实验误差.为了更准确地计算，本文采用滑动窗口确定能量最大的距离单元.选择窗长为N个点，滑窗的步进长度为M个点（慢时间采样频率为M Hz），可用此方法来每秒计算目标实时距离单元.

每个滑窗内的数据为 $k \times N$ 的矩阵，其中，k为距离单元个数.由式（7）可计算出该滑窗内能量最大的距离门，即求此滑窗内所有慢时间采样值的平方和，R_n，i 是经过MTI后数据矩阵的第n行i列值，n为距离门位置，i为窗内慢时间采样点位置.如此可计算滑窗内每个距离门的能量值，取最大值为滑窗能量最大距离门，此距离门即被测对象距离雷达系统距离.

E_{n} = \sum_{i = 1}^{N} R_{n, i}^{2}

(7)

对所有数据进行步进滑窗，可得能量最大距离门组.取距离门组中的众数，看作整个实验中能量最大距离门，若目标不动，可将此距离门中信息作为后续携带人体生命体征信号.

2.3　相位信号处理

对经过距离门切换得到目标信号并提取其相位信息，如图3a所示.由于相位值在 $[- π, π]$ ，因此对相位进行解缠绕得到图3b的波形.心跳信号比呼吸信号的能量小，因此需对解缠绕后的相位进行差分运算来增强心跳信号.从图3c可明显看出存在脉冲噪声，中值滤波算法简单，可有效抑制脉冲噪声，保护信号边缘不被模糊，因此可用于消除脉冲噪声^［9-10］，结果如图3d所示.

图3

新窗口打开| 下载原图ZIP| 生成PPT

图3 相位处理波形

(a)相位展开前；(b)相位展开后；(c)相位展开差分后；

(d)相位展开差分中值滤波后

Fig.3 Phase processing waveform

3 生命体征信号分离

3.1　VMD方法

EMD方法在研究瞬时频率时，提出了本征模态函数（Intrinsic Mode Function，IMF）的概念，是一种将信号分解为一组IMF的方法，其核心是得到一组IMF分量^［11］.EMD算法存在很多局限性，总结如下：

（1）EMD是一种根据经验分解信号的算法，在计算过程中没有严格的数学理论基础；

（2）抗噪性能较差，噪声鲁棒性较差；

（3）容易出现模式混叠.

针对EMD方法存在的上述缺点，提出VMD算法代替EMD算法处理信号.该算法的目的是将原始数据分解为k个模态分量 $u_{k} (t)$ ，每个模态分量都在中心频率ω_k 附近，存在一定的带宽.这样就将问题转化为使所有模态分量的带宽和最小，VMD算法根据经验规定IMF个数.设需要分解模态分量个数为k，分解具体步骤如下^［12］：

（1）对每个模态信号 $u_{k} (t)$ 进行希尔伯特变换：

(δ (t) + \frac{j}{π t}) * u_{k} (t)

(8)

（2）将每个模态经过希尔伯特变换后的频谱搬移至基带：

[(δ (t) + \frac{j}{π t}) * u_{k} (t)] e^{- j ω_{k} t}

(9)

（3）计算式（9）信号梯度平方的L²范数来估计每个IMF的带宽，受约束的变分问题如下式：

\begin{array}{l} \underset{\{u_{k}\}, \{ω_{k}\}}{m i n} \{{\sum_{k} ‖\partial t [(δ (t) + \frac{j}{π t}) * u_{k} (t)] e^{- j ω_{k} t}‖}_{2}^{2}\} \\ s . t . \sum_{k} u_{k} = f \end{array}

(10)

其中， $u_{k} (t)$ 是所有模态，ω_k 是模态对应的中心频率.

为不受约束，使用二次惩罚因子α和拉格朗日乘子 $λ (t)$ .增广拉格朗日如下：

\begin{array}{l} L (\{u_{k}\}, \{ω_{k}\}, λ) : = \\ α \sum_{k} {‖t [(δ (t) + \frac{j}{π t}) * u_{k} (t)] e^{- j ω_{k} t}‖}_{2}^{2} + \\ {‖f (t) - \sum_{k} u_{k} (t)‖}_{2}^{2} + 〈λ (t), f (t) - \sum_{k} u_{k} (t)〉 \end{array}

(11)

利用交替方向乘子法（ADMM）解决变分问题，通过更新 $u_{k}^{n + 1}$ ， $ω_{k}^{n + 1}$ 以及 $λ^{n + 1}$ 寻找拉格朗日鞍点.获得 $u_{k}^{n + 1} (ω)$ 和 $ω_{k}^{n + 1} (ω)$ 最优解，代入ADMM得到完整的VMD算法，流程如下：

（1）初始化参数：初始化 $｛u_{k}^{1}｝｛ω_{k}^{1}｝｛λ^{1}｝$ ，并将n初始化为0；

（2）利用下式对IMF从1到k更新u_k 和ω_k ：

\begin{array}{l} u_{k}^{n + 1} (ω) = \\ \frac{f (ω) - \sum_{i < k} {u_{i}}^{n + 1} (ω) - \sum_{i > k} {u_{i}}^{n} (ω) + \frac{λ^{n} (ω)}{2}}{1 + 2 α {(ω - {ω_{n}}^{k})}^{2}} \end{array}

(12)

ω_{k}^{n + 1} = \frac{\int_{0}^{\infty} ω {|u_{k}^{n + 1} (ω)|}^{2} d ω}{\int_{0}^{\infty} {|u_{k}^{n + 1} (ω)|}^{2} d ω}

(13)

（3）更新λ：

λ^{n + 1} (ω) = λ^{n} (ω) + τ [f (ω) - \sum_{k} u_{k}^{n + 1} (ω)]

(14)

（4）判断是否收敛，条件如下：

\sum_{k} \frac{{‖u_{k}^{n + 1} - u_{k}^{n}‖}_{2}^{2}}{{‖u_{k}^{n}‖}_{2}^{2}} < ε

(15)

其中， $ε = 10^{- 6}$ .若收敛则循环停止，得到IMF，若不收敛，则返回第二步，并重复上述步骤.

3.2　改进VMD

与EMD算法分解过程中IMF个数的不确定性比较，VMD算法可人为规定IMF的个数，该个数根据经验规定.当规定的分解模态数量过小，会导致欠分解；若规定的分解模态数量过大，则会导致过分解.因此本节采用改进VMD算法确定分解模态的数量^［13］.

互相关系数是表征两个信号相关程度的量.当两信号相关程度越高，则它们的互相关系数值越大^［14-15］.当分解的IMF数量最优时，利用VMD对信号分解后的余量是噪声，余量与IMF互不相关^［16］.因此，本文采用余量与IMF叠加信号的互相关系数确定最优IMF的数量，互相关系数值越小，选取的值越优.设信号经过VMD分解后的余量为 $m (n)$ ，分解后所有的IMF重构的信号为 $f (n)$ ，则两者归一化互相关函数为：

ρ (k) = \frac{R_{m f} (k)}{\sqrt[]{R_{f f} (0) R_{m m} (0)}} = \frac{\sum_{n = 0}^{N} f (n) c (n - k)}{\sqrt[]{R_{f f} (0) R_{m m} (0)}}

(16)

其中， $R_{m f} (k)$ 为两者的互相关函数，N为序列长度，k为时间差， $R_{m m} (0)$ 为余量的自相关函数， $R_{f f} (0)$ 为模态重构信号的自相关函数.当 $k = 1$ 时，信号间的相关性最大，用下式确定最优分解模式的数量：

ρ_{m} = \frac{R_{m f} (1)}{\sqrt[]{R_{f f} (0) R_{m m} (0)}}

(17)

4 结果与分析

4.1　仿真信号分离算法比较

由1.1节人体体征信号特点，可建模人体生命信号.呼吸和心跳信号近似看作正弦信号，为便于分析，本节采用正弦函数建模呼吸和心跳信号.设定心跳频率为1.1 Hz、呼吸频率为0.3 Hz.由于心跳信号幅度比呼吸信号幅度微弱很多，所以设心跳信号振幅为0.001 m，呼吸信号振幅为0.01 m，即：

呼吸信号建模为

y_{b} (t) = A_{b} s i n (2 π \cdot f_{b} \cdot t)

(18)

其中，A_b=0.01 m，f_b=0.3 Hz.

心跳信号建模为

y_{h} (t) = A_{h} s i n (2 π \cdot f_{h} \cdot t)

(19)

其中，A_h=0.001 m，f_h=1.1 Hz.

实际测量中存在噪声，因此人为加入信噪比为10 dB的高斯白噪声.仿真建模信号时域波形如图4所示；添加信噪比为10 dB的高斯白噪声的合成信号时域波形如图5所示.

图4

新窗口打开| 下载原图ZIP| 生成PPT

图4 建模呼吸信号时域波形图

(a)建模呼吸信号；(b)建模心跳信号

Fig.4 Modeling respiratory signal time domain waveform

图5

新窗口打开| 下载原图ZIP| 生成PPT

图5 合成信号时域图

Fig.5 Synthetic signal time domain plot

从图4和图5可看出，合成信号波形接近呼吸信号波形，是由于呼吸信号振幅远大于心跳信号振幅，对合成信号添加高斯白噪声，使仿真信号更接近真实值.

对上面建模得到的单人合成生命体征信号用EMD和VMD算法分离信号，得到的结果如图6和图7所示.

图6

新窗口打开| 下载原图ZIP| 生成PPT

图6 EMD分离结果

Fig.6 EMD separation results

图7

新窗口打开| 下载原图ZIP| 生成PPT

图7 VMD分离结果

Fig.7 VMD separation results

EMD算法分离结果如图6所示，由图可见，IMF1是无用的噪声信号，IMF2中包含了心跳信号频率，IMF3包含了呼吸信号频率.

下面讨论用VMD算法分离呼吸和心跳信号.因为建模信号仅有呼吸和心跳信号频率，因此将所需分离模态个数 $k$ 人为设置为2， $τ = 0$ ， $α = 900$ ， $ε = 10^{- 7}$ ，分离结果如图7所示.图7a和图7b分别为IMF1的时域波形和频域波形，由频域波形可看出分离的信号频率为0.3 Hz，即呼吸信号；图7c和图7d为IMF2的时域波形和频域波形，可看到IMF2信号为频率1.1 Hz的心跳信号.由图7的VMD分离结果可明显看出，仿真信号中的呼吸和心跳信号被成功分离，时域波形振幅虽与原信号存在差异，但不影响结果，分离出的呼吸和心跳信号频率正确.

对比图6和图7，VMD算法和EMD算法相比，可以分离出更真实的呼吸和心跳信号.EMD算法分离的IMF数量是不定的，需要迭代至条件不满足退出循环，由于IMF数量的不确定性，对分解后的数据处理影响较大，若不观察分解结果，无法确定呼吸和心跳信号与IMF分量的对应关系，VMD算法可以人为设定分解IMF的个数，在该仿真实验中IMF数量设置为2，呼吸和心跳分量分别对应IMF1和IMF2.

4.2　实测信号测试

本节通过实际测量信号验证改进VMD算法对确定分解模态数量的可行性.采用德州仪器（TI）公司的AWR1642毫米波雷达开发平台，被测目标为一名成年女性呈坐姿平静状态接受检测.雷达发射FMCW波，遇到目标后反射回波信号，对获得信号进行距离FFT，通过距离门切换确定目标所在距离单元，进行相位提取、解缠绕、差分运算及中值滤波，处理后的相位信号幅频图见图8，AWR1642参数见表2.

图8

新窗口打开| 下载原图ZIP| 生成PPT

图8 处理后的相位信号幅频图

Fig.8 Amplitude⁃frequency diagram of the processed phase signal

表2 AWR1642参数

Table 2 AWR1642 parameters

参数	值	备注
开始距离（m）	0.3	被预测目标的开始和结束距离，程序会搜索范围内最大峰值，并假设该峰值对应受试者位置
结束距离（m）	0.9	被预测目标的开始和结束距离，程序会搜索范围内最大峰值，并假设该峰值对应受试者位置
距离单元长度（m）	0.0337	一个距离单元的间隔长度
接收天线	1	要处理的接收天线数量，当前实现来自单个接收天线数据
采样频率 $F_{s}$ （Hz）	20	每1 s在慢时间采样20个点

新窗口打开| 下载CSV

对实际测量的目标信号进行VMD分解，根据生命体征信号特点，设定分解IMF数量k依次为2，3，4，不同k值的余量与模式重构信号归一化互相关系数如表3所示.从表中可看出k=3时，归一化互相关系数为最小值，从图8中可看出实测信号有三个谱峰，因此采用互相关系数法确定的模式数量与实际结果一致.

表3 不同分解模式数量下余量与模式重构信号的互相关系数

Table 3 The cross⁃correlation coefficient between the residual and the mode⁃reconstructed signal under different number of decomposition modes

模式数量	k=2	k=3	k=4
归一化互相关系数	0.16	0.10	0.13

新窗口打开| 下载CSV

对该实测信号进行 $k = 3$ 的VMD模态分解，结果见图9.VMD对信号进行分解后，IMF1频率为0.28 Hz，对应呼吸信号频率，IMF2频率为0.5 Hz，对应呼吸谐波信号频率，IMF3频率为1.29 Hz，对应心跳信号频率，实验结果表明改进VMD方法可有效估计模态分解数量.

图9

新窗口打开| 下载原图ZIP| 生成PPT

图9 VMD分离后各模式时域及频域波形

Fig.9 Time⁃domain and frequency⁃domain waveforms of each mode after VMD separation

5 结论

本文基于毫米波雷达对生命体征进行检测，针对VMD算法分解模态个数需人为设定的缺点，提出了改进VMD分解算法提取呼吸和心跳信号频率.通过仿真和实测实验结果验证，结论如下：

（1）仿真数据实验表明：在已知信号成分数量时，EMD算法模态分解个数不确定，容易产生模态混叠，VMD算法可以避免EMD方法发生的模态混叠，可以分离出更真实的呼吸和心跳信号，得到人体呼吸和心跳频率.

（2）提出一种基于互相关函数法确定模态分解个数的方法，对实际测量的数据进行算法验证，结果表明：本文方法避免了VMD方法按经验选取模态个数产生结果误差的缺点.

参考文献

原文顺序

文献年度倒序

文中引用次数倒序

被引期刊影响因子

[1]

Chen

K M

, Misra

, Wang

,et al.

An x⁃band microwave life⁃detection system

IEEE Transactions on Biomedical Engineering,1986,BME⁃33(7)：697-701.