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

1 長安大學地質(zhì)工程與測繪學院地球物理系， 西安 710054 2 中國石油大學(北京)油氣資源與探測國家重點實驗室， 北京 102249

0 引言

有限差分法具備簡單易行、計算效率和靈活性高等優(yōu)點，因此在地震波場數(shù)值模擬中廣受青睞(Virieux，1986；Carcione et al.，2002；Etgen and O′Brien，2007；馮英杰等，2007；劉洋，2014；劉立彬等，2020；姜占東等，2021).發(fā)展至今，為了提高有限差分在離散精度、穩(wěn)定性和效率等方面的表現(xiàn)，多種改進版本的差分方法得到發(fā)展，例如：規(guī)則網(wǎng)格和交錯網(wǎng)格差分(Dablain，1986；Virieux，1986；董良國等，2000；Yan et al.，2016)、空間域和時空域差分(裴正林，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ù)，Dablain(1986)和Chen(2011)推導了具有時間高階精度的Lax-Wendroff方法來模擬聲波傳播.與Lax-Wendroff方法相比，基于改進模板的時間高階差分方法能夠有效提高時間模擬精度.Liu和Sen(2013)發(fā)展了基于菱形模板的規(guī)則網(wǎng)格差分法來計算二維聲波方程中的拉普拉斯算子，該方法能夠同時獲得時間和空間高階模擬精度.為了提高菱形差分模板的運算效率，Wang等(2016)有效組合菱形算子與傳統(tǒng)十字形算子，發(fā)展了一種時間、空間差分算子長度相互獨立的時間高階差分方法.張保慶等(2016)發(fā)展了時間四階和六階，空間高階精度規(guī)則網(wǎng)格差分法求解二維聲波方程.在Liu和Sen(2013)的工作基礎上，Tan和Huang(2014)針對一階聲波方程發(fā)展了兩種改進的交錯網(wǎng)格差分方法，通過在傳統(tǒng)差分模板中引入額外的網(wǎng)格節(jié)點來共同近似偏導數(shù)，該方案能夠同時獲得空間高階，時間四階和六階模擬精度.Chen等(2017)采用波數(shù)域算子進一步補償時間高階差分的精度.Ren和Li(2017)將Wang等(2016)發(fā)展的時間高階規(guī)則網(wǎng)格差分模板推廣到彈性波交錯網(wǎng)格數(shù)值模擬中.然而，上述時間高階差分方法主要屬于顯式差分，即對波動方程中的偏導數(shù)直接顯式求解.與顯格式差分相比，隱格式差分在算子長度一致時能夠獲得更高的近似精度.因此，應用具有時間高階精度的隱式差分求解地震波動方程逐漸受到更多研究.針對聲波方程，Wang和Liu(2018)發(fā)展了時間高階，空間隱式規(guī)則網(wǎng)格差分，有效壓制了頻散誤差.其后，Ren和Li(2019)發(fā)展了時間高階，空間隱式交錯網(wǎng)格差分算子.

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

1 理論與方法

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

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

(1)

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

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

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

(2a)

(2b)

(2c)

(2d)

(2e)

其中：

(3a)

(3b)

(3c)

(3d)

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

(4)

將方程(4)代入到方程(2)和(3)中并進一步化簡，可以推導出與P波和S波相關的兩種頻散關系式為(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)

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

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

(7a)

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

(7b)

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

(8)

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

fp-ξ(khξ,θ)，

(9)

其中：

(10)

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

(11)

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

(12)

其中κmax為預先設置的最大誤差.

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

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

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

(13)

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

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

(14)

利用泰勒級數(shù)展開式(14)中的三角函數(shù)，高階差分系數(shù)可以通過對比等式兩端多項式系數(shù)獲得.具體地，高階隱式差分系數(shù)可以通過式(15)進行求取(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é)采用相速度相對誤差來分析不同差分方法的數(shù)值精度.為了簡單起見，在表1中對本文涉及的差分方法進行了簡寫.

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

基于平面波理論，可以推導出P波頻散關系式為：

(19a)

其中：

(19b)

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

(20)

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

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

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

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

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

圖4 不同差分方法計算的S波頻散曲線傳播角度分別為(a) θ=0°，(b) θ=22.5°，(c) θ=45°. 相關參數(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)定性分析

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

(21)

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

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

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

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

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

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

(23)

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

(24)

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

(25)

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

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

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

(3)重復上述兩個步驟直到迭代時間中止.

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

4 數(shù)值模擬算例

4.1 均勻模型

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

圖6 均勻模型中不同差分方法計算的0.9 s時刻波場快照局部顯示

4.2 復雜模型

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

圖7 非標準Marmousi速度模型

圖8 非標準Marmousi速度模型參考地震記錄圖

圖9 非標準Marmousi速度模型中不同差分方法計算的地震記錄圖局部顯示

5 結(jié)論

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