CN108956416A - 一种基于Matlab分析致密砂岩储层孔隙表征的方法 - Google Patents
一种基于Matlab分析致密砂岩储层孔隙表征的方法 Download PDFInfo
- Publication number
- CN108956416A CN108956416A CN201810571899.0A CN201810571899A CN108956416A CN 108956416 A CN108956416 A CN 108956416A CN 201810571899 A CN201810571899 A CN 201810571899A CN 108956416 A CN108956416 A CN 108956416A
- Authority
- CN
- China
- Prior art keywords
- image
- hole
- matlab
- function
- sandstone reservoir
- 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 53
- 238000012512 characterization method Methods 0.000 title claims abstract description 30
- 238000004458 analytical method Methods 0.000 title claims abstract description 27
- 239000011148 porous material Substances 0.000 claims abstract description 77
- 238000009826 distribution Methods 0.000 claims abstract description 35
- 238000001878 scanning electron micrograph Methods 0.000 claims abstract description 20
- 238000005266 casting Methods 0.000 claims abstract description 19
- 238000006243 chemical reaction Methods 0.000 claims abstract description 4
- QSHDDOUJBYECFT-UHFFFAOYSA-N mercury Chemical compound [Hg] QSHDDOUJBYECFT-UHFFFAOYSA-N 0.000 claims description 5
- 229910052753 mercury Inorganic materials 0.000 claims description 5
- 238000000605 extraction Methods 0.000 claims description 3
- 238000004364 calculation method Methods 0.000 claims description 2
- 238000004445 quantitative analysis Methods 0.000 abstract description 9
- 238000011160 research Methods 0.000 abstract description 5
- 239000004576 sand Substances 0.000 description 12
- 239000000243 solution Substances 0.000 description 12
- 239000011800 void material Substances 0.000 description 9
- NLYAJNPCOHFWQQ-UHFFFAOYSA-N kaolin Chemical compound O.O.O=[Al]O[Si](=O)O[Si](=O)O[Al]=O NLYAJNPCOHFWQQ-UHFFFAOYSA-N 0.000 description 8
- 229910052622 kaolinite Inorganic materials 0.000 description 8
- 239000002245 particle Substances 0.000 description 8
- 238000010586 diagram Methods 0.000 description 7
- 230000007797 corrosion Effects 0.000 description 5
- 238000005260 corrosion Methods 0.000 description 5
- 229910052900 illite Inorganic materials 0.000 description 5
- VGIBGUSAECPPNB-UHFFFAOYSA-L nonaaluminum;magnesium;tripotassium;1,3-dioxido-2,4,5-trioxa-1,3-disilabicyclo[1.1.1]pentane;iron(2+);oxygen(2-);fluoride;hydroxide Chemical compound [OH-].[O-2].[O-2].[O-2].[O-2].[O-2].[F-].[Mg+2].[Al+3].[Al+3].[Al+3].[Al+3].[Al+3].[Al+3].[Al+3].[Al+3].[Al+3].[K+].[K+].[K+].[Fe+2].O1[Si]2([O-])O[Si]1([O-])O2.O1[Si]2([O-])O[Si]1([O-])O2.O1[Si]2([O-])O[Si]1([O-])O2.O1[Si]2([O-])O[Si]1([O-])O2.O1[Si]2([O-])O[Si]1([O-])O2.O1[Si]2([O-])O[Si]1([O-])O2.O1[Si]2([O-])O[Si]1([O-])O2 VGIBGUSAECPPNB-UHFFFAOYSA-L 0.000 description 5
- 238000011161 development Methods 0.000 description 4
- 238000002474 experimental method Methods 0.000 description 4
- 230000035699 permeability Effects 0.000 description 4
- 238000002347 injection Methods 0.000 description 3
- 239000007924 injection Substances 0.000 description 3
- 239000011435 rock Substances 0.000 description 3
- 241000208340 Araliaceae Species 0.000 description 2
- 238000005481 NMR spectroscopy Methods 0.000 description 2
- 235000005035 Panax pseudoginseng ssp. pseudoginseng Nutrition 0.000 description 2
- 235000003140 Panax quinquefolius Nutrition 0.000 description 2
- 239000002734 clay mineral Substances 0.000 description 2
- 239000003245 coal Substances 0.000 description 2
- 239000012141 concentrate Substances 0.000 description 2
- 238000012790 confirmation Methods 0.000 description 2
- 230000000694 effects Effects 0.000 description 2
- 238000005516 engineering process Methods 0.000 description 2
- 239000012530 fluid Substances 0.000 description 2
- 235000008434 ginseng Nutrition 0.000 description 2
- 239000007788 liquid Substances 0.000 description 2
- 244000132059 Carica parviflora Species 0.000 description 1
- 235000014653 Carica parviflora Nutrition 0.000 description 1
- 208000035126 Facies Diseases 0.000 description 1
- PPTYJKAXVCCBDU-UHFFFAOYSA-N Rohypnol Chemical compound N=1CC(=O)N(C)C2=CC=C([N+]([O-])=O)C=C2C=1C1=CC=CC=C1F PPTYJKAXVCCBDU-UHFFFAOYSA-N 0.000 description 1
- 239000000853 adhesive Substances 0.000 description 1
- 230000001070 adhesive effect Effects 0.000 description 1
- 238000013459 approach Methods 0.000 description 1
- 230000009286 beneficial effect Effects 0.000 description 1
- 239000004568 cement Substances 0.000 description 1
- 229910001919 chlorite Inorganic materials 0.000 description 1
- 229910052619 chlorite group Inorganic materials 0.000 description 1
- QBWCMBCROVPCKQ-UHFFFAOYSA-N chlorous acid Chemical compound OCl=O QBWCMBCROVPCKQ-UHFFFAOYSA-N 0.000 description 1
- 239000004927 clay Substances 0.000 description 1
- 239000003086 colorant Substances 0.000 description 1
- 238000005520 cutting process Methods 0.000 description 1
- 230000007812 deficiency Effects 0.000 description 1
- 230000008021 deposition Effects 0.000 description 1
- 230000005611 electricity Effects 0.000 description 1
- 239000004744 fabric Substances 0.000 description 1
- 239000003292 glue Substances 0.000 description 1
- 239000004615 ingredient Substances 0.000 description 1
- 229910052500 inorganic mineral Inorganic materials 0.000 description 1
- 238000005259 measurement Methods 0.000 description 1
- 239000011707 mineral Substances 0.000 description 1
- 238000007781 pre-processing Methods 0.000 description 1
- 238000002360 preparation method Methods 0.000 description 1
- 238000004321 preservation Methods 0.000 description 1
- 238000003860 storage Methods 0.000 description 1
Classifications
-
- G—PHYSICS
- G01—MEASURING; TESTING
- G01N—INVESTIGATING OR ANALYSING MATERIALS BY DETERMINING THEIR CHEMICAL OR PHYSICAL PROPERTIES
- G01N15/00—Investigating characteristics of particles; Investigating permeability, pore-volume or surface-area of porous materials
- G01N15/08—Investigating permeability, pore-volume, or surface area of porous materials
- G01N15/088—Investigating volume, surface area, size or distribution of pores; Porosimetry
-
- G—PHYSICS
- G01—MEASURING; TESTING
- G01N—INVESTIGATING OR ANALYSING MATERIALS BY DETERMINING THEIR CHEMICAL OR PHYSICAL PROPERTIES
- G01N15/00—Investigating characteristics of particles; Investigating permeability, pore-volume or surface-area of porous materials
- G01N15/08—Investigating permeability, pore-volume, or surface area of porous materials
Landscapes
- Chemical & Material Sciences (AREA)
- Dispersion Chemistry (AREA)
- Physics & Mathematics (AREA)
- Health & Medical Sciences (AREA)
- Life Sciences & Earth Sciences (AREA)
- Analytical Chemistry (AREA)
- Biochemistry (AREA)
- General Health & Medical Sciences (AREA)
- General Physics & Mathematics (AREA)
- Immunology (AREA)
- Pathology (AREA)
- Analysing Materials By The Use Of Radiation (AREA)
Abstract
本发明公开了一种基于Matlab分析致密砂岩储层孔隙表征的方法，首先读取致密砂岩储层的铸体薄片图像和SEM图像；然后调节图像阈值，将铸体薄片图像和SEM图像进行二值化处理得到二值图像；再通过识别二值图像的黑色像素点，计算孔隙度；并求取致密砂岩储层的孔隙半径分布；最后求取致密砂岩储层的孔洞平均直径及方差，并输出二值图像和求取的数据，应用软件Matlab对致密砂岩储层孔隙表征研究主要分为两种方法，第一种利用内置函数对二值图像进行定量分析，适用于纳米级孔隙与微米级孔隙，第二种利用应用linecut即切片法直接对铸体薄片、SEM图像等进行分析，适用于纳米级孔隙。
Description
技术领域
本发明涉及致密砂岩储层孔隙表征分析技术领域，具体为一种基于Matlab分析致密砂岩储层孔隙表征的方法。
背景技术
储集层微观孔隙结构是指储集岩石中孔隙和吼道的几何形状、大小、分布及其相互连通关系，直接影响储层孔隙空间的发育状况。储层储集空间分为毫米孔隙、微米孔隙、纳米孔隙三种级别，才能等通过应用场发射扫描电子显微镜与纳米CT重构技术对中国非常规储层研究中首次发现纳米孔，国内外多名学者对致密砂岩储层孔隙结构、表征等进行研究，多种手段被应用于研究致密砂岩孔隙，以核磁共振NMR测井技术、流体注入实验法和图像观测法这三大类为主。
基于岩石电阻率参数的核磁共振测井技术，以测定储层孔径分布为主。以压汞法为代表的流体注入实验法能够间接地通过压汞曲线获得储集空间分布、孔径大小等参数，但并不能像以由于SEM及纳米CT技术等图像观测法直接地观测到储集空间分布、孔喉特征等特征。
然而以扫描电镜及纳米CT为实验指导，研究结果大多为对孔喉的几何形状及联通关系进行定性描述，对孔隙表征的定量研究仍有不足。
发明内容
为了克服现有技术方案的不足,本发明提供一种基于Matlab分析致密砂岩储层孔隙表征的方法，应用软件Matlab对致密砂岩储层孔隙表征研究主要分为两种方法，第一种利用内置函数对二值图像进行定量分析，适用于纳米级孔隙与微米级孔隙。第二种利用应用linecut即切片法直接对铸体薄片、SEM图像等进行分析，适用于纳米级孔隙，能有效的解决背景技术提出的问题。
本发明解决其技术问题所采用的技术方案是：
一种基于Matlab分析致密砂岩储层孔隙表征的方法，包括如下步骤：
步骤100、读取致密砂岩储层的铸体薄片图像和SEM图像；
步骤200、调节图像阈值，将铸体薄片图像和SEM图像进行二值化处理得到二值图像；
步骤300、识别二值图像的黑色像素点，计算孔隙度；
步骤400、求取致密砂岩储层的孔隙半径分布；
步骤500、求取致密砂岩储层的孔洞平均直径及方差，并输出二值图像和求取的数据。
作为本发明一种优选的技术方案，所述步骤200中图像阈值的调节方法包括：通过软件Photoshop调整图像阈值，以及利用Matlab中imread函数对预处理图像的像素进行读取，在调用Matlab内置应用threshold将图像转为二值图像，进入YCbCr模式。
作为本发明一种优选的技术方案，所述步骤300中孔隙度的计算方法包括：
首先通过Matlab中的size函数对二值图像中的像素进行识别提取，通过函数可以分别获得图像的总像点数Na；
然后调用能够识别0像素值的函数bwarea获取孔隙像素点数Np；
孔隙度P即黑色像素点与总像素点的比值，即为：
作为本发明一种优选的技术方案，所述步骤400中孔隙半径分布的求取方法包括两种：
方法一，利用Matlab中内置函数bwlabel进行求取；
方法二，利用Matlab内置应用GUI_linecut，即切片法进行求取。
作为本发明一种优选的技术方案，所述方法一的具体步骤如下：
步骤411、利用Matlab内置函数bwlabel可对相同像素值的区域进行识别统计，同时可以确定区域的边界，基于此原理确定孔隙区域；
步骤412、对每一个二值图像的孔洞调用bwlabel函数，求出每个区域X轴和Y轴的最大值以及最小值即xmax、xmin、ymax、ymin，所以每一个孔洞的直径大小d为X和Y方向的差值的平均值，即孔径大小d：
步骤413、调用Matlab内部函数imhist对所有孔隙进行叠加并绘制孔径分布直方图；
作为本发明一种优选的技术方案，所述方法二包括：
步骤421、调用Matlab内部程序linecut，将预处理图像进行切片；
步骤422、在每条切线上手动增加控制点，并且设置比例尺，最终将控制点导出成直方图；
步骤423、对切线控制区域内的孔隙进行识别输入，可得到切线上被识别的线段，即为孔径大小；
作为本发明一种优选的技术方案，所述步骤500中孔洞平均直径及方差的求取是通过分别调用函数mean以及函数var进行求取，即为平均孔隙半径以及方差S。
作为本发明一种优选的技术方案，所述方法一中孔洞平均直径及方差为：
作为本发明一种优选的技术方案，所述图像阈值是通过函数imhist建立二值图像直方图,人为的控制阈值大小，把误差控制到最低。
作为本发明一种优选的技术方案，所述孔隙度与孔径分布分析结果通过应用软件Photoshop与压汞曲线对实验结果准确性进行验证，通过软件Photoshop套索及色彩范围功能对Y488致密储层孔隙度像素进行分析，利用色彩范围功能对孔隙进行选择。
与现有技术相比，本发明的有益效果是：本发明应用软件Matlab对致密砂岩储层孔隙表征研究主要分为两种方法，第一种利用内置函数对二值图像进行定量分析，适用于纳米级孔隙与微米级孔隙。第二种利用应用linecut即切片法直接对铸体薄片、SEM图像等进行分析，适用于纳米级孔隙。
附图说明
图1为本发明Matlab中图像二值化界面图；
图2为本发明利用二值图像定量分析微米级孔隙的示意图；
图3为本发明Matlab中GUI_linecut切片法操作界面图；
图4为本发明的流程图；
图5为本发明实施方式中盒8段致密砂岩微米级孔隙的示意图；
图6为本发明实施方式中盒8段致密砂岩纳米级孔隙的示意图；
图7为本发明实施方式中盒8段伊利石晶间孔定量分析图；
图8为本发明不同阈值控制的孔径大小分布示意图；
图9为本发明利用软件Photoshop进行对孔隙进行识别示意图；
图10为本发明函数法与切片法方法比较示意图。
具体实施方式
下面将结合本发明实施例中的附图，对本发明实施例中的技术方案进行清楚、完整地描述，显然，所描述的实施例仅仅是本发明一部分实施例，而不是全部的实施例。基于本发明中的实施例，本领域普通技术人员在没有做出创造性劳动前提下所获得的所有其他实施例，都属于本发明保护的范围。
如图4所示，本发明提供了一种基于Matlab分析致密砂岩储层孔隙表征的方法，包括如下步骤：
步骤一、读取致密砂岩储层的铸体薄片图像和SEM图像；
铸体薄片图像和SEM图像在这里作为分析储层孔隙相关表征的手段，其中，铸体薄片是通过将有色液态胶在真空加压下注入岩石孔隙空间，待液态胶固化后磨制成的岩石薄片，由于岩石孔隙被有色胶充填，故在显微镜下十分醒目，容易辨认；SEM图像是通过扫描电子显微镜读取的致密砂岩储层图像。
步骤二、调节图像阈值，将铸体薄片图像和SEM图像进行二值化处理得到二值图像；
在该步骤中，由于计算机图像的分析处理，需要将图像预处理为二值图像，而图像二值化就是将图像上的像素点的灰度值设置为0或255，也就是将整个图像呈现出明显的黑白效果，二值图像最大幅度地显示图像内部结构。
对应在密砂岩储层的图像时，像素点为0的点显示黑色，代表孔隙区域，像素点为255的点显示为白色，代表非孔隙区域(包括颗粒、杂基及胶结物)。
将图像二值化的手段也主要有两种，一种即通过软件Photoshop调整图像阈值；另一种利用Matlab中imread函数对预处理图像的像素进行读取，在调用内置应用threshold将图像转为二值图像，进入YCbCr模式，如图1所示。
图1中，右上角为控制二值图像的阈值范围，Y为颜色亮度成分，Cb与Cr分别为蓝色和红色的浓度偏移，通过界面上功能将转为二值图像，Y起着阈值作用，调节Y值还原图像孔隙空间，阈值的调节对致密砂岩储层表征起着重要的影响。
通过预处理获取的二值图像如图2所示，根据二值图像可以进一步计算图像的参数。
步骤三、通过Matlab中的size函数对二值图像中的像素进行识别提取，通过函数可以分别获得图像的总像点数Na(所有黑白像素点总和)；然后调用能够识别0像素值的函数bwarea获取孔隙像素点数Np(黑色像素点)；孔隙度P即黑色像素点与总像素点的比值，即为：
步骤400、求取致密砂岩储层的孔隙半径分布，由于储层致密砂岩孔隙的复杂性、非均质性以及相对低孔隙度(<12％)，对于致密砂岩储层的孔隙半径分布求取时有两种途径。
其一：利用Matlab中内置函数bwlabel进行求取；
利用Matlab内置函数bwlabel可对相同像素值的区域进行识别统计，同时可以确定区域的边界，基于此原理确定孔隙区域(黑色像素值连通区域)。对每一个孔洞调用bwlabel函数，求出每个区域的X和Y轴的最大值以及最小值即xmax、xmin、ymax、ymin，所以对于每一个孔洞的直径大小d为X和Y方向的差值的平均值，即孔径大小d：
进一步调用内部函数imhist对所有孔隙进行叠加绘制孔径分布直方图，直方图可清晰地展示孔径直方图。最终分别调用函数mean以及函数var求取平均孔隙半径以及方差S:
其二：利用Matlab内置应用GUI_linecut，即切片法进行求取，如图3所示：
将预处理图像进行切片，在每条且线上手动增加控制点，并且设置比例尺，最终将控制点导出成直方图。
该种方法适用性较高，对图像要求较少，因为只有一个X方向控制，所以精确度较差，由于切片法是对图像进行分割，相比第一种方法因计算机自动识别孔隙坐标，切片法适用于孔喉分布相对均匀，弯曲度较低的图像。
虽然切片法的应用有其自身的局限性，但由于此方法操作简单，可对铸体薄片图像和SEM图像直接进行分析。
以SEM图像分析为例，加载图像后需调节图像比例尺级别，进一步控制图像切线数量，并对切线控制区域内的孔隙进行识别输入，切线上被识别的线段即为孔径大小。
上述两种途径最终可以获得孔径分布直方图，进一步对孔隙平均值、非均质性以及连通性进行分析，通过调用Matlab内置函数mean函数可对平均孔隙半径进行求取。在获取每个孔洞半径的基础上调用内置函数var函数对所有孔洞的方差进行求取，平均孔隙半径及其方差S。
对于致密砂岩储层非均质性来说，方差越大，表现整个储层孔隙分布越复杂，非均质性越强。利用方差S来定量判断储层非均质性需要借助微米级图像，对于纳米级图像如SEM图像，由于图像尺寸相对较小，图像所选取的储层区域缺少代表性。通过对致密砂岩储层图像二值化预处理后，通过Matlab相关内置函数可以定量求取储层砂岩孔隙度、孔隙半径分布、平均孔隙度、非均质性等表征参数。
步骤五、输出二值图像和求取的数据。
以下以鄂尔多斯盆地东部上古生界致密砂岩为例进行说明。
鄂尔多斯盆地资源丰富，煤系地层发育，具有丰富的煤系气资源，上古生界为典型的海陆过渡相沉积，盒8致密砂岩作为上古生界储集空间，发育致密砂岩气储层，选区井位为Y488井。致密砂岩储层孔隙度较低，孔隙结构复杂，微米孔隙与纳米孔隙并存，微米孔隙以次生溶蚀孔、较大的粒间孔、微裂缝等为主，充当储层主要的储集空间，孔隙直径相对较大。纳米孔隙以岩屑粒内孔、自生矿物晶间孔如高岭石晶间孔等为主，作用于联通微米孔隙的孔喉，主导着储层渗透性，孔隙结构复杂。
通过铸体薄片与SEM图片，对储层微米级孔隙及纳米级孔隙表征进行定量研究，分析储层在不同尺度上孔隙大小分布及联通性等相关参数。由于微米级孔隙与纳米级孔隙具有不同的尺寸级别，所展现的孔隙空间、结构以及大小也大不一样，分别以两种不同尺度上储层孔隙表征进行说明。
第一种，微米级孔隙。
研究区微米级空隙以次生溶蚀孔和粒间孔为主，选区Y488盒8段目标层位铸体薄片一张，首先利用内置函数imread读取图像，读取图像后通过Matlab内置应用threshold，将图像转化为二值图像同时调节图像阈值，所调节的阈值应精准地显示完整孔隙，生成的二值图像如图5(b)。建立二值图像后，根据公式(1)及函数bwarea求解得到孔隙度，如图5(b)。基于Matlab所识别的孔径分布图，见图5(c)，引用mean函数以及var函数对图像所有孔径进行均质及方差的求解。
由图5(b)测得储层孔隙度为9.5825％，图5(c)显示了Matlab所识别的储集空间范围，由内置函数bwlabel计算所得，
孔隙直径主要分布在200μm以下，平均孔隙半径为158.4357μm，储层储集空间主要以次生溶蚀孔为主，集中在孔径大小200μm以上这部分孔隙，200μm以下这部分微米孔隙较小，集中在储层吼道部位，控制着储层渗透率。孔径方差为118.1767，孔隙分布较均匀，联通性较差。应用函数bwlabel计算致密砂岩孔隙表征时，由于微米级图像二值化时引起的噪点误差，使得大量岩屑粒内孔被误认为颗粒，粘土矿物晶间孔等纳米级孔隙会直接与粘土矿物一同被系统误认为孔隙，使得整个图像微米级孔隙消失，整体孔径大小偏大。因此要借助SEM图像对储层中纳米级孔隙进行定量表征。
第二种，纳米级孔隙。
据邹才能等学者提出的鄂尔多斯盆地致密砂岩储层孔喉直径主要集中在40-700nm之间，与微米孔隙相比，纳米孔隙结构相对复杂，非均质性也更强，针对能够展示微米级孔隙的SEM图像，由于SEM图像所呈现的储集空间为纳米级孔隙，其复杂性及非均质性也远大于微米孔隙，因此纳米级图像所呈现的孔隙度也并不具有代表性，这里主要研究孔径分布，平均孔喉大小以及非均质性。研究区纳米级孔隙以高岭石晶间孔及岩屑粒内孔为代表，局部发育少量伊利石/绿泥石混，分别应用切片法以及内置函数bwlabel对两种孔隙类型进行研究。
高岭石晶间孔：高岭石呈书珊状赋存于储层颗粒之间，高岭石晶间孔分布相对均匀如图6(d)，孔隙叠加区较少，所以应用切片法定量分析高岭石晶间孔。高岭石晶间孔大小主要集中于50～400nm，平均大小为163.2961nm，方差为68.6932，见图6(f)，孔隙分布相对均匀，对于纳米级图像，由于比例尺非常小，图像非均质性非常强，所以对于纳米级孔隙只分析其孔径大小及分布状况。
岩屑粒内孔：与高岭石晶间孔相比，岩屑粒内孔孔径分布非均质性较强，由Matlab内置函数分析，岩屑平均孔隙大小为153.0664nm，标准差为232.8941，孔径分布范围较广。
伊利石晶间孔：储层中纤维状伊利石晶间孔结构非常复杂，大量纳米孔隙叠合在一起，结构复杂，切片法并不适用。因此应用Matlab内置函数定量分析伊利石孔隙，二值图像如图7(b)所示，利用内置函数size、bwlabel等对储层孔隙进行计算，最终孔隙图像输出如图7(c)，可以发现大量直径较小孔隙的未被识别出，这是由于图像分辨率以及比例尺所决定的，因此所测的孔径分布、平均孔隙直径稍微偏大。
整个盒8段致密砂岩储集空间孔隙由微米级孔隙与纳米级孔隙共同构成，综上，目标层位主要储集空间集中在100～300μm，以次生溶蚀孔为主，吼道主要集中在150nm左右，以粘土矿物晶间孔为主，储层孔隙度为11.984％，整个储层孔隙分布较均匀，但联通性差，吼道直径较低。目标层位致密储层砂岩储集空间由微米级孔洞及纳米级孔喉构成，前者控制这储层空间发育，后者影响储层渗透率，两者皆对储层起着至关重要的作用，研究储层储集空间离不开微米孔隙，研究储层渗透性也建立在纳米级孔隙的基础之上，两者相互依托，共同作用于致密砂岩储层。
在本实施方式中，二值图像的阈值是用来控制图像细节的“阀门”，阈值越高，图像细节越多，这种现象使得Matlab计算图像孔隙度时，代表致密砂岩中颗粒与杂基的白色像素也越多，因阈值升高而产生的噪点代替原来颗粒与杂基的位置，产生很多代表噪点的黑色像素值，使得Matlab计算孔隙度时会产生大量孔径半径较小的孔隙，孔喉结构变得复杂，所测的的孔隙度、孔径大小、联通性等参宿均小于理论值。相反地，阈值越低，图像细节越少，孔隙空间被代表颗粒及杂基的白色像素所代替，孔隙空间大幅度减少，因此所测的孔隙度、孔径大小、颗粒联通性均等参数均大于理论值。综上所述，无论阈值偏大或偏小都会使得所测孔隙度、孔径半径偏小，所以通过函数imhist建立二值图像直方图，人为的控制阈值大小非常有必要，把误差控制到最低。
以目标层位致密砂岩储层铸体薄片为例，阈值为112时，见图8(a)，图像孔隙分布与铸体图像最为接近。二值图像阈值为75时，见图8(b)，储层孔隙度为6.2967％，与标准值偏差较多，平均空隙大小为156.4357μm，与标准孔径大小相比偏小，孔孔方差偏大，孔径分布更差。二值图像阈值为166时，见图8(c)，图像细节增多，二值图像孔隙图明显偏大，如图，同时孔径边缘产生大量噪点，对应出现大量假性“孔洞”，孔径方差为131.3525更加杂乱，其分析结果如下：
在本实施方式中，应用Matlab对致密砂岩储层孔隙的定量化研究中，对实验对象进行计算时，误差是一个不可避免的影响，影响误差的因素包括二值图像阈值的确认、应用函数bwlabel计算孔洞大小、应用切片法测量孔隙半径等步骤。对于二值图像阈值的确认已在上文中解释，需通过二值图像直方图人为调节阈值。而在应用Matlab内置函数bwlabel对孔洞大小进行计算时，函数对孔洞xmin、xmax、ymin、ymax这四个参数进行控制计算，对于一些规则储集空间如颗粒溶蚀孔非常准确，但对于一些条带状、片状等歪度较为复杂的孔隙会产生误差，使得测得孔隙半径偏大一些。对于条带状、片状、三角状等复杂的孔隙而言，切片法所测结果会精准一些，但由于切片法是对图像某一方向(x轴或y轴)进行测量，同时图像切割条数也会影响到孔隙测量次数，所以切片法结果具有一些偶然性，需要进行多次统计。
综上，实验对于孔隙度与孔径分布分析结果具有一定偏差，应用Photoshop与压汞曲线对实验结果准确性进行验证。多位学者应用Photoshop对定量分析储层孔隙度[24]，分析结果误差非常低。本文依靠软件Photoshop套索及色彩范围功能对Y488致密储层孔隙度像素进行分析，利用色彩范围功能对“红色”孔隙进行选择，如图9，根据直方图辨别出，孔隙所占像素值164211，总图像像素值1783545，因此孔隙度为9.2721％。而利用Matlab所测孔隙度为9.5825％，偏差为3.2392％，误差相对较低，并且可进行批量加载计算，适用范围广。
函数法与切片法分析图像比较如图10所示，致密砂岩孔隙空间由纳米级孔隙及微米孔隙共同构成，与纳米级孔隙相比，微米级图像所显示的孔隙结构更为复杂，能够展现出更多细小的孔隙，分析孔径分布主要用到Matlab内置函数bwlabel与GUI切片法。利用函数bwlabel需建立在二值图像上，由于形成二值图像会产生误差，所以生成的结果也会相应偏离原孔隙表征参数，但由于其适用性高，可应用于一切二值图像。切片法是利用切线将图像分割成多个区域，每个孔径大小由平行于切线方向孔洞最大宽度所决定，越复杂的孔隙结构，在每个区域孔隙叠合率高，并且部分小尺寸孔隙并未被切线相切，但其不需建立在二值图像上便可进行操作，避免了生成二值图像所带来的误差，适用于微米级图像。
应用软件Matlab对致密砂岩储层孔隙表征研究主要分为两种方法，第一种利用内置函数对二值图像进行定量分析，适用于纳米级孔隙与微米级孔隙。第二种利用应用linecut即切片法直接对铸体薄片、SEM图像等进行分析，适用于纳米级孔隙。
对于本领域技术人员而言，显然本发明不限于上述示范性实施例的细节，而且在不背离本发明的精神或基本特征的情况下，能够以其他的具体形式实现本发明。因此，无论从哪一点来看，均应将实施例看作是示范性的，而且是非限制性的，本发明的范围由所附权利要求而不是上述说明限定，因此旨在将落在权利要求的等同要件的含义和范围内的所有变化囊括在本发明内。不应将权利要求中的任何附图标记视为限制所涉及的权利要求。
Claims (10)
1.一种基于Matlab分析致密砂岩储层孔隙表征的方法，其特征在于，包括如下步骤：
步骤100、读取致密砂岩储层的铸体薄片图像和SEM图像；
步骤200、调节图像阈值，将铸体薄片图像和SEM图像进行二值化处理得到二值图像；
步骤300、识别二值图像的黑色像素点，计算孔隙度；
步骤400、求取致密砂岩储层的孔隙半径分布；
步骤500、求取致密砂岩储层的孔洞平均直径及方差，并输出二值图像和求取的数据。
2.根据权利要求1所述的一种基于Matlab分析致密砂岩储层孔隙表征的方法，其特征在于：所述步骤200中图像阈值的调节方法包括：通过软件Photoshop调整图像阈值，以及利用Matlab中imread函数对预处理图像的像素进行读取，在调用Matlab内置应用threshold将图像转为二值图像，进入YCbCr模式。
3.根据权利要求1所述的一种基于Matlab分析致密砂岩储层孔隙表征的方法，其特征在于：所述步骤300中孔隙度的计算方法包括：
首先通过Matlab中的size函数对二值图像中的像素进行识别提取，通过函数可以分别获得图像的总像点数Na；
然后调用能够识别0像素值的函数bwarea获取孔隙像素点数Np；
孔隙度P即黑色像素点与总像素点的比值，即为：
4.根据权利要求1所述的一种基于Matlab分析致密砂岩储层孔隙表征的方法，其特征在于：所述步骤400中孔隙半径分布的求取方法包括两种：
方法一，利用Matlab中内置函数bwlabel进行求取；
方法二，利用Matlab内置应用GUI_linecut，即切片法进行求取。
5.根据权利要求4所述的一种基于Matlab分析致密砂岩储层孔隙表征的方法，其特征在于：所述方法一的具体步骤如下：
步骤411、利用Matlab内置函数bwlabel可对相同像素值的区域进行识别统计，同时可以确定区域的边界，基于此原理确定孔隙区域；
步骤412、对每一个二值图像的孔洞调用bwlabel函数，求出每个区域X轴和Y轴的最大值以及最小值即xmax、xmin、ymax、ymin，所以每一个孔洞的直径大小d为X和Y方向的差值的平均值，即孔径大小d：
步骤413、调用Matlab内部函数imhist对所有孔隙进行叠加并绘制孔径分布直方图；
6.根据权利要求1所述的一种基于Matlab分析致密砂岩储层孔隙表征的方法，其特征在于：所述方法二包括：
步骤421、调用Matlab内部程序linecut，将预处理图像进行切片；
步骤422、在每条切线上手动增加控制点，并且设置比例尺，最终将控制点导出成直方图；
步骤423、对切线控制区域内的孔隙进行识别输入，可得到切线上被识别的线段，即为孔径大小；
7.根据权利要求1所述的一种基于Matlab分析致密砂岩储层孔隙表征的方法，其特征在于：所述步骤500中孔洞平均直径及方差的求取是通过分别调用函数mean以及函数var进行求取，即为平均孔隙半径以及方差S。
8.根据权利要求5所述的一种基于Matlab分析致密砂岩储层孔隙表征的方法，其特征在于：所述方法一中孔洞平均直径及方差为：
9.根据权利要求1所述的一种基于Matlab分析致密砂岩储层孔隙表征的方法，其特征在于：所述图像阈值是通过函数imhist建立二值图像直方图,人为的控制阈值大小，把误差控制到最低。
10.根据权利要求1所述的一种基于Matlab分析致密砂岩储层孔隙表征的方法，其特征在于：所述孔隙度与孔径分布分析结果通过应用软件Photoshop与压汞曲线对实验结果准确性进行验证，通过软件Photoshop套索及色彩范围功能对Y488致密储层孔隙度像素进行分析，利用色彩范围功能对孔隙进行选择。
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN201810571899.0A CN108956416A (zh) | 2018-06-06 | 2018-06-06 | 一种基于Matlab分析致密砂岩储层孔隙表征的方法 |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN201810571899.0A CN108956416A (zh) | 2018-06-06 | 2018-06-06 | 一种基于Matlab分析致密砂岩储层孔隙表征的方法 |
Publications (1)
| Publication Number | Publication Date |
|---|---|
| CN108956416A true CN108956416A (zh) | 2018-12-07 |
Family
ID=64493279
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN201810571899.0A Pending CN108956416A (zh) | 2018-06-06 | 2018-06-06 | 一种基于Matlab分析致密砂岩储层孔隙表征的方法 |
Country Status (1)
| Country | Link |
|---|---|
| CN (1) | CN108956416A (zh) |
Cited By (16)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN110644979A (zh) * | 2019-09-03 | 2020-01-03 | 中国石油大学(北京) | 孔隙流体赋存状态的获取方法及其装置 |
| CN110766742A (zh) * | 2019-09-20 | 2020-02-07 | 中国地质大学（武汉） | 基于图像识别软件准确测定白云岩面孔率的方法 |
| CN110927194A (zh) * | 2019-12-11 | 2020-03-27 | 中国科学院地质与地球物理研究所 | 确定泥页岩有机孔含量和孔径分布的方法 |
| CN111485850A (zh) * | 2019-01-09 | 2020-08-04 | 中国石油天然气股份有限公司 | 基于大孔道数据的油井堵水方法及装置 |
| CN111504875A (zh) * | 2020-04-28 | 2020-08-07 | 中国地质大学（北京） | 一种致密砂岩孔喉参数的提取及计算方法及提取装置 |
| CN111505038A (zh) * | 2020-04-28 | 2020-08-07 | 中国地质大学（北京） | 一种基于阴极发光技术定量分析砂岩胶结的实现方法 |
| CN111721683A (zh) * | 2019-03-19 | 2020-09-29 | 中国石油天然气股份有限公司 | 岩样微观孔隙结构定量检测方法及装置 |
| CN111982743A (zh) * | 2020-08-31 | 2020-11-24 | 长春工程学院 | 一种基于物联网的火山岩鉴定方法、系统、终端及介质 |
| CN112051200A (zh) * | 2020-08-26 | 2020-12-08 | 中国石油大学(北京) | 一种致密砂岩储层孔隙结构的定量评价方法和装置 |
| CN113533158A (zh) * | 2021-07-06 | 2021-10-22 | 中国地质大学（北京） | 一种基于sem图像的煤储层孔隙结构参数定量分析方法 |
| CN113570651A (zh) * | 2021-07-06 | 2021-10-29 | 中国地质大学（北京） | 基于sem图像的碳酸盐岩储层孔隙半径分布定量方法 |
| CN113570652A (zh) * | 2021-07-06 | 2021-10-29 | 中国地质大学（北京） | 基于sem图像的砂岩储层矿物晶间孔的定量分析方法 |
| CN113670791A (zh) * | 2021-08-03 | 2021-11-19 | 中国地质大学（北京） | 一种对储层孔隙及颗粒表面积的定量分析方法 |
| CN113870200A (zh) * | 2021-09-16 | 2021-12-31 | 中国地质大学（北京） | 一种定量分析砂岩与碳酸盐岩储层中孔隙配位数的方法 |
| CN114441405A (zh) * | 2021-12-22 | 2022-05-06 | 中国地质大学（北京） | 基于压实和胶结减孔趋势的次生增孔幅度定量评价方法 |
| CN117409408A (zh) * | 2023-12-15 | 2024-01-16 | 北京大学 | 层理缝参数获取方法、装置、设备及可读存储介质 |
Citations (7)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN101639434A (zh) * | 2009-08-27 | 2010-02-03 | 太原理工大学 | 基于显微图像分析固体材料孔隙结构的方法 |
| CN104535475A (zh) * | 2015-01-08 | 2015-04-22 | 中国石油大学(北京) | 碳酸盐岩微观结构的确定方法及装置 |
| US20150153326A1 (en) * | 2013-11-26 | 2015-06-04 | Bergen Teknologioverforing AS, BTO | Quantitative analysis of contact-dependent cell-to-cell transfer and disease transmission |
| CN104933760A (zh) * | 2015-06-18 | 2015-09-23 | 中国地质大学（北京） | 一种重构土壤ct图片三维重建及土壤孔隙搜索方法 |
| CN104931400A (zh) * | 2015-05-29 | 2015-09-23 | 宁夏大学 | 颗粒材料孔隙组构的定量测试与图像分析方法 |
| CN105352873A (zh) * | 2015-11-26 | 2016-02-24 | 中国石油大学(北京) | 页岩孔隙结构的表征方法 |
| CN107133630A (zh) * | 2016-02-29 | 2017-09-05 | 中国石油化工股份有限公司 | 一种基于扫描图像判断碳酸盐岩孔隙类型的方法 |
-
2018
- 2018-06-06 CN CN201810571899.0A patent/CN108956416A/zh active Pending
Patent Citations (7)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN101639434A (zh) * | 2009-08-27 | 2010-02-03 | 太原理工大学 | 基于显微图像分析固体材料孔隙结构的方法 |
| US20150153326A1 (en) * | 2013-11-26 | 2015-06-04 | Bergen Teknologioverforing AS, BTO | Quantitative analysis of contact-dependent cell-to-cell transfer and disease transmission |
| CN104535475A (zh) * | 2015-01-08 | 2015-04-22 | 中国石油大学(北京) | 碳酸盐岩微观结构的确定方法及装置 |
| CN104931400A (zh) * | 2015-05-29 | 2015-09-23 | 宁夏大学 | 颗粒材料孔隙组构的定量测试与图像分析方法 |
| CN104933760A (zh) * | 2015-06-18 | 2015-09-23 | 中国地质大学（北京） | 一种重构土壤ct图片三维重建及土壤孔隙搜索方法 |
| CN105352873A (zh) * | 2015-11-26 | 2016-02-24 | 中国石油大学(北京) | 页岩孔隙结构的表征方法 |
| CN107133630A (zh) * | 2016-02-29 | 2017-09-05 | 中国石油化工股份有限公司 | 一种基于扫描图像判断碳酸盐岩孔隙类型的方法 |
Non-Patent Citations (2)
| Title |
|---|
| SVEN MEISTER: "grain and particle analysis with line intersection method", 《MATHWORKS》 * |
| 雷惊雷 等: "应用MATLAB提取纳米模板特征几何参数", 《功能材料》 * |
Cited By (22)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN111485850A (zh) * | 2019-01-09 | 2020-08-04 | 中国石油天然气股份有限公司 | 基于大孔道数据的油井堵水方法及装置 |
| CN111721683A (zh) * | 2019-03-19 | 2020-09-29 | 中国石油天然气股份有限公司 | 岩样微观孔隙结构定量检测方法及装置 |
| CN110644979A (zh) * | 2019-09-03 | 2020-01-03 | 中国石油大学(北京) | 孔隙流体赋存状态的获取方法及其装置 |
| CN110766742A (zh) * | 2019-09-20 | 2020-02-07 | 中国地质大学（武汉） | 基于图像识别软件准确测定白云岩面孔率的方法 |
| CN110766742B (zh) * | 2019-09-20 | 2022-05-03 | 中国地质大学（武汉） | 基于图像识别软件准确测定白云岩面孔率的方法 |
| CN110927194A (zh) * | 2019-12-11 | 2020-03-27 | 中国科学院地质与地球物理研究所 | 确定泥页岩有机孔含量和孔径分布的方法 |
| CN110927194B (zh) * | 2019-12-11 | 2020-08-18 | 中国科学院地质与地球物理研究所 | 确定泥页岩有机孔含量和孔径分布的方法 |
| CN111504875B (zh) * | 2020-04-28 | 2021-06-04 | 中国地质大学（北京） | 一种致密砂岩孔喉参数的提取及计算方法及提取装置 |
| CN111505038A (zh) * | 2020-04-28 | 2020-08-07 | 中国地质大学（北京） | 一种基于阴极发光技术定量分析砂岩胶结的实现方法 |
| CN111504875A (zh) * | 2020-04-28 | 2020-08-07 | 中国地质大学（北京） | 一种致密砂岩孔喉参数的提取及计算方法及提取装置 |
| CN112051200B (zh) * | 2020-08-26 | 2022-02-15 | 中国石油大学(北京) | 一种致密砂岩储层孔隙结构的定量评价方法和装置 |
| CN112051200A (zh) * | 2020-08-26 | 2020-12-08 | 中国石油大学(北京) | 一种致密砂岩储层孔隙结构的定量评价方法和装置 |
| CN111982743A (zh) * | 2020-08-31 | 2020-11-24 | 长春工程学院 | 一种基于物联网的火山岩鉴定方法、系统、终端及介质 |
| CN113533158A (zh) * | 2021-07-06 | 2021-10-22 | 中国地质大学（北京） | 一种基于sem图像的煤储层孔隙结构参数定量分析方法 |
| CN113570652A (zh) * | 2021-07-06 | 2021-10-29 | 中国地质大学（北京） | 基于sem图像的砂岩储层矿物晶间孔的定量分析方法 |
| CN113570651A (zh) * | 2021-07-06 | 2021-10-29 | 中国地质大学（北京） | 基于sem图像的碳酸盐岩储层孔隙半径分布定量方法 |
| CN113570651B (zh) * | 2021-07-06 | 2023-02-07 | 中国地质大学（北京） | 基于sem图像的碳酸盐岩储层孔隙半径分布定量方法 |
| CN113670791A (zh) * | 2021-08-03 | 2021-11-19 | 中国地质大学（北京） | 一种对储层孔隙及颗粒表面积的定量分析方法 |
| CN113870200A (zh) * | 2021-09-16 | 2021-12-31 | 中国地质大学（北京） | 一种定量分析砂岩与碳酸盐岩储层中孔隙配位数的方法 |
| CN114441405A (zh) * | 2021-12-22 | 2022-05-06 | 中国地质大学（北京） | 基于压实和胶结减孔趋势的次生增孔幅度定量评价方法 |
| CN117409408A (zh) * | 2023-12-15 | 2024-01-16 | 北京大学 | 层理缝参数获取方法、装置、设备及可读存储介质 |
| CN117409408B (zh) * | 2023-12-15 | 2024-03-08 | 北京大学 | 层理缝参数获取方法、装置、设备及可读存储介质 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| CN108956416A (zh) | 一种基于Matlab分析致密砂岩储层孔隙表征的方法 | |
| CN104535475B (zh) | 碳酸盐岩微观结构的确定方法及装置 | |
| US11506650B2 (en) | Method for automatic quantitative statistical distribution characterization of dendrite structures in a full view field of metal materials | |
| CN104619952B (zh) | 具可靠多相渗透性决定数字岩石分析系统及方法 | |
| AU2017239499B2 (en) | Digital rock analysis systems and methods that reliably predict a porosity-permeability trend | |
| CN109087274A (zh) | 基于多维融合及语义分割的电子器件缺陷检测方法及装置 | |
| CN107817199B (zh) | 一种致密砂岩多尺度孔隙模型的构建方法和应用 | |
| CN104819923B (zh) | 基于核磁共振的低渗透砂岩储层孔隙结构定量反演方法 | |
| CN106053315B (zh) | 一种生物碎屑灰岩储层孔隙结构分类方法 | |
| CN104751473B (zh) | 碳酸岩多尺度孔洞特征的确定方法及装置 | |
| CN108802073A (zh) | 一种基于数字岩心的岩石电性参数获取方法及装置 | |
| CN103969168B (zh) | 一种疏松矿物横截面孔隙率的定量测定方法 | |
| CN110021030A (zh) | 一种岩土体材料数字图像的分割阈值确定方法 | |
| CN108734714A (zh) | 一种基于Matlab分析碳酸盐岩结构的方法 | |
| CN110320137A (zh) | 一种基于数字岩心的多尺度融合方法 | |
| CN108956420A (zh) | 一种混凝土孔隙特征提取方法及系统 | |
| CN112051200B (zh) | 一种致密砂岩储层孔隙结构的定量评价方法和装置 | |
| CN108061697B (zh) | 土体三维孔隙率计算方法 | |
| CN113074656A (zh) | 工件孔洞测量方法 | |
| CN110441209A (zh) | 一种基于致密储层数字岩心计算岩石渗透率的方法 | |
| CN110906887B (zh) | 一种服装缝纫平整度检测方法 | |
| Liu et al. | Multicomponent digital core construction and three-dimensional micro-pore structure characterization of shale | |
| CN115221792A (zh) | 一种基于机器学习的天然气水合物类型划分装置及方法 | |
| CN114612578A (zh) | 一种有效储层的识别方法及装置 | |
| WO2024108689A1 (zh) | 一种花斑状碳酸盐岩的储层孔隙度和/或储层渗透率确定方法 |
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: 20181207 |