基于線性化巖石物理反演的物性參數(shù)預(yù)測方法

張佳佳，印興耀，張廣智，谷一鵬，樊祥剛

（1.中國石油大學(xué)（華東）地球科學(xué)與技術(shù)學(xué)院，山東青島 266580；2.深層油氣地質(zhì)與勘探教育部重點(diǎn)實(shí)驗(yàn)室，山東青島 266580；3.海洋國家實(shí)驗(yàn)室海洋礦產(chǎn)資源評價(jià)與探測技術(shù)功能實(shí)驗(yàn)室，山東青島 266071）

0 引言

地震數(shù)據(jù)中蘊(yùn)含著豐富的地震彈性參數(shù)和儲(chǔ)集層物性參數(shù)信息[1]，儲(chǔ)集層物性參數(shù)預(yù)測主要分兩步：第1步是利用常規(guī)地震反演方法由地震數(shù)據(jù)獲取地震彈性參數(shù)，現(xiàn)有的地震反演技術(shù)已經(jīng)非常成熟[2-5]，能夠反演得到較為準(zhǔn)確的彈性參數(shù)，例如縱橫波速度、密度等；第 2步是利用巖石物理反演由地震彈性參數(shù)預(yù)測儲(chǔ)集層參數(shù)，如孔隙度、飽和度、泥質(zhì)含量等[6-8]。巖石物理反演首先需要構(gòu)建一個(gè)理論的或經(jīng)驗(yàn)的巖石物理模型，即建立地震彈性參數(shù)與儲(chǔ)集層物性參數(shù)之間的關(guān)系[9-10]。其次就是巖石物理反演采用的反演算法，現(xiàn)階段絕大多數(shù)采用了概率統(tǒng)計(jì)學(xué)反演方法等[11-12]。

巖石物理模型和巖石物理反演算法都直接影響儲(chǔ)集層物性參數(shù)反演的精度。例如在許多巖石物理模型中，接觸膠結(jié)模型[13]數(shù)學(xué)形式相對簡單，被廣泛應(yīng)用于巖石物理反演中[14-15]。但是實(shí)際生產(chǎn)中彈性參數(shù)與物性參數(shù)之間的關(guān)系非常復(fù)雜，例如利用等效介質(zhì)理論[16]來構(gòu)建巖石物理模型。如果選擇復(fù)雜巖石物理模型則需要采用概率統(tǒng)計(jì)學(xué)反演算法，如貝葉斯方法或蒙特卡洛算法等。國內(nèi)外很多學(xué)者在這兩方面開展了大量研究工作，如Bachrach等建立了物性參數(shù)與彈性參數(shù)之間的統(tǒng)計(jì)關(guān)系，采用貝葉斯反演方法實(shí)現(xiàn)了孔隙度和含水飽和度的聯(lián)合反演[17]。Spikes等利用接觸膠結(jié)模型建立了物性參數(shù)與彈性參數(shù)之間的關(guān)系，采用貝葉斯反演方法實(shí)現(xiàn)了孔隙度、泥質(zhì)含量和含水飽和度的預(yù)測[18]。Bosch利用Wyllie時(shí)間平均方程[14]對巖石進(jìn)行建模，采用馬爾可夫鏈蒙特卡洛抽樣方法和最小二乘迭代反演物性參數(shù)[19]。鄧?yán)^新等提出了基于隨機(jī)巖石物理模型的貝葉斯反演方法，實(shí)現(xiàn)對孔隙度與含水飽和度進(jìn)行聯(lián)合反演[20]。Grana和Rossa利用接觸膠結(jié)模型[13]進(jìn)行巖石物理建模，在貝葉斯理論框架下進(jìn)行儲(chǔ)集層參數(shù)后驗(yàn)概率分布估計(jì)，實(shí)現(xiàn)儲(chǔ)集層參數(shù)反演[21]。印興耀等建立了彈性阻抗與儲(chǔ)集層物性參數(shù)之間的統(tǒng)計(jì)關(guān)系，利用貝葉斯理論估計(jì)物性參數(shù)的后驗(yàn)概率分布，實(shí)現(xiàn)儲(chǔ)集層物性參數(shù)反演[11]。Grana提出了對膠結(jié)接觸模型線性近似的方法，利用貝葉斯理論求解巖石物理反演問題，估計(jì)儲(chǔ)集層物性參數(shù)[22]。林恬等利用包含物模型構(gòu)建巖石物理模型，提出了約束最小二乘法與信賴域的儲(chǔ)集層參數(shù)地震反演方法[12]。Li等利用等效介質(zhì)理論[23]構(gòu)建三維巖石物理模版，采用最小二乘原理尋找三維模版網(wǎng)格點(diǎn)的方法預(yù)測物性參數(shù)。楊培杰建立了砂泥巖儲(chǔ)集層物性參數(shù)和地震彈性參數(shù)之間的定量關(guān)系，基于貝葉斯理論同步反演儲(chǔ)集層孔隙度和含水飽和度[24]。上述巖石物理反演方法選用的巖石物理模型比較簡單，例如統(tǒng)計(jì)經(jīng)驗(yàn)關(guān)系和膠結(jié)接觸模型，另一個(gè)問題就是很多反演算法利用的是概率統(tǒng)計(jì)學(xué)反演算法，需要進(jìn)行全局尋優(yōu)或者隨機(jī)抽樣，求解巖石物理反演問題的計(jì)算效率較低。因此，本文提出了一種巖石物理模型線性化方法，對巖石物理模型進(jìn)行泰勒展開，得到巖石物理反演問題的一階近似表達(dá)式。然后再利用阻尼最小二乘算法直接求解巖石物理反演問題，獲得巖石物理反演問題的解析解，不需要全局尋優(yōu)或者隨機(jī)抽樣，而是直接求逆運(yùn)算，計(jì)算效率較高。本文利用 Gassmann方程[25]和臨界孔隙度模型[26]對含泥質(zhì)砂巖儲(chǔ)集層進(jìn)行巖石物理建模，建模過程比較復(fù)雜，構(gòu)建的地震彈性參數(shù)和儲(chǔ)集層物性參數(shù)之間的關(guān)系是非線性的。實(shí)際測井?dāng)?shù)據(jù)和地震數(shù)據(jù)驗(yàn)證表明，線性化巖石物理反演方法可以獲得較為準(zhǔn)確的物性參數(shù)，可以用來近似實(shí)際巖石物理模型。

1 方法原理

1.1 線性化巖石物理反演

巖石物理反演問題一般可寫為：

（1）式中，d為彈性參數(shù)，m為待預(yù)測的物性參數(shù)，f為彈性參數(shù)與物性參數(shù)之間的關(guān)系式，即巖石物理模型。為簡單起見，首先假設(shè)（1）式中所有變量均為單一參數(shù)的標(biāo)量，然后再把它們推廣到多參數(shù)，例如d為縱波速度，m為孔隙度，f為Wyllie時(shí)間平均方程[19]或者Raymer改進(jìn)公式[27]。

實(shí)際生產(chǎn)中，巖石物理模型f并不像Wyllie時(shí)間平均方程或者Raymer改進(jìn)公式這樣簡單，往往是非常復(fù)雜而且是非線性的。地球物理中經(jīng)常利用泰勒展開來近似復(fù)雜的非線性函數(shù)，特別是在地震成像中[28]。為了求解巖石物理反演問題，本文也采用了基于泰勒展開的線性化巖石物理反演方法[22]。

利用泰勒展開對巖石物理反演問題進(jìn)行線性化，這里使用一階泰勒展開，那么（1）式可以改寫為：

對（2）式中的多項(xiàng)式進(jìn)行展開和合并，巖石物理反演問題可以改寫為線性形式：

再把這種線性化巖石物理反演方法推廣到多個(gè)參數(shù)，那么彈性參數(shù)d和物性參數(shù)m為包含多參數(shù)變量的矢量形式。例如，d為縱波速度vp、橫波速度vs和密度ρ等彈性參數(shù)組成的向量，m為孔隙度φ、含水飽和度Sw和泥質(zhì)含量C等物性參數(shù)組成的向量。同樣巖石物理模型f也寫成矢量形式，那么巖石物理反演問題可以寫為：

對公式（4）進(jìn)行一階泰勒展開，變?yōu)?/p>

其中，Jm0為函數(shù)f在自變量m0處的Jacobian矩陣。再將其改寫成線性化形式則變?yōu)椋?/p>

1.2 巖石物理模型建模

本文選用 Gassmann方程[25]和臨界孔隙度模型[26]對含泥質(zhì)砂巖進(jìn)行巖石物理建模，具體建模過程如下：

①巖石基質(zhì)的體積模量Km、剪切模量μm和密度ρm可以利用Voigt-Reuss-Hill平均方法[29]計(jì)算：

②飽和流體的體積模量Kfl和密度ρfl計(jì)算公式為：

這里利用流體均勻混合公式[9]來計(jì)算飽和流體的體積模量。

③巖石骨架的體積模量Kdry和剪切模量μdry可以利用臨界孔隙度模型[26]計(jì)算：

這里臨界孔隙度值φc雖然是經(jīng)驗(yàn)值，譬如砂巖取40%、灰?guī)r取 60%等，但是臨界孔隙度值包含了巖石的孔隙類型或者孔隙形狀對縱橫波速度的影響[30]。

④飽和巖石的體積模量Ksat和剪切模量μsat可以用Gassmann方程[25]計(jì)算：

⑤飽和巖石的縱波速度vP、橫波速度vs和密度ρ計(jì)算公式為：

由于上述巖石物理建模過程即（7）—（18）式是非線性的，所以可以對其進(jìn)行線性化處理。將彈性參數(shù)d（即縱波速度vP、橫波速度vS和密度ρ）分別對物性參數(shù)m（孔隙度φ、泥質(zhì)含量C和含水飽和度Sw）求導(dǎo)，計(jì)算其對應(yīng)的Jacobian矩陣：

上式中，縱波速度vP對孔隙度φ的偏導(dǎo)數(shù)為：

橫波速度vs對孔隙度φ的偏導(dǎo)數(shù)為：

密度ρ對孔隙度φ的偏導(dǎo)數(shù)為：

縱波速度vP對含水飽和度Sw的偏導(dǎo)數(shù)為：

橫波速度vs對含水飽和度Sw的偏導(dǎo)數(shù)為：

密度ρ對含水飽和度Sw的偏導(dǎo)數(shù)為：

縱波速度vP、橫波速度vs和密度ρ對泥質(zhì)含量C的偏導(dǎo)數(shù)太復(fù)雜，故本文沒有給出具體表達(dá)式。

再將 Jacobian矩陣代入到線性化巖石物理反演問題（6）式中。求解線性化之后的巖石物理反演問題可以有很多方法，例如最小二乘算法，阻尼最小二乘算法或者貝葉斯方法。這里選用阻尼最小二乘方法來求解，可得：

1.3 線性化巖石物理模型精度分析

為了檢驗(yàn)線性化巖石物理反演的精度，分別利用實(shí)際巖石物理模型以及線性化巖石物理模型進(jìn)行含泥質(zhì)砂巖的巖石物理正演模擬分析。假設(shè)含泥質(zhì)砂巖的礦物成分為石英和黏土，孔隙流體為水和氣混合，相應(yīng)的彈性模量和密度分別為：Kq= 3 6GPa ，μq= 3 6GPa ，ρq= 2 .65g/cm3，Kc= 2 1GPa ，μc=15GPa，ρc=2.45 g/cm3，Khc=0.0208GPa，ρhc= 0 .001g/cm3，Kw=2.25 GPa，ρw= 1 .03g/cm3。線性化巖石物理模型均在孔隙度、泥質(zhì)含量和飽和度的均值處進(jìn)行一階泰勒展開。

圖1分別展示了彈性參數(shù)（縱波速度、橫波速度及密度）與孔隙度之間的變化關(guān)系，圖中含泥質(zhì)砂巖的泥質(zhì)含量C以10%為間隔從0到50%變化，含水飽和度Sw為常數(shù)30%，孔隙度φ從0到30%均勻變化。由圖1a和1b可以看出，線性化巖石物理模型與實(shí)際巖石物理模型在孔隙度為 15%左右處吻合最好，在孔隙度較低和較高時(shí)有一定的誤差，這是因?yàn)樵诳紫抖染?5%處進(jìn)行的一階泰勒展開。由圖1c可以看出，線性化巖石物理模型與實(shí)際巖石物理模型的密度完全一樣，這是因?yàn)槊芏扔?jì)算公式本身就是線性的，線性化巖石物理模型與實(shí)際巖石物理模型是完全相同的。

圖1 實(shí)際巖石物理模型及線性化巖石物理模型計(jì)算的彈性參數(shù)與孔隙度之間的變化關(guān)系

圖2分別展示了彈性參數(shù)（縱波速度、橫波速度及密度）與泥質(zhì)含量之間的變化關(guān)系。圖中含泥質(zhì)砂巖的孔隙度φ從以10%間隔從0到30%變化，含水飽和度Sw為常數(shù)50%，泥質(zhì)含量C從0到50%均勻變化。由圖2a和 2b可以看出，線性化巖石物理模型與實(shí)際巖石物理模型在 25%處吻合最好，在泥質(zhì)含量較低和較高時(shí)有一定的誤差，這是因?yàn)樵谀噘|(zhì)含量均值處（25%）進(jìn)行的一階泰勒展開。由圖2c同樣可以看出，線性化巖石物理模型與實(shí)際巖石物理模型的密度是完全一樣的，這是因?yàn)槊芏扔?jì)算公式本身就是線性的，線性化巖石物理模型與實(shí)際巖石物理模型是完全相同的。

圖3分別展示了彈性參數(shù)（縱波速度、橫波速度及密度）與含水飽和度之間的變化關(guān)系。圖中含泥質(zhì)砂巖的孔隙度φ以10%間隔從0到30%變化，泥質(zhì)含量C為常數(shù)25%，含水飽和度Sw從0到100%均勻變化。由圖3a和3b可以看到，利用均勻飽和方式計(jì)算飽和流體的體積模量，實(shí)際巖石物理模型和線性化巖石物理模型計(jì)算的縱波速度吻合很好，僅在含水飽和度Sw接近 100%左右誤差較大，是符合近似要求的。由圖3c同樣可以看出，線性化巖石物理模型與實(shí)際巖石物理模型的密度是完全一樣的，這是因?yàn)槊芏扔?jì)算公式本身就是線性的，線性化巖石物理模型與實(shí)際巖石物理模型是完全相同的。

圖2 實(shí)際巖石物理模型及線性化巖石物理模型計(jì)算的彈性參數(shù)與泥質(zhì)含量之間的變化關(guān)系

圖3 實(shí)際巖石物理模型及線性化巖石物理模型計(jì)算的彈性參數(shù)與含水飽和度之間的變化關(guān)系

2 實(shí)際應(yīng)用對比

2.1 測井?dāng)?shù)據(jù)結(jié)果對比

首先將線性化巖石物理反演方法應(yīng)用于實(shí)際測井?dāng)?shù)據(jù)。A井位于中國東部海上油田，儲(chǔ)集層為含氣砂巖，測井曲線包括孔隙度、泥質(zhì)含量和含水飽和度等物性參數(shù)，以及縱波速度、橫波速度和密度等彈性參數(shù)（見圖4）。測井?dāng)?shù)據(jù)經(jīng)過井震標(biāo)定后由深度域轉(zhuǎn)換成時(shí)間域，并按照2 ms進(jìn)行重采樣。主力產(chǎn)氣層為時(shí)間深度為2 760～2 800 ms的一套厚砂巖儲(chǔ)集層，其上覆2 718 ms左右有一層較薄的泥巖層，2 810～2 822 ms有一套較薄的砂巖儲(chǔ)集層，其下覆為一套較厚泥巖層。

圖4 A井實(shí)際測井?dāng)?shù)據(jù)

在進(jìn)行線性化巖石物理反演之前，首先進(jìn)行正演模擬分析。含泥質(zhì)砂巖的礦物成分為石英和黏土，孔隙流體為水和氣混合，相應(yīng)的彈性模量和密度分別：Kq= 3 8GPa ，μq= 4 0GPa ，ρq= 2 .65g/cm3，Kc=15 GPa，μc= 7 GPa，ρc= 2 .55g/cm3，Khc= 0 .0208GPa ，ρhc= 0 .001g/cm3，Kw=2.25GPa ，ρw= 1 .03g/cm3。

利用精確巖石物理模型以及線性化巖石物理模型，由圖4中的孔隙度、泥質(zhì)含量和含水飽和度等物性參數(shù)來計(jì)算縱波速度、橫波速度和密度，如圖5所示。圖中藍(lán)色實(shí)線代表實(shí)際測井?dāng)?shù)據(jù)，黑色實(shí)線代表實(shí)際巖石物理模型正演計(jì)算結(jié)果，紅色虛線代表線性化巖石物理模型正演計(jì)算結(jié)果。由圖5可見，線性化巖石物理模型與實(shí)際巖石物理模型正演計(jì)算結(jié)果基本一致，與實(shí)際測井曲線有一定誤差，但整體變化趨勢相同。線性化巖石物理模型與實(shí)際彈性參數(shù)之間的平均誤差，縱波速度約為4.0%，橫波速度約為3.0%，密度約為1.2%。在2 820～2 860 ms泥巖段誤差較大，主要是由于泥巖的彈性模量變化較大，很難獲取精確值。

圖5 實(shí)際巖石物理模型與線性化巖石物理模型正演計(jì)算的彈性參數(shù)以及實(shí)際測井彈性參數(shù)

其次，再對測井曲線進(jìn)行線性化巖石物理反演，即利用圖4d—4f中縱波速度、橫波速度和密度來反演圖4a—4c中的孔隙度、泥質(zhì)含量和含水飽和度。使用的彈性模量和密度等參數(shù)與圖5中一致，反演結(jié)果如圖6所示。由圖6可見，線性化巖石物理模型反演結(jié)果與實(shí)際測井曲線基本一致，雖然幅值沒有完全吻合，但整體變化趨勢是相同的。這里給出了線性化巖石物理模型反演結(jié)果與實(shí)際測井曲線之間的相關(guān)系數(shù)和平均誤差，其中孔隙度的相關(guān)系數(shù)為64.8%，平均誤差約為28%，含水飽和度的相關(guān)系數(shù)為85.4%，平均誤差約為17%，泥質(zhì)含量的相關(guān)系數(shù)為90.2%，平均誤差約為10%。

圖6 線性化巖石物理模型反演與實(shí)際測井物性參數(shù)

2.2 地震數(shù)據(jù)結(jié)果對比

在進(jìn)行測井曲線的巖石物理反演之后，再對地震數(shù)據(jù)進(jìn)行線性化巖石物理反演。選擇一條過A井的二維地震剖面，由疊前同時(shí)反演[31]分別獲得了縱波速度、橫波速度和密度等彈性參數(shù)的二維剖面（見圖7）。由二維彈性參數(shù)剖面上可以看到，無論是縱波速度、橫波速度還是密度在2 800 ms左右都有一套較厚的低值區(qū)域，這是主要的目的層段。再利用圖7中的彈性參數(shù)進(jìn)行物性參數(shù)反演（見圖8），圖8a—8c分別展示了線性化巖石物理反演得到的孔隙度、泥質(zhì)含量和含氣飽和度，圖8d—8f分別展示了利用彈性參數(shù)與物性參數(shù)之間的經(jīng)驗(yàn)關(guān)系擬合得到的孔隙度、泥質(zhì)含量和含氣飽和度。由反演結(jié)果可見，兩種反演結(jié)果得到的主要目的層段物性參數(shù)與A井上的物性參數(shù)吻合較好，并且能夠獲取主要目的層段的物性參數(shù)空間展布情況，特別是線性化巖石物理反演含氣飽和度效果要優(yōu)于經(jīng)驗(yàn)關(guān)系擬合結(jié)果。

圖7 地震彈性參數(shù)反演剖面

3 討論

圖8 線性化巖石物理反演的物性參數(shù)

本文提出了一種線性化巖石物理反演方法，由于巖石物理模型一般并不是線性的，所以本文方法較適用于非線性化不是特別強(qiáng)的巖石物理模型，尤其適用于孔隙結(jié)構(gòu)較簡單的儲(chǔ)集層，如碎屑巖儲(chǔ)集層；對于孔隙結(jié)構(gòu)復(fù)雜的儲(chǔ)集層，如碳酸鹽巖儲(chǔ)集層并不太適用。另外值得注意的是，不同的巖石類型應(yīng)該選用不同的巖石物理模型。本文是針對含泥質(zhì)砂巖儲(chǔ)集層進(jìn)行建模，選擇了臨界孔隙度模型和Gassmann方程是合理的，如果針對其他類型儲(chǔ)集層則需要選用其他的巖石物理模型。文中給出的巖石物理正演模型表明，實(shí)際巖石物理模型與線性化巖石物理模型是可以近似的。利用線性化巖石物理模型的優(yōu)勢就是可以獲取巖石物理反演問題的解析解。本文中選擇阻尼最小二乘方法來求解線性反演問題，求解速度較快，但對巖石物理建模中使用的彈性模量和密度等參數(shù)較為敏感。另外，不同于其他巖石物理反演方法，線性化巖石物理反演方法獲得的泥質(zhì)含量和飽和度反演精度要優(yōu)于孔隙度的反演精度。主要是由于在 Jacobian矩陣的逆矩陣中與泥質(zhì)含量以及飽和度有關(guān)的系數(shù)項(xiàng)較大，而與孔隙度有關(guān)的系數(shù)項(xiàng)較小，所以導(dǎo)致泥質(zhì)含量與飽和度的反演效果要優(yōu)于孔隙度反演效果。當(dāng)然，本文中的彈性參數(shù)是由AVO疊前地震反演獲得，下一步還可以直接利用地震數(shù)據(jù)進(jìn)行巖石物理反演。

4 結(jié)論

提出了一種線性化巖石物理反演方法，即通過對實(shí)際巖石物理模型進(jìn)行泰勒級數(shù)展開，得到一階近似的線性化巖石物理模型，再利用阻尼最小二乘方法求解線性化后的巖石物理反演問題，獲得精確的解析解。利用了 Gassmann方程和臨界孔隙度模型對含泥質(zhì)砂巖進(jìn)行實(shí)際巖石物理建模，并且給出了其線性化巖石物理模型表達(dá)式。本方法適用于線性或輕微非線性的巖石物理模型，對于高度非線性巖石物理模型可能不適應(yīng)。實(shí)際測井?dāng)?shù)據(jù)和地震數(shù)據(jù)應(yīng)用表明，線性化巖石物理反演方法能夠反演得到比較好的物性參數(shù)結(jié)果。

符號注釋：

C——泥質(zhì)含量，%；C0——泥質(zhì)含量均值，%；d——彈性參數(shù)；d——彈性參數(shù)矢量形式；f——巖石物理模型；f′——巖石物理模型f的一階導(dǎo)數(shù)；f——巖石物理模型矢量形式；f（m0）——巖石物理模型在m0（即φ0、C0和Sw0）處的值；I——單位矩陣；Jm0——Jacobian矩陣在m0（即φ0、C0和Sw0）處的值；Kc——泥質(zhì)礦物的體積模量，GPa；Kdry——巖石骨架的體積模量，GPa；Kfl——飽和流體的體積模量，GPa；Km——巖石基質(zhì)的體積模量，GPa；Khc——油或氣等碳?xì)浠衔锏捏w積模量，GPa；Kq——砂質(zhì)礦物的體積模量，GPa；Ksat——飽和巖石的體積模量，GPa；Kw——水的體積模量，GPa；m——待預(yù)測的物性參數(shù)；m0——物性參數(shù)m某一具體值，一般選擇m0為待預(yù)測模型參數(shù)的均值；m——待預(yù)測的物性參數(shù)矢量形式；m0——物性參數(shù)矢量形式m某一具體值，即φ0、C0和Sw0；Sw——含水飽和度，%；Sw0——飽和度均值，%；vP——飽和巖石的縱波速度，km/s；vs——飽和巖石的橫波速度，km/s；ρ——飽和巖石的密度，g/cm3；ρc——礦物的密度，g/cm3；ρfl——飽和流體的密度，g/cm3；ρhc——油或氣等碳?xì)浠衔锏拿芏?，g/cm3；ρm——巖石基質(zhì)的密度，g/cm3；ρq——砂質(zhì)礦物的密度，g/cm3；ρw——水的密度，g/cm3；φ——孔隙度，%；φ0——孔隙度均值，%；φc——巖石的臨界孔隙度，%；μc——泥質(zhì)礦物的剪切模量，GPa；μdry——巖石骨架的剪切模量，GPa；μm——巖石基質(zhì)的剪切模量，GPa；μq——砂質(zhì)礦物的剪切模量，GPa；μsat——飽和巖石的剪切模量，GPa；ε——阻尼因子。