張 櫨， 林龍年

（華東師范大學 腦功能基因組學教育部、上海市重點實驗室，上海 200062）

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

張 櫨， 林龍年

（華東師范大學 腦功能基因組學教育部、上海市重點實驗室，上海 200062）

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

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

0 引 言

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

1 方 法

1.1 功率譜密度

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

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

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

其中E（t）為時不相關(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〉為時間延遲q的協(xié)方差矩陣.

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

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

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

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

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

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

2 結(jié) 果

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

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

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

2.2 基于卡爾曼濾波器的實時頻譜分析對ripple的偵測

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

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

10組不同小鼠的場電位記錄（每組3段，各5 s時長）數(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模型的實時頻譜分析在偵測場電位ripple段的能量變化時，較STFT方法更為精確.

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

3 討 論

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

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

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

［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

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

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

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