基于二維三次卷積插值算法的辛幾何射線追蹤

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

1桂林理工大學地球科學學院，桂林 541004

2桂林理工大學廣西礦冶與環境科學實驗中心，桂林 541004

3中國科學院廣州地球化學研究所，廣州 510640

1 引言

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

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

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

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

2.1 射線追蹤系統及其初始化

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

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

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

本文中哈密爾頓函數H 選用可分系統哈密爾頓函數形式：

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

其中，ζ＝1／v，du為時間步長，dt為地震波走時，他們之間滿足關系

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

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

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

2.2 辛幾何算法

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

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

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

此外，對于式（7）的第三項，我們用數值積分的方法，可以計算得到該地震波射線的走時.

2.3 二維三次卷積插值

二維三次卷積插值是利用插值點g（x，z）周圍16個節點進行加權求和（圖2），該方法能夠更好地描述介質模型中該點的速度值，并且可以用中心差分格式計算出插值點的一階導數（韓復興等，2008）.

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

根據二維三次卷積插值方法，插值點g（x，z）處的值為

同時，利用中心差分方法可以計算得到插值點g（x，z）的x、z方向的導數為

其中

3 數值計算與分析

3.1 SAM與FMM算法的走時誤差分析

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

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

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

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

圖4 SAM的地震波射線空間誤差小圓表示射線的實際空間位置，實線為SAM計算得到的射線位置.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 辛幾何算法與龍格—庫塔算法的射線路徑分析

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

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

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

4 計算實例

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

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

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

根據具有不同初始入射角（φ0）射線計算結果，其哈密爾頓系統函數值H的誤差圖如圖7所示.

由圖7可知哈密爾頓系統函數的數值（H）保持在-0.0254左右，隨著計算步數的增大沒有大幅度的波動，比同階的龍格—庫塔算法效果更好.由此可見，SAM具有長時間運算保持哈密爾頓系統穩定性的特點.

5 結論

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

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

（2）對于同樣模型的計算結果，SAM的運算速度快，從而大幅度提高了射線追蹤的效率；

（3）SAM能夠保持哈密爾頓系統穩定性，同時也驗證了SAM在射線追蹤中的有效性和準確性.

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，Fedoriuk 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.

附中文參考文獻

陳景波，秦孟兆.2000.辛幾何算法在射線追蹤中的應用.數值計算與計算機應用，21（4），255-265.

馮康，秦孟兆.2003.哈密頓系統的辛幾何算法.杭州：浙江科學技術出版社.

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

韓復興，孫建國，楊昊.2008.基于二維三次卷積插值算法的波前構建射線追蹤.吉林大學學報：地球科學版，38（2）：336-340，346.

李世雄.2001.波動方程的高頻近似與辛幾何.北京：科學出版社.

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

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

王秀明，張海瀾，王東.2003.利用高階交錯網格有限差分法模擬地震波在非均勻孔隙介質中的傳播.地球物理學報，46（6）：842-849.

吳立明，許云，烏達巴拉.1995.高斯束射線法在二維非均勻介質復雜構造中的應用.地球物理學報，38（S1）：144-152.

徐濤，徐果明，高爾根等.2004.三維復雜介質的塊狀建模和試射射線追蹤.地球物理學報，47（6）：1118-1126.

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

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