基于等效介質(zhì)理論的孔隙縱橫比分布反演

歐陽芳， 趙建國， 李智， 肖增佳， 賀艷曉， 趙皓， 任靜

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

0 引言

巖石的彈性性質(zhì)受孔隙度、孔隙結(jié)構(gòu)、礦物成分、及流體性質(zhì)等多種因素的綜合影響，其中孔隙結(jié)構(gòu)是影響巖石特性的重要因素.由于巖石孔隙空間分布十分復(fù)雜，孔隙結(jié)構(gòu)難以定量表征，因此需要借助理想的幾何形狀對復(fù)雜的孔隙結(jié)構(gòu)進(jìn)行近似處理.目前，最常用的孔隙結(jié)構(gòu)表征方法是將巖石中的孔隙包含物近似為橢球體，并采用橢球體的短長軸之比(即孔隙縱橫比)作為描述孔隙結(jié)構(gòu)特征的主要參數(shù).為了研究橢球狀孔隙包含物的彈性性質(zhì)，Eshelby(1957)基于彈性理論和假想實(shí)驗(yàn)建立了橢球包含物彈性質(zhì)及其幾何形狀與巖石基質(zhì)彈性模量之間的內(nèi)在聯(lián)系，并給出了橢球包含物的彈性場表達(dá)式.在Eshelby理論的基礎(chǔ)上，Wu(1966)與Dunn和Ledbetter(1995)進(jìn)一步推導(dǎo)了橢球孔隙壓縮系數(shù)P和剪切柔度Q的解析表達(dá)式，其中P和Q又稱為孔隙幾何因數(shù)(Mavko et al., 2009)，它們均為孔隙包含物彈性模量、孔隙縱橫比、以及其所在背景基質(zhì)彈性模量的函數(shù)，是描述橢球孔隙結(jié)構(gòu)特征的重要參數(shù).這一系列工作為基于橢球孔隙包含物假設(shè)的各種巖石物理模型的建立和發(fā)展奠定了重要理論基礎(chǔ).

等效介質(zhì)理論是聯(lián)系巖石微觀孔隙結(jié)構(gòu)與巖石等效彈性性質(zhì)的重要橋梁.目前，運(yùn)用較為廣泛的等效介質(zhì)理論包括：Kuster-Toks?z(KT)模型(Kuster and Toks?z，1974)、Mori-Tanaka(MT)模型(Mori and Tanaka，1973；Benveniste，1987； Ferrari，1994)、自相容近似(Self-Consistent Approximation, SCA)模型(Budiansky，1965；Berryman，1980;Berryman and Wang, 1995；Wu，1966)和微分等效介質(zhì)(Differential effective medium, DEM)模型(Zimmerman，1984；Norris，1985；Berryman et al.，2002；李宏兵和張佳佳，2014).KT模型由Kuster和Toks?z(1974)在長波長假設(shè)條件下基于散射理論推導(dǎo)得到，具有計(jì)算巖石等效彈性性質(zhì)的顯式表達(dá)式.與KT模型類似，MT模型也具有顯式的解析理論公式，因此KT模型和MT模型的數(shù)值實(shí)現(xiàn)十分簡單，且更便于實(shí)際應(yīng)用.但已有研究表明，KT模型和MT模型僅適用于包含物稀疏分布的情形(Kuster and Toks?z，1974；Berryman and Bergr，1996；Mavko et al., 2009)，并且在某些條件下，KT模型和MT模型的等效彈性模量可能會超出Hashin-Shtrikman界限(Berryman，1980；Norris，1989；Ferrari，1991).Berryman和Bergr(1996)指出，只有當(dāng)巖石基質(zhì)占巖石總體積的比例高于70～80%時(shí)，KT模型和MT模型才能得到比較合理的結(jié)果.因此，在運(yùn)用KT和MT模型時(shí)需十分謹(jǐn)慎.與顯式的KT 模型和MT模型不同，SCA模型和DEM模型屬于隱式耦合模型，這些模型均涉及一系列迭代或微分運(yùn)算.盡管SCA模型和DEM模型的計(jì)算相對復(fù)雜，但它們考慮了包含物之間的彈性相互作用，從而可以將孔隙包含物推廣到含量稍高的情況(Mavko et al., 2009).在上述等效介質(zhì)模型中，巖石等效彈性模量均為孔隙縱橫比和孔隙度的非線性復(fù)雜函數(shù).由此可見，孔隙結(jié)構(gòu)參數(shù)的獲取對等效介質(zhì)理論的有效運(yùn)用十分重要.

一直以來，許多學(xué)者嘗試采用各種方式從實(shí)際數(shù)據(jù)中提取孔隙縱橫比，根據(jù)這些方法的反演結(jié)果，大致可以分為兩類：等效縱橫比反演法和縱橫比分布反演法.等效縱橫比反演法通常假設(shè)巖石僅含一組或兩組孔隙，并利用實(shí)際數(shù)據(jù)(如巖芯數(shù)據(jù)、測井?dāng)?shù)據(jù)等)反演計(jì)算各組孔隙的等效縱橫比，比如：Kumar和Han (2005) 假設(shè)巖石孔隙空間僅由兩種孔隙構(gòu)成，即背景孔和球形孔或背景孔和微裂隙，利用Wyllie時(shí)間平均方程和DEM模型，由含水飽和碳酸鹽巖的實(shí)測速度反演計(jì)算兩種孔隙的等效孔隙縱橫比及體積含量；之后， Xu和Payne(2009)在此基礎(chǔ)上提出了適用于碳酸鹽巖的Xu-Payne模型；Zhao等(2013)進(jìn)一步將該建模思路應(yīng)用于測井和地震數(shù)據(jù)的孔隙縱橫比反演；Fortin等(2007)與Adelinet等(2011)假設(shè)巖石由球形硬孔和單一縱橫比的微裂隙構(gòu)成，在各個(gè)壓力下分別對砂巖和玄武巖進(jìn)行了孔隙縱橫比反演；De Paula等(2012)提出了一種由干燥巖石縱橫波速度-壓力數(shù)據(jù)獲取中孔(即縱橫比介于0.01和0.2之間的孔隙)和軟孔(即縱橫比小于0.01的孔隙)等效孔隙縱橫比的反演方法；郭繼亮等(2017)采用單重孔隙巖石物理模型對碳酸鹽巖孔隙形態(tài)參數(shù)(或孔隙縱橫比)和孔隙度進(jìn)行了反演；李宏兵等(2019)基于單一孔隙縱橫比的DEM模型建立了復(fù)雜孔隙儲層的三維巖石物理量版.

盡管等效縱橫比法實(shí)現(xiàn)簡單且容易應(yīng)用，但球狀硬孔隙的假設(shè)通常難以滿足實(shí)際情形；此外，微裂隙對壓力的變化十分敏感，假設(shè)所有微裂隙僅具有一個(gè)“等效縱橫比”而不考慮壓力對孔隙形態(tài)的影響往往是不夠的.與等效縱橫比反演不同，縱橫比分布反演法認(rèn)為巖石的孔隙空間是由一系列縱橫比不同的孔隙構(gòu)成的，因而更符合實(shí)際條件.早在20世紀(jì)60年代，Walsh(1965)便指出，巖石縱橫波速度的壓力依賴性是由巖石中孔隙的閉合引起的；不久，Morlier(1971)和Toks?z等(1976)根據(jù)速度與壓力之間的變化關(guān)系，提出了孔隙縱橫比分布的反演思路；之后，Cheng和Toks?z(1979)、Zimmerman(1991)與David和Zimmerman(2012)進(jìn)一步將孔隙縱橫比分布反演法與等效介質(zhì)理論結(jié)合起來.在這些孔隙縱橫比反演方法中，David和Zimmerman(2012)的孔隙結(jié)構(gòu)模型(簡稱D-Z模型)最為經(jīng)典，該模型假設(shè)巖石由固體礦物基質(zhì)、一組縱橫比相等的硬孔隙以及多組縱橫比不等的微裂隙構(gòu)成，并認(rèn)為固體礦物基質(zhì)和硬孔隙均不受壓力影響，然后以累積裂隙密度為橋梁，借助等效介質(zhì)理論將巖石彈性模量和孔隙縱橫比聯(lián)系起來；并在此基礎(chǔ)上，利用超聲縱橫波速度的壓力依賴性反演巖石硬孔隙和各組微裂隙的孔隙縱橫比及孔隙度.盡管D-Z模型能夠獲得巖石的完整孔隙縱橫比分布，但在該模型中，多重孔隙巖石的累積裂隙密度是由單重孔隙條件下的裂隙密度公式近似計(jì)算得到的，而這種近似處理導(dǎo)致D-Z模型在許多情況下無法獲得良好的反演精度.此外，D-Z模型僅考慮了DEM和MT理論，而在許多情況下，KT和SCA等其他等效介質(zhì)理論的應(yīng)用也十分廣泛.

為此，本文提出一種基于虛擬降壓的孔隙縱橫比分布反演方法，將基于DEM和MT的經(jīng)典D-Z模型推廣到KT和SCA中，結(jié)合四種等效介質(zhì)理論建立了一套完整的反演流程，該流程通過模擬假想降壓過程實(shí)現(xiàn)累積裂隙密度的準(zhǔn)確反演，從而有效克服經(jīng)典D-Z模型的精度不足.本文首先簡要回顧了四種等效介質(zhì)模型的基本理論；然后，根據(jù)壓力對孔隙形態(tài)的影響規(guī)律，推導(dǎo)了微裂隙孔隙縱橫比和裂隙密度的計(jì)算公式；在此基礎(chǔ)上，給出了基于虛擬降壓的孔隙縱橫比分布反演流程；最后，應(yīng)用實(shí)際巖芯數(shù)據(jù)對該反演方法進(jìn)行了測試和分析.

1 等效介質(zhì)理論的基本公式

等效介質(zhì)理論是聯(lián)系微觀孔隙結(jié)構(gòu)與巖石彈性性質(zhì)的重要工具.本文假設(shè)巖石孔隙形狀均為扁圓橢球體狀(見圖1，圖中α表示孔隙縱橫比，a和c分別為橢球孔的長軸和短軸)，并利用KT模型、MT模型、DEM模型和SCA模型四種不同的等效介質(zhì)理論實(shí)現(xiàn)硬孔隙和微裂隙孔隙縱橫比的反演計(jì)算.這里我們首先回顧四種等效介質(zhì)理論的基本公式.

圖1 扁圓橢球狀孔隙Fig.1 The oblate-ellipsoid pore

1.1 KT模型

基于散射理論，Kuster和Toks?z(1974)推導(dǎo)了多重孔隙巖石的Kuster-Toksoz表達(dá)式，其一種推廣形式可以表示為

(1)

以及

(2)

對于單重孔隙的干燥巖石，即N=1且Ki=Gi=0，式(1)可以簡化為

(3)

1.2 MT模型

根據(jù)Mori-Tanaka理論(Mori and Tanaka, 1973)，多重孔隙巖石的等效彈性模量滿足如下方程(Berryman and Berge, 1996)：

(4)

如果巖石僅由體積百分比分別為1-φ和φ的礦物基質(zhì)和單重孔隙構(gòu)成，則式(4)可以退化為

(5)

式中，(Km,Gm)為巖石礦物基質(zhì)的體積模量和剪切模量；(Ki,Gi)為孔隙包含物的體積模量和剪切模量.上述過程運(yùn)用了關(guān)系Pmm=Qmm=1，此式表明當(dāng)孔隙包含物材料與主相材料相同時(shí)，幾何因數(shù)為1.在干燥條件下，式(5)還可以進(jìn)一步簡化為

(6)

1.3 DEM模型

在多重孔隙情形下，Norris形式的DEM公式可以寫成(Norris, 1985; 李宏兵和張佳佳, 2014)

(7)

上述微分方程滿足初始條件：

(8)

在干燥條件下，上述耦合微分方程具有解析近似式，即(李宏兵和張佳佳, 2014)

(9)

式中，a,b,S0i,S1i,S2i和S3i的具體表達(dá)式詳見附錄B.

1.4 SCA模型

Berryman(1980; Berryman and Wang,1995)給出了N相混合物的自相容近似(SCA或CPA)的一般形式：

(10)

對于由礦物基質(zhì)和單重孔隙構(gòu)成的兩相混合物，Wu(1966)給出了如下自相容模量的計(jì)算公式：

(11)

在干燥條件下，上式變?yōu)?/p>

(12)

這里需要注意的是，盡管式(10)和式(11)均為SCA模型公式，但兩者來源卻完全不同.通常，如果一個(gè)模型既能由準(zhǔn)靜態(tài)分析得到，也能通過散射近似獲得，則稱該模型為“物理可實(shí)現(xiàn)”，否則稱“物理不可實(shí)現(xiàn)”.在上述公式中，Berryman提出的SCA模型為“物理可實(shí)現(xiàn)”模型，而Wu提出的SCA模型為“物理不可實(shí)現(xiàn)”模型.事實(shí)上，KT模型和MT模型均為“物理不可實(shí)現(xiàn)”模型，而DEM模型為“物理可實(shí)現(xiàn)”模型(Berryman and Bergr, 1996).雖然“物理不可實(shí)現(xiàn)”模型無法適用于所有類型的孔隙形狀，但對于許多形狀(如球形)，它們?nèi)阅芙o出十分準(zhǔn)確的估計(jì)結(jié)果.

2 微裂隙孔隙縱橫比及裂隙密度

圖2 受壓后巖石內(nèi)部硬孔隙和微裂隙的形態(tài)變化示意圖Fig.2 Schematic diagram of stiff pores and micro-cracks in rock at different hydrostatic pressures

(13)

2.1 孔隙縱橫比的理論公式

下面我們從微裂隙隨有效壓力的變化規(guī)律出發(fā)，推導(dǎo)微裂隙孔隙縱橫比的計(jì)算公式.考慮壓力P下某一個(gè)縱橫比為α的干燥硬幣形微裂隙(即近似于孔隙縱橫比很小的扁橢球孔隙)，由壓縮定律可知，其相對體積變化為

dVc(P)/Vc(P)=-Cc(P,α)dP,

(14)

式中，Cc(P,α)表示壓力P下縱橫比為α的微裂隙的壓縮系數(shù)；Vc(P)表示壓力P下微裂隙的體積，dVc(P)為相應(yīng)的體積變化量.在微裂隙的閉合過程中，通常假設(shè)其長軸a保持不變，即dVc/Vc=dα/α.因此，式(14)變?yōu)?/p>

dα=-[αCc(P,α)]dP,

(15)

當(dāng)α很小(即α→0)時(shí)，有(Zimmerman, 1991)

(16)

即αCc(P,α)與α無關(guān).式中，K和ν分別為體積模量和泊松比.

將式(15)從0到P進(jìn)行積分，得到

(17)

式中，α0表示零壓力下的微裂隙初始縱橫比；第二項(xiàng)表示壓力由0變化至P的過程中微裂隙縱橫比的變化量.此外，硬幣形微裂隙的縱橫比通常較小(<0.01)，因此第二項(xiàng)幾乎是與α無關(guān)的(參見式(16))，這便意味著，在相同的壓力變化區(qū)間內(nèi)，所有微裂隙的縱橫比變化量是相同的.

假設(shè)該微裂隙的閉合壓力p等于P，即α(P)=0，則式(17)變?yōu)?/p>

(18)

式中，α0(p)表示閉合壓力為p的微裂隙的初始孔隙縱橫比，在后文中，閉合壓力均用p表示，以與P相區(qū)分.上式表明，當(dāng)壓力由0增至該微裂隙的閉合壓力p時(shí)，其縱橫比變化量剛好等于其初始縱橫比大小，即當(dāng)施加的有效壓力P等于微裂隙的閉合壓力p時(shí)，微裂隙剛好完全閉合.注意，利用式(18)為計(jì)算單一微裂隙的孔隙縱橫比公式，需要已知該微裂隙的壓縮系數(shù)Cc，而此參數(shù)通常難以測量，因此，我們還需要進(jìn)一步將式(18)與實(shí)驗(yàn)可測量的物理量(如巖石壓縮系數(shù))聯(lián)系起來.

設(shè)巖石的總孔隙空間Ωp由N組縱橫比為αi的微裂隙組成，其中i=1,2,…,N，取Cp(P)為壓力P下總孔隙空間的壓縮系數(shù)，于是有

dΩp(P)/Ωp(P)=-Cp(P)dP,

Ωp(P)=∑idVc(P,αi).

(19)

聯(lián)立式(14)和式(19)，得到

φc(P)Cp(P)=∑iφi(P)Cc(P,αi),

(20)

式中，φc=Ωp/Ω為微裂隙的總孔隙度，φi=Vc(αi)/Ω為第i組微裂隙的孔隙度，其中Ω表示巖石體積.由于巖石中充滿了各種各樣孔隙縱橫比的微裂隙，因此，式(20)中的求和符號可以轉(zhuǎn)換為積分符號，且積分區(qū)間應(yīng)包含壓力P下保持開孔的所有微裂隙的孔隙縱橫比.為此，不妨令

φi=c(α)dα,

(21)

則式(20)變?yōu)?/p>

(22)

式中，α表示保持張開狀態(tài)微裂隙的最小縱橫比；c(α)稱為孔隙度分布函數(shù)，根據(jù)孔隙度與裂隙密度Γ之間的關(guān)系，即φ=4παΓ/3，c(α)可以類似地表示為

(23)

(24)

當(dāng)微裂隙縱橫比較小時(shí)，αCc(P,α)與α幾乎無關(guān)，可以提到積分符號之外.于是，有

(25)

利用式(24)，式(18)可以重新寫成

(26)

(27)

2.2 裂隙密度的理論公式

假設(shè)巖石僅包含一組縱橫比為α=c/a的硬幣形微裂隙，其中a和c分別為微裂隙的長軸和短軸，則該組微裂隙的裂隙密度Γ可以表示為

(28)

式中，Ω為巖石體積；L為巖石中微裂隙的數(shù)量.于是，微裂隙的總孔隙度φc可以寫成

(29)

對于縱橫比較小的硬幣形微裂隙，其幾何因數(shù)P和Q具有顯式解析表達(dá)式(David and Zimmerman, 2011)：

(30)

(31)

式中，ν為巖石基質(zhì)的泊松比，而巖石基質(zhì)的選取與等效介質(zhì)理論有關(guān)，比如，對于KT模型和MT模型，ν表示背景基質(zhì)的泊松比，而對于DEM模型和SCA模型，ν則表示DEM或SCA等效介質(zhì)的泊松比.另外，由式(29)—(31)可以發(fā)現(xiàn)，φP和φQ與α無關(guān)，且均為裂隙密度Γ的函數(shù).由此可見，將φP和φQ引入等效介質(zhì)模型中便可以將巖石彈性模量與裂隙密度Γ聯(lián)系起來.

對于KT模型，我們可以將φP和φQ代入式(3)中得到裂隙密度的計(jì)算公式：

(32)

式中，Kb，Gb和νb分別表示背景基質(zhì)的體積模量、剪切模量和泊松比.注意，這里的背景基質(zhì)可以是等效礦物基質(zhì)m，也可以是其它等效介質(zhì)，其選擇視具體問題而定.

對于MT模型，將φP和φQ代入式(6)中得到

(33)

考慮到微裂隙孔隙度φc?1，我們在上述公式中運(yùn)用了近似關(guān)系1-φc≈1.

對于DEM模型，David和Zimmerman(2011)推導(dǎo)了如下裂隙密度的解析近似式：

(34)

此近似式的精度略依賴于泊松比vb，但在裂縫密度很高的情況下，其計(jì)算誤差仍小于2%(David and Zimmerman, 2011).

對于SCA模型，將φP和φQ代入式(12)中有

(35)

3 基于虛擬降壓的反演策略

3.1 微裂隙的壓力響應(yīng)描述

設(shè)巖石中所有的微裂隙在壓力Pc下剛好全部閉合，且壓力區(qū)間[0,Pc]由相等的N個(gè)子區(qū)間和N+1個(gè)壓力點(diǎn)構(gòu)成，其中壓力點(diǎn)包括0,ΔP,2ΔP,…,NΔP，此處ΔP=Pc/N表示壓力增量.假設(shè)巖石中N組微裂隙的閉合壓力剛好為p1=ΔP,p2=2ΔP,…,pN=NΔP，且ΔP足夠小時(shí)，各壓力子區(qū)間內(nèi)的孔隙縱橫比可以認(rèn)為是恒定的.當(dāng)有效壓力Pe∈[pi,pi+1)時(shí)，巖石中初始孔隙縱橫比小于α0(pi+1)的所有微裂隙全部閉合，而α0>α0(pi)的微裂隙仍處于張開狀態(tài).由式(16)—(17)可知，當(dāng)壓力增加ΔP時(shí)，所有微裂隙的縱橫比變化量是相同的，這意味著，在壓力Pe下巖石中仍處于開孔狀態(tài)的各組微裂隙的孔隙縱橫比均應(yīng)相對于其初始值減小α0(pi).此時(shí)，各組開孔微裂隙的孔隙縱橫比應(yīng)為

α0(pi+1)-α0(pi),α0(pi+2)-α0(pi),…,α0(pN)-α0(pi),

(36)

圖3展示了不同壓力下巖石中各組微裂隙的開孔情況：當(dāng)實(shí)驗(yàn)壓力為0時(shí)，所有微裂隙均處于初始狀態(tài)；當(dāng)實(shí)驗(yàn)壓力增至p1，初始縱橫比小于或等于α0(p1)的微裂隙閉合，而其它微裂隙的縱橫比均減小α0(p1)；當(dāng)實(shí)驗(yàn)壓力繼續(xù)增至p2，滿足α0≤α0(p2)的微裂隙全部閉合，而仍處于開孔狀態(tài)微裂隙的縱橫比相對于其初始值減小α0(p2)；以此類推下去，當(dāng)實(shí)驗(yàn)壓力增至壓力Pc=NΔP時(shí)，巖石中的所有微裂隙全部閉合.由此可見，一旦獲得各組微裂隙的初始縱橫比(即微裂隙的初始縱橫比分布)，利用上述過程便能求出不同壓力下巖石的微裂隙縱橫比分布.

圖3 增壓過程中巖石內(nèi)各組微裂隙的形態(tài)變化示意圖Fig.3 The change of each group of micro-cracks in rock during pressure increasing

3.2 孔隙縱橫比反演流程

為了獲得微裂隙的初始縱橫比，首先需要求解微裂隙的累積裂隙密度(見式(27))，這里我們借助假想降壓過程進(jìn)行計(jì)算.注意，這里的假想降壓過程是指增壓逆過程，即對應(yīng)于圖3中由右往左、由下至上的過程，而并不等同于實(shí)際降壓過程.事實(shí)上，吳春燕等(2016)的實(shí)驗(yàn)研究表明，在實(shí)際降壓過程中巖石孔隙空間通常難以恢復(fù)到原始狀態(tài).

基于虛擬降壓的反演流程主要包括三部分，即微裂隙的累積裂隙密度反演、微裂隙孔隙縱橫比分布計(jì)算以及硬孔縱橫比反演，工作流程見圖4.

圖4 基于虛擬降壓的孔隙縱橫比反演流程圖Fig.4 Inversion workflow of aspect ratio distribution based on a thought unloading method

(1)累積裂隙密度反演

根據(jù)假想降壓思路，我們將降壓過程分為N個(gè)變化狀態(tài)，即pN→pN-1,…,p2→p1,p1→0，其中pN等于最高實(shí)驗(yàn)壓力Pc.在每個(gè)變化狀態(tài)，干燥巖石均可以視為單重孔隙結(jié)構(gòu).這意味著，我們可以直接運(yùn)用第2.2節(jié)推導(dǎo)的單重孔隙巖石的等效介質(zhì)理論公式(32)—(35)計(jì)算微裂隙的裂隙密度，具體過程如下：

④ 利用各組微裂隙的裂隙密度計(jì)算累積裂隙密度：

(37)

在上述過程中，我們自然而然地認(rèn)為實(shí)驗(yàn)壓力點(diǎn)是等間隔分布的.然而，在很多情況下，實(shí)際數(shù)據(jù)并不滿足這些要求，此時(shí)我們可以利用關(guān)系式(13)對實(shí)驗(yàn)數(shù)據(jù)點(diǎn)進(jìn)行擬合，然后利用擴(kuò)展后的數(shù)據(jù)完成上述過程.

(38)

(2)孔隙縱橫比分布函數(shù)計(jì)算

(39)

(40)

(41)

φi(α0)=c(α0)dα0.

(42)

(43)

由此可見，對于α0>α0(pi)的微裂隙，其壓力Pe下的裂隙密度分布函數(shù)與初始分布函數(shù)γ(α0)相同，其原因在于：裂隙密度只與長軸a有關(guān)(見式(28))，而a不隨壓力變化，這意味著在微裂隙沒有完全閉合的情況下，其裂隙密度不隨外部壓力而改變，而只有當(dāng)微裂隙完全閉合后，該微裂隙對巖石裂隙密度的貢獻(xiàn)方降為零.

一般而言，壓力對巖石體積Ω的影響可以忽略不計(jì)，因此由微裂隙孔隙度的定義φi=4πa3α/(3Ω)可知，孔隙度φi與縱橫比α的比值等于常量.于是，在有效壓力Pe下，φi(α)可由初始孔隙縱橫比α0和初始孔隙度分布函數(shù)φi(α0)表示為

φi(α)=φi(α0)α/α0,

(44)

此外，聯(lián)立式(41)、式(42)和式(44)還可以求出有效壓力Pe下的孔隙度分布函數(shù)c(α)以及累積孔隙度分布函數(shù)C(α).

(3)硬孔縱橫比反演

硬孔縱橫比反演是基于高壓狀態(tài)下的實(shí)驗(yàn)數(shù)據(jù)實(shí)現(xiàn)的.在高壓條件下，干燥巖石可以視為由等效礦物基質(zhì)和單重硬孔構(gòu)成的兩相混合物(見圖2d)，此時(shí)利用等效介質(zhì)理論便可以從高壓數(shù)據(jù)中提取出硬孔縱橫比.圖5給出了硬孔縱橫比的反演流程，具體步驟如下：

圖5 硬孔縱橫比反演流程圖Fig.5 Inversion workflow for aspect ratio of the stiff pores

① 獲取巖芯樣品的基本信息，包括樣品的總孔隙度φ、礦物組分及其含量、干燥巖石的密度ρD以及各有效壓力下的超聲縱橫波速度VP,D(P)和VS,D(P)，并根據(jù)反演得到的微裂隙的總孔隙度φc，計(jì)算硬孔的孔隙度φs=φ-φc.對于φc?φs的情況，也可以由φs≈φ來估算硬孔的孔隙度；

② 利用Voigt-Reuss-Hill平均公式(簡稱V-R-H平均)計(jì)算等效礦物基質(zhì)模量(Km,Gm)，即(Mavko et al., 2009)

(45)

式中，fi為各礦物組分的體積含量；Mi表示各礦物組分的體積模量或剪切模量.

(46)

4 實(shí)際應(yīng)用分析

為了分析上述反演流程的可靠性和適用性，我們分別選取了具有代表性的2塊砂巖樣品S1和S2與2塊致密碳酸鹽巖樣品C1和C2，并采用基于虛擬降壓的反演方法對其孔隙結(jié)構(gòu)參數(shù)進(jìn)行了提取.表1列出了這些巖芯樣品的基本參數(shù)和數(shù)據(jù)來源，其中樣品S1為純石英砂巖，孔隙度為4%，巖芯數(shù)據(jù)來自于文獻(xiàn)David (2012)；樣品S2為泥質(zhì)砂巖，孔隙度為23.9%，主要礦物成分為石英(66%)、長石(17%)、粘土(15%)以及極少量的黃鐵礦和菱鐵礦，其孔隙結(jié)構(gòu)特征見圖6(a—b)；樣品C1和C2為純凈的致密碳酸鹽巖，孔隙度分別為3%和0.6%，其礦物基質(zhì)均由100%的方解石構(gòu)成，兩塊碳酸鹽巖樣品均為致密灰?guī)r，其中樣品C1主要發(fā)育晶間孔和粒內(nèi)孔，而樣品C2具有較強(qiáng)的非均質(zhì)特征，孔隙以局部發(fā)育的平直微裂縫為主，孔隙結(jié)構(gòu)特征見圖6(c—f).

表1 砂巖和碳酸鹽巖樣品的基本參數(shù)Table 1 Physical properties of the sandstone and carbonate samples

圖7 干燥巖芯樣品的縱橫波速度-壓力曲線(a) 純凈砂巖樣品S1; (b) 泥質(zhì)砂巖樣品S2； (c) 孔洞型致密碳酸鹽巖樣品C1； (d) 裂縫型致密碳酸鹽巖樣品C2.Fig.7 Velocity versus confining pressure for dry rock samples(a) Pure sandstone sample S1； (b) Shaley sandstone sample S2; (c) Tight carbonate sample C1 with inter-granular pores； (d) Tight carbonate sample C2 with cracks.

圖7給出了不同實(shí)驗(yàn)壓力下干燥砂巖樣品S1，S2與碳酸鹽巖樣品C1，C2的縱橫波速度超聲實(shí)驗(yàn)測量結(jié)果.由圖7可見，砂巖樣品和碳酸鹽巖樣品對壓力的彈性響應(yīng)存在顯著差異.當(dāng)有效壓力由0 MPa增至20 MPa時(shí)，砂巖樣品的縱橫波速度迅速增大；但隨著有效壓力的繼續(xù)增大，其縱橫波速度的增長趨勢逐漸減緩，達(dá)到最大值后趨于平穩(wěn).這表明砂巖樣品中存在大量小縱橫比的微裂隙(如粒間孔)以及少量縱橫比較大的微裂隙，由于小縱橫比微裂隙的閉合壓力較小且數(shù)量極多，所以當(dāng)壓力略微增加時(shí)大量的微裂隙閉合，從而使巖石的縱橫波速度迅速增加.從圖7a—b可以推測出，砂巖樣品中的大部分微裂隙在有效壓力到達(dá)20 MPa時(shí)便已全部閉合，而僅剩余少數(shù)縱橫比較大的微孔處于開孔狀態(tài)；當(dāng)壓力繼續(xù)增大并超過剩余微孔的閉合壓力后，這些縱橫比較大的孔隙才相繼閉合，但由于其數(shù)量較少，巖石速度隨壓力的增加趨勢也逐漸減緩.對于致密碳酸鹽巖樣品而言，縱橫波速度隨壓力的變化則相對平緩(見圖7c—d).這意味著，與小縱橫比微裂隙占主導(dǎo)的砂巖樣品相比，碳酸鹽巖樣品中各種微裂隙的數(shù)量差異相對較小.因此，在有效壓力較低的范圍內(nèi)，其縱橫波速度隨壓力的增加趨勢也相對平緩.

(47)

由表2可見，對于純凈砂巖樣品S1和非均質(zhì)程度較低的孔洞型碳酸鹽巖樣品C1，KT、MT、DEM和SCA模型均能取得很好的反演效果，其反演誤差均小于0.5%，且這些模型估計(jì)的硬孔隙縱橫比十分接近(樣品S1：～0.4，樣品C1：～0.18).對于泥質(zhì)砂巖樣品S2，四種模型的反演精度略低，其中DEM模型的效果最佳，其反演誤差為2.05%，KT模型的效果最差，其反演誤差為4.31%.由此可見，泥質(zhì)及其他摻雜礦物對砂巖硬孔縱橫比的反演精度影響較大，其原因在于：巖石彈性模量的主要貢獻(xiàn)來源于礦物成分，而對于包含有多種礦物組分的非純凈砂巖，V-R-H平均僅能提供一個(gè)粗略的等效礦物基質(zhì)模量估計(jì)值，因此存在一定的誤差.另外，對于裂縫型碳酸鹽巖樣品C2，四種等效介質(zhì)理論的反演結(jié)果均不理想，最大反演誤差接近5%.由于C2為純凈碳酸鹽巖樣品，礦物組分的影響幾乎可以忽略，因此，該誤差主要來源于碳酸鹽巖的復(fù)雜孔隙結(jié)構(gòu)分布.由圖6f可見，樣品C2存在局部發(fā)育的裂縫，具有較強(qiáng)的非均質(zhì)性，而本文采用的反演算法理論上更適合非均質(zhì)性較弱的巖石.

表2 基于不同等效介質(zhì)理論反演的硬孔隙縱橫比Table 2 Inversion results of the characteristic aspect ratio for the stiff pore obtained from different effective medium theories

圖8 基于4種等效介質(zhì)理論反演的不同有效壓力下砂巖(a—b)和碳酸鹽巖(c—d)樣品的累積裂隙密度分布Fig.8 Cumulative crack density distribution as a function of the aspect ratio calculated by different effective medium theories at different pressures for sandstone (a—b) and carbonate samples (c—d)

圖9 基于4種等效介質(zhì)理論反演的不同有效壓力下砂巖(a—b)和碳酸鹽巖(c—d)樣品的微裂隙孔隙度Fig.9 Crack porosity as a function of the aspect ratio calculated by different effective medium theories at different pressures for sandstone (a—b) and carbonate samples (c—d)

下面我們從數(shù)據(jù)統(tǒng)計(jì)的角度分析四種等效介質(zhì)理論的反演效果.我們以31塊干燥泥質(zhì)砂巖和40塊干燥碳酸鹽巖的超聲縱橫波速度作為輸入數(shù)據(jù)，其中砂巖樣品的主要礦物成分為石英、長石、方解石、白云石和粘土，孔隙度分布范圍為3.5%～30%；碳酸鹽巖樣品主要由方解石構(gòu)成，其孔隙度分布范圍為0.3%～3% .圖10給出了四種理論的模擬結(jié)果及超聲實(shí)驗(yàn)數(shù)據(jù).從圖中可以看出，無論是泥質(zhì)砂巖還是碳酸鹽巖，四種等效介質(zhì)理論的模擬速度與實(shí)測速度均具有很好的一致性.結(jié)合圖8和圖9可見，這些理論的反演精度和反演結(jié)果均十分接近.由此可見，本文提出的孔隙結(jié)構(gòu)反演方法并不十分依賴于等效介質(zhì)理論的選擇，四種等效介質(zhì)理論均能作為反演孔隙結(jié)構(gòu)參數(shù)的有效工具.但是，在計(jì)算效率方面，由于KT和MT模型具有顯式的理論表達(dá)式，其計(jì)算速度比隱式的DEM和SCA模型更快.

為了對比本文方法與現(xiàn)有方法，我們采用基于虛擬降壓的反演方法(見第4.2節(jié))和D-Z方法(Zimmerman, 1999; David and Zimmerman, 2012;鄧?yán)^新等, 2015)計(jì)算了不同有效壓力下干燥巖石的彈性模量，以此評估兩種方法的反演精度.本文方法與D-Z方法的主要區(qū)別在于，D-Z方法直接采用單重孔隙巖石的裂隙密度公式(式(33)和式(34))近似反演多重孔隙巖石的累積裂隙密度；而本文方法則借助假想降壓過程實(shí)現(xiàn)多重孔隙巖石的累積裂隙密度計(jì)算.以DEM模型為例，圖11和圖12分別從單塊巖芯與統(tǒng)計(jì)分析的角度對比了兩種方法的反演結(jié)果.由圖11可見，對于純凈砂巖樣品S1，本文方法和D-Z方法均能取得很好的反演效果；但對于泥質(zhì)砂巖和碳酸鹽巖，由D-Z方法計(jì)算的體積模量與實(shí)測數(shù)據(jù)存在較大誤差，且計(jì)算誤差隨著有效壓力的減小而增大，這表明D-Z模型所采用的近似做法(即直接利用單重孔隙公式(33)—(34)近似計(jì)算多重孔隙巖石的累積裂隙密度)存在一定局限性.另一方面，從圖11—12中可以看出，本文方法相較于D-Z方法具有更高的精度，無論是砂巖還是碳酸鹽巖，由本文方法計(jì)算得到的體積模量和剪切模量與實(shí)際數(shù)據(jù)吻合更好.由此可見，本文方法比D-Z方法更具優(yōu)勢，能夠有效克服經(jīng)典D-Z反演方法的精度不足.

圖10 基于不同等效介質(zhì)理論的模型反演結(jié)果與實(shí)際測量結(jié)果對比Fig.10 Comparison of the measured velocities and the results inverted using different effective medium theories

圖11 本文方法和D-Z方法基于DEM理論反演的干燥巖石彈性模量對比(a) 純凈砂巖樣品S1; (b) 泥質(zhì)砂巖樣品S2； (c) 孔洞型致密碳酸鹽巖樣品C1； (d) 裂縫型致密碳酸鹽巖樣品C2.Fig.11 Elastic Moduli of dry rock calculated with DEM by using our method and D-Z method(a) Pure sandstone sample S1； (b) Shaley sandstone sample S2; (c) Tight carbonate sample C1 with inter-granular pores； (d) Tight carbonate sample C2 with cracks.

圖12 本文方法和D-Z模型反演結(jié)果對比Fig.12 Comparison of our method and D-Z method

5 結(jié)論

本文對經(jīng)典D-Z方法進(jìn)行了改進(jìn)和完善，并將其推廣至KT、MT、DEM和SCA四種等效介質(zhì)模型的情形，建立了一套完整的孔隙縱橫比反演流程.以一系列砂巖樣品和碳酸鹽巖樣品為例，分析了本文方法在實(shí)際巖芯孔隙縱橫比提取中的應(yīng)用效果，數(shù)值分析結(jié)果表明：基于本文方法從干燥巖芯速度-壓力數(shù)據(jù)中提取的孔隙縱橫比分布、累積裂隙密度分布和微裂隙孔隙度分布可以很好地解釋實(shí)際巖芯中微裂隙的壓力響應(yīng)規(guī)律；此外，本文方法并不十分依賴于等效介質(zhì)理論的選擇，即四種理論均能取得十分好的反演精度，且反演的孔隙結(jié)構(gòu)參數(shù)差異很小；與經(jīng)典D-Z方法相比，本文方法有效改善了D-Z方法中利用單重孔隙公式近似計(jì)算多重孔隙巖石累積裂隙密度所引起的精度不足.在實(shí)際應(yīng)用方面，本文所建立的反演流程可以為巖石物理實(shí)驗(yàn)和理論研究提供重要的孔隙結(jié)構(gòu)參數(shù)提取工具；另外，結(jié)合深度和壓力的對應(yīng)關(guān)系，此方法也有望應(yīng)用于測井?dāng)?shù)據(jù)的巖石孔隙結(jié)構(gòu)參數(shù)反演中.

附錄A 橢球狀孔隙的幾何因數(shù)

對于縱橫比小于1的橢球狀孔隙，其幾何因數(shù)P和Q可以表示為(Mavko et al.，2009)

×(F2F4)-1],

式中，

F1=1+A[1.5(fi+θi)-R(1.5fi+2.5θi-4/3)],

F2=1+A[1+1.5(fi+θi)-0.5R(3fi+5θi)]

+B(3-4R)+0.5A(A+3B)(3-4R)[fi+

F3=1+A[1-(fi+1.5θi)+R(fi+θi)],

F4=1+0.25A[fi+3θi-R(fi-θi)],

F5=A[-fi+R(fi+θi-4/3)]+Bθi(3-4R),

F6=1+A[1+fi-R(fi+θi)]

+B(1-θi)(3-4R),

F7=2+0.25A[3fi+9θi-R(3fi+5θi)]

+Bθi(3-4R),

F8=A[1-2R+0.5fi(R-1)+0.5θi(5R-3)]

+B(1-θi)(3-4R),

F9=A[(R-1)fi-Rθi]+Bθi(3-4R),

其中：

A=Gi/Gm-1,

R=(Km/Gm+4/3)-1

附錄B DEM解析近似式的系數(shù)

在DEM的解析近似式中，a,b,S0i,S1i,S2i和S3i的具體表達(dá)式分別為

S1i=(θi-fi)/(2-3fi-3θi),

S2i=4/3,

式中，

-0.2F2F4(F′4F5+F4F′5+F′6F7+F6F′7

F′1=-A(1.5fi+2.5θi-4/3)R′,

F′2={-0.5A(3fi+5θi)-4B-0.5A(A+3B)

F′3=A(fi+θi)R′,

F′4=-0.25A(fi-θi)R′,

F′5=[A(fi+θi-4/3)-4Bθi]R′,

F′6=[-A(fi+θi)-4B(1-θi)]R′,

F′7=[-0.25A(3fi+5θi)-4Bθi]R′,

F′8=[A(-2+0.5fi+2.5θi)-4B(1-θi)]R′,

F′9=[A(fi-θi)-4Bθi]R′,

R′=-R2.