各向異性雙相介質(zhì)彈性波場(chǎng)褶積算法數(shù)值模擬①

宮 猛，李信富

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

0 引言

雙相彈性介質(zhì)理論認(rèn)為實(shí)際的地下介質(zhì)是由固相、液相組成的。固相的多孔隙骨架是均勻的、各向同性的彈性固體；液相的充滿孔隙空間的物質(zhì)是具有粘彈性的、不可壓縮的流體。特別是含油儲(chǔ)層具有較大的孔隙度，表現(xiàn)出明顯的雙相介質(zhì)性質(zhì)。雙相介質(zhì)理論與單相介質(zhì)理論不同，它充分地考慮了介質(zhì)的結(jié)構(gòu)、流體與氣體的特殊性質(zhì)、局部特性與整體效應(yīng)的關(guān)系，因此更能準(zhǔn)確地描述實(shí)際地層結(jié)構(gòu)和地層性質(zhì)，自然也就更能適應(yīng)越來(lái)越復(fù)雜的油氣儲(chǔ)藏勘探的實(shí)際需要，從而引起了國(guó)內(nèi)外地震學(xué)家和勘探地震學(xué)家們的高度重視，由此而發(fā)展起來(lái)的正演和反演研究具有更好的應(yīng)用前景。1951年，Gassmann提出了關(guān)于彈性波在多孔介質(zhì)中的傳播理論［1］，并建立了著名的Gassmann方程（反映了速度與孔隙度之間的定量關(guān)系）。之后，Biot根據(jù)潮濕土壤的電位特性和聲學(xué)中聲波的吸收特性，發(fā)展了Gassmann的流體飽和多孔隙雙相介質(zhì)理論［2－3］，奠定了雙相介質(zhì)波動(dòng)理論的基礎(chǔ)。Biot理論充分考慮了孔隙介質(zhì)的雙相特性，發(fā)現(xiàn)了第二類縱波，并指出粘滯力控制孔隙流體的相對(duì)運(yùn)動(dòng)是彈性波在孔隙介質(zhì)傳播過(guò)程中發(fā)生衰減的重要機(jī)理，并于后人的工作中得到了驗(yàn)證、發(fā)展和應(yīng)用。

由于Biot雙相介質(zhì)波動(dòng)方程在復(fù)雜地質(zhì)環(huán)境下沒(méi)有確定的解析解，只能通過(guò)數(shù)值方法求得，所以人們對(duì)Biot雙相介質(zhì)的數(shù)值模擬做了大量研究。Zhu和McMechan用有限差分法模擬了雙相介質(zhì)［4］；Hassanzadeh用有限差分法模擬了雙相介質(zhì)中縱波的傳播［5］；牛濱華等討論了裂隙含流體、氣體各向異性介質(zhì)一維二、三分量波動(dòng)方程［6］，然后用有限元方法進(jìn)行了波場(chǎng)數(shù)值模擬；N.Dai等用一階應(yīng)力－速度波動(dòng)方程模擬了各向異性雙相介質(zhì)［7］；牟永光基于Biot理論對(duì)雙相PTL、雙相EDA、及雙相PTL十EDA介質(zhì)中彈性波問(wèn)題進(jìn)行了深入研究［8］，給出了雙相各向異性介質(zhì)中彈性波方程的有限差分方法；邵秀民和藍(lán)志凌推導(dǎo)了有人工邊界時(shí)流體飽和多孔介質(zhì)波動(dòng)方程的有限元計(jì)算公式［9］；Arntsen和Carcione對(duì)雙相介質(zhì)中的慢縱波進(jìn)行了模擬研究［10］；楊寬德等從同時(shí)包含兩種力學(xué)機(jī)制的孔隙彈性波方程出發(fā)，利用有限差分法對(duì)含流體孔隙各向同性介質(zhì)中的地震波和聲波進(jìn)行了數(shù)值模擬［11］，并與基于Biot流動(dòng)的Biot理論之模擬結(jié)果進(jìn)行比較；楊頂輝利用有限元方法對(duì)雙相PTL介質(zhì)和雙相各向同性介質(zhì)中的彈性波傳播進(jìn)行了數(shù)值模擬［12］；劉洋等推導(dǎo)了各向異性雙相介質(zhì)的偽譜法數(shù)值解法和有限元解法［13－14］；孫衛(wèi)濤和楊慧珠采用交錯(cuò)網(wǎng)格技術(shù)建立了各向異性孔隙介質(zhì)波動(dòng)方程的高精度差分格式，并對(duì)這類差分格式的頻散特性和穩(wěn)定性作了詳細(xì)分析討論，解決了計(jì)算穩(wěn)定性和邊界反射問(wèn)題［15］；王東等利用基于Biot理論的孔隙彈性介質(zhì)的高階交錯(cuò)網(wǎng)格有限差分算法，模擬了具有隨機(jī)分布特征的多種流體飽和巖石中聲波［16］；裴正林基于Biot理論給出了三維雙相各向異性介質(zhì)應(yīng)力－速度彈性波方程交錯(cuò)網(wǎng)格任意偶數(shù)階精度有限差分解法，對(duì)三維雙相橫向各向異性介質(zhì)中彈性波場(chǎng)進(jìn)行了模擬［17－18］。

雙相各向異性介質(zhì)是非常復(fù)雜的介質(zhì)，也是目前地震學(xué)研究的前沿、熱點(diǎn)課題和難點(diǎn)課題之一。雙相各向異性理論是以雙相各向同性理論和地震各向異性理論為基礎(chǔ)，產(chǎn)生和發(fā)展的重要原因在于這種理論能夠更真實(shí)地描述地下介質(zhì)問(wèn)題。Biot建立了孔隙各向異性介質(zhì)理論，是研究雙相各向異性問(wèn)題的基礎(chǔ)。

本文通過(guò)對(duì)李信富等提出的方法進(jìn)行改進(jìn)［19］，推出各向異性雙相介質(zhì)中地震波傳播數(shù)值模擬計(jì)算的褶積微分算法，并對(duì)之進(jìn)行數(shù)值計(jì)算檢驗(yàn)。

1 雙相介質(zhì)波場(chǎng)模擬基本方程

Biot雙相介質(zhì)的運(yùn)動(dòng)平衡方程為

式中，ρ11，ρ12，ρ22為質(zhì)量系數(shù)，ρ11表示單元體中固體相對(duì)流體運(yùn)動(dòng)時(shí)固體部分總的等效質(zhì)量，ρ22為流體相對(duì)固體運(yùn)動(dòng)時(shí)流體部分總的等效質(zhì)量，ρ12表示流體和固體之間的質(zhì)量耦合系數(shù)，是粘滯、摩擦等效應(yīng)的綜合反映，又稱為視質(zhì)量。

為了得到一階速度－應(yīng)力方程，首先要將固相位移分量和流相位移分量從各自滿足的方程分離開(kāi)來(lái)。令式 （1）×ρ22－式（2）×ρ12，并整理得到固相位移分量公式：

令式 （1）×ρ12－式（2）×ρ11，并整理得到流相位移分量公式：

設(shè)固相速度矢量為v3×1＝ ［vx，vy，vz］，流相速度矢量為V3×1＝ ［Vx，Vy，Vz］，將速度對(duì)時(shí)間的一階導(dǎo)數(shù)替換式（3）與式（4）中位移對(duì)時(shí)間的二階導(dǎo)數(shù)，我們可以進(jìn)一步得到固相速度和流相速度的一階時(shí)間導(dǎo)數(shù)方程形式：

B3×3是耗散系數(shù)，對(duì)于各向同性Biot雙相介質(zhì)bx＝by＝bz；對(duì)于垂向橫向各向同性Biot雙相介質(zhì)bx＝by。

進(jìn)而，可以建立一階速度－應(yīng)力波動(dòng)方程各分量形式如下：

固相速度分量

流相速度分量

固相應(yīng)力分量

流相應(yīng)力

2 雙相介質(zhì)波場(chǎng)模擬數(shù)值檢驗(yàn)

為了檢驗(yàn)本文方法的有效性，利用該方法對(duì)均勻橫向各向同性雙相介質(zhì)和分層雙相介質(zhì)中彈性波場(chǎng)傳播進(jìn)行了數(shù)值模擬計(jì)算。

2.1 橫向各向同性雙相介質(zhì)褶積算法波場(chǎng)特征

本模型考察橫向各向同性雙相介質(zhì)中的地震波波場(chǎng)特征。網(wǎng)格點(diǎn)數(shù)為200×200；空間步長(zhǎng)為10 m；時(shí)間步長(zhǎng)為0.001s；介質(zhì)方位角為零。震源位于模型中心位置，采用傾斜集中力源激發(fā)，Richer子波主頻為15Hz，所用彈性參數(shù)如表1所示。圖1是不同分量的波場(chǎng)快照。

從波場(chǎng)快照中可以看出，（1）在方位角為零的橫向各向同性雙相介質(zhì)中，二維三分量數(shù)值模擬只能觀測(cè)到三種波，即快縱波qP1、快橫波qS1和慢縱波qP2，這三種波具有各向異性特性，且慢縱波屬于第二類波，其余為的一類波。（2）由于不同方向上地震波的傳播速度不同，地震波的波前面不再是圓形，而成為橢圓形，其曲率的大小由各向異性參數(shù)決定。（3）從波場(chǎng)快照中還可以看出，在二維三分量橫向各向同性雙相介質(zhì)波場(chǎng)模擬中x分量和z分量的波場(chǎng)特征相似，y分量波場(chǎng)有所區(qū)別，y分量是由固體骨架各向異性引起的。

表1 均勻雙相TIV介質(zhì)物性參數(shù)

圖1 橫向各向同性雙相介質(zhì)模型彈性波場(chǎng)快照（介質(zhì)傾角為0°；方位角為0°；t＝280ms）Fig.1 Snapshots of elastic wavefield in transverse homogeneous two－phase media（the dip is zero degree and the azimuths is zero degree.The propagation time of the wave is 280ms）.

2.2 兩層橫向各向同性雙相介質(zhì)波場(chǎng)模擬

本模型考察兩層橫向各向同性雙相介質(zhì)水平分界面上的的地震波波場(chǎng)特征。使用二維二分量格式進(jìn)行模擬，介質(zhì)方位角為零。模型如圖2所示，網(wǎng)格點(diǎn)數(shù)為200×200；空間步長(zhǎng)為10m；時(shí)間步長(zhǎng)為0.001s；震源位于（100，60）處，采用傾斜集中力源激發(fā)，Richer子波頻率為15Hz，主對(duì)稱軸為z軸。分界面在z＝80處，接收排列在z＝50處，與界面平行。兩層介質(zhì)的彈性參數(shù)見(jiàn)表2。圖3是各分量的波場(chǎng)快照。圖3中可見(jiàn)，快縱波、慢縱波和快橫波在界面上產(chǎn)生透射或反射后可以相互轉(zhuǎn)換，其速度大小依次為：快縱波、快橫波和慢縱波。由于各向異性介質(zhì)中地震波傳播的特殊性，只有在三維空間模擬才能全面認(rèn)識(shí)波場(chǎng)的空間變化。

圖2 兩層水平分界面模型Fig.2 Two－layered model with a horizontal subsurface.

從圖中可以看出：（1）雙相TIV介質(zhì)中的反射域彈性波有三類：直達(dá)類波：直達(dá)快縱波（qP1）、直達(dá)橫波（qS1）、直達(dá)慢縱波（qP2）；反射波類：快縱波的反射波（qP1qP1）、橫波的反射波（qS1qS1）、慢縱波的反射波（qP2qP2）；轉(zhuǎn)換波類：由快縱波形成的轉(zhuǎn)換橫波（qP1qS1）、由快縱波形成的轉(zhuǎn)換慢縱波（qP1qP2）、由慢縱波形成的轉(zhuǎn)換快縱波（qP2qP1）和由橫波形成的轉(zhuǎn)換慢縱波（qS1qP2）。（2）從第一類波轉(zhuǎn)換為第二類波后，該波型具有第二類波的性質(zhì)；同樣，從第二類波轉(zhuǎn)換為第一類波后，該波型具有第一類波的性質(zhì)。

表2 兩層雙相TIV介質(zhì)物性參數(shù)

圖3 兩層雙相TIV介質(zhì)中彈性波場(chǎng)快照（介質(zhì)傾角為0°；方位角為0°；t＝350ms）Fig.3 Snapshots of elastic wavefield in two－layered transverse homogeneous two－phase media（the dip is zero degree and the azimuths is zero degree.The propagation time of the wave is 350ms）.

3 結(jié)論與討論

本文利用時(shí)間錯(cuò)格有限差分－空間褶積微分算子法來(lái)計(jì)算非均勻介質(zhì)中的地震波場(chǎng)。時(shí)間錯(cuò)格有限差分－空間褶積微分算子法應(yīng)用于地震波場(chǎng)數(shù)值模擬的主導(dǎo)思想是：利用基于Forsyte廣義正交多項(xiàng)式褶積微分算子的有效表示計(jì)算波場(chǎng)對(duì)空間的偏導(dǎo)數(shù)，采用錯(cuò)格有限差分法計(jì)算對(duì)時(shí)間的偏導(dǎo)數(shù)，其計(jì)算思路類似于偽譜法。該方法是一種全新的、快速的、高精度的地震波場(chǎng)模擬方法，它相對(duì)于有限元法而言，可使用較大的網(wǎng)格；相對(duì)于有限差分法而言，頻散效應(yīng)小，計(jì)算精度高；相對(duì)于偽譜法而言，計(jì)算速度較快。

數(shù)值模擬結(jié)果表明，各向異性介質(zhì)中波場(chǎng)數(shù)值模擬的褶積算法非常準(zhǔn)確地模擬了各向異性介質(zhì)中的波場(chǎng)傳播過(guò)程，對(duì)各種波的刻畫非常清楚，為研究復(fù)雜介質(zhì)中地震波的傳播特征提供了一種新的方法選擇。

［1］Gassmann F.Elastic waves through a packing of spheres［J］.Geophysics，1951（16）：673－685.

［2］Biot M A.Theory of propagation of elastic waves in a fluidsaturated porous solid. Ⅰ.low－frequency range［J］.J.Acoust.Soc.Am.，1956（28）：168－178.

［3］Biot M A.Theory of propagation of elastic waves in a fluidsaturated porous solid.Ⅱ.higher－frequency range［J］.J.Acoust.Soc.Am.，1956（28）：179－191.

［4］Zhu X，G A McMechan.Numerical simulation of seismic responses of poroelastic reservoirs using Biot theory［J］.Geophysics，1991（56）：328－339.

［5］Hassanzadeh S.Acoustic modeling in fluid－saturated porous media［J］.Geophysics，1991（56）：424－435.

［6］牛濱華，吳有校，孫春巖.裂隙含流體氣體各向異性介質(zhì)波場(chǎng)數(shù)值模擬［J］.長(zhǎng)春地質(zhì)學(xué)院學(xué)報(bào)，1994，24（4）：454－460.

［7］N Dai，Vafidis A，Kanasewieh E R.Wave propagation in heterogeneous，porous media：A velocity－stress，finite difference method［J］.Geophysics，1995，60（2）：327－34.

［8］牟永光.儲(chǔ)層地球物理學(xué)［M］.北京：石油工業(yè)出版社，1996：5.

［9］邵秀民，藍(lán)志凌.流體飽和多孔介質(zhì)波動(dòng)方程的有限元解法［J］.地球物理學(xué)報(bào)，2000，43（2）：264－278.

［10］Arntsen B，Carcione J M.Numerical simulation of the Biot slow wave in water－saturated Nivelsteiner sandstone［J］.Geophysics，2001（66）：890－896.

［11］楊寬德，楊頂輝，王書強(qiáng).基于Biot－Squirt方程的波場(chǎng)模擬［J］.地球物理學(xué)報(bào)，2002，45（6）：853－861.

［12］楊頂輝.雙相各向異性介質(zhì)中彈性波方程的有限元解法及波場(chǎng)模擬［J］.地球物理學(xué)報(bào)，2002，45（4）：575－583.

［13］劉洋，李承楚.雙相各向異性介質(zhì)中彈性波傳播偽譜法數(shù)值模擬研究［J］.地震學(xué)報(bào)，2000，22（2）：132－138.

［14］劉洋，魏修成.雙相各向異性介質(zhì)中彈性波傳播有限元方程及數(shù)值模擬［J］.地震學(xué)報(bào)，2003，25（2）：154－162.

［15］孫衛(wèi)濤，楊慧珠.雙相各向異性介質(zhì)彈性波場(chǎng)有限差分正演模擬［J］.固體力學(xué)學(xué)報(bào)，2004，25（l）：21－28.

［16］王東，張海瀾，王秀明.部分飽和孔隙巖石中聲波傳播數(shù)值研究［J］.地球物理學(xué)報(bào)，2006，49（2）：524－532.

［17］裴正林.二維雙相各向異性介質(zhì)彈性波方程交錯(cuò)網(wǎng)格高階有限差分法模擬［J］.中國(guó)石油大學(xué)學(xué)報(bào)（自然科學(xué)版），2006，30（2）：16－20.

［18］裴正林.三維各向同性介質(zhì)彈性波方程交錯(cuò)網(wǎng)格高階有限差分法模擬［J］.石油物探，2006，41（2）：308－315.

［19］李信富，李小凡.地震波傳播的褶積微分算子法數(shù)值模擬［J］.地球科學(xué)－中國(guó)地質(zhì)大學(xué)學(xué)報(bào)，2008，33（6）：861－866.