基于Beamlet變換的數字地質露頭裂縫線性特征提取

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

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

0 引 言

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

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

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

1 算法設計

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

1.1 RGB轉灰度圖像

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

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

(1)

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

1.2 多尺度自適應增強算法

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

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

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

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

rh=Gmean+(Gmax-Gmean)×f

(2)

rl=Gmean-(Gmean-Gmin)×f

(3)

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

I′=I×f′

(4)

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

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

1.3 大津閾值分割算法

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

假設圖像有n個灰度級別，根據灰度區間[0,k]和[k+1,n]將圖像分成A、B兩部分，分別計算A、B區域的平均灰度mA、mB和全圖的平均灰度m，以及A、B區域出現的概率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字典是以不同尺度、方向和位置的小線段為基構成的集合，對一幅分辨率為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基構成的集合就稱為Beamlet詞典。在尺度j上，Beamlet基的數量為6(2K+J-j-1)2+8(2K+J-j-1)+2。圖3為不同尺度的Beamlet基，圖4為用一個連續的Beamlet鏈來逼近一條線段，即一條Beam。

2.2 Beamlet變換

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

(6)

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

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

(7)

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

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

fi1i2=Ave{f|Pixel(i1i2)}

(8)

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

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

(9)

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

2.3 離散Beamlet基連接

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

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

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

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

4) 循環考慮所有“未檢測”狀態的小矩形。

3 Beamlet算法實現具體步驟

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

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

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

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

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

4 實驗結果與分析

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

4.1 模擬圖像分析

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

4.2 裂縫檢測

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

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

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

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

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

四組測試結果如表1所示。從精度評價結果可以看出，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 結 語

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