基于路徑形態(tài)學(xué)和正弦函數(shù)族匹配的電成像測(cè)井縫洞自動(dòng)識(shí)別與分離方法研究

王磊 ，沈金松 *，衡海亮 ，魏帥帥

1 中國(guó)石油大學(xué)(北京)油氣資源與探測(cè)國(guó)家重點(diǎn)實(shí)驗(yàn)室，北京 102249

2 中國(guó)石油大學(xué)(北京)CNPC物探重點(diǎn)實(shí)驗(yàn)室，北京 102249

0 引言

碳酸巖鹽儲(chǔ)層在各種類(lèi)型的巖性?xún)?chǔ)層中的占比較大，從已公布數(shù)據(jù)來(lái)看，碳酸巖鹽儲(chǔ)層的開(kāi)采產(chǎn)出可達(dá)全球范圍內(nèi)油氣產(chǎn)出的六成左右[1-4]。該類(lèi)型油氣藏裂縫較發(fā)育，因此對(duì)縫洞型碳酸鹽巖儲(chǔ)層的研究至關(guān)重要。該類(lèi)型儲(chǔ)層的油氣主要存在交錯(cuò)的裂縫、經(jīng)風(fēng)化腐蝕和流體經(jīng)過(guò)形成的孔洞之中，其構(gòu)成成分有原生基質(zhì)孔隙及次生裂縫、溶孔、溶洞[5]?？紫缎?、溶洞型、裂縫型為基本儲(chǔ)集類(lèi)型，孔洞型、縫洞型、孔縫型、孔縫洞復(fù)合型為過(guò)渡類(lèi)型[5-6]。正是由于縫洞型儲(chǔ)層的復(fù)雜性，縫洞的測(cè)井識(shí)別與評(píng)價(jià)仍存在多方面的難題，例如，缺乏儲(chǔ)層類(lèi)型準(zhǔn)確識(shí)別和劃分有效方法，缺少儲(chǔ)集空間定量表征的高精度方法。因此，實(shí)際生產(chǎn)亟需研究縫洞自動(dòng)識(shí)別與分離算法。

縫洞型儲(chǔ)層并非今天才進(jìn)入人們的視線(xiàn)。李瑞、陸云龍等通過(guò)井口直徑的計(jì)算方法結(jié)合RLLD/RLLS理論，歸納出裂縫出現(xiàn)的可能性的值，并以此作為儲(chǔ)層含油性評(píng)價(jià)的標(biāo)準(zhǔn)[7-8]。楊輝廷等利用測(cè)井及鉆井?dāng)?shù)據(jù)對(duì)塔里木盆地裂縫和溶蝕孔洞信息進(jìn)行了較為準(zhǔn)確的計(jì)算，對(duì)與水平地層呈不同角度的裂縫進(jìn)行了提取，對(duì)孔洞飽含油、氣、水的分類(lèi)通過(guò)自然伽馬算法進(jìn)行了劃分[9]。陳科貴等從巖層內(nèi)電流經(jīng)過(guò)縫洞和非縫洞部分的串并特征進(jìn)行分析，試圖以此來(lái)計(jì)算巖層的滲透率，并認(rèn)為巖層的滲透率不是一成不變的[10]。

電成像測(cè)井對(duì)碳酸鹽巖縫洞型儲(chǔ)層的評(píng)價(jià)作用在油氣業(yè)界逐漸獲得越來(lái)越多的關(guān)注，國(guó)內(nèi)外學(xué)者和油氣服務(wù)公司發(fā)表了許多成果。為了將裂縫和孔洞從基質(zhì)中區(qū)分開(kāi)來(lái)并計(jì)算獲得參數(shù)，斯倫貝謝公司的Delhomme首先發(fā)表了成像的閾值法，這可以將縫洞的外形顯示出來(lái)，同時(shí)指出裂縫所在深度段[11]；Hall等嘗試使用Hough變換計(jì)算正弦型的裂縫信息和地質(zhì)走向數(shù)據(jù)，就符合正弦分布的連續(xù)裂縫取得了較好效果[12]。景建恩等針對(duì)塔河油田碳酸鹽巖儲(chǔ)層的裂縫，通過(guò)對(duì)FMI測(cè)井圖像數(shù)據(jù)使用R型聚類(lèi)算法，建立了裂縫指標(biāo)數(shù)據(jù)集合，并對(duì)裂縫出現(xiàn)可能性做出了客觀(guān)評(píng)價(jià)[13]。李曦寧、李振苓和沈金松等提出了基于二值圖的完備路徑形態(tài)學(xué)算法，用于分離地層層理與裂縫[14-15]。除此之外，人工智能算法在測(cè)井領(lǐng)域的應(yīng)用實(shí)例也越來(lái)越多[16-22]。針對(duì)電成像測(cè)井，XUEYongchao等[19]利用遺傳算法與神經(jīng)網(wǎng)絡(luò)結(jié)合對(duì)多種測(cè)井資料進(jìn)行處理進(jìn)而實(shí)現(xiàn)裂縫識(shí)別。Azim[20]提出一種魯棒性高的神經(jīng)網(wǎng)絡(luò)算法(ANN)，利用Fortran語(yǔ)言建立以后向傳播為基礎(chǔ)的訓(xùn)練過(guò)程，引入5口井的微電阻率成像及傳統(tǒng)測(cè)井曲線(xiàn)數(shù)據(jù)的80%作為訓(xùn)練集，20%作為測(cè)試集，以此來(lái)估計(jì)裂縫密度和分形兩個(gè)維度上的裂縫性質(zhì)，預(yù)測(cè)儲(chǔ)層段井壁附近的三維裂縫分布圖。李冰濤等[21]提出基于TensorFlow的圖像語(yǔ)義分割模型DeepLabv3+來(lái)對(duì)經(jīng)由Labelme工具標(biāo)定后的裂縫數(shù)據(jù)進(jìn)行像素分割和提取，效果較好。鄒文波對(duì)人工智能在巖性識(shí)別、低電阻率油層識(shí)別、儲(chǔ)層參數(shù)評(píng)價(jià)、儲(chǔ)層裂縫孔洞識(shí)別等4個(gè)測(cè)井重點(diǎn)方向的研究現(xiàn)狀進(jìn)行了分析和總結(jié)[22]。然而，由于受復(fù)雜地質(zhì)環(huán)境和井孔干擾的影響，電導(dǎo)率圖像的縫洞自動(dòng)識(shí)別和定量評(píng)價(jià)仍存在諸多困難，例如受地層或極板切割的裂縫或孔洞因?yàn)槠洳贿B續(xù)性，對(duì)其進(jìn)行定量評(píng)價(jià)仍較為困難。

本文根據(jù)數(shù)學(xué)形態(tài)學(xué)基本理論和不完備路徑算子適應(yīng)于彎曲線(xiàn)性結(jié)構(gòu)追蹤以及對(duì)間斷有較強(qiáng)容忍度的特點(diǎn)，將不完備路徑形態(tài)學(xué)算法、正弦函數(shù)族算法和二階矩橢圓擬合算法應(yīng)用于電導(dǎo)率測(cè)井圖像，通過(guò)給定長(zhǎng)度的路徑算子、容忍算子和正弦函數(shù)族數(shù)據(jù)集完成對(duì)裂縫參數(shù)的準(zhǔn)確計(jì)算，形成了一種不完備路徑形態(tài)學(xué)算法與正弦函數(shù)族結(jié)合計(jì)算裂縫參數(shù)信息以及用二階矩自動(dòng)對(duì)溶蝕孔洞進(jìn)行橢圓擬合并計(jì)算孔洞參數(shù)信息的方法。為測(cè)試和驗(yàn)證該方法的有效性和適應(yīng)性，將所研究方法應(yīng)用于某油田多口井的實(shí)測(cè)電成像數(shù)據(jù)，初步證實(shí)了算法的高精度和有效性。

1 用于電導(dǎo)率圖像基質(zhì)和縫洞電導(dǎo)率異常自動(dòng)分離的Otsu算法

在這一部分，應(yīng)用大津法(Otsu)從井壁采集獲得的FMI圖像中自動(dòng)分離裂縫和孔洞[23]。Otsu算法的基本思想是指定一個(gè)閾值，該閾值可以參考圖像灰度圖的像素分布柱狀圖以便將圖像分成兩部分，且兩部分的方差達(dá)到最大。

如圖1所示，設(shè)計(jì)了一個(gè)圖像G，它是Q×R像素大小的圖像，

設(shè)g(x,y)是分布在(x,y)范圍內(nèi)的圖像的灰度函數(shù)。通過(guò)在灰度圖的像素值區(qū)間里測(cè)試不同的閾值來(lái)獲得最優(yōu)的圖像分割閾值thresh，該閾值可以讓被分割成兩部分的方差σ2(thresh)達(dá)到最大，該算法的具體實(shí)現(xiàn)如圖1所示。

圖1 Otsu算法示意圖Fig.1 A schematic diagram of Otsu method

對(duì)于二維電導(dǎo)率圖像來(lái)說(shuō)，可以把閾值thresh作為分割的參照。高于thresh的所有像素點(diǎn)設(shè)定為1，該像素點(diǎn)對(duì)應(yīng)于井壁巖石的基質(zhì)部分，低于thresh的所有像素點(diǎn)設(shè)定為0，該像素點(diǎn)對(duì)應(yīng)于井壁巖石的裂縫、溶蝕孔或洞部分。通過(guò)這一過(guò)程，灰度電導(dǎo)率測(cè)井圖像矩陣變成了0-1二值矩陣。圖2a是原始灰度縫洞成像圖，圖2b是基于小波變換與快速行進(jìn)算法插值后的圖像[24]；圖2c是Otsu分割后的裂縫、溶蝕孔或洞圖像。圖中可以看出，裂縫、溶蝕孔或洞可以較好地與基質(zhì)部分實(shí)現(xiàn)分離。

2 不完備路徑形態(tài)學(xué)縫洞識(shí)別方法

圖2所示的電成像測(cè)井電導(dǎo)率圖像中暗色或黑色斑塊或條帶顯示了可能存在的裂縫或溶蝕孔洞，如何從這些圖像中識(shí)別和分離出該類(lèi)電導(dǎo)率異常是電成像測(cè)井縫洞儲(chǔ)層識(shí)別與評(píng)價(jià)的首要任務(wù)。下文首先簡(jiǎn)單介紹從電導(dǎo)率圖像中提取高電導(dǎo)率異常并分離縫洞的方法原理。

圖2 (a)原始灰度縫洞成像圖；(b)基于小波變換與快速行進(jìn)算法插值后的圖像；(c)Otsu分割后的縫洞圖像Fig.2 (a) Original electric imaging logging image; (b)Electric imaging image interpolated by means of wavelet transform and fast matching method; (c) Fracture-vug image separated by the Otsu method

2.1 基于二值圖路徑開(kāi)的裂縫識(shí)別

路徑形態(tài)學(xué)與傳統(tǒng)的數(shù)學(xué)形態(tài)學(xué)開(kāi)閉運(yùn)算不同，包含了諸多與路徑相關(guān)的概念。首先是帶方向的臨近關(guān)系的概念。設(shè)E是一個(gè)已知節(jié)點(diǎn)的集合，代表像素點(diǎn)的位置。使用二元的臨近關(guān)系“?”來(lái)定義一個(gè)在這些節(jié)點(diǎn)之上的有向圖。具體地，x?y表示存在一條從x到y(tǒng)的路徑。如果x?y，就稱(chēng)y是x的后繼節(jié)點(diǎn)，x是y的前驅(qū)節(jié)點(diǎn)。這個(gè)概念在圖3中有展示。b1、b2、b3是a的前驅(qū)節(jié)點(diǎn)，a1、a2、a3是b的后繼節(jié)點(diǎn)。圖3中可以看到，箭頭表示的臨近關(guān)系規(guī)定了節(jié)點(diǎn)與節(jié)點(diǎn)之間的繼承方向是從左至右，分別為45°，0°，-45°，而臨近關(guān)系的角度并不局限于45°的倍數(shù)，可以任意選擇。在FMI圖像中應(yīng)用該方法時(shí)，通過(guò)事先定義臨近關(guān)系，可以把結(jié)點(diǎn)間以不同的方向串連起來(lái)，提取與實(shí)際測(cè)井資料相符的圖像特征。

圖3 鄰接關(guān)系示意圖，b1、b2、b3是a的三個(gè)前驅(qū)節(jié)點(diǎn)，a1、a2、a3是b的三個(gè)后繼節(jié)點(diǎn)Fig.3 Adjacency relation of set a and b, where a is the successor of b1, b2, b3, and b is the predecessor of a1, a2, a3.

如果存在一個(gè)空間分布結(jié)構(gòu)恒定不變的二維節(jié)點(diǎn)圖，當(dāng)若干節(jié)點(diǎn)滿(mǎn)足圖4所示的臨近關(guān)系，那么其中會(huì)包括多條不同形狀的長(zhǎng)度為L(zhǎng)的路徑。由于每個(gè)點(diǎn)的前驅(qū)結(jié)點(diǎn)與后繼結(jié)點(diǎn)的數(shù)量均為3個(gè)，以任一結(jié)點(diǎn)為起始點(diǎn)的長(zhǎng)度為L(zhǎng)的路徑會(huì)有3L-1種情況，因此，路徑形態(tài)學(xué)算法的耗時(shí)與L呈指數(shù)型上升，同時(shí)對(duì)內(nèi)存的要求也會(huì)逐步提高。針對(duì)這一問(wèn)題，Talbot等提出了一種新的快速算法[25-28]。假設(shè)存在東西方向的二維臨近關(guān)系節(jié)點(diǎn)集合E，X是節(jié)點(diǎn)圖E中的任意子集。給定函數(shù)b：

對(duì)于二維臨近關(guān)系節(jié)點(diǎn)集合E，在每個(gè)像素點(diǎn)p存儲(chǔ)兩個(gè)值：起始點(diǎn)為p且方向?yàn)橛业淖铋L(zhǎng)路徑的長(zhǎng)度λ+[p](該長(zhǎng)度不包含點(diǎn)p自身)，以及起始點(diǎn)為p且方向?yàn)樽蟮淖铋L(zhǎng)路徑的長(zhǎng)度λ-[p](該長(zhǎng)度不包含點(diǎn) 自身)，其計(jì)算公式如下：

其中，p1，p2表示二維節(jié)點(diǎn)圖中點(diǎn)p的橫縱坐標(biāo)。

當(dāng)二維臨近關(guān)系節(jié)點(diǎn)圖滿(mǎn)足圖4(a)且b[p]=真時(shí)，可以通過(guò)對(duì)每個(gè)符合條件的路徑上的節(jié)點(diǎn)重復(fù)利用式(3)和(4)的關(guān)系，分別計(jì)算得到節(jié)點(diǎn)圖中每個(gè)像素點(diǎn)p(p1,p2)的向右最長(zhǎng)路徑節(jié)點(diǎn)數(shù)λ+[p]、向左最長(zhǎng)路徑節(jié)點(diǎn)數(shù)λ-[p]，進(jìn)而根據(jù)式(5)求得包含節(jié)點(diǎn)p的最長(zhǎng)路徑總的節(jié)點(diǎn)數(shù)λ[p]：

圖4 四種鄰接關(guān)系的組合Fig.4 Four kinds of adjacency relations.

如果點(diǎn)p自身是空缺態(tài)，即b(p)=假，那么穿過(guò)點(diǎn)p的最長(zhǎng)路徑不存在，設(shè)為λ[p]=0。

圖5a～e是一個(gè)路徑形態(tài)學(xué)算法在一個(gè)8×9的二維臨近關(guān)系節(jié)點(diǎn)矩陣中的示意圖。圖中黑色節(jié)點(diǎn)表示b[p]=真，灰色節(jié)點(diǎn)表示b[p]=假，臨近關(guān)系滿(mǎn)足圖4a。圖5a表示從右側(cè)第一列開(kāi)始計(jì)算經(jīng)過(guò)該列節(jié)點(diǎn)的最長(zhǎng)路徑長(zhǎng)度，并逐列向左推進(jìn)直至最左列(如綠框所示)，所有符合臨近關(guān)系的黑色節(jié)點(diǎn)都用綠色箭頭連接起來(lái)并在圖5c中統(tǒng)計(jì)經(jīng)過(guò)各節(jié)點(diǎn)的右向最長(zhǎng)路徑長(zhǎng)度λ-[p]。圖5b表示從左向右計(jì)算圖5d中λ+[p]的過(guò)程(如藍(lán)框所示)。圖5e按式(5)計(jì)算路徑開(kāi)的值λ[p]。此時(shí)如想獲得長(zhǎng)度大于L的路徑，只需從圖5e中找出符合條件的節(jié)點(diǎn)，這就是路徑開(kāi)運(yùn)算追蹤裂縫分布的原理。

圖5 二維臨近關(guān)系路徑開(kāi)運(yùn)算示意圖Fig.5 The diagram of the path opening algorithm in a 2-D adjacency relations

2.2 基于二維節(jié)點(diǎn)圖的不完備路徑開(kāi)的裂縫識(shí)別運(yùn)算

在二維節(jié)點(diǎn)圖中使用長(zhǎng)度因子L進(jìn)行路徑開(kāi)運(yùn)算時(shí)，各節(jié)點(diǎn)間的臨近關(guān)系是逐點(diǎn)相連的，路徑中不存在節(jié)點(diǎn)缺失的情況，而實(shí)際的裂縫常因地層運(yùn)動(dòng)、儀器與井壁耦合性不良或含有不規(guī)則噪聲導(dǎo)致測(cè)井圖像中某一條裂縫從中間斷開(kāi)，形成多條路徑長(zhǎng)度較短的小裂縫。因此若只使用長(zhǎng)度因子進(jìn)行裂縫的提取，當(dāng)L取值較大時(shí)，會(huì)大大降低對(duì)有間斷裂縫的提取效果。

針對(duì)這一問(wèn)題，可以從Talbot等改進(jìn)的不完備路徑形態(tài)學(xué)算法獲得借鑒，該算法引入了容忍因子K，當(dāng)某一條路徑上有K個(gè)節(jié)點(diǎn)丟失時(shí)(b[p]=假)，使用路徑長(zhǎng)度因子L仍可以對(duì)間斷路徑進(jìn)行有效識(shí)別和提取，尤其是對(duì)背景包含較多點(diǎn)狀噪聲的路徑提取效果較好。容忍因子的使用示例如圖6所示。圖6a是包含間斷的裂縫，可以看到裂縫右側(cè)有2個(gè)節(jié)點(diǎn)丟失；圖6b是L=6時(shí)的完備路徑開(kāi)的提取結(jié)果，右側(cè)裂縫的路徑長(zhǎng)度為5，不符合路徑長(zhǎng)度條件，無(wú)法有效提??；圖6c是L=6，K=3時(shí)的不完備路徑開(kāi)的提取結(jié)果，由于容忍因子K的存在，丟失的2個(gè)節(jié)點(diǎn)沒(méi)有影響到裂縫右側(cè)部分的提取。由此可見(jiàn)路徑長(zhǎng)度因子L與容忍因子K同時(shí)使用可以提高間斷裂縫的識(shí)別效果。

圖6 基于不完備路徑的裂縫提取示意圖Fig.6 The diagram of the fracture extraction based on the incomplete path opening method

用λ[p,k]表示節(jié)點(diǎn)p在缺失路徑長(zhǎng)度為k即路徑容忍因子K=k情況下通過(guò)該點(diǎn)的最長(zhǎng)路徑長(zhǎng)度。需要注意的是，缺失路徑長(zhǎng)度k表示多個(gè)連續(xù)節(jié)點(diǎn)的缺失，而不是某一條路徑中多處缺失節(jié)點(diǎn)的個(gè)數(shù)之和，如圖6c中缺失的2個(gè)節(jié)點(diǎn)必須連在一起。類(lèi)似完備路徑開(kāi)運(yùn)算的過(guò)程，不完備路徑開(kāi)同樣是先左向逐列計(jì)算λ-[p,k1]，再右向逐列計(jì)算λ+[p,k2]。當(dāng)b[p]=真時(shí)，k1+k2=k；當(dāng)b[p]=假時(shí)，k1+k2+1=k。如果經(jīng)過(guò)該節(jié)點(diǎn)的路徑在這里存在長(zhǎng)度為k的間斷，需要利用上一個(gè)間斷長(zhǎng)度為k-1的符合臨近關(guān)系的最長(zhǎng)路徑來(lái)輔助計(jì)算，詳細(xì)公式在下面給出。

若b[p]=假：

若b[p]=真：

根據(jù)公式(6)～(10)對(duì)每個(gè)節(jié)點(diǎn)p(p1,p2)遞歸計(jì)算λ-[p,k]和λ+[p,k]。

2.3 不完備路徑形態(tài)學(xué)對(duì)模擬含噪聲電成像數(shù)據(jù)的處理分析

為了考查不完備路徑形態(tài)學(xué)與Otsu算法對(duì)模擬含噪聲電成像數(shù)據(jù)的適應(yīng)能力和處理效果，設(shè)計(jì)了含噪且溶蝕孔洞與多條相交裂縫共生的縫洞儲(chǔ)層電成像響應(yīng)模型。下文將測(cè)試不同參數(shù)對(duì)去噪和縫洞分離的效果。

如圖7a所示，設(shè)計(jì)兩條相交裂縫和兩個(gè)不同尺度溶蝕孔洞的組合模型，大小是552×314像素。圖中可以看到兩個(gè)正弦裂縫均被斷成了多個(gè)部分，還有兩個(gè)線(xiàn)段表示直線(xiàn)型裂縫或者雁陣型誘導(dǎo)縫，兩個(gè)橢圓表示不同直徑和縱橫比的孔洞，并給圖像加上方差為0.005的二維高斯噪聲。圖7b針對(duì)圖7a進(jìn)行了彩色圖轉(zhuǎn)換成灰度圖的操作，圖7c利用了Otsu算法對(duì)圖7b進(jìn)行了二值化分割，其中白色部分像素值是1，黑色部分是0。可以看到圖像的背景都被去掉了，僅剩裂縫、孔洞和噪聲，正弦型裂縫中存在間斷的最長(zhǎng)距離小于15個(gè)像素長(zhǎng)度。

圖8是利用不同的不完備路徑形態(tài)學(xué)參數(shù)組合對(duì)圖7c兩條裂縫圖像提取的效果。其中，圖8a是L=90，K=0時(shí)的提取效果，即完備路徑形態(tài)學(xué)下L=90的效果。此時(shí)因?yàn)楸尘霸朦c(diǎn)的路徑長(zhǎng)度都小于90，所以噪點(diǎn)被去除了，同時(shí)被去除的還有一條直線(xiàn)型小裂縫和兩個(gè)橢圓狀孔洞。圖8b是L=90，K=15時(shí)的提取效果，此時(shí)圖7c中一個(gè)較大的孔洞較長(zhǎng)的一部分再次出現(xiàn)，另一個(gè)直線(xiàn)型小裂縫亦然。圖8c是L=150，K=0時(shí)的提取效果。此時(shí)小的直線(xiàn)裂縫和孔洞都被去掉了，但是因?yàn)榇藭r(shí)是完備路徑形態(tài)學(xué)提取，所以?xún)蓷l正弦型裂縫中間斷型存在的部分就被略去了，這無(wú)法滿(mǎn)足要求。

圖7 兩條正弦型裂縫示意圖Fig.7 A schematic diagram of two sinusoidal curves

圖8 不完備路徑形態(tài)學(xué)算法不同參數(shù)組合的提取效果Fig.8 Different extraction results of combinations of parameters by incomplete path opening operation

由于正弦型裂縫中存在間斷的最長(zhǎng)距離小于15個(gè)像素長(zhǎng)度，所以選擇容忍閾值15。圖8d是L=150，K=15時(shí)的提取效果。在增加了容忍閾值的設(shè)置后，兩條正弦型裂縫被完整提取出來(lái)，對(duì)比圖8c，可以看到不完備路徑算法的優(yōu)勢(shì)。

繼續(xù)增大L，圖8e是L=200，K=0時(shí)的提取效果，此時(shí)兩條正弦型裂縫中的大部分被略去了。圖8f是L=200，K=15時(shí)的提取效果，增加容忍閾值的設(shè)置，正弦型裂縫中仍然有一小部分缺失，這表明L設(shè)置過(guò)大。圖8g～h是L=280時(shí)的提取效果，作為L(zhǎng)設(shè)置過(guò)大的極端情況，無(wú)論K值是否是0，裂縫都無(wú)法被完整提取。綜上，可以看出使用不完備路徑形態(tài)學(xué)在參數(shù)組合L=150，K=15的時(shí)候，兩條正弦型裂縫的提取效果最好。

3 密度聚類(lèi)算法以及正弦函數(shù)族算法對(duì)相交裂縫的分離

在實(shí)際的測(cè)井資料中多是如圖7a那樣兩條甚至多條交織在一起的正弦曲線(xiàn)，為了把這兩條正弦曲線(xiàn)分開(kāi)，本文在使用不完備路徑形態(tài)學(xué)算法的基礎(chǔ)上結(jié)合密度聚類(lèi)算法以及正弦函數(shù)族匹配的算法來(lái)完成。

3.1 密度聚類(lèi)算法對(duì)多條相交裂縫的提取

圖9a是對(duì)圖8c使用L=150，K=15的路徑形態(tài)學(xué)算子提取的結(jié)果，之后進(jìn)行去噪，采用長(zhǎng)和寬均是4的disk型結(jié)構(gòu)元素對(duì)圖像進(jìn)行開(kāi)運(yùn)算操作，即對(duì)原圖像先腐蝕后膨脹，圖9b是去噪后的結(jié)果，圖中可以看到白色噪點(diǎn)被去除，而兩條正弦曲線(xiàn)幾乎沒(méi)有受到破壞，去噪效果良好。下邊分析對(duì)圖9b進(jìn)行兩個(gè)方向的路徑形態(tài)學(xué)分離。

仍然使用容忍閾值L=150，K=15的不完備路徑形態(tài)學(xué)算子對(duì)圖9b進(jìn)行提取，但不再是如圖4所示的四個(gè)方向的提取，而是選取了其中的兩個(gè)方向，以方便后續(xù)討論。圖10a是算子“下-右下-右”的示意圖，圖10b是用該算子提取的結(jié)果；圖10c是算子“上-右上-右”的示意圖，圖10d是用該算子提取的結(jié)果?？梢钥吹綀D10b和圖10d分別提取出了帶有間斷的裂縫的一部分。雖然在每一條裂縫上還留著與之接近垂直方向的裂縫的“刺”，但它對(duì)下一步的分離效果幾乎沒(méi)有影響。為了完成纏繞裂縫的分離，選用聚類(lèi)方法來(lái)實(shí)現(xiàn)。

圖9 開(kāi)運(yùn)算去噪效果圖Fig.9 A denoising result by open mathematical morphology method

圖10 使用兩個(gè)不同方向的不完備路徑形態(tài)參數(shù)組合提取的結(jié)果Fig.10 Results by using two different directions of incomplete path opening parameter combinations

基于聚類(lèi)的算法是機(jī)器學(xué)習(xí)中涉及對(duì)數(shù)據(jù)進(jìn)行分組的一種算法，通過(guò)分析給定數(shù)據(jù)集在不同域的分布規(guī)律，可以把數(shù)據(jù)分成多組，同一組內(nèi)的數(shù)據(jù)具有相同或類(lèi)似的特征，而不同組之間的數(shù)據(jù)特征不一致。常見(jiàn)的聚類(lèi)算法有K-means聚類(lèi)、Mean-Shift聚類(lèi)、基于密度的帶噪聲的空間聚類(lèi)(DBSCAN)等[29]。本文使用的是DBSCAN算法，它包含兩個(gè)重要參數(shù)：ε(鄰域距離參數(shù))和minPts(鄰域最少樣本點(diǎn)參數(shù))，該算法用這兩個(gè)參數(shù)來(lái)尋找數(shù)據(jù)集內(nèi)各點(diǎn)在某一坐標(biāo)域內(nèi)具有較高密度的數(shù)據(jù)點(diǎn)的組合，并把組合劃歸成不同的簇，具體運(yùn)算關(guān)系可以用下列式子表示：

樣本數(shù)據(jù)集D= {x1,x2,...,xm}。

定義以下概念：

ε-鄰接距離

設(shè)xj∈D，其ε-鄰接距離指的是在數(shù)據(jù)集D中與點(diǎn)xj不遠(yuǎn)于ε的數(shù)據(jù)，用Nε(xj) = {xi∈D|dist(xi,xj)≤ε}表示；

核心有效數(shù)據(jù)點(diǎn)

若xj的鄰域具有不少于minPts個(gè)數(shù)據(jù)點(diǎn)，滿(mǎn)足Nε(xj)≥minPts的條件，那么xj被認(rèn)為是核心有效數(shù)據(jù)點(diǎn)，minPts數(shù)包含了xj自身。

直接觸及關(guān)系

對(duì)xi與xj，若xi是核心有效數(shù)據(jù)點(diǎn)且xj與xi的距離不大于ε-鄰接距離，則稱(chēng)xj由xi構(gòu)成直接觸及關(guān)系；

觸及關(guān)系

對(duì)xi與xj，若存在xk使得xi與xj均與xk構(gòu)成直接觸及關(guān)系，則稱(chēng)xj由xi構(gòu)成觸及關(guān)系；

圖11是一個(gè)簡(jiǎn)單的算法示意圖，圖中minPts=3。圖中各個(gè)圓環(huán)表示ε-鄰接距離，核心有效數(shù)據(jù)點(diǎn)是x1，可以看到x2、x4與x1均構(gòu)成直接觸及關(guān)系，而x3與x1構(gòu)成觸及關(guān)系，因此x3與x4也構(gòu)成觸及關(guān)系。

圖11 DBSCAN算法示意圖Fig.11 A schematic diagram of DBSCAN algorithm.

從上述定義出發(fā)，基于密度的帶噪聲的空間聚類(lèi)的應(yīng)用算法給出了“簇”的定義，即一定區(qū)域內(nèi)最大數(shù)目的多個(gè)互相構(gòu)成觸及關(guān)系的數(shù)據(jù)點(diǎn)連接到一起的集合。若x為核心有效數(shù)據(jù)點(diǎn)，那么找到所有與之有觸及關(guān)系的集合為X={x'∈D|x'由x密度可達(dá)}，X滿(mǎn)足簇的連接性與最大性的性質(zhì)。

按(ε，minPts)參數(shù)重復(fù)上面的過(guò)程，直到數(shù)據(jù)集內(nèi)所有點(diǎn)都被該算法掃描完畢。

對(duì)于有間斷的兩條正弦曲線(xiàn)，選用不同的聚類(lèi)參數(shù)組合會(huì)得到不同的結(jié)果。圖12是在ε=3，minPts=3的情況下獲得的結(jié)果。

圖12b和e分別是用兩條帶有間隔的正弦曲線(xiàn)通過(guò)不完備路徑形態(tài)學(xué)在圖12a和圖12d方向的提取結(jié)果。對(duì)圖12b和e進(jìn)行密度聚類(lèi)，參數(shù)設(shè)定為ε=3，minPts=3，這時(shí)如果3個(gè)及以上的像素點(diǎn)每?jī)蓚€(gè)像素點(diǎn)周?chē)淮笥?，它們都會(huì)被劃入同一個(gè)聚類(lèi)簇，得到的結(jié)果如圖12c和f所示。由于間隔的存在，同一條正弦曲線(xiàn)中間隔都大于3個(gè)像素距離，可以看到圖12c和f均被DBSCAN算法自動(dòng)分成了7種顏色對(duì)應(yīng)的7部分，本應(yīng)該通過(guò)密度聚類(lèi)算法分到一起的正弦曲線(xiàn)段因?yàn)殚g斷的存在而被分割成多個(gè)部分，這沒(méi)有達(dá)到做密度聚類(lèi)的初衷。

圖12 兩條帶有間隔的正弦曲線(xiàn)使用不完備路徑形態(tài)學(xué)提取后密度聚類(lèi)的結(jié)果。(a)不完備路徑形態(tài)算子的提取方向一；(b)使用a獲得的提取結(jié)果；(c)對(duì)圖b的密度聚類(lèi)結(jié)果；(d)不完備路徑形態(tài)算子的提取方向二；(e)使用d獲得的提取結(jié)果；(f)對(duì)圖e的密度聚類(lèi)結(jié)果Fig.12 Results by DBSCAN after using incomplete operation on two sinusoidal curves with several intervals.(a)First direction of incomplete path opening operator; (b)The result by using operator a; (c) The result of DBSCAN on b; (d)Second direction of incomplete path opening operator; (e)The result by using operator d; (f) The result of DBSCAN on e

嘗試使用另一組聚類(lèi)參數(shù)組合。圖13是在如果ε=30，minPts=3的情況下獲得的結(jié)果。

圖13b和e分別是用兩條帶有間隔的正弦曲線(xiàn)通過(guò)不完備路徑形態(tài)學(xué)算法在圖13a和d方向上的提取結(jié)果。對(duì)圖13b和e進(jìn)行密度聚類(lèi)，參數(shù)設(shè)定為ε=30，minPts=3，這時(shí)如果3個(gè)及以上的像素點(diǎn)每?jī)蓚€(gè)像素點(diǎn)周?chē)淮笥?0，它們都會(huì)被劃入同一個(gè)聚類(lèi)簇，得到的結(jié)果如圖13c和f所示。由于同一條正弦曲線(xiàn)中間隔都小于30個(gè)像素距離，所以可以看到圖13c和f被分成了3種顏色對(duì)應(yīng)的3部分，原來(lái)屬于一條正弦曲線(xiàn)的部分被歸為一類(lèi)，因此密度聚類(lèi)的結(jié)果是正確的。

圖13 兩條帶有間隔的正弦曲線(xiàn)使用不完備路徑形態(tài)學(xué)提取后密度聚類(lèi)的結(jié)果(a)不完備路徑形態(tài)算子的提取方向一；(b)使用a獲得的提取結(jié)果；(c)對(duì)圖b的密度聚類(lèi)結(jié)果；(d)不完備路徑形態(tài)算子的提取方向二；(e) 使用d獲得的提取結(jié)果；(f)對(duì)圖e的密度聚類(lèi)結(jié)果Fig.13 Results by DBSCAN after using incomplete operation on two sinusoidal curves with several intervals (a)First direction of incomplete path opening operator; (b)The result by using operator a; (c) The result of DBSCAN on b;(d) Second direction of incomplete path opening operator; (e)The result by using operator d; (f) The result of DBSCAN on e

3.2 基于正弦函數(shù)族的裂縫參數(shù)提取算法

由電導(dǎo)率圖像提取和分離裂縫與溶蝕孔洞后，下一個(gè)關(guān)鍵任務(wù)是裂縫和溶蝕孔洞參數(shù)的定量評(píng)價(jià)。對(duì)于規(guī)則裂縫與井孔相交在電導(dǎo)率圖像上呈一條有初始相位角的正弦曲線(xiàn)，非規(guī)則裂縫則呈現(xiàn)變形的曲線(xiàn)。溶蝕孔洞則在電導(dǎo)率圖像上大多呈橢圓形，也有不規(guī)則形狀。因此，根據(jù)分離的裂縫與溶蝕孔洞分別擬合，求取裂縫開(kāi)度、方位傾角以及溶蝕孔洞的長(zhǎng)短軸、縱橫比和面孔率等參數(shù)。下面針對(duì)不同縫洞的組合形狀，分析相應(yīng)的擬合方法。

Hough變換可以針對(duì)電成像測(cè)井圖像中的正弦型裂縫進(jìn)行提取，方法是把圖像變換到參數(shù)域后計(jì)算出裂縫的振幅和相位參數(shù)[30]。在轉(zhuǎn)換的過(guò)程中，如果圖像中的裂縫曲線(xiàn)歸屬于同一個(gè)正弦曲線(xiàn)，那么在參數(shù)域會(huì)疊加在一點(diǎn)上。在測(cè)井解釋中，Hough變換參數(shù)域是作為一個(gè)投票器來(lái)對(duì)圖像中的各點(diǎn)變換后的信息進(jìn)行投票。在投票結(jié)束后，參數(shù)域中較高的那一點(diǎn)會(huì)分別顯示出圖像中對(duì)應(yīng)曲線(xiàn)的振幅和相位。

Hough變換的優(yōu)勢(shì)是能快速提取微電導(dǎo)率測(cè)井成像中的裂縫信息，但它也有缺點(diǎn)。第一，噪音的存在影響了正弦曲線(xiàn)的提??；第二，對(duì)高角度裂縫擬合容易出現(xiàn)較大誤差。

從前面的描述中，知道電導(dǎo)率測(cè)井成像中的裂縫信息可以用正弦函數(shù)來(lái)表示，因此可以嘗試使用模式識(shí)別的方式來(lái)提取裂縫。在本文中，構(gòu)建了一個(gè)正弦函數(shù)族來(lái)進(jìn)行實(shí)際圖像的裂縫提取。所有接近正弦型的裂縫理論上都可以用正弦函數(shù)族來(lái)進(jìn)行擬合。擬合的原理為2-D相關(guān)算法[11]：

其中，r表示2-D相關(guān)系數(shù)，A表示實(shí)際的裂縫成像數(shù)據(jù)，B是正弦函數(shù)族。和分別表示A和B的平均值。

為了對(duì)比Hough變換和正弦函數(shù)族擬合提取裂縫的優(yōu)勢(shì)，這里分別考慮低角度和高角度的裂縫(水平裂縫和豎直裂縫也可以按照類(lèi)似的方法處理)。正弦函數(shù)族的基本參數(shù)設(shè)置如下：傾角區(qū)間設(shè)為-75o～75o，相位區(qū)間設(shè)為0o～360o。在井徑已知的情況下，正弦函數(shù)的振幅可以通過(guò)傾角計(jì)算出來(lái)，深度可以通過(guò)實(shí)際的FMI成像資料獲得?；谶@樣的設(shè)定，正弦函數(shù)族就可以被建立起來(lái)。在實(shí)際的計(jì)算過(guò)程中，傾角區(qū)間和相位區(qū)間的步長(zhǎng)應(yīng)該根據(jù)測(cè)井?dāng)?shù)據(jù)的大小和粗細(xì)等情況斟酌，因?yàn)楫?dāng)單位步長(zhǎng)太小時(shí)，正弦函數(shù)族的計(jì)算代價(jià)會(huì)非常高；當(dāng)單位步長(zhǎng)太大時(shí)，最后的匹配計(jì)算精度會(huì)變差，因此應(yīng)尋找一個(gè)計(jì)算代價(jià)與精度的平衡點(diǎn)。

圖14a展示了三條殘縫，本節(jié)試圖用不同的方法來(lái)重建原始的正弦曲線(xiàn)。從圖14b～c重建結(jié)果和圖14e輪盤(pán)的對(duì)比結(jié)果可以看到正弦函數(shù)族匹配的結(jié)果精度更高，更接近殘縫原始的振幅，而使用Hough變換重建的正弦曲線(xiàn)振幅超過(guò)了原始模型的振幅，同時(shí)從圖14d看到Hough變換域也并不是全部收斂到某一個(gè)“亮點(diǎn)”上。因此當(dāng)原始?xì)埧p的信息不夠多時(shí)，使用Hough變換提取的結(jié)果精度較低。

圖15是使用正弦函數(shù)族對(duì)圖13b和13e聚類(lèi)后的一共6部分分別進(jìn)行重建恢復(fù)的結(jié)果。文中使用了32GB的內(nèi)存RAM來(lái)進(jìn)行重建工作，匹配精度較高。可以看到各部分的重建裂縫與模型最初的正弦型很接近，但有一些振幅出現(xiàn)了偏差，這是由于殘縫提供的信息太少導(dǎo)致的。該匹配過(guò)程可以在更高內(nèi)存的工作站上使用相位、振幅和深度間隔步長(zhǎng)更小的正弦函數(shù)族來(lái)進(jìn)行匹配，精度更高，但會(huì)花費(fèi)更多的時(shí)間和內(nèi)存計(jì)算成本。

圖15 使用正弦函數(shù)族對(duì)圖13b和圖13e聚類(lèi)后的各部分分別進(jìn)行重建恢復(fù)的結(jié)果Fig.15 The reconstruction results of every part after DBSCAN in Fig.13b and Fig.13e

接下來(lái)對(duì)恢復(fù)重建的正弦函數(shù)進(jìn)行歸類(lèi)，假定恢復(fù)后的兩個(gè)正弦函數(shù)的相位差及深度差滿(mǎn)足：

(1)φi-φj＜π/4,1 ≤i≤ 6 ,1 ≤j≤ 6且i≠j；

(2)Depthi-Depthj＜ 0.2,1 ≤i≤ 6 ,1 ≤j≤6且i≠j。

即可以認(rèn)定這兩個(gè)正弦函數(shù)屬于同一條裂縫，那么它們對(duì)應(yīng)的原來(lái)的殘縫也隸屬于同一條裂縫。這是因?yàn)樵趯?shí)際地層中，深度間隔小于0.2 m的任意兩條殘縫的相位如果小于π/4，那通常是由于地層擠壓或拉伸造成同一個(gè)地層發(fā)生了錯(cuò)位導(dǎo)致的。基于該思想，把圖15a～f中所有隸屬于一條裂縫的殘縫疊加到一起，可以得到如圖16的結(jié)果，圖中可以看到，兩條正弦曲線(xiàn)被較好地分離開(kāi)來(lái)。

圖16 (a)使用形態(tài)學(xué)算法去噪后的電導(dǎo)率測(cè)井圖像；(b)從a中分離出來(lái)的第一條裂縫；(c)從a中分離出來(lái)的第二條裂縫Fig.16 (a)The electric imaging logging image after denoising by mathematical morphology; (b)The first separated fracture from a; (c) The second separated fracture from a

4 基于二階中心矩的溶蝕孔洞擬合與參數(shù)提取算法

上一節(jié)對(duì)碳酸鹽巖縫洞型儲(chǔ)層的電導(dǎo)率圖像利用不完備路徑形態(tài)學(xué)和正弦函數(shù)族算法實(shí)現(xiàn)了纏繞裂縫的分離，但僅分離了裂縫，為了后續(xù)測(cè)井處理時(shí)計(jì)算滲透率、飽和度，同深度段的溶蝕孔洞信息需要被進(jìn)一步計(jì)算，例如孔隙縱橫比、準(zhǔn)確的深度分布等。本節(jié)將圍繞該主題，使用基于二階中心矩橢圓的孔洞擬合和參數(shù)提取算法得到溶蝕孔洞在空間上的分布及特征。本節(jié)先對(duì)二階矩?cái)M合算法進(jìn)行公式推導(dǎo)，然后展示模型數(shù)據(jù)的計(jì)算結(jié)果。

4.1 二階中心矩溶蝕孔洞橢圓擬合算法

傳統(tǒng)文獻(xiàn)關(guān)于橢圓擬合的思路分為兩個(gè)路線(xiàn)：最小平方擬合[31]和聚類(lèi)算法[32]。與傳統(tǒng)的基于散點(diǎn)圖的橢圓擬合不同，孔洞在圖像上是一塊區(qū)域，不是多個(gè)邊緣散落點(diǎn)的組合，提取更為復(fù)雜，涉及的參數(shù)也更多，且孔洞結(jié)構(gòu)的復(fù)雜性又增大了其參數(shù)提取的難度。為了便于實(shí)際應(yīng)用，文中用不同尺度與偏心率的二階矩橢圓形函數(shù)來(lái)擬合孔洞。

一般來(lái)說(shuō)，可以對(duì)一個(gè)實(shí)際物體在任意維度上通過(guò)定義一個(gè)矩來(lái)獲取它的有效信息。對(duì)實(shí)際物體可以認(rèn)為它有一個(gè)總的質(zhì)心，這對(duì)應(yīng)的是一階矩。在圖像中，二階矩可以被用來(lái)決定與該物體對(duì)應(yīng)的橢圓，從中可以找出長(zhǎng)軸和短軸的長(zhǎng)度及方向。

圖像的矩通過(guò)圖像中像素的密集程度來(lái)定義。在一個(gè)像素密度為I(i,j)的灰度圖中，簡(jiǎn)單的二維(p+q)矩Mp,q按如下給出：

在下面的討論中，灰度圖像矩的計(jì)算對(duì)圖像的局部對(duì)比度有較強(qiáng)依賴(lài)，因此處理起來(lái)更加困難。所以在進(jìn)行計(jì)算前，先把圖像變成二值圖，即圖像中物體的像素值變成1，而背景值變成0。對(duì)于一個(gè)二值圖，(p,q)矩Mp,q可以被定義如下：

零階矩M0,0給出了物體的總像素點(diǎn)數(shù)，可以理解為其面積。一階矩M1,0和M0,1可以理解成該物體分別在水平和垂直方向的重心：

從這里開(kāi)始，算法就變得復(fù)雜起來(lái)，從二階矩M2,0，M1,0和M0,2中提取等效橢圓的參數(shù)不是最簡(jiǎn)便的方法，因?yàn)楸仨毷褂枚A中心矩，這與計(jì)算物體的剛度二階矩不同。它們的表達(dá)式如下：

經(jīng)過(guò)上述計(jì)算，二值圖像中物體的協(xié)方差矩陣變成

該協(xié)方差矩陣的特征向量對(duì)應(yīng)了物體等效橢圓的長(zhǎng)軸和短軸方向。如圖17所示，等效橢圓度的長(zhǎng)軸與x軸的夾角θ、長(zhǎng)軸長(zhǎng)度l和短軸長(zhǎng)度w的最終形式為：

圖17 擬合后的橢圓示意圖Fig.17 The schematic diagram of ellipse approximation

確定橢圓中心位置后，還需統(tǒng)計(jì)目標(biāo)區(qū)橢圓的面積，以便定量計(jì)算面孔率。橢圓面積計(jì)算公式為：

橢圓對(duì)應(yīng)的三維空間的橢球體積計(jì)算公式為：

其中，a、b分別對(duì)應(yīng)了橢圓的長(zhǎng)、短軸長(zhǎng)度的一半。

橢圓的縱橫比計(jì)算公式為：

4.2 模擬數(shù)據(jù)測(cè)試

為了測(cè)試前文橢圓擬合算法的有效性，本文選取了散列的硬幣圖像進(jìn)行二階矩橢圓擬合，并給出了橢圓的中心點(diǎn)，效果如圖18所示，可以看到二值圖中每枚硬幣使用二階中心矩橢圓擬合后找到了對(duì)應(yīng)的橢圓及中心。

圖18 (a)原始硬幣圖像；(b)用二階中心矩?cái)M合各硬幣的橢圓Fig.18 (a)The original images of coins; (b)Enclosing ellipses of these coins based on the second-order centric moments method

5 測(cè)井?dāng)?shù)據(jù)計(jì)算與縫洞儲(chǔ)層評(píng)價(jià)應(yīng)用

5.1 不完備路徑形態(tài)學(xué)和正弦函數(shù)族算法在提取裂縫參數(shù)過(guò)程中的應(yīng)用

圖19是對(duì)正弦函數(shù)族方法與改進(jìn)的Hough變換方法在實(shí)際電成像測(cè)井資料中的效果進(jìn)行比較。

圖19a中包含三條斷裂的裂縫，其中上下兩條縫發(fā)育較好，但中間的裂縫被巖層的基質(zhì)切斷。圖19b是使用不完備路徑形態(tài)學(xué)算子L=40，K=5對(duì)二值化后的圖像提取的結(jié)果，基質(zhì)、孔洞等被剔除，只剩下三條較大的裂縫。圖19c是用改進(jìn)的Hough變換重建的裂縫，圖中上下兩條裂縫被較好重建出來(lái)，但中間的裂縫由于原來(lái)電導(dǎo)率圖像中只有部分信息，因此用Hough變換重建出現(xiàn)了錯(cuò)誤。圖19d是用正弦函數(shù)族進(jìn)行模式識(shí)別的結(jié)果，可以看到三條裂縫都獲得了較好的恢復(fù)效果，中間殘縫的重建結(jié)果與原電導(dǎo)率圖像吻合。圖19e為實(shí)際電成像測(cè)井資料用兩種方法得到的參數(shù)的羅盤(pán)圖，可以清晰看到正弦函數(shù)族得到的參數(shù)和模型參數(shù)基本吻合，而Hough變換則存在較大誤差。因此可以斷定，當(dāng)初始模型有效信息較少，裂縫不完整的情況下，正弦函數(shù)族提取的裂縫參數(shù)比Hough變換得到的結(jié)果更準(zhǔn)確。

圖19 (a)原始的電導(dǎo)率測(cè)井圖像；(b)使用不完備路徑形態(tài)學(xué)算子L=40，K=5提取的空白條帶插值后的結(jié)果；(c)使用Hough變換重建的裂縫；(d)使用正弦函數(shù)族重建的裂縫；(e)Hough變換和正弦函數(shù)族重建的結(jié)果對(duì)比Fig.19 (a) The original electric imaging logging data; (b) The extraction by incomplete path opening operation with L=40 and K=5 of the blank area interpolated data; (c) The reconstruction fractures using Hough transform; (d) The reconstruction fractures using the cluster of sinusoid; (e) The comparison of Hough transform and the cluster of sinusoid results

接下來(lái)使用正弦函數(shù)族與不完備路徑形態(tài)學(xué)方法結(jié)合，對(duì)實(shí)際的微電導(dǎo)率成像中存在纏繞現(xiàn)象的裂縫進(jìn)行分離。

首先使用不完備路徑形態(tài)學(xué)路徑算子和容忍算子來(lái)進(jìn)行裂縫的提取。

圖20a是原圖，大小是735×441，圖中可以看到有兩條相交裂縫，在裂縫周?chē)植贾鞣N溶蝕孔洞。圖20b是原圖先使用Otsu算法進(jìn)行二值化之后用開(kāi)運(yùn)算處理的結(jié)果。圖20c是L=30K=5的提取效果，可以看出因?yàn)槁窂介L(zhǎng)度因子較小，且容忍因子K非0，所以一部分溶蝕孔洞仍然得以保留。圖20d是L=56，K=0的提取效果，可以看出路徑長(zhǎng)度因子變大使一部分溶蝕孔洞被剔除，但構(gòu)成裂縫的長(zhǎng)度較短的部分由于容忍因子為0，沒(méi)有被提取出來(lái)。圖20e是L=56，K=5時(shí)的提取效果，可以看出連通圖20d遺漏部分的兩條裂縫被較好提取出來(lái)。圖20f是L=56，K=15時(shí)的提取效果，由于容忍因子K過(guò)大，導(dǎo)致裂縫周邊的孔洞信息重新出現(xiàn)。綜上，提取目標(biāo)裂縫最好的參數(shù)是圖20c中的L=56，K=5。

圖20 (a)原始電阻率測(cè)井圖像；(b)使用Otsu算法二值化后的圖像；(c)L=30 K=5的不完備路徑形態(tài)學(xué)提取結(jié)果；(d)L=56 K=0的不完備路徑形態(tài)學(xué)提取結(jié)果；(e)L=56 K=5的不完備路徑形態(tài)學(xué)提取結(jié)果；(f)L=56 K=15的不完備路徑形態(tài)學(xué)提取結(jié)果Fig.20 (a) The original electric imaging logging image; (b)The binary image after using Otsu method; (c) The extraction result by incomplete path opening operation with L=30 and K=5; (d) The extraction result by incomplete path opening operation with L=56 and K=0; (e) The extraction result by incomplete path opening operation with L=56 and K=5; (f)The extraction result by incomplete path opening operation with L=56 and K=15

下面分別使用“下-右下-右”方向的路徑形態(tài)學(xué)算子和“上-右上-右”方向的路徑形態(tài)學(xué)算子對(duì)圖20c采用L=56，K=5的參數(shù)進(jìn)行提取，得到如圖21(b)和21(d)的結(jié)果，然后對(duì)這兩圖使用ε=30，minPts=3的密度聚類(lèi)參數(shù)組合，可以獲得圖21c和21e的聚類(lèi)分組，圖中可以看到，對(duì)用兩個(gè)方向的路徑形態(tài)學(xué)算子提取的數(shù)據(jù)分別進(jìn)行聚類(lèi)后，數(shù)據(jù)被分成了4部分，每部分用不同的顏色表示，其中兩圖右下角的部分是無(wú)效數(shù)據(jù)。

對(duì)圖21c和21e有效數(shù)據(jù)使用正弦函數(shù)族進(jìn)行提取和擬合。每部分對(duì)應(yīng)的正弦曲線(xiàn)如圖22所示，且每部分對(duì)應(yīng)的正弦曲線(xiàn)的參數(shù)列在表1中。

表1 6條擬合得到的正弦曲線(xiàn)的參數(shù)Table 1 The parameters of 6 sine curves.

圖21 兩條裂縫使用不完備路徑形態(tài)學(xué)提取后密度聚類(lèi)的結(jié)果 (a)不完備路徑形態(tài)算子的提取方向一；(b)使用a獲得的提取結(jié)果；(c)圖b的密度聚類(lèi)結(jié)果；(d)不完備路徑形態(tài)算子的提取方向二；(e) 使用d獲得的提取結(jié)果；(f)圖e的密度聚類(lèi)結(jié)果Fig; 21 Results by DBSCAN after using incomplete operation on two fractures (a)First direction of incomplete path opening operator; (b)The result by using operator a;(c) The result of DBSCAN on b; (d) Second direction of incomplete path opening operator; (e)The result by using operator d; (f) The result of DBSCAN on e

圖22 使用正弦函數(shù)族對(duì)圖21b和圖21e聚類(lèi)后的各部分分別進(jìn)行重建恢復(fù)的結(jié)果Fig.22 The reconstruction results of every part after DBSCAN in Fig.21b and Fig.21e

上文中在對(duì)恢復(fù)重建的正弦函數(shù)進(jìn)行歸類(lèi)時(shí)，假定兩個(gè)正弦函數(shù)的相位差和深度差滿(mǎn)足一定條件即可以認(rèn)定這兩個(gè)正弦函數(shù)屬于同一條裂縫，那么它們對(duì)應(yīng)的原來(lái)的殘縫也隸屬于同一條裂縫?；谠摷僭O(shè)，把六部分殘縫進(jìn)行疊加，獲得圖23的結(jié)果。圖中可以看到，兩條裂縫被較好地分離開(kāi)來(lái)，且兩條裂縫的深度、振幅及相位如表2所示。

表2 兩條分離后的裂縫參數(shù)Table 2 The parameters of 20 separated fractures.

圖23 (a)使用不完備路徑形態(tài)學(xué)提取結(jié)果；(b)從a中分離出來(lái)的第一條裂縫；(c)從a中分離出來(lái)的第二條裂縫Fig.23 (a) The extraction result by incomplete path opening operation; (b) The first separated fracture from a; (c) The second separated fracture from a

5.2 二階矩橢圓擬合算法在提取裂縫參數(shù)過(guò)程中的應(yīng)用

對(duì)上節(jié)中圖21a的微電導(dǎo)率測(cè)井圖像，使用不完備路徑形態(tài)學(xué)、密度聚類(lèi)和正弦函數(shù)族算法分離出裂縫后，對(duì)剩余的孔洞部分使用二階矩橢圓擬合算法，效果分別如圖24所示。

圖24 (a)對(duì)圖19a使用Otsu算法二值化后的圖像；(b)從a中分離出來(lái)的裂縫；(c)從a中分離出來(lái)的孔洞；(d)孔洞基于二階矩橢圓擬合的結(jié)果Fig.24 (a) The binary image to Fig.19 after using Otsu method; (b) Fracture separated from a; (c) Vugs separated from a; (d) The equivalent ellipses of the vugs based on the second-order moments method

對(duì)圖24中的縫洞參數(shù)進(jìn)行縱橫比和面積(像素意義)統(tǒng)計(jì)，得到的結(jié)果顯示在表3中。

表3 圖23孔洞參數(shù)統(tǒng)計(jì)Table 3 The parameter statistics of vugs in Fig.23

6 結(jié)論

本文提出了一種基于不完備路徑形態(tài)學(xué)、正弦函數(shù)族和二階矩橢圓擬合算法由電成像測(cè)井電導(dǎo)率圖像自動(dòng)提取和分離裂縫和溶蝕孔洞的方法，通過(guò)縫洞模型的模擬數(shù)據(jù)和實(shí)測(cè)數(shù)據(jù)初步驗(yàn)證了方法的有效性，并得到了如下認(rèn)識(shí)：

(1) 路徑形態(tài)學(xué)算法應(yīng)用于電成像測(cè)井電導(dǎo)率圖像的去噪、縫洞邊緣點(diǎn)檢測(cè)與分離，可有效地提取裂縫和溶蝕孔洞的信息，實(shí)現(xiàn)縫洞自動(dòng)識(shí)別、分離和定量表征。

(2) 文中研究的不完備路徑形態(tài)學(xué)處理算法完全由數(shù)據(jù)驅(qū)動(dòng)，算法自身具有去噪能力，可以通過(guò)掃描路徑算子長(zhǎng)度和容忍度參數(shù)，構(gòu)造適合處理數(shù)據(jù)體噪聲特點(diǎn)的路徑算子，從而可直接對(duì)有噪聲的成像測(cè)井資料進(jìn)行處理，檢測(cè)縫洞發(fā)育帶；對(duì)孔洞區(qū)域進(jìn)行二階矩橢圓擬合的算法可以獲得孔洞在像素意義下的面積和縱橫比，可以為后續(xù)估計(jì)對(duì)應(yīng)深度段的孔洞的孔隙度和連通性、滲透率提供支撐，且有助于估計(jì)目的層段的含油儲(chǔ)量。

(3) 文中算法在成像測(cè)井資料處理中的應(yīng)用仍處于探索階段，例如，不同形態(tài)和尺度的縫洞，在邊緣檢測(cè)時(shí)要用不同的路徑算子尺度和容忍度。因此，如何在成像測(cè)井資料處理中根據(jù)地質(zhì)特征自適應(yīng)確定合理的算子長(zhǎng)度與容忍度參數(shù)，仍是一個(gè)需要深入研究的課題。