图像处理技术在岩矿薄片鉴定分析中的应用已经成为地质行业的常用手段. 矿物颗粒分割是薄片鉴定和矿物成分识别的基础, 其效果直接影响后续鉴定分析的准确性. 利用岩矿薄片在固定视域下的正交偏光序列图像进行矿物颗粒的提取是一种新的研究方法, 已取得一定成果. 在实际生产应用中主要采用基于灰度阈值的分割算法^[1]、基于统计区域合并(SRM)的算法^[2]以及基于边缘流^[3]的算法.

超像素分割算法^[4]的最大优势在于边缘精准定位以及算法复杂度低. 本文使用基于熵率的超像素算法^[5,6]对偏光序列融合图像进行初步分割, 通过快速区域合并算法^[7,8]合并具有相似性的过分割区域, 同时, 基于矿物颗粒在正交偏光序列图像中的变化规律, 对结果进行二次合并, 改善超像素算法过分割严重的问题, 实现了对岩矿薄片图像中颗粒的分割提取. 本文分割算法的流程图如图1所示.

1 基于偏光序列图像的超像素分割

超像素^[9]是指具备相似颜色、灰度、纹理等特征的相邻像素组成的图像块, 超像素分割实质上就是依据各种特征将相似或相同的像素聚集到一起的聚类问题. 超像素分割算法在国内外已广泛应用, 主要分为基于图论和梯度下降的两类方法, 本文采用基于图论的熵率超像素分割算法对岩矿薄片偏光序列图进行矿物颗粒分割.

1.1 熵率超像素分割

熵率法超像素^[5]提出了将图上随机游走熵率和平衡项相结合的目标函数, 并在聚类过程中不断优化目标函数实现对图像的分割. 该方法分割出的超像素块相对规则、均匀, 其主要思想分为以下四个方面.

(1)图的构造: 采取映射的方式, 将图片映射至无向图 $G = (V,E)$ , 其中 $V$ 表示像素点集, $E$ 为边集, 并将像素间的相似性定义为边的权值 $\omega $ . 目标是寻找到集合 $A \subseteq E$ , 使生成的无向图 $G = (V,A)$ 中包括 $K$ 个连通的子图.

(2)随机游走熵率: 熵率用来度量随机过程的不确定性, 可以很好的表征信息的冗余程度、内部结构的紧凑性, 随机游走模型的加入能够获取紧凑、均匀的超像素区域块. 无向图 $G = (V,A)$ 上的随机游走熵率计算公式如下:  

$H'(A) = - \sum\limits_i {{\mu _i}} \sum\limits_j {{p_{i,j}}(A)\log ({p_{i,j}}(A))} $

(1)

其中, ${\mu _i} = {\omega _i}/{\omega _T}$ , 且 ${\omega _T}{\rm{ = }}\sum\nolimits_{i = 1}^{\left| V \right|} {{\omega _i}} $ , ${\omega _i}$ 为连接节点 $i$ 所有边的权值和; ${p_{i,j}}$ 为随机游走模型转移概率, 其定义如式(2)所示, 式(2)中 ${e_{i,j}}$ 为连接点i、j间的边, ${\omega _{i,j}}$ 为边 ${e_{i,j}}$ 的权值.

${p_{i,j}}(A) = \left\{ {\begin{array}{*{20}{l}} \frac{{{\omega _{i,j}}}}{{{\omega _i}}}& if\; i \ne j\; and\; {e_{i,j}} \in A \\ 0& if\; i \ne j\; and \;{e_{i,j}} \notin A \\ 1 - \frac{{\sum {_{i:{e_{i,j}} \in A}{\omega _{i,j}}} }}{{{\omega _i}}} &if \;i = j \\ \end{array}} \right.$

(2)

(3)平衡条件: 在目标函数中添加平衡条件, 可以约束限制熵率的随机性. 令 $A$ 为已选择边集, 划分出 $A$ 中超像素块集为 ${S_A} = \left\{ {{S_1},{S_2}, \cdots ,{S_{{N_A}}}} \right\}$ , ${N_A}$ 为 $G = (V,A)$ 中连通子集的数目, ${Z_A}$ 为子图聚类的概率分布:

${p_{{Z_A}}}(i) = \frac{{\left| {{S_i}} \right|}}{{\left| V \right|}},i = \{ 1,2, \cdots ,{N_A}\} $

(3)

由此, 定义平衡函数为:

$B(A) = H({Z_A}) - {N_A} = - \sum\limits_i {{p_{{Z_A}}}(i)\log ({p_{{Z_A}}}(i)) - {N_A}} $

(4)

其中, ${Z_A}$ 的熵 $H({Z_A})$ 可以使聚类结果有相似的尺寸, 同时 ${N_A}$ 约束集群的数目. 图2中用不同粗细的线条连接代表不同的集群, 可见图2(a)中的平衡项大, 集群分布更均匀.

(4)目标函数: 组合随机游走熵率和平衡条件, 得到聚类最优化的目标函数如式(5)所示, 式中集合 $A$ 满足 $A \subseteq E$ , 且 $\lambda $ 为正实数. 聚类过程中通过最大化目标函数得到图像的最佳分割结果.

$\mathop {\max }\limits_A [H'(A) + \lambda B(A)]$

(5)

1.2 基于偏光序列图像的熵率超像素分割

岩矿薄片在单偏光下获取的图像中, 矿物颗粒之间粘连比较严重、缺少清晰的边界, 当颗粒颜色相差不大时, 很难辨认出具体是几个颗粒, 因此, 给颗粒分割带来极大的困难. 在正交偏光镜下, 矿物颗粒在视域下会表现出一系列特有的光学特性^[10], 如消光特性、干涉色等. 正交偏光下的消光特性^[10]在非均质矿物中普遍存在, 当正交偏光镜旋转一周, 矿物共会出现四次消光现象, 并且, 不同矿物颗粒在不同偏光角度下消光状态也可能不一样, 加上矿物颗粒在正交偏光下的干涉色, 由此显现出明显的边缘. 本文所需的序列图像由自主研制的偏光图像采集与拼接系统获取, 先在单偏光下采集视域的一张图像, 然后保持显微镜载物台固定, 每旋转10°正交偏光镜采集一张图像, 共旋转120°, 共采集13张固定视域下的正交偏光图像, 图3为部分序列图像.

考虑到矿物颗粒本身的复杂性, 以及单个正交偏光角度下矿物颗粒的不完整性给分割带来的影响, 首先选取一个消光周期(90°)内的正交偏光序列图像, 即拍摄的前十张正交偏光序列图像, 依据式(6)进行图像差分^[11]融合, 形成一张拥有视域中全部颗粒的正交偏光图像, 如图4(a)所示.

$g(x,y) = \left[ \begin{gathered} {g_R}(x,y) \\ {g_G}(x,y) \\ {g_B}(x,y) \\ \end{gathered} \right] = \sum\limits_{i = 1}^9 {\left[ \begin{gathered} {f_{(i + 1)R}}(x,y) - {f_{iR}}(x,y) \\ {f_{(i + 1)G}}(x,y) - {f_{iG}}(x,y) \\ {f_{(i + 1)B}}(x,y) - {f_{iB}}(x,y) \\ \end{gathered} \right]} /9$

(6)

其中, $g(x,y)$ 为融合后图像, ${f_i}(x,y)$ 是序列图中第 $i$ 张正交偏光图像, 且 $g(x,y)$ 的R、G、B分量 $ {g_R}(x,y)$ 、 ${g_G}(x,y)$ 、 ${g_B}(x,y)$ 均已映射至 $0 \sim 255$ .

为了减弱R、G、B分量间的相互影响, 提高矿物颗粒分割的效果, 本文将融合后的图像 $g(x,y)$ 转化至Lab空间. 对Lab空间下的融合图像采用1.1节中的算法进行图像分割, 其分割结果如图4(b)所示.

2 矿物颗粒目标的区域合并 2.1 基于融合图像与单偏光图像的快速区域合并

由图4(b)的结果可以看出, 利用超像素分割算法可以完整的保留矿物颗粒的边缘信息, 但是使矿物颗粒出现过分割现象. 一般而言, 属于同一个矿物颗粒的区域具有相似性, 利用其相似性将同质区域合并. 本文首先采用文献[7]中提出的基于RAG(区域邻接图)和NNG(最近邻接图)数据结构的快速区域合并算法对同质区域进行合并, 其中合并判定准则结合了目标区域融合图像的差分直方图和单偏光图像的灰度直方图信息.

(1)相似性合并准则的确定

① 融合图像的差分直方图: 差分直方图反映了正交偏光下矿物颗粒图像在一定方向、一定距离上相邻点灰度之间差的概率分布^[12,13]. 在正交偏光的融合图像中, 矿物颗粒本身的形态、纹理比较复杂, 较颜色直方图而言, 差分直方图具有颗粒目标区域图像的综合纹理信息, 因此, 提取偏光序列融合图像的差分直方图作为相似度主要判定准则.

设一幅具有 $M(M = 1,2, \cdots ,L)$ 灰度级的图像 $F$ , 区域 $D$ 内的两个像素 ${f_1}{\rm{ = }}f(x,y)$ 、 ${f_2} = f(x + {d_x},y + {d_y})$ , $({d_x},{d_y}) \in D$ , 对于偏移 $({d_x},{d_y})$ 的像素差分按式(7)计算如下:

$d(x,y) = {f_1} - {f_2}$

(7)

由此可得归一化的差分直方图如下:

${P_d}(i) = {h_d}(i)/N,i = - L, - L + 1, \cdots ,0,1, \cdots ,L$

(8)

式中, ${h_d}(i)$ 为区域 $D$ 中 $d(x,y)$ 的分布, $N$ 为区域 $D$ 的像素总数. 在计算每一个区域的差分直方图矩阵时, 本文在每个像素的8邻域中选取四个方向 $({d_x},{d_y}) = (d,0),(d,d), $ $(0,d),( - d,d)$ 进行合并计算.

统计学中, 可以用Bhattacharyya系数^[14]计算两个统计样本的重叠量, 从而进行两组样本的相关性测量, 其计算公式如式(9)所示. 在进行直方图的相似性计算时, 该计算方法可以得到最好效果.

$BC(p,q) = \sum\limits_{x \in X} {\sqrt {p(x)q(x)} } $

(9)

式中, p(x)、q(x)为两组不同的概率分布, 且 $0 < BC < 1$ .

因此, 根据式(9), 定义融合图像差分直方图的相似度合并准则 ${h_C}$ 为:

${h_C} = \sum\limits_{i = {\rm{ - }}L}^L {\sqrt {{P_{dp}}(i) \times {P_{dq}}(i)} } $

(10)

其中, P_dp、P_dq分别为两不同区域p、q的差分直方图.

② 单偏光灰度直方图: 单偏光图像中, 同一矿物颗粒所在区域的灰度值差异变化不大, 可将区域的灰度值作为次要判定准则. 直方图统计是提取区域灰度特征的常用方法, 根据式(9)可得单偏光下灰度相似性合并准则 ${h_A}$ 为:

${h_A} = \sum\limits_{i = 0}^{255} {\sqrt {{H_p}(i) \times {H_q}(i)} } $

(11)

其中, H_p、H_q分别为两不同区域p、q的灰度直方图.

综合①、②定义本文方法中使用的相似性合并准则如式(12)所示:

$h = \alpha {h_C} + (1 - \alpha ){h_A}$

(12)

(2)快速区域合并算法

区域邻接图(RAG)的数据结构可以表示图像中全部区域之间的邻近关系, 利用该数据结构并依据(1)中讨论的相似性准则实现区域合并, 同时, 使用最近邻接图(NNG)实现基于RAG区域合并算法的加速优化. 具体操作过程如下:

① 初始化数据结构: 计算所有超像素区域下对应融合图像的差分直方图以及单偏光图像的灰度直方图. 将全部的超像素区域作为单独的节点, 按照区域邻接图(RAG)的规则生成一个无向邻接图, 并据此生成对应的最近邻接图(NNG);

② 搜索邻接情况: 采用基于NNG的搜索方法, 对于每一个区域确定其所有的邻接区域;

③ 合并: 依据(1)中的相似性合并准则计算两邻接区域的相似性, 若达到阈值, 将两邻接区域合并成一个区域;

④ 更新数据结构: 合并之后计算该区域相应的差分直方图以及灰度直方图数据, 并更新RAG和NNG;

⑤ 重复②~④, 直至没有区域需要合并为止.

图5为合并前后对比图. 可以看出, 使用该方法基本可以实现大部分矿物颗粒的过分割合并, 但由于矿物颗粒本身的复杂性, 仍有部分矿物颗粒的过分割区域未实现合并. 本文在此结果基础上, 利用矿物颗粒在正交偏光序列图像中的变化特性对未实现合并的矿物颗粒再次合并, 最终达到较为理想的结果.

2.2 基于正交偏光序列图像的二次区域合并 2.2.1 颗粒目标图层的初步提取

采用熵率超像素算法可以得到良好的颗粒边界, 其中碎屑岩颗粒的边界也较为理想. 由于最终需要得到颗粒目标图形层进行粒度分析、成分识别, 因此, 利用矿物颗粒在序列图中的亮度特性采用如下步骤初步获取颗粒目标图层:

① 对2.1节的边缘结果进行目标反转操作得到全部目标区域;

② 依据式(13)对①中的所有区域进行消光性筛选, 将亮度差低于阈值 $T$ 的区域置为背景色, 获得矿物颗粒区域;

${I_{\max }} - {I_{\min }} > T$

(13)

其中, I_max、I_min分别为目标区域在正交偏光的前后序列图中的最大亮度均值与最小亮度均值, T为亮度筛选阈值条件. 经大量实验, T取25筛选效果最佳;

③ 对②中的结果进行去噪、填孔、数字形态学、边界平滑等操作.

根据上述步骤初步得到矿物颗粒目标的图形层图像如图6所示.

2.2.2 颗粒目标图层的消光性合并

经过区域亮度筛选可基本获取颗粒目标区域, 但仍存在部分颗粒被分为若干个, 因此, 基于图6的结果, 利用矿物颗粒目标在序列图中的消光变化特性继续对仍存在的过分割颗粒目标区域进行相似性合并.

对于图7(a)中的完整颗粒区域W以及其在图7(b)中的过分割区域A、B, 每个目标区域在正交偏光序列图像下的亮度均值变化趋势如图8所示.

虽然矿物颗粒出现了过分割现象, 由于其存在共同的消光特性, 亮度均值在序列图像中具有相同的变化趋势. 在统计学中, 相关系数用来描述两个序列之间的相关程度, 设两个亮度均值序列为 $a = ({a_1},{a_2}, \cdots ,{a_N}), $ $b = ({b_1},{b_2}, \cdots ,{b_N})$ , 则以式(14)计算序列的a、b的相关系数 ${\beta _{ab}}$ :

${\beta _{ab}} = \frac{{\operatorname{cov} (a,b)}}{{{\sigma _a}{\sigma _b}}}$

(14)

其中, $\operatorname{cov} (a,b)$ 为两个序列的协方差, σ_a、σ_b分别为序列a、b的标准差.

因此, 本文以两区域在正交偏光序列图下亮度均值序列的相关系数作为区域消光特性的合并准则, 当两区域的相关系数 $\beta $ 小于阈值 ${\beta _T}$ 时, 判定两区域为同一颗粒目标, 将两目标区域进行合并. 对初步提取的颗粒图层中所有颗粒目标进行合并预测, 直至无颗粒目标需要合并为止, 得到最终的分割结果.

3 实验结果分析

本文基于熵率超像素算法, 针对矿物颗粒在正交偏光前后序列图像中的变化特性完成矿物颗粒图层的提取. 本文实验程序使用C++语言, 在VS2015平台下基于MFC框架编写实现.

经过大量薄片图像的测试, 本文的分割算法能够取得较好的分割效果. 以图3中矿物薄片的原始偏光序列图像为例, 将本文算法与文献[2]、文献[3]中的矿物颗粒分割算法的矿物颗粒提取结果进行对比分析, 如图9所示. 其中, 图9(b)为本文算法得到的颗粒图层, 图9(c)、图9(d) 分别为SRM和边缘流算法得到的颗粒目标图层.

结合该薄片的原始偏光序列图像(图3)以及原始融合图像(图9(a))可以看出, 文献[2]中利用SRM算法得到的分割结果容易出现大量的欠分割现象, 使很多矿物颗粒目标粘连在一起, 并且一些矿物颗粒的边缘定位不够准确; 文献[3]中利用边缘流算法得到的分割结果大部分比较理想, 由于边缘流算法对边界定位的敏感性, 对于表面纹理特征较为复杂的碱性长石易出现过分割情况且不易融合, 使提取区域不够完整, 如图9(d)中标注出的颗粒目标. 本文的分割算法可以较为完整地提取视域中的矿物颗粒, 改善了碱性长石的过分割现象, 同时, 对于粒径较大的碎屑岩颗粒也有了较为完整的提取结果, 符合实际生产需求.

4 结论

本文针对固定视域下岩矿薄片的偏光序列图像, 提出一种矿物颗粒目标提取分割的方法: 首先将岩矿薄片在正交偏光下的序列图像进行差分融合后再基于熵率超像素算法分割矿物颗粒, 然后提取分割区域在融合图像下的差分直方图以及单偏光图像下的灰度直方图作为相似性特征, 对过分割区域进行快速区域合并, 最后通过矿物颗粒在正交偏光下的消光特性对分割区域筛选后再次进行颗粒融合, 从而得到最终的颗粒分割结果.

实验结果表明, 该算法综合了视域下偏光序列图的特征信息, 能够较为准确的提取矿物颗粒目标, 并且, 对碎屑岩颗粒的提取效果有一定的提升, 为矿物颗粒的分割提取提供了一个新的方向.