CN105354600B - 一种砂岩显微薄片的自动分类方法 - Google Patents

Abstract

本发明公开了一种砂岩显微薄片的自动分类方法，其步骤是：1)对砂岩的显微薄片图像进行预处理，分割矿物颗粒；2)以颗粒为单位分别提取纹理和形状特征，选择代表颗粒，构成薄片图像的特征向量；3)对目标地域的砂岩薄片，基于欧氏距离采用近邻法选择相似的砂岩薄片样本，构建训练集合；4)精化训练集数据，包括类间平衡和特征降维；5)训练随机森林分类器，预测砂岩薄片的种类。本方法运用图像处理技术和机器学习方法，自动鉴别来自不同地域的砂岩薄片种类，能够解决由于地域不同导致的砂岩薄片结构差异较大而难于互相借鉴的问题。

Description

一种砂岩显微薄片的自动分类方法

技术领域

本发明涉及不同地域砂岩显微薄片的自动分类问题，特别针对地域不同导致的砂岩薄片结构差异较大而难于互相借鉴，运用图像处理技术和机器学习方法，实现不同地域砂岩薄片的自动分类。

背景技术

砂岩是一种主要由砂粒胶结而成的陆源沉积岩，主要由石英和长石组成。按照其成分和成因可划分为：石英砂岩、长石砂岩和岩屑砂岩三大类。砂岩是一种使用最广泛的建筑用石材，还是石油和天然气的储集层，砂岩的分类在油气和矿物勘探等领域有重要意义。

砂岩薄片分类根据砂岩的纹理和形状特征，确定砂岩薄片的种类，进而确定其结构、来源和成因。现有砂岩薄片分类方法多采用人工，费时费力，且过度依赖个人经验，不可重复。另外，由于形成的地质地域或环境不同，不同地域砂岩的结构差异较大，给砂岩薄片的分类带来了很大困难。

发明内容

本发明所要解决的技术问题是提供一种砂岩显微薄片的自动分类方法，该方法运用图像处理技术和机器学习方法实现对不同地域砂岩显微薄片的自动分类。

为达到上述目的，本方法采用如下的步骤：

1)读取砂岩显微薄片图像，对其进行滤波、去噪、平滑等预处理；采用图像分割技术，对砂岩显微薄片图像进行区域分割，获得矿物颗粒图像；

2)基于砂岩显微薄片图像中的矿物颗粒，分别计算纹理和形状特征指标，构建薄片图像的特征向量；

3)对目标地域的砂岩薄片图像，计算其特征向量，采用近邻法选择若干个相似的有标注的砂岩薄片图像样本，构建训练集合；

4)精化训练集数据：采用随机取样方法对初始训练集中少数类的样本进行重采样，实现训练集数据的类间平衡；分别计算每个特征指标的ANOVAF值，设定百分比，排除不必要的特征指标；

5)训练随机森林分类器，以目标地域砂岩薄片图像为输入，自动预测其种类。

上述步骤2)中基于砂岩显微薄片图像的矿物颗粒构建薄片图像的特征向量的处理过程是：首先针对薄片图像中的每个矿物颗粒，分别计算纹理特征和形状特征指标，组成矿物颗粒的特征向量X；然后对所有矿物颗粒的特征向量加权平均，得到均值特征向量X_m，计算公式如下：

其中k为薄片图像中包含的矿物颗粒总数；再后选取前K个面积最大的矿物颗粒的特征向量X_i(i∈[1,K])；最后将X₁,X₂,...,X_K,X_m串接，构成薄片图像的特征向量。

上述步骤2)中计算矿物颗粒的纹理特征的处理过程是：首先将颗粒图像转换为灰度图像，计算位移矢量d分别为(1,0)、(1,1)、(0,1)以及(-1,1)的4个归一化灰度共生矩阵P：矩阵P是一个256×256的矩阵，其第i行第j列的矩阵元素p[i,j](0≤i≤255,0≤j≤255)代表按矢量d分离的像素对的个数除以颗粒图像的总像素数量，像素对中第一个像素的灰度值为i，另一个像素的灰度值为j；分别按照4个位移矢量计算，共得到4个灰度共生矩阵；

然后针对每一个灰度共生矩阵P，计算能量Energy、熵Entropy、对比度Contrast以及相关系数Correlation指标，计算公式如下，其中p[i,j](0≤i≤255,0≤j≤255)代表矩阵P中元素：

基于4个灰度共生矩阵，每个矿物颗粒获得16个纹理特征值。

上述步骤2)中计算矿物颗粒的形状特征的处理过程是：首先计算颗粒图像中包含的像素点总数，即颗粒的区域面积A；统计其边界像素点的个数，即区域周长C；然后求颗粒区域的凸包，指包含颗粒的最小凸多边形；计算凸包面积H，即凸包所包含的像素点总数；据以计算凸包差CD和密制度Solidity指标，计算公式如下：

CD＝H-A

Solidity＝A/H

接下来求包含颗粒区域的最小外接矩形，计算其面积ERA，即外接矩形包含的像素点数量；从而获得每个矿物颗粒的5个形状特征值：区域面积A、周长C、凸包差CD、密致度Solidity、以及外接矩形面积ERA；合并形状特征和纹理特征，每个矿物颗粒转换为包含21个特征值的特征向量X。

上述步骤3)中基于目标地域砂岩薄片图像的特征向量，采用近邻法构建训练集的计算过程是：对目标地域的砂岩薄片图像，首先计算其特征向量与现有样本集中各砂岩薄片图像特征向量的欧氏距离；给定两个特征向量X₁＝<x₁₁,x₁₂,...,x_1m>和X₂＝<x₂₁,x₂₂,...,x_2m>，其欧氏距离的计算公式如下：

其中x_1i和x_2i分别表示特征向量X₁和X₂的第i个特征指标值，m为特征向量中特征变量的总数；

然后按照距离递增顺序对样本集中各砂岩薄片图像排序，选取前s个相似的带标注砂岩薄片图像，构建训练集；如果有多张目标砂岩薄片图像，则分别针对每一张薄片图像，从样本集中选取s个近邻薄片图像，取其合集构成训练集。

上述步骤4)中采用重采样技术完成训练集类间平衡的处理过程是：首先统计训练集中各类砂岩薄片图像的数量；选择数量最多的砂岩种类，以其数量为基准值；如果其他种类的薄片图像数量少于基准值的80％，则采用随机重采样技术，创建该种类薄片图像的复本，直至数量达到基准值的80％。

上述步骤4)中计算每个特征指标的ANOVAF值，排除不必要特征指标的处理过程是：首先基于训练集中每一个特征指标，计算其ANOVAF值，公式如下：

其中i表示特征指标的序号，m表示特征指标的数量，即薄片图像特征向量的长度；n表示训练集中薄片图像实例的数量；μ_i表示特征指标i的均值；μ表示训练集中所有指标值的均值；x_ij表示特征指标i在第j个图像实例中的指标值；

然后设定百分比，保留30％的特征指标，计算要保留的特征指标数量M；

随后删除ANOVAF值最小的特征指标，如果剩余的特征指标数量大于M，则重新计算特征指标的ANOVAF值，重复删除过程；最后得到精化后的训练集。

上述步骤5)中采用随机森林分类器自动预测砂岩薄片种类的计算过程是：首先在训练集上训练随机森林分类器；然后基于目标地域砂岩显微薄片图像的特征向量，删除其中不必要的特征指标，即训练集精化过程中删除的特征指标集，作为输入由随机森林分类器预测砂岩薄片的种类。

本发明方法基于砂岩显微薄片图像，运用图像处理技术和机器学习方法，充分利用砂岩纹理清晰，颗粒形状独特的特性，提取纹理和形状特征指标，适用于砂岩的自动分类；特别针对由于地域不同导致的砂岩薄片结构差异较大而难于互相借鉴的问题，本发明提出基于欧氏距离采用近邻法选择与目标砂岩薄片相似的薄片图像样本构建训练集，并对训练集数据进行精化，训练随机森林分类器，从而达到理想的分类效果。本发明方法计算简单，具有扩展性，分类效率高，可有效降低砂岩分类时间和成本，可随着砂岩薄片数据储备的增加而提高自动分类的准确性；在油气勘探、砂岩分类和矿物鉴别中具有应用价值。

附图说明

图1是不同地域砂岩显微薄片自动分类方法的总体框架；

图2是岩石学领域砂岩分类的三角图；

图3是来自不同地域的石英砂岩和长石砂岩显微薄片示意图；

图4是基于砂岩显微薄片构建并精化训练集的处理流程图。

具体实施方式

本发明的主要目的是实现对不同地域砂岩显微薄片的自动分类，运用图像处理技术和机器学习方法，以矿物颗粒为单位统计纹理及形状特征，构成薄片图像的特征向量，基于欧氏距离采用近邻法选择与目标砂岩薄片样本相似的薄片样本构建训练集，并对训练集数据进行精化，训练随机森林分类器，自动预测目标地域砂岩薄片的种类。

图1所示为不同地域砂岩显微薄片自动分类方法的技术框架。方法的输入是砂岩显微薄片图像；方法的输出是砂岩薄片的种类(石英砂岩、长石砂岩或者岩屑砂岩)。为保证方法的正确应用，需要预先制备砂岩显微薄片图像的样本数据集，并在数据集中标注每张薄片图像属于哪一类砂岩。技术框架分为5个步骤：首先对砂岩的显微薄片图像进行预处理，分割矿物颗粒；再以颗粒为单位分别提取纹理和形状特征，构成薄片图像的特征向量；对目标地域的砂岩薄片，基于欧氏距离采用近邻法构建训练集；然后应用类间平衡、特征降维方法精化训练集数据；最后训练随机森林分类器，预测目标砂岩薄片的种类。

图2所示为岩石学领域砂岩分类的三角图。砂岩是一种陆源沉积岩，主要由砂粒胶结而成。绝大部分砂岩由石英和长石组成。砂岩按照其成分和成因可划分为：石英砂岩、长石砂岩和岩屑砂岩三大类。石英砂岩中，石英和各种硅质岩屑的含量占砂级岩屑总量的90％以上，仅含少量或不含长石、岩屑和重矿物；长石砂岩中，碎屑成分主要是石英和长石，其中长石碎屑含量占砂级岩屑总量的25％以上，石英含量低于75％，可能包含较多的云母和重矿物；岩屑砂岩中，岩屑含量占砂级碎屑总量的25％以上，石英含量低于75％，并可含少量长石(低于25％)，重矿物含量较高且种类复杂。砂岩是一种使用最广泛的建筑用石材，还是石油和天然气的储集层，砂岩的分类在油气和矿物勘探等方面有重要意义。

图3所示为来自不同地域的石英砂岩和长石砂岩显微薄片的对比图。从图中可以看出，不同地域不同种类砂岩由于其组成成分相同，颜色相似，不容易从颜色上进行区分；但是，砂岩通常具有较为清晰和丰富的纹理和颗粒形状，适于通过纹理和形状特征予以分类；由于形成的地质地域或环境不同，不同地域砂岩的结构差异较大，相互之间难于借鉴；例如图3中山南地区的长石砂岩与日喀则地区长石砂岩有明显的不同，反而在结构上与日喀则地区石英砂岩类似，而日喀则地区石英砂岩与那曲地区石英砂岩也完全不同；给砂岩显微薄片的自动分类带来了很大困难。

因此，本发明提出基于欧氏距离采用近邻法选择与目标砂岩薄片样本相似的薄片样本构建训练集，并对训练集数据进行精化，训练随机森林分类器，从而达到更好的分类效果。本发明采用的步骤如下：

1)读取砂岩显微薄片图像，对其进行滤波、去噪、平滑等预处理；采用图像分割技术，对砂岩显微薄片图像进行区域分割，获得矿物颗粒图像；

2)基于砂岩显微薄片图像中的矿物颗粒，分别计算纹理和形状特征指标，构建薄片图像的特征向量；

3)对目标地域的砂岩薄片图像，计算其特征向量，采用近邻法选择若干个相似的有标注的砂岩薄片图像样本，构建训练集合；

4)精化训练集数据：采用随机取样方法对初始训练集中少数类的样本进行重采样，实现训练集数据的类间平衡；分别计算每个特征指标的ANOVAF值，设定百分比，排除不必要的特征指标；

5)训练随机森林分类器，以目标地域砂岩薄片图像为输入，自动预测其种类。

上述步骤2)中基于砂岩显微薄片图像的矿物颗粒构建薄片图像的特征向量的处理过程是：首先针对薄片图像中的每个矿物颗粒，分别计算纹理特征和形状特征指标，组成矿物颗粒的特征向量X；然后对所有矿物颗粒的特征向量加权平均，得到均值特征向量X_m，计算公式如下：

其中k为薄片图像中包含的矿物颗粒总数；再后选取前K个面积最大的矿物颗粒的特征向量X_i(i∈[1,K])；最后将X₁,X₂,...,X_K,X_m串接，构成薄片图像的特征向量。

上述步骤2)中计算矿物颗粒的纹理特征的处理过程是：首先将颗粒图像转换为灰度图像，计算位移矢量d分别为(1,0)、(1,1)、(0,1)以及(-1,1)的4个归一化灰度共生矩阵P：矩阵P是一个256×256的矩阵，其第i行第j列的矩阵元素p[i,j](0≤i≤255,0≤j≤255)代表按矢量d分离的像素对的个数除以颗粒图像的总像素数量，其中第一个像素的灰度值为i，另一个像素的灰度值为j；分别按照4个位移矢量计算，可以得到4个灰度共生矩阵；

然后针对每一个灰度共生矩阵P，计算能量Energy、熵Entropy、对比度Contrast以及相关系数Correlation指标，计算公式如下，其中p[i,j](0≤i≤255,0≤j≤255)代表矩阵P中元素：

基于4个灰度共生矩阵，每个矿物颗粒获得16个纹理特征值。

上述步骤2)中计算矿物颗粒的形状特征的处理过程是：首先计算颗粒图像中包含的像素点总数，即颗粒的区域面积A；统计其边界像素点的个数，即区域周长C；然后求颗粒区域的凸包，指包含颗粒的最小凸多边形；计算凸包面积H，即凸包所包含的像素点总数；据以计算凸包差CD和密制度Solidity指标，计算公式如下：

CD＝H-A (10)

Solidity＝A/H (11)

接下来求包含颗粒区域的最小外接矩形，计算其面积ERA，即外接矩形包含的像素点数量；从而获得每个矿物颗粒的5个形状特征值：区域面积A、周长C、凸包差CD、密致度Solidity、以及外接矩形面积ERA；合并形状特征和纹理特征，每个矿物颗粒转换为包含21个特征值的特征向量X。

上述步骤3)中基于目标地域砂岩薄片图像的特征向量，采用近邻法构建训练集的计算过程是：按图4所示，对目标地域的砂岩薄片图像，首先计算其特征向量与现有样本集中各砂岩薄片图像特征向量的欧氏距离，给定两个特征向量X₁＝<x₁₁,x₁₂,...,x_1m>和X₂＝<x₂₁,x₂₂,...,x_2m>，其欧氏距离的计算公式如下：

其中x_1i和x_2i分别表示特征向量X₁和X₂的第i个特征指标值，m为特征向量中特征变量的总数；

然后按照距离递增顺序对样本集中各砂岩薄片图像排序，选取前s个相似的带标注砂岩薄片图像，构建训练集；如果有多张目标砂岩薄片图像，则分别针对每一张薄片图像，从样本集中选取s个近邻薄片图像，取其合集构成训练集。

上述步骤4)中采用重采样技术完成训练集类间平衡的处理过程是：首先统计训练集中各类砂岩薄片图像的数量；选择数量最多的砂岩种类，以其数量为基准值；如果其他种类的薄片图像数量少于基准值的80％，则采用随机重采样技术，创建该种类薄片图像的复本，直至数量达到基准值的80％。

上述步骤4)中计算每个特征指标的ANOVAF值，排除不必要特征指标的处理过程是：如图4所示，首先基于训练集中每一个特征指标，计算其ANOVAF值，公式如下：

其中i表示特征指标的序号，m表示特征指标的数量，即薄片图像特征向量的长度；n表示训练集中薄片图像实例的数量；μ_i表示特征指标i的均值；μ表示训练集中所有指标值的均值；x_ij表示特征指标i在第j个图像实例中的指标值；

然后设定百分比，保留30％的特征指标，计算要保留的特征指标数量M；

随后删除ANOVAF值最小的特征指标，如果剩余的特征指标数量大于M，则重新计算特征指标的ANOVAF值，重复删除过程；最后得到精化后的训练集。

上述步骤5)中采用随机森林分类器自动预测砂岩薄片种类的计算过程是：首先在训练集上训练随机森林分类器；然后基于目标地域砂岩显微薄片图像的特征向量，删除其中不必要的特征指标，即训练集精化过程中删除的特征指标集，作为输入由随机森林分类器预测砂岩薄片的种类。

本发明方法充分利用砂岩纹理清晰，颗粒形状独特的特性，提取纹理和形状特征指标，适用于砂岩的自动分类；特别针对由于地域不同导致的砂岩薄片结构差异较大而难于互相借鉴的问题，本发明提出基于欧氏距离采用近邻法选择与目标砂岩薄片相似的薄片图像样本构建训练集，并对训练集数据进行精化，训练随机森林分类器，从而达到理想的分类效果。通过收集西藏地区不同地域的砂岩薄片图像作为实验数据，实验结果表明本发明方法对不同地域砂岩薄片分类具有较高的准确性：精度达到89.3％，能够达到地质勘探的基本要求。另外，本发明方法具有较好的扩展性：其一可以采用其他的机器学习技术；其二通过砂岩薄片数据的不断储备，可以进一步提高砂岩分类的准确性。

本发明方法的具体应用途径很多，以上所述仅是本发明的优选实施方式。应当指出，对于本技术领域的普通技术人员来说，在不脱离本发明原理的前提下，还可以做出若干改进，这些改进也应视为本发明的保护范围。

Claims (6)

1.一种砂岩显微薄片的自动分类方法，其特征在于包括以下步骤：

1)读取砂岩显微薄片图像，对其进行预处理；之后采用图像分割技术对砂岩显微薄片图像进行区域分割，获得矿物颗粒图像；

2)基于砂岩显微薄片图像中的矿物颗粒，分别计算纹理和形状特征指标，构建薄片图像的特征向量；

3)对目标地域的砂岩薄片图像，计算其特征向量，采用近邻法选择若干个相似的有标注的砂岩薄片图像样本，构建训练集合；

4)精化训练集数据：采用随机取样方法对初始训练集中少数类的样本进行重采样，实现训练集数据的类间平衡；分别计算每个特征指标的ANOVA F值，设定百分比，排除不必要的特征指标；

5)训练随机森林分类器，以目标地域砂岩薄片图像为输入，自动预测其种类；

上述步骤2)基于砂岩显微薄片图像的矿物颗粒构建薄片图像的特征向量的处理过程是：首先针对薄片图像中的每个矿物颗粒，分别计算纹理特征和形状特征指标，组成矿物颗粒的特征向量X；然后对所有矿物颗粒的特征向量加权平均，得到均值特征向量X_m，计算公式如下：

其中k为薄片图像中包含的矿物颗粒总数；再后选取前K个面积最大的矿物颗粒的特征向量X_i(i∈[1,K])；最后将X₁,X₂,...,X_K,X_m串接，构成薄片图像的特征向量；

上述步骤2)中计算矿物颗粒的纹理特征指标的处理过程是：首先将颗粒图像转换为灰度图像，计算位移矢量d分别为(1,0)、(1,1)、(0,1)以及(-1,1)的4个归一化灰度共生矩阵P：矩阵P是一个256×256的矩阵，其第i行第j列的矩阵元素p[i,j](0≤i≤255,0≤j≤255)代表按矢量d分离的像素对的个数除以颗粒图像的总像素数量，像素对中第一个像素的灰度值为i，另一个像素的灰度值为j；分别按照4个位移矢量计算，共得到4个灰度共生矩阵；

然后针对每一个灰度共生矩阵P，计算能量Energy、熵Entropy、对比度Contrast以及相关系数Correlation指标，计算公式如下，其中p[i,j](0≤i≤255,0≤j≤255)代表矩阵P中元素：

基于4个灰度共生矩阵，每个矿物颗粒获得16个纹理特征值。

2.根据权利要求1所述的砂岩显微薄片的自动分类方法，其特征在于，

上述步骤2)中计算矿物颗粒的形状特征指标的处理过程是：首先计算颗粒图像中包含的像素点总数，即颗粒的区域面积A；统计其边界像素点的个数，即区域周长C；然后求颗粒区域的凸包，指包含颗粒的最小凸多边形；计算凸包面积H，即凸包所包含的像素点总数；据以计算凸包差CD和密制度Solidity指标，计算公式如下：

CD＝H-A

Solidity＝A/H

接下来求包含颗粒区域的最小外接矩形，计算其面积ERA，即外接矩形包含的像素点数量；从而获得每个矿物颗粒的5个形状特征值：区域面积A、周长C、凸包差CD、密致度Solidity、以及外接矩形面积ERA；合并形状特征和纹理特征，每个矿物颗粒转换为包含21个特征值的特征向量X。

3.根据权利要求2所述的砂岩显微薄片的自动分类方法，其特征在于，

上述步骤3)的过程是：对目标地域的砂岩薄片图像，首先计算其特征向量与现有样本集中各砂岩薄片图像特征向量的欧氏距离；给定两个特征向量X₁＝<x₁₁,x₁₂,...,x_1m>和X₂＝<x₂₁,x₂₂,...,x_2m>，其欧氏距离的计算公式如下：

其中x_1i和x_2i分别表示特征向量X₁和X₂的第i个特征指标值，m为特征向量中特征变量的总数；

然后按照距离递增顺序对样本集中各砂岩薄片图像排序，选取前s个相似的带标注砂岩薄片图像，构建训练集；如果有多张目标砂岩薄片图像，则分别针对每一张薄片图像，从样本集中选取s个近邻薄片图像，取其合集构成训练集。

4.根据权利要求3所述的砂岩显微薄片的自动分类方法，其特征在于，

上述步骤4)中采用重采样技术完成训练集类间平衡的处理过程是：首先统计训练集中各类砂岩薄片图像的数量；选择数量最多的砂岩种类，以其数量为基准值；如果其他种类的薄片图像数量少于基准值的80％，则采用随机重采样技术，创建该种类薄片图像的复本，直至数量达到基准值的80％。

5.根据权利要求4所述的砂岩显微薄片的自动分类方法，其特征在于，

上述步骤4)中计算每个特征指标的ANOVA F值，排除不必要特征指标的处理过程是：首先基于训练集中每一个特征指标，计算其ANOVA F值，公式如下：

其中i表示特征指标的序号，m表示特征指标的数量，即薄片图像特征向量的长度；n表示训练集中薄片图像实例的数量；μ_i表示特征指标i的均值；μ表示训练集中所有指标值的均值；x_ij表示特征指标i在第j个图像实例中的指标值；

然后设定百分比，保留30％的特征指标，计算要保留的特征指标数量M；

随后删除ANOVA F值最小的特征指标，如果剩余的特征指标数量大于M，则重新计算特征指标的ANOVA F值，重复删除过程；最后得到精化后的训练集。

6.根据权利要求5所述的砂岩显微薄片的自动分类方法，其特征在于，

上述步骤5)中采用随机森林分类器自动预测砂岩薄片种类的计算过程是：首先在训练集上训练随机森林分类器；然后基于目标地域砂岩显微薄片图像的特征向量，删除其中不必要的特征指标，即训练集精化过程中删除的特征指标集，作为输入由随机森林分类器预测砂岩薄片的种类。