起伏層狀介質(zhì)中曲線網(wǎng)格有限差分法與射線法波場模擬對(duì)比研究*

楊尚倍 白超英 何雷宇

1) 中國西安710054長安大學(xué)地質(zhì)工程與測繪學(xué)院2) 中國西安710054長安大學(xué)計(jì)算地球物理研究所

?

起伏層狀介質(zhì)中曲線網(wǎng)格有限差分法與射線法波場模擬對(duì)比研究*

楊尚倍1)白超英1,2),*何雷宇1)

1) 中國西安710054長安大學(xué)地質(zhì)工程與測繪學(xué)院2) 中國西安710054長安大學(xué)計(jì)算地球物理研究所

針對(duì)處理起伏地表(或含地下不規(guī)則波阻抗界面)條件下發(fā)展起來的地震波場數(shù)值模擬算法的模擬結(jié)果與解析解(大多數(shù)情形下無法得到)無法進(jìn)行對(duì)比， 且其有效性和正確性難以驗(yàn)證的情況， 本文提出了一種可以相互驗(yàn)證波場數(shù)值模擬結(jié)果與射線追蹤數(shù)值模擬結(jié)果的正確性和有效性的佐證方法， 驗(yàn)證了參考射線追蹤法. 其中， 波場數(shù)值模擬中采用曲線網(wǎng)格DRP/opt MacCormack有限差分法， 射線追蹤模擬則采用分區(qū)多步三角網(wǎng)格最短路徑算法． 通過系統(tǒng)對(duì)比上述兩種方法得到的波場快照、 單炮地震記錄， 以及合成理論地震圖的結(jié)果顯示， 本方法相互作證了兩種方法所得結(jié)果的正確性和有效性. 雙層和三層起伏層狀模型的對(duì)比分析結(jié)果表明， 這種方法不但可以加深理解地震波在復(fù)雜介質(zhì)中的傳播規(guī)律， 同時(shí)射線法的引入為清晰識(shí)別和標(biāo)定地震波場數(shù)值模擬中各種不同震相提供了一種便捷的途徑.

分區(qū)多步三角網(wǎng)格最短路徑法 曲線網(wǎng)格有限差分法 地震波傳播 數(shù)值模擬 起伏層狀介質(zhì)

引言

地震波場數(shù)值模擬技術(shù)是指假定地下介質(zhì)和物性參數(shù)已知， 模擬地震波在地下介質(zhì)中的傳播規(guī)律， 并計(jì)算地面或地下各觀測點(diǎn)所觀測到的數(shù)值地震記錄的一種地震模擬方法(Carcioneetal， 2002)． 該方法是一種先將特定的地質(zhì)、 地球物理問題作合理簡化后抽象形成數(shù)學(xué)模型， 再采用數(shù)值計(jì)算方法來獲取地震響應(yīng)的過程． 地震波場數(shù)值模擬的理論基礎(chǔ)主要有射線理論和波動(dòng)方程理論(Carcioneetal, 2002; 朱多林， 白超英， 2011)．

因此， 深入研究這兩種數(shù)值模擬算法的共性和差異， 將為分析地震波在復(fù)雜介質(zhì)中傳播的細(xì)節(jié)提供理論依據(jù)， 同時(shí)射線理論又是波動(dòng)方程高頻近似下的理論， 所以二者在高頻下應(yīng)有相同的模擬結(jié)果， 兩者本應(yīng)相互佐證和共同發(fā)展完善， 但是實(shí)際上卻是平行獨(dú)立的． 王濤等(2014， 2016)分別對(duì)各向同性、 各向異性介質(zhì)中射線理論與高階交錯(cuò)網(wǎng)格有限差分波場數(shù)值模擬的結(jié)果作了系統(tǒng)的對(duì)比研究； 作為補(bǔ)充和完善， 本文擬對(duì)模型較為復(fù)雜的情形， 即地表起伏(含地下不規(guī)則波阻抗界面)時(shí)， 采用曲線網(wǎng)格DRP/opt MacCormack有限差分法和三角網(wǎng)最短路徑射線法， 通過波場快照、 單炮地震記錄、 合成理論地震圖的系統(tǒng)對(duì)比， 以期相互驗(yàn)證兩種不同數(shù)值模擬方法的有效性和正確性．

1 兩種數(shù)值模擬算法

1.1 分區(qū)多步三角網(wǎng)格最短路徑射線追蹤算法基本原理

采用規(guī)則網(wǎng)格劃分介質(zhì)模型時(shí)(Moser， 1991; 趙愛華， 張中杰， 2004; 趙愛華， 丁志峰， 2005; 王童奎等， 2007)， 為了模擬起伏界面和提高計(jì)算精度， 一般通過不斷縮小網(wǎng)格節(jié)點(diǎn)間距或加密單元內(nèi)次生節(jié)點(diǎn)的個(gè)數(shù)來實(shí)現(xiàn)， 這種方法會(huì)導(dǎo)致網(wǎng)格節(jié)點(diǎn)或單元內(nèi)次生節(jié)點(diǎn)總數(shù)過多， 造成計(jì)算量過大和對(duì)于復(fù)雜模型適應(yīng)性差等缺點(diǎn)． 鑒于此， 本文采用三角網(wǎng)格單元剖分下的分區(qū)多步最短路徑射線追蹤算法(multistage triangular shortest-path method， 簡寫為MTSPM)(Baietal， 2012)， 來實(shí)現(xiàn)二維各向同性復(fù)雜層狀介質(zhì)中的初至波、 反射和轉(zhuǎn)換波的射線追蹤．

1.1.1 模型分區(qū)

為了實(shí)現(xiàn)分區(qū)多步計(jì)算， 按照模型的地表和速度界面情況將模型分成相對(duì)獨(dú)立的層狀區(qū)域， 這些區(qū)域之間由速度界面相連接(黃國嬌等， 2015)． 在射線追蹤時(shí)根據(jù)多次波的類型來分配速度模型和計(jì)算步驟， 通過獨(dú)立計(jì)算單個(gè)區(qū)域內(nèi)各節(jié)點(diǎn)的走時(shí)， 可以避免計(jì)算射線不需要通過的區(qū)域的走時(shí)， 節(jié)約計(jì)算時(shí)間(李曉玲等， 2013)， 這在復(fù)雜地學(xué)模型中是很有必要的．

1.1.2 模型參數(shù)化

二維模型如圖1所示， 采用不規(guī)則三角網(wǎng)格單元剖分對(duì)該模型進(jìn)行參數(shù)化． 三角單元的3個(gè)角點(diǎn)定為主節(jié)點(diǎn)， 在主節(jié)點(diǎn)之間等間距插入次級(jí)節(jié)點(diǎn)以確保計(jì)算精度． 三角單元的內(nèi)部無節(jié)點(diǎn)， 次級(jí)節(jié)點(diǎn)僅分布在網(wǎng)格單元的邊界上(圖1b)， 但炮點(diǎn)和檢波器可位于其內(nèi)(李曉玲等， 2013)．

圖1 含起伏界面的層狀模型示意圖

1.1.3 速度界面的離散化

實(shí)際問題中， 需要預(yù)先給定速度界面． 根據(jù)速度不連續(xù)界面的分布情況將界面離散化， 為保證同等計(jì)算精度， 界面上離散點(diǎn)的間距應(yīng)不低于三角網(wǎng)格單元上次級(jí)節(jié)點(diǎn)的間距， 離散點(diǎn)速度值可由單元主節(jié)點(diǎn)采樣速度通過插值的方法獲得(Baietal， 2012)．

1.1.4 波前計(jì)算和射線追蹤原理

在MTSPM算法中， 射線路徑和波前掃描是同時(shí)進(jìn)行的． 按照所需射線類型， 從震源所在區(qū)域開始逐區(qū)掃描(黃國嬌， 2014)， 即從震源所在的單元開始計(jì)算， 先計(jì)算當(dāng)前單元所有節(jié)點(diǎn)的走時(shí)， 然后進(jìn)行波前擴(kuò)展， 依次計(jì)算當(dāng)前區(qū)域所有網(wǎng)格單元的節(jié)點(diǎn)后， 使波前停留在該區(qū)的速度離散界面上， 此時(shí)保存速度界面上離散點(diǎn)的走時(shí)和相應(yīng)的射線路徑； 下一步根據(jù)所要追蹤的射線類型選擇下一個(gè)計(jì)算區(qū)域， 直至追蹤到檢波器為止(李曉玲等， 2013)． 本文中約定上標(biāo)代表上行波， 下標(biāo)代表下行波， 數(shù)字表示分區(qū)號(hào)， 例如P1P2P2P1表示第一區(qū)激發(fā)的P波下行經(jīng)界面1透射至第二區(qū)， 經(jīng)界面2反射進(jìn)入第二區(qū)后， 又經(jīng)界面1透射至第一區(qū)．

如圖1所示， 以速度不連續(xù)(間斷)面為界將模型分為3個(gè)區(qū)域． 當(dāng)計(jì)算界面1處的透射震相或反射震相時(shí)， 需要從震源所在的單元開始進(jìn)行第一區(qū)的波前擴(kuò)展計(jì)算， 第一區(qū)內(nèi)所有節(jié)點(diǎn)的走時(shí)追蹤完成后， 保存震源到界面離散點(diǎn)的最小走時(shí)和相應(yīng)的射線路徑， 此時(shí)波前停留在第一區(qū)的離散界面上． 接著， 若計(jì)算透射波P1P2， 則需要通過調(diào)用第二區(qū)的 P 波速度來進(jìn)行第二區(qū)向下的波前擴(kuò)展掃描計(jì)算， 從界面1上走時(shí)最小的點(diǎn)開始計(jì)算； 若計(jì)算反射波P1P1， 則需要通過調(diào)用第一區(qū)的P波速度進(jìn)行第一區(qū)向上的波前擴(kuò)展掃描計(jì)算． 波前擴(kuò)展計(jì)算從節(jié)點(diǎn)i到節(jié)點(diǎn)j的最小走時(shí)表達(dá)式為

(1)

式中，tij為當(dāng)前次級(jí)震源節(jié)點(diǎn)i到計(jì)算網(wǎng)格中的節(jié)點(diǎn)j的最小走時(shí)，Nj為當(dāng)前波前面上的節(jié)點(diǎn)總數(shù)，v(xi)為節(jié)點(diǎn)i的速度值，v(xj)為節(jié)點(diǎn)j的速度值，D(xi，xj)為節(jié)點(diǎn)i與j之間的距離．

多次波是由入射、 透射、 反射和轉(zhuǎn)換波按照一定的規(guī)律組合而成的， 同理可實(shí)現(xiàn)任意多震相地震射線的追蹤計(jì)算(王濤等， 2016)． 以二維起伏層狀模型(圖2)為例， 采用MTSPM算法進(jìn)行多次波射線追蹤． 該介質(zhì)模型的長度x為500 m， 深度z為600 m， 地表起伏； 模型內(nèi)有兩個(gè)起伏反射界面， 平均深度分別為250 m和500 m． 界面將模型分為3個(gè)區(qū)域， 各區(qū)域的彈性參數(shù)分別為： 第一區(qū)，vP=3.5 km/s，vS=2.02 km/s，ρ= 1.0 g/cm3； 第二區(qū)，vP=4.0 km/s，vS=2.31 km/s，ρ=1.0 g/cm3； 第三區(qū)，vP=4.5 km/s，vS=2.60 km/s，ρ=1.0 g/cm3． 炮點(diǎn)位于(250 m， 100 m)處， 34個(gè)檢波器等間距分布在地表． 圖2a表示經(jīng)由兩個(gè)界面的一次反射P波的射線路徑， 其震相為P1P1， P1P2P2P1； 圖2b表示兩個(gè)界面的一次反射S波的射線路徑， 其震相為S1S1， S1S2S2S1； 圖2c--d分別表示經(jīng)由兩個(gè)界面的一次反射轉(zhuǎn)換波的射線路徑， 其震相分別為P1S1， P1P2S2S1和S1P1， S1S2P2P1． 遵循上述原理可以計(jì)算追蹤更為復(fù)雜的多次波(多震相)射線路徑．

1.2 曲線網(wǎng)格DRP/opt MacCormack有限差分法原理

1.2.1 曲線坐標(biāo)系下的波動(dòng)方程

曲線網(wǎng)格是一種適合復(fù)雜地表介質(zhì)的網(wǎng)格離散方法， 網(wǎng)格生成的原則是使離散后的網(wǎng)格邊界與地表形態(tài)吻合， 以避免人為產(chǎn)生的階梯邊界引起的虛假散射(蔣麗麗， 孫建國， 2008)． 該方法要求界面光滑可導(dǎo)， 因此不適用于諸如斷層、 界面拐點(diǎn)、 尖點(diǎn)等不可導(dǎo)界面． 曲線網(wǎng)格可以通過計(jì)算空間到物理空間的坐標(biāo)變換來獲得， 如圖3所示．

圖2 起伏層狀均勻介質(zhì)中多震相射線路徑

圖3 曲線網(wǎng)格和計(jì)算網(wǎng)格的網(wǎng)格變換關(guān)系

曲線網(wǎng)格生成后， 建立直角坐標(biāo)系與曲線坐標(biāo)系下網(wǎng)格點(diǎn)的一一對(duì)應(yīng)關(guān)系為

(2)

(3)

由式(3)可得

(4)

其中，J=x,ξz,η-x,ηz,ξ.

直角坐標(biāo)系下， 各向同性介質(zhì)的二維彈性波一階速度-應(yīng)力方程為

(5)

其中：vx和vz為質(zhì)點(diǎn)速度；τxx，τzz和τxz為應(yīng)力分量；ρ為密度；λ和μ為拉梅常數(shù)．

利用鏈鎖法則和式(4)—(5)， 可推出曲線坐標(biāo)系下的波動(dòng)方程為

(6)

用矩陣形式可表示為

(7)

其中，

其中，χ=λ+2μ．

1.2.2 差分格式

本文采用DRP/opt MacCormack有限差分格式， Hixon和Turkel(2000)采用Tam和Webb(1993)提出的DRP (dispersion relation preserving)方法對(duì)MacCormack格式的系數(shù)進(jìn)行了優(yōu)化， 得到了低頻散、 低耗散的DRP/opt MacCormack格式． 該格式為向前差分和向后差分兩個(gè)偏心差分格式， 通過其中一個(gè)對(duì)波場進(jìn)行預(yù)估， 然后再用另一個(gè)對(duì)其矯正， 最后恢復(fù)為一個(gè)中心差分格式． 該格式通過偏心差分格式引入耗散項(xiàng)， 通過兩個(gè)偏心差分交替進(jìn)行以獲得中心差分的截?cái)嗾`差精度．

DRP/opt MacCormack格式的兩個(gè)偏心差分格式----向前差分和向后差分格式分別為

(8)

(9)

其中， 差分系數(shù)a-1=-0.30874，a0=-0.63260，a1=1.23300，a2=-0.33340，a3=0.04168.

式(8)—(9)定義了適用于一維的單邊差分算子， 對(duì)于式(6)的二維問題， Zhang和Chen(2006)采用的混合單邊差分算子分別為

(10)

在時(shí)間方向上， 本文采用Hixon(1998)， Hixon和Turkel(2000)提出的4—6階低色散和低損耗的龍格庫塔(Runge-Kutta， 簡寫為RK)法， 即使用四階RK法和六階RK法混合使用的高階時(shí)間積分， 具體實(shí)施步驟如下：

在時(shí)間循環(huán)中： 第一步采用四階RK時(shí)間積分， 表達(dá)式為

(11)

第二步采用六階RK時(shí)間積分， 表達(dá)式為

(12)

第三步再采用四階RK時(shí)間積分， 表達(dá)式為

(13)

第四步再采用六階RK時(shí)間積分， 表達(dá)式分別為

(14)

其中： Δt表示時(shí)間步長； h為中間變量； 四階RK時(shí)間積分公式中， α2=0.5， α3=0.5， α4=1， β1=1/6， β2=1/3， β3=1/3， β4=1/6； 六階RK時(shí)間積分公式中， α2=0.353323， α3=0.999597， α4=0.152188， α5=0.534216， α6=0.603907， β1=0.467621， β2=0.137286， β3=0.170975， β4=0.197572， β5=0.282263， β6=0.165142．

1.2.3 自由邊界條件

為了適應(yīng)復(fù)雜地形起伏模型中的地震波模擬， 采用曲線網(wǎng)格來離散計(jì)算區(qū)域． 在曲線坐標(biāo)系中， 本文釆用Zhang等(2012)提出的牽引力鏡像法來處理彎曲地表的自由表面條件．

考慮到自由地表處的自由表面條件要求牽引力T在地表處為零， 即

(15)

式中， σ為應(yīng)力張量， n為法向矢量， 表達(dá)式為

(16)

將式(6)代入式(15)， 得

(17)

將式(6)代入式(17)， 得

(18)

其中,

以上所述的三角網(wǎng)格單元剖分下的分區(qū)多步最短路徑射線追蹤算法(MTSPM)和曲線網(wǎng)格中求解一階速度-應(yīng)力方程的DRP/opt MacCormack有限差分算法均對(duì)起伏界面有很強(qiáng)的適應(yīng)能力， 為了驗(yàn)證這兩種方法所得結(jié)果的正確性和有效性， 下面將通過數(shù)值模擬實(shí)例進(jìn)行對(duì)比分析.

2 數(shù)值模擬實(shí)例對(duì)比研究

在數(shù)值模擬實(shí)例中， 以雙層和三層各向同性介質(zhì)模型為研究對(duì)象． 兩個(gè)模型大小均為 3.2 km×3.2 km; 震源位于介質(zhì)模型(1.6 km， 1.4 km)處， 震源采用的主頻為60 Hz， 并加載在垂直方向的瑞克(Ricker)子波上； 120個(gè)檢波器均勻地分布在1 km深度的水平面上． 在曲線網(wǎng)格DRP/opt MacCormack差分法波場數(shù)值模擬中， 網(wǎng)格間距dx=dz=1 m， 時(shí)間采樣間隔dt=1 ms， 地表為自由地表邊界條件， 人工邊界采用最佳匹配層(perfectly matched layer， 簡寫為PML)吸收邊界． 為了圖像清晰起見， 僅考慮較早時(shí)刻的直達(dá)波、 反射波和透射波的情形．

圖4 雙層介質(zhì)模型Fig.4 The two-layer medium model

2.1 雙層介質(zhì)模型

雙層介質(zhì)模型如圖4所示， 模型的界面平均深度為1.8 km， 上層介質(zhì)速度參數(shù)vP=3.5 km/s， ρ=2×103kg/m3； 下層介質(zhì)速度參數(shù)vP=4.5 km/s，ρ=3×103kg/m3． 各層vS由vS=vP/1.732計(jì)算而得.

圖5給出了雙層介質(zhì)模型中t=0.6 s時(shí)的x分量波場快照， 分別為曲線網(wǎng)格DRP/opt MacCormack差分法模擬的x分量波場快照(圖5a)和MTSPM射線追蹤理論預(yù)測的多次波波前(圖5b)． 可以看出： 隨著時(shí)間的增大， 波逐漸向外傳播， 波前面也越來越大， 在波前到達(dá)地下界面時(shí)發(fā)生了反射、 轉(zhuǎn)換和透射， 并且在地表出現(xiàn)了反射波和轉(zhuǎn)換波； 在不考慮波場傳播中能量損耗的條件下， 有限差分法計(jì)算的波前和MTSPM射線追蹤算法預(yù)測的波前具有相同的形狀和到時(shí)．

圖6給出了雙層介質(zhì)模型在上述炮檢排列下用曲線網(wǎng)格DRP/opt MacCormack差分法模擬的x分量單炮地震記錄剖面(圖6a)和MTSPM射線追蹤理論計(jì)算的相應(yīng)地震記錄剖面(圖6b， 無能量衰減)， 僅考慮初至波、 界面和自由地表的反射或反射轉(zhuǎn)換波情形． 可以看出， 在不考慮波場模擬中出現(xiàn)的能量衰減時(shí)， 兩種方法計(jì)算的地震記錄剖面有高度的相似性， 圖6b中幾乎所有出現(xiàn)的多次波震相均能從圖6a中的地震記錄中找到相應(yīng)的震相． 然而， 由于射線法是不考慮能量衰減的一種高頻近似方法， 所以一些震相， 特別是較晚時(shí)間記錄的多次波會(huì)出現(xiàn)不完全一致現(xiàn)象． 盡管如此， 但從到時(shí)和波形完全可以識(shí)別出這兩種地震記錄中相對(duì)應(yīng)的震相， 所以在地震記錄復(fù)雜的情況下， 可以參考射線追蹤法預(yù)測的多次波來識(shí)別波場數(shù)值模擬所得到的地震記錄震相．

圖5 雙層介質(zhì)模型中t=0.6 s時(shí)x分量波場快照. 虛線表示模型速度界面

圖6 雙層模型中x分量單炮地震記錄剖面

圖7給出了雙層各向同性介質(zhì)模型中曲線網(wǎng)格DRP/opt MacCormack差分法模擬合成的x分量(圖7a)和z分量(圖7b)的合成理論地震圖(檢波器位于x=1.6 km，z=1.0 km處)以及MTSPM射線追蹤算法預(yù)測可能出現(xiàn)地震震相的理論初至?xí)r間， 可以看出這兩種方法預(yù)測得到的初至?xí)r間基本吻合．

圖7 基于曲線網(wǎng)格DRP/opt MacCormack差分法得到的雙層介質(zhì)

圖8 三層介質(zhì)模型Fig.8 The three-layer medium model

2.2 三層介質(zhì)模型

三層介質(zhì)模型如圖8所示， 界面1的平均深度為1.8 km， 界面2的平均深度為2.8 km， 第一層介質(zhì)速度參數(shù)為vP=3.5 km/s，ρ=1.5×103kg/m3； 第二層介質(zhì)速度參數(shù)為vP=4.5 km/s，ρ=2.5×103kg/m3； 第三層介質(zhì)速度參數(shù)為vP=5.5 km/s，ρ=3.5×103kg/m3． 各層vS由vS=vP/1.732計(jì)算而得.

圖9給出了三層介質(zhì)模型中曲線網(wǎng)格DRP/opt MacCormack差分法模擬的x分量波場快照(圖9a)和MTSPM射線追蹤理論預(yù)測的多次波波前(圖9b)， 可以看出： 與雙層模型一樣， 隨著時(shí)間的增大， 波逐漸向外傳播， 波前面也越來越大， 在波前到達(dá)地下界面時(shí)發(fā)生了反射、 轉(zhuǎn)換和透射， 并且在地表出現(xiàn)了反射波和轉(zhuǎn)換波； 在不考慮波場傳播中能量損耗的條件下， 有限差分法計(jì)算的波前與MTSPM射線追蹤算法預(yù)測的波前具有相同的形狀和到時(shí)． 與雙層模型的波場快照相比， 三層模型除了有來自自由地表和界面1的反射外， 還有來自于界面2的反射， 因此其震相更為豐富．

圖10給出了三層介質(zhì)模型在上述炮檢排列下用曲線網(wǎng)格DRP/opt MacCormack差分法模擬得到的x分量單炮地震記錄剖面(圖10a)和MTSPM射線追蹤理論計(jì)算的相應(yīng)地震記錄剖面(圖10b， 無能量衰減)， 僅考慮了初至波、 界面和自由地表的反射或反射轉(zhuǎn)換波情形． 可以看出： 在不考慮波場模擬中出現(xiàn)的能量衰減時(shí)， 兩種方法計(jì)算得到的單炮地震記錄剖面具有高度的相似性； 從圖10a的地震記錄可見介質(zhì)模型較復(fù)雜時(shí)， 很難從地震記錄上完全分辨出不同的震相， 所以在地震記錄復(fù)雜的情況下， 可以參考射線追蹤法預(yù)測的多次波來識(shí)別波場數(shù)值模擬得到的地震記錄震相．

圖9 三層均勻介質(zhì)模型中t=0.6 s時(shí)x分量波場快照. 虛線表示模型速度界面

圖10 三層模型中x分量單炮地震記錄剖面

圖11給出了三層均勻各向同性介質(zhì)模型中曲線網(wǎng)格DRP/opt MacCormack差分法模擬合成的x分量(圖11a)和z分量(圖11b)的合成理論地震圖(檢波器位于x=1.6 km，z=1.0 km處)以及MTSPM射線追蹤算法預(yù)測可能出現(xiàn)地震震相的初至?xí)r間， 可以看出， 三層介質(zhì)模型較雙層介質(zhì)模型震相豐富， 兩種方法預(yù)測得到的初至?xí)r間基本吻合．

圖11 基于曲線網(wǎng)格DRP/opt MacCormack差分法得到的三層介質(zhì)模型中

3 討論與結(jié)論

本文分別采用曲線網(wǎng)格DRP/opt MacCormack有限差分算法和分區(qū)多步三角網(wǎng)格最短路徑射線追蹤算法對(duì)含起伏界面模型進(jìn)行了模擬計(jì)算， 通過波場快照圖、 單炮地震記錄， 以及合成理論地震圖的系統(tǒng)對(duì)比研究， 得到以下結(jié)論：

1) 在速度模型中采用不規(guī)則三角單元進(jìn)行參數(shù)化時(shí)， 可以擬合任意不規(guī)則起伏地表和地下不規(guī)則速度界面， 利用MTSPM算法可以實(shí)現(xiàn)多震相地震波走時(shí)及相應(yīng)射線路徑的追蹤計(jì)算； 而曲線網(wǎng)格有限差方法同樣可以高效地、 精確地解決地震波在包含起伏地表且各界面光滑可導(dǎo)的各向同性介質(zhì)中的波場傳播問題．

2) 由于射線追蹤僅僅追蹤了各種震相中最快的初至波， 因而無法追蹤多值射線(如逆向波、 回轉(zhuǎn)波等)， 故射線法不能全面地刻畫完全彈性波場； 同時(shí)射線法是不考慮能量衰減的一種高頻近似方法， 而波場模擬方法存在能量衰減， 因此， 在較晚時(shí)刻的地震記錄(或波場快照)中， 由于波傳播能量的不斷衰減使得震相很弱， 無法與射線法結(jié)果一一對(duì)應(yīng)．

3) 王濤等(2016)對(duì)這兩種方法在各向異性介質(zhì)中的應(yīng)用進(jìn)行的對(duì)比研究， 僅限于模型單元規(guī)則參數(shù)化的情形， 而未考慮地表起伏(或地下界面不規(guī)則)的情況． 當(dāng)?shù)乇泶嬖谄鸱?或地下界面不規(guī)則)時(shí)， 所用的兩種方法也與之前的截然不同． 本文通過對(duì)比分析驗(yàn)證了參考射線追蹤法預(yù)測的多次波， 可以簡單易行地幫助識(shí)別和標(biāo)定地震波場數(shù)值模擬得到的波場快照圖、 單炮地震記錄剖面， 以及合成理論地震圖中的主要震相． 這對(duì)于后續(xù)的地震偏移成像、 全波形反演等實(shí)際應(yīng)用具有重要的參考價(jià)值， 同時(shí)也是一種可以驗(yàn)證波動(dòng)方程數(shù)值模擬或射線法模擬正確與否的有效方法．

黃國嬌. 2014. 各向同性/各向異性介質(zhì)中多震相走時(shí)同時(shí)反演方法技術(shù)研究[D]. 西安: 長安大學(xué)地質(zhì)工程與測繪學(xué)院： 16--19.

Huang G J. 2014.SimultaneousTraveltimeInversioninComplexIsotropic/AnisotropicMediumUsingMulti-PhaseArrivals[D]. Xi’an: College of Geological Engineering and Surveying， Chang’an University： 16--19 (in Chinese).

黃國嬌， 白超英， 錢衛(wèi). 2015. 球坐標(biāo)系下多震相走時(shí)三參數(shù)同時(shí)反演成像[J]. 地球物理學(xué)報(bào)， 58(10): 3627--3638.

Huang G J， Bai C Y， Qian W. 2015. Simultaneous inversion of three model parameters with multiple phases of arrival times in spherical coordinates[J].ChineseJournalofGeophysics， 58(10): 3627--3638 (in Chinese).

蔣麗麗， 孫建國. 2008. 基于Poisson方程的曲網(wǎng)格生成技術(shù)[J]. 世界地質(zhì)， 27(3): 298--305.

Jiang L L， Sun J G. 2008. Source terms of elliptic system in grid generation[J].GlobalGeology， 27(3): 298--305 (in Chinese).

李曉玲， 白超英， 胡光義. 2013. 起伏層狀TI介質(zhì)中多次波射線追蹤[J]. 石油地球物理勘探， 48(6): 924--931.

Li X L， Bai C Y， Hu G Y. 2013. Multiple ray tracing in an undulated layered TI media[J].OilGeophysicalProspect-ing， 48(6): 924--931 (in Chinese).

王濤， 張延騰， 白超英. 2014. 2D各向異性TTI介質(zhì)中地震波場的射線法和波動(dòng)方程數(shù)值模擬對(duì)比研究[J]. 地球物理學(xué)進(jìn)展， 29(2): 580--587.

Wang T， Zhang Y T， Bai C Y. 2014. Seismic wavefield propagation on 2D anisotropic TTI media through ray theory and wave-equation simulation and comparison[J].ProgressinGeophysics， 29(2): 580--587 (in Chinese).

王濤， 白超英， 楊尚倍. 2016. 各向同性介質(zhì)中射線理論與有限差分地震波場模擬方法比較研究[J]. 地球物理學(xué)進(jìn)展， 31(2): 606--613.

Wang T， Bai C Y， Yang S B. 2016. Comparison of ray theory and FDM for simulating seismic wavefield in isotropic media[J].ProgressinGeophysics， 31(2): 606--613 (in Chinese).

王童奎， 張美根， 李小凡， 李瑞華， 龍桂華. 2007. PS轉(zhuǎn)換波界面二次源法射線追蹤[J]. 地球物理學(xué)進(jìn)展， 22(1): 165--170.

Wang T K， Zhang M G， Li X F， Li R H， Long G H. 2007. Interface points as secondary sources method raytracing PS converted waves[J].ProgressinGeophysics， 22(1): 165--170 (in Chinese).

趙愛華， 張中杰. 2004. 三維復(fù)雜介質(zhì)中轉(zhuǎn)換波走時(shí)快速計(jì)算[J]. 地球物理學(xué)報(bào)， 47(4): 702--707.

Zhao A H， Zhang Z J. 2004. Fast calculation of converted wave traveltime in 3-D complex media[J].ChineseJournalofGeophysics， 47(4): 702--707 (in Chinese).

趙愛華， 丁志峰. 2005. 寬角反射地震波走時(shí)模擬的雙重網(wǎng)格法[J]. 地球物理學(xué)報(bào)， 48(5): 1141--1147.

Zhao A H， Ding Z F. 2005. A double-grid algorithm for calculating traveltimes of wide-angle reflection waves[J].ChineseJournalofGeophysics， 48(5): 1141--1147 (in Chinese).

朱多林， 白超英. 2011. 基于波動(dòng)方程理論的地震波場數(shù)值模擬方法綜述[J]. 地球物理學(xué)進(jìn)展， 26(5): 1588--1599.

Zhu D L， Bai C Y. 2011. Review on the seismic wavefield forward modelling[J].ProgressinGeophysics， 26(5): 1588--1599 (in Chinese).

Bai C Y， Li X L， Wang Q L， Peng J B. 2012. Multiple arrival tracking within irregular triangular or tetrahedral cell model[J].JGeophysEng， 9(1): 29--38.

Carcione J M， Herman G C， Kroode A P E. 2002. Seismic modeling[J].Geophysics， 67(4): 1304--1325.

Crampin S. 1984. Effective anisotropic elastic constants for wave propagation through cracked solids[J].GeophysJInt， 76(1): 135--145.

Hixon R. 1998. Evaluation of a high-accuracy MacCormack-type scheme using benchmark problems[J].JCompAcoust， 6(3): 291--305.

Hixon R， Turkel E. 2000. Compact implicit MacCormack-type schemes with high accuracy[J].JCompPhys， 158(1): 51--70.

Igel H， Mora P， Riollet B. 1995. Anisotropic wave propagation through finite-difference grids[J].Geophysics， 60(4): 1203--1216.

Komatitsch D， Barnes C， Tromp J. 2000. Simulation of anisotropic wave propagation based upon a spectral element method[J].Geophysics， 65(4): 1251--1260.

Moser T J. 1991. Shortest path calculation of seismic rays[J].Geophysics， 56(1): 59--67.

Tam C K W， Webb J C. 1993. Dispersion-relation-preserving finite difference schemes for computational acoustics[J].JCompPhys， 107(2): 262--281.

Zhang W， Chen X F. 2006. Traction image method for irregular free surface boundaries in finite difference seismic wave simulation[J].GeophysJInt， 167(1): 337--353.

Zhang W， Zhang Z G， Chen X F. 2012. Three-dimensional elastic wave numerical modelling in the presence of surface topography by a collocated-grid finite-difference method on curvilinear grids[J].GeophysJInt， 190(1): 358--378.

Zhou H， Chen X F. 2008. Localized boundary integral equation-discrete wave number method for simulating wave propagation in irregular multiple layers， part I： Theory[J].BullSeismolSocAm， 99(3): 1984--1994.

Zhu J L， Dorman J. 2000. Two-dimensional， three-component wave propagation in a transversely isotropic medium with arbitrary orientation-finite-element modeling[J].Geophysics， 65(3): 934--942.

Comparison of seismic wavefield simulation between the curvilinear-grid finite difference method and ray method in undulated layered medium

Yang Shangbei1)Bai Chaoying1，2),*He Leiyu1)

1)CollegeofGeologyEngineeringandGeomatics，Chang’anUniversity，Xi’an710054，China2)InstituteofComputingGeophysics，Chang’anUniversity，Xi’an710054，China

Commonly it is not possible to verify and test the accuracy and validity of the seismic wavefield modeling method， especially when the model surface (or subsurface) is undulated. For this reason， this paper tries an alternative way to mutually verify and test the computational accuracy and the solution correctness of both the ray theory (the multistage triangular shortest-path method) and the wave-equation simulation method (the curvilinear-grid finite difference method of the DRP/opt MacCormack format) in layered isotropic medium. Through the analysis and comparison of wave field snapshot， common source gather profile and synthetic seismogram， it is not only able to verify the accuracy and correctness of each of the methods at least for kinematic features， but also to thoroughly understand the kinematic and dynamic features of the wave propagation in isotropic media. The results show that the curvilinear-grid finite difference method is able to yield the correct results even for complex isotropic media; the multistage triangular shortest-path method is capable of predicting similar kinematic features as the wave-equation simulation method does， which can be used to mutually test each other for methodology accuracy and solution correctness. In addition， with the aid of the ray tracing results， it is easy to identify the multi-phases (or multiples) in the wave field snapshot， common source point gather seismic section and synthetic seismogram predicted by the wave-equation simulation method， which is a key issue for later seismic application.

multistage triangular shortest-path method； curvilinear-grid finite difference method； seismic wave propagation； numerical simulation; undulated layered medium

國家科技重大專項(xiàng)子課題(2011ZX05024-001-03)資助．

2016-03-26收到初稿， 2016-05-23決定采用修改稿.

10.11939/jass.2016.06.005

P315.3+1

A

楊尚倍， 白超英， 何雷宇. 2016. 起伏層狀介質(zhì)中曲線網(wǎng)格有限差分法與射線法波場模擬對(duì)比研究. 地震學(xué)報(bào), 38(6): 854--868. doi:10.11939/jass.2016.06.005.

Yang S B, Bai C Y, He L Y. 2016. Comparison of seismic wavefield simulation between the curvilinear-grid finite difference method and ray method in undulated layered medium.ActaSeismologicaSinica, 38(6): 854--868. doi:10.11939/jass.2016.06.005.

*通訊作者 e-mail: baicy@chd.edu.cn