劉向君，朱洪林，梁利喜

油氣藏地質(zhì)及開(kāi)發(fā)工程國(guó)家重點(diǎn)實(shí)驗(yàn)室，西南石油大學(xué)，成都 610500

1 引言

巖石是一種天然的多孔介質(zhì)，其內(nèi)部除了固體基質(zhì)還分布有大量不規(guī)則的孔隙以及孔隙空間流體，這些組分的物理性質(zhì)以及微觀孔隙結(jié)構(gòu)特征直接影響著宏觀巖石物理屬性，如強(qiáng)度、彈性模量、滲透率、電阻率、聲波速度等.模擬孔隙尺度的物理現(xiàn)象、理解微觀作用機(jī)理是準(zhǔn)確獲取巖石物理性質(zhì)的關(guān)鍵，探明巖石微觀組構(gòu)與宏觀物性之間的內(nèi)在聯(lián)系，對(duì)于解決石油、地質(zhì)等地球物理領(lǐng)域中的實(shí)際工程問(wèn)題具有十分重要的意義，而這一切僅靠傳統(tǒng)巖石物理研究手段是實(shí)現(xiàn)不了的.近年來(lái)，國(guó)內(nèi)外已有學(xué)者（?ren and Bakke，2002；Arns etal.，2004a，2004b；Hu，2007；Okabe and Blunt，2005；Zhao etal.，2007）通過(guò)多種方法建立了能夠刻畫(huà)孔隙空間分布的三維數(shù)字巖芯，在此基礎(chǔ)上開(kāi)展數(shù)值模擬，從而計(jì)算巖石物性參數(shù).這種方法被稱為數(shù)字巖石物理，由于研究是基于數(shù)字化平臺(tái)的虛擬實(shí)驗(yàn)，因而具有可重復(fù)性，可同時(shí)模擬多重物理響應(yīng)并探討相互關(guān)系，且微觀影響因素可控，還能模擬傳統(tǒng)巖石物理實(shí)驗(yàn)難以測(cè)量的物理量，并節(jié)省大量人力物力資源.數(shù)字巖石物理的這些優(yōu)勢(shì)，使其逐漸成為地球物理學(xué)的研究熱點(diǎn).

盡管如此，現(xiàn)有的主要研究成果還集中在國(guó)外，主要來(lái)自于挪威的Numerical Rock團(tuán)隊(duì)（2002，2011）、澳大利亞國(guó)立大學(xué)的Arns（2004）、Knackstedt等人（2002）、英國(guó)帝國(guó)理工學(xué)院的Okabe等人（2005）、Hu等人（2007），美國(guó)斯坦福大學(xué)的Keehm（2003）、Sain（2010），德 國(guó) 卡 爾 斯 魯 厄 大 學(xué) 的Saenger等人（2004，2008）以及美國(guó)的數(shù)字巖石物理公司Ingrain（2010）.而國(guó)內(nèi)目前還處于方興未艾的階段，中國(guó)石油大學(xué)的Zhao等人（2007）、Liu等人（2009a，2009b）、陶果等人（2005）、岳文正等人（2004），以及西南石油大學(xué)的Su等人（2010）開(kāi)展了相應(yīng)研究，取得了一定的成果.由于滲流機(jī)理在提高采收率中的重要地位，前面的研究大多集中于滲流特性模擬，且均采用格子玻茲曼方法或基于帝國(guó)理工的兩相流代碼；對(duì)巖石聲、電、彈性性質(zhì)的研究還比較零散，且其中的數(shù)值模擬大都基于Garboczi教授的開(kāi)源代碼；在數(shù)字巖芯建模方法的選擇上，由于微CT成本太高，多數(shù)學(xué)者基于二維薄片信息采用數(shù)學(xué)方法進(jìn)行三維重構(gòu)，而這會(huì)導(dǎo)致微觀孔隙結(jié)構(gòu)過(guò)于理想化或隨機(jī)化，無(wú)法反映真實(shí)；在研究對(duì)象的選擇上，幾乎都以澳大利亞國(guó)立大學(xué)提供的楓丹白露砂巖或Bera砂巖數(shù)字巖芯為載體（兩者可視為均質(zhì)純砂巖，結(jié)構(gòu)簡(jiǎn)單），而對(duì)復(fù)雜巖石研究缺乏.但總的來(lái)說(shuō)，前面的研究無(wú)論是側(cè)重于數(shù)字巖芯的三維重構(gòu)或孔隙網(wǎng)絡(luò)模型構(gòu)建，還是后期的巖石物理數(shù)值模擬分析，都為推動(dòng)數(shù)字巖石物理這一新技術(shù)的發(fā)展做出了不可磨滅的貢獻(xiàn)，只是研究手段還可以再豐富些、研究?jī)?nèi)容還應(yīng)該更為全面、系統(tǒng)化.綜合上述分析，本文以常規(guī)砂巖樣品為例，通過(guò)微CT掃描結(jié)合先進(jìn)的三維可視化軟件Avizo建立了具有真實(shí)孔隙結(jié)構(gòu)特征的三維數(shù)字巖芯模型，在此基礎(chǔ)上，利用Avizo強(qiáng)大的幾何模型前處理、后處理功能，將Avizo與多場(chǎng)耦合有限元軟件Comsol完美對(duì)接，實(shí)現(xiàn)了多種巖石物理參數(shù)的數(shù)值模擬，從而在避免繁瑣的算法研究和程序開(kāi)發(fā)的同時(shí)，為數(shù)字巖石物理的大規(guī)模發(fā)展提供了一條新的途徑，也為該領(lǐng)域的研究人員提供了一套可借鑒的研究思路.

2 三維數(shù)字巖芯建模

2.1 微CT成像

微CT掃描作為一種無(wú)損檢測(cè)物體內(nèi)部結(jié)構(gòu)的技術(shù)，是當(dāng)前建立三維數(shù)字巖芯最直接和最準(zhǔn)確的方法，其原理是根據(jù)巖石中不同密度的成分對(duì)X射線吸收系數(shù)不同以達(dá)到區(qū)分孔隙和骨架的目的.本研究中巖芯三維圖像的采集均在美國(guó)Xradia公司生產(chǎn)的MicroXCT-400（圖1a）上完成，其最高采樣分辨率可達(dá)1μm.實(shí)驗(yàn)樣品為直徑約8mm的圓柱體砂巖（圖1b），一個(gè)樣品可獲取983張980×1005像素的二維CT切片圖，空間分辨率為2.1μm／體素，將這些二維切片圖依次疊加組合便得到巖樣的三維灰度圖像.圖1c為其中一張切片的灰度圖，灰色、白色的巖石骨架（高密度）和黑色的孔隙（低密度）在圖像中清晰可辨.

2.2 圖像處理

微CT掃描獲得的巖芯灰度圖像中存在各種類(lèi)型的系統(tǒng)噪聲，降低圖像質(zhì)量的同時(shí)也不利于后續(xù)的定量分析，因此圖像處理第一步是通過(guò)濾波算法增強(qiáng)信噪比.針對(duì)三維圖像，比較常用的濾波算法有低通線性濾波、高斯平滑濾波及中值濾波，通過(guò)綜合對(duì)比三種算法的濾波效果，本研究中選用中值濾波器.巖芯灰度圖像經(jīng)中值濾波器進(jìn)行濾波處理之后，孔隙和巖石骨架之間的過(guò)渡變得自然，邊界也變得平滑，同時(shí)也盡可能地保留了圖像重要特征信息（圖1d）.但為了更好地區(qū)分及量化孔隙和骨架，還需采用圖像分割方法對(duì)灰度圖像進(jìn)行合理的二值劃分.圖像二值化的關(guān)鍵在于分割閾值的選取，鑒于本文用于微CT掃描的巖芯已知實(shí)測(cè)孔隙度，所以可采用基于巖芯孔隙度尋求到的最佳分割閾值來(lái)對(duì)圖像進(jìn)行分割.以實(shí)測(cè)孔隙度為約束尋求分割閾值k＊的公式如下：

圖1 微CT技術(shù)流程（a）微CT儀器；（b）巖樣；（c）CT切片；（d）濾波后切片；（e）二值化結(jié)果.Fig.1 The process of micro-CT technology（a）MicroXCT-400；（b）Rock sample；（c）One of CT slices；（d）The slice after filtering；（e）The result of binarization.

式中，巖芯孔隙度為φ，灰度閾值為k，圖像的最大、最小灰度值分別為Imax、Imin，灰度值為i的體素?cái)?shù)為p（i），灰度低于閾值的體素表征孔隙，其余代表骨架.以最終搜尋到的k＊作為分割閾值，得到分割后的二值圖像（圖1e），其中黑色為孔隙，白色為骨架.在此基礎(chǔ)上，還可根據(jù)實(shí)際需要，采用數(shù)學(xué)形態(tài)學(xué)算法對(duì)其作進(jìn)一步精細(xì)處理，即通過(guò)開(kāi)運(yùn)算移除孤立體素，通過(guò)閉運(yùn)算填充細(xì)小孔洞，連接鄰近體素.

2.3 3D體表面重建

理論上數(shù)字巖芯尺寸越大，就越能準(zhǔn)確表征巖石的微觀孔隙結(jié)構(gòu)和宏觀特性，然而數(shù)字巖芯尺寸越大，對(duì)計(jì)算機(jī)存儲(chǔ)和運(yùn)算能力要求就越高，因此折衷方案是選取代表元體積（REV），姜黎明等（2012）通過(guò)多次試驗(yàn)表明當(dāng)數(shù)字巖芯大小為200×200×200體素時(shí)，其物理性質(zhì)（比如孔隙度、彈性模量等）幾乎不再受尺寸的影響.在本文研究中，出于計(jì)算存儲(chǔ)和計(jì)算速度的考慮，選取代表元體積為200×200×200體素.

采用Marching Cube算法從圖像處理結(jié)果的REV三維數(shù)據(jù)體中提取表面的三角面片集，再用光照模型對(duì)三角面片進(jìn)行渲染，進(jìn)而形成巖芯的三維體表面圖像，至此三維數(shù)字巖芯建模工作完成（圖2）.

3 數(shù)字巖石物理實(shí)驗(yàn)

在數(shù)字巖芯的基礎(chǔ)上，通過(guò)各種形態(tài)學(xué)算法及數(shù)值模擬手段，可以統(tǒng)計(jì)、計(jì)算多種巖石物理參數(shù)，這即是所謂的數(shù)字巖石物理實(shí)驗(yàn).

3.1 孔隙結(jié)構(gòu)量化及表征

基于上述步驟所建數(shù)字巖芯的孔隙模型中（圖2c），大部分孔隙與孔隙之間接觸緊密，很難區(qū)分單個(gè)孔隙的邊界，這不利于后期定量統(tǒng)計(jì)孔隙體積分布及孔徑分布.為此，需要識(shí)別出每個(gè)孔隙的邊界，并對(duì)其進(jìn)行標(biāo)記.本文在研究中采用快速分水嶺算法進(jìn)行孔隙邊緣檢測(cè)，其基本原理是把圖像看作地學(xué)上的拓?fù)涞孛?，圖像上每一像素點(diǎn)的灰度值表示該點(diǎn)海拔高度，每一個(gè)局部極小值及其影響區(qū)域稱為集水盆地，集水盆地的邊界則形成分水嶺.通過(guò)該算法每個(gè)孔隙都能清楚地識(shí)別，類(lèi)似于都貼上了獨(dú)有的標(biāo)簽（圖3a），可以很方便地對(duì)號(hào)提取以進(jìn)行定量分析.一旦每個(gè)孔隙體積確定，可以統(tǒng)計(jì)出孔隙體積的分布直方圖（圖3c），還可以根據(jù)下面公式計(jì)算孔隙度：

式中，孔隙度φ為小數(shù)，Vpore為單個(gè)孔隙體積，單位pix3，Vvoxel為總體素的體積，單位pix3，pix是指一個(gè)像素，在本文研究中為2.1μm.

由表1可見(jiàn)，計(jì)算所得孔隙度略低于實(shí)測(cè)孔隙度，分析誤差來(lái)源，主要是圖像處理平滑造成，剔除掉的一部分小孔對(duì)計(jì)算孔隙度應(yīng)有所貢獻(xiàn).

圖2 數(shù)字巖芯模型（a）孔隙和骨架；（b）骨架（孔隙透明）；（c）孔隙（骨架透明）.Fig.2 Digital core model（a）Pore and frame；（b）Frame（with pore transparent）；（c）Pore（with frame transparent）.

圖3 孔隙結(jié)構(gòu)量化及表征（a）孔隙標(biāo)記圖；（b）孔隙網(wǎng)絡(luò)模型；（c）孔隙體積分布；（d）孔隙直徑分布.Fig.3 Quantification and characterization of pore structure（a）Label image of pore；（b）Pore network model；（c）The distribution of pore volume；（d）The distribution of pore diameter.

表1 孔隙度計(jì)算結(jié)果Table 1 The computation result of porosity

為了更加簡(jiǎn)明直觀地展示孔隙空間的拓?fù)浣Y(jié)構(gòu)，本文在數(shù)字巖芯的基礎(chǔ)上，采用形態(tài)學(xué)細(xì)化算法獲取孔隙空間中軸線，并將中軸線節(jié)點(diǎn)定義為孔隙，節(jié)點(diǎn)之間的連接線定義為喉道，由此建立了能夠簡(jiǎn)化表征孔隙空間拓?fù)浣Y(jié)構(gòu)的等價(jià)孔隙網(wǎng)絡(luò)模型（圖3b），圖中球體表征孔隙，管束表征喉道.球體體積與相應(yīng)位置的孔隙體積近似相等，每個(gè)孔隙的等效孔徑則可通過(guò)公式（3）確定，最終統(tǒng)計(jì)得到孔徑分布直方圖（圖3d）.

式中，等效孔隙直徑Deq單位為pix.

3.2 絕對(duì)滲透率數(shù)值模擬

巖石的絕對(duì)滲透率衡量的是飽和單相流體通過(guò)其孔隙空間的能力，這就要求巖石內(nèi)部必須存在相互連通的有效孔隙，才能提供相應(yīng)滲流路徑.因此，在絕對(duì)滲透率的數(shù)值模擬中，為保證數(shù)模能順利進(jìn)行并較快收斂，首先需對(duì)數(shù)字巖芯的孔隙空間進(jìn)行連通性測(cè)試，移除孤立“死孔”，然后再對(duì)孔隙空間進(jìn)行四面體網(wǎng)格剖分及優(yōu)化，最后通過(guò)有限元求解器實(shí)現(xiàn)數(shù)值模擬.本文采用Comsol軟件的不可壓縮Navier-Stokes方程模塊來(lái)完成孔隙空間的微流動(dòng)模擬，流體基本屬性按常態(tài)下水的參數(shù)賦值，模型中相對(duì)立的兩面分別作為速度入口及壓力出口邊界，其余流動(dòng)邊界及孔壁視為無(wú)滑移壁面（流速為0）.據(jù)此分別模擬了X、Y、Z三個(gè)方向的滲流特性，模擬得到三個(gè)方向的速度場(chǎng)分布及流線圖如圖4.

在數(shù)值模擬結(jié)果中，由出口或入口邊界上對(duì)流動(dòng)速度進(jìn)行積分，可以得到通過(guò)巖樣的體積流量，再代入達(dá)西定律公式（4）中即可求得絕對(duì)滲透率：

式中，流量Q單位cm3·s-1，巖芯截面積A單位cm2，巖芯長(zhǎng)度L單位cm，流體黏度μ單位為mPa·s，壓差ΔP單位MPa，計(jì)算所得滲透率K單位為μm2.三個(gè)方向的滲透率模擬結(jié)果見(jiàn)表2.分別計(jì)算三者的算術(shù)平均值、幾何平均值、調(diào)和平均值并與實(shí)驗(yàn)室氣測(cè)、液測(cè)結(jié)果進(jìn)行對(duì)比，發(fā)現(xiàn)計(jì)算結(jié)果均低于氣測(cè)值，且看不出有明顯的聯(lián)系，這是由于氣體滑脫效應(yīng)的存在，同一巖石的氣測(cè)滲透率為液測(cè)結(jié)果的3～5倍不等，而通過(guò)與液測(cè)結(jié)果的對(duì)比可看出：三個(gè)方向滲透率的幾何平均與實(shí)驗(yàn)結(jié)果較為接近.

圖4 速度場(chǎng)分布（顏色越亮，流速越大）Fig.4 The distribution of velocity field in the X，Y，Zdirections（the brighter the colors，the higher the velocity）

圖5 流線圖Fig.5 The streamline in the X，Y，Zdirections

表2 滲透率結(jié)果對(duì)比Table 2 The comparison of permeability result

3.3 彈性參數(shù)數(shù)值模擬

巖石的彈性參數(shù)（體積模量、剪切模量等）在地球物理勘探與測(cè)井領(lǐng)域發(fā)揮著重要作用.從結(jié)構(gòu)上看，巖石是由骨架和孔隙流體組成的復(fù)合介質(zhì)，巖石的彈性實(shí)則是各組分彈性性質(zhì)綜合而成的有效彈性.因此，多孔巖石的彈性參數(shù)不僅取決于固體骨架的彈性性質(zhì)，巖石中孔隙的大小、幾何形狀以及孔隙流體性質(zhì)都會(huì)對(duì)巖石總體彈性參數(shù)產(chǎn)生一定的影響，定量研究這之間的關(guān)系也一直是地球物理領(lǐng)域關(guān)注的熱點(diǎn)，對(duì)于油氣勘探開(kāi)發(fā)中儲(chǔ)層的預(yù)測(cè)具有重要指導(dǎo)意義.時(shí)至今日，國(guó)際上比較有代表性的經(jīng)典理論及經(jīng)驗(yàn)公式包括：微分有效介質(zhì)理論、Voigt-Reuss-Hill模型、Hashin-Shtrikman邊界方程、自洽理論、Gassmann方程、Wyllie公式、Raymer方程、百靈方程等.其中，流體替換Gassmann方程是研究孔隙飽和流體對(duì)巖石聲波速度影響最常用的理論，國(guó)內(nèi)外已有學(xué)者（Saenger，2008；姜黎明等，2012）基于三維數(shù)字巖芯分別通過(guò)有限元／旋轉(zhuǎn)-交錯(cuò)網(wǎng)格有限差分技術(shù)模擬了流體替換對(duì)巖石彈性性質(zhì)的影響，并與Gassmann方程計(jì)算結(jié)果進(jìn)行了對(duì)比驗(yàn)證，結(jié)果吻合度較高，證實(shí)了數(shù)值模擬復(fù)雜孔隙巖石有效彈性參數(shù)的可靠性.然而該方程的基本假設(shè)之一是孔隙空間飽和無(wú)摩擦流體，對(duì)于孔隙充填物為黏滯性稠油、瀝青的情況不再適用.為此，Ciz和Shapiro（2007）提出了針對(duì)高黏度物質(zhì)或固相充填孔隙的近似Gassmann公式：

式中，Ksat、KB、KA、Kdry分別指有效體積模量、基質(zhì)礦物體積模量、孔隙充填相體積模量、干巖石體積模量，單位均為GPa，孔隙度φ為小數(shù).公式（5）是一個(gè)近似固相替換方程，其準(zhǔn)確度依賴于巖石的孔隙結(jié)構(gòu)，為評(píng)價(jià)其有效性，采用傳統(tǒng)物理實(shí)驗(yàn)方法無(wú)法實(shí)現(xiàn)，因此，本文采用Comsol軟件的結(jié)構(gòu)力學(xué)模塊開(kāi)展數(shù)值模擬研究.基于三維數(shù)字巖芯模型用有限元求解孔隙尺度的線彈性方程，其理論基礎(chǔ)是最小勢(shì)能原理，對(duì)于給定的數(shù)字巖芯，施加一個(gè)宏觀的體積應(yīng)變，在周期性邊界條件的控制下，利用共軛梯度法通過(guò)把體系的彈性勢(shì)能最小化來(lái)求取由這個(gè)外加應(yīng)變引起的平均應(yīng)力，進(jìn)而求得巖石的有效彈性參數(shù).

數(shù)值模擬中，礦物基質(zhì)體積模量KB根據(jù)Hill（1963）的研究取為36GPa，體積應(yīng)變賦為0.001，基于簡(jiǎn)化后的網(wǎng)格（圖6），模擬了六種不同孔隙充填相（KA分別取0.1GPa、5.1GPa、10.1GPa、15.1 GPa、20.1GPa、25.1GPa）的巖石有效體積模量.其中，公式（5）中干巖石體積模量Kdry不同于氣體飽和巖石的體積模量，它對(duì)應(yīng)于孔隙相體積模量為0的情況，而氣體具有不可忽略的體積模量，因此通過(guò)實(shí)驗(yàn)獲得理想的Kdry比較困難，一般通過(guò)經(jīng)驗(yàn)公式求得.而本文采用數(shù)值模擬計(jì)算Kdry時(shí)，為保障數(shù)模運(yùn)算能夠完成且又不影響結(jié)果，取孔隙充填相KA為0.0001GPa，計(jì)算得到Kdry約為21GPa.圖7為KA取25.1GPa時(shí)模擬結(jié)果應(yīng)力分布圖，由圖可見(jiàn)局部高應(yīng)力區(qū)出現(xiàn)在孔隙充填相的邊緣附近.

將數(shù)值模擬所得到的巖石有效體積模量與公式（5）計(jì)算出的結(jié)果進(jìn)行對(duì)比（圖8），發(fā)現(xiàn)兩者基本吻合，再次證明：無(wú)論是孔隙中飽和流體還是固相充填，基于巖石的微觀結(jié)構(gòu)用有限元的方法模擬復(fù)合巖石的彈性性質(zhì)都是可行且可靠的.

4 結(jié)論與展望

（1）利用微CT掃描結(jié)合先進(jìn)的三維可視化軟件Avizo建立的數(shù)字巖芯可以精細(xì)地捕捉巖石真實(shí)孔隙結(jié)構(gòu)特征，基于數(shù)字巖芯模型開(kāi)展數(shù)字巖石物理實(shí)驗(yàn)，可以準(zhǔn)確地獲取孔隙度并統(tǒng)計(jì)得到孔隙體積分布、孔徑分布特征，以及能夠簡(jiǎn)化表征孔隙空間拓?fù)浣Y(jié)構(gòu)的孔隙網(wǎng)絡(luò)模型；

（2）將Avizo與多場(chǎng)耦合有限元軟件Comsol對(duì)接，可以形象直觀地模擬流體在真實(shí)孔隙空間的滲流過(guò)程，并計(jì)算獲得絕對(duì)滲透率；

（3）運(yùn)用Comsol軟件模擬計(jì)算了固相充填孔隙情況下巖石的有效彈性參數(shù)，數(shù)模結(jié)果與固相替換近似Gassmann方程能夠很好的相互驗(yàn)證.

圖6 簡(jiǎn)化網(wǎng)格Fig.6 The simplified grid

圖7 應(yīng)力分布圖（KA＝25.1GPa）Fig.7 The stress distribution（KA＝25.1GPa）

圖8 有效體積模量結(jié)果對(duì)比Fig.8 The results contrast of effective bulk modulus

隨著計(jì)算機(jī)技術(shù)的發(fā)展，數(shù)字巖石物理必將成為一項(xiàng)重要的技術(shù)手段參與到油氣田的開(kāi)發(fā)建設(shè).本文的研究提供了數(shù)字巖石物理一套新的方法體系，同時(shí)也為進(jìn)一步研究低滲致密砂巖的相滲、巖電、聲波傳播及其相互聯(lián)系奠定了基礎(chǔ).

Arns C H，Bauget F，Ghous A，etal.2004a.Digital Core Laboratory：petrophysical analysis from 3Dimaging of reservoir core fragments.Petrophysics，46（4）：260-277.

Arns C H，Knackstedt M A，Pinczewski W V，etal.2004b.Virtual permeametry on microtomographic images.Journal of Petroleum Science and Engineering，45（1-2）：41-46.

Ciz R，Shapiro S A.2007.Generalization of Gassmann equations for porous media saturated with a solid material.Geophysics，72（6）：A75-A79.

Derzhi N，Dvorkin J，Diaz E，etal.2010.Comparison of traditional and digital rock physics techniques to determine the elastic core parameters in cretaceous formations，Abu Dhabi.SPE138586.

Rezaei-Gomari S，Berg F，Mock A，etal.2011.Electrical and petrophysical properties of siliciclastic reservoir rocks from pore-scale modeling.SCA International Symposium，Texas，18-21.

Hill R.1963.Elastic properties of reinforced solids：some theoretical principles.Journal of the Mechanics and Physics of Solids，11（5）：357-372.

Hu D，Touati M，Blunt M J.2007.Pore network modeling：analysis of pore size distribution of Arabian core samples.SPE105156.

Hu D.2007.Micro-CT imaging and pore network extraction［Ph.D.thesis］.London：Imperial College.

Jiang L M，Sun J M，Liu X F，etal.2012.Numerical study of the effect of natural gas saturation on the reservoir rocks′elastic parameters.Well Logging Technology（in Chinese），36（3）：239-243.

Keehm Y.2003.Computational rock physics：transport properties in porous media and applications［Ph.D.thesis］.Stanford：Stanford University.

Knackstedt M A，Arns C H，Pinczewski W V，etal.2002.Computation of linear elastic properties from microtomographic images：methodology and match to theory and experiment.Geophysics，67（5）：1396-1405.

Liu X F，Sun J M，Wang H T.2009a.Numerical simulation of rock electrical properties based on digital cores.Applied Geophysics，6（1）：1-7.

Liu X F，Sun J M，Wang H T.2009b.Reconstruction of 3-D digital cores using a hybrid method.Applied Geophysics，6（2）：105-112.

Okabe H，Blunt M J.2005.Pore space reconstruction using multiple-point statistics.Journal of Petroleum Science and Engineering，46（1-2）：121-137.

Φren P E，Bakke S.2002.Process based reconstruction of sandstones and prediction of transport properties.Transport in Porous Media，46（2-3）：311-343.

Saenger E H，Bohlen T.2004.Finite-difference modeling of viscoelastic and anisotropic wave propagation using the rotated staggered grid.Geophysics，69（2）：583-591.

Saenger E H.2008.Numerical methods to determine effective elastic properties.International Journal of Engineering Science，46（6）：598-605.

Sain R.2010.Numerical simulation of pore-scale heterogeneity and its effects on elastic，electrical and transport properties［Ph.D.thesis］.Stanford：Stanford University.

Su N，Duan Y G，Chen W，etal.2010.Three-dimensional reconstruction of micro pore structure.International Conference on Computational and Information Sciences（ICCIS）.120-123.

Tao G，Yue W Z，Xie R H，etal.2005.A new method for theoretical modeling and numerical experiments on petrophysical studies.Progress in Geophysics（in Chinese），20（1）：4-11.

Wang K W，Sun J M，Geng S C，etal.2006.Percolation network study of shale effects on rock electrical properties under different salinity.Chinese J.Geophys.（in Chinese），49（6）：1867-1872.

Yue W Z，Tao G，Zhu K Q.2004.Simulation of electrical properties of rock saturated with multi-phase fluids using lattice gas automation.Chinese J.Geophys.（in Chinese），47（5）：905-910.

Zhao X C，Yao J，Yi Y J.2007.A new stochastic method of reconstructing porous media.Transport in Porous Media，69（1）：1-11.

附中文參考文獻(xiàn)

姜黎明，孫建孟，劉學(xué)鋒等.2012.天然氣飽和度對(duì)巖石彈性參數(shù)影響的數(shù)值研究.測(cè)井技術(shù)，36（3）：239-243.

陶果，岳文正，謝然紅等.2005.巖石物理的理論模擬和數(shù)值實(shí)驗(yàn)新方法.地球物理學(xué)進(jìn)展，20（1）：4-11.

王克文，孫建孟，耿生臣等.2006.不同礦化度下泥質(zhì)對(duì)巖石電性影響的逾滲網(wǎng)絡(luò)研究.地球物理學(xué)報(bào)，49（6）：1867-1872.

岳文正，陶果，朱克勤.2004.飽和多相流體巖石電性的格子氣模擬.地球物理學(xué)報(bào)，47（5）：905-910.