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

?

搜索引擎方法在實(shí)時(shí)地震監(jiān)測(cè)中的應(yīng)用＊

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

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

關(guān)鍵詞搜索引擎方法；實(shí)時(shí)地震監(jiān)測(cè)；震源機(jī)制解

引言

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

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

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

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

1　計(jì)算結(jié)果

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

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

圖1　地震臺(tái)站和虛擬試驗(yàn)區(qū)的位置。

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

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

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

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

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

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

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

1.2 快速搜索結(jié)果

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

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

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

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

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

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

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

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

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

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

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

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

圖6　震源深度不確定性。

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

圖7　三次地震測(cè)試搜索結(jié)果的比較。

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

2　討論

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

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

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

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

3　方法

3.1 計(jì)算機(jī)快速搜索技術(shù)

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

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

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

3.2 利用MRKD樹對(duì)地震圖進(jìn)行索引和分類

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

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

圖8　快速搜索法標(biāo)示地震波形。

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

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

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

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

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

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

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

圖10　不同方法下的計(jì)算機(jī)搜索時(shí)間對(duì)比。

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

資料來源：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

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

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

參考文獻(xiàn)

［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

學(xué)術(shù)論文

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

doi：10.3969／j.issn.0235－4975.2015.10.003

文獻(xiàn)標(biāo)識(shí)碼：A；

中圖分類號(hào)：P315.63；