一种基于多重卷积计算的数字岩心微孔区体素网格粗化方法

1.本发明属于数字岩心分析技术，应用于电化学、地下油气开采、岩土工程等多方面，重点用于油气流体矿藏的数值模拟技术。


背景技术：

2.基于数字岩心的数值模拟作为一种成本低，周期短的研究手段，对油气储层渗流特性的确定和后续开采方案的优化具有重要意义。然而，对于页岩、致密砂岩和碳酸盐岩等这些具有多尺度孔隙结构的岩心，受扫描设备视域和分辨率相互制约的影响，很难在单一分辨率下提取到完整的孔隙结构：当采用较低的分辨率进行扫描时，所提取的孔隙结构会忽略低于图像分辨率的孔隙，影响模拟的准确性；当提高图像的分辨率来获取微孔结构时，数字岩心的真实尺寸远小于用像素
·
表征的体元大小，导致模拟的结果不具有代表性。基于数字岩心的多尺度数值模型是通过多阈值分割将三维扫描图像体素划分为固体区、大孔区和微孔区，生成多尺度数字岩心，并在大孔区和微孔区采用不同尺度的控制方程进行求解，以实现多尺度数字岩心的高精度数值模拟。所述多尺度是指有些岩石的孔隙分布呈现双峰的形状，其岩石孔隙的孔径有数个量级尺寸变化。
3.目前，针对多尺度孔隙结构数字岩心所开发的流体渗流数值模型主要有3类：双重孔隙网络模型、微观-连续介质模型和孔隙网络-达西模型。双重孔隙网络模型对大孔区和微孔区同时进行了简化，计算速度快，但模拟精度较低。微观-连续介质模型则基于原始图像，对大孔区和微孔区体素进行标记，通过求解darcy-brinkman-stoke方程和质量守恒方程来得到体素单元中的压力和速度分布，该方法模拟精度较高，但计算效率较低，所模拟的数字岩心尺寸通常小于3003体素；孔隙网络-达西模型结合了孔隙网络模型和微观-连续介质模型的优点，其将大孔区简化为孔隙网络，而微孔则视为连续介质，通过hagen
–
poiseuille方程和darcy方程分别控制流体在大孔区和微孔区的流动，在保证计算精度前提下提高了计算效率，但随着模拟岩心的增大，微孔区体素数量会成指数倍增加，当存在数千万至几亿微孔区体素时，与darcy方程相关的计算量会非常巨大，使得数值模拟无法进行下去。
4.综上可知，现有多尺度数字岩心模型在保证计算精度的前提下，受限于计算网格单元数量和计算效率，往往只能对较小的区域进行模拟研究。因此，如何在保证计算结果精度的前提下，极大地减少计算孔隙网格的数量，提高计算效率，是多尺度数字岩心数值模型所要解决的重要问题。


技术实现要素：

5.本发明的目的是为了解决现有数值模型中计算效率、计算精度与计算域大小相互制约的问题，提供一种基于多重卷积计算的数字岩心微孔区体素网格粗化方法，它能在确保大孔区与微孔区交界面处网格精度和考虑微孔区非均质性的前提下，实现从微孔区中心向外的多级体素网格粗化，极大地减少了计算网格单元数量，提高了计算效率。
6.本发明的技术方案是，一种基于多重卷积计算的数字岩心微孔区体素网格粗化方法，包括以下步骤：
7.步骤1、导入原始三维扫描图像，对图像进行滤波除噪声、多点阈值分割，获得大孔区和微孔区的分布信息；
8.步骤2、对步骤1获得的大孔区体素，利用图像处理中的分水岭算法对大孔区进行分割，基于分割图提取孔隙网络，并获取孔喉(指孔隙和喉道)的相关几何信息，包括体积、比表面积、内切圆半径；
9.步骤3、对步骤1获得的微孔区体素数据进行多重卷积计算，依次遍历卷积矩阵，对满足粗化条件的卷积值所映射的微孔区体素进行粗化标记；
10.步骤4、根据步骤3获得粗化标记的微孔区矩阵和步骤2获得的分割矩阵，生成微孔区粗化后的网格单元信息和连接关系；
11.步骤5、将大孔区的孔隙网络与微孔区的网格单元合并，生成数字岩心混合模型的计算网格文件。
12.进一步，还包括步骤6、使用步骤5生成的数字岩心混合模型的计算网格文件进行数值模拟计算。
13.本发明的技术效果是：
14.在保留微孔区几何分布特征和非均质性前提下，克服了现有数字岩心直接数值模拟方法计算网格单元多，难以开展全岩心数值模拟的缺点，在保证模拟结果准确性的前提下，极大地减少了微孔区的计算网格数量，显著提高计算效率。且本方法发明适用范围比较广泛，对于任何以图像体素为计算网格单元进行的数值模拟，均可用本发明实施体素网格的可控多级粗化。
附图说明
15.本发明的附图说明如下：
16.图1为本方法发明的步骤框图；
17.图2为实施例的岩心图像；
18.图3为实施例步骤1阈值分割提取的大孔区可视化图像；
19.图4为实施例步骤1阈值分割提取的微孔区可视化图像；
20.图5为实施例中大孔区提取的孔隙网络；
21.图6为未粗化和粗化的微孔区显示图；
22.图7为未粗化与粗化的混合模型计算网格单元数量对比图；
23.图8为未粗化与粗化的混合模型模拟结果对比图。
具体实施方式
24.下面结合附图和实施例对本发明作进一步说明：
25.如图1所示，本方法发明包括以下步骤：
26.步骤1、导入三维岩心图像(如微纳米ct扫描图像)进行滤波除噪、多阈值分割，获取固体区、大孔区和微孔区的分布信息；
27.多阈值分割与常规的阈值分割相同，只是使用了两个阈值来分别提取图像中的大
孔体素和微孔体素；具体为：
28.根据图像的ct值分布，设置阈值a、b，其中b＞a＞0，其与各相关系如下：
29.当像素ct值小于a时，其为大孔；
30.当像素ct值大于a小于b时，其为微孔；
31.当像素ct值大于b时，其为固体。
32.提取大孔区和微孔区信息，分别存储在两个与原始图像尺寸相同的体素矩阵中。在大孔区矩阵中大孔区体素值为1，其余(指微孔和固体)体素值为0；在微孔区矩阵中微孔区体素值为1，其它(指大孔和固体)体素值为0。
33.微孔区矩阵和大孔矩阵区用于提取微孔网格单元和大孔孔隙网络。
34.步骤2、对步骤1获得的大孔区体素，利用图像处理中的分水岭算法对大孔区进行分割，基于分割图提取孔隙网络，并获取孔喉的相关几何信息，包括体积、比表面积、内切圆半径。
35.有关图像分割的分水岭算法文献有：gostick,jefft，versatile and efficientpore network extraction method using marker-based watershed segmentation[j].physical review e,2017,96(2):023307(“基于标记分水岭分割的通用高效孔隙网络提取方法”，gostick,jefft，物理评论e，2017,96(2):023307)。
[0036]
步骤3、对步骤1获得的微孔区体素矩阵数据进行多重卷积计算，并微孔区体素进行粗化标记。
[0037]
步骤31、由步骤1获得的微孔区矩阵进行第一次卷积计算，卷积核大小为2
×2×
2，权重系数为1，步长为2，得到三维矩阵a1，然后，对a1矩阵做阈值分割提取满足粗化条件的区域，对阈值分割后的矩阵做腐蚀计算得到第一级粗化矩阵b1；
[0038]
步骤32、由步骤31获得的矩阵b1，进行第二次卷积计算，卷积核大小设为2
×2×
2，权重系数为1，步长为2，得到三维矩阵a2，然后，对a2依此进行阈值分割和腐蚀计算得到第二级粗化矩阵b2；
[0039]
步骤33、由步骤32获得的矩阵b2，进行第三次卷积计算，卷积核大小设为2
×2×
2，权重系数为1，步长为1，得到三维矩阵a3，然后，对a3依此进行阈值分割和欧式距离计算得的第三级粗化矩阵b3；
[0040]
步骤31、步骤32和步骤33中的阈值分割算法式为：
[0041][0042]
式(1)中，ai(x,y,z)为第i次卷积计算得到的单位矩阵在坐标点(x,y,z)位置的卷积值，ti为判断ai矩阵中卷积值是否满足粗化条件所设置的阈值，i＝1、2、3。
[0043]
步骤31和步骤32中，腐蚀算法将卷积矩阵中值不为0的区域向内收缩，所述腐蚀算法所采用的结构元素为3
×3×
3的全置矩阵，其腐蚀算式为：
[0044][0045]
式(2)中，bj(x,y,z)为第j次卷积得到的单位矩阵在坐标点(x,y,z)位置对应的卷积值，tj为判断aj矩阵中卷积值是否满足腐蚀条件所设置的阈值，j＝1、2。
[0046]
步骤33中，对矩阵a3中满足粗化条件的体素做欧式距离计算，其计算式为：
[0047][0048]
式(3)中，xc、yc、zc为矩阵中值不为0的体素坐标点，xb、yb、zb为距离该非零体素点最近的值为零的体素点。
[0049]
步骤34、对步骤1获得的微孔区矩阵进行第四次、第五次、第六次卷积计算，卷积核大小分别为2
×2×
1、2
×1×
2、1
×2×
2，权重系数均为1，步长为1，得到三维矩阵a4、a5和a6。
[0050]
每个卷积矩阵的值都会映射微孔区矩阵中的一个区域，卷积层数越高的卷积值所映射的区域越大，b3中的每个卷积值在微孔矩阵中映射的区域为8
×8×
8个像素。
[0051]
步骤35、由步骤33获得的矩阵b3，检索矩阵中的每个卷积值，当不为0时，对卷积值所映射的微孔区矩阵中的体素进行标记，并将b3、a2、a1、a4、a5、a6中映射区域包含已标记体素的卷积值修改为0，重复判断b3卷积值、微孔矩阵标记、修改所述b3、a2、a1、a4、a5、a6卷积值的过程，直到遍历b3矩阵。
[0052]
标记是通过修改微孔矩阵值来完成的：在微孔区矩阵中微孔体素值为1，在遍历卷积值的过程中，当找到第一个满足粗化条件的卷积值时，将其在原始图像中所映射的体素值都修改为k 1(k为所提取的孔隙网络中大孔数量，在分水岭算法所得到的分割图像中，大孔体素值是从1～k)；第二个则标记为k 2，依此类推。从k 1开始标记微孔单元，是为了提取大孔与微孔的连接关系。
[0053]
将卷积值修改为0的作用是：当再遍历到该卷积值时，就会判断其不满足粗化条件，不在参与粗化。
[0054]
遍历a2、a1、a4、a5、a6矩阵，在a2中卷积值为64时，在a1中为8时，在a4、a5、a6中为4时，分别相应地对卷积值所映射的微孔区矩阵的体素进行标记，并将a2、a1、a4、a5、a6中映射区域包含已标记体素的卷积值修改为0；重复判断、微孔矩阵标记、修改卷积值的过程，直到所有矩阵遍历完成。
[0055]
对b3矩阵所采用的遍历方法是对非零值从大到小进行检索，对矩阵a2、a1、a4、a5、a6采用沿x、y、z三个方向顺序检索。b3、a2、a1、a4、a5、a6中被修改为0的卷积值不参与之后的粗化判断。
[0056]
遍历粗化标记完的微孔矩阵，依次对该矩阵体素值为1，也就是未参与粗化的微孔体素进行标记。
[0057]
步骤4、根据步骤3获得粗化标记的微孔区矩阵和步骤2获得的分割矩阵，生成微孔区粗化后的网格单元信息和连接关系；
[0058]
粗化后微孔区的网格单元信息包括网格单元尺寸、灰度值、网格单元间的接触面积和中心点间距。
[0059]
粗化后微孔区网格单元的灰度值为在原始图像中包含的所有体素灰度值的平均值，通过该值可估计微孔区网格单元的孔隙率、渗透率等材料参数，对非均质微孔区数字岩心进行数值模拟。
[0060]
步骤5、将大孔区的孔隙网络与微孔区的网格单元合并，生成数字岩心混合模型的计算网格文件；
[0061]
所述孔隙网络是将孔隙结构简化为一种用球(孔隙)和棍(喉道)表示的孔隙网络，参见图5中的结构。所述网格单元是有限元和有限体积法模拟所使用的连续网格。
[0062]
步骤6、用步骤5生成的数字岩心混合模型的计算网格文件进行数值模拟计算。
[0063]
实施例
[0064]
岩心图像如图2所示，该数字岩心尺寸为200
×
200
×
200体素，分辨率为5.35μm，其中固体灰度值为0，微孔灰度值为1，大孔固体值为2。基于该数字岩心图像，依照以下步骤：
[0065]
步骤1、导入数字岩心图像，通过多点阈值分割得到大孔区和微孔区体素信息。图3和图4分别为大孔区和微孔区的可视化图像。从图3和图4看出，大孔连通性较好，而微孔区则成不连通的块状，填充在大孔之间。
[0066]
提取大孔区和微孔区信息，分别存储在两个与原始图像尺寸相同的体素矩阵中，构成大孔区矩阵和微孔区矩阵。
[0067]
步骤2、采用分水岭分割算法提取大孔区的孔隙网络，如图5所示，与图3对比可以看出，分水岭算法所提取的孔隙网络能很好的表征孔隙结构。
[0068]
获得孔喉的几何信息包括体积、比表面积、内切圆半径等，这些信息是通过分水岭分割算法提取到。
[0069]
步骤3、对步骤1获得的微孔区体素矩阵数据进行多重卷积计算，并微孔区体素进行粗化标记。
[0070]
1)、对微孔区矩阵中的微孔区体素进行第一次卷积计算，卷积核大小为2
×2×
2，权重系数为1，步长为2，得到三维矩阵a1，然后，对a1矩阵做阈值分割提取满足粗化条件的区域，对阈值分割后的矩阵做腐蚀计算得到第一级粗化矩阵b1；
[0071]
2)、对步骤1)获得的矩阵b1进行第二次卷积计算，卷积核大小设为2
×2×
2，权重系数为1，步长为2，得到三维矩阵a2，然后，对a2依此进行阈值分割和腐蚀计算得到第二级粗化矩阵b2；
[0072]
3)、对步骤2)获得的矩阵b2，进行第三次卷积计算，卷积核大小设为2
×2×
2，权重系数为1，步长为1，得到三维矩阵a3，然后，对a3依此进行阈值分割和欧式距离计算得的第三级粗化矩阵b3；
[0073]
4)、对微孔区矩阵的微孔区体素进行第四次-第六次卷积计算，卷积核大小分别为2
×2×
1、2
×1×
2，1
×2×
2，权重系数均为1，步长为1，得到三维矩阵a4、a5、a6；
[0074]
在未粗化的微孔区中网格单元的大小为1
×1×
1(也就是1个体素单元)，b3用于确定8
×8×
8(512个体素合并而来)的粗化网格，b2用于确定4
×4×
4的粗化网格，b1用于确定2
×2×
2的粗化网格，a4、a5、a6用于确定2
×2×
1、2
×1×
2、1
×2×
2的粗化网格。
[0075]
5)、对矩阵b3中的卷积值进行统计，按从大到小的方式进行遍历，当值不为0时，对卷积值所映射的微孔区矩阵的体素进行标记，并将b3、a2、a1、a4、a5、a6中映射区域包含已标记体素的卷积值修改为0，重复判断b3卷积值、微孔矩阵标记、修改所述b3、a2、a1、a4、a5、a6卷积值的过程，直到遍历b3矩阵；
[0076]
采用沿x，y，z方向顺序检索的方式遍历a2、a1、a4、a5、a6矩阵，在a2中卷积值为64时，在a1中卷积值为8时，在a4、a5、a6中卷积值为4时，分别相应地对卷积值所映射的微孔区矩阵体素进行标记，若a2、a1、a4、a5、a6中卷积值所映射的微孔区矩阵中包含已标记体素，则并将a2、a1、a4、a5、a6中映射区域的卷积值修改为0；重复判断、微孔矩阵标记、修改卷积
值的过程，直到所有矩阵遍历完成。
[0077]
所谓对微孔区矩阵体素进行标记，其实就是在微孔区矩阵中对矩阵值重新赋值。在遍历卷积值的过程中，当找到第一个微孔满足粗化条件的卷积值时，将其在原始图像中所映射的体素值都修改为k 1，第二卷积值的映射区域则标记为k 2，依此类推。
[0078]
遍历粗化标记完的微孔区矩阵，依次对该矩阵体素值为1(也就是未参与粗化的微孔体素)进行标记。体素值是接着粗化标记的值继续累加的。
[0079]
步骤4、对微孔粗化标记矩阵和大孔分水岭分割矩阵做并集操作，遍历并集生成的矩阵，生成粗化微孔区的网格单元信息和连接关系文件；
[0080]
并集操作为：当粗化标记矩阵和大孔分水岭分割矩阵在同一位置的值都为0时，在合并矩阵中该位置值为0，否则为两矩阵中的非零值。
[0081]
连接关系的获取是：在粗化标记完成的原始矩阵中遍历所有微孔体素，查找其与上下左右前后六个方向体素的值，当值为0时，说明相邻体素为固体，不存储任何信息。当值为1～k时，说明其与大孔相连，大于k时，则是微孔与微孔相连。此外，当两相邻体素的值相同时，说明两体素位于同一个粗化单元中，不对其连接关系进行存储。这些连接关系会储存在一个矩阵中；通过对连接关系进行统计，可得到网格单元间的接触面积。
[0082]
网格单元信息的获取是：根据网格单元在矩阵中的位置可得到单元的坐标；将粗化标记与ct图像中的灰度相结合可得到粗化后的单元灰度。
[0083]
最后生成两个矩阵，一个矩阵存储微孔网格内部信息，如单元位置、单元尺寸、单元灰度等；另一个存储微孔网格间的信息，如连接关系、单元间的接触面积等。
[0084]
步骤5、将大孔的孔网与微孔的网格单元合并，生成粗化混合模型数据文件。
[0085]
网格的合并是通过合并存储网格信息的矩阵来实现的，在合并前需确保合并矩阵的列数相同，且每列所存储的信息类型相同，然后沿列方向合并矩阵。
[0086]
粗化前和粗化后的对比如图6所示，从图6中可以看出，粗化模型中的微孔区出现了由中心向外尺寸逐级减小的网格单元，其最小网格单元尺寸均位于微孔区表面，大小为1个体素。
[0087]
图6的可视化图只能展示单元的大小及其空间关系，其详细数据采用txt文本进行保存。
[0088]
说明：步骤4的合并，是为了遍历合并的矩阵，去提取微孔与微孔，微孔与大孔的连接关系以及其他模拟需要的单元位置、单元ct值等信息。步骤5的合并是将采用分水岭算法提取的大孔与大孔的关系以及其它模拟需要的参数信息与步骤4最终得到的微孔与微孔、微孔与大孔的关系合并以及其它参数信息，最终得到一个多尺度的几何模型。
[0089]
步骤6、用生成的粗化混合模型数据文件进行单相不可压缩数值模拟计算，计算该样品的渗透率。
[0090]
单相不可压缩数值模拟计算渗透率参见文献“四川盆地典型富有机质页岩孔隙结构特征及页岩气渗流机理研究”[d]，张莉等，中国科学院大学(中国科学院广州地球化学研究所)，2021，p145。
[0091]
为了验证本方法发明的准确性，统计本发明生成粗化混合模型与初始的未粗化混合模型的计算单元数量对比如图7所示，与未粗化混合模型相比，本方法发明的微孔区网格单元减少了80％以上，而与研究区域的总体素相比，粗化混合模型的网格单元减少了98％。
经步骤6渗透率模拟计算，如图8所示，在不可压缩流动的情况下，微孔区计算网格的粗化对压力分布和渗透率的影响都非常小。