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

（1.中國石油大學（華東）地球科學與技術學院，山東青島 266580；2.深層油氣地質與勘探教育部重點實驗室，山東青島 266580；3.海洋國家實驗室海洋礦產資源評價與探測技術功能實驗室，山東青島 266071）

0 引言

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

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

1 方法原理

1.1 線性化巖石物理反演

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

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

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

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

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

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

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

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

1.2 巖石物理模型建模

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

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

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

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

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

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

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

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

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

上式中，縱波速度vP對孔隙度φ的偏導數為：

橫波速度vs對孔隙度φ的偏導數為：

密度ρ對孔隙度φ的偏導數為：

縱波速度vP對含水飽和度Sw的偏導數為：

橫波速度vs對含水飽和度Sw的偏導數為：

密度ρ對含水飽和度Sw的偏導數為：

縱波速度vP、橫波速度vs和密度ρ對泥質含量C的偏導數太復雜，故本文沒有給出具體表達式。

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

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

為了檢驗線性化巖石物理反演的精度，分別利用實際巖石物理模型以及線性化巖石物理模型進行含泥質砂巖的巖石物理正演模擬分析。假設含泥質砂巖的礦物成分為石英和黏土，孔隙流體為水和氣混合，相應的彈性模量和密度分別為：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。線性化巖石物理模型均在孔隙度、泥質含量和飽和度的均值處進行一階泰勒展開。

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

圖1 實際巖石物理模型及線性化巖石物理模型計算的彈性參數與孔隙度之間的變化關系

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

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

圖2 實際巖石物理模型及線性化巖石物理模型計算的彈性參數與泥質含量之間的變化關系

圖3 實際巖石物理模型及線性化巖石物理模型計算的彈性參數與含水飽和度之間的變化關系

2 實際應用對比

2.1 測井數據結果對比

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

圖4 A井實際測井數據

在進行線性化巖石物理反演之前，首先進行正演模擬分析。含泥質砂巖的礦物成分為石英和黏土，孔隙流體為水和氣混合，相應的彈性模量和密度分別：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中的孔隙度、泥質含量和含水飽和度等物性參數來計算縱波速度、橫波速度和密度，如圖5所示。圖中藍色實線代表實際測井數據，黑色實線代表實際巖石物理模型正演計算結果，紅色虛線代表線性化巖石物理模型正演計算結果。由圖5可見，線性化巖石物理模型與實際巖石物理模型正演計算結果基本一致，與實際測井曲線有一定誤差，但整體變化趨勢相同。線性化巖石物理模型與實際彈性參數之間的平均誤差，縱波速度約為4.0%，橫波速度約為3.0%，密度約為1.2%。在2 820～2 860 ms泥巖段誤差較大，主要是由于泥巖的彈性模量變化較大，很難獲取精確值。

圖5 實際巖石物理模型與線性化巖石物理模型正演計算的彈性參數以及實際測井彈性參數

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

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

2.2 地震數據結果對比

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

圖7 地震彈性參數反演剖面

3 討論

圖8 線性化巖石物理反演的物性參數

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

4 結論

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

符號注釋：

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