基于matlab的岩石表面裂纹统计方法

简文星，余锦风，任 佳，熊亚萍（.中国地质大学（武汉）工程学院，湖北　武汉430074；2.中国地质大学（武汉）教育部长江三峡库区地质灾害研究中心，湖北　武汉430074）

基于matlab的岩石表面裂纹统计方法

简文星1，2，余锦风1，任 佳1，熊亚萍1
（1.中国地质大学（武汉）工程学院，湖北武汉430074；2.中国地质大学（武汉）教育部长江三峡库区地质灾害研究中心，湖北武汉430074）

利用matlab图像处理工具箱设计开发了岩石裂纹检测系统，实现了对岩石表面裂纹的识别检测以及获取裂纹面积和裂隙率两个特征参数。检测系统中图像处理分为图像预处理、裂纹提取和计算三个阶段。首先对裂纹照片进行灰度处理并滤波；然后利用扩充Sobel算子进行边缘检测，突出裂纹与背景的差别，再进行图像分割，将灰度图像变为二值图像；最后统计计算特征参数。采用该检测系统对“浸泡—风干”循环作用下巴东组第四段泥质粉砂岩裂纹扩展情况进行统计分析，结果表明:随着泥质粉砂岩“浸泡—风干”循环次数的增加，裂纹的面积和裂隙率明显增大；初期阶段裂纹的扩展速率较快，后期逐渐减缓。岩石裂纹检测系统能够直接反映裂纹的扩展情况，且操作简单、结果合理，具有较强的实用性，为裂纹扩展规律研究提供了新的手段。

matlab数字图像处理；岩石裂纹检测系统；“浸泡—风干”循环作用；裂纹面积；裂隙率；裂纹扩展规律

三峡库区库水位周期性升降作用使得库岸岩体长期处于“浸泡—风干”循环环境中，“浸泡—风干”循环作用使得库岸岩体结构发生变化，从而导致力学性质劣化，直接影响了库岸的稳定性。“浸泡—风干”循环作用对岩石来说是一种“疲劳作用”［1］，由于岩石内部存在微小裂隙，经过多次“浸泡—风干”循环后，会发展成宏观的裂纹。同时，“浸泡—风干”循环作用对砂岩的次生孔隙的增长有重要的促进作用［2］。裂纹的扩展是导致岩石性质劣化的直接因素之一。库水作用下岩体原有裂隙的扩张及产生新的裂隙都会增加滑坡的变形［3］，因此分析岩石裂纹的扩展规律具有重要的现实意义。

针对岩石裂纹的扩展规律，许多学者进行了大量深入的研究。如陈四利等［4］利用CT扫描研究了三轴加载条件下的岩石的微裂纹萌生、扩展和破裂过程；刘冬梅等［5］利用实时全息干涉法、数字摄像机和图像处理系统对单轴和压剪荷载作用下岩石进行了连续动态观测，描述了岩石裂纹的起裂、扩展和变形过程，有效地判识了岩石裂纹的力学性状；刘涛影等［6］总结了压剪岩石裂纹在高渗透压作用下的起裂规律和分支裂纹尖端强度因子的变化规律，并提出了高渗压作用下有关压剪岩石裂纹渐进破坏的损伤方程；师访等［7］利用有限元法模拟了单轴压缩条件下裂纹的扩展过程，并提出了可以预测张拉、压剪交互作用下岩石的裂纹的萌生和开裂扩展的方案；赵兴东等［8］利用声发射定位技术监测岩石裂纹扩展演化过程，预测了岩样内部裂纹在空间的位置、扩展方向和曲面形态；郑贵平等［9］利用波速对各类岩石进行测试，结果显示裂纹的扩展影响了波速随应力变化的规律。

综上研究可见，学者们采用了不同的方法（如数字模拟、CT扫描、波速测试、声发射等）研究了岩石裂纹的扩展规律，但这些方法需要一定的技术手段和实验条件。近年来，在国外一些研究者采用数字图像处理技术研发了一些裂纹检测系统，如道路表面裂纹自动检测系统［10］、混凝砖表面裂纹的提取与测量系统［11］和地下管道表面裂纹检测系统［12］；在国内也有很多研究者提出了一些效果较好的裂纹检测算法，如程骏等［13］在机器视觉技术基础上，提出了检测耐火砖裂纹的方法；张洪光等［14］利用二值化算法和Agent算法相结合进行路面图像的裂纹检测；郑嵘等［15］则利用形态学和小波变换相结合检测板坯表面图像上的纵向裂纹。

本文基于matlab数字图像处理设计开发了岩石裂纹检测系统，相对以前的技术方法可更加直观地反映裂隙的发展情况，并可将裂隙的微小变化通过图像数据变化直接展现出来。通过实例将该检测系统应用到“浸泡—风干”循环作用下岩石表面裂纹扩展规律研究上，结果表明其操作简单、结果合理，为裂纹扩展规律的研究提供了一种新方法和新思路。

1　岩石裂纹检测系统

1.1岩石裂纹检测系统简介

岩石裂纹检测系统的界面设计比较简单，且操作方便，并可以直观地反映裂纹的扩展情况。该检测系统主要由图形显示界面、图像导入界面和计算界面三个模块组成，其界面视图如图1所示。

图1　岩石裂纹检测系统界面Eig.1 Detection system of rock cracks

岩石裂纹检测系统最核心的部分就是图像的各种处理操作，其包括图像预处理、裂纹提取和计算三个阶段（见图2）:①图像预处理，首先对裂纹照片（RGB图像）进行灰度处理，然后滤波，消除部分噪声，为后续处理提供基础；②裂纹提取，利用扩充Sobel算子进行边缘检测，突出裂纹与背景的差别，再进行图像分割，将灰度图像变为二值图像；③计算，对裂纹的面积和裂隙率两个特征参数进行统计计算。

图2　裂纹图像处理基本流程Eig.2 Basic process of crack image processing

1.2图像预处理

岩石图像在获取过程中，光线强度和方向不一，加上拍摄设备可能发生微小的抖动，会造成图像一定程度失真，需要对图像进行滤波处理，消除噪声。图像处理常见的噪声分为乘性噪声、加性噪声、量化噪声、盐和胡椒噪声。本次采集的岩石图像的噪声主要来源于图像的获取及传输，绝大部分噪声属于加性噪声，因此在图像预处理阶段，选取了消除加性噪声效果好的中值滤波器。

中值滤波是一种基于排序统计理论的非线性平滑滤波信号处理技术。中值滤波主要思路为:①将滤波模板在待处理图像中漫游，并将模板的中心与图像某个像素点重合；②读取模板中各对应像素的灰度值；③将上一步骤中灰度值从小到大排序；④取上一步骤中的中间数据赋值到对应模板中心位置的像素。若模板中有奇数个元素，中值取排序后的中间元素灰度值，否则取中间两个元素灰度的平均值。在一定条件下，中值滤波可以克服线性滤波器处理图像时细节模糊的问题。

利用检测系统的预处理功能处理原始图像，可得到去除噪声后的灰度图像，见图3。

图3　去噪声后的灰度图像Eig.3 Grayscale image after eliminating noise

1.3裂纹提取

裂纹提取是图像处理的关键步骤，主要是把岩石的裂纹从背景图像中提取出来。然而岩石裂纹主要特征表现在裂纹的边缘，因此首先要对岩石裂纹图像进行精确的边缘检测，以便于图像处理后续步骤的顺利进行。

1.3.1边缘检测

图像的边缘是其最基本的特征之一。边缘表明了图像局部特征的不连续性，灰度或结构信息的突变（如灰度级、颜色、纹理结构）就是边缘。matlab数字图像处理中定义了很多边缘检测算子，包括Robert算子、Sobel算子、Prewitt算子、LOG算子和Laplacian算子等。由于岩石裂纹图像噪声多，因此本文选取Sobel算子对图像边缘进行检测。Sobel算子的两个模板如下:

传统的Sobel算子只有水平和垂直方向两个模板，它们对水平和垂直方向的裂纹图像产生最大响应，但并不表示只能检测水平裂纹和垂直裂纹，它们也能检测非水平或非垂直裂纹。但是由于岩石裂纹的发展往往是不规则的，Sobel算子对裂纹的检测具有一定的局限性，对于其他方向裂纹检测效果不是很好，因此需要对Sobel算子进行改进。具体改进方法如下:

首先，基于边缘检测的基本理论构建45°和-45°的Sobel算子模板，两个模板如下:

然后，根据模板计算边缘检测的响应值，即进行卷积运算，选用通用模板S，再利用下式计算出响应值R:

式中:wi代表模板系数；zi代表中心像素点邻域点的灰度值。

在具体编制边缘检测程序时，根据上述方法在中心像素点（x，y）处分别计算水平、垂直、45°和-45°Sobel算子模板的响应值R1、R2、R3和R4。响应值有正有负，但我们只对响应幅值感兴趣。若|Ri|＞|Rj|（i≠j），则取Ri作为中心像素点（x，y）处的输出。边级检测后典型输出图像见图4，图像显示的裂纹与背景之间有较大的视觉差距，便于后续处理运算。

1.3.2图像分割

阈值分割方法，也叫图像二值化，是一种最简单的图像分割方法，它的关键在于找到合适的阈值。阈值分割法基本原理可用下式表示:

式中:Z表示阈值，为图像灰度级内的任意一个灰度级的集合；ZE和ZB分别表示对象和背景的灰度级。

合适的阈值是图像分割成功的关键，对象和背景之间的差距越大，图像分割的效果越好。前面使用扩充Sobel算子提取岩石裂纹增大了裂纹与背景之间灰度级的差距，因此可以利用阈值分割法来处理图像。

图4　边缘检测后的图像Eig.4 Image after edge detection

迭代式阈值分割法比较常用，图像分割效果也较好，因此本检测系统采用了迭代阈值分割法。迭代阈值分割法是将灰度图像中最大灰度值Zmax和最小灰度值Zmin的平均值作为初始的阈值T0［T0= （Zmax+Zmin）/2］，根据阈值Tk将灰度图像分割成两部分区域即对象区域R1和背景区域R2，分别计算两个区域的灰度均值μ1和μ2，再重新计算新阈值Tk=（μ1+μ2）/2，重复上述步骤直到Tk不再变化为止。迭代阈值分割法实现算法见图5。

图5　迭代阈值分割法Eig.5 Iterative threshold segmentation algorithm

利用检测系统的图像分割功能处理边缘检测图像得到的二值图像见图6。由图6可见，图像中很清楚地显示了原始图片中的裂纹，这为统计裂纹的特征参数提供了良好的基础。

图6　图像分割后二值图像Eig.6 Binary image after image segmentation

1.4计算

图像的特性研究主要针对研究所关注的对象。对图像进行定量的参数研究对于准确把握图像特性有重要的意义。检测系统将裂纹的面积（像素尺度）和裂隙率作为图像的定量特征。

统计裂纹面积的方法有很多，包括逐像素提取法、垂直方向数据提取法和13点提取法等。本检测系统选取了逐像素提取法。经图像分割所得的二值化图像，裂纹用暗点表示即为“1”，背景用亮点表示即为“0”，因此逐像素提取法可以针对这一特征，对亮点进行统计，其计算公式如下:

式中:A为裂纹的面积；R为图像区域的点；（x，y）为像素点的坐标；f（x，y）为像素点的值。

对于二值图像来说，裂纹像素点的值f（x，y）为1，背景像素点的值f（x，y）为0，所以可以用公式（6）统计裂纹的面积。

为了研究裂纹的发展规律，我们在像素的尺度定义一个裂隙率的参量，来反映裂纹的相对大小。裂隙率指的是裂纹的面积（像素）与整个图片的面积（像素）之比，其表达式为

式中:A0为裂纹的面积（像素）；A为整个图像的面积（像素）。

2　实例应用分析

三峡库区库水位周期性升降作用使得库岸岩体长期处于“浸泡—风干”循环环境中，“浸泡—风干”循环作用使得库岸岩体的结构发生变化，从而导致力学性质劣化，直接影响了库岸的稳定性。本文采用上述岩石裂纹检测系统对“浸泡—风干”循环作用下巴东组第四段泥质粉砂岩的表面裂纹扩散情况进行了统计分析，以探究岩石裂纹的扩展规律。

2.1岩样的制作

本次岩样采用三峡库区巴东组第四段泥质粉砂岩，经过中国地质大学（武汉）国家重点实验室X射线衍射分析，确定其主要成分为石英、绿泥石、伊利石和方解石等，其中黏土矿物含量高达40%以上，具体成分含量见表1。为了便于统计试验结果，将岩石在长江科学研究院岩石样品加工实验室统一制作成直径为50 mm、高度为100 mm的圆柱形标准试样，见图7。

表1　岩样X射线衍射分析结果Table 1 Test results of X-ray diffraction

图7　加工后的岩样Eig.7 Processed rock samples

2.2试验方案

三峡水库正常运营过程中，库岸消落带岩石处于“浸泡—风干”循环状态，这种状态会导致岩石的劣化，直接体现在岩石裂纹的扩展。为了模拟库水位升降作用，本文采取了“浸泡—风干”循环试验进行研究。具体试验步骤为:①在水中浸泡24 h；②在空气中自然风干24 h。以上两个过程为一个周期，具体试验流程见图8。

将制作的泥质粉砂岩标准试样进行筛选，剔除残缺和带裂纹样品，并选择同一块岩石切割成的试件5块，编号为A1、A2、A3、A4、A5。

2.3岩石裂纹特征参数统计与分析

在开发裂纹检测系统过程中，为了定量地确定岩石裂纹的特征，本文选定两个特征参数即裂纹的面积和裂隙率作为依据，将裂纹的面积（像素尺度）和裂隙率作为图像的定量特征。

按照图8试验流程，利用检测系统对试件A1上截面的裂纹特征参数进行了统计，其统计结果见表2。典型检测结果界面视图见图9。

图8　试验流程图Eig.8 Experiment process

表2　岩石裂纹特征参数统计结果Table 2 Characteristic parameters statistics of rock cracks

为了更好地了解“浸泡—风干”循环作用下裂纹的阶段增长情况，本文依据表2数据绘制了裂隙率随“浸泡—风干”循环次数的关系曲线，见图10。

综合表2和图10分析可知:

（1）随着“浸泡—风干”循环次数的增加，裂纹面积和裂隙率都有所增加，且在第十一次循环时达到最大值，裂纹面积为18 409像素，裂隙率为0.41%，说明“浸泡—风干”循环作用对于岩石裂纹的扩展起着不可忽略的作用，然而裂纹的扩展势必会带来岩石强度的降低，直接导致了岩石的劣化。

图9　岩样裂纹检测系统检测结果界面视图Eig.9 Interface of the detection results from the rock crack detection system

图10　岩石裂纹在“浸泡—风干”循环作用下裂隙率的变化曲线Eig.10 Change curve of fissure percentage under the“immersion-air dry”circulation

（2）在“浸泡—风干”循环作用的不同时期，裂纹的扩展速率不同。经历前三次“浸泡—风干”循环的样品，裂隙率阶段增长量分别达到了0.07%、0.06%和0.08%，说明在“浸泡—风干”循环的初期阶段，岩石裂纹的扩展速率较快；而在第四次“浸泡—风干”循环时，裂隙率阶段增长量下降到了0. 01%，在后期最高的裂隙率阶段增长量为0.04%。从整个过程来看，岩石裂纹的发展呈现初期阶段扩展速率较快，后期逐渐趋于平缓的趋势。

3　结 论

（1）基于matlab数字图像处理开发设计了岩石裂纹检测系统，该检测系统界面设计简单、操作方便，能够将背景和裂纹严格区分，并将裂纹的微小变化通过定量的图像数据变化进行统计，可获取岩石表面裂纹的基本特征参数，能够分析裂纹扩展趋势，具有一定的实用价值。

（2）实例应用表明:随着泥质粉砂岩“浸泡—风干”循环次数的增加，裂纹面积和裂隙率明显增大，岩石裂纹的发展呈现初期阶段扩展速率较快、后期逐渐减缓的趋势。

［1］徐千军，陆杨.干湿交替对边坡长期安全性的影响［J］.地下空间与工程学报，2005，1（7）:1021-1024.

［2］邓华锋，李建林，王孔伟，等.饱和-风干循环过程中砂岩次生孔隙率变化规律研究［J］.岩土力学，2012，33（2）:483-488.

［3］胡畅，牛瑞卿.三峡库区树坪滑坡变形特征及其诱发因素研究［J］.安全与环境工程，2013，20（2）:41-45.

［4］陈四利，冯夏庭，周辉.化学腐蚀下砂岩三轴细观损伤机理及损伤变量分析［J］.岩土力学，2004，25（9）:1363-1367.

［5］刘冬梅，蔡美峰，周玉斌，等.岩石裂纹扩展过程的动态监测研究［J］.岩石力学与工程学报，2006，25（3）:467-472.

［6］刘涛影，曹平，章立峰，等.高渗压条件下压剪岩石裂纹断裂损伤演化机制研究［J］.岩土力学，2012，33（6）:1801-1808.

［7］师访，高峰，李玺茹，等.模拟岩石压剪状态下主次裂纹萌生开裂的扩展有限元法［J］.岩土力学，2014，35（6）:1809-1817.

［8］赵兴东，李元辉，袁瑞甫，等.基于声发射定位的岩石裂纹动态演化过程研究［J］.岩石力学与工程学报，2007，26（5）:944-950.

［9］郑贵平，赵兴东，刘建坡，等.岩石加载过程声波波速变化规律实验研究［J］.东北大学学报（自然科学版），2009，30（8）:1197-1200.

［10］Pynn J，Wright A，Lodge R.Automatic identification of cracks in road surfaces［C］//Seventh International Conference on Image Processing and Its Applications［S.l.］:IEEE，1999，2:671-675.

［11］Ito A，Aoki Y，Hashimoto S.Accurate extraction and measurement of fine cracks from concrete block surface image［C］// IEEE 28th Annual Conference of the Industrial Electronics Society（IECON 02）［S.l.］:IEEE，2002，3:2202-2207.

［12］Eieguth P W，Sinha S K.Automated analysis and detection of cracks in underground scanned pipes［C］//Proceedings of 1999 International Conference on Image Processing［S.l.］:IEEE，1999，4:395-399.

［13］程骏，贺俊吉，肖建立，等.基于机器视觉的耐火砖裂纹检测［J］.微计算机信息，2010（35）:233-235.

［14］张洪光，王祁.基于人工种群和Agent的路面裂纹检测算法［J］.哈尔滨工业大学学报，2007，39（1）:128-132.

［15］郑嵘，方志宏，夏勇.板坯纵向裂纹的图像检测方法［C］//中国计量协会冶金分会.2012年会暨能源计量与节能降耗经验交流会论文集.北京:《冶金自动化》杂志社，2012.

The Statistical Method of Rock Surface Crack Based on Matlab

JIAN Wenxing1，2，YU Jinfeng1，REN Jia1，XIONG Yaping1
（1.Faculty of Engineering，China University of Geoscience，Wuhan 430074，China；2.Three Gorges Research Center for Geo-hazards，Ministry of Education，China University of Geosciences，Wuhan 430074，China）

Based on matlab image processing toolbox，this study develops the rock crack detection system to achieve the surface crack identification and detection of rock，and obtain the characteristic parameters of the crack including area and fissure percentage.Image processing of the detection system includes three stages:image preprocessing，cracks extraction and calculation.Eirstly，the photos of rock surface cracks（RGB images）is transformed into the gray image and filtered；Secondly，the crack edge is detected with Sobel operator to highlight the difference between cracks and background；And then image segmentation is carried out to make the gray-scale image become a binary image；Lastly，characteristic parameters are counted.This paper applies the system to summarize the pelitic siltstone crack propagation behavior in the 4th member section of Badong formationgroup under the conditions of“immersion-air dry”circulation.The results show that with the number of“immersion-air dry”circulation increasing，the crack area and fissure percentage tend to increase significantly，and the growth rate of cracks in the initial stage is fairly fast while that in the later period gradually slow down.Rock crack detection system is practical and reasonable to reflect cracks expanding.It can provide a new means for the study of crack propagation behavior.

digital image processing in matlab；rock crack detection system；immersion-air dry circulation；crack area；fissure percentage；crack propagation behavior

X93；TU45

A

10.13578/j.cnki.issn.1671-1556.2015.05.023

1671-1556（2015）05-0128-06

2015-04-03

2015-05-30

国家自然科学基金项目（41272306）

简文星（1967—），男，博士，教授，主要从事工程地质与岩土工程等方面的研究。E-mail:wxjian@cug.edu.cn