三維彈性波方程的修正時(shí)空優(yōu)化保辛數(shù)值求解方法

王健， 賀茜君， 董興朋， 楊頂輝*， 李靜爽， 黃雪源， 周艷杰

1 清華大學(xué)數(shù)學(xué)科學(xué)系， 北京 100084 2 北京工商大學(xué)數(shù)學(xué)與統(tǒng)計(jì)學(xué)院應(yīng)用統(tǒng)計(jì)系， 北京 100048 3 中國(guó)礦業(yè)大學(xué)(北京)理學(xué)院， 北京 100083

0 引言

地震波場(chǎng)正演模擬是指利用數(shù)值方法求解地震波在已知地下介質(zhì)中傳播的技術(shù)，它是基于波動(dòng)方程反問題研究的基礎(chǔ).此外，正演模擬對(duì)于解釋地震信息，評(píng)估觀測(cè)系統(tǒng)等也具有十分重要的意義.地震波正演模擬有很多實(shí)現(xiàn)方法，包括有限差分方法(finitie difference, FD, Dablain, 1986)、有限元法(finite element, FE, Lysmer and Drake, 1972; Smith, 1975; He et al., 2020)、偽譜法(pseudo-spectral, PS, Kosloff and Baysal, 1982)、譜元法(spectral-element, SEM, Priolo and Seriani, 1991)等.這些方法各有優(yōu)缺點(diǎn)及適用范圍，其中，有限差分方法是目前應(yīng)用最為廣泛的正演模擬算法.

有限差分方法將計(jì)算區(qū)域劃分為有限個(gè)規(guī)則的或不規(guī)則的網(wǎng)格點(diǎn)，利用泰勒展開的方式將波動(dòng)方程中的微分算子轉(zhuǎn)化為網(wǎng)格點(diǎn)上函數(shù)值的差分形式，通過求解網(wǎng)格點(diǎn)上的函數(shù)值近似得到原波動(dòng)方程的解.不同的差分形式對(duì)應(yīng)著不同的有限差分方法，其具有易于實(shí)現(xiàn)、存儲(chǔ)量小、運(yùn)算速度快、易于并行等優(yōu)點(diǎn).Alterman 和Karal(1969)最早將該方法引入到地震波場(chǎng)的正演模擬中，并在后續(xù)得到了廣泛應(yīng)用(如, Moczo et al., 2002).然而，有限差分方法面臨著可能產(chǎn)生嚴(yán)重?cái)?shù)值頻散的問題(Fei and Larner, 1995).這是由于利用差分形式近似微分算子時(shí)，不同頻率波的相速度產(chǎn)生了差別，特別是當(dāng)網(wǎng)格步長(zhǎng)很大或波場(chǎng)梯度很大時(shí)，這種現(xiàn)象尤為明顯.因此當(dāng)波傳播時(shí)間較長(zhǎng)時(shí)，不同頻率的波由于相速度不同，會(huì)發(fā)生分離.這顯然與實(shí)際物理定律相違背，故而產(chǎn)生了虛假信息，影響了正演模擬的精度.

為解決這一問題，Yang 等(2003)將近似解析離散化(nearly analytic discrete, NAD)算子引入波動(dòng)方程的空間偏導(dǎo)數(shù)離散中，同時(shí)利用網(wǎng)格點(diǎn)上的位移及梯度的差分來(lái)近似微分算子.由于利用了網(wǎng)格點(diǎn)上的梯度值，因此與相同數(shù)值精度的傳統(tǒng)有限差分方法相比，NAD算子長(zhǎng)度更小、算子更為緊致.同時(shí)，數(shù)值實(shí)驗(yàn)表明，在相同網(wǎng)格步長(zhǎng)的情況下，NAD算子的數(shù)值頻散明顯低于傳統(tǒng)的有限差分方法，能夠在粗網(wǎng)格情形下有效壓制數(shù)值頻散(Yang et al., 2003).近年來(lái)，基于NAD算子的思想，一系列改進(jìn)算法被提出，包括優(yōu)化近似解析離散化方法(optimal nearly analytic discrete, ONAD, Yang et al., 2006)、帶權(quán)重的近似解析離散化方法(weighted nearly analytic discrete, WNAD, 王磊等, 2009)、近似解析保辛分部龍格-庫(kù)塔方法(nearly analytic symplectic partitioned Runge-Kutta, NSPRK, 馬嘯等, 2010; Ma et al., 2011)、改善近似解析離散方法(refined nearly analytic discrete, RNAD, Yang et al., 2010)等，并已被應(yīng)用于全波形成像(Wang et al., 2019).

壓制數(shù)值頻散的另一種思路是提高離散格式相速度的準(zhǔn)確性.為達(dá)此目的，可以構(gòu)造優(yōu)化有限差分格式，將數(shù)值波數(shù)和真實(shí)波數(shù)的誤差作為優(yōu)化函數(shù)，從而得到最優(yōu)的差分離散系數(shù)(Haras and Ta'asan, 1994; Kim and Lee, 1996; Lee and Seo, 2002).其中，Zhang和Yao(2013)將數(shù)值波數(shù)和真實(shí)波數(shù)的無(wú)窮范數(shù)誤差作為目標(biāo)函數(shù)，通過模擬退火法(Kirkpatrick et al., 1983; Sen and Stoffa, 2013)優(yōu)化目標(biāo)函數(shù)，最終得到最優(yōu)的差分離散系數(shù).數(shù)值結(jié)果表明，使用優(yōu)化的高階有限差分方法可以大大節(jié)省內(nèi)存需求和計(jì)算代價(jià).

波動(dòng)方程可以被視為一個(gè)無(wú)限維線性可分的哈密爾頓系統(tǒng)(Hamilton system)(Schmid, 1987)，其相空間存在著辛幾何結(jié)構(gòu)，可以表示物體的運(yùn)動(dòng)規(guī)律.在哈密爾頓體系下構(gòu)造的辛算法(symplectic method)可以保持哈密爾頓系統(tǒng)的辛幾何結(jié)構(gòu)，具有重要意義.Feng和 Qin(1987)根據(jù)生成函數(shù)理論，給出了任意階精度的隱式保辛數(shù)值格式的構(gòu)造方法，并據(jù)此提出了多種不同類型的辛格式，包括保辛龍格-庫(kù)塔(Runge-Kutta, RK)格式、蛙跳格式等(馮康等, 2003).Okunbor和Skeel(1992)給出了對(duì)于可分哈密爾頓系統(tǒng)的分部龍格-庫(kù)塔(patitioned Runge-Kutta, PRK)格式的構(gòu)造方法.Qin和Zhang(1990)構(gòu)造了一至四階精度、一至四級(jí)的PRK格式.Ma 等(2011)將保辛格式與NAD算子相結(jié)合，給出了NSPRK方法，其數(shù)值頻散誤差更小.Liu 等(2017)提出了一種修正的保辛分部龍格-庫(kù)塔(symplectic partitioned Runge-Kutta, SPRK)格式，使一個(gè)二級(jí)格式達(dá)到了三階時(shí)間精度.Ma和Yang(2017)提出了優(yōu)化SPRK格式，具有保相位的優(yōu)點(diǎn).Ma等(2018)提出了時(shí)空優(yōu)化保辛(time-space optimized symplectic method, TSOS)方法，將優(yōu)化SPRK格式與優(yōu)化有限差分格式結(jié)合，同時(shí)具有保相位和低數(shù)值頻散的優(yōu)點(diǎn)，更易于處理非均勻介質(zhì).

本文將修正的SPRK方法和優(yōu)化有限差分方法思想結(jié)合，發(fā)展了用于求解三維非均勻介質(zhì)中彈性波方程的修正時(shí)空優(yōu)化保辛方法(modified time-space optimized symplectic method, MTSOS)，并分析了該方法的穩(wěn)定性及數(shù)值頻散特性，數(shù)值實(shí)驗(yàn)表明該方法能有效壓制數(shù)值頻散，同時(shí)可以準(zhǔn)確、高效地模擬彈性波在地下介質(zhì)中的傳播.

1 修正時(shí)空優(yōu)化保辛(MTSOS)方法

1.1 時(shí)間格式

一般的2n維哈密爾頓系統(tǒng)可表示為：

(1)

(2)

且f和g分別是關(guān)于u和v的線性函數(shù)，則稱系統(tǒng)(1)為線性可分的哈密爾頓系統(tǒng).Qin和Zhang(1990)指出，波動(dòng)方程是一個(gè)線性可分的哈密爾頓系統(tǒng).對(duì)于彈性波方程：

(3)

可將其表示為線性可分的哈密爾頓系統(tǒng)的形式：

(4)

其中

(5)

(6)

(7)

空間算子L的其他項(xiàng)類似可得.

由Sanz-Serna(1988),s階顯式SPRK格式可表示為：

(8)

其中，un和vn表示該系統(tǒng)在第n個(gè)時(shí)間步的數(shù)值解，c=(c1,c2,…,cs)和d=(d1,d2,…,ds)是該格式的系數(shù)向量，Δt表示該格式的時(shí)間步長(zhǎng).

為保證該格式具有s階精度，需通過泰勒展開的方式確定系數(shù)向量c和d.但要想使格式達(dá)到s階精度，至少需要s級(jí)格式才能滿足要求.Liu 等(2017)提出了一種修正SPRK格式，具體形式為：

(9)

因?yàn)楦袷?9)中引入了修正項(xiàng)Lv1，可以通過泰勒展開的方式確定系數(shù)向量c和d，使一個(gè)二級(jí)格式具有三階時(shí)間精度.

經(jīng)計(jì)算，修正SPRK格式的系數(shù)向量為：

(10)

1.2 空間格式

為了利用修正SPRK格式得到波動(dòng)方程的數(shù)值解，還需要離散式(4)中的空間算子L.有限差分法方法易于編程及并行，是一種常見的數(shù)值離散方法.一般的有限差分格式對(duì)函數(shù)的空間高階偏導(dǎo)數(shù)利用該點(diǎn)及周圍點(diǎn)的函數(shù)值進(jìn)行離散，即

(11)

(12)

(13)

(14)

(15)

(16)

1.3 等效介質(zhì)參數(shù)與邊界條件

(17)

(18)

可以看出，在均勻介質(zhì)中，密度和彈性參數(shù)的等效介質(zhì)參數(shù)值即為其真實(shí)值.在非均勻介質(zhì)中，只需在網(wǎng)格附近利用算術(shù)或調(diào)和平均即可得到等效介質(zhì)參數(shù)值.Ma 等(2018)通過數(shù)值算例證明了等效介質(zhì)參數(shù)處理可以更好地模擬間斷面處地震波的傳播.

自由地表邊界條件是指地表處的應(yīng)力為0(蘭海強(qiáng)等, 2012).為了模擬自由地表邊界條件，Nilsson等(2007)提出了邊界修正法，該方法是一種顯式方法，具有二階空間精度，并具有很好的數(shù)值穩(wěn)定性.本文以三維彈性波方程為例說(shuō)明該方法.

假設(shè)z=0處為自由地表，滿足自由地表邊界條件.為下文表述方便，記(u1,u2,u3)T=(u,v,w)T.在三維彈性波方程中，自由地表邊界條件可表示為：

(19)

設(shè)z=0對(duì)應(yīng)k=1，邊界修正法首先引入一個(gè)虛擬層k=0，通過數(shù)值離散自由地表邊界條件(19)來(lái)得到虛擬層的波場(chǎng).具體來(lái)說(shuō)，用如下格式來(lái)離散自由地表邊界條件(Nilsson et al., 2007)：

(20)

(21)

(22)

對(duì)于其他的k=1處關(guān)于z方向的混合偏導(dǎo)數(shù)，可以類似地得到其數(shù)值離散格式.Nilsson 等(2007)指出，上述的邊界修正法是一種二階方法，且當(dāng)格式滿足CFL條件(Courant-Friedrichs-Lewy condition)時(shí)是穩(wěn)定的.

需要注意的是，由于邊界修正法是一種二階方法，所以當(dāng)計(jì)算區(qū)域內(nèi)部的空間精度超過二階時(shí)，會(huì)造成邊界與計(jì)算區(qū)域內(nèi)精度不匹配的現(xiàn)象，影響計(jì)算區(qū)域內(nèi)的空間精度.為解決這一問題，可以在計(jì)算區(qū)域內(nèi)從邊界層開始逐層提高空間精度，直至達(dá)到計(jì)算區(qū)域內(nèi)的空間精度，以保證數(shù)值格式的穩(wěn)定性.

由于計(jì)算時(shí)間和規(guī)模的限制，一般只能模擬地震波在給定區(qū)域內(nèi)的傳播情況.在這種情況下，需要在計(jì)算區(qū)域的人工邊界添加吸收邊界條件(absorbing boundary conditions)，來(lái)消除人工邊界產(chǎn)生的反射波對(duì)計(jì)算區(qū)域的影響.近年來(lái)很多吸收邊界條件方法被提出，其中包括旁軸近似方法(Clayton and Engquist, 1977)，完全匹配層吸收邊界條件(Berenger, 1994)等.本文采用Martin等(2010)提出的輔助微分方程-完全匹配層(auxiliary differential equation-perfectly matched layer, ADE-PML)吸收邊界條件來(lái)處理人工邊界問題.

將1.1節(jié)提出的修正SPRK時(shí)間格式與1.2節(jié)提出的優(yōu)化有限差分格式相結(jié)合，將空間格式拓展到三維，同時(shí)利用等效介質(zhì)參數(shù)，邊界修正法和完全匹配層吸收邊界條件，即獲得了三維非均勻介質(zhì)彈性波方程正演模擬的修正時(shí)空優(yōu)化保辛(MTSOS)方法.該方法是一種保辛方法，且由于使用了優(yōu)化有限差分格式和等效介質(zhì)參數(shù)，更適用于非均勻介質(zhì)中彈性波方程的求解.在第2節(jié)，將對(duì)三維MTSOS方法的穩(wěn)定性和數(shù)值頻散進(jìn)行分析.

2 三維MTSOS方法的穩(wěn)定性和數(shù)值頻散分析

2.1 穩(wěn)定性條件

本節(jié)首先利用傅里葉分析的方法，求出時(shí)間離散算子的譜半徑，進(jìn)而得到MTSOS方法的穩(wěn)定性條件.對(duì)于三維均勻介質(zhì)中的聲波方程，將諧波解

(23)

代入聲波方程，其中ωnum表示數(shù)值角頻率，k=(kx,ky,kz)表示波矢量，Δt表示時(shí)間步長(zhǎng)，h=Δx=Δy=Δz表示空間步長(zhǎng).可以得到：

(24)

其中增長(zhǎng)矩陣G為：

(25)

I表示單位算子，L表示空間離散算子.設(shè)矩陣G的特征值為ψ，算子L的特征值為.由于波動(dòng)方程是雙曲型偏微分方程，故其空間離散算子<0.由式(25)，ψ的特征多項(xiàng)式f(ψ)和滿足關(guān)系：

f(ψ)=ψ2-tr(G)ψ+det(G)=ψ2-(2+Δt2

(26)

為保證格式穩(wěn)定，需要滿足|ψ|≤1，即：

(27)

解上述不等式，可以得到：

-12≤Δt2≤0.

(28)

(29)

式(29)對(duì)算子L的任意特征值均成立，因此要保證格式穩(wěn)定，最大庫(kù)朗數(shù)αmax需滿足：

(30)

式(30)即為三維MTSOS方法的穩(wěn)定性條件，可以根據(jù)不同精度的空間算子L的最小特征值得到與之對(duì)應(yīng)的最大庫(kù)朗數(shù).在表1中給出了三維情形下不同精度的MTSOS及TSOS方法的最大庫(kù)朗數(shù).可以看出，MTSOS方法的穩(wěn)定性相比TSOS方法有了明顯的提升.

表1 三維情形下不同空間精度的MTSOS方法與TSOS方法的最大庫(kù)朗數(shù)Table 1 Maximum Courant numbers of MTSOS method and TSOS method with different spatial accuracies in 3D case

2.2 數(shù)值頻散分析

為了對(duì)格式進(jìn)行數(shù)值頻散分析，首先定義數(shù)值頻散比率為如下形式(Yang et al., 2006)：

(31)

將諧波解式(23)代入式(24)，可以得到：

(32)

由于方程必有非零解，因此有：

det(exp(iωnumΔt)I-G)=0,

(33)

即：

ωnumΔt=arccos(Re(ψ)).

(34)

又由式(26)可知：

(35)

因此，可得：

(36)

式(36)被稱為數(shù)值頻散關(guān)系式，它表示了數(shù)值頻散比率R與庫(kù)朗數(shù)、采樣率及空間算子的特征值之間的關(guān)系.R越接近于1，說(shuō)明數(shù)值頻散越小.接下來(lái)，本小節(jié)將比較三維情形下MTSOS方法與SPRK方法的數(shù)值頻散比率，以說(shuō)明MTSOS方法的優(yōu)勢(shì).

在圖1和圖2中分別給出了在α=0.1,Sp=0.4時(shí)不同空間精度的MTSOS方法和SPRK方法在不同方向上的頻散誤差|1-R|.從圖中可以看出，MTSOS方法與SPRK方法的數(shù)值頻散誤差都表現(xiàn)出了各向異性的特點(diǎn)：數(shù)值頻散誤差的最大值都出現(xiàn)在沿坐標(biāo)軸的方向，而在與坐標(biāo)軸呈45°角的方向，數(shù)值頻散誤差最小.同時(shí)，隨著空間精度的提升，數(shù)值頻散誤差也呈現(xiàn)出下降的趨勢(shì).

圖1 不同空間精度的MTSOS方法在α=0.1,Sp=0.4時(shí)不同方向的數(shù)值頻散誤差(a) 六階； (b) 八階； (c) 十階； (d) 十二階.Fig.1 The numerical dispersion error of the MTSOS method with different spatial accuracies in different directions when α=0.1,Sp=0.4(a) Sixth order; (b) Eighth order; (c) Tenth order; (d) Twelfth order.

圖2 不同空間精度的SPRK方法在α=0.1,Sp=0.4時(shí)不同方向的數(shù)值頻散誤差(a) 六階； (b) 八階； (c) 十階； (d) 十二階.Fig.2 The numerical dispersion error of the SPRK method with different spatial accuracies in different directions when α=0.1,Sp=0.4(a) Sixth order; (b) Eighth order; (c) Tenth order; (d) Twelfth order.

從圖中還可以看出，從總體上看，MTSOS方法的數(shù)值頻散誤差要小于SPRK方法.為了量化比較兩種方法的數(shù)值頻散誤差，在表2中分別給出了不同空間精度的MTSOS方法與SPRK方法沿各個(gè)方向的最大數(shù)值頻散誤差.可以看出，不同空間精度的MTSOS方法的最大頻散誤差均小于相同空間精度的SPRK方法的最大頻散誤差，而且隨著空間精度的提高，MTSOS方法對(duì)于最大頻散誤差的改進(jìn)更加明顯.

表2 不同空間精度的MTSOS方法及SPRK方法在α=0.1,Sp=0.4時(shí)的最大數(shù)值頻散誤差

3 數(shù)值算例

本節(jié)將MTSOS方法與邊界修正法、吸收邊界條件等結(jié)合，來(lái)計(jì)算三個(gè)三維彈性波模型，以驗(yàn)證MTSOS方法的有效性，以及邊界修正法和吸收邊界條件對(duì)計(jì)算區(qū)域邊界的處理效果.

3.1 均勻介質(zhì)模型

首先，選取三維均勻介質(zhì)模型，利用MTSOS方法來(lái)進(jìn)行波場(chǎng)模擬.計(jì)算區(qū)域?yàn)?0 km×40 km×40 km，介質(zhì)密度為2660 kg·m-3，P波速度為5.8 km·s-1，S波速度為3.199 km·s-1.震源位于(20 km,20 km,18 km)處.震源取為主頻2 Hz的點(diǎn)震源，震源函數(shù)為雷克子波

×exp[-8(0.6f0t-1)2],

(37)

震源位于區(qū)域中心.時(shí)間步長(zhǎng)為2.5 ms，空間步長(zhǎng)為0.2 km.計(jì)算區(qū)域的上方為自由地表邊界條件，其他五個(gè)方向采用吸收邊界條件，吸收層取為15層.采用六階空間精度的MTSOS方法進(jìn)行波場(chǎng)模擬.

在圖3中給出了在t=3.0 s時(shí)刻位移場(chǎng)的不同分量分別在震源所在的XY、XZ和YZ平面的波場(chǎng)快照.從圖中可以清晰的發(fā)現(xiàn)直達(dá)P波及直達(dá)S波，且沒有可見的數(shù)值頻散.

圖3 MTSOS方法計(jì)算得到的t=3.0 s時(shí)刻位移場(chǎng)的不同分量分別在震源所在的XY、XZ和YZ平面的波場(chǎng)快照Fig.3 Snapshots of different components of the displacement field at the time t=3.0 s calculated by the MTSOS method in the XY, XZ and YZ planes where the seismic source is located

為了驗(yàn)證MTSOS方法的準(zhǔn)確性，在(23 km,24 km,18 km)處設(shè)置一個(gè)接收器，并將接收器處的波形記錄與解析解進(jìn)行對(duì)比.圖4是在0～5 s內(nèi)位移的三分量與解析解的比較圖，其中實(shí)線表示解析解，虛線表示數(shù)值解.從圖中可以看出，接收器處的波形記錄與解析解能夠很好的吻合，說(shuō)明數(shù)值模擬結(jié)果是可信的.

圖4 六階MTSOS方法計(jì)算得到的位于(23 km,24 km,18 km)處的接收器的波形圖(a) x分量； (b) y分量； (c) z分量.Fig.4 Waveforms of the receiver located at (23 km,24 km,18 km) calculated by the sixth-order MTSOS method(a) x component; (b) y component; (c) z component.

3.2 雙層介質(zhì)模型

選取一個(gè)三維雙層介質(zhì)模型，來(lái)檢驗(yàn)MTSOS方法在利用等效介質(zhì)參數(shù)法處理模型間斷以及吸收邊界條件的有效性.計(jì)算區(qū)域?yàn)?0 km×20 km×20 km，上層介質(zhì)密度為2600 kg·m-3，P波速度為5.8 km·s-1，S波速度為3.2 km·s-1；下層介質(zhì)密度為3723 kg·m-3，P波速度為9.134 km·s-1，S波速度為4.932 km·s-1，間斷面位于z=20 km處.震源取為主頻8 Hz的點(diǎn)震源，震源函數(shù)為雷克子波.震源位于(10 km,10 km,8 km)處.時(shí)間步長(zhǎng)為1 ms，空間步長(zhǎng)為0.08 km.計(jì)算區(qū)域的上方為自由地表邊界條件，其他五個(gè)方向采用吸收邊界條件，吸收層取為15層.采用六階空間精度的MTSOS方法進(jìn)行波場(chǎng)模擬.

在圖5中給出了在t=1.75 s時(shí)刻位移場(chǎng)的不同分量分別在震源所在的XY、XZ和YZ平面的波場(chǎng)快照.從圖中可以清晰的發(fā)現(xiàn)直達(dá)波、反射波及透射波的各個(gè)震相以及產(chǎn)生的轉(zhuǎn)換波.在XY平面，由于不存在間斷面，波場(chǎng)呈同心圓結(jié)構(gòu)，不同的圓表示了以不同速度傳播的波場(chǎng).其中，最外層為直達(dá)P波(P)，向內(nèi)依次為反射P波(PP)、S波反射后轉(zhuǎn)換成的P波(SP)、直達(dá)S波(S)、P波反射后轉(zhuǎn)換成的S波(PS)、反射S波(SS)等.在XZ及YZ平面，可以看到經(jīng)間斷面得到的反射、透射和轉(zhuǎn)換波，而且沒有可見的數(shù)值頻散.圖6中選取了t=1.75 s時(shí)位移場(chǎng)的x分量在XZ平面的波場(chǎng)快照，并標(biāo)注出所有產(chǎn)生的震相.

圖5 六階MTSOS方法計(jì)算得到的t=1.75 s時(shí)刻位移場(chǎng)的不同分量分別在震源所在的XY、XZ和YZ平面的波場(chǎng)快照Fig.5 The wave field snapshots of different components of the displacement field at the time t=1.75 s calculated by the sixth-order MTSOS method in the XY, XZ, and YZ planes where the seismic source is located

圖6 六階MTSOS方法計(jì)算得到的t=1.75 s時(shí)刻位移場(chǎng)的x分量在震源所在的XZ平面的波場(chǎng)快照Fig.6 The wave field snapshots of the x component of the displacement field at the time t=1.75 s calculated by the sixth-order MTSOS method in the XZ plane where the seismic source is located

在圖7中以位移場(chǎng)的x分量在XZ平面為例給出了從t=0.5 s至t=2.5 s時(shí)刻的波場(chǎng)快照.從圖中可以看出波場(chǎng)隨時(shí)間穩(wěn)定的傳播：最初時(shí)刻僅有直達(dá)P波及直達(dá)S波.當(dāng)波傳播到間斷面，產(chǎn)生了各種反射波、透射波及轉(zhuǎn)換波.各種震相均穩(wěn)定且無(wú)可見的數(shù)值頻散.約t=1.5 s時(shí)，透射P波(Pp)傳播到邊界；約t=2.0 s時(shí)，其他波的震相也陸續(xù)傳播到邊界.其中，自由地表處的直達(dá)P波經(jīng)過自由地表的反射形成反射P波和轉(zhuǎn)換S波，而在其他的邊界處，波場(chǎng)均被吸收邊界條件很好地吸收，未產(chǎn)生可見的反射波.該數(shù)值實(shí)驗(yàn)很好地說(shuō)明等效介質(zhì)參數(shù)可以準(zhǔn)確地處理模型的間斷，從而更適用于求解非均勻介質(zhì)彈性波方程，而邊界修正法與吸收邊界條件可以很好地與數(shù)值算法結(jié)合，給出準(zhǔn)確的邊界條件.

圖7 六階MTSOS方法計(jì)算得到的位移場(chǎng)的x分量在XZ平面上從t=0.5 s至t=2.5 s時(shí)刻的波場(chǎng)快照Fig.7 The wave field snapshots of the x component of the displacement field on the XZ plane from t=0.5 s to t=2.5 s calculated by the sixth-order MTSOS method

3.3 含水裂隙介質(zhì)模型

選取一個(gè)三維含水裂隙介質(zhì)模型來(lái)考察MTSOS方法對(duì)于細(xì)微的非均勻結(jié)構(gòu)的分辨能力.含水裂隙介質(zhì)是在石油勘探領(lǐng)域經(jīng)常會(huì)接觸到的地質(zhì)結(jié)構(gòu).這是由于石油公司會(huì)向井中注水，用以驅(qū)油并減少開采時(shí)間.因此，模擬彈性波在含水裂隙介質(zhì)模型中的傳播對(duì)于石油勘探具有重要的意義.

在這個(gè)例子中，將計(jì)算區(qū)域取為2 km×2 km×2 km，背景介質(zhì)密度為2170 kg·m-3，P波速度為2.618 km·s-1，S波速度為1.421 km·s-1.震源取為主頻40 Hz的點(diǎn)震源，震源函數(shù)為雷克子波.震源位于(1 km,1 km,1 km)處.時(shí)間步長(zhǎng)為0.2 ms，空間步長(zhǎng)為5 m.計(jì)算區(qū)域的上方為自由地表邊界條件，其他五個(gè)方向采用吸收邊界條件，吸收層取為15層.在(0.9 km,0.9 km,1 km)至(1 km,1 km,1.1 km)有一個(gè)寬度為5 m的含水裂隙，裂隙的兩端恰好分別位于震源所在的XY和YZ、XZ平面內(nèi).在含水裂隙內(nèi)部，密度為1000 kg·m-3，P波速度為1.5 km·s-1，S波速度為0 km·s-1.采用六階空間精度的MTSOS方法進(jìn)行波場(chǎng)模擬.在t=0.25 s時(shí)刻位移場(chǎng)不同分量在震源所在的XY、XZ和YZ平面的波場(chǎng)快照如圖8所示.在圖8中，可以清晰地發(fā)現(xiàn)由于含水裂隙產(chǎn)生的散射波.為了更準(zhǔn)確地顯示散射波的形成過程，分別在(0.925 km,0.92 km,1 km)和(0.885 km,0.88 km,1 km)各放置一個(gè)接收器.注意到第一個(gè)接收器和震源位于含水裂隙的同側(cè)，而第二個(gè)接收器和震源分別位于含水裂隙的兩側(cè).圖9和圖10分別給出了波形圖，其中實(shí)線為不含含水裂隙時(shí)的波形圖，虛線為含含水裂隙時(shí)的波形圖.右側(cè)為波形圖的局部放大圖.可以看出，對(duì)于位于震源同側(cè)的接收器，由于含水裂隙的存在，在直達(dá)S波后可清晰地發(fā)現(xiàn)散射波所形成的一系列尾波.而對(duì)于位于震源異側(cè)的接收器，散射P波與直達(dá)S波幾乎同時(shí)到達(dá)，使得直達(dá)S波的振幅與沒有含水裂隙時(shí)的情形有了明顯的區(qū)別.該數(shù)值實(shí)驗(yàn)說(shuō)明MTSOS方法可以準(zhǔn)確模擬勘探尺度下模型內(nèi)的細(xì)微結(jié)構(gòu)變化對(duì)于波形的影響.

圖8 六階MTSOS方法計(jì)算得到的t=0.25 s時(shí)刻位移場(chǎng)的不同分量分別在震源所在的XY、XZ及YZ平面的波場(chǎng)快照Fig.8 The wave field snapshots of the different components of the displacement field at the time t=0.25 s calculated by the sixth-order MTSOS method in the XY, XZ, and YZ planes where the seismic source is located

圖9 六階MTSOS方法計(jì)算得到的位于(0.925 km,0.92 km,1 km)接收器的波形圖(a) x分量； (b) x分量局部放大圖； (c) y分量； (d) y分量局部放大圖； (e) z分量； (f) z分量局部放大圖.Fig.9 Waveforms at the receiver located at (0.925 km,0.92 km,1 km) calculated by the sixth-order MTSOS method(a) x component; (b) Partial enlarged view of x component; (c) y component; (d) Partial enlarged view of y component; (e) z component; (f) Partial enlarged view of z component.

圖10 六階MTSOS方法計(jì)算得到的位于(0.885 km,0.88 km,1 km)接收器的波形圖(a) x分量； (b) x分量局部放大圖； (c) y分量； (d) y分量局部放大圖； (e) z分量； (f) z分量局部放大圖.Fig.10 Waveforms at the receiver located at (0.885 km,0.88 km,1 km) calculated by the sixth-order MTSOS method(a) x component; (b)Partial enlarged view of x component; (c) y component; (d) Partial enlarged view of y component; (e) z component; (f) Partial enlarged view of z component.

4 結(jié)論

本文發(fā)展了適用于三維非均勻介質(zhì)彈性波方程的正演模擬新方法.該方法在時(shí)間上源自Liu等(2017)提出的修正SPRK格式，利用二級(jí)格式達(dá)到了三階時(shí)間精度；在空間上采用Ma 等(2018)提出的優(yōu)化有限差分格式，并將其推廣到三維情況，同時(shí)與等效介質(zhì)參數(shù)、吸收邊界條件等處理技術(shù)相結(jié)合，獲得了用于求解三維非均勻介質(zhì)彈性波方程的MTSOS新方法.該方法尤其適用于求解非均勻介質(zhì)彈性波方程.我們進(jìn)一步給出了MTSOS方法的穩(wěn)定性條件，并進(jìn)行了數(shù)值頻散分析.理論分析和數(shù)值結(jié)果均表明，新方法的數(shù)值頻散誤差低于相同空間精度的SPRK方法，而且隨著空間精度階數(shù)的提高，這種優(yōu)勢(shì)也越來(lái)越明顯；同時(shí)新方法具有較高的數(shù)值穩(wěn)定性，能有效提高波場(chǎng)模擬的計(jì)算效率.波場(chǎng)模擬結(jié)果進(jìn)一步驗(yàn)證了MTSOS方法能夠有效壓制數(shù)值頻散，清晰地給出彈性波傳播過程中產(chǎn)生的各種波震相，可以模擬非均勻結(jié)構(gòu)對(duì)于波場(chǎng)的影響，證實(shí)了MTSOS新方法的有效性和準(zhǔn)確性.對(duì)于含水裂隙模型，MTSOS方法能夠準(zhǔn)確地反映彈性波經(jīng)過含水裂隙所形成的散射波，接收器處接收到的波形能充分顯示含水裂隙的存在對(duì)波形的影響，表明MTSOS方法能有效模擬非均勻結(jié)構(gòu)對(duì)波場(chǎng)的影響.進(jìn)一步地，在未來(lái)的研究中，可將本研究提出的MTSOS新方法與三維全波形反演方法相結(jié)合，得到基于三維彈性波方程的新的全波形反演算法，并應(yīng)用于實(shí)際地震數(shù)據(jù)的地震成像研究中，以獲得地球內(nèi)部高分辨率的成像結(jié)果.

附錄 優(yōu)化有限差分格式二階空間偏導(dǎo)數(shù)離散格式

(A1)

(A2)

(A3)

(A4)

(A5)

(A6)

(A7)

(A8)