基于VMD-NLMS的運(yùn)動(dòng)狀態(tài)血氧監(jiān)測(cè)算法

張 林,余成波,譚 擁,米有玉

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

1 引言

血氧飽和度(SpO2)是一項(xiàng)表征人體心肺功能及循環(huán)系統(tǒng)健康程度的重要生理參數(shù)[1]。近年來(lái)推出的可穿戴的脈搏血氧監(jiān)測(cè)設(shè)備,可以通過(guò)采集PPG信號(hào)實(shí)時(shí)的監(jiān)測(cè)用戶的血氧飽和度,有效地避免一些突發(fā)的心血管疾病。但可穿戴設(shè)備在為用戶帶來(lái)便捷體驗(yàn)的同時(shí),也引入了更大的運(yùn)動(dòng)干擾,而常見的血氧監(jiān)測(cè)算法的魯棒性較差,并不能對(duì)抗較強(qiáng)的運(yùn)動(dòng)噪聲,使測(cè)量結(jié)果往往存在較大誤差。因此,研究一種魯棒性更強(qiáng)的血氧監(jiān)測(cè)算法有很大的意義和價(jià)值。

目前去除PPG信號(hào)中的運(yùn)動(dòng)偽影的方法主要有自適應(yīng)濾波和經(jīng)驗(yàn)?zāi)B(tài)分解及其改進(jìn)算法等。有陳真誠(chéng)等人[2-4]利用EMD處理PPG信號(hào)時(shí),由于其理論存在缺陷,分解過(guò)程易出現(xiàn)模態(tài)混疊、端點(diǎn)效應(yīng)等問(wèn)題[5],去噪效果并不理想。因此,有Kang等人[6-7]采用VMD處理PPG信號(hào),克服了模態(tài)混疊和邊界效應(yīng)的影響,但重構(gòu)時(shí)僅使用加速度信號(hào)作為參考選擇IMF,使得處理后的PPG信號(hào)仍有噪聲殘留。此外,自適應(yīng)濾波方法也常用于PPG信號(hào)處理。Han等人[8]將同步測(cè)量的加速度信號(hào)作為自適應(yīng)濾波器的參考信號(hào)來(lái)濾除運(yùn)動(dòng)偽影,對(duì)比EMD降低了時(shí)間復(fù)雜度,但其去噪效果過(guò)于依賴參考信號(hào)的質(zhì)量,因此常與其它方法結(jié)合使用[9-11],但結(jié)合算法仍使用加速度信號(hào)作為噪聲參考,不僅增加了硬件成本,且構(gòu)造的運(yùn)動(dòng)噪聲參考信號(hào)較為片面,去噪效果雖有提高,但并不是十分理想。

為了有效地去除PPG信號(hào)中的運(yùn)動(dòng)偽影,本文提出了一種基于變分模態(tài)分解和歸一化變步長(zhǎng)的自適應(yīng)濾波算法來(lái)消除運(yùn)動(dòng)干擾。首先通過(guò)VMD算法將預(yù)處理后的PPG信號(hào)分解為若干個(gè)本征模態(tài)分量,然后通過(guò)計(jì)算每個(gè)模態(tài)的多尺度排列熵值,選擇合適的模態(tài)重構(gòu)噪聲參考信號(hào),最后使用歸一化的最小均方誤差自適應(yīng)算法濾除PPG信號(hào)中的運(yùn)動(dòng)偽影。

2 相關(guān)算法原理

2.1 VMD算法

變分模態(tài)分解[12]是Konstantin等人在2014年提出的一種自適應(yīng)的信號(hào)處理方法,通過(guò)迭代搜尋變分模態(tài)的最優(yōu)解,不斷更新各模態(tài)函數(shù)及中心頻率,得到若干具有一定帶寬的模態(tài)函數(shù)。變分模態(tài)分解是EMD算法的改進(jìn),其優(yōu)點(diǎn)是更好的避免了模態(tài)混疊和端點(diǎn)效應(yīng)等問(wèn)題。

VMD方法將信號(hào)分解成預(yù)先定義的K個(gè)本證模態(tài)分量IMF,這里的IMF即 被重新定義為一個(gè)AM-FM信號(hào)

uk(t)=Ak(t)cos(φk(t))

(1)

其中Ak(t)為瞬時(shí)幅值,Ak(t)≥0;φk(t)為相位;φk(t)對(duì)t求一階微分,得到uk(t)的瞬時(shí)頻率ωk(t)=φk(t)′。

假設(shè)每個(gè)IMF在信號(hào)中心頻率附近具有有限的帶寬,由此構(gòu)造了如下的約束變分問(wèn)題

(2)

為了求解式(2),引入二次懲罰項(xiàng)α和Lagrange乘數(shù)λ(t),把問(wèn)題轉(zhuǎn)換為了無(wú)約束的變分問(wèn)題

L({uk},{ωk},λ)

(3)

利用交替方向乘子(Alternating Direction Method of Multipliers,ADMM)對(duì)式(3)求解得到的K個(gè)IMF的頻域結(jié)果進(jìn)行反傅里葉變換,得到其時(shí)域信號(hào)。

2.2 NLMS算法原理

歸一化最小均方誤差(NLMS)算法采用變步長(zhǎng)的方法來(lái)縮短自適應(yīng)收斂過(guò)程,解決了最小均方誤差(Least Mean Square,LMS)算法在收斂速度和穩(wěn)態(tài)誤差對(duì)步長(zhǎng)要求相互矛盾的問(wèn)題[13],是一種改進(jìn)的LMS算法。NLMS算法具體原理如下

(4)

e(n)=d(n)-y(n)=d(n)-wT(n)x(n)

(5)

(6)

x(n)輸入信號(hào),d(n)為參考信號(hào),w(n)為濾波器抽頭系數(shù),e(n)為誤差信號(hào)。μ為步長(zhǎng)因子,為了保證自適應(yīng)濾波器的性能,取值范圍應(yīng)滿足0≤μ≤2;γ是為了避免xT(n)x(n)過(guò)小導(dǎo)致步長(zhǎng)太大而預(yù)先設(shè)置的矯正值,一般取0≤γ≤1。

3 VMD-NLMS改進(jìn)算法

本文提出的改進(jìn)算法主要分為兩個(gè)部分:基于VMD的運(yùn)動(dòng)噪聲參考信號(hào)重構(gòu)和基于NLMS的PPG信號(hào)的運(yùn)動(dòng)偽影消除。

3.1 基于VMD的參考噪聲重構(gòu)

用VMD算法對(duì)預(yù)處理后的PPG信號(hào)進(jìn)行分解,得到K個(gè)中心頻率由低到高的模態(tài)分量,而利用VMD方法進(jìn)行參考噪聲重構(gòu)的算法思想就是從K個(gè)本征模態(tài)分量中選出包含運(yùn)動(dòng)偽影的模態(tài),然后進(jìn)行重構(gòu)。因此,重構(gòu)的首要問(wèn)題就是建立相關(guān)的準(zhǔn)則判斷并選取出運(yùn)動(dòng)偽影對(duì)應(yīng)的模態(tài)分量。本文通過(guò)計(jì)算IMF的多尺度排列熵MPE來(lái)選擇重構(gòu)信號(hào)的IMF。

多尺度排列熵是一種衡量時(shí)間序列復(fù)雜度的非線性動(dòng)力學(xué)特征指標(biāo)[14-15],是在排列熵基礎(chǔ)上的改進(jìn)和優(yōu)化,通過(guò)對(duì)時(shí)間序列的粗?；?得到不同時(shí)間尺度的熵值信息,文獻(xiàn)[16]中描述了其原理及計(jì)算方法。而對(duì)IMF分量來(lái)說(shuō),計(jì)算得到的MPE值越小,說(shuō)明該尺度下的時(shí)間序列越規(guī)則,越接近PPG信號(hào)的原始成分,因?yàn)闊o(wú)噪的PPG信號(hào)為一個(gè)周期信號(hào);MPE值越大,時(shí)間序列越隨機(jī),越接近噪聲成分。計(jì)算得到各IMF分量的MPE值后,選擇合適的尺度因子得到判斷IMF是否為噪聲分量的閾值,然后根據(jù)設(shè)定的閾值選取合適的IMF分量進(jìn)行重構(gòu)即可得到噪聲參考信號(hào)。

3.2 基于NLMS的運(yùn)動(dòng)偽影消除

將重構(gòu)的噪聲參考信號(hào)和原PPG信號(hào)共同作為NLMS自適應(yīng)濾波器的輸入,得到去除運(yùn)動(dòng)偽影后的PPG信號(hào)。整個(gè)算法流程圖如圖1所示。

圖1 VMD-NLMS算法流程圖

4 數(shù)據(jù)的獲取與預(yù)處理

4.1 PPG信號(hào)采集設(shè)備的搭建

本文利用SFH7050傳感器、MPU6050六軸加速度傳感器、AFE4404模擬前端、主控芯片STM32F103及藍(lán)牙傳輸模塊等模塊搭建了PPG信號(hào)的采集設(shè)備如圖2所示。其中SFH7050傳感器搭載了血氧檢測(cè)所需的波長(zhǎng)為940 nm的紅外光LED和波長(zhǎng)為660 nm的紅光LED[17],而模擬前端AFE4404高度集成了LED驅(qū)動(dòng)電路、放大電路、濾波電路、采樣電路等血氧前端采集必需的功能模塊。

圖2 PPG信號(hào)采集設(shè)備

為了探究不同運(yùn)動(dòng)狀態(tài)下的PPG信號(hào),本文用設(shè)計(jì)的PPG信號(hào)采集裝置采集了常見的三種運(yùn)動(dòng)狀態(tài):步行、跑步和騎行時(shí)腕部所產(chǎn)生的PPG信號(hào),且在采集過(guò)程中,盡量保證手腕自然擺動(dòng),采樣的頻率設(shè)置為256 Hz。如圖3為采集的行走狀態(tài)下的雙路PPG信號(hào)。

圖3 步行狀態(tài)下的雙路PPG信號(hào)

由圖3可知,采集到的未經(jīng)處理的PPG信號(hào)中噪聲較多。因此,若使用未處理的PPG信號(hào)進(jìn)行血氧飽和度值的計(jì)算,結(jié)果與真實(shí)值會(huì)出現(xiàn)很大的誤差,甚至出現(xiàn)錯(cuò)誤的結(jié)果。

4.2 PPG信號(hào)預(yù)處理

由于人體每分鐘的心率正常為40～200 次/分鐘,對(duì)應(yīng)的頻率范圍大約為0.7～3.5 Hz。而運(yùn)動(dòng)偽影也在這個(gè)頻率范圍附近。為了充分保留有用的信息并去掉其它一些不感興趣的噪聲如基線漂移和高頻噪聲等,首先對(duì)PPG信號(hào)進(jìn)行預(yù)處理,采用長(zhǎng)度為8 秒,步長(zhǎng)為2 秒的滑動(dòng)時(shí)間窗口截取PPG信號(hào);然后使用截止頻率為0.7～3.5 Hz的帶通濾波器對(duì)信號(hào)進(jìn)行處理。由于巴特沃斯帶通濾波器在其通頻帶內(nèi)有平穩(wěn)的幅頻特性,濾波器本文選擇了巴特沃斯帶通濾波器,設(shè)置通帶頻率為0.7 Hz～3.5 Hz,阻帶截止頻率為0.5 Hz～5.0 Hz,通帶最大衰減系數(shù)為3 dB,阻帶最小衰減系數(shù)為18 dB,PPG信號(hào)的預(yù)處理結(jié)果如圖4所示。

圖4 預(yù)處理前后的PPG信號(hào)

由圖4可以看到,預(yù)處理前的PPG信號(hào)的上面有呼吸諧波等因素造成的表現(xiàn)為基線漂移的低頻噪聲的干擾,有工頻干擾等其它一些因素造成的表現(xiàn)為“毛刺”狀的高頻噪聲的干擾,還有由于采集時(shí)的運(yùn)動(dòng)干擾帶來(lái)的表現(xiàn)為幅值的大幅增大或減小的運(yùn)動(dòng)噪聲的干擾。而經(jīng)過(guò)預(yù)處理后的PPG信號(hào),從時(shí)域特征來(lái)看,高頻噪聲帶來(lái)的信號(hào)的“毛刺”得到了很好的抑制,基線的漂移也有了一定程度的抑制,但運(yùn)動(dòng)噪聲并未得到有效地抑制。因此,預(yù)處理后的PPG信號(hào)需要進(jìn)一步的去除信號(hào)中的噪聲,才能進(jìn)行血氧飽和度的計(jì)算。

5 實(shí)驗(yàn)結(jié)果及分析

5.1 去噪結(jié)果及對(duì)比分析

為了進(jìn)一步的去除PPG信號(hào)中的運(yùn)動(dòng)噪聲,本文對(duì)預(yù)處理后的PPG信號(hào)進(jìn)行VMD分解,而VMD算法的分解效果的好壞,主要取決于模態(tài)個(gè)數(shù)K和懲罰因子α的取值。K取值過(guò)大,會(huì)造成信號(hào)的過(guò)分解;而K取值過(guò)小,則會(huì)造成模態(tài)混疊,即不同中心頻率的信號(hào)成分經(jīng)過(guò)VMD分解后處在了同一個(gè)本征模態(tài)分量中,產(chǎn)生“虛假”分量。而α的取值主要影響VMD分解后的得到各模態(tài)的帶寬大小模態(tài)混疊程度。如果懲罰因子α的值太小,則每個(gè)頻率的帶寬太大,并且將相似頻率誤認(rèn)為同一中心頻率的大帶寬部分,從而將其分解為同一分量,出現(xiàn)模態(tài)混疊的問(wèn)題。當(dāng)該值太大時(shí),部分頻帶寬度減小,可能會(huì)導(dǎo)致信號(hào)過(guò)分解,出現(xiàn)虛假分量。K的取值可由待處理信號(hào)的頻譜中峰值的個(gè)數(shù)來(lái)確定,因此根據(jù)預(yù)處理后信號(hào)的頻譜中峰的個(gè)數(shù),本文中K設(shè)置為7。預(yù)處理后的信號(hào)經(jīng)過(guò)變分模態(tài)分解后得到了IMF1-IMF7共7個(gè)本證模態(tài)分量如圖5所示,各個(gè)本征模態(tài)分量的中心頻率與處理前信號(hào)的峰值所對(duì)應(yīng)的頻率分別如表1所示,對(duì)比分解前后各模態(tài)的中心頻率(cf),與分解前信號(hào)的主要頻率(mf),發(fā)現(xiàn)頻率值基本對(duì)應(yīng),可見VMD的結(jié)果可以較好反應(yīng)原信號(hào)的頻率特征。

表1 PPG信號(hào)的主要頻率與IMF的中心頻率

圖5 PPG信號(hào)VMD分解得到的IMF分量

懲罰因子α的取值可通過(guò)信號(hào)分解后各模態(tài)在頻域的混疊程度來(lái)判斷。如圖6所示為K取7,α取3300時(shí),PPG信號(hào)經(jīng)分解后各個(gè)IMF的頻率-幅值圖,由圖可得,模態(tài)直接并沒(méi)有發(fā)生混疊,所以,α可按經(jīng)驗(yàn)取值為3300。

圖6 各個(gè)IMF分量的頻率-幅值圖

PPG信號(hào)經(jīng)過(guò)VMD分解后得到7個(gè)IMF。本文通過(guò)計(jì)算各IMF分量的MPE來(lái)確定各IMF分量所代表的頻率成分。計(jì)算MPE時(shí),嵌入維數(shù)m、時(shí)間延遲τ、尺度因子s這些參數(shù)的選擇會(huì)對(duì)熵值大小造成影響。

時(shí)間延遲τ的選擇對(duì)熵值的計(jì)算影響較小,本文計(jì)算了τ為1～5時(shí),單個(gè)模態(tài)的排列熵值如圖7所示?？梢钥吹?τ的變化對(duì)熵值的大小的影響較小,因此本文中,時(shí)間延遲τ設(shè)置為1。嵌入維數(shù)m的取值一般為3～7[18]。m過(guò)大時(shí),信號(hào)在相空間中重構(gòu)時(shí)會(huì)花費(fèi)更多的時(shí)間,計(jì)算時(shí)間會(huì)因此而大大增加;m過(guò)小又會(huì)使信號(hào)在相空間重構(gòu)時(shí)信息太少而丟失信號(hào)的細(xì)節(jié)特征。本文統(tǒng)計(jì)了當(dāng)時(shí)間延遲τ取1,m取3～7時(shí)計(jì)算單個(gè)排列熵所用的時(shí)間如表2?？紤]到計(jì)算效率和特征損失的平衡,本文的嵌入維數(shù)設(shè)置為6。尺度因子s的值對(duì)熵值的計(jì)算有較大的影響。若s取值過(guò)小,則不能最大限度的提取信號(hào)的特征信息;若s取值過(guò)大,可能抹除信號(hào)之間的復(fù)雜度差異,還會(huì)導(dǎo)致熵值計(jì)算不穩(wěn)定。本文計(jì)算了s從1到20時(shí),7個(gè)IMF的熵值如圖8所示,當(dāng)s大于7后各個(gè)IMF的排列熵值開始重合,無(wú)法反映出IMF的表征特征;當(dāng)s為4～6時(shí),各IMF分量的排列熵值區(qū)別較為明顯,其熵值如表3所示。但為了在保證熵值穩(wěn)定性的同時(shí),最大限度的提取信號(hào)的特征信息,將本文的尺度因子s設(shè)置為6。

表2 不同嵌入維數(shù)下排列熵的計(jì)算時(shí)間

表3 IMF在尺度因子不同時(shí)的多尺度排列熵

圖7 單個(gè)IMF在不同時(shí)間延遲下的排列熵

圖8 各IMF不同尺度下的排列熵值

本文通過(guò)驗(yàn)證最終的去噪效果,來(lái)確定MPE的最佳閾值。經(jīng)50次試驗(yàn)后,用來(lái)構(gòu)造噪聲參考信號(hào)的IMF排列熵閾值范圍分布在(0,0.2320)∪(0.3450,+∞)的次數(shù)為43次,因此將閾值設(shè)置為σmin≤0.2320和σmax≥0.3450,即當(dāng)MPE≤0.2320 時(shí),為低頻噪聲成分,當(dāng)MPE≥0.3450 時(shí),為高頻噪聲成分,此時(shí)重構(gòu)的噪聲參考信號(hào)最好。將重構(gòu)的參考信號(hào)和預(yù)處理后的PPG信號(hào)共同輸入NLMS自適應(yīng)濾波器中去噪得到的干凈的PPG信號(hào)。

此外,本文選擇了用加速度信號(hào)作為參考選擇IMF進(jìn)行PPG信號(hào)重構(gòu)的EMD算法(A-EMD)和VMD算法(A-VMD)對(duì)PPG信號(hào)進(jìn)行處理,結(jié)果如圖9所示。并計(jì)算了三種算法去噪前后的SNR和MSE如表4。

表4 三種方法處理PPG信號(hào)的SNR和MSE

圖9 三種不同方法處理PPG信號(hào)的結(jié)果

由圖9知,使用A-EMD方法處理PPG信號(hào)后,一些較強(qiáng)的運(yùn)動(dòng)偽影得到了濾除,但由于發(fā)生了模態(tài)混疊,導(dǎo)致在信號(hào)重構(gòu)時(shí)丟失了部分有用的信息,在4 s～8 s出現(xiàn)較大的波形失真;使用A-VMD處理PPG信號(hào)后,運(yùn)動(dòng)偽影得到了較好的抑制,但信號(hào)丟失了一些峰值特征,且運(yùn)動(dòng)偽影較強(qiáng)時(shí),抑制效果不夠明顯;使用本文方法處理PPG信號(hào)后,不僅有效的抑制了運(yùn)動(dòng)偽影,而且還保留了信號(hào)在一些細(xì)節(jié)處的特征。由表4可知,本文方法對(duì)比另外兩種算法SNR提高了1 dB～2 dB,MSE也更小,因此去噪效果更好。圖10為不同運(yùn)動(dòng)狀態(tài)下通過(guò)本文方法對(duì)PPG信號(hào)的處理結(jié)果。

圖10 VMD-NLMS處理前后的PPG信號(hào)和生成的參考信號(hào)

圖10中,紅線為紅外光PPG信號(hào),黑色線為紅光PPG信號(hào),可以看出,經(jīng)過(guò)本文算法處理后,三種運(yùn)動(dòng)狀態(tài)的PPG信號(hào)在保持脈搏信號(hào)時(shí)域特征的同時(shí),較好的抑制了運(yùn)動(dòng)偽影。

5.2 血氧飽和度計(jì)算結(jié)果及對(duì)比分析

本文通過(guò)提取去噪后雙路PPG信號(hào)的直流分量和交流分量計(jì)算血氧飽和度值。實(shí)驗(yàn)對(duì)象為5名身體健康的人,分別采集1分鐘內(nèi)的步行、跑步、騎行三種狀態(tài)下左手腕部的PPG信號(hào),然后使用三種算法處理PPG信號(hào)并計(jì)算血氧飽和度值,同時(shí)選取標(biāo)準(zhǔn)的脈搏血氧儀同步檢測(cè)(同步采集時(shí)右手盡可能保持靜止)的血氧飽和度值作為標(biāo)準(zhǔn)值,結(jié)果如表5-7所示。由此可計(jì)算各方法測(cè)量血氧值的平均絕對(duì)誤差值(Mean Absolute Deviation,MAE)和平均絕對(duì)誤差百分比(Mean Absolute Error Percentage, MAEP)如表8所示。

表5 步行狀態(tài)下的血氧飽和度計(jì)算值

表6 跑步狀態(tài)下的血氧飽和度計(jì)算值

表7 騎行狀態(tài)下的血氧飽和度計(jì)算值

表8 血氧飽和度值測(cè)量結(jié)果的誤差….(%)

由表8可知相比與A-EMD和A-VMD兩種方法,本文提出的算法平均絕對(duì)誤差和平均絕對(duì)誤差百分比減小了2.0%～3.0%,說(shuō)明使用本文算法得到的血氧值更接近于真實(shí)值。

6 結(jié)論

本文提出了一種VMD-NLMS的方法來(lái)去除PPG信號(hào)中的運(yùn)動(dòng)偽影。通過(guò)VMD方法構(gòu)造了噪聲參考信號(hào),作為自適應(yīng)濾波器的輸入來(lái)去除PPG信號(hào)中的運(yùn)動(dòng)偽影。該方法的優(yōu)點(diǎn)在于自適應(yīng)濾波器的參考信號(hào)是利用VMD分解從PPG信號(hào)內(nèi)部產(chǎn)生的而不是使用三軸加速度信號(hào)作為噪聲參考,不僅節(jié)省了硬件成本,而且重構(gòu)的噪聲參考信號(hào)比用單一加速度作為參考信號(hào)更有運(yùn)動(dòng)噪聲代表性。通過(guò)對(duì)常見的3種運(yùn)動(dòng)狀態(tài)下PPG信號(hào)的分析和處理,計(jì)算相應(yīng)的血氧飽和度,證明了本算法更強(qiáng)的魯棒性和準(zhǔn)確性;通過(guò)與基于加速度的EMD算法和基于加速度的VMD算法去噪結(jié)果的對(duì)比,證明了本算法的優(yōu)越性。因此,本文提出的算法可以很好的應(yīng)用于脈搏血氧值的計(jì)算。