卷積神經(jīng)網(wǎng)絡(luò)在遠(yuǎn)-近地震震相拾取中的應(yīng)用及模型解釋＊

申中寅吳慶舉

（中國北京 100081 中國地震局地球物理研究所）

引言

隨著地震觀測規(guī)模與密度的提高，大量涌現(xiàn)的地震數(shù)據(jù)對震相拾取算法的性能及可靠性提出了更高要求.基于地震信號對特定物理量的擾動(dòng)，震相拾取包括震相識別與到時(shí)測量兩部分.這些物理量包括能量（如長短窗均值比）（Allen，1982）、峰度（Saragiotiset al，2004）、赤池信息準(zhǔn)則（Sleeman，van Eck，1999）等指標(biāo).相較于傳統(tǒng)方法，卷積神經(jīng)網(wǎng)絡(luò)（convolutional neural network，縮寫為CNN）能充分利用地震波形的全局特征，具有良好的魯棒性和泛化能力（ 趙明等，2019a，b；Zhuet al，2019）.根據(jù) CNN 在走時(shí)測量的參與程度，相關(guān)算法依次為：① CNN識別震相，利用傳統(tǒng)方法測量到時(shí)（如CNN+dbshear，趙明等，2019a）；② CNN識別震相，利用得分函數(shù)測量到時(shí)（Zhuet al，2019）；③ 地震波形經(jīng)CNN直接輸出震相類型及其到時(shí)，如 U-net （Ronnebergeret al，2015；趙明等，2019b）和多任務(wù) CNN （李健等，2020）.從①到③，整個(gè)震相拾取流程趨近端對端模式（輸入原始波形，直接輸出震相到時(shí)），并伴隨模型和標(biāo)簽復(fù)雜度的增加.與此同時(shí)，對CNN結(jié)構(gòu)與訓(xùn)練的探索（于子葉等，2020；周本偉等，2020）為方法優(yōu)化及其內(nèi)在原理闡釋開辟了道路.

另一方面，對CNN透明化的嘗試推進(jìn)了相關(guān)解釋算法在查漏優(yōu)化、可靠性評估及逆向?qū)W習(xí)等環(huán)節(jié)的應(yīng)用（Selvarajuet al，2016）.在震相識別領(lǐng)域，包括CNN在內(nèi)的機(jī)器學(xué)習(xí)方法，當(dāng)前大致處于復(fù)現(xiàn)人工挑揀的水平，因此有必要深入了解模型的潛在缺陷及其震相判別機(jī)制.在眾多解釋算法中，平滑GradCAM++和類模型可視化（class model visualisation，縮寫為CMV）（Simonyanet al，2014）分別從具體個(gè)案和模型整體出發(fā)，定量評估了CNN的決策敏感區(qū)域以及特征復(fù)現(xiàn)能力.其中 GradCAM 系列（Selvarajuet al，2016；Chattopadhyayet al，2018；Omeizaet al，2019）整合了類激活映射（class activation mapping，縮寫為 CAM ）的類別敏感性（Zhouet al，2016）與導(dǎo)向背傳播的像素級分辨率（guided back-propagation，縮寫為 GBP）（Springenberget al，2014），而平滑 GradCAM++ （Omeizaet al，2019）又進(jìn)一步優(yōu)化了最末卷積層的梯度權(quán)重及其導(dǎo)向背傳播過程.CMV通過反向傳播損失函數(shù)合成各類別的高分輸入圖形，是評估模型特征提取能力的重要依據(jù).在不觸及模型及數(shù)據(jù)結(jié)構(gòu)的前提下，上述算法已成功用于CNN圖像處理模型的解釋.然而在震相識別領(lǐng)域，CNN解釋算法的應(yīng)用與實(shí)踐尚處起步階段.為此，本研究采用CMV與平滑GradCAM++分析訓(xùn)練得到的CNN模型，以考察其震相識別的可靠程度與內(nèi)在機(jī)制.

鑒于此，本文擬采用CNN識別震相的得分函數(shù)測量到時(shí)來拾取北京國家地球觀象臺(tái)（以下簡稱北京臺(tái)）記錄的遠(yuǎn)震（P，S）和近震（Pg，Sg）波形，以期考察輕量CNN（參數(shù)不超過1萬）在小樣本集（少于1萬）的表現(xiàn).此外，模型解釋算法也被用于評估訓(xùn)練所得CNN的性能及其決策機(jī)制.

通過有放回隨機(jī)抽樣構(gòu)造訓(xùn)練樣本：隨機(jī)抽取固定數(shù)目的樣本并確保各震相數(shù)量相同（每個(gè)震相取2 000個(gè)，樣本子集大小為5×2000＝10 000）；重復(fù)上述抽樣過程，生成子樣本（10次采樣生成10個(gè)子樣本）.子樣本集合作為整體構(gòu)成一個(gè)完整的輪次參與訓(xùn)練.統(tǒng)計(jì)結(jié)果表明，子樣本集合囊括了幾乎全部震相及大部分“噪聲”.較小規(guī)模（5×500＝2500）的測試樣本同樣由10個(gè)子集構(gòu)成.上述的隨機(jī)采樣方法在充分利用現(xiàn)有數(shù)據(jù)的同時(shí)，還有助于訓(xùn)練過程中局部最優(yōu)點(diǎn)的遍歷.

2 模型搭建與訓(xùn)練

2.1 CNN模型結(jié)構(gòu)

鑒于不同震相持續(xù)時(shí)間的差異，為模型輸入長（40 s）?短（10 s）窗三分量記錄，以應(yīng)對不同時(shí)間尺度的波形特征.模型的輸出向量依次由噪聲N和震相P，S，Pg和Sg的概率組成.隱藏層結(jié)構(gòu)參考Zhu和Beroza （2019），隨卷積-池化依次倍縮（倍增）特征層的維度（數(shù)目）.由于激活函數(shù)選用ReLU，卷積層權(quán)重的隨機(jī)初始化采用he_normal分布，以避免梯度傳播的不穩(wěn)定現(xiàn)象（Heet al，2015）.相關(guān)論述可參閱附件2.研究還進(jìn)一步考察殘差塊添加對CNN性能的影響.殘差塊采用“瓶頸（bottle-neck）”結(jié)構(gòu)以控制參數(shù)規(guī)模（Heet al，2016）.

CNN和瓶頸式殘差塊結(jié)構(gòu)如圖2所示.圖2a展示了不同深度的卷積-池化序列.以4層CNN 為例，輸入張量（1000采樣點(diǎn)，6道）在雙邊鏡像延長至（1024，6）后，歷經(jīng) 4次卷積-池化成為（8，64），最后展平并全連接至輸出層.其它深度的網(wǎng)絡(luò)結(jié)構(gòu)依次類推.由于卷積層尺寸在深度為7時(shí)倍縮至1，模型深度搜索止于7層.圖2b展示了瓶頸式殘差塊的內(nèi)部構(gòu)造：數(shù)據(jù)流在依次卷積減-增特征層后恢復(fù)初始形狀，并與自身對位相加輸出結(jié)果.

圖2 卷積神經(jīng)網(wǎng)絡(luò)（a）和瓶頸式殘差塊（b）結(jié)構(gòu)示意圖紅色虛線對應(yīng)于層深為4的CNN結(jié)構(gòu)Fig.2 The structure of the convolutional neural network （a） and a bottle-neck residual block （b）The red dotted line corresponds to the CNN structure with a depth of 4 layers

2.2 模型訓(xùn)練

模型訓(xùn)練采用早停及L2正則化控制過擬合.其中早停用于動(dòng)態(tài)判定訓(xùn)練迭代的終止：以驗(yàn)證集的準(zhǔn)確率（損失函數(shù)）為準(zhǔn)，當(dāng)驗(yàn)證集準(zhǔn)確率（損失函數(shù)）連續(xù)3次低于（高于）最優(yōu)值時(shí)，終止訓(xùn)練并輸出最大準(zhǔn)確率（最小損失）模型.研究主要考察最大準(zhǔn)確率模型.L2正則化則懲罰模型復(fù)雜程度，最終損失函數(shù)為

式中：wm和M為模型參數(shù)及其總數(shù)；λ是L2因子的權(quán)重，可由數(shù)值實(shí)驗(yàn)確定合理范圍.考慮到λ搜索和殘差塊添加的運(yùn)算成本，模型最優(yōu)結(jié)構(gòu)的搜索將分步進(jìn)行（結(jié)果見第3節(jié)）：

1） 從1至7依次增加卷積-池化層數(shù)目，考察網(wǎng)絡(luò)深度對震相識別能力的影響.考慮到模型初始化的隨機(jī)效應(yīng)，訓(xùn)練將重復(fù)10次進(jìn)行.不考慮L2正則化（λ＝0）.

2） 在1）最優(yōu)模型的基礎(chǔ)上，搜索合適的λ取值λ*，所得模型及λ*將用于后續(xù)測試.λ搜索由疏到密分兩步進(jìn)行.

3） 在 2）最優(yōu)模型的基礎(chǔ)上，分層依次添加殘差模塊（Heet al，2016）.當(dāng)模型準(zhǔn)確率連續(xù)三次小于同深度最優(yōu)結(jié)果時(shí)，終止當(dāng)前層位的殘差塊添加.

在步驟2）和3）的訓(xùn)練中，在既有最優(yōu)模型的基礎(chǔ)上進(jìn)行了參數(shù)微調(diào).這不僅節(jié)省訓(xùn)練時(shí)間，還充分利用了步驟1）的隨機(jī)搜索結(jié)果，有助于模型參數(shù)的客觀比較.

2.3 性能指標(biāo)

混淆矩陣（confusion matrix）記錄了 CNN 對測試集樣本不同震相的識別能力，是模型性能評估的重要手段.在圖3中，混淆矩陣的行和列序號分別對應(yīng)樣本的真實(shí)標(biāo)簽和預(yù)測結(jié)果：元素Cij即標(biāo)簽為i的樣本被分類為j的數(shù)目.一般而言，模型性能越好，對角線元素越大.為全面客觀地描述模型表現(xiàn)，混淆矩陣派生出包括準(zhǔn)確率、震相精度、召回率和F1得分等指標(biāo).

圖3 混淆矩陣Fig.3 Confusion matrix

準(zhǔn)確率A反映了模型的總體震相識別水平，是混淆矩陣對角線元素之和占樣本總數(shù)的比例.精度P和召回率R描述了模型對特定震相的識別能力.其中震相精度反映了預(yù)測類別的正確比例.召回率則關(guān)注指定標(biāo)簽的成功檢出率.模型準(zhǔn)確率、震相i的精度、震相i的召回率表達(dá)式分別為：

F1得分被用于考察模型的整體性能，該準(zhǔn)則傾向于為指標(biāo)均衡的模型賦予高分.當(dāng)精度與召回率同等權(quán)重時(shí)，震相i的F1得分定義為

而整體F1得分

則作為模型性能的最終衡量指標(biāo)，衡量震相識別的總體表現(xiàn).

3 訓(xùn)練結(jié)果

卷積層深度.未清洗樣本訓(xùn)練時(shí)不同模型深度N的準(zhǔn)確率A、損失函數(shù)和F1t如圖4a，b和附件1的表1所示.受模型初始化等隨機(jī)因素影響，各輪次結(jié)果波動(dòng)顯著.盡管如此，依舊可見以下規(guī)律：單卷積層顯著劣于其它結(jié)構(gòu)；當(dāng)卷積層深度大于等于2時(shí)，CNN準(zhǔn)確率維持較高水平；中等深度（4，5層）模型表現(xiàn)較好.其中4層CNN準(zhǔn)確率最高（A＝0.826，F(xiàn)1t＝82.4），而最低損失值則出現(xiàn)在5層網(wǎng)絡(luò).單層網(wǎng)絡(luò)表現(xiàn)較差主要源于過簡結(jié)構(gòu)的欠擬合，而深層（特別是7層）準(zhǔn)確率低的結(jié)構(gòu)成因有待論證.具體數(shù)值及統(tǒng)計(jì)結(jié)果可參見附件1的表1.

圖4 卷積層深度對模型性能的影響（a） 未清洗數(shù)據(jù)，第4輪訓(xùn)練中不同卷積層深度（線上序號）模型的表現(xiàn)；（b） 未清洗數(shù)據(jù)，各輪次訓(xùn)練中不同卷積層數(shù)模型的最高準(zhǔn)確率和最低損失函數(shù)；（c） 已清洗數(shù)據(jù)，第10輪訓(xùn)練中不同卷積層深度（線上序號）模型的表現(xiàn)；（d） 已清洗數(shù)據(jù)，各輪次訓(xùn)練中不同卷積層數(shù)模型的最高準(zhǔn)確率及最低損失值Fig.4 The influence of convolutional layer depth on model performance（a） Data unwashed，the model perfomance for different depths （marked by numbers） during the 4th training round；（b） Data unwashed，the maximum accuray and minimum loss with different depths in the total 10 training rounds；（c） Data washed，the model perfomance for different depths （marked by numbers） during the 10th training round；（d） Data washed，the maximum accuray and minimum loss with different depths in the total 10 training rounds

L2正則化.在最優(yōu)深度的基礎(chǔ)上（N＝4，A＝0.826），測試了L2正則化系數(shù)（λ）的影響，采樣間隔按疏密依次進(jìn)行（附件1中表2）.疏測試平行進(jìn)行3次，比較λ的不同量級：0.001，0.1，1，10，100.當(dāng)λ為10，100時(shí)，模型準(zhǔn)確率分別為0.842和0.838，明顯高于未正則化的結(jié)果.為獲得λ的最優(yōu)范圍，密測試從10到100進(jìn)行（步長為10）.由于結(jié)果穩(wěn)定性良好（最大偏差不超過0.005），密測試只進(jìn)行1輪.結(jié)果表明，當(dāng)λ＝30時(shí)模型準(zhǔn)確率最高（0.860）.

殘差塊添加.在L2正則化最優(yōu)模型（N＝4，λ＝30，A＝0.860）的基礎(chǔ)上逐層添加殘差塊.相比卷積層深度及L2正則化，殘差塊帶來的改善十分有限（模型準(zhǔn)確率＋0.002至＋0.005不等）.因此，本文不再考察殘差塊對CNN的影響.

數(shù)據(jù)清洗及重新訓(xùn)練.基于最優(yōu)模型（N＝4，λ＝30，A＝0.860），逐一審查訓(xùn)練集和測試集的錯(cuò)判案例，并清洗實(shí)際波形與標(biāo)簽對應(yīng)不佳的樣本.被清洗的數(shù)據(jù)主要來自：噪聲疊加，信噪比過低，震相分析標(biāo)注的偏差，噪聲自動(dòng)截取時(shí)相鄰震相的混入等.由于剔除了不合適的震相和噪聲疊加樣本（特別是P與S），單個(gè)訓(xùn)練集與驗(yàn)證集分別收縮至5×800和5×175.

給定λ＝30，卷積層深度N取1，2，···，7訓(xùn)練10次，結(jié)果見圖4c，d和附件1表1.得益于訓(xùn)練樣本波形的改善以及驗(yàn)證集質(zhì)量的提高，數(shù)據(jù)清洗使模型準(zhǔn)確率大幅提高.雖然最高準(zhǔn)確率出現(xiàn)在5層CNN （A＝0.971，F(xiàn)1t＝97.1），模型性能隨網(wǎng)絡(luò)深度變化的趨勢依然不變：適中的深度（N＝4，5，6）普遍優(yōu)于過淺、過深情形（N＝1，2，7）.與此同時(shí)，卷積層深度為5的CNN依然具有最低的損失值.

4 模型解釋

本節(jié)使用兩種模型解釋算法評估最優(yōu)CNN模型（清洗樣本集，N＝5，λ＝30，A＝0.971）.其中類模型可視化CMV旨在反演各震相的最高得分波形，平滑GradCAM++ （Omeizaet al，2019）則用于勾勒輸入波形的決策敏感區(qū)域.這些梯度算法能在不觸及模型結(jié)構(gòu)的前提下，定量評估目標(biāo)震相的特征提取及其判別機(jī)理.

4.1 類模型可視化（CMV）

作為基于模型本身的解釋方法，CMV根據(jù)得分函數(shù)梯度修正空白輸入（I0＝def0），以獲取使震相c得分yc最大的特征波形.為減少類別混染，yc取“softmax”激活之前的數(shù)值.為控制幅值及波動(dòng)水平，目標(biāo)函數(shù)具有形式

圖5 類模型可視化算法的偽代碼（Nguass＝20）Fig.5 Pseudo code of Class Model Visualisation（Ngauss＝20）

完整的CMV反演結(jié)果可見附件3，圖7a選取了其中的長窗分量，可見如下特征：

1） 短窗振幅顯著小于長窗；

2） P和Pg垂直分量優(yōu)勢明顯，Sg水平分量有所發(fā)育，S水平分量占優(yōu)；

3） 噪聲的較大振幅出現(xiàn)在窗口前半部分；

4） 地震震相波形各異：P與Pg振幅突增明顯，存在明顯后續(xù)能量；S受前驅(qū)干擾影響顯著，后半窗波形模糊；Sg則以清晰的Pg前驅(qū)能量為特征.

其中2）——4）表明訓(xùn)練所得CNN模型能捕捉到目標(biāo)震相的基本特征，但1）暗示短窗數(shù)據(jù)未能有效參與預(yù)測.此外，震相S呈現(xiàn)出較低的波形質(zhì)量，震相Sg波形明顯偏離窗口中心位置.這些都是CNN應(yīng)用中必須注意的問題，也是未來需要改進(jìn)的方向.

4.2 平滑GradCAM++

GradCAM++旨在評估波形不同部位對模型預(yù)測的影響，算法可分解為“Grad”和“CAM”兩部分.其中Grad通過誤差量的導(dǎo)向背傳播（guided back-propagation，縮寫為GBP）描繪圖像細(xì)節(jié)對決策函數(shù)的貢獻(xiàn).而CAM則反映最末卷積層不同位置對分類結(jié)果的影響.CAM結(jié)果經(jīng)上采樣與Grad對位相乘，得到輸入波形對模型決策的敏感度分布.相關(guān)公式可見Chattopadhyay等（2018），圖6為平滑 GradCAM++算法的偽代碼.Omeiza等（2019）的平滑GradCAM++藉由加噪平均進(jìn)一步提高了模型解釋水平，也是本文采用的算法.

圖6 平滑 GradCAM++的偽代碼Fig.6 Pseudo code of smooth GradCAM++

為獲得震相敏感區(qū)域，平滑GradCAM++只選取高分案例（得分＞0.99）.相比二維圖像處理，震相識別更易出現(xiàn)接近1.0的得分.此時(shí)過小的損失函數(shù)會(huì)隨GBP消失為零，導(dǎo)致空白結(jié)果.噪聲疊加可有效防止過高得分引起的GBP消失，這也是選用平滑GradCAM++的重要理由.研究向輸入波形添加高斯噪聲（1/2最大振幅）并平均所得GBP （10次），并與CMV一同構(gòu)成了CNN解釋的依據(jù).

圖7b展示了震相波形（綠色）及其平滑GradCAM++結(jié)果（灰色），完整內(nèi)容可附件3.與CMV類似，平滑GradCAM++的短窗振幅普遍偏低.長窗垂直分量在模型預(yù)測中起著決定性的作用.P與Pg對窗口中心兩側(cè)較為敏感：前側(cè)的平靜與后側(cè)的振動(dòng)均對震相判別不可或缺.Sg也具有前側(cè)（Pg振動(dòng)）和后側(cè)（Sg振動(dòng)）兩個(gè)敏感區(qū)域，只是后者幅度相對較小.對于震相S，窗口后半段并未出現(xiàn)可見的敏感區(qū)域.結(jié)合CMV （4.1節(jié)）結(jié)果可以看出，模型在波形分量平衡及震相S識別等方面尚有較大改進(jìn)空間.

圖7 模型解釋結(jié)果（a） 類模型可視化（CMV）；（b） 平滑 G radCAM++.綠色背景為原始波形Fig.7 Results of model interpretation（a） Class model visualization （CMV）； （b） Smooth GradCAM++.The green background stands for the original waveforms

5 連續(xù)波形下的模型部署

參照Zhu等 （2019）將2020年4月的連續(xù)波形用于震相拾取實(shí)測，以規(guī)避訓(xùn)練集和測試集的影響.只有當(dāng)震相（P，S，Pg和Sg）概率大于閾值時(shí)才認(rèn)定為有效觸發(fā)，得分函數(shù)最大時(shí)刻即為震相測量到時(shí).為平衡掃描效率和精度，研究采用變步長滑動(dòng)窗口進(jìn)行震相的搜索和測量，窗口長度以20 s 長窗為準(zhǔn).特征函數(shù)選用模型的輸出得分.具體流程如下：

1） 長步長（1 s）搜索震相，合并相鄰?fù)鹣啻翱冢?/p>

2） 短步長（0.1 s）計(jì)算目標(biāo)震相的特征函數(shù)，最大值對應(yīng)時(shí)刻即震相的測量到時(shí).

在連續(xù)波形掃描之前（5.2節(jié)），有必要考察震相片段的局部掃描結(jié)果（5.1節(jié)），以便合理選定閾值.

5.1 震相片段的掃描

為確定不同震相閾值，首先對波形片段（基于測震目錄）進(jìn)行短步長掃描，概率得分及走時(shí)偏差如圖8和附件1表3 （前4行）所示.不同震相表現(xiàn)為：震相P，具有大量較低得分，走時(shí)偏差集中在0 s附近；震相S，得分與走時(shí)偏差均顯著分散；震相Pg，得分較高（＞0.999），除過個(gè)別情況（走時(shí)偏差＞10 s，與Sg混淆），走時(shí)測量略有提前；震相Sg，得分較高（＞0.999），走時(shí)整體延遲.其中Pg-Sg走時(shí)偏差印證了CMV反演波形中震相初動(dòng)相對窗口中心位置的偏離.為此，我們以0.999作為短步長測量的得分閾值，以平衡震相拾取的準(zhǔn)確性與召回率.考慮到模型解釋中震相S的較差表現(xiàn)（圖7），震相S的得分閾值被調(diào)高至0.999 999以屏蔽干擾.為避免稀疏采樣引起的漏檢，長步長（1 s）搜索采用了更為寬松的閾值（P＞0.9，S＞0.99，Pg＞0.9，Sg＞0.9）.本研究設(shè)定的閾值（0.999）明顯高于 Zhu等 （2019）的0.5和Ross等（2016）的0.98.這可能來自樣本及訓(xùn)練等環(huán)節(jié)，具體可見6.3節(jié).

圖8 連續(xù)觀測震相片段的走時(shí)測量結(jié)果（a） 各震相走時(shí)偏差及模型得分的分布；（b） 各閾值下諸震相的識別數(shù)目（藍(lán)色方點(diǎn)）與走時(shí)偏差（紅色圓點(diǎn)）Fig.8 Traveltime mearurements for seismic phase clips cut from continuous observation（a） Travetime residuals and CNN-predicted phase scores；（b） Detection number （blue square） and traveltime residuals （red dot） of the phases

5.2 連續(xù)波形掃描

根據(jù)預(yù)設(shè)閾值掃描了2020年4月的連續(xù)波形.鑒于Pg-Sg總成對出現(xiàn)且到時(shí)間隔有限（小于 6 0 s），事先修剪了未成對（或間隔超過 6 0 s）的Pg，Sg拾取記錄.當(dāng)震相拾取到時(shí)與目錄到時(shí)相差小于 5 s時(shí)才定為檢出（S 震相放寬至 1 0 s），結(jié)果見圖9和附件1表3 （最末行）.可以看出，震相間的干擾十分有限（僅有1例Sg誤識別為Pg）.識別偏差主要來自CNN的噪聲過敏感（噪聲被誤判為震相）與震相欠敏感（震相被誤判為噪聲）.正確識別結(jié)果的得分函數(shù)分布如圖10所示.為突出相對變化 趨勢，得分函數(shù)經(jīng)過了單調(diào)變換x→ ? lg（1－x）、去均值和最大值歸一等處理.其中震相P僅表現(xiàn)為高頻成分的陡增，卻依然被準(zhǔn)確拾取.而波形清晰的Pg，Sg卻呈現(xiàn)出明顯的走時(shí)偏差.對于錯(cuò)誤結(jié)果，本節(jié)將從召回率和準(zhǔn)確率兩方面分別討論.

圖9 連續(xù)波形掃描結(jié)果的混淆矩陣Fig.9 The confusion matrix for the continuous waveform scanning

圖10 正確的震相識別水平虛線由低到高依次對應(yīng)長、短步長的得分閾值，原始波形經(jīng)歸一化處理，下同F(xiàn)ig.10 Seismic phases detected correctlyHorizontal dotted lines correspond to thresholds for long（lower） and short （upper） scanning steps，the waveform is normalized，the same below

震相召回率反映了CNN的靈敏程度.震相未能召回的原因主要有：震相概率得分過低，走時(shí)測量偏差過大.不同震相的召回率分別為：P （66.7%），S （10%），Pg （40%），Sg （60%）.對于震相P，較低的信噪比壓低了概率得分，是漏檢的主要原因（大于20例，圖11b）.而波形清晰的漏檢案例常具有豐富的低頻成分，并與較高的得分閾值相關(guān)（圖11a）.除了較高的閾值，震相S的召回率也受波形質(zhì)量與模型性能制約（圖11c?d）.噪聲擾動(dòng)（導(dǎo)致過低得分）及走時(shí)拾取的系統(tǒng)偏差是Pg-Sg漏檢的重要原因（圖11e?f）.

震相精度反映了CNN的抗噪能力.不同震相精度分別為：P （29.3%），S （41.7%），Pg （6.8%），Sg （10.2%）.較低的震相精度表明大量噪聲波形被誤判為地震震相.對于震相P，除過個(gè)別波形失真、疑似地震（圖11g）、震相誤讀（圖11i）外，誤判案例主要來自周期為4——5 s的背景噪聲的擾動(dòng)（總數(shù)不小于150，圖11j）.此外，得分函數(shù)振蕩造成的搜索窗口分裂也是導(dǎo)致誤判的原因之一（圖11h）.震相S較高的得分閾值在造成大量漏檢的同時(shí)，也有效降低了噪聲混入的概率.除過個(gè)別震相誤判（圖11k），震相S的假正例主要來自模型對背景噪聲的過敏感（圖11l）.后者也是Pg-Sg假正例的主要來源（圖11n）.與此同時(shí)，模型亦成功拾取了未錄入測震目錄的Pg-Sg事件（圖11m）.

圖11 錯(cuò)誤震相的識別（e） Pg 走時(shí)測量偏差過大；（f） Pg-Sg 波形不清晰；（g） 疑似地震事件；（h） 搜索窗口分裂；（i） 后續(xù)震相（ScP）干擾；（j） 常見 P震相誤判波形；（k） 誤判為 S 的 S cS 震相；（l） S 震相的常見誤判波形；（m） 未錄入地震目錄的 P g-Sg；（n） Pg-Sg 常 見誤判波形Fig.11 Cases of the wrong detections（e） Too enormous Pg travetime residual；（f） Unclear Pg-Sg waveform；（g） Suspisious earthquake；（h） Splitting in searching window；（i） Later-coming phase （ScP）；（j） Common false P detection；（k） ScS identified as S；（l） Common false S detection； （m） Pg-Sg event not in catalogue；（n） Common false Pg-Sg detection

圖11 錯(cuò)誤震相的識別（a） P 震 相波形清晰；（b） P 震 相波形不清晰；（c） S 震 相波形清晰；（d） S 震 相波形不清晰Fig.11 Cases of the wrong detections（a） Clear P waveform；（b） Unclear P waveform；（c） Clear S waveform；（d） Unclear S waveform

6 討論

上文介紹了樣本制備、訓(xùn)練以及模型解釋、部署的整個(gè)流程，此節(jié)將重點(diǎn)說明其中涉及的一些問題.

6.1 訓(xùn)練樣本規(guī)模

雖然模型的參數(shù)量被控制在較低的量級，其復(fù)雜度依然高于訓(xùn)練樣本數(shù)目.為此我們對比了震相樣本規(guī)模從200到800時(shí)各深度模型的訓(xùn)練結(jié)果，訓(xùn)練重復(fù)5次.

圖12a展示了CNN準(zhǔn)確率隨訓(xùn)練樣本規(guī)模的變化趨勢.4層CNN的準(zhǔn)確率普遍較高，5層次之.隨著樣本規(guī)模的增加，深層（6，7）模型的表現(xiàn)逐漸優(yōu)于淺層模型（2，3）.圖12b對應(yīng)樣本實(shí)際使用量與訓(xùn)練樣本規(guī)模之間的關(guān)系.受樣本池規(guī)模所限，樣本使用量的增速隨樣本規(guī)模的增加持續(xù)減小.這雖然會(huì)影響模型性能的提高，但無法解釋圖12a呈現(xiàn)的復(fù)雜趨勢.

圖12 樣本規(guī)模對模型準(zhǔn)確率及樣本利用情況的影響（a）模型準(zhǔn)確率與訓(xùn)練樣本規(guī)模的關(guān)系；（b）訓(xùn)練樣本規(guī)模與樣本使用量的關(guān)系Fig.12 Relationship among train data size，model accuracy，and samples used（a） Relationship between train data size and model accuracy （b） Relationship between train data size and the number of samples used

圖13列舉了模型準(zhǔn)確率隨深度的變化，是圖12a的詳細(xì)補(bǔ)充.當(dāng)訓(xùn)練樣本規(guī)模小于500時(shí)，模型準(zhǔn)確率與樣本規(guī)模正向相關(guān).而當(dāng)樣本規(guī)模超過500時(shí)，模型準(zhǔn)確率增速明顯放緩（3，4，6層）乃至下降（2，5，7層）.除淺層（2，3層）外，模型最高準(zhǔn)確率均未對應(yīng)最大樣本規(guī)模.此外，即使樣本規(guī)模僅為200，4層CNN的準(zhǔn)確率依然超過了0.96，反映了輕量CNN對小規(guī)模訓(xùn)練樣本的較強(qiáng)適應(yīng)能力（Zhuet al，2019）.

圖13 各深度模型準(zhǔn)確率與訓(xùn)練樣本規(guī)模的關(guān)系Fig.13 Relationship between the raining data size and model accuracy

由此可見，擴(kuò)充樣本有助于模型準(zhǔn)確率的提高，當(dāng)前樣本規(guī)模尚足以滿足模型訓(xùn)練需要.同時(shí)，訓(xùn)練樣本規(guī)模亦無法解釋深層（6，7層）CNN準(zhǔn)確率的下降.

6.2 正則化系數(shù)

在CNN的深度探索階段，研究采用了固定的正則化系數(shù)λ.為考察λ取值對結(jié)果的影響，本節(jié)將系統(tǒng)考察最優(yōu)λ隨CNN深度的分布趨勢.

圖14a展示了模型準(zhǔn)確率隨λ的變化趨勢，訓(xùn)練進(jìn)行3次并選取最高準(zhǔn)確率結(jié)果.隨深度增加，最優(yōu)λ也相應(yīng)增大.即便考慮高λ值對深層CNN準(zhǔn)確率的提高，5層模型依舊具有最好的表現(xiàn).圖14b展示了λ與（模型系數(shù)的平方均值）的關(guān)系.由于正則項(xiàng)約束的增強(qiáng)，模型均方隨λ增加持續(xù)下降.而對于給定λ，模型系數(shù)的平方均值則與層數(shù)反向相關(guān)，體現(xiàn)了模型總體激活程度的減弱.最小二乘擬合結(jié)果揭示了最優(yōu)λ與模型系數(shù)均方滿足

表明不同深度的最優(yōu)模型具有大致相同的正則項(xiàng)（圖14c），其內(nèi)在機(jī)制有待進(jìn)一步澄清.

最優(yōu)λ隨模型深度的分布可解釋兩次深度搜索結(jié)果的不同（3.1節(jié)與3.4節(jié)）.根據(jù)圖14a，當(dāng)λ＜1時(shí)4層CNN具有明顯優(yōu)勢，對應(yīng)了第一輪深度搜索的結(jié)果（無正則項(xiàng)，相當(dāng)于λ＝0）.而當(dāng)進(jìn)行第二輪搜索時(shí)，λ＝30 （第一輪4層CNN的λ搜索結(jié)果）使5層CNN的性能大幅提升，并使后者成為最優(yōu)搜索模型.據(jù)此可知，當(dāng)λ＝80時(shí)（當(dāng)前取值30），5層網(wǎng)絡(luò)的性能還有望進(jìn)一步提升.

圖14 正則化對各深度模型的影響（a） 正則化系數(shù)λ、模型層數(shù)與模型準(zhǔn)確率之間的關(guān)系，圓點(diǎn)對應(yīng)最高模型準(zhǔn)確率；（b） 正則化系數(shù)λλ與模型均方之間的關(guān)系，其中空心圖形對應(yīng)最高準(zhǔn)確率模型（注意2層數(shù)據(jù)點(diǎn)疊覆于3層之下）；（c） 模型損失函數(shù)的構(gòu)成Fig.14 The effects of regularization on CNN with different depths（a） Relationship among regularization factor λ，CNN depth，and accuracy；（b） Relationship between regularization factor λ and the squaremean of model weights，with hollow patterns standing for CNN with highest accuracy （2-layer dot is overlapped by that of 3-layer）；（c） Variation of loss function with CNN depth

6.3 連續(xù)波形掃描的過高閾值

本文的震相得分閾值明顯高于同類研究的結(jié)果（Rosset al，2016；Zhuet al，2019），這雖然出于屏蔽噪聲干擾的需要，但也與震相本身的過高得分密不可分.對于后者，本節(jié)將從數(shù)據(jù)泄露和模型過擬合兩方面進(jìn)行討論.

數(shù)據(jù)泄露會(huì)造成虛高的訓(xùn)練精度，通常來自模型對個(gè)別特征過度依賴以及驗(yàn)證集同訓(xùn)練集的混染.前者常見于聚類分析的不當(dāng)參數(shù)選取，表現(xiàn)為個(gè)別指標(biāo)完全控制模型分類.作為基于波形的分類算法，CNN對震相特定結(jié)構(gòu)的過度敏感會(huì)導(dǎo)致CMV和平滑GradCAM++結(jié)果的異常，這已為模型解釋（第4節(jié)）所排除.此外，訓(xùn)練集與驗(yàn)證集截取時(shí)段的嚴(yán)格分隔也確保了二者的有效分離.然而，調(diào)參與訓(xùn)練階段對樣本的共享作為潛在的數(shù)據(jù)泄露機(jī)制，有待于樣本的擴(kuò)充與細(xì)化.

過擬合也會(huì)導(dǎo)致模型虛高的準(zhǔn)確率，具體表現(xiàn)為測試集的準(zhǔn)確率相對訓(xùn)練集顯著偏低.訓(xùn)練樣本欠缺和模型過高的復(fù)雜程度是過擬合的重要誘因.研究采用L2正則化、最大池化，以及早停降低過擬合風(fēng)險(xiǎn).本節(jié)將從模型選擇、訓(xùn)練樣本數(shù)量、L2正則化等三方面進(jìn)行討論.

受隨機(jī)性影響，模型選取會(huì)以多次訓(xùn)練的最高準(zhǔn)確率結(jié)果為準(zhǔn).由于CNN訓(xùn)練往往收斂于局部最優(yōu)點(diǎn)而全局最優(yōu)點(diǎn)常為過擬合點(diǎn)，有必要考察最高準(zhǔn)確率模型（更接近全局最優(yōu)點(diǎn)）的過擬合情況.圖15a對比了最高準(zhǔn)確率模型相對其它模型（5層第4輪訓(xùn)練）的震相得分，可見前者的P，Sg得分反而偏低，暗示模型收斂點(diǎn)對過擬合的影響十分有限.

訓(xùn)練樣本數(shù)量也是影響模型過擬合的重要因素.為方便對比，選取準(zhǔn)確率相近的4層模型（A=0.962±0.03）.圖15b展示了不同訓(xùn)練樣本規(guī)模下的震相得分，未見震相得分的系統(tǒng)變化.

L2正則化通過平滑CNN參數(shù)來控制模型的復(fù)雜程度，其最優(yōu)系數(shù)隨模型深度增加（6.2節(jié)）.連續(xù)波形掃描所用5層模型在訓(xùn)練過程中用到的正則系數(shù)（30）顯著小于該深度的最優(yōu)取值（80）.圖15c對比了不同正則化系數(shù)訓(xùn)練下5層模型的震相得分（模型準(zhǔn)確率均為0.971），可見L2正則化對震相得分的抑制作用.

圖15 不同模型的震相得分對比（a） 模型選擇；（b） 樣本規(guī)模；（c） 正則化系數(shù).黑色實(shí)線對應(yīng)y ＝ x ，震相得分經(jīng)單調(diào)變換 x → ? lg（1－x）以便于展示Fig.15 Phase scores for the selected models（a） Model selection；（b） Training sample size；（c） L2 regularization factor.The black solid line corresponds to y＝ x ，with phase scores transformed monotonously by x → ? lg（1－x） for better view

綜合上述，正則化系數(shù)的不當(dāng)選取是導(dǎo)致震相得分虛高的重要原因，不排除數(shù)據(jù)泄露的可能影響.盡管研究涉及的樣本規(guī)模不足以造成過擬合程度的系統(tǒng)變化，高質(zhì)量訓(xùn)練樣本的積累依然是CNN訓(xùn)練質(zhì)量的前提和保證.

7 結(jié)論

本文探討了卷積層深度、正則化、殘差塊添加和數(shù)據(jù)清洗等因素對CNN震相（P，S，Pg，Sg）拾取能力的影響.結(jié)果表明，適中的卷積層深往往具有更優(yōu)的性能.深度為4，5的CNN 網(wǎng)絡(luò)分別具有（8，64），（4，64）的最末卷積層（圖2a）.在提取足量特征（64）的同時(shí)，上述層位依舊具有一定的局部分辨能力，從而能有效地處理震相序列沿時(shí)間軸的分布.與此同時(shí)，合適的正則化系數(shù)及數(shù)據(jù)清洗也能有效提高震相識別的準(zhǔn)確率.相比而言，殘差塊的有限作用可能與地震數(shù)據(jù)（相對二維圖像識別）較低的復(fù)雜度（于子葉等，2020）有關(guān).

為了客觀評估模型的震相識別機(jī)制，采用CMV和平滑GradCAM++解釋所得CNN.其中CMV反演波形基本復(fù)現(xiàn)了各震相的主要特征，包括不同震相垂直與水平分量的相對大小.而平滑GradCAM++則揭示了前方平靜段（P與Pg）或前驅(qū)能量（Sg）對震相識別的重要作用.與此同時(shí)，模型解釋也揭示了模型及訓(xùn)練的問題.對于CMV，S震相波形噪聲過大，Sg波形初動(dòng)偏離過大；對于平滑GradCAM++，S震相的波形區(qū)段對得分幾乎不作任何貢獻(xiàn).而短窗分量的過低參與度導(dǎo)致長窗垂直分量主導(dǎo)CNN預(yù)測，不利于多分量地震圖的綜合利用.

最后，將CNN模型用于連續(xù)波形的掃描.掃描按長、短步長依次識別、拾取震相.結(jié)果表明模型具有初步的遠(yuǎn)-近震相拾取能力，但靈敏度和抗噪性亟需提升.在討論章節(jié)，我們初步分析了訓(xùn)練樣本規(guī)模、正則化參數(shù)、數(shù)據(jù)泄露和過擬合等因素對所得結(jié)果的影響，為后續(xù)深入研究提供了一定參考.

切實(shí)提高訓(xùn)練樣本的質(zhì)量與數(shù)量，改善S震相訓(xùn)練水平并平衡各分量的模型參與度，是未來深入研究的重要方向.為此，后續(xù)研究可從以下幾個(gè)方面著手：

1） 適當(dāng)延長震相S的選取時(shí)段，或從理論地震圖生成訓(xùn)練樣本，從源頭緩解震相的不足.

2） 在模型架構(gòu)上，輸入的長短窗分量流入各自的卷積?池化序列，以壓制長窗垂直分量的優(yōu)勢地位（附件1中圖1a）.

3） 采用軟標(biāo)簽編碼訓(xùn)練樣本，適當(dāng)平移震相波形，在豐富樣本的同時(shí)優(yōu)化得分函數(shù)對走時(shí)的敏感（附件1中圖1b）.考慮到短窗長度，震相得分的非零區(qū)間半徑取5秒.對于Pg與Sg的中間時(shí)段，該值參考Pg-Sg的到時(shí)差值（中點(diǎn)處Pg，Sg標(biāo)簽均為0.5）.

4） 隨機(jī)抽樣生成的訓(xùn)練樣本集合單獨(dú)參與CNN的訓(xùn)練，得到各自對應(yīng)的子模型，最終結(jié)果取決于每個(gè)子模型得分的“投票”結(jié)果.這種集成學(xué)習(xí)策略有望利用隨機(jī)性消除局部極值點(diǎn)的影響，壓制過擬合造成的“極端”預(yù)測結(jié)果.

本文模型訓(xùn)練采用的操作系統(tǒng)為Fedora 30，中央處理器為Intel Core i7-6700CPU@3.40 GHz×8，硬盤為 ST1000LM044 HN-M101SAD （1 TB，7200 轉(zhuǎn)），1 個(gè) epoch （以清洗后樣本集，深度為5的CNN為例）耗時(shí)約10分鐘.模型的搭建與訓(xùn)練在python3.7下的tensorflow2.2環(huán)境完成，震相數(shù)據(jù)處理采用obspy包，圖件繪制采用matplotlib和cartopy包.中國地震局地球物理研究所李麗研究員為本文提出了寶貴建議.北京國家觀象臺(tái)朱戰(zhàn)斌高級工程師與周江林高級工程師為本文提供了數(shù)據(jù)支持，審稿專家為本文提出的寶貴修改意見，作者在此一并表感謝.