蘇 波 庹先國(guó) 劉知貴

(①中國(guó)工程物理研究院研究生院，四川綿陽(yáng) 621900； ②四川理工學(xué)院，四川自貢 643002；③西南科技大學(xué)計(jì)算機(jī)科學(xué)與技術(shù)學(xué)院，四川綿陽(yáng) 621010)

1 引言

地震波波動(dòng)方程的直接數(shù)值求解方法需要對(duì)空間和時(shí)間進(jìn)行數(shù)值離散[1]。針對(duì)空間離散，常用的方法有有限差分法(FDM)[2,3]、有限元法(FEM)[4-6]、偽譜法(PSM)[7-9]、譜元法(SEM)[10,11]和近似解析離散法(NADM)[12-14]等。近年來(lái)，基于粒子的方法作為傳統(tǒng)基于連續(xù)波動(dòng)方程求解方法的一個(gè)替代選項(xiàng)而被用于地震波數(shù)值模擬。Takekawa等[15]將Suzuki等[16]提出的Hamiltonian particle method(HPM)用于地震波數(shù)值模擬。

相比以上在空間離散方面的大量研究，對(duì)時(shí)間離散方面的研究相對(duì)較少[17]。時(shí)間離散在顯式格式上大多采用二階中心差分(second-order central difference，CD)，然而該算法在時(shí)間步長(zhǎng)較大時(shí)將導(dǎo)致嚴(yán)重的數(shù)值頻散，且隨著大規(guī)模運(yùn)算的進(jìn)行，數(shù)值頻散效應(yīng)造成累積誤差逐漸增大?？紤]到隱式格式在不考慮庫(kù)朗數(shù)的情況下可以取得相對(duì)較大的時(shí)間步長(zhǎng)，羅明秋等[18]對(duì)隱式辛格式做了相關(guān)研究，但隱式算法需要求解方程組，計(jì)算量較大。Chen[19]討論了三種高階時(shí)間離散策略(Lax-Wendroff方法、三級(jí)四階辛格式RKN方法和三階分步Runge-kutta方法)，大量數(shù)值實(shí)驗(yàn)表明，辛算法對(duì)時(shí)間誤差壓制能力強(qiáng)于Lax-Wendroff方法。實(shí)際計(jì)算表明，二階格式精度往往不夠高[11,20]，雖然采用傳統(tǒng)三階、四階辛格式提高了時(shí)間精度，但是它們具有負(fù)的辛系數(shù)，這與時(shí)間演進(jìn)方向不符，正的辛系數(shù)才符合實(shí)際情況，更具有長(zhǎng)時(shí)間穩(wěn)定性[21]。

本文首先基于Hamilton力學(xué)，并借鑒分子動(dòng)力學(xué)的基本原理，介紹一種簡(jiǎn)單、直觀、物理意義明確的空間準(zhǔn)粒子離散體系。然后利用Liu等[17,22]提出的構(gòu)造高階辛算法策略，也即在傳統(tǒng)辛格式中額外加入空間算子，構(gòu)造出一種新的具有時(shí)間三階精度的辛格式，達(dá)到修正精度的目的。最后，為了測(cè)試新的辛格式和空間準(zhǔn)粒子體系的準(zhǔn)確性，用均勻模型的解析解與本文方法的數(shù)值解進(jìn)行對(duì)比。為了測(cè)試穩(wěn)定性，進(jìn)一步選取Sigsbee 2B速度模型進(jìn)行驗(yàn)證。

2 空間離散—準(zhǔn)粒子體系

地震波的傳播過(guò)程是能量逐步耗散殆盡的過(guò)程，但在實(shí)際應(yīng)用中常常將其假設(shè)為無(wú)能量損耗過(guò)程[18]。馮康[23]指出，所有保守?zé)o耗散系統(tǒng)均可以表示成Hamilton形式。本文利用彈性波方程描述地震波傳播現(xiàn)象，為了討論方便，只考慮各向同性彈性介質(zhì)。

2.1 Hamilton體系的線性形式

針對(duì)Hamilton體系的動(dòng)力演化特性，由分子動(dòng)力學(xué)及能量守恒法則，可將系統(tǒng)看成由N個(gè)粒子組成(圖1)。

圖1 準(zhǔn)粒子體系編號(hào)規(guī)則

假設(shè)每個(gè)粒子的質(zhì)量為mi，系統(tǒng)總能量由動(dòng)能和勢(shì)能構(gòu)成。動(dòng)能由各個(gè)質(zhì)點(diǎn)i在x、y方向上對(duì)應(yīng)的動(dòng)量pxi、pyi決定；勢(shì)能E由分子離開(kāi)平衡位置的位移量qxi、qyi決定，Hamilton量可表示為

(1)

為簡(jiǎn)化起見(jiàn)，在該體系中各個(gè)粒子僅與上、下、左、右和4個(gè)對(duì)角粒子發(fā)生相互作用，且認(rèn)為其作用力和相對(duì)位移呈近似線性關(guān)系。假設(shè)微觀粒子具有不可辨性,體系勢(shì)能在平衡位置的Taylor展開(kāi)為

(2)

式中粒子在平衡位置處的勢(shì)能E0為常數(shù)，E的一階導(dǎo)數(shù)在平衡位置為零，而其二階偏導(dǎo)數(shù)正定。令E0=0，由式(1)可以進(jìn)一步得到Hamilton體系的線性形式為

(3)

式中αij、βij和γij分別表示在x、y方向的粒子間相互作用系數(shù)。

2.2 粒子間的相互作用系數(shù)

由于彈性波是矢量波，在各個(gè)方向上粒子間的相互作用關(guān)系不同，但在相同軸上對(duì)中心準(zhǔn)粒子的作用力大小相同。準(zhǔn)粒子間的相互作用系數(shù)為

(4)

(5)

(6)

令二維網(wǎng)格粒子間距為d，根據(jù)連續(xù)偏導(dǎo)數(shù)與離散量的關(guān)系，由式(6)可得

(7)

(8)

比較。式中：u為位移；λ、μ為拉梅常數(shù)，滿足以下關(guān)系

(9)

式中c0為泊松比。由式(7)～式(9)可以得到一組準(zhǔn)粒子體系相互作用系數(shù)與P波速度vP、介質(zhì)密度ρ以及泊松比c0之間的相互關(guān)系[24]

(10)

將式(10)作為空間離散系數(shù)代入式(5)，實(shí)現(xiàn)對(duì)波動(dòng)方程的空間離散。

3 時(shí)間離散—修正辛格式

文中通過(guò)在二階辛格式[25]的基礎(chǔ)上構(gòu)建一個(gè)新的修正辛格式，該格式具有三階時(shí)間精度，并在理論上推導(dǎo)了其數(shù)值頻散和穩(wěn)定性。為了方便討論，記修正辛格式為M1、二級(jí)二階辛格式[21]為M2、三級(jí)三階辛格式[26]為M3。

3.1 構(gòu)造修正辛格式

(11)

式(11)進(jìn)一步變換為Hamilton形式

(12)

其中

為了求解半離散方程組(式(12))，常用二階方法[23]

(13)

(14)

對(duì)于式(14)可以基于Taylor原理得到三階精度條件

(15)

3.2 穩(wěn)定性分析

為了驗(yàn)證修正辛格式的穩(wěn)定性，設(shè)在t=nΔt時(shí)刻式(8)的簡(jiǎn)諧解為

(16)

式中：ωnum為角頻率；kx=kcosφ和ky=ksinφ分別為x,y方向的波數(shù),φ為波數(shù)矢量k=(kx,ky)與x軸間的夾角。

將式(16)代入式(14)得到形如

(υn+1,un+1)T=G(υn,un)T

的時(shí)間演進(jìn)方程。式中：G為增長(zhǎng)矩陣。設(shè)G和空間算子L的特征值分別為ψ和ξ，可得G的特征多項(xiàng)式為

(17)

使式(14)的修正辛格式穩(wěn)定的充要條件為|ψ|≤1[27]。由式(17)易得

|≤2?

-5.308≤Δt2ξi≤0

(18)

由于得到穩(wěn)定性解析表達(dá)式較困難，借由Liu等[17]的討論可以得到

(19)

3.3 頻散分析

數(shù)值頻散是影響數(shù)值模擬結(jié)果的重要因素之一[28]。為了定量分析頻散效應(yīng)，定義數(shù)值頻散比率

和采樣率

式中λ為波長(zhǎng)。則頻散關(guān)系可進(jìn)一步表示為[27]

(20)

式中ξS為由式(17)得到的特征值的最小絕對(duì)值，其對(duì)應(yīng)于S波的特征值。

圖2為數(shù)值頻散誤差曲線

(|R-1|cosφ, |R-1|sinφ)

在(x1,x2)平面上的投影，其中φ為波傳播方向與x1軸的夾角，α=0.5。由圖可見(jiàn)：在兩種采樣率的情況下，M1與M3的頻散誤差曲線交織在一起，幾乎相同；在采用較小采樣率時(shí)M1和M3的誤差均小于M2(圖2a)；當(dāng)采用較大采樣率時(shí)，在軸向上M2的誤差略小于另外兩種格式(圖2b)。

穩(wěn)定性參數(shù)方 法M1M2M3|Δt2ξi|max5.30847.107a0.44760.44460.4429b0.8380.72740.9695

4 數(shù)值實(shí)驗(yàn)

4.1 拉梅問(wèn)題

利用拉梅問(wèn)題驗(yàn)證本文方法對(duì)彈性波數(shù)值模擬的計(jì)算精度和效率。圖3為數(shù)值解與解析解[29]在水平與垂直方向的波形對(duì)比圖。由圖可見(jiàn)：由本文方法所得的數(shù)值解與解析解擬合較好(圖3a、圖3b)；局部放大圖顯示M1與M3相互疊加在一起，且與解析解擬合較好，而M2與解析解擬合略差，說(shuō)明M1與M3的精度相當(dāng)，且高于M2(圖3c、圖3d)。表2為三種方法的計(jì)算效率和精度對(duì)比。由表可見(jiàn)：M2用時(shí)最少(159.428s)，但是水平和垂直方向的誤差總量較M1與M3偏大,與前文對(duì)二階辛格式的分析認(rèn)識(shí)一致；M1占用的計(jì)算內(nèi)存最少，雖然其計(jì)算時(shí)間比M2多，但在計(jì)算時(shí)間少的情況下幾乎與M3的計(jì)算精度相同。

圖3 數(shù)值解與解析解在水平與垂直方向的波形對(duì)比圖 (a)水平分量； (b)垂直分量； (c)圖a的局部放大； (d)圖b的局部放大

表2 三種方法的計(jì)算效率和精度對(duì)比

4.2 Sigsbee 2B速度模型

為了檢驗(yàn)修正辛—準(zhǔn)粒子法在非均勻介質(zhì)模型中彈性波計(jì)算的穩(wěn)定性，選取Sigsbee 2B速度模型(圖4)進(jìn)行測(cè)試。模擬結(jié)果如圖5和圖6所示。

圖5為不同方法得到的Sigsbee 2B模型的水平位移分量波場(chǎng)快照。由圖可見(jiàn)，當(dāng)M1的時(shí)間步長(zhǎng)大于M2、M3時(shí)，三種方法的波場(chǎng)模擬結(jié)果相似。

由于在非均勻介質(zhì)模型中沒(méi)有解析解，考慮到較小的時(shí)間步長(zhǎng)可以取得更精確的結(jié)果，因此選擇具有高階時(shí)間精度的M3作為參考解，且令其一直取較小的時(shí)間步長(zhǎng)不變。圖6為檢波點(diǎn)(3800m，1400m)處的水平位移分量記錄對(duì)比圖。由圖可見(jiàn)，M1與M2取兩種時(shí)間步長(zhǎng)的結(jié)果與參考解相近，相互疊加在一起(圖6a、圖6b)，得到的M1、M2在時(shí)間步長(zhǎng)為1、2ms時(shí)(與參考解)的誤差曲線(圖6c～圖6f)表明，M1不論取較小的時(shí)間步長(zhǎng)還是較大的時(shí)間步長(zhǎng)，其誤差均小于M2。數(shù)值結(jié)果證實(shí)了修正—準(zhǔn)粒子方法在模擬彈性波方面的優(yōu)越能力。

圖4 Sigsbee 2B速度模型

模型速度變化范圍為1437～4511m/s，選取模型的400×150個(gè)速度網(wǎng)格點(diǎn)，網(wǎng)格間距為20m。爆炸源為主頻10Hz的Ricker子波,位于模型正中心，PML吸收層為20個(gè)網(wǎng)格厚度

圖5 不同方法得到的Sigsbee 2B模型的水平位移分量波場(chǎng)快照 (a)M1(Δt=2ms,t=0.6s)； (b)M1(Δt=2ms,t=0.8s)； (c)M2(Δt=1ms,t=0.6s)； (d)M2(Δt=1ms,t=0.8s)； (e)M3(Δt=1ms,t=0.6s)； (f)M3(Δt=1ms,t=0.8s)

圖6 檢波點(diǎn)(3800m，1400m)處的水平位移分量記錄對(duì)比圖

(a)M1(Δt=1ms)、M2(Δt=1ms)、M3(Δt=1ms)的水平位移分量記錄對(duì)比；(b)M1(Δt=2ms)、M2(Δt=2ms)、M3(Δt=1ms)的水平位移分量記錄對(duì)比；(c)M1(Δt=1ms)誤差曲線；(d)M1(Δt=2ms)誤差曲線；(e)M2(Δt=1ms)誤差曲線； (f)M2(Δt=2ms)誤差曲線

在圖a、圖b中M3標(biāo)注為參考解

5 結(jié)束語(yǔ)

時(shí)間離散方面，本文在二階辛格式基礎(chǔ)給出一種新的修正辛格式，并與空間準(zhǔn)粒子離散方法相結(jié)合得到了地震波模擬的修正辛—準(zhǔn)粒子法。結(jié)合理論推導(dǎo)與數(shù)值算例發(fā)現(xiàn)，修正辛—準(zhǔn)粒子法在數(shù)值頻散壓制和數(shù)值穩(wěn)定性提升等方面具有明顯的優(yōu)勢(shì)。尚需指出，在空間離散方面，修正辛—準(zhǔn)粒子法所采用的離散點(diǎn)較少，只考慮了粒子與其周?chē)藗€(gè)點(diǎn)的相互作用，需研究引入更多粒子參與計(jì)算的情況，并根據(jù)修正策略構(gòu)造其他修正辛算法，并將算法與其他空間離散法相結(jié)合做進(jìn)一步研究。

[1] 曹禹,楊孔慶.對(duì)聲波和彈性波傳播模擬的Hamilton系統(tǒng)方法.物理學(xué)報(bào),2003,52(8):1984-1992. Cao Yu,Yang Kongqing.Hamiltonian system approach for simulation of acoustic and elastic wave propagation.Acta Physica Sinica,2003,52(8):1984-1992.

[2] Alford R,Kelly K,Boore D M.Accuracy of finite-diffe-rence.Geophysics,1974,39(6):834-842.

[3] Virieux J.P-SV wave propagation in heterogeneous media:velocity-stress finite-difference method.Geophysics,1986,51(4):889-901.

[4] Lysmer J,Drake L A.A finite element method for seismology.Methods in Computational Physics Advances in Research & Applications,1972,11:181-216.

[5] Moser F,Jacobs L J,Qu J.Modeling elastic wave propagation in waveguides with the finite element method.Ndt & E International,1999,32(4):225-234.

[6] Liu S,Li X,Wang W et al.A mixed-grid finite element method with PML absorbing boundary conditions for seismic wave modelling.Journal of Geophysics and Engineering,2014,11(5):055009.

[7] Fornberg B.High-order finite difference and the pseudospectral method on staggered grids.Siam Journal on Numerical Analysis,1990,27(4):904-918.

[8] 馬德堂,朱光明.有限元法與偽譜法混合求解彈性波動(dòng)方程.地球科學(xué)與環(huán)境學(xué)報(bào),2004,26(1):61-64. Ma Detang,Zhu Guangming.Hybrid method combining finite difference and pseudospectral method for solving the elastic wave equation.Journal of Earth Sciences and Environment,2004,26(1):61-64.

[9] 唐懷谷,何兵壽.一階聲波方程時(shí)間四階精度差分格式的偽譜法求解.石油地球物理勘探,2017,52(1):71-80. Tang Huaigu,He Bingshou.Pseudo spectrum method of first-order acoustic wave equation finite-difference schemes with fourth-order time difference accuracy.OGP,2017,52(1):71-80.

[10] Komatitsch D,Tromp J.The spectral-element me-thod,Beowulf computing,and global seismology.Science,2002,298(5599):1737-1742.

[11] 汪文帥.基于辛—譜元方法的地震波場(chǎng)模擬研究[學(xué)位論文].北京：中國(guó)科學(xué)院大學(xué),2013.

[12] Yang D H,Teng J,Zhang Z et al.A nearly analytic discrete method for acoustic and elastic wave equations in anisotropic media.BSSA,2003,93(2):882-890.

[13] Chen Y,Song G,Xue Z et al.An improved optimal nearly analytical discretized method for 2D scalar wave equation in heterogeneous media based on the modified nearly analytical discrete operator.Geophy-sics,2014,79(6):T349-T362.

[14] 汪勇,段焱文,王婷等.優(yōu)化近似解析離散化方法的二維彈性波波場(chǎng)分離模擬.石油地球物理勘探,2017,52(3):458-467. Wang Yong,Duan Yanwen,Wang Ting et al.Numerical simulation of elactic wave separation in 2D isotropic medium with the optimal nearly-analytic discretization.OGP,2017,52(3):458-467.

[15] Takekawa J,Madariage R,Mikada H et al.Numerical simulation of seismic wave propagation produced by earthquake by using a particle method.Geophysical Journal International,2012,191(3):1305-1316.

[16] Suzuki Y,Koshizuka S,Oka Y.Hamiltonian moving-particle semi-implicit (HMPS) method for incompressible fluid flows.Computer Methods in Applied Mechanics and Engineering,2007,196(29-30):2876-2894.

[17] Liu S,Yang D,Ma J.A modified symplectic PRK scheme for seismic wave modeling.Computers & Geosciences,2017,99:28-36.

[18] 羅明秋,劉洪,李幼銘.地震波傳播的哈密頓表述及辛幾何算法.地球物理學(xué)報(bào),2001,44(1):120-128. Luo Mingqiu,Liu Hong,Li Youming.Hamiltonian description and symplectic method of seismic wave propagation.Chinese Journal of Geophysics,2001,44(1):120-128.

[19] Chen J B.High-order time discretizations in seismic modeling.Geophysics,2007,72(5):SM115-SM122.

[20] Liu S,Li X,Wang W et al.A new kind of optimal se-cond-order symplectic scheme for seismic wave simulations.Science in China:Earth Sciences,2014,57(4):751-758.

[21] Liu S,Yang D,Lang C et al.Modified symplectic schemes with nearly-analytic discrete operators for acoustic wave simulations.Computer Physics Communications,2017,213:52-63.

[22] Liu S,Li X,Wang W et al.A modified symplectic scheme for seismic wave modeling.Journal of Applied Geophysics,2015,116(5):110-120.

[23] 馮康.哈密爾頓系統(tǒng)的辛幾何算法.浙江杭州：浙江科技出版社,2003.

[24] 馬毅恒,井西利,邢小寧等.彈性波的三維Hamilton系統(tǒng)方法模擬.石油地球物理勘探,2012,47(4):524-531. Ma Yiheng,Jing Xili,Xing Xiaoning et al.3D Hamiltonian system approach simulation for elastic wave.OGP,2012,47(4):524-531.

[25] 劉少林,李小凡,汪文帥等.求解彈性波方程的辛RKN格式.地球物理學(xué)報(bào),2015,58(4):1355-1366. Liu Shaolin,Li Xiaofan,Wang Wenshuai.A symplectic RKN scheme for solving elastic wave equations.Chinese Journal of Geophysics,2015,58(4):1355-1366.

[26] Li X,Li Y,Zhang M et al.Scalar seismic-wave equation modeling by a multisymplectic discrete singular convolution differentiator method.BSSA,2011,101(4):1710-1718.

[27] Basabe J D D,Sen M K.Grid dispersion and stability criteria of some common finite-element methods for acoustic and elastic wave equations.Geophysics,2007,72(6):T81-T95.

[28] 梁文全,王彥飛,楊長(zhǎng)春.聲波方程數(shù)值模擬矩形網(wǎng)格有限差分系數(shù)確定法.石油地球物理勘探,2017,52(1):56-62. Liang Wenquan,Wang Yanfei,Yang Changchun.Acoustic wave equation modeling with rectangle grid finite difference operator and its linear time space domain solution.OGP,2017,52(1):56-62.

[29] De Hoop A.A modification of Cagniard's method for solving seismic pulse problems.Applied Scientific Research,Section B,1960,8(1):349-356.