一种基于EM与孔隙度的数字岩心孔隙分割方法

一种基于em与孔隙度的数字岩心孔隙分割方法
技术领域
1.本发明属于数字岩心研究领域，具体涉及一种基于em与孔隙度的数字岩心孔隙分割方法。


背景技术：

2.随着石油工业的不断发展，对油气资源的需求也在不断增加，目前对油气层的开采已经逐渐由常规储层向非常规储层转化。因此，岩石物理研究在油气藏评价中所占据的地位也越来越重要。岩石物理实验分析物性时在定量表征和微观参数分析上具有局限性，同时，这些实验开展起来耗时长，难度大，并且非常规储层结构更加复杂会导致取心成功率低。
3.数字岩心技术是基于ct技术和信息处理技术发展起来的，近几年来逐渐应用到岩石物理中，比较热门的研究方向有采用数值模拟进行孔隙度、粒度分布、渗透率和饱和度等的计算和分析。数字岩心技术能够在无损情况下提取岩心的内部孔隙结构，并提供可视化的三维结果，对储层物性的模拟和分析具有重要帮助，对油藏勘探开发也具有重要意义。相比于传统的岩石物理实验，数字岩心技术具有速度快、效率高、可重复性强的优点，相比于薄片分析，数字岩心技术能够重建岩心内的三维结构，能够更加深入地分析岩石的微观特征，进而有利于建立储层参数和岩石物理特性的关系。因此，ct技术越来越成为储层物性分析的有利手段，对数字岩心的研究具有重要的意义。
4.数字岩心由借助于ct扫描技术采集的图像序列构成，一般为灰度图像，其中孔隙对应图像中较暗即灰度值较低的部分，单张数字岩心中内容简单，孔隙特征比较明显，但是孔隙结构复杂无固定形态和特征、孔径范围跨域微米到纳米、不同ct数据中其对比度与饱和度也不相同。数字岩心孔隙的提取与重建效果依赖于孔隙的分割技术，也影响基于孔隙的后续模拟环节。因此，精准的孔隙分割是数字岩心研究的首要问题。
5.传统的数字岩心分割方法主要基于像素的亮度、形态学特征和纹理特征等进行，主要方法有灰度阈值法、多阈值法、区域生长法等。并且，学者们在此基础上也发展出了更先进的处理方法，比如，zhipeng xu等人针对孔隙ct扫描三维图像提出了一种新的测定方法，对图像孔隙结构的几何和拓扑参数进行计算，并利用常用的数学形态学运算进行孔隙检测。deng等人通过裂缝扩展过程生成定制掩模，在强度直方图上放大信号产生用于局部阈值提取的裂缝，邱茂鑫提出采用图像灰度以及协方差建立两个阈值，对ct图像中灰度值进行分割以提取孔隙。
6.随着信息处理技术和计算机技术的发展，学者们也开始尝试用基于机器学习的方法实现岩石孔隙的分割。比如，程国建等对薄片采用均值聚类与bp神经网络相结合的方法提取岩石图像中的孔隙。刘春等改进了种子算法，也是从铸体薄片图像上自动分割出孔隙，并计算孔隙的几何特征。姜枫对铸体薄片完善了整个基于机器学习的孔隙识别过程，使用迁移学习和多角度图像扩增进行网络参数学习，使用rock net网络自动提取特征。
7.上述这些方法对数字岩心的孔隙分割、识别、检测和分类起到了有效的帮助，但也
存在不同程度的弊端。主要体现在，传统方法简单易实现，也是目前工程应用中的常用方法，但该种方法需要与人交互完成，分割精度也较低，且受ct扫描效果的影响，基于孔隙分割后计算出的孔隙度往往与实验偏差较大。而基于机器学习的方法能够自动完成分割过程，但其单张ct图像的计算量较大，目前也主要针对铸体薄片实现，对ct图像进行孔隙提取的较少。同时，不同岩心、不同扫描过程获得的数字岩心灰度效果不同，目前尚无统一的训练样本库，方法的泛化能力也不高。
8.因此，如何通过分析数字岩心自身的特点，设计出一种具有较强通用性，且分割精度与分割效率较高的孔隙分割方法，是本领域内一个亟待解决的问题。


技术实现要素：

9.本发明实施例的目的在于提供一种基于em与孔隙度的数字岩心孔隙分割方法，以实现提高孔隙分割方法的通用性，以及提高分割精度与分割效率的目的。具体技术方案如下：
10.从待孔隙分割的数字岩心数据体中选择满足预设要求的预设数量个样本；其中，所述预设要求为每个样本的均值和方差分别对应小于所述数字岩心数据体的均值和方差；
11.获取em算法在当前次迭代中初始的概率参数；其中，所述概率参数包括孔隙和背景的先验概率，以及孔隙和背景的概率分布参数；
12.基于所述预设数量个样本和所述当前次迭代中初始的概率参数，得到当前次迭代中更新的概率参数；其中，所述当前次迭代中更新的概率参数包括当前次迭代中更新的先验概率和更新的概率分布参数；
13.判断所述当前次迭代是否满足迭代停止条件；
14.若否，利用所述当前次迭代中更新的概率分布参数，计算当前次迭代的孔隙度；基于所述当前次迭代的孔隙度、上一次迭代的孔隙度，以及对所述数字岩心数据体对应的岩心样本进行实测得到的实测孔隙度，确定迭代速度的待变化类别；根据所述待变化类别对所述当前次迭代中更新的概率分布参数进行修正，以得到下一次迭代中初始的概率参数继续进行迭代；其中，所述待变化类别包括加速和减速；
15.若是，利用所述当前次迭代中更新的概率分布参数，获得所述数字岩心数据体中各像素的孔隙分割结果。
16.在本发明的一个实施例中，所述从待孔隙分割的数字岩心数据体中选择满足预设要求的预设数量个样本，包括：
17.从所述数字岩心数据体中随机选取多个数字岩心切片图像；
18.对所述多个数字岩心切片图像中的每一个，分别计算均值和方差；
19.将获得的所有均值求取平均值的结果作为所述数字岩心数据体的均值；并将获得的所有方差求取平均值的结果作为所述数字岩心数据体的方差；
20.每次从所述数字岩心数据体中选取一个未使用过的数字岩心切片图像，并利用预设尺寸从中选取一部分子图像作为候选样本图像，计算该候选样本图像的均值和方差，在确定该候选样本图像的均值和方差满足所述预设要求时将该候选样本图像确定为样本，并重复该样本选取过程直至样本数量达到所述预设数量。
21.在本发明的一个实施例中，所述孔隙和背景的概率分布参数，包括：
22.孔隙类在高斯概率分布下的均值、方差，以及背景类在高斯概率分布下的均值、方差。
23.在本发明的一个实施例中，所述基于所述预设数量个样本和所述当前次迭代中初始的概率参数，得到当前次迭代中更新的概率参数，包括：
24.根据所述预设数量个样本中各像素的灰度值、所述当前次迭代中初始的概率参数，计算表征后验概率的各灰度值对应的孔隙和背景的分类软标签；
25.利用所述各灰度值对应的孔隙和背景的分类软标签，重新计算孔隙和背景的先验概率，得到当前次迭代中更新的先验概率；
26.利用所述当前次迭代中更新的先验概率，重新计算孔隙和背景的概率分布参数，得到当前次迭代中更新的概率分布参数。
27.在本发明的一个实施例中，所述利用所述当前次迭代中更新的概率分布参数，计算当前次迭代的孔隙度，包括：
28.根据所述当前次迭代中更新的概率分布参数，计算所述预设数量个样本中每个像素属于孔隙和背景的概率；
29.利用所述预设数量个样本中每个像素属于孔隙和背景的概率，对所述数字岩心数据体进行孔隙分割，得到孔隙分割的二值结果；
30.根据所述孔隙分割的二值结果，计算出所述数字岩心数据体在当前次迭代的孔隙度。
31.在本发明的一个实施例中，所述基于所述当前次迭代的孔隙度、上一次迭代的孔隙度，以及对所述数字岩心数据体对应的岩心样本进行实测得到的实测孔隙度，确定迭代速度的待变化类别，包括：
32.利用所述当前次迭代的孔隙度、所述上一次迭代的孔隙度，以及所述实测孔隙度，计算得到当前次迭代的待变化类别初始判别值，并对所述待变化类别初始判别值进行归一化处理，得到当前次迭代的待变化类别判别值；
33.若所述待变化类别判别值大于0.5，则确定迭代速度的待变化类别为加速；反之，则确定迭代速度的待变化类别为减速。
34.在本发明的一个实施例中，所述当前次迭代的待变化类别初始判别值的计算公式，包括：
[0035][0036]
其中，δpsi表示所述当前次迭代的待变化类别初始判别值；psi表示所述当前次迭代的孔隙度；ps
i-1
表示所述上一次迭代的孔隙度；ps
test
表示所述实测孔隙度；
[0037]
所述对所述待变化类别初始判别值进行归一化处理，得到当前次迭代的待变化类别判别值的过程中，所采用的计算公式包括：
[0038][0039]
其中，δps'i表示当前次迭代的待变化类别判别值。
[0040]
在本发明的一个实施例中，所述根据所述待变化类别对所述当前次迭代中更新的
概率分布参数进行修正，包括：
[0041]
根据所述当前次迭代的待变化类别判别值计算对应的修正权值；
[0042]
基于所述修正权值，利用待变化类别对应的调节公式，对所述当前次迭代中更新的概率分布参数中的孔隙类和背景类在高斯概率分布下的均值进行修正。
[0043]
在本发明的一个实施例中，所述根据所述当前次迭代的待变化类别判别值计算对应的修正权值的过程中，所采用的计算公式，包括：
[0044][0045]
其中，ω表示所述修正权值。
[0046]
在本发明的一个实施例中，所述待变化类别对应的调节公式，包括：
[0047]mcj
'＝(1 ω)
×mcj
[0048]
其中，m
cj
表示所述当前次迭代中更新的概率分布参数中的孔隙类和背景类在高斯概率分布下的均值；j＝1,2分别表示孔隙类和背景类；m
cj
'表示对m
cj
进行修正后的数值。
[0049]
本发明实施例所提供的基于em与孔隙度的数字岩心孔隙分割方法中，首先，为避免高灰度值矿物的影响，基于灰度分布特征从待孔隙分割的数字岩心数据体中选择满足预设要求的预设数量个样本，使得样本的均值和方差分别对应小于所述数字岩心数据体的均值和方差。同时，该种样本筛选方式使得em迭代过程所需要的样本从待分割的数岩心数据体自身采集，避免了样本库数据不匹配的问题，能够提高算法的通用性。然后，对预设数量个样本采用em算法进行无监督的学习，针对每次迭代，利用初始的概率参数得到当前次迭代中更新的概率参数，在满足迭代停止条件时，利用当前次迭代中更新的概率分布参数，获得所述数字岩心数据体中各像素的孔隙分割结果。在不满足迭代停止条件时，利用所述当前次迭代中更新的概率分布参数，计算当前次迭代的孔隙度，并结合实验获得的实测孔隙度确定迭代加速还是减速，并相应对所述当前次迭代中更新的概率分布参数进行修正，以得到下一次迭代中初始的概率参数继续进行迭代。本发明实施例通过在参数学习过程中结合实测孔隙度精细调节概率分布参数的变化，以确定最终的概率分布参数，能够避免参数调节过度，并能够加速迭代进程，降低运算耗时。因此在分割数字岩心的孔隙时能够提高孔隙自动分割的精度和收敛的速度。
附图说明
[0050]
图1为本发明实施例所提供的一种基于em与孔隙度的数字岩心孔隙分割方法的流程示意图；
[0051]
图2为本发明实施例的修正权值的数学波形图；
[0052]
图3为本发明实施例实验中对相同的待分割的数字岩心切片图像用不同样本学习分割时的影响对比结果图；其中，图3(a)为待分割的数字岩心切片图像；图3(b)从左至右依次为样本1、样本2、样本3和样本4对应的图像；图3(c)～图3(g)分别是样本1、样本2、样本3和样本4经过经典em算法迭代1次、6次、11次、16次和21次的分割结果；
[0053]
图4为本发明实施例实验给出的样本筛选示例；其中，图4(a)为本发明实施例方法筛选出的满足预设要求的3个典型学习样本；图4(b)为本发明实施例方法筛选出的3个不合格样本；
[0054]
图5为本发明实施例实验针对相同样本，采用otus方法、经典em算法和本发明实施例方法进行孔隙提取的效果图；其中，图5(a)为样本a-1、a-2和a-3对应的图像。图5(b)为样本a-1、a-2和a-3采用otus方法进行孔隙提取的效果图；图5(c)为样本a-1、a-2和a-3采用经典em算法进行孔隙提取的效果图；图5(d)为样本a-1、a-2和a-3采用本发明实施例方法进行孔隙提取的效果图。
具体实施方式
[0055]
下面将结合本发明实施例中的附图，对本发明实施例中的技术方案进行清楚、完整地描述，显然，所描述的实施例仅仅是本发明一部分实施例，而不是全部的实施例。基于本发明中的实施例，本领域普通技术人员在没有做出创造性劳动前提下所获得的所有其他实施例，都属于本发明保护的范围。
[0056]
如何分析数字岩心自身的特点设计出具有较强通用性，且能够综合分割精度与计算耗时的孔隙分割方法是本发明的研究目的，为此，本发明实施例提供了一种基于em与孔隙度的数字岩心孔隙分割方法。
[0057]
需要说明的是，本发明实施例所提供的一种基于em与孔隙度的数字岩心孔隙分割方法的执行主体可以为一种基于em与孔隙度的数字岩心孔隙分割装置，所述装置可以运行于电子设备中。其中，该电子设备可以为一服务器或终端设备，当然并不局限于此。
[0058]
如图1所示，本发明实施例所提供的一种基于em与孔隙度的数字岩心孔隙分割方法，可以包括如下步骤：
[0059]
s1，从待孔隙分割的数字岩心数据体中选择满足预设要求的预设数量个样本。
[0060]
该步骤是从待分割的数字岩心数据体中采集无标签的学习样本，作为学习和训练的依据。但可以理解的是，数字岩心数据体由多层数字岩心切片图像构成，而数字岩心切片图像中一般含有低灰度值的孔隙和其他灰度值的矿物这两种成分。由于矿物的种类复杂，在数字岩心切片图像中体现出的灰度值也不均匀，如果学习样本中混有较多高亮度值的矿物，会提高样本整体的灰度均值，因此会导致概率分布模型偏向高亮度值，容易使灰度值相对较低的矿物被划分成孔隙，这样就难以对孔隙进行有效的分割。因此在选择样本时应尽量避开混有显示较高灰度值的矿物。
[0061]
并且，发明人在研究过程中发现，当数字岩心切片图像中含有较多灰度值较高的矿物时，数字岩心数据体整体的均值和方差均较大，样本的质量会影响概率分布参数的选取，因此本发明实施例在随机选择样本的基础上基于样本的统计特征进行样本筛选，以保障样本的质量和学习的效果。
[0062]
基于上述原因，本发明实施例从待孔隙分割的数字岩心数据体中选择满足预设要求的预设数量个样本。其中，所述预设要求为每个样本的均值和方差分别对应小于所述数字岩心数据体的均值和方差。
[0063]
可选的一种实施方式中，s1可以包括以下步骤：
[0064]
s11，从所述数字岩心数据体中随机选取多个数字岩心切片图像。
[0065]
具体的，从所述数字岩心数据体的若干个数字岩心切片图像中，随机选取m张，其中，m为大于0的自然数。
[0066]
s12，对所述多个数字岩心切片图像中的每一个，分别计算均值和方差。
[0067]
具体的，对随机选出的m个数字岩心切片图像，分别计算各自的均值和方差。其中，第i个数字岩心切片图像的均值和方差分别以和表示。
[0068]
关于均值和方差的计算过程，请参见相关技术，在此不做具体说明。
[0069]
s13，将获得的所有均值求取平均值的结果作为所述数字岩心数据体的均值；并将获得的所有方差求取平均值的结果作为所述数字岩心数据体的方差。
[0070]
在该步骤中，将m个数字岩心切片图像获得的所有均值求和后再取平均值，作为所述数字岩心数据体的均值，以及将m个数字岩心切片图像获得的所有方差求和后再取平均值作为所述数字岩心数据体的方差，是以m个数字岩心切片图像的统计结果近似计算样本整体的灰度均值和方差。具体的，所述数字岩心数据体的均值和方差分别以和表示，计算方式如下：
[0071][0072]
s14，每次从所述数字岩心数据体中选取一个未使用过的数字岩心切片图像，并利用预设尺寸从中选取一部分子图像作为候选样本图像，计算该候选样本图像的均值和方差，在确定该候选样本图像的均值和方差满足所述预设要求时将该候选样本图像确定为样本，并重复该样本选取过程直至样本数量达到所述预设数量。
[0073]
其中，所述预设尺寸可以表示为n
×
m，n
×
m小于单张数字岩心切片图像的尺寸。因此，可以理解的是，从所述数字岩心数据体的同一个数字岩心切片图像中可以获取多个候选样本图像。
[0074]
针对该步骤，可以从所述数字岩心数据体的各层数字岩心切片图像中，每次按照一定的顺序选取未使用过的，也就是不重复的一个数字岩心切片图像，然后从该数字岩心切片图像中选取候选样本图像，或者每次也可以从所述数字岩心数据体的各层数字岩心切片图像随机选取不重复的一个数字岩心切片图像进而从中选取候选样本图像。
[0075]
针对选取出的每个候选样本图像，计算其均值和方差，如果该候选样本图像的均值小于所述数字岩心数据体的均值且该候选样本图像的方差小于所述数字岩心数据体的方差则确定该候选样本图像可以作为样本；反之，则继续选取下一个候选样本图像进行判断，直至选取出k个样本，其中k表示所述预设数量，可以根据需要设定一个大于0的具体数值。所选取出的k个样本作为无监督学习中样本在后续步骤使用。
[0076]
可选的一种实施方式中，每次可以从所述数字岩心数据体的各层数字岩心切片图像中，一次性随机选取出q个未使用过的候选样本图像，作为一批候选样本图像，然后针对该批候选样本图像中的每个分别计算均值和方差，判断其是否能够作为样本保留，如果该批候选样本图像所得到的样本数量不满足k个，则从所述数字岩心数据体未筛选到的位置再随机选取下一批的q个候选样本图像再次进行筛选，直至所得到的样本数量累计满足k个。其中q≥k。
[0077]
s2，获取em算法在当前次迭代中初始的概率参数。
[0078]
ct图像进行孔隙分割时，也是一种根据样本分布特点推理出某一个像素属于孔隙的概率是否大于属于背景的概率的问题，因此如果能够确定孔隙和背景的概率分布，也就能够实现分割。学者们在贝叶斯理论基础上，引入最大期望，建立了em算法(expectation-maximization algorithm，最大期望算法或称为期望最大化算法)，通过迭代计算出概率分布的参数以确定具体模型，但是该过程存在一定随机性，且迭代过程运算量较大。因此，本发明实施例以ct扫描的数字岩心特点出发，从效果和收敛性两方面提出新的数字岩心孔隙分割方法。
[0079]
本发明实施例中，所述概率参数包括孔隙和背景的先验概率，以及孔隙和背景的概率分布参数。
[0080]
基于em算法实现孔隙分割在初次进行迭代时，某个灰度值下的像素属于孔隙或者属于背景是未知的，因此需要初始化一个携带参数的孔隙-灰度的概率分布和一个背景-灰度的概率分布。可选的一种实施方式中，均可以采用高斯分布实现。
[0081]
为了便于描述，以c1和c2分别代表孔隙类和背景类。所述概率参数中孔隙和背景的先验概率可以分别用p
c1
和p
c2
表示。所述孔隙和背景的概率分布参数，包括：孔隙类在高斯概率分布下的均值m
c1
、方差s
c1
，以及背景类在高斯概率分布下的均值m
c2
、方差s
c2
。也就是说，针对k个样本，在em算法的每次迭代中，初始的概率参数包括6个数值：m
c1
，s
c1
，m
c2
，s
c2
，p
c1
，p
c2
。
[0082]
针对首次迭代，初始的概率参数可以根据需要设定预设值，比如可以为m
c1
＝0.05；s
c1
＝0.1；m
c2
＝0.5；s
c2
＝0.3；p
c1
＝0.5；p
c2
＝0.5。
[0083]
s3，基于所述预设数量个样本和所述当前次迭代中初始的概率参数，得到当前次迭代中更新的概率参数。
[0084]
其中，所述当前次迭代中更新的概率参数包括当前次迭代中更新的先验概率和更新的概率分布参数。
[0085]
可选的一种实施方式中，s3可以包括以下步骤：
[0086]
s31，根据所述预设数量个样本中各像素的灰度值、所述当前次迭代中初始的概率参数，计算表征后验概率的各灰度值对应的孔隙和背景的分类软标签。
[0087]
本发明实施例在贝叶斯理论基础上根据公式(2)和(3)获得样本中灰度值对应的孔隙和背景的分类软标签。
[0088][0089]
p
c2_soft
(i)＝1-p
c1_soft(i)ꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀ
(3)
[0090]
其中，i为像素点，x(i)表示像素灰度值，取值在[0-255]之间；pdf()表示高斯概率分布函数；p
c1_soft
(i)表示灰度值对应的孔隙的分类软标签；p
c2_soft
(i)表示灰度值对应的背景的分类软标签；p
c1_soft
(i)和p
c2_soft
(i)体现该灰度值属于孔隙或样本的概率，也就是后验概率。j＝1，2分别表示孔隙和背景。
[0091]
s32，利用所述各灰度值对应的孔隙和背景的分类软标签，重新计算孔隙和背景的先验概率，得到当前次迭代中更新的先验概率。
[0092]
该步骤所采用的计算公式如下：
[0093][0094]
其中，n表示样本中所含像素点个数，j＝1，2分别表示孔隙和背景。
[0095]
s33，利用所述当前次迭代中更新的先验概率，重新计算孔隙和背景的概率分布参数，得到当前次迭代中更新的概率分布参数。
[0096]
该步骤所采用的计算公式如下：
[0097][0098]
其中，公式(5)中等号左侧的m
cj
、s
cj
表示更新的概率分布参数；j＝1，2分别表示孔隙和背景；n
cj
表示样本中属于j类的像素点个数。
[0099]
s4，判断所述当前次迭代是否满足迭代停止条件。
[0100]
其中，迭代停止条件包括所述当前次迭代达到预设迭代次数，或者所述当前次迭代中更新的概率分布参数实现收敛。
[0101]
预设迭代次数可以根据需要设定，比如为30次等。所述当前次迭代中更新的概率分布参数实现收敛是指相比于上一次迭代中更新的概率分布参数，对应参数的差值小于一定的预设值，比如0.0001等。
[0102]
具体的，若满足上述任一个条件，则表示所述当前次迭代满足迭代停止条件，则执行s6；若两个条件都不满足，则表示所述当前次迭代不满足迭代停止条件，则执行s5。
[0103]
若否，执行s5，利用所述当前次迭代中更新的概率分布参数，计算当前次迭代的孔隙度；基于所述当前次迭代的孔隙度、上一次迭代的孔隙度，以及对所述数字岩心数据体对应的岩心样本进行实测得到的实测孔隙度，确定迭代速度的待变化类别；根据所述待变化类别对所述当前次迭代中更新的概率分布参数进行修正，以得到下一次迭代中初始的概率参数继续进行迭代。
[0104]
其中，所述待变化类别包括加速和减速。
[0105]
岩心的孔隙空间是油气赋集和运输的通道，孔隙度是反映岩心物性的重要参数，是孔隙体积与岩心体积的比值。常规方法是通过双室孔隙测试仪测量岩心样本获得，测量精度也较高，但这种方法仅仅能得到岩心的总孔隙度，也就是总孔隙体积与岩心体积的比值，无法获得岩心的微观孔隙形态特点，而微观孔隙形态会影响油气的运移等。通过对ct扫描得到的数字岩心切片图像进行数字岩心提取孔隙后能够建立可视化的孔隙结构，显示出岩心内所有的孔隙形态，并定量计算孔隙形态参数，如孔隙度等，但前提是孔隙提取的精度够高，也就是孔隙分割要准确，因此本发明实施例考虑将实测孔隙度与数字岩心的孔隙提取过程相结合，通过与实验测得孔隙度进行比较判断，以提高分割的准确性。
[0106]
可选的一种实施方式中，所述利用所述当前次迭代中更新的概率分布参数，计算当前次迭代的孔隙度，包括：
[0107]
a1，根据所述当前次迭代中更新的概率分布参数，计算所述预设数量个样本中每
个像素属于孔隙和背景的概率。
[0108]
具体的，根据s33公式(5)计算得到的更新的概率分布参数m
cj
和s
cj
，j＝1,2，可以计算k个样本中每个像素分别属于孔隙和背景的概率，比如一个像素属于孔隙和背景的概率分别为0.2、0.8等。该步骤的具体过程请结合相关技术理解，在此不做详细说明。
[0109]
a2，利用所述预设数量个样本中每个像素属于孔隙和背景的概率，对所述数字岩心数据体进行孔隙分割，得到孔隙分割的二值结果。
[0110]
针对k个样本中每个像素，利用该像素分别属于孔隙和背景的概率，确定两概率值较大者对应的类别作为该像素最终的类别，即最终为孔隙或者背景，从而对所述数字岩心数据体进行孔隙分割，得到孔隙分割的二值结果。该过程属于现有技术，在此不做具体说明。可以理解的是，所得到的孔隙分割的二值结果中孔隙可以用0表示，显示为黑色，背景可以用1表示，显示为白色。
[0111]
a3，根据所述孔隙分割的二值结果，计算出所述数字岩心数据体在当前次迭代的孔隙度。
[0112]
具体的，可以根据所述孔隙分割的二值结果，利用所述数字岩心数据体中被划分为孔隙的像素数/所述数字岩心数据体中像素总数，计算出所述数字岩心数据体在当前次迭代的孔隙度。
[0113]
可以理解的是，每次迭代后得到的概率分布参数都不相同，因此计算出的孔隙度也不同。
[0114]
可选的一种实施方式中，所述基于所述当前次迭代的孔隙度、上一次迭代的孔隙度，以及对所述数字岩心数据体对应的岩心样本进行实测得到的实测孔隙度，确定迭代速度的待变化类别，包括：
[0115]
b1，利用所述当前次迭代的孔隙度、所述上一次迭代的孔隙度，以及所述实测孔隙度，计算得到当前次迭代的待变化类别初始判别值，并对所述待变化类别初始判别值进行归一化处理，得到当前次迭代的待变化类别判别值。
[0116]
其中，所述当前次迭代为第一次迭代时，所述上一次迭代的孔隙度可以利用em算法在第一次迭代中初始的概率参数计算出一个数值，或者也可以给定一个预设值。
[0117]
该步骤中，所述当前次迭代的待变化类别初始判别值的计算公式，包括：
[0118][0119]
其中，δpsi表示所述当前次迭代的待变化类别初始判别值；psi表示所述当前次迭代的孔隙度；ps
i-1
表示所述上一次迭代的孔隙度；其中，式(6)中i表示当前次迭代的迭代序号；ps
test
表示所述实测孔隙度。其中，所述实测孔隙度是对待孔隙分割的数字岩心数据体所对应的岩心样本预先利用双室孔隙测试仪测量获得的。
[0120]
以下对构思设计出公式(6)的过程予以简要说明。
[0121]
本发明实施例发明人在使用em算法的时候发现，迭代过程较为耗时，导致数字岩心整体处理的效率不高，并且容易出现分割过度情况，影响孔隙分割的精度。因此本发明实施例结合孔隙度，提出基于孔隙度参数调节的方法，提高孔隙分割的精度，加速收敛的速度。由于每一次迭代的过程得到的概率分布参数都能体现出不同的分割效果，进而体现出分割出的孔隙度的变化情况。具体的，在迭代未结束的情况下，与上一次迭代相比，概率分
布参数的变化如果导致待分割切片孔隙度变化较小，且与实测孔隙度差别较大，则说明迭代收敛过程应当加速，概率分布参数的变化幅度应当加大，反之，则说明迭代收敛过程应当减速，概率分布参数的变化幅度应当减小。简而言之，在迭代中基于孔隙度进行概率分布参数调整，从而调整迭代收敛的速度。
[0122]
具体的，考量待分割切片孔隙度变化情况以及与实测孔隙度差别情况的过程中，需要考虑到所述当前次迭代的孔隙度psi、所述上一次迭代的孔隙度ps
i-1
以及所述实测孔隙度ps
test
这三者间的关系，因此，本发明实施例首先建立三者间的孔隙度的变化关系如式(7)和(8)。
[0123][0124][0125]
式(7)表示第i次迭代与第i-1次迭代的孔隙度变化率；式(8)表示第i次迭代后孔隙度与实测孔隙度的逼近程度。当式(7)的结果大于式(8)的结果时，说明孔隙度变化速率高于孔隙度趋近于实测值的速率，应当加速收敛过程；当式(7)的结果小于式(8)的结果时，说明孔隙度变化速率低于孔隙度趋近于实测值的速率，应当减缓收敛过程；当然，如果当式(7)的结果等于式(8)的结果时，说明可以不改变收敛过程的速度。
[0126]
因此，将式(7)和式(8)进行比较得到为了防止分母出现0的情况，为分母加1，最终得到式(6)表示的所述当前次迭代的待变化类别初始判别值δpsi。
[0127]
当|δpsi|》1时，提高概率分布参数的调节速度；当|δpsi|《1时，降低概率分布参数的调节速度。但由于不同岩心的扫描数据中，δpsi的范围并不相同，因此，对所述待变化类别初始判别值δpsi进行归一化处理，得到当前次迭代的待变化类别判别值δpsi'。
[0128]
具体的，所述对所述待变化类别初始判别值进行归一化处理，得到当前次迭代的待变化类别判别值的过程中，所采用的计算公式包括：
[0129][0130]
其中，δps'i表示当前次迭代的待变化类别判别值。
[0131]
可以理解的是，经过式(9)，当前次迭代的待变化类别判别值δpsi'的数值范围被归一化在[0.1,0.9]之间。
[0132]
b2，若所述待变化类别判别值大于0.5，则确定迭代速度的待变化类别为加速；反之，则确定迭代速度的待变化类别为减速。
[0133]
具体的，如果当前次迭代的待变化类别判别值δpsi'》0.5，则确定迭代速度的待变化类别为加速；如果当前次迭代的待变化类别判别值δpsi'《0.5，则确定迭代速度的待变化类别为减速。
[0134]
可选的一种实施方式中，根据所述待变化类别对所述当前次迭代中更新的概率分布参数进行修正，包括：
[0135]
c1，根据所述当前次迭代的待变化类别判别值计算对应的修正权值。
[0136]
可选的一种实施方式中，该步骤采用的计算公式包括：
[0137][0138]
其中，ω表示所述修正权值。
[0139]
式(10)的修正权值是由sigmoid函数坐标变换后推导出来的，其数学波形如图2所示。图2中横轴为δpsi'，纵轴为ω。
[0140]
c2，基于所述修正权值，利用待变化类别对应的调节公式，对所述当前次迭代中更新的概率分布参数中的孔隙类和背景类在高斯概率分布下的均值进行修正。
[0141]
其中，所述待变化类别对应的调节公式，包括：
[0142]mcj
'＝(1 ω)
×mcj
ꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀ
(11)
[0143]
其中，m
cj
表示所述当前次迭代中更新的概率分布参数中的孔隙类和背景类在高斯概率分布下的均值；j＝1,2分别表示孔隙类和背景类；m
cj
'表示对m
cj
进行修正后的数值。
[0144]
可见，根据所述待变化类别对所述当前次迭代中更新的概率分布参数进行修正后，将修正后的结果作为下一次迭代中初始的概率参数，返回s2继续进行迭代。其中，修正后的结果包括所述当前次迭代中更新的概率分布参数中的孔隙类和背景类在高斯概率分布下的方差s
cj
、修正得到的m
cj
'，以及当前次迭代中孔隙和背景的先验概率p
c1
,p
c2
。
[0145]
可见，本发明实施例基于孔隙度变化率和分割孔隙度与实测孔隙度的逼近程度，建立了待变化类别对应的调节公式，能够调整em算法的收敛速度，提高分割精度。
[0146]
若是，执行s6，利用所述当前次迭代中更新的概率分布参数，获得所述数字岩心数据体中各像素的孔隙分割结果。
[0147]
具体的，得到当前次迭代中更新的所述孔隙和背景的概率分布参数m
c1
、s
c1
、m
c2
、s
c2
后，可以利用上述参数计算所述数字岩心数据体每个切片中的像素属于孔隙和背景的概率。针对每个像素，以确定出的较高概率值确定其像素类型，实现孔隙分割。具体过程请结合相关现有技术理解，在此不做详细说明。
[0148]
数字岩心一般为灰度图像，其中孔隙对应图像中较暗即灰度值较低的部分，单张数字岩心中，内容简单，孔隙特征比较明显，但是孔隙结构复杂无固定形态和特征、孔径范围跨域微米到纳米、不同ct数据中其对比度与饱和度也不相同，同时，考虑到实际工程问题，对整个数字岩心数据的处理耗时也是一个必要的考虑因素。
[0149]
本发明实施例所提供的基于em与孔隙度的数字岩心孔隙分割方法中，首先，为避免高灰度值矿物的影响，基于灰度分布特征从待孔隙分割的数字岩心数据体中选择满足预设要求的预设数量个样本，使得样本的均值和方差分别对应小于所述数字岩心数据体的均值和方差。同时，该种样本筛选方式使得em迭代过程所需要的样本从待分割的数岩心数据体自身采集，避免了样本库数据不匹配的问题，能够提高算法的通用性。然后，对预设数量个样本采用em算法进行无监督的学习，针对每次迭代，利用初始的概率参数得到当前次迭代中更新的概率参数，在满足迭代停止条件时，利用当前次迭代中更新的概率分布参数，获得所述数字岩心数据体中各像素的孔隙分割结果。在不满足迭代停止条件时，利用所述当前次迭代中更新的概率分布参数，计算当前次迭代的孔隙度，并结合实验获得的实测孔隙度确定迭代加速还是减速，并相应对所述当前次迭代中更新的概率分布参数进行修正，以
得到下一次迭代中初始的概率参数继续进行迭代。本发明实施例通过在参数学习过程中结合实测孔隙度精细调节概率分布参数的变化，以确定最终的概率分布参数，能够避免参数调节过度，并能够加速迭代进程，降低运算耗时。因此在分割数字岩心的孔隙时能够提高孔隙自动分割的精度和收敛的速度。
[0150]
为了验证本发明实施例方法的有效性，以下以实验结果进行说明。
[0151]
具体采用本发明实施例方法对采集到的数字岩心进行孔隙的提取，并与em、otus算法进行对比。otus(称为大津法/最大类间方差法)是图像分割领域的经典算法，也是目前在数字岩心工程处理中最常采用的孔隙提取方法之一。该实验环节从工程可实现化的角度出发，将本发明实施例方法与两种典型算法进行比较分析。
[0152]
1)学习样本采集效果对孔隙提取的影响
[0153]
本发明实施例方法从待孔隙分割的数字岩心数据体中采集学习样本，挑选不含有较高灰度值的矿物的样本作为学习样本。当样本选取不合适时会对分割结果造成影响。如图3为对相同的待分割的数字岩心切片图像用不同样本学习分割时的影响对比。其中图3(a)为待分割的数字岩心切片图像；图3(b)从左至右依次为样本1、样本2、样本3和样本4对应的图像。其中，样本1、样本2是本发明实施例选出的满足预设要求的样本；样本3中存在高灰度值的矿物且孔隙较少；样本4局部存在中等灰度值的明显矿物且孔隙较少。
[0154]
图3(c)～图3(g)分别是这四个样本经过经典em算法迭代1次、6次、11次、16次和21次的分割结果。图3(c)～图3(g)中每行中从左至右分别对应样本1、样本2、样本3和样本4在相应迭代次数下的分割结果。
[0155]
从图3可以看出，训练样本为样本1和样本2时，经过经典em算法迭代运算后，分割结果是比较相似的。样本3中存在高灰度值的矿物，且孔隙比较少，因此导致当经典em算法迭代到第6次时孔隙分割已经严重失真，但是根据经典em算法理论，即使结果错误但是未满足迭代停止条件仍会继续迭代运算下去。样本4局部存在中等灰度值的明显矿物，且孔隙较少，从分割的效果看，进行到第11次迭代后，分割结果也开始出现明显失真，接下来的迭代运算也并未修正错误，反而加重了失真。
[0156]
可见，上述实验验证了本发明实施例从待孔隙分割的数字岩心数据体中选择满足预设要求的样本的必要性，通过本发明实施例的样本筛选能够提高分割效果。
[0157]
采用本发明实施例方法从待孔隙分割的数字岩心数据体中自动提取样本，并筛选出真正有效的样本，图4为本发明实施例实验给出的样本筛选示例；其中图4(a)为本发明实施例方法筛选出的满足预设要求的3个典型学习样本；图4(b)为本发明实施例方法筛选出的3个不合格样本。
[0158]
2)孔隙提取结果
[0159]
针对本发明实施例选取出的3个样本，如图5(a)所示，分别采用otus、经典em算法和本发明实施例方法进行处理，得到的孔隙提取效果，分别如图5(b)、图5(c)和图5(d)所示。
[0160]
其中，图5(a)为3个样本，也就是3个待分割切片，分别以a-1、a-2和a-3表示。图5(b)为样本a-1、a-2和a-3采用otus方法进行孔隙提取的效果图；图5(c)为样本a-1、a-2和a-3采用经典em算法进行孔隙提取的效果图；图5(d)为样本a-1、a-2和a-3采用本发明实施例方法进行孔隙提取的效果图。
[0161]
从图5可以看出，图5(a)的a-1含有较多的特征明显的大尺寸孔隙，高亮度矿物分布相对比较均匀，三种方法均能够比较好地提取孔隙；图5(a)的a-2含有明显的孔隙和分布不均匀的高亮度矿物，otus方法提取结果中包含较多微小孔隙和噪声，相比之下经典em方法和本发明实施例方法能更好地够克服噪声的影响；图5(a)的a-3含有较少孔隙且形态特征复杂，矿物的亮度和分布也极不均匀，从提取效果看，经典em和本发明实施例方法则能有效地实现孔隙的分割，而otus方法则对高亮度矿物进行了识别，孔隙提取失败。
[0162]
也就是说，当切片中出现高灰度值矿物，且孔隙相对较少时采用otus方法会导致分割明显失真，而经典em算法和本发明实施例方法针对这类样本则可有效避免这个问题。
[0163]
对三个岩心样本采用上述三种方法进行处理，得到单个切片的孔隙度如表1所示。
[0164]
表1.孔隙度对比
[0165][0166]
表1中采用实测孔隙度为参考标准。实测孔隙度是对岩心样本整体进行实验测试得到的，用于处理的数字岩心是在扫描数据体的基础上去除扫描影响后截取的部分，因此数字岩心孔隙提取后的孔隙度本来就和岩心的实验孔隙度存在偏差，但多个样本进行比较时整体趋势是一致的。从表1中可以看出采用otus方法计算出的孔隙度变化率过高，与em算法和本发明实施例方法对比可以看出，后两种方法分割后的趋势一致，能吻合实际孔隙的测试结果。
[0167]
关于em算法和本发明实施例方法迭代过程的对比可以参见表2～表4。表2～表4分别是针对三个数字岩心切片图像，采用em算法和本发明实施例方法迭代时，迭代次数以及对应的孔隙度，其中，表格中迭代次数下的小数数值是该次迭代得到的孔隙度。\表示未进行对应次迭代。
[0168]
表2 slice1迭代过程对比
[0169][0170]
表3 slice2迭代过程对比
[0171][0172]
表4 slice3迭代过程对比
[0173][0174]
从表2～表4可以看出，本发明实施例方法能够在保障分割效果的同时以更快的速度进行收敛，从对三个slice的处理数据来看，相比em算法，本发明实施例方法能够平均减少2次迭代，处理速率平均提高了28％。
[0175]
通过将本发明实施例方法与经典em算法和otus算法的处理效果进行比较，表明本发明实施例基于统计特征的样本采集方法不需要提前建立样本库，能有效提高分割的准确性，避免高亮矿物对样本的影响，且对不同类型样本都能取得良好的分割效果。同时，本发明实施例相比于传统em算法，迭代效率提升近似28％，非常有利于数字岩心分析的工程应用，对数字岩心的分割与基于数字岩心的其余研究均有促进意义。
[0176]
以上所述仅为本发明的较佳实施例而已，并非用于限定本发明的保护范围。凡在本发明的精神和原则之内所作的任何修改、等同替换、改进等，均包含在本发明的保护范围内。