劉致水 包乾宗 劉俊州 時(shí) 磊

(①長(zhǎng)安大學(xué)地質(zhì)工程與測(cè)繪學(xué)院，陜西西安 710054； ②中國(guó)石化石油勘探開(kāi)發(fā)研究院，北京 100083)

0 引言

巖石物理模型是描述儲(chǔ)層特征與彈性參數(shù)之間關(guān)系的數(shù)學(xué)函數(shù)。針對(duì)不同地質(zhì)目標(biāo)，基于等效介質(zhì)理論，人們提出多種巖石物理模型以描述巖石中礦物顆粒與一定形狀巖石孔隙之間的關(guān)系[1-4]。

目前，考慮孔隙結(jié)構(gòu)對(duì)巖石彈性參數(shù)影響的巖石物理模型主要有2D孔隙和3D孔隙模型兩種：2D孔隙形狀包括橢圓形、裂縫形、多邊形； 3D孔隙形狀有球形、橢球形、針形和裂縫形。

在2D孔隙研究中，Kuster等[2]推導(dǎo)了含多種橢圓形孔隙巖石的彈性參數(shù)公式。Zimmerman[5]系統(tǒng)研究了多種2D孔隙形狀對(duì)巖石彈性參數(shù)的影響。Mauge等[6]研究了各向同性與各向異性情況下裂縫孔隙對(duì)巖石彈性參數(shù)的影響。Thorpe等[7]在自適應(yīng)巖石物理模型框架下推算了巖石中含方向隨機(jī)的橢圓形孔隙的彈性模量公式。Jasiuk等[8-9]導(dǎo)出了含多邊形孔隙介質(zhì)的彈性模量公式。Kachanov等[4]推導(dǎo)了含2D規(guī)則多邊形孔隙巖石的巖石物理模型公式，使用3個(gè)參數(shù)描述規(guī)則形狀孔隙對(duì)巖石彈性參數(shù)的影響； 且利用Savin方法[10]求得三角形、正方形等幾種2D規(guī)則多邊形孔隙的孔隙形狀因子。

在3D孔隙研究中，Mackenzie[11]構(gòu)建了含3D球形孔隙的巖石物理模型。Zimmerman[12]推算出含橢球體孔隙巖石的剛度張量。Eshelby[13]求得含低孔隙度球形孔隙巖石的彈性張量。Berryman利用自適應(yīng)巖石物理模型分析了含隨機(jī)分布的橢球體包含物巖石的彈性張量[1]； 又在自適應(yīng)巖石物理模型基礎(chǔ)上給定了多種3D孔隙的形狀因子，最常用的包括球形、針形及裂縫形[14]。

這些經(jīng)典模型基本都只適用于礦物簡(jiǎn)單、孔隙類型單一的情況，且都僅適用于高頻條件。Xu等[15]將微分Kuster-Toks?z模型[2]、Gassmann方程[16]、Wyllie時(shí)間平均公式[17]組合起來(lái)，得到一種兼顧兩種礦物、兩種2D橢圓形孔隙、適用于低頻條件的巖石物理模型。其后眾多學(xué)者基于不同需求，利用不同的基礎(chǔ)巖石物理模型進(jìn)行組合，獲得多個(gè)在業(yè)界應(yīng)用效果良好的巖石物理模型。這些模型使用的基礎(chǔ)巖石物理模型主要為自適應(yīng)模型和微分等效介質(zhì)(DEM)理論，所考慮的孔隙形狀也基本是2D橢圓形孔隙[15，18-22]和3D橢球體孔隙[23-25]，且鮮見(jiàn)其他孔隙類型的巖石物理模型的研究與應(yīng)用。

2D橢圓形模型和3D橢球體模型應(yīng)用較多的原因在于其描述孔隙形狀變化的參數(shù)只有孔隙縱橫比，研究者可反演該參數(shù)[21，25]，并賦予其具體地質(zhì)意義。雖然Kachanov等[4]推導(dǎo)出2D規(guī)則多邊形孔隙巖石物理模型，但該模型的孔隙形狀參數(shù)有三個(gè)，不同多邊形孔隙的形狀因子具體數(shù)值不同，應(yīng)用時(shí)需根據(jù)Savin方法[10]求取孔隙形狀參數(shù)具體數(shù)值； 且?guī)r石孔隙形狀極為復(fù)雜，這幾種規(guī)則多邊形孔隙未能囊括巖石所有復(fù)雜孔隙形狀； 當(dāng)每次使用時(shí)，研究者尚不能預(yù)知巖石孔隙形狀并確定孔隙形狀因子； 三種參數(shù)不利于巖石物理研究中經(jīng)常使用的數(shù)學(xué)方法的應(yīng)用。因此，限制了該模型的推廣和應(yīng)用。

為此，引入一個(gè)等效孔隙形狀因子g，替代Kachanove 2D模型中的多個(gè)形狀因子參數(shù)，得到簡(jiǎn)化的2D規(guī)則多邊形孔隙巖石物理模型。文中厘定了等效孔隙形狀因子的理論參數(shù)范圍，并通過(guò)數(shù)值模擬、實(shí)驗(yàn)室測(cè)試及實(shí)際測(cè)井資料應(yīng)用，證實(shí)了本文方法對(duì)砂巖儲(chǔ)層的適用性。

1 基本原理

1.1 Kachanove巖石物理模型

經(jīng)典Mori-Tanaka巖石物理模型是基于平均應(yīng)力場(chǎng)近似理論推導(dǎo)的含孔隙干巖石模型。在此基礎(chǔ)上Kachanove等[4]給出含2D規(guī)則多邊形孔隙干巖石的楊氏模量Edry和泊松比σdry計(jì)算公式

(1)

(2)

式中：Em和σm是巖石基質(zhì)的楊氏模量和泊松比；φ是孔隙度；h2和h3是孔隙形狀因子，部分情況下h3與另一因子h1有關(guān)系，具體細(xì)節(jié)見(jiàn)表1。

表1是三角形、四邊形、五邊形的孔隙形狀示意圖及其對(duì)應(yīng)的形狀因子h1、h2和h3。利用Savin方法[10]求取2D規(guī)則多邊形孔隙形狀因子。在應(yīng)用該模型時(shí)，需先由巖石鏡下薄片觀察巖石的孔隙形狀，確定孔隙形狀因子，再計(jì)算巖石彈性模量。此外，該模型只能應(yīng)用于表1所列幾種2D規(guī)則多邊形孔隙，其他形狀孔隙則據(jù)Savin方法[10]推導(dǎo)，且此模型中孔隙形狀因子參數(shù)較多，難以與自適應(yīng)算法等數(shù)學(xué)方法結(jié)合，限制了該模型的可應(yīng)用性。

1.2 簡(jiǎn)化的多邊形孔隙巖石物理模型

引入一個(gè)等效孔隙形狀因子

g=2(h3-h2)-1

(3)

替代式(1)和式(2)中多個(gè)孔隙形狀因子，即簡(jiǎn)化Kachanove模型。表1列舉了幾種多邊形孔隙的g值。

將式(3)代入式(1)和式(2)，可得

(4)

(5)

由楊氏模量E、泊松比σ與體積模量K、剪切模量G之間的關(guān)系[26]

經(jīng)過(guò)數(shù)學(xué)運(yùn)算，得到

(6)

(7)

式中：Kdry、Gdry分別是干巖石體積模量和剪切模量；Km和Gm分別是巖石基質(zhì)的體積模量和剪切模量。

為明確g的取值范圍，對(duì)上兩式取g的導(dǎo)數(shù)

(8)

(9)

若式(8)和式(9)有意義，則式中分母部分不為0。即要使式(8)有意義，則g的取值需滿足

(10)

gφ+1≠0

(11)

要使式(9)有意義，則g取值需滿足式(11)和下式

(12)

因?yàn)棣铡蔥0，1]，巖石基質(zhì)泊松比σm∈(0，0.5)，經(jīng)推導(dǎo)可知：g∈(1,+∞)為式(6)定義域，g∈(-1,+∞)為式(7)定義域。在計(jì)算縱波速度時(shí)需同時(shí)使用式(6)和式(7)，此時(shí)形狀因子g∈(1,+∞)。

圖1為不同孔隙度的砂泥巖體積模量(式(6)，圖1a)及其對(duì)g的導(dǎo)數(shù)(式(8)，圖1b)隨g變化的正演曲線。圖中不同的曲線表示孔隙度在0.01～0.34區(qū)間以0.03的步長(zhǎng)變化。其中砂泥巖基質(zhì)的體積和剪切模量分別為39.0、32.8GPa； 孔隙含水，水的體積模量為2.2GPa，g在1～600區(qū)間變化。

圖1 體積模量(a)和體積模量的導(dǎo)函數(shù)(b)隨等效孔隙形狀因子的變化曲線

從該圖可見(jiàn)： ①形狀因子g越接近1，體積模量越大；隨著g增大，巖石體積模量降低。 ②g較小時(shí)，模量降低速率快； 隨著g增大，模量降低速率變慢； 當(dāng)g大于某個(gè)數(shù)時(shí)，彈性模量的變化率基本可忽略，因此實(shí)際應(yīng)用中g(shù)的取值并非無(wú)限大，而是一個(gè)有限的數(shù)，圖中當(dāng)g大于約500時(shí)，彈性模量的變化率即可忽略，故本文設(shè)定取值范圍是g∈(1,500]。 ③彈性模量隨g的變化率與孔隙度有關(guān)，當(dāng)孔隙度較小時(shí)，模量隨g的變化率較小。進(jìn)一步研究得知，剪切模量變化規(guī)律與體積模量一致。

2 模型測(cè)試與應(yīng)用

在改進(jìn)的多邊形孔隙模型中，彈性模量只與一個(gè)參數(shù)(等效孔隙形狀因子g)有關(guān)。因此，在已知孔隙度、巖石基質(zhì)彈性模量和g時(shí)，結(jié)合2D-Kachanove模型和Gassmann方程[16]，可求取飽和流體巖石的縱、橫波速度[26]

(13)

(14)

式中：VPc、VSc分別是飽和流體巖石的縱、橫波速度；Ksat、Gsat和ρsat分別是飽和流體巖石的體積模量、剪切模量和密度，且有ρsat=ρm(1-φ)+ρfφ，其中ρm、ρf分別為巖石基質(zhì)、孔隙流體的密度，對(duì)于干巖石，ρf取0。需要說(shuō)明的是，若巖石基質(zhì)礦物較復(fù)雜，可通過(guò)VRH平均公式[26]計(jì)算多種礦物混合情形下的巖石基質(zhì)體積模量和剪切模量。

2.1 初始和簡(jiǎn)化的Kachanove 2D模型對(duì)比

假設(shè)巖石基質(zhì)是含黏土的石英，其體積模量和剪切模量分別為39.0、32.8GPa，密度為2.65g/cm3； 孔隙含水，其體積模量和密度分別為2.2GPa和0.99g/cm3。利用式(1)和式(2)可計(jì)算含幾種規(guī)則形狀孔隙巖石的速度—孔隙度關(guān)系曲線(圖2)。由該圖可見(jiàn)：

(1)相同孔隙度下含不同形狀孔隙巖石的速度隨孔隙度的變化規(guī)律不同，具有圓角的凸邊孔隙速度大于具有尖角的凹邊孔隙。如在三角形孔隙中，速度大小排序?yàn)棰?②； 四邊形孔隙中，速度排序?yàn)棰?④； 五邊形孔隙中，速度排序?yàn)棰?⑥。

(2)分別針對(duì)凸邊和凹邊孔隙，對(duì)比三角形、四邊形、五邊形孔隙，可見(jiàn)五邊形孔隙>三邊形孔隙>四邊形孔隙，其中三邊形孔隙與四邊形孔隙近似。

結(jié)合表1、圖2可見(jiàn)，對(duì)于不同形狀的孔隙，當(dāng)g增大時(shí)巖石速度降低，g降低時(shí)巖石速度增高，通過(guò)g的變化，可達(dá)到代替原始模型中多個(gè)孔隙形狀因子的目的。

圖3為由改進(jìn)模型計(jì)算的巖石速度隨孔隙度、等效孔隙形狀因子變化的曲線圖。圖中巖石基質(zhì)、孔隙流體參數(shù)與圖2一致，綠色曲線部分是g在1.01～1.91之間以0.10為步長(zhǎng)進(jìn)行變化，黑色曲線部分是g在2～100之間以1為步長(zhǎng)進(jìn)行變化。對(duì)比圖2和圖3可知，改進(jìn)模型除能表達(dá)圖2中幾種孔隙類型外，g在理論定義域內(nèi)的其他取值提供了描述更多孔隙形狀的可能性，從而極大地?cái)U(kuò)大了原始模型描述孔隙類型的可能性。

圖2 含表1中幾種孔隙的巖石縱波速度(a)、橫波速度(b)隨孔隙度變化曲線

圖3 由改進(jìn)模型計(jì)算的巖石縱波速度(a)、橫波速度(b)隨孔隙度、等效孔隙形狀因子的變化曲線HS-Up和HS-Low分別指Hashin-Shtrikman上、下限

圖3中還給出了Hashin-Shtrikman上、下限[26]，可見(jiàn)g=1.61時(shí)計(jì)算的縱波速度隨孔隙度的變化曲線與Hashin-Shtrikman上限基本一致，g=100時(shí)計(jì)算的縱波速度隨孔隙度的變化曲線接近Hashin-Shtrikman下限且在其上方。當(dāng)g在1.01～100取值時(shí)，橫波速度曲線總在Hashin-Shtrikman上、下限之間。數(shù)值實(shí)驗(yàn)還發(fā)現(xiàn)，g取較大數(shù)值范圍(如1～10000)時(shí)，所得結(jié)果也基本一致。需強(qiáng)調(diào)的是，本文改進(jìn)模型與經(jīng)典DEM模型呈現(xiàn)的規(guī)律類似。

圖4為2D橢圓形孔隙的DEM模型計(jì)算的縱、橫波速度[25]與孔隙度關(guān)系圖。對(duì)比圖4與圖3可知，改進(jìn)模型通過(guò)g的變化描述孔隙結(jié)構(gòu)對(duì)速度的影響，而DEM模型則通過(guò)橢圓形孔隙的孔隙縱橫比(橢圓短軸與長(zhǎng)軸之比)α描述孔隙結(jié)構(gòu)對(duì)速度的影響，兩者效果基本一致。

圖4 通過(guò)2D橢圓形孔隙的DEM模型計(jì)算的縱波速度(a)和橫波速度(b)與孔隙度關(guān)系圖中曲線是橢圓形孔隙的孔隙縱橫比α從1.00到0.01之間按照步長(zhǎng)0.01進(jìn)行變化

2.2 實(shí)驗(yàn)數(shù)據(jù)測(cè)試

利用Han等[27]在超聲實(shí)驗(yàn)室40MPa條件下針對(duì)砂泥巖測(cè)量的75個(gè)數(shù)據(jù)點(diǎn)測(cè)試本文方法。Han數(shù)據(jù)的孔隙度范圍是0.02～0.31，黏土含量為0～0.5。將本文方法數(shù)值正演結(jié)果與實(shí)測(cè)數(shù)據(jù)進(jìn)行疊合(圖5)。所用參數(shù)如表2所示。該圖顯示，隨著孔隙度增加，巖石速度降低； 孔隙度不變時(shí)，隨著g變大，巖石速度降低； 孔隙度及形狀因子g相同時(shí)，泥巖速度低于砂巖。通過(guò)孔隙度、黏土含量、g的變化，可包括所有數(shù)據(jù)點(diǎn)，即改進(jìn)模型可解釋測(cè)量數(shù)據(jù)的孔隙度—黏土含量—速度三者之間關(guān)系。

表2 巖石組成成分的彈性參數(shù)和密度[26-28]

在巖石物理研究過(guò)程中，理想情況是通過(guò)巖心等實(shí)物觀測(cè)巖石的孔隙結(jié)構(gòu)，并將實(shí)際觀測(cè)參數(shù)代入巖石物理模型。然而，實(shí)際應(yīng)用中很難精確確定待研究巖石的真實(shí)孔隙形態(tài)，因此需將巖石物理模型與數(shù)學(xué)工具結(jié)合，采取迭代或自適應(yīng)方法求取等效孔隙結(jié)構(gòu)參數(shù)。從圖5可知，當(dāng)以較小步長(zhǎng)變化孔隙形狀因子，總能找到一條曲線穿過(guò)圖5a中每個(gè)實(shí)測(cè)樣點(diǎn)，即將該樣點(diǎn)的孔隙形狀等效為該曲線所對(duì)應(yīng)的孔隙形狀，再以此孔隙形狀因子求算橫波速度。還可看出，穿過(guò)圖5a中每個(gè)樣點(diǎn)的曲線不一定與圖5b中曲線吻合，但差距不太大，其原因是在測(cè)量縱、橫波速度時(shí)，人、儀器因素必然導(dǎo)致針對(duì)同一目標(biāo)的測(cè)量數(shù)據(jù)有誤差。

圖5 基于改進(jìn)模型的縱波速度(a)和橫波速度(b)正演結(jié)果黑線代表往純石英基質(zhì)加入0～0.35孔隙度(干孔隙空腔)的變化； 紅線代表以50%石英+50%黏土混合物作為基質(zhì)(用VRH平均公式計(jì)算混合礦物體積模量和剪切模量)，往其中加入0～0.35孔隙度(干孔隙空腔)的變化； 從上到下的黑線和紅線均代表等效孔隙形狀因子g在1.01～30.01區(qū)間以步長(zhǎng)1變化。樣點(diǎn)為一組砂巖巖心測(cè)試數(shù)據(jù)[27]，其顏色對(duì)應(yīng)不同的黏土含量

本文利用速度預(yù)測(cè)效果驗(yàn)證模型的可應(yīng)用性?，F(xiàn)實(shí)中的巖石都是非常復(fù)雜的混合物，影響巖石彈性模量的因素包括巖石的組成礦物及其彈性模量、孔隙度、孔隙流體、等效孔隙形狀因子等，而本文改進(jìn)的巖石物理模型是考慮純粹的理想情況。因此，參考Xu-White砂泥巖速度預(yù)測(cè)流程[15]的思路，采取多個(gè)模型結(jié)合以構(gòu)建等效介質(zhì)模型的速度計(jì)算流程可簡(jiǎn)述為： ①使用VRH平均公式[26]計(jì)算巖石基質(zhì)體積模量和剪切模量； ②利用式(6)和式(7)計(jì)算含干孔隙巖石的彈性模量； ③利用Wood方程[26]計(jì)算混合流體的體積模量，并用Gassmann方程[16]計(jì)算飽和流體巖石的彈性模量； ④利用式(13)計(jì)算飽和流體巖石縱波速度； ⑤利用計(jì)算的縱波速度VPc與實(shí)測(cè)縱波速度VPm構(gòu)建目標(biāo)函數(shù)ε=|VPm-VPc(g)|/VPm，并使ε取極小值，以求取等效孔隙形狀因子g； ⑥利用g最終求得橫波速度。

圖6為利用文獻(xiàn)[27]數(shù)據(jù)在縱波速度約束下得到的計(jì)算速度與實(shí)測(cè)速度的對(duì)比圖。從圖6a可知計(jì)算縱波速度與實(shí)測(cè)縱波速度一致，這是由于等效孔隙形狀因子是由縱波速度計(jì)算得到的，基于此孔隙形狀因子計(jì)算的縱波速度與實(shí)測(cè)縱波速度吻合度較高(基本是相等的)。圖6b中數(shù)據(jù)點(diǎn)均勻分布于對(duì)角線(紅線)附近，說(shuō)明預(yù)測(cè)橫波速度與實(shí)測(cè)數(shù)據(jù)吻合度較高。

圖6 預(yù)測(cè)縱波(a)、橫波(b)速度與實(shí)測(cè)速度交會(huì)圖

選取計(jì)算與實(shí)測(cè)速度之間相對(duì)誤差的平均值A(chǔ)E、均方根誤差RMSE、相關(guān)系數(shù)R2等參數(shù)[28]定量評(píng)價(jià)預(yù)測(cè)結(jié)果。由表3可見(jiàn)，計(jì)算結(jié)果與實(shí)測(cè)數(shù)據(jù)之間誤差較小，即利用改進(jìn)模型能獲得較準(zhǔn)確的砂巖儲(chǔ)層橫波速度，這說(shuō)明本文改進(jìn)模型適用于砂泥巖儲(chǔ)層。

表3 實(shí)測(cè)數(shù)據(jù)與預(yù)測(cè)數(shù)據(jù)之間的的誤差統(tǒng)計(jì)表

由縱波速度求得的孔隙形狀因子如圖7所示。可見(jiàn)等效孔隙形狀因子隨孔隙度增大呈“先降后略增”趨勢(shì)，且隨黏土含量的增大而增大。需要說(shuō)明的是，由縱波速度雖可求取每個(gè)樣點(diǎn)的孔隙形狀因子，但由于測(cè)量、測(cè)井解釋過(guò)程中的誤差，以及孔隙流體、礦物等參數(shù)影響，所求得的孔隙形狀因子并不一定能精確表征孔隙的具體形狀，其中包含了由其他因素造成的誤差，但通過(guò)該參數(shù)的自適應(yīng)可幫助研究者獲得接近實(shí)測(cè)數(shù)據(jù)的橫波速度。

圖7 等效孔隙形狀因子與孔隙度(a)和黏土含量(b)的交會(huì)圖

2.3 測(cè)井?dāng)?shù)據(jù)測(cè)試

將本文方法應(yīng)用于四川盆地西部地區(qū)三疊系沙河街組砂巖儲(chǔ)層段。圖8展示了該區(qū)A井目的層段的測(cè)井解釋、速度預(yù)測(cè)及孔隙形狀因子計(jì)算等綜合結(jié)果。利用多礦物測(cè)井解釋方法[29]估算孔隙度、礦物體積含量，解釋的礦物包括黏土、石英、方解石、長(zhǎng)石。經(jīng)過(guò)測(cè)試、標(biāo)定，計(jì)算過(guò)程中使用的礦物體積模量和剪切模量如表2所示。將礦物的體積分?jǐn)?shù)、孔隙度數(shù)據(jù)引入巖石物理模型，在縱波速度數(shù)據(jù)約束下預(yù)測(cè)橫波速度。由圖8可見(jiàn)，本文方法預(yù)測(cè)結(jié)果與實(shí)測(cè)數(shù)據(jù)吻合度較高，縱、橫波速度的相對(duì)誤差基本控制在8%以內(nèi)。表3為實(shí)測(cè)數(shù)據(jù)與預(yù)測(cè)數(shù)據(jù)之間的誤差統(tǒng)計(jì)，可見(jiàn)測(cè)井?dāng)?shù)據(jù)預(yù)測(cè)的橫波速度的誤差平均值為0.008、均方根誤差為0.137、相關(guān)系數(shù)為0.853。圖8和表3均說(shuō)明新方法在砂泥巖測(cè)井資料中的應(yīng)用效果良好。

圖8 四川盆地西部地區(qū)A井三疊系沙河街組砂巖儲(chǔ)層段速度預(yù)測(cè)與實(shí)測(cè)數(shù)據(jù)對(duì)比圖

3 結(jié)論

(1)針對(duì)規(guī)則形狀孔隙的2D-Kachanove模型孔隙形狀因子的參數(shù)多、僅能表示幾種孔隙形狀、難以與自適應(yīng)等數(shù)學(xué)算法結(jié)合等問(wèn)題，本文引入新等效孔隙形狀因子以替換原始模型的多個(gè)孔隙形狀因子，得到一個(gè)簡(jiǎn)化的2D-Kachanove巖石物理模型。等效孔隙形狀因子的理論取值范圍是(1,+∞)。

(2)數(shù)值算例表明： 隨著g增大，巖石彈性模量降低，且降低速率逐漸變小； 當(dāng)g大于某個(gè)數(shù)時(shí)，彈性模量的變化率基本可忽略，因此實(shí)際應(yīng)用中g(shù)的取值并非無(wú)限大，而是一個(gè)有限的數(shù)。相對(duì)于原始模型僅能考慮幾種孔隙形狀，改進(jìn)模型可考慮g在理論取值范圍內(nèi)的任意值，從而擴(kuò)大了原始2D-Kachanove模型的應(yīng)用范圍并使其更易于推廣。

(3)利用簡(jiǎn)化巖石物理模型對(duì)實(shí)驗(yàn)室測(cè)試數(shù)據(jù)及測(cè)井資料進(jìn)行橫波速度預(yù)測(cè)，結(jié)果證明簡(jiǎn)化模型在砂(泥)巖儲(chǔ)層中有較好應(yīng)用效果。