一種基于插值技術(shù)高精度計(jì)算稀疏網(wǎng)格地震定位中震源軌跡的方法

趙愛(ài)華

中國(guó)地震局地球物理研究所, 北京 100081

0 引言

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

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

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

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

1 震源軌跡方程

1.1 到時(shí)差約束的震源軌跡方程

(1)

由方程(1)可得到時(shí)差殘差(RDT)場(chǎng)FRDT：

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

(2)

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

1.2 到時(shí)約束的震源軌跡方程

(3)

(4)

到時(shí)殘差(RAT)場(chǎng)FRAT為

(5)

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

值得指出的是，到時(shí)約束的震源軌跡不僅與作為約束條件的觀測(cè)到時(shí)有關(guān)，還與構(gòu)建發(fā)震時(shí)間場(chǎng)的觀測(cè)到時(shí)有關(guān).發(fā)震時(shí)間場(chǎng)不同，同一到時(shí)約束的震源軌跡也不同，甚至相差很大(趙愛(ài)華等，2015).

2 震源軌跡參考點(diǎn)的計(jì)算

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

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

圖1 均勻介質(zhì)中事件以同一臺(tái)站P波與S波到時(shí)差約束的震源軌跡及其殘差場(chǎng)棕色點(diǎn)為震源軌跡所經(jīng)過(guò)模型單元的中心點(diǎn)；紅色十字為使用插值方法計(jì)算的震源軌跡參考點(diǎn)；藍(lán)色實(shí)線是震源軌跡的理論解,綠色虛線為±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)組成點(diǎn)對(duì)

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

(2)確定法線點(diǎn)對(duì)

選取同一模型節(jié)點(diǎn)其全部點(diǎn)對(duì)中絕對(duì)梯度最大者為法線點(diǎn)對(duì).法線點(diǎn)對(duì)中絕對(duì)殘差較小的點(diǎn)更靠近震源軌跡，震源軌跡經(jīng)過(guò)其所在模型單元，即該節(jié)點(diǎn)為震源軌跡節(jié)點(diǎn).

(3)插值求解

使用法線點(diǎn)對(duì)及其兩側(cè)延長(zhǎng)線上的模型節(jié)點(diǎn)構(gòu)建殘差的插值函數(shù)，選取插值函數(shù)中殘差為零的點(diǎn)作為軌跡參考點(diǎn).

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

(6)

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

3 射線路徑的計(jì)算

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

4 模型算例

4.1 速度模型

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

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

圖2 以云南地區(qū)為背景的速度模型(修改自白志明和王椿鏞,2004)紅色十字為震源位置；黑色倒三角形為地震臺(tái)站；灰線為地震反射/折射波射線路徑.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 到時(shí)差約束的震源軌跡

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

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

圖3 圖2中事件以到時(shí)差約束的震源軌跡(a) 不同臺(tái)站PmP波到時(shí)差； (b) 相同臺(tái)站PmP波與SmS波到時(shí)差; (c) 相同臺(tái)站PmP波與Pg波到時(shí)差.不同顏色和線型的曲線表示不同震源軌跡，其旁邊的數(shù)字“i”或數(shù)字對(duì)“i-j”表示震源軌跡以臺(tái)站Ri或臺(tái)站對(duì)Ri-Rj的到時(shí)差約束.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.

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

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

4.3 到時(shí)約束的震源軌跡

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

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

圖5 圖2中事件不同觀測(cè)情況以到時(shí)約束的震源軌跡(a)和(b) 所有臺(tái)； (c)和(d) 近臺(tái)； (e)和(f)遠(yuǎn)臺(tái)； (g)和(h) 右側(cè)臺(tái)； (i)和(j) 擾動(dòng)到時(shí). 左邊圖僅使用PmP波到時(shí)；右邊圖使用PmP和Pg波到時(shí)，其中藍(lán)線和綠線分別對(duì)應(yīng)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.

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

在區(qū)域臺(tái)網(wǎng)的地震定位中，事件常常位于臺(tái)網(wǎng)之外.為模擬該種情形，僅保留右側(cè)5個(gè)臺(tái)站R5-R9.PmP/Pn波到時(shí)約束的震源軌跡如圖5g所示.震源軌跡在震源處以和水平方向成約135°(以逆時(shí)針?lè)较驗(yàn)檎?的方位成束分布，對(duì)震源位置約束較弱.增加Pg波到時(shí)信息后，震源軌跡(圖5h)在震源處的方位分布得到顯著改善，對(duì)震源位置特別是震源深度的約束顯著加強(qiáng).

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

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

5 討論

5.1 計(jì)算誤差

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

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

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

圖6 圖1中到時(shí)差殘差沿震源軌跡法線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 計(jì)算效率

由表1可以看出，粗網(wǎng)格時(shí)使用插值法計(jì)算軌跡參考點(diǎn)可達(dá)到細(xì)網(wǎng)格時(shí)震源軌跡節(jié)點(diǎn)法的精度.以最簡(jiǎn)單的線性插值為例，d=9 km時(shí)其軌跡參考點(diǎn)的精度比震源軌跡節(jié)點(diǎn)法d=1 km時(shí)的精度還高.對(duì)于給定模型，d=1 km時(shí)的模型單元數(shù)是d=9 km時(shí)的81倍(二維)或729倍(三維).在震源軌跡計(jì)算中，地震波走時(shí)場(chǎng)的計(jì)算最為費(fèi)時(shí)，計(jì)算時(shí)間和模型單元數(shù)成正比.這意味著，使用新的軌跡參考點(diǎn)計(jì)算方法可使用于計(jì)算地震波走時(shí)場(chǎng)的時(shí)間大大減少，盡管后者的計(jì)算與地震事件無(wú)關(guān)，對(duì)于穩(wěn)定區(qū)域僅需1次，但對(duì)于尺度較大的三維模型或當(dāng)?shù)卣鹋_(tái)站較多時(shí)，這種減少的效果還是非常顯著的.

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

表2 震源軌跡參考點(diǎn)計(jì)算對(duì)計(jì)算效率的影響Table 2 Effect of the calculation of FLRPs on the efficiency

在射線路徑計(jì)算部分，需要將路徑中模型節(jié)點(diǎn)替換為計(jì)算的軌跡參考點(diǎn).由于操作簡(jiǎn)單，而且當(dāng)使用稀疏網(wǎng)格時(shí)，需要替換的節(jié)點(diǎn)較少，因此這部分所花費(fèi)的時(shí)間微乎其微，對(duì)軌跡的計(jì)算效率影響不大.

5.3 算法的實(shí)用性

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

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

圖7 華北地區(qū)2012年8月26日ML4.2地震的震源軌跡(a) P波到時(shí)差約束； (b) P和S波到時(shí)差約束； (c) P波到時(shí)約束； (d) P波和S波到時(shí)約束. 黑色線為趙愛(ài)華等(2015)計(jì)算結(jié)果， 彩色線為本文計(jì)算結(jié)果(詳情見(jiàn)正文). 倒三角形為地震臺(tái)站, 紅色十字為北京臺(tái)網(wǎng)所定震源位置在臺(tái)站連線上的投影.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.

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

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

圖8 圖7d中綠色實(shí)線震源軌跡(對(duì)應(yīng)CAD臺(tái)站P波走時(shí))的殘差場(chǎng)(a)及其局部放大圖(b)棕色點(diǎn)為使用線性插值計(jì)算的震源軌跡參考點(diǎn);菱形為震源軌跡兩側(cè)相鄰的模型節(jié)點(diǎn)(間距0.2 km),其中藍(lán)色的表示殘差小于0,紅色的表示殘差不小于0, 空心的表示殘差絕對(duì)值大于0.2 s,實(shí)心的表示殘差絕對(duì)值不大于0.2 s. 綠色線為到時(shí)殘差等值線(單位為秒).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).

真實(shí)事件的計(jì)算結(jié)果表明：本文算法具有良好的穩(wěn)健性，即使震源軌跡所在區(qū)域殘差變化劇烈，也能得到可靠結(jié)果；對(duì)震源軌跡組成的段數(shù)有較強(qiáng)的自適應(yīng)性，即使分開(kāi)的兩段軌跡之間的區(qū)域絕對(duì)殘差很小也不會(huì)產(chǎn)生虛假軌跡，可很好地同時(shí)兼顧計(jì)算軌跡的完整性與精細(xì)性；計(jì)算參數(shù)易于設(shè)置且不必隨震源軌跡幾何形態(tài)而頻繁修改, 適于大批量事件的自動(dòng)處理.由此可知，本文基于插值技術(shù)計(jì)算震源軌跡的方法具有較強(qiáng)的實(shí)用性.

6 結(jié)論

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

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

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