Jie ZhangHaijiang ZhangEnhong ChenYi ZhengWenhuan KuangXiong Zhang

1）Laboratory of Seismology and Physics of Earth’s Interior，School of Earth and Space Sciences，University of Science and Technology of China，Hefei 230026，China

2）School of Computer Science and Technology，University of Science and Technology of China，Hefei 230027，China

?

搜索引擎方法在實時地震監(jiān)測中的應用＊

Jie Zhang1）Haijiang Zhang1）Enhong Chen2）Yi Zheng2）Wenhuan Kuang1）Xiong Zhang1）

1）Laboratory of Seismology and Physics of Earth’s Interior，School of Earth and Space Sciences，University of Science and Technology of China，Hefei 230026，China

2）School of Computer Science and Technology，University of Science and Technology of China，Hefei 230027，China

摘 要當?shù)卣鸢l(fā)生時，地震學家需要利用記錄的地震圖盡快去推斷震中、震級和震源機制解。如果能快速判定這些信息，就能進行及時疏散，并采取應急措施以減輕地震災害?，F(xiàn)在有一些先進的方法能在地震后幾秒內(nèi)判斷出最初的震中和震級，但可能需要數(shù)分鐘到幾小時才能估算出震源機制解。這里我們提出針對大型地震圖數(shù)據(jù)庫，利用一種類似于互聯(lián)網(wǎng)搜索引擎的地震搜索引擎方法，來尋找輸入地震波數(shù)據(jù)的最佳匹配。該方法比當前的精確查找法快幾千倍。以2012年3月8日新疆MW5.9地震為例，該搜索引擎方法能在長周期地震面波到達后1s內(nèi)推算出地震震源機制參數(shù)。

關鍵詞搜索引擎方法；實時地震監(jiān)測；震源機制解

引言

在地震學界，實時地震信息發(fā)布早已是一項重要的研究工作［1－3］。近幾年，地震預警系統(tǒng)的發(fā)展進一步推動了此項工作的研究步伐，其能在地震發(fā)生后數(shù)秒到1min左右向公眾發(fā)布預警信息［4］。全世界已有一些地震預警系統(tǒng)正在實施應用中，包括日本的緊急地震速報系統(tǒng)（REIS）［5］，墨西哥的地震預警系統(tǒng)（SAS）［6］，中國臺灣的地震預警系統(tǒng)（VSN）［7］和土耳其的地震預警系統(tǒng)（IERREWS）［8］。在美國加州，開發(fā)的ElarmS預警系統(tǒng)已被列入一項重要的研究計劃，該系統(tǒng)已進行了離線測試，但尚未完全部署［4］。地震學家已研制出可靠的算法用于自動估算地震的震源信息［9－10］。例如，日本REIS系統(tǒng)通過利用來自密集監(jiān)測網(wǎng)的數(shù)據(jù)，在P波到達后5s內(nèi)推算出震中和震級［5］。然而，即使是利用最近發(fā)表的新方法，在地震可能發(fā)生的位置網(wǎng)格點上提前計算好格林函數(shù)，然后利用格林函數(shù)反演矩張量，仍然需要幾分鐘甚至更長時間才能獲取震源機制解［11－13］。利用GPS數(shù)據(jù)進行相似的研究工作也只能在幾分鐘內(nèi)確定大地震的矩張量［14－15］。

除了震中和震級外，實時獲取震源機制解也同樣重要。例如，海嘯預報需要完整的震源參數(shù)，包括震源深度、震級和斷層滑動走向等［16］。2010年10月5日發(fā)生在蘇門答臘島西海岸的MW7.7淺源地震，導致當?shù)匕l(fā)生了海嘯，浪高達到3m，幾分鐘內(nèi)影響了整個島嶼。據(jù)報道，死亡人數(shù)超過400人［17］。震源機制解研究表明，這次地震屬于逆斷層運動，從而造成海水運動。相比較，2012年4月11日在印度洋發(fā)生了MW8.6大地震和隨后的MW8.2余震，盡管在印度洋海域發(fā)布了海嘯預警，但最終并未形成海嘯。通過震源機制解研究表明，這兩次地震都是順著走向滑動造成的［18］，因此垂直向運動位移相對較小，不足以形成海嘯。震源機制解的實時估算對于監(jiān)測斷層活動也非常重要。例如，Bouchon等［19］對1999年伊茲米特MW7.6地震進行了研究，發(fā)現(xiàn)其一系列前震的震源機制解表明在斷層加速運動至破裂過程中存在相似的斷層滑動。實時獲取震群的震源機制解可能有助于識別斷層活動特征。這類信息能讓我們快速發(fā)現(xiàn)地震活動區(qū)并實時跟蹤斷層運動。

當前的挑戰(zhàn)為，當一些臺站收到地震數(shù)據(jù)后幾秒內(nèi)必須自動快速估算出震源機制解。我們開發(fā)了一種基于圖像的地震搜索引擎，類似于互聯(lián)網(wǎng)搜索引擎，在一個單一的AMD工作站的大型數(shù)據(jù)庫中搜索相似的地震波形圖，在1s內(nèi)估算出地震震源機制參數(shù)。計算機搜索技術的重要發(fā)展有助于促進文字、圖像、視頻和來自網(wǎng)絡尺度數(shù)據(jù)集的音頻等搜索行業(yè)的發(fā)展［20－32］。與錄音和一維圖像相似，一幅地震波形圖可認為是地震臺根據(jù)時間記錄的地面運動情況，其包含的信息包括震源和地下介質(zhì)。通過假設已知地球速度模型，我們采用正演模擬方法，在離散的網(wǎng)格點上為可能的地震震源機制和震中位置構建波形數(shù)據(jù)庫。我們的目標是，對于任何一個新的地震記錄，都能從數(shù)據(jù)庫中找出最佳匹配。這個方法是完全自動的，不需要輸入任何參數(shù)或人工干預。因此，它可以用于通常的地震參數(shù)報告。

我們利用在新疆試驗區(qū)中的3個地震實例測試我們的地震搜索引擎法。我們選擇了一個經(jīng)緯度5°×5°的區(qū)域構建一個快速搜索數(shù)據(jù)庫。3個地震事件發(fā)生后，通過搜索引擎法在收到長周期面波數(shù)據(jù)后1s內(nèi)，都可以得到搜索結果。在數(shù)據(jù)庫覆蓋范圍內(nèi)的任意單個地震事件，都可以得到足夠精確的震源參數(shù)。如果單個地震事件發(fā)生在數(shù)據(jù)庫覆蓋范圍以外，或者是多個事件相繼發(fā)生在鄰近區(qū)域，搜索引擎能自動通過預設的互相關閾值排除無效的結果。

1　計算結果

1.1 用于新疆試驗區(qū)的搜索引擎數(shù)據(jù)庫

圖1顯示在南疆一個5°×5°的試驗區(qū)，3個固定地震臺（暗紅色三角）能記錄5°～15°范圍內(nèi)的地震事件。這些臺站包括MAKZ （46.8°N，82.0°E，哈薩克斯坦馬坎奇），KBL（34.5°N，69.0°E，阿富汗喀布爾）和LSA（29.7°N，91.1°E，中國拉薩）。從2000 年1月1日至今，根據(jù)美國地質(zhì)調(diào)查局（USGS）國家地震信息中心（NEIC）的記錄，在試驗區(qū)內(nèi)共發(fā)生MW4.0以上地震51次。

圖1　地震臺站和虛擬試驗區(qū)的位置。

在3個以暗紅色三角形所標示的地震臺站區(qū)域范圍內(nèi)（MAKZ：46.8°N，82.0°E，哈薩克斯坦馬坎奇；KBL：34.5°N，69.0°E，阿富汗喀布爾；LSA：29.7°N，91.1°E，中國拉薩），由虛擬震源構建大型合成數(shù)據(jù)集，創(chuàng)建一個快速搜索數(shù)據(jù)庫。3次地震由星形標示（其定位信息均由美國地質(zhì)調(diào)查局的國家地震信息中心提供），用于測試不同的狀況。紫色星代表2012年3月8日發(fā)生的MW5.9地震；藍色星代表2011年9月15日發(fā)生的MW5.3地震；綠色星代表2014 年4月30日發(fā)生的MW5.3地震在本研究中，我們著重于對研究區(qū)內(nèi)小震和中震的實時分析，以確定它們的震源和雙力偶震源機制解。我們的方法包括以下幾步：①構建格林函數(shù)綜合數(shù)據(jù)庫，在橫向均勻介質(zhì)中設置離散網(wǎng)格節(jié)點，計算每個節(jié)點作為虛擬震源的格林函數(shù)，計算時考慮所有可能的雙力偶源解；②在數(shù)據(jù)庫中利用主成分分析法（principal components analysis，PCA）減少時間采樣點個數(shù)；③利用互聯(lián)網(wǎng)工業(yè)發(fā)展的快速搜索算法，快速找出最佳匹配的虛擬源；④對震源機制解的確認和搜索結果分辨率的定量化。

為了設計一個用于監(jiān)測試驗區(qū)的地震搜索引擎，我們構建了一個包含大量合成地震記錄的數(shù)據(jù)庫，這些記錄與試驗區(qū)內(nèi)3D網(wǎng)格上的每個虛擬源相一致。我們的三分量地震圖計算采用Thompson－Haskell傳播矩陣技術在多層半空間的地震點源彈性波模型［33］。在這個研究中，我們首先關注試驗區(qū)內(nèi)5°～15°范圍內(nèi)的淺源地震，通過模擬一些研究區(qū)的歷史地震測試PREM地球模型［34］。在區(qū)域內(nèi)利用該模型模擬長周期地震波（0.01～0.05Hz）的傳播是合理的。我們也利用試驗區(qū)內(nèi)和試驗區(qū)外的真實地震事件測試了搜索引擎，同時也測試了2個真實地震構造的假設雙震事件。

如圖1所示，試驗區(qū)為36°N～41°N，79°E～84°E，間隔為0.2°。深度網(wǎng)格從5～60km，間隔為5km。因此，在這個3D網(wǎng)格中大約有8 112個虛擬源。每個網(wǎng)格點的震源機制解離散如下：走向范圍從10°～350°，間隔為20°；傾角范圍從5°～80°，間隔為15°；滑動角范圍從－170°～170°，間隔為20°。這樣，在每個虛擬源共形成1 944個不同的事件。因此對于一個地震臺來說，我們都能產(chǎn)生15 789 168（8 122×1 944）三分量地震圖用于3D網(wǎng)格上的地震事件，也就是說，總共有47 367 504個三分量地震圖用于這3個地震臺。為構建搜索數(shù)據(jù)庫，我們將3個地震臺的所有三分量地震圖進行融合，并為每個虛擬地震構建一個長的超級道（supertrace）。因此，在每個搜索數(shù)據(jù)庫中共有15 789 168條超級道。

對于任意虛擬源位置，計算量最大的是在卷積震源機制解之前震源和接收器之間的格林函數(shù)的計算。對于一維地球模型來說，應該有9組格林函數(shù)，其中8組能用于計算雙力偶源的三分量地震圖［13，33］。幸運的是，格林函數(shù)不依賴于震源機制解，因此，我們僅需要計算每個網(wǎng)格點的一組格林函數(shù)。

此外，我們提出一個內(nèi)插方法，用于有效地計算離地震臺任一給定距離位置的理論地震圖。這個方法可以大大降低搜索數(shù)據(jù)庫構建的計算量，且精度可以接受。圖2對內(nèi)插方法進行了圖解說明。以地震臺為中心，以一定的插值間隔畫實線圓，同時，以實線圓為中心在兩邊畫虛線，標記插值區(qū)域。在地球變平換算后的任何一個插值間隔內(nèi)，我們對與實線圓一致的震中距上的唯一一個虛擬源計算格林函數(shù)。對于在3D網(wǎng)格上相同

圖2　一種用于準備搜索數(shù)據(jù)庫的插值方法。

針對一個地震臺，定義了具有一定間隔的實線圓。對于實線圓上的虛擬源，利用一維地球模型計算合成格林函數(shù)。對于落在兩條相鄰虛線之間相同間隔的網(wǎng)格點，利用與那些位于中央實線的點相同的格林函數(shù)作為近似值，如此可減少計算時間深度的任何虛擬源，如果它位于相同的間隔，那么它將被分配相同的格林函數(shù)作為近似值。在這個研究中，我們的震中距范圍是5°～15°。如果我們選擇0.2°作為間隔，這意味著我們僅必須為600個（50間隔×12深度）虛擬震源計算地震圖，而不必考慮臺站數(shù)量，與24 366（8 122×3臺站）個虛擬源相比，最高能降低40%計算工作量。在單個工作站上，在新疆試驗區(qū)進行地震圖計算構建數(shù)據(jù)庫需要花費大約13h，產(chǎn)出15 789 168個超級道并構建樹形結構用于快速搜索還需要再用30min。因此，如果有高性能的計算資源，創(chuàng)建一個智能搜索數(shù)據(jù)庫是一個非常有效的方法。

1.2 快速搜索結果

在新疆試驗區(qū)，2012年3月8日發(fā)生了MW5.9地震（圖1中紫色星所示）。根據(jù)全球矩張量解（CMT）記錄［35］，該地震事件發(fā)生在39.49°N，81.47°E，深度44.4km。我們利用此事件作為輸入，用我們提出的地震搜索引擎來確定其震中和震源機制解。此外，我們選擇2個其他地震用于特殊事件的系統(tǒng)測試。這些事件包括2014年4月30日發(fā)生的MW5.3地震，位于43.02°N，94.26°E，深度10.0km（圖1中綠色星所示；根據(jù)NEIC地震目錄，地震位于試驗區(qū)外［36］），以及發(fā)生于2011年9月15日的MW5.3地震，位于36.32°N，82.50°E，深度11.6km（圖1中藍色星所示；根據(jù)NEIC地震目錄，地震位于試驗區(qū)內(nèi)［36］）。對于發(fā)生在2011年9月15日的地震來說，我們將整個事件的發(fā)生時刻進行平移，來模擬其在2012年3月8日地震后40s發(fā)生的情況。通過兩個地震重疊構建人為的雙震記錄。我們的搜索數(shù)據(jù)庫并未專門考慮雙震或多震事件。然而，人造的雙震記錄可用于測試我們的搜索引擎方法如何處理此類事件。

圖3顯示了2012年3月8日地震輸入數(shù)據(jù)的超級道（上部紅線）以及數(shù)據(jù)庫中從1 到1 000，每間隔100的最佳匹配。由于每個超級道是由3個地震臺的三分量數(shù)據(jù)融合而成，而我們的搜索是基于從相似的事件中匹配9個可用的全波形。在圖3中，最佳搜索結果顯示的震源機制解，與全球矩張量解相

圖3　輸入數(shù)據(jù)和前1 000個搜索結果。

利用全球矩張量結果（紅色）將3個地震臺的數(shù)據(jù)與前1 000個搜索結果中每間隔100的結果（黑色）進行比較。左邊為震源機制解序列索引和最大互相關值，右邊為震源深度、經(jīng)緯度和震源機制解

利用全球矩張量結果（紅色）將3個地震臺的數(shù)據(jù)與前10個搜索結果（黑色）進行比較。左邊為震源機制解序列索引和最大互相關值，右邊為震源深度、經(jīng)緯度和震源機制解似，震源位置與全球矩張量解記錄相比，水平偏差15km，深度偏差0.6km。前200個結果表明，與全球矩張量解相比震源位置偏差在25km內(nèi)，震源深度偏差在5km內(nèi)。圖4顯示了在1 000個結果中的前10個最佳匹配結果。在所有10個搜索結果中的震源機制解都與全球矩張量解結果非常接近。

在10個最佳結果中，最大的互相關系數(shù)從0.869 5下降到0.862 6，變化很小，意味著這些結果存在非唯一性，這可能是受到有限的數(shù)據(jù)和觀測站幾何分布的影響。由于它定義了該結果的信度水平，因此了解這種情況非常重要。無論如何，這些結果表明震中位置、震源深度和震源機制解都能合理地包含在這個事件中。在這一特定的地震事件中，震源深度和震源機制解的估算可能比震中位置的估計更為可靠。

圖5顯示了在測試網(wǎng)格中震源深度為45km的最佳匹配結果，同時給出了全球矩張量解定位的輸入事件（紫色星）。在每個網(wǎng)格點，圖中表明了與輸入事件最佳匹配的震源機制解，且具有最大的互相關系數(shù)。我們的最佳搜索位置與全球矩張量解的經(jīng)緯度相差約半個網(wǎng)格（0.1°）。繪制像圖5中不同深度的震源機制解，能幫助我們理解在三維空間中解的不確定性。

圖6表明在最佳匹配位置下，矩張量解和最大互相關系數(shù)與深度的關系。虛線為我們的最佳估計震源深度是45.0km，而全球矩張量解的深度為44.4km。這兩個結果相當接近，意味著我們的結果與全球矩張量解一致。

圖7顯示了最佳的1 000個搜索結果的最大互相關系數(shù)曲線。該曲線給出了搜索結果非唯一性的直接說明。如果曲線隨著結果快速下降，意味著能很好地得到最佳結果。否則，表明存在太多的非唯一的或近似的結果。在這個研究中，褐色曲線快速下降，因此，該結果的信度水平很高。

我們也利用了2014年4月30日發(fā)生在試驗區(qū)外的地震測試搜索引擎。圖7中顯示，藍綠色曲線代表該事件中最佳匹配的1 000個搜索結果的最大互相關系數(shù)，所有的最大互相關系數(shù)都小于0.40，這是一個非

圖5　地震橫向定位不確定性。

對于震源深度為45 km，輸入值和搜索結果之間的最大互相關值的分布，以及最優(yōu)匹配結果的震源機制解。紫色星標記為全球矩張量解的定位結果。彩色比例尺表示最大互相關值的范圍

圖6　震源深度不確定性。

曲線表示在最優(yōu)匹配結果的地面點，震源深度在5～60 km，輸入數(shù)據(jù)和搜索結果之間的最大互相關值。實線表示全球矩張量結果震源深度為44.4 km的值，虛線表示快速搜索法的最優(yōu)匹配結果為震源深度45 km常低的值。這主要是因為在3個地震臺的搜索數(shù)據(jù)集中不能找到與發(fā)生在試驗區(qū)外地震相匹配的波形。同樣地，圖7中的橘色曲線表明，雙震的搜索也出現(xiàn)很低的互相關系數(shù)。通過利用人工和真實數(shù)據(jù)的實驗，我們發(fā)現(xiàn)互相關系數(shù)0.7是該系統(tǒng)在特定區(qū)域判別無效結果的閾值，包括地震發(fā)生在搜索數(shù)據(jù)集覆蓋范圍以外、多個地震疊加以及復雜破裂過程的大地震等情況。這意味著地震搜索引擎可以自動判別結果的有效性，避免錯誤報告。在這些測試事件中，當輸入值與搜索結果的最大互相關系數(shù)低于0.7時，這些搜索結果被認為是無效的。該研究結果意味著任何的已有搜索數(shù)據(jù)庫都應該進行更新，或者采用不同方法來更新數(shù)據(jù)。對于不同的區(qū)域和監(jiān)測臺網(wǎng)，在應用該系統(tǒng)前，都必須

圖7　三次地震測試搜索結果的比較。

實線表示輸入數(shù)據(jù)和前1 000個搜索結果之間的最大互相關值。褐色曲線表示試驗區(qū)內(nèi)2012年3 月8日發(fā)生的MW5.9地震的結果。藍綠色曲線表示試驗區(qū)外2014年4月30日發(fā)生的MW5.3地震的結果。橘色曲線表示由2012 年3月8日地震和2011年9月15日MW5.3地震構建出虛擬雙震的結果。虛線表示用于驗證搜索結果的互相關閾值進行重新估算和設置閾值。如果搜索數(shù)據(jù)庫可以包含特殊事件，搜索引擎方法應該能用于處理不同事件。我們應該探討許多可能性來改善數(shù)據(jù)庫。

2　討論

通過新疆測試區(qū)我們已經(jīng)表明，利用地震搜索引擎，在收到長周期面波后1s內(nèi)，我們能給出地震震中位置、震源深度和震源機制解。但不只是一個解，搜索引擎真正得到的是近似結果的子集，它也能幫助我們在1s內(nèi)來評估結果的信度水平。這項研究對于震源機制解的快速估算是一個重要的改進，它對于開展海嘯預警和監(jiān)測大型斷層活動更加有效。在這個研究過程中，我們并未探討震級參數(shù)。這是因為震級可以直接利用振幅信息進行估算［37］。

這個新方法需要在局部或者監(jiān)測區(qū)域范圍內(nèi)為試驗區(qū)構建一個搜索數(shù)據(jù)庫。我們的正演模擬計算僅限于點源，因此這個研究以震級小于7.5級地震為對象。然而，我們的搜索引擎并不會局限于特定震級大小或波譜類型。如果我們能構建一個包含模擬的大地震信息的數(shù)據(jù)庫，那么搜索引擎就能找出最佳波形匹配值，并能獲取地震震源機制參數(shù)。對于大地震的動態(tài)破裂過程來說，構建這樣一個搜索數(shù)據(jù)庫可以說是一個挑戰(zhàn)，也是很復雜的工作?？赡艿慕鉀Q方案包括，考慮震源規(guī)模對中低頻數(shù)據(jù)予以適當權重［2］，近似有限源格林函數(shù)合成計算［12］，以及利用最大方差縮減方法從搜索結果中定義震中位置［15］。一些地震確實在歷史上重復出現(xiàn)，在地震活躍區(qū)，利用歷史地震的震源信息構建一個數(shù)據(jù)庫，也將是一個值得關注的途徑。

我們的研究主要采用低頻面波作為輸入數(shù)據(jù)。對于局部或整個區(qū)域范圍內(nèi)，利用相對準確的三維速度模型，例如南加州，地震搜索引擎在該地區(qū)就可利用高頻體波。在研究中，由于體波能合理地利用格林函數(shù)進行估算，地震搜索引擎能在收到地震信息后的1s內(nèi)利用P波初動，快速得到震源機制解（包括震中和震級）。對于寬頻帶或高頻數(shù)據(jù)，可考慮利用人工合成輸入數(shù)據(jù)校準P波到時信息。幸運的是，在數(shù)據(jù)庫中校準無噪聲地震合成記錄相對容易。對于輸入數(shù)據(jù)，我們簡單地對于P波到時采用一系列較小的正負偏移，并同時利用所有原始記錄和偏移記錄進行快速搜索。我們選取與最小誤差值相關的結果作為最后的結果。這個方法有助于消除P波到時不準確的影響。

在這個研究中，我們在新疆選擇了一個5°×5°的范圍，震源深度在5～60km，在單一的CPU計算機上測試我們的地震搜索引擎。如果我們在空間上擴大研究區(qū)范圍，或對虛擬源和局部結果減小網(wǎng)格插值間隔，搜索數(shù)據(jù)庫也將增加更多信息。我們用于計算格林函數(shù)的插值方法表明，對于合成數(shù)據(jù)集來說，計算時間僅僅與震中與臺站間的距離相關；因此，密集臺站或者監(jiān)測網(wǎng)并不會顯著增加計算時間。考慮到不同頻率體波和面波的相關數(shù)據(jù)信息的重要性，我們根據(jù)同一地震不同數(shù)據(jù)信息，設計了多線程搜索方式，以確保結果的可靠性。我們也通過采用搜索引擎法實時對波形數(shù)據(jù)進行分組，并在數(shù)據(jù)接收過程中對結果進行不斷更新。我們將進行進一步的研究工作。

3　方法

3.1 計算機快速搜索技術

在本研究中，地震搜索引擎利用類似于網(wǎng)絡搜索引擎的快速搜索方法，后者主要用于在網(wǎng)絡上查找以文字為主的信息，大量的搜索結果以索引和歸類的方式存放于一個大型數(shù)據(jù)庫中［20］。圖像搜索引擎與其類似，但它是基于圖像的，并對內(nèi)容以及其他信息進行相似性評估，用于獲取與記錄圖像的最優(yōu)匹配值［21－23］。網(wǎng)絡技術已經(jīng)激發(fā)其他領域的相關研究，例如，Aguiar和Beroza［38］應用網(wǎng)絡排序方法分析地震波的相似性。

在計算機搜索領域，數(shù)據(jù)中的大量采樣點被作為它的“維度”。地震圖的相似性搜索是一個高維問題，需要努力地從大型數(shù)據(jù)集中搜索真正的結果。降維法能在保持數(shù)據(jù)本質(zhì)特征的同時，有效地減少時間樣本的數(shù)量［39］。除了降維以外，一個重要的方法是近似臨近搜索法，它能比精確查找快很多，而且能提供近似最優(yōu)精確度［24］。

已有的快速搜索方法多是基于散列結構或樹結構，例如區(qū)域敏感散列結構［25］和多重隨機K維樹方法（MRKD）［26－27］。通過利用這兩種方法處理地震圖，我們發(fā)現(xiàn)MRKD樹結構方法用于處理高維地震圖數(shù)據(jù)集的變化始終很快。Silpa－Anan和Hartley［26］及Muja和Lowe［24］認為，MRKD樹方法能很好地處理計算機視覺中的高維數(shù)據(jù)集問題。這個方法通過在樹結構的每一層上將振幅變化最大的波形數(shù)據(jù)拆分成兩半，從數(shù)據(jù)集構建多重樹結構。數(shù)據(jù)庫利用有限樣本獲取地震數(shù)據(jù)中最重要的特征，然后當有記錄到達時，我們根據(jù)樹結構搜索最優(yōu)匹配結果。這個方法需要通過log2N次計算比較，從一個樹結構中找出第一個候選作為最優(yōu)匹配結果，N代表數(shù)據(jù)庫中地震圖的數(shù)量［28］。如果需要更多的搜索結果，則需要附加的回溯工作。

3.2 利用MRKD樹對地震圖進行索引和分類

圖8舉例說明了如何利用數(shù)據(jù)庫中虛擬的4個地震圖構建K維樹。4個地震圖的初至波被記錄下來并對齊排列。在每個時間樣本中，計算4個地震圖的振幅平均值。圖8a的最下面一條曲線顯示了每個時間樣本振幅與平均值的方差。最大的振幅方差為0.143，平均值為0.176，從樹結構中挑選出的是第864個采樣點（維）。選擇相同的時間樣本，振幅比平均值小的地震圖（ID＝1，4）放置于維節(jié)點的左側，而那些振幅大于平均值的地震圖（ID＝2，3）位于右側。上述進程被遞歸式地用于左右單獨排列的地震圖，到任意維節(jié)點下僅留下一個地震圖（葉節(jié)點）為止。

圖8b顯示了在以上進程后形成的二級樹。實際上，應該利用同一數(shù)據(jù)集構建多重隨機K－D樹并用于搜索。為構建多重樹，我們可以利用Muja和Lowe的方法［24］，從前m個具有最大方差的維中隨機地選擇一個分離維。在我們的應用程序中，利用固定值m＝5，從一個大型地震數(shù)據(jù)集中構建了128個樹。當需要查詢的地震圖輸入時，相同維（時間采樣點）的數(shù)據(jù)集相對于平均值的振幅方差將與分區(qū)值進行比較，并用于確定哪一半數(shù)據(jù)屬于查詢數(shù)據(jù)。這個搜索過程是在所有樹中分別進行的。在每個維節(jié)點，計算輸入振幅值和數(shù)據(jù)集平均值間的累積L2距離值，同時通過增加L2距離值，保留所有隨機樹中單一的優(yōu)先隊列。在地震圖ID有效的每個葉節(jié)點，計算輸入和合成地震圖之間

圖8　快速搜索法標示地震波形。

MRKD樹法用于構建快速搜索法數(shù)據(jù)庫。（a）假設4個地震圖的數(shù)據(jù)集，對于所有的道，計算每個數(shù)據(jù)維的平均值和RMS方差；（b）在樹頂上，最大方差（0.143）出現(xiàn)在第864個樣本，左邊標示該樣本的振幅值小于平均值（0.176），右邊為那些大于平均值的樣本。這個進程迭代地應用于每一層，直到每一邊僅留下一個地震圖為止的L2標準距離，并根據(jù)距離遞增順序保存在候選隊列中。該算法需要迭代搜索來查找許多近似解。優(yōu)先隊列確定了用于下一次迭代對比的地震圖，并且候選隊列返回大量按照與輸入地震圖相似性遞減順序排列的結果［26］。該搜索方法是基于利用最高概率的理念，選擇輸入值的最相鄰近似值。其原理在地震圖匹配中是可行的，因為本質(zhì)上對于真實地震圖的精確匹配是不存在的，這是由于我們目前對于地球三維結構認知不完整，以及準確反演真實波傳播效果的能力有限。

3.3 優(yōu)化快速搜索性能

降維是加速搜索進程的另一個重要方法。其目標是識別數(shù)據(jù)模式及壓縮模式下的數(shù)據(jù)表達，以突出數(shù)據(jù)的相似性和差異性。在各種降維方法中，主成分分析法（PCA）被認為是能可靠地將數(shù)據(jù)嵌入到低維線性子空間的好方法［39］。為了應用主成分分析法，我們首先計算數(shù)據(jù)集的協(xié)方差矩陣，然后求解其特征向量和特征值。對于n維的數(shù)據(jù)集來說，它應該有n個特征向量和特征值。如果我們僅選擇p個具有最大特征值的特征向量作為投影基，那么最后數(shù)據(jù)集僅有p維，并不會損失太多信息。實際上，我們需要檢查特征值的分布并確定p的最優(yōu)值。

選擇最優(yōu)參數(shù)用于快速搜索方法非常重要。圖9說明了我們?nèi)绾芜x擇樹的數(shù)量，搜索結果的數(shù)量和主成分分析法的輸入維，用于處理一個數(shù)據(jù)集中超過1 500萬的合成地震圖。圖9a顯示了在MRKD樹中，返回的1 000個結果中搜索精度與樹的數(shù)量的關系。搜索精度被定義為MRKD樹與精確查找1 000個搜索結果的比值。對于128樹來說，其精確度大約是76%，意味著在給出的1 000個結果中，760種快速搜索法的結果與該搜索法的精確搜索一致。該精確度應該足以用于選擇最優(yōu)波形匹配。圖9b表示搜索前1 000個解的精度與搜索候選解個數(shù)之間的關系，其中用了5個不同的合成地震波形作為輸入。在前1 000個結果中，更多的搜索返回值產(chǎn)出更多的精確結果，但其需要更多的計算量?？雌饋磉x擇1 000個結果和128樹能在精度和計算量之間達到一個好的平衡。對于以上的測試數(shù)據(jù)來說，原始記錄包含3 072個時間樣本（維）。圖9c表明特征值與從數(shù)據(jù)協(xié)方差矩陣獲取的維數(shù)的關系。該曲線可用于確定最優(yōu)p值，并用于上述探討的主成分分析法降維。該曲線表明，在100維以上的特征值非常接近于0。因此，可

圖9　快速搜索法的最優(yōu)參數(shù)估計。

搜索精度被定義為MRKD樹與精確查找超過1 000個搜索結果的比值。虛線表示參數(shù)的選擇。（a）為1 000個搜索結果的搜索精度與樹的數(shù)量的關系。選擇128樹的精度約為76%。（b）為利用128樹，前1 000個結果的搜索精度與基于5個不同合成地震記錄測試搜索候選方案數(shù)量的關系。總共1 000個輸出結果，足以達到約76%的精度。（c）為利用主成分分析降維后，數(shù)據(jù)協(xié)方差矩陣特征值與維數(shù)的關系。該曲線表明，大于100維的特征值非常接近于0，因此可以去除。新數(shù)據(jù)集的維數(shù)可降至100，且不會造成太多信息損失以去除這些與小的特征值相關的特征向量，而不會造成許多信息損失。選取100作為最優(yōu)p值，新數(shù)據(jù)集的維數(shù)可降至100。

基于數(shù)據(jù)集概率的計算機搜索算法的性能差異很大［24］。對于三分量地震臺站來說，我們首先使波形歸一化，然后融合三分量數(shù)據(jù)生成每個事件的超級道。圖10表明MRKD樹法在未利用主成分分析法時，大約

圖10　不同方法下的計算機搜索時間對比。

對于返準L2需要超過17min來完成相同的工作。利用互相關查找最優(yōu)匹配需要27h。這些測試都通過單一的AMD6136處理器進行。如上所述，MRKD樹法在單一樹需要log2N次對比來查找第一個候選用于最優(yōu)匹配，N表示數(shù)據(jù)庫中地震圖的數(shù)量。然而，在圖10中，MRKD樹搜索時間與地震數(shù)量的關系在對數(shù)坐標系中幾乎是線性的。這是因為我們實際上將數(shù)據(jù)集分割到多個更小的數(shù)據(jù)集，用于獨立順序檢索，然后再合成結果。將來，數(shù)據(jù)集將構建成并行計算的形式?；氐? 000個最佳匹配結果，在單一CPU上利用3種不同的方法進行測試。這些方法包括確定性線性搜索（圓圈標示），未進行降維的MRKD樹法（星號標示）和利用主成分分析法降維的MRKD樹法（三角形標示）需要0.2s在3 300萬超級道輸入值中查找1 000個最優(yōu)匹配，而利用主成分分析法，僅需要大約0.06s。精確查找計算不匹配的標

資料來源：Zhang J，Zhang H，Chen E，et al.Real－time earthquake monitoring using a search engine method.Nature Communications，2014，5：5664.doi：10.1038／ncomms6664

（福建省地震局 王林 譯，黃宏生 校；中國地震局地球物理研究所 王輝 復校）

（譯者電子信箱，王林：wl＿0117＠163.com）

參考文獻

［1］Allen R M.Seconds before the big one.Sci.Am.，2011，304：75－79

［2］Ekstrom G，Nettles M，Dziewonski A M.The global CMT project 2004—2010：centroid－moment tensors for 13017earthquakes.Phys.Earth Planet.Inter.，2012，200－201：1－9

［3］Satriano C，Lomax A，Zollo A.Real－time evolutionary earthquake location for seismic early warning.Bull.Seis.Soc.Amer.，1998，98：1482－1494

［4］Allen R M，Gasparini P，Kamigaichi O，et al.The status of earthquake early warning around the world：an introductory overview.Seismol.Res.Lett.，2009，80：682－693

［5］Nakamura H，Horiuchi S，Wu C，et al.Evaluation of the real－time earthquake information system in Japan.Geophys.Res.Lett.，2008，36：L00B01

［6］Espinosa－Aranda J M，Rodríquez F H.Mexico City seismic alert system.Seismol.Res.Lett.，1995，66：42－53

［7］Wu Y M，Teng T L.A virtual sub－network approach to earthquake early warning.Bull.Seis.Soc.Amer.，2002，92：2008－2018

［8］Erdik M，Aydinoglu N，F(xiàn)ahjan Y，et al.Istanbul earthquake rapid response and the early warning system.Bull.Earthquake Eng.，2003，1：157－163

［9］Horiuchi S，Negishi H，Abe K，et al.An automatic processing system for broadcasting earthquake alarms.Bull.Seis.Soc.Amer.，2005，95：708－718

［10］Kanamori H.Real－time seismology and earthquake damage mitigation.Annu.Rev.Earth Planet.Sci.，2005，33：195－214

［11］Lee J，F(xiàn)riederich W，Meier T.Real time monitoring of moment magnitude by waveform inversion.Geophys.Res.Lett.，2012，39：L02309

［12］Guilhem A，Dreger D S.Rapid detection and characterization of large earthquakes using quasi－finitesource Green’s functions in continuous moment tensor inversion.Geophys.Res.Lett.，2011，38：L13318

［13］Tsuruoka H，Kawakatsu H，Urabe T.GRiD MT（grid－based real－time determination of moment tensors）monitoring the log－period seismic wavefield.Phys.Earth Planet.Inter.，2009，175：8－16

［14］Melgar D，Bock Y，Crowell B W.Real－time centroid moment tensor determination for large earthquakes from local and regional displacement records.Geophys.J.Int.，2012，188：703－718

［15］Crowell B W，Bock Y，Melgar D.Real－time inversion of GPS data for finite fault modeling and rapid hazard assessment.Geophys.Res.Lett.，2012，39：L09305

［16］Papazachos B C，Dimitriu P P.Tsunamis in and near Greece and their relation to the earthquake focal mechanisms.Natural Hazards，1991，4：161－170

［17］Ismoyo B.Indonesia battles disasters on two fronts.The Jakarta Globe 1.2010.http：∥thejakartaglobe.beritasatu.com／archive／indonesia－battles－disasters－ontwo－fronts／404050／

［18］Yue H，Lay T，Koper K D.En e′chelon and orthogonal fault ruptures of the 11April 2012great intraplate earthquakes.Nature，2012，490：245－249

［19］Bouchon M，Karabulut H，Aktar M，et al.Extended nucleation of the 1999 Mw 7.6Izmit earthquake.Science，2011，331：877－880

［20］Henzinger M.Search technologies for the Internet.Science，2007，317：468－471

［21］Bach J R，F(xiàn)uller C，Gupta A，et al.Virage image search engine：an open framework for image management.SPIE，1996，2670：76－87

［22］Sahami M，Mittal V，Baluja S，et al.The happy searcher：challenges in web information retrieval.Proc.8th Pacific Rim Int.Conf.Artificial Intelligence，Springer，2004，3157：3－12

［23］Zhou W，Li H，Lu Y，et al.Large scale image retrieval with geometric coding.ACM International Conference on Multimedia（MM），2011：1349－1352

［24］Muja M，Lowe D G.Fast approximate nearest neighbors with automatic algorithm configuration.Proc.International Conf.Computer Vision Theory and Applications（VISAPP’09），2009，1：331

［25］Slaney M，Casey M.Locality－sensitive hashing for finding nearest neighbors.IEEE Signal Process.Mag.，2008，25：128－131

［26］Silpa－Anan C，Hartley R.Optimised KD－trees for fast image descriptor matching.Computer Vision and Pattern Recognition 1－8，Institute of Electrical and Electronics Engineers（IEEE），2008

［27］Wu P，Hoi S C H，Nguyen D D，et al.Randomly projected KD－trees with distance metric learning for image retrieval.Proc.17th Int.Conf.Adv.Multi.Model.Springer，2011：371－382

［28］Friedman J H，Bentley J L，F(xiàn)inkel R A.An algorithm for finding best matches in logarithmic expected time.ACM Trans.Math.Softw.，1977，3：209－226

［29］Berchtold S，Bhm C，Keim D A，et al.A cost model for nearest neighbor search in high－dimensional data space.Proc.16th ACM SIGACTSIGMOD－SIGART Symp.Principles of Database Systems，ACM Press，1997：78－86

［30］Yianilos P N.Data structure and algorithms for nearest neighbor search in general metric spaces.Proc.5th ACM－SIAM Symp.Discr.Algorithms，Society for Industrial and Applied Mathematics，1993：311－321

［31］Seidl T，Kriegel H P.Optimal multi－step k－nearest neighbor search.Proc.ACM SIGMOD Int.Conf.on Management of Data 27，SIGMOD，1998：154－165

［32］Arya S，Mount D M，Netanyahu N S，et al.Proc.of the Fifth Annual ACM－STAM Symp.on Discrete Algorithms.Association for Computing Machinery（ACM），1994：573－582

［33］Zhu L，Rivera L A.A note on the dynamic and static displacements from a point source in multilayered media.Geophys.J.Int.，2002，148：619－627

［34］Dziewonski A M，Anderson D L.Preliminary reference Earth model.Phys.Earth Planet Inter.，1981，25：297－356

［35］Ekstrom G，Nettles M.Global CMT Web Page.2014.http：∥www.globalcmt.org／

［36］National Earthquake Information Center—NEIC.2014.http：∥earthquake.usgs.gov／regional／neic／

［37］Shearer P M.Introduction to Seismology.Cambridge Univ.Press，1999

［38］Aguiar A C，Beroza G C.Page rank for earthquakes.Seismol.Res.Lett.，2014，85：344－350

［39］Jolliffe I T.Principal Component Analysis.Springer－Verlag，1986

學術論文

收稿日期：(＊)2015－01－26；采用日期：2015－07－10。

doi：10.3969／j.issn.0235－4975.2015.10.003

文獻標識碼：A；

中圖分類號：P315.63；