CN113854992A - 基于77GHz毫米雷达的非接触式精确心率检测方法 - Google Patents
基于77GHz毫米雷达的非接触式精确心率检测方法 Download PDFInfo
- Publication number
- CN113854992A CN113854992A CN202111190201.9A CN202111190201A CN113854992A CN 113854992 A CN113854992 A CN 113854992A CN 202111190201 A CN202111190201 A CN 202111190201A CN 113854992 A CN113854992 A CN 113854992A
- Authority
- CN
- China
- Prior art keywords
- signal
- radar
- frequency
- algorithm
- heartbeat
- Prior art date
- Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
- Pending
Links
- 238000001514 detection method Methods 0.000 title claims abstract description 33
- 238000001228 spectrum Methods 0.000 claims abstract description 35
- 210000000779 thoracic wall Anatomy 0.000 claims abstract description 31
- 238000006073 displacement reaction Methods 0.000 claims abstract description 19
- 230000001121 heart beat frequency Effects 0.000 claims abstract description 8
- 230000033001 locomotion Effects 0.000 claims abstract description 6
- 238000000354 decomposition reaction Methods 0.000 claims description 34
- 238000000034 method Methods 0.000 claims description 33
- 238000005070 sampling Methods 0.000 claims description 22
- 230000008569 process Effects 0.000 claims description 12
- 230000029058 respiratory gaseous exchange Effects 0.000 claims description 10
- 238000007781 pre-processing Methods 0.000 claims description 9
- 230000003044 adaptive effect Effects 0.000 claims description 7
- 230000003190 augmentative effect Effects 0.000 claims description 6
- 230000003595 spectral effect Effects 0.000 claims description 5
- 238000005516 engineering process Methods 0.000 claims description 4
- 239000000284 extract Substances 0.000 claims description 4
- 230000003068 static effect Effects 0.000 claims description 4
- 230000008859 change Effects 0.000 claims description 3
- 230000003111 delayed effect Effects 0.000 claims description 3
- 238000001914 filtration Methods 0.000 claims description 3
- 238000005259 measurement Methods 0.000 claims description 3
- 230000010363 phase shift Effects 0.000 claims description 3
- 230000009467 reduction Effects 0.000 claims description 3
- 230000036391 respiratory frequency Effects 0.000 claims description 3
- 238000012216 screening Methods 0.000 claims description 3
- 230000000276 sedentary effect Effects 0.000 claims description 3
- 239000000126 substance Substances 0.000 claims description 3
- 230000009466 transformation Effects 0.000 claims description 2
- 238000000605 extraction Methods 0.000 abstract description 5
- 238000004364 calculation method Methods 0.000 abstract description 3
- 238000010009 beating Methods 0.000 abstract description 2
- 210000000038 chest Anatomy 0.000 abstract description 2
- 238000004804 winding Methods 0.000 abstract description 2
- 238000004458 analytical method Methods 0.000 description 7
- 238000012544 monitoring process Methods 0.000 description 6
- 238000002474 experimental method Methods 0.000 description 5
- 230000036541 health Effects 0.000 description 5
- 238000011156 evaluation Methods 0.000 description 3
- 238000012545 processing Methods 0.000 description 3
- 230000035945 sensitivity Effects 0.000 description 3
- 238000003745 diagnosis Methods 0.000 description 2
- 238000010586 diagram Methods 0.000 description 2
- 230000003321 amplification Effects 0.000 description 1
- 230000002612 cardiopulmonary effect Effects 0.000 description 1
- 238000006243 chemical reaction Methods 0.000 description 1
- 238000010835 comparative analysis Methods 0.000 description 1
- 238000011161 development Methods 0.000 description 1
- 230000003862 health status Effects 0.000 description 1
- 208000019622 heart disease Diseases 0.000 description 1
- 230000006872 improvement Effects 0.000 description 1
- 230000035772 mutation Effects 0.000 description 1
- 238000003199 nucleic acid amplification method Methods 0.000 description 1
- 230000036285 pathological change Effects 0.000 description 1
- 231100000915 pathological change Toxicity 0.000 description 1
- 230000002093 peripheral effect Effects 0.000 description 1
- 230000000284 resting effect Effects 0.000 description 1
- 230000033764 rhythmic process Effects 0.000 description 1
- 238000000926 separation method Methods 0.000 description 1
- 210000000115 thoracic cavity Anatomy 0.000 description 1
- 238000012795 verification Methods 0.000 description 1
- 210000001835 viscera Anatomy 0.000 description 1
- 210000000707 wrist Anatomy 0.000 description 1
Images
Classifications
-
- A—HUMAN NECESSITIES
- A61—MEDICAL OR VETERINARY SCIENCE; HYGIENE
- A61B—DIAGNOSIS; SURGERY; IDENTIFICATION
- A61B5/00—Measuring for diagnostic purposes; Identification of persons
- A61B5/02—Detecting, measuring or recording pulse, heart rate, blood pressure or blood flow; Combined pulse/heart-rate/blood pressure determination; Evaluating a cardiovascular condition not otherwise provided for, e.g. using combinations of techniques provided for in this group with electrocardiography or electroauscultation; Heart catheters for measuring blood pressure
- A61B5/024—Detecting, measuring or recording pulse rate or heart rate
-
- A—HUMAN NECESSITIES
- A61—MEDICAL OR VETERINARY SCIENCE; HYGIENE
- A61B—DIAGNOSIS; SURGERY; IDENTIFICATION
- A61B5/00—Measuring for diagnostic purposes; Identification of persons
- A61B5/05—Detecting, measuring or recording for diagnosis by means of electric currents or magnetic fields; Measuring using microwaves or radio waves
- A61B5/0507—Detecting, measuring or recording for diagnosis by means of electric currents or magnetic fields; Measuring using microwaves or radio waves using microwaves or terahertz waves
-
- A—HUMAN NECESSITIES
- A61—MEDICAL OR VETERINARY SCIENCE; HYGIENE
- A61B—DIAGNOSIS; SURGERY; IDENTIFICATION
- A61B5/00—Measuring for diagnostic purposes; Identification of persons
- A61B5/72—Signal processing specially adapted for physiological signals or for diagnostic purposes
- A61B5/7203—Signal processing specially adapted for physiological signals or for diagnostic purposes for noise prevention, reduction or removal
-
- A—HUMAN NECESSITIES
- A61—MEDICAL OR VETERINARY SCIENCE; HYGIENE
- A61B—DIAGNOSIS; SURGERY; IDENTIFICATION
- A61B5/00—Measuring for diagnostic purposes; Identification of persons
- A61B5/72—Signal processing specially adapted for physiological signals or for diagnostic purposes
- A61B5/7225—Details of analog processing, e.g. isolation amplifier, gain or sensitivity adjustment, filtering, baseline or drift compensation
-
- A—HUMAN NECESSITIES
- A61—MEDICAL OR VETERINARY SCIENCE; HYGIENE
- A61B—DIAGNOSIS; SURGERY; IDENTIFICATION
- A61B5/00—Measuring for diagnostic purposes; Identification of persons
- A61B5/72—Signal processing specially adapted for physiological signals or for diagnostic purposes
- A61B5/7235—Details of waveform analysis
-
- A—HUMAN NECESSITIES
- A61—MEDICAL OR VETERINARY SCIENCE; HYGIENE
- A61B—DIAGNOSIS; SURGERY; IDENTIFICATION
- A61B5/00—Measuring for diagnostic purposes; Identification of persons
- A61B5/72—Signal processing specially adapted for physiological signals or for diagnostic purposes
- A61B5/7235—Details of waveform analysis
- A61B5/725—Details of waveform analysis using specific filters therefor, e.g. Kalman or adaptive filters
-
- A—HUMAN NECESSITIES
- A61—MEDICAL OR VETERINARY SCIENCE; HYGIENE
- A61B—DIAGNOSIS; SURGERY; IDENTIFICATION
- A61B5/00—Measuring for diagnostic purposes; Identification of persons
- A61B5/72—Signal processing specially adapted for physiological signals or for diagnostic purposes
- A61B5/7235—Details of waveform analysis
- A61B5/7253—Details of waveform analysis characterised by using transforms
Landscapes
- Health & Medical Sciences (AREA)
- Life Sciences & Earth Sciences (AREA)
- Engineering & Computer Science (AREA)
- General Health & Medical Sciences (AREA)
- Medical Informatics (AREA)
- Veterinary Medicine (AREA)
- Public Health (AREA)
- Animal Behavior & Ethology (AREA)
- Physics & Mathematics (AREA)
- Signal Processing (AREA)
- Biophysics (AREA)
- Pathology (AREA)
- Biomedical Technology (AREA)
- Heart & Thoracic Surgery (AREA)
- Surgery (AREA)
- Molecular Biology (AREA)
- Physiology (AREA)
- Artificial Intelligence (AREA)
- Psychiatry (AREA)
- Computer Vision & Pattern Recognition (AREA)
- Cardiology (AREA)
- Power Engineering (AREA)
- Nuclear Medicine, Radiotherapy & Molecular Imaging (AREA)
- Radiology & Medical Imaging (AREA)
- Measuring Pulse, Heart Rate, Blood Pressure Or Blood Flow (AREA)
Abstract
本发明提出了一种基于77GHz毫米波雷达的非接触式精确心跳检测方法。首先，采用连续波多普勒雷达体征检测的原理捕捉人体心脏跳动引起的微弱胸腔位移引起的频移，即雷达向受试者通过发射端发射T(t)。随后T(t)信号的相位被受试者胸壁运动x(t)所调制，反射信号被雷达接收端所捕获；其次，反射信号经过前端的低噪声放大器后，得到基带信号，采用正交混频器解调基带信号，再经过增益放大器后，由反正切解调和接卷绕计算可得到胸壁信号；最后，胸壁模拟信号经过模数转换器得到数字信号x(t)，在数字信号处理后在PC端经过心跳信号提取算法VMD和心跳频谱估计算法CZT后得到心跳频谱结果。
Description
技术领域
本发明涉及心率检测领域，特别是涉及一种基于77GHz毫米雷达的非接触式精确心跳检测方法。
背景技术
随着生活质量的提高，人们对生命健康越来越重视，而心率检测对于日常健康监测和临床医疗诊断有着重大作用。一方面，心率能反映基本的健康状态；另一方面，心跳的节律和强度在很大程度上能体现内脏器官的病理变化。相比传统的接触式心跳监测，如ECG、PPG、腕带式血氧仪等，基于毫米波雷达的非接触式心跳检测可实现远程监护，在某些特殊场合有着显著优势，如灾后生命检测、婴儿监护、疲劳监测等。
毫米波雷达具有体积小、空间分辨率高、大带宽等优点，毫米波波长较短，对胸壁位移这类微小运动具有更高灵敏度。Obeid等人通过2.4GHz、5.8GHz和60GHz三种多普勒雷达对比，验证了一个事实：雷达工作频率越高，对小位移检测灵敏度越高。因此，77GHz毫米波雷达已被广泛应用与车载雷达、人体姿态识别、体征检测等领域。
雷达心跳检测的挑战在于呼吸及其谐波强干扰以及频谱分辨率不足。前者由于呼吸引起的胸壁位移(1～12mm)远大于心跳引起的胸壁位移(0.1～0.5mm)，导致心跳信号被淹没。Hu等人提出连续小波变换与经验模态分解相结合的方法，实现了心肺信号分离与提取，但该方法对信噪要求较高，需预先基带信号校准。随后，变分模态分解算法被提出，在变分框架中自适应分解信号，解决了EMD模态混叠问题。后者，传统心率估计采用超过10s的时间窗口来提高频率分辨率。但时间窗口较长，丢失的心跳细节越多，实时性越差。对于短时窗下的心率估计，将频率时间相位回归技术用于估计心率，可以显著提高了频率分辨率，但短时窗下会出现严重主频估计偏差。然而，实时心率估计往往只需要对心率范围内的窄带频谱区间进行细微观察与分析，频谱细化算法能对信号频谱中感兴趣的频带进行局部放大，有利于心率的准确估计。
发明内容
本发明提出一种基于77GHz毫米雷达的非接触式精确心率检测的方法，用于估计实时心率。
一种基于77GHz毫米雷达的非接触式精确心率检测的方法，其步骤如下：
步骤一：采用77GHz毫米波雷达发射雷达信号，获取回波信号；
步骤二：通过前端的低噪声放大器和正交混频器得到基带信号；
步骤三：经过增益放大器增强和反正切解调以及解卷绕技术得到胸壁位移信号；
步骤四：胸壁信号经过模数转换器得到数字信号；
步骤五：得到的数字信号选择巴特沃斯滤波器进行雷达信号预处理；
步骤六：设计自适应变分模态分解方法(AVMD)得到本征模函数分量(IMF)进行心跳信号提取；
步骤七：设计并采用线性调频Z变换算法(CZT)估计心跳频谱；
所述获取回波信号的方法是：通过IWR1443芯片上的雷达模块发射端发射频率为77GHz的毫米波T(t)到目标特征体处，雷达信号可表示为
由频域中相移的变化来捕捉由呼吸和心跳引起的胸壁位移信号，反射信号R(t)被雷达接收端所捕获，反射信号可表示为：
其中c为无线电的波速，λ＝c/f表示雷达波长，
我们可以把R(t)近似看作经过2d0/c时间延迟并相位调制后的T(t)，反射信号R(t)再经过前端的低噪声放大器后，得到基带信号B(t)，
其中，
为常数相位偏移，与接收机本身参数以及d0相关；残余相位噪声简化为
然后，由正交混频器解调基带信号B(t)，同相和正交信号分别表示为：
再经过增益放大器后，由呼吸和心跳引起的胸壁信号x(t)由反正切解调和解卷绕技术计算可得，表达式如下：
对得到的雷达所捕获的胸壁位移信号进行预处理，预处理步骤包括去除背景噪声和多项式趋势、带通滤波，由于雷达测量范围内除受试者胸壁运动外还存在大量静止物体，通过减去原始雷达回波数据的平均值来去除静止背景噪声，另外，由于热噪声和时间漂移引起的发射单元幅值的不稳定导致回波数据呈现严重的多项式趋势，通过减去最小二乘法拟合的曲线来去除多项式趋势，此外，基于呼吸频率典型范围为0.1～0.3Hz、HR范围为1～3Hz这一事实，巴特沃斯滤波器，其带通范围设置为0.8～4Hz，并应用于去噪后的胸壁位移信号x(t)，就得到了预处理后的胸壁位移信号
基于变分模态分解(VMD)设计自适应变分模态分解算法(AVMD)：
VMD算法中本征模函数表达式为
其中，Ak(t)、
分别为瞬时幅值和相位，则uk(t)的瞬时频率
在VMD算法中，为获取IMF分量，VMD摒弃EMD的循环筛选的方式，而在变分框架下分解，在变分模型的迭代求解过程中，更新各个IMD分量的中心频率和带宽，根据信号频域特征自适应分解信号频带，最终得到窄带IMF分量，假设将原始信号
分解为K个IMF分量，相应的约束变分模型可表示为：
其中，uk代表VMD分解的K个IMF分量；wk表示IMF分量的中心频率；
为求解上述约束变分问题最优解，引入二次惩罚因子α和拉格朗日算法λ构造增广拉格朗日函数：
最后利用交替方向乘子算法求解上述增广拉格朗日函数的鞍点，得到约束模型最优解，然后得到
的分解结果；
分解的模态数K和二次惩罚因子α的选择控制原始信号和分解信号的重构误差，参数的经验选择可能导致较高的重构误差，甚至丢失原始信号的某些重要信息，K过小，会出现模态混叠的现象；K过大，会出现过分解，产生虚假分量，α越小，IMF分量带宽越大，越容易模态混叠；α越大，IMF分量带宽越小，避免模态混叠，但计算量增大，实际应用中α＞fs/2,fs为采样频率，一般可取α＝fs；
自适应变分模态分解避免了VMD分解参数经验选择导致严重重构误差；AVMD算法提取心跳过程如下：对
进行VMD分解并计算每个IMF分量与
的相关系数，设置相关系数阈值δ，若相关系数最小值小于δ则停止分解，得到IMF分量，其中相关系数最大的即为心跳信号xh(t)；
采用CZT算法对心跳频谱估计的具体流程是：
CZT本质是计算沿Z平面上一段螺旋线周线等角度间隔采样的有限采样点的Z变换值，其定义可表示为：
其中，x(n)为采用序列，n＝1，2，...，N-1；A0和θ0分别表示起始采样点的矢量半径和相角；W0表示螺旋线伸展趋势；
表示相邻采样点间的夹角；M表示目标频带的采样点数，m＝0，1，...，M-1；
线性调频Z变换算法(CZT)首先根据静坐时正常心率范围确定目标频带(fmin，fmax)与采样频率fz，然后利用公式:
计算所需要的频谱分辨率Δf以及M、θ0、
令A0＝W0＝1，代入相应参数至公式(4a)中可得目标频带的频谱：
在上述目标频谱X(z)中寻找最大频谱峰值，谱峰对应的频率即为当前心率HR。
本发明的优点在于：
1)使用77HGz毫米波雷达，提高检测灵敏度。
2)AVMD算法实现了在呼吸及其谐波的强干扰下对心跳信号的精确提取。
3)将频谱细化算法引入心跳频谱估计，该方法实现了快速、精确的心率估计，并解决了短时窗下频谱分辨率不足的难题。
4)本发明的算法在短时间内也可以精确且快速地提取心率特征。
附图说明
图1是本发明77GHz毫米波雷达心跳检测系统模型图。
图2是本发明CZT算法与传统FFT算法心跳信号的局部频谱图。
图3(a)是本发明在T＝10s时CZT算法与ECG算法和FTPR算法心率估计曲线实验比较结果。
图3(b)是本发明在T＝3s时CZT算法与ECG算法和FTPR算法心率估计曲线实验比较结果。
具体实施方式
下面结合附图对本发明作进一步详细的阐述。
本发明提出了一种基于77GHz毫米波雷达的非接触式精确心跳检测方法。首先，采用连续波多普勒雷达体征检测的原理捕捉人体心脏跳动引起的微弱胸腔位移引起的频移，即雷达向受试者通过发射端发射T(t)。随后T(t)信号的相位被受试者胸壁运动x(t)所调制，反射信号被雷达接收端所捕获；其次，反射信号经过前端的低噪声放大器后，得到基带信号，采用正交混频器解调基带信号，再经过增益放大器后，由反正切解调和接卷绕计算可得到胸壁信号；最后，胸壁模拟信号经过模数转换器得到数字信号x(t)，在数字信号处理后在PC端经过心跳信号提取算法VMD和心跳频谱估计算法CZT后得到心跳频谱结果。
如图1所示：一种基于AVMD和CZT组合的77GHz毫米雷达的非接触式精确心率检测的方法，其步骤如下：
步骤一：采用77GHz毫米波雷达发射雷达信号，获取回波信号；
步骤二：通过前端的低噪声放大器和正交混频器得到基带信号；
步骤三：经过增益放大器增强和反正切解调以及解卷绕技术得到胸壁位移信号；
步骤四：胸壁信号经过模数转换器得到数字信号；
步骤五：得到的数字信号选择巴特沃斯滤波器进行雷达信号预处理；
步骤六：设计自适应变分模态分解方法(AVMD)得到本征模函数分量(IMF)进行心跳信号提取；
步骤七：设计并采用线性调频Z变换算法(CZT)估计心跳频谱。
所述的采用相位采样获取回波信号的方法是：作为优选，本发明主控芯片采用IWR1443芯片，运算速度快，拥有丰富的外设结构。通过IWR1443芯片上的雷达模块发射端发射频率为77GHz的毫米波T(t)到目标特征体处(雷达信号可表示为
)中，f为雷达工作频率，
为初始相位。由频域中相移的变化来捕捉由呼吸和心跳引起的胸壁位移信号，反射信号R(t)被雷达接收端所捕获,反射信号可表示为：
其中c为无线电的波速，λ＝c/f表示雷达波长。
我们可以把R(t)近似看作经过2d0/c时间延迟并相位调制后的T(t)。反射信号R(t)再经过前端的低噪声放大器后，得到基带信号B(t)。
其中，
为常数相位偏移，与接收机本身参数以及d0相关；残余相位噪声简化为
然后，由正交混频器(I/Q Mixer)解调基带信号B(t)，同相和正交信号分别表示为：
再经过增益放大器(Gain)后，由呼吸和心跳引起的胸壁信号x(t)由反正切解调和解卷绕技术计算可得，表达式如下：
最后，模拟信号x(t)经过模数转换器(ADC)得到数学信号x(t)，用于数字信号处理(DSP)以及PC端心跳检测算法的开发。
采用心跳检测算法获取心跳频谱的方法是：为实现复杂环境下多普勒雷达对心跳信号的精确估计，提出了一种雷达体征数据处理的具体实现方案，主要包括三部分：雷达信号预处理、心跳信号提取、心跳频谱估计。雷达捕获的受试者胸壁位移信号x(t)作为本发明算法的输入。
预处理步骤包括去除背景噪声和多项式趋势、带通滤波。由于雷达测量范围内除受试者胸壁运动外还存在大量静止物体，通过减去原始雷达回波数据的平均值来去除静止背景噪声。另外，由于热噪声和时间漂移引起的发射单元幅值的不稳定导致回波数据呈现严重的多项式趋势，通过减去最小二乘法拟合的曲线来去除多项式趋势。
此外，基于呼吸频率典型范围为0.1～0.3Hz、HR范围为1～3Hz这一事实，带通滤波器选择巴特沃斯滤波器，其通带范围设置为0.8～4Hz，并应用于去噪后的x(t)，得到预处理后的信号
基于变分模态分解(VMD)设计自适应变分模态分解算法(AVMD)：
VMD算法中本征模函数表达式为
其中，Ak(t)、
分别为瞬时幅值和相位，则uk(t)的瞬时频率
在VMD算法中，为获取IMF分量，VMD摒弃EMD的循环筛选的方式，而在变分框架下分解。在变分模型的迭代求解过程中，更新各个IMD分量的中心频率和带宽，根据信号频域特征自适应分解信号频带，最终得到窄带IMF分量。假设将原始信号
分解为K个IMF分量，相应的约束变分模型可表示为：
其中，uk代表VMD分解的K个IMF分量；wk表示IMF分量的中心频率。
为求解上述约束变分问题最优解，引入二次惩罚因子α和拉格朗日算法λ构造增广拉格朗日函数：
最后利用交替方向乘子算法求解上述增广拉格朗日函数的鞍点，得到约束模型最优解，然后得到
的分解结果。
分解的模态数K和二次惩罚因子α的选择控制原始信号和分解信号的重构误差，参数的经验选择可能导致较高的重构误差，甚至丢失原始信号的某些重要信息。K过小，会出现模态混叠的现象；K过大，会出现过分解，产生虚假分量。α越小，IMF分量带宽越大，越容易模态混叠；α越大，IMF分量带宽越小，避免模态混叠，但计算量增大。实际应用中α＞fs/2，fs为采样频率，一般可取α＝fs。
自适应变分模态分解避免了VMD分解参数经验选择导致严重重构误差。AVMD算法提取心跳过程如下：首先根据采样频率设定二次惩罚因子α，随后对
进行VMD分解并计算每个IMF分量与
的相关系数，设置相关系数阈值δ，若相关系数最小值小于δ则停止分解，得到IMF分量，其中相关系数最大的即为心跳信号xh(t)。
表1的算法1总结了AVMD算法提取心跳信号的处理流程。
表1
采用CZT算法对心跳频谱估计的具体流程是；
线性调频Z变换相比传统的FFT只能得到粗略的“全体频谱”，CZT能实现信号的窄带分析，显著提高频率估计准确性。CZT本质是计算沿Z平面上一段螺旋线周线等角度间隔采样的有限采样点的Z变换值，其定义可表示为：
其中，x(n)为采用序列，n＝1，2，...，N-1；A0和θ0分别表示起始采样点的矢量半径和相角；W0表示螺旋线伸展趋势；
表示相邻采样点间的夹角；M表示目标频带的采样点数，m＝0，1，...，M-1；
线性调频Z变换算法(CZT)首先根据静坐时正常心率范围确定目标频带(fmin，fmax)与采样频率fz，然后利用公式:
计算所需要的频谱分辨率Δf以及M、θ0、
令A0＝W0＝1，代入相应参数至公式(4a)中可得目标频带的频谱：
在上述目标频谱X(z)中寻找最大频谱峰值，谱峰对应的频率即为当前心率HR。表2是算法2详细阐述CZT算法心跳频谱估计流程。
表2
为分析CZT算法心率估计的准确性，与传统的FFT进行了对比分析。以雷达信号处理过程得到的心跳信号xh(t)为例，ECG测量的心率作为参考，信号采样频率fz＝32Hz,时间窗口t＝3s，静坐心率频带为1～2Hz，其频谱计算结果如图2所示，FFT和CZT估计的心率分别为1.333Hz和1.427Hz，参考心率为1.418Hz。对比心率估计误差，FFT是CZT的10倍左右。此外，表3列举了四组不同时间段下两种算法的心率估计误差，CZT心率估计误差明显低于FFT，验证了CZT频谱细化算法心率估计的准确性。
表3
实施例
本发明体征检测实验在真实的办公环境下实施，实验的参考心跳数据由传统的接触式ECG模块测量，ECG电极贴于受试者心脏附近，采样频率设为500Hz。另外，实验采用德州仪器公司设计的IWR1443Boost毫米波雷达平台来测量体征信息，其相关实验参数如表4所示。实验过程中，受试者静坐与雷达正前方0.5m处，胸腔与雷达基本处于同一水平线，保持正常呼吸，记录两分钟实验数据。
此外，在实际应用中，需要选择合适的时间分辨率与频率分辨率。一般而言，时间窗口越短，时间分辨率越好，呈现更多心跳细节，但同时频率分辨率越低，检测精度越差。作为优选，本发明选择3s和10s两种不同时间窗口长度用于实时心率检测。
表4
性能分析：
为验证本发明提出的心跳检测算法的有效性与鲁棒性，我们在不同的实验条件下测量了多组实验数据，并采用多个评价指标对本文方法的性能进行评估。
此外，在实际应用中，需要选择合适的时间分辨率与频率分辨率。一般而言，时间窗口越短，时间分辨率越好，呈现更多心跳细节，但同时频率分辨率越低，检测精度越差。本发明基于77GHz毫米波雷达非接触式检测心跳方法具有普适性。
心率(HR)和心率变异性(HRV)是人们在心脏跳动方面衡量身体健康的两个重要指标。
心率(HR)是指每分钟的平均心跳次数，常被用做健康检测的基本指标。心率准确度(HRA)定义为多普勒雷达检测到的心率与ECG检测到的参考心率偏差在±2％以内的时间百分比，我们也将其作为心率指标。
心率变异性(HRV)是指连续心跳之间的时间间隔，揭示了每次心跳的规律性。心率变异性分析对心脏疾病的诊断和治疗具有重要意义。不同于在时域中通过检测相邻心跳的时间间隔来确定心率，本发明采用在频域中搜索心跳频率的心率检测方法，其中连续心跳时间间隔BBIs可用BPM来表示：
另外，平均相对误差(MRE)能被表示为：
HRV分析评价指标：RR间期标准差(SDNN)、连续RR间期的均方根差(RMSSD)，其中，NBBI为连续心跳间隔数，可分别表示为：
表5(a)和表5(b)为CZT算法与FTPR算法的HRV分析结果。如表5(a)所示，CZT算法的MRE误差范围为0.41％～0.60％，表明雷达与ECG间一致性很好。另外，所有的RMSSD都小于3ms，这表明长时窗时心率细节严重丢失，与HR分析结果一致。如表5(b)所示，FTPR、CZT算法的平均MRE分别为2.34％、1.01％。两种算法与ECG间HRV评价指标差值的平均值分别为SDNN(4.60ms、1.68ms)、RMSSD(4.93ms、1.53ms)。对比上述误差，CZT检测性能明显优于FTPR，验证了短时窗下本发明算法的优越性。
表5(a)
表5(b)
为验证本发明CZT算法的可行性，在相同检测环境标准下与ECG算法和FTPR算法进行比较验证。
通过分析3s和10s两个不同时窗下HR和HRV相关指标，来对比FTPR、CZT算法的心率检测性能。
心率(HR)是指每分钟的平均心跳次数，常被用作健康监测的基本指标。心率准确度(HRA)定义为多普勒雷达检测到的心率与ECG检测到的参考心率偏差在±2％以内的时间百分比。
图3(a)和图3(b)展示了不同时窗下不同算法的HR估计曲线，ECG作为参考。如图3(a)所示，在±2％误差范围内，FTPR、CZT算法的HR准确率高达95.5％、97.3％。然而，当心率出现较大突变时FTPR心率曲线中存在严重的偏离点。如图3(b)所示，由于3s时窗较短，频谱分辨率较低，FTPR算法的HR准确率仅为45.76％，但CZT算法的HR准确率仍超过了90％，且呈现了更多的心跳细节。此外，当短时心率发生较大突变时CZT算法估计的心率值与参考心率一致性仍较好，进一步验证了CZT算法的有效性与鲁棒性。
Claims (1)
1.基于77GHz毫米雷达的非接触式精确心率检测的方法，其特征在于，其步骤如下：
步骤一：采用77GHz毫米波雷达发射雷达信号，获取回波信号；
步骤二：通过前端的低噪声放大器和正交混频器得到基带信号；
步骤三：经过增益放大器增强和反正切解调以及解卷绕技术得到胸壁位移信号；
步骤四：胸壁信号经过模数转换器得到数字信号；
步骤五：得到的数字信号选择巴特沃斯滤波器进行雷达信号预处理；
步骤六：设计自适应变分模态分解方法(AVMD)得到本征模函数分量(IMF)进行心跳信号提取；
步骤七：设计并采用线性调频Z变换算法(CZT)估计心跳频谱；
所述获取回波信号的方法是：通过IWR1443芯片上的雷达模块发射端发射频率为77GHz的毫米波T(t)到目标特征体处，雷达信号可表示为
由频域中相移的变化来捕捉由呼吸和心跳引起的胸壁位移信号，反射信号R(t)被雷达接收端所捕获，反射信号可表示为：
其中c为无线电的波速，λ＝c/f表示雷达波长，
我们可以把R(t)近似看作经过2d0/c时间延迟并相位调制后的T(t)，反射信号R(t)再经过前端的低噪声放大器后，得到基带信号B(t)，
其中，
为常数相位偏移，与接收机本身参数以及d0相关；残余相位噪声简化为
然后，由正交混频器解调基带信号B(t)，同相和正交信号分别表示为：
再经过增益放大器后，由呼吸和心跳引起的胸壁信号x(t)由反正切解调和解卷绕技术计算可得，表达式如下：
对得到的雷达所捕获的胸壁位移信号进行预处理，预处理步骤包括去除背景噪声和多项式趋势、带通滤波，由于雷达测量范围内除受试者胸壁运动外还存在大量静止物体，通过减去原始雷达回波数据的平均值来去除静止背景噪声，另外，由于热噪声和时间漂移引起的发射单元幅值的不稳定导致回波数据呈现严重的多项式趋势，通过减去最小二乘法拟合的曲线来去除多项式趋势，此外，基于呼吸频率典型范围为0.1～0.3Hz、HR范围为1～3Hz这一事实，巴特沃斯滤波器，其带通范围设置为0.8～4Hz，并应用于去噪后的胸壁位移信号x(t)，就得到了预处理后的胸壁位移信号
基于变分模态分解(VMD)设计自适应变分模态分解算法(AVMD)：
VMD算法中本征模函数表达式为
其中，Ak(t)、
分别为瞬时幅值和相位，则uk(t)的瞬时频率
在VMD算法中，为获取IMF分量，VMD摒弃EMD的循环筛选的方式，而在变分框架下分解，在变分模型的迭代求解过程中，更新各个IMD分量的中心频率和带宽，根据信号频域特征自适应分解信号频带，最终得到窄带IMF分量，假设将原始信号
分解为K个IMF分量，相应的约束变分模型可表示为：
其中，uk代表VMD分解的K个IMF分量；wk表示IMF分量的中心频率；
为求解上述约束变分问题最优解，引入二次惩罚因子α和拉格朗日算法λ构造增广拉格朗日函数：
最后利用交替方向乘子算法求解上述增广拉格朗日函数的鞍点，得到约束模型最优解，然后得到
的分解结果；
分解的模态数K和二次惩罚因子α的选择控制原始信号和分解信号的重构误差，参数的经验选择可能导致较高的重构误差，甚至丢失原始信号的某些重要信息，K过小，会出现模态混叠的现象；K过大，会出现过分解，产生虚假分量，α越小，IMF分量带宽越大，越容易模态混叠；α越大，IMF分量带宽越小，避免模态混叠，但计算量增大，实际应用中α＞fs/2，fs为采样频率，一般可取α＝fs；
自适应变分模态分解避免了VMD分解参数经验选择导致严重重构误差；AVMD算法提取心跳过程如下：对
进行VMD分解并计算每个IMF分量与
的相关系数，设置相关系数阈值δ，若相关系数最小值小于δ则停止分解，得到IMF分量，其中相关系数最大的即为心跳信号xh(t)；
采用CZT算法对心跳频谱估计的具体流程是：
CZT本质是计算沿Z平面上一段螺旋线周线等角度间隔采样的有限采样点的Z变换值，其定义可表示为：
其中，x(n)为采用序列，n＝1，2，...，N-1；A0和θ0分别表示起始采样点的矢量半径和相角；W0表示螺旋线伸展趋势；
表示相邻采样点间的夹角；M表示目标频带的采样点数，m＝0，1，...，M-1；
线性调频Z变换算法(CZT)首先根据静坐时正常心率范围确定目标频带(fmin，fmax)与采样频率fz，然后利用公式：
计算所需要的频谱分辨率Δf以及M、θ0、
令A0＝W0＝1，代入相应参数至公式(4a)中可得目标频带的频谱：
在上述目标频谱X(z)中寻找最大频谱峰值，谱峰对应的频率即为当前心率HR。
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202111190201.9A CN113854992A (zh) | 2021-10-13 | 2021-10-13 | 基于77GHz毫米雷达的非接触式精确心率检测方法 |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202111190201.9A CN113854992A (zh) | 2021-10-13 | 2021-10-13 | 基于77GHz毫米雷达的非接触式精确心率检测方法 |
Publications (1)
| Publication Number | Publication Date |
|---|---|
| CN113854992A true CN113854992A (zh) | 2021-12-31 |
Family
ID=78999376
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN202111190201.9A Pending CN113854992A (zh) | 2021-10-13 | 2021-10-13 | 基于77GHz毫米雷达的非接触式精确心率检测方法 |
Country Status (1)
| Country | Link |
|---|---|
| CN (1) | CN113854992A (zh) |
Cited By (5)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN114403820A (zh) * | 2022-01-13 | 2022-04-29 | 武汉烽火凯卓科技有限公司 | 一种移动目标的生命体征检测方法及系统 |
| CN115089143A (zh) * | 2022-08-24 | 2022-09-23 | 江苏亿连通信技术有限公司 | 毫米波雷达生命体征信号提取和测量方法 |
| CN115312195A (zh) * | 2022-10-10 | 2022-11-08 | 安徽交欣科技股份有限公司 | 一种基于情绪数据计算个体心理异常的健康评估方法 |
| CN115969339A (zh) * | 2023-02-13 | 2023-04-18 | 山东师范大学 | 非接触的心率监测方法、系统、存储介质及设备 |
| CN117579136A (zh) * | 2024-01-17 | 2024-02-20 | 南京控维通信科技有限公司 | Tdma中网控系统对反向突发的aupc及acm控制方法 |
Citations (5)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN108614259A (zh) * | 2018-05-02 | 2018-10-02 | 电子科技大学 | 一种基于超宽带雷达传感器的心跳呼吸特征监测方法 |
| CN112303504A (zh) * | 2020-11-09 | 2021-02-02 | 吉林大学 | 一种基于改进的变分模式分解算法的供水管道泄漏位置检测方法 |
| CN112617773A (zh) * | 2019-09-24 | 2021-04-09 | 新加坡国立大学 | 一种用于健康监测的信号处理方法及信号处理装置 |
| CN112754441A (zh) * | 2021-01-08 | 2021-05-07 | 杭州环木信息科技有限责任公司 | 一种基于毫米波的非接触式心跳检测方法 |
| CN113273978A (zh) * | 2021-05-21 | 2021-08-20 | 电子科技大学 | 一种基于超宽带雷达的人体呼吸和心跳频率的检测方法 |
-
2021
- 2021-10-13 CN CN202111190201.9A patent/CN113854992A/zh active Pending
Patent Citations (5)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN108614259A (zh) * | 2018-05-02 | 2018-10-02 | 电子科技大学 | 一种基于超宽带雷达传感器的心跳呼吸特征监测方法 |
| CN112617773A (zh) * | 2019-09-24 | 2021-04-09 | 新加坡国立大学 | 一种用于健康监测的信号处理方法及信号处理装置 |
| CN112303504A (zh) * | 2020-11-09 | 2021-02-02 | 吉林大学 | 一种基于改进的变分模式分解算法的供水管道泄漏位置检测方法 |
| CN112754441A (zh) * | 2021-01-08 | 2021-05-07 | 杭州环木信息科技有限责任公司 | 一种基于毫米波的非接触式心跳检测方法 |
| CN113273978A (zh) * | 2021-05-21 | 2021-08-20 | 电子科技大学 | 一种基于超宽带雷达的人体呼吸和心跳频率的检测方法 |
Non-Patent Citations (4)
| Title |
|---|
| 孙曙光等: "基于变分模态分解与Teager能量算子的谐波/间谐波检测方法", 电测与仪表, vol. 57, no. 02, 30 July 2019 (2019-07-30), pages 109 - 115 * |
| 孙曙光等: "基于变分模态分解与Teager能量算子的谐波/间谐波检测方法", 电测与仪表, vol. 57, no. 2, 25 January 2020 (2020-01-25), pages 109 - 115 * |
| 蒋腾等: "基于连续小波变换的多数据心率提取方法", 现代雷达, vol. 41, no. 05, 15 May 2019 (2019-05-15), pages 22 - 26 * |
| 马可等: "CZT和ZFFT频谱细化性能分析及FPGA实现", 计算机测量与控制, vol. 24, no. 02, 25 February 2016 (2016-02-25), pages 288 - 298 * |
Cited By (6)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN114403820A (zh) * | 2022-01-13 | 2022-04-29 | 武汉烽火凯卓科技有限公司 | 一种移动目标的生命体征检测方法及系统 |
| CN115089143A (zh) * | 2022-08-24 | 2022-09-23 | 江苏亿连通信技术有限公司 | 毫米波雷达生命体征信号提取和测量方法 |
| CN115312195A (zh) * | 2022-10-10 | 2022-11-08 | 安徽交欣科技股份有限公司 | 一种基于情绪数据计算个体心理异常的健康评估方法 |
| CN115969339A (zh) * | 2023-02-13 | 2023-04-18 | 山东师范大学 | 非接触的心率监测方法、系统、存储介质及设备 |
| CN117579136A (zh) * | 2024-01-17 | 2024-02-20 | 南京控维通信科技有限公司 | Tdma中网控系统对反向突发的aupc及acm控制方法 |
| CN117579136B (zh) * | 2024-01-17 | 2024-04-02 | 南京控维通信科技有限公司 | Tdma中网控系统对反向突发的aupc及acm控制方法 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| CN113854992A (zh) | 基于77GHz毫米雷达的非接触式精确心率检测方法 | |
| EP2757944B1 (en) | Systems and methods for determining respiration information from a photoplethysmograph | |
| EP2757943B1 (en) | Systems and methods for determining respiration information from a photoplethysmograph | |
| US20150150515A1 (en) | Respiration rate extraction from cardiac signals | |
| CN113951856A (zh) | 一种基于多普勒雷达心跳检测的频谱估计方法 | |
| CN114305364B (zh) | 基于毫米波雷达的血压检测方法、系统及设备 | |
| Han et al. | UWB radar for non-contact heart rate variability monitoring and mental state classification | |
| CN113273978B (zh) | 一种基于超宽带雷达的人体呼吸和心跳频率的检测方法 | |
| TW201225912A (en) | Method for measuring physiological parameters | |
| CN115363547B (zh) | 一种基于超宽带雷达相参积累的人体生命体征检测方法 | |
| CN108814579A (zh) | 一种基于emd分解的心电、呼吸联合计算心率变异性的方法 | |
| TWI505816B (zh) | 血氧飽和度檢測方法及裝置 | |
| Aardal et al. | Radar cross section of the human heartbeat and respiration in the 500MHz to 3GHz band | |
| CN111685760B (zh) | 一种基于雷达测量的人体呼吸频率计算方法 | |
| Pan et al. | A spectrum estimation approach for accurate heartbeat detection using Doppler radar based on combination of FTPR and TWV | |
| Will et al. | Intelligent signal processing routine for instantaneous heart rate detection using a Six-Port microwave interferometer | |
| Zhao et al. | Multi-target vital signs remote monitoring using mmWave FMCW radar | |
| CN115813363A (zh) | 一种基于心跳二次谐波重构的心率估计方法 | |
| JP7170147B2 (ja) | 非一時的なコンピュータ可読記憶媒体、コンピュータ実装方法及びシステム | |
| US20220323023A1 (en) | Method for determining respiratory rate | |
| Sameera et al. | Effect of respiration harmonics on beat-to-beat analysis of heart signal | |
| CN112043256A (zh) | 一种基于雷达的多目标心率实时测量方法 | |
| Zamrath et al. | Robust and computationally efficient approach for Heart Rate monitoring using photoplethysmographic signals during intensive physical excercise | |
| CN115590489B (zh) | 一种基于调频连续波雷达的无接触血压监测方法 | |
| Chen et al. | A developed algorithm for oscillometric blood pressure measurement |
Legal Events
| Date | Code | Title | Description |
|---|---|---|---|
| PB01 | Publication | ||
| PB01 | Publication | ||
| SE01 | Entry into force of request for substantive examination | ||
| SE01 | Entry into force of request for substantive examination |