CN105139444B - 基于岩心二维颗粒图像的三维颗粒结构重建方法 - Google Patents

Abstract

本发明涉及一种基于岩心二维颗粒图像的三维颗粒结构重建方法。即基于岩心二维图像的三维建模问题，利用二维图像中包含的矿物颗粒信息，推断与二维颗粒结构相对应的三维颗粒结构信息。具体提出基于模拟退火算法与颗粒生长算法相结合的岩心颗粒微观结构重建算法，使重建的三维颗粒结构与二维颗粒结构参考图像具有相似的形态分布，并能够很好的捕获岩心颗粒结构特征；其微观结构中矿物颗粒的三维分布不仅展现了矿物颗粒空间分布，且在三维颗粒结构基础上量化不同矿物成分对其结构性能的影响；重建的三维微观结构为真实岩心微观结构提供了一个很好解释；可用于岩心微观结构电学特性及渗流特性的研究，具有实用价值。

Description

基于岩心二维颗粒图像的三维颗粒结构重建方法

技术领域

本发明涉及一种基于二维图像的三维建模方法，尤其涉及一种基于岩心二维颗粒图像的三维颗粒结构重建方法，属于三维图像重建技术领域。

背景技术

在石油地质研究中，岩心三维微观结构是研究岩心宏观物理特性的基础。以岩心二维颗粒图像为基础进行岩心三维颗粒结构的重建，是获取岩心三维颗粒结构信息的主要手段之一。目前，基于岩心二维图像重建与其对应的三维结构的方法，大多侧重于孔隙的三维重建，而对岩心中矿物颗粒的三维空间分布的研究相对较少。在目前缺少对岩心三维颗粒微观结构研究的情况下，研究者基于岩心二维颗粒图像的三维颗粒结构重建，重建的三维颗粒微观结构为真实岩心微观结构提供了一个很好的解释。目前对岩心二维颗粒图像重建与其对应的三维颗粒结构的重建算法主要有基于过程的重建算法和最大熵原理重建算法等。

所述基于过程的重建算法的基本思想是模拟岩心的沉积过程。用该方法构建岩心三维颗粒结构主要包括三个过程：颗粒沉积过程、压实过程和成岩过程。在构建数字岩心三维颗粒结构过程中，假定岩心颗粒粒径分布满足正态分布，且岩心颗粒形状均为球形。此时，颗粒在沉积过程中遵循重力势能最小原理，主要利用下降与滚动算法确定三维数字岩心中每一岩心颗粒球心的稳定位置；然后将所有沉积岩心颗粒向下移动来模拟岩心的压实和颗粒重排过程；最后模拟岩心矿物的地质形成过程，在成岩过程模拟中，考虑了石英胶结质的生长和粘土物质的填充作用。

1984年美国电气工程师C.E.Shannon在其《通信的数学理论》一文中首次提到了信息熵的概念，把熵作为一个随机事件的“不确定性或信息量的度量。”而最大熵算法是一种选择随机变量统计特性最符合客观情况的准则，在这种情况下，概率分布最均匀，预测的风险最小。因为这时概率分布的信息熵最大，所以人们称这种算法叫“最大熵算法”。所述最大熵原理重建算法构建微观结构的核心思想是把微观结构看成随机场并利用最大熵原理构造随机场的概率分布函数，即在已知二维图像随机场的概率分布函数的状况下估计重建三维结构最大可能的概率分布，使重建三维结构的概率分布函数接近原始二维图像的微观结构；并能够将任意多的概率信息并入到重建后的微观结构模型中，为匹配目标提供最优选择是最大熵原理重建算法的优点。但其最大缺点是需要大规模的计算，因此重建微观三维颗粒结构非常耗时且所重建效果并不理想。应用最大熵原理研究微观结构的重建主要是针对平面微观结构，而目前基于最大熵原理的岩心三维颗粒结构重建工作未见报道。因此，如何解决这一问题，是所属技术领域科技工作者面临急待解决的问题。基于此，本课题组潜心研究，提出一种基于岩心二维颗粒图像的三维颗粒结构重建方法；该研究项目受国家自然科学基金项目《岩石微观非均质结构三维图像重建及分辨率提升技术研究》(61372174)资助。

发明内容

本发明的目的在于为了解决上述现有技术中所存在的技术问题，而提供一种基于岩心二维颗粒图像的三维颗粒结构重建方法。该方法是基于岩心二维图像的三维建模问题，利用二维颗粒图像中所包含的矿物颗粒信息，推断与二维矿物颗粒结构信息相对应的三维颗粒结构信息，采用重建算法，使重建的三维颗粒结构与二维颗粒结构参考图像具有相似的形态分布。

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

本发明提供的一种基于岩心二维颗粒图像的三维颗粒结构重建方法，包括以下步骤：

(1)提取原始岩心二维颗粒图像中矿物颗粒的个数，对不同矿物颗粒赋予不同的灰度级，将处于不同灰度级的矿物颗粒进行标号，以区分出不同的颗粒；

(2)统计步骤(1)提取的岩心二维颗粒图像中每个矿物颗粒的面积；采用遍历整幅岩心二维颗粒图像，用计算机记录每个矿物颗粒标号下二维矿物颗粒像素值的个数，即统计出岩心二维矿物颗粒的面积；

(3)将步骤(2)统计出的每个矿物颗粒的面积，按其大小给出二维矿物颗粒面积分布直方图；

(4)对步骤(3)给出的每个二维矿物颗粒面积分布直方图区间的二维矿物颗粒，采用第一阶段模拟退火算法来预测其直方图区间中每个岩心二维矿物颗粒的三维颗粒结构；

(5)将步骤(4)预测的每个岩心二维矿物颗粒的三维颗粒结构放入重建的三维空间中；整个重建过程当满足预测终止条件时结束，即获得每个岩心二维矿物颗粒的三维颗粒结构；

(6)对步骤(5)在第一阶段预测重建所得每个岩心二维矿物颗粒三维结构的基础上采用第二阶段颗粒生长算法预测重建，即在余下的未被颗粒占据的空间中，应用颗粒生长算法模拟岩心颗粒的生长过程，整个生长过程直到满足预测终止条件时结束，即完成岩心二维颗粒图像的三维颗粒结构重建。

上述方案中，步骤(3)中所述统计出的每个矿物颗粒的面积，按其大小给出二维矿物颗粒面积分布直方图，即找出其中的最大面积和最小面积，通过最大面积与最小面积的差值除以10，即将二维矿物颗粒面积分成10等份，计算出二维矿物颗粒面积分布直方图每个区间的大小，再将每个颗粒按其面积大小存放到不同的各区间中。

上述方案中，步骤(4)中所述通过第一阶段模拟退火算法预测直方图区间中每个岩心二维矿物颗粒的三维颗粒结构时，首先提取该岩心二维矿物颗粒，然后在三维空间中以该岩心二维矿物颗粒所占的体积分数随机布点作为初始结构，并以岩心二维矿物颗粒的两点概率函数作为重建约束条件，预测三维颗粒结构。

上述方案中，步骤(5)中所述将预测的岩心二维矿物颗粒的三维颗粒结构放入重建的三维空间中，是在重建每一个三维颗粒后，在存储颗粒的三维空间中随机寻找没有被颗粒占据且能够完整存放三维颗粒的空间点放置重建的三维颗粒；当重建三维颗粒结构的直方图与二维颗粒图像相匹配或重建三维颗粒结构占整个三维空间的百分比与二维颗粒图像相同时，重建达到终止条件。

上述方案中，步骤(6)中所述的第二阶段颗粒生长算法是在步骤(5)中重建结束后，在余下的三维空间中，随机选取一个没有被矿物颗粒占据的空间点作为矿物颗粒生长核，以该生长核为中心向其18邻域生长；该生长核生长后的结构作为新的矿物颗粒生长核，重复上述过程，直到重建三维颗粒结构的直方图与二维颗粒图像相匹配，或重建三维颗粒结构占整个三维空间的百分比与二维颗粒图像相同时，重建即结束。

上述方案中，所述步骤(6)中所述的颗粒生长算法，是在颗粒生长过程中，如果生长的颗粒和其他颗粒结构相接触则停止该方向上颗粒的生长，完成一个颗粒的生长之后，随机再选取一个没有被矿物颗粒占据的空间点作为新的颗粒生长核，重复上述过程，直到颗粒生长的体积与预测的颗粒结构的体积相同为止，即完成所述三维颗粒结构的预测。

上述方案中，所述以生长核为中心向其18邻域生长，它的最近邻域为6邻域和12次近邻，其生长概率的比值为4:1，生长核生长后的结构作为新的生长核，重复上述过程，直到颗粒生长的体积与预测的颗粒结构的体积相同为止，即完成所述三维颗粒结构的预测。

上述方案中，所述生长概率的设置是基于矿物颗粒结构是一种聚合物的形状，因而生长核的最近6邻域的生长概率高于12次近邻的生长概率，本发明中将生长核的最近6邻域和12次近邻的生长概率的比值设置为4:1。

所述模拟退火算法是一种求解优化组合问题的通用算法，其原理是将可行解和目标函数看作是物质体系的状态和能量，通过模拟物质体系的退火过程来求优化问题的最优解的过程。所述颗粒生长算法借鉴自多晶材料的晶粒生长过程。随着材料科学的发展，人们逐渐认识到多晶材料的物理、化学性能不仅依赖于材料的成分，而且在很大程度上也依赖于晶粒生长过程中的晶粒尺寸和晶粒组织结构。晶粒生长过程主要模拟多晶材料的晶粒演化过程。对晶粒长大现象的研究可分为两大类：1.晶粒长大动力学的研究，即研究晶粒尺寸的变化；2,晶粒长大过程中晶粒组织形态的研究，如晶界形状、晶界夹角等。

本发明基于岩心二维颗粒图像重建与其对应的三维颗粒结构的方法的基本原理如下：

本发明的方法所述基于岩心二维颗粒图像的三维颗粒结构建模方法，关键在于依据给定的岩心二维颗粒图像中矿物颗粒的结构进行预测，通过模拟退火算法重建与岩心二维颗粒图像中每个矿物颗粒相对应的三维颗粒结构，然后将预测的矿物颗粒的三维颗粒结构放入重建的三维空间中，并在余下的未被颗粒占据的三维空间中，进行颗粒生长重建。即以岩心二维颗粒图像为参考图像，预测岩心二维矿物颗粒的三维颗粒结构及其空间分布。所述预测包含两个阶段，第一阶段：通过模拟退火算法预测每一岩心颗粒的三维颗粒结构，并将预测的三维颗粒结构置于重建的三维空间中，随着放置的岩心颗粒越来越多，三维空间中可放置岩心颗粒的空间越来越小，直到三维空间中再也不能够放下预测的三维颗粒结构时，重建的第一阶段结束；第二阶段：在模拟退火算法预测每一岩心颗粒三维颗粒结构基础上应用颗粒生长算法，在余下的未被颗粒占据的三维空间中，随机选取一个没有被矿物颗粒占据的空间点作为矿物颗粒生长核，然后以该矿物颗粒生长核为中心向其18邻域生长，矿物颗粒生长核生长后的结构作为新的矿物颗粒生长核，重复上述过程，直到重建三维颗粒结构的直方图与二维颗粒图像相匹配或重建三维颗粒结构占整个三维空间的百分比与二维颗粒图像相同时，重建达到终止条件。所述以生长核为中心向其18邻域生长，它的最近邻域为6邻域和12次近邻，其生长概率的比值为4:1，生长核生长后的结构作为新的生长核，重复上述过程，直到颗粒生长的体积与预测的颗粒结构的体积相同为止。

本发明与现有技术相比具有以下的优点及有益的技术效果：

本发明提出的基于岩心二维颗粒图像的三维颗粒结构重建方法，该方法采用基于模拟退火算法与颗粒生长算法相结合的岩心二维矿物颗粒结构的预测方法。在缺少岩心三维颗粒结构研究的情况下，基于岩心二维颗粒图像的三维颗粒结构重建，重建的三维颗粒微观结构为真实岩心微观结构提供了一个很好的解释。本发明以岩心二维颗粒图像为基础，获取其颗粒结构的二维分布，并以获取的颗粒结构二维分布为基础，重建与二维颗粒结构分布相一致的三维颗粒结构。通过比较重建的三维颗粒结构不同切面与二维参考图像中颗粒的粒径分布、形状因子、岩心定名等特性参数，表明本发明所提出的基于模拟退火算法与颗粒生长算法相结合的岩心颗粒微观结构重建算法能够很好的捕获岩心颗粒的结构特征。通过实施例，岩心微观结构中矿物颗粒的三维分布不仅展现了矿物颗粒的空间分布，而且能在三维颗粒结构基础上量化不同矿物成分对其结构性能的影响；其重建的三维颗粒微观结构为真实岩心微观结构提供了一个很好的解释；可以应用于岩心微观结构电学特性及渗流特性的研究，为该研究提供了基础数据和理论依据，具有岩心研究的实用价值。

附图说明

图1-1是本发明实施例中给定的原始岩心二维矿物颗粒图像；

图1-2是本发明实施例中从图1-1中提取的二维矿物颗粒图；

图1-3是本发明实施例中原始岩心二维颗粒图像的孔隙分布图；

图1-4本发明实施例中统计的二维矿物颗粒面积分布直方图；

图2-1是本发明实施例中提取的岩心二维矿物颗粒参考图；用于说明模拟退火算法与其对应的三维颗粒的重建过程；

图2-2是本发明实施例中依据矿物颗粒随机产生的三维颗粒结构的二维切面图；

图2-3是本发明实施例中模拟退火算法重建三维颗粒结构的二维切面图；

图2-4是本发明实施例中通过模拟退火算法重建三维颗粒结构进行去噪和复原的二维切面结构；

图3-1是本发明实施例中提取的岩心二维矿物颗粒参考图，用于说明颗粒生长算法与其对应的三维颗粒的重建过程；

图3-2是本发明实施例中随机选取一个没有被矿物颗粒占据的空间点作为矿物颗粒生长核；

图3-3是本发明实施例中颗粒生长核在二维切面中的八个生长方向；

图3-4是本发明实施例中通过颗粒生长算法生成的颗粒某一切面的最终结果图；

图4为图1-1原始岩心矿物颗粒经过模拟退火算法和颗粒生长算法后重建的三维颗粒结果图；

图5-1、图5-2、图5-3分别为从重建三维颗粒结构图4中选取的三张切面图，通过比较这三张切面图像与原始岩心二维颗粒图像的粒径分布特性来说明本发明所提出的方法的有效性；

图6-1为图5-1中孔隙的分布图；图6-2为图5-2中孔隙的分布图；图6-3为图5-3中孔隙的分布图，通过比较这三维孔隙分布图与原始岩心二维颗粒图像中孔隙结构的孔隙度来说明所提出的算法的有效性。

具体实施方式

下面用具体实施例并结合附图对本发明作进一步详细说明，但所述实施例只是对本发明的实现方法作一个具体的详细说明，而不应理解为是对本发明保护内容的任何限制。

实施例：

为了使本发明所述基于岩心二维颗粒图像重建与其对应的三维颗粒结构的方法更加便于理解和接近于真实应用，下面对提取原始岩心二维颗粒图像中矿物颗粒的个数，通过统计二维矿物颗粒的面积，并按统计面积的大小给出面积分布直方图，再通过模拟退火算法和颗粒生长算法两个阶段的预测获得三维颗粒结构，将三维颗粒结构放入重建的三维空间，最终完成重建二维矿物颗粒对应的三维颗粒结构，以及对重建的三维颗粒结构进行分析的一系列过程的操作流程作整体说明。其中包括本发明的模拟退火算法所重建的核心颗粒预测方法。

其具体操作步骤如下：

(1)对给定的原始岩心二维矿物颗粒图像如图1-1所示，其图像大小为200×200个像素，点长度为6.7μm，提取其中的二维矿物颗粒如图1-2所示，其中黑色表示孔隙，其余不同的颜色代表不同的岩心矿物颗粒，对不同矿物颗粒赋予不同的灰度级，将处于不同灰度级的不同矿物颗粒进行标号；具体标号过程为：遍历图像的所有像素，读取其像素值及该像素的左邻域像素与上邻域像素的像素值，如果该像素左邻域和上邻域像素值都与该像素值不同，判断该像素属于一个新的颗粒，因此给该像素一个新的标号。如果该像素左邻域和上邻域像素值其一与该像素值相同，则判断该像素和这一相同邻域属于同一个颗粒，此时将这一邻域像素值的标号赋予该像素。如果该像素左邻域和上邻域像素值都与该像素值相同，则判断该像素和这两邻域属于同一个颗粒，此时将左邻域的标号赋予该像素和上邻域的像素。此时1-2图中颗粒的体积分数为φ₁＝0.937，孔隙的体积分数为φ₂＝0.063，岩心二维颗粒图像中矿物颗粒之间形成的孔隙分布显示在图1-3中所示。

(2)统计步骤(1)提取的岩心二维颗粒图像中每个矿物颗粒的面积大小；具体计算方法及计算过程为：采用遍历整幅岩心二维颗粒图像，在计算机内存空间中记录统计每一个标号下的二维矿物颗粒的像素值的个数，而每个标号下二维矿物颗粒的像素值的个数即代表其面积，即可统计出岩心二维矿物颗粒的面积。

(3)将步骤(2)计算出的每个矿物颗粒的面积，按其大小给出二维矿物颗粒面积分布直方图；其二维矿物颗粒面积分布直方图如图1-4所示；并依据整幅图像颗粒面积的分布，将面积划分为不同的区间，每个颗粒依据面积大小被归到不同的区间；对于二维颗粒图像来说，其直方图指的是每一区间中颗粒的面积占整体图像的百分比；而对于三维颗粒结构，其直方图指的是每一区间中颗粒的体积占整个三维空间的体积百分比；本实施例在研究岩心颗粒重建时，将二维矿物颗粒的面积划分为十个不同区间。即用最大面积与最小面积的差值除以10，将二维矿物颗粒的面积分成10等份，就是划分为十个不同区间。

(4)对步骤(3)给出的每个二维矿物颗粒面积分布直方图区间的二维矿物颗粒，采用模拟退火算法来预测每个直方图区间中岩心二维矿物颗粒相对应的三维颗粒结构；

在重建岩心矿物颗粒三维颗粒结构时，先从直方图中包含二维颗粒面积最大的区间开始重建，完成该区间颗粒的重建后，依次按其面积减小的直方图区间来重建三维空间中矿物颗粒的三维颗粒结构，直到完成最小尺寸颗粒的重建；如果在直方图某一区间重建过程中，重建该区间的每个颗粒都不能在三维空间中找到一个合适的位置来存放，此时认为该区间颗粒重建已完成，继续进行下一区间颗粒的重建；如果从最小尺寸开始重建，放置在三维空间中的小颗粒把重建空间弄的很零碎，容易导致后面重建的大颗粒在三维空间中找不到一个合适的存放位置。

(5)所述模拟退火算法作为预测矿物颗粒三维颗粒结构第一阶段预测重建：将步骤(4)预测的每个岩心二维矿物颗粒的三维颗粒结构放入重建的三维空间，通过模拟退火算法预测每一岩心颗粒的三维颗粒结构，直到满足预测终止条件为止，即获得每个岩心二维矿物颗粒的三维颗粒结构；

所述模拟退火算法第一阶段预测三维颗粒结构的终止条件为直方图匹配和拥挤度，该阶段重建结束条件为：直方图匹配和拥挤度。对于给定的岩心二维颗粒图像，统计其中每一个二维颗粒所占的面积，并依据整幅图像颗粒大小的分布，将其面积划分为不同的区间，每个颗粒依据其面积大小被归为不同的区间；对于二维颗粒图像来说，直方图指的是每一区间中颗粒的面积占整体图像的百分比；而对于三维颗粒结构则指的是每一区间中颗粒的体积所占整个三维空间的百分比；拥挤度指的是岩心矿物颗粒成分占整个结构的百分比，即岩石相像素点的个数占整个二维颗粒图像或三维颗粒结构像素点个数的比例。拥挤度用来保证重建的颗粒所占的百分比整体上与岩心二维颗粒图像中各种颗粒成分所占的比重相一致性。这两个参数共同作为矿物颗粒重建的约束条件，只要满足其中一个条件，就能完成整个矿物颗粒的三维颗粒结构重建。在预测单个颗粒三维颗粒结构时，依据该颗粒的结构特征和面积，在三维空间中以该颗粒所占的百分比随机布点作为初始结构，并以颗粒的两点概率函数作为重建约束条件，预测与其对应的三维颗粒结构，然后将重建的颗粒结构置于三维空间中。所述模拟退火算法在预测颗粒的三维颗粒结构时，以预测三维颗粒结构的两点概率函数与原始二维矿物颗粒的两点概率函数之间的误差来表征能量，其表达式如下：

(P₁₁)_orig代表原始二维矿物颗粒结构的两点概率函数，(P₁₁)_rec代表预测三维颗粒结构的两点概率函数，并以它们之间的误差作为重建过程中的能量E。

上式所表示的能量E是度量预测颗粒的三维颗粒结构与原始二维矿物颗粒结构统计的相似性，预测颗粒的统计分布特征与二维矿物颗粒结构的统计分布特征越接近，其能量值E越小，表征预测结构与原始颗粒越相似。将重建完成的单个颗粒放置在岩心的三维空间中，在重建过程中，应用周期性边界约束条件，周期性边界条件(Periodic BoundaryConditions,PBC)是边界条件的一种，反映的是如何利用边界条件替代所选部分系统受到周边环境的影响。可以看作是如果去掉周边环境，保持该系统不变应该附加的条件，也可以看作是由部分的性质来推广表达全局的性质。此约束条件重建的单颗颗粒三维颗粒结构可能被三维空间分割为几个颗粒块，并且在重建结构中会产生少量的噪点，此时需按周期性边界条件将重建结构中的几个颗粒块重新恢复成一个完整的颗粒，即将处于二维切面结构两侧的被分割开的颗粒拼接为一个颗粒后重新放入二维切面结构中并去除其中的噪点，然后将颗粒结构置于重建的三维空间中。图2-1到图2-4给出了三维颗粒重建过程中的二维切面示意图的变化过程：图2-1为矿物颗粒参考结构；图2-2为随机产生的三维颗粒结构的二维切面图；图2-3为重建结构的二维切面图；图2-4为重建结构进行去噪和复原的二维切面结构。完成一个颗粒的预测之后，再选取一个矿物颗粒结构重复上述颗粒重建过程，直到满足该阶段颗粒重建的终止约束条件为止。由于岩心二维颗粒图像中颗粒所占的百分比通常比较高，模拟退火算法预测的矿物结构通常不能完全放置在岩心的三维空间中。因此，在模拟退火算法预测三维结构基础上进一步进行第二阶段预测三维颗粒结构。

(6)所述第二阶段采用颗粒生长算法预测矿物颗粒三维结构：在步骤(5)第一阶段预测重建所得每个岩心二维矿物颗粒三维颗粒结构基础上进行第二阶段预测，即在余下的未被矿物颗粒占据的空间基础上，应用颗粒生长算法模拟岩心颗粒的生长过程，整个生长过程直到满足预测终止条件时结束，即完成岩心二维颗粒图像相对应的三维颗粒结构重建。

在第一阶段模拟退火算法预测结构基础上，即在余下的未被矿物颗粒占据的空间基础上，模拟岩心颗粒的生长过程，颗粒生长算法终止条件仍为：直方图匹配和拥挤度。在余下的三维空间中，随机选取一个没有被矿物颗粒占据的空间点作为矿物颗粒生长核，然后以该生长核为中心向其18邻域生长，为了得到各相同性的矿物颗粒结构，设置生长核的最近6邻域和12次近邻的生长概率的比值为4:1。生长概率的设置是基于矿物颗粒结构是一种聚合物的形状，因而生长核的最近6邻域的生长概率高于12次近邻的生长概率，因此，本实施例中将生长核的最近6邻域和12次近邻的生长概率的比值设置为4:1。生长核生长后的结构作为新的生长核，重复上述过程，直到颗粒生长的体积与预测颗粒结构的体积相同为止，即完成该颗粒的预测。在颗粒生长过程中如果生长的颗粒和其他颗粒结构相接触则停止该方向上颗粒的生长；在完成一个颗粒的生长之后，随机再选取一个没有被矿物颗粒占据的空间点作为新的颗粒生长核，重复上述过程，直到满足颗粒重建的终止约束条件，即完成整个矿物颗粒的三维重建。图3-1到图3-4给出了一个三维颗粒生长过程中的二维切面示意图的变化过程：图3-1为矿物颗粒参考结构；图3-2为随机选取的颗粒生长核；图3-3为颗粒生长核在二维切面中的八个生长方向；图3-4为该切面中颗粒生长的最终结果。

(7)完成上述两个阶段的重建即完成岩心二维矿物颗粒的三维颗粒结构重建，如图4所示，其中不同的灰度级表示不同的矿物颗粒。

所述岩心矿物颗粒与岩心孔隙相比，它们的结构特征不一样，矿物颗粒具有自己的结构特征。基于岩心二维颗粒图像的重建，用与矿物颗粒特征相一致的参数来描述重建结果。所述三维颗粒的主要参数包括面积、周长、粒径、形状因子、体积等。而对整个颗粒结构来说，则用这些参数的统计平均来表征。本实施例中，用岩心颗粒的粒径分布、形状因子、岩心定名等来说明本发明的有效性。本实施例对重建结构和原始结构中矿物颗粒的粒径分布和依据三维矿物颗粒粒径分布给出的岩心粒度定名及矿物颗粒的平均形状因子的分析来说明所提出的重建算法的有效性。从重建的三维颗粒结构中随机选取三个相互垂直的切面，比较这三个切面与参考图像中颗粒的结构、分布特性和孔隙度。图5-1、图5-2和图5-3分别为从重建结构中选取的切面图，并对其中的颗粒赋为不同的灰度级，将它们定名为切面1、切面2、切面3，切面中不同的灰度级代表不同的颗粒，图6-1、图6-2和图6-3分别为这三个切面图对应的孔隙分布。从三幅图5和三幅图6中可以看出，重建结构切面图中颗粒的分布与参考图像中颗粒的分布类似，具体的颗粒性能参数如下面的表1和表2所示。

表1 切面1矿物颗粒粒度分类

表2 切面2和切面3矿物颗粒粒度分类

从表1可知，给出了原始结构图1和切面1中矿物颗粒的粒径分布。表2为切面2和切面3的矿物颗粒粒径分布。从表1和表2中矿物颗粒的粒度分布可以得出原始结构图1和切面1、切面2、切面3的岩心粒度定名都为含粗粉砂极细砂岩。原始结构图1中颗粒的平均形状因子和孔隙度分别为0.57和0.063；重建结构图中切面1、切面2和切面3中颗粒平均形状因子分别为0.62、0.62和0.60，孔隙度分别为0.062、0.047、0.063。

颗粒的形状因子表征颗粒结构与圆的接近程度，用下式来计算：

F＝4πA_p/L² (2)

F表示颗粒的性质因子，A_p为颗粒的面积，L为颗粒的周长，对于岩心二维颗粒图像，一般用颗粒的平均形状因子来表示其结构形状。

颗粒的孔隙度为岩样中所有孔隙空间体积之和与该岩样体积的比值，以百分数表示。

从对重建结构切面图与原始图像的颗粒粒径分布、岩石粒度定名、平均形状因子及孔隙度的分析可以看出重建结构不同切面与原始图像中颗粒和孔隙的统计参数比较接近，从而反映了三维重建结构不同切面都捕获了参考图像中矿物颗粒主要的分布和结构特征，但各切面与参考图像中颗粒的分布和结构特征之间仍然存在一定的偏差，这也反映了重建三维颗粒结构不同切面之间颗粒的变化性。即说明本发明的岩心矿物颗粒三维重建算法的有效性。

本实施例中，提取岩心二维图像的颗粒结构，并统计每个矿物颗粒的面积，按颗粒面积大小给出二维矿物颗粒面积分布直方图，对给定的每个直方图区间的二维矿物颗粒，预测与该区间颗粒相对应的三维颗粒结构，并将预测的三维颗粒结构放入重建的三维空间中。通过比较重建三维颗粒结构不同切面与二维颗粒图像中颗粒的粒径分布、形状因子、岩心定名等特性参数，从而验证了本发明提出的基于岩心二维颗粒图像的三维颗粒结构重建方法的正确性。

上述实施例只是本发明的优选实施例，并不是对本发明所述技术方案的限定，只要是不经过创造性劳动即可在上述实施例的基础上实现的技术方案，均应视为落入本发明内容的保护范围内。

Claims (7)

1.一种基于岩心二维颗粒图像的三维颗粒结构重建方法，其特征在于包括以下步骤：

(1)提取原始岩心二维矿物颗粒图像中矿物颗粒的个数，对不同矿物颗粒赋予不同的灰度级，将处于不同灰度级的矿物颗粒进行标号，以区分出不同的颗粒；

(2)统计步骤(1)提取的原始岩心二维矿物颗粒图像中每个矿物颗粒的面积；采用遍历整幅岩心二维矿物颗粒图像，用计算机记录每个矿物颗粒标号下的二维颗粒像素值的个数，而每个标号下二维矿物颗粒的像素值个数代表其面积，并统计出岩心二维颗粒的面积；

(3)将步骤(2)统计出的每个矿物颗粒的面积，按其大小给出二维矿物颗粒面积分布直方图；

(4)对步骤(3)给出的每个二维矿物颗粒面积分布直方图区间的二维矿物颗粒，采用第一阶段模拟退火算法来预测其直方图区间中每个岩心二维矿物颗粒的三维结构；

(5)将步骤(4)预测的每个岩心二维矿物颗粒的三维颗粒结构放入重建的三维空间中；整个重建过程当满足预测终止条件时结束，则能获得每个岩心二维矿物颗粒的三维颗粒结构；

(6)对步骤(5)在第一阶段预测重建所得每个岩心二维矿物颗粒的三维结构的基础上采用第二阶段颗粒生长算法预测重建，在余下的未被颗粒占据的空间中，应用颗粒生长算法模拟岩心颗粒的生长过程，整个生长过程直到满足预测终止条件时结束，至此完成岩心二维颗粒图像的三维颗粒结构重建。

2.根据权利要求1所述的方法，其特征在于步骤(3)中所述按其大小给出二维矿物颗粒面积分布直方图，找出其中的最大面积和最小面积，用最大面积与最小面积的差值除以10，并将二维矿物颗粒面积分成10等份，计算出二维颗粒面积分布直方图每个区间的大小，再将每个颗粒按其面积大小存放到不同的各区间中。

3.根据权利要求1所述的方法，其特征在于步骤(4)中所述第一阶段模拟退火算法预测直方图区间中每个岩心二维颗粒的三维颗粒结构时，首先提取该岩心二维颗粒，然后在三维空间中以该岩心二维颗粒所占的体积分数随机布点作为初始结构，并以岩心二维颗粒的两点概率函数作为重建约束条件，预测三维颗粒结构。

4.根据权利要求1所述的方法，其特征在于步骤(5)中所述将预测的岩心二维矿物颗粒的三维颗粒结构放入重建的三维空间中，是在重建每一个三维颗粒后，在存储颗粒的三维空间中随机寻找没有被颗粒占据且能够完整存放三维颗粒的空间点放置重建的三维颗粒；当重建的三维颗粒结构的直方图与二维颗粒图像相匹配或重建三维颗粒结构占整个三维空间的百分比与二维颗粒图像相同时，重建达到终止条件。

5.根据权利要求1所述的方法，其特征在于步骤(6)中所述的第二阶段颗粒生长算法为：在步骤(5)重建结束后，在余下的三维空间中，随机选取一个没有被矿物颗粒占据的空间点作为矿物颗粒生长核，然后以该生长核为中心向其18邻域生长；该生长核生长后的结构作为新的矿物颗粒生长核，重复上述过程，直到重建三维颗粒结构的直方图与二维颗粒图像相匹配，或重建三维颗粒结构占整个三维空间的百分比与二维颗粒图像相同时，重建结束。

6.根据权利要求1或5所述的方法，其特征在于所述以生长核为中心向其18邻域生长，它的最近邻域为6邻域和12次近邻，其生长概率的比值为4:1，生长核生长后的结构作为新的生长核，重复上述过程，直到颗粒生长的体积与预测的颗粒结构的体积相同为止，完成所述三维颗粒结构的预测。

7.根据权利要求1所述的方法，其特征在于步骤(6)中所述的颗粒生长算法，是在颗粒生长过程中，如果生长的颗粒和其他颗粒结构相接触则停止该方向上颗粒的生长，完成一个颗粒的生长之后，随机再选取一个没有被矿物颗粒占据的空间点作为新的颗粒生长核，重复上述过程，直到颗粒生长的体积与预测的颗粒结构的体积相同为止，完成所述三维颗粒结构的预测。