基于Canny優(yōu)化算法物理模擬實(shí)驗(yàn)覆巖裂隙場(chǎng)自動(dòng)識(shí)別研究

王順陽(yáng)，賀易航，韓 博

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

大規(guī)模煤炭資源開(kāi)采破壞了圍巖的應(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)我國(guó)煤炭工業(yè)健康穩(wěn)步發(fā)展具有現(xiàn)實(shí)意義。

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

基于此，本文以某礦15235工作面為工程研究背景，開(kāi)展物理模擬實(shí)驗(yàn)，反演煤層開(kāi)挖過(guò)程中上覆巖層的運(yùn)移過(guò)程，就采動(dòng)過(guò)程中發(fā)育的裂隙場(chǎng)開(kāi)展識(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算法利用一階微分來(lái)計(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)來(lái)表達(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)，否則該像素被抑制；然后通過(guò)閾值來(lái)識(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)化算法來(lái)自動(dòng)選取最優(yōu)閾值，實(shí)現(xiàn)圖像前景信息(裂隙)與背景的分割。

現(xiàn)假設(shè)一幅圖像的像素大小為M×N，其灰度級(jí)為L(zhǎng)，圖像上某點(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類(lèi)，即背景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ù)的類(lèi)間方差函數(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值下的類(lèi)間方差，當(dāng)δ2(t)最大時(shí)的t值即為所求最優(yōu)閾值，即:

(13)

2 工程驗(yàn)證

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

15235工作面開(kāi)采2號(hào)煤層，平均開(kāi)采深度約780 m，工作面走向長(zhǎng)度870 m，傾向長(zhǎng)度140 m，緩傾斜煤層，開(kāi)采厚度5.5 m，走向長(zhǎng)壁后退式放頂煤開(kāi)采。實(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í)別過(guò)程

下面以峰峰礦區(qū)15235工作面為工程背景，分別對(duì)相似材料實(shí)驗(yàn)中的第16次和第29次開(kāi)挖(工作面推進(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)可以看出，覆巖裂隙主要集中在工作面兩端(開(kāi)采邊界的內(nèi)側(cè))，根據(jù)礦山開(kāi)采沉陷學(xué)“上三帶”劃分方法，當(dāng)工作面推進(jìn)度為255 cm時(shí)，上覆巖層中垮落帶高度為13.31 cm，斷裂帶高度為63.36 cm(由于只模擬部分巖層，彎曲帶高度無(wú)法確定)，左右垮落角分別為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)移過(guò)程，反饋上覆巖層的受損狀態(tài)。開(kāi)展采動(dòng)覆巖裂隙定量提取和分析研究對(duì)于井下煤層氣抽放、水體(含水層)下開(kāi)采具有一定的實(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é)果見(jiàn)表1。由表1可知，算法優(yōu)化后較優(yōu)化前上覆巖層的裂隙率和裂隙增長(zhǎng)率均有不同程度的提高，且裂隙率和裂隙增長(zhǎng)率與工作面開(kāi)采尺寸正相關(guān)，開(kāi)采尺寸越大裂隙率越大。當(dāng)工作面推進(jìn)至255 cm時(shí)算法優(yōu)化前后的裂隙率分別為5.32%、7.41%，裂隙率增大了21.53%。這是因?yàn)樗惴▋?yōu)化前部分裂隙模糊不清，類(lèi)間方差函數(shù)的引入使原來(lái)未能識(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ī)則、雜亂無(wú)章的復(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)行二值化處理，再把圖像分割成邊長(zhǎng)為s的小正方形，然后將包含裂隙的小正方形記作1，其余的記作0。若包含裂隙的小正方形數(shù)量用N(s)表示，則改變s的大小，N(s)的大小也隨之變化。實(shí)際上分形維數(shù)計(jì)算就是繪制lg(s)-lgN(s)之間的關(guān)系曲線，而擬合直線斜率的絕對(duì)值即為分形維數(shù)。計(jì)算公式見(jiàn)式(14)：

(14)

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

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

根據(jù)分形維數(shù)理論，對(duì)工作面不同推進(jìn)度三幅覆巖裂隙場(chǎng)影像進(jìn)行了圖像分形維數(shù)計(jì)算，結(jié)果見(jiàn)表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)性，表明相似材料模型在開(kāi)挖過(guò)程中，上覆巖層裂紋生長(zhǎng)具有自相似的分形特性，反映了不同推進(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ù)最大，說(shuō)明此狀態(tài)下圖像表觀裂隙被等分后，所占小正方形數(shù)量最多，即裂隙率最高，這與像素法定量分析的結(jié)果一致。

圖5 工作面不同推進(jìn)進(jìn)度時(shí)計(jì)盒數(shù)和小正方形邊長(zhǎng)之間的擬合關(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)Ｐ?，通過(guò)對(duì)定焦相機(jī)批量獲取的多基線近景攝影測(cè)量影像進(jìn)行拼接和校正，獲得了目標(biāo)區(qū)的正射影像。

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

(3)基于Matlab研發(fā)平臺(tái)，采用Canny優(yōu)化算法，得到兩組圖像類(lèi)間方差權(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)行了定量分析，上覆巖層裂紋生長(zhǎng)具有自相似的分形特性，反映了不同推進(jìn)狀態(tài)下上覆巖層裂隙發(fā)育的表觀結(jié)構(gòu)特征。