分段光滑曲線邊界波動方程數(shù)值模擬研究

王忠仁，劉瑞，趙博雄

吉林大學(xué)儀器科學(xué)與電氣工程學(xué)院地球信息探測儀器教育部重點實驗室，長春 130026

1 引言

隨著人類對礦產(chǎn)資源需求的不斷增加，地表和淺層的金屬礦產(chǎn)資源越來越少，尋找深部隱伏金屬礦產(chǎn)已成為拓展新的找礦空間的必由之路.在勘探深度和精度方面都有明顯優(yōu)勢且在油氣勘探領(lǐng)域取得廣泛成功的地震勘探方法有希望在尋找深部隱伏金屬礦產(chǎn)方面發(fā)揮作用.由于金屬礦多產(chǎn)出在造山帶，因而通常表現(xiàn)為起伏地表的復(fù)雜特征.不僅如此，金屬礦區(qū)往往經(jīng)歷強(qiáng)烈的構(gòu)造變動和巖漿活動，礦體形狀復(fù)雜，在空間尺度上要遠(yuǎn)遠(yuǎn)小于油氣藏.因而金屬礦地震勘探具有地表起伏地下地質(zhì)體復(fù)雜的顯著特征，其地震波傳播規(guī)律要比水平地表和層狀介質(zhì)情形復(fù)雜得多.采用地震數(shù)值模擬技術(shù)是揭示起伏地表復(fù)雜地質(zhì)體情形下地震波傳播規(guī)律的有效手段.

地震波方程求解的數(shù)值方法主要有有限元法（楊頂輝，2002）、邊界元法（Fu and Wu，2000）和有限差分法（王德利等，2005；鄭海山和張中杰，2005）.有限元法可有效處理復(fù)雜邊界，邊界點與內(nèi)點具有相同的離散處理手段，但采用有限元法求解與時間變量有關(guān)的發(fā)展方程，每個時間層都要求解一個大型線性代數(shù)方程組，計算速度很慢.邊界元法比有限元法計算量小，處理復(fù)雜外邊界比較有效，但是當(dāng)內(nèi)邊界較多時，該方法計算效率明顯降低.常規(guī)的矩形網(wǎng)格有限差分法具有計算速度快的顯著特點，但該方法處理復(fù)雜邊界比較困難.

在復(fù)雜邊界地震波方程數(shù)值模擬方面，國內(nèi)外許多學(xué)者開展了大量富有探索性的研究工作.薛東川等（2007）采用有限元法對起伏地表復(fù)雜介質(zhì)波動方程進(jìn)行正演模擬，Ge和Chen（2008）采用邊界元方法就復(fù)雜地下分界面情形進(jìn)行地震波場數(shù)值模擬，李緒宣等（2011）采用邊界元方法分析復(fù)雜海底的地震散射特征.采用矩形網(wǎng)格有限差分方法進(jìn)行復(fù)雜邊界地震波數(shù)值模擬仍然是國內(nèi)外研究的一個主流方向，Jih等（1988）對自由表面為折線邊界情形區(qū)分六種情況給出了邊界點處理的一種局部旋轉(zhuǎn)坐標(biāo)方法，Pérez-Ruiz等（2005）在自由表面之上引入了虛介質(zhì)的概念并區(qū)分13種情形對邊界點進(jìn)行處理，Appel?和Petersson（2009）采用曲線坐標(biāo)系到直角坐標(biāo)系映射方法處理復(fù)雜自由邊界，Lan和Zhang（2011）就三維情形利用坐標(biāo)映射方法處理非規(guī)則自由表面，裴正林（2004）采用零速度法和廣義虛像法相結(jié)合的辦法處理自由邊界，董良國等（2007）將自由邊界上的點劃分為22類并分別采用不同的差分格式處理邊界點.為了克服單一數(shù)值方法的局限性，又相繼出現(xiàn)了幾種組合方法.黃自萍等（2004）采用有限元與有限差分耦合的區(qū)域分裂法，克服了單純用有限差分法對區(qū)域的依賴性.馬德堂和朱光明（2004）采用有限元與偽譜法混合求解方法，并采用過渡區(qū)域解決了兩種算法的銜接問題.管西竹等（2011）采用邊界元和體積元結(jié)合的方法在頻域求解含有面積分和體積分的積分方程.此外，由于三角網(wǎng)差分如同有限元一樣可適應(yīng)復(fù)雜邊界，褚春雷和王修田（2005）采用三角網(wǎng)有限差分法求解復(fù)雜邊界地震波方程，孫小東等（2012）研究了基于三角網(wǎng)有限差分法的疊前逆時偏移技術(shù).在這些方法中，組合方法與三角網(wǎng)差分也都存在計算效率低的弱點.鑒于矩形網(wǎng)有限差分法在計算效率上具有突出的優(yōu)點，而一般的二維地球物理模型邊界都可由分段光滑曲線來逼近，本文研究二維分段光滑曲線邊界地震波方程矩形網(wǎng)格有限差分的有效算法和數(shù)值模擬技術(shù).

2 方法原理

采用矩形網(wǎng)格有限差分法求解地震波傳播的波動方程數(shù)學(xué)模型，有兩類邊界條件最難處理，其一是地表Γ上的自由邊界條件

其二是地下分界面Γs上的銜接條件

其中u為波場值（位移），v為介質(zhì)波速，n表示外法線方向，下腳標(biāo)U和D所標(biāo)注的分別是分界面上部和下部的相應(yīng)物理量.

矩形網(wǎng)格有限差分法處理復(fù)雜邊界問題的關(guān)鍵有兩點，其一是對復(fù)雜邊界進(jìn)行合理的網(wǎng)格化形成網(wǎng)格邊界點，其二是利用網(wǎng)格邊界點附近的波場值計算自由邊界條件（1）式和銜接條件（2）式中的方向?qū)?shù).

2.1 邊界網(wǎng)格化

如果邊界Γ是封閉曲線，可將Γ分解成多支依次連接的曲線，每支曲線都可由單值函數(shù)來表示.不妨假定邊界?？捎蓡沃岛瘮?shù)z＝φ（x），（xl≤x≤xr）來表示，其中xl和xr分別為邊界Γ左右端點的x坐標(biāo).按如下的就近遷移原則將實際邊界點遷移到最近的網(wǎng)格點形成網(wǎng)格邊界點：

（1）選定包含計算區(qū)域R在內(nèi)的一個矩形區(qū)域［0，X］×［0，Z］，其中X和Z分別為矩形區(qū)域的長度和高度，取定x和z方向的空間步長分別為Δx，Δz，設(shè)X／Δx＝M，Z／Δz＝N，則網(wǎng)格點坐標(biāo)為

（2） 設(shè)函數(shù)nint表示取整數(shù)運(yùn)算，計算Γ左右端點x坐標(biāo)對應(yīng)的水平網(wǎng)格點的序號為

（3）計算

其中jφi是邊界Γ上的φ（xi）對應(yīng)水平網(wǎng)格點序號i的網(wǎng)格化近似值.

如果沒有特殊需要，通常情況下都采用Δx＝Δz的正方形網(wǎng)格，本文計算過程中采用的就是正方形網(wǎng)格.

2.2 法線方向上波場值的插值計算

對于自由表面邊界條件來說，計算方向?qū)?shù)需要自由面以下介質(zhì)內(nèi)部沿邊界法線方向上的波場值；對于內(nèi)邊界銜接條件而言，計算方向?qū)?shù)需要內(nèi)部邊界兩側(cè)沿邊界法線方向上的波場值.一般情況下，法線與網(wǎng)格節(jié)點并不相交，因而法線上的波場值通常是未知的.本文給出了一種插值方法來計算法線方向上的波場值.

該插值方法在正方形網(wǎng)格條件下要求將角度區(qū)間 ［0°，180°］以45°為間隔分四種情況來討論，這里我們僅以［0°，45°］情況為例進(jìn)行論述.

首先假設(shè)邊界Γ是光滑曲線z＝φ（x），曲線在（xi，φ（xi））處的切線L與ox軸夾角為α，0≤α≤45°.將邊界點就近遷移到網(wǎng)格邊界點（i，j），過點（i，j）做平行于切線L的直線則直線可作為網(wǎng)格邊界點（i，j）處的切線.為敘述方便，不妨假定邊界Γ上的點恰好位于網(wǎng)格邊界點 （i，j）處，直線恰為點）處的切線，對應(yīng)的法線為與ox軸的夾角為α（如圖1）.

圖1 插值計算示意圖Fig.1 Sketch of interpolation calculation

設(shè)點 （i-1，j+1），（i，j+1）分別為D和C（如

同理：

在正方形網(wǎng)格條件下，由于Δz＝Δx，于是有

2.3 方向?qū)?shù)與邊界條件計算

在圖1中記（i，j）點為E，首先假設(shè)Γ是內(nèi)邊界.當(dāng)Γ作為上部介質(zhì)的邊界時，邊界Γ上點E處的外法線方向?qū)?shù)可按下式計算

當(dāng)Γ作為下部介質(zhì)的邊界時，邊界Γ上點E處的外法線方向?qū)?shù)的計算公式為

于是，銜接條件為

可寫成

如果假設(shè)圖1中Γ為自由面，且Γ下部為地球介質(zhì)，則自由面邊界條件為

差分后得

于是有

在上述公式（6）—（9）中，φ′（xi）僅僅是空間變量x的函數(shù)，與時間變量無關(guān)，在差分法的波場時間遞推計算之前，可一次性計算出φ′（xi）.因此，這種算法明顯具有復(fù)雜度低與計算量小的特點.

2.4 正則導(dǎo)數(shù)定義

在前面的討論中，公式（4）—（7）的計算要求φ′（xi）處處存在，也就是要求邊界Γ處處光滑，而實際地層中常出現(xiàn)斷層、尖滅等非光滑邊界情形，討論分段光滑邊界更具普遍意義.由于分段光滑邊界尖點處的導(dǎo)數(shù)不存在，要想使前述算法仍然可用，就必須給出尖點處的某種廣義導(dǎo)數(shù).

首先考慮分段光滑曲線邊界的特例——折線邊界，在折線邊界的尖點處我們引入一種正則導(dǎo)數(shù)，具體定義如下：

計算向量的模為

構(gòu)造單位向量為

圖2 尖點正則導(dǎo)數(shù)定義Fig.2 Definition of regular derivative for a cusp

構(gòu)造差向量為

對于一般的分段光滑曲線，可根據(jù)尖點兩側(cè)的常規(guī)左右導(dǎo)數(shù)的變化趨勢做出兩條在尖點相交的切線，設(shè)交點為p1，在兩條切線上分別取點p0和p2，與折線尖點正則導(dǎo)數(shù)類似可定義分段光滑曲線尖點處的正則導(dǎo)數(shù).

3 正演模擬計算與波場特征分析

3.1 正演模擬計算

構(gòu)造圖3所示的地球介質(zhì)模型I.地表Γ0為曲線自由邊界.區(qū)域內(nèi)部邊界Γ1為斜線邊界，內(nèi)部邊界Γ2為封閉折線邊界，Pk（k＝1，2，…，7）為折線上的尖點，區(qū)域內(nèi)部邊界都采用銜接條件.左邊界x＝0、右邊界x＝1000m和下邊界z＝1000m為人工截斷邊界，處理為吸收邊界.圖中v1是上層介質(zhì)的波速，v2是Γ1以下Γ2外部介質(zhì)的波速，v3是Γ2所圍區(qū)域內(nèi)部介質(zhì)的波速.

依據(jù)以上邊界條件建立如下的二維地震縱波數(shù)學(xué)模型為

圖3 地球介質(zhì)模型IFig.3 Earth medium model I

在由（10）—（17）式給出的數(shù)學(xué)模型中，（10）式是內(nèi)點的泛定方程，時間和空間的二階偏導(dǎo)數(shù)均采用二階中心差商進(jìn)行離散；（12）—（14）式分別是左邊界、右邊界和底邊界上的吸收邊界條件，由于吸收邊界是與坐標(biāo)軸平行的直線邊界，可用常規(guī)差分法直接離散；（15）式是地表Γ0上的自由邊界條件，（16）—（17）式分別是內(nèi)邊界Γ1和Γ2上的銜接條件，（15）—（17）式中的方向?qū)?shù)采用本文給出的插值法計算.

采用矩形網(wǎng)有限差分法對上述模型進(jìn)行地震正演數(shù)值模擬.震源子波采用主頻為62.5Hz的雷克子波，炮點坐標(biāo)為（500，5），地面接收排列的起點為x＝100，道間距為10m，設(shè)置80道，計算中取空間步長為Δz＝Δx＝1.25m，時間步長為Δt＝0.25ms，記錄長度為0.9s，模擬得到的波場系列快照和地震記錄見圖4和圖5.

3.2 波場特征分析

在圖4中的220ms波場快照中，T1是地震波經(jīng)過界面Γ1向下方介質(zhì)透射的波前，R1是對應(yīng)的反射波前，A1、B1位于內(nèi)邊界Γ1上.A、B兩點在曲線地表Γ0上，地表檢波器接收到的A、B兩點的波場是直達(dá)波，在圖5的地震記錄中對應(yīng)的是直達(dá)波同相軸D.在260ms波場快照中，弧與是介質(zhì)V2中的入射波前，T2是透射到介質(zhì)V3中的透射波前，R2是對應(yīng)的反射波前，A2位于內(nèi)邊界Γ2的線段上，B2位于內(nèi)邊界Γ2的線段上.

300ms波場快照中新出現(xiàn)的波前R3是Γ2上段產(chǎn)生的反射，R4是Γ2上P2P3段產(chǎn)生的反射，R4與R2形成的“蝴蝶結(jié)”是凹陷反射的結(jié)果.320ms波場快照中新出現(xiàn)的波前R5是Γ2上段產(chǎn)生的反射，從該快照中還可看出，當(dāng)反射波前R2向上透過Γ1時，波前的形狀、曲率以及能量分布的均勻性都發(fā)生了改變.350ms波場快照中新出現(xiàn)的波前R6是Γ2上段產(chǎn)生的反射，而反射波R1即將到達(dá)地表Γ0并將被地表檢波器所接收，在圖5的地震記錄中在350ms以后將開始出現(xiàn)反射波同相軸R1.

在400ms以后的波場快照中，除了反射波之外還將討論自由表面引起的多次波（包括二次入射波、二次反射波和二次透射波）.400ms波場快照中反射波R1經(jīng)過自由表面反射后產(chǎn)生了二次入射波DS1.在480ms和560ms的波場快照中，DS2是反射波R2經(jīng)過自由表面反射后生成的二次入射波.而且在此期間反射波前R2完整經(jīng)過了地表的檢波器排列，在圖5的地震記錄中出現(xiàn)了對應(yīng)的反射波同相軸R2；反射波前R3按道號由大到小的規(guī)律也被地表檢波器依次接收，在圖5中出現(xiàn)了對應(yīng)的反射波同相軸R3，并將繼續(xù)向小道號方向擴(kuò)展；反射波前R4只有波震面的一小部分被少量小道號地表檢波器所接收，在圖5中對應(yīng)的是同相軸R4，對應(yīng)時間是在550ms到600ms之間.在700ms的波場快照中，TS1是二次入射波DS1經(jīng)過介質(zhì)分界面產(chǎn)生的下行透射波前（這里可包含多次下行透射），M1是DS1入射到界面Γ1產(chǎn)生的二次反射波前；TS2是二次入射波DS2經(jīng)過介質(zhì)分界面產(chǎn)生的下行透射波前，M2是DS2入射到界面Γ1產(chǎn)生的二次反射波前.在600ms與700ms期間反射波前R5、R6也依次到達(dá)地表檢波器排列，在圖5的反射記錄中，反射波R5、R6表現(xiàn)為交叉“蝴蝶結(jié)”狀同相軸.而二次反射波前M1和M2將在720ms至850ms期間先后到達(dá)地表檢波器排列，在圖5的地震記錄上對應(yīng)兩個二次波同相軸M1和M2.

圖4 波場系列快照Fig.4 Snapshots of wave fields

3.3 邊界條件處理對比計算

我們研究的法向?qū)?shù)插值計算方法，自由表面外邊界處理只是內(nèi)邊界處理的一個簡單特例.為此我們選擇了另一種處理內(nèi)邊界的方法進(jìn)行比較，就是內(nèi)邊界不給定銜接條件（也就不涉及方向?qū)?shù)的計算），只用介質(zhì)參數(shù)（這里是波速）的不同分區(qū)來表征內(nèi)邊界.圖6是參數(shù)分區(qū)表征邊界情形模擬的地震記錄.

粗略地看，圖5與圖6給出的兩種計算結(jié)果似乎沒有太大差別.但仔細(xì)觀察就會發(fā)現(xiàn)，給定銜接條件計算法向?qū)?shù)的情形下，反射波的相位與入射波（或直達(dá)波）的相位是相反的，符合真實的物理規(guī)律；而在不給定銜接條件的情形下，反射波的相位與入射波（或直達(dá)波）的相位是相同的，這與真實的物理規(guī)律是相悖的.通過圖7給出的波場快照可以更明顯地看出兩種情形反射波相位的差別.

為了衡量插值算法的計算效率，我們采用節(jié)3.1給出的計算參數(shù)，在雙核處理器為AMD Athlon（t）II×2250的普通微機(jī)上進(jìn)行對比計算，銜接條件邊界情形的計算時間為16.5min，無銜接條件的參數(shù)分區(qū)表征邊界情形的計算時間為15.9min，二者的時間差就是處理內(nèi)部銜接條件邊界Γ1和Γ2所耗費(fèi)的機(jī)時，這表明對本文給定的模型而言，其內(nèi)部銜接條件邊界插值處理的計算量只占全部模擬過程所用計算量的3.6%，因此插值處理銜接邊界條件的計算效率是可以接受的.當(dāng)然對于不同模型不同總長度的內(nèi)邊界，其計算效率是不同的.

3.4 尖點導(dǎo)數(shù)對波場繞射特征的影響

圖5 正演模擬的地震記錄Fig.5 Seismic record by forward simulation

圖6 參數(shù)分區(qū)表征邊界情形模擬的地震記錄Fig.6 Seismic record at boundary condition characterized by parameter partition

圖7 220ms波場快照（a）銜接邊界情形；（b）參數(shù)分區(qū)表征邊界情形.Fig.7 Wave field snapshots at 220ms（a）The case at join condition boundary；（b）The case at boundary of parameter partition characterization.

在尖點正則導(dǎo)數(shù)的定義下，我們實現(xiàn)了對地球介質(zhì)模型I的數(shù)值模擬計算.為了研究尖點導(dǎo)數(shù)對波場特征的影響，就要給出另一種尖點導(dǎo)數(shù)，通過對比兩種不同的尖點導(dǎo)數(shù)對應(yīng)的波場及波場差，揭示不同的尖點導(dǎo)數(shù)處理技術(shù)對波場特征的影響規(guī)律.為此，我們又研究了另一種尖點導(dǎo)數(shù).對分段光滑曲線而言，尖點兩側(cè)以尖點為中心的去心鄰域內(nèi)曲線上的導(dǎo)數(shù)存在而且連續(xù)，在曲線上由尖點一側(cè)的導(dǎo)數(shù)向尖點處取單側(cè)極限，將尖點處的這個極限值作為尖點的一種廣義導(dǎo)數(shù)，我們稱之為單側(cè)連續(xù)導(dǎo)數(shù).

圖8 地球介質(zhì)模型ⅡFig.8 Earth medium modelⅡ

由于模型I過于復(fù)雜，掩蓋了尖點繞射的微弱特征，我們建立了模型II（圖8）用來討論尖點問題，除了介質(zhì)模型結(jié)構(gòu)不同之外，其余的計算參數(shù)與模型I相同.分別采用正則導(dǎo)數(shù)和單側(cè)連續(xù)導(dǎo)數(shù)兩種定義對模型II進(jìn)行數(shù)值計算，圖9a是采用正則導(dǎo)數(shù)計算得到的160、200ms和230ms的波場快照，圖9b是采用單側(cè)連續(xù)導(dǎo)數(shù)計算得到的對應(yīng)波場快照，圖9c是二者的波場差.圖9a中的C1、C2和C3分別是尖點P2、P3和P4產(chǎn)生的繞射波.繞射波前在均勻介質(zhì)中都是以尖點為圓心的圓周，與反射波相比繞射波能量相對較弱.由圖9a和b可見，兩種方法的模擬效果在宏觀上幾乎看不出差別.然而由圖9c的波場差可見，兩者存在著微小差別，這種差別產(chǎn)生的原因僅僅是兩種方法在尖點導(dǎo)數(shù)的處理技術(shù)上有所不同，而這種不同并沒有引起傳播規(guī)律的改變，改變的只是尖點對繞射波能量大小的貢獻(xiàn)略有不同.

3.5 穩(wěn)定性與數(shù)值頻散問題討論

泛定方程（10）所對應(yīng)的齊次方程的二階顯式差分格式為

其中Δx＝Δz＝h，p＝vΔt／h為網(wǎng)比，i，j，k是使得x＝iΔx，z＝j(luò)Δz，t＝kΔt的整數(shù).

圖9 模型Ⅱ波場快照（a）正則導(dǎo)數(shù)法；（b）單側(cè)連續(xù)導(dǎo)數(shù)法；（c）（a）與（b）波場差.Fig.9 Wave field snapshots for modelⅡ（a）Wavefields of regular derivative method；（b）Wavefields of unilateral continuous derivative method；（c）Wavefields difference between（a）and（b）.

差分格式（18）的穩(wěn)定性條件為（Alford etal.，1974）通過大量的數(shù)值計算表明，在采用差分格式（18）對內(nèi)點進(jìn)行離散的前提下，本文提出的處理邊界的法向?qū)?shù)插值算法在網(wǎng)比滿足（19）式的限制下是穩(wěn)定的，穩(wěn)定性條件（19）式?jīng)Q定了時間步長Δt與空間步長h的相互制約關(guān)系.

利用有限差分法數(shù)值求解波動方程僅靠穩(wěn)定性條件限制是不能有效解決數(shù)值頻散問題的.為了有效消除數(shù)值頻散誤差，通常要求在一個地震波長范圍內(nèi)涵蓋足夠多的空間網(wǎng)格點.本文正演模擬計算中取空間步長h＝1.25m，對給定的震源子波使得在一個最小地震波長內(nèi)涵蓋了23個空間網(wǎng)格點，有效地消除了數(shù)值頻散現(xiàn)象.通過減小空間步長的辦法雖然可以減弱數(shù)值頻散誤差，但同時也大大增加了計算量.近年來，（Ma etal.，2011；Yang etal.，2012a，2012b；Tong etal.，2013）研究了幾種具有弱數(shù)值頻散特性的數(shù)值方法，可以放寬空間步長的限制，并有效壓制數(shù)值頻散誤差.

4 結(jié)論

給出的邊界網(wǎng)格化方法可方便地處理橫向變化劇烈的內(nèi)外邊界.提出的法向?qū)?shù)插值算法具有復(fù)雜度低和計算量小的顯著特點.通過補(bǔ)充定義分段光滑曲線在尖點處的廣義導(dǎo)數(shù)（正則導(dǎo)數(shù)或單側(cè)連續(xù)導(dǎo)數(shù)），可將法向?qū)?shù)插值算法拓展到分段光滑曲線邊界情形.該算法彌補(bǔ)了矩形網(wǎng)格有限差分法處理復(fù)雜邊界問題的缺陷.模型試算結(jié)果驗證了該算法的有效性.

根據(jù)波場快照中各個波前到達(dá)檢波器排列位置所對應(yīng)的快照時間，可以準(zhǔn)確地分析與識別地震記錄中的各個同相軸.利用波場系列快照技術(shù)，可揭示地震波場在復(fù)雜邊界（包括地表自由邊界和內(nèi)部邊界）控制下的傳播規(guī)律.

Alford R M，Kelly K R，Boore D M.1974.Accuracy of finitedifference modeling of the acoustic wave equation.Geophysics，39（6）：834-842.

Appel?D，Petersson N A.2009.A stable Finite difference method for the elastic wave equation on complex geometries with free surfaces.Commun.Comput.Phys.，5（1）：84-107.

Chu C L，Wang X T.2005.Seismic Modeling with a Finite-Difference Method on Irregular Triangular Grids.Periodical of Ocean University of China（in Chinese），35（1）：43-48.

Dong L G，Guo X L，Wu X F，etal.2007.Finite difference numerical simulation for the elastic wave propagation in rugged topography.Natural Gas Industry（in Chinese），27（10）：38-41.

Fu L Y，Wu R S.2000.Infinite boundary element absorbing boundary for wave propagation simulations.Geophysics，65（2）：596-602.

Ge Z X，Chen X F.2008.An efficient approach for simulating wave propagation with the boundary element method in multilayered media with irregular interfaces.Bulletin of the Seismological Society of America，98（6）：3007-3016.

Guan X Z，F(xiàn)u L Y，Tao Y，etal.2011.Boundary-volume Integral equation numerical modeling for complex near surface.Chinese J.Geophys.（in Chinese），54（9）：2357-2367.

Huang Z P，Zhang M，Wu W Q，etal.2004.A domain decomposition method for numerical simulation of the elastic wave propagation.Chinese J.Geophys.（in Chinese），47（6）：1094-1100.

Jih R S，Mclaughlin K L，Der Z A.1988.Free-boundary conditions of arbitrary polygonal topography in a two-dimensional explicit elastic finite-difference scheme.Geophysics，53（8）：1045-1055.

Lan H Q，Zhang Z J.2011.Three-dimensional wave-field simulation in heterogeneous transversely isotropic medium with irregular free surface.Bulletin of the Seismological Society of America，101（3）：1354-1370.

Li X X，Yu G X，F(xiàn)u L Y，etal.2011.Analysing seismic scattering characteristics of complex seabed by using the boundary-element simulation method.China Offshore Oil and Gas（in Chinese），23（6）：357-361，368.

Ma D T，Zhu G M.2004.Hybrid method combining finite difference and pseudospectral method for solving the elastic wave equation.Journal of Earth Sciences and Environment（in Chinese），26（1）：61-64.

Ma X，Yang D H，and Liu F Q.2011.A nearly analytic symplectically partitioned Runge-Kutta method for 2-D elastic wave equations.Geophysical Journal International，187（1）：480-496.

Pei Z L.2004.Numerical modeling using staggered-grid high-order finite-difference of elastic wave equation on arbitrary relief surface.Oil Geophysical Prospecting（in Chinese），39（6）：629-634.

Pérez-Ruiz J A，Luzón F，García-Jerez A.2005.Simulation of an irregular free surface with a displacement finite-difference scheme.Bulletin of the Seismological Society of America，95（6）：2216-2231.

Sun X D，Li Z C，Wang X L.2012.Pre-stack reverse-time migration using a finite difference method based on triangular grids.Progress in Geophysics（in Chinese），27（5）：2077-2083.>

Tong P，Yang D H，Hua B L，etal.2013.A high-order stereomodeling method for solving wave equations.Bulletin of the Seismological Society of America，103（2A）：811-833.

Wang D L，He Q D，Han L G.2005.Multi-azimuth threecomponent surface seismic modeling for cracked monoclinic media.Chinese J.Geophys.（in Chinese），48（2）：386-393.

Xue D C，Wang S X，Jiao S J.2007.Wave equation finite-element modeling including rugged topography and complicated medium.Progress in Geophysics（in Chinese），22（2）：522-529.

Yang D H.2002.Finite element method of the elastic wave equation and wavefield simulation in two-phase anisotropic media.Chinese J.Geophys.（in Chinese），45（4）：575-583.

Yang D H，Tong P，Deng X Y.2012a.A central difference method with low numerical dispersion for solving the scalar wave equation.Geophysical Prospecting，60（5）：885-905.

Yang D H，Wang N，Liu F Q.2012b.A strong stability-preserving predictor-corrector method for the simulation of elastic wave propagation in anisotropic media.Communications in Computational Physics，12（4）：1006-1032.

Zheng H S，Zhang Z J.2005.Synthetic seismograms of nonlinear seismic waves in anisotropic（VTI）media.Chinese J.Geophys.（in Chinese），48（3）：660-671.

附中文參考文獻(xiàn)

褚春雷，王修田.2005.非規(guī)則三角網(wǎng)格有限差分法地震正演模擬.中國海洋大學(xué)學(xué)報，35（1）：43-48.

董良國，郭曉玲，吳曉豐等.2007.起伏地表彈性波傳播有限差分法數(shù)值模擬.天然氣工業(yè)，27（10）：38-41.

管西竹，符力耘，陶毅等.2011.復(fù)雜地表邊界元-體積元波動方程數(shù)值模擬.地球物理學(xué)報，54（9）：2357-2367.

黃自萍，張銘，吳文青等.2004.彈性波傳播數(shù)值模擬的區(qū)域分裂法.地球物理學(xué)報，47（6）：1094-1100.

李緒宣，于更新，符力耘等.2011.應(yīng)用邊界元模擬方法分析復(fù)雜海底地震散射特征.中國海上油氣，23（6）：357-361，368.

馬德堂，朱光明.2004.有限元法與偽譜法混合求解彈性波動方程.地球科學(xué)與環(huán)境學(xué)報，26（1）：61-64.

裴正林.2004.任意起伏地表彈性波方程交錯網(wǎng)格高階有限差分法數(shù)值模擬.石油地球物理勘探，39（6）：629-634.

孫小東，李振春，王小六.2012.三角網(wǎng)格有限差分法疊前逆時偏移方法研究.地球物理學(xué)進(jìn)展，27（5）：2077-2083.

王德利，何樵登，韓立國.2005.裂隙型單斜介質(zhì)中多方位地面三分量記錄模擬.地球物理學(xué)報，48（2）：386-393.

薛東川，王尚旭，焦淑靜.2007.起伏地表復(fù)雜介質(zhì)波動方程有限元數(shù)值模擬方法.地球物理學(xué)進(jìn)展，22（2）：522-529.

楊頂輝.2002.雙相各向異性介質(zhì)中彈性波方程的有限元解法及波場模擬.地球物理學(xué)報，45（4）：575-583.

鄭海山，張中杰.2005.橫向各向同性（VTI）介質(zhì)中非線性地震波場模擬.地球物理學(xué)報，48（3）：660-671.