CN113706413B - 一种获取眼底光电体积描记信号的方法 - Google Patents
一种获取眼底光电体积描记信号的方法 Download PDFInfo
- Publication number
- CN113706413B CN113706413B CN202110987651.4A CN202110987651A CN113706413B CN 113706413 B CN113706413 B CN 113706413B CN 202110987651 A CN202110987651 A CN 202110987651A CN 113706413 B CN113706413 B CN 113706413B
- Authority
- CN
- China
- Prior art keywords
- signal
- rppg
- matrix
- noise
- fundus
- Prior art date
- Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
- Active
Links
- 238000000034 method Methods 0.000 title claims abstract description 35
- 238000013186 photoplethysmography Methods 0.000 title claims abstract description 32
- 239000011159 matrix material Substances 0.000 claims abstract description 61
- 238000001228 spectrum Methods 0.000 claims abstract description 18
- 238000012847 principal component analysis method Methods 0.000 claims abstract description 9
- 230000004927 fusion Effects 0.000 claims description 13
- 238000005457 optimization Methods 0.000 claims description 8
- 238000005316 response function Methods 0.000 claims description 8
- 208000006550 Mydriasis Diseases 0.000 claims description 4
- 230000008602 contraction Effects 0.000 claims description 3
- 238000012935 Averaging Methods 0.000 claims description 2
- 230000017105 transposition Effects 0.000 claims description 2
- 210000003491 skin Anatomy 0.000 description 7
- 210000004204 blood vessel Anatomy 0.000 description 6
- 238000005259 measurement Methods 0.000 description 5
- 239000008280 blood Substances 0.000 description 4
- 238000012545 processing Methods 0.000 description 4
- 210000004369 blood Anatomy 0.000 description 3
- 230000000737 periodic effect Effects 0.000 description 3
- 210000001525 retina Anatomy 0.000 description 3
- 230000002123 temporal effect Effects 0.000 description 3
- 210000001519 tissue Anatomy 0.000 description 3
- 206010052428 Wound Diseases 0.000 description 2
- 208000027418 Wounds and injury Diseases 0.000 description 2
- 230000009286 beneficial effect Effects 0.000 description 2
- 230000036772 blood pressure Effects 0.000 description 2
- 210000001508 eye Anatomy 0.000 description 2
- 238000005286 illumination Methods 0.000 description 2
- 230000002911 mydriatic effect Effects 0.000 description 2
- 238000000513 principal component analysis Methods 0.000 description 2
- 208000025865 Ulcer Diseases 0.000 description 1
- 210000005252 bulbus oculi Anatomy 0.000 description 1
- 230000000747 cardiac effect Effects 0.000 description 1
- 238000000354 decomposition reaction Methods 0.000 description 1
- 230000007547 defect Effects 0.000 description 1
- 210000004207 dermis Anatomy 0.000 description 1
- 238000003745 diagnosis Methods 0.000 description 1
- 201000010099 disease Diseases 0.000 description 1
- 208000037265 diseases, disorders, signs and symptoms Diseases 0.000 description 1
- 239000003814 drug Substances 0.000 description 1
- 229940079593 drug Drugs 0.000 description 1
- 230000002526 effect on cardiovascular system Effects 0.000 description 1
- 230000000694 effects Effects 0.000 description 1
- 210000002615 epidermis Anatomy 0.000 description 1
- 238000000605 extraction Methods 0.000 description 1
- 230000035876 healing Effects 0.000 description 1
- 230000031700 light absorption Effects 0.000 description 1
- 238000000691 measurement method Methods 0.000 description 1
- 238000012986 modification Methods 0.000 description 1
- 230000004048 modification Effects 0.000 description 1
- 238000012544 monitoring process Methods 0.000 description 1
- 210000001328 optic nerve Anatomy 0.000 description 1
- 230000003287 optical effect Effects 0.000 description 1
- 210000001747 pupil Anatomy 0.000 description 1
- 230000010344 pupil dilation Effects 0.000 description 1
- 238000011160 research Methods 0.000 description 1
- 230000004256 retinal image Effects 0.000 description 1
- 238000004088 simulation Methods 0.000 description 1
- 238000006467 substitution reaction Methods 0.000 description 1
- 231100000397 ulcer Toxicity 0.000 description 1
Classifications
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
- G06T5/00—Image enhancement or restoration
- G06T5/70—Denoising; Smoothing
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
- G06T2207/00—Indexing scheme for image analysis or image enhancement
- G06T2207/10—Image acquisition modality
- G06T2207/10016—Video; Image sequence
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
- G06T2207/00—Indexing scheme for image analysis or image enhancement
- G06T2207/20—Special algorithmic details
- G06T2207/20092—Interactive image processing based on input by user
- G06T2207/20104—Interactive definition of region of interest [ROI]
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
- G06T2207/00—Indexing scheme for image analysis or image enhancement
- G06T2207/30—Subject of image; Context of image processing
- G06T2207/30004—Biomedical image processing
- G06T2207/30041—Eye; Retina; Ophthalmic
Landscapes
- Physics & Mathematics (AREA)
- General Physics & Mathematics (AREA)
- Engineering & Computer Science (AREA)
- Theoretical Computer Science (AREA)
- Eye Examination Apparatus (AREA)
- Measuring Pulse, Heart Rate, Blood Pressure Or Blood Flow (AREA)
Abstract
本发明公开了一种获取眼底光电体积描记信号的方法,包括如下步骤:S1、从采集的眼底视频中提取原始rPPG信号,并确定rPPG信号模型;S2、对原始rPPG信号,根据信噪比或者信号幅度预先进行第一次去噪;S3、对第一次去噪后的rPPG信号进行第二次去噪:使用稳健主成分分析法对去噪问题进行建模,利用矩阵因子分解去除外部噪声;S4、对步骤S3去除外部噪声后的rPPG信号,通过估计心跳信号的稀疏频谱来去除内部噪声;S5、对步骤S4去除内部噪声后的rPPG信号融合时间窗,得到随时间变化的rPPG信号。
Description
技术领域
本发明涉及计算机视觉以及图像处理领域,具体涉及一种从采集的眼 底图像中获取眼底光电体积描记信号的方法。
背景技术
PPG(Photoplethysmography,光电体积描记)是一种检测组织中微血 管血容量变化的非侵入性光学技术,利用此技术获得的PPG信号可用来计算人 体某些重要生理指标例如血压、心率等。此种方法相比于传统方法(例如使用袖 带血压计测量血压,使用心电图仪测量心率等)更加地便捷、高效,而且能够实现实时连续监测,有助于尽早发现与心血管有关的疾病,使病人得到及时治疗。
目前PPG信号测量技术可大体上分为接触式和无接触式两类。接触式 方法为了进行精确测量,需要将光电探测器直接置于皮肤表面,当光源照亮组织 后,光电探测器测量一段时间内由于血容量改变导致的光吸收的变化。非接触式 方法也称为远程光电体积描记法(remote photoplethysmography,rPPG),其基 于血液在血管中流动导致皮肤反射的光强发生改变,从捕获的被采集对象的视频中提取出PPG信号。
接触式PPG信号采集方法由于需要将光电探测器直接置于皮肤表面, 这限制了其在伤口诊断(例如烧伤、溃疡、创伤等)和皮肤愈合评估或需要自由运动情况下的实用性。而非接触式方法采集的是光照射于皮肤后反射的光的强度 变化,这种改变会受到来自皮肤表皮层和真皮层的干扰,导致测量结果不准确。
发明内容
本发明的主要目的在于克服上述现有技术的不足,提出一种非接触式 采集眼底图像并从中提取出PPG信号的方法,解决现有的非接触式方法测量结 果不准确的技术问题。
为解决上述技术问题,本发明采用以下技术方案:
一种获取眼底光电体积描记信号的方法,包括如下步骤:S1、从采集 的眼底视频中提取原始rPPG信号,并确定rPPG信号模型;S2、对原始rPPG信 号,根据信噪比或者信号幅度预先进行第一次去噪;S3、对第一次去噪后的rPPG信号进行第二次去噪:使用稳健主成分分析法对去噪问题进行建模,利用矩阵因子分解去除外部噪声;S4、对步骤S3去除外部噪声后的rPPG信号,通过估计 心跳信号的稀疏频谱来去除内部噪声;S5、对步骤S4去除内部噪声后的rPPG信号融合时间窗,得到随时间变化的rPPG信号。
进一步地,前述获取眼底光电体积描记信号的方法还包括:利用眼底 相机在非散瞳情况下采集眼底视频。
进一步地,步骤S1具体包括如下步骤:S11、对采集的眼底视频的每 一帧图像进行区域划分,并对感兴趣区域进行跟踪;在每一帧图像中,对每个区 域的像素强度进行空间平均,以获得原始rPPG信号;S12、计算每个区域随时 间变化的平均像素强度,并在原始rPPG信号中减去该平均像素强度,然后使用带通滤波器将信号的频率范围限制在30~300bpm,该频率范围包含了感兴趣的心 脏信号;S13、对每个区域随时间变化的平均像素强度,利用心跳信号、通道响 应函数和通道噪声建立关于时间的信号模型,并将该信号模型改写为向量形式,得到所述rPPG信号模型。
更进一步地,步骤S11中采用KLT跟踪器和RANSAC算法对每一帧眼 底图像进行感兴趣区域跟踪。
更进一步地,步骤S13包括:对第j个区域随时间t变化的平均像素强 度pj(t)建模如下:
pj(t)=hj(t)*yj(t)+nj(t) (1)
其中,*是线性卷积算子,yj(t)表示在第j个区域观测到的心跳信号,hj(t)和nj(t)分别表示通道响应函数和通道噪声;由于已知心跳信号在频域中是稀疏的,将式 (1)改写为向量形式,得:
pj=hj*F-1xj+nj (2)
其中,F-1表示傅里叶变换F的逆变换,F是大小为W的一维离散傅里叶变换, W为时间窗长度;xj为心跳信号yj的稀疏频谱,且有xj∈CW,yj∈RW,C表示 复数集合,R表示实数集合;忽略通道响应函数的影响,得到所述rPPG信号模 型:
pj=F-1xj+nj (3)。
进一步地,步骤S2具体包括:对原始rPPG信号,将信噪比低于预设 信噪比阈值的区域或者信号最大幅值高于预设幅度阈值的区域视为噪声过多区 域并予以丢弃,实现所述第一次去噪。
更进一步地,所述预设信噪比阈值为θSNR=0.2,所述预设幅度阈值为 θamp=16。
进一步地,步骤S3具体包括如下步骤:
S31、将第一次去噪后的rPPG信号堆叠成矩阵,得到rPPG矩阵P;
S32、将rPPG矩阵P建模为:
P=Y+N=Y+E+S (4)
其中,Y为包含心跳信号的低秩矩阵;N为噪声矩阵,包含内部噪声E和稀疏外 部噪声S;
S33、使用稳健主成分分析法建立如下优化问题:
其中,||Z||*=∑kσk(Y)表示矩阵Z的核范数,其等于矩阵Z的奇异值σk之和; 矩阵S的l1范数||S||1=∑t,j|S(t,j)|,其等于矩阵S各项元素S(t,j)的绝对值之和; 参数γ控制被划分到噪声分量中的信号能量的相对比例,取值0.02~0.1;
S34、将矩阵Z分解为两个因子,令Z=LMT,其中L和M为两个低 秩因子,L∈RW×r,M∈RK×r,W表示时间窗的长度,K表示每一帧图像被划分成 K个区域,r表示秩估计参数且r<W,K;T表示矩阵转置;之后||Z||*可被表示成 如下形式:
其中||·||F表示对其内部矩阵求弗罗贝尼乌斯范数;将式(6)代入式(5)进行求解, 实现从rPPG矩阵P中消除稀疏外部噪声S。
更进一步地,步骤S4具体包括:
S41、将消除稀疏外部噪声之后的输出信号zj建模为zj=F-1xj+ej,对 应的矩阵形式为:
其中,ej为内部噪声信号,I为单位矩阵;
S42、定义如下优化问题:
其中,λ、μ均为正则项系数,λ的取值范围为0.02~0.3,μ的取值范 围为0.15~1.0;矩阵X的l2,1范数定义为:
S43、使用迭代阈值收缩算法求解式(8),以计算出矩阵X和矩阵E, 然后去除内部噪声E。
更进一步地,步骤S5融合时间窗包括:
S51、构建加权平均时间融合窗口其中,/>为滑动窗口,Qm表示前一个窗口中不存在的新的rPPG信号,Qo表示当前窗口 中与前一个窗口重叠的rPPG信号的比例,/>X表示前一窗口的滤波输出,/>是/>存在于当前窗口的部分,α为融合时间窗的阈值,取值范围为0.01~0.05;
S52、再次使用稳健主成分分析法去除时间融合窗口中的噪声,并估 计出所有时间窗口的稀疏频谱;
S53、用步骤S52估计出的稀疏频谱,结合时间窗信息,获得与采集的 眼底视频时间等长的随时间变化的rPPG信号。
本发明前述技术方案提供的获取眼底光电体积描记信号的方法,其有 益效果包括:本发明利用眼球是由透明组织构成且血管靠近主要的反射层(即视 网膜和视神经)的表面这一特点,采集眼底视频图像作为提取眼底光电体积描记 信号的数据基础,可避免皮肤层对光电体积描记信号的干扰,为后续提取光电体 积描记信号提供更可靠的数据基础;在此基础上,本发明可直接从眼底图像中提取到更加准确的眼底光电体积描记信号。无需将光电探测器直接置于皮肤,因此 适用领域更广,信号采集更加便利。
附图说明
图1是本发明实施例获取眼底光电体积描记信号的方法流程图;
图2是本发明实施例利用稳健主成分分析法对原始rPPG信号进行去噪 的流程图。
具体实施方式
下面结合附图和具体的实施方式对本发明作进一步说明。
本发明实施例提供一种获取眼底光电体积描记信号的方法,旨在从采 集的眼底视频图像中提取出眼底光电体积描记信号。首先,可利用一款能够在非 散瞳情况下采集眼底视频的眼底相机,采集一定时长的眼底视频,然后以采集的 眼底视频作为数据基础来从中提取眼底光电体积描记信号。其中,能够在非散瞳 情况下采集眼底视频的眼底相机,其可以在45°视野下拍摄视网膜的非散瞳图 像。为了在非散瞳情况下拍摄到眼底的高质量图像,本发明一实施例采用近红外 光LED照明的相机附件,这些LED以环形排布安装在红外摄像机附件上。在该 眼底相机附件中,使用940nm的红外光源,在此基础上,可为摄像机安装带宽 为10nm的带通滤波器。然后,使用Zemax软件进行射线路径的模拟,以确保 LED被放在恰当的位置。这样做可以减少角膜反射,同时提供高视野。采集眼 底视频时,相机的红外(940nm)LED照亮视网膜而不收缩瞳孔,这使设备与眼 睛对齐并聚焦于视网膜图像,该眼底相机不需使用药物进行瞳孔扩张,即可采集到清晰的眼底视频图像。
本发明实施例提供了从采集的眼底视频图像中提取出眼底光电体积描 记信号的方法,图1为该方法的流程图。参考图1,该方法包括如下步骤S1~S5:
步骤S1、从采集的眼底视频中提取原始rPPG信号,并确定rPPG(remote photoplethysmography,远程PPG)信号模型。具体包括S11~S13:
S11、对采集的眼底视频的每一帧图像进行区域划分,并对感兴趣区域 进行跟踪;在每一帧图像中,对每个区域的像素强度进行空间平均,以获得原始rPPG信号。具体而言,对于每一帧眼底图像,可以划分为48个区域,然后使用Kanade-Lucas-Tomasi(KLT)跟踪器和RANSAC算法来对感兴趣的区域进行跟踪,以观察每个区域的像素强度随时间的变化规律,在此基础上,以对某个区域j在 某一时刻t的pj(t)、hj(t)、nj(t)、yj(t)等值进行研究。应当理解的是,每一帧图像 划分的区域数量不限于48,也可为其它数值,本发明对此不作限制,48仅为一 种参考的实施方式。
S12、计算每个区域随时间变化的平均像素强度,并在原始rPPG信号 中减去该平均像素强度,然后使用带通滤波器将信号的频率范围限制在 30~300bpm,该频率范围包含了感兴趣的心脏信号。在此我们使用了一个长度为 10秒的时间窗,并使用了一个只有10帧来自新的时间窗的时间窗的重叠(0.33 秒用于每秒30帧(fps)录制的视频)。
S13、对于每个区域j∈{1,…,48},所测量的平均像素强度pj(t)是一个一 维时间序列信号,其中t∈{1,…,10}是长度为10的时间窗口内的时间视频帧索引。 对每个区域随时间变化的平均像素强度pj(t),利用心跳信号、通道响应函数hj(t) 和通道噪声nj(t)建立关于时间的信号模型:
pj(t)=hj(t)*yj(t)+nj(t) (1)
然后将该信号模型改写为向量形式:
pj=hj*F-1xj+nj (2)
其中:*是线性卷积算子;yj(t)表示在第j个区域观测到的心跳信号;hj(t)和nj(t)分别表示通道响应函数和通道噪声;F-1表示傅里叶变换F的逆变换,F是大小 为W的一维离散傅里叶变换,W为时间窗长度;xj为心跳信号yj的稀疏频谱, 且有xj∈CW,yj∈RW,C表示复数集合,R表示实数集合。忽略通道响应函数 的影响,得到如下的rPPG信号模型:
pj=F-1xj+nj (3)。
步骤S2、对原始rPPG信号,根据信噪比或者信号幅度预先进行第一 次去噪。所述的第一次去噪,主要目的在于去除噪声过多区域,由于眼底血管某 些区域可能含有更好的rPPG信号,而某些区域则含有更多的噪声,因此我们需 要确定哪些区域可能包含较多的噪声,并且在任何处理之前将其移除,以便它们不会影响信号估计。具体操作方法为:对原始rPPG信号,将信噪比低于预设信 噪比阈值的区域或者信号最大幅值高于预设幅度阈值的区域视为噪声过多区域 并予以丢弃,实现所述第一次去噪。在一种优选的实施例中,预设信噪比阈值为 θSNR=0.2,预设幅度阈值为θamp=16。θamp=16是因为采集的眼底视频中,rPPG信号的平均幅度约为4,而这个值为rPPG信号平均幅度的4倍,可以最优过滤 掉包含噪声较多的区域。为了测量SNR,我们用频谱中最大峰值周围[-15bpm, 15bpm]区域的功率谱曲线下面积除以[30bpm,300bpm]频谱中剩余范围的曲线下 面积,所得到的比值即为SNR。在每个时间窗口内,我们可能会移除不同的区 域。为了后续进行时间窗口的融合,可先通过从相邻区域内插来合成缺失区域中 的信号。
步骤S3、对第一次去噪后的rPPG信号进行第二次去噪:使用稳健主 成分分析法对去噪问题进行建模,利用矩阵因子分解去除外部噪声。对于每个时 间窗(长度为10s),我们将48个区域的rPPG信号堆叠到一个大小为10×48 的rPPG矩阵P中。rPPG矩阵中包含原始的rPPG信号,其受到由于不准确的运 动对准、突然的照明变化以及rPPG信号强度跨区域的变化而产生的大量噪声的 污染。然而,所有的区域都应该表达由心脏周期产生的相同的周期性生理信号。此外,当噪声被去除时,在时间窗持续时间上的底层心跳信号的周期性产生了低 秩矩阵,因此,将rPPG矩阵P建模为:
P=Y+N=Y+E+S (4)
其中,Y为包含心跳信号的低秩矩阵;N为噪声矩阵,包含内部噪声E和稀疏外 部噪声S;
具体而言,外部噪声是由突然的照明变化和区域跟踪误差引起的。这 些通常发生在相对于时间处理窗口较短的持续时间内,并且影响少量区域。内部 噪声表征了心跳信号不明显的眼底血管区域。在这种情况下,我们希望在估计心 跳信号时抑制这些区域。为了从P中提取Y的估计值并抑制外部噪声,我们使 用了稳健主成分分析法(robust principalcomponent analysis,RPCA),并建立了 如下优化问题:
其中,||Z||*=∑kσk(Y)表示矩阵Z的核范数,其等于矩阵Z的奇异值σk之和; 矩阵S的l1范数||S||1=∑t,j|S(t,j)|,其等于矩阵S各项元素S(t,j)的绝对值之和; 参数γ控制被划分到噪声分量中的信号能量的相对比例,取值0.02~0.1。较小的 γ值允许更多的信号被视为噪声,本例中优选γ为0.05。
为了求解式(5),可以采用Mansour等人于2016年提出的方法(Mansour等人.Factorized robust matrix completion.Handbook of Robust Low-Rank and SparseMatrix Decomposition:Applications in Image and Video Processing,2016.),即基于因子分解的鲁棒矩阵补全方法。我们将低秩矩阵Z分解成两个因子,令 Z=LMT,其中L和M为两个低秩因子,L∈RW×r,M∈RK×r,W表示时间窗的 长度,K表示每一帧图像被划分成K个区域,r表示秩估计参数且r<W,N;T表 示矩阵转置。这样一来,||Z||*可被表示成如下形式:
其中||·||F表示对其内部矩阵求弗罗贝尼乌斯范数;由于RPCA模型能够消除rPPG测量中的稀疏外部噪声,因此将式(6)代入式(5)进行求解,即可实现从rPPG 矩阵P中消除稀疏外部噪声S。图2中显示了使用RPCA方法对rPPG信号进行去噪的示例。
然而,在某些情况下,可能会发生来自眼底血管区域的信号在整个时 间窗口内都有噪声的情况。这样的噪声分布仍可以被建模为低秩,因此不会被 RPCA去除,鉴于此,我们将在步骤S4中解决这些噪声伪影。
步骤S4、对步骤S3去除外部噪声后的rPPG信号,通过估计心跳信号 的稀疏频谱来去除内部噪声。在一个短的时间窗内,心跳信号是近似周期性的, 其由主频和谐波两部分组成。因此,心跳信号的频谱应该是稀疏的。此外,相同 的心跳信号驱动所有眼底血管区域的rPPG信号的周期性行为。所以,来自所有区域j的信号yj的无噪声频谱xj应该具有相同的支持。
再次考虑式(4)的模型,将消除稀疏外部噪声之后的输出信号zj建模为 zj=F-1xj+ej,并改写成如下的矩阵形式:
其中,ej为内部噪声信号,I为单位矩阵;
如果一个区域有噪声,我们希望在整个时间窗中,该区域都被包含到 矩阵E里,这可以通过强制E中该区域所对应的完整列为零或非零来实现。另 一方面,由于X中的频率分量应为稀疏的,并且在所有区域具有相同的支持, 我们希望X中的列是联合稀疏的,也就是说,我们希望X的整行要么全为零, 要么非零。因此,定义如下优化问题:
其中,λ、μ均为正则项系数,λ的取值范围为0.02~0.3,μ的取值范 围为0.15~1.0;优选地,λ=0.2、μ=1。矩阵X的l2,1范数定义为:
上述优化问题的解决方案可以使用迭代阈值收缩算法例如FISTA算法, 其中收缩函数应适当地应用于X的行范数和E的列范数。图2中显示了恢复的 信号X和噪声E的频谱。使用迭代阈值收缩算法求解式(8),可以计算出矩阵X 和矩阵E,然后去除内部噪声E。FISTA是一种迭代求解如式(8)这种优化问 题的方法,当式(8)即目标函数确定后,调用已写好的FISTA包就可对其进行 求解,具体求解步骤不再详述。
步骤S5、对步骤S4去除内部噪声后的rPPG信号融合时间窗,得到随 时间变化的rPPG信号。通过步骤S4,我们得到了每个时间窗对应的稀疏频谱, 在步骤S5,我们需要进行时间窗的融合,以获得所有时间窗即采集视频的完整 rPPG信号的稀疏频谱。
由于心跳信号随时间变化缓慢,我们可以认为rPPG观测是来自几乎平 稳过程的多通道测量。所以,我们使用滑动窗口其中Qm表示前一个 窗口中不存在的新的rPPG数据,Qo表示在前一个窗口也在当前窗口中的rPPG 数据的比例即当前窗口中与前一个窗口重叠的rPPG信号的比例。为了更好地抑 制噪声,我们构建了加权平均时间融合窗口其中/>是前一窗口的滤波输出,/>是/>存在于当前窗口的部分,α为融合时间窗的阈值, 取值范围为0.01~0.05,优选地α=0.03。然后我们再次使用RPCA过程进一步去 除时间融合窗口/>中的噪声,并估计出所有时间窗口的稀疏频谱,该稀疏频谱 揭示了经过去噪后的rPPG信号的频率、幅度等特征,由这些特征结合时间窗信 息,利用式/>对每个时间窗的稀疏频谱X进行时间窗融合, 即可得到与采集的眼底视频时间等长的随时间变化的rPPG信号。
以上内容是结合具体的优选实施方式对本发明所作的进一步详细说明, 不能认定本发明的具体实施只局限于这些说明。对于本发明所属技术领域的技术 人员来说,在不脱离本发明构思的前提下,还可以做出若干等同替代或明显变型, 而且性能或用途相同,都应当视为属于本发明的保护范围。
Claims (6)
1.一种获取眼底光电体积描记信号的方法,其特征在于,包括如下步骤:
S1、从采集的眼底视频中提取原始rPPG信号,并确定rPPG信号模型;
步骤S1具体包括:
S11、对采集的眼底视频的每一帧图像进行区域划分,并对感兴趣区域进行跟踪;在每一帧图像中,对每个区域的像素强度进行空间平均,以获得原始rPPG信号;
S12、计算每个区域随时间变化的平均像素强度,并在原始rPPG信号中减去该平均像素强度,然后使用带通滤波器将信号的频率范围限制在30~300bpm,该频率范围包含了感兴趣的心脏信号;
S13、对每个区域随时间变化的平均像素强度,利用心跳信号、通道响应函数和通道噪声建立关于时间的信号模型,并将该信号模型改写为向量形式,得到所述rPPG信号模型;
其中,S13具体包括:
对第j个区域随时间t变化的平均像素强度pj(t)建模如下:
pj(t)=hj(t)*yj(t)+nj(t) (1)
其中,*是线性卷积算子,yj(t)表示在第j个区域观测到的心跳信号,hj(t)和nj(t)分别表示通道响应函数和通道噪声;
由于已知心跳信号在频域中是稀疏的,将式(1)改写为向量形式,得:
pj=hj*F-1xj+nj (2)
其中,F-1表示傅里叶变换F的逆变换,F是大小为W的一维离散傅里叶变换,W为时间窗长度;xj为心跳信号yj的稀疏频谱,且有xj∈CW,yj∈RW,C表示复数集合,R表示实数集合;
忽略通道响应函数的影响,得到所述rPPG信号模型:
pj=F-1xj+nj (3);
S2、对原始rPPG信号,根据信噪比或者信号幅度预先进行第一次去噪,以将包含较多噪声的区域移除,并且,在每个时间窗移除不同的区域;为了后续进行时间窗口的融合,先通过从相邻区域内插来合成缺失区域中的信号;
S3、对第一次去噪后的rPPG信号进行第二次去噪:使用稳健主成分分析法对去噪问题进行建模,利用矩阵因子分解去除外部噪声;
步骤S3具体包括:
S31、将第一次去噪后的rPPG信号堆叠成矩阵,得到rPPG矩阵P;
S32、将rPPG矩阵P建模为:
P=Y+N=Y+E+S (4)
Y为包含心跳信号的低秩矩阵;N为噪声矩阵,包含内部噪声E和稀疏外部噪声S;
S33、使用稳健主成分分析法建立如下优化问题:
其中,||Z||*=∑kσk(Y)表示矩阵Z的核范数,其等于矩阵Z的奇异值σk之和;矩阵S的l1范数||S||1=∑t,j|S(t,j)|,其等于矩阵S各项元素S(t,j)的绝对值之和;参数γ控制被划分到噪声分量中的信号能量的相对比例,取值0.02~0.1;
S34、将矩阵Z分解为两个因子,令Z=LMT,其中L和M为两个低秩因子,L∈RW×r,M∈RK×r,W表示时间窗的长度,K表示每一帧图像被划分成K个区域,r表示秩估计参数且r<W,K;T表示矩阵转置;之后||Z||*可被表示成如下形式:
其中||·||F表示对其内部矩阵求弗罗贝尼乌斯范数;将式(6)代入式(5)进行求解,实现从rPPG矩阵P中消除稀疏外部噪声S;
S4、对步骤S3去除外部噪声后的rPPG信号,通过估计心跳信号的稀疏频谱来去除内部噪声;
步骤S4具体包括:
S41、将消除稀疏外部噪声之后的输出信号zj建模为zj=F-1xj+ej,对应的矩阵形式为:
其中,ej为内部噪声信号,I为单位矩阵;
S42、定义如下优化问题:
其中,A:=[F-1I],λ、μ均为正则项系数,λ的取值范围为0.02~0.3,μ的取值范围为0.15~1.0;矩阵X的l2,1范数定义为:
S43、使用迭代阈值收缩算法求解式(8),以计算出矩阵X和矩阵E,然后去除内部噪声E;
S5、对步骤S4去除内部噪声后的rPPG信号融合时间窗,得到随时间变化的rPPG信号。
2.如权利要求1所述的获取眼底光电体积描记信号的方法,其特征在于,还包括:利用眼底相机在非散瞳情况下采集眼底视频。
3.如权利要求1所述的获取眼底光电体积描记信号的方法,其特征在于,步骤S11中采用KLT跟踪器和RANSAC算法对每一帧眼底图像进行感兴趣区域跟踪。
4.如权利要求1所述的获取眼底光电体积描记信号的方法,其特征在于,步骤S2具体包括:对原始rPPG信号,将信噪比低于预设信噪比阈值的区域或者信号最大幅值高于预设幅度阈值的区域视为噪声过多区域并予以丢弃,实现所述第一次去噪。
5.如权利要求4所述的获取眼底光电体积描记信号的方法,其特征在于,所述预设信噪比阈值为θSNR=0.2,所述预设幅度阈值为θamp=16。
6.如权利要求1所述的获取眼底光电体积描记信号的方法,其特征在于,步骤S5融合时间窗包括:
S51、构建加权平均时间融合窗口其中,/>为滑动窗口,Qm表示前一个窗口中不存在的新的rPPG信号,Qo表示当前窗口中与前一个窗口重叠的rPPG信号的比例,/>表示前一窗口的滤波输出,/>是/>存在于当前窗口的部分,α为融合时间窗的阈值,取值范围为0.01~0.05;
S52、再次使用稳健主成分分析法去除时间融合窗口中的噪声,并估计出所有时间窗口的稀疏频谱;
S53、用步骤S52估计出的稀疏频谱,结合时间窗信息,获得与采集的眼底视频时间等长的随时间变化的rPPG信号。
Priority Applications (1)
Application Number | Priority Date | Filing Date | Title |
---|---|---|---|
CN202110987651.4A CN113706413B (zh) | 2021-08-26 | 2021-08-26 | 一种获取眼底光电体积描记信号的方法 |
Applications Claiming Priority (1)
Application Number | Priority Date | Filing Date | Title |
---|---|---|---|
CN202110987651.4A CN113706413B (zh) | 2021-08-26 | 2021-08-26 | 一种获取眼底光电体积描记信号的方法 |
Publications (2)
Publication Number | Publication Date |
---|---|
CN113706413A CN113706413A (zh) | 2021-11-26 |
CN113706413B true CN113706413B (zh) | 2023-08-18 |
Family
ID=78655157
Family Applications (1)
Application Number | Title | Priority Date | Filing Date |
---|---|---|---|
CN202110987651.4A Active CN113706413B (zh) | 2021-08-26 | 2021-08-26 | 一种获取眼底光电体积描记信号的方法 |
Country Status (1)
Country | Link |
---|---|
CN (1) | CN113706413B (zh) |
Citations (2)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
CN106846268A (zh) * | 2017-01-04 | 2017-06-13 | 温州大学 | 一种高斯‑脉冲混合图像噪声去除方法 |
CN108652605A (zh) * | 2018-03-27 | 2018-10-16 | 上海交通大学 | 基于单路ppg信号的实时血压监测装置 |
-
2021
- 2021-08-26 CN CN202110987651.4A patent/CN113706413B/zh active Active
Patent Citations (2)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
CN106846268A (zh) * | 2017-01-04 | 2017-06-13 | 温州大学 | 一种高斯‑脉冲混合图像噪声去除方法 |
CN108652605A (zh) * | 2018-03-27 | 2018-10-16 | 上海交通大学 | 基于单路ppg信号的实时血压监测装置 |
Non-Patent Citations (1)
Title |
---|
A dynamic opto-physiological model to effectively interpret retinal microvascular circulation;Harnani Hassan et al;《ResearchGate》;第1-14页 * |
Also Published As
Publication number | Publication date |
---|---|
CN113706413A (zh) | 2021-11-26 |
Similar Documents
Publication | Publication Date | Title |
---|---|---|
Fan et al. | Robust blood pressure estimation using an RGB camera | |
Sun et al. | Motion-compensated noncontact imaging photoplethysmography to monitor cardiorespiratory status during exercise | |
CN110662480A (zh) | 用于测量和处理对象的生理信号的设备、系统和方法 | |
CN105813564A (zh) | 用于确定对象的生命体征的设备和方法 | |
CN106793962A (zh) | 用于使用视频图像来连续估计人体血压的方法和装置 | |
Yu et al. | Video-based heart rate measurement using short-time Fourier transform | |
JP6620999B2 (ja) | 生体情報計測装置、生体情報計測プログラム、及び生体情報計測方法 | |
CN111387959A (zh) | 一种基于ippg的非接触式生理参数检测方法 | |
Allen et al. | Photoplethysmography (PPG): state-of-the-art methods and applications | |
CA3079625A1 (en) | System and method for camera-based stress determination | |
CN114387479A (zh) | 一种基于人脸视频的非接触式心率测量方法及系统 | |
JP2019097757A5 (zh) | ||
Qiao et al. | Revise: Remote vital signs measurement using smartphone camera | |
JP6866399B2 (ja) | カメラ及びパルスオキシメータを使用した高解像度血液灌流撮像 | |
CN113706413B (zh) | 一种获取眼底光电体积描记信号的方法 | |
Medina et al. | Continuous blood pressure estimation in wearable devices using photoplethysmography: A Review | |
Liu et al. | Lightweight and interpretable convolutional neural network for real-time heart rate monitoring using low-cost video camera under realistic conditions | |
Nair et al. | Illumination invariant non-invasive heart rate and blood pressure estimation from facial thermal images using deep learning | |
Panigrahi et al. | Video-based HR measurement using adaptive facial regions with multiple color spaces | |
Alqudah et al. | Multiple time and spectral analysis techniques for comparing the PhotoPlethysmography to PiezoelectricPlethysmography with electrocardiography | |
Balaraman et al. | Recent Innovations and Improvements in Remote Heart Rate and Heart Disease Measuring Methods Using RGB Camera | |
CN114271800A (zh) | 一种办公环境下的非侵扰式连续血压监测方法及应用 | |
Xiao et al. | Combination of denoising algorithms for video-based non-contact heart rate measurement | |
Gambi et al. | Sensitivity of the contactless videoplethysmography-based heart rate detection to different measurement conditions | |
Ellington | Elucidation of speckleplethysmography as an alternative to photoplethysmography for deriving blood pressure |
Legal Events
Date | Code | Title | Description |
---|---|---|---|
PB01 | Publication | ||
PB01 | Publication | ||
SE01 | Entry into force of request for substantive examination | ||
SE01 | Entry into force of request for substantive examination | ||
GR01 | Patent grant | ||
GR01 | Patent grant |