段友祥，李寧寧，孫歧峰，李 鈺

(中國石油大學(華東) 計算機與通信工程學院，山東 青島 266580)

0 引 言

地質剖面圖是用規(guī)定的符號、花紋和顏色,按一定的比例，沿一定的方向，表示一定的距離內地下一定深度內地質現(xiàn)象的圖件,一般要表示出地層分層、斷層、巖性等地層結構和巖體屬性特征。它形象直觀地表達了地層的結構構造和地層的沉積規(guī)律，是地層在垂向上最直觀、有效的表達方式。油藏地質剖面圖能為系統(tǒng)分析區(qū)域或局部的油藏地質條件、正確指導油藏資源開發(fā)利用和優(yōu)化管理提供依據。

現(xiàn)在很多老油田已經進入了后期開發(fā)，小油藏、隱蔽油藏等復雜油藏成為主要的開發(fā)對象，運用的技術手段也越來越多樣化，要求的基礎地質數(shù)據資料也越來越多。但是，由于受早期條件的限制，這些老油田的很多地質資料不全或不符合現(xiàn)在開發(fā)技術手段的要求。例如，作為重要地質成果圖件的地質剖面圖，過去大多采用手工繪制(很容易轉化為位圖圖像)或軟件繪制位圖圖像(又稱為柵格圖像)。它們雖然能夠逼真展現(xiàn)地層剖面，但是難以運用現(xiàn)代計算機技術進行處理，圖像的數(shù)據量大，不能進行保真縮放，難以滿足現(xiàn)代的鉆井導向、油藏精確描述識別等先進技術的要求。因此，對過去陳舊的地質資料進行深度加工和處理是老油田有效開發(fā)的重要基礎，矢量化就是其中重要的內容之一。矢量化是將油藏地質剖面圖等柵格圖像轉換為矢量圖像的過程，其實質是將圖像數(shù)據轉換為矢量圖形數(shù)據，從而更方便地對圖形進行處理，并精確地表達空間對象的位置、長度等信息[1]。

1 油藏地質剖面圖矢量化方法研究

隨著計算機技術的發(fā)展以及圖像處理技術應用的深入，圖像矢量化對象由早期的二值位圖(即黑白圖片)變?yōu)楝F(xiàn)在的彩色圖像，矢量化算法也在不斷發(fā)展。圖像矢量化的流程一般包括圖像分割、圖像輪廓提取、曲線特征點提取、曲線擬合等主要環(huán)節(jié)。為了更好地實現(xiàn)矢量化，針對各個環(huán)節(jié)存在的突出問題，人們研究和提出了各種不同的解決思路，從而產生了很多不同的矢量化算法，主要包括基于輪廓的矢量化算法、基于邊緣的矢量化算法、基于細化的算法、基于梯度網格的矢量化算法等。

錢曉峰等[2]提出了一種輪廓點列快速矢量化算法,極大地減少了輪廓需存儲的像素點個數(shù),節(jié)省了內存空間,并為進一步處理如形狀匹配、編碼等提供了基礎。黃雪優(yōu)[3]采用Canny算法和基于閾值的算法提取邊緣,對邊緣進行編輯獲取閉合的輪廓，根據輪廓點矢量變化率判斷輪廓的特征點并完成輪廓分段,對分段輪廓進行擬合初步獲得圖像矢量化邊緣。最后采用矢量圖手動編輯的方式對獲得的矢量圖邊緣做出修正,減少特征點個數(shù),提高矢量圖的精度。石文瑩[4]提出了采用基于細化和自適應網格相結合的矢量化方法，解決了傳統(tǒng)基于細化方法容易在交叉處產生畸變的情況。楊云等[5]提出了一種線目標半自動矢量化的方法。該方法通過在用戶選擇的線目標上加一個滑動窗,采用顏色空間轉換、k-均值聚類和區(qū)域生長相結合的方法對窗口內的當前線目標進行自適應分割;再通過不斷地沿前進方向移動窗口和對窗口內的線目標進行分割、細化及序貫跟蹤,來完成線目標的矢量化。王紅巖等[1]針對煤礦地質剖面圖提出了采用不同的算法刪除各種巖性符號，只保留各個地層的邊界線，最后自動矢量化的算法。這些算法在實際問題解決中都得到了應用，一些矢量化軟件采用了其中的部分算法，但是對于一些特殊問題，這些算法還存在不足。

文中在深入研究已有矢量化算法的基礎上，針對油藏地質剖面圖的特點，提出了一種基于顏色聚類結合改進的Canny邊緣檢測的油藏地質剖面圖矢量化方法。

1.1 油藏地質剖面圖矢量化流程

根據矢量化的一般流程，確定地質剖面圖矢量化的流程，如圖1所示。

圖1 地質剖面圖矢量化流程

1.2 顏色聚類分割

油藏地質剖面圖的地層層位之間一般采用不同的圖例填充，其中圖例由顏色和圖案組成。對圖像首先進行顏色空間轉換，然后可以采用顏色聚類消除層位之間的圖例。

1.2.1 RGB轉Lab

常見的顏色模式有RGB模式、HIS模式、HSV模式和Lab模式等。由于采用RGB模式時，各種顏色在數(shù)值上的差距不能很好地反應其色彩差異，于是先將圖像由RGB模式轉換到Lab模式。Lab模式是國際照明委員會于1976年制定的一種色彩模式，是一種均勻的顏色空間，也是經常用來描述人眼可見的所有顏色的最完備的色彩模型。Lab顏色空間由L、a和b三個坐標軸組成。L表示亮度，值域是從0到100逐漸由黑變白；a和b的取值范圍都是-128～+127，a值從小到大的變化是指從綠色變?yōu)榧t色，b值從小到大的變化是從藍色變?yōu)辄S色。一般先將RGB轉化到XYZ空間，然后再由XYZ空間轉化到Lab色彩空間，見式1～3。

(1)

(2)

(3)

設在Lab色彩空間模式下的兩種顏色分別為C1和C2，則C1與C2的顏色色差定義為在L、a、b三個通道的歐氏距離，即：

ED(C1,C2)=

(4)

1.2.2 確定初始聚類中心

假設油藏地質剖面圖的圖例庫為Library(其中Library是由圖例填充時的顏色值L、a、b組成，L、a、b分別為Lab顏色模型的三個分量)，根據圖例庫Library分析油藏地質剖面圖地層輪廓的填充顏色，從而確定初始聚類中心的顏色值。對于給定的油藏地質剖面圖Image(RGB轉Lab后的圖像)中的每一個像素點p,計算其與每種圖例顏色Library(C)(C為圖例庫中的一種顏色)的歐氏距離，每個像素點取對應最小歐式距離的點的顏色值建立顏色量化圖像Image2,即：Image2=Library(C)。其中任意C1≠C，滿足：ED(Library(C),Image(p))

1.2.3 顏色聚類

已知k個初始顏色聚類中心,循環(huán)執(zhí)行以下兩步：(1)對每個像素點p，計算其與每個初始聚類中心的歐氏距離，并將該像素劃分到距離最近的聚類中心上；得到該像素的類別Classp∈{1,2,…,k},即將像素點p劃分到第Classp類。(2)重新計算聚類中心。新的聚類中心為圖例庫Library中與類別平均色最接近的顏色。對于每個類別,先計算該類所包含的所有像素的平均色，然后計算該平均色與Library中的每一種顏色的距離，將距離最小值設為新的聚類中心；直到聚類中心未發(fā)生變化或達到最大迭代次數(shù)。

1.2.4 平 滑

顏色聚類后圖像中還存在一些孤立的顏色塊。采用圖的最小割算法進行平滑處理，通過最小化能量函數(shù)[6]實現(xiàn)，重新分配各像素的類別。最終得到顏色聚類分割后的油藏地質剖面圖。

1.3 改進的自適應Canny邊緣檢測算法

邊緣檢測是圖像處理和計算機視覺中的重要組成部分。邊緣檢測的目的是識別包括圖像特征信息和形狀質量的大部分的亮度變化。眾所周知的經典邊緣檢測算子包括Log、Robert、Prewitt、Laplace、Sobel和Canny[7]等。不同的算子有不同的優(yōu)勢。由于Canny算子具有更好的信噪比和檢測精度，因此，Canny算子成為評估其他方法的基準，其基本步驟如下：

(1)利用一維高斯函數(shù)，對圖像進行低通平滑濾波；

(2)用一階偏導有限差分計算平滑后圖像中各點的梯度值和梯度方向；

(3)對梯度幅值進行非極大值抑制；

(4)用雙閾值算法檢測和連接邊緣。

但Canny算子也存在如下缺點：

(1)傳統(tǒng)的Canny邊緣檢測利用高斯函數(shù)進行平滑濾波時會過度平滑圖像，使圖像變得模糊，可能扭曲和丟失弱邊緣，特別是邊緣的連接和拐角處；同時，高斯函數(shù)的方差需要人為設定；

(2)在計算像素點的梯度幅值時，對噪聲較敏感，容易檢測出假邊緣或丟失一些真實邊緣的細節(jié)[8]；

(3)在高低閾值的選取時，由于高低閾值不是由圖像邊緣的特征信息決定的，而是需要人為事先設定好的，因此不具有自適應性。

針對以上缺點，Cheng Y[9]用非線性擴散濾波器代替高斯函數(shù)平滑圖像，求梯度幅值時考慮了像素對角線方向，最后用Otsu算法自動獲取閾值,在邊緣檢測中取得了更好的精度。Zhao M等[10]針對傳統(tǒng)Canny方法使用高斯函數(shù)去噪時，平滑度取決于高斯核的大小且只對于高斯噪聲效果好的缺點，引進了DCT壓縮域的概念，提出通過DCT變換將2D離散像素轉換為連續(xù)DCT域，得到DCT系數(shù)和估計噪聲，對DCT系數(shù)進行修正，然后通過IDCT獲得平滑圖像的方法，保留了更多有效的邊緣信息。Ma X等[11]為解決傳統(tǒng)Canny邊緣檢測算子過度平滑圖像和適應性弱的問題，改進了參數(shù)Sigma和高低閾值的獲取方法。Biswas R等[12]針對Canny算法對于模糊圖像高低閾值難于確定的缺點，提出了基于type-2模糊集概念的算法來自動確定閾值。尚長春等[13]采用一種新的四階偏微分方程的降噪算法對圖像去噪,進一步提高降噪效果；采用自適應閾值的方法對圖像邊緣進行檢測,實現(xiàn)了雙閾值的自適應提??；基于模糊判決的理論,提出了一種有效的邊緣連接方法。溫陽東等[14]提出了一種基于改進非極大值抑制(NMS)過程和雙閾值求取方法的Canny邊緣檢測算子，能夠獲得更好的圖像邊緣。

在研究以上改進算法的基礎上，結合油藏地質剖面圖矢量化的特點，即對地層水平層位輪廓矢量化，提出了對Canny算法改進的思路，步驟如下：

1.利用一維高斯函數(shù)，對圖像進行低通平滑濾波(采用Sigma為1；Gaussian mask size 為5)。

2.用Sobel算子對圖像進行卷積，求像素點的梯度，采用水平、45°和135°方向上的一階梯度模板，如圖2所示。

圖2 一階梯度模板

三個方向上的一階梯度分量分別為Px(i,j)、P45°(i,j)、P135°(i,j)，由圖2中對應的一階梯度模板與圖像進行卷積得到(其中i,j分別表示圖像像素點的x坐標和y坐標)。每個像素點的梯度幅值和梯度方向計算分別為：

(5)

(6)

但是求完梯度后依然會存在一些噪聲，采用如下方法消除部分噪聲。

(2)Matrix為梯度圖像中每一點的梯度平均值累加和矩陣，其中每一點存放了該點所經過的窗口的梯度平均值的累加和，Matrix1為梯度平均值的累加次數(shù)，計算Matrix中每一點的平均值為Matrix2。

(3)判斷如果M中某一點的梯度值小于Matrix2中對應點梯度值的b倍(b為倍數(shù)因子)，則設該點的梯度值為0，否則該點的梯度值不變。(采用鄰域窗口20行20列，step為5，倍數(shù)因子為2)。

3.對梯度幅值進行非極大值抑制。

4.用雙閾值算法檢測和連接邊緣，其中的高低閾值選取方法采用改進的Otsu自動閾值選取方法[15]，連接邊緣后，刪除邊緣長度小于等于2的邊緣。

1.4 邊緣輪廓矢量化

經過邊緣輪廓提取后的油藏地質剖面圖只包含水平的地層層位輪廓，需要將其矢量化，以便于后期的編輯和應用。采用Two-pass算法先對輪廓信息進行連通區(qū)域標記，然后對標記的每個區(qū)域采用Douglas-Peucker矢量數(shù)據壓縮算法提取特征點(文中采用的DP閾值為1)，最后采用三次貝塞爾曲線對提取的特征點進行擬合，從而完成油藏地質剖面圖的矢量化。

2 實驗結果及分析

2.1 實驗1

為了驗證提出算法的可行性，實驗在MATLAB中以Lena圖像為例，并分別添加椒鹽噪聲(噪聲密度為0.05)、高斯噪聲(均值為0、方差為0．01)以及泊松噪聲(噪聲密度為0.05)來進行驗證，并與傳統(tǒng)經典算法Log、Canny邊緣檢測算法的檢測效果作比較。從定量的角度分析，分別計算各邊緣檢測算子所對應的峰值信噪比(peak signal-to-noise ratio，PSNR)，如表1和表2所示。

表1 不同噪聲環(huán)境下各算法對應的PSNR dB

表2 不同噪聲密度下各算法對應的PSNR dB

上述實驗結果表明，針對不同類型的噪聲以及不同程度的高斯噪聲，文中算法具有較好的抗噪性和準確性。

2.2 實驗2

利用提出的油藏地質剖面圖矢量化流程和改進算法，在MATLAB中用實際地質剖面圖進行了矢量化實驗，如圖3和圖4所示。

圖3 實驗結果1

實驗過程和結果如下：

(1)采用傳統(tǒng)的Canny邊緣檢測算法，對地質剖面圖進行檢測，結果見圖3(b)和圖3(e)。可以看出，圖3(b)噪聲很多，效果很差。采用改進的Canny邊緣檢測算法，對地質剖面圖進行檢測，結果見圖3(c)和圖3(f)?？梢钥闯鲈肼暽?，較好地保留了水平的層位輪廓，效果好。

(2)對原始圖先進行顏色聚類，結果見圖4(a)和圖4(e)；然后再采用改進的Canny邊緣檢測算法進行邊緣檢測，得到輪廓檢測結果，見圖4(b)和圖4(f)；最后再進行矢量化，結果見圖4(c)和圖4(g)。

(3)對圖4(b)和圖4(f)采用美國斯坦福大學的Vector magic軟件進行矢量化，結果見圖4(d)和圖4(h)。可以看出，文中矢量化結果與Vector magic基本相同，但是文中矢量化方法速度更快，得到的曲線特征點更少，更適合矢量化結果的存儲和編輯。與文獻[1]中針對每一種圖例進行刪除的算法相比，文中算法實用性更強。

圖4 實驗結果2

3 結束語

油藏地質剖面圖中的巖性復雜多樣，針對其特點提出了油藏地質剖面圖的矢量化流程，增加了顏色聚類分割，改進了Canny邊緣檢測算法。實驗結果表明，該方法很好地解決了油藏地質剖面的矢量化問題，矢量化速度快、質量高，為在隨鉆地質導向綜合決策中應用油藏地質剖面奠定了基礎。下一步包含斷層的油藏地質剖面的矢量化將是研究的重點。