1 引言

心脏病是人类死亡率最高的疾病, 提高早期诊断的效率、有效地进行防范和治疗, 能够提高人们的身体健康水平. 心率是临床心脏疾病检测的重要指标之一, 保持心率的日常稳定对于心脏疾病的防范治疗有很重要的意义^[1]. 测量心率的方法有很多种, 传统中医通过摸桡动脉的搏动, 用来间接测定心率情况. 医生通过听诊器测量, 在心前区胸骨左缘第4肋间左右每分钟听到的心音数, 就是心率^[2]. 更为精准的测量方法有, 动态心电图和心脏彩超等检查, 需要利用电极和辅助带固定, 常用于医院重大疾病的检测. 普通情况下人们采用手环, 指夹式脉搏血氧仪等设备检测心率. 上述方法都是接触式心率检测, 在舒适性和便捷性方面存在缺点. 近年来针对心率的非接触测量技术被深入研究.

目前主流的非接触式心率测量方法是成像式光电容积脉搏波描记法(Imaging PhotoPlethysmoGraphy, IPPG)^[3]. 该方法认为心脏跳动带来的血液循环导致人体血管中血液容积的周期性变化, 在光源照射皮肤时, 血管中不同容积的血液对于光的吸收也随之发生周期性变化, 相机成像接收到的反射光强度保留了血容量变化脉冲信号(Blood Volume Pulse, BVP). 对摄像头采集的人脸视频进行处理, 可以提取出血液容积的变化信息, 进而通过信号处理算法得到心率等人体生理体征信号. 2008年, Verkruysse等人^[4]首次指出使用环境光作为视频采集环境提取脉搏波信号, 并认为图像绿色通道相对于其他通道对于BVP信号有更好地反映. 2010年, Poh等人^[5]提出了基于盲源分离的心率检测方法, 利用独立成分分析(Independent Component Analysis, ICA)将感兴趣区域(Region Of Interest, ROI)的图像红绿蓝3个通道亮度变化解混分离成3个独立的信号, 选择第2个源信号作为脉搏波信号来估计心率. 2014年, Guo等人^[6]应用独立矢量分析(Independent Vector Analysis, IVA)算法处理面部区域的颜色信号, 实验结果显示, 与基于ICA的方法相比, 测量得到的心率更加准确. 2016年, Chwyl等人^[7]提出了一种基于贝叶斯估计的心率检测方法, 通过贝叶斯最小化框架估计脉搏波信号波形. 上述方法仍存在不足之处, 对于采集过程中的运动伪差和光照不均等干扰问题, 仍需要研究者们进一步深入研究.

本文基于人脸视频图像测量人体心率, 主要工作包括: (1)基于IPPG提出了一种多特征区域结合盲源分离算法的非接触式心率测量方法; (2)通过人脸特征点检测算法确定多个特征区域进行信号提取. 降低了视频拍摄过程中运动伪差对心率测量结果的干扰. (3)通过多特征区域结合快速独立成分分析(FastICA)算法提取脉搏波信号, 降低了视频中存在的人脸面部光照不均问题对于心率测量结果的影响; (4)在公开的大型数据集DEAP^[8]进行实验, 对提出的非接触式测量方法进行了分析和验证.

2 基于人脸视频的心率测量方法

本文提出了一种新的基于IPPG的非接触式心率测量方法. 该方法通过多特征区域结合FastICA算法提取脉搏波信号, 检测心率, 结合了图像绿色通道保留更多BVP信息和盲源分离算法的双重优势. 方法流程图如图1所示, 包括视频帧处理, 特征区域处理, 信号处理, 以及心率计算4个部分, 下文具体对相应的流程进行详细说明.

2.1 视频帧及特征区域处理

基于人脸视频图像的非接触式心率信号检测, 为了提高BVP信号的信噪比和减小系统误差, 需要进行特征区域的选取. 传统的选取方法采用额头^[9,10], 脸颊部位^[11], 整个人脸^[12]作为信号提取的感兴趣区域. 相比之下, 额头部位很容易有遮挡情况出现, 而选择整个人脸区域误差较大. 目前基于IPPG原理的非接触心率测量研究中多为单特征区域, 在视频采集过程中很容易受到光照和运动伪差的影响, 单特征区域提取的BVP信号在后续处理中难以滤除光照和运动带来的干扰. 而选择多个特征区域. 提取的不同信号之间就可以相互补偿, 有效地避免运动伪差对提取的脉搏波信号质量的影响, 基于上述考虑, 本文利用人脸特征点检测结合区域中心定位的方法在两侧脸颊各选取3块矩形部分重叠的特征区域作为最终的信号提取区域.

2.1.1 区域中心定位

本文提出的非接触式实时心率检测方法, 人脸检测及关键点定位算法采用landmark算法, 该算法可以提取到面部68个特征点, 特征点示意如图2所示.

至少选择3个以上特征点作为一个特征点组, 计算特征点组构成的区域质心作为特征区域中心点. 本文方法中每个特征点组选择3个特征点, 三角形的质心计算方法如式(1)所示, 计算三角形3个顶点横纵坐标的和除以3即为三角形质心点坐标.

$x = \frac{{\left( {{x_1} + {x_2} + {x_3}} \right)}}{3},y = \frac{{\left( {{y_1} + {y_2} + {y_3}} \right)}}{3}$

(1)

对3个特征点计算得到的区域中心稳定性进行分析, 如图3, 以鼻尖特征点 $O$ 为参考点, 计算区域中心点A, 特征点组的两个特征点B, C, 相对于鼻尖特征点 $O$ 的向量 $\overrightarrow {O{{A}}} $ , $\overrightarrow {O{{B}}} $ , $\overrightarrow {O{{C}}} $ .

选择一个视频进行测试, 该视频中存在较大姿态变动, 分别计算3个向量在视频帧间X坐标和Y坐标变化幅度的方差, 如表1所示: 可以看到区域中心对应的方差相对于特征点都要更小.

因此根据特征点确定的区域中心点相对于人脸检测的特征点更加稳定, 根据区域中心坐标截取矩形区域作为BVP信号提取的特征区域, 提高了视频帧间的定位和跟踪准确度, 降低了视频采集过程中运动伪差对信号的干扰.

2.1.2 多特征区域选择

根据landmark算法对包含人脸的视频图像处理得到人脸检测框和68个人脸特征点的位置坐标, 为了保证特征区域包含尽可能多的人脸区域, 通过上文提到的区域中心定位对20个特征点组构成的特征区域的中心点稳定性进行分析, 最终确定了6个三点的特征点组如图4所示, R表示右侧脸颊, L表示左侧脸颊. 每一个特征点组构成一个三角区域, 以三角区域的质心作为特征区域的中心, 截取矩形特征区域. 每个特征区域的长宽为人脸检测框的10%, 单侧脸颊的3个特征区域重合度为40%到60%.

对视频中每一帧图像都进行上述的特征区域选择处理, 每一帧图像得到6个特征区域, 如图5所示, 单侧脸颊3个特征区域存在部分重叠, 两侧脸颊特征区域相互对称. 一方面保证了特征区域覆盖更多的脸颊区域, 在处理过程中更容易去除噪声, 保留主要信号. 另一方面, 提取到的多路信号存在更多的相似性, 信号在后续处理过程中相互补偿,可以得到信噪比更好的脉搏波信号.

2.1.3 G-BVP信号提取

通过对视频帧图像的多特征区域选取处理, 得到了6路特征区域的图像序列, 需要对特征区域进行处理提取BVP信号. BVP信号的提取处理包括视频帧图像特征区域的平滑处理和通道均值化. 本文在信号提取部分对6路图像序列进行处理提取到6路G-BVP信号.

平滑处理包括卷积核大小为5的均值滤波和两层金字塔降采样. 通过平滑处理, 去除视频帧图像像素点异常导致的噪声干扰. 平滑之后对特征区域进行通道分离. 由于在图像灰度处理过程中, 图像绿色通道的分量占比权重更大, 反射光强度信息更多. 通过图像绿色通道提取的BVP信号相对于图像红色, 蓝色通道提取的BVP信号更加能够反映人体血液容积的变化信息^[13]. 所以通道分离后选取图像绿色通道进行均值化, 对6个特征区域的图像序列逐帧进行上述处理得到最终的G-BVP信号. 本文对一个测试视频进行上述处理得到的6路G-BVP信号如图6所示, 可以看出图中信号干扰比较严重.

2.2 信号处理

基于人脸视频图像的非接触式心率信号提取, 脉搏波信号质量在视频采集过程中很容易受到光照的影响. 现有研究都是寻找更好的信号分离算法. 在非接触式生理信号检测的研究中, ICA被用于信号的分离以去除噪声干扰, 但ICA通常被用来处理同一个特征区域得到的图像RGB三通道的BVP信号. 本文借鉴该方法, 从视频序列中提取6路G-BVP信号, 再对6路G-BVP信号选用基于负熵最大的FastICA算法进行盲信号分离, 降低光照对脉搏波信号的干扰, 提高信号的信噪比.

2.2.1 脉搏波信号的盲源分离模型

盲源分离是指在信号的理论模型和源信号无法精确获知的情况下, 如何从混迭信号(观测信号)中分离出各源信号的过程^[14]. 本文针对脉搏波信号的盲源分离原理图如图7所示.

图7中 $S(t)$ 为 $n$ 维源信号矢量, 针对本文问题, 假定源信号矢量中一个为脉搏波信号, 其他源信号矢量为运动, 光照等其他干扰信号. 脉搏波信号和干扰信号经过未知的混合矩阵A混合为多维包含干扰的观测信号矩阵X, 每个 $x(t)$ 都是脉搏波信号和干扰信号的组合. X即本文通过多特征区域的选择和处理得到的G-BVP信号. 我们建立数学模型如式(2)所示:

$X(t) = AS(t)$

(2)

由于混合矩阵A和源信号矢量未知, 所以需要通过算法求解一个分离矩阵 $W$ , 得到估计矢量信号 $Y(t)$ , 使 $Y(t)$ 近似于 $S(t)$ , 估计矢量信号中一个为脉搏波信号, 其余信号为干扰信号. 估计矢量用源信号矢量表示, 如式(3)所示.

$Y(t) = WX(t) = WAS(t) = PS(t)$

(3)

式中, 估计信号矢量 $Y(t)$ 和源信号矢量 $S(t)$ 的近似程度用矩阵P来表示. 明显P是一个单位矩阵. 分离矩阵W与混合矩阵A互为逆矩阵.

2.2.2 基于负熵最大的FastICA算法

FastICA是由芬兰赫尔辛基大学Hyvarien等人提出来的盲源分离算法. FastICA算法主要有基于负熵和基于峭度两大类^[15], 由于高斯分布的随机变量具有最大的微分熵, 因此可以利用微分熵最大来判断变量之间的相互独立性^[16]. 本文通过基于负熵的FastICA算法, 从G-BVP信号中分离脉搏波信号计算心率.

FastICA算法首先对观测信号进行预处理以降低输入信号的相关性, 预处理包括去均值和白化处理. 去均值和白化处理可通过式(4)和式(5)实现, 白化矩阵 $V$ 的求解通过输入信号矩阵的协方差矩阵的特征值分解来实现.

$\tilde x(i){\rm{ = }}\tilde x(i) - \frac{1}{n}\sum\nolimits_{i = 1}^n {x(i)} $

(4)

${x_1}(i) = V\tilde x(i)$

(5)

白化处理之后, 以负熵作为判断条件, 按照数学里最优化和矩阵分析的相关知识迭代求解分离矩阵 $W$ . 假设输入 $X(i)$ 的概率密度函数为 ${p_1}$ , 高斯概率密度函数为 ${p_g}$ . 则负熵的定义如式(6)所示:

$J(x) = \int {{p_1}(x)\lg \left[ {\frac{{{p_1}(x)}}{{{p_g}(x)}}} \right]} {\rm{d}}x$

(6)

通过式(6)可以得到 ${p_1}(x)$ 与高斯分布的相似程度. 由于其计算困难, 用一个非二次函数对一个标准化的随机向量的负熵简化近似为式(7):

$J(y) \approx \sum {{k_i}} {\{ E[{G_i}(y)] - E[{G_i}(v)]\} ^2}$

(7)

其中, $v$ 是标准化的高斯变量, 假设 ${{y = }}{{{w}}^{\rm T}}{{p}}$ 为需要估计的源信号, 则式(7)为目标函数. 定义g为二次函数G的导函数, 在y标准化前提下, 由Kunhn-Tucker条件, 目标函数的最大化问题转化为在约束 $\left\| {{{{w}}^2}} \right\| = 1$ 的条件下求解式(8):

$E\left[ {pg({w^{\rm T}}p)} \right] - \beta w = 0$

(8)

其中, $\beta = E\left[ {{w^{\rm T}}pg({w^{\rm T}}p)} \right]$ , 由牛顿法可得:

${x_{n + 1}} = {x_n} - \frac{{f({x_n})}}{{{f^/}({x_n})}}$

(9)

则FastICA算法中分离矩阵的迭代更新形式如式(10):

$w \leftarrow E\left[ {pg({w^{\rm T}}p)} \right] - E\left[ {{g^/}({w^{\rm T}}p)} \right]w$

(10)

对满足收敛条件的 $w$ 进行归一化, 得到最终的分离矩阵 ${W_{}}$ , ${W_{}}$ 左乘 ${x_1}(i)$ 即可得到最终的估计矢量信号.

2.3 心率计算

通过FastICA算法处理6路G-BVP信号, 得到盲源分离后的6路信号, 本文通过频谱分析确定脉搏波信号求得心率值. 首先通过快速傅里叶变换求得信号频谱, 计算信号在心率通频带范围内的谱峰数量及谱峰峰值, 幅值超过6路信号频谱中最大幅值80%的频率点确定为谱峰, 选择6路信号中在心率通频带范围内谱峰数最少的信号为脉搏波信号, 如果同时有两路以上信号满足谱峰数最少且相等, 则判断谱峰峰值大的信号为最终的脉搏波信号.

图8为FastICA算法处理得到的6路信号在心率通频带范围内的频谱, 图中横轴为频率点, 纵轴为幅值, ch为分离通道, 红色线条表示最大幅值的80%处, 红色圆圈表示峰值. 算法处理的G-BVP信号由一个15 s的测试视频处理得出, 心率信号通频带范围设定为0.8 Hz到4 Hz, 则频谱图中心率通频带范围对应横轴坐标12到60处. 以图8为例确定脉搏波信号, 可以看到ch1, ch3, ch4, ch5四个通道的信号存在谱峰, 由于ch4通道的信号通频带范围内频谱只有一个谱峰, 所以图8所示的分离信号中我们选择第4个通道作为脉搏波信号.

脉搏波信号求解心率有两种方法, 一种是计算时域信号峰值数即为心率值, 另一种是计算频域峰值对应的频率点, 再进一步换算成心率值. 本文选择频域计算的方法, 对脉搏波信号进行频域理想带通滤波, 滤除因为采集设备及环境引起的高频噪声. 选取0.8 Hz到4 Hz的通频带范围, 对应心率值为每分钟40次到240次. 再通过峰值检测求得信号频谱图能量最大的峰值处坐标.

$HR = 60 \times \left( {\frac{{Mr}}{t}} \right)$

(11)

根据式(11)可以求得心率值 , 其中 $t$ 为视频采集时间, $Mr$ 对应信号谱峰位置的序号. 以图8为例, 脉搏波信号为ch4, 视频采集时间为15 s, Mr为18, 按照公式计算得到心率值为72 bpm. 根据此方法计算心率简单, 节省时间, 对于噪声较多的信号, 信号峰值与实际心率值对应峰值容易有偏离, 通过本文上述步骤进行处理得到的心率信号信噪比高, 计算得到的心率值精度较高, 降低了固有的系统误差.

3 实验与结果分析 3.1 .实验数据获取

为了验证本文提出的方法, 本文在国外公开数据集DEAP上进行了实验. DEAP数据集是一个使用脑电信号, 生理信号, 以及视频进行情绪分析的多模态数据集. 该数据集记录了32个参与者受刺激情况下的脑电信号和其他的生理信号, 数据集的基本信息如表2和表3所示. 32个参与者中有22名参与者提供了视频数据, 其中有11名女性参与者, 11名男性参与者, 每名参与者包含约40个60秒的视频数据, 共计874个视频数据. 本文的实验中使用了这些视频数据用于测量心率值. 数据集原始的数据记录, 有32个.bdf文件, 每个文件有48个录制频道, 其中32个EEG通道, 12个外周通道, 3个未使用通道, 1个状态通道. 其中第39个通道包含了由接触式设备记录的BVP信号, 本文使用了第39个通道的BVP信号计算心率作为测量结果的参照标准.

对于非接触式心率测量来说, DEAP数据集并不是一个理想的数据集, 数据集提供的视频有3点对于测量结果造成较大的影响.

(1)如图9所示, 参与者面部存在电极, 电线, 眼镜的遮挡, 如图遮挡增加了特征区域检测算法的跟踪难度, 非接触式心率测量的实验要求特征区域尽可能稳定.

(2)参与者在相对较暗的环境下观看彩色视频, 进一步加大了人脸面部光照不均对检测结果的影响, 而非接触式测量要求面部光照尽可能稳定和均匀.

(3)参与者在录制过程中并不被要求静止不动, 参与者在实验过程中可以自由转动头部, 甚至可以喝水等大幅度动作, 这导致录制视频中存在较大的运动伪差问题, 而非接触式心率测量要求视频主体尽可能保持不动.

3.2 实验评价指标

为了评价本文提出的非接触式心率测量方法, 我们在DEAP数据集的874个视频中选择了300个视频进行实验测量心率值 $H{R_{\rm {nc}}}$ , 采用数据集提供的BVP数据计算心率作为真实心率值 $H{R_{\rm{gt}}}$ , 则单次实验的误差 $H{R_{\rm{err}}}$ 如式(12)所示. 本文采用了统计学指标评价所提出方法的准确性和相关性.

$H{R_{\rm {err}}} = H{R_{\rm{nc}}} - H{R_{\rm {gt}}}$

(12)

采用均方根误差 $RMS{E_e}$ 以及平均误差率 ${M_r}$ 对本文方法测量结果与真实心率值的准确性进行分析, 均方根误差和平均误差率的计算方法如式(13), 式(14)所示.

$RMS{E_e} = \sqrt {\frac{1}{N}{{\sum\limits_{n = 1}^N {\left| {H{R_{\rm {err}}}\left( n \right)} \right|} }^2}} $

(13)

${M_r} = \frac{1}{N}\sum\limits_{n = 1}^N {\frac{{\left| {H{R_{\rm {err}}}\left( n \right)} \right|}}{{H{R_{\rm{gt}}}\left( n \right)}}} \times 100\% $

(14)

采用相关性分析来非接触式测量方法与真实心率数据一致性的评价, 相关性分析是指对两个或多个具备相关性的变量元素进行分析, 从而衡量两个变量因素的相关密切程度. 采用Pearson相关系数进行计算, 如式(15)所示:

${\rho _{x \cdot y}} = \frac{{\operatorname{cov} \left( {x,y} \right)}}{{{\sigma _x} \cdot {\sigma _y}}}$

(15)

式中, $\operatorname{cov} \left( {x,y} \right)$ 表示系统测试数据和血氧仪测试数据的协方差, ${\sigma _x}$ 和 ${\sigma _y}$ 分别表示测量心率数据和真实心率数据的标准差,

3.3 实验对比

本文使用DEAP数据集的300个视频进行实验, 使用Verkruysse等人提出的单一通道结合频谱分析的方法, Poh等人提出的基于ICA的三通道盲源分离方法, 以及基于独立矢量分解的方法, 分别对300个视频进行实验, 和本文提出的多区域单通道的FastICA方法进行对比, 实验结果指标如表4所示.

从表4中可以看到, 本文所提出的非接触式心率测量方法的误差相对于Verkruysse, Poh提出的方法, 以及IVA的方法都要更小, 在相关性方面也优于上述3种方法. 表中Verkruysse和Poh的实验结果可以看到, Poh的方法的平均误差及均方根误差都要优于Verkruysse的方法, 误差率反而更低, 此处是由于Poh的方法采用盲源分离的算法处理多通道, 由于盲源分离算法输出通道没有顺序, Poh的方法在脉搏波信号的选择上采用了固定通道导致的. 本文方法在使用快速独立成分分析的基础上, 加以频谱分析确定脉搏波信号, 提高了测量结果的准确性, 有效的避免了因盲源分离算法输出通道没有一定的顺序造成的干扰.

为了验证本文方法在实际环境下的性能, 本文连续五天在室内自然光环境下进行测量, 每天测量一组数据, 每组数据包含10次测量, 时间为一分钟, 测试数据和指夹式血氧仪数据对比分析. 实验结果如表5所示.

本文采用指夹式血氧仪测量心率作为参照心率, 采用平均误差以及相关系数对每组测量数据进行分析, 可以看到每组的平均误差最大为3.8%, 表现良好, 由于每组测量数据为连续测量, 本文方法测量的心率值总在参照心率值上下浮动, 所以单组测量相关性表现较差, 相关系数最小为0.329. 本文对于5组实测数据的总体一致性进行分析, 相关系数为0.913. 验证了本文方法在实际环境下测量心率的有效性.

此外, 为了验证本文方法在弱光照, 光照不均匀的实测环境下的性能, 本文在实验过程中控制测试者房间顶部的LED灯管的打开数量, 模拟人脸受光强度, 左右脸颊受光不均匀的情况. 如图10所示, 测试者头顶共有6个灯管, 图10(a)为正常光照环境, 测试过程中6个灯管全部打开; 图10(b)为弱光照环境, 测试过程中关闭位于测试者头部前方的两个灯管, 模拟人脸受光强度较低的情况; 图10(c)和图10(d)为人脸受光不均匀环境, 测试过程中分别关闭测试者头部左前方和右前方灯管, 模拟人脸左右两侧受光不均匀的情况.

分别在正常光照, 弱光照, 光照不均3种情况下进行了实验, 每种光照条件下进行了连续20次测量, 时间为两分钟, 测试数据和指夹式血氧仪数据对比分析, 实验结果如表6所示, 可以看到3种不同光照条件下测试数据的平均误差均比较小, 最大为1.13%, 相关系数的数值大小也比较相近. 说明了本文方法在不同光照条件下的测试结果具有一定的鲁棒性, 验证了本文方法对于弱光照, 光照不均情况的光照干扰问题起到了一定的抑制作用.

4 结论与展望

本文我们提出了一种新的基于人脸视频的非接触式心率测量方法, 针对视频采集过程中运动伪差和光照不均对提取到的脉搏波信号质量的影响, 提出了多特征区域结合单一通道FastICA提取脉搏波信号的解决方案, 并在公共数据集DEAP上进行了实验分析, 结果表明本文方法在准确性和相关性方面均优于基于ICA, IVA的非接触式心率方法. 在实际非受控环境下采集数据进行测量验证, 结果表现良好.