基于時間高精度-空間隱式矩形交錯網格有限差分的彈性波數值模擬

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

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

0 引言

有限差分法具備簡單易行、計算效率和靈活性高等優點，因此在地震波場數值模擬中廣受青睞(Virieux，1986；Carcione et al.，2002；Etgen and O′Brien，2007；馮英杰等，2007；劉洋，2014；劉立彬等，2020；姜占東等，2021).發展至今，為了提高有限差分在離散精度、穩定性和效率等方面的表現，多種改進版本的差分方法得到發展，例如：規則網格和交錯網格差分(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)、泰勒展開和優化類差分方法(杜啟振等, 2010；Chu and Stoffa，2012；Liu，2014；王之洋等，2015；印興耀等，2015；李振春等，2016)等.總體而言，上述方法中多數可以通過增加差分算子長度來獲得空間高階模擬精度，但時間模擬精度仍為二階.

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

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

1 理論與方法

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

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

(1)

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

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

圖1 交錯網格差分模板沿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波庫朗數，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).根據方程(5)和(6)，可以采用泰勒級數展開法求取高階差分系數.基于前人的推導思路(Chen et al.，2017；Ren and Li，2017，2019)，此處直接給出N=2和3時關于P波的泰勒展開差分系數.

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

(7a)

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

(7b)

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

(8)

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

fp-ξ(khξ,θ)，

(9)

其中：

(10)

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

(11)

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

(12)

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

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

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

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

(13)

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

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

(14)

利用泰勒級數展開式(14)中的三角函數，高階差分系數可以通過對比等式兩端多項式系數獲得.具體地，高階隱式差分系數可以通過式(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)

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

(16)

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

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

(17)

其中：

(18)

2 數值精度分析

2.1 頻散分析

本節采用相速度相對誤差來分析不同差分方法的數值精度.為了簡單起見，在表1中對本文涉及的差分方法進行了簡寫.

表1 不同交錯網格差分方法的簡寫

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

(19a)

其中：

(19b)

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

(20)

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

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

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

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

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

2.2 穩定性分析

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

(21)

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

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

圖5 不同網格比和不同長度差分算子M所計算的穩定性曲線

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

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

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

(23)

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

(24)

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

(25)

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

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

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

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

在計算過程中，通過結合P波和S波單獨的波動方程與改進差分方案能夠實現彈性波場的高精度外推.

4 數值模擬算例

4.1 均勻模型

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

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

4.2 復雜模型

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

圖7 非標準Marmousi速度模型

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

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

5 結論

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