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

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

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

0 引言

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

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

類似于時(shí)間域情形，頻率域成像過程的主要計(jì)算量也集中于對地震波傳播規(guī)律的正演模擬.因此，正演算法的優(yōu)劣直接決定著全波形成像結(jié)果的精度和整體計(jì)算效率.目前用于地震波場模擬的主流正演算法包括有限差分方法(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)等.有限差分方法因其原理直觀、格式構(gòu)造簡單、實(shí)現(xiàn)容易、便于并行等優(yōu)點(diǎn)，已成為地震波模擬和成像的有力工具之一(Yang et al., 2006; 劉少林等，2013).然而傳統(tǒng)有限差分方法容易引起數(shù)值頻散效應(yīng)，當(dāng)遇到復(fù)雜地質(zhì)結(jié)構(gòu)(例如強(qiáng)間斷面)時(shí)，往往需要加密離散網(wǎng)格以提升數(shù)值精度，這將導(dǎo)致計(jì)算效率降低(Moczo et al., 2000; Liu et al., 2017c).近似解析離散化(Nearly Analytic Discrete，NAD)方法作為一種新型差分方法，由楊頂輝等(Yang et al., 2003)首先引入到時(shí)間域地震波場數(shù)值計(jì)算中.該方法的基本構(gòu)造思想是同時(shí)采用波場值以及波場梯度值來近似高階偏導(dǎo)數(shù)項(xiàng)(Yang et al., 2004)，這不僅繼承了普通有限差分方法的優(yōu)點(diǎn)，同時(shí)也能在粗網(wǎng)格下有效地壓制數(shù)值頻散，進(jìn)而得到準(zhǔn)確的數(shù)值結(jié)果，因此該方法可以提高地震波數(shù)值模擬的計(jì)算效率(Tong et al., 2013; Liu et al., 2015).

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

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

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

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

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

(1)

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

(2)

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

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

(4)

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

(5)

(6)

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

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

(8)

(9)

(7c)式離散為

(10)

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

若將二維計(jì)算區(qū)域中的網(wǎng)格節(jié)點(diǎn)按行排列，在每一行中按照從左向右的次序并且同一節(jié)點(diǎn)處的波場及其梯度值按照如下規(guī)則放置：

(11)

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

Cx=b,

(12)

阻抗矩陣C具有如下稀疏分塊結(jié)構(gòu)：

(13)

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

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

2 不精確旋轉(zhuǎn)分塊三角預(yù)處理迭代方法

2.1 不精確旋轉(zhuǎn)分塊三角預(yù)處理子的構(gòu)造過程

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

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

(14)

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

(15)

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

(16)

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

(17)

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

(18)

其中

(19)

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

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

2.2 不精確旋轉(zhuǎn)分塊三角預(yù)處理子的數(shù)值評估

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

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

(20)

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

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

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

3 波場模擬

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

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

(21)

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

3.1 均勻模型

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

表2 對各種介質(zhì)模型進(jìn)行波場模擬的基本參數(shù)表Table 2 Parameter table for wave-field simulation for various media models

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

圖2 均勻介質(zhì)中由四階NAD方法計(jì)算得到的單頻波波場快照(a),(c),(e) 水平分量; (b),(d),(f) 垂直分量; (a)和(b), (c)和(d), (e)和(f)分別對應(yīng)于頻率取值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 均勻介質(zhì)中由四階NAD (a&b)、OFD (c&d)、SG (e&f)方法計(jì)算t=0.25 s時(shí)刻的時(shí)間域波場快照(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.

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

圖4 由四階NAD方法(a&b)、四階OFD方法(c&d)和四階SG方法(e&f)計(jì)算的正則化數(shù)值解與解析解的波形對比(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方法能夠壓制數(shù)值頻散，需要加密相應(yīng)的離散網(wǎng)格.當(dāng)OFD方法的網(wǎng)格參數(shù)設(shè)置為h=0.015 km,nx=nz=267，SG方法的參數(shù)設(shè)置為h=0.017 km,nx=nz=233，相應(yīng)的波場快照如圖5所示，其中位于左列的子圖對應(yīng)于OFD方法，右列對應(yīng)于SG方法，第一行為水平分量，第二行為垂直分量，通過觀察未發(fā)現(xiàn)明顯的數(shù)值頻散.比較三種方法在成功壓制數(shù)值頻散時(shí)所需的計(jì)算時(shí)間(見表3)，可以看到NAD方法用時(shí)最少，相比于次之的OFD方法可以節(jié)省8.9%，比SG方法減少20.4%.

圖5 均勻介質(zhì)中在細(xì)網(wǎng)格下由OFD(a&b)和SG(c&d)方法計(jì)算t=0.25 s時(shí)刻的時(shí)間域波場快照(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 各種數(shù)值方法獲得準(zhǔn)確解所需的最少計(jì)算時(shí)間(單位：min)Table 3 Least computing time to obtain accurate solutions by various numerical schemes (unit：min)

3.2 數(shù)值頻散分析

圖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 雙層介質(zhì)模型

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

圖7 雙層介質(zhì)中由四階NAD方法計(jì)算得到的單頻波波場快照(a),(c),(e) 水平分量; (b),(d),(f) 垂直分量，(a)和(b), (c)和(d), (e)和(f)分別對應(yīng)于頻率取值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 雙層介質(zhì)中由四階NAD(a&b)、OFD(c&d)、SG(e&f)方法計(jì)算t=0.25 s時(shí)刻的時(shí)間域波場快照(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方法能夠成功地壓制數(shù)值頻散，分別加密離散網(wǎng)格至參數(shù)設(shè)置為h=0.0147 km,nx=nz=272以及h=0.0169 km,nx=nz=237, 相應(yīng)的波場快照見圖9，其中左右兩列子圖分別對應(yīng)于OFD和SG方法，第一行子圖為水平分量，第二行為垂直分量，通過觀察未看到明顯的數(shù)值頻散.比較三種方法在得到準(zhǔn)確波場結(jié)果的前提下所花費(fèi)的最少計(jì)算時(shí)間(如表3所列)，可以看出NAD方法用時(shí)最少，相比于OFD方法和SG方法，可以分別加速9.1%和19.6%.

圖9 雙層介質(zhì)中在細(xì)網(wǎng)格下由OFD(a&b)和SG(c&d)方法計(jì)算t=0.25 s時(shí)刻的時(shí)間域波場快照(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波速度結(jié)構(gòu)Fig.10 S-wave velocity of Marmousi model

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

圖11 Marmousi模型中由四階NAD方法計(jì)算頻率取值分別為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方法計(jì)算的理論地震圖(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方法計(jì)算t=0.3 s (a&b), 0.6 s (c&d), 1.2 s (e&f)時(shí)刻的時(shí)間域波場快照(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 結(jié)論

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

附錄A 四階NAD方法網(wǎng)格差分模板的具體構(gòu)造過程

以x方向?yàn)槔?，若對波場u在(i,j)點(diǎn)處的二階和三階偏導(dǎo)數(shù)項(xiàng)進(jìn)行離散，NAD格式所對應(yīng)的網(wǎng)格差分模板可表示為

(A1)

(A2)

和

(A3)

求解(A2)和(A3)式，可得網(wǎng)格差分模板系數(shù)

(A4)

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

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

(B1)

(B2)

(B3)

(B4)

(B5)

(B6)

(B7)

(B8)

(B9)

其中j為分塊行號，對應(yīng)于計(jì)算區(qū)域內(nèi)網(wǎng)格節(jié)點(diǎn)的不同位置，上述表達(dá)式中dx和dz的取值也隨之變化.

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

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

(C1)

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

(C2)

矩陣

(C3)

的元素可以表示為

(C4)

(C5)

(C6)

(C7)

(C8)

(C9)

(C10)

(C11)

(C12)

(C13)

(C14)

(C15)

事實(shí)上，對于P波的頻散分析完全類似，所得結(jié)果只相差一個(gè)因子r2.