萬 紅 張 超 劉新玉 尚志剛

(鄭州大學(xué)電氣工程學(xué)院，鄭州 450001)

波形特征的高斯混合模型鋒電位分類算法

萬 紅*張 超 劉新玉 尚志剛

(鄭州大學(xué)電氣工程學(xué)院，鄭州 450001)

鋒電位分類是進(jìn)行大腦信息處理機(jī)制研究的基本步驟之一。針對(duì)鋒電位信號(hào)的復(fù)雜性和非平穩(wěn)性，從統(tǒng)計(jì)聚類的角度出發(fā)，采用高斯混合模型描述鋒電位信號(hào)的概率密度函數(shù)，提出了一種新的基于波形變化率特征和高斯混合模型聚類的鋒電位分類方法。首先計(jì)算鋒電位的波形變化率，然后利用最大差異方法獲得鋒電位波形的低維特征，最后采用高斯混合模型算法對(duì)特征進(jìn)行聚類，實(shí)現(xiàn)鋒電位分類。采用開放的仿真數(shù)據(jù)分析了該算法的分類精度和可行性，然后分別利用來自5只大鼠和1只恒河猴初級(jí)視覺皮層的實(shí)測(cè)數(shù)據(jù)驗(yàn)證了該算法的實(shí)用性，并與主成分分析特征的高斯混合模型聚類、幅值特征的高斯混合模型聚類和變化率特征的k均值聚類等3種方法進(jìn)行了比較。仿真實(shí)驗(yàn)中，在噪聲水平為0.05、0.10、0.15、0.20時(shí)，誤分率分別為1.18%±1.18%、1.41%±1.06%、2.27%±1.51%、2.98%±2.06%, 低于其他3種方法；實(shí)測(cè)實(shí)驗(yàn)中，恒河猴數(shù)據(jù)的J3準(zhǔn)則值為13.50±5.26，大鼠數(shù)據(jù)的J3準(zhǔn)則值為26.43±10.46。與其他3種方法相比，平均J3準(zhǔn)則值較大，且顯著高于幅值特征的高斯混合模型聚類算法。所提出的方法表現(xiàn)出較高的分類精度和較好的類可分性，為實(shí)現(xiàn)神經(jīng)元鋒電位的可靠分類提供了一種有效的手段。

鋒電位分類；波形變化率；最大差異算法；高斯混合模型

引言

在神經(jīng)科學(xué)中，鋒電位(spike)是神經(jīng)信息接收、編碼和傳遞的載體，具有重要的研究意義。在微電極陣列胞外記錄方式中，每根微電極尖端可能被多個(gè)神經(jīng)元包圍，即一個(gè)電極可能同時(shí)記錄到多個(gè)神經(jīng)元發(fā)放的鋒電位。因此需要首先把微電極陣列記錄到的不同神經(jīng)元的放電活動(dòng)區(qū)分開來，即實(shí)現(xiàn)spike分類。然而在神經(jīng)電生理實(shí)驗(yàn)中，神經(jīng)信號(hào)會(huì)受到背景噪聲的污染，使spike波形發(fā)生改變，增加了spike分類的難度。尋求簡單高效的特征提取方法，進(jìn)而實(shí)現(xiàn)精確的spike分類是目前神經(jīng)信息編碼機(jī)制研究一直致力解決的問題。

Spike分類一般包括帶通濾波、檢測(cè)、波形對(duì)齊、特征提取、降維和聚類等階段。其中，特征提取是spike分類中至關(guān)重要的一步，特征提取的質(zhì)量對(duì)分類結(jié)果有決定性作用。經(jīng)典的特征提取方法主要提取spike的波形特征，如幅值、波寬、能量或相位[1-2]等，這種方法計(jì)算簡單、效率較高，尤其對(duì)于波形差異較大的spike分類效果較好，但是對(duì)于波形比較接近的spike則不甚理想；此外，主成分分析(principal component analysis，PCA)[3-4]、小波變換[5-6]等特征提取方法應(yīng)用也比較廣泛，雖然這些方法提高了特征的分類性能，但是PCA難以完整刻畫出spike波形，而且小波變換對(duì)小波基的選擇并沒有有效的準(zhǔn)則。

此外，特征聚類也是spike分類中的關(guān)鍵步驟，目前常用的方法有k均值聚類[7]、超順磁聚類[8]等，雖然這些方法都取得了較好的聚類效果，但是k均值需要首先根據(jù)初始聚類中心確定一個(gè)初始劃分，而且這個(gè)初始聚類中心對(duì)聚類結(jié)果有很大影響；超順磁聚類雖然是無監(jiān)督聚類方法，但計(jì)算繁瑣、速度較慢。封洲燕等利用PCA結(jié)合最大最小距離聚類的spike分類方法，在大鼠海馬區(qū)spike信號(hào)的分類中取得了較好的分類效果[9]。但是這種方法只適用于多通道的spike分類。

針對(duì)以上問題，本研究提出一種新的spike分類方法，即首先計(jì)算spike波形變化率，并基于spike的偶極子模型從理論上探討了波形變化率在描述spike特征時(shí)的優(yōu)勢(shì)，然后利用最大差異方法獲得spike的低維特征，最后采用高斯混合模型實(shí)現(xiàn)特征聚類。波形變化率，即一階導(dǎo)數(shù)，能夠有效度量不同神經(jīng)元spike波形之間的差異，而最大差異算法是一種自適應(yīng)的特征提取方法，可以根據(jù)spike波形自動(dòng)尋找到差異最顯著的特征。因此，本研究提出的spike分類方法可以根據(jù)神經(jīng)信號(hào)的不同自適應(yīng)的提取出具有顯著差異的spike波形特征，提高了spike的分類性能。實(shí)驗(yàn)中利用仿真數(shù)據(jù)和大鼠初級(jí)視覺皮層(又稱V1區(qū))和恒河猴V1區(qū)神經(jīng)信號(hào)對(duì)算法進(jìn)行了驗(yàn)證。

1 材料與方法

1.1 鋒電位分類方法

假設(shè)一次實(shí)驗(yàn)(即trial)中，對(duì)于給定通道，微電極共記錄到了N個(gè)spike，分別來自于C個(gè)神經(jīng)元，而且每個(gè)spike有M個(gè)采樣點(diǎn)組成，即Xi={xi,1,xi,2,…,xi,M}(i=1,2,…,N)，則這N個(gè)spike可以表示為

(1)

式中，Z為N×M矩陣，表示微電極陣列在一次trial中記錄到的所有spike波形的幅值。

Spike分類就是利用模式識(shí)別方法將這些無標(biāo)簽的spike信號(hào)分為不同類別，以實(shí)現(xiàn)對(duì)單個(gè)神經(jīng)元活動(dòng)特性進(jìn)行研究的目的。Spike分類過程如圖1所示。

圖1 Spike分類過程Fig.1 Spike sorting procedure

1.1.1 鋒電位特征提取

特征提取是spike分類中的關(guān)鍵步驟，特征提取的質(zhì)量直接關(guān)系到分類的性能。從理論上講，同一神經(jīng)元發(fā)放的鋒電位波形相同，而不同神經(jīng)元的鋒電位在波形上存在一定差異，然而，在信號(hào)采集過程中，背景噪聲模糊了這種差異。變化率，又稱陡度[7]，表示幅值隨時(shí)間變化的改變程度，可以表征波形變化的細(xì)微差別，是波的一個(gè)重要特征。

本研究使用spike變化率特征的靈感最初是來自于spike的偶極子模型[7]，偶極子p在距偶極子中心r處產(chǎn)生的電勢(shì)可表示為

(2)

式中，p為偶極子算子，r為電極與神經(jīng)元之間的距離，θ為向量p與向量r之間的夾角，ε0為電介質(zhì)常數(shù)。

偶極子是由細(xì)胞外的電流形成的，該電流是由Na+從細(xì)胞外空間流向細(xì)胞內(nèi)空間，以平衡細(xì)胞膜而形成的。

從式(2)可以看出，距離神經(jīng)元遠(yuǎn)近不同的兩個(gè)電極，假如第1個(gè)電極與神經(jīng)元間的距離小于第2個(gè)電極與神經(jīng)元間的距離。在這種情況下，從這兩個(gè)電極記錄到的spike會(huì)表現(xiàn)出相同的形狀但幅值不同。這是因?yàn)閟pike幅值反比例于r2，由于兩個(gè)記錄到的spike具有不同的幅值，所以spike幅值可以作為區(qū)分兩個(gè)不同神經(jīng)元spike的特征。而變化率是波形的一階導(dǎo)數(shù)，刻畫的是波形在不同時(shí)刻的波形變化趨勢(shì)，能明顯反映出幅值變化急劇程度，是波形變化的顯著數(shù)字特征。因此，利用spike波形的變化率作為分類特征，并利用最大差異法[10]對(duì)特征進(jìn)行降維。

(3)

式中，i=1,2,…,M-1，Zk,i表示矩陣Z中的第k行、第i列的元素值，Δt為采樣間隔。

(4)

式中，Ld為局部差異向量。

1.1.2 高斯混合模型分類算法

假設(shè)spike數(shù)據(jù)集X={x1,x2,…,xN}服從由K個(gè)高斯成分(即K類)組成的高斯混合模型(gaussianmixturemodel，GMM)。對(duì)于單個(gè)spike樣本xi，其高斯混合分布的概率密度函數(shù)為

(5)

(6)

GMM的訓(xùn)練過程需要利用最大期望(expectationmaximization，EM)算法進(jìn)行非線性概率函數(shù)的優(yōu)化，首先對(duì)確定的分類數(shù)目進(jìn)行模型初始化，并給出θk的初始估計(jì)值θ0k。待求出參數(shù)λ=(πk,θk)后，就可以得到每個(gè)spike屬于每一類的概率，根據(jù)最大后驗(yàn)概率準(zhǔn)則，將該spike判為概率最大的那一類。EM算法的計(jì)算過程如下：

EM算法有兩步，即期望歩(E-step)和最大值歩(M-step)，通過更新參數(shù)πk，μk，σk在E步和M步之間進(jìn)行迭代，使似然函數(shù)達(dá)到最大。

1)E-step：通過似然函數(shù)計(jì)算當(dāng)前估計(jì)和上次估計(jì)的差值Q(p,pold)[5]，表示為

(7)

式中，N和K分別表示spike總數(shù)和混合模型中高斯成分的個(gè)數(shù)。

2)M-step：使式(7)中的Q(p,pold)達(dá)到最大，得到新的參數(shù)πk，new、μk，new、σk，new，更新過程可由下列式子得出

(8)

這些新的參數(shù)可以在又一次的迭代中得到估計(jì)和修正，直到使Q(p,pold)達(dá)到預(yù)定的收斂條件，至此就可以獲得參數(shù)λ，完成spike的分類。

1.2 評(píng)價(jià)指標(biāo)

1.2.1 誤分率

對(duì)spike分類而言，由于分類而造成的誤分樣本的個(gè)數(shù)是一個(gè)重要的評(píng)價(jià)指標(biāo)。因此，使用誤分率來定量衡量分類算法的性能。假設(shè)待分?jǐn)?shù)據(jù)集X={x1,x2,…,xN}可分為K類，每類的誤分樣本數(shù)為wi(i=1,…,K)，則誤分率R可表示為[11]

(9)

1.2.2J3準(zhǔn)則值

對(duì)于spike分類而言，從模式識(shí)別的角度看，提取的特征應(yīng)當(dāng)滿足兩個(gè)條件：一是要使每一類的特征值緊密的分布在某個(gè)值附近，即類內(nèi)方差要??；二是不同類之間的特征值分布的距離要盡可能遠(yuǎn)，即類間距離要大。由類可分性準(zhǔn)則可知，特征的類可分性越好，其J3準(zhǔn)則值就會(huì)越大。J3準(zhǔn)則表征了特征的類可分性，定義為[12]

J3=tr{S-1Sm}

(10)

1.3 數(shù)據(jù)采集

為了驗(yàn)證所提出方法的有效性和實(shí)用性，采用仿真數(shù)據(jù)研究了該算法的可行性和分類精度，用實(shí)測(cè)數(shù)據(jù)驗(yàn)證了實(shí)際分類效果，并與多種分類算法進(jìn)行了比較。關(guān)于仿真數(shù)據(jù)和實(shí)測(cè)數(shù)據(jù)的來源分別簡述如下：

為了說明該算法的有效性，仿真數(shù)據(jù)選用了公開的spike分類測(cè)試數(shù)據(jù)，來源于文獻(xiàn)[13]。該數(shù)據(jù)共包含有4組仿真信號(hào)，每組信號(hào)中由噪聲水平分別為0.05、0.10、0.15、0.20的4段數(shù)據(jù)組成, 噪聲水平定義為信號(hào)背景噪聲的標(biāo)準(zhǔn)差。這4組數(shù)據(jù)分別記為data1、data2、data3和data4。每段數(shù)據(jù)時(shí)間長度為60 s，采樣頻率為24 kHz，共包含有來自3個(gè)神經(jīng)元發(fā)放的大約3 400個(gè)spike，每個(gè)spike由64個(gè)采樣點(diǎn)組成。

為了驗(yàn)證該算法的實(shí)用性，采用的實(shí)測(cè)數(shù)據(jù)包含兩部分，實(shí)測(cè)數(shù)據(jù)1為從CRCNS.org神經(jīng)數(shù)據(jù)共享平臺(tái)—恒河猴V1區(qū)采集得到，選擇了其中5組典型的數(shù)據(jù)，分別記為data5、data6、data7、data8和data9；實(shí)測(cè)數(shù)據(jù)2為本課題組從大鼠V1區(qū)采集得到，選用其中5組來自不同大鼠且具有代表性的數(shù)據(jù)，分別記為data10、data11、data12、data13和data14。

實(shí)測(cè)數(shù)據(jù)1：來源于Collaborative Research in Computational Neuroscience-Data sharing(CRCNS)的數(shù)據(jù)共享平臺(tái)[14]。實(shí)驗(yàn)動(dòng)物為體重為4～5 kg的恒河猴1只，利用氯胺酮麻醉后在V1區(qū)植入8×8微電極陣列(Neuronexus，電極間距為200 μm，電極直徑30 μm)采集得到，詳細(xì)實(shí)驗(yàn)過程請(qǐng)參考文獻(xiàn)[15]。所選用的5組數(shù)據(jù)分別包含有32 042、54 749、41 097、31 151和28 116個(gè)spike，采樣頻率為24.4 kHz，每個(gè)spike由12個(gè)采樣點(diǎn)組成。

實(shí)測(cè)數(shù)據(jù)2：實(shí)驗(yàn)動(dòng)物為體重200～250 g的成年Long Evans(LE)大鼠5只，雌雄不拘。利用10%水合氯醛麻醉后在V1區(qū)植入2×8微電極陣列(Microprobe，鉑銥合金，電極間距250 μm，電極直徑125 μm，尖端直徑15 μm，電極阻抗在0.5～1.0 MΩ之間)通過數(shù)據(jù)采集系統(tǒng)(CerberusTM，Blackrock Inc.)采集得到，詳細(xì)實(shí)驗(yàn)過程請(qǐng)參考文獻(xiàn)[16]。所選用的5組數(shù)據(jù)分別包含有747、1 329、810、1 130和1 417個(gè)spike，采樣頻率為30 kHz，每個(gè)spike由44個(gè)采樣點(diǎn)組成。

2 結(jié)果

圖2 3個(gè)典型的spike模板及其變化率。(a)～(c)3個(gè)典型的spike模板波形；(d)～(f)對(duì)應(yīng)的變化率Fig.2 Three spike templates and its waveform changing rate. (a)～(c) Waveforms of the three spike templates；(d)～(f) Waveforms changing rate of the three spike templates

圖3 3種不同方法對(duì)相同數(shù)據(jù)特征提取的統(tǒng)計(jì)分布對(duì)比。(a)波形變化率最大差異法；(b)PCA方法；(c)原始波形最大差異法Fig.3 Statistical distribution of features based on three extraction methods for the same data. (a) The maximum-difference of waveform changing rate method; (b) The PCA method; (c) The maximum-difference of original waveform method

圖2為3個(gè)不同的spike模板波形圖及其變化率。由圖可知，3個(gè)spike模板的變化率差異要比原始信號(hào)的波形差異更顯著。圖3為波形變化率最大差異法、PCA法和原始波形最大差異法3種算法提取的第一維特征值統(tǒng)計(jì)分布，其中的粗實(shí)線為高斯曲線擬合結(jié)果。由圖可知，除PCA法呈單峰分布外，本方法和原始波形最大差異法所提取特征的分布均呈雙峰形狀，但本方法提取的特征兩峰之間距離較大，而且沒有重疊。由此可知本方法提取的特征要優(yōu)于PCA和原始波形最大差異特征提取方法。

2.1 仿真結(jié)果

圖4為其中一組仿真數(shù)據(jù)的聚類結(jié)果。圖中(a)是4個(gè)不同噪聲水平(分別是0.05、0.10、0.15、0.20)的仿真數(shù)據(jù)波形，從中可以看出，隨著噪聲水平的增加，三類spike之間的差異越來越模糊；(b)～(e)分別為4種不同聚類效果方法的結(jié)果，即變化率特征的高斯混合模型聚類(變化率特征+GMM)、PCA特征的高斯混合模型聚類(PCA+GMM)、幅值特征的高斯混合模型聚類(幅值特征+GMM)、變化率特征的k均值聚類(變化率特征+k均值)。由圖可知，隨著噪聲水平的增加，4種分類方法中不同spike類之間的重疊程度都逐漸增加，差異度變小，但是與其他3種方法相比，本方法類間距變化較小，依然保持著較好的分類效果，有著良好的分類性能。

為了進(jìn)一步度量本算法的分類性能，表1給出了4種不同分類方法對(duì)4個(gè)仿真數(shù)據(jù)集的誤分率對(duì)比結(jié)果。從表中可以看出，隨著噪聲水平的增加，4種方法的誤分率都呈逐漸升高趨勢(shì)，但與其他3種方法相比，在不同噪聲水平時(shí)本方法的誤分率都是最低的，分別為1.18%、1.41%、2.27%和2.98%。與圖4的分類結(jié)果一致。綜合圖4和表1的結(jié)果可知，從特征提取的角度看，與PCA特征提取方法相比，最大差異特征提取方法能夠有效選擇出差異最顯著的特征，而且與原始波形相比，變化率受噪聲的影響相對(duì)較小，具有一定的抗噪能力；從聚類的角度看，GMM算法對(duì)spike波形特征的非高斯分布可以進(jìn)行有效擬合，聚類效果也優(yōu)于經(jīng)典的k均值算法。

圖4 不同特征提取方法分類結(jié)果對(duì)比(圖中的不同顏色代表不同類別的鋒電位特征值分布，最左側(cè)圖中數(shù)字表示噪聲水平)。(a)Spike波形；(b)變化率特征+GMM算法；(c)PCA+GMM算法；(d)幅值特征+GMM算法；(e)變化率特征+k均值算法Fig.4 Results for different spike sorting methods from simulated data(The different colors denote the features of spike from the different clusters and the number in left figures denote noise level). (a) Spike waveforms; (b) Spike sorting method based on changing rate features and GMM cluster; (c) Spike sorting method based on PCA features and GMM cluster; (d) Spike sorting method based on waveform features and GMM cluster; (e) Spike sorting method based on waveform changing rate features and k-means cluster

表1 仿真數(shù)據(jù)的不同方法誤分率對(duì)比表Tab.1 Comparison of misclassification rate of different methods in simulated data

注：平均誤分率用(均值±標(biāo)準(zhǔn)差)表示(n=4)；*表示P<0.05，**表示P<0.01。

Note: The average misclassification rate is presented by mean ± standard deviation (n=4),*denotesP<0.05 and**denotesP<0.01.

2.2 實(shí)測(cè)結(jié)果

實(shí)測(cè)數(shù)據(jù)的分類結(jié)果如圖5所示，其中前兩組數(shù)據(jù)(data5和data6)來自于恒河猴V1區(qū)，后兩組數(shù)據(jù)(data10和data11)來自于大鼠V1區(qū)。圖中(a)是實(shí)測(cè)數(shù)據(jù)的spike波形；(b)為本方法的分類結(jié)果，(c)～(e)是其他3種不同方法的分類結(jié)果，與仿真數(shù)據(jù)的分類方法一致。從圖(b)和(e)可以看出，當(dāng)特征提取方法相同時(shí)，分類結(jié)果相似，方框標(biāo)出了兩種方法的分類結(jié)果不同之處。但是，從分類結(jié)果來看，本方法的聚類結(jié)果要明顯好于其他方法，類間距離明顯，類與類之間沒有重疊，與仿真數(shù)據(jù)的結(jié)果一致。

圖5 實(shí)測(cè)數(shù)據(jù)分類結(jié)果對(duì)比(圖中的不同顏色代表不同類別的鋒電位特征值分布)。(a)Spike波形；(b)變化率特征+GMM算法；(c)PCA+GMM算法；(d)幅值特征+GMM算法；(e)變化率特征+k均值算法Fig.5 Results for different spike sorting methods from real data (The different colors denote the features of spike from the different clusters). (a) Spike waveforms; (b) Spike sorting method based on changing rate features and GMM cluster; (c) Spike sorting method based on PCA features and GMM cluster; (d) Spike sorting method based on amplitude features and GMM cluster; (e) Spike sorting method based on changing rate features and k-means cluster

由于實(shí)測(cè)數(shù)據(jù)中每個(gè)spike信號(hào)的類別是未知的，難以用上述的誤分率進(jìn)行衡量，為了能定量度量本算法在實(shí)測(cè)數(shù)據(jù)中的分類性能，此處選擇類可分性度量指標(biāo)，J3準(zhǔn)則值，對(duì)上述分類算法進(jìn)行定量評(píng)估。表2給出了實(shí)測(cè)數(shù)據(jù)1和實(shí)測(cè)數(shù)據(jù)2的分類J3準(zhǔn)則值。恒河猴的數(shù)據(jù)結(jié)果顯示(實(shí)測(cè)數(shù)據(jù)1)，本算法的類可分性顯著優(yōu)于幅值特征+GMM算法(P<0.05)。雖然與PCA+GMM和變化率特征+k均值算法相比，并沒有顯著統(tǒng)計(jì)學(xué)差異，但是在平均水平上，本算法J3準(zhǔn)則值高于這兩種算法；大鼠的數(shù)據(jù)結(jié)果顯示(實(shí)測(cè)數(shù)據(jù)2)，本算法顯著優(yōu)于PCA+GMM(P<0.05)和幅值特征+GMM(P<0.01)，而且平均J3準(zhǔn)則值也高于變化率特征+k均值算法。綜上可知，無論是恒河猴數(shù)據(jù)還是大鼠數(shù)據(jù)，本算法的平均J3準(zhǔn)則值都大于其他3種算法，表現(xiàn)出了較好的類可分性。

表2 實(shí)測(cè)數(shù)據(jù)不同方法的J3準(zhǔn)則值對(duì)比表Tab.2 Comparison of J3 criterion values of different methods in real data

注：J3平均值用(均值±標(biāo)準(zhǔn)差)表示(n=5)；*表示與本算法相比P<0.05，**表示與本算法相比P<0.01。

Note: The averageJ3value is presented by mean ± standard deviation (n=5),*denotesP<0.05 compared with the proposed algorithm,**denotesP<0.01 compared with the proposed algorithm.

3 討論

Spike分類被認(rèn)為是解析大腦信息處理機(jī)制的第1步[17]。Spike的產(chǎn)生機(jī)制決定了波形變化對(duì)于spike分類的重要性，而偶極子模型描述了影響spike波形變化的因子。本研究針對(duì)波形特征給出了一種新的spike分類方法，即首先提取出spike波形的變化率特征，然后利用最大差異法對(duì)波形變化率特征進(jìn)行降維，最后利用高斯混合模型聚類算法對(duì)特征進(jìn)行聚類，實(shí)現(xiàn)spike分類。算法驗(yàn)證方面采用來自公開的spike分類仿真測(cè)試數(shù)據(jù)和來自不同動(dòng)物種類V1區(qū)的實(shí)測(cè)數(shù)據(jù)，保證了實(shí)驗(yàn)結(jié)果的客觀性和通用性。

特征提取和聚類對(duì)于spike分類而言是兩個(gè)至關(guān)重要的步驟。在特征提取上，通過前期大量的實(shí)驗(yàn)數(shù)據(jù)分析，本研究提出了一種新的spike分類特征提取方法，即基于spike的波形變化率，利用最大差異算法獲得差異最顯著的spike分類特征。本研究提出的spike分類特征提取方法最大的優(yōu)勢(shì)在于其是一種自適應(yīng)的特征提取方法，可以根據(jù)spike波形自適應(yīng)的提取出差異最顯著的特征，而且對(duì)于spike波形相近的神經(jīng)元，也具有較好的分類效果。在聚類上，由于神經(jīng)元活動(dòng)的復(fù)雜性和非平穩(wěn)性，以及特征選擇的自適應(yīng)性，在spike特征分布上往往表現(xiàn)出非高斯特性，因此，選用GMM算法對(duì)上述spike分類特征進(jìn)行聚類，GMM分類算法是一種統(tǒng)計(jì)聚類方法，可以有效表征spike特征的非高斯分布。

在實(shí)驗(yàn)分析上，本研究選用了公開的spike分類仿真數(shù)據(jù)和多種實(shí)測(cè)數(shù)據(jù)，并與目前常用的spike分類算法進(jìn)行了對(duì)比。由仿真數(shù)據(jù)結(jié)果可知，本方法的平均誤分率相對(duì)較低，對(duì)于不同的噪聲水平，即0.05、0.10、0.15、0.20，其spike誤分率分別為1.18%±1.18%、1.41%±1.06%、2.27%±1.51%、2.98%±2.06%(見表1)，表明本算法具有較高的分類精度；從對(duì)實(shí)測(cè)數(shù)據(jù)結(jié)果來看，本方法的平均J3準(zhǔn)則值是4種方法中最大的，恒河猴數(shù)據(jù)為13.50±5.26，大鼠數(shù)據(jù)為26.43±10.46(如表2)，說明本算法得到的結(jié)果具有最大類可分性。仿真和實(shí)測(cè)數(shù)據(jù)的結(jié)果具有一致性，表明了所提出方法的有效性。

綜上可知，不論從特征提取的角度，還是從聚類的角度，本方法都體現(xiàn)了較好的分類性能。波形變化率能夠有效度量不同神經(jīng)元spike波形之間的差異，GMM算法對(duì)spike波形特征的非高斯分布可以進(jìn)行有效擬合，聚類效果也優(yōu)于經(jīng)典的聚類算法。但是，也不可否認(rèn)的是，與原始波形相比，波形變化率可能會(huì)放大噪聲，不過在嵌入式高放大倍數(shù)的spike信號(hào)采集方式下，加之特征提取前spike經(jīng)過了濾波、閾值檢測(cè)、波形對(duì)齊等步驟進(jìn)行預(yù)處理，已將噪聲干擾減小到最小，信噪比相對(duì)較高，故一般來說這個(gè)問題對(duì)spike分類而言影響較小。

4 結(jié)論

本研究提出了一種基于波形變化率特征和高斯混合模型聚類的spike分類方法，即首先提取出spike波形的變化率特征，然后利用最大差異法對(duì)波形變化率特征進(jìn)行降維，最后利用高斯混合模型聚類算法對(duì)特征進(jìn)行聚類，實(shí)現(xiàn)spike分類。并利用仿真信號(hào)和來自大鼠V1區(qū)及恒河猴V1區(qū)神經(jīng)信號(hào)進(jìn)行了驗(yàn)證。仿真和實(shí)測(cè)數(shù)據(jù)結(jié)果表明，本方法不僅可以自適應(yīng)提取出具有顯著差異的spike特征，而且具有較高的分類精度和較好的類可分性，為基于spike信號(hào)的神經(jīng)信息編碼和解碼研究奠定了基礎(chǔ)。

[1] Lewicki MS. A review of methods for spike sorting: the detection and classification of neural action potentials [J]. Network: Comp Neural, 1998, 9(4): R53-R78.

[2] Kim KH, Kim SJ. Neural spike sorting under nearly 0-dB signal-to-noise ratio using nonlinear energy operator and artificial neural-network classifier [J]. IEEE Trans Biomed Eng, 2000, 47(10): 1406-1411.

[3] Adamos DA, Kosmidis EK, Theophilidis G. Performance evaluation of PCA-based spike sorting algorithms [J]. Comput Meth Prog Bio, 2008, 91(3): 232-244.

[4] 丁偉東，袁景淇，梁培基. 多電極鋒電位信號(hào)檢測(cè)與分類方法研究 [J]. 儀器儀表學(xué)報(bào), 2006, 27(12): 1636-1640.

[5] Bestel R, Daus AW, Thielemann C. A novel automated spike sorting algorithm with adaptable feature extraction [J]. J Neurosci Methods, 2012, 211(1): 168-178.

[6] 丁穎，范影樂，楊勇. 基于小波特征與動(dòng)態(tài)高斯混合模型的動(dòng)作電位分類研究 [J]. 儀器儀表學(xué)報(bào), 2011, 32(2): 475-480.

[7] Paraskevopoulou SE, Barsakcioglu DY. Feature extraction using first and second derivative extrema (FSDE) for real-time and hardware-efficient spike sorting [J]. J Neurosci Methods, 2013, 215(1): 29-37.

[8] Wild J, Prekopcsak Z. Performance comparison of extracellular spike sorting algorithms for single-channel recordings [J]. J Neurosci Methods, 2012, 203(2): 369-376.

[9] 王靜，封洲燕. 多通道神經(jīng)元鋒電位檢測(cè)和分類的新方法 [J]. 生物化學(xué)與生物物理進(jìn)展, 2009, 36(5): 641-647.

[10] Gibson S, Judy JW. Technology-aware algorithm design for neural spike detection, feature extraction, and dimensionality reduction [J]. IEEE Trans Neur Sys Reh, 2010, 18(5): 469-478.

[11] Geng Xinling, Hu Guangshu. Unsupervised feature selection by kernel density estimation in wavelet-based spike sorting [J]. Biomed Signal Proces, 2012, 7(2): 112-117.

[12] Theodoridis S, Koutroumbas K, 著. 李晶皎, 王愛俠, 王驕，等譯. 模式識(shí)別 [M]. 北京: 電子工業(yè)出版社, 2010: 188-190.

[13] Martinez J, Pedreira C, Ison MJ, et al. Realistic simulation of extracellular recordings [J]. J Neurosci Methods, 2009, 184(2): 285-293.

[14] Chu CCJ, Chien PF, Hung CP. Multi-electrode recordings of ongoing activity and responses to parametric stimuli in macaque V1 [DB/OL]. CRCNS.org, http://dx.doi.org/10.6080/K0J1012K, 2014-03-06/2014- 10-10.

[15] Chu CCJ, Chien PF, Hung CP. Tuning dissimilarity explains short distance decline of spontaneous spike correlation in macaque V1 [J]. Vision Res, 2014, 96(1): 113-132.

[16] 萬紅，張曉娜，劉新玉，等. 基于局部場(chǎng)電位固有模態(tài)分量的響應(yīng)調(diào)諧特性研究 [J]. 中國生物醫(yī)學(xué)工程學(xué)報(bào), 2013, 32(3): 292-298.

[17] Gibson S, Judy JW, Markovic D. Spike sorting: the first step in decoding the brain [J]. IEEE Signal Proc Mag, 2012, 29(1): 124-14.

Spike Sorting Using Gaussian Mixture Model of Waveform Feature

Wan Hong*Zhang Chao Liu Xinyu Shang Zhigang

(SchoolofElectricalEngineering,ZhengzhouUniversity,Zhengzhou450001,China)

Spike sorting is one of basic steps to study brain information processing mechanism. Regarding to the complexity and non-stationary of spike signals, a new spike sorting method based on waveform changing rate features and Gaussian mixture model (GMM) cluster was proposed in this work, in which the GMM described the probability density function of spike features from statistical clustering viewpoint. In the proposed method, firstly, the changing rate of spike waveforms was calculated, then low-dimensional features were obtained using maximum-difference method, and the features were finally sorted by the GMM. The accuracy and feasibility of the algorithm was measured by the simulated data from the published database, and the practicality was validated by the real data from primary visual cortex of five rats and one macaque monkey. And its performance was compared to other three sorting methods that based on principal component features and GMM cluster, based on waveform features and GMM cluster as well as based on waveform changing rate features andk-meancluster. For the simulated data, the misclassification rate of proposed method below the other three methods, that is 1.18%±1.18%, 1.41%±1.06%, 2.27%±1.51% and 2.98%±2.06% when the noise level is 0.05, 0.10, 0.15 and 0.20, respectively. For the real data, theJ3value of proposed method is 13.50±5.26 with monkey data and 13.50±5.26 with rat data. Compared with the other three methods, the proposed method gives the maximumJ3value and is higher significantly than the sorting method based on waveform features and GMM cluster. It has higher precision and classification performance as well as provides an effective approach to achieve reliable spike sorting.

spike sorting; waveform changing rate; maximum-difference algorithm; gaussian mixture model

10.3969/j.issn.0258-8021. 2016. 04.003

2014-10-28， 錄用日期:2016-05-10

國家自然科學(xué)基金(U1304602)；河南省重點(diǎn)科技攻關(guān)項(xiàng)目(122102210102, 162102310167)

R318.04

A

0258-8021(2016) 04-0402-09

*通信作者(Corresponding author)， E-mail: wanhong@zzu.edu.cn