孫章慶，孫建國＊，岳玉波，江兆南

1吉林大學地球探測科學與技術學院，長春 130026

2國土資源部應用地球物理綜合解釋理論開放實驗室—波動理論與成像技術實驗室，長春 130026

3東方地球物理公司物探技術研究中心，河北涿州 072751

4重慶市地質(zhì)礦產(chǎn)勘查開發(fā)局208水文地質(zhì)工程地質(zhì)隊，重慶 408300

1 引言

地震波走時是描述地震波運動學特征的一個很重要的物理參數(shù)，其刻畫了地震波由震源點激發(fā)后到達地下介質(zhì)空間中各個位置坐標點所需耗消的走時，而該走時和地下介質(zhì)的速度參數(shù)的分布情況密切相關，所以地震波走時計算方法常常被應用于走時反演、層析成像、偏移成像、速度分析等一些重要的地震數(shù)據(jù)處理技術中.地震波走時計算的精度、效率、穩(wěn)定性等性質(zhì)直接影響著這些地震數(shù)據(jù)處理的效果和效率（Sun，1998，1999，2000，2004；井西利等，2007；瞿辰等，2007），同時考慮到三維地震勘探正在全面大范圍展開的大背景，所以對三維地震波走時計算進行研究具有很重要的理論意義和實際價值.

在走時計算方面，目前采用的主要方法有兩點射線追蹤法（Julian and Gubbins，1977）、最短路徑射線追蹤法（Moser，1991；張美根等，2006）、波前構建法（Vinje et al.，1993）、走時插值法（Asakawa and Kawanaka，1993；Wang and Ma，1999）、有限差分法（Vidale，1988；Sethian and Popovici，1999；Sun et al.，2011）、辛幾何算法（秦孟兆和陳景波，2000；高亮等，2000）等.其中，兩點射線追蹤法在早期天然地震數(shù)據(jù)處理中起著非常重要的作用，但在現(xiàn)代三維地震數(shù)據(jù)處理中其計算效率相對較低.最短路徑射線追蹤法具有不錯的走時計算精度，但由于其需要設置很多網(wǎng)格線上的計算節(jié)點，所以在求解三維問題時其需要耗費相對更大的計算和存儲成本.波前構建法能同時計算走時和射線路徑，并且其計算精度也相對不錯，不過在三維算法時，其中的網(wǎng)格定位、插值問題等都相對非常復雜.實際上，走時插值法和有限差分法是目前被廣泛采用的三維走時算法，其中有限差分法具有計算效率高、局部走時算法簡潔等優(yōu)勢，但其計算精度有限，而具有更高計算精度的高階有限差分需要花費相對更多的計算成本.走時插值法具有不錯的計算精度，但是其算法的整體實現(xiàn)策略的效率相對有限差分法低很多.因此，本文將采用走時插值法作為局部走時算法，而在整體實現(xiàn)策略上采用有限差分算法中的快速算法.

目前，三維走時插值法主要采用B樣條插值法（張東等，2013）、雙線性插值法（李培明等，2013；劉鋒等，2012；梅勝全等，2010；張東等，2009）等算法.在雙線插值法中，有兩方面關鍵算法.一方面是局部走時算法，其包括計算公式和局部計算策略等技術環(huán)節(jié)；另一方面是算法的整體實現(xiàn)策略，也即算法的整體實現(xiàn)步驟.在局部走時算法方面，目前普遍認為雙線性插值法求解極小值的方程為一個超越方程，無法求出解析解（其可能解分布在一條雙曲線上），所以一些學者分別采用網(wǎng)格界面剖分法（張東等，2009）、快速插值法（李培明等，2013；梅勝全等，2010）、最速下降法（劉鋒等，2012）來獲得一個約束條件下接近于真實解的一個近似解.而在走時計算的整體實現(xiàn)步驟方面，目前普遍采用“向前處理”的按行列掃描計算（張東等，2009）、最短路徑射線追蹤算法的實現(xiàn)步驟（李培明等，2013；梅勝全等，2010）等.針對以上兩個方面的關鍵技術，筆者充分利用平面波雙線性假設的結論，推導出了雙線性法求解極小值超越方程的解析計算公式，同時構造性地證明了該解析計算公式的正確性和相對應的地震波傳播的物理規(guī)律.此外，筆者還借鑒快速推進法中的迎風差分格式的構建思想，提出了迎風雙線性插值的局部走時計算的實現(xiàn)策略，并以快速推進法中的窄帶技術（這是一種滿足地震波傳播基本規(guī)律的、計算效率高且靈活的波前擴展算法）作為算法的整體實現(xiàn)策略.最后的精度分析和計算實例表明了本文算法的正確性和在面對復雜介質(zhì)時的穩(wěn)定性和有效性.

2 算法的描述

為了獲得一種精度高且兼顧效率的三維走時計算方法，本文將提出一種快速推進迎風雙線性插值法的三維地震波走時算法.如圖1所示，這是一種網(wǎng)格算法，和以往基于網(wǎng)格的雙線插值法、有限差分法、最短路徑追蹤法等類似，該算法主要包括三個方面的核心內(nèi)容：局部計算公式、局部實現(xiàn)策略以及整體實現(xiàn)策略.

圖1 算法的描述Fig.1 The description of the algorithm

如圖1所示，局部計算公式的建立可描述為如下問題：在計算的某一個時刻一個平面上的點A、B、C、D的走時值已知，且該平面內(nèi)的走時值為雙線性分布，則如何構建計算與該平面鄰近的網(wǎng)格節(jié)點F的走時值的計算公式；局部實現(xiàn)策略則可以描述為如下問題：在計算過程中，與走時值待求的F點相鄰近的、類似于平面ABCD的、走時值已知且分布滿足雙線性分布的平面可能會有很多種情況，而在不同情況時應該分別采用怎樣對應的不同的處理方法；最后，整體實現(xiàn)策略可以描述為如下問題：從源點S處給定的初始條件出發(fā)，應該采用怎樣的實現(xiàn)步驟或過程來逐次地計算整個計算空間內(nèi)所有網(wǎng)格節(jié)點的走時值.

與如上描述的三個問題相對應，這實際上涉及到算法的局部計算公式、算法的局部實現(xiàn)策略、算法的整體實現(xiàn)策略，所以接下來本文將分別詳細闡述這三方面問題.最后，還將通過算法的精度分析和計算實例來驗證算法的精度、正確性和有效性.

3 算法的局部計算公式：雙線性插值解析公式

為了計算三維地震波走時，首先需要推導建立局部走時計算公式.如上所述，局部走時計算公式的推導可描述為：如圖2所示，設點A、B、C、D的走時值已知，分別為：tA、tB、tC、tD；正方體網(wǎng)格的網(wǎng)格間距為h；假設平面ABCD上的走時值為雙線性分布；求F點的走時值tF的計算公式.

圖2 算法的局部計算公式：雙線性插值公式Fig.2 The local computation formulas of algorithm：the bilinear interpolation formulas

為了推導計算tF的公式，首先以點A為坐標原點，向量AB方向為x軸正方向，向量AD方向為y軸正方向，向量AF方向為z軸正方向，建立局部笛卡爾坐標系，則有點 A、B、C、D、F的坐標分別為A（0，0，0），B（h，0，0），C（h，h，0），D（0，h，0），F(xiàn)（0，0，h）；其次，假定到達F點的射線經(jīng)過了走時值分布為已知的平面ABCD內(nèi)的任意一點E，設E點的坐標為 E（x，y，0），因為 E 點被限定位于平面ABCD內(nèi)，則有0≤x≤h，0≤y≤h.過E點作平行于y軸的直線，該直線與線段AB、CD分別相交于E1和E2點，則點E1、E2的坐標分別為E1（x，0，0）、E2（x，h，0）.

根據(jù)平面ABCD上的走時值為雙線性分布的假設，線段AB、CD及E1E2上的走時值均為線性分布，則有

對上述（1）式進行整理可得

其中k1＝ （tB－tA）/h，k2＝ （tD－tA）/h，

k3＝ （tA＋tC－tB－tD）/h2.

設在局部正方體網(wǎng)格單元內(nèi)地震波傳播速度為均勻，并基于上述到達F點的射線經(jīng)過平面ABCD內(nèi)的E（x，y，0）點的假定，可得

其中s為當前正方體網(wǎng)格單元內(nèi)的地震波傳播速度的倒數(shù)（即慢度）.將式（2）代入（3）中，則有

根據(jù)Fermat原理，只要確定能使得（4）式取最小值的E（x，y，0）點的位置，即可獲得計算tF的公式，所以對（4）式tF分別關于變量x，y求偏導數(shù)，并令該偏導數(shù)為零，即

分析方程組（5）可知：只要求出該方程組的解，即可獲得E（x，y，0）點的位置坐標，然后再將其帶回到（4）式即可獲得計算tF的公式.但是，一些研究已表明方程組（5）實際上是一個超越方程組，無法求出相應的解析解（其可能解分布在一條雙曲線上），所以一些學者分別采用網(wǎng)格界面剖分法（張東等，2009）、快速插值法（李培明等，2013；梅勝全等，2010）、最速下降法（劉鋒等，2012）來獲得一個約束條件下接近于真實解的一個近似解.

與以往的雙線插值法的近似處理不同，在此筆者將充分利用雙線性假設的結論，來獲取方程組（5）的精確解析解，并構造性地證明計算tF的公式的正確性以及其背后滿足的地震波傳播波的物理事實.根據(jù)走時插值法的提出者Asakawa的線性假設（Asakawa and Kawanaka，1993），實際上即是平面波的假設，其闡述的是如果當前傳播的為平面波，則該波傳播經(jīng)過的二維空間中的走時分布在二維空間中的任意一條網(wǎng)格線上均為線性，而在三維空間中的一個平面上的走時分布則均為雙線性.同樣基于平面波的假設，并假定在局部的相對較小的一個網(wǎng)格單元內(nèi)地震波的傳播速度為勻速.如圖3所示，則無論從那個方向（假定箭頭方向為平面波的任意一個傳播方向）途經(jīng)平面ABCD均有tD－tA＝tC－tB，因為如圖3所示，設該平面波與平面ABCD相交于如圖3所示的一系列平行的虛線，其中經(jīng)過點A、B、C、D的虛線分別位于平面波tA、tB、tC、tD時刻的波前上，因為這些波前是相互平行的，并且網(wǎng)格四邊形ABCD為一正方形，所以必有tA時刻的波前傳到tD時刻花費的走時等于tB時刻的波前傳到tC時刻花費的走時，也即是tD－tA＝tC－tB.

圖3 平面波假設分析Fig.3 The analysis of plane wave assumption

基于上述雙線性假設，即平面波假設的進一步結論tD－tA＝tC－tB，可以得出k3＝0.將該結論帶回到（5）式中，則（5）可簡化為

再分析（6）式，若在k1、k2均不為零的情況下，則用下式除以上式再整理有：y＝k2x/k1，即方程組（6）的解實際上是位于一條直線上（而并非一條雙曲線上），并且將該直線代回方程組（6）中的任意一個方程，均可以獲得解析的x、y的解.不過要獲得各種情況下的方程組（6）的解，還得分情況討論k1、k2取不同值時，方程組（6）的解的情況：

①k1＝0且k2＝0時

此時，將k1、k2代回方程組（6）可得：x＝y(tǒng)＝0，再將該結論代回（4）式中可得計算tF的公式為

②k1＝0且k2＜0時

在此種情況下，將k1＝0代入方程組（6）可得x＝0，再將該結論代回方程組（6）中可將其簡化為

融入0≤y≤h的限定條件，求解（8）式可得

其中Δt＝tD－tA，且需滿足：（此是基于0≤y≤h的限定條件）.綜合x＝0的結論，可知 E點的坐標為0），再將其代入（4）式中即可獲得此種情況下計算tF的公式：

③k2＝0且k1＜0時

這種情況與情況②相似，不同之處僅在于交換了變量x與y的角色，所以同理情況②的推導過程可得：

其中：Δt＝tB－tA，且需滿足：（此是基于0≤x≤h的限定條件）.同時E點的坐標為計算tF的公式為：

④k1＜0且k2＜0

這種情況下，直接求解方程組（6），即將上述結論

代入方程組（6）的上式，并融入0≤x≤h、0≤y≤h的限定條件，則可得方程（6）的解：

其中，Δt1＝tB－tA、Δt2＝tD－tA，且需滿足：0≤（此是基于0≤x≤h、0≤y≤h的限定條件）.將（14）式代回（5）式中，即可獲得計算tF的公式：同時，此時E點的坐標為：

⑤k1＞0或k2＞0

這種情況下，因為上述將E點限定于平面ABCD內(nèi)的限定條件有0≤x≤h、0≤y≤h，所以很容易分析出，此時方程組（6）無解.這種無解的情況表面上表明：上述計算公式可能出現(xiàn)不穩(wěn)定的情況.不過，這種所謂的不穩(wěn)定的情況，在將上述算法融入一定的局部實現(xiàn)策略和算法的整體實現(xiàn)策略后，是不可能出現(xiàn)的，關于該問題將在后續(xù)內(nèi)容中闡述.

經(jīng)過以上的推導，實際上已經(jīng)獲得了當k1、k2取不同值的情況時，點E（x，y，0）的位置坐標和tF的計算公式.如圖4所示，下面將構造性地證明這些公式的正確性，以及它們對應的不同情況下的地震波傳播的物理事實.

圖4 雙線性插值公式的分析Fig.4 The analysis of bilinear interpolation formulas

上述情況①：當k1＝0且k2＝0時，即有tB＝tA，tD＝tA，也即tA＝tB＝tD，該式對應的物理事實為點A、B、D位于地震波傳播過程中的同一個等時面上.此外，根據(jù)不在同一條直線上的三點確定一個平面的數(shù)學事實和上述基于的平面波假設，則有與A、B、D點位于同一平面上的C點也位于該等時面上，所以可以進一步得出tA＝tB＝tC＝tD，也即平面ABCD為平面波的一個等時面.同時，根據(jù)局部笛卡爾坐標系可知，直線AF與平面ABCD垂直，再根據(jù)射線路徑和等時面相互垂直的物理規(guī)律，可以得出在這種情況下平面波沿著z軸正方向傳播，而經(jīng)過平面ABCD到達F點的射線與平面ABCD相交于A點，所以此時自然有：x＝y(tǒng)＝0，tF＝tA＋sh，也就是說上述情況①得的方程組（6）的解和tF的計算公式符合在這種情況下平面波的傳播規(guī)律，即是正確的.

上述情況②：當k1＝0且k2＜0時，即有tB＝tA，tD＜tA，此時再根據(jù)雙線性假設，因為tB＝tA，所以平面ABCD上的走時分布沿x方向沒有變化，而僅沿著y方向有變化，所以此時的插值問題實際上已簡化為二維插值問題，而E點則被固定到距F更近的線段AD上，即有x＝0.此外，基于二維空間中Asakawa提出的走時線性插值法同樣可以獲得其中Δt＝tD－tA.再分析tF的計算公式，實際上其可以變形為該方程實際上是二維Eikonal方程 （?t/?z）2＋（?t/?y）2＝s2在網(wǎng)格上以A點為中心的差分離散方程，也就是說此處得到的tF的計算公式，在局部上滿足Eikonal方程，所以其是正確的.同時上述0的限定條件還能保證tF的計算公式不會出現(xiàn)負數(shù)開平方的不穩(wěn)定情況.

上述情況③：當k2＝0且k1＜0時，此時的情況與情況②是類似的，也可以簡化為一個二維插值問題，不同之處僅在于坐標方向x與y交換了角色，所以在此不再重復闡述.

上述情況④：當k1＜0且k2＜0時，即有tB＜tA，tD＜tA，在這種情況下插值問題才真真變?yōu)橐粋€三維雙線性插值問題，在此條件下方程組（6）有解.求解該方程組即可獲得E點的坐標為x＝－h(huán)Δt1/同時將它們代回（5）式得到tF的計算公式：tF＝tA＋其中Δt1＝tB－tA、Δt2＝tD－tA. 同 樣 其 可 以 變 形 為該方程實際上是 三 維 Eikonal方 程 （?t/?x）2＋ （?t/?y）2＋（?t/?z）2＝s2在網(wǎng)格上以A點為中心的差分離散方程，也就是說此處得出的tF的計算公式在局部上是滿足Eikonal方程的，所以其是正確的.此外，再仔細分析E點的坐標可得：y＝kx，其中k＝Δt2/Δt1，也就是說，實際上E點是位于坐標平面xAy內(nèi)的一條直線上，該直線的斜率為k.k的取值實際上和已知條件tA、tB、tD的取值情況有關：k＝Δt2/Δt1＝ （tD－tA）/（tB－tA）.如圖4所示，tD＜tB時有k＞1，此時直線途徑ACD半?yún)^(qū)域，也即E點位于走時分布較小（因為雙線假設和tD＜tB）的區(qū)域.這實際上與Fermat原理和地震波傳播的基本規(guī)律是相符合的，即地震波總是從走時值相對更小的區(qū)域傳向下一個區(qū)域.而當tB＜tD時有k＜1，此時有如上同樣的結論.最后當tB＝tD時，則k＝1，此時E點位于正方形ABCD的對角線AC上.綜上所述，在情況④下得出的E點的坐標位置是滿足地震波的傳播規(guī)律和Fermat原理的，tF的計算公式在局部上滿足Eikonal方程.同時0≤2s2h2/3的限定條件，還能保證tF的計算公式不會出現(xiàn)負數(shù)開平方的不穩(wěn)定情況.

上述情況⑤：當k1＞0或k2＞0時，即有tB＞tA，tD＞tA，此時上述的結論為無解.實際上，再回顧方程組（6），此時并非方程組（6）無解，只是解的情況為x＜0，y＜0，即E點不在平面ABCD內(nèi).實際上，為了公式推導的方便和局部網(wǎng)格計算的限定，將E點限定于走時值已知的平面ABCD內(nèi)，而在該限定條件下x＜0，y＜0的解不能滿足該限定條件，所以致使方程組（6）無解.事實上，x＜0，y＜0的解的情況說明：在這種情況下E點位于其他的鄰近的網(wǎng)格單元內(nèi).同時tA＜tB、tA＜tD，再加上雙線性假設，暗含了其他相鄰的網(wǎng)格單元上的分布的走時值會比網(wǎng)格單元ABCD上的更小.而根據(jù)Fermat原理，此時計算tF的公式則應該通過其他網(wǎng)格單元上推導獲得，也即在計算時應該采用其他網(wǎng)格單元進行計算.該問題實際上引伸出了局部走時計算策略的問題：如圖5所示，在進行三維局部走時計算時，包圍被算點F周圍的類似于ABCD的網(wǎng)格單元實際上有24個，從幾何的角度講，到達F點的射線必與用于封閉F點的類似于ABCD的24個網(wǎng)格單元中的一個相交，而該交點正是我們上述需要求解的E點.當然要是在24個網(wǎng)格單元中逐次搜尋計算E點，必定是一種繁瑣且計算效率低下的算法.為了解決該問題，此處借鑒迎風差分的思想，提出了一種非常簡潔的迎風雙線性插值的局部實現(xiàn)策略，下面將詳細闡述該策略.

4 算法的局部實現(xiàn)策略：迎風雙線性插值法

如圖5所示，在三維走時計算的某一時刻，當前需要計算的是點F的走時值tF.如圖2所示，上述推導了網(wǎng)格單元ABCD走時分布為已知的情況下點F的走時值tF的計算公式.但實際上，在被算點F周圍的24個類似于ABCD的網(wǎng)格單元中，走時分布為已知的網(wǎng)格單元可能有很多個，那么究竟應該采用那個網(wǎng)格單元作為已知條件是合理的和穩(wěn)定的？

在此，引入迎風差分的基本思想（Sethian and Popovici，1999），該思想是指地震波總是從地震波走時值分布更小的區(qū)域向走時值未知的區(qū)域傳播，而在局部數(shù)值計算實現(xiàn)時就應該迎著走時值分布更小的區(qū)域去做差分代替微分的近似.迎風差分格式的基本思想實際上是Huygens原理和Fermat原理的綜合利用，主要體現(xiàn)在其隱含的將被算點周圍的走時值最小的網(wǎng)格節(jié)點當作子震源點，也正因為此迎風差分格式具有無條件穩(wěn)定性.同樣借鑒該思想，在此提出迎風雙線性插值法，該方法的核心借用迎風差分格式選取最小鄰近走時點的方法，選取被算點周圍的最小走時分布網(wǎng)格單元，然后再采用類似于圖2所示的推導的公式進行局部走時計算.總體上講，實際就是要構建相應的迎風雙線性插值公式.

如圖5a所示的局部正方體網(wǎng)格，以點F為中心的大正方體的六個面分別由以 A1、A2、A3、A4、A5、A6為中心的4個類似于圖2中的ABCD的網(wǎng)格單元組成，這樣就總共有24個類似于ABCD的網(wǎng)格單元.在當前計算點F的走時值tF的時刻，這24個網(wǎng)格單元上的26個網(wǎng)格節(jié)點的走時值可能為已知，可能為未知.由后續(xù)將闡述的算法的整體實現(xiàn)策略的初始化步驟（詳見后續(xù)闡述內(nèi)容的“算法的整體實現(xiàn)步驟：快速推進法”部分）可知，若網(wǎng)格節(jié)點的走時值為未知（即沒有被計算出來），其值為固定的一個非常大的數(shù)（這個數(shù)為人為規(guī)定的，其值遠大于所有網(wǎng)格節(jié)點最終可能算出的走時值），若網(wǎng)格節(jié)點的走時值被計算完成，其值即為最終的走時計算結果.在一般情況下，當前被算點F周圍的26個網(wǎng)格節(jié)點的走時值為一部分被計算完成，一部分未開始計算（如圖6a所示）.而在被計算完成的部分里面，先完成的走時值一定是小于后完成的走時值（該結論由后續(xù)內(nèi)容的窄帶技術的闡述可知），也就是說與點F鄰近的26個網(wǎng)格節(jié)點均有一個走時值，而用于計算tF的網(wǎng)格單元的走時分布應為最小.具體挑選最小走時分布區(qū)域的方法為：①確定與點F直接相鄰的6個網(wǎng)格節(jié)點的走時值最小的網(wǎng)格節(jié)點Amin，其中有 tAmin＝ min（tA1，tA2，tA3，tA4，tA5，tA6）；② 在 以Amin為中心的平面上分別確定其他兩個方向上的最小走時網(wǎng)格節(jié)點Bmin，Dmin，若設Amin＝A5則分別有tBmin＝min（tB51，tB52），tDmin＝ min（tD51，tD52）；③Amin，Bmin，Dmin所在的網(wǎng)格單元即為最終的類似于公式推導過程中的網(wǎng)格單元ABCD，與該網(wǎng)格單元相對應，tF的計算公式相應地變?yōu)?/p>

圖5 算法的局部實現(xiàn)策略：迎風雙線性插值Fig.5 The local implementation strategy of algorithm：the upwind bilinear interpolation

其中Δt1＝tBmin－tAmin，Δt2＝tDmin－tAmin.

公式（16）即為最終的迎風雙線性插值公式，與常規(guī)的雙線性插值公式相比，該公式實際上是一種優(yōu)化的公式.該公式隱含著一個優(yōu)化的局部實現(xiàn)策略，而該優(yōu)化的局部實現(xiàn)策略將Huygens原理和Fermat原理融入到算法的局部計算過程中，使得算法是在滿足地震波傳播的基本規(guī)律的基礎上實現(xiàn)的.與迎風差分格式類似，此處的迎風雙線性插值法也具有無條件穩(wěn)定性.

5 算法的整體實現(xiàn)策略：窄帶技術

上述內(nèi)容分別建立了算法的局部走時計算公式和局部實現(xiàn)策略，實際上要想完成整個計算區(qū)域的地震波走時計算，還需要擬定一個算法的整體實現(xiàn)步驟，即算法的整體實現(xiàn)策略.算法的整體實現(xiàn)策略往往決定著算法的穩(wěn)定性和計算效率.基于對計算效率、無條件穩(wěn)定性、對迎風雙線性插值的適應性等因素的綜合考慮，在此引入快速推進法中的窄帶技術（Sethian and Popovici，1999）作為算法的整體實現(xiàn)策略.

如圖6b所示，窄帶技術的核心思想是采用一個包含當前所有波前點的窄帶來近似模擬地震波前，通過窄帶的擴展演化來模擬地震波前的傳播過程，其充分考慮了地震波前傳播過程中的熵守恒理論，同時此處因為迎風差分思想的融入，所以窄帶技術具有無條件穩(wěn)定性.窄帶技術在實現(xiàn)過程中將計算空間中的網(wǎng)格節(jié)點分為三種類型：完成點，其走時計算已經(jīng)完成的點，在圖6b中用黑色填充的點表示；窄帶點，當前窄帶內(nèi)的點，也即近似波前上的網(wǎng)格節(jié)點，在圖6b中用灰色填充的點表示；遠離點，窄帶擴展還未到達的點，也即還未進行走時計算的網(wǎng)格節(jié)點，在圖6b中用白色填充的點表示.

圖6 算法的整體實現(xiàn)策略：窄帶技術Fig.6 The global implementation strategy of algorithm：the narrow band technique

窄帶技術的實現(xiàn)過程可以分為初始化和擴展計算兩個部分.在初始化時，震源點為唯一的完成點，其走時值為零；與震源點直接相鄰的正方體網(wǎng)格單元內(nèi)的網(wǎng)格節(jié)點為初始窄帶點，它們構成初始窄帶，即初始波前，它們的走時值通過局部的解析公式求取；其它所有剩下的計算空間中的網(wǎng)格節(jié)點為遠離點，其初始值設置為一個相對很大的數(shù)（遠大于最終可能計算出來的最大的走時值即可，這主要是為了實現(xiàn)上述迎風雙線性插值對最小走時分布區(qū)域的選定而設定的，這也是為什么此處的初始值不設為0的原因）.

擴展計算的步驟為：①通過比較找出窄帶內(nèi)走時值最小的窄帶點作為子震源點，并將其類型由窄帶點改為完成點；②對子震源直接鄰近的6個網(wǎng)格節(jié)點的類型進行判斷，如果其類型為遠離點則采用公式（16）計算它的走時值，并將其類型由遠離點改為窄帶點；如果其類型為窄帶點則再次計算它的走時值，并與其原來的走時值比較，較小者被保留；如果其類型為完成點則不對其作任何處理；③判斷窄帶內(nèi)的窄帶點是否還存在，是，則跳回步驟①重復循環(huán)計算，否則，計算結束.

上述窄帶技術的實現(xiàn)過程，與常規(guī)快速推進法中的窄帶技術的不同之處僅在于局部擴展計算采用的計算公式不一樣，常規(guī)快速推進法采用的是有限差分公式，而此處采用的是雙線性插值公式.由于將迎風差分的基本思想融入到了此處的雙線性插值公式中，提出了迎風雙線性插值法，所以此處的局部計算公式和整體實現(xiàn)策略的窄帶技術是完全兼容的.不過此處的迎風線性插值法和有限差分法在計算精度和效率方面卻有很大程度的差別，所以下面將對算法進行精度和效率分析.

6 算法的精度和效率分析

6.1 均勻介質(zhì)模型

為了驗證本文算法的計算精度和效率，選用一階精度有限差分法作對比分析.對比分析時，采用的模型為三維均勻介質(zhì)模型，模型大小為1.0km×1.0km×1.0km，計算時采用的網(wǎng)格間距分別為20.0、10.0、5.0m，震源點置于（0.5，0.5，0.0km）點處.如圖7a所示，誤差分析時分別提取三維模型的三個剖面：Z＝0.0km（圖7a的左圖）、Y＝0.5km（圖7a的中圖）、X－Y＝0.0km（圖7a的右圖）上的計算結果.

圖7b、圖7c、圖7d分別給出了當網(wǎng)格間距取10.0m時，剖面Z＝0.0km、Y＝0.5km、X－Y＝0.0km上，本文方法與一階精度有限差分法的誤差對比.分析這些誤差分布圖表明，在均勻介質(zhì)中本文方法能獲得正確的計算結果，并且本文方法具有很好的計算精度，其精度要遠高于一階精度有限差分法.

表1給出了均勻介質(zhì)中取不同網(wǎng)格間距時，一階精度有限差分法和本文算法在剖面Z＝0.0km、Y＝0.5km、X－Y＝0.0km上計算結果的平均相對誤差和整個模型計算的CPU耗時的對比，分析表1可以發(fā)現(xiàn)：①兩種方法的計算精度均受網(wǎng)格間距的影響，網(wǎng)格間距越小其計算精度越高；②當采用相同網(wǎng)格間距時，在所有剖面上本文的算法均有著比一階精度有限差分法更高的計算精度；③本文算法比一階精度有限差分法的計算效率高很多.通過仔細分析兩種算法實現(xiàn)過程可得出促使本文算法計算效率更高的主要原因是：在進行第5節(jié)“算法的整體實現(xiàn)策略”中闡述的步驟②的各網(wǎng)格節(jié)點走時的更新迭代計算時，本文算法精度更高，其更易快速獲取收斂于精確解附近的計算結果，所以其大幅度減少了更新重復計算的次數(shù)，進而大幅提高了計算效率.

表1 均勻介質(zhì)中本文方法與一階精度有限差分法的精度和效率對比Table 1 Comparing the accuracy and CPU time of our method with the 1st finite－difference method in homogeneous media

6.2 層狀介質(zhì)模型

為了分析地下速度不連續(xù)界面對走時計算的影響程度，將如圖7a所示的均勻介質(zhì)模型改為一個兩層的水平層狀介質(zhì)模型，其中分界面的深度為0.5km，上下兩層地震波傳播速度分別為1500.0m·s－1、1000.0m·s－1，其他模型參數(shù)和計算參數(shù)與圖7a所示的均勻介質(zhì)模型完全相同.表2給出了層狀介質(zhì)中取不同網(wǎng)格間距時，一階精度有限差分法和本文算法在剖面Z＝0.0km、Y＝0.5km、X－Y＝0.0km上計算結果的平均相對誤差和整個模型計算的CPU耗時的對比，對比分析表2和表1可以發(fā)現(xiàn)，算法在如上采用的層狀介質(zhì)模型中的計算精度略高于均勻介質(zhì)模型，而計算效率大體略低于均勻介質(zhì)模型，但總體上兩者只有細微的數(shù)值差別.該計算結果表明，無論采用一階精度有限差分法，還是采用本文算法進行走時計算，地下速度不連續(xù)界面對走時計算的影響程度很小.此外，基于該結論，在此不再重復顯示水平層狀介質(zhì)模型對應的與如圖7b—d類似的剖面上的誤差分布，因為它們基本上與圖7b—d相同，僅有極細微的差別.

表2 層狀介質(zhì)中本文方法與一階精度有限差分法的精度和效率對比Table 1 Comparing the accuracy and CPU time of our method with the 1st finite－difference method in layered media

圖7 算法的精度和效率分析Fig.7 The analysis of computational accuracy and efficiency for our method

圖8 三維復雜介質(zhì)模型中地震波走時計算結果Fig.8 The traveltime computational results in the 3D complex model

7 計算實例

為了驗證本文算法在面對復雜介質(zhì)時的穩(wěn)定性和有效性，引入一個由二維Marmousi模型生成的三維模型，該模型地下介質(zhì)分布非常復雜.如圖8a所示，模型的尺度為3.83km×3.83km×1.21km，計算時采用的網(wǎng)格間距為Δx＝Δy＝Δz＝10.0m，震源置于（0.0，0.0，0.0km）處.圖8b—e分別給出了Z＝0.0km、Y＝0.0km、X＝0.0km、X－Y＝0.0km這4個剖面上的地震波走時值的等時線分布，從這4個圖可以看出，本文算法在面對地下復雜介質(zhì)分布時，能夠有效地、無條件穩(wěn)定地獲得符合地震波傳播規(guī)律的計算結果.

在此需要特別提出的是，Z＝0.0km剖面上的計算結果.如圖8a和圖8b所示，Z＝0.0km剖面上的地震波傳播速度的分布實際上是均勻的，即均勻介質(zhì).按照地震波傳播規(guī)律，在均勻介質(zhì)中地震波走時的等時線應該為一系列的同心圓弧，但是圖8b顯示的計算結果并非如此，難道是本文算法在面對復雜介質(zhì)時出現(xiàn)了不適應和錯誤嗎？其實并非如此！實際上，其更能說明本文算法的有效性和正確性，因為本文采用的迎風雙線性插值技術和窄帶技術是嚴格以Fermat原理和Huygens原理為理論基礎的，同時其計算出來的走時值為初至波走時值.也正因為此，Z＝0.0km剖面上，遠離震源的區(qū)域開始接收到了來自下覆高速地層的首波，而這類型波的走時小于Z＝0.0km剖面上的直達波的走時，所以它們被作為最終計算結果.其實，這種現(xiàn)象的出現(xiàn)在折射波法勘探和三維地震觀測系統(tǒng)設計方面有重要的指導意義，有利于分析接收數(shù)據(jù)的特征和指導有針對性勘探的觀測系統(tǒng)設計.

8 結論與討論

為了獲得簡潔且精度高的三維地震波走時計算結果，筆者提出了一種新的快速推進迎風雙線性插值法，該方法包括如下核心技術環(huán)節(jié)：①充分利用平面波雙線插值的假設，采用分情況進行參數(shù)討論的方式，推導出了不同已知條件下雙線插值極小值方程的解析解，進而獲得了解析的局部走時計算公式；②引入迎風差分思想提出了迎風雙線插值的局部走時計算策略；③將上述局部走時計算公式和策略與常規(guī)快速推進法中的窄帶技術相結合提出了算法的整體實現(xiàn)步驟.由于采用了上述核心技術，新算法具有如下特征：①通過構造性的證明，新推導的雙線插值極小值方程的解析計算公式滿足Fermat原理和Eikonal方程；②迎風雙線插值的局部走時計算策略算法簡潔、易于編程實現(xiàn)且無條件穩(wěn)定；③在算法的整體實現(xiàn)過程中，因為引入了迎風差分的思想，所以迎風雙線性插值的局部走時計算策略能夠無條件穩(wěn)定地與窄帶技術兼容，進而有效地適應任意復雜介質(zhì).最后，算法的精度和效率分析驗證了上述提出的雙線插值極小值方程的解析公式是正確有效的，并且其計算結果具有比一階精度有限差分法高很多的計算精度和效率；同時計算實例還表明新提出的三維地震波走時計算方法能夠無條件穩(wěn)定地適應復雜介質(zhì).

Asakawa E，Kawanaka T.1993.Seismic ray tracing using linear traveltime interpolation.Geophysical Prospecting，41（1）：99－111.

Gao L，Li Y M，Chen X R，et al.2000.An attempt to seismic ray tracing with symplectic algorithm.Chinese J.Geophys.（in Chinese），43（3）：402－410，doi：10.3321/j.issn：0001－5733.2000.03.014.

Jing X L，Yang C C，Wang S Q.2007.A improved seismic reflection tomographic method.Chinese J.Geophys.（in Chinese），50（6）：1831－1836.

Julian B R，Gubbins D.1977.Three－dimensional seismic ray tracing.J.Geophys.，43：95－113.

Li P M，Mei S Q，Ma Q B.2013.An improved bilinear interpolation travel－time ray tracing method.Oil Geophysical Prospecting （in Chinese），48（4）：553－558.

Liu F，Zhang D，Yang Y，et al.2012.A fast numerical solution to minimum equation in 3－D seismic LTI ray tracing.Journal of Wuhan University （Natural Science Edition）（in Chinese），58（5）：395－400.

Mei S Q，Deng F，Zhong B S，et al.2010.The 3Dray tracing method base on the improved bilinear traveltime interpolation.Computing Techniques for Geophysical and Geochemical Exploration （in Chinese），32（2）：152－157.

Moser T J.1991.Shortest path calculation of seismic rays.Geophysics，56（1）：59－67.

Qin M Z，Chen J B.2000.Maslov asymptotic theory and symplectic algorithm.Chinese J.Geophys.（in Chinese），43（4）：522－533，doi：10.3321/j.issn：0001－5733.2000.04.013.

Qu C，Zhou H L，Zhao D P.2007.Deep structure beneath the west margin of Philippine Sea Plate and South China Sea from P and S wave travel time tomography.Chinese J.Geophys.（in Chinese），50（6）：1757－1768.

Sethian J A，Popovici A M.1999.3－D traveltime computation using the fast marching method.Geophysics，64（2）：516－523.

Sun J G.1998.On the limited aperture migration in two dimensions.Geophysics，63（3）：984－994.

Sun J G.1999.On the aperture effect in 3D Kirchhoff－type migration.Geophysical Prospecting，47（6）：1045－1076.

Sun J G.2000.Limited－aperture migration.Geophysics，65（2）：584－595.

Sun J G.2004.True－amplitude weight functions in 3D limitedaperture migration revisited.Geophysics，69（4）：1025－1036.

Sun J G，Sun Z Q，Han F X.2011.A finite difference scheme for solving the eikonal equation including surface topography.Geophysics，76（4）：T53－T63.

Vidale J E.1988.Finite－difference calculation of traveltimes.Bull.Seism.Soc.Am.，78（6）：2062－2076.

Vinje V，Iversen E，Gj?ystdal H.1993.Traveltime and amplitude estimation using wavefront construction.Geophysics，58（8）：1157－1166.

Wang H Z，Ma Z T.1999.Traveltime calculation in 3D media with arbitrary velocity distribution.∥ 69th Annual International Meeting，SEG，Expanded Abstracts，1778－1781.

Zhang D，F(xiàn)u X R，Yang Y，et al.2009.3D seismic ray tracing algorithm based on LTI and partition of grid interface.Chinese J.Geophys.（in Chinese），52（9）：2370－2376，doi：10.3969/j.issn.0001－5733.2009.09.023.

Zhang D，Zhang T T，Qiao Y F，et al.2013.A 3－D ray tracing method based on B－spline traveltime interpolation.Oil Geophysical Prospecting （in Chinese），48（4）：559－566.

Zhang M G，Jia Y G，Wang M Y，et al.2006.A global minimum traveltime ray tracing algorithm of wavefront expanding with interface points as secondary sources.Chinese J.Geophys.（in Chinese），49（4）：1169－1175.

附中文參考文獻

高亮，李幼銘，陳旭榮等.2000.地震射線辛幾何算法初探，地球物理 學 報，43（3）：402－410，doi：10.3321/j.issn：0001－5733.2000.03.014.

井西利，楊長春，王世清.2007.一種改進的地震反射層析成像方法.地球物理學報，50（6）：1831－1836.

李培明，梅勝全，馬青坡.2013.一種改進的雙線性插值射線追蹤方法.石油地球物理勘探，48（4）：553－558.

劉鋒，張東，楊艷等.2012.三維LTI射線追蹤極小值方程的快速數(shù)值解法.武漢大學學報（理學版），58（5）：395－400.

梅勝全，鄧飛，鐘本善等.2010.基于改進的雙線性走時插值的三維射線追蹤.物探化探計算技術，32（2）：152－157.

秦孟兆，陳景波.2000.Maslov漸近理論與辛幾何算法.地球物理學報，43（4）：522－533，doi：10.3321/j.issn：0001－5733.2000.04.013.

瞿辰，周蕙蘭，趙大鵬.2007.使用縱波和橫波走時層析成像研究菲律賓海板塊西邊緣帶和南海地區(qū)的深部結構.地球物理學報，50（6）：1757－1768.

張東，傅相如，楊艷等.2009.基于LTI和網(wǎng)格界面剖分的三維地震射線追蹤算法.地球物理學報，52（9）：2370－2376，doi：10.3969/j.issn.0001－5733.2009.09.023.

張東，張婷婷，喬友鋒等.2013.三維走時場B樣條插值射線追蹤方法.石油地球物理勘探，48（4）：559－566.

張美根，賈豫葛，王妙月等.2006.界面二次源波前擴展法全局最小走時射線追蹤技術.地球物理學報，49（4）：1169－1175.