一種基于插值技術高精度計算稀疏網格地震定位中震源軌跡的方法

趙愛華

中國地震局地球物理研究所, 北京 100081

0 引言

地震定位在防震減災中起著十分重要的作用，不僅關乎大震救災，而且是進行大震預測的重要基礎(張國民等，2001；陳運泰，2009).大震之前，異常的地震活動范圍往往很大，例如1975年的海城地震(朱鳳鳴等，1982)；大震之后，余震活動的范圍可達數百公里，例如2011年的日本MW9.0地震(Asano et al., 2011)；數千公里之外的事件也可能觸發本地的中強震，例如2008年9月11月印度尼西亞的MW6.6地震觸發了約4750 km之外的日本MW6.9地震(Lin, 2012)，因此，為準確判斷震情(盡管很難)，常常需要對大尺度范圍的事件進行定位.高精度的地震定位不僅是分析地震孕育環境和發生機理的基礎，而且對區域、全球深部結構的研究(曾融生,1991；滕吉文等，2004；Rawlinson et al.，2010)及核爆監測(Bowers and Selby,2009)都具有重要意義.和過去相比，現代地震臺網的定位能力有了很大改善，但定位精度，特別是區域和全球尺度地震定位中震源的深度精度仍有待提高(Schweitzer, 2006; Bondár and Storchak, 2011).

區域和全球尺度的地震定位(例如Xie et al.,1996；Engdahl et al., 1998)和煤炭、油氣開采中的地震定位(例如Boltz et al., 2014; Bao and Eaton, 2016; 常旭和王一博，2019)，在原理上是一樣的，即在地下空間中尋求一點使理論計算的和實際觀測的到時(或到時差)之差最小(Thurber,1986; Moser et al.,1992)，所不同的是前者的速度模型比后者在尺度上大得多.為減小理論到時誤差以提高定位精度，區域和全球的地震定位也越來越多地使用更接近實際的三維速度模型(陳棋福等，2001；Lin et al., 2007; De Kool and Kennett, 2014; Myers et al., 2015).由于三維速度模型的地震波走時難以表示成解析表達式，因此通常使用射線追蹤技術進行數值求解，而模型則以網格剖分成模型單元.區域和全球地震定位中的速度模型，尺度較大，不宜以較小網格剖分，否則剖分的單元過多、走時計算量太大.再者，速度模型主要來自地震層析成像結果，而目前區域和全球尺度的地震層析成像，分辨率還比較低(例如，Kohler et al.，2003; Lin et al., 2010; Xin et al., 2019; Simmons et al., 2019)，這樣，即使以較小網格剖分模型，也不會在根本上提高理論走時的計算精度.當模型以稀疏網格剖分時，由于剖分的單元尺寸較大，走時計算將變得較為復雜，例如，模型界面難以用模型單元的中心點連線近似表示，否則，會人為地產生較大的齒狀斷層；相鄰單元中心點之間的射線路徑用直線代替將可能導致較大誤差.因此，對于大多數的地震定位方法而言，例如，蓋革法(Geiger, 1912)、主事件法(Spence, 1980)、雙差法(Waldhauser and Ellsworth, 2000)、搜索法(Prugger and Gendzwill, 1988； Kennett and Sambridge,1992；Billings，1994；Shearer,1997；Sambridge and Kennett, 2001)以及震源參數與速度結構聯合反演法(Crosson, 1976; Zhang and Thurber, 2003)等，區域和全球地震定位的困難主要在于理論走時的計算，但對于地震定位交切法則不盡然.交切法基于幾何思想定位，將震源軌跡交會最集中的點作為震源位置，具有穩健直觀的優點(Garza et al., 1977,1979；傅淑芳和劉寶誠，1991；Pujol, 2004; 廉超等，2006)；但當速度結構較為復雜時，震源軌跡難以找到解析解(Zhao and Ding, 2007, 2009; 周建超和趙愛華，2012).因此，將傳統交切法應用于高精度的區域和全球地震定位，還需要解決大尺度復雜模型中震源軌跡的求解問題.

對于復雜速度模型中的震源軌跡，目前主要有兩種求解方法.第一種，以殘差(理論計算的與實際觀測的到時或到時差之差)絕對值較小的單元節點近似表示震源軌跡，單元節點通常取為模型單元的中心點(Zhou，1994；Font et al., 2004；Theunissen et al., 2012).這種方法，簡單易行，但當模型以稀疏網格剖分時，由于模型節點的間距較大，得到的震源軌跡有較大誤差.為減小誤差，在稀疏網格中初步確定震源位置后，再以細網格對初步確定的震源其周圍區域進行剖分，更精確地確定震源位置.這樣做，雖然在一定程度上解決了問題，但無疑增加了工作量，而且實現起來也較煩瑣.另外，以絕對殘差較小的點表示的震源軌跡不是理論上的曲面而是厚度不均的板帶，軌跡所交會的為一個區域而非一個點，即使速度模型和到時都準確.針對該問題，本文作者與合作者提出了另外一種求解震源軌跡的方法，其基本思想是：選取震源軌跡所經過模型單元的節點為軌跡參考點，以絕對殘差場中連接參考點的射線路徑表示震源軌跡(趙愛華等，2008，2015；趙愛華，2018).相比于第一種方法，第二種方法得到的震源軌跡更為精細，根據震源軌跡的交會情況可直觀地評估定位效果，便于對觀測到時進行優化組合、設計最佳的地震定位方案；不過，該種方法同樣存在當模型以稀疏網格剖分時所得到的震源軌跡誤差較大的問題.高精度的區域和全球地震定位，不僅范圍廣，而且臺站多、事件量大，要求震源軌跡的計算有較高的精度和效率.

為此，本文嘗試對基于射線追蹤技術的震源軌跡計算方法(趙愛華，2018)進行改進，使之滿足區域和全球事件定位的需要：當模型以稀疏網格剖分時，計算的震源軌跡也具有較高的精度，并且計算效率較高.雖然地震定位是一個三維問題，但為便于敘述，方法說明仍以二維模型為例.

1 震源軌跡方程

1.1 到時差約束的震源軌跡方程

(1)

由方程(1)可得到時差殘差(RDT)場FRDT：

FRDT(x;Rj1,Wj1;Rj2,Wj2)=[T(x;Rj1,Wj1)

(2)

其中，I(x)為空間單位場，其值處處為1.

1.2 到時約束的震源軌跡方程

(3)

(4)

到時殘差(RAT)場FRAT為

(5)

其中，I(x)為空間單位場.

值得指出的是，到時約束的震源軌跡不僅與作為約束條件的觀測到時有關，還與構建發震時間場的觀測到時有關.發震時間場不同，同一到時約束的震源軌跡也不同，甚至相差很大(趙愛華等，2015).

2 震源軌跡參考點的計算

對于速度結構復雜模型，震源軌跡難以準確計算，不便于進行誤差分析.為此，考慮一個均勻介質模型，模型大小為150 km×90 km，介質P波和S波的速度分別為6.0 km·s-1和3.6 km·s-1(趙愛華等，2015).為考察網格稀疏時的震源軌跡計算誤差，將模型剖分成較大單元，例如5 km×5 km的正方形，在其淺部布設5個地震臺站Ri(13.5+(i-1)×30.0 km,-2.5 km)(i=1,2,3,…,5).設在點x0(37.5 km,-31.5 km)處，時間t0=0 s時發生地震.圖1顯示了臺站R3處P波和S波到時差約束的震源軌跡(藍色實線)及其殘差場，圖中棕色點為軌跡所經過模型單元的中心點(為便于敘述，將其稱為震源軌跡節點).可以看出，在均勻各向同性介質中，同一臺站P波與S波到時差約束的震源軌跡為圓形；震源軌跡節點和軌跡的理論解，有明顯的偏差，這一點在圖1中局部放大的圖中顯示得更為清楚.可以推知，當模型單元較大時，選取震源軌跡節點為軌跡參考點、以連接軌跡參考點的線作為震源軌跡會有較大誤差.

易知，選取的參考點越接近實際震源軌跡，計算的震源軌跡越精確，最理想的情況是選取震源軌跡上的點即殘差場中的零值點為軌跡參考點.觀察圖1可知，殘差從震源軌跡一側到另一側，從正變為負或反之；沿軌跡法線方向，殘差變化最快.據此，可以假定：在震源軌跡兩側附近，殘差沿軌跡法線方向的變化可表征為簡單的函數.若使用震源軌跡兩側法線上鄰近點構建插值函數，則插值函數的零點應離實際震源軌跡較近.在計算震源軌跡之前，準確的軌跡法線方向是難以確定的.不過，根據震源軌跡在殘差場中的特點，和軌跡法線相近的方向，例如法線點對方向(趙愛華，2018)卻是可以確定的.法線點對由位于震源軌跡兩側的鄰點組成，其方向和軌跡法線方向相差不超過45°.這樣，基于插值技術計算法線點對間殘差為零的點作為震源軌跡參考點，也應具有較高精度.根據這種認識，本文提出了如下計算震源軌跡參考點的算法：

圖1 均勻介質中事件以同一臺站P波與S波到時差約束的震源軌跡及其殘差場棕色點為震源軌跡所經過模型單元的中心點；紅色十字為使用插值方法計算的震源軌跡參考點；藍色實線是震源軌跡的理論解,綠色虛線為±1.0 s的殘差等值線.Fig.1 A hypocentral locus constrained with the travel time difference between the P - and S- wave at the same station and its residual field of an event in a homogeneous mediumBrown dots are nodes of the mode elements traversed by the hypocentral locus; red crosses indicate the focal locus reference points by linear interpolation. Blue lines show the theoretical solutions of the hypocentral locus; the green lines are the residual contours of ±1.0 s.

(1)組成點對

對每個模型節點，考慮其上下、左右4個鄰點(三維，再加前后2個鄰點)； 若存在和節點殘差正負極性不同的鄰點，則將這樣的鄰點和節點組成點對.

(2)確定法線點對

選取同一模型節點其全部點對中絕對梯度最大者為法線點對.法線點對中絕對殘差較小的點更靠近震源軌跡，震源軌跡經過其所在模型單元，即該節點為震源軌跡節點.

(3)插值求解

使用法線點對及其兩側延長線上的模型節點構建殘差的插值函數，選取插值函數中殘差為零的點作為軌跡參考點.

需要指出的是，盡管模型節點的坐標是二維或三維的，但殘差在法線點對方向上的變化是一維的，因此構建插值函數較為簡單.對于法線點對A(xA,zA,RA)和B(xB,zB,RB)(其中，R為殘差值)，若殘差在點對AB間線性變化，則AB間殘差為零的點P,其坐標為

(6)

若殘差在點對AB間非線性變化，則可以用二分法求取殘差為零的點.圖1中，紅色十字是使用線性插值方法計算的震源軌跡參考點.可以看出，計算的軌跡參考點和軌跡的理論解(藍色實線)幾乎完全重合，具有較高精度.

3 射線路徑的計算

震源軌跡的參考點是離散的點，雖然可直接用于地震定位，但以其表示的震源軌跡不太形象和細密；若用線將離散的軌跡參考點連接起來，則有利于更直觀地考察震源軌跡在空間的展布.由式(2)和(5)可知，震源軌跡在絕對殘差場中具有近于零的低值.這樣，若將絕對殘差看作地震波速度的倒數(即地震波慢度)，在絕對殘差最小處激發地震波，則地震波將沿震源軌跡快速傳播，而在軌跡兩側傳播較慢.因此，震源軌跡參考點可在絕對殘差場中以射線路徑連接(趙愛華等，2008).絕對殘差場，非均勻性較強(參見圖1)，因此，連接軌跡參考點的射線路徑，應選用穩健性好的射線追蹤方法計算，例如最小走時樹法(Nakanishi and Yamaguchi, 1986).最小走時樹射線追蹤方法以構建最小走時樹的方式計算地震波傳播路徑和走時(Dijkstra, 1959；Cao and Greenhalgh, 1993)，沒有陰影區問題，計算的走時全局最小，適于復雜模型的計算；經過不斷地改進和完善(Moser, 1991;王輝和常旭, 2000; Zhao et al., 2004; 趙愛華和徐濤, 2012)，計算效率也較高.當震源軌跡由多段組成時，僅將絕對殘差場中全局最小的那個節點設置為射線路徑初始點，計算結果將包含虛假震源軌跡；若將射線路徑計算區域限制在絕對殘差較小的區域，計算區域依震源軌跡的段數組成自適應地劃分為若干連通區域，選取每個連通區域中絕對殘差最小的點為震源點，則無虛假軌跡問題(趙愛華等，2015).鑒于根據殘差極性可選取全部或幾乎全部的震源軌跡節點(趙愛華，2018)，本文將選出的震源軌跡節點所在的模型單元作為射線路徑的計算區域.常規的最小走時樹法僅計算模型節點之間的射線路徑，而不能處理非模型節點的點.注意到震源軌跡節點和計算的軌跡參考點存在對應關系，因而可采取先計算連接震源軌跡節點的射線路徑、再以計算的軌跡參考點代替射線路徑中震源軌跡節點的辦法來計算震源軌跡.

4 模型算例

4.1 速度模型

為驗證方法的適用性，將改進的震源軌跡計算方法應用于一個接近實際的速度模型.該模型以云南地區為背景，修改自遮放—賓川速度剖面(白志明和王椿鏞，2004).模型速度結構如圖2所示，S波與P波速度之比vS:vP為常數1∶1.73.可以看出，速度的縱向和橫向非均勻性都較強；地殼與地幔之間的界面，即莫霍(Moho)面，深度變化顯著.參照文獻(趙愛華等，2015)做法，將模型剖分成邊長為1 km的正方形單元，模型單元內速度設為常數；布設于地表的9個地震臺站Ri(i=1, 2, 3…9)，深度為0.5 km,臺站間距為20或30 km；將地震事件設置在點x0(130.5 km, -27.5 km)處，發震時間為t0=0 s.震源至地震臺站6種震相(Pg、Sg、PmP、SmS、Pn和Sn)的到時，使用子波傳播區域(控制計算精度)動態調整、效率較高的最小走時樹射線追蹤方法(趙愛華等，2003)計算.PmP/Pn波的射線路徑，如圖2中灰線所示.

對于圖2中模型的事件，作者及其合作者曾計算了Pg波和Sg波其到時差和到時約束的震源軌跡(趙愛華等，2015).已有的研究表明，使用莫霍面反射波(PmP和SmS)或折射波(Pn和Sn)有利于震源深度的確定(白玲等，2003；Husebye et al, 2013; Wagner et al., 2013).關于反射波和折射波的震源軌跡，目前研究得還較少.因此，本文將重點考察這類震相到時或到時差約束的震源軌跡.

圖2 以云南地區為背景的速度模型(修改自白志明和王椿鏞,2004)紅色十字為震源位置；黑色倒三角形為地震臺站；灰線為地震反射/折射波射線路徑.Fig.2 A velocity model set in Yunnan area (modified from Bai and Wang, 2004)The red letter “+” indicates the hypocenter; black inverse triangles denote the seismic stations; gray lines show the ray paths of reflected or refracted seismic waves.

4.2 到時差約束的震源軌跡

全部9個臺站兩兩組合，對于某個震相，總計可產生36個到時差.由Pg波到時差約束的震源軌跡可知，其中有些到時差約束的震源軌跡在方位分布上是非常接近的.因此，計算所有臺站組合的震源軌跡是不必要的.為便于和Pg波對比，本文計算了相同臺站組合、PmP/Pn波到時差約束的震源軌跡，計算結果如圖3a所示.圖例線上的數字對“i-j”表示震源軌跡所對應的臺站組合為Ri和Rj.可以看出，震源軌跡在震源附近有良好的方位分布，能很好地約束震源位置，這和Pg波到時差約束的震源軌跡(趙愛華等，2015)相似；但和后者僅交會于震源點不同，在莫霍面之下、約在震源點的對稱位置還有較集中的交會.震源軌跡在莫霍面發生急劇方向變化可能源于界面兩側(上：反射波/折射波；下：直達波)的走時場具有不同的分布特征.

相同臺站PmP/Pn波與SmS/Sn波到時差約束的震源軌跡顯示在圖3b中，軌跡旁邊的數字“i”表示軌跡所對應的臺站為Ri.可以看出，軌跡在震源附近有很好的方位分布；在震源處，近臺(例如R3)的軌跡以水平分布為主，對震源深度有較好的約束，而遠臺(例如R1)的軌跡則以垂直分布為主，可較好地約束震源的水平位置.值得注意的是，近臺(例如R3和R4)的震源軌跡為在地下閉合的曲線.可以證明，若莫霍面水平、地殼為均勻各向同性介質，則同一臺站(位于地表)反射縱、橫波到時差約束的震源軌跡(反射界面之上部分)為圓形，圓心在臺站之下2倍地殼厚度處.因此，對于近臺，震源軌跡可能不會出露于地表，這點和直達縱、橫波到時差約束的震源軌跡不同.

圖3 圖2中事件以到時差約束的震源軌跡(a) 不同臺站PmP波到時差； (b) 相同臺站PmP波與SmS波到時差; (c) 相同臺站PmP波與Pg波到時差.不同顏色和線型的曲線表示不同震源軌跡，其旁邊的數字“i”或數字對“i-j”表示震源軌跡以臺站Ri或臺站對Ri-Rj的到時差約束.Fig.3 Hypocentral loci constrained with arrival time differences for the event in Fig.2(a) Traveltime differences between PmP waves at different stations; (b) Traveltime differences between PmP and SmS waves at the same stations; (c) Those between PmP and Pg waves at the same stations. Curves in different colors and styles represent different hypocentral loci and their sideward figures “i” or figure pairs “i-j” imply that they are constrained with the arrival time differences from stations Ri or station pairs Ri-Rj.

相同臺站直達縱波(Pg)和反射/折射波(PmP/Pn)到時差約束的震源軌跡如圖3c所示，軌跡旁邊的數字“i”表示軌跡所對應的臺站為Ri.可以看出，無論是近臺(例如R3、R4和R5)還是遠臺(例如R1、R2和R9)，軌跡在震源處均以水平分布為主，似乎對震源深度都有很好的約束.根據以前的研究(趙愛華等，2015)可知，震源軌跡約束震源位置的能力，不僅與軌跡在震源附近的方位分布有關，還與軌跡的穩定性有關.軌跡越穩定，受擾動影響越小，對震源位置的約束越強.在絕對殘差場中，震源軌跡越穩定，其所在的“峽谷”(將絕對殘差看作高程)越窄陡；反之亦然.圖4顯示了圖3c中近臺R3和遠臺R9所對應的震源軌跡及其絕對殘差場.為清晰顯示震源軌跡，對絕對殘差做了取對數處理.易見，相同臺站PmP/Pn波與Pg波到時差約束的震源軌跡其所在的“峽谷”，近臺的為陡窄的“V”形，而遠臺的為半側寬緩的反“L”形.這說明，前者軌跡較為穩定，對震源深度有較強約束；后者軌跡穩定性較弱，較小的擾動就可使其形態或位置發生顯著變化，雖然在震源處以水平分布為主，但對震源深度的實際約束則較弱.這與地震定位中“近臺約束強、遠臺約束弱”的經驗相符.

圖4 圖3c中近臺(上)和遠臺(下)震源軌跡的殘差場Fig.4 Residual fields of hypocentral loci for near (up) and far (down) stations in Fig.3c

4.3 到時約束的震源軌跡

全部臺站PmP/Pn波到時約束的震源軌跡如圖5a所示.和到時差約束的震源軌跡(圖3a)相比，到時約束的震源軌跡在震源附近的方位分布略窄，不過有兩條軌跡在震源點具有較大的方位夾角，因此仍能較好地約束震源位置.當以全部臺站PmP/Pn波和Pg波到時構建發震場時，PmP/Pn波(藍線)和Pg波(綠線)到時約束的震源軌跡顯示在5b中.可以看出，增加Pg波到時后，震源軌跡與之前的相比，方位分布大大改善.較為全面的方位覆蓋使震源位置得到有力約束；更多條震源軌跡約束震源，有利于減少隨機干擾對定位結果的影響.

在地震定位實踐中，事件被近臺和遠臺都記錄到是較理想的情況.微小地震，激發的能量較弱，往往難以被遠臺記錄到而僅能用近臺定位.研究表明，近臺在約束震源水平位置方面起著關鍵作用(Bondár et al, 2004).大震之后在震源附近沿發震構造布設流動臺陣其目的即是用近臺對余震活動進行更好地監測.圖5c和5d顯示了近臺(R2-R6)PmP/Pn波及其與Pg波到時約束的震源軌跡.可以看出，當僅有近臺時，PmP/Pn波單獨約束的震源軌跡在震源處的分布方位上更加不均，除了有一條近水平分布外，其余則都集中地分布于垂直方向很窄的方位范圍內；PmP/Pn波和Pg波共同約束的震源軌跡，則和全部臺站的情形類似，在震源處的方位分布全面且較均勻，對震源位置有很好的約束.

圖5 圖2中事件不同觀測情況以到時約束的震源軌跡(a)和(b) 所有臺； (c)和(d) 近臺； (e)和(f)遠臺； (g)和(h) 右側臺； (i)和(j) 擾動到時. 左邊圖僅使用PmP波到時；右邊圖使用PmP和Pg波到時，其中藍線和綠線分別對應PmP和Pg波. 紅色“+”為震源位置.Fig.5 Arrival time constrained hypocentral loci of the event in Fig.2 for different observations(a) and (b) complete observation; (c) and (d) near observation; (e) and (f) far observation; (g) and (h) right-side observation; (i) and (j) disturbed arrivaltimes. Hypocentral loci in the left subfigures involve only PmP-wave data; those in the right ones involve both PmP- and Pg-waves, where they are shown in blue and green respectively for PmP and Pg waves. The red “+” indicate the true hypocentral position.

在地震臺網稀疏地區，很多時候僅有遠臺可用于定位.為模擬缺少近臺的定位情形，僅保留距震源較遠的5個臺站：R1、R2、R7—R9.PmP/Pn波及其與Pg波到時約束的震源軌跡分別如圖5(e，f)所示.PmP/Pn波單獨約束的震源軌跡在震源處均近垂直分布，可很好地約束震中位置，但難以約束震源深度；增加Pg波到時后，軌跡在震源處的分布方位：范圍略有加寬，但仍以垂直方向為主.這表明，缺少近臺，即使聯合使用PmP/Pn波和Pg波也難以很好地約束震源深度.

在區域臺網的地震定位中，事件常常位于臺網之外.為模擬該種情形，僅保留右側5個臺站R5-R9.PmP/Pn波到時約束的震源軌跡如圖5g所示.震源軌跡在震源處以和水平方向成約135°(以逆時針方向為正)的方位成束分布，對震源位置約束較弱.增加Pg波到時信息后，震源軌跡(圖5h)在震源處的方位分布得到顯著改善，對震源位置特別是震源深度的約束顯著加強.

震相到時誤差是影響定位精度的關鍵因素之一.采用互相關技術(例如Akram and Eaton, 2016)或人工智能技術(例如Ross et al., 2018;Wang et al., 2019)可大大提高到時拾取精度，但若環境噪聲較大，準確地拾取到時仍較困難.為模擬到時拾取存在誤差情況，使用FORTRAN編譯器中的隨機函數對PmP/Pn波和Pg波的理論到時加以擾動，兩種波的擾動幅度均為±0.5 s，擾動均方值分別為0.321 s和0.356 s.PmP/Pn波單獨約束、PmP/Pn波和Pg波聯合約束的震源軌跡分別如圖5(i,j)所示.可以看出，當觀測到時有誤差時，震源軌跡不再交會于震源點，而是交會成一個小的區域.圖5i中閉合的軌跡(箭頭所指)明顯偏離震源點，而且形態發生了較大變化(參見圖5a中箭頭所指軌跡), 表明其穩定性可能不是很強.和單獨震相到時約束的震源軌跡相比，兩種震相到時聯合約束的震源軌跡其交會最密集的點更接近實際震源位置.

和Pg波到時約束的震源相比(參見趙愛華等，2015)，PmP/Pn波到時約束的震源軌跡在震源處的方位分布是類似的；不同的是，后者除了在震源處有較集中的交會外，在莫霍面之下也有密集交會的點.軌跡在莫霍面發生急劇方向變化可能主要源于界面上下的地震波走時場具有不同分布特征.

5 討論

5.1 計算誤差

由計算過程可知，震源軌跡的計算精度決定于軌跡參考點的精度.軌跡參考點越接近實際的震源軌跡，以連接參考點之射線路徑表示的震源軌跡越精確.對于復雜介質模型(例如圖2中模型)，難以找到震源軌跡的理論解.為此，本文以圖1中均勻介質模型中事件的震源軌跡(其解析解為標準的圓形)來例，分析模型單元尺寸及插值方法對震源軌跡參考點計算精度的影響.模型單元邊長d=1、3、5、7、9 km時，使用不同方法計算軌跡參考點的結果如表1所示.可以看出：對于選用震源軌跡所經過模型單元中心點作為軌跡參考點即震源軌跡節點法，軌跡參考點的誤差和單元邊長d基本成正比，最大值和平均值分別約為d/2和d/4；插值法軌跡參考點的精度雖然總體上也是隨著單元邊長d的增大而降低，但比震源軌跡節點法的高幾十(線性插值)到幾百倍(拋物插值和三次樣條插值)；當模型單元較大時，插值函數次數越高，計算的軌跡參考點誤差越小：以d=9 km為例，線性(1次)插值法、拋物插值(2次)法和三次樣條(3次)法軌跡參考點的誤差最大值分別為0.110、0.027和0.015 km；但這并不意味著使用更高次的插值函數計算軌跡參考點的效果會更好.高次插值函數存在震蕩現象，求取函數法線點對間的零值點(可能有多個)較為復雜，并且可能由于插值函數對震源軌跡附近區間殘差的變化表征得更差而使得其零值點離實際震源軌跡更遠.

在表1中，三次樣條差值法d=3 km時的軌跡參考點誤差為0.004 km,比其d=7 km時0.003 km的誤差還要大，似乎不合理.為解釋該現象，考察了圖1中殘差沿震源軌跡法線方向的變化，如圖6所示.軌跡法線為水平，顯示的區間從A(103.5 km,-31.5 km)至B(112.5,-31.5 km).圖中紅色虛線為連接點A和B的直線.可以看出，殘差在震源軌跡附近沿法線方向近似為線性變化,但有波動.這種波動有可能使得以震源軌跡兩側較近的點(對應較小模型單元)比以軌跡兩側較遠的點(對應較大模型單元)構建三次樣條插值函數、函數的零值點離實際軌跡更遠.

表1 不同方法計算圖1中震源軌跡參考點的精度比較Table 1 Accuracy comparison of different methods for calculating the hypocentral locus in Fig.1

圖6 圖1中到時差殘差沿震源軌跡法線A(103.5,-31.5)B(112.5,-31.5)方向的變化Fig.6 Variation of traveltime difference residual along the normal A(103.5,-31.5)B(112.5,-31.5) of the hypocentral locus in Fig.1

5.2 計算效率

由表1可以看出，粗網格時使用插值法計算軌跡參考點可達到細網格時震源軌跡節點法的精度.以最簡單的線性插值為例，d=9 km時其軌跡參考點的精度比震源軌跡節點法d=1 km時的精度還高.對于給定模型，d=1 km時的模型單元數是d=9 km時的81倍(二維)或729倍(三維).在震源軌跡計算中，地震波走時場的計算最為費時，計算時間和模型單元數成正比.這意味著，使用新的軌跡參考點計算方法可使用于計算地震波走時場的時間大大減少，盡管后者的計算與地震事件無關，對于穩定區域僅需1次，但對于尺度較大的三維模型或當地震臺站較多時，這種減少的效果還是非常顯著的.

在選取震源軌跡節點的基礎上，基于插值技術求取法線點對間殘差為零的點作為軌跡參考點需要花費額外的計算時間.表2列出了不同插值法計算軌跡參考點所額外花費的CPU時間，CPU時間使用編程語言FORTRAN中的函數timef計量.需要指出的是，程序連續重復運行多次，每次得到的CPU時間并不完全一致，表中所列是出現次數最多的結果.盡管數據不是很精確,但仍可以看出，使用插值法計算軌跡參考點額外花費的時間不是很多，以d=1 km時為例，線性、拋物、三次樣條三種插值法計算軌跡參考點所額外花費的時間分別為0.007 s、0.039 s和0.045 s，和常規方法(以震源軌跡節點為參考點)計算震源軌跡的時間(0.117 s)相比，分別增加5.98%、33.33%和38.46%，和地震波走時場的計算時間(單次為3.09 s)相比，增加的比例則更低.隨著單元尺寸的增大，由于震源軌跡節點個數減少，使用插值法計算軌跡參考點所額外花費的時間越來越少.

表2 震源軌跡參考點計算對計算效率的影響Table 2 Effect of the calculation of FLRPs on the efficiency

在射線路徑計算部分，需要將路徑中模型節點替換為計算的軌跡參考點.由于操作簡單，而且當使用稀疏網格時，需要替換的節點較少，因此這部分所花費的時間微乎其微，對軌跡的計算效率影響不大.

5.3 算法的實用性

和理論實驗相比，實際的地震定位有更多的不確定性因素，例如地下的速度結構難以準確獲知、拾取的震相到時存在或多或少誤差等.計算真實事件的震源軌跡可較好地檢驗本文算法的實用性.實際的地震存在于三維空間.對于三維速度模型，根據殘差正負極性可選取全部或幾乎全部的震源軌跡法線點對(趙愛華，2018)；法線點對雖然在三維空間分布，但其上的殘差變化是一維的，這樣，使用一元的插值方法仍可高精度地計算出三維殘差場中的零值點作為震源軌跡參考點.因此，本文所提出的震源軌跡計算方法在理論上是適用于三維速度模型的.盡管如此，為便于展示和分析，本文仍以二維模型中真實事件的計算為例，說明算法的實用性.

在二維空間中考察震源軌跡的分布，震中和臺站應位于同條直線.文獻(趙愛華等，2015)提供了較為符合這一要求的實例.地震發生在華北地區，北京數字遙測地震臺網(北京臺網)的測定結果為：震中位于117.42°E、39.61°N，震源深13 km，發震時間為2012年8月26日7時13分34.3秒.該地震具有良好的觀測方位覆蓋，并且有近臺(CAD臺，震中距僅3km)控制，定位到時殘差較小，可以認為定位結果可靠.震中和三個臺站(從南至北依次為EWZ臺、CAD臺和XAZ臺)基本位于同條直線(見趙愛華等，2015中圖6).根據三個臺站的Pg波和Sg波走時，參照華北東部地殼模型(滕吉文等，1979)，趙愛華等(2015)構建了一個水平層狀均勻的速度模型，模型由三層介質組成、尺度為60 km×35 km.該模型對于北京臺網確定的震源位置和震相到時，EWZ臺、CAD臺和XAZ臺的Pg波和Sg波走時殘差(理論走時與觀測走時之差，單位為s)分別為：(-0.08，0.40)、(-0.06，-0.35)和(0.05，0.05).他們使用0.1 km×0.1 km的網格計算了三個臺站Pg波和Sg波的走時場.和地震波走時場相比，震源軌跡的計算需要更多內存；地震波走時場的空間采樣率對計算軌跡的完整性與精細性影響不大.為此，當計算震源軌跡的程序可能因內存超限而不能運行時，他們對得到的地震波走時場進行了粗網格(0.2 km×0.2 km)重采樣以完成計算.使用其粗網格采樣的走時場和北京臺網確定的震相到時，本文計算了真實事件的震源軌跡，計算結果如圖7所示，對應不同臺站或臺站組合的震源軌跡以不同彩色線表示.可以看出：不同臺站Pg波到時差約束的震源軌跡(圖7a)交匯得較為集中，交匯點和震源位置僅略有偏差；同一臺站Pg波與Sg波到時差約束的震源軌跡(圖7b)交匯成一個區域，走時殘差越小的到時其約束的震源軌跡離震源點越近；Pg波到時約束的震源軌跡(圖7c)不僅交會得較集中，而且對震源位置有較好的方位覆蓋，其中臺站CAD所對應的軌跡(綠色線)在震源處近水平分布，對震源深度有較強約束；Pg波與Sg波到時約束的震源軌跡(7d, 彩色實線對應Pg波，彩色虛線對應Sg波)和圖7b中的軌跡類似，交匯得不太集中，但交匯區包含震源位置，可能意味著使用較多的震源軌跡進行定位有助于減少隨機擾動對定位結果的影響.Pg波單獨約束的震源軌跡(圖7a和7c)比Pg波與Sg波共同約束的震源軌跡(圖7b和7d)交匯得更為集中，這可能與Pg波走時殘差小、Sg波走時殘差較大有關.當震源軌跡交匯成區域時如何確定震源位置，是使用震源軌跡進行地震定位需要解決的另外一個關鍵問題.

圖7 華北地區2012年8月26日ML4.2地震的震源軌跡(a) P波到時差約束； (b) P和S波到時差約束； (c) P波到時約束； (d) P波和S波到時約束. 黑色線為趙愛華等(2015)計算結果， 彩色線為本文計算結果(詳情見正文). 倒三角形為地震臺站, 紅色十字為北京臺網所定震源位置在臺站連線上的投影.Fig.7 Hypocentral loci of the 26 August 2012 ML4.2 earthquake in North China(a) Constraint of arrival time differences between P waves; (b) Constraint of those between P and S waves; (c) Constraint of P-wave arrival times; (d) Constraint of P- and S-wave ones. Hypocentral loci in black come from Zhao et al. (2015) and those in colors (red, green and blue) are results in this study (see text for details). The inverse triangles represent the seismic stations and the red crosses indicate the projection of the hypocenter determined by Beijing Seismic Network to the station line.

對比文獻(趙愛華等，2015)的計算結果(圖7中黑色線)，可以發現本文計算的震源軌跡更為完整，如圖7中水平箭頭所指示部分.軌跡完整性的差異與軌跡參考點的選取或計算方法不同有關.趙愛華等(2015)在絕對殘差場中使用“削皮”法選取軌跡參考點，選取的參考點可能不都是震源軌跡節點，為此去掉較短路徑使計算軌跡精細的修飾會損害軌跡的完整性.本文計算的軌跡參考點與選取的軌跡節點一一對應，在殘差場中基于殘差正負極性選取軌跡節點，選出的點不但全是軌跡節點而且包含全部或幾乎全部的震源軌跡節點，因而計算的震源軌跡相當完整.

圖7中兩種計算結果另外一個顯著不同之處是早先計算的震源軌跡存在間斷，如垂直箭頭所指示部分.為分析計算軌跡中斷原因，圖8顯示了中斷最嚴重的震源軌跡(圖7d中對應CAD臺Pg波到時、以綠色實線表示的軌跡)的走時殘差場.可以看出，在殘差變化緩慢的區域，以線性插值計算的軌跡參考點(棕色點)和繪圖軟件以克里克插值法(空間采樣率亦為0.2 km×0.2 km)得到的0.0 s殘差等值線(綠色線)符合得相當好; 在殘差變化較劇烈的區域，兩者不太相符(見圖8a).將圖8a中計算軌跡左邊中斷附近的殘差場(紫色箭頭所指區域)放大(圖8b)，易見計算的軌跡參考點位于正殘差節點(紅色菱形)和負殘差節點(藍色菱形)之間，并且靠近殘差絕對值較小一側的節點，這說明以線性插值方法計算的殘差零值點較為準確.軌跡參考點不間斷的分布表明相應區段的震源軌跡為連續分布.計算軌跡中斷部分，總有某側、有時甚至雙側的相鄰節點絕對殘差值大于0.2 s(見圖8b中以空心菱形表示的節點).這樣，按照文獻(趙愛華等，2015)的方法，為保證計算軌跡的完整性，軌跡參考點選取區域的閾值RFLRP和射線路徑計算區域的閾值RFL均要大于0.2 s.若RFL>0.2 s，則射線路徑的計算區域將包括圖8a中部兩段軌跡之間的區域，使得單個連通區域內包含兩段軌跡、產生虛假軌跡(連接兩段軌跡之間的射線路徑).為防止產生虛假軌跡，趙愛華等(2015)將RFLRP和RFL取為較小的值(0.05 s)，使計算的震源軌跡產生中斷.本文算法可在整個模型區域中選取震源軌跡節點，只要將RFLRP取得足夠大(例如0.5 s)，就可選出全部或幾乎全部的軌跡節點，將選出的軌跡節點所在單元設置為射線路徑的計算區域，使得每個連通區域僅包含一段震源軌跡，即使兩段軌跡之間區域的絕對殘差很小也不會產生虛假軌跡，從而可很好地同時兼顧計算軌跡的完整性與精細性.

圖8 圖7d中綠色實線震源軌跡(對應CAD臺站P波走時)的殘差場(a)及其局部放大圖(b)棕色點為使用線性插值計算的震源軌跡參考點;菱形為震源軌跡兩側相鄰的模型節點(間距0.2 km),其中藍色的表示殘差小于0,紅色的表示殘差不小于0, 空心的表示殘差絕對值大于0.2 s,實心的表示殘差絕對值不大于0.2 s. 綠色線為到時殘差等值線(單位為秒).Fig.8 Entire (a) and enlarged partial (b) residual fields of the hypocentral locus (responding to the P- wave arrival time at station CAD) shown as solid green lines in Fig.7dBrown dots are focal locus reference points by linear interpolation. Diamonds indicate model nodes (spaced 0.2 km apart) adjacent to the hypocentral locus, where the blue are negative and the red are nonnegative residual values; the empty are bigger and the solid are not bigger than 0.2 s of the absolute residual. Green lines are contours of arrival time residual (in second).

真實事件的計算結果表明：本文算法具有良好的穩健性，即使震源軌跡所在區域殘差變化劇烈，也能得到可靠結果；對震源軌跡組成的段數有較強的自適應性，即使分開的兩段軌跡之間的區域絕對殘差很小也不會產生虛假軌跡，可很好地同時兼顧計算軌跡的完整性與精細性；計算參數易于設置且不必隨震源軌跡幾何形態而頻繁修改, 適于大批量事件的自動處理.由此可知，本文基于插值技術計算震源軌跡的方法具有較強的實用性.

6 結論

本文提出的基于插值技術計算殘差零值點為參考點的震源軌跡計算方法，不僅保持了常規方法(以震源軌跡節點為參考點)的優點，如適用于復雜模型計算、對震源軌跡的組成段數及穩定性沒有限制等，而且即使剖分模型的網格較為稀疏，計算的震源軌跡也具有很高精度.以稀疏網格剖分模型，可大大減少模型單元數量，從而相應地減少了震源軌跡計算(特別是其中最為費時的地震波場計算)所需的計算機內存和CPU時間，對于大尺度模型和臺站較多的臺網，效果更為顯著.稀疏網格中震源軌跡快速高精度計算的實現，使利用震源軌跡對區域和全球范圍內的地震事件進行快速高精度的定位成為可能.

數值實驗表明：莫霍面反射縱波(PmP)對地殼中事件震源的約束和直達縱波(Pg)相似，但與后者不同的是，前者的震源軌跡除了較集中地交會于震源點附近外，在莫霍面之下也有較集中的交會；同一臺站PmP-Pg波到時差約束的震源軌跡在震源附近以水平分布為主，對震源深度有較好約束；聯合使用PmP-Pg波到時可更好地約束震源位置、減少隨機干擾的影響；相比于遠臺，近臺對震源位置特別是震源深度有更好的約束.

致謝中國地震局地球物理研究所的吳慶舉研究員對本項工作給予了支持、房立華研究員在真實事件數據獲取方面給予了幫助；美國普林斯頓大學的劉前程博士提出了建設性意見；作者同時也感謝評審專家的中肯評述.