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

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

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

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

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

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

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

(1)

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

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

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

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

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

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

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

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

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

(2)

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

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

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

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

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

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

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

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

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

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

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

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

圖7 砂礫巖巖性精細(xì)評(píng)價(jià)成果圖及與NRA巖性識(shí)別方法對(duì)比

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

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

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

7 結(jié)束語(yǔ)

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

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