基于新的差分結構的時-空域高階有限差分波動方程數值模擬方法

張保慶， 周輝， 陳漢明， 盛善波

1 油氣資源與探測國家重點實驗室， CNPC物探重點實驗室， 中國石油大學(北京)， 北京　102249 2 東方地球物理公司研究院大港分院， 天津大港　300280 3 中國石油天然氣勘探開發(fā)公司， 北京　100034

基于新的差分結構的時-空域高階有限差分波動方程數值模擬方法

張保慶1, 2， 周輝1*， 陳漢明1， 盛善波3

1 油氣資源與探測國家重點實驗室， CNPC物探重點實驗室， 中國石油大學(北京)， 北京102249 2 東方地球物理公司研究院大港分院， 天津大港300280 3 中國石油天然氣勘探開發(fā)公司， 北京100034

摘要傳統(tǒng)的高階有限差分波動方程數值模擬方法采用高階差分算子近似空間偏導數，能有效抑制空間頻散.然而，傳統(tǒng)的有限差分法僅采用二階差分算子近似時間偏導數，這使得地震波場沿時間外推的精度較低.當采用較大的時間采樣間隔，傳統(tǒng)的有限差分法模擬波場會出現明顯的時間頻散，甚至不穩(wěn)定.本文基于新的差分結構和中心網格剖分，發(fā)展了一種空間任意偶數階精度、時間四階和六階精度的時空域有限差分方法.基于對離散后的頻散關系進行泰勒展開，本文推導了時空域高階有限差分算子的差分系數.相速度分析表明時間四階、六階精度的差分方法能顯著地減小傳統(tǒng)時間二階精度差分方法的時間頻散.在相同的精度下與傳統(tǒng)差分法比較，本文發(fā)展的時間四階、六階有限差分方法的計算效率比傳統(tǒng)方法高.均勻和非勻均介質中的波場數值模擬實驗進一步證實本文研究的時空高階有限差分方法的優(yōu)越性.

關鍵詞波動方程； 時空域； 高階有限差分； 數值模擬

1引言

波動方程數值模擬是當前熱門的逆時偏移、最小二乘偏移和全波形反演技術的基礎.波動方程數值模擬的精度和效率在很大程度上決定了波動方程偏移和反演的精度和效率.當前，波動方程數值模擬的方法主要包括基于傅里葉變換的偽譜法，基于非規(guī)則網格剖分的有限元類方法以及基于差分近似的有限差分方法.偽譜法采用傅里葉變換計算空間偏導數，能完全壓制由于空間離散造成的數值頻散，但其缺點是計算效率較低(謝桂生等，2005)，并且由于傅里葉變換的全局性，吸收邊界的處理較為麻煩.有限元類方法剖分網格更靈活，能準確地模擬各種不規(guī)則邊界，從而避免矩形網格剖分造成的階梯繞射.然而，受制于巨大的計算量，目前有限元類方法在逆時偏移和全波形反演中很少使用.相比之下，有限差分法無疑是一種更靈活簡便的數值算法，已被廣泛用于各類偏微分方程的數值求解，在波動方程偏移和反演中也大受歡迎.

Alterman and Karal(1968)最早采用有限差分法模擬二維層狀介質中的彈性波場，Kelly et al.(1976)則討論了彈性波動方程數值模擬的兩種有限差分方法，分別對應中心網格有限差分法和交錯網格有限差分法.Virieux(1984, 1986)正式將Yee(1966)的交錯網格離散方法引入到波動方程數值模擬中，并求解一階速度-應力形式的聲波和彈性波動方程.早期的交錯網格有限差分法僅采用二階差分算子分別近似空間和時間偏導數項，當采用較大的網格間距時，數值頻散嚴重，需要校正(吳國忱和王華忠，2005).為了提高有限差分數值模擬的精度，Levander(1988)采用四階精度的差分算子近似空間偏導數，顯著地改善了二階差分算子的計算精度，形成了空間四階精度、時間二階精度的差分方法，本文將其簡記為(4,2)差分法.與二階差分方法相比，(4,2)差分方法更好地平衡了計算精度和計算效率，因此被廣泛用于地震波場的數值模擬(Graves, 1996；Moczo et al., 2000；Shi et al., 2012).為了進一步提高空間差分算子的精度，許多學者進一步討論了高階有限差分數值模擬方法(Holberg, 1987；Fornberg, 1987；Kristek et al. 2010；Dong et al., 2013；杜啟振等，2015).通過增加差分算子的長度，空間精度可以達到任意偶數(2M)階(劉洋等，1988)，本文將傳統(tǒng)的高階有限差分方法記為(2M,2).

盡管傳統(tǒng)的高階有限差分方法能有效地抑制空間頻散，但長期以來很少有學者關注時間頻散.事實上，地震波在時間和空間里傳播，因網格離散造成的數值頻散同時包括空間頻散和時間頻散，傳統(tǒng)的高階差分(2M,2)方法時間外推的精度只有二階，當采用較大的時間步長時會出現明顯的時間頻散，甚至不穩(wěn)定.近些年，許多旨在提高地震波場時間外推精度的數值模擬方法得以發(fā)展.這些方法主要包括基于傅立葉變換的矩陣低秩分解方法(Fomel et al., 2013；Song et al., 2013；Fang et al., 2014)和基于頻散關系的時-空域有限差分方法(Finkelstein and Kastner, 2007；Tan and Huang, 2014).基于矩陣低秩分解的數值模擬方法在波場外推的每個時間步至少需要兩次反傅里葉變換，計算量較大.因此，Song等(2013)和Fang等(2014)在矩陣低秩分解的基礎上分別發(fā)展了基于中心網格和基于交錯網格的有限差分方法，進一步提高了計算效率.Tan和Huang(2014)基于新的差分結構發(fā)展了時間四階和六階精度、空間任意偶數階精度的交錯網格有限差分方法，分別記為(2M,4)和(2M,6).頻散分析表明(2M,4)和(2M,6)交錯網格有限差分方法能顯著地提高傳統(tǒng)(2M,2)交錯網格有限差分方法的精度，并且比傳統(tǒng)方法具有更好的穩(wěn)定性.國內學者董良國等(2000)在求解彈性波動方程的一階速度-應力形式時將速度分量對時間的奇數階偏導數轉化成應力分量對空間的偏導數，發(fā)展了時間四階精度、空間任意偶數階精度的交錯網格有限差分法.楊慶節(jié)等(2014)推導了可導函數任意次導數的任意階精度的差分近似及相應的差分系數，進一步完善了董良國等(2000)的高階差分法.

高階差分算子的差分系數可利用泰勒多項式展開進行確定，也可采用一些優(yōu)化的方法確定差分系數，如基于最小二乘的有限差分方法(Liu et al., 2014; 張杰和吳國忱，2015)，基于正則化的有限差分方法(Wang et al., 2014)以及基于非線性優(yōu)化的有限差分方法(梁文全等，2015).時-空域有限差分方法的差分系數隨速度而變換，當速度劇烈變化時，基于優(yōu)化的方法需要對每個不同的速度值進行一次優(yōu)化，引入了額外的計算量.相比之下，基于泰勒展開的方法給出了有限差分系數的解析表達式，使得計算差分系數的計算量基本可以忽略.

本文基于新的差分結構，發(fā)展了基于中心網格剖分的空間任意2M階，時間四階和六階精度的有限差分方法，與Tan和Huang(2014)討論的(2M,4)和(2M,6)交錯網格有限差分法相比，本文的差分方法適用于離散二階偏導數，因此適用于模擬二階波動方程.當不考慮介質的密度變化時，本文研究的(2M,4)和(2M,6)中心網格有限差分方法具有更高的計算效率.

2波動方程的差分離散格式

2.1時間四階精度的差分方法

地震波在二維聲學介質中傳播的控制方程為

(1)

其中t表示時間，xoz組成二維(2D)笛卡爾坐標系，v為地震波的傳播速度，p為壓力.采用如下差分格式離散方程(1)，公式為

(2)

其中上標表示時間網格點的標號，下標表示空間網格點的標號，a0,0,a1,1和am,0(m=1,2,…,M)表示差分系數，h為網格大小，Δt為時間步長，M為空間差分算子長度.方程(2)為一種新的差分結構，與傳統(tǒng)的差分結構相比，額外包含了非軸向的四個空間網格節(jié)點.Tan和Huang(2014)的研究表明，額外包含非軸向的四個網格節(jié)點剛好使得差分近似的截斷誤差為o(Δt4)+o(h2M)(符號o表示高階無窮小)，因此時間差分能達到四階精度，空間差分能達到2M階精度.本文將基于新的差分結構(2)并采用泰勒展開的方法確定方程(2)中的差分系數.對方程(2)兩端做傅里葉變換，并利用傅里葉變換的時移性質可得到如下頻散關系為

(3)

其中γ=vΔt/h為網格比，kx,kz為2D離散波數分量，ω為角頻率.利用ω=vk替換ω并對方程(3)中的余弦函數做泰勒展開得到：

(4)

(5)

將方程(5)代入(4)并按h2j項整理得到

(6)

對比方程(6)兩端h2j項的系數可得到：

(7)

(8)

(9)

(10)

將方程(10)代入方程(7)和(9)求得差分系數為

(11)

利用方程(11)中的差分系數，差分格式(2)能達到空間2M階精度，時間四階精度.

2.2時間六階精度的差分方法

Tan and Huang (2014)發(fā)展了時間六階精度的交錯網格有限差分法，與其時間四階精度的差分方法相比，其差分結構里又額外包含了八個網格節(jié)點.采用泰勒展開確定差分系數時，額外的網格節(jié)點使得求取差分系數的方程適定(即方程的個數剛好等于差分系數的個數)，并且差分截斷誤差為o(Δt6)+o(h2M)，因此時間差分能達到六階精度，空間差分能達到2M階精度.類似地，為發(fā)展時間六階精度的中心網格差分方法，本文采用如下的差分格式為

(12)

其中d0,0,d1,2和dm，0(m=1,2,…,M)為差分系數.類似地，對方程(12)兩端做傅里葉變換，并利用其時移性質得到如下頻散關系為

(13)

將ω=vk代入方程(13)，并對方程(13)兩端的余弦函數做泰勒展開得到

(14)

對比(14)式兩端h2j項的系數可得到：

(15)

(16)

γ2j-2, (j=1,2,…,M).

(17)

為獲得時間外推的六階精度，方程(16)兩端混

(18)

聯立方程(15)、(17)和(18)得到空間任意偶數階，時間六階精度的差分系數為

(19)

3頻散分析

根據頻散關系(3)和(13)求得角頻率ω，進而得到數值相速度為vnum=ω/k，定義數值相速度和真實地震波傳播速度之間的比值δ=vnum/v，用該比值來描述頻散程度.該比值越接近于1，頻散越小，越遠離1頻散越嚴重.另外，空間頻散造成地震波的相位滯后，即δ<1；而時間頻散則造成地震波的相位超前，即δ>1(Levander，1988；Tan and Huang，2014).據此可大致區(qū)分時間頻散和空間頻散.根據以上思路，時間四階差分算子的δ的表達式為

(20)

而時間六階差分算子的δ的表達式為

(21)

圖1顯示了當網格比取不同值時，傳統(tǒng)時間二階差分方法和本文推導的時間四階、六階差分方法的頻散誤差示意圖，空間差分算子長度統(tǒng)一取為2M=16.當網格比取較小值γ=0.2時，傳統(tǒng)的時間二階差分方法仍然存在較明顯的時間頻散，如圖1a所示.當網格比增加到γ=0.4(實際上是時間采樣步長增加一倍)，傳統(tǒng)的方法頻散誤差顯著增加，如圖1b所示.相比之下，時間四階精度的差分方法在γ=0.2時頻散誤差很小，當時間步長增加一倍時，頻散誤差也沒有顯著增加，如圖1c和圖1d所示.對比圖1d和圖1a也可得出結論，當傳統(tǒng)方法的時間步長取為時間四階精度差分方法的一半時，其頻散誤差仍然比時間四階精度的差分方法大.圖1e和圖1f則分別顯示了γ=0.2和γ=0.4時時間六階精度差分方法的頻散誤差，很明顯，時間六階精度的差分方法壓制時間頻散的效果最好.圖1g和圖1h則分別顯示了γ=0.6時本文研究的時間四階和六階精度差分算子的頻散誤差，在高波數部分四階差分算子出現了較明顯的頻散，而六階差分算子的精度更高，誤差更小.對比圖1h和圖1a可以得出結論，即使傳統(tǒng)時間二階差分方法所取的時間步長為時間六階差分方法的1/3，其頻散誤差仍然比后者嚴重.

4穩(wěn)定性分析

(22)

取尼奎斯特波數(kx,kz)=(π/h,π/h)代入方程(22)得到：

(23)

將方程(11)中a0,0的表達式代入(23)式，經過整理得到(2M, 4)差分方法的穩(wěn)定性條件為

(24)

其中int表示直接截取浮點數的整數部分.同理可得到(2M,6)差分方法的穩(wěn)定性條件為

圖1　頻散示意圖(a) 傳統(tǒng)(16,2)方法，γ=0.2; (b) 傳統(tǒng)(16,2)方法，γ=0.4; (c) 新的(16,4)方法，γ=0.2; (d) 新的(16,4)方法γ=0.4; (e) 新的(16,6)方法，γ=0.2; (f) 新的(16,6)方法γ=0.4; (g) 新的(16,4)方法γ=0.6; (h) 新的(16,6)方法γ=0.6.Fig.1　Dispersion relations(a) The traditional (16,2) scheme with γ=0.2 and (b) with γ=0.4; (c) The new (16,4) scheme with γ=0.2 and (d) with γ=0.4; (e) The new (16,6) scheme with γ=0.2 and (f) with γ=0.4, (g) The new (16,4) scheme with γ=0.6, and the new (16,6) scheme with γ=0.6.

(25)

圖2　傳統(tǒng)時間二階差分方法，新的時間四階、六階差分方法取不同差分算子長度時對應的最大網格比Fig.2　The allowed maximum Courant-Friedrichs-Lewy (CFL) number with ranging stencil length from 4 to 32, for the traditional temporal 2nd-order, and our temporal 4th- and sixth-order FD schemes.

圖2進一步比較了取不同差分算子長度時，傳統(tǒng)的時間二階差分，本文所研究的時間四階和六階差分方法能取得的最大網格比.從圖中可以看出，時間四階和六階差分方法顯著地改善了傳統(tǒng)二階差分方法的穩(wěn)定性，時間六階差分算子的穩(wěn)定性比時間四階差分算子的穩(wěn)定性略好.

5計算效率分析

時間四階和六階精度差分算子基于新的差分結構，該差分結構除了包含軸向2M+1個離散點外，還包含少數的非軸向離散點，與傳統(tǒng)的差分格式相比，引入了額外的計算量.然而，時間高階差分方法改善了傳統(tǒng)二階方法的精度，在相同的精度下比較，傳統(tǒng)方法需要使用更小的時間步長.整體而言，時間高階差分方法能提高計算效率.通過一系列的頻散分析可以得到，本文討論的時間四階差分算子采用約2.7倍于傳統(tǒng)方法的時間步長時，仍然達到相似的精度；而時間六階差分算子則可采用約4倍于傳統(tǒng)方法的時間步長時，仍然能達到可比擬的精度.表1列出了傳統(tǒng)的時間二階差分算子、本文討論的時間四階、六階差分算子包含的離散點數，每個時間步內計算空間偏導數的浮點運算次數以及達到相似精度時所采用的時間步長.當差分算子長度2M=16時，與傳統(tǒng)的時間二階差分方法相比，本文發(fā)展的時間四階差分算子和時間六階差分算子的理論加速比分別為2.4和3.0.

表1　傳統(tǒng)時間二階差分算子與時間四階、

6數值例子

本文的第一個數值模擬例子模擬均勻介質中的聲波波場，介質的地震波傳播速度為2000 m·s-1，網格大小為20 m，采用主頻為12.5 Hz的雷克子波作為震源，置于模型的中間，差分算子的長度2M=16.圖3顯示了采用不同差分方法模擬得到的3.0 s的波場切片.圖3a采用傳統(tǒng)的時間二階差分算子計算，時間步長為5.0 ms，模擬結果出現了顯著的相位扭曲，說明傳統(tǒng)的時間二階差分方法的精度低；圖3b是時間步長減半為2.5 ms時傳統(tǒng)方法模擬的結果，其數值頻散仍然比較明顯；進一步減小時間步長為2 ms，傳統(tǒng)的時間二階精度差分法的精度有所提高，但仍然可見微弱的時間頻散，如圖3c所示；圖3d顯示的波場切片為傳統(tǒng)的二階方法模擬，采用的時間步長已經縮小4倍達到1.25 ms，未見明顯的數值頻散.圖3e和3f是分別采用本文發(fā)展的時間四階、六階差分方法模擬得到的波場切片，采用時間步長為5.0 ms.對比圖3b和3e可以發(fā)現，盡管傳統(tǒng)的二階差分方法采用的時間步長為時間四階差分方法的一半，但其數值頻散仍然比四階方法嚴重；對比圖3c和3e可以發(fā)現，當四階差分方法采用的時間步長為傳統(tǒng)方法的2.5倍時，其壓制數值頻散的效果仍然更好，這與前文的頻散分析結論是一致的.進一步對比圖3d和3f，兩圖中的波場均未見明顯的數值頻散，而此時時間六階差分方法采用的時間步長為傳統(tǒng)二階方法的4倍，這與前文的頻散分析結論也是一致的.

本文的第二個數值算例模擬在二維BP鹽丘中傳播的聲波波場，圖4為相應的速度模型，該模型的最大速度差異比約為3.4.離散采用的網格大小為15 m，采用主頻為12.5 Hz的雷克子波作為震源，位于(x,z)=(7500,300) m處，數值模擬采用的差分長度為2M=16.分別采用傳統(tǒng)的時間二階差分方法、本文發(fā)展的時間四階、六階差分方法模擬聲波的傳播.為對比各種差分方法模擬的結果，將傳統(tǒng)時間二階差分方法采用微小的時間步長0.15 ms模擬的地震記錄作為參考解(見圖5).

圖3　采用不同差分方法模擬均勻介質中3 s時的波場切片，(a)傳統(tǒng)時間二階差分方法的模擬結果，時間步長5.0 ms，(b)傳統(tǒng)時間二階差分方法的模擬結果，時間步長2.5 ms，(c)傳統(tǒng)時間二階差分方法的模擬結果，時間步長 2.0 ms，(d)傳統(tǒng)時間二階差分方法的模擬結果，時間步長1.25 ms，(e)時間四階差分方法模擬的結果，時間步長5.0 ms，(f)時間六階差分方法模擬的結果，時間步長5.0 msFig.3　Acoustic snapshots at 3.0 s computed by (a) the traditional temporal 2nd-order FD scheme with the time step of 5.0 ms, (b) with the time step of 2.5 ms, (c) with the time step of 2.0 ms, and (d) with the time step of 1.25 ms, (e) by our temporal 4th-order and (f) 6th-order FD schemes with the time step of 5.0 ms

圖4　BP鹽丘模型Fig.4　BP salt model

圖5　傳統(tǒng)時間二階差分方法采用微小時間步長0.15 ms生成的參考地震記錄Fig.5　Reference seismic recording generated by traditional temporal 2nd-order FD scheme with a tiny time step of 0.15 ms

圖6　(x,z)=(12,0) km處的單道地震記錄對比，時間段4.8～5.2 s，(a)時間二階、四階和六階差分方法的模擬結果，采用相同的時間步長1.5 ms，(b)時間二階差分方法采用0.75 ms，0.5 ms時間采樣間隔的模擬結果，以及時間四階、六階差分方法采用1.5 ms采樣間隔的模擬結果Fig.6　Recorded seismograms at (x,z)=(12,0) km during 4.8～5.2 s, (a) computed by temporal 2nd-order, 4th-order and 6th-order FD schemes with the same time step of 1.5 ms, (b) computed by temporal 2nd-order FD scheme with the time step of 0.75 ms, 0.5 ms, and computed by the temporal 4th-order and 6th-order FD schemes with the time step of 1.5 ms

圖6中對比了遠偏移距接收點(x,z)=(12,0) km處的振動圖，時間區(qū)段為4.8～5.2 s.圖6a上、中、下三個子圖分別對比了傳統(tǒng)時間二階差分方法、本文發(fā)展的時間四階、六階差分方法模擬的地震道與參考道的波形差異，三種方法均采用1.5 ms的時間步長.由圖可知，傳統(tǒng)的時間二階差分方法模擬的結果與參考解存在明顯的差異，而本文發(fā)展的時間四階、六階差分方法與參考解吻合較好.圖6b上、中兩子圖顯示了傳統(tǒng)時間二階差分方法模擬的地震道，時間步長分別采用0.75 ms和0.5 ms，分別為圖6a中時間步長的1/2倍和1/3倍.圖6b中的第三個子圖則顯示了時間四階、六階差分方法模擬的地震道.仔細對比圖6b中的三個子圖，可發(fā)現采用1.5 ms時間步長的四階、六階差分方法模擬的結果幾乎一致，并且與參考解的吻合程度略好于采用0.75 ms時間步長的二階差分方法.然而，與理論的頻散分析結論有些差異的是，當傳統(tǒng)的時間二階差分方法采用的時間步長為0.5 ms時(見圖6b中第二個子圖)，其吻合參考解的效果略微勝過采用1.5 ms時間步長的六階差分方法.我們也對比了其他地震道，發(fā)現了一致的規(guī)律.這說明模擬非均勻介質中的地震波場，時間四階精度搭配空間高階精度的差分方法已經足夠，而時間六階差分算子的優(yōu)勢并未體現.

7結論

本文基于中心網格剖分發(fā)展了適用于二階聲波方程的時間四階、六階，空間任意偶數階的有限差分數值模擬方法，并基于泰勒多項式展開推導了差分系數的解析表達式.與傳統(tǒng)的時間二階差分方法相比，本文發(fā)展的時間高階差分方法具有更好的穩(wěn)定性，允許更大的時間采樣步長.頻散分析表明本文發(fā)展的時間四階、六階精度的差分方法能顯著地改善傳統(tǒng)的時間二階精度差分方法的精度，在相同的頻散誤差限下，本文發(fā)展的時間四階、六階精度差分方法允許更大的時間步長，因此提高了波場數值模擬的計算效率.在相似的頻散條件下，理論的計算量分析表明本文發(fā)展的(16,4)和(16,6)差分算子的計算效率分別是傳統(tǒng)(16,2)方法的2.4倍和3.0倍.然而，非均勻介質中波場模擬的數值實驗證實，時間六階精度差分方法模擬的結果與時間四階精度差分方法模擬的結果幾乎完全一致，考慮到時間四階差分方法包含的浮點運算更少，因此，時間四階精度搭配空間高階精度更可取.此外，本文發(fā)展的時-空高階差分方法可以很容易地推廣至三維情況.

References

Alterman Z, Karal F. 1968. Propagation of elastic waves in layered media by finite difference methods.Bull.Seism.Soc.Am., 58(1): 367-398.Dong C H, Wang S X, Zhao J G, et al. 2013. Numerical experiment and analysis of the differential acoustic resonance spectroscopy for elastic property measurements.J.Geophys.Eng., 10(5): 1-10.

Dong L G, Ma Z T, Cao J Z, et al. 2000. The staggered-grid high-order difference method of one-order elastic equation.ChineseJ.Geophys. (in Chinese), 43(6): 856-864.

Du Q Z, Guo C F, Gong X F. 2015, Hybrid PS/FD numerical simulation and stability analysis of pure P-wave propagation in VTI media.ChineseJ.Geophys. (in Chinese), 58(4): 1290-1304.Fang G, Fomel S, Du Q Z, et al. 2014. Lowrank seismic-wave extrapolation on a staggered grid.Geophysics, 79(3): T157-T168.

Finkelstein B, Kastner R. 2007. Finite difference time domain dispersion reduction schemes.J.Comput.Phys., 221(1): 422-438.Fomel S. Ying L X, Song X L. 2013. Seismic wave extrapolation using lowrank symbol approximation.Geophys.Prospect., 61(3): 526-536.

Fornberg B. 1987. The pseudospectral method: comparisons with finite differences for the elastic wave equation.Geophysics, 52(4): 483-501.

Graves R W. 1996. Simulating seismic wave propagation in 3D elastic media using staggered-grid finite differences.Bull.Seism.Soc.Am., 86(4): 1091-1106.Holberg O. 1987. Computational aspects of the choice of operator and sampling interval for numerical differentiation in large-scale simulation of wave phenomena.Geophys.Prospect., 35(6): 629-655.

Kelly K, Ward R, Treitel S, et al. 1976. Synthetic seismograms: A finite-difference approach.Geophysics, 41(1): 2-27.

Kristek J, Moczo P, Galis M. 2010. Stable discontinuous staggered grid in the finite-difference modeling of seismic motion.Geophy.J.Int., 183(3): 1401-1407.Levander A R. 1988. Fourth-order finite-difference P-SV seismograms.Geophysics, 53(11): 1425-1436.Liang W Q, Wang Y F, Yang C C. 2015. Determination on the implicit finite-difference operator based on optimization method in time-space domain.Geophy.Prosp.Petroleum, 54(3): 254-259.

Liu H F, Dai N X, Niu F, et al. 2014. An explicit time evolution method for acoustic wave propagation.Geophysics, 79(3): T117-T124.

Liu Y, Li C C, Mu Y G. 1998. Finite-difference numerical modeling of any even-order accuracy.OilGeophys.Prosp., 33(1): 1-10.Moczo P, Kristek J, Halada L.2000. 3D fourth-order staggered-grid finite-difference schemes: Stability and grid dispersion.Bull.Seism.Soc.Am., 90(3): 587-603.

Shi R Q, Wang S X, Zhao J G. 2012. An unsplit complex-frequency-shifted PML based on matched Z-transform for FDTD modelling of seismic wave equations.J.Geophys.Eng. 9(2): 218-229.

Song X L, Fomel S, Ying L X. 2013. Lowrank finite-differences and lowrank Fourier finite-differences for seismic wave extrapolation.Geophy.J.Int., 193(2): 960-969.

Tan S R, Huang L J. 2014. An efficient finite-difference method with high-order accuracy in both time and space domains for modelling scalar-wave propagation.Geophy.J.Int., 197(2): 1250-1267.

Virieux J. 1984. SH-wave propagation in heterogeneous media: Velocity-stress finite-difference method.Geophysics, 49(11): 1933-1957.

Virieux J. 1986. P-SV wave propagation in heterogeneous media: Velocity stress finite-difference method.Geophysics, 51(4): 889-901.

Wang Y F, Liang W, Nashed Z, et al. 2014. Seismic modeling by optimizing regularized staggered-grid finite-difference operators using a time-space-domain dispersion-relationship-preserving method.Geophysics, 79(5): T277-T285.

Wu G C, Wang H Z. 2005. Analysis of numerical dispersion in wave-field simulation.ProgressinGeophysics. (in Chinese), 20(1), 58-65.Xie G S, Liu H, Zhao L G. 2005. Parallel Algorithm based on the multithread technique for pseudospectal modeling of seismic wave.ProgressinGeophysics. (in Chinese), 20(1), 17-23.

Yang Q J, Liu C, Geng M X, et al. 2014. Staggered-grid finite-difference scheme and coefficients deduction of any number of derivatives.J.JilinUniv. (Earthscienceedition), 44(1): 375-385.

Yee K S. 1966. Numerical solution of initial boundary value problems involving Maxwell′s equations in isotropic media.IEEETrans.AntennasPropagat, 14(3): 302-307.

Zhang J, Wu G C. 2015, Seismic modeling method of two phase medium utilizing staggered grid with optimal difference coefficients.ProgressinGeophysics. (in Chinese), 30(5), 2301-2311.

附中文參考文獻

董良國，馬在田，曹景中等. 2000. 一階彈性波方程交錯網格高階差分解法. 地球物理學報. 43(6): 856-864.

杜啟振，郭成鋒，公緒飛. 2015. VTI介質純P波混合法正演模擬及穩(wěn)定性分析. 地球物理學報. 58(4): 1290-1304.

梁文全，王彥飛，楊長春. 2015. 基于優(yōu)化方法的時間空間域隱格式有限差分算子確定方法. 石油物探. 54(3): 254-259.

謝桂生，劉洪，趙連功. 2005. 偽譜法地震波正演模擬的多線程并行計算. 地球物理學進展. 20(1): 17-23.

劉洋，李承楚，牟永光. 1998. 任意偶數階精度有限差分法數值模擬. 石油地球物理勘探. 33(1): 1-10.

吳國忱，王華忠. 2005. 波場模擬中的數值頻散分析與校正策略. 地球物理學進展. 20(1): 58-65.

楊慶節(jié)，劉財，耿美霞等. 2014. 交錯網格任意階導數有限差分格式及差分系數推導. 吉林大學學報(地球科學版). 44(1): 375-385.

張杰，吳國忱. 2015. 基于優(yōu)化差分系數的雙相介質交錯網格正演模擬. 地球物理學進展. 30(5): 2301-2311.

(本文編輯劉少華)

Time-space domain high-order finite-difference methods for seismic wave numerical simulation based on new stencils

ZHANG Bao-Qing1，2，ZHOU Hui1*， CHEN Han-Ming1， SHENG Shan-Bo3

1StateKeyLaboratoryofPetroleumResourcesandProspecting,CNPCKeyLabofGeophysicalExploration,ChinaUniversityofPetroleum,Beijing102249,China2DagangBranchofGRI,BGPInc.,CNPC,TianjinDagang300280,China3ChinaNationalOil&GasExplorationandDevelopmentCorporation,Beijing100034，China

AbstractThe traditional finite-difference (FD) seismic wave simulation scheme adopts high-order FD operators to discretize the spatial derivatives, and the second-order FD operator to discretize the temporal derivative. Therefore, the traditional high-order FD method only achieves high-order accuracy in space, but exhibits low-order accuracy in time. When a relatively large time step is applied, the traditional FD method suffers from visible temporal dispersion and even instability. This paper develops new time-space domain high-order FD methods that attain arbitrary even-order accuracy in space, fourth- and sixth-order accuracy in time. The new FD methods are developed based on new FD stencils and a centered-grid. The FD coefficients are determined from the discrete dispersion relation using the Taylor-series expansion (TE) approach. Dispersion analysis indicates that our temporal fourth- and sixth-order FD methods improve the accuracy of the traditional temporal second-order FD method significantly. Computational cost analysis demonstrates that our temporal high-order FD methods are more efficient than the traditional temporal second-order method. Numerical simulation of seismic waves in homogeneous and heterogeneous media further validates the effectiveness of our high-order FD methods.KeywordsSeismic wave equation; Time-space domain; High-order finite-difference; Numerical simulation

基金項目國家重點基礎研究發(fā)展計劃(973計劃)(2013CB228603)，國家自然科學基金(41174119)和中石油物探新方法新技術研究(2014A-3609)聯合資助.

作者簡介張保慶，男，1973年生，高級工程師，主要從事地震資料處理等工作. E-mail:zhangbaoqing1640@sina.com *通訊作者周輝，教授，博士生導師，主要從事地震勘探和地質雷達探測等地球物理資料的處理和正反演研究. E-mail: huizhou@cup.edu.cn

doi:10.6038/cjg20160523 中圖分類號P631

收稿日期2015-08-04，2016-04-22收修定稿

張保慶， 周輝， 陳漢明等. 2016. 基于新的差分結構的時-空域高階有限差分波動方程數值模擬方法. 地球物理學報，59(5):1804-1814,doi:10.6038/cjg20160523.

Zhang B Q, Zhou H, Chen H M, et al. 2016. Time-space domain high-order finite-difference methods for seismic wave numerical simulation based on new stencils.ChineseJ.Geophys. (in Chinese),59(5):1804-1814,doi:10.6038/cjg20160523.