基于模糊C 均值和SLIC 的納米孔阻斷事件的識別與研究

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

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

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

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

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

1 小波閾值去噪原理

1.1 納米孔電流模型

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

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

其中：i 為t 時刻實際的噪聲值；高斯白噪聲時μ 為零；σ 為噪聲標準差，σ 由納米孔狀態以及數據采集裝置的噪聲分布決定。

1.2 小波分析基本理論

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

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

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

1.3 小波閾值去噪步驟

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

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

(3)利用處理過的小波系數進行信號重構，恢復信號。

2 納米孔阻斷事件檢測

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

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

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

2.1 納米孔電流的灰度化

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

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

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

2.2 基線電流檢測

為了將納米孔電流信號對應的灰度圖二值化，即分離出非事件電流以及阻斷事件電流，需要先獲得基線電流，再依據基線電流選擇合適的閾值。滑動平均的方法是基線電流檢測中常用的方法，但是該方法受滑動窗口大小的影響較大，不同的窗口大小帶來的差異性較大，且窗口大小的確定依賴于試探以及經驗。為了避免窗口的選擇對基線檢測的影響，本文采用小波分解的方法實現對基線電流的檢測。

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

圖3 示出了采用小波分解方法以及滑動平均方法對圖2 仿真信號中的基線電流進行計算后的對比圖，橫坐標中每1 μs 都代表一個點。考慮到阻斷事件會對滑動平均法計算出的基線電流值產生影響，需要選擇較大的窗口來減小這種不利的影響。作為對比，滑動窗口的大小分別選為2 000 和3 000 以及與實驗數據長度一致的窗口大小。

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

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

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

2.3 模糊C 均值聚類算法

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

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

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

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

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

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

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

當

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

(2)采用回溯的方式將起始位置向前回溯一定的點，回溯的長度按以下規則選取：選取回溯過程中第1 個與檢測出的基線電流相交的點的位置，計算出該位置和起始位置相距點的個數N，再以相交點作為起始點繼續向前回溯N 個點，以2N 個點作為最終回溯的點的數目；

(3)對回溯的數據點進行模糊C 均值聚類，類簇數目設為2，分別代表非事件電流的點以及短阻斷事件的點。

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

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

3 納米孔多級事件

實際的納米孔電流信號中，不僅包含短阻斷事件，還包含多級事件。多級事件即包含多個臺階電流的信號，每一個臺階電流都可以看作阻斷電流信號達到穩定狀態，每一個穩定的狀態稱為一個臺階。典型的二級阻斷事件的示意圖如圖6 所示。

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

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

通過觀察幅值直方圖可以發現，屬于同一個臺階的數據點之間距離相近，而屬于不同的臺階之間的數據點相距較遠，對應到相應的灰度圖上，可以認為距離相近的數據點之間的像素顏色相近，而相距較遠的數據點之間的像素顏色差別較大。

圖7 示出了典型的二級事件及其對應的灰度圖。可以發現，二級事件對應的灰度圖共有3 個顏色差異較大的色塊，其中黑色色塊對應基線電流，其余色塊可以認為都對應著一個穩態電流幅值，即一個臺階。為了實現對臺階數目的統計，只需統計不同的色塊個數，而計算每個臺階所包含的阻斷電流幅值以及阻斷時間，則需準確獲得每個臺階的起始位置、終止位置和臺階內部的像素值。

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

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

3.1 SLIC 超像素算法

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

SLIC 算法的流程如下：

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

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

由于每個像素都會被多個中心點搜索到，所以每個像素點都會與周圍的中心點存在多個距離，在進行聚類時，選擇D'最小的中心點作為該像素的聚類中心。

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

圖8 示出了給SLIC 算法指定不同的超像素個數時的圖像分割效果。圖8（a）為使用4 個聚類中心對多級事件的灰度圖進行切割的結果（綠色虛線）。可以發現，當指定正確的切割的超像素個數時，SLIC 算法能夠實現對臺階的準確切割。然而在自動處理實際的納米孔多級事件時，其臺階數無法事先確定，也就無法為SLIC 算法指定分割的超像素個數。

圖 8 利用SLIC 算法對不同的超像素個數的分割結果Fig. 8 Splitting results using SLIC with different superpixels

從圖8 可以發現，當給SLIC 算法指定不同的超像素個數時，會得到不同的分割結果。圖8（a）指定了4 個超像素，最終灰度圖被分成了4 個超像素。圖8（b）和圖8（c）灰度圖的分割結果也和指定的超像素個數一致。這種需要指定超像素個數的方式不利于多級事件灰度圖的準確分割。

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

3.2 超像素融合

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

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

4 仿真實驗與結果分析

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

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

4.1 納米孔阻斷事件檢測

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

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

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

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

4.2 短阻斷事件的修正

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

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

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

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

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

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

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

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

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

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

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

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

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

4.3 多級事件的臺階檢測

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

圖14～圖17 的分割結果表明，利用SLIC 算法能夠實現多級事件中臺階的有效分割，并且利用超像素融合的方式能夠解決SLIC 算法過分割的問題。實驗結果表明，在多級事件的檢測上，SLIC 算法初始化聚類中心數目的選取并不是越大越好，過大的聚類中心數目會產生復雜且不規則的分割圖像，不利于分割后的超像素融合。

4.4 SLIC 算法效率

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

表 3 不同模糊參數m 下的阻斷事件修正的結果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 典型單級阻斷事件的灰度圖分割結果Fig. 14 Splitting result on grayscale of one-level blockade event

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

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

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

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

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

5 結束語

本文提出了基于納米孔電流信號的灰度圖實現對阻斷事件的檢測，為納米孔阻斷事件的檢測提供了新的方法。針對閾值法對短阻斷事件的影響，采用FCM 方法對短阻斷事件進行長度修正，通過利用數據點之間的相似性來進行數據點的分類，即屬于開孔電流或者阻斷信號，使其效果要優于常用的DBC 以及FWHM 兩類方法。在多級事件的檢測問題上，利用超像素算法SLIC 實現對其對應灰度圖色塊的分割，能夠有效地利用電流幅值以及對應的位置信息，使得具有相近電流幅值而位置上不連續的點不會被歸于同一個超像素。超像素融合有效地解決了SLIC 算法帶來的灰度圖被過分割問題，能夠將屬于同一臺階且位置上相鄰的超像素進行融合，防止屬于同一臺階的超像素被分割成多個，影響多級事件的臺階檢測的準確性。短阻斷事件的準確提取以及SLIC 算法對多級事件臺階的準確切割，表明在納米孔數據的分析上采用圖像處理的方法的可行性，同時也能夠為后續的納米孔電流阻斷事件的其他分析提供可靠的數據。