李同宇 張建中

(①海底科學(xué)與探測技術(shù)教育部重點(diǎn)實(shí)驗室，山東青島 266100； ②青島海洋科學(xué)與技術(shù)試點(diǎn)國家實(shí)驗室海洋礦產(chǎn)資源評價與探測技術(shù)功能實(shí)驗室，山東青島 266061； ③中國海洋大學(xué)海洋地球科學(xué)學(xué)院，山東青島 266100)

1 引言

地震射線追蹤及走時計算是地震學(xué)的基本問題之一，被應(yīng)用于層析成像、地震定位、偏移以及地震數(shù)據(jù)采集設(shè)計等領(lǐng)域[1，2]。傳統(tǒng)的射線追蹤方法主要針對兩點(diǎn)射線追蹤問題，可分為試射法[3]和彎曲法[4]。這兩種方法都存在一定局限性： 前者需不斷試驗射線入射角，計算效率較低； 后者易于陷入局部最優(yōu)解。20世紀(jì)80年代后期，出現(xiàn)了基于地震波前走時的射線追蹤方法，如程函方程的有限差分法[5-8]、最短路徑法[9-13]和線性插值算法[14-24]等。與傳統(tǒng)方法相比，該類方法具有較高計算精度和效率，得到廣泛研究和應(yīng)用。

Asakawa等[14]提出基于線性插值的二維走時計算和射線追蹤方法(LTI)，指出由于可通過對LTI求導(dǎo)得到程函方程有限差分[5]公式，因而LTI方法可被視作程函方程有限差分法的更高版本，并對比說明了LTI方法具有更高的計算精度和效率。但是，LTI算法從震源開始逐步向外圍推進(jìn)的過程中，考慮的波傳播方向有限，計算的節(jié)點(diǎn)處走時不一定是最小走時，也不能正確追蹤逆向傳播的射線。為了克服這些問題，該方法被不斷的修正和改進(jìn)。張建中等[15，16]結(jié)合LTI與最短路徑算法，提出了動態(tài)網(wǎng)絡(luò)射線追蹤方法。黃月琴等[17]將二維LTI算法推廣至三維情形。Zhang等[18]結(jié)合雙線性走時插值和波前擴(kuò)展法，提出了基于規(guī)則單元的三維射線追蹤方法，隨后Huang等[19]將該方法推廣至適用于不規(guī)則六面體單元。李培明等[20]提出了一種改進(jìn)的雙線性插值射線追蹤方法，并結(jié)合動態(tài)變密度網(wǎng)格剖分，提高了旅行時的計算精度和效率。黃翼堅等[21]通過對傳統(tǒng)初至波LTI算法的射線方向進(jìn)行約束，得到了一種直達(dá)波旅行時LTI迭代算法。張東等[22]通過使用LTI計算向前過程的離散旅行時場，并用B樣條插值梯度反向追蹤射線路徑，提出了初至波走時梯度射線追蹤算法(MTG)，較大程度消除了計算網(wǎng)格離散化帶來的誤差。在此基礎(chǔ)上，張婷婷等[23]結(jié)合逐次網(wǎng)格剖分技術(shù)及MTG方法，實(shí)現(xiàn)了三維層狀介質(zhì)中的多波走時計算及射線追蹤，隨后針對在速度突變區(qū)域用B樣條插值計算旅行時梯度會帶來誤差的問題，張婷婷等[24]在三維MTG射線追蹤算法基礎(chǔ)上提出了一種改進(jìn)的B樣條/線性聯(lián)合插值的三維射線追蹤(MMTG)算法。

LTI方法假設(shè)地震走時沿單元邊界線性變化，單元邊界上任意一點(diǎn)的走時可由該邊界上已知點(diǎn)走時的線性插值函數(shù)表示。但實(shí)際上，單元邊界上的走時呈非線性變化，因而當(dāng)離散單元較大時，LTI方法的線性假設(shè)會導(dǎo)致較大的走時計算和射線追蹤誤差。針對這一問題，Zhang等[25]提出了基于規(guī)則單元的線性走時擾動插值(LTPI)方法，提高了波前走時的計算精度。為了更好模擬地形和地層界面的起伏及各地層速度的變化，本文將復(fù)雜介質(zhì)模型離散成不規(guī)則單元，推導(dǎo)了二維情況下適用于不規(guī)則單元的LTPI方程，并結(jié)合波前擴(kuò)展算法，提出了一種基于線性走時擾動插值的射線追蹤方法，數(shù)值實(shí)驗證明該方法具有較高的計算精度、效率以及很強(qiáng)的對復(fù)雜介質(zhì)的適應(yīng)性。

2 方法原理

基于LTPI的射線追蹤算法包括兩步：波前走時計算和射線追蹤。在波前走時計算過程中，根據(jù)惠更斯原理，結(jié)合群波前擴(kuò)展算法(GMM)和LTPI方法，從震源開始逐步向外計算離散模型中各網(wǎng)格節(jié)點(diǎn)的波前走時；追蹤射線時，根據(jù)費(fèi)馬原理，從接收點(diǎn)開始反向逐步確定射線與相關(guān)單元邊界的交點(diǎn)，直至追蹤至炮點(diǎn)，將炮點(diǎn)、各個交點(diǎn)和接收點(diǎn)順次連接，則得到從炮點(diǎn)到接收點(diǎn)的射線路徑。

2.1 模型離散

為模擬地形與地層界面的起伏，復(fù)雜模型被離散成如圖1所示的一系列不規(guī)則四邊形單元[19]。具體離散過程如下。

首先，在x方向等間距剖分模型。記模型在x方向的起始坐標(biāo)為xc，剖分間距為Δx，剖分網(wǎng)格數(shù)目為n。離散模型的橫坐標(biāo)可表示為

x(i)=xc+iΔxi=0，1,2,…，n

(1)

然后，由淺至深分別將每個地層劃分成多個薄層。在同一地層中，薄層數(shù)目相等，并且在同一x坐標(biāo)處，各薄層在z方向上的厚度相等。不同地層的薄層數(shù)目可以不同。假設(shè)第k個地層中劃分的薄層數(shù)目為dk，則第k個地層中的第m個薄層在x=x(i)處的底面深度值可表示為

(2)

式中l(wèi)k(i)和lk+1(i)分別表示第k個地層的頂面和底面在x=x(i)處的深度值。

這種離散方式可很好地擬合起伏地形和地層界面，不僅考慮了每個地層內(nèi)速度的縱向和橫向變化以及地層速度與界面的耦合問題，且易于自動實(shí)現(xiàn)。

2.2 基于走時擾動插值的走時計算

在波前走時計算過程中，需計算與各級震源點(diǎn)相鄰的網(wǎng)格節(jié)點(diǎn)處的波傳播時間，即在一個單元內(nèi)由已知節(jié)點(diǎn)走時計算未知節(jié)點(diǎn)走時。在單元較大時，LTI方法的線性走時假設(shè)會導(dǎo)致較大的計算誤差，因此本文引入線性走時擾動插值(LTPI)方法[25]，并將它擴(kuò)展到適于不規(guī)則單元的情況。

在圖2中，假設(shè)射線穿過不規(guī)則單元的上頂邊界。S點(diǎn)表示震源點(diǎn)，P1和P2為單元上頂邊界兩個走時已知節(jié)點(diǎn)，坐標(biāo)分別為(xi，zi)(i=1，2)，其波前走時分別為t1和t2。需計算射線從震源點(diǎn)S穿過單元上頂邊界傳至該單元邊界節(jié)點(diǎn)E的走時。

設(shè)從震源傳至E點(diǎn)的射線與單元上頂邊界的交點(diǎn)為P點(diǎn)，設(shè)上頂邊界的中點(diǎn)為P0(x0，z0)，P點(diǎn)坐標(biāo)可表示為

(3)

式中：r為P點(diǎn)與P0點(diǎn)橫坐標(biāo)之差，當(dāng)P點(diǎn)在P0點(diǎn)右側(cè)時r大于零，否則r小于零；a0為常數(shù)，有

(4)

若P點(diǎn)走時為tP，則E點(diǎn)走時tE可表示為

tE=tP+slPE

(5)

式中：lPE表示P點(diǎn)至E點(diǎn)的射線長度；s表示該單元內(nèi)的平均慢度。

現(xiàn)在采用LTPI方法[25]計算P點(diǎn)走時和坐標(biāo)。將圖2中震源點(diǎn)S至單元上頂邊界之間的介質(zhì)等效成勻速介質(zhì)。在等效勻速介質(zhì)中，從S點(diǎn)至單元上頂邊界的射線路徑為直線。若S點(diǎn)至P1、P2點(diǎn)的直線長度分別為l1、l2，則等效勻速介質(zhì)的平均慢度seq可表示為

(6)

由此可求得等效勻速介質(zhì)中從震源點(diǎn)S至點(diǎn)Pi(i=1，2)的射線走時為

(7)

Δti=ti-seqlii=1,2

(8)

對于均勻模型，等效介質(zhì)射線與實(shí)際射線重合，走時擾動為零。對于非均勻模型，等效勻速介質(zhì)中的直射線與實(shí)際射線不同，走時擾動一般為非零值，且遠(yuǎn)小于參考走時。

據(jù)式(6)～式(8)，P點(diǎn)處走時擾動可表示為

ΔtP=tP-seqlSP

(9)

式中l(wèi)SP為等效勻速介質(zhì)中S點(diǎn)至P點(diǎn)的射線長度。將式(9)代入式(5)，S點(diǎn)至E點(diǎn)的走時可表示為

tE=ΔtP+seqlSP+slPE

(10)

式中P點(diǎn)處的波前走時被分解為參考走時seqlSP和走時擾動ΔtP。其中參考走時由定義求得，而ΔtP可由單元邊界上已知點(diǎn)的走時擾動表示。

假定走時擾動沿單元邊界線性變化，P點(diǎn)的走時擾動ΔtP可用線性方程表示為

ΔtP=ar+b

(11)

式中a和b為常系數(shù)，可由P1和P2點(diǎn)處的已知走時擾動Δt1和Δt2表示

(12)

設(shè)炮點(diǎn)S坐標(biāo)為(xS，zS)，E點(diǎn)坐標(biāo)為(xE，zE)，據(jù)式(3)有

(13)

為簡化計算， 將lSP=f1(r)和lPE=f2(r)在r=0處分別進(jìn)行泰勒級數(shù)展開，保留至二次方項，并把f1(0)和f2(0)分別簡化記為f1和f2，則有

(14)

其中

將式(6)、式(11)和式(14)代入式(10)中，由費(fèi)馬原理可知，E點(diǎn)的波前走時滿足公式

(15)

求解式(15)可得

(16)

其中

由式(16)求得r值后，根據(jù)式(3)可知射線與單元上頂邊界的交點(diǎn)P的坐標(biāo)，將r代入式(10)即可求得射線經(jīng)過單元上頂邊界到達(dá)E點(diǎn)的走時。射線穿過單元其他邊界時的走時計算方法與此類似。

若不規(guī)則單元的其他邊界節(jié)點(diǎn)的走時也已知，同樣也可用上述方法求出射線穿過其他邊界到達(dá)E點(diǎn)的走時，然后將各個走時中的最小者作為所求E點(diǎn)處的波前走時。

2.3 波前擴(kuò)展

本文采用群波前擴(kuò)展算法(GMM)擴(kuò)展波前。復(fù)雜介質(zhì)離散后，每一個單元都被認(rèn)為是均勻的，射線在單元內(nèi)沿直線傳播。波前擴(kuò)展的二維情形如圖1 所示。從震源點(diǎn)開始，利用LTPI計算出與其相鄰的各個網(wǎng)格節(jié)點(diǎn)上的波傳播時間，這些節(jié)點(diǎn)組成波前窄帶。在波前窄帶中需要按照一定規(guī)則找出次級震源點(diǎn)，再從次級震源點(diǎn)處開始計算與其相鄰的網(wǎng)格節(jié)點(diǎn)處的波傳播時間，這些網(wǎng)格節(jié)點(diǎn)便組成了新的波前窄帶[19]。波前窄帶需要不斷向外推進(jìn)直至獲得介質(zhì)中全部網(wǎng)格節(jié)點(diǎn)處的波前走時。

記當(dāng)前波前窄帶為N，令

(17)

式中：sN，min表示波前窄帶各個網(wǎng)格節(jié)點(diǎn)中波傳播的最小慢度，即最大波傳播速度； Δx和Δz為波前窄帶各網(wǎng)格在x和z方向上的離散網(wǎng)格距。

由GMM算法可知若波前窄帶中兩個相鄰網(wǎng)格節(jié)點(diǎn)的波傳播時間之差小于δtN，則這兩點(diǎn)處的波傳播能量不會互相影響。次級震源集合G的選擇需遵循如下規(guī)則

G={(i,j)∈N,ti,j≤tN,min+δtN}

(18)

式中tN,min表示波前窄帶N中的最小波傳播時間。

波前走時計算步驟如下：

(1)將震源點(diǎn)的波前走時記為0，并且記M=2，表示該處已獲得波前走時； 其余網(wǎng)格節(jié)點(diǎn)的波前時間均記為無窮大，并令其M=0，表示該網(wǎng)格節(jié)點(diǎn)未計算波前走時；

(2)計算與震源相鄰的網(wǎng)格節(jié)點(diǎn)處的波的傳播時間后，將時間保存到波前窄帶N中，并令M=1，表示該處網(wǎng)格節(jié)點(diǎn)已計算出波傳播時間，但還未確定波前時間；

(3)波前窄帶中滿足式(18)的網(wǎng)格節(jié)點(diǎn)處的波傳播能量不會相互干擾，因此將其作為次級震源，將滿足條件的網(wǎng)格節(jié)點(diǎn)從波前窄帶N中移至次級震源集合G中，并令M=2；

(4)對次級震源集合G中的網(wǎng)格節(jié)點(diǎn)，更新與其相鄰的網(wǎng)格節(jié)點(diǎn)處波傳播時間，要注意相鄰網(wǎng)格節(jié)點(diǎn)需滿足M≠2，若節(jié)點(diǎn)處M=0，則將該節(jié)點(diǎn)移至波前窄帶N中并令M=1；

(5)若集合N非空，則跳轉(zhuǎn)至步驟(3)；若集合N為空，則波前擴(kuò)展結(jié)束。

2.4 射線追蹤

在計算出復(fù)雜介質(zhì)全部網(wǎng)格節(jié)點(diǎn)處的波前走時后，根據(jù)費(fèi)馬原理和互換原理，可從接收點(diǎn)出發(fā)反向逐步確定出射線與相應(yīng)單元邊界的交點(diǎn)，直至追蹤至炮點(diǎn)； 然后順次連接該對炮點(diǎn)和接收點(diǎn)之間確定的所有交點(diǎn)，便得到具有最小走時的直達(dá)波射線路徑。具體射線追蹤步驟如下。

(1)在接收點(diǎn)R所在單元內(nèi)，利用LTPI方法計算射線路徑與該單元各邊界的交點(diǎn)P(圖3)。若接收點(diǎn)位于單元節(jié)點(diǎn)或單元邊界上，則需要在接收點(diǎn)所在的所有單元中進(jìn)行計算。

(2)將具有最小走時的交點(diǎn)P作為新的接收點(diǎn)，判斷接收點(diǎn)是否位于炮點(diǎn)所在單元，若是，執(zhí)行步驟(3)； 否則返回步驟(1)。

(3)順次連接炮點(diǎn)、射線與相應(yīng)單元邊界的各個交點(diǎn)和接收點(diǎn)，就得到了該對炮點(diǎn)和接收點(diǎn)對應(yīng)的直達(dá)波射線路徑。

圖3 一個單元內(nèi)射線路徑追蹤示意圖

3 數(shù)值實(shí)驗

為測試本文方法，設(shè)計了四種模型：均勻介質(zhì)速度模型(模型1)、垂向漸變速度模型(模型2)、傾斜層狀速度模型(模型3)和含異常體速度模型(模型4)。前兩種模型的波前走時和射線路徑有解析解，可使用均方根誤差

(19)

3.1 均勻介質(zhì)速度模型

模型1是尺寸為2000m×2000m的均勻介質(zhì)速度模型，離散為一系列80m×80m的規(guī)則單元。模型速度為4000m/s，震源坐標(biāo)為(0，1000m)。圖4a和圖4b分別為用LTI和LTPI方法計算得到的波前走時(紅色虛線)與理論值(灰色實(shí)線)的對比圖，用相應(yīng)方法計算的波前走時與理論值之間的絕對誤差分布分別如圖4c和圖4d所示?？梢奓TPI方法計算的波前走時更接近理論值，絕對誤差更小。

圖4 模型1波前走時等值線(單位：ms)及絕對誤差(a)LTI方法走時； (b)LTPI方法走時； (c)LTI方法誤差； (d)LTPI方法誤差

為說明計算誤差和計算時間隨離散單元大小的變化，對模型1分別以10、20、40、50、80、100和200m的網(wǎng)格尺寸進(jìn)行離散，然后計算了LTI和LTPI方法的波前走時均方根誤差，并記錄了CPU時間(圖5)?？梢娫谙嗤x散網(wǎng)格距下，LTI和LTPI方法所用CPU時間相近，而LTPI方法的均方根誤差值要遠(yuǎn)小于LTI方法。隨著網(wǎng)格距的增大，LTI方法的均方根誤差值增長迅速，而LTPI方法的均方根誤差值則變化較為平緩。因此在一定誤差要求下，LTPI方法可以使用比LTI方法更少的單元數(shù)目，從而提高計算效率。

將模型1以80m的網(wǎng)格尺寸進(jìn)行離散，分別使用LTI和LTPI方法計算射線路徑，與理論路徑的對比如圖6a和圖6b所示。使用相應(yīng)方法計算的各個接收點(diǎn)對應(yīng)的射線走時及其相對誤差對比分別如圖6c和圖6d所示。相比LTI方法，LTPI方法計算的射線路徑和射線走時更接近理論值。在不同接收點(diǎn)處使用LTPI方法得到的走時相對誤差穩(wěn)定在1.0×10-3附近，小于LTI方法對應(yīng)的誤差值?？梢姡琇TPI方法具有更高的計算精度。

圖5 LTI和LTPI方法在不同離散網(wǎng)格距下的均方根誤差和CPU時間

3.2 垂向漸變速度模型

速度隨深度線性變化介質(zhì)的波前走時和射線路徑具有解析解[26]。為將本文方法結(jié)果與變速介質(zhì)的理論值進(jìn)行對比，建立了速度隨深度線性增加的模型2，其尺寸為2000m×2000m，離散網(wǎng)格尺寸為40m×40m，速度隨深度線性變化函數(shù)為v=1500+2z(深度單位為m，速度單位為m/s)。炮點(diǎn)位于(1000m，0)。圖7a和圖7b分別為用LTI和LTPI方法計算出的波前走時(紅虛線)與理論值(灰實(shí)線)的對比圖，圖7c和圖7d分別是用相應(yīng)方法計算的波前走時絕對誤差分布圖?？梢奓TPI方法比LTI方法計算的波前走時更貼近理論值，計算精度更高。

圖6 模型1的射線路徑及走時對比圖(a)LTI方法射線； (b)LTPI方法射線； (c)計算走時與理論走時； (d)計算走時相對誤差

圖7 模型2波前走時等值線(單位：ms)和絕對誤差分布圖(a)LTI方法走時； (b)LTPI方法走時； (c)LTI方法誤差； (d)LTPI方法誤差

為說明本文方法對離散單元大小的響應(yīng)，將模型2分別以10、20、40、50、80、100和200m網(wǎng)格尺寸做離散，分別用LTI和LTPI方法計算波前走時均方根誤差并記錄了CPU時間(圖8)?？梢娫谕痪W(wǎng)格距下，LTI與LTPI方法所用的CPU時間相近，但LTI方法的均方根誤差比LTPI方法大一倍以上。隨著網(wǎng)格距的增大，LTI方法的均方根誤差值快速增大，而LTPI方法的均方根誤差值變化相對平緩。因此在一定誤差要求下，LTPI方法可使用比LTI方法更少的單元數(shù)目以提高計算效率。

將模型2以40m的網(wǎng)格尺寸離散，炮點(diǎn)坐標(biāo)為(1000m，2000m)，圖9a和圖9b分別為使用LTI和LTPI方法計算的射線路徑與理論路徑的對比圖，不同接收點(diǎn)處分別使用這兩種方法得到的射線走時及其相對誤差分別如圖9c和圖9d所示。可見在速度隨深度線性增加的模型中，使用LTPI方法得到的射線路徑和射線走時更加貼近理論值，并且對應(yīng)的走時相對誤差要小于LTI方法。

圖8 不同網(wǎng)格距下LTI和LTPI方法計算的均方根誤差和CPU時間

圖9 模型2中射線路徑及走時對比圖(a)LTI方法射線； (b)LTPI方法射線； (c)計算走時與理論走時； (d)計算走時相對誤差

3.3 傾斜層狀速度模型

模型3由三層傾斜層狀介質(zhì)組成，各層速度從上到下分別為2000、3000和4000m/s。模型尺寸為2000m×2000m，震源位于(1000m，0)。分別使用LTI和LTPI方法計算離散網(wǎng)格距為20m(灰色實(shí)線)和80m(紅色虛線)時的波前走時(圖10)?？梢婋x散網(wǎng)格尺寸變化時，使用LTPI方法得到的波前走時更加一致，與LTI方法相比，LTPI方法具有良好的穩(wěn)定性，能更好地適應(yīng)大網(wǎng)格距離散模型。

用網(wǎng)格距為40m的網(wǎng)格離散模型3。為獲取傾斜層狀介質(zhì)中的理論射線路徑，先設(shè)定炮點(diǎn)位置為(1000m，1800m)，再從炮點(diǎn)處以一定的出射角出射射線，根據(jù)Snell定律確定射線路徑，將射線在地表的出射點(diǎn)作為接收點(diǎn)。然后分別使用LTI和LTPI方法計算了這些炮檢點(diǎn)所對應(yīng)的射線路徑(圖11a和圖11b所示)?？梢娪肔TPI方法計算的射線路徑更接近理論射線路徑。不同接收點(diǎn)處對應(yīng)的射線走時及其相對誤差分別如圖11c和圖11d所示，使用LTPI方法計算得到的走時相對誤差穩(wěn)定在1.0×10-3以下，小于LTI方法的相對誤差。實(shí)驗結(jié)果說明LTPI方法比LTI方法具有更高的計算精度。

圖10 離散網(wǎng)格距分別為20m和80m時，模型3波前走時等值線圖(單位ms)(a)LTI方法； (b)LTPI方法

圖11 模型3中射線路徑與走時對比圖(a)LTI方法射線； (b)LTPI方法射線； (c)計算走時與理論走時； (d)計算走時相對誤差

3.4 含異常體速度模型

模型4尺寸為2000m×2000m，由三層起伏層狀介質(zhì)組成，從上到下速度分別為1700、2800和1800m/s。在上層介質(zhì)中存在速度為2100m/s的高速體，中間層介質(zhì)含一速度為2300m/s的低速體。震源位于(1000m，2000m)。分別用LTI和LTPI方法計算了離散網(wǎng)格距為20m(灰色實(shí)線)和80m(紅色虛線)時波前走時(圖12)?？梢妰煞N方法計算的波前走時趨勢基本一致，反映出遇到低速體、高速體和地層界面時的波前變化特征。相比LTI方法，LTPI方法在網(wǎng)格距增大時走時變化更小，穩(wěn)定性更高。

圖12 離散網(wǎng)格距分別為20m和80m時模型4波前走時等值線(單位ms)圖(a)LTI方法波前走時； (b)LTPI方法波前走時

將模型4分別用20m和80m的網(wǎng)格距離散，并分別使用LTI和LTPI方法計算射線路徑(圖13)?？梢奓TI和LTPI兩種方法計算的射線路徑都表現(xiàn)出盡量避開低速體和趨向高速體的特征。但是，當(dāng)模型離散網(wǎng)格距發(fā)生變化時，LTI方法計算的射線路徑變化較大，而LTPI方法計算得到的射線路徑更為合理且變化更小，反映出LTPI方法更具優(yōu)勢。

圖13 模型4射線路徑圖(a)LTI方法； (b)LTPI方法

4 結(jié)論

線性走時擾動插值(LTPI)方法將實(shí)際波前走時分解為等效勻速介質(zhì)中的參考走時與走時擾動。參考走時沿單元邊界非線性變化，雖然假設(shè)走時擾動在單元邊界上線性變化，但與參考走時相比，走時擾動很小。因此，參考走時與走時擾動之和仍保持著沿單元邊界非線性變化的特征，從而克服了常用的線性走時插值(LTI)方法在離散單元較大時計算走時誤差大的問題，提高了波前走時的計算精度。

采用不規(guī)則單元離散復(fù)雜介質(zhì)模型，能夠更好地模擬地形和地層界面的起伏形態(tài)以及各地層內(nèi)速度的縱橫向變化。針對不規(guī)則單元，建立了基于LTPI的走時計算方程，并結(jié)合GMM波前擴(kuò)展算法，提出了一種基于線性走時擾動插值的射線追蹤方法，對復(fù)雜介質(zhì)模型具有很強(qiáng)的適應(yīng)性。

不同模型的測試結(jié)果表明， 相比LTI方法，LTPI方法在波前走時和射線路徑方面具有更高的計算精度，并且隨著離散模型單元尺寸的增加，LTPI方法的計算誤差增加更為緩慢。在離散單元數(shù)目相同的情況下，LTPI與LTI方法的計算效率相當(dāng)，但LTPI的計算精度比LTI方法高出一倍以上；在一定精度要求下，LTPI方法可使用離散單元數(shù)目更少的模型，從而提高計算效率。