倪宗軍,陳 輝,張 昀,蘇 敏,鄭秀娟
1.四川大學(xué) 電氣工程學(xué)院,成都610065
2.西安交通大學(xué) 電子與信息學(xué)部,西安710049
成像式光電容積描記法(imaging Photoplethysmography,iPPG)是一種非接觸式生理參數(shù)檢測技術(shù)[1],其基本原理是由于皮下淺層血管血流灌注使得皮膚的顏色發(fā)生輕微變化,因此可以通過成像設(shè)備采集的人體表面部位的視頻信號獲得血容量脈沖(Blood Volume Pulse,BVP)信號[2]。通過測量BVP波形的兩個連續(xù)峰(或谷底)之間的時間間隔,可以從BVP信號中提取心率(Heart Rate,HR)和呼吸率(Respiratory Rate,RR),甚至可以進一步得到人體血氧、血壓、心率變異性、血管微循環(huán)等生理參數(shù),對人體健康的監(jiān)護具有重要意義。但由于人體脈搏波信號十分微弱,因此利用成像式光電手段獲得BVP信號極易受到干擾,噪聲的來源主要包含低頻基線漂移和高頻噪聲,分別對應(yīng)運動噪聲和光照變化產(chǎn)生的噪聲,而人體脈搏波信號的主要信息處于低頻和高頻之間部分。
iPPG技術(shù)可以追溯到2008年,Verkruysse等人首次提出通過iPPG測量生命體征的通用方法。首先,確定暴露皮膚區(qū)域為感興趣區(qū)域(Region of Interest,ROI),通常可以選取為面部的前額區(qū)域,接著取ROI中所有像素值的平均值作為原始BVP信號,而后對原始BVP信號進行帶通濾波,保留與生理信號相對應(yīng)的頻率分量,然后采用快速傅里葉變換(FFT)的頻譜分析方法提取心率和呼吸率對應(yīng)的頻率分量[3]。該通用方法對于運動及光照變化的穩(wěn)健性較差。2010年,Poh等人在該通用方法基礎(chǔ)上引入獨立分量分析(ICA)方法處理面部視頻的RGB三通道信號,從而降低運動噪聲帶來的影響,得出視頻中最強的BVP分量并用于估計HR[4]。接下來,Poh等人又將時域濾波器集成到ICA方法中,從而穩(wěn)健地提取BVP信號,并實現(xiàn)心率、呼吸率及心率變異性等相關(guān)生理參數(shù)的估計[5]。Lewandowska等人提出僅使用兩個通道和應(yīng)用主成分分析來計算心率,該方法可以較為有效地提取符合心肺頻率的信號從而降低光照等噪聲的影響[6]。在標準膚色假設(shè)下,CHROM方法被用來消除運動偽影帶來的影響,該方法在一定程度上提高了iPPG方法的運動魯棒性[7]。雷恒波等人利用經(jīng)驗?zāi)B(tài)分解法(Empirical Mode Decomposition,EMD)將BVP信號分解成一系列不同頻率成分的本征模態(tài)函數(shù)(Intrinsic Mode Function,IMF),最后通過頻域分析法得到心率[8]。EMD去噪的基本原理是從分解的IMF分量中去除高頻噪聲引起的IMF分量。該分量基本不含有效信號成分,去除該分量后,用剩余的IMF分量重構(gòu)原信號,即達到了去噪的效果。但是,該方法只是去除了部分噪聲,其他IMF分量中仍殘留部分噪聲,而且高頻分量中也會含有有效信息。另外EMD還存在著其他問題,如混疊模態(tài)、端點效應(yīng)和停止條件等等。而另外一種信號分析方法-變分模態(tài)分解(Variational Mode Decomposition,VMD)方法,模態(tài)分量數(shù)目K是可以根據(jù)模態(tài)分量瞬時頻率的均值來人為設(shè)定最優(yōu)的,在頻域上能自適應(yīng)地分解出各中心頻率對應(yīng)的有效成分,對低頻段的特征提取具有更好的準確性和穩(wěn)定性,它克服了EMD、經(jīng)驗小波變換等模式中存在的模態(tài)混合、不能正確地消除附加噪聲以及確定本征模態(tài)函數(shù)個數(shù)的問題,具有很好的噪聲魯棒性和降噪效果。但該方法仍然存在各IMF中含有殘留噪聲的問題,如何去除各IMF中存在的噪聲,又保留有效信號是精確提取BVP信號的關(guān)鍵,而且上述方法獲取的生理參數(shù)估計的準確性都受到重構(gòu)的BVP信號質(zhì)量的影響。
因此,本文基于變分模態(tài)分解法提出了一種新型的聯(lián)合能量熵和信號質(zhì)量檢測的自適應(yīng)生理參數(shù)檢測方法,可以很好地分解出BVP信號的高頻噪聲和低頻基線漂移噪聲,實現(xiàn)在自然條件下的人臉視頻中穩(wěn)健地估計心率及呼吸率。
本文提出了一種基于VMD的新型生理參數(shù)測量的方法,共由三個主要步驟構(gòu)成,分別為原始BVP信號獲取、基于變分模態(tài)分解的能量熵閾值自適應(yīng)去噪以及基于信號質(zhì)量檢測的自適應(yīng)生理參數(shù)估計。
獲取原始BVP信號,首先需要選擇受試者臉上暴露的皮膚區(qū)域,即面部ROI。首先使用MTCNN[9]方法檢測人臉,該人臉檢測器采用CNN方法,利用非最大抑制[10]方法尋找最佳的人臉區(qū)域,與目前流行的其他人臉檢測方法相比,具有更高的準確性和魯棒性。在檢測到人臉后,使用CE-CLM[11]方法來形成68個人臉特征點,接著通過選取合適的人臉特征點來包含鼻子和額頭(剔除嘴角和眼睛區(qū)域)的區(qū)域作為面部ROI。
為了有效地提高運動魯棒性,降低運動帶來的噪聲影響,在獲得的面部ROI基礎(chǔ)上,計算視頻RGB三通道的平均值,并采用CHROM[7]方法,由RGB三通道的信號組合得到原始的BVP信號。
1.2.1 變分模態(tài)分解的基本原理
Dragomiretskiy等提出了一種可變尺度的信號分解方法[12],該方法將信號f(t)分解為K個本征模態(tài)分量uk(t),定義本征模態(tài)函數(shù)uk(t)為一個調(diào)幅-調(diào)頻信號,如式(1)所示:
(1)通過Hilbert變換,構(gòu)造每個uk(t)的解析函數(shù)以獲得相應(yīng)的邊際譜。
(2)通過指數(shù)混合調(diào)制到以wk為中心頻率的頻帶上,將每個模型uk(t)的頻譜轉(zhuǎn)移到“基帶”。
(3)由解調(diào)信號的梯度平方和(L2范數(shù))最小來估計各IMF分量的帶寬。即:
(4)通過引入懲罰因子α和Lagrange算子λ(t),將有約束的變分問題轉(zhuǎn)化為無約束的變分問題:
(5)采用乘法算子交替方向法(ADMM)解決上述變分問題,迭代優(yōu)化uk+1、ωk+1、λk+1求得擴展拉格朗日表達式的“鞍點”。采用傅里葉變換,求取二次優(yōu)化問題的頻域解更新公式為:
1.2.2 能量熵閾值自適應(yīng)去噪
VMD分解算法將信號分解為具有一定帶寬的、頻率由高到低的K個模態(tài)分量。由于BVP信號中光照隨機噪聲主要分布在高頻部分,運動偽影噪聲則主要分布在低頻部分,而心肺信號主要信息處于低頻和高頻之間部分,因此需要將K個模態(tài)分量中的低頻部分運動偽影的噪聲去除,高頻部分的光照噪聲去除,將剩余的IMF分量重構(gòu),即可得到去噪后的信號。
首先,依據(jù)各IMF總體信息與BVP信號的相關(guān)性剔除相關(guān)性差的分量,接著根據(jù)剩余IMF中有效信息和噪聲對BVP信號的相關(guān)性不同的特性,對各IMF分量中的有效信息進行檢測與定位,將有效信息引起的系數(shù)剔除,噪聲引起的系數(shù)保留,得到新的分量系數(shù);再將保留的噪聲IMF分量系數(shù)分成若干子區(qū)間,分別計算各子區(qū)間的噪聲能量熵,能量熵最大的區(qū)域,可以認為是全部由噪聲引起的,只要計算出該區(qū)域系數(shù)的能量,再利用噪聲的能量構(gòu)造閾值,對各原IMF分量系數(shù)進行閾值處理,剔除每個IMF分量中噪聲的影響。下面給出改進算法的關(guān)鍵計算步驟。
先計算各IMF分量與原始信號的相關(guān)系數(shù),定義rk(n)為第k個IMF分量與f(n)的相關(guān)系數(shù):
接著將剩余的各IMF分量uk(n)分成t等份,定義rk(i)為第k個IMF分量uk(t)與f(t)的相關(guān)系數(shù):
為使相關(guān)系數(shù)與IMF分量系數(shù)具有可比性,定義規(guī)范化相關(guān)系數(shù)Rk(n)為:
其中,uk(i)為第k個IMF分量中的第i段信號,通過各IMF分量的規(guī)范化相關(guān)系數(shù)與各IMF分量中的每段成分的相關(guān)系數(shù)的比較,得到各IMF分量中噪聲信號的系數(shù)位置。
將IMF中的各個噪聲區(qū)間提取出來,根據(jù)信號在時頻域具有能量守恒特性,設(shè)每個噪聲區(qū)間長度為M,定義第k個IMF分量第p個噪聲區(qū)間的能量為:
將uk(p)分成l等份,每個小區(qū)間的能量記為Eki(p),采樣點數(shù)為則有:
第k個IMF分量中第p個噪聲區(qū)間的第i個子區(qū)間的能量Eki(p)在該分量的總能量Ek(p)中的概率為:
則第i個子區(qū)間對應(yīng)的能量熵為:
得到第k個IMF分量中第p個噪聲區(qū)間的能量熵序列:
搜索第k個IMF分量的熵值最大子區(qū)間,為:
求第k個IMF分量第p個噪聲區(qū)間的第m個子區(qū)間的IMF分量系數(shù)的平均值:
定義第k個IMF分量第p個噪聲區(qū)間的閾值Tk(p)為:
其中,σk(p)為第k個IMF分量的第p個噪聲區(qū)間的熵值最大子區(qū)間的系數(shù)的平均值作為噪聲方差。根據(jù)得到的閾值Tk(p),利用軟閾值函數(shù)分別對各IMF分量中的每個噪聲區(qū)間的分量系數(shù)進行閾值處理,得到新的IMF分量系數(shù),接著利用新的IMF分量系數(shù)重構(gòu)信號,得到去噪后的信號。如圖1所示,將原始的BVP信號經(jīng)過變分模態(tài)分解后,得到5個IMF分量,然后通過相關(guān)性檢測,剔除第一個明顯屬于運動偽影噪聲的IMF分量,起到了去趨勢化的效果,接著通過對每個保存下來的IMF分量進行更細致的相關(guān)性分析,對屬于噪聲區(qū)間的分量進行噪聲能量熵閾值分析,將得到的能量熵閾值重構(gòu)噪聲區(qū)間的分量,從而得到更準確的IMF分量,最后對重構(gòu)的BVP信號進行相應(yīng)范圍的濾波即可得到理想的心率和呼吸率的時域信號。故該噪聲能量熵閾值的確定方法可根據(jù)信號中噪聲的能量特征自適應(yīng)地去除IMF分量的噪聲,重構(gòu)出理想的BVP信號。
1.3.1 基于方差表征序列的信號質(zhì)量檢測
從原始BVP信號中提取得到去噪后的心率信息后,為了繼續(xù)檢測該信號的質(zhì)量,采用Pang等人提出的方差特征序列(Variance Characterization Series,VCS)檢測方法[13],進行心率信號的質(zhì)量檢測。步驟如下:
步驟1找出心率信號中的所有局部最大值Mi,i=1,2,…,以及局部最小值mi,i=1,2,…。
步驟2計算每一個最大局部最大值Mi和接下來七個局部最大值的方差σMi,然后再計算每一個局部最小值mi和接下來七個局部最小值的方差σmi。
步驟3對于每個方差,設(shè):
對于每個最大值和最小值,可以得到時間序列x(t)的兩個VCS。如圖2所示,當不同心率信號的波形變化時,VCS的變化情況。綠色三角代表心率信號的極大值,紅色圓圈代表心率信號的極小值。對于每個最大值的VCS,基于兩個準則來判斷心率信號的質(zhì)量:
(1)δMi的值遠大于THAM。
(2)(i+4)th的最大值和ith的最大值之間的距離遠小于distFAR或者(i+2)th的最大值和ith的最大值之間的距離遠大于distNEAR。
THAM的設(shè)定由經(jīng)驗給定,而distFAR和distNEAR的設(shè)定由之前估計的心率Epre和采樣頻率p給定,即:
圖1 提取心率和呼吸率成分的時間序列步驟
圖2 不同心率成分信號以及相對應(yīng)的方差表征序列
1.3.2 自適應(yīng)生理參數(shù)估計
VCS質(zhì)量檢測合格表明已經(jīng)得到較理想的BVP信號,該信號中存在強烈的心肺頻率成分,接著分別對符合心率和呼吸頻率的心肺信號進行三階的巴特沃斯帶通濾波提取相應(yīng)的頻率成分,濾波范圍分別為(0.7~4 Hz)和(0.1~0.7 Hz),然后使用快速傅里葉(FFT)轉(zhuǎn)換到頻域上分析,最高的頻率成分即是對應(yīng)的心肺頻率fhr、fbr。最終心率(HR)和呼吸值(RR)可估計為:
當VCS質(zhì)量檢測不合格時,表明心肺信號中還存在較強烈的噪聲頻率,這個時候選擇時域的分析方法來獲取心率和呼吸率值,在心肺信號上使用數(shù)峰值的方法來近似得到心肺速率的估計值etime。為了更進一步提高準確性,如圖3所示,接著進行了頻域跟蹤的方法。接下來是具體的步驟:
(1)將心肺信號進行三次樣條插值處理。由于動作或光照帶來的影響是短暫的,進行插值處理后,可以加強心肺頻率成分,從而削弱噪聲帶來的影響。
(2)使用FFT得到插值后的心肺信號的頻譜分布,找到前3個最大頻率成分。
(3)將這3個估計的頻率fest與前一次估計的fpre進行比較,選擇與前一估計最接近的頻率為頻域跟蹤結(jié)果的估計epre。然后與時域分析的估計值etime比較,如果兩種方法的估計值相差在10 beat/min(每分鐘心跳次數(shù))以內(nèi),則取兩者的平均值作為最終的心肺頻率估計efinal。如果相差大于10 beat/min,則說明這段心肺信號的質(zhì)量太差,還存在較多噪聲,則選用機器學(xué)習(xí)的方法來預(yù)測結(jié)果eML。本文選用了背景ROI和人臉ROI的80個特征,其中包括頻域中前5個最大頻率成分、時域信號中波峰到波谷的距離、相鄰波峰、波谷的間距、ICA變換后的結(jié)果等等,最終將80個特征用支持向量機做回歸預(yù)測作為最終的估計efinal。
圖3 頻譜跟蹤算法
選用公開的DEAP數(shù)據(jù)集[14]進行該方法性能驗證并與目前已有方法進行比較。DEAP數(shù)據(jù)集是一個公共多模態(tài)數(shù)據(jù)庫,包括22名參與者(22~35歲)的面部視頻。各參與者均拍攝了40段長度為1 min的面部視頻,同時還通過指夾式設(shè)備和呼吸帶記錄了心率及呼吸率數(shù)據(jù),作為方法性能評估的金標準。視頻拍攝的分辨率為720×576,幀率為50 f/s。參與者被要求靜坐電腦前觀看“情緒誘導(dǎo)實驗”中的各種視頻,隨著視頻的變化會改變參與者臉部的光照,模擬了自然條件下的光照變化,因此DEAP數(shù)據(jù)集上的實驗結(jié)果可以反映各方法在光照變化時的性能效果。實驗從22名參與者中選取金標準數(shù)據(jù)質(zhì)量高(沒有過多的動作偽跡干擾,以保證金標準的準確性)的數(shù)據(jù)進行后續(xù)的對比實驗。
為了進一步進行實際場景中的頭部運動和光照變化的影響實驗比較,本文還采集了各種不同實際場景下的實驗數(shù)據(jù)進行方法的性能比較。使用羅技C920網(wǎng)絡(luò)攝像頭默認設(shè)置拍攝受試面部視頻,使用消費級魚躍血氧儀(YX303)心率結(jié)果作為金標準(由于血氧儀只含心率結(jié)果,故沒有進行呼吸率結(jié)果的比較)。采集了總共20個受試者在5種不同場景下的各5 min視頻。場景1,受試者在室內(nèi)電腦前坐立,距離攝像頭前0.5 m處,可進行點頭、搖頭、微笑等動作,光源為普通LED燈,面部位置光強為300 lux;場景2,受試者在室內(nèi)坐立,距離攝像頭前0.5 m,可進行點頭、搖頭、微笑等動作,光源為普通LED燈,調(diào)節(jié)面部處光強從10 lux到750 lux變化;場景3,受試者在室內(nèi)站立,距離攝像頭前0.5 m處,進行身體左右搖晃、慢速走動等動作,光源為普通LED燈,面部位置處光強為300 lux;場景4,受試者站立在戶外街道,距離攝像頭前0.5 m處,進行點頭、搖頭、微笑等動作;場景5,受試者站立在戶外廣告牌下,距離攝像頭前0.5 m處,進行點頭、搖頭、微笑等動作。
本文中實現(xiàn)了基于面部視頻的實時生理信息采集系統(tǒng)軟件,所有實驗在個人筆記本電腦(微星GP63,Win10操作系統(tǒng),內(nèi)存8 GB,處理器i7-8750H,主頻2.2 GHz)上運行,本文方法運行速度可以平均穩(wěn)定在25 f/s。
BVP信號處理中最常用的方法有離散時間傅里葉變換(DTFT)、離散小波變換(DWT)、峰值計數(shù)(峰值計數(shù))、獨立分量分析(ICA)、自回歸分析(AR)和經(jīng)驗?zāi)B(tài)分解(EMD)。快速傅里葉變換(FFT)復(fù)雜性較低,獨立分量分析(ICA)方法的使用頻率很高,自回歸分析(AR)方法比較新穎,經(jīng)驗?zāi)B(tài)分解(EMD)與變分模態(tài)分解原理相似,所以本文一共選擇了FFT[2]、ICA[4]、AR[15]和EMD[8]方法來進行性能比較。本文在整體方法流程上做了很多改進,因此為了說明其他非BVP因素(如ROI獲取、人臉跟蹤方法、信號平滑以及頻率跟蹤算法)的改進帶來的提升,將整個方法與其他論文提出的整體方法進行比較,本文選取的ROI區(qū)域如圖4所示,其中綠色空心圓表示選取的部分特征點,藍色區(qū)域是選取的粗略ROI,紅色矩形是精細ROI,包括鼻子和前額區(qū)域。
圖4 面部ROI選擇結(jié)果示例
為了公平比較,ICA、FFT、AR和EMD中用到的參數(shù)均參照原論文設(shè)置為最優(yōu)值,在處理不同視頻時,所有參數(shù)也保持不變。視頻時間窗都選擇30 s(包括25 s的重疊窗),最終比較結(jié)果如表1~3所示。
由表1和表2結(jié)果可見,ICA方法的性能優(yōu)于FFT方法,這是因為ICA方法增加了多重時域濾波,使得時間序列更加穩(wěn)定。而傅里葉變換是一種經(jīng)典的頻譜分析方法,在采樣頻率穩(wěn)定,且數(shù)據(jù)長度足夠的情況下,生理參數(shù)估計可以取得較好的精確度。然而在使用傅里葉變換分析低頻段頻譜時,其頻率分辨率是不足的,而且實際應(yīng)用中通常含有光照、表情動作等因素的影響,會在獲取的原始BVP信號中引入了非高斯白噪聲,將帶來多個偽峰值,尤其是運動噪聲帶來的頻率峰值。通常運動干擾信號的頻率成分會超過BVP信號中心肺頻率成分。在這種條件下,單純地將最高峰值的頻率點作為心肺頻率,其估計精度將很難保證。AR[15]方法雖然可以減少人工光照的影響,但是對于一些運動偽影的干擾則不能很好地抑制。而ICA方法是假定觀測到的信號為盲源的線性混合,而實際的環(huán)境光噪聲和其他運動噪聲可能是非線性和時變的,特別是在運動過程中,因為運動引起的血容量的生理變化也可能是非線性的。此外,ICA假設(shè)源信號的數(shù)量等于RGB通道的數(shù)量,但潛在噪聲源的數(shù)量是不同的,且在實際應(yīng)用難以確定。因此,ICA并不總是能正確地分解RGB通道來提取所需的心肺信號的頻率。而EMD方法是通過去除屬于噪聲的IMF分量進行降噪,但存在著模態(tài)混疊和分解數(shù)目不可控的缺陷,如何選擇去除合適的IMF分量是重構(gòu)出精確的BVP信號的關(guān)鍵,而且之前討論過,去除的IMF分量里也包含著有效的信息,而剩余的IMF分量中還存在著其他噪聲,這是EMD方法準確度不是很理想的關(guān)鍵原因。從表3可知,當在更復(fù)雜的現(xiàn)實場景下(包含變化的光照、頭部運動),這些方法達不到滿意的估計結(jié)果。主要原因是,已有方法在人臉ROI獲取時對動作的跟隨性能不高,隨著人臉轉(zhuǎn)動,會出現(xiàn)人臉跟蹤失敗,進而無法得到理想的BVP信號。而且,在已有方法中也缺乏對復(fù)雜噪聲的處理,如圖5所示,原始BVP信號中包含較強烈的偽影噪聲時,可能會出現(xiàn)運動偽影帶來的峰值高于心率成分的峰值,其中因AR方法與FFT方法一樣使用綠色通道,故未標出。ICA方法在一定程度上減弱了運動偽影帶來的影響,而EMD方法通過選取合適的本征模態(tài)函數(shù)也可以較大程度減少運動偽影的噪聲,但都還殘留著部分噪聲。本文的方法簡化了環(huán)境光的過濾,著重在ROI的精確選取上,以及應(yīng)用CHROM方法可以對運動噪聲起到一定的抑制作用。如圖6所示,當同時出現(xiàn)強烈的運動偽影及光照變化的噪聲時,其他方法的頻譜中會殘留大量的偽峰,而通過應(yīng)用基于VMD方法的能量熵閾值自適應(yīng)去噪方法,具有更好的抗噪性,能夠重構(gòu)出更精確的BVP信號,而且根據(jù)重構(gòu)出的BVP信號的質(zhì)量通過VCS信號質(zhì)量分析可以自適應(yīng)地選擇不同的分析方法以及提出使用的頻譜跟蹤算法中機器學(xué)習(xí)的方法利用了背景ROI和面部ROI的信號得到光照和運動噪聲特征,可以更準確地估計出心肺信號頻率。然而本文方法還不能有效地去除因劇烈運動帶來的噪聲影響。因此,在未來的工作中,將提高本文方法的健壯性,優(yōu)化方法流程,使其適應(yīng)更多的實際場景,從而提高方法在不同情況下的準確性。
表1 DEAP數(shù)據(jù)集下不同方法的心率比較結(jié)果
表2 DEAP數(shù)據(jù)集下不同方法的呼吸率比較結(jié)果
表3 自采數(shù)據(jù)集下不同方法的心率比較結(jié)果
本文提出了一種基于人臉視頻測量生理信息的方法,通過改進ROI檢測方法,可以自適應(yīng)地精準選擇ROI。針對直接去除低頻或高頻分量的EMD去噪方法中,不能有效區(qū)分各IMF分量中的有效信息和噪聲信息的問題,本文提出了基于變分模態(tài)分解的相關(guān)能量熵閾值自適應(yīng)去噪方法,該方法根據(jù)噪聲的能量熵自適應(yīng)地確定各IMF分量的閾值,重構(gòu)出更精確的BVP信號。相比傳統(tǒng)單一分析方法如FFT、數(shù)峰值等等,本文還使用了基于信號質(zhì)量的自適應(yīng)分析方法,有效地提升了應(yīng)用范圍,最后為了提高測量的精度。同時,通過采用頻譜跟蹤算法,可以更準確地估算心肺頻率值。在公開的數(shù)據(jù)集和自采數(shù)據(jù)集上進行的性能驗證實驗,證明了本文方法在生理參數(shù)估計精度較已有方法有較大提升,并且可以更好地適用于實際環(huán)境中。
圖5 不同算法對含較強運動偽影噪聲的BVP信號處理結(jié)果示例
圖6 不同算法對含強烈的運動偽影噪聲和光照噪聲的BVP信號處理結(jié)果示例