迎風(fēng)型組合緊致差分的二維地震波場數(shù)值模擬

周誠堯,汪 勇,蔡偉祥,桂志先，陳驁卓

(1. 長江大學(xué) 油氣資源與勘探技術(shù)教育部重點實驗室，湖北 武漢 430100;2. 中石化重慶涪陵頁巖氣勘探開發(fā)有限公司，重慶 408014)

0 前 言

地震正演模擬研究是當(dāng)前研究地震波傳播性質(zhì)的重要的方法。隨著地震正演模擬的發(fā)展，人們對于數(shù)值結(jié)果的精度要求越來越高[1]。常用的數(shù)值模擬方法主要有射線追蹤法和波動方程法，屬于有限差分法。波動方程法然而傳統(tǒng)的有限差分?jǐn)?shù)值模擬方法有著兩個較為明顯的缺點：一是需要用來進(jìn)行數(shù)值模擬的網(wǎng)格點數(shù)較多進(jìn)而增大了計算時間；二是傳統(tǒng)的顯式差分格式具有較大的數(shù)值頻散[2]。緊致差分方法是為應(yīng)對這種情形應(yīng)運而生的。Lele于1992年提出了對稱型緊致差分格式。傅德薰和馬延文在緊致差分格式中引入迎風(fēng)機(jī)制[3]，提出了五格點五階精度的迎風(fēng)緊致格式[4]。引入迎風(fēng)機(jī)制理論上可以提高緊致差分格式的穩(wěn)定性。在相同的計算網(wǎng)格中，超緊致差分格式有著更高的精度和分辨率[5-6]，而超緊致差分格式的穩(wěn)定性相比于普通緊致差分的穩(wěn)定性有一定幅度的下降。本文所使用的迎風(fēng)超緊致有限差分?jǐn)?shù)值模擬法是在Chu的三點六階超緊致差分的基礎(chǔ)上引入迎風(fēng)機(jī)制而推導(dǎo)出的三點五階迎風(fēng)超緊致差分格式，其穩(wěn)定性近一步提高[7]。在地震數(shù)值模擬應(yīng)用方面，王書強將緊致差分格式用于彈性波方程的數(shù)值模擬[8]，本文所提及方法鮮有文獻(xiàn)報道，本文屬首次將迎風(fēng)超緊致差分格式用于二維地震波場數(shù)值模擬，而相比于傳統(tǒng)的有限差分以及普通緊致差分?jǐn)?shù)值模擬，本文所采用的迎風(fēng)超緊致有限差分?jǐn)?shù)值模擬方法擁有高精度、高效率以及低頻散的優(yōu)點，相較于原本的三點六階超緊致差分格式，三點五階迎風(fēng)超緊致差分格式具有更高的穩(wěn)定性。鑒于上述優(yōu)點，本文將利用迎風(fēng)型超緊致有限差分格式，對聲波方程時間四階離散格式的位移場空間偏導(dǎo)數(shù)進(jìn)行離散，實現(xiàn)地震波場的數(shù)值模擬。

1 迎風(fēng)超緊致差分格式的原理及應(yīng)用

迎風(fēng)超緊致有限差分格式是將迎風(fēng)機(jī)制引入到超緊致有限差分格式中來。Krishnan提出了如下的緊致差分格式[9]：

(1)

其中a1、a0、a2、b0、b1、b2、c1、c2、c0、c3、c4為差分公式系數(shù)，h表示網(wǎng)格間距。該公式可以達(dá)到六階精度，相比于一般的緊致差分格式，所用到的網(wǎng)格點數(shù)更少。本文中所選取需要引入迎風(fēng)機(jī)制的超緊致格式為三點六階超緊致差分格式，如下所示：

(2)

上述公式引入迎風(fēng)機(jī)制后得到本文所研究的三點五階迎風(fēng)超緊致差分格式如下所示：

(3)

接下來則是將上述迎風(fēng)超緊致差分方法運用于求解聲波方程初值問題中來[10]。在各向均勻同性介質(zhì)中聲波方程為如下公式：

(4)

其中的u表示聲波在x與z方向上的位移。v則是聲波在各向均勻同性介質(zhì)中的速度。對于上式利用二階中心差分對時間導(dǎo)數(shù)離散可得到如下公式：

(5)

其中un+1為n+1時刻位移un-1為n-1時刻位移un為n時刻位移。

迎風(fēng)型組合緊致差分格式的應(yīng)用關(guān)鍵點在于對其空間導(dǎo)數(shù)的計算。本文利用公式(3)求公式(5)中u的z方向上的空間偏導(dǎo)數(shù)，矩陣表示為AF=Bu，其中A為公式(3)左端的差分系數(shù)矩陣，B為公式(3)右端的差分系數(shù)矩陣，u為位移場矩陣。F為位移場空間一階和二階導(dǎo)數(shù)矩陣。這些矩陣分別如下所示，系數(shù)矩陣A：

(6)

F與u分別為：

(7)

(8)

系數(shù)矩陣B如下所示：

(9)

因此對于聲波方程的空間導(dǎo)數(shù)利用迎風(fēng)超緊致差分格式逼近后得到如下關(guān)系式：AF=Bu，得F=A-1·Bu。若將F定義為聲波方程空間z方向上的一二階導(dǎo)數(shù)的集合，那么設(shè)FT為x方向上的一二階導(dǎo)數(shù)集合，可知FT=u·(B)T·(A-1)T。而聲波方程所需要的為二階導(dǎo)數(shù)部分。因此取z、x方向上的二階導(dǎo)數(shù)部分并命名為F1與FT1，如下所示：

同理對上述結(jié)果再次進(jìn)行矩陣變化可得到z與x方向上的三四階導(dǎo)數(shù)的集合F′及F′T。對x、z依次求導(dǎo)的集合為G，取z、x方向上的四階導(dǎo)數(shù)部分及對x對z依次二階導(dǎo)數(shù)部分，并命名為F′1、F′T1與G1?？芍?/p>

(12)

公式(5)變?yōu)槿缦滦问剑?/p>

un+1=2un-un-1+(Δtv)2(FT1+F1)+(Δtv)2(F′1+2G1+F′T1)

(13)

利用關(guān)系式(12)及該聲波方程所設(shè)定初始條件即可進(jìn)行聲波方程的數(shù)值模擬。

2 用于聲波方程中的頻散分析

頻散現(xiàn)象是指波的傳播相速度與群速度不一致而引起的現(xiàn)象[11-12]。這里我們利用Fourier分析。由于聲波方程只需要離散二階空間導(dǎo)數(shù)部分，因此我們只需得知二階導(dǎo)數(shù)的修正波數(shù)即可。頻散分析大體步驟如下。

1)根據(jù)修正波數(shù)的定義利用x及z方向上的修正波數(shù)kxkz與空間步長h乘積的二次方公式求出差分格式二階導(dǎo)數(shù)x方向上的修正波數(shù)和z方向上的修正波數(shù)，其中kxh=khcosθ，kzh=khsinθ,θ即波的傳播方向與x軸的夾角。

2)再將上述格式代入聲波方程中得到如下格式：

2cos(kv′Δt)=2-(vΔt/h)2[(kx″h)2+(kz″h)2]+(vΔt/h)4[(kx″h)4+2(kx″h·kz″h)2+(kz″h)4]/12

(14)

kx″h與kz″h分別為差分格式二階導(dǎo)數(shù)x方向上的修正波數(shù)和z方向上的修正波數(shù)與空間步長的乘積，v′為數(shù)值波速，v為波的真速度。

4)做出頻散曲線圖并與幾種常見的緊致差分格式的頻散曲線圖做對比并分析。

根據(jù)頻散壓制原則，如果不存在數(shù)值頻散則速度比恒等于1。若偏離1越大，則說明該方法的數(shù)值頻散越嚴(yán)重。反之，則說明該方法能更好地壓制數(shù)值頻散。圖1顯示了迎風(fēng)超緊致格式與其它兩種相同精度得差分格式在不同θ時候的速度比曲線，其中α為courant數(shù)。

1 USCD; 2 UCD; 3 CD

如上圖所示，橫坐標(biāo)φ為波數(shù)與空間步長的乘積，縱坐標(biāo)γ為速度比，其中數(shù)值1表示修正波速與理論波速一致。曲線1為迎風(fēng)超緊致差分方法速度比曲線，曲線2為五階迎風(fēng)緊致差分方法速度比曲線，曲線3為五階一般緊致差分速度比曲線。從圖中的速度比曲線可以看出：在不同的角度下，3種方法的的偏移趨勢沒有明顯的變化；從3種曲線的對比看出迎風(fēng)超緊致方法在偏離1時所對應(yīng)的橫坐標(biāo)數(shù)值，因此根據(jù)圖中數(shù)值可以得到上述4種方法在保證壓制頻散的時候所需要的最低網(wǎng)格點數(shù)為4.0、5.6、6.7個。幾種對比方法中，迎風(fēng)超緊致三點五階方法對于頻散的壓制效果最好。

3 用于聲波方程中的穩(wěn)定性分析

對于穩(wěn)定性分析，本文采用ui,j=ξneIh(kxi+kzj)代入位移場及位移場空間一階和二階導(dǎo)數(shù)中得到如下表達(dá)式：

(15)

α=vmaxΔt/h<0.81

這里和未引入迎風(fēng)機(jī)制的三點六階超緊致差分格式的穩(wěn)定性進(jìn)行一個對比，如表1所示。

表1 不同方法courant數(shù)對比

迎風(fēng)超緊致格式具有較高的courant數(shù)，即穩(wěn)定性高，普通緊致六階差分方法的courant數(shù)(即穩(wěn)定性)明顯高于三點六階超緊致差分格式，而與未引入迎風(fēng)機(jī)制的三點六階超緊致格式相比，三點五階迎風(fēng)超緊致格式的穩(wěn)定性有著明顯的提高。

4 截斷誤差分析

對于精度的分析本文采用截斷誤差的對比分析。表2為由Taylor截斷誤差分析所得，內(nèi)容為幾種方法的二階導(dǎo)數(shù)的截斷誤差?？梢缘玫揭韵陆Y(jié)論：迎風(fēng)超緊致格式不僅具有超出中心差分格式的精度，和標(biāo)準(zhǔn)緊致差分格式相比，在同樣利用三點的情況下具有更高的精度。

表2 二階導(dǎo)數(shù)不同方法截斷誤差

5 模型試算

5.1 均勻介質(zhì)模型

圖2波場快照所選取均勻介質(zhì)模型的聲波速度設(shè)定為4 000 m/s，空間步長均取為25 m，激發(fā)位置位于模型的中心處，震源為20 Hz雷克子波。圖中x軸表示為模型的橫向長度，y軸為縱向長度。小節(jié)將三點五階迎風(fēng)超緊致格式與另外3種方法對均勻介質(zhì)進(jìn)行數(shù)值模擬后的結(jié)果進(jìn)行對比，圖2的時間步長Δt=0.001 s，空間步長則盡量取較大的數(shù)值以對比4種方法哪種在頻散壓制方面具有優(yōu)勢。圖2的均勻介質(zhì)模型的長度及深度均取4.8 km這個較大數(shù)值。

圖2 空間步長為25 m,400 ms時各方法波場快照

從圖2可知普通緊致差分所模擬得到的結(jié)果存在頻散現(xiàn)象。當(dāng)空間步長取到25 m時，七點六階差分格式與普通六階緊致差分格式均有強烈的頻散現(xiàn)象，三點六階超緊致差分格式有輕微的頻散現(xiàn)象。而三點五階迎風(fēng)超緊致差分格式的結(jié)果幾乎沒有頻散現(xiàn)象；將上述幾種方法分別對應(yīng)比較可知，三點五階迎風(fēng)超緊致差分格式優(yōu)勢相當(dāng)明顯。

5.2 水平介質(zhì)模型

水平介質(zhì)模型設(shè)置為兩層，每層的厚度為2 km，長度設(shè)置與深度相同取4 km。其中水平介質(zhì)模型速度設(shè)置為：第1層波速為4 000 m/s；第2層波速設(shè)為4 500 m/s。以地表中心進(jìn)行激發(fā)，震源為20 Hz雷克子波，時間與空間步長分別取0.001 s與10 m。將迎風(fēng)超緊致差分格式方法應(yīng)用于該模型中進(jìn)行數(shù)值模擬。結(jié)果如圖3所示，圖中x軸為模型橫向長度，y軸為聲波傳播時間。

圖4 水平層狀介質(zhì)模型聲波方程數(shù)值模擬地面地震記錄

圖3中波場快照清晰得表現(xiàn)出聲波在水平層狀介質(zhì)中傳播效果，沒有數(shù)值頻散和不穩(wěn)定數(shù)值結(jié)果，圖3中所得到的地面地震記錄反映清晰，直達(dá)波、反射波顯示清楚，說明該方法能夠有效而清晰地模擬聲波在多層各向同性介質(zhì)模型中的傳播。

6 結(jié) 論

從地震正演模擬的聲波方程數(shù)值模擬入手，簡單介紹了傳統(tǒng)的有限差分，由于精度的要求問題，往往不能滿足需求，因此引出了三點五階迎風(fēng)超緊致差分方法，也是本文主要介紹的內(nèi)容，它具有精度高、分辨率高、高穩(wěn)定性、網(wǎng)格節(jié)點數(shù)需求較少等優(yōu)點。在介紹了三點五階超緊致差分格式的原理后，對比分析了該差分格式的截斷誤差、穩(wěn)定性和頻散等問題，隨后做了均勻介質(zhì)的聲波方程數(shù)值模擬，并與其它方法進(jìn)行對比分析，得出的結(jié)果是該方法不僅在精度上有著較為明顯的優(yōu)勢，而且對頻散的壓制方面也有很好的效果。最后，將這種超緊致差分格式應(yīng)用于對層狀介質(zhì)的聲波方程數(shù)值模擬中，取得了良好的效果，表明該方法能夠有效地模擬聲波在多層各向同性介質(zhì)模型中的傳播。相比原來的三點六階超緊致差分格式，三點五階迎風(fēng)超緊致差分格式的穩(wěn)定性有明顯的提高。

本文在進(jìn)行聲波方程數(shù)值模擬時，由于該方法使用的是均勻網(wǎng)格，而且對于地震正演模擬來說，只是對均勻介質(zhì)及水平層狀介質(zhì)的模擬還是不夠的，迎風(fēng)超緊致格式的優(yōu)點沒有完全得體現(xiàn)出來。今后的研究中，可以應(yīng)用于其他非均勻的網(wǎng)格中，比如交錯網(wǎng)格，亦可以結(jié)合其他復(fù)雜介質(zhì)進(jìn)行地震正演模擬分析。

參考文獻(xiàn)：

[1] Komatissch D, Barnes C, Tromp J. Simulation of anisotropic wave propagation based upon a spectral element method[J]. Geophysics, 2000, 65(4): 1251-1260.

[2] Fu D X, Ma Y W, Liu H. Proceedings, 5th International symposium on CFD. Sendai. 1993 1:104.

[3] Fu D X, Ma Y W. A high order accuracy difference scheme for complex flow fields[J]. J ComputPhys,1997(134): 1-15.

[4] 王芳芳. 緊致差分格式的理論與分析[D]. 遼寧：東北大學(xué)理學(xué)院，2010.

[5] Lele S K. Compact finite difference scheme with spectral-like resolution[J]. Journal of Computational Physics, 1992,103(1):16-42.

[6] 林東, 詹杰民. 淺水方程組合型緊致差分格式[J]. 計算力學(xué)學(xué)報, 2008, 25(6):791-796.

[7] 王書強，楊頂輝，楊寬德. 彈性波方程的緊致差分方法[J]. 清華大學(xué)學(xué)報(自然科學(xué)版)，2002，4(8)：1128-1131.

[8] 梁賢，田振夫. 求解Navier-Stokes 方程組的組合緊致迎風(fēng)格式[J]. 計算物理，2008,25(6)：660-665.

[9] Krishnan Mahesh. A family of high order finite difference schemes with good spectral resolution[J]. Journal of Computational Physics, 1998(145): 332-358.

[10] Ehaterinaris J A. Implicit, high-resolution, compact schemes for gas dynamics and aeroacoustics[J]. Journal of Computational Physics, 1999,156(2):272-299.

[11] Blayo E. Compact finite difference schemes for ocean models[J]. Journal of Computational Physics, 2000,164(2):241-257.

[12] Sengupta T K, Ganeriwal G, De S. Analysis of central and upwind compact schemes[J]. Journal of Computational Physics,2003,192(2):677-694.