在運(yùn)動(dòng)狀態(tài)下脈搏信號(hào)的級(jí)聯(lián)濾波算法

文 武，孫再敏

(1.重慶郵電大學(xué) 通信與信息工程學(xué)院，重慶 400065；2.重慶郵電大學(xué) 通信新技術(shù)應(yīng)用研究中心，重慶 400065；3.重慶信科設(shè)計(jì)有限公司，重慶 401121)

引 言

脈搏信號(hào)是人體一種生理信號(hào)，包含著人體常見的心血管系統(tǒng)病理信息[1]。脈搏波蘊(yùn)藏著人體大量的健康狀態(tài)信息狀態(tài)，可以為診斷疾病或者人體亞健康狀態(tài)提供重要的科學(xué)依據(jù)。傳統(tǒng)的檢測(cè)儀器用粘性金屬電極和導(dǎo)電凝膠與患者身體接觸，給患者帶來不便。早在20世紀(jì)30年代，HERTZMAN就首次提出了光電容積脈搏波信號(hào)描記法[2]。光電容積脈搏波描記術(shù)(photoplethysmography，PPG)是一種眾所周知的用于監(jiān)測(cè)生理信息的非侵入性方法,可以使用光電容積脈搏描記法來監(jiān)測(cè)受測(cè)試者的健康狀況。生理信息監(jiān)測(cè)還可以作為一種運(yùn)動(dòng)監(jiān)測(cè)手段,有助于人們合理安排運(yùn)動(dòng)強(qiáng)度和運(yùn)動(dòng)量[3]。

在劇烈的運(yùn)動(dòng)狀態(tài)下，運(yùn)動(dòng)偽跡的存在使得脈搏波不能以理想的波形存在，為后續(xù)心率提取帶來困難。為了減小運(yùn)動(dòng)偽跡，參考文獻(xiàn)[4]中提出了一種新的方法使得PPG信號(hào)中的偽影(motion artifact,MA)去除結(jié)合了準(zhǔn)確的心率(heart rate,HR)頻率估計(jì)。參考文獻(xiàn)[5]中采用了一種時(shí)間序列方法，該方法通過修正個(gè)體來修正運(yùn)動(dòng)偽影PPG脈沖，所提出的方法將每個(gè)脈沖替換為全局平均脈沖和特定脈沖的凸組合局部平均脈沖而不會(huì)丟失時(shí)間信息的脈沖。參考文獻(xiàn)[6]和參考文獻(xiàn)[7]中使用小波變換去除運(yùn)動(dòng)偽跡，這種方法雖然可以將脈搏信號(hào)進(jìn)行分解，但是需要分析噪聲在哪一層，而且對(duì)于需要確定小波基及設(shè)定閾值。參考文獻(xiàn)[8]中利用遞歸最小二乘(recursive least square,RLS)和歸一化最小均方(normalizd least mean square,NLMS)自適應(yīng)濾波器用于對(duì)PPG信號(hào)進(jìn)行降噪。參考文獻(xiàn)[9]中使用奇異譜分析將脈搏波分解，將分解后分量的主頻與加速度信號(hào)的主頻進(jìn)行分析，去除加速度中噪聲頻率成分，重構(gòu)后實(shí)現(xiàn)去噪。參考文獻(xiàn)[10]中提出了一種基于多參考自適應(yīng)噪聲抵消技術(shù)的腕式PPG信號(hào)的魯棒HR估計(jì)算法，其主要挑戰(zhàn)是設(shè)計(jì)一個(gè)合格的參考噪聲信號(hào)到自適應(yīng)濾波器。參考文獻(xiàn)[11]中提出3種方法去除運(yùn)動(dòng)偽跡,包括3個(gè)部分：(1)PPG和加速度計(jì)信號(hào)的時(shí)頻譜估計(jì)；(2)通過從PPG信號(hào)中減去加速度計(jì)信號(hào)的時(shí)頻譜來去除運(yùn)動(dòng)偽影;(3)用三次樣條插值法剔除影響心率的殘余運(yùn)動(dòng)偽影。

自適應(yīng)濾波對(duì)于參考信號(hào)有要求，需要選擇合適的參考信號(hào)?；趩螀⒖夹盘?hào)，本文中加入了三軸陀螺儀信號(hào)，并且運(yùn)用奇異譜分析和快速橫向遞歸最小二乘(fast transverse recursive least square,FTRLS)算法有效去除了頻域上的運(yùn)動(dòng)干擾譜峰。

1 理論研究

1.1 奇異譜分析

奇異譜分析(singular spectrum analysis，SSA)先將1-D時(shí)域信號(hào)首先構(gòu)造成高維的軌跡矩陣，然后進(jìn)行奇異值分解(singular value decomposition,SVD)，將奇異分量歸類分組，最后重構(gòu)矩陣[12]。計(jì)算軌跡矩陣和奇異值分解稱為時(shí)間序列分解，分組和重構(gòu)步驟稱為重建[13]。奇異譜分析具體步驟如下。

(1)計(jì)算信號(hào)的軌跡矩陣X。將1維信號(hào)YN=(y1,…,yN)進(jìn)入多維的序列X1,…,XK。L是窗口長(zhǎng)度,N為信號(hào)長(zhǎng)度，2≤L≤N，K=N-L+1。軌跡矩陣X=[X1,…,XK]。

(2)計(jì)算矩陣S=XXT。將矩陣S的奇異值分解，計(jì)算矩陣S的特征值和特征向量。λ1,…,λL是S的逐漸遞減的特征值，U1,…,UL是這些特征值相應(yīng)的特征向量。

(1)

式中，d=max{i}。軌跡X可以寫成X=X1+…+Xd。則：

(2)

(3)分組。基本矩陣Xi被分割成m個(gè)不相交子集I1,…,Im。將對(duì)應(yīng)于I子集的矩陣XI定義為XI=XI(1)+…XI(p)。軌跡矩陣X現(xiàn)在可以被寫成：

X=XI(1)+…+XI(m)

(3)

(4)對(duì)角線平均。每個(gè)矩陣XI被轉(zhuǎn)換成一個(gè)新的長(zhǎng)度N的時(shí)間序列。設(shè)X是元素xij的L×K矩陣，1≤i≤L，1≤j≤K，L′=min(L,K)，K′=max(L,K)，N=L+K-1。若L

(4)

1.2 FTRLS

FTRLS算法[14]可視為兩個(gè)橫向?yàn)V波器，聯(lián)合過程估計(jì)器和輔助濾波器并行工作。FTRLS算法大致可分解為前向預(yù)測(cè)、后向預(yù)測(cè)和聯(lián)合估計(jì)3個(gè)部分,在收斂過程中實(shí)現(xiàn)3個(gè)部分之間的參量的互相交換、相互作用，從而實(shí)現(xiàn)FTRLS算法[15]。下文中，e(k,N)為先驗(yàn)誤差，ε(k,N)為后驗(yàn)誤差，下標(biāo)f和b分別表示前向和后向。3個(gè)部分的大致過程如下。

(1)前向預(yù)測(cè)。在前向預(yù)測(cè)關(guān)系式中，瞬時(shí)前向后驗(yàn)誤差為：

εf(k,N)=x(k)-wfT(k,N)x(k-1,N)=

(5)

式中,x(k)為PPG信號(hào)，wf(k,N)為濾波器在k時(shí)刻的前向預(yù)測(cè)權(quán)重向量。

前向預(yù)測(cè)權(quán)重向量[16]更新方程為:

wf(k,N)=wf(k-1,N)+

φ(k-1,N)εf(k,N)

(6)

式中,φ(k-1,N)為中間變量。

(2)后向預(yù)測(cè)。在后向預(yù)測(cè)關(guān)系式中，后向的后驗(yàn)和先驗(yàn)預(yù)測(cè)誤差關(guān)系表達(dá)式為：

(7)

式中，γ(k-1,N)為k-1時(shí)刻后驗(yàn)誤差與前驗(yàn)誤差的轉(zhuǎn)換因子。后向預(yù)測(cè)權(quán)重向量更新方程[17]為：

wb(k,N)=wb(k-1,N)+

φ(k-1,N)εb(k,N)

(8)

(3)聯(lián)合估計(jì)。其先驗(yàn)誤差可寫成：

e(k,N)=d(k)-wT(k-1,N)x(k,N)

(9)

式中，d(k)為期望信號(hào)。聯(lián)合過程估計(jì)濾波器權(quán)重向量為：

w(k,N)=w(k-1,N)+φ(k,N)ε(k,N)

(10)

2 濾波算法

正常脈搏波形和存在運(yùn)動(dòng)偽跡干擾脈搏波波形如圖1所示。在自行車上高速騎行時(shí)，運(yùn)動(dòng)偽跡已經(jīng)明顯嚴(yán)重破壞了脈搏波波形，在強(qiáng)烈運(yùn)動(dòng)條件下脈搏波已不是正常波形，所以去除運(yùn)動(dòng)偽跡顯得尤為重要。

Fig.1 Motion artifact pulse wave and normal pulse wave

由于加速度計(jì)適用于長(zhǎng)時(shí)間追蹤角度變化，加速度計(jì)計(jì)算角度沒有累積誤差，三軸陀螺儀由角速度積分得到的角度由于溫漂、單次迭代等原因,往往偏差較大[18]。所以用互補(bǔ)濾波法使得它們的優(yōu)勢(shì)融合起來，進(jìn)行角度矯正[19]。

FTRLS算法是通過快速橫向?yàn)V波算法代替?zhèn)鹘y(tǒng)的RLS算法[20]。比起傳統(tǒng)的RLS算法,FTRLS降低了計(jì)算復(fù)雜度，減少了算法收斂時(shí)所花費(fèi)的時(shí)間。本文中采用的FTRLS算法框架如圖2所示。利用PPG信號(hào)中的運(yùn)動(dòng)偽跡分量具有同時(shí)加速度信號(hào)的頻率，首先將校正后的三軸加速度計(jì)通過SSA，經(jīng)過嵌入生成矩陣、SVD、歸類分組和重構(gòu)，分組時(shí)將相同頻率成分的分量分為一組；其次，將奇異譜分析分組后的信號(hào)為參考信號(hào)，通過快速橫向?yàn)V波RLS濾波器，將構(gòu)造出的新分量依次輸入到3級(jí)的FTRLS自適應(yīng)濾波器中,每個(gè)FTRLS自適應(yīng)濾波器接收前級(jí)產(chǎn)生的剩余信號(hào)和本級(jí)對(duì)應(yīng)的輸入的參考信號(hào)分量。如圖3所示，PPG信號(hào)和x軸加速度參考信號(hào)sx輸入到FTRLS自適應(yīng)濾波器，輸出信號(hào)作為2級(jí)輸入和y軸加速度信號(hào)sy進(jìn)行自適應(yīng)濾波，同理輸出作為3級(jí)輸入與z軸加速度信號(hào)sz自適應(yīng)濾波。

Fig.2 Filtering framework

Fig.3 Three-stage adaptive filtering

3 數(shù)據(jù)庫(kù)

采用mimic數(shù)據(jù)庫(kù)，采集了8名參與者(3名男性和5名女性)的腕部PPG數(shù)據(jù)，年齡為22歲～32歲(平均26.5歲)。參與者進(jìn)行了一種或多種不同類型的鍛煉，數(shù)據(jù)庫(kù)是在加速度計(jì)和陀螺儀運(yùn)動(dòng)測(cè)量?jī)煞N情況下進(jìn)行體育鍛煉時(shí)記錄的信號(hào)。此數(shù)據(jù)庫(kù)與以前的公共數(shù)據(jù)庫(kù)不同，它包括在胸部同時(shí)采集的心電圖(electrocardiography,ECG)，三軸陀螺儀信號(hào)、±2g和±16g三軸加速度信號(hào)、以及ECG采樣時(shí)間和R峰時(shí)間，因?yàn)樗ㄔ谛凶?、跑步、輕松騎自行車和高阻力下騎自行車[21]。例如在高阻力下騎自行車時(shí)采集的12s到14s脈搏數(shù)據(jù)，如圖4所示。

表1表示每個(gè)數(shù)據(jù)記錄的持續(xù)時(shí)間?！?”表示數(shù)據(jù)庫(kù)中參與者沒有參與活動(dòng)。

Table 1 Duration of each data record

Fig.4 The 12s～14s record information

4 實(shí) 驗(yàn)

4.1 奇異譜分析

本文中采用互補(bǔ)濾波方法，利用三軸加速度計(jì)和三軸陀螺儀得到真實(shí)角度，利用陀螺儀角度來矯正加速度信號(hào)。三軸速度計(jì)需要濾掉高頻信號(hào)，三軸陀螺儀需要濾掉低頻信號(hào)，所以互補(bǔ)系數(shù)正好相加為1，調(diào)整它們所占的比重進(jìn)行濾波。如圖5所示,z軸角度矯正，加速度波形受到的干擾被濾掉，陀螺儀低頻信號(hào)也被濾去。

Fig.5 z-axis angle correction

矯正后的加速度信號(hào)，利用SSA算法分組為不同頻率成分的信號(hào)重構(gòu)。SSA算法可以對(duì)窗口長(zhǎng)度和分組兩個(gè)參量進(jìn)行優(yōu)化。對(duì)于窗口長(zhǎng)度，最好選擇小于將要分解信號(hào)長(zhǎng)度一半的值，可以平衡誤差。在測(cè)試中，選擇的窗口長(zhǎng)度為350,重建效果最好。將加速度數(shù)據(jù)矯正后經(jīng)過SSA分析，得到不同頻率的信號(hào)以及主要成分。

SSA算法相當(dāng)于壓縮重構(gòu)，去除冗余信息，減小信息大量存儲(chǔ)問題，提高后續(xù)計(jì)算速度問題[22]?？紤]不同特征值進(jìn)行重建時(shí)，原始信號(hào)與重構(gòu)信號(hào)互相關(guān)結(jié)果如表2所示。當(dāng)信息被丟失過多時(shí)，信號(hào)將無法得

Table 2 Cross-correlation results of original signal and reconstructed signal

到正確的重構(gòu)，結(jié)果表明，只考慮前8個(gè)特征三元組就可以在不損失太多信息情況下重建信號(hào)。

4.2 FTRLS

采用信號(hào)為1000采樣點(diǎn)，采樣率為256。將分組以后的加速度信號(hào)作為FTRLS的參考信號(hào)。結(jié)果如圖6所示,自行車高速騎行下s1_high_resistance_bike的原始信號(hào)和FTRLS濾波去除運(yùn)動(dòng)偽跡以后的信號(hào)，以及對(duì)應(yīng)的功率譜密度。FIR濾波器階數(shù)為10，定義遺忘因子為1。去噪前和去噪后的功率譜對(duì)比，可以看到在FTRLS濾波以后頻率窗口在500左右的第二大波峰去除掉，從而運(yùn)動(dòng)干擾帶來的波峰去除，為后續(xù)心率計(jì)算帶來便利，從濾波后的脈搏波的圖中可以看到還保留了重博波信息。

Fig.6 Pulse wave to remove motion artifacts

4.3 濾波效果評(píng)價(jià)

為了驗(yàn)證提出的基于3級(jí)FTRLS濾波方法，將3軸加速度角度矯正后利用奇異譜方法分組重建，采用最小均方算法(least mean square,LMS)和RLS濾波，結(jié)果如圖7所示。原始信號(hào)是s1_high_resistance_bike的數(shù)據(jù)，LMS算法、RLS算法都有一定的濾波效果，但是頻率窗口在500左右時(shí)，只有FTRLS濾除了功率譜中第二大峰。

Fig.7 Three methods of s1_high_resistance_bike to remove motion artifacts

本文中采用均方根誤差(root mean square error，RMSE)ERMSE和信噪比(signal-to-noise ratio，SNR)RSNR對(duì)評(píng)價(jià)濾波效果評(píng)估。 若RSNR值越大，ERMSE值越小，表明濾波效果越好[23]。均方根誤差定義為：

(11)

式中，s(k)為信號(hào)，y(k)為濾波后信號(hào)。

信噪比計(jì)算公式為：

RSNR=10lg(Ps/Pn)

(12)

式中，Pn為噪聲方差，Ps為信號(hào)方差。

從表3中可知，本文中提出的濾波算法去除運(yùn)動(dòng)偽跡比LMS,RLS的信噪比分別提升了12.2%,6.7%。同時(shí)，均方根誤差分別減少了30%,11%。相比LMS,RLS方法，本文中處理脈搏波的速度更快，完全滿足脈搏波在實(shí)時(shí)處理的需求。

Table 3 Comparison of three filtering methods

5 結(jié) 論

本文中提出的方法打破了自適應(yīng)濾波單一的三軸加速度參考信號(hào)，采用了陀螺儀信號(hào)作為參考信號(hào)，結(jié)合了奇異譜分析和FTRLS算法，能有效去除噪聲,還能保留重博波信息。仿真實(shí)驗(yàn)驗(yàn)證了這種去除運(yùn)動(dòng)偽跡算法的有效性，在高速自行車騎行下脈搏波中去除運(yùn)動(dòng)偽跡取得了很好的效果。