孫呈霞，吳向東，劉 仙,*

(1.燕山大學(xué) 電氣工程學(xué)院，河北 秦皇島 066004；2.北京理工大學(xué) 自動(dòng)化學(xué)院，北京 100081)

0 引言

由于直接記錄的腦節(jié)律可能會(huì)因?yàn)樯窠?jīng)元和放大器中的噪聲干擾以及記錄裝置中的不確定性而變得不準(zhǔn)確，所以狀態(tài)估計(jì)和參數(shù)辨識(shí)在腦科學(xué)研究中起著十分重要的作用[1]。Deng等人[2]將無(wú)跡卡爾曼濾波器(Unscented Kalman Filter, UKF)算法與一種基于同步的方法相結(jié)合，從被噪聲嚴(yán)重破壞的活動(dòng)電位中估計(jì)出了FitzHugh-Nagumo模型和Hidmarsh-Rose模型的所有未知參數(shù)。Liu等人[3]和Shan等人[4]證明，UKF也可以有效估計(jì)神經(jīng)群模型的不可測(cè)狀態(tài)以及參數(shù)。Lankarany等人[5]提出雙擴(kuò)展卡爾曼濾波器算法和雙無(wú)跡卡爾曼濾波器算法來(lái)重構(gòu)Hodgkin-Huxley模型的動(dòng)力學(xué)，并通過(guò)記錄的噪聲干擾下的膜電壓數(shù)據(jù)辨識(shí)出了模型的未知參數(shù)。隨后，Liao等人[6]證實(shí)，容積卡爾曼濾波器(Cubature Kalman Filter, CKF)及其改進(jìn)濾波算法可以用于估計(jì)Hodgkin-Huxley模型的狀態(tài)和辨識(shí)模型的參數(shù)，并且統(tǒng)計(jì)結(jié)果表明雙容積卡爾曼濾波器方法的辨識(shí)精度更高。Armando等人[7]發(fā)現(xiàn)CKF也適用于估計(jì)癲癇持續(xù)狀態(tài)下的一類神經(jīng)群模型的狀態(tài)和參數(shù)。

如上所述，在神經(jīng)科學(xué)領(lǐng)域，大部分用于辨識(shí)參數(shù)和估計(jì)狀態(tài)的方法是基于模型提出的，所以數(shù)學(xué)模型的選擇也是一個(gè)關(guān)鍵點(diǎn)。完善的數(shù)學(xué)模型可以準(zhǔn)確地模擬真實(shí)的腦節(jié)律，這是我們研究辨識(shí)和估計(jì)方法的前提[8-9]。目前用于模擬腦電圖(Electroencephalogram, EEG)信號(hào)的模型主要有兩大類。一類是從微觀層面上模擬單個(gè)神經(jīng)元的活動(dòng)，如前面提到的FitzHugh-Nagumo模型[10]、Hidmarsh-Rose模型[11]和Hodgkin-Huxley模型[12]。大腦不同區(qū)域的信息往往要通過(guò)耦合許多神經(jīng)元模型來(lái)模擬。這種建模方式的缺點(diǎn)是它需要很強(qiáng)的計(jì)算能力來(lái)求解模擬狀態(tài)。另一類是從宏觀層面上模擬大量神經(jīng)元的平均活動(dòng)，例如神經(jīng)群模型[13-14]。類似于Jansen-Rit模型這樣的宏觀神經(jīng)群模型既簡(jiǎn)單又更具生理意義，是神經(jīng)科學(xué)研究中最常用的模型。之前提到的那些卡爾曼濾波方法在神經(jīng)科學(xué)領(lǐng)域的局限性在于它們無(wú)法有效辨識(shí)突變參數(shù)。已有研究指出，腦功能障礙可能是由于大腦系統(tǒng)內(nèi)部參數(shù)之間的平衡關(guān)系被破壞而引起的內(nèi)部失穩(wěn)所致[15-16]。當(dāng)局部參數(shù)突然發(fā)生變化且超出了系統(tǒng)的自我調(diào)整范疇時(shí)，該區(qū)域可能會(huì)成為一個(gè)病灶源，產(chǎn)生的異常節(jié)律在耦合時(shí)延的作用下傳播到其他區(qū)域，引發(fā)嚴(yán)重的神經(jīng)疾病。因此突變參數(shù)的辨識(shí)對(duì)于神經(jīng)疾病的預(yù)防和治療來(lái)說(shuō)意義重大。Li等人[17]利用強(qiáng)跟蹤容積卡爾曼濾波器(Strong Tracking Cubature Kalman Filter, STCKF)來(lái)完成脈沖機(jī)動(dòng)衛(wèi)星系統(tǒng)中的實(shí)時(shí)估計(jì)和辨識(shí)任務(wù)。STCKF算法通過(guò)在狀態(tài)誤差協(xié)方差矩陣中引入漸消因子來(lái)調(diào)整卡爾曼增益矩陣，從而解決了非線性系統(tǒng)中突變參數(shù)的跟蹤問(wèn)題。STCKF為有效辨識(shí)腦系統(tǒng)中的突變參數(shù)創(chuàng)造了可能。

基于以上原因，本文嘗試?yán)肧TCKF算法來(lái)估計(jì)一類耦合Jansen-Rit模型中的突變興奮性增益參數(shù)。這項(xiàng)工作有望為腦部疾病病灶源的確定提供可能的方法。

1 Jansen-Rit模型

本文選用Jansen和Rit[18]建立的神經(jīng)群模型作為基礎(chǔ)模型來(lái)模擬實(shí)際腦活動(dòng)背景下的不同節(jié)律。由于大多數(shù)實(shí)驗(yàn)數(shù)據(jù)(如由腦電圖記錄的EEG信號(hào))反映的是大量神經(jīng)元的集群活動(dòng)，所以這類完全來(lái)源于現(xiàn)象學(xué)的宏觀模型對(duì)于研究大腦機(jī)制等非常有效[19]。此外，與微觀模型相比，這類模型對(duì)計(jì)算能力的要求不是很高。Jansen-Rit模型可以由8個(gè)一階微分方程描述：

(1)

表1 模型參數(shù)的生理意義和標(biāo)準(zhǔn)值

Tab.1 Physiological meanings and standard values of the model parameters

參數(shù)生理意義標(biāo)準(zhǔn)值He平均興奮性突觸增益3.25mVHi平均抑制性突觸增益22mVτe興奮性膜平均時(shí)間常數(shù)與樹(shù)突平均時(shí)延0.0108sτi抑制性膜平均時(shí)間常數(shù)與樹(shù)突平均時(shí)延0.02sτd來(lái)源于某些群落的傳出連接平均時(shí)延0.0303sC1主體細(xì)胞對(duì)興奮反饋回路中中間神經(jīng)元樹(shù)突產(chǎn)生的平均突觸連接數(shù)135C2中間神經(jīng)元對(duì)興奮反饋回路中主體細(xì)胞樹(shù)突產(chǎn)生的平均突觸連接數(shù)108C3主體細(xì)胞對(duì)抑制反饋回路中中間神經(jīng)元樹(shù)突產(chǎn)生的平均突觸連接數(shù)33.75C4中間神經(jīng)元對(duì)抑制反饋回路中主體細(xì)胞樹(shù)突產(chǎn)生的平均突觸連接數(shù)33.752e0最大點(diǎn)燃率5s-1v0對(duì)應(yīng)于點(diǎn)燃率的突觸后電位6mVrSigmoid變換的斜率0.56mV-1

已有研究表明，異常腦節(jié)律可能是某些神經(jīng)群落中興奮性參數(shù)He突然增大引發(fā)的，然后在耦合延遲的作用下棘波信號(hào)傳播到其他群落。這表明，局部He的突然變化可能是腦功能障礙發(fā)生的根源。本文希望通過(guò)及時(shí)捕捉He的突然變化，為尋找病灶源提供一種可能的方法。

2 突變參數(shù)辨識(shí)方法

一類受到測(cè)量噪聲干擾的離散非線性狀態(tài)空間模型可以描述為[20]

(2)

其中，下標(biāo)k表示采樣時(shí)間點(diǎn)；xk∈RNx是采樣時(shí)間點(diǎn)k處的狀態(tài)向量，Nx表示該狀態(tài)向量的維度；pk-1∈RNP是一個(gè)已知的隨機(jī)輸入，NP表示向量維度；yok∈RNy是測(cè)量向量，Ny表示向量維度；f(·):RNx+Nu→RNy和h(·):RNx→RNy分別表示非線性動(dòng)力學(xué)函數(shù)和測(cè)量函數(shù)；wk表示測(cè)量噪聲，服從均值為零方差為WK的高斯分布。

1)預(yù)測(cè)過(guò)程

該預(yù)測(cè)過(guò)程與CKF算法[22]的預(yù)測(cè)過(guò)程相同。具體描述為

(3)

其中上標(biāo)*表示引入漸消因子和遺忘因子前的情況；Uk-1表示采樣時(shí)間點(diǎn)k-1處的后驗(yàn)估計(jì)狀態(tài)誤差協(xié)方差矩陣Pk|k-1的喬利斯基分解因子；Chol(·)表示喬利斯基分解函數(shù)；φi為求容積點(diǎn)；ζ是一組標(biāo)準(zhǔn)正交點(diǎn)集，集合中的元素表示為

2)計(jì)算漸消因子λk

受強(qiáng)跟蹤濾波器(Strong Tracking Filter, STF)算法[23]的啟發(fā)，STCKF算法通過(guò)在CKF算法中引入漸消因子k和遺忘因子來(lái)提高對(duì)突變狀態(tài)的跟蹤能力。計(jì)算漸消因子λk的過(guò)程可以描述為

(4)

計(jì)算漸消因子λk:

(5)

其中，εk表示殘差序列；yok和Wk與前文一致，分別表示測(cè)量向量和測(cè)量噪聲的方差；Vk是一個(gè)與殘差相關(guān)的中間函數(shù)；Nk和Mk是用于計(jì)算漸消因子λk的矩陣，而tr[Nk]和tr[Mk]分別為矩陣Nk和Mk的跡；ρ是殘差序列的遺忘因子，并且0≤ρ≤1，通常取ρ=0.95；β為弱化因子，通常β≥1。

3)校正過(guò)程

該過(guò)程涉及到一個(gè)在計(jì)算出漸消因子λk之后重新計(jì)算先驗(yàn)估計(jì)狀態(tài)誤差協(xié)方差矩陣Pk|k-1的步驟：

(6)

校正過(guò)程的其余公式與CKF相同：

(7)

除了估計(jì)系統(tǒng)狀態(tài)之外，CKF的另一個(gè)重要特性是，當(dāng)與增廣狀態(tài)向量法[20]相結(jié)合時(shí)可以用于辨識(shí)未知參數(shù)，STCKF亦是如此。增廣向量法通過(guò)將未知參數(shù)視為虛擬狀態(tài)，把系統(tǒng)表達(dá)式改寫為

(8)

其中qk∈RNq表示參數(shù)向量，狀態(tài)向量由常數(shù)參數(shù)向量增廣。在濾波器的更新作用下，即使沒(méi)有相應(yīng)的隨機(jī)終止參數(shù)方程，辨識(shí)得到的參數(shù)向量qk仍會(huì)隨時(shí)間變化。

通過(guò)將衰減因子λk和遺忘因子ρ引入到CKF算法，來(lái)近似一個(gè)給定非線性系統(tǒng)的狀態(tài)誤差協(xié)方差矩陣，可以有效克服STF迭代過(guò)程中Jacobi矩陣的逼近精度低和計(jì)算繁瑣的問(wèn)題。同時(shí)，STCKF算法繼承了STF對(duì)突變狀態(tài)的強(qiáng)跟蹤能力。

3 數(shù)值結(jié)果

本文利用一類3個(gè)節(jié)點(diǎn)耦合的Jansen-Rit模型來(lái)模擬真實(shí)的EEG信號(hào)，其連接結(jié)構(gòu)如圖1所示。這類模型中每個(gè)群落的數(shù)學(xué)表達(dá)如式(1)所示，其中N=3。受到測(cè)量噪聲干擾的3個(gè)節(jié)點(diǎn)耦合的Jansen-Rit模型輸出可以描述為

y(t)=h[x(t)]+w(t),

(9)

圖1 3個(gè)節(jié)點(diǎn)耦合的Jansen-Rit模型的連接圖

Fig.1 Connection diagram of a three-node coupled Jansen-Rit mode

第一組仿真實(shí)驗(yàn)中，群落1的興奮性參數(shù)He的變化過(guò)程設(shè)置為3.6 mV→4 mV→3.4 mV→3.5 mV。在實(shí)際的腦系統(tǒng)中，引起腦疾病發(fā)作的局部參數(shù)突變可能不會(huì)像這樣連續(xù)發(fā)生，而通常是從標(biāo)準(zhǔn)值變?yōu)橐粋€(gè)非標(biāo)準(zhǔn)值。設(shè)置多次突變的目的是驗(yàn)證STCKF對(duì)突變參數(shù)的強(qiáng)跟蹤性能。由于本文的主要目的是辨識(shí)突變參數(shù)，并且群落2和3的狀態(tài)估計(jì)結(jié)果與群落1相似，本節(jié)只給出了群落1的相關(guān)結(jié)果。群落1的動(dòng)力學(xué)特性如圖2所示，可以看出模型輸出的峰值和頻率也會(huì)隨著興奮性增益值的變化而變化。

圖2 第一組模擬實(shí)驗(yàn)中群落1的動(dòng)力學(xué)

Fig.2 The dynamics of population 1 in the first set of simulation experiments

這里引入均方根誤差(Root Mean Square Error,RMSE)[24]的概念來(lái)檢驗(yàn)濾波算法的精度。計(jì)算第i個(gè)狀態(tài)分量或局部參數(shù)的RMSE，公式為

(10)

圖3顯示的是測(cè)量噪聲影響下，經(jīng)過(guò)CKF和STCKF的濾波處理后得到的輸出估計(jì)和參數(shù)辨識(shí)結(jié)果，其中實(shí)線表示實(shí)際值，虛線表示輸出估計(jì)結(jié)果，點(diǎn)劃線表示參數(shù)辨識(shí)結(jié)果。圖3(a)和(b)給出了兩種濾波算法的輸出估計(jì)結(jié)果，可以看出這兩種濾波算法都可以準(zhǔn)確估計(jì)群落1的輸出。為了更直觀地評(píng)價(jià)這兩種濾波算法的精度，計(jì)算輸出的RMSE值，得到R(y1)CKF=0.255 7和R(y1)STCKF=0.211 9，說(shuō)明STCKF濾波算法的估計(jì)精度高于CKF。圖3(c)和(d)給出了兩種濾波算法的參數(shù)辨識(shí)結(jié)果，可以看出當(dāng)辨識(shí)結(jié)果已經(jīng)達(dá)到穩(wěn)定狀態(tài)時(shí)，CKF算法幾乎失去對(duì)突變參數(shù)的跟蹤能力，而STCKF算法能夠訊速而準(zhǔn)確地辨識(shí)出新的參數(shù)值。相關(guān)結(jié)果表明，STCKF算法不僅可以辨識(shí)出耦合Jansen-Rit模型中的突變興奮性增益參數(shù)He，而且在估計(jì)精度方面優(yōu)于CKF算法。

圖3 第一組仿真實(shí)驗(yàn)得到的輸出估計(jì)和參數(shù)辨識(shí)的結(jié)果

Fig.3 Output estimation and parameter identification results from the first set of simulation experiments

第二組仿真實(shí)驗(yàn)中，將群落1的興奮性參數(shù)He的變化過(guò)程設(shè)置為3.25 mV→3.5 mV→3.25 mV→3.6 mV→3.25 mV→3.4 mV→3.25 mV。此時(shí)，群落1的動(dòng)態(tài)特性如圖4所示。圖5是經(jīng)CKF和STCKF算法處理后得到的輸出估計(jì)和參數(shù)辨識(shí)結(jié)果。與圖3一致，實(shí)線表示實(shí)際值，虛線表示輸出估計(jì)結(jié)果，點(diǎn)劃線表示參數(shù)辨識(shí)結(jié)果。計(jì)算得到R(y1)CKF=0.195 0以及R(y1)STCKF=0.174 1。因此，可以得出與第一組實(shí)驗(yàn)一致的結(jié)論。也進(jìn)一步證明了STCKF算法在突變參數(shù)辨識(shí)和估計(jì)精度方面的優(yōu)越性。

圖4 第二組模擬實(shí)驗(yàn)中群落1的動(dòng)力學(xué)

Fig.4 The dynamics of population 1 in the second set of simulation experiments

圖5 第二組仿真實(shí)驗(yàn)得到的輸出估計(jì)和參數(shù)辨識(shí)的結(jié)果

Fig.5 Output estimation and parameter identification results from the second set of simulation experiments

4 結(jié)論

本文針對(duì)耦合Jansen-Rit模型中參數(shù)突變的特點(diǎn)，給出了一種結(jié)合增廣向量法的STCKF方法。該方法不僅可以估計(jì)模型的輸出，也可以辨識(shí)模型中的突變參數(shù)。為了突出所給方法的優(yōu)越性，將其與CKF方法進(jìn)行比較。實(shí)驗(yàn)結(jié)果表明：1)CKF和STCKF算法都可以準(zhǔn)確估計(jì)耦合Jansen-Rit模型的輸出，并且通過(guò)計(jì)算輸出的RMSE值，證明STCKF算法的估計(jì)精度高于CKF；2)CKF和STCKF算法都可以用于辨識(shí)一類耦合Jansen-Rit模型的參數(shù)，但是當(dāng)參數(shù)突然變化時(shí)CKF算法幾乎失去對(duì)參數(shù)的跟蹤能力，而STCKF算法仍能夠迅速而準(zhǔn)確地辨識(shí)出發(fā)生了變化的參數(shù)值。本文工作有望成為進(jìn)一步研究腦疾病抑制方法的良好開(kāi)端。但是，仍有一些問(wèn)題需要進(jìn)一步探討。例如，如何將這項(xiàng)研究與控制方案結(jié)合起來(lái)，以系統(tǒng)地減輕甚至消除神經(jīng)疾病。這也是我們未來(lái)的工作重點(diǎn)。