楊鳳英，印興耀,劉 博

(中國石油大學(xué)(華東)地球科學(xué)與技術(shù)學(xué)院，山東青島266580)

基于巖石物理理論，結(jié)合地質(zhì)、地震資料構(gòu)建巖石物理模型，用測井資料作為約束可以進(jìn)行儲層參數(shù)預(yù)測[1-3]。Gassmann[4]給出了低頻條件下飽和巖石的等效模量與巖石基質(zhì)模量、干巖石骨架模量、流體模量以及孔隙度之間的關(guān)系。巖石基質(zhì)模量可以用等效介質(zhì)理論求解[5]。干巖石骨架彈性模量的計算方法常見的有以下幾類：①實驗室?guī)r石樣品直接測量干巖石骨架模量，雖然精度高但是成本也相對昂貴，不具有普遍適用性[6]。②經(jīng)驗公式法，如：Krief等[7]基于Raymer等[8]得到的實驗數(shù)據(jù)給出了干巖石骨架模量隨孔隙度變化的關(guān)系；Nur[9]提出了臨界孔隙度概念，認(rèn)為干巖石骨架模量隨孔隙度呈線性關(guān)系；Mavko等[10]指出Nur的臨界孔隙度公式與Krief等的公式基本一致。③理論計算法，如：Eshelby[11]和Walsh[12]認(rèn)為巖石骨架模型不僅與孔隙度有關(guān)還與包含物縱橫比有關(guān)；Kuster和Toksoz[13]根據(jù)長波一階散射理論得到了含多種包含物形狀的等效模量，即Kuster-Toksoz(K-T)理論，但該理論只適用于孔隙度與孔隙縱橫比之比遠(yuǎn)小于1的情形；Pride等[14]給出了含有固結(jié)系數(shù)的巖石等效模量；Lee[15]對Pride等人公式中的剪切模量做了修改；張佳佳等[16]研究了幾種巖石骨架模型的適用性，并發(fā)現(xiàn)Pride模型比Krief模型和Nur模型適用范圍更廣，且Krief模型和Nur模型對低壓的低孔巖石不適用；微分等效介質(zhì)(DEM)理論[17-19]通過往固體礦物相中逐漸加入包含物相來模擬雙相混合物得到了耦合的微分方程組，由于方程耦合很難直接對方程積分得到精確的解析解，通常只能得到近似的數(shù)值解；Keys和Xu[20]結(jié)合K-T理論和DEM理論，在干巖石泊松比不隨孔隙度變化的假設(shè)下得到了干巖石骨架彈性模量的Keys-Xu近似公式，但用他們給出的干巖石彈性模量近似式計算得到的干巖石泊松比卻并不是常數(shù)。

在前人研究的基礎(chǔ)上，我們將含有可變參數(shù)的干巖石骨架等效模量代入耦合的DEM微分方程，使得耦合的微分方程解耦為常微分方程，經(jīng)一系列數(shù)學(xué)運(yùn)算得到一種新的可變干巖石骨架等效模型；進(jìn)一步分析可變參數(shù)的干巖石骨架等效模量的影響因素，并用實驗室測量數(shù)據(jù)對可變干巖石骨架等效模型及其它干巖石等效模型進(jìn)行對比測試。

1 方法原理

1.1 可變干巖石骨架模型

微分等效介質(zhì)(DEM)理論可以用來計算干巖石或飽和巖石的等效彈性模量：

式中：Ki,μi為包含物相的體積模量和剪切模量；K*,μ*為加入包含物相后的體積模量和剪切模量；y為包含物含量；P,Q是給定的幾何因數(shù)(Mavko等[10])；上標(biāo)(*i)指的是此幾何因數(shù)是針對具有等效模量K*,μ*的背景介質(zhì)中的包含物材料i。

對于干巖石骨架，包含物的模量Ki=μi=0，則DEM方程變?yōu)?/p>

DEM模型之所以是耦合的是因為幾何因數(shù)P,Q中含有要求解的模量信息。我們引入含有可變參數(shù)m,n的干巖石骨架等效模量：

(5)

(6)

式中：Km,μm為基質(zhì)礦物體積模量和剪切模量；Kd,μd為干巖石骨架體積模量和剪切模量；m≥0,n≥0為可變參數(shù)。

將含有可變參數(shù)的干巖石骨架等效模量代替要求的模量代入文獻(xiàn)[10]給出的P,Q表達(dá)式，就得到含有可變參數(shù)的P(m,n),Q(m,n)。這樣耦合的DEM微分方程變成了常微分方程，可以對其進(jìn)行求解。通過數(shù)學(xué)運(yùn)算(見附錄A)可以得到含有可變參數(shù)的干巖石等效彈性模量表達(dá)式，即本文導(dǎo)出的可變干巖石骨架等效模型(Variable Dry Rock Frame Equivalent Model,VDEM)：

其中：P(0),Q(0)與可變參數(shù)無關(guān)；P(1),Q(1)與可變參數(shù)有關(guān)，且當(dāng)m=n時，P(1)=Q(1)=0，此時得到的干巖石等效彈性模量與Keys-Xu近似式的結(jié)果一致；φ為孔隙度。

1.2 可變參數(shù)對干巖石等效彈性模量的影響

以石英作為基質(zhì)，模量取Km=37GPa,μm=44GPa，孔隙縱橫比α=0.1，給定孔隙度范圍，分析可變參數(shù)m,n對干巖石等效彈性模量的影響。首先分析m對干巖石等效模量的影響：令n=0，m從1變化到5，分析干巖石等效體積模量和剪切模量隨孔隙度變化，如圖1所示。再令m=0，讓n從1到7變化，分析不同n時干巖石等效模量隨孔隙度的變化，如圖2所示。

由圖1和圖2可知，干巖石體積模量Kd和剪切模量μd都隨孔隙度增加而減小；固定n時，在同一孔隙度條件下干巖石體積模量隨m的增加而增大，剪切模量隨m的增加而減小，且m對Kd影響比對μd影響更明顯；固定m時，在同一孔隙度值下，干巖石體積模量隨n增大而減小，剪切模量隨n的增加而增大。

1.3 孔隙縱橫比對干巖石等效模量的影響

同樣以石英作為基質(zhì)礦物，可變參數(shù)為m=4,n=1，孔隙度變化范圍為0～40%，取不同孔隙縱橫比α=[0.01,0.1,0.5,0.9]，分析干巖石骨架等效模量隨孔隙縱橫比的變化。由圖3可知，在同一孔隙度條件下，干巖石體積模量和剪切模量都隨孔隙縱橫比的增加而增大。

圖1 干巖石等效模量在不同參數(shù)m下隨孔隙度的變化

圖2 干巖石等效模量在不同參數(shù)n下隨孔隙度的變化

圖3 干巖石等效模量隨孔隙縱橫比α的變化

2 實驗室數(shù)據(jù)對比分析

用Han等[21]在40MPa壓力下測得的砂巖數(shù)據(jù)對本文導(dǎo)出的可變干巖石骨架等效模型進(jìn)行測試。首先用Han等給出的方法計算巖石基質(zhì)模量，然后用測量得到的飽和巖石的縱波速度、橫波速度以及密度計算得到飽和巖石模量，最后根據(jù)Gassmann方程計算干巖石骨架模量。將計算結(jié)果與Nur公式，Lee公式以及本文導(dǎo)出的可變干巖石骨架等效模型得到的干巖石骨架彈性模量進(jìn)行對比，所得結(jié)果如圖4和圖5所示。表1給出了所需公式及參數(shù)。

圖4 基于文獻(xiàn)[21]在40MPa壓力下實測砂巖數(shù)據(jù)的干巖石骨架體積模量計算結(jié)果對比

圖5 基于文獻(xiàn)[21]在40MPa壓力下實測砂巖數(shù)據(jù)的干巖石骨架剪切模量計算結(jié)果對比

參數(shù)Kdμd其它參數(shù)Nur公式Kd=Km1-φφc()μd=μm1-φφc()臨界孔隙度φc=0.4Lee公式Kd=Km1-φ1+cφμd=μm1-φ1+rcφ,r=1+2c1+c固結(jié)系數(shù)c=4VDEMKd=Km(1-φ)P(0)+P(1)eφP(1)μd=μm(1-φ)Q(0)+Q(1)eφQ(1)m=4,n=1α=0.2

從圖4中可以看到，含有固結(jié)系數(shù)的Lee公式得到的干巖石骨架體積模量相對誤差主要為負(fù)值，說明得到的干巖石體積模量偏??；從干巖石骨架體積模量相對誤差上來看，可變干巖石骨架等效模型(VDEM)和Nur公式得到的結(jié)果優(yōu)于Lee公式的結(jié)果，但VDEM所得干巖石骨架體積模量的相對誤差比Nur公式得到的相對誤差更為集中且分布在0附近，說明VDEM得到的干巖石等效體積模量和Han測得的結(jié)果更為接近。

由圖5可知，Nur公式得到的干巖石骨架剪切模量相對誤差大部分為正，說明其計算結(jié)果整體偏大；VDEM得到的干巖石骨架剪切模量與Lee公式得到的結(jié)果相比，VDEM所得結(jié)果的相對誤差小且集中分布在0附近，說明VDEM得到的干巖石等效剪切模量精度更高。

3 結(jié)束語

通過引入含有可變參數(shù)的干巖石骨架等效模量將耦合的微分等效介質(zhì)方程進(jìn)行解耦，進(jìn)一步化簡得到可變干巖石骨架等效模型(VDEM)?？勺兏蓭r石骨架等效模型避免了求解耦合的微分方程，使問題得到簡化，提高了干巖石骨架等效彈性模量的計算效率。且該模型得到的干巖石骨架彈性模量不僅與礦物基質(zhì)有關(guān)，還與孔隙度以及表征孔隙形狀的孔隙縱橫比有關(guān)，與其他經(jīng)驗公式相比提高了模型的適用性。基于實驗室測量數(shù)據(jù)的試算結(jié)果表明，與其他干巖石骨架彈性模量公式的計算結(jié)果相比，可變干巖石骨架等效模型提高了干巖石骨架等效彈性模量的預(yù)測精度。

附錄A 可變干巖石骨架模型公式推導(dǎo)

(A1)

式中：

其中：a=3Km+4μm，b=3nKm+4mμm，c=3μm，d=3mμm；f，θ是包含物縱橫比的函數(shù)(具體形式請參見文獻(xiàn)[10])。

對于干巖石骨架有Ki=μi=0，DEM方程變?yōu)?/p>

用上述的P替換方程(A8)中的P(*i)(y)可得到：

(A10)

將(10)中的y從0到φ積分可以得到：

(A11)

由此可以得到含有孔隙度為φ的干巖石等效體積模量K(φ)：

(A12)

式中：P(0)與可變參數(shù)m,n無關(guān)，而P(1)與兩個可變參數(shù)有關(guān)，且當(dāng)m,n相等或都為0時，可以得到P(1)=0，此時干巖石等效體積模量為K(φ)=Km(1-φ)P(0)。這個結(jié)果與Keys等[20]給出的干巖石體積模量近似式一致。

同理可得幾何因數(shù)Q的表達(dá)式：

(A13)

其中：

B7=5(D3D10+D4D9+D5D8)

(A30)

B8=5(D4D10+D5D9)

(A31)

B9=5D5D10(A32)D1至D10分別為

(A33)

將(A33)對y從0到φ進(jìn)行積分，可以得到

(A34)

因此含有孔隙度φ的干巖石等效剪切模量為

(A35)

式中：Q(0)與可變參數(shù)m,n無關(guān)，Q(1)和兩個調(diào)節(jié)參數(shù)有關(guān)，且當(dāng)m,n相等或都為0時，可以得到Q(1)=0，此時的干巖石等效剪切模量為μ(φ)=μm(1-φ)Q(0)，這個結(jié)果也與Keys-Xu近似式給出的干巖石剪切模量近似結(jié)果一致。

參 考 文 獻(xiàn)

[1] 郭棟,印興耀,吳國忱.橫波速度計算方法與應(yīng)用[J].石油地球物理勘探,2007,42(5):535-538

Guo D,Yin X Y,Wu G C.Methods of S-wave velocity computation and its applications[J].Oil Geophysical Prospecting,2007,42(5):535-538

[2] 周中彪.基于巖石物理模型的測井約束橫波速度計算方法研究[J].物探化探計算技術(shù),2010,32(5):536-541

Zhou Z B.The study on S-wave velocity calculation method constrained by logging based on rock physics models[J].Computing Techniques for Geophysical and Geochemical Exploration,2010,32(5):536-541

[3] 劉欣欣,印興耀,張峰.一種碳酸鹽巖儲層橫波速度估算方法[J].中國石油大學(xué)學(xué)報(自然科學(xué)版),2013,37(1):42-49

Liu X X,Yin X Y,Zhang F.An approach to predict S-wave velocity of carbonate rocks[J].Journal of China University of Petroleum,2013,37(1):42-49

[4] Gassmann F.über die elastizit?t poroser medien[J].Vierteljahrsschrift der Naturforschenden Gesellschaft in Zürich Vier，1951,96(1):1-23

[5] 印興耀,劉欣欣,曹丹平.基于Biot相洽理論的致密砂巖彈性參數(shù)計算方法[J].石油物探,2013,52(5):445-451

Yin X Y,Liu X X,Cao D P.Tight sandstone elastic parameters calculation method based on Boit-Consistent theory[J].Geophysical Prospecting for Petroleum,2013,52(5):445-451

[6] Murphy W F.Acoustic measures of partial gas saturation in tight sandstones[J].Journal of Geophysical Research,1984,89(B13):11549-11559

[7] Krief N,Garat J,Stellingwerff J,et al.A petrophysical interpretation using the velocities of P and S waves (full-waveform sonic)[J].The Log Analysis,1990,31(6):355-369

[8] Raymer L L,Hunt E R,Gardner J S.An improved sonic transit time-to-porosity transform[J].SPWLA 21stAnnual Logging Symposium,1980:1-12

[9] Nur A.Critical porosity and the seismic velocities in rocks[J].Eos,Transactions American Geophysical Union,1992,73(1):43-66

[10] Mavko G,Mukerji J,Dvorkin J.The rock physics handbook:tools for seismic analysis in porous media[M].NewYork:Cambridge University Press,2003:1-329

[11] Eshelby J D.The determination of the elastic field of an ellipsoidal inclusion,and related problems[J].Proceedings of the Royal Society of London,1957,241(1226):376-396

[12] Walsh J B.The effective of cracks on the compressibility of rock[J].Journal of Geophysical Research,1965,70(2):381-389

[13] Kuster G T,Toksoz M N.Velocity and attenuation of seismic waves in two-phase media:Part I,theoretical formulations[J].Geophysics,1974,39(5):587-606

[14] Pride S R,Berryman J G,Harris J M.Seismic attenuation due to wave-induced flow[J].Journal of Geophysical Research,2004,109:B01201

[15] Lee M W.A simple method of predicting S-wave velocity[J].Geophysics,2006,69(5):161-164

[16] 張佳佳,李宏兵,劉懷山,等.幾種巖石骨架模型的適用性研究[J].地球物理學(xué)進(jìn)展,2010,25(5):1697-1702

Zhang J J,Li H B,Liu H S,et al.Accuracy of dry frame models in the study of rock physics[J].Process of Geophysics,2010,25(5):1697-1702

[17] Cleary M P,Lee S M,Chen I W.Self-consistent techniques for the heterogeneous media[J].Journal of the Engineering Mechanics Division,1980,106(5):861-887

[18] Norris A.A differential scheme for the effective moduli of composites[J].Mechanics of Materials,1985,4(1):1-16

[19] Zimmerman R W.Compressibility of sandstones[M].New York:Elsevier,1991:1-173

[20] Keys R G,Xu S.An approximation for the Xu-White velocity model[J].Geophysics,2002,67(5):1406-1414

[21] Han D H,Nur A,Morgan D.Effects of porosity and clay content on wave velocities in sandstones[J].Geophysics,1986,51(11):2093-2107