涂秉宇，熊淑华，滕奇志，田刚

（1.四川大学电子信息学院图像信息研究所，成都 610065；2.新疆油田公司行政事务中心，克拉玛依 834000）

0 引言

在储集层岩石微观结构中，由三个或三个以上岩石颗粒所包围的空间被称为孔腔，相邻两个孔腔之间的最窄连接部分称为喉道，孔腔和连接它的部分喉道的总体称为孔隙[1]。喉道是相邻岩石孔隙之间最窄的部分，其大小和分布都是预测多孔介质的渗流性的重要因素[2]。因此对于喉道的研究是岩石孔隙结构分析的重要一环。

目前，微观孔隙结构的分析过程是：利用铸体薄片在光学显微镜下选取典型视域成像，再在获取的图像上根据铸体薄片中孔隙的颜色特性将孔隙部分作为目标提取出来，然后对孔隙部分进行喉道的划分。而喉道一般是处在目标区域凹陷或凸出最窄处，因此，对于喉道的分割可以通过对目标图像拐点处的检测、筛选与匹配来实现。最后再根据石油天然气行业标准[1]进行参数计算，从而得出微观孔隙结构的数据。

本文提出了一种基于拐点检测和匹配的自动喉道分割方法，该方法通过检测拐点信息，得到拐点的位置，再通过两两匹配拐点，形成拐点对，以此来确定喉道位置，从而实现对岩石孔隙连通区域的分割。

1 喉道分割算法

对于喉道分割给出以下示例，图1是一幅孔隙原图及其分割结果图，（a）是孔隙原图，（b）是孔隙提取结果图，（c）是将孔隙部分叠加到原图的效果图，（d）是喉道分割结果图，其中红色分割线为喉道。由此可以看出，喉道分割的本质在图像处理中的意义是找出孔腔与孔腔之间最窄的部分。由于喉道往往处在边界的凹陷处，如果对目标边界进行拐点检测，再根据拐点的位置便能准确找出边界凹陷处，基于此就能准确地进行喉道分割。所以本文提出了基于拐点的喉道分割算法，该算法首先通过八邻域法获得孔隙二值图的边界，引入边界链码代替曲率，通过链码的变化找出孔隙边界上的拐点。对找出的拐点进行筛选，过滤掉伪拐点，得出有效拐点。然后根据岩石孔隙结构的物理特性对有效拐点两两匹配，并处理拐点数为奇数的情况，从而形成喉道，最终得出喉道分割位置。

图1 孔隙示例图

1.1 获得目标区域的边界

孔隙提取后的效果图一般为二值图，而对于二值图目标区域的边界判定可以使用八邻域连通标记法。所以本文对孔隙区域采用八邻域法进行标记判别，具体步骤为：首先针对孔隙图片的提取结果进行处理，其顺序是从左上方的第一个像素点开始，从左到右，从上到下依次遍历，并按照图2所示掩膜的方法对孔隙区域进行标记，直至整幅图搜索结束。最后根据八邻域判别边界，对标记的孔隙区域的所有像素点进行迭代，判断其邻域是否属于该标记区域，若其不属于，则断定为该区域的边界，否则判断为该区域所属的点。

图2 八邻域掩膜

1.2 边界曲线的拐点检测

（1）曲线的点曲率

边界曲线进行拐点检测需对曲线上的点进行曲率计算，然后根据曲率判断拐点位置。

如图3所示，假设函数y=f（x）代表的是一条曲线，P1和P2是曲线上任意两点，曲线上P1点的曲率k定义为点沿着曲线的切线方向与水平X轴的夹角的变化[4-5]。

α1为点P1的切线夹角，α2为点P2的切线夹角，△l为曲线P1P2的弧长，△α为切线变化的角度。

图3 曲率示意图

（2）链码

由曲率的定义可知，其计算需要求极限和乘除运算，相对较复杂[6]。本文引入计算更为简便的链码来表示边界曲率，避免乘除法运算，可以有效地节省计算时间。

所谓链码是用来反映像素点和其邻近像素点方向的一种编码[7-8]，其编码方式如图4，用当前像素点指向它的八邻域方向来表示。

图4 八连通域链码

链码可以表示边界的斜率、曲率，即当前边界点的切线方向。链码差则是和曲率成正比的量[3]，其计算公式为：

式中D（xi）为边缘曲线上第i个像素点的链码值。但是由于只有0～7八个方向，细化程度差，导致精度缺失。文献[9]算法根据链码特性利用其计算得到的曲率进行局部平均从而得以提升算法效果。而文献[10]引入了N点链码差，利用进入和离开点Xi的N个链码平均值之差，提升算法精确度。以上算法都只是针对八邻域进行方向编码，而八邻域方向编码将2π平面分为8个方向，每个方向精度的偏差达到π/4，这存在量化精度不足的问题。因此本文采用了K邻域链码中的16邻域链码方法[6]。16邻域链码在8邻域链码的基础上向外扩展一圈，将2π平面分为16个方向，这使得16邻域链码相较于8邻域链码可以在不增加过多运算量的情况下提升方向精度和增加边界平滑[10]，这对于边界走势复杂的孔隙区域而言具有实用性和可行性。

（3）拐点检测算法

本文基于链码的拐点检测算法如下：

①对提取得到的边缘轮廓曲线进行16邻域链码编码。

②按照曲率定义，使用差分替代微分，则边缘曲线第i点的切线角度变化的差分表示值θi，可通过公式（3）求得：

其中，C（16,i）表示第i个像素的16邻域链码值。

由于θi是拐角的差分表示，所以应该将其角度归一化到[0,π]。而在16邻域链码编码中，[0,π]的区间对应的链码区间为[0,8]，所以点Pi处的曲率计算为：

为了进一步扩大拐点和非拐点处曲率的差距，所以使用以下公式，得到最终曲率值：

③根据曲率值得到曲率局部峰值的位置。即确定一个阈值 t，当 ei≥t>ei-1时曲率曲线为上升，当 ei>t≥ei+1时，曲线为下降，那么在此范围内与波峰位置左右相对应的边缘点即是待选拐点。其中t值应为边界曲率值的平均值，计算公式如下：

④判断凹凸性。以相距拐点前后相邻n点的中点为判别条件，如果中点像素值属于目标区域则为凸点，否则为凹点。

孔隙连通区域内部存在内部孔洞时，内部孔洞目标的边界为内轮廓，孔隙连通区域边界为外轮廓。对内轮廓取其凸点，经由以下筛选和匹配等步骤后与外轮廓上的凹点形成喉道分割线，用以进行有效的喉道分割。

1.3 筛选算法

由于孔隙区域边界凹陷处呈块状，而凹陷块中曲率的变化会导致检测出多个拐点，同时检测出的拐点中也可能有不属于喉道位置的噪声点，所以经由拐点检测后得出的拐点仅是待选拐点。还应对待选拐点进行筛选，去除相邻点和噪声点的干扰。为了达到消除噪声和冗余点的目的，本文采用的筛选算法如下：

（1）对于孔隙连通区域目标边缘凹陷块的情况，根据孔隙形状特征，本文通过某个距离阈值来判断待选凹点是否相邻，在阈值范围内的待选凹点属于相邻凹点。相邻的待选凹点，将其看作一个凹点群，取其中曲率最大的点作为该组凹点的有效凹点，代表这组凹点。而对于孔隙连通区域内部的孔洞，边缘检测形成的内轮廓上会有凸起块，通过拐点检测后找到的相邻凸点即为一个凸点群，取其中曲率最大的点作为有效凸点。因此，检测到的前后两个拐点之间的距离小于某个阈值的时候，我们认为这两个拐点是相邻的，应舍弃曲率较小的点。通过实验发现，两点间距离阈值一般设定为4～6个像素。

（2）由于孔隙边界情况的复杂性，所以由第一步筛选之后所得的拐点，容易出现拐点比较“浅”，即如图5所示的情况：通过第一步筛选算法后，得出A为有效拐点，但是A到直线BC的距离却很短，即A到直线的距离值非常的小。根据喉道的定义可知，喉道为两孔腔之间最短的地方，所以如图5所示的点不应该计入有效拐点。

图5 由于太“浅”不能作为拐点

因此本文提出确定一个有效的d值对此类拐点进行判别。由于全薄片图像中孔隙区域的边缘情况复杂，如果d值的选取过小，会导致一些原本不是匹配点的拐点被误判；取值过大，则会导致漏掉一些真正的拐点。所以d的取值不能为固定值。本文算法采用了自适应d值，即使用所有孔隙目标根据如下公式计算得到平均面积大小：

最后根据该面积大小的等效圆半径确定适中的d值。等效圆半径r计算如式（8）：

（3）对于所有已得的拐点进行以上两步筛选迭代，直至所有拐点筛选完毕。

1.4 匹配算法

喉道处在孔隙边界两端凹陷处，而要找到相应喉道则必须对边界两端凹陷处的拐点进行配对，如图1（d）所示。因此通过筛选算法得出的有效拐点还需进行拐点对匹配算法，从而形成喉道线对孔隙连通区域进行有效的分割。

本文采用的匹配算法如下：

（1）喉道线的有效凹点对应该是位于孔隙连通区域两侧，如图6所示。

两个有效凹点待匹配的示意图，其中的E1和E2为有效凹点，D1、F1和 D2、F2分别是 E1点和 E2点的等距前继点与后继点。将两个待匹配凹点的连线向两端进行延伸，如果E1一侧的延伸线位于E1D1向量与E1F1向量形成的锐角区域内，E2一侧的延伸线位于E2D2向量与E2F2向量形成的锐角区域内，且E2位于E1D1和E1F1的反向延长线所覆盖的边界区域内。只有上述条件均满足的情况下，满足匹配条件，该组有效凹点是相对应匹配的，否则该组有效凹点不对应匹配。

（2）通过判别步骤1进行匹配后，如果出现当前待匹配的有效凹点对应多个待匹配点的情况。根据孔隙喉道的定义，喉道是孔腔与孔腔之间最短的部分，所以需引入距离判别进一步判别正确的凹点匹配对，即只要分割线穿越目标区域的欧氏距离短即可满足此条件，如图7所示。

图7 最短距离匹配

A1点对应候选的匹配点存在A2、B2两个点，计算A1A2和A1B2的欧氏距离，取距离最短的A2进行匹配，欧氏距离公式如下：

（3）由于喉道是最窄处，而根据图1中喉道的宽度和其所属孔隙连通区域的等效圆直径比较可知，分割线线长必定满足不大于某个D值。即D应满足公式（10）：

其中C为孔隙连通区域等效圆周长，D为两匹配凹点间距离。

（4）由于分割线必须位于孔隙内部，所以喉道分割线上的像素点应该全部属于孔隙连通区域。因此对匹配点对连线上的点进行抽样分析，如果存在像素点属于背景区域则舍弃该喉道，否则该连线可以作为喉道。

（5）对于孔隙连通区域内部存在孔洞目标的情况，应该将内轮廓上检测到的凸点进行筛选，再与外轮廓上的凹点根据以上四点匹配准则进行有效匹配，形成分割线，如图8所示。

图8 最短距离匹配

1.5 奇数拐点的处理方法

上文的匹配算法是在一个孔隙连通区域中恰好筛选出了偶数个有效拐点的情况下进行一一匹配的。但是当遇到奇数个拐点的情况时，会在孔隙局部区域留下未做匹配的有效拐点。因此，根据喉道定义，本文通过选取对角边缘延长线所覆盖的边界区域进行最短距离判断，找到距离有效拐点最近的边界点与此拐点进行匹配，如图9所示。

图9 匹配非拐点

反向延长线范围可以找出边界上距离A、C点最近的B、D点，同时通过上文的匹配算法可以有效的判断AB、CD是否能够成为分割线。

通过此法，我们可以处理拐点数为奇数的情况，也防止了欠分割的情况，使得基于拐点的喉道分割算法更具有实际应用的价值。

2 测试与分析

为了验证本文算法的可行性，我们采用实际生产中的铸体孔隙图片，对不同情况下的孔隙进行实验。图10分别给出了不同情况下的孔隙原图及其对应的分割效果图，左边是原图，右边是对应的分割效果图。

图10 实验结果图

其中（a）（b）为孔隙中有一个孔洞情况的示例。图（c）（d）为孔隙中有多个孔洞情况的示例。图（e）（f）为孔隙中无孔洞但是区域形状复杂的情况示例。从图11中可看出，本文算法能适应各种情况的喉道分割。

同时本文针对基于形态学流域分割目标方法和腐蚀膨胀分割目标方法的结果进行实验比对。图11给出了不同算法的分割结果对比图，其中（a）为原图。（b）为基于流域算法的分割结果，其能有效地分割出一部分喉道，但是由于其形成分水岭的过程中受到孔隙部分局部极值影响，从而导致过分割现象。（c）（d）（e）分别为基于腐蚀膨胀算法17、20、25次后的分割结果图，图（c）在17次的腐蚀膨胀后会出现欠分割现象。图（d）进行了20次腐蚀后会出现小颗粒被完全腐蚀掉，存在欠分割的现象。图（e）腐蚀膨胀25次导致粘连性较大的区域得不到分割，而不是喉道区域的位置被分割。（f）为本文算法分割结果图，可以看到孔隙中的喉道被全部分割出，没有过分割和欠分割的现象。因为基于拐点的喉道分割算法从喉道定义出发，在分割喉道时只考虑目标的边界信息，避免了孔隙区域局部极值和算法循环次数不同带来极大误差的影响，因此本文算法能有效地进行喉道分割。

图11 分割对比图

3 结语

本文提出了一种基于拐点的喉道分割算法，该算法主要通过检测出图像中孔隙连通区域的边界，然后对边界进行拐点检测，并通过筛选算法确定有效拐点。最后进行有效拐点对的匹配，从而形成喉道分割线。通过理论分析和实验表明，该算法具有以下特性：①算法具有很强的鲁棒性，可以分割不同形态下的孔隙目标，例如内部存在孔洞情况。②算法准确度高，无论在外轮廓的凹陷处还是内轮廓的凸出处，均能准确地找出喉道位置。③算法计算量小，通过引入链码表示曲率，可以有效减少计算量。