朱　強(qiáng),姜蘆倩,張　偉

(1.中國(guó)科學(xué)技術(shù)大學(xué)地球和空間科學(xué)學(xué)院,安徽合肥230026;2.南方科技大學(xué)地球與空間科學(xué)系,廣東深圳518055)

地震波場(chǎng)數(shù)值模擬是研究地震波傳播規(guī)律的主要手段,也是基于波動(dòng)方程進(jìn)行地震反演和成像的基礎(chǔ)。有限差分方法具有實(shí)現(xiàn)簡(jiǎn)單、應(yīng)用范圍廣、內(nèi)存需求小、易于并行的優(yōu)點(diǎn),被廣泛應(yīng)用于理論地震學(xué)和勘探地震研究中。目前的有限差分算法可以根據(jù)變量在計(jì)算網(wǎng)格上的分布不同而分成同位網(wǎng)格和交錯(cuò)網(wǎng)格兩種網(wǎng)格系統(tǒng)。交錯(cuò)網(wǎng)格布局又可以根據(jù)波場(chǎng)分量在計(jì)算網(wǎng)格上的相對(duì)位置關(guān)系被大致分為三種格式:標(biāo)準(zhǔn)交錯(cuò)網(wǎng)格[1-2]、完全交錯(cuò)網(wǎng)格[3-4]和旋轉(zhuǎn)交錯(cuò)網(wǎng)格[5]。交錯(cuò)網(wǎng)格有限差分方法通常采用中心差分的方式求解一階速度-應(yīng)力方程組,可以避免同位網(wǎng)格簡(jiǎn)單中心差分格式的奇偶失聯(lián)問(wèn)題[6],是目前求解波動(dòng)方程數(shù)值解的主流方法之一。

有限差分波動(dòng)方程數(shù)值解存在不同類型的誤差,其中與地震波頻率有關(guān)的色散/耗散誤差可以通過(guò)使用高階差分算子或者優(yōu)化格式來(lái)壓制[7-10]。而包含內(nèi)部速度間斷面的復(fù)雜速度模型采用有限差分網(wǎng)格離散時(shí),由于網(wǎng)格與速度界面不重疊,同一物理模型采用不同的介質(zhì)離散方法,因此會(huì)導(dǎo)致不同的網(wǎng)格離散模型,進(jìn)而引起正演結(jié)果的不同。ALTERMAN等[11]采用虛擬網(wǎng)格點(diǎn)和顯式實(shí)施內(nèi)部邊界條件的方式處理二階位移方程中的水平間斷面問(wèn)題。VIRIEUX[1,12]發(fā)展了基于一階速度-應(yīng)力方程的交錯(cuò)網(wǎng)格有限差分方法,通過(guò)網(wǎng)格點(diǎn)離散取值表征速度模型,并且展示了這種方法在大泊松比情況下的穩(wěn)定性。GRAVES[13]給出了交錯(cuò)網(wǎng)格有限差分算法處理三維彈性介質(zhì)參數(shù)的格點(diǎn)平均方法。MUIR等[14]首次將等效介質(zhì)參數(shù)理論[15]應(yīng)用到有限差分方法中。MOCZO等[16]提出了體積分平均法,用來(lái)處理交錯(cuò)網(wǎng)格有限差分算法中的介質(zhì)參數(shù)間斷面,解決了物理界面在網(wǎng)格格點(diǎn)之間變化時(shí)的表征問(wèn)題。MOCZO等[17]、KRISTEK等[18]將離散后界面附近的介質(zhì)近似為正交各向異性介質(zhì)來(lái)實(shí)施間斷面的邊界條件,從而有效地壓制了界面誤差。LISITSA等[19]通過(guò)平面波方程分析了交錯(cuò)網(wǎng)格有限差分方法處理界面的理論誤差。VISHNEVSKY等[20]對(duì)標(biāo)準(zhǔn)交錯(cuò)網(wǎng)格、完全交錯(cuò)網(wǎng)格和旋轉(zhuǎn)交錯(cuò)網(wǎng)格在層狀界面下的誤差收斂情況進(jìn)行了總結(jié)。但是,上述這些研究都沒(méi)有定量評(píng)估和給出在任意復(fù)雜模型中準(zhǔn)確計(jì)算(滿足給定誤差要求)反射/透射波的網(wǎng)格大小要求,尤其是沒(méi)有考慮準(zhǔn)確計(jì)算P-S轉(zhuǎn)換波所要求的網(wǎng)格大小,所以目前進(jìn)行有限差分地震模擬通常僅考慮格式本身的色散和耗散誤差對(duì)網(wǎng)格大小的限制,而沒(méi)有考慮介質(zhì)離散方法不同導(dǎo)致的與介質(zhì)界面有關(guān)的反射/透射波計(jì)算準(zhǔn)確性的影響。

本文首先介紹了交錯(cuò)網(wǎng)格有限差分方法原理,然后總結(jié)了目前有限差分地震模擬中所使用的介質(zhì)離散方法,最后定量分析了不同介質(zhì)離散方法對(duì)不同地震波場(chǎng)模擬準(zhǔn)確性的影響。

1　交錯(cuò)網(wǎng)格有限差分方法

在二維笛卡爾坐標(biāo)系下,波動(dòng)方程的速度-應(yīng)力形式表示為:

式中:vx與vz代表速度分量;τxx,τzz與τxz代表應(yīng)力分量;ρ表示介質(zhì)密度;t表示時(shí)間;Cij(i=1,2,3;j=1,2,3)為彈性剛度矩陣,各向同性介質(zhì)情況下,Cij由兩個(gè)獨(dú)立的分量構(gòu)成,即拉梅常數(shù)λ與μ。

為了分析問(wèn)題方便,假定計(jì)算網(wǎng)格是矩形,顯式M(M=2,4,6,8,…)階交錯(cuò)網(wǎng)格差分算子Dx定義為:

(3)

式中:h表示網(wǎng)格步長(zhǎng),αm表示差分系數(shù)。則交錯(cuò)網(wǎng)格格式如下:

(4)

2　介質(zhì)參數(shù)離散方法

對(duì)于介質(zhì)參數(shù)不連續(xù)的內(nèi)部界面,其邊界條件應(yīng)該滿足位移連續(xù)與法向牽引力連續(xù)的條件。目前較為常用的介質(zhì)參數(shù)離散方法有直接取值法、格點(diǎn)平均法[13]、體積分平均法[16]與正交各向異性等效介質(zhì)法[17-18]。本文以一個(gè)傾斜界面模型為例,說(shuō)明不同介質(zhì)離散方法的差異(圖1)。如圖1a所示,不同介質(zhì)的分界面傾斜穿過(guò)網(wǎng)格,傾角為30°,上層介質(zhì)參數(shù)為vP1=2000m/s,vS1=1410m/s,ρ1=1160kg/m3;下層介質(zhì)參數(shù)為vP2=4000m/s,vS2=2310m/s,ρ2=2050kg/m3。圖1b 是采用廣義反射-透射系數(shù)方法[21-23]按圖1a中的觀測(cè)系統(tǒng)計(jì)算的接收波形(通過(guò)水平層計(jì)算后旋轉(zhuǎn)得到),震源為炸藥震源。其中第一個(gè)事件(第一個(gè)檢波器0～1.0s之間)是直達(dá)P波,第二個(gè)事件(第一個(gè)檢波器1.5～2.0s之間)是反射P波,第三個(gè)事件是P-S轉(zhuǎn)換波(第一個(gè)檢波器2.3～2.6s之間)。以下介紹不同的介質(zhì)離散方法、圖1a所示模型用不同介質(zhì)離散方法得到的不同離散模型,以及與圖1b對(duì)應(yīng)的不同事件計(jì)算結(jié)果的差異(圖2,圖3)。

圖1　斜層模型a 模型示意; b 道集

圖2　修正介質(zhì)參數(shù)示意a 直接取值法; b 格點(diǎn)平均法; c 體積分平均法; d 正交各向異性等效介質(zhì)法

圖3　斜層模型局部道集展示(圖1b中2～3s波形)a 直接取值法; b 格點(diǎn)平均法; c 體積分平均法; d 正交各向異性等效介質(zhì)法

2.1　直接取值法

直接取值法即直接將物理模型根據(jù)計(jì)算域劃分的網(wǎng)格予以離散,每個(gè)格點(diǎn)取其與物理模型相對(duì)應(yīng)點(diǎn)的介質(zhì)參數(shù)。界面傾斜穿過(guò)網(wǎng)格時(shí),該方法會(huì)形成階梯狀的介質(zhì)參數(shù)間斷界面,見(jiàn)圖2a。根據(jù)惠更斯原理,此時(shí)階梯上的每一個(gè)點(diǎn)相當(dāng)于一個(gè)點(diǎn)震源,從而導(dǎo)致界面上產(chǎn)生人工的虛假波。圖3a展示了采用直接取值法離散介質(zhì)所產(chǎn)生的散射波形,可以看到在每道地震記錄的P-S轉(zhuǎn)換波后都伴隨著明顯的波形扭曲(第一個(gè)檢波器2.6s之后)。因此直接取值法不能精確模擬界面的實(shí)際位置,誤差較為明顯。

2.2　格點(diǎn)平均法

格點(diǎn)平均法對(duì)界面附近的介質(zhì)做了簡(jiǎn)單的平均化處理,避免了介質(zhì)參數(shù)的劇烈間斷(圖2b),但這種方法本身無(wú)法對(duì)單位網(wǎng)格之內(nèi)的界面變化進(jìn)行細(xì)致的描述,因此當(dāng)界面位置與網(wǎng)格位置沒(méi)有耦合的時(shí)候,其模擬精度較差(圖3b)。

2.3　體積分平均法

當(dāng)界面在單位網(wǎng)格間變化的時(shí)候,簡(jiǎn)單的格點(diǎn)平均法無(wú)法準(zhǔn)確反映界面位置。MOCZO等[16]提出的體積分平均法可以更好地描述界面位置,提高反射/透射波的模擬精度。在體積分平均法中,某個(gè)點(diǎn)在計(jì)算域的介質(zhì)參數(shù)不僅僅由這個(gè)點(diǎn)本身所在位置的介質(zhì)參數(shù)決定,還受到以它為中心的單位網(wǎng)格區(qū)域內(nèi)的介質(zhì)參數(shù)的影響。

設(shè)界面位于x=0處,在一維情況下對(duì)波動(dòng)方程在單位網(wǎng)格內(nèi)做空間積分:

(7)

對(duì)(7)式應(yīng)用均值定理,得到:

(8)

則可以定義密度分量ρ的積分平均為:

(9)

ρA反映了這個(gè)被積分單元的平均響應(yīng)。將其推廣到二維情況下,同理有:

如公式(10)至公式(12)所示,對(duì)于網(wǎng)格上的某一點(diǎn),應(yīng)用體積分平均法對(duì)其周圍單位空間的介質(zhì)參數(shù)進(jìn)行積分然后取平均,即可得到修正后的彈性介質(zhì)參數(shù)。圖2c展示了體積分平均法處理后的介質(zhì)參數(shù),與直接取值法和格點(diǎn)平均法相比,這種方法能夠反映界面在單位網(wǎng)格間的具體變化。同時(shí),體積分平均法對(duì)傾斜界面做了合適的處理,與直接取值法相比有效壓制了因傾斜界面介質(zhì)離散不準(zhǔn)確而產(chǎn)生的散射波(圖3c)。

2.4　正交各向異性等效介質(zhì)法

層狀介質(zhì)研究表明,由五個(gè)彈性系數(shù)描述的均勻橫向各向同性介質(zhì)是層狀各向同性介質(zhì)疊加物的長(zhǎng)波等效[15,24]。其等效介質(zhì)的彈性參數(shù)可以通過(guò)對(duì)層狀介質(zhì)拉梅常數(shù)的合理平均來(lái)近似獲得。正交各向異性等效介質(zhì)法基于這一等效理論,通過(guò)將介質(zhì)界面附近的彈性參數(shù)等效為一種正交各向異性介質(zhì)彈性參數(shù),模擬界面在單位網(wǎng)格內(nèi)的不同位置。

二維各向同性介質(zhì)情況下:

(13)

式中:εxx,εxz與εzz表示應(yīng)變分量。將(13)式應(yīng)力-應(yīng)變關(guān)系應(yīng)用于一個(gè)水平介質(zhì)界面,τzz,τxz與εxx在界面兩端連續(xù),τxx,εxz與εzz在界面兩端不連續(xù),用上標(biāo)“+”表示在界面兩端具有連續(xù)性的分量,上標(biāo)“-”表示在界面兩端間斷的分量,則(13)式可以改寫為:

(14)

其中,

(15)

(14)式中彈性剛度矩陣被橫線與豎線分為4個(gè)部分,將每一個(gè)部分看作一個(gè)整體變量,解出界面兩端間斷的分量τ-與ε-,同時(shí)對(duì)等式兩端取平均可得:

(16)

式中,“〈〉”表示對(duì)單位網(wǎng)格求積分平均。(16)式反映了位于界面上的微小單元的應(yīng)力-應(yīng)變關(guān)系,由于〈τ+〉=τi,j,〈ε+〉=εi,j,因此可改寫為:

(17)

(18)

設(shè)M=λ+2μ,則:

(19)

式中:〈〉H表示調(diào)和平均,通過(guò)這樣的平均化處理將界面附近的介質(zhì)視為等效的正交各向異性介質(zhì),將其彈性剛度矩陣代入應(yīng)力-應(yīng)變關(guān)系以模擬界面上的邊界條件。與體積分平均法一樣,正交各向異性等效介質(zhì)法能夠準(zhǔn)確地描述界面在單位網(wǎng)格之間的實(shí)際位置(圖2d),從而能夠獲得更加準(zhǔn)確的地震記錄(圖3d)。

3　不同介質(zhì)離散方法對(duì)不同事件模擬準(zhǔn)確性的影響

本文通過(guò)改變計(jì)算域網(wǎng)格的每波長(zhǎng)網(wǎng)格數(shù)(points per wavelength,PPW)來(lái)觀察界面誤差的變化情況,分析了當(dāng)每波長(zhǎng)網(wǎng)格數(shù)逐漸增大時(shí),由界面計(jì)算產(chǎn)生的不同反射/透射波的相對(duì)誤差。為了定量評(píng)估交錯(cuò)網(wǎng)格有限差分算法對(duì)界面計(jì)算的準(zhǔn)確性,首先定義相對(duì)誤差為:

(20)

式中:‖u‖2表示u的二范數(shù),ufd表示有限差分解,uref表示參考解。本文通過(guò)截取相應(yīng)事件窗口的波形計(jì)算相對(duì)誤差,使用PPW來(lái)評(píng)估網(wǎng)格的分辨率。對(duì)

于Ricker子波,其PPW定義為:

(21)

式中:fc為Ricker子波中心頻率,Δh為空間網(wǎng)格步長(zhǎng),vmin表示P波或S波的最小速度。

在具體探索各種介質(zhì)參數(shù)處理方法在不同波形下的相對(duì)誤差PPW約束之前,首先觀察不同反射/透射角度對(duì)相對(duì)誤差的影響。如圖4所示設(shè)置一個(gè)兩層模型,上層介質(zhì)參數(shù)為vP1=2000m/s,vS1=1000m/s,ρ1=1200kg/m3,下層介質(zhì)參數(shù)為vP2=4000m/s,vS2=2000m/s,ρ2=2100kg/m3,震源位于(1000m,1500m)處。在距離界面上方400m與下方300m深度,分別沿著圖4中的虛線布置兩排檢波器,相對(duì)誤差隨不同反射/透射角的變化如圖5所示。由于某些位置反射/透射波與P-S轉(zhuǎn)換波發(fā)生了耦合,因此在展示相對(duì)誤差隨反射/透射角度的變化時(shí),使用整體波形計(jì)算其相對(duì)誤差。

圖4　不同PPW相對(duì)誤差測(cè)試模型

3.1　透射波和反射波相對(duì)誤差隨網(wǎng)格的變化

在圖5中相對(duì)誤差較大的位置附近,同時(shí)在時(shí)間上能夠區(qū)分反射/透射波與P-S轉(zhuǎn)換波的前提下,將反射波接收點(diǎn)置于(1930m,2060m)處,透射波接收點(diǎn)置于(2260m,3160m)處。通過(guò)保持震源、檢波器與界面相對(duì)位置不變,逐漸縮小網(wǎng)格空間步長(zhǎng)以得到不同的PPW。為了使獲得的結(jié)論更加一般化,將實(shí)際物理界面位置置于兩個(gè)網(wǎng)格之間(不與任何網(wǎng)格格點(diǎn)發(fā)生耦合),使用主頻為10Hz的Ricker子波作為震源子波,分別對(duì)比不同反射/透射波在使用不同介質(zhì)離散方法后的相對(duì)誤差。為了避免色散誤差的影響,本文采用基于泰勒展開(kāi)的10階交錯(cuò)網(wǎng)格格式,該格式滿足格式的色散誤差要求的PPW為5左右,而下面測(cè)試的PPW最小是6,因此所有測(cè)試都滿足格式的色散誤差要求。

圖5　相對(duì)誤差隨偏移距的變化a 反射波; b 透射波

圖6分別展示了反射波與透射波相對(duì)誤差隨PPW的變化,可以看到不同介質(zhì)離散方法對(duì)P波影響很大。直接取值法具有極大的不穩(wěn)定性,其對(duì)于波形的誤差顯著大于其它方法。就不同波形而言,相同情況下透射波的準(zhǔn)確性普遍高于反射波。圖7展示了反射波與透射波分別在相對(duì)誤差約為10%與5%時(shí)的波形對(duì)比。

3.2　P-S轉(zhuǎn)換波誤差隨網(wǎng)格的變化

圖8展示了P-S轉(zhuǎn)換波相對(duì)誤差隨PPW的變化。由于界面上發(fā)生了P-S波的轉(zhuǎn)換,因此轉(zhuǎn)換波(反射)的精度顯著低于其它波形。由于直接取值法與格點(diǎn)平均法不能精確描述界面在網(wǎng)格間的具體位置,因此它們?cè)诖志W(wǎng)格情況下誤差較大,而通過(guò)增大PPW提高模擬精度的計(jì)算代價(jià)高昂。P-S轉(zhuǎn)換波在相對(duì)誤差約為10%與5%時(shí)的波形對(duì)比見(jiàn)圖9。

圖6　相對(duì)誤差隨PPW的變化a 反射波; b 透射波

圖7　相對(duì)誤差波形對(duì)比(10%與5%)a 反射波波形; b 透射波波形

圖8　P-S轉(zhuǎn)換波相對(duì)誤差隨PPW的變化a 反射波; b 透射波

圖9　P-S轉(zhuǎn)換波相對(duì)誤差波形對(duì)比(10%和5%)a 反射波波形; b 透射波波形

3.3　多次反射波誤差隨網(wǎng)格的變化

為了觀察多次反射波相對(duì)誤差隨網(wǎng)格的變化,設(shè)置一個(gè)平層模型,模型大小為nx×nz=2500m×1500m,包括3層介質(zhì):頂層與底層縱波速度為4000m/s,橫波速度為2000m/s,密度2050kg/m3;中間低速層縱波速度為2000m/s,橫波速度為1000m/s,密度1410kg/m3。震源(700m,800m)與檢波點(diǎn)(1330m,770m)均位于模型中間的低速層。圖10展示了多次反射波在不同PPW情況下,使用不同介質(zhì)離散方法得到的相對(duì)誤差。圖11和圖12展示了誤差約10%和5%時(shí)的波形對(duì)比(實(shí)際誤差分別為9.61%與4.79%),可以看到在相對(duì)誤差約為5%時(shí),波形的擬合已經(jīng)相對(duì)較好。

圖10　多次反射波相對(duì)誤差隨PPW的變化

圖11　多次反射波相對(duì)誤差波形對(duì)比

表1與表2分別展示了4種介質(zhì)離散方法達(dá)到10%與5%的相對(duì)誤差時(shí),不同類型反射/透射波各自的PPW要求。

圖12　多次反射波相對(duì)誤差波形對(duì)比局部(圖11中1.0～1.6s波形)

個(gè)

表2　不同反射/透射波PPW約束(相對(duì)誤差5%)　個(gè)

4　結(jié)束語(yǔ)

本文針對(duì)介質(zhì)的彈性系數(shù)處理對(duì)于交錯(cuò)網(wǎng)格波場(chǎng)模擬準(zhǔn)確性造成的影響進(jìn)行了分析,給出了不同介質(zhì)離散方式下各個(gè)反射/透射波相對(duì)誤差的PPW約束參考,得到以下結(jié)論：

1) 在使用交錯(cuò)網(wǎng)格有限差分進(jìn)行正演計(jì)算時(shí),不僅要考慮格式本身的色散誤差,還需要考慮介質(zhì)離散方式對(duì)模擬精度的影響,特別是在復(fù)雜模型中,即使在滿足色散誤差的情況下,仍然存在較為嚴(yán)重的界面誤差,這類誤差對(duì)整體波形有著不可忽視的影響。

2) 界面誤差對(duì)PPW的要求普遍大于4階格式的色散誤差,因此,在選取交錯(cuò)網(wǎng)格的網(wǎng)格步長(zhǎng)時(shí),為了保證誤差在一定范圍之內(nèi),應(yīng)該以橫波最小速度對(duì)應(yīng)的PPW為約束,在可以接受的誤差范圍內(nèi)首先確定網(wǎng)格的最大空間步長(zhǎng)。

在實(shí)際交錯(cuò)網(wǎng)格計(jì)算中,首先應(yīng)該確定計(jì)算精度要求,根據(jù)所采取的介質(zhì)參數(shù)處理算法,參考界面誤差隨PPW的變化關(guān)系,以最小橫波速度為基準(zhǔn),確定網(wǎng)格的最大空間步長(zhǎng),再根據(jù)柯朗-弗里德里希斯-列維(CFL)穩(wěn)定性條件確定時(shí)間步長(zhǎng),從而在確保計(jì)算準(zhǔn)確性的前提下,減少計(jì)算資源需求。

[1]VIRIEUX J.P-SV wave propagation in heterogeneous media:velocity-stress finite-difference method[J].Geophysics,1986,51(4):889-901

[2]LEVANDER A R.Fourth-order finite-difference P-SV seismograms[J].Geophysics,1988,53(11):1425-1436

[3]MCGARRY R,PASALIC D,ONG C.Anisotropic elastic modeling on a Lebedev grid:dispersion reduction and grid decoupling[J].Expanded Abstracts of 81stAnnual Internat SEG Mtg,2011:2829-2833

[4]LISITSA V,VISHNEVSKIY D.Lebedev scheme for the numerical simulation of wave propagation in 3D anisotropic elasticity[J].Geophysical Prospecting,2010,58(4):619-635

[5]SAENGER E H,GOLD N,SHAPIRO S A.Modeling the propagation of elastic waves using a modified finite-difference grid[J].Wave Motion,2000,31(1):77-92

[6]裴正林.三維各向同性介質(zhì)彈性波方程交錯(cuò)網(wǎng)格高階有限差分法模擬[J].石油物探,2005,44(4):308-315

PEI Z L.Numerical simulation of elastic wave equation in 3-D isotropic media with staggered-grid high-order difference method[J].Geophysical Prospecting for Petroleum,2005,44(4):308-315

[7]李勝軍,孫成禹,高建虎,等.地震波數(shù)值模擬中的頻散壓制方法分析[J].石油物探,2008,47(5):444-448

LI S J,SUN C Y,GAO J H,et al.Analysis of dispersion suppression in wave equation numerical simulation[J].Geophysical Prospecting for Petroleum,2008,47(5):444-448

[8]陳可洋.地震波數(shù)值模擬中差分近似的各向異性分析[J].石油物探,2010,49(1):19-22

CHEN K Y.Anisotropic analysis of difference approximation in seismic wave numerical modeling[J].Geophysical Prospecting for Petroleum,2010,49(1):19-22

[9]梁文全,王彥飛,楊長(zhǎng)春.基于優(yōu)化方法的時(shí)間-空間域隱格式有限差分算子確定方法[J].石油物探,2015,54(3):254-259

LIANG W Q,WANG Y F,YANG C C.Determination on the implicit finite-difference operator based on optimization method in time-space domain[J].Geophysical Prospecting for Petroleum,2015,54(3):254-259

[10]陳東,梁文全,辛維.適用于聲波方程數(shù)值模擬的時(shí)間-空間域隱式有限差分算子優(yōu)化方法[J].地球物理學(xué)報(bào),2016,59(4):1491-1497

CHEN D,LIANG W Q,XIN W.Acoustic wave equation modeling based on implicit finite difference operators in the time-space domain[J].Chinese Journal of Geophysics,59(4):1491-1497

[11]ALTERMAN Z,KARAL F C.Propagation of elastic waves in layered media by finite difference methods[J].Bulletin of the Seismological Society of America,1968,58(1):367-398

[12]VIRIEUX J.SH-wave propagation in heterogeneous media:velocity-stress finite-difference method[J].Geophysics,1984,49(11):1933-1942

[13]GRAVES R W.Simulating seismic wave propagation in 3D elastic media using staggered-grid finite differences[J].Bulletin of the Seismological Society of America,1996,86(4):1091-1106

[14]MUIR F,DELLINGER J,ETGEN J,et al.Modeling elastic fields across irregular boundaries[J].Geophysics,1992,57(9):1189-1193

[15]BACKUS G E.Long-wave elastic anisotropy produced by horizontal layering[J].Journal of Geophysical Research,1962,67(11):4427-4440

[17]MOCZO P,KRISTEK J,GLIS M.The finite-difference modelling of earthquake motions:waves and ruptures[M].Cambridge:Cambridge University Press,2014:199-226

[18]KRISTEK J,MOCZO P,CHALJUB E,et al.An orthorhombic representation of a heterogeneous medium for the finite-difference modelling of seismic wave propagation[J].Geophysical Journal International,2017,208(2):1250-1264

[19]LISITSA V,PODGORNOVA O,TCHEVERDA V.On the interface error analysis for finite difference wave simulation[J].Computational Geosciences,2010,14(4):769-778

[20]VISHNEVSKY D,LISITSA V,TCHEVERDA V,et al.Numerical study of the interface errors of finite-difference simulations of seismic waves[J].Geophysics,2014,79(4):T219-T232

[21]CHEN X.Seismogram synthesis for multi-layered media with irregular interfaces by global generalized reflection/transmission matrices method I,theory of two-dimensional SH case[J].Bulletin of the Seismological Society of America,1990,80(6A):1696-1724

[22]CHEN X.Seismogram synthesis for multi-layered media with irregular interfaces by global generalized reflection/transmission matrices method Ⅱ,applications for 2D SH case[J].Bulletin of the Seismological Society of America,1995,85(4):1094-1106

[23]CHEN X.Seismogram synthesis for multi-layered media with irregular interfaces by global generalized reflection/transmission matrices method III,theory of 2D P-SV case[J].Bulletin of the Seismological Society of America,1996,86(2):389-405

[24]SCHOENBERG M,MUIR F.A calculus for finely layered anisotropic media[J].Geophysics,1989,54(5):581-589