各向異性雙相介質彈性波場褶積算法數值模擬①

宮 猛，李信富

（中國地質大學（北京）地球物理與信息技術學院，北京 100083）

0 引言

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

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

雙相各向異性介質是非常復雜的介質，也是目前地震學研究的前沿、熱點課題和難點課題之一。雙相各向異性理論是以雙相各向同性理論和地震各向異性理論為基礎，產生和發展的重要原因在于這種理論能夠更真實地描述地下介質問題。Biot建立了孔隙各向異性介質理論，是研究雙相各向異性問題的基礎。

本文通過對李信富等提出的方法進行改進［19］，推出各向異性雙相介質中地震波傳播數值模擬計算的褶積微分算法，并對之進行數值計算檢驗。

1 雙相介質波場模擬基本方程

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

式中，ρ11，ρ12，ρ22為質量系數，ρ11表示單元體中固體相對流體運動時固體部分總的等效質量，ρ22為流體相對固體運動時流體部分總的等效質量，ρ12表示流體和固體之間的質量耦合系數，是粘滯、摩擦等效應的綜合反映，又稱為視質量。

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

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

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

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

進而，可以建立一階速度－應力波動方程各分量形式如下：

固相速度分量

流相速度分量

固相應力分量

流相應力

2 雙相介質波場模擬數值檢驗

為了檢驗本文方法的有效性，利用該方法對均勻橫向各向同性雙相介質和分層雙相介質中彈性波場傳播進行了數值模擬計算。

2.1 橫向各向同性雙相介質褶積算法波場特征

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

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

表1 均勻雙相TIV介質物性參數

圖1 橫向各向同性雙相介質模型彈性波場快照（介質傾角為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 兩層橫向各向同性雙相介質波場模擬

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

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

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

表2 兩層雙相TIV介質物性參數

圖3 兩層雙相TIV介質中彈性波場快照（介質傾角為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 結論與討論

本文利用時間錯格有限差分－空間褶積微分算子法來計算非均勻介質中的地震波場。時間錯格有限差分－空間褶積微分算子法應用于地震波場數值模擬的主導思想是：利用基于Forsyte廣義正交多項式褶積微分算子的有效表示計算波場對空間的偏導數，采用錯格有限差分法計算對時間的偏導數，其計算思路類似于偽譜法。該方法是一種全新的、快速的、高精度的地震波場模擬方法，它相對于有限元法而言，可使用較大的網格；相對于有限差分法而言，頻散效應小，計算精度高；相對于偽譜法而言，計算速度較快。

數值模擬結果表明，各向異性介質中波場數值模擬的褶積算法非常準確地模擬了各向異性介質中的波場傳播過程，對各種波的刻畫非常清楚，為研究復雜介質中地震波的傳播特征提供了一種新的方法選擇。

［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］牛濱華，吳有校，孫春巖.裂隙含流體氣體各向異性介質波場數值模擬［J］.長春地質學院學報，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］牟永光.儲層地球物理學［M］.北京：石油工業出版社，1996：5.

［9］邵秀民，藍志凌.流體飽和多孔介質波動方程的有限元解法［J］.地球物理學報，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］楊寬德，楊頂輝，王書強.基于Biot－Squirt方程的波場模擬［J］.地球物理學報，2002，45（6）：853－861.

［12］楊頂輝.雙相各向異性介質中彈性波方程的有限元解法及波場模擬［J］.地球物理學報，2002，45（4）：575－583.

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

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

［15］孫衛濤，楊慧珠.雙相各向異性介質彈性波場有限差分正演模擬［J］.固體力學學報，2004，25（l）：21－28.

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

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

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

［19］李信富，李小凡.地震波傳播的褶積微分算子法數值模擬［J］.地球科學－中國地質大學學報，2008，33（6）：861－866.