王 潔，郭天翔，盧云山，趙 冰，熊 鵬，杜海曼

1.河北大學(xué) 電子信息工程學(xué)院，河北 保定 071000

2.康泰醫(yī)學(xué)系統(tǒng)（秦皇島）股份有限公司，河北 秦皇島 066004

現(xiàn)代多模態(tài)顱內(nèi)記錄設(shè)備逐漸具備植入式腦機(jī)接口系統(tǒng)高時(shí)空分辨率的要求[1]，測(cè)量細(xì)胞活動(dòng)的最流行技術(shù)之一是使用細(xì)胞外電極陣列。現(xiàn)有多電極技術(shù)已經(jīng)能夠同時(shí)記錄數(shù)千個(gè)神經(jīng)元鋒電位（spike）的電活動(dòng)[2-3]。與細(xì)胞內(nèi)活動(dòng)記錄不同，這些細(xì)胞外記錄不是直接接觸神經(jīng)元的，而是需要處理電極記錄的信號(hào)來提取電極周圍不同細(xì)胞發(fā)出的鋒電位信號(hào)，鋒電位信號(hào)的檢測(cè)性能是隨后信號(hào)特征檢測(cè)及聚類等智能分析的基礎(chǔ)和前提，在大腦探索及腦類疾病的預(yù)防檢測(cè)中意義重大，應(yīng)用范圍從神經(jīng)編碼、腦類認(rèn)知到軍事醫(yī)學(xué)等[1，4]。

鋒電位檢測(cè)的準(zhǔn)確性影響著后續(xù)腦電信號(hào)特征提取和聚類精確度[5]。目前已有多名研究人員和團(tuán)隊(duì)在此領(lǐng)域做出貢獻(xiàn)。上個(gè)世紀(jì)由Turin[6]最早提出的檢測(cè)到鋒電位的方法為匹配濾波，它能從單電極采集的信號(hào)中識(shí)別出神經(jīng)元，匹配濾波的思想是將觀察到的信號(hào)匹配到神經(jīng)元鋒電位模板，在誤差范圍內(nèi)匹配的信號(hào)將被認(rèn)為是鋒電位。然而，該方法不適應(yīng)現(xiàn)代多電極采集的鋒電位信號(hào)檢測(cè)，且檢測(cè)精度較低。為了減少計(jì)算復(fù)雜性、節(jié)約時(shí)間并適應(yīng)多電極針信號(hào)采集，文獻(xiàn)[7]使用了振幅閾值法，即只有當(dāng)振幅越過預(yù)定義的閾值時(shí)檢測(cè)到的信號(hào)為鋒電位信號(hào)，但是根據(jù)數(shù)據(jù)標(biāo)準(zhǔn)差公式或人工手動(dòng)選擇的閾值雖在高信比噪數(shù)據(jù)中有較高檢測(cè)準(zhǔn)確性，但其對(duì)于低信噪比數(shù)據(jù)檢測(cè)精度仍有待提高。隨著小波在鋒電位檢測(cè)中的使用，低信噪比下的鋒電位檢測(cè)精確度明顯提升，但其在計(jì)算上要求很高，此外，小波檢測(cè)算法的精度在很大程度上取決于其對(duì)母小波的選擇。近年來，具有瞬時(shí)特性和較低計(jì)算資源需求優(yōu)勢(shì)的非線性能量檢測(cè)算法（non-linear energy，NEO）在電信號(hào)的檢測(cè)中也得到了廣泛應(yīng)用[8-9]，但當(dāng)NEO應(yīng)用于低信噪比植入式腦電數(shù)據(jù)時(shí)，檢測(cè)精度不能滿足信號(hào)分析的需求。文獻(xiàn)[10]在NEO基礎(chǔ)上提出了一種混合的神經(jīng)元峰值檢測(cè)算法——縮放性能量檢測(cè)算法（scaled energy operators，SEO），雖降低了功耗，但其鋒電位檢測(cè)精度有待進(jìn)一步提高。文獻(xiàn)[11]提出了一種基于稀疏表示和形態(tài)成分分析的鋒電位檢測(cè)算法，加強(qiáng)了該算法的抗噪能力以及適應(yīng)能力，但耗時(shí)長(zhǎng)，需要提高求解速度。在文獻(xiàn)[12]中，為提高鋒電位信號(hào)的檢測(cè)能力，提出了基于CUDA的實(shí)時(shí)鋒電位分析系統(tǒng)框架[12]，但是其文獻(xiàn)中假設(shè)的神經(jīng)元信號(hào)和噪聲都是高斯分布，對(duì)檢測(cè)含不確定背景噪聲的信號(hào)可能會(huì)有較高的錯(cuò)誤率。文獻(xiàn)[13]針多通道微電極陣列記錄的鋒電位的微弱性以及易受干擾特性，結(jié)合主成分分析（PCA）小波分析，提出PCA-小波（PCAW）與整體平均經(jīng)驗(yàn)?zāi)B(tài)分解（EEMD）聯(lián)合的去噪新方法（PCWE），但是該方法在低信噪比信號(hào)中表現(xiàn)較差，需加強(qiáng)該算法的魯棒性。由此可見，根據(jù)以上文獻(xiàn)所總結(jié)，鋒電位檢測(cè)的難點(diǎn)在于所提檢測(cè)算法既要適應(yīng)現(xiàn)代多電極針數(shù)據(jù)采集，鋒電位波形隨電極針位置或細(xì)胞內(nèi)電流的類型以及空間分布造成的低幅值特點(diǎn)，也要在因背景環(huán)境干擾造成含不確定背景噪聲數(shù)據(jù)中保持較高的檢測(cè)精度。

綜上所述，為了解決鋒電位檢測(cè)過程中因幅值較低和不確定背景噪聲而出現(xiàn)的漏檢、誤檢情況，本文提出了一種基于啟發(fā)式閾值的鋒電位檢測(cè)算法。在該方法中首先對(duì)于一些低幅值鋒電位信號(hào)，優(yōu)化橢圓濾波器參數(shù)，降低有用信號(hào)衰減程度，保留原始信號(hào)中幅值較小的鋒電位，然后使用啟發(fā)式閾值檢測(cè)公式，有效降低復(fù)雜采集環(huán)境中引入的混雜噪聲干擾，實(shí)現(xiàn)腦電信號(hào)鋒電位自動(dòng)檢測(cè)。

1 方法

本文在植入式腦電采集數(shù)據(jù)的鋒電位的提取過程中，優(yōu)化橢圓濾波器，保留了低幅值鋒電位信號(hào)，同時(shí)結(jié)合針對(duì)不確定背景噪聲的啟發(fā)式閾值檢測(cè)公式提出了一種新型鋒電位檢測(cè)算法，且對(duì)相關(guān)參數(shù)不斷調(diào)整優(yōu)化，將其應(yīng)用在不確定背景噪聲的獼猴肢體伸展抓握運(yùn)動(dòng)范式下獲得的真實(shí)數(shù)據(jù)，以及具有低幅值鋒電位信號(hào)的細(xì)胞外模擬記錄數(shù)據(jù)中的5個(gè)不同信噪比數(shù)據(jù)中進(jìn)行檢測(cè)驗(yàn)證。

1.1 數(shù)據(jù)來源

實(shí)驗(yàn)數(shù)據(jù)部分主要由兩部分組成：真實(shí)數(shù)據(jù)和模擬數(shù)據(jù)。真實(shí)數(shù)據(jù)來自于獼猴肢體伸展抓握運(yùn)動(dòng)范式，該范式中，獼猴需要伸手到指定目標(biāo)位置處對(duì)目標(biāo)物體進(jìn)行抓握，在本次實(shí)驗(yàn)中，使用前2個(gè)通道的記錄數(shù)據(jù)做相關(guān)分析。模擬數(shù)據(jù)來自英國(guó)萊斯特大學(xué)神經(jīng)工程實(shí)驗(yàn)室提供的細(xì)胞外模擬記錄數(shù)據(jù)[14-15]，該數(shù)據(jù)集源于真實(shí)的數(shù)據(jù)庫(kù)及真實(shí)的模擬細(xì)胞外記錄，團(tuán)隊(duì)根據(jù)噪聲的振幅和頻譜分布對(duì)背景噪聲的模擬源于疊加大量不同振幅與頻譜的鋒電位信號(hào)，另外再在每個(gè)數(shù)據(jù)集中添加1個(gè)低信噪比、低幅值和高方差特性的多單元鋒電位和2個(gè)頻率較低但振幅相對(duì)較大的單單元鋒電位。數(shù)據(jù)集中包含5個(gè)不同信噪比數(shù)據(jù)集，其信噪比按照給出的數(shù)據(jù)合成方式在該文獻(xiàn)中分別標(biāo)明為4，4，2，2，3[14]。

1.2 濾波

信號(hào)通常都有自身的幅頻特性，且經(jīng)由電極采集到的原始信號(hào)伴隨噪聲，考慮到鋒電位信號(hào)本身特性，對(duì)于采集到的腦電信號(hào)要進(jìn)行濾波操作。一般神經(jīng)學(xué)定義突觸間的電活動(dòng)是低頻慢波，稱為局部場(chǎng)電位，通常在300 Hz以下。過高頻率的波形都是噪聲，一般在3 000 Hz以上。研究顯示，數(shù)字濾波器的非線性相位響應(yīng)可很大程度上扭曲鋒電位形狀[16]，且會(huì)改變信號(hào)偽影的外觀，使其看起來與真實(shí)的脈沖類似，胞體的神經(jīng)沖動(dòng)都是比較高頻的快波，故采用零相位濾波，使用300～3 000 Hz的帶通橢圓濾波器對(duì)原始信號(hào)進(jìn)行濾波，橢圓濾波器在通帶和阻帶內(nèi)都是等波紋的逼近方式，相同階數(shù)下相比于巴特沃斯、切比雪夫?yàn)V波器通帶與阻帶的逼近特性良好[17]。MATLAB的信號(hào)處理工具箱提供了設(shè)計(jì)橢圓濾波器的函數(shù)：ellip函數(shù)。

利用其可以設(shè)計(jì)高通（當(dāng)ftype=high時(shí)）、或帶阻濾波器（當(dāng)ftype=stop，且Wn=[W1,W2]時(shí)）。參數(shù)n表示最小階數(shù)，參數(shù)Rp表示通帶最大波紋度（單位dB），參數(shù)Rs表示阻帶最小衰減（單位dB），對(duì)于文中所評(píng)估的模擬數(shù)據(jù)，使用文獻(xiàn)[18]所給的濾波器參數(shù)會(huì)導(dǎo)致一些在峰值較低且在閾值線的鋒電位信號(hào)不能被檢測(cè)，故將濾波器的Rp至1×10-6dB，有用信號(hào)盡可能無衰減的通過，保留原始信號(hào)中幅值較低信號(hào)的波形。

對(duì)來自萊斯特大學(xué)神經(jīng)工程實(shí)驗(yàn)室的模擬細(xì)胞外記錄數(shù)據(jù)集中的第3個(gè)數(shù)據(jù)集進(jìn)行鋒電位檢測(cè)實(shí)驗(yàn)，此數(shù)據(jù)集所包含的2個(gè)單單元鋒電位信號(hào)具備低幅值特點(diǎn)，如上圖1所示，橫坐標(biāo)表示Rp值，縱坐標(biāo)表示精確度，圖1中的曲線則為各自對(duì)應(yīng)的鋒電位檢測(cè)精度，從圖中可以看出隨著橢圓濾波器中的最大波紋度的數(shù)值的減小，鋒電位檢測(cè)精確度增大并逐漸趨于平穩(wěn)，由于考慮到濾波器的引入會(huì)對(duì)原有信號(hào)帶來?yè)p耗，所以選取相對(duì)較小的Rp值。

圖1 不同Rp檢測(cè)精度圖Fig.1 Different Rp test accuracy diagram

1.3 閾值檢測(cè)方法

本文使用四種閾值檢測(cè)方法，分別對(duì)模擬數(shù)據(jù)和真實(shí)數(shù)據(jù)進(jìn)行鋒電位閾值檢測(cè)并對(duì)比。在振幅閾值法中通常使用標(biāo)準(zhǔn)差閾值檢測(cè)算法，表述為以下表達(dá)式，根據(jù)背景噪聲對(duì)信號(hào)的影響，使用設(shè)置為4倍的噪聲標(biāo)準(zhǔn)差為閾值來進(jìn)行鋒電位的檢測(cè)[11]，表達(dá)式中的x是經(jīng)過濾波器得到的濾波信號(hào)。

能量算法主要測(cè)量信號(hào)與其導(dǎo)數(shù)之間的交叉能量，這使得它們適合于峰值檢測(cè)，因?yàn)殇h電位峰值是神經(jīng)信號(hào)的最高能量成分。Mukhopadhyay和Ray提出將非線性能量算法（NEO），如公式（3）所示用于鋒電位信號(hào)的閾值檢測(cè)。公式（4）中離散加速能量算法DEAO（discrete energy acceleration operator）在文獻(xiàn)[10]中也被用于鋒電位的檢測(cè)，且在此文獻(xiàn)中針對(duì)能量算法做了進(jìn)一步改進(jìn)，提出了公式（5）并將其用于低信噪比模擬數(shù)據(jù)，使其鋒電位檢測(cè)精度有了一定程度的提高。

式中，增加了一個(gè)參數(shù)p，經(jīng)過經(jīng)驗(yàn)取值，當(dāng)p=1.6時(shí)，鋒電位檢測(cè)的精確度最高。

2 結(jié)果與討論

2.1 模擬數(shù)據(jù)

2.1.1 多種信噪比下鋒電位檢測(cè)精度

文中采用精確度（accuracy，ACC）作為衡量鋒電位檢測(cè)準(zhǔn)確度的標(biāo)準(zhǔn)：

其中上述公式中正檢數(shù)（truly detected spikes，TDS）為鋒電位檢測(cè)算法檢出的正確的鋒電位信號(hào)，漏檢數(shù)（missed spikes，MS）為沒有被檢測(cè)出來的鋒電位信號(hào)，誤檢數(shù)（false alarms，F(xiàn)A）為檢測(cè)出來的不是鋒電位信號(hào)的誤檢個(gè)數(shù)。

表1是萊斯特大學(xué)神經(jīng)工程實(shí)驗(yàn)室提供的細(xì)胞外模擬記錄數(shù)據(jù)的5個(gè)數(shù)據(jù)集的鋒電位總數(shù)及多單元、單單元鋒電位所占百分比，在所提供細(xì)胞外模擬記錄數(shù)據(jù)集中，模擬數(shù)據(jù)3和4中擁有兩個(gè)低幅值且形狀相似的單單元鋒電位信號(hào)。表2是3種不同的鋒電位檢測(cè)算法和本文所提檢測(cè)算法分別在含不同信噪比、不同幅值大小鋒電位信號(hào)的模擬數(shù)據(jù)集上的TDS、MS、FA及ACC對(duì)比。

表1 模擬數(shù)據(jù)集中各種峰值單位的百分比Table 1 Percentages of various peak units in simulated dataset

表2 不同鋒電位檢測(cè)算法檢測(cè)精度Table 2 Detection accuracy of different spike detection algorithms

表1中的S_1到S_5是萊斯特大學(xué)提供的5個(gè)不同信噪比的細(xì)胞外模擬記錄數(shù)據(jù)，其中，表格第2列代表各個(gè)數(shù)據(jù)集中的鋒電位總個(gè)數(shù)，多單元與兩個(gè)單單元鋒電位所占百分比分別占據(jù)表的第3到5列，不難看出，每個(gè)數(shù)據(jù)集中兩個(gè)單單元鋒電位信號(hào)所占比例相似，而在數(shù)據(jù)集1與數(shù)據(jù)集5中多單元鋒電位信號(hào)所占比例較大，在90%以上，且鋒電位的總體數(shù)目相對(duì)較少。數(shù)據(jù)集2、3和4中鋒電位總數(shù)在3 000以上，多單元鋒電位信號(hào)所占比例為67%左右。

從表2可以看出，在鋒電位檢測(cè)精確度上，文獻(xiàn)[10]中所提到的SEO檢測(cè)算法要高于文獻(xiàn)[9]中的NEO和DEAO，但是在本文所提算法中，鋒電位檢測(cè)精度相較于SEO算法在每個(gè)細(xì)胞外模擬數(shù)據(jù)集上都有一定的提高，這其中對(duì)于低信噪比的數(shù)據(jù)5來說，本文提出的鋒電位檢測(cè)算法對(duì)鋒電位信號(hào)的檢測(cè)精確度相對(duì)于給出的NEO、DEAO和SEO算法精度都要高，證明該方法可用于低信噪比信號(hào)檢測(cè)。另外，在上表中還可觀察到本算法在擁有兩個(gè)較低幅值單單元鋒電位信號(hào)的數(shù)據(jù)3中檢測(cè)精確度也相對(duì)較高。為了更直觀地表現(xiàn)出不同閾值檢測(cè)算法的檢測(cè)精確度，將其繪制成如圖2所示。

圖2 4種閾值檢測(cè)算法精確度圖Fig.2 Accuracy of four threshold detection algorithms

圖2中圓圈、方框和三角形分別表示5個(gè)模擬數(shù)據(jù)集分別經(jīng)過NEO、DEAO和SEO檢測(cè)算法的鋒電位信號(hào)檢測(cè)精度，五角星代表本文所提算法在不同信噪比下細(xì)胞外模擬數(shù)據(jù)集的鋒電位信號(hào)檢測(cè)精度。從圖中的折線圖中可以看出在不同信噪比下，數(shù)據(jù)中鋒電位信號(hào)的檢測(cè)精度是不一致的，且使用本文所提鋒電位檢測(cè)算法的精確度比其余能量算法都要高。

2.1.2 不同信噪比單單元鋒電位檢測(cè)精度

表3是這兩個(gè)數(shù)據(jù)集分別經(jīng)過兩種不同鋒電位檢測(cè)算法計(jì)算并統(tǒng)計(jì)得出的漏檢數(shù)和誤檢數(shù)。

表3 數(shù)據(jù)3、4單元鋒電位信號(hào)MS與FA數(shù)統(tǒng)計(jì)表Table 3 Data 3，4 single spike signal MS and FA statistical table

在表3中，第1列和第3列為經(jīng)過文獻(xiàn)[15]中的鋒電位檢測(cè)算法得出的漏檢數(shù)和誤檢數(shù)，剩余兩列為本文提出的檢測(cè)算法得出的漏檢數(shù)和誤檢數(shù)。在表3中可以看出，經(jīng)過所提檢測(cè)算法檢測(cè)后的漏檢數(shù)為0，且誤檢數(shù)較中值檢測(cè)公式較低，因此本文方法對(duì)于模糊在閾值線附近的低幅值鋒電位信號(hào)具有較優(yōu)越的檢測(cè)性能。

2.2 真實(shí)數(shù)據(jù)

在真實(shí)的獼猴肢體伸展抓握運(yùn)動(dòng)范式腦電信號(hào)采集過程中，會(huì)不可避免地受到噪聲影響，造成信號(hào)背景噪聲的不確定性，同時(shí)，真實(shí)的鋒電位形態(tài)學(xué)特征也大有不同，幅值不盡相同，故而將提出的鋒電位檢測(cè)算法應(yīng)用在上述范式所采集的真實(shí)數(shù)據(jù)中進(jìn)行鋒電位信號(hào)檢測(cè)，并繪制出實(shí)驗(yàn)結(jié)果如圖3所示。

圖3 真實(shí)數(shù)據(jù)鋒電位閾值檢測(cè)結(jié)果圖Fig.3 Result diagram of realdata spike threshold detection

從圖3中可以看見一條橫亙于縱坐標(biāo)的直線，這便是經(jīng)過鋒電位檢測(cè)算法計(jì)算得出的閾值線，同時(shí)在圖中也可觀察到一些低幅值的在閾值線附近的鋒電位信號(hào)被檢測(cè)出。

為進(jìn)一步驗(yàn)證由所提算法檢測(cè)到的鋒電位信號(hào)，對(duì)獼猴肢體伸展抓握運(yùn)動(dòng)范式下采集得到的真實(shí)植入式腦電數(shù)據(jù)分別使用主成分分析以及K多元均值方法進(jìn)行了特征提取和聚類操作。這其中，數(shù)據(jù)具有不確定背景噪聲，同時(shí)在進(jìn)行鋒電位信號(hào)檢測(cè)前幅值大小未知。

在圖4（a）和（b）中，分別使用了獼猴肢體伸展抓握運(yùn)動(dòng)范式下采集得到的腦電信號(hào)的前兩個(gè)通道數(shù)據(jù)，經(jīng)過本文鋒電位檢測(cè)算法檢測(cè)后，再經(jīng)上述特征提取并聚類后的鋒電位信號(hào)結(jié)果圖。從圖中清晰地可以看出每個(gè)通道包含的不同的鋒電位信號(hào)，說明本文鋒電位檢測(cè)算法可用于真實(shí)采集的含不確定背景噪聲的植入式腦電數(shù)據(jù)，且從兩個(gè)圖中可以觀察到，真實(shí)數(shù)據(jù)的幅值較低，從側(cè)面驗(yàn)證了本文算法對(duì)低幅值鋒電位信號(hào)的提取能力。

圖4 真實(shí)數(shù)據(jù)鋒電位分選結(jié)果圖Fig.4 Result diagram of realdata spike sorting

3 結(jié)束語

本文主要針對(duì)含不確定背景噪聲且低幅值的腦電數(shù)據(jù)采集信號(hào)提出一種基于啟發(fā)式閾值的鋒電位檢測(cè)算法。將其用于低幅值且含不確定背景噪聲的真實(shí)數(shù)據(jù)中進(jìn)行驗(yàn)證，取得了良好的實(shí)驗(yàn)效果，同時(shí)也在多種信噪比的模擬細(xì)胞記錄數(shù)據(jù)集下進(jìn)行了相關(guān)鋒電位檢測(cè)實(shí)驗(yàn)并統(tǒng)計(jì)。實(shí)驗(yàn)結(jié)果表明，該鋒電位檢測(cè)算法對(duì)于低幅值的鋒電位信號(hào)檢測(cè)精度相對(duì)于一般或者改進(jìn)的能量算法的檢測(cè)精度要高一些，最優(yōu)閾值下NEO的平均檢測(cè)準(zhǔn)確度為48.07%，DEAO為58.19%，SEO的檢測(cè)準(zhǔn)確率平均為63.11%，而本文所提算法達(dá)到了65.21%，在多種信噪比數(shù)據(jù)中檢測(cè)精度相對(duì)突出。算法的開發(fā)也可同時(shí)減少專家手動(dòng)標(biāo)注閾值的復(fù)雜性，節(jié)省人力資源。但值得注意的是，即使新提出算法表現(xiàn)良好，但平均精確度仍待提升，所以要求研發(fā)者繼續(xù)開發(fā)高效且精確度較高的鋒電位檢測(cè)算法。

致謝 本文實(shí)驗(yàn)中所使用的獼猴伸展抓握運(yùn)動(dòng)范式數(shù)據(jù)由軍事醫(yī)學(xué)研究院軍事認(rèn)知與腦科學(xué)研究所前沿交叉學(xué)科研究室提供，特別感謝該團(tuán)隊(duì)對(duì)本項(xiàng)工作的支持。