基于数字图像处理的非均质岩石材料细观尺度应力分析

2017-01-17周汉国刘思萌郭建春彭成乐侯江朋

中国石油大学学报（自然科学版） 2016年6期

李静, 周汉国, 刘思萌, 郭建春, 彭成乐, 侯江朋

(1.中国石油大学储运与建筑工程学院,山东青岛 266580; 2.西南石油大学油气藏地质及开发工程国家重点实验室,四川成都 610500; 3.中国石化胜利油田分公司投资发展处,山东东营 257000; 4.中国石油天然气管道局,河北廊坊 065000)

李静1, 周汉国2,3, 刘思萌4, 郭建春2, 彭成乐1, 侯江朋1

以冀中坳陷碳酸盐岩储层为例,利用碳酸盐岩微观结构的图像信息,确定地应力数值模拟计算网格的合适尺寸,将数字图像处理技术引入非均质碳酸盐岩岩石材料细观尺度的应力分布计算中,将表征岩石细观结构特征的量与岩石物理力学参数等宏观的量有机结合,精确预测碳酸盐岩非均质储层裂缝发育及分布的区域。结果表明:研究区域最大水平主应力为-59.51 ～-95.09 MPa,最小水平主应力为-34.05 ～-52.58 MPa;Mises等效应力为56.84～83.15 MPa;研究区域最大总应变为2.097×10-4,最小总应变为1.392×10-4;岩心样品在左上区域和右下区域均出现应力集中和变形不协调的情况,预测为裂缝主要扩展和分布的区域;基于数字图像处理技术能够更真实、有效地表征储层岩石的非均质性。

碳酸盐岩; 非均质; 岩石断面; 数字图像处理; 数值模拟

岩石内部细观结构有着不同的物理、化学及力学性质,表现出高度的非均质性,决定了岩石在外力作用下的应力—应变响应、裂缝开展以及破坏模式[1-4]，给工程建设、矿产和油气勘探开发等带来大量难题,为此,国内外学者进行了一系列岩石细观结构的研究。Blair等[5-6]基于统计学和随机理论对岩石非均质性进行描述,研究了其破裂行为;徐文杰等[7]应用数字图像技术对土石混合体结构进行了研究,并分析了大型直剪试验的土石混合体块石含量与抗剪强度的关系; Chen等[8-9]利用数字图像处理技术,分析了非均质岩石的破坏机制;陈从新等[10]建立了岩石的细观模型,研究了岩石流—固耦合特性。将数字图像处理技术应用于油气储层岩石地应力分析和裂缝预测尚未见报道,因此,笔者以冀中坳陷碳酸盐岩储层岩石为例,基于数字图像技术结合数值模拟,进行非均质岩石材料细观尺度的地应力计算和裂缝预测研究。

1 数字图像处理基本原理

1.1 数字图像表示

任何图像都可以表示成二维函数[11]f(x,y),其中x,y是任一点的平面坐标,f为任意一点x,y处的幅度值,称为该点的亮度,也即灰度(灰度值包含了图像的基本信息,是对图像进行进一步处理的基础),此时的坐标值x、y和幅值f是连续的。将图像中连续的坐标值x、y和幅值f变为离散有限的x、y和f的过程,也即将连续的二维函数f(x,y)转化为离散的二维矩阵的过程称为图像数字化。

1.2 数字图像处理

数字图像处理的主要过程是:首先用CT机、X-Ray、扫描装置、照相机、数码相机等成像设备采集所研究物体的图像,并将其转化为储存在计算机中的数字图像(二维矩阵);然后使用计算机对其进行分析处理;最后得到所需要的信息。

数字图像处理技术有很多,其中最为基础的是空间域处理技术,空间域是指图像平面本身[12-13]。空间域处理技术的基本原理是对图像本身像素直接进行处理,其基本原理数学表述如下:

g(x,y)=T[f(x,y)].

(1)

式中,f(x,y)为输入的已知原图像;g(x,y)为处理后的结果;T为在点(x,y)的一个领域上定义的对原图像f进行处理的算子。

由于岩石中含有不同的细观介质,而这些介质一般具有不同的光学特性和光谱性质[14],因此,图像中目标对象的不同组成成分、内部结构以及缺陷就表现出不同的亮度、灰度或颜色。使用数码相机或显微镜获取岩石新鲜断面的二维数码照片,以及光谱成像仪得到光谱成像立方体,采用数字图像处理技术将这些图像进行数字化处理,数字化后的岩石图像中每个像素点的灰度f都是一个唯一的整数,式(1)中定义的图像处理算子T直接作用在灰度f上,可以将岩石中不同的细观介质清晰地辨别、提取出来[11]。

2 岩心样品的图像处理

实验用岩心样品来自任丘潜山任28井,层位属于蓟县系雾迷山组(Jxw),井深为3 243.25 m,岩性为碳酸盐岩。由于碳酸盐岩具有强烈的非均质性[15],地应力数值模拟建立模型时无法准确判断岩石介质各自的分布区域;而在岩石样品横截面(二维)和块体(三维)的数字图像中,包含了建模时所需的大量岩石介质细观(物化特性和几何结构特性)信息,如介质的种类、数量及其空间分布,因此,采用图像分割技术对岩心光片的数字图像进行处理,进一步分离出岩石中的不同介质,判断岩石样品中石灰岩与白云岩的空间分布,最终建立能够精确表征岩石不均匀性的计算模型。

2.1 图像灰度处理

通过数码相机对岩心样品进行扫描,得到如图1所示的物理图像,横向扫描线和纵向扫描线相互交叉形成具有一定宽度和高度的小区域,这些小区域是组成图像的像素点[12]。图1为宽为7.829 mm、高为17.372 mm的岩心物理图像,其横向、纵向扫描线分别为357条和351条,每条线的横向和纵向间距分别为0.022 3 mm和0.048 7 mm。

将图1导入MATLAB软件,按照式(1)的原则,利用函数rgb2gray将图像变换成灰度图像,灰度图像是一个矩阵,矩阵的值表示灰度的浓淡,由于本幅图像像素是unit8类型,所以灰度的范围是[0,255];然后利用MATLAB中自带的图像处理函数对转化成的灰度图像进行处理;最后利用imhist函数绘制出灰度直方图,如图2所示,图2表示了[0,255]256个灰度值分别在该幅图像中出现的频率。

图1 岩心的表面图像Fig.1 Surface image of cores

图2 岩心表面灰度值变量分布直方图Fig.2 Core image grey value variable distribution histogram

2.2 图像分割

图像分割是指将图像中的不同组成成分单独分割提取出来,分割的精确程度取决于实际问题的需要。就相对比较简单的单色图像而言,灰度值的不连续性和相似性[11]是对单色图像进行分割的两个主要依据。其中,若分割主要依据灰度值的不连续性,则算法主要是基于其突变进行的,选取突变处的灰度值作为阀值,对图像进行分割操作;若分割基于灰度值的相似性,则算法是根据一组事先定义好的规则,将图像分割成相似的区域。由图2可知,岩心样品图像的灰度直方图有明显的突变波峰和波谷,因此,可以采取第一类算法进行图像分割。

图像的阈值分割方法是一种对灰度值有突变或者不连续程度很大的图像进行分割的简单有效的方法[16]。灰色图像的灰度值定义为图像分割中的阈值,灰度值愈小,染色愈深,异染性颗粒含量愈多,即灰度值与异染性颗粒的含量呈反比,研究区域白云岩的灰度值小于石灰岩的灰度值。由图2可知,灰度值在120处有一个峰值,说明在120处图像中有许多像素点在此处集中,这些像素点为含有机化合物的白云岩,选择115可以把白云岩方解石脉和石灰岩区域分开。同样选择205作为第二阈值可以把白云岩从图像中分割出来,应用MATLAB软件编程,对图1进行图像分割后的结果如图3所示,两种介质在图像中被清晰地表征出来。

图3 阈值分割结果Fig.3 Threshold region segmentation method results

在图3中,共有37 800个像素点,其中白云岩和石灰岩的像素数量分别为20 105和17 695,一个像素点的实际面积为0.001 1 mm2,因此,白云岩和石灰岩的实际面积分别为22.115 5和19.464 5 mm2。图3表明,数字图像经过分割处理后能够比较精确地反映出该岩石截面中不同介质的数量和空间分布。由于在获取图像的过程中存在杂质和噪声,导致处理后的图像中白云岩区域存在部分石灰岩颗粒,石灰岩区域中存在部分白云岩的颗粒,为后面的数值模拟模型的建立带来困难,故还需要对已经处理过的图3所示的结果进行再次处理修正。去除白云岩区域中的石灰岩和石灰岩区域中的白云岩后,再对白云岩区域进行轮廓跟踪,以一个像素宽的多边形近似作为这两种介质的边界。最终的处理结果如图4所示。

将图4进行矢量化处理,应用有限元软件ANSYS建立模型,模型按照图像阈值分割方法的处理结果进行网格划分,据此进行地应力数值模拟计算。

图4 图像处理的最终结果Fig.4 Final results of digital image processing

3 地应力数值模拟计算

地应力是地壳内的岩石受到地质构造运动等作用而产生的力,岩石本身是地应力的载体,地应力通过岩石传递,因此,岩石的力学特性对地应力在地层中的传递、衰减等产生巨大影响。利用任丘28井的测井资料计算出岩石力学参数[17-19],如表1所示。岩石总体表现为脆性,故将其看成线弹性材料。

表1 任丘28井岩石力学参数Table 1 Rock mechanics parameters of well Renqiu28

3.1 计算模型

为了确定与实际情况相符合的模型边界条件,同时考虑简化计算,假定X轴所指方向代表实际的西-东方向,Y轴所指方向代表南-北方向,在西部界面施加X向位移约束,在南部界面施加Y向位移约束[20]。模型边界承受井点围压作用,根据测井资料选取井深为3 243.25 m,围压分别为120.60和50.50 MPa,模型的荷载施加情况、边界约束情况如图5(a)所示。

模型采用1∶1比例进行建模,模型X轴方向取6.04 mm,Y轴方向取7.1 mm。根据图像分割结果,应用三角形单元进行网格划分,采用四节点PLANE182单元,共生成计算单元5 361个,节点个数总共为6 283个,模型划分网格后的结果如图5(b)所示。

图5 模型荷载施加、边界约束及有限元模型Fig.5 Model loaded situation image and schematic diagram of model meshing

3.2 计算结果及其分析

对该研究区域施加均布荷载,进行地应力数值模拟计算,得到如图6、7所示的最大、最小水平主应力云图及应力方向图。

由图6(a)可知,研究区域最大水平主应力为-59.51～-95.09MPa,呈压应力状态。在左上区域和右下区域均出现应力集中情况。右下区域介质应力达到最大值,大于91.14MPa,为可能出现裂缝扩展和分布的区域。右上和左下区域的应力分布较为均匀,应力范围在-75.32 ～-83.23MPa,表明该区域不易产生裂缝。

由图6(b)可知,研究区域最小水平主应力范围在-34.05 ～-52.58MPa,呈压应力状态。应力在右上区域石灰岩处分布均匀,应力集中在-40.23 ～-42.28MPa。左上区域和右下区域处应力值较大,在左上区域产生了应力集中现象,应力值大于50.52MPa。白云岩区域应力呈带状分布,自中央区域向两侧应力值逐渐增大。岩石介质交界处最小水平主应力变化较为明显,应力值跨度范围大。

如图7所示,绿色线所指方向代表岩心截面最大主应力方向,蓝色线所指方向代表岩心截面最小主应力方向,两者呈垂直相交状态。其中右上区域和中部区域最大主应力方向为北东方向,左上区域和右下区域处最大主应力方向为北西方向,向左下区域逐渐变为南东方向。

Mises应力是一种基于剪切应变能的等效应力,其大小为

式中,σ1,σ2,σ3分别为第一、第二、第三主应力。

第四强度理论认为变形能密度是引起材料屈服破坏的原因,当Mises应力值达到材料本身的屈服应力值时,材料屈服破坏。研究区域的Mises等效应力如图8 所示。

图6 最大、最小水平主应力云图Fig.6 The maximum and minimum horizontal principal stress nephogram

图7 最大、最小水平主应力方向示意图Fig.7 Direction image of the maximum principal stress and the minimum horizontal principal stress

图8 Mises等效应力云图Fig.8 Mises stress nephogram

由图8可得,研究区域Mises等效应力为56.84～83.15 MPa,其中白云岩区域的Mises应力为69.99～83.15 MPa,明显大于石灰岩区域56.84～66.71 MPa。左上区域和右下区域的Mises应力变化范围较大,尤其是白云岩和石灰岩的交界面处。

应变云图反映了研究区域不同位置不同介质的变形情况,研究区域的总应变云图如图9所示。

图9 研究区总应变云图Fig.9 Strain nephogram of study area

从图9可以看出,整个研究区域应变为(1.392～2.097)×10-4,其中最大应变(2.097×10-4)主要发生在右下区域,最小应变(1.392×10-4)主要发生在左上区域,白云岩区域的应变明显小于石灰岩区域。

4 结论

(1)岩石内部介质的实际空间分布对其应力分布影响很大。冀中坳陷碳酸盐岩储层具有强烈的非均质性,研究区域最大主应力范围为-59.51～-95.09 MPa,最小主应力范围为-34.05 ～-52.58 MPa,应力分布不均匀,在左上区域和右下区域均出现应力集中情况,为可能出现裂缝扩展和分布的区域。

(2)研究区域Mises等效应力为56.84 ～83.15 MPa,其中白云岩区域的Mises应力为69.99～83.15 MPa,石灰岩区域的Mises应力为56.84～66.71 MPa。岩心样品左上区域和右下区域的Mises应力变化范围较大,说明该两区域会比其他区域先达到屈服应力而发生破坏形成裂缝。

(3)研究区域最大总应变(2.097×10-4)主要发生在右下区域,最小总应变(1.392×10-4)主要发生在左上区域,石灰岩和白云岩交界区域应变明显不同,变化较大,尤其是在左上区域和右下区域,故此两区域很容易产生因变形不协调而形成的裂缝。

(4)地层中岩石内部裂缝的发育程度与岩石受到的应力大体呈正比关系,岩石应力集中的区域可能是裂缝发育程度较高的区域。岩心样品左上区域和右下区域出现应力集中、变形不协调现象,预测该两区域为可能出现裂缝扩展和分布的区域。

[1] 王立恩,姜复东.碳酸盐岩储层非均质性定量表征方法[J].天然气技术,2009,3(1):27-29. WANG Lien, JIANG Fudong. Quantitatively charactering heterogeneity in carbonate rock reservoirs [J]. Natural Gas Technology, 2009,3(1):27-29.

[2] 李静,查明,刘震.基于声波测井资料的地应力分布研究:以饶阳凹陷任北奥陶系潜山为例[J].岩土力学,2011,32(9):2765-2770. LI Jing, ZHA Ming, LIU Zhen. Research on crustal stress distribution based on acoustic logging data: taking north region of Renqiu Ordovician buried hill of Raoyang Depression for example[J]. Rock and Soil Mechanics, 2011,32(9):2765-2770.

[3] 朱泽奇,肖培伟,盛谦,等.基于数字图像处理的非均质岩石材料破坏过程模拟[J].岩土力学,2011,32(12):3780-3786. ZHU Zeqi, XIAO Peiwei, SHENG Qian, et al. Numerical-simulation of fracture propagation of heterogeneous rock material based on digital image processing [J]. Rock and Soil Mechanics, 2011,32(12):3780-3786.

[4] 邹飞,李海波,周青春,等.基于数字图像灰度相关性的类岩石材料损伤分形特征研究[J].岩土力学,2012,33(3):731-738. ZOU Fei, LI Haibo, ZHOU Qingchun, et al. Fractal features study of rock-like material damage based on gray correlation of digital image [J]. Rock and Soil Mechanics, 2012,33(3):731-738.

[5] BLAIR S C, COOK N G W. Analysis of compressive fracture in rock using statistical techniques(p I): a no-linear rule based model[J]. Int J Rock Mech Min Sci, 1998,35:837-848.

[6] 唐春安,赵文.岩石破裂全过程分析软件系统RFPA2D[J].岩石力学与工程学报,1997,16(5):507-508. TANG Chunan, ZHAO Wen. RFPA2Dsystem for rock failure process analysis [J]. Chinese Journal of Rock Mechanics and Engineering, 1997,16(5):507-508.

[7] 徐文杰,胡瑞林,岳中琦.基于数字图像处理的土石混合体细观结构[J].辽宁工程技术大学学报(自然科学版),2008,27(1):51-53. XU Wenjie, HU Ruilin,YUE Zhongqi. Meso-structure character of soil-rock mixtures based on digital image [J]. Journal of Liaoning Technical University(Natural Science), 2008,27(1):51-53.

[8] CHEN S, YUE Z Q, THAM L G. Digital image-based numerical modeling method for prediction of inhomogeneous rock failure[J]. International Journal of Rock Mechanics & Mining Sciences, 2004,41(6):939-957.

[9] 于庆磊,杨天鸿,刘洪磊,等.中低围压花岗岩细观破坏机制的数值模拟[J].东北大学学报(自然科学版),2009,30(7):1026-1029. YU Qinglei, YANG Tianhong, LIU Honglei, et al. Numerical simulation on granite failure mechanism at meso-level under moderate or low confining pressure[J]. Journal of Northeastern University (Natural Science), 2009,30(7):1026-1029.

[10] 陈从新,刘秀敏,刘才华.数字图像技术在岩石细观力学研究中的应用[J].岩土力学,2010,31(增1):53-61. CHEN Congxin, LIU Xiumin, LIU Caihua. Application of digital image processing to rock mesomechanics [J]. Rock and Soil Mechanics, 2010,31(sup1):53-61.

[11] RAFAEL C G, RICHARD E W, STEVEN L E. Digital image processing using MATLAB [M]. 2nd ed.Beijing: Publishing House of Electronics Industry, 2014.

[12] 胡学龙.数字图像处理[M].北京:电子工业出版社,2014.

[13] 杨松,邵龙潭,郭晓霞,等.基于混凝土裂纹数字图像的有限元网格生成[J].计算力学学报,2012,29(4):635-640. YANG Song, SHAO Longtan, GUO Xiaoxia, et al. Finite element mesh generation based on digital image for concrete crack[J]. Chinese Journal of Computational Mechanics, 2012,29(4):635-640.

[14] 李静,查明,郭元岭,等.光谱成像信息的数据融合技术在储层表征中的应用[J].光谱学与光谱分析,2011,31(10):2639-2642. LI Jing, ZHA Ming, GUO Yuanling,et al. Application of data fusion of microscopic spectral imaging in reservoir characterization [J]. Spectroscopy and Spectral Analysis,2011, 31(10):2639-2642.

[15] 王子振,王瑞和,单殉,等.常规方法预测碳酸盐岩地层压力的偏差分析[J].中国石油大学学报(自然科学版),2014,38(5):96-101. WANG Zizhen, WANG Ruihe, SHAN Xun, et al. Uncertainty analysis of pore pressure prediction in carbonate formation using conventional methods[J]. Journal of China University of Petroleum (Edition of Natura1 Science), 2014,38(5):96-101.

[16] YUE Z Q, CHEN S, THAM L G. Finite element modeling of geomaterials using digital image processing[J]. Computers and Geotechnics, 2003,30(5):375-397.

[17] 胡国忠,王宏图,贾剑青.岩石的动静弹性模量的关系[J].重庆大学学报(自然科学版),2005,28(3):102-105. HU Guozhong, WANG Hongtu, JIA Jianqing. Relationship between dynamic and static value of elastic modulus in rock[J]. Journal of Chongqing University (Natural Science Edition), 2005,28(3):102-105.

[18] 盛金昌,刘继山,速宝玉.基于图像数字化技术的裂隙岩石多场耦合分析[J].工程力学,2007,24(10):30-35. SHENG Jinchang, LIU Jishan, SU Baoyu. Coupled multiphysics analysis in fractured rock masses based on digital image processing technique[J]. Engineering Mechanics, 2007,24(10):30-35.

[19] 孙建孟,闫国亮,姜黎明,等.基于数字岩心研究流体性质对裂缝性低渗透储层弹性参数的影响规律[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.

[20] 戴俊生,刘敬寿,杨海盟,等.铜城断裂带阜二段储层应力场数值模拟及开发建议[J].中国石油大学学报(自然科学版),2016,40(1):1-9. DAI Junsheng, LIU Jingshou, YANG Haimeng, et al. Numerical simulation of stress field of Fu-2 member in Tongcheng fault zone and development suggestions[J]. Journal of China University of Petroleum (Edition of Natural Science), 2016,40(1):1-9.

(编辑沈玉英)

Stress analysis of heterogeneous rock material at mesoscopic scale based on digital image processing technology

LI Jing1, ZHOU Hanguo2,3, LIU Simeng4, GUO Jianchun2, PENG Chengle1, HOU Jiangpeng1

(1.College of Pipeline and Civil Engineering in China University of Petroleum, Qingdao 266580, China;2.StateKeyLaboratoryofOilandGasReservoirGeologyandExploitation,SouthwestPetroleumUniversity,Chengdu610500,China;3.InvestmentDevelopmentDepartmentofSINOPECShengliOilfieldCompany,Dongying257000,China;4.ChinaPetroleumPipelineBureau,Langfang065000,China)

Taking the carbonate reservoir of Jizhong depression as an example in this paper, the image information of carbonate rock microscopic structure was used to determine the appropriate grid size for the numerical simulation of ground stress field. The digital image processing technology was introduced into the distribution calculation of meso-scale stress of heterogeneous carbonate rock material. Combining those quantities characterizing the microscopic structure feature of rock with the physical and mechanical parameters of rock and other macroscopic quantities, the fracture developing and distributing area in carbonate heterogeneous reservoirs can be predicted accurately. The results show that the maximum horizontal principal stress range is -59.51--95.09 MPa, and the minimum horizontal principal stress is -34.05--52.58 MPa. Mises equivalent stress is 56.84-83.15 MPa. The maximum total strain of the studied area is 2.097×10-4, and the minimum total strain is 1.392×10-4. The stress concentration and incompatible deformation occur in the up-left and the down-right regions of the core samples, which is predicted to be the main developing and distributing area of the fractures. The numerical calculation model based on digital image processing technology can reflect the real microscopic structure of reservoir rock material, and characterize the heterogeneity of reservoir rock more reliably and effectively.

carbonate; heterogeneous; rock section; digital image processing; numerical simulation

2016-03-13

国家自然科学基金项目(41272141)；国家科技重大专项(2016ZX05002-002)

李静(1967-)，女，教授，博士，硕士生导师，主要研究方向为岩石力学与储层岩石裂缝预测。E-mail：lijing0681@163.com。

1673-5005(2016)06-0143-07

10.3969/j.issn.1673-5005.2016.06.018

TE 122.2

李静,周汉国,刘思萌,等. 基于数字图像处理的非均质岩石材料细观尺度应力分析[J]. 中国石油大学学报(自然科学版), 2016,40(6):143-149.

LI Jing, ZHOU Hanguo, LIU Simeng, et al. Stress analysis of heterogeneous rock material at mesoscopic scale based on digital image processing technology [J]. Journal of China University of Petroleum (Edition of Natural Science), 2016,40(6):143-149.