三維波動(dòng)方程時(shí)空域混合網(wǎng)格有限差分?jǐn)?shù)值模擬方法

胡自多， 劉威， 雍學(xué)善， 王小衛(wèi)， 韓令賀， 田彥燦

1 中國(guó)石油勘探開(kāi)發(fā)研究院西北分院， 蘭州 730020 2 中國(guó)石油天然氣集團(tuán)有限公司油藏描述重點(diǎn)實(shí)驗(yàn)室， 蘭州 730020

0 引言

波動(dòng)方程數(shù)值模擬是逆時(shí)偏移(Baysal et al., 1983; Virieux et al., 2011; 陳生昌和周華敏, 2018)和全波形反演(Tarantola, 1984; Pratt et al., 1998; 董良國(guó)等, 2015)的重要基礎(chǔ).有限元法(Marfurt, 1984; Moczo et al., 2010, 2011)、虛譜法(Reshef et al., 1988; 黃建平等, 2016)和有限差分法(Alterman and Karal, 1968; Liu and Sen, 2011a; 梁文全等, 2013)是求解波動(dòng)方程的三大主流數(shù)值算法(皮紅梅等，2009；胡恒山，2018).相比其他兩類(lèi)方法，有限差分法占用內(nèi)存更小，計(jì)算效率更高，因而成為應(yīng)用最普遍的波動(dòng)方程數(shù)值模擬算法(蔣韜等，2008；張志禹等，2017；姜占東等，2021).但是，有限差分法采用差分算子近似波動(dòng)方程中的時(shí)間和空間偏微分算子，導(dǎo)致時(shí)間和空間數(shù)值頻散，降低了模擬精度(Alford et al., 1974; Dablain, 1986)；同時(shí)，逆時(shí)偏移和全波形反演巨大的計(jì)算量，發(fā)展高精度和高效率的有限差分模擬算法具有重要意義.

Dablain(1986)研究指出，時(shí)間和空間高階差分格式能夠有效壓制數(shù)值頻散、提高模擬精度.然而，時(shí)間高階差分格式占用內(nèi)存過(guò)大，穩(wěn)定性降低(Chen, 2007, 2011).針對(duì)標(biāo)量波動(dòng)方程，常規(guī)高階差分方法采用時(shí)間2階和空間高階(2M)差分格式，并利用空間域頻散關(guān)系和泰勒級(jí)數(shù)展開(kāi)求解差分系數(shù)，本文稱(chēng)為常規(guī)空間域高階有限差分法(簡(jiǎn)稱(chēng)為CSD-FDM).CSD-FDM的差分系數(shù)算法僅衡量了空間差分算子(Laplace差分算子)的精度，而沒(méi)有考慮差分離散波動(dòng)方程的精度，空間差分算子雖然能達(dá)到2M階精度，但差分離散波動(dòng)方程僅具有2階精度.實(shí)際上，有限差分法通過(guò)迭代求解差分離散波動(dòng)方程實(shí)現(xiàn)波動(dòng)方程數(shù)值模擬，差分離散波動(dòng)方程的差分精度才能更準(zhǔn)確地描述有限差分法的模擬精度.Liu和Sen (2009)保持CSD-FDM的差分格式不變，提出利用時(shí)空域頻散關(guān)系改進(jìn)差分系數(shù)算法，本文稱(chēng)為時(shí)空域高階有限差分法(簡(jiǎn)稱(chēng)為T(mén)SD-FDM).TSD-FDM比CSD-FDM具有更高的模擬精度和更強(qiáng)的穩(wěn)定性.TSD-FDM在二維和三維標(biāo)量波動(dòng)方程模擬中，差分離散波動(dòng)方程分別沿8個(gè)和48個(gè)傳播方向具有2M階差分精度，其他傳播方向只能達(dá)到2階差分精度，存在數(shù)值各向異性.基于相同思路，Liu和Sen(2011b)提出了面向應(yīng)力-速度聲波方程的時(shí)空域高階交錯(cuò)網(wǎng)格差分法，同樣有效提高了模擬精度和穩(wěn)定性.嚴(yán)紅勇等進(jìn)一步推廣應(yīng)用于聲波、VTI介質(zhì)、黏聲介質(zhì)逆時(shí)偏移中，有效壓制了數(shù)值頻散造成的成像假象(嚴(yán)紅勇和劉洋, 2013；Yan and Liu，2013a，b).基于泰勒級(jí)數(shù)展開(kāi)的差分系數(shù)算法計(jì)算效率高，通常低波數(shù)成分具有較高的模擬精度，但高波數(shù)模擬精度迅速降低.為了提高高波數(shù)成分的模擬精度，Zhang和Yao (2012,2013)、Liu(2013,2014)通過(guò)最小二乘優(yōu)化頻散誤差，改進(jìn)差分系數(shù)算法提高模擬精度.

除了改進(jìn)差分系數(shù)算法，設(shè)計(jì)更合理的差分格式是提高模擬精度的另一重要途徑.針對(duì)二維和三維應(yīng)力-速度聲波方程，Tan和Huang(2014a)利用坐標(biāo)軸網(wǎng)格點(diǎn)和非坐標(biāo)軸網(wǎng)格點(diǎn)構(gòu)建空間差分算子近似一階空間微分算子，提出了混合交錯(cuò)網(wǎng)格差分方法，并利用時(shí)空域頻散關(guān)系和泰勒展開(kāi)求解差分系數(shù)，差分離散波動(dòng)方程沿任意傳播方向可達(dá)到4階、6階差分精度，并具有更強(qiáng)的穩(wěn)定性.Tan和Huang(2014b)利用最小二乘法優(yōu)化差分系數(shù)，模擬精度進(jìn)一步提高.Ren和Li(2017)針對(duì)二維彈性波方程，構(gòu)建了更具一般性的混合交錯(cuò)網(wǎng)格差分格式，彈性波差分離散方程沿任意傳播方向可達(dá)到2N(N<5)階差分精度.混合交錯(cuò)網(wǎng)格差分格式，應(yīng)力和速度的空間一階偏導(dǎo)算子，各坐標(biāo)軸相互獨(dú)立計(jì)算，因此2D可較易推廣到3D.

針對(duì)二維標(biāo)量波動(dòng)方程，Liu和Sen(2013)、張保慶等(2016)、胡自多等(2016)和Wang等(2016)先后提出將二維Laplace微分算子近似為常規(guī)笛卡爾坐標(biāo)系中坐標(biāo)軸網(wǎng)格點(diǎn)構(gòu)建的M個(gè)Laplace差分算子和旋轉(zhuǎn)笛卡爾坐標(biāo)系中坐標(biāo)軸網(wǎng)格點(diǎn)構(gòu)建的N個(gè)Laplace差分算子的加權(quán)平均，構(gòu)建了不同的混合網(wǎng)格差分格式，利用時(shí)空域頻散關(guān)系和泰勒展開(kāi)求解差分系數(shù)，差分離散波動(dòng)方程沿任意傳播方向理論上可達(dá)到任意偶數(shù)階差分精度，明顯提高了二維標(biāo)量波動(dòng)方程有限差分法的模擬精度，穩(wěn)定性進(jìn)一步增強(qiáng)，還能采用更大的時(shí)間采樣間隔以獲得更高的計(jì)算效率.不同于混合交錯(cuò)網(wǎng)格差分格式，標(biāo)量波動(dòng)方程壓力對(duì)空間的2階偏導(dǎo)(Laplace算子)作為一個(gè)整體計(jì)算，2D混合網(wǎng)格差分格式不能直接推廣到3D.而地震勘探的對(duì)象通常是三維介質(zhì)，發(fā)展三維混合網(wǎng)格差分格式更具實(shí)際意義.

如何利用三維笛卡爾坐標(biāo)系中非坐標(biāo)軸網(wǎng)格點(diǎn)構(gòu)建Laplace差分算子是建立三維混合網(wǎng)格差分格式需要解決的關(guān)鍵問(wèn)題.本文將三維笛卡爾坐標(biāo)系中的非坐標(biāo)軸網(wǎng)格點(diǎn)分為兩類(lèi)：坐標(biāo)平面內(nèi)的非坐標(biāo)軸網(wǎng)格點(diǎn)和坐標(biāo)平面外的非坐標(biāo)軸網(wǎng)格點(diǎn)，并系統(tǒng)推導(dǎo)出了這兩類(lèi)非坐標(biāo)軸網(wǎng)格點(diǎn)構(gòu)建三維Laplace差分算子的方法.在此基礎(chǔ)上，將三維Laplace微分算子近似為坐標(biāo)軸網(wǎng)格點(diǎn)構(gòu)建的M個(gè)Laplace差分算子和非坐標(biāo)軸網(wǎng)格點(diǎn)構(gòu)建的N個(gè)Laplace差分算子的加權(quán)平均，提出了一種面向三維標(biāo)量波動(dòng)方程的混合網(wǎng)格M+N型差分方法(簡(jiǎn)稱(chēng)為MG-FDM)，并根據(jù)時(shí)空域頻散關(guān)系和泰勒級(jí)數(shù)展開(kāi)建立了差分系數(shù)求解方程組，推導(dǎo)出差分系數(shù)的通解.增大N的取值，三維MG-FDM給出的差分離散波動(dòng)方程沿任意傳播方向能夠達(dá)到4階、6階、甚至任意偶數(shù)階差分精度.頻散分析表明，相比三維CSD-FDM和TSD-FDM，計(jì)算效率基本相同時(shí)，MG-FDM的數(shù)值頻散壓制效果更好，模擬精度更高；模擬精度基本相當(dāng)時(shí)，MG-FDM能采用更大的時(shí)間采樣間隔，計(jì)算效率更高.穩(wěn)定性分析表明，MG-FDM穩(wěn)定性更強(qiáng)，為其采用更大的時(shí)間采樣間隔以提高計(jì)算效率奠定了基礎(chǔ).三維數(shù)值模擬實(shí)驗(yàn)進(jìn)一步證實(shí)了MG-FDM在提高模擬精度和計(jì)算效率方面的優(yōu)越性.

1 三維混合網(wǎng)格差分格式構(gòu)建方法

1.1 三維常規(guī)空間域高階(CSD-FDM)和時(shí)空域高階(TSD-FDM)差分格式

常密度介質(zhì)中，三維標(biāo)量波動(dòng)方程為：

(1)

有限差分法數(shù)值求解波動(dòng)方程普遍采用2階時(shí)間差分方案，式(1)右邊的時(shí)間偏微分可差分表示為：

(2)

三維常規(guī)空間域高階差分(CSD-FDM)和三維時(shí)空域高階差分(TSD-FDM)，均采用圖1所示的差分格式，波動(dòng)方程(1)中的三維Laplace微分算子可差分近似為：

(3)

圖1 三維CSD-FDM和TSD-FDM的差分格式Fig.1 FD stencil of 3D CSD-FDM and TSD-FDM

式(3)、(2)代入式(1)：

(4)

1.2 三維混合網(wǎng)格M+N型差分格式構(gòu)建

構(gòu)建三維混合網(wǎng)格M+N型差分格式(MG-FDM)的基本思想是聯(lián)合利用坐標(biāo)軸網(wǎng)格點(diǎn)和非坐標(biāo)軸網(wǎng)格點(diǎn)一起差分近似三維Laplace算子.如何利用非坐標(biāo)軸網(wǎng)格點(diǎn)構(gòu)建Laplace差分算子是建立三維MG-FDM需要解決的關(guān)鍵問(wèn)題.

圖2a—c給出了三維笛卡爾坐標(biāo)系中與差分中心點(diǎn)距離相等的三組非坐標(biāo)軸網(wǎng)格點(diǎn)，這三組網(wǎng)格點(diǎn)與差分中心點(diǎn)的距離依次增大.由于無(wú)法將任意一組非坐標(biāo)軸網(wǎng)格點(diǎn)置于坐標(biāo)原點(diǎn)位于差分中心點(diǎn)的三維旋轉(zhuǎn)笛卡爾坐標(biāo)系的坐標(biāo)軸上，導(dǎo)致不能直接利用非坐標(biāo)軸網(wǎng)格點(diǎn)構(gòu)建三維Laplace差分算子.為了利用非坐標(biāo)軸網(wǎng)格點(diǎn)構(gòu)建三維Laplace差分算子，本文將非坐標(biāo)軸網(wǎng)格點(diǎn)分成兩類(lèi)：坐標(biāo)平面內(nèi)的非坐標(biāo)軸網(wǎng)格點(diǎn)(如圖2a、c)和坐標(biāo)平面外的非坐標(biāo)軸網(wǎng)格點(diǎn)(如圖2b).

(5)

(6)

將三維坐標(biāo)系中坐標(biāo)軸網(wǎng)格點(diǎn)差分格式(圖1)與非坐標(biāo)軸網(wǎng)格點(diǎn)差分格式(圖2)進(jìn)行組合，構(gòu)建出三維MG-FDM的差分格式(圖3).圖1與圖2a中的差分格式組合，可構(gòu)建出三維MG-FDM(N=1)的差分格式(圖3a).式(3)與式(5)加權(quán)求和：

(7)

其中cm(m=1,2,…,M)和c1,1,0為權(quán)系數(shù).

圖2 三維坐標(biāo)系中非坐標(biāo)軸網(wǎng)格點(diǎn)構(gòu)建的差分格式(a) 坐標(biāo)平面內(nèi)的非坐標(biāo)軸網(wǎng)格點(diǎn)構(gòu)建的差分格式(網(wǎng)格點(diǎn)與中心點(diǎn)的距離是 坐標(biāo)平面外的非坐標(biāo)軸網(wǎng)格點(diǎn)構(gòu)建的差分格式(網(wǎng)格點(diǎn)與中心點(diǎn)的距離是 坐標(biāo)平面內(nèi)的非坐標(biāo)軸網(wǎng)格點(diǎn)構(gòu)建的差分格式(網(wǎng)格點(diǎn)與中心點(diǎn)的距離是Fig.2 FD stencils constructed by off-axial grid nodes in the 3D Cartesian coordinate system(a) FD stencil constructed by the off-axial grid nodes in the coordinate plane (distance between grid node and center grid node is (b) FD stencil constructed by the off-axial grid nodes outside the coordinate plane (distance between grid node and center grid node is (c) FD stencil constructed by the off-axial grid nodes in the coordinate plane (distance between grid node and center grid node is

圖3 三維MG-FDM的差分格式(a) MG-FDM(N=1)； (b) MG-FDM(N=2)； (c) MG-FDM(N=3).Fig.3 FD stencils of 3D MG-FDM

式(7)表明，三維MG-FDM(N=1)是將三維Laplace微分算子近似為坐標(biāo)軸網(wǎng)格點(diǎn)構(gòu)建的M個(gè)Laplace差分算子和非坐標(biāo)軸網(wǎng)格點(diǎn)構(gòu)建的1個(gè)Laplace差分算子的加權(quán)平均.式(7)和(2)代入(1)：

(8)

式(8)是三維MG-FDM(N=1)的差分離散波動(dòng)方程，附錄B給出了三維MG-FDM(N=2,3)對(duì)應(yīng)的差分離散波動(dòng)方程，同理可導(dǎo)出N取任意值時(shí)三維MG-FDM的差分離散波動(dòng)方程.MG-FDM實(shí)質(zhì)上是將三維Laplace微分算子近似為坐標(biāo)軸網(wǎng)格點(diǎn)構(gòu)建的M個(gè)Laplace差分算子和非坐標(biāo)軸網(wǎng)格點(diǎn)構(gòu)建的N個(gè)Laplace差分算子的加權(quán)平均.

2 差分系數(shù)計(jì)算方法和差分精度分析

三維CSD-FDM和TSD-FDM的差分系數(shù)計(jì)算，分別利用空間域頻散關(guān)系和時(shí)空域頻散關(guān)系.為了闡述三維MG-FDM差分系數(shù)計(jì)算方法和差分精度對(duì)比，我們首先列出三維CSD-FDM和TSD-FDM的差分系數(shù)算法，進(jìn)而推導(dǎo)出三維TSD-FDM和MG-FDM差分系數(shù)的通解.

2.1 三維CSD-FDM的差分系數(shù)算法

在均勻介質(zhì)中，三維標(biāo)量波動(dòng)方程具有如下離散形式的平面波解：

(9)

其中kx=ksinφcosθ,ky=ksinφsinθ,kz=kcosφ，ω為角頻率，k為波數(shù)，φ和θ為平面波的傳播方向角.

將式(9)代入式(3)中，得到CSD-FDM 的Laplace差分算子頻散關(guān)系，即空間域頻散關(guān)系：

+2cos(mkzh)-6].

(10)

對(duì)余弦函數(shù)進(jìn)行泰勒展開(kāi)，可導(dǎo)出三維CSD-FDM的差分系數(shù)通解為：

(m=1,2,…,M).

(11)

三維CSD-FDM差分系數(shù)計(jì)算，利用空間域頻散關(guān)系，僅考慮空間差分算子(Laplace差分算子)的差分精度.

2.2 三維TSD-FDM的差分系數(shù)算法

三維TSD-FDM差分系數(shù)計(jì)算，考慮了差分離散波動(dòng)方程的差分精度.將平面波解式(9)代入差分離散方程式(4)，得到差分離散波動(dòng)方程的頻散關(guān)系，即時(shí)空域頻散關(guān)系：

(12)

(13)

其中r=υτ/h為Courant條件數(shù)，表示單位時(shí)間采樣間隔(τ)內(nèi)波傳播的距離與空間采樣間隔(h)之比.

(m=1,2,…,M)，

(14)

其中f(j,φ,θ)=sin2jφ(cos2jθ+sin2jθ)+cos2jφ.可以看出，三維TSD-FDM的差分系數(shù)am(m=0,1,2,…,M)是地震波傳播方向角φ和θ的函數(shù).Liu和Sen(2009)計(jì)算差分系數(shù)時(shí)，取φ=π/2,θ=π/8，分析出三維TSD-FDM的差分離散波動(dòng)方程沿48個(gè)傳播方向具有2M階差分精度.

取φ=0,θ=0，式(14)可簡(jiǎn)化得到TSD-FDM差分系數(shù)通解的一個(gè)特殊形式：

(m=1,2,…,M)，

(15)

取φ=0,θ=0，三維TSD-FDM的差分離散波動(dòng)方程僅沿6個(gè)傳播方向(3個(gè)坐標(biāo)軸的正負(fù)方向)具有2M階差分精度.因此，取φ=π/2,θ=π/8計(jì)算的差分系數(shù)，TSD-FDM的差分離散波動(dòng)方程能達(dá)到更高的模擬精度.

2.3 三維MG-FDM的差分系數(shù)算法

借鑒時(shí)空域頻散關(guān)系和泰勒展開(kāi)求解差分系數(shù)的方法，本文推導(dǎo)出了MG-FDM的差分系數(shù)通解.

將平面波解式(9)代入三維MG-FDM(N=1)的差分離散波動(dòng)方程式(8)，得到：

+cos(kxh+kyh)+cos(kyh-kzh)+cos(kyh+kzh)+cos(kzh-kxh)+cos(kzh+kxh)-6]，

(16)

式(16)描述了三維MG-FDM(N=1)的差分離散波動(dòng)方程的時(shí)空域頻散關(guān)系，泰勒展開(kāi)余弦函數(shù)可得：

+(ky-kz)2j+(ky+kz)2j+(kz-kx)2j+(kz+kx)2j]h2j-2.

(17)

(j=1,2,…,M)，

(18)

(19)

2.4 差分精度分析和對(duì)比

差分精度是一種定量描述有限差分法模擬精度的傳統(tǒng)方法，通常將誤差函數(shù)關(guān)于時(shí)間采樣間隔τ或者空間采樣間隔h的最小冪指數(shù)稱(chēng)為差分精度階數(shù).根據(jù)頻散關(guān)系式(10)，定義三維CSD-FDM的Laplace差分算子的誤差函數(shù)εCSD-FDM為：

(20)

根據(jù)頻散關(guān)系式(12)，定義三維CSD-FDM的差分離散波動(dòng)方程的誤差函數(shù)ECSD-FDM為：

(21)

根據(jù)三維TSD-FDM差分離散波動(dòng)方程的頻散關(guān)系式(12)，結(jié)合式(13)，TSD-FDM的差分離散波動(dòng)方程的誤差函數(shù)：

(22)

計(jì)算差分系數(shù)時(shí)，取φ=π/2,θ=π/8，三維TSD-FDM的差分離散波動(dòng)方程，沿特定的48個(gè)傳播方向可達(dá)到2M階差分精度，但沿其他傳播方向僅具有2階差分精度.整體來(lái)講，三維TSD-FDM的差分離散波動(dòng)方程僅具有2階差分精度.

根據(jù)三維MG-FDM(N=1)的差分離散波動(dòng)方程的頻散關(guān)系式(16)，定義MG-FDM(N=1)的差分離散波動(dòng)方程的誤差函數(shù)EMG-FDM(N=1)為：

(23)

結(jié)合式(17)、(18)，誤差函數(shù)EMG-FDM(N=1)可表示為：

(24)

同理，可分析出三維MG-FDM(N=2,3)的差分精度.表1給出了三維CSD-FDM、TSD-FDM和MG-FDM(N=1,2,3)的差分精度，可以看出，相比三維CSD-FDM和TSD-FDM，MG-FDM能夠有效提高差分離散波動(dòng)方程的差分精度.增大N的取值，MG-FDM的差分離散波動(dòng)方程理論上可達(dá)到任意偶數(shù)階差分精度.

表1 三維CSD-FDM、TSD-FDM和MG-FDM(N=1,2,3)的差分離散波動(dòng)方程的差分精度Table 1 Difference accuracy of the discrete difference wave equation for 3D CSD-FDM, TSD-FDM and MG-FDM (N=1,2,3)

3 數(shù)值頻散分析和穩(wěn)定性分析

3.1 數(shù)值頻散分析

數(shù)值頻散使得相速度υph與地震波的真實(shí)傳播速度不相等，并且隨地震波的傳播方向變化，呈現(xiàn)數(shù)值各向異性特征.本文采用歸一化相速度δ(φ,θ)=υph/υ描述相速度的數(shù)值頻散特性，根據(jù)相速度定義υph=ω/k和三維MG-FDM(N=1)的差分離散波動(dòng)方程的頻散關(guān)系式(16)，可得出MG-FDM(N=1)的歸一化相速度δ(φ,θ)的表達(dá)式：

(25)

(26)

其中，G=λ/h，λ為波長(zhǎng)，G為波長(zhǎng)與空間采樣間隔之比，則G值越大，1/G越小，表示空間采樣越密，反之，空間采樣越稀疏.

歸一化相速度δ(φ,θ)的值越接近1，相速度數(shù)值頻散越??；δ(φ,θ)>1，相速度偏大，稱(chēng)為時(shí)間數(shù)值頻散；δ(φ,θ)<1，相速度偏小，稱(chēng)為空間數(shù)值頻散.另外,δ(φ,θ)隨地震波傳播方向角φ和θ變化越小，數(shù)值各向異性越弱，反之，數(shù)值各向異性越強(qiáng).相速度數(shù)值頻散越小，數(shù)值各向異性越弱，則數(shù)值模擬精度越高；相反，相速度數(shù)值頻散越大，數(shù)值各向異性越強(qiáng)，則數(shù)值模擬精度越低.

同樣地，可導(dǎo)出三維CSD-FDM、TSD-FDM和MG-FDM(N=2,3)的歸一化相速度δ(φ,θ)的表達(dá)式.δ(φ,θ)是傳播方向角φ和θ的函數(shù)，分析相速度頻散特性時(shí)給出了8個(gè)傳播角度(φ=π/2;θ=0,π/16,2π/16,3π/16,4π/16)和(φ=π/12,2π/12,3π/12;θ=2π/16)的相速度頻散曲線(xiàn).

圖4給出了三維CSD-FDM(M=2,6,10)、TSD-FDM(M=2,6,10)和MG-FDM(M=2,6,10;N=1)的相速度頻散曲線(xiàn)，Courant條件數(shù)r的取值均為r=0.3.分析對(duì)比相速度頻散曲線(xiàn)特征，可以得出：

(1)三維CSD-FDM(M=2)、TSD-FDM(M=2)和MG-FDM(M=2;N=1)的相速度頻散特征相似，均具有明顯的空間數(shù)值頻散，模擬精度都很低.

(2)三維CSD-FDM(M=6,10)具有明顯的時(shí)間數(shù)值頻散；TSD-FDM(M=6,10)的相速度頻散曲線(xiàn)較發(fā)散，同時(shí)存在空間和時(shí)間數(shù)值頻散，表現(xiàn)出明顯的數(shù)值各向異性(數(shù)值頻散特征隨地震波的傳播方向變化)；MG-FDM(M=6,10;N=1)的相速度頻散曲線(xiàn)收斂，數(shù)值各向異性明顯減弱，數(shù)值頻散幅值也明顯減小.

圖4 三維CSD-FDM、TSD-FDM和MG-FDM(N=1)的相速度頻散曲線(xiàn)(r=0.3)(a)(b)(c) CSD-FDM(M=2,6,10); (d)(e)(f) TSD-FDM(M=2,6,10); (g)(h)(i) MG-FDM(M=2,6,10;N=1)Fig.4 Phase velocity dispersion curves of 3D CSD-FDM, TSD-FDM and MG-FDM (N=1) with r=0.3

綜合(1)和(2)可以看出：相比三維CSD-FDM和TSD-FDM，M取值較小時(shí)(M=2左右)，MG-FDM在壓制相速度數(shù)值頻散方面無(wú)明顯優(yōu)勢(shì)；M取值較大時(shí)(M≥6)，MG-FDM的相速度數(shù)值頻散和數(shù)值各向異性均明顯減小，因而具有更高的模擬精度.

圖5給出了三維MG-FDM(M=8,15;N=1,2,3)的相速度頻散曲線(xiàn)(r=0.3).需要注意，圖5a、b、c和圖5d、e、f兩組頻散曲線(xiàn)具有不同的縱軸刻度.

圖5 三維MG-FDM的相速度頻散曲線(xiàn)(r=0.3)(a)(b)(c) MG-FDM(M=8;N=1,2,3)； (d)(e)(f) MG-FDM(M=15;N=1,2,3).Fig.5 Phase velocity dispersion curves of 3D MG-FDM with r=0.3

(1)三維MG-FDM(M=8;N=1,2,3)的相速度數(shù)值頻散特征基本相同；相比MG-FDM(M=15;N=1,2)，MG-FDM(M=15;N=3)的相速度頻散曲線(xiàn)更收斂，數(shù)值頻散幅值和數(shù)值各向異性均明顯減小.因此，三維標(biāo)量波動(dòng)方程數(shù)值模擬對(duì)精度要求較高時(shí)，建議采用MG-FDM(N=1)，且M取值為8左右；對(duì)模擬精度要求極其苛刻時(shí)，可采用MG-FDM(N=3)，且M取值為15左右；N取值大于4，同時(shí)M取值更大的三維MG-FDM因?yàn)閿?shù)值計(jì)算量巨大，計(jì)算效率極低而很少采用.三維MG-FDM可根據(jù)模擬精度要求選擇適當(dāng)?shù)腗和N值.

(2)三維MG-FDM(M=8;N=1,2)的相速度數(shù)值頻散特征基本相同，同樣地，MG-FDM(M=15;N=1,2)的相速度數(shù)值頻散特征也基本相同，和三維MG-FDM(N=1,2)的差分離散波動(dòng)方程均具有4階差分精度的結(jié)論一致(見(jiàn)表1).

圖6給出了r取值分別為r=0.1,0.2,0.4時(shí)，三維CSD-FDM(M=10)，TSD-FDM(M=10)和MG-FDM(M=8;N=1)的相速度頻散曲線(xiàn).根據(jù)r=υτ/h，空間采樣間隔h保持不變時(shí)，r取值從0.1增大至0.2和0.4，相應(yīng)的時(shí)間采樣間隔τ將分別增大至2倍和4倍.

(1)r取值為0.1，三維CSD-FDM表現(xiàn)出輕微的時(shí)間頻散，具有較高的模擬精度；r取值增大至0.2和0.4，時(shí)間頻散顯著增強(qiáng)，模擬精度低.

(2)r取值為0.1和0.2時(shí)，三維TSD-FDM的相速度頻散曲線(xiàn)較收斂，表現(xiàn)出輕微的空間和時(shí)間頻散，具有較高的模擬精度；r取值增大至0.4時(shí)，頻散曲線(xiàn)嚴(yán)重發(fā)散，數(shù)值頻散的幅值明顯增大，數(shù)值各向異性顯著增強(qiáng)，模擬精度低.

(3)當(dāng)r取值為0.1、0.2和0.4時(shí)，三維MG-FDM的相速度頻散曲線(xiàn)收斂性均較好，數(shù)值頻散微弱，模擬精度高.

圖6 三維CSD-FDM(M=10)，TSD-FDM(M=10)和MG-FDM(M=8;N=1)的相速度頻散曲線(xiàn)(a)(b)(c) 分別為r=0.1,0.2,0.4時(shí)CSD-FDM(M=10)的相速度頻散曲線(xiàn)； (d)(e)(f) 分別為r=0.1,0.2,0.4時(shí)TSD-FDM(M=10)的相速度頻散曲線(xiàn)； (g)(h)(i) 分別為r=0.1,0.2,0.4時(shí)MG-FDM(M=8;N=1)的相速度頻散曲線(xiàn).Fig.6 Phase velocity dispersion curves of 3D CSD-FDM (M=10), TSD-FDM (M=10) and MG-FDM (M=8; N=1)(a),(b) and (c) CSD-FDM (M=10) with r=0.1,0.2,0.4; (d),(e) and (f) TSD-FDM (M=10) with r=0.1,0.2,0.4; (g),(h) and (i) MG-FDM (M=8;N=1) with r=0.1,0.2,0.4.

綜合可知，三維MG-FDM(M=8;N=1)可以比CSD-FDM(M=10)和TSD-FDM(M=10)取更大的r值(速度模型υ和空間采樣間隔h相同時(shí)，r取值越大等價(jià)于時(shí)間采樣間隔τ越大)，同時(shí)數(shù)值頻散大小基本相當(dāng).因此，MG-FDM(M=8;N=1)能采用較大的時(shí)間采樣間隔以提高計(jì)算效率，同時(shí)保持較高的模擬精度.

3.2 穩(wěn)定性分析

有限差分法通過(guò)迭代求解差分離散波動(dòng)方程模擬地震波的傳播過(guò)程，必須確保迭代過(guò)程穩(wěn)定.根據(jù)式(16)可得出：

-cos(kxh+kyh)-cos(kyh-kzh)-cos(kyh+kzh)-cos(kzh-kxh)-cos(kzh+kxh)]，

(27)

取空間波數(shù)kx=ky=kz=π/h，并且0≤1-cos(ωτ)≤2，則可得到：

(28)

圖7給出了三維CSD-FDM、TSD-FDM和MG-FDM(N=1,2,3)的穩(wěn)定性曲線(xiàn)，描述了穩(wěn)定性條件約束下的最大r取值隨M的變化關(guān)系.穩(wěn)定性條件約束前提下，r取值越小，穩(wěn)定性越弱；相反，r取值越大，穩(wěn)定性越強(qiáng).分析圖中的穩(wěn)定性曲線(xiàn)可以看出：

(1)M的取值增大，三維CSD-FDM、TSD-FDM和MG-FDM(N=1,2,3)的穩(wěn)定性均下降.

(2)N的取值增大，MG-FDM的穩(wěn)定性增強(qiáng)，但差分離散波動(dòng)方程的差分精度相同時(shí)，穩(wěn)定性基本相同，MG-FDM(N=1,2)的差分離散波動(dòng)方程均為4階差分精度，穩(wěn)定性基本一致.

(3)三維CSD-FDM的穩(wěn)定性最弱，TSD-FDM的穩(wěn)定性明顯增強(qiáng)，MG-FDM(N=1,2,3)的穩(wěn)定性進(jìn)一步增強(qiáng).MG-FDM(N=1,2,3)的強(qiáng)穩(wěn)定性為數(shù)值模擬采用更大的時(shí)間采樣進(jìn)而提高計(jì)算效率奠定了基礎(chǔ).

圖7 三維CSD-FDM，TSD-FDM和MG-FDM(N=1,2,3)的穩(wěn)定性曲線(xiàn)Fig.7 Stability curves of 3D CSD-FDM, TSD-FDM and MG-FDM (N=1,2,3)

4 計(jì)算效率分析

計(jì)算效率是衡量有限差分格式優(yōu)劣的另一項(xiàng)重要指標(biāo).圖6中的相速度頻散曲線(xiàn)表明，三維CSD-FDM和TSD-FDM均可采用較小的時(shí)間采樣間隔來(lái)減小數(shù)值頻散，以獲得較高的模擬精度，但計(jì)算量會(huì)顯著增加.MG-FDM的相速度頻散曲線(xiàn)收斂性好，模擬精度高，穩(wěn)定性強(qiáng)，可采用更大的時(shí)間采樣間隔以提高計(jì)算效率，同時(shí)保持較高的模擬精度.

下面以三維CSD-FDM(M=10)、TSD-FDM(M=10)和MG-FDM(M=8;N=1)為例分析三種差分格式具有基本相同的模擬精度條件下的計(jì)算效率.這三種差分格式均包含61個(gè)網(wǎng)格點(diǎn)，單次時(shí)間迭代的浮點(diǎn)運(yùn)算量基本相同，模擬時(shí)間長(zhǎng)度相等時(shí)，計(jì)算效率之比約等于時(shí)間采樣間隔之比.

分析圖6中的相速度頻散曲線(xiàn)，可以認(rèn)為CSD-FDM(M=10)取r=0.1，TSD-FDM(M=10)取r=0.2和MG-FDM(M=8;N=1)取r=0.4基本能保持同樣高的模擬精度.根據(jù)r=vτ/h,在速度v和空間采樣間隔h保持不變的條件下，時(shí)間采樣間隔τ之比等價(jià)于r之比.此分析表明，模擬精度基本相當(dāng)?shù)那疤嵯拢琓SD-FDM可采用2倍于CSD-FDM的時(shí)間采樣間隔，MG-FDM可采用4倍于CSD-FDM的時(shí)間采樣間隔.因此，TSD-FDM的計(jì)算效率是CSD-FDM的2倍；MG-FDM的計(jì)算效率是CSD-FDM的4倍，是TSD-FDM的2倍.表2列出了模擬精度基本相當(dāng)時(shí)，三種差分格式采用的時(shí)間采樣間隔和理論加速比.

表2 三維CSD-FDM、TSD-FDM和MG-FDM三種差分格式的計(jì)算效率分析Table 2 Analysis of computational efficiency for 3D CSD-FDM, TSD-FDM and MG-FDM

5 數(shù)值模擬實(shí)例

5.1 層狀介質(zhì)模型

圖8a給出了一個(gè)8 km×3 km×6 km的三層模型：v1=2000 m·s-1,h1=3.0 km;v2=2450 m·s-1,h2=1.2 km;v3=2680 m·s-1,h3=1.8 km.空間采樣間隔Δx=Δy=Δz=h=10 m，模型網(wǎng)格數(shù)為nx×ny×nz=801×301×601，主頻25 Hz的Ricker子波作為震源，位于網(wǎng)格點(diǎn)(51,21,3).三維CSD-FDM(M=10)、TSD-FDM(M=10)和MG-FDM(M=8;N=1)分別采用不同的時(shí)間采樣間隔進(jìn)行數(shù)值模擬.圖8b給出了MG-FDM(M=8;N=1)采用時(shí)間采樣間隔τ=1.5 ms進(jìn)行數(shù)值模擬得到的炮集，為了便于對(duì)比分析，圖9給出了三維CSD-FDM(M=10)采用時(shí)間采樣間隔τ=0.5 ms(圖9a)和τ=1.0 ms(圖9b)，TSD-FDM(M=10)采用τ=1.0 ms(圖9c)和τ=1.5 ms(圖9d)，MG-FDM(M=8;N=1)采用τ=1.0 ms(圖9e)和τ=1.5 ms(圖9f)進(jìn)行數(shù)值模擬得到的局部炮集，局部區(qū)域范圍由坐標(biāo)標(biāo)出.

圖8 三維層狀速度模型及數(shù)值模擬炮集(a) 層狀速度模型； (b) MG-FDM(M=8;N=1)數(shù)值模擬炮集，時(shí)間采樣間隔τ=1.5 ms.Fig.8 3D layered velocity model and numerical modeling seismogram(a) 3D layered velocity model； (b) Numerical modeling seismogram using MG-FDM(M=8;N=1) with τ=1.5 ms.

圖9 三維層狀介質(zhì)模型數(shù)值模擬的局部炮集(a)(b) CSD-FDM(M=10), 時(shí)間采樣間隔τ=0.5 ms和τ=1.0 ms; (c)(d) TSD-FDM(M=10), 時(shí)間采樣間隔τ=1.0 ms和τ=1.5 ms; (e)(f) MG-FDM(M=8;N=1), 時(shí)間采樣間隔τ=1.0 ms和τ=1.5 ms.Fig.9 Local parts of the numerical modeling seismograms on 3D layered model(a) and (b) CSD-FDM(M=10) with τ=0.5 ms and τ=1.0 ms; (c) and (d) TSD-FDM(M=10) with τ=1.0 ms and τ=1.5 ms; (e) and (f) MG-FDM(M=8;N=1) with τ=1.0 ms and τ=1.5 ms.

圖10給出了三維CSD-FDM(M=10)采用τ=0.5 ms,0.75 ms,1.0 ms，TSD-FDM(M=10)采用τ=0.5 ms,1.0 ms,1.25 ms,1.5 ms和MG-FDM(M=8;N=1)采用τ=1.0 ms,1.25 ms,1.5 ms進(jìn)行數(shù)值模擬得到的單道波形圖，檢波器位于網(wǎng)格點(diǎn)(650,3,3).

圖10 三維層狀介質(zhì)模型數(shù)值模擬單道波形對(duì)比圖檢波器位于(650,3,3)，藍(lán)色為參考波形,CSD-FDM(M=10)采用極小時(shí)間采樣間隔τ=0.1 ms得到；紅色為不同差分格式模擬波形. ①②③分別為CSD-FDM(M=10), 時(shí)間采樣間隔τ=0.5 ms,0.75 ms,1.0 ms；④⑤⑥⑦分別為T(mén)SD-FDM(M=10)，時(shí)間采樣間隔τ=0.5 ms,1.0 ms,1.25 ms,1.5 ms; ⑧⑨⑩分別為M2M+N-FD(M=8;N=1), 時(shí)間采樣間隔τ=1.0 ms,1.25 ms,1.5 ms.Fig.10 Comparison of numerical modeling single trace waveforms on 3D layered model The geophone is located at (650,3,3), and the blue line is the reference waveform which is obtained by CSD-FDM (M=10) with a very small time sampling interval τ=0.1 ms; The red lines are the modeling waveform with different FD schemes. ①②③CSD-FDM(M=10) with τ=0.5 ms,0.75 ms,1.0 ms；④⑤⑥⑦TSD-FDM(M=10) with τ=0.5 ms,1.0 ms,1.25 ms,1.5 ms； ⑧⑨⑩ M2M+N-FD(M=8; N=1) with τ=1.0 ms,1.25 ms,1.5 ms.

從圖9和圖10可以看出：CSD-FDM采用時(shí)間采樣間隔τ=0.5 ms，未出現(xiàn)明顯數(shù)值頻散，采用τ=0.75 ms,1.0 ms，出現(xiàn)明顯的時(shí)間頻散；TSD-FDM采用τ=0.5 ms,1.0 ms，未出現(xiàn)明顯的數(shù)值頻散，采用τ=1.25 ms,1.5 ms，出現(xiàn)較明顯的數(shù)值頻散；MG-FDM采用τ=1.0 ms和τ=1.5 ms，均未出現(xiàn)明顯的數(shù)值頻散.表3給出了層狀介質(zhì)模型單炮模擬時(shí)，三種差分格式采用不同時(shí)間采樣間隔的計(jì)算機(jī)耗時(shí)和加速比，數(shù)值模擬均采用Inter Xeon CPU E5-2670處理器.

層狀模型數(shù)值模擬結(jié)果表明：確保不出現(xiàn)明顯數(shù)值頻散，實(shí)現(xiàn)高精度數(shù)值模擬，MG-FDM(M=8;N=1)可達(dá)到CSD-FDM(M=10)的計(jì)算效率的3.31倍，可達(dá)到TSD-FDM(M=10)的計(jì)算效率的1.5倍.

表3 三維層狀介質(zhì)模型CSD-FDM、TSD-FDM和MG-FDM數(shù)值模擬計(jì)算效率對(duì)比Table 3 Comparison of computational efficiency for numerical modeling on 3D layered model with CSD-FDM, TSD-FDM and MG-FDM

5.2 鹽丘模型

圖11b給出了MG-FDM(M=8;N=1)采用時(shí)間采樣間隔τ=1.5 ms進(jìn)行數(shù)值模擬得到的炮集，同樣為了分析對(duì)比方便，圖12僅給出三種差分格式數(shù)值模擬的兩個(gè)局部炮集.圖12a、b給出了CSD-FDM采用時(shí)間采樣間隔τ=1.0 ms，圖12c、d給出TSD-FDM采用τ=1.25 ms，圖12e、f給出MG-FDM采用τ=1.5 ms進(jìn)行數(shù)值模擬得到的兩個(gè)局部炮集.圖13給出了CSD-FDM采用τ=0.5 ms,1.0 ms，TSD-FDM(M=10)采用τ=1.0 ms,1.25 ms和MG-FDM采用τ=1.25 ms,1.5 ms的兩個(gè)單道波形記錄，兩個(gè)檢波器分別位于網(wǎng)格點(diǎn)(3,650,3)和(560,650,3).

圖11 修改的三維SEG/EAGE鹽丘速度模型及數(shù)值模擬炮集(a) 三維鹽丘速度模型； (b) MG-FDM(M=8;N=1)數(shù)值模擬炮集，時(shí)間采樣間隔τ=1.5 ms，頂部顯示了3.0 s時(shí)刻的切片.Fig.11 Revised 3D SEG/EAGE salt model and numerical modeling seismogram(a) Revised 3D SEG/EAGE salt model；(b) Numerical modeling seismogram using MG-FDM(M=8;N=1) with τ=1.5 ms. The top shows the slice at 3.0 s.

圖12 修改的三維SEG/EAGE鹽丘模型數(shù)值模擬的兩個(gè)局部炮集(a)(b) CSD-FDM(M=10), 時(shí)間采樣間隔τ=1.0 ms; (c)(d) TSD-FDM(M=10), 時(shí)間采樣間隔τ=1.25 ms; (e)(f) MG-FDM(M=8;N=1), 時(shí)間采樣間隔τ=1.5 ms.Fig.12 Two local parts of the numerical modeling seismogram on the revised 3D SEG/EAGE salt model(a) and (b) CSD-FDM(M=10) with τ=1.0 ms； (c) and (d) TSD-FDM(M=10) with τ=1.25 ms； (e) and (f) MG-FDM(M=8;N=1) with τ=1.5 ms.

從圖12和圖13可以看出：CSD-FDM采用τ=1.0 ms，表現(xiàn)出明顯的時(shí)間數(shù)值頻散；TSD-FDM采用τ=1.25 ms，表現(xiàn)出較明顯的數(shù)值頻散；MG-FDM采用τ=1.5 ms，未出現(xiàn)明顯的數(shù)值頻散.表4給出了鹽丘模型單炮模擬時(shí)，三種差分格式采用不同時(shí)間采樣間隔的計(jì)算機(jī)耗時(shí)和加速比.

圖13 修改的三維SEG/EAGE鹽丘模型數(shù)值模擬單道波形對(duì)比圖(a)(b)檢波器分別位于(3,650,3)和(560,650,3)，藍(lán)色為參考波形,CSD-FDM(M=10)采用極小時(shí)間采樣間隔τ=0.1 ms得到；紅色為不同差分格式模擬波形. ①②分別為CSD-FDM(M=10), 時(shí)間采樣間隔τ=0.5 ms,1.0 ms；③④分別為T(mén)SD-FDM(M=10)，時(shí)間采樣間隔τ=1.0 ms,1.25 ms;⑤⑥分別為M2M+N-FD(M=8;N=1), 時(shí)間采樣間隔τ=1.25 ms,1.5 ms.Fig.13 Numerical modeling single trace waveform on the revised 3D SEG/EAGE salt model(a) and (b) Geophone at (3,650,3) and (560,650,3). ① and ② CSD-FDM(M=10) with τ=0.5 ms,1.0 ms；③ and ④ TSD-FDM(M=10) with τ=1.0 ms,1.25 ms； ⑤and ⑥ M2M+N-FD(M=8;N=1) with τ=1.25 ms,1.5 ms.

表4 修改的三維SEG/EAGE鹽丘模型CSD-FDM、TSD-FDM和MG-FDM數(shù)值模擬計(jì)算效率對(duì)比Table 4 Comparison of computational efficiency for numerical modeling on revised 3D SEG/EAGE Salt model with CSD-FDM, TSD-FDM and MG-FDM

鹽丘模型數(shù)值模擬結(jié)果表明，三維MG-FDM(M=8;N=1)能采用更大的時(shí)間采樣間隔以提高計(jì)算效率，同時(shí)能更有效地減小數(shù)值頻散，獲得更高的模擬精度.

6 結(jié)論與討論

本文提出了一種面向三維標(biāo)量波動(dòng)方程數(shù)值模擬的混合網(wǎng)格有限差分方法(MG-FDM)，將三維Laplace微分算子近似表示為坐標(biāo)軸網(wǎng)格點(diǎn)構(gòu)建的M個(gè)Laplace差分算子和非坐標(biāo)軸網(wǎng)格點(diǎn)構(gòu)建的N個(gè)Laplace差分算子的加權(quán)平均，充分利用了與差分中心點(diǎn)距離更近的非坐標(biāo)軸網(wǎng)點(diǎn)，有效減小數(shù)值頻散和數(shù)值各向異性，獲得更高的模擬精度.通過(guò)對(duì)比分析，得出以下結(jié)論：

(1)三維CSD-FDM的空間差分算子(Laplace 差分算子)雖然能達(dá)到2M階差分精度，但其差分離散波動(dòng)方程僅具有2階差分精度；TSD-FDM的差分離散波動(dòng)方程沿特定的48個(gè)傳播方向可達(dá)到2M階差分精度，但沿其他傳播方向僅具有2階差分精度；MG-FDM的差分離散波動(dòng)方程沿任意傳播方向可達(dá)到4階、6階、甚至任意偶數(shù)階差分精度.

(2)相比三維CSD-FDM和TSD-FDM，計(jì)算效率基本相同時(shí)，MG-FDM的數(shù)值頻散壓制效果更好，模擬精度更高；模擬精度基本相同時(shí)，MG-FDM可以采用更大的時(shí)間采樣間隔以提高計(jì)算效率；穩(wěn)定性分析表明，MG-FDM穩(wěn)定性更強(qiáng)，為采用更大的時(shí)間采樣間隔以提高計(jì)算效率奠定了基礎(chǔ).

(3)數(shù)值模擬實(shí)例進(jìn)一步證實(shí)了三維MG-FDM在模擬精度和計(jì)算效率方面的優(yōu)越性，也驗(yàn)證了其普遍適用性.

但是，本文提出的三維MG-FDM也具有一定的局限性，推導(dǎo)過(guò)程隱含了三個(gè)方向的空間采樣間隔相等，即MG-FDM要求采用立方體網(wǎng)格.適用于長(zhǎng)方體網(wǎng)格的更具一般性的MG-FDM，有待進(jìn)一步研究.

附錄A 三維坐標(biāo)系中非坐標(biāo)軸網(wǎng)格點(diǎn)構(gòu)建Laplace差分算子

(A1)

(A2)

式(A2)給出了圖2a中12個(gè)坐標(biāo)平面內(nèi)的非坐標(biāo)軸網(wǎng)格點(diǎn)構(gòu)建三維Laplace差分算子的表達(dá)式.

(A3)

(A4)

(A5)

(A6)

式(A6)給出了圖2c中24個(gè)坐標(biāo)平面內(nèi)的非坐標(biāo)軸網(wǎng)格點(diǎn)構(gòu)建三維Laplace差分算子的表達(dá)式.

≈P(x,y,z,t)+Px(x,y,z,t)h+Py(x,y,z,t)h+Pz(x,y,z,t)h

+Pxy(x,y,z,t)h2+Pyz(x,y,z,t)h2+Pzx(x,y,z,t)h2.

(A7)

≈8P(x,y,z,t)+4Pxx(x,y,z,t)h2+4Pyy(x,y,z,t)h2+4Pzz(x,y,z,t)h2，

(A8)

(A9)

式(A9)給出了圖2b中8個(gè)坐標(biāo)平面外的非坐標(biāo)軸網(wǎng)格點(diǎn)構(gòu)建三維Laplace差分算子的表達(dá)式.

附錄B 三維MG-FDM(N=2,3)的差分離散波動(dòng)方程和差分系數(shù)通解

圖3b中三維MG-FDM(N=2)的差分格式由圖1和圖2a、b中的差分格式組合構(gòu)建，三維MG-FDM(N=2)的差分離散波動(dòng)方程可表示為：

(B1)

三維MG-FDM(N=2)的差分系數(shù)通解為：

(B2)

圖3c中三維MG-FDM(N=3)的差分格式由圖1和圖2a、b、c中的差分格式組合構(gòu)建，三維MG-FDM(N=3)的差分離散波動(dòng)方程可表示為：

(B3)

(B4)