王順陽，賀易航，韓 博

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

大規(guī)模煤炭資源開采破壞了圍巖的應(yīng)力平衡狀態(tài)，導(dǎo)致上覆巖層發(fā)生移動變形進(jìn)而形成裂隙，影響煤礦的安全生產(chǎn)[1-2]。因此，研究采動覆巖裂隙場發(fā)育規(guī)律對于保障礦井安全生產(chǎn)，促進(jìn)我國煤炭工業(yè)健康穩(wěn)步發(fā)展具有現(xiàn)實意義。

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

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

1 Canny算法裂隙識別的不足及其改進(jìn)

1.1 基于Canny算法的裂隙識別

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

(1)

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

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

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)

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

(4)

(5)

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

1.2 Canny算法裂隙識別的不足及優(yōu)化

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

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

(6)

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

背景C0和目標(biāo)C1出現(xiàn)的概率ω0和ω1可分別表達(dá)為式(7)和式(8)：

(7)

(8)

背景C0和目標(biāo)C1的平均灰度值分別為：

(9)

(10)

圖像的總平均灰度值為：

(11)

為了確定圖像分割的最優(yōu)閾值，平衡背景方差和目標(biāo)方差值，構(gòu)建了顧及權(quán)重系數(shù)的類間方差函數(shù)δ2(t)：

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

(12)

式中，a為權(quán)重系數(shù)，表征圖像目標(biāo)像素個數(shù)與整幅圖像像素個數(shù)之比。

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

(13)

2 工程驗證

2.1 相似材料模型鋪設(shè)

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

2.2 裂隙識別過程

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

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

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

圖2 算法優(yōu)化前后的閾值分割結(jié)果Fig.2 Threshold segmentation results before and after algorithm optimization

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

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

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

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

2.3 裂隙定量分析

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

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

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

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

(14)

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

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

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

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

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

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

3 結(jié)論

(1)以15235工作面為工程背景，構(gòu)建了幾何相似比為1/100的相似材料實驗?zāi)Ｐ?，通過對定焦相機(jī)批量獲取的多基線近景攝影測量影像進(jìn)行拼接和校正，獲得了目標(biāo)區(qū)的正射影像。

(2)針對復(fù)雜環(huán)境下Canny算法裂隙識別的不足，構(gòu)建了顧及權(quán)重系數(shù)的類間方差函數(shù)，確定了影像信息分割最優(yōu)閾值表達(dá)式，實現(xiàn)目標(biāo)影像前景和背景的完整分割以及裂隙信息的完整識別。

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

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