基于時(shí)間高精度-空間隱式矩形交錯(cuò)網(wǎng)格有限差分的彈性波數(shù)值模擬

徐世剛， 包乾宗， 任志明， 劉洋

1 長安大學(xué)地質(zhì)工程與測(cè)繪學(xué)院地球物理系， 西安 710054 2 中國石油大學(xué)(北京)油氣資源與探測(cè)國家重點(diǎn)實(shí)驗(yàn)室， 北京 102249

0 引言

有限差分法具備簡(jiǎn)單易行、計(jì)算效率和靈活性高等優(yōu)點(diǎn)，因此在地震波場(chǎng)數(shù)值模擬中廣受青睞(Virieux，1986；Carcione et al.，2002；Etgen and O′Brien，2007；馮英杰等，2007；劉洋，2014；劉立彬等，2020；姜占東等，2021).發(fā)展至今，為了提高有限差分在離散精度、穩(wěn)定性和效率等方面的表現(xiàn)，多種改進(jìn)版本的差分方法得到發(fā)展，例如：規(guī)則網(wǎng)格和交錯(cuò)網(wǎng)格差分(Dablain，1986；Virieux，1986；董良國等，2000；Yan et al.，2016)、空間域和時(shí)空域差分(裴正林，2004；Finkelstein and Kastner，2007；Liu and Sen, 2009a；梁文全等，2013；Wang et al.，2014)、顯格式和隱格式差分(Lele，1992；Liu and Sen, 2009b)、泰勒展開和優(yōu)化類差分方法(杜啟振等, 2010；Chu and Stoffa，2012；Liu，2014；王之洋等，2015；印興耀等，2015；李振春等，2016)等.總體而言，上述方法中多數(shù)可以通過增加差分算子長度來獲得空間高階模擬精度，但時(shí)間模擬精度仍為二階.

在保證空間高階精度的前提下如何提高有限差分在時(shí)間域的模擬精度成為重要的研究方向.通過采用不同的空間導(dǎo)數(shù)組合來替換高階時(shí)間導(dǎo)數(shù)，Dablain(1986)和Chen(2011)推導(dǎo)了具有時(shí)間高階精度的Lax-Wendroff方法來模擬聲波傳播.與Lax-Wendroff方法相比，基于改進(jìn)模板的時(shí)間高階差分方法能夠有效提高時(shí)間模擬精度.Liu和Sen(2013)發(fā)展了基于菱形模板的規(guī)則網(wǎng)格差分法來計(jì)算二維聲波方程中的拉普拉斯算子，該方法能夠同時(shí)獲得時(shí)間和空間高階模擬精度.為了提高菱形差分模板的運(yùn)算效率，Wang等(2016)有效組合菱形算子與傳統(tǒng)十字形算子，發(fā)展了一種時(shí)間、空間差分算子長度相互獨(dú)立的時(shí)間高階差分方法.張保慶等(2016)發(fā)展了時(shí)間四階和六階，空間高階精度規(guī)則網(wǎng)格差分法求解二維聲波方程.在Liu和Sen(2013)的工作基礎(chǔ)上，Tan和Huang(2014)針對(duì)一階聲波方程發(fā)展了兩種改進(jìn)的交錯(cuò)網(wǎng)格差分方法，通過在傳統(tǒng)差分模板中引入額外的網(wǎng)格節(jié)點(diǎn)來共同近似偏導(dǎo)數(shù)，該方案能夠同時(shí)獲得空間高階，時(shí)間四階和六階模擬精度.Chen等(2017)采用波數(shù)域算子進(jìn)一步補(bǔ)償時(shí)間高階差分的精度.Ren和Li(2017)將Wang等(2016)發(fā)展的時(shí)間高階規(guī)則網(wǎng)格差分模板推廣到彈性波交錯(cuò)網(wǎng)格數(shù)值模擬中.然而，上述時(shí)間高階差分方法主要屬于顯式差分，即對(duì)波動(dòng)方程中的偏導(dǎo)數(shù)直接顯式求解.與顯格式差分相比，隱格式差分在算子長度一致時(shí)能夠獲得更高的近似精度.因此，應(yīng)用具有時(shí)間高階精度的隱式差分求解地震波動(dòng)方程逐漸受到更多研究.針對(duì)聲波方程，Wang和Liu(2018)發(fā)展了時(shí)間高階，空間隱式規(guī)則網(wǎng)格差分，有效壓制了頻散誤差.其后，Ren和Li(2019)發(fā)展了時(shí)間高階，空間隱式交錯(cuò)網(wǎng)格差分算子.

到目前為止，時(shí)間高階隱式差分方案主要針對(duì)聲波方程設(shè)計(jì)的，針對(duì)彈性波方程的較少，此外在離散過程中所采用的網(wǎng)格單元主要為正方形，靈活性較低.為了提高彈性波數(shù)值模擬的精度和靈活性，本文發(fā)展了基于矩形網(wǎng)格單元(離散網(wǎng)格單元的長度和寬度在空間不同方向不相等)的改進(jìn)時(shí)間高精度-空間隱式交錯(cuò)網(wǎng)格差分法.本文的核心思想是將改進(jìn)模板與時(shí)間二階差分算子聯(lián)合用于求解彈性波方程中的時(shí)間導(dǎo)數(shù)，同時(shí)采用隱式差分算子求解空間導(dǎo)數(shù)，提供了泰勒展開和最小二乘兩種算法求取的差分系數(shù)，將波場(chǎng)分離引入到改進(jìn)差分方案中以獲得高精度彈性波場(chǎng).數(shù)值分析和模型算例驗(yàn)證了改進(jìn)方法的有效性.

1 理論與方法

1.1 基于矩形單元的時(shí)間高精度交錯(cuò)網(wǎng)格離散格式

各向同性介質(zhì)中的二維彈性波方程可以表示為下述一階形式(Virieux，1986；董良國等，2000)：

(1)

其中，vx和vz為速度分量，τxx、τzz和τxz為應(yīng)力分量，ρ為密度，λ和μ為拉梅常數(shù).

時(shí)間二階，空間高階交錯(cuò)網(wǎng)格差分法通常被用于離散方程(1)中的時(shí)間和空間導(dǎo)數(shù)，模擬精度在空間域可以達(dá)到高階，但時(shí)間精度僅為二階.為了進(jìn)一步提高精度，在前人的研究基礎(chǔ)上(Chen et al., 2017; Ren and Li，2017，2019)，本文設(shè)計(jì)了一種改進(jìn)的交錯(cuò)網(wǎng)格差分模板，如圖1a所示.與現(xiàn)有基于正方形網(wǎng)格單元的差分模板相比(圖1b)(Ren and Li，2017, 2019)，改進(jìn)模板所采用的長方形網(wǎng)格單元更具靈活性.結(jié)合改進(jìn)差分模板和傳統(tǒng)時(shí)間二階差分算子來近似方程(1)中的時(shí)間導(dǎo)數(shù)，離散格式可以修改為：

圖1 交錯(cuò)網(wǎng)格差分模板沿x軸示意圖

(2a)

(2b)

(2c)

(2d)

(2e)

其中：

(3a)

(3b)

(3c)

(3d)

基于平面波理論，方程(1)中的波場(chǎng)分量可以表示為：

(4)

將方程(4)代入到方程(2)和(3)中并進(jìn)一步化簡(jiǎn)，可以推導(dǎo)出與P波和S波相關(guān)的兩種頻散關(guān)系式為(Ren and Li, 2017, 2019):

(5a)

(5b)

其中rp-x=vpΔt/hx，rp-z=vpΔt/hz，rs-x=vsΔt/hx，rs-z=vsΔt/hz分別為沿x軸和z軸的P波和S波庫朗數(shù)，A和B為：

×sin[(m-1/2)kxhx]cos(nkzhz)，

(6a)

×sin[(m-1/2)kzhz]cos(nkxhx).

(6b)

相關(guān)模擬實(shí)驗(yàn)表明，算子N的長度取到3能夠同時(shí)兼顧計(jì)算精度和效率(Tan and Huang，2014；Chen et al., 2017; Ren and Li, 2017, 2019).根據(jù)方程(5)和(6)，可以采用泰勒級(jí)數(shù)展開法求取高階差分系數(shù).基于前人的推導(dǎo)思路(Chen et al.，2017；Ren and Li，2017，2019)，此處直接給出N=2和3時(shí)關(guān)于P波的泰勒展開差分系數(shù).

(1)N= 2，差分系數(shù)的表達(dá)式為：

(7a)

(2)N= 3，差分系數(shù)的表達(dá)式為：

(7b)

其中ξ,ψ∈{x,z}且有ξ≠ψ，cm(m=1,2,…,M)為傳統(tǒng)空間域差分系數(shù)，表達(dá)式為：

(8)

與S波相關(guān)的差分系數(shù)可以采用類似方法進(jìn)行求解.

fp-ξ(khξ,θ)，

(9)

其中：

(10)

×dθdkhξ，(g=1,2,…,M)，

(11)

其中β為khξ的上限，可由式(12)決定：

(12)

其中κmax為預(yù)先設(shè)置的最大誤差.

1.2 基于矩形單元的空間高階隱式交錯(cuò)網(wǎng)格離散格式

與顯式差分相比，在相同的算子長度下，隱式差分能夠產(chǎn)生更高的模擬精度(Liu and Sen，2009b；Liu，2014；Ren and Li，2019).為了提高有限差分空間模擬精度，本文采用高階隱式交錯(cuò)網(wǎng)格有限差分法求解方程(1)中的空間偏導(dǎo)數(shù).

以偏導(dǎo)數(shù)?τxx/?x為例，其對(duì)應(yīng)的隱式差分格式可以表示為(Liu and Sen，2009b；Liu，2014)：

(13)

其中，q=dτxx/dx，dm和a為隱式差分系數(shù)，q可以通過求解三對(duì)角矩陣獲得(Liu and Sen，2009b)，方程(1)中的其余空間導(dǎo)數(shù)可以采用類似的差分格式求取.

基于平面波解，方程(13)可以化簡(jiǎn)為(Liu，2014)：

(14)

利用泰勒級(jí)數(shù)展開式(14)中的三角函數(shù)，高階差分系數(shù)可以通過對(duì)比等式兩端多項(xiàng)式系數(shù)獲得.具體地，高階隱式差分系數(shù)可以通過式(15)進(jìn)行求取(Liu and Sen，2009b；Liu，2014)：

a=(1-2d′M+1)/d′M+1，dm=d′m/d′M+1，

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

(15)

其中，系數(shù)d′m(m=1,2,…,M-1,M)通過求解方程(16)組獲得(Liu and Sen, 2009b)：

(16)

此外，基于最小二乘算法的高階差分系數(shù)通過式(17)求取(Liu，2014)：

(n=1,2,...,M-1,M)，

(17)

其中：

(18)

2 數(shù)值精度分析

2.1 頻散分析

本節(jié)采用相速度相對(duì)誤差來分析不同差分方法的數(shù)值精度.為了簡(jiǎn)單起見，在表1中對(duì)本文涉及的差分方法進(jìn)行了簡(jiǎn)寫.

表1 不同交錯(cuò)網(wǎng)格差分方法的簡(jiǎn)寫

基于平面波理論，可以推導(dǎo)出P波頻散關(guān)系式為：

(19a)

其中：

(19b)

基于方程(19)，定義式(20)為：

(20)

其中δp代表P波相速度的相對(duì)誤差.如果采用vs替換方程(20)中的vp，S波相速度相對(duì)誤差δs可以采用類似方法求取.當(dāng)δp和δs越接近1，頻散誤差越小.當(dāng)Ix=Iz=0時(shí)，方程(19)和(20)可以用來計(jì)算顯式差分法的頻散誤差.

采用均勻模型測(cè)試不同方法的數(shù)值頻散，圖2展示了沿三個(gè)固定傳播角度的P波頻散曲線.可以觀察到，傳統(tǒng)基于泰勒展開和最小二乘優(yōu)化的時(shí)間二階顯式和隱式差分法計(jì)算的頻散曲線偏離1，表明其存在較強(qiáng)頻散誤差.相比之下，時(shí)間四階和六階TE-THR-ESFD(12,2N)(N=2，3)的計(jì)算結(jié)果更接近于1，能夠有效壓制數(shù)值頻散.此外，隱式差分TE-和LS-THR-ISFD(12,2N)(N=2，3)具有更高模擬精度，尤其是優(yōu)化隱式差分方案.關(guān)于S波的頻散誤差函數(shù)δs，可以得出類似的模擬結(jié)果.

圖2 不同差分方法計(jì)算的P波頻散曲線傳播角度分別為(a)θ=0°，(b)θ=22.5°，(c)θ=45°. 相關(guān)參數(shù)為vp=4500 m·s-1，Δt=1 ms，hx=10 m，hz=7.5 m.

圖3和圖4分別對(duì)比了不同時(shí)間高精度差分方法計(jì)算的P波和S波頻散曲線.從圖中可以觀察到，相比于泰勒展開法，最小二乘優(yōu)化能夠獲得更高的精度和更寬的有效波數(shù)范圍.基于P波頻散關(guān)系的差分法能夠產(chǎn)生較小的P波數(shù)值頻散和嚴(yán)重的S波數(shù)值頻散.反之，基于S波頻散關(guān)系的差分法能夠產(chǎn)生明顯的P波數(shù)值頻散和輕微的S波數(shù)值頻散.這些計(jì)算結(jié)果表明在彈性波數(shù)值模擬過程中應(yīng)該采用相互獨(dú)立的P波和S波差分方案去分別實(shí)現(xiàn)P波和S波的波場(chǎng)外推.

圖3 不同差分方法計(jì)算的P波頻散曲線傳播角度分別為(a) θ=0°，(b) θ=22.5°，(c) θ=45°. 相關(guān)參數(shù)為vp=4500 m·s-1，vs=2250 m·s-1，Δt=1 ms，hx=10 m，hz=7.5 m.

圖4 不同差分方法計(jì)算的S波頻散曲線傳播角度分別為(a) θ=0°，(b) θ=22.5°，(c) θ=45°. 相關(guān)參數(shù)為vp=4500 m·s-1，vs=2250 m·s-1，Δt=1 ms，hx=10 m，hz=7.5 m.

2.2 穩(wěn)定性分析

對(duì)于給定的彈性波模擬，rp-x通常大于rs-x，因此本文僅利用P波穩(wěn)定性因子來進(jìn)行穩(wěn)定性分析.應(yīng)用特征值分析方法(Chen，2011；Liu and Sen，2013；Ren and Li，2017，2019)，本文差分方案的穩(wěn)定性條件可以表示為：

(21)

其中sp-x代表沿x方向的P波最大庫朗數(shù)：

圖5展示了不同差分方法計(jì)算的兩組穩(wěn)定性曲線圖，其中網(wǎng)格單元的長度和寬度之比不一致.可以觀察到，隨著差分算子M的增加，穩(wěn)定性因子呈遞減趨勢(shì)，表明穩(wěn)定性逐漸變差.此外，傳統(tǒng)泰勒展開和最小二乘差分方法具有最小的穩(wěn)定性因子，表明這些方法容易遭受不穩(wěn)定.相比于顯式方案，隱式方案會(huì)降低穩(wěn)定性.隨著參數(shù)N的增加，穩(wěn)定性呈遞增趨勢(shì).

圖5 不同網(wǎng)格比和不同長度差分算子M所計(jì)算的穩(wěn)定性曲線

3 結(jié)合時(shí)間高精度隱式差分方案和波場(chǎng)分離的彈性波數(shù)值模擬

如前文所述，P波差分系數(shù)僅能壓制P波頻散，S波差分系數(shù)也僅能獲得高精度的S波模擬結(jié)果.對(duì)于轉(zhuǎn)換波，上述差分方案不能有效壓制其頻散誤差.因此，為了獲得高精度的P波、S波和轉(zhuǎn)換波波場(chǎng)，本文將結(jié)合時(shí)間高精度隱式差分模板和波場(chǎng)分離技術(shù)進(jìn)行彈性波數(shù)值模擬(Chen et al.，2017；Ren and Li，2017).間接分離算法通過在時(shí)間-空間域求解相互獨(dú)立的P波和S波方程來實(shí)現(xiàn)彈性波場(chǎng)分離，該方法簡(jiǎn)單易行，且具有較高計(jì)算效率(馬德堂和朱光明，2003；李振春等，2007；Xiao and Leaney，2010；Wang et al.，2015；Chen et al.，2017；Ren and Li，2017)，因此本文將該方法直接推廣到改進(jìn)差分方案中.

方程(1)中的全波場(chǎng)vx和vz可以表示為下述波場(chǎng)分量之和(馬德堂和朱光明，2003；李振春等，2007)

(23)

采用統(tǒng)一的正應(yīng)力τp替換方程(1)中的τxx和τzz，P波分量可以表示為：

(24)

相應(yīng)的，S波分量可以將方程(23)和(24)代入方程(1)中整理化簡(jiǎn)得到：

(25)

基于改進(jìn)差分方法和波場(chǎng)分離技術(shù)的彈性波數(shù)值模擬流程可以劃分為以下步驟(Ren and Li，2017)：

(1)利用P波和S波差分系數(shù)分別計(jì)算方程(24)和(25)，實(shí)現(xiàn)P波和S波的波場(chǎng)外推.

(2)基于方程(23)獲得每個(gè)迭代時(shí)間內(nèi)的全波場(chǎng)vx和vz.

(3)重復(fù)上述兩個(gè)步驟直到迭代時(shí)間中止.

在計(jì)算過程中，通過結(jié)合P波和S波單獨(dú)的波動(dòng)方程與改進(jìn)差分方案能夠?qū)崿F(xiàn)彈性波場(chǎng)的高精度外推.

4 數(shù)值模擬算例

4.1 均勻模型

第一個(gè)算例采用均勻模型，模型尺寸為7500 m×6000 m，采用矩形網(wǎng)格單元對(duì)其進(jìn)行剖分，網(wǎng)格大小為12.5 m×10 m.P波和S波的速度分別為3400 m·s-1和1900 m·s-1，密度為2300 kg·m-3.震源選用主頻30 Hz的雷克子波，置于模型中央產(chǎn)生振動(dòng).時(shí)間采樣間隔為1.5 ms.圖6對(duì)不同方法計(jì)算的波場(chǎng)快照進(jìn)行了局部展示，其中圖6a為優(yōu)化時(shí)間二階顯式差分采用較長算子長度和較小時(shí)間步長計(jì)算的參考解.對(duì)比不同波場(chǎng)切片，可以觀察到，與參考解相比，泰勒展開時(shí)間二階顯式差分(TE-ESFD [16,2])出現(xiàn)較強(qiáng)數(shù)值頻散，如圖6b所示.相比之下，圖6c中的優(yōu)化時(shí)間二階差分(LS-ISFD [8,2])，圖6d、e中的泰勒展開時(shí)間高精度差分方法(TE-THR-ESFD [16,4]和TE-THR-ISFD [8,4])僅能在一定程度上抑制數(shù)值頻散.與上述幾種差分方法相比，改進(jìn)的LS-THR-ISFD (8,4)能夠有效壓制數(shù)值頻散，產(chǎn)生高精度的模擬結(jié)果，如圖6f所示.此外，相比于其他方法，改進(jìn)的優(yōu)化隱式差分方案可以采用較小的算子長度獲得相似精度的模擬結(jié)果，能夠有效保證計(jì)算效率.

圖6 均勻模型中不同差分方法計(jì)算的0.9 s時(shí)刻波場(chǎng)快照局部顯示

4.2 復(fù)雜模型

第二個(gè)算例采用非標(biāo)準(zhǔn)的Marmousi模型，圖7僅展示了P波速度模型，S波和P波速度之比設(shè)置為0.5，密度為常數(shù)2000 kg·m-3.將計(jì)算區(qū)域剖分為501×353個(gè)網(wǎng)格單元，每個(gè)網(wǎng)格單元的尺寸為15 m×10 m.時(shí)間采樣間隔為1.2 ms.主頻15 Hz的雷克子波置于點(diǎn)(3750 m, 10 m)處產(chǎn)生P波震源.接收排列的深度為150 m.采用10層厚度的一階彈性波混合吸收邊界條件壓制人工截?cái)噙吔绶瓷?任志明和劉洋，2014).圖8為時(shí)間二階優(yōu)化顯式差分法計(jì)算的參考地震記錄圖.對(duì)圖8進(jìn)行局部放大，并將不同差分方法計(jì)算的局部地震記錄在圖9中進(jìn)行顯示.從圖中可以觀察到，相比于參考解，泰勒展開類差分方法(TE-ESFD [16, 2]，TE-THR-ESFD [16,4]和TE-THR-ISFD [8,4])存在較強(qiáng)數(shù)值頻散干擾，如圖9b—d中箭頭所示.圖9e、f中的時(shí)間二階優(yōu)化差分和時(shí)間高精度優(yōu)化差分(LS-ESFD [16,2]和LS-THR-ISFD [8,4])能夠顯著壓制數(shù)值頻散，提高模擬精度.此外，與采用較長算子長度的顯式差分相比，隱式差分可以采用較小的算子長度獲得與之相近的模擬精度，可以減少浮點(diǎn)數(shù)運(yùn)算，有效提高計(jì)算效率，尤其是采用優(yōu)化差分系數(shù)的隱式差分方案.

圖7 非標(biāo)準(zhǔn)Marmousi速度模型

圖8 非標(biāo)準(zhǔn)Marmousi速度模型參考地震記錄圖

圖9 非標(biāo)準(zhǔn)Marmousi速度模型中不同差分方法計(jì)算的地震記錄圖局部顯示

5 結(jié)論

本文通過結(jié)合十字形和菱形差分模板，發(fā)展了一種基于矩形網(wǎng)格單元的時(shí)間高精度，空間隱式交錯(cuò)網(wǎng)格差分模板來模擬彈性波傳播.改進(jìn)方案結(jié)合時(shí)間二階離散格式和改進(jìn)差分模板來求解彈性波方程中的時(shí)間導(dǎo)數(shù)，采用隱式差分算子近似空間導(dǎo)數(shù).應(yīng)用泰勒級(jí)數(shù)展開和最小二乘優(yōu)化兩種方法計(jì)算了高階差分系數(shù).為了產(chǎn)生高精度的P波、S波和轉(zhuǎn)換波波場(chǎng)，本文采用與P波和S波相關(guān)的差分系數(shù)來分別求解獨(dú)立的P波和S波方程.數(shù)值分析和模型算例證明本文隱式差分方案能夠有效壓制數(shù)值頻散干擾，產(chǎn)生高精度的模擬結(jié)果.