張佳佳， 曾慶才， 印興耀， 李紅兵， 姜仁，黃家強， 張連群， 張廣智， 谷一鵬

1 中國石油大學(華東)地球科學與技術學院， 山東青島 266580 2 深層油氣地質(zhì)與勘探教育部重點實驗室， 山東青島 266580 3 中國石油勘探開發(fā)研究院， 北京 100083

0 引言

隨著油氣勘探開發(fā)程度的不斷提高，碳酸鹽巖、致密砂巖和頁巖等儲層已經(jīng)成為油氣增儲上產(chǎn)的重點領域.這些儲層孔隙類型多樣，孔隙結構復雜，非均質(zhì)性很強，屬于典型的多重孔隙儲層(含有多種孔隙類型，如孔隙、裂縫和溶洞等).多重孔隙儲層預測的難點就是尋找高可采儲量和高經(jīng)濟價值的優(yōu)質(zhì)儲層，而巖石孔隙結構是影響多重孔隙儲層油氣儲量和產(chǎn)能的重要因素(Anselmetti and Eberli，1993；Lucia，1995；Baechle et al., 2005；Li and Chen，2013).因此，巖石孔隙結構表征是多重孔隙儲層預測和流體識別的關鍵.

實驗室測量或者測井數(shù)據(jù)都有豐富的數(shù)據(jù)或直觀的手段來表征儲層孔隙結構，例如巖石薄片分析、成像測井等，但是實驗室測量或者測井數(shù)據(jù)很難進行橫向預測.而善于進行橫向預測的地震數(shù)據(jù)只有縱波速度、橫波速度以及密度三種彈性參數(shù)，無法直觀的描述儲層孔隙結構.更重要的是，多重孔隙儲層的彈性參數(shù)與孔隙結構之間的關系非常復雜，例如相同孔隙度和流體成分、不同孔隙結構的碳酸鹽巖之間的速度差異可以達到2500 m·s-1(L?n?y，2006；Sayers，2008；Baechle et al.，2005，2008，2009).利用地震數(shù)據(jù)表征儲層孔隙結構存在兩個核心問題：一個是如何描述儲層孔隙結構，另一個是如何建立孔隙結構與彈性參數(shù)之間的關系.很多學者在這兩方面做了嘗試和探索.

巖石物理中多利用孔隙縱橫比(等效橢圓的短軸與長軸之比)描述儲層孔隙形狀，因此最常見的方法是利用孔隙縱橫比表征儲層孔隙結構，再利用巖石物理模型建立彈性參數(shù)與孔隙縱橫比之間的關系.其中一類方法是利用總的等效孔隙縱橫比定性的描述儲層孔隙結構，例如總的等效孔隙縱橫比偏大，認為鑄?？椎扔部紫遁^多，反之偏小則認為裂縫等軟孔隙較多.例如Sun等(2012)根據(jù)實際測量的縱波速度反演碳酸鹽巖儲層的平均孔隙縱橫比，再對平均孔隙縱橫比進行剖分表征儲層孔隙結構.李宏兵等(2013)提出了利用Gassmann方程(Gassmann，1951)和DEM解析模型(Li and Zhang，2010，2011，2012)由縱波速度反演等效孔隙縱橫比預測橫波速度以及描述儲層孔隙結構.

另一類方法是利用不同的孔隙縱橫比來表征不同的孔隙類型(如裂縫的孔隙縱橫比值偏小，為0.01左右，基質(zhì)孔隙的孔隙縱橫比值居中，往往在0.10左右，鑄?？椎目紫犊v橫比偏大，往往在0.60左右)，再由彈性參數(shù)等反演不同孔隙類型的體積含量，相比于總的等效縱橫比方法來說更加定量化表征儲層孔隙結構.例如Kumar和Han(2005)提出了利用DEM模型(Berryman，1980)和Gassmann方程估算巖石不同孔隙類型的孔隙縱橫比以及它們各自體積含量的方法.Xu和Payne(2009)將適用于砂泥巖的Xu和White模型(1995)推廣到碳酸鹽巖，提出了一種適用于碳酸鹽巖的Xu-Payne模型，然后建立巖石縱波速度與巖石不同孔隙類型的孔隙縱橫比以及體積含量之間的關系圖版.Zhao等(2013)參照Kumar和Han(2005)以及Xu和Payne(2009)的方法，提出了利用地震數(shù)據(jù)進行巖石物理反演表征巖石的孔隙結構的方法，以此推測碳酸鹽巖儲層的地質(zhì)演化歷史.

除了利用孔隙縱橫比描述孔隙結構之外，很多學者還提出了構建一個新參數(shù)表征儲層孔隙結構，實際上這些構建的新的孔隙結構參數(shù)往往都是與孔隙縱橫比有直接或者間接的關系.例如蔣煉等(2011)運用Eshelby(1957)橢球包體裂縫理論和Gassmann方程，推導了碳酸鹽巖孔隙度與巖石孔隙縱橫比以及巖石彈性參數(shù)之間的關系式，用來模擬非均質(zhì)性較強的生物礁灘儲層的孔隙結構，該關系式包含了孔隙縱橫比的影響.Sun(2004)提出了一種孔隙結構參數(shù)用來描述巖石孔隙結構和孔隙大小對其彈性參數(shù)的影響，該孔隙結構參數(shù)與孔隙縱橫比也存在間接的關系.

在巖石物理模型除了傳統(tǒng)的經(jīng)驗模型、理論模型之外，還有一類綜合經(jīng)驗模型和理論模型的啟發(fā)式模型(Heuristic model)在行業(yè)中應用廣泛(Zhang et al.，2020)，因為啟發(fā)式模型綜合了二者的優(yōu)點.臨界孔隙度模型(Nur，1992)是一種常見的經(jīng)驗模型，因為其具有一定的物理意義和地質(zhì)含義，并且簡單實用，所以在巖石物理建模過程經(jīng)常使用.很多專家學者也對該模型進行了發(fā)展，Markov等(2012，2013)提出的GDEM(Generalized Differential Effective-Medium)模型，在傳統(tǒng)的DEM模型的基礎上，考慮了巖石的臨界孔隙度對巖石彈性特征的影響.臨界孔隙度可以用來建立巖石骨架與巖石基質(zhì)之間的關系，那么能否利用巖石的臨界孔隙度來描述儲層孔隙結構，是需要解答的問題.

針對這一問題，本文利用臨界孔隙度模型和等效介質(zhì)理論的Kuster-Toks?z方程(Kuster and Toks?z，1974)推導了臨界孔隙度與孔隙縱橫比之間的關系，進而利用極化(形狀)因子建立臨界孔隙度與彈性參數(shù)之間的關系，構建了能夠包含多種孔隙類型的臨界孔隙度模型——多孔可變臨界孔隙度模型.基于多孔可變臨界孔隙度模型，本文提出了一種利用儲層的彈性參數(shù)反演儲層不同孔隙類型的體積含量的方法.

1 多孔可變臨界孔隙度模型

1.1 臨界孔隙度與孔隙縱橫比之間的關系

等效介質(zhì)理論或者經(jīng)驗模型通常用來描述復雜儲層的孔隙結構與彈性參數(shù)之間的關系.Kuster和Toks?z(1974)基于長波長一階散射理論，使得彈性波通過等效介質(zhì)產(chǎn)生位移場和波經(jīng)每個包含物(孔隙)散射引起的位移場相同，從而推導出包含N種孔隙類型的巖石的等效彈性模量的表達式.

(1)

(2)

Nur(1992)提出了臨界孔隙度的概念，所謂臨界孔隙度就是當巖石的孔隙度超過一定臨界值時，巖石的組成礦物之間相互分散就不再是巖石了，這時候?qū)目紫抖染褪桥R界孔隙度，如圖1所示，利用臨界孔隙度建立了巖石骨架和巖石基質(zhì)彈性模量之間的關系：

(3)

(4)

式中，Kdry和μdry分別是巖石骨架的體積模量和剪切模量，φ為巖石的孔隙度，φc為巖石的臨界孔隙度.

圖1 巖石臨界孔隙度示意圖(據(jù)Nur，1992)Fig.1 Schematic diagram of critical porosity of rock (from Nur, 1992)

Nur(1992)認為同一巖性巖石的臨界孔隙度都是相同的，譬如說砂巖取40%，灰?guī)r60%等，巖石的臨界孔隙度大小取決于巖石內(nèi)部結構，它對于顆粒巖石可能居中，對于裂縫巖石可能非常小，而對于泡沫狀的巖石可能非常能夠大.

巖石孔隙結構通常用孔隙縱橫比來表征，即模擬孔隙的等效橢圓的短軸與長軸之比，如圖2所示.

結合臨界孔隙度模型和Kuster-Toks?z方程可以推導出巖石的臨界孔隙度與巖石孔隙結構(孔隙縱橫比)之間的關系，具體推導過程詳見附錄B.

(5)

(6)

式中，φcK為體積模量的臨界孔隙度值，φcμ為剪切模量的臨界孔隙度值，由式(5)和(6)可以看到巖石的臨界孔隙度值是巖石孔隙結構(孔隙縱橫比)以及巖石基質(zhì)體積模量和剪切模量的函數(shù).

假設純凈砂巖(巖石基質(zhì)體積模量Km=38 GPa和剪切模量μm=44 GPa)含有單一孔隙類型的孔隙，孔隙縱橫比分別為α=0.12，孔隙度為φ=0.30，那么巖石的臨界孔隙度可以利用式(5)和(6)計算可得到：φcK=0.36，φcμ=0.45，如圖3所示.

假設純凈碳酸鹽巖(巖石基質(zhì)體積模量Km=77 GPa和剪切模量μm=32 GPa)含有三種不同類型的孔隙，孔隙縱橫比分別為α1=0.01，α2=0.12和α3=0.8，對應為裂縫、基質(zhì)孔隙、鑄模孔，各自的孔隙體積含量分別為x1=0.04，x2=0.2和x3=0.06，那么巖石的臨界孔隙度可以利用式(5)和(6)計算可得到：φcK1=0.03，φcK2=0.36和φcK3=0.99，φcμ1=0.05，φcμ2=0.45和φcμ3=0.99，如圖4所示.

圖4 三種孔隙類型碳酸鹽巖的臨界孔隙度值Fig.4 Critical porosity values for carbonate with three pore types

1.2 多孔可變臨界孔隙度模型

建立巖石臨界孔隙度與孔隙結構之間的關系，再將其代入到Kuster-Toks?z理論模型中就可以獲得包含多種孔隙類型的臨界孔隙度模型表達式.如果考慮包含多種孔隙類型的巖石，各種孔隙類型的體積含量為xi，φcKi分別為各種孔隙的體積含量的臨界孔隙度值，φcμi分別為各種孔隙的剪切含量的臨界孔隙度值，那么就可以得到包含多種孔隙類型的臨界孔隙度模型的體積模量和剪切模量表達式(具體推導過程見附錄C)：

(7)

(8)

其中求和是對多種孔隙類型，N為巖石包含的孔隙類型的數(shù)目.

式(7)和(8)可以稱為多孔可變臨界孔隙度模型，新模型可以包含多種孔隙類型，孔隙類型可以用該孔隙類型對應的臨界孔隙度表征.不同孔隙類型的臨界孔隙度值是不同的，取決于孔隙縱橫比以及巖石基質(zhì)體積模量與剪切模量之比.值得注意的是，由于要滿足臨界孔隙度的物理意義，即巖石的孔隙度要小于臨界孔隙度，所以對于多重孔隙儲層來說，每種孔隙類型的體積含量要小于其對應的臨界孔隙度，即xi<φci，因此計算過程中需要進行迭代求解，具體迭代過程同與Xu和White模型(Xu and White，1995)一樣.

選取Han(1986)75個干燥砂巖中的10個純凈砂巖，利用單孔可變臨界孔隙度模型(即式(7)和(8)中的孔隙類型N=1)進行巖石物理建模，其中巖石基質(zhì)體積模量Km=38 GPa和剪切模量μm=44 GPa，見圖5(a和b).從圖中可以看，單孔可變臨界孔隙度模型只能用單一的孔隙縱橫比表征儲層孔隙結構，對于儲層孔隙結構變化則需要選取不同的孔隙縱橫比值，才能夠模擬得到較好的巖石彈性模量，并且只是定性描述孔隙結構，無法定量描述孔隙類型的體積含量等參數(shù).

選取Anselmetti和Eberli(1993)77個碳酸鹽巖樣品，利用雙孔可變臨界孔隙度模型進行巖石物理建模，即式(7)和(8)中的孔隙類型N=2.這些碳酸鹽巖巖石樣品的孔隙類型包括致密巖石、粒內(nèi)孔、粒間孔、鑄模孔(>300 μm)、鑄?？?<300 μm)和微孔隙等六種，如圖5所示.這里假設有三種孔隙類型在這些碳酸鹽巖樣品中，硬孔隙、基質(zhì)孔隙和軟孔隙，它們的孔隙縱橫比分別為0.6、0.15和0.03，巖石基質(zhì)體積模量Km=77 GPa和剪切模量μm=32 GPa.這里雖然劃分為三種孔隙類型，但是仍然要利用雙孔可變臨界孔隙度模型，主要原因是目前實際巖石物理模擬過程大多數(shù)都假設巖石中包含兩種孔隙類型，即同時包含硬孔隙和基質(zhì)孔隙，或者同時包含基質(zhì)孔隙和軟孔隙(Kumar and Han，2005； Xu and Payne，2009).孔隙類型確定好以后，就需要模擬不同孔隙類型的體積含量的變化，圖6顯示了77個碳酸鹽巖樣品的縱波速度隨孔隙度以及硬孔體積含量或軟孔體積含量變化曲線，從圖中可以看到碳酸鹽巖的縱波速度隨著孔隙度增加以及軟孔含量增加或者硬孔含量減小而減小.這就表明利用雙孔可變臨界孔隙度模型同時包含兩種孔隙類型(軟孔隙與基質(zhì)孔隙組合或者硬孔隙與基質(zhì)孔隙組合)，改變不同孔隙類型的體積含量，就能夠模擬得到準確的碳酸鹽巖速度.相對于單孔可變臨界孔隙度模型，多孔可變臨界孔隙度模型可以定量描述儲層孔隙類型的體積含量等.

圖5 利用單孔可變臨界孔隙度模型模擬砂巖彈性模量(a) 體積模量; (b) 剪切模量.Fig.5 Simulating sandstone elastic modulus with single-porosity variable critical porosity model(a) Bulk modulus; (b) Shear modulus.

圖6 利用雙孔可變臨界孔隙度模型模擬碳酸鹽巖縱波速度Fig.6 Simulating P-wave velocity of carbonate rocks with the dual-porosity variable critical porosity model

2 基于多孔可變臨界孔隙度模型的儲層孔隙結構表征

多孔可變臨界孔隙度模型可以包含多種孔隙類型，如孔隙、裂縫和溶洞等，不同的孔隙類型的孔隙結構不相同，臨界孔隙度值也不相同.因此，可以利用多孔可變臨界孔隙度模型進行孔隙結構表征.當儲層巖石含有多種孔隙類型的時候，如碳酸鹽巖，可以仿照Kumar和Han的思想，利用多孔可變臨界孔隙度模型反演不同孔隙類型的體積含量.具體實現(xiàn)步驟為：

(1)選取巖石基質(zhì)礦物組分的彈性模量和孔隙流體組成成分的彈性模量和密度，計算巖石基質(zhì)的體積模量Km、剪切模量μm和縱波速度VPm和孔隙流體的體積模量Kfl和縱波速度VPfl：

(9)

(10)

(11)

(12)

(13)

ρfl=Soilρoil+Sgasρgas+Swaterρwater,

(14)

(15)

(16)

(17)

(18)

(19)

(20)

(21)

(22)

(23)

(7)采用非線性全局尋優(yōu)算法修改步驟(3)給定的不同孔隙類型的體積含量，計算新的碳酸鹽巖儲層巖石骨架的彈性模量，重復步驟(3)—(6)，計算相應的誤差，比較誤差之間的大小，最后得到最優(yōu)的不同孔隙類型的體積含量.

圖7為利用多孔可變臨界孔隙度模型對某碳酸鹽巖儲層不同孔隙類型體積含量反演的結果.這里假設有三種孔隙類型，硬孔、基質(zhì)孔隙和軟孔，它們的孔隙縱橫比α分別為0.8、0.12和0.01.利用多孔可變臨界孔隙度模型反演得到三種不同孔隙類型的體積含量，圖7f為硬孔孔隙度，圖7g為軟孔孔隙度.

圖8為利用多孔可變臨界孔隙度模型對某致密砂巖儲層軟孔隙體積含量反演的結果.這里假設有兩種孔隙類型，基質(zhì)孔隙和軟孔，它們的孔隙縱橫比α分別為0.12和0.01.利用多孔可變臨界孔隙度模型反演得到兩種不同孔隙類型的體積含量，圖8g為軟孔隙的體積含量，可以看到在3775～3880 m軟孔隙的體積含量占總的孔隙度的5%以上，指示了軟孔隙含量較多，實際鉆井表明該層段裂縫含量較多，滲透性很好，是非常好的“甜點”儲層.

3 討論

需要指出的是，本文的提出的多孔可變臨界孔隙度模型對傳統(tǒng)的臨界孔隙度模型的發(fā)展和修正，建立了巖石臨界孔隙度與儲層孔隙結構(孔隙縱橫比)之間的理論關系，是一種綜合經(jīng)驗模型和理論模型的啟發(fā)式模型(Heuristic model).多孔可變臨界孔隙度模型是基于等效介質(zhì)理論Kuster-Toks?z方程推導出來的，雖然Kuster-Toks?z方程要求孔隙彼此之間是孤立的，近似高頻假設條件，但是多孔可變臨界孔隙度模型計算過程中采用了類似Xu-White模型的計算過程，每一步計算得到的等效介質(zhì)又被作為新的背景介質(zhì)，因此這種迭代計算方法既考慮了巖石孔隙之間的相互作用，同時也可以一定程度上解決巖石物理中不同的尺度問題(頻率問題).另外，針對巖石不同頻率(超聲波高頻與地震低頻)的彈性特征差異還可以利用臨界孔隙度值差異來進行等效模擬.

圖7 利用多孔可變臨界孔隙度模型反演碳酸鹽巖儲層硬孔孔隙度和軟孔孔隙度(a) 密度； (b) 泥質(zhì)含量； (c) 孔隙度； (d) 含水飽和度； (e) 縱波速度； (f) 反演得到的硬孔孔隙度； (g) 反演得到的軟孔孔隙度.Fig.7 Inversion of porosity of different pore types of carbonate reservoir based on multiple-porosity variable critical porosity model(a) Density； (b) Shale volume content； (c) Porosity； (d) Water saturation ；(e) P-wave velocity； (f) Inverted porosity of stiff pore； (g) Inverted porosity of soft pore.

圖8 利用多孔可變臨界孔隙度模型反演致密砂巖軟孔隙的體積含量(a) 密度; (b) 泥質(zhì)含量; (c) 孔隙度; (d) 含水飽和度; (e) 縱波速度; (f) 橫波速度; (g) 反演得到的軟孔體積含量.Fig.8 Inversion of volume content of soft pore of tight sandstone based on multiple-porosity variable critical porosity model(a) Density; (b) Shale volume content; (c) Porosity; (d) Water saturation; (e) P-wave velocity; (f) S-wave velocity; (g) Inverted volume content of soft pore.

4 結論

本文結合經(jīng)驗模型臨界孔隙度模型和等效介質(zhì)理論Kuster-Toks?z方程，推導了巖石的臨界孔隙度與巖石孔隙結構的關系，并且推導了能夠包含多種孔隙類型的臨界孔隙度模型，并將之定義為多孔可變臨界孔隙度模型.基于多孔可變臨界孔隙度模型可以利用儲層的彈性參數(shù)反演不同孔隙類型的體積含量.實驗室測量數(shù)據(jù)和實際測井數(shù)據(jù)表明，多孔可變臨界孔隙度模型能夠適用于多重孔隙儲層巖石物理建模和孔隙結構表征.

致謝感謝匿名評閱人對本文提出的寶貴修改意見.

附錄A 極化因子P和Q表達式

具有任意孔隙縱橫比的橢球狀包含物的P和Q系數(shù)可寫為(Kuster和Toks?z，1974)：

(A1)

(A2)

其中的張量Tijkl將均勻遠場應變場與橢球包含物的應變相聯(lián)系(Wu，1966)，Beryman(1980)給出了計算P和Q所需要的相關的標量：

Tiijj=3F1/F2,

(A3)

(A4)

其中

+B(3-4R)+(A/2)(A+3B)(3-4R)

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

F4=1+(A/4)[f+3θ-R(f-θ)],

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

F7=2+(A/4)[3f+9θ-R(3f+5θ)]

+Bθ(3-4R)

F8=A[1-2R+(f/2)(R-1)+(θ/2)(5R-3)]

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

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

其中A，B和R為

A=μi/μm-1

且

R=(1-2νm)/2(1-νm)

函數(shù)θ為

以上分別針對扁長和扁圓橢球體，且

注意：對扁圓橢球體α<1，對扁長橢球體α>1.

附錄B 臨界孔隙度與孔隙縱橫比之間的關系

(B1)

當巖石包含單一孔隙類型的時候，則孔隙的體積含量就等于巖石的孔隙度，即x=φ，那么式(B1)變?yōu)?/p>

(B2)

因此，當φ=φcK時，φcK定義為體積模量的臨界孔隙度值，可以得到：

(B3)

式(B3)可改寫為

(B4)

這就是體積模量的臨界孔隙度值φcK與巖石孔隙結構α之間的關系表達式.

(B5)

當巖石包含單一孔隙類型時，孔隙的體積含量就等于巖石的孔隙度，即x=φ，那么式(B5)變?yōu)?/p>

(B6)

因此，當φ=φcμ時，φcμ定義為剪切模量的臨界孔隙度值，可以得到：

(B7)

式(B7)可以改寫為

(B8)

這就是剪切模量的臨界孔隙度值φcμ與巖石孔隙結構α之間的關系表達式.

附錄C 多孔可變臨界孔隙度模型

如果考慮包含多種孔隙類型的巖石，各種孔隙類型的體積含量為xi，φcKi分別為各種孔隙的體積含量的臨界孔隙度值，那么有

(C1)

因此

(C2)

(C3)

再將式(C3)代入式(B1)中，就可以得到包含多種孔隙類型的臨界孔隙度模型的體積模量表達式

(C4)

其中求和是對多種孔隙類型，N為巖石包含的孔隙類型的數(shù)目.

同樣的，如果考慮包含多種孔隙類型的巖石，各種孔隙類型的體積含量為xi，φcμi分別為各種孔隙的剪切含量的臨界孔隙度值，那么有

(C5)

因此

(C6)

(C7)

再將式(C7)代入式(B5)中，就可以得到包含多種孔隙類型的臨界孔隙度模型的剪切模量表達式

(C8)

其中求和是對多種孔隙類型，N為巖石包含的孔隙類型的數(shù)目.