基于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