CN114732390A - 一种基于调频连续波雷达的非接触式心率变异性监测方法 - Google Patents

一种基于调频连续波雷达的非接触式心率变异性监测方法 Download PDF

Info

Publication number
CN114732390A
CN114732390A CN202210539306.9A CN202210539306A CN114732390A CN 114732390 A CN114732390 A CN 114732390A CN 202210539306 A CN202210539306 A CN 202210539306A CN 114732390 A CN114732390 A CN 114732390A
Authority
CN
China
Prior art keywords
signal
acceleration
sequence
heartbeat
frequency
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
Application number
CN202210539306.9A
Other languages
English (en)
Inventor
刘涛
韩翔宇
翟潜
韩梅梅
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
Zhejiang University ZJU
Original Assignee
Zhejiang University ZJU
Priority date (The priority date 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 date listed.)
Filing date
Publication date
Application filed by Zhejiang University ZJU filed Critical Zhejiang University ZJU
Priority to CN202210539306.9A priority Critical patent/CN114732390A/zh
Publication of CN114732390A publication Critical patent/CN114732390A/zh
Pending legal-status Critical Current

Links

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
    • A61B5/02405—Determining heart rate variability
    • 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/24—Detecting, measuring or recording bioelectric or biomagnetic signals of the body or parts thereof
    • A61B5/316—Modalities, i.e. specific diagnostic methods
    • A61B5/318—Heart-related electrical modalities, e.g. electrocardiography [ECG]
    • A61B5/346—Analysis of electrocardiograms
    • A61B5/349—Detecting specific parameters of the electrocardiograph cycle
    • 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
    • A61B5/7257—Details of waveform analysis characterised by using transforms using Fourier transforms
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F17/00—Digital computing or data processing equipment or methods, specially adapted for specific functions
    • G06F17/10—Complex mathematical operations
    • G06F17/14—Fourier, Walsh or analogous domain transformations, e.g. Laplace, Hilbert, Karhunen-Loeve, transforms
    • G06F17/141—Discrete Fourier transforms
    • G06F17/142—Fast Fourier transforms, e.g. using a Cooley-Tukey type algorithm
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F17/00—Digital computing or data processing equipment or methods, specially adapted for specific functions
    • G06F17/10—Complex mathematical operations
    • G06F17/18—Complex mathematical operations for evaluating statistical data, e.g. average values, frequency distributions, probability functions, regression analysis
    • Y—GENERAL TAGGING OF NEW TECHNOLOGICAL DEVELOPMENTS; GENERAL TAGGING OF CROSS-SECTIONAL TECHNOLOGIES SPANNING OVER SEVERAL SECTIONS OF THE IPC; TECHNICAL SUBJECTS COVERED BY FORMER USPC CROSS-REFERENCE ART COLLECTIONS [XRACs] AND DIGESTS
    • Y02—TECHNOLOGIES OR APPLICATIONS FOR MITIGATION OR ADAPTATION AGAINST CLIMATE CHANGE
    • Y02A—TECHNOLOGIES FOR ADAPTATION TO CLIMATE CHANGE
    • Y02A90/00—Technologies having an indirect contribution to adaptation to climate change
    • Y02A90/10—Information and communication technologies [ICT] supporting adaptation to climate change, e.g. for weather forecasting or climate simulation

Landscapes

  • Health & Medical Sciences (AREA)
  • Physics & Mathematics (AREA)
  • Engineering & Computer Science (AREA)
  • Life Sciences & Earth Sciences (AREA)
  • Mathematical Physics (AREA)
  • General Physics & Mathematics (AREA)
  • Cardiology (AREA)
  • Data Mining & Analysis (AREA)
  • Veterinary Medicine (AREA)
  • Computational Mathematics (AREA)
  • Biophysics (AREA)
  • Pathology (AREA)
  • Biomedical Technology (AREA)
  • Heart & Thoracic Surgery (AREA)
  • Medical Informatics (AREA)
  • Molecular Biology (AREA)
  • Surgery (AREA)
  • Animal Behavior & Ethology (AREA)
  • General Health & Medical Sciences (AREA)
  • Public Health (AREA)
  • Theoretical Computer Science (AREA)
  • Pure & Applied Mathematics (AREA)
  • Mathematical Optimization (AREA)
  • Mathematical Analysis (AREA)
  • Software Systems (AREA)
  • Physiology (AREA)
  • General Engineering & Computer Science (AREA)
  • Databases & Information Systems (AREA)
  • Algebra (AREA)
  • Signal Processing (AREA)
  • Nuclear Medicine, Radiotherapy & Molecular Imaging (AREA)
  • Artificial Intelligence (AREA)
  • Radiology & Medical Imaging (AREA)
  • Computer Vision & Pattern Recognition (AREA)
  • Discrete Mathematics (AREA)
  • Psychiatry (AREA)
  • Bioinformatics & Cheminformatics (AREA)
  • Bioinformatics & Computational Biology (AREA)
  • Evolutionary Biology (AREA)
  • Operations Research (AREA)

Abstract

本发明公开了一种基于调频连续波雷达的非接触式心率变异性监测方法，具体方法包括向监测对象所在区域发射周期性的线性调频波，并采集雷达回波信号；依次采用混频和FFT对雷达回波信号进行预处理；随后提取监测对象心脏部位所在的位置；采用相位相关的方法，提取监测对象心脏部位的运动信息，并据此计算其加速度信号；对所得加速度信号进行平滑，并通过峰值检测的方法对各次心跳之间分割点的位置进行估计；生成单次心跳的加速度模板信号并对加速度信号进行精密分割，得到各次心跳的心跳间期，计算心率变异性指标。本发明采用微波作为探测媒介，具有抗干扰能力强，对非金属障碍物穿透力强，集成度高，功耗低，测量精度高，鲁棒性强等诸多优点。

Description

一种基于调频连续波雷达的非接触式心率变异性监测方法
技术领域
本发明涉及高新技术改造传统产业与非接触式生命体征监测领域，特别涉及一种具体的基于调频连续波雷达的非接触式心率变异性监测算法的设计。
背景技术
近年来基于非接触式生理特征监测的研究方兴未艾，各种非接触式生理特征监测方法不断涌现。它们的监测媒介多种多样，包括红外，声波和光学等等。然而基于上述媒介的非接触式生理特征监测技术都具有一些明显的缺点。例如：基于红外传感器的监测手段会受到监测对象附近热源的影响；基于超声传感器的监测手段会受到空气湿度、以及棉纱等能够吸收声波的材料的影响；基于光学的传感器会受到光照条件、烟雾以及监测对象身体遮挡等因素的影响。上述缺点限制了这些非接触式生理特征监测方法的推广与应用。
另一方面，心率变异性(heart rate variability，HRV)是人体一项十分重要的生命体征指标，包含着神经体液因素对心血管系统调节的信息，与情绪唤醒度，交感神经张力及平衡性等生理体征有着十分密切的联系。在现有技术中，已经有多种监测手段能够对人体的心跳信号进行监测，从而对心率变异性指标进行分析。然而，心率变异性指标分析对测量精度要求较高；且与呼吸等生理特征相比，心跳信号强度十分微弱，极易受到干扰。因此现有的心率变异性监测手段大多仍基于接触式的心电监测仪或指夹式脉搏波传感器，测量过程中佩戴的指夹，电极等会明显降低使用者的使用体验。
发明内容
为了解决现有非接触式生理特征监测方法测量精度不高，使用条件受限，易受环境因素干扰等方面的问题，同时通过非接触式的方式优化传统心率变异性监测方法的使用体验，本发明提出一种基于调频连续波雷达的非接触式心率变异性监测方法。通过对雷达回波信号处理算法的创新设计，提高回波信号中心跳相关信号的信噪比，最终实现较好的心率变异性监测效果。
本发明所采用的具体技术方案如下：
本发明提供了一种基于调频连续波雷达的非接触式心率变异性监测方法，具体如下：
S1：向监测对象所在区域发射周期性的线性调频波，并使用接收天线采集相应的雷达回波信号；
S2：依次采用混频和快速傅里叶变换的方式对所采集的雷达回波信号进行预处理；
S3：从预处理后的雷达回波信号中提取监测对象心脏部位所在的位置；
S4：采用相位相关的方法，从步骤S3所得位置中提取监测对象心脏部位的运动信息，并据此计算其加速度信号；
S5：通过计算短时平均功率以对所得加速度信号进行平滑，并通过峰值检测的方法对各次心跳之间分割点的位置进行估计；
S6：基于步骤S5中估计得到的各次心跳之间分割点的位置，生成单次心跳的加速度模板信号；随后利用该加速度模板信号对步骤S4中所得加速度信号进行精密分割，得到各次心跳的心跳间期；
S7：基于步骤S6所得的心跳间期，计算心率变异性指标，进而实现非接触式心率变异性监测。
作为优选，所述步骤S1中，线性调频波的脉冲调频范围为77GHz～81GHz，单个脉冲持续周期为50us；快时间轴采样频率为4MHz，采样点数为128；慢时间轴采样频率为100Hz。
作为优选，所述步骤S2中，经混频方式处理后的输出为一个中频信号；该中频信号的幅值为一常数，频率和相位均正比于反射物体相对接收天线阵列的距离。
作为优选，所述步骤S3具体如下：
S3-1：对于接收到的第n个雷达回波信号，经过步骤S2预处理后提取其第k个元素的相位
Figure BDA0003647563420000021
；
S3-2：选取长度为T的时间窗，并依次按照步骤S3-1中的方法处理，得到序列Φk：
Figure BDA0003647563420000031
S3-3：对所得序列Φk解缠绕并通过快速傅里叶变换计算其频谱，统计其在心跳频段[0.8Hz,2Hz]的能量Ek；
S3-4：遍历k以找到Ek的最大值，并据此确定监测对象心脏部位相对于雷达天线的距离，即为监测对象心脏部位所在的位置。
作为优选，所述步骤S4具体如下：
根据公式
Figure BDA0003647563420000032
提取反映监测对象心脏部位与接收天线阵列间距离变化的相位序列
Figure BDA0003647563420000033
解缠绕后获得对应的体动信号序列{Rm}；其中，λ为所发射调频连续脉冲的波长；
随后通过公式
Figure BDA0003647563420000034
求得对应的加速度信号序列{am}；其中，Δt为慢时间轴的采样间隔，Rm+1为体动信号序列{Rm}中的第m+1项，Rm为体动信号序列{Rm}中的第m项，Rm-1为体动信号序列{Rm}的第m-1项。
作为优选，所述步骤S5具体如下：
S5-1：定义在n0时刻长度为L邻域内加速度信号的平均功率PL(n0)为
Figure BDA0003647563420000035
其中，L为计算平均功率选取时间窗长度的一半，ai为加速度信号序列{am}中的第i项；
取L1为单次心跳周期时长一半对应的采样点数，取L2＝0.4L1；按照下式，逐项处理加速度信号序列{am}，以得到平滑后的心跳信号序列{Ym}；
Figure BDA0003647563420000036
其中，L1和L2表示在计算平均功率时选取了不同的时间窗长度，
Figure BDA0003647563420000041
表示在m时刻L1长度邻域内加速度信号的平均功率；
Figure BDA0003647563420000042
表示在m时刻L2长度邻域内加速度信号的平均功率；
S5-2：对于平滑后的心跳信号序列{Ym}，通过峰值检测的方式找出其所有的极大值点在序列中对应的下标t
t＝[t1,t2,…ti…,tI]，
所得ti即是各次心跳之间的分割点在序列中对应的下标；其中I表示找到极大值点的总个数。
作为优选，所述步骤S6具体如下：
S6-1：基于步骤S5中所得各次心跳周期分割点的位置的估计值t，将加速度信号序列{am}进行分割，得到如下(I-1)个数据段：
Figure BDA0003647563420000043
S6-2：对步骤S6-1所得数据段进行平均，得到心跳加速度模板信号Template：
Figure BDA0003647563420000044
计算时通过插值保证S6-1中各个加速度信号的片段
Figure BDA0003647563420000045
长度一致；
S6-3：选取一组新的分割点集S：
S＝[s1,s2,……,sI]，
得到一组对加速度信号序列{am}的新的截断：
Figure BDA0003647563420000046
S6-4：通过如下公式定义步骤S6-3所得新的截断与加速度模板信号的匹配程度：
Figure BDA0003647563420000047
在计算时通过插值保证每个
Figure BDA0003647563420000051
与Template长度相同；
S6-5：重复步骤S6-3～S6-4以计算每组分割点集对应的截断与加速度模板信号的匹配程度，找出D最小的一组分割点集Ss作为当前模板信号下的最优分割；基于Ss按照步骤S6-2以生成新的模板信号；
S6-6：重复S6-5所述过程，直到某次重复过程中通过步骤S6-4计算得到的D小于事先设定的阈值，此时的分割点集S即为最终得到各次心跳之间分割点的位置。
本发明相对于现有技术而言，具有以下有益效果：
本发明采用微波作为探测媒介，可以有效避免传统接触式手段使用方法复杂，使用体验差且使用人群受限等问题，具有抗干扰能力强，对非金属障碍物穿透力强，集成度高，功耗低，测量精度高，鲁棒性强等诸多优点。
附图说明
图1为本发明监测方法的技术流程图；
图2为本发明监测方法的具体步骤流程图；
图3为实施例中各次心跳间期测量效果展示图。
具体实施方式
下面结合附图和具体实施方式对本发明做进一步阐述和说明。本发明中各个实施方式的技术特征在没有相互冲突的前提下，均可进行相应组合。
如图1所示，本发明提供了一种基于调频连续波雷达的非接触式心率变异性监测方法，该方法主要包括以下过程：对雷达回波信号的预处理，监测对象心脏位置的确定，心跳信号的提取，各次心跳之间的分割和心率变异性指标的计算。下面将根据图2对该方法的具体步骤进行说明。
S1：利用发射天线向监测对象所在区域发射周期性的线性调频波，并使用接收天线采集相应的雷达回波信号。
在实际应用时，线性调频波的脉冲调频范围为77GHz～81GHz，单个脉冲持续周期为50us；快时间轴采样频率为4MHz，采样点数为128；慢时间轴采样频率为100Hz。
S2：依次采用混频和快速傅里叶变换(FFT)的方式对所采集的雷达回波信号进行预处理。其中，经混频方式处理后的输出为一个中频信号；该中频信号的幅值为一常数，频率和相位均正比于反射物体相对接收天线阵列的距离。该步骤具体如下：
混频输出为一个中频信号xout，
Figure BDA0003647563420000061
该中频信号的幅值A为常数，频率和相位均正比于反射面和天线的距离R，如下式所示：
Figure BDA0003647563420000062
其中S为线性调频连续波频率随时间的变化率，c为光速，λ为所用调频连续波的波长，均为常数。最后对混频输出的中频信号进行FFT，得到频谱中不同频率的成分即代表不同距离处的物体反射的雷达回波信号。
S3：从预处理后的雷达回波信号中提取监测对象心脏部位所在的位置，步骤具体如下：
S3-1：对于接收到的第n个雷达回波信号，经过步骤S2预处理后提取其第k个元素的相位
Figure BDA0003647563420000063
S3-2：选取长度为T的时间窗，并依次按照步骤S3-1中的方法处理，得到序列Φk：
Figure BDA0003647563420000064
S3-3：对所得序列Φk解缠绕并通过快速傅里叶变换计算其频谱，统计其在心跳频段[0.8Hz,2Hz]的能量Ek；
S3-4：遍历k以找到Ek的最大值，并据此确定监测对象心脏部位相对于雷达天线的距离，即为监测对象心脏部位所在的位置。
S4：采用相位相关的方法，从步骤S3所得位置中提取监测对象心脏部位的运动信息，并据此计算其加速度信号。该步骤具体如下：
根据公式
Figure BDA0003647563420000065
提取反映监测对象心脏部位与接收天线阵列间距离变化的相位序列
Figure BDA0003647563420000071
解缠绕后获得对应的体动信号序列{Rm}；其中，λ为所发射调频连续脉冲的波长；
随后通过公式
Figure BDA0003647563420000072
求得对应的加速度信号序列{am}；其中，Δt为慢时间轴的采样间隔。
上述步骤S3和S4统一为提取监测对象心跳信息，在实际应用时，例如，接收天线接收到第n个雷达回波信号经过上述预处理后得到的频谱，该频谱为一个复数序列，记此序列为xn，该序列中的第k个元素记为xn,k，通过反正切函数可以求得该复数元素的相位
Figure BDA0003647563420000073
Figure BDA0003647563420000074
其中Im(xn,k)表示复数xn,k的虚部，Re(xn,k)表示复数xn,k的实部。
选取时间窗长度为T，按上述预处理方法计算在时间窗长度T内的每个雷达回波信号的频谱
xn+1,xn+2,……,xn+T
并分别选取每个频谱序列中的第k个元素，按前述方法计算其相位
Figure BDA0003647563420000076
后组成一个新的序列Φk
Figure BDA0003647563420000075
对该序列进行FFT求得其频谱，并据此计算其在心跳频段[0.8Hz,2Hz]的总能量Ek；遍历上述过程中的k的取值即可得到监测范围内不同距离处物体在心跳频段能量值E
E＝[E1,E2,……,EN]
其中N表示预处理后得到频谱的点数。找出E中最大的一项(例如第p项)，即代表在该距离处心跳频段的能量最大，认为此最大项Ep对应的下标p即监测对象心脏所在的位置。选取每个预处理后的频谱中的第p项并按照上述反正切函数提取其相位，即可得到一个相位序列
Figure BDA0003647563420000081
对上述相位序列进行解缠绕后，即可得到代表监测对象体动信息的体动信号序列{Rm},m＝1,2,……,T；
所述监测对象心跳信息的提取是通过对监测对象的体动信号序列{Rm}计算二阶差分信号实现的。对于上述得到的体动信号序列{Rm}，按照下式逐项计算其加速度信号序列{am}
Figure BDA0003647563420000082
其中Δt为慢时间轴的采样间隔。
S5：由于体动信号中的呼吸信号与心跳信号相比虽然幅值更大，但是其变化过程相对缓慢。因此通过计算体动信号的加速度信号并对呼吸信号进行抑制同时对心跳信号进行加强。通过计算短时平均功率以对所得加速度信号进行平滑，并通过峰值检测的方法对各次心跳之间分割点的位置进行估计。
该步骤具体如下：
对心跳信号的分割首先通过计算功率对上述得到的加速度信号序列{am}进行平滑。定义在n0时刻的平均功率为
Figure BDA0003647563420000083
其中L为计算平均功率选取时间窗长度的一半。
按照下式，逐项处理加速度信号序列{am}，得到平滑后的心跳信号序列{Ym}：
Figure BDA0003647563420000084
其中L1，L2表示在计算平均功率时选取了不同的时间窗长度，取L1为单次心跳周期时长(例如0.8～1.2s)一半对应的采样点数，取L2＝0.4L1；
平滑后的心跳信号序列{Ym}是一个类正弦信号，通过寻找序列{Ym}中的峰值，并记录每个峰值点在序列中的下标：
t＝[t1,t2,……,tI]
其中I表示序列{Ym}中峰值点的个数。ti即估计得到的各次心跳之间分割点在序列中的下标。据此即可得到各次心跳的心跳间期。
S6：然而在进行上述平滑处理时，对于一段时间内平均功率的计算会引起加速度信号序列{am}中部分细节的丢失，导致基于上述心跳分割点求得的心跳间期误差较大，不能满足心率变异性指标的监测要求。因此，需要对加速度信号序列{am}进行进一步精细化分割处理。即基于步骤S5中估计得到的各次心跳之间分割点的位置，生成单次心跳的加速度模板信号；随后利用该加速度模板信号对步骤S4中所得加速度信号进行精密分割，得到各次心跳的心跳间期。
该步骤具体如下：
基于上述估计得到的各次心跳分割点所在位置t，对加速度信号序列{am}进行截断，分成如下片段
Figure BDA0003647563420000091
其中
Figure BDA0003647563420000092
表示[t1,t2]两个分割点之间加速度信号序列{am}对应的序列。
对上述加速度信号片段进行平均，得到心跳加速度信号的模板信号Template
Figure BDA0003647563420000093
对于长度不相同的加速度信号片段，需要通过插值的方式统一其长度。
用该模板信号Template对原加速度信号序列{am}重新进行分割。在原加速度信号序列{am}中选取一组分割点S＝[s1,s2,……,sI]。得到一组对加速度信号序列{am}的新的截断
Figure BDA0003647563420000094
通过下式定义新的截断与模板的匹配程度
Figure BDA0003647563420000101
当某一截断
Figure BDA0003647563420000102
与模板Template长度不同时仍然通过插值使该截断的长度与模板的长度保持一致。
反复调整分割点S直到在当前模板信号下D取值达到最小；记此时的分割点集为S1，基于S1按照上述方法生成新的模板信号，记为Template1，基于Template1重新对原加速度信号序列{am}进行分割。直到计算得到的D小于事先设定的阈值，即认为完成了对当前加速度信号序列{am}的分割。此时的分割点集
S＝{s1,s2,……sI}
即为最终各次心跳之间分割点的位置。
S7：基于步骤S6所得的心跳间期，计算心率变异性指标，进而实现非接触式心率变异性监测。
具体的，据此计算各次心跳间期RR
RR＝[s2-s1,s3-s2,……,sI-sI-1]
基于上述心跳间期RR，按照下式计算心率变异性相关指标MEAN，SDNN等
Figure BDA0003647563420000103
Figure BDA0003647563420000104
实施例
为使本发明实现的技术手段、创作特征、达成目的与功效易于明白了解，下面结合具体实施方式，进一步阐述本发明。
本实施例中提出了一种基于调频连续波雷达的非接触式心率变异性监测方法，具体如下：
S1向监测对象所在区域发射周期性的线性调频波，并使用天线采集相应的回波信号。所发射的线性调频波的线性调频范围为77-81GHz，单个脉冲持续周期为50us；快时间轴采样频率为4MHz，采样点数为128；慢时间轴采样频率为100Hz。
S2采用混频和FFT的方式对雷达回波信号进行预处理，混频输出为一个中频信号xout。
Figure BDA0003647563420000115
该中频信号的幅值A为常数，频率和相位均正比于反射面和天线的距离R，如下式所示：
Figure BDA0003647563420000111
其中S为线性调频连续波频率随时间的变化率，c为光速，λ为所用调频连续波的波长，在本例中，S＝70MHz/μs，λ＝4mm。
最后对混频输出的中频信号进行FFT，得到频谱中不同频率的成分即代表不同距离处的物体反射的雷达回波信号。
S3从雷达的回波信号中提取监测对象心脏部位所在的位置，具体指相对于天线的距离
S3-1对于接收到的第n个雷达回波信号，经过预处理后通过如下反正切函数的方式提取其第k个元素xn,k的相位
Figure BDA0003647563420000112
Figure BDA0003647563420000113
其中Im(xn,k)表示复数xn,k的虚部，Re(xn,k)表示复数xn,k的实部。
S3-2选取长度为T的时间窗并依次按照S3-1中的方法处理时间窗中的所有雷达回波信号，得到一个新的序列Φk
Figure BDA0003647563420000114
S3-3对S3-2中的序列Φk解缠绕并通过FFT计算其频谱，统计其在心跳频段[0.8Hz,2Hz]的能量Ek
S3-4遍历k并找到Ek的最大值，并据此确定监测对象心脏部位相对于雷达天线的距离。
S4根据
Figure BDA0003647563420000121
提取反映监测对象心脏部位与天线阵列间距离变化的相位序列
Figure BDA0003647563420000122
解缠绕后获得对应的运动信息序列{Rm}；并通过
Figure BDA0003647563420000123
求得对应的加速度信号序列{am}
其中Δt是慢时间轴的采样间隔
S5通过计算短时平均功率对S4中的加速度信号进行平滑，并通过峰值检测的方法对各次心跳之间分割点的位置进行估计
S5-1取L1＝50,L2＝20，通过
Figure BDA0003647563420000124
Figure BDA0003647563420000125
对S4中得到的加速度信号序列{am}进行平滑，得到平滑后的心跳信号序列{Ym}，{Ym}应为一个类正弦信号。
S5-2对于平滑后的心跳信号序列{Ym},通过峰值检测的方式找出其所有的极大值点在序列中对应的下标t
t＝[t1,t2,……,tI]
认为ti即是各次心跳之间的分割点对应的下标，其中I表示序列{Ym}中峰值点的个数
S6基于S5中估计得到的各次心跳之间分割点生成单次心跳的加速度模板信号，利用该模板信号对S4中的加速度信号进行精密分割，得到各次心跳的心跳间期。
S6-1基于S5中得到的各次心跳周期分割点的估计值t，将加速度信号序列{am}进行分割，得到如下I-1个数据段
Figure BDA0003647563420000131
S6-2对S6-1中加速度信号片段进行平均，得到心跳加速度信号的模板信号Template
Figure BDA0003647563420000132
通过插值保证S6-1中各个加速度信号序列{am}的片段长度一致。
S6-3选取一组新的分割点S
S＝[s1,s2,……,sI]
得到一组对加速度信号序列{am}的新的截断
Figure BDA0003647563420000133
S6-4通过下式定义新的截断与模板的匹配程度
Figure BDA0003647563420000134
通过插值保证每个截断
Figure BDA0003647563420000135
与模板Template长度相同。
S6-5：
重复S6-3～S6-4所述过程计算每组分割点集对应的截断与加速度模板信号的匹配程度，找出D最小的一组分割点集Ss作为当前模板信号下的最优分割。基于Ss按照S6-2所述过程生成新的模板信号
S6-6：重复S6-5所述过程，直到某次重复过程中S6-4计算得到的D小于事先设定的阈值，此时的分割点集S即为最终得到各次心跳之间分割点的位置。基于该分割点集S计算各次心跳的心跳间期。
S7基于S6中得到的心跳间期按照心率变异性定义公式，计算心率变异性指标。以MEAN指标为例
Figure BDA0003647563420000141
采用上述实施例的方法对监测对象的心率变异性进行监测，得到如下结果
图3为各次心跳间期测量效果展示，图中红色曲线表示通过电脑化心电图机采集的心电信号得到的各次心跳间期，为基准信号；蓝色曲线表示基于本专利提出的方法，通过调频连续波雷达测量得到的各次心跳间期RR
RR＝[s2-s1,s3-s2,……,sI-sI-1]
表1为在6组时长为5min的实验数据中心率变异性MEAN指标的测量结果。
表1
Figure BDA0003647563420000142
以上所述的实施例只是本发明的一种较佳的方案，然其并非用以限制本发明。有关技术领域的普通技术人员，在不脱离本发明的精神和范围的情况下，还可以做出各种变化和变型。因此凡采取等同替换或等效变换的方式所获得的技术方案，均落在本发明的保护范围内。

Claims (7)

1.一种基于调频连续波雷达的非接触式心率变异性监测方法，其特征在于，具体如下：
S1：向监测对象所在区域发射周期性的线性调频波，并使用接收天线采集相应的雷达回波信号；
S2：依次采用混频和快速傅里叶变换的方式对所采集的雷达回波信号进行预处理；
S3：从预处理后的雷达回波信号中提取监测对象心脏部位所在的位置；
S4：采用相位相关的方法，从步骤S3所得位置中提取监测对象心脏部位的运动信息，并据此计算其加速度信号；
S5：通过计算短时平均功率以对所得加速度信号进行平滑，并通过峰值检测的方法对各次心跳之间分割点的位置进行估计；
S6：基于步骤S5中估计得到的各次心跳之间分割点的位置，生成单次心跳的加速度模板信号；随后利用该加速度模板信号对步骤S4中所得加速度信号进行精密分割，得到各次心跳的心跳间期；
S7：基于步骤S6所得的心跳间期，计算心率变异性指标，进而实现非接触式心率变异性监测。
2.根据权利要求1所述的基于调频连续波雷达的非接触式心率变异性监测方法，其特征在于，所述步骤S1中，线性调频波的脉冲调频范围为77GHz～81GHz，单个脉冲持续周期为50us；快时间轴采样频率为4MHz，采样点数为128；慢时间轴采样频率为100Hz。
3.根据权利要求1所述的基于调频连续波雷达的非接触式心率变异性监测方法，其特征在于，所述步骤S2中，经混频方式处理后的输出为一个中频信号；该中频信号的幅值为一常数，频率和相位均正比于反射物体相对接收天线阵列的距离。
4.根据权利要求1所述的基于调频连续波雷达的非接触式心率变异性监测方法，其特征在于，所述步骤S3具体如下：
S3-1：对于接收到的第n个雷达回波信号，经过步骤S2预处理后提取其第k个元素的相位
Figure FDA0003647563410000011
S3-2：选取长度为T的时间窗，并依次按照步骤S3-1中的方法处理，得到序列Φk：
Figure FDA0003647563410000021
S3-3：对所得序列Φk解缠绕并通过快速傅里叶变换计算其频谱，统计其在心跳频段[0.8Hz,2Hz]的能量Ek；
S3-4：遍历k以找到Ek的最大值，并据此确定监测对象心脏部位相对于雷达天线的距离，即为监测对象心脏部位所在的位置。
5.根据权利要求1所述的基于调频连续波雷达的非接触式心率变异性监测方法，其特征在于，所述步骤S4具体如下：
根据公式
Figure FDA0003647563410000022
提取反映监测对象心脏部位与接收天线阵列间距离变化的相位序列
Figure FDA0003647563410000023
解缠绕后获得对应的体动信号序列{Rm}；其中，λ为所发射调频连续脉冲的波长；
随后通过公式
Figure FDA0003647563410000024
求得对应的加速度信号序列{am}；其中，Δt为慢时间轴的采样间隔，Rm+1为体动信号序列{Rm}的第(m+1)项，Rm为体动信号序列{Rm}的第m项，Rm-1为体动信号序列{Rm}的第(m-1)项。
6.根据权利要求1所述的基于调频连续波雷达的非接触式心率变异性监测方法，其特征在于，所述步骤S5具体如下：
S5-1：定义在n0时刻长度为L邻域内加速度信号的平均功率PL(n0)为
Figure FDA0003647563410000025
其中，L为计算平均功率选取时间窗长度的一半，ai为加速度信号序列{am}中的第i项；
取L1为单次心跳周期时长一半对应的采样点数，取L2＝0.4L1；按照下式，逐项处理加速度信号序列{am}，以得到平滑后的心跳信号序列{Ym}；
Figure FDA0003647563410000031
其中，L1和L2表示在计算平均功率时选取了不同的时间窗长度，
Figure FDA0003647563410000032
表示在m时刻L1长度邻域内加速度信号的平均功率；
Figure FDA0003647563410000033
表示在m时刻L2长度邻域内加速度信号的平均功率；
S5-2：对于平滑后的心跳信号序列{Ym}，通过峰值检测的方式找出其所有的极大值点在序列中对应的下标t
t＝[t1,t2,…ti…,tI]，
所得ti即是各次心跳之间的分割点在序列中对应的下标；其中I表示找到极大值点的总个数。
7.根据权利要求1所述的基于调频连续波雷达的非接触式心率变异性监测方法，其特征在于，所述步骤S6具体如下：
S6-1：基于步骤S5中所得各次心跳周期分割点的位置的估计值t，将加速度信号序列{am}进行分割，得到如下(I-1)个数据段：
Figure FDA0003647563410000034
S6-2：对步骤S6-1所得数据段进行平均，得到心跳加速度模板信号Template：
Figure FDA0003647563410000035
计算时通过插值保证S6-1中各个加速度信号的片段
Figure FDA0003647563410000036
长度一致；
S6-3：选取一组新的分割点集S：
S＝[s1,s2,……,sI]，
得到一组对加速度信号序列{am}的新的截断：
Figure FDA0003647563410000041
S6-4：通过如下公式定义步骤S6-3所得新的截断与加速度模板信号的匹配程度：
Figure FDA0003647563410000042
在计算时通过插值保证每个
Figure FDA0003647563410000043
与Template长度相同；
S6-5：重复步骤S6-3～S6-4以计算每组分割点集对应的截断与加速度模板信号的匹配程度，找出D最小的一组分割点集Ss作为当前模板信号下的最优分割；基于Ss按照步骤S6-2以生成新的模板信号；
S6-6：重复S6-5所述过程，直到某次重复过程中通过步骤S6-4计算得到的D小于事先设定的阈值，此时的分割点集S即为最终得到各次心跳之间分割点的位置。
CN202210539306.9A 2022-05-17 2022-05-17 一种基于调频连续波雷达的非接触式心率变异性监测方法 Pending CN114732390A (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202210539306.9A CN114732390A (zh) 2022-05-17 2022-05-17 一种基于调频连续波雷达的非接触式心率变异性监测方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202210539306.9A CN114732390A (zh) 2022-05-17 2022-05-17 一种基于调频连续波雷达的非接触式心率变异性监测方法

Publications (1)

Publication Number Publication Date
CN114732390A true CN114732390A (zh) 2022-07-12

Family

ID=82287764

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202210539306.9A Pending CN114732390A (zh) 2022-05-17 2022-05-17 一种基于调频连续波雷达的非接触式心率变异性监测方法

Country Status (1)

Country Link
CN (1) CN114732390A (zh)

Cited By (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN117731298A (zh) * 2024-01-05 2024-03-22 长春理工大学 一种基于fmcw雷达的非接触式心电信号反演方法

Cited By (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN117731298A (zh) * 2024-01-05 2024-03-22 长春理工大学 一种基于fmcw雷达的非接触式心电信号反演方法

Similar Documents

Publication Publication Date Title
Jezewski et al. A novel technique for fetal heart rate estimation from Doppler ultrasound signal
CN110584631A (zh) 一种基于fmcw雷达的静态人体心跳和呼吸信号提取方法
EP3446248A2 (en) Extraction of features from physiological signals
CN109009124B (zh) 基于超宽带雷达的呼吸频率测量及目标定位方法
CN114732390A (zh) 一种基于调频连续波雷达的非接触式心率变异性监测方法
CN114847911A (zh) 一种基于毫米波雷达的多人生命体征监测方法
CN114052744B (zh) 基于脉冲神经网络的心电信号分类方法
CN113273978B (zh) 一种基于超宽带雷达的人体呼吸和心跳频率的检测方法
WO2023226223A1 (zh) Ppg信号质量评估方法及装置以及ppg信号处理方法及系统
CN106510674A (zh) 血压信号去干扰的方法和装置、血压检测系统
WO2023004688A1 (zh) 一种基于多普勒雷达的非接触式呼吸监测方法
CN103690169A (zh) 呼吸信息检测方法及系统
Liu et al. A coarse-to-fine robust estimation of FMCW radar signal for vital sign detection
CN112168152A (zh) 呼吸心跳的检测方法、检测装置与计算机可读存储介质
Li et al. A novel approach to phase space reconstruction of single lead ECG for QRS complex detection
CN115736854A (zh) 一种基于毫米波雷达的呼吸心跳监测系统
CN117838083A (zh) 一种基于毫米波雷达的体征快速精确检测方法
CN113116320A (zh) 一种基于vmd的fmcw雷达生命信号检测方法
CN117100241A (zh) 一种心跳间期测量方法及装置
CN117530666A (zh) 呼吸异常识别模型训练方法、呼吸异常识别方法及设备
CN115015867B (zh) 一种基于超宽带雷达的身份识别和跌倒检测方法
CN114145725B (zh) 一种基于无创连续血压测量的ppg采样率估算方法
CN111693990B (zh) 一种基于24GHz雷达的简易手势识别方法
Lu et al. Accurate heart beat detection with Doppler radar using bidirectional GRU network
Zhang et al. Radar-Beat: Contactless beat-by-beat heart rate monitoring for life scenes

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