基于歐拉放大的橈動脈脈搏波提取方法

李江山, 楊學志, 吳 秀, 金 兢

(1.合肥工業(yè)大學 計算機與信息學院,安徽 合肥 230601; 2.工業(yè)安全與應急技術(shù)安徽省重點實驗室,安徽 合肥 230601)

0 引 言

脈搏波是反映人體心血管系統(tǒng)健康狀況的重要生理參數(shù)。心臟泵血運動促進血液流經(jīng)全身,這是整個動脈系統(tǒng)循環(huán)不斷的動力源泉。動脈系統(tǒng)隨心臟運動而周期性地運動,并引起皮膚表面動脈組織的變化,從而形成搏動特征,脈搏波就是記錄這一搏動特征的參數(shù)。心血管功能的變化或者病變必然會導致脈搏博動特征的改變,因此,從脈搏波波形中可以獲取心血管系統(tǒng)的特征信息[1]。

血液傳播引起的皮膚表面變化可以分為2類:① 皮膚表面顏色明暗的變化;② 由血液隨著心臟泵血造成皮膚表面運動的變化[2]。這些變化難以被肉眼察覺,而通過視頻放大算法可以將之可視化并從中提取相關(guān)生理信息,因此,基于視頻法的非接觸式脈搏波檢測算法的核心工作就是從含有人體皮膚表面組織(面部、腕部等)的視頻中提取出脈搏波。

文獻[3]最早提出了通過普通攝像頭拍攝人臉視頻來檢測脈搏波的方法,該方法對一段10 s左右含有人臉信息的視頻進行檢測。首先對每幀視頻進行人臉檢測選出含有人臉信息的區(qū)域作為感興趣區(qū)域(region of interest, ROI);然后計算ROI內(nèi)部的像素平均值,得到RGB彩色空間下的3個顏色通道;再利用獨立成分分析(independent component analysis,ICA)算法將3個顏色通道作為3個獨立的源信號,并選擇其中的G通道作為待處理信號;最后通過時域分析和功率譜估計得到脈搏波和心率值。但是在后來的研究中發(fā)現(xiàn)該方法具有一定的缺陷,血容量脈沖(blood volume pulse,BVP)信號不一定是ICA分量的第1個基向量,有時候是其他的基向量[4]。

文獻[5]提出了基于歐拉視角的視頻放大(Eulerian video magnification,EVM)算法,該算法能夠?qū)⑷祟惾庋鄄豢梢姷奈⑿∽兓梢暬?為非接觸式脈搏波提取的研究提供了新的思路。文獻[6]將EVM算法應用于非接觸式脈搏波檢測算法中,提出了對由血液傳播引起的面部微弱明暗變化進行顏色放大,并從中提取出脈搏波、心率等生理信息;該方法通過心電圖儀進行驗證,具有較高的準確性,但需要在受測者極度配合的情況下才能準確地測量心率。文獻[7]將EVM應用于人臉視頻上,對人體面部血液明暗變化進行顏色放大估算心率值；但是該方法在非合作情況下準確性較低。文獻[8]使用特征點跟蹤和旋轉(zhuǎn)校正等方法降低人臉運動產(chǎn)生的運動干擾,在非合作情況下具有一定優(yōu)勢；但該方法使用理想帶通濾波做時域處理,所提取的脈搏波易損失局部細節(jié)。文獻[9]將EVM方法與ICA方法進行對比,驗證了EVM方法比ICA方法具有更好的檢測精度。文獻[10]將EVM算法應用于腕部橈動脈脈搏波的提取檢測中,對腕部橈動脈微小搏動變化進行運動放大,從而得到脈搏波信息;該方法相比于面部顏色放大,具有更多的脈搏博動信息,橈動脈脈搏波保留了更豐富的局部細節(jié)[11];但該方法對環(huán)境中存在的運動干擾十分敏感,需要受測者十分配合,且成功率較低。

針對以上問題,本文在歐拉視角框架下,利用視頻序列的亮度變化與位置變化的對應關(guān)系,建立基于亮度變化的尺動脈跳動的描述,通過對變化的亮度信號進行提取與處理,實現(xiàn)對橈動脈脈搏波的檢測。

1 本文算法

本文提出了一種基于運動放大的脈搏波提取方法,通過對腕部橈動脈進行視頻運動放大提取脈搏波;然后對ROI進行最優(yōu)化處理,選取信噪比(signal to noise ratio,SNR)最高的區(qū)域作為待處理源信號;最后選取ROI中3個子區(qū)域作為觀測信號,使用盲源分離算法提取脈搏波,進一步降低運動干擾。

1.1 視頻運動放大

橈動脈搏動信號十分微弱,人類的肉眼無法直接觀察這一現(xiàn)象[11],直接使用計算機或者普通攝像頭設(shè)備也難以直接獲取視頻中的脈搏信息。為了精準地獲取視頻中的橈動脈信號,首先,對腕部的橈動脈區(qū)域進行視頻放大,這里使用EVM算法對含有橈動脈的視頻從空間和時間2個角度進行處理,結(jié)合橈動脈的搏動特征來提取脈搏波。

本文算法流程如下:首先,在空間域上,利用多尺度的拉普拉斯金字塔將視頻幀序列分解為不同空間頻率的子帶;然后，使用無限脈沖響應(infinite impulse response,IIR)濾波器做時域濾波,對視頻序列使用傅里葉變換算法做時域變換,使用IIR濾波器使得視頻信號保持在一定范圍的帶通內(nèi),用于后續(xù)的最優(yōu)ROI區(qū)域選取。

(1) 空間濾波。圖像處理中的金字塔分解處理是一種多尺度、多分辨率的方法,拉普拉斯金字塔的作用是將圖像在不同的空間頻帶上,針對不同的特征與信息進行處理[12]。拉普拉斯金字塔是基于高斯金字塔基礎(chǔ)上進行運算的,因此首先對原圖像進行高斯金字塔分解。以視頻幀F(xiàn)0作為高斯圖像金字塔第0層,對圖像進行高斯濾波和隔行隔列的降采樣,得到高斯金字塔的第1層金字塔F1;再對第1層執(zhí)行高斯濾波和降采樣,得到第2層F2。重復以上過程,得到F1、F2……FN構(gòu)成高斯圖像金字塔，N為所求解的最大金字塔層次。Fl的表達式為:

(1)

(2)

至此,高斯金字塔構(gòu)建過程完成。

(3)

由此,L1、L2、L3……LN構(gòu)成了拉普拉斯金字塔。金字塔每層圖像分別表示不同尺度下的頻率空間帶。通過將不同層次的灰度層拉開,有助于顯現(xiàn)微弱的細節(jié)變化。

(2) 時域濾波。IIR濾波器在其截止頻率附近的某一頻率范圍內(nèi),存在一個信號衰減程度逐漸增加的過程,即常見的濾波器滾降現(xiàn)象，該特性更符合人體橈動脈的頻域特性。因此,使用IIR濾波器對拉普拉斯金字塔圖像分解后的圖像序列進行時域濾波,通過設(shè)置IIR濾波器的上下帶通,獲取脈搏波所在頻率范圍內(nèi)的橈動脈搏動信號。

(3) 視頻重構(gòu)。對時域濾波后的信號進行放大,將R通道的幅值保持不變,G通道和B通道分別放大100倍,即對其通道內(nèi)的每個像素點值乘以100。每幀放大之后都擴展到原始視頻幀大小,并與原始視頻幀圖像相加,得到放大后的視頻幀,使運動信號的變化情況能被肉眼觀察到。腕部動脈視頻放大前、后效果如圖1所示。

圖1 腕部動脈視頻放大前、后效果對比

1.2 最優(yōu)ROI選取

視頻運動放大之后,對視頻幀進行分塊,將640×480的視頻幀分成40×40的小窗口,該窗口大小能夠保證完整覆蓋橈動脈搏動區(qū)域,視頻將被分為192個窗口。計算每個窗口在時域上的SNR,然后求出的最大信噪比作為最優(yōu)ROI。因為橈動脈的搏動信息比環(huán)境噪聲具有更強烈的運動信號,含有橈動脈的區(qū)域具有更高的SNR,因此選出最高SNR的區(qū)域必定含有橈動脈的搏動信號。執(zhí)行步驟如下所述。

(1) 沿著圖像的x軸與y軸方向,從左到右、從上到下,以40為步長進行分塊,共計192個分塊。分塊在時間維度上保持不變,令Fmn(s,k,t)為分塊像素值的函數(shù),其中,m、n為分塊的位置坐標,s、k、t分別為分塊內(nèi)部像素的位置以及時間序列。

(2) 求第(m,n)個ROI子區(qū)域的像素平均值,計算公式如下:

(4)

其中,ROImn(t)為第(m,n)個ROI子區(qū)域的第t幀像素平均值,顯然ROImn表示第(m,n)個ROI子區(qū)域塊在時間序列上的信號值,是一個1×300的一維信號。

(3) 求每個子ROI的SNR,并求得最大SNR的子區(qū)域,計算公式如下:

SNRmax(u,v)=max{SNR(ROImn)}

(5)

其中,SNG()為求信噪比的函數(shù);max為求最大值的函數(shù)。

ROI的SNR分布如圖2a所示,通過對視頻求出的最佳ROI子區(qū)域如圖2b所示,對照視頻幀選出橈動脈區(qū)域。最優(yōu)ROI的選取能夠從一定程度上降低環(huán)境噪聲的干擾,提高脈搏波的提取精度,所求的ROI區(qū)域可以作為下一步ICA的源信號。

圖2 ROI的SNR分布

1.3 脈搏波提取

選出最優(yōu)ROI之后,采用ICA算法[13]對ROI區(qū)域進行信號分離,進一步分離噪聲混疊信號。手臂皮膚表面存在多個運動源信號,其中一部分是橈動脈的搏動信號,另一部分是環(huán)境中存在的噪聲干擾信號。這些噪聲可以分別視為一種源信號。使用ICA算法能夠?qū)⑦@些信號進一步地分開,從混合信號中分離出脈搏波信號,進一步降低其他運動信號對脈搏波的干擾。

假設(shè)x(t)=[x1(t)x2(t) …xm(t)]T為ROI區(qū)域內(nèi)m個時間長度為n的運動信號,其中t=1,2,…,n;設(shè)s(t)=[s1(t)s2(t) …sm(t)]T為m個未知的相互獨立的源信號,其中t=1,2,…,n。則觀測信號與源信號存在如下關(guān)系:

x=As

(6)

其中,A為未知時不變的n×m混疊矩陣。使用ICA算法從觀測信號矢量中通過矩陣運算得到相互獨立的源信號矢量。當m=n時,可以構(gòu)造混合矩陣A的逆矩陣A-1求解源信號,則有:

y=Wx

(7)

其中:y為系統(tǒng)分離的信號矢量,也是源信號s(t)的近似值;W為解混矩陣,即W=A-1。

基于概率論和統(tǒng)計學原理,如果源信號大多數(shù)服從高斯分布,那么經(jīng)過線性混疊后的信號依然服從高斯分布。本文取觀測信號與源信號的個數(shù)為3,即m=n=3,根據(jù)概率學理論可得解混矩陣W,從而求得源信號的近似值s(t)。

基于以上分析,從ROI區(qū)域中選取3個區(qū)域作為觀測信號,對每個區(qū)域做像素平均得到平均值,該平均值作為該區(qū)域在某一時間點上的觀測值。對所有視頻幀序列上的該區(qū)域做像素平均,得到該區(qū)域在時間序列上的一維觀測信號,該觀測信號是橈動脈信號與噪聲的混合。通過ICA可以從觀測信號中分離出橈動脈運動信號和噪聲,得到較為純凈的橈動脈運動信號。脈搏波信號提取結(jié)果如圖3所示。

圖3 脈搏波信號

在提取脈搏波之后,通過快速傅里葉變換(fast Fourier transform,FFT)分析脈搏波的功率譜,然后從功率譜中估算心率值,功率譜中的最大值對應的頻率即心率值,心率可以作為評價脈搏波提取質(zhì)量的定量指標。FFT在小樣本數(shù)目的前提下具有較低的精度,因此本文采用一種傅里葉系數(shù)迭代插值法對輸入的脈搏波信號進行增量差值,在保證不改變數(shù)據(jù)內(nèi)容的情況下提高FFT運算精度,實現(xiàn)算法精度與復雜度的平衡[14]。

從脈搏波信號估算心率的3個步驟如下所述。

(1) 計算樣本數(shù)N=300的脈搏波信號y的功率譜P(t)和功率譜最大值位置對應的頻率fmax,即

P(t)=|FFT(y)|2,t=0,1,…,N-1;

fmax=max{P(t)}

(8)

(2) 對功率譜P(t)進行傅里葉系數(shù)迭代插值。首先初始化參數(shù),令迭代方程中的傅里葉系數(shù)起始偏差e0=0和迭代次數(shù)Q=50;然后進行1≤k≤Q的for循環(huán)迭代。實驗驗證可得,迭代次數(shù)為50時,傅里葉系數(shù)偏差基本保持收斂[14]。循環(huán)部分的公式如下:

ek=ek-1+r(ek-1)

(9)

其中:d=±0.5;r(ek-1)=0.5(|S+0.5|-|S-0.5|)/(|S+0.5|+|S-0.5|)。

(3) 估算心率。將(8)式求得的最大值對應頻率和(9)式求得的迭代法偏差值代入下式,得到估算的心率值,即

(10)

其中:HRy為脈搏波信號y的心率值;eQ為迭代了Q次之后的傅里葉系數(shù)偏差;fs為視頻采樣率;N為樣本數(shù)。所求心率作為脈搏波信號定量分析的參考指標。

2 實驗結(jié)果與分析

本文設(shè)計了2組實驗,驗證本文方法的有效性和可行性。首先,使用傳統(tǒng)硬件脈搏波檢測設(shè)備作為實驗真值,檢測本文算法的準確性;然后,與文獻[8]和文獻[10]提出的非接觸式的脈搏波檢測方法對比,以心率作為定量分析指標,在不同的光照條件下驗證本文方法具有的優(yōu)勢。

2.1 實驗設(shè)計

實驗裝置設(shè)計如圖4所示。使用一個普通網(wǎng)絡攝像頭連接筆記本電腦,攝像頭對準手腕橈動脈位置,置于手腕上方0.5 m左右距離,同時,受測者需要主動配合,將手腕平放在桌面上,避免較大幅度的手臂運動。

圖4 實驗裝置示意圖

電腦通過MatLab控制攝像頭錄制視頻。設(shè)置攝像設(shè)備采集圖像為RGB彩色空間,幀率為30幀/s,視頻分辨率為640×480,拍攝時長為10 s,總計300幀圖像。10名成年人參加了本次準確性測試,包括8名男性和2名女性,在每次測試時,采集受測者的手腕視頻信息,同時使受測者佩戴基于指尖的PWS-20D脈搏波檢測設(shè)備。該設(shè)備采用光學體積描記術(shù)原理,基于指尖進行脈搏波的采集,可以精確地采集受測者的脈搏波信息,在拍攝視頻時,使用該設(shè)備測量脈搏波作為參考值。同時,從定量分析的角度來評價算法的有效性,采用平均誤差ME、標準差S、均方誤差RMSE、準確率平均值HRac和皮爾森相關(guān)系數(shù)r作為評價指標。

2.2 實驗結(jié)果分析

為了驗證本文方法的有效性和準確性,按照本節(jié)所述方法,對實驗的10名志愿者進行非接觸式手腕脈搏波提取,同時使用PWS-20S脈搏波檢測儀同步測量分析。10名受測者每人分別做3組測試,每組測試包括非接觸式和同步接觸式脈搏波采集各1次,共30對數(shù)據(jù)。

采用Bland-Altman法進行一致性評估,分析結(jié)果如圖5所示。

圖5 Bland-Altman法估算2種測量方法的一致性

圖5中散點K的計算方法如下:

K(x,y)=((HRr+HRv)/2,HRr-HRv)

(11)

其中:HRr為PWS-20D脈搏波采集儀的測量心率;HRv為本文方法的估算心率;x坐標(橫坐標)為一次實驗中2種方法測量值的平均值;y坐標(縱坐標)為同一次實驗中2種方法測量值的差值。置信率取95%(S=1.96)。通過Bland-Altman方法對實驗數(shù)據(jù)進行一致性估計,誤差均值為0.160 0,標準方差為2.451 2,95%的置信區(qū)間為[-4.507 2, 4.624 1],說明本文方法與接觸式脈搏波測量設(shè)備具有很好的一致性。

為了充分驗證本文方法對光照影響的抗干擾能力,實驗將在自然光、日光燈與弱光3種不同的光照強度場景下分別與文獻[8]、文獻[10]方法進行對比。其中,文獻[8]采用了歐拉顏色放大方法,基于面部色度信息提取脈搏波;文獻[10]采用了歐拉運動放大方法,基于橈動脈提取脈搏波。本次實驗的視頻時長為10 s,幀率為30幀/s。

不同光照場景下心率檢測結(jié)果見表1所列。實驗結(jié)果表明,在不同光照場景下,本文方法均具有較好的性能,且在室內(nèi)弱光場景下,本文方法準確率仍能達到90.69%,平均誤差只有2.054 5,具有明顯優(yōu)勢。

實驗結(jié)果的對比如圖6所示。由圖6可知,本文方法在各方面指標均優(yōu)于文獻[8]和文獻[10]的方法。其主要原因如下:

(1) 本文方法在ROI的選取上進行了精確劃分,排除了環(huán)境噪聲的干擾。通過計算每個子ROI區(qū)域的SNR,選出含有橈動脈運動信息量最高的ROI區(qū)域作為源信號;再通過ICA算法將選出的ROI區(qū)域中含有的混疊局部運動噪聲進一步分離,使提取的脈搏波受環(huán)境噪聲的干燥較小。

(2) 文獻[8]方法對面部顏色明暗變化進行放大以提取脈搏波,顏色變化對光照變化比較敏感,因此該方法隨著光照條件的下降檢測精度也隨之降低。

(3) 文獻[10]方法同樣采用了基于橈動脈運動放大的方法,但采用了Lab顏色空間作為皮膚分割和ROI選取的依據(jù),Lab空間的L通道受光強影響較大,環(huán)境光照對放大效果有較大的影響。

表1 不同光照場景下心率檢測結(jié)果

圖6 實驗結(jié)果的對比

3 結(jié) 論

本文提出了一種非接觸式脈搏波的提取方法。該方法基于歐拉運動放大原理,首先對腕部橈動脈的搏動進行視頻放大；再對放大后的視頻進行ROI區(qū)域劃分,對每個子區(qū)域求解SNR,選取橈動脈運動信號占比最高的區(qū)域作為最優(yōu)ROI;然后使用ICA算法對ROI內(nèi)的脈搏波信號進一步做信號分離,分離出脈搏波運動信號和環(huán)境噪聲,使脈搏波信號更加純凈,處理后的脈搏波信號具有更好的局部細節(jié)信息。分別在自然光、日光燈、室內(nèi)弱光場景下,在多組實驗數(shù)據(jù)集上進行了對比實驗,實驗結(jié)果驗證了本文方法的有效性和穩(wěn)定性,與基于色度信息分析的脈搏波提取方法、基于人臉視頻的脈搏波提取方法相比,本文方法更具有優(yōu)勢。