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

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

0 引言

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

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

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

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

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

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

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

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

(1)

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

(2)

三維常規(guī)空間域高階差分(CSD-FDM)和三維時空域高階差分(TSD-FDM)，均采用圖1所示的差分格式，波動方程(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)鍵問題.

圖2a—c給出了三維笛卡爾坐標(biāo)系中與差分中心點(diǎn)距離相等的三組非坐標(biāo)軸網(wǎng)格點(diǎn)，這三組網(wǎng)格點(diǎn)與差分中心點(diǎn)的距離依次增大.由于無法將任意一組非坐標(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)分成兩類：坐標(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個Laplace差分算子和非坐標(biāo)軸網(wǎng)格點(diǎn)構(gòu)建的1個Laplace差分算子的加權(quán)平均.式(7)和(2)代入(1)：

(8)

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

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

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

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

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

(9)

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

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

+2cos(mkzh)-6].

(10)

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

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

(11)

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

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

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

(12)

(13)

其中r=υτ/h為Courant條件數(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)計算差分系數(shù)時，取φ=π/2,θ=π/8，分析出三維TSD-FDM的差分離散波動方程沿48個傳播方向具有2M階差分精度.

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

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

(15)

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

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

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

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

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

(16)

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

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

(17)

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

(18)

(19)

2.4 差分精度分析和對比

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

(20)

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

(21)

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

(22)

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

根據(jù)三維MG-FDM(N=1)的差分離散波動方程的頻散關(guān)系式(16)，定義MG-FDM(N=1)的差分離散波動方程的誤差函數(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能夠有效提高差分離散波動方程的差分精度.增大N的取值，MG-FDM的差分離散波動方程理論上可達(dá)到任意偶數(shù)階差分精度.

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

(25)

(26)

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

歸一化相速度δ(φ,θ)的值越接近1，相速度數(shù)值頻散越小；δ(φ,θ)>1，相速度偏大，稱為時間數(shù)值頻散；δ(φ,θ)<1，相速度偏小，稱為空間數(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ù)，分析相速度頻散特性時給出了8個傳播角度(φ=π/2;θ=0,π/16,2π/16,3π/16,4π/16)和(φ=π/12,2π/12,3π/12;θ=2π/16)的相速度頻散曲線.

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

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

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

圖4 三維CSD-FDM、TSD-FDM和MG-FDM(N=1)的相速度頻散曲線(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取值較小時(M=2左右)，MG-FDM在壓制相速度數(shù)值頻散方面無明顯優(yōu)勢；M取值較大時(M≥6)，MG-FDM的相速度數(shù)值頻散和數(shù)值各向異性均明顯減小，因而具有更高的模擬精度.

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

圖5 三維MG-FDM的相速度頻散曲線(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)的相速度頻散曲線更收斂，數(shù)值頻散幅值和數(shù)值各向異性均明顯減小.因此，三維標(biāo)量波動方程數(shù)值模擬對精度要求較高時，建議采用MG-FDM(N=1)，且M取值為8左右；對模擬精度要求極其苛刻時，可采用MG-FDM(N=3)，且M取值為15左右；N取值大于4，同時M取值更大的三維MG-FDM因為數(shù)值計算量巨大，計算效率極低而很少采用.三維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)的差分離散波動方程均具有4階差分精度的結(jié)論一致(見表1).

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

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

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

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

圖6 三維CSD-FDM(M=10)，TSD-FDM(M=10)和MG-FDM(M=8;N=1)的相速度頻散曲線(a)(b)(c) 分別為r=0.1,0.2,0.4時CSD-FDM(M=10)的相速度頻散曲線； (d)(e)(f) 分別為r=0.1,0.2,0.4時TSD-FDM(M=10)的相速度頻散曲線； (g)(h)(i) 分別為r=0.1,0.2,0.4時MG-FDM(M=8;N=1)的相速度頻散曲線.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相同時，r取值越大等價于時間采樣間隔τ越大)，同時數(shù)值頻散大小基本相當(dāng).因此，MG-FDM(M=8;N=1)能采用較大的時間采樣間隔以提高計算效率，同時保持較高的模擬精度.

3.2 穩(wěn)定性分析

有限差分法通過迭代求解差分離散波動方程模擬地震波的傳播過程，必須確保迭代過程穩(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)定性曲線，描述了穩(wěn)定性條件約束下的最大r取值隨M的變化關(guān)系.穩(wěn)定性條件約束前提下，r取值越小，穩(wěn)定性越弱；相反，r取值越大，穩(wěn)定性越強(qiáng).分析圖中的穩(wěn)定性曲線可以看出：

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

(2)N的取值增大，MG-FDM的穩(wěn)定性增強(qiáng)，但差分離散波動方程的差分精度相同時，穩(wěn)定性基本相同，MG-FDM(N=1,2)的差分離散波動方程均為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ù)值模擬采用更大的時間采樣進(jìn)而提高計算效率奠定了基礎(chǔ).

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

4 計算效率分析

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

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

分析圖6中的相速度頻散曲線，可以認(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保持不變的條件下，時間采樣間隔τ之比等價于r之比.此分析表明，模擬精度基本相當(dāng)?shù)那疤嵯?，TSD-FDM可采用2倍于CSD-FDM的時間采樣間隔，MG-FDM可采用4倍于CSD-FDM的時間采樣間隔.因此，TSD-FDM的計算效率是CSD-FDM的2倍；MG-FDM的計算效率是CSD-FDM的4倍，是TSD-FDM的2倍.表2列出了模擬精度基本相當(dāng)時，三種差分格式采用的時間采樣間隔和理論加速比.

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

5 數(shù)值模擬實例

5.1 層狀介質(zhì)模型

圖8a給出了一個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)分別采用不同的時間采樣間隔進(jìn)行數(shù)值模擬.圖8b給出了MG-FDM(M=8;N=1)采用時間采樣間隔τ=1.5 ms進(jìn)行數(shù)值模擬得到的炮集，為了便于對比分析，圖9給出了三維CSD-FDM(M=10)采用時間采樣間隔τ=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ù)值模擬炮集，時間采樣間隔τ=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), 時間采樣間隔τ=0.5 ms和τ=1.0 ms; (c)(d) TSD-FDM(M=10), 時間采樣間隔τ=1.0 ms和τ=1.5 ms; (e)(f) MG-FDM(M=8;N=1), 時間采樣間隔τ=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ù)值模擬單道波形對比圖檢波器位于(650,3,3)，藍(lán)色為參考波形,CSD-FDM(M=10)采用極小時間采樣間隔τ=0.1 ms得到；紅色為不同差分格式模擬波形. ①②③分別為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; ⑧⑨⑩分別為M2M+N-FD(M=8;N=1), 時間采樣間隔τ=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采用時間采樣間隔τ=0.5 ms，未出現(xiàn)明顯數(shù)值頻散，采用τ=0.75 ms,1.0 ms，出現(xiàn)明顯的時間頻散；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ì)模型單炮模擬時，三種差分格式采用不同時間采樣間隔的計算機(jī)耗時和加速比，數(shù)值模擬均采用Inter Xeon CPU E5-2670處理器.

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

表3 三維層狀介質(zhì)模型CSD-FDM、TSD-FDM和MG-FDM數(shù)值模擬計算效率對比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)采用時間采樣間隔τ=1.5 ms進(jìn)行數(shù)值模擬得到的炮集，同樣為了分析對比方便，圖12僅給出三種差分格式數(shù)值模擬的兩個局部炮集.圖12a、b給出了CSD-FDM采用時間采樣間隔τ=1.0 ms，圖12c、d給出TSD-FDM采用τ=1.25 ms，圖12e、f給出MG-FDM采用τ=1.5 ms進(jìn)行數(shù)值模擬得到的兩個局部炮集.圖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的兩個單道波形記錄，兩個檢波器分別位于網(wǎng)格點(diǎn)(3,650,3)和(560,650,3).

圖11 修改的三維SEG/EAGE鹽丘速度模型及數(shù)值模擬炮集(a) 三維鹽丘速度模型； (b) MG-FDM(M=8;N=1)數(shù)值模擬炮集，時間采樣間隔τ=1.5 ms，頂部顯示了3.0 s時刻的切片.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ù)值模擬的兩個局部炮集(a)(b) CSD-FDM(M=10), 時間采樣間隔τ=1.0 ms; (c)(d) TSD-FDM(M=10), 時間采樣間隔τ=1.25 ms; (e)(f) MG-FDM(M=8;N=1), 時間采樣間隔τ=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ù)值頻散；TSD-FDM采用τ=1.25 ms，表現(xiàn)出較明顯的數(shù)值頻散；MG-FDM采用τ=1.5 ms，未出現(xiàn)明顯的數(shù)值頻散.表4給出了鹽丘模型單炮模擬時，三種差分格式采用不同時間采樣間隔的計算機(jī)耗時和加速比.

圖13 修改的三維SEG/EAGE鹽丘模型數(shù)值模擬單道波形對比圖(a)(b)檢波器分別位于(3,650,3)和(560,650,3)，藍(lán)色為參考波形,CSD-FDM(M=10)采用極小時間采樣間隔τ=0.1 ms得到；紅色為不同差分格式模擬波形. ①②分別為CSD-FDM(M=10), 時間采樣間隔τ=0.5 ms,1.0 ms；③④分別為TSD-FDM(M=10)，時間采樣間隔τ=1.0 ms,1.25 ms;⑤⑥分別為M2M+N-FD(M=8;N=1), 時間采樣間隔τ=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ù)值模擬計算效率對比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ù)值頻散，獲得更高的模擬精度.

6 結(jié)論與討論

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

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

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

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

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

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

(A1)

(A2)

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

(A3)

(A4)

(A5)

(A6)

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

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

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

(B1)

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

(B2)

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

(B3)

(B4)