基于時空域交錯網(wǎng)格有限差分法的應(yīng)力速度聲波方程數(shù)值模擬

彭更新,劉 威,郭念民,胡自多,徐凱馳,裴廣平

(1.中國石油天然氣集團公司塔里木油田分公司,新疆庫爾勒841000;2.中國石油勘探開發(fā)研究院西北分院,甘肅蘭州730020)

有限差分法具有占用內(nèi)存小和編程實現(xiàn)簡單等優(yōu)點,是應(yīng)用最廣泛的一種波動方程數(shù)值模擬方法[1-3]。有限差分法利用差分算子近似表示波動方程中的時間和空間偏微分算子,進(jìn)而通過迭代求解差分離散波動方程實現(xiàn)波動方程數(shù)值模擬。差分算子近似微分算子,會導(dǎo)致地震波的傳播速度與真實速度不相等,并且不同頻率成分的傳播速度不相同,即出現(xiàn)數(shù)值頻散現(xiàn)象[4-5]。固有的數(shù)值頻散嚴(yán)重影響有限差分法的數(shù)值模擬精度;同時,相比基于射線理論[6]和高斯束理論[7]的數(shù)值模擬方法,有限差分法的計算效率相對較低。因此,提高模擬精度,同時保持甚至提高計算效率是有限差分?jǐn)?shù)值模擬方法面臨的重要挑戰(zhàn)。

波動方程數(shù)值模擬廣泛采用的有限差分法可以分成兩類,即常規(guī)網(wǎng)格有限差分法和交錯網(wǎng)格有限差分法。常規(guī)網(wǎng)格有限差分法主要應(yīng)用于二階標(biāo)量波動方程數(shù)值模擬[8-10];交錯網(wǎng)格有限差分法主要應(yīng)用于一階應(yīng)力-速度聲波[11-13]或彈性波動方程數(shù)值模擬[14-16]。相比常規(guī)網(wǎng)格差分法,交錯網(wǎng)格差分法通常具有更高的模擬精度,并且對速度對比度(介質(zhì)最大速度與最小速度之比)較大的介質(zhì)數(shù)值模擬穩(wěn)定性更高[17],同時交錯網(wǎng)格差分法更便于施加應(yīng)力和位移震源。

KANE[18]首次將交錯網(wǎng)格有限差分應(yīng)用于電磁波方程數(shù)值模擬,VIRIEUX[14-15]將時間二階和空間二階交錯網(wǎng)格有限差分法應(yīng)用于地震波場的數(shù)值模擬。為了提高交錯網(wǎng)格差分法的模擬精度,LEVANDER[19]在VIRIEUX[14-15]研究成果的基礎(chǔ)上,進(jìn)一步發(fā)展了時間二階和空間四階交錯網(wǎng)格差分法。通過增大空間差分算子長度構(gòu)建出空間任意偶數(shù)階交錯網(wǎng)格差分法,相應(yīng)的差分系數(shù)通?；诳臻g域頻散關(guān)系和泰勒展開求解[20]。本文將這種時間二階、空間2M階并基于空間域頻散關(guān)系和泰勒展開計算差分系數(shù)的交錯網(wǎng)格差分法稱為空間域交錯網(wǎng)格有限差分法(SD-SFD)。差分離散波動方程在時空域同時迭代求解,而空間域交錯網(wǎng)格有限差分法僅基于空間域頻散關(guān)系計算差分系數(shù),導(dǎo)致空間域交錯網(wǎng)格有限差分法容易出現(xiàn)數(shù)值頻散,模擬精度低。

本文旨在改進(jìn)空間域交錯網(wǎng)格有限差分法的差分系數(shù)算法,進(jìn)而提高模擬精度。在保持時間二階、空間2M階差分格式不變的情況下,提出利用時空域頻散關(guān)系和泰勒展開計算差分系數(shù),稱這種改進(jìn)差分系數(shù)算法的交錯網(wǎng)格差分法為時空域交錯網(wǎng)格有限差分法(TSD-SFD)。首先,給出時間二階、空間2M階交錯網(wǎng)格差分法對二維應(yīng)力-速度聲波方程的差分離散方程;其次,推導(dǎo)基于時空域頻散關(guān)系和泰勒展開的差分系數(shù)算法;再次,分析和對比時空域交錯網(wǎng)格有限差分法和空間域交錯網(wǎng)格有限差分法的數(shù)值頻散特性和穩(wěn)定性特性;最后,利用層狀介質(zhì)模型和塔里木盆地典型復(fù)雜構(gòu)造模型進(jìn)行數(shù)值模擬,對比兩種交錯網(wǎng)格差分法的模擬精度和計算效率。

1 交錯網(wǎng)格有限差分法

二維應(yīng)力速度聲波方程可以表示為:

(1)

式中:P=P(x,z,t)為壓力場,u=u(x,z,t)和w=w(x,z,t)分別為質(zhì)點震動速度場的x和z分量,κ為體積模量,ρ為介質(zhì)密度。

應(yīng)力速度聲波方程數(shù)值模擬普遍采用交錯網(wǎng)格有限差分法。交錯網(wǎng)格有限差分法的基本思想是將壓力場P、質(zhì)點震動速度場u和w定義在交錯的網(wǎng)格系統(tǒng)中,如圖1所示。

圖1 二維交錯網(wǎng)格差分法示意

交錯網(wǎng)格有限差分法通常采用時間二階、空間2M階差分法。對波動方程中的一階時間偏導(dǎo)數(shù)?P/?t、?u/?t和?w/?t采用二階差分近似,可以得到:

(2)

(3)

(4)

對波動方程中的一階空間偏導(dǎo)數(shù)采用2M階差分近似,可以得到:

(5)

(6)

(7)

(8)

式中:am(m=1,2,…,M)為差分系數(shù)。

將方程(2)～(8)代入方程(1),可以得到:

(9)

方程(9)為時間二階、空間2M階交錯網(wǎng)格差分法對應(yīng)力速度聲波方程的差分離散波動方程?？臻g域交錯網(wǎng)格有限差分法和時空域交錯網(wǎng)格有限差分法均通過迭代求解方程(9)實現(xiàn)波動方程數(shù)值模擬,只是二者采用的差分系數(shù)算法不同。

2 差分系數(shù)計算

差分系數(shù)直接影響交錯網(wǎng)格有限差分法的數(shù)值模擬精度和穩(wěn)定性,因此,差分系數(shù)算法是交錯網(wǎng)格有限差分法的重要研究內(nèi)容。

空間域交錯網(wǎng)格有限差分法僅在空間域計算差分系數(shù),相應(yīng)的差分系數(shù)算法可概括為:首先,根據(jù)平面波理論和空間差分算子表達(dá)式推導(dǎo)空間域頻散關(guān)系;其次,對空間域頻散關(guān)系中的三角函數(shù)進(jìn)行泰勒級數(shù)展開,建立差分系數(shù)求解方程組,再通過解方程組計算差分系數(shù)。

時空域交錯網(wǎng)格有限差分法在時空域計算差分系數(shù),相應(yīng)的差分系數(shù)算法可概括為:首先,根據(jù)平面波理論和差分離散波動方程推導(dǎo)時空域頻散關(guān)系;其次,對時空域頻散關(guān)系中的三角函數(shù)進(jìn)行泰勒級數(shù)展開,建立差分系數(shù)求解方程組,通過解方程組計算差分系數(shù)。

空間域交錯網(wǎng)格有限差分法和時空域交錯網(wǎng)格有限差分法的差分系數(shù)計算過程均采用了平面波理論。在均勻介質(zhì)中,應(yīng)力速度聲波方程具有如下形式的離散平面波解:

(10)

(11)

(12)

kx=kcosθkz=ksinθ

(13)

式中:AP、Au和Aw為平面波的振幅,k為波數(shù),θ為平面波傳播方向與x軸的夾角,ω為角頻率。

2.1 空間域差分系數(shù)算法

空間域交錯網(wǎng)格有限差分法僅在空間域計算差分系數(shù),將離散平面波解方程(10)代入空間差分算子方程(5)得到:

(14)

方程(14)為空間2M階交錯網(wǎng)格差分法的空間域頻散關(guān)系,泰勒展開其中的正弦函數(shù)得到:

(15)

(16)

解方程(16)得到:

(17)

m=1,2,…,M

方程(17)為空間域交錯網(wǎng)格有限差分法的差分系數(shù)通解,可以看出其差分系數(shù)僅與M的取值有關(guān),而與地震波的傳播速度無關(guān)。

2.2 時空域差分系數(shù)算法

時空域交錯網(wǎng)格有限差分法在時空域計算差分系數(shù),將平面波解方程(10)～(13)代入差分離散波動方程(9)得到:

(18a)

(18b)

(18c)

消除方程18a～18c中的AP,Au和Aw,并且考慮到ω=υk,κ=ρυ2得到:

(19)

式中:r=υτ/h為Courant條件數(shù),表示單位時間步長內(nèi)地震波的傳播距離與空間采樣步長之比。

方程(19)為時間二階和空間2M階交錯網(wǎng)格差分法的時空域頻散關(guān)系。利用時空域頻散關(guān)系和泰勒展開求解差分系數(shù),需要選擇一個特定θ值。為了盡可能簡化差分系數(shù)計算過程,取θ=0或θ=π/2,可以得到:

(20)

泰勒展開式(20)中的正弦函數(shù)得到:

(21)

使方程(21)左右兩邊k2j-1h2j-2的系數(shù)對應(yīng)相等得到:

(22)

j=1,2,…,M

解方程(22)得到:

(23)

m=1,2,…,M

方程(23)為時空域交錯網(wǎng)格有限差分法的差分系數(shù)通解。取r=0時,方程(23)與方程(17)完全等價,因此,空間域交錯網(wǎng)格有限差分法的差分系數(shù)是時空域交錯網(wǎng)格有限差分法的差分系數(shù)的一個特解。同時,方程(23)表明時空域交錯網(wǎng)格有限差分法的差分系數(shù)不僅與M的取值相關(guān),還與r=υτ/h取值相關(guān),數(shù)值模擬過程中時間步長τ和空間步長h取值固定,差分系數(shù)隨速度υ自適應(yīng)變化,這是時空域交錯網(wǎng)格有限差分法比空間域交錯網(wǎng)格有限差分法具有更高模擬精度的根本原因。

3 數(shù)值頻散分析

數(shù)值頻散的嚴(yán)重程度直接影響交錯網(wǎng)格有限差分法的數(shù)值模擬精度。本文采用歸一化相速度誤差εph(θ)度量數(shù)值頻散的大小。根據(jù)時間二階、空間2M階交錯網(wǎng)格差分法的時空域頻散關(guān)系方程(19),可導(dǎo)出εph(θ)的表達(dá)式為:

(24)

(25)

當(dāng)εph(θ)=0時,相速度等于真實速度,無數(shù)值頻散;εph(θ)>0時,相速度大于真實速度,表現(xiàn)出時間頻散;εph(θ)<0時,相速度小于真實速度,表現(xiàn)出空間頻散。

空間域交錯網(wǎng)格有限差分法和時空域交錯網(wǎng)格有限差分法的歸一化相速度誤差εph(θ)的表達(dá)式完全相同,均由方程(24)和(25)表示,但是二者的差分系數(shù)am(m=1,2,…,M)不同,從而表現(xiàn)出不同的數(shù)值頻散特性。

根據(jù)εph(θ)的方程(24)和方程(25),結(jié)合兩種差分法各自的差分系數(shù),可繪制兩種差分法的相速度頻散曲線,進(jìn)而分析數(shù)值頻散特性。

圖2給出了r=0.2和r=0.3時,空間域交錯網(wǎng)格有限差分法(M=10)和時空域交錯網(wǎng)格有限差分法(M=10)的相速度數(shù)值頻散曲線,分析對比可以看出:

圖2 相速度數(shù)值頻散曲線

1)r取值增大,空間域交錯網(wǎng)格有限差分法(M=10)和時空域交錯網(wǎng)格有限差分法(M=10)的相速度數(shù)值頻散均明顯增大;

2)在波動方程有限差分法數(shù)值模擬中,通常保證單位波長內(nèi)有4個采樣點,對應(yīng)kh/π的取值為0.5,在0≤kh/π≤0.5范圍內(nèi),r取值相同時,時空域交錯網(wǎng)格有限差分法(M=10)的相速度數(shù)值頻散明顯小于空間域交錯網(wǎng)格有限差分法(M=10),因此,時空域交錯網(wǎng)格有限差分法(M=10)壓制數(shù)值頻散的效果更好,模擬精度更高;

3)在0≤kh/π≤0.5范圍內(nèi),空間域交錯網(wǎng)格有限差分法(M=10,r=0.2)和時空域交錯網(wǎng)格有限差分法(M=10,r=0.3)的相速度數(shù)值頻散大小基本相當(dāng),但時空域交錯網(wǎng)格有限差分法(M=10)采用的r取值更大,意味著能采用更大的時間步長,從而提高計算效率,同時不損失模擬精度;

4)當(dāng)kh/π≥0.6時,時空域交錯網(wǎng)格有限差分法(M=10)存在一定的空間數(shù)值頻散,這時相比空間域交錯網(wǎng)格有限差分法在壓制數(shù)值頻散,提高模擬精度方面的優(yōu)勢會降低。

4 穩(wěn)定性分析

根據(jù)時間二階、空間2M階交錯網(wǎng)格差分法的時空域頻散方程得到:

(26)

取空間波數(shù)kx=kz=π/h,并且根據(jù)三角函數(shù)有界性sin2(rkh/2)≤1,可得到:

(27)

式中:S為穩(wěn)定性因子。

空間域交錯網(wǎng)格有限差分法和時空域交錯網(wǎng)格有限差分法的穩(wěn)定性條件均由方程(27)所描述,但兩種方法的差分系數(shù)不同。

圖3a和圖3b分別給出了兩種差分法的穩(wěn)定性因子S隨r取值變化的曲線。空間域的差分系數(shù)與r取值無關(guān),其穩(wěn)定性因子S不隨r變化,呈現(xiàn)為一條平行于r坐標(biāo)軸的直線;時空域的差分系數(shù)為r的函數(shù),其穩(wěn)定性因子S表現(xiàn)為隨r取值變化的曲線。

根據(jù)穩(wěn)定性條件不等式r≤S,穩(wěn)定性因子S曲線(直線)與直線S=r的交點為穩(wěn)定性條件約束下的最大r取值。圖3c給出了兩種差分法穩(wěn)定性條件不等式約束下的最大r取值隨M的變化曲線,稱為穩(wěn)定性曲線。穩(wěn)定性條件不等式約束下的最大r取值越大,穩(wěn)定性越強,反之,則穩(wěn)定性越弱。圖3c表明,空間域交錯網(wǎng)格有限差分法和時空域交錯網(wǎng)格有限差分法的穩(wěn)定性均隨M增大而降低;當(dāng)M≥1時,時空域交錯網(wǎng)格有限差分法的穩(wěn)定性明顯強于空間域交錯網(wǎng)格有限差分法,因此利用時空域交錯網(wǎng)格有限差分法進(jìn)行數(shù)值模擬時,能采用更大的r取值,即采用更大的時間步長,從而提高計算效率。

圖3 穩(wěn)定性曲線

5 數(shù)值模擬實例

利用層狀介質(zhì)模型和塔里木盆地典型復(fù)雜構(gòu)造模型開展數(shù)值模擬實驗,進(jìn)一步分析、對比兩種方法的數(shù)值模擬精度和計算效率。

5.1 層狀介質(zhì)模型

圖4給出了一個網(wǎng)格尺寸為Δx=Δz=h=15 m,網(wǎng)格數(shù)為nz×nx=601×801的層狀介質(zhì)模型。主頻為25 Hz的Ricker子波作為震源,位于網(wǎng)格點(401,3)。空間域交錯網(wǎng)格有限差分法(M=10)和時空域交錯網(wǎng)格有限差分法(M=10)分別采用時間步長τ=1.5 ms和τ=1.0 ms進(jìn)行數(shù)值模擬。圖5和圖6分別給出了利用兩種差分法(M=10)進(jìn)行數(shù)值模擬得到的2.1 s時刻壓力場P和質(zhì)點震動速度場x分量u的波場快照。

圖4 層狀介質(zhì)速度模型

圖5 層狀介質(zhì)模型數(shù)值模擬2.1 s時刻波場快照(時間步長τ=1.5 ms)

圖6 層狀介質(zhì)模型數(shù)值模擬2.1 s時刻波場快照(時間步長τ=1.0 ms)

對比兩組波場快照,可以看出:采用時間步長τ=1.5 ms時,空間域交錯網(wǎng)格有限差分法(M=10)存在嚴(yán)重的數(shù)值頻散,而時空域交錯網(wǎng)格有限差分法(M=10)僅表現(xiàn)輕微的數(shù)值頻散;減小時間步長至τ=1.0 ms,空間域交錯網(wǎng)格有限差分法(M=10)的數(shù)值頻散明顯減弱,但仍存在較明顯的數(shù)值頻散,時空域交錯網(wǎng)格有限差分法(M=10)的數(shù)值頻散進(jìn)一步減弱,無明顯的數(shù)值頻散。因此,采用相同時間步長時,即計算效率基本相同時,時空域交錯網(wǎng)格有限差分法比空間域交錯網(wǎng)格有限差分法能更有效地壓制數(shù)值頻散,獲得更高的模擬精度。

進(jìn)一步分析可以看出,時空域交錯網(wǎng)格有限差分法(M=10)采用時間步長τ=1.5 ms比空間域交錯網(wǎng)格有限差分法(M=10)采用時間步長τ=1.0 ms數(shù)值頻散更弱,模擬精度更高。因此,相比空間域交錯網(wǎng)格有限差分法,時空域交錯網(wǎng)格有限差分法能采用更大的時間步長以提高計算效率,同時達(dá)到更高的模擬精度。

5.2 塔里木盆地典型復(fù)雜構(gòu)造模型

圖7為塔里木盆地典型復(fù)雜構(gòu)造模型,模型網(wǎng)格規(guī)模為nz×nx=1 051×1 201,網(wǎng)格尺寸為Δx=Δz=h=15 m。主頻25 Hz的Ricker子波作為震源,位于網(wǎng)格點(501,3)。兩種差分法均采用時間步長τ=1.0 ms進(jìn)行數(shù)值模擬。圖8為兩種差分法數(shù)值模擬3.0 s時刻壓力場P的波場快照,圖9為兩種差分法數(shù)值模擬得到的壓力場P的炮集記錄。

圖7 塔里木盆地典型復(fù)雜構(gòu)造速度模型

圖8 塔里木盆地復(fù)雜構(gòu)造模型數(shù)值模擬3.0 s時刻壓力場P的波場快照(時間步長τ=1.0 ms)

圖9 塔里木盆地復(fù)雜構(gòu)造模型數(shù)值模擬壓力場P的炮集記錄(時間步長τ=1.0 ms)

對比圖8和圖9可以看出,空間域交錯網(wǎng)格有限差分法(M=10)數(shù)值模擬得到的波場快照和炮集記錄均表現(xiàn)出明顯的數(shù)值頻散,模擬精度低;時空域交錯網(wǎng)格有限差分法(M=10)數(shù)值模擬的波場快照和炮集記錄均無明顯的數(shù)值頻散,模擬精度高。塔里木盆地典型復(fù)雜構(gòu)造模型模擬結(jié)果進(jìn)一步證明了時空域交錯網(wǎng)格有限差分法能更有效地壓制數(shù)值頻散,從而提高模擬精度。

6 結(jié)論

鑒于波動方程有限差分?jǐn)?shù)值模擬在時間域和空間域迭代求解,本文針對時間二階、空間2M階交錯網(wǎng)格差分法,推導(dǎo)出基于時空域頻散關(guān)系和泰勒展開的差分系數(shù)算法,提出時空域交錯網(wǎng)格有限差分法,并進(jìn)行數(shù)值頻散分析、穩(wěn)定性分析和數(shù)值模擬實驗,結(jié)果如下。

1)相比空間域交錯網(wǎng)格有限差分法,時空域交錯網(wǎng)格有限差分法能更有效地壓制數(shù)值頻散,模擬精度更高;同時,時空域交錯網(wǎng)格有限差分法還具有更強的穩(wěn)定性,能夠采用更大的時間步長提高計算效率。

2)層狀介質(zhì)模型和塔里木盆地典型復(fù)雜構(gòu)造模型數(shù)值模擬結(jié)果進(jìn)一步驗證了時空域交錯網(wǎng)格有限差分法在提高模擬精度和計算效率方法面的優(yōu)越性,也驗證了該方法的普遍適用性。

本文提出時空域交錯網(wǎng)格有限差分法僅針對二維應(yīng)力速度聲波方程開展了研究,仍存在以下不足。

1)文中僅開展二維應(yīng)力速度聲波方程數(shù)值模擬,針對三維情況有待進(jìn)一步研究。

2)相比空間域交錯網(wǎng)格有限差分法,時空域交錯網(wǎng)格有限差分格式的數(shù)值頻散明顯減小,但是其頻散曲線仍然較發(fā)散,下一步研究通過構(gòu)建更合理的交錯網(wǎng)格差分法改善相速度頻散曲線的收斂性,從而進(jìn)一步提高其模擬精度。

3)文中未討論時空域交錯網(wǎng)格有限差分法對彈性波方程的適用性,將進(jìn)一步開展深入研究。