基于运动特征预估的呼吸率检测快速算法

朱明扬,陈 鲸,杨学志,沈 晶,吴克伟

(1.合肥工业大学 计算机与信息学院,安徽 合肥 230601; 2.工业安全与应急技术安徽省重点实验室,安徽 合肥 230601; 3.合肥工业大学 软件学院,安徽 合肥 230601; 4.合肥师范学院 电子信息与电气工程学院,安徽 合肥 230061; 5.安徽微波与通信工程技术研究中心,安徽 合肥 230061)

随着计算机科学技术的进步,视频图像处理技术快速发展,并与医疗领域相结合,在医学诊断和人体健康监测等方面成为现代医学的重要助力。在医学领域,呼吸率是重要的健康指标,对呼吸率的监测在呼吸疾病的预防方面具有重大意义。传统的呼吸率检测多为接触式方法,通过传感器探测人体因呼吸运动引起的变化来获取呼吸率,该类方法需要在被测者身上粘贴传感器,会引起诸多不便,且对于长期癫痫患者、烧伤患者、睡眠中患者难以实施。基于视频的呼吸率检测方法是一种远距离、非接触式的生理测量方法,可以避免带来患者身体上的不适,并提供长时间、自动化的呼吸率检测,相比于接触式检测方法更具优势,是未来呼吸率检测领域的发展方向。

近年来,基于视频的呼吸率检测方法得到了广泛的研究。其中,从测量方法上可分为直接测量法和间接测量法。从心电信号、光电容积脉搏波或动脉血压信号等信号中提取呼吸信号属于间接测量。文献[1]通过小波分析和经验模态分解从光电容积脉搏波中提取呼吸率；文献[2]通过连续小波变换从脉搏波信号中提取呼吸率,这些方法都是从已有的完整信号进行再研究,实时性不强。基于机器视觉记录胸部起伏的呼吸率检测方法属于直接测量。文献[3]提出了基于最大似然的呼吸率检测方法,但该方法没有对呼吸区域进行预估,对全图进行最大似然估计定位呼吸区域,因此算法时间复杂度较高;文献[4]提出了基于相位的呼吸率检测算法,该算法具有良好的抗噪性能,但是算法定位的呼吸区域过大,且后续在使用复可控金字塔分解获取呼吸信号时没有选择呼吸运动方向,影响了算法整体速度;文献[5]利用光流法加欧拉算法进行呼吸率检测,实现对人体胸腹区域的准确定位,但对于整个视频的放大处理和光流编码需要庞大的计算开销,最终影响整体算法实时性。

针对以上问题,本文提出了一种基于运动特征预估的呼吸率快速检测方法。通过对视频中呼吸运动区域和运动方向的预估,实现基于相位的呼吸率快速检测。首先通过人脸检测确定人脸区域,并结合人体头身比选出胸部区域,然后建立呼吸信号模型,在人体呼吸率有效范围内利用最大似然法精准定位呼吸区域,并且通过梯度信息确定呼吸区域内主要运动方向。在获取有效的呼吸区域后,结合呼吸运动方向,采用基于相位的复可控金字塔算法提取呼吸信号,最后根据峰值检测获取呼吸率。

呼吸率检测方法由呼吸区域选取和呼吸信号分析两部分组成。精确的呼吸区域可以减少后续呼吸信号处理的数据量。本文初步选取到胸部区域后,进一步通过最大似然估计选取更精确的呼吸区域。人脸检测仅需一帧图像即可实现,因此相较于文献[3]、文献[5]能够更快选取呼吸区域,后续经最大似然估计选取的呼吸区域较文献[4]更为精确,提升了实时性。此外,本文提出使用方向梯度直方图算法[6](histograms of oriented gradients,HOG)获取运动的主要方向,在后续使用复可控金字塔分解获取呼吸相位信号时仅沿主要方向进行分解,相较文献[4]有效减少了信号处理的计算量,提升了实时性。本文通过获取更精确的呼吸区域和呼吸运动主要方向减少了计算量,提升了实时性。

1 本文算法

本文提出的呼吸率检测算法以基于运动特征的预估为基础,通过人脸检测初步选出胸部区域,并建立信号模型,通过最大似然法精确定位胸部呼吸区域。然后通过梯度信息确定主要运动方向,为复可控金字塔的方向选取提供依据。最后采用基于相位的处理方法,提取呼吸信号,并根据峰值检测获得呼吸率。呼吸率检测的流程图如图1所示,图1包括4个主要步骤:① 基于最大似然的呼吸区域预估;② 基于方向梯度直方图的运动方向预估;③ 呼吸相位信号提取;④ 基于峰值检测的呼吸率估计。

1.1 基于最大似然的呼吸区域预估

基于视觉的人体呼吸率检测方法通过检测人体胸部一段时间内起伏的次数来实现。由于胸部特征的提取存在困难,没有直接的胸部识别算法。因此,本文采用人脸检测算法首先定位人脸位置,然后根据人体比例关系,由人脸区域初步确定胸部区域位置。此时的呼吸区域的选取没有针对具体呼吸情况,选取区域过大,会造成检测算法计算量大,实时性低。为了提高算法的实时性,对区域内的像素建立信号模型,通过最大似然估计计算信号基频和运动幅度,根据运动幅度较大的点所在位置选取较为准确的呼吸区域。

1.1.1 呼吸区域初步定位

人脸检测是胸部定位的第1步,本文使用基于肤色模型[7]的人脸检测方法来获取人脸位置。通过人脸检测器进行人脸定位后,使用矩形框框选人脸范围并获取4个顶点的坐标,为胸部区域初定位提供位置信息。

图1 算法流程图

在获取人脸范围后,根据头身比获取胸部区域。成人身高一般是头高的8倍,结合实验结果本文选择第6头高和第7头高之间作为胸部区域。具体区域为:从头顶开始分别向下移动1.5倍头高和2倍头高作为纵坐标,横坐标保持不变,将4点相连初步获得胸部区域,然后进行基于最大似然的呼吸区域定位。呼吸区域选取示意图如图2所示。

图2 呼吸区域选取示意图

图2中：红色区域为人脸区域;蓝色区域为经头身比选取的胸部区域;绿色区域为本文选取的精确呼吸区域。

1.1.2 呼吸区域精确定位

假设所选胸部区域的像素亮度的变化仅由呼吸运动引起,将表示这种变化的周期信号建模为:

X[n]=C+Acos(2πf0TSn+Φ)+W[n]

(1)

其中:矩阵大小为所选胸部区域大小,即M1×M2;C为视频中不变的量;A为幅度;f0为基频;TS为采样周期;n为帧序列;Φ为初始相位;W[n]为高斯噪声。本文使用最大似然对基频f0和幅度A进行估计,首先对基频f0进行估计,然后根据基频f0估计幅度A,最后根据幅度A来选取包含呼吸运动的区域。

为简化符号,引入表示视频帧X[n]的矩阵向量化算子:

XV[n]=vec(X[n])=

[x(1,1,n) …x(1,M2,n)

x(2,1,n) …x(2,M2,n)

⋮

x(M1,1,n) …x(M1,M2,n)]T

(2)

其中:vec为矢量化函数,将矩阵转化为列向量;列向量XV[n]的大小为M1M2×1;x(i,j,n)为视频帧X[n]矩阵中第n帧中坐标(i,j)处的亮度值。相应的逆算子被定义为

X[n]=vec-1(XV[n])

(3)

av[n]cos(2πf0TSn+φv[p])]2

(4)

(5)

由于呼吸运动是周期运动,设呼吸率范围为[fmin,fmax],其中fmin、fmax分别为呼吸率的最低值和最高值。(5)式可以改写为:

(6)

其中:kmin和kmax分别为能够取到的最接近fmin和fmax的频率值;Sp[k]为离散的周期图谱[8]。第p个像素位置[9]的离散值可表示为:

(7)

(8)

经处理后,人体呼吸区域如图2b所示,绿色矩形区域为所选精确呼吸区域。由图2可知,经过最大似然估计,呼吸区域显著减小,算法所需计算量减少,整体算法速度明显提升。以下是基于最大似然的呼吸区域选取算法。

输入:由1.1.1节初步选取的呼吸区域内视频帧X[n];呼吸区域内像素点集合Points;人体呼吸率最高值和最低值fmax、fmin。

1.矢量化处理。简化符号,引入表示视频帧X[n]的矩阵向量化算子XV[n]。

for each row ofX[n]

point ofXV[n]=row ofX[n];

end

for eachfbetween [fmin,fmax]

for each point of Points

for each frame ofX[n]

calculate the discrete periodogramSp[k];

end

end

calculate Maximum likelihood estimate off;

end

for each point of Points

for each frame ofX[n]

end

end

for each point of Points

end

find biggest element ofAasA[p];

select region aroundA[p] asXV[p]

1.2 基于HOG的运动方向预估

通过方向梯度直方图算法HOG[6]获取运动的主要方向,为复可控金字塔中滤波器的方向选择提供依据。

HOG算法通过灰度化将视频帧转化为二维矩阵,并将矩阵分解为一些大小为16×16的块,且相邻块重叠50%区域。由以下公式计算块中每点的梯度信息:

(9)

(10)

(11)

(12)

(13)

图3 方向梯度效果图

图3表示图2b中呼吸区域对应的S矩阵效果图。由图3可知,此时胸部呼吸的主要运动方向为2方向,对应45°。后续使用复可控金字塔分解视频图像时,仅沿45°方向分解,通过减少信号通道数,提升检测算法实时性。以下是基于HOG的运动方向预估算法。

输入:由1.1.2节选取的呼吸区域內视频帧灰度矩阵X[n];呼吸区域內像素点集合Points。

1.建立每帧对应串联矩阵C。将X[n]每帧分解为若干等大小的块Blocks;计算每块中像素点梯度大小和方向。

Decompose each frame into Blocks;

for each point in each block of Blocks

calculate gradient magnitudeG(x,y);

calculate gradient directionΘ(x,y);

end

计算每块对应的向量,进而获取串联矩阵C。

for each block of Blocks

calculate the vector corresponding to block;

end

concatenate all vectors to get the MatrixC

2.计算和矩阵S。将所有C矩阵取平均计算平均矩阵M;使用C矩阵和M矩阵计算和矩阵S。

take average of all MatrixCto get MatrixM。

for each frame

MatrixS=S+|C-M|;

end

for eachθbetween [0,π]

calculate Maximum likelihood estimate ofθ;

end

1.3 呼吸运动相位信号提取

完成对胸部运动特征的预估后,接下来需要提取呼吸信号。如图 1所示,首先在空间上,使用复数可控金字塔对视频序列沿主要运动方向进行不同尺度的子带分解,获取相位信息;然后在时域上,对子带进行带通滤波、将获取到的序列取平均得到表征呼吸运动的子带序列;最后将所有序列求和得到基于相位的呼吸信号。

1.3.1 空间分解

胸部起伏在视频中表现为像素点亮度的周期性变化。为了减少计算量并突出亮度信息,本文方法将色彩空间由RGB域转换到Gray域,将三通道整合为单一通道。RGB与Gray-Scale转换公式为:

Igray=0.298 9R+0.587 0G+0.114 0B

(14)

其中:R、G、B分别为红、绿、蓝3个颜色通道矩阵;Igray为转换后的灰度矩阵。

为了在频率空间中提取相位信息,本文通过傅里叶变换实现图像从时域空间到频域空间的转换,并在频域空间进行相位信息提取、理想带通滤波等操作用于提取呼吸信号。根据文献[10],本文采用复可控金字塔沿主要运动方向将空间域拓展至不同相位域,通过多次迭代将图像分解为不同尺度的子带信息,图像的相位信息蕴含在子带信息中。以下是复可控金字塔的实现过程。

(1) 根据输入图像尺寸,确定金字塔分解层数。计算公式为：

(15)

其中,h、w为图像的高和宽。由(15)式可知,选取更加精确的呼吸区域后,金字塔分解层数变少,算法时间复杂度降低。

(2) 滤波器方向选取。首先设置方向参数N,即金字塔每层共有N个方向(本文N=4),然后根据1.2节结果选取主要运动方向,在后续金字塔分解中仅沿主要运动方向进行不同尺度分解。

(3) 算法首先将输入图像分解为高通子带H0与低通子带L0,其中低通子带包含图像全局信息,高通子带包含图像细节。然后按照参数N与1.2节结果进行主要方向带通滤波,获得主要方向的带通子带序列Bk(k=0,1,…,N)和子带部分L1。

(4) 将子带部分L1进行二抽样后,按步骤(3)中低通子带分解方法继续分解,不断重复这一步骤,直到尺度分解层数达到M。

操作完成后,共获得M+2(H0+M+1)个子带序列,其中M个带通子带序列B中包含呼吸相位信息。至此,方向可控金字塔分解完成。相较于原有N方向分解,仅沿主要方向分解后的通道数显著减少,算法时间复杂度降低。

1.3.2 呼吸信号提取与呼吸率计算

成人呼吸频率一般分为3种情况:呼吸过速(超过24次/min)、呼吸正常(12～24次/min)、呼吸过缓(低于12次/min)。因此本文采用理想带通滤波器对子带序列进行滤波,通频带为0.1～1.0 Hz(对应呼吸率6～60 次/min)。

获取呼吸运动序列后,根据文献[11],将序列转换为平均相位信号来表征呼吸运动,平均相位信号为:

(16)

其中,Wl、Hl分别为各级金字塔分解后的图像尺寸。

将多通道信号按帧相加,计算呼吸区域内平均相位随时间的变化序列得到呼吸信号，视频提取的呼吸波形如图4所示。

图4 视频提取的呼吸波形

呼吸频率的计算公式为:

(17)

其中:N为视频总帧数;FS为帧率;P为波峰检测算法检测到的呼吸波形中波谷数目。

2 实验结果与分析

2.1 实验设备及数据

实验通过MATLAB软件调用摄像机采集完成的视频,并对采集视频进行呼吸率估计。实验中被测者与摄像机水平相距0.5 m,并拍摄人体上半身图像,如图5所示。实验在室内稳定白炽灯环境下使用普通光学摄像机采集视频,在分辨率1 280×720像素,帧率30帧/s条件下,在RGB空间下保存成AVI格式。

图5 呼吸率检测设备场景

由于成人呼吸频率可分为呼吸过速、正常与呼吸过缓3种状态,考虑呼吸状态与测量对象的差异性[12],实验要求测试者在呼吸频率6～36 次/min范围内进行周期呼吸,对测试者进行多次视频采集和呼吸率检测,验证算法的准确性。测试者按照显示屏提示进行特定频次的呼吸作为真实呼吸率进行参考。

2.2 实验结果

实验复现了文献[3-5]的视频呼吸率检测方法进行本文方法的性能验证。本文采用3种指标进行性能评估。第1个指标为平均误差Me,即真实值与检测值的偏差,即

(18)

其中:N为视频的数目;Rev为呼吸率估计值;Rtrue为真实呼吸率值。

第2个指标是平均准确率,记做Rac,根据相关医学指导建议[13],以误差±3 次/min作为正确组,其余为错误组。

(19)

第3个指标为视频时长与算法用时的比例,记做Rt,用来反映实时性能。

Rt=tv/td

(20)

其中：tv为视频时长;td为获取呼吸率所需时间。Rt越大表示实时性越好。

实验在室内稳定白炽灯条件下对测试者进行呼吸率测量,并与其他算法进行实时性和稳定性的比较,最后进行各步骤贡献度验证实验。

2.2.1 检测性能比较

本次实验在视频时长60 s,采样率30帧/s,室内稳定白炽灯条件下进行。将采集的视频采用文献[3-5]方法和本文方法检测呼吸率,并根据检测结果进行算法性能比较。在上述实验条件下,呼吸率检测性能的记录见表1所列。为直观解释表格数据,本文选取30组实验数据制成散点图,如图6所示。

表1 理想环境下呼吸率检测性能对比

由表1和图6可知,本文方法与文献[3-5]方法在呼吸率检测准确率上均保持了很好的性能。文献[3]方法直接对全图建立呼吸模型并进行最大似然估计,能够准确地定位到呼吸区域,从而保证了呼吸率检测的准确性,然而算法对全图进行的呼吸信号建模和最大似然估计都需要计算开销,影响了算法的运行速度;文献[4]的方法虽然进行了胸口定位和基于相位的图像处理,但选取的呼吸区域过大,并且由于没有选取主要方向,增加了信号通道数,导致计算开销增大,降低了算法实时性;文献[5]的方法通过人脸检测加光流编码处理,精确定位了运动区域,且使用欧拉视频放大[14](Eulerian video magnification,EVM)算法进行预处理,但对于全图的EVM处理和光流编码在视频质量较高、分辨率较大时,影响算法整体的实时性。

图6 30组呼吸率散点图

本文方法在提取呼吸信号前用人脸检测加最大似然估计法精确定位呼吸区域,减少了需要处理的像素数。然后通过确定主要运动方向,在呼吸信号提取阶段时,仅沿主要运动方向进行时空处理,有效地减少了金字塔层数与通道数,降低了计算复杂度,在保持检测精度的同时提高了算法实时性。

2.2.2 算法稳定性验证

稳定性验证可以用于评估算法是否适用于真实场景。本节通过不同时长、不同衣服特征量、不同呼吸区域大小、不同帧率4种情况下分别进行呼吸率检测,并与文献[3-5]方法进行对比,验证方法的稳定性。

对不同视频时长条件下的呼吸率检测进行稳定性验证。在光照条件为普通白炽灯稳定光源环境,呼吸区域大小64×64,普通带花纹T恤条件下,分别采集10、20、40、60 s时长的视频进行验证。计算4种情况下的平均准确率Rac和Rt,见表2、表3所列。由于文献[3]和文献[5]的方法在短时长下失去稳定性,仅对文献[4]方法和本文方法进行实时性的比较。

由表2可知,随着视频时长的减少,4种方法呼吸率检测的平均准确率均逐渐下降,其中文献[3]方法、文献[5]方法在10 s视频时长条件下已经失去稳定性。文献[4]方法和本文方法由于通过金字塔分解获得了呼吸运动多通道的子带信息,在视频时长较短的情况下也能弥补数据量较少产生的估计偏差问题。由表3可知,随着视频时长的减少,本文方法虽然实时性一直优于文献[4]方法,但也在下降，这是因为对于任意长度视频,本文算法均花费一定时间选取呼吸区域与运动方向,当时长变短,后续呼吸信号处理减少的时间不足以弥补预处理的时间花费，从而导致对短视频进行呼吸率检测时,本文方法实时性有所下降。

表2 不同时长下的呼吸率检测平均准确率 单位：%

表3 不同时长下的呼吸率检测用时比

对不同衣服特征量条件下的呼吸率检测进行稳定性验证。在光照条件为普通白炽灯稳定光源环境,呼吸区域大小为64×64,视频采集时间为60 s条件下,分别对衣服类型为纯色T恤、花纹T恤2种特征量不同的衣服进行测试,具体衣服类型如图7所示,并以平均准确率Rac作为衡量标准进行稳定性验证,结果见表4所列。

图7 不同衣服特征量效果图

表4 不同衣服特征量下的呼吸率检测平均准确率 单位：%

从表4数据可以看出,当衣服特征量较少时,文献[4]方法、文献[5]方法已经失去了稳定性,而文献[3]方法、本文方法仍保持良好的性能。主要因为文献[4]和文献[5]的方法均使用了阈值进行二值化处理,当所选呼吸区域都缺少特征时,无论是基于亮度的差值[15]还是基于相位的的差值均可能会因为小于阈值而被忽略,导致提取不到运动。在相同视频质量的条件下,文献[3]方法和本文方法没有使用二值化,通过最大似然估计法仍能精确获取呼吸运动信号,因此较其他方法准确率会有一定程度提升。

在不同呼吸区域尺寸条件下进行呼吸率检测稳定性验证。经最大似然估计后选取到胸口区域振幅最大的点的位置,以此位置为中心可选取不同尺寸的呼吸区域。在光照条件为普通白炽灯稳定光源环境,普通带花纹T恤,视频采集时间为60 s条件下进行测试。分别采用128×128、64×64、32×32、16×16进行测试,具体区域大小如图8所示。在2.2.1节性能实验的条件下,分别以50、30、15帧率进行测试。最终以平均准确率Rac、实时性Rt作为衡量标准进行稳定性验证,结果见表5所列。

从表5可以看出,随着呼吸区域逐渐减小,Rt逐渐增大,实时性增强。当呼吸区域过小时,平均准确率Rac显著下降,本文算法失去稳定性。这是因为呼吸区域内包含的数据量太少,影响后续呼吸率提取。视频帧率方面,当帧率为15帧时,所要处理数据量变小,实时性提高,但平均误差和平均准确率均会下降。当帧率为50帧时,数据量增大,但相较30帧率准确率没有显著提高且实时性显著下降。

图8 不同呼吸区域尺寸下的显示效果图

表5 不同呼吸区域尺寸下的呼吸率检测准确性和实时性

2.2.3 逐步测试

为了验证呼吸区域选取和主要运动方向选取对本文算法的贡献度,本文进行逐步测试。在2.2.1节性能实验的条件下,分别在仅选取呼吸区域、仅选取主要运动方向、本文完整算法条件下进行测试。并以平均准确率Rac、实时性Rt作为衡量标准进行验证,结果见表6所列。

表6 不同条件下的呼吸率检测准确性和实时性

从表6中可以看出,在实时性方面呼吸区域选取贡献了34.6%,主要运动方向选取贡献了65.4%,通过精确选取呼吸区域,可有效减少所需计算的像素点数量,进而降低算法时间复杂度;通过确定主要运动方向,在后续金字塔分解中可减少3/4的通道数,有效降低了时间复杂度。在准确性方面,3种方法都保持了良好的性能,说明基于最大似然的呼吸区域和主要运动方向的选取能够有效反映人体呼吸时胸腔的主要运动情况。

3 结 论

本文提出了一种基于运动特征预估的呼吸率检测快速方法。通过基于最大似然的呼吸区域预估和基于HOG的呼吸方向预估为后续视频信号处理降低时间复杂度;最后采用基于相位的欧拉视频处理方法对呼吸区域进行快速时空处理,获取呼吸信号,提取呼吸频率。

在真实场景中,当被测者出现大的运动或定位不到呼吸区域时,本文方法无法准确地检测呼吸率。另外,由于呼吸率的精度要求,当视频时长小于10 s时,会因为频率分辨率太低影响检测结果。在未来的工作中,可针对呼吸运动中存在的大运动干扰、各种情况下的呼吸区域定位和减少检测所需视频时长等方面进行改进,进一步扩大算法的应用范围。