改進BISQ模型的雙相介質地震波場數值模擬及頻散校正

周竹生，唐磊

(中南大學 地球科學與信息物理學院，湖南 長沙，410083)

相對單相介質而言，雙相介質能夠更好地描述實際地層結構和性質，因此，雙相介質的數值模擬更有利于實際介質中地震波傳播規律的認識和研究。在雙相介質的早期研究工作中，Biot假設孔隙介質理論中的固體與流體的相對運動屬于Poiseuille型流動[1-2]，然而，隨著雙相介質研究的不斷深入，研究發現孔隙介質中的流體還存在一種對地震波衰減和頻散起關鍵作用的流動方式即噴射流動。其中，Biot流動力學機制描述的是宏觀現象，而噴射流動機制描述的是局部特征。Dvorkin等[3]將這 2種機制結合起來，提出了BISQ模型，使得模型結構更貼近實際。Mamadou[4]提出改進BISQ模型，其優勢在于在描述流體壓力的表達式中不再含有不便于實際描述的微觀參量(如特征噴射流動長度)，更加方便于地震波波場模擬的實現。本研究采用交錯網格有限差分法對改進BISQ模型的雙相介質地震波場進行數值模擬。數值模擬方法的共性是需要在時空域進行網格離散，以差分代替微分，由此所得到的地震波場必然存在數值頻散現象。為了壓制數值頻散，本研究利用通量校正傳輸(FCT)技術修正了交錯網格有限差分算法，確保在滿足一定計算精度的前提下，可以采用較低的差分階數，以提高計算效率。

1 基本方程

1.1 改進BISQ模型的雙相介質彈性波動方程

在時空域內，基于改進BISQ模型的二維各向同性雙相介質彈性波動方程可以表示為[4]:

1.2 一階速度-應力方程

在二維情況下，基于改進BISQ模型的雙相介質運動平衡方程為：

進一步推導，將位移分量對時間的一階導數換成質點振動速度分量，得到一階速度-應力方程為：

2 一階彈性波方程高階差分近似

2.1 時間域的2M階差分近似

在使用交錯網格法求解一階速度-應力方程時，速度和應力分別是在t+Δt/2和t時刻進行計算的，對于速度與應力分量Vi和σij，利用Taylor公式得到2M階精度的時間差分近似格式：

式中：Δt為時間步長；O(Δt2m)為高階無窮小量，當M=1時，上式即為傳統的二階差分精度近似。

2.2 空間域的2N階差分近似

在計算空間導數時采用下式所示的 2N階差分近似表達式：

2.3 交錯網格高階有限差分格式的構建

對一階速度-應力方程進行高階差分近似，構建了空間2N階和時間2階精度的有限差分格式，其中：i，j和k分別表示x，z和t的離散序號；U，W，P，Q 和 S 分別表示固相 vx，vz，σxx，σzz和 σxz；u，w，F和R分別表示流相的Vx，Vz，S和P。

3 一階方程FCT有限差分算法

Boris和 Book發展的通量校正傳輸方法(FCT)成功應用于聲波方程[15]，FCT的基本思路是將修正的漫射應用于頻散傳輸格式。在進行校正過程中分為漫射和反漫射2個階段，計算的波場被改變以壓制數值頻散引起的偽波動。對于一階速度-應力波動方程(8)，通過交錯網格差分計算得到 k+1時刻的速度和應力值，為了提高計算效率，只需對差分計算得到的速度分量，，和的值進行校正。

(1) 計算k-1時刻的彌散通量：

式中：η1為漫射參數，它是常量或一個線性函數，其取值取決于有限差分階段頻散誤差。

(2) 計算k+1時刻的彌散通量：

式中，η2的取值與η1不同。這是因為振幅和分辨率的損失主要有2方面的原因：一是傳統有限差分運算引起的頻散，二是人為加入的漫射。反漫射運算不僅要補償人為加入漫射引起的損失，而且要補償傳統有限差分運算帶來的振幅損失。

4 波場數值模擬及頻散校正

4.1 波場模擬

為檢驗模擬算法的正確性，建立一個各向同性雙相介質均勻半空間模型，模型網格點數為200×200，模擬過程中為消除人工邊界的影響。采用Cerjan等[8]提出的吸收邊界條件，為保證計算過程中的穩定性，選取空間網格間距為5 m，時間間隔為5 ms，震源位于模型中心，采用主頻為30 Hz的雷克子波，模型的固體參數為：拉梅系數 12.6 GPa，剪切模量為6.72 GPa，固體體積模量為56.7 GPa，密度為2 822 kg/m3，耦合附加密度為85 kg/m3，縱波速度為3 km/s，橫波速度2 000 m/s；流體參數：體積模量為1.51 GPa，密度為 554 kg/m3，孔隙率為 0.248 8，黏滯系數為1×10-6Pa·s，滲透率為 15 md。

采用精度為 O(Δt2+Δx4)和 O(Δt2+Δx16)的交錯網格有限差分法，對固體和流體的水平分量和垂直分量分別進行計算，結果分別如圖1和圖2所示。

從圖1和圖2可以看出：在雙相介質中存在第二類縱波，從內到外依次為慢縱波(p)，橫波(S)和快縱波(P)，其衰減規律明顯不同于第一類縱波，第一類縱波在固相和流相都引起較強的振動，第二類縱波引起較弱的固相振動，較強的流相振動，這種關系取決于雙相介質的性質。其中第一類縱波的相位在固相和流相部分是相同的，第二類縱波的相位在固相和流相部分則相反。通過圖1與圖2對比可以發現：在時間差分精度一定時，空間差分精度階數較低時，頻散嚴重，隨著階數的增加，頻散降低，波場模擬的精度逐漸提高。說明采用交錯網格高階差分能夠對數值頻散起到一定的壓制。

4.2 加FCT修正的交錯網格差分法進行波場模擬及頻散校正

圖3和圖4所示分別為利用加FCT修正的不同精度的交錯網格差分法進行波場模擬得到的波場快照圖。與未進行FCT修正的圖1和圖2對比，可以看到頻散現象得到明顯地壓制，說明FCT方法能夠有效壓制在計算過程中產生的數值頻散，采用低階精度的差分方程加FCT法修正后也可以得到較高階精度的正演模擬結果。交錯網格高階差分與FCT方法相結合，在保持波場特征情況下能有效提高波場模擬的精度和計算效率。

圖1 100 ms時波場快照(時間2階，空間4階)Fig.1 Snapshot of 100 ms (two-order time and four-order space)

圖2 100 ms時波場快照(時間2階，空間16階)Fig.2 Snapshot of 100 ms (two-order time and sixteen-order space)

圖3 FCT修正后的100 ms波場快照(時間2階，空間4階)Fig.3 Snapshot of 100 ms (two-order time and four-order space) by using FCT

圖4 FCT修正后的100 ms波場快照(時間2階，空間16階)Fig.4 Snapshot of 100 ms (two-order time and sixteen-order space) by using FCT

5 結論

(1) 將交錯網格差分法運用到各向同性雙相介質中，給出不同時間和空間精度的一階速度-應力波動方程交錯網格差分算法格式，并提出帶FCT修正的交錯網格差分算法，運用交錯網格差分算法實現各向同性雙相介質地震波場正演數值計算。

(2) 雙相各向同性介質中存在快縱波、慢縱波和橫波，得到3種波在不同相介質中表現的波場特征。

(3) 對比不同差分精度計算得到的波場快照圖，當時間差分精度固定時，空間差分精度隨著階數的增加，頻散降低，差分精度的提高對頻散起到壓制作用。

(4) FCT方法能有效壓制在計算過程中產生的數值頻散，用低階精度的差分方程也可以得到高階精度的正演模擬效果，并且計算效率也得到了很大提高。

[1] Biot M A. Theory of propagation of elastic waves in a fluid-saturate porous solid, low-frequency range[J]. J Acoust Soc Amer, 1956, 28: 168-178.

[2] Biot M A. Theory of propagation of elastic waves in a fluid-saturate porous solid, higher-frequency range[J]. J Acoust Soc Amer, 1956, 28: 179-191.

[3] Dvorkin J, Nur A. Dynamic poroelasticity: A unified model with the squirt and the Boit mechanisms[J]. Geophysics, 1993, 58(4):524.

[4] Mamadou S. Acoustic wave propagation in saturated porous media: Reformulation of the Biot/Squirt flow theory[J]. Journal of Applied Geophysics, 2000, 44(4): 313-325.

[5] 孟慶生, 何樵登, 朱建偉, 等. 基于 BISQ模型雙相各向同性介質中地震波場數值模擬[J]. 吉林大學學報: 地球科學版,2003, 33(2): 217-196.MENG Qing-sheng, HE Qiao-deng, ZHU Jian-wei, et al.Seismic modeling in isotropic porous media based on BISQ model[J]. Journal of Jilin University: Earth Science Edition,2003, 33(2): 217-196.

[6] 劉洋, 李承楚, 牟永光. 任意偶數階精度有限差分法數值模擬[J]. 石油地球物理勘探, 1998, 33(1): 275-284.LIU Yang, LI Cheng-chu, MOU Yong-guang. Finite-difference numerical modeling of any even-order accuracy[J]. OGP, 1998,33(1): 1-10.

[7] 董良國, 馬在田, 曹景忠, 等. 一階彈性波動方程交錯網格高階差分解法穩定性研究[J]. 地球物理學報, 2000, 43(6):853-861.DONG Liang-guo, MA Zai-tian, CAO Jing-zhong, et al. A study on stability of the staggered-grid high-order difference method of first-order elastic wave equation[J]. Chinese Journal of Geophysics, 2000, 43(6): 853-861.

[8] Cerjan C, Kosloff D, Kosloff R, et al. A nonreflecting boundary condition for discrete acoustic and elastic wave equation[J].Geohysics, 1985, 50(4): 705.

[9] 牟永光. 儲層地球物理學[M]. 北京: 石油工業出版社, 1996 27-38.MOU Yong-guang. Reservoir geophysics[M]. Beijing: Oil Industry Press, 1996: 27-38.

[10] 楊寬德, 楊頂輝, 王書強. 基于Biot-Squirt方程的波場模擬[J].地球物理學報, 2002, 45(6): 853-861.YANG Kuan-de, YANG Ding-hui, WANG Shu-qiang. Wavefield simulation based on the Biot-Squirt equation[J]. Chinese Journal of Geophysics, 2002, 45(6): 853-861.

[11] 吳國忱, 王華忠. 波場模擬中的數值頻散分析與校正策略[J].地球物理學進展, 2005, 20(1): 58-65.WU Guo-chen, WAN Hua-zhong. Analysis of numerical dispersion in wave-field simulation[J]. Progress in Geophysics,2005, 20(1): 58-65.

[12] Dablain M A. The application of high differencing to the scalar wave equation[J]. Geophysics, 1985, 51(1): 54-66.

[13] XIAO Shao-ping. An FE-FCT method with implicit functions for the study of shock wave propagation in solids[J]. Wave Motion,2004, 40: 263-276.

[14] TONG Fei. Finite-difference modeling and depth-migration via flux-corrected transport[J]. SEG Expanded Abstracts 1993, 12(5):1074-1077.

[15] Book D L, Boris J P, Hain K. Flux-corrected transport II:Generalization of the method[J]. J Comput Phys, 1975, 18:248-283.