岩石薄片正交偏光序列图像的颗粒分割

李周，滕奇志＊，张余强

（1.四川大学电子信息学院图像信息研究所，成都610065；2.成都西图科技公司，成都 610024）

0 引言

在石油地质行业中，确定岩石中各矿物成分比例是薄片矿物鉴定的重要内容[1]，对岩石颗粒进行准确分割则是确定矿物成分的前提。在岩石薄片图像中，颗粒表面纹理复杂，形状各异，颗粒间差异不明显，导致岩石颗粒图像自动分割的难度较大。

不少学者对岩石颗粒图像的自动分割进行了研究，如文献[2]中B.Obara将岩石图像转换到不同的颜色空间，然后再进行阈值分割。文献[3]中Gorsevski等使用了基于元胞自动机演化的边缘检测方法对岩石图像进行分割。虽然已有较多的颗粒分割算法研究，但是能够实际应用的却比较少。因为在自动分割之后，仍需要进行大量的人工交互，对岩石颗粒的分割结果进行修正，自动化程度不高。

文献[4]中提到，岩石样品中晶体存在消光特性，在不同正交偏光角度下，颗粒的亮度会发生变化，颗粒之间的边界也能清晰地显示出来，所以可以利用这种特性来辅助岩石颗粒的分割。以此为依据，本文通过对采集到的不同正交偏光角度下的岩石薄片序列图的分析处理，最终实现颗粒的自动分割。

1 偏光序列图采集处理

采集岩石薄片图像一般是通过在显微镜下利用投射单偏光来获取图像，但是由于单偏光下岩石薄片图像中的颗粒缺乏清晰的轮廓，颗粒之间没有良好的区分度，如图 1（a）所示。

本文使用正交偏光系统下采集的序列图进行图像分割。正交偏光系统是在单偏光系统的基础上再加上一个单偏光镜，上、下偏光镜振动方向相互垂直，将岩石薄片放在上、下偏光镜之间的载物台上，视域中将呈现一系列光学现象，如消光、干涉色等[6]。采用自主研制的偏光图像采集装置，每旋转正交偏光镜15°获取一幅图像，总共旋转120°，一共采集9张序列图像。图1（b）（c）（d）为其中 0°，60°，90°角度的 3 张图像，可以看出，正交偏光图像中的岩石颗粒清晰，但是由于颗粒的消光现象，在单张正交偏光图中并不能把所有颗粒都显示出来。

2 颗粒分割

图1 单偏光图像和不同角度下的正交偏光图像

利用一组正交偏光序列图像进行颗粒分割，比较直观的做法是对每一幅偏光序列图下的图像进行分割，然后将多幅图像的分割结果进行合并，得到最后的分割结果。但是通过实验发现，该做法存在这样的问题：由于岩石图像本身的复杂性和多样性，分割中存在许多误分割，对多幅图像的处理图进行叠加处理，容易将多幅图像的错误分割结果相叠加。因此本文未选择上述方法，而采用先将序列图像叠加到一起，再做分割处理的方法。

2.1 偏光序列图像的叠加

不同的颗粒在正交偏光镜下的明暗变化的趋势不同，这种趋势如果相同，则可以合理的认为是一个颗粒的概率较大[3]，本文将偏光镜下采集的不同序列图像之间的变化差值累加后得到叠加图像，如式（1）、式（2）所示：

2.2 标记分水岭分割

本文对叠加后的图像采用标记分水岭方法进行颗粒分割，分水岭分割[7]是一类目前使用最为广泛的颗粒形态学分割方法。其中H-minima[8]是一种常见的标记选择方法，采用梯度图像中的区域极小值作为种子点[9]，实现目标分割。本文根据文献[10]选取方向测度的极小点作为标记种子点，将每个像素点的邻域扩充到7×7，有效抑制了因为噪声或者颗粒形状不规则等因素而产生的多个虚假种子点，在一定程度上解决了基于梯度标记分水岭分割方法产生的严重过分割现象。

如图3为采用分水岭分割后用随机颜色填充颗粒，并将背景色置为黑色后的图，由此可以清晰看出，基本所有的颗粒都被提取出来。但在部分区域，如图3中的标记区域，仍存在过分割的现象。

图2 正交偏光下采集的序列图像的叠加图像

图3 分水岭分割后的图像

3 过分割区域聚类

图4 过分割的粘连区域

由图3中标记可见，对于相毗邻的颗粒，可能会出现将其中的颗粒过分割成为若干个小区域的现象。图4是将这部分区域放大后以一定透明度附加到原本的颗粒上的图像。如图4（a）所示，由于右上角的颗粒被过分割，导致原本只有两个颗粒的该区域被分割成多个区域。针对这种现象，本文首先对过分割的图像区域进行角度，面积，位置特征提取，然后采用“自底向上”的AGNES层次聚类算法进行聚类，将该区域聚合为符合实际情况的分割效果。

3.1 图像区域特征表示

以图5中的两个颗粒为例，标记上方颗粒为Grain-U,下方颗粒为Grain-D,其中Grain-U由于过分割导致被分成若干小的区域。

图5 被过分割的两个颗粒

本文从以下几个方面来描述区域的特征：

（1）区域的角度特征

由于岩石颗粒正交偏光的周期消光性，正交偏光镜下旋转360°，颗粒会出现4次消光现象，每90°为一个周期[11]。如图 6（a）（b）（c）（d）：分别是在正交偏光0°，45°，60°，105°下的颗粒图像，可见上下两个颗粒在不同偏光角度下呈现出来的亮度是不一样的，每个区域在一个周期内都存在亮度的最小值和最大值，且不同颗粒的亮度峰值所在角度不同。图6（e）是Grian-U中的区域亮度变化折线图，由于共同的消光性，虽然被过分割为多个不同的区域，但其亮度变化的趋势一样，均在75°达到峰值；而未被过分割的Grian-D颗粒在60°时亮度达到峰值。

因此可以用该峰值的角度来表征该区域的角度特征，如式（3）所示：

其中θLmax表示该区域亮度达到峰值时的角度，以π/2为底对数据进行标准化操作，θ'为标准化之后的角度值。

（2）区域的面积特征

由于分水岭的过分割现象，部分颗粒被分割成了许多小的区域。而这些小的区域的面积相对于其毗邻颗粒的面积较小，具有相似性，因此面积可以作为聚类的一个特征，如式（4）所示：

S'代表该区域的相对面积的大小，Ssum表示过分割区域的总面积，S为当前区域的面积。

（3）区域的位置特征

由图5（b）可以看出，Grian-U过分割区域的重心坐标相对集中，而Grian-D的重心坐标距离其他区域的相对较远。因此可以用区域的重心点坐标来表征该区域所在的位置，如式（5）所示：

(X',Y')代表当前区域重心坐标的相对位置。其中W、H分别表示图像的宽和高，(X,Y)为当前区域的重心坐标。

图6 同一颗粒不同角度的明暗变化

根据公式（3）（4）（5）计算后的特征量值都被限定到[0,1]之间，从而组成一个四维的特征向量V={θ',S',X',Y'}.

3.2 过分割区域聚类

本文采用“自底向上”AGNES层次聚类算法，在不同层次对数据进行划分，从而形成树形的聚类结构[12]。先将过分割区域内所有的单个区域看成一个初始聚类簇，然后找出距离最近的两个簇进行合并，该过程不断重复，直至达到终止聚类个数。本文计算簇间距离使用average linkage距离，簇间距离等于两组对象之间的平均距离。给定聚类簇Ci和Cj，其距离定义为：

其中d(Ci,Cj)表示不同簇Ci和Cj之间的距离，dist(Vp,Vq)是两个区域p和q之间的欧氏加权距离，Vpi表示这个区域在i分量的特征值，ωi表示在i分量的权值。经大量试验，区域的角度特征和面积特征在聚类作用中比重较大，取亮度权重和面积权重ω1=ω2=5；坐标位置权重ω3=ω4=3。

为了得到符合实际的聚类个数，在使用层次聚类的过程中，每一次两个相近簇合并后，需要对该次合并的有效性进行判断，即衡量当前聚类的有效性。首先定义层次聚类中的成本函数为：

其中，k是簇的个数，x表示簇中的数据，Ti是第i个簇中元素的集合；Ci表示第i个簇的中心，m为当前簇中的元素个数。成本函数是各个簇畸变程度之和，若簇内部的成员彼此间越紧凑则畸变程度越小，反之，若簇内部的成员彼此间越分散则畸变程度越大。

根据文献[13]，可以通过肘部法则来估计聚类数量。肘部法会把不同k值的成本函数值反映在曲线上。随着k值的增大，每个簇中包含的样本数会减少，样本离其中心会更近，因此平均畸变程度会减小。但是，随着k值继续增大，平均畸变程度的改善效果会不断减低。当Jk减少缓慢时，就认为即使进一步增大聚类数，畸变程度的改善效果也并不能增强，则此时的这个“肘点”对应的k值就是最佳聚类数目，如图7（a）。上述出现明显“肘部”的情况比较常见，但是Jk值的变化也可能比较平缓没有明显的拐点，如图7（b），这时可以认为该区域最终聚为一类。

图7 J随k值的变化曲线

本文设定当聚类数小于设定数kmax时，即当k取k=1，2，…，kmax，计算其成本函数。根据计算的结果得到Jk值变化对应的k值，则此k值就为最终的聚类个数，如果随着聚类数目的减小，Jk的变化幅度始终相当，变化曲线上没有明显的拐点,则把这个区域归为一类。

聚类过程为：首先设定一个过分割区域内的所有区域各为一簇，根据公式（6）计算所有簇间距离，将距离最近的两个簇合并，然后更新整个簇的个数；再计算新的簇之间的距离，根据距离最近原则再合并两个簇，如此反复。当簇个数小于等于设定值kmax后，开始计算其成本函数 Jk，并保存每次更新的簇的结果，直到最终聚为一类。

每合并一簇后计算每次变化的幅度：

计算ηk的平均值kavg，比较每次聚类后当前幅度值与平均值的大小，若满足ηk≥δ×kavg，则判断其为“肘部”，对应的k值即为最佳聚类个数，k值对应的聚类结果即为最终结果。如果均不满足，即不存在明显“肘部”，此时k取1，聚为一类。经过试验，当设定δ=2时，聚类结果较符合实际情况。

将同一类区域填充为同一种颜色，根据上述算法，图5的最终聚类图如图8所示：

图8 Grain-U和Grain-D最终聚类图

再分别对图4中的过分割区域进行聚类，结果如图9所示，可以看出，本文算法成功将过分区域聚类为符合实际情况的结果。

4 实验结果与分析

本文对二组实际拍摄的偏光序列图进行颗粒分割，并与文献[6]中的颗粒分割方法算法进行对比。

颗粒分割结果如图10所示，其中：（a）（d）分别为两组偏光序列图进行差分叠加后的图像；（b）（e）分别为通过文献[6]中SRM算法分割后的效果图；（c）（f）分别为本文方法对颗粒的分割效果图。文献[6]中先对每一张图片进行分割后再合并，提取到的颗粒比较多，但是过分割现象严重。本文中先将序列图合并后再进行分割，并且对过分割区域进行聚类合并，过分割现象得到很好的抑制，使得结果较于文献[6]的方法更符合实际需求。

5 结语

本文讨论了一种正交偏光图像序列图的颗粒分割算法，首先采集不同偏光角度下的序列图，对序列图进行叠加融合后再基于分水岭算法分割将目标颗粒提取；然后对过分割区域提取角度、面积、位置等特征；最后通过AGNES聚类算法将过分割的颗粒聚类到一起，从而完成颗粒的分割提取。

图9 过分区域聚类后效果图

图10 文献[3]算法以及本文算法对比效果图

实验结果显示，本文算法通过先分离后聚合的方式，成功将过分割的区域融合。但对于没有在分水岭阶段中没有提取到的颗粒，本文算法不能将其分割出来，可在这方面继续加以研究。

[1]徐耀鉴.岩石学[M].北京:地质出版社,2007.

[2]Obara B.A New Algorithm Using Image Colour System Transformation for Rock Grain Segmentation[J].Mineralogy and Petrology,2007,91（3-4）:271-285.

[3]Gorsevski P V,Onasch C M,Farver J R,et al.Detecting Grain Boundaries in Deformed Rocks Using a Cellular Automata Approach[J].Computers&Geosciences,2012,42（3）:136-142.

[4]杨宗瑞.基于偏光序列图像的岩石颗粒分割及矿物分[D].四川大学,2014

[5]Gorsevski P V,Onasch C M,Farver J R,et al.Detecting Grain Boundaries in Deformed Rocks Using a Cellular Automata Approach[J].Computers&Geosciences,2012,42（3）:136-142.

[6]吴拥,苏桂芬,滕奇志,等.岩石薄片正交偏光图像的颗粒分割方法[J].科学技术与工程,2013,13（31）:9201-9206.

[7]Beucher S,Meyer F.Segmentation:The Watershed Transformation.Mathematical Morphology in Image Processing[J].Optical Engineering,1993,34（5）:433-481.

[8]Chanho J,Changick K.Segmenting Clustered Nuclei Using H-minima Transform-based Marker Extraction and Contour Parameterization[J].IEEE Transactions on Bio-medical Engineering,2010,57（10）:2600-4.

[9]Li B,Pan M,Wu Z.An Improved Segmentation of High Spatial Resolution Remote Sensing Image Using Marker-based Watershed Algorithm[C].International Conference on Geoinformatics.IEEE,2012:1-5.

[10]杨海军,梁德群,毕胜.基于图象方向性信息测度的图象象素分类[J].中国图象图形学报,2001,6（5）:429-433.

[11]Fueten F,Goodchild J S.Quartz C-axes Orientation Determination Using The Rotating Polarizer Microscope[J].Journal of Structural Geology,2001,23（6-7）:895-902.

[12]周志华.机器学习[M].清华大学出版社,2016.

[13]Bholowalia P,Kumar A.EBK-Means:A Clustering Technique based on Elbow Method and K-means in WSN[J].International Journal of Computer Applications,2014.