基于深度學習的接收函數自動挑選方法

李志強， 田有,2,3* ， 趙鵬飛,2， 劉財,2， 李洪麗,2

1 吉林大學地球探測科學與技術學院， 長春 130026 2 吉林大學地球信息探測儀器教育部重點實驗室， 長春 130026 3 長白山火山綜合地球物理教育部野外科學觀測研究站， 長春 130026

0 引言

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

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

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

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

1 數據和方法

1.1 接收函數提取

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

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

1.2 神經網絡搭建

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

圖1 深度學習自動挑選接收函數網絡結構每個幾何圖形代表網絡中的一個節點. 網絡下方列出了每一層的類型以及對應的節點數、通道數.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.

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

1.3 網絡訓練

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

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

兩個臺站的網絡結構從8000次訓練左右開始收斂穩定，如圖2所示的Loss變化曲線.二分類問題常用的評價標準包括準確率和召回率，準確率表示正負樣本被正確分類的比例，召回率表示在實際為正的樣本中被預測為正樣本的比例，二者計算公式如公式(1)、(2)所示：

(1)

(2)

其中ACC表示準確率，TPR表示召回率，NTP表示正類樣本被正確分類的數目，NTN表示負類樣本被正確分類的數目，NFP表示負類樣本被分為正類的數目，NFN表示正類樣本被分為負類的數目.

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

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

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

的接收函數標記為保留，將質量較好的接收函數標記為舍棄，甚至會有某些接收函數在挑選時保留亦或舍棄都可以的情況.根據大多數正確分類接收函數規律訓練出的模型，在測試時也會受到人為因素導致的錯誤標簽的影響，這些錯誤分類的接收函數具有一定的“模糊性”，即使人工挑選也可留可棄，對后續處理影響較小.自動挑選而未被人工挑選的接收函數見圖4(d、h)，可以看到這些接收函數分布沒有明顯規律，同時按人工挑選經驗，這些接收函數也是可以保留的.

圖4 人工挑選與深度學習自動挑選的接收函數結果對比(a) MDJ臺原始未挑選的接收函數； (b) MDJ臺自動挑選結果(MDJ_DL)； (c) MDJ臺手工挑選結果(MDJ_MA)； (d) MDJ臺自動挑選而未被人工挑選的接收函數(MDJ_FP)； (e—h)與(a—d)相同，表示BJT臺對應結果.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 地殼結構與各向異性估計

為驗證自動挑選的接收函數在實際處理分析中的有效性，我們分別使用人工挑選的接收函數和自動挑選的接收函數對兩個臺站下方的地殼厚度、波速比以及方位各向異性進行計算，并對對比結果進行分析.

2.1 地殼厚度與波速比

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

+ω3Ai(tPpSs+PsPs)}.

(3)

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

圖5 H-κ疊加獲得的兩臺站下方地殼厚度和平均波速比(a) MDJ臺自動挑選數據H-κ疊加結果； (b) MDJ臺人工挑選數據H-κ疊加結果；(c)MDJ臺全部人工挑選數據H-κ疊加結果.(d—f)與(a—c)相同，表示BJT臺結果.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.

可以看出兩個臺站經深度學習挑選的接收函數H-κ疊加結果和人工結果基本一致，說明神經網絡用于接收函數的自動挑選具有較好的適用性.

2.2 地殼各向異性

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

(4)

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

(5)

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

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

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

圖6 兩臺站下方地殼方位各向異性結構(a) MDJ臺自動挑選出的接收函數，采用遺傳算法(GA)估計各向異性時的適應度函數以及各向異性參數變化； (b) 與各向異性參數估計結果相對應的Ps波理論到時差； (c—d)與(a—b)相同，表示MDJ臺人工挑選的接收函數計算結果； (e—h)與(a—d)相同，表示BJT臺自動、人工挑選出的接收函數估計各向異性結果.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 臺站所在位置及圖6中獲取的各向異性結果對比藍色實線為自動挑選的接收函數估計的結果，紅色實線為人工挑選的接收函數估計的結果.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 訓練集數目的影響

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

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

表1 MDJ臺不同訓練集規模情況下的測試集準確率和召回率Table 1 The influence of testing data size on accuracy and recall of MDJ station

3.2 深度學習方法在<5.5級地震接收函數挑選的作用

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

圖8 不同震級范圍計算的接收函數統計分析(a) MDJ臺; (b) BJT臺.角標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.

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

3.3 多臺聯合訓練

對于永久臺站，可以選取每個臺站幾年的數據進行訓練，將所有模型建立一個臺網字典供以后直接使用.但是對于大量的流動臺站，其觀測周期較短，單臺接收函數訓練結果泛化能力有較大的局限性，此時采用多臺數據聯合訓練，建立一個公共的模型是有必要的.

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

圖9 雙臺聯合訓練自動挑選的接收函數H-κ疊加結果(a) MDJ臺; (b) BJT臺.Fig.9 H-κ stacking results of two stations using joint training sets(a) MDJ station; (b) BJT station.

4 結論

本文發展了一種基于深度學習的接收函數自動挑選方法，參照LeNet5模型構建CNN神經網絡，并使用中國地震局MDJ臺和BJT臺的事件波形提取接收函數，構建訓練集和測試集進行驗證.實際數據應用結果表明，利用本文提出的深度學習自動挑選接收函數方法具有較高的準確率.對于自動挑選的接收函數，我們使用H-κ疊加方法和Ps震相到時差計算了兩個臺站下方的地殼厚度、平均波速比以及地殼方位各向異性等地殼參數，得出以下結論：

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

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

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

(4)對于流動臺站，同時使用多個臺站的數據構建訓練集以提高泛化能力，最終模型在測試集的表現也讓人滿意.

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