基于数字岩心电性数值模拟新方法的研究

孔强夫，胡 松，王晓畅，李 军，熊培祺，王 琬，刘 秘

(1.中国石化石油勘探开发研究院，北京 100083；2.中国石油大学，北京 102249；3.加拿大卡尔加里大学，加拿大阿尔伯塔，T2N4V7)



基于数字岩心电性数值模拟新方法的研究

孔强夫1，胡 松1，王晓畅1，李 军1，熊培祺2，王 琬3，刘 秘2

(1.中国石化石油勘探开发研究院，北京 100083；2.中国石油大学，北京 102249；3.加拿大卡尔加里大学，加拿大阿尔伯塔，T2N4V7)

基于数字岩心的岩石电性数值模拟技术可以弥补传统岩石物理实验的诸多不足，而准确厘定孔隙空间流体分布是电性数值模拟的基础。在对常用模拟方法进行了系统阐述的基础上，分析了利用孔隙网络模型、数学形态学及与其他模拟方法结合开展电性数值模拟研究的优缺点，同时提出了将数学形态学与随机游走结合的新方法，并在图像分析基础上提出利用“等效”处理办法表征水膜。在考虑了数字岩心的连通性、各向异性、水膜等因素的基础上，对高孔、高渗砂岩岩心进行了数字模拟，模拟结果与理论计算结果吻合较好，证实了新方法的可行性，为进一步开展电性数值模拟研究奠定了基础。

数字岩心；电性数值模拟；数学形态学；随机游走

基于数字岩心的岩石电性数值模拟研究引起了国内外学者的广泛关注。澳大利亚国立大学、帝国理工大学、法国石油研究院和斯坦福大学自21世纪初开展了相应的研究工作，中国石油大学也开展了相应的研究工作[1-5]，但尚处于起步阶段。

岩石电性数值模拟的首要任务就是准确确定孔隙空间流体的分布，进而计算不同含水饱和度下岩石的电性参数。模拟不同含水饱和度下孔隙流体分布的方法主要有孔隙网络模型和数学形态学方法[5-13]。孔隙网络模型不仅可以方便研究流体在孔隙中的分布规律，计算流体的传导性能，而且还可以克服计算时间长的缺点。但该方法对孔隙空间进行简化，势必造成岩心真实孔隙信息的丢失，同时难以处理复杂的孔隙空间(如裂缝等)。数学形态学基于岩心真实图像，采用一般为球形的结构元素对孔隙空间做开运算来模拟不同驱替压力下的流体分布，这种方法能够保留岩心的真实孔隙信息，使岩石电性的数值模拟结果更为准确。

确定孔隙空间流体分布后，采用某种数值模拟方法开展电性模拟研究，这些方法主要有格子玻尔兹曼方法、有限元法、随机游走方法、基尔霍夫电路节点法[4-5，20-33]。格子玻尔兹曼方法是基于流体微观模型和细观动理论方程的方法，该方法易于处理复杂的边界条件，具有完全并行性，但对于多相流的模型需进一步完善，且速度分量多时模拟速度变慢；有限元法是利用分块逼近，根据三维数字图像中能量最小原理，确定每一个节点上的电压分布，是电性数值模拟较为经典的方法，但该方法求解比较困难，对于网格的剖分具有一定的随机性；基尔霍夫电路节点法利用基尔霍夫定律在提取孔隙格架的基础上求解每一个节点的电阻率，原理简单，容易实现，但该方法本身要在孔隙网络模型的基础上开展模拟，结果会存在一定的误差，另外运算量非常庞大；随机游走方法则是利用稳定状态下扩散和电流传导之间解相同的拉普拉斯方程，通过粒子的随机游走计算孔隙空间的迂曲度，进而利用迂曲度来计算岩石的电性参数，该方法原理简单易于实现。

本文在前人研究基础上将数学形态学算法与随机游走方法结合进行岩石电性参数的模拟尝试，通过分析数字岩心的连通性，采用“等效”处理方法表征水膜的存在，计算了高孔、高渗砂岩三维数字岩心在岩心轴向上(z方向)的电性参数。岩石地层因素F与实验结果吻合；在整个含水饱和度区间内岩石电阻率指数与实验结果吻合，证实了新方法的可行性，为进一步开展岩石电性数字模拟研究奠定了坚实的基础。

1 方法原理

1.1 数学形态学原理

数学形态学是由一组形态学的代数运算组成的，包括4个基本运算：膨胀(或扩张)、(腐蚀)或侵蚀、开运算和闭运算。膨胀运算和腐蚀运算是数学形态学处理的基础，许多数学形态学算法都是以这两种基本运算作为基础的。

设C和D是中Z2的集合，Φ为空集，C被D腐蚀，记为腐蚀算子⊖。腐蚀运算定义为：

C⊖D={x|(D)x⊆C≠Φ}

(1)

式中x——集合中的任意一点 。

式(1)说明集合C被集合D的腐蚀结果是所有使D被x平移后包含于C的点x的集合。图1表示了一个腐蚀过程，图中阴影部分为腐蚀之后的结果。

腐蚀运算的作用是去除图像中的孤立点和毛刺，缩减区域边界，增大内部孔洞，断开某些区域间的联系。

设C和D是Z2中的集合，Φ为空集，∩为两个集合的交集，x为集合中的任意元素，C被D的膨胀，记为膨胀算子⨁。膨胀运算的定义为：

C⨁D={x|(D)x∩C≠Φ}

(2)

式(2)说明膨胀运算过程是集合D首先做关于其原点的映射，然后平移x。C被D的膨胀运算就是D被所有的x平移后与集合C至少存在一个非零公共元素。

与其他数学形态学算法一样，集合D在膨胀运算中被称为结构元素。图2表示一个膨胀过程，图中各参数含义与图1相同；图中阴影部分为膨胀之后的结果。

膨胀运算的作用是扩大图像区域的边界，减小区域内的孔洞。

取结构元素C作开运算，记为C·D,定义为：

C·D=(C⊖D)⨁D

(3)

即集合C被集合D开运算就是C被D腐蚀后的结果再被D膨胀。

数学形态学中开运算有一个简单的几何解释，假设把圆盘结构元素D看作一个平面内的“滚动球”，C·D的边界为集合D在C中滚动所能达到的最远处的D的边界所构成。集合C中所有朝外的突出部分保持不变，突出且不能容下结构元素的部分被去掉(图3)。

开运算能够在保持总位置和形状不变条件下消除比结构元素小的特定的图像细节、平滑图像的轮廓、扩大缝隙、分离由细线相连的物体、消除细小的突起和孤立点。

取结构元素D作闭运算，记为C·D,定义为：

C·D=(C⨁D)⊖D

(4)

即集合C被集合D闭运算就是C被D膨胀后的结果再被D腐蚀。

闭运算解释为在边界外滚动球，所有朝内的突出部分均被圆滑了，而朝外的则没有影响。集合C最左边的凹入部分被大幅减弱了(图4)。

闭运算能在保持总位置和形状不变条件下填充目标中的孔洞，连接相邻目标和小的断点，从外部平滑边界、填充图像中的缝隙、融合相邻的图像小块。

1.2 随机游走模拟电学参数原理

(5)

式中 τ——相应的迂曲度；

t——扩散时间；

D——溶质的扩散系数；

Dbulk——溶液的扩散系数。

式中比值可以通过可溶物质的均方位移与时间的比值在时间足够长时取极限来表示，即：

(6)

式中 r(t)——t时刻溶质的位置；

r(0) ——初始时刻溶质的位置。

(7)

式中 τe——电流迂曲度；

σ——饱和流体时物质的电导率；

σbulk——导电流体的电导率。

Clennell(1997)[30]详细论述了扩散迂曲度与电流迂曲度之间存在一个修正系数1/φ,即：

(8)

式中 F——地层因素；

ρ(t→∞)——t很大时饱和流体物质的电阻率；

ρbulk——导电流体的电阻率；

φ——连通孔隙度；

τ——含水饱和度100%时孔隙空间的迂曲度。

如果考虑导电的水的有效扩散系数Dw，则F可以表示为：

(9)

后面扩散项的比值可通过随机步行者的布朗运动模拟水分子的扩散来计算。FR即为阿尔奇公式中的地层因素F，是孔隙度与迂曲度的函数,描述纯净的、无泥质且100%含水的砂岩的电阻率与孔隙水的电阻率成正比。

同样地，电阻率指数RI，其描述的是对含水饱和度小于1的纯砂岩而言，其电阻率与同种砂岩在100%含水时的电阻率成正比，是含水饱和度与迂曲度的函数，即：

(10)

式中 τ′——含水饱和度小于100%时孔隙空间的迂曲度；

τ——含水饱和度100%时孔隙空间的迂曲度；

RI——电阻增大系数；

Sw=100%——含水饱和度(100%)。

2 三维孔隙空间流体分布表征

图5a为某砂岩X射线CT扫描的原始二维图像，尺寸为600×600像素，其中黑色代表岩石骨架，白色代表岩石孔隙。原始图像中用0表示岩石骨架，1表示孔隙空间。选取半径R为15个像素的圆作为结构元素，对图5a中的孔隙空间(白色区域)分别进行腐蚀、膨胀和开运算，运算结果如图5b—d所示，红色部分表示孔隙空间经过相应运算后的结果。从图中可以看出膨胀运算扩大图像，腐蚀运算收缩图像。开运算的另一种解释为结构元素B在X内滚动所能达到的边界构成的空间，计算结果如图5d中红色区域所示。从开运算的定义可以看出对孔隙空间做开运算其结果为保留所有大于结构元素半径R的孔隙空间。

在水润湿岩石中，非润湿的油相首先占据孔隙空间中的大孔隙，随着驱替压力的增大，油按照孔隙半径由大到小依次侵入。对于油润湿岩石，可采用类似方法，结构元素半径由大到小依次变化，开运算的结果表征的孔隙空间被地层水占据，其余空间则被油占据。假设孔隙空间开运算的结果表征油驱水过程中的油占据的空间，其余孔隙空间表征地层水占据的空间，则该过程与水润湿岩石的排驱过程相似。因此岩石孔隙空间的开运算可以模拟水润湿岩石的排驱过程，进而确定在不同含水饱和度下孔隙空间中油和水的分布。图6显示的是不同含水饱和度下孔隙流体的分布，其中蓝色表示骨架，绿色表示水，红色表示油，可以看出随着结构元素逐渐减小，含水饱和度逐渐降低。

3 方法验证

在对数字岩心进行连通性分析时，通常是沿着某个方向进行分析(图7)，对于正方形数据体的数字岩心，x、y、z方向连通孔隙的提取过程如下。

(1)对三维数据体进行26邻域扫描，在扫描过程中对每一个独立的连通区域进行标记，按次序标记为1，2，3，…，n，则此三维数据体有n个连通域。

(2)判断x方向上两个边界面是否有相同的标记符，即1，2，3，…，n，假如存在相同的标记符i，j，…，m(1≤i，j，m ≤n)，则此数字岩心在x方向具有连通的孔隙。将标记为i，j，…，m所代表的子数据体存放到与数据体D相同大小的三维矩阵Dx中。

(3)在y、z方向进行同样操作，提取出Dy、Dz。

如图7所示，对三维数据体扫描后标记有1、2、3、4共4个独立的连通域，在x方向的前后两个面都有标记为3的像素点，则Dx就只有标记3的这个连通域被记录下来，Dy只包含2，Dz只包含1。标记为4的连通域不参与后续模拟计算。为了验证连通性算法分析的有效性，以NO2岩心的300×300×300数据体为例，孔隙连通性分析如图8所示。

本次研究共选取5块南堡油田高孔、高渗砂岩岩心(NO1、NO2、NO15、NO17、NO19)做X射线CT扫描构建三维数字岩心，其图像扫描分辨率均为12μm,分割像素尺寸均为300×300×300。在分析数字岩心连通性的基础上利用随机游走法公式(5)、公式(6)分别模拟计算了这5块岩心在完全饱和水状态下3个轴向上的连通孔隙度和迂曲度，结果如表1所示。

表1 5块岩心孔隙度与迂曲度计算结果表

在利用数学形态学方法确定不同含水饱和度下孔隙空间油、水的分布以后，对于水润湿岩石，孔隙空间被油侵入以后，会在孔隙壁上存在一层水膜，这些水膜的存在为电流的传导提供了一条附加路径。在中、高含水饱和度时这种传导作用并不是十分明显；但是在低含水饱和度时，这种传导作用变的十分重要。因此，要在整个三维数字岩心含水饱和度范围内准确模拟岩石的电性参数，水膜的传导作用必须考虑。但由于一般实验测量的水膜厚度有限，在现有的CT扫描分辨率下无法被识别，因此不能直接将油侵孔隙的最外面一层像素视为水膜，在图像分析的基础上采用一种“等效”处理方法，将油相的最外一层像素与水膜的厚度相乘累加之后在油相与骨架的接触位置设置水膜。水膜的表征图如图10所示。本次研究水膜的厚度参考向阳[34]实验室的测量结果，给定的水膜厚度为0.057μm，用0.6μm的岩心扫描图像分辨率得到含水饱和度为26.9%。

4 结 论

(1)在图像分析基础上提出利用“等效”处理办法表征水膜。

(2)数学形态学方法在保证图像真实孔隙信息方面具有一定优势，随机游走模拟电性参数原理简单，易于实现，本文利用两种方法的结合进行了电性数值模拟的探索。

(3)对高孔、高渗岩心的模拟结果与理论计算结果基本吻合，证实了新方法的可行性，为进一步开展电性数字模拟研究奠定了坚实的基础。

[1] 陶果, 岳文正, 谢然红, 等. 岩石物理的理论模拟和数值实验新方法[J]. 地球物理学进展, 2005, 20(1): 4-11.

[2] 朱益华, 陶果. 顺序指示模拟技术及其在 3D 数字岩心建模中的应用[J].测井技术,2007, 31(2): 112-115.

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

[4] 姚军, 赵秀才, 衣艳静, 等. 储层岩石微观结构性质的分析方法[J]. 中国石油大学学报 (自然科学版), 2007, 31(1): 80-86.

[5] 刘学峰.基于数字岩心的岩石声电特性微观数值模拟研究[D].东营：中国石油大学(华东)，2010.

[6] Wang Y, Sharma M M. A network model for the resistivity behavior of partially saturated rocks[C].SPWLA 29th Annual Logging Symposium. Society of Petrophysicists and Well-Log Analysts, 1988.

[7] 高楚桥, 毛志强. 润湿性对岩石电性的影响[J]. 地球物理学进展, 1998, 13(1): 60-72.

[8] 陈家军, 杨建, 田亮. 基于孔隙网络模型的非水溶相液体运移实验研究进展[J]. 地球科学进展, 2007, 22(10): 997-1004.

[9] 陈民锋, 姜汉桥. 基于孔隙网络模型的微观水驱油驱替特征变化规律研究[J]. 石油天然气学报, 2007, 28(5): 91-95.

[10] Suman R J, Knight R J. Effects of pore structure and wettability on the electrical resistivity of partially saturated rocks-A network study[J]. Geophysics, 1997, 62(4): 1151-1162.

[11] Hilpert M, Miller C T. Pore-morphology-based simulation of drainage in totally wetting porous media[J]. Advances in Water Resources, 2001, 24(3): 243-255.

[12] 唐常青. 数学形态学方法及其应用[M]. 科学出版社, 1990.

[13] 李学民，曹俊兴，何晓燕. 用格子玻尔兹曼方法模拟非均匀介质中的电场响应[J]. 地球物理学报, 2004，47(2): 349-353.

[14] Xu You Sheng，Li Hua Bing，Fang Hai Ping，etal. Lattice Boltzmann simulation for nonlinear flow in porous media with coupling reaction[J]. Acta Physica Sinica, 2004, 53(3):773-777.

[15] 岳文正，陶果，朱克勤. 饱和多相流体岩石电性的格子气模拟[J].地球物理学报，2004，4(5):905-910.

[16] 张武生，杨燕华，徐济鋆.格子波尔兹曼方法及其应用[J]. 现代机械, 2003 (4): 4-6.

[17] Martys N S, Chen H. Simulation of multicomponent fluids in complex three-dimensional geometries by the lattice Boltzmann method.[J]. Physical Review E Statistical Physics Plasmas Fluids & Related Interdisciplinary Topics, 1996, 53(1):743-750.

[18] 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.

[19] Mora P. The lattice Boltzmann phononic lattice solid[J]. Journal of statistical physics, 1992, 68(3-4): 591-609.

[20] Adler P M, Jacquin C G, Thovert J F.The formation factor of reconstructed porous media[J]. Water Resources Research, 1992，28(6): 1571-1576.

[21] Amabeoku M O, Al-Ghamdi T M, Toelke J,etal.Evaluation and Application of Digital Rock Physics (DRP) for Special Core Analysis in Carbonate Formations[C].International Petroleum Technology Conference.2013.

[22] Arns C H, Knackstedt M A, Pinczewski M V,etal. Accurate estimation of transport properties from microtomographic images[J]. Geophysical Research Letters, 2001，28(17): 3361-3364.

[23] Arns C H, Knackstedt M A, Pinczewski M V,etal. Accurate estimation of transport properties from microtomographic images[J]. Geophysical Research Letters, 2001，28(17): 3361-3364.

[24] Martys N, Garboczi E J. Length scales relating the fluid permeability and electrical conductivity in random two-dimensional model porous media[J]. Physical Review B, 1992，46(10): 60-80.

[25] Adler P M, Jacquin C G, Thovert J F. The formation factor of reconstructed porous media[J]. Water Resources Research, 1992，28(6): 1571-1576.

[26] Toumelin E. Pore-scale petrophysical models for the simulation and combined interpretation of nuclear magnetic resonance and wide-band electromagnetic measurements of saturated rocks[J]. Dissertations & Theses - Gradworks, 2006, DAI/B 67-12.

[27] Jin G, Torres-Verdin C, Devarajan S,etal. Pore-scale analysis of the Waxman-Smits shaly-sand conductivity model[J]. Petrophysics, 2007，48(2): 104-120.

[28] Latour L L, Mitra P P, Kleinberg R L,etal. Time-dependent diffusion coefficient of fluids in porous media as a probe of surface-to-volume ratio[J]. Journal of Magnetic Resonance, Series A, 1993，101(3): 342-346.

[29] Rasmus J C. Summary of the effects of various pore geometries and their wettabilies on measured and in-situ values of cementation and saturation exponents[J]. Log Analyst, 1986, 28(1):201-206.

[30] Clennell M B.Tortuosity:a guide through the maze[J].Geological Society of London, Special Publications, 1997，122(1): 299-344.

[31] Gurevich A E, Chiligarian G V. Theory, measurements, and interpretation of well logs[J]. Journal of Petroleum Science & Engineering, 1995, 13(3-4):260-261.

[32] Thomas E C, Petrophysics B, Thomas E C,etal.pore-scale analysis of the waxman-smits shaly sand conductivity model[J]. Petrophysics, 2006, 48(2):104-120.

[33] 周灿灿，李潮流，王昌学，等.复杂碎屑岩测井岩石物理与处理评价[M].北京：石油工业出版社, 2013：29-31.

[34] 向阳，张朝举. 致密砂岩气藏水驱动态采收率及水膜厚度研究[J]. 成都理工学院学报, 1999, 26(4): 389-391.

New Numerical Simulation Method of Electrical Properties of Rock Based on Digital Cores

Kong Qiangfu1,Hu Song1，Wang Xiaochang1,Li Jun1, Xiong Peiqi2,Wang Wan3，Liu Mi2

(1.PetroleumExplorationandProductionResearchInstitute,SINOPEC,Beijing100083,China；2.ChinaUniversityofPetroleum,Beijing,102249,China；3.UniversityofCalgary,Calgary,AlbertaT2N4V7,Canada)

Based on the digital core, rock electrical numerical simulation technology can make up for the deficiencies of traditional rock physics experiment. Determining the fluid distribution in the pore accurately is the basis of the numerical simulation of the electrical. After carry out a wide range of literature research for the various methods, the advantages and disadvantages of combining pore-network model or mathematical morphology with other simulation methods to do the electrical numerical simulation research are analyzed. A new method that combined mathematical morphology with random walk was been put forward, and on the basis of image analysis, method of “equivalent treatment” was to characterize water-film. Connectivity of digital core, anisotropy, water film and other factors were considered, the simulation results of core with high porosity-permeability was match well with the experiment, which confirmed the feasibility of new method, laid a foundation for the electrical digital simulation study.

Digital core; Numerical simulation of electrical properties; Mathematical morphology; Random walk

国家科技重大专项“不同类型缝洞储集体测井精细评价”(2016ZX05014-002-001)。

孔强夫(1989年生)，男，硕士，主要从事测井岩石物理数值模拟研究。邮箱：kongqf.syky@sinopec.com。

TE319

A