岩心图像自适应分区三维重建方法与流程

本发明涉及一种分区重建方法，尤其涉及一种岩心图像自适应分区三维重建方法，属于三维图像重建技术领域。

背景技术：

岩心多孔介质作为油气能源在自然界中最重要的存储空间之一，其孔隙空间的形态(几何形状和连通性)强烈地影响着储层中油气的储存和迁移。因此，对于岩心孔隙微观结构的研究逐渐成为人们关注的焦点。

目前获取岩心孔隙结构的方法主要有两种：一是通过计算机断层扫描(ct)以及聚焦离子束扫描电镜(fib-sem)等设备直接构建真实的三维数字岩心；二是利用二维图像进行三维数学建模。第一种方法由于扫描设备比较昂贵，还不能普遍使用。此外，由于高分辨率图像能够揭示有关孔隙网络和孔隙连通性方面的重要细节，对于岩心特性研究具有重要意义，但受扫描设备成像机理的限制，扫描图像的分辨率与图像尺寸是相矛盾的，要获取微米或纳米级的大视域3d孔隙结构(如页岩等致密岩)的清晰图像，是昂贵且耗时的。而2d高分辨率图像则能在较大规模下以更低的成本和更高的效率获取。鉴于此，基于二维图像的岩心三维重建仍然是一个重要的研究方向。

现有的三维重建算法大都基于三维结构是均质且各向同性的假设，然而在实际应用中，岩心微观结构大部分都是非均质、各向异性的。要想获得岩样较完整的孔隙分布情况，需要得到大视域的高分辨率图像，而这类图像通常整体呈现较强的非均质性，且图像尺寸较大，以1cm×1cm的岩样为例，图像尺寸在分辨率为1μm/像素的情况下高达10000×10000，而在更高的nm级别分辨率条件下，扫描得到的图像尺寸更大。而图像尺寸越大，对重建时间与计算机内存的要求越高，直接对这类图像进行三维重建仍是三维重建研究中的一大难点，相关研究较少。

技术实现要素：

本发明的目的就在于解决大视域岩心图像难以直接进行三维重建，难以获取岩样完整结构信息这一难题，而提供的一种对岩样包含的各类结构信息分别进行重建并整合的岩心图像自适应分区三维重建方法。

本发明通过以下技术方案来实现上述目的：

1、本发明所述岩心图像自适应分区三维重建方法包括以下步骤：

(1)计算待划分区域的最佳局部孔隙度模板尺寸；

(2)使用最佳模板对待划分区域所有像素点计算其对应的局部孔隙度；

(3)以待划分区域的孔隙度p为依据，将局部孔隙度大于p的像素点分为一个区域，剩余的像素点归为另一个区域；

(4)根据各区域的面积以及其中最大最小局部孔隙度差异，判断现有的区域是否还需要进行划分；

(5)对还需要划分的区域重复步骤(1)～(4)，直到没有区域需要划分为止；

(6)计算每个区域各自的孔隙度及相对于整幅图像的面积占比，并以这些孔隙度为依据提取每个区域的典型代表单元体rev；

(7)使用传统的三维重建算法对所有rev进行三维重建；

(8)计算所有重建三维孔隙结构的孔喉参数，并将这些参数按照(6)中计算的面积占比进行整合，从而得到整个岩样的三维结构参数。

上述方法的基本原理如下：

实际应用中，岩心图像虽然整体呈现的非均质性较强，但图像不同区域的内部往往孔隙尺寸与分布较为均匀，具有均质性。因此可以将孔隙的尺寸与分布作为区域划分的关键。而局部孔隙度分布曲线是用来判断岩心图像均质特性的重要方法之一，它通过使用固定模板扫描整幅训练图像，计算模板内的孔隙度，即图像局部区域的孔隙度，并对局部孔隙度出现的概率以及分布情况进行统计，从而获取图像的均质或非均质特性。一般来说，图像的局部孔隙度分布范围越窄，表明岩心的均质性越强。因此，局部孔隙度可以在一定程度上反映孔隙的尺寸与分布情况。通过计算以每个像素点为中心的局部区域孔隙度，可以表征该像素点周围的孔隙特征，将具有相近局部孔隙度的像素点划分为同一个区域，则区域的局部孔隙度曲线分布范围较窄，表明区域内部呈现均质性。且考虑到划分后各区域形状不规则，可以根据区域的孔隙度来获取每个区域规则的典型代表单元体rev，将它们作为各区域的训练图像分别进行重建，最后根据每个区域的面积占比对重建的三维结构相关参数进行整合，将整合结果作为整幅图像对应三维结构的相关参数。

具体地，所述步骤(1)中，局部孔隙度定义如下：

局部孔隙度在计算时首先设定一个测量单元测量单元在二维图像中是边长为l的正方形，在三维图像中是边长为l的正方体，其中，n(g)表示g相在测量单元中占的比例，对于二维图像，n(g)表示g相的面积，对于三维图像n(g)表示g相的体积，m为测量单元的数量，δ(x)为狄拉克函数；

最佳局部孔隙度模板的边长根据待划分区域的孔隙尺寸与分布决定，具体将图像中x、y方向上，孔隙与背景弦长的最大值作为最佳模板尺寸m_size，以保证模板在遍历图像过程中不会出现全为孔点或全为背景点的情况，从而更好的反映孔隙的分布情况；

所述步骤(2)中，利用最佳模板计算待划分区域每个像素点的局部孔隙度p_local，该值通过将模板中心与像素点重合，计算模板内的孔隙度来得到；

所述步骤(3)中，以整个待划分区域的孔隙度p_area为阈值，若待划分区域中像素点的p_local<＝p_area，将满足该条件的像素点划分为一区，满足p_local>p_area的像素点自动归为另一区；

所述步骤(4)中，首先判断现有的这些区域面积是否足够大，对于面积够大的区域，计算该区域像素点最大局部孔隙度与最小局部孔隙度之间的差异，若差异较大，表明该区域具有非均质性，需要继续进行划分；

所述步骤(5)中，对于判断出还需要继续划分的区域，因为区域的孔隙尺寸与分布同之前的区域会有差异，最佳模板尺寸也会变化，因此需重新计算各区域的最佳局部孔隙度模板尺寸，重复步骤(1)～(4)，直到没有区域需要划分；

所述步骤(6)中，首先需要计算划分后这些区域相对于整幅图像的面积占比，由于各区域的孔尺寸与分布都较为均匀，但是区域形状不规则，需要从这些区域中选取规则的典型代表单元体rev，为了使选取的rev更具代表性，选择将区域的孔隙度作为标准，使rev与它代表的区域孔隙度相近；

所述步骤(7)中，将各区域的rev作为训练图像，使用现有的三维重建算法，比如模拟退火算法对训练图像进行三维重建；

所述步骤(8)中，计算所有重建三维结构的孔喉参数，并将这些参数按照其代表区域的面积占比进行整合，将整合结果作为整幅图像对应三维结构的结构参数。

本发明的有益效果在于：

本发明充分考虑了真实岩心图像的孔隙分布情况，不需要直接对整幅岩心图像进行重建，而是利用孔隙特征对图像进行区域划分，使区域内部呈现均质性。由于同一区域内的孔隙尺寸与分布都具有相似性，可以从中选择部分具有代表性的区域进行重建，从而大大减少重建数据量，解决岩心图像尺寸过大时难以直接进行三维重建的难题。通过对各区域的重建结果按照区域面积占比进行整合，能够在一定程度上反映整幅岩心图像对应三维结构的孔隙特性。

附图说明

图1是真实岩心图像；

图2为图1真实岩心图像对应的孔隙图形层；

图3是图像进行区域划分后的分区情况；

图4是各区域选取的典型代表单元体rev；

图5为各区域rev重建后的三维结构；

具体实施方式

下面结合附图对本发明作进一步说明：

(1)图1为原始岩心图像，该图像的尺寸为5984×5878，分辨率为2.277μm/像素，首先对其进行孔隙提取，获取二值的岩心孔隙图形层，如图2所示。

(2)计算待划分区域的孔隙度及其对应的最佳模板尺寸，使用最佳模板对图像进行遍历，从而获取图像每个像素点对应的局部孔隙度。第一次计算得到整幅图像的孔隙度为1.705％，最佳局部孔隙度模板尺寸为321。

(3)将待划分区域内局部孔隙度值小于区域孔隙度的像素点分为一区，其余像素点为另一区。

(4)判断划分后的区域是否还需要进行划分。由于图像尺寸较大，这里设置当区域面积大于1000×1000，且区域内像素点的最大最小局部孔隙度差异大于20％时，则该区域还需要继续进行划分。

(5)重复步骤(2)～(4)，直到没有区域需要进行划分，此时整幅图像被划分为3区，为了便于观察，将不同的区域赋予不同的颜色并输出着色后的图像，如图3所示。

(6)计算不同区域的孔隙度以及它们相对于整幅图像的面积占比，将结果展示如表1所示。

表1

(7)根据各区域的孔隙度，从其中选取孔隙度与区域孔隙度相近的规则区域，作为各区域的典型代表单元体rev，即训练图像，各区域选取的训练图像如图4所示，从左到右依次为浅灰区域、中灰区域以及深灰区域的训练图像，尺寸均为512×512，孔隙度分别为0.263％、1.051％、3.365％。

(8)使用稳相分级模拟退火算法对这些训练图像进行三维重建，结果如图5所示，它们的展示顺序与图4的训练图像一致。计算每个结构的孔喉参数，结果如表2所示。

表2

(9)对步骤(8)中计算得到的各区域三维结构孔喉参数按照区域的面积占比进行整合，结果如表3所示。

表3

通过观察岩心图像孔隙图形层中孔隙的大小与分布，与自适应分区后的情况比较，可以看到具有相似孔隙尺寸与分布的区域被划分为同一区，说明本方法能够较好的实现岩心图像的自适应分区，使整体呈现非均质性的图像，经过分区后区域内部具备均质性，从而可以使用局部孔隙信息代表更大区域的孔隙特征，大大减少重建的计算量以及重建时间。

上述实施例只是本发明的较佳实施例，并不是对本发明技术方案的限制，只要是不经过创造性劳动即在上述实施例的基础上实现的技术方案，均应视为落入本发明专利的权利保护范围内。