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

1 中國石油大學(xué)(北京)油氣資源與探測國家重點實驗室，北京 102249

2 中國石油大學(xué)(北京)CNPC物探重點實驗室，北京 102249

0 引言

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

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

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

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

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

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

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

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

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

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

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

圖2所示的電成像測井電導(dǎo)率圖像中暗色或黑色斑塊或條帶顯示了可能存在的裂縫或溶蝕孔洞，如何從這些圖像中識別和分離出該類電導(dǎo)率異常是電成像測井縫洞儲層識別與評價的首要任務(wù)。下文首先簡單介紹從電導(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 基于二值圖路徑開的裂縫識別

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

圖3 鄰接關(guān)系示意圖，b1、b2、b3是a的三個前驅(qū)節(jié)點，a1、a2、a3是b的三個后繼節(jié)點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.

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

對于二維臨近關(guān)系節(jié)點集合E，在每個像素點p存儲兩個值：起始點為p且方向為右的最長路徑的長度λ+[p](該長度不包含點p自身)，以及起始點為p且方向為左的最長路徑的長度λ-[p](該長度不包含點 自身)，其計算公式如下：

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

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

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

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

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

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

2.2 基于二維節(jié)點圖的不完備路徑開的裂縫識別運算

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

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

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

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

若b[p]=假：

若b[p]=真：

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

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

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

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

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

圖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

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

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

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

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

3.1 密度聚類算法對多條相交裂縫的提取

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

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

圖9 開運算去噪效果圖Fig.9 A denoising result by open mathematical morphology method

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

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

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

定義以下概念：

ε-鄰接距離

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

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

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

直接觸及關(guān)系

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

觸及關(guān)系

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

圖11是一個簡單的算法示意圖，圖中minPts=3。圖中各個圓環(huán)表示ε-鄰接距離，核心有效數(shù)據(jù)點是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ā)，基于密度的帶噪聲的空間聚類的應(yīng)用算法給出了“簇”的定義，即一定區(qū)域內(nèi)最大數(shù)目的多個互相構(gòu)成觸及關(guān)系的數(shù)據(jù)點連接到一起的集合。若x為核心有效數(shù)據(jù)點，那么找到所有與之有觸及關(guān)系的集合為X={x'∈D|x'由x密度可達(dá)}，X滿足簇的連接性與最大性的性質(zhì)。

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

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

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

圖12 兩條帶有間隔的正弦曲線使用不完備路徑形態(tài)學(xué)提取后密度聚類的結(jié)果。(a)不完備路徑形態(tài)算子的提取方向一；(b)使用a獲得的提取結(jié)果；(c)對圖b的密度聚類結(jié)果；(d)不完備路徑形態(tài)算子的提取方向二；(e)使用d獲得的提取結(jié)果；(f)對圖e的密度聚類結(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

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

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

圖13 兩條帶有間隔的正弦曲線使用不完備路徑形態(tài)學(xué)提取后密度聚類的結(jié)果(a)不完備路徑形態(tài)算子的提取方向一；(b)使用a獲得的提取結(jié)果；(c)對圖b的密度聚類結(jié)果；(d)不完備路徑形態(tài)算子的提取方向二；(e) 使用d獲得的提取結(jié)果；(f)對圖e的密度聚類結(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)率圖像提取和分離裂縫與溶蝕孔洞后，下一個關(guān)鍵任務(wù)是裂縫和溶蝕孔洞參數(shù)的定量評價。對于規(guī)則裂縫與井孔相交在電導(dǎo)率圖像上呈一條有初始相位角的正弦曲線，非規(guī)則裂縫則呈現(xiàn)變形的曲線。溶蝕孔洞則在電導(dǎo)率圖像上大多呈橢圓形，也有不規(guī)則形狀。因此，根據(jù)分離的裂縫與溶蝕孔洞分別擬合，求取裂縫開度、方位傾角以及溶蝕孔洞的長短軸、縱橫比和面孔率等參數(shù)。下面針對不同縫洞的組合形狀，分析相應(yīng)的擬合方法。

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

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

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

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

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

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

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

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

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

(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)定這兩個正弦函數(shù)屬于同一條裂縫，那么它們對應(yīng)的原來的殘縫也隸屬于同一條裂縫。這是因為在實際地層中，深度間隔小于0.2 m的任意兩條殘縫的相位如果小于π/4，那通常是由于地層擠壓或拉伸造成同一個地層發(fā)生了錯位導(dǎo)致的。基于該思想，把圖15a～f中所有隸屬于一條裂縫的殘縫疊加到一起，可以得到如圖16的結(jié)果，圖中可以看到，兩條正弦曲線被較好地分離開來。

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

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

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

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

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

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

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

從這里開始，算法就變得復(fù)雜起來，從二階矩M2,0，M1,0和M0,2中提取等效橢圓的參數(shù)不是最簡便的方法，因為必須使用二階中心矩，這與計算物體的剛度二階矩不同。它們的表達(dá)式如下：

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

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

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

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

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

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

橢圓的縱橫比計算公式為：

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

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

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

5 測井?dāng)?shù)據(jù)計算與縫洞儲層評價應(yīng)用

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

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

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

圖19 (a)原始的電導(dǎo)率測井圖像；(b)使用不完備路徑形態(tài)學(xué)算子L=40，K=5提取的空白條帶插值后的結(jié)果；(c)使用Hough變換重建的裂縫；(d)使用正弦函數(shù)族重建的裂縫；(e)Hough變換和正弦函數(shù)族重建的結(jié)果對比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

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

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

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

圖20 (a)原始電阻率測井圖像；(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é)算子對圖20c采用L=56，K=5的參數(shù)進(jìn)行提取，得到如圖21(b)和21(d)的結(jié)果，然后對這兩圖使用ε=30，minPts=3的密度聚類參數(shù)組合，可以獲得圖21c和21e的聚類分組，圖中可以看到，對用兩個方向的路徑形態(tài)學(xué)算子提取的數(shù)據(jù)分別進(jìn)行聚類后，數(shù)據(jù)被分成了4部分，每部分用不同的顏色表示，其中兩圖右下角的部分是無效數(shù)據(jù)。

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

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

圖21 兩條裂縫使用不完備路徑形態(tài)學(xué)提取后密度聚類的結(jié)果 (a)不完備路徑形態(tài)算子的提取方向一；(b)使用a獲得的提取結(jié)果；(c)圖b的密度聚類結(jié)果；(d)不完備路徑形態(tài)算子的提取方向二；(e) 使用d獲得的提取結(jié)果；(f)圖e的密度聚類結(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ù)族對圖21b和圖21e聚類后的各部分分別進(jìn)行重建恢復(fù)的結(jié)果Fig.22 The reconstruction results of every part after DBSCAN in Fig.21b and Fig.21e

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

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

圖23 (a)使用不完備路徑形態(tài)學(xué)提取結(jié)果；(b)從a中分離出來的第一條裂縫；(c)從a中分離出來的第二條裂縫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ù)過程中的應(yīng)用

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

圖24 (a)對圖19a使用Otsu算法二值化后的圖像；(b)從a中分離出來的裂縫；(c)從a中分離出來的孔洞；(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

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

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

6 結(jié)論

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

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

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

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