劉欣欣，印興耀，張 峰

(1.中國(guó)石油大學(xué)地球科學(xué)與技術(shù)學(xué)院，山東青島 266580;2.中國(guó)石油大學(xué)地球物理與信息工程學(xué)院，北京 102249)

碳酸鹽巖儲(chǔ)層是重要的油氣勘探目標(biāo)。該類儲(chǔ)層非均質(zhì)性很強(qiáng)，孔隙度、滲透率和儲(chǔ)層其他特性之間的關(guān)系非常復(fù)雜，導(dǎo)致儲(chǔ)層特征描述困難［1］。橫波速度有利于儲(chǔ)層巖性、物性和含油氣性的分析與識(shí)別，在疊前地震反演、AVO屬性分析、測(cè)井巖性解釋以及多波資料解釋等過(guò)程中都有非常重要的作用。由于橫波測(cè)井成本高，解釋難度大等原因，實(shí)際生產(chǎn)中許多井資料缺少橫波測(cè)井信息［2］。基于巖石物理理論的方法是計(jì)算地震橫波速度的重要手段，許多學(xué)者通過(guò)巖石物理試驗(yàn)和分析建立了多種經(jīng)驗(yàn)關(guān)系［3-4］及理論模型［5-7］。目前，國(guó)內(nèi)相關(guān)研究主要側(cè)重于使用測(cè)井或其他資料約束提高橫波速度計(jì)算精度或者改善實(shí)際應(yīng)用效果［8-12］。針對(duì)碳酸鹽巖儲(chǔ)層的特殊性和復(fù)雜性，Kumar等［13-14］研究了孔隙幾何形狀對(duì)橫波速度的影響，提出了估算不同類型孔隙的縱橫比和孔隙度的方法。Vega等［15］研究了非均質(zhì)碳酸鹽巖中流體置換模型的適用性，指出斑塊飽和模型能預(yù)測(cè)更大的速度變化。筆者根據(jù)巖石物理理論和模型計(jì)算飽和巖石的橫波速度，討論碳酸鹽巖中孔隙微結(jié)構(gòu)參數(shù)和礦物組分彈性模量的求取方法，進(jìn)而確定巖石的孔隙結(jié)構(gòu)和巖性參數(shù)，并在此基礎(chǔ)上計(jì)算橫波速度。

1 碳酸鹽巖儲(chǔ)層橫波速度計(jì)算原理

巖石是由固體的巖石骨架和流動(dòng)的孔隙流體組成的雙相體，在波長(zhǎng)遠(yuǎn)遠(yuǎn)大于巖石的不均勻尺度時(shí)，巖石的彈性性質(zhì)等效于雙相體的等效彈性性質(zhì)［16］。根據(jù)巖石物理多孔介質(zhì)理論，針對(duì)碳酸鹽巖儲(chǔ)層巖性和孔隙結(jié)構(gòu)的復(fù)雜性，考慮巖石礦物組分、孔隙微結(jié)構(gòu)、孔隙流體的影響，采用不同的巖石物理理論分別計(jì)算碳酸鹽巖儲(chǔ)層的巖石基質(zhì)、干燥巖石骨架以及飽和巖石的等效彈性模量，進(jìn)而計(jì)算儲(chǔ)層的橫波速度。

1.1 巖石基質(zhì)等效彈性模量計(jì)算

Hashin和Shtrikman(1963)提出不同組分混合彈性體模量的最低上限和最高下限(HS界限)［17］。圖1是用不同理論計(jì)算的方解石和泥質(zhì)組成的巖石基質(zhì)的體積和剪切模量。HS上下限比Voigt和Reuss界限范圍要窄。在沒(méi)有足夠的信息的情況下，近似使用HS上下限的平均來(lái)計(jì)算巖石基質(zhì)的彈性模量。

圖1 不同等效介質(zhì)模型計(jì)算的巖石基質(zhì)體積模量和剪切模量Fig.1 Bulk modulus and shear modulus of rock matrix calculated by different effective medium models

1.2 巖石干燥骨架等效彈性模量計(jì)算

Kuster等［18］假設(shè)孔隙為橢球形，基于散射理論推導(dǎo)出巖石的體積模量K、剪切模量G與孔隙度和孔隙縱橫比的聯(lián)系(KT模型)，令M=(K，G)，T=(P，Q)，則有

式中，φi為第 i種孔隙的體積分?jǐn)?shù);Mi=(Ki，Gi)，Ki和Gi分別為孔隙內(nèi)含物的體積模量和剪切模量;Mm=(Km，Gm)，Km和Gm分別為背景基質(zhì)的體積模量和剪切模量;系數(shù) Tmi=(Pmi，Qmi)，Pmi和 Qmi反映第i種孔隙對(duì)巖石基質(zhì)的影響，是孔隙縱橫比的函數(shù)。

微分有效介質(zhì)理論(DEM模型)是通過(guò)往固體礦物相中逐漸加入包含物相來(lái)模擬雙相混合物。彈性模量的耦合的微分方程［19］為

Berryman［20］基于彈性波散射理論，推導(dǎo)出含橢圓形包含體的巖石等效彈性模量，提出自相容近似模型(SCA模型):

使用以上三種模型計(jì)算干燥巖石骨架的彈性模量，結(jié)果見(jiàn)圖2(孔隙縱橫比為0.5)?？紫抖容^小時(shí)，三種模型的計(jì)算結(jié)果相當(dāng)。由于KT模型要求孔隙度與孔隙縱橫比的比值要遠(yuǎn)小于1，因此適用于計(jì)算較高縱橫比或較小孔隙度巖石的有效彈性模量;DEM模型沒(méi)有考慮孔隙包含物之間散射的影響，而且包含物的加入順序會(huì)影響微分彈性模量的計(jì)算結(jié)果;SCA模型假設(shè)孔隙是相互孤立的，而且波長(zhǎng)遠(yuǎn)大于包含物尺寸，對(duì)孔隙較大的巖石是適用的，但是如果孔隙為針狀、盤狀和硬幣狀，那么SCA模型的使用必須謹(jǐn)慎。對(duì)于縱橫比較大(＞0.1)以及縱橫比較小(＜0.1)的孔隙，分別使用KT模型和DEM模型計(jì)算巖石骨架的彈性模量。

圖2 不同模型計(jì)算的干燥巖石骨架體積模量(左)和剪切模量(右)Fig.2 Bulk modulus and shear modulus of dry rock frame calculated by different effective medium models

1.3 飽和巖石的等效彈性模量的計(jì)算

Gassmann方程假設(shè)巖石宏觀上均勻分布，孔隙之間相互連通，但在碳酸鹽巖地層中，膠結(jié)或者溶蝕作用使得孔隙形狀非常不規(guī)則，或者填隙物的存在堵塞局部流體的流動(dòng)通道，導(dǎo)致孔隙流體部分飽和并且不均勻分布，因此Gaussman方程的適用性受到一定的限制。White給出了斑塊飽和模型的定量計(jì)算［21］。在不同斑塊內(nèi)的流體飽和狀態(tài)不同，當(dāng)巖石滲透率較低或者含水飽和度增大時(shí)，誘發(fā)的孔隙壓力增量不能在完全飽和的斑塊和未完全飽和的斑塊之間達(dá)到平衡，從而使巖石骨架硬化，此時(shí)巖石體積模量增大而剪切模量保持不變［22］。假設(shè)巖石孔隙中填充油、氣、水三相流體的混合物，各相流體飽和度分別為Sf1、Sf2、Sf3，斑塊飽和模型將不均勻分布的各相流體分開處理，首先分別計(jì)算各相流體單獨(dú)分布時(shí)巖石的彈性模量:

使用上述模量的等效平均計(jì)算斑塊飽和巖石的等效彈性模量:

2 關(guān)鍵參數(shù)的求取

2.1 碳酸鹽巖孔隙微觀結(jié)構(gòu)參數(shù)計(jì)算

孔隙形狀是影響碳酸鹽巖石彈性的重要參數(shù)［23-24］。國(guó)內(nèi)外關(guān)于碳酸鹽巖孔隙出現(xiàn)了多種分類方法(Lucia，1995;L?n?y，2006)［25-26］?？紤]到碳酸鹽巖孔隙形狀，本文將孔隙分為孔洞(縱橫比0.5～1.0)、粒間孔隙(縱橫比0.1～0.5)以及裂隙(縱橫比0.01～0.1)?？锥创硖妓猁}巖中堅(jiān)硬的近圓形孔隙或溶洞，裂隙代表巖石中最柔軟的孔隙，受到外力作用時(shí)易閉合，對(duì)巖石速度的影響較大。在含有泥質(zhì)的情況下，考慮孔隙的連通性再細(xì)分出第四種孔隙:泥質(zhì)孔隙，主要代表含束縛水或者不連通孔隙。則總孔隙度φ可以表示為

式中，φvu、φin、φcr、φcl分別表示孔洞、粒間孔隙、裂隙、泥質(zhì)孔隙的孔隙度。假設(shè)四種孔隙的縱橫比分別為 αvu、αin、αcr、αcl。

對(duì)于泥質(zhì)孔隙，根據(jù)Xu-White模型對(duì)其進(jìn)行描述。碳酸鹽巖中只含有粒間或者晶間孔隙(縱橫比約為0.1)時(shí)，根據(jù)DEM理論計(jì)算的縱波速度與Wyllie時(shí)間平均方程的計(jì)算結(jié)果非常接近;對(duì)于鑄?？紫痘蛘呖锥吹葓A形孔隙(縱橫比約為1)，DEM理論計(jì)算的縱波速度接近于HS上邊界;而對(duì)于裂隙(縱橫比大約為0.01)，DEM理論計(jì)算的縱波速度接近于HS下邊界［23］。為此使用下面的方法計(jì)算這三種孔隙的縱橫比:

(1)給定初值條件。假設(shè)孔洞、粒間孔隙、裂隙的孔隙縱橫比α分別為0.8、0.1、0.01;

(2)假設(shè)巖石的全部孔隙為孔洞，分別據(jù)DEM理論、HS上限計(jì)算巖石縱波速度，得兩速度的差ΔvP;

假設(shè)巖石的全部孔隙為粒間孔隙，使用DEM理論、Wyllie時(shí)間平均方程計(jì)算縱波速度，得ΔvP;

假設(shè)巖石的全部孔隙為裂隙，使用DEM理論、HS下限計(jì)算縱波速度，得ΔvP;

將實(shí)測(cè)縱波速度與Wyllie時(shí)間平均方程計(jì)算的縱波速度之差作為縱波速度偏差，該偏差是由于不同的孔隙類型引起的，正偏差代表鑄?？紫兜容^大縱橫比孔隙的影響，零偏差代表粒間孔隙的影響，而負(fù)偏差表示裂縫或者高含氣量的影響［27］。統(tǒng)計(jì)正、零、負(fù)速度偏差三者的比例，作為孔洞、粒間孔隙、裂隙三種孔隙分配的初始比例。根據(jù)實(shí)測(cè)縱波速度ΔvPm，計(jì)算這三種孔隙的孔隙度:

(1)按照初始比例對(duì)除泥質(zhì)孔隙以外的孔隙進(jìn)行分配，得到三種孔隙的初始孔隙度φvu、φin、φcr;

(2)計(jì)算巖石的縱波速度，得ΔvPe，計(jì)算其與實(shí)測(cè)縱波速度ΔvPm的差ΔvP;

若ΔvP＞ε，則φcr=φcr+Δφ，且① 若 φvu＞0，φvu=φvu-Δφ;② 若φvu=0，則φin=φin-Δφ。轉(zhuǎn)步驟(3);

若ΔvP＜ -ε，則 φvu=φvu+Δφ，且① 若 φcr＞0，φcr=φcr-Δφ;② 若φcr=0，則φin=φin-Δφ。轉(zhuǎn)步驟(3)。

使用上述方法可求取四種孔隙度各自的縱橫比和孔隙度，也就得到了碳酸鹽巖孔隙微結(jié)構(gòu)參數(shù)，然后對(duì)于不同縱橫比的孔隙，使用KT模型和DEM模型計(jì)算干燥巖石骨架的等效彈性模量。

2.2 巖石礦物組分彈性模量反演

碳酸鹽巖礦物組成復(fù)雜，包括方解石或白云石，經(jīng)?；旌仙倭康氖?，泥質(zhì)或無(wú)水石膏(硬石膏)，以及其他非碳酸鹽自生礦物和陸源碎屑礦物。本文中將碳酸鹽巖基質(zhì)等效為白云石、方解石和泥質(zhì)等礦物組分的混合物。考慮到隨著深度變化，地層的溫度和壓力、巖石的巖性等都會(huì)變化，礦物組分的模量可能隨之發(fā)生變化，而且在不同地區(qū)礦物模量也不是一成不變的。儲(chǔ)層縱橫波速度在一定程度上反映巖石礦物彈性性質(zhì)的變化，但是縱橫波速度受多種因素的影響，其與礦物組分彈性模量的顯函數(shù)關(guān)系式很難建立。為此，根據(jù)縱橫波速度計(jì)算方法，建立目標(biāo)函數(shù)，通過(guò)使其取最小值反演求取各種礦物組分的體積模量K和剪切模量U:

式中，K=(K1，K2，…，KN)，U=(U1，U2，…，UN)分別為N種礦物組分的體積模量和剪切模量;F=(f1，f2，…，fN)為各組分的體積分?jǐn)?shù);Φ =(φvu，φin，φcr，φcl)為四種孔隙的孔隙度，A=(αvu，αin，αcr，αcl)為四種孔隙的縱橫比;Sf為流體飽和度;ρs為飽和巖石的密度;vPm、vPe分別為實(shí)測(cè)和計(jì)算的縱波速度。

目標(biāo)函數(shù)式(9)取最小值的過(guò)程是一個(gè)非線性多元函數(shù)的求解問(wèn)題。本文采用自適應(yīng)遺傳算法［28］，利用目標(biāo)函數(shù)本身的信息建立尋優(yōu)方向，具備全局最優(yōu)搜索性，且具有內(nèi)在的并行計(jì)算特性。具體的計(jì)算流程如圖3所示。為了減少多解性，使用礦物組分彈性模量的經(jīng)驗(yàn)值范圍對(duì)待求參數(shù)進(jìn)行約束。

使用上述方法可以反演得到隨深度變化的礦物組分彈性模量。但是，該過(guò)程是建立在孔隙微結(jié)構(gòu)參數(shù)已知的基礎(chǔ)上的，而確定孔隙微結(jié)構(gòu)參數(shù)時(shí)，又需要礦物組分的彈性模量值作為輸入，這樣導(dǎo)致了礦物組分彈性模量和孔隙微結(jié)構(gòu)參數(shù)互為前提條件而又互相影響。為此，首先使用礦物組分彈性模量的經(jīng)驗(yàn)值來(lái)計(jì)算孔隙微結(jié)構(gòu)參數(shù)，據(jù)此計(jì)算礦物組分的彈性模量。最后利用計(jì)算得到的礦物組分彈性模量值和孔隙微結(jié)構(gòu)參數(shù)計(jì)算儲(chǔ)層的橫波速度。

圖3 礦物組分彈性模量反演流程Fig.3 Flow chart for mineral elastic modulus inversion

3 實(shí)際資料應(yīng)用效果分析

對(duì)某碳酸鹽巖研究區(qū)的實(shí)際測(cè)井資料進(jìn)行橫波速度計(jì)算。該研究區(qū)內(nèi)的碳酸鹽巖礦物主要是方解石和白云石，含有少量的泥質(zhì)。儲(chǔ)層內(nèi)的孔隙流體主要是水和天然氣。A井的已知測(cè)井?dāng)?shù)據(jù)如圖4所示，該井的實(shí)測(cè)橫波速度用于檢驗(yàn)橫波速度估算結(jié)果的質(zhì)量。

首先計(jì)算巖石的孔隙微結(jié)構(gòu)參數(shù)。用表1所示礦物組分彈性模量的經(jīng)驗(yàn)值，計(jì)算得到的孔洞、粒間孔隙、裂隙以及泥質(zhì)孔隙的縱橫比和孔隙度分別如圖5和圖6所示。計(jì)算的孔隙微結(jié)構(gòu)是實(shí)際孔隙的等效近似，可以在一定程度上反映巖石中不同形狀孔隙的組成情況。從圖中可以看出，5.34 km深度附近的地層，縱波速度降低，計(jì)算的粒間孔隙和裂隙縱橫比減小，而孔洞的孔隙度降低，粒間孔隙和裂隙的孔隙度增大，符合低縱橫比孔隙導(dǎo)致巖石縱波速度較低的理論趨勢(shì)。計(jì)算得到的橫波速度(圖7左圖)的低頻趨勢(shì)與實(shí)測(cè)縱橫波速度基本吻合，但是高頻細(xì)節(jié)上存在不可忽視的誤差。

圖4 A井已知測(cè)井?dāng)?shù)據(jù)Fig.4 Well logging data of well A

表1 礦物組分經(jīng)驗(yàn)彈性參數(shù)Table 1 Empirical elastic parameters of mineral components

圖5 不同孔隙的孔隙縱橫比Fig.5 Pore aspect ratio of different pores

根據(jù)本文中的方法反演方解石、白云石以及泥質(zhì)的彈性模量，結(jié)果如圖8所示。由于方解石含量較高，因此相對(duì)于白云石和泥質(zhì)來(lái)說(shuō)，方解石的彈性模量對(duì)縱波速度的影響較大。在不含白云石的深度上，無(wú)法反演求取白云石的彈性模量，為此假設(shè)其等于表1中所示的經(jīng)驗(yàn)值。從圖中可以看出，在5.46 km深度附近的地層，計(jì)算的橫波速度比實(shí)測(cè)橫波速度高，而反演得到的礦物組分彈性模量值比經(jīng)驗(yàn)值小，這與理論趨勢(shì)相吻合。最終計(jì)算的橫波速度如圖7中右圖所示，與圖7左圖所示的橫波速度計(jì)算結(jié)果相比，大大減小了由于礦物組分彈性模量不準(zhǔn)確引入的誤差，尤其是在5.46 km深度附近，橫波速度的低值得到了很好的反映。

圖6 不同孔隙的孔隙度Fig.6 Porosity of different pores

計(jì)算的縱波速度與實(shí)測(cè)值的相關(guān)系數(shù)為0.99，估算的橫波速度與實(shí)測(cè)值的相關(guān)系數(shù)為0.94，平均相對(duì)誤差為1.9%，平均絕對(duì)誤差為59.9 m/s，方差為77.1 m/s。計(jì)算得到的縱橫波速度以及泊松比與實(shí)測(cè)值的吻合程度較高且誤差接近于正態(tài)分布(圖9、10，圖10中曲線為高斯分布擬合曲線)證明了方法的有效性。如果結(jié)合測(cè)井解釋信息，可以進(jìn)一步進(jìn)行巖性或者流體巖石物理分析，為地震反演或者解釋提供有利的幫助和指導(dǎo)。

圖10 縱橫波速度以及泊松比估算結(jié)果誤差分布直方圖Fig.10 Error distribution histogram for P-wave velocity，S-wave velocity and Poission ratio

4 結(jié)束語(yǔ)

碳酸鹽巖儲(chǔ)層的非均質(zhì)性很強(qiáng)，進(jìn)行巖石物理分析和計(jì)算時(shí)，需要明確各種巖石物理理論模型的適用性，建立適用于復(fù)雜碳酸鹽巖儲(chǔ)層的巖石物理理論和方法。對(duì)碳酸鹽巖孔隙微結(jié)構(gòu)進(jìn)行合理的描述是計(jì)算橫波速度的關(guān)鍵。在缺少巖心數(shù)據(jù)或者試驗(yàn)數(shù)據(jù)的情況下，根據(jù)孔隙形狀和連通性對(duì)碳酸鹽巖孔隙進(jìn)行劃分，得到的孔隙微結(jié)構(gòu)是實(shí)際巖石孔隙系統(tǒng)的等效近似，可以在一定程度上反映巖石中不同形狀孔隙的構(gòu)成情況。碳酸鹽巖儲(chǔ)層的礦物組分非常復(fù)雜，使用巖石速度信息作為約束條件反演計(jì)算得到的不同深度點(diǎn)上的礦物組分彈性模量，可以在一定程度上反映儲(chǔ)層巖性的縱向變化。實(shí)際地層中的碳酸鹽巖由于裂隙的定向排列可能表現(xiàn)出各向異性特征，由于缺少相關(guān)的試驗(yàn)及測(cè)井信息，無(wú)法對(duì)裂隙的方位和密度等作出合理的描述，因此沒(méi)有考慮裂隙誘導(dǎo)各向異性的影響。相關(guān)的各向異性理論和方法有待進(jìn)一步研究。

［1］ AHR W M，ALLEN D，BOYD A，et.al.Confornting the carbonate conundrum［J］.Oilfield Review，2006，18-29.

［2］ 楚澤涵，陳豐，劉祝萍，等.估算地震橫波速度的新方法［J］.測(cè)井技術(shù)，1995，19(5):313-318.CHU Ze-han，CHEN Feng，LIU Zhu-ping，et al.Estimation of formation Vs［J］.Well Logging Technology，1995，19(5):313-318.

［3］ CASTAGNA J P，BATZLE M L，EASTWOOD R L.Relationships between compressional-wave and shear-wave velocities in clastic silicate rocks［J］.Geophysics，1985，50(4):571-581.

［4］ HAN D H，NUR A，MORGAN D.Effects of porosity and clay content on wave velocities in sandstones［J］.Geophysics，1986，51(11):2093-2107.

［5］ GASSMANN F. über die Elastizit?t poroser Medien［J］.Veirteljahrsschrift der Naturforschenden Gesellschaft in Zürich，1951，96:1-23.

［6］ GREENBERG M L，CASTAGNA J P.Shear-wave velocity estimation in porous rocks:theoretical formulation，prelimulation，preliminary verification and applications［J］.Geophysical Prospecting，1992，40(2):195-209.

［7］ XU S Y，WHITE R E.A new velocity model for claysand mixtures［J］.Geophysical Prospecting，1995，43(1):91-118.

［8］ 郭棟，王興謀，印興耀，等.橫波速度計(jì)算方法與應(yīng)用［J］.油氣地球物理，2007，5(3):18-22.GUO Dong，WANG Xing-mou，YIN Xing-yao，et al.The calculation method and application of transverse wave［J］.Petroleum Geophysics，2007，5(3):18-22.

［9］ 邵才瑞，印興耀，張福明，等.利用常規(guī)測(cè)井資料基于巖石物理和多礦物分析反演橫波速度［J］.地球科學(xué)——中國(guó)地質(zhì)大學(xué)學(xué)報(bào)，2009，34(4):699-707.SHAO Cai-rui，YIN Xing-yao，ZHANG Fu-ming，et al.Shear wave velocity inversion with routine well logs based on rock physics and multi-minerals nalysis［J］.Earth Science—Journal of China University of Geosciences，2009，34(4):699-707.

［10］ 張璐.基于巖石物理的地震儲(chǔ)層預(yù)測(cè)方法應(yīng)用研究［D］.青島:中國(guó)石油大學(xué)(華東)地球科學(xué)與技術(shù)學(xué)院，2009.ZHANG Lu.Application of rock physics theory in seismic reservior discrimination［D］.Qingdao:School of Geosciences in China University of Petroleum，2009.

［11］ 李呈呈.碳酸鹽巖巖石物理模型的研究及應(yīng)用［D］.青島:中國(guó)石油大學(xué)(華東)地球科學(xué)與技術(shù)學(xué)院，2011.LI Cheng-cheng.The research and application of rock physics model for carbonate rocks［D］.Qingdao:School of Geosciences in China University of Petroleum，2011.

［12］ 周中彪.基于巖石物理模型的測(cè)井約束橫波速度計(jì)算方法研究［J］.物探化探計(jì)算技術(shù)，2010，32(5):536-541.ZHOU Zhong-biao.The study on S-wave velocity calculation method constrained by logging based on rock physics models［J］.Computing Techniques for Geophysical and Geochemical Exploration，2010，32(5):536-541.

［13］ KUMAR M，HAN D H.Pore shape effect on elastic properties of carbonate rocks［C/OL］.Expanded Abstract of 75th Annual International SEG Meeting，2005，1477-1480.［2010-09-22］ http://dx.doi.org/10.1190/1.2147969.

［14］ BAECHLE G T，COLPAERT A，EBERLI G P et.al.Modeling velocity in carbonates using a dual porosity DEM model［C/OL］.Expanded Abstract of 77th AnnualInternationalSEG Meeting， 2007，1589-1593.［2011-06-23］ http://dx.doi.org/10.1190/1.2792799.

［15］ VEGA S，BERTEUSSEN K，Yue-feng SUN.Is Gassmann1the best model for fluid substitution in heterogeneous carbonates?［C/OL］.Expanded Abstract of 77th Annual International SEG Meeting，2007，1575-1578.［2012-03-29］ http://dx.doi.org/10.1190/1.2792796.

［16］ 陳颙，黃庭芳，劉恩儒.巖石物理學(xué)［M］.合肥:中國(guó)科學(xué)技術(shù)大學(xué)出版社，2009:17-20.

［17］ BERRYMAN J G.Mixture theories for rock properties［M］//AHRENS T J.A handbook of physical constants.Washington D C:American Geophysical U-nion，1995:205-228.

［18］ KUSTER G T，TOKS?Z M N.Velocity and attenuation of seismic waves in two-phase media［J］.Geophysics，1974，39(5):587-618.

［19］ BERRYMAN J G.Single-scattering approximations for coefficients in Biot1s equations of poroelasticity［J］.Journal of the Acoustical Society of America，1992，91:551-571.

［20］ BERRYMAN J G.Long-wavelength propagation in composite elastic media［J］.Acoustical Society of America，1980，68(6):1801-1831.

［21］ WHITE J E.Computed seismic speeds and attenuation in rocks with partial gas saturation［J］.Geophysics，2002，67:1406-1414.

［22］ 王炳章，朱曄，王丹.多孔介質(zhì)的流體機(jī)制模型及其頻散機(jī)理［J］.勘探地球物理進(jìn)展，2008，31(6):405-413.WANG Bing-zhang， ZHU Ye， WANG Dan.Fluid mechanism models and their velocity dispersions in porous media［J］.Progress in Exploration Geophysics，2008，31(6):405-413.

［23］ ANSELMETTI F S，EBERLI G P.The velocity-deviation log:a tool to predict pore type and permeability trends in carbonate drill holes from sonic and porosity or density logs［J］.AAPG Bulletin，1999，83(3):450-466.

［24］ WEGER R J，BAECHLE G T，MASAFERRO J L，et al.Effects of porestructure on sonic velocity in carbonates［C/OL］.Expanded Abstract of 74th Annual International SEG Meeting，2004，23(1):1774-1777.［2011-10-29］ http://dx.doi.org/10.1190/1.1845169.

［25］ LUCIA F J.Rock-fabric/petrophysical classification of carbonate pore space for reservoir characterization［J］.AAPG Bulletin，1995，79(9):1275-1300.

［26］ L?N?Y A.Making sense of carbonate pore systems［J］.AAPG Bulletin，2006，90(9):1381-1405.

［27］ ANSELMETTI F S，EBERLI G P.Controls on sonic velocity in carbonates［J］.Pure and Applied Geophysics，1993，141(2/4):287-323.

［28］ SIVANANDAM S N，DEEPA S N.Introduction to genetic algorithms［M］.New York:Springer Publishing Company，2007.