一种改进型多孔介质三维重建方法

技术领域

本发明属于计算机技术领域，尤其涉及一种改进型多孔介质三维重建方法。

背景技术

根据有限张多孔介质图像的统计信息，基于数学理论的随机重建法，是多孔介质的三维数字岩心重建最有吸引力的选择之一。模拟退火算法作为随机重建法的主流技术，因其可以根据材料特性采用不同类型的统计函数来重建多孔介质数字岩心而被广泛应用。传统的模拟退火重建中，重建的孔隙结构模型对比表明，它们在相关函数、几何结构和拓扑结构与参考模型具有高度的一致性，限制传统模拟退火算法(SAA)应用的主要障碍仍然是耗时问题。在传统SAA重建过程中，通过交换系统中的像素进行扰动更新，需要相当多的时间才能形成孔隙的基本结构。每一次像素交换后，相关函数统计值的更新也耗费了大量的时间。像素体的交换方式以及相关函数的统计值是导致SAA重建效率低下的主要因素。近年来，一些学者在这两个环节做出了重要的贡献，Tang、Zhao等采用不同相邻(DPNS)像素选择方案，降低重建时像素的冗余运算来提高像素交换的效率[6,7]；Ju等采用多线程的并行运算方案，充分调动计算机的性能，提高相关函数值的统计效率来改善SAA的重建效率[8]；Song等提出了基于增量的模拟退火算法[9]，计算像素交换前后的变化量作为贡献值，替代传统线性路径函数和两点概率函数的统计值，提高了相关函数值的更新效率，使得SAA重建大规模多孔介质的成为可能，可以说是目前SAA孔结构数字岩心重建最高效的算法之一。尽管如此，正如Song，Ju等人研究所指出的，单一环节的改善效果有限，重建大尺寸模型仍然存在困难，耗时长的问题依然需要被解决。

发明内容

本发明实施例的目的在于提供一种改进型多孔介质三维重建方法，旨在解决上述背景技术中提出的问题。

本发明实施例是这样实现的，一种改进型多孔介质三维重建方法，所述方法包括以下步骤：

包括以下步骤：

步骤1，选择需要重建的多孔介质材料二维切片图像，进行孔隙分割和提取，生成仅含有孔隙和基质两相是二值图像，并以此作为三维重建时的参考图像。

步骤2，基于参考图像中的孔隙像素数量，计算孔隙率。基于孔隙率，生成随机分布的孔隙与基质两相的多孔介质初始重建岩心。

步骤3，采用随机方向的抽样统计法进行相关函数(两点概率函数，线性路径函数)的统计计算。两点概率函数和线性路径函数的计算形式如下：

/>

其中，S(r)，L(r)分别表示两点概率函数和线性路径函数，N(I(x

统计时，随机选择重建系统中的孔隙像素点N个，并以这些点为起始点，产生[0,180°]范围内与水平方向(三维时为水平面)的夹角，随后在系统内的生成统计距离r的线段，线段两端的像素点都为孔隙相的数量，并计算频率值作为两点概率函数的统计值。线段上的所有像素点都为孔隙相的频率值作为线性路径函数的统计值。

分别统计参考图像和初始重建岩心的相关函数统计值，根据下列计算式，得到模拟退火系统的初始能量。

其中，E

步骤4，搭建元胞自动机的像素演化模块，来进行三维重建时的系统扰动。其中元胞自动机详细过程见下文具体实施步骤4.

步骤5，设定算法基础参数，运行步骤3像素交换模块，系统进入迭代环节，像素交换后的重建系统的能量将被重新计算，系统第k次扰动时的能量根据下式所确定：

重复进行系统的循环扰动，根据模拟退火判定准则以确保元胞自动机的演化是否符合重建要求，直至模拟退火三维重建系统达到能量降温条件，终止循环迭代，输出重建的三维模型。

本发明实施例提供的一种改进型多孔介质三维重建方法，其相比于现有技术，具有以下有益效果：重建效率相比于传统的重建算法有较大的改善。相关函数、孔隙的连通性、孔径分布与真实样品具有很好的一致性，重建质量优于传统的模拟退火三维重建算法。

附图说明

图1是一种改进型多孔介质三维重建方法的流程图；

图2是一种改进型多孔介质三维重建方法中对参考图像的预处理示意图；

图3是一种改进型多孔介质三维重建方法中随机生成的初始重建岩心；

图4是一种改进型多孔介质三维重建方法中两点概率函数和线性路径函数的统计示意图；

图5是一种改进型多孔介质三维重建方法中元胞自动机的像素邻域；

图6是一种改进型多孔介质三维重建方法中像素演化规则示意图；

图7是一种改进型多孔介质三维重建方法中元胞自动机模块示意图；

图8是一种改进型多孔介质三维重建方法中重建的数字岩心是重建的数字岩心示意图。

具体实施方式

为了使本发明的目的、技术方案及优点更加清楚明白，以下结合附图及实施例，对本发明进行进一步详细说明。应当理解，此处所描述的具体实施例仅仅用以解释本发明，并不用于限定本发明。

以下结合具体实施例对本发明的具体实现进行详细描述。

本发明提供的一种改进型多孔介质三维重建方法，解决了背景技术中的技术问题。

如图1所示，为本发明的一个实施例提供的一种改进型多孔介质三维重建方法的主流程图，所述一种改进型多孔介质三维重建方法包括：

包括5个步骤：①确定参考图像以及图像的预处理；②初始重建岩心生成；③相关函数的统计方法确定以及初始能量的计算；④元胞自动机像素扰动模型搭建；⑤重建系统扰动更新，温度冷却及重建结束。其具体流程如图1所示。

步骤1，参考图像的确定以及图像预处理，为了获取更准确的孔隙结构，对二维图像进行重建前的预处理。如图2所示，基于Micro-CT二维切片图像，采用Retinx图像增强技术，突出图像中的孔隙像素，利用阈值分割法获得三维重建的二值图片，以此作为三维重建时的参考系统。其中孔隙为白色，基质为黑色。

步骤2，生成初始重建岩心，根据参考系统中的孔隙率，随机生成孔隙与基质像素随机分布的三维数据体，作为初始重建系统。如图3所示，基质为灰色，孔隙为透明色。其中孔隙率比例17.73％

步骤3，两点概率函数和线性路径函数的计算形式如下：

其中，S(r)，L(r)分别表示两点概率函数和线性路径函数，N(I(x

采用随机方向的抽样统计法进行相关函数(两点概率函数，线性路径函数)的统计计算。具体如下：首先，随机选择系统中的体素点N个，以此作为统计线段的起始点。随后随机分配一个[0,180°]的角度，根据起始点和角度生成统计距离r的统计线段(如图4中S′

分别统计参考图像和初始重建岩心的相关函数统计值，根据下列计算式，得到模拟退火系统的初始能量。

/>

其中，E

步骤4.元胞自动机模型搭建，包括输入状态、输出状态、演化规则三个部分的简易自动机模型。将重建模型中的像素视为一系列的元胞个体，像素与其周围的空间组合形成一系列的离散子集(简称为元胞)。在对称性的情况下，元胞由中心和邻域组成。图1中给出了邻域半径为1的Von Neumann neighbourhood(VNN)和Moore neighbourhood(MN)两种常见的CA邻域形式(Maerivoet and De Moor,2005)。图5(a)是二维邻域，每一个格子代表二维图像的像素(孔隙或基质)，红色标记为有限自动机的中心。图5(b)是三维邻域，每个方块代表三维模型的体素。从图中可以看出，二维CA邻域中，VNN由邻域中心和四个相邻单元组成，MN由邻域中心和8个单元组成。三维的VNN有6个相邻单元，MN邻域有26个相邻单元。

确定交换像素的邻域范围后，元胞自动机在进行像素交换时，根据像素邻域的属性判别中心像素是否作为扰动时的像素。具体为：首先，随机选择重建系统中以孔隙和基质像素为中心的邻域，对邻域的属性信息进行统计，并记录当前的状态信息，作为输入状态，其定义为：

X

其中s

随后，定义中心像素的演化规则。设定一个允许中心像素进行演化的判定阈值，如邻域中心像素为孔隙时，若输入状态中基质像素数量大于设定的阈值时，则下一刻邻域中心像素演化为基质像素。反之，中心像素则保持不变。同理，若邻域中心像素为基质像素时，输入状态中的孔隙像素数大于设定的阈值时，下一刻邻域中心基质像素变化为孔隙。

自动机的演化规则为输入状态到输出状态的单值映射，表示为f:X

Y

完成上述元胞自动机的搭建后，系统在扰动时，同步进行中心像素的转化，即当某个元胞的中心体素由基质转化为孔隙，同时另一元胞的中心体素将由孔隙转化为基质。如图6所示，图6(a)的中心像素由基质转化为孔隙，同时如图6(b)所示中心像素由孔隙转化为基质。

完成上述元胞自动机模型的设定，如图7所示，搭建一个对重建系统进行扰动算法模块。算法可以归纳为三个部分，首先随机选择重建系统中孔隙和基质像素，以此为中心像素。随后记录相邻单元的像素组成，生成自动机的输入状态。最后进行像素演化判别，符合像素演化条件的进行像素演化。

步骤5，三维重建系统进入扰动更新阶段，当系统中的像素交换后，系统能量将被重新计算，系统第k次扰动后的能量计算如下式：

根据Metropolis准则判断新的重建模型是否被接受，判定准则定义如下：

/>

P(ΔE)≥rand[0,1) (8)

其中，

为了能够加载上述方法和系统能够顺利运行，该系统除了包括上述各种模块之外，还可以包括比上述描述更多或更少的部件，或者组合某些部件，或者不同的部件，例如可以包括输入输出设备、网络接入设备、总线、处理器和存储器等。

本应该理解的是，虽然本发明各实施例的流程图中的各个步骤按照箭头的指示依次显示，但是这些步骤并不是必然按照箭头指示的顺序依次执行。除非本文中有明确的说明，这些步骤的执行并没有严格的顺序限制，这些步骤可以以其它的顺序执行。而且，各实施例中的至少一部分步骤可以包括多个子步骤或者多个阶段，这些子步骤或者阶段并不必然是在同一时刻执行完成，而是可以在不同的时刻执行，这些子步骤或者阶段的执行顺序也不必然是依次进行，而是可以与其它步骤或者其它步骤的子步骤或者阶段的至少一部分轮流或者交替地执行。

以上所述实施例的各技术特征可以进行任意的组合，为使描述简洁，未对上述实施例中的各个技术特征所有可能的组合都进行描述，然而，只要这些技术特征的组合不存在矛盾，都应当认为是本说明书记载的范围。

以上所述实施例仅表达了本发明的几种实施方式，其描述较为具体和详细，但并不能因此而理解为对本发明专利范围的限制。应当指出的是，对于本领域的普通技术人员来说，在不脱离本发明构思的前提下，还可以做出若干变形和改进，这些都属于本发明的保护范围。因此，本发明专利的保护范围应以所附权利要求为准。

以上所述仅为本发明的较佳实施例而已，并不用以限制本发明，凡在本发明的精神和原则之内所作的任何修改、等同替换和改进等，均应包含在本发明的保护范围之内。