基于Canny優化算法物理模擬實驗覆巖裂隙場自動識別研究

王順陽，賀易航，韓 博

(1.河南理工大學 測繪與國土信息工程學院，河南 焦作 454000； 2.河南省巖石礦物測試中心，河南 鄭州 450014)

大規模煤炭資源開采破壞了圍巖的應力平衡狀態，導致上覆巖層發生移動變形進而形成裂隙，影響煤礦的安全生產[1-2]。因此，研究采動覆巖裂隙場發育規律對于保障礦井安全生產，促進我國煤炭工業健康穩步發展具有現實意義。

通過閱讀國內外參考文獻發現，針對瀝青路面裂隙、巖土體裂隙等開展識別研究的比較多[3-5]。針對相似材料模擬實驗開展的研究主要可分為2個方面：①覆巖裂隙發育特征和時空演化規律研究。Wang等[6]采用相似材料模擬實驗和顆粒流數值模擬實驗，對孔隙變化和裂隙發育規律進行了研究，并與現場原位測試的結果進行了對比。Behrooz等[7-8]針對上組煤和下組煤存在部分重疊的2個長壁開采工作面采動覆巖形變特征進行了分析，采用包括近景攝影測量和激光掃描儀在內的監測手段對覆巖位移量進行了監測，并將監測結果與有限元法所得結果進行了對比。郭昌貴等[9]利用FLAC軟件建立煤層開采數值模擬模型，研究了不同高度的導水裂隙帶的發育規律。柴敬等[10]采用光纖光柵傳感器對采場覆巖的內部應變進行了監測，以此來研究關鍵層的初次破斷和周期破斷特征。②覆巖導水裂隙帶發育高度預測研究。Wang[11]基于關鍵層沉降的分析和評價提出了采動覆巖裂隙率的理論分布模型，并采用UDEC數值模擬方法證實了裂隙的“拱形”發育特征。此外，有關學者基于圖像分割對采礦過程中裂縫等問題進行研究，張昊[12]針對模糊C均值圖像分割方法存在大量噪聲干擾問題提出了一種融合空間噪聲信息的改進模糊C均值圖像分割方法，結果表明新的分割方法能夠實現圖像的精準分割。徐平等[13]針對采用閾值分割算法分割圖像時，背景存在干擾的問題，提出利用分割算法對圖像進行初始分割和二次分割，研究結果表明分割出的圖像邊界清晰能夠較好地抑制噪聲。鄒盛等[14]針對現階段煤礦帶式輸送機機輥缺乏有效的監測，傳統人工巡檢方法效率低、自動化程度低等問題，提出基于Hough變換的帶式輸送機托輥圖像分割方法，實驗結果表明該算法可以實現不同位置托輥的自動分割，能夠為煤礦安全開采帶來一定的技術基礎。邵佳園等[15]針對油基泥漿裂縫識別問題，選用具有各向異性的曲波變換對超聲波圖像進行了圖像分割，此方法能夠獲取圖像邊緣信息，對圖像分割有一定的作用。綜上，前人從不同行業背景和專業視角研究了采動覆巖裂隙的時空演化規律，構建了形式各異、功能多樣的數學模型，豐富了裂隙/裂隙場的識別理論。但針對相似材料模擬實驗在開挖過程中形成的裂隙場，開展自動識別研究的還不多。

基于此，本文以某礦15235工作面為工程研究背景，開展物理模擬實驗，反演煤層開挖過程中上覆巖層的運移過程，就采動過程中發育的裂隙場開展識別和定量分析研究。

1 Canny算法裂隙識別的不足及其改進

1.1 基于Canny算法的裂隙識別

Canny邊緣檢測算法是一種數字圖像裂隙識別方法，它是基于變分原理導出的一種用高斯模板導數逼近的最優算子。裂隙識別包括4個步驟：高斯濾波、梯度及梯度方向計算、非極大值抑制和邊緣連接。Canny算法首要步驟就是利用高斯函數對圖像進行濾波，圖像上任意點(x，y)二維高斯函數可表達為：

(1)

式中，σ為高斯濾波器參數，控制圖像濾波的平滑程度。

關于梯度及梯度方向計算，Canny算法利用一階微分來計算圖像f各點處的梯度幅值及梯度方向，進而得到相應的梯度幅值圖像C和梯度方向。任意點(i，j)處沿x和y兩個方向的偏導數Cx(i，j)和Cy(i，j)可用式(2)和式(3)來表達：

Cx(i，j)=[f(i，j+1)-f(i，j)+f(i+1，j+1)-f(i+1，j)]/2

(2)

Cy(i，j)=[f(i，j)-f(i+1，j)+f(i+1，j+1)]

(3)

根據式(2)和式(3)可解算出任意點(i，j)處的梯度及其方向：

(4)

(5)

將當前像素的梯度強度與沿正負梯度方向上的2個像素的梯度強度進行比較，如果前者與后者相比較大，則將當前像素保留為邊緣點，否則該像素被抑制；然后通過閾值來識別邊緣像素點，最后對識別出的邊緣像素點進行連接，從而達到對裂隙邊緣的自動提取。

1.2 Canny算法裂隙識別的不足及優化

Canny算法在裂隙識別方面具有一定的優勢，但需要人為設定閾值，一方面費時費力，導致計算資源浪費，另一方面不可避免地影響裂隙識別的準確度。為了克服這種缺陷，論文提出一種優化算法來自動選取最優閾值，實現圖像前景信息(裂隙)與背景的分割。

現假設一幅圖像的像素大小為M×N，其灰度級為L，圖像上某點(x，y)處的灰度值為f(x，y)，由于f(x，y)∈[0，L-1]，如果圖像中灰度級 包含的像素個數為fi，則該級灰度出現的概率P(i)可表達為式(6)：

(6)

將圖像中的像素按灰度級用閾值t(灰度值)劃分為2類，即背景C0和目標C1，背景C0的灰度級為0～t-1，目標C1的灰度級為t～L-1，則背景C0和目標C1對應像素的灰度值滿足條件：{f(x，y)

背景C0和目標C1出現的概率ω0和ω1可分別表達為式(7)和式(8)：

(7)

(8)

背景C0和目標C1的平均灰度值分別為：

(9)

(10)

圖像的總平均灰度值為：

(11)

為了確定圖像分割的最優閾值，平衡背景方差和目標方差值，構建了顧及權重系數的類間方差函數δ2(t)：

δ2(t)=ω0[μ-μ0(t)]2+aω1[μ-μ1(t)]2

(12)

式中，a為權重系數，表征圖像目標像素個數與整幅圖像像素個數之比。

將式(9)、式(10)和式(11)代入式(12)即可計算不同t值下的類間方差，當δ2(t)最大時的t值即為所求最優閾值，即:

(13)

2 工程驗證

2.1 相似材料模型鋪設

15235工作面開采2號煤層，平均開采深度約780 m，工作面走向長度870 m，傾向長度140 m，緩傾斜煤層，開采厚度5.5 m，走向長壁后退式放頂煤開采。實驗采用CM400/22平面應力試驗臺，模型架尺寸4 000 mm×220 mm×1 500 mm，模型相似比為1/100，模型重度之比為0.6，時間相似準數為24。相似模型中骨料為河砂、石膏、膩子粉，緩凝劑為硼砂，模擬煤層上方一部分巖層，其余部分加等效載荷。模型鋪設高度133.5 cm，上覆巖層厚度120 cm，下伏地層厚度8 cm。根據工作面2-2鉆孔柱狀圖，相似材料模型共鋪設41層。為了能采用數字散斑法獲取模型表面任意點的形變量，需要制作人工紋理(膩子襯底)，人工紋理用毛筆蘸墨汁點涂完成，以增強形變的量測精度。

2.2 裂隙識別過程

下面以峰峰礦區15235工作面為工程背景，分別對相似材料實驗中的第16次和第29次開挖(工作面推進距離分別為118、255 cm)覆巖應力穩定后模型表面出現的裂隙進行識別。首先采用Agisoft Mateshape圖像處理平臺，對現場采集的影像進行拼接、校正，生成正射影像(圖1)。由圖1可知，擬處理區域環境較為復雜，模型表面光線偏暗，而且模型架、護板、網格線、編碼點和非編碼點等都可能影響裂隙的識別效果。因此，首先需要對圖2進行預處理(剔除模型架、護板等)，消除潛在的外部影響因素，再進行后續的圖像處理。

圖1 工作面不同推進度時影像對比度增強處理Fig.1 Image contrast enhancement processing under different advancing degrees of working face

為了對比Canny算法優化前后的閾值分割效果，作者基于Matlab 研發平臺，編制了計算程序，對比分析算法優化前后兩幅影像數據的閾值分割效果以及模型表觀裂隙識別結果，驗證算法的有效性。由于實驗獲取的模型圖像整體較暗，識別目標與背景的灰度值接近，因此需要對選取的各幅影像進行對比度增強預處理，處理后的結果如圖2所示。基于式(13)，對圖2(a)、圖2(b)進行閾值分割，分割結果如圖2所示。根據式(13)，確定的兩幅影像權重系數分別為5.73%、7.68%。對比優化前后各幅影像可知，算法優化前，各影像均不同程度地存在裂隙邊緣噪聲、裂隙輪廓模糊不清，嚴重影像裂隙識別的精準度。算法優化后，各幅影像表觀裂隙紋理清晰，對比度明顯增強，閾值分割效果顯著提升。

圖2 算法優化前后的閾值分割結果Fig.2 Threshold segmentation results before and after algorithm optimization

完成圖像閾值分割之后，即可進行相似材料表面裂隙的邊緣識別與連接，工作面推進至118 cm時算法優化前后的裂隙識別結果如圖3所示，工作面推進至255 cm時算法優化前后的裂隙識別如圖4所示。

圖3 工作面推進至118 cm時算法優化前后的裂隙識別結果Fig.3 Fracture identification results before and after algorithm optimization when the working face advances to 118 cm

對比圖4(a)和圖4(b)可以看出，整個裂隙區劃分為Ⅰ、Ⅱ、Ⅲ三個部分，使用Canny常規算法對模型全景圖像進行處理后，第Ⅰ部分識別的裂隙模糊不清，部分裂隙未能有效識別，第Ⅱ部分主要裂隙未能識別，第Ⅲ部分裂隙模糊不清。算法優化后，裂隙的整體識別率較高，裂隙邊緣輪廓較為清晰。此外，從圖4(b)可以看出，覆巖裂隙主要集中在工作面兩端(開采邊界的內側)，根據礦山開采沉陷學“上三帶”劃分方法，當工作面推進度為255 cm時，上覆巖層中垮落帶高度為13.31 cm，斷裂帶高度為63.36 cm(由于只模擬部分巖層，彎曲帶高度無法確定)，左右垮落角分別為52°和55°，裂隙場整體呈“梯形”。Canny優化算法實現了相似材料模型表面主要裂隙和細微裂隙的準確識別，裂隙損失小，裂隙場輪廓較為清晰，證實了算法的有效性和實用性。

圖4 工作面推進至255 cm時算法優化前后的裂隙識別結果Fig.4 Fracture identification results before and after algorithm optimization when the working face advances to 255 cm

2.3 裂隙定量分析

采動覆巖裂隙能夠響應裂紋的萌生、擴展和動態運移過程，反饋上覆巖層的受損狀態。開展采動覆巖裂隙定量提取和分析研究對于井下煤層氣抽放、水體(含水層)下開采具有一定的實用價值和意義。在裂隙識別的基礎上，作者擬采用像素法和計盒維數法進一步從微觀角度對裂隙進行定量提取、分析和相互驗證。

(1)基于像素法裂隙定量分析。像素法是對采動覆巖裂隙進行定量分析的常用方法之一，評價指標采用裂隙率。裂隙率是擬識別圖像裂隙區域面積與識別區域總面積之比，面積用區域內像素個數之和表示。采用像素法可分別對工作面不同推進度的采動覆巖表觀裂隙進行定量分析，結果見表1。由表1可知，算法優化后較優化前上覆巖層的裂隙率和裂隙增長率均有不同程度的提高，且裂隙率和裂隙增長率與工作面開采尺寸正相關，開采尺寸越大裂隙率越大。當工作面推進至255 cm時算法優化前后的裂隙率分別為5.32%、7.41%，裂隙率增大了21.53%。這是因為算法優化前部分裂隙模糊不清，類間方差函數的引入使原來未能識別的裂隙得以識別，裂隙識別率得到了不同程度的提升。

表1 基于像素法相似材料模型采動覆巖表觀裂隙定量分析結果Tab.1 Quantitative analysis results of apparent fractures of mining overburden based on pixel method similar material mod

(2)基于分形維數的裂隙定量分析。分形維數是一種分形的定量表征，能定量化描述自然界不規則、雜亂無章的復雜現象。原理是利用計算機存儲技術基于點覆蓋法，計算圖像各個狀態的分形維數，從而反映圖像裂隙的微觀狀態。本部分擬采用計盒維數法，基于Matlab平臺編制計算程序，實現工作面不同推進度時裂隙分形維數的自動計算和分析。計盒維數法首先將圖像進行二值化處理，再把圖像分割成邊長為s的小正方形，然后將包含裂隙的小正方形記作1，其余的記作0。若包含裂隙的小正方形數量用N(s)表示，則改變s的大小，N(s)的大小也隨之變化。實際上分形維數計算就是繪制lg(s)-lgN(s)之間的關系曲線，而擬合直線斜率的絕對值即為分形維數。計算公式見式(14)：

(14)

式中，D為裂隙的分形維數；s為小正方形的邊長；

N(s)為包含裂隙的小正方形數量。

根據分形維數理論，對工作面不同推進度三幅覆巖裂隙場影像進行了圖像分形維數計算，結果見表2和圖5所示。

表2 工作面不同推進進度裂隙場分形維數計算Tab.2 Calculation of fractal dimension fracture field in different advancing progress of working face

由表2可知，兩幅影像分形維數計算相關系數均較高，擬合誤差較小，具有良好的線性相關性，表明相似材料模型在開挖過程中，上覆巖層裂紋生長具有自相似的分形特性，反映了不同推進狀態下上覆巖層裂隙發育的表觀結構特征。當工作面推進度分別為118 cm和255 cm時，采動覆巖裂隙場的分形維數分別為1.811和1.885，分形維數為1～2。分形理論認為，二維條件下若分形維數位于0～2，則分形維數值符合客觀規律。因此，采用優化Canny算法提取的3個裂隙場均具有較好的分形特征。其中，工作面推進度為255 cm時的裂隙場分形維數最大，說明此狀態下圖像表觀裂隙被等分后，所占小正方形數量最多，即裂隙率最高，這與像素法定量分析的結果一致。

圖5 工作面不同推進進度時計盒數和小正方形邊長之間的擬合關系Fig.5 Fitting between the number of meter boxes and the side length of small square at different progress of working face

3 結論

(1)以15235工作面為工程背景，構建了幾何相似比為1/100的相似材料實驗模型，通過對定焦相機批量獲取的多基線近景攝影測量影像進行拼接和校正，獲得了目標區的正射影像。

(2)針對復雜環境下Canny算法裂隙識別的不足，構建了顧及權重系數的類間方差函數，確定了影像信息分割最優閾值表達式，實現目標影像前景和背景的完整分割以及裂隙信息的完整識別。

(3)基于Matlab研發平臺，采用Canny優化算法，得到兩組圖像類間方差權重系數分別為5.73%、7.68%，影像信息分割最優閾值分別為115、128，裂隙場信息提取完整，左右垮落角分別為52°和55°，覆巖垮落帶和斷裂帶高度分別為13.31 cm和63.36 cm，算法優化后裂隙率相對增量明顯提升。

(4)采用像素法和分形理論分別對工作面推進度為118 cm和255 cm時的覆巖裂隙場進行了定量分析，上覆巖層裂紋生長具有自相似的分形特性，反映了不同推進狀態下上覆巖層裂隙發育的表觀結構特征。