卡爾曼濾波器在海馬場(chǎng)電位ripple節(jié)律分析中的應(yīng)用

張 櫨， 林龍年

（華東師范大學(xué) 腦功能基因組學(xué)教育部、上海市重點(diǎn)實(shí)驗(yàn)室，上海 200062）

卡爾曼濾波器在海馬場(chǎng)電位ripple節(jié)律分析中的應(yīng)用

張 櫨， 林龍年

（華東師范大學(xué) 腦功能基因組學(xué)教育部、上海市重點(diǎn)實(shí)驗(yàn)室，上海 200062）

利用自適應(yīng)自回歸（adaptive autoregressive，AAR）模型和卡爾曼濾波器算法，分析小鼠海馬CA1區(qū)場(chǎng)電位ripple高頻振蕩的時(shí)頻特性.研究發(fā)現(xiàn)，與傳統(tǒng)的基于短時(shí)傅立葉變換的實(shí)時(shí)頻譜分析方法相比，利用AAR模型以及卡爾曼濾波器算法的參數(shù)化方法在對(duì)ripple高頻振蕩信號(hào)進(jìn)行實(shí)時(shí)頻譜分析時(shí)，具有更高的時(shí)域和頻域分辨率.因此，基于卡爾曼濾波器得到的ripple能量變化，可更為準(zhǔn)確、實(shí)時(shí)地反映ripple高頻振蕩的發(fā)生與變化過程.

卡爾曼濾波器； 海馬CA1區(qū)； ripple； 實(shí)時(shí)頻譜分析

0 引 言

在中樞神經(jīng)系統(tǒng)中記錄到的局部場(chǎng)電位活動(dòng)，具有多種行為相關(guān)性的節(jié)律性振蕩，而在大腦的海馬結(jié)構(gòu)中，這些節(jié)律性振蕩尤為明顯.在大鼠海馬一般可記錄到4種腦電（EEG）節(jié)律振蕩，按頻段依次為theta節(jié)律（4～12 Hz）、beta節(jié)律（12～30 Hz）、gamma節(jié)律（30～100 Hz）和ripple節(jié)律（100～200 Hz）［1］.這些節(jié)律性電活動(dòng)隨著動(dòng)物行為狀態(tài)的不同而強(qiáng)弱交替，構(gòu)成了一種動(dòng)態(tài)平衡.當(dāng)動(dòng)物處于行走、探索狀態(tài)以及快波睡眠時(shí)，場(chǎng)電位（local field potential，LFP）呈現(xiàn)出大幅度且有規(guī)律的theta振蕩并且伴隨著gamma節(jié)律；而當(dāng)動(dòng)物處于慢波睡眠（slow wave sleep）、安靜站立等狀態(tài)時(shí)，在海馬CA1區(qū)則可記錄到高頻（120～200 Hz）的場(chǎng)電位振蕩［2，3］.這樣的高頻場(chǎng)電位節(jié)律振蕩被稱作ripple，一般在海馬CA1區(qū)的錐體層幅度最大.在功能上，ripple高頻節(jié)律被認(rèn)為與記憶的鞏固過程有關(guān).最新的研究發(fā)現(xiàn)，海馬中存在大量的位置細(xì)胞，一個(gè)位置細(xì)胞會(huì)在動(dòng)物接近或穿過某一特定位置的時(shí)候增強(qiáng)活動(dòng)，因此新環(huán)境中的空間信息由位置細(xì)胞的活動(dòng)狀態(tài)所表征［4］，而具有相交或相鄰位置域的多個(gè)位置細(xì)胞，會(huì)出現(xiàn)放電的相關(guān)性以及有序性，這種關(guān)系在動(dòng)物慢波睡眠階段可以重現(xiàn)，尤其在ripple節(jié)律發(fā)生的時(shí)候［5－8］.每個(gè)ripple的時(shí)程為50～100 ms，而一系列位置細(xì)胞的有序放電又在如此短的時(shí)間內(nèi)重演，因此ripple的精確偵測(cè)以及對(duì)其動(dòng)力學(xué)變化的分析，對(duì)于進(jìn)一步理解與研究海馬多個(gè)神經(jīng)元在記憶鞏固時(shí)期編碼的機(jī)制有著重要的意義.實(shí)時(shí)頻譜分析可以同時(shí)給出時(shí)間序列信號(hào)在時(shí)域與頻域的動(dòng)力學(xué)描述，因此，可以作為研究ripple高頻振蕩及其他場(chǎng)電位節(jié)律的有力工具.傳統(tǒng)的時(shí)頻譜分析方法為譜圖（spectrogram），即將信號(hào)分做很多重疊的時(shí)段，對(duì)每個(gè)時(shí)段分別應(yīng)用短時(shí)傅立葉變換（short-time Fourier transform，STFT）.其原理在于傅立葉變換將信號(hào)放到一個(gè)以正弦與余弦函數(shù)為基函數(shù)的空間，而這些基函數(shù)彼此正交，于是將原信號(hào)分解成不同頻率成分.但這種方法卻面臨著如何在時(shí)域與頻域分辨率上的協(xié)調(diào)問題（不確定性原理）［9］.因?yàn)樘岣邥r(shí)間分辨率則必然導(dǎo)致低頻率分辨率，反之亦然.而應(yīng)用自適應(yīng)自回歸模型（adaptive autoregressive（AAR）model）來描述信號(hào)，則可以得到更高的頻域分辨率.再應(yīng)用卡爾曼濾波器算法來構(gòu)建AAR模型進(jìn)行實(shí)時(shí)頻譜分析［10－13］，則可進(jìn)一步在保持頻域高分辨率的同時(shí)，得到與信號(hào)采樣頻率一樣精度的時(shí)間分辨率，且無需像STFT一樣提取大量數(shù)據(jù)以得到自相關(guān)序列.本文運(yùn)用卡爾曼濾波算法，對(duì)海馬CA1場(chǎng)電位的ripple高頻振蕩進(jìn)行實(shí)時(shí)頻譜分析，并與傳統(tǒng)的基于STFT的實(shí)時(shí)頻譜分析方法進(jìn)行比較.

1 方 法

1.1 功率譜密度

根據(jù)Cramer、Kolmogoroff和 Wiener等人的結(jié)果［14，15］，協(xié)方差穩(wěn)定序列的自協(xié)方差可以表示為

式（1）為自協(xié)方差函數(shù)的頻率表示，F(xiàn)（ω）為頻譜分布函數(shù).如果Xt為非確定性序列（即不可被過去已知值所準(zhǔn)確預(yù)測(cè)），則有dF（ω）＝f（ω）.f（ω）絕對(duì)連續(xù)，被稱作是功率譜密度.f（ω）dω刻畫了信號(hào)在（ω，ω＋dω）頻段對(duì)整個(gè)信號(hào)方差的貢獻(xiàn).因此，f（ω）出現(xiàn)某個(gè)峰值的時(shí)候，意味著該信號(hào)的相應(yīng)頻段可能非常重要.

記Xt的階數(shù)m的自回歸（AR）模型為

其中E（t）為時(shí)不相關(guān)偏差，其協(xié)方差矩陣為∑.

在方程兩端使用傅利葉變換（Fourier transform），可以得到

其中，H（ω）＝A－1（ω）為傳遞函數(shù)，＊為共軛轉(zhuǎn)制.

1.2 自回歸模型各參數(shù)的確立

自回歸（AR）模型（2）中，右乘Xt－k，其中k＝1，2，…，m，然后在兩端取期望，于是得到Y(jié)ule-Walker方程

其中，μ（q）＝〈XtXt＋q〉為時(shí)間延遲q的協(xié)方差矩陣.

其中N為樣本數(shù).

1.3 自適應(yīng)自回歸模型和卡爾曼濾波器

注意AR模型（2）中的系數(shù)Ak是獨(dú)立于時(shí)間的，當(dāng)應(yīng)用在場(chǎng)電位的theta振蕩中時(shí)是合適的，因?yàn)閠heta振蕩持續(xù)時(shí)間長(zhǎng)且是連續(xù)性的，往往在10 s以上.但當(dāng)應(yīng)用于轉(zhuǎn)瞬性很強(qiáng)的ripple振蕩的時(shí)候，就會(huì)出現(xiàn)問題.這是因?yàn)閞ipple往往發(fā)生在大幅尖波之上，而這些尖波的突然產(chǎn)生與維持并沒有連續(xù)性，ripple持續(xù)的時(shí)間往往只有50～100 ms，模型所要求的平穩(wěn)性條件往往達(dá)不到.在這種情況下，有必要考慮隨時(shí)間變化的模型參數(shù).于是考慮如下自適應(yīng)自回歸模型（AAR）模型，

其中，E（t）為零均值高斯噪聲過程，方差為∑t，此模型與自回歸模型（2）的根本區(qū)別在于前者的模型參數(shù)隨時(shí)間變化而變化.

在不給出證明的情況下，這里給出卡爾曼濾波器方程如下，其理論推導(dǎo)的細(xì)節(jié)詳見文獻(xiàn)［12］.

這樣的估計(jì)更加光滑，也被稱作平滑器（smoother）.據(jù)此可以得到新的參數(shù)估計(jì)［19］，

2 結(jié) 果

2.1 不同行為狀態(tài)下小鼠海馬CA1區(qū)場(chǎng)電位的基本特征

對(duì)不同行為狀態(tài)下，小鼠海馬CA1區(qū)記錄到的場(chǎng)電位進(jìn)行功率譜密度分析，結(jié)果顯示，在探索活動(dòng)以及REM睡眠狀態(tài)下，場(chǎng)電位表現(xiàn)出規(guī)律的theta節(jié)律振蕩，同時(shí)在theta波峰處，出現(xiàn)gamma節(jié)律快速振蕩（見圖1A和C）.另外，在REM睡眠狀態(tài)下記錄到的theta節(jié)律相對(duì)更為規(guī)則（見圖1C）.而在動(dòng)物探索活動(dòng)的時(shí)候，一般還有伴有beta節(jié)律（見圖1A）.在慢波睡眠狀態(tài)，則以在大幅尖波上出現(xiàn)的ripple節(jié)律高頻振蕩為主 （見圖1B）.

圖1 海馬CA1區(qū)場(chǎng)電位及其功率譜密度分布Fig.1 LFP signals recorded in hippocampal CA1 and power spectral density disctribution

2.2 基于卡爾曼濾波器的實(shí)時(shí)頻譜分析對(duì)ripple的偵測(cè)

實(shí)時(shí)頻譜分析可以準(zhǔn)確地描述場(chǎng)電位ripple高頻頻段的發(fā)生以及在頻域的動(dòng)態(tài)變化，當(dāng)ripple發(fā)生時(shí)，在100～200 Hz的頻域上可以看到能量相應(yīng)增強(qiáng)（見圖2 B和C）.本文對(duì)傳統(tǒng)的基于短時(shí)傅立葉變換（STFT）的譜圖（spectrogram）（見圖2B）以及基于卡爾曼濾波器算法的實(shí)時(shí)頻譜分析圖進(jìn)行了對(duì)比分析（見圖2C）.對(duì)比相應(yīng)的原始場(chǎng)電位（見圖2A），可以發(fā)現(xiàn)基于卡爾曼濾波器算法的實(shí)時(shí)頻譜分析在頻域以及時(shí)域上具有更高的分辨率.同時(shí)，傳統(tǒng)的基于短時(shí)傅立葉變換（STFT）的譜圖由于必須對(duì)信號(hào)進(jìn)行分段加窗處理（此處為50 ms的hanning窗，1 ms的時(shí)間步長(zhǎng)），并采用時(shí)間窗中的中位時(shí)間作為該窗的時(shí)刻對(duì)應(yīng)點(diǎn)因此，因此對(duì)ripple波的偵測(cè)有一個(gè)開始時(shí)的時(shí)間提前量和結(jié)束時(shí)的時(shí)間延遲量.而卡爾曼濾波器本身是AAR模型的一種以采樣時(shí)間為時(shí)間步長(zhǎng)的自適應(yīng)迭代算法，再應(yīng)用平滑器對(duì)迭代進(jìn)行修正，從而在偵測(cè)ripple波時(shí)，開始的時(shí)間提前量以及結(jié)束的時(shí)間延遲量均很小.

為了進(jìn)一步定量比較這兩種方法對(duì)ripple波偵測(cè)準(zhǔn)確性的差異，將每個(gè)時(shí)間點(diǎn)ripple段的能量求和，作為ripple能量隨時(shí)間變化的表示（見圖3C）.同樣的，基于SFFT的實(shí)時(shí)頻譜分析也可以得到ripple能量隨時(shí)間的變化過程（見圖3D）.將得到的這兩組ripple能量變化曲線，分別與場(chǎng)電位ripple濾波信號(hào)（見圖3B細(xì)線）的幅度（見圖3B粗線）相比較，從而得到兩組相關(guān)系數(shù)ρAAR與ρSFFT.由于ripple節(jié)律隨時(shí)間的強(qiáng)弱變化可以由場(chǎng)電位ripple濾波信號(hào)的幅度（見圖3B粗線）進(jìn)行描述，因此該相關(guān)系數(shù)越接近1，說明相應(yīng)的實(shí)時(shí)頻譜分析越能準(zhǔn)確表達(dá)場(chǎng)電位ripple段的能量變化.

10組不同小鼠的場(chǎng)電位記錄（每組3段，各5 s時(shí)長(zhǎng)）數(shù)據(jù)分析顯示，ρAAR（mean＝0.677，median＝0.708）高于ρSFFT（mean＝0.630，median＝0.671），（見圖4A，p?0.001，Wilcoxon signed rank test）.相關(guān)系數(shù)指數(shù)（見圖4B）亦顯示ρAAR＞ρSFFT（p?0.001，Student test）.以上結(jié)果表明，基于ARR模型的實(shí)時(shí)頻譜分析在偵測(cè)場(chǎng)電位ripple段的能量變化時(shí)，較STFT方法更為精確.

圖4 STFT方法與卡爾曼濾波器方法偵測(cè)ripple振蕩的效果比較Fig.4 Compartment of ripple detection through STFT and Kalman filter

3 討 論

AAR模型及卡爾曼濾波器算法中涉及到許多參數(shù)的設(shè)置.首先，AAR模型應(yīng)用于非平穩(wěn)序列的階數(shù)，不同于AR模型對(duì)平穩(wěn)序列可以利用AIC等準(zhǔn)則進(jìn)行確定［13］.其次，模型描述數(shù)據(jù)的準(zhǔn)確性受到狀態(tài)更新噪音為以及觀測(cè)噪音的影響，而卡爾曼濾波器是一個(gè)迭代算法，因此各參數(shù)的初始值的不同設(shè)置將影響到模型估計(jì)的結(jié)果.對(duì)于此，并沒有完整的解決方案［13，20，21］，參考已有的研究工作，以及考慮到場(chǎng)電位數(shù)據(jù)本身的特性，我們將模型階數(shù)p設(shè)為9；初始觀測(cè)噪音為0.001，初始狀態(tài)更新噪音為0.000 1.

研究結(jié)果顯示，基于STFT的實(shí)時(shí)頻譜分析，無法在時(shí)域和頻域同時(shí)高精度地表示場(chǎng)電位的ripple節(jié)律.而基于ARR模型以及卡爾曼濾波器算法得到的信號(hào)實(shí)時(shí)頻譜分析，提供了一個(gè)比傳統(tǒng)STFT更為有力的工具.該法在腦電（EEG）的分析方法中已有一些應(yīng)用［22，23］，而平滑器的應(yīng)用不但讓得到的頻譜更加光滑，而且消除了對(duì)信號(hào)表征的時(shí)間延遲［13］；另一方面，平滑器是基于信號(hào)的離線分析為前提，當(dāng)進(jìn)行信號(hào)的在線解碼時(shí)，便無法應(yīng)用.Wilson的研究小組，基于二階ARR模型與卡爾曼濾波器的算法，偵測(cè)海馬區(qū)場(chǎng)電位ripple信號(hào)的發(fā)生與起始，得到了比傳統(tǒng)的運(yùn)用ripple濾波信號(hào)的幅度（通過希爾伯特變換）偵測(cè)方法更好的效果［21］.而我們將原始場(chǎng)電位直接用ARR模型直接描述，由此可以得到被Nyquist frequency限制的大范圍的頻譜，更準(zhǔn)確地描述了場(chǎng)電位ripple頻段頻率以及能量大小的變化.

AR模型同ARR模型一樣也是參數(shù)化方法，但前者中各個(gè)參數(shù)不隨時(shí)間變化而變化，并要求信號(hào)是平穩(wěn)信號(hào).一般認(rèn)為12 s以上的EEG信號(hào)是平穩(wěn)的［24］，而發(fā)生在尖波之上持續(xù)時(shí)間僅為50～100 ms的ripple信號(hào)則難以達(dá)到這一要求.對(duì)于低頻的theta振蕩，由于其持續(xù)時(shí)間較長(zhǎng)（往往達(dá)到數(shù)十秒以上），而且其節(jié)律振蕩的子周期也在120 ms左右.因此不但可以應(yīng)用AR模型來進(jìn)行實(shí)時(shí)頻譜分析，傳統(tǒng)的SFFT也可以應(yīng)用較大的時(shí)間窗來得到較高的頻域分辨率，從而使方法不同造成的相對(duì)較低的時(shí)間分辨率可以被忽略.

［1］ ANDERSEN P.The Hippocampus Book［M］.Oxford：Oxford University Press，2007.

［2］ BUZSAKI G，HORVATH Z，URIOSTE R，et al.High-frequency network oscillation in the hippocampus［J］.Science，1992，256：1025-1027.

［3］ O′KEEFE J，NADEL L.The Hippocampus as a Cognitive Map［M］.Oxford：Clarendon Press，1978.

［4］ O′KEEFE J，DOSTROYSKY J.The hippocampus as a spatial map.Preliminary evidence from unit activity in the freely-moving rat［J］.Brain Res，1971，34：171-175.

［5］ WILSON M A，MCNAUGHTON B L.Reactivation of hippocampal ensemble memories during sleep［J］.Science，1994，265：676-679.

［6］ LEE A K，WILSON M A.Memory of sequential experience in the hippocampus during slow wave sleep［J］.Neuron，2002，36：1183-1194.

［7］ FOSTER D J，WILSON M A.Reverse replay of behavioural sequences in hippocampal place cells during the awake state［J］.Nature，2006，440：680-683.

［8］ O＇NEILL J，SENIOR T J，ALLEN K，et al.Reactivation of experience-dependent cell assembly patterns in the hippocampus［J］.Nat Neurosci，2008（11）：209-215.

［9］ BENEDICKS M.On fourier transforms of functions supported on sets of finite Lebesgue measure［J］.Journal of Mathematical Analysis and Applications，1985，106：180-183.

［10］ KALMAN R E.A new approach to linear filtering and prediction problems［J］.Transactions of the ASME，1960：35-45.

［11］ KALMAN R E，BUCY R S.New results in linear filtering and prediction theory［J］.Transactions of the ASME Series D，Journal of Basic Engineering，1961，83：95-107.

［12］ ANDERSON B D O，MOORE J B.Optimal Filtering［M］.New Jersey：Prentice Hall，1979.

［13］ TARVAINEN M P，HILTUNEN J K，RANTA-AHO P O，et al.Estimation of nonstationary EEG with Kalman smoother approach：an application to event-related synchronization （ERS）［J］.Biomedical Engineering，IEEE Transactions on，2004，51：516-524.

［14］ CRAMER H.On the theory of stationary random processes［J］.The Annals of Mathematics，1940，41：215-230.

［15］ WIENER N.Extrapolation，Interpolation，and Smoothing of Stationary Time Series，with Engineering Applications［M］.Cambridge：Technology Press of the Massachusetts Institute of Technology，1949.

［16］ MORF M，VIEIRA A，LEE D T L，et al.Recursive multichannel maximum entropy spectral estimation［J］.Geoscience Electronics，IEEE Transactions on，1978，16：85-94.

［17］ DING M，BRESSLER S L，YANG W，et al.Short-window spectral analysis of cortical event-related potentials by adaptive multivariate autoregressive modeling：data preprocessing，model validation，and variability assessment［J］.Biological Cybernetics，2000，83：35-45.

［18］ AKAIKE H.A new look at the statistical model identification［J］.Automatic Control，IEEE Transactions on，2003，19：716-723.

［19］ HAMITON J D.Time Series Analysis［M］.Princeton：Princeton University Press，1994.

［20］ ALOIS S.The electroencephalogram and the adaptive autoregressive model：theory and applications［D］.Technischen Universitat Graz，2000.

［21］ NGUYEN D P，KLOOSTERMAN F，BARBIERI R，et al.Characterizing the dynamic frequency structure of fast oscillations in the rodent hippocampus［J］.Front Integr Neurosci，2009（3）：11-11.

［22］ BOHLIN T.Analysis of EEG signals with changing spectra using a short-word Kalman estimator［J］.Mathematical Biosciences，1977，35：221-259.

［23］ AMOLD M，MILNER X H R，WITTE H，et al.Adaptive AR modeling of nonstationary time series by means of Kalman filtering［J］.Biomedical Engineering，IEEE Transactions on，1998，45：553-562.

［24］ COHEN B A，SANCES A.Stationarity of the human electroendephalogram［J］.Med Biol Eng Comput，1977，15：513-518.

Analysis of hippocampal ripple osillations by application of Kalman filter

ZHANG Lu， LIN Long-nian

（Key Laboratory of Brain Functional Genomics，Ministry of Education and Shanghai，Institute of Brain Functional Genomics，East China Normal University，Shanghai 200062，China）

This paper studied high frequency ripple（100～200 Hz）oscillations in hippocanpal CA1 area by applications of adaptive autoregressive（AAR）model and Kalman filter.Compared with traditional real time frequency analysis of time seies based on short term Fourier transfrom（STFT），improved time and frequency resolutions in time-frequency representation could be achieved by parametric methold obtained by AAR model and Kalman filter algorithms.Thus，the occurance of ripple oscillations and the variation of ripple power could be addresed more accurate by Kalman filter than that of STFT.

Kalman filter； hippocampal CA1； ripple； real-time frequency analysis

Q6

A

10.3969／j.issn.1000-5641.2011.06.010

1000-5641（2011）06-0081-08

2010-11

上海市教育委員會(huì)科研創(chuàng)新項(xiàng)目（09ZZ44）；國家自然科學(xué)基金重大項(xiàng)目（30990262）

張櫨，男，博士研究生.E-mail：math2437＠hotmail.com.

林龍年，男，教授.E-mail：lnlin＠brain.ecnu.edu.cn.