時(shí)間-空間高階精度矩形交錯(cuò)網(wǎng)格隱式有限差分聲波正演模擬

王靜， 劉洋,3*， 周泓宇

1 中國(guó)石油大學(xué)(北京)油氣資源與探測(cè)國(guó)家重點(diǎn)實(shí)驗(yàn)室， 北京 102249 2 中國(guó)石油大學(xué)(北京)CNPC物探重點(diǎn)實(shí)驗(yàn)室， 北京 102249 3 中國(guó)石油大學(xué)(北京)克拉瑪依校區(qū)， 克拉瑪依 834000

0 引言

有限差分方法(Kindelan et al.，1990；Robertsson et al.，1994；劉洋等，1998；Carcione et al.，2002；Moczo et al.，2002，2014；Du et al.，2009；Di Bartolo et al.，2012；Song et al.，2013；Fang et al.，2014；Etemadsaeed et al.，2016；Koene et al，2020；胡自多等，2021；趙明哲等，2022)因其簡(jiǎn)單易行、計(jì)算效率高和內(nèi)存小的特點(diǎn)而被廣泛應(yīng)用于地震波正演、成像及反演中.有限差分方法的本質(zhì)是采用多個(gè)臨近點(diǎn)的加權(quán)求和來近似中心點(diǎn)的時(shí)間和空間導(dǎo)數(shù).通常情況下，用離散差分算子數(shù)值逼近波動(dòng)方程中的連續(xù)微分算子會(huì)引起時(shí)間及空間頻散.兩個(gè)重要因素將會(huì)影響數(shù)值頻散的大小：差分格式和差分系數(shù)的求取方法，前者決定了用于近似中心差分點(diǎn)的時(shí)間及空間導(dǎo)數(shù)的網(wǎng)格點(diǎn)分布情況；后者給出了不同網(wǎng)格點(diǎn)的權(quán)重.規(guī)則網(wǎng)格(Dablain，1986)和交錯(cuò)網(wǎng)格(Yee，1966；Virieux，1984，1986；Graves，1996；董良國(guó)等，2000；Moczo et al.，2000；裴正林，2004；李振春等，2007；Dong et al.，2013；姜占東等，2021)是有限差分方法中最常用的兩種用于求解一階偏導(dǎo)數(shù)的差分格式.與規(guī)則網(wǎng)格相比，交錯(cuò)網(wǎng)格因其具有更高的精度和更好的穩(wěn)定性而被廣泛應(yīng)用于地震波正演模擬中.

有限差分模擬的空間精度可以通過增加差分算子長(zhǎng)度來提高，如高階顯式和隱式(Lele，1992；Du et al.，2009；Liu and Sen，2009a；Chu and Stoffa，2012a，c；Liu，2014；陳東等，2016；Yang et al.，2017；Wang and Liu，2018；Wang et al.，2018，2021)交錯(cuò)網(wǎng)格有限差分方法，前者通過增加空間算子長(zhǎng)度很容易實(shí)現(xiàn)空間高階精度模擬，但存在“飽和效應(yīng)”(Kosloff et al.，2010)的問題，即空間精度的提高率隨著算子長(zhǎng)度的增加而降低，并且增加了大量的計(jì)算消耗.隱式有限差分方法的發(fā)展很好地解決了這一問題，隱式差分方法可以使用較短的算子長(zhǎng)度達(dá)到與顯式方法相同的精度，且有效降低了計(jì)算耗時(shí).然而上述顯式和隱式有限差分方法的時(shí)間導(dǎo)數(shù)通常由二階中心差分格式進(jìn)行離散，導(dǎo)致在長(zhǎng)時(shí)間的波場(chǎng)傳播過程中，時(shí)間頻散變得嚴(yán)重，特別是當(dāng)使用大的時(shí)間步長(zhǎng)時(shí).因此，更多學(xué)者致力于發(fā)展時(shí)間高精度波場(chǎng)模擬方法(Lax and Wendroff，1964；Dablain，1986；Tam and Webb，1993；Etgen and Brandsberg-Dahl，2009；Fowler et al.，2010；Alkhalifah，2013；梁文全等，2013；Tan and Huang，2014；Ren et al.，2017；Koene et al.，2018；Ren and Li，2019；Zhou et al.，2021).

有限差分模擬的時(shí)間精度可以通過采用空間導(dǎo)數(shù)代替時(shí)間導(dǎo)數(shù)的Lax-Wendroff方法(Lax and Wendroff，1964)來實(shí)現(xiàn)，然而精確的時(shí)間高階算法依賴于使用譜方法計(jì)算空間導(dǎo)數(shù)，因此計(jì)算量大.發(fā)展時(shí)間高階有限差分方法能夠有效提高計(jì)算效率.Finkelstein和Kastner(2007，2008)給出了在固定頻率下滿足時(shí)空域頻散關(guān)系的有限差分系數(shù)，實(shí)現(xiàn)了對(duì)特定頻率下時(shí)間頻散的有效壓制.Liu和Sen(2009b，2011)進(jìn)一步發(fā)展了連續(xù)頻率范圍內(nèi)時(shí)空域規(guī)則網(wǎng)格和交錯(cuò)網(wǎng)格有限差分系數(shù)的求取方法，相比于傳統(tǒng)空間域方法的時(shí)間二階精度，該方法在二維情況下可將8個(gè)波傳播方向上的時(shí)間精度提高到任意偶數(shù)階.為了進(jìn)一步提高模擬精度，Liu和Sen(2013)提出了菱形差分格式，該方法在任意方向上均可達(dá)到時(shí)間及空間高階精度，但相比于傳統(tǒng)十字形差分格式，增加了一定的計(jì)算成本.Tan和Huang(2014)之后發(fā)展了時(shí)間4階和6階精度、空間均為任意偶數(shù)階精度的兩種時(shí)空域交錯(cuò)網(wǎng)格有限差分格式，Chen等(2016)將其推廣為更加適用于某些特殊地質(zhì)情況的矩形網(wǎng)格離散形式的差分方法.Ren等(2017)采用十字形和菱形結(jié)合的差分格式(Wang et al.，2016)實(shí)現(xiàn)了交錯(cuò)網(wǎng)格時(shí)間高階有限差分聲波正演模擬，緊接著類似的工作被拓展到彈性波模擬(Ren and Li，2017).

差分系數(shù)的求取方法也直接影響數(shù)值模擬結(jié)果，泰勒級(jí)數(shù)展開方法可以在低波數(shù)區(qū)域達(dá)到較高的模擬精度，但隨著波數(shù)的增加，精度會(huì)迅速下降.為了在更大的波數(shù)范圍獲得高精度模擬結(jié)果，許多學(xué)者對(duì)空間差分算子優(yōu)化和時(shí)空差分算子同時(shí)優(yōu)化策略進(jìn)行了研究.Chu和Stoffa(2012b)發(fā)現(xiàn)，通過在波數(shù)域設(shè)置一些二項(xiàng)式窗口推導(dǎo)出的優(yōu)化差分系數(shù)可以有效壓制數(shù)值頻散.Zhang和Yao(2013)利用模擬退火算法求解由最大范數(shù)構(gòu)造的目標(biāo)函數(shù)，有效提高了模擬精度和效率.Liu(2014)發(fā)展了全局優(yōu)化的顯式和隱式交錯(cuò)網(wǎng)格有限差分方法用于計(jì)算波動(dòng)方程的一階偏導(dǎo)數(shù)，該方法通過最小化選定區(qū)間內(nèi)的誤差來獲得優(yōu)化的空間差分系數(shù).對(duì)于多孔彈性介質(zhì)，Itzá等(2016)基于全局最優(yōu)策略的隱式交錯(cuò)網(wǎng)格差分格式實(shí)現(xiàn)了二維低頻介質(zhì)中傳播的地震波模擬.Ren和Li(2019)提出了基于最小二乘優(yōu)化算法的時(shí)間高階空間隱式有限差分方法，但該方法沒有依據(jù)整個(gè)離散波動(dòng)方程的頻散關(guān)系來近似時(shí)間和空間導(dǎo)數(shù)，并不是真正意義上的時(shí)空域方法.因此，時(shí)空域時(shí)間高階空間隱式交錯(cuò)網(wǎng)格有限差分方法還具有很大的研究潛力以提高地震波正演模擬精度.

現(xiàn)階段發(fā)展的時(shí)間高階有限差分方法主要為空間顯式或隱式方法，但通常為基于正方形網(wǎng)格離散形式的，將低了其對(duì)實(shí)際地質(zhì)問題的靈活適應(yīng)性.本文發(fā)展了一種基于矩形網(wǎng)格離散形式的時(shí)間高階空間隱式有限差分格式用于精確且高效地?cái)?shù)值求解一階變密度聲波方程，該差分格式在傳統(tǒng)十字形矩形交錯(cuò)網(wǎng)格隱式差分格式的基礎(chǔ)上增加了一些額外的離散點(diǎn)，達(dá)到有效壓制時(shí)間頻散的目的.基于本文新差分格式的時(shí)空域頻散關(guān)系及變量替換的思想，我們分別采用泰勒級(jí)數(shù)展開和最小二乘優(yōu)化兩種方法來計(jì)算相應(yīng)的有限差分系數(shù).通過頻散和穩(wěn)定性分析、數(shù)值算例來驗(yàn)證本文基于泰勒級(jí)數(shù)展開和最小二乘優(yōu)化的新差分方法相較于傳統(tǒng)時(shí)間二階空間域隱式差分方法在模擬精度、穩(wěn)定性和效率方面的優(yōu)勢(shì).

1 矩形交錯(cuò)網(wǎng)格隱式有限差分格式

二維介質(zhì)中，一階應(yīng)力-速度變密度聲波方程可表述為：

(1)

其中，p(x,z,t)為標(biāo)量壓力場(chǎng)，vx(x,z,t)和vz(x,z,t)分別代表x和z方向的質(zhì)點(diǎn)偏振速度.K(x,z)=ρ(x,z)c2(x,z)為體積模量，ρ(x,z)和c(x,z)分別為介質(zhì)密度和地震波傳播速度.

1.1 傳統(tǒng)時(shí)間二階隱式有限差分格式

(2a)

(2b)

(2c)

其中：

(3a)

(3b)

(4a)

(4b)

(4c)

(5)

據(jù)Liu和Sen(2009a)的分析可知，上述基于泰勒級(jí)數(shù)展開的空間域隱式差分格式可以達(dá)到2M+2階空間精度，具有與傳統(tǒng)4M階顯式方法(當(dāng)bα=0時(shí))基本相同的精度.相較于顯式有限差分方法，隱式方法雖然能夠顯著提高空間模擬精度，但時(shí)間精度依舊只能夠達(dá)到二階.為方便敘述，我們將基于泰勒級(jí)數(shù)展開求取差分系數(shù)的空間任意偶數(shù)階、時(shí)間二階精度的矩形交錯(cuò)網(wǎng)格隱式有限差分方法計(jì)為泰勒(2M+2，2).

1.2 新時(shí)間高階空間隱式有限差分格式

圖1 本文基于矩形交錯(cuò)網(wǎng)格離散形式的時(shí)間高階空間隱式有限差分格式離散?/?x示意圖

(6a)

(6b)

(7a)

(7b)

(7c)

該差分格式具有時(shí)間2N階精度和空間2M+2階精度，適用于任意矩形網(wǎng)格剖分的一階波動(dòng)方程的時(shí)空域隱式有限差分方法.

2 新時(shí)間高階空間隱式有限差分格式時(shí)空域差分系數(shù)求取方法

基于平面波理論和變量替換的思想，本文分別采用泰勒級(jí)數(shù)展開和線性優(yōu)化兩種方法求取相應(yīng)的差分系數(shù)并給出時(shí)間及空間精度具體階數(shù)的推導(dǎo)證明.為方便敘述，將本文具有空間2M+2階和時(shí)間2N階精度的矩形交錯(cuò)網(wǎng)格時(shí)間高階空間隱式有限差分格式記為(2M+2，2N)，則基于該格式的采用泰勒級(jí)數(shù)展開和線性優(yōu)化方法求取差分系數(shù)的兩種有限差分方法可分別記為泰勒(2M+2，2N)和優(yōu)化(2M+2，2N).

2.1 泰勒級(jí)數(shù)展開方法

假設(shè)本文(2M+2，2N)差分格式(7)具有如下形式的平面波解：

(8a)

(8b)

(8c)

(9)

其中：

(10a)

(10b)

(11a)

(11b)

以及rx=cτ/hx、rz=cτ/hz分別為x、z方向的網(wǎng)格比.

(12)

(13)

其中：

(14a)

(14b)

(15)

并將獨(dú)立波數(shù)項(xiàng)(kαhα)2l(α=x或z)和混合波數(shù)項(xiàng)(kxhx)2l-2ξ(kzhz)2ξ前面的系數(shù)進(jìn)行整合可以得到：

(16)

其中：

(17a)

(17b)

(18a)

(18b)

通過分別對(duì)比公式(16)等號(hào)兩端的獨(dú)立波數(shù)項(xiàng)(kαhα)2l和混合波數(shù)項(xiàng)(kxhx)2l-2ξ(kzhz)2ξ前面的系數(shù)，可以得到如下表達(dá)式：

(19)

(20)

(21)

將其進(jìn)一步整理為：

(22)

(23)

2.2 線性優(yōu)化方法

分解出公式(9)中的僅包含獨(dú)立波數(shù)項(xiàng)kα的表達(dá)式并對(duì)等號(hào)左側(cè)余弦函數(shù)應(yīng)用二倍角公式可得：

×sin[(u-0.5)β]，

(24)

(25)

其中：

(26)

(27)

(w=1,2,…,M+1).

(28)

為了獲得求解公式(28)所需的Bmax，我們定義優(yōu)化差分系數(shù)求取過程的最大誤差限定條件為：

(29)

其中：

(30)

3 頻散分析

本節(jié)通過對(duì)比分析本文泰勒(2M+2，2N)和優(yōu)化(2M+2，2N)有限差分方法與傳統(tǒng)泰勒(2M+2，2)有限差分方法的頻散曲線來表明本文兩種差分方法在模擬精度方面的優(yōu)勢(shì).

數(shù)值模擬速度與真實(shí)速度的比值δ可以衡量不同有限差分方法的數(shù)值頻散，基于公式(9)，可推導(dǎo)出本文(2M+2，2N)差分格式的頻散關(guān)系為：

(31)

其中，cFD和c分別為數(shù)值模擬速度和真實(shí)速度.δ>1，表現(xiàn)為時(shí)間頻散；δ<1，表現(xiàn)為空間頻散；δ=1，則無數(shù)值頻散.

圖2展示了本文泰勒(2M+2，2N)和優(yōu)化(2M+2，2N)有限差分方法的頻散曲線隨著空間算子長(zhǎng)度參數(shù)M和時(shí)間算子長(zhǎng)度參數(shù)N的變化.x和z方向的網(wǎng)格比分別為rx=0.35、rz=0.42，設(shè)置優(yōu)化(2M+2，2N)差分方法的差分系數(shù)求取時(shí)的最大允許誤差為10-4，其他相關(guān)參數(shù)已在圖中注明.由圖2a(θ=0)可見：隨著M的增大，本文兩種差分方法的頻散曲線均趨近于δ=1的零頻散參考線；與泰勒(2M+2，2N)差分方法相比，優(yōu)化(2M+2，2N)差分方法能夠在更大的波數(shù)區(qū)間內(nèi)壓制空間頻散，即取得更高的空間模擬精度.由圖2b(θ=π/4)可見：在給定相同M的情況下，隨著N的增大，本文兩種差分方法的頻散曲線形態(tài)相近(相同顏色的頻散曲線基本重合)，同時(shí)時(shí)間頻散都得到了很好的壓制.因此，本文(2M+2，2N)差分格式的空間和時(shí)間模擬精度的提高可以通過采用更大的M和N來實(shí)現(xiàn).

圖2 本文兩種矩形交錯(cuò)網(wǎng)格時(shí)間高階空間隱式有限差分方法的頻散曲線隨M和N的變化

圖3給出了本文泰勒(2M+2，2N)和優(yōu)化(2M+2，2N)差分方法與傳統(tǒng)泰勒(2M+2，2)差分方法在三種時(shí)間步長(zhǎng)的情況下沿不同傳播方向(0-π/4)的頻散曲線.設(shè)置空間算子長(zhǎng)度參數(shù)M= 6，速度c= 2200 m·s-1，x、z方向網(wǎng)格離散間距分別為22 m、18 m，其余參數(shù)可見圖上注釋.如圖3a—c所示：當(dāng)采用較小的時(shí)間步長(zhǎng)τ=1 ms(網(wǎng)格比rx=1、rz=1.23)時(shí)，三種矩形交錯(cuò)網(wǎng)格隱式有限差分方法的頻散曲線差異不大，均能獲得較為滿意的模擬結(jié)果.當(dāng)時(shí)間步長(zhǎng)τ增大到2 ms(網(wǎng)格比rx=2、rz=2.44)甚至3 ms(網(wǎng)格比rx=3、rz=3.67)時(shí)，傳統(tǒng)泰勒(2M+2，2)差分方法的頻散曲線逐漸向上遠(yuǎn)離δ=1的零頻散參考線，出現(xiàn)了嚴(yán)重的時(shí)間頻散(圖3d、g).然而，通過增大時(shí)間算子長(zhǎng)度參數(shù)N，本文的泰勒(2M+2，2N)和優(yōu)化(2M+2，2N)有限差分方法均未表現(xiàn)出明顯的時(shí)間頻散，并且優(yōu)化(2M+2，2N)差分方法在更大的波數(shù)區(qū)間內(nèi)進(jìn)一步壓制了空間頻散，是三種矩形交錯(cuò)網(wǎng)格隱式有限差分方法中模擬精度最高的.通過以上對(duì)比分析可以發(fā)現(xiàn)，當(dāng)使用大時(shí)間步長(zhǎng)進(jìn)行波場(chǎng)模擬時(shí)，傳統(tǒng)泰勒(2M+2，2)差分方法的時(shí)間頻散會(huì)顯著增加，模擬精度較低，而本文兩種基于(2M+2，2N)差分格式的時(shí)間高階空間隱式有限差分方法可通過增加時(shí)間算子長(zhǎng)度參數(shù)N來有效壓制時(shí)間頻散，因此并不會(huì)導(dǎo)致模擬精度的降低，并且優(yōu)化(2M+2，2N)差分方法又能在此基礎(chǔ)上進(jìn)一步提高空間模擬精度.本文兩種矩形交錯(cuò)網(wǎng)格時(shí)間高階空間隱式有限差分方法更適用于求解大時(shí)間步長(zhǎng)延拓時(shí)的地震波場(chǎng)，在獲得高精度模擬結(jié)果的同時(shí)能有效地節(jié)省計(jì)算時(shí)間.

圖3 三種矩形交錯(cuò)網(wǎng)格隱式有限差分方法的頻散曲線隨時(shí)間步長(zhǎng)τ的變化

4 穩(wěn)定性分析

本節(jié)采用傳統(tǒng)的特征值方法來推導(dǎo)本文(2M+2，2N)差分格式的穩(wěn)定性條件.對(duì)公式(9)進(jìn)行空間傅里葉變換，可得到以下表達(dá)式：

(32)

(33)

考慮到頻散誤差會(huì)隨波數(shù)增加而增加，因此將最大Nyquist波數(shù)kxhx=kzhz=π代入公式(33)中并進(jìn)行整理可得：

(34)

圖4展示了本文泰勒(2M+2，2N)和優(yōu)化(2M+2，2N)兩種矩形交錯(cuò)網(wǎng)格時(shí)間高階空間隱式有限差分方法隨空間和時(shí)間算子長(zhǎng)度參數(shù)M和N的穩(wěn)定性變化曲線，并與傳統(tǒng)泰勒(2M+2，2)差分方法的穩(wěn)定性結(jié)果進(jìn)行對(duì)比.相關(guān)參數(shù)設(shè)置為：rx=0.39、rz=0.44，基于最小二乘優(yōu)化方法求取差分系數(shù)時(shí)的最大允許誤差為10-4.如圖4所示：三種矩形交錯(cuò)網(wǎng)格隱式有限差分方法的穩(wěn)定性均隨著空間算子長(zhǎng)度參數(shù)M的增加而降低，其中傳統(tǒng)泰勒(2M+2，2)差分方法具有最嚴(yán)格的穩(wěn)定性條件；本文泰勒(2M+2，2N)差分方法具有最寬松的穩(wěn)定性條件.當(dāng)時(shí)間算子長(zhǎng)度參數(shù)N增加時(shí)，無論是泰勒(2M+2，2N)還是優(yōu)化(2M+2，2N)差分方法的穩(wěn)定性條件均得到明顯改善，表明時(shí)間精度的提高可以有效增強(qiáng)地震波場(chǎng)模擬的穩(wěn)定性.

圖4 三種矩形交錯(cuò)網(wǎng)格隱式有限差分方法的穩(wěn)定性曲線

5 模型算例

本節(jié)將采用簡(jiǎn)單均勻介質(zhì)模型和二維SEG/EAGE鹽丘模型來驗(yàn)證本文兩種矩形交錯(cuò)網(wǎng)格時(shí)間高階空間隱式有限差分方法相較于傳統(tǒng)空間域時(shí)間二階空間高階隱式差分方法具有更高的模擬精度和計(jì)算效率.

5.1 簡(jiǎn)單模型

第一個(gè)數(shù)值算例是模擬速度為2500 m·s-1的地震波在6 km × 6 km均勻介質(zhì)模型中傳播的地震響應(yīng).時(shí)間步長(zhǎng)為2 ms，模型縱、橫向的網(wǎng)格比分別0.33和0.25.三種矩形交錯(cuò)網(wǎng)格隱式有限差分方法的空間算子長(zhǎng)度參數(shù)M為5，其中本文兩種矩形交錯(cuò)網(wǎng)格時(shí)間高階空間隱式有限差分方法的時(shí)間算子長(zhǎng)度參數(shù)N為3.震源選擇主頻為20 Hz的雷克子波并在模型的中心點(diǎn)處激發(fā).

圖5給出了三種矩形交錯(cuò)網(wǎng)格隱式有限差分方法1.9 s時(shí)的波場(chǎng)快照，參考解(圖5a)為傳統(tǒng)泰勒(2M+2，2)差分方法取M=12，τ=0.5 ms時(shí)的模擬結(jié)果.為了更加精確地對(duì)比波場(chǎng)模擬結(jié)果，我們將每個(gè)波場(chǎng)快照中白色虛線長(zhǎng)方形框內(nèi)圖像進(jìn)行放大顯示并置于相應(yīng)波場(chǎng)圖下方.如圖5b所示：傳統(tǒng)泰勒(2M+2，2)差分方法會(huì)產(chǎn)生明顯的時(shí)間(白色箭頭標(biāo)識(shí)處)和空間(白色實(shí)線框標(biāo)識(shí)處)頻散，因橫向空間采樣間隔(20 m)大于縱向空間采樣間隔(15 m)導(dǎo)致波場(chǎng)快照橫向上的空間頻散更加嚴(yán)重.本文泰勒(2M+2，2N)差分方法的波場(chǎng)快照中(圖5c)雖然顯示出與傳統(tǒng)泰勒(2M+2，2)差分方法的波場(chǎng)快照中(圖5b)相近的空間頻散(白色實(shí)線框標(biāo)識(shí)處)，然而并無可觀察到的時(shí)間頻散，表明本文泰勒(2M+2，2N)差分方法能夠有效壓制時(shí)間頻散，提高時(shí)間模擬精度.圖5c中可觀察到的空間頻散可采用最小二乘優(yōu)化算法來進(jìn)一步壓制，獲得如圖5d所示的與參考解最為接近的波場(chǎng)快照結(jié)果.三種方法中，本文優(yōu)化(2M+2，2N)差分方法能夠有效減少時(shí)間和空間頻散，獲得最為精確的波場(chǎng)模擬結(jié)果.

圖5 不同矩形交錯(cuò)網(wǎng)格隱式有限差分方法1.9 s時(shí)的波場(chǎng)快照及其局部放大圖

5.2 二維SEG/EAGE鹽丘模型

為了進(jìn)一步驗(yàn)證本文兩種矩形交錯(cuò)網(wǎng)格時(shí)間高階空間隱式有限差分方法在模擬精度和效率方面的優(yōu)勢(shì)，下面對(duì)二維SEG/EAGE鹽丘模型(圖6)進(jìn)行一階變密度聲波方程正演模擬.該模型的速度變化范圍為1.5 km·s-1到4.5 km·s-1，密度由經(jīng)驗(yàn)公式(Castagna et al., 1993)獲得，計(jì)算區(qū)域大小為6 km×2 km且具有18 m×10 m的空間采樣間隔，時(shí)間步長(zhǎng)為1.5 ms.震源為27 Hz主頻的雷克子波并在(x，z)=(6 km，0.1 km)的位置處激發(fā).接收點(diǎn)均勻排布在地表，接收總時(shí)長(zhǎng)為4 s.圖7展示了采用三種矩形交錯(cuò)網(wǎng)格隱式有限差分方法對(duì)圖6中SEG/EAGE鹽丘模型進(jìn)行正演模擬2.3 s時(shí)的波場(chǎng)快照和參考解(傳統(tǒng)泰勒(2M+2，2)差分方法，M=15，τ=0.75 ms)，以及相應(yīng)的局部放大圖，圖8為最終接收到地震記錄的局部放大圖.由圖7和圖8可見：

圖6 二維SEG/EAGE鹽丘速度(a)及密度(b)模型

(1)傳統(tǒng)泰勒(2M+2，2)差分方法的正演模擬結(jié)果會(huì)產(chǎn)生明顯的空間和時(shí)間頻散，如圖7b和圖8b中白色實(shí)線框和白色箭頭標(biāo)識(shí)處所示.本文的泰勒(2M+2，2N)差分方法和優(yōu)化(2M+2，2N)差分方法的波場(chǎng)快照(圖7c、d)及地震記錄(圖8c、d)中基本觀察不到時(shí)間頻散，表明本文的兩種時(shí)間高階空間隱式有限差分方法能顯著提高地震波場(chǎng)正演模擬的時(shí)間精度.

圖7 SEG/EAGE鹽丘模型中不同矩形交錯(cuò)網(wǎng)格隱式有限差分方法2.3 s時(shí)的波場(chǎng)快照

(2)本文泰勒(2M+2，2N)差分方法的波場(chǎng)模擬結(jié)果(圖7c和圖8c)在有效壓制時(shí)間頻散的同時(shí)仍殘留一些空間頻散，優(yōu)化(2M+2，2N)差分方法則可以很好地解決這一問題，獲得與參考解基本相近的時(shí)間和空間高精度正演模擬結(jié)果(圖7d和圖8d).

圖8 SEG/EAGE鹽丘模型中三種矩形交錯(cuò)網(wǎng)格隱式有限差分方法計(jì)算所得地震記錄

接下來，仍然采用二維SEG/EAGE鹽丘模型來說明本文兩種矩形交錯(cuò)網(wǎng)格時(shí)間高階空間隱式有限差分方法的效率優(yōu)勢(shì).為了保證對(duì)比的公平性，我們基于統(tǒng)一計(jì)算平臺(tái)(i7-4790內(nèi)核，3.60 GHz處理器，8 GB內(nèi)存)，并保證三種差分方法之間具有相同精度的前提下比較計(jì)算耗時(shí).定義數(shù)值模擬解與參考解之間的歸一化L2范數(shù)誤差(Bohlen and Saenger，2006)來定量衡量每種差分方法的模擬精度：

(35)

其中，Iσ和Jσ分別代表數(shù)值模擬結(jié)果和參考解，L是空間采樣總點(diǎn)數(shù).

圖9為三種矩形交錯(cuò)網(wǎng)格隱式有限差分方法在(3450 m, 500 m)處的波形圖.由圖可知：當(dāng)時(shí)間步長(zhǎng)為1 ms時(shí)，傳統(tǒng)泰勒(2M+2，2)差分方法計(jì)算出的波形與參考解差異較大，表現(xiàn)出明顯的相位偏移(時(shí)間頻散)，由公式(35)計(jì)算所得誤差為0.236.然而本文兩種矩形交錯(cuò)網(wǎng)格時(shí)間高階空間隱式有限差分方法在相同時(shí)間步長(zhǎng)的情況下模擬出的波形與參考解基本匹配，其中泰勒(2M+2，2N)差分方法的誤差為0.081；優(yōu)化(2M+2，2N)差分方法能進(jìn)一步提高模擬精度，其與參考解的誤差僅為0.013.為了使傳統(tǒng)泰勒(2M+2，2)差分方法能夠達(dá)到與本文泰勒(2M+2，2N)和優(yōu)化(2M+2，2N)差分方法基本相同的精度，我們縮短用于傳統(tǒng)泰勒(2M+2，2)差分方法正演模擬的時(shí)間步長(zhǎng)至原來的一半(τ=0.5 ms).據(jù)此，傳統(tǒng)泰勒(2M+2，2)差分方法與參考解匹配良好，誤差降為0.089.由此可知：本文兩種(2M+2，2N)差分方法即便采用了相較于傳統(tǒng)泰勒(2M+2，2)差分方法兩倍的時(shí)間步長(zhǎng)，其模擬精度仍然高于傳統(tǒng)泰勒(2M+2，2)方法.在保證三種矩形交錯(cuò)網(wǎng)格隱式有限差分方法基本達(dá)到了相同精度的前提下，表1列出了SEG/EAGE鹽丘模型(圖6)中三種差分方法正演模擬所需的計(jì)算時(shí)間.本文泰勒(2M+2，2N)和優(yōu)化(2M+2，2N)差分方法的計(jì)算耗時(shí)分別為662.57 s和653.26 s，相較于傳統(tǒng)泰勒(2M+2，2)差分方法879.91 s的計(jì)算時(shí)間，分別節(jié)省了24.7%和25.8%.值得注意的是，盡管傳統(tǒng)泰勒(2M+2，2)差分方法以增加計(jì)算成本為代價(jià)(時(shí)間步長(zhǎng)減半)，其模擬精度仍然低于本文兩種(2M+2，2N)差分方法.因此，本文(2M+2，2N)有限差分格式在保證高精度正演模擬的同時(shí)可以采用更大的時(shí)間步長(zhǎng)，從而提高計(jì)算效率.

圖9 SEG/EAGE鹽丘模型中三種矩形交錯(cuò)網(wǎng)格隱式有限差分方法計(jì)算所得波形圖

表1 SEG/EAGE鹽丘模型中三種矩形交錯(cuò)網(wǎng)格隱式有限差分方法地震波正演模擬計(jì)算時(shí)間(CPU)

6 結(jié)論

本文發(fā)展的適用于矩形網(wǎng)格剖分形式的時(shí)間高階空間隱式有限差分格式可以顯著提高一階變密度聲波方程的模擬精度.基于變量替換的思想，使得本文時(shí)間高階空間隱式有限差分格式的時(shí)空域頻散關(guān)系與差分系數(shù)之間的非線性關(guān)系線性化，從而降低了差分系數(shù)的求解難度.相應(yīng)的有限差分系數(shù)的求取采用了泰勒級(jí)數(shù)展開和線性優(yōu)化兩種方法，前者可以達(dá)到時(shí)間2N階、空間2M+2階模擬精度，后者通過基于泰勒級(jí)數(shù)展開和最小二乘結(jié)合的線性優(yōu)化方法，進(jìn)一步提高了本文矩形交錯(cuò)網(wǎng)格時(shí)間高階空間隱式有限差分格式的空間模擬精度.通過與傳統(tǒng)空間域矩形交錯(cuò)網(wǎng)格隱式有限差分格式的多方面對(duì)比，表明了本文差分格式在模擬精度和穩(wěn)定性方面的優(yōu)越性.此外，在相同精度的情況下，基于本文差分格式的泰勒(2M+2，2N)和優(yōu)化(2M+2，2N)差分方法可以使用比傳統(tǒng)(2M+2，2)差分方法更大的時(shí)間步長(zhǎng)，因此降低了計(jì)算耗時(shí)，更適用于大尺度模型的地震波場(chǎng)模擬.