CN112971771B - 一种基于毫米波重建心电图的方法 - Google Patents

Abstract

本发明公开了一种基于毫米波重建心电图的方法。本发明包括以下步骤：步骤1)毫米波雷达朝人体胸口位置发射调频连续波进行探测，并对探测信号进行锯齿调制；步骤2)根据弗里斯传播方程，将胸腔表面参数化为时变的非刚性表面；步骤3)获取胸腔表面纯净的中频信号，利用中频信号中的相位项进行细粒度胸部传感得到多维心脏传感信号；步骤4)重建心电图：将解调后的胸腔振动输入到训练好的时序神经网络中，重建多维度的心电图。采用本发明重建的心电图具有良好的保真度以及完整的心脏临床特征，在促进心脏监测应用领域拥有巨大潜力。

Description

一种基于毫米波重建心电图的方法

技术领域

本发明涉及无线感知领域，尤其涉及一种基于毫米波重建心电图的方法。

背景技术

随着人们生活水平的提高，无接触式的健康检测的需求也越来越大，尤其是对心率的检测，心率的变化可以反映人心肺健康状况。除了健康检测的应用，无接触式的心电计算也能应用到远距离人体探测方面。采用毫米波探测人体，具有成本低，距离远和不易受环境影响的优势。而目前，检测心电的技术都是接触式的，因此应用环境十分受限。

发明内容

为了解决背景技术中的问题，本发明提供了一种基于毫米波重建心电图的方法。

本发明采用的技术方案包括以下步骤：

步骤1)毫米波雷达朝人体胸口位置发射调频连续波进行探测，并对探测信号进行锯齿调制：

步骤2)根据弗里斯传播方程，将胸腔表面参数化为时变的非刚性表面：

步骤3)获取胸腔表面纯净的中频信号；

步骤4)利用中频信号中的相位项进行细粒度胸部传感得到多维心脏传感信号R；

步骤5)重建心电图：

5.1)选取不同活动中不同对象的心脏传感信号以及心脏传感信号对应的心电图数据作为样本数据；

5.2)构建时序神经网络，时序神经网络由长短期记忆网络(LSTM)和注意力机制组成；

5.3)将样本数据输入时序神经网络进行训练，训练时使用离散小波变换(DWT)系数之间的距离作为损失函数；

5.4)将待测对象的心脏传感信号R输入步骤4.2)训练后的时序神经网络中，输出对应的心电图数据，从而完成心电图的重建。

所述步骤1)中对探测信号采用下述计算公式进行锯齿调制：

其中，x为在时间t+kT时的信号幅度；T是周期；k表示周期的计数器，k＝0,1,2…；j为调制阶数，j＝1，2,3…；f_start表示探测信号的起始频率；B是带宽。

所述步骤2)中将胸腔表面参数化为时变的非刚性表面的计算公式为：

其中，y(t)表示时变的非刚性表面对调频连续波的反射函数；S表示胸腔表面，dS表示将胸腔表面网格化后分成的每个单元面；ρ_s(t)表示单元面的反射率；r_s(t)表示在t时刻单元面dS和雷达之间的距离。

所述步骤3)具体为：

3.1)探测信号被胸腔表面的单元面反射，不考虑信号传感过程中的杂波，通过计算探测信号和胸腔表面反射信号的差异得到纯净的中频信号；

胸腔表面反射信号表示如下：

其中，t_tr是反射信号的运行时间，δ是反射信号的振幅；

通过计算探测信号和胸腔表面反射信号的差异得到纯净的中频信号的计算公式如下：

δ’表示雷达中混频器转换损耗的振幅；

3.2)考虑到信号传感过程中的杂波，将含杂波的中频信号表示为旋转相量形式：

IF’_k＝A_kexp(jθ_k)+A_cexp(jθ_c)

其中，旋转相量中的第一项为纯净的中频信号IF(t+kT)，第二项为静态杂波；A_k表示第k个周期的振幅；θ_k表示第k个周期的相位，A_c和θ_c表示初始振幅和初始相位；由于调频连续波非常短，旋转相量重写时略去了周期T和时间参数t；

根据含杂波中频信号IF′_k的旋转相量形式判断出多个周期的含杂波中频信号分布在相量空间中中心位于A_c∠θ_c的圆周上，通过三个同周期的含杂波中频信号估计A_c∠θ_c，即求解

的外接圆圆心：

其中，‖·‖是一范数；

将A_c∠θ_c进行坐标转换获得A_cexp(jθ_c)，含杂波的中频信号IF′_k通过减去A_cexp(jθ_c)，得到抑制杂波后的中频信号；然后使用0.8Hz的高通滤波器对抑制杂波后的中频信号进行滤波后得到纯净的中频信号IF(t+kT)。

所述步骤4)具体为：

4.1)根据

对步骤3)中纯净的中频信号进行简化：

其中，f是

表示胸腔单元面dS的相位，并将

设为φ_s；r_s表示单元面dS和雷达之间的距离；δ’_s表示雷达中混频器转换损耗的振幅；cl表示毫米波在空气中的传播光速；由于简化，将中频信号IF(t+kT)中的t+kT换元为t；

4.2)将多个胸腔单元面的中频信号进行叠加：

其中，A、B分别表示叠加后正弦和余弦分量的幅度；s₀表示单元面的序号，IF_s(t)表示叠加后的中频信号；

4.3)将步骤4.2)叠加后的中频信号的相位值φ经傅里叶展开和自由空间传播模型得到心脏传感信号R，心脏传感信号R为不同单元面和雷达之间距离变化值的集合，

根据傅里叶展开级数确定n为16，

表示第n个单元面和雷达之间的距离。

所述步骤5.2)中，长短期记忆网络中输入层的输入为不同时刻下的心脏传感信号R，输出层的输出为t个中间变量β₀，…β_t，输入层输入和输出层输出之间的关系如下：

β₀，…β_t＝FoldR(CELL，[0，R₀，…，R_t])

其中，R_t是t时刻的心脏传感信号；CELL是LSTM单元；FoldR组合函数将LSTM单元和心脏传感信号依照CELL(…CELL(CELL(0,R0),R1)…,R_t)的方式组合后，β₀，…β_t为输出的t个中间变量。

所述的FoldR组合函数的组合方式具体为：LSTM单元接收前一时刻的CELL输出和当前时刻的心脏传感信号R并输出当前时刻的中间变量β，遍历所有时刻下的心脏传感信号R，对应输出t个中间变量β₀，…β_t；

其中，前一时刻的CELL输出为前一时刻中间变量的CELL输出；当前时刻的心脏传感信号为R₀时，前一时刻的CELL输出为0，即CELL(0,R0)。

所述步骤5.2)中，注意力机制定义为：

其中，Q是队列向量，M是密钥向量，V是值向量；d_M是M向量的模长，M^T表示M向量的转置；

将长短期记忆网络输出的t个中间变量β₀，…β_t输入注意力机制，即设置Q＝K＝V＝β₀，…β_t，以将注意力机制应用至LSTM提取的胸腔振动的时序信息上。

所述步骤5.3)中的损失函数具体采用symlet小波，定义为：

其中，x_i是第i个心脏传感信号，y_i是第i个心电图信号；τ是DWT系数；f是高通或低通滤波器；w是symlet小波中的参数，我们经验性地选择n＝1,…,6以到达性能和精度的最优平衡。

本发明的有益效果：

本发明具有可扩展性：本发明可执行远程无接触式心脏监测。因此，用户可以摆脱现有方案中笨重的电极。相对现有技术，皮肤状况不会影响监视质量，从而使更多的潜在应用成为可能。

本发明具有高精度：复杂的胸腔振动是心脏机械活动的重要生命体征，本发明可以感应整个胸部的细微振动。相对现有基于智能穿戴设备的心脏监测，本发明保留了心脏电活动的大部分信息，实现了多维心脏监测。

附图说明

图1为本发明的工作系统示意图。

图2为本发明方法的流程图。

具体实施方式

下面结合附图及具体实施例对本发明作进一步详细说明。

一、心脏机械活动与心电图之间是通过心肌细胞的兴奋收缩耦合效应建立联系的，而毫米波能测量心肌细胞的拉力的变化，下述过程证明了从胸部采集毫米波进而重建多通道的心电图的可行性：

1)心脏兴奋是有组织的心肌细胞动作电位变化，心肌动作电位共分为三种状态：静息状态、去极化状态、复极化状态。心肌细胞的动作电位主要受细胞内外的钠离子、钾离子、氯离子浓度控制。

从医学上讲，戈德曼方程可以描述动作电位与三种离子浓度的关系：

其中，E_m(t)是心肌细胞动作电位，R是气体常量，T是绝对温度，z是离子化合价，F是法拉第常量，P是离子的渗透性，a_o，a_i是细胞内外的离子浓度，K、Na、Cl分别表示钾离子、钠离子、氯离子。

2)所有部位心肌动作电位的总和构成了心电图的曲线，可以被描述为：

其中

为与心肌细胞动作电位的梯度成正比的量；

为心肌细胞动作电位的梯度；

表示心肌细胞的位置和心电图测量点的距离；(i，j，k)表示心电图上每个点的坐标。

3)细胞内钙离子浓度的模型通过下述公式计算：

其中η_c表示心肌动作电位控制的通道效率。c表示心脏膜通道。

4)心脏机械活动的实质是心肌细胞的收缩，而钙离子浓度会影响心肌细胞的收缩。当粗肌丝和细肌丝发生滑动时，心肌细胞会出现收缩现象。具体来说，细肌丝最初被肌钙蛋白复合物“锁定”，所以不会发生滑动，随着越来越多的钙离子进入细胞内部，细肌丝被释放，在这种情况下，粗肌丝附着在细肌丝上开始收缩。在收缩过程中，粗肌丝消耗能量将细肌丝拉向细胞中心，心肌细胞的拉力是所有肌丝的总和：

表示心肌细胞的拉力，τ表示一条肌丝产生的细胞拉力分量

随后在动作电位的控制下，细胞内钙离子浓度降低，肌丝再次被“锁定”。在心脏周期中，心房的心肌细胞首先填满心室，然后心室收缩向肺动脉和主动脉传送血液。所以心脏周期诱导的振动可以通过胸腔传播，导致整个胸部表面微小和复杂的皮肤振动。

5)实验证明通过同时采集毫米波和心电图，发现二者具有相同的周期，并且毫米波包含四倍的心电图的样本熵，即毫米波可以采集更多的信息。

基于毫米波重建心电图的可行性，多维胸部传感和完整的心电图重建，具体过程如图2所示：

二、多维胸部传感获取心脏传感信号

1)胸腔表面振动是心脏收缩的结果，因此我们根据弗里斯传播方程，将胸腔表面参数化为时变的非刚性表面。

其中，S代表胸腔表面；ρ_s(t)表示单元面的反射率，ρ_s(t)被毫米波的方向和频率影响，它典型的取值是在34％-42％之间；r_s(t)表示在t时刻单元面dS和雷达之间的距离；我们的目标是感知和时间相关的胸腔表面的距离分布R_s(t)

2)为了确定非刚性胸腔表面的距离分布，我们选择频率调制连续波(FMCW)进行探测。用于探测的信号被锯齿调制为：

其中，f_start表示起始频率，B是带宽，T是周期，k＝0,1,2…表示周期的计数器。

如果探测信号被胸腔表面元素反射，那么反射信号可以表示为：

其中，t_tr是反射信号的运行时间，δ是反射信号的振幅；

根据公式(5)转换得到：

根据速度、距离和传播时间的关系可以转换得到：

其中，cl表示毫米波在空气中的传播光速。

3)反射信号可以视作原信号的延迟版本。我们发现它们都携带距离分布的信息。但如弗里斯传播方程中所描述的，信号的传播路径不理想，我们利用反射信号的指数部分做进一步的处理。特别的，我们通过计算探测信号和胸腔表面反射信号的差异来得到中频信号：

δ’表示雷达中混频器转换损耗的振幅。公式中相比于频率分量，相位项对细小的距离分布具有更强的敏感性，所以放射性心电图可以利用相位项进行细粒度的胸腔传感。相位项对距离分布更敏感同时也会导致容易产生杂波。

4)在实践中，简短的时间内获得的中频信号可以用向量重写：

IF′_k＝A_kexp(jθ_k)+A_cexp(jθ_c) (11)

第一项为纯净的中频信号，第二项为静态杂波；A_k表示第k个周期的振幅；θ_k表示第k个周期的相位，A_c和θ_c表示初始振幅和初始相位；由于调频连续波非常短，向量重写时略去了周期T和时间参数t；

我们的目标是消除静态杂波。从上边的等式中，我们可以发现真实的中频信号分布在一个中心是A_c∠θ的弧形上，因此我们使用以下的等式通过三个中频信号数据点来估计A_c∠θ：

其中，‖·‖是一范数。通过减去A_cexp(jθ)，杂波得到抑制。然后应用0.8Hz的高通滤波器来移除呼吸的干扰。处理完毕后，我们得到了纯净的中频信号。

5)理想情况下我们可以对来自单一源头的中频信号进行建模，但是实际中，心脏收缩引起的振动会传递到胸腔表面的各个部位，形成各种振动。为了获得心脏活动的积分信息，我们研究了胸腔表面多个来源的相位历史，即多个维度。这种情况下，胸腔表面的中频信号可以表示为表面积分，由于这个积分理论上很难展开，且存在有限的反射源，所以我们简化了这个问题以更好地描述振动叠加过程，首先我们简化中频信号方程：

其中，f是

是胸腔单元面dS的相位φ_s。考虑到胸腔表面有多个反射信号，叠加之后的结果为：

方程表明叠加之后依然保持正弦性质，实际应用中，由于内部传输的非线性失真和可能的多路径效应，信号会是不完美的正弦信号。我们利用基于方程的傅里叶展开和自由空间传播模型，将中频信号分解为多个维度，得到心脏传感信号R，心脏传感信号R为不同距离变化值的集合，

根据傅里叶展开级数确定n为16；

三、心电图重建

由于心脏活动涉及到两百万以上的心脏细胞及其参数，我们提出基于时序神经网络来以较少的参数量从分解的胸部振动中检索完整的心脏活动信息，从而重建心电图。心电图重建分为两步，首先使用同步采集的毫米波多维心脏传感信号和心电图信号训练该时序神经网络。在此过程中，时序神经网络可以学习到心脏传感信号和心电图之间的隐式映射。这一步骤只进行一次，得到的神经网络作为心电图重建模块可以重复使用。在之后的心电图重建阶段，只需将心脏传感信号输入到训练好的时序神经网络，重建心电图。

1)选取不同活动中不同对象的心脏传感信号R以及心脏传感信号R对应的心电图数据作为样本数据。

2)构建时序神经网络，在时序神经网络的架构设计方面，我们重点考虑心电图重建基于的心脏活动的两个性质：高阶马尔可夫性和周期稳定性。

2.1)我们使用长短期记忆网络结构(LSTM)来捕捉胸部振动的高阶马尔可夫性。

β₀,…β_t＝FoldR(CELL,[0,R₀,…,R_t]) (15)

R_t是t时刻胸部振动，CELL是LSTM单元，该单元接收一个另CELL输出和一个时刻的胸部振动并输出一个中间变量β，FoldR组合函数把LSTM单元和胸部振动按照CELL(…CELL(CELL(0,R0),R1)…,R_t)的方式组合，得到t个中间变量β₀,…β_t。

2.2)虽然长短期记忆网络能够利用高阶马尔可夫特性挖掘胸部振动中的信息，但它不能处理周期性信息。由于毫米波多维传感和心电图都是高采样率的信号，需要大量的数据点才能显示周期性。在这种情况下，长短期记忆网络在训练过程中存在梯度消失的问题。因此，我们利用注意力机制，通过周期稳定性来检索长期心脏活动信息：

其中Q是队列向量，M是密钥向量，V是值向量，d_M是M向量的模长。我们将LSTM结构中得到的t个中间变量输入到Q中，并设置Q＝K＝V＝β₀,…β_t以将注意力机制应用到胸部振动的周期性信息上。

3)除架构设计之外，目标函数在神经网络构建中也起着重要的作用。在心电图重建中，我们不仅关注心脏传感信号与真实心电图之间的全局平均距离，而且关注局部临床特征。由于最广泛使用的平均平方误差通常会使具有重要临床参考价值的信号峰谷变平，我们建议使用离散小波变换(DWT)系数之间的距离作为损失函数。具体来说，我们考虑了symlet小波，它的波形接近心电图信号，能够有效地为峰值和低谷分配更大的权重。此外，symlet小波还有助于避免噪声对心电信号的影响。因此，心电图重建中的目标函数被定义为：

其中，x_i是第i个心脏传感信号，y_i是第i个心电图信号，τ是DWT系数。f可以是高通或低通滤波器；我们经验性地选择n＝1,…,6以到达性能和精度的最优平衡。利用所提出的优化目标函数进行训练得到的时序神经网络很好地建模了毫米波心脏多维传感信号到心电图信号的映射，可在后续心电图重建中作为心电图重建模块重复使用。

如图1所示，由本发明构成的无线心电图检测设备包括毫米波发射器、毫米波传感器、信号处理模块和心电图重建模块。毫米波发射器用于发射调频连续波进行探测；毫米波传感器用于接收反射的毫米波并产生中频信号；信号处理模块用于利用中频信号中的相位项进行细粒度胸部传感；心电图重建模块通过利用信号处理模块获得的心脏传感信号进行心电图重建，从而获得多维心电图，通过心电图获取心率等指标。

Claims (7)

1.一种基于毫米波重建心电图的方法，其特征在于，包括以下步骤：

步骤1)毫米波雷达朝人体胸口位置发射调频连续波进行探测，并对探测信号进行锯齿调制；

步骤2)根据弗里斯传播方程，将胸腔表面参数化为时变的非刚性表面；

步骤3)获取胸腔表面纯净的中频信号；

所述步骤3)具体为：

3.1)探测信号被胸腔表面的单元面反射，不考虑信号传感过程中的杂波，通过计算探测信号和胸腔表面反射信号的差异得到纯净的中频信号；

胸腔表面反射信号表示如下：

其中，t_tr是反射信号的运行时间，δ是反射信号的振幅；x为在时间t+kT时的信号幅度；T是周期；k表示周期的计数器，k＝0,1,2…；j为调制阶数，j＝1,2,3…；f_srart表示探测信号的起始频率；B是带宽；

通过计算探测信号和胸腔表面反射信号的差异得到纯净的中频信号IF(t+kT)的计算公式如下：

δ’表示雷达中混频器转换损耗的振幅；

3.2)考虑到信号传感过程中的杂波，将含杂波的中频信号IF′_k表示为旋转相量形式：

IF′_k＝A_kexp(jθ_k)+A_cexp(jθ_c)

其中，旋转相量中的第一项为纯净的中频信号IF(t+kT)，第二项为静态杂波；A_k表示第k个周期的振幅；θ_k表示第k个周期的相位，A_c和θ_c表示初始振幅和初始相位；

根据含杂波中频信号IF′_k的旋转相量形式判断出多个周期的含杂波中频信号分布在相量空间中中心位于A_c∠θ_c的圆周上，通过三个同周期的含杂波中频信号估计A_c∠θ_c，即求解

的外接圆圆心：

其中，‖·‖是一范数；

将A_c∠θ_c进行坐标转换获得A_cexp(jθ_c)，含杂波的中频信号IF′_k通过减去A_cexp(jθ_c)，得到抑制杂波后的中频信号；然后使用0.8Hz的高通滤波器对抑制杂波后的中频信号进行滤波后得到纯净的中频信号IF(t+kT)；

步骤4)利用中频信号中的相位项进行细粒度胸部传感得到多维心脏传感信号R；

步骤5)重建心电图：

5.1)选取不同活动中不同对象的心脏传感信号以及心脏传感信号对应的心电图数据作为样本数据；

5.2)构建时序神经网络，时序神经网络由长短期记忆网络和注意力机制组成；

5.3)将样本数据输入时序神经网络进行训练，训练时使用离散小波变换系数之间的距离作为损失函数；

5.4)将待测对象的心脏传感信号R输入步骤5.3)训练后的时序神经网络中，输出对应的心电图数据，从而完成心电图的重建。

2.根据权利要求1所述的一种基于毫米波重建心电图的方法，其特征在于，所述步骤1)中对探测信号采用下述计算公式进行锯齿调制：

其中，x为在时间t+kT时的信号幅度；T是周期；k表示周期的计数器，k＝0,1,2…；j为调制阶数，j＝1,2,3…；f_start表示探测信号的起始频率；B是带宽。

3.根据权利要求1所述的一种基于毫米波重建心电图的方法，其特征在于，所述步骤2)中将胸腔表面参数化为时变的非刚性表面的计算公式为：

其中，y(t)表示时变的非刚性表面对调频连续波的反射函数；S表示胸腔表面，dS表示将胸腔表面网格化后分成的每个单元面；ρ_s(t)表示单元面的反射率；r_s(t)表示在t时刻单元面dS和雷达之间的距离。

4.根据权利要求1所述的一种基于毫米波重建心电图的方法，其特征在于，所述步骤4)具体为：

4.1)对步骤3)中纯净的中频信号进行简化：

其中，f是

S表示胸腔表面；

表示胸腔单元面dS的相位，并将

设为φ_s；r_s表示单元面dS和雷达之间的距离；δ’_s表示雷达中混频器转换损耗的振幅；cl表示毫米波在空气中的传播光速；

4.2)将多个胸腔单元面的中频信号进行叠加：

其中，A₀、B₀分别表示叠加后正弦和余弦分量的幅度；s₀表示单元面的序号，IF_s(t)表示叠加后的中频信号；

4.3)将步骤4.2)叠加后的中频信号的相位值φ经傅里叶展开和自由空间传播模型得到心脏传感信号R，心脏传感信号R为不同单元面和雷达之间距离变化值的集合，

5.根据权利要求1所述的一种基于毫米波重建心电图的方法，其特征在于，所述步骤5.2)中，长短期记忆网络中输入层的输入为不同时刻下的心脏传感信号R，输出层的输出为t个中间变量β₀,…β_t，输入层输入和输出层输出之间的关系如下：

β₀,…β_t＝FoldR(CELL,[0,R₀,…,R_t])

其中，R_t是t时刻的心脏传感信号；CELL是LSTM单元；FoldR组合函数将LSTM单元和心脏传感信号依照CELL(…CELL(CELL(0,R0),R1)…,R_t)的方式组合后，β₀,…β_t为输出的t个中间变量。

6.根据权利要求1所述的一种基于毫米波重建心电图的方法，其特征在于，所述步骤5.2)中，注意力机制定义为：

其中，Q是队列向量，M是密钥向量，V是值向量；d_M是M向量的模长，M^T表示M向量的转置；

将长短期记忆网络输出的t个中间变量β₀,…β_t输入注意力机制，即设置Q＝K＝V＝β₀,…β_t。

7.根据权利要求1所述的一种基于毫米波重建心电图的方法，其特征在于，所述步骤5.3)中的损失函数具体采用symlet小波，定义为：

其中，x_i是第i个心脏传感信号，y_i是第i个心电图信号；τ是DWT系数；f[]是高通或低通滤波器；w是symlet小波中的参数。

Images