基于修正加速度的對數(shù)歸一化變步長自適應(yīng)濾波的心率估計(jì)算法

譚 擁， 余成波, 張 林

(重慶理工大學(xué)電氣與電子工程學(xué)院， 重慶 401320)

光電容積脈搏波(photoplethysmography，PPG)描記圖法的原理是利用心臟泵血導(dǎo)致外周動(dòng)脈發(fā)生血液容積變化，通過特定波長的光去檢測這種變化，將這種變化以反射光強(qiáng)的方式記錄下來。與心電圖(electrocardiogram，ECG)不同，PPG測量設(shè)備不需要緊貼皮膚，減少了佩戴者的不適感，但這也使得PPG信號易受到因人體運(yùn)動(dòng)而產(chǎn)生的運(yùn)動(dòng)偽影(motion artifact，MA)的影響，致使PPG信號質(zhì)量變差，提取相關(guān)信息困難，這也是目前測量PPG信號面臨最大的挑戰(zhàn)。

針對上述問題，研究者們提出許多去運(yùn)動(dòng)偽影方法。傳統(tǒng)的頻域分析法，如快速傅里葉變換[1]，通過能量分析來濾除運(yùn)動(dòng)噪聲，繼而獲得心率。但該方法只適用于靜止或噪聲較小的情況。Zhang等[2]提出了TRIOKA框架，該方法利用帶通濾波、奇異值分解，對奇異值聚類分組，以加速度作為參考信號去除相關(guān)噪聲，從而使信號稀疏化，最后利用稀疏信號重構(gòu)(sparse signal reconstruction，SSR)獲得一個(gè)高分辨率的PPG信號，對一般運(yùn)動(dòng)下的PPG信號能夠獲得較為準(zhǔn)確的心率，然而該方法復(fù)雜度較高，易造成冗余,資源浪費(fèi)，不適用于長時(shí)間監(jiān)測。孫斌等[3]學(xué)者統(tǒng)計(jì)靜止?fàn)顟B(tài)下PPG信號的峰度、偏度等指標(biāo)，以此判斷PPG信號是否被污染，并將污染的PPG信號片段剪切，以達(dá)到去除運(yùn)動(dòng)偽影的目的，雖然降低了復(fù)雜度，但去除運(yùn)動(dòng)偽影的同時(shí)也把有用信息剪切掉。Dragomiretskiy等[4]學(xué)者在經(jīng)驗(yàn)?zāi)B(tài)分解的基礎(chǔ)上提出一個(gè)完全非遞歸的變分模態(tài)分解(variational mode decom-position，VMD)模型，有效的解決了模態(tài)混疊的問題。Sharma[5]使用VMD對PPG信號進(jìn)行分解、重構(gòu)，但只對靜止?fàn)顟B(tài)下的PPG信號進(jìn)行了研究。此外，還有其他抗運(yùn)動(dòng)偽影方法，如集合經(jīng)驗(yàn)?zāi)B(tài)分解、自適應(yīng)濾波法[6-7]等，均未在復(fù)雜度和準(zhǔn)確率之間尋找到折中點(diǎn)。

運(yùn)動(dòng)偽影主要由外部運(yùn)動(dòng)所造成的，因此研究者將加速度[2,6-7]作為去除MA的參考信號。然而加速度計(jì)是基于慣性參考系[8-9]測量，得到的加速度信號會(huì)受到重力分量的影響。Jarchi等[10]在采集多種運(yùn)動(dòng)下PPG信號時(shí)，還采集了加速度以及三軸角速度。利用陀螺儀測出的三軸角速度將加速度映射到傳感器參考系[8]，從而修正加速度，獲得更準(zhǔn)確地心率估計(jì)值。

結(jié)合加速度和運(yùn)動(dòng)偽影的特點(diǎn)，提出一種基于修正加速度的對數(shù)歸一化變步長自適應(yīng)濾波算法。該算法對PPG信號進(jìn)行0.4～5 Hz的帶通濾波，用三軸角速度對三軸加速度進(jìn)行卡爾曼濾波修正，然后將修正后的加速度作為Log-NLMS自適應(yīng)濾波器的參考信號，濾除運(yùn)動(dòng)偽影，從而得到高質(zhì)量的PPG信號，最后，通過頻域的譜峰準(zhǔn)追蹤來估計(jì)心率值。此外，從同步采集的ECG信號中獲取真實(shí)心率，來驗(yàn)證本文算法的抗運(yùn)動(dòng)干擾能力與心率估計(jì)的準(zhǔn)確性。

1 心率估計(jì)算法

所提出的基于修正的加速度的運(yùn)動(dòng)偽影去除算法主要包含：數(shù)據(jù)預(yù)處理、加速度的修正、PPG信號的自適應(yīng)濾波和譜峰追蹤，其流程圖如圖1所示。

1.1 數(shù)據(jù)預(yù)處理

根據(jù)人體正常心率范圍，首先對原PPG信號進(jìn)行0.4～5 Hz的帶通濾波[2,11]。去除該范圍以外的運(yùn)動(dòng)噪聲、工頻干擾以及基線漂移，不僅達(dá)到使信號稀疏化的目的，而且還消除了多余的噪聲干擾，為接下來的信號自適應(yīng)濾波創(chuàng)造更有利的去噪條件。

1.2 加速度修正

在運(yùn)動(dòng)狀態(tài)下，加速度計(jì)測量的慣性參考系下加速度信號是由重力加速度和外部加速度的共同作用的結(jié)果，如式(1)所示，這樣會(huì)導(dǎo)致外部加速度測量誤差較大。

YA=Sg+Sa+nA

(1)

式(1)中：上標(biāo)S表示慣性參考系；下標(biāo)A表示加速度計(jì)；YA、nA分別表示加速度計(jì)測量值和測量誤差；Sg為重力分量；Sa為外部加速度分量。

根據(jù)文獻(xiàn)[8]基于旋轉(zhuǎn)的過程，將慣性參考系下的重力映射到傳感器參考系下，消除重力分量影響，修正加速度，可分別表示為

Sg=g×SZ

(2)

SZ=[-sinβcosβsinγcosβcosγ]

(3)

式中：β、γ分別為俯仰角和翻滾角，可由陀螺儀獲得；g為重力加速度，9.8 m/s2；SZ為旋轉(zhuǎn)向量。

利用六軸加速度傳感器獲得的加速度和角速度。然而，陀螺儀測量的角速度需要通過積分才能獲得俯仰角和翻滾角，由測量誤差積分所產(chǎn)生的誤差是令人無法接受的，因此陀螺儀不能用于長時(shí)間的角度測量，需要隔一段時(shí)間對陀螺儀測量姿態(tài)角進(jìn)行修正。這種誤差可以通過卡爾曼濾波的方法來減小[10,12]，從而有效將傳感器參考系下的加速度中的重力分量去除。采用線性卡爾曼濾波，其定義為

Xk=AXk-1+wk

(4)

zk=HXk-1+vk

(5)

式中：A為狀態(tài)轉(zhuǎn)移矩陣；H為測量矩陣；Xk為k時(shí)刻的狀態(tài)向量；zk為測量值；wk為均值為0、協(xié)方差矩陣為Q、且服從正態(tài)分布的過程噪聲；vk為均值為0，協(xié)方差矩陣為R，且服從正態(tài)分布的測量噪聲。利用捷聯(lián)式積分，式(1)、式(2)可將狀態(tài)轉(zhuǎn)移方程和測量方程描述為

SZt=[I-dt(YG,t-1×)]SZt-1-dt(SZt-1×)nG

(6)

(7)

(8)

(9)

計(jì)算卡爾曼增益Kk可表示為

(10)

應(yīng)用卡爾曼增益對當(dāng)前估計(jì)值進(jìn)行修正：

(11)

(12)

圖2 原始PPG信號、原始加速度及修正加速度FFT頻對比

1.3 Log-NLMS變步長自適應(yīng)濾波

歸一化最小均方差算法(normalized least mean square，NLMS)是一種改進(jìn)的自適應(yīng)濾波算法，它有效解決了最小均方差(least mean square，LMS)算法在收斂速度和誤差穩(wěn)態(tài)相互矛盾的問題[6]。該算法用步長和輸入信號的功率之比來對原有的固定步長進(jìn)行歸一化處理，在迭代計(jì)算中調(diào)整步長以及抽頭系數(shù)，使均方誤差快速達(dá)到最小，這樣既保證了收斂速度，又使得收斂誤差穩(wěn)定在可接受范圍內(nèi)。NLMS算法的步長更新公式為

(13)

(14)

式(14)中：c為一個(gè)極小的常數(shù)，其作用是避免步長分母為零；W(n)為n時(shí)刻的抽頭系數(shù)；e(n)為n時(shí)刻期望信號與輸入信號之間誤差。

對數(shù)函數(shù)的最小均方算法(log-least mean square，Log-LMS)是以減小穩(wěn)態(tài)誤差，提高收斂速為前提，對傳統(tǒng)的sigmoid函數(shù)[13]、反正切函數(shù)[14]的變步長LMS進(jìn)行改進(jìn)的算法。該算法將步長和誤差設(shè)置為對數(shù)關(guān)系，有效地解決了收斂速度和誤差穩(wěn)態(tài)之間的矛盾[15]。其步長更新表達(dá)式為

μ(n)′=blg[a|e(n)|m]

(15)

式(15)中：a表示控制曲線的整體形狀走向；m表示控制曲線底部變化速度；b表示控制曲線的幅度大小。

NLMS和Log-LMS雖然在一定程度改善了傳統(tǒng)LMS算法的缺點(diǎn)，但NLMS使用了全局步長因子，在一定意義上仍會(huì)存在傳統(tǒng)LMS算法同樣的問題，而Log-LMS算法對輸入信號比較敏感[16]。將兩種方法的優(yōu)點(diǎn)進(jìn)行融合，其步長表達(dá)式為

(16)

將修正的加速度信號作為自適應(yīng)濾波輸入， 設(shè)置參數(shù)濾波器階數(shù)L=3，m=2，a=1 000，b=0.001。如圖3所示，Log-NLMS分別將修正前后的加速度作為參考信號所得到的濾波后PPG頻譜進(jìn)行對比。圖3(a)為原始PPG信號的頻譜圖，清楚顯示了心率和MA兩個(gè)主峰，圖3(b)為原始加速度信號作為參考信號的Log-NLMS濾波。由于原始加速度主頻與圖3(a)中的MA主頻有所偏差，無法去除與加速度有關(guān)的MA，以至于接下來的譜峰追蹤無法定位真實(shí)心率位置。如圖3(c)所示，修正后的加速度作為參考信號，MA主頻與圖3(a)的MA主頻重合，能夠較為有效地消除MA，得到準(zhǔn)確地心率估計(jì)值。

圖3 PPG經(jīng)兩種算法濾波后的FFT頻對比

1.4 譜峰追蹤

經(jīng)過濾波處理后的PPG信號在時(shí)域上的波形如圖4(c)所示，雖然在前半段波峰不太突出，但整體上還是能和ECG波峰一一對應(yīng)，而且接下來的譜峰追蹤是從頻域上對圖3(c)進(jìn)行心率估計(jì)，就更加保證心率估計(jì)的正確性。

在心率估計(jì)算法中，譜峰追蹤是算法中一個(gè)重要組成部分。其要求主要包含兩點(diǎn)：①Q(mào)1：人類的心率短時(shí)間內(nèi)不會(huì)發(fā)生突變；②Q2：PPG頻譜中存在與心率對應(yīng)的諧波峰，尤其是二次諧波，與基波形成峰值對。譜峰追蹤具體如下。

(1)初始化：在該階段，首先要求測試者保持靜止不動(dòng)，以獲得靜止時(shí)的心率來作為下一個(gè)時(shí)間窗口心率估計(jì)的先驗(yàn)信息。

(2)峰值選擇：用Nprev表示上個(gè)時(shí)間窗口心率在頻譜中波峰所對應(yīng)的索引值。根據(jù)Q1可得，存在一個(gè)Δs使得當(dāng)前時(shí)間窗口的心率索引在R0=[Nprev-Δs,Nprev-Δs+1,…,Nprev+Δs]這個(gè)范圍內(nèi)，而對應(yīng)的二次諧波則在R1=[2(Nprev-Δs-1)+1,2(Nprev-Δs)+1,…, 2(Nprev+Δs-1)+1]內(nèi)，在(R0,R1)搜索峰值對。注意，在對應(yīng)區(qū)域中只尋找該區(qū)域內(nèi)大于該區(qū)域內(nèi)最大峰值的30%的峰[2]。

(3)驗(yàn)證：為了防止錯(cuò)誤追蹤，導(dǎo)致心率估計(jì)錯(cuò)誤，因此驗(yàn)證結(jié)果可靠性是有必要的。根據(jù)Q1，規(guī)定相鄰時(shí)間窗口心率不允許相差10 BPM(BPM為每秒鐘心臟搏動(dòng)次數(shù))。因此需要對心率進(jìn)行驗(yàn)證，其表達(dá)式為

(17)

式(17)中：N-表示當(dāng)前時(shí)間窗口選擇的峰值索引；Ncur表示驗(yàn)證后的峰值索引；α=6是根據(jù)頻譜劃分網(wǎng)格數(shù)N=5 120和采樣頻率Fs=256 Hz所決定。參數(shù)a根據(jù)具體情況設(shè)置，設(shè)a=2。

圖4 ECG信號與自適應(yīng)濾波后的PPG信號的對比

2 結(jié)果與討論

2.1 實(shí)驗(yàn)數(shù)據(jù)

圖5 PPG采集裝置實(shí)物圖

主要以模塊化的方式采集數(shù)據(jù)。如圖5所示，PPG采集裝置以STM32F103為主控芯片，控制PPG信號采集、六軸加速度信號采集、心電信號采集以及藍(lán)牙數(shù)據(jù)發(fā)送。模擬手環(huán)戴在手腕上，使用波長535 nm的綠光從手腕采集PPG信號，同時(shí)，從六軸加速度采集模塊和心電采集模塊同步獲取六軸加速度信號和心電信號。共采集了常見3種狀態(tài)：步行、低速騎行、高速騎行。每一種狀態(tài)都是先保持30 s靜止，接下的3 min實(shí)施相應(yīng)的動(dòng)作，最后保持30 s靜止，每實(shí)施一種動(dòng)作需要間隔10 min。實(shí)驗(yàn)分別對年齡在20～50歲的3名和2名女性測試者進(jìn)行了數(shù)據(jù)采集，每名測試者采集3組數(shù)據(jù)，分別對應(yīng)3種不同運(yùn)動(dòng)狀態(tài)。

2.2 評價(jià)指標(biāo)

實(shí)驗(yàn)采用平均絕對誤差(average absolute error，AAE)和平均相對誤差(average relative error，ARE)作為估計(jì)心率與真實(shí)心率的評價(jià)指標(biāo)，同時(shí)還采用Bland-Altman[17]圖來驗(yàn)證兩者之間的一致性。平均絕對誤差定義為

(18)

式(18)中：M表示總的時(shí)間窗口個(gè)數(shù)；HRECG(i)、HRPPG(i)分別表示第i個(gè)窗口的真實(shí)心率和估計(jì)心率。平均相對誤差定義為

(19)

Bland-Altman圖用來檢驗(yàn)估計(jì)值與真實(shí)值之間的一致性，其橫坐標(biāo)為兩者的均值，縱坐標(biāo)為兩者之間的差值，以此來展示每個(gè)估計(jì)值和對應(yīng)的真實(shí)值之間的差異。Bland-Altman圖劃定了一個(gè)置信區(qū)[Mean-1.96σ, Mean+1.96σ]，其中Mean表示所有估計(jì)心率與真實(shí)心率差值的均值，σ表示估計(jì)心率與真實(shí)心率差值的標(biāo)準(zhǔn)差。如果95%及以上數(shù)據(jù)點(diǎn)的在這個(gè)區(qū)域內(nèi)，表示置信度高，一致性好。

2.3 實(shí)驗(yàn)結(jié)果

PPG采集裝置的采樣頻率統(tǒng)一設(shè)置為256 Hz。用一個(gè)長為T=10 s的時(shí)間窗口以步長S=2 s在每組信號上滑動(dòng)。采用原始加速度作為參考信號的Log-NLMS(對比算法)和修正加速度作為參考信號的Log-NLMS(本文算法)進(jìn)行實(shí)驗(yàn)，得到的AAE和ARE如表1所示，步行、低速騎行、高速騎行3種狀態(tài)下，經(jīng)過本文算法去噪后，測試者的AAE和ARE整體上都比對比算法要低，尤其是高速騎行狀態(tài)下，本文算法5名測試者的AAE均值2.62 BPM要遠(yuǎn)低于對比算法AAE均值5.66 BPM。結(jié)果表明修正后的加速度作為自適應(yīng)濾波的參考信號比原始加速度作為自適應(yīng)濾波的參考信號去噪效果要好，心率估計(jì)更加準(zhǔn)確。如Bland-Altman圖(圖6)所示，大部分?jǐn)?shù)據(jù)點(diǎn)都在劃定的95%置信區(qū)內(nèi)，經(jīng)計(jì)算，本文算法計(jì)算的心率與真實(shí)心率的誤差有93.33%的概率在置信區(qū)內(nèi)，表明本文算法去運(yùn)動(dòng)偽影效果以及心率估計(jì)都有較高的準(zhǔn)確性。

圖6 真實(shí)心率與本文估計(jì)心率的Bland-Altman圖

表1 步行、低速騎行和高速騎行狀態(tài)下兩種算法心率提取的AAE、ARE

3 結(jié)論

為消除PPG信號中運(yùn)動(dòng)偽影的影響，提高心率估計(jì)的準(zhǔn)確性，設(shè)計(jì)了一種基于修正加速度變步長自適應(yīng)濾波的心率估計(jì)算法，得到以下結(jié)論。

(1)由于加速度計(jì)測量是基于慣性參考系，所以測量的加速度包含重力分量，為了得到外部加度，采用卡爾曼濾波對加速度進(jìn)行修正，消除重力分量的影響，相較于以原始加速度信號為參考信號的心率估計(jì)算法，去除MA能力效果更好，心率估計(jì)更加準(zhǔn)確。

(2)穿戴設(shè)備要求的計(jì)算準(zhǔn)確、功耗低，大部分抗MA算法都未能解決復(fù)雜度和精確度兩者之間的矛盾。Log-NLMS自適應(yīng)濾波算法簡單，收斂速度塊，誤差穩(wěn)態(tài)高，符合可穿戴設(shè)備的需求。