基于微CT技術(shù)的砂巖數(shù)字巖石物理實驗

劉向君，朱洪林，梁利喜

油氣藏地質(zhì)及開發(fā)工程國家重點實驗室，西南石油大學(xué)，成都 610500

1 引言

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

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

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

2.1 微CT成像

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

2.2 圖像處理

微CT掃描獲得的巖芯灰度圖像中存在各種類型的系統(tǒng)噪聲，降低圖像質(zhì)量的同時也不利于后續(xù)的定量分析，因此圖像處理第一步是通過濾波算法增強信噪比.針對三維圖像，比較常用的濾波算法有低通線性濾波、高斯平滑濾波及中值濾波，通過綜合對比三種算法的濾波效果，本研究中選用中值濾波器.巖芯灰度圖像經(jīng)中值濾波器進行濾波處理之后，孔隙和巖石骨架之間的過渡變得自然，邊界也變得平滑，同時也盡可能地保留了圖像重要特征信息（圖1d）.但為了更好地區(qū)分及量化孔隙和骨架，還需采用圖像分割方法對灰度圖像進行合理的二值劃分.圖像二值化的關(guān)鍵在于分割閾值的選取，鑒于本文用于微CT掃描的巖芯已知實測孔隙度，所以可采用基于巖芯孔隙度尋求到的最佳分割閾值來對圖像進行分割.以實測孔隙度為約束尋求分割閾值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的體素數(shù)為p（i），灰度低于閾值的體素表征孔隙，其余代表骨架.以最終搜尋到的k＊作為分割閾值，得到分割后的二值圖像（圖1e），其中黑色為孔隙，白色為骨架.在此基礎(chǔ)上，還可根據(jù)實際需要，采用數(shù)學(xué)形態(tài)學(xué)算法對其作進一步精細處理，即通過開運算移除孤立體素，通過閉運算填充細小孔洞，連接鄰近體素.

2.3 3D體表面重建

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

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

3 數(shù)字巖石物理實驗

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

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

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

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

由表1可見，計算所得孔隙度略低于實測孔隙度，分析誤差來源，主要是圖像處理平滑造成，剔除掉的一部分小孔對計算孔隙度應(yīng)有所貢獻.

圖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 孔隙度計算結(jié)果Table 1 The computation result of porosity

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

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

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

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

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

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

圖4 速度場分布（顏色越亮，流速越大）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é)果對比Table 2 The comparison of permeability result

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

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

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

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

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

4 結(jié)論與展望

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

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

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

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

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

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

隨著計算機技術(shù)的發(fā)展，數(shù)字巖石物理必將成為一項重要的技術(shù)手段參與到油氣田的開發(fā)建設(shè).本文的研究提供了數(shù)字巖石物理一套新的方法體系，同時也為進一步研究低滲致密砂巖的相滲、巖電、聲波傳播及其相互聯(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.

附中文參考文獻

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

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

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

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