印興耀， 秦秋萍， 宗兆云

1 中國石油大學(華東)， 青島　266580 2 海洋國家實驗室海洋礦產(chǎn)資源評價與探測技術(shù)功能實驗室， 青島　266071

?

數(shù)字巖石物理中彈性參數(shù)的有限差分計算方法

印興耀1,2， 秦秋萍1,2， 宗兆云1,2

1 中國石油大學(華東)， 青島266580 2 海洋國家實驗室海洋礦產(chǎn)資源評價與探測技術(shù)功能實驗室， 青島266071

數(shù)字巖石物理利用三維成像技術(shù)和數(shù)學方法，建立數(shù)字巖心，開展多場物理響應(yīng)模擬，計算巖石的等效彈性參數(shù)，為巖石物理學研究開拓了新的領(lǐng)域.本文發(fā)展了一種基于高階有限差分的巖石模量數(shù)值計算的新方法，該方法便于理解，容易實現(xiàn)，占用內(nèi)存相對較少，計算效率高，結(jié)果合理，彌補了常規(guī)巖石物理實驗周期長，成本高，誤差大等不足.該方法將三維數(shù)字巖心樣本嵌入一個具有與巖心樣本骨架顆粒相同彈性性質(zhì)的區(qū)域中擴充成一個新模型，測量由數(shù)字巖石樣本非均勻結(jié)構(gòu)帶來的縱波或橫波峰值振幅的時間差，利用該時間差(相對參考模型)估算縱橫波的等效速度，進而求取等效彈性模量.理論模型和實際巖心的計算結(jié)果表明，數(shù)值模擬結(jié)果與實驗結(jié)果有較高的吻合度，驗證了該方法的合理性.

數(shù)字巖石物理； 高階有限差分； 彈性模量； 速度

1　引言

巖石的宏觀物理性質(zhì)是微觀結(jié)構(gòu)骨架、孔隙等各組成成分的綜合效應(yīng)，對于體積模量、剪切模量、速度等，這些彈性參數(shù)在地球物理勘探領(lǐng)域等實際工作中發(fā)揮著重要作用(劉向君等，2014).巖石物理是連接儲層彈性參數(shù)與物性參數(shù)的橋梁，也是地球物理反演的基礎(chǔ)(宗兆云等，2012a；印興耀等，2015).隨著計算機技術(shù)與成像技術(shù)的提高，數(shù)字巖石物理技術(shù)得到了迅速發(fā)展，根據(jù)巖石微觀結(jié)構(gòu)信息重建反映巖石真實孔隙空間的三維數(shù)字模型，以實際地質(zhì)資料為約束，在計算機上運用數(shù)值計算方法研究孔隙結(jié)構(gòu)、裂縫特征及其與巖石物理特性(如電阻率、彈性模量、速度與波場、滲透率和核磁共振)之間的關(guān)系(孫建孟等，2012；朱偉和單蕊，2014).基于數(shù)字巖心的物理性質(zhì)研究是巖石物理的重要發(fā)展方向，特別是隨著頁巖氣、致密油氣等開發(fā)，發(fā)展數(shù)字巖心巖石物理方法對其具有重要的指導(dǎo)意義.相比于傳統(tǒng)的巖石物理實驗，數(shù)字巖石物理在巖石彈性、物性多尺度模擬方面有較大的優(yōu)勢，針對不同類型儲層(如疏松砂巖、碳酸鹽巖、致密砂巖、頁巖氣、頁巖油等)巖樣，利用數(shù)字化成像技術(shù)，通過巖心整體掃描建立多尺度巖石模型，明確不同類型儲層巖石骨架及孔隙的多尺度特征并將其參數(shù)化，抽象出不同類型儲層巖石數(shù)學模型，發(fā)展基于模型驅(qū)動的疊前地震反演方法，為地震反演提供理論指導(dǎo)(宗兆云等，2012b；宗兆云，2013).與傳統(tǒng)的巖石物理實驗相比，數(shù)字巖石物理實驗操作簡單，速度快而且成本低，對實驗的條件要求低，同時避免了人為操作、實驗環(huán)境及儀器的影響所帶來的誤差，節(jié)省大量人力物力.基于數(shù)字巖心的巖石物理數(shù)值模擬對巖心不具有破壞性，巖心可重復(fù)使用.為了建立多物理場屬性的內(nèi)在關(guān)系，可以在同一塊數(shù)字巖心上進行不同物理場的模擬，如電阻率、彈性模量、速度與波場、滲透率等，并且可以測量傳統(tǒng)巖石物理實驗難以直接測量的物理量，如三相相對滲透率(孫建孟等，2012).此外傳統(tǒng)的巖石物理實驗對于低孔低滲等復(fù)雜油氣儲層，巖心驅(qū)替測試困難，因此對于裂縫發(fā)育的碳酸鹽巖儲層巖石物理參數(shù)的準確獲得也具有重要意義(陳懷震等，2014a).基于數(shù)字巖心可以在微觀上模擬巖石彈性模量隨裂縫密度、裂縫發(fā)育方向的變化，同時可以模擬裂縫介質(zhì)的波場特征，有助于裂縫的識別，建立各向異性巖石物理模型為AVAZ反演提供理論基礎(chǔ)(陳懷震等，2014b).對于碳酸鹽巖及頁巖、油砂等難以取到代表性的巖心，取心費用高、成功率低，對于疏松巖心保存困難，巖石物理實驗難以開展，無法定量研究微觀參數(shù)對巖石宏觀物理屬性的影響，利用數(shù)字巖石物理實驗可代替?zhèn)鹘y(tǒng)巖石物理實驗彌補這些不足.除了上述優(yōu)勢，通過調(diào)整數(shù)字巖心的微觀參數(shù)，利用數(shù)字巖石物理實驗有利于認識儲層參數(shù)對巖石物理屬性的影響，由于數(shù)字巖心的微觀結(jié)構(gòu)直觀可見，為直接計算幾何特征參數(shù)提供了便利，數(shù)值模擬的邊界條件靈活，可連續(xù)呈現(xiàn)物理場的變化過程(朱偉和單蕊，2014).中國的地層結(jié)構(gòu)復(fù)雜，在三維數(shù)字巖心方面的研究起步較晚，為了縮小國際差距，數(shù)字巖石物理是發(fā)展趨勢，地球物理工作者應(yīng)該實施數(shù)字化智能化的戰(zhàn)略.

在油氣勘探與開發(fā)領(lǐng)域數(shù)字巖石物理主要在測井及油氣開發(fā)滲流模擬領(lǐng)域得到研究及應(yīng)用，屈樂(2014)基于低滲透儲層建立三維數(shù)字巖心模型及電性模擬研究.數(shù)字巖石物理在地震勘探領(lǐng)域以彈性模擬為主，目前處于起步階段，研究及應(yīng)用范圍有一定的局限性.Saenger等(2000)在標準交錯網(wǎng)格基礎(chǔ)上提出了旋轉(zhuǎn)交錯網(wǎng)格有限差分方法.Saenger等(2004)假設(shè)數(shù)字巖心模型孔隙中飽含高密度流體，計算等效彈性參數(shù)，這一研究為數(shù)字巖石物理實驗在地震勘探方面的研究奠定了基礎(chǔ).Saenger等(2007)利用旋轉(zhuǎn)交錯網(wǎng)格有限差分方法研究巖石孔隙空間氣、水分布對波傳播的影響，模擬結(jié)果與實驗測量結(jié)果接近，表明該方法的高精度性與魯棒性.基于三維數(shù)字巖石物理，假設(shè)介質(zhì)含有裂縫，利用旋轉(zhuǎn)交錯網(wǎng)格有限差分進行波場動力學模擬(Saenger and Shapiro, 2002).通過廣義Maxwell型和旋轉(zhuǎn)交錯網(wǎng)格有限差分，對黏滯性孔隙流體進行數(shù)值模擬，得到縱橫波速度(Saenger et al., 2011).Saxena和Mavko(2015)基于數(shù)字巖心研究雙相介質(zhì)中流體或固體替代下剪切模量和體積模量的計算方法，討論了孔隙填充物性質(zhì)對彈性模量的影響.對于孔隙聯(lián)通性良好的多礦物巖石，Brown和Korringa(1975)擴展Gassmann流體替換方程用于彈性模量的計算，但是實際應(yīng)用中很少采用此方程， Saxena等(2015)基于三維數(shù)字巖心提出了流體替換的新方法，該方法在計算巖石基質(zhì)模量時，不依賴Gassmann流體替換中的Voigt-Reuss-Hill平均模量或Hashin-Shtrikman界限.除此之外，利用數(shù)字巖心開展準靜態(tài)模擬，建立波速、衰減系數(shù)與孔隙流體黏滯度的關(guān)系，建立計算數(shù)字巖心(黏)彈性參數(shù)的新方法.

巖石彈性微觀數(shù)值模擬方法主要有聲格子、有限元以及旋轉(zhuǎn)交錯網(wǎng)格有限差分方法等(Andr? et al，2013a).一些學者應(yīng)用有限元方法進行了材料的靜力學模擬，計算了復(fù)合材料的等效彈性參數(shù)及巖石的彈性參數(shù)(Garboczi and Day, 1995; Bohn and Garboczi, 2003)，但有限元方法運行速度慢，計算效率低.Del Valle-García和Sánchez-Sesma(2003)利用聲格子方法模擬了聲波在孔隙介質(zhì)中的傳播，但由于該方法的計算量較大，目前還只限于二維數(shù)值模擬研究.本文直接在三維數(shù)字巖心基礎(chǔ)上利用有限差分方法來模擬波的傳播過程，測量由數(shù)字巖石樣本非均勻結(jié)構(gòu)帶來的峰值振幅的時間差，完成彈性模量的計算，驗證了該方法的有效性與合理性，高階有限差分方法占用內(nèi)存相對較少，計算速度快，模擬結(jié)果合理，為數(shù)字巖石物理實驗提供了一套新的研究思路.

2　方法原理

三維標量波方程為

(1)

其中，u(x,y,z,t)為波場，v(x,y,z)為介質(zhì)速度.

對方程(1)進行泰勒級數(shù)展開,得到方程(2)和(3)為

(3)

相加式(2)和式(3)得：

(4)

化簡式(4)可寫為

(5)

取等號右端的前三項，得到關(guān)于時間的二階差分格式為

+O(Δt4)，

(6)

式(6)在時間方向上的截斷誤差為O(Δt4)，是二階精度的差分格式.

關(guān)于x、y、z方向空間二階偏導(dǎo)數(shù)的高階差分格式推導(dǎo)過程是類似的，在此僅對X方向展開討論.假設(shè)X方向上的截斷誤差為O(Δx2N)，N是大于2的偶數(shù).

將u(x+nΔx)和u(x-nΔx)進行泰勒級數(shù)展開，有

(7)

(8)

相加式(7)和式(8)，得：

+O(Δx2N+2).

(9)

將上述n個方程同時乘以ω1，ω2，…，ωn然后相加并化簡整理，寫成矩陣形式為

(10)

其系數(shù)矩陣即為范德蒙矩陣.

(11)

(12)

(13)

根據(jù)前面的推導(dǎo)，截斷誤差為O(Δx2N,Δy2N,Δz2N,Δt2)，方程(1)的高階差分格式可表示為

(14)

圖1　嵌入數(shù)字巖石樣本的新模型(a)和參考模型(b)Fig.1　New model of embedded digital rock samples (a) and reference model (b)

該方法的基本思想是利用有限差分來模擬波的傳播過程，在孔隙尺度遠遠小于波長的情況下，研究非均質(zhì)三維數(shù)字巖心樣本的速度.為了計算巖石的有效彈性參數(shù)，將數(shù)字巖石樣本、參考樣本分別嵌入均勻彈性區(qū)域形成新模型和參考模型如圖1所示，假設(shè)該彈性區(qū)域和參考樣本具有與數(shù)字巖石樣本骨架顆粒相同的彈性性質(zhì).選擇最佳匹配層(Perfectly Matched Layer，簡稱PML)吸收邊界條件來消除邊界效應(yīng).波源在模型頂部激發(fā)，縱波或橫波在模型中傳播，檢波器在底部接收，測量由數(shù)字巖石樣本非均勻結(jié)構(gòu)帶來的峰值振幅的時間差，利用時間差(相對于參考模型)估算等效速度，進而求取等效彈性模量.具體的計算公式為：

圖2　計算流程與方法框圖Fig.2　Chart of computation process and the method

(15)

(16)

(17)

(18)

(19)

(20)

(21)

(22)

(23)

其中h表示被嵌入均勻彈性區(qū)域的高度，H表示三維數(shù)字巖心樣本的高度，V0表示均勻彈性區(qū)域的速度，Veff表示三維數(shù)字巖心樣本的等效速度，t1表示波源激發(fā)平面波經(jīng)過新模型的時間，t2表示波源激發(fā)平面波經(jīng)過參考模型的時間，α、β分別表示估算得到的縱、橫波速度，Keff、Geff分別表示三維數(shù)字巖心樣本的等效體積模量、剪切模量，ρf、ρs、ρeff分別表示三維數(shù)字巖心樣本的孔隙流體密度、骨架密度、等效密度，φ為三維數(shù)字巖心樣本的孔隙度.利用得到的縱、橫波速度根據(jù)公式(19)和(20)即可得到三維數(shù)字巖心樣本的等效彈性模量，如公式(22)和(23)所示.

用于本文研究的理論模型是細層層狀介質(zhì)模型，三維數(shù)字巖心模型為Fontainebleau(楓丹白露)和Berea(貝雷)砂巖模型.以四塊Fontainebleau砂巖及一塊Berea砂巖的三維數(shù)字巖心樣本為例，在三維數(shù)字巖心的基礎(chǔ)上利用高階有限差分計算巖石的體積模量和剪切模量.三維數(shù)字巖心的分辨率分別為 5.68 μm/像素，5.345 μm/像素，巖心尺寸分別為200×200×200像素，400×400×400像素.Fontainebleau砂巖為研究提供了一種理想的三維數(shù)字巖心.首先，F(xiàn)ontainebleau砂巖的彈性模量實驗數(shù)據(jù)充分，且?guī)r石孔隙度的變化范圍廣，為對比實驗結(jié)果和數(shù)值模擬結(jié)果提供了方便.其次，F(xiàn)ontainebleau砂巖均質(zhì)性好，巖石骨架礦物成分單一，只有石英顆粒，而且該類砂巖孔隙結(jié)構(gòu)簡單，只存在粒間孔隙.在進行數(shù)字巖心數(shù)值模擬的研究中巖心的掃描，圖像的獲得是前提.劉學鋒等(2013)對三維數(shù)字巖心建模做了大量研究.如何利用X射線CT建立三維數(shù)字巖心及掃描后的圖像處理劉學峰(2010)做了詳細介紹.除此之外，Andr?等(2013b)對不同巖石類型的數(shù)字巖心的成像做了大量研究.

3　模型試算

3.1細層層狀介質(zhì)-Backus平均

該模型上部介質(zhì)(綠色)縱波速度Vp(1)=5200 m·s-1，橫波速度Vs(1)=2700 m·s-1，密度ρ1=2450 kg·m-3，厚度d1=0.75 m；下部介質(zhì)(黃色)縱波速度Vp(2)=2900 m·s-1，橫波速度Vs(2)=1400 m·s-1，密度ρ2=2340 kg·m-3，厚度d2=0.5 m.根據(jù)Backus平均計算公式可以得出理論垂向等效縱波速度為3761.0 m·s-1，理論垂向等效橫波速度為1854.8 m·s-1；數(shù)值模擬結(jié)果分別為3675.3 m·s-1、1753.2 m·s-1，如圖4所示.誤差分別為2.3%、5.3%.由此可見，數(shù)值模擬結(jié)果與理論計算結(jié)果的吻合程度較高.說明在三維數(shù)字巖心基礎(chǔ)上，采用高階有限差分方法計算三維數(shù)字巖心樣本的速度是可行的，為模量計算打下了基礎(chǔ)，驗證了該方法的有效性.

3.2Fontainebleau(楓丹白露)和Berea(貝雷)砂巖

數(shù)值模擬中，設(shè)巖石骨架的體積模量K=37 GPa，剪切模量G=44 GPa，密度ρs=2.65 g·mm-3.孔隙中飽含流體水，體積模量K=1.02 GPa，剪切模量G=0 GPa，密度ρf=1.00 g·mm-3.砂巖孔隙度φ分別為5%、10%、15%和20%，Berea砂巖孔隙度φ=19.6%.Fontainebleau砂巖及Berea砂巖三維數(shù)字巖心的體積模量和剪切模量的模擬結(jié)果如圖7和圖8所示，圖中綠色正方形數(shù)據(jù)點為實驗結(jié)果(Han et al., 1986)，黃色圓形數(shù)據(jù)點表示Fontainebleau砂巖數(shù)值模擬結(jié)果，紅色圓形數(shù)據(jù)點表示Berea砂巖數(shù)值模擬結(jié)果.巖石彈性模量的數(shù)值模擬結(jié)果與實驗結(jié)果的吻合程度較高.說明在三維數(shù)字巖心基礎(chǔ)上，采用高階有限差分方法計算巖石的彈性模量是可行的，驗證了該方法的有效性.

圖3　細層層狀介質(zhì)模型Fig.3　Schematic model of thin-layered media

圖4　數(shù)值模擬結(jié)果與理論計算結(jié)果比較Fig.4　Comparison of numerical simulation and theoretical results

圖5　三維Fontainebleau模型Fig.5　The 3D model of Fontainebleau

圖6　Berea模型Fig.6　The 3D model of Berea

圖7　體積模量數(shù)值模擬結(jié)果Fig.7　Simulation results of bulk moduli

圖8　剪切模量數(shù)值模擬結(jié)果Fig.8　Simulation results of shear moduli

4　有限差分與有限元方法對比

圖9　有限差分法與有限元法的體積模量計算比較Fig.9　Comparison between the finite element method and finite difference method in simulation results of bulk moduli

圖10　有限差分法與有限元法的剪切模量計算比較Fig.10　Comparison between the finite element method and finite difference method in simulation results of shear moduli

三維Fontainebleau砂巖巖心的分辨率為 5.68 μm/像素，巖心尺寸為200×200×200像素.以孔隙度分別為5%、10%、15%和20%的三維Fontainebleau 砂巖的數(shù)值模擬結(jié)果為例，從圖9和圖10有限元及有限差分的計算結(jié)果不難發(fā)現(xiàn).孔隙度在10%左右的時候有限元的計算結(jié)果要好于有限差分，孔隙度為5%、15%和20%時有限差分的計算結(jié)果與實驗結(jié)果更接近.利用有限差分法和有限元法分別對孔隙度為20%的三維Fontainebleau數(shù)字巖心進行數(shù)值模擬計算，計算時間分別是11757 s、12694 s，如圖11.從而說明了有限差分方法計算精度高，速度快，結(jié)果合理.

圖11　有限差分法與有限元法的計算效率比較Fig.11　Comparison between the finite element method and finite difference method in computational efficiency

5　結(jié)論

基于三維數(shù)字巖心，利用高階有限差分方法對現(xiàn)有的巖石物理等效模型進行驗證，計算巖石的體積模量和剪切模量，數(shù)值模擬結(jié)果與實驗結(jié)果相符，驗證了該方法的可行性與合理性.建立了模量與孔隙流體的關(guān)系，形成計算數(shù)字巖心彈性參數(shù)的新方法.與有限元計算方法相比有限差分計算速度快，精度較高，結(jié)果合理.基于三維數(shù)字巖心利用有限差分方法來模擬波傳播過程計算巖石的等效模量應(yīng)用前景廣泛，為數(shù)字巖石物理的研究提供了一條新的途徑，為地球物理勘探人員提供了一套可借鑒的研究思路.

Andr? H, Combaret N, Dvorkin J, et al. 2013a. Digital rock physics benchmarks—Part II: Computing effective properties.Computers&Geosciences, 50: 33-43.

Andr? H, Combaret N, Dvorkin J, et al. 2013b. Digital rock physics benchmarks—Part I: Imaging and segmentation.Computers&Geosciences, 50: 25-32.

Bohn R B, Garboczi E J. 2003. User manual for finite element and finite difference programs: A parallel version of NISTIR 6269. National Institute of Standards and Technology.

Brown R J S, Korringa J. 1975. On the dependence of the elastic properties of a porous rock on the compressibility of the pore fluid.Geophysics, 40(4): 608-616, doi: 10.1190/1.1440551.

Chen H Z, Yin X Y, Gao C G, et al. 2014b. AVAZ inversion for fluid factor based on fracture anisotropic rock physics theory.ChineseJ.Geophys. (in Chinese), 57(3): 968-978, doi: 10.6038/cjg20140326.

Chen H Z, Yin X Y, Zhang J Q, et al. 2014a. Seismic inversion for

fracture rock physics parameters using azimuthally anisotropic elastic impedance.ChineseJ.Geophys. (in Chinese), 57(10): 3431-3441, doi: 10.6038/cjg20141029.

Del Valle-García R, Sánchez-Sesma F J. 2003. Rayleigh waves modeling using an elastic lattice model.GeophysicalResearchLetters, 30(16): 1866, doi: 10.1029/2003GL017600.

Garboczi E J, Day A R. 1995. An algorithm for computing the effective linear elastic properties of heterogeneous materials: Three-dimensional results for composites with equal phase Poisson ratios.JournaloftheMechanicsandPhysicsofSolids, 43(9): 1349-1362.

Han D H, Nur A, Morgan D. 1986. Effects of porosity and clay content on wave velocities in sandstones.Geophysics, 51(11): 2093-2107.

Liu X F. 2010. Numerical simulation of elastic and electrical properties of rock based on digital cores [Ph. D. thesis] (in Chinese). Qingdao: China University of Petroleum.

Liu X F, Zhang W W, Sun J M. 2013. Methods of constructing 3-D digital cores: A review.ProgressinGeophys. (in Chinese), 28(6): 3066-3072, doi: 10.6038/pg20130630.

Liu X J, Zhu H L, Liang L X. 2014. Digital rock physics of sandstone based on micro-CT technology.ChineseJ.Geophys. (in Chinese), 57(4): 1133-1140, doi: 10.6038/cjg20140411.

Qu L. 2014. Modeling and application of three-dimensional digital core based on low permeability reservoir [Ph. D. thesis] (in Chinese). Xi′an: Northwest University.

Saenger E H, Ciz R, Krüger O S, et al. 2007. Finite-difference modeling of wave propagation on microscale: A snapshot of the work in progress.Geophysics, 72(5): SM293-SM300.

Saenger E H, Enzmann F, Keehm Y, et al. 2011. Digital rock physics: Effect of fluid viscosity on effective elastic properties.JournalofAppliedGeophysics, 74(4): 236-241.

Saenger E H, Gold N, Shapiro S A. 2000. Modeling the propagation of elastic waves using a modified finite-difference grid.WaveMotion, 31(1): 77-92.

Saenger E H, Krüger O S, Shapiro S A. 2004. Numerical considerations of fluid effects on wave propagation: Influence of the tortuosity.GeophysicalResearchLetters, 31(21): L21613.

Saenger E H, Shapiro S A. 2002. Effective velocities in fractured media: a numerical study using the rotated staggered finite-difference grid.GeophysicalProspecting, 50(2): 183-194.

Saxena N, Mavko G. 2015. The embedded-bound method for estimating the change in rock moduli under pore fill and mineral phase substitution.Geophysics, 80(3): L1-L10.

Saxena N, Mavko G, Mukerji T. 2015. Fluid substitution in multimineralic rocks with large mineral stiffness contrast.Geophysics, 80(3): L11-L33.

Sun J M, Jiang L M, Liu X F, et al. 2012. Log application and prospect of digital core technology.WellLoggingTechnology

(in Chinese), 36(1): 1-7.

Yin X Y, Zong Z Y, Wu G C. 2015. Research on seismic fluid identification driven by rock physics.ScienceChinaEarthSciences, 58(2): 159-171, doi: 10.1007/s11430-014-4992-3.

Zhu W, Shan R. 2014. Progress of digital rock physics.OilGeophysicalProspecting(in Chinese), 49(6): 1138-1146.

Zong Z Y, Yin X Y, Wu G C. 2012a. Fluid identification method based on compressional and shear modulus direct inversion.ChineseJ.Geophys. (in Chinese), 55(1): 284-292, doi: 10.6038/j.issn.0001-5733.2012.01.028.

Zong Z Y, Yin X Y, Zhang F, et al. 2012b. Reflection coefficient equation and pre-stack seismic inversion with Young’s modulus and Poisson ratio.ChineseJ.Geophys. (in Chinese), 55(11): 3786-3794, doi: 10.6038/j.issn.0001-5733.2012.11.025.

Zong Z Y. 2013. Methodologies of model driven inversion with Pre-stack seismic data [Ph. D. thesis] (in Chinese). Qingdao: China University of Petroleum.

附中文參考文獻

陳懷震, 印興耀, 張金強等. 2014a. 基于方位各向異性彈性阻抗的裂縫巖石物理參數(shù)反演方法研究. 地球物理學報, 57(10): 3431-3441, doi: 10.6038/cjg20141029.

陳懷震, 印興耀, 高成國等. 2014b. 基于各向異性巖石物理的縫隙流體因子AVAZ反演. 地球物理學報, 57(3): 968-978, doi: 10.6038/cjg20140326.

劉學鋒. 2010. 基于數(shù)字巖心的巖石聲電特性微觀數(shù)值模擬研究[博士論文]. 青島: 中國石油大學.

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

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

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

孫建孟, 姜黎明, 劉學鋒等. 2012. 數(shù)字巖心技術(shù)測井應(yīng)用與展望. 測井技術(shù), 36(1): 1-7.

印興耀, 宗兆云, 吳國忱. 2015. 巖石物理驅(qū)動下地震流體識別研究. 中國科學: 地球科學, 45(1): 8-21.

朱偉, 單蕊. 2014. 虛擬巖石物理研究進展. 石油地球物理勘探, 49(6): 1138-1146.

宗兆云, 印興耀, 吳國忱. 2012a. 基于疊前地震縱橫波模量直接反演的流體檢測方法. 地球物理學報, 55(1): 284-292, doi: 10.6038/j.issn.0001-5733.2012.01.028.

宗兆云, 印興耀, 張峰等. 2012b. 楊氏模量和泊松比反射系數(shù)近似方程及疊前地震反演. 地球物理學報, 55(11): 3786-3794, doi: 10.6038/j.issn.0001-5733.2012.11.025.

宗兆云. 2013. 基于模型驅(qū)動的疊前地震反演方法研究[博士論文]. 青島: 中國石油大學(華東).

(本文編輯張正峰)

Simulation of elastic parameters based on the finite difference method in digital rock physics

YIN Xing-Yao1,2, QIN Qiu-Ping1,2, ZONG Zhao-Yun1,2

1ChinaUniversityofPetroleum,Qingdao266580,China2LaboratoryforMarineMineralResources,QingdaoNationalLaboratoryforMarineScienceandTechnology,Qingdao266071,China

Digital core is a new field of petrophysics experiments, which simulates responses of multiple physics fields and calculates equivalently elastic parameters of rocks through three-dimensional imaging techniques and mathematical methods. In this paper, a new approach is proposed to acquire modulus of rocks based on the high-order finite difference method, which is easy to understand and implement, taking up less memory and making up the weakness of the conventional rock physics experiment with long periods, high costs and big errors. For the purpose of estimating elastic moduli available, the digital rock sample is embedded in a homogeneous elastic region assigned elastic properties of the grain material. With receivers at the bottom of the model, it is possible to measure the time-delay of the peak amplitude of the plane P-or S-waves caused by the inhomogeneous structure of the digital rock sample. With the time-delay (compared to the reference model) we estimate the effective velocity of both the compressional and the shear waves. The study on a theoretical model and digital rock samples indicates that simulation results agree well with the experimental solution, which verifies the validity of the method.

Digital rock physics; High-order finite difference; Elastic modulus; Velocity

10.6038/cjg20161031.

國家自然科學基金-石油化工基金聯(lián)合重點項目(U1562215)、國家油氣重大專項課題(2016ZX05024-004)、中央高?；究蒲袠I(yè)務(wù)費專項資金和中石化地球物理重點實驗室開放基金(33550006-15-FW2099-0027)聯(lián)合資助.

印興耀，男，中國石油大學(華東)教授、博士生導(dǎo)師，從事地球物理理論與方法方面的工作. E-mail：xyyin@upc.edu.cn

10.6038/cjg20161031

P631

2015-09-18，2016-05-13收修定稿

印興耀， 秦秋萍， 宗兆云. 2016. 數(shù)字巖石物理中彈性參數(shù)的有限差分計算方法. 地球物理學報，59(10):3883-3890,

Yin X Y, Qin Q P, Zong Z Y. 2016. Simulation of elastic parameters based on the finite difference method in digital rock physics.ChineseJ.Geophys. (in Chinese),59(10):3883-3890,doi:10.6038/cjg20161031.