基于Beamlet變換的數(shù)字地質(zhì)露頭裂縫線性特征提取

付含聰 鄧 帆* 曾齊紅 邵燕林

1(長江大學(xué)地球科學(xué)學(xué)院 湖北 武漢 430100) 2(中國石油勘探開發(fā)研究院 北京 100083)

0 引 言

碳酸鹽巖油氣藏作為世界上最重要的油氣藏類型之一，廣泛分布在各種陸相盆地。油氣勘探開發(fā)的成功很大程度上取決于根據(jù)巖性充分表征裂縫體系，縫洞識(shí)別與表征是進(jìn)行碳酸鹽巖縫洞系統(tǒng)的研究和儲(chǔ)層預(yù)測的基礎(chǔ)工作[1]。但由于地質(zhì)現(xiàn)象及裂縫成因的復(fù)雜性，包括巖石巖性、巖心粗糙程度、光照陰影噪聲的干擾，很難精確識(shí)別和預(yù)測各種尺度的裂縫，對(duì)數(shù)字地質(zhì)露頭裂縫的線性特征自動(dòng)提取一直是個(gè)難題。目前，基于圖像處理技術(shù)的露頭裂縫提取方法主要有基于微分算子和基于變換兩大類。前者是利用邊緣檢測算子，如Canny算子[2]、Sobel算子[3]、Roberts算子、Prewitt算子[4]等直接對(duì)圖像進(jìn)行邊緣檢測，計(jì)算簡單、運(yùn)算速度快，但對(duì)噪聲比較敏感，特別是在信噪比很低的情況下，容易丟失變化緩慢的邊緣；后者如Radon變換、Hough變換[5]、小波變換等,雖然對(duì)噪聲不敏感，但存在無法確定線段長度及起始點(diǎn)信息、定位不準(zhǔn)確等問題，且對(duì)曲線的檢測較復(fù)雜。

針對(duì)以上提取問題，基于圖像多尺度分析的線特征提取有了很大的發(fā)展。其不僅可以獲得圖像的主要特征，還能以不同的細(xì)節(jié)程度構(gòu)造圖像。Beamlet變換理論是由Donoho等[6]和Huo等[7]提出的一種多尺度圖像分析框架，以不同尺度、不同方向的小線段作為基，最初用于從圖像中提取直線和曲線，由于其抗噪性能強(qiáng)的優(yōu)勢(shì)，許多學(xué)者對(duì)算法進(jìn)行了改進(jìn)。如陳雨等[8]對(duì)基于Beamlet變換的線特征提取算法進(jìn)行改進(jìn)，提出了一種圖像邊緣檢測的新算法；肖進(jìn)勝等[9]結(jié)合Beamlet與K-Means算法有效地解決了車道線識(shí)別等問題；曾接賢等[10]針對(duì)傳統(tǒng)Beamlet無結(jié)構(gòu)算法提取圖像線性特征時(shí)存在重疊模糊的缺陷,提出了改進(jìn)的Beamlet無結(jié)構(gòu)算法與Canny算子相結(jié)合的方法，有效提高了圖像線特征提取的準(zhǔn)確性和連續(xù)性。Beamlet變換算法已經(jīng)在許多重要領(lǐng)域得到了廣泛的應(yīng)用。梅小明等[11]利用Beamlet變換提取了遙感圖像線性特征；Ying等[12]通過Beamlet算法，利用數(shù)字圖像對(duì)道路表面裂縫進(jìn)行了提取和分類；荊智輝等[13]利用Beamlet變換對(duì)地震小斷層和裂縫發(fā)育帶進(jìn)行了提??；Sahli等[14]運(yùn)用Beamlet變換對(duì)機(jī)場跑道進(jìn)行了檢測。這些算法都取得了不錯(cuò)的效果，很大程度上解決了背景噪聲復(fù)雜、線條豐富的圖像線特征提取問題。

本文針對(duì)碳酸鹽巖油氣藏地質(zhì)露頭裂縫難以提取的特點(diǎn)，將Beamlet變換算法引入到數(shù)字露頭裂縫自動(dòng)提取中，設(shè)計(jì)并實(shí)現(xiàn)了一種數(shù)字露頭裂縫自動(dòng)提取算法。首先利用多尺度自適應(yīng)增強(qiáng)算法減弱背景噪聲和光照等影響；然后基于Beamlet算法提取裂縫；最后針對(duì)出現(xiàn)的誤檢、斷裂問題設(shè)計(jì)了離散Beamlet基連接算法，最終實(shí)現(xiàn)裂縫的精確提取。

1 算法設(shè)計(jì)

露頭裂縫提取方法包含圖像增強(qiáng)處理、Beamlet裂縫提取、離散裂縫線性特征連接三個(gè)部分(如圖1所示)。野外采集到的數(shù)字圖像由于巖性、光照、陰影等各種因素的影響，會(huì)使圖像包含復(fù)雜的背景噪聲，為了提高裂縫識(shí)別精度、可靠性、實(shí)用性，在進(jìn)行Beamlet變換之前有必要對(duì)圖像進(jìn)行前期處理。

1.1 RGB轉(zhuǎn)灰度圖像

野外采集的露頭裂縫圖像通常是多通道影像，為了提高運(yùn)行速度，降低算法復(fù)雜度，將多通道圖像轉(zhuǎn)化為單通道圖像[15]。本文利用加權(quán)平均法對(duì)露頭裂縫圖像進(jìn)行灰度轉(zhuǎn)換，計(jì)算公式為：

Gray=0.229×R+0.587×G+0.114×B

(1)

式中：Gray為灰度圖像；R為紅色通道；G為綠色通道；B為藍(lán)色通道。

1.2 多尺度自適應(yīng)增強(qiáng)算法

采用多尺度自適應(yīng)增強(qiáng)算法對(duì)灰度圖像進(jìn)行增強(qiáng)處理。該算法能均勻裂縫光照背景，有效抑制噪聲干擾，且運(yùn)算速度快[16-17]，算法實(shí)現(xiàn)過程主要分為以下五步：

(1) 對(duì)一幅n×n(n=2J)的數(shù)字圖像，如果定義圖像的邊長為單位長度，那么每個(gè)像素可以看成是邊長1/n的方塊，所有像素都位于[0,1]2中。在一定的尺度上將數(shù)字圖像分割成許多小的子正方形，尺度的選擇為2的冪次方，一個(gè)2J×2J的圖片最小可分解到像素級(jí)，最大可分解到2J×2J級(jí)。

(2) 分別計(jì)算每個(gè)子正方形像素灰度平均值Gmean、最小值Gmin、最大值Gmax。

(3) 對(duì)每個(gè)子正方形窗體，設(shè)置一個(gè)閾值上限r(nóng)h和一個(gè)閾值下限r(nóng)l，當(dāng)像素點(diǎn)的灰度值在閾值外時(shí)則被自動(dòng)認(rèn)為是噪聲點(diǎn)，將該噪聲點(diǎn)的灰度值變成該子正方形像素灰度平均值Gmean，當(dāng)像素點(diǎn)灰度值在閾值內(nèi)時(shí)，則不變。閾值[rh,rl]的范圍如下：

rh=Gmean+(Gmax-Gmean)×f

(2)

rl=Gmean-(Gmean-Gmin)×f

(3)

式中：f是一個(gè)限定性因子，其取值根據(jù)不同圖像而定，在本文中f取值為0.5。

I′=I×f′

(4)

式中：I′為增強(qiáng)后的影像；I為原始影像；f′為校正系數(shù)

圖2(b)和圖2(c)為分別采用8×8和16×16的尺度對(duì)灰度圖像圖2(a)進(jìn)行自適應(yīng)增強(qiáng)處理的結(jié)果，結(jié)果表明，在8×8的尺度下，裂縫信息部分丟失，造成圖像失真，但背景均勻；在16×16的尺度下，裂縫比較突出，但背景均勻效果下降。一般裂縫越粗，圖像越大，選擇的尺度就越大。

1.3 大津閾值分割算法

對(duì)數(shù)字圖像進(jìn)行增強(qiáng)處理后，利用大津閾值[18]分割算法分離裂縫和背景噪聲得到二值圖像，具體原理如下：

假設(shè)圖像有n個(gè)灰度級(jí)別，根據(jù)灰度區(qū)間[0,k]和[k+1,n]將圖像分成A、B兩部分，分別計(jì)算A、B區(qū)域的平均灰度mA、mB和全圖的平均灰度m，以及A、B區(qū)域出現(xiàn)的概率PA、PB，則類間方差σ2可表示為：

σ2=PA(mA-m)2+PB(mB-m)2

(5)

通過遍歷得到最大類間方差σ2，即為最佳灰度分割閾值。

2 Beamlet變換理論

多尺度幾何分析工具Beamlet變換主要包括Beamlet字典、Beamlet變換、Beamlet金字塔、Beamlet圖和Beamlet算法五部分。

2.1 Beamlet字典

Beamlet字典是以不同尺度、方向和位置的小線段為基構(gòu)成的集合，對(duì)一幅分辨率為n×n(n=2J)的圖像，所有像素都位于[0,1]2中。在尺度0≤j≤J上，可以將圖像塊[0,1]2分為大小為2-j×2-j的2j×2j個(gè)小方塊。在所有尺度下取固定分辨率σ=2-j-k(k≤0,J=log2n),對(duì)每個(gè)小方塊以其邊界σ為間隔取標(biāo)記點(diǎn)，可以得到Mj=4×2K×2J-j個(gè)標(biāo)記點(diǎn)，Beamlet基就是指任意兩個(gè)標(biāo)記點(diǎn)的連線。所有尺度下的Beamlet基構(gòu)成的集合就稱為Beamlet詞典。在尺度j上，Beamlet基的數(shù)量為6(2K+J-j-1)2+8(2K+J-j-1)+2。圖3為不同尺度的Beamlet基，圖4為用一個(gè)連續(xù)的Beamlet鏈來逼近一條線段，即一條Beam。

2.2 Beamlet變換

假設(shè)f(x1,x2)為[0,1]2上的連續(xù)函數(shù)，則函數(shù)f的連續(xù)Beamlet變換可以理解為線積分的集合。

(6)

式中：Bnp指對(duì)n×n圖像，分辨率為ρ時(shí)，不同尺度下所有Beamlet的集合；f(x,y)為b上圖像函數(shù)；Tf(b)為對(duì)應(yīng)b的Beamlet變換系數(shù)。

對(duì)于實(shí)際的一幅n×n數(shù)字圖像是由n2個(gè)像素組成的離散函數(shù)，若圖像函數(shù)為fi1i2，則fi1i2的Beamlet變換必須首先通過插值得到對(duì)應(yīng)的連續(xù)函數(shù)f(x,y),再利用式(6)得到。連續(xù)函數(shù)f(x,y)的求解公式如下:

(7)

式中：φi1i2(x,y)表示連續(xù)插值函數(shù)。

連續(xù)插值函數(shù)有多種選擇[7]，為提高算法效率，采用平均插值的方法。即fi1i2可以看作是對(duì)連續(xù)函數(shù)f(x,y)上的像素值的平均，則式(7)可改寫為：

fi1i2=Ave{f|Pixel(i1i2)}

(8)

式中：Ave表示平均值函數(shù)；Pixel表示像素點(diǎn)。

在此意義下，一條Beamlet基的Beamlet變換實(shí)際上就是Beamlet基上包含的所有像素的灰度值加權(quán)和，因此式(8)可寫成：

(9)

式中：Wi1i2為Beamlet基b在圖像內(nèi)所包含的線段的長度。

2.3 離散Beamlet基連接

Beamlet變換能有效提取露頭裂縫圖像線性特征，但提取的線段有斷裂，因此需要對(duì)離散的線段進(jìn)行連接，流程如圖5所示。具體步驟如下:

1) 在Beamlet變換選取尺度j的基礎(chǔ)上，首先將含有Beamlet基的小矩陣狀態(tài)定義為“未檢測”，不包含Beamlet基的小矩陣狀態(tài)定義為“空”，計(jì)算所有“未檢測”狀態(tài)小矩陣中Beamlet基的角度。

2) 對(duì)每個(gè)Beamlet基分別計(jì)算在八鄰域內(nèi)相鄰Beamlet基的角度差，若角度差小于設(shè)定的閾值(本文取45°)，則進(jìn)行步驟3)，若不符合，則認(rèn)為不是同一條裂縫，不進(jìn)行連接，直接考慮下一鄰域。

3) 在上一步的基礎(chǔ)上，考慮兩個(gè)Beamlet基的相鄰端點(diǎn)距離是否也在設(shè)定的閾值內(nèi)(本文取10個(gè)像素點(diǎn)的距離)，若符合則連接兩個(gè)Beam基，反之不連接，考慮一條Beam基所有八連通區(qū)域后，將Beam基所在的小矩陣狀態(tài)改為“檢測”。

4) 循環(huán)考慮所有“未檢測”狀態(tài)的小矩形。

3 Beamlet算法實(shí)現(xiàn)具體步驟

1) 首先利用雙線性插值算法將圖像的大小轉(zhuǎn)換為n×n(n=2J)。將圖像灰度化后采用自適應(yīng)增強(qiáng)算法處理，用大津閾值分割算法進(jìn)行分割，得到裂縫二值圖像。

2) 根據(jù)裂縫二值圖像大小選擇要?jiǎng)澐值某叨萰(j=2m,j

3) 按照尺度j將圖像劃分成許多二進(jìn)方塊，根據(jù)二進(jìn)方塊的大小構(gòu)建所有Beamlet基并保存在一個(gè)元組中，即Beamlet字典。用構(gòu)建的Beamlet字典分別對(duì)每個(gè)小方塊進(jìn)行Beamlet變換，計(jì)算所有二進(jìn)方塊的Beamlet變換系數(shù)Tf(b)和相應(yīng)的Beamlet基長度L(b)。

4) 設(shè)定一個(gè)閾值f，當(dāng)方塊的最大Beamlet系數(shù)Tf(b)小于f時(shí)，認(rèn)為該方塊能量不足，將其自動(dòng)識(shí)別為噪聲，并將該方塊置零作為背景。反之當(dāng)方塊的Beamlet系數(shù)大于f時(shí)，選取該方塊最大的Beamlet系數(shù)所對(duì)應(yīng)的Beamlet基長度L(b)作為該子方塊的最優(yōu)基，并顯示出來。這樣既能保證每個(gè)子方塊中只含一條最優(yōu)基，又消除了圖像中非線性特征的點(diǎn)狀、塊狀噪聲。

5) 輸出所有符合條件的Beamlet基，并通過離散Beamlet基連接算法對(duì)斷裂處進(jìn)行連接，最后得到裂縫線性特征的提取結(jié)果圖像。

4 實(shí)驗(yàn)結(jié)果與分析

設(shè)計(jì)了兩組實(shí)驗(yàn)對(duì)算法可行性進(jìn)行了驗(yàn)證，并與Canny邊緣檢測算法和Hough變換算法提取結(jié)果進(jìn)行了比較。第一組為人工生成的帶噪聲模擬圖像，第二組為挑選的局部露頭影像。

4.1 模擬圖像分析

為了驗(yàn)證算法的可行性，本文運(yùn)用Beamlet變換對(duì)自動(dòng)生成的含高斯噪聲的圖像進(jìn)行分析。如圖6所示，圖6(a)為原始圖像，圖像大小為256×256；圖6(b)為高斯噪聲處理(均值為0.4，方差為0.8)的圖像；圖6(c)是采用自適應(yīng)增強(qiáng)算法對(duì)含噪聲圖像進(jìn)行增強(qiáng)后的圖像；圖6(d)是用大津閾值算法分割得到的二值圖像；圖6(e)是利用Beamlet變換對(duì)二值圖像提取線性特征后的結(jié)果。由結(jié)果圖像可以得出，基于像素亮度的圖像處理方式容易受光照和噪音的影響，導(dǎo)致提取結(jié)果出現(xiàn)噪聲和漏檢(圖6(d))，本文方法提取的線性特征準(zhǔn)確，不存在噪聲，對(duì)噪聲不敏感，適用于復(fù)雜背景下線性特征的提取，但存在部分線段斷裂的情況，需要進(jìn)行后處理。

4.2 裂縫檢測

數(shù)字露頭裂縫類型有很多種[19]包括填充縫、非填充縫等。選取這兩類裂縫進(jìn)行了四組實(shí)驗(yàn)。

圖7第一、第二行為填充縫圖像實(shí)驗(yàn)，裂縫為淺色調(diào)，圖上還存在其他淺色區(qū)塊，容易對(duì)裂縫提取造成干擾。圖7第三、第四行為非填充逢，裂縫為深色調(diào)，圖上存在其他深色調(diào)區(qū)塊會(huì)對(duì)裂縫提取造成干擾。同時(shí)圖中還存在明暗不一的情況，會(huì)對(duì)裂縫提取產(chǎn)生負(fù)面影響。圖7(c)-圖7(e)分別為Canny變換、Hough變換和Beamlet變換提取結(jié)果。Canny變換是邊緣提取算法，圖中亮度突變區(qū)域均被很好地檢測出來,但各種非線性特征也被檢測出來，需要大量的后處理。Hough變換對(duì)線性特征方向性有很好的指示，基本不存在非線性特征噪音，但提取的線段不連續(xù)，對(duì)線性特征位置、形態(tài)的表示不精確。Beamlet變換提取的結(jié)果相對(duì)最優(yōu)，準(zhǔn)確地提取出了裂縫的形態(tài)，沒有噪聲，但是存在一些線段中斷的問題，提取結(jié)果需要后續(xù)進(jìn)行連接處理。通過離散Beamlet基連接算法進(jìn)行連接后，與原始圖像進(jìn)行疊加顯示，提取結(jié)果更加精確。

通過使用如下兩種評(píng)價(jià)指標(biāo)對(duì)算法精確度進(jìn)行定性評(píng)價(jià)：

準(zhǔn)確率=正確提取的裂縫長度/真實(shí)裂縫長度

誤檢率=錯(cuò)誤提取的裂縫長度/總的提取裂縫長度

四組測試結(jié)果如表1所示。從精度評(píng)價(jià)結(jié)果可以看出，Canny算法提取準(zhǔn)確度雖然較高，四次實(shí)驗(yàn)準(zhǔn)確率分別為0.84、0.79、0.75、0.72，但存在較高的誤檢率，分別為0.40、0.55、0.49、0.46。Hough算法準(zhǔn)確率較低，四次實(shí)驗(yàn)準(zhǔn)確率分別為0.61、0.39、0.60、0.51，誤檢率相比Canny算法有所降低，分別為0.19、0.22、0.15、0.20。而Beanlet算法四次實(shí)驗(yàn)的準(zhǔn)確率為0.81、0.83、0.89、0.85，誤檢率為0、0.04、0.09、0.13，能在具有較高的準(zhǔn)確率前提下，保證較低的誤檢率。

表1 不同算法提取精度比較

5 結(jié) 語

針對(duì)數(shù)字地質(zhì)露頭裂縫線性特征自動(dòng)提取問題，為了減少人工解譯工作量大、效率低、可靠性差等問題，本文提出一種基于多尺度Beamlet變換算法的精確數(shù)字地質(zhì)露頭裂縫提取方法。在傳統(tǒng)Beamlet變換算法上引入了圖像增強(qiáng)和離散線段連接兩部分。實(shí)驗(yàn)結(jié)果表明該方法能在裂縫表面粗糙度、顏色、光照、陰影等造成的復(fù)雜噪聲下準(zhǔn)確提取裂縫線性特征，能在較高的準(zhǔn)確率前提下，保證較低的誤檢率，具有良好的抗噪能力。與傳統(tǒng)Canny算法和Hough算法相比，提取結(jié)果更加完整、準(zhǔn)確。下一步可考慮采用樣條曲線進(jìn)一步增加提取裂縫的連續(xù)性。