求解彈性波動方程的頻率域近似解析離散化波場模擬方法

郎超， 仇楚鈞， 劉少林， 申文豪， 李小凡， 徐錫偉*

1 北京信息科技大學理學院， 北京 100192 2 清華大學數學科學系， 北京 100084 3 應急管理部國家自然災害防治研究院， 北京 100085 4 中國地質大學(武漢)地球物理與空間信息學院， 武漢 430074

0 引言

隨著計算機硬件水平和地震觀測技術的不斷發展，越來越多的研究者開始關注基于波動方程的高分辨率地震成像方法(Tape et al., 2009; Lei et al., 2020; Liu et al., 2018a)，例如雙程波逆時偏移(劉洪等，2009；Yan et al., 2015)和波動方程伴隨層析成像(Tarantola, 1984; Pratt and Worthington, 1990; Tromp et al., 2004).相比于傳統射線走時地震成像技術，一方面，基于波動方程的成像方法可同時利用地震觀測數據中的走時信息和波形信息，這為地下結構反演提供更多的有力約束(Tarantola, 1986; Tape et al., 2007)；另一方面，波動方程成像方法以聲波或者彈性波方程作為數學模型，能夠更加貼切地模擬實際地震波在地下傳播的動力學特征(Liu et al., 2014b).

基于波動方程的地震成像既可在“時間-空間”域進行，也可在“頻率-空間”域開展(Pratt et al., 1998; Chen and Cao, 2018).相比于時間域情形，頻率域計算主要具有以下四方面優勢：(1)由于頻率域波動方程只涉及空間偏導數項并且關于各個頻率點是相互解耦的，所以頻率域計算具有天然的并行特性；(2)可以避免數值離散誤差的累積(Pratt, 1999)；(3)通過引入復速度或者復頻率(Song and Williamson, 1995)，能夠更方便地刻畫衰減效應(Wang and Rao, 2009)；(4)頻率域反演成像只需選擇少量頻率點計算(Sirgue and Pratt, 2004)，消耗的計算代價較小.因此，研究頻率域全波形成像方法具有重要的科學意義與實際應用價值(Chen and Cao, 2016).

類似于時間域情形，頻率域成像過程的主要計算量也集中于對地震波傳播規律的正演模擬.因此，正演算法的優劣直接決定著全波形成像結果的精度和整體計算效率.目前用于地震波場模擬的主流正演算法包括有限差分方法(Alford et al.,1974; Virieux, 1984; Moczo et al., 2007)、有限元方法(Marfurt, 1984; Liu et al., 2014a； Liu et al., 2017c)、偽譜法(Fornberg, 1975; Kosloff and Baysal, 1982; 黃繼偉和劉洪，2020)、譜元法(Komatitsch et al., 2005; Liu et al.,2017a)等.有限差分方法因其原理直觀、格式構造簡單、實現容易、便于并行等優點，已成為地震波模擬和成像的有力工具之一(Yang et al., 2006; 劉少林等，2013).然而傳統有限差分方法容易引起數值頻散效應，當遇到復雜地質結構(例如強間斷面)時，往往需要加密離散網格以提升數值精度，這將導致計算效率降低(Moczo et al., 2000; Liu et al., 2017c).近似解析離散化(Nearly Analytic Discrete，NAD)方法作為一種新型差分方法，由楊頂輝等(Yang et al., 2003)首先引入到時間域地震波場數值計算中.該方法的基本構造思想是同時采用波場值以及波場梯度值來近似高階偏導數項(Yang et al., 2004)，這不僅繼承了普通有限差分方法的優點，同時也能在粗網格下有效地壓制數值頻散，進而得到準確的數值結果，因此該方法可以提高地震波數值模擬的計算效率(Tong et al., 2013; Liu et al., 2015).

求解離散頻率域波動方程形成的大型稀疏線性代數方程組是頻率域全波形成像的核心問題之一.該方程組的系數矩陣(也叫做“阻抗矩陣”(Pratt, 1999))通常情況下具有不定性，并且在某些常用地震波頻率取值下病態程度嚴重，這為有效求解相應線性方程組帶來了挑戰.目前主流的求解算法有兩種：直接法和迭代法(Ernst and Gander, 2012).直接法是基于對阻抗矩陣進行帶有排序(George, 1977)的三角分解，它的優點是對于多震源情形，分解過程只需進行一次，可以減少計算代價.然而直接法通常只適用于二維小尺度區域地震波模擬與波形成像計算.這是由于該方法不能充分利用阻抗矩陣的稀疏分塊結構，導致分解矩陣產生大量非零元素，引起內存和計算量的急劇增加，計算效率偏低.對于大尺度區域模型頻率域地震波場模擬，迭代法是較為合適的選擇.這類方法主要包括基于代數預處理子的Krylov子空間方法(Saad, 2003)、多重網格(Multigrid，MG)方法(Plessix, 2007)以及區域分解方法(Domain Decomposition Methods，DDM)(Gander et al., 2007)等.迭代法在計算過程中只涉及若干次稀疏矩陣與向量乘積，采用稀疏處理技術可以避免內存和計算操作數的大量消耗(Bai, 2015)，從而更加適用于大規模計算.

為了降低阻抗矩陣的病態程度從而達到快速穩定求解線性方程組的目的，Lang和Ren(2015)提出采用不精確旋轉分塊三角預處理子(Inexact Rotated Block Triangular Preconditioners, IRBTP)結合廣義極小殘量(Generalized Minimum Residues, GMRES)迭代方法(Saad and Schultz, 1986).其中“不精確”是指采用不精確預處理的思想，在保證準確性的前提下可充分利用相關矩陣的稀疏特性，節約計算成本.近期的研究成果(Lang et al., 2020)進一步完善了IRBTP預處理矩陣的特征值性質，可以從理論上保證相應預處理迭代方法的收斂特性.通過采用數值試驗與其他經典迭代方法進行比較，IRBTP預處理迭代方法在加速頻率域聲波波場模擬方面的優勢得到了充分體現(Lang and Yang, 2017).然而IRBTP預處理方法對于頻率域彈性波場模擬的數值效率尚不明確，因此有必要進一步研究.

由于彈性波方程比聲波方程更能準確地刻畫地震波在地下介質中的傳播規律，為了貼近實際情形，通常采用彈性波動方程作為模擬地震波傳播的數學模型.然而彈性波方程的結構較為復雜，相應求解過程比聲波方程消耗的計算代價更多，因此構造高效的頻率域彈性波動方程數值模擬方法至關重要.本文主要介紹頻率域彈性波方程四階NAD方法的詳細離散過程，同時也分別采用四階普通有限差分(Ordinary Finite Difference，OFD)方法(Charl-Hyun et al., 1996)和四階交錯網格(Staggered Grid，SG)(Moczo et al., 2000)方法進行數值離散，在檢驗IRBTP預處理迭代方法的數值效率之后將其運用于相應的彈性波場計算中.分別通過對各種介質模型進行波場模擬、數值頻散分析以及與解析解的波形對比，驗證NAD方法較其他兩種典型數值方法在提高波場求解準確度和計算效率方面的優勢.

1 頻率域彈性波方程的NAD離散過程

考慮各向同性二維頻率域彈性波動方程：

(1)

其中u,v分別表示彈性波場的水平和垂直分量，ρ表示介質密度，λ和μ是拉梅系數，ω=2πf表示角頻率(f是頻率)，s1和s2分別是震源的水平和垂直分量.若引入矩陣和向量記號

(2)

則(1)式可以表示為矩陣分塊形式

當使用NAD方法進行數值離散時，需要對(3)式分別沿x和z方向求偏導數(Liu et al., 2017b, 2018b)，則可得到微分方程組

(4)

與此同時，為了消除人工邊界處產生的反射波，需要采用吸收邊界條件.本文選取完美匹配層(Perfect Matched Layer, PML)邊界條件(Komatitsch and Tromp, 2003).具體地，首先引入復坐標(以x方向為例)

(5)

(6)

z方向也完全類似.將(4)式中的實坐標換成復坐標，再代入(5)和(6)式，則形成帶有PML吸收邊界條件的微分方程組

若令Δx和Δz分別為x和z方向上的空間離散步長，采用四階NAD方法所對應的網格差分模板(具體推導過程可見附錄A)進行離散，則(7a)式可以離散為

(8)

(9)

(7c)式離散為

(10)

其中I2表示二階單位矩陣.

若將二維計算區域中的網格節點按行排列，在每一行中按照從左向右的次序并且同一節點處的波場及其梯度值按照如下規則放置：

(11)

根據(8)—(10)式，可形成線性方程組

Cx=b,

(12)

阻抗矩陣C具有如下稀疏分塊結構：

(13)

值得注意的是，C是復值矩陣并且元素取值依賴于網格規模N，空間離散步長Δx,Δz，以及頻率f等因素.若設置參數N=41×41，Δx=Δz=0.02 km，f=15 Hz，則矩陣C的特征值分布如圖1所示(圖中虛線分別表示實軸和虛軸所在位置)，可看到C的特征值位于復平面內虛軸的兩側(實部有正有負)并且某些特征值非常接近于0.從數值代數的觀點來看，C是不定矩陣并且接近于奇異矩陣，病態程度較高.對于地震波常用頻率取值范圍內的其他參數情形，阻抗矩陣C的這些特點依然存在.這就為快速而穩定地求解線性方程組(12)帶來了困難，進而影響頻率域反演的整體計算效率.因此，構造高效的求解算法對于頻率域彈性波場模擬和全波形反演至關重要.

圖1 阻抗矩陣C的特征值分布Fig.1 Eigenvalue distribution of impedance matrix C

2 不精確旋轉分塊三角預處理迭代方法

2.1 不精確旋轉分塊三角預處理子的構造過程

為了快速有效地求解復線性方程組(12)，本文采用一種基于代數預處理子的預處理Krylov子空間迭代方法.首先分別提取方程(12)中矩陣和向量的實部與虛部，則有

Cx=(Cr+iCi)(xr+ixi)=b=br+ibi,

(14)

其中Cr,xr,br分別表示C,x,b的實部，Ci,xi,bi表示虛部.通過基本代數運算，將等式兩端的實部與虛部分別對應，則方程組(14)可以等價地轉化為實值線性方程組

(15)

針對(15)式中實矩陣A所具有的特殊分塊‘2×2’結構，可以構造預處理子M，使得預處理矩陣AM-1接近于單位矩陣并且條件數遠遠小于阻抗矩陣C，從而關于AM-1的線性方程組求解過程穩定并且需要較小的計算代價，達到加速求解的目的.此外，M-1在迭代計算過程中不必顯式求出，而只需求解以M為系數矩陣的廣義殘量方程組，這也稱作“預處理過程”(Saad, 2003).根據以上分析，預處理子M應當是A的一個好的近似，并且具備某種特殊結構，由此可利用A的子塊矩陣Cr和Ci分別構造“不精確旋轉分塊下三角(Inexact Rotated Block Lower Triangular, IRBLT)”預處理子

(16)

“不精確旋轉分塊上三角(Inexact Rotated Block Upper Triangular, IRBUT)”預處理子

(17)

以及“不精確旋轉分塊三角因子(Inexact Rotated Block Triangular Factor, IRBTF)”預處理子

(18)

其中

(19)

表示Givens正交旋轉矩陣，IN是N階單位矩陣，α>0為常數(通?？梢匀ˇ?1).(αCr+Ci)approx表示子塊矩陣αCr+Ci的某個可逆稀疏近似，可以通過Chebyshev半加速迭代公式(Bai et al., 2013)或者不完全LU分解(ILU)(Saad, 2003)來構造，預處理子命名中“不精確”一詞也是由此產生.

對于系數矩陣AM-1非對稱的一般型線性方程組，能夠有效求解的Krylov子空間迭代方法主要包括GMRES方法和BiCGSTAB方法(Sleijpen and Fokkema, 1993).相比于后者，GMRES方法只針對(15)式中系數矩陣A本身而不涉及法方程組，計算過程中無需AT的信息，從而具有更廣泛的適用性.此外，在這類三角分塊預處理子中，IRBTF預處理子所對應的預處理迭代方法具有最好的收斂特性，獲得給定精度解所需的迭代步數最少(Lang et al., 2020).然而IRBTF預處理子結構較為復雜，每一步預處理過程消耗過多的計算時間，從而導致求解所需的總時間高于IRBLT與IRBUT預處理子，實際計算效率偏低；通過數值算例測試，對于頻率域波動方程求解問題而言，IRBLT的計算效率略優于IRBUT預處理子.因此，在下文的數值試驗中統一采用IRBLT預處理子結合GMRES方法(簡記為IRBLT-GMRES)求解所對應的線性方程組.關于這類預處理迭代方法的具體實現步驟以及“不精確”預處理過程的處理方法可參見文獻(Lang and Ren, 2015).

2.2 不精確旋轉分塊三角預處理子的數值評估

我們采用IRBLT-GMRES方法求解線性方程組(15)，再與基于復轉換(complex-shifted)預處理子M0(Laird, 2001)以及不完全LU分解(Gander and Nataf, 2005)的預處理GMRES方法(分別簡記為M0-GMRES和ILU-GMRES)直接求解方程組(14)進行對比，考查各種方法在頻率域彈性波場模擬方面的數值效率.

考慮計算區域0≤x,z≤5 km，網格規模nx×nz=201×201, PML層數為20層.在此區域內P波速度vP=5 km·s-1, S波速度vS=3.0 km·s-1.實驗過程中，沿水平和垂直方向設置相等的空間離散步長，記作h=Δx=Δz=0.025 km;同時統一采用標準Ricker子波在垂直方向上施加震源作用力(劉少林等，2014)，它的頻率域表達式為(Lang and Yang, 2017)

(20)

其中Amp表示振幅，f0表示震源主頻，可取作20 Hz. 根據以上參數設置，可對阻抗矩陣以及震源右端項取值進行計算.在線性方程組的迭代求解過程中，初始迭代解均設置為0.停止準則是一旦當前迭代步殘量的歐氏范數小于或者等于初始殘量范數的10-6，即達到迭代終止條件.

表1列出了三種迭代方法在不同頻率取值下計算彈性波場所需要的時間，其中后兩列括號中列出的加速比定義為M0-GMRES或ILU-GMRES方法所需計算時間與IRBLT-GMRES計算時間的比值.可以看出所有的加速比均大于1，IRBLT-GMRES方法消耗的計算時間最少.相對于M0-GMRES和ILU-GMRES方法，IRBLT-GMRES方法最多可分別加速3.98倍和15.79倍.此外，通過按列觀察M0-GMRES方法和ILU-GMRES方法的計算時間，可以看到隨著頻率取值的增加，這兩種方法所需計算時間急劇增長，而關于ILU-GMRES方法的計算時間只隨著頻率增長而緩慢增加，計算過程最為穩定.

表1 三種迭代方法計算彈性波場所需時間(單位：s)Table 1 Computing time of three iteration methods for computing elastic wave-fields (unit：s)

3 波場模擬

我們運用四階NAD方法、四階OFD方法以及四階SG方法分別在均勻介質模型、雙層介質模型和Marmousi模型中進行波場模擬來對比各種方法的數值效率.由于從頻率域波場結果中較難判斷所得數值解的精確程度，需要對單頻波波場作離散Fourier逆變換(IDFT)合成時間域波場，通過觀察數值頻散現象來衡量這三種數值方法的準確性.根據對Ricker子波震源的頻譜分析(Liu et al., 2018b)，可采用0～2.5f0這一頻段的波場數據進行時間域波場計算.此外，本文還通過與解析解的波形對比和數值頻散分析進一步考查NAD方法在壓制數值頻散方面的優勢.

為了定量刻畫離散網格的疏密程度，這里定義網格頻率(grid frequency,Gf)為每個波長WL中的最小網格點數，可表示為(Liu, 1997)

(21)

在實際地震波的傳播過程中，最小波速vmin通常取作S波速度vS.

3.1 均勻模型

首先進行均勻介質中的波場模擬，基本計算參數設置見表2，P波速度vP=5.6 km·s-1, S波速度vS=3.0 km·s-1. 在此參數設置下，可計算出網格頻率Gf=3.

表2 對各種介質模型進行波場模擬的基本參數表Table 2 Parameter table for wave-field simulation for various media models

圖2顯示由四階NAD方法計算得到的單頻波波場快照，其中第一行子圖表示水平分量，第二行是垂直分量，每列子圖分別對應于頻率f=10,20,30 Hz. 由此可以合成時間域波場，相應的計算過程也對四階OFD方法和四階SG方法同樣進行.在t=0.25 s時刻由三種方法計算得到的時間域波場結果如圖3所示，可以看到NAD方法計算的波場快照沒有明顯數值頻散，而對于另外兩種方法，則有比較清楚的頻散現象，其中OFD方法對應的頻散情況最為嚴重.

圖2 均勻介質中由四階NAD方法計算得到的單頻波波場快照(a),(c),(e) 水平分量; (b),(d),(f) 垂直分量; (a)和(b), (c)和(d), (e)和(f)分別對應于頻率取值f=10, 20, 30 HzFig.2 Mono-frequency wave-field snapshots by fourth-order NAD method in homogeneous mediumWhere (a), (c), (e) are the horizontal components, (b), (d), (f) are the vertical components, and (a) & (b), (c) & (d), (e) & (f) correspond to the frequency value f=10, 20, 30 Hz, respectively.

圖3 均勻介質中由四階NAD (a&b)、OFD (c&d)、SG (e&f)方法計算t=0.25 s時刻的時間域波場快照(a),(c),(e) 水平分量; (b),(d),(f) 垂直分量.Fig.3 Time-domain wave-field snapshots at t=0.25 s by fourth-order NAD (a&b), OFD (c&d), and SG (e&f) methods in homogeneous mediumWhere (a), (c), (e) are the horizontal components, and (b), (d), (f) are the vertical components.

為了進一步驗證所討論數值方法的準確性，我們放置接收器于(1.0 km,2.0 km)處，合成對應于三種方法的正則化解與解析解(王美霞等，2012)的波形對比如圖4所示.可以看出由NAD方法計算的數值解與解析解的波形匹配程度最高，而對于另外兩種方法，則存在比較明顯的偏差.

圖4 由四階NAD方法(a&b)、四階OFD方法(c&d)和四階SG方法(e&f)計算的正則化數值解與解析解的波形對比(a),(c),(e) 水平分量; (b),(d),(f) 垂直分量.Fig.4 Waveform comparisons of normalized numerical solutions by fourth-order NAD (a&b), OFD (c&d), and SG (e&f) methods with analytic solutions(a),(c),(e) Horizontal components; (b), (d), (f) Vertical components.

若使OFD和SG方法能夠壓制數值頻散，需要加密相應的離散網格.當OFD方法的網格參數設置為h=0.015 km,nx=nz=267，SG方法的參數設置為h=0.017 km,nx=nz=233，相應的波場快照如圖5所示，其中位于左列的子圖對應于OFD方法，右列對應于SG方法，第一行為水平分量，第二行為垂直分量，通過觀察未發現明顯的數值頻散.比較三種方法在成功壓制數值頻散時所需的計算時間(見表3)，可以看到NAD方法用時最少，相比于次之的OFD方法可以節省8.9%，比SG方法減少20.4%.

圖5 均勻介質中在細網格下由OFD(a&b)和SG(c&d)方法計算t=0.25 s時刻的時間域波場快照(a)和(c) 水平分量; (b)和(d) 垂直分量.Fig.5 Time-domain wave-field snapshots at t=0.25 s computed by OFD (a&b) and SG (c&d) methods on finer grids in homogeneous medium(a) (c) The horizontal components; (b) (d) The vertical components.

表3 各種數值方法獲得準確解所需的最少計算時間(單位：min)Table 3 Least computing time to obtain accurate solutions by various numerical schemes (unit：min)

3.2 數值頻散分析

圖6 縱橫波速比r=2(a,c,e)和r=3(b,d,f),波傳播角θ分別為30°(a&b)，45°(c&d)，60°(e&f)下的S波頻散曲線對比Fig.6 S-wave dispersion curve comparison for r=2 (a,c,e) & r=3 (b,d,f), and the wave propagation angle θ valuing 30° (a&b), 45° (c&d), and 60° (e&f)

3.3 雙層介質模型

圖7顯示了雙層介質中由四階NAD方法計算得到的單頻波波場快照，其中第一行子圖為水平分量，第二行為垂直分量，每列子圖分別對應頻率取值為f=15,25,35 Hz.由此合成對應于三種方法在t=0.25 s時刻的時間域波場快照如圖8所示.類似于均勻介質情形，四階NAD方法的波場快照結果最為準確，而另外兩種方法則存在明顯的數值頻散，SG方法的頻散情況略好于OFD方法.

圖7 雙層介質中由四階NAD方法計算得到的單頻波波場快照(a),(c),(e) 水平分量; (b),(d),(f) 垂直分量，(a)和(b), (c)和(d), (e)和(f)分別對應于頻率取值f=15,25,35 Hz.Fig.7 Mono-frequency wave-field snapshots by fourth-order NAD method in two-layer medium(a),(c),(e) The horizontal components; (b),(d),(f) The vertical components, and (a)&(b), (c)&(d), (e)&(f) correspond to the frequency values f=15,25,35 Hz, respectively.

圖8 雙層介質中由四階NAD(a&b)、OFD(c&d)、SG(e&f)方法計算t=0.25 s時刻的時間域波場快照(a),(c),(e) 水平分量; (b),(d),(f) 垂直分量.Fig.8 Time-domain wave-field snapshots at t=0.25 s by fourth-order NAD (a&b), OFD (c&d), and SG (e&f) methods in two-layer medium(a),(c),(e) The horizontal components; (b),(d),(f) The vertical components.

同樣地，為了使四階OFD方法和SG方法能夠成功地壓制數值頻散，分別加密離散網格至參數設置為h=0.0147 km,nx=nz=272以及h=0.0169 km,nx=nz=237, 相應的波場快照見圖9，其中左右兩列子圖分別對應于OFD和SG方法，第一行子圖為水平分量，第二行為垂直分量，通過觀察未看到明顯的數值頻散.比較三種方法在得到準確波場結果的前提下所花費的最少計算時間(如表3所列)，可以看出NAD方法用時最少，相比于OFD方法和SG方法，可以分別加速9.1%和19.6%.

圖9 雙層介質中在細網格下由OFD(a&b)和SG(c&d)方法計算t=0.25 s時刻的時間域波場快照(a)和(c) 水平分量; (b)和(d) 垂直分量.Fig.9 Time-domain wave-field snapshots at t=0.25 s computed by OFD (a&b) and SG (c&d) methods on finer grids in two-layer medium(a)&(c) The horizontal components; (b)&(d) The vertical components.

3.4 Marmousi模型

圖10 Marmousi模型S波速度結構Fig.10 S-wave velocity of Marmousi model

根據四階NAD方法計算得到的單頻波波場快照如圖11所示，其中左列子圖為水平分量，右列為垂直分量，每一行分別對應頻率取值為6 Hz、12 Hz、18 Hz. 若在地表層每個網格節點處放置接收器，所得理論地震圖如圖12所示.此外，在t=0.4 s、0.8 s、1.2 s時刻的時間域波場快照見圖13.通過觀察圖12和13，可以看到四階NAD方法成功壓制數值頻散，波場結果符合實際情形.

圖11 Marmousi模型中由四階NAD方法計算頻率取值分別為6 Hz(a&b)，12 Hz(c&d)，18 Hz(e&f)的單頻波波場快照(a),(c),(e) 水平分量; (b),(d),(f) 垂直分量.Fig.11 Mono-frequency wave-field snapshots at 6 Hz (a&b), 12 Hz (c&d), 18 Hz (e&f) by fourth-order NAD method in Marmousi model(a),(c),(e) The horizontal components; (b),(d),(f) The vertical components.

圖12 Marmousi 模型由四階NAD方法計算的理論地震圖(a) 水平分量; (b) 垂直分量.Fig.12 Theoretical seismograms by fourth-order NAD method in Marmousi model(a) and (b) The horizontal and vertical components, respectively.

圖13 Marmousi模型由四階NAD方法計算t=0.3 s (a&b), 0.6 s (c&d), 1.2 s (e&f)時刻的時間域波場快照(a),(c),(e) 水平分量; (b),(d),(f) 垂直分量.Fig.13 Time-domain wave-field snapshots at t=0.3 s (a&b), 0.6 s (c&d), 1.2 s (e&f) by fourth-order NAD method in Marmousi model(a),(c),(e) are the horizontal components, and (b),(d),(f) are the vertical components.

4 結論

在NAD方法成功運用于頻率域聲波方程數值模擬以及全波形反演之后，本文進一步研究了這種方法在頻率域彈性波場模擬方面的數值效率.對于各向同性介質中的頻率域彈性波動方程，我們詳細推導了四階NAD方法的離散過程并得到大型線性代數方程組.針對阻抗矩陣稀疏分塊結構與數學性質的深入分析結果，揭示了對該線性方程組進行快速有效求解的本質困難.為此，本文引入IRBLT-GMRES預處理迭代求解算法，并將其運用于彈性波場計算.相比于經典的M0-GMRES和ILU-GMRES預處理迭代方法，IRBLT-GMRES方法最高可分別加速頻率域彈性波場模擬大致3.98倍和15.79倍.接著我們采用四階NAD方法分別在均勻介質和雙層介質中進行波場模擬、數值頻散分析以及與解析解的波形對比，并與四階OFD和四階SG方法進行比較.各種數值結果均顯示了NAD方法在壓制數值頻散和提高計算效率方面的優勢，相比于OFD方法和SG方法，可分別加速約9%和20%.此外，四階NAD方法還用于Marmousi模型的彈性波場模擬并得出無明顯數值頻散的地表觀測記錄與波場快照，這進一步說明了NAD方法在復雜介質模型數值模擬中依然可靠.由此可見，IRBLT-GMRES迭代方法配合NAD方法可作為頻率域彈性波數值模擬與全波形反演的有力工具.

附錄A 四階NAD方法網格差分模板的具體構造過程

以x方向為例，若對波場u在(i,j)點處的二階和三階偏導數項進行離散，NAD格式所對應的網格差分模板可表示為

(A1)

(A2)

和

(A3)

求解(A2)和(A3)式，可得網格差分模板系數

(A4)

將其代入(A1)式，即可得到四階NAD方法的離散格式，z方向也可類似推導.

附錄B 阻抗矩陣的非零子塊元素表達式

(B1)

(B2)

(B3)

(B4)

(B5)

(B6)

(B7)

(B8)

(B9)

其中j為分塊行號，對應于計算區域內網格節點的不同位置，上述表達式中dx和dz的取值也隨之變化.

附錄C 關于頻率域彈性波方程四階NAD方法的頻散分析

類似于聲波情形(Lang and Yang, 2017)，考慮齊次頻率域彈性波動方程及其偏導形式

(C1)

采用四階NAD方法所對應的網格差分模板對(C1)式進行數值離散，再代入平面諧波解U=U0e-i(kx·x+kz·z)，其中U0=[u0,v0]T表示常向量，kx,kz分別是波數k在x和z方向上的分量，則有

(C2)

矩陣

(C3)

的元素可以表示為

(C4)

(C5)

(C6)

(C7)

(C8)

(C9)

(C10)

(C11)

(C12)

(C13)

(C14)

(C15)

事實上，對于P波的頻散分析完全類似，所得結果只相差一個因子r2.