李晓媛,武 鹏,刘 允,司红玉,王振龙
(1.郑州大学 电气工程学院,河南 郑州 450001;2.郑州大学 体育学院,河南 郑州 450001;3.郑州大学 生命科学学院,河南 郑州 450001)
心率是指每分钟心脏跳动的次数,它作为人体最重要的生理参数之一,已经被广泛应用于心律失常、心肌缺血、高血压病、慢性心率衰竭等心血管疾病的诊断和情绪检测当中[1-5]。目前,心率测量的金标准是心电图法(Electro Cardio Gram,ECG),这种方法测量的准确度虽然高,但是需要专业医护人员在测试者皮肤表面粘贴电极,给测试者造成一定的不适,也会增加皮肤感染的风险,不适合长时间测量。此外,接触式心率测量也不适用皮肤受损患者,新生婴儿以及隐蔽式刑侦测谎等情形。
通过非接触的方式来测量心率参数是近年来提出的一项新技术。2007年,Pavlidis等[6]通过分析面部热红外图像成功提取到被试者的心率参数。2010年,路国华等[7]使用多普勒雷达测得了被试者的心率参数。但这些设备价格昂贵,而且需要复杂的硬件支持,很难推广到实际应用。Poh等[8]通过低成本的普通摄像头配合独立成分分析法(Independent Component Analysis,ICA)成功提取到了被试者的心率信号,证明了用低成本的摄像头来提取心率信号的可行性。其原理是脸部皮肤表面的血液容积在心脏的作用下呈现周期性变化,因此可以通过提取脸部的光电容积脉搏波(Photoplethysmography,PPG)信号来获得心率参数。2011年,Poh等[9]在原来研究的基础上,利用计算波峰间隔(Interbeat Intervals,IBIs)算法来获取心率参数。2013年,Bousefsaf等[10]利用Morlet小波变换算法重构原始心率波形,并通过IBIs算法来获得心率参数。2014年,Monkaresi等[11]将机器学习算法应用于心率信号的提取,提出了ICA+KNN(K-Nearest Neighbor)的心率波形重构算法,并使用快速傅里叶变换(Fast Fourier Transform,FFT)算法从重构的波形中提取心率参数。2017年,吴庆甜等[12]用联合近似特征对角化(Joint Approximate Diagnoalization of Eigenmatrices,JADE)算法和FFT算法提取驾驶员的心率信息,用于评估驾驶员的疲劳状态。2019年,Martinez等[13]利用红外摄像头提取脸部PPG信号,并用短时傅里叶变换算法提取心率参数。
目前,研究人员已经尝试用多种方法从摄像头录制的视频中提取高质量的PPG信号;但是通过PPG信号计算心率参数的方法却很少,大体可以分为两类,一类是IBIs算法,一类是傅里叶变换算法。IBIs算法通过计算PPG信号波峰的时间间隔来获取心率参数,但由于噪声干扰的存在,会出现波峰偏移、漏检、多检等情形,这会严重影响心率参数测量的准确性。而傅里叶变换只适合计算一段时间内占主要成分的心率信息,不能够得到心率变化的信息,即使利用短时傅里叶变换,在时域特性与频域特性之间权衡时间窗的大小也是一项难题。针对上述问题,本文从YCbCrCg颜色空间Cg通路中提取高质量的PPG信号,用Morlet复小波(Complex Morlet,CMOR)作为母波绘制小波能量谱图,并根据心率参数的生理特性去除伪点噪声,提取随时间变化的心率参数。最后,通过计算静息状态和头部运动状态下非接触式方法测量结果同标准仪器测量结果的平均绝对值误差|Me|,误差的标准差SDe和均方根误差(Root Mean Square Error,RMSE),并绘制两种测量方法的Bland-Altman散点图,验证该方法的准确性。
非接触式心率测量原理在于当心脏跳动时,心室的收缩和舒张会引起血管内血液容积的周期性变化,其反射光的强度会呈现周期性变化,人体脸部皮肤表面含有丰富的毛细血管,通过摄像头记录和特定算法提取分析脸部PPG信号的周期性变化,可以测得相应的心率参数。
本实验招募8名测试者对设计的非接触式心率测量系统进行验证,测试者的性别、年龄、肤色信息如表1所示,其中肤色根据Fitzpatrick肤色分型法[14]分为Ⅰ~Ⅵ型,Ⅰ型对应于白色皮肤,Ⅵ型对应于黑色皮肤。本实验分别在静息状态和头部运动状态下使用普通的网络摄像头(Logitech C920)对每位测试者录制60 s视频。静息状态时,测试者头部面朝摄像头尽量保持静止。头部运动状态时,测试者头部面朝摄像头左右倾斜,倾斜角度约为45°。一个完整的运动周期为:头部居中-倾斜至一侧约45°-倾斜至另一侧约45°-头部居中,运动周期约为6 s。为保证摄像头读入帧率维持在20 frame/s,摄像头的分辨率设为800×600 pixel。实验需关闭摄像头的自动白平衡功能来避免其自动调节局部颜色从而引入噪声信号。数据分析与处理在一款中档台式电脑上进行(Intel core i5 处理器,运行内存为 8 G,系统为Windows 7),所有实验均使用指夹式脉搏血氧仪同步记录心率参数。
表1 测试者基本信息
注:M代表男性,F代表女性。
非接触式心率测量系统分为视频图像处理和信号分析两部分。
2.3.1 视频图像处理
视频图像处理中,先通过摄像头录制一段包含人脸的视频图像(图1(a)),并使用Kanade Lucas Tomasi(KLT)算法[15]对人脸进行识别与跟踪(图1(b)),KLT算法返回矩形人脸框4个顶点的坐标。如果头部倾斜,则检测到的矩形人脸框也是倾斜的,本文设计了相应的角度转换算法来提取倾斜的矩形人脸数据并对它进行重构。人体脸部皮肤表面含有丰富的毛细血管,其反射光构成原始的PPG信号,因此将重构的人脸图像转换到YCbCr颜色空间来进行皮肤检测(图1(c)),转换公式如式(1)所示:
(1)
其中:Y为像素的亮度分量,Cb和Cr分别为蓝色和红色的浓度偏移量成分。为提取信噪比较高的PPG信号,进一步计算了Cg通道即绿色的浓度偏移量成分,如式(2)所示:
(2)
其中:Y′=Kr×R′+Kg×G′+Kb×B′,R′,G′,B′表示摄像头录制视频的三个原始通道,Kb,Kr,Kg为权重因子,且Kb+Kr+Kg=1。参考ITU-R的BT.601协议,取Kb=0.114,Kr=0.299,则Kg=0.587,Y′=0.299×R′+0.587×G′+0.144×B′,代入式(1)和式(2)计算可得:
(3)
(4)
根据亚洲人的皮肤颜色特点,Y,Cb,Cr3个通道的皮肤颜色检测阈值设置如下:
(5)
脸部皮肤表面的血液容积在心跳的作用下呈现周期性变化,这种变化是PPG信号交流分量的重要组成部分。研究表明,血液中红细胞内氧和脱氧血红蛋白对510~590 nm光谱段最为敏感[16],这段光谱对应于可见光的绿/黄光部分,因此本文选择光谱范围较为接近的Cg颜色通道来提取心率信号。通过式(5)确定脸部皮肤位置(图1(d)),提取Cg颜色通道中皮肤位置对应点的像素(图1(e)),对像素强度进行平均得到原始PPG信号x(t)(图1(f))。
图1 通过图像处理算法提取原始PPG信号.(a)摄像头采集一段包含人脸的视频;(b)从视频中提取人脸部分图像;(c)将提取到的图像转换到YCbCr颜色空间;(d)在YCbCr颜色空间进行人脸皮肤检测;(e)将Cg颜色通道的人脸图像和检测到的人脸皮肤位置进行逻辑‘与’操作,提取Cg颜色通道中人脸皮肤位置对应点的像素;(f)对像素强度进行平均得到原始PPG信号
Fig.1 Original PPG signal extracted by image processing algorithm. (a)Camera captures a video containing the face;(b)Face images extracted from video;(c)Extracted images converted to YCbCr color space;(d)Face skin detection performed in YCbCr color space;(e)Logical ‘AND’ operation performed on the face image of theCgcolor channel and the detected face skin position, extracting the pixels corresponding to the position of the face skin in theCgcolor channel;(f)Original PPG signal obtained by averaging the pixel intensity
2.3.2 基于CMOR的心率信号分析
通过原始PPG信号x(t)提取心率数据时,首先对原始信号进行带通滤波,通带频率设置为[40/60,200/60] Hz,对应于心率40~200 bpm(beats per minute),滤除通带外的噪声信号。
小波变换能够同时分析信号的时域分量和频域分量,而傅里叶变换只能分析信号的频域分量。短时傅里叶变换虽然能够通过加窗操作同时观察到信号的时域分量和频域分量;但是短时傅里叶变换窗函数宽度一旦确定就不能改变,窗函数宽度选择较窄,则频域分辨率较差,窗函数宽度选择较宽,则时域分辨率较差,因此窗函数宽度的选择也是一项难题。小波变换由于其窗口形状可变,可以同时在时域和频域分析信号的细节变化,因而广泛应用于生物信号分析中[17-21]。小波变换中,子小波ψτ,s共轭后对PPG信号x(t)进行卷积得到小波变换系数:
(6)
其中:*表示共轭,子小波ψτ,s(t)是通过母小波或基小波ψ(t)伸缩和平移得到的,如式(7)所示:
(7)
其中:s表示伸缩因子或尺度因子,对应分析频率,增大s时,分析频率降低,减少s时,分析频率升高;τ表示平移因子,对应位移信息。小波分析是一种强有力的信号分析工具,想要得到理想的分析效果,根据输入信号选取合适的母小波是一项至关重要的操作。本研究通过大量的实验对比,最终选定CMOR作为母小波来分析PPG信号x(t),CMOR母小波的表达式为:
(8)
其中:fd为带宽参数,fc为小波函数的中心频率。根据PPG信号x(t)的生物特性,本研究中选择fd=5,fc=3,代入式(8)可得:
(9)
它对应的子小波函数为:
(10)
图2(a)是提取的一段原始PPG信号x(t),横轴表示测试时间,纵轴表示像素强度。图2(b)是对它进行CMOR小波变换后生成的能量谱图。本研究开发了相应的算法来提取随时间变化的心率参数,第一步提取能量谱矩阵中每列能量值最大的位置,第二步根据心率参数变化连续的生理特征去除伪点,并绘制随时间变化的心率参数曲线,如图2(c)所示。
(a)原始PPG信号x(t) (a)Original PPG signal x(t)
(b)小波变换能量谱 (b)Energy spectrum of wavelet transform
(c)心率变化曲线 (c)Heart rate curve图2 由PPG信号x(t)中提取心率信号Fig.2 Heart rate signal extracted from PPG signal x(t)
本文利用提出的方法对所有的测试者进行非接触式心率信号检测,并用标准指夹式脉搏血氧仪同步记录心率参数进行对比。其中,第三位测试者在静息状态下的测量结果如图3所示(彩图见期刊电子版),红色表示非接触式算法的测量结果,蓝色表示指夹式脉搏血氧仪的测量结果。
图3 静息状态下两种方法测量结果对比Fig.3 Comparison of results of two measurement methods in resting state
对于每一位测试者,每隔0.5 s取一次两种方法的测量结果,计算静息状态下非接触式算法测量结果同标准接触式仪器测量结果的平均绝对值误差|Me|,误差的标准差SDe和RMSE,结果如表2所示。结果表明,对于所有测试者,两种测量方法的|Me|均小于2 bpm,远低于中华人民共和国医药行业规定的误差标准(误差≤5 bpm),SDe均小于2.5 bpm,RMSE均小于2.6 bpm,表明静息状态下本文提出的非接触式算法同标准仪器的测量结果高度一致。
表2 静息状态下两种测量方法的|Me|,SDe和RMSE 统计结果
图4 静息状态下Bland-Altman一致性分析Fig.4 Bland-Altman consistency analysis in resting state
对于每一位测试者,分别计算头部运动状态下非接触式算法测量结果同标准接触式仪器测量结果的|Me|,SDe和RMSE,结果如表3所示。结果表明,对于所有测试者,头部运动状态下两种测量方法的|Me|均小于2.3 bpm,SDe均小于2.9 bpm,RMSE均小于2.9 bpm。与静息状态相比,头部运动状态下由于运动伪影等噪声干扰的存在,非接触式算法的测量误差有所升高,但依然保持较高的测量精度。
表3 头部运动状态下两种测量方法的|Me|,SDe和RMSE统计结果
Tab.3 Statistical results of |Me|, SDe and RMSE of two measurement methods in head moving state (bpm)
图5 头部运动状态下Bland-Altman一致性分析Fig.5 Bland-Altman consistency analysis in head moving state
心率是一项极为重要的生理参数,被广泛应用于心血管疾病的诊断和情绪检测中,本文提出了一种基于CMOR小波变换的非接触式心率测量方法,用低成本的电脑摄像头拍摄人脸来准确检测被试者的心率参数。并通过计算静息状态和头部运动状态下非接触式算法测量结果同标准接触式仪器测量结果的平均绝对值误差|Me|,误差的标准差SDe和RMSE,绘制两种测量方法的Bland-Altman散点图,证明本文提出的非接触式方法测量结果同标准仪器测量结果具有较强的一致性。
在视频图像处理中,人脸皮肤检测是一项重要的操作,KLT算法跟踪提取到的矩形人脸像素中,存在头发、眼睛、背景等像素引起的噪声信号,直接提取到的PPG信号信噪比较低,如图6(a)所示,经过皮肤检测后提取到PPG信号的信噪比较高,如图6(b)所示。而目前的皮肤检测算法基本都是通过在YCbCr颜色空间中设定Y,Cb,Cr 3个通道的阈值来确定的,但对于不同的肤色类型,如I型白色皮肤,Ⅵ型黑色皮肤,Y,Cb,Cr 3个通道的阈值需要重新设定才能够检测到皮肤像素,未来准备开发针对脸部的自适应肤色检测算法,根据测试者脸部皮肤类型,算法自动调整Y,Cb,Cr 3个通道的阈值,从而进一步提高算法的准确性和鲁棒性。
(a)未经过皮肤检测提取到的PPG信号 (a)PPG signal extracted without skin detection
(b)经过皮肤检测提取到的PPG信号 (b)PPG signal extracted through skin detection
图6 未经皮肤检测与经过皮肤检测提取到的PPG信号对比
Fig.6 Comparison of PPG signals extracted without skin detection with those extracted through skin detection
本文使用Cg颜色通路作为信号源来减少运动伪影等噪声干扰信号的影响。例如,当测试者头部出现自然微小运动时(5 s和20 s),原始RGB视频中G颜色通路的PPG信号会有较大波动,如图7(a)所示;而Cg颜色通路的PPG信号基本不受影响,如图7(b)所示,证明Cg颜色通路作为信号源能够有效减少运动伪影等噪声干扰的影响,这是本文提出的方法能够在头部运动状态下保持较高精度的一个重要原因。其中,G颜色通路已经被验证是R,G,B 3个通路中PPG信号信噪比最高的一个通路[23]。
(a)G颜色通路提取到的PPG信号 (a)PPG signal extracted from G channel
(b)Cg颜色通路提取到的PPG信号 (b)PPG signal extracted from Cg channel
图7 G颜色通路和Cg颜色通路提取到的PPG信号对比
Fig.7 Comparison of PPG signals extracted from G channel with those extracted from Cg channel
在之前的研究中,大多数学者将滤波得到的PPG信号进行傅里叶变换[24-26],找到频率域中幅值最大点,进而确定该时间段的心率参数,如图8(a)所示。本文使用CMOR小波算法来生成PPG信号的能量谱图,同时分析信号的时域分量和频域分量,并根据心率参数变化连续的生理特征去除伪点噪声,提取随时间变化的心率参数,测量精度高,克服了以往傅里叶变换只能提取时间段内单一主心率参数的缺陷,如图8(b)所示。前一种方法只适用于在心率平稳时间段内提取主心率参数,如果该时间段心率有一定波动,则该方法也只能检测到频率域中幅值最大点对应的单一心率参数;而本文提出的方法能够得到随时间变化的心率数据,且实验结果表明,本文方法的测量结果和标准仪器的测量结果具有较强的一致性。
(a)傅里叶变换提取到的心率参数 (a)Heart rate parameters extracted by Fourier transform
(b)CMOR小波能量谱图提取到的心率参数 (b)Heart rate parameters extracted by energy spectrum of CMOR wavelet
图8 傅里叶变换和CMOR小波能量谱图提取到的心率参数对比
Fig.8 Comparison of heart rate parameters extracted by Fourier transform with those extracted by energy spectrum of CMOR wavelet