趙愛華, 丁志峰, 白志明

1 中國地震局地球物理研究所, 北京 100081 2 中國科學(xué)院地質(zhì)與地球物理研究所， 北京 100029

?

基于射線追蹤技術(shù)計算地震定位中震源軌跡的改進方法

趙愛華1, 丁志峰1, 白志明2

1 中國地震局地球物理研究所, 北京 100081 2 中國科學(xué)院地質(zhì)與地球物理研究所， 北京 100029

使用震源軌跡確定震源位置不僅穩(wěn)健而且直觀，但當(dāng)介質(zhì)復(fù)雜時震源軌跡難以給出解析解.基于最小走時樹射線追蹤技術(shù)計算震源軌跡的方法(以軌跡所在的殘差場中殘差最小的點(初始點)至殘差較小的點(震源軌跡代表點)的射線路徑表示震源軌跡)適用于復(fù)雜速度模型，但尚不能正確計算由多段組成的震源軌跡，同時兼顧計算軌跡的完整性和精細(xì)性較為困難，計算參數(shù)設(shè)置煩瑣不適于大批量數(shù)據(jù)的自動處理.針對該方法存在的問題，本文對其進行了改進：(1)采用一種“削皮”算法選取震源軌跡所經(jīng)過的模型單元的節(jié)點作為軌跡代表點；(2)將殘差較小的區(qū)域作為震源軌跡計算區(qū)域(該區(qū)域依軌跡分布自適應(yīng)地劃分為若干個連通區(qū)域)，從未計算的軌跡代表點中選取殘差最小者作為射線路徑初始點，利用最小走時樹算法依次計算所有連通區(qū)域內(nèi)的震源軌跡；(3)通過去掉較短的不再分叉的射線路徑使震源軌跡更為精細(xì).虛擬和真實事件的算例表明，改進方法有效克服了原方法的不足，可便捷地計算復(fù)雜速度模型中事件的震源軌跡，計算的軌跡精細(xì)且較完整.

地震定位； 震源軌跡； 射線追蹤； 最小走時樹方法； “削皮”算法

1 引言

地震定位在地震活動性分析、核爆監(jiān)測及地球深部結(jié)構(gòu)的研究等中起著十分重要的作用(曾融生, 1991；丁志峰等，1999；Richards et al., 2006；陳運泰，2009).為更好地確定震源位置與發(fā)震時間參數(shù)，地震學(xué)家發(fā)展了多種定位方法(田玥和陳曉非，2002；蔡明軍等，2004；楊文東等，2005)，例如蓋革法(Geiger, 1912)、雙差法(Waldhauser and Ellsworth, 2000)、網(wǎng)格搜索法(Shearer, 1997)、交切法(Pujol, 2004; 廉超等, 2006；周建超和趙愛華，2012)以及逆時成像法(許力生等，2013a, 2013b)等.這些方法各有特點，其中交切法不僅原理簡單，而且穩(wěn)健直觀，在地震臺網(wǎng)中有廣泛的應(yīng)用(孟玉梅等, 2001; 宋銳等, 2001).交切法使用震源軌跡確定震源位置，即使僅有少量的觀測數(shù)據(jù)，也能得到有價值的震源位置信息(趙愛華等, 2008).交切法的缺點是定位精度較低，特別是震源深度(Udías, 1999).交切法定位精度低的主要原因是它基于均勻或橫向均勻介質(zhì)模型(此時震源軌跡可準(zhǔn)確求解)，而實際上地球內(nèi)部存在著較強的橫向非均勻性.使用遠(yuǎn)離實際的速度模型必然導(dǎo)致或大或小的定位誤差(Billings et al., 1994; Chen et al., 2006).因此，提高交切法定位精度的關(guān)鍵是使之基于更接近實際的速度模型.此外，得益于地震臺網(wǎng)建設(shè)和反演技術(shù)發(fā)展，愈加廣泛的區(qū)域建立了較精細(xì)的三維結(jié)構(gòu)模型，全球尺度圖像的分辨率也不斷提高，使基于三維速度模型進行高精度定位成為可能.

將傳統(tǒng)交切法擴展至三維速度模型需要解決的首要問題是震源軌跡計算.當(dāng)速度結(jié)構(gòu)復(fù)雜時，震源軌跡難以表達成解析形式.對于復(fù)雜介質(zhì)中震源軌跡的計算問題，目前主要有兩種解決方案.一是以近似滿足震源軌跡方程的點表示震源軌跡，近似的程度以方程容差表征(Zhou，1994；Font et al., 2004；Theunissen et al., 2012).這種方法簡單易行，自適應(yīng)性強，對震源軌跡的幾何形態(tài)沒有限制，其缺點是得到的震源軌跡較為粗略，特別是當(dāng)方程容差取得較大時：即使到時和速度模型都準(zhǔn)確，震源軌跡也不會交于一點，而是交會成一個區(qū)域.另外一種方案基于最小走時樹射線追蹤技術(shù)：震源軌跡以其所在殘差場中殘差最小的點(初始點)至殘差較小的點(作為軌跡代表點)的射線路徑進行表示(趙愛華等，2008).這種方法可精細(xì)計算震源軌跡，根據(jù)震源軌跡的方位分布可直觀地評估定位效果，分析速度模型、到時拾取精度及臺站分布等因素對定位結(jié)果的影響.其缺點是：1)不能處理由多段組成的震源軌跡.由于僅選取殘差場中最小的點為射線路徑初始點，所有軌跡代表點均與該點相連，當(dāng)震源軌跡由分開的多段組成時會產(chǎn)生虛假軌跡(即穿過兩段軌跡之間區(qū)域的射線路徑)；2)選取軌跡代表點煩瑣費時，不適于大批量數(shù)據(jù)的自動處理.震源軌跡代表點按全部模型節(jié)點的百分比選取，軌跡代表點的百分比應(yīng)依軌跡長度而設(shè)，而軌跡長度很難預(yù)知；3)難以同時兼顧計算軌跡的完整性與精細(xì)性.即使準(zhǔn)確知道軌跡長度，由于軌跡所經(jīng)過模型單元的節(jié)點處的殘差并不相同，單純根據(jù)殘差大小選取的軌跡代表點要么不能完整表征軌跡要么冗余使軌跡較為粗略.現(xiàn)代地震臺網(wǎng)記錄的地震事件及觀測數(shù)據(jù)數(shù)量都是巨大的，大量事件的高精度定位，要求震源軌跡的計算不僅準(zhǔn)確而且效率要高.

由以上可知，復(fù)雜介質(zhì)中震源軌跡的計算方法還需完善.本文嘗試對基于射線追蹤技術(shù)的計算方法(趙愛華等，2008)進行改進以克服其不足.為敘述方便，方法說明以二維模型為例.

2 震源軌跡方程

2.1 到時差約束的震源軌跡

(1)

即考慮到計算誤差，震源軌跡在到時差殘差場

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

(2)

中具有接近于0的低值，其中I(x)=1.

2.2 到時約束的震源軌跡

根據(jù)觀測到時可構(gòu)建事件的發(fā)震時間場(趙愛華, 2014)

(3)

其中，Gj和Wj分別為第j個觀測到時的臺站和震相類型(P或S波)，K為觀測到時個數(shù).這樣，第j個到時約束的震源軌跡滿足如下方程：

(4)

即考慮到計算誤差，震源軌跡在到時殘差場

FRAT(x,Gj,Wj)=

(5)

中具有接近于0的低值，其中I(x)=1.

方程(1)和(4)，僅當(dāng)介質(zhì)速度結(jié)構(gòu)較為簡單時才存在解析解.與到時差約束的震源軌跡不同，到時約束的震源軌跡即使對于均勻介質(zhì)也比較復(fù)雜，特別是當(dāng)發(fā)震時間場使用較多到時構(gòu)建時.考慮P波和S波速度分別為6.0km·s-1和3.6km·s-1、尺度為150km×90km的均勻介質(zhì)模型.將模型剖分成3km×3km的正方形單元.設(shè)在點x0(37.5km,-31 .5km)處，時間t0=0 s時發(fā)生地震，記錄到該事件P和S波初至的臺站有5個，分別為Ri(13.5+(i-1)× 30.0 km,-1.5 km)(i=1,2,3,…,5).該事件的發(fā)震時間場及臺站R3處P波的到時殘差場灰度圖像如圖1所示.可以看到：1)P波到時單獨構(gòu)建的發(fā)震時間場(圖1a)與P和S波到時聯(lián)合構(gòu)建的發(fā)震時間場(圖1b)具有相似的結(jié)構(gòu)，這是因為P波與S波的速度之比為常數(shù)；2)單獨P波到時約束的震源軌跡(圖1c)、P和S波到時聯(lián)合約束的震源軌跡(圖1d)均較為復(fù)雜，既非圓亦非雙曲線，特別是后者由分開的兩段組成；3)對于同一到時，發(fā)震時間場不同，其震源軌跡也不同，甚至相差較大.需要指出的是，為突出顯示震源軌跡，到時殘差進行了取對數(shù)處理.

3 軌跡計算

由圖1(c,d)可以看出，殘差較小(例如小于0.2 s)的區(qū)域依震源軌跡的分布劃分為1個(圖1c)或多個(圖1d)連通區(qū)域.由于震源軌跡在殘差場中表現(xiàn)為低值，因此若在每個連通區(qū)域內(nèi)殘差最小點處激發(fā)地震波，地震波將沿該連通區(qū)域內(nèi)的震源軌跡快速傳播，而在軌跡兩側(cè)特別是連通區(qū)域外傳播較慢，此時將殘差看作地震波慢度.這樣，若選取部分單元節(jié)點作為震源軌跡的代表點，則代表點至與其同一連通區(qū)域內(nèi)地震波初始點的射線路徑可近似地表示震源軌跡，且不存在虛假震源軌跡問題.為處理方便，通常以殘差小于閾值的點作為震源軌跡的代表點(Zhou,1994;趙愛華等，2008)，由于震源軌跡所經(jīng)過節(jié)點其殘差并不相同(在圖1(c，d)中表現(xiàn)為灰度不同)，因而難以同時兼顧計算軌跡的完整性和精確性.與之相比，若以震源軌跡所經(jīng)過節(jié)點為其代表點，則計算震源軌跡即較完整也較精細(xì).基于以上認(rèn)識，我們提出改進的震源軌跡計算方法.

圖1 均勻介質(zhì)中事件的發(fā)震時間場(a和b)和到時殘差場灰度圖像(c和d)(a) P波到時構(gòu)建; (b) P與S波到時構(gòu)建; (c)和(d) 分別對應(yīng)(a)和(b)，同一P波到時約束. 發(fā)震時間的等值線間隔為2.0 s，紅色“+”為震源；到時殘差場中，綠線為0.2 s的到時殘差等值線，藍色實線和紅色虛線分別為震源軌跡的計算結(jié)果與理論解，黃色點為選定的震源軌跡代表點.Fig.1 Origin time fields (a, b) and arrival time residual fields (c, d) of an event in a homogeneous medium(a) P-wave arrival times; (b) P- and S-wave data; (c) and (d) one P-wave arrival time jointed with the origin time fields (a) and (b), respectively. The origin time contour interval equals 2.0 s, and the red crosses indicate the hypocenter. In the residual fields, the green lines are contours of 0.2 s. The blue solid and red dashed lines respectively show the calculated and theoretical focal loci. The yellow dots indicate the selected hypocentral locus reference points.

3.1 震源軌跡代表點的選取

在圖1(c，d)中，震源軌跡位于殘差較小連通區(qū)域內(nèi)的中間，猶如位于果中之核.震源軌跡所經(jīng)過的單元節(jié)點，其殘差盡管有大有小，但與軌跡兩側(cè)節(jié)點相比為最小.據(jù)此特點可采用類似層層削去果皮得到果核的方法求取這些節(jié)點.為便于敘述，將殘差小于RFLRP的單元節(jié)點稱為震源軌跡代表點的初選點，其集合記為Θ.從初選點集合Θ中選取震源軌跡代表點的步驟如下：

(1)分類

對于節(jié)點i∈Θ，若其鄰點均在集合Θ內(nèi)，則i為Θ的內(nèi)部點，內(nèi)部點集合記為I；否則i為Θ的外表點，外表點集合記為O.

(2)更新

對于節(jié)點i∈O，若其鄰點中沒有內(nèi)部點或有內(nèi)部點但殘差都較i點大，則保留該點；否則將i點從集合O中移出.之后，將集合Θ更新為O與I之并集.

(3)迭代條件檢測

若內(nèi)部點集合I為空集，則轉(zhuǎn)到步驟(4)；否則回到步驟(1).

(4)精簡

集合O中有時包括兩列沿震源軌跡并排分布的外表點.由于單列點足以表征震源軌跡，因此還需對集合O進行精簡.對于節(jié)點i∈O及與之相鄰的外表點，當(dāng)這些點的個數(shù)超過3個時，則僅保留其中殘差最小的3個.

(5)純化

當(dāng)殘差場的結(jié)構(gòu)較為復(fù)雜時，集合O中可能會有少數(shù)異常點，并非震源軌跡所經(jīng)過節(jié)點.這些異常點殘差相對較大.這里，將殘差超過集合O中點平均值2倍均方差的點視為異常點，并將其從集合O中移出.

在上述算法中，步驟(1)—(3)是用類似層層削皮的方式去掉初選點集合中的外表點；步驟(4)是去掉多余的代表點；步驟(5)是去掉假的代表點.最后保留的集合O中的點作為震源軌跡代表點.圖1(c，d)中黃色點是從殘差小于0.2s的初選點中選取的震源軌跡代表點.可以看出，選取的代表點僅為震源軌跡所經(jīng)過節(jié)點其中的一部分，但很好地表征了震源軌跡的形態(tài)，這對計算震源軌跡來說已足夠.實際上，以射線路徑表示震源軌跡，若能僅選取每段震源軌跡上殘差最小的點及兩端節(jié)點作為代表點更為理想，但實現(xiàn)較為困難.

3.2 射線路徑計算

表示震源軌跡的射線路徑利用改進的最小走時樹方法(Zhao et al., 2004)計算.最小走時樹射線追蹤方法(Nakanishi and Yamaguchi, 1986; Moser, 1991)基于地震學(xué)中的惠更斯原理和費瑪原理，具有穩(wěn)健性好、適用復(fù)雜模型的優(yōu)點.該方法通過構(gòu)建最小走時樹的方式計算地震波從震源至所有模型單元節(jié)點的走時和射線路徑(Cao and Greenhalgh., 1993).對于由分開的多段組成的震源軌跡，為防止產(chǎn)生虛假軌跡，需從包含每段軌跡的連通區(qū)域內(nèi)選出殘差最小的點作為地震波初始點.震源軌跡的形態(tài)分布難以預(yù)知，同時選出其每段軌跡的地震波初始點是困難的.由最小走時樹算法的終止條件可知其運行區(qū)域為連通區(qū)域.這樣，若將射線路徑的計算區(qū)域限制為殘差較小的區(qū)域，則可通過依次為其每個連通區(qū)域指定地震波初始點的方式，利用最小走時樹算法逐段計算出全部震源軌跡.注意到震源軌跡代表點包含了每個連通區(qū)域殘差最小的點，因此地震波初始點可從軌跡代表點中選取.

將射線路徑計算的區(qū)域限定為殘差小于RFL的區(qū)域Ω(這可通過將Ω外的節(jié)點設(shè)置為已得到最小走時的點來實現(xiàn)).從殘差小于RFLRP的單元節(jié)點中選取的震源軌跡代表點集合，這里記為Z.參數(shù)RFL和RFLRP可相同也可不同.表示震源軌跡的射線路徑計算步驟如下：

(1)查找地震波初始點

在震源軌跡代表點集合Z中查詢殘差最小的點s,將s從Z中移出；

(2)射線路徑計算

以s為地震波初始點，在殘差場中計算地震波走時及射線路徑直至最小走時樹算法終止；若節(jié)點j∈Z被計算過，則將其從Z中移出.

(3)迭代條件檢測

如果Z成為空集則結(jié)束；否則回到步驟(1).

3.3 震源軌跡修飾

震源軌跡代表點大多僅近似滿足軌跡方程.當(dāng)其數(shù)量過多時，計算的震源軌跡將不很精細(xì)，需要做進一步的修飾處理.

圖2 以射線路徑表示的震源軌跡示意圖Fig.2 A diagram of one hypocentral locus represented with ray paths

震源軌跡代表點至其各自地震波初始點的射線路徑形成樹狀結(jié)構(gòu)，如圖2所示.若將射線路徑(即圖2中有向線段)看作枝條，則射線經(jīng)過的節(jié)點則為枝條生發(fā)點，當(dāng)其生發(fā)的枝條數(shù)(圓圈內(nèi)數(shù)字)大于或等于2時該點為分杈點.軌跡代表點沿震源軌跡分布，因此射線路徑樹的主干(即射線路徑經(jīng)過次數(shù)nFLP較多的代表點至地震波初始點射線路徑，例如線段BCD)應(yīng)在實際震源軌跡上.如果僅取該部分(趙愛華 等，2008)，將可能會損害軌跡的完整性.例如，取nFLP=2，則將遺失圖2中的震源軌跡段GHI.為此，在對震源軌跡修飾時，我們是去掉長度Ltwig較小的不再分杈的枝條，例如圖2中虛線所示部分.

圖1(c，d)中的震源軌跡，Ltwig=2(以模型單元的邊長為單位)時的計算結(jié)果如圖中藍線所示.可以看出，計算的震源軌跡和理論結(jié)果(圖中紅色虛線)符合得很好.值得指出的是，若將模型剖分為更小的正方形單元(例如1 km×1 km)，則兩者將更為一致.

4 算例

4.1 虛擬事件

4.1.1 速度模型

虛擬事件以云南地區(qū)為背景.云南是我國中強地震多發(fā)地區(qū)，為此在該區(qū)開展了大量深部結(jié)構(gòu)探測活動(Kan et al., 1986; 王椿鏞等，2002).根據(jù)深地震測深資料，獲得了多條高分辨率地殼速度結(jié)構(gòu)剖面(白志明和王椿鏞, 2003; Zhang et al., 2011; Teng et al., 2013).對其中的遮放—賓川速度剖面(白志明和王椿鏞，2004)略做平滑處理，得到如圖3所示模型，S波與P波速度之比vS∶vP設(shè)為常數(shù)1∶1.73.可以看到，該模型在深度上和橫向上都存在較強的非均勻性.將模型剖分成1km×1km的正方形單元，模型單元內(nèi)速度為常數(shù).在地表布設(shè)9個地震臺站Ri(i=1, 2, 3,…,9)，其位置如表1所示.假設(shè)在點x0(130.5 km, -27.5 km)處發(fā)生地震，發(fā)震時間t0=0 s.震源至地震臺站的P和S波到時與射線路徑利用改進最小走時樹算法(Zhao et al., 2004)計算，射線路徑如圖3中黑線所示.

對圖3中事件，利用改進方法計算其分別以到時差和到時約束的兩種震源軌跡.到時差約束的震源軌跡每條僅涉及兩個觀測到時，具有較強的獨立性，研究得較多(白玲等，2003；廉超等，2006；Zhao and Ding, 2007, 2009; 趙愛華等，2008；Theunissen et al., 2012)；到時約束的震源軌跡不僅與每條對應(yīng)的到時有關(guān)，而且還與構(gòu)建發(fā)震時間場的到時有關(guān)，目前研究得還較少.因此，對于前種軌跡僅考慮全部臺站觀測、速度模型與到時均準(zhǔn)確的情況；對于后種軌跡則考慮不同臺站觀測(包括全部臺、稀疏臺、近臺、遠(yuǎn)臺和單側(cè)臺)、速度模型與到時均準(zhǔn)確以及全部臺站觀測、速度模型準(zhǔn)確、到時有擾動的情況.震源軌跡計算參數(shù)RFLRP=0.1 s,RFL=0.2 s,Ltwig=5.

圖3 復(fù)雜模型(修改自白志明等，2004) 震源(紅色十字)至臺站(黑色倒三角形)的聯(lián)線為地震波傳播路徑.Fig.3 A complex model (from Bai et al., 2004)The black lines are ray paths of seismic waves from the hypocenter (cross) to stations (inverse triangles).

表1 復(fù)雜速度模型中地震臺站位置Table 1 Positions of the seismic receivers in the complex model

4.1.2 到時差約束的震源軌跡

當(dāng)僅有部分臺站觀測時，例如只保留臺站R5～R9，由于較多分布方位缺失，上述兩類震源軌跡各自對震源的約束大大減弱.不過由于這兩類軌跡在分布方位上存在互補性，將兩者結(jié)合可在一定程度上彌補觀測上的不足.當(dāng)速度模型或到時存在擾動時，可以預(yù)計軌跡將不會恰好交匯于一點而是交匯成一個區(qū)域(參見趙愛華等，2008).

4.1.3 到時約束的震源軌跡

全部臺站P波到時約束的9條震源軌跡顯示在圖5a中.和到時差約束的震源軌跡(圖4a)相比，到時約束的震源軌跡在震源處的方位分布不是很全面和均勻，大部分軌跡靠近垂直方向分布，不過有兩條軌跡與垂直方向的夾角較大，因此仍能較好地約束震源位置.當(dāng)以全部臺站P波和S波到時構(gòu)建發(fā)震時間場時，P波(藍線)、S波(綠線)到時約束的震源軌跡如圖5b所示.增加S波到時信息后，震源軌跡與之前的方位分布非常相似，僅略有改善，不過有更多軌跡約束震源，可更好地減少抗隨機干擾的影響.

圖4 復(fù)雜速度模型中事件以到時差約束的震源軌跡(a) 不同臺站P波到時差約束的震源軌跡; (b) 相同臺站P波與S波到時差約束的震源軌跡.Fig.4 Hypocentral loci constrained with arrival time differences for the event in the complex model(a) Focal loci constrained with arrival time differences between P waves at different stations; (b) Those with ones between P and S waves at the same stations.

圖5 復(fù)雜速度模型中事件以到時約束的震源軌跡(a)和(b)所有臺; (c)和(d)稀疏臺; (e)和(f) 近臺; (g)和(h) 遠(yuǎn)臺; (i)和(j) 右側(cè)臺; (k)和(l)擾動到時. 左邊圖僅使用P波到時; 右邊圖使用P和S波兩種到時，其中藍線和紅線分別對應(yīng)P波、S波. 紅色“+”為震源位置.Fig.5 Hypocentral loci constrained with arrival times for the event in the complex model(a) and (b) all stations; (c) and (d) sparse stations; (e) and (f) near stations; (g) and (h) far stations; (i) and (j) right-side stations; (k) and (l) noisy arrival times. Focal loci in the left panels are constrained with only P- wave arrival times; those in the right ones with both P- and S- wave data, where the blue and green lines respond to P and S waves, respectively. The red “+” indicate the true hypocenter.

在全部臺站情況中，臺站間距為30或20 km, 而實際臺站間距通常要比之大1倍或更多.為此，僅保留序號為奇數(shù)的臺站R2i-1(i=1, 2, 3，…，5).稀疏臺P波及其與S波到時約束的震源軌跡分別顯示在圖5(c，d)中.可以看出，當(dāng)臺站間距增大后，震源軌跡在震源處的方位分布范圍略有減小，但仍能較好地約束震源位置.與之相似的是近臺(R2—R6)約束的震源軌跡(5(e，f)).當(dāng)僅有遠(yuǎn)臺(R1、R2、R7、R8和R9)時，無論單獨P波還是其與S波約束的震源軌跡(圖5(g，h))在震源處都近垂直分布，難以約束震源深度.

在定位實踐中，事件常常位于地震臺網(wǎng)之外.為模擬該種情形，僅保留右側(cè)5個臺站R5～R9.P波到時約束的震源軌跡如圖5i所示.震源軌跡在震源處以水平為主的方位分布，似乎顯示可很好地約束震源深度，而對震中位置約束較弱.增加S波到時信息后，震源軌跡在震源處的方位分布變得不僅全面而且均勻，如圖5j所示，這與前幾種情況顯著不同.這表明，在局部觀測條件下對地震進行定位時，聯(lián)合使用P波和S波到時尤其重要.此外，還可以看到震源軌跡除了在震源處交會外，還有幾條軌跡在震源的右上地方交會.這意味著總的到時殘差函數(shù)存在局部極小值，當(dāng)以蓋革法(以總的到時殘差為目標(biāo)函數(shù))定位時，所給初始震源參數(shù)須離真值較近.

震相到時受環(huán)境、儀器及人為等因素影響總會或多或少存在誤差.為考察到時拾取誤差對震源軌跡分布的影響，對全部9個臺站的理論到時加幅度為±0.5 s的隨機擾動(擾動利用FORTRAN90編譯器中的隨機函數(shù)生成)，P、S波的均方到時擾動量分別為0.356 s和0.321 s.P波到時單獨約束、P和S波到時聯(lián)合約束的震源軌跡分別如圖5(k，l)所示.可以看出，當(dāng)觀測到時有誤差時，震源軌跡不再交會于震源點，而是交會成一個小的區(qū)域.和圖5k相比，圖5l中震源軌跡交會最密集的點更接近實際震源位置.

4.2 真實事件

4.2.1 事件與速度模型

真實事件來自華北地區(qū).華北地區(qū)是我國重點地震監(jiān)測地區(qū)，該區(qū)地震臺站分布密集且較均勻.為在二維剖面上較好地考察震源軌跡分布，事件震中和臺站應(yīng)成線分布.2012年8月26日(北京時間)發(fā)生在華北東部的ML4.2級地震較符合這一要求.據(jù)北京數(shù)字遙測地震臺網(wǎng)(北京臺網(wǎng))測定，該事件震中位于117.42°E、39.61°N, 深13 km, 發(fā)震時間為7時13分34.3秒.圖6顯示了事件震中及其周圍臺站.需要指出的是，在更大的范圍內(nèi)還有更多觀測臺站.可以看出，該事件不僅有良好的方位觀測覆蓋，而且有近臺控制(CAD臺的震中距僅3 km).基于臺網(wǎng)分布分析(Bondár et al., 2004), 可認(rèn)為北京臺網(wǎng)的定位結(jié)果是可靠的.在由“重復(fù)地震”得到的北京臺網(wǎng)定位精度圖(蔣長勝等，2008)上，震中處的定位精度可達2km或更高.北京臺網(wǎng)基于走時表、采用單純形方法定位，觀測走時及走時殘差列于表2中.由表2中較小的走時殘差也可推斷北京臺網(wǎng)的定位結(jié)果較為可靠.此外，值得一提的是北京臺網(wǎng)和中國地震臺網(wǎng)中心的定位結(jié)果相當(dāng)一致: 震中相差不到1.5 km、震源深度相差3 km，發(fā)震時間相差0.1 s.這樣，可以認(rèn)為事件震中和臺站EWZ、CAD和XAZ基本位于直線AB上.

華北地區(qū)地殼三維速度結(jié)構(gòu)有較多的研究成果，例如嘉世旭等(2005)、齊誠等(2006)、于湘?zhèn)サ?2010).根據(jù)已有結(jié)果可知，沿直線AB方向，介質(zhì)橫向非均勻性較弱.另外，地震臺網(wǎng)采用1D模型進行定位仍較普遍.為此，本文仍采用1D模型.模型由三層介質(zhì)組成，從上至下各層厚度(H)、P波速度(vP)和S波速度(vS)為：H1=5 km、vP1=5.0 km·s-1、vS1=2.6 km·s-1;H2=15 km、vP2=6.1 km·s-1、vS2=3.5 km·s-1;H3=15 km、vP3=6.8 km·s-1、vS3=3.9 km·s-1.模型參照華北東部模型(滕吉文等，1979)、根據(jù)觀測走時構(gòu)建.模型使用0.1 km×0.1 km網(wǎng)格剖分.以北京臺網(wǎng)測定的震源位置為標(biāo)準(zhǔn)，本文模型的走時殘差如表2所示.由走時殘差來看，所構(gòu)建的模型是可以接受的.

表2 真實事件的觀測走時及走時殘差(單位：s)Table 2 Observed traveltimes and traveltime residuals for the real (unit: s)

圖6 真實事件的震中(紅色“+”)及觀測臺站(實三角形)Fig.6 Epicenter (red “+”) of a real event and the seismic stations (solid triangles)

4.2.2 震源軌跡

三個臺站EWZ、CAD和XAZ兩兩組合、P波到時差約束的震源軌跡如圖7a所示.不同軌跡以不同顏色表示.可以看出，P波到時差約束的震源軌跡近似為雙曲線；三條軌跡較好地交匯于一點.和北京臺網(wǎng)確定的震源位置(圖中紅色“+”)相比，交匯點略向右下偏.同一臺站P、S波到時差約束的震源軌跡顯示在圖7b中.從左至右三個臺站對應(yīng)的軌跡分別以紅、綠和藍色標(biāo)示.不難看出，盡管存在縱向非均勻性，軌跡仍大體為圓形，但交匯得不是很集中，特別是EWZ臺軌跡(紅色)偏離臺網(wǎng)所定震源較遠(yuǎn)，這可能主要與該臺Sg震相走時殘差較大有關(guān).

圖7 真實事件的震源軌跡(a) P波到時差約束； (b) P-S波到時差約束； (c) P波到時約束； (d) P和S波到時約束. 紅色“+”為北京地震臺網(wǎng)確定的震源位置.Fig.7 Hypocentral loci of the real event(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. The red signs of “+” indicate the hypocenter located by Beijing seismic network.

圖8 復(fù)雜模型中事件到時殘差(RAT)場中的震源軌跡(a)臺站R6的P波到時約束的震源軌跡, 發(fā)震時刻場使用所有臺站的P波到時構(gòu)建; (b)和(c) 臺站R9的P波到時約束的震源軌跡, 其發(fā)震時間場分別使用臺站R7—R9的P波到時單獨、P波和S波到時共同構(gòu)建. 紅色“+”為震源位置.Fig.8 Hypocentral loci in the fields of RAT (residual of arrival time) for the event in the complex model(a) The hypocentral locus is constrained with the P wave arrival time at station R6, and its origin time field is constructed with P-wave data from all the nine stations. (b) and (c) The hypocentral loci are constrained with the same P-wave arrival time at station R9, but their origin time fields are constructed with only P-wave, and both P-and S-wave data from stations R7—R9, respectively. The red “+” indicate the true hypocenter.

EWZ臺和XAZ臺軌跡存在中斷現(xiàn)象(中斷位于第一界面處)，這源于在兩個臺站處激發(fā)的地震波沿界面折射滑行，走時在界面處發(fā)生躍變.

P波到時約束的震源軌跡顯示在圖7c中.紅、綠和藍色軌跡分別對應(yīng)臺站EWZ、CAD和XAZ.軌跡2條似雙曲線、1條接近圓形.3條軌跡不僅較好地交匯于一點，而且方位分布較為全面，特別是對震源深度有很好地約束.據(jù)此推斷，實際震源可能要比臺網(wǎng)所定位置略深并偏向東北.P、S波到時聯(lián)合約束的震源軌跡如圖7d所示.紅、綠和藍色軌跡分別對應(yīng)臺站EWZ、CAD和XAZ，實線和虛線分別對應(yīng)P、S波.和圖7b類似，軌跡交匯得也不集中，但交匯區(qū)包含了臺網(wǎng)所定震源位置.和單純的P波相比，P、S波約束的震源軌跡交匯得較分散，這可能與S波的速度模型及/或拾取的震相到時誤差較大有關(guān).

5 討論

5.1 震源軌跡穩(wěn)定性

復(fù)雜模型的算例顯示，從震源軌跡的方位分布來看，遠(yuǎn)臺似乎有時更能約束震源深度(見圖4a和圖5i)，與常識明顯不符.實際上，震源軌跡對震源位置的約束，不僅要看軌跡在震源處的分布方位，而且還要考慮其穩(wěn)定性.若將震源軌跡所在的殘差場看作地形圖，則震源軌跡位于山谷中.山谷越窄陡，軌跡越穩(wěn)定；山谷越寬緩，軌跡穩(wěn)定性越差.

圖8a顯示了圖5a中標(biāo)有紅圈的震源軌跡的到時殘差場，該軌跡對應(yīng)臺站R6.可以看出，該軌跡由左、右兩段組成.和左段相比，右段殘差較小(例如小于0.1 s)的條帶要寬得多，穩(wěn)定性較弱.殘差場中藍線為RFLRP=0.01 s時計算的震源軌跡.和圖5a相比，紅圈標(biāo)示部位的軌跡略有不同，在形態(tài)上更接近實際.當(dāng)觀測到時有誤差時，穩(wěn)定性較弱的右段形狀發(fā)生很大變化，見圖5k紅圈標(biāo)示部分.與之相比，穩(wěn)定性強的左段形狀改變較小.這表明，穩(wěn)定性弱的震源軌跡易受干擾影響.

為更清楚地揭示遠(yuǎn)臺對震源深度的實際約束能力，考慮臺站R7、R8和R9到時約束的震源軌跡.圖8b顯示了發(fā)震時間場以P波到時構(gòu)建，臺站R9處P波到時約束的震源軌跡及其所在的殘差場.震源軌跡盡管在震源處近水平分布，但由于穩(wěn)定性差，特別是在震源附近，即使很小的干擾，其位置和形狀都可能發(fā)生很大變化，因此對震源深度的實際約束并不是很強.當(dāng)構(gòu)建發(fā)震時間場時增加S波到時信息，同一到時對應(yīng)的震源軌跡則穩(wěn)定得多，如圖8c所示.臺站R7和R8對應(yīng)的軌跡情況與之類似.因此，增加使用S波到時對震源位置約束的強化，不僅在于震源軌跡方位分布的改善，還在于軌跡穩(wěn)定性的提高.

5.2 震源軌跡計算參數(shù)的選取

震源軌跡代表點的選取參數(shù)RFLRP控制軌跡代表點選取范圍.為保證軌跡完整性，初選點應(yīng)包含軌跡各段端點所在模型單元的節(jié)點，即這些節(jié)點的最大殘差值為RFLRP的下限，通常地震波速度越小、模型單元尺寸越大，該值越大.另一方面，由于是通過不斷“削去”初選點集合外皮、保留其核心作為震源軌跡代表點，因此RFLRP不能大到初選點組成的連通區(qū)域有多個“核”，即每個連通區(qū)域應(yīng)僅包含1段震源軌跡.實際上，滿足上述要求的RFLRP取值范圍是較寬的.圖9顯示了均勻模型算例中P波到時約束的震源軌跡其初選(虛線)、實際(實線)代表點占所有模型節(jié)點百分比隨參數(shù)RFLRP的變化，RFLRP的步長取0.005s.可以看出，RFLRP在約0.03～0.5s的較寬范圍內(nèi)，震源軌跡代表點的選取結(jié)果都很穩(wěn)定，這種穩(wěn)定性意味著參數(shù)RFLRP不需要隨震源軌跡形態(tài)變化而調(diào)整，因而有利于大批量數(shù)據(jù)的自動處理.在保證軌跡完整性的前提下，RFLRP應(yīng)盡量取小，這樣不僅可減少計算量，而且對于穩(wěn)定性較弱的震源軌跡，可保證選取的軌跡代表點有較好的代表性.

圖9 圖1c中震源軌跡代表點(實線)及其初選點(虛線)百分比Fig.9 Percentages of primary (dashed line) and final (solid line) reference points for the focal locus in Fig.1c

射線路徑的計算參數(shù)RFL控制射線路徑的追蹤區(qū)域.為保證震源軌跡的完整性、且不出現(xiàn)虛假軌跡，殘差小于RFL的連通區(qū)域應(yīng)完整包含且僅包含1段軌跡，即軌跡所經(jīng)過模型單元節(jié)點的殘差最大值為RFL的下限，但RFL不能大到有連通區(qū)域包含多于1段震源軌跡.對于只有1段的震源軌跡，RFL的上限可取至包括整個模型空間；對于多段的震源軌跡，RFL的取值范圍通常也是比較寬的，以圖8a為例，RFL在0.1～0.5s范圍內(nèi)取值均可.對于較小的RFL，射線追蹤區(qū)域為1個或多個狹窄條帶，這些條帶占模型區(qū)域的比例很小，因之震源軌跡計算較快.

震源軌跡的修飾參數(shù)Ltwig控制所要去掉的射線路徑的長度，其取值大小與選取的軌跡代表點有關(guān).當(dāng)代表點均沿軌跡分布時，不再分叉的射線路徑較短，Ltwig可取得較小.對于穩(wěn)定性較弱的震源軌跡，當(dāng)參數(shù)RFLRP取得較大時，選取的軌跡代表點會包含代表性較差的點，此時Ltwig應(yīng)取得大些以去掉距軌跡較遠(yuǎn)點的射線路徑.

5.3 算法到三維介質(zhì)的擴展性

地震定位是個三維問題.在三維殘差場中，震源軌跡所經(jīng)過模型單元的節(jié)點殘差局部最小，這些節(jié)點位于殘差較小區(qū)域的中心，更類似于實際的果中之核.因此，選取軌跡代表點的“削皮”算法同樣適用于三維模型，只是在消除多余代表點的細(xì)節(jié)上會略有不同.射線路徑的計算基于最小走時樹算法，該算法對于二維和三維模型，具有相同的框架.在三維模型中，計算的射線路徑同樣具有樹狀結(jié)構(gòu)，因此仍然可以采用類似修剪樹枝的方法來修飾震源軌跡使之精細(xì).這樣，利用三維最小走時樹算法(Klimes and Kvasnicka, 1994; 王輝等，2000；趙愛華和張中杰, 2004)可較容易地將本文提出的改進算法推廣到三維模型.和二維模型相比，三維模型震源軌跡的計算量要大得多，這主要是計算三維走時場較為費時，不過這部分僅需計算一次.因此，對于速度結(jié)構(gòu)穩(wěn)定的地區(qū)，當(dāng)處理同一地震臺網(wǎng)記錄的多個事件時，震源軌跡的計算效率還是較高的.

6 結(jié)論

本文提出的計算震源軌跡的改進方法不僅保持了原方法適于復(fù)雜速度模型的優(yōu)點，而且可計算任意多段的震源軌跡，能較好地兼顧計算軌跡的完整性與精細(xì)性.射線追蹤僅限于殘差較小的區(qū)域，顯著減少了計算量.震源軌跡計算參數(shù)設(shè)置簡單，不需要隨軌跡形態(tài)分布反復(fù)調(diào)整，有利于大批量數(shù)據(jù)的自動處理.改進方法易于擴展到三維介質(zhì)模型.由于相對費時的走時場計算僅需一次，因此對于速度結(jié)構(gòu)穩(wěn)定地區(qū)、同一臺網(wǎng)記錄的多個地震，震源軌跡計算效率較高.

數(shù)值模型計算表明：到時差和到時約束的震源軌跡具有不同的形態(tài)分布，即使介質(zhì)均勻后者也較復(fù)雜；震源軌跡對震源位置的約束，不僅要看軌跡在震源處的方位分布，而且要考慮軌跡的穩(wěn)定性；近臺到時對震源深度有較強的約束；在速度模型及震相到時較準(zhǔn)的情況下增加S波到時可改善對震源位置的約束，對于臺網(wǎng)外地震尤為顯著；使用多條軌跡，有助于減少隨機因素對定位的影響.

致謝 中國科學(xué)院地質(zhì)與地球物理研究所張中杰研究員對本項工作給予了支持和指導(dǎo)；真實事件數(shù)據(jù)來自北京數(shù)字遙測地震臺網(wǎng)，田寶峰博士和房立華博士給予了熱情幫助.作者在此表示衷心感謝.

Bai L, Zhang T Z, Zhang H Z. 2003. Multiplet relative location and wave correlation correction and their application.ActaSeismologicaSinica(in Chinese), 25(6): 591-600.Bai Z M, Wang C Y. 2003. Tomographic investigation of the upper crustal structure and seismotectonic environments in Yunnan Province.ActaSeismologicaSinica(in Chinese), 25(2): 117-127.

Bai Z M, Wang C Y. 2004. Tomography research of the Zhefang-Binchuan and Menglian-Malong wide-angle seismic profiles in Yunnan province.ChineseJ.Geophys. (in Chinese), 47(2): 257-267.

Billings S D, Sambridge M S, Kennett B L N. 1994. Errors in hypocenter location: picking, model, and magnitude dependence.Bull.Seism.Soc.Am., 84(6): 1978-1990.Bondár I, Myers S C, Engdahl E R, et al. 2004. Epicentre accuracy based on seismic network criteria.Geophys.J.Int., 156(3): 483-496.

Cai M J, Shan X M, Xu Y, et al. 2004. Review of earthquake-locating methods from error.JournalofSeismologicalResearch(in Chinese), 27(4): 314-317.

Cao S H, Greenhalgh S. 1993. Calculation of the seismic first-break time field and its ray path distribution using a minimum traveltime tree algorithm.Geophys.J.Int., 114(3): 593-600.

Chen H, Chiu J M, Pujol J, et al. 2006. A simple algorithm for local earthquake location using 3DVPandVSModels: test examples in the central United States and in central eastern Taiwan.Bull.Seismol.Soc.Am., 96(1): 288-305.

Chen Y T. 2009. Earthquake prediction: Retrospect and prospect. Sci.ChinaSer.D-EarthSci. (in Chinese), 39(12): 1633-1658.Ding Z F, He Z Q, Sun W G, et al. 1999. 3-D crust and upper mantle velocity structure in eastern Tibetan plateau and its surrounding areas.ChineseJ.Geophys. (in Chinese), 42(2): 197-205.

Font Y, Kao H, Lallemand S et al. 2004. Hypocentre determination offshore of eastern Taiwan using the Maximum Intersection method.Geophys.J.Int., 158: 655-675.

Geiger L. 1912. Probability method for the determination of earthquake epicenters from the arrival time only.Bull.St.LouisUniv., 8: 60-71.

Jia S X, Qi C, Wang F Y, et al. 2005. Three-dimensional crustal gridded structure of the Capital area.ChineseJ.Geophys. (in Chinese), 48(6): 1316-1324.

Jiang C S, Wu Z L, Li Y T. 2008. Estimating the location accuracy of the Beijing Capital Digital Seismograph Network using repeating events.ChineseJ.Geophys. (in Chinese), 51(3): 817-827.

Kan R J, Hu H X, Zeng R S, et al. 1986. Crustal structure of Yunnan province, People′s Republic of China, from seismic refraction profiles.Science, 234(4775): 433-437.

Lian C, Li S L, Dong M, et al. 2006. Method of spherical surface shearing in earthquake location.JournalofGeodesyandGeodynamics(in Chinese), 26(2): 99-103.

Meng Y M, Zhao Y, Wang B. 2001. The analysis on deviation of rapid location by China Seismic observational Network.Earthquake(in Chinese), 21(3): 65-69.

Moser T J. 1991. Shortest path calculation of seismic rays.Geophysics, 56(1): 59-67.

Nakanishi I, Yamaguchi K. 1986. A numerical experiment on nonlinear image reconstruction from first-arrival times for two-dimensional island arc structure.J.Phys.Earth, 34(2): 195-201.

Pujol J. 2004. Earthquake location tutorial: graphical approach and approximate epicentral location techniques.Seismol.Res.Lett., 75(1): 63-74.Qi C, Zhao D P, Chen Y, et al. 2006. 3-D P and S wave velocity structures and their relationship to strong earthquakes in the Chinese capital region.ChineseJ.Geophys. (in Chinese), 49(3): 805-815.Richards P G, Waldhauser F, Schaff D, et al. 2006. The applicability of modern methods of earthquake location.PureandAppliedGeophysics, 163(2-3): 351-372.

Shearer P M. 1997. Improving local earthquake locations using the L1 norm and waveform cross correlation: Application to the Whittier Narrows, California, aftershock sequence.J.Geophys.Res., 102(B4): 8269-8283.Song R, Gu X H, Wang Y L, et al. 2001. The software system for real-time processing and large earthquake rapid determination of NCDSN.Earthquake(in Chinese), 21(4): 47-59.

Teng J W, Yao H, Zhou H N. 1979. Crustal structure in the Beijing-Tianjin-Tangshan-Zhangjiakou region.ActaGeophysicaSinica(in Chinese), 22(3): 218-236.

Teng J W, Zhang Z J, Zhang X K, et al. 2013. Investigation of the Moho discontinuity beneath the Chinese mainland using deep seismic sounding profiles.Tectonophysics, 609: 202-216.

Theunissen T, Font Y, Lallemand S, et al. 2012. Improvements of the Maximum Intersection Method for 3D Absolute Earthquake Locations.Bull.Seismol.Soc.Am., 102(4): 1764-1785.

Tian Y, Chen X F. 2002. Review of seismic location study.ProgressinGeophysics(in Chinese), 17(1): 147-155.

Udías A. 1999. Principles of Seismology. Cambridge: Cambridge University Press.Waldhauser F, Ellsworth W L. 2000. A double-difference earthquake location algorithm: method and application to the northern Hayward fault, California.Bull.Seismol.Soc.Am., 90(6): 1353-1368.

Wang C Y, Mooney W D, Wang X L, et al. 2002. Study on 3-D velocity structure of crust and upper mantle in Sichuan-Yunnan Region, China.ActaSeismologicaSinica(in Chinese), 24(1): 1-16.

Wang H, Chang X. 2000. 3-D ray tracing method based on graphic structure.ChineseJ.Geophys. (in Chinese), 43(4): 534-541.

Xu L S, Du H L, Yan C, et al. 2013a. A method for determination of earthquake hypocentroid: time-reversal imaging technique I- Principle and numerical tests.ChineseJ.Geophys. (in Chinese), 56(4): 1190-1206, doi: 10.6038/cjg20130414.Xu L S, Yan C, Zhang X, et al. 2013b. A method for determination of earthquake hypocentroid: time-reversal imaging technique II—An examination based on people-made earthquakes.ChineseJ.Geophys. (in Chinese), 56(12): 4009-4027, doi: 10.6038/cjg20131207.

Yang W D, Jin X, Li S Y, et al. 2005. Study of seismic location methods.EarthquakeEngineeringandEngineeringVibration(in Chinese), 25(1): 14-20.

Yu X W, Chen Y T, Zhang H. 2010. Three-dimensional crustal P-wave velocity structure and seismicity analysis in Beijing-Tianjin-Tangshan Region.ChineseJ.Geophys. (in Chinese), 53(8): 1817-1828, doi: 10.3969/j.issn.0001-5733.2010.08.007.

Zeng R S. 1991. The problems of lithospheric structure and geodynamics.AdvanceinEarthSciences(in Chinese), 6(2):1-10.

Zhang Z J, Yang L Q, Teng J W, et al. 2011. An overview of the earth crust under China.Earth-ScienceReviews, 104(1-3): 143-166.

Zhao A H, Zhang Z J, Peng S P. 2003. Fast ray tracing method for converted waves in complex media.JournalofChinaUniversityofMining&Technology(in Chinese), 32(5): 513-516.

Zhao A H, Zhang Z J. 2004. Fast calculation of converted wave traveltime in 3-D complex media.ChineseJ.Geophys. (in Chinese), 47(4): 702-707.Zhao A H, Zhang Z J, Teng J W. 2004. Minimum travel time tree algorithm for seismic ray tracing: improvement in efficiency.J.Geophy.Eng., 1(4): 245-251.

Zhao A H, Ding Z F. 2007. An intersection method for locating earthquakes in complex velocity models.AppliedGeophysics, 4(4): 294-300.

Zhao A H, Ding Z F, Sun W G, et al. 2008. Calculation of focal loci for earthquake location in complex media.ChineseJ.Geophys. (in Chinese), 51(4): 1188-1195.

Zhao A H, Ding Z F. 2009. Earthquake location in transversely isotropic media with a tilted symmetry axis.J.Seismol., 13(2): 301-311.

Zhao A H. 2104. Hypocentral loci in earthquake location under the constraint of arrival time residual norm. ∥ Chen Y T, Jin Z M, Shi Y L, et al. A Study of Interior Physics and Dynamics of Continental China — in Honor of Academician Teng Jiwen′s 60th Anniversary of being Engaged in Geophysical Research (in Chinese). Beijing: Science Press, 265-280.Zhou H W. 1994. Rapid three-dimensional hypocentral determination using a master station method.J.Geophys.Res., 99(B8): 15439-15455.Zhou J C, Zhao A H. 2012. An intersection method for locating earthquakes in 3-D complex velocity models.ChineseJ.Geophys. (in Chinese), 55(10): 3347-3354, doi: 10.6038/j.issn.0001-5733.2012.10.017.

附中文參考文獻

白玲, 張?zhí)熘? 張宏智. 2003. 多重相對定位法和波形相關(guān)校正及其應(yīng)用. 地震學(xué)報, 25(6): 591-600.

白志明, 王椿鏞. 2003. 云南地區(qū)上部地殼結(jié)構(gòu)和地震構(gòu)造環(huán)境的層析成像研究. 地震學(xué)報, 25(2): 117-127.

白志明, 王椿鏞. 2004. 云南遮放—賓川和孟連—馬龍寬角地震剖面的層析成像研究. 地球物理學(xué)報, 47(2): 257-267.

蔡明軍, 山秀明, 徐彥等. 2004. 從誤差觀點綜述分析地震定位方法. 地震研究, 27(4): 314-317.

陳運泰. 2009. 地震預(yù)測: 回顧與展望. 中國科學(xué)D輯: 地球科學(xué), 39 (12): 1633-1658.

丁志峰, 何正勤, 孫為國等. 1999. 青藏高原東部及其邊緣地區(qū)的地殼上地幔三維速度結(jié)構(gòu). 地球物理學(xué)報, 42(2): 197-205.

嘉世旭, 齊誠, 王夫運等. 2005. 首都圈地殼網(wǎng)格化三維結(jié)構(gòu). 地球物理學(xué)報, 48(6): 1316-1324.

蔣長勝, 吳忠良, 李宇彤. 2008. 首都圈地區(qū)“重復(fù)地震”及其在區(qū)域地震臺網(wǎng)定位精度評價中的應(yīng)用. 地球物理學(xué)報, 51(3): 817-827.

廉超, 李勝樂, 董曼等. 2006. 球面交切法地震定位. 大地測量與地球動力學(xué), 26(2): 99-103.

孟玉梅, 趙永, 王斌. 2001. 中國地震觀測臺網(wǎng)地震速報定位偏差的分析. 地震, 21(3): 65-69.

齊誠, 趙大鵬, 陳颙等. 2006. 首都圈地區(qū)地殼P波和S波三維速度結(jié)構(gòu)及其與大地震的關(guān)系. 地球物理學(xué)報, 49(3): 805-815.

宋銳, 顧小紅, 王永力等. 2001. 國家數(shù)字地震臺網(wǎng)中心實時處理及大地震速報軟件系統(tǒng). 地震, 21(4): 47-59.

滕吉文, 姚虹, 周海南. 1979. 北京、天津、唐山和張家口地區(qū)的地殼結(jié)構(gòu). 地球物理學(xué)報, 22(3): 218-236.

田玥, 陳曉非. 2002. 地震定位研究綜述. 地球物理學(xué)進展, 17(1): 147-155.

王椿鏞, Mooney W D, 王溪莉等. 2002. 川滇地區(qū)地殼上地幔三維速度結(jié)構(gòu)研究. 地震學(xué)報, 24(1): 1-16.

王輝, 常旭. 2000. 基于圖形結(jié)構(gòu)的三維射線追蹤方法. 地球物理學(xué)報, 43(4): 534-541.

許力生, 杜海林, 嚴(yán)川等. 2013a. 一種確定震源中心的方法: 逆時成像技術(shù)(一)—原理與數(shù)值實驗. 地球物理學(xué)報, 56(4): 1190-1206, doi: 10.6038/cjg20130414.

許力生, 嚴(yán)川, 張旭等. 2013b. 一種確定震源中心的方法: 逆時成像技術(shù)(二)—基于人工地震的檢驗. 地球物理學(xué)報, 56(12):4009-4027, doi: 10.6038/cjg20131207.

楊文東, 金星, 李山有等. 2005. 地震定位研究及應(yīng)用綜述. 地震工程與工程振動, 25(1): 14-20.

于湘?zhèn)? 陳運泰, 張懷. 2010. 京津唐地區(qū)地殼三維P波速度結(jié)構(gòu)與地震活動性分析. 地球物理學(xué)報, 53(8): 1817-1828, doi: 10.3969/j.issn.0001-5733.2010.08.007.

曾融生. 1991. 大陸巖石圈構(gòu)造與地球動力學(xué). 地球科學(xué)進展, 6(2): 1-10.

趙愛華, 張中杰, 彭蘇萍. 2003. 復(fù)雜地質(zhì)模型轉(zhuǎn)換波快速射線追蹤方法. 中國礦業(yè)大學(xué)學(xué)報, 32(5): 513-516.

趙愛華, 張中杰. 2004. 三維復(fù)雜介質(zhì)中轉(zhuǎn)換波走時快速計算. 地球物理學(xué)報, 47(4): 702-707.

趙愛華, 丁志峰, 孫為國等. 2008. 復(fù)雜介質(zhì)地震定位中震源軌跡的計算. 地球物理學(xué)報, 51(4): 1188-1195.

趙愛華. 2014. 地震定位中到時殘差范數(shù)下的震源軌跡. ∥ 陳運泰, 金振民, 石耀霖等. 中國大陸地球內(nèi)部物理學(xué)與動力學(xué)研究—慶賀滕吉文院士從事地球物理學(xué)研究60周年. 北京: 科學(xué)出版社, 265-280.

周建超, 趙愛華. 2012. 三維復(fù)雜速度模型的交切法地震定位. 地球物理學(xué)報, 55(10): 3347-3354, doi: 10.6038/j.issn.0001-5733.2012.10.017.

(本文編輯 汪海英)

Improvement of the ray-tracing based method calculating hypocentral loci for earthquake location

ZHAO Ai-Hua1, DING Zhi-Feng1, BAI Zhi-Ming2

1InstituteofGeophysics,ChinaEarthquakeAdministration,Beijing100081,China2InstituteofGeologyandGeophysics,ChineseAcademyofSciences,Beijing100029,China

Hypocentral loci are very useful for reliable and visual earthquake location, but they can hardly be analytically expressed when the velocity model is complex. There is a ray-tracing based numerical method to calculate them, in which a focal locus is represented in terms of ray paths from the minimum point (namely initial point) to low residual points (referred as hypocentral locus reference points, HLRPs) in its residual field. It has no restrictions on the complexity of the velocity model and can produce quite fine results. However, this method is incapable of addressing multi-segment loci and inadequate for processing large quantity of data. Additionally, it is rather laborious and difficult to set its controlling parameters for obtaining both fine and complete hypocentral loci. This paper presents an improvement of the ray-tracing based method.The method for calculating hypocentral loci is based on the minimum traveltime tree algorithm for tracing rays. It consists of three steps: (1) HLRPs are selected from nodes of the model cells that the hypocentral locus runs through by means of a so-called peeling method. (2) The calculation domain of a hypocentral locus is defined as such a low residual area that its connected regions each include one segment of the locus and then all the focal locus segments are respectively calculated with the minimum traveltime tree algorithm for tracing rays by repeatedly assigning the minimal point among those HLRPs that have not been traced as an initial point. (3) Short ray paths without branching are removed to make the calculated locus finer.The improved method is applied to a virtual seismic event and a real earthquake. The virtual event takes place in a 300 km×60 km complex velocity model with a background of Yunnan area. It is covered well by 9 seismic stations. The real earthquake occurred in North China and its epicenter and three seismic stations are located in an approximately straight line. Considering the lateral heterogeneity along the line is weak, a horizontally homogeneous crustal model is employed for calculating theoretical travel times. The differences between the observed and the calculated traveltimes are no more than 0.1 s for P waves and 0.4 s for S waves. For the two events, we calculate their hypocentral loci constrained with arrival times and those with arrival time differences when only P waves are used and both P and S waves are used. As for the arrival time constrained hypocentral loci of the virtual event, we investigate them in six different cases: (a) all stations, (b) sparse stations; (c) near stations; (d) far stations; (e) right-side stations; and (f) noisy arrival times. The numerical tests show that one-segment and multi-segment hypocentral loci are calculated correctly and their intactness is maintained well. For the same event, resultant fine hypocentral loci intersect exactly at the hypocenter when the velocity model and the observed arrival times are accurate. The parameters controlling the quality of outputting results can be chosen within a quite wide range and require little adjusting once they are set properly. Hypocentral loci associated with near stations have strong constraints on hypocenters. Additional use of S-wave data can improve the azimuth distribution and the stability of hypocentral loci.The improved method is capable of efficiently calculating hypocentral loci with good completeness and fineness for earthquakes in a complex model. Arrival times from near stations are critically important for determining hypocenters. S-wave data are helpful in strengthening the constraint on hypocenters, especially when the observation is incomplete.

Earthquake location; Hypocentral locus; Ray tracing; Minimum traveltime tree algorithm; Peeling method

趙愛華, 丁志峰, 白志明. 2015. 基于射線追蹤技術(shù)計算地震定位中震源軌跡的改進方法.地球物理學(xué)報，58(9):3272-3285,

10.6038/cjg20150922.

Zhao A H, Ding Z F, Bai Z M. 2015. Improvement of the ray-tracing based method calculating hypocentral loci for earthquake location.ChineseJ.Geophys. (in Chinese),58(9):3272-3285,doi:10.6038/cjg20150922.

10.6038/cjg20150922

P315

2014-09-17，2015-08-26收修定稿

國家自然科學(xué)基金項目(41374098,40974050),中央國家級公益事業(yè)單位基本科研業(yè)務(wù)費專項(DQJB11C03)資助.

趙愛華，男，1970年出生，研究員，2001年在中國科學(xué)院地質(zhì)與地球物理研究所獲理學(xué)博士學(xué)位，主要從事復(fù)雜介質(zhì)射線追蹤方法及其應(yīng)用研究.E-mail:ahzhao123@yahoo.com