基于多变量奇异谱分析的非接触式心率检测方法

1.本发明属于生物医学信号处理的技术领域，尤其涉及一种基于多变量奇异谱分析的非接触式心率检测方法。


背景技术：

2.目前心率的测量方法主要分为接触式与非接触式两类。接触式检测方法通过传统的测量仪器如心电图机、脉搏血氧仪等监测血容量脉冲。但由于测量时需要直接接触到受试者的皮肤，一定程度上会限制人体活动，且长时间的测量可能会引起受试者的不适，也不适用于婴儿及烧伤患者等特殊人群。非接触式测量方法也称远程光电容积脉搏波描记法(remote photoplethysmography，rppg)，利用相机捕捉因血液流动而导致的皮肤颜色的变化，从而提取心率。能够克服传统接触式心率监测的缺点，具有无创、便携、易于实现等特点。
3.在rppg技术中主要存在两大难题—光照变化和运动噪声。目前大多数研究场景都处于室内光源或自然光源条件下，光照变化比较小的情况下让受试者保持静止。rencheng song等人采用集合经验模态分解对筛选的最佳感兴趣区域的绿色通道进行分解，然后提取共同信息的方法(eemd-mcca)，该方法能够有效地减小环境光变化对视频心率提取的影响。但是eemd-mcca是两步法，eemd的使用主要是构造mcca的多通道输入集。但是在每个roi中获得的imf的数量是不同的，需要对imf进行选择和填充，以保持每个集合的通道数量相同。imf的选择过程是启发式的，可能会因为模态混合而失去有用的项。这将降低mcca算法的效果。并且在实际场景中，受试者难免会有运动，因此如何在光照变化的情况下又能消除运动噪声的干扰，从而得到准确心率也成为一个挑战。


技术实现要素：

4.本发明是为解决上述技术所存在的不足，提出了一种基于多变量奇异谱分析的非接触式心率检测方法，以期能去除环境光变化的影响和运动噪声的干扰，从而能提高非接触式视频心率检测的准确性。
5.本发明为解决技术问题采用如下技术方案如下：
6.本发明一种基于多元奇异谱分解的非接触式视频心率检测方法的特点是按照如下步骤进行：
7.步骤1：获取受试者的t帧视频数据并有重叠的划分为l份数据，且每份数据包含n帧视频图像；
8.步骤2：采用人脸检测及面部追踪法确定第l份n帧视频图像的面部感兴趣区域，将每帧视频图像中的面部感兴趣区域划分成q块子区域，1≤l≤l；
9.步骤3：逐帧计算每块子区域在第l份n帧视频图像的像素均值，并根据每块子区域的像素均值，计算每块子区域的光照强度、光照变化和信噪比作为判断的质量指标，从而根据每块子区域的质量指标，从q块子区域中选取p个最佳子区域，再从p个最佳子区域中提取
rgb通道均值信号并转化为hsv信号，提取所述hsv信号中的色度信号，记为h
l
＝{h
l,1
,h
l,2
,...,h
l,p
,...,h
l,p
}，h
l,p
为第l份n帧视频图像的第p个最佳子区域的色度信号，1＜p＜p；
10.步骤4：采用openface法逐帧检测第l份n帧视频图像的人脸特征点，并将第l份n帧视频图像的鼻尖位置特征点的运动轨迹信号记为v
l
＝{v
l,1
,v
l,2
,...,v
l,j
,...,v
l,n
}；其中，v
l,j
为第l份第j帧视频图像中鼻尖位置特征点的位置坐标；
11.步骤5：利用第l份运动轨迹信号v
l
对第l份第p个最佳子区域的色度信号h
l,p
进行自适应滤波lms处理，以去除第l份色度信号h
l,p
中的运动伪影，得到滤波后的第l份第p个色度信号且且表示第l份第j帧视频图像中滤波后的第p个色度信号的像素值，从而得到第l份滤波后的p个色度信号并作为第l个输入信号数据集；
12.步骤6：采用奇异谱分解方法将所述第l个输入信号数据集x
l
分解成若干个分量：
13.步骤6.1：设置窗口的长度为m，且m＜n/2，定义参数k＝n-m+1，按照所述窗口的长度m和参数k从所述滤波后的第l份第p个色度信号中生成维度为m
×
k的第l份第p个轨迹矩阵然后将第l份p个轨迹矩阵拼接后构成维度为pm
×
k的汉克尔轨迹矩阵
14.步骤6.2：计算第l份轨迹矩阵y
lylt
的特征值，记为λ
l,1
,...,λ
l,i
,...,λ
l,pm
，其中，λ
l,i
表示第l份轨迹矩阵y
lylt
的第i个特征值；计算第l份轨迹矩阵y
lylt
的特征值λ
l,1
,...,λ
l,i
,...,λ
l,pm
的标准正交向量，记为u
l,1
,...,u
l,i
,...,u
l,pm
，其中，u
l,i
表示特征值λ
l,i
对应的标准正交向量；计算r＝rank(y
l
)，r表示矩阵y
l
的秩，计算主成分
15.将第l份汉克尔轨迹矩阵y
l
分解为y
l
＝y
l,1
+...+y
l,i
+...+y
l,r
，其中，y
l,i
表示第l份轨迹矩阵y
lylt
第i个特征值λi对应的分解矩阵，且
16.步骤6.3：将第l份索引集{1,...,r}划分为r个不相交的子集i
l,i
＝{i},i＝1,2,...,r，i
l,i
表示第l份第i个子集；根据所述子集，对所述汉克尔轨迹矩阵y
l
进行分组，得到第l份r组矩阵{y
l,i
|1≤i≤r}，y
l,i
由步骤6.2得到；
17.步骤6.4：对第l份r组矩阵y
l,i
进行对角线求平均处理，得到第l份重构信号其中，表示第l份分组矩阵y
l,i
对角线平均处理后的一维信号；
18.步骤7：将第l份重构信号分组，每组包含四个信号，将第一组信号
作为候选脉冲集，计算第l份第一组第i个候选心率信号的主频率以及二次谐波频率的能量之比，从而选择能量之比最大的候选心率信号作为脉冲信号，利用快速傅里叶变换将所述脉冲信号转换成频域形式，从而得到所述脉冲信号的主频率f
main
，计算出第l份n帧视频图像的受试者平均心率值，从而得到第l个心率值hr
l
＝f
main
×
60；
19.步骤8：对按照步骤2-步骤7的过程，得到l份n帧视频图像的l个心率值集合{hr
l
|l＝1,2,
…
,l}；
20.步骤9：计算第l份的n帧视频图像的心率值hr
l
分别与同一受试者在第k份的n帧视频图像的心率值hrk之间的绝对误差是否小于所设定的阈值th1，若小于，则令计数值s+1；否则计数值s保持不变，其中，计数值s的初始值为0；
21.步骤10：判断计数值s是否大于所设置的阈值th2，若大于，则将第l个心率值hr
l
视为有效值并计入心率候选集hr，否则，判断下一个心率值；从而得到最终的心率候选集hr，并计算其平均值作为参考值hr
ref
；
22.步骤11：根据参考值hr
ref
与l个心率值集合{hr
l
}
(l＝1,2,...,l)
，计算第l个心率值hr
l
与参考hr
ref
之间的绝对误差hr
error,l
，如果绝对误差hr
error,l
小于等于设定的阈值th1，则第l个心率值hr
l
作为目标hr
t
，反之，则被作为异常值；
23.若第l个心率值hr
l
为异常值，则从剩余的分解信号中计算最接近参考值hr
ref
的心率值hr
n,new
作为目标hr
t
。
24.与已有技术相比，本发明有益效果体现在：
25.1、本发明选用色度信号：大多数研究选取绿色通道信号，它的优势在于所包含的心率信号强度相比于红色和蓝色通道更大。但是在光线变化的情况下，绿色通道的颜色空间不仅取决于物体的颜色，还取决于来自表面的反射光的强度。相反，色度信号不依赖于亮度，这意味着色度信号相比于绿色信号更能抵抗环境光的变化，从而降低了环境光变化对提取心率的影响。
26.2、本发明采用自适应滤波的方法去除运动噪声：采用openface方法逐帧检测人脸特征点，将鼻尖位置的运动轨迹作为运动信号，利用自适应滤波能够有效去除色度信号中由于头部运动造成的运动伪影，提升了心率测量的准确性。
27.3、本发明将面部感兴趣区域划分为多个子区域并进行最优筛选：不同的面部感兴趣区域所包含的心脏搏动的信息是相同的，相比于单个感兴趣区域提取心率信号源，该方法更强调多个感兴趣区域所包含的共同的信号源成分，对多个区域的联合分析能够更准确的提取心率。
28.4、本发明采用多元奇异谱分解的方法提取心率：多元奇异谱分解避免了传统的经验模态分解方法的模态混叠问题，可以直接分解多个信号，同时考虑它们之间的内在相关性，能够有效地减少心率信号失真率。
附图说明
29.图1为本发明方法的流程图；
30.图2a为本发明面部感兴趣区域示意图；
31.图2b为本发明筛选面部最佳感兴趣区域示意图；
32.图3a为本发明多元奇异谱分析经过筛选得到的候选心率信号1；
33.图3b为本发明多元奇异谱分析经过筛选得到的候选心率信号2；
34.图3c为本发明多元奇异谱分析经过筛选得到的候选心率信号3；
35.图3d为本发明多元奇异谱分析经过筛选得到的候选心率信号4；
36.图4a为本发明多元奇异谱分析经过筛选得到的候选心率信号1频谱图；
37.图4b为本发明多元奇异谱分析经过筛选得到的候选心率信号2频谱图；
38.图4c为本发明多元奇异谱分析经过筛选得到的候选心率信号3频谱图；
39.图4d为本发明多元奇异谱分析经过筛选得到的候选心率信号4频谱图；
40.图5a为本发明心率计算流程图；
41.图5b为本发明筛选异常值流程图。
具体实施方式
42.本实施例中，一种基于多元奇异谱分解的非接触式视频心率检测方法，如图1所示，首先获取人脸面部视频图像序列，并确定面部感兴趣区域；然后将面部感兴趣区域化分为若干子区域，根据光照强度、光照变化和信噪比筛选最优子区域，提取最优子区域的rgb信号转化为hsv信号，将鼻尖部位的运动轨迹作为运动信号；将色度信号与运动信号做自适应滤波，滤波后的色度信号作为每块子区域的输入信号；之后采用多元奇异谱分解处理所有子区域的输入信号并得到所有子区域分解信号数据集；筛选每个子区域第一个分解信号，记为候选心率信号，然后计算所有候选心率信号主频率及二次谐波频率能量占比，筛选能量占比最大的候选心率信号作为脉冲信号，并根据该信号主频值计算出心率值，然后根据心率连续性将异常值挑选出来并用正确的心率值替换。具体的说，是按照如下步骤进行：
43.步骤1：获取受试者的60s共t帧视频数据，分为l份，l＝7，且每份数据包含n帧视频图像。
44.步骤2：对于第l份n帧视频图像，采用viola-jones面部检测器得到一个矩形边框，缩放到60％用以去除非皮肤区域，如背景和头发。然后将整个roi区域划分q＝4
×
4块子区域，如图2a所示，利用kanade-lucas-tomasi算法沿着视频帧跟踪每块子区域的四个顶点。
45.步骤3：逐帧计算每块子区域在第l份n帧视频图像的像素均值，并根据每块子区域的像素均值，计算每块子区域的光照强度、光照变化和信噪比作为判断的质量指标，从q块子区域中选取p＝4个最佳子区域，如图2b所示。再从p个最佳子区域中提取rgb通道均值信号并转化为hsv信号，相比于rgb图像域到hsv图像域的转换相比，信号之间的转换速度更快，可以实时应用。提取所述hsv信号中的色度信号，记为h
l
＝{h
l,1
,h
l,2
,...,h
l,p
,...,h
l,p
}，h
l,p
为第l份n帧视频图像的第p个最佳子区域的色度信号，1＜p＜p。选择色度信号代替绿色信号是因为色度信号更能抵抗光线的变化。在光线变化的情况下，绿色通道的颜色空间不仅取决于物体的颜色，还取决于来自表面的反射光的强度。相反，色度信号不依赖于亮度。这意味着色度信号更能容忍环境光的变化。
46.步骤4：采用openface法逐帧检测第l份n帧视频图像的人脸特征点，并将第l份n帧视频图像的鼻尖位置特征点的运动轨迹信号记为v
l
＝{v
l,1
,v
l,2
,...,v
l,j
,...,v
l,n
}；其中，v
l,j
为第l份第j帧视频图像中鼻尖位置特征点的位置坐标；
47.步骤5：利用第l份运动轨迹信号v
l
对第l份第p个最佳子区域的色度信号h
l,p
进行
自适应滤波lms处理，以去除第l份色度信号h
l,p
中由于头部运动带来的运动伪影。自适应滤波在没有关于待提取信息的先验统计知识的条件下，直接利用观测数据根据某种判据在观测过程中不断的递归更新处理参数，自动跟踪统计性质的变化，以逐步逼近某一最优处理结果。这样的处理方法更符合非平稳情况的需要，也被证明是一种有效的噪声消除方法。自适应滤波后得到第l份第p个色度信号且且表示第l份第j帧视频图像中滤波后的第p个色度信号的像素值，从而得到第l份滤波后的p个色度信号并作为第l个输入信号数据集；
48.步骤6：采用奇异谱分解方法将所述第l个输入信号数据集x
l
分解成若干个分量：
49.步骤6.1：设置窗口的长度为m，且m＜n/2，定义参数k＝n-m+1，按照所述窗口的长度m和参数k从所述滤波后的第l份第p个色度信号中生成维度为m
×
k的第l份第p个轨迹矩阵然后将第l份p个轨迹矩阵拼接后构成维度为pm
×
k的汉克尔轨迹矩阵
50.步骤6.2：计算第l份轨迹矩阵y
lylt
的特征值，记为λ
l,1
,...,λ
l,i
,...,λ
l,pm
，其中，λ
l,i
表示第l份轨迹矩阵y
lylt
的第i个特征值；计算第l份轨迹矩阵y
lylt
的特征值λ
l,1
,...,λ
l,i
,...,λ
l,pm
的标准正交向量，记为u
l,1
,...,u
l,i
,...,u
l,pm
，其中，u
l,i
表示特征值λ
l,i
对应的标准正交向量。计算r＝rank(y
l
)，r表示矩阵y
l
的秩，计算主成分
51.将第l份汉克尔轨迹矩阵y
l
分解为y
l
＝y
l,1
+...+y
l,i
+...+y
l,r
，其中，y
l,i
表示第l份轨迹矩阵y
lylt
第i个特征值λi对应的分解矩阵，且
52.步骤6.3：将第l份索引集{1,...,r}划分为r个不相交的子集i
l,i
＝{i},i＝1,2,...,r，i
l,i
表示第l份第i个子集。根据所述子集，对所述汉克尔轨迹矩阵y
l
进行分组，得到第l份r组矩阵{y
l,i
|1≤i≤r}，y
l,i
由步骤6.2得到；
53.步骤6.4：为了将分组后的矩阵转换为时间序列，对第l份r组矩阵y
l,i
进行对角线求平均处理，得到第l份重构信号其中，表示第l份分组矩阵y
l,i
对角线平均处理后的一维信号；本实例中，采用上述的多元奇异谱分解方法将输入信号分解成若干分量。该方法能够充分利用不同区域输入信号的关联性并且克服了传统经验模态分解方法模式混叠的影响。
54.步骤7：将第l份重构信号分组，每组包含四个信号，将第一组信号
作为候选脉冲集，如图3a、图3b、图3c、图3d所示，本实例中提取了四个分解信号构成了候选脉冲集。计算第l份第一组第i个候选心率信号的主频率以及二次谐波频率的能量之比，，从而选择能量之比最大的候选心率信号作为脉冲信号。如图4a、图4b、图4c、图4d所示，每张图分别对应单个候选心率信号频率谱图，黑色圆圈对应峰值频率点也即主频率点。利用快速傅里叶变换将所述脉冲信号转换成频域形式，从而得到所述脉冲信号的主频率f
main
，计算出第l份n帧视频图像的受试者平均心率值，从而得到第l个心率值hr
l
＝f
main
×
60。选择第一组分解信号当作候选脉冲及是因为在前面的操作步骤中已经去除色度通道中的大部分噪声。因此，心跳信号是色度通道中最显著且相关的振荡信号。
55.步骤8：如图5a所示，按照步骤2-步骤7的过程，得到l份n帧视频图像的l个心率值集合{hr
l
|l＝1,2,
…
,l}；
56.步骤9：计算第l份的n帧视频图像的心率值hr
l
分别与同一受试者在第k份的n帧视频图像的心率值hrk之间的绝对误差，是否小于所设定的阈值th1，若小于，则令计数值s+1；否则计数值s保持不变，其中，计数值s的初始值为0；
57.步骤10：判断数值s是否大于所设置的阈值th2，若大于，则当前的hr
l
将被视为有效值并计入心率候选集hr，否则，当前的hr
l
不计入心率候选集hr，最后，确定参考值hr
ref
为心率候选集hr的平均值；
58.步骤11：获取参考值hr
ref
与l个心率值集合{hr
l
}
(l＝1,2,...,l)
，计算当前hr
l
与参考hr
ref
之间的绝对误差hr
error,l
，如果绝对误差hr
error,l
小于等于设定的阈值th1，则当前的hr
l
被视为目标hr
t
，反之则被视为异常值，若为异常值，从剩余的分解信号计算最接近参考值hr
ref
的心率值hr
n,new
作为目标hr
t
。如图5b所示，根据心率连续性得到去除异常值的目标hr
t
。
59.为了验证本发明提出的视频心率算法的鲁棒性，本发明采用公开数据集cohface和自采数据集bsipl进行算法验证。本实例通过比较两个数据集采集视频的真实心率与待测算法测量心率的误差来分析实验结果，并采用均方根误差(root mean square error，rmse)，平均绝对误差(mean absolute error，mae)，标准差(standarddeviation，sd)和相关系数这四种评价指标来评估算法鲁棒性。将本发明提出的方法与eemd-mcca方法进行了对比，结果如表1和表2所示。
60.表1 cohface两种方法得到的心率检测结果分析
61.评价指标eemd-mccamssa-h均方根误差(bpm)4.812.19平均绝对误差(bpm)2.081.25标准差(bpm)4.331.81相关系数0.910.98
62.表2 bsipl两种方法得到的心率检测结果分析
63.评价指标eemd-mccamssa-h均方根误差(bpm)6.144.00平均绝对误差(bpm)2.951.81
标准差(bpm)5.393.56相关系数0.860.94
64.从表1可以看出，在均方根误差，平均绝对误差以及标准差这三项指标上，mssa-h方法相比较于eemd-mcca有明显的改进，分别是2.62bpm，0.83bpm和2.52bpm，并且相关系数也由0.91提升到0.98，从表2结果可以所提出的mssa-h四项指标均优于eemd-mcca。从得到的结果可以看出，本方法计算得到的心率值和心率真值更接近，相对于eemd-mcca有更好的鲁棒性。
65.综上所述，本发明中提出的基于多元奇异谱分析的非接触式心率检测方法能够较为准确的准确地从视频中提取出人体心率并获取视频心率检测结果，具有良好的鲁棒性。