李川，王有學(xué)，2＊，何曉玲，劉榮平，贠鵬，熊彬，許繼峰，3

1桂林理工大學(xué)地球科學(xué)學(xué)院，桂林 541004

2桂林理工大學(xué)廣西礦冶與環(huán)境科學(xué)實(shí)驗(yàn)中心，桂林 541004

3中國(guó)科學(xué)院廣州地球化學(xué)研究所，廣州 510640

1 引言

射線追蹤是地震學(xué)反演問題的基礎(chǔ)，地震波的走時(shí)和射線路徑廣泛用于走時(shí)反演、層析成像、反射偏移及地震定位等領(lǐng)域中，而地震波射線追蹤的速度及精度決定了反演的速度與質(zhì)量.目前，廣泛使用的射線追蹤方法包括基于射線理論的打靶法（Julian and Gubbins，1977；徐濤等，2004）、彎曲法（Moser et al.，2004）、高 斯 射 線 束 算 法 （and；吳立明等，1995）等，基于網(wǎng)格單元擴(kuò)展的有限差分解程函方程法（Vidale，1988；Cao Shun-hua et al.，1991；王秀明，2003）、最短路徑算法（Moser，1991；張美根等，2006），以及結(jié)合射線和網(wǎng)格單元擴(kuò)展的波前構(gòu)造法（Vinje et al.，1993；趙連鋒等，2003），求解程函方程的快速行進(jìn)法（FMMFast Marching Method）（Sethian，1999；Rawlinson，2003）等.在射線追蹤的過程中，地震波的傳播滿足波動(dòng)方程，我們?cè)谝话闱闆r下很難得到精確解.最早針對(duì)波動(dòng)方程求解的方法是一種高頻近似法即幾何光學(xué)法，幾何光學(xué)法將波動(dòng)現(xiàn)象轉(zhuǎn)化成為射線理論，并用射線管的概念來解釋射線的傳播機(jī)制.但是，對(duì)于地下速度結(jié)構(gòu)存在微小擾動(dòng)的成像過程中，普通的射線理論會(huì)在焦散處出現(xiàn)奇點(diǎn)，這使得幾何漸近射線理論失效（Chapman and Drummond 1982，Chapman 1985）.前蘇聯(lián)學(xué)者M(jìn)aslov根據(jù)方程變換和Fourier積分算子理論，引入具有相同維數(shù)的慢度向量空間，然后再利用Fourier逆變換回到原來的空間，并且引進(jìn)正則（canonical）算子，在相空間中引入適當(dāng)?shù)男羶?nèi)積后成為辛空間，從而可以得到焦散處附近有效地高頻近似解（Maslov and Fedoriuk，1981；李世雄，2001）.

哈密爾頓（Hamiltonian）系統(tǒng)是描述各種守恒的物理和力學(xué)過程的三種基本形式（牛頓力學(xué)、拉格朗日力學(xué)及哈密爾頓力學(xué)）之一，而辛幾何則是該系統(tǒng)的數(shù)學(xué)基礎(chǔ)，并且哈密爾頓系統(tǒng)具有辛結(jié)構(gòu)不變性和能量守恒性.作為彈性力學(xué)中的地震波射線追蹤問題，當(dāng)然也可以用哈密爾頓系統(tǒng)進(jìn)行描述（羅明秋等，2001）.早在1984年馮康先生首先提出了求解哈密爾頓動(dòng)力學(xué)體系的辛幾何算法（SAMSymplectic Algorithm Method）（馮康和秦孟兆，2003），并應(yīng)用于哈密爾頓系統(tǒng)的求解（秦孟兆和陳景波，2000）.陳景波 與 秦孟兆 （2000，2001）在Maslov研究的基礎(chǔ)上，采用辛幾何算法這一專門針對(duì)Hamilton系統(tǒng)的數(shù)值方法，利用Maslov漸近理論對(duì)地震波波場(chǎng)進(jìn)行了數(shù)值模擬，從SAM觀點(diǎn)看來，射線是相空間中的特征線在物理空間上的投影，能有效地彌補(bǔ)傳統(tǒng)方法上的不足（李世雄，2001）.

本文利用SAM算法，結(jié)合二維三次卷積（韓復(fù)興等，2008；Keys，1981），對(duì)復(fù)雜速度介質(zhì)模型的地震波射線追蹤進(jìn)行研究，并同龍格—庫(kù)塔數(shù)值微分算法相比較，SAM可以保證哈密爾頓不變量守恒，具有不隨運(yùn)算時(shí)間增加而減弱的特點(diǎn).

2 辛幾何算法及二維三次卷積插值

2.1 射線追蹤系統(tǒng)及其初始化

程函方程（1）式也寫成如下形式：

在笛卡爾坐標(biāo)系中，我們用下標(biāo)i＝（1，2，3）分別表示（x，y，z）.式（2）中t＝t（xi）表示地震波走時(shí)（程函），pi是介質(zhì)的慢度向量，Δ

程函方程（2）是t（xi）的一階非線性的偏微分方程，并且滿足如下的哈密爾頓系統(tǒng)形式：

本文中哈密爾頓函數(shù)H 選用可分系統(tǒng)哈密爾頓函數(shù)形式：

并且偏微分方程式（3）的在三維坐標(biāo)系下的特征向量解形式為

其中，ζ＝1／v，du為時(shí)間步長(zhǎng)，dt為地震波走時(shí)，他們之間滿足關(guān)系

將式（4）帶入式（5），并用變量τ代替時(shí)間步長(zhǎng)變量u，便可得到需要求解的7個(gè)射線方程組：

在各向同性的三維介質(zhì)中，地震波射線路徑和它的走時(shí)都滿足方程（7），其初始條件的取值由震源S確定.設(shè)i0表示震源處射線與垂直方向之間的夾角，φ0為震源處射線在水平面投影與x1方向的夾角（圖1），那么，由式（2）可得震源位置、變量處pi0及地震波走時(shí)的初始值為

圖1 震源S處射線參數(shù)示意圖Fig.1 The diagram of seismic ray parameters at source S

2.2 辛幾何算法

為了求解哈密爾頓系統(tǒng)式（4）中射線的空間位置，我們采用辛幾何數(shù)值積分方法（Feng and Qin，2010；Thijssen，1999）.

對(duì)于形如H＝U（p）+W（x）形式的可分哈密爾頓系統(tǒng)，式（7）的k階SAM下的解為

待定系數(shù)ak、bk（k＝1，2，3，…）的選取與射線追蹤的精度有很大的關(guān)系，本文選取k＝4（4階），其顯式的辛格式系數(shù)為

此外，對(duì)于式（7）的第三項(xiàng)，我們用數(shù)值積分的方法，可以計(jì)算得到該地震波射線的走時(shí).

2.3 二維三次卷積插值

二維三次卷積插值是利用插值點(diǎn)g（x，z）周圍16個(gè)節(jié)點(diǎn)進(jìn)行加權(quán)求和（圖2），該方法能夠更好地描述介質(zhì)模型中該點(diǎn)的速度值，并且可以用中心差分格式計(jì)算出插值點(diǎn)的一階導(dǎo)數(shù)（韓復(fù)興等，2008）.

圖2 二維三次卷積插值示意圖Fig.2 The diagram of bi-cubic interpolation

根據(jù)二維三次卷積插值方法，插值點(diǎn)g（x，z）處的值為

同時(shí)，利用中心差分方法可以計(jì)算得到插值點(diǎn)g（x，z）的x、z方向的導(dǎo)數(shù)為

其中

3 數(shù)值計(jì)算與分析

3.1 SAM與FMM算法的走時(shí)誤差分析

在速度為2km／s的均勻各向同性二維（x-z）模型中，設(shè)震源S的坐標(biāo)為（50，50），分別用FMM和SAM進(jìn)行地震波射線追蹤，在SAM中，初始角i0＝90°，并且φ0以步長(zhǎng)為18°，時(shí)間步長(zhǎng)τ＝0.01s，計(jì)算得到φ0＝0°～360°的20條射線.然后，利用有限差分方法求解式（7）的第三項(xiàng)，計(jì)算各條地震波射線的走時(shí)，并將其與射線精確走時(shí)的差值取絕對(duì)值，即得到地震波射線走時(shí)的誤差（Δt），其結(jié)果如圖3所示.

圖3 射線追蹤的走時(shí)誤差分析.（a）FMM；（b）SAMFig.3 The traveltimes error in raytracing.（a）FMM；（b）SAM

由圖3可知，用SAM得到的地震波射線走時(shí)最大誤差為1.4×10-6s，比FMM的精度要高出5個(gè)數(shù)量級(jí).

如果將震源S置于點(diǎn)（50，0），初始角i0＝90°，并且φ0以步長(zhǎng)20°、時(shí)間步長(zhǎng)τ＝0.01s，用SAM計(jì)算得到φ0＝10°～170°的地震波射線的空間誤差（見圖4）.

圖4 SAM的地震波射線空間誤差小圓表示射線的實(shí)際空間位置，實(shí)線為SAM計(jì)算得到的射線位置.Fig.4 The special error of seismic rays by using SAM Circles are the accurate rays，solid lines are the rays computed by SAM.

3.2 辛幾何算法與龍格—庫(kù)塔算法的射線路徑分析

在速度隨深度線性變化的各向同性介質(zhì)中，設(shè)v＝1.8+0.3 z，震源點(diǎn)S位于（0，4），初始入射角i0＝90°，時(shí)間步長(zhǎng)τ＝0.01s，φ0以步長(zhǎng)5°，分別用SAM和龍格—庫(kù)塔算法計(jì)算得到φ0＝-60°～-30°的地震波射線路徑，結(jié)果如圖5所示.

圖5 SAM算法（小圓）與龍格—庫(kù)塔算法（實(shí)線）計(jì)算得到的射線路徑對(duì)比Fig.5 The raypaths computed by SAM（circle）and Runge-Kutta algorithm（solid line）

從圖5中我們可以得出，SAM用于地震射線追蹤是可行的，其中龍格—庫(kù)塔算法CPU耗時(shí)1.2s，而用SAM計(jì)算機(jī)只需0.0624s，對(duì)于運(yùn)算量巨大的地震層析成像，SAM會(huì)大幅度地減少計(jì)算時(shí)間.同時(shí)，圖5也說明在計(jì)算步數(shù)不太多的情況下，采用龍格—庫(kù)塔算法求解哈密爾頓方程組（5）的初值問題，已可滿足要求.

4 計(jì)算實(shí)例

對(duì)于圖6所示的二維Marmousi速度模型，設(shè)震源點(diǎn)S位于（250，0），初始入射角i0＝90°，時(shí)間步長(zhǎng)τ＝0.01s，φ0以步長(zhǎng)10°，將四階（k＝4）顯式SAM與二維三次卷積插值算法相結(jié)合，對(duì)φ0＝0°～180°之間的地震波進(jìn)行射線追蹤，其結(jié)果如圖6所示.

圖6 二維Marmousi模型及射線路徑Fig.6 2-D Marmousi model and seismic raypaths

圖7 四階顯式SAM計(jì)算的沿射線軌跡Hamilton函數(shù)值Fig.7 The Hamiltonian（H）value along raypaths by 4-order explicit SAM

根據(jù)具有不同初始入射角（φ0）射線計(jì)算結(jié)果，其哈密爾頓系統(tǒng)函數(shù)值H的誤差圖如圖7所示.

由圖7可知哈密爾頓系統(tǒng)函數(shù)的數(shù)值（H）保持在-0.0254左右，隨著計(jì)算步數(shù)的增大沒有大幅度的波動(dòng)，比同階的龍格—庫(kù)塔算法效果更好.由此可見，SAM具有長(zhǎng)時(shí)間運(yùn)算保持哈密爾頓系統(tǒng)穩(wěn)定性的特點(diǎn).

5 結(jié)論

在利用SAM進(jìn)行地震波射線追蹤的過程中，由于采用了辛幾何射線追蹤算法和二維三次卷積插值，有效地保證了射線追蹤的準(zhǔn)確性和穩(wěn)定性.根據(jù)數(shù)值計(jì)算及分析，可以得出如下結(jié)論：

（1）運(yùn)用SAM進(jìn)行地震波射線追蹤，可以獲得高精度的波前，進(jìn)一步提高了射線空間位置的準(zhǔn)確性；

（2）對(duì)于同樣模型的計(jì)算結(jié)果，SAM的運(yùn)算速度快，從而大幅度提高了射線追蹤的效率；

（3）SAM能夠保持哈密爾頓系統(tǒng)穩(wěn)定性，同時(shí)也驗(yàn)證了SAM在射線追蹤中的有效性和準(zhǔn)確性.

Cao S H，Greenhalgh S.1994.Finite-difference solution of the eikonal equation using an efficient，first-arrival，wavefront tracking scheme.Geophysics，59（4）：632-643.

Chapman C H，Drummond R.1982.Body-wave seismograms in inhomogeneous media using Maslov asymptotic theory.Bull.Seism.Soc.Am.，72（6B）：S277-S317.

Chapman C.1985.Ray theory and its extensions：WKBJ and Maslov seismograms.J.Geophys.，58（1）：27-43.

Chen J B，Qing M Z.2000.Ray tracing by symplectic algorithm.Journal On Numerical Methods And Computer Applications（in Chinese），21（4）：255-265.

Comer R P.1984.Rapid seismic ray tracing in a spherically symmetric earth via interpolation of rays.Bull.seism.Soc.Am.，74（2）：479-492.

Feng K，Qin M Z.2003.Symplectic Geometric Algorithm for Hamiltonian Systems（in Chinese）.Hangzhou：Zhejiang Science and Technology Press.

Feng K，Qin M Z.2010.Symplectic Geometric Algorithms for Hamiltonian Systems.New York：Springer.

Gao L，Li Y M，Chen X R，etal.2000.An attempt to seismic ray tracing with symplectic algorithm.Chinese Journal of Geophysics（in Chinese），43（3）：402-410.

Han F X，Sun J G，Yang H.2008.Ray-tracing of wavefront construction by bicubic convolution interpolation.Journal of Jilin University：Earth Science Edition（in Chinese），38（2）：336-340，346.

Yoshida H.1990.Construction of higher order symplectic integrators.Physics Letters A，150（5-7）：262-268.

Julian B R，Gubbins D.1977.Three-dimensional seismic ray tracing.J.Geophys.，43（1）：95-113.

Li S X.2001.High-frequency Approximation of Wave Equation and Symplectic Algorithm.Beijing：Science Press.

Luo M Q，Liu H，Li Y M.2001.Hamiltonian description and symplectic method of seismic wave propagation.Chinese J.Geophys.（in Chinese），44（1）：120-128.

Maslov V P，F(xiàn)edoriuk M V.1981.Semi-Classical Approximation in Quantum Mechanics.Hingham，Mass：D Reidel Publishing.

Moser T J，Nolet G，Snieder R.1992.Ray bending revisited.Bulletin of the Seismological Society of America，82（1）：259-288.

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

Qing 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

Rawlinson N.2003.Seismic Ray Tracing and Wavefront Tracking in Laterally Heterogeneous Media.Canberra：Australian National University.

Keys R G.1981.Cubic convolution interpolation for digital image processing.IEEE Transaction on Acoustics，Speech，and Signal Processing，29（6）：1153-1160.

Sethian J A.1999.Fast marching methods.SIAM Review，41（2）：199-235.

Thijssen J M.1999.Computational Physics.Cambridge：Cambridge University Press.

Vidale J E.1988.Finite-difference calculations 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 X M，Zhang H L，Wang D.2003.Modelling of seismic wave propagation in heterogeneous poroelastic media using a highorder staggered finite-difference method.ChineseJournalof Geophysics（in Chinese），46（6）：842-849.

Wu L M.Xu Y，Wudabala.Applied study of Gaussian beam method in 2-D inhomogeneous media and laterally varying layered structures.Chinese J.Geophys.（in Chinese），38（S1）：144-152.

Xu T，Xu G M，Gao E G，etal.2004.Block modeling and shooting ray tracing in complex 3-D media.Chinese J.Geophys.（in Chinese），47（6）：1118-1126.

Zhang L B，Rector J W，Michael Hoversten G.2005.Eikonal solver in the celerity domain.Geophys.J.Int.（in Chinese），162（1）：1-8.

Zhang M G，Cheng B J，Li X F，etal.2006.A fast algorithm of shortest path ray tracing.Chinese J.Geophys.（in Chinese），49（5）：1467-1474.

Zhao L F，Zhu J S，Cao J X，etal.2003.Ray tracing using ordinal wave-front reconstruction method.Chinese J.Geophys.（in Chinese），46（3）：415-420.

附中文參考文獻(xiàn)

陳景波，秦孟兆.2000.辛幾何算法在射線追蹤中的應(yīng)用.數(shù)值計(jì)算與計(jì)算機(jī)應(yīng)用，21（4），255-265.

馮康，秦孟兆.2003.哈密頓系統(tǒng)的辛幾何算法.杭州：浙江科學(xué)技術(shù)出版社.

高亮，李幼銘，陳旭榮等.2000.地震射線辛幾何算法初探.地球物理學(xué)報(bào)，43（3）：402-410.

韓復(fù)興，孫建國(guó)，楊昊.2008.基于二維三次卷積插值算法的波前構(gòu)建射線追蹤.吉林大學(xué)學(xué)報(bào)：地球科學(xué)版，38（2）：336-340，346.

李世雄.2001.波動(dòng)方程的高頻近似與辛幾何.北京：科學(xué)出版社.

羅明秋，劉洪，李幼明.2001.地震波傳播的哈密頓表述及辛幾何算法.地球物理學(xué)報(bào)，44（1）：120-128.

秦孟兆，陳景波.2000.Maslov漸近理論與辛幾何算法.地球物理學(xué)報(bào)，43（4）：522-533，doi：10.3321／j.issn：0001-5733.2000.04.013.

王秀明，張海瀾，王東.2003.利用高階交錯(cuò)網(wǎng)格有限差分法模擬地震波在非均勻孔隙介質(zhì)中的傳播.地球物理學(xué)報(bào)，46（6）：842-849.

吳立明，許云，烏達(dá)巴拉.1995.高斯束射線法在二維非均勻介質(zhì)復(fù)雜構(gòu)造中的應(yīng)用.地球物理學(xué)報(bào)，38（S1）：144-152.

徐濤，徐果明，高爾根等.2004.三維復(fù)雜介質(zhì)的塊狀建模和試射射線追蹤.地球物理學(xué)報(bào)，47（6）：1118-1126.

張美根，程冰潔，李小凡等.2006.一種最短路徑射線追蹤的快速算法.地球物理學(xué)報(bào)，49（5）：1467-1474.

趙連鋒，朱介壽，曹俊興等.2003.有序波前重建法的射線追蹤.地球物理學(xué)報(bào)，46（3）：415-420.