李志強(qiáng)， 田有,2,3* ， 趙鵬飛,2， 劉財(cái),2， 李洪麗,2

1 吉林大學(xué)地球探測科學(xué)與技術(shù)學(xué)院， 長春 130026 2 吉林大學(xué)地球信息探測儀器教育部重點(diǎn)實(shí)驗(yàn)室， 長春 130026 3 長白山火山綜合地球物理教育部野外科學(xué)觀測研究站， 長春 130026

0 引言

海量地震數(shù)據(jù)的取得為人類認(rèn)識地球帶來了可能.隨著計(jì)算機(jī)硬件水平的快速提高，應(yīng)用計(jì)算機(jī)自動(dòng)處理地震數(shù)據(jù)，替代繁復(fù)機(jī)械的人工流程進(jìn)行諸如震相識別、走時(shí)初至拾取、接收函數(shù)挑選、余震序列識別等工作是十分必要的.傳統(tǒng)的自動(dòng)處理算法包括：利用短長時(shí)平均比值的STA/LTA方法(Allen, 1978; Withers et al., 1998)進(jìn)行震相識別，利用“V”字分形曲線首個(gè)突變點(diǎn)拾取震相的分形分維法(Sleeman and Van Eck, 1999；趙大鵬等, 2013)，利用不同波震相信噪比挑選接收函數(shù)等(Gao and Liu, 2014).然而這些方法在處理精度和效率上各有局限，仍需發(fā)展更加全面快捷的計(jì)算機(jī)自動(dòng)處理方案.

自Hinton等(2006)的開創(chuàng)性工作以來，深度學(xué)習(xí)蓬勃發(fā)展，并在地震大數(shù)據(jù)的自動(dòng)處理和油氣藏研究中嶄露頭角.付超等(2018)將卷積神經(jīng)網(wǎng)絡(luò)與支持向量機(jī)方法結(jié)合，提出了多波地震油氣儲(chǔ)層分布預(yù)測的深度學(xué)習(xí)方法；林年添等(2018)利用卷積神經(jīng)網(wǎng)絡(luò)預(yù)測地震油氣儲(chǔ)層；趙明等(2019a,b)利用U型卷積神經(jīng)網(wǎng)絡(luò)進(jìn)行Pg和Sg震相的識別和到時(shí)提取，常規(guī)的U網(wǎng)絡(luò)用來處理二維或者三維數(shù)據(jù)，將其降維設(shè)計(jì)即可用來處理一維地震波形；蔣一然和寧杰遠(yuǎn)(2019)利用支持向量機(jī)設(shè)計(jì)了地震信號探測器SSD和震相分離器SPS用來震相的分類和到時(shí)的拾??；李健等(2020)利用美國南加州地震臺(tái)網(wǎng)的200萬條震相和噪聲數(shù)據(jù)訓(xùn)練的卷積神經(jīng)網(wǎng)絡(luò)模型應(yīng)用到實(shí)際地震數(shù)據(jù)的震相拾取當(dāng)中取得了較好的效果；地震緊急預(yù)警中的應(yīng)用可見胡安冬和張海明(2020)的文章；奚先和黃江清(2020)應(yīng)用卷積神經(jīng)網(wǎng)絡(luò)進(jìn)行地震剖面中散射體的定位和成像.許多優(yōu)秀的模型和程序包也被分享到開源社區(qū)當(dāng)中供研究人員學(xué)習(xí)和使用，如用于識別震相的PhasePApy (Chen and Holland, 2016)；泛化能力較強(qiáng)的ConvNet (Ross et al., 2018)；斷層自動(dòng)識別軟件FaultSeg3D (Wu et al., 2018)；用以地震走時(shí)初至拾取的PickNet等(Wang et al., 2019; Geng and Wang, 2020).

深度學(xué)習(xí)在地震震相識別、走時(shí)拾取、斷層識別等領(lǐng)域的發(fā)展和應(yīng)用較快并取得了較好的成果，但是在接收函數(shù)研究中仍少有應(yīng)用.接收函數(shù)和地震波一樣都是一種時(shí)間序列，常規(guī)處理中同樣需要人工重復(fù)機(jī)械的挑選.將應(yīng)用于地震波中的深度學(xué)習(xí)方法遷移到接收函數(shù)挑選是可行且必要的.接收函數(shù)在特定的區(qū)域具有較為一致的形態(tài)，通常會(huì)根據(jù)P波、Ps轉(zhuǎn)換波以及PpSs+PsPs的震相是否清晰一致選擇保留或舍棄，這實(shí)質(zhì)上構(gòu)成了一種二分類問題，且特征數(shù)量較少，使得接收函數(shù)的挑選訓(xùn)練不需要巨量的數(shù)據(jù)，也不需要過于復(fù)雜的網(wǎng)絡(luò)即可達(dá)到令人滿意的準(zhǔn)確度.這一前提保證了即使是流動(dòng)臺(tái)站或其他觀測周期較短的臺(tái)站，也可以通過手動(dòng)挑選少量接收函數(shù)來構(gòu)建訓(xùn)練集，使用訓(xùn)練后的模型對剩余數(shù)據(jù)以及小于5.5級的小震數(shù)據(jù)進(jìn)行挑選，以實(shí)現(xiàn)節(jié)省人工、提高整體接收函數(shù)質(zhì)量的目的.而對于永久臺(tái)站，訓(xùn)練一次模型就可以供以后本臺(tái)所有數(shù)據(jù)的挑選使用，若能建立起區(qū)域臺(tái)網(wǎng)模型字典，將具有較高的應(yīng)用價(jià)值.

卷積神經(jīng)網(wǎng)絡(luò)(Convolutional Neural Network, CNN)是深度學(xué)習(xí)中具有代表性的一種，其人工神經(jīng)元可以使鄰近的神經(jīng)元之間互相產(chǎn)生影響.常規(guī)的卷積神經(jīng)網(wǎng)絡(luò)包括卷積層(Convolutional Neural Layers)、池化層(Max Pooling)和全連接層(Fully Connected Neural Layers，F(xiàn)CN).相比其他深度學(xué)習(xí)結(jié)構(gòu)，CNN在圖像識別和語音識別方面具有一定優(yōu)勢，且已經(jīng)有很成熟的開源框架，可以較容易的搭建出面向接收函數(shù)挑選問題的網(wǎng)絡(luò)結(jié)構(gòu).據(jù)此本文設(shè)計(jì)了一種CNN網(wǎng)絡(luò)來挑選接收函數(shù)，并采用實(shí)際觀測數(shù)據(jù)來驗(yàn)證方法的可行性.

1 數(shù)據(jù)和方法

1.1 接收函數(shù)提取

本文所使用的數(shù)據(jù)取自中國地震局的兩臺(tái)永久臺(tái)站，黑龍江省牡丹江地震臺(tái)(MDJ)和北京地震臺(tái)(BJT).兩個(gè)地震臺(tái)都具有較長的地震觀測周期，能夠提取大量接收函數(shù)應(yīng)用模型訓(xùn)練.兩個(gè)臺(tái)站使用的數(shù)據(jù)下載自美國地震學(xué)研究聯(lián)合會(huì)(Incorporated Research Institutions for Seismology, IRIS)網(wǎng)站(http:∥ds.iris.edu)，選取2000年到2019年的震級大于5.1級、震中距30°～ 90°的三分量事件波形數(shù)據(jù).根據(jù)P波理論到時(shí)截取波形P波前20 s，后60 s.之后進(jìn)行去均值、尖刺、傾斜操作，并對數(shù)據(jù)進(jìn)行0.05～10 Hz的帶通濾波，并由ZNE(垂、北、西)坐標(biāo)系旋轉(zhuǎn)至ZRT(垂、徑向、切向)坐標(biāo)系.采用時(shí)間域迭代反褶積方法(Ligorria and Ammon, 1999)計(jì)算接收函數(shù)，高斯因子設(shè)置為2.5.

對得到的MDJ臺(tái)15329條接收函數(shù)、BJT臺(tái)14524條接收函數(shù)進(jìn)行人工挑選和標(biāo)簽的制作.接收函數(shù)保存為sac文件格式，截取P波前5 s，后25 s共30 s的數(shù)據(jù)，采樣率20 Hz，保留的接收函數(shù)頭文件中t1參數(shù)設(shè)置為1，舍棄的t1參數(shù)設(shè)置為0.

1.2 神經(jīng)網(wǎng)絡(luò)搭建

本研究搭建的神經(jīng)網(wǎng)絡(luò)基于LeNet-5(LeCun et al., 1998, 2015)手寫數(shù)字識別模型.如圖1所示，整個(gè)網(wǎng)絡(luò)有7層，前4層為卷積層和池化層交替，輸入為一維序列，對其進(jìn)行卷積，使用大小為1×5的一維卷積核，滑動(dòng)步長1，深度16，選用ReLU函數(shù)作為激活函數(shù)，此時(shí)神經(jīng)元數(shù)量為600×16=9600個(gè).后接一層池化層，設(shè)定一個(gè)1×2的滑動(dòng)窗口，滑動(dòng)步長為2，選取其中的最大值輸出到下一層，神經(jīng)元數(shù)量減少到4800個(gè).重復(fù)進(jìn)行一次卷積、池化操作，神經(jīng)元數(shù)目進(jìn)一步減少到2400個(gè).后3層為全連接層，大小設(shè)置為256、60和2，采用Dropout(Srivastava et al., 2014)方法，按照0.5的概率舍棄部分神經(jīng)元以防止訓(xùn)練過程中的過擬合情況發(fā)生.

圖1 深度學(xué)習(xí)自動(dòng)挑選接收函數(shù)網(wǎng)絡(luò)結(jié)構(gòu)每個(gè)幾何圖形代表網(wǎng)絡(luò)中的一個(gè)節(jié)點(diǎn). 網(wǎng)絡(luò)下方列出了每一層的類型以及對應(yīng)的節(jié)點(diǎn)數(shù)、通道數(shù).Fig.1 Architecture of the deep learning network used for receiver functions auto-pickingEach shape corresponds to a node in the net. The characteristics of each layer, number of nodes and channels are also listed at the bottom.

對于接收函數(shù)挑選，只有兩種結(jié)果：保留、舍棄.這本質(zhì)上構(gòu)成了一種二分類問題，我們可以對輸出結(jié)果進(jìn)行One-hot編碼，保留編碼為[1,0]，舍棄編碼為[0,1].基于TensorFlow(ver:1.15.0, https: ∥www.tensorflow.org)軟件平臺(tái)和Python(ver:3.7.3)語言環(huán)境搭建深度學(xué)習(xí)模型，每次輸入100個(gè)帶標(biāo)簽的接收函數(shù)作為一個(gè)訓(xùn)練批次，采用加入L2正則化的交叉熵?fù)p失函數(shù)作為Loss函數(shù)，利用Adam優(yōu)化算法(Kingma and Ba, 2014)訓(xùn)練網(wǎng)絡(luò)參數(shù)，最大訓(xùn)練次數(shù)設(shè)置為30000次.訓(xùn)練采用一塊NVIDIA Quadro P5200顯卡，訓(xùn)練耗時(shí)140 s左右.

1.3 網(wǎng)絡(luò)訓(xùn)練

接收函數(shù)形態(tài)受地形影響較大，不同的地殼結(jié)構(gòu)獲得的接收函數(shù)具有一定的差異.例如高原地區(qū)的接收函數(shù)由于地殼過于復(fù)雜，各個(gè)震相會(huì)難以辨別(Shi et al., 2015)，而沉積盆地地區(qū)的接收函數(shù)直達(dá)P波、轉(zhuǎn)換波和多次反射波會(huì)被沉積層混響所覆蓋(朱洪翔等, 2018).因此采用單臺(tái)計(jì)算出的接收函數(shù)，每臺(tái)訓(xùn)練一個(gè)屬于自己的模型，將在一定程度上保證訓(xùn)練和測試樣本的相關(guān)性，從而降低學(xué)習(xí)難度，提高識別準(zhǔn)確率.

首先要構(gòu)建訓(xùn)練集和測試集，我們將兩個(gè)臺(tái)站的所有數(shù)據(jù)進(jìn)行人工挑選，并按照上文所述對t1參數(shù)進(jìn)行標(biāo)記.訓(xùn)練集使用2000年到2016年的數(shù)據(jù)，測試集使用2017年到2019年的數(shù)據(jù).MDJ訓(xùn)練集13531個(gè)接收函數(shù)中帶有保留標(biāo)簽的有1969個(gè)，占比14.5%；測試集1798個(gè)接收函數(shù)中帶有保留標(biāo)簽的有283個(gè)，占比15.7%；BJT訓(xùn)練集12730個(gè)接收函數(shù)中帶有保留標(biāo)簽的有3285個(gè)，占比25.8%；測試集1812個(gè)接收函數(shù)中帶有保留標(biāo)簽的有496個(gè)，占比27.3%.

兩個(gè)臺(tái)站的網(wǎng)絡(luò)結(jié)構(gòu)從8000次訓(xùn)練左右開始收斂穩(wěn)定，如圖2所示的Loss變化曲線.二分類問題常用的評價(jià)標(biāo)準(zhǔn)包括準(zhǔn)確率和召回率，準(zhǔn)確率表示正負(fù)樣本被正確分類的比例，召回率表示在實(shí)際為正的樣本中被預(yù)測為正樣本的比例，二者計(jì)算公式如公式(1)、(2)所示：

(1)

(2)

其中ACC表示準(zhǔn)確率，TPR表示召回率，NTP表示正類樣本被正確分類的數(shù)目，NTN表示負(fù)類樣本被正確分類的數(shù)目，NFP表示負(fù)類樣本被分為正類的數(shù)目，NFN表示正類樣本被分為負(fù)類的數(shù)目.

使用最終訓(xùn)練模型測試兩個(gè)臺(tái)站的測試集，MDJ臺(tái)的準(zhǔn)確率達(dá)到92.3%，召回率達(dá)到71.0%，BJT臺(tái)的準(zhǔn)確率達(dá)到93.1%，召回率達(dá)到84.8%.兩臺(tái)測試集的挑選情況以混淆矩陣熱力圖的形式在圖3給出.MDJ臺(tái)召回率較BJT臺(tái)低，經(jīng)分析是MDJ臺(tái)測試集正類樣本所占比例較低所致.兩臺(tái)站原始接收函數(shù)見圖4(a、e)，自動(dòng)挑選結(jié)果見圖4(b、f),人工挑選結(jié)果見圖4(c、g),可以看出自動(dòng)挑選結(jié)果與人工挑選結(jié)果具有較高的相似性.由于標(biāo)簽的制作具有比較大的主觀因素，因此存在將質(zhì)量稍差

圖2 Loss函數(shù)隨訓(xùn)練次數(shù)變化曲線(a) MDJ臺(tái); (b) BJT臺(tái).Fig.2 Variation curve of loss function during the training processing(a) MDJ station; (b) BJT station.

圖3 測試集混淆矩陣熱力圖(a) MDJ臺(tái); (b) BJT臺(tái).Fig.3 Heatmaps showing confusion matrix of testing datasets(a) MDJ station; (b) BJT station.

的接收函數(shù)標(biāo)記為保留，將質(zhì)量較好的接收函數(shù)標(biāo)記為舍棄，甚至?xí)心承┙邮蘸瘮?shù)在挑選時(shí)保留亦或舍棄都可以的情況.根據(jù)大多數(shù)正確分類接收函數(shù)規(guī)律訓(xùn)練出的模型，在測試時(shí)也會(huì)受到人為因素導(dǎo)致的錯(cuò)誤標(biāo)簽的影響，這些錯(cuò)誤分類的接收函數(shù)具有一定的“模糊性”，即使人工挑選也可留可棄，對后續(xù)處理影響較小.自動(dòng)挑選而未被人工挑選的接收函數(shù)見圖4(d、h)，可以看到這些接收函數(shù)分布沒有明顯規(guī)律，同時(shí)按人工挑選經(jīng)驗(yàn)，這些接收函數(shù)也是可以保留的.

圖4 人工挑選與深度學(xué)習(xí)自動(dòng)挑選的接收函數(shù)結(jié)果對比(a) MDJ臺(tái)原始未挑選的接收函數(shù)； (b) MDJ臺(tái)自動(dòng)挑選結(jié)果(MDJ_DL)； (c) MDJ臺(tái)手工挑選結(jié)果(MDJ_MA)； (d) MDJ臺(tái)自動(dòng)挑選而未被人工挑選的接收函數(shù)(MDJ_FP)； (e—h)與(a—d)相同，表示BJT臺(tái)對應(yīng)結(jié)果.Fig.4 Contrast of manual picking receiver functions and deep learning auto-picking receiver functions(a) All receiver functions in testing dataset of MDJ station; (b) Receiver functions of MDJ station after deep learning auto-picking (MDJ_DL); (c) Receiver functions of MDJ station after manual picking(MDJ_MA); (d) Receiver functions of MDJ station remained by deep learning but not manually remained; (e—h) are same as (a—d) but for BJT station.

2 地殼結(jié)構(gòu)與各向異性估計(jì)

為驗(yàn)證自動(dòng)挑選的接收函數(shù)在實(shí)際處理分析中的有效性，我們分別使用人工挑選的接收函數(shù)和自動(dòng)挑選的接收函數(shù)對兩個(gè)臺(tái)站下方的地殼厚度、波速比以及方位各向異性進(jìn)行計(jì)算，并對對比結(jié)果進(jìn)行分析.

2.1 地殼厚度與波速比

為估計(jì)兩個(gè)臺(tái)站下方的地殼厚度以及波速比，我們采用接收函數(shù)處理中常用的“H-κ疊加”(Zhu and Kanamori，2000)方法.Ps波、PpPs波以及PpSs+PsPs波相對于P波的時(shí)差(tPs、tPpPs、tPpSs+PsPs)可以表示為地殼厚度、地殼縱橫波速度比和射線參數(shù)的函數(shù).將一個(gè)臺(tái)站的所有接收函數(shù)按照公式(3)進(jìn)行疊加，即可估計(jì)臺(tái)站下方地殼厚度H和平均波速比κ.式(3)中，Ai為第i個(gè)接收函數(shù)對應(yīng)時(shí)差的幅值，ω1、ω2、ω3分別為三種震相疊加的權(quán)系數(shù)，MDJ臺(tái)權(quán)系數(shù)設(shè)為0.7、0.2和0.1，BJT設(shè)為0.8、0.1和0.1，N為該臺(tái)站參與疊加的接收函數(shù)的數(shù)目.

+ω3Ai(tPpSs+PsPs)}.

(3)

圖5(a、d)分別為MDJ臺(tái)和BJT臺(tái)測試集通過深度學(xué)習(xí)算法挑選的接收函數(shù)的H-κ疊加結(jié)果；圖5(b、e)分別為MDJ臺(tái)和BJT臺(tái)測試集通過人工挑選的接收函數(shù)H-κ疊加結(jié)果；圖5(c、f)分別為MDJ臺(tái)和BJT臺(tái)所有人工挑選的接收函數(shù)H-κ疊加結(jié)果.圖5(a、b、c)對比可以看出不同挑選方案得到的MDJ臺(tái)下地殼厚度和波速比分別為35.2 km/1.83、34.9 km/1.84、34.8 km/1.84；圖5(d、e、f)對比可以看出不同挑選方案得到的BJT臺(tái)下地殼厚度和波速比分別為38.4 km/1.78、37.9 km/1.79、37.8 km/1.80；

圖5 H-κ疊加獲得的兩臺(tái)站下方地殼厚度和平均波速比(a) MDJ臺(tái)自動(dòng)挑選數(shù)據(jù)H-κ疊加結(jié)果； (b) MDJ臺(tái)人工挑選數(shù)據(jù)H-κ疊加結(jié)果；(c)MDJ臺(tái)全部人工挑選數(shù)據(jù)H-κ疊加結(jié)果.(d—f)與(a—c)相同，表示BJT臺(tái)結(jié)果.Fig.5 Crustal thickness and average VP/VS beneath the two stations obtained by H-κ stacking(a) H-κ stacking results using deep learning auto-picking data of MDJ station; (b) H-κ stacking results of manual picking data of MDJ station; (c) H-κ stacking results of all manual picking data of MDJ station; (d—f) are same as (a—c) but for BJT station.

可以看出兩個(gè)臺(tái)站經(jīng)深度學(xué)習(xí)挑選的接收函數(shù)H-κ疊加結(jié)果和人工結(jié)果基本一致，說明神經(jīng)網(wǎng)絡(luò)用于接收函數(shù)的自動(dòng)挑選具有較好的適用性.

2.2 地殼各向異性

我們進(jìn)一步使用接收函數(shù)Ps震相估計(jì)臺(tái)站下方的地殼方位各向異性.在單層水平地殼假設(shè)下，接收函數(shù)Ps轉(zhuǎn)換波與直達(dá)P波的到時(shí)差可以由反方位角余弦表示(Liu and Niu, 2012; Zheng et al., 2018)，公式如下：

(4)

式(4)中，t為接收函數(shù)Ps轉(zhuǎn)換波與直達(dá)P波的到時(shí)差，t0為各向同性介質(zhì)假設(shè)下的到時(shí)差，Δt為各向異性介質(zhì)產(chǎn)生的到時(shí)差，dt是表征地殼各向異性強(qiáng)度的橫波分裂后快慢橫波到時(shí)差，Φ是表征地殼各向異性快軸方向的快橫波極化方向，Baz為反方位角.我們需要估計(jì)的即是t0、dt、Φ三個(gè)參數(shù).在本研究中，參數(shù)估計(jì)利用遺傳算法(Holland, 1992)進(jìn)行，采用二進(jìn)制編碼，每個(gè)參數(shù)由八位二進(jìn)制表示，總基因長度為24，t0、dt、Φ三個(gè)參數(shù)范圍分別設(shè)置為4～8 s、0～0.8 s、-90°～90°，所有接收函數(shù)Ps轉(zhuǎn)換波振幅值之和作為適應(yīng)度函數(shù)，用公式表示為：

(5)

式中F為適應(yīng)度函數(shù)，RFA(t)表示Ps轉(zhuǎn)換波振幅值，t由式(4)表示.種群設(shè)置為50，遺傳100代截止，交叉概率和變異概率分別設(shè)置為0.6和0.05.

接收函數(shù)預(yù)處理采用Zheng等(2018)提出的方法，為消除不同震中距對到時(shí)差的影響，Ps震相經(jīng)過了動(dòng)校正處理.接下來對接收函數(shù)進(jìn)行疊加平均，使用平均后的接收函數(shù)估計(jì)各個(gè)參數(shù).

圖6(a、e)分別為MDJ臺(tái)和BJT臺(tái)通過深度學(xué)習(xí)方法從測試集中挑選出的接收函數(shù)通過遺傳算法估計(jì)t0、dt、Φ三個(gè)參數(shù)的進(jìn)化圖；圖6(b、f)為圖6(a、e)各向異性估計(jì)結(jié)果所對應(yīng)的接收函數(shù)Ps波理論到時(shí)差曲線；圖6(c、g)分別為MDJ臺(tái)和BJT臺(tái)人工從測試集中挑選出的接收函數(shù)通過遺傳算法估計(jì)t0、dt、Φ三個(gè)參數(shù)的進(jìn)化圖；圖6(d、h)為圖6(c、g)各向異性估計(jì)結(jié)果所對應(yīng)的接收函數(shù)Ps波理論到時(shí)差曲線.將兩臺(tái)各向異性對比結(jié)果繪制在圖7所示的地形圖上，可以看出人工挑選的接收函數(shù)和深度學(xué)習(xí)挑選的接收函數(shù)估計(jì)出的各向異性結(jié)果比較一致.

圖6 兩臺(tái)站下方地殼方位各向異性結(jié)構(gòu)(a) MDJ臺(tái)自動(dòng)挑選出的接收函數(shù)，采用遺傳算法(GA)估計(jì)各向異性時(shí)的適應(yīng)度函數(shù)以及各向異性參數(shù)變化； (b) 與各向異性參數(shù)估計(jì)結(jié)果相對應(yīng)的Ps波理論到時(shí)差； (c—d)與(a—b)相同，表示MDJ臺(tái)人工挑選的接收函數(shù)計(jì)算結(jié)果； (e—h)與(a—d)相同，表示BJT臺(tái)自動(dòng)、人工挑選出的接收函數(shù)估計(jì)各向異性結(jié)果.Fig.6 Crustal azimuth-anisotropy structure beneath the two stations(a) The variation of fitness, anisotropic parameters in genetic algorithm (GA) using deep learning auto-picking receiver functions of MDJ station; (b) Theoretical Ps wave delay time corresponds to estimated anisotropic parameters; (c—d) are same as (a—b) but use all of manual picking receiver functions of MDJ station; (e—h) are same as (a—d) but for BJT station.

圖7 臺(tái)站所在位置及圖6中獲取的各向異性結(jié)果對比藍(lán)色實(shí)線為自動(dòng)挑選的接收函數(shù)估計(jì)的結(jié)果，紅色實(shí)線為人工挑選的接收函數(shù)估計(jì)的結(jié)果.Fig.7 Locations of the two stations and the contrast of anisotropic results obtained from Fig.6Blue lines represent the estimations from deep learning auto-picking receiver functions and the red lines are anisotropy results from all of manual selection receiver functions.

3 討論

3.1 訓(xùn)練集數(shù)目的影響

前文我們使用了2000—2016年16年的接收函數(shù)作為訓(xùn)練集，2017—2019兩年的接收函數(shù)作為測試集，測試集與訓(xùn)練集之比接近1∶ 7，最終模型自動(dòng)挑選接收函數(shù)可以達(dá)到較好的精度.然而，對于布設(shè)時(shí)間較短的流動(dòng)臺(tái)站或其他臨時(shí)臺(tái)站，無法提供大量的數(shù)據(jù)以供訓(xùn)練.為此我們探究了只應(yīng)用較小規(guī)模的訓(xùn)練集是否也能達(dá)到較好的效果.我們重新劃分訓(xùn)練集規(guī)模進(jìn)行訓(xùn)練，剩余數(shù)據(jù)用作測試集進(jìn)行試驗(yàn).

表1為MDJ臺(tái)不同訓(xùn)練集大小的訓(xùn)練結(jié)果，訓(xùn)練集規(guī)模從552條接收函數(shù)遞增，最終增加到13531條.最終訓(xùn)練模型應(yīng)用于測試集后準(zhǔn)確率均達(dá)到了90%以上，召回率70%以上.同時(shí)我們發(fā)現(xiàn)，對于召回率，測試集規(guī)模過大反而導(dǎo)致了召回率呈下降趨勢，僅用一年的數(shù)據(jù)訓(xùn)練的模型，準(zhǔn)確率保證的前提下，召回率反而得到了提升.經(jīng)過分析，我們發(fā)現(xiàn)導(dǎo)致這一現(xiàn)象的原因是對于較小規(guī)模的訓(xùn)練集，其正類樣本所占比例較高，表1 中Remained一列表示訓(xùn)練集中正類樣本所占的比例.如此一來，小規(guī)模訓(xùn)練集能夠獲得相對更好的挑選效果，這一前提為流動(dòng)臺(tái)站采取此類方法提供了可能.

表1 MDJ臺(tái)不同訓(xùn)練集規(guī)模情況下的測試集準(zhǔn)確率和召回率Table 1 The influence of testing data size on accuracy and recall of MDJ station

3.2 深度學(xué)習(xí)方法在<5.5級地震接收函數(shù)挑選的作用

在本研究中使用了近20年時(shí)間跨度的地震數(shù)據(jù)，使得我們可以通過統(tǒng)計(jì)學(xué)對接收函數(shù)挑選中所使用的地震震級情況進(jìn)行分析.如圖8所示，MDJ和BJT臺(tái)站接收到的震級范圍M5.1～5.5的事件占比分別為64.5%和64.7%，占比明顯高于震級≥M5.5的事件，但是從震級范圍M5.1～5.5的事件中人工挑選出的接收函數(shù)分別只占到該震級范圍內(nèi)總事件數(shù)的2.5%和7.07%，而從震級≥M5.5的事件中挑選出的接收函數(shù)占比達(dá)到了23.21%和16.29%，意味著在人工挑選過程中，是沒有必要花費(fèi)大量精力挑選

圖8 不同震級范圍計(jì)算的接收函數(shù)統(tǒng)計(jì)分析(a) MDJ臺(tái); (b) BJT臺(tái).角標(biāo)Y代表保留，N代表舍棄.Fig.8 Statistical analysis of receiver functions in different magnitude range(a) MDJ station; (b) BJT station. The subscript “Y” denotes being remained, “N” denotes being abandoned.

由于人工挑選數(shù)據(jù)的局限性舍棄掉的小震級接收函數(shù)數(shù)據(jù)，在使用自動(dòng)挑選時(shí)則可以得以保留.采用深度學(xué)習(xí)自動(dòng)挑選接收函數(shù)時(shí)，不需要擔(dān)心人工成本的問題.應(yīng)用訓(xùn)練好的模型可以快速從小于5.5級的地震事件中挑選出少量較高質(zhì)量的接收函數(shù)，從而提高數(shù)據(jù)的利用率.

3.3 多臺(tái)聯(lián)合訓(xùn)練

對于永久臺(tái)站，可以選取每個(gè)臺(tái)站幾年的數(shù)據(jù)進(jìn)行訓(xùn)練，將所有模型建立一個(gè)臺(tái)網(wǎng)字典供以后直接使用.但是對于大量的流動(dòng)臺(tái)站，其觀測周期較短，單臺(tái)接收函數(shù)訓(xùn)練結(jié)果泛化能力有較大的局限性，此時(shí)采用多臺(tái)數(shù)據(jù)聯(lián)合訓(xùn)練，建立一個(gè)公共的模型是有必要的.

受臺(tái)站數(shù)目限制，本次研究僅對上文使用的兩個(gè)臺(tái)站進(jìn)行聯(lián)合訓(xùn)練測試.我們將MDJ臺(tái)和BJT臺(tái)的訓(xùn)練集合并，測試集合并，同樣經(jīng)歷30000次訓(xùn)練，最終得到的模型應(yīng)用于測試集準(zhǔn)確性達(dá)到了92.5%，召回率達(dá)到81.3%，再次使用H-κ疊加方法對臺(tái)站下方的地殼厚度和各向異性進(jìn)行估計(jì)，如圖9所示，得到MDJ臺(tái)下地殼厚度為35.2 km，平均波速比為1.83；BJT臺(tái)下地殼厚度為37.9 km，平均波速比為1.79.兩臺(tái)最終計(jì)算結(jié)果與人工挑選的接收函數(shù)計(jì)算結(jié)果或單臺(tái)自動(dòng)挑選的接收函數(shù)計(jì)算結(jié)果均較為一致，表明建立區(qū)域多臺(tái)聯(lián)合接收函數(shù)挑選模型具有一定可行性.更多臺(tái)站準(zhǔn)確分析結(jié)果需要后續(xù)工作繼續(xù)跟進(jìn).

圖9 雙臺(tái)聯(lián)合訓(xùn)練自動(dòng)挑選的接收函數(shù)H-κ疊加結(jié)果(a) MDJ臺(tái); (b) BJT臺(tái).Fig.9 H-κ stacking results of two stations using joint training sets(a) MDJ station; (b) BJT station.

4 結(jié)論

本文發(fā)展了一種基于深度學(xué)習(xí)的接收函數(shù)自動(dòng)挑選方法，參照LeNet5模型構(gòu)建CNN神經(jīng)網(wǎng)絡(luò)，并使用中國地震局MDJ臺(tái)和BJT臺(tái)的事件波形提取接收函數(shù)，構(gòu)建訓(xùn)練集和測試集進(jìn)行驗(yàn)證.實(shí)際數(shù)據(jù)應(yīng)用結(jié)果表明，利用本文提出的深度學(xué)習(xí)自動(dòng)挑選接收函數(shù)方法具有較高的準(zhǔn)確率.對于自動(dòng)挑選的接收函數(shù)，我們使用H-κ疊加方法和Ps震相到時(shí)差計(jì)算了兩個(gè)臺(tái)站下方的地殼厚度、平均波速比以及地殼方位各向異性等地殼參數(shù)，得出以下結(jié)論：

(1)本文提出的深度學(xué)習(xí)方法挑選的接收函數(shù)和人工挑選的接收函數(shù)估計(jì)的地殼厚度、平均波速比和方位各向異性結(jié)果基本一致，證明了本文提出的自動(dòng)挑選方法的可行性.

(2)本研究提出的方法在構(gòu)建訓(xùn)練集過程中，對訓(xùn)練集數(shù)據(jù)數(shù)量有較低的依賴，試驗(yàn)中盡管只應(yīng)用幾百條接收函數(shù)訓(xùn)練，應(yīng)用最終模型挑選上萬條接收函數(shù)也可達(dá)到90%以上準(zhǔn)確率和80%以上的召回率.

(3)本研究提出的自動(dòng)挑選方法可以應(yīng)用于5.1～5.4級遠(yuǎn)震，有效地提高了觀測數(shù)據(jù)的利用率，為后續(xù)的分析、研究提供更好的約束.

(4)對于流動(dòng)臺(tái)站，同時(shí)使用多個(gè)臺(tái)站的數(shù)據(jù)構(gòu)建訓(xùn)練集以提高泛化能力，最終模型在測試集的表現(xiàn)也讓人滿意.

致謝感謝IRIS網(wǎng)站提供的地震數(shù)據(jù)(http:∥ds.iris.edu)，感謝吉林大學(xué)朱洪翔博士提出的寶貴意見，感謝審稿專家給出的修改意見.