一种砂岩薄片显微图像的自动分割方法与流程

本发明涉及砂岩薄片显微图像的自动分割问题，结合正交偏光图像和单偏光图像，运用图像处理技术和机器学习方法，通过对图像块的分类处理，实现将砂岩薄片显微图像中所包含的矿物颗粒逐个分割到单独区域。

背景技术：

砂岩在自然界的分布极广，是构成石油、天然气和地下水的主要储集层。砂岩薄片鉴定能够分析砂岩矿物组分及矿物含量，获取砂岩的类型和结构特征，在油气储集层探测和评估方面具有重要的应用价值。砂岩薄片鉴定的基础是对砂岩薄片显微图像进行图像分割。

现有的砂岩薄片显微图像分割以人工分割为主，耗时耗力，对个人能力和经验有较强依赖性，且其结果具有不可重复性。在砂岩薄片显微图像中，包含大量的矿物颗粒，各个颗粒之间的边界特征模糊，给图像分割带来了难度。此外，砂岩颗粒内部包含复杂多样的晶体微结构，常被误认为是颗粒边界，也对图像分割的结果产生影响。

技术实现要素：

本发明所要解决的技术问题是提供一种砂岩薄片显微图像的自动分割方法，该方法结合正交偏光图像和单偏光图像，运用图像处理技术和机器学习方法，实现砂岩薄片显微图像的自动分割。

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

1)读取砂岩正交偏光显微薄片图像，采用超像素分割技术将其预分割为图像块；

2)基于正交偏光显微薄片图像，提取图像块颜色特征和纹理特征；

3)读取砂岩单偏光显微薄片图像，采用边界检测技术提取图像边界；结合预分割图像块，提取图像块边界特征；

4)基于砂岩颗粒样本数据集训练支持向量机分类器；

5)使用训练过的分类器，计算每个图像块属于石英、长石、岩屑的概率；通过预设阈值确定图像块类别，高于阈值的图像块标注类型，其余图像块不标注类型；

6)根据图像块的颜色特征和纹理特征构造相似矩阵，使用标签传播算法为步骤5)未标注类型的图像块标注类型；

7)将类型相同、边界特征小于阈值的相邻图像块进行合并，形成完整的矿物颗粒。

上述步骤1)中将图像预分割为图像块的具体过程是：首先读入砂岩正交偏光显微图像，沿图像x轴和y轴两个方向以为间隔(其中n图像中像素数目)，将图像分成k个正方形区域，每个区域视为一个初始聚类；然后根据像素颜色特征和坐标值，迭代地计算其属聚类；迭代终止后，每个像素点构成的聚类即为一个图像块。

计算每个像素所属聚类的方法如下：根据k个初始聚类，计算每个聚类的中心{ck}k∈[1,k]，并以ck为中心搜索其周围2h×2h范围内的所有像素的集合ps；对ps内每个像素点i，计算其距离中心ck的距离，如果该距离小于上次迭代的距离则更新该距离，并设置该像素点所属聚类为k；对所有聚类中心迭代一次后，重新计算每个聚类中心以及新旧聚类中心的欧氏距离，如果距离小于预设阈值则停止算法，否则进入下一次迭代。

上述计算两个像素点i和ck之间距离的方法如下：

其中rc、gc和bc为ck的rgb颜色特征值，xc和yc为ck的x轴和y轴坐标；ri、gi和bi为像素点i的rgb颜色特征值，xi和yi为像素点i的x轴和y轴坐标，m优选为8。

上述步骤2)中计算图像块颜色特征的具体过程是：首先将图像块转换为灰度图像，计算基于直方图的颜色特征hist，以及基于统计的特征，包括：均值mean、方差variance、中值median、众数mode及平均绝对偏差mad五项统计指标。

基于直方图的颜色特征hist计算方法如下：将0～255灰度分为16个级别(每个级别包含16个灰度)，计算图像块内每个灰度级别所包含的像素数目{cnt1,cnt2,...,cnt16}，并将该向量归一化得到hist＝{hist1,hist2,...,hist16}；

基于统计的特征计算公式如下，其中n代表图像块中包含的像素数目，ip为像素p的灰度值：

median＝mid(i1,i2,…,in)，其中mid()是中值函数

通过上述计算，每个图像块的颜色特征为20维的特征向量。

上述步骤2)中计算图像块纹理特征的具体过程是：首先将图像块转换为灰度图像，并计算灰度共生矩阵；根据灰度共生矩阵计算纹理特征。

计算图像灰度共生矩阵的方法如下：在图像中任意取一点p1＝(x,y)及偏离它的一点p2＝(x+a,y+b)，计算点对(p1,p2)的灰度值(f1,f2)；让点p1在整幅图像上移动，得到不同的(f1,f2)值；统计每一种(f1,f2)值出现的次数；构建矩阵p，p(i,j)定义为灰度值(i,j)出现的次数；将矩阵p进行归一化，即将矩阵中每个元素除以所有灰度值出现的总次数。定义(a,b)分别为(0,1),(1,0),(1,1)以及(-1,1)，因此得到4个灰度共生矩阵。

根据灰度共生矩阵计算纹理特征的方法如下：分别计算能量energy、熵entropy、最大概率maxprobability、对比度contrast、倒数差分矩idm指标，公式如下：

每个灰度矩阵得到的纹理特征包含5个元素，因此每个图像块的纹理特征为20维的特征向量。

上述步骤3)中采用边界检测技术提取图像边界的具体过程是：首先读入砂岩单偏光显微薄片图像，将之转换为灰度图像，并对图像进行预处理；接着使用canny算子检测图像中的边界；最后，根据图像边界计算图像块的边界特征。

对图像进行预处理的方法如下：首先使用高斯滤波器对图像进行平滑处理；接着对图像进行伽马校正；最后进行直方图均衡化。

对图像进行伽马校正的方法是：对于像素点p，令其灰度值为ip，首先对ip归一化，公式为：x＝(ip+0.5)/256；接着对x进行预补偿，公式为：y＝x^1/γ，选择γ＝2.2；最后进行反归一化，公式为：z＝y*256-0.5，z即为伽马校正后像素点p的灰度值。

对图像进行直方图均衡化的方法是：令图像中有n个像素点，将图像灰度分为l个灰度级，l优选为16，计算灰度级为rk的像素的数目nk，则第k个灰度级出现的概率为p(rk)＝nk/n,k＝0,1,...,l-1，直方图均衡化公式如下：

其中sk是灰度级为k的像素经过变换以后的灰度。

上述步骤4)中训练支持向量机分类器的具体过程是：首先，收集石英、长石、岩屑颗粒的图像样本数据集；统计每类砂岩颗粒图像的数量，对数量较少的颗粒类型实施图像增强的方法，通过平移、翻转、旋转、缩放等方法生成额外的图像样本，使三种类型颗粒图像达到类间平衡；对于每一张砂岩颗粒样本图像，按照步骤2)的方法计算颜色特征和纹理特征，并构成特征向量；按照颗粒图像类别标注，构建分类器训练集；最后，基于图像训练集，训练支持向量机分类器。

上述步骤5)中对图像块类型进行标注的具体过程是：使用步骤4)训练好的支持向量机分类器对步骤2)计算得到的颜色特征向量和纹理特征向量进行预测，得到概率向量p＝{p1,p2,p3}，分别是该图像块属于石英、长石和岩屑的概率；令pmax＝max(p1,p2,p3)，当pmax≥0.7时，标注图像块的类型为pmax所对应的砂岩颗粒的种类。

上述步骤6)中对未标注图像块的类型进行预测的处理过程是：首先将步骤1)得到的每个图像块看作一个顶点，构建相似矩阵；接着使用标签传播算法对未标注类型的图像块进行类型标注。

构建相似矩阵的方法是：将步骤5)有标注类型的图像块构成已标注数据，记为(x1,y1),(x2,y2),…,(xl,yl)，其中xi表示第i个图像块的颜色特征向量和纹理特征向量构成的特征向量，yi∈{1,2,3}表示第i个图像块的类型；未标注类型的图像块记为(xl+1,yl+1),(xl+2,yl+2),…,(xk,yk)，这些图像块的类型未知；将每个图像块看作一个顶点，构建完全图，图中连接顶点i和j的边的权重定义为：

其中，α是超参，优选为1.5。根据wij构建相似矩阵w，并保留每行最大的8个值，其余值全部设置为0。

使用标签传播算法给图像块标记类型的方法是：首先根据相似矩阵w构建概率转移矩阵p，其公式为：

根据有标注数据(x1,y1),(x2,y2),…,(xl,yl)构建二维矩阵m1，其大小为l×3，该矩阵的第i行表示第i个数据样本的类型向量，即：如果第i个数据样本所属类型为j，则第j列元素值为1，其余列值为0；根据无标注数据(xl+1,yl+1),(xl+2,yl+2),…,(xk,yk)构建二维矩阵m2，其大小为(k-l)×3，该矩阵的每个位置元素初始值均设置为1/3；将m1和m2两个矩阵合并得到矩阵f＝{m1；m2}，其大小为k×3；最后进行标签传播，方法如下：

a)执行标签传播：f＝pf；

b)重置f中的有标注数据样本的标签：即将f中前l行重置为m1；

c)重复步骤a)和b)直至收敛。

f中后k-l行的构成的子矩阵即表示未标注图像块的类型。

上述步骤7)中合并图像块的具体过程是：首先计算所有的相邻图像块对d；对于图像块{bi,bj}∈d，根据步骤8)得到bi和bj的类型，如果类型不同则不处理，否则根据步骤3)计算得到bi和bj之间的边界特征，如果边界特征小于预设阈值则合并两个图像块，否则不做处理。

本发明方法基于砂岩正交偏光显微薄片图像和单偏光显微薄片图像，运用图像处理技术和机器学习方法，将图像预分割为图像块，在正交偏光图像中提取颜色特征和纹理特征，在单偏光图像中提取边界特征；训练支持向量机用于分类石英、长石和岩屑颗粒，对预分割图像块进行类型预测，并使用标签传播算法进一步处理类型模糊的图像块；最后根据图像块类型进行区域合并，实现砂岩薄片显微图像的自动分割。本发明方法分割效果好，扩展性强，计算简便，效率高，能够帮助地质学从业人员有效降低砂岩薄片显微图像分割时间并提高结果的准确率；在砂岩薄片鉴定、油气储集层探测和评估方面具有应用价值。

附图说明

图1是砂岩薄片显微图像自动分割方法的总体框架；

图2是砂岩正交偏光显微薄片图像和单偏光显微薄片图像；

图3是对砂岩薄片显微图像进行预分割的处理流程图。

具体实施方式

本发明的主要目的是实现对砂岩薄片显微图像的自动分割，运用图像处理技术和机器学习方法，以超像素算法将图像预分割为图像块，从正交偏光图像提取颜色和纹理特征，从单偏光图像提取边界特征；以石英、长石和岩屑三种颗粒图像训练支持向量机，对图像块类型进行预测，并用标签传播算法进行进一步处理；通过图像块的相邻关系、类型和边界特征进行图像块合并，实现砂岩薄片显微图像自动分割。

图1所示为砂岩薄片显微图像自动分类方法的技术框架。方法的输入是砂岩薄片显微图像(包括正交偏光图像和单偏光图像)；方法的输出是砂岩薄片显微图像的分割结果。为保证方法的正确应用，需要预先制备带标注的矿物颗粒图像，包括石英颗粒图像、长石颗粒图像和岩屑颗粒图像，作为训练样本集。技术框架分为7个步骤：采用超像素分割技术，将输入的正交偏光显微薄片图像预分割为图像块；提取图像块颜色特征和纹理特征；基于单偏光显微薄片图像，提取图像边界特征；基于矿物颗粒样本数据集训练支持向量机分类器；计算每个图像块属于石英、长石、岩屑的概率，通过预设条件确定图像块类型；使用标签传播算法为未标注类型的图像块标注类型；将类型相同、边界特征小于阈值的相邻图像块进行合并，形成完整的矿物颗粒。

图2所示为同一张砂岩薄片所制作而成的正交偏光显微图像和单偏光显微图像。正交偏光显微图像中，矿物颗粒呈现不同颜色的干涉色，部分相邻矿物颗粒由于颜色比较接近，因此之间的边界模糊，难以划分；此外，长石颗粒中存在较复杂的晶体微结构，如：双晶、解理和裂理等，岩屑颗粒中存在大量琐碎的矿物碎屑，这些微结构往往被误认为是颗粒边界而产生错误的分割。单偏光显微图像中，矿物颗粒大多为无色或透明的，能够从中提取部分矿物颗粒边界，并且长石颗粒和岩屑颗粒内部由于晶体微结构产生的伪边界现象不明显。因此，结合正交偏光图像和单偏光图像能够提高砂岩薄片显微图像的准确率。

因此，本发明提出基于超像素算法预分割砂岩图像，结合正交偏光图像和单偏光图像提取图像块特征，并训练分类器识别图像块类型，对未能识别图像块使用标签传播算法进一步处理，从而达到更好的分割效果。本发明采用的步骤如下：

1)读取砂岩正交偏光显微薄片图像，采用超像素分割技术将其预分割为图像块；

2)基于正交偏光显微薄片图像，提取图像块颜色特征和纹理特征；

3)读取砂岩单偏光显微薄片图像，采用边界检测技术提取图像边界；结合预分割图像块，提取图像块边界特征；

4)基于砂岩颗粒样本数据集训练支持向量机分类器；

5)使用训练过的分类器，计算每个图像块属于石英、长石、岩屑的概率；通过预设阈值确定图像块类别，高于阈值的图像块标注类型，其余图像块不标注类型；

6)根据图像块的颜色特征和纹理特征构造相似矩阵，使用标签传播算法为步骤5)未标注类型的图像块标注类型；

7)将类型相同、边界特征小于阈值的相邻图像块进行合并，形成完整的矿物颗粒。

上述步骤1)中将图像预分割为图像块的具体过程是：首先读入砂岩正交偏光显微图像，沿图像x轴和y轴两个方向以为间隔(其中n图像中像素数目)，将图像分成k个正方形区域，每个区域视为一个初始聚类；然后根据像素颜色特征和坐标值，迭代地计算其属聚类；迭代终止后，每个像素点构成的聚类即为一个图像块。

计算每个像素所属聚类的方法如下：根据k个初始聚类，计算每个聚类的中心{ck}k∈[1,k]，并以ck为中心搜索其周围2h×2h范围内的所有像素的集合ps；对ps内每个像素点i，计算其距离中心ck的距离，如果该距离小于上次迭代的距离则更新该距离，并设置该像素点所属聚类为k；对所有聚类中心迭代一次后，重新计算每个聚类中心以及新旧聚类中心的欧氏距离，如果距离小于预设阈值则停止算法，否则进入下一次迭代。

上述计算两个像素点i和ck之间距离的方法如下：

其中rc、gc和bc为ck的rgb颜色特征值，xc和yc为ck的x轴和y轴坐标；ri、gi和bi为像素点i的rgb颜色特征值，xi和yi为像素点i的x轴和y轴坐标，m优选为8。

上述步骤2)中计算图像块颜色特征的具体过程是：首先将图像块转换为灰度图像，计算基于直方图的颜色特征hist，以及基于统计的特征，包括：均值mean、方差variance、中值median、众数mode及平均绝对偏差mad五项统计指标。

基于直方图的颜色特征hist计算方法如下：将0～255灰度分为16个级别(每个级别包含16个灰度)，计算图像块内每个灰度级别所包含的像素数目{cnt1,cnt2,...,cnt16}，并将该向量归一化得到hist＝{hist1,hist2,...,hist16}；

基于统计的特征计算公式如下，其中n代表图像块中包含的像素数目，ip为像素p的灰度值：

median＝mid(i1,i2,…,in)，其中mid()是中值函数

通过上述计算，每个图像块的颜色特征为20维的特征向量。

上述步骤2)中计算图像块纹理特征的具体过程是：首先将图像块转换为灰度图像，并计算灰度共生矩阵；根据灰度共生矩阵计算纹理特征。

计算图像灰度共生矩阵的方法如下：在图像中任意取一点p1＝(x,y)及偏离它的一点p2＝(x+a,y+b)，计算点对(p1,p2)的灰度值(f1,f2)；让点p1在整幅图像上移动，得到不同的(f1,f2)值；统计每一种(f1,f2)值出现的次数；构建矩阵p，p(i,j)定义为灰度值(i,j)出现的次数；将矩阵p进行归一化，即将矩阵中每个元素除以所有灰度值出现的总次数。定义(a,b)分别为(0,1),(1,0),(1,1)以及(-1,1)，因此得到4个灰度共生矩阵。

根据灰度共生矩阵计算纹理特征的方法如下：分别计算能量energy、熵entropy、最大概率maxprobability、对比度contrast、倒数差分矩idm指标，公式如下：

每个灰度矩阵得到的纹理特征包含5个元素，因此每个图像块的纹理特征为20维的特征向量。

上述步骤3)中采用边界检测技术提取图像边界的具体过程是：首先读入砂岩单偏光显微薄片图像，将之转换为灰度图像，并对图像进行预处理；接着使用canny算子检测图像中的边界；最后，根据图像边界计算图像块的边界特征。

对图像进行预处理的方法如下：首先使用高斯滤波器对图像进行平滑处理；接着对图像进行伽马校正；最后进行直方图均衡化。

对图像进行伽马校正的方法是：对于像素点p，令其灰度值为ip，首先对ip归一化，公式为：x＝(ip+0.5)/256；接着对x进行预补偿，公式为：y＝x^1/γ，选择γ＝2.2；最后进行反归一化，公式为：z＝y*256-0.5，z即为伽马校正后像素点p的灰度值。

对图像进行直方图均衡化的方法是：令图像中有n个像素点，将图像灰度分为l个灰度级，l优选为16，计算灰度级为rk的像素的数目nk，则第k个灰度级出现的概率为p(rk)＝nk/n,k＝0,1,...,l-1，直方图均衡化公式如下：

其中sk是灰度级为k的像素经过变换以后的灰度。

上述步骤4)中训练支持向量机分类器的具体过程是：首先，收集石英、长石、岩屑颗粒的图像样本数据集；统计每类砂岩颗粒图像的数量，对数量较少的颗粒类型实施图像增强的方法，通过平移、翻转、旋转、缩放等方法生成额外的图像样本，使三种类型颗粒图像达到类间平衡；对于每一张砂岩颗粒样本图像，按照步骤2)的方法计算颜色特征和纹理特征，并构成特征向量；按照颗粒图像类别标注，构建分类器训练集；最后，基于图像训练集，训练支持向量机分类器。

上述步骤5)中对图像块类型进行标注的具体过程是：使用步骤4)训练好的支持向量机分类器对步骤2)计算得到的颜色特征向量和纹理特征向量进行预测，得到概率向量p＝{p1,p2,p3}，分别是该图像块属于石英、长石和岩屑的概率；令pmax＝max(p1,p2,p3)，当pmax≥0.7时，标注图像块的类型为pmax所对应的砂岩颗粒的种类。

上述步骤6)中对未标注图像块的类型进行预测的处理过程是：首先将步骤1)得到的每个图像块看作一个顶点，构建相似矩阵；接着使用标签传播算法对未标注类型的图像块进行类型标注。

构建相似矩阵的方法是：将步骤5)有标注类型的图像块构成已标注数据，记为(x1,y1),(x2,y2),…,(xl,yl)，其中xi表示第i个图像块的颜色特征向量和纹理特征向量构成的特征向量，yi∈{1,2,3}表示第i个图像块的类型；未标注类型的图像块记为(xl+1,yl+1),(xl+2,yl+2),…,(xk,yk)，这些图像块的类型未知；将每个图像块看作一个顶点，构建完全图，图中连接顶点i和j的边的权重定义为：

其中，α是超参，优选为1.5。根据wij构建相似矩阵w，并保留每行最大的8个值，其余值全部设置为0。

使用标签传播算法给图像块标记类型的方法是：首先根据相似矩阵w构建概率转移矩阵p，其公式为：

根据有标注数据(x1,y1),(x2,y2),…,(xl,yl)构建二维矩阵m1，其大小为l×3，该矩阵的第i行表示第i个数据样本的类型向量，即：如果第i个数据样本所属类型为j，则第j列元素值为1，其余列值为0；根据无标注数据(xl+1,yl+1),(xl+2,yl+2),…,(xk,yk)构建二维矩阵m2，其大小为(k-l)×3，该矩阵的每个位置元素初始值均设置为1/3；将m1和m2两个矩阵合并得到矩阵f＝{m1；m2}，其大小为k×3；最后进行标签传播，方法如下：

a)执行标签传播：f＝pf；

b)重置f中的有标注数据样本的标签：即将f中前l行重置为m1；

c)重复步骤a)和b)直至收敛。

f中后k-l行的构成的子矩阵即表示未标注图像块的类型。

上述步骤7)中合并图像块的具体过程是：首先计算所有的相邻图像块对d；对于图像块{bi,bj}∈d，根据步骤8)得到bi和bj的类型，如果类型不同则不处理，否则根据步骤3)计算得到bi和bj之间的边界特征，如果边界特征小于预设阈值则合并两个图像块，否则不做处理。

本发明方法充分利用砂岩正交偏光图像和单偏光图像的特点，分别提取其颜色特征、纹理特征及边界特征，应用图像处理技术、机器学习方法，实现砂岩薄片显微图像的自动分割；针对砂岩薄片图像中相邻颗粒特征接近、边界模糊的问题，本发明结合单偏光图像辅助边界检测；针对长石和岩屑颗粒中晶体微结构较复杂，影响分割的问题，本发明通过提取正交偏光图像颜色和纹理特征对图像块进行分类，并结合标签传播算法进一步确定图像块所属矿物颗粒类型，避免将晶体微结构误分为边界的发生，从而达到理想的分割效果。通过从西藏等地采集和制作的砂岩薄片显微图像作为实验数据，实验结果表明本发明方法对各种不同类型砂岩薄片的分割具有较高的准确率，能够达到薄片鉴定的基本要求。另外，本发明方法可以应用到其他类型的沉积岩，如：页岩、石灰岩等薄片显微图像的分割，具有较好的扩展性。

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