基于歐拉影像放大的非接觸式心率測(cè)量方法

蘇培權(quán)，許 亮，梁永堅(jiān)

(廣東工業(yè)大學(xué) 自動(dòng)化學(xué)院，廣州 510006)

0 引言

心率作為人體重要生命體征，是評(píng)估人體健康重要參數(shù)，其中，靜息心率是鑒別心臟健康狀況及其他臨床診斷風(fēng)險(xiǎn)預(yù)測(cè)重要因子[1]，因此對(duì)靜息心率日常測(cè)量將有助于對(duì)人體身體狀況監(jiān)測(cè)以及相關(guān)疾病預(yù)防，特別是對(duì)老年人身體功能性下降和疾病風(fēng)險(xiǎn)的預(yù)測(cè)有著重要的意義[2]。目前，心電圖法是心率測(cè)量中應(yīng)用最廣、準(zhǔn)確度高的一種方法，但是由于這類接觸式測(cè)量法中傳感器設(shè)備需要與人體皮膚接觸，造成使用上不便。

早期非接觸式心率測(cè)量法有微波多普勒雷達(dá)[3]、超寬帶脈沖雷達(dá)分析[4-5]以及熱成像處理[6]，然而這些方法成本昂貴，且需要專業(yè)操作。另一方面，家庭式醫(yī)療診斷測(cè)量和遠(yuǎn)程健康狀況監(jiān)控越來越被重視，并且處于快速發(fā)展階段[7]，所以尋找一種低成本、操作方便且實(shí)時(shí)性較強(qiáng)非接觸式心率測(cè)量法具有重要意義。文獻(xiàn)[8] 對(duì)臉部感興趣皮膚區(qū)域采集30 s視頻，根據(jù)皮膚亮度變化采用自回歸頻譜實(shí)現(xiàn)脈搏測(cè)量，是早期通過視頻采集實(shí)現(xiàn)脈搏測(cè)量的方法。近年來，隨著計(jì)算機(jī)視覺技術(shù)發(fā)展，有研究人員利用光電容積描記法(Photo Plethysmo Graphy, PPG)原理對(duì)人體皮膚采集到視頻進(jìn)行處理[9-10]，從而實(shí)現(xiàn)對(duì)心率非接觸式測(cè)量，該方法在降低視頻噪聲時(shí)大多采用線性濾波器，無法降低與心率同頻段信號(hào)噪聲，造成心率提取準(zhǔn)確度較低。Poh等[11-12]和Pursche等[13]采用獨(dú)立成分分析(Independent Component Analysis, ICA)算法進(jìn)行盲源分離去噪，對(duì)于測(cè)量對(duì)象人臉中非皮膚區(qū)域，特別是眼睛眨動(dòng)會(huì)對(duì)信號(hào)波形產(chǎn)生較大干擾，且干擾頻率和心率相接近，不易被濾除[14]。另外，盲源分離并沒有提供一個(gè)判斷方法來區(qū)分分離后哪一組獨(dú)立信號(hào)相對(duì)于實(shí)際心率信號(hào)最具有真實(shí)性，從而影響到心率計(jì)算精度[15]。Xu等[16]針對(duì)盲源分離缺點(diǎn)，受光學(xué)血氧儀啟發(fā)，提出了針對(duì)人體皮膚基于脈搏變化交流信號(hào)成分分析的簡(jiǎn)單數(shù)學(xué)模型心率計(jì)算，該方法在測(cè)量過程中要持續(xù)45～90 s，而且整個(gè)過程要求被測(cè)對(duì)象保持不動(dòng)。短時(shí)傅里葉測(cè)量心率[17]，固定變換窗口寬度限制了心率估計(jì)，不正確窗口大小直接導(dǎo)致心率錯(cuò)誤測(cè)量。文獻(xiàn)[15,18]在Poh等[11]提出方法基礎(chǔ)上對(duì)盲源分離算法改進(jìn)，一定程度上提高了算法魯棒性。針對(duì)測(cè)量過程中被測(cè)對(duì)象移動(dòng)和環(huán)境亮度變化對(duì)測(cè)量精度的影響，文獻(xiàn)[19] 通過臉部跟蹤和對(duì)臉部變化亮度矯正有效提高算法對(duì)移動(dòng)和環(huán)境變化的魯棒性。文獻(xiàn)[20]根據(jù)建立皮膚模型和線性盲源分離從臉部多個(gè)局部區(qū)域提取的體積描記器(PlethysmoGraph, PG)信號(hào)，并使用多數(shù)投票法估計(jì)心率，進(jìn)一步提高對(duì)被測(cè)對(duì)象移動(dòng)和環(huán)境亮度變化的魯棒性，一定程度上克服Li等[19]在復(fù)雜測(cè)量背景臉部變化亮度矯正可靠性低的不足。人臉視頻PPG原理是皮膚對(duì)光吸收和反射隨著皮膚下血管血液灌注量變化而變化，人臉視頻通過捕獲這種變化來估計(jì)心率。非接觸式視頻心率估計(jì)方法，幾乎都在臉部皮膚上根據(jù)PPG原理演變而來，光電容積描記法相關(guān)醫(yī)學(xué)研究[21-22]表明，人體周圍環(huán)境溫度和表皮組織溫度直接影響到皮膚血流灌注，對(duì)PPG信號(hào)產(chǎn)生較大干擾，導(dǎo)致血氧飽和度及心率估計(jì)精度下降。人臉直接與周圍環(huán)境接觸，受環(huán)境溫度影響較大，所以造成基于光電容積描記法原理人臉心率測(cè)量精度和穩(wěn)定性下降。因此，上述心率測(cè)量方法存在心率同頻段噪聲干擾大、受環(huán)境溫度影響較大、對(duì)分離信號(hào)真實(shí)性缺乏判斷等問題。

針對(duì)上述問題，并根據(jù)靜息心率可通過橈動(dòng)脈一分鐘跳動(dòng)次數(shù)來估算，本文提出了一種基于歐拉影像放大技術(shù)[23]的心率測(cè)量新方法。該方法運(yùn)用歐拉影像放大技術(shù)實(shí)現(xiàn)橈動(dòng)脈跳動(dòng)微小變化放大，對(duì)放大后的視頻幀作亮度方差統(tǒng)計(jì)，并根據(jù)統(tǒng)計(jì)結(jié)果分割提取橈動(dòng)脈跳動(dòng)區(qū)域視頻圖像，對(duì)所提取橈動(dòng)脈部位亮度通道進(jìn)行時(shí)頻分析，實(shí)現(xiàn)心率非接觸式自動(dòng)測(cè)量。

1 測(cè)量原理

由于心率測(cè)量可通過對(duì)橈動(dòng)脈跳動(dòng)次數(shù)的計(jì)數(shù)間接實(shí)現(xiàn)。而在影像中橈動(dòng)脈很細(xì)微、跳動(dòng)難以用人肉眼觀察出來，所以本文所提方法先采用歐拉影像放大技術(shù)對(duì)影像中橈動(dòng)脈微小跳動(dòng)進(jìn)行放大，使脈搏跳動(dòng)可視化，同時(shí)便于橈動(dòng)脈跳動(dòng)部位定位處理，如圖1所示。在歐拉影像放大技術(shù)中，對(duì)于微小運(yùn)動(dòng)對(duì)象采用空間多尺度方法放大像素亮度變化值，實(shí)現(xiàn)橈動(dòng)脈跳動(dòng)信號(hào)放大。對(duì)經(jīng)過放大處理視頻幀，由于亮度包含了圖像微小運(yùn)動(dòng)信息，圖像變化主要在亮度上表達(dá)出來，因此亮度變化可用來解析橈動(dòng)脈跳動(dòng)信息。

圖1 心率測(cè)量原理

Fig. 1 Schematic diagram of heart rate measurement

橈動(dòng)脈微小跳動(dòng)信息主要集中在視頻幀像素點(diǎn)亮度通道的變化，所以本文對(duì)采集橈動(dòng)脈區(qū)域視頻進(jìn)行歐拉影像放大處理。歐拉影像放大處理中建立視頻幀拉普拉斯金字塔對(duì)視頻幀空間多分辨率分解，對(duì)視頻幀空間分辨率一樣的圖像層進(jìn)行時(shí)域處理。時(shí)域處理算法，以及與運(yùn)動(dòng)放大之間關(guān)系可根據(jù)光流分析中基于一階泰勒級(jí)數(shù)展開微分近似法進(jìn)行分析，一維運(yùn)動(dòng)信號(hào)具體分析說明如下所示。

函數(shù)I(x,t)表示在圖像空間點(diǎn)位置為x、時(shí)間為t時(shí)亮度，令I(lǐng)(x,0)=f(x)，經(jīng)過時(shí)間t之后產(chǎn)生動(dòng)作位移函數(shù)為δ(t)[23]，則：

I(x,t)=f(x+δ(t))

(1)

對(duì)于微小運(yùn)動(dòng)信號(hào)給予適當(dāng)放大因子放大，可對(duì)式(1)進(jìn)行關(guān)于x的一階泰勒級(jí)數(shù)展開分析，如式(2)所示：

(2)

令經(jīng)過時(shí)域帶通濾波所要放大動(dòng)作信號(hào)為：

(3)

將該信號(hào)放大α(放大因子)倍后，疊加到原信號(hào)I(x,t)，如式(4)所示：

(4)

由此可得：

(5)

將式(5)作關(guān)于一階泰勒級(jí)數(shù)展開逆方向運(yùn)算，得到逼近運(yùn)動(dòng)放大后亮度信號(hào)：

(6)

由式(1)與式(6)比較，可知經(jīng)過運(yùn)動(dòng)放大后，運(yùn)動(dòng)幅度是沒經(jīng)過放大的(1+α)倍。為了防止放大失真，對(duì)所構(gòu)造金字塔中不同空間頻率基帶應(yīng)使用合理放大因子進(jìn)行限制，假設(shè)信號(hào)空間波長(zhǎng)為λ=2π/ω，其中ω為信號(hào)空間頻率，則放大因子α限制公式[23]如式(7)所示：

(1+α)δ(t)<λ/8

(7)

經(jīng)過歐拉影像放大處理，橈動(dòng)脈跳動(dòng)部位對(duì)應(yīng)像素點(diǎn)的亮度值在時(shí)域上波動(dòng)較大，所以在時(shí)域上對(duì)視頻幀像素點(diǎn)亮度通道作亮度方差統(tǒng)計(jì)，以灰度值表示像素亮度方差值大小，灰度值越大代表該像素點(diǎn)亮度變化幅度越大。由于手臂背景中可能存在運(yùn)動(dòng)物體，在亮度方差統(tǒng)計(jì)中就會(huì)產(chǎn)生較大方差值，這會(huì)對(duì)后續(xù)步驟——橈動(dòng)脈定位產(chǎn)生很大影響，所以在進(jìn)行視頻幀亮度方差統(tǒng)計(jì)同時(shí)，對(duì)經(jīng)過歐拉影像放大的視頻幀作皮膚分割，提取視頻中皮膚區(qū)域，并將提取結(jié)果與亮度方差統(tǒng)計(jì)結(jié)果作邏輯與運(yùn)算，去除背景方差。去除背景后僅剩皮膚區(qū)域亮度方差，所以對(duì)該區(qū)域作形態(tài)學(xué)處理實(shí)現(xiàn)橈動(dòng)脈跳動(dòng)定位，并對(duì)所定位區(qū)域像素點(diǎn)根據(jù)其亮度方差大小計(jì)算該區(qū)域亮度加權(quán)平均值，獲取生理信號(hào)，對(duì)所獲信號(hào)進(jìn)行時(shí)頻轉(zhuǎn)換，實(shí)現(xiàn)心率非接觸式測(cè)量。

2 測(cè)量方法

本文運(yùn)用歐拉影像放大技術(shù)對(duì)橈動(dòng)脈微小跳動(dòng)動(dòng)作幅度放大，通過圖像處理技術(shù)對(duì)視頻幀皮膚分割以及方差統(tǒng)計(jì)，實(shí)現(xiàn)對(duì)橈動(dòng)脈跳動(dòng)部位提取，然后對(duì)視頻幀所提取區(qū)域亮度通道時(shí)頻分析，實(shí)現(xiàn)心率自動(dòng)測(cè)量，具體流程包括以下幾個(gè)部分：視頻輸入，橈動(dòng)脈跳動(dòng)微小變化放大，橈動(dòng)脈跳動(dòng)區(qū)域提取和心率計(jì)算。

2.1 橈動(dòng)脈跳動(dòng)微小變化放大

歐拉影像放大技術(shù)主要先由圖像RGB空間轉(zhuǎn)換到Lab空間，并在Lab空間中通過對(duì)視頻幀空間域和時(shí)域上處理，實(shí)現(xiàn)放大視頻中人眼難以察覺到的微小運(yùn)動(dòng)變化，其中包括手腕附近橈動(dòng)脈跳動(dòng)。對(duì)于微小運(yùn)動(dòng)放大空間域上處理主要是建立視頻幀拉普拉斯金字塔，原圖像經(jīng)過高斯低通濾波并以步長(zhǎng)為2隔行隔列采樣即可得到高斯金字塔第1層圖像(第0層為原圖像)，對(duì)所得圖像同樣濾波和采樣，迭代若干次可得到高斯金字塔，所以高斯金字塔把原圖像分解為低通子帶，在低通子帶依次對(duì)高斯金字塔上一層上采樣，再將采樣結(jié)果高斯濾波之后與高斯金字塔下層作差運(yùn)算即可得到拉普拉斯金字塔。在空間多分辨率分解過程中，結(jié)合采集視頻分辨率提高金字塔層數(shù)，可實(shí)現(xiàn)偽影抑制和運(yùn)動(dòng)放大，同時(shí)提高視頻信噪比。在時(shí)域上，運(yùn)用時(shí)域?yàn)V波器對(duì)每一視頻幀金字塔分辨率相同的圖像進(jìn)行時(shí)域?yàn)V波，提取視頻中感興趣微小運(yùn)動(dòng)信號(hào)。時(shí)域?yàn)V波對(duì)象是視頻幀中像素點(diǎn)，視頻中脈搏跳動(dòng)引起對(duì)應(yīng)像素值作同樣周期變化，所以可根據(jù)人體可能出現(xiàn)脈搏頻率范圍設(shè)定所需濾波器截止頻率。基于以上分析，對(duì)橈動(dòng)脈視頻幀運(yùn)用歐拉影像放大技術(shù)，實(shí)現(xiàn)橈動(dòng)脈跳動(dòng)幅度放大，算法流程如圖2所示。

圖2 歐拉影像放大方法流程 Fig. 2 Flowchart of Eulerian video magnification method

橈動(dòng)脈跳動(dòng)微小變化放大具體步驟如下：

步驟1 對(duì)輸入視頻每幀利用前文所述拉普拉斯金字塔進(jìn)行空間多分辨率分解，每幀對(duì)應(yīng)一個(gè)拉普拉斯金字塔，多分辨率金字塔中層次越高，空間頻率越低，信噪比越高。

步驟2 對(duì)每個(gè)拉普拉斯金字塔中空間分辨率相同的圖像作時(shí)域帶通濾波，時(shí)域?yàn)V波結(jié)果即是提取視頻L通道中脈搏信號(hào)，視頻幀時(shí)域?yàn)V波截止頻率由脈搏頻率范圍(0.3～3 Hz)決定。

步驟3 根據(jù)放大因子對(duì)步驟2濾波結(jié)果提取信號(hào)進(jìn)行放大，其中放大因子隨著空間頻率增大(信噪比的下降)而減小，即由金字塔頂層到底層逐漸減小，然后將放大后多分辨率金字塔合成圖像。

步驟4 物體微小位移主要與亮度有關(guān)，因此在圖像的Lab空間中L亮度通道放大脈搏跳動(dòng)信號(hào)時(shí)，抑制a、b通道信號(hào)，可提高輸出視頻信噪比。根據(jù)式(4)，放大后信號(hào)與原信號(hào)相加，將放大后圖像疊加到原視頻幀，完成圖像最后合成并轉(zhuǎn)換回RGB空間，重構(gòu)視頻。

步驟5 輸出視頻，可見視頻幀標(biāo)記位置時(shí)域上紋理變化比放大前明顯很多。

2.2 橈動(dòng)脈跳動(dòng)區(qū)域提取

橈動(dòng)脈跳動(dòng)微小變化放大之后，橈動(dòng)脈跳動(dòng)部分時(shí)域上亮度值變化幅度增大，根據(jù)放大后視頻幀像素點(diǎn)時(shí)域上亮度值的變化，可實(shí)現(xiàn)對(duì)橈動(dòng)脈跳動(dòng)區(qū)域提取。物體在運(yùn)動(dòng)過程中伴隨著視頻幀對(duì)應(yīng)位置像素點(diǎn)亮度通道上變化，所以經(jīng)過橈動(dòng)脈跳動(dòng)微小變化放大后，相對(duì)于皮膚其他區(qū)域，橈動(dòng)脈跳動(dòng)部位所對(duì)應(yīng)像素亮度值變化幅度明顯較大。因此對(duì)視頻橈動(dòng)脈跳動(dòng)前1～3個(gè)周期像素亮度值，根據(jù)方差計(jì)算公式作時(shí)域上像素點(diǎn)亮度方差統(tǒng)計(jì)，統(tǒng)計(jì)結(jié)果以亮度方差統(tǒng)計(jì)灰度圖表示，其中灰度值較高表示該像素點(diǎn)亮度變化幅度較大。同時(shí)，對(duì)橈動(dòng)脈跳動(dòng)微小變化放大后視頻幀進(jìn)行皮膚區(qū)域提取，提取結(jié)果與亮度方差統(tǒng)計(jì)灰度圖進(jìn)行邏輯與運(yùn)算，去除手臂背景方差，避免可能存在復(fù)雜運(yùn)動(dòng)背景和皮膚邊緣部分，對(duì)后續(xù)處理產(chǎn)生影響，以便進(jìn)一步確定橈動(dòng)脈跳動(dòng)位置。最后，根據(jù)去除背景方差圖，利用圖像形態(tài)學(xué)處理，進(jìn)行橈動(dòng)脈跳動(dòng)區(qū)域定位，從而實(shí)現(xiàn)橈動(dòng)脈跳動(dòng)區(qū)域提取。橈動(dòng)脈跳動(dòng)區(qū)域提取流程如圖3所示。

圖3 橈動(dòng)脈跳動(dòng)部位提取流程 Fig. 3 Flowchart of radial artery pulsing region extraction

2.2.1 統(tǒng)計(jì)方差

橈動(dòng)脈細(xì)微跳動(dòng)放大視頻幀中，結(jié)合人體可能出現(xiàn)的心率(20～180 beat/min)和采集視頻幀速率(n幀/s)，用橈動(dòng)脈至少一個(gè)跳動(dòng)周期對(duì)應(yīng)的視頻幀計(jì)算像素點(diǎn)時(shí)域上亮度方差，用到視頻幀數(shù)N=60n/20，根據(jù)式(8)計(jì)算像素點(diǎn)(x,y)時(shí)域上亮度方差σ2(x,y)：

)2

(8)

其中：Li(x,y)代表時(shí)域上像素點(diǎn)的亮度。

利用前N幀計(jì)算出每個(gè)像素點(diǎn)時(shí)域上亮度方差，令方差灰度圖矩陣為A，則A(x,y)=σ2(x,y)，A(x,y)越大表示像素亮度幅度波動(dòng)越大。

2.2.2 皮膚分割

對(duì)橈動(dòng)脈跳動(dòng)微小變化放大后視頻幀進(jìn)行皮膚分割，提取皮膚區(qū)域，可去除亮度方差統(tǒng)計(jì)灰度圖中運(yùn)動(dòng)背景影響，降低皮膚邊緣存在亮度變化干擾，使后續(xù)提取脈搏信號(hào)信噪比更高，因此有必要將視頻畫面中皮膚區(qū)域進(jìn)行分割提取。在圖像顏色空間YCrCb空間中，Y是亮度，Cr、Cb分別是紅色、藍(lán)色的色度，為了對(duì)亮度變化具有較高適應(yīng)性，YCrCb空間對(duì)皮膚檢測(cè)具有較高效率[24]，所以通過在Cr、Cb通道設(shè)置閾值實(shí)現(xiàn)皮膚分割，符合式(9)認(rèn)為是皮膚區(qū)域[24]。對(duì)不同光照亮度皮膚分割結(jié)果如圖4所示，其中(d)、(e)、( f)是不同光照亮度圖像(a)、(b)、(c)對(duì)應(yīng)的皮膚二值區(qū)域，即感興趣區(qū)域(Region Of Interest, ROI)，可見YCrCb空間對(duì)皮膚檢測(cè)具有較高魯棒性。

(9)

本文參考文獻(xiàn)[24]方法，并針對(duì)背景中可能出現(xiàn)與皮膚像素值接近像素的誤判，結(jié)合本文具體應(yīng)用改進(jìn)。經(jīng)過式(9)閾值后產(chǎn)生二值圖像，對(duì)該二值圖像在閾值分割后進(jìn)行連通域大小計(jì)算，并選擇面積最大連通域，對(duì)最大連通域中可能存在孔洞經(jīng)行檢測(cè)和填充，并進(jìn)行邊緣腐蝕從而分割提取出皮膚區(qū)域。皮膚分割結(jié)果RGB圖像和二值圖像如圖5所示。

圖4 不同光照下皮膚區(qū)域初步提取結(jié)果 Fig. 4 Initial extraction results of skin area under different light conditions

圖5 皮膚分割結(jié)果 Fig. 5 Results with skin segmentation

2.2.3 橈動(dòng)脈跳動(dòng)區(qū)域定位

亮度方差灰度圖與皮膚分割二值圖像進(jìn)行邏輯與運(yùn)算后，得到僅有皮膚區(qū)域亮度方差的方差灰度圖。通過對(duì)所提取的皮膚區(qū)域方差灰度圖作閾值分割、孔洞填充，腐蝕、最大連通域選擇和膨脹，最后，提取最大連通域外接矩形區(qū)域作為生理信號(hào)區(qū)域，其中閾值分割中閾值由經(jīng)驗(yàn)值和皮膚區(qū)域方差值范圍作動(dòng)態(tài)設(shè)置，具體步驟如下。

步驟1 將統(tǒng)計(jì)方差灰度圖矩陣與皮膚分割二值圖像矩陣作點(diǎn)乘運(yùn)算，消除皮膚背景方差灰度，避免背景中出現(xiàn)較高方差灰度值對(duì)橈動(dòng)脈跳動(dòng)區(qū)域定位造成影響。

步驟2 對(duì)皮膚區(qū)域方差灰度圖作動(dòng)態(tài)閾值分割，閾值取值如式(10)所示：

(10)

步驟3 對(duì)步驟2 閾值分割保留下來區(qū)域?qū)?yīng)的灰度進(jìn)行灰度腐蝕，使灰度最高連通域相對(duì)最大化，灰度較低連通域最小化；然后進(jìn)行連通域大小計(jì)算，并選擇最大連通域作為膨脹對(duì)象，將其恢復(fù)到腐蝕前大小；最后在亮度通道對(duì)該連通域取外接矩形，并將此矩形中像素點(diǎn)作為橈動(dòng)脈跳動(dòng)部位提取結(jié)果。

2.3 心率計(jì)算

經(jīng)過上述步驟處理得到橈動(dòng)脈跳動(dòng)矩形區(qū)域，心率計(jì)算將對(duì)該區(qū)域中像素亮度求加權(quán)平均值，并將這些加權(quán)平均值構(gòu)成信號(hào)作快速傅里葉變換得到傅里葉能量譜，將能量譜中能量最高頻率乘于60，即可實(shí)現(xiàn)心率計(jì)算。

由于最后提取橈動(dòng)脈矩形區(qū)域中方差值大小不一，即橈動(dòng)脈跳動(dòng)部對(duì)應(yīng)像素點(diǎn)亮度波動(dòng)幅度存在差異，直接對(duì)矩形中像素點(diǎn)亮度取均值作為該幀脈搏亮度值，會(huì)造成生理信號(hào)信噪比下降。因此，根據(jù)橈動(dòng)脈跳動(dòng)矩形部位方差灰度圖中灰度值大小，對(duì)采集視頻幀的橈動(dòng)脈跳動(dòng)部位矩陣像素點(diǎn)亮度值求加權(quán)平均值，加權(quán)平均值在時(shí)域上形成計(jì)算心率所需脈搏波。亮度加權(quán)平均值計(jì)算和基于傅里葉能量譜分析心率計(jì)算具體步驟如下。

步驟1 橈動(dòng)脈跳動(dòng)矩形區(qū)域中像素亮度求加權(quán)平均值，權(quán)值wi如式(11)所示：

(11)

步驟2 對(duì)脈搏波信號(hào)作快速傅里葉變換，結(jié)合人體可能出現(xiàn)的心率，選取頻率帶寬為0.3～3 Hz(20～180 beat/min)作為分析頻段，并在此頻段范圍內(nèi)選擇能量最高點(diǎn)對(duì)應(yīng)的頻率作為心率測(cè)量結(jié)果，生理信號(hào)頻域能量譜如圖6(b)所示。

步驟3 由步驟3中傅里葉能量譜可得能量最高所對(duì)應(yīng)的頻率為fh=1.400 Hz，心率HR可由下式求出：HR=fh×60=1.400×60=84 beat/min。

3 實(shí)驗(yàn)與分析

3.1 數(shù)據(jù)收集與處理

本文實(shí)驗(yàn)中，利用攝像頭能清晰采集到橈動(dòng)脈部位圖像，要求拍攝過程中手避免大幅度運(yùn)動(dòng)。設(shè)定采集原始視頻圖像空間為RGB，幀率為30幀/s，拍攝時(shí)長(zhǎng)10～15 s，獲得300～450幀圖像。將采集到視頻直接作歐拉影像放大處理，實(shí)現(xiàn)橈動(dòng)脈跳動(dòng)微小變化放大。如圖7，橈動(dòng)脈跳動(dòng)微小變化放大前人眼難以觀察到的橈動(dòng)脈微小跳動(dòng)，脈動(dòng)放大后顯示出明顯凹陷和凸起，如圖7中圓圈區(qū)域所示。

圖6 原始信號(hào)與傅里葉能量譜 Fig. 6 Original signal and Fourier energy spectrum

圖7 歐拉影像放大前后對(duì)比 Fig. 7 Comparison of before and after Eulerian video magnification

3.2 橈動(dòng)脈跳動(dòng)區(qū)域定位

圖8 橈動(dòng)脈跳動(dòng)區(qū)域定位 Fig. 8 Positioning for pulsing region of radial artery

3.3 心率計(jì)算

對(duì)所提取橈動(dòng)脈區(qū)域，在亮度通道上根據(jù)式(11)計(jì)算每幀對(duì)應(yīng)加權(quán)平均值，數(shù)據(jù)歸一化后可得到如圖9(a)的脈搏波，并對(duì)該生理信號(hào)作快速傅里葉變換，脈搏波傅里葉能量譜如圖9(b)所示，然后選取能量譜中最高能量所對(duì)應(yīng)頻率計(jì)算心率：心率HR可由下式求出：HR=fh×60=1.166 7×60≈70 beat/min。

圖9 脈搏波與能量譜 Fig. 9 Pulse wave and energy spectrum

3.4 結(jié)果分析

由于本文方法主要是基于橈動(dòng)脈視頻提取心率，所以在不同強(qiáng)度光照環(huán)境中橈動(dòng)脈所在位置光照度對(duì)實(shí)驗(yàn)結(jié)果精確度造成一定影響。為了驗(yàn)證本文方法，借助光照度計(jì)，在日常室內(nèi)不同光照度條件下獲取橈動(dòng)脈視頻進(jìn)行心率計(jì)算精度(日常室內(nèi)光照度一般范圍為100～1 000 lux)，按照上述實(shí)驗(yàn)方法步驟，在光照度為20～1 000 lux范圍內(nèi)由光照度計(jì)測(cè)得光照度值并選取10個(gè)不同光照度作為實(shí)驗(yàn)光照環(huán)境，對(duì)同一測(cè)試志愿者進(jìn)行非接觸式心率測(cè)量；同時(shí)進(jìn)行脈搏血氧儀測(cè)量，其中每一相同光照度進(jìn)行5組實(shí)驗(yàn)，并以脈搏血氧儀測(cè)量結(jié)果作為對(duì)比標(biāo)準(zhǔn)，計(jì)算在不同光照度下心率測(cè)量平均絕對(duì)誤差和最大絕對(duì)誤差，不同光照度(E)實(shí)驗(yàn)誤差曲線如圖10所示。當(dāng)光照度E≤80 lux時(shí)，由于采集的視頻質(zhì)量較差，直接導(dǎo)致心率測(cè)量結(jié)果誤差比較大，當(dāng)光照度逐漸充足時(shí)，心率測(cè)量結(jié)果誤差明顯下降并趨于穩(wěn)定，符合中華人民共和國(guó)醫(yī)藥行業(yè)標(biāo)準(zhǔn)(誤差 ≤5 beat/min)的要求。

圖10 光照度實(shí)驗(yàn)誤差曲線 Fig. 10 Error curve of difference lux

為了進(jìn)一步驗(yàn)證本文方法準(zhǔn)確性，按照本文實(shí)驗(yàn)方法步驟，在本實(shí)驗(yàn)對(duì)10個(gè)測(cè)試志愿者進(jìn)行非接觸式心率測(cè)量，并將測(cè)量結(jié)果與脈搏血氧儀同步測(cè)量結(jié)果對(duì)比分析。其中，10個(gè)測(cè)試者分為10組，即每組1人，每人非接觸式和同步接觸式心率測(cè)量各進(jìn)行5次。對(duì)得到50對(duì)數(shù)據(jù)采用Bland-Altman 法對(duì)本次實(shí)驗(yàn)與脈搏血氧儀測(cè)量進(jìn)行一致性評(píng)估。

如圖11所示，坐標(biāo)上點(diǎn)S的數(shù)學(xué)表達(dá)式表示為：

S(x,y)=((s1+s2)/2,(s1-s2))

(13)

其中s1，s2分別為同一次實(shí)驗(yàn)中用本文方法和脈搏血氧儀測(cè)量心率值，橫坐標(biāo)x表示同一次測(cè)量中兩種方法測(cè)量平均值，縱坐標(biāo)y表示同一次測(cè)量中兩種方法測(cè)量值差值，置信率取95%(1.96標(biāo)準(zhǔn)差)。

圖11 Bland-Altman法估計(jì)兩種測(cè)量方法結(jié)果的一致性 Consistency estimation of results of two methods by Bland-Altman method

圖11的Bland-Altman 法一致性估計(jì)，對(duì)50對(duì)測(cè)量數(shù)據(jù)作符合度分析，誤差均值為0.120 0，標(biāo)準(zhǔn)方差為2.291 4，95%的置信區(qū)間為[-4.371,4.611]，可說明本文心率測(cè)量方法與脈搏血氧儀測(cè)量法具有很好一致性。對(duì)50對(duì)測(cè)量數(shù)據(jù)從平均誤差、標(biāo)準(zhǔn)方差和相關(guān)系數(shù)將本文與同樣采用Bland-Altman 法一致性估計(jì)的文獻(xiàn)[11,16]方法比較，對(duì)比結(jié)果如表1所示。同樣是將測(cè)量結(jié)果與脈搏血氧儀結(jié)果對(duì)比，與前兩種方法相比，本文測(cè)量方法的平均誤差絕對(duì)值和標(biāo)準(zhǔn)方差相對(duì)較小，標(biāo)準(zhǔn)方差相對(duì)于前兩者分別下降50.5%和32.6%，說明本文方法心率測(cè)量誤差明顯下降。另外，本文方法測(cè)量結(jié)果與脈搏血氧儀結(jié)果測(cè)量結(jié)果相關(guān)系數(shù)0.963 1，介于前兩種算法之間，經(jīng)實(shí)驗(yàn)分析，影響誤差以及相關(guān)系數(shù)的主要原因是測(cè)量過程中手腕部位移動(dòng)幅度較大。實(shí)驗(yàn)的10組數(shù)據(jù)的誤差具體如表2所示。

表1 各類方法測(cè)量結(jié)果對(duì)比Tab. 1 Measurement results comparison of various algorithms

表2 本文方法實(shí)驗(yàn)測(cè)量誤差 beat/minTab. 2 Error of experimental measurement for the proposed method beat/min

實(shí)驗(yàn)計(jì)算得出10組誤差平均絕對(duì)誤差為1.76 beat/min，平均最大絕對(duì)誤差為3.7 beat/min，較近國(guó)內(nèi)相關(guān)研究文獻(xiàn)[14]平均絕對(duì)誤差2.0 beat/min，平均最大絕對(duì)誤差為4.0 beat/min，與之相比，本文測(cè)量方法平均絕對(duì)誤差下降12%，平均最大絕對(duì)誤差下降7.5%。根據(jù)中華人民共和國(guó)醫(yī)藥行業(yè)標(biāo)準(zhǔn)(誤差≤5 beat/min)的要求，本文方法滿足心率測(cè)量應(yīng)用，適用于個(gè)人家庭應(yīng)用或遠(yuǎn)程醫(yī)療應(yīng)用。

4 結(jié)語(yǔ)

本文針對(duì)接觸式心率測(cè)量所帶來不便，盲分離信號(hào)準(zhǔn)確性缺乏判斷，以及基于PPG原理非接觸式心率測(cè)量受外部環(huán)境溫度影響較大等不足，提出一種新的心率非接觸式測(cè)量方法。實(shí)驗(yàn)結(jié)果表明，本文心率測(cè)量方法在日常室內(nèi)光照環(huán)境中具有較好的適應(yīng)性，且與脈搏血氧儀測(cè)量法具有很好一致性，與基于PPG原理的測(cè)量方法相比，有效地降低了測(cè)量誤差，符合中華人民共和國(guó)醫(yī)藥行業(yè)標(biāo)準(zhǔn)，可達(dá)到預(yù)期心率測(cè)量效果，在日常家庭保健和遠(yuǎn)程醫(yī)療上應(yīng)用具有較大發(fā)展前景。經(jīng)過實(shí)驗(yàn)分析，本文測(cè)量誤差主要是測(cè)量過程中受測(cè)對(duì)象手腕發(fā)生幅度較大位移，下一步的主要工作將從受測(cè)對(duì)象手腕脈搏處提取相關(guān)紋理特征，并對(duì)其進(jìn)行識(shí)別跟蹤，有效降低心率測(cè)量誤差。

參考文獻(xiàn)(References)

[1] STONE M L, TATUM P M, WEITKAMP J-H, et al. Abnormal heart rate characteristics before clinical diagnosis of necrotizing enterocolitis [J]. Journal of Perinatology Official Journal of the California Perinatal Association, 2013, 33(11): 847-850.

[2] OGLIARI G, MAHINRAD S, STOTT D J, et al. Resting heart rate, heart rate variability and functional decline in old age [J]. Canadian Medical Association Journal, 2015, 187(15): E442-E449.

[3] GRENEKER E F. Radar sensing of heartbeat and respiration at a distance with applications of the technology [C]// RADAR 97: Proceedings of the 1997 Radar Systems. [S.l.]: The Institution of Engineering and Technology, 1997: 150-154.

[4] HIGASHIKATURAGI K, NAKAHATA Y, MATSUNAMI I, et al. Non-invasive respiration monitoring sensor using UWB-IR [C]// ICUWB 2008: Proceedings of the 2008 IEEE International Conference on Ultra-Wideband. Piscataway, NJ: IEEE, 2008: 29-34.

[5] LAZARO A, GIRBAU D, VILLARINO R. Analysis of vital signs monitoring using an Ir-UWB radar [J]. Progress in Electromagnetics Research, 2010, 100: 265-284.

[6] GARBEY M, SUN N, MERLA A, et al. Contact-free measurement of cardiac pulse based on the analysis of thermal imagery [J]. IEEE Transactions on Biomedical Engineering, 2007, 54(8):1418-1426.

[7] TRESKES R W, van de VELDE V, ATSMA D E, et al. Redesigning healthcare: The 2.4 billion euro question?: Connecting smart technology to improve outcome of patients [J]. Netherlands Heart Journal: Monthly Journal of the Netherlands Society of Cardiology & the Netherlands Heart Foundation, 2016, 24(7/8): 441-446.

[8] TAKANO C, OHTA Y. Heart rate measurement based on a time-lapse image [J]. Medical Engineering & Physics, 2007, 29(8): 853-857.

[9] ZHENG J, HU S, CHOULIARAS V, et al. Feasibility of imaging photoplethysmography [C]// BMEI ’08: Proceedings of the 2008 International Conference on Biomedical Engineering and Informatics. Washington, DC: IEEE Computer Society, 2008, 2: 72-75.

[10] VERKRUYSSE W, SVAASAND L O, NELSON J S. Remote plethysmographic imaging using ambient light [J]. Optics Express, 2008, 16(26):21434-21445.

[11] POH M Z, McDUFF D J, PICARD R W. Non-contact, automated cardiac pulse measurements using video imaging and blind source separation [J]. Optics Express, 2010, 18(10): 10762-10774.

[12] POH M Z, McDUFF D J, PICARD R W. Advancements in noncontact, multiparameter physiological measurements using a webcam [J]. IEEE Transactions on Bio-medical Engineering, 2011, 58(1):7-11.

[13] PURSCHE T, KRAJEWSKI J, MOELLER R. Video-based heart rate measurement from human faces [C]// ICCE 2012: Proceedings of the 2012 IEEE International Conference on Consumer Electronics. Piscataway, NJ: IEEE, 2012: 544-545.

[14] 劉祎,歐陽(yáng)健飛,閆勇剛.基于普通攝像頭的心率測(cè)量方法研究[J].計(jì)算機(jī)工程與應(yīng)用,2016,52(7):210-214. (LIU Y, OUYANG J F, YAN Y G. Method of measuring heart rate using a webcam [J]. Computer Engineering and Applications, 2016, 52(7)：210-214.)

[15] ZHAO F, LI M, QIAN Y, et al. Remote measurements of heart and respiration rates for telemedicine [J]. PLOS ONE, 2013, 8(10): e71384.

[16] XU S, SUN L, ROHDE G K. Robust efficient estimation of heart rate pulse from video [J]. Biomedical Optics Express, 2014, 5(4): 1124-35.

[17] YU Y-P, KWAN B-H, LIM C-L, et al. Video-based heart rate measurement using short-time Fourier transform [C]// ISPACS 2013: Proceedings of the 2013 International Symposium on Intelligent Signal Processing and Communications Systems. Piscataway, NJ: IEEE, 2013: 704-707.

[18] YU Y-P, RAVEENDRAN P, LIM C-L. Dynamic heart rate measurements from video sequences [J]. Biomedical Optics Express, 2015, 6(7): 2466-2480.

[19] LI X, CHEN J, ZHAO G, et al. Remote heart rate measurement from face videos under realistic situations [C]// CVPR ’14: Proceedings of the 2014 IEEE Conference on Computer Vision and Pattern Recognition. Piscataway, NJ: IEEE, 2014: 4264-4271.

[20] LAM A, KUNO Y. Robust heart rate measurement from video using select random patches [C]// ICCV ’15: Proceedings of the 2015 IEEE International Conference on Computer Vision. Washington, DC: IEEE Computer Society, 2015: 3640-3648.

[21] ZHANG X Y, ZHANG Y T. The effect of local mild cold exposure on pulse transit time [J]. Physiological Measurement, 2006, 27(7): 649-660.

[22] KHAN M, PRETTY C G, AMIES A C, et al. Analysing the effects of cold, normal, and warm digits on transmittance pulse oximetry [J]. Biomedical Signal Processing & Control, 2016, 26: 34-41.

[23] WU H-Y, RUBINSTEIN M, SHIH E, et al. Eulerian video magnification for revealing subtle changes in the world [J]. ACM Transactions on Graphics, 2012, 31(4): Article No. 65.

[24] MAHMOUD T M. A new fast skin color detection technique [J]. World Academy of Science: Engineering & Technolog, International Journal of Computer and Information Engineering, 2008, 2(7): 2354-2358.

This work is partially supported by the National Natural Science Foundation of China (21376091), the Guangdong Provincial Science and Technology Plan Project (2015A030401089, 2015B090922004, 2016B090927004).

SUPeiquan, born in 1991, M. S. candidate. His research interests include machine vision, machine learning.

XULiang, born in 1971, Ph.D., senior engineer. His research interests include machine vision, machine learning, wireless sensor network.

LIANGYongjian, born in 1991, M. S. candidate. His research interests include machine vision, machine learning.