基于相位处理的颈动脉脉搏波提取方法

高冠群, 杨学志, 刘雪南, 张 龙, 王定良

(1.合肥工业大学 计算机与信息学院,安徽 合肥 230601; 2.工业安全与应急技术安徽省重点实验室,安徽 合肥 230601)

0 引 言

心脏作为人体心血管系统的发动机,负责给整个心血管系统提供动力,促进血液在全身周而复始地流动,在皮肤表面引发颜色变化的同时对血管壁引发压力变化,进而形成搏动的特征[1],脉搏波便是记录这些波动参数的曲线,同时也反映了心血管系统的健康状况[2],是评价健康水平的重要指标。

血液在整个心血管系统中的流动,会使皮肤表面产生2种变化:颜色变化[3]和运动变化[4]。颜色变化是由于血液容积的变化导致对自然光的吸收率不同,进而引发皮肤周期性的明暗变化。运动变化则是由于血液在流经各级动脉系统时会对各级动脉血管壁产生不同的压力波,压力波随血液传播在皮肤表面引发皮肤的周期性微小运动。通过图像光电容积脉搏波描记法[5](image photo plethysmography,IPPG),可以利用远距离非接触式拍摄的视频获取脉搏波信号。

文献[6]使用人脸视频作为输入进行脉搏波检测,并成功计算心率。该方法在人脸中划分感兴趣区域(region of interest,ROI),在ROI区域中通过计算像素值的变化来表征面部的颜色变化,通过独立成分分析(independent component analysis,ICA)模型[7]从R、G、B 3个信号中分离出3个独立的基向量并选用第1个向量估算心率。文献[8]改进了上述算法,对每一个向量进行功率谱分析并选用最大的1个向量作为血液容积脉搏波(blood volumn pulse,BVP)信号。

文献[9]基于欧拉视频放大(Eulerian video magnification,EVM)算法[10]提出了一种非接触式脉搏波检测方法,将脸部颜色即血液引起的变化进行放大,而后取全部像素值进行平均得出BVP,统计该BVP波形的极值点计算心率,该方法仅使用像素平均值进行计算,受噪声干扰严重,准确率低。

文献[11]提出了根据头部运动提取脉搏波并计算心率的算法,由于血液流经颈动脉到达头部,促使头部产生小幅度运动。算法通过跟踪人头部的特征点,并根据感兴趣的时域频带进行过滤,而后使用主成分分析法(principal component analysis,PCA)[12]提取由脉冲引起的周期性信号。最后通过检测信号的频谱来提取平均脉冲率,并使用峰值检测算法获得精确的心率,该方法受头部微小晃动噪声影响较大,影响检测精度。

文献[13]提出了一种抗干扰的心率检测算法,该算法针对运动干扰问题,使用检测人脸特征点以及角点跟踪(kanade lucas tomasi tracking,KLT)算法[14]跟踪人脸并实施矫正。该算法提高了在理想情况下的检测精度,但是因使用理想带通滤波损失了很多检测细节,并且计算量庞大。

文献[15]基于EVM算法,提出了对放大腕部桡动脉微小变化提取脉搏波的算法,桡动脉的搏动变化可以比面部颜色的变化提供更多的信息,波形更加完备,但是该算法使用颜色空间进行分割和选取,抗光照干扰的能力较差,成功率较低。

针对上述问题,本文提出一种基于相位变化的颈动脉脉搏波提取算法,利用相位变化和运动变化的对应关系,建立基于相位变化[16]的颈动脉搏动信号的描述实现对脉搏波的检测,不仅有效保存了波形的细节信息,而且提高了检测精度和抗噪性能。

1 基于相位的脉搏波提取方法

本文提出了一种基于相位变化的脉搏波提取方法,算法主要分为以下4个部分:① 使用EVM进行视频放大,并将放大后的视频分块计算信噪比(signal noise ratio,SNR),选择信噪比最高的分块,即为ROI区域;② 对ROI区域进行复可控金字塔分解[17],得到局部频谱和相位谱,对相位谱进行计算得到相位差;③ 通过时域带通滤波选择感兴趣频段内的相位差信号,并使用主成分分析算法[12]提取主要成分;④ 通过信号选择分离出纯净的脉搏波信号,最后估算心率。

整体算法框架如图1所示。

图1 整体算法框架

1.1 选择最优ROI区域

选择最优ROI区域的方法是通过EVM算法放大颈动脉的微小运动,然后对每幅图像分块,计算每块的信噪比,取信噪比最大的块作为感兴趣区域。

输入的视频每一帧的大小为640×480,分块后选择子块进行计算可以滤除环境噪声,同时节省运算时间和降低计算量,并降低 PCA提取主成分的难度。

EVM在空间域上利用拉普拉斯金字塔将视频分解为不同空间分辨率的帧序列;然后对视频帧序列使用快速傅里叶变换(fast Fourier transform,FFT)变换转换到频域进行滤波,使用有限脉冲响应(infinite impulse response,IIR)滤波器滤除不相关频段的信号,而后进行视频放大和重构。

将放大后的视频重新输入,并将每一帧图分成大小为40×40的块,共计192块,均匀分割整幅图像的同时保证每个分块中包含足够的颈动脉搏动信息。计算每个子块的信噪比,选出信噪比最大的块作为最优ROI区域。

本节输入视频为640×480×600矩阵,令输入矩阵表示为G(i,j,t)。其中:t为视频的总帧数;i、j分别为在第t帧内像素点的横坐标和纵坐标。选取最大信噪比过程如图2所示。

图2 最优ROI区域选取示意图

选取最优ROI的具体步骤如下:

(1) 对图像分块,沿x轴与y轴方向,以40为边长进行分块,每一帧图像分为192个块,令Gmn(s,k,t)表示分块函数。其中:m、n分别为每一块的坐标位置；s、k、t为分块内部的像素位置和时间序列。

(2) 求第(m,n)个分块的像素平均值,计算公式如下:

(1)

其中,ROImn(t)为第(m,n)个分块第t帧的像素平均值,是1个1×600的时间序列。

(3) 求每个子块的信噪比,并求出信噪比最大的分块坐标,计算公式如下:

(2)

其中,(u,v)为所求最大信噪比分块所在位置的坐标。

1.2 复可控金字塔分解与计算相位差

金字塔分解及计算相位差流程如图3所示。

图3 金字塔分解与计算相位差流程

对视频中ROI区域的视频帧序列做处理。使用复可控金字塔,对图像进行不同尺度和不同方向的分解。不同方向滤波器的设计公式如下:

(3)

(3)式中,k表示分解的总方向个数,每个方向b下的方向模板是anglemaskb,b取值范围为1，2，3，…，k。例如,当k=4时,可以得到4个方向的滤波器,即使用复可控金字塔对4个方向的运动进行分解。由于颈动脉是1条垂直方向的大动脉,视频中计算公式可以被捕捉到的运动是水平方向运动,因此选择水平方向相位谱。

假设视频输入信号在时刻t,以视频中ROI区域内的1个像素点为例,位置x处的像素,一维信号的强度表示为l(x,t)。随着时间推移,像素的位置产生移动后的强度计算公式如下:

(4)

其中,δ(t)为位移函数中较小的幅度。首先,利用傅里叶级数分解,将位移图像轮廓f(x+δ(t))表示成多个复杂的正弦曲线之和,计算公式如下[16]:

(5)

其中:每一个子带都对应着1个频率ω；φω为相位信号。

根据FFT变换的平移性质,可以计算出I(x,t)的傅里叶变换结果，公式如下:

I(x,t)=f(x-δ(t))=

(6)

对于(6)式中的每一个频率ω,其复指数形式如下:

Sω(x,t)=Aωei(φω+ωδ(t))e-iωx

(7)

当t=0时,相位为φω,其他时刻的相位为φω=φδ(t)。将t≠0时刻的相位与t=0时刻的相位相减,得相位差公式如下:

Bω(x,t)=ωδ(t)

(8)

则点x处的相位差是1个1×600的一维信号,600是视频帧的个数。这里需要计算ROI分块中所有像素点的相位差信号,然后进行下一步的时域滤波和PCA提取主成分。

1.3 PCA提取主成分与计算心率

得到相位差值后,进行时域带通滤波,通过设置合适的带通上下阈值,提取到感兴趣频段内的运动,这里使用经典的带通滤波器IIR带通滤波器完成滤波。

将每一个点的相位差信号作为1个样本数据,组成1个输入数据矩阵,即

mt=[B1(t)B2(t) …BN(t)]

(9)

其中：t为视频帧数；N为ROI分块中的像素个数；N=40×40,共计1 600个像素点。

求输入矩阵的协方差矩阵，即

(10)

(11)

协方差矩阵的特征向量和特征值为:

(12)

其中:Λm为对角矩阵,对应求出的特征为λ1,λ2,…,λN;Φm为对应的特征向量;Φ1,Φ2,…,ΦN组成特征向量矩阵。

将特征向量按特征值的大小进行从上到下的排列,取前D个特征值所对应的特征向量,并组成矩阵QL。这里,D个特征值对应着D个主要成分,其中参数D是通过计算主成分的贡献率从而得到的。D个主成分一一对应的信号都可以通过对特征向量矩阵进行反投影得到,公式如下:

(13)

PCA方法提取了D个主要成分,但还要从中找出最合适的1个主成分代表脉搏波。虽然主成分按照方差大小依次排列,第1个主成分解释了大部分方差,但S1可能无法最清晰地分析脉搏波,应该选择周期性最强的Si。实验发现前5个主成分包含了原始信号中98%以上的信息,因此只在前5个主成分中选择即可[11],根据功率谱的最大值占总功率谱值的百分比来计算,选择前5个中最大的1个,即为提取的脉搏波。所提取的脉搏波如图4所示。

图4 脉搏波信号示意图

提取到脉搏波之后,通过FFT变换得到脉搏波对应的功率谱,并通过功率谱计算心率值,功率谱中功率最大的值所对应的频率即为心跳的频率。估算心率流程如下:

输入总帧数为N=600的脉搏波信号y以及功率谱P(f),计算功率谱中峰值所对应的频率fmax。

(14)

根据(14)式求得的最大频率,计算心率值。

HRy=60fmax

(15)

(15)式中,HRy为心率值,通过计算得到的心率值可以作为脉搏波提取算法的定量分析的重要参考指标。

2 实验结果与分析

本节将设计2组对比试验验证该算法的有效性和准确性。

(1) 使用接触式检测设备PWS-20D的值作为真值,与算法的实验值对比,验证本文提出算法的准确性和稳定性。

(2) 在不同光照条件下,以心率值作为定量分析的标准,分别与文献[13]和文献[15]提出的非接触式脉搏波检测算法进行对比,验证本文提出算法在准确度和抗光照干扰方面的优势。

2.1 实验设计

实验装置如图5所示。使用笔记本电脑自带的摄像头作为采集设备,PWS-20D设备也同时连接在电脑上,笔记本摄像头正对受测者颈部,正对距离大概0.5 m左右。采集视频时应保证受测者上半身(尤其是颈部和头部)避免较大幅度的运动。在采集颈部视频的同时,也使用PWS-20D设备通过接触式方法采集指尖的脉搏波,该设备基于PPG原理,对指尖发出强光并通过传感器接收信号,可以准确采集脉搏波波形并自动保存成一维数组,能够作为视频处理所得脉搏波的参考。

图5 实验装置示意图

实验过程通过surface笔记本电脑控制,程序运行环境是Matlab,视频录制采用RGB颜色空间,视频帧率为30帧/s,设置视频分辨率为640×480,摄制时长为20 s,总计600帧。本次实验中共15名受测者参与实验,包括10名男性和5名女性,年龄为22～68岁,每次采集均是同步进行颈部视频的录制和PWS-20D设备的波形采集。最后,从定量分析的角度来评价算法的有效性,采用平均误差Me、均方误差RMSE和准确率HRac作为本节的评价指标。

2.2 实验结果分析

第1组实验,为了验证所提出算法的准确性和有效性,通过本文算法计算得出的心率实验值和PWS-20D设备的心率真值进行对比。对15名志愿者进行颈部颈总动脉脉搏波的提取,同步使用非接触式脉搏波检测设备PWS-20D采集脉搏波,计算心率值。

为了验证算法的准确性,对通过多组共同测量的数据得到的散点图进行评估。实验中,分别对15名受测者每人做3组实验,每组实验均包含接触式检测和非接触式检测,共45对数据。受测者包括术后康复期的病人(术后康复期的病人心率一般较低),以及静息状态下的年轻人和轻微运动后的年轻人。使用散点图进行一致性分析，如图6所示。

图6 一元线性回归分析检测结果的一致性

图6中横坐标是心率参考值,纵坐标是心率实验值,每组输入数据包含1组心率实验值和真实值的数据,对应了坐标系中的1个点。参考值的点会呈现出斜率为1的直线,本实验中参考值如图6中的蓝线,横贯坐标轴,实验点则以(真实值,实验值)的形式分布在直线附近。图6中星号所代表的每一个点的横坐标即为接触式设备PWS-20D的值,纵坐标为通过本文方法测量出来的实验值。

从图6可以看出,所有的点都分布在直线附近,表明本文提出的非接触式脉搏波检测算法具有很好的准确性和一致性。

为了验证本文算法抗光照干扰的能力,设计第2组实验,在自然光、弱光和LED灯3种不同光源下,将本文算法分别与文献[13]和文献[15]进行对比。

文献[13]使用面部作为检测部位,通过颜色放大脸部毛细血管中血液颜色的微小变化提取脉搏波;文献[15]使用腕部作为检测部位,通过欧拉运动放大的方法对腕部桡动脉进行微小运动的放大提取脉搏波。

本次实验输入视频时长为20 s,帧率为30帧/s。不同环境光下3种方法的心率检测结果见表1所列。实验结果对比如图7所示。

图7 实验结果对比图

表1 不同光照场景下心率检测结果

实验结果表明,在3种光源下本文方法均具有较好的准确性和稳定性,且在弱光场景下,本文算法的准确率依旧达到91.45%,平均误差只有3.208,具有明显优势。

为了验证本文算法抗运动干扰的能力,设计第3组实验,实验中共10名受测者,在自然光背景下,5名受测者进行受检测部位的小幅运动(包括说话和点头),3名受测者保持不动,2名受测者进行受测部位的大运动(包括转动脖子,左右晃动)。录制共计30个视频,将本文算法分别和文献[13]及文献[15]进行运动干扰的对比,结果见表2所列。

表2 运动干扰下心率检测结果

从表1可以看出,本文方法各项指标比文献[13]和文献[15]的方法在不同光照干扰下均有显著提升。

其主要原因如下:

(1) 文献[13]使用脸部视频,通过颜色放大检测表面毛细血管中血液的明暗变化,以此检测脉搏波并计算心率。该方法将视频由RGB通道转到YIQ通道进行计算,虽然能够较好地观测亮度变化,但是大大降低了抗环境光亮度变化干扰的能力。利用复可控金字塔进行方向分解并选择水平运动方向的相位谱计算相位差,同时,通过相位差进行计算在原理上便可以较好地减少光照干扰。

(2) 文献[15]采用腕部桡动脉视频作为输入信号,但是该方法将整个手臂的区域整体进行计算得出脉搏波,易受到环境光亮度变化的干扰,会受到手臂运动噪声的干扰,导致检测结果不够准确。

本文方法通过视频放大和计算信噪比选择颈动脉搏动区域作为ROI区域,滤除大量环境噪声,因此比文献[15]的准确度高。

从表2可以看出,本文方法各项指标较文献[13]和文献[15]的方法在运动干扰下也有显著提升。其主要原因如下:

(1) 文献[13]通过旋转矫正解决面部运动干扰,但矫正会减少数据量,同时矫正对于头部转动等运动效果很差。

(2) 文献[15]对检测手臂整体区域进行处理,容易受到环境运动噪声的干扰从而影响准确率。

本文方法首先通过ROI选择感兴趣区域来滤除大量环境噪声的干扰,同时利用复可控金字塔进行方向分解并选择水平运动方向的相位谱计算相位差,可以直接滤除水平方向之外其他方向的运动噪声,因此比文献[13]和文献[15]的准确率高。

3 结 论

本文提出了一种基于相位处理的脉搏波提取算法,该算法基于IPPG原理,通过录制颈部颈总动脉的搏动视频,首先选取最优 ROI区域,然后对ROI区域进行复可控金字塔分解,取水平方向相位谱计算相位差,然后进行时域滤波,通过PCA算法提取主成分,即为所求脉搏波信息,使得处理后的波形更加纯净且具有更多的细节信息。本文设计了3组实验对算法准确性和稳定性进行验证,通过与PWS-20D设备进行准确率对比;与文献[13]和文献[15]所提的方法在不同光照条件下进行对比实验,验证了本文算法具有较高的准确度和较强的抗光噪声干扰能力。