郭繼亮，牛濱華，孫春巖，劉 暢

1 中國地質(zhì)大學(xué)（北京），地球物理與信息技術(shù)學(xué)院，北京 100083

2 中國地質(zhì)大學(xué)（北京），工程技術(shù)學(xué)院，北京 100083

1 引 言

本文使用彈性固體模型，彈性流體模型和臨界孔隙度介質(zhì)模型，通過運用等效用力學(xué)模型理論，建立流體飽和線性彈性孔隙模型.并結(jié)合波場模擬，研究該介質(zhì)模型中地震波傳播的特征和規(guī)律.對孔隙介質(zhì)中地震波傳播的研究由來已久.Biot理論被視為研究流體飽和孔隙介質(zhì)彈性波傳播的基礎(chǔ).Biot，M.A.承襲了 Von Terzaghi[1］的理念，將流體飽和孔隙介質(zhì)作為一個整體來處理.Biot理論涵蓋了高頻波和低頻波在流體飽和孔隙介質(zhì)中的轉(zhuǎn)播[2-3］.一些學(xué)者應(yīng)用Biot理論對地震波在孔隙介質(zhì)中的反射和透射進行了研究[4-5］.Zhu等[6］根據(jù)Biot理論，使用二維有限差分方法模擬了孔隙彈性儲層的地震響應(yīng).Berryman和 Wang[7-8］發(fā)展了雙孔介質(zhì)理.Pride等[9］結(jié)合Biot理論，討論了三種情況下P波在沉積巖中的衰減情況.Li等[10］基于Biot理論給出了二維在飽和圓弧形沉積河谷SV波散射問題的解析解.Chand等[11］比較了包括三相Biot模型的不同模型對含水合物沉積巖速度的描述.Dong等[12］結(jié)合Biot理論推導(dǎo)了三維半球形凹陷飽和土場地對平面P波散射問題的解析解.Wei等[13］以Biot理論為基礎(chǔ)，討論了黏滯流體飽和各向異性孔隙介質(zhì)中彈性波的頻散和衰減.Biot理論和以其為基礎(chǔ)的研究取得了很大的成就，但這套理論的不利因素是參數(shù)較多[14］.

等效力學(xué)模型通過基本力學(xué)模型不同方式的組合，得到目標介質(zhì)的應(yīng)力應(yīng)變關(guān)系，即本構(gòu)關(guān)系.被廣泛的應(yīng)用于描述黏彈性介質(zhì)中的應(yīng)力應(yīng)變關(guān)系，例如 Maxwell模型[15］和 Kelvin-Voigt模型[16-17］.臨界孔隙度模型出現(xiàn)于20世紀90年代，Nur通過對大量巖石樣品物性測量和分析提出了臨界孔隙度和臨界孔隙度模型[18-19］.此后 Mavko等[20］，Dvorkin等[21］，Mukerji等[22］和 Chen等[23］學(xué)者對臨界孔隙度模型開展了進一步的研究.在此基礎(chǔ)上牛濱華等[24-26］提出并實現(xiàn)了基于臨界孔隙度模型各組分彈性參數(shù)的數(shù)值計算方法；并以此為基礎(chǔ)進行了孔隙度介質(zhì)對彈性模量的研究，提出了基于臨界孔隙度模型的地震波傳播研究的思路和方法[27-32］.本文承襲這個思路繼續(xù)探討孔隙度對介質(zhì)彈性模量和地震波傳播的影響.

本文研究工作的前提為：基本固體模型，流體模型，臨界孔隙介質(zhì)模型和目標介質(zhì)模型為均勻彈性各向同性，流體與骨架完全耦合；基本力學(xué)模型和目標介質(zhì)力學(xué)模型中的應(yīng)力應(yīng)變都是線性變化的.

2 基本理論

利用等效力學(xué)模型和臨界孔隙度模型，建立流體飽和均勻彈性各向同性孔隙介質(zhì)中的本構(gòu)關(guān)系，進而建立該介質(zhì)模型中的波動方程，研究孔隙和孔隙流體對介質(zhì)彈性性質(zhì)的影響.

2.1 基本力學(xué)模型和目標力學(xué)模型

在孔隙度從0到1的區(qū)間內(nèi)，分布著三個基本力學(xué)模型和一個目標力學(xué)模型，如圖1所示.兩個固體型矢量彈性力學(xué)模型分別是孔隙度為0的S模型和孔隙度φE=φ的目標模型E；兩個流體型標量彈性孔隙模型分別是孔隙度為φC的臨界孔隙度模型C和孔隙度為1的流體模型F.其中S和F為骨架礦物和孔隙流體的成分，C代表臨界孔隙度條件下的介質(zhì)模型.

4個力學(xué)模型分布在兩個孔隙度范圍內(nèi)：一，固體型矢量彈性區(qū)間，即孔隙度分布于 [0～φC）（模型S和C之間）的介質(zhì).這時介質(zhì)整體呈固相，由礦物骨架承受壓力起支撐作用.二，流體型標量彈性區(qū)間，即孔隙度分布于（φC～1］（模型C與F之間）的介質(zhì)，這時介質(zhì)整體呈流相，骨架顆粒在流體中呈懸浮態(tài)，流體承受壓力起支撐作用.孔隙度為φC的C模型是固體型和流體型兩類彈性孔隙流體介質(zhì)的臨界點，或稱其為轉(zhuǎn)換點.四種力學(xué)模型的彈性參數(shù)如表1所示.

圖1 基本力學(xué)模型與目標模型之間關(guān)系示意圖Fig.1 Relation between characteristic elements S、C，F(xiàn) and target element E

表1 四個力學(xué)模型的彈性屬性參數(shù)表Table 1 Elastic parameters of the four elements

2.2 目標模型本構(gòu)關(guān)系

運用等效力學(xué)模型，建立基本力學(xué)模型與目標力學(xué)模型的彈性參數(shù)之間的聯(lián)系，即用基本力學(xué)模型的應(yīng)力應(yīng)變描述目標模型的應(yīng)力應(yīng)變關(guān)系.這時基本力學(xué)模型代表的骨架礦物、孔隙流體為流體飽和孔隙介質(zhì)的組成部分，因此基本力學(xué)模型的彈性參數(shù)稱為組分彈性參數(shù).等效力學(xué)模型中的經(jīng)典方法是：彈性模型與黏性模型的并聯(lián)網(wǎng)絡(luò)是等應(yīng)變系統(tǒng)，該網(wǎng)絡(luò)系統(tǒng)生成黏彈性固體型Voigt-Kelvin介質(zhì)；彈性模型與黏性模型的串聯(lián)網(wǎng)絡(luò)是等應(yīng)力系統(tǒng)，該網(wǎng)絡(luò)系統(tǒng)生成黏彈性流體型Maxwell介質(zhì).文中利用在不同相態(tài)區(qū)間內(nèi)的孔隙介質(zhì)組分之間的等應(yīng)力和等應(yīng)變關(guān)系，通過等效力學(xué)模型，建立流體飽和孔隙介質(zhì)中的本構(gòu)關(guān)系.

流體飽和孔隙介質(zhì)中，當介質(zhì)處于矢量彈性區(qū)間內(nèi)時，礦物骨架承壓，礦物骨架模型和孔隙流體模型的應(yīng)變相同，類似Voigt-Kelvin介質(zhì)中的固體彈性模型和黏性流體模型.當介質(zhì)模型是流體處在標量彈性區(qū)間內(nèi)時，骨架礦物懸浮于流體中，流體承壓，礦物骨架模型與流體模型所受應(yīng)力相同，類似Maxwell介質(zhì)中的固體彈性模型和黏性流體模型.可以視骨架礦物和孔隙流體為對應(yīng)的單相固體矢量彈性介質(zhì)和單相流體標量彈性介質(zhì)，結(jié)合臨界孔隙度模型，運用等效力學(xué)模型方法構(gòu)建目標孔隙流體介質(zhì)模型.

3 目標孔隙介質(zhì)模型構(gòu)建

目標介質(zhì)模型E的體積模量和剪切模量的構(gòu)建都包含了一次模型等應(yīng)力組合和一次模型等應(yīng)變組合.如圖1所示，通過兩步建立模型E的彈性模量組分表達式.第一步：依據(jù)等應(yīng)力條件，將S和F彈性模量加權(quán)組合得到C的彈性模量；第二步：依據(jù)等應(yīng)變條件，將C與S的彈性模量加權(quán)組合得到E的彈性模量.

3.1 力學(xué)模型的應(yīng)力與應(yīng)變關(guān)系

下面介紹四種介質(zhì)中的應(yīng)力應(yīng)變關(guān)系.模型S和E為固體型矢量彈性介質(zhì)，其體積模量由下述式子表示：

式中，模型E的壓力PE與三個正應(yīng)力的關(guān)系是體積應(yīng)變系數(shù)θE與三個正應(yīng)變的關(guān)系是同理，對于模型S有，模型S和E的剪切模量由下述式子表示：

模型C和F為流體型標量彈性介質(zhì)，其體積模量可以寫成：

模型C和F的剪切模量均為0，即，

上述（1）至（4）式的應(yīng)力與應(yīng)變關(guān)系都可以視為本構(gòu)關(guān)系的一種表達形式.

3.2 體積模量的組分表達形式

在均勻各向同性條件下，體積模量、剪切模量和密度可以完整表現(xiàn)介質(zhì)的彈性屬性.體積模量和剪切模量的定義表達式是本構(gòu)關(guān)系的一種表達形式，因此可以用這樣的表達式刻畫相關(guān)模型的彈性.

首先求取臨界體體積模量的組分表達形式.彈性流體模型C可以由彈性固體模型S和彈性流體模型F的等應(yīng)力組合表示.這種等應(yīng)力組合表明，C的體積應(yīng)變等于S和F兩者體積應(yīng)變的加權(quán)組合.等應(yīng)力組合的應(yīng)變加權(quán)系數(shù)分別為（1-φC）和φC，于是有

兩個力學(xué)模型的串聯(lián)系統(tǒng)是等應(yīng)力系統(tǒng)，即PC=則有

上式是KC關(guān)于KS和KF的組分方程.方程（6）的形式與Nur給出的臨界點KC的體積模量表達形式一致，這也從等效力學(xué)模型角度，對臨界理論提供了一種解釋.

下面推導(dǎo)目標模型體積模量的組分表達形式.彈性固體型E可以由彈性固體模型S和彈性流體模型C的等應(yīng)變組合表示.這種等應(yīng)變組合表明，E的正應(yīng)力等于S和C兩者正應(yīng)力的加權(quán)組合.應(yīng)力加權(quán)系數(shù)分別為于是有，

式中孔隙度φ滿足0≤φ＜φC.S和C兩個力學(xué)模型等應(yīng)變組合，即θE=θS=θC；則有

式（8）稱為KE的組分表達形式，即KE的組分方程.

3.3 其它參數(shù)的組分表達形式

運用類似的方法得到E的剪切模量的組分表述式.E的切應(yīng)力等于S和C兩者切應(yīng)力的加權(quán)組合.應(yīng)力組合的加權(quán)系數(shù)分別為于是有，

式中孔隙度φ滿足0≤φ＜φC.把（2）式和（4）式代入（9）式，得到μE的組分表達形式

目標介質(zhì)模型E的密度ρE表達式就是眾所周知的密度組分方程即，

由目標介質(zhì)模型彈性參數(shù)的組分表達形式可得到這種介質(zhì)中縱波和橫波的相速度組分表達形式.縱波和橫波的相速度組分表達式為

上述過程得到的體積模量和煎切模量的組分表達式與Nur等提出的臨界孔隙介質(zhì)模型中的彈性參數(shù)表達式完全一致；這也驗證了Nur等臨界孔隙度模型的正確性.

3.4 本構(gòu)方程和彈性矩陣

得到介質(zhì)模型E彈性參數(shù)的組分表達形式后，討論該介質(zhì)模型的本構(gòu)關(guān)系.由于介質(zhì)為均勻彈性各向同性，本構(gòu)方程的組分形式可以通過求取該介質(zhì)彈性矩陣中彈性參數(shù)的組分表達形式的方法得到，也可直接由推導(dǎo)該模型中介質(zhì)應(yīng)力應(yīng)變關(guān)系得到.

3.4.1 彈性矩陣組分表達形式

基于前面得到的模型E的體積模量KE、剪切模量μE的表達式，根據(jù)均勻彈性各向同性介質(zhì)彈性矩陣的結(jié)構(gòu)特點，可以得到E的彈性矩陣CE6×6的表達式：

模型E的本構(gòu)方程可以寫成：

3.4.2 應(yīng)力應(yīng)變關(guān)系的組分表達形式

針對給的4個力學(xué)模型的本構(gòu)關(guān)系，運用等效力學(xué)原理，討論目標介質(zhì)E的本構(gòu)關(guān)系和彈性矩陣.介質(zhì)模型E是通過模型S和C通過等應(yīng)變組合構(gòu)建而成，E所受的應(yīng)力為S和C的應(yīng)力線形組合，所以存在以下關(guān)系

由上面的本構(gòu)關(guān)系可以得到均勻彈性流體飽和孔隙介質(zhì)中的波動方程，其具體形式為

4 數(shù)值模擬

為檢驗孔隙對地震波傳播的影響，基于矢量彈性流體飽和孔隙介質(zhì)的波動方程，運用旋轉(zhuǎn)交錯網(wǎng)格有限差分方法對二維孔隙介質(zhì)模型（XOZ面）進行彈性波場數(shù)值模擬.波場模擬中模型組分的彈性參數(shù)見表2.飽水模型的彈性參數(shù)參考了Berger在Lawrence Livermore National Laboratory進行的流體飽和人工合成砂巖孔隙度-波速關(guān)系的物理模型試驗[33］.

表2 飽水介質(zhì)在組分的彈性參數(shù)Table 2 Elastic parameters of components of porous media

根據(jù)表2中的骨架和流體的彈性信息以及臨界孔隙度值，通過臨界孔隙度模型計算，可以得到不同孔隙度條件下目標模型的彈性參數(shù)，如圖2所示，根據(jù)方程（8），（10），（11）和（12），隨著孔隙度的增長，介質(zhì)的整體彈性參數(shù)隨之下降，相應(yīng)的整體速度也隨之減小.當孔隙度小于在臨界值時，隨著孔隙度增大整體體積模量（除臨界點）、整體剪切模量和整體密度的變化正比于孔隙度變化；而速度的變化在靠近臨界孔隙度的區(qū)域更加明顯，即減小得更快.這符合臨界孔隙度模型對介質(zhì)的描述.整體縱波波速從孔隙度為0時的6000m／s降到孔隙度為臨界值時的1597.4m／s，減少了73.4%；相應(yīng)的整體橫波速度由3464.2m／s降到0m／s，減少了100%；相應(yīng)的整體密度值由2500kg／m3降到1772.9kg／m3，減少了29.1%.

圖2 目標模型孔隙度與彈性參數(shù)關(guān)系圖實線圖對應(yīng)模量坐標，其中BM代表體積摸量，SM代表剪切模量，虛線圖對應(yīng)速度密度坐標，其中D代表密度，VP代表縱波速度VS代表橫波速度.Fig.2 Relation between porosity and elastic parameters of E Solid lines show modulus；BMmeans bulk modulus and SM means shear modulus，dash lines show velocities and densities；D，VPand VSrespectively mean density Pwave velocities and S wave velocities.

如圖3a所示，模型被離散為網(wǎng)格距為2m的方形網(wǎng)格，共512×1024個網(wǎng)格點.模型分為兩層，層位置為z=512m，在模型的上半部分，也就是I區(qū)域中，介質(zhì)的孔隙度為0，在模型的下半部分，即區(qū)域II中，介質(zhì)的孔隙度分別為0.1，0.2，0.3，0.4和0.4847，孔隙中飽水.位于坐標（1024，10）的震源為主頻60Hz的Ricker應(yīng)力子波，采用不同垂向應(yīng)力和切向應(yīng)力組合的縱橫波偶合震源.模擬的時間步長為0.2ms.兩個觀測點的坐標為（600，0）和（600，600）.采用時間二階精度，空間四階精度的算法.

分別對以上5個不同孔隙度的模型進行數(shù)值模擬，t=190ms時X分量的波場快照如圖3（b-f）所示.

由圖3可以看出，波在介質(zhì)I和介質(zhì)II的分界面發(fā)生反射和透射.隨著孔隙度的增加，介質(zhì)II中的整體縱橫波速度變小.另外，孔隙度為0，0.1，0.2，0.3，0.4對應(yīng)的波場3（b-e）中II區(qū)域中同時存在縱波和橫波，因為這時的介質(zhì)都是矢量彈性的；而孔隙度為臨界值φC=0.4847的圖3f的波場中僅有縱波成分，因為這時的介質(zhì)是標量彈性的.孔隙度增大過程中反射波的能量逐漸增加，透射波的能量逐漸減小.

兩個觀測點的地震記錄如圖4所示，地表接收點（坐標（600，0））位于區(qū)域I中，接收反射記錄；地下接收點（坐標（600，600））位于區(qū)域II中，接收透射記錄.反射記錄中三個波形從左至右依次為PP波反射，PS+SP波反射和SS波反射記錄，由于反射波在介質(zhì)區(qū)域I中傳播，不同模型中反射波的走時相同，但由于孔隙度導(dǎo)致的分界面波阻抗差異，反射波振幅隨區(qū)域II中介質(zhì)孔隙度增大而增大如圖4a所示.透射記錄中除臨界孔隙介質(zhì)外各個記錄中均存在PP，PS，SP和SS四種波的波形，臨界孔隙介質(zhì)中只含PP和SP兩種波形，彈性波走時隨區(qū)域II中介質(zhì)孔隙度增大而增大，如圖4b所示.

通過對孔隙和孔隙流體導(dǎo)致的空間彈性參數(shù)差異的模型進行數(shù)值模擬，直觀的反映了孔隙和孔隙流體對彈性波傳播的影響.也驗證了臨界孔隙度模型地震建模的效果.

5 結(jié) 論

本文利用等效力學(xué)模型原理，將彈性固體模型、彈性流體模型和臨界孔隙度模型進行線性組合，得到了彈性（即無黏性）飽和流體孔隙介質(zhì)模型，主要結(jié)果如下：

（1）運用等效力學(xué)模型構(gòu)建流體飽和組分彈性孔隙介質(zhì)的本構(gòu)關(guān)系.本文基于孔隙介質(zhì)組分的彈性關(guān)系，利用臨界孔隙度模型，通過應(yīng)用等效力學(xué)模型，建立了流體飽和組分彈性孔隙的本構(gòu)關(guān)系.本構(gòu)關(guān)系中彈性矩陣的具體形式可以通過求取該介質(zhì)彈性矩陣中彈性參數(shù)的組分表達形式的方法得到，也可直接由推導(dǎo)該模型中介質(zhì)應(yīng)力應(yīng)變關(guān)系得到.

（2）在構(gòu)建模型的過程中引入了臨界孔隙度模型，并將臨界孔隙介質(zhì)模型作為一個基本力學(xué)模型單元.臨界孔隙模型的引進不僅給出了孔隙介質(zhì)的相態(tài)變化點，也將礦物骨架模型和孔隙流體模型有機的聯(lián)系在一起，將臨界孔隙介質(zhì)模型作為目標模型的直接組成部分.利用等效力學(xué)模型得到的臨界孔隙介質(zhì)的彈性參數(shù)表達其形式與Nur給出的經(jīng)驗公式形式一致，這也從等效力學(xué)模型角度，證明了線性應(yīng)力應(yīng)變前提下，臨界孔隙模型的正確性.

（3）通過數(shù)值模擬，在組分模型的前提下，初步探討了孔隙和孔隙流體對地震波傳播的影響.通過數(shù)值模擬初步表現(xiàn)了介質(zhì)中孔隙對反射波和透射波傳播速度和振幅速度的影響.

6 討 論

實際介質(zhì)中孔隙和孔隙流體對地震波傳播有多方面影響.本文通過臨界孔隙度模型討論完全彈性條件下孔隙和孔隙流體對地震波傳播的影響.未考慮孔隙流體的黏性，流體與骨架耦合關(guān)系和孔隙的散射作用等因素對地震波的影響，這也是我們下一步研究的重點.本文中的均勻彈性各項同性前提下的工作不但是研究流體黏性對地震波研究的前提，也是研究各向異性孔隙介質(zhì)地震波研究的基礎(chǔ).從本構(gòu)關(guān)系角度講，組分模型地震波傳播的工作可以很好的和經(jīng)典黏彈性介質(zhì)和各向異性介質(zhì)地震波傳播的工作相結(jié)合.

（References）

[1]Terzaghi K.Auftrieb und Kapillardruck an betonierten Talsperren.Die Wasserwirtschaft，1933，（26）：397-399.

[2]Biot M A.Theory of propagation of elastic waves in a fluidsaturated porous solid：I.Low-frequency range.Acoust.Soc.Am.，1956，28（2）：168-178.

[3]Biot M A.Theory of propagation of elastic waves in a fluidsaturated porous solid：II.higher frequency range.Acoust.Soc.Am.，1956，28（2）：179-191.

[4]Hajra S，Mukhopadhyay A.Reflection and refraction of seismic waves incident obliquely at the boundary of a Liquidsaturated porous solid.Bulletin of the Seismological Society of America，1982，72（5）：1509-1533.

[5]Deresiewicz H.The effect of boundaries on wave propagation in a Liquid-filled porous solid：IV.surface waves in a halfspace.Bull.Seism.Soc.Am，1962，52（3）：627-638.

[6]Zhu X，McMechan G A.Numerical simulation of seismic responses of poroelastic reservoirs using biot theory.Geophysics，1991，56（3）：328-339.

[7]Berryman J G，Wang H F.The elastic coefficients of doubleporosity models for fluid transport in jointed rock.Journal of Geophysical Research-Solid Earth，1995，100（B12）：24611-24627.

[8]Berryman J G，Wang H F.Elastic wave propagation and attenuation in a double-porosity dual-permeability medium.International Journal of Rock Mechanics and Mining Sciences，2000，37（1-2）：63-78.

[9]Pride S R，Berryman J G，Harris J M.Seismic attenuation due to wave-induced flow.Journal of Geophysical Research-Solid Earth，2004，109：B01201，doi：10.1029／2003JB002639.

[10]李偉華，趙成剛.飽和土沉積谷場地對平面SV波的散射問題的解析解.地球物理學(xué)報，2004，47（5）：911-919.Li W H，Zhao C G.Scattering of plane SV waves by circular arc alluvial valleys with saturated soil deposits.Chinese J.Geophys.，2004，47（5）：911-919.

[11]Chand S，Minshull T A，Gei D，et al.Elastic velocity models for gas-hydrate-bearing sediments—a comparison.Geophysical Journal International，2004，159（2）：573-590.

[12]董俊，趙成剛.三維半球形凹陷飽和土場地對平面P波散射.地球物理學(xué)報，2005，48（3）：680-688.Dong J，Zhao C G.An analytic solution for the diffraction of plane P-wave by three-dimensional hemispherical canyons in a fluid-saturated porous media half space.Chinese J.Geophys.（in Chinese），2005，48（3）：680-688.

[13]魏修成，盧明輝，巴晶等.含黏滯流體各向異性孔隙介質(zhì)中彈性波的頻散和衰減.地球物理學(xué)報，2008，51（1）：213-220.Wei X C，Lu M H，Ba J，et al.Dispersion and attenuation of elastic waves in a viscous fluid-saturated anisotropic porous solid.Chinese J.Geophys.（in Chinese），2008，51（1）：213-220.

[14]陳颙，黃庭芳，劉恩儒.巖石物理學(xué).安徽：中國科學(xué)技術(shù)大學(xué)出版社，2009.Chen Y，Huang T F，Liu E R.Rock Physics（in Chinese）.Anhui：Press of University of Science and Technology of China，2009.

[15]Maxwell J C.On the dynamical theory of gases.Phil.Trans.Roy.Soc.London，1867，157：49-88.

[16]Kelvin L.Vortex statics.Collected works，1875，IV：115-128.

[17]Voigt W.Uber innere reibung fester korper，insbesondere der metalle.Annalen der Physik，1892，283（12）：671-693.

[18]Nur A.Critical porosity and the seismic velocities in rocks.EOS.Trans.Am.Geophys Union，1992，73：43-66.

[19]Nur A，Marion D，Yin H.Wave velocites in sediments.／／Hovem J M，Richardson M D，Stoll R D，eds.Shear Wave in Marine Sediments.Dordrecht：Kluwer Academic Publisher，1991：131-140.

[20]Mavko G， Mukerji T，Dvorkin J.The Rock Physics Handbook.New York：Cambridge University Press，1998.

[21]Dvorkin J，Prasad M，Sakai A，et al.Elasticity of marine sediments：Rock physics modeling.Geophysical Research Letters，1999，26（12）：1781-1784.

[22]Mukerji T，Berryman J， Mavko G，et al.Differential effective medium modeling of rock elastic moduli with critical porosity constraints.Geophysical Research Letters，1995，22（5）：555-558.

[23]Chen Q，Nur A.Critical concentration models for porous materials.／／ Corapcioglu M Y，ed.Advance in Porous Media.New York：Elsevier.1994：169-308.

[24]牛濱華，孫春巖，閆國英等.含氣介質(zhì)臨界點，流體和骨架彈性參數(shù)的數(shù)值計算方法.地球物理學(xué)報，2010，53（6）：1495-1501.Niu B H，Sun C Y，Yan G Y，et al.Numerical calculation method for getting critical point，pore fluid and framework parameters of gas-bearing media.Chinese J.Geophys.（in Chinese），2010，53（6）：1495-1501.

[25]牛濱華，郭繼亮，孫春巖.CLYJ飽水石灰?guī)r骨架和流體彈性參數(shù)估算.地球物理學(xué)報.2011，54（10）：2652-2664.Niu B H，Guo J L，Sun C Y，et al.Estimating elastic parameters of framework and pore fluid from CLYJ saturated limestone samples.Chinese J.Geophys.，2011，54（10）：2652-2664.

[26]Niu B，Sun C Y，Yan G Y，et al.Linear numerical calculation method for obtaining critical point，pore fluid，and framework parameters of gas-bearing media.Applied Geophysics，2009，6（4）：319-326.

[27]賈冀輝，牛濱華，孫春巖等.水合物密度的數(shù)值計算及儲層識別.現(xiàn)代地質(zhì)，2011，25（6）：1180-1186.Jia J H，Niu B H，Sun C Y，et al.Density calculation and reservoir identification of gas hydrates.Geoscience（in Chinese），2011，25（6）：1180-1186.

[28]牛濱華，孫春巖，李明.半空間介質(zhì)與地震波傳播.北京：石油工業(yè)出版社，2002.Niu B H，Sun C Y，Li M.Seismic Wave Propagation in Half Space Media （in Chinese）.Beijing：Petroleum Industry Press，2002.

[29]牛濱華，孫春巖.固體彈性介質(zhì)與地震波傳播：半空間均勻各向同性單相.北京：地質(zhì)出版社，2005.Niu B H，Sun C Y.Seismic Wave Propagation in Single Phase Homogeneous Isotropic Half Space Media（in Chinese）.Beijing：Geological Puldishing House，2005.

[30]牛濱華，孫晟，孫春巖等.Biot介質(zhì)密度參數(shù)的容差密度表達.現(xiàn)代地質(zhì)，2007，21（3）：551-556.Niu B H，Sun S，Sun C Y，et al.Expressions of density Parameters of Biot medium with differential density.Geoscience（in Chinese），2007，21（3）：551-556.

[31]鄒荃，牛濱華，孫春巖等.飽水石灰?guī)r骨架和流體彈性參數(shù)的估算.石油地球物理勘探，2011，46（5）：763-768.Zou Q，Niu B H，Sun C Y，et al.Estimation of component elastic parameters from saturated limestone samples.Oil Geophysical Prospecting （in Chinese），2011，46（5）：763-768.

[32]張艷，牛濱華，孫春巖等.基于典型飽水砂巖樣品的多維空間多元統(tǒng)計密度反演方法.現(xiàn)代地質(zhì)，2011，25（6）：1174-1179.Zhang Y，Niu B H，Sun C Y，et al.Density inversion of typical water-saturated sandstones with multivariate statistical analysis in multidimensional space.Geoscience（in Chinese），2011，25（6）：1174-1179.

[33]Berge P A，Bonner B P，Berryman J G.Ultrasonic velocityporosity relationships for sandstone analogs made from fused glass beads.Geophysics，1995，60（1）：108-119.