旋轉網格和常規(guī)網格混合的時空域聲波有限差分正演

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

1 成都理工大學油氣藏地質及開發(fā)工程國家重點實驗室， 成都　610059 2 中國石油勘探開發(fā)研究院西北分院， 蘭州　730020

?

旋轉網格和常規(guī)網格混合的時空域聲波有限差分正演

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

1 成都理工大學油氣藏地質及開發(fā)工程國家重點實驗室， 成都610059 2 中國石油勘探開發(fā)研究院西北分院， 蘭州730020

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

有限差分正演模擬方法具有計算效率高、內存占用小、簡單易實現(xiàn)的優(yōu)點，因而廣泛應用于地震波正演模擬(Alterman and Karal, 1968; Alford et al., 1974; Kelly et al., 1976)、逆時偏移(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ù)進行網格差分近似造成的，時間數(shù)值頻散通常使相速度變大而出現(xiàn)“波至超前”，空間數(shù)值頻散使相速度變小而出現(xiàn)“波至拖尾”，并且頻率越大，數(shù)值頻散誤差越大(Dablain, 1986)，嚴重影響了地震波高頻成分的模擬精度.因此，壓制數(shù)值頻散一直是有限差分正演模擬的重要研究內容.

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

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

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

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

吸收邊界是地震波數(shù)值模擬中為了消除人工邊界反射而引入的.一階應力速度聲波方程正演模擬時，普遍采用分裂形式的PML吸收邊界(Collino and Tsogka, 2001).二階標量聲波方程進行正演模擬時，采用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吸收邊界具有輔助變量少，易于實現(xiàn)的優(yōu)點.本文提出的混合網格2M+N型有限差分法，基于二階標量聲波方程，因此采用NPML吸收邊界.

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

2　方法原理

2.1聲波方程混合網格2M+N型有限差分格式構建

二維時間-空間域常密度標量聲波方程可寫為：

(1)

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

(2)

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

(3)

圖1　傳統(tǒng)2M階和時空域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階和時空域2M階聲波方程有限差分法的離散化方程，但兩種方法求取差分系數(shù)的方法不同.傳統(tǒng)2M階有限差分法僅基于空間域頻散關系計算差分系數(shù)，而時空域2M階有限差分法則基于時空域頻散關系計算差分系數(shù).

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

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

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

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

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

(6)

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

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

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

圖2　混合網格2M+N型有限差分格式 (a)—(f)分別為混合網格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階(時空域2M階)混合網格2M+1型混合網格2M+2型混合網格2M+3型混合網格2M+4型混合網格2M+5型混合網格2M+6型拉普拉斯算子長度4M+14M+54M+94M+134M+174M+214M+25

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

方程(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ù)進行泰勒級數(shù)展開，ωτ=kvτ=khr，可以得到：

(10)

(11)

(12)

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

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

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

表2　傳統(tǒng)2M階，時空域2M階和混合網格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頻散分析

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

(13)

其中vph為相速度，v為地震波在介質傳播中的真實速度.

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

(14)

(15)

其中G=λ/h，λ為波長，則G表示單位波長內網格點的個數(shù).

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

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

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

圖3　頻散關系曲線(a) 傳統(tǒng)2M階； (b) 時空域2M階； (c)—(h) 混合網格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.

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

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

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

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

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

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

圖4　時空域2M階和混合網格2M+1型有限差分格式頻散曲線(a)、(c) 分別為時空域2M階(M=6,7)；(b)、(d) 分別為混合網格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　混合網格2M+6型和2M+1型有限差分格式頻散曲線(a)、(c)分別為混合網格2M+6(M=4,5)型；(b)、(d)分別為混合網格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，可以看出，混合網格2M+1(M=5,6)型分別比混合網格2M+6(M=4,5)型有限差分格式數(shù)值頻散更小，正演模擬精度更高.因此，與混合網格2M+6型有限差分格式相比，混合網格2M+1型有限差分格式通過取相對較大的M值，可以達到計算效率高，同時數(shù)值頻散小、模擬精度高的效果.圖5再次表明，在M值不是足夠大的情況下，應該采用混合網格2M+1型有限差分格式；菱形網格不是最佳差分格式.

3.2穩(wěn)定性分析

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

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

(16)

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

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

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

圖6　穩(wěn)定性曲線(a) 傳統(tǒng)2M階； (b) 時空域2M階； (c)—(f) 混合網格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)定性最弱；時空域2M階和混合網格2M+1型有限差分格式穩(wěn)定性基本相當；混合網格2M+3，2M+4，2M+6型有限差分格式，比時空域2M階有限差分格式穩(wěn)定性略強.因此，混合網格2M+N型有限差分格式，增大N的取值，會增強穩(wěn)定性，但同時也會增大正演模擬的計算量.

4　數(shù)值模擬

4.1均勻介質模型正演模擬

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

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

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

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

4.2層狀介質模型正演模擬

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

圖7　正演模擬6s時刻波場快照及單道波形圖(a) 傳統(tǒng)2M (M=6)階; (b) 時空域2M (M=6)階; (c) 混合網格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　層狀介質模型及正演炮集 (a) 傳統(tǒng)2M (M=6)階；(b)時空域2M (M=6)階；(c)混合網格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中實線和虛線方框部分(a)—(c)分別對應圖8b—8d中實線方框部分；(d)—(f)分別對應圖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) 時空域2M (M=12)階； (d) 混合網格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　二維鹽丘模型基于混合網格2M+1 (M=11)型有限差分法逆時偏移結果Fig.12　2D salt model RTM result based on mixed 2M+1 (M=11) style FD method

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

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

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

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

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

4.3二維鹽丘模型正演模擬及逆時偏移

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

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

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

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

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

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

5　結論與討論

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

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

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

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

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

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

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

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

(A1)

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

(A2)

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

(A3)

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

(A4)

(A5)

(A6)

(A7)

方程(A5)(A6)和(A7)總共M+2個方程，剛好可以求解出M+2個差分系數(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.

附中文參考文獻

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

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

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

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

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

嚴紅勇, 劉洋. 2012. 黏彈TTI介質中旋轉交錯網格高階有限差分數(shù)值模擬. 地球物理學報, 55(4): 1354-1365, doi: 10.6038/j.issn.0001-5733.2012.04.031.

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

張衡, 劉洪, 唐祥德等. 2015. 基于平均導數(shù)優(yōu)化方法的VTI介質頻率空間域正演. 地球物理學報, 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.

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

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

E-mail：huzd@petrochina.com.cn

10.6038/cjg20161027

P631

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

胡自多， 賀振華， 劉威等. 2016. 旋轉網格和常規(guī)網格混合的時空域聲波有限差分正演. 地球物理學報，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.