劉少林，李小凡，劉有山，陳世仲

1 中國科學(xué)院地質(zhì)與地球物理研究所，中國科學(xué)院地球深部研究重點(diǎn)實(shí)驗(yàn)室，北京 100029

2 中國科學(xué)院大學(xué)，北京 100049

1 引 言

地震波正演模擬技術(shù)一直是地震學(xué)研究的熱點(diǎn)，無論是中小尺度的地震波逆時(shí)偏移［1－2］、全波形反演［3－4］，還是區(qū)域性地震動研究［5－6］以及全球地震波模擬［7－9］和地球介質(zhì)非均勻性反演等方面，正演模擬都起著關(guān)鍵作用.經(jīng)過幾十年的發(fā)展，地震波模擬方法已逐漸完善，總體而言，這些方法可以分為三大類，即射線法、積分方程法以及波動方程直接解法，Carcione等［10］和 Yang等［11］對這些方法做了詳細(xì)的回顧與討論.本文主要討論聲波方程的直接解法.

無論是經(jīng)典的有限差分法［12－14］，還是偽譜法［15－16］、有限元法［4，6］、譜元法［17－18］等眾多方法，在時(shí)間離散上常常使用二階中心差分或Newmark格式，雖然計(jì)算效率較高，但低階的時(shí)間近似與高階的空間離散格式往往不匹配，以致影響了地震波模擬的精度.針對時(shí)間精度較低等問題，Dablain［19］將Lax等［20］提出的Lax－Wendroff方法運(yùn)用至地震波模擬之中，其思想是運(yùn)用高次空間微分修正時(shí)間精度.為了計(jì)算方便，一般得到時(shí)間四階精度.但是直接運(yùn)用有限長度的網(wǎng)格點(diǎn)波場值修正時(shí)間高次項(xiàng)，精度往往不夠，針對此問題，Tong等［21］利用Yang等［22］發(fā)展的 Nearly Anayltic Discrete Method（NAD）近似空間高階微分，得到了時(shí)間四階、空間八階的兩步Stereo－Modeling Method（STEM）地震波模擬方法，兩步STEM法在實(shí)際地震波模擬中精度的提高、數(shù)值頻散壓制和數(shù)值方向各向異性控制等方面均明顯優(yōu)于Lax－Wendroff方法.

Chen在研究Lax－wendroff方法之后發(fā)現(xiàn)該方法在實(shí)際運(yùn)用中誤差增長較快［23］，針對此缺點(diǎn)，Chen將高階的辛 Runge－Kutta－Nystr?m（RKN）與辛Partitioned－Runge－Kutta（PRK）格式運(yùn)用至聲波模擬之中［24－25］，因辛格式嚴(yán)格保持 Hamiltion系統(tǒng)微分二形式不變的特性［26］，有效地解決了地震波長時(shí)間計(jì)算中的誤差累積等問題.Li等［27－28］在時(shí)間離散上采用三級三階的辛PRK格式、空間離散上采用褶積微分算子法，形成了一套較為高效的地震波模擬方法，這些格式都是顯式辛格式.與顯式辛格式不同的隱式辛格式，其優(yōu)點(diǎn)為無條件穩(wěn)定性，Luo等［29－30］主要討論了二階隱式辛格式，但隱式辛格式涉及到微分算子的求逆運(yùn)算，直接LU分解消耗大量的時(shí)間，譜因式分解雖計(jì)算速度快，但面臨著精度損失、邊界處誤差過大等弊端，改進(jìn)的混合算法雖然精度和計(jì)算效率上有所提高，但仍無法滿足高精度與高效率地震波模擬的要求.Ma等［31］在分析總結(jié)前人工作基礎(chǔ)上，提出在時(shí)間離散上采用二階的辛PRK格式［32］、空間離散上采用NAD算子，得到了一種高效、高精度的地震波模擬方法.

本文將偽譜法離散后的半離散聲波方程變換至Hamilton系統(tǒng)，在分析不同時(shí)間積分格式面臨的諸多問題之后，提出在時(shí)間離散上使用二級的辛RKN格式，以保證較高的計(jì)算效率，運(yùn)用根數(shù)理論得到辛條件方程組.針對兩個(gè)自由度的方程組，為了實(shí)際模擬需要，從精度、頻散、穩(wěn)定性三方面優(yōu)化系數(shù)，最后得到了三組優(yōu)化系數(shù)，理論上分析了新得到的優(yōu)化格式在求解聲波方程時(shí)的優(yōu)良特性.在實(shí)際運(yùn)用之中，可以根據(jù)不同的實(shí)際需要選擇不同的格式以得到更好的模擬結(jié)果.

2 辛RKN格式

2.1 聲波方程的偽譜法離散

在地震波模擬與成像中有限差分算子得到了廣泛的應(yīng)用［19］，但有限長度的有限差分算子精度較低、數(shù)值頻散明顯［11］.偽譜法在精度和數(shù)值頻散壓制方面明顯優(yōu)于有限差分法［33］，隨著計(jì)算機(jī)硬件的發(fā)展，偽譜法有望運(yùn)用至三維高精度地震波模擬與成像之中［16］.

考慮如下形式的地震聲波方程：

其中，u（x，y，z，t）為聲波波場值，c（x，y，z）為聲速.在空間上采用均勻網(wǎng)格對（1）式離散，即uijk≈u（iΔx，jΔy，kΔz），i＝1，2，…，Ni，j＝1，2，…，Nj，k＝1，2，…，Nk，Δx，Δy，Δz為空間網(wǎng)格間距.若采用偽譜法計(jì)算（1）式中的空間微分，得到的半離散方程如下：表示空間節(jié)點(diǎn)的離散波數(shù)值，記f（u）＝c2F－1［w2·F（u）］.

其中，u為離散后的波場向量，F(xiàn)，F(xiàn)－1分別表示Fourier正變換和逆變換，w＝ ［w1，1，1，…，wNxNyNz］，

2.2 辛RKN格式及條件方程

對于（2）式可以采用多種時(shí)間離散格式，如三階的 Runge－Kutta方 法 （RK）［34］、Newmark 格 式［5，7，8］、二階辛 PRK 格式［31－32］、三階辛 PRK 格式［23，27－28，34］、四階的辛RKN格式［24－26］以及最近Yang等發(fā)展的顯式化后的 RK格式與 Adams格式［35－36］.這些方法中二階的辛PRK格式效率最高，每一個(gè)時(shí)間步只需兩次正反Fourier變換，高階格式雖然有較高的精度，但每個(gè)時(shí)間步計(jì)算量過大以致嚴(yán)重影響計(jì)算效率.在大尺度以及全球尺度地震波模擬時(shí)，由于空間微分計(jì)算量過大，每個(gè)時(shí)間步無法承受過多中間步，所以本文使用二階辛RKN格式對（2）式進(jìn)行離散.

引入變元v＝u·，則方程（2）可降階為

考慮（3）式為 Hamilton系統(tǒng)情形，即（3）式可表示為如下形式：

其中H（u，v）為一標(biāo)量函數(shù)，由（3）、（4）式知H滿足如下等式：

因此H滿足：當(dāng)f（u）為某一標(biāo)量函數(shù)V的梯度時(shí)才可考慮（4）式的辛算法.對于一般的s級RKN格式可以寫為［26］：

其中，h為時(shí)間步長，ci，aij，ˉbj，bj為辛系數(shù).可以證明，如果（6）式是顯式辛的，（6）式應(yīng)滿足如下方程［26，37］：

在地震波模擬時(shí)，希望格式（6）的計(jì)算效率與二階的PRK格式相近［31－32］，當(dāng)s＝2時(shí)即滿足高效地震波模擬的要求.在精度上希望格式（6）精度盡可能提高，對于非約化的辛RKN格式［26，37］（即格式中無冗余級，bj≠0），滿足（7）式的格式（6）是p階的，可以用圖形表示微分關(guān)系的根數(shù)理論［37－38］，即任何s－樹sτ，都有一個(gè)有根的s－樹ρsτ由sτ的一胖節(jié)點(diǎn)提升得到，ρsτ的權(quán)與密度滿足如下等式：

其中φ為ρsτ的權(quán)，γ為ρsτ的密度.如果二級的辛RKN格式是三階的，有4個(gè)非多余有根s－樹需滿足（8）式，將辛條件（7）帶入（8）式化簡得到如下等式：

2.3 一種誤差最小辛RKN格式

求解（9）式發(fā)現(xiàn)，方程無實(shí)數(shù)解，即二級的辛RKN格式無法使精度達(dá)到三階，我們可以根據(jù)截?cái)嗾`差最小原理得到精度接近三階的誤差最小辛RKN格式，即（9）式中的前兩式滿足的同時(shí)，構(gòu)造目標(biāo)函數(shù)，使目標(biāo)函數(shù)盡量最小，

運(yùn)用非線性優(yōu)化方法，可以得到如下一種誤差最小辛RKN格式（記為M1）：

2.4 一種強(qiáng)穩(wěn)定的辛RKN格式

運(yùn)用二級辛RKN格式進(jìn)行地震波模擬時(shí)，時(shí)間步長與空間步長應(yīng)滿足一定的比例關(guān)系，否則計(jì)算失穩(wěn)以至無法進(jìn)行.為了使分析過程中量綱一致，引入變量~v＝vh，考慮如下的簡諧解，

由（6）式得到的時(shí)間演進(jìn)方程為

其中，θ＝ckh，e1＝b1c21＋b1c22，e2＝b1b2（c2－c1）.由（13）式知時(shí)間演進(jìn)方程的特征值為

其中

選擇不同的c1，c2，滿足（9）、（15）式的同時(shí)力圖使θ取到最大值.最后得到的一種強(qiáng)穩(wěn)定的辛RKN格式（記為 M2）：

2.5 一種最小頻散辛RKN格式

由于用時(shí)間和空間的離散點(diǎn)值近似替代時(shí)間空間連續(xù)變化的函數(shù)，必然導(dǎo)致真實(shí)信號中的部分信號無法拾取而導(dǎo)致數(shù)值頻散.本文用偽譜法計(jì)算空間微分，在Nyquist波數(shù)范圍內(nèi)偽譜法對波數(shù)完全覆蓋，用格式（6）進(jìn)行地震波模擬時(shí)，數(shù)值頻散將主要來源于時(shí)間離散.由（13）、（14）式可知，（13）式左端的exp（iωh）應(yīng)與λ相等，由此可以得到相速度c與真實(shí)速度c之比滿足如下等式：

其中，θ又可以表示為θ＝rkΔ，庫朗數(shù)

選擇適當(dāng)?shù)膸炖蕯?shù)（如r＝1），當(dāng)kΔ由0變化至π，選擇不同的系數(shù)使得相速度與真實(shí)速度最接近，為此構(gòu)造如下目標(biāo)函數(shù)：

用非線性最優(yōu)化方法得到一種最小頻散辛RKN格式（記為 M3），

3 三種辛RKN格式特性分析

第2節(jié)從精度、穩(wěn)定性和數(shù)值頻散優(yōu)化三方面得到了3種辛RKN格式，本節(jié)就這三方面與常見格式對比，以凸顯本文格式的特性.對比選用的格式包括：二階辛 PRK 格式［31－32］（記為 M4）、二階優(yōu)化PRK格式［39］（記為 M5）和三級三階PRK格式［27－28］（記為 M6）.

3.1 精度分析

將（10）式中的E1視為截?cái)嗾`差，對比M1—M6 6種格式的E1.對比結(jié)果如表1所示.就精度而言，三階格式M6精度最高；M1精度其次，其精度逼近三階格式M6；M2與M5精度接近，M3和M4精度最差.

表1 6種辛格式截?cái)嗾`差對比Table 1 Error truncation of six symplectic schemes

3.2 穩(wěn)定性分析

按照2.4節(jié)分析方法，可以得到M1—M6的穩(wěn)定性條件.表2給出了6種格式從一維到三維的穩(wěn)定性極限（即庫朗數(shù)r所取的極大值）.從表2可以看出，理論上M3穩(wěn)定性最好.在常見格式中三階格式M6穩(wěn)定性最好，二階的M4格式穩(wěn)定性最差，優(yōu)化的二階格式M5介于兩者之間，整體而言常見格式穩(wěn)定性要遜色于本文3種格式.

表2 6種格式一維到三維穩(wěn)定性比較Table 2 Stability limits of six symplectic schemes for 1D，2D，and 3Dscalar wave simulations

3.3 頻散分析

按照2.3節(jié)的分析方法，可以得到M1—M6的頻散關(guān)系，圖1（a、b）為r＝0.5和r＝0.8時(shí)6種格式1D情形下的頻散曲線.由圖1a可知，r取較小值時(shí)M6頻散最小，M3頻散其次，M4與M5頻散較大；在圖1b中當(dāng)r＝0.8時(shí)M4與M5已經(jīng)失穩(wěn)，故未在圖中顯示，當(dāng)kΔx取較小值時(shí)，M6頻散最小，但kΔx取較大值時(shí)，M6頻散將超過M3.由圖對比可知，M3頻散一直維持在較低范圍.當(dāng)r＝0.8時(shí)，按照Basabe和Sen［40］頻散限定標(biāo)準(zhǔn)，即頻散小于0.01，1D情況下，M3一個(gè)波長的最小采樣點(diǎn)數(shù)為2.38；同理2D情形下，M3一個(gè)波場內(nèi)的最小采樣點(diǎn)數(shù)為3.36.

圖1 6種格式r＝0.5和r＝0.8時(shí)的頻散關(guān)系對比Fig.1 Comparison of numerical dispersion of six symplectic schemes when r＝0.5（a）and r＝0.8（b）

4 數(shù)值試驗(yàn)

4.1 半無限均勻介質(zhì)模型

為了測試這組辛RKN格式對聲波的模擬精度，設(shè)計(jì)了2D均勻介質(zhì)模型，模型大小為2550m×2550m，震源位于模型中心，接收點(diǎn)在震源左側(cè)400m處.震源由主頻為30Hz的Ricker子波激發(fā)，空間網(wǎng)格間隔為10m，時(shí)間步長為1ms.聲波波速為3000m／s.定義接收點(diǎn)的總體誤差其中，ui為t＝ih時(shí)刻的數(shù)值解，ur為解析解.模擬結(jié)果如圖2，計(jì)算效率與總體誤差如表3所示.

由圖2可知，6種方法模擬結(jié)果與解析解高度一致，M1、M2和 M6與解析解最靠近，M3、M4和M5偏離解析解較遠(yuǎn)，由于壓制高波數(shù)處的數(shù)值頻散，M3對低波數(shù)模擬的精度影響較大，表現(xiàn)在圖2中為相位較為超前.由表3可以看出，M6誤差最小，M1與M2誤差其次，其中M1與M2精度足夠接近；M3—M5誤差較大；就計(jì)算效率而言M6的計(jì)算時(shí)間為M1—M5的1.5倍，在計(jì)算效率相同的情況下，M1、M2比M4、M5有較大的精度提升.

4.2 分層介質(zhì)模型

為了檢驗(yàn)本文方法的數(shù)值頻散壓制能力，設(shè)計(jì)如圖3的分層介質(zhì)模型，模型大小為3825m×3825m，分界面位于1920m深度處，上層介質(zhì)的波速為2000m／s，下層介質(zhì)波速為3000m／s，震源位于（1920m，－1710m）處，由主頻為60Hz的Ricker子波激發(fā)產(chǎn)生，時(shí)間間隔為1ms，空間間隔為15m.模擬結(jié)果如圖4所示.

圖2 6種格式計(jì)算得到的聲波波形與解析解對比Fig.2 Comparison of waveforms generated by six sympletic methods and analytic solution

圖3 分層均勻介質(zhì)模型以及模型參數(shù)Fig.3 Two layer homogeneous media and model parameters

圖4（a—d）分別為M4—M6和M3四種方法計(jì)算得到的0.5s時(shí)刻的波場快照.從圖4a可以看到除了直達(dá)波、反射波、透射波以及首波外，還在直達(dá)波波前、透射波波前附近出現(xiàn)了強(qiáng)烈的數(shù)值頻散，圖

圖4 4種方法模擬分層聲波介質(zhì)中地震波傳播0.5s時(shí)的波場快照（a—c）M4—M6方法；（d）M3方法.Fig.4 The snapshots of scalar wave propagation in two layered medium at time t＝0.5s

表3 6種辛格式計(jì)算效率和總體誤差對比Table 3 Computational efficiency and total error comparison of six symplectic methods

4b在透射波波前也出現(xiàn)了強(qiáng)烈的數(shù)值頻散；圖4（c，d）中只在透射波前出現(xiàn)了輕微的頻散，M3和M6計(jì)算結(jié)果較為一致.M3—M5計(jì)算效率基本相同，M6計(jì)算耗時(shí)為M3—M5的1.5倍，在高頻地震波模擬時(shí)使用M3，在保證計(jì)算效率的同時(shí)，能較好地壓制數(shù)值頻散.

4.3 非均勻介質(zhì)模型——SEG／EAGE鹽丘模型

為了測試本文方法在非均勻介質(zhì)中波場計(jì)算的有效性和穩(wěn)定性，選擇SEG／EAGE模型做測試，模型速度結(jié)構(gòu)如圖5所示.模型中波速變化范圍為1524～4480m／s，除了若干個(gè)起伏分層界面外，還存在鹽丘高阻體，使得介質(zhì)橫向速度變化極為強(qiáng)烈.選擇震源位于模型中心，其為主頻是40Hz的Ricker子波，空間網(wǎng)格間距為20m，時(shí)間間隔為1ms時(shí)M2、M5和M6三種方法的模擬結(jié)果如圖6所示，時(shí)間間隔為2ms時(shí)的模擬結(jié)果如圖7所示.

當(dāng)時(shí)間間隔為1ms時(shí)，M2、M5和M6可以得到幾乎相同的波場快照，從圖6中可以看到，由于介質(zhì)的非均勻性，在速度突變的界面上產(chǎn)生了強(qiáng)烈的反射和多次反射，以及在速度突變的角點(diǎn)處產(chǎn)生了明顯的繞射和散射，三種方法均是穩(wěn)定的，并得到了可靠的結(jié)果.當(dāng)時(shí)間間隔增加一倍，即為2ms時(shí)，M5已經(jīng)失穩(wěn)，而M2和M6依然穩(wěn)定，從圖7中可以看到兩種方法得到的快照和圖6中的快照無論是整體還是細(xì)節(jié)上都極其一致，說明本文M2使用大時(shí)間步長的模擬結(jié)果依然是可靠的.

圖5 SEG／EAGE鹽丘模型Fig.5 The velocity profile of SEG／EAGE salt model

圖6 三種方法模擬得到的非均勻介質(zhì)中0.5s時(shí)的波場快照（a）M2；（b）M5；（c）M6，時(shí)間間隔為1ms.Fig.6 The snapshots of wavefields in heterogeneous media when t＝0.5s（a—c）are calculated by M2，M5，and M6，respectively.Time interval is 1ms.

圖7 兩種方法模擬得到的非均勻介質(zhì)中0.5s時(shí)的波場快照（a）M2；（b）M6，時(shí)間間隔為2ms.Fig.7 The snapshots of wavefields in heterogeneous media when t＝0.5s（a—b）are calculated by M2，and M6，respectively.Time interval is 2ms.

5 討 論

本文對地震聲波方程空間上使用偽譜法離散，

時(shí)間離散上選用高效的二階辛RKN格式，通過以誤差最小、穩(wěn)定域最大和頻散最小為依據(jù)構(gòu)造了三種辛RKN格式.在理論分析中，和常見格式對比，

理論上論證了本文方法在精度、穩(wěn)定性和數(shù)值頻散等方面的優(yōu)勢；在數(shù)值實(shí)驗(yàn)中，通過與解析解、分層介質(zhì)和非均勻介質(zhì)中不同方法波場模擬結(jié)果對比，

數(shù)值結(jié)果進(jìn)一步佐證了本文方法在保證計(jì)算效率的同時(shí)，在精度提高、穩(wěn)定域增加和數(shù)值頻散壓制等方面均具有明顯改進(jìn).可以根據(jù)不同實(shí)際需要選擇不同方法，如在高精度地震波模擬時(shí)選擇M1，在高頻地震波模擬時(shí)選擇M3，在強(qiáng)烈非均勻介質(zhì)中地震波模擬時(shí)選擇M2.本文方法為高效地震波模擬、成像提供了一種可靠的選擇.

（

）

［1］ Whitmore N D.Iterative depth migration by backward time propagation.53rd Annual International meeting，SEG，Expanded Abstracts，1983：827－830.

［2］ 劉紅偉，李博，劉洪等.地震疊前逆時(shí)偏移高階有限差分算法及GPU實(shí)現(xiàn).地球物理學(xué)報(bào)，2010，53（7）：1725－1733.Liu H W，Li B，Liu H，et al.The algorithm of high order finite difference pre－stack reverse time migration and GPU implementation.ChineseJ.Geophys.（in Chinese），2010，53（7）：1725－1733.

［3］ Song Z M，Williamson P R，Pratt R G.Frequency－domain acoustic－wave modelling and inversion of cross－h(huán)ole data，PartⅡ：Inversion method，synthetic experiments and realdata results.Geophysics，1995，60（3）：796－809.

［4］ 張美根，王妙月，李小凡等.時(shí)間域全波場各向異性彈性參數(shù)反演.地球物理學(xué)報(bào)，2003，46（1）：94－100.Zhang M G，Wang M Y，Li X F，et al.Full wavefield inversion of anisotropic elastic parameters in the time domain.ChineseJ.Geophys.（in Chinese），2004，46（1）：94－100.

［5］ Komatitsch D，Liu Q，Tromp J.Simulations of ground motion in the Los Angeles basin based upon the spectralelement method.Bull.Seismol.Am.Soc.，2004，94（1）：187－206.

［6］ 張懷，周元澤，吳忠良等.福州盆地強(qiáng)地面運(yùn)動特征的有限元數(shù)值模擬.地球物理學(xué)報(bào)，2009，52（5）：1270－1279.Zhang H，Zhou Y Z，Wu Z L，et al.Finite element analysis of seismic wave propagation characteristics in Fuzhou basin.ChineseJ.Geophys.（in Chinese），2009，52（5）：1270－1279.

［7］ Komatitsch D，Tromp J.Spectral－element simulations of global seismic wave propagation—I.Validation.Geophys.J.Int.，2002，149（2）：390－412.

［8］ Komatitsch D，Tromp J.Spectral－element simulations of global seismic wave propagation—Ⅱ.Three－dimensional models，oceans，rotation and self－gravitation.Geophys.J.Int.，2002，150（1）：303－318.

［9］ Yan Z Z，Zhang H，Yang C C，et al.Spectral element analysis on the characteristics of seismic wave propagation triggered by WenchuanMs8.0earthquake.ScienceinChina SeriesD：EarthSciences，2009，52（6）：764－773.

［10］ Carcione J M，Herman G C，ten Kroode A P E.Seismic modeling.Geophysics，2002，76（4）：1304－1325.

［11］ Yang D H，Tong P，Deng X Y.A central difference method with low numerical dispersion for solving the scalar wave equation.GeophysicalProspecting，2012，60（5）：885－905.

［12］ Virieux J.SH－wave propagation in heterogeneous media：Velocity－stress finite－difference method.Geophysics，1984，49（11）：1933－1942.

［13］ Virieux J.P－SV wave propagation in heterogeneous media：Velocity－stress finite－difference method.Geophysics，1986，51（4）：889－901.

［14］ Yang D H，Liu E，Zhang Z J，et al.Finite－difference modelling in two－dimensional anisotropic media using a fluxcorrected transport technique.Geophys.J.Int.，2002，148（2）：320－328.

［15］ Kolsoff D D，Baysal E.Forward modeling by a Fourier method.Geophysics，47（10）：1402－1412.

［16］ 龍桂華，李小凡，江東輝.基于交錯(cuò)網(wǎng)格Fourier偽譜微分矩陣算子的地震波場模擬GPU加速方案.地球物理學(xué)報(bào)，2010，53（12）：2964－2971.Long G H，Li X F，Jiang D H.Accelerating seismic modeling with staggered－grid Fourier Pseudo－spectral differentiation matrix operator method on graphics processing unit.ChineseJ.Geophys.（in Chinese），2010，53（12）：2964－2971.

［17］ Komatitsch D，Barnes C，Tromp J.Wave propagation near a fluid－solid interface：A spectral－element approach.Geophysics，2000，65（2）：623－631.

［18］ Komatitsch D，Barnes C，Tromp J.Simulation of anisotropic wave propagation based upon a spectral element method.Geophysics，2000，65（4）：1251－1260.

［19］ Dablain M A.The application of high－order differencing to the scalar wave equation.Geophysics，1986，51（1）：54－66.

［20］ Lax P D，Wendroff B.Difference schemes for hyperbolic equations with high order of accuracy.Communicationson PureandAppliedMathematics，1964，17（3）：381－398.

［21］ Tong P，Yang D H，Wang M X.A high－order stereomodeling method for solving wave equations.Bull.Seism.Am.Soc.，2013，103（2A）：811－833.

［22］ Yang D H，Teng J W，Zhang J F，et al.A nearly analytic discrete method for acoustic and elastic wave equations in anisotropic media.Bull.Seismol.Soc.Am.，2003，93（2）：882－890.

［23］ Chen J B. High－order time discretizations in seismic modeling.Geophysics，2007，72（5）：115－122.

［24］ Chen J B.Modeling the scalar wave equation with Nystr?m methods.Geophysics，2006，71（5）：151－158.

［25］ Chen J B.Lax－Wendroff and Nystr?m methods for seismic modelling.GeophysicalProspecting，2009，57（6）：931－941.

［26］ Okunbor D，Skeel R D.Explicit canonical methods for Hamiltonian systems.Math.Comp.，1992，59（200）：439－455.

［27］ Li X F，Li Y Q，Zhang M G，et al.Scalar seismic－wave equation modeling by a multisymplectic discrete singular convolution differentiator method.Bull.Seismol.Am.Soc.，2011，101（4）：1710－1718.

［28］ Li X F，Wang W S，Lu M W，et al.Structure－preserving modelling of elastic waves：a symplectic discrete singular convolution differentiator method.Geophys.J.Int.，188（3）：1382－1392.

［29］ Luo M Q，Liu H，Li Y M.LU decomposition with spectral factorization in seismic imaging.ChineseJ.Geophys.，2003，46（3）：602－612.

［30］ Luo M Q，Zhu G T，Liu H，et al.A hybrid matrix inversion method for 3－D implicit prestack depth migration.ChineseJ.Geophys.，2003，46（5）：978－987.

［31］ Ma X，Yang D H，Liu F Q.A nearly analytic symplectically partitioned Runge－Kutta method for 2－D seismic wave equations.Geophys.J.Int.，2011，187（1）：480－496.

［32］ 孫耿.波動方程的一類顯式辛格式.計(jì)算數(shù)學(xué)，1997，（1）：1－10.Sun G.A class of explicitly symplectic schemes for wave equations.Comput.Math.（in Chinese），1997，（1）：1－10.

［33］ Fornberg B.High－order finite differences and the pseudospectral method on staggered grids.SIAMJ.Numer.Anal.，1990，27（4）：904－918.

［34］ 汪文帥，李小凡，魯明文等.基于多辛結(jié)構(gòu)譜元法的保結(jié)構(gòu)地震波場模擬.地球物理學(xué)報(bào)，2012，55（10）：3427－3439.Wang W S，Li X F，Lu M W，et al.Structure－preserving modeling for seismic wavefields based upon a multisymplectic spectral element method.ChineseJ.Geophys.（in Chinese），2012，55（10）：3427－3439.

［35］ Yang D H，Wang N，Chen S，et al.An explicit method based on the implicit Runge－Kutta algorithm for solving wave equations.Bull.Seismol.Soc.Am.，2009，99（6）：3340－3354.

［36］ Yang D H，Wang Lei，Deng X Y.An explicit split－step algorithm of the implicit Adams method for solving 2－D acoustic and elastic wave equations.Geophys.J.Int.，2010，180（1）：291－310.

［37］ 馮康，秦孟兆.哈密爾頓系統(tǒng)的辛幾何算法.杭州：浙江科學(xué)出版社，2003.Feng K，Qin M Z.Symplectic Geometric Algorithms for Hamiltionian Systems （in Chinese）.Hangzhou：Zhejiang Science ＆ Technology Press，2003.

［38］ Hairer E.Geometric Numerical Intergration I （2nd ed）.Berlin and New York：Springer－Verlag.

［39］ MaLachlan R I， Atela P. The accuracy of symplectic integrators.Nonlinearity，1992，5（2）：541－562.

［40］ Basabe J D D，Sen K M.Grid dispersion and stability criteria of some common finite－element methods for acoustic and elastic wave equations.Geophysics，2007，72（6）：T81－T95.