本發(fā)明涉及一種油藏地質(zhì)剖面圖矢量化方法，屬于計(jì)算機(jī)輔助設(shè)計(jì)技術(shù)領(lǐng)域，具體地說是一種基于改進(jìn)的自適應(yīng)canny(canny，johnf.canny于1986年開發(fā)出來的一個(gè)多級邊緣檢測算法)邊緣檢測的油藏地質(zhì)剖面圖矢量化方法。

背景技術(shù)：

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

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

例如文獻(xiàn)1：selingerp.potrace:apolygon-basedtracingalgorithm[j].2003.中提出了基于多邊形的位圖輪廓矢量化算法potrace，該算法的思想是把位圖分解為一些路徑，將每條路徑都近似為一個(gè)最優(yōu)多邊形，由多邊形轉(zhuǎn)化為光滑的輪廓，最后由貝塞爾曲線進(jìn)行擬合產(chǎn)生矢量化圖像。

文獻(xiàn)2：石文瑩.圖像特征閾值分割及其矢量化的研究[d].合肥工業(yè)大學(xué),2012.中提出了采用基于細(xì)化和自適應(yīng)網(wǎng)格相結(jié)合的矢量化方法，解決了傳統(tǒng)基于細(xì)化方法容易在交叉處產(chǎn)生畸變的情況。

文獻(xiàn)3：王紅巖,毛善君.煤礦地質(zhì)剖面圖的計(jì)算機(jī)矢量化方法探討[j].計(jì)算機(jī)工程與科學(xué),2005,(10):35-36+68.針對于煤礦地質(zhì)剖面圖提出了采用不同的算法刪除各種巖性符號，只保留各個(gè)地層的邊界線，最后自動(dòng)矢量化的方法，該方法需要采用不同的算法刪除不同的巖性符號，但是地層巖性符號眾多，不可能處理每一種巖性符號，實(shí)用性差。

專利cn201610051086公布了一種基于色彩分割的控規(guī)圖紙自動(dòng)矢量化方法,該方法將圖像轉(zhuǎn)換為hsv空間進(jìn)行匹配分割，將分割后的結(jié)果進(jìn)行二值化，最后進(jìn)行連通體抽取，得到矢量坐標(biāo)。該方法可以快速分析圖紙,完成整個(gè)矢量圖紙的轉(zhuǎn)換過程,但是對于噪聲多的圖像的輪廓信息存在丟失現(xiàn)象。上述方法在實(shí)際問題解決中都得到了應(yīng)用，一些矢量化軟件采用了其中的部分算法，但是對于地層剖面圖的矢量化，這些算法還存在不足，不能滿足地層剖面圖矢量化的要求。

技術(shù)實(shí)現(xiàn)要素：

本發(fā)明目的:針對油藏地質(zhì)剖面圖的特點(diǎn)，提出了顏色聚類結(jié)合改進(jìn)的自適應(yīng)canny邊緣檢測算法的一種油藏地質(zhì)剖面圖矢量化方法，用于支持油藏地質(zhì)剖面圖的矢量化。

一種油藏地質(zhì)剖面圖矢量化方法，包括以下步驟：

步驟s1：運(yùn)用matlab(matlab，美國mathworks公司出品的商業(yè)數(shù)學(xué)軟件)軟件，通過將direct(direct，即“數(shù)字化油藏表征軟件系統(tǒng)”，direct為其英文全稱“digitalreservoircharacterizationtool”的縮寫)數(shù)字化油藏表征軟件系統(tǒng)輸出的地質(zhì)剖面圖作為輸入，將該油藏地質(zhì)剖面圖根據(jù)本文的矢量化方法獲得矢量化結(jié)果；

步驟s2：對步驟s1中的輸入進(jìn)行預(yù)處理，采用高斯濾波對圖像進(jìn)行平滑濾波，然后在平滑的基礎(chǔ)上采用拉普拉斯銳化方法來增強(qiáng)圖像；

步驟s3：對預(yù)處理后的油藏地質(zhì)剖面圖由rgb(rgb，代表紅、綠、藍(lán)三個(gè)通道的顏色)顏色空間模型轉(zhuǎn)為lab(lab，lab色彩模型中明度l和有關(guān)色彩的三個(gè)要素)顏色空間模型，再進(jìn)行顏色聚類分割，得出初始油藏地質(zhì)剖面圖；

步驟s4：采用改進(jìn)的自適應(yīng)canny邊緣檢測算法對所述初始油藏地質(zhì)剖面圖進(jìn)行邊緣檢測，獲取所述初始油藏地質(zhì)剖面圖的邊緣輪廓；

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

步驟s42：用sobel(sobel，由irwinsobel提出的一種一階梯度算子)算子對圖像進(jìn)行卷積求像素點(diǎn)的梯度，采用x軸、45°和135°方向上的一階梯度模板，模板如下所示：；三個(gè)方向上的一階梯度分量分別為gx(i,j)、g45°(i,j)、g135°(i,j)，由三個(gè)對應(yīng)的一階梯度模板與圖像進(jìn)行卷積得到，其中(i,j)表示圖像中i行j列的像素點(diǎn)，每個(gè)像素點(diǎn)的梯度幅值為：梯度方向?yàn)椋?imgfile="bda0001320939990000033.gif"wi="518"he="136"img-content="drawing"img-format="gif"orientation="portrait"inline="no"/>gy(i,j)為y軸方向的一階梯度分量；

步驟s43：依據(jù)大小為n行m列的梯度圖像m，計(jì)算以(i,j)為起點(diǎn)，大小為(i+a)行(j+a)列的窗口的梯度的累加和sum,avg為點(diǎn)(i,j)在該窗口的梯度累加和的平均值，則其中i,j分別表示圖像i行j列的像素點(diǎn),i∈[0，n-1-a],j∈[0，m-1-a]，其中a為取某一像素點(diǎn)的鄰域窗口的大小，鄰域窗口的移動(dòng)步長設(shè)為step；

步驟s44：矩陣matrix為梯度圖像m中每一點(diǎn)的梯度平均值累加和矩陣，其中每一點(diǎn)存放了該點(diǎn)所經(jīng)過的窗口的梯度平均值的累加和，矩陣matrix1為每一點(diǎn)的梯度平均值的累加次數(shù)，計(jì)算matrix中每一點(diǎn)的平均值為矩陣matrix2；

步驟s45：判斷如果m中某一點(diǎn)的梯度值小于matrix2中對應(yīng)點(diǎn)梯度值的b倍，b為倍數(shù)因子，則該點(diǎn)為梯度值為0，否則該點(diǎn)的梯度值保持原值不變；

步驟s46：對梯度幅值進(jìn)行非極大值抑制；

步驟s47：用雙閾值算法檢測和連接邊緣，得到最終的邊緣圖像；

步驟s5：對所述始油藏地質(zhì)剖面圖的邊緣輪廓進(jìn)行矢量化，得到目標(biāo)油藏地質(zhì)剖面圖；

步驟s51：對邊緣檢測的結(jié)果進(jìn)行連通區(qū)域標(biāo)記；

步驟s52：對標(biāo)記的每個(gè)區(qū)域采用douglas–peucker(douglas–peucker，亦稱為拉默-道格拉斯-普克算法、迭代適應(yīng)點(diǎn)算法、分裂與合并算法)矢量數(shù)據(jù)壓縮算法提取特征點(diǎn)；

步驟s53：由三次貝塞爾曲線對提取的特征點(diǎn)進(jìn)行擬合，完成油藏地質(zhì)剖面圖的矢量化；

步驟s6：輸出矢量化結(jié)果。

附圖說明

附圖1為本發(fā)明的處理流程示意圖。

附圖2為輸入的油藏地質(zhì)剖面圖。

附圖3為附圖2顏色聚類的結(jié)果。

附圖4為附圖3經(jīng)改進(jìn)的自適應(yīng)canny檢測的結(jié)果。

附圖5為附圖4矢量化的結(jié)果。

具體實(shí)施方式

下面結(jié)合附圖1-5，運(yùn)用matlab軟件，對本發(fā)明作進(jìn)一步的描述：一種采用改進(jìn)的自適應(yīng)canny邊緣檢測算法矢量化油藏地質(zhì)剖面圖的矢量化方法，具體實(shí)施方案包括如下步驟，步驟如圖1所示：

(1)位圖圖像輸入。在matlab中讀入圖2作為初始的輸入；

(2)圖像預(yù)處理。輸入的油藏地質(zhì)剖面圖圖2中除了代表地層界線的線段外，還有大量充填于地層界線間的巖性符號，而且這些巖性符號形狀多樣(如充填區(qū)域、短線、圓、實(shí)心圓、小三角形等)且含有部分噪聲，對圖2進(jìn)行預(yù)處理，首先采用高斯濾波對圖像進(jìn)行平滑操作，然后在平滑的基礎(chǔ)上采用拉普拉斯銳化方法來增強(qiáng)圖像，通過圖像預(yù)處理可以去除圖像中的部分噪聲。

(3)圖像分割。油藏地質(zhì)剖面圖的矢量化是針對于圖中的水平層位輪廓信息的矢量化，而圖中的各種圖例不需要保留，首先對(2)中預(yù)處理后的圖像進(jìn)行顏色空間轉(zhuǎn)換，由rgb顏色空間模型轉(zhuǎn)為lab顏色空間模型，lab顏色空間模型有l(wèi)、a和b三個(gè)坐標(biāo)軸組成。l表示亮度，值域是從0到100逐漸由黑變白；a和b的取值范圍都是-128到+127，a值從小到大的變化是指從綠色變?yōu)榧t色，b值從小到大的變化是從藍(lán)色變?yōu)辄S色。一般先將rgb空間轉(zhuǎn)化到xyz(xyz，xyz顏色模型中的三個(gè)分量)空間然后再由xyz空間轉(zhuǎn)化到lab色彩空間然后輸入圖像聚類的初始聚類中心數(shù)k為7，假設(shè)油藏地質(zhì)剖面圖的圖例庫為library，其中l(wèi)ibrary是由圖例填充時(shí)的顏色值l、a、b組成，l、a、b分別為lab顏色模型的三個(gè)分量，根據(jù)圖例庫library分析油藏地質(zhì)剖面圖地層輪廓的填充顏色，從而確定k個(gè)初始聚類中心的顏色值。對于顏色空間轉(zhuǎn)換后的油藏地質(zhì)剖面圖中的每一個(gè)像素點(diǎn)p,計(jì)算其與每種圖例顏色library(c)的歐式距離，其中c為圖例庫中的一種顏色，每個(gè)像素點(diǎn)取對應(yīng)最小歐式距離的點(diǎn)的顏色值建立顏色量化圖像image2,即image2＝library(c)其中滿足ed(library(c),image(p))<ed(library(c1),image(p)),其中設(shè)在lab色彩空間模式下的兩種顏色分別為c1，c2，則c1與c2的顏色色差定義為在l、a、b三個(gè)通道的歐式距離。統(tǒng)計(jì)image2中各種顏色所占的比例，并按升序排列，得到image2的顏色直方圖histogram。取出前n種顏色，n為滿足的最小n值，histi表示直方圖中的第i種顏色的占比，i取值1～n，即保證所有取出的顏色總占比不小于90％的同時(shí)，使取出的顏色最少。設(shè)聚類數(shù)為k,對得到的n種顏色采用層次聚類算法進(jìn)行聚類，得到k類的中心色即為圖像的k個(gè)初始的顏色聚類中心color{1,…k}。

已知k個(gè)初始顏色聚類中心。重復(fù)以下步驟：1)每個(gè)像素點(diǎn)p，計(jì)算其與每個(gè)初始聚類中心的歐式距離，并將該像素劃分到距離最近的聚類中心上；得到該像素的類別colorp{1,2,…,k},即將像素點(diǎn)p劃分到第colorp類；2)重新計(jì)算聚類中心。新的聚類中心為圖例庫library中與類別平均色最接近的顏色。對于每個(gè)類別,先計(jì)算該類所包含的所有像素的平均色，然后計(jì)算該平均色與library中的每一種顏色的距離，將距離最小值設(shè)為新的聚類中心；直到聚類中心未發(fā)生變化或達(dá)到最大迭代次數(shù)時(shí)終止。得到如圖3所示的聚類結(jié)果。

從圖3的聚類結(jié)果中可以得出本文采用的圖像分割方法較好地消除了地層之間的圖例，得到了準(zhǔn)確的地層層位輪廓信息。

(4)邊緣檢測。采用改進(jìn)的自適應(yīng)canny邊緣檢測算法，對圖3進(jìn)行邊緣檢測，設(shè)圖3為f(x,y)；首先利用一維高斯函數(shù):對圖像進(jìn)行低通平滑濾波處理其中表示卷積操作，卷積核大小為5行5列，σ取2，x，y表示x坐標(biāo)和y坐標(biāo)，得到處理后的圖像g(x,y)；然后采用sobel算子對圖像g(x,y)根據(jù)步驟s42中一階梯度模板進(jìn)行卷積求像素點(diǎn)的梯度分量，然后將梯度分量帶入步驟s42所述的梯度公式中計(jì)算梯度，得到梯度圖像m(x,y),計(jì)算完梯度后，仍然會(huì)存在噪聲，根據(jù)步驟s43、步驟s44和步驟s45所述的方法消除部分梯度噪聲，其中a為20，step為5，b為2；然后對梯度幅值進(jìn)行非極大值抑制，首先根據(jù)步驟s42中所述的梯度方向公式計(jì)算每個(gè)坐標(biāo)點(diǎn)的梯度方向，以待判斷像素點(diǎn)為中心分為四個(gè)區(qū)域，代表梯度方向近似的四個(gè)可能角度分別為0°，45°，90°，135°，對每一個(gè)像素點(diǎn)的梯度進(jìn)行非極大值抑制，首先將待判斷像素梯度值m(x,y)與沿著梯度方向的2個(gè)8鄰域像素的梯度m(x±1,y±1)進(jìn)行比較，如果該點(diǎn)像素點(diǎn)位置的梯度幅值沒有沿梯度方向的兩個(gè)相鄰像素的梯度值大，則說明該點(diǎn)的梯度值不是局部最大值需被抑制，令m(x,y)＝0，即該點(diǎn)不是邊緣點(diǎn)。用此方法進(jìn)行抑制全局梯度幅值中的非極大值，這一步驟會(huì)排除非邊緣點(diǎn)，保留候選的邊緣；最后用雙閾值算法檢測和連接邊緣，對非極大值抑制后的圖像與高低閾值進(jìn)行判斷，如果某一像素位置的梯度幅值超過高閾值，則該像素被保留為邊緣像素；如果某一像素位置的梯度幅值小于低閾值，則該像素被排除；如果某一像素位置的梯度幅值介于高、低閾值之間，則判斷該像素8鄰域空間的像素是否存在高于高閾值的像素，如果存在，則該像素將被保留。其中高低閾值通過otsu(otsu，一種使類間方差最大的自動(dòng)確定閾值方法)算法計(jì)算，邊緣檢測的結(jié)果如圖4所示，其中高閾值為144，低閾值為77；

從圖4可以看出本文采用的圖像邊緣檢測方法對分割后的圖像進(jìn)行檢測的結(jié)果保留了所有的邊緣輪廓，未出現(xiàn)邊緣丟失和斷裂的現(xiàn)象。

(5)層位線矢量化。首先采用two-pass(two-pass，一種連通域標(biāo)記算法)算法對圖4進(jìn)行連通區(qū)域標(biāo)記，得到8個(gè)連通區(qū)域；然后對標(biāo)記的每個(gè)區(qū)域采用douglas–peucker矢量數(shù)據(jù)壓縮算法提取特征點(diǎn)，其中閾值threshold為1；最后采用三次貝塞爾曲線對提取的特征點(diǎn)進(jìn)行擬合；圖5為矢量化結(jié)果，圖5在保留了原有的層位輪廓變化趨勢的同時(shí)，得到的邊緣輪廓的特征點(diǎn)更少，便于矢量化數(shù)據(jù)的存儲和編輯。

(6)矢量圖輸出。將擬合的結(jié)果保存并輸出，從而完成了整個(gè)矢量化流程，結(jié)果如圖5所示。

本發(fā)明的有益效果是：該發(fā)明較好的分割了地層剖面圖，保留了原有的地層中各層的層位信息；改進(jìn)的自適應(yīng)canny邊緣檢測算法提取的輪廓線準(zhǔn)確，保留了原圖中的層位趨勢；得到的矢量化結(jié)果在保留了各層的層位信息的基礎(chǔ)上，得到的層位曲線的特征點(diǎn)少，易于矢量化數(shù)據(jù)的存儲和編輯。

本發(fā)明不局限于上述具體實(shí)施方式，根據(jù)上述內(nèi)容，按照本領(lǐng)域的普通技術(shù)知識和慣用手段，在不脫離本發(fā)明上述基本技術(shù)思想前提下，本發(fā)明還可以做出其它多種形式的等效修改、替換或變更，均屬于本發(fā)明的保護(hù)范圍。