王 娟 劉俊民 李 健 邱宏茂 王曉明 商 杰 蓋 磊
(禁核試北京國家數據中心和北京放射性核素實驗室,北京 100085)
地震監測是全面禁止核試驗條約規定的4種核查手段中遠區監測最有效的方法。條約生效后,CTBTO國際監測系統(IMS)通過與全球通訊基礎設施(GCI)的安全衛星鏈接傳送位于全球170個地震臺站的數據至國際數據中心(IDC),其中50個基本地震臺站1周7天1天24小時不間斷地傳送數據,120個輔助地震臺站在需要時,由程序自動按照既定的規則申請傳送相應時間段的數據[1]。地震監測數據進入IDC的數據處理系統后,首先進行臺站處理:即檢測信號、識別震相、計算震相特性(到時、振幅、方位角、慢度等)。單臺臺站處理完成后,根據計算出的震相特征,利用臺網關聯算法將來自不同臺站的信號進行關聯,推測這些信號可能來自同一事件,最終形成自動處理公報。分析員以自動處理公報為基礎,經過詳細分析,添加自動關聯漏掉的事件,刪除虛假關聯的事件,最終形成審定事件公報(REB)。目前,IDC的REB公報也已經成為全球地震領域比較認可的公報之一[2]。想要提高自動處理事件關聯的有效率,降低事件漏檢率,提高自動處理公報質量,減小分析員的工作量,是IDC一直以來致力研究的方向之一。
目前在IDC數據處理系統中,臺網關聯算法采用全局關聯算法(GA),其采用了基于全球格點的搜索算法來關聯震相形成事件[3]。GA首先將全球按緯度、經度分為覆蓋整個地球表面相互重疊的若干個圓形格點,格點文件包括了格點位置信息以及與臺站間的關系,然后對地球上給定位置(格點中心點)的所有可能檢測組合進行全局搜索,組合形成許多假設的事件,隨后合并相同事件,剔除冗余和虛假事件,最終形成公報產品。GA在IDC運行20多年來,在波形數據的實時處理中表現出了良好的處理性能,但是GA還是漏掉了感興趣的震級范圍內將近三分之一的地震事件,且自動處理公報中大約一半的事件為虛假事件[4-5]。
隨著計算機性能的提高、大數據分析、數據挖掘及機器學習等先進數據處理技術的廣泛應用,如何利用海量的歷史地震事件的特征來提高事件關聯性能成為學者們研究的熱點和IDC關注的重點之一。2009年開始,美國加州大學伯克利分校(University of California,Berkeley)的科學家與IDC開始合作研究基于震相特征模型的網絡處理-垂直綜合地震分析(Network processing-Vertically Integrated Seismic Analysis,NET-VISA)軟件[4-5],旨在替代GA成為IDC下一代全球震相自動關聯和事件檢測定位軟件。本文介紹了基于貝葉斯推理方法的NET-VISA算法的基本原理及在IDC歷年測試中的進展,并對禁核試北京國家數據中心(NDC)參加IDC組織的阿爾法測試(Alpha Tester Group,ATG)中朝鮮 6 次核試驗的 NET-VISA關聯結果作了分析對比。
NET-VISA核心是基于貝葉斯框架的關聯算法,由基于物理的生成概率模型(GM)和啟發式推理算法(IA)組成。該方法首先建立一個綜合考慮了信號到時、方位角、慢度、事件震級、正確震相、誤檢信號和尾波震相等事件特征的全球尺度的概率模型,然后由基于貝葉斯方法設計的推理算法(IA)搜索具有最大后驗概率(MAP)的事件集,并在關聯到達的事件之間進行迭代,最終確定事件的時間、位置和大小[5-6]。
假設X為覆蓋所有可能地震事件集合的一個隨機變量,每個事件由時間、位置、深度和震級來定義,Y為臺網內全部地震臺站的所有可能的信號觀測值,Pθ(X)描述了參數化的事件先驗概率,P?(Y|X)描述了信號的傳播和測量概率(包括走時、吸收與散射、噪聲等)。對于給定的觀測值Y=y,事件的后驗分布P(X|Y=y)即為事件參數。則NET-VISA的事件模型可以由最基本的貝葉斯公式來描述:

在NET-VISA中,一個到達信號最終分為關聯事件的震相 Λ ,誤檢信號 ξ和尾波信號 κ,在給定觀測信號A的情況下,NET-VISA結合了事件假設e,震相Λ,誤檢信號ξ 和尾波信號 κ得到的聯合概率模型[3]:

推理算法的最終目的是得到使后驗概率達到最大的事件參數的解:

NET-VISA生成概率模型由先驗模型和似然函數模型構成,主要包括事件先驗模型、震相似然函數模型、誤檢信號模型和尾波信號模型等部分。NETVISA中只考慮了14種震相:即P、Pn、PKP、S、PKPbc、PcP、pP、Lg、PKPab、PKKPbc、Pg、Rg、Sn和 ScP的震相模型。事件先驗模型反映了事件頻度、地理分布、深度分布和震級分布等;震相似然函數模型根據統計的事件發生時的震相檢測概率和震相屬性特征而建立;誤檢信號模型描述了由噪聲產生的誤檢信號的特征在各臺站的分布;尾波信號模型描述了由于信號傳播的散射而帶來的尾波信號,減少其對事件關聯的影響,通過尾波信號與其觸發震相的延遲、尾波信號的方位角、慢度和振幅特征來構建模型[3-5]。
NET-VISA的輸入從信號檢測開始,并基于所有標為P震相的檢測建立種子事件。通過構建事件、改善到達、改善事件、刪除事件等步驟查找最佳事件集以解釋觀察到的檢測結果,NET-VISA在改善到達和改善事件之間進行迭代以優化事件。圖1為NETVISA生成事件的示意圖。

圖1 NET-VISA 事件生成窗口Fig.1 NET-VISA windows
由圖1可知,NET-VISA將連續的到達數據流劃分為不同的窗口。在構建事件窗口中利用每個獨立震相的相關特征參數構建種子事件,圖1中構建事件窗口的窗長為W,窗長從t0=0開始到t1=t0+W(實際應用中W取 30 min)。到達窗口從t0?MT到t2=t0+W+MT,其中MT是震相的最大走時(33.3 min 或2 000 s),在此窗口中的每個到達都要進行與事件窗口中的事件相關聯。對于每個檢測在檢測窗口內通過“檢測改善算法”后,與事件窗口中的事件相關聯,進行迭代以優化事件。NET-VISA采用定義的事件度量標準和檢測評分標準來計算檢測和事件的評分,事件得分小于1的事件將被刪除。構建事件窗口向前移動步長S后,在輸出窗口中的最終確定的事件將寫入到事件公報中。
2009年3月,美國加州大學伯克利分校的斯圖亞特·羅素(Stuart Russell)教授在其主持的 ISS09(International Scientific Studies Project 2009,國際科學研究項目2009)項目中首次提出了NET-VISA算法。后來Arora博士在2009年6月的ISS09大會上報告了有限數據集的測試結果。2009年11月,在ISS09技術會議上討論了基于3個月的地震數據研究結果,當時相對于SEL3雖有改善,但還沒有量的體現。2010年2月,Russell教授在籌委會第34次核查工作組會議上報告了NET-VISA原型系統在以3個月數據集為基礎的測試結果,相對于SEL3,系統性能提高了5.5%。表1為當時的統計數據,表中SEL3代表用全局關聯算法(GA)產生的自動處理公報事件,LEB是分析員審定后的事件。

表1 NET-VISA 原型系統 3 個月測試結果Table 1 The three months test results of NET-VISA prototype
隨后Arora博士于2012—2013年針對地震數據處理研究開發了軟件,2014年集成了水聲處理流程,2015年開始了對次聲數據的處理。在PTS的組織下,Arora博士與IDC的科學家在2015年以離線方式進行了1年的地震數據處理測試,2016—2017年以離線運行方式針對包含次聲數據處理的NETVISA進行了測試,2018年1月—2019年2月在IDC測試平臺上對NET-VISA進行了在線測試。測試結果表明,在一般情況下,NET-VISA產生的公報要比IDC自動處理產生的公報SEL3更加完整和精準,對于REB的覆蓋率可以提高約15%(80%—90%與65%—75%),不一致率降低約5%(45%—55%與50%—60%)。
2019年PTS針對NET-VISA軟件的測試成立了ATG,NDC參加了阿爾法測試全過程,順利完成了全部4項測試,其中第2項測試是對2017年9月3日朝鮮核試驗事件進行數據處理,我們進行了擴展測試,對朝鮮6次核試驗逐個進行了類似處理,對NET-VISA關聯結果進行了分析對比。
利用阿爾法測試中IDC提供的NET-VISA程序對朝鮮6次核試驗數據逐個進行了重關聯。為了方便在同等條件下與IDC/GA關聯算法比較,在臺站列表方面,選擇IDC/SEL3朝鮮核試驗事件關聯臺站表中全部地震臺站;在數據段方面,以REB事件前約10分鐘整分為起始時間,持續時長2小時,以保證全部可能震相可被關聯。經NET-VISA關聯后的朝鮮地區地震事件基本參數列于表2。其中第4、5、6次核試驗事件均被NET-VISA分裂成兩個事件;第6次核試驗事件過后約8分半鐘還有一次小震級事件。

表2 NET-VISA定位的歷次朝鮮核試驗相關事件基本參數列表Table 2 Basic parameters of events related to several North Korean nuclear tests located by NET-VISA
以IDC/REB為參考,考察SEL3和VSEL(NETVISA產生的公報)事件定時和定位精準程度和存在的問題。定位和定時誤差按以下公式計算:

其中,E代表誤差,MSE代表均方誤差。
SEL3和VSEL定時誤差分別列于表3,圖2是對應表3的定時誤差直方圖。定位誤差列于表4,圖3是對應表4的定位誤差直方圖。表3和表4中VSEL2是VSEL的分裂事件,直方圖中NetVISA1對應于VSEL,NetVISA2對應于VSEL2,用2017-09-03(1)和 2017-09-03(2)區別 2017 年 9月第 1 次事件與第2次事件。定時均方誤差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 事件定時誤差直方圖Fig.2 Histogram of timing error of SEL3 and VSEL events

表3 以 REB 為參考的 SEL3 和 VSEL 定時誤差(s)Table 3 SEL3 and VSEL timing errors (s) with REB as reference
從測試結果來看,NET-VISA軟件模塊運行基本平穩、可靠。NetVISA檢出了全部核試驗事件,但是把第4、5、6次核試驗事件均分裂成了兩次事件;事件定時和定位方面,NetVISA比IDC的GA稍差一點。
NET-VISA作為IDC下一代臺網關聯軟件,從2009年概念提出、立項、研發、測試,前后已經過了10年,目前正在對NET-VISA軟件進行發布前的用戶體驗測試,雖然IDC的在線和離線測試均表明NET-VISA產生的公報要比IDC目前使用的全球關聯方法生成的自動處理公報SEL3更加完整和精準,但是ATG測試組的報告仍然反應出不少問題,這些問題反饋到PTS,有助于對NET-VISA軟件進行改進及模型優化。下一步,IDC將結合波形整體特征及更豐富的歷史信號特征完善NET-VISA軟件,進一步提高事件關聯的質量。