王 娟 劉俊民 李 健 邱宏茂 王曉明 商 杰 蓋 磊

（禁核試北京國家數(shù)據(jù)中心和北京放射性核素實(shí)驗(yàn)室，北京 100085）

引言

地震監(jiān)測(cè)是全面禁止核試驗(yàn)條約規(guī)定的4種核查手段中遠(yuǎn)區(qū)監(jiān)測(cè)最有效的方法。條約生效后，CTBTO國際監(jiān)測(cè)系統(tǒng)（IMS）通過與全球通訊基礎(chǔ)設(shè)施（GCI）的安全衛(wèi)星鏈接傳送位于全球170個(gè)地震臺(tái)站的數(shù)據(jù)至國際數(shù)據(jù)中心（IDC），其中50個(gè)基本地震臺(tái)站1周7天1天24小時(shí)不間斷地傳送數(shù)據(jù)，120個(gè)輔助地震臺(tái)站在需要時(shí)，由程序自動(dòng)按照既定的規(guī)則申請(qǐng)傳送相應(yīng)時(shí)間段的數(shù)據(jù)[1]。地震監(jiān)測(cè)數(shù)據(jù)進(jìn)入IDC的數(shù)據(jù)處理系統(tǒng)后，首先進(jìn)行臺(tái)站處理：即檢測(cè)信號(hào)、識(shí)別震相、計(jì)算震相特性（到時(shí)、振幅、方位角、慢度等）。單臺(tái)臺(tái)站處理完成后，根據(jù)計(jì)算出的震相特征，利用臺(tái)網(wǎng)關(guān)聯(lián)算法將來自不同臺(tái)站的信號(hào)進(jìn)行關(guān)聯(lián)，推測(cè)這些信號(hào)可能來自同一事件，最終形成自動(dòng)處理公報(bào)。分析員以自動(dòng)處理公報(bào)為基礎(chǔ)，經(jīng)過詳細(xì)分析，添加自動(dòng)關(guān)聯(lián)漏掉的事件，刪除虛假關(guān)聯(lián)的事件，最終形成審定事件公報(bào)（REB）。目前，IDC的REB公報(bào)也已經(jīng)成為全球地震領(lǐng)域比較認(rèn)可的公報(bào)之一[2]。想要提高自動(dòng)處理事件關(guān)聯(lián)的有效率，降低事件漏檢率，提高自動(dòng)處理公報(bào)質(zhì)量，減小分析員的工作量，是IDC一直以來致力研究的方向之一。

目前在IDC數(shù)據(jù)處理系統(tǒng)中，臺(tái)網(wǎng)關(guān)聯(lián)算法采用全局關(guān)聯(lián)算法（GA），其采用了基于全球格點(diǎn)的搜索算法來關(guān)聯(lián)震相形成事件[3]。GA首先將全球按緯度、經(jīng)度分為覆蓋整個(gè)地球表面相互重疊的若干個(gè)圓形格點(diǎn)，格點(diǎn)文件包括了格點(diǎn)位置信息以及與臺(tái)站間的關(guān)系，然后對(duì)地球上給定位置（格點(diǎn)中心點(diǎn)）的所有可能檢測(cè)組合進(jìn)行全局搜索，組合形成許多假設(shè)的事件，隨后合并相同事件，剔除冗余和虛假事件，最終形成公報(bào)產(chǎn)品。GA在IDC運(yùn)行20多年來，在波形數(shù)據(jù)的實(shí)時(shí)處理中表現(xiàn)出了良好的處理性能，但是GA還是漏掉了感興趣的震級(jí)范圍內(nèi)將近三分之一的地震事件，且自動(dòng)處理公報(bào)中大約一半的事件為虛假事件[4-5]。

隨著計(jì)算機(jī)性能的提高、大數(shù)據(jù)分析、數(shù)據(jù)挖掘及機(jī)器學(xué)習(xí)等先進(jìn)數(shù)據(jù)處理技術(shù)的廣泛應(yīng)用，如何利用海量的歷史地震事件的特征來提高事件關(guān)聯(lián)性能成為學(xué)者們研究的熱點(diǎn)和IDC關(guān)注的重點(diǎn)之一。2009年開始，美國加州大學(xué)伯克利分校（University of California，Berkeley）的科學(xué)家與IDC開始合作研究基于震相特征模型的網(wǎng)絡(luò)處理-垂直綜合地震分析（Network processing-Vertically Integrated Seismic Analysis，NET-VISA）軟件[4-5]，旨在替代GA成為IDC下一代全球震相自動(dòng)關(guān)聯(lián)和事件檢測(cè)定位軟件。本文介紹了基于貝葉斯推理方法的NET-VISA算法的基本原理及在IDC歷年測(cè)試中的進(jìn)展，并對(duì)禁核試北京國家數(shù)據(jù)中心（NDC）參加IDC組織的阿爾法測(cè)試（Alpha Tester Group，ATG）中朝鮮 6 次核試驗(yàn)的 NET-VISA關(guān)聯(lián)結(jié)果作了分析對(duì)比。

1 NET-VISA 基本原理

NET-VISA核心是基于貝葉斯框架的關(guān)聯(lián)算法，由基于物理的生成概率模型（GM）和啟發(fā)式推理算法（IA）組成。該方法首先建立一個(gè)綜合考慮了信號(hào)到時(shí)、方位角、慢度、事件震級(jí)、正確震相、誤檢信號(hào)和尾波震相等事件特征的全球尺度的概率模型，然后由基于貝葉斯方法設(shè)計(jì)的推理算法（IA）搜索具有最大后驗(yàn)概率（MAP）的事件集，并在關(guān)聯(lián)到達(dá)的事件之間進(jìn)行迭代，最終確定事件的時(shí)間、位置和大小[5-6]。

假設(shè)X為覆蓋所有可能地震事件集合的一個(gè)隨機(jī)變量，每個(gè)事件由時(shí)間、位置、深度和震級(jí)來定義，Y為臺(tái)網(wǎng)內(nèi)全部地震臺(tái)站的所有可能的信號(hào)觀測(cè)值，Pθ(X)描述了參數(shù)化的事件先驗(yàn)概率，P?(Y|X)描述了信號(hào)的傳播和測(cè)量概率（包括走時(shí)、吸收與散射、噪聲等）。對(duì)于給定的觀測(cè)值Y=y，事件的后驗(yàn)分布P(X|Y=y)即為事件參數(shù)。則NET-VISA的事件模型可以由最基本的貝葉斯公式來描述：

在NET-VISA中，一個(gè)到達(dá)信號(hào)最終分為關(guān)聯(lián)事件的震相 Λ ，誤檢信號(hào) ξ和尾波信號(hào) κ，在給定觀測(cè)信號(hào)A的情況下，NET-VISA結(jié)合了事件假設(shè)e，震相Λ，誤檢信號(hào)ξ 和尾波信號(hào) κ得到的聯(lián)合概率模型[3]：

推理算法的最終目的是得到使后驗(yàn)概率達(dá)到最大的事件參數(shù)的解：

NET-VISA生成概率模型由先驗(yàn)?zāi)Ｐ秃退迫缓瘮?shù)模型構(gòu)成，主要包括事件先驗(yàn)?zāi)Ｐ?、震相似然函?shù)模型、誤檢信號(hào)模型和尾波信號(hào)模型等部分。NETVISA中只考慮了14種震相：即P、Pn、PKP、S、PKPbc、PcP、pP、Lg、PKPab、PKKPbc、Pg、Rg、Sn和 ScP的震相模型。事件先驗(yàn)?zāi)Ｐ头从沉耸录l度、地理分布、深度分布和震級(jí)分布等；震相似然函數(shù)模型根據(jù)統(tǒng)計(jì)的事件發(fā)生時(shí)的震相檢測(cè)概率和震相屬性特征而建立；誤檢信號(hào)模型描述了由噪聲產(chǎn)生的誤檢信號(hào)的特征在各臺(tái)站的分布；尾波信號(hào)模型描述了由于信號(hào)傳播的散射而帶來的尾波信號(hào)，減少其對(duì)事件關(guān)聯(lián)的影響，通過尾波信號(hào)與其觸發(fā)震相的延遲、尾波信號(hào)的方位角、慢度和振幅特征來構(gòu)建模型[3-5]。

NET-VISA的輸入從信號(hào)檢測(cè)開始，并基于所有標(biāo)為P震相的檢測(cè)建立種子事件。通過構(gòu)建事件、改善到達(dá)、改善事件、刪除事件等步驟查找最佳事件集以解釋觀察到的檢測(cè)結(jié)果，NET-VISA在改善到達(dá)和改善事件之間進(jìn)行迭代以優(yōu)化事件。圖1為NETVISA生成事件的示意圖。

圖1 NET-VISA 事件生成窗口Fig.1 NET-VISA windows

由圖1可知，NET-VISA將連續(xù)的到達(dá)數(shù)據(jù)流劃分為不同的窗口。在構(gòu)建事件窗口中利用每個(gè)獨(dú)立震相的相關(guān)特征參數(shù)構(gòu)建種子事件，圖1中構(gòu)建事件窗口的窗長(zhǎng)為W，窗長(zhǎng)從t0=0開始到t1=t0+W（實(shí)際應(yīng)用中W取 30 min）。到達(dá)窗口從t0?MT到t2=t0+W+MT，其中MT是震相的最大走時(shí)（33.3 min 或2 000 s），在此窗口中的每個(gè)到達(dá)都要進(jìn)行與事件窗口中的事件相關(guān)聯(lián)。對(duì)于每個(gè)檢測(cè)在檢測(cè)窗口內(nèi)通過“檢測(cè)改善算法”后，與事件窗口中的事件相關(guān)聯(lián)，進(jìn)行迭代以優(yōu)化事件。NET-VISA采用定義的事件度量標(biāo)準(zhǔn)和檢測(cè)評(píng)分標(biāo)準(zhǔn)來計(jì)算檢測(cè)和事件的評(píng)分，事件得分小于1的事件將被刪除。構(gòu)建事件窗口向前移動(dòng)步長(zhǎng)S后，在輸出窗口中的最終確定的事件將寫入到事件公報(bào)中。

2 NET-VISA 測(cè)試進(jìn)展[7-11]

2009年3月，美國加州大學(xué)伯克利分校的斯圖亞特·羅素（Stuart Russell）教授在其主持的 ISS09（International Scientific Studies Project 2009，國際科學(xué)研究項(xiàng)目2009）項(xiàng)目中首次提出了NET-VISA算法。后來Arora博士在2009年6月的ISS09大會(huì)上報(bào)告了有限數(shù)據(jù)集的測(cè)試結(jié)果。2009年11月，在ISS09技術(shù)會(huì)議上討論了基于3個(gè)月的地震數(shù)據(jù)研究結(jié)果，當(dāng)時(shí)相對(duì)于SEL3雖有改善，但還沒有量的體現(xiàn)。2010年2月，Russell教授在籌委會(huì)第34次核查工作組會(huì)議上報(bào)告了NET-VISA原型系統(tǒng)在以3個(gè)月數(shù)據(jù)集為基礎(chǔ)的測(cè)試結(jié)果，相對(duì)于SEL3，系統(tǒng)性能提高了5.5%。表1為當(dāng)時(shí)的統(tǒng)計(jì)數(shù)據(jù)，表中SEL3代表用全局關(guān)聯(lián)算法（GA）產(chǎn)生的自動(dòng)處理公報(bào)事件，LEB是分析員審定后的事件。

表1 NET-VISA 原型系統(tǒng) 3 個(gè)月測(cè)試結(jié)果Table 1 The three months test results of NET-VISA prototype

隨后Arora博士于2012—2013年針對(duì)地震數(shù)據(jù)處理研究開發(fā)了軟件，2014年集成了水聲處理流程，2015年開始了對(duì)次聲數(shù)據(jù)的處理。在PTS的組織下，Arora博士與IDC的科學(xué)家在2015年以離線方式進(jìn)行了1年的地震數(shù)據(jù)處理測(cè)試，2016—2017年以離線運(yùn)行方式針對(duì)包含次聲數(shù)據(jù)處理的NETVISA進(jìn)行了測(cè)試，2018年1月—2019年2月在IDC測(cè)試平臺(tái)上對(duì)NET-VISA進(jìn)行了在線測(cè)試。測(cè)試結(jié)果表明，在一般情況下，NET-VISA產(chǎn)生的公報(bào)要比IDC自動(dòng)處理產(chǎn)生的公報(bào)SEL3更加完整和精準(zhǔn)，對(duì)于REB的覆蓋率可以提高約15%（80%—90%與65%—75%），不一致率降低約5%（45%—55%與50%—60%）。

2019年P(guān)TS針對(duì)NET-VISA軟件的測(cè)試成立了ATG，NDC參加了阿爾法測(cè)試全過程，順利完成了全部4項(xiàng)測(cè)試，其中第2項(xiàng)測(cè)試是對(duì)2017年9月3日朝鮮核試驗(yàn)事件進(jìn)行數(shù)據(jù)處理，我們進(jìn)行了擴(kuò)展測(cè)試，對(duì)朝鮮6次核試驗(yàn)逐個(gè)進(jìn)行了類似處理，對(duì)NET-VISA關(guān)聯(lián)結(jié)果進(jìn)行了分析對(duì)比。

3 朝鮮核試驗(yàn) NET-VISA 測(cè)試結(jié)果分析

利用阿爾法測(cè)試中IDC提供的NET-VISA程序?qū)Τr6次核試驗(yàn)數(shù)據(jù)逐個(gè)進(jìn)行了重關(guān)聯(lián)。為了方便在同等條件下與IDC/GA關(guān)聯(lián)算法比較，在臺(tái)站列表方面，選擇IDC/SEL3朝鮮核試驗(yàn)事件關(guān)聯(lián)臺(tái)站表中全部地震臺(tái)站；在數(shù)據(jù)段方面，以REB事件前約10分鐘整分為起始時(shí)間，持續(xù)時(shí)長(zhǎng)2小時(shí)，以保證全部可能震相可被關(guān)聯(lián)。經(jīng)NET-VISA關(guān)聯(lián)后的朝鮮地區(qū)地震事件基本參數(shù)列于表2。其中第4、5、6次核試驗(yàn)事件均被NET-VISA分裂成兩個(gè)事件；第6次核試驗(yàn)事件過后約8分半鐘還有一次小震級(jí)事件。

表2 NET-VISA定位的歷次朝鮮核試驗(yàn)相關(guān)事件基本參數(shù)列表Table 2 Basic parameters of events related to several North Korean nuclear tests located by NET-VISA

以IDC/REB為參考，考察SEL3和VSEL（NETVISA產(chǎn)生的公報(bào)）事件定時(shí)和定位精準(zhǔn)程度和存在的問題。定位和定時(shí)誤差按以下公式計(jì)算：

其中，E代表誤差，MSE代表均方誤差。

SEL3和VSEL定時(shí)誤差分別列于表3，圖2是對(duì)應(yīng)表3的定時(shí)誤差直方圖。定位誤差列于表4，圖3是對(duì)應(yīng)表4的定位誤差直方圖。表3和表4中VSEL2是VSEL的分裂事件，直方圖中NetVISA1對(duì)應(yīng)于VSEL，NetVISA2對(duì)應(yīng)于VSEL2，用2017-09-03（1）和 2017-09-03（2）區(qū)別 2017 年 9月第 1 次事件與第2次事件。定時(shí)均方誤差MSE（SEL3，time）=0.796 1 s、MSE（VSEL，time）=2.006 7 s，定位均方誤差MSE（SEL3，dis.）=0.10051°、MSE（VSEL，dis.）=0.44713°。

圖3 SEL3 和 VSEL 定位誤差直方圖Fig.3 Histogram of SEL3 and VSEL location errors

表4 以 REB 為參考的 SEL3 和 VSEL 定位誤差（°）Table 4 SEL3 and VSEL location errors （°） with REB as reference

圖2 SEL3 和 VSEL 事件定時(shí)誤差直方圖Fig.2 Histogram of timing error of SEL3 and VSEL events

表3 以 REB 為參考的 SEL3 和 VSEL 定時(shí)誤差（s）Table 3 SEL3 and VSEL timing errors （s） with REB as reference

從測(cè)試結(jié)果來看，NET-VISA軟件模塊運(yùn)行基本平穩(wěn)、可靠。NetVISA檢出了全部核試驗(yàn)事件，但是把第4、5、6次核試驗(yàn)事件均分裂成了兩次事件；事件定時(shí)和定位方面，NetVISA比IDC的GA稍差一點(diǎn)。

4 結(jié)論

NET-VISA作為IDC下一代臺(tái)網(wǎng)關(guān)聯(lián)軟件，從2009年概念提出、立項(xiàng)、研發(fā)、測(cè)試，前后已經(jīng)過了10年，目前正在對(duì)NET-VISA軟件進(jìn)行發(fā)布前的用戶體驗(yàn)測(cè)試，雖然IDC的在線和離線測(cè)試均表明NET-VISA產(chǎn)生的公報(bào)要比IDC目前使用的全球關(guān)聯(lián)方法生成的自動(dòng)處理公報(bào)SEL3更加完整和精準(zhǔn)，但是ATG測(cè)試組的報(bào)告仍然反應(yīng)出不少問題，這些問題反饋到PTS，有助于對(duì)NET-VISA軟件進(jìn)行改進(jìn)及模型優(yōu)化。下一步，IDC將結(jié)合波形整體特征及更豐富的歷史信號(hào)特征完善NET-VISA軟件，進(jìn)一步提高事件關(guān)聯(lián)的質(zhì)量。