雙相TI介質中彈性波交錯網格高階有限差分法數值模擬

尹學愛，邱光輝

（1.濱州學院，山東 濱州 256603；2.中化地質礦山總局 山東地質勘查院，山東 濟南 250013）

雙相TI介質中彈性波交錯網格高階有限差分法數值模擬

尹學愛1，邱光輝2

（1.濱州學院，山東 濱州 256603；2.中化地質礦山總局 山東地質勘查院，山東 濟南 250013）

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

雙相橫向各向同性介質；交錯網格；有限差分法；彈性波

0 前言

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

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

作者在Biot理論的基礎上，建立了Biot雙相介質波動方程交錯網格差分格式，通過數值模擬對雙相橫向各向同性介質中的彈性波傳播特征進行了觀測分析。

1 原理

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

（1）波長遠大于顆粒尺寸，顆粒粒徑又大于孔隙尺寸。

（2）固體與流體的質點運動相對較小，在線性彈性理論范圍內。

（3）流體相在整個介質中是連續的，而不連通的孔道可以視為固體骨架。

（4）固體骨架認為是均勻和統計各向同性的，并且忽略了重力的影響。

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

（6）介質是完全飽和的。

基于上述假設，可以得到下面相關方程。

1.1 應變～位移關系

若用u3×1＝ ［uxuyuz］T表示固相的平均位移矢量，U3×1＝ ［UxUyUz］T表示流相的平均位移矢量。則固相應變分量可以用位移分量表示，固相正應變為：

固相切應變為：

固體部份的體應變θ，可由位移矢量場的散度表示：

同樣地，流相部份的體應變ε為：

1.2 應力～應變關系

根據廣義Hooke定律，雙相TI介質中的應力～應變關系（即本構方程）表示為

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

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

對于橫向各向同性雙相介質，

其中 d66＝0.5（d11－d12），在D6×6中共有五個獨立的物性參數，Q1×6為：

其中 Q1＝Q2，共有二個獨立的物性參數，R1×1＝R中只有一個獨立物性參數。

1.3 應力～位移關系

根據分析力學中的Lagrange方程，可導出雙相TI介質中當有耗散存在時的平衡運動方程式：

式中

對于各向同性雙相介質b11＝b22＝b33，對于垂向橫向各向同性雙相介質b11＝b22。在式（8）中ρ11、ρ12、ρ22為質量系數，ρ11表示單元體中固相相對于流相運動時固相部份總的等效質量；ρ22為流相相對于固體相運動時流體部份總的等效質量；ρ12表示流相和固相之間的質量耦合系數，是粘滯、摩擦等效應的綜合反映，又稱為視質量。

2 有限差分算法的實現

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

（1）對單變量的二階波動方程直接利用時空二階中心差分離散求解（規格網格差分）。

（2）將以位移表示的二階波動方程化為以應力和質點速度表示的一階雙曲線方程組，用應力和速度的交錯網格求解，交錯網格是一種較為先進的差分格式。

作者在本文利用交錯網格有限差分法原理，實現了雙相介質相關微分方程向差分方程的轉化，為下一步數值模擬打下了基礎。

雙相介質的運動平衡方程式為：

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

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

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

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

下面推導應力分量的一階方程形式，雙相介質的本構方程為：

將方程式（16）兩邊同時對時間求導，我們就可以得到固相應力和流相應力的一階時間導數方程形式。對于固相有式（17）。

對于流相有式（18）。

其中

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

表1 彈性波場分量和彈性參數位置Tab.1 Elastic wave field components and elastic parameter

圖1 三維雙相介質交錯網格示意圖Fig.1 3Dtwo－phase medium staggered－grid

（1）固相位移各分量的交錯網格差分離散格式為式（20）。

（2）流相位移各分量的交錯網格差分離散格式為式（21）。

其中

（3）固相應力各分量的交錯網格差分格式為式（22）。

（4）流相應力分量的交錯網格差分格式為式（23）。

其中

上面離散格式中，DΔx、DΔy、DΔz是空間求導算子，分別代表了對x、y、z的導數；n表示時間上的離散點；i表示x方向上的離散點；j表示y方向上的離散點；k表達z方向上的離散點。

若令y方向的空間導數為零，則差分格式退化為x～z平面上的二維交錯網格差分格式；若令y分量上的波場值都為零，則退化為只有x、z兩個分量的差分格式。作者在本文的模擬都是在二維二分量上求取的。

3 數值模擬結果及分析

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

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

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

表2 均勻雙相介質物性參數Tab.2 Parameters of medium physical

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

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

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

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

根據物性參數與速度的關系，得到速度，用速度驗證波前面，得知本模擬結果是正確的。

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

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

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

根據物性參數與速度的關系，可得到速度，用速度驗證波前面，得知本模擬結果是正確的。

4 結論

（1）根據Biot理論，作者對雙相介質理論進行了討論，推導了雙相橫向各向同性介質波動方程和雙相橫向各向同性速度～應力彈性波方程，通過交錯網格高階有限差分數值模擬，研究了彈性波在雙相橫向各向同性介質中的傳播情況。

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

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

從數值模擬的結果可以看出：交錯網格高階有限差分法是一種模擬精度高、計算效率好，實用性強的彈性波場數值模擬方法。

［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］ 孫衛濤，楊慧珠.雙相各向異性介質彈性波場有限差分正演模擬［J］.固體力學學報，2004，25（1）：21.

［8］ 裴正林.三維雙相各向異性介質彈性波方程交錯網格高階有限差分法模擬［J］.中國石油大學學報，2006，30（2）：16.

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

［10］王秀明，張海瀾，王東.利用高價交錯網格有限差分法模擬地震波在非均勻孔隙介質中的傳播［J］.地球物理學報，2003，46（6）：842.

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

P 631.4

A

10.3969／j.issn.1001－1749.2012.05.06

1001—1749（2012）05—0533—08

濱州學院科研基金（BZXYL1107）

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

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