基于深度學(xué)習(xí)到時(shí)拾取自動(dòng)構(gòu)建長(zhǎng)寧地震前震目錄

趙明， 唐淋， 陳石， 蘇金蓉， 張淼

1 中國(guó)地震局地球物理研究所， 北京 100081 2 北京白家疃國(guó)家地球科學(xué)野外觀測(cè)研究站， 北京 100095 3 四川省地震局， 成都 610041 4 加拿大戴爾豪斯大學(xué)地球和環(huán)境科學(xué)系， 新斯科舍省 哈利法克斯 B3H 4R2

0 引言

微震的識(shí)別和定位一直以來(lái)是地震研究的熱點(diǎn)和難點(diǎn).其中，有效地識(shí)別前震活動(dòng)和建立完備的前震目錄，對(duì)于理解前震和震群活動(dòng)歷程，地震觸發(fā)和地震成核機(jī)制等十分關(guān)鍵(Ross et al., 2019a).近年來(lái)，隨著觀測(cè)技術(shù)和儀器精度的提高，密集臺(tái)陣的大規(guī)模布設(shè)，地震監(jiān)測(cè)能力也有了很大提升，但與此同時(shí)需要處理的地震數(shù)據(jù)量不斷增長(zhǎng)，對(duì)高精度、可實(shí)用的自動(dòng)化數(shù)據(jù)處理方法和手段的需求越來(lái)越迫切.

地震數(shù)據(jù)處理一般包含挑取震相到時(shí)、多臺(tái)震相關(guān)聯(lián)、事件定位、測(cè)定震級(jí)等幾個(gè)環(huán)節(jié)，最終以地震目錄的形式發(fā)布.挑取震相、震相關(guān)聯(lián)等是后續(xù)地震定位和測(cè)定震級(jí)的基礎(chǔ)和前提，目前很大程度上仍依賴于人工處理.人工處理的主要挑戰(zhàn)在小震震相的識(shí)別和到時(shí)挑取上，小震事件一是數(shù)量多，根據(jù)Gutenberg-Richter經(jīng)驗(yàn)法則(Gutenberg and Richter, 1944)，地震每下降一個(gè)震級(jí)，數(shù)量會(huì)增加約10倍；二是信噪比相對(duì)較弱，很容易被誤檢和漏檢.多年以來(lái)人們發(fā)展了很多算法進(jìn)行小震的自動(dòng)識(shí)別與檢測(cè).其中，基于波形相似性的模板匹配算法在小震的檢測(cè)識(shí)別上十分有效(Peng and Zhao, 2009; Shelly et al., 2007; Yang et al., 2009; Zhang and Wen, 2015)，但其計(jì)算量較大，不能滿足實(shí)際應(yīng)用中高效快速處理大量連續(xù)波形數(shù)據(jù)的需求.

近年來(lái)，深度學(xué)習(xí)和人工智能技術(shù)的快速發(fā)展，在地震學(xué)領(lǐng)域也得到一系列成功應(yīng)用，主要體現(xiàn)在小震自動(dòng)識(shí)別(趙明等, 2019a; Perol et al., 2018; Yang et al., 2020)、震相挑取(趙明等, 2019b; Ross et al., 2018; Wang et al., 2019; 劉芳等, 2020)、震相關(guān)聯(lián)(Ross et al., 2019b)、定位(Zhang et al., 2020)等幾個(gè)方面，都有比較明顯的進(jìn)展.這些研究的共同點(diǎn)在于，將深度神經(jīng)網(wǎng)絡(luò)用于復(fù)雜地震波形抽象特征的自動(dòng)提取，在使用大量標(biāo)簽數(shù)據(jù)進(jìn)行監(jiān)督學(xué)習(xí)之后，訓(xùn)練好的模型通常在一些與訓(xùn)練數(shù)據(jù)相似的測(cè)試數(shù)據(jù)集或研究區(qū)域具有可媲美甚至超過(guò)傳統(tǒng)方法的效果.不過(guò)，深度神經(jīng)網(wǎng)絡(luò)方法實(shí)用化的關(guān)鍵在于模型能否達(dá)到比較高的泛化能力，即是否能在訓(xùn)練樣本以外的數(shù)據(jù)上也取得不錯(cuò)的效果.得益于海量高質(zhì)量訓(xùn)練數(shù)據(jù)樣本的收集與整理，一部分深度學(xué)習(xí)算法模型取得了較高的泛化能力.例如使用北加州地震臺(tái)網(wǎng)889個(gè)臺(tái)站近30年的事件波形數(shù)據(jù)(超過(guò)700000個(gè)帶P、S震相標(biāo)簽的樣本)訓(xùn)練得到的U型卷積神經(jīng)網(wǎng)絡(luò)P、S震相挑取算法(簡(jiǎn)稱PhaseNet)，適用于100 km震中距范圍內(nèi)的近震事件震相挑取(Zhu and Beroza, 2019).目前PhaseNet與傳統(tǒng)地震定位流程結(jié)合已經(jīng)被用于直接處理多個(gè)地區(qū)的連續(xù)地震波形，表現(xiàn)出出色的泛化能力，比如：2019年加州Ridgecrest地震序列(Liu et al., 2020)，2010年Guy-Greenbrier地震序列的誘發(fā)地震(Park et al., 2020)，以及短周期密集臺(tái)陣記錄的誘發(fā)地震(Wang et al., 2020a).

本研究借鑒了Liu等(2020)處理2019年加州Ridgecrest地震序列的方法和思路，將深度神經(jīng)網(wǎng)絡(luò)震相自動(dòng)挑取技術(shù)、最新發(fā)展的震相自動(dòng)關(guān)聯(lián)技術(shù)與發(fā)展相對(duì)成熟的傳統(tǒng)地震定位方法結(jié)合起來(lái)，構(gòu)建了一套直接從處理連續(xù)地震波形出發(fā)到產(chǎn)出地震目錄的自動(dòng)化流程.其具體步驟包括：(1)通過(guò)遷移學(xué)習(xí)將具備較高泛化能力的深度學(xué)習(xí)模型應(yīng)用至目標(biāo)研究區(qū)域的P、S震相自動(dòng)挑取；(2)應(yīng)用震相關(guān)聯(lián)技術(shù)自動(dòng)關(guān)聯(lián)P、S震相并進(jìn)行初步定位；(3)依次利用傳統(tǒng)絕對(duì)定位法和相對(duì)定位法在初步定位基礎(chǔ)上進(jìn)一步優(yōu)化定位結(jié)果，得到事件的發(fā)震時(shí)刻、經(jīng)緯度和深度; (4)自動(dòng)仿真并量取關(guān)聯(lián)好的S波形振幅并計(jì)算ML地方震級(jí).

基于上述自動(dòng)化數(shù)據(jù)處理流程，我們以2019年6月17日四川宜賓市長(zhǎng)寧縣發(fā)生的MS6.0地震震前小震活動(dòng)為研究對(duì)象，構(gòu)建對(duì)應(yīng)的自動(dòng)識(shí)別地震目錄.根據(jù)中國(guó)地震臺(tái)網(wǎng)發(fā)布的公開(kāi)地震目錄，長(zhǎng)寧地震震前半個(gè)月只有小于3級(jí)的小震發(fā)生，所有小震均集中沿長(zhǎng)寧背斜構(gòu)造線發(fā)生(長(zhǎng)約25 km，寬約5 km)，并且在長(zhǎng)寧地震震中50 km范圍內(nèi)有相對(duì)密集的臺(tái)站覆蓋(易桂喜等,2019; Lei et al.,2020)，這些條件對(duì)于驗(yàn)證和應(yīng)用我們的自動(dòng)化數(shù)據(jù)處理與地震目錄構(gòu)建流程是非常理想的.我們將得到的自動(dòng)識(shí)別目錄與公開(kāi)目錄進(jìn)行了對(duì)比，驗(yàn)證了自動(dòng)識(shí)別目錄的精度和可靠性，以及在小震識(shí)別方面的突出優(yōu)勢(shì).研究成果可以為長(zhǎng)寧地震震前微震活動(dòng)的快速特征分析提供參考.

1 臺(tái)站和數(shù)據(jù)

本研究所使用的連續(xù)波形數(shù)據(jù)來(lái)自2019年6月1日00時(shí)00分—6月17日22時(shí)55分期間，長(zhǎng)寧MS6.0地震震中75 km范圍內(nèi)21個(gè)臺(tái)站記錄的原始數(shù)據(jù)(圖1)，其數(shù)據(jù)主要有50 Hz(10臺(tái)短周期儀器)和100 Hz(11臺(tái)寬頻帶儀器)兩種，總計(jì)約28 GB.臺(tái)站對(duì)震區(qū)形成了較好的方位覆蓋，超過(guò)一半的臺(tái)(12個(gè))位于離震中30 km范圍內(nèi)，保證了定位結(jié)果的可靠性.按照深度學(xué)習(xí)模型對(duì)輸入數(shù)據(jù)的要求，我們對(duì)數(shù)據(jù)進(jìn)行了以下統(tǒng)一預(yù)處理：(1)去均值和去線性化；(2)統(tǒng)一采樣至100 Hz；(3)使用1 Hz的高通濾波器進(jìn)行濾波(提高信噪比).

本研究還收集了2019年四川臺(tái)網(wǎng)5—8月間的人工拾取震相報(bào)告，包括與16595個(gè)地震事件關(guān)聯(lián)的120233對(duì)P、S震相.我們根據(jù)震相報(bào)告和對(duì)應(yīng)的事件波形制作了用于遷移學(xué)習(xí)的訓(xùn)練集和驗(yàn)證集，其中以我們所研究的長(zhǎng)寧鹽礦區(qū)19個(gè)臺(tái)站(27.5°N—29.4°N，104°E—106°E)在2019年6月1日—30日期間記錄的32535條波形為驗(yàn)證集，其余103638條波形為訓(xùn)練集.

此外，長(zhǎng)寧MS6.0地震震前(2019年6月1日—6月17日22時(shí)55分)在研究區(qū)域內(nèi)人工檢測(cè)到了101個(gè)事件(ML0～2.9)，震相報(bào)告中與之關(guān)聯(lián)的P震相有474個(gè)，S震相有553個(gè)，這部分?jǐn)?shù)據(jù)在下文用于不同結(jié)果的對(duì)比.以上所有數(shù)據(jù)均由四川省地震局提供.

2 方法

本文自動(dòng)構(gòu)建地震目錄的流程分為四步，最終得到的是地震事件的發(fā)震時(shí)刻，經(jīng)緯度及深度位置，震級(jí)，以及與事件相關(guān)聯(lián)的P、S震相到時(shí)信息.具體流程細(xì)節(jié)如下.

2.1 基于深度學(xué)習(xí)模型的震相自動(dòng)識(shí)別

首先在PhaseNet算法以及利用北加州數(shù)據(jù)訓(xùn)練的模型(以下簡(jiǎn)稱北加州模型)基礎(chǔ)上，加上第1節(jié)所述的訓(xùn)練集和驗(yàn)證集進(jìn)行遷移學(xué)習(xí)訓(xùn)練，得到更適用于四川地區(qū)數(shù)據(jù)的檢測(cè)模型.由于新加入的數(shù)據(jù)相比北加州模型所用的訓(xùn)練數(shù)據(jù)是比較少的，所以我們?cè)谟?xùn)練過(guò)程中采用了一些數(shù)據(jù)增強(qiáng)的技巧，包括隨機(jī)平移，加噪、事件疊加、噪聲疊加、重采樣等(Zhu et al., 2020).為了證明PhaseNet算法在長(zhǎng)寧區(qū)域的適用性，我們分別使用遷移學(xué)習(xí)訓(xùn)練之后的模型、北加州模型、RNN模型(Zhou et al., 2019)、STA/LTA方法對(duì)同一測(cè)試集進(jìn)行到時(shí)拾取，并給出了每種算法預(yù)測(cè)值與其標(biāo)簽值的誤差直方圖，如圖2所示，我們定義與標(biāo)簽的絕對(duì)值誤差在0.5 s以內(nèi)的拾取震例為真正例(TP)，由此可以計(jì)算出精度(Pr)、召回率(Re)、F1綜合評(píng)估值等深度學(xué)習(xí)算法評(píng)價(jià)指標(biāo)，用于衡量算法的震相檢測(cè)能力，此外誤差直方圖的均值(μ)、方差(σ)、絕對(duì)平均誤差(MAE)則能夠在統(tǒng)計(jì)意義上反映算法到時(shí)拾取的準(zhǔn)確性.對(duì)比結(jié)果表明，遷移學(xué)習(xí)模型和北加州模型在所有六項(xiàng)指標(biāo)上都顯著優(yōu)于RNN模型和STA/LTA算法，而遷移學(xué)習(xí)訓(xùn)練的模型則在北加州模型基礎(chǔ)上有所提升，尤其在S震相的拾取上，同時(shí)修正了北加州模型直接用于四川地區(qū)數(shù)據(jù)所帶來(lái)的系統(tǒng)誤差(圖2a，c中誤差分布明顯偏向正值區(qū)域，這說(shuō)明網(wǎng)絡(luò)更傾向于給出稍滯后于人工標(biāo)注的到時(shí)預(yù)測(cè)，而遷移學(xué)習(xí)之后模型的誤差分布更為均勻，沒(méi)有這一傾向性).

圖1 (a) 研究區(qū)域及臺(tái)站分布.方框?yàn)檠芯繀^(qū)域，也是長(zhǎng)寧地震及其余震序列發(fā)生地點(diǎn)；白色圓圈表示所用臺(tái)站的震中距最大范圍，即以方框中心為圓心半徑 75 km范圍內(nèi)；藍(lán)線為人工勘察斷層(雷興林等(Lei et al., 2019))，綠線為構(gòu)造線； (b) 圖(a)方框的放大圖，圓點(diǎn)為我們的HYPODD重定位目錄，其顏色深淺代表事件的發(fā)生時(shí)間(從2019年6月1日至17日)，紅色同心圓為鹽礦注水井，圖中兩個(gè)紅色虛線方框?yàn)榈卣鸺袇^(qū)域，標(biāo)記為1號(hào)和2號(hào)區(qū)域以 便于正文中討論.Fig.1 (a) Study area and station distribution. Black rectangle represents the zoomed-in area in (b), as well as the seismicity zone of the Changning earthquake sequence; White circle with radius 75 km denotes the largest station distance; Blue lines and green lines are the surveyed faults (Lei et al., 2019) and mapped tectonic lines, respectively; (b) Zoomed-in region in (a). Dots denote earthquake locations in our HYPODD catalog, which are colored by their origin time (from June 1 to 17, 2019). Red concentric circle indicates the water injection well for salt mining. Red rectangles represent two regions with large number of earthquakes, which will be discussed in the main text as regions #1 and #2.

圖2 不同算法(從上到下依次為：遷移學(xué)習(xí)PhaseNet模型，北加州PhaseNet模型，RNN模型， STA/LTA)的自動(dòng)震相拾取誤差分布圖，左邊四圖為P波，右邊四圖為S波Fig.2 Performances of different picking algorithms (from up to down: Transfer learning PhaseNet model, PhaseNet model- Northern California, RNN model, STA/LTA): the left four figures are for P picks and the right four are for S picks

圖4 (a) P波走時(shí)-震中距關(guān)系圖，其中紅點(diǎn)表示REAL關(guān)聯(lián)震相，藍(lán)點(diǎn)表示人工關(guān)聯(lián)震相； (b) S波走時(shí)-震中距關(guān)系圖；(1)(2)為(a)(b)中虛線圈標(biāo)注的“離群點(diǎn)”所對(duì)應(yīng)的波形，其中藍(lán)色代表P到時(shí)，紅色代表S到時(shí)，實(shí)線為機(jī)器 拾取，虛線為人工拾取Fig.4 (a) Travel time to Hypocenter distance curves for the associated P phases. Red and blue dots indicate associated P phase picks in our REAL catalog and the routine catalog, respectively; (b) Same as (a), except for the S picks; (1)(2) are the waveform examples corresponding to “outliers” marked in (a)(b), where blue represents P arrival, red represents S arrival, the solid line is machine learning picks, and the dashed line is manual picks

圖5 (a)人工地震目錄，(b) REAL初步定位目錄，(c) VELEST目錄和(d)HYPODD目錄 圓點(diǎn)顏色代表震源深度，大小代表震級(jí)，其中(a)圖中虛線圓圈內(nèi)是3.2節(jié)討論的定位有誤的ML2.0人工目錄事件， 紅色五角星代表長(zhǎng)寧MS6.0地震，紫線為人工勘察斷層，綠線為構(gòu)造線，紅色同心圓代表鹽礦注水井.Fig.5 (a) Routine catalog, (b) REAL catalog, (c) VELEST catalog, and (d) HYPODD catalog Dots represent earthquakes in different catalogs, which are colored by their depth and scaled by their magnitude. The dot within the dotted circle in (a) is the ML2.0 event that discussed in Section 3.2, which was located wrongly in routine catalog. Red pentagram represents the MS6.0 Changning mainshock; Purple and green lines denote the surveyed faults and mapped tectonic lines, respectively. Red concentric circle represents the water injection well for salt mining.

圖6 人工目錄與VELEST目錄的震源參數(shù)對(duì)比圖 (a) VELEST目錄震級(jí)(紅色)與人工計(jì)算的震級(jí)(藍(lán)色)對(duì)比； (b) 兩種目錄共同事件的震級(jí)差； (c) 兩種目錄共同事件的發(fā)震時(shí)刻差； (d) 兩種目錄共同事件的水平位置差.Fig.6 The comparison of source parameters between routine catalog and the VELEST catalog (a) The number distribution of different magnitudes in the routine catalog (blue) and the VELEST catalog (red); (b) The magnitude difference of the common events in the two catalogs; (c) The origin time difference of the common events in the two catalogs; (d) The horizontal location difference of the common events in the two catalogs.

圖7 VELEST和人工目錄中差異較大的共同事件的波形圖 其中(a)和(b)是兩個(gè)目錄中的一組共同事件，(c)和(d)是另外一組共同事件.(a)和(c)為人工目錄中位置和關(guān)聯(lián)震相， (b)和(d)為VELEST目錄中的位置和關(guān)聯(lián)震相.Fig.7 Waveforms of common events with large differences between the VELEST and routine catalog (a) and (b) are one group of common events, while (c) and (d) are another. (a) and (c) correspond to the locations and associated phases in routine catalog; (b) and (d) correspond to the locations and associated phases in VELEST catalog.

我們將驗(yàn)證集上表現(xiàn)最優(yōu)的模型(即遷移學(xué)習(xí)后的PhaseNet)用于連續(xù)地震波形的P、S震相自動(dòng)拾取：首先將預(yù)處理之后的連續(xù)波形每間隔15 s依次切割成30 s長(zhǎng)(3000采樣點(diǎn))的三分量地震波形片段，然后將這些部分重疊的波形片段輸入PhaseNet模型進(jìn)行預(yù)測(cè)，預(yù)測(cè)曲線如圖3a所示，其中紅藍(lán)色高斯概率分布的峰值分別對(duì)應(yīng)于P、S到時(shí).在配備有一塊NVIDIA GTX 1080ti GPU加速卡的DELL Precision-7920-Tower工作站上，對(duì) 21個(gè)臺(tái)站的連續(xù)波形數(shù)據(jù)進(jìn)行了掃描，一共自動(dòng)挑取了131976個(gè)P到時(shí)，124357個(gè)S候選到時(shí)，單臺(tái)一天的連續(xù)波形掃描約需37 s，總共耗時(shí)約3.6 h.

2.2 自動(dòng)關(guān)聯(lián)與初步定位

多臺(tái)震相關(guān)聯(lián)是定位的基礎(chǔ)，同時(shí)也能很大程度上排除一部分誤識(shí)別震相，以及不屬于研究范圍內(nèi)的震相.這一步采用了Zhang等(2019)開(kāi)發(fā)的REAL算法.REAL算法是一種三維網(wǎng)格搜索算法，可以將按時(shí)間順序排列的多臺(tái)P、S到時(shí)與特定事件相關(guān)聯(lián)，同時(shí)通過(guò)震相到時(shí)計(jì)數(shù)和與理論到時(shí)的走時(shí)殘差計(jì)算來(lái)定位(圖3b).其搜索范圍為以P波最早記錄的臺(tái)站為中心，水平向27.75 km，深度方向15 km，網(wǎng)格大小為2.2 km×2.2 km×2 km.只有關(guān)聯(lián)到的P、S到時(shí)個(gè)數(shù)達(dá)到一定閾值才被認(rèn)為是一個(gè)事件，同時(shí)到時(shí)計(jì)數(shù)最多的網(wǎng)格所在位置作為該事件的初步定位結(jié)果；如果有多個(gè)網(wǎng)格擁有著共同最大數(shù)目關(guān)聯(lián)震相，則其中走時(shí)殘差最小的網(wǎng)格為最優(yōu)位置.本研究選取了P到時(shí)不少于3個(gè)，P+S不少于6個(gè)作為事件觸發(fā)閾值.在定位模型方面，我們用到了易桂喜等(2019)文章中表1的長(zhǎng)寧地區(qū)一維速度模型.這一步在研究區(qū)域范圍內(nèi)初步檢測(cè)并定位了1067個(gè)事件，以及與事件相關(guān)聯(lián)的4358個(gè)P震相，5022個(gè)S震相.

2.3 絕對(duì)定位與相對(duì)定位

在REAL關(guān)聯(lián)與初步定位基礎(chǔ)上，我們先后采用VELEST(Kissling et al.,1994)和HYPODD(Waldhauser and Ellsworth, 2000)進(jìn)一步優(yōu)化定位結(jié)果.VELEST是一種基于最小二乘的高效絕對(duì)定位算法.HYPODD 通過(guò)引入地震對(duì)走時(shí)殘差的差能夠得到高精度的震源位置參數(shù).二者都是地震事件定位常用的算法.在VELEST絕對(duì)定位后，我們使用空隙角≤220°對(duì)定位結(jié)果進(jìn)行篩選，并排除結(jié)果中走時(shí)殘差≥0.3 s的事件，得到870個(gè)絕對(duì)定位事件.這一嚴(yán)格設(shè)置雖然會(huì)減少有效地震的數(shù)目，但是能保證我們自動(dòng)檢測(cè)的可靠性和定位的精度.然后進(jìn)一步使用HYPODD進(jìn)行高精度相對(duì)定位.在這一步，我們使用了嚴(yán)格的參數(shù)篩選，即：事件-臺(tái)站對(duì)距離不超過(guò)50 km，而且定位震相不低于8對(duì).總共產(chǎn)生416個(gè)相對(duì)定位事件.

2.4 計(jì)算ML震級(jí)

通常測(cè)定地方性震級(jí)時(shí)要將寬頻帶數(shù)字地震記錄仿真成短周期位移記錄，并利用S波最大振幅計(jì)算，其計(jì)算公式為ML=logAμ+R(Δ)，其中Aμ為水平分量的最大振幅，R為量規(guī)函數(shù)，可查表獲得(見(jiàn)唐淋等(2018)).我們使用四川省地震局提供的實(shí)時(shí)仿真程序?qū)?shù)據(jù)進(jìn)行了Wood-Anderson型儀器實(shí)時(shí)仿真，并自動(dòng)量取S波在N、E分量上的最大振幅并取平均，代入公式中即可自動(dòng)計(jì)算得到震級(jí).

3 結(jié)果分析

3.1 拾取到時(shí)可靠性

在震相關(guān)聯(lián)方面，如圖4所示，與事件成功關(guān)聯(lián)的P、S震相數(shù)目均是人工的9倍以上.人工和自動(dòng)關(guān)聯(lián)到時(shí)主要集中在震中距<60 km范圍，二者的走時(shí)-震中距曲線均沿一條直線分布，并且能夠很好的匹配.我們注意到<10 km范圍人工拾取到時(shí)偏離直線的“離群點(diǎn)”更多一些，這種情形基本都是由于人工拾取的S震相不夠準(zhǔn)確所致，如圖4(a,b)中圈出的兩對(duì)“離群點(diǎn)”所對(duì)應(yīng)的波形(1)(2)，可以看到機(jī)器和人工在P震相的拾取上基本一致，但S震相拾取差異較大.當(dāng)S-P到時(shí)差不到2 s時(shí)，即使很小的到時(shí)拾取誤差也會(huì)造成較大的定位誤差，從而表現(xiàn)為走時(shí)-震中距曲線上的“離群點(diǎn)”，而深度學(xué)習(xí)對(duì)這些<10 km的事件自動(dòng)挑取的到時(shí)更為準(zhǔn)確.

3.2 目錄對(duì)比

圖5展示了101個(gè)人工定位，1067個(gè)REAL初步定位，870個(gè)VELEST絕對(duì)定位以及416個(gè)HYPODD相對(duì)定位結(jié)果.VELEST目錄更注重完整性，包含了所有三臺(tái)及三臺(tái)以上的定位事件，而HYPODD定位結(jié)果則利用嚴(yán)格的參數(shù)篩選保留了高精度相對(duì)定位事件.

首先我們與人工目錄對(duì)比了發(fā)震時(shí)刻、水平位置、震級(jí)方面的差異，分別計(jì)算了三者的均值和標(biāo)準(zhǔn)差.其中，由于人工目錄的深度并不準(zhǔn)確，所以我們只比較了水平定位誤差.計(jì)算結(jié)果表明，總共101個(gè)人工目錄事件，可以在VELEST目錄找到85個(gè)事件與之匹配，這85個(gè)事件與人工目錄的平均發(fā)震時(shí)刻誤差為0.36±0.07 s，震級(jí)差為0.15±0.024級(jí)，位置差為1.45±0.028 km.對(duì)16個(gè)未能與人工目錄匹配的事件，經(jīng)過(guò)檢查發(fā)現(xiàn)，10個(gè)是因?yàn)樽邥r(shí)殘差過(guò)大(0.3～0.5 s)而被VELEST目錄排除，4個(gè)是因?yàn)槎ㄎ晃恢寐湓趫D5所示研究區(qū)域的邊界外而被排除，1個(gè)因?yàn)榕_(tái)站間隙角過(guò)大(260°)被排除，1個(gè)因?yàn)榘l(fā)震時(shí)刻與另外一個(gè)大震級(jí)事件過(guò)近(2 s)被排除.

相匹配事件的震級(jí)差主要分布在[-0.4級(jí),0.2級(jí)]區(qū)間(圖6b)，超過(guò)這一范圍事件不超過(guò)10個(gè)，經(jīng)過(guò)檢查發(fā)現(xiàn)一是由于這些事件所關(guān)聯(lián)的震相到時(shí)和臺(tái)站數(shù)目不一致，通常自動(dòng)關(guān)聯(lián)的震相和臺(tái)站更多一些，二是算法自動(dòng)量取的最大振幅值與臺(tái)網(wǎng)工作人員手動(dòng)量取的振幅不一致.發(fā)震時(shí)刻差95%以上分布在[-0.6 s,0.4 s]區(qū)間(圖6c)，水平位置差95%以上在4 km范圍以內(nèi)(圖6d).對(duì)少數(shù)發(fā)震時(shí)刻、水平位置差超過(guò)區(qū)間的事件，我們將其對(duì)應(yīng)的多臺(tái)關(guān)聯(lián)波形按震中距排列并標(biāo)注到時(shí)和發(fā)震時(shí)刻進(jìn)行檢查，發(fā)現(xiàn)大都是人工目錄的結(jié)果存在偏差，例如圖7a和7b分別為其中一個(gè)2.0級(jí)事件的人工和VELEST定位事件波形，該事件人工和VELEST的發(fā)震時(shí)刻相差1.1 s，水平位置相差有5 km，從單臺(tái)看大部分的P、S到時(shí)挑取是準(zhǔn)確的，但將這些到時(shí)按震中距排列以后，圖7a的到時(shí)與理論到時(shí)曲線不符，且離震中更遠(yuǎn)的XWE臺(tái)P波到時(shí)比離震中更近的SHB臺(tái)P波到時(shí)更早，明顯是不合理的，而在我們關(guān)聯(lián)定位結(jié)果中近臺(tái)多了MET、TLO，并且去掉了XWE臺(tái)，觀測(cè)到時(shí)與理論到時(shí)曲線符合得很好(圖7b)，這說(shuō)明VELEST確定的發(fā)震時(shí)刻和位置更準(zhǔn)確，而人工結(jié)果出現(xiàn)了較大的發(fā)震時(shí)刻和定位偏差.圖7c和7d所示事件位置差達(dá)到了10 km，但人工定位只有三臺(tái)六個(gè)震相，而VELEST關(guān)聯(lián)的臺(tái)站和震相數(shù)都更多，其定位結(jié)果更可靠.以上人工目錄中的異常事件表明，基于人工拾取到時(shí)和手動(dòng)定位很難完全避免主觀性因素，而VELEST、HYPODD目錄的自動(dòng)構(gòu)建過(guò)程遵循統(tǒng)一的標(biāo)準(zhǔn)，從而盡可能避免了人為不確定因素.

在震級(jí)分布和數(shù)量上，我們分<0級(jí)，0～0.9級(jí)，1.0～1.9級(jí)，≥2.0級(jí)四個(gè)區(qū)間統(tǒng)計(jì)，如圖6a,c所示，其中VELEST目錄的震級(jí)范圍為ML-1.1～2.3，而人工目錄的范圍是ML0～2.9.VELEST目錄在0～1.0級(jí)區(qū)間的事件數(shù)目是人工目錄的8倍，并且包含更多低于0級(jí)的事件，將微震檢測(cè)能力下限提升至-1級(jí)左右.由于我們的方法檢測(cè)出了遠(yuǎn)多于人工目錄的1級(jí)以下事件，為了評(píng)估這些新增事件的真實(shí)性，我們按定位走時(shí)殘差、關(guān)聯(lián)到時(shí)數(shù)、震級(jí)對(duì)這些事件進(jìn)行排序，并檢查那些定位走時(shí)殘差最大(0.2～0.29 s)、關(guān)聯(lián)到時(shí)數(shù)最少(6～7個(gè))、震級(jí)低的事件(<0.3級(jí))，這樣的事件一共有48個(gè).通過(guò)對(duì)事件關(guān)聯(lián)波形的人工檢查，我們發(fā)現(xiàn)只有不超過(guò)5個(gè)疑似事件因?yàn)閿?shù)據(jù)信噪比較低，肉眼只能確認(rèn)其最近的臺(tái)站記錄為地震信號(hào)，其他臺(tái)站難以判斷，其余事件從波形和關(guān)聯(lián)情況看均可自信地判定為真實(shí)事件.

3.3 前震時(shí)空活動(dòng)性

如圖8所示，從地震活動(dòng)的每日頻次統(tǒng)計(jì)來(lái)看，人工目錄的地震活動(dòng)性從月初的每日20余次逐漸降低為震前的10余次，2級(jí)以上的地震也是月初(1—9日)發(fā)生比較頻繁，震前一周僅有一次.VELEST目錄微震數(shù)目更多，也呈現(xiàn)相同趨勢(shì)，即隨時(shí)間推移地震活動(dòng)性逐步減弱，震前一周都比較平靜.不過(guò)值得注意的有兩點(diǎn)：一是16日07時(shí)27分48秒在[28.319°N ,104.872°E]發(fā)生一次2.2級(jí)地震(人工目錄測(cè)定為2.3級(jí))，而且這次地震無(wú)論時(shí)間和空間上均十分靠近長(zhǎng)寧MS6.0主震，有可能是一次前震；二是我們統(tǒng)計(jì)了每日不同時(shí)間段的地震發(fā)生頻次，如圖8a所示灰色背景標(biāo)注的是每日06∶00—18∶00的時(shí)間段(大致相當(dāng)于白天)，白色背景則對(duì)應(yīng)于18∶00到次日06∶00的時(shí)間段(晚上)，分別統(tǒng)計(jì)這兩個(gè)時(shí)間段的地震發(fā)生頻次可以看出，晚上的地震事件數(shù)量普遍高于白天，這可能是因?yàn)橥砩吓_(tái)站的背景噪聲小，間接提升了臺(tái)站的微震監(jiān)測(cè)能力.類似的地震活動(dòng)周期性現(xiàn)象在加拿大Raton盆地的誘發(fā)地震監(jiān)測(cè)上也有發(fā)現(xiàn)，作者分析后認(rèn)為可能是風(fēng)噪引起(Wang et al.,2019).

圖8 2019年6月1日—17日VELEST目錄(棕色)與人工地震目錄(青色)的對(duì)比 (a) 每日頻次圖; 其中白色和灰色背景分別對(duì)應(yīng)于時(shí)間段：00∶00—06∶00和18∶00—24∶00 (白)，06∶00—18∶00(灰)； (b)累積數(shù)圖; 五角星代表震級(jí)相對(duì)較大的事件(ML≥1.9).Fig.8 The comparison between the VELEST catalog (brown) and the routine catalog (cyan) from June 1 to June 17, 2019 (a) Daily frequency diagram where the white and gray backgrounds respectively correspond to the time period: 00∶00—06∶00 and 18∶00— 24∶00 (white), 06∶00—18∶00 (gray) and (b) the corresponding cumulative number. Pentagrams denote earthquakes with ML≥1.9.

從前震的空間分布特征看，無(wú)論是VELEST還是HYPODD目錄，都可以看到前震主要沿長(zhǎng)寧—雙河復(fù)式大背斜的東南端(圖1b中2號(hào)區(qū)域)、西北端(圖1b中1號(hào)區(qū)域)聚集，而長(zhǎng)寧MS6.0主震正是發(fā)生在2號(hào)區(qū)域的中心位置.從震源深度看(圖5c,d)，2號(hào)區(qū)域內(nèi)前震深度較淺，而1號(hào)區(qū)域深度更深，這與長(zhǎng)寧MS6.0地震發(fā)生以后的余震深度分布是一致的(如易桂喜等，2019).此外,根據(jù)震相關(guān)聯(lián)與定位結(jié)果的統(tǒng)計(jì)，在2號(hào)區(qū)域，離長(zhǎng)寧鹽礦區(qū)的注水井(圖1b中紅色同心圓)最近的四個(gè)臺(tái)，其中CJW臺(tái)參與了750個(gè)事件的定位，CNI臺(tái)695個(gè)，MET臺(tái)612個(gè)，GXA臺(tái)608個(gè)，這些事件絕大部分均為<1.0級(jí)的事件，揭示出主震震前小震活動(dòng)已經(jīng)非常頻繁.相比之下人工目錄不具備檢測(cè)這類小震事件的能力，因此顯得震前鹽礦區(qū)的地震活動(dòng)相對(duì)平靜.我們可以合理推斷，這些小震活動(dòng)可能與鹽井注水活動(dòng)緊密相關(guān).Wang等利用InSAR數(shù)據(jù)的反演結(jié)果也揭示了鹽礦的注水采鹽活動(dòng)與本次長(zhǎng)寧MS6.0地震序列存在明顯的因果關(guān)系(Wang et al., 2020b).

4 結(jié)論

隨著地震數(shù)據(jù)量的快速增長(zhǎng)，自動(dòng)化和智能化處理是未來(lái)發(fā)展的趨勢(shì)，其中最近幾年方興未艾的人工智能和大數(shù)據(jù)挖掘是驅(qū)動(dòng)這一變革的關(guān)鍵技術(shù)(Bergen et al., 2019).深度學(xué)習(xí)目前仍是人工智能發(fā)展最快、潛力最大的一個(gè)分支研究領(lǐng)域，已經(jīng)廣泛用于社會(huì)生產(chǎn)的方方面面；在引入防震減災(zāi)領(lǐng)域后，它在地震學(xué)大數(shù)據(jù)分析方面顯示出很強(qiáng)的普適性、高精度與高效率.同時(shí)，傳統(tǒng)基于到時(shí)的定位方法發(fā)展相對(duì)成熟，其精度比較可靠，計(jì)算量較小.因此，將深度學(xué)習(xí)到時(shí)拾取、震相關(guān)聯(lián)技術(shù)與傳統(tǒng)定位方法聯(lián)系起來(lái)，構(gòu)建一套自動(dòng)化的地震定位和目錄產(chǎn)出流程，是非常值得探索的方向.本研究應(yīng)用遷移學(xué)習(xí)技術(shù)，將目前已具有較強(qiáng)泛化能力的PhaseNet震相自動(dòng)識(shí)別算法利用本地?cái)?shù)據(jù)訓(xùn)練之后，成功應(yīng)用于長(zhǎng)寧鹽礦區(qū)MS6.0地震震前的連續(xù)波形震相識(shí)別，并結(jié)合REAL震相自動(dòng)關(guān)聯(lián)技術(shù)、VELEST和HYPODD傳統(tǒng)定位技術(shù)，得到了長(zhǎng)寧地震前震活動(dòng)的自動(dòng)化絕對(duì)定位與相對(duì)定位目錄.自動(dòng)化產(chǎn)出地震目錄的流程的時(shí)間消耗主要集中在前期參數(shù)和閾值的測(cè)試調(diào)整，其后期應(yīng)用避免了人為主觀因素，提高了數(shù)據(jù)處理的效率，還提高了微震定位精度和目錄的完整性.在本研究中，自動(dòng)化絕對(duì)定位目錄不僅能夠識(shí)別和定位85%以上的人工目錄事件，同時(shí)在<1.0級(jí)小震事件識(shí)別數(shù)是人工的8倍以上.在準(zhǔn)確度方面，對(duì)照人工結(jié)果，自動(dòng)目錄在發(fā)震時(shí)刻、震級(jí)、水平位置等主要指標(biāo)上均已非常接近人工目錄，可以用于長(zhǎng)寧地震震前微震活動(dòng)時(shí)空特征的快速分析.

完備的高精度地震目錄對(duì)于長(zhǎng)寧區(qū)域地震活動(dòng)性和長(zhǎng)寧MS6.0主震的發(fā)震機(jī)理研究至關(guān)重要.未來(lái)將當(dāng)?shù)貥?gòu)造地質(zhì)、注水活動(dòng)數(shù)據(jù)，和我們自動(dòng)確定的前震時(shí)空分布相結(jié)合，有望綜合理解這次強(qiáng)震的發(fā)震機(jī)制.

致謝四川省地震局為本研究提供了人工震相數(shù)據(jù)、人工目錄和連續(xù)波形數(shù)據(jù)，朱尉強(qiáng)博士分享了PhaseNet算法和在北加州數(shù)據(jù)集上訓(xùn)練得到的模型(https:∥github.com/wayneweiqiang/PhaseNet)，周一劍博士提供了用于對(duì)比的RNN模型和識(shí)別程序(https:∥github.com/YijianZhou/CDRP_TF)，中國(guó)地震局地球物理研究所房立華研究員對(duì)論文初稿提出了寶貴意見(jiàn)，并提供了長(zhǎng)寧地區(qū)斷層結(jié)構(gòu)數(shù)據(jù)，在此表示衷心感謝.