文 武，孫再敏

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

引 言

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

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

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

1 理論研究

1.1 奇異譜分析

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

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

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

(1)

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

(2)

(3)分組?；揪仃嘪i被分割成m個不相交子集I1,…,Im。將對應(yīng)于I子集的矩陣XI定義為XI=XI(1)+…XI(p)。軌跡矩陣X現(xiàn)在可以被寫成：

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

(3)

(4)對角線平均。每個矩陣XI被轉(zhuǎn)換成一個新的長度N的時間序列。設(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]可視為兩個橫向濾波器，聯(lián)合過程估計器和輔助濾波器并行工作。FTRLS算法大致可分解為前向預(yù)測、后向預(yù)測和聯(lián)合估計3個部分,在收斂過程中實現(xiàn)3個部分之間的參量的互相交換、相互作用，從而實現(xiàn)FTRLS算法[15]。下文中，e(k,N)為先驗誤差，ε(k,N)為后驗誤差，下標(biāo)f和b分別表示前向和后向。3個部分的大致過程如下。

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

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

(5)

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

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

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

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

(6)

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

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

(7)

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

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

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

(8)

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

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

(9)

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

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

(10)

2 濾波算法

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

Fig.1 Motion artifact pulse wave and normal pulse wave

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

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

Fig.2 Filtering framework

Fig.3 Three-stage adaptive filtering

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

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

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

Table 1 Duration of each data record

Fig.4 The 12s～14s record information

4 實 驗

4.1 奇異譜分析

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

Fig.5 z-axis angle correction

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

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

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

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

4.2 FTRLS

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

Fig.6 Pulse wave to remove motion artifacts

4.3 濾波效果評價

為了驗證提出的基于3級FTRLS濾波方法，將3軸加速度角度矯正后利用奇異譜方法分組重建，采用最小均方算法(least mean square,LMS)和RLS濾波，結(jié)果如圖7所示。原始信號是s1_high_resistance_bike的數(shù)據(jù)，LMS算法、RLS算法都有一定的濾波效果，但是頻率窗口在500左右時，只有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對評價濾波效果評估。 若RSNR值越大，ERMSE值越小，表明濾波效果越好[23]。均方根誤差定義為：

(11)

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

信噪比計算公式為：

RSNR=10lg(Ps/Pn)

(12)

式中，Pn為噪聲方差，Ps為信號方差。

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

Table 3 Comparison of three filtering methods

5 結(jié) 論

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