基于模糊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 算法对多级事件台阶的准确切割，表明在纳米孔数据的分析上采用图像处理的方法的可行性，同时也能够为后续的纳米孔电流阻断事件的其他分析提供可靠的数据。