一種修正的Kuster-Toks?z巖石物理模型及應(yīng)用

劉致水 孫贊東 董 寧 劉俊州 劉致秀 王亞靜

(①長(zhǎng)安大學(xué)地質(zhì)工程與測(cè)繪學(xué)院，陜西西安 710054； ②中國(guó)石油大學(xué)(北京)，北京 102249； ③中國(guó)石化石油勘探開(kāi)發(fā)研究院，北京100083； ④中國(guó)石油東方地球物理公司，河北涿州 072751； ⑤慶陽(yáng)職業(yè)技術(shù)學(xué)院，甘肅慶陽(yáng) 745000)

1 引言

巖石物理模型是描述彈性參數(shù)與儲(chǔ)層特征之間關(guān)系的數(shù)學(xué)函數(shù)[1-3]。針對(duì)不同的地質(zhì)目標(biāo)和問(wèn)題，前人提出了多種巖石物理理論和模型，包括基礎(chǔ)理論[2-7]，針對(duì)碎屑巖[8-12]、碳酸鹽巖[13-18]、頁(yè)巖[19-22]等的應(yīng)用模型；其中Kuster-Toks?z(KT)模型[4]、微分等效介質(zhì)模型[2，3](DEM)比較經(jīng)典，應(yīng)用廣泛。KT模型能夠同時(shí)考慮多種具有不同形狀的包含物對(duì)彈性參數(shù)的影響，其中不同包含物之間的關(guān)系平等；它的缺點(diǎn)是認(rèn)為巖石中的包含物是“稀疏”存在的，也就是包含物含量不能太高[6-8，14，15]。DEM模型通過(guò)將包含物微分并多次加入巖石基質(zhì)中，從而不受“稀疏”含量包含物的限制。但是DEM模型只能考慮單一包含物類型，在處理多種類型包含物巖石的問(wèn)題時(shí)，只能多次重復(fù)使用DEM公式[14]，且計(jì)算結(jié)果依賴于包含物的添加順序[1]。

針對(duì)上述問(wèn)題，Kuster等[4]、Xu等[8]分別提出用微分等效的思想通過(guò)“迭代”的方式實(shí)現(xiàn)模擬，但是他們并未給出相應(yīng)的公式。尤其需要注意的是，這種模擬不是簡(jiǎn)單地用KT公式迭代就可以得到真實(shí)值，在迭代過(guò)程中由于物質(zhì)的替換，參與計(jì)算的包含物含量和真實(shí)的包含物含量數(shù)值是不同的，如果僅僅用KT公式迭代計(jì)算，則得到的結(jié)果必然是錯(cuò)誤的。李宏兵等[6]基于KT模型推導(dǎo)了一種微分格式的微分等效多孔隙巖石物理模型，進(jìn)而得到了相應(yīng)的解析公式，但是其微分格式與DEM格式一致，計(jì)算過(guò)程中需要求解微分方程；且解析公式是針對(duì)空孔隙的，不能應(yīng)用于包含物體積模量和剪切模量不為零的情況，即不能用于計(jì)算固體包含物的情況，這極大地縮小了模型的應(yīng)用范圍，例如該模型不能將有機(jī)質(zhì)顆粒(干酪根)作為包含物考慮，從而使其很難應(yīng)用于富有機(jī)質(zhì)巖石。在富有機(jī)質(zhì)巖石中，干酪根的體積模量和剪切模量都較小[1,20]。一些學(xué)者[19-22]在構(gòu)建富有機(jī)質(zhì)巖石的巖石物理模型時(shí)將干酪根等效為固體顆粒，使用包含物模型表征干酪根對(duì)巖石速度的影響，這就需要能夠同時(shí)考慮多種不同形狀的固體和流體包含物、又不受“稀疏”含量包含物的限制的巖石物理模型。針對(duì)速度預(yù)測(cè)，Xu等[8]在KT模型基礎(chǔ)上，將巖石中的孔隙分為泥質(zhì)孔隙和砂質(zhì)孔隙，提高了速度預(yù)測(cè)的精度，但是假設(shè)泥質(zhì)孔隙、砂質(zhì)孔隙與巖石中的泥巖和砂巖含量成正比，這種假設(shè)在大多數(shù)情況下是沒(méi)有根據(jù)的。Xu等[15]針對(duì)碳酸鹽巖，將巖石中的孔隙細(xì)分為泥質(zhì)孔隙、印模孔、粒間孔、微裂隙，使用這種復(fù)雜的孔隙類型進(jìn)行速度預(yù)測(cè)得到較好的效果，但是每種孔隙在孔隙度中的體積分?jǐn)?shù)很難準(zhǔn)確求取，使該方法缺乏嚴(yán)謹(jǐn)?shù)臄?shù)學(xué)基礎(chǔ)。

本文將微分等效思維引入KT巖石物理模型，得到一個(gè)適用于高濃度包含物的修正Kuster-Toks?z巖石物理模型，并給出了計(jì)算過(guò)程中所使用的“計(jì)算包含物體積含量”與巖石中真實(shí)包含物體積含量之間的數(shù)學(xué)關(guān)系。利用所建立的新巖石物理模型，構(gòu)建了考慮干酪根與孔隙兩種包含物的富有機(jī)質(zhì)巖石速度預(yù)測(cè)流程，其中干酪根分為球形顆粒與硬幣狀顆粒，兩種顆粒之比固定；包含流體孔隙分為球形孔隙與硬幣狀孔隙，計(jì)算過(guò)程中兩種孔隙類型的體積比可變(自適應(yīng))，在縱波速度的約束下可求得孔隙類型的含量，最后通過(guò)孔隙類型體積比求取橫波速度。

2 基本原理

(1)

(2)

(3)

(4)

硬幣狀裂縫孔隙或包含物的形狀因子為

(5)

(6)

在KT模型中，假設(shè)包含物隨機(jī)分布，以使其等效為各向同性。通過(guò)把包含物模量設(shè)為零模擬干巖石，把包含物的剪切模量設(shè)為零模擬流體飽和巖石。

KT模型有兩個(gè)限制條件： 其一為模擬的是高頻條件下飽和巖石的屬性，適用于超聲實(shí)驗(yàn)室條件； 其二為要求包含物是“稀疏”的，也就是包含物體積含量不能太大。對(duì)于第一個(gè)限制，公認(rèn)且有效的做法是通過(guò)引入Gassmann方程[23]以得到低頻條件下的巖石彈性模量[14]； 對(duì)于第二個(gè)限制，本文引入微分等效的思想，其基本過(guò)程(圖1)為： 將包含物劃分為N小份(極小量)，分為N次加入巖石基質(zhì)中，每一小量代表不同的成分與形狀，每次加入后包含物相替代同體積分?jǐn)?shù)的基質(zhì)相，從而形成混合相，將這個(gè)混合相作為新的基質(zhì)相以進(jìn)行下一次迭代，重復(fù)N次直到加入的包含物相體積分?jǐn)?shù)達(dá)到真實(shí)值。

式(1)和式(2)可以被重寫為

(7)

(8)

圖1 多種包含物微分等效介質(zhì)巖石物理模型的等效過(guò)程示意圖

(9)

(10)

進(jìn)而可寫為

(11)

(12)

每次迭代都可以看作使用小量的包含物替代巖石背景基質(zhì)，迭代完成后，包含物隨機(jī)分散在巖石基質(zhì)中。除第一次迭代外，新加入的包含物所替代的不只是巖石基質(zhì)，也有包含物的一部分，而本質(zhì)上新加入包含物需要替代的物質(zhì)只能是巖石基質(zhì)。因此，每次加入基質(zhì)中的包含物y不是真實(shí)的包含物體積含量分量(dVt)，僅有(1-Vt)y對(duì)迭代過(guò)程中的真實(shí)包含物體積含量分量dVt有貢獻(xiàn)，即每次迭代過(guò)程中的真實(shí)包含物體積含量的變化量為

dVt=(1-Vt)y

(13)

那么，巖石中的真實(shí)包含物體積含量與y之間的關(guān)系為

Vt=1-e-∑y

(14)

式中∑y稱為“計(jì)算包含物體積含量”。如果迭代次數(shù)為N，那么上式為

Vt=1-e-Ny

(15)

此時(shí)，“計(jì)算包含物體積含量”為Ny，y可寫為

(16)

圖2為真實(shí)包含物體積含量(Vt)與“計(jì)算包含物體積含量(Ny)”之間的關(guān)系曲線。由圖可見(jiàn)，當(dāng)包含物含量小于0.1時(shí)，兩者差距較小；隨著真實(shí)包含物體積含量的增大，所需要的計(jì)算包含物體積含量數(shù)值越大。

圖2 真實(shí)包含物與計(jì)算包含物體積含量關(guān)系曲線

表1 不同孔隙度與孔隙形狀條件下體積模 量計(jì)算精度達(dá)到0.001所需迭代次數(shù)

表2 不同孔隙度與孔隙形狀條件下剪切模 量計(jì)算精度達(dá)到0.001所需迭代次數(shù)

3 模型對(duì)比分析

使用兩組數(shù)值進(jìn)行計(jì)算，對(duì)比新模型與KT模型之間的差別。

第一組：往砂巖基質(zhì)孔隙中加入氣體，其中砂巖基質(zhì)體積模量為37GPa，剪切模量為44GPa；孔隙流體的體積模量Kf為0.01GPa；分別對(duì)球形孔、硬幣狀裂縫孔隙進(jìn)行模擬，硬幣狀裂縫的孔隙縱橫比分別設(shè)為0.01、0.05、0.70、1.00，此時(shí)式(9)～式(12)中M=1，式(16)中真實(shí)包含物體積含量為孔隙度φ； 數(shù)值模擬中假設(shè)孔隙度為0～100%，需要注意的是，真實(shí)巖石中的孔隙度不能大于臨界孔隙度[5]。數(shù)值模擬結(jié)果如圖3所示：當(dāng)孔隙度較小時(shí)，新模型與KT模型一致，例如針對(duì)孔隙縱橫比為0.01的硬幣狀裂縫孔隙，在孔隙度小于5%時(shí)兩者基本一致；隨著孔隙度增大，新模型與KT模型的差距變大，與KT模型相比，新模型在高孔隙度時(shí)可以得到更合理的值。一種極端情況是當(dāng)孔隙度為100%時(shí)，真實(shí)值應(yīng)該是孔隙包含物的參數(shù)，新模型可以得到真實(shí)值，而KT模型所計(jì)算的結(jié)果不等于包含物的參數(shù)，這顯然是不對(duì)的。

第二組：往頁(yè)巖基質(zhì)中加入干酪根顆粒，頁(yè)巖基質(zhì)體積模量為39.54GPa，剪切模量為25.68GPa;干酪根顆粒體積模量為2.9GPa，剪切模量為2.7GPa[1]；分別對(duì)球形顆粒、硬幣狀顆粒進(jìn)行模擬，硬幣狀顆粒的縱橫比分別為0.01、0.05、0.70、1.00，此時(shí)式(9)～式(12)中M=1，式(16)中Vt是干酪根體積含量；數(shù)值模擬中假設(shè)干酪根體積含量為0～100%。數(shù)值模擬結(jié)果如圖4所示，其結(jié)果與圖3一致。

圖3 不同孔隙形狀時(shí)KT模型與新模型正演結(jié)果對(duì)比 (a)體積模量； (b)剪切模量

圖4 不同干酪根顆粒形狀時(shí)KT模型與新模型正演結(jié)果對(duì)比 (a)體積模量； (b)剪切模量

4 模型在富有機(jī)質(zhì)巖石中的應(yīng)用

富有機(jī)質(zhì)巖石的成分包括礦物、干酪根、飽和流體孔隙。干酪根的物理參數(shù)較為特殊，如表3所示：干酪根的體積模量、密度與石英、方解石、黏土等礦物相差較大，與水接近；但是它具有剪切模量，又與流體不同。本文將干酪根等效為巖石中的一種包含物，與飽和流體孔隙一同加入巖石基質(zhì)中，從而進(jìn)行富有機(jī)質(zhì)巖石的速度預(yù)測(cè)。

利用Vernik等[24]在超聲實(shí)驗(yàn)室針對(duì)干巖石測(cè)量的數(shù)據(jù)(70MPa圍壓環(huán)境)說(shuō)明干酪根、孔隙對(duì)富有機(jī)質(zhì)巖石速度的影響并測(cè)試新模型。本文使用了其測(cè)量的59個(gè)干巖石樣品數(shù)據(jù)中的51個(gè)，未使用縱、橫波速度較高的8個(gè)樣品。所使用的51個(gè)樣品的干酪根含量為1.4%～42.3%，孔隙度為0.6%～30.9%。圖5是樣品的縱波速度、橫波速度與孔隙度、干酪根含量交會(huì)圖。由圖可見(jiàn)：在干酪根含量近似不變時(shí)(例如深藍(lán)色時(shí))，隨著孔隙度的增加，巖石的縱、橫波速度整體呈減小趨勢(shì)；在孔隙度近似不變時(shí)(例如孔隙度接近0)，可非常明顯地看到隨干酪根含量的增加，巖石的縱、橫波速度減小。此外，縱、橫波速度數(shù)據(jù)點(diǎn)較分散，例如圖中用直線段標(biāo)注的五個(gè)數(shù)據(jù)點(diǎn)，其孔隙度近10.0%，而縱波速度在3.2～5.1km/s范圍內(nèi)變化，橫波速度在2.1km/s～3.1km/s范圍內(nèi)變化，本文認(rèn)為這種速度分散是由孔隙形狀及干酪根顆粒形狀的變化所引起。

表3 巖石組成成分的彈性參數(shù)和密度[1]

Xu等[8]將孔隙分為柔性孔隙與剛性孔隙，使用Wu[25]的橢球體模型表征孔隙類型(孔隙縱橫比的大小表征孔隙的形狀)，其中剛性孔隙的孔隙縱橫比為0.12～0.15，柔性孔隙的孔隙縱橫比為0.02～0.03。剛性孔隙與柔性孔隙的體積比與巖石中砂巖/泥巖體積比一致。本文將頁(yè)巖中的干酪根和孔隙等效為具有一定形狀的包含物，其中干酪根體積含量(Vk)等效為由剛性顆粒(Vks)與柔性顆粒(Vkc)組成，將孔隙度(φ)等效為由剛性孔隙度φs與柔性孔隙度φc組成，其中剛性顆粒和孔隙使用多種孔隙理論[1]中的球型包含物表征，柔性顆粒和孔隙使用硬幣形狀包含物表征(顆粒及孔隙縱橫比取0.02)。理論上通過(guò)新模型可以考慮多種包含物及孔隙類型，但是難點(diǎn)為需要知道每種類型包含物各自的體積分?jǐn)?shù)。

圖5 富有機(jī)質(zhì)巖石的縱波速度(a)、橫波速度(b)與孔隙度、干酪根含量的交會(huì)圖(數(shù)據(jù)來(lái)自文獻(xiàn)[24])

由圖5可見(jiàn)，干酪根含量較高(暖色表征)時(shí)，縱波與橫波速度的分散性較低，因此假定速度分散主要由孔隙形狀的變化引起，在計(jì)算中固定剛性干酪根顆粒Vks與柔性干酪根顆粒Vkc的比例，而剛性孔隙與柔性孔隙的體積比是可變(自適應(yīng))的。經(jīng)過(guò)多次試算并求取預(yù)測(cè)值相對(duì)于實(shí)際值的均方根誤差(RMSE)，本文最終選擇剛性干酪根顆粒與柔性干酪根顆粒的比例為1∶1。圖6為縱、橫波的均方根誤差隨剛性干酪根顆粒在干酪根總量中的占比Vks的變化(間隔為-0.05)。由圖6可見(jiàn)，在1～0之間，縱波的RMSE有兩個(gè)極小值：當(dāng)Vks=0.75時(shí)縱波的均方根誤差為0.0606km/s，其所對(duì)應(yīng)的橫波均方根誤差為0.0812km/s；當(dāng)Vks=0.5時(shí)縱波的均方根誤差為0.0607km/s，其所對(duì)應(yīng)的橫波均方根誤差為0.0772km/s。因此本文選擇Vks=0.5(即剛性干酪根顆粒與柔性干酪根顆粒相等)進(jìn)行計(jì)算。

圖6 縱、橫波速度的均方根誤差隨剛性 干酪根顆粒的體積比Vks的變化

此時(shí)，式(15)中Vt指孔隙度φ與干酪根體積含量Vk之和；式(9)～式(12)中M=4，包含物體積分?jǐn)?shù)Vi有四個(gè)：球形干酪根顆粒V1、硬幣狀干酪根顆粒V2、球形孔隙V3、硬幣狀孔隙V4，其具體表達(dá)式為

(17)

(18)

(19)

(20)

式中：Vks+Vkc=1，本文Vks、Vkc取值均為0.5；φc+φs=1，φc與φs體積比可變。

利用式(11)、式(12)、式(16)及式(17)～式(20)，將干酪根顆粒和孔隙(含流體孔隙或干孔隙)加入巖石基質(zhì)中，得到體積模量K和剪切模量G，則巖石的縱波速度和橫波速度可以表示為

(21)

(22)

式中ρ=ρm(1-φ-Vk)+ρfφ+ρkVk，ρm、ρk、ρf分別指巖石基質(zhì)、干酪根、孔隙中流體的密度，在干巖石中ρf取為零。需要說(shuō)明的是，通過(guò)上述過(guò)程得到的是超聲波條件下的高頻速度，要得到低頻條件下的飽和流體巖石速度，則須利用式(11)～式(16)以及式(17)～式(20)獲得干燥空腔的等效模量，再利用Gassmann理論往空腔中加入流體。

圖7為利用上述富有機(jī)質(zhì)巖石速度計(jì)算流程正演結(jié)果與實(shí)測(cè)數(shù)據(jù)的對(duì)比。由圖可知：當(dāng)孔隙度不變、干酪根含量增加時(shí)，巖石的速度降低；在干酪根含量不變、孔隙度增加時(shí)，巖石的速度降低；隨著球形孔隙含量的降低，巖石的速度降低。通過(guò)孔隙度、干酪根含量的變化，可以覆蓋所有數(shù)據(jù)點(diǎn)，即新模型可以解釋測(cè)量數(shù)據(jù)的孔隙度—干酪根含量—速度之間的關(guān)系。

這樣，就建立了縱、橫波速度與可變的孔隙參數(shù)(φc、φs)之間的關(guān)系，通過(guò)使計(jì)算的縱波速度vPcal與實(shí)測(cè)的縱波速度vPmea之差達(dá)到最小來(lái)建立目標(biāo)函數(shù)

ε=|vPcal-vPmea|

(23)

采用迭代算法使該式取極小值以反演得到孔隙參數(shù)(φc和φs)。計(jì)算過(guò)程中需要輸入的數(shù)據(jù)包括干酪根體積含量Vk、孔隙度φ、基質(zhì)體積含量1-Vk-φ； 輸入的參數(shù)還包括： 頁(yè)巖基質(zhì)體積模量為39.54GPa、剪切模量為25.68GPa、密度為2.64g/cm3，干酪根顆粒體積模量為2.9GPa、剪切模量為2.7GPa、密度為1.3g/cm3，孔隙中氣體的體積模量為0.18GPa、密度為0.26g/cm3。包含物總數(shù)M設(shè)為4，迭代次數(shù)N設(shè)為200。最終，將這些參數(shù)代入式(9)～式(16)以計(jì)算橫波速度。計(jì)算過(guò)程中取ε=0.005為約束量。圖8為預(yù)測(cè)縱波與實(shí)測(cè)縱波交會(huì)圖，圖中位于對(duì)角線上的數(shù)據(jù)點(diǎn)代表預(yù)測(cè)速度與實(shí)測(cè)速度之間吻合度高。由圖可見(jiàn)，新模型與KT模型所得縱波差別較小，絕大部分?jǐn)?shù)據(jù)點(diǎn)實(shí)測(cè)縱波速度與計(jì)算縱波速度吻合度高。有幾個(gè)點(diǎn)的實(shí)測(cè)數(shù)據(jù)與預(yù)測(cè)數(shù)據(jù)不吻合，超出了約束量范圍，其原因是所取的基質(zhì)參數(shù)、干酪根顆粒的形狀與該數(shù)據(jù)點(diǎn)的情況差別較大，孔隙度、孔隙類型、干酪根含量的變化不足以消除這種影響，說(shuō)明在實(shí)際應(yīng)用中，需要根據(jù)具體的研究對(duì)象選擇基質(zhì)參數(shù)以及干酪根顆粒形狀。

圖9顯示了新模型和KT模型所得到的球形孔隙和硬幣狀孔隙在包含物中所占的體積分?jǐn)?shù)。由圖可見(jiàn)：①KT模型所得結(jié)果中球形孔隙體積分?jǐn)?shù)或硬幣狀孔隙體積分?jǐn)?shù)為零的情況多于新模型，說(shuō)明KT模型在孔隙變化范圍內(nèi)不能達(dá)到縱波約束量的情況多于新模型，這與圖8結(jié)果一致；②小孔隙巖樣所得的球形孔隙體積分?jǐn)?shù)或硬幣狀孔隙體積分?jǐn)?shù)差別不大，對(duì)于大孔隙度巖樣，新模型所得球形孔隙體積分?jǐn)?shù)小于KT模型結(jié)果，所得硬幣狀孔隙體積分?jǐn)?shù)略大于KT模型結(jié)果，說(shuō)明新模型在剛性孔隙含量較少的情況下得到了準(zhǔn)確的速度預(yù)測(cè)結(jié)果，這與圖3、圖4正演結(jié)果一致。

圖10為新模型和KT模型預(yù)測(cè)橫波與實(shí)測(cè)橫波結(jié)果的對(duì)比。圖中的點(diǎn)離對(duì)角線越近表明預(yù)測(cè)速度與實(shí)測(cè)速度之間的吻合度越高。顯然，新模型所得預(yù)測(cè)速度與實(shí)測(cè)速度吻合度高(圖10a)。

圖7 基于新模型的縱波速度(a)和橫波速度(b)正演結(jié)果

其中紅線表示基質(zhì)中加入5%的干酪根顆粒，縱、橫波速度隨孔隙度(0～40%，干孔隙空腔)的變化，黑線表示基質(zhì)中加入40%的干酪根顆粒，縱、橫波速度隨孔隙度(0～40%)的變化。從上到下的六條紅線和黑線分別代表球型孔隙體積比(φs)為100%、80%、60%、40%、20%、0

圖8 縱波速度約束計(jì)算的結(jié)果與實(shí)測(cè)結(jié)果對(duì)比 (a)新模型； (b)KT模型

圖9 新模型與KT模型計(jì)算的硬幣狀裂縫孔隙與球形孔隙占包含物總量的體積分?jǐn)?shù) (a)新模型硬幣狀孔隙； (b)新模型球形孔隙； (c)KT模型硬幣狀孔隙； (d)KT模型球形孔隙

圖10 預(yù)測(cè)的橫波速度與實(shí)測(cè)橫波速度交會(huì)圖 (a)新模型； (b)KT模型

5 結(jié)論

本文將微分等效思維引入KT模型，提出了一種適用于高濃度包含物的修正Kuster -Toks?z巖石物理模型。該巖石物理模型保留了KT模型能夠同時(shí)考慮多種包含物類型對(duì)彈性參數(shù)的影響的優(yōu)點(diǎn)，又適合于高濃度包含物的情況。將該模型應(yīng)用于富有機(jī)質(zhì)巖石中，將孔隙和干酪根固體顆粒作為包含物同時(shí)加入巖石基質(zhì)，結(jié)果證明新模型比KT模型更適合于包含物(孔隙和干酪根固體)濃度較高的情況；利用建立的巖石物理模型對(duì)實(shí)驗(yàn)室測(cè)試的富有機(jī)質(zhì)巖石資料進(jìn)行橫波速度預(yù)測(cè)，干酪根顆粒采取固定比例的球形與硬幣狀組合，孔隙采取可變的球形與硬幣形狀組合，在縱波速度的約束下預(yù)測(cè)橫波速度，結(jié)果證明，與KT模型相比，新模型的預(yù)測(cè)誤差更小。

[1] Mavko G，Mukerji T and Dvorikin J.The Rock Phy-sics Handbook：Tools for Seismic Analysis in Porous Media.Cambridge University Press，1998.

[2] Norris A N.A differential scheme for the effective moduli of composites.Mechanics of Materials，1985，4(6)：1-16.

[3] Zimmerman R W.Compressibility of Sandstones.Elsevier，New York，1991.

[4] Kuster G and Toks?z M.Velocity and attenuation of seismic waves in two phase media.Geophysics，1974，39(5)：587-618.

[5] Nur A.Critical porosity and the seismic velocities in rocks.Eos Transactions American Geophysical Union，1992，73(1)：43-66.

[6] 李宏兵，張佳佳.多重孔巖石微分等效介質(zhì)模型及其干燥情形下的解析近似式.地球物理學(xué)報(bào)，2014，57(10)：3422-3430. Li Hongbing，Zhang Jiajia.A differential effective medium model of multiple porosity rock and its analytical approximations for dry rock.Chinese Journal of Geophysics，2014，57(10)：3422-3430.

[7] 印興耀，化世榜，宗兆云.基于線性近似的微分等效介質(zhì)方程解耦方法.石油地球物理勘探，2016，51(2)：281-287. Yin Xingyao，Hua Shibang，Zong Zhaoyun.A decoupling approach for differential equivalent equations based on linear approximation.OGP，2016，51(2)：281-287.

[8] Xu S and White R E.A new velocity model for clay-sand mixtures.Geophysical Prospecting，1995，43(1)：91-118.

[9] 劉雅杰，李生杰，王永剛等.橫波預(yù)測(cè)技術(shù)在蘇里格氣田儲(chǔ)層預(yù)測(cè)中的應(yīng)用.石油地球物理勘探，2016，51(1)：165-172. Liu Yajie，Li Shengjie，Wang Yonggang et al.Reservoir prediction based on shear wave in Sulige Gas Field.OGP，2016，51(1)：165-172.

[10] 唐杰，王浩，姚振岸等.基于巖石物理診斷的橫波速度估算方法.石油地球物理勘探，2016，51(3)：537-543. Tang Jie，Wang Hao，Yao Zhen’an et al.Shear wave velocity estimation based on rock physics diagnosis.OGP，2016，51(3)：537-543.

[11] 侯波，陳小宏，張孝珍等.臨界孔隙度Pride模型及其應(yīng)用.石油地球物理勘探，2012，47(2)：277-281. Hou Bo，Chen Xiaohong，Zhang Xiaozhen et al.Critical porosity Pride model and its application.OGP，2012，47(2)：277-281.

[12] 郭棟，印興耀，吳國(guó)忱等.橫波速度計(jì)算方法與應(yīng)用.石油地球物理勘探，2007，42(5)：535-538. Guo Dong，Yin Xingyao，Wu Guochen et al.Computational approach of S-wave velocity and application.OGP，2007，42(5)：535-538.

[13] Sun S Z，Stretch S R and Brown R J.Comparison of borehole velocity-prediction models and estimation of fluid saturation effects：From rock physics to exploration problem.Journal of Canadian Petroleum Techno-logy，2004，43(3)：18-26.

[14] Sun S Z，Wang H Y，Liu Z S et al.The theory and application of DEM-Gassmann rock physics model for complex carbonate reservoirs.The Leading Edge，2012，31(2)：152-158.

[15] Xu S and Payne M.Modeling elastic properties in car-

bonate rocks.The Leading Edge，2009，28(1)：66-74.

[16] Wang H Y，Sun S Z，Yang H J et al.The influence of pore structure on P- and S-wave velocities in complex carbonate reservoirs with secondary storage space.Petroleum Science，2011，8(4)：394-405.

[17] 張秉銘，劉致水，劉俊州等.鄂爾多斯盆地北部復(fù)雜碳酸鹽巖橫波速度預(yù)測(cè)研究.石油物探，2017，56(3)：328-337. Zhang Bingming，Liu Zhishui，Liu Junzhou et al.An improved S-wave velocity prediction method for complex carbonate reservoir in North Ordos Basin，China.GPP，2017，56(3)：328-337.

[18] 張廣智，李呈呈，印興耀等.基于修正Xu-White模型的碳酸鹽巖橫波速度估算方法.石油地球物理勘探，2012，47(5)：717-722. Zhang Guangzhi，Li Chengcheng，Yin Xingyao et al.A shear velocity estimation method for carbonate rocks based on the improved Xu-White model.OGP，2012，47(5)：717-722.

[19] 董寧，霍志周，孫贊東等.泥頁(yè)巖巖石物理建模研究.地球物理學(xué)報(bào)，2014，57(6)：1990-1998. Dong Ning，Huo Zhizhou，Sun Zandong et al.An investigation of a new rock physics model for shale.Chinese Journal of Geophysics，2014，57(6)：1990-1998.

[20] 劉致水，孫贊東.新型脆性因子及其在泥頁(yè)巖儲(chǔ)層預(yù)測(cè)中的應(yīng)用.石油勘探與開(kāi)發(fā)，2015，42(1)：117-124. Liu Zhishui,Sun Zandong.New brittleness indexes and their application in shale/clay gas reservoir prediction.Petroleum Exploration and Development，2015，42(1)：117-124.

[21] 胡起，陳小宏，李景葉.基于單孔隙度縱橫比模型的有機(jī)頁(yè)巖橫波速度預(yù)測(cè)方法.地球物理學(xué)進(jìn)展，2014，29(5)：2388-2394. Hu Qi，Chen Xiaohong，Li Jingye.Shear velocity prediction for organic shales based on the single aspect ratio model.Progress in Geophysics，2014，29(5)：2388-2394.

[22] Guo Z，Chapman M，Li X.A shale rock physics model and its application in the prediction of brittleness index，mineralogy，and porosity of the Barnett Shale.SEG Technical Program Expanded Abstracts，2012，31.

[23] Gassmann F.Elastic waves through a packing of spheres.Geophysics，1951，16(4)：673-685.

[24] Vernik L and Liu X Z.Velocity anisotropy in shales：a petrophysical study.Geophysics，1997，62(2)：521-532.

[25] Wu T T.The effect of inclusion shape on the elastic moduli of a two phase material.International Journal of Solids & Structures，1966，2(1)：1-8.