盛 宙， 顏秉勇， 周家樂， 王慧鋒

（華東理工大學(xué)信息科學(xué)與工程學(xué)院，上海 200237）

近年來，納米孔電流檢測技術(shù)的發(fā)展取得了顯著的進(jìn)步，并且已經(jīng)成為生物分子研究領(lǐng)域的重要方法之一，被廣泛應(yīng)用于DNA[1-3]、RNA[4-5]、蛋白質(zhì)[6-7]以及其他納米顆粒[8-9]的分析。納米孔電流檢測技術(shù)基于庫爾特計(jì)數(shù)器(Coulter counter)原理：當(dāng)分析物通過兩端施加偏置電壓的納米孔時(shí)，會阻塞納米孔的部分體積，使得納米孔中的離子電流發(fā)生變化，從而形成阻斷電流；當(dāng)分析物穿過納米孔后，離子電流信號得以恢復(fù)，對這些阻斷電流事件進(jìn)行統(tǒng)計(jì)，就可以對通過納米孔的粒子進(jìn)行計(jì)數(shù)并且獲得和粒子相關(guān)的信息。然而，在實(shí)際過程中，采集到的納米孔電流信號十分微弱，通常是pA 級別，且容易被噪聲淹沒，因此，微弱的電流信號給納米孔電流信號的分析帶來了巨大的挑戰(zhàn)。

目前，對納米孔電流信號的研究主要集中在對阻斷事件的研究，包括阻斷電流幅值和阻斷事件的持續(xù)時(shí)間，這些特征能夠反映通過納米孔的分析物的一些特性，包括生物特性以及化學(xué)特性等[10]。從事件持續(xù)的時(shí)間上看主要分為短阻斷事件以及多級阻斷事件。短阻斷事件表現(xiàn)為阻斷電流幅值較小以及持續(xù)的時(shí)間較短，持續(xù)時(shí)間一般小于500 μs[11]，且由于其幅值一般較小，易被噪聲淹沒。在對短阻斷事件的處理上，需要預(yù)先對信號進(jìn)行濾波處理，但會使部分短阻斷信號發(fā)生失真。針對濾波帶來的短阻斷信號失真的問題，可以采用基于二階差分校正算法（Second-Order Differential Method，DBC）[10,12]、基于斜坡的修正方法[11]、半幅全寬（Full Width Half Maximum,FWHM）[13]以及等效電路的方法[14]對失真的短阻斷事件進(jìn)行修正，實(shí)驗(yàn)結(jié)果表明這些方法均能夠在一定程度上修正失真的短阻斷信號。多級事件具有較大的電流幅值以及較長的持續(xù)時(shí)間，使其相對于短阻斷事件更易提取，且受濾波的影響較小。此外，多級事件內(nèi)部包含了更多的具有差異化的電流幅值，所承載的信息相對于短阻斷事件更多。通過對不同類型的事件進(jìn)行分析可以獲得和分析物相關(guān)的一些信息。Arvind 等[14]建立了納米孔通道的等效電路模型，能夠有效提取短阻斷事件的信息，但是對多級事件的處理效果不佳。Raillon 等[15]將累加和算法用于固體納米孔電流多級事件信息的提取，但是該方法并不適用于短阻斷事件信息的提取。Zhang 等[16-17]對短阻斷事件以及多級事件分別構(gòu)建不同的隱馬爾科夫模型來提取阻斷事件的信息。

上述方法都是對短阻斷事件或者多級事件進(jìn)行針對性的處理，無法實(shí)現(xiàn)對同時(shí)包含短阻斷事件以及多級事件的真實(shí)納米孔電流信號進(jìn)行統(tǒng)計(jì)分析。對此，本文提出了一種基于灰度圖的方法實(shí)現(xiàn)對同時(shí)包含短阻斷事件以及多級事件的納米孔電流信號的處理，并最終獲得相應(yīng)的事件統(tǒng)計(jì)信息，包括事件類型、幅值以及事件持續(xù)時(shí)間。在數(shù)據(jù)預(yù)處理階段，為避免噪聲淹沒短阻斷事件，采用小波閾值去噪方法抑制噪聲。研究表明，小波閾值去噪能夠在去噪的同時(shí)保留信號的細(xì)節(jié)信息，且不會造成信號的嚴(yán)重失真[18-19]。在事件檢測階段，將預(yù)處理后的納米孔電流信號轉(zhuǎn)換成相應(yīng)的灰度圖，并選擇合適的閾值將灰度圖進(jìn)行二值化處理，以提取出阻斷事件。對于閾值法自身的缺陷對檢測出的短阻斷事件持續(xù)時(shí)間的影響，利用模糊C 均值（Fuzzy C-Means，F(xiàn)CM）[20-22]方法實(shí)現(xiàn)對短阻斷事件持續(xù)時(shí)間的修正，使修正后的短阻斷事件的持續(xù)時(shí)間更加接近實(shí)際的事件持續(xù)時(shí)間。對于檢測出的多級事件，利用簡單的線性迭代聚類（Simple Linear Iterative Clustering，SLIC）方法對多級事件電流信號對應(yīng)的灰度圖分割以檢測出臺階數(shù)以及定位出每一個(gè)臺階的起始與終止位置，再根據(jù)每個(gè)臺階內(nèi)部對應(yīng)的像素點(diǎn)值計(jì)算相應(yīng)臺階的平均電流，得到最終的多級事件信號。仿真結(jié)果表明，本文方法能夠有效地檢測出短事件以及多級事件，且用去噪后的信號進(jìn)行阻斷事件檢測的準(zhǔn)確率要高于直接對信號進(jìn)行阻斷事件檢測。在短阻斷事件的修正上，基于FCM 的方法比DBC 和FWHM 方法有著更好的修正效果。SLIC 超像素算法能夠?qū)Χ嗉壥录呐_階進(jìn)行有效分割，實(shí)現(xiàn)多級事件臺階的自動檢測。

1 小波閾值去噪原理

1.1 納米孔電流模型

納米孔電流的模型可以等效為[15]

式中： i (t) 為 納米孔電流； is(t) 為有效信號電流；ie(t)為噪聲電流；t 為時(shí)間序列，t= 1 ,2,3,4,··· ,N，N 為信號長度，且噪聲為高斯白噪聲，滿足式（2）給出的高斯正態(tài)分布模型：

其中：i 為t 時(shí)刻實(shí)際的噪聲值；高斯白噪聲時(shí)μ 為零；σ 為噪聲標(biāo)準(zhǔn)差，σ 由納米孔狀態(tài)以及數(shù)據(jù)采集裝置的噪聲分布決定。

1.2 小波分析基本理論

設(shè) Ψ (t)∈L2(t) ，其中Ψ(t)表示一個(gè)平方可積的函數(shù)，即能量有限的空間，若其傅里葉變換滿足允許條件，

此時(shí)，稱Ψ(t)為一個(gè)基本小波或者母小波。將母小波Ψ(t)經(jīng)伸縮和平移，得到小波序列，

式中：a 為伸縮因子，且a≠0；b 為平移因子。對于任意函數(shù) f (t)∈L2(R) ，其連續(xù)小波變換為

1.3 小波閾值去噪步驟

(1)選擇合適的小波基、分解尺度對含有噪聲的信號進(jìn)行離散小波變換，得到不同分解尺度下的小波系數(shù)；

(2)選擇合適的閾值，利用閾值函數(shù)分別對不同尺度下的小波系數(shù)進(jìn)行處理；

(3)利用處理過的小波系數(shù)進(jìn)行信號重構(gòu)，恢復(fù)信號。

2 納米孔阻斷事件檢測

納米孔電流信號的分析依賴于阻斷事件的檢測，且事件檢測的準(zhǔn)確程度會影響最終的分析結(jié)果。閾值法由于其計(jì)算簡便的優(yōu)點(diǎn)，在阻斷事件的檢測中得到了廣泛的應(yīng)用，但是該方法受噪聲影響較大。圖1 示出了閾值法檢測阻斷電流事件的原理。選擇一個(gè)合適的閾值N×σ，N 通常取為5，當(dāng)電流偏離基線Ibase的距離超過該閾值時(shí)，即認(rèn)為事件發(fā)生，把該點(diǎn)作為事件的起始位置（S），并把電流值超過閾值的第1 個(gè)點(diǎn)作為事件的結(jié)束位置（E）。起止點(diǎn)與終止點(diǎn)之間的時(shí)間差為事件的持續(xù)時(shí)間。

圖 1 閾值法檢測阻斷事件的示意圖Fig. 1 Diagram of blockade event detected by threshold method

實(shí)驗(yàn)數(shù)據(jù)表明，實(shí)際的納米孔電流信號的基線并不是一個(gè)固定的常數(shù)，而是伴隨著低頻波動。為了減小基線波動對事件檢測的影響，Gu 等[10]采用滑動窗口平均的方法，利用局部的閾值來檢測事件，從而避免了基線波動對事件檢測的影響。Raillon 等[15]采用自適應(yīng)閾值的方法減小基線波動引入的事件檢測誤差。Plesa 等[23]采用迭代檢測算法，迭代計(jì)算基線電流，但是迭代的方式會引入較大的計(jì)算量。

2.1 納米孔電流的灰度化

實(shí)際上，阻斷事件可以認(rèn)為是不同大小的阻斷電流深度表征，而阻斷深度可以認(rèn)為是事件達(dá)到穩(wěn)態(tài)時(shí)對應(yīng)的電流或者未達(dá)穩(wěn)態(tài)時(shí)的峰值電流相對于基線電流的差值，因此只要是阻斷事件，都和基線電流的幅值存在一定的差值，且差值的大小表征了阻斷電流的深度，同時(shí)可以認(rèn)為基線電流的阻斷深度為零。利用這種具有差異性的阻斷電流深度可以生成具有不同亮度的灰度圖。以圖2（a）的仿真信號為例，其灰度圖如圖2（b）所示，其中顏色的深淺代表了不同電流幅值的大小，顏色越深表明阻斷深度越小，顏色越淺表明阻斷深度越大，相近的、連續(xù)的顏色所占據(jù)的時(shí)間寬度表示相應(yīng)事件的持續(xù)時(shí)間。此外，利用電流信號的灰度圖方法可以在后續(xù)的多級事件的臺階檢測中有效地利用數(shù)據(jù)點(diǎn)的位置信息。

圖 2 納米孔電流仿真信號（a）、灰度圖（b）以及二值化后的灰度圖（c）Fig. 2 Simulated nanopore current （a ） , grayscale （b） and thresholding grayscale（c）

從圖2（b）可以看出，在灰度圖下，阻斷事件電流和基線電流可以通過顏色有效區(qū)分，且顏色越深表明當(dāng)前位置越靠近基線電流。在對阻斷事件的檢測上，只需將事件的起始位置與終止位置檢測出來，而不需要關(guān)注阻斷事件內(nèi)部所攜帶的差異化信息。在這種情況下，可以認(rèn)為同一事件內(nèi)部所對應(yīng)的顏色可以由同一種顏色表示。為了實(shí)現(xiàn)這樣的目標(biāo)，采用二值化方法將生成的灰度圖進(jìn)行二值化處理，結(jié)果如圖2（c）所示。由于二值化方法依賴于閾值的選取，因此需要對納米孔電流信號的基線電流進(jìn)行檢測，從而得到二值化的閾值。

2.2 基線電流檢測

為了將納米孔電流信號對應(yīng)的灰度圖二值化，即分離出非事件電流以及阻斷事件電流，需要先獲得基線電流，再依據(jù)基線電流選擇合適的閾值?；瑒悠骄姆椒ㄊ腔€電流檢測中常用的方法，但是該方法受滑動窗口大小的影響較大，不同的窗口大小帶來的差異性較大，且窗口大小的確定依賴于試探以及經(jīng)驗(yàn)。為了避免窗口的選擇對基線檢測的影響，本文采用小波分解的方法實(shí)現(xiàn)對基線電流的檢測。

小波分解的小波系數(shù)可以分為細(xì)節(jié)系數(shù)與近似系數(shù)，細(xì)節(jié)系數(shù)反映了原始信號中包含的高頻部分的信息，近似系數(shù)反映了原始信號中包含的低頻部分的信息。 實(shí)際的納米孔電流信號中含有低頻波動的基線電流可以認(rèn)為是原始信號中的低頻部分所包含的某一信號成分，因此對原信號進(jìn)行小波分解后，原始信號中的低頻部分的相關(guān)信息也包含在分解后的近似系數(shù)中，因此可以利用近似系數(shù)來近似地計(jì)算原始信號中的基線電流值。

圖3 示出了采用小波分解方法以及滑動平均方法對圖2 仿真信號中的基線電流進(jìn)行計(jì)算后的對比圖，橫坐標(biāo)中每1 μs 都代表一個(gè)點(diǎn)?？紤]到阻斷事件會對滑動平均法計(jì)算出的基線電流值產(chǎn)生影響，需要選擇較大的窗口來減小這種不利的影響。作為對比，滑動窗口的大小分別選為2 000 和3 000 以及與實(shí)驗(yàn)數(shù)據(jù)長度一致的窗口大小。

圖 3 滑動平均及小波分解方法獲得的基線電流對比Fig. 3 Comparison of baselines currents acquired by wavelet decomposition and moving average

可以發(fā)現(xiàn)，兩類方法都能計(jì)算出基線電流，但是采用滑動平均方法的誤差更大。從波形上來看，小波分解方法獲得的基線電流更加平滑，而滑動窗口方法得到的基線電流波動較大，無法反映出真實(shí)的基線電流。雖然采用和實(shí)驗(yàn)數(shù)據(jù)長度相一致的窗口大小也能獲得平滑的基線電流，但是無法反映實(shí)際的基線電流的波動，且始終和實(shí)際的基線電流之間存在較大的偏差，該誤差要比基于小波分解方法獲得的基線電流的誤差更大。實(shí)際的電流信號中包含著阻斷深度較小的短阻斷事件，滑動窗口方法帶來的大誤差會使得部分短阻斷事件無法被檢測出來或者存在較大誤檢測的可能性，且計(jì)算的基線電流值易受阻斷事件的影響。而小波分解方法得到的基線電流相對平滑且誤差較小，受阻斷事件的影響較小，對短阻斷事件識檢測的影響也較小。此外，為了加快納米孔電流信號實(shí)際的處理速度，可以在利用小波閾值去噪方法對信號進(jìn)行預(yù)處理的過程中，同時(shí)計(jì)算信號中的基線電流值，從而加快整個(gè)計(jì)算過程，使整個(gè)處理過程在大數(shù)據(jù)量時(shí)更加高效。

為了實(shí)現(xiàn)事件的自動檢測，在獲得二值化后的圖像后還需自動確定事件的起始位置與終止位置。由于二值化后的灰度圖是黑白相間的樣條，分別代表非事件電流以及事件電流，且只有黑白兩種像素，通過遍歷像素，比較像素值的差異，即可近似地得出事件的起始位置與終止位置，從而檢測出阻斷事件。

2.3 模糊C 均值聚類算法

盡管通過小波分解方法能夠獲得基線電流值，但是由于閾值法自身的缺陷，在用閾值法進(jìn)行阻斷事件的檢測時(shí)會帶來不可避免的誤差，主要表現(xiàn)為檢測出的阻斷事件的持續(xù)時(shí)間要比實(shí)際的阻斷事件的持續(xù)時(shí)間短。圖4 示出了采用閾值法檢測出的短阻斷事件的起始與終止位置以及實(shí)際的短阻斷事件的起始與終止位置的對比，其中S1、E1 表示利用閾值法檢測到的阻斷事件起始與終止位置，S2、E2 表示實(shí)際的阻斷事件的起始與終止位置?？梢园l(fā)現(xiàn)采用閾值法檢測出的短阻斷事件的長度要比阻斷事件的真實(shí)長度短。

圖 4 基于閾值法的阻斷事件的檢測Fig. 4 Blockade event detected based on threshold

為了減小閾值法自身的缺陷給短阻斷事件檢測帶來的誤差，通常會對檢測出的阻斷事件的位置進(jìn)行回溯處理。Gu 等[10,12]采用二階差分方法對檢測出的短阻斷事件長度進(jìn)行修正，結(jié)果表明，該方法在一定程度上能夠減小閾值法帶來的誤差，但是二階差分方法依賴于參數(shù)的選取，不同的參數(shù)會帶來不同的擬合效果，對短阻斷事件的修正效果也會存在差異[10,12,24]。由于實(shí)際的短阻斷事件的阻斷深度、持續(xù)時(shí)間以及形狀均存在差異，利用該方法時(shí)，很難找到一組通用的二階微分項(xiàng)系數(shù)對所有類型的短阻斷事件都能進(jìn)行有效的修正。

考慮到修正方法需要具備泛化能力，本文采用模糊C 均值聚類方法實(shí)現(xiàn)對短阻斷事件的修正，利用數(shù)據(jù)之間的相似度來劃分?jǐn)?shù)據(jù)，達(dá)到修正短阻斷事件的目的。其原理在于：屬于阻斷事件的數(shù)據(jù)點(diǎn)之間存在較大的相似性，而與開孔電流的點(diǎn)之間的相似性較小，利用這種相似性區(qū)別出阻斷事件電流以及開孔電流。作為無監(jiān)督機(jī)器學(xué)習(xí)的重要方法之一，模糊C 均值聚類方法利用某種相似性度量來衡量數(shù)據(jù)之間的相似性，并將這些數(shù)據(jù)劃分為若干子集，使相似的數(shù)據(jù)點(diǎn)盡可能歸于一類，而把不相似的數(shù)據(jù)點(diǎn)劃分到不同的類中。模糊C 均值由K-means聚類算法發(fā)展而來，兩者的區(qū)別在于：使用K-means對數(shù)據(jù)聚類時(shí)，每一個(gè)數(shù)據(jù)都只屬于某一確定的類簇，這種硬劃分的方式無法有效處理數(shù)據(jù)間有重疊的情況，而模糊C 均值聚類方法以[0, 1]的概率來劃分?jǐn)?shù)據(jù)，某一數(shù)據(jù)通常并不只屬于某一個(gè)類，而是可能屬于多個(gè)類，數(shù)據(jù)屬于某一類的程度用概率來衡量，最終將數(shù)據(jù)點(diǎn)歸入概率最大的一類，實(shí)現(xiàn)數(shù)據(jù)的聚類。研究表明，當(dāng)數(shù)據(jù)間存在較多的重疊時(shí)，模糊C 均值聚類方法相比于K-means 和其他聚類算法可以獲得更好的聚類效果[25]。

設(shè) X ={x1,x2,x3,···,xN}?Rs是一個(gè)待分類數(shù)據(jù)集，模糊C 均值聚類通過最小化式（7）所示的目標(biāo)函數(shù)將數(shù)據(jù)樣本X 劃分為若干類：

其中：m 為模糊參數(shù)；N 為樣本點(diǎn)數(shù)目；C 為聚類的類別數(shù)目；μi為第i 個(gè)類別的均值；wij∈[0, 1]是樣本xj屬于類別i 的隸屬度；是衡量樣本與聚類中心的相似度的距離度量。

模糊C 均值通過迭代的方式更新模糊隸屬度wij以及聚類中心μi，更新公式如式（8）和式（9）所示。

當(dāng)

(1)使用二值化后的灰度圖初步定位出短阻斷事件的起始位置以及終止位置；

(2)采用回溯的方式將起始位置向前回溯一定的點(diǎn)，回溯的長度按以下規(guī)則選?。哼x取回溯過程中第1 個(gè)與檢測出的基線電流相交的點(diǎn)的位置，計(jì)算出該位置和起始位置相距點(diǎn)的個(gè)數(shù)N，再以相交點(diǎn)作為起始點(diǎn)繼續(xù)向前回溯N 個(gè)點(diǎn)，以2N 個(gè)點(diǎn)作為最終回溯的點(diǎn)的數(shù)目；

(3)對回溯的數(shù)據(jù)點(diǎn)進(jìn)行模糊C 均值聚類，類簇?cái)?shù)目設(shè)為2，分別代表非事件電流的點(diǎn)以及短阻斷事件的點(diǎn)。

圖5 示出了對從圖2（a）所示的仿真信號中選取的持續(xù)時(shí)間為123 μs 的短阻斷事件采用模糊C 均值的方法修正之后的事件位置。可以發(fā)現(xiàn)，采用基于閾值的方法在定位事件的起始位置（S1）以及終止位置（E1）時(shí)，引入了較大的誤差。當(dāng)采用小波分解的方法得到基線電流時(shí)，采用將第1 個(gè)超出基線電流的點(diǎn)作為事件的終止位置的方法，可以使獲得的事件終止位置（E3）更加趨向于真實(shí)的阻斷事件終止位置（E2）。對于起始點(diǎn)S1 來說，其位置離真實(shí)的起始位置（S2）較遠(yuǎn)，而修正后的起始點(diǎn)（S3）縮短了S1 和S2 之間的距離，使得修正后的起始點(diǎn)位置更加地接近真實(shí)的事件的起始位置（S2）。

圖 5 檢測出的短阻斷事件的位置以及采用模糊C 均值修正后的短阻斷事件的位置Fig. 5 Locations of detected short event and locations of short event calibrated by FCM

3 納米孔多級事件

實(shí)際的納米孔電流信號中，不僅包含短阻斷事件，還包含多級事件。多級事件即包含多個(gè)臺階電流的信號，每一個(gè)臺階電流都可以看作阻斷電流信號達(dá)到穩(wěn)定狀態(tài)，每一個(gè)穩(wěn)定的狀態(tài)稱為一個(gè)臺階。典型的二級阻斷事件的示意圖如圖6 所示。

圖6（a） 示出了多級事件對應(yīng)的幅值直方圖，共含有3 個(gè)峰值，包括一個(gè)基線電流峰值，除去基線電流峰值后的峰值數(shù)目后與阻斷電流的臺階數(shù)目一致，橫坐標(biāo)表示某一電流值出現(xiàn)的頻數(shù)。圖6（b）示出的納米孔電流是二級阻斷事件。在判斷阻斷事件的臺階數(shù)時(shí)，雖然可以通過繪制幅值直方圖的方式直觀地獲得臺階數(shù)，但是無法實(shí)現(xiàn)自動判斷臺階數(shù)，尤其在數(shù)據(jù)量巨大時(shí)，無法高效地實(shí)現(xiàn)對阻斷事件臺階數(shù)的準(zhǔn)確檢測。

圖 6 二級阻斷事件示意圖Fig. 6 Diagram of two-level blockage event

通過觀察幅值直方圖可以發(fā)現(xiàn)，屬于同一個(gè)臺階的數(shù)據(jù)點(diǎn)之間距離相近，而屬于不同的臺階之間的數(shù)據(jù)點(diǎn)相距較遠(yuǎn)，對應(yīng)到相應(yīng)的灰度圖上，可以認(rèn)為距離相近的數(shù)據(jù)點(diǎn)之間的像素顏色相近，而相距較遠(yuǎn)的數(shù)據(jù)點(diǎn)之間的像素顏色差別較大。

圖7 示出了典型的二級事件及其對應(yīng)的灰度圖?？梢园l(fā)現(xiàn)，二級事件對應(yīng)的灰度圖共有3 個(gè)顏色差異較大的色塊，其中黑色色塊對應(yīng)基線電流，其余色塊可以認(rèn)為都對應(yīng)著一個(gè)穩(wěn)態(tài)電流幅值，即一個(gè)臺階。為了實(shí)現(xiàn)對臺階數(shù)目的統(tǒng)計(jì)，只需統(tǒng)計(jì)不同的色塊個(gè)數(shù)，而計(jì)算每個(gè)臺階所包含的阻斷電流幅值以及阻斷時(shí)間，則需準(zhǔn)確獲得每個(gè)臺階的起始位置、終止位置和臺階內(nèi)部的像素值。

由于每個(gè)臺階內(nèi)部的像素值相近，可以在某種程度上認(rèn)為是同一個(gè)像素值，而不同的臺階之間的像素值差異較大，可以認(rèn)為整個(gè)圖像是由幾種特定的像素值構(gòu)成。對這些像素值進(jìn)行統(tǒng)計(jì)即可獲得多級事件的臺階數(shù)，像素的類別即為事件臺階數(shù)，除去基線電流對應(yīng)的像素外，每種類別像素值占據(jù)的時(shí)間寬度即為每個(gè)臺階的持續(xù)時(shí)間。如圖7（b）所示，紅色豎線表示相鄰臺階的邊界。

圖 7 二級阻斷事件以及相應(yīng)的灰度圖Fig. 7 Two-level blockade event and grayscale

3.1 SLIC 超像素算法

超像素[26]是指具有相似紋理、顏色、亮度等特征的相鄰像素構(gòu)成的有一定視覺意義的像素塊，通過以像素分組的方式減少待處理像素的數(shù)目，加快后續(xù)圖像處理的速度，目前已經(jīng)成為圖像分割領(lǐng)域重要的手段之一。簡單線性迭代聚類（SLIC）[27]是實(shí)現(xiàn)超像素的算法之一，其思想是利用像素之間的相似性對中心點(diǎn)附近的像素進(jìn)行聚類處理，當(dāng)聚類達(dá)到穩(wěn)定狀態(tài)時(shí)，可以認(rèn)為同一類簇下的像素形成了一個(gè)較大的超像素，從而實(shí)現(xiàn)圖像分割。與其他超像素算法相比[28]，SLIC 算法具有計(jì)算效率高、速度快且能有效地保留圖像邊緣信息等優(yōu)點(diǎn)。在多級事件數(shù)據(jù)量大的處理上，高速處理是必要的，因此選擇SLIC 算法能夠保證多級事件檢測的處理效率。

SLIC 算法的流程如下：

(1)選擇聚類的類簇?cái)?shù)C，即最終分割的超像素的個(gè)數(shù)。按照設(shè)定的超像素個(gè)數(shù)，在圖像內(nèi)均勻地分配聚類中心，在初始情況下，每個(gè)超像素包含的原始像素個(gè)數(shù)S=N/C，N 為原始像素的數(shù)目。為避免初始化的中心點(diǎn)落在圖像梯度值較大的輪廓邊界上，需在中心點(diǎn)鄰域內(nèi)重新選擇中心點(diǎn)，鄰域大小一般選為3×3。計(jì)算鄰域內(nèi)所有像素點(diǎn)的梯度值，將新的中心點(diǎn)移動到梯度值最小的像素點(diǎn)位置；

(2)分別計(jì)算中心點(diǎn)2S×2S 范圍內(nèi)的像素點(diǎn)與中心點(diǎn)之間的距離D'，D'的定義如下：

由于每個(gè)像素都會被多個(gè)中心點(diǎn)搜索到，所以每個(gè)像素點(diǎn)都會與周圍的中心點(diǎn)存在多個(gè)距離，在進(jìn)行聚類時(shí)，選擇D'最小的中心點(diǎn)作為該像素的聚類中心。

(3)根據(jù)每次的聚類結(jié)果更新聚類中心，直至兩次迭代過程中聚類中心不再發(fā)生變化，結(jié)束迭代。

圖8 示出了給SLIC 算法指定不同的超像素個(gè)數(shù)時(shí)的圖像分割效果。圖8（a）為使用4 個(gè)聚類中心對多級事件的灰度圖進(jìn)行切割的結(jié)果（綠色虛線）?？梢园l(fā)現(xiàn)，當(dāng)指定正確的切割的超像素個(gè)數(shù)時(shí)，SLIC 算法能夠?qū)崿F(xiàn)對臺階的準(zhǔn)確切割。然而在自動處理實(shí)際的納米孔多級事件時(shí)，其臺階數(shù)無法事先確定，也就無法為SLIC 算法指定分割的超像素個(gè)數(shù)。

圖 8 利用SLIC 算法對不同的超像素個(gè)數(shù)的分割結(jié)果Fig. 8 Splitting results using SLIC with different superpixels

從圖8 可以發(fā)現(xiàn)，當(dāng)給SLIC 算法指定不同的超像素個(gè)數(shù)時(shí)，會得到不同的分割結(jié)果。圖8（a）指定了4 個(gè)超像素，最終灰度圖被分成了4 個(gè)超像素。圖8（b）和圖8（c）灰度圖的分割結(jié)果也和指定的超像素個(gè)數(shù)一致。這種需要指定超像素個(gè)數(shù)的方式不利于多級事件灰度圖的準(zhǔn)確分割。

從圖8（b）和圖8（c）的分割結(jié)果可以看出，雖然超參數(shù)的指定改變了最終的分割結(jié)果，但是每一個(gè)臺階之間均被有效分割，如綠色虛線所示。超出臺階數(shù)目的超像素在臺階內(nèi)部進(jìn)行了像素分割，如紅色實(shí)線所示。因此，在實(shí)際使用SLIC 時(shí)，可以指定一個(gè)較大的超像素?cái)?shù)目產(chǎn)生過分割，再對分割后的超像素進(jìn)行超像素融合，以實(shí)現(xiàn)對屬于同一個(gè)臺階的超像素進(jìn)行合并。

3.2 超像素融合

對SLIC 算法帶來的圖像過度分割的問題，需要使用超像素融合的方法對分割的、屬于同一個(gè)臺階的超像素進(jìn)行融合，以合并完整的臺階信號。典型的超像素融合方法是利用超像素的特征信息計(jì)算超像素之間的相似性，并將相似性高的超像素合并，達(dá)到超像素融合的目的。Hu 等[29]采用譜聚類的方式實(shí)現(xiàn)對SLIC 算法得到的超像素進(jìn)行高效的聚類融合，但是該方法對復(fù)雜場景下的超像素聚類效果不理想。

結(jié)合多級事件的超像素分割的特征來看，過度的分割問題主要發(fā)生在每個(gè)臺階事件內(nèi)部，如圖8（b）和圖8（c）所示，因此需要進(jìn)行超像素融合的超像素在每個(gè)臺階內(nèi)部相鄰的超像素之間，且超像素在時(shí)間序列上具有連續(xù)的特征。為了加快融合過程以適應(yīng)大數(shù)量的處理要求，本文直接利用臺階內(nèi)部點(diǎn)之間像素值相近的特性，通過計(jì)算超像素內(nèi)部的平均像素值將具有相似像素值的超像素進(jìn)行融合。具體思路為：選擇一個(gè)合適的閾值ε，在時(shí)間序列上遍歷每一個(gè)超像素，當(dāng)相鄰的兩個(gè)超像素的平均像素值的差小于該閾值時(shí)，則認(rèn)為兩個(gè)超像素屬于同一個(gè)臺階，并將這兩個(gè)超像素融合，否則認(rèn)為這兩個(gè)超像素屬于不同的臺階，將兩個(gè)超像素之間的邊界當(dāng)作兩個(gè)臺階的邊界。此外，選取合適的ε 還可以在一定程度上抑制噪聲對超像素融合的影響。

4 仿真實(shí)驗(yàn)與結(jié)果分析

圖 9 納米孔電流信號的處理流程Fig. 9 Diagram of processing nanopore signal

納米孔電流的處理過程如圖9 所示。

4.1 納米孔阻斷事件檢測

在阻斷事件的檢測上采用圖10 (a)所示的仿真信號，分別對仿真信號疊加均值為0，標(biāo)準(zhǔn)差分別為0.4、0.8、1.2、1.6、2.0 的高斯白噪聲。可以看出，隨著噪聲標(biāo)準(zhǔn)差的逐漸增大，阻斷深度越小的短事件越不容易觀測到，給事件的提取帶來了嚴(yán)重的阻礙。為了有效地檢測出每一個(gè)有效事件，尤其是短阻斷事件，有必要對信號進(jìn)行去噪處理。圖10 仿真信號中的有效事件數(shù)為12 個(gè)，為了量化事件檢測的誤差率，將該仿真信號周期性延長，使得最終的仿真信號中包含192 個(gè)事件，其中短阻斷事件112 個(gè)，多級事件80 個(gè)。針對含噪聲信號和去噪后的信號，分別對兩類信號的灰度圖進(jìn)行二值化處理，再利用修正基線的閾值法進(jìn)行阻斷事件的自動檢測，包含短阻斷事件以及多級事件。為了比較未去噪和去噪的信號在阻斷事件檢測方面的差異，定義事件數(shù)目檢測誤差率如式（13）所示：

其中：Nd為本文方法檢測出的事件數(shù)目；Nr為實(shí)際的事件數(shù)目。表1 示出了短阻斷事件以及多級事件檢測數(shù)目的誤差率。

表 1 阻斷事件的檢測數(shù)目誤差率Table 1 Error rate of blockage events number

從表1 可以發(fā)現(xiàn)，在含有噪聲的情況下，短阻斷事件無法被準(zhǔn)確地檢測出來，且平均誤差檢測率隨著噪聲水平的增大而呈現(xiàn)增大的趨勢，表明噪聲的存在會影響短事件的檢測。在噪聲標(biāo)準(zhǔn)差σ＜0.4 pA時(shí)，本文方法能夠準(zhǔn)確檢測出未去噪信號中的多級事件，但是當(dāng)0.4 pA＜σ＜2.0 pA 時(shí)，多級事件的檢測雖然出現(xiàn)了誤差，但是噪聲標(biāo)準(zhǔn)差的增大并沒有使得檢測的多級事件數(shù)目的誤差呈現(xiàn)增大的趨勢，表明當(dāng)噪聲在一定范圍內(nèi)時(shí)，對多級事件檢測數(shù)目的影響較小。對于去噪后的信號，無論是短事件還是多級事件，其數(shù)目都能有效地獲得。

4.2 短阻斷事件的修正

為了比較檢測出的事件的位置與實(shí)際的事件的位置之間的誤差，將檢測出的事件的起始與終止時(shí)刻仿真信號中實(shí)際的起始與終止位置進(jìn)行比較。為了便于觀察，圖11 只示出了前8 200 個(gè)點(diǎn)中包含12 個(gè)阻斷事件的一段信號、事件的實(shí)際位置以及本文方法檢測出的事件位置。

圖 11 基于電流灰度圖的阻斷事件檢測Fig. 11 Blockade events of detected based on current ’s grayscale

從圖11 可以發(fā)現(xiàn)，對于短事件和多級事件，本文方法都能夠檢測出來，但是在事件的持續(xù)時(shí)間上，本文方法檢測出的事件持續(xù)時(shí)間要比實(shí)際的事件持續(xù)時(shí)間短，主要原因在于判斷阻斷事件發(fā)生時(shí)，需要電流值偏離基線電流超過某一閾值，使得檢測出的事件起始位置要后于事件實(shí)際的位置。

表2 對比了本文方法檢測出的事件長度、修正的阻斷事件長度以及阻斷事件實(shí)際長度之間的關(guān)系。由于用于事件檢測的仿真信號是對前8 200 個(gè)點(diǎn)的周期性延長，且前8 200 個(gè)點(diǎn)中包含12 個(gè)事件，因此，對前8 200 個(gè)點(diǎn)的信號中的12 個(gè)事件來說，在8 200 個(gè)點(diǎn)之后的信號中都有15 個(gè)與之完全相同的事件?？梢詫⑦@16 個(gè)相同的事件看成一組，則總計(jì)12 組。對每一組事件檢測出的事件長度分別取平均值，則可以得到本文方法檢測出的每個(gè)事件的平均長度，并將其作為檢測出的事件長度。表2 中前7 組為短阻斷事件，第8 組～第12 組為多級事件（組號與圖11 所示阻斷事件的順序并對應(yīng)），m 為1.8。

表 2 原始事件長度與本文方法事件長度的對比Table 2 Comparison between real events’ length and events’length acquired by proposed method

從表2 可以看出，對短事件來說，其檢測出的長度誤差較大，且大于長事件檢測長度的誤差。當(dāng)采用模糊C 均值方法對短事件檢測長度修正后，其長度誤差相較于未修正的短阻斷事件檢測長度的誤差得到了明顯的改善，表明使用模糊C 均值方法能夠降低短阻斷事件檢測長度的誤差。而對多級事件來說，其檢測出的長度的誤差較小，平均在10%左右。

圖12 示出了前8 200 個(gè)點(diǎn)中包含的6 個(gè)短阻斷事件（無先后順序，且與圖11 所示短阻斷事件的順序無關(guān)），采用模糊C 均值方法對檢測出的短阻斷事件進(jìn)行修正的結(jié)果。其中，藍(lán)色圓點(diǎn)表示事件實(shí)際的位置，紅色菱形表示閾值法檢測出的短阻斷事件的位置，黑色三角形表示采用模糊C 均值修正之后短阻斷事件的位置。模糊參數(shù)m 的取值為1.8。

從圖12 中可以發(fā)現(xiàn)，對檢測出的短阻斷事件用模糊C 均值進(jìn)行修正后，事件的位置與實(shí)際的阻斷事件的位置更加接近。此外，計(jì)算出的基線值使得短阻斷事件的結(jié)束位置的判斷也更加接近實(shí)際的阻斷事件結(jié)束的位置。

為了對比不同的模糊參數(shù)m 對短阻斷事件長度修正的影響，在m 最佳選取范圍[1.5,2.5]內(nèi)分別選擇1.6、1.8、2.0、2.2、2.4 作為對比，其修正結(jié)果的對比如表3 所示。從表3 可知，在m 的最佳選取范圍內(nèi)，m 值的變化對最終的短阻斷事件長度修正的影響較小。此外，由于各個(gè)短阻斷事件的幅值以及持續(xù)時(shí)間之間存在較大的差異，且其檢測易受多種因素的干擾，如基線電流是否存在較大誤差以及噪聲是否濾除等，因此無法找到一個(gè)通用的m 值對所有類型的短阻斷事件的修正都能夠達(dá)到最佳的效果，而從表3 可以看出，本文選擇m=1.8 時(shí)能夠獲得比較好的效果。

圖13 示出了使用模糊C 均值方法（FCM）、二階差分修正方法（DBC）以及半幅全寬的方法（FWHM）對短阻斷事件進(jìn)行修正后事件長度的相對誤差比較結(jié)果。相對誤差的定義如式（14）所示：

其中：Li為第i 組事件的修正后的平均長度；Lri為第i 組事件的實(shí)際長度；Lerr表示事件長度的相對誤差，其正負(fù)表示修正后的事件是否存在過修正的問題。

圖 12 短阻斷事件的長度修正的結(jié)果Fig. 12 Results of short events calibrated by FCM

圖13 中黑色水平虛線表示事件長度誤差為0 的位置?？梢园l(fā)現(xiàn)，采用FWHM 方法在短阻斷事件長度的修正上會引入較大的誤差，而DBC 和FCM 方法均可以實(shí)現(xiàn)對短阻斷事件的修正，且誤差相對較小。DBC 方法雖然能獲得比較好的修正效果，但同時(shí)也引入了過修正的問題，使得部分噪聲信號被當(dāng)作有效信號，影響短阻斷事件的分析。FCM 方法能夠在獲得較好的修正效果的同時(shí)，還有效避免了將噪聲信號當(dāng)作有效信號的問題。

4.3 多級事件的臺階檢測

多級事件的檢測主要是獲得事件電流臺階數(shù)以及獲得每個(gè)臺階之間的間隔點(diǎn)，即需要獲得每個(gè)臺階的起始與終止位置。圖14～圖17 分別示出了對4 類典型的多級事件分割后的最終結(jié)果，即先采用SLIC 方法對相應(yīng)的灰度圖進(jìn)行過度分割后，再采用超像素融合的方法將相鄰的、相似的超像素進(jìn)行融合，SLIC 方法的聚類中心C 的取值均為10。

圖14～圖17 的分割結(jié)果表明，利用SLIC 算法能夠?qū)崿F(xiàn)多級事件中臺階的有效分割，并且利用超像素融合的方式能夠解決SLIC 算法過分割的問題。實(shí)驗(yàn)結(jié)果表明，在多級事件的檢測上，SLIC 算法初始化聚類中心數(shù)目的選取并不是越大越好，過大的聚類中心數(shù)目會產(chǎn)生復(fù)雜且不規(guī)則的分割圖像，不利于分割后的超像素融合。

4.4 SLIC 算法效率

由于實(shí)際的納米孔數(shù)據(jù)包含有大量的多級事件，因此使用SLIC 算法對多級事件臺階進(jìn)行自動檢測時(shí)，需要考慮算法的計(jì)算效率問題。第一，當(dāng)多級事件包含的數(shù)據(jù)點(diǎn)逐漸增大時(shí)，采用SLIC 算法時(shí)的計(jì)算效率，即數(shù)據(jù)點(diǎn)的長度和計(jì)算時(shí)間的關(guān)系；第二，多級事件的臺階數(shù)是否對SLIC 算法效率有影響。圖18 示出了不同數(shù)據(jù)長度、不同的臺階數(shù)目下SLIC 算法的執(zhí)行時(shí)間，橫坐標(biāo)中每1 ms 包含1 000 個(gè)數(shù)據(jù)點(diǎn)。所測試數(shù)據(jù)均為已知臺階數(shù)的多級事件信號，以保證最終分割結(jié)果的準(zhǔn)確性，最終的時(shí)間包含了超像素融合所需的時(shí)間。測試環(huán)境為：Intel(R)Core(TM) i5-6300HQ CPU @2.30 GHz, 8 GB RAM, 64 bit Win10 企業(yè)版操作系統(tǒng)，軟件為Matlab2017b。

表 3 不同模糊參數(shù)m 下的阻斷事件修正的結(jié)果Table 3 Results of short events calibrated by FCM with different m

圖 13 FCM、FWHM、DBC 的修正效果對比Fig. 13 Comparison of calibrating effect on short events among FCM, FWHM and DBC

圖 14 典型單級阻斷事件的灰度圖分割結(jié)果Fig. 14 Splitting result on grayscale of one-level blockade event

圖 15 典型2 級阻斷事件的灰度圖分割結(jié)果Fig. 15 Splitting result on grayscale of two-level blockade event

圖 16 典型3 級阻斷事件的灰度圖分割結(jié)果Fig. 16 Splitting result on grayscale of three-level blockade event

圖 17 典型4 級阻斷事件的灰度圖分割結(jié)果Fig. 17 Splitting result on grayscale of four-level blockade event

圖 18 SLIC 算法的效率Fig. 18 Efficiency of SLIC

從圖18 可以發(fā)現(xiàn)，隨著數(shù)據(jù)量的增大，SLIC 算法所需的計(jì)算時(shí)間與數(shù)據(jù)點(diǎn)的長度近似呈現(xiàn)一種線性關(guān)系。當(dāng)多級事件中包含大量的數(shù)據(jù)點(diǎn)時(shí)，會帶來計(jì)算時(shí)間上的巨大開銷。針對大數(shù)據(jù)量的多級事件來說，可以利用采樣方式減少數(shù)據(jù)點(diǎn)的數(shù)目，如分別對奇數(shù)序列和偶數(shù)序列進(jìn)行分割，再將分割結(jié)果融合。此外，可以發(fā)現(xiàn)多級事件的臺階數(shù)與SLIC 算法的效率并沒有直接的關(guān)系。

5 結(jié)束語

本文提出了基于納米孔電流信號的灰度圖實(shí)現(xiàn)對阻斷事件的檢測，為納米孔阻斷事件的檢測提供了新的方法。針對閾值法對短阻斷事件的影響，采用FCM 方法對短阻斷事件進(jìn)行長度修正，通過利用數(shù)據(jù)點(diǎn)之間的相似性來進(jìn)行數(shù)據(jù)點(diǎn)的分類，即屬于開孔電流或者阻斷信號，使其效果要優(yōu)于常用的DBC 以及FWHM 兩類方法。在多級事件的檢測問題上，利用超像素算法SLIC 實(shí)現(xiàn)對其對應(yīng)灰度圖色塊的分割，能夠有效地利用電流幅值以及對應(yīng)的位置信息，使得具有相近電流幅值而位置上不連續(xù)的點(diǎn)不會被歸于同一個(gè)超像素。超像素融合有效地解決了SLIC 算法帶來的灰度圖被過分割問題，能夠?qū)儆谕慌_階且位置上相鄰的超像素進(jìn)行融合，防止屬于同一臺階的超像素被分割成多個(gè)，影響多級事件的臺階檢測的準(zhǔn)確性。短阻斷事件的準(zhǔn)確提取以及SLIC 算法對多級事件臺階的準(zhǔn)確切割，表明在納米孔數(shù)據(jù)的分析上采用圖像處理的方法的可行性，同時(shí)也能夠?yàn)楹罄m(xù)的納米孔電流阻斷事件的其他分析提供可靠的數(shù)據(jù)。