修正辛—準粒子法地震波數(shù)值模擬

蘇 波 庹先國 劉知貴

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

1 引言

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

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

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

2 空間離散—準粒子體系

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

2.1 Hamilton體系的線性形式

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

圖1 準粒子體系編號規(guī)則

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

(1)

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

(2)

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

(3)

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

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

由于彈性波是矢量波，在各個方向上粒子間的相互作用關(guā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)可以得到一組準粒子體系相互作用系數(shù)與P波速度vP、介質(zhì)密度ρ以及泊松比c0之間的相互關(guān)系[24]

(10)

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

3 時間離散—修正辛格式

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

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

(11)

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

(12)

其中

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

(13)

(14)

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

(15)

3.2 穩(wěn)定性分析

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

(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

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

(17)

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

|≤2?

-5.308≤Δt2ξi≤0

(18)

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

(19)

3.3 頻散分析

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

和采樣率

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

(20)

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

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

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

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

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

4 數(shù)值實驗

4.1 拉梅問題

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

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

表2 三種方法的計算效率和精度對比

4.2 Sigsbee 2B速度模型

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

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

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

圖4 Sigsbee 2B速度模型

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

圖5 不同方法得到的Sigsbee 2B模型的水平位移分量波場快照 (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 檢波點(3800m，1400m)處的水平位移分量記錄對比圖

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

在圖a、圖b中M3標注為參考解

5 結(jié)束語

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

[1] 曹禹,楊孔慶.對聲波和彈性波傳播模擬的Hamilton系統(tǒng)方法.物理學(xué)報,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] 馬德堂,朱光明.有限元法與偽譜法混合求解彈性波動方程.地球科學(xué)與環(huán)境學(xué)報,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] 唐懷谷,何兵壽.一階聲波方程時間四階精度差分格式的偽譜法求解.石油地球物理勘探,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] 汪文帥.基于辛—譜元方法的地震波場模擬研究[學(xué)位論文].北京：中國科學(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)化近似解析離散化方法的二維彈性波波場分離模擬.石油地球物理勘探,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é)報,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é)報,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] 梁文全,王彥飛,楊長春.聲波方程數(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.