雙相TI介質(zhì)中彈性波交錯(cuò)網(wǎng)格高階有限差分法數(shù)值模擬

尹學(xué)愛，邱光輝

（1.濱州學(xué)院，山東 濱州 256603；2.中化地質(zhì)礦山總局 山東地質(zhì)勘查院，山東 濟(jì)南 250013）

雙相TI介質(zhì)中彈性波交錯(cuò)網(wǎng)格高階有限差分法數(shù)值模擬

尹學(xué)愛1，邱光輝2

（1.濱州學(xué)院，山東 濱州 256603；2.中化地質(zhì)礦山總局 山東地質(zhì)勘查院，山東 濟(jì)南 250013）

應(yīng)用雙相介質(zhì)波動(dòng)方程，推導(dǎo)了雙相橫向各向同性介質(zhì)（TI）中波動(dòng)方程的有限差分格式，對(duì)雙相TI介質(zhì)中彈性波有限差分?jǐn)?shù)值進(jìn)行了模擬。結(jié)果表明，彈性波在雙相TI介質(zhì)中傳播時(shí)，除了存在常規(guī)的快縱波（qP1）和橫波以外，還存在慢縱波（qP2）。并且慢縱波的速度明顯小于快縱波，而且受耗散系數(shù)的影響衰減地很快，所以在實(shí)際中很難觀測(cè)到慢縱波?？炜v波在固相和流相中相位相同，而慢縱波在固相和流相中的相位相反。慢縱波在流相中振幅大，而在固相中的振幅較小。

雙相橫向各向同性介質(zhì)；交錯(cuò)網(wǎng)格；有限差分法；彈性波

0 前言

經(jīng)典的地震波理論只適合于研究固體或流體單相介質(zhì)中地震波傳播規(guī)律。然而，無論是砂巖儲(chǔ)層還是碳酸鹽儲(chǔ)層，或者是海底沉積物，都是由固體和流體兩種部份組成，即由固體和流體組成的雙相介質(zhì)或多相介質(zhì)。由于地震波在雙相或多相介質(zhì)傳播時(shí)，其傳播規(guī)律與在單相介質(zhì)中不同，因此，根據(jù)雙相介質(zhì)地震波理論，通過正、反演數(shù)值模擬，才能可靠準(zhǔn)確地確定儲(chǔ)層的厚度、空間展布，以及儲(chǔ)層的孔隙度、滲透率和油氣飽和度等儲(chǔ)層參數(shù)，所以對(duì)油氣勘探開發(fā)等實(shí)際工作具有十分重要的意義。

但是，該領(lǐng)域的研究還存在大量凾待解決的問題，其中最突出的是計(jì)算速度和穩(wěn)定性，國(guó)內(nèi)、外學(xué)者仍在進(jìn)行對(duì)于雙相介質(zhì)模型和高性能算法的研究。Gassmann［1］提出了關(guān)于彈性波在多孔介質(zhì)中的傳播理論，并建立了著名的Gassmann方程；Biot［2～4］根據(jù)潮濕土壤的電位特性和聲學(xué)中聲波的吸收特性，發(fā)展了Gassmann的流體飽和多孔隙雙相介質(zhì)理論，奠定了雙相介質(zhì)波動(dòng)理論的基礎(chǔ)，但由于Biot雙相介質(zhì)波動(dòng)方程在復(fù)雜地質(zhì)環(huán)境下沒有確定的解析解，只能通過數(shù)值方法求得，所以國(guó)內(nèi)、外對(duì)Biot雙相介質(zhì)的數(shù)值模擬做了大量研究；Zhu和McMechan［5］用有限差分法模擬了雙相介質(zhì)；N.Dai［6］等用一階應(yīng)力～速度波動(dòng)方程模擬了各向異性雙相介質(zhì)；孫衛(wèi)濤和楊慧珠［7］采用交錯(cuò)網(wǎng)格技術(shù)，建立了各向異性孔隙介質(zhì)波動(dòng)方程的高精度差分格式，并對(duì)這類差分格式的頻散特性和穩(wěn)定性作了詳細(xì)分析討論，解決了計(jì)算穩(wěn)定性和邊界反射問題；裴正林［8、9］基于Biot理論給出了三維雙相各向異性介質(zhì)應(yīng)力～速度彈性波方程交錯(cuò)網(wǎng)格任意偶階精度有限差分解法，對(duì)三維雙相各向異性介質(zhì)中彈性波場(chǎng)進(jìn)行了模擬；王秀明等［10、11］利用基于Biot理論的非均勻孔隙彈性介質(zhì)的高階交錯(cuò)網(wǎng)格有限差分算法，模擬了地震波在其中的傳播。

作者在Biot理論的基礎(chǔ)上，建立了Biot雙相介質(zhì)波動(dòng)方程交錯(cuò)網(wǎng)格差分格式，通過數(shù)值模擬對(duì)雙相橫向各向同性介質(zhì)中的彈性波傳播特征進(jìn)行了觀測(cè)分析。

1 原理

由飽和流體孔隙介質(zhì)的Biot理論可知，聲波傳播的Biot線性理論基于以下六種假設(shè)：

（1）波長(zhǎng)遠(yuǎn)大于顆粒尺寸，顆粒粒徑又大于孔隙尺寸。

（2）固體與流體的質(zhì)點(diǎn)運(yùn)動(dòng)相對(duì)較小，在線性彈性理論范圍內(nèi)。

（3）流體相在整個(gè)介質(zhì)中是連續(xù)的，而不連通的孔道可以視為固體骨架。

（4）固體骨架認(rèn)為是均勻和統(tǒng)計(jì)各向同性的，并且忽略了重力的影響。

（5）流體流動(dòng)是Poiseuille類型，巖石骨架與孔隙流體之間存在相對(duì)運(yùn)動(dòng)，且流體的流動(dòng)滿足廣義達(dá)西定律（這一限制只適用于低頻情況），并忽略熱效應(yīng)。

（6）介質(zhì)是完全飽和的。

基于上述假設(shè)，可以得到下面相關(guān)方程。

1.1 應(yīng)變～位移關(guān)系

若用u3×1＝ ［uxuyuz］T表示固相的平均位移矢量，U3×1＝ ［UxUyUz］T表示流相的平均位移矢量。則固相應(yīng)變分量可以用位移分量表示，固相正應(yīng)變?yōu)椋?/p>

固相切應(yīng)變?yōu)椋?/p>

固體部份的體應(yīng)變?chǔ)?，可由位移矢量?chǎng)的散度表示：

同樣地，流相部份的體應(yīng)變?chǔ)艦椋?/p>

1.2 應(yīng)力～應(yīng)變關(guān)系

根據(jù)廣義Hooke定律，雙相TI介質(zhì)中的應(yīng)力～應(yīng)變關(guān)系（即本構(gòu)方程）表示為

其中 σ6×1＝［σxxσyyσzzτyzτxzτxy］T是應(yīng)力矩陣的列矩陣表達(dá)式；e6×1＝［exxeyyezzeyz

exzexy］T是應(yīng)變矩陣的列矩陣表達(dá)式；S是流相應(yīng)力；ε是流相體應(yīng)變；D6×6是固相的物性矩陣，它包含與固相有關(guān)的彈性參數(shù)；R1×1是與流相有關(guān)的彈性參數(shù)；Q1×6是與固相和流相之間耦合有關(guān)的彈性參數(shù)，且Q6×1＝Q1×6T。

對(duì)于橫向各向同性雙相介質(zhì)，

其中 d66＝0.5（d11－d12），在D6×6中共有五個(gè)獨(dú)立的物性參數(shù)，Q1×6為：

其中 Q1＝Q2，共有二個(gè)獨(dú)立的物性參數(shù)，R1×1＝R中只有一個(gè)獨(dú)立物性參數(shù)。

1.3 應(yīng)力～位移關(guān)系

根據(jù)分析力學(xué)中的Lagrange方程，可導(dǎo)出雙相TI介質(zhì)中當(dāng)有耗散存在時(shí)的平衡運(yùn)動(dòng)方程式：

式中

對(duì)于各向同性雙相介質(zhì)b11＝b22＝b33，對(duì)于垂向橫向各向同性雙相介質(zhì)b11＝b22。在式（8）中ρ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ì)量。

2 有限差分算法的實(shí)現(xiàn)

目前有限差分法在數(shù)值模擬中具有其它方法不可替代的優(yōu)勢(shì)，差分離散的途徑主要有兩種：

（1）對(duì)單變量的二階波動(dòng)方程直接利用時(shí)空二階中心差分離散求解（規(guī)格網(wǎng)格差分）。

（2）將以位移表示的二階波動(dòng)方程化為以應(yīng)力和質(zhì)點(diǎn)速度表示的一階雙曲線方程組，用應(yīng)力和速度的交錯(cuò)網(wǎng)格求解，交錯(cuò)網(wǎng)格是一種較為先進(jìn)的差分格式。

作者在本文利用交錯(cuò)網(wǎng)格有限差分法原理，實(shí)現(xiàn)了雙相介質(zhì)相關(guān)微分方程向差分方程的轉(zhuǎn)化，為下一步數(shù)值模擬打下了基礎(chǔ)。

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

為了得到一階速度～應(yīng)力方程，首先要將固相位移分量和流相位移分量各自滿足的方程分離開來。

令式 （11）×ρ22－式（12）×ρ12，并整理得到固相位移分量公式為式（13）。

令式（11）×ρ12－式（12）×ρ11，并整理得到流相位移分量公式為式（14）。

設(shè)固相速度矢量為v3×1＝ ［vxvyvz］，流相速度矢量為V3×1＝ ［VxVyVz］，將速度對(duì)時(shí)間的一階導(dǎo)數(shù)替換成公式（13）和式（14）中位移對(duì)時(shí)間的二階導(dǎo)數(shù)，我們可以進(jìn)一步得到固相速度和流相速度的一階時(shí)間導(dǎo)數(shù)方程形式。

下面推導(dǎo)應(yīng)力分量的一階方程形式，雙相介質(zhì)的本構(gòu)方程為：

將方程式（16）兩邊同時(shí)對(duì)時(shí)間求導(dǎo)，我們就可以得到固相應(yīng)力和流相應(yīng)力的一階時(shí)間導(dǎo)數(shù)方程形式。對(duì)于固相有式（17）。

對(duì)于流相有式（18）。

其中

我們建立了雙相介質(zhì)一階速度～應(yīng)力波動(dòng)方程，接下來就可以建立雙相介質(zhì)交錯(cuò)網(wǎng)格差分離散格式［21］。首先看交錯(cuò)網(wǎng)格離散格式中雙相介質(zhì)相應(yīng)的波場(chǎng)分量和彈性參數(shù)的空間位置，見下頁表1和下頁圖1。

表1 彈性波場(chǎng)分量和彈性參數(shù)位置Tab.1 Elastic wave field components and elastic parameter

圖1 三維雙相介質(zhì)交錯(cuò)網(wǎng)格示意圖Fig.1 3Dtwo－phase medium staggered－grid

（1）固相位移各分量的交錯(cuò)網(wǎng)格差分離散格式為式（20）。

（2）流相位移各分量的交錯(cuò)網(wǎng)格差分離散格式為式（21）。

其中

（3）固相應(yīng)力各分量的交錯(cuò)網(wǎng)格差分格式為式（22）。

（4）流相應(yīng)力分量的交錯(cuò)網(wǎng)格差分格式為式（23）。

其中

上面離散格式中，DΔx、DΔy、DΔz是空間求導(dǎo)算子，分別代表了對(duì)x、y、z的導(dǎo)數(shù)；n表示時(shí)間上的離散點(diǎn)；i表示x方向上的離散點(diǎn)；j表示y方向上的離散點(diǎn)；k表達(dá)z方向上的離散點(diǎn)。

若令y方向的空間導(dǎo)數(shù)為零，則差分格式退化為x～z平面上的二維交錯(cuò)網(wǎng)格差分格式；若令y分量上的波場(chǎng)值都為零，則退化為只有x、z兩個(gè)分量的差分格式。作者在本文的模擬都是在二維二分量上求取的。

3 數(shù)值模擬結(jié)果及分析

二維均勻雙相橫向各向同性介質(zhì)模型彈性參數(shù)見表2。模型大小為400m×400m，震源為爆炸點(diǎn)源，高斯子波主頻為30Hz，震源位于模型中心，時(shí)間步長(zhǎng)為10－4s，網(wǎng)格大小為2.5m，交錯(cuò)網(wǎng)格差分精度為o（Δt2，Δx6）。

圖2和圖3（見下頁）是利用表2中模型2所得到的二維均勻雙相橫向各向同性（VTI）介質(zhì)彈性波波場(chǎng)模擬結(jié)果。

下頁圖4和后面的圖5是利用表2中模型1所得到的二維均勻雙相橫向各向同性（HTI）介質(zhì)彈性波波場(chǎng)模擬結(jié)果。

表2 均勻雙相介質(zhì)物性參數(shù)Tab.2 Parameters of medium physical

圖2 均勻雙相介質(zhì)中彈性波快照Fig.2 Snapshots of TI media

圖5 彈性波場(chǎng)快照與彈性波場(chǎng)記錄（t＝0.09s）Fig.5 Snapshots and time series for elastic wave field

彈性波在雙相各向異性介質(zhì)中傳播時(shí)存在四種波，即快縱波、快橫波、慢橫波和慢縱波。由于模擬為二維，且傾角和方位角為零度或者90°時(shí)，只能同時(shí)看到三種波，所以從圖2中可以看到快縱波qP1、快橫波qS1（SH波）和慢縱波qP2。從上頁圖4中可以看到快縱波qP1、快橫波qS2（SV波）和慢縱波qP2。這三種波具有各向異性特性，只有慢縱波屬于第二類波，其余為第一類波。

圖2是利用表2中模型2所得到的，二維均勻雙相橫向各向同性（VTI）介質(zhì)彈性波波場(chǎng)模擬結(jié)果。從圖2中可以看到，彈性波在雙相橫向各向同性（VTI）介質(zhì)中傳播時(shí)存在三種波，即快縱波qP1、快橫波qS1（SH）、慢縱波qP2。從圖2中可知，當(dāng)傾角和方位角都為零度時(shí)，快橫波和慢橫波在固流相中形成橫波分裂盲點(diǎn)，速度一樣，無法區(qū)分。在固相x、z分量上含有快縱波qP1、快橫波qS1和慢縱波qP2；在流相x、z分量上也含有快縱波qP1、快橫波qS1和慢縱波qP2，其中快縱波qP1和快橫波qS1在z分量的能量比x分量要強(qiáng)一些。從圖2中還可以看出，快縱波qP1和快橫波qS1在方位上各向異性，并且在一個(gè)方位上是耦合的，同時(shí)可以看出橫波波前產(chǎn)生尖角和三分叉現(xiàn)象。這一特殊現(xiàn)象使波場(chǎng)復(fù)雜化，并且容易引起誤解。

根據(jù)物性參數(shù)與速度的關(guān)系，得到速度，用速度驗(yàn)證波前面，得知本模擬結(jié)果是正確的。

圖3是z方向上雙相橫向各向同性（VTI）介質(zhì)中彈性波波場(chǎng)快照與彈性波場(chǎng)記錄的對(duì)比。其中檢波器道間距為5m和最小偏移距為零，檢波器排列在震源上方37.5m處水平接收。從圖3中可以看出，各種波對(duì)應(yīng)良好，只是在流相中的快橫波有頻散現(xiàn)象，這應(yīng)該是快縱波與快橫波耦合造成的。

圖4是利用表2中模型1所得到的二維均勻雙相橫向各向同性（HTI）介質(zhì)彈性波波場(chǎng)模擬結(jié)果。由圖4可知，在固相x、z分量上含有快縱波qP1、快橫波qS2（SV波）和慢縱波qP2，但是慢縱波qP2在x分量上能量很弱，在z分量上基本看不到。在流相x、z分量上也含有快縱波qP1、快橫波qS2和慢縱波qP2，快橫波qS2和慢縱波qP2是耦合的。從圖4整體上看，快縱波qP1能量在z方向上要比x方向上要小一些。

圖5是z方向上雙相橫向各向同性（HTI）介質(zhì)中彈性波波場(chǎng)快照與彈性波波場(chǎng)記錄的對(duì)比。其中檢波器道間距為5m和最小偏移距為零，在震源深度水平接收。從圖5中可以看出，波場(chǎng)快照與波場(chǎng)紀(jì)錄對(duì)應(yīng)良好。從圖5中可知，當(dāng)傾角和方位角都為零度時(shí)，快橫波和慢橫波在固流相中形成橫波分裂盲點(diǎn)，速度一樣，無法區(qū)分。慢縱波qP2在x、z方向上能量均較固相中要強(qiáng)很多，這進(jìn)一步證明了在流相中更容易觀測(cè)到慢縱波qP2。慢縱波qP2的振幅大小，主要受耗散系數(shù)b大小的影響。耗散系數(shù)越小，慢縱波越明顯，反之，則越不明顯。

根據(jù)物性參數(shù)與速度的關(guān)系，可得到速度，用速度驗(yàn)證波前面，得知本模擬結(jié)果是正確的。

4 結(jié)論

（1）根據(jù)Biot理論，作者對(duì)雙相介質(zhì)理論進(jìn)行了討論，推導(dǎo)了雙相橫向各向同性介質(zhì)波動(dòng)方程和雙相橫向各向同性速度～應(yīng)力彈性波方程，通過交錯(cuò)網(wǎng)格高階有限差分?jǐn)?shù)值模擬，研究了彈性波在雙相橫向各向同性介質(zhì)中的傳播情況。

（2）在雙相介質(zhì)VTI和HTI中，除了存在常規(guī)的縱波（快縱波）和橫波以外，還存在第二類縱波（慢縱波），慢縱波的速度明顯小于快縱波，而且受耗散系數(shù)的影響衰減地很快，所以在實(shí)際中很難觀測(cè)到第二類縱波。

（3）在雙相介質(zhì)中，快縱波在固相和流相中相位相同，而慢縱波在固相和流相中的相位相反。慢縱波在流相中振幅大，而在固相中的振幅較小，即在相同條件下，從流相中更容易觀測(cè)到慢縱波。

從數(shù)值模擬的結(jié)果可以看出：交錯(cuò)網(wǎng)格高階有限差分法是一種模擬精度高、計(jì)算效率好，實(shí)用性強(qiáng)的彈性波場(chǎng)數(shù)值模擬方法。

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

［2］ BIOT M A.Theory of propagation of elastic waves in fluid filled porous solid I.higher requency range［J］.J.Acou.Soc.Am.，1956，28（2）：179.

［3］ BIOT M A.Mechanics of deformation and acoustic propagation in porous media［J］.J.appl.Phys.，1962，33（4）：1482.

［4］ BIOT M A.Generalized theory of acoustic propagation in porous dissipative media［J］.J.Acoust.Soc.Am.，1962，34（5）：1254.

［5］ ZHU X，MCMECHAN G A.McMechan.Numerical simulation of seismic responses of poroelastic reservoirs using Biot theory［J］.Geophysics，1991，56（3）：328.

［6］ DAI N，VAFIDIS A，KANASEWICH E R.Wave propagation in heterogeneous，porous media：avelocity－stress，finite－difference method［J］.Geophysics，1995，60（2）：327.

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

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

［9］ 裴正林.雙相各向異性介質(zhì)彈性波傳播交錯(cuò)網(wǎng)格高階有限差分法模擬［J］.石油地球物理勘探，2006，41（2）：137.

［10］王秀明，張海瀾，王東.利用高價(jià)交錯(cuò)網(wǎng)格有限差分法模擬地震波在非均勻孔隙介質(zhì)中的傳播［J］.地球物理學(xué)報(bào)，2003，46（6）：842.

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

P 631.4

A

10.3969／j.issn.1001－1749.2012.05.06

1001—1749（2012）05—0533—08

濱州學(xué)院科研基金（BZXYL1107）

2011－10－21 改回日期：2012－05－30

尹學(xué)愛（1979－），女，漢，山東鄒平人，教師，目前主要從事彈性波數(shù)值模擬方面的研究。