蔡振宇 蓋增喜

?

人工智能在拾取地震 P 波初至中的應(yīng)用——以汶川地震余震序列為例

蔡振宇 蓋增喜?

北京大學(xué)地球與空間科學(xué)學(xué)院, 北京 100871

為了準(zhǔn)確而迅速地拾取大量地震事件的 P 波初至, 將深度學(xué)習(xí)方法引入微地震 P 波初至到時拾取研究中, 對卷積神經(jīng)網(wǎng)絡(luò)的結(jié)構(gòu)進(jìn)行改造, 以便適應(yīng)地震波形數(shù)據(jù)的特點和 P 波初至拾取的要求。該算法只需要輸入 10s 窗口的三分量地震波形數(shù)據(jù), 就可以自動地判定 P 波初至?xí)r刻, 無需掃描連續(xù)波形, 運算時間遠(yuǎn)遠(yuǎn)小于長短窗、模板匹配等傳統(tǒng)方法。使用該算法訓(xùn)練汶川地震主震后 2008 年 7—8 月 7467 條人工拾取的余震 P 波初至到時, 將得到的模型對測試集中 1867 條數(shù)據(jù)的計算結(jié)果與人工拾取結(jié)果對比,誤差小于 0.5s 者占比達(dá)到 98.9%。在低信噪比條件下, 該方法仍能保持較好的拾取能力。

人工智能; 機器學(xué)習(xí); 深度學(xué)習(xí); 小波變換; 初至拾取

汶川地震主震發(fā)生后, 不斷地有余震發(fā)生。為了更好地開展抗震救災(zāi)工作, 需要更快速和精準(zhǔn)的手段來識別余震的初至?xí)r刻并構(gòu)建余震地震目錄。如果能夠自動地識別這些微弱的地震信號和遺漏的早期余震, 找到余震震相的初至?xí)r刻, 將對完善地震目錄、地震預(yù)警、全面分析地震活動性以及監(jiān)測地震破裂區(qū)域的震后形變等具有重要意義。地震震相到時拾取是地震學(xué)研究的重要課題, 是地震預(yù)警、地震定位及地球內(nèi)部結(jié)構(gòu)等研究的基礎(chǔ)。地震預(yù)警技術(shù)[1]是近二十年發(fā)展起來的一種減輕地震損失、減少人員傷亡、降低地震次生災(zāi)害的有效手段。通常情況下, 由于地震 P 波波速大于破壞性較大的 S 波和面波, 且地震波速度遠(yuǎn)遠(yuǎn)小于電磁波速度, 因此地震預(yù)警系統(tǒng)可在破壞性地震發(fā)生后, 基于近震源地區(qū)地震臺站監(jiān)測記錄得到的地震波形數(shù)據(jù), 快速地偵測地震, 在破壞性地震波到達(dá)特定目標(biāo)區(qū)前, 發(fā)布地震警報。

從技術(shù)層面上看, 地震震相的自動拾取是實現(xiàn)基于地震臺網(wǎng)地震預(yù)警的首要條件。在傳統(tǒng)地震學(xué)中, 一般是在地震圖上人工拾取震相。震相的自動識別開始于對海量波形數(shù)據(jù)的自動處理, 隨著實時地震學(xué)的發(fā)展, 對震相的自動精確拾取受到越來越多的重視。對地震預(yù)警而言, 目前主要考慮的震相為容易自動拾取的 P 波震相和 S 波震相。在各種地震震相自動拾取方法中, 廣泛應(yīng)用的主要有兩種, 一種是反映幅值瞬時變化的長短時窗比方法, 另一種是基于互相關(guān)技術(shù)的模板匹配濾波方法。這兩種方法或者依賴于閾值的設(shè)定, 或者依賴于模板的選取, 在一定程度上限制了應(yīng)用范圍。人工智能領(lǐng)域中近年來廣泛使用的神經(jīng)網(wǎng)絡(luò)模型可以根據(jù)輸入的訓(xùn)練數(shù)據(jù)訓(xùn)練參數(shù), 而不用人工設(shè)定參數(shù), 可能更適用于不同地區(qū)不同性質(zhì)的微地震震相到時拾取, 特別是在信噪比較低的情形下。

長短時窗比(STA/LTA)方法由 Stevenson[2]提出, Earle 等[3]完善了該方法, 并且將其應(yīng)用于地震研究中。該方法利用有效信號與噪聲的振幅差異, 通過兩個長度不同的滑動時窗來判斷初至?xí)r刻, 目前廣泛應(yīng)用于以微地震為代表的初至拾取研究中[4]。其中, 短時窗用于截取微弱的有用信號, 長時窗內(nèi)信號的平均值則反映背景噪聲水平, 包絡(luò)線上短時窗內(nèi)信號的平均值與長時窗內(nèi)信號平均值的比值用來表征信號振幅的變化。與其他震相拾取方法相比, 該方法具有算法簡單、計算速度快的特點, 對不同的震相均可拾取。但是, 該方法依賴于閾值的選取, 對于特定的情形, 需要人工設(shè)定閾值; 另一方面, 長短時窗比方法對信噪比的要求較高, 對信噪比較低的微地震拾取效果欠佳。

模板匹配濾波技術(shù)[5]是在滑動窗互相關(guān)檢測技術(shù)的基礎(chǔ)上發(fā)展起來的, 通過掃描連續(xù)波形來自動拾取震相, 是在信噪比較低的情況下提取微弱信號的一種有效方法。該方法將現(xiàn)有地震目錄中記錄的事件作為模板, 掃描連續(xù)波形, 從中識別出與模板事件波形有一定相似度的地震事件。與長短時窗比方法相比, 該方法具有在信噪比較低情況下仍適用的優(yōu)點。然而, 該方法也有其局限性, 一方面檢測出地震的數(shù)目依賴于所選取的模板數(shù)目, 另一方面, 只能識別重復(fù)的地震信號。與長短窗比方法類似, 該方法也需要設(shè)定互相關(guān)系數(shù)的閾值。近年來的新方法主要圍繞這兩種方法進(jìn)行改進(jìn)和優(yōu)化, 選取不同的特征函數(shù)來計算長短時窗比, 如 Zhang 等[6]2017年將分形引入長短時窗比方法中。但是, 這些改進(jìn)不能從本質(zhì)上解決閾值的人工設(shè)定以及低信噪比情形下效果差的問題。

以上兩類方法需要掃描連續(xù)波形, 存在運行速度慢的特點。如果能夠?qū)⑦B續(xù)數(shù)據(jù)切割成相互獨立的窗口作為輸入, 將大大減少運算時間。此外, 長短窗比方法中的特征函數(shù)和閾值參數(shù)均需要人工選取, 若能避開人工的局限, 直接從數(shù)據(jù)中挖掘出特征和參數(shù), 將大大簡化初至拾取的流程。

隨著大數(shù)據(jù)時代的來臨, 人工智能領(lǐng)域取得飛速的發(fā)展, 一方面數(shù)據(jù)量呈指數(shù)式增長, 另一方面計算機硬件計算能力不斷提升, 特別是隨著 GPU廣泛應(yīng)用于深度學(xué)習(xí)計算中, 過去應(yīng)用范圍受局限的深度學(xué)習(xí)方法開始應(yīng)用到各個領(lǐng)域中。深度學(xué)習(xí)方法具有結(jié)構(gòu)更復(fù)雜、可訓(xùn)練參數(shù)更多以及無需人工選取特征等優(yōu)點, 在眾多領(lǐng)域逐漸超越傳統(tǒng)機器學(xué)習(xí)方法, 得到更廣泛的應(yīng)用。

人工智能在地震學(xué)領(lǐng)域的應(yīng)用目前處于起步階段。過去幾十年來, 隨著寬頻帶地震儀的普及, 地震數(shù)據(jù)的質(zhì)和量都得到大幅度的提升, 為人工智能方法在地震學(xué)領(lǐng)域的應(yīng)用提供了數(shù)據(jù)量方面的支持。地震波和聲波同屬機械波, 地震波數(shù)據(jù)與聲波數(shù)據(jù)也有一些共同的特性, 這為將人工智能語音識別的一些方法拓展到對地震波數(shù)據(jù)的分析處理中提供了可能性。人工神經(jīng)網(wǎng)絡(luò)模型的雛形在 20 世紀(jì)40 年代就已提出, 20 世紀(jì) 90 年代就有人試圖將其應(yīng)用于震相的識別中[7], 但由于當(dāng)時數(shù)據(jù)量和計算能力等多方面的限制, 沒能取得好的效果。

鑒于上述背景, 本文另辟蹊徑, 繞開上述兩類方法的窠臼, 將人工智能領(lǐng)域的卷積神經(jīng)網(wǎng)絡(luò)結(jié)構(gòu)引入地震波 P 波初至?xí)r刻拾取中, 并根據(jù)地震波波形數(shù)據(jù)的特性和 P 波初至拾取問題的特殊性, 對傳統(tǒng)的卷積神經(jīng)網(wǎng)絡(luò)結(jié)構(gòu)進(jìn)行一定程度的調(diào)整, 充分利用深度學(xué)習(xí)方法在大數(shù)據(jù)量情形下的優(yōu)勢。我們將 P 波初至?xí)r刻拾取問題轉(zhuǎn)化為一個監(jiān)督學(xué)習(xí)中的回歸問題, 通過對大量已經(jīng)標(biāo)記好到時的連續(xù)波形數(shù)據(jù)進(jìn)行訓(xùn)練, 使機器學(xué)習(xí)得到一個可以重復(fù)使用的卷積神經(jīng)網(wǎng)絡(luò)模型。該模型一旦訓(xùn)練結(jié)束, 在實時拾取 P 波初至?xí)r, 只需要將連續(xù)波形數(shù)據(jù)切割成10s 的數(shù)據(jù)段作為輸入, 即可迅速得到拾取到時結(jié)果, 無須對原始波形進(jìn)行逐點掃描, 大大節(jié)省運行時間。此外, 在模型的訓(xùn)練過程中, 利用 GPU 加速, 使得訓(xùn)練時間也大大縮短。以本文具體使用的模型為例, 從輸入數(shù)據(jù)到訓(xùn)練出最終模型只需要 4~ 5 小時。與長短時窗比方法相比, 該模型的構(gòu)建既不需要人工選取振幅或能量等具有實際物理意義的特征函數(shù), 也不需要選取閾值。該方法也不需要選取模板, 而是直接輸入三分量連續(xù)波形數(shù)據(jù), 輸入數(shù)據(jù)被簡單地歸一化后, 即可通過訓(xùn)練好的模型得到拾取的到時時刻。在訓(xùn)練的過程中, 卷積神經(jīng)網(wǎng)絡(luò)通過一系列濾波器實現(xiàn)參數(shù)最優(yōu)化, 從而從波形數(shù)據(jù)中找到識別初至的特征, 這一切都是機器通過訓(xùn)練習(xí)得的, 不需要人工挑選特征或模板。

1 數(shù)據(jù)分析與預(yù)處理

根據(jù)中國地震臺網(wǎng)中心的統(tǒng)計, 汶川地震主震發(fā)生后, 2018 年 5 月共有 6.0 級以上余震 5 次, 5.0~ 5.9 級余震 25 次, 4.0~4.9 級余震 161 次; 6 月共有 6.0級以上余震 0 次, 5.0~5.9 級余震 3 次, 4.0~4.9 級余震31 次; 7 月共有 6.0 級以上余震 1 次, 5.0~5.9 級余震2 次, 4.0~4.9 級余震 15 次。當(dāng)震源周邊的地震臺站監(jiān)測到這些余震信號時, 通過提取 P 波和 S 波的初至?xí)r刻, 可以確定余震的震中位置, 為抗震救災(zāi)工作提供重要信息。然而, 在現(xiàn)實情況中, 由于受到背景噪聲水平、記錄儀器分布等多種因素的影響, 許多微弱的地震信號被淹沒在噪聲中, 無法直接觀測到。另一方面, 在大地震發(fā)生后的短時間內(nèi), 由于強震尾波的影響以及大量余震事件的發(fā)生, 波形相互重疊, 難以區(qū)分。

本文采用汶川地震后(2008 年 7 月 1—31 日)四川及鄰省 16 個地震臺站的三分量連續(xù)波形數(shù)據(jù)[8], 由國家數(shù)字測震臺網(wǎng)數(shù)據(jù)備份中心提供。

本文將三分量連續(xù)原始波形分割成時間窗口長度為 10s 的片段, 數(shù)據(jù)原始采樣率為 100Hz, 故數(shù)據(jù)格式為 1000×3 的二維矩陣。對每段數(shù)據(jù)進(jìn)行 1Hz 高通濾波, 濾掉長周期背景波動。然后, 對每段數(shù)據(jù)的每一道分量單獨進(jìn)行歸一化操作, 即將數(shù)據(jù)除以整段數(shù)據(jù)幅值絕對值的最大值, 使其最大值的絕對值為 1。本文采用人工標(biāo)記好 P 波到時的三分量數(shù)據(jù) 9334 條, 從中隨機抽取數(shù)據(jù)分為訓(xùn)練集、驗證集和測試集。其中, 5974 條(64%)作為訓(xùn)練集, 用于訓(xùn)練模型參數(shù); 1493 條(16%)作為驗證集, 用于判斷訓(xùn)練停止時刻; 1867 條(20%)作為測試集, 用于驗證模型樣本外的泛化能力, 各樣本集之間數(shù)據(jù)量的比例按照慣例[9], 測試集與訓(xùn)練集和驗證集之和的比例以及驗證集與訓(xùn)練集的比例均取 1:4。

為了準(zhǔn)確地拾取不同頻率特性的地震 P 波初至, 本文使用連續(xù)小波變換, 將每一分量的原始數(shù)據(jù)轉(zhuǎn)換為小波時頻譜。小波基采用中心頻率為 3 Hz, 頻率帶寬為 3Hz 的 Morlet 小波, 縮放因子取 10, 得到時間域為 0~10s, 采樣率為 0.01s, 頻率域為 1~ 10Hz, 頻率間隔為 1Hz 的 1000×10×3 三維矩陣。分別將時間序列數(shù)據(jù)和小波變換時頻譜數(shù)據(jù)作為輸入進(jìn)行訓(xùn)練, 并比較所得模型的識別能力。

2 模型構(gòu)建

本文以語音和圖像識別領(lǐng)域常用的卷積神經(jīng)網(wǎng)絡(luò)模型結(jié)構(gòu)為基礎(chǔ), 并針對地震波形數(shù)據(jù)的特點和到時拾取問題的特性, 對部分神經(jīng)網(wǎng)絡(luò)中的部分層級結(jié)構(gòu)做了一定的修改。當(dāng)以小波變換時頻譜作為輸入數(shù)據(jù)構(gòu)建模型時, 每一個頻率采樣點使用與以上神經(jīng)網(wǎng)絡(luò)結(jié)構(gòu)相同的模型, 最后通過一定方法計算不同頻率采樣點的回歸值, 得到最終的到時結(jié)果。

2.1 時間域波形數(shù)據(jù)卷積神經(jīng)網(wǎng)絡(luò)模型

McCulloch 等[10]首先將生物神經(jīng)網(wǎng)絡(luò)中的神經(jīng)元模型引入數(shù)學(xué)領(lǐng)域, 提出沿用至今的“M-P 神經(jīng)元模型”。在該模型中, 神經(jīng)元接收到來自個其他神經(jīng)元傳遞的輸入信號, 這些輸入信號通過具有權(quán)重的連接進(jìn)行傳遞, 神經(jīng)元接收到的總輸入值會與神經(jīng)元的閾值進(jìn)行比較, 然后通過激活函數(shù)處理來產(chǎn)生神經(jīng)元的輸出。將諸如邏輯斯諦函數(shù)這樣的非線性函數(shù)作為激活函數(shù), 當(dāng)神經(jīng)網(wǎng)絡(luò)的層數(shù)和結(jié)點數(shù)增加時, 理論上可以無限逼近任意非線性函數(shù)。因此, 神經(jīng)網(wǎng)絡(luò)模型具有極大的模型表達(dá)能力。

從理論上講, 參數(shù)越多的模型復(fù)雜度越高, 這意味著它能完成更復(fù)雜的學(xué)習(xí)任務(wù)。但是, 這種模型通常訓(xùn)練效率低, 易陷入過擬合, 樣本外測試效果差。隨著數(shù)據(jù)量的提升, 深度學(xué)習(xí)模型(層數(shù)很多的神經(jīng)網(wǎng)絡(luò)模型)的表達(dá)能力大大加強。在深度學(xué)習(xí)領(lǐng)域, 卷積神經(jīng)網(wǎng)絡(luò)自提出以來, 就在圖像識別和語音識別等領(lǐng)域有了廣泛的應(yīng)用[11]。鑒于地震波形數(shù)據(jù)與語音數(shù)據(jù)的相似性, 已經(jīng)有人將卷積神經(jīng)網(wǎng)絡(luò)應(yīng)用于地震信號和噪音的分類以及地震區(qū)域的定位工作中[12], 但將卷積神經(jīng)網(wǎng)絡(luò)應(yīng)用于到時拾取問題的研究中還無人問津。

一個典型的神經(jīng)網(wǎng)絡(luò)結(jié)構(gòu)主要由以下 5 個部分組成。

1)輸入層: 在到時拾取問題中, 輸入層為 1000 ×3×1 的三維矩陣, 1000 表示時間域的采樣點數(shù), 3 表示,,三個分量, 1 表示輸入層的深度。從輸入層開始, 卷積神經(jīng)網(wǎng)絡(luò)通過不同的神經(jīng)網(wǎng)絡(luò)結(jié)構(gòu)將上一層的三維矩陣轉(zhuǎn)化為下一層的三維矩陣, 直到最后的全連接層。

2)卷積層: 卷積神經(jīng)網(wǎng)絡(luò)中最重要的部分。與傳統(tǒng)的全連接層不同, 卷積層中每一個節(jié)點的輸入只是上一層神經(jīng)網(wǎng)絡(luò)的一小塊, 即卷積核。在圖像識別領(lǐng)域, 常用的卷積核大小為 3×3 或 5×5[13]。由于到時拾取問題的輸入層長寬比非常大, 本文中每個卷積層的卷積核大小均采用 2×1。卷積層試圖將神經(jīng)網(wǎng)絡(luò)中的每一小塊進(jìn)行更深入的分析, 從而得到抽象程度更高的特征。一般來說, 通過卷積層處理過的節(jié)點矩陣變得更深, 深度取決于過濾器深度。深度為的過濾器表示通過卷積層將原始數(shù)據(jù)的一張圖卷積成張圖。雖然每一層卷積核只利用臨近兩個采樣點數(shù)據(jù)之間的關(guān)系, 但通過疊加多個卷積層, 模型可以捕捉到更多臨近采樣點之間的關(guān)系來構(gòu)造特征。本文模型共疊加 5 個卷積層。

3)池化層: 池化層神經(jīng)網(wǎng)絡(luò)不改變?nèi)S矩陣的深度, 但可以縮小矩陣的大小?？梢哉J(rèn)為池化操作是將一張分辨率較高的圖片轉(zhuǎn)化為分辨率較低的圖片。通過池化層, 可以進(jìn)一步縮小最后全連接層中節(jié)點的個數(shù), 從而減少整個神經(jīng)網(wǎng)絡(luò)的參數(shù)。使用池化層可以加快計算速度, 并防止過擬合, 在一定程度上避免模型對訓(xùn)練集擬合效果較好, 但對測試集擬合誤差較大的問題。但是, 也有研究指出池化層對模型效果的影響不大[14]。通常情況下, 在將卷積神經(jīng)網(wǎng)絡(luò)應(yīng)用于分類問題(如圖像的分類)時, 池化層起到提取低頻抽象特征的作用。由于到時拾取問題的關(guān)鍵不在于提取抽象特征, 而是找到數(shù)據(jù)特征的突變時間點, 所以使用池化層會喪失高頻信息, 導(dǎo)致拾取結(jié)果的精度下降。因此, 本文模型結(jié)構(gòu)中不使用池化層, 直接將多層卷積層連接起來。

4)全連接層: 卷積神經(jīng)網(wǎng)絡(luò)一般由 1~2 個全連接層給出最后的分類結(jié)果。經(jīng)過多輪卷積層的處理后, 圖像中的信息已經(jīng)被抽象為信息含量更高的特征, 可以將卷積層看成自動圖像特征提取的過程。完成特征提取后, 使用全連接層完成分類任務(wù)。Perol 等[12]在處理地震和噪音的分類問題時, 只使用一層全連接層, 共兩個節(jié)點, 用于區(qū)分地震和噪音兩類結(jié)果。他們在處理地震定位分類問題時, 也只用一層全連接層, 共 10 個節(jié)點(文獻(xiàn)[12]中將分類結(jié)果分為 10 個地區(qū)塊)。到時拾取是回歸問題而不是分類問題, 要求輸出連續(xù)的浮點數(shù)值, 而不是表示類別的整數(shù), 因此需要更多的全連接層和節(jié)點來提高模型的表達(dá)能力。本文模型包含兩個全連接層, 第一個有 1024 個節(jié)點, 第二個只有 1 個節(jié)點, 用于輸出最后到時拾取結(jié)果。

5)柔性最大激活函數(shù)層: 主要用于分類問題, 計算得到當(dāng)前樣例屬于不同種類的概率。由于到時拾取并非分類問題, 故本文模型不使用該層。

綜上所述, 本文模型輸入層為 1000×3×1 的三維矩陣, 輸入層后連續(xù)疊加 5 個卷積層, 其過濾器深度分別為 6, 16, 16, 32 和 32, 提取的特征越來越深。除過濾器深度外, 5 個卷積層的其他參數(shù)均相同, 卷積核尺寸為 2×1, 卷積核的移動步長為 2, 激活函數(shù)為 ReLU 函數(shù)。ReLU 函數(shù)是人工神經(jīng)網(wǎng)絡(luò)中常用的激活函數(shù), 當(dāng)輸入信號為負(fù)值時, 輸出都是 0; 當(dāng)輸入信號為正值時, 輸出等于輸入。其表達(dá)式為

ReLU() = max(0,), (1)

其中,為神經(jīng)元的輸入。ReLU 函數(shù)被認(rèn)為有一定的仿生學(xué)原理, 并且在實踐中通常比其他常用激活函數(shù)(如邏輯斯蒂函數(shù))的效果更好, 因此廣泛使用于計算機視覺[15]人工智能領(lǐng)域。

在卷積層后設(shè)置一個隨機丟失層, 隨機丟失率為 0.5 (訓(xùn)練過程中每次更新參數(shù)時, 隨機地斷開50%輸入神經(jīng)元, 用于防止過擬合)。隨機丟失層后設(shè)置一個壓平層, 用來將輸入“壓平”, 即將多維輸入一維化, 完成卷積層到全連接層的過渡。壓平層后設(shè)置一個節(jié)點數(shù)為 1024, 激活函數(shù)為 ReLU 函數(shù)的全連接層。全連接層之后再設(shè)置一個隨機丟失率為 0.5 的隨機丟失層。最后一層是節(jié)點數(shù)為 1 的全連接層。

2.2 小波時頻譜卷積神經(jīng)網(wǎng)絡(luò)模型

考慮到不同震級的地震具有不同的頻率特性, 為了使模型對不同頻率特性的地震波形數(shù)據(jù)都具有良好的識別能力, 本文進(jìn)行基于小波變換的帶通濾波。小波變換是一種具有較好時頻局部化特性的方法, 具有可在時間域和頻率域同時進(jìn)行分析的特點, 在地震數(shù)據(jù)處理中廣泛應(yīng)用, 也應(yīng)用于震相拾取工作中[16]。小波變換將信號分解成一系列小波函數(shù)的疊加, 這些小波函數(shù)均由一個小波基函數(shù)經(jīng)過平移與尺度伸縮得到。利用這種不規(guī)則的小波函數(shù), 可以逼近非穩(wěn)態(tài)信號中尖銳變化的部分, 也可以逼近離散、不連續(xù)、具有局部特性的信號, 從而更真實地反映原信號在某一時間尺度上的變化。小波變換可以分為連續(xù)小波變換和離散小波變換。離散小波變換常用于降噪與數(shù)據(jù)壓縮, 連續(xù)小波變化更適用于信號特征的提取。作為時間序列間歇式波動特征提取的工具, 連續(xù)小波變化廣泛地應(yīng)用于地球物理學(xué)研究中[17]。

在連續(xù)小波變換中, 移動帶通濾波器的帶寬(頻率域窗口)隨著濾波器中心頻率的增加而擴(kuò)大。通常定義一個實信號()的連續(xù)小波變換為如下卷積:

其中, 窗口函數(shù)()為內(nèi)核小波, *號表示共軛。參數(shù)和分別為小波變換的尺度和平移量。對于每個尺度, 小波核按 1/的比例縮放, 按平移得到小波系數(shù)(a,b)。

為了有效地提取原始波形中的時頻信息, 小波基的選取十分重要。Farge[18]指出小波基選擇時需要考慮的因素, 例如正交與非正交、負(fù)值與實值以及小波基的寬度與圖形等。正交小波函數(shù)一般用于離散小波變換, 非正交小波函數(shù)既可用于離散小波變換, 也可用于連續(xù)小波變換[19]。在對時間序列進(jìn)行分析時, 通常希望得到平滑連續(xù)的小波振幅, 因此非正交小波函數(shù)較合適。此外, 要得到時間系列振幅和相位兩方面的信息, 就要選擇復(fù)值小波, 因為復(fù)值小波具有虛部, 可以對相位進(jìn)行很好的表達(dá)。Morlet 小波是高斯包絡(luò)下的單頻率正弦函數(shù), 沒有尺度函數(shù), 是非正交分解, 并且是由 Gaussian調(diào)節(jié)的指數(shù)復(fù)值小波[20]。Morlet 小波的圖形(圖 1)和地震波信號也具有一定的相似性, 因此本文選擇 Morlet 小波作為小波基函數(shù)。

小波時頻譜卷積神經(jīng)網(wǎng)絡(luò)模型的輸入數(shù)據(jù)為1000×10×3 的三維矩陣。使用中心頻率為 3Hz, 頻率帶寬為 3Hz, 縮放因子為 10 的Morlet 小波基, 經(jīng)過連續(xù)小波變換, 將每一分量的原始波形數(shù)據(jù)分別轉(zhuǎn)換為小波時頻譜, 其中頻率域范圍為 1~10Hz, 頻率間隔為 1Hz。對頻率域的每一個分量, 分別構(gòu)造一個如 2.1 節(jié)所述的卷積神經(jīng)網(wǎng)絡(luò)模型, 通過訓(xùn)練得到 10 個結(jié)構(gòu)相同、參數(shù)不同的卷積神經(jīng)網(wǎng)絡(luò)模型。將每個頻率上的 1000×3×1 三維矩陣作為輸入數(shù)據(jù), 通過模型計算得到每個頻率上的到時拾取值。對每個頻率, 計算到時拾取值與該頻率到時拾取值之差小于 1s 的頻率個數(shù)。以該頻率個數(shù)最大的頻率作為基準(zhǔn)頻率, 計算到時拾取值與該頻率到時拾取值之差小于 1s 的到時拾取值的平均值作為最終到時拾取值。

圖1 Morlet小波基(據(jù)文獻(xiàn)[20]修改)

2.3 卷積神經(jīng)網(wǎng)絡(luò)模型的訓(xùn)練過程

我們使用誤差逆?zhèn)鞑ニ惴ㄟM(jìn)行訓(xùn)練, 基于梯度下降策略, 按目標(biāo)的負(fù)梯度方向?qū)?shù)進(jìn)行調(diào)整。對每個訓(xùn)練樣例, 首先將輸入示例提供給輸入層神經(jīng)元; 然后逐層將信號前傳, 直到產(chǎn)生輸出層的結(jié)果; 再計算輸出層的誤差, 將誤差逆向傳播至隱層神經(jīng)元; 最后根據(jù)隱層神經(jīng)元的誤差, 對連接權(quán)和閾值進(jìn)行調(diào)整。該迭代過程循環(huán)進(jìn)行, 直到達(dá)到預(yù)設(shè)的停止條件。本文使用的停止條件為早停策略, 即將樣本內(nèi)數(shù)據(jù)分成訓(xùn)練集和驗證集, 訓(xùn)練集用來計算梯度, 更新連接權(quán)和閾值, 驗證集用來估計誤差。若最近的 1000 次訓(xùn)練中訓(xùn)練集誤差降低而驗證集誤差升高, 則停止訓(xùn)練, 同時返回具有最小驗證集誤差的連接權(quán)和閾值。訓(xùn)練過程中采用最小二乘誤差作為損失函數(shù)的計算損失, 使用學(xué)習(xí)速率為0.001的自適應(yīng)矩估計法來優(yōu)化參數(shù)。

3 模型結(jié)果

3.1 時間域波形數(shù)據(jù)卷積神經(jīng)網(wǎng)絡(luò)模型結(jié)果

基于時間域波形數(shù)據(jù)卷積神經(jīng)網(wǎng)絡(luò)模型, 對5974 條訓(xùn)練集數(shù)據(jù)進(jìn)行訓(xùn)練, 調(diào)整參數(shù), 將整個數(shù)據(jù)集中的 1867 條數(shù)據(jù)作為測試集, 當(dāng)測試集數(shù)據(jù)結(jié)果的最小二乘誤差在最近的 1000 次訓(xùn)練中都未減小時, 停止訓(xùn)練。為了便于評價模型的泛化能力, 我們將訓(xùn)練集和測試集結(jié)果的誤差分布做成直方圖(圖2)。

圖 2(a)顯示 7467 條訓(xùn)練集+驗證集樣本波形數(shù)據(jù)通過時間域波形數(shù)據(jù)卷積神經(jīng)網(wǎng)絡(luò)模型計算得到的到時拾取結(jié)果與人工拾取到時之間的誤差分布, 人工拾取到時的標(biāo)注精度為 0.1s, 65.8%的誤差小于 0.1s, 89.6%的誤差小于 0.2s, 98.4%的誤差小于0.5 s。

圖 2(b)顯示 1867 條測試集樣本波形數(shù)據(jù)通過時間域波形數(shù)據(jù)卷積神經(jīng)網(wǎng)絡(luò)模型計算得到的到時拾取結(jié)果與人工拾取到時之間的誤差分布, 42.8%的誤差小于 0.1s, 67.4%的誤差小于 0.2s, 90.6%的誤差小于 0.5s?？梢? 該模型在測試集的表現(xiàn)比訓(xùn)練集有較大的下降, 誤差大于 1.0s 的仍有 2.1%。因此, 我們對其中誤差小于 0.1s 和大于 1.0s 的數(shù)據(jù)分別畫出波形進(jìn)行分析, 結(jié)果如圖3所示。

從圖 3 看出, 對信噪比較高的數(shù)據(jù), 該模型有不錯的表現(xiàn)。圖 3(a)顯示使用時間域波形數(shù)據(jù)卷積神經(jīng)網(wǎng)絡(luò)模型對測試集中一條三分量波形數(shù)據(jù)拾取P 波初至到時的結(jié)果, 模型拾取到時比人工拾取到時滯后 0.03s。由于人工拾取到時的標(biāo)注精度為 0.1s, 故可以認(rèn)為對該條數(shù)據(jù)的拾取結(jié)果是準(zhǔn)確的。

對于信噪比較低的數(shù)據(jù), 該模型的識別能力較弱。圖 3(b)顯示使用時間域波形數(shù)據(jù)卷積神經(jīng)網(wǎng)絡(luò)模型對測試集中某條低信噪比數(shù)據(jù)拾取 P 波初至到時的結(jié)果, 模型拾取到時比人工拾取到時滯后 1.1 s, 可見模型的表現(xiàn)不理想。

3.2 小波時頻譜卷積神經(jīng)網(wǎng)絡(luò)模型結(jié)果

不同的微地震可能具有不同的頻率特性, 在不了解震源機制的情況下, 很難確定濾波頻率。如果能夠先對波形數(shù)據(jù)進(jìn)行小波變換, 將不同頻段的特征輸入提取到卷積神經(jīng)網(wǎng)絡(luò)中, 就可以讓機器去學(xué)習(xí)在不同頻段拾取到時而得到不同的模型, 從而將頻率因子納入模型中。

在新的模型中, 我們首先對數(shù)據(jù)做小波變換, 然后分別訓(xùn)練不同頻率的神經(jīng)網(wǎng)絡(luò)模型參數(shù)。對于輸入數(shù)據(jù), 每一個頻段由機器給出一個返回結(jié)果, 如圖 4 中黃線所示。對于每個頻段, 計算所有其他頻段的機器返回結(jié)果與該頻段的機器返回結(jié)果之差小于 1.0s 的頻段數(shù), 選取最大者對應(yīng)的頻段為基準(zhǔn)頻段, 最終輸出結(jié)果為與該基準(zhǔn)頻段結(jié)果之差小于1.0s 的所有頻段輸出結(jié)果的平均值。對于同樣的數(shù)據(jù), 該模型只比人工拾取結(jié)果滯后 0.03s, 在 0.1s的精度范圍內(nèi), 比不做小波變換的模型表現(xiàn)更好。

與 3.1 節(jié)的模型同樣, 我們統(tǒng)計了該模型在訓(xùn)練集和測試集誤差的分布情況。圖 5(a)顯示使用時間域波形數(shù)據(jù)卷積神經(jīng)網(wǎng)絡(luò)模型對 7467 條訓(xùn)練集+驗證集樣本波形數(shù)據(jù)計算得到的到時拾取結(jié)果與人工拾取到時之間的誤差分布, 人工拾取到時的標(biāo)注精度為 0.1s, 78.9%的誤差小于 0.1s, 97.1%的誤差小于 0.2s, 99.9%的誤差小于 0.5s。圖 5(b)顯示使用時間域波形數(shù)據(jù)卷積神經(jīng)網(wǎng)絡(luò)模型對 1867 條測試集樣本波形數(shù)據(jù)計算得到的到時拾取結(jié)果與人工拾取到時之間的誤差分布, 人工拾取到時的標(biāo)注精度為0.1s, 68.5%的誤差小于 0.1s, 90.0%的誤差小于 0.2s, 98.9%的誤差小于 0.5s。表 1 列出兩個模型的誤差統(tǒng)計分布, 可以看出, 小波變換明顯地提高了模型的識別精度。

圖2 時間域波形數(shù)據(jù)卷積神經(jīng)網(wǎng)絡(luò)模型拾取到時與人工拾取到時誤差分布

紅線表示人工拾取的到時, 藍(lán)線表示模型拾取的到時。下同

黃色豎線表示不同頻率模型拾取的到時, 綠色豎線表示模型最終給出的到時, 紅色豎線表示人工拾取的到時

圖5 小波時頻譜卷積神經(jīng)網(wǎng)絡(luò)模型拾取到時與人工拾取到時誤差分布

3.3 數(shù)據(jù)疊加高斯背景噪聲后模型的拾取能力

對比 3.1 節(jié)與 3.2 節(jié)的結(jié)果可知, 信噪比是模型拾取精度的重要影響因素。本文訓(xùn)練集和測試集的數(shù)據(jù)是隨機劃分的, 未對信噪比做任何要求。為了測試模型對信噪比的穩(wěn)定性, 選取其中信噪比較高的數(shù)據(jù), 通過疊加不同比例的白噪聲來對比信噪比對模型拾取結(jié)果的影響程度。圖 6(a)為使用小波時頻譜卷積神經(jīng)網(wǎng)絡(luò)模型對某條原始波形數(shù)據(jù)計算得到的到時拾取結(jié)果(圖 6 只展示分量波形圖), 模型拾取到時比人工拾取到時提前 0.06s。由于人工拾取到時的標(biāo)注精度為 0.1s, 故可以認(rèn)為該條數(shù)據(jù)拾取準(zhǔn)確。圖 6(b)~(h)依次為圖 6(a)中原始波形疊加其振幅最大值的 10%, 20%, 30%, 40%, 50%, 60%和 70%白噪聲后的疊加波形以及使用小波時頻譜卷積神經(jīng)網(wǎng)絡(luò)模型計算得到的到時拾取結(jié)果, 拾取誤差分別為?0.02, 0.09, 0.06, ?0.04, 0.03, 0.40 和0.59s。由于人工拾取到時的標(biāo)注精度為 0.1s, 故可以認(rèn)為疊加 50%振幅白噪聲時的拾取結(jié)果較準(zhǔn)確; 當(dāng)白噪聲大于 50%振幅后, 隨著信噪比降低, 拾取精度逐漸下降。在信噪比不是特別低的情況下, 該模型對于信噪比的穩(wěn)定性較好, 拾取結(jié)果沒有受到很大的影響。

表1 時間域與小波時頻譜卷積神經(jīng)網(wǎng)絡(luò)模型在測試集的誤差分布對比

3 結(jié)論

本文分別構(gòu)造兩個卷積神經(jīng)網(wǎng)絡(luò)模型來拾取 P波初至, 一個直接使用時間域波形數(shù)據(jù)輸出結(jié)果, 另一個先使用 Morlet 小波基對原始數(shù)據(jù)進(jìn)行小波變換, 然后對不同尺度下的小波變換波形分別構(gòu)造數(shù)個結(jié)構(gòu)相同的卷積神經(jīng)網(wǎng)絡(luò)模型, 綜合不同尺度下的模型輸出值而得到最終拾取結(jié)果。對比兩種模型的測試集結(jié)果誤差, 經(jīng)過小波變換的模型對 P 波初至的拾取準(zhǔn)確度更高, 誤差小于 0.1s 的占 68.5%, 小于 0.2s 的占 90.0%, 小于 0.5s的占 98.9%。

通過對可以準(zhǔn)確拾取的原始數(shù)據(jù)疊加不同比例的白噪聲, 分析該小波變換卷積神經(jīng)網(wǎng)絡(luò)模型在不同信噪比下的拾取準(zhǔn)確度, 結(jié)果表明, 當(dāng)白噪音與原始波形振幅之比小于 50%時, 拾取誤差仍然小于0.1s, 可見在較低信噪比的情況下, 該模型仍能保證準(zhǔn)確地拾取 P 波初至到時。當(dāng)白噪音與原始波形振幅之比達(dá)到 70%時, 該模型基本上失去拾取 P 波初至到時的能力。

除在低信噪比情況下仍能適用的優(yōu)點外, 本模型還具有運算速度快的特點。運算的主要耗時在于模型的訓(xùn)練, 一旦模型訓(xùn)練完成, 本文中 9334 條10s 長度的波形數(shù)據(jù)在不到 1s 的時間內(nèi)就能全部輸出拾取結(jié)果。中國地震局地球物理研究所與阿里云計算有限公司聯(lián)合主辦的阿里云天池“余震捕捉 AI大賽”要求震相拾取精度為 0.6s, 在此賽事中, 我們的模型對測試集中 1867 條數(shù)據(jù)的計算誤差小于 0.6s 的占 99.1%, 取得相當(dāng)高的樣本外準(zhǔn)確率。該模型不僅可以準(zhǔn)確地識別訓(xùn)練過的樣本, 對沒有經(jīng)過訓(xùn)練的新樣本的識別也有很高的準(zhǔn)確率。因此可以預(yù)期, 將該模型應(yīng)用于新的地震事件自動地拾取 P 波初至, 在大大節(jié)省運行時間的前提下, 也能取得高準(zhǔn)確率。

圖6 不同信噪比下數(shù)據(jù)小波時頻譜卷積神經(jīng)網(wǎng)絡(luò)模型與人工拾取到時對比

[1]Satriano C, Wu Y M, Zollo A, et al. Earthquake early warning: concepts, methods and physical grounds. Soil Dynamics and Earthquake Engineering, 2011, 31(2): 106?118

[2]Stevenson P R. Microearthquakes at Flathead Lake, Montana: a study using automatic earthquake process-ing. Bulletin of the Seismological Society of America, 1976, 66(1): 61?80

[3]Earle P S, Shearer P M. Characterization of global seismograms using an automatic-picking algorithm. Bulletin of the Seismological Society of America, 1994, 84(2): 366?376

[4]Liu H, Zhang J. STA/LTA algorithm analysis and improvement of microseismic signal automatic detec-tion. Progress in Geophysics, 2014, 29(4): 1708?1714

[5]Gibbons S J, Ringdal F. The detection of low mag-nitude seismic events using array-based wave-form correlation. Geophysical Journal International, 2006, 165(1): 149?166

[6]Zhang J, Tang Y, Li H. STA/LTA fractal dimension algorithm of detecting the P-wave arrival. Bulletin of the Seismological Society of America, 2017, 108 (1): 230?237

[7]張范民, 李清河. 利用人工神經(jīng)網(wǎng)絡(luò)理論對地震信號及地震震相進(jìn)行識別. 地震工程學(xué)報, 1998, 20 (4): 43?49

[8]鄭秀芬, 歐陽飚, 張東寧, 等. “國家數(shù)字測震臺網(wǎng)數(shù)據(jù)備份中心”技術(shù)系統(tǒng)建設(shè)及其對汶川大地震研究的數(shù)據(jù)支撐. 地球物理學(xué)報, 2009, 52(5): 1412? 1417

[9]Ripley B D. Pattern recognition and neural networks. Technometrics, 1999, 39(2): 233?234

[10]McCulloch W S, Pitts W. A logical calculus of the ideas immanent in nervous activity. Bulletin of Math-ematical Biology, 1990, 52(1/2): 99?115

[11]LeCun Y, Bengio Y. The handbook of brain theory and neural networks. Cambridge: MIT Press, 1995, 255? 258

[12]Perol T, Gharbi M, Denolle M. Convolutional neural network for earthquake detection and location. Scien-ce Advances, 2018, 4(2): e1700578

[13]LeCun Y, Bottou L, Bengio Y, et al. Gradient-based learning applied to document recognition. Procee-dings of the IEEE, 1998, 86(11): 2278?2324

[14]Springenberg J T, Dosovitskiy A, Brox T, et al. Striving for simplicity: the all convolutional net [EB/ OL]. (2014?12?21) [2018?04?13]. https://arxiv.org/ pdf/1412.6806.pdf

[15]Glorot X, Bordes A, Bengio Y. Deep sparse rectifier neural networks // Proceedings of the Fourteenth In-ternational Conference on Artificial Intelligence and Statistics. Crete, 2011: 315?323

[16]劉希強, 周蕙蘭. 用于三分向記錄震相識別的小波變換方法. 地震學(xué)報, 2000, 22(2): 125?131

[17]Chakraborty A, Okaya D. Frequency-time decomposi-tion of seismic data using wavelet-based methods. Geophisics, 1995, 60(6): 1906?1916

[18]Farge M. Wavelet transforms and their applications to turbulence. Annual Review of Fluid Mechanics, 1992, 24(1): 395?458

[19]Grinsted A, Moore J C, Jevrejeva S. Application of the cross wavelet transform and wavelet coherence to geophysical time series. Nonlinear Processes in Geophysics, 2004, 11(5/6): 561?566

[20]Torrence C, Compo G P. A practical guide to wavelet analysis. Bulletin of the American Meteorological Society, 1998, 79(1): 61?78

Using Artificial Intelligence to Pick P-Wave First-Arrival of the Microseisms: Taking the Aftershock Sequence of Wenchuan Earthquake as an Example

CAI Zhenyu, GE Zengxi?

School of Earth and Space Sciences, Peking University, Beijing 100871

In order to accurately and quickly pick up P-wave first-arrival of a large number of seismic events, deep learning method is introduced into the micro seismic P-wave first-arrival picking problem. The structure of convolution neural network is adjusted to apply to the characteristics of the seismic waveform data and first-arrival picking problem. The algorithm takes a 10s-window three-component seismic waveform data as input instead of scanning the continuous waveform. So the running time is far less than traditional methods such as STA/LTA and template matching. The algorithm is applied to aftershocks of 2008 Wenchuan earthquake in July and August, using 7467 manual picked first-arrival data as training dataset. Among the 1867 testing data, 98.9% of the P arrival times picked using this algorithm have an error less than 0.5 s compare to the results picked manually. This method can still maintain good pick-up capability under the condition of low signal-to-noise ratio.

artificial intelligence; machine learning; deep learning; wavelet transform; first-arrival picking

, E-mail: zge@pku.edu.cn

10.13209/j.0479-8023.2018.036

, E-mail: zge@pku.edu.cn

國家重點研發(fā)項目(2018YFC1504203)、國家自然科學(xué)基金(41774047)和中國地質(zhì)調(diào)查局地質(zhì)調(diào)查項目(DD20160082)資助

2018?05?18;

2018?06?24;

2019?04?03