CN110021030B - 一种岩土体材料数字图像的分割阈值确定方法 - Google Patents

一种岩土体材料数字图像的分割阈值确定方法 Download PDF

Info

Publication number
CN110021030B
CN110021030B CN201910174541.9A CN201910174541A CN110021030B CN 110021030 B CN110021030 B CN 110021030B CN 201910174541 A CN201910174541 A CN 201910174541A CN 110021030 B CN110021030 B CN 110021030B
Authority
CN
China
Prior art keywords
gray level
gray
rock
segmentation threshold
derivative
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.)
Active
Application number
CN201910174541.9A
Other languages
English (en)
Other versions
CN110021030A (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 Mining and Technology CUMT
Original Assignee
China University of Mining and Technology CUMT
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 Mining and Technology CUMT filed Critical China University of Mining and Technology CUMT
Priority to CN201910174541.9A priority Critical patent/CN110021030B/zh
Priority to US16/965,307 priority patent/US20210201497A1/en
Priority to PCT/CN2019/086147 priority patent/WO2020177215A1/zh
Publication of CN110021030A publication Critical patent/CN110021030A/zh
Application granted granted Critical
Publication of CN110021030B publication Critical patent/CN110021030B/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/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/11—Region-based segmentation
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T7/00—Image analysis
    • G06T7/40—Analysis of texture
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06V—IMAGE OR VIDEO RECOGNITION OR UNDERSTANDING
    • G06V10/00—Arrangements for image or video recognition or understanding
    • G06V10/40—Extraction of image or video features
    • G06V10/44—Local feature extraction by analysis of parts of the pattern, e.g. by detecting edges, contours, loops, corners, strokes or intersections; Connectivity analysis, e.g. of connected components
    • 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
    • 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/30—Subject of image; Context of image processing
    • G06T2207/30181—Earth observation
    • 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
    • Y02P—CLIMATE CHANGE MITIGATION TECHNOLOGIES IN THE PRODUCTION OR PROCESSING OF GOODS
    • Y02P90/00—Enabling technologies with a potential contribution to greenhouse gas [GHG] emissions mitigation
    • Y02P90/30—Computing systems specially adapted for manufacturing

Landscapes

  • Engineering & Computer Science (AREA)
  • Computer Vision & Pattern Recognition (AREA)
  • Physics & Mathematics (AREA)
  • General Physics & Mathematics (AREA)
  • Theoretical Computer Science (AREA)
  • Multimedia (AREA)
  • Image Analysis (AREA)
  • Analysing Materials By The Use Of Radiation (AREA)

Abstract

本发明公开了一种岩土体材料数字图像的分割阈值确定方法，所述方法包括如下步骤：S1：获取岩土体材料的SEM图像的灰度直方图曲线；S2：根据所述灰度直方图曲线，确定分割阈值T的取值范围；S3：获取所述灰度直方图曲线的二阶导数；S4：根据所述灰度直方图曲线的二阶导数和分割阈值T的取值范围，确定所述分割阈值T的大小。本发明能够快速准确的确定岩土体材料的数字图像的分割阈值，准确地从数字图像中将岩土体材料孔隙或裂缝结构与表面土骨架结构区分开来，为后续基于数字图像的岩土体材料的深入研究提供了准确的分割阈值，为准确的提取孔隙或裂缝结构，提供了有效的技术支持。

Description

一种岩土体材料数字图像的分割阈值确定方法
技术领域
本发明涉及图像分割技术领域，尤其涉及一种岩土体材料数字图像的分割阈值确定方法。
背景技术
随着微、细观成像技术的发展，研究人员可以利用该技术对岩土体的孔隙结构形态特征进行直接观测。对于SEM扫描电镜技术，由于其较低的测试成本和良好的成像效果，被逐渐应用于各个领域。
使用SEM扫描电镜技术对各种岩土体材料的表面微观孔隙结构进行直接观测，通过定量化的分析，对岩土体材料渗透率进行预测。在对SEM图像定量化分析的过程中，最关键的一步就是对图像进行二值化。二值化的主要目的是将SEM图像中的岩土体材料的孔隙结构与表面土骨架结构区分开来，提取出其中的孔隙结构，而孔隙结构是决定岩土体材料渗透率大小的主要因素，因此如何准确的提取孔隙结构就是数字图像法中二值化过程最为关键的技术。而二值化处理过程中，图像的分割阈值确定决定着孔隙结构提取二值化的准确性，因此，确定图像的分割阈值确定显得极为重要。
近几十年来，国内外众多学者已经提出了大量的图像二值化分割阈值确定计算方法。但由于数字图像的复杂性，以及图像分割问题依赖于具体应用领域的具体情况，因此，直到现在为止没有一种通用的二值化分割阈值确定算法。
中国专利公布号：CN 102841220 A；公布日：2012年12月26日，公开了一种基于孔隙率的黏土扫描电镜照片图像分割方法，包括步骤有：采用土力学常规实验方法测得黏土试样的干密度，算出干密度的三分之二次幂；实验室制备黏土试样的微结构样品；采用ZD-A3冻干仪将制作的样品冻干并抽真空，真空度达到到10-6Pa时进行扫描电镜，观察样品的微结构并获得扫描照片；采用图像分析软件提取两个阈值对应的二值化数据；对比两个阈值条件下抽取到颗粒体和孔隙率的参数特征，算出孔隙体面积比；当孔隙体面积比等于干密度的三分之二次幂时，停止逼近，此时的阈值即为与所测干密度相对应的图像分割阈值确定。有益效果是采用该方法实现了三维参数向二维化的科学转变，排除了人为主观因数的干扰，使得土体微结构的研究更加准确。但是该发明的不足之处在于：该发明中Leica QWin图像分析软件对SEM照片提取出的阈值，虽然进行了人为确认，但是在进行人为确认的过程中，涉及多次运算，过程繁琐的同时还易出错，从而不利于长期使用。
发明内容
发明目的：针对现有岩土体材料数字图像二值化中分割阈值的确定过程复杂易出错的问题，本发明提出一种岩土体材料数字图像的分割阈值确定方法。
技术方案：为实现本发明的目的，本发明所采用的技术方案是：
一种岩土体材料数字图像的分割阈值确定方法，所述方法包括如下步骤：
S1：获取岩土体材料的SEM图像的灰度直方图曲线；
S2：根据所述灰度直方图曲线，确定分割阈值T的取值范围；
S3：获取所述灰度直方图曲线的二阶导数；
S4：根据所述灰度直方图曲线的二阶导数和分割阈值T的取值范围，确定所述分割阈值T的大小。
更进一步地，所述步骤S1获取图像的灰度直方图曲线之前还包括：读取所述岩土体材料的SEM图像，获取所述SEM图像中每个像素的灰度级i、所述每个像素的灰度级i所对应的总像素数量ni。
更进一步地，所述步骤S1获取图像的灰度直方图曲线具体如下：
S1.1：确定所述岩土体材料灰度图像中每个像素的灰度级i；
S1.2：根据如下公式，获取所述岩土体材料灰度图像的灰度直方图曲线上各个点：
Figure BDA0001989096900000021
其中，i表示灰度级，N表示图像像素总数，ni表示图像中所有灰度级为i的像素的总像素个数，L表示灰度级的种类数；
S1.3：根据所述P(i)，拟合获取岩土体材料灰度图像的灰度直方图曲线。
更进一步地，所述步骤S2确定分割阈值T的取值范围具体如下：
步骤S2.1：根据所述灰度直方图曲线，确定所述灰度直方图曲线中的峰值数目；
步骤S2.2：根据所述峰值数目，确定所述岩土体材料的结构；
步骤S2.3：获取所述峰值对应像素的灰度级imax；
步骤S2.4：根据所述岩土体材料的结构和峰值对应像素的灰度级imax，确定所述分割阈值T的取值范围。
更进一步地，步骤S2.2确定所述岩土体材料的结构具体如下：
如果所述灰度直方图曲线只有一个峰值，所述岩土体材料的结构为含孔隙结构；
如果所述灰度直方图曲线有两个峰值，所述岩土体材料的结构为含裂缝结构。
更进一步地，步骤S3获取所述灰度直方图曲线的二阶导数之前包括：获取所述灰度直方图曲线的一阶导数，所述一阶导数为:
Figure BDA0001989096900000022
其中，i表示灰度级，ni表示图像中每个像素的灰度级i所对应的总像素数量。
更进一步地，确定所述孔隙结构的分割阈值T具体如下：
SA4.1：所述灰度直方图曲线只有一个峰值，所述分割阈值T的取值范围为：
T＜imax
其中，imax为峰值对应像素的灰度级；
SA4.2：获取所述灰度直方图曲线的二阶导数，所述二阶导数为：
Figure BDA0001989096900000031
其中，i表示灰度级，ni表示图像中每个像素的灰度级i所对应的总像素数量；
SA4.3：根据所述二阶导数，确定在分割阈值T的取值范围内二阶导数的最大值aimax；
SA4.4：确定所述二阶导数最大值aimax所对应像素的灰度级iT,所述分割阈值T为：
T＝iT
其中，iT为二阶导数最大值aimax所对应像素的灰度级。
更进一步地，确定所述裂缝结构的分割阈值T具体如下：
SB4.1：所述灰度直方图曲线有两个峰值，所述分割阈值T的取值范围为：
imax1＜T＜imax2
其中，imax1为第一个峰值对应像素的灰度级，imax2为第二个峰值对应像素的灰度级；
SB4.2：获取所述灰度直方图曲线的二阶导数，所述二阶导数为：
Figure BDA0001989096900000032
其中，ni表示图像中所有灰度级为i的像素的总像素个数，i表示灰度级；
SB4.3：根据所述二阶导数，确定在分割阈值T的取值范围内二阶导数的最大值aimax；
SB4.4：确定所述二阶导数最大值aimax所对应像素的灰度级iT,所述分割阈值T为：
T＝iT
其中，iT为二阶导数最大值aimax所对应像素的灰度级。
有益效果：与现有技术相比，本发明的技术方案具有以下有益技术效果：
(1)本发明先是通过待测岩土体材料的SEM图像的灰度直方图曲线，确定分割阈值的取值范围，再通过灰度直方图曲线的二阶导数确定分割阈值的具体大小，通过不断的缩小范围直至精确取值，能够更进一步地保证分割阈值取值的准确性；
(2)本发明根据待测岩土体材料自身的SEM图像进行分析，能够确保分割阈值更加满足待测岩土体材料数字图像二值化的需求，从而准确的从数字图像中将岩土体材料的孔隙或裂缝结构与表面土骨架结构区分开来；
(3)本发明为后续基于数字图像的岩土体材料的深入研究提供了准确的分割阈值，为准确的提取岩土体材料的孔隙或裂缝结构，提供了有效的技术支持。
附图说明
图1是本发明的流程示意图；
图2是膨润土SEM图像；
图3是膨润土颗粒基本单元的示意图；
图4是SEM扫描岩土体材料原理剖面示意图；
图5是灰度直方图曲线对应岩土体结构示意图；
图6是不同分割阈值压实膨润土二值图；
图7是SEM图像随分割阈值增大的二值化提取过程；
图8是压裂后煤样SEM图像；
图9是压裂后煤样灰度直方图曲线；
图10是膨润土的灰度直方图曲线；
图11是膨润土的灰度直方图曲线二阶导数曲线；
图12是膨润土的最终二值化图像；
图13是压裂后的煤样的灰度直方图曲线二阶导数曲线；
图14是压裂后的煤样的最终二值化图像。
具体实施方式
为使本发明实施例的目的、技术方案和优点更加清楚，下面将结合本发明实施例中的附图，对本发明实施例中的技术方案进行清楚、完整地描述。其中，所描述的实施例是本发明一部分实施例，而不是全部的实施例。因此，以下对在附图中提供的本发明的实施例的详细描述并非旨在限制要求保护的本发明的范围，而是仅仅表示本发明的选定实施例。
参考图1，本实施例提供了一种岩土体材料数字图像的分割阈值确定方法，本方法计算得到的分割阈值能够有效地提取出SEM图像中的孔隙结构或裂缝结构，同时也符合岩土体材料孔隙及裂缝的实际分布情况。从而为基于数字图像法的岩土体材料细观机理研究提供了良好的技术支持。
参考图4，图4为扫描电子显微镜通过用聚焦电子束扫描样品的表面来产生样品表面的图像，其中扫描电子显微镜英文为：Scanning Electron Microscope，可以缩写为SEM，简称扫描电镜，是一种电子显微镜。
使用SEM电镜扫描得到的岩土体材料表面微观结构图像为灰度图像，其中灰度图像的灰度级从0到255，总共有256种取值，同时SEM图像中每个像素点颜色的深浅都代表着一个灰度级。
由于岩土体材料的孔隙结构或裂缝结构距离电子束发射口较远，因此最终的SEM图像中代表孔隙或裂隙结构的像素点呈现出的灰度级一般在0～90；而岩土体材料的表面土骨架结构上的凸起颗粒距离较近，因此最终图像中代表岩土体材料表面土骨架结构的像素点呈现出的灰度级一般在150～255；同时其中岩土体材料的表面土骨架结构占大多数，具体大于50％，且处于同一平面中，因此其灰度级一般在90～150，对应像素数量最多。
在本实施例中，具体地讲，本分割阈值确定方法具体包括如下步骤：
步骤S1：通过MATALB代码读取待测岩土体材料的SEM图像，获取SEM图像中每个像素的灰度级i以及每个像素的灰度级i所对应的总像素数量ni。
步骤S2：获取待测岩土体材料的SEM图像的灰度直方图曲线，具体过程如下：
步骤S2.1：根据每个像素的灰度级i、每个像素的灰度级i所对应的总像素数量ni以及如下公式，获取待测岩土体材料灰度图像的灰度直方图曲线上的各个点：
Figure BDA0001989096900000051
其中，i表示灰度级，N表示图像像素总数，ni表示图像中所有灰度级为i的像素的总像素个数，L表示灰度级的种类数；
步骤S2.2：根据步骤S2.1中的各个P(i)，拟合获取岩土体材料灰度图像的灰度直方图曲线。
步骤S3：根据灰度直方图曲线，确定分割阈值T的取值范围，具体过程如下：
步骤S3.1：根据灰度直方图曲线，确定灰度直方图曲线中的峰值数目；
步骤S3.2：根据峰值数目，确定待测岩土体材料的结构，具体为：
如果灰度直方图曲线中只有一个峰值，待测岩土体材料的结构为含孔隙结构；
如果灰度直方图曲线中有两个峰值，待测岩土体材料的结构为含裂缝结构；
步骤S3.3：根据峰值数目，获取相应峰值对应像素的灰度级imax；
步骤S3.4：根据待测岩土体材料的结构、相应峰值对应像素的灰度级imax，确定分割阈值T的取值范围。
值得注意的是，由于待测岩土体材料结构的不确定，分割阈值T的取值范围的确定也是不固定的，由结构对应峰值的数目决定。
步骤S4：获取灰度直方图曲线的二阶导数，具体过程如下：
步骤S4.1：获取灰度直方图曲线的一阶导数，一阶导数为:
Figure BDA0001989096900000052
其中，i表示灰度级，ni表示图像中每个像素的灰度级i所对应的总像素数量；
步骤S4.2：获取灰度直方图曲线的二阶导数，二阶导数为：
Figure BDA0001989096900000053
其中，i表示灰度级，ni表示图像中每个像素的灰度级i所对应的总像素数量。
步骤S5：根据灰度直方图曲线的二阶导数和分割阈值T的取值范围，确定分割阈值T的大小，具体过程如下：
步骤S5.1：根据步骤S4.2中的各个ai，拟合获取二阶导数曲线；
步骤S5.2：确定在二阶导数曲线中，分割阈值T取值范围内的最大值aimax，同时明确最大值aimax所对应像素的灰度级iT,分割阈值T为：
T＝iT
其中，iT为二阶导数最大值所对应像素的灰度级。
实施例1含孔隙结构的岩土体材料
本实施例提供了一种岩土体材料数字图像的分割阈值确定方法，在本实施例中，具体地讲，待测岩土体材料选择膨润土。其中膨润土数字图像的分割阈值确定方法，具体如下：
步骤SA1：参考图2，图2为膨润土的SEM图像，通过MATALB代码读取膨润土SEM图像，获取膨润土SEM图像中每个像素的灰度级i以及每个像素的灰度级i所对应的总像素数量ni。
步骤SA2：获取膨润土SEM图像的灰度直方图曲线，具体过程如下：
步骤SA2.1：根据每个像素的灰度级i、每个像素的灰度级i所对应的总像素数量ni以及如下公式，获取膨润土灰度图像的灰度直方图曲线上的各个点：
Figure BDA0001989096900000061
其中，i表示灰度级，N表示图像像素总数，ni表示图像中所有灰度级为i的像素的总像素个数，L表示灰度级的种类数。
步骤SA2.2：根据步骤SA2.1中的各个P(i)，拟合获取膨润土灰度图像的灰度直方图曲线。参考图10，图10为膨润土的灰度直方图曲线，从膨润土的灰度直方图曲线中，能够直观的看到膨润土SEM图像中每个灰度级i的像素数量的占比情况。
步骤SA3：根据灰度直方图曲线，确定分割阈值T的取值范围，具体过程如下：
步骤SA3.1：根据灰度直方图曲线，确定灰度直方图曲线中的峰值数目，在本实施例中，参考图10可知，灰度直方图曲线中的峰值只有一个，因此膨润土的结构为孔隙结构。
根据压汞法及岩心覆压孔渗自动测试系统测试得到岩土体的孔隙率小于30％，其与表面土骨架相比来说，孔隙结构占小部分。在本实施例中，具体地讲，根据压汞法及岩心覆压孔渗自动测试系统测试得到干燥压实膨润土的孔隙率为20％～30％，对于膨润土SEM图像来说，其中的孔隙结构应该在25％左右。参考图3可知，膨润土的孔隙结构又分有几个不同的层次，通过SEM电子镜扫描观察到的膨润土孔隙结构属于表观孔隙结构，故通过数字图像法得到的孔隙率应该小于20％。
步骤SA3.2：根据步骤SA3.1可知，膨润土的灰度直方图曲线中只有一个峰值，因此在本实施例中，只需要确定一个峰值对应像素的灰度级imax。
步骤SA3.3：根据膨润土的结构、相应峰值对应像素的灰度级imax，确定分割阈值T的取值范围，因为膨润土的结构为孔隙结构，且其灰度直方图曲线只有一个峰值，故膨润土数字图像的分割阈值T的取值范围为：
T＜imax
其中，imax为峰值对应像素的灰度级。
参考图5，由于所有含孔隙结构的岩土体材料的灰度直方图都是单峰曲线，所以膨润土的结构为孔隙结构。其中由于表面土骨架占大多数，因此峰值处对应像素一定属于表面土骨架。值得注意的是，通过压汞法及岩心覆压孔渗自动测试系统测试得到干燥压实膨润土的孔隙率小于30％，因此，膨润土的数字图像中代表孔隙结构的像素也肯定小于30％，故膨润土中的表面土骨架结构超过70％，因此灰度直方图中峰值部分对应像素肯定是代表表面土骨架。显然地，孔隙结构对应的灰度级一般在0～90，而剩下的在90～225的灰度级对应的则为表面土骨架结构，故代表孔隙结构像素的灰度级小于代表表面土骨架像素的灰度级。
参考图6，图6为随着分割阈值T从0逐渐增大到255时，压实膨润土二值图的变化过程。其中，图6中黑色部分就是根据不同的分割阈值提取出的孔隙结构，显然有一些是非常不合理的。但是通过这个变化过程，可以观察到随着分割阈值T的逐渐增大，孔隙率也在逐渐增大。
参考图7，图7为二值化图随着分割阈值T的变化过程，其中二值化图像孔隙率的增大实际上就是图中黑色像素数量随着分割阈值T的增大而增大。当分割阈值T从孔隙结构的灰度级到达表面土骨架灰度级时，由于表面土骨架结构占大多数，故图像中每个像素的灰度级i所对应的总像素数量ni会存在一个突变，因此，找到这个从孔隙结构到表面土骨架结构突变的这个点对应的灰度级i，即确定分割阈值T的大小。
步骤SA4：获取灰度直方图曲线的二阶导数，具体过程如下：
步骤SA4.1：获取灰度直方图曲线的一阶导数，一阶导数为:
Figure BDA0001989096900000071
其中，i表示灰度级，ni表示图像中每个像素的灰度级i所对应的总像素数量。
步骤SA4.2：获取灰度直方图曲线的二阶导数，二阶导数为：
Figure BDA0001989096900000072
其中，i表示灰度级，ni表示图像中每个像素的灰度级i所对应的总像素数量。
步骤SA5：根据灰度直方图曲线的二阶导数ai和分割阈值T的取值范围，确定分割阈值T的大小，具体过程如下：
步骤SA5.1：根据步骤SA4.2中的各个ai，拟合获取二阶导数曲线，参考图11，图11为膨润土的灰度直方图曲线的二阶导数曲线。
步骤SA5.2：确定在分割阈值T取值范围内，二阶导数曲线中的最大值aimax，同时确定最大值aimax所对应像素的灰度级iT,分割阈值T为：
T＝iT
其中，iT为二阶导数最大值所对应像素的灰度级。
参考图12，图12为膨润土最终的二值化图像。其中分割阈值T＝iT，故在膨润土的灰度直方图曲线中，灰度级i小于分割阈值T处所对应的范围为膨润土的孔隙结构，灰度级i大于分割阈值T处所对应的范围为膨润土的表面土骨架结构。
实施例2含裂缝结构的岩土体材料
本实施例提供了一种岩土体材料数字图像的分割阈值确定方法，在本实施例中，具体地讲，待测岩土体材料选择压裂后的煤样。其中压裂后煤样的数字图像的分割阈值确定方法，具体如下：
步骤SB1：参考图8，图8为压裂后煤样的SEM图像，通过MATALB代码读取压裂后煤样的SEM图像，获取压裂后煤样SEM图像中每个像素的灰度级i以及每个像素的灰度级i所对应的总像素数量ni。
步骤SB2：获取压裂后煤样SEM图像的灰度直方图曲线，具体过程如下：
步骤SB2.1：根据每个像素的灰度级i、每个像素的灰度级i所对应的总像素数量ni以及如下公式，获取压裂后煤样灰度图像的灰度直方图曲线上的各个点：
Figure BDA0001989096900000081
其中，i表示灰度级，N表示图像像素总数，ni表示图像中所有灰度级为i的像素的总像素个数，L表示灰度级的种类数。
步骤SB2.2：根据步骤SB2.1中的各个P(i)，拟合获取压裂后煤样灰度图像的灰度直方图曲线。参考图9，图9为压裂后煤样SEM图像的灰度直方图曲线。
步骤SB3：根据灰度直方图曲线，确定分割阈值T的取值范围，具体过程如下：
步骤SB3.1：根据灰度直方图曲线，确定灰度直方图曲线中的峰值数目，在本实施例中，参考图9可知，灰度直方图曲线中的峰值有两个，因此压裂后煤样的结构为裂缝结构。
步骤SB3.2：根据步骤SB3.1可知，压裂后煤样的灰度直方图曲线中有两个峰值，因此在本实施例中，需要确定两个峰值分别对应像素的灰度级imax1和imax2，其中，imax1为第一个峰值对应像素的灰度级，imax2为第二个峰值对应像素的灰度级。
步骤SB3.3：根据压裂后煤样的结构、相应峰值对应像素的灰度级imax1和imax2，确定分割阈值T的取值范围，因为压裂后煤样的结构为孔隙结构，且其灰度直方图曲线有两个峰值，故压裂后煤样数字图像的分割阈值T的取值范围为：
imax1＜T＜imax2
其中，imax1为第一个峰值对应像素的灰度级，imax2为第二个峰值对应像素的灰度级。
参考图9，含裂缝结构的岩土体材料的灰度直方图曲线均为双峰曲线，故压裂后煤样的结构为含裂缝结构的岩土体材料，在灰度级i从小到大变化的过程中，图像中第一个峰值处的像素代表着裂缝处的像素，而第二个峰值处的像素代表着表面骨架处的像素，因此含裂缝结构的岩土体材料的分割阈值T介于两个峰值之间。
步骤SB4：获取灰度直方图曲线的二阶导数，具体过程如下：
步骤SB4.1：获取灰度直方图曲线的一阶导数，一阶导数为:
Figure BDA0001989096900000091
其中，i表示灰度级，ni表示图像中每个像素的灰度级i所对应的总像素数量。
步骤SB4.2：获取灰度直方图曲线的二阶导数，二阶导数为：
Figure BDA0001989096900000092
其中，i表示灰度级，ni表示图像中每个像素的灰度级i所对应的总像素数量。
步骤SB5：根据灰度直方图曲线的二阶导数ai和分割阈值T的取值范围，确定分割阈值T的大小，具体过程如下：
步骤SB5.1：根据步骤SB4.2中的各个ai，拟合获取二阶导数曲线，参考图13，图13为压裂后煤样的灰度直方图曲线的二阶导数曲线。
步骤SB5.2：确定在分割阈值T取值范围内，二阶导数曲线中的最大值aimax，同时确定最大值aimax所对应像素的灰度级iT,分割阈值T为：
T＝iT
其中，iT为二阶导数最大值所对应像素的灰度级。
参考图14，图14为压裂后煤样最终的二值化图像。其中分割阈值T＝iT，故在压裂后煤样的灰度直方图曲线中，灰度级i小于分割阈值T处所对应的范围为压裂后煤样的裂缝结构，而灰度级i大于分割阈值T处所对应的范围则为压裂后煤样的表面土骨架结构。
以上示意性的对本发明及其实施方式进行了描述，该描述没有限制性，附图中所示的也只是本发明的实施方式之一，实际的结构和方法并不局限于此。所以，如果本领域的普通技术人员受其启示，在不脱离本发明创造宗旨的情况下，不经创造性的设计出与该技术方案相似的结构方式及实施例，均属于本发明的保护范围。

Claims (5)

1.一种岩土体材料数字图像的分割阈值确定方法，其特征在于，所述方法包括如下步骤：
S1：获取岩土体材料的SEM图像的灰度直方图曲线；
S2：根据所述灰度直方图曲线，确定分割阈值T的取值范围，具体如下：
步骤S2.1：根据所述灰度直方图曲线，确定所述灰度直方图曲线中的峰值数目；
步骤S2.2：根据所述峰值数目，确定所述岩土体材料的结构；
步骤S2.3：获取所述峰值对应像素的灰度级imax；
步骤S2.4：根据所述岩土体材料的结构和峰值对应像素的灰度级imax，确定所述分割阈值T的取值范围；
S3：获取所述灰度直方图曲线的二阶导数；
S4：根据所述灰度直方图曲线的二阶导数和分割阈值T的取值范围，确定所述分割阈值T的大小；其中
(1)当所述灰度直方图曲线只有一个峰值，步骤S4具体如下：
SA4.1：所述分割阈值T的取值范围为：
T<imax
其中，imax为峰值对应像素的灰度级；
SA4.2：获取所述灰度直方图曲线的二阶导数，所述二阶导数为：
其中，i表示灰度级，ni表示图像中每个像素的灰度级i所对应的总像素数量；
SA4.3：根据所述二阶导数，确定在分割阈值T的取值范围内二阶导数的最大值aimax；
SA4.4：确定所述二阶导数最大值aimax所对应像素的灰度级iT,所述分割阈值T为：
T＝iT
其中，iT为二阶导数最大值aimax所对应像素的灰度级；
(2)当所述灰度直方图曲线有两个峰值，步骤S4具体如下：
SB4.1：所述分割阈值T的取值范围为：
imax1<T<imax2
其中，imax1为第一个峰值对应像素的灰度级，imax2为第二个峰值对应像素的灰度级；
SB4.2：获取所述灰度直方图曲线的二阶导数，所述二阶导数为：
其中，ni表示图像中所有灰度级为i的像素的总像素个数，i表示灰度级；
SB4.3：根据所述二阶导数，确定在分割阈值T的取值范围内二阶导数的最大值aimax；
SB4.4：确定所述二阶导数最大值aimax所对应像素的灰度级iT,所述分割阈值T为：
T＝iT
其中，iT为二阶导数最大值aimax所对应像素的灰度级。
2.根据权利要求1所述的一种岩土体材料数字图像的分割阈值确定方法，其特征在于，所述步骤S1获取图像的灰度直方图曲线之前还包括：读取所述岩土体材料的SEM图像，获取所述SEM图像中每个像素的灰度级i、所述每个像素的灰度级i所对应的总像素数量ni。
3.根据权利要求1或2所述一种岩土体材料数字图像的分割阈值确定方法，其特征在于，所述步骤S1获取图像的灰度直方图曲线具体如下：
S1.1：确定所述岩土体材料灰度图像中每个像素的灰度级i；
S1.2：根据如下公式，获取所述岩土体材料灰度图像的灰度直方图曲线上各个点：
其中，i表示灰度级，N表示图像像素总数，ni表示图像中所有灰度级为i的像素的总像素个数，L表示灰度级的种类数；
S1.3：根据所述P(i)，拟合获取岩土体材料灰度图像的灰度直方图曲线。
4.根据权利要求1所述的一种岩土体材料数字图像的分割阈值确定方法，其特征在于，步骤S2.2确定所述岩土体材料的结构具体如下：
如果所述灰度直方图曲线只有一个峰值，所述岩土体材料的结构为含孔隙结构；
如果所述灰度直方图曲线有两个峰值，所述岩土体材料的结构为含裂缝结构。
5.根据权利要求4所述的一种岩土体材料数字图像的分割阈值确定方法，其特征在于，步骤S3获取所述灰度直方图曲线的二阶导数之前包括：获取所述灰度直方图曲线的一阶导数，所述一阶导数为:
其中，i表示灰度级，ni表示图像中每个像素的灰度级i所对应的总像素数量。
CN201910174541.9A 2019-03-05 2019-03-05 一种岩土体材料数字图像的分割阈值确定方法 Active CN110021030B (zh)

Priority Applications (3)

Application Number Priority Date Filing Date Title
CN201910174541.9A CN110021030B (zh) 2019-03-05 2019-03-05 一种岩土体材料数字图像的分割阈值确定方法
US16/965,307 US20210201497A1 (en) 2019-03-05 2019-05-09 Method for determining segmentation threshold of digital image of rock-soil material
PCT/CN2019/086147 WO2020177215A1 (zh) 2019-03-05 2019-05-09 一种岩土体材料数字图像的分割阈值确定方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN201910174541.9A CN110021030B (zh) 2019-03-05 2019-03-05 一种岩土体材料数字图像的分割阈值确定方法

Publications (2)

Publication Number Publication Date
CN110021030A CN110021030A (zh) 2019-07-16
CN110021030B true CN110021030B (zh) 2023-04-25

Family

ID=67189412

Family Applications (1)

Application Number Title Priority Date Filing Date
CN201910174541.9A Active CN110021030B (zh) 2019-03-05 2019-03-05 一种岩土体材料数字图像的分割阈值确定方法

Country Status (3)

Country Link
US (1) US20210201497A1 (zh)
CN (1) CN110021030B (zh)
WO (1) WO2020177215A1 (zh)

Families Citing this family (12)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN110458816B (zh) * 2019-08-06 2021-04-27 北京工商大学 一种基于阈值回归的纤维材料孔隙率分析方法
CN112085693B (zh) * 2020-06-24 2022-09-20 中国科学院武汉岩土力学研究所 土石混合体内部结构的孔隙比评估及形态重建方法及系统
CN113112453B (zh) * 2021-03-22 2022-03-22 深圳市华启生物科技有限公司 胶体金检测卡识别方法、系统、电子设备及存储介质
CN113702259A (zh) * 2021-08-19 2021-11-26 国家烟草质量监督检验中心 一种卷烟整体孔隙均匀度的检测方法
CN113702258B (zh) * 2021-08-19 2024-01-19 国家烟草质量监督检验中心 一种卷烟轴向孔隙分布的检测方法
CN114168884B (zh) * 2021-12-02 2024-08-02 成都工业学院 一种岩土样本的概率分布获取方法
CN114705606B (zh) * 2022-04-06 2024-07-02 同济大学 一种基于网络化分析的岩石内部关键渗流节点的封堵方法
CN115049566B (zh) * 2022-08-15 2022-10-25 聊城扬帆田一机械有限公司 一种平板夯激振模式智能调节系统
CN116433663B (zh) * 2023-06-13 2023-08-18 肥城恒丰塑业有限公司 一种土工格室质量智能检测方法
CN116503394B (zh) * 2023-06-26 2023-09-08 济南奥盛包装科技有限公司 基于图像的印刷制品表面粗糙度检测方法
CN116977230B (zh) * 2023-09-22 2024-01-02 济宁市质量计量检验检测研究院(济宁半导体及显示产品质量监督检验中心、济宁市纤维质量监测中心) 一种扫描电子显微镜图像优化增强方法
CN117291945B (zh) * 2023-11-24 2024-02-13 山东省济宁生态环境监测中心(山东省南四湖东平湖流域生态环境监测中心) 基于图像数据的土壤腐蚀污染检测预警方法

Citations (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN101710425A (zh) * 2009-12-25 2010-05-19 南京航空航天大学 基于图像灰度梯度和灰度统计直方图的自适应预分割方法

Family Cites Families (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN105654501B (zh) * 2016-02-22 2019-07-09 北方工业大学 基于模糊阈值的自适应图像分割方法
CN106340029A (zh) * 2016-08-23 2017-01-18 湖南文理学院 基于Beta‑Gamma散度的灰度图像阈值分割方法
CN107590815A (zh) * 2017-09-07 2018-01-16 陕西师范大学 基于萤火虫群优化法的图像多阈值分割方法

Patent Citations (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN101710425A (zh) * 2009-12-25 2010-05-19 南京航空航天大学 基于图像灰度梯度和灰度统计直方图的自适应预分割方法

Non-Patent Citations (2)

* Cited by examiner, † Cited by third party
Title
Analysis of SEM digital images to quantify crack network pattern area in chromium electrodeposits;M.Vidal等;《Surface and Coatings Technology》;20160115;全文 *
Characterization of the effect of normal load on the discontinuity morphology in direct shear specimens using X-ray micro-CT;Bryan S. A. Tatone等;《Acta Geotechnica》;20140506;全文 *

Also Published As

Publication number Publication date
CN110021030A (zh) 2019-07-16
US20210201497A1 (en) 2021-07-01
WO2020177215A1 (zh) 2020-09-10

Similar Documents

Publication Publication Date Title
CN110021030B (zh) 一种岩土体材料数字图像的分割阈值确定方法
Xu et al. An investigation into the relationship between saturated permeability and microstructure of remolded loess: A case study from Chinese Loess Plateau
CN108763711B (zh) 一种基于岩心扫描图像分块数值模拟的渗透率预测方法
CN105352873B (zh) 页岩孔隙结构的表征方法
Liu et al. Multi-scale fractal analysis of pores in shale rocks
Kane et al. Microstructural characterization and pore structure analysis of nuclear graphite
US9558588B2 (en) Method for building a 3D model of a rock sample
CN105445160B (zh) 一种沥青混合料的空隙特征及其提取方法
CN113706603B (zh) 页岩有机质孔隙连通性分类表征的方法
Liu et al. A new method for threshold determination of gray image
CN105649615A (zh) Ct定量、三维可视化测试储层致密油赋存状态的方法
CN108061697B (zh) 土体三维孔隙率计算方法
Al Ibrahim et al. An automated petrographic image analysis system: Capillary pressure curves using confocal microscopy
CN104268830A (zh) 基于数字图像确定不均匀岩土材料渗透系数的方法
CN102841220A (zh) 一种基于孔隙率的黏土扫描电镜照片图像分割方法
CN113295720B (zh) 一种利用ct扫描进行微米级矿物识别的装置及方法
Anangsha et al. A new autonomous program customized for computing surface cracks in an unsaturated soil in a 1-D column
Shi et al. A comprehensive assessment of image processing variability in pore structural investigations: conventional thresholding vs. machine learning approaches
CN104217424A (zh) 冰湖终碛垄颗粒的分析方法
CN110988004B (zh) 一种粒间压溶成因石英胶结物含量的评价方法
CN114518309A (zh) 一种页岩油储层中夹层孔隙有效三维表征的方法及装置
Moraru et al. Euler number: a method for statistical analysis of ancient pottery porosity
Plaisted et al. Testing of expansive clays in a centrifuge permeameter
Bolshakov et al. Investigation of the pore space structure by a scanning electron microscope using the computer program collector
Moraru et al. Ancient pottery analysis using SEM image processing

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