一種簡化的二維規(guī)則多邊形孔隙巖石物理模型

劉致水 包乾宗 劉俊州 時 磊

(①長安大學地質(zhì)工程與測繪學院，陜西西安 710054； ②中國石化石油勘探開發(fā)研究院，北京 100083)

0 引言

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

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

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

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

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

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

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

1 基本原理

1.1 Kachanove巖石物理模型

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

(1)

(2)

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

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

1.2 簡化的多邊形孔隙巖石物理模型

引入一個等效孔隙形狀因子

g=2(h3-h2)-1

(3)

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

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

(4)

(5)

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

經(jīng)過數(shù)學運算，得到

(6)

(7)

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

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

(8)

(9)

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

(10)

gφ+1≠0

(11)

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

(12)

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

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

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

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

2 模型測試與應用

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

(13)

(14)

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

2.1 初始和簡化的Kachanove 2D模型對比

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

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

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

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

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

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

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

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

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

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

2.2 實驗數(shù)據(jù)測試

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

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

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

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

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

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

圖6 預測縱波(a)、橫波(b)速度與實測速度交會圖

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

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

由縱波速度求得的孔隙形狀因子如圖7所示?？梢姷刃Э紫缎螤钜蜃与S孔隙度增大呈“先降后略增”趨勢，且隨黏土含量的增大而增大。需要說明的是，由縱波速度雖可求取每個樣點的孔隙形狀因子，但由于測量、測井解釋過程中的誤差，以及孔隙流體、礦物等參數(shù)影響，所求得的孔隙形狀因子并不一定能精確表征孔隙的具體形狀，其中包含了由其他因素造成的誤差，但通過該參數(shù)的自適應可幫助研究者獲得接近實測數(shù)據(jù)的橫波速度。

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

2.3 測井數(shù)據(jù)測試

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

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

3 結(jié)論

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

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

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