基于裂缝平移的岩心CT序列图像裂缝批量提取

许康，滕奇志

（四川大学电子信息学院图像信息研究所，成都 610065）

0 引言

数字岩心由于其独特的优势，近些年来在石油和天然气行业已经得到了广泛地运用［1］。建立一个准确客观的三维数字岩心模型是应用数字岩心技术进行各种数值模拟实验的基础［2］，本文所使用到的一种建模的方法是，先通过CT逐层扫描岩心得到序列图，然后再使用相关软件批处理序列图得到处理后的二值图像，最后再逐层叠加图层即得到三维的岩心结构。

岩心CT序列图像的二维批量处理过程是构建岩心三维模型过程中的关键步骤，其目的就在于从复杂的岩心结构中提取出感兴趣的目标，如岩心的孔隙颗粒和裂缝等，以便于三维模型的构建以及三维参数的计算。岩心当中的微观裂缝对于研究岩心的连通性和渗流特性具有重要意义［3］［5］。研究某岩心裂缝的相关参数的一般步骤是，先对二维CT序列图像进行岩心裂缝目标的提取，再通过一系列算法实现数字岩心裂缝三维结构的重建，进而计算裂缝的各类参数。然而目前CT成像的对比度往往不够理想，且易受到噪声影响，导致二维批处理提取裂缝困难，且一些细小的裂缝在提取过程中容易断裂，较难做到完整的提取。

常见的数字岩心裂缝分割提取的办法一般是先进行图像增强，再通过阈值分割或区域生长等算法来提取裂缝。但是这类方法往往难以解决CT图像对比度低和噪声难以去除等问题，同时也会存在裂缝提取不完整、易断裂等情况。如果要获得好的分割效果，通常需要进行大量的人工修补工作。而在实际应用中，岩心CT序列图的数量通常较大，在批处理过程中如果每一张都需要人为修补，将耗费大量的人力和时间。文献［3］将CT序列看做一个连续变化的视频帧序列，利用区域聚类和光流阈值分割进行帧间像素运动的预测来分割裂缝，然而其时间复杂度较高，所以在工程运用中存在局限性。文献［4］采用相关性阈值分割，先对首帧序列图计算阈值，再对之后的每张序列图自动调整分割阈值，以获得较好的分割效果，然而其主要应用于孔隙提取，对裂缝目标仍然难以做到完整提取。

本文的方法同样考虑到CT序列帧间较强的相关性［7］，且经过进一步观察，裂缝在序列图逐帧递进的过程中变化较缓慢，裂缝整体呈方向性移动趋势，因此提出一种基于裂缝平移的岩心CT序列图像裂缝提取的方法。以第一帧阈值分割预提取的裂缝为模板，将其人为手工补全并去噪，然后分别对图像中的每条裂缝进行逐帧平移提取，同时根据序列图中该裂缝的真实形态变化情况来微调平移后的裂缝形态。逐帧平移变化至最后一帧，从而实现序列图裂缝的批量提取。

1 算法原理

岩心CT序列图像二维批量处理的过程中存在两个重要的图层，其中一个是CT设备扫描岩心之后获得的原始图像层，另一个就是用于保存提取之后的孔隙、颗粒和裂缝等目标用于三维重建的二值图像层，我们大部分的阈值分割和目标提取操作都在该层完成。

本文的算法是将序列图中的裂缝目标以裂缝平移的方式逐条提取，具体流程如图1所示。

图1 裂缝平移提取流程

算法的具体原理如下：针对这种存在方向性移动趋势的裂缝，通过计算出序列图像每帧之间裂缝移动的方向以及距离，就可以将上一帧二值图像层中提取到的裂缝目标移动到下一帧图层的对应位置上，以实现裂缝目标的粗提取。且由于裂缝帧间差异较小，因此可以模拟图像层裂缝的真实形态，使粗提取裂缝在保证不断裂的前提下发生细微的变化。为使变化后的裂缝与裂缝真实情况尽可能相近，因此采取一种部分区域二值化的方法对裂缝周边区域的像素点进行进一步的提取划分，同时针对区域二值化过程中可能发生的细微裂缝断裂情况，需要再结合Dijkstra算法以及区域生长算法进行裂缝的部分区域补全。再对CT图像中的每条裂缝都执行本文的算法，即可实现CT序列图像裂缝的批量提取。

1.1 预提取裂缝

先取序列图的第一幅图像进行裂缝目标的初步提取，采用常见的阈值分割方法提取即可，通过阈值算法计算出合适的阈值作为二值化图片的界限，将裂缝目标在二值图层中置为目标色，其余为背景色。具体情况如图2 所示，深色部分即为提取到的裂缝目标，它叠加在原始图像层的上方。

图2 阈值分割后裂缝断裂情况

由于裂缝相对背景对比度不高且存在噪点的干扰导致裂缝的提取效果不够理想，具体表现为：①噪点和裂缝杂糅在一起难以去除。②一些裂缝在提取的过程中由于颜色较浅而出现了断裂的情况，这会导致裂缝参数的计算出现较大误差。

因此在阈值分割之后必须对裂缝目标进行人为的修补和去噪。人工交互处理的过程主要是人为手动将断裂的裂缝连接好，以及将噪点、岩石、孔隙等非裂缝目标置为背景色。

1.2 裂缝平移

由于本文的算法采取逐条提取裂缝的模式，因此可以只针对感兴趣的裂缝目标做分割提取，进而排除一些细微的有干扰的裂缝项。

下面是裂缝单条平移提取的具体过程。人工选择覆盖其中一条裂缝某段的部分矩形区域如图3 所示，注意选择的区域块中只能包含一条裂缝否则将会出现干扰影响该裂缝的平移提取。选择的部位尽量为裂缝目标中较直的部分以更好的拟合直线计算方向及距离。

图3 裂缝平移方向及距离计算

图3 中的粗线裂缝代表上一帧处理好的裂缝目标，细线裂缝表示裂缝需要平移到下一帧的目标位置，这是由第二帧裂缝的真实位置决定的。由于CT序列图帧间较强的相关性，同时裂缝存在一定的厚度，所以在实践中前后两帧裂缝的位置偏移通常较小且存在重合的部分。

（1）在划定的矩形区域内将前一帧中粗线表示的裂缝目标像素点提取成一个集合，便于后续的计算操作。先扫描矩形框内的点寻找裂缝目标，以第一个目标色点做为起始点，执行八邻域的聚类操作。具体过程如下。

定义一个队列Q，设起始点为p(x,y)，（x,y为其对应的坐标）。并将其放入队列中，设待提取目标点的集合为P。

1）将p存入队列Q中。

2）取队列Q中的首元素，设为x，将该点存入点集P。弹出首元素，同时把x的八邻域点中为目标色且未访问过的点存入Q中。

3）当Q不为空时循环执行第二步，直到没有点满足条件此时聚类完成。

图4 裂缝目标的八领域聚类操作

执行聚类操作后裂缝将以坐标数据的形式存储在集合P中，通过索引坐标值可以获得该裂缝点位置在岩心原始图像层的具体灰度值与其在二值图层的状态，即目标或背景。

计算该集合像素点的灰度平均值E（P）。其中灰度平均值的计算公式如下：

（2）在当前矩形选框内，以E（P）为阈值进一步分割后一帧图像，达到粗提取浅色部分裂缝的目标，设后一帧图像粗提取后的集合为B。

对集合P中的像素点做线性回归拟合成一条直线，该直线反映了集合P总的数据趋势［6］。直线方程的最终表达式为hθ(x) =θ0+θ1x, 为了让拟合成的直线与实际值误差最小，求得其损失函数为：

解得：

（3）计算集合B到直线hθ(x) =θ0+θ1x在x、y方向的平均移动距离，即可将前一帧的裂缝平行移动相应的距离到下一帧。

每次平移之后的裂缝在位置上基本与前一帧裂缝的位置相一致，但是考虑到CT序列图像中的真实裂缝在每帧之间的形态肯定会存在细微的差异，因此为了模拟这种形态上的差异，考虑采用一种区域二值化的算法进行局部阈值分割。

1.3 改变裂缝平移形态

由于裂缝存在一定的宽度，所以在线性拟合时生成的直线更趋近于裂缝区域的中轴线。因此在平移裂缝的过程中也会将裂缝移动到下一帧裂缝目标的中轴线附近。考虑在中轴线吻合的情形下，裂缝的形态变化就基本只集中在其外层的部分区域，即二维点集的周边区域。因此可以保持裂缝目标中心区域不变化，只改变周边区域就可以模拟出裂缝逐帧形态的变化过程。

一条裂缝在岩心中的形态变化具体包括向外围区域扩张或者向内部区域收缩（裂缝中途断裂的情况暂不考虑），因此将裂缝目标区域内靠外一层的点集与裂缝区域最外围一圈点集作为其大概率发生形态变化的区域。

假定裂缝的主要形态变化只集中在这些区域，如图5所示。

图5 裂缝周围部分二值化区域

最靠近中心的一圈白色区域为裂缝的中心部分，中间层的灰色部分为裂缝外圈的点集，最外圈一层为背景区域。因此我们保留中心部分的状态，同时合并中间层与最外层，将其作为局部阈值分割的目标区域，对移动后的裂缝进行调整优化。

模拟的灰度图如图5，其中灰色与黑色代表裂缝区域。黑色为其中心区域，灰色为裂缝最外圈的部分，白色为背景中紧贴裂缝的一圈区域。

获得目标区域的具体做法如下：遍历裂缝目标，若其八邻域内的点存在不属于裂缝目标的点，则其为灰色区域的部分，同时其八邻域中不属于裂缝目标的点即属于白色区域。

采用最大类间方差法（OSTU 法）计算目标区域的分割阈值［7］该算法能得到类间分离的最佳分割阈值，具体算法的流程如下：

设图像的最大灰度级为L-1，统计区域内不同灰度级下像素点的个数n，区域内像素点总数为MN=n0+n1+ …+nL-1，归一化直方图分量以及不同灰度等级k下的累积和分别为：

求出使得σ2B(k)取得最大值的阈值k*，此时的k*即为最大类间方差分割阈值。通过该分割阈值对目标区域进行局部二值化分割，可以实现裂缝平移目标的形态变化，具体表现为：将平移而来但是实际不属于该帧裂缝的像素点置为背景，同时将该帧图像中裂缝扩展的像素点置为目标。

1.4 基于Dijkstra算法的裂缝自动补全

当裂缝宽度小于3 个像素时，在划分裂缝区域时就会缺少上文所述的中心区域。在这种情形下，如果该部分目标区域的原始灰度值又较高时，此时经过局部阈值分割之后，裂缝就有可能会出现小部分断裂的情况。

针对这些小部分断裂处，需要采取一种基于Dijkstra最短路径算法［10-12］的裂缝自动补全方法。

设裂缝的两处断裂处的端点分别为start 和end，以这两点为矩形对角线构建一个矩形框，并对该矩形框构造带权图G，表示如下：

其中，V代表矩形框中的像素点，即带权图中的顶点；E代表图中的一个边，图中的每个顶点（除边缘处）都有8 条边与之相交；W代表带权图中每条边的权值，具体对应原始灰度图像中顶点所代表的像素点的平均灰度值。

Dijkstra最短路径算法的具体流程如下：

（1）设置断裂裂缝的两头端点为start 和end，以两点为对角线构造矩形带权图G。

（2）创建两个集合P，Q，其中P为已标识点的最短路径集合，初始为start 点，对应的最短路径d为0，Q为未标识点的最短路径集合，包含除start 点外的所有顶点，若与start 点相邻则初始值为权值边w，否则为无穷大。

（3）从Q集合中选取d最小的点，将该点设为x并划为已标识点，确认点x的所有处于矩形框内的邻域点yi，使用已知的x至yi的权重边w[x][yi]，重新计算yi到起始点的最短距离并更新P、Q集合，计算规则是：

d[yi]= min(d[yi],d[x]+w[x][yi])

（4）重复迭代至Q集合为空时结束，此时P集合中end与起始点start的距离即最短路径。

（5）经过回溯找到start 点和end 点的最短路径路线，算法结束。

此时将位于最短路径路线上的点在图形层中置为裂缝目标，这样裂缝的小部分断裂处就实现了自动连接工作。

但是Dijkstra算法只能将断裂部分以单像素宽的路径连接成裂缝骨架，连接处与真实裂缝目标仍存在一定的差异，因此需要再通过一种区域生长算法将单像素宽的裂缝骨架横向扩展，使其与真实的裂缝厚度相一致。

具体过程如下：统计单像素宽的裂缝连接处的所有像素点的平均灰度值G，再将裂缝骨架上的所有目标点当做区域生长的初始种子点，并加入区域生长点集合A中。扫描集合A中的每一个点，若该点的八邻域范围内像素点的灰度值g（x，y）满足如下的关系表达式：

其中T为人为设定的阈值，则将该点在二值图层中置为裂缝目标，并将该点加入到A中，继续下一个点的判定。重复循环至A集合为空时，即实现了单像素裂缝的横向区域生长，此时的裂缝的小部分断裂处就实现了横向与纵向的总体自动补全。

2 实验结果分析

本文对存在裂缝平移趋势的多组CT序列图进行了实验，展示其中的两组序列图二值分割的结果如图6、图7所示。

图6 第一组CT序列图分割结果对比（续）

图6 第一组CT序列图分割结果对比

图7 第二组CT序列图分割结果对比

将本文算法与固定阈值分割算法与文献［4］所提及的相关性阈值分割算法进行对比。其中的交互式阈值分割算法即通过人为交互的方式手动设定二值分割图片的阈值，以实现肉眼可见的最佳分割效果，之后对CT序列图的每一帧都采取第一帧的固定阈值来进行分割提取。相关性阈值分割则采取首帧固定阈值，之后再对相邻帧的实际图像改变分割的阈值，使分割的效果更佳。

从对比图可以明显看出，在CT序列图中不论是采用固定阈值分割法还是相关性分割方法，其提取到的裂缝都会存在较多断裂处，且裂缝杂糅、形态复杂。

对于第一组序列图来说，其裂缝本身就较为完整，但是仍然存在部分细微裂缝的情况，本文算法相比固定阈值分割与相关性二值分割来说有较好的提取效果。对于该序列图中长裂缝目标的局部断裂处，本文算法成功地实现了完整提取，并延续到之后的每一帧提取上。

对于第二组裂缝存在较多微小裂缝干扰的情况，本文算法先通过人工交互的方式修补连接第一帧断裂裂缝，去除噪点并删去干扰的微小裂缝项，简化裂缝结构，随后平移提取目标裂缝项，成功实现了裂缝的针对性提取。因此采取本文的方法既可以选择性地排除微小裂缝的干扰，同时还能够保持了裂缝提取的完整性。

第二组提取结果的裂缝断裂处放大图如图8所示，可以看到在序列图100 帧时作为目标的裂缝仍然提取成功且没有发生断裂，同时排除了一些较浅或不明显的其他裂缝的干扰。

图8 序列图100帧断裂裂缝放大处对比

为了验证CT 序列图裂缝整体提取的正确性，将采用本文算法提取出的裂缝目标生成三维模型［9］［13］，如图9所示。

图9 本文方法三维重建模型结果

可以进一步总结出使用本文方法存在的两大优势：

（1）可以针对性地对感兴趣的裂缝目标进行提取，进一步排除细小裂缝和干扰裂缝的影响，使重建后的三维模型结构较为清晰；

（2）维持了二维裂缝目标的连续性，使得在批处理提取裂缝不产生局部断裂，最终重建的三维裂缝模型就会更加地密集且完整。

3 结语

本文依据岩心CT序列图像中裂缝目标在帧间存在较强的相关性这一特点，提出了一种基于裂缝平移的裂缝批量提取算法。在序列图首帧进行裂缝的预提取，再通过人工交互的方式补全裂缝，之后再逐条平移提取裂缝，同时还根据裂缝在图像层的真实灰度值来改变平移裂缝的形态以适应真实的裂缝变化，并实现了基于Dijkstra算法的小部分断裂裂缝的自动补全，做到了针对性地批处理提取裂缝。

通过对比二维提取结果和三维重建结果，证明了本文算法批量提取裂缝的优越性，在有效去除干扰裂缝及各类噪点的同时还能够提高裂缝批量提取的完整度，在大批量岩心CT序列图像的裂缝提取过程中有较好的应用。

本文的方法仅适用于部分裂缝帧间相关性较强的岩心CT图像。如果单条裂缝在序列图递进的过程中分裂成了两条，或者层数较多裂缝发生改变较大时，采用本文方法就难以做到一次性全部正确分割，需要考虑以第一张裂缝分裂的序列图为节点前后分别平移提取分割裂缝，这增多了人工交互环节。另外，本文算法默认裂缝最外围一圈为其大概率发生形态变化的区域，然而当裂缝帧间相差较多时或者裂缝扩展范围较大时，本文算法就没有做到很好地还原真实裂缝形态，在这种情形下就应该适当增加局部二值化区域的大小，实现裂缝目标在宽度上的完整提取。