采用MUSIC算法的按需心率估算

心率(HR)监测是许多现有可穿戴产品和临床设备的关键特性。这些设备一般测量光电容积脉搏波(PPG)信号，为获得该信号，须利用LED照射人体皮肤，然后用光电二极管测量血流引起的反射光强度变化。PPG信号形态与动脉血压(ABP)波形相似，这使得该信号成为受科学界欢迎的非侵入性心率监测工具。PPG信号的周期性与心脏节律相对应。因此，可以根据PPG信号估算心率。然而，受血液灌流不良、环境光线以及最重要的运动伪像(MA) 1的影响，心率估算性能会降低。业界已提出许多信号处理技术来消除MA噪声，包括ADI公司的运动抑制和频率跟踪算法，通过使用一个靠近PPG传感器放置的三轴加速度传感器来实现。当没有运动时，最好能有一个按需算法来向跟踪算法提供快速且更精确的心率估算。本文改造了多信号分类(MUSIC)频率估计算法，以利用ADI医疗健康手表平台，根据手腕上的PPG信号实现高精度按需心率估算，图1所示为其框图。该图的细节将在后面的内容中说明。
图1. 利用腕上PPG信号的基于MUSIC的按需心率估计算法

ADI医疗健康手表提供的PPG信号概述
当LED发光时，血液和组织会吸收不同数量的光子，导致光电检测器检测到不同的结果。光电检测器测量血液脉动的变化并输出一个电流，该电流随后经放大和滤波以供进一步分析。 图2a显示了一个由交流(ac)和直流(dc)分量组成的一般PPG信号。PPG波形的直流分量检测组织、骨骼和肌肉反射的光信号，以及动脉和静脉血液的平均血容量。交流分量则表示心动周期的收缩期和舒张期之间发生的血容量变化，交流分量的基频取决于心率。图2b是来自 ADPD107 手表的PPG信号，这在之前的《模拟对话》文章中已介绍过。ADI多感知手表的目标是测量人体手腕上的多项生命体征。ADI手表有PPG、心电图(ECG)、皮肤电活动(EDA)、加速度(ACC)和温度传感器。本文仅关注PPG和ACC传感器。
现在我们仔细看看PPG和ABP波形的相似之处。ABP波形是由于左心室射出血液造成的。主压力沿全身血管网流动并到达多个部位，动脉阻力和顺应性的显著变化引起反射。第一个部位是胸主动脉和腹主动脉之间的接合处，其引起第一次反射，通常称为收缩晚期波。第二个反射部位是腹主动脉和髂总动脉之间的接合处。主波被再次反射回来，产生一个很小的下降，称为重搏切迹，这可以在第一次和第二次反射之间观察到。还有其他较小的反射，这些反射在PPG信号中被平滑掉2。本文的重点是心率估计，其仅取决于PPG信号的周期性，此算法不考虑PPG的确切形态。
图2a. 含交流和直流部分的典型PPG信号
图2b. ADI医疗保健手表PPG信号

PPG信号预处理
PPG信号易受周边组织的不良血液灌流和运动伪像的影响是众所周知的1。为将这些因素的影响降至最小，以免干扰随后的PPG分析和心率估计，须有一个预处理阶段。需要一个带通滤波器来消除PPG信号的高频成分（如电源）和低频成分（如毛细血管密度和静脉血容量的变化、温度变化等等）。图3a显示了滤波之后的PPG信号。使用一组信号质量指标来找到适合于按需算法的PPG信号第一个窗口。第一次检查涉及ACC数据和PPG信号，以确定是否能检测到一段无运动的数据，然后衡量其他信号质量指标。如果三个方向上存在高于ACC数据绝对值的特定阈值的运动，则按需算法会拒绝根据这样的数据窗口进行估计。下一信号质量检查是基于数据段特征的某种自相关。图3b显示了经滤波的PPG信号的一个自相关例子。可接受信号段的自相关表现出如下特性：具有至少一个局部峰值，并且对应于最高可能心率的峰值不超过某一数量；局部峰值从高到低递减，间隔时间递增；以及其他一些特性。仅计算与有意义的心率（位于30 bpm到220 bpm范围内）相对应的间隔时间的自相关。
当有足够的数据段连续通过质量检查时，算法的第二阶段就会使用基于MUSIC的算法算出准确的心率。
图3a. 经过带通滤波的图1b中PPG信号
图3b. 图2a中信号图的自相关

基于MUSIC的按需心率估计算法
MUSIC是一种基于子空间的方法，使用谐波信号模型，可以高精度地估算频率3。对于受到噪声破坏的PPG信号，傅立叶变换(FT)可能表现不佳，因为我们需要的是高分辨率心率估计算法。此外，FT将时域噪声均匀分布到整个频域中，限制了估算的确定性。使用FT很难在较大峰值附近观察到较小峰值4。因此，在本研究中，我们使用基于MUSIC的算法进行心率的频率估计。MUSIC背后的关键思想是噪声子空间与信号子空间正交，所以噪声子空间的零点会指示信号频率。下面的步骤是这种心率估计算法的总结：

• 从数据中删除平均和线性趋势
  
• 计算数据的协方差矩阵
  
• 对协方差矩阵应用奇异值分解(SVD)
  
• 计算信号子空间阶数
  
• 形成信号或噪声子空间的伪谱
  
• 找出MUSIC伪谱的峰值作为心率估计值
  
  

MUSIC必须应用奇异值分解，并且必须在整个频率范围内搜索频谱峰值。我们来看一些数学算式，以使上述步骤更清晰。假设经滤波的PPG信号有一个长度为m的窗口，表示为 xm 且 m ≤ L （其中L为给定窗口中经滤波PPG信号的总样本数）。那么，第一步是形成样本协方差矩阵，如下所示：
然后对样本协方差矩阵应用SVD，如下所示：
其中，U为协方差矩阵的左特征向量， Λ 为特征值的对角矩阵， V 为右特征向量。下标 s 和 n 分别代表信号和噪声子空间。正如之前提到的，使用信号已经通过信号质量检查阶段的先备知识，对基于MUSIC的算法进行修改以用于心率估计，因此预处理步骤之后，信号中唯一存在的频率成分是心率频率。接下来形成信号和噪声子空间，假设模型阶数只包含一个单音，如下所示：
其中 p = 2 为模型数。仅考虑有意义心率限值内的频率。这会大 大减少计算量，使嵌入式算法的实时实现成为可能。搜索频率 向量定义为：
其中，k为心率目标频率范围内的频点， L 为 xm(t) 中数据的窗口长度。然后，下面的伪谱使用噪声子空间特征向量找出MUSIC的峰值，如下所示。

这里使用伪谱一词，是因为它表明所研究信号中存在正弦分量，但它不是一个真正的功率谱密度。图4显示了基于MUSIC的算法处理5秒数据窗口得到的示例结果，在1.96 Hz处有一个很陡的峰值，换算为心率是117.6 bpm。
图4. 使用PPG数据的基于MUSIC估计的一个示例

基于MUSIC的按需心率估计算法的结果
我们已经在一个包含1289个测试案例(data1)的数据集上测试了该算法的性能，并且在数据开始时，测试对象被要求静止。表1给出了基于MUSIC算法的结果，并指出估计的心率是否在参考(ECG)的2 bpm和5 bpm精度范围内，以及估计时间的第50百分位数（中值）和第75百分位数。表1中的第二行显示了对于一个包含298个测试案例(data2)的数据集，存在周期性运动（如步行、慢跑、跑步）时该算法的性能。通过检测运动，如果任一数据被视为不可靠而遭到拒绝，或者是认为不受运动影响而精确估算得到心率，则认为该算法是成功的。在内存使用方面，假设缓冲区大小为500（即100 Hz时为5秒），对于目标频率范围（30 bpm至220 bpm），所需总内存约为3.4 kB，每次调用花费2.83周期。
表1. 基于MUSIC的按需心率估计算法的性能数值
结语
基于MUSIC的按需算法是ADI公司医疗保健业务部门生命体征监测小组提出的众多算法之一。在我们医疗健康手表中使用的按需算法与这里讨论的基于MUSIC的方法不同，前者的计算成本较低。ADI公司为传感器（嵌入式）和边缘节点提供软件和算法功能，使其从数据中获取有价值的信息，仅将最重要的数据发送到云端，让我们的客户和合作伙伴可以在本地做出决策。我们选择应用的标准是，其成果对于我们的客户来说非常重要，并且我们拥有独特的测量专业技术。本文只是对ADI公司研发的算法的简单介绍。凭借我们在传感器设计方面的现有专业知识，以及我们在生物医学算法开发（包括嵌入式和云计算）方面的努力，ADI公司将拥有独特的优势来为全球医疗健康市场提供最先进的算法和软件。
参考电路
1 Tamura, Toshiyo Tamura, Yuka Maeda, Masaki Sekine, 和 Masaki Yoshida. “可穿戴光电容积脉搏波传感器——过去和现在.” Electronics, 第3卷第2期，2014年。
2 R. Couceiro, P. Carvalho, R.P. Paiya, J. Henriques, I. Quintal, M. Antunes, J. Muehlsteff, C. Eickholt, C. Brinkmeyer, M. Kelm, 和 C. Meyer. “根据手指光电血管容积图的多高斯拟合评估心血管功能.” Physiological Measurement ，第36卷第9期，2015年。
3 Petre Stoica and Randolph L. Moses. 信号频谱分析 . Pearson Prentice Hall，2005年
4Steven W. Smith。 面向科学家和工程师的数字信号处理指南 。California Technical Publishing，1997年。
致谢
作者感谢ADI公司的Sefa Demirtas、Bob Adams和Tony Akl对此算法的开发所提供的帮助和宝贵意见。
作者
Foroohar Foroozan是信号处理科学家。Foroohar于2015年8月加入ADI公司，领导医疗健康业务部面向生命体征和家用监测系统的多伦多算法团队。同时，她就护理点超声成像系统与ADI成像团队展开合作。加入ADI之前，她担任加拿大Geotech Ltd.的研发科学家，致力于新一代机载电磁地球物理测量系统的智能滤波。她是Sunnybrook Research Institute博士后，在2012年至2013年期间致力于脑血管病图的3D、超高分辨率超声成像，正在申请PCT专利。她于2011年获得加拿大多伦多约克大学-拉松德工学院计算机科学博士学位。她对生物医学系统中的信号处理和算法感兴趣，主要致力于生命体征系统和生物医学成像。她是安省专业工程师协会(P.Eng.)成员和IEEE高级成员。