起伏地表地震波旅行時混合網格線性插值射線追蹤計算方法

王 琦 朱 盼 葉 佩 李 勤 李慶春*

(①長安大學地質工程與測繪學院，陜西西安 710054； ②中國國土資源航空物探遙感中心，北京 100083； ③西安科技大學地質與環境學院，陜西西安 710054)

1 引言

旅行時線性插值(Linear traveltime interpolation 簡稱LTI)是Asakawa等[1]提出的基于線性假設網格單元擴展射線追蹤算法。由于該方法計算速度快、原理簡單，得到了很大發展。趙改善等[2]結合界面二次源法將該方法推廣于追蹤反射波旅行時；聶建新等[3]將旅行時二次插值與線性插值方法聯合，降低了累積誤差；張賽民等[4]用拋物線插值取代線性插值，改善了因線性引起的誤差；張東等[5]為了求得最小走時，在向前處理過程中采用了多方向循環的計算方法。

近年來多次波越來越引起人們關注，對多次波的正確模擬和認識是壓制多次波干擾的前提。目前對多次波的射線追蹤以基于網格單元擴展的方法為主。Rawlinson等[6]利用快速步進法結合分區多步計算技術實現了層狀介質中多次透射波和多次反射波的計算，De Kool等[7]將此方法擴展到三維球面坐標系中；唐小平等[8]實現了分區多步改進的最短路徑算法，可在二維和三維層狀介質中追蹤透射、反射(或折射)和轉換的多次波；趙瑞等[9]提出了分區多步不規則最短路徑算法，也可實現多次波的追蹤計算，并與分區多步FMM算法進行對比，認為此方法比分區多步FMM算法具有高精度和高效率的優勢。本文為了說明分區多步LTI算法的優勢，與文獻[9]中的分區多步FMM算法和分區多步ISPM算法的計算精度和效率亦做了比對。

關于起伏地表地震波傳播問題，近年來國內外學者開展了相關研究。岳玉波等[10]用插值函數對起伏界面進行擬合，利用兩點射線追蹤方法實現了均勻速度模型中初值射線路徑計算，但這種方法計算效率較低； 孫章慶等[11]在矩形網格中針對起伏地表改進了快速步進法(FMM)，實現了起伏地表條件下的旅行時計算，此方法的網格仍然是單一的矩形網格，計算過程中如果網格過于稀疏會使界面呈階梯狀，引起較大誤差，加密網格又會使計算效率大打折扣且仍不能很好地逼近起伏地表或復雜構造形態； 孫章慶等[12]用FMM法中窄帶技術代替傳統的LTI法的波前擴展方式實現了起伏地表射線追蹤，該方法采用三角網與矩形網相結合的方式對模型進行剖分，使起伏地表下地震波走時的計算精度得到了提高，但并未涉及到多次波和轉換波走時的計算； Lelièvre 等[13]利用基于不規則網格的FMM法計算三維起伏地表模型的初至波旅行時，但該方法無論是在網格剖分還是在旅行時計算階段，計算量都非常大，而且此基于不規則網格的方法僅可以處理剖分三角形為銳角的情況。Sun等[14]提出在起伏地表附近采用間距不均一的網格剖分模型，用FMM法實現了起伏地表走時計算，但該方法需要對起伏地表的走時進行插值處理，這樣勢必會引起誤差； Lan等[15]通過求解曲線坐標系下的程函方程，實現了起伏地表下各向異性介質的射線追蹤，該方法為精確求解，需要對求解區域進行精細剖分，會嚴重影響計算效率； 趙后越等[16]利用規則網格和非規則網格相結合的方式進行網格剖分，結合最短路徑法實現了各向異性起伏地表條件下初至波射線追蹤，此方法為了提高計算精度會在網格邊界加入新的節點，但同時降低了運算效率。此外國內外許多學者也對基于網格單元擴展的射線追蹤算法做了有益的研究[17-26]。本文從網格剖分方式和局部走時計算兩個方面出發開展研究，提出一種能更好地逼近真實界面和更加精確高效地計算起伏地表下地震波旅行時的方法，文中討論了該方法的穩定性，用理論模型試算說明該方法具有良好的適應性。

2 矩形網格LTI方法的基本原理

圖1 矩形網格局部旅行時計算示意圖 (a)射線從矩形網格的橫邊穿過； (b)射線從矩形網格的縱邊穿過

(1)

其中tb≥ta。根據幾何關系，C點的旅行時可以表示為

(2)

(3)

可解出

(4)

進而可以求得D點的坐標

(5)

則C點的旅行時可表示為

(6)

在矩形網格中所有邊界不是垂直就是水平的，當點A、B位于橫邊時，za=zb，如圖1a所示，對應的公式可簡化為

(7)

(8)

(9)

當點A、B位于縱邊時xa=xb，如圖1b所示，對應的公式可簡化為

(10)

(11)

(12)

僅考慮以上兩種情況，即可得到矩形網格LTI方法關于D點的坐標與tc的計算公式。

LTI算法的實現過程可分為向前和向后兩個步驟進行。

第一步(向前處理)：從炮點開始，據式(6)計算炮點到網格邊界上各個節點的最小旅行時;

第二步(向后處理)：從檢波點開始，按式(5)計算所有射線路徑與各個網格邊界的交點。

以上每一步的實現過程均在文獻[1]中有詳細敘述，這里不再贅述。

3 混合網格LTI射線追蹤方法

混合網格LTI方法的實現過程可以分為三個步驟：①模型剖分：②計算局部旅行時；③計算目標類型地震波旅行時。

3.1 模型剖分

為了更接近實際地表或速度界面，同時又能兼顧計算效率，本文采用矩形網格與不規則四邊形網格相結合的方法對模型進行剖分，如圖2所示。首先對整個模型區域用矩形網格剖分，對地表和速度界面函數分別計算與網格縱邊的交點，得到離散點地表點和速度界面點； 然后將相鄰的界面離散點或網格節點選擇性連接，構成一條由多條折線段組成近似于實際地表或速度界面的折線，把這些折線段與原矩形網格邊界和節點組合，就構成了不規則四邊形網格； 最后對除地表以上的各個網格節點賦予相應速度值，就完成了對整個模型的剖分。

上述混合網格在對模型剖分時會遇到三種情況：①地表或速度界面3個相鄰的離散點都在原矩形網格的同一層內；②地表或速度界面3個相鄰離散點中左端或右端的離散點與其余兩個不在原矩形網格的同一層內；③地表或速度界面3個相鄰離散點中兩端的離散點與中間的離散點不在原矩形網格的同一層內。以上所述的三種情況分別對應于圖2中編號為①、②、③的陰影部分，為了便于觀察，將陰影部分按順序放大，如圖3所示。情況①時直接連接三個相鄰的離散點即可構成不規則四邊形網格(圖3中的①)；情況②時首先連接在原矩形網格同一層的兩個離散點，然后連接地表或界面與縱邊界無交點的原矩形網格的左下(或右下或左上或右上)角節點與中間離散點，構成非規則四邊形網格(圖3中的②)；情況③時分別連接地表或界面與縱邊界無交點的原矩形網格的左下(或右下或左上或右上)角節點、右下(或左下或左上或右上)角節點與中間離散點，構成不規則四邊形網格(圖3中的③)。圖3中各個顏色不同的區域組成了利用上述方法建立的混合網格。

圖2和圖3中藍色實線代表單一矩形網格剖分時對地表或界面的逼近，紅色實線代表在混合網格剖分時逼近的地表或界面。可以看出后者對地表或界面的擬合程度更高，剖分方法也更加合理。

圖2 混合網中模型剖分與局部旅行時計算示意圖

3.2 混合網格LTI局部旅行時計算公式的建立

通過前兩節的分析得知，適用于矩形網格的局部旅行時計算公式并不適用于混合網格，因此首先要建立針對不規則四邊形的局部旅行時計算公式。

圖4 不規則四邊形網格中局部旅行時計算示意圖

設A、B點所在直線的方程為Ax+Bz+C=0，分別將A、B點的坐標代入即可得

(za-zb)x-(xa-xb)z+za(xa-xb)-

xa(za-zb)=0

(13)

上式中tb>ta， Δt=tb-ta，求得td后，由幾何關系，C點的旅行時可以表示為

(14)

(15)

(16)

(17)

將式(16)、式(17)代入式(15)得

(18)

再將式(16)、式(17)、式(18)代入式(14)便得到了適用于不規則四邊形網格局部走時的計算公式

(19)

D點的坐標可由A、B點的坐標線性插值得到

(20)

以上即為不規則四邊形網格下求解tc和(xd,zd)的基本公式。

3.3 穩定性分析

在3.2的推導過程中，從tc關于r的一階偏導數為零，即

(21)

可以得出： Δt的符號必須與d2-r的符號一致，否則在實數域內無極值點存在。此外，D點必須介于A、B點之間或與B點(A點)重合，否則計算得到的結果無實際物理意義，或與前提假設不符。記點D與點E之間的距離為l，則

(22)

理論上l的取值存在7種情況：

(1)xd∈(xa,xe)，且zd∈(za,ze)，則l∈(0,d2)；

(2)xd∈(xe,xb)，且zd∈(ze,zb)，則l∈(-d1+d2,0)；

(3)xd=xe，且zd=ze，則l=0；

(4)xd=xa，且zd=za，則l=d2；

(5)xd=xb，且zd=zb，則l=-(d1-d2)；

(23)

或

(24)

(25)

當l∈(-∞,-d1+d2)時，求得

有

(26)

故后兩種情況無實際物理意義，在這里不予考慮。在式(19)、式(20)的基礎上針對以上前5種情況單獨考慮，可得到每種情況的表達式

(27)

(28)

(3)xd=xe，且zd=ze時，r=d2，此時Δt=0，代入式(19)、式(20)，得

tc=ta+sd3或tc=tb+sd3

(29)

(30)

(31)

(32)

(33)

(34)

以上敘述建立了不規則四邊形網格局部旅行時的計算公式，并且對其穩定性做了分析。針對每種情況分別利用相應的公式就可以穩定且高效地計算出tc和(xd,zd)。

3.4 目標波型的旅行時計算

本文提出的方法與分區多步計算技術[6]結合，就可以追蹤任意多次反射(或反射轉換)或透射(或透射轉換)波。旅行時的計算分為向前處理和向后處理兩個過程，在向前處理過程中計算每個網格節點的最小旅行時，向后處理過程是對射線路徑的追蹤以及計算整條射線路徑旅行時之和。在計算反射波或透射波的旅行時時需要引入分區多步的思想，即將整個模型區域按所需追蹤的波的類型和速度特征分成幾個獨立的計算區域，然后按照所分計算區域，在每個區域內獨立進行計算。要運用分區多步計算技術，首先要建立模型并網格化(這里采用混合網格進行剖分)，速度界面也要離散成不連續的點。

3.4.1 模型分區

圖5a所示為多次反射波在二維層狀介質模型中的分區情況，圖中展示了一條在界面折返三次的多次反射波。根據分區多步的思想，該多次反射波旅行時的計算過程可以在四個獨立的計算區域(圖中1、2、3、4代表4個分區)內進行。具體實現過程如下。該多次反射波由四段組成： 第一段是來自地表到第二個反射界面的的下行波，在計算時將第一、二層作為第一個分區來計算第一段下行波的波前； 第二段是自第二個反射界面上的反射點到第一個反射界面的上行波，在計算時將第二層作為第二個分區來計算第二段上行波的波前；第三段是自第一個反射界面反射點到第二個反射界面的下行波，在計算時將第二層作為第三個分區來計算第三段下行波的波前；第四段是自第二個反射界面反射點到地表的上行波，在計算時將第一、二層作為第四個分區來計算第四段上行波的波前。

圖5b所示為多次透射轉換波在二維層狀介質模型中的分區情況，圖中展示了一條在界面兩次透射、一次反射、三次轉換(射線每透過一次界面或每反射一次就發生一次轉換)的轉換波射線。按照分區多步思想，該多次透射轉換波旅行時的計算過程可以在四個獨立的計算區域(圖中1、2、3、4代表4個分區)內進行。假設該射線初始入射為P波，則第一段是自地面到第一個透射界面的下行P波，在計算時將第一層作為第一個分區來計算第一段下行P波的波前，此時計算所用速度為第一層的P波速度；第二段是自第一層透射點到第二層反射界面的下行S波，在計算時將第二層作為第二個分區來計算第二段下行S波的波前，此時計算所用速度為第二層的S波速度；第三段是自第二層反射點到第一層透射界面的上行P波，在計算時將第二層作為第三個分區來計算第三段上行P波的波前，此時計算所用速度為第二層的P波速度；第四段是自第一層透射點到地表的上行S波，在計算時將第一層作為第四個分區來計算第四段上行S波的波前，此時計算所用速度為第一層的S波速度。這樣就可以追蹤P-P波、P-S波、S-S波、S-P波。

圖5 二維層狀介質模型中多次反射、透射轉換波分區示意圖 (a)多次反射波； (b)多次透射轉換波

3.4.2 多步計算

圖6所示為反射(或反射轉換)或透射(或透射轉換)波的波前擴展示意圖，具體實現過程如下。

圖6 波前分區擴展示意圖

(a)波前由震源擴展到速度界面離散點； (b)若反射，由界面離散點旅行時計算上行波波前； (c)若透射，由界面離散點旅行時計算下行波波前

(1)計算由震源到第一個分區下界面的下行波的波前，并找出界面離散點上走時最小的點，如果是P波入射，則計算時用P波速度，如果是S波入射，則計算時用S波速度，如圖6a所示。

(2)將步驟(1)中計算出的走時最小的界面離散點作為新的震源。如果需要計算反射波的波前：由界面上新的震源向上計算上行波在該分區內的波前，如圖6b所示，如果需要追蹤轉換波，則采用轉換波的速度即可計算其相應的波前。如果需要計算透射波的波前：由界面上新的震源向下計算下行波在該分區內的波前，如圖6c所示，如果需要追蹤轉換波，則采用轉換波的速度即可計算其相應的波前。

(3)結合步驟(1)、步驟(2)計算波前，可以追蹤任意多次反射(或反射轉換)或透射(或透射轉換)波。

(4)由以上三個步驟計算波前，再結合LTI的向后處理就可以追蹤各類型波的射線路徑。

4 數值算例

為了驗證本文提出的計算方法的計算精度和效率，首先將分區多步LTI算法分別與分區多步FMM算法和分區多步ISPM算法的計算精度及效率進行比對，再選取三個速度模型，分別為起伏地表一維線性增加速度模型、帶有垂直斷層的層狀介質模型和加上起伏地表與速度界面的的等價Marmousi模型，進行試算。對第一個模型分別用分區多步矩形網格和混合網格LTI法計算多次反射波的旅行時，對比分析算法的計算精度。對第二個模型分別追蹤層間多次反射波和多次透射轉換波，對第三個模型追蹤其層間多次反射波，合成各個目標波形的地震記錄。

4.1 方法對比及精度分析

目前基于網格單元波前擴展的多次波波前追蹤方法中較為成熟的主要有分區多步FMM算法和分區多步ISPM算法。為了驗證分區多步LTI算法的計算精度和計算效率，對同一個速度模型分別用以上三種方法，在4種網格間距下計算反射波走時。模型尺寸為100km×40km，速度為4000m/s，地表水平，在30km處有一條水平反射界面。在模型參數化時，分區多步FMM算法中正方形網格尺寸分為四種尺度，分別為1000m×1000m、500m×500m、250m×250m和125m×125m；分區多步ISPM算法中網格尺寸為4000m×4000m，并在網格單元邊界上相應加入3、7、15和31個次級節點；分區多步LTI算法中網格大小和分區多步FMM算法中設置相同，這樣保證了以上三種方法中網格節點的間距相同。炮點在模型的左上角坐標為(0，0)，100個檢波器等間距(1000m)排列在地面，最小炮檢距為1000m。圖7和圖8分別給出了隨炮檢距變化的分區多步FMM算法(圖7a)，分區多步ISPM算法(圖7b)和分區多步LTI算法(圖8)相對于解析解在4種不同節點距下反射波走時的相對誤差。表1給出了三種方法分別所用CPU耗時。從圖7、圖8和表1中可見，除了網格間距為1000m×1000m的情形外分區多步LTI算法的精度均優于分區多步FMM和分區多步ISPM算法，在計算效率方面，分區多步LTI算法明顯優于分區多步FMM算法而低于分區多步ISPM算法。

圖7 隨著炮檢距變化的FMM算法
(a)和ISPM算法(b)數值解相對于解析解在4種不同網格距下反射波旅行時的相對誤差(引自文獻[9])

圖8 隨著炮檢距變化的分區多步LTI算法數值解相對于解析解在4種不同網格距下反射波旅行時的相對誤差

表1 均勻速度模型中FMM、ISPM和分區多步LTI算法 在4種不同網格間距下CPU耗時對比表

注：表中FMM和ISPM算法數據引自文獻[9]

圖9a是一個大小為400m×400m、每一層的速度皆為一維線性增加的三層速度模型，其速度表達式如圖9a中所列，其中h1、h2、h3表示每一層的厚度。炮點置于(200m，-20m)處，在地表上布設101個檢波點，橫向道間距均為4m，模型中剖分網格尺寸(在混合網格中僅指橫向網格間距)分別為0.1m×0.1m、0.2m×0.2m、0.4m×0.4m、1.0m×1.0m、2.0m×2.0m和4.0m×4.0m， 分別利用本文提出的方法和完全矩形網剖分的方法計算自炮點經界面三次折返到起伏地表檢波點的射線路徑的旅行時。

Kim[27]通過改變節點間距的方法來做誤差分析，此方法的前提是誤差收斂。認為隨著網格間距的縮小，誤差也會隨之減小，最終收斂于最小值；反之，隨著網格間距增大，誤差也會隨之向上浮動。

本文采用此方法進行誤差分析。為了驗證本文所采用的方法比矩形網LTI法在進行多次反射波射線追蹤時的數值穩定性以及射線路徑的可靠性更好，本文通過改變網格單元大小進行誤差分析，以最小網格單元尺寸為0.1m×0.1m的結果作為參考值，計算出其他五種網格單元尺寸相對于單元尺寸為0.1m×0.1m時各檢波點的多次反射波旅行時的相對誤差。

從圖9b(采用分區多步混合網LTI算法)、圖9c(采用分區多步矩形網LTI算法)和表2中可以看出：兩種算法計算多次反射波旅行時所得的誤差都隨網格大小縮小而逐漸收斂。但相比較而言，利用分區多步混合網格LTI算法所計算的旅行時的誤差總體上要低于分區多步矩形網格LTI算法，除網格尺寸為4.0m×4.0m以外，其余網格尺寸下計算的相對誤差基本在0.3%以內，且在相同網格尺寸下，分區多步混合網格LTI算法的計算精度高于分區多步矩形網格LTI算法，在保證精度相當的條件下，分區多步混合網格LTI算法的計算效率高于分區多步矩形網格LTI算法。據此說明分區多步混合網格LTI算法是穩定的，進行多次反射波模擬的射線路徑是可靠的。

4.2 二維層狀介質多次反射波模擬

圖10為在(200m，-20m)處單炮激發、11道接收的多次反射波波前在每個計算區域內的擴展過程和多次波的射線路徑。從圖10b～圖10e波前在各個計算區域中的擴展過程和方式可以看出，在速度均勻的同一區域內波前大致呈相互平行的弧線。當波前穿過速度分界面時，就會隨速度的改變發生變化，這是由速度分界面兩邊的速度變化引起的，這種現象符合地震波在層狀介質中的傳播規律。圖10g為根據圖10a模型模擬計算的單炮合成地震記錄，在合成記錄中分別顯示了直達波、來自Ⅰ與Ⅱ波阻抗界面的一次反射P波以及多次反射P波。從圖10g中可以看出， ②與②*相比完全不一致，如A、B兩處相差較大，這是由于檢波點分別位于水平地表和起伏地表引起的，即如果將反射波靜校正到虛擬的水平地表(浮動基準面)，地震記錄就會出現上述問題，不能客觀地反映界面的實際位置。

表2 圖9a所示模型由分區多步混合網格LTI和分區多步矩形網LTI計算所得的 多次反射波旅行時的誤差、相對誤差和CPU在5種不同網格尺寸下對比表

圖10 層間多次波波前在各分區中傳播過程及射線路徑(波前間隔為0.01s)

(a)模型參數和分區編號； (b)下行P波波前在第一個分區內傳播； (c)上行P波波前在第二個分區內傳播； (d)下行P波波前在第三個分區內傳播； (e)上行P波波前在第四個分區內傳播； (f)層間多次波射線路徑； (g)單炮合成記錄

4.3 二維層狀介質多次透射、轉換、一次反射波模擬

圖11為在(200m，-20m)處單炮激發、11道接收的多次透射轉換波波前在各個分區中的傳播過程以及多次透射轉換波的射線路徑。從圖11b～圖11e中波前在各個分區中的傳播過程可以看出，當波前以橫波速度傳播時，波前的等時線較密(如圖11c、圖11e)，當波前以縱波速度傳播時波前等時線相對較為稀疏(如圖11b、圖11d)，這是因為橫波傳播速度相對于縱波傳播速度較慢，符合地震波在介質中的傳播規律。從圖11f中可以看出：在各個速度均勻的分層中，射線路徑是直線，但是通過速度分界面時射線就會發生偏折。若入射的是S波，透射或反射為P波時，相應的出射角或反射角會增大；若入射的是P波，透射或反射為S波，相應的出射角或反射角會減小，這符合斯奈爾定律。可以說明在計算時分區徹底，計算方法正確，方法靈活性較強，可以追蹤到層狀介質中各種透射、轉換、反射波的射線路徑和旅行時。圖11g為在圖11a中模擬的單炮合成記錄。記錄中分別顯示了各道直達波、分別來自Ⅰ與Ⅱ波阻抗界面的P-P波、P-S波、圖10f中所示的多次反射P波以及圖11f中所示的多次透射轉換波。

4.4 Marmousi模型起伏檢波點反射波及多次波模擬

圖12為在人為劃定劇烈起伏地表和兩條反射界面(圖12a中加粗的黑色實線)后的Marmousi模型中，在(2500m，-500m)處單炮激發、11道接收的多次反射P波波前在各個分區中的傳播及射線路徑。由圖12b～圖12e中可以看出，波前在擴展過程中在速度較高的區域等時線被拉伸，在速度較低的區域等時線被壓縮，且波前等時線連續，滿足地震波的傳播規律，說明基于混合網格剖分的LTI方法對構造復雜的速度模型有很強的適應能力，適用于復雜速度模型的多次波模擬。圖12g給出了在圖12a的模型中模擬的單炮合成記錄。在合成記錄中可以看到初至波、分別來自于第Ⅰ與第Ⅱ波阻抗界面的一次反射P波和圖12f所示的多次反射P波。

圖11 多次透射轉換波波前在各分區中的傳播過程及射線路徑(波前間隔為0.01s) (a)模型參數和分區編號； (b)下行P波波前在第一個分區內傳播； (c)下行S波波前在第二個分區內傳播；(d)上行P波 波前在第三個分區內傳播； (e)上行S波波前在第四個分區內傳播； (f)多次透射轉換波射線路徑； (g)單炮合成記錄

(a)模型參數和分區編號； (b)下行P波波前在第一個分區內傳播； (c)上行P波波前在第二個分區內傳播； (d)下行P波波前在第三個分區內傳播； (e)上行P波波前在第四個分區內傳播； (f)層間多次波射線路徑; (g)帶地形的Marmousi模型單炮合成記錄

5 結論

本文提出了用在矩形中加入不規則四邊形網格后形成的混合網格來取代單一矩形網格對模型進行剖分，提高了對起伏地表和速度界面的處理能力，在矩形網格局部旅行時計算公式的基礎上，建立了適用于混合網格局部旅行時計算的公式，并證明了該公式的穩定性。利用分區多步混合網格LTI法實現了相關類型地震波(初至波、一次反射波、多次反射波、多次透射轉換波)旅行時的計算和射線路徑的追蹤。

通過對帶地形的層狀模型和速度不均勻模型的試算，不但說明了該方法的計算結果精度高于同等網格間距下的矩形網格LTI方法，而且證明了該方法對復雜速度模型具有良好的適應能力。為起伏地表和復雜構造情況下地震波旅行時的計算提供了一種新的方法。

[1] Asakawa E and Kawanaka T.Seismic ray tracing using linear traveltime interpolation.Geophysical Prospecting,1993,41(4):99-111．

[2] 趙改善,郝守玲.基于旅行時線性插值的地震射線追蹤算法．石油物探,1998,32(2):14-24. Zhao Gaishan,Hao Shouling.Seismic ray tracing algorithm based on the linear traveltime interpolation.GPP,1998,32(2):14-24．

[3] 聶建新,楊慧珠.地震波旅行時二次/線性聯合插值法．清華大學學報,2003,43(11):1495-1498. Nie Jianxin,Yang Huizhu.Quadrarical linear travel time interpolation of seismic ray-tracing.Journal of Tsinghua University,2003,43(11):1495-1498．

[4] 張賽民,周竹生,陳靈君等.對旅行時進行拋物型插值的地震射線追蹤方法．地球物理學進展,2007,22(1):43-48. Zhang Saimin,Zhou Zhusheng,Chen Lingjun et al.Seismic ray tracing method of applying parabolic interpolation to travel time.Progress in Geophysics,2007,22(1):43-48．

[5] 張東,謝寶蓮,楊艷等.一種改進的線性旅行時插值射線追蹤算法．地球物理學報,2009,52(1):200-205. Zhang Dong,Xie Baolian,Yang Yan et al.A ray tra-cing method based on improved linear traveltime interpolation.Chinese Journal of Geophysics,2009,52(1):200-205．

[6] Rawlinson N,Sambridge M.Multiple reflection and transmission phases in complex layered media using multistage fast marching method.Geophysics,2004,69(5):1338-1350．

[7] De Kool M,Rawlinson N,Sambridge M.A practical grid based method for tracking multiple refraction and reflection phases in three-dimensional heterogeneous media.Geophysical Journal International,2006,167(1):253-270．

[8] 唐小平,白超英.最短路徑算法下三維層狀介質中多次波追蹤．地球物理學報,2009,52(10):2635-2643. Tang Xiaoping,Bai Chaoying.Multiple ray tracing within 3-d layered media with the shortest path me-thod.Chinese Journal of Geophysics,2009,52(10):2635-2643．

[9] 趙瑞,白超英.復雜層狀模型中多次波快速追蹤——一種基于非規則網格的最短路徑算法．地震學報,2010,32(4):433-444. Zhao Rui,Bai Chaoying.Fast multiple ray tracing within complex layered media: the shortest path method based on irregular grid cells.Acta Seismologica Sinica,2010,32(4):433-444．

[10] 岳玉波.起伏地表下初值射線追蹤的實現．勘探地球物理進展,2007,30(5):388-391. Yue Yubo.Realization of initial value ray tracing for rugged topography.Progress in Exploration Geophy-sics,2007,30(5):388-391．

[11] 孫章慶,孫建國.波前快速推進法起伏地表地震波走時計算．勘探地球物理進展,2007,30(5):92-97. Sun Zhangqing,Sun Jianguo.Traveltime computation using fast marching method from rugged topography.Progress in Exploration Geophysics,2007,30(5):92-97．

[12] 孫章慶,孫建國,韓復興.復雜地表條件下基于線性插值和窄帶技術的地震波走時計算．地球物理學報,2009,52(11):2846-2853. Sun Zhangqing,Sun Jianguo and Han Fuxing.Traveltimes computation using linear interpolation and narrow band technique under complex topographical conditions.Chinese Journal of Geophysics,2009,52(11):2846-2853．

[13] Lelièvre P G,Farquharson C G and Hurich C A.Computing first-arrival seismic traveltimes on unstructured 3-D tetrahedral grids using the fast marching method.Geophysical Journal International,2011,184(2)：885-896．

[14] Sun Jianguo,Sun Zhangqing and Han Fuxing.A finite difference scheme for solving the eikonal equation including surface to topography.Geophysics,2011,76(4):53-56．

[15] Lan Haiqing and Zhang Zhongjie.Topography-dependent eikonal equation and its solver for calculating first-arrival traveltimes with an irregular surface.Geophysical Journal International,2013,193(2):1010-1026．

[16] 趙后越,張美根.起伏地表條件下各向異性地震波最短路徑射線追蹤．地球物理學報,2014,57(9):2910-2917. Zhao Houyue,Zhang Meigen.Tracing seismic shortest path rays in anisotropic medium with rolling surface.Chinese Journal of Geophysics,2014,57(9):2910-2917．

[17] 孫章慶.起伏地表條件下的地震波走時與射線路徑計算[學位論文].吉林長春:吉林大學,2011. Sun Zhangqing.The Seismic Traveltimes and Ray Path Computation under Undulating Earth’s Surface Condition[D].Jilin University,Changchun,Jilin,2011．

[18] 李永博,李慶春,吳瓊等.快速行進法射線追蹤提高旅行時計算精度和效率的改進措施．石油地球物理勘探,2016,51(3):467-473. Li Yongbo,Li Qingchun,Wu Qiong et al.Improved fast marching method for higher calculation accuracy and efficiency of traveltime.OGP,2016,51(3):467-473．

[19] 王珺,李永慶．遺傳算法和正交時頻原子相結合的地震記錄快速匹配追蹤．石油地球物理勘探,2016,51(5):881-888,893. Wang Jun,Li Yongqing.Seismic trace fast matching pursuit based on genetic algorithm and orthogonal time-frequency atom.OGP,2016,51(5):881-888,893．

[20] 秦寧,王延光,楊曉東等．非水平地表高斯束疊前深度偏移及山前帶應用實例．石油地球物理勘探,2017,52(1):81-86. Qi Ning,Wang Yanguang,Yang Xiaodong et al.Gaussian bean prestack depth migration for undula-ting-surface area in piedmont zone.OGP,2017,52(1):81-86．

[21] 黃國嬌,白超英.二維復雜層狀介質中地震多波旅行時聯合反演成像．地球物理學報,2010,53(12): 2972-2981. Huang Guojiao,Bai Chaoying.Simultaneous inversion with multiple traveltimes within 2-d complex layered media.Chinese Journal of Geophysics,2010,53(12):2972-2981．

[22] 桑運云,孫曉軍.起伏地表下基于拋物插值的最短路徑射線追蹤．石油物探,2014,53(2):1441-1448. Sang Yunyun,Sun Xiaojun.Shortest path raytracing based on parabolic traveltime interpolation in irregular topography.GPP,2014,53(2):1441-1448．

[23] 侯爵,張忠杰,蘭海強等.起伏地表下地震波傳播數值模擬方法研究進展．地球物理學進展,2014,29(2):488-497. Hou Jue,Zhang Zhongjie,Lan Haiqiang et al.Progress in numerical simulation of seismic wave propagation under an undulating surface.Progress in Geophysics,2014,29(2)：488-497．

[24] 孫章慶,孫建國,韓復興.基于線性插值和窄帶技術的走時計算方法．石油地球物理勘探,2009,44(4):436-441. Sun Zhangqing,Sun Jianguo and Han Fuxing.Travel-time computation based on linear interpolation and narrow-band technique.OGP,2009,44(4):436-441．

[25] 葉佩,李慶春.旅行時線性插值射線追蹤提高計算精度和效率的改進方法．吉林大學學報(地球科學版),2013,43(1):291-298. Ye Pei,Li Qingchun.Improvements of linear traveltime interpolation ray tracing for the accuracy and efficiency.Journal of Jilin University (Earth Science Edition),2013,43(1):291-298．

[26] 宋御杰.基于混合網格剖分的起伏地表旅行時計算方法研究[學位論文].陜西西安: 長安大學,2015. Song Yujie.The Travel-time Computation Method Study Under Undulating Earth’s Surface Based on Hybrid mesh[D].Chang’an University,Xi’an,Shaanxi,2015．

[27] Kim S.3-D eikonal solvers: First-arrival traveltimes.Geophysics,2002,67(4):1225-1231.