國際地震中心改進(jìn)的定位程序

Istvn Bondr　Dmitry Storchak

?

國際地震中心改進(jìn)的定位程序

國際地震中心(ISC)是一個非政府的、非營利性機構(gòu)，它的首要任務(wù)是產(chǎn)出權(quán)威的地球地震活動性報告?！秶H地震中心公報》報道了50年(1960～2011年)的地震活動性。近年來由于全球臺站日益增加，使得報告的地震事件尤其是事件的震相數(shù)量急劇增加。相似的射線路徑會引起相關(guān)的走時預(yù)測誤差，這些誤差歸因于：模型中沒有加入地球非均勻性，這導(dǎo)致了對定位不確定性的低估和不利的臺網(wǎng)幾何結(jié)構(gòu)以及定位偏差。因此，全球的臺站分布越密集和越不均衡，大多數(shù)定位算法中所做出的假設(shè)(觀測數(shù)據(jù)是獨立的)越不成立。

計算地震學(xué)理論地震學(xué)

0　引言

國際地震中心，由于它的非官方身份，意味著在地震學(xué)上開展獨一無二的國際合作。國際地震中心是一個國際性的非營利機構(gòu)，它由美國國家自然科學(xué)基金會(NSF)、英國皇家學(xué)會、日本氣象廳(JMA)、俄羅斯科學(xué)院、中國地震局(CEA)、印度氣象局(IMD)、加拿大自然資源部和全球50個左右其他負(fù)責(zé)產(chǎn)出國家地震公報的單位如地震觀測臺、大學(xué)、官方的研究所和業(yè)務(wù)機構(gòu)支持。

國際地震中心的首要任務(wù)是產(chǎn)出權(quán)威的地球地震活動性報告。為達(dá)到該目標(biāo)，國際地震中心從全世界范圍的地震機構(gòu)收集地震公報(包括震源、斷層面解、矩張量、震級、烈度報告、震相拾取和振幅讀取)，來獲取每個地震事件最完整的一套臺站報告(Willemann and Storchak，2001)。國際地震中心的分析人員合并公報中的數(shù)據(jù)，并且重新定位和仔細(xì)修訂符合一定篩選標(biāo)準(zhǔn)的地震事件(例如，只有一個機構(gòu)上報的小事件我們不予修訂)來出版《國際地震中心公報》。即便如此，我們只修訂了將近20%的地震事件，每個上報的事件均存儲在國際地震中心的數(shù)據(jù)庫中，可從國際地震中心的網(wǎng)站(www.isc.ac.uk)上獲取。1999年之前上報的小事件，由于操作原因沒有被包含進(jìn)去。

修訂后的《國際地震中心公報》比實際時間晚24個月，因為國際地震中心要等到所有上報和臺站報告都收集齊(例如，有些機構(gòu)以年為單位發(fā)布它們的公報)。圖1表明近年來報告的地震數(shù)量和相關(guān)的事件震相呈指數(shù)級增加。數(shù)據(jù)量增長的現(xiàn)象，歸因于在過去的20年期間地方臺站、國家和區(qū)域臺網(wǎng)數(shù)量的激增。國際監(jiān)測系統(tǒng)(www.ctbto.org)臺網(wǎng)在報告震相方面也非常多產(chǎn)。最近非常密集的臺網(wǎng)(例如，美國臺陣(USArray，IberArray)已經(jīng)部署好并且開始向國際地震中心上報震相了。因此，國際地震中心要處理持續(xù)增長的數(shù)據(jù)量，同時還要保持《國際地震中心公報》的準(zhǔn)確度。

為了應(yīng)對不均衡和密集的臺網(wǎng)帶來的挑戰(zhàn)，我們考慮美國臺陣對智利—玻利維亞邊境區(qū)域遠(yuǎn)震事件定位的影響。圖2a展示了我們感興趣的區(qū)域和美國臺陣臺站的定位。由于美國臺陣臺站都位于一個感興趣區(qū)域很窄的方位帶，很多射線路徑是相似的，所以在地球三維速度結(jié)構(gòu)模型沒有給出的情況下，產(chǎn)生了相關(guān)的走時預(yù)測誤差。我們對發(fā)生在2007～2008年的20個大到足以記錄到的遠(yuǎn)震事件進(jìn)行了重定位。圖2b展示了采用美國臺陣和不采用美國臺陣的定位結(jié)果。我們假定地方臺網(wǎng)定位更準(zhǔn)確，盡管它們不像地面真實事件一樣準(zhǔn)確。因此，公平地講，使用從美國臺陣的臺站拾取的震相會使定位結(jié)果變得更差。如果忽略相關(guān)系統(tǒng)誤差，將導(dǎo)致誤差橢圓被低估(真實的定位沒有落在誤差橢圓中)并出現(xiàn)定位偏差。所以，一旦忽視了相關(guān)的模型誤差結(jié)構(gòu)，將會使自己陷入困境。事實上，正如圖2c所示，當(dāng)本文介紹的新的國際地震中心定位程序考慮了相關(guān)誤差，定位偏差就會顯著地減少。

國際地震中心從1964年成立至今，一直致力于為科學(xué)研究提供同質(zhì)量的地震公報。因此，國際地震中心采用的定位算法除了一些小的改進(jìn)外，40年來幾乎沒有改變(Bolt，1960；Adamsetal，1982)。然而，特別是在近十年來，引入了一些小的改進(jìn)。從2001年國際地震中心開始在定位中采用S，Sg，Sn和Sb震相(Storchaketal，2000)。據(jù)國際地震學(xué)與地球內(nèi)部物理學(xué)協(xié)會(Schweitzer and Storchak，2006)組織的幾個專題討論會的推薦，國際地震中心開始使用ak135速度模型(Kennettetal，1995)，用樣條曲線算法(Buland and Chapman，1983)計算ak135走時預(yù)測結(jié)果，來代替杰弗里斯-布倫走時表(Jeffreys and Bullen，1940)。這就緩解了P波、S波和PKP波的走時在杰弗里斯-布倫走時表基線上不同的問題(Dziewonski and Anderson，1981；Engdahletal，1998)。同時，為了達(dá)到震相識別的目標(biāo)，國際地震中心開始采用與ak135速度模型一致的國際地震學(xué)與地球內(nèi)部物理學(xué)協(xié)會標(biāo)準(zhǔn)震相表(Storchaketal，2003)。然而，為了與杰弗里斯-布倫走時表進(jìn)行比較，在本文中提出的定位算法改變之前，國際地震中心還沒有在定位算法中加入后續(xù)震相(除了震中距到100°的初至P波和S波震相)。

圖1　(a)每年報告的地震事件數(shù)量(灰色)和國際地震中心每年處理的地震事件數(shù)量(黑色)。(b)每年使用的震相數(shù)。報告的事件數(shù)和相關(guān)事件震相數(shù)都隨時間呈指數(shù)級增長。國際地震中心震源解數(shù)量的顯著減少反映了國際地震中心在程序上為避免小事件的重定位做出的改變，國際地震中心的分析在地方臺網(wǎng)操作員報告的定位上沒有實質(zhì)性改進(jìn)

為應(yīng)對從非常不均衡的臺網(wǎng)得到的持續(xù)增長的數(shù)據(jù)量帶來的挑戰(zhàn)，我們引入了新的國際地震中心定位程序以確保高效處理數(shù)據(jù)，并且進(jìn)一步改進(jìn)國際地震中心修訂事件的定位精度。新定位程序的核心是新的國際地震中心定位算法。本文中我們描述了新的國際地震中心定位算法的主要部分，并且通過重新定位一組分布在全球的地面真實事件的方式(Bondr and McLaughlin，2009b)做了驗證測試，同時重定位了整個《國際地震中心公報》。注意，討論國際地震中心的分析修訂程序超出了本文的范圍，本文主要關(guān)注新的國際地震中心定位算法。

1　新的國際地震中心定位算法

1.1目前國際地震中心采用的定位算法

國際地震中心當(dāng)前采用的定位算法的核心是杰弗里斯的“一致約化”法(Jeffreys，1932)，將走時殘差分布模擬為兩個混合的高斯分布，第二個高斯分布由一個明顯較大的方差來描述，這個方差解釋了由離群值(不好的拾取、震相識別錯誤等)造成的觀測粗尾。該方法是迭代式重新加權(quán)線性最小二乘算法(Anderson，1982)，在計算時通過減小離群值的權(quán)重來估計震源參數(shù)。每一步迭代都重新計算依賴于殘差的加權(quán)函數(shù)。Buland(1986)曾指出：對非線性問題引入非線性加權(quán)函數(shù)使得一個原本就非線性的問題變得更非線性了。

在該定位算法中我們只考慮Pg，Pb，Pn，Sg，Sb和Sn波(到8°)以及P波和S波(到100°，當(dāng)S波超過20°時，我們將權(quán)值以2為因子降低)。我們識別和計算一些其他震相(PKP分支震相、深度震相等)的殘差，但它們不用于定位。形式不確定性用95%的置信水平來衡量。模型誤差(由沒有被模型描述的地球速度非均勻性導(dǎo)致的走時預(yù)測值的系統(tǒng)誤差)沒有做出規(guī)定，并且我們認(rèn)為觀測值是獨立的且服從正態(tài)分布。

上報的震源按其質(zhì)量劃分等級，該質(zhì)量高低基于依次增加的方位間隔，而這又取決于定位所用的臺站分布。在定位過程中，定位算法首先試圖用各個按等級劃分后的上報震源來找到自由深度解，作為初始值。如果算法找到了一個收斂的解，算法就會轉(zhuǎn)而求解體波和面波震級。若無法收斂，則重復(fù)該過程，但這時我們會用一個深度固定的試驗震源。如果沒有找到解，再次重復(fù)該過程，且這時的深度會限定為缺省值(10km或35km)，由各自的Flinn-Engdahl(Youngetal，1996)的區(qū)域編號指定。

國際地震中心計算體波和面波震級的條件是：一個事件有足夠數(shù)量的上報的振幅和周期讀取資料。對于體波，臺站的mb是利用Gutenberg和Richter(1956)公式，由上報的震中距在21°～100°之間、周期3s的P波的振幅計算得到的。臺站的MS是利用布拉格公式(Vaneketal，1962)，由周期在10～60s，事件的震源深度60km，震中距在20°～100°讀取的面波振幅計算得到。臺網(wǎng)的震級是相應(yīng)的臺站震級取平均得到的。注意，甚至用一個單獨臺站的震級就可以產(chǎn)出一個臺網(wǎng)的震級而且不計算誤差估計。

下面闡述的驗證測試中，我們利用當(dāng)前的國際地震中心定位算法作為衡量新定位算法性能的基準(zhǔn)。

1.2國際地震中心新的定位算法

雖然當(dāng)前的國際地震中心的定位方法已經(jīng)為科學(xué)界提供了50年的服務(wù)，但它確實需要改進(jìn)了。杰弗里斯的“一致約化”方法在只利用初至P波(和可能的S波)震相時是最好的。因為它平等地對待所有震相。為保證后續(xù)震相(觀測可靠性低)權(quán)重降低，結(jié)合后續(xù)震相時需要進(jìn)一步調(diào)整。對觀測加入一個相對加權(quán)結(jié)構(gòu)，不僅會使算法復(fù)雜化，還會牽制加權(quán)迭代算法的效率。

線性定位算法對定位的初始點非常敏感。當(dāng)前的程序做出了這樣的假設(shè)：在上報的震源中有一個好的初始震源可以使用，雖然上報的任何一個震源在搜索空間中接近全局最小是沒有保障的。而且，當(dāng)數(shù)據(jù)在深度上沒有分辨力時(例如，當(dāng)初至沒有在P波走時曲線的拐點之內(nèi)時)，試圖找到一個自由深度解是徒勞的。當(dāng)沒有深度解時，算法會簡單地從發(fā)震時間-深度的權(quán)衡曲線上選一個點。目前的國際地震中心定位算法假設(shè)觀測誤差是獨立的。正如我們之前指出的一樣，隨著臺網(wǎng)變得越來越密集，如果我們想提高(或簡單地維持)定位精度，就不可避免地要考慮相關(guān)的走時預(yù)測誤差。最終，發(fā)布的臺網(wǎng)震級可能來源于對一個單獨臺站的測量，這樣更容易導(dǎo)致產(chǎn)出錯誤的事件震級估計。

為解決這些問題，我們開發(fā)了一套定位算法：

(1)定位中采用所有ak135(Kennettetal，1995)預(yù)測的震相(包括深度震相)，并對高程、橢率(Dziewonski and Gilbert，1976；Kennett and Gudmundsson，1996；Engdahletal，1998)和深度震相的反射點(Engdahletal，1998)進(jìn)行校正。

(2)通過鄰域算法得到初始震源估計(NA；Sambridge 1999；Sambridge and Kennett，2001)。

(3)采用一種全部數(shù)據(jù)協(xié)方差矩陣來解釋相關(guān)模型誤差的先驗估計，并執(zhí)行迭代線性反演(Bondr and McLaughlin，2009a)。

(4)當(dāng)且僅當(dāng)有深度解時，嘗試一個自由深度解，否則將深度固定為一個取決于區(qū)域的缺省深度。

(5)將不確定性衡量到90%的置信水平，并且以各種距離范圍計算定位質(zhì)量。

(6)基于對上報的地表反射波通過深度震相疊加(Murphy and Barker，2006)來獲得深度震相的深度估計。

(7)提供可靠的帶有不確定性的臺網(wǎng)震級估計。

1.2.1震相

利用ak135走時預(yù)測的最大優(yōu)勢之一是：與杰弗里斯-布倫走時表(Jeffreys and Bullen，1940)不同的是，它不受限于P波、S波和PKP震相之間基線不一致的問題。而且，ak135提供了大量的可用到定位中的國際地震學(xué)與地球內(nèi)部物理學(xué)協(xié)會標(biāo)準(zhǔn)地震表中的震相，最著名的是PKP波分支和深度敏感震相。在本文中，我們提到了定位中要用到的作為時間限定的震相。因為我們通常處理許多時間限定的震相，所以定位算法中不必包含慢度和方位角測量值。根據(jù)上報的震源參數(shù)的中位值來識別最初的震相。

圖2　(a)在定位智利—玻利維亞邊境地區(qū)的大小適度的地震時，美國臺陣臺站控制了臺網(wǎng)幾何結(jié)構(gòu)。(b)目標(biāo)區(qū)域放大圖。當(dāng)美國臺陣臺站(三角形)被包括在定位中時，相對于地方臺網(wǎng)(圓形)，解的定位偏差增加了10km。倒三角表示定位時沒有利用美國臺陣。(c)相對于地方臺網(wǎng)解的定位偏差的累積分布。在假設(shè)誤差獨立時，采用美國臺陣的臺站(虛線)使定位變得更糟糕。當(dāng)考慮相關(guān)誤差(實線)時，定位偏差明顯減小

利用WG84橢球參數(shù)做高程和橢率校正(Dziewonski and Gilbert，1976；Kennett and Gudmundsson，1996；Engdahletal，1998)，并將這兩種校正加入ak135走時預(yù)測。我們用Engdahl等(1998)的算法進(jìn)行深度震相、反射點(在地表的反射點高程校正)和水深(針對pwP震相)校正。同時用ETOPO1全球地形模型(Amante and Eakins，2009)獲得高程或反射點的水深。

震相拾取誤差被描述為一個先驗的測量誤差估計，它來源于對地面真實事件的殘差分布(殘差是通過對地面真實事件的定位計算得到)的觀測，該地面真實事件來自國際地震學(xué)與地球內(nèi)部物理學(xué)協(xié)會的參考事件表(Bondr and McLaughlin，2009b)。圖3的直方圖是地面真實事件相對于Pb，Pn，P和PKPdf震相的殘差。

對于那些沒有足夠數(shù)量地面真實事件的觀測震相，我們確定了一種先驗的測量誤差，以便保持相對加權(quán)模式的一致性。初至P波類型震相(P，Pn，Pb，Pg)的拾取比后續(xù)震相精確，所以它們的測量誤差最小，為0.8s。對于初至S震相(S，Sn，Sb，Sg波)，測量誤差被設(shè)定為1.5s。穿過或從地球內(nèi)/外核反射的震相的測量誤差估計稍大(PKP，PKS，PKKP，PKKS等P′P′分支，PKiKP，PcP和PcS為1.3s，SKP，SKS，SKKP，SKKS等S′S′分支，SKiKP，ScP，ScS為1.8s)是為了解釋不同震相分支中可能的識別誤差。自由表面的反射和轉(zhuǎn)換震相(PnPn，PbPb，PgPg，PS，PnS，PgS和SnSn，SbSb，SgSg，SP，SPn和SPg波)不常被觀測到，而且具有更大的不確定性，所以它們的測量誤差大，為2.5s。同樣，測量誤差為2.5s分配給了更長周期典型的突發(fā)衍射震相(Pdif，Sdif和PKPdif)。通常將觀測到的深度震相(pP，sP，pS，sS和pwP)的先驗測量誤差設(shè)為1.3s，其余的深度震相(pPKP，sPKP，pSKS，sSKS分支和pPb，sPb，sSb，pPn，sPn，sSn)的測量誤差設(shè)為1.8s。我們將不太可靠的深度震相(pPg，sPg，sSg，pPdif，pSdif、sPdif和sSdif)的測量誤差設(shè)為2.5s。注意，我們同樣允許距離依存的測量誤差。例如，考慮到震中距在15°～28°之間的遠(yuǎn)區(qū)域距離上可能的震相識別誤差，對于Pn波和P波我們將先驗測量誤差從0.8s提高為1.2s；Sn波和S波從1.5s提高為1.8s。對于清晰的PP波和SS波，震中距在40°～80°的測量誤差分別設(shè)定為1.3s和1.8s。但震中距在25°～40°時測量誤差增加為1.8s和2.5s。

圖3　地面真實事件相對于(a)Pb。(b)Pn。(c)P波和(d)PKPdf震相的殘差分布(相對于國際地震學(xué)與地球內(nèi)部物理學(xué)協(xié)會參考事件表中地面真實事件計算的殘差)。實線表示最佳擬合高斯分布。smad：標(biāo)準(zhǔn)中位值絕對偏差

上文描述的相對加權(quán)模式，在算法中降低了拾取震相的不確定性或更容易對震相識別誤差降低權(quán)重。最佳的是：測量誤差應(yīng)該取決于觀測信噪比，但除了IMS臺站外，其他上報的很難保證。自從國際地震中心利用上報的質(zhì)量參差不齊的參數(shù)數(shù)據(jù)以來，我們寧愿選擇一套更為保守的先驗測量誤差估計。

1.2.2相關(guān)的模型誤差結(jié)構(gòu)

絕大多數(shù)定位算法，線性的或非線性的，都假設(shè)所有觀測誤差都是獨立的。該假設(shè)在兩個臺站間的距離小于地方速度非均勻的長度尺度時不成立。當(dāng)出現(xiàn)相關(guān)走時預(yù)測誤差時，數(shù)據(jù)的協(xié)方差矩陣就不再是對角陣，并且觀測中的冗余減少了自由度的有效數(shù)量。因此，忽視了相關(guān)誤差結(jié)構(gòu)，不可避免地導(dǎo)致定位不確定性的低估。用不均衡的臺網(wǎng)定位的事件或許也會導(dǎo)致定位估計的偏差。

Chang等(1983)證明了在線性定位算法中，當(dāng)可以得到非對角數(shù)據(jù)協(xié)方差矩陣估計時，解釋相關(guān)誤差結(jié)構(gòu)是相對直接的。為確定數(shù)據(jù)協(xié)方差矩陣，我們采用了Bondr和McLaughlin(2009a)描述的方法。他們假設(shè)：臺站之間的距離可以由射線路徑的相似性來近似。這種簡化假設(shè)允許從通用的來自地面真實殘差的P波變差函數(shù)模型(圖4)對臺站對之間的協(xié)方差進(jìn)行估計。由于在《國際地震中心公報》中絕大多數(shù)的震相是遠(yuǎn)震的P波震相，我們期盼通用變差函數(shù)模型在全球任何地方都表現(xiàn)得相當(dāng)好。

數(shù)據(jù)協(xié)方差矩陣的元素的估計是這樣獲得的：

(1)

注意，置信誤差橢圓給出了震中估計的精度，但不指明定位的準(zhǔn)確性(偏差)。人們能獲得非常精確的定位參數(shù)估計，但由于速度模型中三維速度結(jié)構(gòu)的不清楚而會導(dǎo)致有偏差。定位估計的準(zhǔn)確性只能在可以獲得地面真實(GT)事件信息的情況下測量。我們定義位置范圍值為：

(2)式中x和y是笛卡爾坐標(biāo)系，在坐標(biāo)系中GTx震中(x代表了地面真實準(zhǔn)確性，0～5km)由誤差橢圓的半軸定義，誤差橢圓按比例縮放到所需的置信水平，并集中在地震定位上。因為地面真實震中也有誤差項(GTx)，

圖4　來自地面真實事件殘差的遠(yuǎn)震P波震相的通用變差函數(shù)模型。用各向同性、固定的變差函數(shù)模型產(chǎn)生先驗的數(shù)據(jù)協(xié)方差矩陣估計

為了解釋這個不確定性，在計算事件的可能位置時，我們擴大置信誤差橢圓。我們將地面真實震中的坐標(biāo)插入誤差橢圓(從定位算法中獲得，并根據(jù)地面真實不確定性擴大)的標(biāo)準(zhǔn)方程中，來計算定位覆蓋參數(shù)。如果覆蓋參數(shù)小于1，那么地面真實震中在誤差橢圓中；如果大于1，那么誤差橢圓就沒有覆蓋真實位置。

因為這種表示中協(xié)方差只依賴于臺站距離，協(xié)方差矩陣(和它的逆)只需計算一次。我們假設(shè)不同的震相由于不同的射線路徑是不相關(guān)的，距離大于1 000km的臺站對也不相關(guān)。因此數(shù)據(jù)的協(xié)方差矩陣是一個稀疏的分塊對角矩陣。而且，如果在各個震相分塊的臺站以它們最近的相鄰距離排序，那么震相分塊就會變?yōu)榉謮K對角的。為減少計算大矩陣逆的時間，我們開發(fā)了具有分塊對角的結(jié)構(gòu)，通過一塊一塊的方式求協(xié)方差矩陣。

協(xié)方差矩陣的奇異值分解為：

(3)

(4)

使數(shù)據(jù)集正交化并且把冗余觀測投影到零空間。

1.2.3深度的分辨

原則上，如果有從源發(fā)出的上行波和下行波的組合，也即，如果有臺站覆蓋了初至震相走時的垂向分量偏導(dǎo)數(shù)改變符號的距離范圍(地方臺網(wǎng))，或者如果震相有相反符號的垂向慢度(深度震相)，就可以確定出深度。地核反射震相，例如PcP波，以及在更小范圍內(nèi)，二次震相(尤其是S波)也能有助于分辨深度。

我們擬定了一些標(biāo)準(zhǔn)來測試一個地震事件的上報數(shù)據(jù)是否能充分地分辨深度：

(1)地方臺網(wǎng)：在0.2°以內(nèi)有一個或多個有確定時間的震相的臺站。

(2)深度震相：至少有兩個機構(gòu)(為了降低由沒有經(jīng)驗的分析人員導(dǎo)致的錯誤的幾率)上報5或5個以上有確定時間的深度震相。

(3)核反射震相：至少有兩個機構(gòu)上報5或5個以上有確定時間的核反射震相(PcP，ScS波)。

(4)地方/附近區(qū)域的S波震相：在5°以內(nèi)有5或5個以上有確定時間的S波和P波對。

如果滿足以上任一標(biāo)準(zhǔn)，我們就嘗試獲得一個讓深度自由變化的解答；否則我們將深度固定在一個依賴震中定位的缺省深度。注意，一個地方/區(qū)域震相和PKP波觀測之間的組合或許也有助于分辨深度，但我們還沒有應(yīng)用這個標(biāo)準(zhǔn)。已知的人工震源事件和滑坡的深度被固定在地表。

1.2.4深度震相疊加

我們將深度震相直接用在定位中，深度震相疊加方法(Murphy and Barker，2006)提供了一種獲得穩(wěn)定深度估計的獨立手段。圖5闡釋了該方法。對每個深度震相——初至P波對都計算時差曲線(以震中距為參數(shù))，并且畫出了以深度作為函數(shù)的圖。間隔由深度震相的測量誤差估計決定。該誤差估計再投影到深度上觀測的時差周圍，產(chǎn)生一條在深度上與觀測時差一致的矩形函數(shù)道。道的疊加表示所有觀測到的深度震相對應(yīng)的深度分布。由深度震相限制的深度由疊加中位值決定，并且它的不確定性用標(biāo)準(zhǔn)中位值絕對偏差(SMAD，L1范數(shù)等價的標(biāo)準(zhǔn)偏差)來刻畫。需要注意的是，如果深度不確定性與震中誤差有關(guān)，也即，震源誤差橢圓有很大的傾伏成分時，深度不確定性可能會被低估。

由于深度是從深度震相疊加方法中獲得的，這就隱含了深度依賴震中本身，我們只執(zhí)行兩次深度震相疊加：第一次，對初始定位，得到合理的下文中將介紹的網(wǎng)格搜索中的深度起始點；第二次，對最終定位，獲得由深度震相限制的深度的最終估計。

1.2.5初始震源

對于記錄不充分的事件，上報的震源可能呈現(xiàn)出大的發(fā)散并且會有較大的定位誤差，尤其是它們只作為遠(yuǎn)震被記錄到時。為了獲得好的初始震源估計，我們對線性定位算法采用了鄰域算法(NA)(Sambridge，1999；Sambridge and Kennett，2001)。鄰域算法是一個非線性的網(wǎng)格搜索方法，具有探索巨大搜索空間的能力，并能快速地趨近全局最優(yōu)值。Kennett(2006)討論了鄰域算法的細(xì)節(jié)和在定位地震中的應(yīng)用。

我們在上報的震源參數(shù)的中位值周圍采用寬泛定義的搜索區(qū)域——在震中中位值周圍以2°為半徑的圓內(nèi)執(zhí)行搜索，即在發(fā)震時刻中位值周圍10s，上報深度中位值上下150km。這些缺省的搜索參數(shù)是通過試錯法在上報的震源參數(shù)中位值中達(dá)到完成時間和粗差修正量之間的折衷獲得的。注意，如果我們對深度解的測試失敗，我們就將深度固定為區(qū)域依存的缺省深度。初始的震源估計會是一個在鄰域算法的實驗震源中有最小L1范數(shù)失配的值。一旦接近全局最優(yōu)值，我們便開始用線性定位算法來得到最終解和相應(yīng)的形式不確定性。

圖5　Murphy和Barke(2006)的深度震相疊加方法的說明。(a)每個臺站生成深度函數(shù)的(它們也依賴震中距)預(yù)測時差曲線(深度震相—初至P波走時)。(b)相對于每個觀測時差，將矩形函數(shù)放在相應(yīng)的深度上產(chǎn)生深度道。矩形函數(shù)的寬度由深度震相的先驗測量誤差定義；集中在觀測時差上然后投影到x軸(深度)。(c)疊加深度道。深度震相深度及其不確定性分別以中位值和疊加的標(biāo)準(zhǔn)中位值絕對偏差來定義

1.2.6迭代線性化定位算法

Gwm=WGm=Wd=dW

(5)

式中G是(N×M)矩陣，包括N數(shù)據(jù)對M模型參數(shù)的偏微分，m是(M×1)模型調(diào)解矢量[ΔT，Δx，Δy，Δz]T，d是時間殘差的(N×1)矢量，W是(N×N)方程(4)的投影矩陣。所以，我們在本征坐標(biāo)系下求解這個反演問題，轉(zhuǎn)換的觀測值是獨立的，也即，dW代表觀測殘差的線性組合，即“本征殘差”。方程(5)可由奇異值分解求解，產(chǎn)生廣義逆：

(6)

和模型調(diào)節(jié)：

(7)

第j次迭代后，模型矢量調(diào)整為mj+1=mj+mest。一旦獲得收斂解，定位不確定性由后驗的模型協(xié)方差矩陣定義：

(8)

模型協(xié)方差矩陣產(chǎn)生四維誤差橢球，它的推算提供二維誤差橢圓和一維的深度以及發(fā)震時刻的誤差。誤差橢圓包圍的給定百分之α水平的置信區(qū)間被定義為：

(9)

(10)

這里方差標(biāo)定因子s2被定義為：

(11)

Fα是一個在α=90%，具有M和K+N-M自由度、M=2和N個獨立觀測值的F統(tǒng)計值，也即，觀測值總數(shù)減去投影到零空間的觀測值的數(shù)量。K被設(shè)為很大的值(99 999)，以使形式不確定性估計近似“覆蓋”誤差橢圓(Evernden，1969)，其假設(shè)了已知精確的先驗誤差估計，因而降低了后驗殘差的作用。這也是國際數(shù)據(jù)中心產(chǎn)出修訂事件報告所沿用的規(guī)則。注意，由于我們把方程組投影到本征坐標(biāo)系下，獨立觀測值的數(shù)量少于總觀測值數(shù)量。所以，估計的定位誤差橢圓必然變大，提供了定位不確定性的更真實的表征。該方法主要的進(jìn)步是每個用于事件定位的投影矩陣只計算一次。

2　驗證測試

為了對比由新定位程序帶來的改進(jìn)，我們同時使用目前的國際地震中心定位程序(構(gòu)成基線)和新定位算法，重定位了國際地震學(xué)與地球內(nèi)部物理學(xué)協(xié)會參考事件表(Bondr and McLaughlin，2009b)中約7 200個GT0-5的地震事件。由于兩個定位程序均使用了ak135模型，我們預(yù)計定位改進(jìn)較小。然而，由于我們考慮了相關(guān)模型誤差結(jié)構(gòu)，我們預(yù)計會在定位覆蓋范圍(置信誤差橢圓是否包含真實震中)上有顯著改進(jìn)。由于在定位程序中對深度解進(jìn)行了測試，并直接使用了深度震相，我們還預(yù)計在深度確定上有所改進(jìn)。在隨后的圖中，點線表示基線(目前國際地震中心定位程序)，實線表示新定位程序的結(jié)果。

2.1爆炸

GT0-5數(shù)據(jù)集中的爆炸包括核爆炸和化學(xué)爆炸。一些化學(xué)爆炸(折射剖面截圖)被具有大方位間隙的很少的臺站記錄到，因此易于產(chǎn)生較大的定位誤差。我們將這些化學(xué)爆炸列入驗證測試中，因為對于定位程序，它們代表了極端情況(稀疏的、不均衡的臺網(wǎng))。圖6展示了基線和新定位算法用于爆炸(化學(xué)爆炸及核爆炸)的定位偏差和定位覆蓋范圍的累積分布?？傮w而言，兩種算法在定位精度上并沒有顯著差異，因為90%的定位都與真實定位相距在17km以內(nèi)。然而，當(dāng)考慮了相關(guān)誤差后，定位覆蓋范圍得到顯著改進(jìn)。注意，如果適應(yīng)于地面真實精度的誤差橢圓包含真實定位，那么由公式(2)描述的定位覆蓋范圍參數(shù)小于1，并且當(dāng)真實定位落在誤差橢圓外時，參數(shù)大于1。因此，實際覆蓋范圍是百分比水平，而覆蓋范圍參數(shù)是1。基線定位程序95%的置信誤差橢圓僅僅包含了當(dāng)時45%左右的真實震中。當(dāng)考慮了相關(guān)誤差后，從新定位程序中獲得的90%的置信誤差橢圓包含了當(dāng)時70%左右的真實震中。

圖7展示了法屬玻利尼西亞和新地島核試驗場的結(jié)果，它們代表了頻譜中的極端情況。在法屬玻利尼西亞的土阿莫土群島和方加陶法環(huán)礁的核爆炸由非常不均衡的遠(yuǎn)震臺網(wǎng)記錄到，其中日本、歐洲和美國加利福尼亞的高密度臺網(wǎng)占主導(dǎo)。因此，法屬玻利尼西亞核試驗場為證明相關(guān)誤差的效果給出了一個教科書式的案例。的確，正如圖7a所示，當(dāng)考慮相關(guān)誤差后，我們減少了定位偏差并實現(xiàn)90%置信的真實覆蓋，雖然基線算法提供了一個極端的25%的真實覆蓋。另一方面，在新地島試驗場(圖7b)進(jìn)行的核爆炸由遠(yuǎn)震臺站和區(qū)域臺站共同記錄到，并且我們的定位和基線算法相比變差。這可以因為一旦臺站的協(xié)同作用被去除，區(qū)域數(shù)據(jù)的影響是增加的，并且簡單的ak135模型對該區(qū)域的區(qū)域走時預(yù)測影響(例如，Krementskayaetal，2001；Hicksetal，2004；Levshinetal，2007)使得解更糟。不過，我們還是在定位覆蓋上獲得了改進(jìn)。

2.2地震：深度不約束和用深度震相約束的深度解答

占國際地震學(xué)與地球內(nèi)部物理學(xué)協(xié)會參考事件表3/4的某些地震有足夠的深度解來得到自由深度解。圖8展示了來自地面真實震源的定位偏差和深度差異以及利用基線和新定位算法的定位覆蓋的累積分布。對于自由深度解，利用新定位算法，我們在震中和深度上都獲得了系統(tǒng)的提高。特別是，我們達(dá)到了90%的真實定位覆蓋。而且，95%的事件與真實深度偏離小于10km。對于可以獲得深度震相深度解(圖9)的地震，對震中和深度的定位提高更為顯著。

圖6　爆炸的定位偏差和位置范圍值的累積分布。點線表示目前國際地震中心定位程序的結(jié)果；實線表示新定位程序的結(jié)果。位置覆蓋范圍由公式(2)計算出。如果位置范圍參數(shù)小于或等于1，那么誤差橢圓包含真實定位。當(dāng)前定位程序95%的置信誤差橢圓的實際覆蓋范圍只有45%左右。當(dāng)考慮了相關(guān)誤差后，新定位程序90%的置信橢圓的實際位置范圍包含了大于70%的真實定位

圖7　(a)法屬玻利尼西亞核爆炸的定位偏差和定位覆蓋范圍的累積分布。當(dāng)考慮了相關(guān)誤差后(實線)，我們在定位和覆蓋范圍上均取得了顯著的改進(jìn)。(b)新地島核試驗場核爆炸的定位偏差和位置范圍值的累積分布。相關(guān)誤差去除了遠(yuǎn)震臺站的協(xié)同作用，但顯示出區(qū)域性ak135模型預(yù)測對這一區(qū)域的不足，這導(dǎo)致定位變差了

2.3地震：固定深度解

正如前面所述，如果有深度解，我們只嘗試得到自由深度解。如果數(shù)據(jù)沒有深度解，我們就把深度固定為區(qū)域依存缺省深度。與將深度固定到常見深度值(35km或10km)不同，我們根據(jù)Bolton等(2006)的方法，并且利用EHB公報(Engdahletal，1998)在全球網(wǎng)格選擇生成缺省深度值，EHB公報的深度值優(yōu)于國際地震中心的深度估計值。我們只考慮EHB自由深度解并產(chǎn)出缺省深度的1°×1°網(wǎng)格。缺省深度被解釋為網(wǎng)格單元中所有深度的中位值，條件是在單元中至少有5個事件，并且75%～25%的四分位距小于100km。后面強加的約束是為了避免這個區(qū)域同時有深的和淺層地震活動性。缺省深度網(wǎng)格沒有給出深度值的定位，我們將深度固定在CRUST2.0(Bassinetal，2000)地殼的中點。

圖10a展示了從EHB自由深度解得到的缺省深度網(wǎng)格。EHB缺省深度解網(wǎng)格很好地表示了主要地震區(qū)域的深度分布。然而，由于EHB事件選取標(biāo)準(zhǔn)，EHB公報只有5級以上地震是完整的，并且缺省深度網(wǎng)格漏掉了像北美洲西部和地中海這樣的地震活動性以相對小的事件為主的重要區(qū)域。

最初，我們用上文所述的缺省深度網(wǎng)格來得到區(qū)域依存的缺省深度值，但我們對固定深度事件的結(jié)果不滿意。所以，受到通過新定位算法獲得的自由深度解的持續(xù)改善的鼓舞，我們著手創(chuàng)建新的缺省深度網(wǎng)格，提供地震活動性很完善的缺省深度。我們用新定位程序重新定位整個《國際地震中心公報》，并同時采用自由深度解和EHB自由深度解，包括標(biāo)記為有可靠深度估計的(Bob Engdahl，私人通訊，2010)固定深度EHB地震，來生成新的缺省深度網(wǎng)格。顯著增加的事件數(shù)量(從90 000到815 000)允許我們把網(wǎng)格細(xì)化為0.5°×0.5°分辨率并且

圖9　定位偏差、相對于地面真實情況的深度差異以及地震的深度震相深度解的位置范圍值的累積分布。新的定位算法(實線)在震源和位置范圍上提供了明顯的改進(jìn)

填補沒有被EHB覆蓋上的空白。圖10b展示了更新的缺省深度網(wǎng)格。對于缺省深度網(wǎng)格沒有提供缺省深度值的位置，現(xiàn)在我們不采用CRUST2.0，為了獲得更可靠的深度缺省估計，而把深度固定為上報深度的中位值。我們對如何獲得每個在《國際地震中心公報》中的事件的深度進(jìn)行描述。

圖11展示了定位偏差、相對于地面真實的深度差異以及來自基線算法(點線)和利用EHB缺省深度網(wǎng)格重新定位算法獲得的(虛線)和更新的缺省深度網(wǎng)格(實線)的固定深度解的位置范圍的累積分布。更新缺省深度網(wǎng)格在有地震活動性的區(qū)域提供了合理的缺省深度。回想目前的國際地震中心定位程序?qū)⑸疃裙潭ㄔ谝粋€機構(gòu)上報的試驗震源上，并且對于記錄良好的地面真實事件，至少有一家上報機構(gòu)得到的深度是對的。雖然我們對固定深度事件沒有獲得定位改進(jìn)，但更新缺省深度網(wǎng)格為地震活躍區(qū)域提供了適當(dāng)?shù)娜笔∩疃取?/p>

2.4成對比較

在前面的章節(jié)里，我們討論了新舊定位程序產(chǎn)出的自由和固定深度解結(jié)果。以另外一種方式看結(jié)果，是用兩種定位程序定位的事件的累積分布。這可以讓我們對改進(jìn)進(jìn)行量化。圖12展示了定位偏差、深度差異和用基線以及新定位算法定位國際地震學(xué)與地球內(nèi)部物理學(xué)協(xié)會參考事件表中1 180個爆炸和6 003個地震的所定位置范圍的累積分布。由于大多數(shù)在國際地震學(xué)與地球內(nèi)部物理學(xué)協(xié)會參考事件表中的GT0-5的事件記錄良好，大多數(shù)定位偏差是不明確的，也即，相對于地面真實的定位偏差在地面真實的精度內(nèi)。然而，更多的事件得到了改善而不是變差了。對于我們致力于改進(jìn)深度的地震和爆炸，在確定位置范圍方面獲得了顯著改進(jìn)。

圖13展示了地震群的更多細(xì)節(jié)。圖13a展示由基線算法得到的5 243個地震的自由深度解的累積分布。圖13b展示由基線算法得到的760個地震的固定深度解的累積分布。注意，由于深度解測試，新的定位程序相對于自由深度事件能夠獲得自由深度解，老的定位程序不能獲得自由深度解，反之亦然。對于這兩種情況，圖13表明，在老的定位程序獲得自由深度解或固定深度解時，我們用新定位程序改進(jìn)了震源參數(shù)和位置范圍的確定。

圖10　(a)由EHB自由深度解得到的1°×1°網(wǎng)格上的缺省深度。(b)由EHB自由深度解獲得的更新缺省深度網(wǎng)格和標(biāo)記為可靠深度的EHB事件，以及由整個國際地震中心公報用新的定位算法重新定位的自由深度解。事件的數(shù)目非常龐大，這允許我們將網(wǎng)格變?yōu)楦?xì)的0.5°×0.5°的分辨率，并允許我們將占主導(dǎo)的如在北美洲和地中海的更小震級事件填寫進(jìn)地震活動區(qū)，而且勾畫出澳洲板塊的南邊界。注意彩色點的大小代表了網(wǎng)格單元的大小(原圖為彩色圖——譯注)

圖11　定位偏差、相對于地面真實的深度差異以及固定深度解的位置范圍值的累積分布。點線代表來自基線算法的解，虛線是通過EHB缺省深度網(wǎng)格獲得的解，實線表示的解利用了更新的缺省深度網(wǎng)格。更新的缺省深度網(wǎng)格在有地震活動性的區(qū)域提供了合理的缺省深度

2.5相關(guān)誤差的影響

如前所述，考慮相關(guān)誤差會產(chǎn)生更多可靠的非確定性估計，并且對于嚴(yán)重不均衡的臺網(wǎng)會降低偏差。所以，我們期待在確定位置范圍上有全面改善，并且對于少數(shù)事件在震源參數(shù)上有改進(jìn)。換句話說，我們期待大部分改進(jìn)在分布的尾部上。

圖14顯示中位值和90%百分比定位偏差、發(fā)震時刻和來自地面真實的深度差異，也顯示了誤差橢圓范圍、位置范圍和作為臺站數(shù)量的函數(shù)，為7 183個用老的和新的定位程序一般定位的事件的誤差橢圓離心率。在中位值水平(細(xì)虛線和實線)，震源參數(shù)幾乎沒有區(qū)別。雖然稀疏的臺網(wǎng)也會遇到相關(guān)模型誤差的問題，但如果相關(guān)誤差結(jié)構(gòu)被忽略，密集的臺網(wǎng)就有足夠數(shù)量的相關(guān)射線路徑可能引入定位偏差。一旦臺站數(shù)超過100個，在90%置信水平(粗虛線和實線)，我們獲得了震中、深度和發(fā)震時刻一致的改進(jìn)。

圖12　定位偏差、相對于地面真實的深度差異以及用基線(點線)和新的定位算法(實線)定位(a)中1 180個爆炸和(b)中6 003個地震的位置范圍的累積分布。更多的事件得到了改進(jìn)而不是變差了

在定位算法中通過估計全部數(shù)據(jù)的協(xié)方差矩陣并結(jié)合相關(guān)誤差結(jié)構(gòu)的主要優(yōu)點之一是：從相關(guān)觀測的數(shù)量中分離出90%的誤差橢圓。的確，當(dāng)相關(guān)誤差被并入定位算法中時，超過100個臺站的誤差橢圓范圍在中位值與90%百分比水平是持平的，因此保持了位置范圍的覆蓋。我們達(dá)到了90%的覆蓋，第90百分位曲線會離開1。由于我們只有約80%的覆蓋，曲線僅僅超過1以上。另一方面，忽略相關(guān)誤差時，誤差橢圓范圍伴隨臺站數(shù)目無限收縮，所以失去了覆蓋?；叵胍幌潞篁?zāi)Ｐ蛥f(xié)方差矩陣依賴于獨立觀測數(shù)量的情況。因為相等的非相關(guān)數(shù)據(jù)(殘差的線性組合)是減少的，由后驗?zāi)Ｐ蛥f(xié)方差矩陣描述的形式定位不確定性變大，導(dǎo)致擴大和更圓的以誤差橢圓離心率(離心率越小，誤差橢圓越圓)指示的誤差橢圓。

2.6與EHB的比較

圖15展示了定位偏差、發(fā)震時刻、對在地面真實數(shù)據(jù)集中的EHB事件(虛線)對于基線(點線)和新定位算法(實線)的相對于地面真實的深度差異的累積分布。圖中表明用新的定位程序符合或超過了EHB定位精度，但不符合EHB的深度精度。

圖13　定位偏差、相對于地面真實的深度差異以及用基線(點線)定位算法定位(a)5 243個自由深度解和(b)760個固定深度解地震的定位覆蓋的累積分布。新的定位程序(實線)改進(jìn)了更多事件

圖14　中位值(細(xì)灰線)和90%百分比(粗黑線)定位偏差、發(fā)震時刻和深度差異、誤差橢圓范圍、定位覆蓋和采用老的(點線)和新的(實線)定位方法并以臺站數(shù)量為函數(shù)的誤差橢圓離心率。當(dāng)我們考慮了相關(guān)誤差時，誤差橢圓的大小趨于平穩(wěn)，提供了更好的位置覆蓋，并且一旦包含臺網(wǎng)的所有信息，震源偏差就減少。因為(本征的)臺網(wǎng)更均衡了，誤差橢圓變得更圓了

2.7臺站和臺網(wǎng)震級

體波或面波臺站震級在老的定位程序中是用同樣的公式計算出來的。大多數(shù)上報機構(gòu)——采用國際地震學(xué)與地球內(nèi)部物理學(xué)協(xié)會推薦的測量振幅和周期，國際地震中心就可以根據(jù)國際地震學(xué)與地球內(nèi)部物理學(xué)協(xié)會的標(biāo)準(zhǔn)(Bormannetal，2009)計算一系列震級了。

對于上報的每個震級，我們從上報的P波周期≤3s，震中距在21°～100°之間的震級利用Gutenberg和Richter(1956)的公式計算mb。我們在周期10～60s，震源深度≤60km，震中距20°～160°的事件用Prague公式(Vaneketal，1962)測定面波震級MS。我們設(shè)臺站震級為在上報臺站的振幅和周期對中最大A/T值的那個臺站的震級。為避免離群臺站震級導(dǎo)致的錯誤，現(xiàn)在我們至少需要三個有效臺站震級并且按α修削中位值來計算臺網(wǎng)震級。臺網(wǎng)震級定為去除了低于和高于20%的臺站震級(α修削過程)后的臺站震級中位值，并且它的不確定性由標(biāo)準(zhǔn)中位值絕對偏差圍繞中位值給出。由于新程序至少要求三個臺站震級來產(chǎn)出臺網(wǎng)震級，這減緩了《國際地震中心公報》中由于引入離群振幅觀測導(dǎo)致顯而易見的震級誤差的風(fēng)險。

3　全部《國際地震中心公報》的自動重定位

如上所述，我們重新定位了全部《國際地震中心公報》(1960～2010年)，包括4年的《國際地震摘要》(ISS，國際地震中心的前身)目錄(Villaseor and Engdahl，2005，2007)，目的是獲得更新后的來源于地震活動性本身的缺省深度網(wǎng)格。重定位《國際地震中心公報》時，我們盡可能地重新模仿國際地震中心自動定位程序(大多數(shù)年份的公報已經(jīng)經(jīng)過國際地震中心分析員修訂)，也即，不管是否存在國際地震中心的結(jié)果，我們都忽略了國際地震中心、EHB和地面真實的震源，并且定位了公報中的所有事件。

圖16a展示了用于事件定位的時間累積分布，在這次應(yīng)用中重定位了2 507 492個事件。在標(biāo)準(zhǔn)雙核Linux機器上，96%的事件在一秒鐘內(nèi)被定位。這些時間包括鄰域算法網(wǎng)格搜索、數(shù)據(jù)協(xié)方差矩陣的反演、線性迭代定位算法和震級計算。圖16b展示了獲得解的平均運行時間是相關(guān)震相數(shù)目的函數(shù)。不出所料，執(zhí)行時間隨觀測數(shù)目的增加呈指數(shù)級增加，并且當(dāng)觀測數(shù)量超過1 500個時，事件定位的平均用時超過了一分鐘。這對于自動定位是可接受的運行。對于國際地震中心的分析員來說，修訂能力可進(jìn)一步提高，因為在那個階段，由于自動定位已提供給編輯人員，所以不必執(zhí)行操作中最慢的一步，即鄰域算法搜索。

最后，我們比較了用新的定位程序獲得的自動定位和修訂的《國際地震中心公報》。

圖15　定位偏差、相對于EHB地震發(fā)震時刻和相對于地面真實的深度差異的累積分布。點線代表用基線定位算法得到的解，虛線是EHB的解，實線表示新定位程序的解。從新定位程序獲得的解符合或超過了EHB定位準(zhǔn)確度，但不符合EHB的深度準(zhǔn)確度

當(dāng)前修訂的《國際地震中心公報》(1960～2008-08)包含約1百萬個地震。新定位程序的自動定位提供了約2.3百萬個地震的震源(我們重定位了沒有被國際地震中心修訂的事件)。修訂的《國際地震中心公報》和重定位目錄中共有的地震數(shù)目為993 685個。下面我們僅考慮一般定位地震的數(shù)據(jù)集。

我們期望通過對后續(xù)震相的使用和對深度解的檢驗，來減少事件定位的分散。地震活動的緊密性被習(xí)慣地用作定位改進(jìn)的定性估計。為了檢驗我們的假設(shè)，我們通過計算震中的熵來量化地震活動性的分散。我們采用Nicholson等(2000)的定義來計算點集的熵。他們的定義使用了沃羅努瓦圖(Shewchuk，1996)。沃羅努瓦單元的面積與事件的密度成比例，并由此可被用于估計空間中地震分布的概率密度函數(shù)。沃羅努瓦圖(德洛奈三角測量的幾何對偶)對應(yīng)一個獨特的區(qū)域劃分，每個單元包含一個事件，并且單元中的任意點離這個事件要比離其他事件近。因此，沃羅努瓦單元定義了自然鄰域——單元越小，事件越密集。如公式(12)，我們定義了點集的熵：

(12)

式中，N是事件的數(shù)目，ai是沃羅努瓦單元的面積，A是凸包的面積(包含點集的最小凸多邊形)。S是非正的，并且當(dāng)事件分布均勻時達(dá)到它的最大值0。因此，如果相對于從修訂的《國際地震中心公報》中獲得的熵減少，我們就接受我們的假設(shè)(也即，新的定位程序使地震活動性變緊密)。為確保正確，我們應(yīng)該測量來自震源點群的三維沃羅努瓦單元體積，然而在球面上那樣做是微不足道的。因此我們采取在地球表面測量二維熵的近似值的方式。

我們計算了修訂的和重定位的《國際地震中心公報》中每個地震區(qū)域震中的熵。圖17a通過地震區(qū)域號碼展示了從修訂的國際地震中心公報中和從使用新定位程序重定位的公報中分別獲得的每個地震區(qū)域里震中的熵。圖17b展示了重定位公報和《國際地震中心公報》地震中熵的差別。50個地震區(qū)域中，有44個新定位程序降低了數(shù)據(jù)集的熵，有些情況下還很顯著，并且總體上，平均超過了-0.5。因此，新定位程序提供了更好的事件群集，使得地震活動性更加緊密。圖18a～g進(jìn)一步說明了這一點。我們標(biāo)繪了來自修訂《國際地震中心公報》(左側(cè)展板)和新定位程序獲得的(右側(cè)展板)不同地震活動區(qū)域事件的震源。即使在相對大比例的地圖中，新定位程序獲得的地震活動的緊密性也很明顯。

4　討論和結(jié)論

我們提出的國際地震中心新定位程序利用了所有具有有效的ak135模型預(yù)測的上報震相，通過鄰域算法獲得初始震源估計，并且在線性迭代最小二乘定位算法中考慮相關(guān)走時預(yù)測誤差結(jié)構(gòu)。

我們已經(jīng)指出，如果在由非對角數(shù)據(jù)協(xié)方差矩陣定義的本征坐標(biāo)系下解方程，就可以考慮相關(guān)模型誤差結(jié)構(gòu)。根據(jù)Bondr和McLaughlin(2009a)，我們用通用變差函數(shù)模型為遠(yuǎn)震P波構(gòu)建了數(shù)據(jù)協(xié)方差矩陣的先驗估計。使用變差函數(shù)的基本假設(shè)是：射線路徑間的相似性可以用兩個臺站之間的距離近似，并且各向同性的變差函數(shù)充分地刻畫了相關(guān)結(jié)構(gòu)。基于這些簡化假設(shè)，數(shù)據(jù)協(xié)方差矩陣和它的逆只需被計算一次。

然而，各向同性假設(shè)忽視了三維速度結(jié)構(gòu)，并且在有更小的平均多種長度尺度的速度擾動的更短射線路徑下變得不符合。不可否認(rèn)，由遠(yuǎn)震震相得到的通用變差函數(shù)模型不能充分捕捉到地方和區(qū)域距離下的所有三維速度變化的細(xì)節(jié)。所以，對于只有地方或近區(qū)域記錄到的較小事件，我們的方法可能會低估或高估相似路徑下走時預(yù)測之間的關(guān)聯(lián)。有人可能嘗試導(dǎo)出區(qū)域和地方距離的變差函數(shù)模型，但那必須是特定區(qū)域的且只提供不好的三維速度模型替代。此外，回想一下由變差函數(shù)模型得到的先驗測量誤差方差加到數(shù)據(jù)協(xié)方差矩陣的對角線。大的測量誤差會支配數(shù)據(jù)協(xié)方差矩陣，并傾向于削弱相關(guān)結(jié)構(gòu)的走時預(yù)測誤差，也即，它們減小了非對角元素間的相關(guān)長度。增加背景值可以平衡這個影響。因此，我們引入近3倍的增

圖16　(a)用于定位事件的時間累積分布。在標(biāo)準(zhǔn)雙核Linux機器上，96%的事件在一秒鐘內(nèi)被定位。(b)事件定位的平均運行時間是有關(guān)震相數(shù)量的函數(shù)

圖17　(a)從修訂的《國際地震中心公報》(灰色)和用新定位程序重定位(黑色)獲得的每個地震區(qū)震中的熵。(b)熵之間的差別。對大多數(shù)區(qū)域，新定位程序降低了數(shù)據(jù)集的熵，因此使得地震活動分布更緊密

加背景值是考慮到在地方和區(qū)域距離下沒有解釋的模型誤差，更重要的是：強制考慮相關(guān)誤差結(jié)構(gòu)忽視了相當(dāng)大的對角元素。顯而易見，全球的三維速度模型能真正解決這個問題。的確，完美的走時預(yù)測將不會有模型誤差，且非常相似的射線路徑的走時預(yù)測可能保持不相關(guān)。不過，我們說明了即使用擬合遠(yuǎn)震觀測的一維模型和變差函數(shù)模型也能獲得像90%置信誤差橢圓覆蓋真實位置80%～85%的現(xiàn)實的非確定性估計。

前面已提到，我們對地面真實事件說明的定位改進(jìn)是一致的，但較小。這并不奇怪，因為國際地震學(xué)與地球內(nèi)部物理學(xué)協(xié)會參考事件表中的大多數(shù)事件是以小的方位間隙和P波類型震相占主導(dǎo)記錄的。在這些情況下，我們期望定位有顯著改善，只針對大量的相關(guān)射線路徑共同引入定位偏差的嚴(yán)重不均衡的臺網(wǎng)。另一方面，《國際地震中心公報》代表了過多的從合理到最不利的臺網(wǎng)幾何結(jié)構(gòu)的臺站配置。因此，我們可以期望在重定位整個《國際地震中心公報》時有更驚人的定位改進(jìn)。雖然在這種情況下由于缺乏地面真實事件信息我們不能衡量定位精度的改進(jìn)，但事件定位更好的聚集表明新的定位程序的確勝過基線算法。

圖18　修訂的《國際地震中心公報》(老定位程序，左側(cè))和重定位的《國際地震中心公報》(新定位程序，右側(cè))不同區(qū)域一般地震活動圖的比較。(a)北安第斯山。(b)南安第斯山。(c)中美洲。(d)阿拉斯加。(e)關(guān)島—本州—琉球群島。(f)興都庫什—帕米爾高原。(g)安的列斯群島南部。當(dāng)使用新的定位程序重定位后，事件變得更加緊密(原圖為彩色圖——譯注)

這些改進(jìn)歸因于新定位算法的幾個功能的組合。它們是：鄰域算法網(wǎng)格搜索為開始線性反演提供了很好的初始震源；考慮相關(guān)走時預(yù)測誤差提供了可靠的形式不確定性估計和減少不均衡的臺網(wǎng)定位偏差；在定位算法中使用深度敏感震相得到相當(dāng)好的深度估計；對于沒有深度解的事件，缺省深度網(wǎng)格給出了合理的缺省深度估計?？坍嫷诙秃罄m(xù)震相的影響更為困難，盡管它們攜帶了非常寶貴的關(guān)于地震震源的信息，但只能謹(jǐn)慎對待它們。舉例來說，定位算法中利用遠(yuǎn)區(qū)域震相時，由于三維非均勻性在一維速度模型中沒有被描述出來，可能使定位更糟。這就是我們?yōu)槭裁唇o后續(xù)震相分配較大先驗測量誤差的原因。

只有較好的走時預(yù)測能提供顯著的定位改進(jìn)，比如全球三維速度模型。盡管過去的幾年我們已經(jīng)看到在研究三維速度模型、給出走時精確預(yù)測方面(例如，Myers and Schultz，2000；Ritzwolleretal，2003；Yangetal，2004；Morozovetal，2005；Murphyetal，2005；Reiteretal，2005；Flanaganetal，2007；Zuccaetal，2009；Myersetal，2010)取得了相當(dāng)大的進(jìn)步，但它們還沒能被應(yīng)用到日常工作中。而且，給國際地震中心提供海量的震相報告的時候，用目前國家最先進(jìn)的射線追蹤方法來追蹤通過三維模型的射線將是昂貴的。因此，我們有更適中的目標(biāo)：用全球一維模型獲得可能最好的震源和相應(yīng)的不確定性估計。

我們計劃使用新的國際地震中心定位程序重建整個《國際地震中心公報》，目的是提供地球地震活動性質(zhì)量一致的公報。注意，當(dāng)前的《國際地震中心公報》不是完全同質(zhì)量的。例如，國際地震中心剛開始在定位中僅使用了初至P波震相，直到2001年才開始在定位程序中引入了初至S波震相；2006年才引入ak135模型；分別于1999年、2005年和2006年決定選擇哪些報告事件用于國際地震中心修訂的程序。此外，國際地震中心僅從1978年才開始發(fā)布MS震級，盡管自1971年開始面波振幅就由各個機構(gòu)定期上報。根據(jù)全球地震風(fēng)險模型(www.globalquakemodel.org)計劃，當(dāng)前從歷史臺站公報中收集面波振幅讀數(shù)的努力，將會進(jìn)一步提高震級評估的連續(xù)性和一致性。因此，重建《國際地震中心公報》將使我們有機會去校正公報中已知的矛盾，去掉假事件和錯誤；為個別具有錯誤時間戳的臺站鑒定和標(biāo)記時段；并且給人為事件重新分配事件類型標(biāo)志。

在重建過程中，我們將引入基礎(chǔ)數(shù)據(jù)集，這在最初的《國際地震中心公報》出版時是沒有的，例如，覆蓋了1960～1963年的《國際地震摘要》數(shù)據(jù)，還有來自固定臺網(wǎng)因遲報而未被包括在國際地震中心修訂中的數(shù)據(jù)，或者不慎沒有被國際地震中心使用的數(shù)據(jù)。我們還將引入在政局動蕩或行政糾紛之后恢復(fù)回顧性的來自固定臺網(wǎng)的數(shù)據(jù)。如果有的話，我們將從包括海底地震計的臨時布署中收集震相拾取。我們預(yù)計這些努力將需要數(shù)年才能完成。

我們已證明新的國際地震中心定位算法提供了小而一致的定位改進(jìn)、相當(dāng)大的深度測定改進(jìn)和顯著的更加精確的形式不確定性估計。在有地震活動的地方，缺省的深度網(wǎng)格提供了合理的深度估計。我們已經(jīng)說明由新算法獲得的定位和深度精度配得上或已超過了EHB的精度。通過后續(xù)震相和深度解測試，以及十分聚集的事件定位，我們也證明了新算法為地球地震活動性提供了改進(jìn)的觀點。

數(shù)據(jù)與來源

本文所有數(shù)據(jù)來自參考文獻(xiàn)中列出的出版資料。國際地震學(xué)與地球內(nèi)部物理學(xué)協(xié)會參考事件表、EHB公報和《國際地震中心公報》均由國際地震中心網(wǎng)站主辦(www.isc.ac.uk)。國際地震中心目前使用的定位程序可以從國際地震中心網(wǎng)站上下載，新的定位程序也即將公布。我們使用了de Hoon等(2004)開發(fā)的開源Cluster3.0庫和Wessel與Smith(1991)開發(fā)的通用繪圖工具(GMT)。

Adams，R.D.，Hughes，A.A.& Gregor，D.M.，1982.Analysis procedures at the International Data Centre，Phys.Earthplanet.Int.，30，85-92.

Amante，C.& Eakins，B.W.，2009.ETOPO1 1 arc-minute global relief model：procedures，data sources and analysis，NOAATechnicalMemorandumNESDISNGDC-24.

Anderson，K.R.，1982.Robust earthquake location usingMestimates，Phys.Earthplanet.Int.，30，119-130.

Bassin，C.，Laske，G.& Masters，G.，2000.The current limits of resolution for surfacewave tomography in North America，EOS，Trans.Am.geophys.Un.，81，F(xiàn)897.

Bolt，B.A.，1960.The revision of earthquake epicentres，focal depths and origin time using a high-speed computer，Geophys.J.R.astr.Soc.，3，434-440.

Bolton，M.K.，Storchak，D.A.& Harris，J.，2006.Updating default depth in the ISC bulletin，Phys.Earthplanet.Inter.，158，27-45.

Bormann，P.，Liu，R.，Xu，Z.，Ren，K.，Zhang，L.& Wendt，S.2009.First application of the new IASPEI teleseismic magnitude standards to data of the China National Seismographic Network，Bull.seism.Soc.Am.，99，1868-1891.

Buland，R.& Chapman，C.H.，1983.The computation of seismic travel times，Bull.seism.Soc.Am.，73，1271-1302.

Buland，R.，1986.Uniform reduction error analysis，Bull.seism.Soc.Am.，76，217-230.

Chang，A.C.，Shumway，R.H.，Blandford，R.R.& Barker，B.W.，1983.Two methods to improve location estimates—preliminary results，Bull.seism.Soc.Am.，73，281-295.Dziewonski，A.M.& Anderson，D.L.，1981.Prelimi-nary reference earth model，Phys.Earthplanet.Inter.，25，297-356.

Dziewonski，A.M.& Gilbert，F(xiàn).，1976.The effect of small，aspherical perturbations on travel times and a re-examination of the correction for ellipticity，Geophys.J.R.astr.Soc.，44，7-17.

Engdahl，E.R.，van der Hilst，R.& Buland，R.，1998.Global teleseismic earthquake relocation with improved travel times and procedures for depth determination，Bull.seism.Soc.Am.，88，722-743.

Evernden，J.，1969.Precision of epicenters obtained by small numbers of world-wide stations，Bull.seism.Soc.Am.，59，1365-1398.

Flanagan，M.P.，Myers，S.C.& Koper，K.D.，2007.Regional travel-time uncertainty and seismic location improvement using a three-dimensional a priori velocity model，Bull.seism.Soc.Am.，97，804-825.

Gutenberg，B.& Richter，C.F.，1956.Magnitude and energy of earthquakes，Ann.Geof.，9，1-15.

Hicks，E.C.，Kvaerna，T.，Mykkeltveit，S.，Schweitzer，J.&Ringdal，F(xiàn).，2004.Travel-times and attenuation relations for regional phases in the Barents Sea region，Pureappl.Geophys.，161，1-19.

de Hoon，M.J.L.，Imoto，S.，Nolan，J.& Miyano，S.，2004.Open source clustering software，Bioinformatics，20，1453-1454.

Jeffreys，H.，1932.An alternative to the rejection of observations，Proc.Roy.Soc.Lond.，187，78-87.

Jeffreys，H.& Bullen，K.E.，1940.Seismological Tables，British Association for the Advancement of Science，Gray-Milne Trust，London.Jordan，T.H.& Sverdrup，K.A.，1981.Teleseismic location techniques and their application to earthquake clusters in the South-central Pacific，Bull.seism.Soc.Am.，71，1105-1130.

Kennett，B.L.N.，2006.Non-linear methods for event location in a global context，Phys.Earthplanet.Int.，158，45-64.

Kennett，B.L.N.& Gudmundsson，O.，1996.Ellipticity corrections for seismic phases，Geophys.J.Int.，127，40-48.

Kennett，B.L.N.，Engdahl，E.R.，& Buland，R.，1995.Constraints on seismic velocities in the Earth from traveltimes，Geophys.J.Int.，122，108-124.

Kremenetskaya，E.，Asming，V.& Ringdal，F(xiàn).，2001.Seismic location calibration of the Euro-pean Arctic，Pureappl.Geophys.，158，117-128.

Levshin，A.L.，Schweitzer，J.，Weidle，C.，Shapiro，N.M.&Ritzwoller，M.H.，2007.Surface wave tomography of the Barents Sea and surrounding regions，Geophys.J.Int.，170，441-459.

Morozov，I.B.，Morozova，E.A.，Smithson，S.B.，Richards，P.G.，Khalturin，V.I.& Solodilov，L.N.，2005.3D first-arrival regional calibration model of Northern Eurasia，Bull.seism.Soc.Am.，95，951-964.

Murphy，J.R.& Barker，B.W.，2006.Improved focal-depth determination through automated identification of the seismic depth phases pP and sP，Bull.seism.Soc.Am.，96，1213-1229.

Murphy，J.R.，et al.，2005.Calibration of International Monitoring System(IMS)stations in Eastern Asia for improved seismic event location，Bull.seism.Soc.Am.，95，1535-1560.

Myers，S.C.& Schultz，C.A.，2000.Improving sparse network seismic location with Bayesian kriging and teleseismically constrained calibration events，Bull.seism.Soc.Am.，90，199-211.

Myers，S.C.，et al.，2010.A crust and upper mantle model for Eurasia and North Africa for Pn travel-time calculation，Bull.seism.Soc.Am.，100，640-656.

Nicholson T.，Sambridge，M.，& Gudmundsson，O.，2000.On entropy and clustering in earthquake hypocenter distributions，Geophys.J.Int.，142，37-51.Reiter，D.，Rodi，W.，& Johnson，M.，2005.Develop-ment of a tomographic upper-mantle velocity model beneath Pakistan and Northern India，Bull.seism.Soc.Am.，95，926-940.

Ritzwoller，M.H.，Shapiro，N.M.，Levshin，E.A.，Bergman，E.A.，& Engdahl，E.R.，2003.Ability of a global three-dimensional model to locate regional events，J.geophys.Res.，108(B7)，2353，doi：10.1029/2002JB002167.

Sambridge，M.，1999.Geophysical inversion with a neighbourhood algorithm.I.Searching the parameter space，Geophys.J.Int.，138，479-494.

Sambridge，M.& Kennett，B.L.N.，2001.Seismic event location：non-linear inversion using a neighbourhood algorithm，Pureappl.Geophys.，158，241-257.

Schweitzer，J.& Storchak，D.A.，2006.Modernizing the ISC location procedures，Editorial，Phys.Earthplanet.Inter.，158，1-3.

Shewchuk，J.R.，1996.Engineering a 2D quality mesh generator and Delaunay triangulator，F(xiàn)irstWorkshoponComputationalGeometry，Philadelphia，ACM，pp.124-133.

Storchak，D.A.，Chen，Q.F.，Willemann，R.J.，& Andrianirina，M.，2000.Improved locations for moderately large earthquakes using regional S and PKP.EOS，Trans.Am.geophys.Un.，81(48)，S71C-01.

Storchak，D.A.，Schweitzer，J.，& Bormann，P.，

2003.The IASPEI standard seismic phase list，Seism.Res.Lett.，74，761-772.

Vanek，J.，Zatopek，A.，Karnik，V.，Kondorskaya，N.V.，Riznichenko，Y.V.，Savarensky，Y.F.，Solovev，S.L.& Shebalin，N.V.，1962.Standardization of magnitude scales，Bull.(Izvest.)Acad.Sci.USSR，Geophys.Ser.，2，108-111.

Wessel，P.& Smith，W.H.F.，1991.Free software helps map and display data，EOS，Trans.Am.geophys.Un.，72，441.

Willemann，R.J.& Storchak，D.A.，2001.Data collection at the International Seismological Centre，Seism.Res.Lett.，72，440-453.

Yang，X.，et al.，2004.Validation of regional and teleseismic traveltime models by relocating GT events，Bull.seism.Soc.Am.，94，897-919.

Young，J.B.，Presgrave，B.W.，Aichele，H.，Wiens，D.A.& Flinn，E.A.，1996.The Flinn-Engdahl regionalization scheme：the 1995 revision，Phys.Earthplanet.Int.，96，223-297.

Zucca，J.J.，et al.，2009.The prospect of using three-dimensional Earth models to improve nuclear explosion monitoring and ground-motion hazard assessment，Seism.Res.Lett.，80，31-39.

張瑩瑩(1986—)，女，中國地震臺網(wǎng)中心助理工程師，固體地球物理學(xué)專業(yè)碩士研究生，主要從事地震監(jiān)測、地殼速度結(jié)構(gòu)和各向異性方面的研究。E-mail：zyycugb1986@163.com。

張瑩瑩 譯.2016.國際地震中心改進(jìn)的定位程序.世界地震譯叢.47(5)：376-400.doi：10.16738/j.cnki.issn.1003-3238.201605002

中國地震臺網(wǎng)中心張瑩瑩譯

中國地震局地球物理研究所吳何珍校

為了應(yīng)對這一挑戰(zhàn)，我們?yōu)閲H地震中心開發(fā)了一種新的定位算法，這種算法考慮了相關(guān)誤差結(jié)構(gòu)。我們采用國際地震學(xué)與地球內(nèi)部物理學(xué)協(xié)會(IASPEI)所有的在ak135模型中有有效走時預(yù)測的標(biāo)準(zhǔn)震相來得到更精確的事件定位。本文中，我們闡述了新的國際地震中心定位方法，為檢驗此方法的可靠性，我們用該方法重新定位在國際地震學(xué)與地球內(nèi)部物理學(xué)協(xié)會參考事件表中的地面真實事件，我們也對整個《國際地震中心公報》做了重新定位。

我們的定位算法雖然只有很小的改進(jìn)，但是對定位結(jié)果的改進(jìn)卻很顯著，尤其是在深度的確定上和更精確的形式不確定性分析上。我們論證了新的定位算法；通過后續(xù)震相的使用和深度解的測試，相當(dāng)多的地震震群變得更緊密了；因此該算法為地球地震活動性提供了一個改進(jìn)的觀點。