擴(kuò)展的卡爾曼濾波在近紅外光譜提取腦血流信號(hào)中的研究

劉頌陽(yáng), 劉光達(dá), 劉卓婭, 邱吉慶， 蔡 靖,朱展鵬， 張 程, 齊 遠(yuǎn)， 張 尚

1. 吉林大學(xué)， 吉林 長(zhǎng)春 130061 2. 吉林大學(xué)珠海學(xué)院， 廣東 珠海 519041 3. 吉林大學(xué)第一醫(yī)院， 吉林 長(zhǎng)春 130061

引 言

大腦作為人體的神經(jīng)中樞， 是人類(lèi)高級(jí)認(rèn)知功能的載體。 近年來(lái)， 用于研究大腦運(yùn)作機(jī)制、 疾病診斷的腦功能成像方法得到了長(zhǎng)足的發(fā)展， 腦功能成像技術(shù)主流的方法有腦電圖(electro encephalo gram, EEG)、 腦磁圖(magneto encephalo graphy, MEG)、 功能磁共振成像(functional magnetic resonance imaging， FMRI)、 功能性近紅外光譜技術(shù)(functional near infrared spectroscopy, FNIRS)等。 在上述方法中， 功能性近紅外光譜技術(shù)(FNIRS)因具有較高的空間分辨率、 適中的時(shí)間分辨率、 無(wú)創(chuàng)費(fèi)用低等優(yōu)點(diǎn)成為腦功能成像技術(shù)研究的熱點(diǎn)。 各種腦成像技術(shù)優(yōu)劣對(duì)比見(jiàn)表1。

功能性近紅外光譜技術(shù)(functional near infrared spectroscopy， fNIRS)探測(cè)腦組織活動(dòng)的方式是間接方式。 大腦的各種神經(jīng)活動(dòng)需要氧氣， 腦部血流(Celebral blood flow)做為大腦供氧的重要渠道， 它對(duì)大腦的神經(jīng)活動(dòng)是十分敏感的， fNIRS正是利用探測(cè)腦血流動(dòng)力學(xué)的參數(shù)變化， 來(lái)判斷腦組織的各種神經(jīng)活動(dòng)的。 顯而易見(jiàn)的是， 大腦作為人的神經(jīng)中樞， 各種生理活動(dòng)如心跳、 呼吸、 脈搏等都會(huì)引起神經(jīng)活動(dòng)， 也自然會(huì)引起腦血流動(dòng)力學(xué)的參數(shù)變化， 因此， 如何從包含各種生理干擾的信號(hào)中提取腦組織神經(jīng)活動(dòng)信號(hào)是近紅外光譜提取腦功能信號(hào)難點(diǎn)。

表1 各種腦成像技術(shù)優(yōu)劣對(duì)比

目前， 常見(jiàn)的提取近紅外光腦信號(hào)的算法有： 脈搏色素譜法[1]、 EEMD-ICA法[3]主成分分析法(PCA)[5]、 獨(dú)立成分分析法(ICA)[6]、 相干平均法[7]、 自適應(yīng)濾波[8]等。 脈搏色素譜法采用注射顯影劑吲哚青綠的方式進(jìn)行示蹤估計(jì)； EEMD-ICA對(duì)信號(hào)進(jìn)行希爾伯特變化， 再進(jìn)行濾波； PCA和ICA這兩種方法適用于連續(xù)光譜的近紅外光譜系統(tǒng)； 自適應(yīng)濾波法可以消除心跳干擾； 相干平均法可以濾除脈搏等生理干擾。 上述算法在對(duì)近紅外腦神經(jīng)活動(dòng)信號(hào)提取時(shí)都取得了一定的進(jìn)展。

但是， 上述算法在提取近紅外光譜腦功能信號(hào)時(shí)， 重視各種生理干擾， 忽視了在提取信號(hào)中的測(cè)量干擾， 如儀器精密度、 信號(hào)傳輸中的串?dāng)_等。 研究[10-11]表明， 近紅外光譜腦信號(hào)用非線性算法提取可以去除測(cè)量干擾。 如Fritson使用非線性算法將信號(hào)Volterra[11]級(jí)數(shù)化， 大大提高了信號(hào)的信噪比。 本工作提出的擴(kuò)展的卡爾曼(Extend Klaman Filter， EKF)算法是一種非線性算法。 根據(jù)近紅外光譜提取的腦信號(hào)中的干擾成分是多維非平穩(wěn)隨機(jī)的， 利用其時(shí)變性、 功率譜不固定， 建立了數(shù)學(xué)模型， 將信號(hào)進(jìn)行泰勒級(jí)數(shù)分解， 取一階函數(shù)， 再進(jìn)行濾波和估計(jì)， 進(jìn)行了高斯(Gauss)濾波實(shí)驗(yàn)、 Valsava實(shí)驗(yàn)、 視覺(jué)誘發(fā)實(shí)驗(yàn)， 并將EKF算法和主流的EEMD算法進(jìn)行比對(duì)。

1 理論基礎(chǔ)

1.1 修正的朗伯比爾定律

朗伯比爾定律(Lambert-Beer law)是分光光度法的基本定律， 是描述物質(zhì)對(duì)某一波長(zhǎng)光吸收的強(qiáng)弱與吸光物質(zhì)的濃度及其厚度間的關(guān)系。 適用于所有的電磁輻射和所有的吸光物質(zhì)， 包括氣體、 固體、 液體、 分子、 原子和離子。 其公式為

(1)

式(1)中，ελ為吸光物質(zhì)在波長(zhǎng)λ(nm)下的消光系數(shù)，c為吸光物質(zhì)的濃度(mol·L-1)，d為光在物質(zhì)中的傳輸距離，I0和I分別代表出射、 入射光強(qiáng)，A為吸光度。

但事實(shí)上， 光在生物組織中會(huì)有散射、 吸收等現(xiàn)象， 圖1給出了光在組織傳播的路徑。 因此， 光在組織中的實(shí)際傳輸距離遠(yuǎn)大于物質(zhì)厚度。 為了消除這一誤差， 1998年， Deply D T等對(duì)朗伯比爾定律進(jìn)行了修正[12]， 修正的朗伯比爾定律見(jiàn)式(2)

(2)

圖1 光在組織中的路徑

1.2 功能性近紅外光譜

近紅外光譜是指波長(zhǎng)在650～950 nm的不可見(jiàn)光， 科學(xué)家Jobsis[2]在1977年發(fā)現(xiàn)近紅外光對(duì)人體組織具有良好的透射性， 首次驗(yàn)證了用近紅外光譜技術(shù)無(wú)創(chuàng)監(jiān)測(cè)組織血氧的可行性。 各種生理活動(dòng)都離不開(kāi)氧氣， 腦組織所需氧氣的輸送靠的是腦血流中的兩種血紅蛋白， 氧合血紅蛋白(Hb02)和還原血紅蛋白(HbR)， 當(dāng)腦血流流經(jīng)腦組織時(shí)， 血液中HbO2釋放氧氣， 供腦組織使用， HbO2轉(zhuǎn)化為HbR; 當(dāng)血液流經(jīng)肺泡時(shí)， HbR又獲得氧氣轉(zhuǎn)化為HbO2。 兩種血紅蛋白的總量基本恒定， 由上述可知， 人腦組織的生理活動(dòng)會(huì)引起HbR和HbO2兩種血紅蛋白含量[6]的變化。

人體大腦組織中， 水的含量占據(jù)了大部分的比重， 約為75%， 除此之外， 還有氧合血紅蛋白(HbO2)、 還原血紅蛋白(HbR)和黑色素等。 在650～950 nm近紅外光波段， 水的吸收率很低， 氧合血紅蛋白(HbO2)和還原血紅蛋白(HbR)的吸收率較高， 不同物質(zhì)對(duì)不同波長(zhǎng)的光的吸收度如圖2， 這一區(qū)間也稱(chēng)為光窗期。

近紅外光譜腦信號(hào)提取技術(shù)正是基于待測(cè)組織中的水、 氧合血紅蛋白(HbO2)、 還原血紅蛋白(HbR)在近紅外波段具有不同吸收光譜特性， 以及近紅外光可以穿過(guò)0.5～3 cm的外層生物組織到達(dá)腦組織這一特性進(jìn)行測(cè)定的。 若分別選取波長(zhǎng)為λ1和λ2的近紅外光源， 利用修正的朗伯比爾定律，可得入射光和出射光的光密度變化， 見(jiàn)式(3)—式(6)

圖2 組織的吸光度圖譜

(3)

(4)

進(jìn)而計(jì)算， 可得。

(5)

(6)

式中， ΔA是光密度變化，ε是某一物質(zhì)在某一波長(zhǎng)的消光系數(shù)，L是光源和光源接收器的距離， DPF和L的乘積可以近似表示光在組織中走過(guò)的實(shí)際路徑， Δc代表的是濃度變化。 進(jìn)而計(jì)算可得氧合血紅蛋白(HbO2)和還原血紅蛋白(HbR)兩種物質(zhì)的濃度變化。

2 擴(kuò)展的卡爾曼濾波算法及近紅外光譜模型的建立

2.1 卡爾曼濾波算法

卡爾曼濾波(Kalman filtering)是1961年由Rudolf Kalman提出， 是一種利用線性系統(tǒng)狀態(tài)方程， 通過(guò)系統(tǒng)輸入輸出觀測(cè)數(shù)據(jù)， 基于誤差平方和最小的原理進(jìn)行遞歸計(jì)算的方法， 其本質(zhì)是參數(shù)化的beiyesi模型， 通過(guò)對(duì)下一時(shí)刻系統(tǒng)的初步狀態(tài)估計(jì)以及測(cè)量得出的反饋相結(jié)合， 得到該時(shí)刻無(wú)限逼近真實(shí)值的狀態(tài)估計(jì)。 在現(xiàn)代隨機(jī)最優(yōu)控制和隨機(jī)信號(hào)處理技術(shù)中， 信號(hào)和噪聲往往是多維非平穩(wěn)隨機(jī)過(guò)程， 對(duì)于離散域線性系統(tǒng)， 其線性系統(tǒng)狀態(tài)預(yù)測(cè)方程式

xk=Axk-1+Bxk-1+wk-1

(7)

p(w)～N(0，Q)

(8)

cov(w)=E(wwT)=Q

(9)

式中，xk和xk-1分別為在k時(shí)刻、k-1時(shí)刻的狀態(tài)值，uk-1為在k-1時(shí)刻的控制輸入，wk-1為過(guò)程噪聲， 并且服從高斯分布；A為狀態(tài)轉(zhuǎn)移系數(shù)矩陣，B為控制輸入的增益矩陣，Q為過(guò)程激勵(lì)噪聲的協(xié)方差矩陣。

測(cè)量方程為

Zk=Hxk+vk

(10)

P(v)～(0,R)

(11)

cov(v)=E(vvT)=R

(12)

式中，Zk為k時(shí)刻的測(cè)量值，Vk為測(cè)量噪聲，H為測(cè)量矩陣，R為測(cè)量噪聲協(xié)方差矩陣。 并且滿足：E(w)=E(v)=0。

2.2 擴(kuò)展的卡爾曼濾波及模型建立

擴(kuò)展的卡爾曼濾波是在卡爾曼濾波的基礎(chǔ)上， 將非線性系統(tǒng)函數(shù)作泰勒(Taylor)級(jí)數(shù)展開(kāi)， 再進(jìn)行線性的卡爾曼濾波， 這就是擴(kuò)展的卡爾曼濾波(extend Kalman filter， EKF)。

考慮離散時(shí)間非線性動(dòng)態(tài)系統(tǒng)

(13)

其中，w(k)是過(guò)程噪聲， 協(xié)方差是Qk；Vk是測(cè)量噪聲， 其方差是Rk+1。f是狀態(tài)變量x和時(shí)間k相關(guān)的非線性函數(shù)，h是狀態(tài)方程中與狀態(tài)變量x和觀測(cè)量z相關(guān)的非線性函數(shù)。 進(jìn)行泰勒(Taylor)展開(kāi)， 取展開(kāi)后的一階線性部分作為非線性模型的近似。

狀態(tài)方程表示為

=Fk-1xk-1+wk-1+uk-1

(14)

式(14)中，f(xk|k)進(jìn)行泰勒展開(kāi)，

(15)

(16)

測(cè)量方程表示為

=Hkxk+Vk

(17)

式(17)中， Gh(xk)表示在xk|kG處進(jìn)行泰勒展開(kāi)，

(18)

參考經(jīng)典卡爾曼濾波的推導(dǎo)過(guò)程， 可以得到擴(kuò)展卡爾曼濾波的五個(gè)方程式(19)—式(23)

xk|k-1=f(xk-1|k-1)+wk-1

(19)

pk|k-1=Fpk-1|k-1FT+Q

(20)

Kk=pk|k-1HT(Hpk|k-1HT+R)-1

(21)

xk|k=xk|k-1+Kk(Zk-h(xk|k-1))

(22)

pk|k=(1-kkH)pk|k-1

(23)

式(19)是狀態(tài)預(yù)測(cè)方程， 式(20)是誤差協(xié)方差方程， 式(21)是卡爾曼增益方程， 式(22)是濾波校正方程， 式(23)是誤差協(xié)方差更新方程。

2.3 EKF模型的建立

簡(jiǎn)而言之， 人腦的神經(jīng)活動(dòng)需要氧氣， 氧氣需要腦血流中的兩種血紅蛋白進(jìn)行輸送， 必然會(huì)引起氧合血紅蛋白(HbO2)和還原血紅蛋白(HbR)濃度的變化。 由擴(kuò)展的朗伯比爾定律及近紅外光譜的特性可知， 這種變化可以由近紅外光譜測(cè)量得到， 此時(shí)測(cè)量的數(shù)據(jù)是包含了生理噪聲和測(cè)量噪聲， 將測(cè)量所得到的數(shù)據(jù)經(jīng)擴(kuò)展的卡爾曼濾波， 取得HbO2和HbR濃度的變化。 圖3表述了上述關(guān)系。

圖3 腦血流動(dòng)力學(xué)數(shù)學(xué)模型圖

3 實(shí)驗(yàn)結(jié)果與討論

3.1 擴(kuò)展的卡爾曼濾波(EKF)對(duì)Gauss噪聲濾波有效性的仿真實(shí)驗(yàn)

近紅外光譜腦信號(hào)采集過(guò)程中的噪聲干擾符合高斯分布， 取干擾信號(hào)w(k)和測(cè)量噪聲信號(hào)v(k)的幅值均為0.01的高斯白噪聲信號(hào)， 輸入信號(hào)幅值為1.0、 頻率為1.5的正弦信號(hào)。 使用EKF實(shí)現(xiàn)信號(hào)的濾波， 取Q=1,R=1。 仿真時(shí)間是5 s， 在MATLAB中編程進(jìn)行仿真， 圖4是在信號(hào)中疊加了高斯(Gauss)噪聲的波形圖， 圖5是經(jīng)EKF濾波后的信號(hào)。 從MATLAB仿真結(jié)果看， 信號(hào)平滑、 無(wú)失真、 沒(méi)有畸變， 該算法在去除高斯(Gauss)白噪聲方面效果明顯。

3.2 EKF提取HbR和HbO2濃度變化

3.2.1 實(shí)驗(yàn)裝置

為了測(cè)量整個(gè)前額葉區(qū)域的血紅蛋白和還原血紅蛋白的變化量， 設(shè)計(jì)了近紅外光譜腦信號(hào)采集裝置。 其中， 直流穩(wěn)壓電源和光源驅(qū)動(dòng)電路給近紅外光光源穩(wěn)定供電， 光源為波長(zhǎng)750和830 nm的二極管近紅外光源， 用同相放大和低通濾波做初步處理， 將處理后的信號(hào)發(fā)送至上位機(jī)， 上位機(jī)用C#編寫(xiě)， 其框圖見(jiàn)圖6。

圖4 疊加噪聲的原始信號(hào)

圖5 EKF濾波的信號(hào)

圖6 裝置設(shè)計(jì)框圖

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

為了獲得HbO2和HbR在血液中的濃度變化曲線， 使用波長(zhǎng)分別為750和830 nm的近紅外光源進(jìn)行Valsava實(shí)驗(yàn)， 近紅外光源和光源探測(cè)器用黑布完全遮起來(lái)， 被試者是25歲成年男性， 身體健康， 無(wú)煙酒嗜好。 耳朵中放置隔音棉， 盡可能減少外界環(huán)境干擾， 進(jìn)行閉氣-呼氣實(shí)驗(yàn)。 在實(shí)驗(yàn)開(kāi)始的0～60 s， 被試者正常呼吸； 第60～120 s， 被試者閉氣； 第120 s之后， 被試者恢復(fù)正常呼吸。 采集到的HbO2和HbR經(jīng)擴(kuò)展的卡爾曼濾波后變化曲線如圖7所示。

圖7 Valsava實(shí)驗(yàn)血紅蛋白變化曲線圖

3.2.3 視覺(jué)誘發(fā)實(shí)驗(yàn)與分析

為了驗(yàn)證本方法在提取氧合血紅蛋白(HbO2)和還原血紅蛋白(HbR)濃度變化的對(duì)外界刺激的敏感性， 進(jìn)行視覺(jué)誘發(fā)實(shí)驗(yàn)。 Arturs發(fā)現(xiàn)， 視覺(jué)誘發(fā)實(shí)驗(yàn)會(huì)引起腦血流參數(shù)的變化[14]， Zaletel發(fā)現(xiàn)視覺(jué)誘發(fā)實(shí)驗(yàn)和腦血流流速變化是正相關(guān)關(guān)系[15]。 本工作采用黑白相間的圖片作為視覺(jué)誘發(fā)源， 如圖8所示。

圖8 視覺(jué)誘發(fā)源

將視覺(jué)誘發(fā)源信號(hào)在電腦屏幕上交替播放， 交替頻率是5 Hz, 近紅外光源和光源探測(cè)器用黑布完全遮起來(lái)， 被試者是25歲成年男性， 身體健康， 無(wú)煙酒嗜好。 耳朵中放置隔音棉， 盡可能減少外界環(huán)境干擾。 讓被試者觀看圖片， 測(cè)得氧合血紅蛋白(HbO2)和還原血紅蛋白(HbR)濃度變化如圖9所示。

3.3 EKF和EEMD算法比對(duì)

EKF算法在提取近紅外光腦信號(hào)的有效性得到了驗(yàn)證， 進(jìn)一步對(duì)其濾波效果進(jìn)行評(píng)價(jià)， 比對(duì)方法是利用均方根誤差RMSE (root mean square error)、 相關(guān)系數(shù)r(correlation coefficient)。 RMSE可以表示測(cè)量值和真值的離散情況， 參數(shù)越小， 表明濾波效果越好。r表示的是相關(guān)系數(shù)，r越接近1，表示效果越好。 對(duì)同一被試者采集10組數(shù)據(jù)， 分別計(jì)算EEMD算法和EKF算法的RMSE值、r值， 結(jié)果如表2。 可以看出， EKF算法對(duì)比EEMD算法， 其RMSE值提高了0.96%，r值提高了0.6%。

圖9 視覺(jué)誘發(fā)實(shí)驗(yàn)血紅蛋白濃度變化曲線圖

表2 EKF和EEMD算法比對(duì)

4 結(jié) 論

設(shè)計(jì)了功能性近紅外光譜(fNIRS)腦血流信號(hào)采集裝置， 進(jìn)行了Valsava和視覺(jué)誘發(fā)實(shí)驗(yàn)， 采用擴(kuò)展的卡爾曼濾波(EKF)建立數(shù)學(xué)模型， 對(duì)所采集信號(hào)進(jìn)行了EKF處理。 通過(guò)實(shí)驗(yàn)表明： 本文所提出的裝置和算法可以有效去除符合高斯分布的干擾， 提取出腦血流中氧合血紅蛋白(HbO2)和還原血紅蛋白(HbR)變化曲線圖； 和主流的EEMD提取腦信號(hào)算法比對(duì)其RMSE值提高了0.96%，r值提高了0.6%。 表明本文所提出的方法有一定的優(yōu)越性。 并可以為相關(guān)腦部疾病診斷如癲癇病灶定位、 抑郁研究等提供了一種有效的腦信號(hào)提取途徑。