一种基于贝叶斯估计的脉搏波峰值检测算法

李思楠,赵 海

1(河北民族师范学院 物理与电子工程系,河北 承德 067000)2(东北大学 计算机科学与工程学院,沈阳 110000)

1 引 言

脉搏信号与人体心血管系统有着紧密的联系,脉搏波主波峰的位置对应着心跳发生的时刻,往往通过定义相邻脉搏波峰之间的间隔来划分完整的心跳周期,进而计算出人体的心率.目前主流的体域网产品,如智能手环、手表等也都具备了采集和记录脉搏波信号的功能.利用光电容积原理[1]在可穿戴产品中采集脉搏波信号来提取峰值点以计算心率等心脏相关的参数是一种更为简单和低成本的方式.同时为使用者瞬时心率的记录和评价提供了更为便利的条件.目前针对脉搏波峰值点提取的研究中,Aniruddha等人在2007年提出一种基于小波变换的脉搏波峰值提取方法[2];Schmidit等人在2010年提出了将脉搏波进行高通滤波处理后并进行归一化阈值筛选的方法来确定脉搏波中的峰值点和波谷点[3],Pachauri等人在2011年研究使用滑动时间窗来计算脉搏波的能量包络,进而自动划分出脉搏波中的特征点[4];Yazdani等人在2018年提出了一种利用长短时间窗结合阈值设定的方式来进行峰值点提取的方法[5].以上方法是在平稳的脉搏波状态下进行的研究,且都需要人为设定阈值进行筛选,未考虑其他状态下的检测过程.

然而,可穿戴产品采集到的脉搏波信号不仅会掺杂周围环境引入的噪声,还会因为人体定向的运动,导致信号中叠加了运动干扰,导致脉搏波形可能会偏离了正常的形态结构.如何从带有噪声和运动干扰的脉搏波信号中准确识别特征点已成为脉搏波研究领域的一个重要研究方向.Raghuram等人提出用小波分解的方法来分析带有运动干扰的脉搏波信号[6],不过小波基函数的选择会对结果产生较大影响,需要根据经验进行选择;孙旭月等人采用经验模态分解的方法来对脉搏波信号进行分析处理[7],但该方法未直接提出峰值点的检测过程;Fukushima等人提出一种以加速度传感器采集的信号作为参考来消除脉搏波信号中运动干扰并提取特征点的方法[8];Yousefi等人采集双通道脉搏波信号,以其中一路作为参考信号,提出了归一化最小方差算法来计算脉搏波特征点[9];以上两种方法需要额外的一路信号作为参考,增加了检测成本;Scholkmann等人针对于带有运动干扰的准周期信号提出了一种自动多尺度峰值检测(AMPD)算法[10],通过计算信号局部最大尺度值的方式来对峰值点进行检测,该方法在心电信号峰值检测中取得了良好效果,但复杂形态下的脉搏波信号主波峰不一定是局部最大值点,该算法会造成峰值点被遗漏.

本文受到文献[10]中方法的启发,提出一种基于贝叶斯估计的脉搏波主波峰值点检测改进算法(B-AMPD),从时域提取出脉搏信号的峰值点(本文提到的峰值点检测均指主波峰值点的检测),该方法无需额外参考信号的引入,利用脉搏波峰峰值间隔的先验概率分布[11]来预测峰值点位置出现的概率,这样降低了传统AMPD算法遗漏峰值点的缺陷造成对检测结果的影响,通过MIT-BIH标准数据库的验证,在灵敏度和阳性预测值性能上比传统的AMPD算法都得到了提升.而且对于掺杂了运动干扰的脉搏波信号,本文的算法在灵敏度检测上取得了更好的效果,大大降低了峰值点被漏判的情况,为脉搏波的周期合理划分和瞬时心率的准确计算提供了基础.

2 基于多尺度峰值检测的脉搏波峰值提取

2.1 自动多尺度峰值检测(AMPD)算法

AMPD算法的目的是通过分析信号的局部最大尺度图来检测信号的峰值.用X=[x1,x2,…,xN]表示采集得到的N个采样点的脉搏波信号,AMPD算法通过一个滑动窗来计算该信号的局部最大尺度值,滑动窗的长度wk(wk=2k|k=1,2,…,L)是一个不断变化的量,其中k被定义为信号的分析尺度,L=[N/2]-1(符号[m]表示不大于m的最小整数).滑动窗大小随着不同的尺度k值而变化,以不同的分辨率来覆盖整个脉搏波信号.通过在信号中每个尺度k处执行公式(1)的标准来搜索局部最大值,对于每个k值,当i=k+2,…,N-k+1时有:

(1)

由此得到一个矩阵M,即:

(2)

其中第k行就包含了对应的滑动窗的长度wk,这个矩阵M中1的位置就是在每个尺度k下的局部最大值(潜在的峰值点的位置),如果对于每个尺度k,矩阵M中某一列i的所有元素值都为1,则这一列所对应的采样点i的位置即是峰值点.由此,可以计算出M中所有列元素全为1的列标i,进而得出这组信号X的所有峰值点位置.

2.2 基于拉普拉斯概率模型的峰值点概率估计

根据AMPD算法,矩阵M每一列中元素1的个数所占的比例,即是该采样点被判断为峰值点的概率,可以表示为:

(3)

从脉搏波信号中计算瞬时心率,首先要利用AMPD算法来检测每个波形周期内的峰值点位置,如图1所示,第j个峰值点所对应时间用tj表示,将相邻的两个峰值点之间的时间间隔定义为:

rj=tj+1-tj

(4)

图1 峰值点时间位置和峰峰值间隔rFig.1 Instants at which a peak is detected,and r,the intervals between peaks(RR intervals)

文献[11]中研究了r中的元素并不是独立且相同分布的,它们很大程度上取决于先验的分布.并且定义了另外一种表示脉搏波信号峰值的方法是定义一个距离d,用它来表示连续记录的峰峰值间隔之间的差异,即:

dj=rj+1-rj

(5)

根据文献[11]中的研究结果,拉普拉斯分布模型可以用来分析d,而采样点i成为峰值点的概率p2可以通过拉普拉斯概率密度分布来计算:

(6)

其中,di跟公式(5)中表示的内容有所不同,它表示为di=ri-ri-1,ri表示为采样点i与其前一个峰值点之间的间隔,而ri-1表示采样点i的上一个峰峰值间隔时间.u是拉普拉斯密度函数中的位置参数,它是di的中位数值,而b作为密度函数中的尺度参数,其计算公式为:

(7)

其中K表示d中元素的个数.计算出p1和p2后,可以通过比较对应位置p1和p2的大小,如公式(8),比较二维数组[p1(i);p2(i)],选取对应列中的最大值位置来进一步确定峰值点的位置,即:

(8)

3 基于贝叶斯模型的峰值点检测

3.1 基于Beta分布模型的峰值点概率估计

由于人的瞬时心率是在不断的变化的动态过程,为了进一步增强峰值点判断的准确性,需要构建一个动态可更新的概率估计模型.本文采用贝叶斯原理来构建此过程,首先,假设任何一个采样点i成为峰值点的概率表示为θi(θi初始值可以通过p2(i)=θi来计算).根据贝叶斯公式,可以推断后验分布概率为:

posteriori=P(θi|AMPDoutput(i))
=P(AMPDoutput(i)|θi)·P(θi)prior

(9)

其中,P(θi)prior是θi的先验概率分布,P(AMPDoutput(i)|θi)代表θi已定条件下的通过AMPD算法输出的采样点i为峰值点的似然估计.同时,需要为P(θi)prior构建一个分布模型,并使其在不断学习过程中进行更新.这里选用Beta分布来构建此模型,如图2所示,为利用Beta方程构建的P(θi)prior的概率分布模型,它是一个偏正态的概率模型,当θi很小的时候(采样点i靠近上一个峰值点)曲线偏向左边;当θi较大时(采样点i靠近下一个峰值点)曲线偏向右边;而采样点靠近两个峰值点中间则曲线更趋近于均匀分布.Beta分布的这种统计学特性可以为P(θi)prior建模提供良好的效果.

(10)

图2 根据Beta方程构建P(θi)prior先验概率分布模型Fig.2 Modelling for the prior probability distribution(P(θi)prior)as Beta function

其中,Γ表示伽玛方程,Γ(α+β)/Γ(α)Γ(β)是为了使得所有概率总和为1的归一化常数.α和β被称为控制分布形状的超参数,需要通过设置α和β这两个参数值来准确的表征出每一个采样点i的先验概率分布.Beta分布的期望值可以表示为:

E[θi-prior]=α/(α+β)

(11)

对于似然估计P(AMPDoutput(i)|θi),如果将AMPDoutput(i)看作是利用AMPD算法计算得出的一组只包含0或1的二进制数组(公式(3)中概率p1(i)为1则AMPDoutput(i)为1,p1(i)不为1则AMPDoutput(i)为0).假设每个采样点的峰值概率θi是独立同分布的,则通过AMPD算法输出的二进制数组是符合伯努利分布的形式.因此,AMPDoutput(i)的似然函数可以表示为:

(12)

将公式(10)和公式(12)带入公式(9),可以得到一个θi的后验概率分布是一个有全新参数的Beta分布形式:B(α′,β′),表示为：

(posteriori=B(α′,β′)=θn+α-1(1-θi)β-n

(13)

后验概率posteriori在每个脉搏波波型周期内进行更新,如果AMPD算法对采样点i的判别输出为1(判别为峰值点),即n=1,则α′=α+1,β′=β.同理,AMPD算法对采样点i的判别输出为0(判别不是峰值点),即n=0,则α′=α+1,β′=β+1.

3.2 贝叶斯学习过程

如果将连续两个峰峰值之间的间隔定义为一个脉搏波周期,从第一个脉搏周期开始,周期内的采样点的先验概率θi通过公式(6)来初始化,而到了下一个周期,θi的后验概率分布通过α和β的参数更新利用公式(13)进行计算.实际学习过程如图3所示,α′=α+1导致了Beta曲线向右移动(增加了采样点i成为峰值的置信度),而β′=β+1导致Beta曲线向左移动(降低了该采样点i成为峰值的置信度).这样的学习过程在每个周期内重复,因为那些最可能成为峰值点的采样点i的后验概率分布曲线向右倾斜,它们将具有较高的期望值.

由此,采样点i被判定为峰值点的置信度可以通过θi的后验概率来进行计算:

(14)

图3 学习过程示意图Fig.3 Diagram of learning procedure

其中阈值threshold通过经验进行设置,本文中设置为0.9.最后,将p3(i)的值与公式(3)计算出的p1(i)进行比较,取二者中的较大值对应的采样点i即是最终认定的峰值点,即:

(15)

4 实验测试和评价

4.1 样本选择

为了检验此算法的性能,本文选择15组MIT-BIH标准数据库里的信号进行分析,它是由美国的麻省理工学院提供的包含有多种人体生理信号的公开数据库,其中有9组数据是从研究慢性阻塞性睡眠呼吸暂停综合症实验中获得,实验者年龄32～56岁,实验者采集过程中处于仰卧放松或睡眠状态下,采样率250Hz,与脉搏波同步采集的还有人体的心电(ECG)、脑电(EEG)和呼吸信号,累计记录时长超过了8小时.其余6组数据源自运动中的脉搏波采集实验,实验者年龄22～32岁,实验者处于健身房中的跑步机上进行采集,可以自由设定跑步机的速度来调整自身步速,采样率256Hz,每组采样过程持续4～6分钟,同步采集的还有人体的心电(ECG)信号.

为了便于统计分析,本文选取时间长度为每组3分钟,并且选择同步采集到的心电(ECG)信号为检测标准,MIT-BIH数据中已经给出了心电信号ECG的检测结果,可以跟据心电信号中每一个R波峰都在脉搏波信号中有一个主波峰相对应的关系,将心电信号的检测结果作为脉搏波峰值点检测算法的正确性检验指标.

4.2 验证方法与评价指标

实验软件环境为Matlab2010b,为了评价所提出的算法的好坏,对该算法进行了两方面的验证.首先,从该算法对给定信号干扰强度下的峰值点检测进行评价,从数据库中选择一组平稳光滑的脉搏波信号,对其分别加入不同信噪比(SNR)强度的干扰,研究其峰值点检测的性能随着引入干扰强度的变化趋势.其次,对比该算法与AMPD算法的检测效果,选择不同人体生理状态下的脉搏波信号,分别利用两种算法进行峰值检测,对性能指标进行对比,观察该算法的提升效果.

对于脉搏波峰值点提取算法的检测效果,可以通过灵敏度(Sensitivity,Se)、阳性预测值(Positive Predictive Valve,PPV)和检错率(Detection Error Rate,DER)[12,13]三个指标进行评价.它们的定义分别为:

(16)

(17)

(18)

其中,TP表示真实被检测到的峰值点,FP表示检测到误判为真的虚假峰值点,FN代表被遗漏的真实峰值点,TN表示成功被排除的虚假峰值点.由此可知,灵敏度(Se)就表示了在所有峰值点中利用算法检测出的真实峰值点的比例;而阳性预测值(PPV)表征被检测出的峰值点中真实的峰值点所占的比例;检错率代(DER)表了检测过程中虚假峰值点被误判为真和被漏掉的真实峰值点占所有峰值点的比例.

5 结果分析与结论

图4所示为MIT-BIH 数据库中选择slp01a记录的时长为3分钟的平稳光滑的脉搏波信号和其参考系(心电信号),以心电信号中峰值点个数为标准,向脉搏波信号中加入均值为0.3、信噪比强度为分别为10、15、20、25、30(dB)的泊松分布随机数来模拟信号受到不同程度的运动干扰,分别统计利用本文算法检测到的峰值点的评价指标.

图4 平稳的脉搏波信号Fig.4 Smooth pulse wave signals

图5 算法的灵敏度、阳性预测值和检错率三项指标随信噪比变化的评价效果图Fig.5 Performance evaluation:Positive Predictive Valve,Sensitivity and Detection Error Rate wtih the change of SNR

图5中横坐标表示信噪比强度,从右向左,信噪比下降,也代表了加入干扰的增强.左侧纵坐标表征了两项性能指标Se(实心三角符号表示)和PPV(实心方格符号表示)的变化趋势,从右向左随着干扰的增加,Se的值有明显的下降,而PPV值下降幅度不大,这说明在已检测出的峰值点中真实的峰值点比例很高.而右侧纵坐标反映了DER(实心圆点表示)的变化趋势,从右向左随着干扰加强,DER也明显上升,这说明DER上升的主要原因还是因为一些峰值点被遗漏导致.

表1 文本算法在MIT-BIH数据库识别峰值点效果
Table 1 Performance of the proposed algorithm in the detection of pulse wave peaks

记录名称Se(%)PPV(%)DER(%)slp02a98.0899.512.40slp02b10099.500.50slp0310099.060.95slp0499.5498.202.29slp1498.6299.082.36slp1610098.151.95slp3298.5299.012.46slp4198.0597.574.39s1_walk97.5495.487.07s2_walk98.3697.094.59s3_walk98.2698.263.47s3_run97.2897.585.12s4_run95.2497.777.30s5_run94.7597.087.58平均值98.1798.153.69

图6 几种生理状态下的脉搏波峰值点检测Fig.6 Pulse wave peak detection in several physiological states

可穿戴产品在实际的脉搏采集过程中,运动干扰引入的大小往往是难以准确衡量的,为了评估本文提出的算法在不同场景下的检测效果,选择了剩余的14组脉搏波数据进行实验.统计结果如表1所示,这14组实验中人体生理状态包含了睡眠、步行和跑步.平均的灵敏度(Se)检测能达到98.17%,阳性预测值(PPV)能达到98.15%,检错率(DER)3.69%.从统计结果可以看出,在睡眠状态下,本文提出的算法具备非常高的Se和PPV,当随着生理状态的改变,运动幅度的加大对本文峰值点的检测产生了影响,出现了峰值点遗漏和误判的情况.

图6中,a表示睡眠状态下脉搏波舒张期幅值跃升的情况,本文的改进算法可以将峰值点进行准确识别;b表示睡眠状态下脉搏波舒张期出现纹波的情况,主波峰形态未发生改变,算法依然可以将峰值点识别;c表示步行状态下的检测,步行状态下脉搏波的波形结构发生了形变,个别舒张期峰值高于主波峰幅度,就如c中左侧所示,真正的峰值点被遗漏而将虚假峰值点误判为主波峰;d表示的为人体处于跑步状态下情况,脉搏波形态发生了明显形变,有些波形的舒张期消失,仅保留了主波峰反映了左心室向主动脉射血的程度[14],本文的算法可以将保留有明显主波形状的峰值点进行识别.

表2 本文提出的算法和AMPD算法在Se、PPV
和DER三项性能指标上的对比
Table 2 Comparison of the proposed algorithm and AMPD algorithm on the performance of Se,PPV and DER

生理状态算法名称平均Se平均PPV平均DER睡眠状态AMPD算法96.35%97.28%4.37%本文算法99.10%98.76%2.16%步行状态AMPD算法87.26%90.08%15.34%本文算法98.05%96.94%5.04%跑步状态AMPD算法80.61%88.37%20.74%本文算法95.76%97.47%6.67%

表2显示本文提出的算法和AMPD算法在不同状态下的性能对比.在睡眠状态下,脉搏波的信号相对平稳,本文的算法只是在Se值上比AMPD算法有明显提升,在PPV值上两者十分接近,说明在平稳状态下传统的AMPD算法也很少将虚假峰值点误判为真实峰值点,性能的落后在于出现漏判峰值点的情况.在步行状态下,运动干扰的加入使得脉搏波波形发生了改变,AMPD算法和本文的算法性能指标都有不同程度的下降,运动干扰的增加使得AMPD算法的Se值明显下降,波形的形变会影响AMPD算法对局部最大值的判断,造成大量峰值点检测的遗漏,这样使得DER明显下降.而本文提出的算法根据了贝叶斯估计的理论对峰值点的漏判进行了较大的弥补,在跑步状态下,脉搏波中的干扰较步行状态下进一步增强,脉搏波的舒张期基本消失,只留下了主波峰的基本形态结构,AMPD算法的检测的Se和PPV值较睡眠状态下分别下降了15.74%和8.91%,DER较睡眠状态下提高了16.37%,而用本文算法计算出的Se和PPV较睡眠状态下分别下降了3.34%和1.29%.这表明在三种生理状态下(睡眠、步行、跑步)改进算法在性能上明显优于AMPD算法.比较步行和跑步两种状态下AMPD算法的PPV值下降幅度远低于Se,说明运动干扰对AMPD算法的主要影响在于使得峰值点遗漏数量增加,改进算法则在此性能上有明显的优势.

从文本进行的实验可以看出,当脉搏波中掺杂的干扰是定向的而非杂乱无章的干扰,未使脉搏波主波形态结构发生根本性改变的前提下,提出的改进算法比AMPD算法在检测脉搏波峰值点性能Se、PPV和DER上都有所提升,在Se上提升最为显著.

6 结束语

本文提出了一种基于贝叶斯估计的改进脉搏波峰值点检测算法,在传统AMPD算法的基础上引入贝叶斯学习理论,对脉搏波主波峰值点检测取得了良好的效果,面对不同生理状态对脉搏波引入的运动干扰,利用MIT-BIH的标准数据库作为检测标准,依然在平均灵敏度(Se)和阳性预测值(PPV)上取得了98%以上的效果.同时,从本文中的实验也可以发现,随着运动干扰的加大,当脉搏波基本形态结构发生明显的形变时,本文的算法也依然会出现逐渐增多的漏判和误判,这对算法的灵敏度影响较大,但比AMPD算法有明显的提升.后续的研究中可以考虑波形在形变状态下对于峰值点位置的判别和检测准确率的提升方法.