楊春雨　王小磊　李棟

射線追蹤方法作為一種快速有效的波場(chǎng)近似計(jì)算方法，不僅對(duì)于地震波理論研究具有重要意義，而且也直接應(yīng)用于地震波反演及偏移成像等過(guò)程。該文在收集、整理國(guó)內(nèi)外有關(guān)研究資料的基礎(chǔ)上，著重介紹了最短路徑算法，并通過(guò)對(duì)Marmousi模型的計(jì)算，驗(yàn)證了最短路徑方法的高精度，及對(duì)復(fù)雜模型的適應(yīng)性。

射線追蹤的方法種類較多。傳統(tǒng)方法有基于初值問(wèn)題的打靶法（Shooting method， i.e.， Julian and Gubbins，1977）和基于邊值問(wèn)題的彎曲法（Bending method， i.e.， Um and Thurber， 1987），但是他們不能處理介質(zhì)中較強(qiáng)的速度變化，有時(shí)無(wú)法求出全局最小走時(shí)，計(jì)算效率較低，陰影區(qū)內(nèi)無(wú)射線。隨著射線追蹤方法的發(fā)展，出現(xiàn)了直接從Huygens原理或Fermat原理出發(fā)，采用等價(jià)波前描述地震波場(chǎng)特征的方法。Vidale（1988，1990）和Podvin（1990）等人則從程函方程出發(fā)，首先求出走時(shí)場(chǎng)分布，再計(jì)算走時(shí)場(chǎng)的最速下降方向的辦法，得到每一條接收點(diǎn)到震源的射線路徑。隨后，Qin（1992）等人對(duì)Vidale的方法作了改進(jìn)，提出了波前擴(kuò)展方法。黃聯(lián)捷、李幼銘、吳如山（1992）基于Huygens 原理提出WFRT方法，根據(jù)計(jì)算精度的要求，逐次細(xì)化網(wǎng)格，而求得最小走時(shí)。Sava 和Fomel（1998）提出了HWT（huygencs wavefront tracing）法。Moser（1991）提出了根據(jù)費(fèi)馬原理的最短路徑法。

自上世紀(jì)八十年代后期發(fā)展起來(lái)的基于網(wǎng)格單元的射線追蹤算法由于其諸多優(yōu)點(diǎn)，因而倍受人們的青睞。與傳統(tǒng)的射線追蹤算法相比而言，基于網(wǎng)格單元的算法具有四大優(yōu)點(diǎn)：①可利用波振面向外擴(kuò)展傳播的原理一次性計(jì)算出速度模型中所有網(wǎng)格節(jié)點(diǎn)的射線走時(shí)及相應(yīng)的路徑，并能正確的追蹤檢波器位于射線陰影區(qū)的衍射波射線；②算法數(shù)值計(jì)算穩(wěn)定，所得到的解總是全局最佳射線路徑及相應(yīng)最小走時(shí)；③多炮多道接收時(shí)表現(xiàn)得更為高效，精度也比傳統(tǒng)射線追蹤方法要高；④在連續(xù)介質(zhì)中網(wǎng)格單元算法始終能找到初至波走時(shí)，而傳統(tǒng)射線方法則只能找出唯一走時(shí)，且很難判別這個(gè)走時(shí)是屬于初至波還是屬于后續(xù)波（De Kool et al，2006）。目前，基于網(wǎng)格單元的射線追蹤算法已成功用于地震定位、地震偏移成像和地震層析成像中。

在眾多基于網(wǎng)格單元的射線追蹤算法中，最短路徑射線追蹤算法，簡(jiǎn)稱SPM以其精度高、計(jì)算穩(wěn)健以及較強(qiáng)的模型適應(yīng)性而廣受推崇。

2 最短路徑算法

最短路徑算法是在網(wǎng)格化節(jié)點(diǎn)速度場(chǎng)中計(jì)算所有節(jié)點(diǎn)走時(shí)的一種十分有效的射線追蹤方法。這種方法是由網(wǎng)格節(jié)點(diǎn)的連線作為具有走時(shí)的射線段，采用類似Dijkstra算法來(lái)求取給定點(diǎn)到網(wǎng)格上所有節(jié)點(diǎn)的最短路徑。根據(jù)費(fèi)馬原理兩點(diǎn)間最短路徑（最小走時(shí)路徑）相應(yīng)于實(shí)際射線路徑。

最短路徑網(wǎng)格通常是由單元或中心節(jié)點(diǎn)來(lái)參數(shù)化。Nakanishi和Yamaguchi （1986）把速度場(chǎng)參數(shù)化為由常速度單元組成，而節(jié)點(diǎn)定義在單元邊界上。該算法的優(yōu)勢(shì)在于每對(duì)節(jié)點(diǎn)間的走時(shí)可以很容易的計(jì)算 ，這里d是兩節(jié)點(diǎn)間的距離，s是含節(jié)點(diǎn)單元的波慢度。計(jì)算精度可通過(guò)縮小單元尺度或增加單元邊界上節(jié)點(diǎn)的個(gè)數(shù)來(lái)提高。

另一種形成網(wǎng)格節(jié)點(diǎn)的方法是利用規(guī)則的速度節(jié)點(diǎn)，并進(jìn)行線性連接（Moser，1991）。兩連接節(jié)點(diǎn)A和B間的走時(shí)可簡(jiǎn)單地由： 得到，這里 和 分別為節(jié)點(diǎn)A和B處的波慢度值，而d則為兩節(jié)點(diǎn)間的距離。

一旦網(wǎng)格節(jié)點(diǎn)結(jié)構(gòu)及節(jié)點(diǎn)間走時(shí)計(jì)算的方式確定之后，接下來(lái)就是計(jì)算整個(gè)走時(shí)場(chǎng)和相應(yīng)射線路徑。Dijkstra （1959）最早提出原始的網(wǎng)格理論算法，該文原載于中國(guó)社會(huì)科學(xué)院文獻(xiàn)信息中心主辦的《環(huán)球市場(chǎng)信息導(dǎo)報(bào)》雜志http：//www.ems86.com總第547期2014年第15期-----轉(zhuǎn)載須注名來(lái)源其中計(jì)算時(shí)間正比于 ，而 為所有節(jié)點(diǎn)總數(shù)。此算法的概念十分簡(jiǎn)單，即：有總數(shù)為Q的未知走時(shí)節(jié)點(diǎn)，起初 Q 含M個(gè)元素，而 P 是空集，將Q集節(jié)點(diǎn)的走時(shí)設(shè)置為任意大的數(shù)。算法將炮點(diǎn)所在節(jié)點(diǎn)加入P集開(kāi)始，然后計(jì)算臨近節(jié)點(diǎn)的走時(shí)。這些組成了可能的走時(shí)，然后算法從中挑選最小走時(shí)，將其加入 P 集直至所有 Q 集內(nèi)節(jié)點(diǎn)計(jì)算完畢為止。如果 Q 集內(nèi)的節(jié)點(diǎn)在上次循環(huán)中已有計(jì)算走時(shí)，則選擇具有最小走時(shí)使之更新。完整的走時(shí)場(chǎng)可通過(guò)M次迭代得到，射線路徑可通過(guò)記錄節(jié)點(diǎn)更新的順序獲得。

3.模型實(shí)例

為了檢驗(yàn)最短路徑算法的有效性、精確性及對(duì)復(fù)雜模型的適應(yīng)性，我們選擇了較為復(fù)雜的marmousi模型。將模型劃分為三個(gè)計(jì)算區(qū)域。炮點(diǎn)位于模型左上角（0m，0m）處，910個(gè)檢波器等間距（100m）的分布在地表，最大炮檢距為9100m。我們分別使用波場(chǎng)模擬方法和最短路徑算法計(jì)算來(lái)自第一、二個(gè)不規(guī)則界面的一次反射波，并計(jì)算兩種算法計(jì)算走時(shí)的相對(duì)誤差Ere。根據(jù)計(jì)算結(jié)果，兩種算法的相對(duì)誤差在0.2%之內(nèi)，而且，最短路徑計(jì)算走時(shí)均小于波場(chǎng)模擬計(jì)算走時(shí)，更接近真實(shí)解，并且走時(shí)差大小與反射界面的起伏呈很好的相關(guān)性。說(shuō)明最短路徑算法能夠更好的模擬界面的起伏特征，保證計(jì)算的精度。

4 結(jié)論

通過(guò)使用最短路徑算法對(duì)marmousi模型的試算，我們可以得出以下結(jié)論：

最短路徑算法是一種十分有效的射線追蹤算法，具有精度高，計(jì)算穩(wěn)健的特點(diǎn)；最短路徑算法能夠適應(yīng)十分復(fù)雜的地質(zhì)模型，具有較強(qiáng)的實(shí)用性；計(jì)算過(guò)程中最短路徑算法需要較多內(nèi)存空間，計(jì)算耗時(shí)較長(zhǎng)，需要改進(jìn)。

（作者單位：1.長(zhǎng)安大學(xué)地質(zhì)工程與測(cè)繪學(xué)院；2.長(zhǎng)安大學(xué)地球科學(xué)與國(guó)土資源學(xué)院）