基于UNet++的地震P波初至拾取研究

劉志文，王進強，b，王廣鑫

(北京科技大學 a.土木與資源學院，b.金屬礦山高效開采與安全教育部重點實驗室，北京 100083)

地震P波到時拾取是地震數據處理中關鍵的一環，拾取效率與精度直接影響震源定位的效率和精度。隨著數字地震臺網的大量建設，海量測震數據的產生，人工拾取P波不僅費時費力，還難以進行快速智能預警，自動化快速精確識別震相對于地震學的研究有著重要的意義[1]。傳統的STA/LTA(長短時窗均值比法)[2]方法雖然過程簡單，但精度不夠，對低信噪比的數據拾取效果不好，且閾值需要自己設定；利用AIC(赤池信息準則)[3]來尋找全局極小值作為P波初至，雖然效果要比STA/LTA好，但是計算復雜，對低信噪比數據的拾取效果不理想。1985年MAEDA提出M-AIC[4]，可跳過AR系數計算，直接用地震數據計算信號的AIC值。1999年SLEEMAN提出AR-AIC[5]模型用于震相的拾取，效果比STA/LTA好，但是依賴于時窗長度和AR系數的設定，也存在對低信噪比數據拾取效果不理想的情況。

語義分割是深度學習的一個研究領域，其將圖片中的每個像素點進行預測分類，并將每個類別的邊緣標注出來用不同的顏色對不同類別進行填充，從而達到邊緣檢測和分類的目的。將地震P波初至自動拾取問題看成一個二分類的問題(是否為P波初至)，需要判斷地震序列中的每個采樣點是否是P波初至。因為這一判斷和語義分割問題類似，因此可以借鑒語義分割模型對P波初至進行自動拾取研究。FCN(全卷積網絡)、U-Net(U型網絡)、UNet++、DeepLabV3都是很好的語義分割模型，在圖像處理方面取得了很好的效果。FCN是語義分割模型的先驅，它的結構全部由卷積層組成，所以叫全卷積網絡，其將來自深層、粗略層的語義信息與來自淺層、精細層的外觀信息相結合，以產生準確和詳細的分割[6]。U-Net形狀類似字母U的形狀，所以叫U型網絡，它通過下采樣抓住上下文的信息得到高級語義特征圖，再通過上采樣將圖片恢復到原分辨率進行輸出[7]。ZHU et al[8]基于U-Net網絡進行改進，提出了PhaseNet，相比于原始的U-Net網絡結構，PhaseNet每個水平層都少了一個卷積層。PhaseNet對P波初至拾取的準確率達到了93.9%，遠高于AR-AIC 55.8%的準確率，其對地震震相的高精度識別驗證了U型網絡可以很好地應用到地震P波初至的自動拾取上。趙明等[9]基于U-Net網絡提出了U型網絡對震相的拾取，將U-Net中加入了Dropout層，能有效防止模型過擬合，但拾取精度不佳。MOUSAVI[10]提出了EQT模型用于地震事件檢測和震相拾取，EQT模型包括多層卷積池化層、殘差網絡層、雙向LSTM層、LSTM層、Transformer層，最后將其分成三部分，用于檢測地震事件、P波到時、S波到時。P波到時預測的這一部分包括一層LSTM、局部注意力機制和多個上采樣、卷積操作，最后經過Sigmoid層輸出P波概率分布。雖然EQT在拾取精度上很高，但是該模型龐大臃腫、參數量太多，耗費計算資源也更大。

UNet++是基于U-Net改進的一種神經網絡結構，在語義分割的任務上的表現比U-Net好，主要是它整合了不同層次的特征，比U-Net增加了更多的特征拼接操作。理論上可以利用UNet++模型進行地震P波拾取，但相關研究未見報道。本文首先對UNet++網絡結構進行改造，然后進行地震波到時拾取研究，并與傳統方法進行對比，驗證本文方法的可行性與優越性。

1 方法和數據

1.1 基于UNet++的模型

UNet++是一個基于嵌套密集的跳躍連接的語義分割模型，結構背后的基本假設是：當來自編碼器網絡逐漸豐富的高分辨率特征圖在與來自解碼器網絡的相應語義豐富的特征圖融合之前，該模型可以更有效地捕捉對象的細粒度細節[11]。

為了使UNet++能應用于地震P波初至拾取，且達到很好的效果，本文基于原始網絡結構在3個方向進行了改進：網絡降維、改變單個Block中的操作和改變網絡的深度。

為了驗證單個Block中的操作對模型效果的影響，首先在保持原始網絡深度不變的情況下，改變單個Block中的操作，有增加或者減少單個Block中一組卷積和樣本歸一化操作的選擇，對比不同情況下在驗證集上漏拾數、正確拾取數、誤拾數的變化，實驗結果見表1.從表1可知增加一組操作，在誤拾取和正確拾取數上的結果優于不變和減少一組操作。雖然漏拾數比減少一組效果差，但從整體上來考慮，增加一組操作的效果明顯更好。結果也可以用理論來解釋，增加一組卷積和樣本歸一化操作，可以使神經網絡輸出更多的特征通道數，學到更多的地震P波特征，因此拾取P波的效果也更好。

確定每個Block中的操作后，再確定網絡結構的深度，分別做三組實驗，在3、4、5層網絡深度的情況下，網絡在驗證數據集上的表現，得到的結果如圖1所示。

表1 改變單個Block中的操作模型在驗證集上的表現Table 1 The validation set upon the change in the operation model in a single Block

圖1 不同網絡深度模型在驗證集上的表現Fig.1 Performance of different network depth models on the validation set

從圖1可以看出隨著網絡深度的增加，正確拾取數逐漸增加，漏拾數逐漸減少，誤拾數不是在網絡深度為5時表現最優。雖然誤拾數是在網絡深度為4時表現最優，但考慮到模型的任務主要是得到更多的正確拾取數，所以網絡深度確定為5層。

最終確定本文模型結構如圖2所示。圖中每個Block中的操作都一樣，Conv1、Conv2、Conv3都是一維卷積，卷積核大小為3，步長為1，填充為1，這樣做是為了使卷積后得到的結果和一維地震序列的長度保持一樣。樣本歸一化(Batch Normalizaion)操作有助于防止模型過擬合，其效果比添加Dropout層要好[12]。ReLU增加了神經網絡各層之間的非線性關系，且可以緩解梯度消失；池化操作(即下采樣)為最大池化，核大小為2，步長為2，經過池化后的數據長度會減半；上采樣操作的縮放因子為2，將下采樣后的樣本恢復到原長度，插值方式為線性插值。跳躍連接將每個Block得到的輸出在特征維度上進行向量的拼接。網絡結構中第一個Block接收地震數據的輸入，輸入的向量維度為(B,3,6 000)，B表示樣本批次大小，本文中取64，最后一個Block的輸出再經過一次一維卷積得到預測結果，輸出結果的向量維度為(B,2,6 000).

圖2 本文模型網絡結構圖Fig.2 Network structure diagram of the model in this paper

綜上，本文模型和UNet++不同之處有：1) 維度不同，UNet++是二維圖像的語義分割模型，地震數據是一維時間序列數據，不能直接用于地震P波的拾取。對其進行降維，卷積、池化、批次歸一化都降為一維操作；2) 網絡結構的深度不同：原始網絡深度為4層，本文模型為5層；3) 每個Block中的操作不同：原始網絡中每個Block包含Conv1(卷積1)、BN1(批次歸一化1)、Conv2(卷積2)、 BN2(批次歸一化2)、ReLU的操作，本文模型多了一組卷積和批次歸一化操作。

1.2 損失函數與優化器的選擇

在本文中將微震P波自動拾取看作一個二分類問題，損失函數選取交叉熵損失(CrossEntropyLoss)，公式如下：

(1)

式中：N為微震序列信號長度，yij是二值化編碼標簽，i=1，2表示有兩個類別分別為P波初至和非P波初至，其表達式為：

(2)

其中，z為最后一層輸出的張量，形狀大小為(m,2,N)，m為樣本批次大小。

深度學習中常用的優化器有SGD(隨機梯度下降)、Momentum、RMSprop[13]、Adam等，本文中的模型使用的是Adam，Adam優化器是Momentum和RMSprop的結合，能很快地找到正確方向并前進，且收斂速度快，可以讓模型更快地到達全局最優值。

1.3 數據預處理

訓練、驗證、測試數據均來源于STEAD[14](斯坦福地震數據集)，數據集的Github地址為https:∥github.com/smousavi05/STEAD，從該地址下載好需要的數據集。從源數據中篩選出三通道信噪比都低于某個值的數據，這里選擇低于20 dB的數據。選擇低信噪比的數據的原因是：高信噪比數據的P波自動拾取傳統算法就能做得很好，對低信噪比數據進行P波自動拾取研究更有意義。選出小于20 dB的數據共56 130條，每條數據長度都為6 000，采樣頻率為100 Hz.原始數據示例見圖3(a)，從圖中可以看出信號信噪比低，噪聲和有效信號混疊在一起，很難快速準確地確定P波初至點。

為了減少噪聲對模型的干擾，讓模型更好地學習P波波形特征，對原始波形數據進行降噪處理。目前，廣泛使用的微震信號降噪方法主要有以下幾類：1) 利用有效波和噪聲的頻率差異去噪；2) 利用有效波和噪聲的傳播方向 (視速度) 差異去噪；3) 利用有效波和噪聲的空間分布差異去噪；4) 利用微地震資料的區域統計差異去噪；5) 利用信號間的相干性區分噪聲；6) 基于神經網絡的微地震去噪方法[15]。本文中用到的地震數據震級大部分比2.5M小，為微震事件。針對微震信號的隨機非平穩性，降噪方法選取第三類中的小波閾值濾波法。小波閾值法是以小波變換理論為中心發展而來的，小波變換在信號的低頻部分具有良好的頻率分辨率，在高頻部分具有良好的時間分辨率[16-17]。利用這一優勢，DONOHO et al[18]提出一種基于小波變換的小波閾值濾波算法，TO et al[19]將其應用于微震信號降噪方面，并與維納濾波法進行對比，證明了這種方法可以用于微震信號降噪。閾值函數選擇軟閾值，硬閾值函數存在不連續點，易出現偽吉布斯點，選擇軟閾值函數使降噪后的信號更平滑；閾值選取為固定閾值，用固定閾值或啟發式閾值進行降噪與選極大極小閾值或無偏風險估計閾值進行降噪相比，前者的降噪效果比較徹底，而后者相對保守，前者降噪效果好。在微震信號降噪時，采用固定閾值或啟發式閾值降噪比較完全，更為有效[20]；小波函數選擇的是Sym8，Sym8小波基函數與地震信號的特征類似，有更好地提取微震信號的能力，降噪效果更好。通常小波分解的頻段范圍與采樣頻率有關，若N層分解，則各個頻段大小為Fs/2/2N[21].當小波函數和采樣周期選定之后，信號多分辨率分解的各層所占的頻帶是一定的，將各層頻帶與實際信號的有效頻帶相比較，就可以確定小波分解層數[22]。經計算分解層數為6，過大則會造成信號失真。最后對經小波閾值去噪后的波形進行歸一化處理，得到的數據如圖3(b)所示，從圖中可以明顯看出P波初至點，去噪效果好，并且對數據進行歸一化后，減少了不同振幅大小對模型特征學習的干擾。

圖3 數據降噪前后對比Fig.3 Comparison before and after data noise reduction

將數據集中實際P波到時的前后25個點都標注為1，其余標注為0.選前后25個點的原因是人工拾取的P波到時存在-0.25～0.25 s的誤差，本文數據采樣頻率為100 Hz，所以是前后25個點。圖4是將圖3(b)中Z分量數據進行標注處理前后的對比圖(為了讓P波初至和非P波初至標簽對比明顯，截取了前20 s的數據)，P波到時為7 s.

圖4 數據標注前后對比Fig.4 Before and after data labeling

將數據預處理好后，按照4∶1的比例劃分訓練集和驗證集，訓練數據有44 904條，驗證數據有11 226條。

1.4 模型訓練與驗證

用Pytorch框架搭建本文模型網絡結構，并在Kaggle(谷歌提供的虛擬機)上訓練和驗證，Kaggle配置為16G的運行內存，73G儲存，GPU為Tesla P100，最大功耗300 W.訓練和驗證樣本都按批次輸入到模型，每個批次的大小為64條原始數據，訓練和驗證的輪數為5，總共花了70 min.喂入一個批次進行訓練就是1個Time step(時間步)，總共有3 510個Time step，每輪有702個Time step，每50個Time step記錄一次損失函數的值。

每輪訓練完再喂入驗證數據集進行驗證，得到模型在驗證集上的表現，繪制其精確度和召回率隨訓練輪數變化的圖像，如圖5(b)所示。由于訓練數據龐大，模型在第一輪訓練后幾乎已經收斂，驗證集上數據的精確率和召回率不隨訓練輪數增加有明顯的增加，并且有時呈現下降的趨勢，但精確率基本維持在93.1%～94.3%之間，召回率在87%～92.2%之間。從圖5(a)中可以看出，繼續增加訓練輪數，訓練集上的損失值還會繼續下降，但下降幅度微小，在一個微小的區間內震蕩，趨于平穩，對精確率和召回率的提升可以忽略。由圖5(b)可知在第一輪訓練后，驗證集上的精確率就已經達到93.1%，第二輪繼續提升，之后三輪保持穩定。召回率在第一輪訓練后達到89.6%，第二輪有所下降，之后三輪緩慢攀升，有趨于穩定的態勢。綜上，選擇第5輪的模型作為最終模型，其在驗證集上精確率和召回率分別為94%、92.2%.

圖5 模型訓練與驗證過程Fig.5 Model training and verification process

2 測試和對比

測試數據集的信噪比也是低于20 dB，且經過和訓練、驗證數據一樣的處理，一共有150條。為了客觀評估本文網絡模型在測試集上的表現，我們選取了2種常用的震相拾取算法進行對比：STA/LTA和AR-AIC，對于3種方法的檢測結果，以人工到時拾取結果為標準，保證結果的可靠性。為評估3種方法的表現，使用平均值(μ)、方差(δ)、精確率(H)、召回率(R)作為評價指標，4個指標計算公式如下：

(3)

(4)

(5)

(6)

分別計算3種方法在測試集上得到的4個指標的值，結果見表2.其中STA/LTA和AR-AIC用的是Obspy(處理微震數據的Python庫)中的算法，代碼在本機中實現。在Obspy中，STA/LTA法拾取P波初至需傳入的參數為長短時窗寬度；AR-AIC法拾取P波初至需傳入的參數為低通頻率、高通頻率、長短時窗寬、AR系數和時窗方差。

表2 本文模型與STA/LTA、AR-AIC算法對比Table 2 Comparison of the model in this paper with STA/LTA and AR-AIC algorithms

從表2中可知， 本文模型在測試集上的精確率是98.00%，比STA/LTA高32%，比AR-AIC高22%，測試集數據信噪比大部分都低，如果數據信噪比高一些，精確率還會更高。平均值反映方法正確拾取P波時與人工拾取誤差絕對值的平均值，方差則反映方法正確拾取時，每次的拾取誤差絕對值和總體拾取誤差絕對值的平均值的分布是否集中，這兩個值都是越小越好。從表2中可知本文模型要優于AR-AIC，AR-AIC要優于STA/LTA.因為測試的樣本比較少，所以3種方法在召回率的表現上都為1，沒有漏拾的情況。

本文模型在150條測試數據上的表現：147條數據正確拾取，3條數據錯誤拾取(誤差分別為0.625 s、0.91 s、-1.06 s)，0條漏拾。為了觀察模型正確拾取時與人工拾取的誤差分布，繪制直方圖如圖6所示，從分布上來看大部分誤差集中在0～0.05 s，說明模型的拾取精度高，完全可以代替人工拾取。

圖7給出了本文模型在不同信噪比下與人工拾取的結果對比。圖7(a)中的數據雖然經過小波閾值去噪信噪比有所提高，但是起跳點依然不明顯，本文模型能夠無誤差地拾取P波到時；圖7(b)中本文模型拾取到時與人工拾取到時相差0.1 s，說明本文模型拾取延后；圖7(c)中本文模型拾取到時與人工拾到時相差-0.215 s，說明本文模型過早拾??；圖7(d)中本文模型拾取到時與人工拾到時相差0.625 s，模型拾取錯誤；圖7(b)、(c)中本文模型拾取結果和人工拾取結果誤差小，從圖中看兩條線重合在一起，將其部分放大才能看出差距，放大后見圖8(a)、(b).

圖6 模型拾取正確時與人工拾取到時的誤差分布Fig.6 Error distribution when the model is picked correctly and when it is picked manually

圖7 本文模型在不同信噪比下與人工拾取的結果對比Fig.7 Comparison of the model in this paper with the results of manual picking under different signal-to-noise ratios

從圖7(d)中可以觀察出對于信噪比太低的數據，P波起跳點被噪聲所淹沒，本文方法依然能夠進行到時拾取，但誤差較大。對于這種數據可以看看其它通道有沒有更好的數據或者改進去噪算法，讓P波起跳點更明顯，達到拾取P波初至的目的。

為了驗證模型的泛化能力，對輸入數據進行加噪處理，得到的預測結果和未加噪的數據進行對比，如圖9所示。圖9中實際P波到時是8 s，模型對未加噪數據的預測結果為8 s，對加噪10%數據的預測結果為8.01 s，可見模型對于加噪后的數據仍然能以很小的誤差拾取P波初至，表明模型泛化能力強。

PhaseNet和EQT是兩種先進的震相拾取深度學習模型，得到了廣泛的關注和初步的應用[23]。在P波初至拾取精度上PhaseNet達到了0.96，EQT為0.99[10]，本文模型為0.98.本文模型和PhaseNet在結構上存在很大的差異，本文模型增加了同層之間的殘差連接，從第二層開始就進行了上采樣，更多的殘差連接有助于避免訓練過程中的梯度消失問題，而從第二層開始上采樣能抓住不同層次之間的特征，讓網絡自己去學習不同深度的特征，因此在P波初至拾取精度上本文模型更好些。EQT模型龐大，有很深的網絡結構，需要訓練的參數量多，可能會造成模型過擬合，本文模型參數量與之相比要少很多，訓練所花費的資源比EQT模型少，但本文模型幾乎能達到和EQT同樣的P波初至拾取精度，因此本文模型具有良好的應用價值。

圖8 圖7(b)、(c)部分放大圖Fig.8 Enlarged view of part (b) and (c) of Fig. 7

圖9 模型泛化能力驗證Fig.9 Validation of model generalization ability

3 結論

1) 三種方法在測試集上的表現證明本文模型在P波初至拾取上明顯優于STA/LTA和AR-AIC方法，可以替代傳統P波自動拾取方法。

2) 模型訓練需要花較長的時間，依據訓練好的模型進行單個數據預測，耗時可以忽略不計(小于1 s)，能夠滿足實際需求。

3) 使用訓練好的模型對加噪后的數據進行測試，模型也能很好地拾取，驗證了本文模型泛化能力較好。

4) 面對日益增多的地震數據，每條數據的P波初至都由人工拾取顯得很不合理。本文將深度學習模型UNet++用于地震數據的P波初至拾取，取得了很好的效果，對推進深度學習在P波初至自動拾取上的應用有較強的現實意義。