修正辛—準粒子法地震波數值模擬

蘇 波 庹先國 劉知貴

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

1 引言

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

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

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

2 空間離散—準粒子體系

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

2.1 Hamilton體系的線性形式

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

圖1 準粒子體系編號規則

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

(1)

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

(2)

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

(3)

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

2.2 粒子間的相互作用系數

由于彈性波是矢量波，在各個方向上粒子間的相互作用關系不同，但在相同軸上對中心準粒子的作用力大小相同。準粒子間的相互作用系數為

(4)

(5)

(6)

令二維網格粒子間距為d，根據連續偏導數與離散量的關系，由式(6)可得

(7)

(8)

比較。式中：u為位移；λ、μ為拉梅常數，滿足以下關系

(9)

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

(10)

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

3 時間離散—修正辛格式

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

3.1 構造修正辛格式

(11)

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

(12)

其中

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

(13)

(14)

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

(15)

3.2 穩定性分析

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

(16)

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

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

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

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

(17)

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

|≤2?

-5.308≤Δt2ξi≤0

(18)

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

(19)

3.3 頻散分析

數值頻散是影響數值模擬結果的重要因素之一[28]。為了定量分析頻散效應，定義數值頻散比率

和采樣率

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

(20)

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

圖2為數值頻散誤差曲線

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

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

穩定性參數方 法M1M2M3|Δt2ξi|max5.30847.107a0.44760.44460.4429b0.8380.72740.9695

4 數值實驗

4.1 拉梅問題

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

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

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

4.2 Sigsbee 2B速度模型

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

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

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

圖4 Sigsbee 2B速度模型

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

圖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 結束語

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

[1] 曹禹,楊孔慶.對聲波和彈性波傳播模擬的Hamilton系統方法.物理學報,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] 馬德堂,朱光明.有限元法與偽譜法混合求解彈性波動方程.地球科學與環境學報,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] 汪文帥.基于辛—譜元方法的地震波場模擬研究[學位論文].北京：中國科學院大學,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] 汪勇,段焱文,王婷等.優化近似解析離散化方法的二維彈性波波場分離模擬.石油地球物理勘探,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] 羅明秋,劉洪,李幼銘.地震波傳播的哈密頓表述及辛幾何算法.地球物理學報,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] 馮康.哈密爾頓系統的辛幾何算法.浙江杭州：浙江科技出版社,2003.

[24] 馬毅恒,井西利,邢小寧等.彈性波的三維Hamilton系統方法模擬.石油地球物理勘探,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格式.地球物理學報,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] 梁文全,王彥飛,楊長春.聲波方程數值模擬矩形網格有限差分系數確定法.石油地球物理勘探,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.