基于CT的数字岩心三维建模

林承焰，王 杨，杨 山，任丽华，由春梅，吴松涛，吴玉其，张依旻

1.中国石油大学(华东)地球科学与技术学院，山东 青岛 266580 2.山东省油藏地质重点实验室，山东 青岛 266580 3.中国石油西南油气田分公司勘探开发研究院，成都 610051 4.大庆油田有限责任公司勘探开发研究院，黑龙江 大庆 163000 5.中石油勘探开发研究院，北京 100083

0 引言

在石油工业的推动下，数字岩心技术在油藏工程以及地学领域得到广泛应用和推广，其中包括数字岩心建模及流动模拟技术等。目前国内外对于数字岩心建模方法的研究主要有两种：一是采用X射线CT(computed tomography)扫描岩心，利用图像处理算法将二维图像构建为三维数字岩心；二是基于岩石二维信息(如铸体薄片和岩石粒度分析资料等)，利用建模算法构建三维数字岩心[1]。美国学者Rosenberg等[2]利用CT扫描技术建立了枫丹白露砂岩的三维数字岩心，Arns[3]也通过三维CT扫描得到了孔隙度分别为8%、13%、15%和22%四块枫丹白露砂岩的三维数字岩心，扫描分辨率为5.68 μm/像素。相比于CT扫描建立数字岩心，基于薄片分析的图像建模方法只需要一定数量的岩石切片扫描图像，获取较为方便、经济，该方法被广泛采用[4]。图像建模算法主要包括随机法和过程法，其中随机法主要包括高斯场法、模拟退火法和序贯指示法等，主要是根据数学算法在岩石二维图像统计特征的约束下构建三维数字岩心，使其与岩石二维图像具有相似的统计特征[5]。Joshi[6]首次利用高斯场法建立三维数字岩心。Hazlett[7]首次利用模拟退火法建立三维数字岩心。Bakke等[8]结合岩石的颗粒粒径分布，通过对沉积类岩石地质形成过程(包括沉积、压实和成岩作用)的模拟来构建数字岩心。Arns等[9]通过建立的数字模型计算得到了渗透率等数值，并与实验室测量所得的结果进行比较来验证模型的准确性。刘洋[10]通过对岩心进行切片分析，获取过程法构建数字岩心所需要的相关储层参数，利用过程法的建模算法，实现了岩心的数字可视化。Youssef等[11]应用图像处理算法对模型孔喉进行了两相分割和孔喉参数统计。Ryazanov等[12]运用微观孔喉网络模型模拟了不同润湿性条件下的剩余油类型及其体积含量。刘向君等[13]通过采用COMSOL软件的不可压缩方程模块完成了孔隙空间的微流动模拟。李易霖等[14]通过Avizo建立了致密砂岩多尺度数字岩心模型，并对储层微观孔隙特征进行了定量表征。20世纪80年代中期发展形成的格子玻尔兹曼方法(lattice Boltzmann method, LBM)，来源于完全离散格子气动机算法，目前主要模型有：RK格子波尔兹曼模型、SC模型、自由能方法模型和基于连续性波尔兹曼方程的热力学一致性多相流体理论模型。Pan等[15]首次使用格子玻尔兹曼法和SC模型模拟两相流体在孔隙吼道中的流动，实现了在介质中微观油驱水的数字化过程。

本文依据二维CT资料，通过Avizo软件将CT生成的切片数据进行精确的三维体重构。建模过程中，采用了更为准确的非局部均值滤波算法(NLM filter)进行滤波降噪处理，通过对骨架和孔喉分割，认识储层空间结构，分析孔喉结构特征。流动模拟研究采用等压面的假设，在微观尺度上建立了微观流体流动模拟模型，以期为数字岩心动态模拟提供一种新方法。

1 数字岩心模型建立

CT扫描建立数字岩心主要是依据岩心内部各部分组成不同，其对X射线的吸收强度就不同，因此通过计算机检测射线经过样品吸收后的强度，可以计算得到这个层面上X射线的吸收部分，根据吸收的射线强度就可以转化成对应灰度值以表征不同的结构组分。主要仪器包括X射线源，固定岩样的载物台和收集扫描后射线的X射线检测器(图1)。X射线源向载物台上的物品发出X射线，与物品发生一系列作用后被接收器接收，接收器处理并将其转化为电信号后返回给计算机，记录不同位置的灰度信息。测量样品的密度越大，对X射线的吸收就越多。利用计算机重建可得到三维灰度图像，该图像中的灰度值与被测样品的组分密度值相互对应[16]。

图1 X射线CT扫描原理图Fig.1 Schematic diagram of X-ray computed tomography

1.1 切片获取与模型构建

X射线CT扫描是建立三维数字岩心最常用的方法，但是在实际应用中也存在一些问题[17-18]：由于同一块岩心在不同分辨率下得到的孔隙网络模型基本结构参数和渗流特征有明显差异，如果选取的分辨率过低，就难以识别出岩心中的微细孔喉，使得建模得到的孔隙度较岩心的实际孔隙度小；如果选取的分辨率过高，扫描所得的数据量偏大，因此该方法对计算机等软硬件条件要求较高。本文选取H152工区的细粉砂岩岩样(表1)，根据储层物性特点，选取分辨率为1.49 μm/像素，扫描得到岩样的二维图像，通过对二维切片进行顺序调整、角度矫正和黑点消除，构建三维数字岩心模型(图2)。

1.2 平滑滤波

岩心扫描后，由于噪声等的干扰，会出现灰度值不同的孤立像素点。对于骨架及孔隙来说，即使是相同的组分，也可能会因为密度及含流体性质的细微差异而呈现出不同的灰度值；因此需要采用适当的滤波算法，去除噪声干扰，同时使相同组分灰度值分布更为集中，不同组分灰度值区分更加明显，便于后期分割阈值的选取。不同的滤波算法中，中值滤波为高强度的平滑滤波，但在滤波过程中没有考虑到样品的结构信息，仅仅是对选择区域内所有样品的灰度值进行大小排序，然后选择中间值直接赋值。而非局部均值滤波(NLM filter)不同于中值滤波的原理，其基本思想是：考虑图像的自相似性，当前像素的预测值可由图像中所有与它结构相似的像素加权平均得到[19]。其具体方法为：在选择区域内对比目标和其他部分的结构信息，根据结构上的相似性赋值权重，相似性越高，在最终灰度值的计算中这部分的权重就越大；反之越小。这样就能充分考虑样品结构上的相似性，保留样品的结构信息，便于不同组分的分割。在实际研究工作中，选取的滤波搜索区域为21像素，目标区域为5像素，通过滤波处理得到灰度值分布直方图(图3)和二维横截面灰度图像滤波效果图(图4)。由于中值滤波的高强度平滑作用，因此滤波后的图像高度平滑但是清晰度大大降低(图4b)；而非局部均值滤波处理后，相同组分的单峰形态更加集中，灰度值双峰对比更加明显(图3)，不同灰度值区域的交界处图像有明显的平滑处理，相同组分内部灰度值保真度较高(图4c)，便于后期相态的识别与划分。

表1 CT样品参数表

a.2D切片；b.3D模型。图2 三维数字岩心模型构建Fig.2 Construction of 3D digital model

图3 NLM滤波前(a)后(b)灰度值分布直方图Fig.3 Histogram of gray value distribution before NLM filtering (a) and after (b)

a.原始图像；b.中值滤波图像；c. NLM 滤波图像。图4 NLM滤波和中值滤波效果对比图Fig.4 Contrast figure of NLM filter and median filter

1.3 孔喉分割

微观模型的准确性关键取决于骨架和孔喉的分割。数字岩心建模中孔喉的识别方法主要有两种：一种是根据已知岩样的实际孔隙度大小进行孔喉划分，使划分后的模型孔隙度与实际相同；另一种是依据最大类间距方法进行分割阈值的选取，通过设置不同的阈值，使得组内方差与组间方差的比值最小。本文采用第二种方法，应用Avizo软件进行方差计算和两相分割，分割效果如图5所示。

经过阈值分割后，大部分的孔隙与喉道都可以区分出来，但是由于扫描分辨率与岩样孔喉结构两方面的原因，当相邻的颗粒之间的间隙非常小，即微孔喉尺寸小于扫描分辨率时，部分微小的孔喉被忽略，颗粒被处理成直接接触的形式(图5b)；因此需要应用分离算法对颗粒进行分割并做一定的优化。分离算法是分水岭算法、距离转换算法和数值重建算法的综合体现，通过计算灰度梯度模量对颗粒进行分割，确定颗粒最佳轮廓，显示微小孔喉分布，识别结果如图5c所示。

分离算法的处理可以将岩样中的颗粒分离，得到部分微孔微喉；但是算法的处理过程相对机械，分割产生的孔喉需要进行模型的优化运算。模型的优化主要是对模型进行开运算(opening)和闭运算(closing)。开运算的思路是用运算元沿着图像的边缘进行腐蚀运算，将图像的狭窄部分剔除，然后再用此运算元沿着模型的边缘进行膨胀运算。经过开运算处理后，可以消除小物体、在纤细点处分离物体、平滑较大物体的边界，同时并不明显改变其面积。闭运算的处理与开运算的处理过程相反，即先用运算元进行膨胀运算，然后再用此运算元进行腐蚀边缘，最终可以起到充填细小空间、平滑边缘、连接孤立图像的作用。通过对比可以发现，模型优化后孤立的像素点明显减少，骨架的边缘也明显平滑(图6)。

a.原始灰度图像；b.两相分割图像(白色为岩石骨架，黑色为孔喉，颗粒直接线接触，孔隙不连通)；c.微孔喉识别图像。图5 二值分割效果图Fig.5 Contrast figure of interactive threshold

图6 开运算闭运算前(a)后(b)对比图Fig.6 Contrast figure before opening and closing (a) and after (b)

1.4 体积元提取

由于计算机计算性能的限制，微观特征分析需要在整个模型中选取一定大小的体积元，以此表征整个模型。本文体积元(REV)的选取以孔隙度为约束条件进行。即在模型中选取3个像素点为中心，分别构建棱长为50像素的立方体并计算孔隙度，然后依次以50像素增大立方体的棱长并计算孔隙度，做出3个立方体孔隙度-棱长散点曲线，如图7中3条曲线所示。由图7可知，当立方体棱长大于300像素时，不同位置立方体的自相关函数趋于一致。因此，该样品的REV尺寸为300像素，即数字岩心的体积元棱长为300像素。多次选取体积元进行的计算研究表明，当数字岩心体积元棱长为300像素时，其物理性质几乎不受尺寸的影响，建立REV模型如图8所示。

图7 孔隙度与立方体边长的关系曲线Fig.7 Relation curve of porosity and cube

2 数字岩心模型分析

2.1 静态参数分析

1)颗粒直径

利用标定分析算法可计算模型内部颗粒的体积分布参数V3d，然后得到等效直径公式：

(1)

利用公式(1)可计算出颗粒等效直径分布情况(表2)，并对每个颗粒的直径进行统计，建立颗粒等效直径分布直方图(图9)。

根据表2和图9可知，粒径主要分布区间为10～70像素。已知扫描该样品所用分辨率为1.49 μm/像素，则数字岩心模型粒径主要分布在14.9～104.3 μm；与研究区细粉砂岩的粒度主峰基本一致。数字岩心模型中最大颗粒直径为129.2像素，约190 μm，这是由于极少部分的微细孔喉未能分割，造成颗粒之间为直接接触，导致粒径明显偏大。

2)孔隙度

统计阈值分割后模型中骨架部分和孔隙部分的像素体积，结果如表3所示。孔隙度为孔隙部分的像素在整个模型中所占的体积比，约为15.54%。而研究区实际孔隙度为22.5%。计算孔隙度明显小于实验室测试值，其主要原因是部分微孔喉尺寸明显小于分辨率1.49 μm，因此未能在切片上明显识别，最终导致孔隙度偏小。

a.体积元渲染(绿色边框为镂空部分边界)；b.体积元正交视图(白色为岩石骨架，黑色为孔喉)。图8 体积元三维渲染及正交视图Fig.8 Volume rendering and ortho view of REV

均值/像素最小值/像素最大值/像素中值/像素55.034.66129.2034.50

图9 颗粒等效直径分布直方图Fig.9 Histogram of diameters of labels’ Deq

模型总体积/(106μm3)骨架体积/(106μm3)孔喉体积/(106μm3)孔隙度/%8.947.551.3915.54

3)孔喉直径分布

选取模型内部的孔喉，计算得到孔隙体积V3d的分布情况，依据公式(1)计算内径大小，得到其孔喉直径分布直方图(图10)。研究区孔喉粒径较小，集中分布在20像素以下，大孔发育较少，储层以中孔低渗为特征。

图10 孔喉等效直径分布直方图Fig.10 Histogram of diameters of pores’ Deq

4)孔喉连通性与曲折度

孔喉连通性和孔喉曲折度是影响流体流动的重要参数，一般来说连通性越好，曲折度越低，流体的流动越容易。曲折度算法原理为计算孔隙边界在某一轴向上不同切面之间的质点距离，z(n)面到z(0)面距离叠加之和与H的比值即为曲折度(图11)。连通性计算是依赖于渗流的一种算法，通过沿着某一轴向进行渗流，得到连通的体积与全体积的比值，即为岩样的连通性。

图11 曲折度计算原理Fig.11 Principle of centroid path of tortuosity

经计算：模型的孔喉连通性为92.16%，连通性较好；曲折度为2.62，即孔喉边界的曲线长度是对应节点间直线长度的2.62倍，孔喉发育相对复杂，由此产生的流动阻力较大。根据模型连通性测试结果，依据数学逻辑算法，可得到模型内部的孤立孔喉分布情况，结果如图12所示。

5)孔喉网络骨架提取

根据燃烧算法，用孔隙中心与喉道中轴线分别代表孔隙和喉道，将复杂不规则的孔喉空间用相互连接的中轴线骨架模型来表征(图13a)。尽管中轴线骨架无法完整地刻画孔隙类型和尺寸，但是它准确地描述了孔喉网络的结构特征[20]。根据孔喉直径计算结果，将等效直径值赋值到骨架模型中，可得到模型的孔喉结构填充模型(图13b)，该模型中不同颜色展示了模型内部孔喉发育特征：模型内部多发育小孔细喉，很少存在大的连通孔隙。图13c展示了某一切面上灰度值图像和孔喉的对应特征。根据图13b、c可得出：模型内部孔喉不均一性较为突出，孔喉大小分布不均，孔喉结构非均质性明显。

a.全孔喉；b.连通孔喉；c.孤立孔喉。图12 孔喉连通性测试结果Fig.12 Results of connectivity

a.孔喉网络骨架模型；b.孔喉结构填充模型；c.切片对比。图13 孔喉模型Fig.13 Pore-throat model

多孔介质中孔隙和喉道的连接关系反映了多孔介质内部的连通性。一个孔隙所连接的喉道数量称为配位数，配位数值的大小反映了网络模型的连通性好坏，同时配位数也是影响岩石渗透率的重要因素[21]。因此，根据建立的孔喉网络骨架模型可分析模型的孔喉配位数并建立频率分布直方图，结果如图14所示。由图14可知：模型中存在一定数量的配位数为1的孔隙，剩余油在这些孔隙中富集；除此之外，配位数主要为2～6，相对较低，孔喉连通性较差，模型渗透性较低，研究区储层具有低渗透特征。

图14 配位数频率分布直方图Fig.14 Frequency distribution histogram of coordination number

2.2 流动特征分析

1)绝对渗透率计算

根据建立的模型，利用软件的有限元模拟模块分别在x,y,z3个方向计算模型的渗透率数值，计算时设置进口压力为130 kPa，出口压力为100 kPa，流体黏度为0.001 Pa·s，得到x,y,z方向的渗透率分别为：33.7×10-3、41.7×10-3、23.8×10-3μm2。渗透率较小，属于低渗透储层，且不同方向渗透率值略有差异。实验室测量该井段的水平渗透率值为19.28×10-3μm2，由于部分连通微孔喉未能识别出来，渗流阻力明显变小，计算值略大于实验室测量值。

2)速度场流线表征

根据计算的渗透率值，以流线的形式表征相态在模型内部的流动路径，左端为入口，右端为出口(图15)。根据流线分布情况，模型内部主流线区域为大孔喉连通部分，流线较为密集，模型角隅部分流线稀疏，存在大量的未波及区，是剩余油潜在富集区域。

3)流动模拟

在多孔介质流动中，当介质内部彼此连通、渗流过程发生得非常缓慢时，即在准静止条件下，每一个等压面上其液-液界面和液-固界面都是平衡的，在不考虑重力和动力情况下，起作用的力就是毛管压力[22]。这样，在任一时刻，流体内部的压力会彼此传递，液面最前端压力值保持相等，即在理想条件下流体液面形态、移动方向与等压面保持一致。因此，压降与相态流动液面等效，即液面前缘为等压面。通过等压面的变化来表征模型中相态的流动，并根据压降的变化来计算流体饱和度大小。根据Stokes方程，设定流体为不可压缩的牛顿流体，液固边界无滑移，在准静止条件下模拟层流渗流过程，结果如图16所示。

a. x轴；b. y轴；c. z轴。蓝色表示模型中连通的孔喉空间，红色流线表示速度场分布。图15 速度场流线表征Fig.15 Streamlines of velocity in the x, y, z axis

a. Sw=0；b. Sw=16.7%；c. Sw=78%。Sw为含水饱和度。红色表示孔喉空间，蓝色表示水相。图16 流动模拟Fig.16 Percolation simulation

由图16可知，在一端注水的情况下：由于孔喉大小的不同，液面会受到毛管压力的作用，大孔喉连通区域液面出现突进现象，在三维上呈现为不平行的液面特征(图16b)；随着持续注水，水占据大部分的孔喉空间，但是部分微孔喉和孤立孔喉未被波及(图16c)。该方法可简化模型内部流动过程，初步建立毛管压力作用下的渗流模型，模型中的未波及区域为剩余油潜在富集区域。

3 结论与认识

1)通过对优选的H152样品进行X射线CT扫描，可以构建三维数字岩心模型，能够模拟孔喉网络结构，建立了微观尺度储层模型，定量表征了粒径、孔喉半径、配位数、孔隙度、渗透率等动静态参数及低渗透储层特征。

2)结合研究区储层样品孔渗实际资料，对所建数字岩心模型进行验证和校对，能够为储层流动模拟提供可靠的地质模型。

3)基于等压面的微观渗流模拟初步揭示了低渗透储层流体渗流特征，提供了一种流动模拟新方法，但是由于只考虑了毛管压力因素，还需进一步分析孔喉的润湿性、流体重力等其他多种参数的影响。

[1] 刘学锋. 基于数字岩心的岩石声电特性微观数值模拟研究[D].青岛：中国石油大学(华东),2010.

Liu Xuefeng. Numerical Simulation of Elastic and Electrical Properties of Rock Based on Digital Cores [D]. Qingdao: China University of Petroleum (East China), 2010.

[2] Rosenberg E, Lynch J, Guéroult P, et al. High Re-solution 3D Reconstructions of Rocks and Composites[J]. Oil & Gas Science & Technology, 1999, 54(4):497-511.

[3] Arns C. H. The Influence of Morphology on Physical Properties of Reservoir Rocks [D]. Sydney: The University of New South Wales, 2002.

[4] 姚军,赵秀才,衣艳静,等. 数字岩心技术现状及展望[J]. 油气地质与采收率,2005,12(6)：52-54.

Yao Jun, Zhao Xiucai, Yi Yanjing, et al. Status and Prospect on Digital Core Technology [J]. Petroleum Geology and Recovery Efficiency, 2005,12(6):52-54.

[5] 谢一婷. 基于数字岩心技术的低渗油藏渗流机理及合理井距研究[D].成都:西南石油大学,2013.

Xie Yiting. The Research of Seepage Flow in Low Permeability Reservoir and the Reasonable Well Spacing Based on Digital Core Technology [D]. Chengdu: Southwest Petroleum University, 2013.

[6]Joshi M. A Class of Stochastic Models for Porous Media[D]. Lawrence Kansas: University of Kansas, 1974.

[7] Hazlett R D. Statistical Characterization and Stochastic Modeling of Pore Networks in Relation to Fluid Flow [J]. Mathematical Geosciences, 1997, 29(6):801-822.

[8]Bakke S, Ren P E. 3-D Pore-Scale Modelling of Sandstones and Flow Simulations in the Pore Networks[J]. Spe Journal, 1997, 2(2):136-149.

[9] Arns C H, Averdunk H, Bauget F, et al. Digital Core Laboratory: Reservoir Core Analysis from 3D Images[C]// The 6th North America Rock Mechanics Symposium (NARMS). Houston: American Rock Mechanics Association, 2004: ARMA-04-498.

[10] 刘洋. 过程法构建数字岩心技术[D].青岛：中国石油大学 (华东),2007.

Liu Yang. Construction of Digital Core by Process-Based Simulation Method [D]. Qingdao: China University of Petroleum (East China), 2007.

[11]Youssef S, Rosenberg E, Gland N F, et al. High Resolution CT And Pore-Network Models To Assess Petrophysical Properties Of Homogeneous And Heterogeneous Carbonates[C]// Reservoir Characterization and Simulation Conference. Abu Dhabi: Society of Petroleum Engineers. doi：10.2118/111427-MS.

[12] Ryazanov A V, Sorbie K S, Dijke M I J V. Structure of Residual Oil as a Function of Wettability Using Pore-Network Modelling[J]. Advances in Water Resources, 2014, 63(2):11-21.

[13] 刘向君,朱洪林,梁利喜. 基于微CT技术的砂岩数字岩石物理实验[J]. 地球物理学报, 2014, 57(4):1133-1140.

Liu Xiangjun, Zhu Honglin, Liang Lixi. Digital Rock Physics of Sandstone Based on Micro-CT Technology [J]. Chinese Journal of Geophysics, 2014, 57(4): 1133-1140.

[14] 李易霖, 张云峰, 丛琳,等. X-CT 扫描成像技术在致密砂岩微观孔隙结构表征中的应用：以大安油田扶余油层为例[J]. 吉林大学学报(地球科学版), 2016, 46(2):379-387.

Li Yilin, Zhang Yunfeng, Cong Lin, et al. Application of X-CT Scanning Technique in the Characterization of Micro Pore Structure of Tight Sandstone Reservoir: An Example from Fuyu Oil Layer in Daan Oil Filed[J]. Journal of Jilin University (Earth Science Edition), 2016, 46(2): 379-387.

[15] Pan C, Hilpert M, Miller C T. Lattice-Boltzmann Simulation of Two-Phase Flow in Porous Media[J]. Water Resources Research, 2004, 40(1):62-74.

[16] 屈乐. 基于低渗透储层的三维数字岩心建模及应用[D].西安:西北大学,2014.

Qu Le. Modeling and Application of Three-Dimensional Digital Cores of Low Permeability Reservoir [D]. Xi’an: Northwestern University, 2014.

[17] 刘学锋,张伟伟,孙建孟. 三维数字岩心建模方法综述[J]. 地球物理学进展,2013,28(6):3066-3072.

Liu Xuefeng, Zhang Weiwei, Sun Jianmeng. Methods of Constructing 3-D Digital Cores: A Review [J]. Progress in Geophysical, 2013,28(6):3066-3072.

[18] 王晨晨,姚军,杨永飞,等. 基于CT扫描法构建数字岩心的分辨率选取研究[J]. 科学技术与工程,2013,13(4):1049-1052.

Wang Chenchen, Yao Jun, Yang Yongfei, et al. The Research of Selected Resolution Based on CT Scanning Method for Constructing Digital Core[J]. Science, Technology and Engineering, 2013,13(4):1049-1052.

[19] 张爽,周慧鑫,牛肖雪,等. 基于非局部均值滤波与时域高通滤波的非均匀性校正算法[J]. 光子学报,2014,43(1):147-150.

Zhang Shuang, Zhou Huixin, Niu Xiaoxue, et al. Temporal High-Pass Filter Nonunifority Correction Algorithm Based on Non-Local Means Filter for Infrared Focal Plane Array[J]. Acta Photonica Sinica, 2014, 43(1): 147-150.

[20] 杨山. 基于驱替实验及数字岩心的微观剩余油研究[D].青岛:中国石油大学(华东),2015.

Yang Shan.The Research of Microscopic Remaining Oil on the Basis of Displacement Experiment and Digital Core[D]. Qingdao: China University of Petroleum (East China), 2015.

[21] 任奕明. 微观水驱油网络模拟研究[D].成都:西南石油大学,2014.

Ren Yiming. The Network Simulation of Microscopic Water Flooding[D]. Chengdu: Southwest Petroleum University, 2014.

[22] Mohammadmoradi P, Kantzas A. Pore-Scale Perme-ability Calculation Using CFD and DSMC Techniques[J]. Journal of Petroleum Science & Engineering, 2016, 146:515-525.