CN113570651A - 基于sem图像的碳酸盐岩储层孔隙半径分布定量方法 - Google Patents

基于sem图像的碳酸盐岩储层孔隙半径分布定量方法 Download PDF

Info

Publication number
CN113570651A
CN113570651A CN202110760736.9A CN202110760736A CN113570651A CN 113570651 A CN113570651 A CN 113570651A CN 202110760736 A CN202110760736 A CN 202110760736A CN 113570651 A CN113570651 A CN 113570651A
Authority
CN
China
Prior art keywords
pixel
pore
image
value
point
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.)
Granted
Application number
CN202110760736.9A
Other languages
English (en)
Other versions
CN113570651B (zh
Inventor
久博
黄文辉
郝睿林
于春兰
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
China University of Geosciences Beijing
Original Assignee
China University of Geosciences Beijing
Priority date (The priority date 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 date listed.)
Filing date
Publication date
Application filed by China University of Geosciences Beijing filed Critical China University of Geosciences Beijing
Priority to CN202110760736.9A priority Critical patent/CN113570651B/zh
Publication of CN113570651A publication Critical patent/CN113570651A/zh
Application granted granted Critical
Publication of CN113570651B publication Critical patent/CN113570651B/zh
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Images

Classifications

    • 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
    • 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
    • G01N23/00—Investigating or analysing materials by the use of wave or particle radiation, e.g. X-rays or neutrons, not covered by groups G01N3/00 – G01N17/00, G01N21/00 or G01N22/00
    • G01N23/22—Investigating or analysing materials by the use of wave or particle radiation, e.g. X-rays or neutrons, not covered by groups G01N3/00 – G01N17/00, G01N21/00 or G01N22/00 by measuring secondary emission from the material
    • G01N23/225—Investigating or analysing materials by the use of wave or particle radiation, e.g. X-rays or neutrons, not covered by groups G01N3/00 – G01N17/00, G01N21/00 or G01N22/00 by measuring secondary emission from the material using electron or ion
    • G01N23/2251—Investigating or analysing materials by the use of wave or particle radiation, e.g. X-rays or neutrons, not covered by groups G01N3/00 – G01N17/00, G01N21/00 or G01N22/00 by measuring secondary emission from the material using electron or ion using incident electron beams, e.g. scanning electron microscopy [SEM]
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T5/00—Image enhancement or restoration
    • G06T5/70—Denoising; Smoothing
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T5/00—Image enhancement or restoration
    • G06T5/77—Retouching; Inpainting; Scratch removal
    • 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/13—Edge detection
    • 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/136—Segmentation; Edge detection involving thresholding
    • 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/181—Segmentation; Edge detection involving edge growing; involving edge linking
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T2207/00—Indexing scheme for image analysis or image enhancement
    • G06T2207/10—Image acquisition modality
    • G06T2207/10056—Microscopic image
    • G06T2207/10061—Microscopic image from scanning electron microscope
    • Y—GENERAL TAGGING OF NEW TECHNOLOGICAL DEVELOPMENTS; GENERAL TAGGING OF CROSS-SECTIONAL TECHNOLOGIES SPANNING OVER SEVERAL SECTIONS OF THE IPC; TECHNICAL SUBJECTS COVERED BY FORMER USPC CROSS-REFERENCE ART COLLECTIONS [XRACs] AND DIGESTS
    • Y02—TECHNOLOGIES OR APPLICATIONS FOR MITIGATION OR ADAPTATION AGAINST CLIMATE CHANGE
    • Y02A—TECHNOLOGIES FOR ADAPTATION TO CLIMATE CHANGE
    • Y02A90/00—Technologies having an indirect contribution to adaptation to climate change
    • Y02A90/30—Assessment of water resources

Landscapes

  • Engineering & Computer Science (AREA)
  • Physics & Mathematics (AREA)
  • General Physics & Mathematics (AREA)
  • Theoretical Computer Science (AREA)
  • Computer Vision & Pattern Recognition (AREA)
  • Chemical & Material Sciences (AREA)
  • Life Sciences & Earth Sciences (AREA)
  • Health & Medical Sciences (AREA)
  • Analytical Chemistry (AREA)
  • Biochemistry (AREA)
  • General Health & Medical Sciences (AREA)
  • Immunology (AREA)
  • Pathology (AREA)
  • Dispersion Chemistry (AREA)
  • Geometry (AREA)
  • Analysing Materials By The Use Of Radiation (AREA)
  • Image Processing (AREA)

Abstract

本发明公开了一种基于SEM图像的碳酸盐岩储层孔隙半径分布定量方法，包括以下步骤：获取每种碳酸盐岩储层在电镜扫描下的SEM图像，SEM图像内按照自适应阈值方式对每种碳酸盐岩储层二值化处理并转化为二值图像；消除将非孔隙特征的像素点，并将两个孔隙单元连接以修复噪点；重新遍历完成噪点消除和噪点修复的二值图像，将像素值突变的像素点作为孔隙边界，并且基于孔隙边缘划分二值图像中的孔隙空间个体；确定每个像素点在二维坐标系内的坐标值计算每个孔隙空间个体的孔隙半径，并输出孔隙半径分布直方图和平均孔隙半径；本发明基于常规扫描电镜图像，试验成本低，且明够对多种碳酸盐岩储层半径进行定量分析，并且能够提供可视化结果。

Description

基于SEM图像的碳酸盐岩储层孔隙半径分布定量方法
技术领域
本发明涉及碳酸盐岩储层孔隙分析技术领域，具体涉及一种基于SEM图像的碳酸盐岩储层孔隙半径分布定量方法。
背景技术
自形程度指组成岩石的矿物的形态特点，根据矿物自形程度可以分为三种结构，分别为自形粒状结构(矿物颗粒按自已的结晶习性发育成被规则的晶面所包围的自形晶)；它形粒状结构(矿物颗粒多呈不规则的形态－它形晶，少有完整规则的晶面)；半自形粒状结构(矿物颗粒按结晶习性发育一部分规则的晶面，其它的晶面发育不好呈不规则的形态)。受颗粒自形程度影响，对于半自形或者它形的碳酸盐岩储层颗粒间孔隙非均质非常强，孔隙间隔非常低，一般为5微米以下，整个岩体较为致密，对应的孔隙半径的定量分析较为困难，常压压汞法与氦气吸附法并不适用。
目前对于这些自形与它形颗粒发育的碳酸盐岩致密储层的孔隙半径的计算方法主要由高压压汞法与3D-CT数字岩芯成像技术。高压压汞法不同于常规的恒速压汞法，因其注入压力非常大，使得注入流体能过够进入一些微孔(小于1微米)与小孔(1-5微米)间，从而获得致密储层的孔隙分布。3D-CT数字岩芯成像技术利用高分辨率X射线的微计算层析成像设备(微CT)技术对岩芯孔隙中纳米级、微米级和毫米级的孔隙进行精准地提取，也是目前致密或者半自形-它形碳酸盐岩储层孔隙分布的定量提取手段精准度最高的方法，寻求一个成本较低，准确度较高的，针对孔隙较为致密的半自形-它形的碳酸盐岩孔隙分析方法对于碳酸盐岩储层孔隙半径分布的研究十分重要。
但是上述两种方法还存在的缺陷如下：
(1)高压压汞法能够对小于5微米以下的致密碳酸盐岩储层孔隙分布级大小进行定量分析，但再高压流体注入后，岩体内压徒增，压力再部分脆弱掩体初释放，会导致一些新的孔隙产生。同时其只能间接测量，无法可视化观察；
(2)3D-CT数字岩芯成像技术测试结果，费时，整个拍摄过程的使用成本非常高。
发明内容
本发明的目的在于提供一种基于SEM图像的碳酸盐岩储层孔隙半径分布定量方法，以解决现有技术中高压汞法引起新的孔隙，造成定量分析误差大，且分析成本高的技术问题。
为解决上述技术问题，本发明具体提供下述技术方案：
一种基于SEM图像的碳酸盐岩储层孔隙半径分布定量方法，包括以下步骤：
步骤100、准备各类自形程度的碳酸盐岩储层，获取每种碳酸盐岩储层在电镜扫描下的SEM图像，所述SEM图像内按照自适应阈值方式对每种碳酸盐岩储层二值化处理并转化为二值图像；
步骤200、将像素值突变的像素点围成的小面积作为非孔隙特征进行消除，计算像素值突变的像素点围成的大面积之间的邻近域包含的像素点以将至少两个大面积连接修复；
步骤300、重新遍历完成噪点消除和噪点修复的所述二值图像，将像素值突变的像素点作为孔隙边界，并且基于孔隙边缘划分所述二值图像中的孔隙空间个体；
步骤400、确定每个像素点在所述二维坐标系内的坐标值，基于坐标值计算所述二值图像中的每个所述孔隙空间个体的孔隙半径，并输出每种所述碳酸盐岩储层的孔隙半径分布直方图和平均孔隙半径。
作为本发明的一种优选方案，在步骤100中，将所述SEM图像转化二值图像的具体实现步骤为：
将所述SEM图像转化为矩阵图像，并对所述矩阵图像平均切割为多个图像碎片；
抽样采集每个所述图像碎片中多个像素点的像素值，且确定对每个所述图像碎片二值化处理对应的分界像素值；
将每个所述图像碎片中低于分界像素值的像素点的像素值设置为0，且将高于分界像素值的像素点的像素值设置为255；
将所有的图像碎片按照行列式对齐，形成二值化处理后的二值图像。
作为本发明的一种优选方案，确定分界像素值的实现方式为：
将抽样采集的列间距和行间距形成抽样网格，并获取所述抽样网格的每个交叉点对应的像素值；
求取所有像素点的平均像素值，且将所述平均像素值作为分界像素值。
作为本发明的一种优选方案，在步骤200中，建立关于所述二值图像的二维坐标系，按照矩阵遍历方式识别所述二值图像中所有像素点对应的像素值，标注像素值突变的像素点，并根据像素值突变的像素点围成的面积选择消除噪点或修复噪点，消除噪点的实现方法为：
以所述二值图像的垂直交叉边缘为坐标轴建立二维坐标系，且按照先行后列的方式遍历所述二值图像以确定每个所述像素点的像素值，选定行列遍历时的遍历间距均为一个像素，确定出每行中像素值突变的像素点，且将像素值突变的像素点保存在一个集合内；
将相邻行的像素值突变的像素点统计在同一个集合内，并且判断每个集合中包含的像素值突变的像素点个数；
将个数少于等于设定值的集合定义为噪点，并重新设定该像素点的像素值为255以消除该噪点。
作为本发明的一种优选方案，修复噪点的实现方式为：
将个数多于设定值的集合定义为待修复噪点，对两个相邻的所述集合进行多次方数字形态学运算，根据运算结果将两个相邻的集合之间的像素点的像素值重置为0以进行孔隙连接。
作为本发明的一种优选方案，对两个相邻的所述集合进行多次方数字形态学运算具体为：计算两个相邻的集合中同行或者同列像素点之间的邻域像素，并在所述邻域像素小于设定值时，将像素值为255的像素点设置为0以将两个所述像素点连接。
作为本发明的一种优选方案，在步骤300中，识别每种所述碳酸盐岩储层对应所述二值图像中的孔隙边缘和孔隙空间个体的实现方法为：
步骤301、将每个集合内的所述像素点的行列坐标转化为树状图，且将同列的像素点设置在所述树状图同一层的节点上；
步骤302、在所述二值图像内重新定义对应集合内的所述像素点的RGB值；
步骤303、将所有重新定义RGB值的像素点首尾连接构成的曲线设定为所述孔隙边缘，且将所有重新定义RGB值的像素点包围的区域设定为孔隙空间个体。
作为本发明的一种优选方案，在步骤302中，确定遍历过程中的每一行的第一个像素值突变的所述像素点，按照遍历顺序划分每一行的像素值突变的所述像素点的优先级，并且将不同行的像素值突变的所述像素点的优先级设定为与遍历顺序一一对应，并且将同一行的像素值突变的所述像素点设置为同一个优先级；
将仅有一个像素值突变的所述像素点按照优先级顺序分别设置为一个所述孔隙空间个体的起点和终点，按照从起点到终点的顺序分为两条测绘曲线，两条所述测绘曲线顺次将不同优先级的所述像素值突变的所述像素点依次连接形成一个所述孔隙空间个体。
作为本发明的一种优选方案，在步骤400中，所述孔隙半径为单个孔隙单位的长轴与短轴的平均值的一半，所述二值图像中每个所述孔隙空间个体的所述孔隙半径利用每个所述孔隙空间个体的像素点在所述二维坐标系的最大X坐标，最小X坐标，最大Y坐标，最小Y坐标进行求取，具体的实现步骤为：
选择每个所述孔隙空间个体内像素值突变的所述像素点的横坐标最大值对应的像素点、横坐标最小值对应的像素点，纵坐标最大值的像素点和坐标表最小值的像素点；
每个所述孔隙空间个体的长轴半径为Max(Xmax-Xmin，Ymax-Ymin)，而每个所述孔隙空间个体的短轴半径为Min(Xmax-Xmin，Ymax-Ymin)；
将通过公式[(Xmax-Xmin)+(Ymax-Ymin)]/2计算每个所述孔隙空间个体的孔隙半径。
作为本发明的一种优选方案，在步骤400中，统计每种碳酸盐岩储层中的每个孔隙空间个体对应的所述孔隙半径，并确定所述孔隙半径相同的所述孔隙空间个体的数量，创建关于孔隙半径和孔隙空间个体数量的孔隙半径直方图，并计算每种碳酸盐岩储层样品的孔隙半径平均值。
本发明与现有技术相比较具有如下有益效果：
(1)本发明能够对小于5微米以下的自形、半自形、它形的碳酸盐岩储层半径进行定量分析，并且能够提供可视化结果，已解决现有技术中高压压汞法无法实现直视化测试，且引起新的孔隙造成测量误差的问题；
(2)本发明基于常规扫描电镜图像，并结合计算机算法，成分非常低。
附图说明
为了更清楚地说明本发明的实施方式或现有技术中的技术方案，下面将对实施方式或现有技术描述中所需要使用的附图作简单地介绍。显而易见地，下面描述中的附图仅仅是示例性的，对于本领域普通技术人员来讲，在不付出创造性劳动的前提下，还可以根据提供的附图引伸获得其它的实施附图。
图1为本发明实施例提供的碳酸盐岩孔隙半径定量分析的流程示意图；
图2为本发明实施例提供的不同自形的碳酸盐岩储层的SEM图像；
图3为本发明实施例提供的不同自形的碳酸盐岩储层的未降噪二值图像；
图4为本发明实施例提供的不同自形的碳酸盐岩储层的已降噪二值图像；
图5为本发明实施例提供的不同自形的碳酸盐岩储层的孔隙半径分布直方图。
具体实施方式
下面将结合本发明实施例中的附图，对本发明实施例中的技术方案进行清楚、完整地描述，显然，所描述的实施例仅仅是本发明一部分实施例，而不是全部的实施例。基于本发明中的实施例，本领域普通技术人员在没有做出创造性劳动前提下所获得的所有其他实施例，都属于本发明保护的范围。
如图1所示，本发明提供了一种基于SEM图像的碳酸盐岩储层孔隙半径分布定量方法，该方案能够对小于5微米以下的自形、半自形、它形的碳酸盐岩储层半径进行定量分析，并且能够提供可视化结果，以解决现有技术中高压压汞法无法实现直视化测试，且引起新的孔隙造成测量误差的问题。
具体包括以下步骤：
步骤100、准备各类自形程度的碳酸盐岩储层，获取每种碳酸盐岩储层在电镜扫描下的SEM图像，SEM图像内按照自适应阈值方式对每种碳酸盐岩储层二值化处理并转化为二值图像。
在步骤100中，将SEM图像转化二值图像的具体实现步骤为：
先将SEM图像转化为矩阵图像，并对矩阵图像平均切割为多个图像碎片。
然后抽样采集每个图像碎片中多个像素点的像素值，且确定对每个图像碎片二值化处理对应的分界像素值，其中确定分界像素值的实现方式为：
将抽样采集的列间距和行间距形成抽样网格，并获取抽样网格的每个交叉点对应的像素值；求取所有像素点的平均像素值，且将平均像素值作为分界像素值。将每个图像碎片中低于分界像素值的像素点的像素值设置为0，且将高于分界像素值的像素点的像素值设置为255。
最后将所有的图像碎片按照行列式对齐，形成二值化处理后的二值图像。
图像的二值化就是将图像上的像素点的灰度值设置为0或255，也就是将整个图像呈现出明显的只有黑和白的视觉效果，二值图像也是一种包含目标特征的一种能被计算机所识别的图像，一般将图像二值化基本原理主要是在0-255之间设定一个阈值，大于该阈值的像素，变为255(白色)，小于该阈值的像素变为O(黑色)。
但是当二值图像的亮度不均时，用一个阈值去处理一系列二值图像势必会造成结果的偏离，因此本实施方式采用自适应阈值的设定，即当同一幅图像上的不同部分的具有不同亮度时，确定每个部分的自适应阈值，二值图像上的每一个小区域对应的阈值均不相同，因此在同一幅图像上的不同区域采用的是不同的阈值，从而能在亮度不同的情况下得到更好的结果。
步骤200、将像素值突变的像素点围成的小面积作为非孔隙特征进行消除，计算像素值突变的像素点围成的大面积之间的邻近域包含的像素点以将至少两个大面积连接修复。
在步骤200中，建立关于二值图像的二维坐标系，按照矩阵遍历方式识别二值图像中所有像素点对应的像素值，标注像素值突变的像素点，并根据像素值突变的像素点围成的面积选择消除噪点或修复噪点，消除噪点的实现方法为：
①以二值图像的垂直交叉边缘为坐标轴建立二维坐标系，且按照先行后列的方式遍历二值图像以确定每个像素点的像素值，选定行列遍历时的遍历间距均为一个像素，确定出每行中像素值突变的像素点，且将像素值突变的像素点保存在一个集合内。
按照行列式遍历二值图像每个像素点时，可得到每个像素点的行列坐标以及对应每个像素点的像素值，将每次像素点的像素值与上一个像素点的像素值进行对比，当遍历到像素点的像素值由255突变到0时，则认为该像素点对应孔隙单元的边缘位置，当遍历到像素点的像素值由0突变到255时，则将像素值为0的像素点作为孔隙单元的边缘位置。
②将相邻行的像素值突变的像素点统计在同一个集合内，并且判断每个集合中包含的像素值突变的像素点个数，统计在同一个集合内的像素值突变的像素点为同一个孔隙单元的孔隙边缘对应的像素点，那么如何将像素值突变的像素点分别归类到同一个集合内，具体的实现方式为：
将每个像素值突变的像素点与上一行像素值突变的像素点的行坐标进行对比，如果此时的像素点的行坐标有两个，且与上一行像素值突变的像素点的行坐标的差距为1，那么将此时遍历到的像素点与上一行的像素点保存在同一个集合内。
③将个数少于等于设定值的集合定义为噪点，并重新设定该像素点的像素值为255以消除该噪点，比如说当一个集合仅包含两个或者少于两个的像素点行列坐标，那么将定义该集合对应的像素点为噪点，将该像素点的像素值重新定义为255(白色)以消除噪点
其中修复噪点的实现方式为：
将个数多于设定值的集合定义为待修复噪点，对两个相邻的集合进行多次方数字形态学运算，根据运算结果将两个相邻的集合之间的像素点的像素值重置为0以进行孔隙连接，即计算两个的集合中同行或者同列像素点之间的邻域像素，并在邻域像素小于设定值时，如邻域像素小于8个，则将两个孔隙单元之间的多行或者多列像素值为255的像素点设置为0以将两个像素点连接。
由于扫描电镜方法具有超强的放大倍数，能够观察至纳米级孔隙，在这个过程中受高压影响，SEM图像在拍摄过程中会产生一定数量的噪点，这对于纳米级与微米级孔隙的影响是不可忽视的。因此需要人为对所获取的二值图像进行降噪，降噪包括对非孔隙特征噪点的消除以及受噪点影响导致孔隙被磨灭的空间进行弥补和修复。消除噪点和修复噪点分别采用的MATLAB中函数clean与函数bridge，函数clean是移除孤立的像素，如被0像素点包围的1像素点，也就是将0值像素置255，函数bridge是用于连接断开的像素。如果他有两个非零的不相连(8邻域)的像素。
步骤300、重新遍历完成噪点消除和噪点修复的二值图像，将像素值突变的像素点作为孔隙边界，并且基于孔隙边缘划分二值图像中的孔隙空间个体。
在步骤300中，识别每种碳酸盐岩储层对应二值图像中的孔隙边缘和孔隙空间个体的实现方法为：
步骤301、将每个集合内的像素点的行列坐标转化为树状图，且将同列的像素点设置在树状图同一层的节点上；
步骤302、在二值图像内重新定义对应集合内的像素点的RGB值，确定遍历过程中的每一行的第一个像素值突变的像素点，按照遍历顺序划分每一行的像素值突变的像素点的优先级，并且将不同行的像素值突变的像素点的优先级设定为与遍历顺序一一对应，并且将同一行的像素值突变的像素点设置为同一个优先级；
将仅有一个像素值突变的像素点按照优先级顺序分别设置为一个孔隙空间个体的起点和终点，按照从起点到终点的顺序分为两条测绘曲线，两条测绘曲线顺次将不同优先级的像素值突变的像素点依次连接形成一个孔隙空间个体。
步骤303、将所有重新定义RGB值的像素点首尾连接构成的曲线设定为孔隙边缘，且将所有重新定义RGB值的像素点包围的区域设定为孔隙空间个体。
单个孔隙空间个体即为一个孔隙单元，从而实现对孔隙单元的分割，便于查看识别，并且通过将分割的孔隙单元尺寸与求出的每个孔隙单元的长轴短轴尺寸对比，以实现对碳酸盐岩储层孔隙半径分析的人工验证工作。
步骤400、确定每个像素点在二维坐标系内的坐标值，基于坐标值计算二值图像中的每个孔隙空间个体的孔隙半径，并输出每种碳酸盐岩储层的孔隙半径分布直方图和平均孔隙半径。
孔隙半径为单个孔隙空间个体的长轴与短轴的平均值的一半，根据上述，按照行列遍历二值图像时，将像素值由255突变到0的像素点作为孔隙边缘的像素点，并将此像素点的行列坐标保存在单个孔隙空间个体对应的集合内，因此当计算每个孔隙空间个体的孔隙半径时，先查找出每个集合内的行坐标最大值、行坐标最小值、列坐标最大值和列坐标最小值，然后将行坐标最大值、行坐标最小值、列坐标最大值和列坐标最小值转化为在二维坐标系的最大X坐标，最小X坐标，最大Y坐标，最小Y坐标，求取单个孔隙空间个体的长轴与短轴，具体的实现步骤为：
选择每个孔隙空间个体对应集合内像素值突变的像素点的横坐标最大值对应的像素点、横坐标最小值对应的像素点，纵坐标最大值的像素点和坐标表最小值的像素点，每个孔隙空间个体的长轴半径为Max[(Xmax-Xmin)，(Ymax-Ymin)]，而每个孔隙空间个体的短轴半径为Min[(Xmax-Xmin)，(Ymax-Ymin)]，将通过公式[(Xmax-Xmin)+(Ymax-Ymin)]/2计算每个孔隙空间个体的孔隙半径。
统计每种碳酸盐岩储层中的每个孔隙空间个体对应的孔隙半径，并确定孔隙半径相同的孔隙空间个体的数量，创建关于孔隙半径和孔隙空间个体数量的孔隙半径直方图，并计算每种碳酸盐岩储层样品的孔隙半径平均值。
本实施方式基于常规扫描电镜图像，并结合计算机算法，成分非常低，非常省时，能够对小于5微米以下的自形、半自形、它形的碳酸盐岩储层半径进行定量分析，并且能够提供可视化结果。
基于上述碳酸盐岩储层孔隙半径的定量方法，本实施方式做出以下实验操作，通过对鄂尔多斯盆地奥陶系马家沟组中自形、半自形、它形致密碳酸盐岩储层的孔隙半径进行图像二值化-降噪处理-孔隙单元划分-自动识别-定量计算这个整套流程，以得到定量分析的实验结果。
如图2所示，分别提供自形碳酸盐岩储层的SEM图像，半自形碳酸盐岩储层的SEM图像，以及它形碳酸盐岩储层的SEM图像，对三类自形程度的碳酸盐岩储层扫描电镜图像二值化处理，先通过将每种SEM图像拆分为多尺寸的灰度图像阵列，并且抽样识别灰度图像阵列的不同灰度大小，便可得到自适应赋值前身的多尺度灰度图像的二值分布情况，进一步对每个拆分的灰度图像阵列进行二值化并依次拼接，便可得到如图3所示的目标二值图像。
通过函数clean与函数bridge对已获取的二值图像进行降噪，包括对非孔隙特征的噪点的消除以及受噪点影响导致孔隙被磨灭的空间进行弥补和修复。最终得到了如图4所示的自形、半自形、它形三类碳酸扫描电镜降噪图像。
对降噪处理后不同自形程度的碳酸盐岩储层中的孔隙半径的计算，先进行孔隙单元的识别与划分，再识别每个孔隙单元的半径，最后输出孔隙半径直方图并计算孔隙半径的平均值，最终如图5所示，自形、半自形、它形三类碳酸盐岩储层平均孔隙半径分别为2.2445μm、2.5229μm、1.8560μm，上述工作过程均可以在MATLAB中实现，完成对鄂尔多斯盆地马家沟组自形、半自形、它形碳酸盐岩储层孔隙识别、分布和大小的定量计算，以及可视化输出结果。
以上实施例仅为本申请的示例性实施例，不用于限制本申请，本申请的保护范围由权利要求书限定。本领域技术人员可以在本申请的实质和保护范围内，对本申请做出各种修改或等同替换，这种修改或等同替换也应视为落在本申请的保护范围内。

Claims (10)

1.一种基于SEM图像的碳酸盐岩储层孔隙半径分布定量方法，其特征在于，包括以下步骤：
步骤100、准备各类自形程度的碳酸盐岩储层，获取每种碳酸盐岩储层在电镜扫描下的SEM图像，所述SEM图像内按照自适应阈值方式对每种碳酸盐岩储层二值化处理并转化为二值图像；
步骤200、将像素值突变的像素点围成的小面积作为非孔隙特征进行消除，计算像素值突变的像素点围成的大面积之间的邻近域包含的像素点以将至少两个大面积连接修复；
步骤300、重新遍历完成噪点消除和噪点修复的所述二值图像，将像素值突变的像素点作为孔隙边界，并且基于孔隙边缘划分所述二值图像中的孔隙空间个体；
步骤400、确定每个像素点在所述二维坐标系内的坐标值，基于坐标值计算所述二值图像中的每个所述孔隙空间个体的孔隙半径，并输出每种所述碳酸盐岩储层的孔隙半径分布直方图和平均孔隙半径。
2.根据权利要求1所述的一种基于SEM图像的碳酸盐岩储层孔隙半径分布定量方法，其特征在于：在步骤100中，将所述SEM图像转化二值图像的具体实现步骤为：
将所述SEM图像转化为矩阵图像，并对所述矩阵图像平均切割为多个图像碎片；
抽样采集每个所述图像碎片中多个像素点的像素值，且确定对每个所述图像碎片二值化处理对应的分界像素值；
将每个所述图像碎片中低于分界像素值的像素点的像素值设置为0，且将高于分界像素值的像素点的像素值设置为255；
将所有的图像碎片按照行列式对齐，形成二值化处理后的二值图像。
3.根据权利要求2所述的一种基于SEM图像的碳酸盐岩储层孔隙半径分布定量方法，其特征在于，确定分界像素值的实现方式为：
将抽样采集的列间距和行间距形成抽样网格，并获取所述抽样网格的每个交叉点对应的像素值；
求取所有像素点的平均像素值，且将所述平均像素值作为分界像素值。
4.根据权利要求1所述的一种基于SEM图像的煤储层孔隙结构参数定量分析方法，其特征在于：在步骤200中，建立关于所述二值图像的二维坐标系，按照矩阵遍历方式识别所述二值图像中所有像素点对应的像素值，标注像素值突变的像素点，并根据像素值突变的像素点围成的面积选择消除噪点或修复噪点，消除噪点的实现方法为：
以所述二值图像的垂直交叉边缘为坐标轴建立二维坐标系，且按照先行后列的方式遍历所述二值图像以确定每个所述像素点的像素值，选定行列遍历时的遍历间距均为一个像素，确定出每行中像素值突变的像素点，且将像素值突变的像素点保存在一个集合内；
将相邻行且与所述集合内上一个像素点的行坐标差值为1的像素值突变的像素点统计在同一个集合内，并且判断每个集合中包含的像素值突变的像素点个数；
将个数少于等于设定值的集合定义为噪点，并重新设定该像素点的像素值为255以消除该噪点。
5.根据权利要求4所述的一种基于SEM图像的碳酸盐岩储层孔隙半径分布定量方法，其特征在于，修复噪点的实现方式为：
将个数多于设定值的集合定义为待修复噪点，对两个所述集合进行多次方数字形态学运算，根据运算结果将两个所述集合之间的像素点的像素值重置为0以进行孔隙连接。
6.根据权利要求5所述的一种基于SEM图像的碳酸盐岩储层孔隙半径分布定量方法，其特征在于，对两个相邻的所述集合进行多次方数字形态学运算具体为：计算两个所述集合中同行或者同列像素点之间的邻域像素，并在所述邻域像素小于设定值时，将所述邻域像素内像素值为255的像素点设置为0以将两个所述像素点连接。
7.根据权利要求6所述的一种基于SEM图像的碳酸盐岩储层孔隙半径分布定量方法，其特征在于：在步骤300中，识别每种所述碳酸盐岩储层对应所述二值图像中的孔隙边缘和孔隙空间个体的实现方法为：
步骤301、将每个集合内的所述像素点的行列坐标转化为树状图，且将同列的像素点设置在所述树状图同一层的节点上；
步骤302、在所述二值图像内重新定义对应集合内的所述像素点的RGB值；
步骤303、将所有重新定义RGB值的像素点首尾连接构成的曲线设定为所述孔隙边缘，且将所有重新定义RGB值的像素点包围的区域设定为孔隙空间个体。
8.根据权利要求7所述的一种基于SEM图像的煤储层孔隙结构参数定量分析方法，其特征在于：在步骤302中，确定遍历过程中的每一行的第一个像素值突变的所述像素点，按照遍历顺序划分每一行的像素值突变的所述像素点的优先级，并且将不同行的像素值突变的所述像素点的优先级设定为与遍历顺序一一对应，并且将同一行的像素值突变的所述像素点设置为同一个优先级；
将仅有一个像素值突变的所述像素点按照优先级顺序分别设置为一个所述孔隙空间个体的起点和终点，按照从起点到终点的顺序分为两条测绘曲线，两条所述测绘曲线顺次将不同优先级的所述像素值突变的所述像素点依次连接形成一个所述孔隙空间个体。
9.根据权利要求1所述的一种基于SEM图像的碳酸盐岩储层孔隙半径分布定量方法，其特征在于：在步骤400中，所述孔隙半径为单个孔隙单位的长轴与短轴的平均值的一半，所述二值图像中每个所述孔隙空间个体的所述孔隙半径利用每个所述孔隙空间个体的像素点在所述二维坐标系的最大X坐标，最小X坐标，最大Y坐标，最小Y坐标进行求取，具体的实现步骤为：
选择每个所述孔隙空间个体内像素值突变的所述像素点的横坐标最大值对应的像素点、横坐标最小值对应的像素点，纵坐标最大值的像素点和坐标表最小值的像素点；
每个所述孔隙空间个体的长轴半径为Max(Xmax-Xmin，Ymax-Ymin)，而每个所述孔隙空间个体的短轴半径为Min(Xmax-Xmin，Ymax-Ymin)；
将通过公式[(Xmax-Xmin)+(Ymax-Ymin)]/2计算每个所述孔隙空间个体的孔隙半径。
10.根据权利要求9所述的一种基于SEM图像的碳酸盐岩储层孔隙半径分布定量方法，其特征在于：在步骤400中，统计每种碳酸盐岩储层中的每个孔隙空间个体对应的所述孔隙半径，并确定所述孔隙半径相同的所述孔隙空间个体的数量，创建关于孔隙半径和孔隙空间个体数量的孔隙半径直方图，并计算每种碳酸盐岩储层样品的孔隙半径平均值。
CN202110760736.9A 2021-07-06 2021-07-06 基于sem图像的碳酸盐岩储层孔隙半径分布定量方法 Active CN113570651B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202110760736.9A CN113570651B (zh) 2021-07-06 2021-07-06 基于sem图像的碳酸盐岩储层孔隙半径分布定量方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202110760736.9A CN113570651B (zh) 2021-07-06 2021-07-06 基于sem图像的碳酸盐岩储层孔隙半径分布定量方法

Publications (2)

Publication Number Publication Date
CN113570651A true CN113570651A (zh) 2021-10-29
CN113570651B CN113570651B (zh) 2023-02-07

Family

ID=78163754

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202110760736.9A Active CN113570651B (zh) 2021-07-06 2021-07-06 基于sem图像的碳酸盐岩储层孔隙半径分布定量方法

Country Status (1)

Country Link
CN (1) CN113570651B (zh)

Cited By (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN115410049A (zh) * 2022-10-31 2022-11-29 中国石油大学(华东) 一种岩体溶蚀程度的分类评估方法及装置
CN117409408A (zh) * 2023-12-15 2024-01-16 北京大学 层理缝参数获取方法、装置、设备及可读存储介质

Citations (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN102496020A (zh) * 2011-10-31 2012-06-13 天津大学 基于累积边缘点可视灰度范围直方图的图像二值化方法
CN108763997A (zh) * 2018-04-13 2018-11-06 广州中大微电子有限公司 一种基于交叉方向快速二值化的二维码降噪方法及系统
CN108956416A (zh) * 2018-06-06 2018-12-07 中国地质大学（北京） 一种基于Matlab分析致密砂岩储层孔隙表征的方法

Patent Citations (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN102496020A (zh) * 2011-10-31 2012-06-13 天津大学 基于累积边缘点可视灰度范围直方图的图像二值化方法
CN108763997A (zh) * 2018-04-13 2018-11-06 广州中大微电子有限公司 一种基于交叉方向快速二值化的二维码降噪方法及系统
CN108956416A (zh) * 2018-06-06 2018-12-07 中国地质大学（北京） 一种基于Matlab分析致密砂岩储层孔隙表征的方法

Non-Patent Citations (3)

* Cited by examiner, † Cited by third party
Title
BO JIU 等: "《Quantitative Analysis of Micron-Scale and Nano-Scale Pore Throat Characteristics of Tight Sandstone Using Matlab》", 《APPLIED SCIENCES》 *
柳林: "《基于OpenCV的数字图像处理技术》", 31 December 2020 *
郭达志: "《空间信息技术与资源环境保护》", 30 November 2007 *

Cited By (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN115410049A (zh) * 2022-10-31 2022-11-29 中国石油大学(华东) 一种岩体溶蚀程度的分类评估方法及装置
CN117409408A (zh) * 2023-12-15 2024-01-16 北京大学 层理缝参数获取方法、装置、设备及可读存储介质
CN117409408B (zh) * 2023-12-15 2024-03-08 北京大学 层理缝参数获取方法、装置、设备及可读存储介质

Also Published As

Publication number Publication date
CN113570651B (zh) 2023-02-07

Similar Documents

Publication Publication Date Title
CN113570651B (zh) 基于sem图像的碳酸盐岩储层孔隙半径分布定量方法
CN108053417B (zh) 一种基于混合粗分割特征的3D U-Net网络的肺分割装置
CN107507173B (zh) 一种全切片图像的无参考清晰度评估方法及系统
CN109978839B (zh) 晶圆低纹理缺陷的检测方法
CN108447050A (zh) 一种基于超像素的工件表面缺陷分割方法
CN108805863B (zh) 深度卷积神经网络结合形态学检测图像变化的方法
CN117078671B (zh) 一种甲状腺超声影像智能分析系统
CN112184725B (zh) 一种沥青路面图像的结构光光条中心提取方法
JP2013238449A (ja) ひび割れ検出方法
CN115272319B (zh) 一种矿石粒度检测方法
CN113570652B (zh) 基于sem图像的砂岩储层矿物晶间孔的定量分析方法
CN114820625A (zh) 一种汽车顶块缺陷检测方法
CN114881965A (zh) 一种基于人工智能和图像处理的木板节子检测方法
CN113222992A (zh) 基于多重分形谱的裂纹特征表征方法及系统
CN114494318B (zh) 基于大津算法的角膜动态形变视频提取角膜轮廓的方法
CN116523898A (zh) 一种基于三维点云的烟草表型性状提取方法
CN110930423B (zh) 一种物体边缘特征识别提取方法
CN112784922A (zh) 智能云端医疗影像的提取和分类方法
CN113591740B (zh) 基于深度学习的复杂河流环境下泥沙颗粒识别方法及装置
CN115690073A (zh) 一种激光增材制造陶瓷微结构局部表征方法、装置及介质
CN114152211B (zh) 一种基于显微图像处理的压裂支撑剂圆度测量方法
CN112862816B (zh) 一种hrtem图像中煤芳香烃晶格条纹的智能提取方法
CN104504702A (zh) 基于方格搜索法的水泥刻槽路面裂缝识别方法
CN111709939B (zh) 一种结构对称性的编织陶瓷基复合材料细观组分分类方法
CN108198245B (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