CN110163826A - 一种黄土微结构图像的处理方法 - Google Patents
一种黄土微结构图像的处理方法 Download PDFInfo
- Publication number
- CN110163826A CN110163826A CN201910435101.4A CN201910435101A CN110163826A CN 110163826 A CN110163826 A CN 110163826A CN 201910435101 A CN201910435101 A CN 201910435101A CN 110163826 A CN110163826 A CN 110163826A
- Authority
- CN
- China
- Prior art keywords
- image
- loess
- microstructure
- hole
- value
- 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.)
- Withdrawn
Links
- 238000003672 processing method Methods 0.000 title claims abstract description 8
- 238000012545 processing Methods 0.000 claims abstract description 11
- 238000004458 analytical method Methods 0.000 claims abstract description 10
- 230000000877 morphologic effect Effects 0.000 claims abstract description 6
- 238000005259 measurement Methods 0.000 claims abstract description 4
- 238000000034 method Methods 0.000 claims description 52
- 238000009826 distribution Methods 0.000 claims description 33
- 239000011148 porous material Substances 0.000 claims description 29
- 239000002689 soil Substances 0.000 claims description 28
- 239000002245 particle Substances 0.000 claims description 24
- 241000220317 Rosa Species 0.000 claims description 15
- 230000009466 transformation Effects 0.000 claims description 10
- 238000013139 quantization Methods 0.000 claims description 5
- 238000002360 preparation method Methods 0.000 claims description 3
- 238000013461 design Methods 0.000 claims description 2
- 238000001878 scanning electron micrograph Methods 0.000 claims 1
- 230000000694 effects Effects 0.000 abstract description 16
- 238000003709 image segmentation Methods 0.000 abstract description 15
- 238000007781 pre-processing Methods 0.000 abstract description 5
- 230000002708 enhancing effect Effects 0.000 abstract description 3
- 230000003321 amplification Effects 0.000 description 11
- 238000003199 nucleic acid amplification method Methods 0.000 description 11
- 230000008569 process Effects 0.000 description 11
- 230000011218 segmentation Effects 0.000 description 7
- 238000001914 filtration Methods 0.000 description 6
- 230000008859 change Effects 0.000 description 4
- 230000003014 reinforcing effect Effects 0.000 description 4
- 238000011160 research Methods 0.000 description 4
- 238000005516 engineering process Methods 0.000 description 3
- 230000006872 improvement Effects 0.000 description 3
- 230000001965 increasing effect Effects 0.000 description 3
- 238000012360 testing method Methods 0.000 description 3
- 230000008901 benefit Effects 0.000 description 2
- 230000005540 biological transmission Effects 0.000 description 2
- 238000005520 cutting process Methods 0.000 description 2
- 238000010894 electron beam technology Methods 0.000 description 2
- 239000012530 fluid Substances 0.000 description 2
- 238000012804 iterative process Methods 0.000 description 2
- 238000012986 modification Methods 0.000 description 2
- 230000004048 modification Effects 0.000 description 2
- 238000005192 partition Methods 0.000 description 2
- 230000009467 reduction Effects 0.000 description 2
- 238000000926 separation method Methods 0.000 description 2
- 230000000007 visual effect Effects 0.000 description 2
- 241001270131 Agaricus moelleri Species 0.000 description 1
- 244000035744 Hura crepitans Species 0.000 description 1
- 244000137852 Petrea volubilis Species 0.000 description 1
- 230000001133 acceleration Effects 0.000 description 1
- 230000009286 beneficial effect Effects 0.000 description 1
- 230000015572 biosynthetic process Effects 0.000 description 1
- 230000000052 comparative effect Effects 0.000 description 1
- 230000007812 deficiency Effects 0.000 description 1
- 238000011161 development Methods 0.000 description 1
- 230000018109 developmental process Effects 0.000 description 1
- 238000009792 diffusion process Methods 0.000 description 1
- 235000013399 edible fruits Nutrition 0.000 description 1
- 230000005611 electricity Effects 0.000 description 1
- 238000011156 evaluation Methods 0.000 description 1
- 238000000605 extraction Methods 0.000 description 1
- 239000004744 fabric Substances 0.000 description 1
- 238000011049 filling Methods 0.000 description 1
- 238000003384 imaging method Methods 0.000 description 1
- 238000009413 insulation Methods 0.000 description 1
- 230000003993 interaction Effects 0.000 description 1
- 230000001788 irregular Effects 0.000 description 1
- 239000000463 material Substances 0.000 description 1
- 230000003287 optical effect Effects 0.000 description 1
- 230000035699 permeability Effects 0.000 description 1
- 239000011435 rock Substances 0.000 description 1
- 210000003296 saliva Anatomy 0.000 description 1
- 238000005070 sampling Methods 0.000 description 1
- 230000005476 size effect Effects 0.000 description 1
- 238000001228 spectrum Methods 0.000 description 1
- 238000010408 sweeping Methods 0.000 description 1
- 238000003786 synthesis reaction Methods 0.000 description 1
- WFKWXMTUELFFGS-UHFFFAOYSA-N tungsten Chemical compound [W] WFKWXMTUELFFGS-UHFFFAOYSA-N 0.000 description 1
- 229910052721 tungsten Inorganic materials 0.000 description 1
- 239000010937 tungsten Substances 0.000 description 1
- 238000001771 vacuum deposition Methods 0.000 description 1
Classifications
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
- G06T5/00—Image enhancement or restoration
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
- G06T5/00—Image enhancement or restoration
- G06T5/40—Image enhancement or restoration using histogram techniques
-
- 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
- G06T7/00—Image analysis
- G06T7/0002—Inspection of images, e.g. flaw 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/155—Segmentation; Edge detection involving morphological operators
Landscapes
- Engineering & Computer Science (AREA)
- Physics & Mathematics (AREA)
- General Physics & Mathematics (AREA)
- Theoretical Computer Science (AREA)
- Computer Vision & Pattern Recognition (AREA)
- Quality & Reliability (AREA)
- Image Processing (AREA)
Abstract
本发明公开了一种黄土微结构图像的处理方法，处理方法包括如下步骤：步骤1，制备黄土样品，采集黄土样品的图像；步骤2，对图像进行预处理；步骤3，对图像进行增强；步骤4，对图像进行分割；步骤5，对图像进行形态学处理。本发明通过average和laplacian滤波器综合使用增强效果明显；通过histeq函数可以扩展微结构图像灰度级，便于后期测量与分析；通过运用信息信息量数理统计能力的最大熵算法，图像分割效果明显。
Description
技术领域
本发明涉及一种黄土微结构图像的处理方法。
背景技术
天然土体是一种非连续、非均质的结构性材料，其宏观物理力学特性受土的结构的影响和控制，因此，对土的结构性的认识和研究是现代土力学的基础和核心。作为土的一种固有特性，结构性通过自身的变化规律来影响土的诸多工程特性，如压缩性、剪切强度、渗透性等，这就使得在工程实践中更难把握土性的变化。而获取高清的土样微结构图像、提取合理的微结构参数和量化参数的具体方法是研究土的结构性的前提。因此可见，提出高效解决土样微结构获取及微结构参数量化的方法具有非常重要的意义。
现代微观结构研究寻求将土体微观结构参数直接用于土体力学本构模型的建立。但是总体来说，关于微观结构的定量研究无论是技术上还是原理方面都还存在一些困难，如信息提取技术问题，土体在外力或其他因素作用下的结构变化，微观结构研究的尺寸效应，微观力学模型的建立等。如何将黄土的微观结构特征与宏观强度和变形特性建立联系仍处于探索阶段，如何建立考虑黄土微观结构特性参数的动本构模型，并应用于解决工程实际，这些研究更是处于探索阶段，对解决实际工程问题还有相当的距离。
发明内容
针对现有技术的不足，本发明提供了一种黄土微结构图像的处理方法，包括如下步骤：
步骤1，制备黄土样品，采集黄土样品的图像；
步骤2，对图像进行预处理；
步骤3，对图像进行增强；
步骤4，对图像进行分割；
步骤5，对图像进行形态学处理，得到黄土微结构图像。
步骤1包括：将干燥的黄土土样制备成10mm×10mm×2mm的薄片，并对要观测的自然断面进行镀金处理，得到黄土样品，使用扫描电子显微镜获取黄土样品图像。
步骤2包括：使用matlab中的二维线性空间滤波器对图像进行去噪。
步骤3包括：
步骤3-1，对图像进行灰度变换；
步骤3-2，对图像对比度进行调整；
步骤3-3，对图像进行直方图处理。
步骤4包括：基于最大熵法设定阈值T，将图像中灰度值大于或等于阈值T的区域的灰度值赋值为1，小于阈值T的区域的灰度值赋值为0，从而将黄土微结构中的颗粒和孔隙分割开。
步骤5包括：使用IPP(Image-Pro Plus)软件滤波组件中流域分割法Watershed和变细Thinning工具，识别并分离图像中相互重叠的部分，得到黄土微结构图像。
本发明还包括步骤6，对黄土微结构图像量化分析，具体包括如下步骤：
步骤a1，根据黄土微结构图像，获取黄土微结构孔隙分布特征；
步骤a2，获取黄土微结构孔隙形态分布特征；
步骤a3，分析黄土微结构孔隙定向分布特征。
步骤a1包括：将黄土微结构孔隙分为四类：微孔隙、小孔隙、孔隙和大孔隙，并统计黄土微结构孔隙随时间和空间的变化规律，统计黄土微结构孔隙所占百分比，得到孔隙分布特征图。
步骤a2包括：
步骤a2-1，采用椭圆率和分维数两个微观参数描述黄土微结构孔隙形态分布特征，根据如下公式计算椭圆率：
K＝L/B
其中，L为黄土微结构孔隙的最长弦，B为垂直于最长弦的弦长；
步骤a2-2，对于规则图形的周长S与测量单位尺寸ε一次方成正比，面积A则与ε二次方成正比，则图形周长S和面积就可以写成S∝A1/2比例关系，对于在二维空间内不规则分形周长与面积关系表达式如下：
[S(ε)]1/D＝a0ε(1-D)/D[A(ε)]1/2
式中a0为与土颗粒形状有关的系数，取值0-1，图形越规则，a0取值越接近1，上式两边取对数得：
测得周长S数据后，绘制如下直线：
由所得直线能够得到：分维数D就是直线部分的斜率的倒数，将上式中面积A和周长S换成等效面积和等效周长，则存在以下关系：
log S(ε)＝D/2logA(ε)+C
式中，D为SEM图像(SEM图像通过扫描电子显微镜scanning electron microscope扫描黄土样品图像得到)的分维数，C为0～1的常数。
步骤a3包括：
步骤a3-1，采用概率熵和玫瑰图的方式分析黄土微结构孔隙定向分布特征，得到结构单元体，并根据如下公式计算概率熵Hm：
其中，Pi为结构单元体在一方位区中出现的概率，即为在一方位区上结构单元体的定向强度；n为在结构单元体排列方向[0～N]中等分的方位区数；N表示方位角，取值0～180°；
步骤a3-2，将黄土微结构孔隙的角均分9个区位，并统计出每个区位所占整体面积的百分比，绘制得到玫瑰图。
本发明首先通过MATLAB图像处理工具箱对图像进行去噪，增强及分割等处理，提高图像质量；然后在IPP软件中对处理后的黄土微结构图像进行测量和分析。
本发明的有益效果是：详细介绍了采用matlab和IPP相结合的方式对黄土微结构处理流程。研究发现：在图像增强方面，通过average和laplacian滤波器综合使用增强效果明显；通过histeq函数可以扩展微结构图像灰度级，便于后期测量与分析；通过运用信息信息量数理统计能力的最大熵算法，图像分割效果明显。
附图说明
下面结合附图和具体实施方式对本发明做更进一步的具体说明，本发明的上述或其他方面的优点将会变得更加清楚。
图1是KYKY-2800B型扫描电子显微镜。
图2a是100倍数下黄土微结构图像。
图2b是200倍数下黄土微结构图像。
图2c是400倍数下黄土微结构图像。
图2d是800倍数下黄土微结构图像。
图2e是1000倍数下黄土微结构图像。
图3a是原始图像。
图3b是average滤波效果。
图3c是laplacian滤波效果。
图4是正常直方图。
图5a是原始图像及与对应的直方图。
图5b是imadjust函数增强效果图及与对应的直方图。
图5c是adapthisteq函数增强效果图及与对应的直方图。
图5d是histeq函数增强效果图及与对应的直方图。
图6a是Oust算法图像分割结果。
图6b是迭代法图像分割结果。
图6c是Bernsen算法图像分割结果。
图6d是最大熵法图像分割结果。
图7a是识别出的图像。
图7b是分离后的图像。
图8a是青海西宁地区的黄土微结构图像。
图8b是甘肃永登地区的黄土微结构图像。
图8c是甘肃兰州地区的黄土微结构图像。
图8d是宁夏西吉地区的黄土微结构图像。
图8e是山西芮城地区的黄土微结构图像。
图9是孔隙分类直方图。
图10a是青海孔隙分布特征图。
图10b是永登孔隙分布特征图。
图10c是兰州孔隙分布特征图。
图10d是西吉孔隙分布特征图。
图10e是芮城孔隙分布特征图。
图11是孔隙椭圆率直方图。
图12a是青海孔隙椭圆率分布图。
图12b是永登孔隙椭圆率分布图。
图12c是兰州孔隙椭圆率分布图。
图12d是西吉孔隙椭圆率分布图。
图12e是芮城孔隙椭圆率分布图。
图13a是青海(4m)孔隙的周长与面积分维关系。
图13b是永登(4m)孔隙的周长与面积分维关系。
图13c是兰州(4m)孔隙的周长与面积分维关系。
图13d是西吉(4m)孔隙的周长与面积分维关系。
图13e是芮城(4m)孔隙的周长与面积分维关系。
图14是颗粒或孔隙方位。
图15a是青海孔隙定向分布玫瑰图。
图15b是永登孔隙定向分布玫瑰图。
图15c是兰州孔隙定向分布玫瑰图。
图15d是西吉孔隙定向分布玫瑰图。
图15e是芮城孔隙定向分布玫瑰图。
具体实施方式
下面结合附图及实施例对本发明做进一步说明。
步骤一：黄土微结构图像获取
(1)试验仪器
仪器采用的是KYKY-2800B型扫描电子显微镜，如图1所示。KYKY-2800B型扫描电镜是一种观测样品微结构的光学仪器，是通过电子束与样品的相互作用能够产生样品表面放大的形貌像，并按被扫描时序将样品的形貌像建立起来的。此显微镜下可以观测到的黄土微观结构有，土的颗粒形态、孔隙特征、接触关系及胶结类型等。具体的技术指标如下所示：
(1)分辨率：4.5nm(钨丝阴极)
(2)放大倍数：15X～250000X
(3)加速电压：0.1kV～30kV
(4)试样尺寸：Φ60mm(最大)
(5)样品台
X Y Z是样品台的几何尺寸、T表示倾斜角度的范围、R表示转动角度；
1)标准样品台：X＝Y＝50mm,Z＝25mm,T＝-5°～+90°,R＝360°
2)大样品台：X＝80mm,Y＝50mm,Z＝30mm,T＝0～+90°,R＝360°
3)高真空高温样品台：最高温度达1000℃
(6)扫描方式：全幅、面、线、点、选区双放大及屏幕分割
(7)真空系统：全自动高速扩散泵真空系统
(2)样品制备
制备电子扫描试样时，试样观察面不应磨制一个平面，这会破坏原始结构，从而失去观测的意义。因此，课题组采用人工掰断的方法，保持断面不受碰撞，保证观察到的图像能够真实反映原状土体的微观状态。具体操作如下：
1)试样制备
根据经验，较干燥的土样在制备扫描电镜试样时较为方便。将原状土样和动三轴试验后的试样晾置一段时间后，用手轻轻掰断原状试样，直至出现较为平整的自然断面为止。由于扫面电镜试样观测是二次成像的原理，试样越薄观察效果越明显，因此，为了便于观测分析，底部用砂纸磨平或用小刀削平，最后将试样制备成10mm×10mm×2mm的薄片。
2)样品充电
导电性良好的样品，可以直接进行观察。而黄土作为一种绝缘土体，直接进行观察，电子束易在土体表面积累电荷而影响到土体的观测效果，因此，将粘结好的黄土扫描电镜试样放置在真空镀膜机中，对将要观测的自然断面进行镀金处理。由于黄土自然断面凹凸不平，一般镀2-3次最佳。
3)图像拍摄
KYKY-2800B型扫描电子显微镜放大倍数在15-25万倍。过大或过小的放大倍数均会使得图像计算值与实际值产生偏差，因此合适的放大倍数选择对于微结构图像测取工作是非常重要的。图2a、图2b、图2c、图2d、图2e是观测试样放大倍数分别为100、200、400、800和1000倍的放大效果。由图可知，放大倍数不易过小，虽然较小的放大倍数观测的到区域较大，但黄土微结构的细节信息模糊甚至观测不到；放大倍数不宜过大，虽然部分细节能够观测清晰，但局部信息往往并不能代表整体的结构特点，会导致计算值与实际值偏差较大。经反复试验，为保证图像信息的丰富和完整，放大倍数在400-600倍之间较为适宜。
步骤二：图像预处理
由于图像采集时光线的强弱、电压频率的波动以及拍摄图像时拍摄视角、仪器使用等原因，常常造成像素空间关系发生变化、图像背景不均匀和对比度不明显的现象，因此需要对图像进行预处理。
(1)去噪
图像的获取和传输过程是数字图像噪声的主要来源。在图像获取过程中主要受图像传感器的工作情况的影响；传输信道受到噪声污染会导致图像噪声的产生。
人们结合实际图像的特点、噪声频谱分布规律以及统计特征，发展了多种多样的去噪方法。常用的滤波方法有：非线性中值滤波、线性均值滤波、小波滤波等。
针对黄土微结构图像自身的特点，本次选用matlab中二维线性空间滤波器，这些滤波器可以通过分发fspecial得到，该函数生成一个滤波模板w，语法为
w＝fspecial(’type’，parameters) (1)
通过fspecial函数，可以调用9种不同的滤波器，并且根据图像处理的实际效果对滤波器的类型、形状、尺寸等参数进行调整，直到产生满意的图像为止。
从图3a、图3b、图3c可知，经过average均值滤波器滤波后发现，原始图像高频部分被去掉，图像中的一些噪声得到了消除，但由于均值滤波器滤波会使图像产生钝化的效果，图像颗粒和孔隙边缘变得模糊不清晰。结合使用laplacian滤波器锐化的功能，提高了图像清晰度，达到了降噪的效果。
步骤三：图像增强
图像增强，顾名思义，就是感兴趣的信息加以增强图像，不感兴趣的部分加以抑制。
(1)灰度变换
灰度即使用黑色调表示物体。每个灰度对象的亮度值都在0％(白色)到100％(黑色)之间。灰度变换主要通过改变图像数据占原始灰度范围，对像素点进行处理，使其在视觉上得到改变。常用到的灰度变换方式有：图像反转、对数变换、伽马变换等。
(2)对比度调整
对比度调整实际是改变原始图像各部分的反差，通过对比度增强或对比度扩展，增强灰度值间的动态范围，将感兴趣部分的特征扩展，使之占据可现实灰度级更大的部分。经常用到的对比度调整有：对比度拉伸变换以及对数变换。
(3)直方图处理
直方图处理是图像增强中使用的最为广泛的方法，直方图展示了图像中像素的灰度分布。在理论上来说，一张曝光良好的照片在不同亮度级别下细节都应该是非常丰富的，各个亮度值上都有像素分布，像一座起伏波荡的小山丘，如图4所示。而实际的黄土微结构图像由于机器及人为等因素，灰度级配往往不是很均匀。在图像直方图处理过程中往往用到的方法有：直方图均衡化、直方图拉伸及直方图匹配。
本发明采用不同的增强图像方法，并对各自的效果进行比较。不难发现，图5a中的所采集到的图像偏暗，直方图的峰值偏向坐标轴左侧，这时较为困难的测量图像微观信息；通过imadjust函数对原始图像进行灰度变化和adapthisteq函数进行对比度调整，从图5b和图5c中不难发现，进过调整后的图像变亮，灰度级范围分别在55～250和45～220之间，较原始图像有明显增强，但直方图分布不均匀，没有达到预期的效果；图5d是经过histeq函数直方图均衡化后的图像，颗粒区域强度增强，孔隙区域强度降低，像素占有全部的灰度级且分布均匀，图像达到预期效果。因此，本发明以直方图处理方式主来增强主黄土微结构图像。
在黄土微结构图像中，本发明较为关心的是黄土颗粒与孔隙之间的界限是否足够清晰，便于分割处理。通过图像去燥和增强等预处理过程可以使得图像便于观测，然而对于图像去燥和增强过程通用的理论是不存在的，图像预处理过程的好坏均由观察者自身判断，图像质量的视觉评价也是一种高度主观的过程。因此，在图像预处理过程中，处理的方法不仅仅拘泥于一种或者某几种，而是多种方法相结合，真正达到图像质量改善的目的。
4、步骤四：图像分割
图像分割是图像处理过程中的关键技术。图像分割技术随着计算机技术的发展，又有了质的飞跃。图像分割过程中使用最为广泛的方法是阈值分割，其分割过程就是图像二值化的过程，主要原理如下：
在一幅灰度图像中给定一个阈值T，灰度值大于或等于阈值T的将其赋值为1，而小于阈值T的灰度值赋值为0，进而通过阈值可以将黄土微结构中的颗粒和孔隙分割开来。阈值分割的优点是速度快，计算简单。但最大的问题是阈值的选择，阈值选择过大会将部分颗粒划分成孔隙，而阈值选之过小又会将孔隙判为颗粒，因此，决定图像处理结果好坏重中之重是选择合理的阈值。常用到阈值分割几种方法有：
(1)Oust算法
Oust算法是由日本大津学者提出来的，他认为图像由前景和背景两部分构成。如果两部分的类间方差越大，则图像的差别越大。当图像的类间方差最大时，前景和背景就会基本分离。因此，通过选择合适的阈值使得两部分的方差最大，继而达到图像准确分割。
(2)迭代法
迭代法是运用逐渐逼近发求出灰度最大值Zmax和最小值Zmin。令初始阈值为其平均值T0，根据T0大小将图象分割为两个部分，分别计算两个区域内平均灰度值Z1和Z2，利用此平均值重复计算阈值T0，直至T0不再发生变化为止。迭代法运算速度快，但在迭代过程中，由于初始值设置不合理会造成迭代过程的失败。
(3)Bernsen算法
Bernsen算法是一种较为原始的算法，以当前像素为中心，计算某邻域内像素的最大值M与最小值N，并令两者的均值为T；如果M-N>S(初设置)，则当令阈值为T；反之，说明窗口所在区域的灰度级差别较小。若T>128，则灰度值为255，反之，灰度值为0。
(4)最大熵法
最大熵是一种简单有效的图像分割方法。T.Pun将熵的概念第一次引入到图像分割的过程中来，对于黄土微观图像来说，阈值分割后将图像分为颗粒和孔隙两部分，对图像信息量的量化进行衡量，分割后颗粒和孔隙细节越丰富，信息量越大，则图像分割效果越明显。
图6a、图6b、图6c、图6d给出了四种阈值分割效果对比图，图6a是Oust算法图像分割结果(T＝102)。图6b是迭代法图像分割结果(T＝109)。图6c是Bernsen算法图像分割结果(T＝135)。图6d是最大熵法图像分割结果(T＝90)。
Oust算法和迭代法计算得到的阈值相近，颗粒和孔隙得到明显分割，但由于这两种算法只考虑像素本身的灰度值，目标物与背景灰度不明显时，使得附着在大颗粒上的粘粒信息明显丢失；Bernsen算法仅较好的提取了颗粒和孔隙边缘信息，颗粒内部大部分被孔隙填充，与实际不符；借助于熵对图像信息量数理统计的能力，分割后得到的信息量大，颗粒和孔隙明显分离。综合对比几种阈值分割结果，借助于最大熵阈值处理方式，图像分割处理效果较好。
5、步骤五：图像形态学处理
在图像处理过程中，经常还用到形态学处理。形态学处理可以识别颗粒和孔隙的形状，比图像骨架、图像边界以及凹凸情况等等。通过数学形态学中开、闭、腐蚀、膨胀等基本变换，可以提取出结构单元体基本骨架结构。在IPP软件中几乎包含了形态学中所有的处理方式，并且可以处理结果进行实时调整，为图像处理带来便利条件。
图7中采用的是IPP软件滤波组件中流域分割法(Watershed)和变细(Thinning)工具，识别并分离结构图像中相互重叠的部分，图7a是识别出的图像，图7b是分离后的图像，与之前图像分割处理结果对比发现，重叠的结构单元体被分离出来，而且颗粒和孔隙之间的边界界限更加明显。
6、步骤六：土微结构量化方法
选取黄土高原五个地区黄土采样，按上述五个步骤得到如图8a、图8b、图8c、图8d、图8e所示的微结构图像；
(1)黄土微结构孔隙分布特征
将五处地区的黄土孔隙分为四类：微孔隙、小孔隙、孔隙和大孔隙，并统计了黄土孔隙随时间和空间的变化规律，其中五处地点的黄土均在4m和8m处取样，此外，兰州和芮城的还在12m深度进行取样，以作对照，统计结果如图9所示。
在图9中，统计了五处典型地区黄土微结构孔隙所占百分比。从上图可以看出，中孔隙在孔隙中占有举足轻重的地位，其所占面积百分比在40％～65％左右。这是因为大孔隙面积虽然较大，但数量较少，而中、小孔隙数量较多，因此，中孔隙孔隙面积较大。从空间上看，自西北向东南方向，大孔隙的含量逐渐增加，中、小孔隙特别是中孔隙的含量显着减小。
从图10a、图10b、图10c、图10d、图10e中可知，五处黄土的孔隙分布特征图均呈反“√”型，在孔隙分布图中的拐点在中孔隙处，表明中孔隙所占百分比最大。随着深度的增加，大孔隙面积逐渐减小，中孔隙面积所占比重逐渐增加。这是由于地应力随着深度增加而增大，孔隙发生了挤压和变形，使得大孔隙数量减少；部分小颗粒散落到大孔隙中，使得大孔隙转化为中孔隙和小孔隙，导致中孔隙的含量增加。
(2)微结构孔隙形态分布特征
为了进一步分析黄土孔隙形态分布特征，采用椭圆率和分维数两个微观参数进行描述。
1、椭圆率
K＝L/B (3)
式中：L为孔隙的最长弦，B为垂直于最长弦的弦长。
椭圆率是反映孔隙形态的一种基本参数。K值越大，表示孔隙形态越狭长，K值越小，表明孔隙形态越趋于圆形；
从图11得出，黄土孔隙的椭圆率在1.695-2.713左右，其中中孔隙和小孔隙的椭圆率相对较大，微孔隙的椭圆率最小。从椭圆率的直方图中可以看出，在孔隙形态方面，孔径较大的孔隙易受到挤压而发生变形，孔隙形态区域狭长形态，椭圆率较大；孔径较小的孔隙尤其是微孔隙，孔隙形态趋于稳定，不易受到外界压力而发生变形，椭圆率较小。
从图12a、图12b、图12c、图12d、图12e可以看出，随着深度的增加，微孔隙和小孔隙的椭圆率基本不变，而中孔隙和大孔隙变小，尤其是大孔隙椭圆率陡然降低。这是由于地应力随深度增加而增加，微孔隙和小孔隙孔隙结构较稳定，不易发生变形；而大孔隙受外力影响较大，大孔隙受力发生变形，向着结构较稳定的微、小孔隙发展。
2、分维数
分形的重要描述参数是分维，用来描述那些不规则物体的几何形态。维数的数值不一定是整数，其中追著名的是Hausdorff维数，用来表示具有分形特征物体或现象的维数。人们在解决实际问题时，主要用到的分维方法有计盒法、Sandbox法、回转半径法、面积～周长法等等。
本实施例在进行分维数计算采用的是面积～周长法。对于规则图形的周长S与测量单位尺寸ε一次方成正比，面积A则与ε二次方成正比，那么图形周长S和面积就可以写成S∝A1/2比例关系。对于在二维空间内不规则分形周长与面积关系表达式如下：
[S(ε)]1/D＝a0ε(1-D)/D[A(ε)]1/2 (4)
式中a0为与土颗粒形状有关的系数，取值0～1，图形越规则，a0取值越接近1，两边取对数得：
测得数据后，绘制如下直线：
由所得直线可知，分维数D就是直线部分的斜率的倒数。将上式中面积A和周长S换成等效面积和等效周长，则存在以下关系：
log S(ε)＝D/2logA(ε)+C公式(7)
式中，D为SEM图像的分维数，C为常数。
图13a、图13b、图13c、图13d、图13e给出了五处地点取样深度为4m的黄土微结构图像分维数图，从图中可以看出取对数后的周长和面积具有较好的相关性。表1对不同地点和不同深度的黄土孔隙分维数进行了统计，经对比发现，分维数在1.270-1.389左右。在空间上，自西北向东南，同一水平深度黄土的孔隙分维数有逐渐增大的趋势。随着深度的增加，分维数逐渐变小。
表1
(3)微结构孔隙定向分布特征
为了研究孔隙定向分布特征，采用概率熵和玫瑰图的方式对黄土微结构孔隙进行分析。
1、概率熵
根据颗粒的几何特征，胡润林等提出采用颗粒或孔隙最大弦长与水平扫描线之间的夹角确定颗粒或孔隙择优方位(如图14所示)，其范围在[0，π]。
概率熵的表达式为：
式中：Pi为结构单元体在某一方位区中出现的概率，即为在某一方位区上单元体的定向强度，结构单元体是指电镜扫描图像中明显为一体的土颗粒，如图2所示；n为在单元体排列方向[0～N]中等分的方位区数。
这里单元体的排列方向为0°～180°，以10°为单位等分，n＝18。Hm的取值为[0，1]，Hm越小，表明结构单元体的排列方向越一致，有序度越高；Hm越大，说明结构单元体排列越混乱，有序性越低。
由表1可以看出，原状土的孔隙排列无序，概率熵在0.927-0.992之间。随着深度的增加，熵逐渐变小。
2、玫瑰图
由于不同孔隙对土体结构的影响程度不同，因此，在对孔隙进行定向分析时孔隙面积的考虑是必不可少的。玫瑰图是通过将孔隙的角均分9个区位，并统计出每个区位所占面积百分比绘制得到的。玫瑰图可以较为直观的反映孔隙的分布特征。
从图15a、图15b、图15c、图15d、图15e可以看出，五个地点的原状黄土孔隙均表现出一定的定向性，在玫瑰图中显示为“尖锐”的特点，但每幅图并不固定分布在某一区域，例如青海的黄土孔隙主要集中在区位四和区位八，即角度分布在60°～80°和140°～160°之间；兰州和永登的黄土孔隙主要集中在区位四和区位九，即角度分布在60°～80°和160°～180°之间；西吉的黄土主要集中在区位三，即角度分布在40°～60°之间；而芮城的黄土主要集中在区位八，即角度分布在140°～160°之间。随着深度的增加，玫瑰图曲线趋于圆滑，说明孔隙在各个方位分布均匀。
本发明提供了一种黄土微结构图像的处理方法，具体实现该技术方案的方法和途径很多，以上所述仅是本发明的优选实施方式，应当指出，对于本技术领域的普通技术人员来说，在不脱离本发明原理的前提下，还可以做出若干改进和润饰，这些改进和润饰也应视为本发明的保护范围。本实施例中未明确的各组成部分均可用现有技术加以实现。
Claims (10)
1.一种黄土微结构图像的处理方法，其特征在于，包括如下步骤：
步骤1，制备黄土样品，采集黄土样品的图像；
步骤2，对图像进行预处理；
步骤3，对图像进行增强；
步骤4，对图像进行分割；
步骤5，对图像进行形态学处理，得到黄土微结构图像。
2.根据权利要求1所述的方法，其特征在于，步骤1包括：将干燥的黄土土样制备成10mm×10mm×2mm的薄片，并对要观测的自然断面进行镀金处理，得到黄土样品，使用扫描电子显微镜获取黄土样品图像。
3.根据权利要求2所述的方法，其特征在于，步骤2包括：使用matlab中的二维线性空间滤波器对图像进行去噪。
4.根据权利要求3所述的方法，其特征在于，步骤3包括：
步骤3-1，对图像进行灰度变换；
步骤3-2，对图像对比度进行调整；
步骤3-3，对图像进行直方图处理。
5.根据权利要求4所述的方法，其特征在于，步骤4包括：基于最大熵法设定阈值T，将图像中灰度值大于或等于阈值T的区域的灰度值赋值为1，小于阈值T的区域的灰度值赋值为0，从而将黄土微结构中的颗粒和孔隙分割开。
6.根据权利要求5所述的方法，其特征在于，步骤5包括：使用IPP软件滤波组件中流域分割法Watershed和变细Thinning工具，识别并分离图像中相互重叠的部分，得到黄土微结构图像。
7.根据权利要求5所述的方法，其特征在于，还包括步骤6，对黄土微结构图像量化分析，具体包括如下步骤：
步骤a1，根据黄土微结构图像，获取黄土微结构孔隙分布特征；
步骤a2，获取黄土微结构孔隙形态分布特征；
步骤a3，分析黄土微结构孔隙定向分布特征。
8.根据权利要求7所述的方法，其特征在于，步骤a1包括：将黄土微结构孔隙分为四类：微孔隙、小孔隙、孔隙和大孔隙，并统计黄土微结构孔隙随时间和空间的变化规律，统计黄土微结构孔隙所占百分比，得到孔隙分布特征图。
9.根据权利要求8所述的方法，其特征在于，步骤a2包括：
步骤a2-1，采用椭圆率和分维数两个微观参数描述黄土微结构孔隙形态分布特征，根据如下公式计算椭圆率：
K＝L/B
其中，L为黄土微结构孔隙的最长弦，B为垂直于最长弦的弦长；
步骤a2-2，对于规则图形的周长S与测量单位尺寸ε一次方成正比，面积A则与ε二次方成正比，则图形周长S和面积就可以写成S∝A1/2比例关系，对于在二维空间内不规则分形周长与面积关系表达式如下：
[S(ε)]1/D＝a0e(1-D)/D[A(ε)]1/2
式中a0为与土颗粒形状有关的系数，上式两边取对数得：
测得周长S的数据后，绘制如下直线：
由所得直线能够得到：分维数D就是直线部分的斜率的倒数，将上式中面积A和周长S换成等效面积和等效周长，则存在以下关系：
logS(ε)＝D/2logA(ε)+C
式中，D为SEM图像的分维数，C为常数。
10.根据权利要求8所述的方法，其特征在于，步骤a3包括：
步骤a3-1，采用概率熵和玫瑰图的方式分析黄土微结构孔隙定向分布特征，得到结构单元体，并根据如下公式计算概率熵Hm：
其中，Pi为结构单元体在一方位区中出现的概率，即为在一方位区上结构单元体的定向强度；n为在结构单元体排列方向[0～N]中等分的方位区数，N表示方位角，取值0～180°；
步骤a3-2，将黄土微结构孔隙的角均分9个区位，并统计出每个区位所占整体面积的百分比，绘制得到玫瑰图。
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN201910435101.4A CN110163826A (zh) | 2019-05-23 | 2019-05-23 | 一种黄土微结构图像的处理方法 |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN201910435101.4A CN110163826A (zh) | 2019-05-23 | 2019-05-23 | 一种黄土微结构图像的处理方法 |
Publications (1)
| Publication Number | Publication Date |
|---|---|
| CN110163826A true CN110163826A (zh) | 2019-08-23 |
Family
ID=67632349
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN201910435101.4A Withdrawn CN110163826A (zh) | 2019-05-23 | 2019-05-23 | 一种黄土微结构图像的处理方法 |
Country Status (1)
| Country | Link |
|---|---|
| CN (1) | CN110163826A (zh) |
Cited By (2)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN110514671A (zh) * | 2019-09-24 | 2019-11-29 | 广西贺州市红星纸业有限公司 | 一种新型浆板洁净度检测方法 |
| CN114235837A (zh) * | 2021-12-20 | 2022-03-25 | 重庆大学 | 基于机器视觉的led封装表面缺陷检测方法、装置、介质及设备 |
Citations (2)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN108564124A (zh) * | 2018-04-13 | 2018-09-21 | 山东农业大学 | 一种基于支持向量机的稻瘟病菌孢子显微图像检测识别方法 |
| CN108830877A (zh) * | 2018-06-08 | 2018-11-16 | 中国石油化工股份有限公司 | 岩心体视显微图像定量描述方法 |
-
2019
- 2019-05-23 CN CN201910435101.4A patent/CN110163826A/zh not_active Withdrawn
Patent Citations (2)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN108564124A (zh) * | 2018-04-13 | 2018-09-21 | 山东农业大学 | 一种基于支持向量机的稻瘟病菌孢子显微图像检测识别方法 |
| CN108830877A (zh) * | 2018-06-08 | 2018-11-16 | 中国石油化工股份有限公司 | 岩心体视显微图像定量描述方法 |
Non-Patent Citations (2)
| Title |
|---|
| 徐世民 等: "基于Matlab和IPP的黄土孔隙微观结构研究", 《地震工程学报》 * |
| 徐世民: "考虑微结构的黄土动力特性试验研究", 《中国优秀硕士学位论文全文数据库 工程科技Ⅱ辑》 * |
Cited By (2)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN110514671A (zh) * | 2019-09-24 | 2019-11-29 | 广西贺州市红星纸业有限公司 | 一种新型浆板洁净度检测方法 |
| CN114235837A (zh) * | 2021-12-20 | 2022-03-25 | 重庆大学 | 基于机器视觉的led封装表面缺陷检测方法、装置、介质及设备 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| Heilbronner et al. | Image analysis in earth sciences: microstructures and textures of earth materials | |
| AU2019101789A4 (en) | Method for quantitatively measuring spatial structure of soil particulate organic matter | |
| CN110443793B (zh) | 一种沥青混合料空隙分布均匀性评价方法 | |
| CN109584286A (zh) | 一种基于广义回归神经网络的沥青路面构造深度计算方法 | |
| CN110163826A (zh) | 一种黄土微结构图像的处理方法 | |
| CN113706603B (zh) | 页岩有机质孔隙连通性分类表征的方法 | |
| Wan et al. | Wavelet-based statistical features for distinguishing mitotic and non-mitotic cells in breast cancer histopathology | |
| Gayathri et al. | A survey of breast cancer detection based on image segmentation techniques | |
| CN104318564A (zh) | 一种矿物颗粒分相的方法 | |
| CN114565658A (zh) | 基于ct技术的孔隙尺寸计算方法及装置 | |
| Lottin et al. | On methods determining the fractal dimension of combustion aerosols and particleclusters | |
| RP et al. | An effective SEM-based image analysis system for quantitative mineralogy | |
| CN108061697A (zh) | 土体三维孔隙率计算方法 | |
| Caron et al. | A method for detecting artificial objects in natural environments | |
| CN108205821A (zh) | 基于计算机视觉的工件表面三维重构方法 | |
| Tan et al. | A review of image processing and quantification analysis for solid oxide fuel cell | |
| Tovey et al. | Quantitative micro-porosity and orientation analysis in soils and sediments | |
| CN113390905B (zh) | 一种基于ct扫描技术量化土壤团聚体孔隙空间结构方法 | |
| Rosioru et al. | Morphometric assessments for the healthy rat hepatocytes | |
| CN116091465A (zh) | 一种基于三维扫描的织物起毛起球等级自动评定方法 | |
| JP3663465B2 (ja) | 画像情報を強調するための、ピクセルのデジタル式精密強度処理法 | |
| Bischin et al. | Computerized morphometric assessment of the human red blood cells treated with cisplatin | |
| Tovey et al. | Practical aspects of automatic orientation analysis of micrographs | |
| Wei et al. | Clustering-oriented multiple convolutional neural networks for optical coherence tomography image denoising | |
| CN114518309A (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 | ||
| WW01 | Invention patent application withdrawn after publication |
Application publication date: 20190823 |
|
| WW01 | Invention patent application withdrawn after publication |