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

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

0 引 言

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

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

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

1 算法設(shè)計

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

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

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

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

(1)

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

1.2 多尺度自適應增強算法

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

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

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

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

rh=Gmean+(Gmax-Gmean)×f

(2)

rl=Gmean-(Gmean-Gmin)×f

(3)

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

I′=I×f′

(4)

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

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

1.3 大津閾值分割算法

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

假設(shè)圖像有n個灰度級別，根據(jù)灰度區(qū)間[0,k]和[k+1,n]將圖像分成A、B兩部分，分別計算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)成的集合，對一幅分辨率為n×n(n=2J)的圖像，所有像素都位于[0,1]2中。在尺度0≤j≤J上，可以將圖像塊[0,1]2分為大小為2-j×2-j的2j×2j個小方塊。在所有尺度下取固定分辨率σ=2-j-k(k≤0,J=log2n),對每個小方塊以其邊界σ為間隔取標記點，可以得到Mj=4×2K×2J-j個標記點，Beamlet基就是指任意兩個標記點的連線。所有尺度下的Beamlet基構(gòu)成的集合就稱為Beamlet詞典。在尺度j上，Beamlet基的數(shù)量為6(2K+J-j-1)2+8(2K+J-j-1)+2。圖3為不同尺度的Beamlet基，圖4為用一個連續(xù)的Beamlet鏈來逼近一條線段，即一條Beam。

2.2 Beamlet變換

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

(6)

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

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

(7)

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

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

fi1i2=Ave{f|Pixel(i1i2)}

(8)

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

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

(9)

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

2.3 離散Beamlet基連接

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

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

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

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

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

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

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

2) 根據(jù)裂縫二值圖像大小選擇要劃分的尺度j(j=2m,j

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

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

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

4 實驗結(jié)果與分析

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

4.1 模擬圖像分析

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

4.2 裂縫檢測

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

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

通過使用如下兩種評價指標對算法精確度進行定性評價：

準確率=正確提取的裂縫長度/真實裂縫長度

誤檢率=錯誤提取的裂縫長度/總的提取裂縫長度

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

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

5 結(jié) 語

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