陳連木，萬(wàn)永革＊，陳 鑫

（1.防災(zāi)科技學(xué)院，河北 燕郊 065201；2.廣西地震局，南寧 530022）

0 引言

地震學(xué)中的射線追蹤法是指給定發(fā)射點(diǎn)和接收點(diǎn)位置及介質(zhì)的波速，求從發(fā)射點(diǎn)到接收點(diǎn)的射線軌跡及其走時(shí)。作為一種快速有效的波場(chǎng)近似計(jì)算方法，地震射線追蹤對(duì)研究波在介質(zhì)中的傳播路徑和層析成像有重要的意義。傳統(tǒng)的射線追蹤方法通常包括邊值問題的彎曲法和初值問題的試射法。試射法在全局搜索接收器方面效果較好，而彎曲法的計(jì)算效率較高。為解決兩點(diǎn)間射線追蹤問題，許多學(xué)者提出了解決方法。高爾根、徐果明［1］提出了逐步迭代射線追蹤方法，但沒考慮對(duì)超低速區(qū)的處理；徐昇等［2］提出了射線追蹤的微變網(wǎng)格方法；張霖斌等［3］提出有限差分法射線追蹤；Langan［4］等在Cassell［5］的基礎(chǔ)上，設(shè)塊內(nèi)速度線性變化，并按Cassell的思路實(shí)現(xiàn)射線追蹤。本文主要是以一維速度模型追蹤水平層狀介質(zhì)2點(diǎn)的方法為原理，通過Matlab程序?qū)崿F(xiàn)該種方法的計(jì)算結(jié)果，并與其他方法進(jìn)行比較，指出其優(yōu)勢(shì)與不足。

1 射線追蹤方法

1.1 “打靶法”

打靶法是已知射線的起始點(diǎn)位置和給定的初始方向，通過擾動(dòng)初射方向來(lái)確定射線路徑。已知地震位置和臺(tái)站位置，算出震中距，利用試射法求得在震源層的θj，試射法的思路是：對(duì)于第1層的震源，不必用試射法，離源角通過可直接求出。對(duì)于位于其它層的地震，由于地下介質(zhì)速度隨深度會(huì)增大，臺(tái)站所在的位置一定在離源角的αmin＝arctan（最小值）和（最大值）所射出的射線到達(dá)地面的位置之間（h為震源深度，dz為震源到該層底部的距離）。我們以這2個(gè)邊界為初值，以它們的平均值為離源角，利用Δ ＝dz·。試算出射到地面的震中距（θj為震源層的入射角，νj為震源層的速度，νl為第1層的速度）。如果此震中距大于臺(tái)站的震中距，則說(shuō)明離源角取得大了，則令αmax＝α；反之，則說(shuō)明離源角取得小了，則令αmin＝α，進(jìn)行下一次迭代；再次以為離源角迭代計(jì)算。這樣反復(fù)迭代直到計(jì)算得到的震中距和實(shí)際臺(tái)站震中距在一定的誤差范圍內(nèi)，則此時(shí)的離源角θj就可直接射到臺(tái)站上［6］。

圖1 打靶法— 直達(dá)波傳播示意圖

圖2 求解p參數(shù)法路徑示意圖

1.2 求解p 參數(shù)法

求解p參數(shù)法是求解滿足震中距的方程式（1）的射線參數(shù)p。由于p參數(shù)為離源角正弦和該層速度的比值，知道了該層速度就知道了射線在這里如何彎曲。震中距的表達(dá)式為

震中距Δ0的表達(dá)式可以寫成關(guān)于p的非線性方程

該方程可采用牛頓迭代法進(jìn)行求解，該算法必須給定初始p參數(shù)作為第一次迭代的初始值，其一階迭代公式為

其中，f（pi）為（2）式的值，f′（pi）為（2）式的1階導(dǎo)數(shù)，i為迭代次數(shù)。這樣逐次迭代，直到得到的震中距和實(shí)際臺(tái)站震中距在一定的誤差范圍內(nèi)。

這種方法在震中距不大的情況下（Δ0＜24.6 km）都可以適用。但是當(dāng)震中距較大時(shí)，由于最大速度層（速度為νm）中的入射角θm趨于π／2 時(shí)，θm＝，即p很小的變化就會(huì)引起θm的巨大變化，從而導(dǎo)致震中距Δ的巨大變化而出現(xiàn)發(fā)散。因此，只能求解較小震中距的地震波走時(shí)。

1.3 改進(jìn)后的求解p 參數(shù)法

由于上面的迭代算法在震中距較大時(shí)發(fā)散，田玥［7］提出了改進(jìn)后的成層介質(zhì)中計(jì)算p參數(shù)的算法，其步驟如下：

改進(jìn)后的方法是通過將射線參數(shù)p轉(zhuǎn)換成別的參數(shù)q，q是最高速度層中射線走過的水平距離，νm是最高速度層中的速度。

p與q的轉(zhuǎn)換關(guān)系為

則q與震中距的關(guān)系為

同理，方程（5）可以寫成關(guān)于q的非線性方程

用牛頓法求解方程（6）的公式為

將得到的新參數(shù)q通過公式（4）轉(zhuǎn)換為射線參數(shù)p。用牛頓法求解過程需要對(duì)q0的初值進(jìn)行估計(jì)［7］。

2 兩點(diǎn)間射線追蹤的Matlab程序流程

編制程序的核心是水平層狀介質(zhì)2點(diǎn)間射線追蹤。我們根據(jù)給出的震中距求得p參數(shù)。為了增加結(jié)果的準(zhǔn)確性，分別運(yùn)用正、逆梯度速度模型以及含有高速層的速度模型進(jìn)行模擬。

程序流程見圖3：

圖3 平層介質(zhì)射線追蹤流程圖

程序運(yùn)行時(shí)應(yīng)注意以下問題：

在進(jìn)行逆梯度速度求解時(shí)應(yīng)注意：對(duì)于打靶法最大離源角αmax的初始值應(yīng)該設(shè)為，其中νj為震源層的速度，νm為最高速度層的速度。這是因?yàn)椋?dāng)時(shí)，Δ將表現(xiàn)為虛數(shù)。

3 三種方法的比較

在計(jì)算精度和速度結(jié)構(gòu)相同 的情況下，比較3種方法的所用時(shí)間。其中，直達(dá)波的震源深度設(shè)為19km，反射波的震源深度設(shè)為13km，計(jì)算精度均設(shè)為1.0×10－6。

首先，利用不同梯度速度結(jié)構(gòu)和假定的震源深度運(yùn)行射線追蹤程序，變換不同的震中距得到一系列追蹤的震中距—走時(shí)數(shù)據(jù)（圖4）?？梢钥闯觯谟?jì)算精度一樣，震中距一樣的情況下，無(wú)論是對(duì)于含有高速層的正速度梯度還是逆速度梯度的直達(dá)波、反射波來(lái)說(shuō)，都有以下特點(diǎn)：

表1 地殼速度模型

圖4 3種方法在不同速度模型下的計(jì)算效率

①于求解p參數(shù)法而言，該種方法只能求解震中距較?。ㄐ∮?4.6km）的p參數(shù)，這是因?yàn)楫?dāng)震中距較大時(shí)，由于最大速度層（速度為νm）中的入射角θm趨于π／2時(shí)，，即p很小的變化就會(huì)引起θm的巨大變化，從而導(dǎo)致震中距Δ的巨大變化而出現(xiàn)發(fā)散。

②對(duì)于直達(dá)波、反射波來(lái)說(shuō)，打靶法所用的時(shí)間總是最短（追蹤速度為改進(jìn)后的求解p參數(shù)法的2～3倍），求解p參數(shù)法次之，而改進(jìn)后的求解p參數(shù)法所用的時(shí)間相對(duì)較長(zhǎng)。這是因?yàn)楦倪M(jìn)后的求解p參數(shù)法在用牛頓迭代法求解q時(shí)需要進(jìn)行較為復(fù)雜的求導(dǎo)運(yùn)算，并且最后還需將參數(shù)q轉(zhuǎn)換為參數(shù)p，故所用時(shí)間相對(duì)較長(zhǎng)。

③3種方法的求解時(shí)間相差較小，而且求解時(shí)間也較快，最長(zhǎng)時(shí)間也不超過0.025s。

4 結(jié)論與討論

本文采用MATLAB 程序編制了打靶法、求解p參數(shù)法和改進(jìn)求解p參數(shù)法實(shí)現(xiàn)了水平成層介質(zhì)中的射線追蹤算法，并對(duì)算法的運(yùn)行效率進(jìn)行比較。可得出，直接求解p參數(shù)法所能追蹤到的距離較短，故該方法只適合小距離的追蹤；改進(jìn)后的求解p參數(shù)法能夠追蹤的距離較遠(yuǎn)，但是所用的時(shí)間相對(duì)較長(zhǎng)；打靶法的追蹤速度較快，并且追蹤的距離也較遠(yuǎn)。

本文是在水平層狀介質(zhì)中進(jìn)行的一維射線追蹤，對(duì)只考慮三維速度結(jié)構(gòu)的情況沒有進(jìn)行嘗試，但目 前 地 震 定 位［8－11］、近 震 地 殼 模 型 的 計(jì) 算［12－16］大 多采用成層速度模型進(jìn)行計(jì)算。因此，進(jìn)行該模型的探討對(duì)于理解地震走時(shí)計(jì)算和反演成層速度結(jié)構(gòu)均有重要的理論意義。

［1］ 高爾根，徐果明.二維速度隨機(jī)分布逐步［J］.地球物理學(xué)報(bào)，1996，39：302－308.

［2］ 徐昇，楊長(zhǎng)春，劉洪，等.射線追蹤的微變網(wǎng)絡(luò)方法［J］.地球物理學(xué)報(bào)，1996，39（1）：97－102.

［3］ 張霖斌，劉迎曦，趙振峰.有限差分法射線追蹤［J］.石油地球物理勘探，1993，28（6）：673－677，648.

［4］ Langan R T，Lerche I，Cutler R T.Tracing of rays though heterogeneous media：an accurate and efficient procedure［J］.Geophysics，1985，50（9）：1456－1465.

［5］ Cassell B R.A method for calculating synthetic seismograms in laterally varying media［J］.Geophysics，1982，69（2）：339－354.

［6］ 萬(wàn)永革.地震學(xué)導(dǎo)論［M］.北京：科學(xué)出版社，2015.（待出版）.

［7］ 田玥，陳曉非.水平層狀介質(zhì)中的快速兩點(diǎn)間射線追蹤方法［J］.地震學(xué)報(bào)，2005，27（02）：147－154.

［8］ 萬(wàn)永革，李鴻吉.遺傳算法在確定震源位置中的應(yīng)用［J］.地震地磁觀測(cè)與研究.1995，16（6）：1－7.

［9］ 萬(wàn)永革，盛書中，程萬(wàn)正，等.考慮到時(shí)誤差的地震定位算法及其在四川地區(qū)2001－2008年地震定位的應(yīng)用［J］.地震地質(zhì)，2012，34（1）：1－10.

［10］ 劉雙慶，謝靜，王熠熙.一種提高振動(dòng)源掃描算法信噪比的方法［J］.華北地震科學(xué)，2015，33（1）：46－51.

［11］ 劉芳.用sPn震相計(jì)算內(nèi)蒙古地震震源深度.大地測(cè)量與地球動(dòng)力學(xué).2010（30）增刊（II）：14－17.

［12］ Kissling E.Geotomography with local earthquake data［J］.Reviews of Geophysics，1988，26（4）：659－698.

［13］ Kissling E，Ellsworth W L，Eberhart－Phillips D，et al.Initial reference models in local earthquake tomography［J］.J Geophys Res，1994，99（B10）：19635－19646.

［14］ 賈麗華，李秀麗，李君，等.利用寬頻帶地震數(shù)據(jù)資料研究遼寧地區(qū)的地殼結(jié)構(gòu)［J］.華北地震科學(xué)，2013，31（3）：1－7.

［15］ 孫譯，賴曉玲.利用斷層圍陷波資料研究汶川MS8.0地震構(gòu)造特征［J］.華北地震科學(xué)，2014，32（3）：1－4.

［16］ 陳立軍，胡奉湘，陳曉逢.全球地震柱的地震層析成像證據(jù)［J］.華南地震，2013，33（4）：1－10.