考慮成巖固結作用的雙重孔隙介質巖石物理模型及其應用

高 強 李紅兵* 潘豪杰

(①中國石油勘探開發研究院油氣地球物理研究所，北京 100083； ②長江大學地球物理與石油資源學院，湖北武漢 430110)

0 引言

縱、橫波速度在儲層物性(孔隙度、飽和度等)參數預測和巖性及孔隙流體識別等方面扮演著不可替代的角色。準確的縱波和橫波速度信息對精確描述儲層特征、消除流體識別假象以及提高疊前地震反演和AVO屬性分析精度等方面具有重要意義。然而，在實際油氣勘探中，由于偶極聲波測井測量施工成本較高，造成大多數測井數據中缺少橫波速度資料。為此，人們廣泛嘗試利用巖石物理理論預測橫波速度[1-4]。

Castagna等[5]給出了水飽和狀態下碎屑硅質巖橫波速度與縱波速度的經驗關系的泥巖基線公式。Han等[6]針對75塊砂巖樣品開展了不同壓力下的超聲實驗測量，擬合出橫波速度與縱波速度的線性經驗關系式。Nur等[7]假設巖石孔隙度小于臨界孔隙度時，干巖石模量與巖石孔隙度線性相關，構建了臨界孔隙度模型。

雖然上述經驗模型使用方便，但其物理含義欠明確且需結合大量實際工區的測井和巖心數據進行校正。Kuster等[8]利用長波一階散射理論推導出彈性模量與孔隙度和孔隙形狀之間的定量關系式。Xu等[9]將巖石孔隙分成砂巖孔隙和泥巖孔隙，引入砂泥巖孔隙縱橫比表征孔隙形狀，綜合應用Wyllie時間方程[10]、KT模型[8]和Gassmann方程[11]建立了砂泥巖巖石物理模型。該模型考慮了孔隙結構及孔隙形狀對巖石彈性模量的影響，具有明確物理意義，在橫波速度預測及巖石物理分析中得到廣泛應用[12-13]。張廣智等[14]提出了改進的Xu-White模型并應用于碳酸鹽巖橫波速度預測。

除了孔隙結構和孔隙類型外，巖石顆粒之間固結作用的強弱對巖石縱、橫波速度也有重要影響。Pride等[15]通過引入固結參數表征巖石顆粒之間固結程度，建立了干巖石模量與基質模量的關系式。Lee[16]修正了Pride模型計算剪切模量的表達式，并將其應用于橫波速度預測。張佳佳等[17]和侯波等[18]都將Pride模型成功地應用于橫波速度預測。

上述經驗公式和理論模型雖可較好地預測巖石縱、橫波速度，但大多數模型僅考慮了成巖固結作用或聚焦于單一孔隙結構，缺少能同時兼顧多重孔隙結構與成巖作用的有效理論與方法，導致現有模型在復雜孔隙介質儲層預測中呈現出適用性較差、穩定性不足和精度較低等特征。

為了兼顧成巖固結作用與雙重孔隙結構對巖石骨架彈性模量的影響，本文提出一種考慮成巖固結作用的雙重孔隙介質巖石物理模型。首先，基于微分等效介質(DEM)理論[19]推導出含有(軟、硬)雙重孔隙DEM解析模型； 再從雙重孔隙DEM模型出發，借鑒Pride模型思想，建立了考慮成巖固結作用的雙重孔隙介質巖石物理模型； 然后，結合Gass-mann方程建立飽和巖石縱、橫波速度計算模型； 最后，將該模型分別應用于實驗數據和測井數據的速度預測，并對其預測結果進行對比分析，驗證了本文所提考慮成巖固結作用的雙重孔隙介質巖石物理模型的適用性和有效性。

1 方法原理

1.1 單孔隙DEM模型及其近似解

經典DEM理論是通過將等體積增量的填充物逐漸添加到背景固體礦物中模擬雙相等效介質，直至添加的填充物體積為真實孔隙度[20-21]。Berryman[22]給出DEM模型的微分方程組耦合形式

(1)

(2)

式中：Ki和μi分別為第i種填充物的體積模量和剪切模量；K*和μ*對應為加入填充物后的體積模量和剪切模量；y為填充物含量；P*和Q*為填充物孔隙形狀幾何因子。

對于干巖石骨架，填充物的體積模量和剪切模量分別為Ki=0和μi=0，則式(1)、式(2)可簡化為

(3)

(4)

對式(3)、式(4)從y=0到y=φ進行積分，并滿足初始條件K*(0)=Kma、μ*(0)=μma，令Kdry=K*、μdry=μ*，得到單孔隙結構干巖石骨架模量

Kdry=Kma(1-φ)P

(5)

μdry=μma(1-φ)Q

(6)

式中：Kdry、μdry分別為干巖石骨架的體積模量和剪切模量；Kma、μma對應為巖石基質的體積模量和剪切模量；P和Q為與孔隙形狀有關的幾何因子；φ為巖石的總孔隙度。

1.2 雙重孔隙DEM模型

碎屑巖儲層中通常發育兩種孔隙類型，即微裂縫和孔隙[23]。本文利用DEM理論模擬雙孔巖石的彈性性質，即是將單孔隙DEM理論拓展到雙重孔隙DEM模型。為此，將上述孔隙對應定義為軟孔和硬孔兩種類型，分別用αso、αst表示軟孔和硬孔的孔隙縱橫比?？紤]孔隙添加順序，首先向巖石基質中加入軟孔，按式(5)和式(6)可得巖石骨架模量

(7)

(8)

然后，將已引入軟孔的基質視作新巖石基質，并向其中加入硬孔，可得最終干巖石骨架的模量

(9)

(10)

式中：φst為硬孔的孔隙度，且有φ=φso+φst；Pst和Qst為硬孔對應的幾何因子，是關于孔隙縱橫比αst的函數，計算方法同上。將式(7)、式(8)代入式(9)和式(10)，可得含有雙重孔隙的干巖石骨架的模量

Kdry=Kma(1-φso)Pso(1-φst)Pst

(11)

μdry=μma(1-φso)Qso(1-φst)Qst

(12)

若先向巖石基質中加入硬孔，按照式(5)和式(6)可得到此時干巖石骨架的模量

(13)

(14)

然后，將已引入硬孔的基質視作新的巖石基質，并向其中加入軟孔，可得最終干巖石骨架的模量

(15)

(16)

將式(13)、式(14)代入式(15)和式(16)中，可得含雙重孔隙的干巖石骨架模量

Kdry=Kma(1-φst)Pst(1-φso)Pso

(17)

μdry=μma(1-φst)Qst(1-φso)Qso

(18)

與式(11)、式(12)對比可知，由雙重孔隙DEM模型計算的彈性模量與兩種孔隙的加入順序無關。

1.3 考慮成巖固結作用的雙重孔隙模型

雙重孔隙DEM模型可直接求解、計算簡便且考慮了孔隙結構特征對彈性波速度的影響，但該模型忽略了巖石顆粒之間的固結壓實作用對巖石彈性模量的影響。Pride等[15]通過引入固結參數表征巖石顆粒之間的固結程度，提出固結巖石的體積模量和剪切模量的簡單表達式

(19)

(20)

式中：WK(c)、Wμ(c)為反映巖石骨架固結程度的函數，且有WK(c)=1+cφ、Wμ(c)=1+1.5cφ；c為固結參數，對于砂巖，通常取2

Pride模型是一個經驗模型，它揭示了巖石顆粒之間的固結壓實作用對巖石彈性模量的影響。當c接近于2時，表明巖石顆粒之間固結程度較大，巖石骨架表現為更堅硬，彈性模量較高； 當c接近于20時，表明巖石顆粒之間固結程度較小，巖石骨架表現為更疏松，彈性模量較低。

侯波等[24]將Pride模型思路應用于Keys-Xu模型，提出綜合成巖作用和孔隙形狀的巖石物理模型

(21)

(22)

式中

由式(19)、式(20)可知，Pride模型對巖石彈性模量的影響主要表現在分母上增加了一個固結參數的函數。本文采用相同做法，將Pride模型思路應用于新導出的雙重孔隙DEM模型，提出考慮成巖固結作用的雙重孔隙巖石物理模型

(23)

(24)

代入式(23)和式(24)，可得

(25)

(26)

當c逐漸減少并接近于0，即巖石完全固結時，式(25)和式(26)就變成式(11)和式(12)，即DEM模型的近似模型。式中φ=φso+φst=φ(fso+fst)，其中fso、fst分別為軟孔、硬孔孔隙所占總孔隙的比例。

為了探究孔隙度、固結參數、孔隙縱橫比和軟孔體積比對干巖石骨架模量和速度的影響，首先對本文模型、Xu-White模型和Pride模型進行數值模擬并對比其結果。以石英作為主介質，其體積模量和剪切模量分別取為38、44GPa。巖石中含有兩種不同形狀的孔隙，其孔隙縱橫比分別為0.12和0.03，其孔隙體積占比分別取為98%和2%。

圖1、圖2分別是干巖石骨架模量、速度隨孔隙度的變化特征。其中藍色、綠色、黑色實線代表本文模型(PDEM)在固結參數分別取c=8，4，0時計算結果； 紅色虛線代表Xu-White模型計算結果； 藍色、綠色、黑色虛線代表Pride模型在固結參數分別取c=8，4，0時計算結果。從這兩圖可看出： 干巖石骨架彈性模量和速度隨總孔隙度的增加而減??； 隨固結參數增加(固結程度降低)而降低。當固結參數c為0時，本文模型計算結果與Xu-White模型計算結果相同； 當c>0時，本文模型預測結果均低于Xu-White模型。這主要是因為在DEM理論推導中，Xu-White模型假設巖石骨架是由一個整體構成，而實際上巖石骨架是由巖石顆粒膠結而成的。膠結程度會影響巖石骨架模量和速度。而Pride模型計算結果均大于另兩種模型，主要因為計算干巖石骨架模量時，Pride模型輸入參數僅有孔隙度和固結參數，卻未考慮孔隙結構的影響。

圖1 干巖石體積模量(a)和剪切模量(b)隨孔隙度的變化特征

圖2 縱(a)、橫(b)波速度隨孔隙度的變化特征

圖3和圖4展示孔隙度為0.2、固結參數為2時，干燥巖石彈性模量隨孔隙縱橫比和軟孔體積比的變化特征。從圖3可看出，干巖石骨架彈性模量隨孔隙縱橫比的增加而增加，當軟孔縱橫比增加到0.008時，其增長率變緩。從圖4可見，隨著軟孔體積比的增加，干巖石骨架模量總體呈下降趨勢。

圖3 干巖石體積模量(a)和剪切模量(b)隨軟孔隙縱橫比的變化特征

圖4 干巖石體積模量(a)和剪切模量(b)隨軟孔體積比的變化特征

1.4 飽和巖石縱橫波速度計算

由式(25)、式(26)求得干巖石骨架的體積模量和剪切模量后，可通過Gassmann方程計算流體飽和巖石的體積模量和剪切模量

(27)

μsat=μdry

(28)

最后，由縱、橫波速度與流體飽和巖石的彈性模量和體密度之間的關系，可得到縱、橫波速度

(29)

(30)

式中：VP、VS分別為縱波和橫波速度；ρsat為飽和巖石的體密度，其計算公式為

ρsat=ρma(1-φ)+[ρwSw+ρg(1-Sw)]φ

(31)

式中ρma、ρw、ρg分別為巖石基質、地層水、天然氣的密度。

2 橫波速度預測

2.1 技術流程

利用本文模型預測橫波速度，需提供兩個關鍵參數：固結參數和軟孔體積比。這兩參數可利用本文優化計算流程(圖5)求得。其主要步驟為：

圖5 固結參數和軟孔體積比反演流程圖

(1)利用式(25)和式(26)推導的巖石物理模型計算干巖石骨架模量，通過Gassmann方程進行流體替代得到飽和巖石彈性模量，進而計算縱波速度；

(2)基于模擬退火算法不斷調整輸入參數，優化固結參數c和軟孔隙體積比fso，使預測縱波速度與測量縱波速度的誤差達到最小，進而得到優化后的固結參數和軟孔體積比。

得到優化的固結參數和軟孔體積比后，即可用“1.4節”飽和流體巖石物理模型(式(30))計算橫波速度。

2.2 實驗數據

為了檢驗方法的適用性，將本文模型分別應用于壓力為5和40MPa條件下的實驗測量數據[6]和蘇里格地區壓力為29MPa條件下的實驗測量數據[26]。通過反演流程(圖5)，確定壓力為5MPa時固結參數c為2.9，軟孔體積占比fso為6%； 壓力為40MPa時固結參數c為2.0，軟孔體積占比fso為1.5%； 壓力為29MPa時固結參數c為2.0，軟孔體積占比fso為2.5%。從圖6和圖7可看出，在不同壓力條件下，模型預測結果與實驗數據絕大多數均吻合較好，均在擬合直線附近。需指出的是，由于蘇里格實驗數據中缺少泥質含量數據，本文對所有巖石樣品采用統一的體積模量和剪切模量，這會對橫波速度的預測精度造成影響。模型預測結果與實驗數據對比表明，該模型對常規碎屑巖和致密砂巖儲層均有較好的適應性。

圖6 壓力為5(a)、40(b)MPa時預測的與實測的橫波速度[6]

圖7 壓力為29MPa時預測的與實測的的橫波速度

2.3 測井數據

圖8顯示蘇里格地區M井測井數據，包括縱波速度、橫波速度、密度及解釋的孔隙度、含水飽和度和泥質含量。該工區目的層處于2900～3007m，為低孔、低滲致密砂巖儲層，非儲層段孔隙度幾乎為零。此外，由于孔隙之間連通性差，氣、水關系復雜，使流體在致密砂巖氣藏中通常呈不均勻分布。本次建模所涉各種礦物的彈性模量和密度值參考表1，將分別采用不同的巖石物理模型預測橫波速度并做誤差對比分析。

表1 各種礦物彈性模量及密度

通過不斷調整所選用巖石物理模型的自由參數，使其預測的橫波速度與實測值之間誤差最小。如臨界孔隙度模型(CMP)[7]的臨界孔隙度為0.39，Pride模型[15]的固結參數為2，本文模型(PDEM)的固結參數和軟孔體積比分別為2和0.05，Xu-White(XW)[9]模型中砂巖孔隙縱橫比為0.09、泥巖孔隙縱橫比為0.056，KT[8]、SC[27]和DEM[19]模型的等效孔隙縱橫比為0.13。

將本文模型和其他七種模型(Krief模型[28]、臨界孔隙度模型、Pride模型、Xu-White模型、KT模型、DEM模型、SC模型)應用于圖8測井數據，分別得到圖9所示的橫波速度預測結果(紅線)，并與實測橫波速度(藍線)進行對比?？梢娫?920～2950m、2960～2974m及2980～3007m三個層段，本文模型預測結果與實際橫波測井數據的吻合程度好于其他模型。這是因為本文模型引入了可變固結參數(不同沉積環境及不同巖性的巖石的固結參數通常有差異)和兩種孔隙類型。而非儲層段的微小差異，可能是由含水飽和度和孔隙度解釋的不準確造成的。

圖8 蘇里格M井的測井數據及解釋的物性參數

圖9 壓力為29MPa時預測與實測的橫波速度(a)Krief模型； (b)臨界孔隙度模型； (c)Pride模型； (d)Xu-White模型； (e)KT模型； (f)DEM模型； (g)SC模型； (h)PDEM模型

圖10為基于八種巖石物理模型預測的橫波速度與實際橫波測井數據之間的均方根誤差直方圖?？梢娀诒疚哪Ｐ偷臋M波速度預測結果的均方根誤差(143.8m/s)遠小于其他模型預測結果的均方根誤差(平均為220.0m/s)，其中Krief模型的均方根誤差更是高達278.0m/s。因此，本文模型對致密砂巖儲層的橫波速度具有更高預測精度。

圖10 八種巖石物理模型預測的橫波速度的均方根誤差

3 討論

本文推導出雙重孔隙DEM模型以表征復雜孔隙儲層的孔隙結構，并研究儲層彈性參數與孔隙結構和孔隙度的定量關系?？紫督Y構的基本元素包括孔隙類型、孔隙形狀、孔隙連通性等，現階段可通過鑄體薄片、SEM觀測、微納米CT掃描成像和數字巖心分析技術等微觀結構成像手段做直接表征。如何應用地球物理方法間接定量表征儲層復雜孔隙結構，一些學者曾開展了探索性研究。李紅兵等[29]應用DEM模型和Gassmann方程由縱橫波速度反演孔隙結構，并將其應用于橫波速度預測。Zhang等[30]利用Sun氏模型從疊前地震數據中反演表征孔隙類型的孔隙結構因子。Li等[31]利用Xu-Payne模型從疊前地震數據定量預測了碳酸鹽巖儲層軟孔、硬孔和參考孔的孔隙度大小。這些理論和方法的成功應用為利用本文模型實現微觀孔隙結構量化表征提供了理論借鑒和技術支撐。

成巖固結作用和孔隙結構及類型是影響巖石彈性性質的兩類重要因素，在構建巖石物理模型時應兼顧考慮。經典Xu-White模型只考慮了孔隙結構，圖11是針對同樣的砂巖樣品預測的橫波速度與實測橫波速度的交會圖[8]。與本文模型預測結果(圖6)對比，可見Xu-White模型預測的橫波速度(圖11)整體偏高。這可能是因為Xu-White模型未考慮成巖固結作用的影響。圖12是基于同樣砂巖樣品采用Pride模型預測的橫波速度與實測橫波速度的交會圖。再與本文模型預測結果(圖6)對比，可見本文模型預測結果更靠近于擬合直線，即誤差更小。從測井數據橫波速度均方根誤差(圖10)也可看出，本文模型預測結果的均方根誤差均小于Xu-White模型和Pride模型。干巖石骨架彈性模量不僅與孔隙度和孔隙形狀有關，且與礦物顆粒之間固結程度關聯。本文模型同時考慮孔隙度、孔隙形狀和固結參數，可從成巖固結作用和孔隙結構兩方面解釋干巖石彈性模量的影響機制，顯著減少了僅考慮單一因素導致的計算誤差，能更全面、精確地描述巖石骨架彈性模量隨孔隙度的變化特征，具有更廣闊應用前景。

圖11 Xu-White模型預測的與實測的橫波速度交會[9]

圖12 Pride模型預測的與實測的橫波速度交會

4 結束語

本文基于DEM理論導出雙重孔隙DEM模型，并與考慮成巖固結作用的Pride模型相結合，提出一種兼顧成巖固結作用和孔隙結構影響的巖石物理模型。通過引入固結參數、軟孔和硬孔縱橫比、體積占比，實現對巖石成巖作用和微觀孔隙結構的聯合表征。不同壓力下的巖石物理實驗數據的測試結果表明，該模型對常規碎屑巖儲層和致密砂巖儲層具有較強適用性。與其他巖石物理模型預測結果相比，本文模型預測的橫波速度精度更高。

需指出的是，應用本文模型預測橫波速度過程中，需優化多個參數，因此會增加一定計算量。

附錄A 孔隙縱橫比表達式

巖石填充物的幾何因子P和Q的計算公式為

式中i、j均取1、2、3，對應x、y、z坐標。且有

其中

[f+θ-R(f-θ+2θ2)]

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

B(1-θ)(3-4R)

B(1-θ)(3-4R)

F9=A[(R-1)f-Rθ]+Bθ(3-4R)

上列各式中

A=μi/μma-1

R=(1-2σma)/[2(1-σma)] (σ為泊松比)

該式中第一、第二行分別對應扁長、扁圓球體。

對于扁圓球體，α<1； 對于扁長球體，α>1。