袁 陽， 侯國強(qiáng)， 馬濤濤， 周 輝， 董漢強(qiáng)， 潘 濤

(1.中國石油集團(tuán)測井有限公司吐哈分公司，新疆維吾爾自治區(qū) 吐魯番 838200；2.延長油田油氣勘探公司 石油勘探開發(fā)部，陜西 延安 716000)

砂礫巖油藏埋藏深，儲集性能非均質(zhì)性強(qiáng)，復(fù)雜的巖性導(dǎo)致砂礫巖地層的油氣儲集空間類型變化多樣。復(fù)雜的巖石骨架削弱了儲集流體的常規(guī)測井響應(yīng)，造成常規(guī)測井信息識別、劃分砂礫巖巖性、評價儲層難度很大。

由于電成像測井資料分辨率高，可得到砂礫巖地層中礫巖顆粒形狀、大小、空間位置及砂巖和泥質(zhì)條帶的直觀描述，已經(jīng)應(yīng)用到砂礫巖儲層的巖性識別研究中。但以往的研究主要是巖性的定性識別，表現(xiàn)為基于電成像圖像的圖像-巖性識別模式，需要解釋人員根據(jù)模式人工做出巖性的判斷、劃分，沒有實現(xiàn)計算機(jī)巖性自動定量識別[1-4]。目前，文獻(xiàn)[5]采用數(shù)字圖像處理技術(shù)對砂礫巖圖像進(jìn)行特征提取，再基于人工神經(jīng)網(wǎng)絡(luò)自動識別巖性；文獻(xiàn)[6]提取了電成像圖像中的均值、方差、分形維數(shù)、包絡(luò)和非均質(zhì)5條特征曲線，分別利用相關(guān)性和神經(jīng)網(wǎng)絡(luò)兩種方法進(jìn)行巖性自動識別；文獻(xiàn)[7]通過圖像灰度轉(zhuǎn)化，構(gòu)建灰度共生矩陣，計算樣本的對比度、相關(guān)度、熵、均勻度和能量5個特征，利用貝葉斯判別分析法對巖性進(jìn)行自動判別；文獻(xiàn)[8]通過計算圖像灰度值的平均值和方差，以及方差曲線的變差方差根，構(gòu)建巖性的定量判別標(biāo)準(zhǔn)模板，進(jìn)行巖性定量識別；文獻(xiàn)[9]和文獻(xiàn)[10]從電成像測井圖像中提取定量表征礫石含量的結(jié)構(gòu)指示參數(shù)，結(jié)合常規(guī)測井資料建立交會圖，定量自動識別巖性。但這類基于圖像模板特征提取的識別方法，未考慮圖像的井眼覆蓋率，圖像中極板間的空白條帶未充填，使得樣本圖像的特征計算結(jié)果不準(zhǔn)確；并且這類方法依賴建立的圖像標(biāo)準(zhǔn)模板，標(biāo)準(zhǔn)模板的有限性導(dǎo)致識別各種多變的未知圖像時，出現(xiàn)準(zhǔn)確率低、泛化能力差等問題。

本文以圖像中礫巖顆粒提取為核心，首先采用MCA(Morphological Component Analysis)算法充填圖像中極板間的空白條帶，生成全井眼圖像,再進(jìn)行基于邊緣流的多尺度礫巖顆粒邊緣提取和基于曲線演化的礫巖顆粒區(qū)域分割，得到初步分割后的礫巖顆粒區(qū)域。但由于電成像測量技術(shù)、成像技術(shù)的限制和復(fù)雜井況，初步分割的顆粒區(qū)域與人工視覺判別結(jié)果有很大區(qū)別，因此在處理井段引入?yún)^(qū)塊地質(zhì)沉積狀況時，特別是在將相應(yīng)沉積條件下的礫巖顆粒磨圓度這一先驗知識作為約束條件時，自頂向下地驅(qū)動砂礫巖顆粒區(qū)域的合并和分裂。最后，構(gòu)建砂礫巖粒度譜和巖性精細(xì)剖面，得到砂礫巖巖性的精細(xì)評價成果。本文方法應(yīng)用于YJ油田深層砂礫巖難動用油藏的測井評價中，為優(yōu)選射孔井段和壓裂層段提供了準(zhǔn)確依據(jù)。

1 電成像圖像空白條帶充填

由于電成像測井在大多數(shù)井眼條件下都存在空白條帶，達(dá)不到全井眼覆蓋，為保證礫巖顆粒提取的準(zhǔn)確度，采用MCA算法[11]對空白條帶進(jìn)行充填。MCA將圖像分解為幾何結(jié)構(gòu)un和紋理ut兩種形態(tài)成分，并假設(shè)兩個形態(tài)成分在某個特定的基庫或過完備字典下是稀疏的，且基庫或過完備字典之間具有不相干性。MCA模型表示如下[12]：

(1)

圖1 MCA算法具體實現(xiàn)步驟

圖2 原始成像圖及處理結(jié)果示例

MCA算法具體實現(xiàn)結(jié)果如圖2(b)所示，可以看到圖2(a)的原始圖像經(jīng)算法充填后，較好地保持了礫巖顆粒邊緣的連續(xù)性。本文以每次迭代得到的充填圖像與原始圖像在掩膜模板上的平均像素灰度誤差作為迭代結(jié)束的條件。表1給出了圖2(a)原始圖像分別在迭代5，10，20，30，40，50和100次時的平均像素灰度誤差值。

由表1看到，算法迭代到20次時，誤差急劇下降，隨著迭代次數(shù)的增加，誤差變化不大，說明算法迭代20次以后，充填結(jié)果基本相同，實驗結(jié)果也證明了這一點(diǎn)。因此，通過適當(dāng)?shù)剡x擇誤差閾值(或兩次迭代的差異值)，可以自適應(yīng)地結(jié)束迭代，得到可以接受的充填結(jié)果，從而避免固定較大迭代次數(shù)的方法帶來的計算冗余。

表1 不同迭代次數(shù)下的平均像素灰度誤差值

2 基于邊緣流的多尺度礫巖顆粒邊緣提取

邊緣流算法[13]在初始尺度計算邊緣流能量，并將各個像素點(diǎn)的邊緣流能量投影到水平、垂直方向，進(jìn)行累加，得到初始邊緣流向量場；隨著尺度提升，當(dāng)初始向量場幅度小于設(shè)定幅度閾值時，當(dāng)前尺度的向量場等效于初始向量場；當(dāng)初始向量場與當(dāng)前向量場角度小于設(shè)定的角度閾值時，則將當(dāng)前向量場疊加上初始向量場。在多個尺度重復(fù)上述過程，直至到達(dá)最大設(shè)定尺度，得到最終的多尺度邊緣流向量場。最后，根據(jù)Possion方程，由最終的多尺度邊緣流向量場的散度計算得到圖像的邊緣流函數(shù)，結(jié)果如圖2(c)所示，邊緣流函數(shù)通過歸一化轉(zhuǎn)換為灰度圖像顯示，可以看到灰度值較小的礫巖顆粒邊緣，但有些邊緣并不連續(xù)。

3 基于曲線演化的礫巖顆粒區(qū)域分割

基于曲線演化的砂礫巖顆粒區(qū)域分割算法[14]，是根據(jù)上節(jié)計算得到的邊緣流函數(shù)，通過求解水平集方程，計算區(qū)域分割單像素曲線，其特點(diǎn)是可以將不閉合的斷裂邊緣進(jìn)行封閉。水平集方程如下所示[14]：

(2)

式中，φ為迭代計算曲線演化結(jié)果圖像；t為時間；S為邊緣能量場；▽為取散度；w為平滑權(quán)重系數(shù)；C為邊緣流函數(shù)；K為一個表示收縮或擴(kuò)張的常數(shù)。最后對曲線演化結(jié)果圖像分別以高斯函數(shù)水平、垂直方向的一階導(dǎo)數(shù)為卷積核，進(jìn)行卷積運(yùn)算，將幅度值大于設(shè)定閾值的像素點(diǎn)選作邊緣點(diǎn)，結(jié)果如圖2(d)所示，可以看到藍(lán)色的單像素寬度曲線封閉了各個礫巖顆粒區(qū)域。

4 區(qū)塊地質(zhì)沉積約束條件下的砂礫巖顆粒區(qū)域合并與分裂

從圖2(d)中可以看到，各個區(qū)域被閉合的藍(lán)色單像素曲線區(qū)分開來，但有的礫巖顆粒區(qū)域明顯與人眼視覺判別結(jié)果不同，出現(xiàn)明顯的陡峭凹陷。圖3為砂礫巖顆粒區(qū)域分裂示例。從圖3(a)的區(qū)域分割結(jié)果可以看出，同樣各個區(qū)域被閉合的藍(lán)色單像素曲線區(qū)分開來，但有的礫巖顆粒區(qū)域存在明顯的粘連情況，兩個鄰近的礫巖顆粒被劃分為一個區(qū)域。這些礫巖顆粒的形狀特征出現(xiàn)的概率在當(dāng)前地質(zhì)沉積條件下很低。

分析當(dāng)前區(qū)塊地質(zhì)沉積狀況，由此得到礫巖顆粒磨圓度的統(tǒng)計規(guī)律，作為約束條件知識驅(qū)動顆粒區(qū)域的合并與分裂過程。由顆粒磨圓度判定凹陷的顆粒分割區(qū)域是否進(jìn)行相鄰區(qū)域的合并，以及明顯粘連的顆粒分割區(qū)域是否進(jìn)行區(qū)域的分裂，由一個區(qū)域分裂為多個區(qū)域，并將每個分裂后的區(qū)域作為單獨(dú)的一個礫巖顆粒。具體方法是采用自適應(yīng)曲率計算技術(shù)，定位礫巖顆粒區(qū)域中邊界上的凹點(diǎn)、凹點(diǎn)對和鄰近凸點(diǎn)。由鄰近的非礫巖區(qū)域的相對位置，合并鄰近區(qū)域與該礫巖區(qū)域。同時，根據(jù)定位的凹點(diǎn)對之間的弦長和計算的對應(yīng)弧長，與弧弦比閾值比較，將該礫巖區(qū)域分裂為多個礫巖顆粒，這樣就解決了初步分割區(qū)域的不完整和粘連情況。

圖4為砂礫巖顆粒區(qū)域合并示例。圖4(a)是圖2(d)使用并查集提取的各個分割區(qū)域，得到各區(qū)域的包含像素點(diǎn)，并以不同顏色表示，根據(jù)礫巖灰度閾值，判定是否為礫巖顆粒。提取礫巖顆粒邊界像素，并以順時針排序。自適應(yīng)計算邊界曲率，定位凹點(diǎn)及其鄰近凸點(diǎn)，凹點(diǎn)和其鄰近凸點(diǎn)形成一黃色矩形區(qū)域，如圖4(b)所示。如果鄰近區(qū)域的質(zhì)心在該矩形區(qū)域內(nèi)，則合并鄰近區(qū)域；如圖4(c)中所示的左下方淺綠色表征的礫巖顆粒區(qū)域，經(jīng)過區(qū)域合并處理后，形成了完整的礫巖顆粒區(qū)域，符合當(dāng)前礫巖顆粒磨圓度的設(shè)定。圖3(b)為圖3(a)使用并查集提取的各個分割區(qū)域，計算礫巖區(qū)域各個凹點(diǎn)對之間的弦長以及對應(yīng)的弧長，并根據(jù)凹點(diǎn)的次數(shù)權(quán)重，使用權(quán)重大于1的凹點(diǎn)對對應(yīng)弦，切分該礫巖區(qū)域。圖3(b)中大面積的深藍(lán)色區(qū)域表明存在明顯的粘連情況；如圖3(c)中所示，該區(qū)域由紅色的凹點(diǎn)對之間相連的弦，分裂為兩個不同的區(qū)域(分別用深紅色與深藍(lán)色表征)，這樣，每個區(qū)域的形狀特征符合當(dāng)前顆粒磨圓度的知識約束。

圖4 砂礫巖顆粒區(qū)域合并示例

5 粒度譜構(gòu)建及精細(xì)巖性剖面計算

砂礫巖粒度譜的構(gòu)建和巖性精細(xì)剖面的計算采用滑動窗口方法，如圖5所示，統(tǒng)計某深度點(diǎn)對應(yīng)滑動窗口內(nèi)，礫巖顆粒區(qū)域面積的大小，并由此劃分區(qū)域面積的BIN值(面積大小的統(tǒng)計分組數(shù)目，每個BIN(Binary In Net)代表了各自的面積大小范圍)，統(tǒng)計屬于每個BIN值范圍內(nèi)的礫巖顆粒的數(shù)目，構(gòu)成粒度譜，作為解釋成果直觀地展現(xiàn)砂礫巖顆粒的大小分布。如圖6所示，精細(xì)巖性剖面計算首先根據(jù)礫巖顆粒大小閾值，將礫巖顆粒細(xì)分為細(xì)礫巖、中礫巖和粗礫巖，并根據(jù)滑動窗口內(nèi)各巖性顆粒面積除以滑動窗口面積，得到精細(xì)巖性剖面中各細(xì)分巖性的百分比，而砂巖巖性百分比由礫巖顆粒占比的殘余值計算，并形成曲線，由解釋軟件進(jìn)行充填繪制，形成直觀的砂礫巖精細(xì)巖性剖面。

圖5 砂礫巖粒度譜構(gòu)建原理

圖6 砂礫巖精細(xì)巖性剖面計算原理

6 應(yīng)用效果分析

砂礫巖巖性精細(xì)評價成果圖及與NRA巖性識別方法對比如圖7所示。

圖7 砂礫巖巖性精細(xì)評價成果圖及與NRA巖性識別方法對比

圖7(b)給出了Y920井3555.517～3557.808 m井段的精細(xì)評價綜合處理成果圖，可以看到構(gòu)建的粒度譜和精細(xì)巖性剖面，與圖7(a)所示的NRA(Neutron Resistivity Acoustic)巖性識別方法[15]成果圖進(jìn)行對比可以看到，在解釋井段，兩者的對應(yīng)性較好；解釋井段上部礫巖總組分偏少，對應(yīng)NRA方法的礫狀砂巖；下部礫巖總組分逐漸增加，對應(yīng)NRA方法的礫巖巖性。粒度譜直觀地展現(xiàn)了相應(yīng)井段礫巖顆粒的大小分布，巖性剖面相對于傳統(tǒng)的礫巖、礫狀砂巖和含礫砂巖的巖性劃分，更加精細(xì)地區(qū)分了礫巖顆粒的大小，而不是籠統(tǒng)地根據(jù)井段中所有礫巖部分所占面積劃分砂礫巖巖性，不考慮礫巖顆粒的獨(dú)立性和特殊性。

圖8為Y936井3676～3677.6 m井段處理成果圖，巖性剖面中以細(xì)礫巖和中礫巖為主，對比常規(guī)解釋結(jié)論(油層)，表明細(xì)礫巖、中礫巖的物性、含油性較好，符合YJ油田砂礫巖體區(qū)塊的儲層“四性”關(guān)系，據(jù)此優(yōu)選壓裂設(shè)計井段，壓裂成功實施后，獲得初產(chǎn)近20 t的產(chǎn)能，目前日穩(wěn)產(chǎn)油10 t。

圖8 Y936井3676～3677.6 m井段處理成果圖

7 結(jié)束語

本文的巖性自動識別方法通過圖像預(yù)處理以及在地質(zhì)沉積約束條件知識驅(qū)動下的精細(xì)圖像處理技術(shù)，有效地提取了圖像上的礫巖顆粒，并據(jù)此構(gòu)建粒度譜和精細(xì)的巖性剖面，為用戶提供了直觀、準(zhǔn)確的砂礫巖評價成果。

本文方法可在今后的老井復(fù)查和新井解釋中發(fā)揮重要作用，為開拓服務(wù)市場和油田的勘探開發(fā)提供重要的技術(shù)支撐。