劉倩, 印興耀, 李超

1 中國(guó)石油大學(xué)(華東)地球科學(xué)與技術(shù)學(xué)院， 青島 266580 2 中海油研究總院，北京 100028

?

基于逆建模理論的儲(chǔ)層特征定量預(yù)測(cè)方法

劉倩1, 印興耀1, 李超2

1 中國(guó)石油大學(xué)(華東)地球科學(xué)與技術(shù)學(xué)院， 青島 266580 2 中海油研究總院，北京 100028

定量描述儲(chǔ)層特征的物理量稱為儲(chǔ)層參數(shù)，在儲(chǔ)層描述工作中儲(chǔ)層參數(shù)預(yù)測(cè)是一個(gè)重要環(huán)節(jié).本文以巖石物理理論為基礎(chǔ)研究了利用巖石物理逆建模理論進(jìn)行儲(chǔ)層參數(shù)預(yù)測(cè)的方法.在合理、有效的巖石物理模型的基礎(chǔ)上，逆建模方法通過(guò)彈性參數(shù)的等值面的空間交會(huì)可以預(yù)測(cè)出儲(chǔ)層參數(shù).在此基礎(chǔ)上提出了一種擴(kuò)展的巖石物理逆建模方法，在兩個(gè)不同的儲(chǔ)層參數(shù)域進(jìn)行三維巖石物理逆建模，通過(guò)逆建模結(jié)果的參數(shù)域轉(zhuǎn)換可以預(yù)測(cè)更多的儲(chǔ)層參數(shù)，為儲(chǔ)層特征描述提供了更加充分的數(shù)據(jù)支持.輸入數(shù)據(jù)的種類對(duì)于預(yù)測(cè)結(jié)果的準(zhǔn)確性有較大影響，利用本文建立的彈性參數(shù)適用性分析方法可以選擇出適用性好的彈性參數(shù)組合作為輸入數(shù)據(jù).將本方法應(yīng)用于模型數(shù)據(jù)和實(shí)際工區(qū)數(shù)據(jù)，取得了較好的預(yù)測(cè)效果，證明方法具有一定的實(shí)用性.

儲(chǔ)層參數(shù)； 巖石物理； 逆建模； 儲(chǔ)層參數(shù)域

1 引言

儲(chǔ)層預(yù)測(cè)的最終目標(biāo)是找到含油氣儲(chǔ)層并描述儲(chǔ)層的物理性質(zhì).一直以來(lái)，地球物理學(xué)家們采用多種不同的思路研究如何將縱波速度、橫波速度、彈性阻抗等數(shù)據(jù)轉(zhuǎn)換為描述儲(chǔ)層物理性質(zhì)的數(shù)據(jù)如孔隙度、泥質(zhì)含量、含水飽和度、孔隙縱橫比等.Doyen(1988)基于地質(zhì)統(tǒng)計(jì)學(xué)理論考慮空間約束條件從地震資料中預(yù)測(cè)孔隙度；Mukerji等(2001)和Eidsvik等(2004)綜合巖石物理理論、統(tǒng)計(jì)模式識(shí)別理論和信息論提出利用統(tǒng)計(jì)巖石物理方法來(lái)反演儲(chǔ)層物性參數(shù)；Avseth和?degaard(2004)發(fā)展了基于巖石物理模板的孔隙度和含水飽和度的估算方法；Bachrach(2006)將巖石物理隨機(jī)建模方法和貝葉斯反演理論相結(jié)合，進(jìn)行孔隙度和含水飽和度的聯(lián)合反演；Spikes等(2007)提出了一種聯(lián)合地震資料和測(cè)井?dāng)?shù)據(jù)的概率化地震反演技術(shù)實(shí)現(xiàn)儲(chǔ)層參數(shù)預(yù)測(cè)；Bosch等(2010)提出一種綜合馬爾科夫鏈蒙特卡羅方法和迭代最小二乘最優(yōu)化方法的地質(zhì)統(tǒng)計(jì)學(xué)物性參數(shù)反演方法；Grana和Rossa(2010)基于貝葉斯框架和統(tǒng)計(jì)巖石物理建模理論提出了儲(chǔ)層物性參數(shù)的反演方法和巖性流體分類方法；楊午陽(yáng)和王叢鑌(2010)利用基于模擬退火的疊前同步反演方法實(shí)現(xiàn)了儲(chǔ)層物性參數(shù)的反演；胡華鋒等(2012)綜合統(tǒng)計(jì)巖石物理理論和貝葉斯分類方法提出了物性參數(shù)的聯(lián)合反演方法；Jiang和Spikes(2012)基于自洽模型建立了一種搜索算法來(lái)預(yù)測(cè)頁(yè)巖儲(chǔ)層的孔隙度和孔隙縱橫比；李宏兵等(2013)基于Gassmann方程和DEM解析模型利用非線性全局尋優(yōu)算法來(lái)預(yù)測(cè)等效孔隙縱橫比；印興耀等(2014)建立了彈性阻抗與儲(chǔ)層物性參數(shù)之間的線性關(guān)系，利用彈性阻抗反演方法實(shí)現(xiàn)了物性參數(shù)的估計(jì)；Yin等(2014)基于FFT-MA模擬方法和GDM擾動(dòng)方法提出了地質(zhì)統(tǒng)計(jì)學(xué)反演的方法來(lái)反演物性參數(shù).

本文基于巖石物理理論建立了儲(chǔ)層參數(shù)域，在合理、有效的巖石物理模型的基礎(chǔ)上，在儲(chǔ)層參數(shù)域內(nèi)計(jì)算得到彈性參數(shù)的數(shù)據(jù)場(chǎng)和等值面，將三個(gè)不同彈性參數(shù)的單值等值面進(jìn)行空間交會(huì)，就可以預(yù)測(cè)出這組彈性參數(shù)所對(duì)應(yīng)的儲(chǔ)層參數(shù).在此基礎(chǔ)上提出了擴(kuò)展的巖石物理逆建模方法，在兩個(gè)不同的三維儲(chǔ)層參數(shù)域分別進(jìn)行多次巖石物理逆建模，然后通過(guò)參數(shù)域轉(zhuǎn)換及空間交會(huì)可以進(jìn)行四個(gè)儲(chǔ)層參數(shù)的預(yù)測(cè)，為儲(chǔ)層特征定量預(yù)測(cè)提供更充分的數(shù)據(jù)支持.基于逆建模方法進(jìn)行儲(chǔ)層參數(shù)預(yù)測(cè)的優(yōu)點(diǎn)在于對(duì)輸入數(shù)據(jù)的選取較為靈活，可以根據(jù)不同工區(qū)的特點(diǎn)選取不同的彈性參數(shù)作為輸入數(shù)據(jù)，而建立的巖石物理模型也是有針對(duì)性的，以此可以保證儲(chǔ)層參數(shù)預(yù)測(cè)的精確程度，為儲(chǔ)層描述提供了堅(jiān)實(shí)的數(shù)據(jù)基礎(chǔ).

2 方法原理

在儲(chǔ)層描述工作中，通常將縱波速度、橫波速度、體積模量、剪切模量、Gassmann流體項(xiàng)等數(shù)據(jù)稱為彈性參數(shù)；將孔隙度、泥質(zhì)含量、含水飽和度、孔隙縱橫比、壓力等數(shù)據(jù)稱為儲(chǔ)層參數(shù).通過(guò)建立彈性參數(shù)和儲(chǔ)層參數(shù)之間的聯(lián)系，可以從彈性參數(shù)中預(yù)測(cè)出儲(chǔ)層參數(shù)，這有助于評(píng)價(jià)儲(chǔ)層巖性和含油氣性.

巖石物理模型是連接儲(chǔ)層參數(shù)和彈性參數(shù)的紐帶，是儲(chǔ)層參數(shù)預(yù)測(cè)的理論基礎(chǔ).一個(gè)合理的巖石物理模型可以建立從儲(chǔ)層參數(shù)到彈性參數(shù)的一一對(duì)應(yīng)關(guān)系，這種對(duì)應(yīng)關(guān)系是巖石物理逆建模的基礎(chǔ)(Moyano et al.,2011).在這個(gè)基礎(chǔ)上，逆建模方法通過(guò)在儲(chǔ)層參數(shù)域進(jìn)行等值面空間交會(huì)建立了從彈性參數(shù)到儲(chǔ)層參數(shù)的逆向?qū)?yīng)關(guān)系，實(shí)現(xiàn)了儲(chǔ)層參數(shù)的預(yù)測(cè).

2.1 巖石物理建模

儲(chǔ)層參數(shù)和彈性參數(shù)之間的關(guān)系可以由下面的方程表示：

(1)

其中，m=(m1,m2,m3,…)是儲(chǔ)層參數(shù)，如孔隙度、含水飽和度等；d=(d1,d2,d3,…)是彈性參數(shù)，如速度、模量等；F是巖石物理模型.

進(jìn)行巖石物理建模的基本思想是：在假設(shè)巖石孔隙的特征尺度遠(yuǎn)遠(yuǎn)小于地震波長(zhǎng)的前提下，將巖石孔隙等效為具有單一縱橫比的理想橢球孔(Berryman,1992; Mukerji et al., 1995).按照彈性介質(zhì)理論，分析巖石基質(zhì)的礦物組分與體積分?jǐn)?shù)，計(jì)算基質(zhì)的彈性模量與密度；通過(guò)等效介質(zhì)理論計(jì)算巖石骨架模量，分析儲(chǔ)層環(huán)境下流體性質(zhì)，結(jié)合流體的分布特征及其與孔隙結(jié)構(gòu)的聯(lián)系，估算含流體巖石的彈性模量.

巖石物理逆建模過(guò)程不依賴特定的巖石物理模型，大多數(shù)巖石物理模型都可以根據(jù)逆建模理論建立相應(yīng)的逆建模方法.在實(shí)際應(yīng)用過(guò)程中要根據(jù)具體的儲(chǔ)層特點(diǎn)建立合適、準(zhǔn)確的模型，使得模型能夠準(zhǔn)確描述儲(chǔ)層的巖石物理特點(diǎn)，并保證利用模型重構(gòu)得到的彈性參數(shù)曲線與實(shí)際測(cè)井結(jié)果一致.本文以下面的巖石物理建模過(guò)程為例，進(jìn)行逆建模方法的說(shuō)明.具體的儲(chǔ)層巖石物理模型的構(gòu)建步驟如下：

首先利用Voigt-Reuss-Hill平均理論(Hill，1952)計(jì)算巖石基質(zhì)的等效彈性模量MVRH：

(2)

然后利用微分等效介質(zhì)模型(DEM模型)在基質(zhì)中添加孔隙，計(jì)算巖石骨架的彈性模量，等效模量K*和μ*的耦合微分方程分別是

(3)

初始條件是K*(0)=K1，μ*(0)=μ1，其中K1、μ1是利用方程(2)計(jì)算得到的基質(zhì)的體積模量和剪切模量；K2、μ2是逐漸加入的包含物的體積模量和剪切模量，y是包含物的含量，P(*2)和Q(*2)是背景介質(zhì)中包含物形狀的表達(dá)式.

最后利用Gassmann方程在巖石孔隙中充填流體(Gassmann，1951)，飽和巖石彈性模量的表達(dá)式如下：

(4)

其中，Ksat和μsat是飽和巖石彈性模量，K*和μ*是干巖石骨架彈性模量；Kf是混合流體的體積模量，使用Wood公式計(jì)算得到.

密度的計(jì)算方程為

ρ= φ((1-Sw)ρh+Swρw)

+(1-φ)(VshρC+(1-Vsh)ρQ)，

(5)

其中，ρh、ρw、ρC和ρQ分別是油(氣)、水、粘土和石英的密度.

根據(jù)前述的巖石物理建模方法，結(jié)合表1中的參數(shù)，可以建立彈性參數(shù)隨儲(chǔ)層參數(shù)的變化關(guān)系，如圖1所示.圖1a—1d分別是巖石體積模量和剪切模量隨孔隙度、泥質(zhì)含量、含水飽和度和孔隙縱橫比的變化圖.從圖中可以看出，當(dāng)其他儲(chǔ)層參數(shù)一定的時(shí)候，隨著泥質(zhì)含量的增加，體積模量和剪切模量都減??；隨著含水飽和度的增加，體積模量增加，剪切模量不受影響；隨著孔隙度的增加，體積模量和剪切模量都減小；隨著孔隙縱橫比的增加，體積模量和剪切模量都增加.巖石物理模型建立了彈性參數(shù)與儲(chǔ)層參數(shù)之間的聯(lián)系，不同的巖石物理模型建立的聯(lián)系不盡相同，根據(jù)這些關(guān)系進(jìn)行的逆建模結(jié)果也不同，因此建立準(zhǔn)確、合理的巖石物理模型至關(guān)重要.

表1 基本組分參數(shù)表Table 1 Basic constituent parameters table

圖1 彈性參數(shù)隨儲(chǔ)層參數(shù)變化圖(a) 體積模量隨泥質(zhì)含量和含水飽和度的變化； (b) 剪切模量隨泥質(zhì)含量和含水飽和度的變化； (c) 體積模量隨孔隙度和孔隙縱橫比的變化； (d) 剪切模量隨孔隙度和孔隙縱橫比的變化.Fig.1 The variation of elastic parameters with reservoir parameters(a) The variation of bulk modulus with clay content and water saturation; (b) The variation of shear modulus with clay content and water saturation; (c) The variation of bulk modulus with porosity and aspect ratio; (d) The variation of shear modulus with porosity and aspect ratio.

2.2 巖石物理逆建模

由巖石物理建模過(guò)程可以看出，巖石物理模型表示的巖石彈性參數(shù)與儲(chǔ)層參數(shù)之間的關(guān)系是強(qiáng)非線性的，要推導(dǎo)其逆運(yùn)算的解析表達(dá)式、直接計(jì)算出儲(chǔ)層參數(shù)是很困難的，因此需要利用逆建模方法間接地完成預(yù)測(cè).

2.2.1 三維巖石物理逆建模

巖石物理逆建模是在儲(chǔ)層參數(shù)域中，根據(jù)巖石物理模型建立彈性參數(shù)的單值等值面，通過(guò)等值面的空間交會(huì)得到交點(diǎn)坐標(biāo)，完成儲(chǔ)層參數(shù)的預(yù)測(cè).在三維儲(chǔ)層參數(shù)域，需要三個(gè)單值等值面的交會(huì)才能得到一個(gè)穩(wěn)定的交點(diǎn).三維巖石物理逆建模的具體步驟如下：

以一組實(shí)測(cè)彈性參數(shù)d0(K0,μ0,ρ0,f0,M0)作為逆建模的輸入數(shù)據(jù)，要保證反問(wèn)題是非欠定的，需要至少有三種輸入數(shù)據(jù)，這里選取的是體積模量K0、Gassmann流體項(xiàng)f0和密度ρ0(Russell et al.，2003).

選取孔隙度φ、含水飽和度Sw和等效孔隙縱橫比α建立一個(gè)三維直角坐標(biāo)系，稱為一個(gè)儲(chǔ)層參數(shù)域.分別以φ,α,Sw為x軸、y軸和z軸，取值范圍分別是(0～0.4)，(0.001～1)和(0～1)，各坐標(biāo)軸均勻采樣.則方程(1)可以表示為

(6)

?φ∈(φmin,φmax),Sw∈(Swmin,Swmax),α∈(αmin,αmax)其中，φmin,φmax、Swmin,Swmax和αmin,αmax分別是孔隙度、含水飽和度和等效孔隙縱橫比的最小值和最大值.

儲(chǔ)層參數(shù)域中每個(gè)點(diǎn)的坐標(biāo)是一組儲(chǔ)層參數(shù)值，利用建立的巖石物理模型計(jì)算出每一個(gè)點(diǎn)對(duì)應(yīng)的彈性參數(shù)值，可得到一個(gè)三維的數(shù)據(jù)場(chǎng).圖2是體積模量的三維數(shù)據(jù)場(chǎng)，這個(gè)數(shù)據(jù)場(chǎng)中包含所有可能的體積模量值.根據(jù)相同的原理可以建立Gassmann流體項(xiàng)和密度的三維數(shù)據(jù)場(chǎng).

圖2 儲(chǔ)層參數(shù)域的體積模量數(shù)據(jù)場(chǎng)Fig.2 The bulk modulus data cube in the reservoir parameter domain

這里采用Marching Cubes算法求取數(shù)據(jù)場(chǎng)中的單值等值面(Lorensen and Cline，1987).算法的主要思想是將三維數(shù)據(jù)場(chǎng)劃分為多個(gè)小的立方體，作為最小搜索單元(稱為體元).逐個(gè)處理數(shù)據(jù)場(chǎng)的體元，選出與等值面交會(huì)的體元，采用線性插值計(jì)算出與等值面交會(huì)的體元的棱邊上的交點(diǎn)，根據(jù)體元8個(gè)頂點(diǎn)與等值面的相對(duì)位置，將這些交點(diǎn)連接成體元內(nèi)的等值面，所有體元中的等值面構(gòu)成整個(gè)數(shù)據(jù)場(chǎng)的等值面.在等值面的空間結(jié)構(gòu)比較簡(jiǎn)單的情況下，經(jīng)典的Marching Cubes算法就可以完成；在等值面的空間構(gòu)型較為復(fù)雜情況下，需要使用改進(jìn)的Marching Cubes算法(梁秀霞等，2006)，以保證曲面擬合的精度和效率.

利用Marching Cubes算法計(jì)算體積模量為K0的單值等值面(這里用ISK0(φ,Sw,α)來(lái)表示)，如圖3a所示.等值面在坐標(biāo)軸的投影寬度代表該彈性參數(shù)對(duì)該儲(chǔ)層參數(shù)的約束能力，投影越窄約束能力越強(qiáng)，預(yù)測(cè)得到的儲(chǔ)層參數(shù)越準(zhǔn)確.同理可以得到密度為ρ0在儲(chǔ)層參數(shù)域的等值面(用ISρ0(φ,Sw,α)來(lái)表示)，如圖3b所示.將體積模量的單值等值面和密度的單值等值面交會(huì)在同一個(gè)儲(chǔ)層參數(shù)域中，即可找到體積模量和密度值均與實(shí)測(cè)值相等的點(diǎn).如圖4a所示，這些點(diǎn)就分布在兩個(gè)等值面的交線上，可以通過(guò)求解方程ISK0(φ,Sw,α)=ISρ0(φ,Sw,α)獲取交點(diǎn)坐標(biāo).利用模擬退火算法求取交線上的點(diǎn)坐標(biāo)，只需建立如方程(7)的目標(biāo)函數(shù)，在兩個(gè)單值等值面的定義域交集內(nèi)求取目標(biāo)函數(shù)的最小值即可.

圖3 儲(chǔ)層參數(shù)域的單值等值面(a) 體積模量； (b) 密度； (c) Gassmann流體項(xiàng).Fig.3 Single isosurfaces in the reservoir parameter domain(a) Bulk modulus； (b) Density； (c) Gassmann fluid term.

圖4 不同單值等值面的交會(huì)(a) 兩個(gè)等值面的交會(huì)； (b) 三個(gè)等值面的交會(huì).Fig.4 Intersections of different isosurfaces(a) 2 isosurfaces；(b) 3 isosurfaces.

J1=(φK0-φρ0)2+(SwK0-Swρ0)2+(αK0-αρ0)2，

(7)

其中，φK0,SwK0,αK0,φρ0,Swρ0,αρ0是單值等值面ISK0(φ,Sw,α)中點(diǎn)的坐標(biāo)和單值等值面ISρ0(φ,Sw,α)中點(diǎn)的坐標(biāo).

計(jì)算出Gassmann流體項(xiàng)為f0在儲(chǔ)層參數(shù)域的單值等值面(用ISf0(φ,Sw,α)來(lái)表示)，如圖3c所示.將其與圖4a中的兩個(gè)單值等值面交會(huì)在一起，結(jié)果如圖4b所示，三個(gè)等值面的交集部分是圖中黑色圓圈標(biāo)出的點(diǎn)，這個(gè)點(diǎn)的體積模量、Gassmann流體項(xiàng)和密度的值均與實(shí)測(cè)值(K0,f0,ρ0)相等，它的坐標(biāo)即所求儲(chǔ)層參數(shù)，可以通過(guò)求解方程ISK0(φ,Sw,α)=ISf0(φ,Sw,α)=ISρ0(φ,Sw,α)獲取交點(diǎn)坐標(biāo).利用模擬退火算法求取交點(diǎn)坐標(biāo)，只需建立如方程(8)的目標(biāo)函數(shù)，在三個(gè)單值等值面的定義域交集內(nèi)求取目標(biāo)函數(shù)的最小值即可，得到的最終解是(φ=0.19,Sw=0.49,α=0.30).

J2=(φK0-φρ0)2+(SwK0-Swρ0)2+(αK0-αρ0)2

+(φK0-φf(shuō)0)2+(SwK0-Swf0)2+(αK0-αf0)2,

(8)

其中，φf(shuō)0,Swf0,αf0是等值面ISf0(φ,Sw,α)中的點(diǎn)的坐標(biāo).

圖5 不同初值情況下得到的解的相對(duì)誤差Fig.5 Relative errors under conditions of different initial values

這里使用的模擬退火算法是非?？焖倌M退火算法(Ingber，1989；陳華根等，2004)，通過(guò)不斷的試驗(yàn)，調(diào)整算法的參數(shù)，使之達(dá)到最優(yōu)，不僅有較高的運(yùn)算效率，而且不易陷入局部極小值.為了研究初始值對(duì)于本方法的影響，以d0(K0,f0,ρ0)輸入數(shù)據(jù)，選取了一系列不同的初始值用于三維巖石物理逆建模中的模擬退火運(yùn)算，進(jìn)行10000次模擬退火運(yùn)算得到的結(jié)果取平均值，在10000次模擬結(jié)果中均未出現(xiàn)局部極值的情況.圖5是計(jì)算結(jié)果的相對(duì)誤差隨初始值的變化關(guān)系，圖中的橫軸是這些點(diǎn)與真實(shí)值所在坐標(biāo)的距離，距離越小說(shuō)明初始值與真實(shí)值更接近，縱軸是估計(jì)得到的儲(chǔ)層參數(shù)的相對(duì)誤差，其中黑色線、粉色線和綠色線分別表示孔隙度、孔隙縱橫比和含水飽和度的相對(duì)誤差.從圖中可以看出，初始值越接近真實(shí)值，誤差越小，因此在實(shí)際應(yīng)用中使用的初始值應(yīng)該盡量準(zhǔn)確，這不僅有助于提高準(zhǔn)確性，還可以減少迭代次數(shù)，提高收斂的效率.在井位置處可以利用測(cè)井得到的儲(chǔ)層參數(shù)曲線經(jīng)過(guò)平滑處理后的結(jié)果作為初始值，在井旁道位置處可以利用測(cè)井得到的儲(chǔ)層參數(shù)曲線進(jìn)行插值外推的結(jié)果作為初始值.

2.2.2 擴(kuò)展的巖石物理逆建模

為了給儲(chǔ)層特征描述工作提供更多的數(shù)據(jù)，基于巖石物理逆建模理論可以進(jìn)行四個(gè)儲(chǔ)層參數(shù)的預(yù)測(cè).四個(gè)參數(shù)的預(yù)測(cè)對(duì)于前述的巖石物理逆建模方法來(lái)說(shuō)是一種欠定的反問(wèn)題，無(wú)法得到唯一解，因此建立擴(kuò)展的巖石物理逆建模方法：選擇要預(yù)測(cè)的四個(gè)儲(chǔ)層參數(shù)：泥質(zhì)含量Vsh、孔隙度φ、含水飽和度Sw和等效孔隙縱橫比α，可以建立四個(gè)不同的三維參數(shù)域(φ,Sw,α)、(φ,Vsh,α)、(φ,Sw,Vsh)和(Vsh,Sw,α)，在其中選擇兩個(gè)參數(shù)域(φ,Sw,α)和(φ,Vsh,α)；Vsh和Sw在各自的定義域內(nèi)均勻取值，針對(duì)每個(gè)Vsh的不同取值在參數(shù)域(φ,Sw,α)內(nèi)、針對(duì)每個(gè)Sw的不同取值在參數(shù)域(φ,Vsh,α)內(nèi)分別進(jìn)行三維巖石物理逆建模，將兩組逆建模結(jié)果綜合交會(huì)進(jìn)行預(yù)測(cè).

這樣通過(guò)逆建模得到了n個(gè)不同泥質(zhì)含量情況下的交點(diǎn).將m1(φ,Sw,α,Vsh)投影到儲(chǔ)層參數(shù)域(φ,Sw,α)中，通過(guò)插值得到曲線m1(φ,Sw,α)，投影方程表示為

m1(φ,Sw,α,Vsh)T=m1(φ,Sw,α)，

(9)

其中，T是4×3的投影矩陣，由于坐標(biāo)軸Vsh與參數(shù)域(φ,Sw,α)是正交的，因此投影矩陣中的元素t11=t22=t33=1，其他元素都為0；m1(φ,Sw,α)表示m1(φ,Sw,α,Vsh)投影以后的曲線，這條線包含最終的解.

(10)

(11)

其中，φm1,φm2,Swm1,Swm2,αm1,αm2分別表示曲線m1(φ,Sw,α)和m2(φ,Sw,α)上的點(diǎn)坐標(biāo).

圖6 不同泥質(zhì)含量情況下的體積模量數(shù)據(jù)場(chǎng)示意圖Fig.6 Data cubes of bulk modulus with different clay contents

圖7 兩個(gè)儲(chǔ)層參數(shù)域的逆建模結(jié)果的投影的交會(huì)圖Fig.7 Intersection of projections of inverse modeling results from two reservoir parameter domains

2.3 輸入數(shù)據(jù)的適用性分析

可以用于巖石物理逆建模的彈性參數(shù)有許多，在進(jìn)行儲(chǔ)層參數(shù)預(yù)測(cè)時(shí)，選取合適的彈性參數(shù)是預(yù)測(cè)成功的關(guān)鍵步驟之一.對(duì)彈性參數(shù)數(shù)據(jù)的選取主要取決于在建立的巖石物理模型框架下參數(shù)的單值等值面的形狀，相互平行的單值等值面沒(méi)有交線，不平行的單值等值面容易得到交線，越趨于正交的單值等值面越有利于獲得準(zhǔn)確的最終解.兩個(gè)單值等值面之間的正交性可以表示為

(12)

圖8 儲(chǔ)層參數(shù)預(yù)測(cè)流程圖Fig.8 Flowchart of reservoir parameter prediction

建立如下的適用性因子來(lái)衡量一組彈性參數(shù)的適用性：

(13)

其中，TOF是一組彈性參數(shù)的適用性因子，OFdidj是彈性參數(shù)組中任意兩種不同彈性參數(shù)的正交性因子，i=1,2,…n，j=1,2,…n，i≠j，N是正交性因子的個(gè)數(shù).

方程(13)是一個(gè)參數(shù)組中各參數(shù)等值面之間的正交性因子的均值，若TOF趨近于1，則這個(gè)參數(shù)組中的各彈性參數(shù)的等值面之間的正交性差，易產(chǎn)生多個(gè)解；若TOF趨近于0，則各彈性參數(shù)的等值面之間的正交性強(qiáng)，易找到唯一解.

圖9 原始模型(a) 泥質(zhì)含量； (b) 孔隙度； (c) 等效孔隙縱橫比； (d) 含水飽和度.Fig.9 Original model(a) Clay content； (b) Porosity； (c) Effective pore aspect ratio； (d) Water saturation.

利用方程(13)計(jì)算彈性參數(shù)組的適用性因子，可以在進(jìn)行巖石物理逆建模之前就了解各組參數(shù)的適用性.整個(gè)基于逆建模的儲(chǔ)層參數(shù)預(yù)測(cè)流程如圖8所示.

3 模型試算

為了檢驗(yàn)基于逆建模的儲(chǔ)層特征定量預(yù)測(cè)方法的有效性，建立如圖9所示的模型，其中有54個(gè)點(diǎn)，對(duì)應(yīng)的泥質(zhì)含量、孔隙度、含水飽和度和等效孔隙縱橫比分別如圖9a—9d所示.然后根據(jù)前述的巖石物理建模方法可以計(jì)算出多個(gè)彈性參數(shù)曲線.為了找到合理的彈性參數(shù)，首先根據(jù)方程(13)計(jì)算不同參數(shù)組合的適用性，分別計(jì)算模型中每個(gè)點(diǎn)的適用性因子，然后取平均值，得到的結(jié)果如圖10所示.圖10中橫坐標(biāo)1到5分別表示M/μ,K,ρ,f、M/μ,Vp,ρ,μ、K,f,μ,Vp/Vs、Ip,Is,K,μ和Vp,Vs,ρ,K這些參數(shù)的組合，縱坐標(biāo)為適用性因子的值.

通過(guò)比較可以看出，體積模量K、Gassmann流體項(xiàng)f、密度ρ和模量比M/μ是適用性較強(qiáng)的組合，因此可以利用這些彈性參數(shù)進(jìn)行儲(chǔ)層參數(shù)預(yù)測(cè)，其曲線如圖11所示.利用這四個(gè)參數(shù)進(jìn)行擴(kuò)展的巖石物理逆建模，得到的預(yù)測(cè)結(jié)果如圖12所示.圖12中空心點(diǎn)是圖9中的原始模型數(shù)據(jù)，實(shí)心點(diǎn)是預(yù)測(cè)出的儲(chǔ)層參數(shù)值，從圖中可以看出，四個(gè)儲(chǔ)層參數(shù)的值都較為準(zhǔn)確，相比而言含水飽和度和等效孔隙縱橫比的精度略差.

圖10 彈性參數(shù)組合的適用性對(duì)比Fig.10 Comparison of applicability of different combinations of elastic parameters

圖11 彈性參數(shù)曲線(a) 體積模量； (b) Gassmann流體項(xiàng)； (c) 密度； (d) 模量比.Fig.11 Logs of elastic parameters(a) Bulk modulus； (b) Gassmann fluid term； (c) Density； (d) Moduli ratio.

圖12 預(yù)測(cè)結(jié)果與原始模型的對(duì)比(a) 泥質(zhì)含量； (b) 孔隙度； (c) 等效孔隙縱橫比； (d) 含水飽和度.Fig.12 Comparison of predictions and original data(a) Clay content； (b) Porosity； (c) Effective pore aspect ratio； (d) Water saturation.

4 實(shí)際應(yīng)用

將本文所述方法應(yīng)用于中國(guó)東部某油田的實(shí)際資料.綜合工區(qū)的測(cè)井、錄井、地質(zhì)和實(shí)驗(yàn)數(shù)據(jù)，建立合理的巖石物理模型，并選取合適的基礎(chǔ)參數(shù).首先利用該工區(qū)的一口井(A井)的數(shù)據(jù)進(jìn)行擴(kuò)展的巖石物理逆建模，結(jié)果如圖13所示.圖13中從左到右依次是泥質(zhì)含量、孔隙度、含水飽和度和等效孔隙縱橫比，其中實(shí)線是測(cè)井?dāng)?shù)據(jù)，虛線是預(yù)測(cè)結(jié)果.通過(guò)對(duì)比可以看出，預(yù)測(cè)結(jié)果較為準(zhǔn)確.

圖13 利用井?dāng)?shù)據(jù)的預(yù)測(cè)結(jié)果與實(shí)際測(cè)量數(shù)據(jù)的對(duì)比Fig.13 Comparison of predictions from well data and the measured data

接下來(lái)進(jìn)行地震反演和基于地震反演結(jié)果的擴(kuò)展的巖石物理逆建模.圖14是過(guò)工區(qū)中A井的地震剖面，進(jìn)行疊前地震反演之前，需要對(duì)地震數(shù)據(jù)進(jìn)行保幅處理.通過(guò)彈性參數(shù)組的適用性分析發(fā)現(xiàn)體積模量、Gassmann流體項(xiàng)、密度和模量比M/μ是適用性較好的參數(shù)組合，可以作為輸入數(shù)據(jù).

隨后利用疊前反演得到以上四種彈性參數(shù)，并進(jìn)行擴(kuò)展的巖石物理逆建模，儲(chǔ)層參數(shù)預(yù)測(cè)結(jié)果如圖15所示.圖15a—15d分別是預(yù)測(cè)得到的泥質(zhì)含量、孔隙度、含水飽和度和等效孔隙縱橫比剖面，在圖15a—15c中井位置處的黑色井曲線分別是對(duì)應(yīng)的實(shí)際測(cè)井曲線.從圖15中看出，2.4 s、2.5 s和2.65 s深度處有三層砂巖，分別稱其為第一層、第二層和第三層砂巖.其中，第一層和第三層砂巖的泥質(zhì)含量較低、等效孔隙縱橫比較高，是較好的砂巖儲(chǔ)層，而第二層砂巖雖然孔隙度高，但是泥質(zhì)含量偏高，不是好的砂巖儲(chǔ)層.第一層和第二層的含水飽和度值偏高，不是含油氣儲(chǔ)層，第三層的含水飽和值較低，發(fā)育的是有利儲(chǔ)層.實(shí)際測(cè)井解釋結(jié)果和鉆井結(jié)果顯示，這三個(gè)位置鉆遇砂巖，第一層和第二層是含水砂巖層，第三層是含油砂巖層.通過(guò)比較可以看到，預(yù)測(cè)結(jié)果和實(shí)際數(shù)據(jù)匹配較好，證明方法具有一定的實(shí)用性.

圖14 過(guò)A井的地震剖面Fig.14 Seismic profile across well A

5 結(jié)論

儲(chǔ)層參數(shù)是我們理解巖石的幾何特性和動(dòng)力學(xué)特性的必備要素，也是油氣勘探中評(píng)價(jià)儲(chǔ)層的重要依據(jù)，因此進(jìn)行儲(chǔ)層參數(shù)預(yù)測(cè)具有重要意義.本文基于巖石物理理論研究了三維巖石物理逆建模方法，通過(guò)在儲(chǔ)層參數(shù)域建立彈性參數(shù)的單值等值面并進(jìn)行等值面的空間交會(huì)，完成了儲(chǔ)層參數(shù)的預(yù)測(cè).在此基礎(chǔ)上，提出了擴(kuò)展的巖石物理逆建模方法預(yù)測(cè)四個(gè)儲(chǔ)層參數(shù).巖石物理逆建模是一種靈活性很強(qiáng)的方法，使用的巖石物理模型和輸入?yún)?shù)都不是固定的，可以根據(jù)具體的目標(biāo)儲(chǔ)層特征有針對(duì)性地選擇，這樣有利于提高預(yù)測(cè)結(jié)果的精確程度.經(jīng)過(guò)模型數(shù)據(jù)的檢驗(yàn)和實(shí)際工區(qū)的應(yīng)用，本方法取得了不錯(cuò)的應(yīng)用效果，說(shuō)明了方法具有一定的實(shí)用性.

圖15 儲(chǔ)層參數(shù)預(yù)測(cè)結(jié)果剖面(a) 泥質(zhì)含量； (b) 孔隙度； (c) 含水飽和度； (d) 等效孔隙縱橫比.Fig.15 Profiles of predicted reservoir parameters(a) Clay content； (b) Porosity； (c) Water saturation； (d) Effective pore aspect ratio.

該方法的前提是有合理的巖石物理模型，因此在實(shí)際應(yīng)用中巖石物理模型的建立至關(guān)重要，模型和參數(shù)的選擇要盡可能的準(zhǔn)確、合理.為了提高預(yù)測(cè)的效率和結(jié)果的精度，建立了彈性參數(shù)組適用性評(píng)價(jià)方法，這種方法可以在儲(chǔ)層參數(shù)預(yù)測(cè)之前就優(yōu)選出適用性強(qiáng)的參數(shù)組合.類似地，在其他類型的儲(chǔ)層的描述中也可以建立類似的儲(chǔ)層參數(shù)預(yù)測(cè)方法.

Avseth P A, ?degaard E. 2004. Well log and seismic data analysis using rock physics templates.FirstBreak, 22(10): 37-43.

Bachrach R. 2006. Joint estimation of porosity and saturation using stochastic rock-physics modeling.Geophysics, 71(5): O53- O63. Berryman J G. 1992. Single-scattering approximations for coefficients in Biot′s equations of poroelasticity.TheJournaloftheAcousticalSocietyofAmerica, 91(2): 551-571. Bosch M, Mukerji T, Gonzalez E F. 2010. Seismic inversion for reservoir properties combining statistical rock physics and geostatistics: A review.Geophysics, 75(5): 75A165-75A176.

Chen H G, Wu J S, Wang J L, et al. 2004. Mechanism study of simulated annealing algorithm.JournalofTongjiUniversity(NaturalScience), 32(6): 802-805.

Doyen P M. 1988. Porosity from seismic data: A geostatistical approach.Geophysics, 53(10): 1263-1275.

Eidsvik J, Avseth P, More H, et al. 2004. Stochastic reservoir characterization using prestack seismic data.Geophysics, 69(4): 978-993. Gassmann F. 1951. Uber die elastizitat, poroser medien. Vierteljahrsschrift der NaturforschendenGesellschaftinZurich, 96: 1-23.Grana D, Rossa E D. 2010. Probabilistic petrophysical-properties estimation integrating statistical rock physics with seismic inversion.Geophysics, 75(3): O21-O37.

Hill R. 1952. The elastic behaviour of a crystalline aggregate. Proceeding of the Physical Society.SectionA, 65(5): 349-354.

Hu H F, Yin X Y, Wu G C. 2012. Joint inversion of petrophysical parameters based on Bayesian classification.GeophysicalProspectingforPetroleum(in Chinese), 51(3): 225-232.

Ingber L. 1989. Very fast simulated re-annealing.MathematicalandComputerModelling, 12(8): 967-973. Jiang M J, Spikes K. 2012. Estimation of the porosity and pore aspect ratio of the Haynesville Shale using the Self-Consistent Model and a Grid Search Method. ∥ 82th Annual International Meeting, Society of Exploration Geophysicists. Expanded Abstracts, 1-5. Li H B, Zhang J J, Yao F C. 2013. Inversion of effective pore aspect ratios for porous rocks and its applications.ChineseJ.Geophys. (in Chinese), 56(2): 608-615, doi: 10.6038/cjg20130224.

Liang X X, Zhang C M, Liu Y, et al. 2006. A topology complexity based method to approximate isosurface with trilinear interpolated triangular patch.JournalofComputerResearchandDevelopment(in Chinese), 43(3): 528-535. Lorensen W E, Cline H E. 1987. Marching cubes: a high resolution 3D surface construction algorithm.ACMSIGGRAPHComputerGraphics, 21(4): 163-169. Moyano B, Jensen E H, Johansen T A. Improved quantitative calibration of rock physics models.PetroleumGeoscience, 2011, 17(4): 345-354.Mukerji T, Berryman J, Mavko G, et al. 1995. Differential effective medium modeling of rock elastic moduli with critical porosity constraints.GeophysicalResearchLetters, 22(5): 555-558.Mukerji T, Avseth P, Mavko G, et al. 2001. Statistical rock physics: Combining rock physics, information theory, and geostatistics to reduce uncertainty in seismic reservoir characterization.TheLeadingEdge, 20(3): 313-319.

Russell B H, Hedlin K, Hilterman F J, et al. 2003. Fluid-property discrimination with AVO: A Biot-Gassmann perspective.Geophysics, 68(1): 29-39. Spikes K, Mukerji T, Dvorkin J, et al. 2007. Probabilistic seismic inversion based on rock-physics models.Geophysics, 72(5): R87-R97.

Yang W Y, Wang C B. 2010. Utilizing pre-stack simultaneous inversion to predict reservoir physical properties.OilGeophysicalProspecting(in Chinese), 45(3): 414-417.Yin X Y, Cui W, Zong Z Y, et al. 2014. Petrophysical property inversion of reservoirs based on elastic impedance.ChineseJ.Geophys. (in Chinese), 57(12): 4132-4140, doi: 10.6038/cjg20141224.

Yin X Y, Sun R Y, Wang B L, et al. 2014. Simultaneous inversion of petrophysical parameters based on geostatistical a priori information.AppliedGeophysics, 11(3): 311-320.

附中文參考文獻(xiàn)

陳華根, 吳健生, 王家林等. 2004. 模擬退火算法機(jī)理研究. 同濟(jì)大學(xué)學(xué)報(bào)(自然科學(xué)版), 32(6): 802-805.

胡華鋒, 印興耀, 吳國(guó)忱. 2012. 基于貝葉斯分類的儲(chǔ)層物性參數(shù)聯(lián)合反演方法. 石油物探, 51(3): 225-232.

李宏兵, 張佳佳, 姚逢昌. 2013. 巖石的等效孔隙縱橫比反演及其應(yīng)用. 地球物理學(xué)報(bào), 56(2): 608-615, doi: 10.6038/cjg20130224. 梁秀霞, 張彩明, 劉毅等. 2006. 拓?fù)浣Y(jié)構(gòu)正確的三線性插值曲面的三角片逼近. 計(jì)算機(jī)研究與發(fā)展, 43(3): 528-535.

楊午陽(yáng), 王叢鑌. 2010. 利用疊前AVA同步反演預(yù)測(cè)儲(chǔ)層物性參數(shù). 石油地球物理勘探, 45(3): 414-417.

印興耀, 崔維, 宗兆云等. 2014. 基于彈性阻抗的儲(chǔ)層物性參數(shù)預(yù)測(cè)方法. 地球物理學(xué)報(bào), 57(12): 4132-4140, doi: 10.6038/cjg20141224.

(本文編輯 胡素芳)

Quantitative prediction of reservoir characteristics based on inverse modeling theory

LIU Qian1, YIN Xing-Yao1, LI Chao2

1SchoolofGeosciences,ChinaUniversityofPetroleum(Huadong)，Qingdao266580,China2CNOOCResearchInstitute,Beijing100028，China

Reservoir parameters are used to describe the reservoir characteristics quantitatively. It′s important to predict reservoir parameters in reservoir characterizations. This paper proposes a method to predict reservoir parameters based on inverse modeling theory. In the framework of rock physics theory, this paper focuses on parameter prediction with inverse modeling in the reservoir parameter domain. First of all, an accurate rock physics model was established according to the characteristic of the target reservoir. Then, in the reservoir parameter domain, elastic parameters at every point were calculated to form elastic data cubes. Isosurfaces of the elastic data cubes can be calculated with the Marching Cubes algorithm. The reservoir parameter prediction was implemented by the isosurface intersection of three different elastic parameters in the same domain. Based on these, an extended inverse modeling method has also been proposed. It contains two main parts-inverse modelings in two different reservoir parameter domains and the domain transformations of the inverse modeling results. Different kinds of elastic parameters make different contributions to the accuracy of inverse modeling. An applicability evaluation method was proposed for the choice of elastic parameters as input data. A model data and a real example have been chosen to test the inverse modeling method.

Reservoir parameter; Rock physics; Inverse modeling; Reservoir parameter domain

10.6038/cjg20160931.

國(guó)家重點(diǎn)基礎(chǔ)研究發(fā)展計(jì)劃(973計(jì)劃)項(xiàng)目(2013CB228604)資助.

劉倩，1988年生，女，博士研究生，主要從事地震巖石物理建模及地震波場(chǎng)模擬方面的研究. E-mail：liuqian88115@163.com

10.6038/cjg20160931

P631

2015-11-21，2016-01-11收修定稿

劉倩, 印興耀, 李超. 2016. 基于逆建模理論的儲(chǔ)層特征定量預(yù)測(cè)方法. 地球物理學(xué)報(bào)，59(9):3491-3502,

Liu Q, Yin X Y, Li C. 2016. Quantitative prediction of reservoir characteristics based on inverse modeling theory.ChineseJ.Geophys. (in Chinese),59(9):3491-3502,doi:10.6038/cjg20160931.

In the model test, there are 54 points with different reservoir properties in the original model. A rock physics model was established and several elastic parameters were calculated. The prediction result was consistent with the original model data. In the real example, the proposed method was applied to a set of data from Eastern China. A site-specific rock physics model was built according to the well data. The inverse modeling method was first validated by the application on the well data. Inverted elastic parameters were also taken as the input data of inverse modeling. The accuracy of the inverse modeling with inverted elastic data was also high enough. Based on the theory of rock physics modeling, an extended inverse modeling method in reservoir parameter domain has been proposed to predict reservoir parameters quantitatively in this paper. The accuracy of the rock physics model is the foundation of the inverse modeling. The choice of input data is very important to the stability of the inverse modeling results. So the applicability of the input data should be evaluated before the inverse modeling. Tests on model data and real data have proved the rationality and applicability of this method in the prediction of reservoir parameters.