王順陽，賀易航，韓 博

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

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

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

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

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

1.1 基于Canny算法的裂隙識(shí)別

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

(1)

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

關(guān)于梯度及梯度方向計(jì)算，Canny算法利用一階微分來計(jì)算圖像f各點(diǎn)處的梯度幅值及梯度方向，進(jìn)而得到相應(yīng)的梯度幅值圖像C和梯度方向。任意點(diǎn)(i，j)處沿x和y兩個(gè)方向的偏導(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)可解算出任意點(diǎn)(i，j)處的梯度及其方向：

(4)

(5)

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

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

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

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

(6)

將圖像中的像素按灰度級(jí)用閾值t(灰度值)劃分為2類，即背景C0和目標(biāo)C1，背景C0的灰度級(jí)為0～t-1，目標(biāo)C1的灰度級(jí)為t～L-1，則背景C0和目標(biāo)C1對(duì)應(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)像素個(gè)數(shù)與整幅圖像像素個(gè)數(shù)之比。

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

(13)

2 工程驗(yàn)證

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

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

2.2 裂隙識(shí)別過程

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

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

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

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

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

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

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

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

2.3 裂隙定量分析

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

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

表1 基于像素法相似材料模型采動(dòng)覆巖表觀裂隙定量分析結(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ì)算機(jī)存儲(chǔ)技術(shù)基于點(diǎn)覆蓋法，計(jì)算圖像各個(gè)狀態(tài)的分形維數(shù)，從而反映圖像裂隙的微觀狀態(tài)。本部分?jǐn)M采用計(jì)盒維數(shù)法，基于Matlab平臺(tái)編制計(jì)算程序，實(shí)現(xiàn)工作面不同推進(jìn)度時(shí)裂隙分形維數(shù)的自動(dòng)計(jì)算和分析。計(jì)盒維數(shù)法首先將圖像進(jìn)行二值化處理，再把圖像分割成邊長為s的小正方形，然后將包含裂隙的小正方形記作1，其余的記作0。若包含裂隙的小正方形數(shù)量用N(s)表示，則改變s的大小，N(s)的大小也隨之變化。實(shí)際上分形維數(shù)計(jì)算就是繪制lg(s)-lgN(s)之間的關(guān)系曲線，而擬合直線斜率的絕對(duì)值即為分形維數(shù)。計(jì)算公式見式(14)：

(14)

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

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

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

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

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

圖5 工作面不同推進(jìn)進(jìn)度時(shí)計(jì)盒數(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的相似材料實(shí)驗(yàn)?zāi)Ｐ?，通過對(duì)定焦相機(jī)批量獲取的多基線近景攝影測(cè)量影像進(jìn)行拼接和校正，獲得了目標(biāo)區(qū)的正射影像。

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

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

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