CN108734714B - 一种基于Matlab分析碳酸盐岩结构的方法 - Google Patents
一种基于Matlab分析碳酸盐岩结构的方法 Download PDFInfo
- Publication number
- CN108734714B CN108734714B CN201810571924.5A CN201810571924A CN108734714B CN 108734714 B CN108734714 B CN 108734714B CN 201810571924 A CN201810571924 A CN 201810571924A CN 108734714 B CN108734714 B CN 108734714B
- Authority
- CN
- China
- Prior art keywords
- carbonate rock
- matlab
- pores
- image
- function
- 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.)
- Expired - Fee Related
Links
- 239000011435 rock Substances 0.000 title claims abstract description 70
- BVKZGUZCCUSVTD-UHFFFAOYSA-L Carbonate Chemical compound [O-]C([O-])=O BVKZGUZCCUSVTD-UHFFFAOYSA-L 0.000 title claims abstract description 68
- 238000000034 method Methods 0.000 title claims abstract description 45
- 239000011148 porous material Substances 0.000 claims abstract description 140
- 239000002245 particle Substances 0.000 claims abstract description 114
- 238000005266 casting Methods 0.000 claims abstract description 14
- 239000008187 granular material Substances 0.000 claims description 15
- 238000010586 diagram Methods 0.000 claims description 12
- 235000019738 Limestone Nutrition 0.000 claims description 11
- 239000006028 limestone Substances 0.000 claims description 11
- 238000003384 imaging method Methods 0.000 claims description 10
- 230000008569 process Effects 0.000 claims description 9
- QSHDDOUJBYECFT-UHFFFAOYSA-N mercury Chemical compound [Hg] QSHDDOUJBYECFT-UHFFFAOYSA-N 0.000 claims description 7
- 229910052753 mercury Inorganic materials 0.000 claims description 7
- 238000004458 analytical method Methods 0.000 claims description 5
- 238000004364 calculation method Methods 0.000 claims description 5
- 238000006243 chemical reaction Methods 0.000 claims description 4
- 238000012360 testing method Methods 0.000 claims description 3
- 238000004445 quantitative analysis Methods 0.000 description 12
- 230000015572 biosynthetic process Effects 0.000 description 7
- 238000005755 formation reaction Methods 0.000 description 7
- 239000013078 crystal Substances 0.000 description 6
- 239000000243 solution Substances 0.000 description 6
- 238000012512 characterization method Methods 0.000 description 4
- 239000010459 dolomite Substances 0.000 description 4
- 229910000514 dolomite Inorganic materials 0.000 description 4
- 230000035699 permeability Effects 0.000 description 4
- 230000009466 transformation Effects 0.000 description 4
- 241000234435 Lilium Species 0.000 description 3
- 238000005481 NMR spectroscopy Methods 0.000 description 3
- 239000004568 cement Substances 0.000 description 3
- 238000005516 engineering process Methods 0.000 description 3
- 239000011499 joint compound Substances 0.000 description 3
- 239000008188 pellet Substances 0.000 description 3
- 238000011160 research Methods 0.000 description 3
- 241000195493 Cryptophyta Species 0.000 description 2
- 239000003086 colorant Substances 0.000 description 2
- 238000004891 communication Methods 0.000 description 2
- 238000004512 die casting Methods 0.000 description 2
- 238000002474 experimental method Methods 0.000 description 2
- 239000012530 fluid Substances 0.000 description 2
- 229930195733 hydrocarbon Natural products 0.000 description 2
- 150000002430 hydrocarbons Chemical class 0.000 description 2
- 238000002347 injection Methods 0.000 description 2
- 239000007924 injection Substances 0.000 description 2
- 238000007781 pre-processing Methods 0.000 description 2
- 238000012545 processing Methods 0.000 description 2
- 238000012935 Averaging Methods 0.000 description 1
- 229910021532 Calcite Inorganic materials 0.000 description 1
- 241000206572 Rhodophyta Species 0.000 description 1
- 238000005299 abrasion Methods 0.000 description 1
- 230000009471 action Effects 0.000 description 1
- 230000009286 beneficial effect Effects 0.000 description 1
- -1 biological debris Substances 0.000 description 1
- 125000005587 carbonate group Chemical group 0.000 description 1
- 229910001748 carbonate mineral Inorganic materials 0.000 description 1
- 150000004649 carbonic acid derivatives Chemical class 0.000 description 1
- 230000008859 change Effects 0.000 description 1
- 239000003245 coal Substances 0.000 description 1
- 230000007547 defect Effects 0.000 description 1
- 230000008021 deposition Effects 0.000 description 1
- 238000011161 development Methods 0.000 description 1
- 230000018109 developmental process Effects 0.000 description 1
- 230000000694 effects Effects 0.000 description 1
- 230000003628 erosive effect Effects 0.000 description 1
- 239000010419 fine particle Substances 0.000 description 1
- 230000005484 gravity Effects 0.000 description 1
- 238000010191 image analysis Methods 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
- 230000000877 morphologic effect Effects 0.000 description 1
- 238000011002 quantification Methods 0.000 description 1
- 230000006798 recombination Effects 0.000 description 1
- 238000005215 recombination Methods 0.000 description 1
- 239000004576 sand Substances 0.000 description 1
- 238000001878 scanning electron micrograph Methods 0.000 description 1
- 239000013049 sediment Substances 0.000 description 1
- 230000011218 segmentation Effects 0.000 description 1
- 239000000126 substance Substances 0.000 description 1
- XLYOFNOQVPJJNP-UHFFFAOYSA-N water Substances O XLYOFNOQVPJJNP-UHFFFAOYSA-N 0.000 description 1
Images
Classifications
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
- G06T7/00—Image analysis
- G06T7/10—Segmentation; Edge detection
- G06T7/11—Region-based segmentation
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
- G06T7/00—Image analysis
- G06T7/60—Analysis of geometric attributes
Landscapes
- Engineering & Computer Science (AREA)
- Physics & Mathematics (AREA)
- Computer Vision & Pattern Recognition (AREA)
- General Physics & Mathematics (AREA)
- Theoretical Computer Science (AREA)
- Geometry (AREA)
- Analysing Materials By The Use Of Radiation (AREA)
Abstract
本发明公开了一种基于Matlab分析碳酸盐岩结构的方法,包括如下步骤:首先读取碳酸盐岩的铸体薄片;然后选择色彩模式,将碳酸盐岩的铸体薄片图像二值化;再识别二值图像的黑色像素点,计算孔隙度,同时识别二值图像的白色像素点,分辨颗粒大小,计算粒度;然后求取碳酸盐岩的粒度与孔径大小频率分布图,最后求取碳酸盐岩粒度与孔径的平均直径及方差,并输出二值图像和求取的参数,通过Matlab定量计算不同碳酸盐孔隙结构,平均误差率为+4.661%,其中粒间孔平均误差率为‑4.370%,铸模孔平均误差率为+6.734%,生屑孔误差率为+6.144%,微裂缝平均误差率为‑0.314%,铸模孔与生屑孔的孔隙结构最为适用。
Description
技术领域
本发明涉及致密砂岩储层孔隙表征分析技术领域,具体为一种基于Matlab分析碳酸盐岩结构的方法。
背景技术
致密砂岩储层通常为储层渗透率低的砂岩储层,致密低渗透储层是一个相对概念,世界上并无一固定的标准和界限,它是由不同国家不同时期的资源和技术经济条件和世界油气价格而决定的。
储集层是指具有连通孔隙、允许油气在其中储存和渗滤的岩层。世界上已发现的油气储量大多数来自沉积岩层,其中以砂岩和碳酸盐岩储集层最为重要,裂缝性泥岩和煤层也可作为储集层;火成岩和变质岩储集层中也有工业性油气发现。
储层(Reservoir)的储集能力是由储集层的岩石物理性质决定的,通常包括其孔隙性、渗透性;孔隙性决定了储层储存能力的大小,渗透性决定了储集物的渗流能力。
碳酸盐岩是主要由方解石和白云石等碳酸盐矿物组成的沉积岩,其岩石结构包括颗粒、胶结物、泥晶、孔隙几部分,颗粒可分为内碎屑、鲕粒、球粒、藻粒、团块这几部分,岩石微观孔隙结构是指岩石中孔隙和吼道的几何形状、大小、分布及其相互连通关系。国内外多名学者对碳酸盐岩孔隙结构、表征等进行研究,多种手段被应用于研究致密砂岩孔隙,以核磁共振NMR测井技术、流体注入实验法和图像观测法这三大类为主。基于岩石电阻率参数的核磁共振测井技术,以测定储层孔径分布为主。以压汞法为代表的流体注入实验法能够间接地通过压汞曲线获得储集空间分布、孔径大小等参数,但并不能像以扫描电镜SEM及纳米CT技术等图像观测法直接地观测到储集空间分布、孔喉特征等特征。然而以扫描电镜及纳米CT为实验指导,研究结果大多为对孔喉的几何形状及联通关系进行定性描述,对孔隙表征的定量研究仍有不足。
发明内容
为了克服现有技术方案的不足,本发明提供一种基于Matlab分析碳酸盐岩结构的方法,能有效的解决背景技术提出的问题。
本发明解决其技术问题所采用的技术方案是:
一种基于Matlab分析碳酸盐岩结构的方法,包括如下步骤:
步骤100、读取碳酸盐岩的铸体薄片;
步骤200、选择色彩模式,将碳酸盐岩的铸体薄片图像二值化;
步骤300、识别二值图像的黑色像素点,计算孔隙度,同时识别二值图像的白色像素点,分辨颗粒大小,计算粒度;
步骤400、求取碳酸盐岩的粒度与孔径大小频率分布图;
步骤500、求取碳酸盐岩粒度与孔径的平均直径及方差,并输出二值图像和求取的参数。
作为本发明一种优选的技术方案,所述步骤200中,铸体薄片图像二值化的方法为:将铸体薄片图像导入Matlab的内置Color Threshod应用中通过HSV成像模式进行图像二值化转换。
作为本发明一种优选的技术方案,不同的成像模式适用于不同图像色彩空间,最终目的均为最大限度地识别碳酸盐岩结构并将其转为二值图像,所述HSV模式因其可以调控图像色调H、饱和度S以及明度V。
作为本发明一种优选的技术方案,所述孔隙度的求取方法为:通过Matlab中的size函数对二值图像中的像素进行识别提取,通过函数可以分别获得图像的总像点数Na,再进一步调用能够识别0像素值的函数bwarea获取孔隙像素点数Np,即黑色像素点Np,孔隙度P即黑色像素点与总像素点的比值即:
作为本发明一种优选的技术方案,所述粒度的求取方法包括:
首先,通过Matlab内置函数bwlabel对相同像素值的区域进行识别统计,同时可以确定区域的边界,基于此原理确定孔隙区域,即黑色或白色像素值连通区域;
然后对碳酸盐岩的每一个颗粒以及孔洞调用bwlabel函数,求出第i个区域的X和Y轴的最大值以及最小值即ximax、ximin、yimax、yimin,所以对于第i个颗粒或孔洞的直径大小di为X和Y方向的差值的平均值,即直径大小di:
作为本发明一种优选的技术方案,所述步骤400中,碳酸盐岩的粒度与孔径大小频率分布图通过调用内部函数imhist对所有颗粒及孔洞进行叠加进行绘制。
作为本发明一种优选的技术方案,所述碳酸盐岩粒度与孔径的平均直径及方差的求取方法为:
作为本发明一种优选的技术方案,该方法适用的颗粒包括内碎屑颗粒、鲕粒颗粒、球粒灰岩及云质灰岩。
作为本发明一种优选的技术方案,所述二值图像成像精度的误差主要是:利用函数bwlabel对二值图像中孔隙与颗粒进行识别会产生误差,通过Matlab对60-100个样品的孔隙度进行测量,其中每个样品取三个铸体薄片进行均值计算,最终结果与高压压汞测试结果进行对比,分析误差。
作为本发明一种优选的技术方案,在应用Matlab内置函数bwlabel对粒径大小及孔径大小进行计算时,函数bwlabel会对所有的孔隙及图像中的噪点进行统计,并对该统计过程加以限定调节即最小识别范围,最小识别范围为2μm。
与现有技术相比,本发明的有益效果是:本发明通过Matlab定量计算不同碳酸盐孔隙结构,平均误差率为+4.661%,其中粒间孔平均误差率为-4.370%,铸模孔平均误差率为+6.734%,生屑孔误差率为+6.144%,微裂缝平均误差率为-0.314%;对于不同碳酸盐岩石孔隙而言,通过Matlab定量分析铸模孔与生屑孔的孔隙结构最为适用。
附图说明
图1为本发明Matlab中ColorThreshold色彩模式选择界面图;
图2为本发明Color Threshold功能界面图;
图3为本发明实施方式中岩石结构定量分析相关图像;
图4为本发明的流程图;
图5为本发明读取铸体薄片/SEM图像的方法;
图6为本发明实施方式中对不同类型颗粒的碳酸盐岩定量分析示意图;
图7为本发明实施方式中对不同碳酸盐岩孔隙类型定量分析示意图;
图8为本发明实施方式中不同孔隙类型孔隙度与高压压汞孔隙度示意图;
图9为本发明实施方式中孔隙度实测值与真实值差值比较示意图;
图10为本发明实施方式中不同孔隙类型误差率示意图;
图11为本发明实施方式中RGB、HSV、L*a*b*三种色彩空间对比示意图。
具体实施方式
下面将结合本发明实施例中的附图,对本发明实施例中的技术方案进行清楚、完整地描述,显然,所描述的实施例仅仅是本发明一部分实施例,而不是全部的实施例。基于本发明中的实施例,本领域普通技术人员在没有做出创造性劳动前提下所获得的所有其他实施例,都属于本发明保护的范围。
如图4所示,本发明提供了一种基于Matlab分析碳酸盐岩结构的方法,包括如下步骤:
步骤一、读取碳酸盐岩的铸体薄片;
步骤二、选择色彩模式,将碳酸盐岩的铸体薄片图像二值化;
步骤三、识别二值图像的黑色像素点,计算孔隙度,同时识别二值图像的白色像素点,分辨颗粒大小,计算粒度;
步骤四、求取碳酸盐岩的粒度与孔径大小频率分布图;
步骤五、求取碳酸盐岩粒度与孔径的平均直径及方差,并输出二值图像和求取的参数。
在本实施方式步骤一中,利用铸体薄片图像分析岩石学结构是常见的手段,但是应用软件Matlab对铸体薄片图像分析进行相比现有技术具有突出的特点,能定量计算不同碳酸盐孔隙结构,平均误差率为+4.661%,其中粒间孔平均误差率为-4.370%,铸模孔平均误差率为+6.734%,生屑孔误差率为+6.144%,微裂缝平均误差率为-0.314%,铸模孔与生屑孔的孔隙结构最为适用,粒间孔与微裂缝相对较差。
应用软件Matlab分析处理图像时,需将图像预处理为软件可识别的二值图像,图像二值化就是将图像上的像素点转化为0或255,二值图像中灰度只有两种,非黑即白,也就是将整个图像呈现出明显的黑白效果,以这种方式来操作图像可以更容易识别出图像的结构特征。
因此将铸体薄片图像导入Matlab的内置Color Threshod应用中进行图像二值化转换,共有成像模式,分别为RGB模式、HSV模式、YCbCr模式、L*a*b*模式,如图1所示。
在本实施方式中,不同成像模式适用于不同图像色彩空间,最终目的均为最大限度地识别岩石学结构并将其转为二值图像,HSV模式因其可以调控图像色调H、饱和度S以及明度V,对于呈现褐黄色的碳酸盐颗粒及蓝色充铸孔隙的色彩空间而言,HSV模式能够准确地识别碳酸盐岩图像岩石学结构,所以,步骤二应用HSV模式将图像转化为二值图像,如图2所示。
图2中,左侧图像中黄色为岩石孔隙空间,右侧H为色彩空间,控制图像整个范围,将呈现呈蓝色的孔隙空间从色度盘中去掉,S与V分别为饱和度与明度,稍作调整,最终调整后的整个图像的色彩空间在右下方呈现,空间中呈现孔隙的蓝色基本消失,通过ShowBinary将其转为二值图像输出至Matlab。
通过预处理获取二值图像(如图3c)后,进一步计算图像孔隙度、粒度、孔喉分布、平均孔隙直径等参数。本实施方式以孔隙度、粒度、平均孔隙直径为例。
(1)孔隙度的参数求取方法为:
通过Matlab中的size函数对二值图像中的像素进行识别提取,通过函数可以分别获得图像的总像点数Na(所有黑白像素点总和),再进一步调用能够识别0像素值的函数bwarea获取孔隙像素点数Np,即黑色像素点Np,如图3d,孔隙度P即黑色像素点与总像素点的比值即:
对于岩石粒度以及孔隙度来说,Matlab内置函数bwlabel可对相同像素值的区域进行识别统计,同时可以确定区域的边界,基于此原理确定孔隙区域(黑色或白色像素值连通区域)。
(2)粒度的求取方法:
首先,通过Matlab内置函数bwlabel对相同像素值的区域进行识别统计,同时可以确定区域的边界,基于此原理确定孔隙区域,即黑色或白色像素值连通区域;
然后对碳酸盐岩的每一个颗粒以及孔洞调用bwlabel函数,求出第i个区域的X和Y轴的最大值以及最小值即ximax、ximin、yimax、yimin,所以对于第i个颗粒或孔洞的直径大小di为X和Y方向的差值的平均值,即直径大小di:
(3)平均孔隙直径的求取方法:
在本实施方式中,通过内置函数,内置函数imhist最终可以获得粒度与孔径分布直方图,进一步对孔隙平均值非均质性以及连通性进行分析,通过调用Matlab内置函数mean函数可对平均孔隙半径进行求取。在获取每个颗粒及孔洞的基础上调用内置函数var函数对所有孔洞的方差进行求取,平均孔隙半径及其方差S。
对于地层非均质性来说,方差越大,表现整个地层孔隙分布越复杂,颗粒变化大,非均质性越强。对碳酸盐岩图像二值化预处理后,通过Matlab相关内置函数可以定量求取储层砂岩孔隙度、孔隙半径分布、平均孔隙度、非均质性等表征参数。具体流程如图4所示。
碳酸盐岩的结构组分主要由颗粒、孔隙、泥晶及亮晶胶结物组成,其中碳酸盐岩颗粒包括内碎屑、鲕粒、藻粒、球粒、团块、生物颗粒,本实施方式分别以内碎屑、鲕粒、球粒及生物颗粒这四种颗粒类型为例,定量分析孔径分布。同时,孔隙也由原生孔隙、次生孔隙两种,分别以原生粒间孔、晶间孔、生物骨架孔、铸模孔、收缩孔这五类孔隙为例,定量分析地层孔隙度、孔径分布及大小、平均孔径大小及方差。
本实施方式应用Matlab对内碎屑、鲕粒、球粒及生物颗粒这四种颗粒类型分析,对颗粒含量、颗粒分布、颗粒粒径大小及分布、方差进行定量求解,并研究其方法的适用性。
其中,内碎屑颗粒主要是沉积盆地中沉积不久的汐水流、风暴流、重力流等的作用,破碎半固结或固结的各种碳酸盐沉积物,受波浪、潮搬运、磨蚀、再沉积而成的。以卡沙干油田KE-4井为例,应用图4流程图,首先通过Threshold将铸体薄片转为二值图像,如图6a-2,通过size函数测得地层颗粒含量为81.631%,进一步利用bwlabel函数对颗粒进行识别计算,所识别孔隙如图6a-3,粒径频率分布直方图如图6a-4,可以看出粒径主要集中在50μm左右,颗粒分布相对集中,引用Matlab内置函数mean及var函数求解平均孔隙度及粒径方差。
结果如表1所示,内碎屑粒径平均值为70.891μm,方差为47.520,数值较低,颗粒分选性较强。铸体薄片中大部分颗粒能够被系统所识别,但部分颗粒会被分割为几小部分或一些粒径较小的颗粒被拼合成较大的颗粒,这些误差主要来源于函数bwlabel对颗粒边缘的识别,而控制颗粒边缘主要受所选择的色彩模式及相关参数的影响,不同的色彩模式及参数会对生成的二值图像产生不同的误差。对于内碎屑颗粒,局部颗粒边缘会产生噪点,被系统误识别为颗粒边缘进而分割为几部分,而同时一些点-线接触的颗粒未被识别出,被函数一同识别,最终造成颗粒的拼合。但总体对内碎屑颗粒而言,误差较小,适用性较高。
鲕粒颗粒是具有核心和同心层结构的球状颗粒,根据鲕粒的结构和形态特征,可把鲕粒划分为以下类型:正常鲕、表皮鲕、复鲕、椭形鲕、放射鲕及单晶鲕这几种类型。本实施方式选取塔里木盆地塔中4井亮晶鲕粒灰岩为例,对其颗粒含量、颗粒分布、颗粒粒径大小及分布、方差进行定量分析。对二值图像进行读取分析,最终颗粒分布图如图6b-3所示,与6b-1图像相比,大部分鲕粒被还原较好,但部分鲕粒被Matlab错误识别为两个以上的颗粒,这会导致粒径计算结果与实际值偏低,平均粒径大小偏低,方差增大,利用Matlab所计算粒径分布图结果如图6b-4,所测的颗粒含量、平均粒径及粒径方差如表1所示,由于鲕粒形状非常规则,系统识别度非常高。
球粒是较细粒的(粗粉砂级或砂级)、由灰泥组成的、不具特殊内部结构的、球形或卵形的、分选较好的颗粒。以塔里木盆地中4井亮晶-生屑球粒灰岩为例,应用图4程序流程图对球粒灰岩进行定量分析,二值图像见图6c-2,与图6c-1相比,部分海百合颗粒由于颜色与亮晶胶结物相似,被Matlab误认为颗粒以外的物质,导致平均粒径偏小,但由于海百合粒径与球粒粒径相差较多,最终粒径方差将比理论值偏低,分选性也会变高。球粒颗粒分布图见图6c-3,图中部分球粒颗粒被系统分割为几小部分,这也是误差来源之一,适用性较高。
生物颗粒以以鄂尔多斯盆地延628井奥陶系云质灰岩为例,地层晶间孔发育,颗粒属于自形-半自形结构,由Matlab内置函数bwlabel对二值图像进行分析,得到图6d-3颗粒分布图,图中左上方部分晶间孔不发育的区域,系统未识别出细小的晶间孔,一些颗粒相联通,颗粒粒径与原铸体薄片相比偏高,利用内置函数imhist得到孔径频率分布直方图,由mean函数计算得到平均粒径为88.752μm。
表1:图6中不同类型地层的颗粒比例/%、平均粒径/μm及方差
如图5所示,应用Matlab对不同颗粒结构定量分析,在整个流程中,会产生误差,包括选择色彩模式、铸体图像的二值转化、颗粒被分割或重组等,每一步地操作都会对下一步产生影响,误差在整个过程中会叠加,因此对于图像的预处理非常重要,包括色彩模式的选择以及二值图像的转化。
对于不同碳酸盐岩孔隙类型定量分析如下:
碳酸盐岩孔隙类型分为原生孔隙、次生孔隙、裂缝三部分,包括原生粒间孔、原生晶间孔、铸模孔、生物孔、微裂缝等几个部分,主要的碳酸盐岩的颗粒类型包括原生粒间孔,所以下主要讨论晶间孔、铸模孔、生屑孔、微裂缝的定量分析及程序适用性。
(1)晶间孔
以鄂尔多斯盆地东缘延1757井奥陶系白云岩晶间孔为例,白云岩晶体自形程度较好,利用图5方法处理过程对白云岩晶间孔进行定量分析,图7a-2为Matlab生成的二值图像,与铸体薄片相比,大部分孔隙被识别出,但同时一些附着在薄片的黑色脏污也被误认为孔隙,这是误差来源之一。此外,一些孔洞也被系统分割成几部分或拼合几部分,这也是另一部误差来源。但对晶间孔的识别非常准确,相关定量分析结果如表2所示。
(2)铸模孔
样品取自卡莎干油田KE-4井,颗粒灰岩早期的强烈的选择性溶蚀形成了大量的铸模孔,如图7b-1所示,铸模孔形状较为规则,利用Matlab进行识别并利用相关内置函数进行定量分析后,孔隙重合度较高,除极少数孔隙别分割为几部分。
(3)生屑骨架孔
铸体薄片泥晶颗粒灰岩取样于卡沙干油田KED04井,颗粒主要为鲕粒、生物碎屑、砂屑等,生物颗粒见有孔虫、红藻、棘皮类,有孔虫体腔内的软体发生腐烂之后形成生物体腔孔如图7c-1,利用Matlab对铸体薄片二值化,孔隙包括有孔虫体腔内孔、粒间孔都能够被识别出,但一些孔隙会被分解为几小部分,影响孔径分布及平均孔隙度,孔径分布相对中,孔隙度为15.37%,程序非常适用。
(4)微裂缝-缝合线
以鄂尔多斯盆地东缘延1758井奥陶系泥晶灰岩为例,Matlab能够很好识别出微裂缝,如图7d-2,利用size函数及公式能够便捷的求出孔隙度,但利用函数bwlalel对孔隙识别时,泥晶中一些细小的微裂缝未被识别出,同时大量微裂缝被分为细小的颗粒,所以应用Matlab定量分析微裂缝-缝合线时,只适用于定量分析孔隙度,对孔径分布、大小并不适用。
表2:图7中不同地层孔隙类型的孔隙度/%、平均粒径/μm及方差
不同孔隙类型 | 晶间孔 | 铸模孔 | 生屑骨架孔 | 微裂缝+缝合线 |
孔隙度/% | 9.13 | 15.37 | 8.1675 | 11.094 |
平均孔径/μm | 33.7323 | 49.1569 | 24.1185 | 22.854 |
孔径方差 | 23.7345 | 36.3214 | 16.6671 | 16.1559 |
关于二值图像的成像精度,首先,在整个流程中,从铸体图像转化为二值图像,到利用函数bwlabel对二值图像中孔隙与颗粒进行识别均会产生一定的误差,本实施方式通过Matlab对60个样品的孔隙度进行测量,其中每个样品取三个铸体薄片进行均值,最终结果与高压压汞测试结果进行对比,分析误差,结果见图8。
图8中,红色代表利用Matlab的测量结果,蓝色代表高压压汞孔隙度,四种被圈定的孔隙类型分别为粒间孔、铸模孔、生屑孔及微裂缝。测量值与真实值存在一定偏差,不同类型的偏差值也不一样,偏差结果如图9所示。
由图9可以看出,铸模孔与生屑孔两种孔隙类型测量值大部分大于真实值,这主要是在图像二值化过程中,由于色彩空间模式的选择及参数调控所产生的噪点引起的,噪点被Matlab误认为孔隙,由此造成孔隙度比真实值有所增加。粒间孔与微裂缝两种孔隙类型的一些颗粒同铸模孔与生屑孔一样,实际测量值有所偏大,但多数测量值低于理论值,这主要归根于粒间孔与微裂缝自身特征,由于一些粒间孔与微裂缝孔隙较窄,半径大小为纳米级孔隙,并未被Matlab识别,因此形成该种现象,具体误差率见图10。
大部分孔隙误差集中在-5%~10%之间,平均误差率为+4.661%,其中粒间孔平均误差率为-4.370%,标准差为8.818,铸模孔平均误差率为+6.734%,标准差为6.913,生屑孔误差率为+6.144%,标准差为6.731,微裂缝平均误差率为-0.314%,标准差为+9.271。所以对于不同碳酸盐岩石孔隙而言,通过Matlab定量分析铸模孔与生屑孔的孔隙结构最为适用,粒间孔与微裂缝相对较差。
在本实施方式中,图像二值化的精准度直接影响颗粒与孔隙分布,对定量分析产生最直接的影响。而图像二值化包括色彩空间的选择及参数的调控,Color Threshold中共有四种色彩空间模式RGB、HSV、YCbCr、L*a*b*,不同色彩空间模式适用于不同的色域的铸体薄片,以最常见的“黄蓝色”铸体薄片为例(黄色为碳酸盐岩干涉色,蓝色为孔隙被染色),分别以RGB、HSV、L*a*b*三种常规色彩空间模式对图像进行二值化处理,见图11。
对于RGB色彩模式而言,r、g、b三个变量并不能很好的消除色彩空间中的蓝色,调节变量,将会对色彩空间进行线性变换,导致图像中杂基、暗色矿物产生大量噪点,L*a*b*模式中L表示亮度(Luminosity),a表示从洋红色至绿色的范围,b表示从黄色至蓝色的范围,虽然控制参数b可以很好地消除图像中的蓝色,但由于同样是线性变换,一些代表颗粒的黄色也会受到波及,导致颗粒中产生大量噪点。对于HSV模式,H控制的色度盘可以轻易地将蓝色去除,为空间中的非线性变换,并不会影响其它色彩,基本无噪点产生,因此也最适用于碳酸盐岩铸体薄片。
在本实施方式中,应用Matlab对碳酸盐岩石结构定量分析的整个流程中,图像噪点是对计算结果最大的影响,因此去除噪点是必要的一部,应用Matlab内置函数bwlabel对粒径大小及孔径大小进行计算时,函数会对所有的孔隙及图像中的噪点进行统计,因此对于这个统计过程加以限定调节即最小识别范围,最小识别范围是根据研究图像而定的,本实施方式中的限定界限为2μm,但在对噪点去除的同时,一些微裂缝及较小的粒间孔也会被剔除,因此,对于粒径较低的碳酸盐岩如泥晶灰岩及微裂缝发育的碳酸盐岩并不适用。
在本实施方式中,通过Matlab定量计算不同碳酸盐孔隙结构,平均误差率为+4.661%,其中粒间孔平均误差率为-4.370%,铸模孔平均误差率为+6.734%,生屑孔误差率为+6.144%,微裂缝平均误差率为-0.314%。
对于不同碳酸盐岩石孔隙而言,通过Matlab定量分析铸模孔与生屑孔的孔隙结构最为适用,粒间孔与微裂缝相对较差。对于粒径较低的碳酸盐岩如泥晶灰岩及微裂缝发育的碳酸盐岩并不适用。
对于本领域技术人员而言,显然本发明不限于上述示范性实施例的细节,而且在不背离本发明的精神或基本特征的情况下,能够以其他的具体形式实现本发明。因此,无论从哪一点来看,均应将实施例看作是示范性的,而且是非限制性的,本发明的范围由所附权利要求而不是上述说明限定,因此旨在将落在权利要求的等同要件的含义和范围内的所有变化囊括在本发明内。不应将权利要求中的任何附图标记视为限制所涉及的权利要求。
Claims (8)
1.一种基于Matlab分析碳酸盐岩结构的方法,其特征在于,包括如下步骤:
步骤100、读取碳酸盐岩的铸体薄片;
步骤200、选择色彩模式,将碳酸盐岩的铸体薄片图像二值化;
步骤300、识别二值图像的黑色像素点,计算孔隙度,同时识别二值图像的白色像素点,分辨颗粒大小,计算粒度;
步骤400、求取碳酸盐岩的粒度与孔径大小频率分布图;
步骤500、求取碳酸盐岩粒度与孔径的平均直径及方差,并输出二值图像和求取的参数;
所述二值图像成像精度的误差分析方法包括:
利用函数bwlabel对二值图像中孔隙与颗粒进行识别会产生误差,通过Matlab对60-100个样品的孔隙度进行测量,其中每个样品取三个铸体薄片进行均值计算,最终结果与高压压汞测试结果进行对比,分析误差;
在应用Matlab内置函数bwlabel对粒径大小及孔径大小进行计算时,函数bwlabel会对所有的孔隙及图像中的噪点进行统计,并对该统计过程加以限定调节即最小识别范围,最小识别范围为2μm。
2.根据权利要求1所述的一种基于Matlab分析碳酸盐岩结构的方法,其特征在于:所述步骤200中,铸体薄片图像二值化的方法为:将铸体薄片图像导入Matlab的内置ColorThreshod应用中通过HSV成像模式进行图像二值化转换。
3.根据权利要求2所述的一种基于Matlab分析碳酸盐岩结构的方法,其特征在于:不同的成像模式适用于不同图像色彩空间,最终目的均为最大限度地识别碳酸盐岩结构并将其转为二值图像,所述HSV成像模式因其可以调控图像色调H、饱和度S以及明度V。
6.根据权利要求1所述的一种基于Matlab分析碳酸盐岩结构的方法,其特征在于:所述步骤400中,碳酸盐岩的粒度与孔径大小频率分布图通过调用内部函数imhist对所有颗粒及孔洞进行叠加进行绘制。
8.根据权利要求1所述的一种基于Matlab分析碳酸盐岩结构的方法,其特征在于:该方法适用的颗粒包括内碎屑颗粒、鲕粒颗粒、球粒灰岩及云质灰岩。
Priority Applications (1)
Application Number | Priority Date | Filing Date | Title |
---|---|---|---|
CN201810571924.5A CN108734714B (zh) | 2018-06-06 | 2018-06-06 | 一种基于Matlab分析碳酸盐岩结构的方法 |
Applications Claiming Priority (1)
Application Number | Priority Date | Filing Date | Title |
---|---|---|---|
CN201810571924.5A CN108734714B (zh) | 2018-06-06 | 2018-06-06 | 一种基于Matlab分析碳酸盐岩结构的方法 |
Publications (2)
Publication Number | Publication Date |
---|---|
CN108734714A CN108734714A (zh) | 2018-11-02 |
CN108734714B true CN108734714B (zh) | 2022-11-25 |
Family
ID=63932055
Family Applications (1)
Application Number | Title | Priority Date | Filing Date |
---|---|---|---|
CN201810571924.5A Expired - Fee Related CN108734714B (zh) | 2018-06-06 | 2018-06-06 | 一种基于Matlab分析碳酸盐岩结构的方法 |
Country Status (1)
Country | Link |
---|---|
CN (1) | CN108734714B (zh) |
Families Citing this family (7)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
CN109856029B (zh) * | 2019-02-01 | 2021-07-30 | 中海石油(中国)有限公司上海分公司 | 一种基于图像分析的孔隙度评价方法 |
CN111504875B (zh) * | 2020-04-28 | 2021-06-04 | 中国地质大学(北京) | 一种致密砂岩孔喉参数的提取及计算方法及提取装置 |
CN111505038A (zh) * | 2020-04-28 | 2020-08-07 | 中国地质大学(北京) | 一种基于阴极发光技术定量分析砂岩胶结的实现方法 |
CN111413284B (zh) * | 2020-05-18 | 2022-03-22 | 攀钢集团攀枝花钢铁研究院有限公司 | 钒钛烧结矿中主要物相的定量检测方法 |
CN112433069A (zh) * | 2020-12-01 | 2021-03-02 | 安徽理工大学 | 一种奥陶系灰岩顶部相对隔水层的微观结构特征判别方法 |
CN113670791A (zh) * | 2021-08-03 | 2021-11-19 | 中国地质大学(北京) | 一种对储层孔隙及颗粒表面积的定量分析方法 |
CN113870200B (zh) * | 2021-09-16 | 2022-05-31 | 中国地质大学(北京) | 一种定量分析砂岩与碳酸盐岩储层中孔隙配位数的方法 |
Citations (2)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
CN102749046A (zh) * | 2012-07-23 | 2012-10-24 | 中国地质大学(武汉) | 岩体结构面直剪试验中剪切面积的测量方法 |
CN106896213A (zh) * | 2017-02-22 | 2017-06-27 | 中国地质大学(武汉) | 一种基于点云数据的岩体结构面智能识别与信息提取方法 |
Family Cites Families (2)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
US6937744B1 (en) * | 2000-06-13 | 2005-08-30 | Microsoft Corporation | System and process for bootstrap initialization of nonparametric color models |
US20160203591A1 (en) * | 2015-01-09 | 2016-07-14 | Umm Al-Qura University | System and process for monitoring the quality of food in a refrigerator |
-
2018
- 2018-06-06 CN CN201810571924.5A patent/CN108734714B/zh not_active Expired - Fee Related
Patent Citations (2)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
CN102749046A (zh) * | 2012-07-23 | 2012-10-24 | 中国地质大学(武汉) | 岩体结构面直剪试验中剪切面积的测量方法 |
CN106896213A (zh) * | 2017-02-22 | 2017-06-27 | 中国地质大学(武汉) | 一种基于点云数据的岩体结构面智能识别与信息提取方法 |
Non-Patent Citations (2)
Title |
---|
分析碳酸盐岩孔隙系统数字图像的新方法;赵永刚 等.;《计算机应用研究》;20061001;第169-171页 * |
应用MATLAB提取纳米模板特征几何参数;雷惊雷;《功能材料》;20111020;第1893-1897页 * |
Also Published As
Publication number | Publication date |
---|---|
CN108734714A (zh) | 2018-11-02 |
Similar Documents
Publication | Publication Date | Title |
---|---|---|
CN108734714B (zh) | 一种基于Matlab分析碳酸盐岩结构的方法 | |
CN110458119B (zh) | 一种非接触式测量的混凝土骨料级配快速识别方法 | |
CN107228816B (zh) | 一种泥页岩中不同类型孔隙演化评价方法 | |
CN107133630B (zh) | 一种基于扫描图像判断碳酸盐岩孔隙类型的方法 | |
CN110264466A (zh) | 一种基于深度卷积神经网络的钢筋检测方法 | |
CN108961246B (zh) | 一种基于人工智能的扫描电镜图像孔隙识别方法 | |
CN102494976B (zh) | 一种超细晶粒钢晶粒的自动测量及其形态分类方法 | |
Quan et al. | The method of the road surface crack detection by the improved Otsu threshold | |
CN104237103A (zh) | 一种孔隙连通性定量表征方法及装置 | |
CN114067114B (zh) | 一种基于藻类细胞形态的面积核分割计数方法 | |
CN112150430B (zh) | 一种利用岩石细观结构数字图像的数值分析方法 | |
CN107314957B (zh) | 一种岩石块度尺寸分布的测量方法 | |
CN112051200B (zh) | 一种致密砂岩储层孔隙结构的定量评价方法和装置 | |
CN104268600B (zh) | 一种基于Minkowski距离的矿物浮选泡沫图像纹理分析及工况识别方法 | |
YANG et al. | Bubble size estimation using interfacial morphological information for mineral flotation process monitoring | |
CN110443793A (zh) | 一种沥青混合料空隙分布均匀性评价方法 | |
CN104318564A (zh) | 一种矿物颗粒分相的方法 | |
CN103971367B (zh) | 水文资料图像分割方法 | |
CN108986077A (zh) | 基于双树复小波域共生增广矩阵的浮选泡沫工况识别方法 | |
CN107220946B (zh) | 一种岩石运输带上不良块度图像的实时剔除方法 | |
CN113570652B (zh) | 基于sem图像的砂岩储层矿物晶间孔的定量分析方法 | |
CN112164029B (zh) | 一种用于加气混凝土孔结构的图像分析方法 | |
CN107977968A (zh) | 基于建筑物阴影信息挖掘的建筑物分层检测方法 | |
CN115359003A (zh) | 两步式隧道灰度图像的裂缝识别方法、系统、介质及设备 | |
CN115640546A (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 | ||
GR01 | Patent grant | ||
GR01 | Patent grant | ||
CF01 | Termination of patent right due to non-payment of annual fee |
Granted publication date: 20221125 |