CN112364504B - 基于CT扫描技术和Gabriel图的非均质多孔介质模型构建方法 - Google Patents

Abstract

本发明涉及基于CT扫描技术和Gabriel图的非均质多孔介质模型构建方法，包括步骤如下：(1)采用CT扫描技术获取孔隙结构以及粗糙壁面信息；(2)绘制Delaunay三角形，更新获得Gabriel图；(3)生成一组孔隙宽度；(4)判断Gabriel图中是否含有无效网格边；(5)基于每条网格边的位置作为孔径分布，拓展孔隙宽度，遍布所有像素点；(6)对多孔介质模型中的光滑孔隙壁面进行粗糙化。本发明利用CT扫描技术获取真实岩心内孔隙壁面的粗糙度分布，利用Gabriel图构造的孔隙网络分布具有连通性差以及可以随机出现盲道的特点，与真实地层的更加贴近。

Description

基于CT扫描技术和Gabriel图的非均质多孔介质模型构建 方法

技术领域

本发明涉及基于CT扫描技术和Gabriel图的非均质多孔介质模型构建方法，属于多孔介质模拟流动机理研究领域。

背景技术

多孔介质在自然界中广泛存在，如常见的岩石储层、大坝坝体等。当前，可以通过构建合理的多孔介质模型来认识多孔介质内流体运动的规律。当前，多孔介质内流体流动研究主要基于两类多孔介质模型，包括人造构造多孔介质模型以及天然多孔介质模型，其中，人为构造多孔介质模型由于具有可以独立改变其孔隙度、孔隙尺寸、配位数、均质性等物性参数的优点，因而在实验及数值模拟研究中得到了广泛的应用。

近年来随着显微技术以及微纳尺度实验的发展，研究人员发现孔隙壁面粗糙性对于流体输运、热传递以及物质扩散等物理化学过程存在显著的影响。传统人为构造的多孔介质模型主要有小球堆积模型、毛管束模型等等，均未考虑壁面粗糙性对于流体流动的影响。而真实多孔介质如地层岩石其孔隙壁面均是无规则粗糙表面，人为构造的多孔介质模型与其存在较大的差异。

现有孔隙网络分布的生成方法中，如Delaunay三角剖分，形成的孔隙网络分布均质性与连通性较好且未能考虑孔隙连通性中“盲道”的存在。通常在真实的多孔介质中许多孔隙之间非均质性强，连通性差，且伴随“盲道”的存在，上述剖分方法构造的孔隙网络分布与真实的情况不符。

发明内容

针对现有技术的不足，本发明提供了基于CT扫描技术和Gabriel图的非均质多孔介质模型构建方法。本发明利用CT扫描技术获取真实岩心内孔隙壁面的粗糙度分布，利用Gabriel图构造的孔隙网络分布具有连通性差以及可以随机出现盲道的特点，与真实地层的更加贴近。

本发明的技术方案为：

基于CT扫描技术和Gabriel图的非均质多孔介质模型构建方法，包括步骤如下：

步骤(1)，采用CT扫描技术扫描岩石样品，获取孔隙结构信息以及粗糙壁面信息，所述孔隙结构信息包括多孔介质孔隙宽度，所述粗糙壁面信息包括多孔介质孔隙的粗糙度；

步骤(2)，绘制Delaunay三角形，更新获得Gabriel图，并保存Gabriel图中网格边编号；

步骤(3)，基于步骤(1)CT扫描所得的孔隙结构信息，分析获得孔隙宽度的概率密度函数，利用孔隙宽度的概率密度函数生成一组孔隙宽度W₁、W₂、W₃…W_k，W_k表示标号为k的网格边对应的孔隙宽度；

步骤(4)，判断步骤(2)中所得的Gabriel图中是否含有无效网格边；

如果存在无效网格边，则将无效网格边去除，并取无效网格边的中点作为新的网格节点，然后连接新的网格节点，得到修正后的Gabriel图；再重新进行网格边编号，并保存新的网格边编号；

如果不存在无效网格边，则无需修正Gabriel图，保存步骤(2)得到的网格边编号；

步骤(5)，基于步骤(4)所得的Gabriel图中的每条网格边的位置作为孔隙分布，遍布目标区域的所有像素点，得到多孔介质模型文件，根据多孔介质模型文件生成光滑孔隙壁面的多孔介质模型；

步骤(6)，基于步骤(1)CT扫描所得的粗糙壁面信息，获得粗糙分布的概率密度函数；并利用所得的粗糙分布的概率密度函数对步骤(5)所得的多孔介质模型中的光滑孔隙壁面进行粗糙化，获得多孔介质孔隙壁面粗糙数据文件，从而完成非均质多孔介质模型的构建。

根据本发明优选的，步骤(5)中，基于步骤(4)所得的Gabriel图中的每条网格边的位置作为孔隙分布，遍布目标区域的所有像素点，得到多孔介质模型文件，根据多孔介质模型文件生成光滑孔隙壁面的多孔介质模型；具体过程为：

基于步骤(4)所得的Gabriel图中的每条网格边的位置作为孔隙分布，遍布目标区域的所有像素点；距离编号k的网格边小于等于0.5W_k的像素点的色彩数值设置为0，其余部分像素点的色彩数值设置为1，当像素点的色彩数值为0时，表示像素点对应位置的几何结构为多孔介质孔隙空间；当像素点的色彩数值为1时，表示像素点对应位置的几何结构为基岩，形成多孔介质模型文件；所述多孔介质模型文件包括每个像素点的色彩数值和像素点对应位置的几何结构；根据多孔介质模型文件生成光滑孔隙壁面的多孔介质模型。

根据本发明优选的，步骤(6)中，所述多孔介质孔隙壁面粗糙数据文件包括粗糙化处理后的二维平面内每个像素点对应位置的几何结构及像素点对应的色彩数值；粗糙化处理后的二维平面内每个像素点位置包括多孔介质孔隙空间和基岩，每个像素点对应的色彩数值包括0和1；当像素点的色彩数值为0时，表示像素点对应位置的几何结构为多孔介质孔隙空间；当像素点的色彩数值为1时，表示像素点对应位置的几何结构为基岩。

根据本发明的优选的，步骤(6)中，利用所得的粗糙分布的概率密度函数对步骤(5)所得的多孔介质模型中的光滑孔隙壁面进行粗糙化，具体步骤为：

a.将步骤(5)所得的多孔介质模型中光滑孔隙壁面所在位置作为基准面；

b.利用所得的粗糙分布的概率密度函数随机产生一组偏离基准面的距离ΔH₁、ΔH₂……ΔH_n；该组偏离基准面的距离服从步骤(1)中分析获得的粗糙分布的概率密度函数；

c.对所获取的偏离基准面距离ΔH进行判断，并进行不同的处理，具体为：

如果偏离基准面的距离ΔH等于0，则不对光滑孔隙壁面做任何处理；

如果偏离基准面的距离ΔH大于0，则对光滑的孔隙壁面进行凹陷处理，凹陷方向朝向多孔介质孔隙空间的外部，处理后粗糙壁面上的像素点距离基准面的距离为ΔH_n，并且将凹陷壁面垂直至基准面的所有像素点的色彩数值设置为0，表示凹陷壁面垂直至基准面的所有像素点对应的位置为孔隙空间；

如果偏离基准面的距离ΔH小于0，则对光滑孔隙壁面也进行凸起处理，凸起方向朝向多孔介质孔隙空间内部，处理后粗糙壁面上的像素点距离基准面的距离为ΔH，并且将凸起壁面垂直至基准面的所有像素点的色彩数值设置为1，表示凹陷壁面垂直至基准面的所有像素点对应的位置为基岩。

经粗糙化处理后的多孔介质孔隙，充分考虑了粗糙度带来的影响，相比较于传统的多孔介质模型无法考虑粗糙度是一大优势。并且粗糙度的分布取样自真实的岩心，从而更加贴近现场真实情况。

根据本发明优选的，所述步骤(4)中，无效网格边的判定方法如下：

根据步骤(3)生成的一组孔隙宽度W₁、W₂、W₃…W_k，取最大宽度W_max作为判断标准，

如果三角形的网格边的长度小于等于W_max，则判定该三角形的网格边为无效网格边；

如果三角形的网格边大于W_max，则判定该三角形的网格边为有效网格边，并进行编号保存。

去除无效网格边的多孔介质模型与真实地层情况更加贴近，符合真实的孔隙分布。

根据本发明优选的，步骤(3)中，基于步骤(1)CT扫描所得的孔隙结构信息，分析获得孔隙宽度的概率密度函数的过程：

步骤(3-1)，通过CT扫描获取的孔隙宽度数据，对获得数据进行统计；

步骤(3-2)，统计一共N个孔隙宽度数据，其中孔隙宽度为w1的有n1个，则出现w1的概率为n1/N；依次类推，计算不同孔隙宽度对应的概率；

步骤(3-3)，将所有孔隙宽度作为X轴数据，对应的概率为Y轴数据，进行曲线拟合得到孔隙宽度的概率密度函数；

进一步优选的，若存在震荡数据，则采用分段函数进行表示。

根据本发明的优选的，步骤(6)中，基于步骤(1)CT扫描所得的粗糙壁面信息，获得粗糙分布的概率密度函数，具体过程包括：

步骤(6-1)，通过步骤(1)CT扫描获取的粗糙壁面信息数据，对获得数据进行统计；

步骤(6-2)，统计一共S个壁面粗糙度数据，其中粗糙度为h1的有a1个，则出现h1的概率为a1/S；依次类推，计算不同壁面粗糙度对应的概率；

步骤(6-3)，将所有粗糙度作为X轴数据，对应的概率为Y轴数据，进行曲线拟合得到孔隙宽度的概率密度函数。

进一步优选的，步骤(6)中，若存在震荡数据，则采用分段函数进行表示。

根据本发明优选的，步骤(2)中，绘制Delaunay三角形，更新获得Gabriel图，并保存Gabriel图中网格边编号；具体过程包括：

步骤(2-1)，设定一个长为L_a、宽为L_b的二维平面，以二维平面的左下角顶点为原点建立直角坐标系，并在二维平面内生成A*B个像素点，A、B为正整数，所有像素点均匀遍布整个二维平面，获取所有像素点数据；所述像素点数据包括各个像素点在直角坐标系中的坐标值；随机选取n个像素点作为离散点，n为正整数，获得离散点数据；所述离散点数据包括各个离散点在直角坐标系中的坐标值；

步骤(2-2)，基于步骤2-1所得的各个离散点在直角坐标系中的坐标值，绘制Delaunay三角形；对Delaunay三角形的每条网格边进行网格边编号L1、L2、L3…Lj；编号完成后进行分组，以每一个三角形的三边为一组，并进行分组编号T1、T2、T3…Ti；并获得目标区域Delaunay三角网格信息，所述Delaunay三角网格信息包括每条网格边的长度以及每个网格节点的坐标信息；网格节点即组成Delaunay三角形的顶点；

步骤(2-3)，基于步骤2-2所得的网格边编号、分组编号以及目标区域Delaunay三角网格信息，绘制以每条网格边为直径的圆；

针对任一组三角形，判断以一个三角形的网格边为直径的圆内是否含有该组三角形的第三点，若含有第三点，则去除这个三角形的网格边；

若不含有第三点，则保留这个三角形的网格边；

完成所有判断后，去除所有的圆形，更新获得Gabriel图；并保存Gabriel图中网格边编号。

本发明的有益效果为：

1、本发明给出的基于CT扫描技术和Gabriel图的非均质多孔介质模型构建方法，利用CT扫描技术获取真实岩心内孔隙壁面的粗糙度分布，用于构造多孔介质模型孔隙空间的粗糙壁面及确定孔隙半径；利用Gabriel图确定孔隙网络分布，可以避免传统构造的多孔介质孔隙结构的分布过于均匀，周期性过强，连通性过好的问题，并且基于Gabriel图构造的孔隙网络分布具有连通性差以及可以随机出现盲道的特点，形成的多孔介质的模型与真实地层的更加贴近。

2、本发明通过CT扫描依据真实的岩石内部的孔隙结构，考虑了真实的粗糙结构，孔隙形状，孔隙大小等因素，与以往的模型相比更加贴近真实的多孔介质。

附图说明

图1为Gabriel图的定义示意图；

图2为发明基于CT扫描技术和Gabriel图的非均质多孔介质模型构建方法的流程图；

图3为实施例步骤(2)在目标区域随机给定离散点示意图；

图4为实施例步骤(3)根据离散点形成Delaunay三角形网格示意图；

图5为实施例步骤(4)更新后形成的Gabriel图示意图；

图6为实施例步骤(6)去除无效网格边，修正形成Gabriel图的过程示意图；

图7为实施例步骤(7)将步骤(6)Gabriel图的网格边扩展为孔隙示意图；

图8为实施例步骤(8)将步骤(7)的光滑壁面粗糙化示意图；

图9为对比例1中提供的传统多孔介质模型示意图。

具体实施方式

下面结合说明书附图和实施例对本发明进一步限定，但不限于此。

实施例1

基于CT扫描技术和Gabriel图的非均质多孔介质模型构建方法，如图2所示，包括步骤如下：

步骤(1)，采用CT扫描技术扫描岩石样品，获取孔隙结构信息以及粗糙壁面信息，所述孔隙结构信息包括多孔介质孔隙孔径大小，所述粗糙壁面信息包括多孔介质孔隙的粗糙分布信息；

CT扫描技术即电子计算机断层扫描，它是利用精确准直的X线束、γ射线、超声波等，具有扫描时间快，图像清晰等特点。

步骤(2)，设定一个长为1cm、宽为1cm的二维平面，并以二维平面的左下角顶点为原点建立直角坐标系，并在设定的二维平面内生成1000*1000个像素点，所有像素点均匀遍布整个二维平面，获取所有像素点数据，随机选取54个像素点作为离散点，离散点如图3所示,获得离散点数据，所述像素点数据包括各个像素点在直角坐标系中的坐标值；所述离散点数据包括各个离散点在直角坐标系中的坐标值；

步骤(3)，基于步骤(2)所得的各个离散点在直角坐标系中的坐标值，绘制Delaunay三角形；如图4所示，对Delaunay三角形的每条网格边进行网格边编号L1、L2、L3…Lj；编号完成后进行分组，以每一个三角形的三边为一组，并进行分组编号T1、T2、T3…Ti；并获得目标区域Delaunay三角网格信息，所述Delaunay三角网格信息包括每条网格边的长度以及每个网格节点的坐标信息；网格节点即组成Delaunay三角形的顶点；

步骤(4)，基于步骤(3)所得的网格边编号、分组编号以及目标区域Delaunay三角网格信息，绘制以每条网格边为直径的圆；

针对任一组三角形，判断以一个三角形的网格边为直径的圆内是否含有该组三角形的第三点，若含有第三点，则去除这个三角形的网格边；

若不含有第三点，则保留这个三角形的网格边；

完成所有判断后，去除所有的圆形，更新获得Gabriel图；并保存Gabriel图中网格边编号；更新后形成的Gabriel图如图5所示。

Gabriel Graph(简称GG)是一种拓扑模型，定义为存在点集E,u、v、w属于E,如果节点u和节点v之间，直径为uv的圆内，不存在其它顶点w，则节点u和节点v为Gabriel图的边(u，v)。图1形象地说明了Grabriel图的定义，节点u和节点v为直径的圆，节点u和节点v都在圆上，即是圆形区域内。若(u，v)是GG中的边，则在节点u和v之间的圆形区域内，不能包含有任何证明节点w，所以图1的uv不是GG边。

步骤(5)，基于步骤(1)CT扫描所得的孔隙结构信息，分析获得孔隙宽度的概率密度函数，利用孔隙宽度的概率密度函数生成一组孔隙宽度W₁、W₂、W₃…W_k；

步骤(5)中，孔隙宽度为随机设定并服从分析步骤(1)中CT扫描所得的孔隙结构信息后获得孔隙宽度的概率密度函数。具体过程为：

步骤(5-1)，通过CT扫描获取的孔隙宽度数据，对获得数据进行统计；

步骤(5-2)，统计一共N个孔隙宽度数据，其中孔隙宽度为w1的有n1个，则出现w1的概率为n1/N；依次类推，计算不同孔隙宽度对应的概率；

步骤(5-3)，将所有孔隙宽度作为X轴数据，对应的概率为Y轴数据，进行曲线拟合得到孔隙宽度的概率密度函数；若存在震荡数据，则采用分段函数进行表示。

具体的,CT扫描后获得一组孔径大小数据，如W₁＝1μm，W₂＝1.2μm，W₃＝1.2μm，W₄＝1μm，W₅＝1.2μm，W6＝1.25μm，W₇＝1.25μm，W₈＝1.25μm，W₉＝1.2μm，W₁₀＝1μm。

其中孔径大小为1μm的有3个，1.2μm的有4个，1.25μm的有3个。所以生成1μm的概率为0.3，生成1.2μm的概率为0.4，生成1.25μm的概率为0.3。那么如果需要生成100个孔径大小，则服从上述的概率分布。即1μm有30个，1.2μm有40个，1.25μm有30个。

上述数据较少近似可以看成是一种均匀分布比较简单。一旦有大量的数据，统计出来的结果会更加准确，更加贴近事实，可以使用同样的方法进行统计。获取在扫描孔径最小值到最大值之间所有孔隙大小的概率。

本实施例在扫描岩心后获得的密度函数如下：

孔隙宽度的概率密度函数f(W)，如式(I)所示：

式(I)中，W为孔隙宽度，W就是孔径的宽度，f(W)表示出现W宽度的概率值。

步骤(6)，判断步骤(4)中所得的Gabriel图中是否含有无效网格边；

如果存在无效网格边，则将无效网格边去除，并取无效网格边的中点作为新的网格节点，然后连接新的网格节点，得到修正后的Gabriel图；最后重新进行网格边编号，并保存新的网格边编号；

如果不存在无效网格边，则无需修正Gabriel图，保存步骤(4)得到的网格边编号信息；

所述步骤(6)中，无效网格边的判定方法如下：

根据步骤(5)生成的一组孔隙宽度W₁、W₂、W₃…W_k，取最大宽度W_max作为判断标准，

如果三角形的网格边的长度小于等于W_max，则判定该三角形的网格边为无效网格边；

如果三角形的网格边大于W_max，则判定该三角形的网格边为有效网格边，并进行编号保存。

去除无效网格边的多孔介质模型与真实地层情况更加贴近，符合真实的孔隙分布。

如图6所示，图6的左图为去除无效网格边前的示意图，左图黑色圈出部分为无效网格边；图6的右图为去除无效网格边后的示意图。

步骤(7)，基于步骤(6)所得的Gabriel图中的每条网格边的位置作为孔径分布，遍布目标区域的所有像素点；距离编号k的网格边小于等于0.5W_k的像素点的色彩数值设置为0，其余部分像素点的色彩数值设置为1，当像素点的色彩数值为0时，表示像素点对应位置的几何结构为多孔介质孔隙空间；当像素点的色彩数值为1时，表示像素点对应位置的几何结构为基岩，形成多孔介质模型文件，所述多孔介质模型文件包括每个像素点的色彩数值和像素点对应位置的几何结构；根据多孔介质模型文件生成光滑孔隙壁面的多孔介质模型；

对图7的上图中黑色圈出的区域的网格边扩展为孔隙，当k取值为a，b，c，d，e时，Gabriel图的网格边La、Lb、Lc、Ld、Le扩展为孔隙空间的示意图如图7的下图中所示，此步骤下的孔隙壁面是光滑的，这与真实的情况是不一致的。

步骤(8)，基于步骤(1)CT扫描所得的粗糙壁面信息，获得粗糙分布的概率密度函数；并利用所得的粗糙分布的概率密度函数对步骤(7)所得的多孔介质模型中的光滑孔隙壁面进行粗糙化，获得多孔介质孔隙壁面粗糙数据文件，所述多孔介质孔隙壁面粗糙数据文件包括粗糙化处理后的二维平面内每个像素点对应位置的几何结构及像素点对应的色彩数值；粗糙化处理后的二维平面内每个像素点位置包括多孔介质孔隙空间和基岩，每个像素点对应的色彩数值包括0和1；当像素点的色彩数值为0时，表示像素点对应位置的几何结构为多孔介质孔隙空间；当像素点的色彩数值为1时，表示像素点对应位置的几何结构为基岩；从而完成非均质多孔介质模型的构建。

步骤(8)中，利用所得的粗糙分布的概率密度函数对步骤(7)所得的多孔介质模型中的光滑孔隙壁面进行粗糙化，具体步骤为：

a.将步骤(7)所得的多孔介质模型中光滑孔隙壁面所在位置作为基准面；

b.利用所得的粗糙分布的概率密度函数随机产生一组偏离基准面的距离ΔH₁、ΔH₂……ΔH_n；该组偏离基准面的距离服从步骤(1)中分析获得的粗糙分布的概率密度函数，具体步骤为：

步骤(8-1)，通过步骤(1)CT扫描获取的粗糙壁面信息数据，对获得数据进行统计；

步骤(8-2)，统计一共S个壁面粗糙度数据，其中粗糙度为h1的有a1个，则出现h1的概率为a1/S；依次类推，计算不同壁面粗糙度对应的概率；

步骤(8-3)，将所有粗糙度作为X轴数据，对应的概率为Y轴数据，进行曲线拟合得到孔隙宽度的概率密度函数。步骤(8)中，若存在震荡数据，则采用分段函数进行表示。

本实施例下的粗糙分布的概率密度函数如下：

式(II)中ΔH为像素点距离基准面的距离，ΔH是偏离基准面的距离，f(ΔH)是出现ΔH的概率；

c.对所获取的偏离基准面距离ΔH进行判断，并进行不同的处理，如图8中的上图所示，具体为：

如果偏离基准面的距离ΔH等于0，则不对光滑孔隙壁面做任何处理；

如果偏离基准面的距离ΔH大于0，则对光滑的孔隙壁面进行凹陷处理，凹陷方向朝向多孔介质孔隙空间的外部，处理后粗糙壁面上的像素点距离基准面的距离为ΔH_n，并且将凹陷壁面垂直至基准面的所有像素点的色彩数值设置为0，表示凹陷壁面垂直至基准面的所有像素点对应的位置为孔隙空间；

如果偏离基准面的距离ΔH小于0，则对光滑孔隙壁面也进行凸起处理，凸起方向朝向多孔介质孔隙空间内部，处理后粗糙壁面上的像素点距离基准面的距离为ΔH，并且将凸起壁面垂直至基准面的所有像素点的色彩数值设置为1，表示凹陷壁面垂直至基准面的所有像素点对应的位置为基岩。图7下图中的光滑孔隙壁面经过粗糙化后，如图8中的下图所示。

经粗糙化处理后的多孔介质孔隙，充分考虑了粗糙度带来的影响，相比较于传统的多孔介质模型无法考虑粗糙度是一大优势。并且粗糙度的分布取样自真实的岩心，从而更加贴近现场真实情况。

本实施例提出了基于CT扫描技术和Gabriel图的粗糙壁面多孔介质模型构建方法。利用CT扫描技术获取真实岩心内孔隙壁面的粗糙度分布及孔径分布，用于构造多孔介质模型孔隙空间的粗糙壁面及确定孔隙半径，利用Gabriel图确定孔隙网络分布。克服了当前微尺度渗流机制研究中多孔介质模型几何结构的均质性过强，孔隙的连通结构与实际孔隙偏差较大，未考虑孔隙壁面的粗糙性等缺陷。将CT扫描获得的真实多孔介质的相关信息用于构建多孔介质模型，既可以保留人为构造多孔介质模型的优点，又能使所构建的多孔介质模型更加符合实际。

本专利采用CT扫描技术得到图像，保存图像信息，基于统计规律，形成对应的概率密度函数，最后在得到的Gabriel图上赋予相应的值。其次采用Gabriel图剖分的孔隙分布会存在一些不连通的结构，称之为“盲道”，Gabriel图中盲道存在更加符合真实的情况。

对比例1

现有技术中，使用人为方法构造二维多孔介质模型，如图9所示，该模型黑色部分为孔隙空间，白色部分为基岩。孔隙空间壁面是光滑的，没有考虑到真实孔隙中的粗糙度问题。与本实施例相比它的孔隙分布更加简单，不能很好的体现孔隙分布的复杂性。这样的模型在进行研究时会产生较大的误差。

与对比例1相比，本发明提出了一种基于CT扫描技术和Gabriel图的粗糙壁面多孔介质模型构建方法。本发明利用CT扫描技术获取真实岩心内孔隙壁面的粗糙度分布及孔径分布，用于构造多孔介质模型孔隙空间的粗糙壁面及确定孔隙半径，利用Gabriel图确定孔隙网络分布。克服了当前微尺度渗流机制研究中多孔介质模型几何结构的均质性过强，孔隙的连通结构与实际孔隙偏差较大，未考虑孔隙壁面的粗糙性等缺陷。

Claims (7)

1.基于CT扫描技术和Gabriel图的非均质多孔介质模型构建方法，其特征在于，包括步骤如下：

步骤(1)，采用CT扫描技术扫描岩石样品，获取孔隙结构信息以及粗糙壁面信息，所述孔隙结构信息包括多孔介质孔隙宽度，所述粗糙壁面信息包括多孔介质孔隙的粗糙度；

步骤(2)，绘制Delaunay三角形，更新获得Gabriel图，并保存Gabriel图中网格边编号；

步骤(3)，基于步骤(1)CT扫描所得的孔隙结构信息，分析获得孔隙宽度的概率密度函数，利用孔隙宽度的概率密度函数生成一组孔隙宽度W₁、W₂、W₃…W_k，W_k表示标号为k的网格边对应的孔隙宽度；

步骤(4)，判断步骤(2)中所得的Gabriel图中是否含有无效网格边；

如果存在无效网格边，则将无效网格边去除，并取无效网格边的中点作为新的网格节点，然后连接新的网格节点，得到修正后的Gabriel图；再重新进行网格边编号，并保存新的网格边编号；

如果不存在无效网格边，则无需修正Gabriel图，保存步骤(2)得到的网格边编号；

步骤(5)，基于步骤(4)所得的Gabriel图中的每条网格边的位置作为孔隙分布，遍布目标区域的所有像素点，得到多孔介质模型文件，根据多孔介质模型文件生成光滑孔隙壁面的多孔介质模型；

步骤(6)，基于步骤(1)CT扫描所得的粗糙壁面信息，获得粗糙分布的概率密度函数；并利用所得的粗糙分布的概率密度函数对步骤(5)所得的多孔介质模型中的光滑孔隙壁面进行粗糙化，获得多孔介质孔隙壁面粗糙数据文件，从而完成非均质多孔介质模型的构建；

其中，步骤(2)中，绘制Delaunay三角形，更新获得Gabriel图，并保存Gabriel图中网格边编号；具体过程包括：

步骤(2-1)，设定一个长为La、宽为Lb的二维平面，以二维平面的左下角顶点为原点建立直角坐标系，并在二维平面内生成A*B个像素点，A、B为正整数，所有像素点均匀遍布整个二维平面，获取所有像素点数据；所述像素点数据包括各个像素点在直角坐标系中的坐标值；随机选取n个像素点作为离散点，n为正整数，获得离散点数据；所述离散点数据包括各个离散点在直角坐标系中的坐标值；

步骤(2-2)，基于步骤(2-1)所得的各个离散点在直角坐标系中的坐标值，绘制Delaunay三角形；对Delaunay三角形的每条网格边进行网格边编号L1、L2、L3…Lj；编号完成后进行分组，以每一个三角形的三边为一组，并进行分组编号T1、T2、T3…Ti；并获得目标区域Delaunay三角网格信息，所述Delaunay三角网格信息包括每条网格边的长度以及每个网格节点的坐标信息；

步骤(2-3)，基于步骤(2-2)所得的网格边编号、分组编号以及目标区域Delaunay三角网格信息，绘制以每条网格边为直径的圆；

针对任一组三角形，判断以一个三角形的网格边为直径的圆内是否含有该组三角形的第三点，若含有第三点，则去除这个三角形的网格边；

若不含有第三点，则保留这个三角形的网格边；

完成所有判断后，去除所有的圆形，更新获得Gabriel图；并保存Gabriel图中网格边编号。

2.根据权利要求1所述的基于CT扫描技术和Gabriel图的非均质多孔介质模型构建方法，其特征在于，步骤(5)中，基于步骤(4)所得的Gabriel图中的每条网格边的位置作为孔隙分布，遍布目标区域的所有像素点，得到多孔介质模型文件，根据多孔介质模型文件生成光滑孔隙壁面的多孔介质模型；具体过程为：

基于步骤(4)所得的Gabriel图中的每条网格边的位置作为孔隙分布，遍布目标区域的所有像素点；距离编号k的网格边小于等于0.5W_k的像素点的色彩数值设置为0，其余部分像素点的色彩数值设置为1，当像素点的色彩数值为0时，表示像素点对应位置的几何结构为多孔介质孔隙空间；当像素点的色彩数值为1时，表示像素点对应位置的几何结构为基岩，形成多孔介质模型文件；所述多孔介质模型文件包括每个像素点的色彩数值和像素点对应位置的几何结构；根据多孔介质模型文件生成光滑孔隙壁面的多孔介质模型。

3.根据权利要求2所述的基于CT扫描技术和Gabriel图的非均质多孔介质模型构建方法，其特征在于，步骤(6)中，所述多孔介质孔隙壁面粗糙数据文件包括粗糙化处理后的二维平面内每个像素点对应位置的几何结构及像素点对应的色彩数值；粗糙化处理后的二维平面内每个像素点位置包括多孔介质孔隙空间和基岩，每个像素点对应的色彩数值包括0和1；当像素点的色彩数值为0时，表示像素点对应位置的几何结构为多孔介质孔隙空间；当像素点的色彩数值为1时，表示像素点对应位置的几何结构为基岩。

4.根据权利要求3所述的基于CT扫描技术和Gabriel图的非均质多孔介质模型构建方法，其特征在于，步骤(6)中，利用所得的粗糙分布的概率密度函数对步骤(5)所得的多孔介质模型中的光滑孔隙壁面进行粗糙化，具体步骤为：

a.将步骤(5)所得的多孔介质模型中光滑孔隙壁面所在位置作为基准面；

b.利用所得的粗糙分布的概率密度函数随机产生一组偏离基准面的距离ΔH₁、ΔH₂……ΔH_n；

c.对所获取的偏离基准面距离ΔH进行判断，并进行不同的处理，具体为：

如果偏离基准面的距离ΔH等于0，则不对光滑孔隙壁面做任何处理；

如果偏离基准面的距离ΔH大于0，则对光滑的孔隙壁面进行凹陷处理，凹陷方向朝向多孔介质孔隙空间的外部，处理后粗糙壁面上的像素点距离基准面的距离为ΔH_n，并且将凹陷壁面垂直至基准面的所有像素点的色彩数值设置为0，表示凹陷壁面垂直至基准面的所有像素点对应的位置为孔隙空间；

如果偏离基准面的距离ΔH小于0，则对光滑孔隙壁面也进行凸起处理，凸起方向朝向多孔介质孔隙空间内部，处理后粗糙壁面上的像素点距离基准面的距离为ΔH，并且将凸起壁面垂直至基准面的所有像素点的色彩数值设置为1，表示凹陷壁面垂直至基准面的所有像素点对应的位置为基岩。

5.根据权利要求1所述的基于CT扫描技术和Gabriel图的非均质多孔介质模型构建方法，其特征在于，所述步骤(4)中，无效网格边的判定方法如下：

根据步骤(3)生成的一组孔隙宽度W₁、W₂、W₃…W_k，取最大宽度W_max作为判断标准，

如果三角形的网格边的长度小于等于W_max，则判定该三角形的网格边为无效网格边；

如果三角形的网格边大于W_max，则判定该三角形的网格边为有效网格边，并进行编号保存。

6.根据权利要求1所述的基于CT扫描技术和Gabriel图的非均质多孔介质模型构建方法，其特征在于，步骤(3)中，基于步骤(1)CT扫描所得的孔隙结构信息，分析获得孔隙宽度的概率密度函数的过程：

步骤(3-1)，通过CT扫描获取的孔隙宽度数据，对获得数据进行统计；

步骤(3-2)，统计一共N个孔隙宽度数据，其中孔隙宽度为w1的有n1个，则出现w1的概率为n1/N；依次类推，计算不同孔隙宽度对应的概率；

步骤(3-3)，将所有孔隙宽度作为X轴数据，对应的概率为Y轴数据，进行曲线拟合得到孔隙宽度的概率密度函数。

7.根据权利要求1所述的基于CT扫描技术和Gabriel图的非均质多孔介质模型构建方法，其特征在于，步骤(6)中，基于步骤(1)CT扫描所得的粗糙壁面信息，获得粗糙分布的概率密度函数，具体过程包括：

步骤(6-1)，通过步骤(1)CT扫描获取的粗糙壁面信息数据，对获得数据进行统计；

步骤(6-2)，统计一共S个壁面粗糙度数据，其中粗糙度为h1的有a1个，则出现h1的概率为a1/S；依次类推，计算不同壁面粗糙度对应的概率；

步骤(6-3)，将所有粗糙度作为X轴数据，对应的概率为Y轴数据，进行曲线拟合得到孔隙宽度的概率密度函数。

Images