射線追蹤方法定位近震震源空間位置

張潛 尹耿 王玉石 林國良

摘要：介紹了一種基于地震波射線理論、不需要先驗速度結構的定位震源空間位置的方法——射線追蹤法。該方法利用地震射線會反向匯聚到震源的幾何性質，對地面介質分層，并用網格搜索反演方法計算分層速度參數，進而確定震源位置；然后逐步增加分層，重復前面的步驟，對震源位置進行校正，最終得到一個較為精確的定位結果。通過虛擬事件對該方法進行了測試，測試結果表明，該方法理論上能夠用于近震空間位置的定位；使用該方法對一個真實震側進行定位，并和云南地震臺網的定位的結果進行了比較，兩者較為接近。

關鍵詞：地震定位；無速度結構；射線理論；P波偏振

中圖分類號：P3153文獻標識碼：A文章編號：1000-0666（2017）02-0203-08

0引言

震源空間位置的確定是地震學的經典問題——地震定位的一部分，也是地震學研究的基礎。地震定位，無論是對于地球物理學研究，如地球內部的環境構造、地震的發生機制，還是對于震后救援工作的需求，如震害范圍和程度的預估、地震趨勢的預測以及地震預警，都有著無法替代的重要性。因此，地震定位的方法以及提高地震定位精度和速度的方法，一直都是地震學家的研究目標。

早期的地震定位主要是通過幾何作圖完成的。Geiger（1912）根據各個觀測點到時差最小化的原則，將地震波走時方程組線性化，然后通過最小二乘法來定位地震。這也為以后使用計算機來定位地震提供了基礎。20世紀70年代后，伴隨著現代計算機技術的發展，地震學家們基于Geiger的理論，給出了一系列地震地位的程序和算法（Lee，Lahr，1975；Klein，1978；Lienert et al，1986；Nelson，Vidale，1990），并根據實際應用中出現的問題，提出了相應的改進方法（Lienert et al，1986；Prugger，Gendzwill，1988）。隨后，多事件定位方法（Douglas，1976；Crosson，1976；Spence，1980；Waldhauser，Ellsworth，2000）、基于距離殘差的空間域定位方法（Romney，1957；Lomnitze，1977）和非線性的定位方法（Prugger，Gendzwill，1988；Thurber，1985；Tarantola，Valette，1982）也相繼出現。震相、振幅，甚至是全波形的信息都被用于震源位置的確定（Engdahl，et al，1998；Warren，Shearer，2005；羅艷等，2013；Zhao，Helmberger，1994；Zhu，Helmberger，1996）。

在地震定位的過程中，影響結果精度的因素大致有：觀測系統的結構、速度結構的選取或假設、震相識別以及觀測數據（包含觀測點的空間位置）的讀取和計算過程中產生的誤差（朱元清，趙仲和，1997；田玥，陳曉非，2002）。絕大多數時間域和空間域的定位方法，都需要使用某一組先驗的速度結構。因此，該種方法中速度結構的選取，將直接影響到結果的誤差水平。而通過波形信息來確定震源位置時，震相識別、震源機制解和地球結構模型的可靠性，將決定定位結果的可靠性。

本文將介紹一種基于地震波射線理論，且不需要先驗速度結構的定位近震震源位置的方法——射線追蹤法，并使用虛擬數據和一個真實震例對該方法進行了檢驗。[HJ2mm]

1基本原理

假設地球介質是水平均勻且各向同性的分層介質。根據惠更斯原理和地震波射線理論，一束地震波在傳播時，射線和分層界面法線的夾角滿足方程[KH*2]

[SX（]v1sinθ1[SX）]=……=[SX（]vjsinθj[SX）]=……=[SX（]vnsinθn[SX）]=const[JY]（1）[KH*3/4D]

其中：vj和θj是射線在第j層的波速和夾角（震源所在位置為第n層，地表為第一層上邊界）。將地表處的夾角和波速作為已知條件，則[KH*2]

θj=arcsin（[SX（]vjv1[SX）]sinθ1）[JY]（2）[KH*3/4D]

從（2）式可以看出，射線與法線的夾角只與地表入射角θn和當前波速與地表波速的比值vj/v1（記為rj）有關。另一方面，對于一束從震源發出的射線i，其在第j層介質中的傳播，符合方程[KH*1]

sij=Δhj/tanθij[JY]（3）[KH*1D]

式中：sij和θij分別代表射線i在第j層介質中的水平位移和與法線的夾角；Δhj表示第j層的厚度。同時，對于射線i，滿足方程[KH*1]

θij=arcsin（rjsinθi1）[JY]（4）[KH*1]

結合式（3），可得：[KH*1]

sij=Δhj/tan（arcsin（rjsinθi1））[JY]（5）[KH*1]

因此，震源到射線接收點（即觀測點）之間的水平位移為：[KH*2]

[JP2]Si[KG-*4]=[KG-*2]∑[DD（]nj=1[DD）]sij[KG-*4]=[KG-*2]∑[DD（]nj=1[DD）]Δhj/tan（arcsin（rjsinθi1））[JY]（6a）[KH*2D][JP]

或

[JP3]Si[KG-*4]=[KG-*4]∫[DD（]H[KG*4]i0[DD）]si（h）dh[KG-*4]=[KG-*4]∫[DD（]H[KG*4]i0[DD）]1/tan（arcsin（r（h）sinθi1））dh[JP][KH-*4][JX*1][JY]（6b）[KH*1D]

其中：Hi代表震源深度。θi1和Pg波的偏振方向相關，可以通過觀測得到；而Si和Hi代表著觀測點和震源空間位置的矢量差。因此，我們可以得到關于震源空間位置F（x，y，h）和r（r1，r2，…rn）的方程（組）：[KH*1]

[JB（{][KF（]（xi-y）2+（yi-y）2[KF）]-Si=0[KG*2]（hi-h）-H=0[JB）][JY]（7）[KH*1D]

其中：（xi，yi，hi）為觀測點i的空間位置。本文所介紹的方法，就是要找到一組（F，r），使（7）式成立，或者使等式兩邊的差值盡量小，從而實現震源的定位。

[KG（0.15mm]由于式（6）和（7）含有三角函數，不能采用線性方程組的方法來求解，而必須尋找一種合適的反演方法來確定（F，r）。網格搜索法是一種不錯的方法，但是其計算量會隨著參數的增加而呈指數增長。如果待求解的參數越多，則計算所花時間越長。因此，我們采用了一種“分步”的網格搜索法來定位震源，待求參數變少，計算效率得到了提升。圖1是對這種方法的一個簡單的展示，由于看起來像追溯射線的路徑回到震源，我們稱該方法為“射線追蹤法”，其操作步驟如下：[KG）]

（1）根據各個觀測點接收到的Pg波的偏振方向，反向生成一組射線；尋找到各條射線距離最短的一點F1，并將其作為初始定位的結果（F1，r1）（r1=1）。

（2）將介質平均分成2層，每層的厚度為h1/2，r1取步驟（1）的結果。使用網格搜索法，尋找第二層的波速比r2，以及相應的點F2，使F2到根據（r1，r2）做出的反向射線的距離和最小。將（F2，r2）作為最新的結果；

（3）重復步驟（2），直到最后得到的結果（Fn，rn）滿足一定的收斂條件。Fn即為最終的定位結果。

通過射線追蹤法，我們一步步將介質的分層數從2層增加至n層，且每一次分層都只用搜索新增層的波速比，從而大大減少了計算量。最終，我們會得到震源的位置Fn以及1組和各層相對應的波速比（r1，r2，…，rn）。

2測試數據

為了驗證此地震定位方法的準確性，本文虛擬了多組數據對該方法進行測試。測試數據生成[FL）][SD1，1]

[JZ][XC張潛1.TIF][KH*2/3D][HT5K][JZ]

圖1射線追蹤法步驟的簡單演示[JZ]

Fig1Demonstration of ray tracing method[KH*1]

的方法如下：

（1）選擇1個一維速度模型作為基準模型，假定震中位置為0，震源深度為10 km。

（2）隨機生成一組出射角為（θ1n，…，θmn）的射線，根據步驟1的速度模型和震源位置，計算各條射線到達地表（深度為0）處的角度（θ11，…，θm1）和震中距（S1，…，Sm）。

（3）為（S1，…，Sm）和（θ11，…，θm1）引入誤差，組成一組數據，分別對應m個虛擬觀測點的坐標和接收到射線的Pg波偏振方向。

對于震中一定范圍以外的觀測點，Pg波的到達往往會晚于Pn波。此時，Pg波的振幅信息會和其它震相的信息混雜在一起，這就對Pg波偏振方向的確定造成了困難。因此，在生成測試數據時，要注意將θ值限定在一定的范圍內。在本研究中，θ的取值范圍分為3檔，分別為≤ 30°、≤ 45°和≤ 60°，本文將對結果進行比較分析，以觀察不同的出射角范圍對定位精度的影響。同理，為了解觀測誤差和觀測點數量對于定位精度的影響，本文也選取了相應的參數來進行測試。為了計算上的方便，我們將三維空間簡化為二維，即去掉了一個水平維度，變成了如圖1所示的情況。這樣的簡化并不會對測試造成本質性的影響。

本文共虛擬了14萬組數據來進行測試，其詳細信息見表1。所有的數據按照不同的觀測點數量、出射角度范圍、誤差范圍等因素分入14個大組，每個大組包含1萬組數據。除觀測點數量以外，生成每一組數據所用到的角度和誤差都是在限定范圍內隨機生成的。顯然，當誤差水平為0，則意味著該組數據沒有引入相應的誤差。

3測試結果及比較

我們虛構了一個震中位置為0，震源深度為10 km的震源，并根據隨機生成的參數生成了觀測數據，最終反演得到了14萬組定位結果。對這些結果進行比較，以觀察不同的條件對于反演結果的影響。圖2～4直觀的表現了定位結果的差異。

圖2展示的是在無誤差情況下，觀測點的數量和分布對于定位結果的影響。以圖4可以看出，在無誤差情況下，所有的結果都很接近虛擬的震源位置，尤其是震中的誤差，被控制在數米之內；而深度的誤差也不超過05 km。觀測點的數量和分布都會對震中的反演造成影響，不同的是，對于深度的確定，觀測點的分布有著明顯的影響，但觀測點數量的影響卻幾乎可以忽略。綜合起來看，觀測點數量越多，分布的越密集，定位的效果就越好；觀測點的分布對于定位結果的影響比觀測點數量更加顯著。需要注意的是，圖2b中，隨著觀測點的分布越來越稀疏，定位結果中深度的概率峰值很明顯的變小。形象的說，觀測點的分布越稀疏，得到的震源深度越淺。同樣，圖2d顯示，隨著觀測點數量的增加，概率密度峰值有輕微的右移，即深度的定位結果輕微的增大。

圖3展示的是不同的誤差對于定位結果的影響。在引入誤差后，反演結果的概率密度曲線與圖2相比，明顯變得平緩，定位結果的誤差水平也顯著地增加。但是概率密度峰值對應的震中位置和深度并沒有偏離虛擬震源的位置。絕大多數的震中位置的誤差在50 m以內，深度的誤差水平較高，但是也集中在2 km以內。綜合起來看，Pg波偏振方向的誤差和觀測點坐標的誤差對于定位結果有很明顯的影響，并且都滿足誤差越小，定位越準確的規律。相比較而言，觀測點坐標誤差的影響要更大一些，但并不顯著。同時，隨著觀測點坐標誤差的增大，深度概率密度峰值也稍稍左移。并且，坐標誤差的增大會導致定位深度的減小。在反演結果中，出現了極少數遠遠偏離虛擬震源位置的異常結果（主要是深度值異常，圖中未顯示），在后面的章節中，將對這種情況進行討論。從反演結果的誤差水平上來看，各種誤差對于定位結果的影響要遠遠大于觀測點的分布和數量對結果的影響。

圖4顯示的是在一定誤差水平下（Pg波偏振方向的正弦值誤差≤ 1%，觀測點的坐標誤差≤ 1%）的定位結果。同圖2相比，定位的精度極大的降低了，同時深度的概率密度峰值隨觀測點分布和數量的變化趨勢都被一定程度的放大了。只是因為誤差的影響太大，把這個變化給掩蓋住了。同樣的，反演結果中也出現了極少的遠遠偏離虛擬震源的異常點。

4震例檢測

本文選取了2014年8月12日05：[KG-*2]55：[KG-*2]2495發生在云南省大理州賓川縣的28級地震對本方法進行檢測。在震中20 km范圍內，共有15個地震記錄儀獲取了本次地震的記錄。定位的結果見表2。

通過表2可看出，本方法的定位結果與云南地震臺網的定位結果很接近，震中位置相差小于2 km，震源深度相差也小于2 km。

5討論與結論

通過對射線追蹤方法的測試，我們認為該方法確實可以確定震源的空間位置。但是，由于需要獲得比較“純凈的”Pg波偏振信息，因此，該方法不適用于遠震定位。

該方法測試時，發現誤差對于定位精確度的影響非常大，遠遠超過了觀測點的分布和數量對于定位精確度的影響。總的來說，降低誤差水平，提高觀測點的數量和密度都能夠提升定位結果的精確度；反之亦然。另外，觀測點的密度降低和坐標誤差的增大會使深度的定位結果變小；而觀測點數量的增加會導致深度的結果變大。

和深度不同，震中的最大概然值并沒有隨著其它條件的變化而變化，造成兩者間區別的原因可能是觀測點在不同方向上的分布不同。由于出射角的生成是隨機的，因此，觀測點在水平方向上的分布應該是比較均勻的，不會產生大部分的觀測點都位于震中某一邊的情況；而在豎直方向上，所有觀測點的值都為0，全部位于震源的上方。在有誤差存在的反演結果中，出現了極少數的遠遠偏離震源的異常值。對于這些異常的產生，我們推測有2種可能的原因：隨機誤差的不均勻；反演程序自身的缺陷。對于震中距不同的觀測點來說，其在反演過程中的權重是不一樣的。這一點可以通過圖2a、b和圖4a、b看出。因此當某個權重較大的觀測點獲得了較大的誤差后，就會對定位過程產生嚴重的影響。另一方面，如果收斂條件的設定不合適，也可能會導致反演程序給出局部極小值點。

從測試中可以看出，震中定位的精確度是很高的。因此，將射線追蹤法（或者其它定位方法）得到的震中作為先驗條件，再對震源深度進行單獨定位，也許會得到更加可靠的結果。另外，我們注意到，在反演的過程中，我們將地表和震源中間的介質分層，并得到了相應的波速比值r（r1，r2，…，rn）。利用這組數值，我們可以通過線性反演獲得“發震時刻”和各層介質的“波速”。但是，由于我們并不清楚這組波速比值是否對應著真實的速度結構，因此，我們也不知道通過該方法得到的“發震時刻”和“波速”是否是正確的。

通過一例真實地震的檢驗，發現本定位方法的結果和其他研究者的結果較為接近，因此，我們認為射線追蹤方法可以用于真實地震的定位。

今后，我們將從5個方面展開工作：改進程序和算法，降低異常結果的產生；提高定位結果，尤其是深度定位結果的精確度；發震時刻的定位；復雜介質條件下，射線追蹤法的應用；用更多真實的地震數據來檢測該方法，并最終將其應用于實際的地震定位中。

[HTK]在本文的寫作和修改過程中，楊潤海高工提供了地震資料和指導，顏其中研究員對本文提出了建議，我的碩士導師倪四道教授和崔建文研究員也給予了指導和幫助，在此一并致謝。同時，也感謝審稿專家的修改意見和編輯的辛勤工作。

參考文獻：

李聰，雷建設2014滇西南地區地殼速度結構的區域波形反演[J].科學通報，59（34）：3398-3415

羅艷，曾祥芳，倪四道2013震源深度測定方法研究進展[J].地球物理學進展，28（5）：2309-2321

田玥，陳曉非2002地震定位研究綜述[J].地球物理學進展，17（1）：147-155

朱元清，趙仲和1997提高地震定位精度新方法的研究[J].地震地磁觀測與研究，18（5）：59-67

CROSSON R S1976Crustal structure modeling of earthquake data，1，Simultaneous least squares estimation of hypocenter and velocity parameters[J].J Geophys Res，81（17）：3036-3046

DOUGLAS A1976Joint epicenter determination[J].Nature，215：45-48

ENGDAHL E R，HILST R V D，RAYMOND B1998Global teleseismic earthquake relocation with improved travel times and procedures for depth determination[J].Bull Seism Soc Am，88（3）：722-743

GEIGER L1912Probability method for the determination of earthquake epicenters from arrival time only[J].Bull St Louis Univ，（8）：60-71

KLEIN F W1978Hypocenter location program HYPOINVERSE Part I：Users guide to versions 1，2，3 and 4[R].U S Geol Sur Open-File Rept，78-694

LEE W H K，LAHR J C1975HYPO71：A computer program for determining hypocenter，magnitude，and first motion pattern of local earthquake[R].U S Geol Surv Open-File Rept，75-311

LIENERT B R，BERG E，FRAZER L N1986Hypocenter：An earthquake location method using centered，scaled，and adaptively damped last squares[J].Bull Seism Soc Am，76（3）：771-783

LOMNITZE C1977A fast epicenter location program[J].Bull Seism Soc Am，67（2）：425-431

NELSON G D， VIDALE J E1990Earthquake locations by 3-D finite-difference travel times[J].Bull Seism Soc Am，80（2）：395-410

PRUGGER A F，GENDZWILL D J1988Microearthquake location：A nonlinear approach that makes use of a simplex stepping procedure[J].Bull Seism Soc Am，78（2）：799-815

ROMNEY C1957Seismic waves from the Dixie Valley Fairview Peak earthquakes[J].Bull Seism Soc Am，47（4）：301-319

SPENCE W1980Relative epicenter determination using P-wave arrival-time differences[J].Bull Seism Soc Am，70（1）：171-183

TARANTOLA A，VALETTE B1982Inverse problem = quest for information[J].J Geophys R，50：159-170

THURBER C H1985Nonlinear earthquake location：theory and examples[J].Bull Seism Soc Am，75（3）：779-790

WALDHAUSER F，ELLSWORTH W L2000A double-difference earthquake location algorithm：method and application to the Northern Hayward Fault，California[J].Bull Seism Soc Am，90（6）：1353-1368

WARREN L，SHEARER P2005Using the effects of depth phases on P-wave spectra to determine earthquake depths[J].Bull Seism Soc Am，95（1）：173-184

ZHAO L S，HELMBERGER D V1994Source estimation from broadband regional seismograms[J].Bull Seism Soc Am，84（1）：91-104

ZHU L P，HELMBERGER D V1996Advancement in source estimation techniques using broadband regional seismograms[J].Bull Seism Soc Am，86（5）：1634-1641[ZK）][HJ][FL）][LM]

[STHZ][WT4HZ][JZ]Near Earthquake Locating Based on Ray Tracing Method

[WT5B1][STBZ][JZ]ZHANG Qian1，2，YIN Geng2，WANG Yushi3，LIN Guoliang2

[WT5"B1X][JZ]（1 University of Science and Technology of China，Hefei 230026，Anhui，China）[JZ]

（2 Earthquake Administration of Yunnan Province，Kunming 650224，Yunnan，China）[JZ]

（3 Institute of Geophysics，CEA，Beijing 100081，China）

[WT5HZ][JZ]Abstract[WTB1]

We demonstrate a ray tracing method for earthquake locating without prior velocity structure，which is based on ray theory As the seismic ray focus on the hypocenter，first we divide the layer into sub-layers，and compute the velocities of the sub-layers and the location of the hypocenter by grid-search；then we increase the number of the sub-layers step by step，and repeat the process before correcting the hypocenter location；at last，we get a more accurate result Simulative data have been used to test the method，the results show that，this method can be used for earthquake location in theory We use this method to locate a natural earthquake In contrast with the location result of Yunnan Seismic Network，we find the two results are close