王 艷 李 晗 陳佳悅 陳衛(wèi)衛(wèi) 王夢菲 閆 迪 李宏霞 王衛(wèi)星

(1.西安航空職業(yè)技術(shù)學院人工智能學院,陜西 西安 710089;2.遼寧省交通高等?？茖W校汽車工程系,遼寧 沈陽 110122;3.倫敦大學學院計算機科學系,倫敦 WC1E 6BT;4.西安航空職業(yè)技術(shù)學院自動化工程學院,陜西 西安 710089;5.長安大學信息工程學院,陜西 西安 710064)

在各種巖石工程應(yīng)用中,巖石節(jié)理裂隙會顯著影響巖體的物理性質(zhì)[1-3]。地質(zhì)科學中已經(jīng)開發(fā)并采用了許多測量方法。同時,巖石裂隙對礦山工程中的巖石爆破質(zhì)量的影響也很大。同時,對于公路、鐵路、隧道和大壩,精確的巖石節(jié)理裂隙數(shù)據(jù)對于評估巖石邊坡的穩(wěn)定性至關(guān)重要。與此相關(guān)的研究有很多,如在混凝土和水泥工程中,裂縫將主導設(shè)施的破壞過程。為了解決由裂縫引起的工程問題,許多研究人員和工程師研究并開發(fā)了用于測量這些裂縫的各種設(shè)備和數(shù)學算法。

工程地質(zhì)學中有一套完整的理論和方法來測量和描述巖體中各種節(jié)理裂隙的方向、密度和其他特征[1]。早些年,1976 年P(guān)riest 和Hudson 引入了巖石不連續(xù)間距[2];然后,在1981 年,他們測量了一個方向上的裂隙間距[3];之后,進行了更精確的節(jié)理裂隙測量[4]。近年來,在野外測量、巖心鉆探、鉆孔電視和攝像機測量中獲得巖石節(jié)理裂隙的1D 和2D 的測量和分析可以用圖像技術(shù)進行。一個典型的例子:巖石質(zhì)量指標(RQD)被定義為由≥0.1 m 的全長組成的鉆孔巖芯或掃描線的百分比,并已被廣泛用作巖體巖石節(jié)理裂隙發(fā)育程度的工程測量[5-6]。

通過從巖石鉆孔中獲取三維RQD 信息,研究人員使用蒙特卡洛模擬方法創(chuàng)建了一個假設(shè)的三維裂縫網(wǎng)絡(luò)。該方法允許對所獲得的數(shù)據(jù)進行測量和分析。以驗證該方法,將計算的RQD值與實際RQD值進行比較,得出了令人滿意的擬合[5],這可能與巖石表面而非鉆孔的3D 裂縫特征有關(guān)。建立了幾個模型中的網(wǎng)絡(luò)可信渲染,用于節(jié)理裂隙網(wǎng)絡(luò)表征。然后,應(yīng)用最佳擬合Ferret 方法自動確定節(jié)理裂隙面積,然后計算裂隙長寬、裂隙線端點及交叉點、顆粒分布等參數(shù)來綜合評價節(jié)理裂隙的分布和發(fā)育狀況。最終得到網(wǎng)絡(luò)的孔隙率、密度、連接性和復雜性。然而,了解巖石裂隙的三維分布要復雜得多。因此,大多數(shù)研究人員都進行了三維建模和模擬,以避免直接在巖體中進行測量。

自20 世紀80 年代以來,圖像分析技術(shù)得到了迅速發(fā)展,并研究了許多用于巖石節(jié)理裂隙測量和分析的算法[6]。在早期對巖石節(jié)理裂隙的圖像處理和分析的研究中,主要工作是對單條裂隙的檢測和識別(如孔徑檢測)[7]。近年來,研究已擴展到多裂隙處理和分析,包括巖石節(jié)理裂隙方位測量、空間分析等[8]。以前的研究是或已經(jīng)是基于與不同巖石工程相關(guān)的巖石特征而進行的。通常,在巖石節(jié)理裂隙圖像處理和分析之前,圖像是通過不同的傳感器而獲取的,所以圖像的類型和設(shè)備有關(guān)。然后根據(jù)圖像的類型和特征,利用計算機軟件對圖像進行預處理、分割和裂縫測量與分析。最后,將裂隙分析數(shù)據(jù)應(yīng)用于估算應(yīng)力、應(yīng)變、位移、滲流計算等。在過去的10 年里,對于3D 巖石節(jié)理裂隙的測量,激光掃描技術(shù)已在現(xiàn)場使用,主要用于整個巖石表面或單個大裂隙的測量和分析[9-10]。而在實驗室中,對于小型的3D 節(jié)理裂隙巖石樣本,可以更詳細地對節(jié)理裂隙參數(shù)進行測量和分析,一般都應(yīng)用了CT 技術(shù)[10]。無論應(yīng)用何種硬件,都可以使用1D 和2D 裂縫信息來估計3D 巖石節(jié)理裂隙情況。通常,在涉及單個或多個巖石節(jié)理裂隙的圖像處理中,首先將圖像增強算法應(yīng)用于巖石節(jié)理裂隙圖像[11-12],然后根據(jù)對象(節(jié)理裂隙)的特征研究提取算法[13-14]。在圖像分析中,巖石節(jié)理裂隙的標準參數(shù)包括2D 孔隙度、孔徑、節(jié)理裂隙密度和方向等。

近5 年來,語義分割在路面裂縫檢測中已經(jīng)成為了一個研究熱點。其分割網(wǎng)絡(luò)是基于深卷積網(wǎng)絡(luò)和神經(jīng)網(wǎng)絡(luò)的模型,語義分割最基本的任務(wù)是對路面圖像中的不同類別圖像像素點進行歸類,分類的目的以區(qū)分圖像中不同的目標對象(如:裂縫和非裂縫)[15-16]。進一步的發(fā)展在 U-net 等模型結(jié)構(gòu)中引進了擴張卷積的運算,可在不降維數(shù)的情況下增大局部感受野信息和收集多尺度信息。在 Deeplabv1 等網(wǎng)絡(luò)中,能夠更進一步地增加展開卷積之比例,并且在相鄰的像素關(guān)系推理中加入條件隨機場模型。當前,Deeplab 系列網(wǎng)絡(luò)仍然在不斷地發(fā)展中,其骨干網(wǎng)絡(luò)也得到了應(yīng)用[17]。但天然的巖石節(jié)理裂隙比人造的混凝土裂縫要復雜得多,所以深度學習等語義分割只能用在特殊場景中,很難有普遍的適應(yīng)性。

綜上所述,巖石節(jié)理裂隙的測量和分析在巖石工程中非常重要,圖像測量手段有一定的優(yōu)越性,但由于節(jié)理裂隙圖像千變?nèi)f化及噪聲太多,導致二維和三維的節(jié)理裂隙測量和分析都會遇到很多問題?；诖?本文研究了一種復雜巖石節(jié)理裂隙邊界線的檢測和跟蹤算法,研究內(nèi)容主要包括三個步驟:① 根據(jù)巖石節(jié)理裂隙圖像的特點,提出一種新的圖像預處理方法來平滑和增強節(jié)理裂隙信息,研究內(nèi)容主要包括圖像縮小,Geodesic Shadow Removal 和基于Hessian 矩陣的圖像增強;② 根據(jù)Live-wire Contour 的算法思想,研究了一種新的圖像分割算法,用于提取節(jié)理裂隙邊界線上的特征點;③ 在點距離和角度差都滿足基本要求的前提下,綜合獲取的特征點方向、亞像素精度定位及強度進行臨近特征點連接。

1 圖像預處理

在通常情況下,巖石表面比較粗糙,噪聲較多,直接進行目標邊緣線的跟蹤比較困難,所以要進行圖像預處理。預處理主要包括3 個方面:圖像縮小,如果圖像尺寸太大,為了減少計算量,要進行圖像縮小,這里要進行的是有選擇地縮小,圖像縮小本著保持邊緣信息及減少噪聲的原則進行;圖像平滑,主要是盡可能多地去除圖像中的各種噪聲,如:黑白點、陰影、三維體等;圖像增強,也就是節(jié)理裂隙邊緣的增強,增強邊緣的原則是盡可能同時抑制噪聲的產(chǎn)生和增強。3個方面的具體敘述如下。

1.1 圖像縮小

為了減少圖像存儲量及后續(xù)算法運行時間,我們需要設(shè)計一種算法,既能縮小圖像尺寸又不丟失節(jié)理裂隙信息,這種算法與多尺度空間理論有關(guān)。我們將在原始圖像中,以每個檢測點為中心,取3×3 模板的3 個最小灰度值的平均灰度值作為縮小圖像對應(yīng)點的值,如圖1 所示,圖1(b)～(e)是圖4(a)4 種不同算法的縮小結(jié)果。

1.2 圖像平滑

進行Geodesic Shadow Removal 平滑[18],該算法主要包括5 個部分:數(shù)學形態(tài)學“閉合”,高斯平滑,分水嶺,直方圖均衡化,對比度增強。

利用形態(tài)學的“閉合”函數(shù)來消除由于裂隙具有與陰影區(qū)域相似的強度而產(chǎn)生的裂隙。在陰影區(qū)域估計中,它可以防止裂隙被計入陰影。如果應(yīng)用高斯平滑濾波器去除噪聲,可進行下面三步驟。

(1)利用分水嶺法將圖像分割成一系列具有不同平均光照的測地線區(qū)域,{Ri|i=1,…,N,…,M},以S={Ri|i=1,2,…,N}區(qū)域為陰影區(qū)域,以B={Ri|i=N+1,N+2,…,M}區(qū)域為背景區(qū)域。這里,N=M。

(2)選擇亮度補償?shù)淖罴褏⒖紖^(qū)域。我們把一個測地區(qū)域設(shè)為Hi=max(Ii)-min(Ii) 和把一個測地區(qū)域的灰度值的跨度設(shè)為H。它是指與半影區(qū)域相對應(yīng)的測地區(qū)域的H值通常大于背景區(qū)域。它通過公式將圖像大致分為陰影區(qū)和背景區(qū),其中HS和HB分別表示低值和高值的類別。我們可以通過公式{HS,HB} =OTSU({Hi|i=1,2,…,M}) 選擇最佳參考區(qū)域其中,Ai=ave(Ri) 和Di=dev(Ri) 分別表示圖像中像素的平均值和標準差。

(3)將像素(i,j)∈S分配為新值I'i,j=qIi,j+e。其中,Ii,j表示像素(i,j)處的灰度值,I'i,j表示亮度補償?shù)慕Y(jié)果,和DB分別表示S和B的標準差,和分別表示S和B的平均灰度值。這種算法可以增強圖像的對比度。

1.3 節(jié)理裂隙增強

Hessian 矩陣I在點p(x,y)定義為

相應(yīng)的節(jié)理裂隙定義如下:

式中,td(p,θ;σ)=gxxcos2θ+gyysin2θ+gxysin2θ.

對于本文提出的一套方向濾波器—前向過濾器tf(p;σ,ψ1) 和后向過濾器tb(p;σ,ψ2) ,其定義如下:

在上述公式中,采用角度參數(shù)ψ1、ψ2檢測相鄰像素之間的裂紋跡象;常數(shù)d是偏移量參數(shù),盡可能設(shè)置為適當?shù)闹?若其值較低,其貢獻就小,將引起圖像分割錯誤,以至于將不是節(jié)理裂隙的像素錯誤地歸為節(jié)理裂隙像素。

圖像通過定向濾波器的響應(yīng)分別為Tf(p;σ,ψ1)=tf(p;σ,ψ1)×I(p)和Tb(p;σ,ψ2)=tb(p;σ,ψ2)×I(p)。增強節(jié)理裂隙的整體響應(yīng)如式(5)所示:

然后在多個方向上搜索最大響應(yīng)作為結(jié)果的輸出,如圖2 所示。

2 巖石節(jié)理裂隙邊界檢測基于Live-wire 算法的改進

Live wire 是在1992 年由Barrett 和Mortensen 利用交互式邊界提取的一種圖像分割方法[19]。其基本思想是利用動態(tài)規(guī)劃來獲取圖像中給定的2 點間最短路徑,然后構(gòu)造合理的代價函數(shù)來檢測目標的邊界(Live-wire contour),但該算法不是全自動的。本文為了增強其自動化程度,首先對掩模(template)及中心線(skeleton)進行圖像的去噪聲處理及平滑處理得到了新掩模(template_pro)及新的中心線(skeleton_pro);然后把Live-wire 的代價函數(shù)寫成:

式中,fZ(q) 是Laplace 的過零點;fG(q) 是梯度的大小值;fD(p,q) 是梯度的方向;ωZ、ωG、ωD是相應(yīng)權(quán)重。

式中,IL(q)≠0 是圖像及通過LOG 邊界檢測算子的卷積后在圖像中的q點處的新值。

式中,IG(q)是原始圖像與梯度算子卷積運算后在點q處的值。

定義D(p)= (Gy(p),-Gx(p))為點p處的梯度方向。

其中,dp(p,q)=D(p)·L(p,q),dq(p,q)=L(p,q)·D(q)。

L(p,q)定義為

因為Canny 邊界檢測器可以在降低噪聲和精準定位目標邊界點之間選擇一種平衡的方案,故本文的代價函數(shù)是基于Canny 檢測器而構(gòu)造的,從而可以降低噪聲的干擾和提高定位的精度。

節(jié)理裂隙邊界種子點搜索原理如圖3 所示,3 點S、P、E為1 條中心線上的種子點;而T、Q、F和L、R、G分別為節(jié)理裂隙邊界種子點左側(cè)和右側(cè)的對應(yīng)種子點。

圖3 邊界種子點的搜索原理Fig.3 Principle of seed point searching on object boundary

在圖3 中可將邊界線的權(quán)重定義為

根據(jù)上述,那么對一幅無方向的加權(quán)G(V,E)連通圖,可定義其最短路徑為G(V,E),設(shè)其最短路徑圖為Gpath(Vpath,Epath)。首先從G(V,E)中,提取出源點u0,然后放進Vpath,這樣,我們有Vpath= {u0},Epath=?;然后從任一的u∈Vpath,v∈V-Vpath路徑當中,我們選擇權(quán)重weight(u,v) 為最小,把頂點記為weight(u,v)為最小的,其頂點記為u*,v*,而把邊記為＜u*,v*＞,Vpath←Vpath∪{v*},Epath←Epath∪ {＜u*,v*＞},修正源點u0到V-Vpath各個頂點的權(quán)值weight(u0,vi),修正方法如下:

最終,一直重復上述過程,當Vpath=V時,停止。這樣通過Epath中的那條路徑就被確定為最短的路徑。具體步驟及檢測階段效果可以圖4 為例。

圖4 節(jié)理裂隙的邊界檢測Fig.4 Fracture edge tracing

3 實驗與分析

上面詳細地介紹了基于改進的Live-wire Contour算法的邊緣檢測算法及圖像預處理算法[20]。為了表述清楚,這里選擇一幅典型的有明顯裂隙寬度的圖像進行測試,為了體現(xiàn)新算法的優(yōu)越性,我們選擇了常用的11 種算法進行對比試驗。如圖5 所示,原始的巖石節(jié)理裂隙圖像是一幅光照分布不均的圖像,根據(jù)試驗數(shù)據(jù)統(tǒng)計,裂隙邊緣提取過程中用到的參數(shù)設(shè)置為:根據(jù)同類100 幅圖像的訓練,得出Frangi 濾波器對圖像增強的過程中估計的σmax=σmin= 3,二值化過程中選取的T=16,節(jié)理裂隙篩選過程中所選取的參數(shù)S=10,r=3。

圖5 節(jié)理裂隙檢測算法的比較Fig.5 Comparison of joint crack detection algorithm

圖5(b) ～圖5(d)為3 種不同的閾值檢測結(jié)果[7],由于光照不均,大津閾值使得近50%的面積為背景(黑色,包括裂隙),但實際上,裂縫的面積很少,估計約有10%;自動4 個閾值分割結(jié)果要好一些,如果考慮黑色為裂隙,大部分裂隙被檢測到,但有些裂隙沒有完全被檢測出來,而非裂隙區(qū)域也有黑色的像素;動態(tài)閾值分割得要比前面的閾值分割稍微好一些,但也存在著類似的問題,并且產(chǎn)生了許多小目標,為后續(xù)的圖像處理造成困難。圖5(e)～圖5(h)是4種邊界掃描結(jié)果[21]。一階導數(shù)的Sobel 和二階導數(shù)的Laplacian 的邊界檢測結(jié)果是梯度圖像,圖像中的假邊緣點太多,而真的邊緣點又沒有連成線,這會給后續(xù)的二值閾值處理造成很大的難度;而Canny-1 是低閾值的邊界檢測,其結(jié)果給出了太多的噪聲邊界;當設(shè)置高閾值時,Canny-2 產(chǎn)生了欠分割的問題。圖5(i)～圖5(l)是4 種較為復雜和精確的圖像分割算法結(jié)果?；趫D論的MST (最小生成樹)[22]顯然有著欠分割的問題,而聚類分析給出的結(jié)果要好一些,但仍有欠分割的問題,在FCM(模糊聚類)[7]的結(jié)果圖中,棕黃色的區(qū)域包括了大部分裂隙,但有過分割和欠分割問題,合并與分離(Merging & Splitting)基本上能夠檢測出大部分裂隙的區(qū)域,但仍有少量的欠分割現(xiàn)象,還有就是設(shè)置參數(shù)時,要做多次試驗。

總之對于這種噪聲較多且光照不均的巖石節(jié)理裂隙圖像,無論是傳統(tǒng)的圖像分割算法(如基于相似性的閾值法及基于不連續(xù)性的邊界掃描法) 還是更復雜的圖像分割算法 (如基于圖論及基于聚類分析的算法) 都很難直接湊效。

如果用本文研究的預處理和邊緣檢測算法對同樣的圖像進行處理,則可以得到比較好的結(jié)果。當對圖5(a)的圖像進行預處理后,可以得到圖6(a)中的圖像,光照度變得均勻,噪聲會減少。在此基礎(chǔ)上進行裂隙邊緣的特征點檢測及連接,得到初步的裂隙邊緣檢測結(jié)果,如圖6(b)所示,最后進行裂隙擴展和過濾,得到裂隙充填的圖像,如圖6(c)所示。

圖6 新算法對圖5(a)的檢測過程和結(jié)果Fig.6 Fracture image (Fig.5(a)) detection procedure and results by new algorithm

4 結(jié) 論

本文研究了一種節(jié)理裂隙邊緣線跟蹤的新算法。在跟蹤之前,若圖像太大,進行有選擇的縮小,圖像縮小本著保持邊緣信息及減少噪聲的原則進行。然后基于改進Geodesic Shadow Removal 和Hessian 矩陣等算法將具有粗糙面的巖石節(jié)理裂隙圖像進行平滑和增強。根據(jù)圖像預處理的結(jié)果進行節(jié)理裂隙邊緣跟蹤。跟蹤分成2 個步驟:① 基于Live-wire 的思想檢測出節(jié)理裂隙邊緣上的特征點;② 根據(jù)特征點之間的距離及夾角來連接成長短不一的節(jié)理裂隙線段,然后再連接線段以形成完整的節(jié)理裂隙邊緣線。為了驗證該算法的優(yōu)缺點,選擇了200 幅圖像進行試驗,均得到了較為滿意的結(jié)果。為了分析該算法的有效性,選擇了11 種傳統(tǒng)的算法和更精細和復雜的算法對200 幅圖像進行處理分析。結(jié)果表明:對于有明顯節(jié)理裂隙寬度且具有較多噪聲的圖像,新算法有著較強的適應(yīng)性,能夠準確地提取節(jié)理裂隙邊緣線。今后進一步的研究是選取更多特性不同的圖像進行試驗,以改進算法來擴展其適應(yīng)范圍。