楊春雨　王小磊　李棟

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

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

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

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

2 最短路徑算法

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

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

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

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

3.模型實例

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

4 結(jié)論

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

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

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