林承焰，王 楊，楊 山，任麗華，由春梅，吳松濤，吳玉其，張依旻

1.中國石油大學(華東)地球科學與技術(shù)學院，山東 青島 266580 2.山東省油藏地質(zhì)重點實驗室，山東 青島 266580 3.中國石油西南油氣田分公司勘探開發(fā)研究院，成都 610051 4.大慶油田有限責任公司勘探開發(fā)研究院，黑龍江 大慶 163000 5.中石油勘探開發(fā)研究院，北京 100083

0 引言

在石油工業(yè)的推動下，數(shù)字巖心技術(shù)在油藏工程以及地學領(lǐng)域得到廣泛應(yīng)用和推廣，其中包括數(shù)字巖心建模及流動模擬技術(shù)等。目前國內(nèi)外對于數(shù)字巖心建模方法的研究主要有兩種：一是采用X射線CT(computed tomography)掃描巖心，利用圖像處理算法將二維圖像構(gòu)建為三維數(shù)字巖心；二是基于巖石二維信息(如鑄體薄片和巖石粒度分析資料等)，利用建模算法構(gòu)建三維數(shù)字巖心[1]。美國學者Rosenberg等[2]利用CT掃描技術(shù)建立了楓丹白露砂巖的三維數(shù)字巖心，Arns[3]也通過三維CT掃描得到了孔隙度分別為8%、13%、15%和22%四塊楓丹白露砂巖的三維數(shù)字巖心，掃描分辨率為5.68 μm/像素。相比于CT掃描建立數(shù)字巖心，基于薄片分析的圖像建模方法只需要一定數(shù)量的巖石切片掃描圖像，獲取較為方便、經(jīng)濟，該方法被廣泛采用[4]。圖像建模算法主要包括隨機法和過程法，其中隨機法主要包括高斯場法、模擬退火法和序貫指示法等，主要是根據(jù)數(shù)學算法在巖石二維圖像統(tǒng)計特征的約束下構(gòu)建三維數(shù)字巖心，使其與巖石二維圖像具有相似的統(tǒng)計特征[5]。Joshi[6]首次利用高斯場法建立三維數(shù)字巖心。Hazlett[7]首次利用模擬退火法建立三維數(shù)字巖心。Bakke等[8]結(jié)合巖石的顆粒粒徑分布，通過對沉積類巖石地質(zhì)形成過程(包括沉積、壓實和成巖作用)的模擬來構(gòu)建數(shù)字巖心。Arns等[9]通過建立的數(shù)字模型計算得到了滲透率等數(shù)值，并與實驗室測量所得的結(jié)果進行比較來驗證模型的準確性。劉洋[10]通過對巖心進行切片分析，獲取過程法構(gòu)建數(shù)字巖心所需要的相關(guān)儲層參數(shù)，利用過程法的建模算法，實現(xiàn)了巖心的數(shù)字可視化。Youssef等[11]應(yīng)用圖像處理算法對模型孔喉進行了兩相分割和孔喉參數(shù)統(tǒng)計。Ryazanov等[12]運用微觀孔喉網(wǎng)絡(luò)模型模擬了不同潤濕性條件下的剩余油類型及其體積含量。劉向君等[13]通過采用COMSOL軟件的不可壓縮方程模塊完成了孔隙空間的微流動模擬。李易霖等[14]通過Avizo建立了致密砂巖多尺度數(shù)字巖心模型，并對儲層微觀孔隙特征進行了定量表征。20世紀80年代中期發(fā)展形成的格子玻爾茲曼方法(lattice Boltzmann method, LBM)，來源于完全離散格子氣動機算法，目前主要模型有：RK格子波爾茲曼模型、SC模型、自由能方法模型和基于連續(xù)性波爾茲曼方程的熱力學一致性多相流體理論模型。Pan等[15]首次使用格子玻爾茲曼法和SC模型模擬兩相流體在孔隙吼道中的流動，實現(xiàn)了在介質(zhì)中微觀油驅(qū)水的數(shù)字化過程。

本文依據(jù)二維CT資料，通過Avizo軟件將CT生成的切片數(shù)據(jù)進行精確的三維體重構(gòu)。建模過程中，采用了更為準確的非局部均值濾波算法(NLM filter)進行濾波降噪處理，通過對骨架和孔喉分割，認識儲層空間結(jié)構(gòu)，分析孔喉結(jié)構(gòu)特征。流動模擬研究采用等壓面的假設(shè)，在微觀尺度上建立了微觀流體流動模擬模型，以期為數(shù)字巖心動態(tài)模擬提供一種新方法。

1 數(shù)字巖心模型建立

CT掃描建立數(shù)字巖心主要是依據(jù)巖心內(nèi)部各部分組成不同，其對X射線的吸收強度就不同，因此通過計算機檢測射線經(jīng)過樣品吸收后的強度，可以計算得到這個層面上X射線的吸收部分，根據(jù)吸收的射線強度就可以轉(zhuǎn)化成對應(yīng)灰度值以表征不同的結(jié)構(gòu)組分。主要儀器包括X射線源，固定巖樣的載物臺和收集掃描后射線的X射線檢測器(圖1)。X射線源向載物臺上的物品發(fā)出X射線，與物品發(fā)生一系列作用后被接收器接收，接收器處理并將其轉(zhuǎn)化為電信號后返回給計算機，記錄不同位置的灰度信息。測量樣品的密度越大，對X射線的吸收就越多。利用計算機重建可得到三維灰度圖像，該圖像中的灰度值與被測樣品的組分密度值相互對應(yīng)[16]。

圖1 X射線CT掃描原理圖Fig.1 Schematic diagram of X-ray computed tomography

1.1 切片獲取與模型構(gòu)建

X射線CT掃描是建立三維數(shù)字巖心最常用的方法，但是在實際應(yīng)用中也存在一些問題[17-18]：由于同一塊巖心在不同分辨率下得到的孔隙網(wǎng)絡(luò)模型基本結(jié)構(gòu)參數(shù)和滲流特征有明顯差異，如果選取的分辨率過低，就難以識別出巖心中的微細孔喉，使得建模得到的孔隙度較巖心的實際孔隙度?。蝗绻x取的分辨率過高，掃描所得的數(shù)據(jù)量偏大，因此該方法對計算機等軟硬件條件要求較高。本文選取H152工區(qū)的細粉砂巖巖樣(表1)，根據(jù)儲層物性特點，選取分辨率為1.49 μm/像素，掃描得到巖樣的二維圖像，通過對二維切片進行順序調(diào)整、角度矯正和黑點消除，構(gòu)建三維數(shù)字巖心模型(圖2)。

1.2 平滑濾波

巖心掃描后，由于噪聲等的干擾，會出現(xiàn)灰度值不同的孤立像素點。對于骨架及孔隙來說，即使是相同的組分，也可能會因為密度及含流體性質(zhì)的細微差異而呈現(xiàn)出不同的灰度值；因此需要采用適當?shù)臑V波算法，去除噪聲干擾，同時使相同組分灰度值分布更為集中，不同組分灰度值區(qū)分更加明顯，便于后期分割閾值的選取。不同的濾波算法中，中值濾波為高強度的平滑濾波，但在濾波過程中沒有考慮到樣品的結(jié)構(gòu)信息，僅僅是對選擇區(qū)域內(nèi)所有樣品的灰度值進行大小排序，然后選擇中間值直接賦值。而非局部均值濾波(NLM filter)不同于中值濾波的原理，其基本思想是：考慮圖像的自相似性，當前像素的預(yù)測值可由圖像中所有與它結(jié)構(gòu)相似的像素加權(quán)平均得到[19]。其具體方法為：在選擇區(qū)域內(nèi)對比目標和其他部分的結(jié)構(gòu)信息，根據(jù)結(jié)構(gòu)上的相似性賦值權(quán)重，相似性越高，在最終灰度值的計算中這部分的權(quán)重就越大；反之越小。這樣就能充分考慮樣品結(jié)構(gòu)上的相似性，保留樣品的結(jié)構(gòu)信息，便于不同組分的分割。在實際研究工作中，選取的濾波搜索區(qū)域為21像素，目標區(qū)域為5像素，通過濾波處理得到灰度值分布直方圖(圖3)和二維橫截面灰度圖像濾波效果圖(圖4)。由于中值濾波的高強度平滑作用，因此濾波后的圖像高度平滑但是清晰度大大降低(圖4b)；而非局部均值濾波處理后，相同組分的單峰形態(tài)更加集中，灰度值雙峰對比更加明顯(圖3)，不同灰度值區(qū)域的交界處圖像有明顯的平滑處理，相同組分內(nèi)部灰度值保真度較高(圖4c)，便于后期相態(tài)的識別與劃分。

表1 CT樣品參數(shù)表

a.2D切片；b.3D模型。圖2 三維數(shù)字巖心模型構(gòu)建Fig.2 Construction of 3D digital model

圖3 NLM濾波前(a)后(b)灰度值分布直方圖Fig.3 Histogram of gray value distribution before NLM filtering (a) and after (b)

a.原始圖像；b.中值濾波圖像；c. NLM 濾波圖像。圖4 NLM濾波和中值濾波效果對比圖Fig.4 Contrast figure of NLM filter and median filter

1.3 孔喉分割

微觀模型的準確性關(guān)鍵取決于骨架和孔喉的分割。數(shù)字巖心建模中孔喉的識別方法主要有兩種：一種是根據(jù)已知巖樣的實際孔隙度大小進行孔喉劃分，使劃分后的模型孔隙度與實際相同；另一種是依據(jù)最大類間距方法進行分割閾值的選取，通過設(shè)置不同的閾值，使得組內(nèi)方差與組間方差的比值最小。本文采用第二種方法，應(yīng)用Avizo軟件進行方差計算和兩相分割，分割效果如圖5所示。

經(jīng)過閾值分割后，大部分的孔隙與喉道都可以區(qū)分出來，但是由于掃描分辨率與巖樣孔喉結(jié)構(gòu)兩方面的原因，當相鄰的顆粒之間的間隙非常小，即微孔喉尺寸小于掃描分辨率時，部分微小的孔喉被忽略，顆粒被處理成直接接觸的形式(圖5b)；因此需要應(yīng)用分離算法對顆粒進行分割并做一定的優(yōu)化。分離算法是分水嶺算法、距離轉(zhuǎn)換算法和數(shù)值重建算法的綜合體現(xiàn)，通過計算灰度梯度模量對顆粒進行分割，確定顆粒最佳輪廓，顯示微小孔喉分布，識別結(jié)果如圖5c所示。

分離算法的處理可以將巖樣中的顆粒分離，得到部分微孔微喉；但是算法的處理過程相對機械，分割產(chǎn)生的孔喉需要進行模型的優(yōu)化運算。模型的優(yōu)化主要是對模型進行開運算(opening)和閉運算(closing)。開運算的思路是用運算元沿著圖像的邊緣進行腐蝕運算，將圖像的狹窄部分剔除，然后再用此運算元沿著模型的邊緣進行膨脹運算。經(jīng)過開運算處理后，可以消除小物體、在纖細點處分離物體、平滑較大物體的邊界，同時并不明顯改變其面積。閉運算的處理與開運算的處理過程相反，即先用運算元進行膨脹運算，然后再用此運算元進行腐蝕邊緣，最終可以起到充填細小空間、平滑邊緣、連接孤立圖像的作用。通過對比可以發(fā)現(xiàn)，模型優(yōu)化后孤立的像素點明顯減少，骨架的邊緣也明顯平滑(圖6)。

a.原始灰度圖像；b.兩相分割圖像(白色為巖石骨架，黑色為孔喉，顆粒直接線接觸，孔隙不連通)；c.微孔喉識別圖像。圖5 二值分割效果圖Fig.5 Contrast figure of interactive threshold

圖6 開運算閉運算前(a)后(b)對比圖Fig.6 Contrast figure before opening and closing (a) and after (b)

1.4 體積元提取

由于計算機計算性能的限制，微觀特征分析需要在整個模型中選取一定大小的體積元，以此表征整個模型。本文體積元(REV)的選取以孔隙度為約束條件進行。即在模型中選取3個像素點為中心，分別構(gòu)建棱長為50像素的立方體并計算孔隙度，然后依次以50像素增大立方體的棱長并計算孔隙度，做出3個立方體孔隙度-棱長散點曲線，如圖7中3條曲線所示。由圖7可知，當立方體棱長大于300像素時，不同位置立方體的自相關(guān)函數(shù)趨于一致。因此，該樣品的REV尺寸為300像素，即數(shù)字巖心的體積元棱長為300像素。多次選取體積元進行的計算研究表明，當數(shù)字巖心體積元棱長為300像素時，其物理性質(zhì)幾乎不受尺寸的影響，建立REV模型如圖8所示。

圖7 孔隙度與立方體邊長的關(guān)系曲線Fig.7 Relation curve of porosity and cube

2 數(shù)字巖心模型分析

2.1 靜態(tài)參數(shù)分析

1)顆粒直徑

利用標定分析算法可計算模型內(nèi)部顆粒的體積分布參數(shù)V3d，然后得到等效直徑公式：

(1)

利用公式(1)可計算出顆粒等效直徑分布情況(表2)，并對每個顆粒的直徑進行統(tǒng)計，建立顆粒等效直徑分布直方圖(圖9)。

根據(jù)表2和圖9可知，粒徑主要分布區(qū)間為10～70像素。已知掃描該樣品所用分辨率為1.49 μm/像素，則數(shù)字巖心模型粒徑主要分布在14.9～104.3 μm；與研究區(qū)細粉砂巖的粒度主峰基本一致。數(shù)字巖心模型中最大顆粒直徑為129.2像素，約190 μm，這是由于極少部分的微細孔喉未能分割，造成顆粒之間為直接接觸，導致粒徑明顯偏大。

2)孔隙度

統(tǒng)計閾值分割后模型中骨架部分和孔隙部分的像素體積，結(jié)果如表3所示?？紫抖葹榭紫恫糠值南袼卦谡麄€模型中所占的體積比，約為15.54%。而研究區(qū)實際孔隙度為22.5%。計算孔隙度明顯小于實驗室測試值，其主要原因是部分微孔喉尺寸明顯小于分辨率1.49 μm，因此未能在切片上明顯識別，最終導致孔隙度偏小。

a.體積元渲染(綠色邊框為鏤空部分邊界)；b.體積元正交視圖(白色為巖石骨架，黑色為孔喉)。圖8 體積元三維渲染及正交視圖Fig.8 Volume rendering and ortho view of REV

均值/像素最小值/像素最大值/像素中值/像素55.034.66129.2034.50

圖9 顆粒等效直徑分布直方圖Fig.9 Histogram of diameters of labels’ Deq

模型總體積/(106μm3)骨架體積/(106μm3)孔喉體積/(106μm3)孔隙度/%8.947.551.3915.54

3)孔喉直徑分布

選取模型內(nèi)部的孔喉，計算得到孔隙體積V3d的分布情況，依據(jù)公式(1)計算內(nèi)徑大小，得到其孔喉直徑分布直方圖(圖10)。研究區(qū)孔喉粒徑較小，集中分布在20像素以下，大孔發(fā)育較少，儲層以中孔低滲為特征。

圖10 孔喉等效直徑分布直方圖Fig.10 Histogram of diameters of pores’ Deq

4)孔喉連通性與曲折度

孔喉連通性和孔喉曲折度是影響流體流動的重要參數(shù)，一般來說連通性越好，曲折度越低，流體的流動越容易。曲折度算法原理為計算孔隙邊界在某一軸向上不同切面之間的質(zhì)點距離，z(n)面到z(0)面距離疊加之和與H的比值即為曲折度(圖11)。連通性計算是依賴于滲流的一種算法，通過沿著某一軸向進行滲流，得到連通的體積與全體積的比值，即為巖樣的連通性。

圖11 曲折度計算原理Fig.11 Principle of centroid path of tortuosity

經(jīng)計算：模型的孔喉連通性為92.16%，連通性較好；曲折度為2.62，即孔喉邊界的曲線長度是對應(yīng)節(jié)點間直線長度的2.62倍，孔喉發(fā)育相對復雜，由此產(chǎn)生的流動阻力較大。根據(jù)模型連通性測試結(jié)果，依據(jù)數(shù)學邏輯算法，可得到模型內(nèi)部的孤立孔喉分布情況，結(jié)果如圖12所示。

5)孔喉網(wǎng)絡(luò)骨架提取

根據(jù)燃燒算法，用孔隙中心與喉道中軸線分別代表孔隙和喉道，將復雜不規(guī)則的孔喉空間用相互連接的中軸線骨架模型來表征(圖13a)。盡管中軸線骨架無法完整地刻畫孔隙類型和尺寸，但是它準確地描述了孔喉網(wǎng)絡(luò)的結(jié)構(gòu)特征[20]。根據(jù)孔喉直徑計算結(jié)果，將等效直徑值賦值到骨架模型中，可得到模型的孔喉結(jié)構(gòu)填充模型(圖13b)，該模型中不同顏色展示了模型內(nèi)部孔喉發(fā)育特征：模型內(nèi)部多發(fā)育小孔細喉，很少存在大的連通孔隙。圖13c展示了某一切面上灰度值圖像和孔喉的對應(yīng)特征。根據(jù)圖13b、c可得出：模型內(nèi)部孔喉不均一性較為突出，孔喉大小分布不均，孔喉結(jié)構(gòu)非均質(zhì)性明顯。

a.全孔喉；b.連通孔喉；c.孤立孔喉。圖12 孔喉連通性測試結(jié)果Fig.12 Results of connectivity

a.孔喉網(wǎng)絡(luò)骨架模型；b.孔喉結(jié)構(gòu)填充模型；c.切片對比。圖13 孔喉模型Fig.13 Pore-throat model

多孔介質(zhì)中孔隙和喉道的連接關(guān)系反映了多孔介質(zhì)內(nèi)部的連通性。一個孔隙所連接的喉道數(shù)量稱為配位數(shù)，配位數(shù)值的大小反映了網(wǎng)絡(luò)模型的連通性好壞，同時配位數(shù)也是影響巖石滲透率的重要因素[21]。因此，根據(jù)建立的孔喉網(wǎng)絡(luò)骨架模型可分析模型的孔喉配位數(shù)并建立頻率分布直方圖，結(jié)果如圖14所示。由圖14可知：模型中存在一定數(shù)量的配位數(shù)為1的孔隙，剩余油在這些孔隙中富集；除此之外，配位數(shù)主要為2～6，相對較低，孔喉連通性較差，模型滲透性較低，研究區(qū)儲層具有低滲透特征。

圖14 配位數(shù)頻率分布直方圖Fig.14 Frequency distribution histogram of coordination number

2.2 流動特征分析

1)絕對滲透率計算

根據(jù)建立的模型，利用軟件的有限元模擬模塊分別在x,y,z3個方向計算模型的滲透率數(shù)值，計算時設(shè)置進口壓力為130 kPa，出口壓力為100 kPa，流體黏度為0.001 Pa·s，得到x,y,z方向的滲透率分別為：33.7×10-3、41.7×10-3、23.8×10-3μm2。滲透率較小，屬于低滲透儲層，且不同方向滲透率值略有差異。實驗室測量該井段的水平滲透率值為19.28×10-3μm2，由于部分連通微孔喉未能識別出來，滲流阻力明顯變小，計算值略大于實驗室測量值。

2)速度場流線表征

根據(jù)計算的滲透率值，以流線的形式表征相態(tài)在模型內(nèi)部的流動路徑，左端為入口，右端為出口(圖15)。根據(jù)流線分布情況，模型內(nèi)部主流線區(qū)域為大孔喉連通部分，流線較為密集，模型角隅部分流線稀疏，存在大量的未波及區(qū)，是剩余油潛在富集區(qū)域。

3)流動模擬

在多孔介質(zhì)流動中，當介質(zhì)內(nèi)部彼此連通、滲流過程發(fā)生得非常緩慢時，即在準靜止條件下，每一個等壓面上其液-液界面和液-固界面都是平衡的，在不考慮重力和動力情況下，起作用的力就是毛管壓力[22]。這樣，在任一時刻，流體內(nèi)部的壓力會彼此傳遞，液面最前端壓力值保持相等，即在理想條件下流體液面形態(tài)、移動方向與等壓面保持一致。因此，壓降與相態(tài)流動液面等效，即液面前緣為等壓面。通過等壓面的變化來表征模型中相態(tài)的流動，并根據(jù)壓降的變化來計算流體飽和度大小。根據(jù)Stokes方程，設(shè)定流體為不可壓縮的牛頓流體，液固邊界無滑移，在準靜止條件下模擬層流滲流過程，結(jié)果如圖16所示。

a. x軸；b. y軸；c. z軸。藍色表示模型中連通的孔喉空間，紅色流線表示速度場分布。圖15 速度場流線表征Fig.15 Streamlines of velocity in the x, y, z axis

a. Sw=0；b. Sw=16.7%；c. Sw=78%。Sw為含水飽和度。紅色表示孔喉空間，藍色表示水相。圖16 流動模擬Fig.16 Percolation simulation

由圖16可知，在一端注水的情況下：由于孔喉大小的不同，液面會受到毛管壓力的作用，大孔喉連通區(qū)域液面出現(xiàn)突進現(xiàn)象，在三維上呈現(xiàn)為不平行的液面特征(圖16b)；隨著持續(xù)注水，水占據(jù)大部分的孔喉空間，但是部分微孔喉和孤立孔喉未被波及(圖16c)。該方法可簡化模型內(nèi)部流動過程，初步建立毛管壓力作用下的滲流模型，模型中的未波及區(qū)域為剩余油潛在富集區(qū)域。

3 結(jié)論與認識

1)通過對優(yōu)選的H152樣品進行X射線CT掃描，可以構(gòu)建三維數(shù)字巖心模型，能夠模擬孔喉網(wǎng)絡(luò)結(jié)構(gòu)，建立了微觀尺度儲層模型，定量表征了粒徑、孔喉半徑、配位數(shù)、孔隙度、滲透率等動靜態(tài)參數(shù)及低滲透儲層特征。

2)結(jié)合研究區(qū)儲層樣品孔滲實際資料，對所建數(shù)字巖心模型進行驗證和校對，能夠為儲層流動模擬提供可靠的地質(zhì)模型。

3)基于等壓面的微觀滲流模擬初步揭示了低滲透儲層流體滲流特征，提供了一種流動模擬新方法，但是由于只考慮了毛管壓力因素，還需進一步分析孔喉的潤濕性、流體重力等其他多種參數(shù)的影響。

[1] 劉學鋒. 基于數(shù)字巖心的巖石聲電特性微觀數(shù)值模擬研究[D].青島：中國石油大學(華東),2010.

Liu Xuefeng. Numerical Simulation of Elastic and Electrical Properties of Rock Based on Digital Cores [D]. Qingdao: China University of Petroleum (East China), 2010.

[2] Rosenberg E, Lynch J, Guéroult P, et al. High Re-solution 3D Reconstructions of Rocks and Composites[J]. Oil & Gas Science & Technology, 1999, 54(4):497-511.

[3] Arns C. H. The Influence of Morphology on Physical Properties of Reservoir Rocks [D]. Sydney: The University of New South Wales, 2002.

[4] 姚軍,趙秀才,衣艷靜,等. 數(shù)字巖心技術(shù)現(xiàn)狀及展望[J]. 油氣地質(zhì)與采收率,2005,12(6)：52-54.

Yao Jun, Zhao Xiucai, Yi Yanjing, et al. Status and Prospect on Digital Core Technology [J]. Petroleum Geology and Recovery Efficiency, 2005,12(6):52-54.

[5] 謝一婷. 基于數(shù)字巖心技術(shù)的低滲油藏滲流機理及合理井距研究[D].成都:西南石油大學,2013.

Xie Yiting. The Research of Seepage Flow in Low Permeability Reservoir and the Reasonable Well Spacing Based on Digital Core Technology [D]. Chengdu: Southwest Petroleum University, 2013.

[6]Joshi M. A Class of Stochastic Models for Porous Media[D]. Lawrence Kansas: University of Kansas, 1974.

[7] Hazlett R D. Statistical Characterization and Stochastic Modeling of Pore Networks in Relation to Fluid Flow [J]. Mathematical Geosciences, 1997, 29(6):801-822.

[8]Bakke S, Ren P E. 3-D Pore-Scale Modelling of Sandstones and Flow Simulations in the Pore Networks[J]. Spe Journal, 1997, 2(2):136-149.

[9] Arns C H, Averdunk H, Bauget F, et al. Digital Core Laboratory: Reservoir Core Analysis from 3D Images[C]// The 6th North America Rock Mechanics Symposium (NARMS). Houston: American Rock Mechanics Association, 2004: ARMA-04-498.

[10] 劉洋. 過程法構(gòu)建數(shù)字巖心技術(shù)[D].青島：中國石油大學 (華東),2007.

Liu Yang. Construction of Digital Core by Process-Based Simulation Method [D]. Qingdao: China University of Petroleum (East China), 2007.

[11]Youssef S, Rosenberg E, Gland N F, et al. High Resolution CT And Pore-Network Models To Assess Petrophysical Properties Of Homogeneous And Heterogeneous Carbonates[C]// Reservoir Characterization and Simulation Conference. Abu Dhabi: Society of Petroleum Engineers. doi：10.2118/111427-MS.

[12] Ryazanov A V, Sorbie K S, Dijke M I J V. Structure of Residual Oil as a Function of Wettability Using Pore-Network Modelling[J]. Advances in Water Resources, 2014, 63(2):11-21.

[13] 劉向君,朱洪林,梁利喜. 基于微CT技術(shù)的砂巖數(shù)字巖石物理實驗[J]. 地球物理學報, 2014, 57(4):1133-1140.

Liu Xiangjun, Zhu Honglin, Liang Lixi. Digital Rock Physics of Sandstone Based on Micro-CT Technology [J]. Chinese Journal of Geophysics, 2014, 57(4): 1133-1140.

[14] 李易霖, 張云峰, 叢琳,等. X-CT 掃描成像技術(shù)在致密砂巖微觀孔隙結(jié)構(gòu)表征中的應(yīng)用：以大安油田扶余油層為例[J]. 吉林大學學報(地球科學版), 2016, 46(2):379-387.

Li Yilin, Zhang Yunfeng, Cong Lin, et al. Application of X-CT Scanning Technique in the Characterization of Micro Pore Structure of Tight Sandstone Reservoir: An Example from Fuyu Oil Layer in Daan Oil Filed[J]. Journal of Jilin University (Earth Science Edition), 2016, 46(2): 379-387.

[15] Pan C, Hilpert M, Miller C T. Lattice-Boltzmann Simulation of Two-Phase Flow in Porous Media[J]. Water Resources Research, 2004, 40(1):62-74.

[16] 屈樂. 基于低滲透儲層的三維數(shù)字巖心建模及應(yīng)用[D].西安:西北大學,2014.

Qu Le. Modeling and Application of Three-Dimensional Digital Cores of Low Permeability Reservoir [D]. Xi’an: Northwestern University, 2014.

[17] 劉學鋒,張偉偉,孫建孟. 三維數(shù)字巖心建模方法綜述[J]. 地球物理學進展,2013,28(6):3066-3072.

Liu Xuefeng, Zhang Weiwei, Sun Jianmeng. Methods of Constructing 3-D Digital Cores: A Review [J]. Progress in Geophysical, 2013,28(6):3066-3072.

[18] 王晨晨,姚軍,楊永飛,等. 基于CT掃描法構(gòu)建數(shù)字巖心的分辨率選取研究[J]. 科學技術(shù)與工程,2013,13(4):1049-1052.

Wang Chenchen, Yao Jun, Yang Yongfei, et al. The Research of Selected Resolution Based on CT Scanning Method for Constructing Digital Core[J]. Science, Technology and Engineering, 2013,13(4):1049-1052.

[19] 張爽,周慧鑫,牛肖雪,等. 基于非局部均值濾波與時域高通濾波的非均勻性校正算法[J]. 光子學報,2014,43(1):147-150.

Zhang Shuang, Zhou Huixin, Niu Xiaoxue, et al. Temporal High-Pass Filter Nonunifority Correction Algorithm Based on Non-Local Means Filter for Infrared Focal Plane Array[J]. Acta Photonica Sinica, 2014, 43(1): 147-150.

[20] 楊山. 基于驅(qū)替實驗及數(shù)字巖心的微觀剩余油研究[D].青島:中國石油大學(華東),2015.

Yang Shan.The Research of Microscopic Remaining Oil on the Basis of Displacement Experiment and Digital Core[D]. Qingdao: China University of Petroleum (East China), 2015.

[21] 任奕明. 微觀水驅(qū)油網(wǎng)絡(luò)模擬研究[D].成都:西南石油大學,2014.

Ren Yiming. The Network Simulation of Microscopic Water Flooding[D]. Chengdu: Southwest Petroleum University, 2014.

[22] Mohammadmoradi P, Kantzas A. Pore-Scale Perme-ability Calculation Using CFD and DSMC Techniques[J]. Journal of Petroleum Science & Engineering, 2016, 146:515-525.