詹躍榮 范影樂(lè)* 楊文偉 楊 勇 李 軼（杭州電子科技大學(xué)自動(dòng)化學(xué)院，杭州 3008）

2（杭州電子科技大學(xué)生物醫(yī)學(xué)工程與儀器研究所，杭州 310018）

引言

動(dòng)作電位（spike）是神經(jīng)元所接收、分析和傳遞的信號(hào)載體，是神經(jīng)元網(wǎng)絡(luò)活動(dòng)的基本表現(xiàn)形式。近年來(lái)，隨著植入式多電極陣列（multi-electrode arrays，MEA）記錄技術(shù)的發(fā)展，用細(xì)胞外神經(jīng)元?jiǎng)幼麟娢贿M(jìn)行記錄和分析已成為神經(jīng)生理學(xué)的重要研究手段。然而，MEA中單個(gè)電極所采集到的信號(hào)往往是電極附近多個(gè)神經(jīng)元發(fā)放的動(dòng)作電位和大量神經(jīng)系統(tǒng)內(nèi)外噪聲的疊加。因此，如何有效地對(duì)spike信號(hào)進(jìn)行檢測(cè)和模式分類，是研究神經(jīng)元?jiǎng)幼麟娢恍蛄校╯pike train）信息特征的重要工作，并可為后續(xù)神經(jīng)元信息編碼研究奠定基礎(chǔ)。

動(dòng)作電位分類過(guò)程主要包括信號(hào)檢測(cè)、特征提取、模式分類三部分。針對(duì)動(dòng)作電位的特征提取，國(guó)內(nèi)外學(xué)者提出了主成分分析法（principal component analysis， PCA）［1］、小 波 分 析（wavelet analysis，WA）［2－3］和基于相空間（phase space）的特征提取［4］等多種分析方法。目前，主要的研究方法為主成分分析法和小波分析的方法。其中，PCA方法通過(guò)樣本協(xié)方差矩陣求取主要特征值，但由于spike信號(hào)具有一定的非線性時(shí)變特性，因此 PCA可能無(wú)法表征spike信號(hào)的完整信息。基于小波分析的方法通過(guò)對(duì)spike信號(hào)進(jìn)行多層分解，在時(shí)頻域上對(duì)spike信號(hào)特征進(jìn)行動(dòng)態(tài)描述。但是，此方法得到的信號(hào)信息是由小波系數(shù)整體合并來(lái)反映的，而在實(shí)際分類中小波系數(shù)是單獨(dú)作用的，很多有用的波形信息可能被丟失。在模式分類階段，通過(guò)對(duì)特征向量的分析處理，利用不同的數(shù)據(jù)特性進(jìn)行統(tǒng)計(jì)聚類。目前，大多數(shù)的神經(jīng)元?jiǎng)幼麟娢环诸惙椒ㄊ腔趧?dòng)作電位的波形差異，包括模板匹配［5］，以及基于神經(jīng)元放電特性的動(dòng)作電位分類方法［6］。由于神經(jīng)元的電生理活動(dòng)存在相互影響，且反映神經(jīng)編碼特性的放電時(shí)間間隔分布較為復(fù)雜，因此使用上述方法具有較大的局限性。

神經(jīng)元電生理活動(dòng)通常被視為處于一個(gè)不穩(wěn)定的混沌狀態(tài)［7］，而描述混沌機(jī)制的一種有效手段就是分形理論；關(guān)聯(lián)維數(shù)是分形維數(shù)的一種，它作為混沌系統(tǒng)的重要特征，反映了系統(tǒng)的復(fù)雜程度，可用于刻畫神經(jīng)元?jiǎng)幼麟娢坏奶卣餍畔?。本研究提出?yīng)用關(guān)聯(lián)維數(shù)，對(duì)動(dòng)作電位的波形特征進(jìn)行無(wú)監(jiān)督自動(dòng)提取。考慮到K均值（K-means）方法具有較好的收斂性，并能夠在一定程度上解決聚類疊加的問(wèn)題，所以采用這種方法進(jìn)行特征聚類，最終實(shí)現(xiàn)動(dòng)作電位的無(wú)監(jiān)督分類。

1 相空間重構(gòu)及其相關(guān)參數(shù)的確定

根據(jù)嵌入理論，單個(gè)時(shí)間序列可以重構(gòu)系統(tǒng)相空間，且重構(gòu)的系統(tǒng)相空間與原始系統(tǒng)拓?fù)涞葍r(jià)。如果嵌入維數(shù)m≥2d+1（d為動(dòng)力系統(tǒng)維數(shù)），則這個(gè)動(dòng)力系統(tǒng)的吸引子空間幾何結(jié)構(gòu)就會(huì)被完全打開，嵌入相空間就可以把有規(guī)律的軌跡恢復(fù)出來(lái)。

設(shè)植入式多電極陣列某通道采集到一個(gè)spike序列，即

式中，xi為采樣值，n為采樣點(diǎn)個(gè)數(shù)。

選擇合適的嵌入延遲時(shí)間τ和嵌入維m，得到一個(gè)m維的嵌入相空間，相空間中的向量可表示為

令N=n－（m－1）τ，則重構(gòu)的多維相空間可表示為

估計(jì)延遲時(shí)間τ可通過(guò)式（4）～式（6）來(lái)確定，有式中，rj=jσ/2，σ為給定時(shí)間序列的標(biāo)準(zhǔn)差，C（m，r，t） 為關(guān)聯(lián)積分。

在嵌入維數(shù)m的求取過(guò)程中，參考偽最近鄰方法（FNN）的思想進(jìn)行計(jì)算，有

式中，xi（m）為m維相空間中的第i個(gè)向量，f（i，m）為第i個(gè)向量的最近鄰點(diǎn)的下標(biāo)。

當(dāng)m大于某個(gè)m0時(shí)，若F（m）不再明顯地發(fā)生變化并接近于1，則此時(shí)m0+1為最小嵌入維數(shù)［9］。

2 關(guān)聯(lián)維數(shù)

關(guān)聯(lián)維數(shù)以相關(guān)性定義，當(dāng)系統(tǒng)狀態(tài)隨時(shí)間變化時(shí)，狀態(tài)變量前后的相關(guān)性可以用來(lái)表征信號(hào)是否有確定的規(guī)律及其程度，所以關(guān)聯(lián)維數(shù)反映了系統(tǒng)的幾何自相似特性。在確定嵌入延遲時(shí)間τ和嵌入維數(shù)m后，在重構(gòu)相空間中利用關(guān)聯(lián)積分計(jì)算關(guān)聯(lián)維數(shù)。

設(shè)Yi是重構(gòu)的相空間中第i個(gè)向量，計(jì)算其余N－1個(gè)向量與Yi的距離，采用最大模表示Euclidean距離，即

rij=d（Yi－Yj）=

定義關(guān)聯(lián)積分，有

3 關(guān)聯(lián)維數(shù)提取神經(jīng)元?jiǎng)幼麟娢惶卣?/h2>

式中，N為相空間代表點(diǎn)（狀態(tài)矢量）的數(shù)目，ε為相空間中給定超小球的半徑，Θ（·）為 Heaviside函數(shù)。

當(dāng)ε充分小時(shí)，關(guān)聯(lián)維數(shù)可定義為

在實(shí)際計(jì)算關(guān)聯(lián)維數(shù)時(shí)，利用 ln（c?（ε））～lnε曲線，通過(guò)線性擬合可得到D2。用高斯函數(shù)u（rij，ε）代替階躍函數(shù)，以改善 ln（c（ε））～lnε關(guān)系的光滑性。高斯函數(shù)定義為

選取一組仿真的典型spike數(shù)據(jù)，如圖1所示。仿真數(shù)據(jù)來(lái)源于加州理工大學(xué)Andersen實(shí)驗(yàn)室提供的C_Easy1_noise01數(shù)據(jù)，其噪聲水平為0.10，采樣頻率為24 kHz，每個(gè)動(dòng)作電位的時(shí)間長(zhǎng)度為8/3 ms。圖2（a）表示 C_Easy1_noise01數(shù)據(jù)中3類不同的spike信號(hào)，每個(gè)子圖為20個(gè)同類spike信號(hào)疊加在一起，這些同類spike信號(hào)在波形上具有類似性，因此狀態(tài)變量的前后相關(guān)性應(yīng)具有類似性。圖2（b）表示對(duì)應(yīng)的3類 spike信號(hào)的 ln（c（ε））～lnε曲線，其中同類 spike信號(hào)的ln（c（ε））～lnε曲線相似度極高，則通過(guò)線性擬合得到的關(guān)聯(lián)維數(shù)也是十分接近的；但不同類型 spike信號(hào)的 ln（c（ε））～lnε曲線存在較大差異，有利于區(qū)分不同類型的 spike信號(hào)。這也表明，利用關(guān)聯(lián)維數(shù)對(duì)動(dòng)作電位模式分類是可行的。

圖1 含噪spike序列Fig.1 Spike train with noise

圖2 3類不同spike信號(hào)及其相應(yīng)的ln（c（ε））～lnε曲線。（a）3類不同spike信號(hào)波形；（b）3類不同spike對(duì)應(yīng)的 ln（c（ε））～lnε曲線Fig.2 Three different types of spike signals and corresponding curve of ln（c（ε））～ lnε.（a）three different types of spike signals；（b）ln（c（ε））～ lnεcurve of three types

在求取關(guān)聯(lián)維數(shù)時(shí)，首先對(duì)采集到的 spike信號(hào)序列進(jìn)行相空間重構(gòu)。以圖2中的3類spike信號(hào)為例，計(jì)算確定嵌入延遲時(shí)間τ=2，嵌入維數(shù)m=12。如圖3所示，當(dāng)嵌入維數(shù)m＞11時(shí)，F(xiàn)（m）不再明顯發(fā)生變化并且接近于1，取最小嵌入維數(shù)m=12。

圖3 嵌入維數(shù)選擇Fig.3 The diagram of choosing embedding dimension

圖4為來(lái)源于上述3類spike信號(hào)的疊加，3類spike信號(hào)總數(shù)為2 990個(gè)。圖5為對(duì)應(yīng)3類 spike信號(hào)的關(guān)聯(lián)維數(shù)分布，可看出3類spike信號(hào)的關(guān)聯(lián)維數(shù)分布具有明顯的分層現(xiàn)象。本研究模式分類算法采用K均值方法，對(duì)計(jì)算得到的spike關(guān)聯(lián)維數(shù)實(shí)現(xiàn)無(wú)監(jiān)督分類。

圖4 3類spike信號(hào)疊加Fig.4 The superposition of three different types of spike

動(dòng)作電位分類算法有下列步驟。

步驟1：采用閾值法，檢測(cè)隱藏在原始信號(hào)中的spike 信號(hào) spi（i=1，2，3，…，n），閾值 Thr確定為［3］

圖5 3類spike信號(hào)關(guān)聯(lián)維數(shù)分布Fig.5 The distribution of correlation dimension of three different types of spike

式中：k為經(jīng)驗(yàn)值，視噪聲而定；E為原始信號(hào)經(jīng)巴特沃斯帶通濾波器（帶寬為300～6 000Hz）后的信號(hào)，median為取中值。

步驟2：確定嵌入延遲時(shí)間τ和合適的嵌入維數(shù)m。

步驟3：計(jì)算每個(gè)spike所對(duì)應(yīng)的關(guān)聯(lián)維數(shù)。

步驟4：利用不同類spike的關(guān)聯(lián)維數(shù)具有明顯的分層分布現(xiàn)象，通過(guò)K均值實(shí)現(xiàn)spike信號(hào)的無(wú)監(jiān)督分類。

4 數(shù)據(jù)說(shuō)明

本算法的驗(yàn)證數(shù)據(jù)分為仿真數(shù)據(jù)和實(shí)驗(yàn)記錄數(shù)據(jù)兩部分。

仿真數(shù)據(jù)來(lái)源于加州理工大學(xué)Andersen實(shí)驗(yàn)室提供的spike含噪仿真數(shù)據(jù)，每段神經(jīng)元spike信號(hào)采樣時(shí)長(zhǎng)為6 000 ms，采樣頻率為24 kHz，每個(gè)spike片段設(shè)置為8/3 ms，即64個(gè)采樣點(diǎn)。仿真數(shù)據(jù)提供每個(gè)spike信號(hào)所屬的神經(jīng)元類別，為本算法的準(zhǔn)確率性能評(píng)價(jià)提供了參照標(biāo)準(zhǔn)。

實(shí)驗(yàn)記錄數(shù)據(jù)來(lái)源于本實(shí)驗(yàn)室采集到的大鼠初級(jí)運(yùn)動(dòng)皮層神經(jīng)元電信號(hào)。通過(guò)手術(shù)將16通道的多電極陣列（MEA）植入大鼠的初級(jí)運(yùn)動(dòng)皮層，然后在術(shù)后第三天采集大鼠該皮層神經(jīng)元的自發(fā)放電序列；MEA采集到的spike模擬信號(hào)經(jīng)前置放大器（PreAmp）初步放大后，進(jìn)入Plexon多通道數(shù)據(jù)采集處理器（MAP），在 MAP提供的40 kHz采樣頻率下轉(zhuǎn)換為spike數(shù)字信號(hào)，并在PC機(jī)SortClient信號(hào)采集軟件配合下，記錄微電極陣列各通道采集的spike序列數(shù)據(jù)。實(shí)驗(yàn)中 MAP系統(tǒng)的增益設(shè)置為1 000，采樣時(shí)長(zhǎng)為8 000 ms，每個(gè) spike片段時(shí)長(zhǎng)為1.4 ms，即54個(gè)采樣點(diǎn)。大鼠初級(jí)運(yùn)動(dòng)皮層神經(jīng)元spike數(shù)據(jù)采集系統(tǒng)如圖6所示。

圖6 大鼠初級(jí)運(yùn)動(dòng)皮層神經(jīng)元spike數(shù)據(jù)采集系統(tǒng)Fig.6 Spike data acquisition system of primary motor cortex in rats

5 結(jié)果及其分析

5.1 仿真數(shù)據(jù)結(jié)果及分析

為了驗(yàn)證關(guān)聯(lián)維數(shù)對(duì)spike信號(hào)分類的普遍有效性，筆者隨機(jī)選擇4組仿真spike信號(hào)進(jìn)行實(shí)驗(yàn)，每組仿真數(shù)據(jù)為3類不同數(shù)目的spike信號(hào)疊加，分類結(jié)果統(tǒng)計(jì)如表1所示。其中，準(zhǔn)確個(gè)數(shù)定義為實(shí)驗(yàn)中每一類正確分類的 spike信號(hào)個(gè)數(shù)，錯(cuò)分個(gè)數(shù)定義為錯(cuò)分到該類但實(shí)屬它類的spike信號(hào)個(gè)數(shù)，準(zhǔn)確率定義為每類正確分類之和與待分類spike信號(hào)的總數(shù)之百分比，1#、2#、3#分別表示每組數(shù)據(jù)中3類不同的spike信號(hào)，M和SD分別定義為每組數(shù)據(jù)3類不同spike信號(hào)分類準(zhǔn)確率的平均值和標(biāo)準(zhǔn)差。從表1中可知，利用關(guān)聯(lián)維數(shù)對(duì)spike信號(hào)模式進(jìn)行分類，準(zhǔn)確率達(dá)到95％以上，每類只有小部分出現(xiàn)了錯(cuò)分。如前所述，由于受微電極陣列技術(shù)的限制，植入式多電極陣列中的各個(gè)電極通常并不能準(zhǔn)確定位到單個(gè)神經(jīng)元。采集的信號(hào)一般來(lái)源于神經(jīng)細(xì)胞間質(zhì)，是電極周圍多個(gè)神經(jīng)元電活動(dòng)的綜合反映；信號(hào)記錄將不可避免地受到外界刺激擾動(dòng)和周圍神經(jīng)細(xì)胞的相互電磁等外干擾的影響，以及細(xì)胞體膜參數(shù)、放電閾值的隨機(jī)起伏變化等因素引起的內(nèi)噪聲影響。因此，對(duì)于部分?jǐn)?shù)據(jù)，關(guān)聯(lián)維數(shù)方法也出現(xiàn)了錯(cuò)分，但從表1可看出關(guān)聯(lián)維數(shù)能有效刻畫不同類spike信號(hào)的特征，并可達(dá)到較為滿意的分類效果。

表1 4組仿真信號(hào)的spike分類結(jié)果統(tǒng)計(jì)Tab.1 Clustering result of four groups of simulated spikes

表2給出了分別采用PCA方法、小波分析方法和基于關(guān)聯(lián)維數(shù)方法的分類結(jié)果對(duì)比，其中PCA方法是利用matlab統(tǒng)計(jì)工具箱里的princomp函數(shù)來(lái)提取 spike的特征信息，而小波分析方法是利用matlab小波分析工具箱里的wavedec函數(shù)進(jìn)行小波分解，進(jìn)而提取spike的特征信息。表2中的1#、2#、3#分別表示每組數(shù)據(jù)中3類不同的spike信號(hào)，M和SD分別定義為每組數(shù)據(jù)3類不同spike信號(hào)分類準(zhǔn)確率的平均值和標(biāo)準(zhǔn)差，其中類間準(zhǔn)確率是每一組數(shù)據(jù)中一類spike聚類正確率。從表2中每組數(shù)據(jù)分類的總準(zhǔn)確率可得出，基于關(guān)聯(lián)維數(shù)的分類方法總體優(yōu)于PCA分類方法。因類間準(zhǔn)確率反映了每類spike信號(hào)的分類效果，根據(jù)表2的類間準(zhǔn)確率可得出，除個(gè)別spike信號(hào)PCA方法的分類效果好于基于關(guān)聯(lián)維數(shù)方法的分類效果外，大部分類間準(zhǔn)確率顯示，基于關(guān)聯(lián)維數(shù)分類方法優(yōu)于PCA分類方法。表2展示針對(duì)不同的數(shù)據(jù)組，小波分析方法和基于關(guān)聯(lián)維數(shù)的分類方法各有長(zhǎng)處，總體分類性能大體相當(dāng)。但小波分析方法針對(duì)不同的數(shù)據(jù)，需要選擇合適的小波基來(lái)分析，如果面對(duì)大數(shù)據(jù)量時(shí)，如何選擇合適的小波基將耗費(fèi)大量時(shí)間，影響分類的效率，因此基于關(guān)聯(lián)維數(shù)的分類方法在可實(shí)現(xiàn)性方面要優(yōu)于小波分析方法。

5.2 實(shí)驗(yàn)記錄數(shù)據(jù)結(jié)果及分析

微電極陣列第11通道采集的大鼠初級(jí)運(yùn)動(dòng)皮層神經(jīng)元?jiǎng)幼麟娢恍蛄袛?shù)據(jù)，如圖7所示。經(jīng)閾值檢測(cè)和濾波后，spike信號(hào)的疊加如圖8所示。在對(duì)spike信號(hào)序列重構(gòu)相空間時(shí)，取嵌入延遲時(shí)間τ=2，取嵌入維數(shù) m=16。

表2 仿真數(shù)據(jù)的不同分類方法的結(jié)果Tab.2 Clustering result of simulated data wi th different sorting techniques

圖7 大鼠初級(jí)運(yùn)動(dòng)皮層神經(jīng)元spike序列Fig.7 Spike train of primary motor cortex in rats

從圖8中可見(jiàn)，雖然采集到的大鼠初級(jí)運(yùn)動(dòng)皮層神經(jīng)元spike信號(hào)的波形比較相似，但仍可以明顯看出存在兩類spike信號(hào)。圖9為兩類spike信號(hào)在相空間中計(jì)算出的關(guān)聯(lián)維數(shù)分布，顯示兩類spike信號(hào)的關(guān)聯(lián)維數(shù)具有明顯的分層分布現(xiàn)象。圖10為經(jīng)K均值分類后的spike信號(hào)疊加，可看出有一部分spike信號(hào)錯(cuò)分，但大部分spike信號(hào)分類是正確的。這說(shuō)明，基于關(guān)聯(lián)維數(shù)的神經(jīng)元?jiǎng)幼麟娢荒Ｊ椒诸惙椒☉?yīng)用于真實(shí)數(shù)據(jù)是有效的。

圖8 大鼠初級(jí)運(yùn)動(dòng)皮層神經(jīng)元spike信號(hào)疊加Fig.8 Spikes of rats’primary motor cortex

圖9 spike信號(hào)關(guān)聯(lián)維數(shù)分布Fig.9 The distribution of spike’s correlation dimension

由于神經(jīng)元?jiǎng)幼麟娢恍盘?hào)是時(shí)變、非平穩(wěn)信號(hào)，傳統(tǒng)的時(shí)域或頻域分析都無(wú)法準(zhǔn)確地表達(dá)神經(jīng)元?jiǎng)幼麟娢坏奶卣餍畔??；陉P(guān)聯(lián)維數(shù)的特征提取優(yōu)勢(shì)在于它刻畫了信號(hào)的非線性特性，通過(guò)將動(dòng)作電位的時(shí)間序列重構(gòu)到高維的相空間中，使得原本在低維空間中無(wú)法區(qū)別的特征信息在高維相空間中可以進(jìn)行有效區(qū)分，增強(qiáng)了非同類動(dòng)作電位的可區(qū)分度。真實(shí)數(shù)據(jù)的實(shí)驗(yàn)結(jié)果表明，本算法框架具有較好的分類效果，代替人工分類具有一定的可行性。

圖10 分類后spike信號(hào)疊加Fig.10 The superposition of classified spike

6 結(jié)語(yǔ)

植入式多電極陣列應(yīng)用中的動(dòng)作電位模式分類技術(shù)，是神經(jīng)編碼研究的前期基礎(chǔ)，因此有著較好的研究意義和應(yīng)用前景。筆者針對(duì)神經(jīng)元?jiǎng)幼麟娢荒Ｊ椒诸?，提出了一種基于關(guān)聯(lián)維數(shù)的特征提取及分類新算法；給出了關(guān)聯(lián)維數(shù)計(jì)算過(guò)程中延遲時(shí)間和嵌入維數(shù)等關(guān)鍵參數(shù)的選擇依據(jù)，以及算法實(shí)現(xiàn)的主要步驟。以仿真實(shí)驗(yàn)數(shù)據(jù)為例，結(jié)合K均值模式聚類方法，實(shí)現(xiàn)了spike信號(hào)的無(wú)監(jiān)督分類。與PCA和小波變換等被廣泛使用的方法進(jìn)行了比較，結(jié)果證明本算法框架對(duì)仿真實(shí)驗(yàn)數(shù)據(jù)模式分類的準(zhǔn)確率達(dá)到了95％以上，表明本算法具有一定的優(yōu)勢(shì)。在實(shí)驗(yàn)室環(huán)境下，進(jìn)行了大鼠初級(jí)運(yùn)動(dòng)皮層神經(jīng)元電活動(dòng)的采集，并將本方法應(yīng)用于真實(shí)的實(shí)驗(yàn)數(shù)據(jù)，取得了良好效果，證實(shí)了非同類神經(jīng)元?jiǎng)幼麟娢坏年P(guān)聯(lián)維數(shù)具有較好的可分性，同時(shí)也驗(yàn)證了本算法的有效性。

［1］AdamosDA， Kosmidis EK， TheophilidisG. Performance evaluation of PCA-based spike sorting algorithms［J］.Comput Methods Prog Biol，2008，91（3）：232 －244.

［2］Geng Xinling，Hu Guangshu，Tian Xin.Neural spike sorting using mathematical morphology，multiwavelets transform and hierarchical clustering［J］.Neurocomputing，2010，73：707 －715.

［3］Quiroga RQ，Nadasdy Z，Ben-Shaul Y.Unsupervised spike detection and sorting with wavelets and superparamagnetic clustering［J］.Neural Comput，2004，16：1661 － 1687.

［4］Chibirovaa OK， Aksenovaa TI， Benabida A， et al.Unsupervised spike sorting of extracelluar electrophysiological recording in subthalamic nucleus of Parkinsonian patients［J］.Biosystems，2005，79：159－171.

［5］Vargas-Irwin C，Donoghue JP.Automated spike sorting using density grid contour clustering and subtractive waveform decomposition［J］.Journal of Neuroscience Methods，2007，164（1）：1－18.

［6］Delesclue M，Pouzat C.Efficient spike-sorting of multi-state neurons using inter-spike intervals information［J］.Journal of Neuroscience Methods，2006，150：16－29.

［7］彭建華，劉延柱.腦科學(xué)中若干非線性動(dòng)力學(xué)問(wèn)題［J］.力學(xué)進(jìn)展，2003，33（3）：325－321.

［8］Kim HS，Eykholt R，Salas JD.Nonlinear dynamics，delay times and embedding windows［J］.Physica D（S0167 － 2789），1999，127：48－60.

［9］Cao Liangyue.Practical method of determining the minimum embedding dimension of a scalar time series［J］.Physica D，1997，110：43－50.