自适应去噪的非接触式生理参数检测方法

2021-03-09倪宗军郑秀娟

计算机工程与应用 2021年5期

倪宗军，陈辉，张昀，苏敏，郑秀娟

1.四川大学电气工程学院，成都610065

2.西安交通大学电子与信息学部，西安710049

成像式光电容积描记法（imaging Photoplethysmography，iPPG）是一种非接触式生理参数检测技术[1]，其基本原理是由于皮下浅层血管血流灌注使得皮肤的颜色发生轻微变化，因此可以通过成像设备采集的人体表面部位的视频信号获得血容量脉冲（Blood Volume Pulse，BVP）信号[2]。通过测量BVP波形的两个连续峰（或谷底）之间的时间间隔，可以从BVP信号中提取心率（Heart Rate，HR）和呼吸率（Respiratory Rate，RR），甚至可以进一步得到人体血氧、血压、心率变异性、血管微循环等生理参数，对人体健康的监护具有重要意义。但由于人体脉搏波信号十分微弱，因此利用成像式光电手段获得BVP信号极易受到干扰，噪声的来源主要包含低频基线漂移和高频噪声，分别对应运动噪声和光照变化产生的噪声，而人体脉搏波信号的主要信息处于低频和高频之间部分。

iPPG技术可以追溯到2008年，Verkruysse等人首次提出通过iPPG测量生命体征的通用方法。首先，确定暴露皮肤区域为感兴趣区域（Region of Interest，ROI），通常可以选取为面部的前额区域，接着取ROI中所有像素值的平均值作为原始BVP信号，而后对原始BVP信号进行带通滤波，保留与生理信号相对应的频率分量，然后采用快速傅里叶变换（FFT）的频谱分析方法提取心率和呼吸率对应的频率分量[3]。该通用方法对于运动及光照变化的稳健性较差。2010年，Poh等人在该通用方法基础上引入独立分量分析（ICA）方法处理面部视频的RGB三通道信号，从而降低运动噪声带来的影响，得出视频中最强的BVP分量并用于估计HR[4]。接下来，Poh等人又将时域滤波器集成到ICA方法中，从而稳健地提取BVP信号，并实现心率、呼吸率及心率变异性等相关生理参数的估计[5]。Lewandowska等人提出仅使用两个通道和应用主成分分析来计算心率，该方法可以较为有效地提取符合心肺频率的信号从而降低光照等噪声的影响[6]。在标准肤色假设下，CHROM方法被用来消除运动伪影带来的影响，该方法在一定程度上提高了iPPG方法的运动鲁棒性[7]。雷恒波等人利用经验模态分解法（Empirical Mode Decomposition，EMD）将BVP信号分解成一系列不同频率成分的本征模态函数（Intrinsic Mode Function，IMF），最后通过频域分析法得到心率[8]。EMD去噪的基本原理是从分解的IMF分量中去除高频噪声引起的IMF分量。该分量基本不含有效信号成分，去除该分量后，用剩余的IMF分量重构原信号，即达到了去噪的效果。但是，该方法只是去除了部分噪声，其他IMF分量中仍残留部分噪声，而且高频分量中也会含有有效信息。另外EMD还存在着其他问题，如混叠模态、端点效应和停止条件等等。而另外一种信号分析方法-变分模态分解（Variational Mode Decomposition，VMD）方法，模态分量数目K是可以根据模态分量瞬时频率的均值来人为设定最优的，在频域上能自适应地分解出各中心频率对应的有效成分，对低频段的特征提取具有更好的准确性和稳定性，它克服了EMD、经验小波变换等模式中存在的模态混合、不能正确地消除附加噪声以及确定本征模态函数个数的问题，具有很好的噪声鲁棒性和降噪效果。但该方法仍然存在各IMF中含有残留噪声的问题，如何去除各IMF中存在的噪声，又保留有效信号是精确提取BVP信号的关键，而且上述方法获取的生理参数估计的准确性都受到重构的BVP信号质量的影响。

因此，本文基于变分模态分解法提出了一种新型的联合能量熵和信号质量检测的自适应生理参数检测方法，可以很好地分解出BVP信号的高频噪声和低频基线漂移噪声，实现在自然条件下的人脸视频中稳健地估计心率及呼吸率。

1 方法

本文提出了一种基于VMD的新型生理参数测量的方法，共由三个主要步骤构成，分别为原始BVP信号获取、基于变分模态分解的能量熵阈值自适应去噪以及基于信号质量检测的自适应生理参数估计。

1.1 原始BVP信号获取

获取原始BVP信号，首先需要选择受试者脸上暴露的皮肤区域，即面部ROI。首先使用MTCNN[9]方法检测人脸，该人脸检测器采用CNN方法，利用非最大抑制[10]方法寻找最佳的人脸区域，与目前流行的其他人脸检测方法相比，具有更高的准确性和鲁棒性。在检测到人脸后，使用CE-CLM[11]方法来形成68个人脸特征点，接着通过选取合适的人脸特征点来包含鼻子和额头（剔除嘴角和眼睛区域）的区域作为面部ROI。

为了有效地提高运动鲁棒性，降低运动带来的噪声影响，在获得的面部ROI基础上，计算视频RGB三通道的平均值，并采用CHROM[7]方法，由RGB三通道的信号组合得到原始的BVP信号。

1.2 变分模态分解与能量熵阈值自适应去噪

1.2.1 变分模态分解的基本原理

Dragomiretskiy等提出了一种可变尺度的信号分解方法[12]，该方法将信号f(t)分解为K个本征模态分量uk(t)，定义本征模态函数uk(t)为一个调幅-调频信号，如式（1）所示：

（1）通过Hilbert变换，构造每个uk(t)的解析函数以获得相应的边际谱。

（2）通过指数混合调制到以wk为中心频率的频带上，将每个模型uk(t)的频谱转移到“基带”。

（3）由解调信号的梯度平方和（L2范数）最小来估计各IMF分量的带宽。即：

（4）通过引入惩罚因子α和Lagrange算子λ(t)，将有约束的变分问题转化为无约束的变分问题：

（5）采用乘法算子交替方向法（ADMM）解决上述变分问题，迭代优化uk+1、ωk+1、λk+1求得扩展拉格朗日表达式的“鞍点”。采用傅里叶变换，求取二次优化问题的频域解更新公式为：

1.2.2 能量熵阈值自适应去噪

VMD分解算法将信号分解为具有一定带宽的、频率由高到低的K个模态分量。由于BVP信号中光照随机噪声主要分布在高频部分，运动伪影噪声则主要分布在低频部分，而心肺信号主要信息处于低频和高频之间部分，因此需要将K个模态分量中的低频部分运动伪影的噪声去除，高频部分的光照噪声去除，将剩余的IMF分量重构，即可得到去噪后的信号。

首先，依据各IMF总体信息与BVP信号的相关性剔除相关性差的分量，接着根据剩余IMF中有效信息和噪声对BVP信号的相关性不同的特性，对各IMF分量中的有效信息进行检测与定位，将有效信息引起的系数剔除，噪声引起的系数保留，得到新的分量系数;再将保留的噪声IMF分量系数分成若干子区间，分别计算各子区间的噪声能量熵，能量熵最大的区域，可以认为是全部由噪声引起的，只要计算出该区域系数的能量，再利用噪声的能量构造阈值，对各原IMF分量系数进行阈值处理，剔除每个IMF分量中噪声的影响。下面给出改进算法的关键计算步骤。

先计算各IMF分量与原始信号的相关系数，定义rk(n)为第k个IMF分量与f(n)的相关系数：

接着将剩余的各IMF分量uk(n)分成t等份，定义rk(i)为第k个IMF分量uk(t)与f(t)的相关系数：

为使相关系数与IMF分量系数具有可比性，定义规范化相关系数Rk(n)为：

其中，uk(i)为第k个IMF分量中的第i段信号，通过各IMF分量的规范化相关系数与各IMF分量中的每段成分的相关系数的比较，得到各IMF分量中噪声信号的系数位置。

将IMF中的各个噪声区间提取出来，根据信号在时频域具有能量守恒特性，设每个噪声区间长度为M，定义第k个IMF分量第p个噪声区间的能量为：

将uk(p)分成l等份，每个小区间的能量记为Eki(p)，采样点数为则有：

第k个IMF分量中第p个噪声区间的第i个子区间的能量Eki(p)在该分量的总能量Ek(p)中的概率为：

则第i个子区间对应的能量熵为：

得到第k个IMF分量中第p个噪声区间的能量熵序列：

搜索第k个IMF分量的熵值最大子区间，为：

求第k个IMF分量第p个噪声区间的第m个子区间的IMF分量系数的平均值：

定义第k个IMF分量第p个噪声区间的阈值Tk(p)为：

其中，σk(p)为第k个IMF分量的第p个噪声区间的熵值最大子区间的系数的平均值作为噪声方差。根据得到的阈值Tk(p)，利用软阈值函数分别对各IMF分量中的每个噪声区间的分量系数进行阈值处理，得到新的IMF分量系数，接着利用新的IMF分量系数重构信号，得到去噪后的信号。如图1所示，将原始的BVP信号经过变分模态分解后，得到5个IMF分量，然后通过相关性检测，剔除第一个明显属于运动伪影噪声的IMF分量，起到了去趋势化的效果，接着通过对每个保存下来的IMF分量进行更细致的相关性分析，对属于噪声区间的分量进行噪声能量熵阈值分析，将得到的能量熵阈值重构噪声区间的分量，从而得到更准确的IMF分量，最后对重构的BVP信号进行相应范围的滤波即可得到理想的心率和呼吸率的时域信号。故该噪声能量熵阈值的确定方法可根据信号中噪声的能量特征自适应地去除IMF分量的噪声，重构出理想的BVP信号。

1.3 自适应生理参数估计

1.3.1 基于方差表征序列的信号质量检测

从原始BVP信号中提取得到去噪后的心率信息后，为了继续检测该信号的质量，采用Pang等人提出的方差特征序列（Variance Characterization Series，VCS）检测方法[13]，进行心率信号的质量检测。步骤如下：

步骤1找出心率信号中的所有局部最大值Mi,i=1,2,…,以及局部最小值mi,i=1,2,…。

步骤2计算每一个最大局部最大值Mi和接下来七个局部最大值的方差σMi，然后再计算每一个局部最小值mi和接下来七个局部最小值的方差σmi。

步骤3对于每个方差，设：

对于每个最大值和最小值，可以得到时间序列x(t)的两个VCS。如图2所示，当不同心率信号的波形变化时，VCS的变化情况。绿色三角代表心率信号的极大值，红色圆圈代表心率信号的极小值。对于每个最大值的VCS，基于两个准则来判断心率信号的质量：

（1）δMi的值远大于THAM。

（2）(i+4)th的最大值和ith的最大值之间的距离远小于distFAR或者(i+2)th的最大值和ith的最大值之间的距离远大于distNEAR。

THAM的设定由经验给定，而distFAR和distNEAR的设定由之前估计的心率Epre和采样频率p给定，即：

图1 提取心率和呼吸率成分的时间序列步骤

图2 不同心率成分信号以及相对应的方差表征序列

1.3.2 自适应生理参数估计

VCS质量检测合格表明已经得到较理想的BVP信号，该信号中存在强烈的心肺频率成分，接着分别对符合心率和呼吸频率的心肺信号进行三阶的巴特沃斯带通滤波提取相应的频率成分，滤波范围分别为（0.7～4 Hz）和（0.1～0.7 Hz），然后使用快速傅里叶（FFT）转换到频域上分析，最高的频率成分即是对应的心肺频率fhr、fbr。最终心率（HR）和呼吸值（RR）可估计为：

当VCS质量检测不合格时，表明心肺信号中还存在较强烈的噪声频率，这个时候选择时域的分析方法来获取心率和呼吸率值，在心肺信号上使用数峰值的方法来近似得到心肺速率的估计值etime。为了更进一步提高准确性，如图3所示，接着进行了频域跟踪的方法。接下来是具体的步骤：

（1）将心肺信号进行三次样条插值处理。由于动作或光照带来的影响是短暂的，进行插值处理后，可以加强心肺频率成分，从而削弱噪声带来的影响。

（2）使用FFT得到插值后的心肺信号的频谱分布，找到前3个最大频率成分。

（3）将这3个估计的频率fest与前一次估计的fpre进行比较，选择与前一估计最接近的频率为频域跟踪结果的估计epre。然后与时域分析的估计值etime比较，如果两种方法的估计值相差在10 beat/min（每分钟心跳次数）以内，则取两者的平均值作为最终的心肺频率估计efinal。如果相差大于10 beat/min，则说明这段心肺信号的质量太差，还存在较多噪声，则选用机器学习的方法来预测结果eML。本文选用了背景ROI和人脸ROI的80个特征，其中包括频域中前5个最大频率成分、时域信号中波峰到波谷的距离、相邻波峰、波谷的间距、ICA变换后的结果等等，最终将80个特征用支持向量机做回归预测作为最终的估计efinal。

图3 频谱跟踪算法

2 实验结果及分析

2.1 实验数据集和设置

选用公开的DEAP数据集[14]进行该方法性能验证并与目前已有方法进行比较。DEAP数据集是一个公共多模态数据库，包括22名参与者（22～35岁）的面部视频。各参与者均拍摄了40段长度为1 min的面部视频，同时还通过指夹式设备和呼吸带记录了心率及呼吸率数据，作为方法性能评估的金标准。视频拍摄的分辨率为720×576，帧率为50 f/s。参与者被要求静坐电脑前观看“情绪诱导实验”中的各种视频，随着视频的变化会改变参与者脸部的光照，模拟了自然条件下的光照变化，因此DEAP数据集上的实验结果可以反映各方法在光照变化时的性能效果。实验从22名参与者中选取金标准数据质量高（没有过多的动作伪迹干扰，以保证金标准的准确性）的数据进行后续的对比实验。

为了进一步进行实际场景中的头部运动和光照变化的影响实验比较，本文还采集了各种不同实际场景下的实验数据进行方法的性能比较。使用罗技C920网络摄像头默认设置拍摄受试面部视频，使用消费级鱼跃血氧仪（YX303）心率结果作为金标准（由于血氧仪只含心率结果，故没有进行呼吸率结果的比较）。采集了总共20个受试者在5种不同场景下的各5 min视频。场景1，受试者在室内电脑前坐立，距离摄像头前0.5 m处，可进行点头、摇头、微笑等动作，光源为普通LED灯，面部位置光强为300 lux；场景2，受试者在室内坐立，距离摄像头前0.5 m，可进行点头、摇头、微笑等动作，光源为普通LED灯，调节面部处光强从10 lux到750 lux变化；场景3，受试者在室内站立，距离摄像头前0.5 m处，进行身体左右摇晃、慢速走动等动作，光源为普通LED灯，面部位置处光强为300 lux；场景4，受试者站立在户外街道，距离摄像头前0.5 m处，进行点头、摇头、微笑等动作；场景5，受试者站立在户外广告牌下，距离摄像头前0.5 m处，进行点头、摇头、微笑等动作。

本文中实现了基于面部视频的实时生理信息采集系统软件，所有实验在个人笔记本电脑（微星GP63，Win10操作系统，内存8 GB，处理器i7-8750H，主频2.2 GHz）上运行，本文方法运行速度可以平均稳定在25 f/s。

2.2 实验结果对比

BVP信号处理中最常用的方法有离散时间傅里叶变换（DTFT）、离散小波变换（DWT）、峰值计数（峰值计数）、独立分量分析（ICA）、自回归分析（AR）和经验模态分解（EMD）。快速傅里叶变换（FFT）复杂性较低，独立分量分析（ICA）方法的使用频率很高，自回归分析（AR）方法比较新颖，经验模态分解（EMD）与变分模态分解原理相似，所以本文一共选择了FFT[2]、ICA[4]、AR[15]和EMD[8]方法来进行性能比较。本文在整体方法流程上做了很多改进，因此为了说明其他非BVP因素（如ROI获取、人脸跟踪方法、信号平滑以及频率跟踪算法）的改进带来的提升，将整个方法与其他论文提出的整体方法进行比较，本文选取的ROI区域如图4所示，其中绿色空心圆表示选取的部分特征点，蓝色区域是选取的粗略ROI，红色矩形是精细ROI，包括鼻子和前额区域。

图4 面部ROI选择结果示例

为了公平比较，ICA、FFT、AR和EMD中用到的参数均参照原论文设置为最优值，在处理不同视频时，所有参数也保持不变。视频时间窗都选择30 s（包括25 s的重叠窗），最终比较结果如表1～3所示。

2.3 结果讨论

由表1和表2结果可见，ICA方法的性能优于FFT方法，这是因为ICA方法增加了多重时域滤波，使得时间序列更加稳定。而傅里叶变换是一种经典的频谱分析方法，在采样频率稳定，且数据长度足够的情况下，生理参数估计可以取得较好的精确度。然而在使用傅里叶变换分析低频段频谱时，其频率分辨率是不足的，而且实际应用中通常含有光照、表情动作等因素的影响，会在获取的原始BVP信号中引入了非高斯白噪声，将带来多个伪峰值，尤其是运动噪声带来的频率峰值。通常运动干扰信号的频率成分会超过BVP信号中心肺频率成分。在这种条件下，单纯地将最高峰值的频率点作为心肺频率，其估计精度将很难保证。AR[15]方法虽然可以减少人工光照的影响，但是对于一些运动伪影的干扰则不能很好地抑制。而ICA方法是假定观测到的信号为盲源的线性混合，而实际的环境光噪声和其他运动噪声可能是非线性和时变的，特别是在运动过程中，因为运动引起的血容量的生理变化也可能是非线性的。此外，ICA假设源信号的数量等于RGB通道的数量，但潜在噪声源的数量是不同的，且在实际应用难以确定。因此，ICA并不总是能正确地分解RGB通道来提取所需的心肺信号的频率。而EMD方法是通过去除属于噪声的IMF分量进行降噪，但存在着模态混叠和分解数目不可控的缺陷，如何选择去除合适的IMF分量是重构出精确的BVP信号的关键，而且之前讨论过，去除的IMF分量里也包含着有效的信息，而剩余的IMF分量中还存在着其他噪声，这是EMD方法准确度不是很理想的关键原因。从表3可知，当在更复杂的现实场景下（包含变化的光照、头部运动），这些方法达不到满意的估计结果。主要原因是，已有方法在人脸ROI获取时对动作的跟随性能不高，随着人脸转动，会出现人脸跟踪失败，进而无法得到理想的BVP信号。而且，在已有方法中也缺乏对复杂噪声的处理，如图5所示，原始BVP信号中包含较强烈的伪影噪声时，可能会出现运动伪影带来的峰值高于心率成分的峰值，其中因AR方法与FFT方法一样使用绿色通道，故未标出。ICA方法在一定程度上减弱了运动伪影带来的影响，而EMD方法通过选取合适的本征模态函数也可以较大程度减少运动伪影的噪声，但都还残留着部分噪声。本文的方法简化了环境光的过滤，着重在ROI的精确选取上，以及应用CHROM方法可以对运动噪声起到一定的抑制作用。如图6所示，当同时出现强烈的运动伪影及光照变化的噪声时，其他方法的频谱中会残留大量的伪峰，而通过应用基于VMD方法的能量熵阈值自适应去噪方法，具有更好的抗噪性，能够重构出更精确的BVP信号，而且根据重构出的BVP信号的质量通过VCS信号质量分析可以自适应地选择不同的分析方法以及提出使用的频谱跟踪算法中机器学习的方法利用了背景ROI和面部ROI的信号得到光照和运动噪声特征，可以更准确地估计出心肺信号频率。然而本文方法还不能有效地去除因剧烈运动带来的噪声影响。因此，在未来的工作中，将提高本文方法的健壮性，优化方法流程，使其适应更多的实际场景，从而提高方法在不同情况下的准确性。

表1 DEAP数据集下不同方法的心率比较结果

表2 DEAP数据集下不同方法的呼吸率比较结果

表3 自采数据集下不同方法的心率比较结果

3 结束语

本文提出了一种基于人脸视频测量生理信息的方法，通过改进ROI检测方法，可以自适应地精准选择ROI。针对直接去除低频或高频分量的EMD去噪方法中，不能有效区分各IMF分量中的有效信息和噪声信息的问题，本文提出了基于变分模态分解的相关能量熵阈值自适应去噪方法，该方法根据噪声的能量熵自适应地确定各IMF分量的阈值，重构出更精确的BVP信号。相比传统单一分析方法如FFT、数峰值等等，本文还使用了基于信号质量的自适应分析方法，有效地提升了应用范围，最后为了提高测量的精度。同时，通过采用频谱跟踪算法，可以更准确地估算心肺频率值。在公开的数据集和自采数据集上进行的性能验证实验，证明了本文方法在生理参数估计精度较已有方法有较大提升，并且可以更好地适用于实际环境中。