一種油藏地質(zhì)剖面圖矢量化方法

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

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

0 引 言

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

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

1 油藏地質(zhì)剖面圖矢量化方法研究

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

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

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

1.1 油藏地質(zhì)剖面圖矢量化流程

根據(jù)矢量化的一般流程，確定地質(zhì)剖面圖矢量化的流程，如圖1所示。

圖1 地質(zhì)剖面圖矢量化流程

1.2 顏色聚類分割

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

1.2.1 RGB轉(zhuǎn)Lab

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

(1)

(2)

(3)

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

ED(C1,C2)=

(4)

1.2.2 確定初始聚類中心

假設(shè)油藏地質(zhì)剖面圖的圖例庫為Library(其中Library是由圖例填充時的顏色值L、a、b組成，L、a、b分別為Lab顏色模型的三個分量)，根據(jù)圖例庫Library分析油藏地質(zhì)剖面圖地層輪廓的填充顏色，從而確定初始聚類中心的顏色值。對于給定的油藏地質(zhì)剖面圖Image(RGB轉(zhuǎn)Lab后的圖像)中的每一個像素點p,計算其與每種圖例顏色Library(C)(C為圖例庫中的一種顏色)的歐氏距離，每個像素點取對應(yīng)最小歐式距離的點的顏色值建立顏色量化圖像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中的每一種顏色的距離，將距離最小值設(shè)為新的聚類中心；直到聚類中心未發(fā)生變化或達(dá)到最大迭代次數(shù)。

1.2.4 平 滑

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

1.3 改進(jìn)的自適應(yīng)Canny邊緣檢測算法

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

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

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

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

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

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

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

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

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

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

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

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

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

圖2 一階梯度模板

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

(5)

(6)

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

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

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

3.對梯度幅值進(jìn)行非極大值抑制。

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

1.4 邊緣輪廓矢量化

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

2 實驗結(jié)果及分析

2.1 實驗1

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

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

表2 不同噪聲密度下各算法對應(yīng)的PSNR dB

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

2.2 實驗2

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

圖3 實驗結(jié)果1

實驗過程和結(jié)果如下：

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

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

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

圖4 實驗結(jié)果2

3 結(jié)束語

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