基于數(shù)字巖心孔喉特征的等效孔隙縱橫比有效性研究

譚開俊， 趙建國， 滕團(tuán)余， 劉欣澤， 閆博鴻

1 中國石油勘探開發(fā)研究院西北分院， 蘭州 730020 2 中國石油大學(xué)(北京)油氣資源與探測(cè)國家重點(diǎn)實(shí)驗(yàn)室, 北京 102249 3 中國石油天然氣股份有限公司玉門油田環(huán)慶分公司， 甘肅酒泉 735019

0 引言

巖石微觀孔隙結(jié)構(gòu)形態(tài)對(duì)巖石物性具有重要的影響，以定量化儲(chǔ)層預(yù)測(cè)為目標(biāo)，油氣勘探中對(duì)巖石微觀孔隙結(jié)構(gòu)研究的要求日益提高.碳酸鹽巖儲(chǔ)層由于成巖后期強(qiáng)膠結(jié)、溶蝕、以及礦物充填等作用次生孔隙非常發(fā)育，這使其孔隙結(jié)構(gòu)極其復(fù)雜(Bathurst, 1972; Lucia, 1983)，常見的有鑄模孔、粒內(nèi)溶孔、粒間溶孔、晶間孔、裂縫等，Anselmetti和Eberli(1993, 1997)發(fā)現(xiàn)這些復(fù)雜多變的孔隙結(jié)構(gòu)對(duì)彈性性質(zhì)的控制作用甚至比孔隙度還大；Zhang等(2015)研究發(fā)現(xiàn)在中國四川碳酸鹽巖儲(chǔ)層中，在相同孔隙度下由于多種孔隙結(jié)構(gòu)的發(fā)育其縱波速度會(huì)出現(xiàn)約有1000 m·s-1的速度差異；Sayers(2008)也指出具有相同飽和度、孔隙度的數(shù)據(jù)點(diǎn)在解釋量版上的分布并不集中.這導(dǎo)致我們?cè)趦?chǔ)層預(yù)測(cè)時(shí)很容易將孔隙結(jié)構(gòu)對(duì)速度、彈性性質(zhì)造成的影響歸結(jié)在孔隙度、孔隙流體等其他因素上，使結(jié)果產(chǎn)生很大的偏差(李宏兵等，2019).另一方面，Anselmetti和Eberli(1999)、Eberli等(2003)還通過分析裂縫與孔洞的分布發(fā)現(xiàn)裂縫發(fā)育會(huì)顯著增強(qiáng)滲透率這一規(guī)律，并試圖用之來近似預(yù)測(cè)巖石滲透率.因此，孔隙結(jié)構(gòu)特征的研究對(duì)碳酸鹽巖等巖性的儲(chǔ)層預(yù)測(cè)、含油性評(píng)價(jià)具有重要意義.

為了盡可能定量化地預(yù)測(cè)與表征孔隙結(jié)構(gòu)，國內(nèi)外許多學(xué)者從多種角度進(jìn)行了大量的實(shí)驗(yàn)與理論研究.從巖石物理實(shí)驗(yàn)研究的角度出發(fā)，Anselmetti等(1998)通過圖像數(shù)字分析技術(shù)(Digital Image Analysis，DIA)提取了17塊白云巖鑄體薄片的孔隙尺寸分布、孔隙形狀(孔隙周長與孔隙面積的比值)等作為孔隙結(jié)構(gòu)表征參數(shù)，隨后又使用Wyllie時(shí)間平均方程計(jì)算速度與實(shí)測(cè)速度的差值為標(biāo)準(zhǔn)，定義了300多塊碳酸鹽巖巖心的孔隙結(jié)構(gòu)變化情況，最后將孔隙結(jié)構(gòu)用于滲透率的預(yù)測(cè)(Anselmetti and Eberli，1999；Eberli et al.，2003).Verwer等(2010)分析指出具有較高孔隙比表面的碳酸鹽巖樣品剪切模量會(huì)明顯降低.目前使用高分辨率的掃描電鏡還能夠更進(jìn)一步研究納米級(jí)的微小孔隙，三維全巖心CT成像技術(shù)可以更準(zhǔn)確地獲得巖石內(nèi)部骨架與孔隙的復(fù)雜幾何形狀(Coenen et al.，2004)，建立真實(shí)的巖石模型.此外隨著實(shí)驗(yàn)手段的進(jìn)步，Volokitin等(2001)將核磁共振技術(shù)應(yīng)用于巖心的孔隙結(jié)構(gòu)描述：首先通過離心過程獲取巖心處于不同飽和情況時(shí)的橫向弛豫時(shí)間T2譜，再由一定的標(biāo)定方法與理論轉(zhuǎn)換得到孔喉半徑分布頻率.還可以對(duì)巖心進(jìn)行恒速或高壓壓汞，通過進(jìn)汞與退汞壓力變化曲線與相應(yīng)理論依據(jù)能夠較準(zhǔn)確地推算孔隙、喉道、孔喉比大小及含量分布，實(shí)現(xiàn)孔隙結(jié)構(gòu)定量化研究(Wang et al.，2019；王永詩等，2021).

在巖石物理理論模型方面，眾多學(xué)者通常采用等效介質(zhì)理論將孔隙結(jié)構(gòu)與彈性性質(zhì)聯(lián)系起來.經(jīng)典的等效介質(zhì)理論主要有K-T模型(Kuster and Toks?z，1974)、Mori-Tanaka(MT)模型(Benveniste，1987)和微分等效介質(zhì)(Differential Equivalent Medium，DEM)模型(Berryman et al.，2002)等，這些模型均將孔隙理想化為縱橫比不同的橢球體，在此基礎(chǔ)上計(jì)算多孔巖石等效彈性模量.Xu和White(1995)考慮碎屑巖中的孔隙在砂巖組分與泥巖組分中具有不同的結(jié)構(gòu)(孔隙縱橫比)，以DEM、K-T模型為基礎(chǔ)發(fā)展了著名的Xu-White模型，用以估算碎屑巖儲(chǔ)層的橫波速度.Kumar和Han(2005)引入了大、中、小三種縱橫比的孔隙，給出了聯(lián)合DEM模型與Wyllie時(shí)間平均方程估算三種孔隙的縱橫比及體積分?jǐn)?shù)的實(shí)施流程.Xu和Payne(2009)應(yīng)用Kumar的工作思路并加入泥巖孔隙，對(duì)碳酸鹽巖測(cè)井曲線進(jìn)行了三種孔隙結(jié)構(gòu)及含量的預(yù)測(cè)，其對(duì)孔隙結(jié)構(gòu)的刻畫通過實(shí)測(cè)橫波速度與預(yù)測(cè)值的對(duì)比得到驗(yàn)證.Zhao等(2013)將Kumar的孔隙分類及計(jì)算方法應(yīng)用于地震數(shù)據(jù)，并通過電阻率成像測(cè)井對(duì)孔隙結(jié)構(gòu)進(jìn)行了驗(yàn)證.Pan等(2015)、趙建國等(2021a，b)基于數(shù)字巖心也對(duì)地震孔隙結(jié)構(gòu)預(yù)測(cè)方面進(jìn)行了研究.而李宏兵等(2013)認(rèn)為在研究孔隙結(jié)構(gòu)對(duì)彈性性質(zhì)的影響時(shí)，即使巖石中含有多種孔隙類型都可將其等效為具有單一縱橫比的孔隙來表示，并通過具有多重孔DEM模型與單一孔隙類型DEM模型之間的對(duì)比證明了這種等效的合理性.除了上述將巖石中孔隙等效為一種或有限的幾種結(jié)構(gòu)外，David和Zimmerman(2012)考慮了巖石中含多種孔隙的情況，并通過變壓力超聲干燥測(cè)量速度數(shù)據(jù)與DEM理論反演得到了具有不同縱橫比孔隙的分布，結(jié)果通過飽和測(cè)量數(shù)據(jù)進(jìn)行驗(yàn)證.國內(nèi)歐陽芳等(2021)也通過超聲實(shí)驗(yàn)與多種等效介質(zhì)理論對(duì)此方面進(jìn)行了分析.然而這種方法很難應(yīng)用于測(cè)井、地震等大尺度的情況.

總體上目前對(duì)孔隙結(jié)構(gòu)的描述可以分成兩大類，一種是從實(shí)驗(yàn)角度獲取，包括薄片、CT掃描數(shù)字巖心、掃描電鏡等直接成像手段，以及通過核磁共振、壓汞等物理實(shí)驗(yàn)方法得到孔隙、喉道半徑等微觀幾何形狀參數(shù)；另外一種是將預(yù)先獲取的彈性模量聯(lián)合相應(yīng)的等效介質(zhì)理論，通過反演的思路計(jì)算等效孔隙縱橫比或類似參數(shù)，相比其他孔隙結(jié)構(gòu)參數(shù)而言基于等效介質(zhì)理論孔隙縱橫比更容易獲得，并且能夠應(yīng)用于測(cè)井、地震等大尺度數(shù)據(jù).然而，等效孔隙縱橫比畢竟是從彈性性質(zhì)出發(fā)通過理論間接計(jì)算得到的，能否確切地表征巖石真實(shí)的復(fù)雜孔隙結(jié)構(gòu)特征，在于其是否能夠與巖石中直接觀測(cè)到的孔隙幾何形狀參數(shù)之間建立起定量的關(guān)系(Baechle et al.，2008).雖然Anselmetti等(1998)、Baechle等(2008)使用DIA技術(shù)能夠提取鑄體薄片的孔隙結(jié)構(gòu)參數(shù)，可以與DEM模型中的孔隙縱橫比建立一定的聯(lián)系，但是鑄體薄片畢竟屬于二維圖像，與三維巖心實(shí)際的孔隙空間相比有很大的差距，且不同位置的薄片結(jié)構(gòu)變化或許很大.此外，使用薄片的孔隙結(jié)構(gòu)參數(shù)提取缺乏明確的孔喉分割理論支撐，可能不夠精準(zhǔn).因而如果能夠定量化地建立孔隙縱橫比與實(shí)際巖石三維圖像的孔隙特征參數(shù)之間的聯(lián)系，那么在利用等效介質(zhì)理論表征孔隙結(jié)構(gòu)時(shí)，其合理性與準(zhǔn)確性就具有了強(qiáng)有力的支撐，CT掃描技術(shù)無疑是最合適的手段.

本文以CT掃描數(shù)字巖心為基礎(chǔ)，從模擬的彈性性質(zhì)和圖像提取的孔隙網(wǎng)絡(luò)幾何特征兩方面研究了孔隙結(jié)構(gòu)的定量化表征方法，建立了等效孔隙縱橫比與孔喉特征參數(shù)間的定量關(guān)系：我們首先分別對(duì)基于數(shù)字巖心的彈性模擬技術(shù)與孔隙網(wǎng)絡(luò)提取方法進(jìn)行介紹；然后，對(duì)砂巖與碳酸鹽巖兩類不同巖性巖心進(jìn)行CT掃描，并對(duì)掃描圖像進(jìn)行圖像處理以保證能精確重構(gòu)三維數(shù)字巖心模型；接下來將數(shù)字巖心分割為子塊，應(yīng)用上述兩種方法獲得每個(gè)子塊的等效彈性模量與多種孔喉幾何特征參數(shù)；最后，構(gòu)建等效介質(zhì)理論與模擬彈性模量之間的損失函數(shù)方程，通過優(yōu)化算法來計(jì)算等效孔隙縱橫比，再分別與各子塊的不同孔喉幾何特征參數(shù)進(jìn)行了交會(huì)分析， 最終發(fā)現(xiàn)平均喉道長度和等效孔隙縱橫比有較為明顯的線性擬合關(guān)系，表明它們對(duì)孔隙結(jié)構(gòu)的表征具有良好的一致性.

1 基于數(shù)字巖心的彈性模擬與孔隙網(wǎng)絡(luò)提取

對(duì)于物體彈性性質(zhì)的數(shù)值模擬，我們使用的是Garboczi和Day(1995)提出的有限元模擬方法，它使用了靜態(tài)線彈性方程的變分形式.在數(shù)值計(jì)算時(shí)首先給定一初始位移場(chǎng)，利用共軛梯度算法對(duì)線彈性方程離散后的有限元控制方程進(jìn)行迭代求解，計(jì)算三維數(shù)字巖心數(shù)據(jù)體內(nèi)任意一體素的位移量，由位移得到應(yīng)力應(yīng)變值并平均得到宏觀量，再使用胡克定律即可計(jì)算巖石等效彈性模量.這種方法以彈性模量會(huì)影響物體中的應(yīng)力應(yīng)變分布為前提，根據(jù)最小位能原理，物體在受到應(yīng)力產(chǎn)生應(yīng)變后，會(huì)向著應(yīng)變位能減小的趨勢(shì)變化，最終在應(yīng)變位能最小時(shí)達(dá)到平衡狀態(tài)，利用此時(shí)物體的應(yīng)力應(yīng)變，計(jì)算得到的是物體處于小應(yīng)變狀態(tài)下的等效彈性模量，因此是一種靜態(tài)的模擬.

從圖形幾何學(xué)的角度，在利用數(shù)字巖心圖像直接確定孔隙結(jié)構(gòu)特征方面，我們使用了孔隙網(wǎng)絡(luò)法.這主要是因?yàn)樵诶硐肭闆r下，我們期望巖心內(nèi)部的孔隙能夠通過聯(lián)通、分割算法(夏良正和李久賢，1999)將每部分連在一起的孔隙空間單獨(dú)標(biāo)定出來，這樣就可以很容易確定孔徑、球度、表比面積等幾何形狀特征參數(shù).然而大多數(shù)聯(lián)通在一起的孔隙空間都是極度不規(guī)則的，很難定量評(píng)價(jià)形狀特征，或者有些巖心中的孔隙空間是整體聯(lián)通的，無法對(duì)其進(jìn)行分割與評(píng)價(jià)，因此需要采用特定的方法將較為復(fù)雜的孔隙空間轉(zhuǎn)化成一種形狀比較規(guī)則的空間來等價(jià)表示.受到滲流性質(zhì)模擬領(lǐng)域的啟發(fā)，在研究中我們采用最大球方法提取了三維數(shù)字巖心的孔隙網(wǎng)絡(luò).這種方法首先將不規(guī)則的孔隙空間轉(zhuǎn)化成了由大量內(nèi)切球體充填的空間，對(duì)所有球體進(jìn)行一定的劃分后將孔隙空間中較寬的區(qū)域定義為孔隙，孔隙之間的窄區(qū)域通過圓柱型的喉道相連接，兩者組合成為孔隙網(wǎng)絡(luò).這種網(wǎng)絡(luò)能夠近似等效地表示巖心的孔隙空間，因此其幾何形狀特征參數(shù)就能夠表示孔隙結(jié)構(gòu)特征(Dong and Blunt, 2009).

(1)

(2)

其中，ra、rb表示孔隙對(duì)應(yīng)的最大球a、b的半徑，rt表示喉道最大球的半徑，k表示孔喉長度分割系數(shù)，喉道長度lt定義為喉道總長度lab減去兩個(gè)孔隙長度la、lb后的值.以上就能夠得到孔隙網(wǎng)絡(luò)所有的幾何特征參數(shù). 從式(1)、(2)中可以發(fā)現(xiàn)，這種孔喉劃分方法是在孔、喉對(duì)應(yīng)的最大球半徑與分割系數(shù)的約束下分別定義了孔隙與喉道的范圍，在喉道長度范圍內(nèi)的所有最大球都?xì)w屬于喉道，而在孔隙長度范圍內(nèi)的所有最大球則歸屬于孔隙.這種方法雖然能夠簡便的獲得喉道與孔隙的長度特征，但是其中分割系數(shù)是一個(gè)人為控制的比例系數(shù)，Dong和Blunt (2009)在計(jì)算時(shí)給定為0.6，在本文中為了分析該參數(shù)對(duì)孔喉分割的影響，我們分別計(jì)算了不同系數(shù)值時(shí)的喉道長度參數(shù)并進(jìn)行對(duì)比.

圖1 孔隙-喉道劃分示意圖Fig.1 Pore-throat segmentation

為了更好地理解孔隙空間轉(zhuǎn)換為孔隙網(wǎng)絡(luò)后孔、喉位置對(duì)應(yīng)的孔隙空間區(qū)域，為孔隙結(jié)構(gòu)的定量分析提供基礎(chǔ)，我們構(gòu)建了由球體規(guī)則排列而成的巖石模型并提取孔隙網(wǎng)絡(luò)，見圖2.孔隙網(wǎng)絡(luò)中8個(gè)較寬的區(qū)域被定義為孔隙(紅色球體)，孔隙間的連接區(qū)域?yàn)楹淼?藍(lán)色圓柱體)，這有效地實(shí)現(xiàn)了不同孔型孔隙空間的定量判別.為進(jìn)一步分析孔隙網(wǎng)絡(luò)模型是否能夠?qū)r石宏觀整體的孔隙結(jié)構(gòu)有比較定量化的表征，我們嘗試提取了有明顯孔隙結(jié)構(gòu)特征巖心的孔隙網(wǎng)絡(luò).圖3a中的數(shù)字巖心樣品為典型的裂縫性孔隙結(jié)構(gòu)，孔隙空間成扁平狀延展，圖3c中為孔洞型的數(shù)字巖心樣品，其孔隙空間多為連續(xù)的類球狀，延展小數(shù)量多，分布較為均勻；圖3b與3d為對(duì)兩者提取的孔隙網(wǎng)絡(luò)模型，可以發(fā)現(xiàn)裂縫型樣品對(duì)應(yīng)的孔隙網(wǎng)絡(luò)模型中，藍(lán)色圓柱狀喉道的長度大多相比孔洞型樣品的喉道要明顯更長一些，裂縫幾乎都由喉道組成.而孔洞型樣品的喉道較為短小，基本由孔隙與喉道連接而組成.統(tǒng)計(jì)喉道長度、半徑頻率分布譜，發(fā)現(xiàn)兩種喉道相關(guān)參數(shù)在裂縫型樣品中都要大于孔洞型樣品，分布范圍也更寬，但是其中喉道長度差別更明顯，見圖4、5，因此在下文中我們重點(diǎn)應(yīng)用喉道長度對(duì)孔隙結(jié)構(gòu)進(jìn)行定量化分析.

2 數(shù)字巖心的圖像預(yù)處理

為了分析與驗(yàn)證孔隙網(wǎng)絡(luò)模型中的孔喉參數(shù)能否定量化描述巖石的孔隙結(jié)構(gòu)，我們分別選擇了孔隙結(jié)構(gòu)較為復(fù)雜的碳酸鹽巖樣品和孔隙結(jié)構(gòu)相對(duì)均一的砂巖樣品，共同來進(jìn)行孔隙結(jié)構(gòu)定量化分析與對(duì)比.所有樣品都為切割好的直徑38 mm的圓柱體，如圖6、7所示，使用本實(shí)驗(yàn)室的CT設(shè)備掃描成像后得到的圖像分辨率可達(dá)每體素20.7678 μm.為了提取后續(xù)運(yùn)算所需的孔隙結(jié)構(gòu)信息，CT掃描原始圖像需要經(jīng)過對(duì)比度增強(qiáng)、濾波去噪、相邊緣增強(qiáng)以及圖像二值化分割幾個(gè)主要的圖像處理步驟，本文以圖6中1號(hào)白云巖巖心為例來展示圖像處理的過程及效果，見圖8，處理僅對(duì)從原始圖像中心剪裁的邊長為1000體素的正方體圖像進(jìn)行.

圖2 最大球法提取簡單球體堆積模型的孔隙網(wǎng)絡(luò)(a) 球體堆積模型骨架空間； (b) 球體堆積模型孔隙空間； (c) 球體堆積模型孔隙網(wǎng)絡(luò).Fig.2 Process of pore-network extracted by maximal ball method for a simple sphere packing model(a) Skeleton of sphere packing model； (b) Pore space of sphere packing model； (c) Pore-network of sphere packing model.

圖3 不同孔隙類型數(shù)字巖心樣品孔隙網(wǎng)絡(luò)提取(a) 典型裂縫型樣品數(shù)字巖心模型； (b) 裂縫型巖心孔隙網(wǎng)絡(luò)模型； (c) 典型孔洞型樣品數(shù)字巖心模型； (d) 孔洞型巖心孔隙網(wǎng)絡(luò)模型.Fig.3 Extraction of pore-network for digital core of different pore structure(a) Digital core model of typical fracture core； (b) Pore-network model of fracture core； (c) Digital core model of typical porous core； (d) Pore-network model of porous core.

圖4 裂縫、孔洞型樣品喉道長度分布譜Fig.4 Histogram of throat length in fracture and porous cores

圖5 裂縫、孔洞型樣品喉道半徑分布譜Fig.5 Histogram of throat radius in fracture and porous cores

對(duì)數(shù)字巖心的圖像處理是進(jìn)行彈性性質(zhì)模擬及孔隙網(wǎng)絡(luò)模型提取的關(guān)鍵步驟，直接決定著孔隙結(jié)構(gòu)定量化表征的準(zhǔn)確性.圖6中第二排展示了7塊碳酸鹽巖樣品CT掃描成像處理后的三維孔隙空間分布情況，可以發(fā)現(xiàn)該類樣品具有比較復(fù)雜的孔隙結(jié)構(gòu)類型、非均質(zhì)性較強(qiáng)、孔隙度都很小，圖7中第二排展示了砂巖巖心的孔隙三維分布，直觀看來砂巖樣品的孔隙結(jié)構(gòu)明顯比較單一且分布均勻，和碳酸鹽巖樣品的區(qū)別較大，孔隙度也遠(yuǎn)大于碳酸鹽巖.

圖6 碳酸鹽巖樣品Fig.6 Carbonate sample

圖7 砂巖樣品Fig.7 Sandstone sample

圖8 CT掃描圖像處理中間過程及結(jié)果(a) 原始掃描圖像切片； (b) 對(duì)比度增強(qiáng)結(jié)果； (c) 濾波； (d) 二值化分割結(jié)果.Fig.8 Result and process of CT image processing(a) Original scanned image slice； (b) Contrast enhancement result； (c) Filtering for image； (d) Binarized image segmentation result.

在后續(xù)的模擬運(yùn)算前，我們還需對(duì)整個(gè)數(shù)字巖心數(shù)據(jù)進(jìn)行子塊分割，這主要有兩方面的原因，首先是受限于計(jì)算機(jī)的處理能力、內(nèi)存和運(yùn)算所需要的時(shí)間，一般較難處理10003個(gè)體素的數(shù)據(jù)量；另外單個(gè)或少量巖心計(jì)算結(jié)果較為單一，無法分析巖心性質(zhì)的統(tǒng)計(jì)性規(guī)律.對(duì)此Dvorkin等(2011)認(rèn)為如果將數(shù)字巖心內(nèi)部子塊數(shù)據(jù)集的模擬數(shù)據(jù)點(diǎn)組成趨勢(shì)線，那么這將在多個(gè)尺度上都比較穩(wěn)定，并且這樣的趨勢(shì)線不會(huì)隨巖心取樣空間位置和尺度的變化而產(chǎn)生劇烈波動(dòng)， 因此我們需要將大正方體的巖心分割成大小合適的子塊數(shù)據(jù).但是分割后的巖心子塊也應(yīng)足夠大，這樣才能包含足夠充分的孔隙信息，與真實(shí)巖心的孔隙結(jié)構(gòu)特征也更相近.姜黎明等(2012)通過代表體積元(REV)方法測(cè)試發(fā)現(xiàn)，巖心塊尺寸選在2003時(shí)通常能夠包含較豐富的孔隙結(jié)構(gòu)且計(jì)算較快，因此本文使用了該尺寸進(jìn)行子塊分割，之后就可以對(duì)子塊巖心進(jìn)行彈性性質(zhì)模擬與孔隙網(wǎng)絡(luò)提取.

在彈性性質(zhì)模擬時(shí)還需給定巖心中骨架與孔隙部分各自的彈性模量， 那么根據(jù)等效介質(zhì)理論的思想，巖石的等效彈性性質(zhì)除孔隙結(jié)構(gòu)、孔隙度外，還會(huì)受到骨架的礦物組成與巖石飽和流體情況兩個(gè)因素的影響.其中碳酸鹽巖的骨架礦物組成較為純凈，常為方解石或白云石.本文樣品經(jīng)XRD(X-ray diffraction，X射線衍射)礦物分析發(fā)現(xiàn)白云石占90%以上，礦物組成單一，因此在比較所有碳酸鹽巖彈性模量時(shí)骨架礦物組成對(duì)各樣品實(shí)際彈性性質(zhì)的影響應(yīng)很小，在此我們給定骨架部分為白云巖的彈性參數(shù)(體積模量為80 GPa，剪切模量為58 GPa，該模量通過對(duì)同井位中、孔隙度0.1%左右的一非常致密巖心測(cè)量得到)；與之相反砂巖骨架的礦物組成通常比較復(fù)雜，因此在對(duì)比砂巖彈性性質(zhì)時(shí)礦物成分是重要的影響因素，但是考慮到我們的研究核心在于孔隙結(jié)構(gòu)的定量化表征，孔隙結(jié)構(gòu)對(duì)彈性性質(zhì)的影響是我們主要考慮的因素，所以仍舊需要假定巖石骨架礦物成分單一，這樣所有樣品的模擬彈性模量才僅會(huì)受到孔隙結(jié)構(gòu)的影響，而不會(huì)受不同樣品礦物組成不同所影響.考慮到較純凈的砂巖(如楓丹白露砂巖)基本為石英，本文砂巖經(jīng)XRD分析石英含量也基本占主導(dǎo)，因此給定砂巖的骨架部分為石英的彈性參數(shù)(體積模量為37 GPa，剪切模量為44 GPa；Mavko et al.，1998).同樣為了排除孔隙流體對(duì)樣品模擬彈性模量的影響，我們假設(shè)兩種巖性樣品的孔隙中都完全飽和空氣.

3 等效孔隙結(jié)構(gòu)定量化表征對(duì)比

根據(jù)研究思路，我們首先聯(lián)合模擬彈性模量與等效介質(zhì)理論估算等效孔隙縱橫比，以該參數(shù)表征數(shù)字巖心整體等效的孔隙結(jié)構(gòu)特征.目前已經(jīng)發(fā)展了多種常用的以孔隙縱橫比來表征孔隙形狀的等效介質(zhì)模型，它們的假設(shè)條件、適用情況以及加入孔隙的方式都有所不同，考慮到計(jì)算效率本文使用了較為簡單的顯式MT模型，但是為了對(duì)比不同理論計(jì)算的等效孔隙縱橫比是否會(huì)對(duì)孔隙結(jié)構(gòu)表征產(chǎn)生顯著的影響， 我們也給出了DEM模型的計(jì)算結(jié)果，見附錄B.通過等效孔隙縱橫比分布的對(duì)比表明，等效介質(zhì)理論不同，對(duì)最后孔隙結(jié)構(gòu)刻畫造成的影響很小.便于計(jì)算的包含N種包裹物的MT模型表示為式(3)(Benveniste，1987)：

(3)

圖9a展示了所有碳酸鹽巖數(shù)字巖心模擬縱波速度與孔隙度交會(huì)特征，可以發(fā)現(xiàn)交會(huì)圖中數(shù)據(jù)點(diǎn)分布總體比較分散，在孔隙度為10%左右時(shí)縱波速度的差異高達(dá)2500 m·s-1左右，這只能是不同的孔隙結(jié)構(gòu)引起的；而在砂巖數(shù)字巖心模擬縱波速度與孔隙度交會(huì)圖中數(shù)據(jù)點(diǎn)非常集中，呈明顯的近線性分布，同一孔隙度下的縱波速度差距很小，見圖9b.在求取等效孔隙縱橫比時(shí)將巖心孔隙度代入式(3)中，計(jì)算的等效模量作為輸出，以模擬彈性模量Kmodeling、μmodeling作為約束值作差，構(gòu)建包含體積模量與剪切模量兩個(gè)約束條件下的方程組(4)，通過優(yōu)化算法求解就可以得到當(dāng)前彈性模量與孔隙度條件下對(duì)應(yīng)的等效孔隙縱橫比.式(4)中OF表示目標(biāo)函數(shù).

(4)

圖9中每個(gè)數(shù)據(jù)點(diǎn)的顏色代表該點(diǎn)等效孔隙縱橫比，顏色越偏暖色值表示縱橫比值越大，反之縱橫比值越小.圖中藍(lán)色虛線表示理論計(jì)算的不同等效孔隙縱橫比時(shí)縱波速度隨孔隙度變化曲線.在碳酸鹽巖樣品中MT模型計(jì)算的等效縱橫比分布在0.03～0.4范圍內(nèi)，分布區(qū)間比較寬(見頻率直方圖10中紅色所示)，同一孔隙度下都表現(xiàn)為速度越快等效孔隙縱橫比越大，總體上所有數(shù)據(jù)點(diǎn)可以根據(jù)顏色分為下部的冷色值與上部的暖色值兩類；砂巖樣品的等效縱橫比分布在0.2～0.5之間，分布區(qū)間都要小于碳酸鹽巖(見頻率直方圖10中藍(lán)色所示)，更加集中.

圖9 數(shù)字巖心縱波速度-孔隙度交會(huì)中等效孔隙縱橫比分布對(duì)比(色標(biāo)表示孔隙縱橫比)(a) 碳酸鹽巖； (b) 砂巖.Fig.9 Distribution of effective pore aspect ratio of digital core in cross-plot between P-velocity and porosity (color represent pore aspect ratio)(a) Carbonate； (b) Sandstone.

圖10 砂巖與碳酸鹽巖等效孔隙縱橫比頻率分布直方圖Fig.10 Histogram of effective pore aspect ratio in sandstone and carbonate

在使用巖心孔隙網(wǎng)絡(luò)模型提取的孔喉參數(shù)對(duì)孔隙結(jié)構(gòu)進(jìn)行表征時(shí)，考慮到圖4、5中的分析，并且喉道此時(shí)更能表現(xiàn)孔隙空間中窄區(qū)域的特征，從而影響孔隙形狀，我們提取了每塊數(shù)字巖心所有喉道的平均長度、長度排序后的中間長度、經(jīng)過喉道體積加權(quán)平均的喉道長度和長度分布頻率直方圖峰值4個(gè)喉道長度相關(guān)的統(tǒng)計(jì)性參數(shù)作為代表每塊巖心的等效喉道長度，并與等效孔隙縱橫比(基于MT模型計(jì)算)進(jìn)行了交會(huì)分析.在圖11中孔喉長度分割系數(shù)為0.6時(shí)可以明顯發(fā)現(xiàn)碳酸鹽巖的不同喉道長度參數(shù)相比砂巖都具有更寬的分布范圍，這種特征和等效孔隙縱橫比的分布特征是一致的；并且碳酸鹽巖各喉道長度參數(shù)與等效孔隙縱橫比交會(huì)都出現(xiàn)了近似線性趨勢(shì)，其中喉道平均長度的相關(guān)系數(shù)明顯更大、相關(guān)性更強(qiáng)，數(shù)據(jù)分布更集中，表現(xiàn)為縱橫比增大對(duì)應(yīng)喉道長度減小的負(fù)相關(guān)性.這可以解釋為縱橫比越小的孔隙越呈扁平狀分布，延展更長，因此在提取的孔隙網(wǎng)絡(luò)中會(huì)有更長的喉道長度，而縱橫比越大的孔隙越接近球體，對(duì)應(yīng)的喉道會(huì)縮短.圖12將兩類巖心單獨(dú)分析發(fā)現(xiàn)砂巖平均喉道長度與等效孔隙縱橫比交會(huì)無明顯相關(guān)性，而碳酸鹽巖的交會(huì)關(guān)系相比圖11c中兩類巖心時(shí)基本不變，因此巖性并不會(huì)影響兩參數(shù)的相關(guān)性，圖11c中的交會(huì)關(guān)系對(duì)兩種巖性都是適用的.隨后為了分析孔喉分割系數(shù)的影響，我們使用0.4與0.8兩個(gè)不同系數(shù)重新計(jì)算了數(shù)字巖心喉道平均長度，并與等效孔隙縱橫比進(jìn)行了交會(huì)分析，見圖13.其中雖然長度值的范圍發(fā)生了變化，但是所有樣品數(shù)據(jù)點(diǎn)的相對(duì)分布規(guī)律并沒有明顯的改變，這表明分割系數(shù)對(duì)我們的孔隙結(jié)構(gòu)表征研究敏感度較小.根據(jù)以上交會(huì)分析我們選擇了平均喉道長度作為孔隙結(jié)構(gòu)的表征參數(shù)，同樣使用縱波速度與孔隙度做交會(huì)，每個(gè)數(shù)據(jù)點(diǎn)顏色表示平均喉道長度值(由于喉道長度非常小，為了便于比較，我們將所有樣品的平均喉道長度以其中的最大值為基準(zhǔn)進(jìn)行了歸一化)，顏色越偏冷色值代表喉道越短，反之喉道越長.孔隙結(jié)構(gòu)表征結(jié)果見圖14所示，藍(lán)色虛線為MT模型計(jì)算的不同等效孔隙縱橫比時(shí)縱波速度隨孔隙度變化曲線，作為以喉道長度來定量表征孔隙結(jié)構(gòu)的參考.我們能夠明顯看出圖14中數(shù)據(jù)點(diǎn)顏色分布與圖9中相似，碳酸鹽巖能夠大致分為喉道較短的冷色點(diǎn)與較長的暖色點(diǎn)兩類，砂巖則都表現(xiàn)為比較單一的喉道長度.

圖11 孔喉長度分割系數(shù)0.6時(shí)各喉道長度參數(shù)-等效孔隙縱橫比交會(huì)(a) 加權(quán)喉道長度-等效孔隙縱橫比交會(huì)； (b) 喉道長度頻率直方圖峰值-等效孔隙縱橫比交會(huì)； (c) 喉道平均長度-等效孔隙縱橫比交會(huì)； (d) 喉道長度中值-等效孔隙縱橫比交會(huì).Fig.11 Cross-plot between effective pore aspect ratio and parameters which described throat length when pore-throat segmentation coefficient is 0.6(a) Effective pore aspect ratio and weighted throat length； (b) Effective pore aspect ratio and peak value of throat length histogram； (c) Effective pore aspect ratio and average value of throat length； (d) Effective pore aspect ratio and mid-value of throat length.

圖12 不同巖心中孔喉長度分割系數(shù)0.6時(shí)喉道平均長度-等效孔隙縱橫比交會(huì)(a) 碳酸鹽巖； (b) 砂巖.Fig.12 Cross-plot between effective pore aspect ratio and average value of throat length when pore-throat segmentation coefficient is 0.6 in different core(a) Carbonate； (b) Sandstone.

圖13 平均喉道長度-等效孔隙縱橫比交會(huì)(a) 孔喉長度分割系數(shù)0.4； (b) 孔喉長度分割系數(shù)0.8.Fig.13 Cross-plot between effective pore aspect ratio and average value of throat length(a) Pore-throat segmentation coefficient is 0.4；(b) Pore-throat segmentation coefficient is 0.8.

圖14 縱波速度-孔隙度交會(huì)圖中數(shù)字巖心喉道平均長度分布(色標(biāo)表示歸一化喉道平均長度)(a) 碳酸鹽巖； (b) 砂巖.Fig.14 Distribution of average length of throat from digital core in cross-plot between P-velocity and porosity (color represent normalized average throat length)(a) Carbonate； (b) Sandstone.

4 討論與誤差分析

經(jīng)過對(duì)以上兩種孔隙結(jié)構(gòu)表征思路結(jié)果的對(duì)比與分析，證明雖然等效孔隙縱橫比是通過彈性性質(zhì)這一間接屬性來定量化描述孔隙結(jié)構(gòu)的，但是其對(duì)孔隙結(jié)構(gòu)的表征效果與直接從巖心圖像中提取的喉道長度的表征效果是相近的.然而目前我們的巖心樣品中等效孔隙縱橫比與平均喉道長度的交會(huì)關(guān)系并不是非常理想，因此需要對(duì)兩參數(shù)可能產(chǎn)生的誤差進(jìn)行分析.這之中等效孔隙縱橫比由等效介質(zhì)理論與彈性模量計(jì)算得到，其誤差應(yīng)來自于兩個(gè)方面：一個(gè)是彈性模擬的結(jié)果是否相對(duì)準(zhǔn)確，另一個(gè)則是等效介質(zhì)理論(MT模型)能否準(zhǔn)確描述彈性模量和等效孔隙縱橫比之間的關(guān)系.對(duì)此我們建立已知孔隙縱橫比的模型介質(zhì)，通過對(duì)比其模擬模量與等效介質(zhì)理論計(jì)算模量間的差距來評(píng)估誤差.首先使用位置隨機(jī)分布、半徑固定的球體作為孔隙，建立理論上孔隙縱橫比為1的多孔介質(zhì)；再使用4參數(shù)隨機(jī)生長法(Quartet Structure Generation Set，QSGS；Wang et al.，2007)，通過控制生長參數(shù)生成長軸與短軸比不同的橢球體作為孔隙(同一個(gè)模型介質(zhì)中每個(gè)橢球體長短軸都相同)，建立不同孔隙縱橫比的多孔介質(zhì)，具體的長、短軸長度可以通過對(duì)圖像進(jìn)行測(cè)量得到.以上多孔介質(zhì)不考慮孔隙間的連通性，各個(gè)孔隙互不相交，見圖15.

模擬使用石英作為骨架，孔隙中充填空氣.圖16中不同標(biāo)記點(diǎn)表示模型介質(zhì)的模擬模量，粉色虛線表示MT模型計(jì)算的彈性模量隨孔隙度變化曲線，兩者所用孔隙縱橫比互相對(duì)應(yīng).對(duì)比發(fā)現(xiàn)兩模量相差較小，這表明對(duì)于簡單的、孔隙結(jié)構(gòu)單一且分布均勻的多孔介質(zhì)，彈性模擬結(jié)果是可靠的，等效介質(zhì)理論也基本能準(zhǔn)確計(jì)算彈性模量.但是在孔隙結(jié)構(gòu)復(fù)雜時(shí)等效介質(zhì)理論的誤差很難評(píng)價(jià)，這是因?yàn)楹茈y像圖15中的簡單模型一樣通過測(cè)量孔隙長短軸定量化提取復(fù)雜孔隙的等效縱橫比，而只能由彈性模量反演計(jì)算一個(gè)理論值.因此裂縫型孔隙定向排列引起的各向異性、孔隙位置與區(qū)域的非均勻分布都可能使巖石不滿足等效介質(zhì)理論的假設(shè)，從而使等效孔隙縱橫比無法反映真實(shí)的孔隙結(jié)構(gòu)特征，導(dǎo)致圖11c中的擬合誤差.另外等效孔隙縱橫比也無法反映孔隙之間的聯(lián)通特性，忽略了孔隙結(jié)構(gòu)對(duì)滲流性質(zhì)的影響.

喉道長度參數(shù)的誤差主要來源于從孔隙空間提取孔隙網(wǎng)絡(luò)的過程.通過圖2對(duì)一球體顆粒堆積形成的簡單孔隙空間提取孔隙網(wǎng)絡(luò)后發(fā)現(xiàn)，雖然最大球法對(duì)孔隙空間中局部寬區(qū)域的孔隙與局部窄區(qū)域的喉道有明確的定義方法，但是對(duì)兩者中間剩余空間的歸屬還是比較模糊的，這也就導(dǎo)致了確定的喉道長度與孔隙長度可能并不準(zhǔn)確.在上文中我們分析了最大球法分割孔喉時(shí)分割系數(shù)對(duì)喉道長度的影響，發(fā)現(xiàn)該參數(shù)基本只會(huì)影響喉道長度值的整體范圍，對(duì)所用樣品喉道長度分布的相對(duì)趨勢(shì)與規(guī)律性沒有影響.但是孔隙網(wǎng)絡(luò)最開始是以研究滲流特性為目的而展開的，為了簡便而將孔隙空間轉(zhuǎn)化為孔隙網(wǎng)絡(luò)，并認(rèn)為孔隙網(wǎng)絡(luò)能夠近似代替實(shí)際的孔隙空間，這種近似可能不會(huì)對(duì)滲流特性產(chǎn)生明顯影響，但是對(duì)原本孔隙結(jié)構(gòu)的表征肯定存在誤差；此外網(wǎng)絡(luò)模型也相當(dāng)于將原本聯(lián)通在一起的孔隙空間分割開來，即使分割過程有最大球法作為理論依據(jù)，分割的位置也可能會(huì)對(duì)實(shí)際孔隙結(jié)構(gòu)的刻畫產(chǎn)生影響，這些因素都可能導(dǎo)致同一喉道長度下出現(xiàn)多個(gè)等效孔隙縱橫比，使兩者間的關(guān)系出現(xiàn)誤差.因此在分析等效孔隙縱橫比定量化表征復(fù)雜孔隙結(jié)構(gòu)的誤差方面，我們對(duì)彈性模擬、圖像形狀的定量化表征、孔隙網(wǎng)絡(luò)提取及孔喉分割還需要更多的研究和嘗試，對(duì)于孔隙結(jié)構(gòu)復(fù)雜、非均質(zhì)性強(qiáng)的碳酸鹽巖也需要更多的研究樣品，才能獲得更加可靠的等效孔隙縱橫比與平均喉道長度間的經(jīng)驗(yàn)關(guān)系.

圖15 不同單一孔隙縱橫比模型(a) 球狀孔隙模型，孔隙縱橫比為1； (b) 橢球狀孔隙模型， 控制孔隙縱橫比定量變化.Fig.15 Model in different single pore aspect ratio(a) Sphere pore model, pore aspect ratio is 1； (b) Ellipsoid pore model, pore aspect ratio is designed quantitatively.

圖16 不同單一孔隙縱橫比模型模擬彈性模量與等效介質(zhì)理論計(jì)算模量對(duì)比(圖中不同符號(hào)點(diǎn)代表 不同單一孔隙縱橫比模型的模擬模量，虛線為MT模型計(jì)算模量，數(shù)字為對(duì)應(yīng)孔隙縱橫比)(a) 體積模量； (b) 剪切模量.Fig.16 Comparison of simulated elastic modulus with calculated modulus by effective medium theory for different single pore aspect ratio model (Different symbolic points represent the simulated modulus of different single pore aspect ratio models, the dotted line represents calculated modulus by MT model, and the number represents the corresponding pore aspect ratio)(a) Bulk modulus； (b) Shear modulus.

5 結(jié)論與展望

本文從兩種不同角度出發(fā)獲取巖石孔隙結(jié)構(gòu)特征的表征與描述結(jié)果，并對(duì)比了兩種方法對(duì)孔隙結(jié)構(gòu)進(jìn)行表征的相似性.其中，基于彈性性質(zhì)并結(jié)合等效介質(zhì)理論計(jì)算等效孔隙縱橫比的方法是目前最常用的孔隙結(jié)構(gòu)描述方法，在地震以及測(cè)井?dāng)?shù)據(jù)的孔隙結(jié)構(gòu)預(yù)測(cè)上已有比較多的應(yīng)用，該方法并不是直接從孔隙空間的幾何特征中求得孔隙結(jié)構(gòu)參數(shù)，而是通過等效介質(zhì)理論將這種幾何特征間接與巖石彈性性質(zhì)聯(lián)系了起來.本研究中，借助于數(shù)字巖心可以提取儲(chǔ)層巖石的孔隙網(wǎng)絡(luò)模型，并由此出發(fā)提取具有明確地質(zhì)意義的喉道長度，作為表征孔隙結(jié)構(gòu)的幾何參數(shù).我們的研究發(fā)現(xiàn)，彈性性質(zhì)聯(lián)合等效介質(zhì)理論估算的等效孔隙縱橫比參數(shù)與直接從數(shù)字巖心提取的喉道長度這一幾何參數(shù)，二者具有很好的相關(guān)性，喉道長度長對(duì)應(yīng)于等效孔隙縱橫比小的扁平裂縫型孔隙，反之對(duì)應(yīng)于等效孔隙縱橫比大的球形孔洞型孔隙.本研究在一定程度上強(qiáng)有力地證明，彈性性質(zhì)聯(lián)合等效介質(zhì)理論這種間接反演等效孔隙縱橫比的思路，對(duì)形狀復(fù)雜的孔隙結(jié)構(gòu)進(jìn)行定量化表征具有準(zhǔn)確性與合理性.相比在地震或測(cè)井孔隙結(jié)構(gòu)預(yù)測(cè)時(shí)，僅使用成像測(cè)井、少量巖心薄片或橫波預(yù)測(cè)誤差等來驗(yàn)證等效孔隙縱橫比對(duì)孔隙結(jié)構(gòu)描述的效果，本文提供了更加定量化、更加直接的證據(jù).此外，由于等效孔隙縱橫比參數(shù)與喉道長度之間良好的相關(guān)性，而喉道特征通常是描述儲(chǔ)層滲流特性的重要參數(shù)，這使得將基于測(cè)井?dāng)?shù)據(jù)或地震數(shù)據(jù)預(yù)測(cè)的等效孔隙縱橫比屬性體轉(zhuǎn)化為更具開發(fā)價(jià)值的滲流屬性體成為可能.

附錄A

對(duì)于孔隙縱橫比α小于1的極化因子，可以表示為(Mavko et al.，1998)：

式中，(Km,Gm)為骨架的彈性模量；(Ki,Gi)為孔隙包含物的彈性模量；α為孔隙縱橫比.

附錄B

包含N種孔隙的DEM模型可以表示為微分方程式(B1)(Berryman et al.，2002)：

(B1)

圖B1與圖B2中每個(gè)數(shù)據(jù)點(diǎn)的顏色代表該點(diǎn)孔隙縱橫比，顏色越偏暖色值表示縱橫比值越大，反之縱橫比值越小，圖中藍(lán)色虛線表示理論計(jì)算的不同孔隙縱橫比時(shí)縱波速度隨孔隙度變化曲線.對(duì)比發(fā)現(xiàn)不同理論計(jì)算的結(jié)果有所區(qū)別，DEM模型計(jì)算的孔隙縱橫比較MT模型偏大一些，但是不同巖性孔隙縱橫比的分布特征是一致的.在碳酸鹽巖樣品中DEM模型計(jì)算的孔隙縱橫比分布在0.04～0.5范圍內(nèi)，MT模型計(jì)算的縱橫比分布在0.03～0.4范圍內(nèi)，不同模型的孔隙縱橫比分布區(qū)間都比較寬(見頻率直方圖B3a、b中紅色所示)，同一孔隙度下都表現(xiàn)為速度越快孔隙縱橫比越大，總體上所有數(shù)據(jù)點(diǎn)可以根據(jù)顏色分為下部的冷色值與上部的暖色值兩類；砂巖樣品DEM模型計(jì)算的縱橫比在0.3～0.5之間，MT模型則在0.2～0.5之間，不同模型孔隙縱橫比的分布區(qū)間都要小于碳酸鹽巖的分布區(qū)間且縱橫比整體較大(見頻率直方圖B3a、b中藍(lán)色所示)，從顏色上所有的數(shù)據(jù)點(diǎn)都為暖色，沒有明顯的冷色值.以上對(duì)比表明即使使用的等效介質(zhì)理論不同，也不會(huì)對(duì)孔隙縱橫比的計(jì)算及最后對(duì)孔隙結(jié)構(gòu)的表征產(chǎn)生顯著影響.

圖B1 縱波速度-孔隙度交會(huì)圖中碳酸鹽巖數(shù)字巖心等效孔隙縱橫比分布(色標(biāo)表示孔隙縱橫比)(a) 基于DEM模型； (b) 基于MT模型.Fig.B1 Distribution of effective pore aspect ratio of carbonate digital core in cross-plot between P-velocity and porosity (color represent pore aspect ratio)(a) Calculated by DEM model； (b) Calculated by MT model.

圖B2 縱波速度-孔隙度交會(huì)圖中砂巖數(shù)字巖心等效孔隙縱橫比分布(色標(biāo)表示孔隙縱橫比)(a) 基于DEM模型； (b) 基于MT模型.Fig.B2 Distribution of effective pore aspect ratio of sandstone digital core in cross-plot between P-velocity and porosity (color represent pore aspect ratio)(a) Calculated by DEM model； (b) Calculated by MT model.

圖B3 砂巖與碳酸鹽巖等效孔隙縱橫比頻率分布直方圖(a) MT模型計(jì)算； (b) DEM模型計(jì)算.Fig.B3 Histogram of effective pore aspect ratio in sandstone and carbonate(a) Calculated by MT model； (b) Calculated by DEM model.