劉頌陽, 劉光達, 劉卓婭, 邱吉慶， 蔡 靖,朱展鵬， 張 程, 齊 遠， 張 尚

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

引 言

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

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

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

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

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

1 理論基礎(chǔ)

1.1 修正的朗伯比爾定律

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

(1)

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

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

(2)

圖1 光在組織中的路徑

1.2 功能性近紅外光譜

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

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

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

圖2 組織的吸光度圖譜

(3)

(4)

進而計算， 可得。

(5)

(6)

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

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

2.1 卡爾曼濾波算法

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

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

(7)

p(w)～N(0，Q)

(8)

cov(w)=E(wwT)=Q

(9)

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

測量方程為

Zk=Hxk+vk

(10)

P(v)～(0,R)

(11)

cov(v)=E(vvT)=R

(12)

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

2.2 擴展的卡爾曼濾波及模型建立

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

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

(13)

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

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

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

(14)

式(14)中，f(xk|k)進行泰勒展開，

(15)

(16)

測量方程表示為

=Hkxk+Vk

(17)

式(17)中， Gh(xk)表示在xk|kG處進行泰勒展開，

(18)

參考經(jīng)典卡爾曼濾波的推導(dǎo)過程， 可以得到擴展卡爾曼濾波的五個方程式(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ù)測方程， 式(20)是誤差協(xié)方差方程， 式(21)是卡爾曼增益方程， 式(22)是濾波校正方程， 式(23)是誤差協(xié)方差更新方程。

2.3 EKF模型的建立

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

圖3 腦血流動力學數(shù)學模型圖

3 實驗結(jié)果與討論

3.1 擴展的卡爾曼濾波(EKF)對Gauss噪聲濾波有效性的仿真實驗

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

3.2 EKF提取HbR和HbO2濃度變化

3.2.1 實驗裝置

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

圖4 疊加噪聲的原始信號

圖5 EKF濾波的信號

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

3.2.2 Valsava實驗與分析

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

圖7 Valsava實驗血紅蛋白變化曲線圖

3.2.3 視覺誘發(fā)實驗與分析

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

圖8 視覺誘發(fā)源

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

3.3 EKF和EEMD算法比對

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

圖9 視覺誘發(fā)實驗血紅蛋白濃度變化曲線圖

表2 EKF和EEMD算法比對

4 結(jié) 論

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