一類耦合Jansen-Rit模型中突變參數(shù)的非線性辨識

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

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

0 引言

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

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

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

1 Jansen-Rit模型

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

(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í)延0.0108sτi抑制性膜平均時(shí)間常數(shù)與樹突平均時(shí)延0.02sτd來源于某些群落的傳出連接平均時(shí)延0.0303sC1主體細(xì)胞對興奮反饋回路中中間神經(jīng)元樹突產(chǎn)生的平均突觸連接數(shù)135C2中間神經(jīng)元對興奮反饋回路中主體細(xì)胞樹突產(chǎn)生的平均突觸連接數(shù)108C3主體細(xì)胞對抑制反饋回路中中間神經(jīng)元樹突產(chǎn)生的平均突觸連接數(shù)33.75C4中間神經(jīng)元對抑制反饋回路中主體細(xì)胞樹突產(chǎn)生的平均突觸連接數(shù)33.752e0最大點(diǎn)燃率5s-1v0對應(yīng)于點(diǎn)燃率的突觸后電位6mVrSigmoid變換的斜率0.56mV-1

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

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

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

(2)

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

1)預(yù)測過程

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

(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算法通過在CKF算法中引入漸消因子k和遺忘因子來提高對突變狀態(tài)的跟蹤能力。計(jì)算漸消因子λk的過程可以描述為

(4)

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

(5)

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

3)校正過程

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

(6)

校正過程的其余公式與CKF相同：

(7)

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

(8)

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

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

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

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

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

(9)

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

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

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

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

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

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

(10)

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

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

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

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

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

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

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

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

4 結(jié)論

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