胡自多， 賀振華， 劉威， 王宇超， 韓令賀， 王述江， 楊哲

1 成都理工大學(xué)油氣藏地質(zhì)及開發(fā)工程國家重點(diǎn)實(shí)驗(yàn)室， 成都　610059 2 中國石油勘探開發(fā)研究院西北分院， 蘭州　730020

?

旋轉(zhuǎn)網(wǎng)格和常規(guī)網(wǎng)格混合的時(shí)空域聲波有限差分正演

胡自多1, 2， 賀振華1， 劉威2， 王宇超2， 韓令賀2， 王述江2， 楊哲2

1 成都理工大學(xué)油氣藏地質(zhì)及開發(fā)工程國家重點(diǎn)實(shí)驗(yàn)室， 成都610059 2 中國石油勘探開發(fā)研究院西北分院， 蘭州730020

壓制數(shù)值頻散，提高正演模擬精度，一直是有限差分正演模擬研究的重要內(nèi)容.基于時(shí)空域頻散關(guān)系的有限差分法，比基于空間域頻散關(guān)系的傳統(tǒng)有限差分法，模擬精度更高.時(shí)空域聲波方程數(shù)值模擬，普遍采用常規(guī)十字交叉型高階有限差分格式.而在頻率-空間域，普遍采用旋轉(zhuǎn)網(wǎng)格和常規(guī)網(wǎng)格混合的有限差分格式，有效提高了模擬精度和計(jì)算效率.本文將頻率-空間域混合網(wǎng)格有限差分的思想引入到時(shí)空域，提出了時(shí)空域混合網(wǎng)格2M+N型聲波方程有限差分方法.首先推導(dǎo)出基于時(shí)空域頻散關(guān)系的混合網(wǎng)格差分系數(shù)計(jì)算方法，然后進(jìn)行頻散分析、穩(wěn)定性分析，并和傳統(tǒng)高階、時(shí)空域高階有限差分法對比，結(jié)果表明：計(jì)算量相同時(shí)，新方法能有效壓制數(shù)值頻散，顯著提高模擬精度；新方法相比傳統(tǒng)2M階有限差分法，穩(wěn)定性增強(qiáng)，與時(shí)空域2M階有限差分法穩(wěn)定性基本相當(dāng).最后利用新方法進(jìn)行均勻介質(zhì)、層狀介質(zhì)、鹽丘模型的數(shù)值模擬和鹽丘模型的逆時(shí)偏移，模擬效果和成像質(zhì)量進(jìn)一步證實(shí)了該方法的有效性和普遍適用性.

混合網(wǎng)格； 時(shí)空域； 有限差分； 頻散分析； 正演模擬

In the frequency-space domain, the mixed-grid finite-difference method is often used, which can effectively improve the forward modeling accuracy. In the time-space domain, the traditional 2Mth-order finite-difference method is commonly used, which essentially has only 2nd-order accuracy. The time-space domain 2Mth-order finite-difference method, in which the difference coefficients are determined by satisfying time-space dispersion relation, has relatively high modeling accuracy, but its dispersion curves are still divergent. Although the rhombus-grid finite-difference method has indeed high modeling accuracy, it requires high computational cost.

In this article, by introducing the mixed-grid strategy from the frequency-space domain to the time-space domain, we propose a new kind of mixed-grid 2M+Nstyle finite-difference method, and derive the approach for calculating the finite-difference coefficients by satisfying time-space domain dispersion relation. In addition, we conduct dispersion analysis, stability analysis, and numerical simulation. The results demonstrate that (1) with almost the same computational cost, the traditional 2Mth-order finite-difference method has seriously numerical dispersion mainly in the time dispersion, and has the lowest modeling accuracy.The time-space domain 2Mth-order finite-difference method has some time dispersion and space dispersion, and has relatively high modeling accuracy. The mixed-grid 2M+Nstyle finite-difference method has no obvious numerical dispersion, and so has the highest modeling accuracy. (2) WhenMis not too big, a biggerNvalue can hardly decrease the numerical dispersion and increase the forward modeling accuracy, but will increase the computational cost. It suggests that we should take use of the mixed-grid 2M+N(N=1) style finite-difference method for general conditions. The rhombus-grid is a special shape of the mixed-grid 2M+Nstyle finite-difference method, in which theNvalue is usually too big, so the rhombus-grid finite-difference method is not the optimal choice. (3) The mixed-grid 2M+N(N=1) style finite-difference method has stronger stability than the traditional 2Mth-order finite-difference method, and has almost the same stability as the time-space domain 2Mth-order finite-difference method. A biggerNvalue will slightly improve the stability, but also increases the computational amount. (4) Numerical modeling on homogeneous model, layer model and salt model further demonstrates the mixed-grid 2M+Nstyle finite-difference method can effectively suppress the numerical dispersion, and improve the modeling accuracy.

In the end, we effectively eliminate most of the energy reflected from the artificial boundary by adopting an NPML (Nearly Perfectly Matched Layer) absorbing boundary, and carry out forward modeling and RTM on the salt model. There is no obvious numerical dispersion and boundary reflection on the shot gathers of the forward modeling, and RTM results have really good quality. All of these prove the validity and applicability of the mixed-grid 2M+Nstyle finite-difference method suggested in this study.

1　引言

有限差分正演模擬方法具有計(jì)算效率高、內(nèi)存占用小、簡單易實(shí)現(xiàn)的優(yōu)點(diǎn)，因而廣泛應(yīng)用于地震波正演模擬(Alterman and Karal, 1968; Alford et al., 1974; Kelly et al., 1976)、逆時(shí)偏移(Sun and McMechan, 2001; Sun et al., 2006; Jones, 2008; Etgen et al., 2009)和全波形反演(Gerhard et al., 1998; Virieux and Operto, 2009).然而，有限差分正演模擬中普遍存在的數(shù)值頻散降低了模擬精度.數(shù)值頻散是對空間和時(shí)間偏導(dǎo)數(shù)進(jìn)行網(wǎng)格差分近似造成的，時(shí)間數(shù)值頻散通常使相速度變大而出現(xiàn)“波至超前”，空間數(shù)值頻散使相速度變小而出現(xiàn)“波至拖尾”，并且頻率越大，數(shù)值頻散誤差越大(Dablain, 1986)，嚴(yán)重影響了地震波高頻成分的模擬精度.因此，壓制數(shù)值頻散一直是有限差分正演模擬的重要研究內(nèi)容.

為了壓制數(shù)值頻散，高階有限差分法、隱式有限差分法、低秩有限差分法相繼被提出.隱式有限差分法需要求解大型線性方程組，內(nèi)存占用大，計(jì)算效率低(Tian and Dai, 2007; Sherer and Scott, 2005; Nihei and Ishii, 2003)；低秩有限差分法也需要大量的快速傅里葉變換，導(dǎo)致計(jì)算效率低(Fomel et al., 2013; Fang et al., 2014)；高階有限差分法被認(rèn)為是一種兼顧效率和精度的有效方法，廣泛應(yīng)用于地震波正演模擬(Dablain, 1986; Levander, 1988; 劉洋等, 1998)和逆時(shí)偏移波場延拓(劉紅偉等，2010a).

波動方程有限差分?jǐn)?shù)值求解在時(shí)間域和空間域同時(shí)進(jìn)行，而傳統(tǒng)高階有限差分法(本文稱為傳統(tǒng)2M階有限差分法)，僅基于空間域頻散關(guān)系求取差分系數(shù)，本質(zhì)上只有二階空間精度.Liu和Sen(2009)提出了基于時(shí)空域頻散關(guān)系的規(guī)則網(wǎng)格高階有限差分法(本文稱為時(shí)空域2M階有限差分法)，并證明該方法在二維正演模擬中沿8個(gè)方向達(dá)到2M階空間精度，三維中沿48個(gè)方向達(dá)到2M階空間精度.與傳統(tǒng)2M階有限差分方法相比，時(shí)空域2M階有限差分法正演模擬精度明顯提高.Liu和Sen(2011)進(jìn)一步把基于時(shí)空域頻散關(guān)系有限差分的思想推廣到高階交錯(cuò)網(wǎng)格，模擬精度也明顯提高.嚴(yán)紅勇和劉洋(2013)將時(shí)空域2M階有限差分法推廣到逆時(shí)偏移，通過自適應(yīng)減小空間差分階數(shù)而提高計(jì)算效率.傳統(tǒng)2M階和時(shí)空域2M階有限差分正演，目前普遍采用十字交叉型規(guī)則網(wǎng)格或交錯(cuò)網(wǎng)格差分格式.

Saenger等(2000)，Saenger和Bohlen(2004)針對彈性、黏彈性和各向異性介質(zhì)提出了旋轉(zhuǎn)交錯(cuò)網(wǎng)格有限差分法，避免對彈性參數(shù)求平均，提高了正演模擬精度.旋轉(zhuǎn)交錯(cuò)網(wǎng)格有限差分正演模擬普遍應(yīng)用于彈性、黏彈性和各向異性參數(shù)存在突變的介質(zhì)(嚴(yán)紅勇和劉洋, 2012).Liu和Sen(2013)改進(jìn)時(shí)空域2M階有限差分格式，由十字交叉網(wǎng)格推廣到菱形網(wǎng)格，提出了二維聲波方程菱形網(wǎng)格高階有限差分方法，沿任意方向能夠達(dá)到2M階空間精度，進(jìn)一步提高了正演模擬精度.但是，當(dāng)差分階數(shù)增大時(shí)計(jì)算量明顯增大，而且較難推廣到交錯(cuò)網(wǎng)格和三維正演.

頻率-空間域，普遍采用旋轉(zhuǎn)網(wǎng)格和常規(guī)網(wǎng)格混合的差分格式，并且有效地提高了正演模擬的精度和計(jì)算效率.Jo等(1996)提出9點(diǎn)混合網(wǎng)格有限差分格式，在增加少量計(jì)算量的情況下，有效提高了正演模擬精度.借鑒混合網(wǎng)格思想，Shin和Sohn(1998)提出25點(diǎn)混合網(wǎng)格有限差分格式，曹書紅和陳景波(2012)提出17點(diǎn)有限差分格式，正演模擬精度進(jìn)一步提高.這些方法僅適用于正方形網(wǎng)格，對于矩形網(wǎng)格，Chen(2012)對9點(diǎn)混合網(wǎng)格有限差分格式進(jìn)行改進(jìn)，提出了平均導(dǎo)數(shù)法9點(diǎn)混合網(wǎng)格有限差分格式，該方法實(shí)用性更強(qiáng).平均導(dǎo)數(shù)法混合網(wǎng)格有限差分格式被廣泛應(yīng)用，并不斷改進(jìn)(劉璐等, 2013; Tang et al., 2015; 唐祥德等, 2015; 張衡等, 2015).

吸收邊界是地震波數(shù)值模擬中為了消除人工邊界反射而引入的.一階應(yīng)力速度聲波方程正演模擬時(shí)，普遍采用分裂形式的PML吸收邊界(Collino and Tsogka, 2001).二階標(biāo)量聲波方程進(jìn)行正演模擬時(shí)，采用CPML(Convolutional Perfectly Matched Layer)吸收邊界(Komatitsch and Martin, 2007; Chen et al., 2014)和NPML(Nearly Perfectly Matched Layer)吸收邊界(Hu et al., 2007; Chen, 2011).NPML吸收邊界具有輔助變量少，易于實(shí)現(xiàn)的優(yōu)點(diǎn).本文提出的混合網(wǎng)格2M+N型有限差分法，基于二階標(biāo)量聲波方程，因此采用NPML吸收邊界.

本文將頻率-空間域混合網(wǎng)格有限差分的思想引入到時(shí)空域，利用常規(guī)直角坐標(biāo)系中M個(gè)拉普拉斯算子和旋轉(zhuǎn)直角坐標(biāo)系中N個(gè)拉普拉斯算子，構(gòu)建聲波方程混合網(wǎng)格2M+N型有限差分格式.借鑒Liu和Sen(2009, 2013)基于時(shí)空域頻散關(guān)系計(jì)算差分系數(shù)的思路，推導(dǎo)出混合網(wǎng)格2M+N型有限差分格式的差分系數(shù)計(jì)算方法.在此基礎(chǔ)上，進(jìn)行頻散分析、穩(wěn)定性分析，并與傳統(tǒng)2M階、時(shí)空域2M階有限差分方法進(jìn)行比較，結(jié)果表明：當(dāng)計(jì)算量基本相同時(shí)，混合網(wǎng)格2M+N型有限差分法能有效壓制數(shù)值頻散，模擬精度顯著提高；混合網(wǎng)格有限差分法，相比傳統(tǒng)2M階有限差分法，穩(wěn)定性增強(qiáng)，與時(shí)空域2M階有限差分法穩(wěn)定性基本相當(dāng).最后進(jìn)一步將混合網(wǎng)格有限差分法應(yīng)用于均勻介質(zhì)、層狀介質(zhì)和鹽丘模型的正演模擬及鹽丘模型的逆時(shí)偏移，采用NPML吸收邊界有效壓制人工邊界反射，正演模擬精度對比分析和逆時(shí)偏移成像質(zhì)量，進(jìn)一步證實(shí)了該方法的有效性和普遍適用性.

2　方法原理

2.1聲波方程混合網(wǎng)格2M+N型有限差分格式構(gòu)建

二維時(shí)間-空間域常密度標(biāo)量聲波方程可寫為：

(1)

對方程(1)進(jìn)行有限差分法正演模擬時(shí)，由于時(shí)間高階差分內(nèi)存占用大、計(jì)算量大，并且穩(wěn)定性降低(Chen, 2007)，一般采用二階時(shí)間精度、2M階空間精度.時(shí)間二階偏導(dǎo)數(shù)進(jìn)行二階有限差分離散近似，可以得到：

(2)

傳統(tǒng)2M階有限差分法(劉洋等, 1998)和時(shí)空域2M階有限差分法(Liu and Sen, 2009)均采用圖1所示的十字交叉型差分格式，對方程(1)中的拉普拉斯算子進(jìn)行2M階有限差分離散近似，可以得到：

(3)

圖1　傳統(tǒng)2M階和時(shí)空域2M階有限差分格式Fig.1　Traditional 2Mth-order and time-space domain 2Mth-order finite-difference scheme

將式(2)和(3)代入式(1)，可以得到：

(4)

其中，r=v τ/h，r為Courant穩(wěn)定性條件數(shù).

方程(4)是傳統(tǒng)2M階和時(shí)空域2M階聲波方程有限差分法的離散化方程，但兩種方法求取差分系數(shù)的方法不同.傳統(tǒng)2M階有限差分法僅基于空間域頻散關(guān)系計(jì)算差分系數(shù)，而時(shí)空域2M階有限差分法則基于時(shí)空域頻散關(guān)系計(jì)算差分系數(shù).

我們將頻率-空間域混合網(wǎng)格有限差分的思想引入到時(shí)空域，提出了時(shí)空域混合網(wǎng)格2M+N型有限差分法，M表示常規(guī)直角坐標(biāo)系中的M個(gè)拉普拉斯算子，N表示旋轉(zhuǎn)直角坐標(biāo)系中的N個(gè)拉普拉斯算子.圖2a—2f給出了6種混合網(wǎng)格2M+N型有限差分格式的示意圖.利用相同的方法，N取更大的數(shù)值，會得到更為豐富的混合網(wǎng)格2M+N型有限差分格式.

此外，與Liu和Sen(2013)給出的菱形網(wǎng)格有限差分法進(jìn)行對比可以發(fā)現(xiàn)，菱形網(wǎng)格是混合網(wǎng)格中M和N取特定值的特殊格式，而混合網(wǎng)格是菱形網(wǎng)格的擴(kuò)展，更具一般性和靈活性.例如，混合網(wǎng)格2M+N(M=4,N=6)型與四階菱形有限差分格式完全等價(jià).

下面以混合網(wǎng)格2M+1型有限差分格式為例，給出聲波方程(1)的有限差分離散形式.按照圖2a所示的差分格式，對拉普拉斯算子進(jìn)行有限差分離散近似，可以得到：

(5)式(5)表示，在混合網(wǎng)格2M+1型有限差分格式中，拉普拉斯算子可以表示為常規(guī)直角坐標(biāo)系中M個(gè)拉普拉斯算子和旋轉(zhuǎn)直角坐標(biāo)系中1個(gè)拉普拉斯算子的加權(quán)平均，其中a1,a2,…,am;a1,1為對應(yīng)的差分系數(shù).

將式(5)和(2)代入聲波方程(1)，可以得到：

(6)

方程(6)給出了混合網(wǎng)格2M+1型有限差分格式進(jìn)行聲波方程正演模擬的離散形式.利用相同的方法，可以推導(dǎo)出混合網(wǎng)格2M+N型有限差分格式的離散形式.

在模型網(wǎng)格化參數(shù)給定的情況下，拉普拉斯算子長度(網(wǎng)格點(diǎn)數(shù))決定了有限差分正演模擬的計(jì)算量.拉普拉斯算子越長，計(jì)算量越大.表1給出了傳統(tǒng)2M階，時(shí)空域2M階，混合網(wǎng)格2M+N型(N=1,2,…,6)有限差分格式的拉普拉斯算子長度.本文用拉普拉斯算子長度來表述不同差分格式的計(jì)算效率.

2.2混合網(wǎng)格2M+N型有限差分格式差分系數(shù)計(jì)算

圖2　混合網(wǎng)格2M+N型有限差分格式 (a)—(f)分別為混合網(wǎng)格2M+N(N=1,2,…,6)型.Fig.2　Mixed-grid 2M+N style finite-difference schemes(a)—(f) represent mixed-grid 2M+N (N=1,2,…,6), respectively.

差分格式傳統(tǒng)2M階(時(shí)空域2M階)混合網(wǎng)格2M+1型混合網(wǎng)格2M+2型混合網(wǎng)格2M+3型混合網(wǎng)格2M+4型混合網(wǎng)格2M+5型混合網(wǎng)格2M+6型拉普拉斯算子長度4M+14M+54M+94M+134M+174M+214M+25

方程(6)中差分系數(shù)a1,a2,…,am;a1,1的求解是混合網(wǎng)格2M+N型有限差分法的關(guān)鍵環(huán)節(jié).我們以混合網(wǎng)格2M+1型有限差分格式為例，闡述混合網(wǎng)格2M+N型差分系數(shù)的計(jì)算方法.

方程(1)的平面波解的離散形式為：

(7)

(8)

其中，k為波數(shù)，kx、kz分別為水平波數(shù)和垂直波數(shù)，ω為圓頻率，θ表示平面波傳播方向與x軸正半軸的夾角，i2=-1.

將平面波解(7)代入式(6)，可以得到：

+a1,1[cos(kxh-kzh)+cos(kxh+kzh)-2]

(9)

對式(9)中的余弦函數(shù)進(jìn)行泰勒級數(shù)展開，ωτ=kvτ=khr，可以得到：

(10)

(11)

(12)

方程(11)和方程(12)共有M+1個(gè)方程，剛好可以求解出M+1個(gè)差分系數(shù)a1,a2,…,am;a1,1.

從式(11)和式(12)可以看出，混合網(wǎng)格差分系數(shù)不僅與M有關(guān)，還與N、r有關(guān).而時(shí)空域2M階有限差分格式的差分系數(shù)僅與M和r有關(guān)，傳統(tǒng)2M階有限差分格式的差分系數(shù)只與M有關(guān).附錄A給出了混合網(wǎng)格2M+2型有限差分格式的差分系數(shù)計(jì)算方法.利用同樣的方法可以推導(dǎo)出其他混合網(wǎng)格2M+N型有限差分格式的差分系數(shù).

表2給出了M=5,v=3000 m·s-1,τ=0.001 s,h=10 m時(shí)的8種差分格式的差分系數(shù).

表2　傳統(tǒng)2M階，時(shí)空域2M階和混合網(wǎng)格2M+N型有限差分格式的差分系數(shù)Table 2　Finite difference coefficients for traditional 2Mth-order, time-space domain 2Mth-order, and mixed-grid 2M+N style finite-difference schemes

3　頻散分析和穩(wěn)定性分析

3.1頻散分析

地震波在非耗散介質(zhì)中傳播時(shí)，不會產(chǎn)生頻散現(xiàn)象，即不同頻率成分的地震波速度相同.然而，利用有限差分法求解波動方程時(shí)，對波動方程進(jìn)行差分離散會導(dǎo)致不同頻率成分的地震波傳播速度不再相等，即出現(xiàn)數(shù)值頻散現(xiàn)象.數(shù)值頻散是利用有限差分法求解波動方程時(shí)固有的本質(zhì)特征，無法避免，只能減小.通常用數(shù)值頻散的大小來衡量有限差分正演模擬的精度.因此可以定義歸一化相速度δ來描述有限差分格式的數(shù)值頻散：

(13)

其中vph為相速度，v為地震波在介質(zhì)傳播中的真實(shí)速度.

結(jié)合相速度的定義vph=ω/k和方程(9)可以得出混合網(wǎng)格2M+1型有限差分格式的頻散關(guān)系式為：

(14)

(15)

其中G=λ/h，λ為波長，則G表示單位波長內(nèi)網(wǎng)格點(diǎn)的個(gè)數(shù).

δ的值與1越接近，表明數(shù)值頻散誤差越??；δ>1表示存在時(shí)間數(shù)值頻散，δ<1表示存在空間數(shù)值頻散.根據(jù)同樣的方法可以計(jì)算出傳統(tǒng)2M階、時(shí)空域2M階和其他混合網(wǎng)格2M+N型有限差分格式的頻散關(guān)系.

圖3給出了傳統(tǒng)2M(M=5)階，時(shí)空域2M(M=5)階和混合網(wǎng)格2M+N型(M=5;N=1,2,…,6)有限差分格式的頻散曲線；計(jì)算參數(shù)均為M=5,v=3000 m·s-1,τ=0.001 s,h=10 m.通過頻散曲線分析和對比可以看出：

① 傳統(tǒng)2M階(M=5)有限差分格式，空間頻散最小，但是時(shí)間頻散最大，相速度與真實(shí)速度之比約等于1時(shí)對應(yīng)1/G的區(qū)間范圍很小，大致為(0,0.075).

圖3　頻散關(guān)系曲線(a) 傳統(tǒng)2M階； (b) 時(shí)空域2M階； (c)—(h) 混合網(wǎng)格2M+N(N=1,2,…,6)型.Fig.3　Dispersion curves(a) Traditional 2Mth-order; (b) Time-space domain 2Mth-order; (c)—(h) Mixed-grid 2M+N (N=1,2,…,6) style.

② 時(shí)空域2M階(M=5)有限差分格式，與傳統(tǒng)2M階(M=5)有限差分格式相比，時(shí)間頻散減小，但空間頻散增大，頻散曲線仍比較發(fā)散，相速度與真實(shí)速度之比約等于1時(shí)對應(yīng)1/G的區(qū)間范圍略有增大，大致為(0,0.125).因此，時(shí)空域2M階有限差分格式比傳統(tǒng)2M階有限差分格式更能有效地壓制數(shù)值頻散，模擬精度更高.

③ 混合網(wǎng)格2M+N型(M=5;N=1,2,…，6)有限差分格式，與傳統(tǒng)2M階(M=5)和時(shí)空域2M階(M=5)有限差分格式相比，頻散曲線較收斂，時(shí)間頻散和空間頻散都減小，相速度與真實(shí)速度之比約等于1時(shí)對應(yīng)1/G的區(qū)間范圍增大，大致為(0,0.25)，這一區(qū)間是時(shí)空域2M階(M=5)有限差分格式的2倍，為傳統(tǒng)2M階(M=5)有限差分格式的5倍，因此，混合網(wǎng)格2M+N型有限差分格式能夠更為有效地壓制數(shù)值頻散，模擬精度更高.

④ 混合網(wǎng)格2M+N型(M=5;N=1,2,…，6)有限差分格式，數(shù)值頻散曲線特征基本相同，模擬精度基本相當(dāng)，但增大N的取值會增大正演模擬的計(jì)算量，降低效率.因此，在M取值不是足夠大的情況下，建議采用混合網(wǎng)格2M+1型有限差分格式，這樣既可以提高模擬精度，又可以保證模擬效率.

圖4a、4c分別給出了時(shí)空域2M(M=6,7)階有限差分格式的頻散曲線，這兩種差分格式的Laplace算子長度分別為25點(diǎn)和29點(diǎn).圖4b、4d分別給出了混合網(wǎng)格2M+1(M=5,6)型有限差分格式的頻散曲線，這兩種差分格式的Laplace算子長度分別為25點(diǎn)和29點(diǎn).即混合網(wǎng)格2M+1(M=5,6)型分別與時(shí)空域2M(M=6,7)階有限差分格式具有基本相同的計(jì)算量.計(jì)算頻散曲線時(shí)選取的參數(shù)均為v=3000 m·s-1,τ=0.001 s,h=10 m.

對比圖4a與4b、圖4c與4d，可以看出，混合網(wǎng)格2M+1(M=5,6)型分別比時(shí)空域2M(M=6,7)階有限差分格式數(shù)值頻散更小.因此，在計(jì)算量相同的情況下，混合網(wǎng)格2M+1型比時(shí)空域2M階有限差分格式更能有效地壓制數(shù)值頻散.

圖5a、5c分別給出了混合網(wǎng)格2M+6型(M=4,5)有限差分格式的頻散曲線，這兩種差分格式的Laplace算子長度分別為41點(diǎn)和45點(diǎn)；混合網(wǎng)格2M+6型(M=4)有限差分格式等價(jià)于四階菱形網(wǎng)格差分格式(Liu and Sen, 2013).圖5b、5d分別給出了混合網(wǎng)格2M+1(M=5,6)型有限差分格式的頻散曲線，這兩種差分格式的Laplace算子長度分別為25點(diǎn)和29點(diǎn).因此，混合網(wǎng)格2M+1(M=5,6)型分別比混合網(wǎng)格2M+6(M=4,5)型有限差分格式計(jì)算量更小，正演模擬效率更高.

圖4　時(shí)空域2M階和混合網(wǎng)格2M+1型有限差分格式頻散曲線(a)、(c) 分別為時(shí)空域2M階(M=6,7)；(b)、(d) 分別為混合網(wǎng)格2M+1型(M=5,6).Fig.4　Dispersion curves for time-space domain 2Mth-order and mixed-grid 2M+1 style finite-difference methods(a) and (c) are for time-space domain 2Mth-order (M=6,7); (b) and (d) are for mixed-grid 2M+1 (M=5,6) style.

圖5　混合網(wǎng)格2M+6型和2M+1型有限差分格式頻散曲線(a)、(c)分別為混合網(wǎng)格2M+6(M=4,5)型；(b)、(d)分別為混合網(wǎng)格2M+1(M=5,6)型.Fig.5　Dispersion curves for mixed-grid 2M+6 style and 2M+1 style finite-difference methods(a) and (c) are for mixed-grid 2M+6 (M=4,5) style; (b) and (d) are for mixed-grid 2M+1 (M=5,6) style.

對比圖5b與5a、圖5d與5c，可以看出，混合網(wǎng)格2M+1(M=5,6)型分別比混合網(wǎng)格2M+6(M=4,5)型有限差分格式數(shù)值頻散更小，正演模擬精度更高.因此，與混合網(wǎng)格2M+6型有限差分格式相比，混合網(wǎng)格2M+1型有限差分格式通過取相對較大的M值，可以達(dá)到計(jì)算效率高，同時(shí)數(shù)值頻散小、模擬精度高的效果.圖5再次表明，在M值不是足夠大的情況下，應(yīng)該采用混合網(wǎng)格2M+1型有限差分格式；菱形網(wǎng)格不是最佳差分格式.

3.2穩(wěn)定性分析

時(shí)間空間域有限差分法是以差分方程的解代替波動偏微分方程的解，只有離散后的差分方程的解收斂和穩(wěn)定，這種替代才是合理的.穩(wěn)定性條件描述了空間采樣間隔和時(shí)間采樣間隔與差分階數(shù)、差分系數(shù)之間必需滿足的條件關(guān)系，也稱為Courant穩(wěn)定性條件或CFL條件(Courant-Friedrichs-Lewy condition).

采用基于特征值的穩(wěn)定性分析方法(Liu and Sen, 2009, 2013)，得出混合網(wǎng)格2M+1型有限差分格式的穩(wěn)定性條件為：

(16)

利用相同的方法，可以給出傳統(tǒng)2M階、時(shí)空域2M階以及其他混合網(wǎng)格2M+N型有限差分格式的穩(wěn)定性條件和穩(wěn)定性因子S的表達(dá)式.

進(jìn)一步分析可知，混合網(wǎng)格2M+2型和2M+3型具有相同的穩(wěn)定性條件；同樣，混合網(wǎng)格2M+5型與2M+6型穩(wěn)定性條件相同.在進(jìn)行穩(wěn)定性分析時(shí)，具有相同穩(wěn)定性條件的混合網(wǎng)格2M+N型有限差分格式只給出其中一種.

圖6a—6f分別給出了傳統(tǒng)2M階、時(shí)空域2M階、混合網(wǎng)格2M+N(N=1,3,4,6)型有限差分格式穩(wěn)定性因子S隨r的變化曲線.根據(jù)式(16)可知，穩(wěn)定性條件要求r≤S，因此，穩(wěn)定性因子S隨r的變化曲線與直線r的交點(diǎn)為該有限差分格式的穩(wěn)定性條件.圖6g給出了傳統(tǒng)2M階，時(shí)空域2M階，混合網(wǎng)格2M+N(N=1,3,4,6)型有限差分格式的穩(wěn)定性條件.

圖6　穩(wěn)定性曲線(a) 傳統(tǒng)2M階； (b) 時(shí)空域2M階； (c)—(f) 混合網(wǎng)格2M+N (N=1,3,4,6)型； (g) 上述各種有限差分格式穩(wěn)定性曲線對比圖.Fig.6　Stability curves(a) Traditional 2Mth-order; (b) Time-space domain 2Mth-order; (c)—(f) Mixed-grid 2M+N (N=1,3,4,6) style; (g) Stability curves for all above finite-difference methods.

對比可以看出，傳統(tǒng)2M階有限差分格式，穩(wěn)定性最弱；時(shí)空域2M階和混合網(wǎng)格2M+1型有限差分格式穩(wěn)定性基本相當(dāng)；混合網(wǎng)格2M+3，2M+4，2M+6型有限差分格式，比時(shí)空域2M階有限差分格式穩(wěn)定性略強(qiáng).因此，混合網(wǎng)格2M+N型有限差分格式，增大N的取值，會增強(qiáng)穩(wěn)定性，但同時(shí)也會增大正演模擬的計(jì)算量.

4　數(shù)值模擬

4.1均勻介質(zhì)模型正演模擬

為了對比傳統(tǒng)2M階，時(shí)空域2M階和混合網(wǎng)格2M+1型有限差分法正演模擬的精度，設(shè)計(jì)了12 km×12 km的均勻介質(zhì)模型，波速為3000 m·s-1，模型網(wǎng)格化參數(shù)為Δx=Δz=10 m.震源位于模型中心(x=6 km,z=6 km)，震源子波為30 Hz的雷克子波，時(shí)間采樣間隔dt=0.001 s.分別用傳統(tǒng)2M(M=6)階，時(shí)空域2M(M=6)階和混合網(wǎng)格2M+1(M=5)型有限差分格式進(jìn)行正演模擬，三種有限差分格式的拉普拉斯算子長度均為25點(diǎn)，即具有基本相同的計(jì)算效率.

考慮到模型的對稱性，繪制波場快照圖時(shí)只畫出左上角1/4部分.圖7a、7b、7c分別給出了傳統(tǒng)2M(M=6)階、時(shí)空域2M(M=6)階、混合網(wǎng)格2M+1(M=5)型有限差分法在6 s時(shí)刻的波場快照.可以看出，傳統(tǒng)2M(M=6)階有限差分法存在嚴(yán)重的時(shí)間數(shù)值頻散(“波至超前”，與圖3a頻散曲線具有一致性)，時(shí)空域2M(M=6)階有限差分法仍然存在一定的數(shù)值頻散(“波至超前”和“波至拖尾”，與圖3b頻散曲線具有一致性)，而混合網(wǎng)格2M+1(M=5)型有限差分法基本不存在數(shù)值頻散，模擬精度顯著提高.

圖7d是分別將圖7a、7b、7c中虛線位置處的波場進(jìn)行波形顯示.可以看出，傳統(tǒng)2M(M=6)階有限差分法存在嚴(yán)重的數(shù)值頻散，波形發(fā)生嚴(yán)重畸變；時(shí)空域2M(M=6)階有限差分法存在一定程度的數(shù)值頻散，波形發(fā)生一定程度的畸變；混合網(wǎng)格2M+1(M=5)型有限差分法基本不存在數(shù)值頻散，波形未發(fā)生畸變.

因此，在計(jì)算效率基本相等的情況下，與傳統(tǒng)2M階和時(shí)空域2M階有限差分法相比，混合網(wǎng)格2M+1型有限差分法更能有效地壓制數(shù)值頻散，提高正演模擬精度，這與頻散分析的結(jié)果一致.

4.2層狀介質(zhì)模型正演模擬

圖8a是一個(gè)13 km×20 km的三層層狀介質(zhì)模型：h1=6 km,v1=3000 m·s-1；h2=3 km,v2=3600 m·s-1;h3=4 km,v3=4300 m·s-1；模型網(wǎng)格化參數(shù)為：Δx=Δz=10 m.震源位于z=100 m,x=11000 m處，震源子波為30 Hz的雷克子波，時(shí)間采樣間隔為dt=0.001 s.

圖7　正演模擬6s時(shí)刻波場快照及單道波形圖(a) 傳統(tǒng)2M (M=6)階; (b) 時(shí)空域2M (M=6)階; (c) 混合網(wǎng)格2M+1 (M=5)型;(d)黑線位置處波形圖.Fig.7　Forward modeling wavefield snap at 6s and single trace waveform(a) Traditional 2Mth-order (M=6); (b) Time-space domain 2Mth-order (M=6); (c) Mixed-grid 2M+1 (M=5) style; (d) Single trace waveform at the black dashed line indicated station.

圖8　層狀介質(zhì)模型及正演炮集 (a) 傳統(tǒng)2M (M=6)階；(b)時(shí)空域2M (M=6)階；(c)混合網(wǎng)格2M+1 (M=5)型.Fig.8　Layering model and forward modeling shot gather(a) Traditional 2Mth-order (M=6); (b) Time-space domain 2Mth-order (M=6); (c) Mixed-grid 2M+1(M=5) style.

圖9　波形放大顯示圖8b—8d中實(shí)線和虛線方框部分(a)—(c)分別對應(yīng)圖8b—8d中實(shí)線方框部分；(d)—(f)分別對應(yīng)圖8b—8d中虛線方框部分.Fig.9　Waveforms show the part in the solid and dashed frames in Figs.8b—8d(a)—(c) Waveforms show part in the solid frames in Figs.8b—8d, respectively. (d)—(f) Waveforms show part in the dashed frames in Figs.8b—8d, respectively.

圖10　二維鹽丘模型及正演炮集(a) 二維鹽丘模型； (b) 傳統(tǒng)2M (M=12)階； (c) 時(shí)空域2M (M=12)階； (d) 混合網(wǎng)格2M+1 (M=11)型.Fig.10　2D salt model and forward modeling shot gathers(a) 2D salt model; (b) Traditional 2Mth-order (M=12); (c) Time-space domain 2Mth-order (M=12); (d) Mixed-grid 2M+1 (M=11) style.

圖11　(a)—(c)分別放大顯示圖10b—10d中黑色方框部分Fig.11　(a)—(c) Zoom images of the part in the black frames in Figs.10b—10d, respectively

圖12　二維鹽丘模型基于混合網(wǎng)格2M+1 (M=11)型有限差分法逆時(shí)偏移結(jié)果Fig.12　2D salt model RTM result based on mixed 2M+1 (M=11) style FD method

分別用傳統(tǒng)2M(M=6)階,時(shí)空域2M(M=6)階和混合網(wǎng)格2M+1(M=5)型有限差分格式進(jìn)行正演模擬.三種有限差分格式的拉普拉斯算子長度均為25點(diǎn).

圖8b、8c、8d分別給出了傳統(tǒng)2M(M=6)階，時(shí)空域2M(M=6)階和混合網(wǎng)格2M+1型(M=5)有限差分法正演模擬的單炮記錄.正演模擬過程中，均采用NPML吸收邊界消除人工邊界反射，從單炮記錄可以看出，NPML吸收邊界有效地吸收了邊界反射能量.

圖9a、9b、9c分別是圖8b、8c、8d中實(shí)線方框部分的放大顯示.可以看出，傳統(tǒng)2M(M=6)階有限差分法，由于存在較強(qiáng)的時(shí)間數(shù)值頻散引起波形嚴(yán)重畸變；時(shí)空域2M(M=6)階有限差分法，由于存在時(shí)間數(shù)值頻散波形發(fā)生畸變；混合網(wǎng)格2M+1(M=5)型有限差分法，波形基本未發(fā)生畸變，數(shù)值頻散最小.

圖9d、9e、9f分別對圖8b、8c、8d中虛線方框部分進(jìn)行放大顯示，可以看出，傳統(tǒng)2M(M=6)階有限差分法，波形存在嚴(yán)重畸變(時(shí)間數(shù)值頻散引起)；時(shí)空域2M(M=6)階有限差分法，波形也存在較嚴(yán)重的畸變(空間數(shù)值頻散引起)；混合網(wǎng)格2M+1(M=5)型有限差分法，波形基本未發(fā)生畸變，數(shù)值頻散最小.

層狀介質(zhì)模型數(shù)值模擬試驗(yàn)表明，傳統(tǒng)2M階有限差分法數(shù)值頻散最嚴(yán)重，而且以時(shí)間頻散為主，正演模擬精度最低；時(shí)空域2M階有限差分法存在一定程度的時(shí)間數(shù)值頻散和空間數(shù)值頻散，正演模擬精度相對較高；混合網(wǎng)格2M+1型有限差分法能夠進(jìn)一步減小時(shí)間頻散和空間頻散，正演模擬精度最高.

4.3二維鹽丘模型正演模擬及逆時(shí)偏移

為了進(jìn)一步驗(yàn)證本文提出的混合網(wǎng)格2M+N有限差分格式對于復(fù)雜介質(zhì)模型的適用性，利用(SEG/EAGE)鹽丘模型進(jìn)行正演模擬，模擬過程中對實(shí)際的鹽丘模型進(jìn)行了適當(dāng)?shù)母脑?，將鹽丘的速度變?yōu)樵瓉淼?.2倍，模型其他部分的速度變?yōu)樵瓉淼?倍.這種改造主要是考慮到，原始的鹽丘模型速度較小，要想觀測出混合網(wǎng)格2M+1型相比時(shí)空域2M型有限差分格式在壓制數(shù)值頻散，提高正演模擬精度方面的優(yōu)勢，需要將模型設(shè)計(jì)得非常大，正演模擬時(shí)間很長；鹽丘部分速度變?yōu)樵瓉淼?.2倍(而不是將整個(gè)速度模型變?yōu)樵瓉淼?倍)主要是為了兼顧有限差分正演模擬的穩(wěn)定性.

圖10a給出了鹽丘速度模型，模型網(wǎng)格化參數(shù)為nz×nx=419×1351，Δx=Δz=20 m；正演模擬震源位于z=60 m,x=6740 m處，震源子波為主頻25 Hz的雷克子波，時(shí)間采樣間隔dt=0.0015 s.

分別用傳統(tǒng)2M(M=12)階,時(shí)空域2M(M=12)階和混合網(wǎng)格2M+1(M=11)型有限差分格式進(jìn)行正演模擬.三種有限差分格式的拉普拉斯算子長度均為49點(diǎn)，因此這三種有限差分格式的計(jì)算量基本相同.正演模擬和逆時(shí)偏移中，均采用NPML吸收邊界.

圖10b、10c、10d分別給出了傳統(tǒng)2M(M=12)階，時(shí)空域2M(M=12)階和混合網(wǎng)格2M+1型(M=11)有限差分法正演模擬的單炮記錄.圖11b、11c、11d分別是圖10b、10c、10d中黑色方框部分的放大顯示.可以看出：傳統(tǒng)2M(M=12)階有限差分法，存在嚴(yán)重的數(shù)值頻散(主要為時(shí)間數(shù)值頻散)；時(shí)空域2M(M=12)階有限差分法，也存在一定的數(shù)值頻散(主要為空間數(shù)值頻散)；混合網(wǎng)格2M+1(M=11)型有限差分法，無明顯數(shù)值頻散.

鹽丘模型的數(shù)值模擬試驗(yàn)進(jìn)一步驗(yàn)證了前面數(shù)值頻散分析和均勻介質(zhì)、層狀介質(zhì)模型數(shù)值模擬得出的結(jié)論：傳統(tǒng)2M階有限差分法數(shù)值頻散最大，正演模擬精度最低；時(shí)空域2M階有限差分法存在一定的數(shù)值頻散，模擬精度中等；混合網(wǎng)格2M+1型有限差分法無明顯數(shù)值頻散，模擬精度最高.

進(jìn)一步將混合網(wǎng)格2M+1(M=11)型有限差分法推廣應(yīng)用于逆時(shí)偏移，正演模擬參數(shù)與上面給出的參數(shù)相同.總共正演225炮，每炮1351道，炮間距120 m，道間距20 m.為了壓制逆時(shí)偏移中的低頻噪聲同時(shí)保護(hù)成像的振幅和相位信息，首先對成像前模擬數(shù)據(jù)在時(shí)間域積分，然后對成像結(jié)果進(jìn)行拉普拉斯濾波(劉紅偉等，2010b; 嚴(yán)紅勇和劉洋, 2013).圖12給出了逆時(shí)偏移成像結(jié)果，成像精度較高.

5　結(jié)論與討論

本文將頻率-空間域混合網(wǎng)格有限差分正演模擬的思想引入到時(shí)空域，提出時(shí)空域混合網(wǎng)格2M+N型有限差分格式，并導(dǎo)出了基于時(shí)空域頻散關(guān)系的差分系數(shù)計(jì)算方法.在此基礎(chǔ)上進(jìn)行頻散分析，穩(wěn)定性分析，正演模擬和逆時(shí)偏移，結(jié)果表明：

(1) 計(jì)算量基本相等時(shí)，傳統(tǒng)2M階有限差分格式數(shù)值頻散嚴(yán)重，主要為時(shí)間頻散，正演模擬精度最低；時(shí)空域2M階有限差分格式存在一定程度的時(shí)間頻散和空間頻散，正演模擬精度相對較高；混合網(wǎng)格2M+N型有限差分格式數(shù)值頻散最小，正演模擬精度最高.

(2) 混合網(wǎng)格2M+N型有限差分格式，M取值不是足夠大時(shí)，增大N的取值，對改善頻散特性，提高正演模擬精度效果微弱，但會增大計(jì)算量，降低正演模擬效率.因此，我們建議采用混合網(wǎng)格2M+1型有限差分格式；菱形網(wǎng)格是混合網(wǎng)格2M+N型(N取值較大)的特殊形式，所以菱形網(wǎng)格有限差分格式并不是最優(yōu)選擇.

(3) 混合網(wǎng)格2M+N型有限差分格式，比傳統(tǒng)2M階有限差分格式穩(wěn)定性強(qiáng)，與時(shí)空域2M階有限差分格式具有基本相同的穩(wěn)定性，增大N的取值，可以在一定程度上增強(qiáng)其穩(wěn)定性，但計(jì)算量也會相應(yīng)增大.

(4) 均勻介質(zhì)、層狀介質(zhì)和復(fù)雜鹽丘模型正演模擬結(jié)果，驗(yàn)證了混合網(wǎng)格2M+N型有限差分格式能有效壓制數(shù)值頻散，提高正演模擬精度.鹽丘模型逆時(shí)偏移進(jìn)一步驗(yàn)證了混合網(wǎng)格2M+N型有限差分格式具有普遍適用性.

本文提出的混合網(wǎng)格2M+N型有限差分格式能夠有效地提高正演模擬精度，但仍存在一定的局限性：(1) 混合網(wǎng)格2M+N型有限差分格式只適用于空間采樣Δx=Δz的情況，對Δx≠Δz情況需要進(jìn)一步研究.(2) 混合網(wǎng)格2M+N型有限差分格式從聲波方程導(dǎo)出，但對各向異性、彈性波方程的適用性未加研究.上述問題也是下一步的研究方向.

附錄A混合網(wǎng)格2M+2型有限差分格式及其差分系數(shù)求解方法

按照圖2b給出的混合網(wǎng)格2M+2型有限差分格式對拉普拉斯算子進(jìn)行差分離散近似，可以得到：

(A1)

將式(A1)和(2)代入聲波方程(1)，可以得到：

(A2)

將平面波解式(7)代入式(A2)，可以得到：

(A3)

對式(A3)中的余弦函數(shù)進(jìn)行泰勒展開，可以得到：

(A4)

(A5)

(A6)

(A7)

方程(A5)(A6)和(A7)總共M+2個(gè)方程，剛好可以求解出M+2個(gè)差分系數(shù)a1,a2,…,am;a1,1,a2,1.

Alford R M, Kelly K R, Boore D M. 1974. Accuracy of finite-difference modeling of the acoustic wave equation.Geophysics, 39(6): 834-842, doi: 10.1190/1.1440470.

Alterman Z, Karal F C. 1968. Propagation of elastic waves in layered media by finite difference methods.BulletinoftheSeismologicalSocietyofAmerica, 58(1): 367-398.

Cao S H, Chen J B. 2012. A 17-point scheme and its numerical implementation for high-accuracy modeling of frequency-domain acoustic equation.ChineseJ.Geophys. (in Chinese), 55(10): 3440-3449, doi: 10.6038/j.issn.0001-5733.2012.10.027.

Chen H M, Zhou H, Li Y Q. 2014. Application of unsplit convolutional perfectly matched layer for scalar arbitrarily wide-angle wave equation.Geophysics, 79(6): T313-T321, doi: 10.1190/geo2014-0103.1.

Chen J B. 2007. High-order time discretizations in seismic modeling.Geophysics, 72(5): SM115-SM122, doi: 10.1190/1.2750424.

Chen J Y. 2011. Application of the nearly perfectly matched layer for seismic wave propagation in 2D homogeneous isotropic media.GeophysicalProspecting, 59(4): 662-672, doi: 10.1111/j.1365-2478.2011.00949.x.

Chen J B. 2012. An average-derivative optimal scheme for frequency-domain scalar wave equation.Geophysics, 77(6): T201-T210, doi: 10.1190/geo2011-0389.1.

Collino F, Tsogka C. 2001. Application of the perfectly matched absorbing layer model to the linear elastodynamic problem in anisotropic heterogeneous media.Geophysics, 66(1): 294-307, doi: 10.1190/1.1444908.

Dablain M A. 1986. The application of high-order differencing to the scalar wave equation.Geophysics, 51(1): 54-66, doi: 10.1190/1.1442040. Etgen J, Gray S H, Zhang Y. 2009. An overview of depth imaging in exploration geophysics.Geophysics, 74(6): WCA5-WCA17, doi: 10.1190/1.3223188.

Fang G, Fomel S, Du Q Z, et al. 2014. Lowrank seismic-wave extrapolation on a staggered grid.Geophysics, 79(3): T157-T168, doi: 10.1190/geo2013-0290.1.

Fomel S, Ying L X, Song X L. 2013. Seismic wave extrapolation using lowrank symbol approximation.GeophysicalProspecting, 61(3): 526-536, doi: 10.1111/j.1365-2478.2012.01064.x. Gerhard P, Shin C, Hicks. 1998. Gauss-Newton and full Newton methods in frequency-space seismic waveform inversion.GeophysicalJournalInternational, 133(2): 341-362, doi: 10.1046/j.1365-246X.1998.00498.x. Hu W Y, Abubakar A, Habashy T M. 2007. Application of the nearly perfectly matched layer in acoustic wave modeling.Geophysics, 72(5): SM169-SM175, doi: 10.1190/1.2738553.

Jo C H, Shin C, Suh J H. 1996. An optimal 9-point, finite-difference, frequency-space, 2-D scalar wave extrapolator.Geophysics, 61(2): 529-537, doi: 10.1190/1.1443979.

Jones I F. 2008. A modeling study of preprocessing considerations for reverse-time migration.Geophysics, 73(6): T99-T106, doi: 10.1190/1.2981183.

Kelly K R, Ward R W, Treitel S, et al. 1976. Synthetic seismograms: a finite-difference approach.Geophysics, 41(1): 2-27, doi: 10.1190/1.1440605.

Komatitsch D, Martin R. 2007. An unsplit convolutional perfectly matched layer improved at grazing incidence for the seismic wave equation.Geophysics, 72(5): SM155-SM167, doi: 10.1190/1.2757586. Levander A R. 1988. Fourth-order finite-difference P-SV seismograms.Geophysics, 53(11): 1425-1436, doi: 10.1190/1.1442422. Liu H W, Li B, Liu H, et al. 2010a. The algorithm of high order finite difference pre-stack reverse time migration and GPU implementation.ChineseJ.Geophys. (in Chinese), 53(7): 1725-1733, doi: 10.3969/j.issn.0001-5733.2010.07.024.

Liu H W, Liu H, Zou Z. 2010b. The problems of denoise and storage in seismic reverse time migration.ChineseJ.Geophys. (in Chinese), 53(9): 2171-2180, doi: 10.3969/j.issn.0001-5733.2010.09.017.

Liu L, Liu H, Liu H W. 2013. Optimal 15-point finite difference forward modeling in frequency-space domain.ChineseJ.Geophys. (in Chinese), 56(2): 644-652, doi: 10.6038/cjg20130228.Liu Y, Li C C, Mou Y G. 1998. Finite-difference numerical modeling of any even-order accuracy.OilGeophysicalProspecting(in Chinese), 33(1): 1-10.

Liu Y, Sen M K. 2009. A new time-space domain high-order finite-difference method for the acoustic wave equation.JournalofComputationalPhysics, 228(23): 8779-8806, doi: 10.1016/j.jcp.2009.08.027. Liu Y, Sen M K. 2011. Scalar wave equation modeling with time-space domain dispersion-relation-based staggered-grid finite-difference schemes.BulletinoftheSeismologicalSocietyofAmerica, 101(1): 141-159.

Liu Y, Sen M K. 2013. Time-space domain dispersion-relation-based finite-difference method with arbitrary even-order accuracy for the 2D acoustic wave equation.JournalofComputationalPhysics, 232(1): 327-345, doi: 10.1016/j.jcp.2012.08.025. Nihei T, Ishii K. 2003. A fast solver of the shallow water equations on a sphere using a combined compact difference scheme.JournalofComputationalPhysics, 187(2): 639-659, doi: 10.1016/S0021-9991(03)00152-9. Saenger E H, Gold N, Shapiro S A. 2000. Modeling the propagation of elastic waves using a modified finite-difference grid.WaveMotion, 31(1): 77-92, doi: 10.1016/S0165-2125(99)00023-2. Saenger E H, Bohlen T. 2004. Finite-difference modeling of viscoelastic and anisotropic wave propagation using the rotated staggered grid.Geophysics, 69(2): 583-591, doi: 10.1190/1.1707078.

Sherer S E, Scott J N. 2005. High-order compact finite-difference methods on general overset grids.JournalofComputationalPhysics, 210(2): 459-496, doi: 10.1016/j.jcp.2005.04.017.

Shin C, Sohn H. 1998. A frequency-space 2-D scalar wave extrapolator using extended 25-point finite-difference operator.Geophysics, 63(1): 289-296, doi: 10.1190/1.1444323.

Sun R, McMechan G A. 2001. Scalar reverse-time depth migration of prestack elastic seismic data.Geophysics, 66(5): 1519-1527, doi: 10.1190/1.1487098.

Sun R, McMechan G A, Lee C S, et al. 2006. Prestack scalar reverse-time depth migration of 3D elastic seismic data.Geophysics, 71(5): S199-S207, doi: 10.1190/1.2227519.

Tang X D, Liu H, Zhang H. 2015. Frequency-space domain acoustic modeling based on an average-derivative and GPU implementation.ChineseJ.Geophys. (in Chinese), 58(4): 1341-1354, doi: 10.6038/cjg20150421.

Tang X D, Liu H, Zhang H, et al. 2015. An adaptable 17-point scheme for high-accuracy frequency-domain acoustic wave modeling in 2D constant density media.Geophysics, 80(6): T211-T221, doi: 10.1190/geo2014-0124.1.

Tian Z F, Dai S Q. 2007. High-order compact exponential finite difference methods for convection-diffusion type problems.JournalofComputationalPhysics, 220(2): 952-974, doi: 10.1016/j.jcp.2006.06.001.

Virieux J, Operto S. 2009. An overview of full-waveform inversion in exploration geophysics.Geophysics, 74(6): WCC1-WCC26, doi: 10.1190/1.3238367.

Yan H Y, Liu Y. 2012. Rotated staggered grid high-order finite-difference numerical modeling for wave propagation in viscoelastic TTI media.ChineseJ.Geophys. (in Chinese), 55(4): 1354-1365, doi: 10.6038/j.issn.0001-5733.2012.04.031. Yan H Y, Liu Y. 2013. Acoustic prestack reverse time migration using the adaptive high-order finite-difference method in time-space domain.ChineseJ.Geophys. (in Chinese), 56(3): 971-984, doi: 10.6038/cjg20130325.

Zhang H, Liu H, Tang X D, et al. 2015. Forward modeling of VTI media in the frequency-space domain based on an average-derivative optimal method.ChineseJ.Geophys. (in Chinese), 58(9): 3306-3316, doi: 10.6038/cjg20150924.

附中文參考文獻(xiàn)

曹書紅, 陳景波. 2012. 聲波方程頻率域高精度正演的17點(diǎn)格式及數(shù)值實(shí)現(xiàn). 地球物理學(xué)報(bào), 55(10): 3440-3449, doi: 10.6038/j.issn.0001-5733.2012.10.027.

劉紅偉, 李博, 劉洪等. 2010a. 地震疊前逆時(shí)偏移高階有限差分算法及GPU實(shí)現(xiàn). 地球物理學(xué)報(bào), 53(7): 1725-1733, doi: 10.3969/j.issn.0001-5733.2010.07.02.

劉紅偉, 劉洪, 鄒振等. 2010b. 地震疊前逆時(shí)偏移中的去噪與存儲. 地球物理學(xué)報(bào), 53(9): 2171-2180, doi: 10.3969/j.issn.0001-5733.2010.09.017.

劉璐, 劉洪, 劉紅偉. 2013. 優(yōu)化15點(diǎn)頻率-空間域有限差分正演模擬. 地球物理學(xué)報(bào), 56(2): 644-652, doi: 10.6038/cjg20130228. 劉洋, 李承楚, 牟永光. 1998. 任意偶數(shù)階精度有限差分法數(shù)值模擬. 石油地球物理勘探, 33(1): 1-10.

唐祥德, 劉洪, 張衡. 2015. 基于加權(quán)平均導(dǎo)數(shù)的頻率-空間域正演模擬及GPU實(shí)現(xiàn). 地球物理學(xué)報(bào), 58(4): 1341-1354, doi: 10.6038/cjg20150421.

嚴(yán)紅勇, 劉洋. 2012. 黏彈TTI介質(zhì)中旋轉(zhuǎn)交錯(cuò)網(wǎng)格高階有限差分?jǐn)?shù)值模擬. 地球物理學(xué)報(bào), 55(4): 1354-1365, doi: 10.6038/j.issn.0001-5733.2012.04.031.

嚴(yán)紅勇, 劉洋. 2013. 基于時(shí)空域自適應(yīng)高階有限差分的聲波疊前逆時(shí)偏移. 地球物理學(xué)報(bào), 56(3): 971-984, doi: 10.6038/cjg20130325.

張衡, 劉洪, 唐祥德等. 2015. 基于平均導(dǎo)數(shù)優(yōu)化方法的VTI介質(zhì)頻率空間域正演. 地球物理學(xué)報(bào), 58(9): 3306-3316, doi: 10.6038/cjg20150924.

(本文編輯何燕)

Scalar wave equation modeling using the mixed-grid finite-difference method in the time-space domain

HU Zi-Duo1,2, HE Zhen-Hua1, LIU Wei2, WANG Yu-Chao2, HAN Ling-He2, WANG Shu-Jiang2, YANG Zhe2

1StateKeyLaboratoryofOilandGasReservoirGeologyandExploitation，ChengduUniversityofTechnology，Chengdu610059,China2ResearchInstituteofPetroleumExploration&Development-Northwest,PetroChina,Lanzhou730020,China

The finite-difference method has been widely used in seismic forward modeling, RTM (Reverse Time Migration) and FWI (Full Waveform Inversion) because of its easy implementation, small memory and low computational cost. Numerical dispersion, due to discretization of time and space derivatives, seriously affects the forward modeling accuracy of the finite-difference method. So suppressing the numerical dispersion to improve the forward modeling accuracy is the key problem for the finite-difference method.

Mixed-grid; Time-space domain; Finite-difference; Dispersion analysis; Forward modeling

10.6038/cjg20161027.

國家油氣專項(xiàng)“天然氣地球物理烴類檢測、評價(jià)技術(shù)及應(yīng)用”(2016ZX05007-006)資助.

胡自多，男，1971年生，成都理工大學(xué)在讀博士，主要從事地震波正演模擬與偏移成像、地震資料處理等方面的研究工作.

E-mail：huzd@petrochina.com.cn

10.6038/cjg20161027

P631

2015-12-25，2016-07-05收修定稿

胡自多， 賀振華， 劉威等. 2016. 旋轉(zhuǎn)網(wǎng)格和常規(guī)網(wǎng)格混合的時(shí)空域聲波有限差分正演. 地球物理學(xué)報(bào)，59(10):3829-3846,

Hu Z D, He Z H, Liu W, et al. 2016. Scalar wave equation modeling using the mixed-grid finite-difference method in the time-space domain.ChineseJ.Geophys. (in Chinese),59(10):3829-3846,doi:10.6038/cjg20161027.