CN108876901A - 一种基于二维图像和多点统计学的数字岩心重建方法 - Google Patents
一种基于二维图像和多点统计学的数字岩心重建方法 Download PDFInfo
- Publication number
- CN108876901A CN108876901A CN201810456085.2A CN201810456085A CN108876901A CN 108876901 A CN108876901 A CN 108876901A CN 201810456085 A CN201810456085 A CN 201810456085A CN 108876901 A CN108876901 A CN 108876901A
- Authority
- CN
- China
- Prior art keywords
- ijk
- image
- point
- data template
- data
- Prior art date
- Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
- Pending
Links
- 238000000034 method Methods 0.000 title claims abstract description 66
- 239000011435 rock Substances 0.000 claims abstract description 98
- 238000012549 training Methods 0.000 claims abstract description 47
- 239000011148 porous material Substances 0.000 claims abstract description 36
- 238000004422 calculation algorithm Methods 0.000 claims abstract description 16
- 238000004088 simulation Methods 0.000 claims description 22
- 238000002591 computed tomography Methods 0.000 claims description 21
- 238000012512 characterization method Methods 0.000 claims description 16
- 238000005315 distribution function Methods 0.000 claims description 10
- 238000001914 filtration Methods 0.000 claims description 10
- 238000000342 Monte Carlo simulation Methods 0.000 claims description 9
- 238000009826 distribution Methods 0.000 claims description 9
- 238000005295 random walk Methods 0.000 claims description 6
- 238000002948 stochastic simulation Methods 0.000 claims description 6
- 238000004364 calculation method Methods 0.000 claims description 4
- 238000000605 extraction Methods 0.000 claims description 4
- 238000005266 casting Methods 0.000 claims description 2
- 238000002474 experimental method Methods 0.000 abstract description 8
- 238000010586 diagram Methods 0.000 description 26
- 239000000523 sample Substances 0.000 description 22
- 230000008569 process Effects 0.000 description 14
- 238000005516 engineering process Methods 0.000 description 7
- 230000006870 function Effects 0.000 description 5
- 230000011218 segmentation Effects 0.000 description 5
- 238000012545 processing Methods 0.000 description 4
- 238000003860 storage Methods 0.000 description 4
- 238000011161 development Methods 0.000 description 3
- 230000018109 developmental process Effects 0.000 description 3
- 239000000284 extract Substances 0.000 description 3
- 239000010410 layer Substances 0.000 description 3
- 230000007246 mechanism Effects 0.000 description 3
- 238000005192 partition Methods 0.000 description 3
- 230000035699 permeability Effects 0.000 description 3
- 238000011160 research Methods 0.000 description 3
- 238000004458 analytical method Methods 0.000 description 2
- 238000004891 communication Methods 0.000 description 2
- 238000005314 correlation function Methods 0.000 description 2
- 230000007547 defect Effects 0.000 description 2
- 230000007812 deficiency Effects 0.000 description 2
- 230000000694 effects Effects 0.000 description 2
- 239000012530 fluid Substances 0.000 description 2
- 238000003709 image segmentation Methods 0.000 description 2
- 239000002245 particle Substances 0.000 description 2
- 230000000704 physical effect Effects 0.000 description 2
- 238000012805 post-processing Methods 0.000 description 2
- 238000002922 simulated annealing Methods 0.000 description 2
- 239000007787 solid Substances 0.000 description 2
- 238000007619 statistical method Methods 0.000 description 2
- 238000012360 testing method Methods 0.000 description 2
- 230000009466 transformation Effects 0.000 description 2
- BVKZGUZCCUSVTD-UHFFFAOYSA-L Carbonate Chemical compound [O-]C([O-])=O BVKZGUZCCUSVTD-UHFFFAOYSA-L 0.000 description 1
- 230000002159 abnormal effect Effects 0.000 description 1
- 238000010521 absorption reaction Methods 0.000 description 1
- 230000009286 beneficial effect Effects 0.000 description 1
- 230000008901 benefit Effects 0.000 description 1
- 210000000988 bone and bone Anatomy 0.000 description 1
- 239000000470 constituent Substances 0.000 description 1
- 239000013068 control sample Substances 0.000 description 1
- 239000004744 fabric Substances 0.000 description 1
- 238000000227 grinding Methods 0.000 description 1
- 238000009499 grossing Methods 0.000 description 1
- 238000010191 image analysis Methods 0.000 description 1
- 238000003384 imaging method Methods 0.000 description 1
- 230000001788 irregular Effects 0.000 description 1
- 238000004519 manufacturing process Methods 0.000 description 1
- 239000011159 matrix material Substances 0.000 description 1
- 238000005259 measurement Methods 0.000 description 1
- 238000013508 migration Methods 0.000 description 1
- 230000005012 migration Effects 0.000 description 1
- 238000003672 processing method Methods 0.000 description 1
- 239000004576 sand Substances 0.000 description 1
- 239000002356 single layer Substances 0.000 description 1
- 241000894007 species Species 0.000 description 1
- 239000004575 stone Substances 0.000 description 1
- 239000000126 substance Substances 0.000 description 1
Classifications
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
- G06T17/00—Three dimensional [3D] modelling, e.g. data description of 3D objects
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
- G06T2200/00—Indexing scheme for image data processing or generation, in general
- G06T2200/08—Indexing scheme for image data processing or generation, in general involving all processing steps from image acquisition to 3D model generation
Landscapes
- Physics & Mathematics (AREA)
- Engineering & Computer Science (AREA)
- Computer Graphics (AREA)
- Geometry (AREA)
- Software Systems (AREA)
- General Physics & Mathematics (AREA)
- Theoretical Computer Science (AREA)
- Image Analysis (AREA)
Abstract
本发明涉及一种基于二维图像和多点统计学的数字岩心重建方法，本发明在物理实验方法获取岩样的二维训练图像，在得到的二维岩心图像的基础上结合多点统计学方法，获取二维岩心图像上的高阶统计信息，从而构建能够反映岩心真实内部孔隙结构特征的三维数字岩心。克服了传统数值重建算法中只应用两点统计函数的缺点。
Description
技术领域
本发明涉及一种基于二维图像和多点统计学的数字岩心重建方法，属于油气田开发工程中的数值模拟技术和数字岩心领域。
背景技术
多孔介质广泛存在于自然界和人类实验活动的许多领域。油气田开发中的储层属于一种典型的多孔介质。岩石多孔介质作为流体储集和流动的基本介质，其流动机理和输运性质源于其内部微观的孔隙结构，如孔喉大小、形状、内部粗糙程度、相互连通情况、孔喉特征及其构成方式等。准确、定量地表征岩石多孔介质微观孔隙结构对于研究岩石多孔介质中微细观渗流机理和输运性质有着重要的理论意义。
岩心常规实验方法借助实验手段来观测岩石表观性质的变化，来间接地反映岩石孔隙结构及其对物理、力学等宏观性质的影响。这种表观上的观测不能定量地解析岩石内部孔隙的连通性，也无法对决定岩石的表观性质的内部机制进行研究。随着岩石物理三维成像测试技术和计算机技术的发展，基于数字岩心的岩石物理数值模拟技术已经发展成为岩石物理研究的一种重要手段。数字岩心准确度决定了岩石物理数值模拟结果的可靠性和适用性，仅当数字岩心准确反映真实岩心微观特征时，岩石物理模拟结果才具有使用价值。
数字岩心建模方法可分为两类：物理实验方法和数值重建方法。物理实验方法借助高倍光学显微镜、扫描电镜、CT成像仪等高精度仪器获取岩心三维图像，之后对图像进行滤波去燥、分割等操作，即可得到数字岩心。物理实验方法得到的数字岩心忠实于岩心的微观孔隙结构，但具有测试时间长、费用高的缺点。相比三维孔隙图像，多孔介质的二维薄片图像较易获取并且可以得到较高的分辨率。此外，扫描分辨率和样品尺寸相互制约，若提高分辨率，样品尺寸必减少，代表性差，如果要想获得更大尺度和范围内的孔隙结构，必须通过结合其他建模方法对其特征的三维重构才能实现。
用以建立数字岩心的数值重构法通常以岩心切片二维图像为基础，借助各种不同的统计方法或模拟岩石的形成过程来建立数字岩心。迄今已发展了多种重构方法，比较典型的有高斯模拟法、模拟退火法、过程模拟法、马尔科夫连蒙特卡洛法和多点统计学方法，但上述方法重建的数字岩心具有一定的缺陷。例如，高斯模拟法重建时仅有孔隙度和两点相关函数作为约束条件，可能无法反映孔隙空间的结构特征；模拟退火法在重建时以两点函数为基础造成孔隙连通性差、孔隙几何异常、建模时间长；过程模拟法无法模拟复杂孔隙的成岩过程；马尔科夫连蒙特卡洛法无法反映孔隙的长连通性；多点统计学方法在重建各向异性的岩石结构时难以在线岩石结构不同方向上的差异性。出现上述这些缺陷的根本原因是岩石的内部过于复杂。
发明内容
针对现有技术的不足，本发明提供了一种基于二维图像和多点统计学的数字岩心重建方法。
本发明的目的是在得到的二维岩心图像的基础上结合多点统计学方法，获取二维岩心图像上的高阶统计信息，从而构建能够反映岩心真实内部孔隙结构特征的三维数字岩心。克服了传统数值重建算法中只应用两点统计函数的缺点。
术语解释：
1、表征单元体，REV，能够有效表征岩心宏观物理特性的最小岩心单元体，小于REV尺度获得岩石物理波动明显，而大于REV尺度岩石物性趋于稳定。数字图像的尺寸越大，对计算机存储能力要求越高，因此合理选取表征单元体，使其既能代表岩石的宏观性质，又适用于现有计算机条件。
2、孔隙度，是指岩样中所有孔隙空间体积之和与该岩样体积的比值，称为该岩石的总孔隙度，以百分数表示。
3、蒙特卡洛方法，即蒙特·卡罗方法(Monte Carlo method)，也称统计模拟方法，是二十世纪四十年代中期由于科学技术的发展和电子计算机的发明，而被提出的一种以概率统计理论为指导的一类非常重要的数值计算方法。是指使用随机数(或更常见的伪随机数)来解决很多计算问题的方法。与它对应的是确定性算法。蒙特·卡罗方法在金融工程学，宏观经济学，计算物理学(如粒子输运计算、量子热力学计算、空气动力学计算)等领域应用广泛。
4、中值滤波算法，是一种非线性平滑技术，将每一像素点的灰度值设置为该点某邻域窗口内的所有像素点灰度值的中值。
本发明的技术方案为：
一种基于二维图像和多点统计学的数字岩心重建方法，包括：
(1)建立二维训练图像，包括相互正交方向的3张二维图像，所述二维图像3为岩心的铸体薄片、扫描电镜或CT扫描图像的任意切面；
若仅利用单张二维图像，由于仅考虑了单层方向上的孔隙连通性，导致重建结果的连通性与真实岩心差别较大，为了保证重建岩石数字岩心的各向异性和非均质性，引入相互正交方向的三个二维图像作为训练图像。训练图像选取的优劣直接决定了重建结果的准确性，因此训练图像的选取需要有足够的代表性。
(2)选取多重数据模板，数据模板是由n个向量组成的几何形态，设为τn＝{hα；α＝1，2，…，n}。设数据模板中心位置为u，数据模板其他位置uα＝u+hα(α＝1，2，…，n)。指训练图像中的某些性质能够被在它上面移动的模板所捕获，数据模板的几何形态、大小均可变化，尺寸大小可以从3×3到100×100不等，多重数据模板在数据模板的基础上利用网格逐渐密集化的多个数据模板来代替一个大而密集的模板来对训练图像进行扫描；通过多重数据模板扫描步骤(1)建立的二维训练图像，捕获数据事件，是指通过已定义的数据模板在训练图像滑动可获得具体数据事件，这些数据事件又被称为模式特征，设为训练数据的结构特征；
数据模板的几何形状、大小均可变化，随着模板尺寸变大，对计算机存储和运算能力要求变高，模拟运行时间变长。重建三维数字岩心时，选取合适数据模板决定了重建结果的准确性。若选取数据模板较小，不能精确反映训练图像中大尺度孔隙分布情况；反之，若选取的数据模板较大，易损失小尺度孔隙信息，且极大地增加了模拟时间。多重数据模板的使用克服了这个问题。在不同方向上根据训练图像的特征选取不同的数据模板。
(3)获取条件概率：对于任一待模拟点u，确定在给定n个条件数据值S(uα)的情况下，属性S(u)取K个状态值中任一个状态值的条件概率P{S(u)＝sk|dn}分布函数如式(Ⅰ)所示：
式(Ⅰ)中，分母为某个模式出现的概率；α＝1,…,n；
分子为模式和待模拟点u取某个状态值的情况同时出现的概率；
假定一种属性S可取K个状态值{sk；k＝1,2,…K}，由数据模板中n个向量uα位置的n个状态值所组成的“数据事件”dn可以定义为：S(uα)表示在uα位置的状态值，dn表示n个向量在uα位置的S(u1),…S(un)，分别在状态值
为了避免反复扫描训练图像，采用“搜索树”的数据结构来存储条件概率以加速重建过程，同时节省内存。在应用三张正交训练图像过程中，分别在不同方向训练图像上提取带模拟点的数据事件，从相对搜索树上检索到不同方向上的不同条件概率分布函数，由于待模拟点的状态值与三个方向上的条件概率分布函数相关，将不同方向上在该处的条件概率分布函数加权平均计算得到待模拟点总条件概率，如式(Ⅱ)所示：
P(xijk|x(N30(ijk)))＝α{P(xijk|x(Ni,12(ijk)))+P(xijk|x(Nj,12(ijk)))+P(xijk|x(Nk,12(ijk)))} (Ⅱ)
式(Ⅱ)中，P(xijk|x(N30(ijk)))表示待模拟点的条件概率；
x为任意的状态值；
xijk为(i，j，k)坐标系的某一状态值；
N30(ijk)为在三维数据模板下对应x状态出现的概率数；
Ni,12(ijk)为在XY平面数据模板下对应x状态出现的概率数；
Nj,12(ijk)为在YZ平面数据模板下对应x状态出现的概率数；
Nk,12(ijk)为在ZX平面数据模板下对应x状态出现的概率数；
i表示X方向位置坐标；j表示Y方向位置坐标；k表示Z方向位置坐标；
P(xijk|x(Ni,12(ijk)))、P(xijk|x(Nj,12(ijk)))、P(xijk|x(Nk,12(ijk)))分别依次表示代表XY、YZ和ZX平面上的条件概率分布；
α为平均参数；进一步优选的，α＝1/3。
(4)定义一条能将所有的待模拟点遍历一次的随机路径，访问待模拟区域的待模拟点；对随机路径上的每一个待模拟点，通过上述步骤(2)-(3)，利用与步骤(2)选取的多重数据模板提取其数据事件，从搜索树中获取该待模拟点的总条件概率；利用蒙特卡洛方法提取该待模拟点的随机模拟值(即孔隙或者骨架)，并将该随机模拟值作为后续模拟新增待模拟点的条件数据；即：当待模拟区域中的某一模拟点确定为随机模拟值时，为了保证后续待模拟点与该点的连通性以及相关性，将该模拟点作为已知条件数据来约束后续的概率分布函数，直至确定待模拟区域的所有待模拟点，得到岩样的三维数字岩心。
根据本发明优选的，建立二维训练图像，包括：
a、选取表征单元体：在岩心二值图像中任意定位一个体素；以该体素为中心点选取一个立方体，该立方体的边长选取的初始值一般没有固定要求，但是不宜过小，也不宜过大，一般选取二值图像边长的十分之一，统计其孔隙度；逐渐扩大立方体边长，并计算每一个边长岩心立方体的孔隙度；作出孔隙度随岩心二值图像边长的变化规律曲线；定位不同的体素中心，重复该步骤a，当所有的孔隙度随岩心二值图像边长的变化规律曲线均趋于稳定时，对应的最小岩心边长即为表征单元体尺寸；
b、选取合理的二维训练图像：绘制不同方向(X、Y、Z)上孔隙度与层序号的关系曲线，选取与表征单元体孔隙度接近的层且孔隙结构平稳、具有代表性的为二维训练图像。
根据本发明优选的，通过多重数据模板扫描步骤(1)建立的二维训练图像，包括：
A、使用稀疏的粗网格数据模板(粗网格的选取与训练图像大小和孔隙结构复杂性相关，要求粗网格尽量覆盖训练图像中大孔隙结构特征，一般选取的粗网格为21×21×21)扫描训练图像，得到粗网格数据模板下的多点统计信息，多点统计信息即数据事件(模式结构)和条件概率分布函数，以搜索树的形式存储；
B、模拟得到粗网格数据模板下的结果图像，即：利用多点统计学方法结合粗网格数据模板下的多点统计信息(数据事件和条件概率分布函数)，重建出粗网格模板下的图像；
C、将粗网格数据模板下的结果图像作为条件数据复制到细网格数据模板(细网格的选取与训练图像大小和孔隙结构复杂性相关，要求细网格尽量获取训练图像中微孔隙结构特征，一般选取的粗网格为5×5×5)上，利用细网格数据模板扫描训练图像，得到细网格数据模板下的多点统计信息；
D、模拟得到细网格数据模板下的结果图像，即：利用多点统计学方法结合细网格数据模板下的多点统计信息(数据事件和条件概率分布函数)重建出细网格模板下的图像；
根据本发明优选的，所述步骤(4)之后，执行以下步骤：对步骤(4)得到的岩样的三维数字岩心，采用中值滤波算法去除噪声。
去除噪音和不合理的孔喉结构、分布后即可得到最终结果。
本发明的有益效果为：
传统数值重建方法中，往往只利用孔隙度、两点相关函数等低阶统计信息，造成重建的数字岩心孔隙结构异常、连通性较差、模拟时间长等缺点，本发明能够弥补在应用物理实验方法获取岩样三维数字岩心过程中分辨率和物理尺寸相互矛盾的不足，本发明突破了物理实验方法获取数字岩心在物理尺寸方面的限制，同时将获取的岩样的高阶统计信息应用到岩样数字岩心重建过程中，该发明能够应用到常规砂岩、非均质性较强的碳酸盐岩、以及孔喉尺寸小的非常规致密砂岩、泥页岩等不同的岩样，为研究油气田开发过程中储层岩石物理属性的数值模拟和流动模拟计算提供了更好的研究平台，因此极具推广价值，在目前公开发表文献和商业应用软件中尚无类似的提出与应用。
附图说明
图1(a)是Berea砂岩某一切面X射线的灰度图像示意图。
图1(b)是Berea砂岩某一切面二值化图像示意图。
图2是Berea砂岩X射线数据处理流程图。
图3是不同方向的训练图像示意图。
图4(a)是数据模板扫描整个训练图像的示意图。
图4(b)是捕获岩样的模式库的示意图。
图5是重建数字岩心示意图。
图6(a)是重建的数字岩心示意图。
图6(b)是X射线CT扫描得到的数字岩心示意图。
图7(a)是提取的重建数字岩心的孔隙网络模型示意图。
图7(b)是提取的X射线CT扫描的孔隙网络模型示意图。
图8(a)是重建和CT扫描数字岩心的孔隙半径对比图。
图8(b)是重建和CT扫描数字岩心的吼道半径对比图。
图8(c)是重建和CT扫描数字岩心的孔喉比对比图。
图8(d)是重建和CT扫描数字岩心的配位数对比图。
图8(e)是重建和CT扫描数字岩心的孔隙形状因子对比图。
图8(f)是重建和CT扫描数字岩心的吼道形状因子对比图。
具体实施方式
下面结合说明书附图和实施例对本发明作进一步限定，但不限于此。
实施例
一种基于二维图像和多点统计学的数字岩心重建方法，该实施例为X射线CT得到的Berea砂岩，包括：
(1)处理Berea岩样，为了保证岩心CT图像能够清晰反映岩样内部微观结构特征，将Berea岩心钻裁为直径最大为厘米级的岩心柱体，随后对岩心柱体进行打磨和抛光，减低样品表面粗糙度，提高成像质量；
(2)对Berea砂岩岩样开展X射线CT扫描，开展扫描实验的设备为Zeiss公司生产的MicroXCT-400型CT机，其最高分辨率为0.75微米，视域像素为20482。
CT扫描的工作过程如下：将样品固定好，开启X射线源，由射线源发出的射线穿过样品，该过程中X射线强度衰减，衰减后的X射线照射到探测器上，该信号被图像获取软件自动捕获并存储。通过测量物质对射线的吸收系数可以判定物质的组成成分，随后，通过控制样品夹持器将样品精确旋转一定角度，重新扫描并记录衰减后的X射线，将样品累计旋转180后结束实验过程。在完成上述扫描过程后，可以获得岩心的投影数据，然后再由投影数据可以重建岩心灰度图像。图像重建的实质是由采集后的数据求解图像矩阵中各像素的吸收系数，然后重新构造图像。
本案例中对Berea砂岩开展X射线CT扫描过程中的分辨率为6.50微米。Berea砂岩某一切面X射线的灰度图像如图1(a)所示。
(3)灰度图像滤波去燥及二值化分割，重建得到的岩心灰度图像后，需要判定孔隙和骨架的具体分布。由于外界干扰导致图像采集过程中所采集的图像出现偏暗或偏亮、对比度不明显、图像模糊等缺点，将给采集灰度图像分析带来较大的误差，所以在不破坏图像所含有的有用信心的前提下有必要对图像进行预处理以改善其质量。一般情况下，图像的预处理是除去无用信息，突出有用信息。针对所采集岩样灰度图像存在偏暗或偏亮的情况，可以通过调整图像的亮度来改变显示效果。判断图像亮度是否合理须借助图像的灰度直方图，灰度直方图给出了图像中所有像素点的灰度值在[0,255]不同灰度区间的概率，图像亮度调整是一种点处理，即对图像中的每一个像素点的灰度增加、减一个常数即可，通过图像亮度调整，灰度直方图只是整体发生了平移而形状保持不变。
针对图像对比度不够明显或是对比度过大的情况，需对图像对比度加以调整，调整时仍需借助灰度直方图。产生这种情况原因主要由图像采集设备曝光不足、灰度值设置不合理等。用以调整图像对比度的方法有多种，如线性变换、非线性变换、直方图的平坦化处理等。
为了突出图像的边缘信息，加强图像的边缘特征，以便于人眼和机器的识别，需要对图像进行锐化处理。用这种方法可以去掉引起图像质量低下的“模糊”现象，并且能够把图像变的轮廓分明。锐化处理可以分为空域处理方法和频域处理方法两种。而空域锐化中常用的有梯度算子锐化、拉普拉斯算子锐化和模板卷积锐化等方法。一般采用拉普拉斯算子锐化法。
经过处理后的灰度图像中存在各种类型的系统噪声，因此需要通过滤波算法来去除这些噪声。常用的滤波算法有高斯滤波、中值滤波、均值滤波、非线性滤波等。均值滤波是最为常用的平滑噪声技术，即用一像素领域内所有像素的灰度平均值代替该像素原来的灰度。
在对灰度图像去噪后，孔隙与岩石骨架的边缘需要借助图像分割方法对其进行合理划分。图像分割的关键在于分割阀值的选取，如果阀值选取不合理，很容易将目标对象理解成为背景图像，或将背景理解为目标对象，这样在对目标图像进行特征分析时就无法有效分析目标特征。常用的图像二值分割方法有迭代阀值法、最大类间距法、指示克里金法、分水岭算法等。
本实施例使用最大类间距法来进行图像分割，该方法分割的思路如下：将待分割图像中所有像素点的灰度值所构成的总体视为一个集合，假定把灰度值的集合用阀值分为两组，当两组的平均值的方差(组间方差)和各组的方差(组内方差)之比达到最大时，其所对应的分割阀值即为合理值。显然，岩心内部的孔隙空间和岩石古交可以各成一组，每一组内部的差异很小而两组间的性质差异很大；当采用合理的阀值将两组分割开后，两组间的性质差异与组内差异之比将达到极大值。最后将阀值分割后的图像进行堆叠得到三维数字岩心，其文件为一个三维数组。本实施例Berea砂岩某一切面二值化图像如图1(b)所示。
上述步骤(1)-(3)Berea砂岩X射线数据处理流程如图2所示。
(4)选取训练图像，能够有效表征岩心宏观物理特性的最小岩心单元体通常被称为表征单元体(REV)，小于REV尺度获得岩石物理波动明显，而大于REV尺度岩石物性趋于稳定。数字图像的尺寸越大，对计算机存储能力要求越高，因此合理选取表征单元体，使其既能代表岩石的宏观性质，又适用于现有计算机条件。选取表征单元体的方法是：首先，在岩心二值图像中任意定位一个体素；以该体素为中心点选取一定边长的立方体，统计其孔隙度；逐渐扩大立方体边长，并计算每一个边长岩心立方体的孔隙度；作出孔隙度随岩心数字图像边长的变化规律曲线；定位不同的体素中心，重复上述步骤。当所有的曲线均趋于稳定时对应的最小岩心边长即为表征单元体尺寸。
在选取的表征单元体的基础上选取合理的训练图像，其选取的方法是：首先，绘制不同方向(X、Y、Z)上孔隙度和层序号的关系曲线，选取与表征单元体孔隙度接近的层且孔隙结构正常、具有代表性的为训练图像。不同方向的训练图像如图3所示。
(5)选取多重数据模板，由于考虑了三个正交方向的训练图像，所以在重建过程中需要在三个不同方向上选取对应的数据模板。数据模板的选取至关重要，需要选择合适的模板尺寸，既体现待模拟层的变化，又能体现孔隙和骨架的空间连续性。本发明采用多重数据模板的思路，利用网格密集多个数据模板来代替一个大而密集的模板对其进行扫描。具体方法是先使用稀疏的粗网格数据模板扫描训练图像，得到粗网格下的多点统计信息，然后可以模拟得到粗网格下的结果图像；将粗网格下的内容作为条件数据复制到细网格上，然后利用细网格模板扫描训练图像，得到细网格下的多点统计信息，最后模拟得到细网格下的结果图像。不同方向的数据模板应该与该方向的训练图像相匹配，不同方向的数据模板大小、形状可以不同。根据大量重建实例发现：一般数据模板选定尺寸为训练图像尺寸的3％～5％左右。数据模板扫描整个训练图像的示意图如图4(a)所示。最终捕获岩样的模式库的示意图如图4(b)所示。
(6)随机构建，定义一条随机路径访问所有待模拟节点，对随机路径上的每一个模拟节点利用蒙特卡洛方法提取一个值作为该处的具体值，由于该点的概率有三个方向上的概率函数融合而成，将该值与概率值进行比较，得到待模拟节点的数值解，随后将该确定的模拟节点作为硬数据更新条件概率。重建数字岩心如图5所示。重建的数字岩心示意图如图6(a)所示，X射线CT扫描得到的数字岩心示意图如图6(b)所示。提取的重建数字岩心的孔隙网络模型如图7(a)所示，提取的X射线CT扫描的孔隙网络模型如图7(b)所示。
(7)重建数字岩心的后处理，从二维训练图像生成的三维数字岩心，不可避免的出现一些噪声的像素点，这样的噪声数据可以分为两类，一类是在孔隙空间中游离的一些小块的固体颗粒，这种数据在真实的多孔介质送是不存在的，另一类是在孔隙空间与固体骨架交界处，存在一些不规整的毛刺。产生噪声数据的原因有两个方面，一方面模拟图像仅仅是真实多孔介质的一个离散的原因，另一方面是在生成图像的过程中，某些采样点的模拟值实际上是使用随机方法生成的。对于重建数字岩心的后处理采用中值滤波算法来去除噪声，去除噪音和不合理的孔喉结构、分布后即可得到最终结果。
(8)评价重建数字岩心的优劣，由于利用X射线CT已经得到了Berea砂岩岩样的真实数字岩心，将重建的数字岩心和真实的数字岩心在孔隙度、渗透率、孔喉分布、拓扑结构、连通性函数以及孔隙网络模型方面进行比较，以判断重建结构的优劣。
本发明采用统计方法来计算数字岩心的孔隙度；对于渗透率的计算采用格子Boltzmann方法，该方法的主要步骤是碰撞和迁移，以描述在微观上连续，宏观上微小的流体流动过程，格子Boltzmann方法是一种不同于传统数值方法的流体计算方法，它本身具有算法简单、并行容易、处理复杂边界条件等有点，本发明中的计算模型为三维19速度的D3Q19模型；孔隙网络模型是指与数字岩心的孔隙空间具有类似的拓扑结构和几何特征的等价模型，对于孔隙网络模型的提取目前有以形态学为基础的中轴线算法和拓扑学为基础的最大球算法，本发明中采用最大球算法来提取数字岩心的孔隙网络模型。最大球算法的思路是：首先，建立孔隙空间的内切球，并且删除冗余球，随后经过成簇算法建立孔喉空间，经过孔喉分割和拓扑结构计算，建立孔隙网络模型。
图8(a)是重建和CT扫描数字岩心的孔隙半径对比图。图8(b)是重建和CT扫描数字岩心的吼道半径对比图。图8(c)是重建和CT扫描数字岩心的孔喉比对比图。图8(d)是重建和CT扫描数字岩心的配位数对比图。图8(e)是重建和CT扫描数字岩心的孔隙形状因子对比图。图8(f)是重建和CT扫描数字岩心的吼道形状因子对比图。
经过比较发现，无论是孔隙度、渗透率还是孔喉尺寸分布的孔隙网络模型特征，利用本发明重建的三维数字岩心与X射线扫描得到的数字岩心具有较强的一致性。
Claims (5)
1.一种基于二维图像和多点统计学的数字岩心重建方法，其特征在于，包括：
(1)建立二维训练图像，包括相互正交方向的3张二维图像，所述二维图像3为岩心的铸体薄片、扫描电镜或CT扫描图像的任意切面；
(2)选取多重数据模板，通过多重数据模板扫描步骤(1)建立的二维训练图像，捕获数据事件，这些数据事件又被称为模式特征，设为训练数据的结构特征；
(3)获取条件概率：对于任一待模拟点u，确定在给定n个条件数据值S(uα)的情况下，属性S(u)取K个状态值中任一个状态值的条件概率P{S(u)＝sk|dn}分布函数如式(Ⅰ)所示：
式(Ⅰ)中，分母为某个模式出现的概率；α＝1,…,n；
分子为模式和待模拟点u取某个状态值的情况同时出现的概率；
将不同方向上在该处的条件概率分布函数加权平均计算得到待模拟点总条件概率，如式(Ⅱ)所示：
P(xijk|x(N30(ijk)))＝α{P(xijk|x(Ni,12(ijk)))+P(xijk|x(Nj,12(ijk)))+P(xijk|x(Nk,12(ijk)))} (Ⅱ)
式(Ⅱ)中，P(xijk|x(N30(ijk)))表示待模拟点的条件概率；
x为任意的状态值；
xijk为(i，j，k)坐标系的某一状态值；
N30(ijk)为在三维数据模板下对应x状态出现的概率数；
Ni,12(ijk)为在XY平面数据模板下对应x状态出现的概率数；
Nj,12(ijk)为在YZ平面数据模板下对应x状态出现的概率数；
Nk,12(ijk)为在ZX平面数据模板下对应x状态出现的概率数；
i表示X方向位置坐标；j表示Y方向位置坐标；k表示Z方向位置坐标；
P(xijk|x(Ni,12(ijk)))、P(xijk|x(Nj,12(ijk)))、P(xijk|x(Nk,12(ijk)))分别表示XY、YZ和ZX平面上的条件概率分布；
α为平均参数；
(4)定义一条能将所有的待模拟点遍历一次的随机路径，访问待模拟区域的待模拟点；对随机路径上的每一个待模拟点，通过上述步骤(2)-(3)，利用与步骤(2)选取的多重数据模板提取其数据事件，从搜索树中获取该待模拟点的总条件概率；利用蒙特卡洛方法提取该待模拟点的随机模拟值，并将该随机模拟值作为后续模拟新增待模拟点的条件数据；即：当待模拟区域中的某一模拟点确定为随机模拟值时，将该模拟点作为已知条件数据来约束后续的概率分布函数，直至确定待模拟区域的所有待模拟点，得到岩样的三维数字岩心。
2.根据权利要求1所述的一种基于二维图像和多点统计学的数字岩心重建方法，其特征在于，α＝1/3。
3.根据权利要求1所述的一种基于二维图像和多点统计学的数字岩心重建方法，其特征在于，建立二维训练图像，包括：
a、选取表征单元体：在岩心二值图像中任意定位一个体素；以该体素为中心点选取一个立方体，统计其孔隙度；逐渐扩大立方体边长，并计算每一个边长岩心立方体的孔隙度；作出孔隙度随岩心二值图像边长的变化规律曲线；定位不同的体素中心，重复该步骤a，当所有的孔隙度随岩心二值图像边长的变化规律曲线均趋于稳定时，对应的最小岩心边长即为表征单元体尺寸；
b、选取合理的二维训练图像：绘制不同方向(X、Y、Z)上孔隙度与层序号的关系曲线，选取与表征单元体孔隙度接近的层且孔隙结构平稳、具有代表性的为二维训练图像。
4.根据权利要求1所述的一种基于二维图像和多点统计学的数字岩心重建方法，其特征在于，通过多重数据模板扫描步骤(1)建立的二维训练图像，包括：
A、使用粗网格数据模板扫描训练图像，得到粗网格数据模板下的多点统计信息，多点统计信息即数据事件和条件概率分布函数，以搜索树的形式存储；
B、模拟得到粗网格数据模板下的结果图像，利用多点统计学方法结合粗网格数据模板下的多点统计信息，重建出粗网格模板下的图像；
C、将粗网格数据模板下的结果图像作为条件数据复制到细网格数据模板上，利用细网格数据模板扫描训练图像，得到细网格数据模板下的多点统计信息；
D、模拟得到细网格数据模板下的结果图像：利用多点统计学方法结合细网格数据模板下的多点统计信息重建出细网格模板下的图像。
5.根据权利要求1-4任一所述的一种基于二维图像和多点统计学的数字岩心重建方法，其特征在于，所述步骤(4)之后，执行以下步骤：对步骤(4)得到的岩样的三维数字岩心，采用中值滤波算法去除噪声。
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN201810456085.2A CN108876901A (zh) | 2018-05-14 | 2018-05-14 | 一种基于二维图像和多点统计学的数字岩心重建方法 |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN201810456085.2A CN108876901A (zh) | 2018-05-14 | 2018-05-14 | 一种基于二维图像和多点统计学的数字岩心重建方法 |
Publications (1)
| Publication Number | Publication Date |
|---|---|
| CN108876901A true CN108876901A (zh) | 2018-11-23 |
Family
ID=64333762
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN201810456085.2A Pending CN108876901A (zh) | 2018-05-14 | 2018-05-14 | 一种基于二维图像和多点统计学的数字岩心重建方法 |
Country Status (1)
| Country | Link |
|---|---|
| CN (1) | CN108876901A (zh) |
Cited By (17)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN110599473A (zh) * | 2019-09-06 | 2019-12-20 | 中国科学院地质与地球物理研究所 | 一种数字岩芯确定方法、装置及设备 |
| CN111415407A (zh) * | 2020-03-27 | 2020-07-14 | 西北民族大学 | 一种采用多模板系统提升三维重建图像性能的方法 |
| CN111724331A (zh) * | 2019-03-22 | 2020-09-29 | 四川大学 | 一种基于生成网络的多孔介质图像重建方法 |
| CN111784839A (zh) * | 2020-04-17 | 2020-10-16 | 中国科学院力学研究所 | 一种rev连通孔隙空间的构建方法、装置 |
| CN111784724A (zh) * | 2020-05-28 | 2020-10-16 | 长安大学 | 改进型马尔科夫链蒙特卡洛二维岩石切片重构方法及系统 |
| CN111950192A (zh) * | 2020-07-15 | 2020-11-17 | 中海油田服务股份有限公司 | 基于卷积神经网络的孔隙网络模型的建模方法及装置 |
| CN112037124A (zh) * | 2019-06-04 | 2020-12-04 | 中国石油化工股份有限公司 | 基于图像纹理合成的特征可调控的数字岩心重构方法 |
| CN112381916A (zh) * | 2020-12-08 | 2021-02-19 | 西南石油大学 | 一种利用二维薄片图像的数字岩心三维结构重建方法 |
| CN112381845A (zh) * | 2020-12-02 | 2021-02-19 | 中国石油大学(华东) | 岩心图像生成方法、模型训练方法及装置 |
| CN112669243A (zh) * | 2020-12-09 | 2021-04-16 | 山东省科学院海洋仪器仪表研究所 | 一种基于连通性与孔隙度的岩石采样方法 |
| CN113029911A (zh) * | 2021-03-31 | 2021-06-25 | 中国科学院武汉岩土力学研究所 | 一种岩石孔隙度计算方法 |
| CN113281239A (zh) * | 2021-06-18 | 2021-08-20 | 中国石油大学(北京) | 多尺度煤岩孔隙网络生成方法和装置 |
| WO2021179558A1 (zh) * | 2020-03-13 | 2021-09-16 | 中国石油大学(华东) | 一种构建数字岩心的方法及系统 |
| CN114078183A (zh) * | 2021-11-01 | 2022-02-22 | 清华大学 | 多孔介质三维结构的重建方法、装置、设备及介质 |
| CN114764843A (zh) * | 2021-01-12 | 2022-07-19 | 四川大学 | 基于邻域块匹配的多孔介质图像超维重建方法 |
| CN115993376A (zh) * | 2022-12-06 | 2023-04-21 | 东北石油大学 | 一种基于随机生长法的页岩基质数字岩心重构方法 |
| RU2812143C1 (ru) * | 2020-07-15 | 2024-01-23 | Чайна Ойлфилд Сервисез Лимитед | Способ и устройство для измерения характеристик колонки породы для создания модели поровой системы |
-
2018
- 2018-05-14 CN CN201810456085.2A patent/CN108876901A/zh active Pending
Non-Patent Citations (3)
| Title |
|---|
| 刁庆雷: "基于多点统计法构建三维数字岩心的研究", 《中国优秀博硕士学位论文全文数据库(硕士) 基础科学辑》 * |
| 刘学锋: "基于数字岩心的岩石声电特性微观数值模拟研究", 《中国优秀博硕士学位论文全文数据库(博士) 基础科学辑》 * |
| 张挺: "基于多点地质统计的多孔介质重构方法及实现", 《中国优秀博硕士学位论文全文数据库(博士) 基础科学辑》 * |
Cited By (29)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN111724331B (zh) * | 2019-03-22 | 2023-05-09 | 四川大学 | 一种基于生成网络的多孔介质图像重建方法 |
| CN111724331A (zh) * | 2019-03-22 | 2020-09-29 | 四川大学 | 一种基于生成网络的多孔介质图像重建方法 |
| CN112037124A (zh) * | 2019-06-04 | 2020-12-04 | 中国石油化工股份有限公司 | 基于图像纹理合成的特征可调控的数字岩心重构方法 |
| CN110599473B (zh) * | 2019-09-06 | 2022-04-15 | 中国科学院地质与地球物理研究所 | 一种数字岩芯确定方法、装置及设备 |
| CN110599473A (zh) * | 2019-09-06 | 2019-12-20 | 中国科学院地质与地球物理研究所 | 一种数字岩芯确定方法、装置及设备 |
| US11934488B2 (en) | 2020-03-13 | 2024-03-19 | China University Of Petroleum (East China) | Method and system for constructing digital rock |
| WO2021179558A1 (zh) * | 2020-03-13 | 2021-09-16 | 中国石油大学(华东) | 一种构建数字岩心的方法及系统 |
| CN111415407B (zh) * | 2020-03-27 | 2023-04-07 | 西北民族大学 | 一种采用多模板系统提升三维重建图像性能的方法 |
| CN111415407A (zh) * | 2020-03-27 | 2020-07-14 | 西北民族大学 | 一种采用多模板系统提升三维重建图像性能的方法 |
| CN111784839A (zh) * | 2020-04-17 | 2020-10-16 | 中国科学院力学研究所 | 一种rev连通孔隙空间的构建方法、装置 |
| CN111784839B (zh) * | 2020-04-17 | 2023-02-07 | 中国科学院力学研究所 | 一种rev连通孔隙空间的构建方法、装置 |
| CN111784724A (zh) * | 2020-05-28 | 2020-10-16 | 长安大学 | 改进型马尔科夫链蒙特卡洛二维岩石切片重构方法及系统 |
| CN111784724B (zh) * | 2020-05-28 | 2023-05-09 | 长安大学 | 改进型马尔科夫链蒙特卡洛二维岩石切片重构方法及系统 |
| RU2812143C1 (ru) * | 2020-07-15 | 2024-01-23 | Чайна Ойлфилд Сервисез Лимитед | Способ и устройство для измерения характеристик колонки породы для создания модели поровой системы |
| WO2022011894A1 (zh) * | 2020-07-15 | 2022-01-20 | 中海油田服务股份有限公司 | 基于卷积神经网络的孔隙网络模型的建模方法及装置 |
| CN111950192B (zh) * | 2020-07-15 | 2023-10-24 | 中海油田服务股份有限公司 | 基于卷积神经网络的孔隙网络模型的建模方法及装置 |
| CN111950192A (zh) * | 2020-07-15 | 2020-11-17 | 中海油田服务股份有限公司 | 基于卷积神经网络的孔隙网络模型的建模方法及装置 |
| CN112381845A (zh) * | 2020-12-02 | 2021-02-19 | 中国石油大学(华东) | 岩心图像生成方法、模型训练方法及装置 |
| CN112381845B (zh) * | 2020-12-02 | 2023-04-18 | 中国石油大学(华东) | 岩心图像生成方法、模型训练方法及装置 |
| CN112381916A (zh) * | 2020-12-08 | 2021-02-19 | 西南石油大学 | 一种利用二维薄片图像的数字岩心三维结构重建方法 |
| CN112669243A (zh) * | 2020-12-09 | 2021-04-16 | 山东省科学院海洋仪器仪表研究所 | 一种基于连通性与孔隙度的岩石采样方法 |
| CN114764843B (zh) * | 2021-01-12 | 2023-04-18 | 四川大学 | 基于邻域块匹配的多孔介质图像超维重建方法 |
| CN114764843A (zh) * | 2021-01-12 | 2022-07-19 | 四川大学 | 基于邻域块匹配的多孔介质图像超维重建方法 |
| CN113029911A (zh) * | 2021-03-31 | 2021-06-25 | 中国科学院武汉岩土力学研究所 | 一种岩石孔隙度计算方法 |
| CN113281239B (zh) * | 2021-06-18 | 2022-06-24 | 中国石油大学(北京) | 多尺度煤岩孔隙网络生成方法和装置 |
| CN113281239A (zh) * | 2021-06-18 | 2021-08-20 | 中国石油大学(北京) | 多尺度煤岩孔隙网络生成方法和装置 |
| CN114078183A (zh) * | 2021-11-01 | 2022-02-22 | 清华大学 | 多孔介质三维结构的重建方法、装置、设备及介质 |
| CN115993376A (zh) * | 2022-12-06 | 2023-04-21 | 东北石油大学 | 一种基于随机生长法的页岩基质数字岩心重构方法 |
| CN115993376B (zh) * | 2022-12-06 | 2023-09-15 | 东北石油大学 | 一种基于随机生长法的页岩基质数字岩心重构方法 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| CN108876901A (zh) | 一种基于二维图像和多点统计学的数字岩心重建方法 | |
| CN107449707B (zh) | 页岩储层中不同尺度孔隙定量的三维表征确定方法和装置 | |
| CN108053417B (zh) | 一种基于混合粗分割特征的3D U-Net网络的肺分割装置 | |
| CN104335042B (zh) | 从多孔介质的数字表征中选择代表性单元体积的高效方法 | |
| US11590708B2 (en) | Three-dimensional fluid micromodels | |
| CN109285222A (zh) | 有机页岩高分辨率数字岩心构建与分析方法 | |
| Zhao et al. | An integrated method for 3D reconstruction model of porous geomaterials through 2D CT images | |
| CN105806765A (zh) | 一种显微ct扫描土体空间孔隙结构的精细化表征方法 | |
| CN103698803A (zh) | 一种岩石孔隙结构表征方法及装置 | |
| Friese et al. | Analysis of tomographic mineralogical data using YaDiV—Overview and practical case study | |
| CN108122221A (zh) | 弥散加权成像图像中脑缺血区域的分割方法及装置 | |
| CN113609696A (zh) | 基于图像融合的多尺度多组分数字岩心构建方法及系统 | |
| CN110660051B (zh) | 一种基于导航金字塔的张量投票处理方法 | |
| CN116403211A (zh) | 一种基于单细胞病理图像细胞核的分割和聚类方法及系统 | |
| Mahanta et al. | Digital rock physics and application of high-resolution micro-CT techniques for geomaterials | |
| CN114705606A (zh) | 一种基于网络化分析的岩石内部关键渗流节点的封堵方法 | |
| Marques et al. | Deep learning-based pore segmentation of thin rock sections for aquifer characterization using color space reduction | |
| CN114428040A (zh) | 一种页岩油储层储渗空间定量表征及参数获取方法 | |
| CN117036635B (zh) | 一种基于图像纹理分类的页岩多尺度数字岩心构建方法 | |
| Wang et al. | Improved skeleton extraction method considering surface feature of natural micro fractures in unconventional shale/tight reservoirs | |
| Zhang et al. | 3D pore space reconstruction using deep residual deconvolution networks | |
| CN111986078A (zh) | 一种基于引导数据的多尺度岩心ct图像融合重建的方法 | |
| Wu et al. | GPU-based adaptive data reconstruction for large-scale statistical visualization | |
| Sundaram | Permeability and electrical conductivity of rocks hosting multimodal pore systems and fractures | |
| Roldan | Reconstruction of porous structures from FIB-SEM data |
Legal Events
| Date | Code | Title | Description |
|---|---|---|---|
| PB01 | Publication | ||
| PB01 | Publication | ||
| SE01 | Entry into force of request for substantive examination | ||
| SE01 | Entry into force of request for substantive examination | ||
| RJ01 | Rejection of invention patent application after publication | ||
| RJ01 | Rejection of invention patent application after publication |
Application publication date: 20181123 |