全文快速搜索:   高级搜索

  中国石油大学学报(自然科学版)  2018, Vol. 42 Issue (3): 12-21  DOI:10.3969/j.issn.1673-5005.2018.03.002
0

引用本文 [复制中英文]

吴玉其, 林承焰, 任丽华, 等. 基于多点地质统计学的数字岩心建模[J]. 中国石油大学学报(自然科学版), 2018, 42(3): 12-21. DOI: 10.3969/j.issn.1673-5005.2018.03.002.
[复制中文]
WU Yuqi, LIN Chengyan, REN Lihua, et al. Digital core modeling based on multiple-point statistics[J]. Journal of China University of Petroleum (Edition of Natural Science), 2018, 42(3): 12-21. DOI: 10.3969/j.issn.1673-5005.2018.03.002.
[复制英文]

基金项目

国家科技重大专项(2016ZX05054012, 2017ZX05009001);中国石油大学(华东)自主创新科研计划项目(18CX06024A)

作者简介

吴玉其(1992-), 男, 博士研究生, 研究方向为数字岩心与微观剩余油分布机理。E-mail:wuyuqi150348@163.com

通讯作者

林承焰(1963-), 男, 教授, 博士生导师, 研究方向为油藏描述与储层地质学。E-mail:ycdzycms2017@126.com

文章历史

收稿日期:2017-10-05
基于多点地质统计学的数字岩心建模
吴玉其1,2 , 林承焰1,2 , 任丽华1,2 , 闫伟超1 , 王杨1,2 , 陈仕臻1,2 , 由春梅3 , 张丽4     
1. 中国石油大学(华东)地球科学与技术学院, 山东青岛 266580;
2. 山东省油藏地质重点实验室, 山东青岛 266580;
3. 大庆油田有限责任公司勘探开发研究院, 黑龙江大庆 163712;
4. 山东科技大学资源与土木工程系, 山东泰安 271019
摘要: 常用的数字岩心建模方法有X射线CT扫描法、模拟退火法和过程法。X射线CT扫描法虽然能够建立准确的三维数字岩心模型, 但是其价格贵、周期长; 模拟退火法和过程法能够降低建模成本提高建模效率, 不过模拟退火法不能有效地表征孔隙结构的长距离连通性, 过程法不适用于模拟成岩作用复杂的岩石。为了解决以上问题, 提出基于多点地质统计学重构数字岩心模型的方法。以密西根盆地Waverly组的Berea砂岩为例, 从X射线CT扫描法构建的Berea砂岩模型中提取两个1503体素的代表体积元, 令个为训练图像, 另个为真实模型; 把真实模型的两张切片作为条件数据, 使用多点地质统计学法重构出三维数字岩心模型。通过以变差函数、孔隙结构参数和单相流渗流参数为评价指标, 对重构模型的准确性及该方法的可行性进行研究。结果表明:重构模型与真实模型的孔喉几何形态和拓扑结构具有较高的致性, 证明了重构模型的准确性和多点地质统计学法建立数字岩心模型的可行性。
关键词: 数字岩心    多点地质统计学    孔隙结构    变差函数    孔隙网络模型    
Digital core modeling based on multiple-point statistics
WU Yuqi1,2 , LIN Chengyan1,2 , REN Lihua1,2 , YAN Weichao1 , WANG Yang1,2 , CHEN Shizhen1,2 , YOU Chunmei3 , ZHANG Li4     
1. School of Geosciences in China University of Petroleum(East China), Qingdao 266580, China;
2. Reservoir Geology Key Laboratory of Shandong Province, Qingdao 266580, China;
3. Exploration and Development Research Institute of Daqing Oilfield Limited Company, Daqing 163712, China;
4. Department of Resources and Civil Engineering in Shandong University of Science and Technology, Tai'an 271019, China
Abstract: There are three commonly used digital core modeling methods:X-ray computed tomography scanning, simulated annealing and process-based method.Although the X-ray computed tomography scanning method can build an accurate 3D digital core model, it is expensive and time-consuming.The simulated annealing and process-based method can reduce the cost and improve the modeling efficiency, but the simulated annealing method cannot effectively characterize the long-range connectivity of the pore structure, and the process-based method is not suitable for simulating the rocks that have undergone complex diageneses.In order to solve these problems, a method of reconstructing digital core model based on multiple-point statistics was proposed.Taking as an example the Berea sandstone of the Waverly group in the Michigan Basin, two representative volume elements with the size of 1503 voxels were extracted from different locations of the Berea sandstone.One was set for the training image, the other for the real model.Two orthogonal two-dimensional slices were chosen from the real model as the conditioning data.The 3D digital rock models were reconstructed using multiple-point statistics.The variogram, pore structure parameters and percolation parameters of single phase flow experiment were applied as the evaluation indexes to validate the accuracy of the reconstructed models.The result shows that the reconstructed models agree well with the real model in the pore-throat geometry and topology features, which verifies the accuracy of the reconstructed models and feasibility of reconstructing digital core model using multiple-point statistics.
Keywords: digital core    multiple-point statistics    pore structure    variogram    pore network model    

对于石油和天然气行业而言, 数字岩心技术作为新兴的岩石物理试验数值模拟方法, 已经广泛地应用于地质、地震、测井和开发以及提高采收率等各个领域。与常规岩石物理试验相比, 数字岩心技术能够同时进行成岩作用、固体力学、声学、电学、流体力学及流固耦合等试验的模拟[1-3]。开展该方向的研究对石油和天然气行业勘探和开发的意义不言而喻[4]。一个准确的三维数字岩心模型是进行数字岩心分析技术的前提和基础。随着试验仪器的创新和新理论的突破, 国内外学者不断提出新的构建数字岩心模型的方法。到目前为止, 数字岩心建模方法包括三大类, 即物理试验法、数值重建法和混合法[5]。物理试验法指利用试验仪器(如高倍光学显微镜或X射线CT扫描仪等)对岩心样品拍摄或扫描以获取大量的岩心二维图片, 通过建模程序或软件把二维图片重构成三维数字岩心的方法, 主要有序列切片成像法[6]、激光扫描共聚焦显微镜法[7]和X射线CT扫描法[4]等。数值重建法是指以少量的二维薄片图像为基础, 利用二维图片中包含的信息, 通过随机模拟法或沉积岩形成过程模拟法来重建三维数字岩心的方法, 包括随机模拟法和过程法。随机模拟法大多是基于地质统计学方法, 利用少量的二维薄片信息模拟生成数字岩心模型, 包括截断高斯随机域法[8]、模拟退火法[9-10]、马尔科夫链蒙特卡洛法[11]、序贯指示模拟法[12]、多点统计学法[13]、相位恢复法[14-15]等。过程法是基于模拟沉积岩形成过程(沉积和成岩过程)的一种建模方法[16-17]。混合法是指数值重建法中两种方法的混合, 如Politis等[18]、刘学锋和孙建孟[19]结合模拟退火法和过程法重建了三维数字岩心。上述方法中, 目前常用的有3种方法:X射线CT扫描法、模拟退火法和过程法。虽然X射线CT扫描法是比较准确的建模方法, 但是价格贵试验周期长, 而且不能解决样品的分辨率和尺寸相互矛盾的问题。相比X射线CT扫描法, 模拟退火法和过程法成本低效率高。不过模拟退火法建模是基于传统的两点地质统计学理论, 在传统地质统计学中计算两点地质变量空间之间相关性的重要工具是变差函数, 该函数一般只能表征空间上两点之间的相关性, 很难精确地模拟具有复杂几何形态和拓扑结构的孔隙空间, 这往往会使重构出的孔隙结构缺少长距离连通性[20-21]。过程法模拟沉积岩形成过程时假设颗粒全为球体, 这与实际颗粒形态相差很大, 而且成岩作用模拟过程比较简单, 不适于成岩作用复杂的岩石。为了解决上述问题, 笔者依据多点地质统计学法能够复制地质变量复杂空间结构特征的优点, 对多点地质统计学法构建三维数字岩心模型及重构模型的准确性评价进行研究。

1 方法原理 1.1 多点地质统计学

多点统计学(multiple-point statistics, MPS)最早被Deutsch提出[22]。现今MPS已经广泛应用于储层建模中[23-26], 部分建模理论已经比较成熟。MPS分为迭代法和非迭代法, 其中迭代法包括模拟退火法、基于吉布斯取样迭代法和马尔科夫链蒙特卡洛法。1993年, Guardiano和Srivastava[27]提出了一种非迭代算法, 直接从“密集训练图像”提取局部条件概率, 并利用序贯指示模拟方法进行模拟实现。但该方法在计算过程中需要对每个数据事件进行设置, 而且每模拟一个节点需要重新扫描一遍训练图像, 这严重制约了模拟的效率, 直到2000年Strebelle[27]提出了单一标准方程模拟(single normal equation simulation, SNESIM)算法, 应用搜索树(search tree)一次性存储训练图像中所有数据事件的条件概率分布, 并保证了在模拟过程中能快速地提取条件概率分布函数, 进而大大降低了模拟所需时间。

2004年, Okabe和Blunt[13]把MPS法用于数字岩心建模中, 他们以岩心的二维图片为训练图像, 通过旋转该图片获得xyz方向上的条件数据, 然后利用优选的搜索模板扫描训练图像, 产生一个随机的三维数字岩心模型。将模拟结果与模拟退火法的重构模型对比, MPS有效地解决了重构模型孔隙空间长距离连通性的问题。但是Okabe和Blunt建模的不足之处在于他们假定模型内部孔隙结构在水平方向上和垂直方向上为各向同性。

1.2 SNESIM算法

SNESIM算法是MPS方法中进行离散型变量模拟最常用的方法。通过多重数据模板扫描训练图像来获取地质变量不同尺寸的结构特征, 将扫描得到的所有地质模式存储于搜索树中; 在多个已知节点为条件数据的约束下, 求取未知节点的条件概率分布。再基于序贯模拟样式, 模拟每个未知节点。已经模拟节点的值会变成下一个待模拟节点的条件数据。该算法通过使用多重数据模板和训练图像能够考虑多尺度和多节点之间的关联性, 更加适用于具有长距离连通性的数字岩心模型的重构。为了便于理解该算法构建三维数字岩心的步骤, 下面对算法中的几个重要概念进行简要介绍。

1.2.1 数据模板与数据事件

数据模板(data template)是指由中心节点un个向量{u+ha, a=1, 2, …, n}组成的数据几何体, 通常也称为搜索模板(search template), 常用τn表示。数据模板中心处u的状态是未知值S(u), 其他网格节点处的状态一般为已知值S(u+ha)=S(ua), a=1, 2, …, n。相对应地, 以u为中心, 尺寸大小为n的数据事件(data event)dn是由搜索模板τnn个向量处的数据值{S(uα), α=1, 2, …, n}共同构成。

1.2.2 训练图像

影响数字岩心重构模型准确性的关键因素之一在于训练图像的选择。训练图像本质上就是一个概念模型, 它要尽可能地包含要模拟的地质变量所有的组合模式, 如进行数字岩心建模时的孔隙结构模式。理论上, 选择较大尺寸的训练图像会使建模效果要更好, 因为训练图像尺寸越大, 它所包含的孔喉组合模式也就越多; 但由于受计算机性能的限制, 尺寸较大的训练图像会大大降低建模效率, 因此选择训练图像时要综合考虑多种因素。

对于数字岩心建模, 训练图像可以来源于X射线CT扫描的真实重构模型, 也可以是二维镜下薄片, 如铸体薄片或者扫描电镜图像等。图 1(a)为二维训练图像, 它是由普通薄片二值化后得到的; 其中黑色为颗粒, 白色为孔隙。图 1(b)为用二维7×7的数据模板扫描训练图像的过程。

图 1 数据模板扫描训练图像 Fig.1 Search template scanning training image
1.2.3 条件概率分布函数

对于数字岩心中的某一节点或体素点, 其可能的状态值为{sk, k=1, 2}, s1为颗粒, s2为孔隙。假设在一个随机模式下, 当搜索模板扫描完训练图像后, 某一数据事件dn={S(uα)=skα, α=1, 2, …, n}重复次数为c(dn), 在此条件上, 中心节点u取值{S(u)=sk, k=1, 2}的重复次数为ck(dn), 则中心节点取值S(u)为孔隙或颗粒的概率, 用条件概率分布函数(conditional probability distribution function, CPDF)表示为

$ {\rm{Prob}}\left\{ {S\left( \mathit{\boldsymbol{u}} \right) = {s_k}{\rm{|}}S\left( {{\mathit{\boldsymbol{u}}_\alpha }} \right) = {s_{{k_\alpha }}};\alpha = 1, 2, \ldots, n} \right\} \approx \frac{{{c_k}\left( {{d_n}} \right)}}{{c\left( {{d_n}} \right)}}. $ (1)

当搜索模板扫描训练图像后, 得到的每个数据事件对应的CPDF都会被存储到一个树状数据结构体(搜索树)中。搜索树的构建保证了在模拟过程中搜索模板只需扫描一次训练图像即可把所以数据事件记录下来, 这大大提高了建模的效率。

1.2.4 多重网格

在数字岩心重构过程中, 搜索模板尺寸的选择决定着它能够捕获训练图像中多大尺度的孔隙结构特征。当搜索模板尺寸较小时, 它不能捕获训练图像中大尺度的孔隙结构, 增大数据模板尺寸虽能捕获大尺度的结构信息, 但这往往会影响建模效率。因此拟采用多重网格(multiple grid)来解决这个问题以达到既能捕获大尺度孔隙结构信息又能提高建模效率的效果[28]。多重网格中一个网格逐渐密集化的多重(多级次)搜索模板会被使用, 以代替一个大而密集的搜索模板扫描训练图像[29], 具体过程是先用粗网格的搜索模板扫描训练图像得到大尺度的孔隙结构特征, 把这些模拟结果作为细网格模拟时的硬数据, 再进行细网格搜索模板扫描细网格以捕获小尺度的孔隙结构特征。定义Gg为第g重网格, 在Gg中, 第g重搜索模板τng 的范围会被缩小为最粗搜索模板τn的(2g-1)-1, 即

$ \tau _n^g = \left\{ {{2^{g-1}}{\mathit{\boldsymbol{h}}_\alpha }, \alpha = 1, 2, \ldots, n} \right\}. $ (2)

式中, hα为最粗搜索模板τn中搜索半径, 为矢量。图 2中为一个三重网格(图 2 (a))和一个三重搜索模板(图 2(b)), 其中黑色代表有效节点, 为待模拟节点; 黄色代表已知节点, 即条件数据; 白色为无效节点, 可被忽略。当用第一重搜索模板扫描训练图像, 得到第一重网格中的黑色节点, 获得较大尺度的孔隙结构特征。第一重搜索模板被使用后, 搜索模板尺寸会缩小为第一重搜索模板的“一半”, 若第一重搜索模板的有效节点间距离为4, 第二重搜索模板中有效节点间隔会变为2。在第二重搜索模板扫描训练图像时, 第一重网格上已经模拟过的节点会变成此次模拟时的硬数据, 即图 2(a)第二重网格中的黄色节点。在条件数据约束下, 用图 2(b)第二重搜索模板扫描训练图像模拟得到图 2(a)第二重网格的黑色节点。第三重搜索模板扫描训练图像时, 第二重网格中黑色节点也自动会变成下一个搜索模板扫描训练图像时的条件数据, 即图 2(a)第三重网格所示的黄色节点。最后再用图 2(b)第三重搜索模板扫描训练图像, 可模拟得到第三重网格上剩余黑色节点的值。尽管搜索模板τngτn具有相同的有效节点数, 但比τn具有更大的搜索范围, 因此即使不增加搜索模板的尺寸, 多重搜索模板也能捕获多尺度的孔隙结构特征。

图 2 三重网格和三重搜索模板 Fig.2 Three multiple grids and three multi-grid search templates
2 重构数字岩心模型 2.1 数据预处理

以密西根盆地Waverly组的河流相Berea砂岩为例。Berea砂岩均质性较好, 颗粒以石英为主, 包括少量的长石、白云石和黏土矿物, 颗粒分选较好。Berea砂岩的数据通过X射线微米CT扫描仪采集[30], 大小为4003体素, 分辨率为5.345 μm, 孔隙度为19.6%。从全岩样中提取两个1503体素的体积元为试验对象, 令其中一个体积元为训练图像(图 3(a)), 孔隙度为20.10%;另一个体积元为真实模型(图 3(b)), 孔隙度为19.95%, 图 3中红色为颗粒, 蓝色为孔隙。从真实模型中选取两张正交的二维切片(图 3(c))作为已知信息, 从中提取一部分像素点作为条件数据。

图 3 训练图像、真实模型、条件数据和重构模型 Fig.3 Training image, real model, conditioning data and reconstructed models
2.2 建模步骤

以X射线CT扫描法建立的数字岩心为训练图像, 以两张正交的二维切片的部分像素点为条件数据, 使用多点地质统计学法重构三维数字岩心的具体步骤为:

(1) 选择训练图像。选择X射线CT扫描获得的三维图像作为训练图像(图 3(a)), 三维训练图像包含更加真实的立体的孔隙结构模式。为了捕获较大尺度的孔隙结构, 设定搜索模板的搜索半径为40×40×40, 搜索角度为0°。

(2) 使用四重搜索模板扫描训练图像建立搜索树。四重搜索模板不仅能保证捕获不同尺度的孔隙结构模式, 而且能提高建模效率。

(3) 将条件数据分配到相应的网格, 选定一条随机路径, 访问每一个待模拟节点。

(4) 模拟未知节点u时, 保留那些在最大搜索模板范围内的条件数据, 假设条件数据的数量为n, 相应的数据事件为dn, 在搜索树中检索数据事件dn的CPDF。如果检索过程中数据事件dn的重复数小于设置的最小重复数, 就把搜索模板中最远的条件数据去掉, 在条件数据数量为n-1的条件下, 再去检索数据事件dn-1的CPDF; 如果数据事件的重复数仍小于设定的最小重复数, 继续重复上述操作。倘若条件数据的数量一直减小到n=1时仍不满足要求, 就把目标边缘概率赋值给CPDF。

(5) 在一定的CPDF下, 选择u的一个状态值S(u)作为下一个待模拟节点的条件数据。

(6) 沿着随机路径重复(4)和(5)步骤, 直到模拟完所有未知节点。

训练图像选择时应该优先考虑三维图像, 因为岩心中的孔隙空间实质上是三维展布的, 具有立体的空间特征, 这类似于沉积相中的分流河道具有立体的空间几何特征。二维训练图像往往只是捕获了某一平面的孔隙空间模式特征, 并不能完全反映孔隙空间三维几何和拓扑结构特征如配位数, 而三维训练图像则能提供更加真实的孔隙空间模型。目标边缘概率可以约束随机生成模型的孔隙和颗粒比例, 一般设置该值为目标孔隙度和颗粒比例。还可以通过修改Servosystem系数控制孔隙所占的比例, Servosystem系数取值范围为0~1, 取值越大, 模拟结果中孔隙和颗粒的比例与训练图像中的两者比例相差越小, 但是该系数很高时建模效果反而会变差[31], 建模时设置Servosystem系数为0.9。为了保证获取概率较低的孔隙结构模式, 数据事件重复次数设置为1。按照以上步骤生成3个1503体素(网格)的随机模型(图 3(d)~(f))。

3 重构模型不确定性评价

对于数字岩心分析技术, 数字岩心模型的渗流特性等传输性质本质上取决于孔隙空间的几何结构和拓扑结构特征, 重构模型的可靠性不是指重构模型孔隙空间如孔隙和喉道与真实模型中相同位置上的孔隙和喉道一定要相同, 而是指整个岩心的孔隙空间几何形态和拓扑结构具有等效性, 相同概率的孔隙、喉道和孔喉配置关系出现即可。通过对比重构模型和真实模型的变差函数、孔隙结构参数和渗流特性可评价重构模型的不确定性。

3.1 变差函数

变差函数可以反映地质变量在空间两点之间的相关性, 可以评价孔隙结构在空间中的变化。对真实模型和3个重构模型在xyz方向上求取变差函数, 得到相应的曲线(图 4(a)~(c))。计算结果表明重构模型和真实模型在xyz方向上变差函数曲线十分接近, 3个方向上的变差函数基台值变化范围均为0.15~0.17, 滞后距变化范围为18~20个体素, 这表明重构模型与真实模型在xyz某一方向上的孔隙空间变化具有很好的一致性。

图 4 重构模型与真实模型在xyz方向的变差函数曲线 Fig.4 Variogram curves of reconstructed models and real model in x, y and z directions
3.2 孔隙结构参数

仅依靠变差函数并不能表征孔隙结构复杂的空间变化, 本文中还利用能够表征孔喉几何形态和拓扑结构的孔隙结构参数对重构模型的准确性进一步评估。

孔隙度是表征孔隙空间的基本参数之一。在目标边缘概率和Servosystem系数的控制下, 重构模型的孔隙度与真实模型的相差很小, 最大差值仅为0.46%(表 1)。根据孔隙是否连通, 可把孔隙空间分为连通孔隙和孤立孔隙。孔隙空间的连通情况可用连通系数表示, 连通系数越大, 孔隙连续性越好。对总孔隙空间做连通性测试可以提取连通孔隙部分(图 5), 再用减法运算得到孤立孔隙(图 5)。真实模型中孤立孔隙所占孔隙比例不到3%(表 1), 连通性较好, 随机生成模型孤立孔隙比例略大, 约为6%~8%, 4个模型的总孔隙和连通孔隙体积相差很小(表 1), 表明重构模型也具有较好的长距离连通性。迂曲度是表征孔隙空间复杂程度和弯曲程度的一个参数, 是指沿着某个方向流体质点在孔隙中流经的实际长度与孔隙两端点间距离的比值, 该值越大表示孔隙结构越复杂。重构模型的迂曲度与真实模型的迂曲度比较接近(表 1), 表明重构模型与真实模型孔隙空间的复杂程度类似。

表 1 真实模型和重构模型的孔隙空间参数对比 Table 1 Comparison of pore space parameters for real model and reconstructed models
图 5 真实模型和重建模型的孔隙空间分布对比 Fig.5 Comparison on pore space distribution of real model and reconstructed models

根据所建模型的拓扑性质, 可以把孔隙网络模型分为两大类, 即为规则拓扑结构孔隙网络模型和真实拓扑结构孔隙网络模型。真实拓扑结构孔隙网络模型是基于数字岩心提取与孔隙空间几何特征和拓扑特征等价的孔隙网络。提取孔隙网络模型的方法主要包括多向切片扫描法、孔喉居中轴线法、Voronoi多面体法和最大球法[30]。利用改进的最大球法[32]提取真实模型和重构模型的孔隙网络模型并计算孔隙结构参数包括孔喉半径、孔喉比、孔喉形状因子和配位数。在孔隙结构参数中, 形状因子是描述多孔介质模型孔隙喉道截面形态的重要参数, 也是影响微观渗流的重要因素。形状因子G的计算公式为

$ G = \frac{A}{{{P^2}}}. $ (3)

式中, A为孔隙或喉道的截面面积, μm2; P为孔隙或喉道截面的周长, μm。形状因子G越大, 孔喉截面越趋于圆滑, 流体渗流的阻力越小。

真实模型和重构模型的孔喉半径、孔喉比、孔喉形状因子和配位数的分布见图 6。尽管重构模型中小孔隙和喉道(半径小于10 μm)比例略高于真实模型中的小孔隙比例, 但重构模型和真实模型的孔喉半径分布直方图和孔喉比分布曲线总体上比较接近(图 6(a)~(c))。重构模型孔隙中孤立孔隙比例较大(图 6(d)), 约为13%;4个模型的配位数分布绝大多数集中在1~4, 总体分布一致。重构模型和真实模型孔喉形状因子分布曲线集中区间均为0.02~0.04(图 6(e)~(f)), 4个模型的孔隙和喉道的几何形态都比较接近。

图 6 真实模型和重构模型孔隙结构参数分布 Fig.6 Distribution of pore structure parameters for real model and reconstructed models

孔隙网络模型能直观地展示孔隙空间展布和孔喉配置关系。采用改进后的最大球法提取4个模型的孔隙网络模型(图 7), 图中球体为孔隙, 管柱体为喉道。球体或管柱体越大, 表示此处的孔隙或喉道半径越大。无管柱体连接的球体表示此处的孔隙为孤立孔隙, 孔隙配位数为零。观察孔隙网络模型能够看出4个模型大多数孔隙被2~4个喉道连接, 总体的孔喉配置关系相似。重构模型和真实模型的孔隙结构参数的相似性进一步验证了重构模型的孔隙空间几何特征和拓扑特征与真实模型的相近性。

图 7 真实模型和重构模型孔隙网络模型 Fig.7 Pore network models of real model and reconstructed models
3.3 渗流特性

基于Navier-Stokes方程和达西定律, 以真实和重构数字岩心模型为平台, 进行单相渗流数值模拟试验。渗流过程为稳态流动, 流体为不可压缩的牛顿流体, 流固界面为无滑移壁面。模型的z轴方向一端设置为入口, 另一端为出口, 入口压力为130 000 Pa, 出口为100 000 Pa, 流体黏度设置为0.001 Pa·s。

3.3.1 绝对渗透率对比

根据单相流渗流试验, 计算出真实模型和重构模型z轴方向上的绝对渗透率, 真实模型、重构模型1、重构模型2、重构模型3测得的绝对渗透率分别为1 032.56×10-3、986.96×10-3、993.34×10-3、975.27×10-3 μm2, 该样品试验渗透率为1 100×10-3 μm2。重构模型的绝对渗透率与真实模型模拟出的和试验测得的绝对渗透率值比较接近, 进一步证明了重构模型在渗流特性上的可靠性。重构模型的渗透率略小于真实模型的渗透率, 原因是重构模型中孔隙的平均配位数略低于真实模型的平均配位数。

3.3.2 流速场对比

为进一步对比重构模型和真实模型单相渗流试验渗流的特性, 应用Navier-Stokes方程计算单相流的渗流速度, 并使其可视化。图 8(a)显示了流线在4个模型孔隙空间内部的三维分布, 流线从z轴负半轴流向正半轴。孔隙空间用蓝色表示, 流线流经的部位表示流体在孔隙内部的渗流路径, 流线的颜色深浅代表渗流的速度, 颜色越红速度越大, 颜色偏绿速度偏小。图 8(a)显示4个模型的流线均集中在连通性较好的大孔喉处, 边界处连通性不好的孔喉流线稀疏, 孤立的孔隙内部无流线; 在半径较大配位数较多的孔隙处, 渗流阻力小, 流体流量大, 流速会相对较大, 流线越红; 4个模型大孔喉处的流线颜色均偏红, 表明渗流速度比较接近。

图 8 真实模型和重构模型的孔隙空间内部流速场和压力场分布 Fig.8 Velocity field and pressure field distribution in pore space
3.3.3 压力场对比

压力场分布范围也能反映孔喉分布。在单相流试验中, 从入口端到出口端连通孔喉内的压力逐渐降低, 同一压降下, 压力场分布范围与孔喉中毛细管力有关, 即是与孔喉半径有关。图 8(b)展示了真实模型和重构模型从入口端(z=0处)开始的3个压降变化过程。颜色越浅表示压力越小。同一压降梯度下, 压力场分布范围越大代表孔隙结构中的毛细管力越小, 孔喉半径越大; 4个模型在相同的压降梯度下压力场范围比较接近, 这进一步证明了重构模型的可靠性。

4 结论

(1) 提出了基于多点地质统计学法构建三维数字岩心模型的建模方法, 创建了一套以三维图像为训练图像, 二维切片为条件数据, 利用多点地质统计学法重构数字岩心模型的流程。

(2) 以变差函数曲线、孔隙结构参数和单相流渗流特性为评价参数, 验证了重构模型与真实模型的孔隙空间具有相似的几何结构和拓扑结构特征。

(3) 针对非均质性较弱的岩石, 多点地质统计学能够构建较理想的三维数字岩心模型, 但是该方法建立的多孔介质模型中会出现较大比例的微小孔隙。该方法并不太适用非均质性较强的岩心如包含孔洞的碳酸盐岩, 原因是复杂岩心中某些特殊的孔隙结构模式如碳酸盐岩中的孔洞占总体孔隙结构模式的比例很小, 使得SNESIM算法往往很难把这些概率极低的孔隙结构模式复制到重构模型中。

(4) 重构模型评价中仅使用了变差函数、孔隙结构参数和单相流渗流特性等评价参数, 还可使用其他函数如多点概率分布函数和相对渗透率曲线等进一步验证重建模型的准确性。如果计算机配置允许, 利用多点地质统计学可以重构更大尺度的数字岩心, 再以重构的数字岩心为平台进行渗流试验数值模拟, 这对研究油气在不同尺度上的渗流机制及建立微纳米孔喉尺度的渗流机制与厘米—米岩心或储层尺度的宏观渗流规律之间的联系具有更加重要的意义。

致谢 感谢帝国理工大学Branko Bijeljic高级研究员提供的Berea砂岩数据。

参考文献
[1]
孙建孟, 闫国亮, 姜黎明, 等. 基于数字岩心研究流体性质对裂缝性低渗透储层弹性参数的影响规律[J]. 中国石油大学学报(自然科学版), 2014, 38(3): 39-44.
SUN Jianmeng, YAN Guoliang, JIANG Liming, et al. Research of influence laws of fluid properties on elastic parameters of fractured low permeability reservoir rocks based on digital core[J]. Journal of China University of Petroleum (Edition of Natural Science), 2014, 38(3): 39-44.
[2]
杨永飞, 姚军, 王晨晨. 水湿油藏油气水三相渗流模拟[J]. 中国石油大学学报(自然科学版), 2010, 34(1): 79-83.
YANG Yongfei, YAO Jun, WANG Chenchen. Oil-gas-water three phase flow simulation in water-wet reservoir[J]. Journal of China University of Petroleum (Edition of Natural Science), 2010, 34(1): 79-83.
[3]
邹友龙, 谢然红, 郭江峰, 等. 致密储层数字岩心重构及核磁共振响应模拟[J]. 中国石油大学学报(自然科学版), 2015, 39(6): 63-71.
ZOU Youlong, XIE Ranhong, GUO Jiangfeng, et al. Reconstruction of digital core of tight reservoir and simulation of NMR response[J]. Journal of China University of Petroleum (Edition of Natural Science), 2015, 39(6): 63-71.
[4]
BLUNT M J, BIJELJIC B, DONG H, et al. Pore-scale imaging and modelling[J]. Advances in Water Resources, 2013, 51(3): 197-216.
[5]
林承焰, 吴玉其, 任丽华, 等. 数字岩心建模方法研究现状及展望[J]. 地球物理学进展, 2018, 33(2): 679-689.
LIN Chengyan, WU Yuqi, REN Lihua, et al. Review of digital core modeling methods[J]. Progress in Geophysics, 2018, 33(2): 679-689.
[6]
VOGEL H J, ROTH K. Quantitative morphology and network representation of soil pore structure[J]. Advances in Water Resources, 2001, 24(3/4): 233-242.
[7]
FREDRICH J T. 3D imaging of porous media using laser scanning confocal microscopy with application to microscale transport processes[J]. Physics and Chemistry of the Earth, Part A:Solid Earth and Geodesy, 1999, 24(7): 551-561. DOI:10.1016/S1464-1895(99)00079-4
[8]
QUIBLIER J A. A new three-dimensional modeling technique for studying porous media[J]. Journal of Colloid and Interface Science, 1984, 98(1): 84-102. DOI:10.1016/0021-9797(84)90481-8
[9]
TALUKDAR M S, TORSAETER O, HOWARD J J. Stochastic reconstruction of chalk samples containing vuggy porosity using a conditional simulated annealing technique[J]. Transport in Porous Media, 2004, 57(1): 1-15. DOI:10.1023/B:TIPM.0000032737.59531.cf
[10]
JU Y, HUANG Y, ZHENG J, et al. Multi-thread parallel algorithm for reconstructing 3D large-scale porous structures[J]. Computers & Geosciences, 2017, 101: 10-20.
[11]
WU K, NUNAN N, CRAWFORD J W, et al. An efficient Markov chain model for the simulation of heterogeneous soil structure[J]. Soil Science Society of America Journal, 2004, 68(2): 346-351. DOI:10.2136/sssaj2004.3460
[12]
KEEHM Y. Permeability prediction from thin sections:3D reconstruction and Lattice-Boltzmann flow simulation[J]. Geophysical Research Letters, 2004, 31(4): 1-4.
[13]
OKABE H, BLUNT M J. Prediction of permeability for porous media reconstructed using multiple-point statistics[J]. Physical Review E, Statistical, Nonlinear, and Soft Matter Physics, 2004, 70(6 Pt 2): 66135.
[14]
HASANABADI A, BANIASSADI M, ABRINIA K, et al. 3D microstructural reconstruction of heterogeneous materials from 2D cross sections:a modified phase-recovery algorithm[J]. Computational Materials Science, 2016, 111: 107-115. DOI:10.1016/j.commatsci.2015.09.015
[15]
HASANABADI A, BANIASSADI M, ABRINIA K, et al. Efficient three-phase reconstruction of heterogeneous material from 2D cross-sections via phase-recovery algorithm[J]. Journal of Microscopy, 2016, 264(3): 384-393. DOI:10.1111/jmi.2016.264.issue-3
[16]
ØREN P L, BAKKE S. Process based reconstruction of sandstones and prediction of transport properties[J]. Transport in Porous Media, 2002, 46(2): 311-343.
[17]
ØREN P L, BAKKE S. Reconstruction of Berea sandstone and pore-scale modelling of wettability effects[J]. Journal of Petroleum Science and Engineering, 2003, 39(3/4): 177-199.
[18]
POLITIS M G, KIKKINIDES E S, KAINOURGIAKIS M E, et al. A hybrid process-based and stochastic reconstruction method of porous media[J]. Microporous and Mesoporous Materials, 2008, 110(1): 92-99. DOI:10.1016/j.micromeso.2007.09.024
[19]
LIU X F, SUN J M, WANG H T. Reconstruction of 3-D digital cores using a hybrid method[J]. Applied Geophysics, 2009, 6(2): 105-112. DOI:10.1007/s11770-009-0017-y
[20]
ZHAO X, YAO J, YI Y. A new stochastic method of reconstructing porous media[J]. Transport in Porous Media, 2007, 69(1): 1-11. DOI:10.1007/s11242-006-9052-9
[21]
HAJIZADEH A, SAFEKORDI A, FARHADPOUR F A. A multiple-point statistics algorithm for 3D pore space reconstruction from 2D images[J]. Advances in Water Resources, 2011, 34(10): 1256-1267. DOI:10.1016/j.advwatres.2011.06.003
[22]
DEUTSCH, CLAYTONV. GSLIB geostatistical software library and user's guide[M]. Oxford: Oxford University Press, 1992, 126.
[23]
吴胜和. 储层表征与建模[M]. 北京: 石油工业出版社, 2010, 372-378.
[24]
TAHMASEBI P, HEZARKHANI A, SAHIMI M. Multiple-point geostatistical modeling based on the cross-correlation functions[J]. Computational Geosciences, 2012, 16(3): 779-797. DOI:10.1007/s10596-012-9287-1
[25]
STRAUBHAAR J, WALGENWITZ A, RENARD P. Parallel multiple-point statistics algorithm based on list and tree structures[J]. Mathematical Geosciences, 2013, 45(2): 131-147. DOI:10.1007/s11004-012-9437-y
[26]
张伟, 林承焰, 董春梅. 多点地质统计学在秘鲁D油田地质建模中的应用[J]. 中国石油大学学报(自然科学版), 2008, 32(4): 24-28.
ZHANG Wei, LIN Chengyan, DONG Chunmei. Application of multiple-point geostatistics in geological modeling of D Oilfield in Peru[J]. Journal of China University of Petroleum (Edition of Natural Science), 2008, 32(4): 24-28.
[27]
STREBELLE S. Conditional simulation of complex geological structures using multiple-point statistics[J]. Mathematical Geology, 2002, 34(1): 1-21. DOI:10.1023/A:1014009426274
[28]
REMY N, BOUCHER A, WU J. Applied geostatistics with SGeMS:a user's guide[M]. Cambridge: Cambridge University Press, 2009, 168-205.
[29]
张挺, 李道伦, 卢德唐, 等. 基于多点地质统计法的多孔介质重构研究[J]. 中国科学(G辑:物理学力学天文学), 2009, 39(9): 1348-1360.
ZHANG Ting, LI Daolun, LU Detang, et al. Research on the reconstruction method of porous media using multiple-point geostatistics[J]. Science in China Series G:Physics, Mechanics & Astronomy, 2009, 39(9): 1348-1360.
[30]
DONG H, BLUNT M J. Pore-network extraction from micro-computerized-tomography images[J]. Physical Review E, Statistical, Nonlinear, and Soft Matter Physics, 2009, 80(3 Pt 2): 36307.
[31]
LIU Y. Using the snesim program for multiple-point statistical simulation[J]. Computers & Geosciences, 2006, 32(10): 1544-1563.
[32]
闫国亮. 基于数字岩心储层渗透率模型研究[D]. 青岛: 中国石油大学(华东), 2013.
YAN Guoliang. Research of permeability models of reservoirs based on digital cores[D]. Qingdao: China University of Petroleum (East China), 2013.