基于毫米波雷达的心音检测

王 浩，张兴敢

（南京大学 电子科学与工程学院，江苏 南京 210023）

0 引 言

随着毫米波雷达的发展，基于毫米波雷达的非接触式生命信号检测逐渐成为研究的热点。毫米波雷达不仅可以用来连续检测呼吸速率和心跳速率[1-2]，也可以对人体进行睡眠监测[3]。心音信号反映了心脏及主血管的机械运动，心脏的病变也会引起心音的病变，因此有效的心音信息可以辅助诊断心脏疾病[4]。心脏经历一次完整的收缩期和舒张期的时间称作心动周期，心音信号伴随着心脏的跳动重复出现。每个心动周期出现的心音信号通常都包含有第一心音（S1）和第二心音（S2），但这两种心音的产生原因完全不同。第一心音的产生是由于心室收缩时二尖瓣、三尖瓣的闭合；第二心音则是由于心室舒张时主、肺动脉瓣的闭合。LFM（Linear Frequency Modulated）毫米波雷达可以根据瓣膜振动引起的胸腔微小位移实现非接触式的心音信号获取。心音信号的精确分割是对S1、S2 时域上的精确定位，是心音分类[5]、异常心音识别[6]、杂音检测[7]等心音信息分析的基础，也是心音特征提取的关键。

心音信号的分割算法一直是心音研究的重点，文献[8]采用基于阈值和K 类均值聚类的心电分析方法识别心音。文献[9]通过心音信号的短时能量谱和自回归参数提取信号的包络，利用多层感知器神经网络和第一心音和第二心音的频谱特征提取基本的心音，结合舒张和收缩间期变化完成儿童的心音分割。这些基于时域包络特征的心音分割算法会受到提取包络方法的影响，例如通过希尔伯特变换提取的包络有较大的毛刺，从而影响心音分割的准确度。之后又有学者把心音分割问题建模为隐马尔可夫模型。文献[10]在隐马尔可夫模型的基础上将心音信号的分割问题建模为连续时段隐马尔可夫模型（Duration-dependent Hidden Markov Model,DHMM），利用测试数据集计算模型参数，完成心音信号的分割。文献[11]将文献[10]的DHMM 进行了部分改进，提出了基于逻辑回归的隐半马尔可夫模型（Logistic Regression-Hidden Semi-Markov Model, LR-HSMM），对心音信号分割的起始时间选择更具有包容性，并采用拓展的维比特算法对LR-HSMM 进行译码。文献[12]提出的解码峰值检测算法与文献[11]的LR-HSMM 相似，解决了在第二心音（S2）的干扰下，对心音信号中的第一心音（S1）的准确定位问题。基于隐马尔可夫模型的研究需要大量的数据样本计算模型参数，导致算法比较复杂，分割效果受到模型参数限制。

本文以LFM 毫米波雷达为基础，首先利用雷达差拍信号的相位信息获得心脏收缩舒张时产生的胸腔位移，采用带通滤波器完成对胸腔位移信息中心音信号的分离；接着利用基于双特征的心音分割方法，以一阶导数和频率包络作为心音信号特征，通过最大类间方差准则计算阈值，完成心音信号中心音部分和非心音部分的初步分割；最后根据收缩期和舒张期的时间间隔的不同，对心音部分中第一心音和第二心音进行类型判决，将一个心动周期内的心音信号分割为第一心音、收缩期（除去第一心音）、第二心音、舒张期（除去第二心音）四部分。本文通过毫米波雷达采集不同对象的心音信号作为实验数据，并与LR-HSMM 算法进行对比。

1 心音信号检测

1.1 LFM 毫米波雷达系统

整个LFM 毫米波雷达系统主要分为两个子系统：一个是负责雷达信号发送和接收射频模拟子系统；另一个是负责对差频信号进行信号处理的数字子系统。射频模拟子系统主要包含发送和接收天线、模数转换器、混频器以及用于生成线性调频脉冲的频率合成器。数字子系统主要包含数字信号处理器、微型控制器。图1为LFM 毫米波雷达系统的一个简易模型，其中LNA 表示低噪声放大器，PA 表示功率放大器。

图1 LFM 毫米波雷达系统框图

LFM 毫米波雷达的发射信号为：

式中：f0为发射信号的起始频率；S为发射信号的调频斜率；T为发射脉冲信号的周期。在不考虑雷达散射截面积的情况下，毫米波雷达接收到的人体反射回波信号为：

式中：td表示雷达发射信号到达接收端的时延，td= 2Rc，c 为光速，R表示人体胸部与雷达之间的距离。由于在毫米波雷达的工作过程中人体是保持静止的，因此R=R0+ ΔR，R0为人体与毫米波雷达的距离，ΔR为胸腔运动位移。

发射信号sT(t)与接收信号sR(t)经过接收端混频器后可以得到差拍信号sb(t)，sb(t)可以表示为：

由式（3）可以得到差拍信号sb(t)的频率fb=Std=2RSc，相位φb为：

式中：λ为毫米波雷达发射的电磁波波长，λ= cf0。由式（5）可知，差拍信号sb(t)的相位变化Δφb和人体胸腔位移ΔR线性相关，差拍信号sb(t)的相位变化可以代表人体的胸腔运动。

1.2 心音信号提取

心音信号同心跳信号一样都属于微弱的生命信号，而心音信号的检测与心跳信号检测的最大区别在于两者振动的频率不同。借助差拍信号sb(t)的相位变化确定人体胸腔的运动位移后，根据心音信号所在的频率区间完成心音信号的提取，心音信号的提取流程如图2所示。

图2 心音信号的提取流程图

接收端混频后得到的差拍信号sb(t)经过A/D 采样和距离快速傅里叶变换后，得到一维距离像向量，代表快时间维。多个差拍信号经过相同处理后得到的多个距离像向量构成距离-时间二维矩阵W，矩阵的列代表每个距离单元的慢时间采样。确定人体所在的距离单元，需要计算二维矩阵W中每个一维距离像向量的功率，再计算二维矩阵每列的平均功率，以最大功率平均值所对应的距离单元作为人体所在距离单元，胸腔位移运动信息就包含在与人体距离单元数据的相位内。如图3 所示（图3 中的功率均归一化），平均功率最大值对应的距离为0.39 m。

图3 人体目标检测

在距离-时间二维矩阵W中的数据都为复数，可以通过反正切解调的方法获得矩阵W中人体目标所在距离列每一个数据点的相位。但是在LFM 毫米波雷达系统中存在一些非线性因素会影响反正切解调，其中最主要的是未知的直流偏移，而直流偏移的主要来源是静态物体反射的电磁波。假设一个复距离单元数据的虚部为I(t) + DCi，实部为R(t) + DCr，那么解调所得的相位φ(t)为：

式中：DCr、DCi分别为实部和虚部的直流项；φ(t) =I(t)R(t)，为真实相位值。为了消除直流项的影响，本文采用Spath 算法[13]估计直流项，并通过补偿的方式消除数据中的直流偏移，直流偏移补偿前后结果如图4所示。

图4 人体距离列数据的直流补偿

当相位大于π 或小于-π 时，会被相位解调函数投影在(-π,π)范围内，致使相位突变。通过相位解卷绕可以实现相位的连续变化，解卷绕后的相位变化便等同于胸腔的运动变化。由于心跳产生的心音信号非常微弱，需要对解卷绕后的相位进行差分运算，增强心音信号，即zk=φk-φk-1。心音信号中第一心音和第二心音的频率范围[14]大致为10～200 Hz，本文选用IIR 带通椭圆滤波器对差分后的相位进行滤波，分离心音信号。为避免人体低频呼吸信号的干扰，同时考虑不同个体之间的普适性，带通滤波器的频带范围定为20～60 Hz，分离出的心音信号如图5 所示。

图5 LFM 毫米波雷达采集的心音信号

2 基于双特征的心音分割

2.1 双特征阈值分割

在时域方面，假设心音信号为s(t)，则s(t)的一阶导数为：

当Δt较小时，可近似为：

式中Ts为心音信号s(t)的采样周期。s′(t)取绝对值之后得到心音信号s(t)的时域特征sp(t)，从图5 心音信号波形图中可以发现，由振动产生的第一心音和第二心音信号波形会产生明显的上下波动，且第一和第二心音的幅度明显高于非心音部分。根据文献[10]，第一心音和第二心音的持续时间小于收缩期（除去第一心音)和舒张期（除去第二心音），对应心音信号中S1 和S2 部分的时域特征sp(t)的值会明显高于非心音部分。在时频方面，通过短时傅里叶变换（Short Time Fourier Transform,STFT）对心音信号s(t)进行时频分析：

式中h(t)为窗函数，窗长的选择需要满足两个条件：一是要保证具有良好的时间分辨率；二是尽量不超过第一心音或是第二心音的持续时间，本文将窗长设定为100 ms。通过心音信号s(t)的短时傅里叶变换，得到心音信号频带范围为[f1,f2]的频率包络FEnv(t)为：

式中：频率f的间隔为Δf；频率包络FEnv(t)代表[f1,f2]内频率能量之和随时间的变化。心音信号中属于S1 和S2 时间段的频率能量会高于非心音时间段。本文设定了三种心音信号的频率包络特征：FEnv1、FEnv2 和FEnv3，它们代表了三种不同频带范围的频率包络，FEnv1 的 频 率 范 围 为20～30 Hz，FEnv2 的 频 率 范 围 为30～40 Hz，FEnv3 的频率范围为40～60 Hz。

心音信号的时域特征sp(t)和时频特征FEnv(t)构成了心音信号的一组特征向量，即gt=[sp,FEnv]。通过设定阈值，将心音信号在时间维度上初步分为两个部分：第一心音或第二心音部分以及非心音部分。令θ表示阈值向量，即θ=[θ1,θ2]，若sp(t) ＞θ1且FEnv(t) ＞θ2，即gt＞θ，则时刻t归为第一心音或第二心音部分；否则，时刻t归为非心音部分，其中θ1为特征向量sp对应的阈值，θ2为特征向量FEnv 对应的阈值。

2.2 计算阈值向量θ

计算阈值向量θ等价于计算阈值θ1和θ2。本文利用最大类间方差法[15]求解阈值θ1和阈值θ2（阈值θ2计算步骤与θ1相同）。θ1的计算过程如下：

1）对向量sp中的元素进行量化，所有元素值根据四舍五入原则只取小数点后三位，量化后的向量sp表示为sp=[sp,1,sp,1,…,sp,L]，其 中L为向量sp的元素个数。向量sp中元素值为sp,i的元素个数为ni，计算元素值为sp,i的元素出现的概率pi为：

初始化阈值θth，令θth= mean(sp)，其中函数mean(·)为求均值函数。

2）根据阈值θth可以将向量sp中的元素分为两部分：一部分为元素值大于θth的S1 或S2 部分；另一部分为元素值小于θth的非心音部分。计算sp中元素属于S1或S2 部分的概率P1以及非心音的概率P2：

3）计算sp中属于S1 或S2 部分的元素平均值μ1以及非心音部分的元素平均值μ2：

阈值为θth时，类间方差为：

4）令θth=θth+ 0.001，若θth大于sp的最大元素值则停止，否则，转至步骤2）。

经过以上步骤得到阈值θ1：

2.3 心音类型识别

图6 为经过2.1 节的方法处理后得到的简易心音信号时域分布模型。模型中的心音信号在时间维度被分为了不同长度的信号片段，其中黑色的心音信号片段代表第一心音或第二心音，白色部分代表收缩期（除去第一心音）或舒张期（除去第二心音）的非心音部分。只有确定黑色心音片段的心音类型之后才能完成心音信号的最终分割。根据文献[16]，一次完整的心动周期中心室收缩期持续时间一般小于舒张期，且第一心音和第二心音的起始分别对应着心室收缩期和舒张期的起始。根据图6 的心音信号时域分布模型，心音类型识别方法为：若tstart,i-tstart,i-1＜tstart,i+1-tstart,i，则shs,i为第一心音；若tstart,i-tstart,i-1＞tstart,i+1-tstart,i，则shs,i为第二心音，其中shs,i为第i个S1 或S2 部分信号片段，tstart,i为第i个S1 或S2 部分信号片段的起始时刻。心音类型识别完成后，一个心动周期内心音信号可以分割为4 种连续的时间阶段：S1、Sys、S2、Dia，它们分别对应第一心音、除去第一心音之外的收缩期、第二心音、除去第二心音之外的舒张期。

图6 心音信号时域分布模型

3 实验结果及分析

本文使用的毫米波雷达为德州仪器公司的AWR1642，工作频段为77～81 GHz。在PC 端完成配置雷达参数后，毫米波雷达差拍信号采样数据会通过DCA1000 采集卡传回电脑，在PC 端通过Matlab R2020b进行信号处理和心音分割。实验通过6 位不同的测试者采集心音信号数据，在实验过程中，测试者在近距离面对雷达而坐，同时毫米波雷达高度与测试者的心脏高度保持一致。毫米波雷达的参数配置为：调频斜率S为60.012 MHz·μs-1，脉冲起始频率f0为77 GHz，ADC 采样率fs为2 MHz·μs-1，脉冲持续时间T为50 μs，每条数据的采集总时长为10 s。

3.1 心音分割结果

对于每一条心音信号数据，根据2.1 节的内容，通过式（8）和式（10）分别提取心音信号s(t)的一阶导数s′(t)和频率包络FEnv(t)，之后对s′(t)取绝对值得到特征向量sp和FEnv，心音信号的时域特征sp(t)和三种频率包络如图7 所示（三种频率包络均归一化）。从图7b）中可以看出，与原始的心音信号相比，sp(t)放大了S1 和S2 部分与非心音部分之间的幅度差距，图7c）中频率包络FEnv(t)的变化也与心音信号的幅度包络变化一致，即S1 和S2 部分的频率能量高于非心音部分，将两种特征结合可以使得S1 和S2 部分与非心音部分之间的分割更加可靠。

图7 心音信号的时域特征和三种频率包络

根据2.2 节的步骤计算特征向量sp和FEnv 对应的阈值θ1和θ2，若sp(t) ＞θ1且FEnv(t) ＞θ2，则时刻t归为第一心音或第二心音部分；否则，时刻t归为非心音部分。利用阈值向量θ完成对心音信号s(t)中S1 或S2 部分和非心音部分在时间轴上的分割，如图8 所示。此时，心音信号被划分为持续时间不等的S1 或S2 部分以及非心音部分的信号片段，由于临床上心室的收缩期一般小于心室的舒张期，所以根据2.3 节心音类型的识别方法，将S1 或S2 部分的信号片段进行第一心音或第二心音的判决，最终将整个心音信号分割为4 种时间阶段：S1、Sys、S2、Dia。频率包络选定为FEnv1 时，心音分割结果如图9 所示。从图9 可以看出，本文所提出的算法对于心音信号的分割定位取得了很好的效果。

图8 双特征阈值分割

图9 心音信号分割

3.2 心音分割效果对比

为了准确评估所提算法对于心音信号的分割效果，以手动分割心音信号的结果作为参考标准，将本文所提算法与文献[11]中的LR-HSMM 算法进行对比，选取评价指标F1分数来评估两种算法，F1分数是衡量测试准确性的常用指标，它同时包含了精确率p和召回率r，计算公式为：

根据本文方法，6 位测试者心音分割的p、r和F1分数结果如表1 所示。两种算法的指标对比结果如表2 所示。表1 和表2 中结果为所有数据指标的平均值，黑体字表示最优F1分数，其中LR-HSMM 算法中的特征提取部分选用同态包络。从表1 中6 位不同测试者的平均F1分数可知，当频率包络选取FEnv1 时，本文方法对于心音信号的分割效果明显优于FEnv2 和FEnv3。由表2 可知，当频率包络为FEnv1时，无论是第一心音（S1）还是第二心音（S2），本文方法在p、r以及F1分数三项指标上均高于文献[11]的LR-HSMM 算法。其中对于第一心音的F1分数，本文方法达到了（92.29±0.94）%，比LR-HSMM 算法提高了7.07%；对于第二心音的F1分数，本文方法达到了（88.71±4.61）%，比LR-HSMM 算法提高了7.76%。

表1 心音分割的F1 分数 %

表2 两种算法的F1 分数对比 %

4 结 论

本文利用LFM 毫米波雷达获得人体胸腔位移信息，通过带通滤波提取人体的心音信号。利用心音信号的一阶导数和频率包络特征，同时结合心音间隔时间的生理特征完成对心音信号的分割定位。实验结果表明，本文提出的算法对于心音信号的分割定位具有较高的精确度，为后续的心音记录自动分析奠定了基础。下一步的研究将包括心音在内诸多人体身体特征，通过毫米波雷达应用于个人身份的识别，从而进一步扩大心音的用途。

注：本文通讯作者为张兴敢。