一种基于Matlab分析致密砂岩储层孔隙表征的方法与流程

本发明涉及致密砂岩储层孔隙表征分析技术领域，具体为一种基于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图像等进行分析，适用于纳米级孔隙。

对于本领域技术人员而言，显然本发明不限于上述示范性实施例的细节，而且在不背离本发明的精神或基本特征的情况下，能够以其他的具体形式实现本发明。因此，无论从哪一点来看，均应将实施例看作是示范性的，而且是非限制性的，本发明的范围由所附权利要求而不是上述说明限定，因此旨在将落在权利要求的等同要件的含义和范围内的所有变化囊括在本发明内。不应将权利要求中的任何附图标记视为限制所涉及的权利要求。