李一瓊，李小凡，張美根

1 中國科學(xué)院地質(zhì)與地球物理研究所，北京 100029

2 中國地震局地球物理研究所，北京 100081

基于辛格式離散奇異褶積微分算子的彈性波場模擬

李一瓊1，2，李小凡1，張美根1

1 中國科學(xué)院地質(zhì)與地球物理研究所，北京 100029

2 中國地震局地球物理研究所，北京 100081

本文發(fā)展了基于辛格式離散奇異褶積微分算子（SDSCD）的保結(jié)構(gòu)方法模擬彈性波場，求解彈性波動方程時，引入辛差分格式進行時間離散，采用離散奇異褶積微分算子進行空間離散.相比于傳統(tǒng)的偽譜方法，該方法提高了計算精度和穩(wěn)定性.數(shù)值結(jié)果表明SDSCD方法可以有效地抑制數(shù)值頻散，為解決大尺度、長時程地震波場模擬問題提供了合適的數(shù)值方法.

辛算法，離散奇異褶積微分算子，彈性波場模擬

1 引 言

研究地震波在地球介質(zhì)中的傳播規(guī)律，是地球物理學(xué)中認識地球內(nèi)部結(jié)構(gòu)的一項重要手段.在地震波場數(shù)值模擬方法中，波動方程數(shù)值模擬方法一直占有重要的地位，波動方程模擬包含了地震波豐富的動力學(xué)信息，為研究地震波的傳播機理提供了更多的依據(jù).描述地球內(nèi)部的地震波場的波動方程多采用聲波方程和彈性波方程.聲波方程應(yīng)用于地震波研究，雖已獲得諸多成果，但聲波模型多適用于勘探地震學(xué).因地球介質(zhì)一般被視為非均勻固體介質(zhì)，在理論地震學(xué)研究中通常須考慮其彈性性質(zhì)及地震波的矢量性質(zhì)，簡單、機械地套用聲波近似，其精度及物理合法性均將受到質(zhì)疑.彈性波近似較為繁復(fù)，更貼近實際、更廣泛地應(yīng)用于地震波研究.本質(zhì)上，聲波近似可視為彈性波近似的一個特例和簡化.因此，多采用彈性波動方程描述地震波.

地震波波動方程的定解問題包括微分算子、震源、初始條件和邊界條件等.常見的正演數(shù)值模擬方法依據(jù)處理微分項的數(shù)值方式而不同，如有限差分法［1］、偽譜法［2］和有限元法［3］等都已被運用多年，較為成熟.有限差分法是基于Taylor展開，用差分代替微分，將微分方程轉(zhuǎn)化為代數(shù)方程，在計算上易于實現(xiàn)，但有些差分格式是根據(jù)微分算子憑經(jīng)驗湊算得到的，缺乏數(shù)學(xué)理論保證［4］.偽譜法采用了快速傅里葉變換計算空間微分，使其具有較快的計算速度，計算精度要高于經(jīng)典有限差分法.但計算空間微分時動用了區(qū)域內(nèi)的所有變量，有悖于微分的局部屬性，對于復(fù)雜介質(zhì)的地震波模擬結(jié)果并不理想［5－6］.以上兩種數(shù)值方法都不可避免地存在差分近似和網(wǎng)格離散化近似，以及無窮級數(shù)的有限截斷近似，將破壞波動方程及解的結(jié)構(gòu)，造成數(shù)值解相對于精確解的偏離.這一數(shù)值方法的固有問題，在一定范圍內(nèi)，顯現(xiàn)得不是很突出，從而被認為是有效的數(shù)值方法，其計算結(jié)果被接納.與前兩種算法相比，有限元法有較強的數(shù)學(xué)理論保障，其理論基礎(chǔ)是變分原理，離散化構(gòu)建在Sobolev函數(shù)空間框架下［7］，通過將單位元上的分片插值函數(shù)和最小位能理論相結(jié)合求解整個區(qū)域的波場信息，但其計算過程復(fù)雜，計算量大，應(yīng)用到復(fù)雜介質(zhì)中還有一定的困難.

在研究大尺度的地震波場模擬時，我們關(guān)注的是較長時程的地震波計算問題，其波動方程在時間域處理時，具有很強的時空繼承性，其積累誤差較為明顯，需要選擇具有高精度的計算方法.辛算法這一保結(jié)構(gòu)算法正是在諸如這樣的需求下提出的，算法構(gòu)建在哈密爾頓體系下，其構(gòu)造之初，就要求在方程進行離散化和有限化近似過程中，盡可能保持運動方程的結(jié)構(gòu)及其相應(yīng)的守恒量不變，或者，把數(shù)值解對于精確解的偏離控制在設(shè)定的精度以內(nèi)，使得辛算法能從根本上提高解的精度和可信度.計算實驗顯示，辛算法優(yōu)異的穩(wěn)定性和長時間跟蹤能力有著重要的應(yīng)用前景.

辛算法在眾多科學(xué)計算領(lǐng)域得到了應(yīng)用和發(fā)展.地震波場在高頻近似條件下，其射線解刻畫了地震波場的射線走向和波陣面形態(tài)，這在地震波場傳播的漸近解研究、模型研究及波場成像的運用中具有廣泛的實用價值.然而，算法本身的一個重大缺陷是在焦散區(qū)失效［8］，1965年前蘇聯(lián)學(xué)者Maslov將射線問題轉(zhuǎn)化為辛空間中的Lagrange子流形問題，成功地解決了焦散問題，從而大大拓廣了幾何光學(xué)應(yīng)用范圍.陳景波［9］和高亮［7］等用辛幾何算法解決地震射線追蹤問題得到了較好的結(jié)果，和四階Runge－Kutta方法相比，精度相當(dāng)，但辛算法的計算速度要快.在基于波動方程的地震波數(shù)值模擬方面，Chen［10］將辛格式和偽譜法相結(jié)合求解了標量波波動方程，和非辛格式的算法相比，辛格式顯示了較小的誤差增長.運用辛算法能很好地解決時間域波動方程的時間繼承性，因此，我們將為波動方程中時間微分項引入辛差分格式的保結(jié)構(gòu)算子.

波動方程中空間微分項的近似計算同樣至關(guān)重要.20世紀70年代早期，就有學(xué)者提出褶積微分算子的概念，但將其應(yīng)用到地震波場模擬中比較晚.Holberg［11］運用Fourier變換，對褶積微分算子進行了設(shè)計，并在三維彈性波數(shù)值模擬中取得了較好的結(jié)果.Zhou等［12］、張中杰等［13］和戴志陽等［14］也對這一方法進行了深入研究和發(fā)展.龍桂華等［15］在前人的工作基礎(chǔ)上，根據(jù)廣義函數(shù)理論，引入奇異核褶積微分算子，經(jīng)過加窗處理的微分算子與波場變量進行褶積，來求取波場變量關(guān)于空間微分.褶積微分算子和偽譜法相比，其優(yōu)勢在于褶積微分算子可以優(yōu)化為短算子，更能突出空間微分本身的局部屬性；而偽譜法在計算微分時須利用全局信息.褶積微分算子和有限差分算子相比，精度更高，能夠彌補有限差分算子描述復(fù)雜介質(zhì)的不足.因此褶積微分算子應(yīng)用于波動方程的空間微分近似是非常有效的.

本文對彈性波動方程引入顯辛差分格式進行時間離散，采用離散奇異褶積微分算子進行空間離散，建立一個新的彈性波場模擬方法，辛格式離散奇異褶積微分算子方法（SDSCD）.該方法盡管會增加一些計算量，但將提高計算精度和穩(wěn)定性.與傳統(tǒng)的偽譜方法相比，SDSCD方法是全局保結(jié)構(gòu)的，具有高精度和較強的長時間跟蹤能力，將為解決大尺度、長時程地震波場模擬問題提供有效的數(shù)值方法.

2 彈性波動方程的辛差分格式

從彈性動力學(xué)的本構(gòu)方程出發(fā)，可以推導(dǎo)出各向同性介質(zhì)中一階速度應(yīng)力方程：其中，vx和vz分別為沿x和z方向的速度場，σxx、σzz和σxz為應(yīng)力分量，λ和μ為表征介質(zhì)屬性的彈性參數(shù)，ρ為介質(zhì)密度.

令v＝（vx，vz），σ＝（σxx，σzz，σxz，σzx），則（v，σ）＝（vx，vz，σxx，σzz，σxz，σzx），得到方程（1）的矩陣形式：

注意到，矩陣P，Q中含有關(guān)于空間變量的微分算子，處理這些空間微分項有賴于第3節(jié)的奇異核褶積微分算子，暫不詳述.

由我們的先前工作［16－17］，可得到（2）式的k級n階顯式辛格式

其中，波場變量的上標（t）表示時間序列中的當(dāng)前時間點，（t＋1）為下一時間點，Δt表示時間步長，矩陣I為與P，Q同階的單位矩陣.

當(dāng)k＝3時，（3）式寫為彈性波動方程（2）的單步三級三階顯辛格式解［18］

式（4）中的空間導(dǎo)數(shù)可采用第3節(jié)的奇異核褶積微分算子計算.

3 離散奇異褶積微分算子

根據(jù)廣義函數(shù)理論，當(dāng)已知某函數(shù)的離散序列｛u（xk）｝，該函數(shù)的微分形式可由奇異核函數(shù)的微分來近似［19］

其中，δ（x）是奇異核Delta函數(shù)，（q）表示關(guān)于x的q階廣義導(dǎo)數(shù)，W為算子半帶寬.W和高斯包絡(luò)σ是可變的，適當(dāng)選取可提高微分的精度.

奇異核Delta函數(shù)常被用于偏微分方程的數(shù)值求解方面，在數(shù)值計算中，需要考慮用近似函數(shù)來逼近物理上很難實現(xiàn)的理論Delta函數(shù).這里選用了歸一化的Shannon奇異核函數(shù)

其中，Δ為空間步長，σ為高斯函數(shù)包絡(luò)線寬度.對于任意給定的σ≠0，當(dāng)Δ→0時，（6）式趨近于Delta函數(shù).

當(dāng)q＝1時，一階微分算子表示為

這個褶積微分算子系數(shù)關(guān)于k＝0是中心對稱的.

在實際計算過程中，對上述算子進行了加窗處理，從而有效地抑制由于算子截斷而引起的Gibbs效應(yīng).選用Hanning窗函數(shù)

得到處理后的微分算子：

在這里，選取Hanning窗函數(shù)的參數(shù)α＝0.54和β＝8.

4 算法的穩(wěn)定性條件

在時間域地震波場數(shù)值模擬中，計算的時間步長與空間步長的離散尺度必須滿足一定的條件，才能保證迭代計算的可行.否則，可能產(chǎn)生隨指數(shù)量級增長的計算誤差，造成嚴重的數(shù)值頻散，最終導(dǎo)致數(shù)值溢出使計算無法進行.對一種數(shù)值解法，需要知道計算穩(wěn)定的離散參數(shù)區(qū)間，即分析解法的穩(wěn)定性.

在第2節(jié)得到波場變量關(guān)于時間迭代的辛差分顯格式（3）式，現(xiàn)僅考慮三級三階的格式.令矩陣Mi＝（I＋ciΔtP）（I＋diΔt Q），i＝1，2，3，矩陣M＝M3M2M1，則（3）式寫為

形如（10）式的差分格式算法穩(wěn)定的充要條件是矩陣M的譜半徑ρ（M）≤1.由于矩陣P，Q中含有關(guān)于空間變量的微分，不能直接考察矩陣M，需要采取Fourier分析法［20－21］.

利用第3節(jié)的褶積微分算子，可得函數(shù)u（x）微分的離散近似表達式

將（11）式兩邊進行Fourier變換，并由算子系數(shù)的中心對稱，得到：

對（10）式進行關(guān)于空間變量的Fourier變換，P，Q中關(guān)于空間變量的微分項由（11）式的算子代替，變換后的矩陣M仍記為M.

考慮如下定理，設(shè)A∈Cn×n，則對Cn×n上的任一矩陣范數(shù)‖·‖，有ρ（A）≤‖A‖.

由此算法穩(wěn)定的充要條件等價為矩陣M的任一范數(shù)‖M‖≤1.

在實際推導(dǎo)過程中發(fā)現(xiàn)，矩陣M不具有特殊性質(zhì)，不能給出‖M‖一個比較簡單的表達式，故將數(shù)值實驗選取的參數(shù)代入矩陣M，采取試算的方式，來保證滿足穩(wěn)定性條件.在數(shù)值實驗中，對于均勻介質(zhì)模型，采用SDSCD方法，最大Courant數(shù)為0.684，而在同一模型中，常規(guī)Fourier偽譜算法的最大Courant數(shù)為0.64.

5 數(shù)值實驗

為了驗證方法的有效性，將本文提出的SDSCD算法結(jié)果與非辛算法的結(jié)果進行比較.這里所說的非辛算法是指空間離散采用Fourier偽譜算子近似，時間離散用向前差分近似的常規(guī)Fourier偽譜算法.下文中除特別說明外，非辛算法均指此偽譜算法.

在驗證算例中，首先采用了一個具有粗糙傾斜界面的分層介質(zhì)模型.模型的物性參數(shù)見圖1，五角星★表示震源所在位置（1280，1080），三角形▲為檢波點所在位置（1180，980）；震源為脹縮源，由主頻為30Hz的Ricker子波激發(fā)產(chǎn)生；時間步長為1ms，模型沿橫向和縱向各剖分256個網(wǎng)格點，空間間隔均為10m.

圖2和圖3分別給出了SDSCD算法和偽譜算法在0.3s時刻的波場快照.這些快照都非常清晰地刻畫了該模型中的直達P波、反射P波、反射P－S波、透射P波、透射P－S波以及散射波.可以看到模型上下兩層的波速對比度較大，散射波信息十分豐富.圖4和圖5給出了檢波點處0.3s內(nèi)兩種算法分別在水平波場和垂直波場的地震記錄.波場快照和地震記錄都說明了SDSCD算法的有效性.

圖1 具有粗糙傾斜界面的分層均勻介質(zhì)模型Fig.1 Homogeneous layered model with rough inclined interface

為了檢驗SDSCD算法的保結(jié)構(gòu)特性，我們對一個均勻介質(zhì)模型的波場情況進行了長時程的跟蹤，模型邊緣沒有加任何邊界條件.此做法相當(dāng)于增大了模型的尺度，通過進行較多次時間迭代步，考察算法是否保持穩(wěn)定.Chen［22］采用同樣的做法對標量波模型的辛Nystr?m算法進行了驗證.

模型的物性參數(shù)為：縱波波速4000m／s，橫波波速2309.3m／s；密度2400kg／m3，模型沿橫向和縱向各剖分出256個網(wǎng)格點，空間間隔均為20m；震源位于模型中央（2560m，2560m）處，震源為集中力源，由主頻為20Hz的Ricker子波激發(fā)產(chǎn)生；時間步長為2ms.圖6和圖7分別給出了SDSCD算法和偽譜算法在300次時間迭代步后，即0.6s時刻的水平波場快照，可以看到，偽譜算法的結(jié)果數(shù)值頻散較為明顯.圖8和圖9分別給出了SDSCD算法和偽譜算法在5000次時間迭代步后，即10s時刻的水平波場快照，SDSCD算法模擬的波前面十分清晰，而偽譜算法的結(jié)果，則很難辨清波前的形狀.由此可知，SDSCD方法具有較明顯的保結(jié)構(gòu)性，能夠高精度處理長時程波場模擬.

6 結(jié)論與討論

本文給出SDSCD算法用于彈性波場模擬，該方法在計算過程中，存儲較少，易于實現(xiàn)，雖然相比本文用于對比的經(jīng)典偽譜方法，增加了一些計算迭代步驟和時間，但在抑制數(shù)值頻散和提高精度方面都有明顯改善，并且具有較強的穩(wěn)定性.可以在一定范圍內(nèi)加大時間步長和空間間隔，作為研究大尺度、長時程的地震波場模擬的有效方法.

作為方法研究，本文還有一些工作尚待完成，如SDSCD算法的穩(wěn)定性條件尚未給出解析表達式，僅采取了試算，以保證數(shù)值實驗的可行.這種方式對于客觀評價本算法的優(yōu)劣略顯不便.這需要進一步的理論研究，使之更為完善.

（References）

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

［2］ Fornberg B.The pseudospectral method：Comparisons with finite differences for the elastic wave equation.Geophysics，1987，52（4）：483－501.

［3］ Smith W D.The application of finite element analysis to body wave propagation problems.Geophys.J.Int.，1975，42（2）：747－768.

［4］ 鐘萬勰.應(yīng)用力學(xué)的辛數(shù)學(xué)方法.北京：高等教育出版社，2006.Zhong W X.Symplectic Solution Methodology in Applied Mechanics（in Chinese）.Beijing：Higher Education Press，2006.

［5］ ?zdenvar T，McMechan G A.Causes and reduction of numerical artifacts in pseudo－spectral wavefield extrapolation.Geophys.J.Int.，1996，126（3）：819－828.

［6］ Li X F，Zhu T，Zhang M G，et al.Seismic scalar wave equation with variable coefficients modeling by a new convolutional differentiator.Computer Physics Communications，2010，181（11）：1850－1858.

［7］ 高亮，李幼銘，陳旭榮等.地震射線辛幾何算法初探.地球物理學(xué)報，2000，43（3）：402－410.Gao L，Li Y M，Chen X R，et al.An attempt to seismic ray tracing with symplectic algorithm.Chinese J.Geophys.（in Chinese），2000，43（3）：402－410.

［8］ 李世雄.波動方程的高頻近似與辛幾何.北京：科學(xué)出版社，2001.Li S X.The High－Frequency Approximation of Wave Equations and Symplectic Algorithm（in Chinese）.Beijing：Science Press，2001.

［9］ 陳景波，秦孟兆.辛幾何算法在射線追蹤中的應(yīng)用.數(shù)值計算與計算機應(yīng)用，2000，21（4）：255－265.Chen J B，Qin M Z.Ray tracing by symplectic algorithm.Numerical Computation and Computer Application（in Chinese），2000，21（4）：255－265.

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

［11］ Holberg O.Computational aspects of the choice of operator and sampling interval for numerical differentiation in largescale simulation of wave phenomena.Geophysical Prospecting，1987，35（6）：629－655.

［12］ Zhou B，Greenhalgh S A.Seismic scalar wave equation modeling by a convolutional differentiator.Bull.Seism.Soc.Am.，1992，82（1）：289－303.

［13］ 張中杰，滕吉文，楊頂輝.聲波與彈性波場數(shù)值模擬中的褶積微分算子法.地震學(xué)報，1996，18（1）：63－69.Zhang Z J，Teng J W，Yang D H.The convolutional differentiator method for numerical modeling of acoustic and elastic wave－field.Acta Seismologica Sinica（in Chinese），1996，18（1）：63－69.

［14］ 戴志陽，孫建國，查顯杰.地震波場模擬中的褶積微分算子法.吉林大學(xué)學(xué)報（地球科學(xué)版），2005，35（4）：520－524.Dai Z Y，Sun J G，Zha X J.Seismic wave field modeling with convolutional differentiator algorithm.Journal of Jilin University（Earth Science Edition）（in Chinese），2005，35（4）：520－524.

［15］ 龍桂華，李小凡，張美根.基于Shannon奇異核理論的褶積微分算子在地震波場模擬中的應(yīng)用.地球物理學(xué)報，2009，52（4）：1014－1024.Long G H，Li X F，Zhang M G.The application of convolutional differentiator in seismic modeling based on Shannon singular kernel theory.Chinese J.Geophys.（in Chinese），2009，52（4）：1014－1024.

［16］ 李一瓊，李小凡，朱童.基于辛格式奇異核褶積微分算子的地震標量波場模擬.地球物理學(xué)報，2011，54（7）：1827－1834.Li Y Q，Li X F，Zhu T.The seismic scalar wave field modeling by symplectic scheme and singular kernel convolutional differentiator.Chinese J.Geophys.（in Chinese），2011，54（7）：1827－1834.

［17］ 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.Seism.Soc.Am.，2011，101（4）：1710－1718.

［18］ Iwatsu R.Two new solutions to the third－order symplectic integration method.Phys.Lett.A.，2009，373（34）：3056－3060.

［19］ Wei G W.Discrete singular convolution for the solution of the Fokker－Planck equations.J.Chem.Phys.，1999，110（18）：8930－8942.

［20］ 黃志祥，沙威，吳先良等.高階辛算法的穩(wěn)定性與數(shù)值色散性分析.計算物理，2010，27（1）：82－88.Huang Z X，Sha W，Wu X L，et al.Stability and numerical dispersion of high order symplectic schemes.Chinese Journal of Computational Physics（in Chinese），2010，27（1）：82－88.

［21］ 董良國，馬在田，曹景忠.一階彈性波方程交錯網(wǎng)格高階差分解法穩(wěn)定性研究.地球物理學(xué)報，2000，43（6）：856－864.Dong L G，Ma Z T，Cao J Z.A study on stability of the staggered－grid high－order difference method of first－order elastic wave equation.Chinese J.Geophys.（in Chinese），2000，43（6）：856－864.

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

The elastic wave fields modeling by symplectic discrete singular convolution differentiator method

LI Yi－Qiong1，2，LI Xiao－Fan1，ZHANG Mei－Gen1
1 Institute of Geology and Geophysics，Chinese Academy of Sciences，Beijing100029，China
2 Institute of Geophysics，China Earthquake Administration，Beijing100081，China

In this paper，we introduce a structure－preserving method based on symplectic discrete singular convolution differentiator（SDSCD）for simulating elastic wave fields.In the method presented for solving elastic wave equations，physical space is discretized by singular convolution differentiator，whereas a symplectic difference scheme is used for the time discretization.The computational accuracy and stability of this method have been greatly improved compared with traditional pseudo－spectral method.Numerical results suggest the SDSCD algorithm can suppress effectively numerical dispersion，and it is suitable for modeling the large－scale and long－term seismic wave propagation.

Symplectic algorithm，Discrete singular convolution differentiator，Elastic wave fields modeling

10.6038／j.issn.0001－5733.2012.05.029

P631

2011－11－03，2012－03－28收修定稿

國家自然科學(xué)基金重點項目（40437018）、國家自然科學(xué)基金項目（40874024，41174047）、國家重點基礎(chǔ)研究發(fā)展計劃“973”計劃（2007CB209603）和中國地震局地球物理研究所基本科研業(yè)務(wù)費專項（DQJB11B10）聯(lián)合資助.

李一瓊，女，滿族，1979年生，2011年于中國科學(xué)院地質(zhì)與地球物理研究所獲博士學(xué)位，現(xiàn)為中國地震局地球物理研究所在站博士后，主要從事地震波傳播理論和數(shù)值模擬的研究工作.E－mail：liyq＠m(xù)ail.igcas.ac.cn，liyq＠cea－igp.ac.cn

李一瓊，李小凡，張美根.基于辛格式離散奇異褶積微分算子的彈性波場模擬.地球物理學(xué)報，2012，55（5）：1725－1731，

10.6038／j.issn.0001－5733.2012.05.029.

Li Y Q，Li X F，Zhang M G.The elastic wave fields modeling by symplectic discrete singular convolution differentiator method.Chinese J.Geophys.（in Chinese），2012，55（5）：1725－1731，doi：10.6038／j.issn.0001－5733.2012.05.029.

（本文編輯 何 燕）