基于压电陶瓷传感器的非接触式精准逐拍心率提取方法研究

方 震 白忠瑞 陈贤祥 夏 攀 何征岭 赵荣建

①(中国科学院空天信息创新研究院 北京100190)

②(中国科学院大学电子电气与通信工程学院 北京100049)

1 引言

心冲击图(BallistoCardioGram,BCG)是一种无接触式的心脏活动监测手段。它来源于心脏泵血过程中血液流动在人体内部产生的一系列机械冲击力[1]。与心电图(ElectroCatdioGram,ECG)等其他生理信息监测技术相比，BCG避免了电极线和电极贴长时间贴在身上带来的束缚和不适，具有无创、无直接接触和检测方便等优势，十分适用于居家养老、病房等需要长时间监测生理信号的场景。近年来随着传感器和数字信号处理的发展，BCG逐渐受到研究人员的重视，研究表明，BCG可以应用于心率、心率变异性(Heart Rate Variability,HRV)监测[2,3]，也可应用于心脏收缩性以及心输出量变化等心血管功能的评估[4]。

但BCG信号在不同测试条件(如测试设备的安置位置与角度、受试者姿势等)下表现出较大的不一致性[5]，此固有特征在一定程度上限制了分析BCG的形状特征在评估心血管功能方面的应用潜力。相比之下，通过检测逐拍心动周期估计心率或者计算HRV指标，是BCG可行性较高的一种应用路径。

HRV定义为逐次心跳周期差异的变化情况，一般通过分析心电信号中R波间期序列获得。HRV反映了心脏本身窦性心律不齐的程度以及神经体液因素与窦房结之间相互作用的平衡关系。已有的研究表明，HRV是心源性猝死、冠心病、高血压病及慢性心力衰竭等心血管疾病及慢性阻塞性肺疾病、糖尿病等疾病预后的预测因子[6,7]，还能反映出睡眠状况、精神压力状况等多种信息[8,9]。因此使用无感式的BCG信号获取心率变异性指标，对于用户方便舒适地掌握自身健康状况具有重要意义。

通过BCG获取HRV的关键在于获得精准可靠的逐拍心率。近年来，已有很多学者在非接触式获取逐拍心率方面做了研究：在BCG信号采集装置方面，有压电薄膜、压电陶瓷、光纤、毫米波雷达、加速度计、压变电阻以及各类采用MEMS工艺的微传感器等[10–12]；在逐拍心动周期的算法方面，目前有自适应阈值法、峰值探测法、倒谱分析法、多示例学习法、模板匹配法、混合判决方法[13–16]。值得注意的是，目前此方面研究更多地关注逐拍心率检测的覆盖率和平均心率的准确度，但事实上，逐拍心动周期计算的精度，即平均绝对误差(Mean Absolute Error,MAE)对于心率变异性的计算也尤为重要。例如，在5 min短时HRV计算中，若BCG逐拍心动周期序列平均绝对误差为10 ms，则由其计算得到的频域特征HF相比ECG真实值的偏差可达15%～25%,SDNN,p NN50等时域指标偏差也可达10%左右，且此偏差具有随机性，这使得BCG所得到的HRV特征指标的可靠性降低；而当逐拍心动周期的MAE降至4 ms时，逐拍心动周期序列所得到的HF偏差降低到8%以内，SDNN和p NN50等指标偏差降低至2%。

针对目前大多数方法在逐拍心率计算精度方面的不足，本文完成以下工作：(1)通过采用优化的传感前端设计，增加传感器的灵敏度和BCG信号的时间分辨率；(2)通过对心冲击信号进行分析，找到BCG中最适合提取逐拍心动周期的关键成分；(3)提出一种采用聚类的自适应模板匹配算法，以准确提取心动周期信息。最后通过实验验证了方法的准确性。

2 方法

2.1 基于压电陶瓷的传感前端

BCG信号采集装置的设计影响逐拍心率提取的准确度。在各种BCG监测方式中，压电陶瓷传感器由于其灵敏度高、价格低廉、设计灵活等优势，具有良好的应用前景。本文通过压电陶瓷非接触式采集BCG，计算逐拍心率。

本文所设计的传感前端外壳为一个ABS材质塑料圆盒，压电陶瓷片被固定在圆盒内上壁形变最大位置。外壳上盖结构设计为一侧凸起形状，使其受到震动时响应更灵敏。实验中将传感前端放置在胸口下方的床垫下，采集到的身体震动信号会经放大、滤波以及模数转换传输到带有低功耗蓝牙(Bluetooth Low Energy,BLE)功能的主控芯片中进行下一步的传输和处理。BCG信号的采样频率也影响逐拍心率计算的精准性，经对比测试，本文设置BCG采样频率为250 Hz，以保证逐拍心率计算的准确程度，同时最小化资源消耗。

2.2 BCG信号关键成分提取

传感前端采集到的原始的压电信号可用式(1)表示

其中，r(n)表示呼吸引起的振动成分，b(n)表示BCG成分，g(n)表示噪声成分。心脏的跳动是一个较为复杂的机械运动过程[17]，BCG中可以反映心脏机械运动和心脏泵血引起的身体机械运动的很多信息。

但对于一般所测得的BCG信号，由心脏原始的泵血引起的信号相对细微，容易被淹没在由身体四肢震动造成的信号的主能量中，通过不同的处理方式可以凸显出BCG中的不同成分。图1是不同方式处理BCG信号的结果，虚线是同步采集的ECG的R波所在位置。图1(a)是经过48 Hz低通滤波器预处理去除高频噪声后的压电信号，为了求得误差最小的逐拍心动周期误差，我们尝试了多种不同的方式对此信号进行进一步处理，以获取BCG信号中最适合进行精准逐拍心率提取的“关键成分”，下面优选几种具有代表性的方法具体说明：

方法1使用3～24 Hz的巴特沃斯带通滤波器进行滤波，得到BCG概貌波形，其结果如图1(b)。

方法2使用与心冲击信号最相符的db6小波基对原始压电信号进行9级小波分解，并使用小波分解后的第3～7层细节分量重建BCG信号，其结果如图1(c)所示。

方法3使用8～24 Hz的巴特沃斯带通滤波器进行滤波，可以得到BCG信号中的由心脏跳动引起的细微振动成分，其结果如图1(d)所示。

方法4对方法1处理的结果进行二次差分操作，求取其加速度信号，结果如图1(e)所示。

对于不同的BCG处理方法，每拍心跳在信号上造成的最高点并不是同一位置，为方便起见，本文中统称这些点为BCG的“J点”。

2.3 逐拍心率提取算法

2.3.1模板学习

图1 BCG不同成分提取结果

BCG形状多变，难以提前预知其波形。我们使用了亲和传播(Affinity Propagation,AP)聚类算法[18]，从10 s的BCG信号中进行模板的自学习。具体步骤如下：

(1)提取2.2节各方法处理过的10秒BCG信号段中所有的极大值点，以这些极大值点为中心，截取其前后各100采样点，共约0.8 s长度的信号，得到数据样本集B={b0b1···bn}，其中bi代表第i个极值点处截取的信号。

(2)使用AP聚类算法对样本集B进行聚类，并得到模板聚类。AP聚类的具体步骤如下：

步骤1初始化吸引度矩阵R和归属度矩阵A为0矩阵。其中R中元素r(i,k)表示数据bk适合作为数据bi的聚类中心的程度，A的元素a(i,k)表示对象bi适合选择数据对象bk作为其聚类中心的程度。

步骤2分别按照式(2)和式(3)迭代更新吸引度矩阵与归属度矩阵。

其中，s(i,k)是相似度矩阵S的元素，为弱化BCG幅值变化的影响，本文没有使用欧氏距离，而是设置s(i,k)为Pearson相关系数与常数1的差，如式(4)所示。特殊地，i=k时的相似度s(k,k)称为参考度，影响聚类的数目，本文将s(k,k)设置为相似度的中位数M(S)，可有效完成聚类。最后得到的相似度矩阵S如式(5)所示。

步骤4重复步骤2、步骤3，直到矩阵稳定。按最大相似度规则将其他点分配到距离它最近的聚类中心点相应的聚类中，完成聚类。

(3)对各聚类，计算每个样本的中点值的平均值，选择此平均值最大的聚类为理想聚类，计算其算术平均作为最后的BCG模板，以BCG表示，如图2所示。

2.3.2心跳位置探测

通过将模板在滤波后BCG数据上不断右滑，形成相关系数函数。对于一段预处理后的BCG信号，采用以下步骤进行模板匹配，确定J点。

(1)通过以下方式构造原BCG信号等长的相关系数函数cor(x)。

(2)对前10 s BCG信号做FFT变换，然后进行频域寻峰确定心率近似值HRe。在最初的2 s的cor(x)信号中选取最大值，作为第1个J点的位置，以此为起点，向后搜索[60/(HRe+20)s,60/(HRe–20)s]范围内的最大值点作为下一个J点。当连续5个所选局部最大值小于阈值θ时，说明匹配质量下降，重新计算模板，本文中阈值θ设为0.75。

最后，根据式(9)和式(10)计算连续逐拍心动周期和逐拍心率。

其中，IndexJ(n)表示第n个J点的索引位置，Fs为采样频率。逐拍心动周期和逐拍心率之间可方便地转换，为更清楚地比较计算结果，本文使用逐拍心动周期进行误差分析。

3 实验与分析

3.1 实验方案

按照如图3所示方法进行实验。本文使用压电陶瓷设备采集BCG信号，同时为验证BCG逐拍心率计算的准确性，在本研究中同步采集单导联心电信号进行对比。实验开始时按单导联方式在受试者胸腹部佩戴一次性心电电极贴，随后让受试者安静平躺于实验床垫上，开启设备进行BCG与ECG的同步采集。通过2.2节和2.3节所述方法处理BCG和计算BCG的逐拍心动周期，再使用PT算法[19]定位ECG的R波，计算ECG逐拍心动周期。随后分析对于2.2节所述的不同的BCG成分，所提取出的逐拍J-J间期相对于R-R间期的覆盖率和精准度。得到每次心跳的R-R间期与J-J间期。所有数据分析使用Py Charm与Python3.5进行。

3.2 测试示例

图4(a)和图4(b)分别展示了受试者1的一段BCG信号在3～24 Hz带通滤波和8～24 Hz带通滤波的预处理方法下，使用本文第2节所述方法进行逐拍心率提取的结果，图中使用红色圆圈标出了所定位的J点的位置，虚线位置为ECG的R波位置。图5则展示了两种预处理方式所计算得到的逐拍心动周期与ECG逐拍R-R间期的对比图。可见经通带为3～24 Hz的滤波器预处理之后计算得到的J点定位在原始信号上即为每次心跳所引起的压电信号最大点位置，但如图5(a)所示，由此所得到的逐拍心率却与逐拍R-R间期有较大差别。而经通带为8～24 Hz的滤波器处理后，所提取的J点位置不再是原始压电信号每拍心跳的最大值，但却能得到与R-R间期更为吻合的J-J间期，如图5(b)所示。我们称这种现象为“关键成分”现象，即8～24 Hz的成分是去除了较低频的肢体振动等成分、更能反映心室泵血、血液冲击主动脉引起的心脏原生振动的“关键成分”，能够获得与R-R间期相近的J-J间期。

图2 BCG模板聚类与BCG模板

图3 非接触式逐拍心动周期检测准确度验证实验示意图

图4 不同频带滤波的BCG模板匹配计算结果

表1为对受试者1的5 min BCG数据，共计351次心跳统计的结果，心跳覆盖率定义为使用本文算法成功定位到的心跳(即定位J波与R波相隔50 ms内)的比例，平均绝对误差定义为逐拍J-J间期与RR间期之差的绝对值的平均值，而平均误差定义为逐拍J-J间期相对于R-R间期的平均误差。可见对于受试者1的5 min安静平躺的BCG数据，在分别采用4种预处理方式后，3～24 Hz滤波以及db6小波重构都取得了100%的心跳检出率，但却有较高的平均绝对误差。经8～24 Hz滤波和二次差分处理后，虽心跳检出率略有不足，但却取得了很低的平均绝对误差。

图6展示了两种不同预处理方式下受试者1的逐拍R-R间期与J-J间期Bland-Altman图，从中可以看出采用宽带宽滤波器所计算的结果(图6(a))在所有心跳间期分布范围内误差都明显高于窄带滤波器计算结果(图6(b))。两种预处理方式的95%置信区间分别为(–32.26 ms,32.26 ms)和(–12.34 ms,12.34 ms)，且两种处理方式所得逐拍心跳间期误差均值都为0，说明在给定时间段内使用BCG和ECG检测出了相同的心跳次数。此处需要说明的是，虽表1中8～24 Hz滤波器的处理方法心跳覆盖率低于100%，但这仅说明BCG有心跳检测位置与R波位置相差较多，总体依然是检测到了应有的心跳次数。

3.3 数据分析

本文采集了15名受试者的BCG和ECG数据，年龄23～34岁，身高156～182 cm，体重45～90 kg。表2展示了所有受试者各5 min安静平躺状态下采集到的数据的处理结果，并且展示了4种BCG成分对逐拍心率提取性能的比较。

由表2可见对于8～24 Hz的带通滤波和二次差分这两种方式处理，在检出心跳覆盖率方面略有不足，这是因为频率较高时波形表现得较为复杂，影响了模板匹配的效果；而同时这两种方式在平均绝对误差(MAE)方面有着较大的优势，这是因为高频成分更能反映心脏跳动的成分而不受整个身体震动的影响。

图5 ECG-RR间期与不同频率滤波的逐拍BCG-JJ间期折线图

表1 受试者1的数据处理结果对比

图7展示了2.1节所述各成分所计算得到的J-J间期与R-R间期误差绝对值的分布，可以看到对于前两种处理方式，75%的误差分别在20 ms和16 ms以下，而后两种处理方式的误差有75%在4 ms以下，且中位数、下四分位数、下限均为0 ms，说明有1/2以上的J-J间期具有与R-R间期相同的采样点数。

图6 受试者1的逐拍RR间期与JJ间期Bland-Altman图

表2 不同预处理方式下逐拍心率提取性能

图7 各成分误差分布箱型图

此外，本文还从系统整体的角度对比了同类研究工作在静息状态下测量BCG提取逐拍心率的主要性能指标，如表3所示。可见，本文方法能够在覆盖率和逐拍心动周期的准确度方面处于领先水平，尤其是本文MAE可达到4 ms以下，这有助于提高BCG计算HRV的可靠性。

4 结论

本文提出了一种基于压电陶瓷传感器采集BCG并提取精准逐拍心率的方法，该方法是通过优化传感器性能提高BCG信号质量，通过对比BCG信号的不同处理方式，找到最适合提取精准逐拍心率的成分，最后通过AP聚类和模板匹配的方法确定J点，计算逐拍心率。

表3 BCG逐拍心率提取方法性能对比

本文的逐拍心动周期平均绝对误差为3.78 ms，心跳检测覆盖率为98.5%。本研究结果表明，通过设计合适的BCG采集装置、采用合适的信号处理方法处理BCG信号以及选择合适的特征点提取算法，可以使BCG提取逐拍心动周期(J-J间期)相对于R-R间期的误差降至可接受的范围，这对于提高使用BCG进行健康监测的效果尤其是提高HRV计算的可靠性具有重要意义。对本文所提出的BCG“关键成分”成因进行实验仿真和验证，以及对BCG提取逐拍心率精度的影响因素作进一步量化分析将是下一步的研究方向。