基于溶解氡數(shù)據(jù)和長短期記憶網(wǎng)絡(luò)的地震預(yù)報(bào)

劉海軍, 單維鋒, 耿貴珍

(1.防災(zāi)科技學(xué)院應(yīng)急管理學(xué)院，三河 065201；2.防災(zāi)科技學(xué)院經(jīng)濟(jì)管理學(xué)院，三河 065201)

幾十年來的研究成果表明，地震發(fā)生前孕震區(qū)及其周圍地區(qū)地下介質(zhì)應(yīng)力積累和介質(zhì)性質(zhì)變化的信息，有可能通過地下流體的動態(tài)變化特征表現(xiàn)出來。地下流體異常作為可靠地震前兆，在中國地震預(yù)報(bào)中起著舉足輕重的作用[1-3]。氡是一種放射性氣體，在巖石的孔隙和裂隙中存在，研究表明，氡反應(yīng)靈敏，當(dāng)受到外界壓力、振動等作用時(shí)，氡很容易從其賦存的介質(zhì)中逃逸出來，導(dǎo)致地下水中氡濃度出現(xiàn)不同程度的變化。由于氡對地下應(yīng)力變化敏感，是地下流體預(yù)報(bào)地震的靈敏組分之一，也是中國開展最早、最為普遍的水化學(xué)前兆觀測項(xiàng)目[4-5]。

氡異常檢測與分析是水氡觀測數(shù)據(jù)處理的主要工作[6-7]。溶解氣氡濃度觀測數(shù)據(jù)屬于典型的時(shí)間序列數(shù)據(jù)，具有短期的、中期的、長期的變化規(guī)律。影響地溶解氣氡濃度變化的因素很多，目前尚無水氡變化背景場模型[8-9]。利用溶解氣氡數(shù)據(jù)時(shí)間序列的特性，根據(jù)歷史觀測數(shù)據(jù)來研究溶解氣氡濃度變化規(guī)律，研究水氡變化模型，進(jìn)而對未來數(shù)據(jù)進(jìn)行預(yù)測，對于快速檢測異常數(shù)據(jù)、處理缺失數(shù)據(jù)等具有重要意義。目前地震前兆數(shù)據(jù)處理中使用的時(shí)間序列預(yù)測的主要方法有滑動平均(moving average，MA)模型、自回歸(auto-regressive，AR)模型、自回歸滑動平均(auto-regressive and moving average，ARMA)模型[1]及其變種等。這類方法均為線性模型，當(dāng)數(shù)據(jù)量較少時(shí)，線性模型效果良好。而數(shù)據(jù)量較大時(shí)，傳統(tǒng)的線性模型擬合能力有限。長短期記憶(long-short-term memory, LSTM)網(wǎng)絡(luò)[10-11]是一種深度循環(huán)神經(jīng)網(wǎng)絡(luò)。LSTM通過精心設(shè)計(jì)的含有輸入門、遺忘門、輸出門的神經(jīng)單元結(jié)構(gòu)，既能記住歷史數(shù)據(jù)中的信息，又能有效的解決數(shù)據(jù)長期依賴問題，在語音識別、視頻分析等領(lǐng)域取得了顯著的成就[12]，是目前最有效的時(shí)間序列分析技術(shù)，然而該技術(shù)還未有在溶解氣氡濃度預(yù)測方向的應(yīng)用?，F(xiàn)將LSTM技術(shù)引入溶解氣體氡濃度預(yù)測中，并將基于LSTM的溶解氣氡預(yù)測方法與傳統(tǒng)的AR方法、ARMA方法進(jìn)行對比實(shí)驗(yàn)。以期擴(kuò)寬溶解氣氡預(yù)測新思路。

1 模型與算法

1.1 圖問題數(shù)學(xué)描述

地震臺站某觀測點(diǎn)每日對溶解氣氡濃度觀測一次并記錄，某觀測點(diǎn)n個(gè)連續(xù)觀測日的觀測數(shù)據(jù)可以用含有n個(gè)數(shù)據(jù)的時(shí)間序列表示，組成集合X={x1，x2，…,xn}，其中xi為第i日的觀測值。本研究的問題可以定義為根據(jù)連續(xù)m日歷史觀測數(shù)據(jù)，預(yù)測其后面連續(xù)t日的觀測值，即根據(jù){x1，x2，…,xm}預(yù)測{xm+1，xm+2，…,xm+t}的值。具體而言，模型的輸入數(shù)據(jù)為連續(xù)m日的時(shí)間序列數(shù)據(jù)，輸出結(jié)果為從第m+1日開始連續(xù)t日的預(yù)測值。

1.2 數(shù)據(jù)預(yù)處理方法

溶解氣氡數(shù)據(jù)預(yù)處理分為兩部分：缺失數(shù)據(jù)處理和平穩(wěn)性處理。

(1)缺失數(shù)據(jù)處理

實(shí)驗(yàn)中的溶解氣氡數(shù)據(jù)來自中國三個(gè)不同地震臺站實(shí)際觀測數(shù)據(jù)。實(shí)際觀測中，由于儀器故障等因素導(dǎo)致不可避免地存在缺失數(shù)據(jù)。論文中選擇的時(shí)間段缺失數(shù)據(jù)較少，因此采用最近鄰方法進(jìn)行處理，即用與缺失數(shù)據(jù)最近的數(shù)據(jù)填充缺失數(shù)據(jù)。

(2)平穩(wěn)性處理

對溶解氣氡濃度進(jìn)行時(shí)間序列分析建模，首先要將數(shù)據(jù)變?yōu)槠椒€(wěn)數(shù)據(jù)。平穩(wěn)數(shù)據(jù)即數(shù)據(jù)的均值和方差在時(shí)間過程上都是常數(shù)，并且在任何兩時(shí)期的協(xié)方差值僅依賴于該兩時(shí)期的距離，不依賴于計(jì)算這個(gè)協(xié)方差的實(shí)際時(shí)間。實(shí)際觀測的溶解氣氡濃度數(shù)據(jù)往往為非平穩(wěn)過程，因此在建模之前首先要進(jìn)行差分處理，將其變?yōu)槠椒€(wěn)過程，進(jìn)而對其進(jìn)行預(yù)測，然后再進(jìn)行反差分運(yùn)算，得到最終的預(yù)測數(shù)據(jù)。差分計(jì)算公式為

Δxt=xt-xt-1

(1)

式(1)中：xt為時(shí)間序列t時(shí)刻的觀測數(shù)據(jù)。

圖1以姑咱地震臺溶解氣氡觀測數(shù)據(jù)為例，顯示了原始數(shù)據(jù)和差分之后數(shù)據(jù)對照圖。其中圖1(a)為某段時(shí)間溶解氣氡濃度原始觀測數(shù)據(jù)，可以看出，原始觀測數(shù)據(jù)具有明顯的周期性變化，隨著時(shí)間的變化，數(shù)據(jù)的均值呈周期性變化，為非平穩(wěn)數(shù)據(jù)。圖1(b)為進(jìn)行差分處理之后的數(shù)據(jù)，可以看出，經(jīng)過差分處理后，數(shù)據(jù)的均值和方差不受時(shí)間的影響，呈現(xiàn)平穩(wěn)狀態(tài)，可以進(jìn)行時(shí)間序列建模分析。

圖1 姑咱地震臺水化溶解氣氡數(shù)據(jù)差分前后對比

1.3 長短期記憶網(wǎng)絡(luò)模型

1.3.1 長短期記憶網(wǎng)絡(luò)基本結(jié)構(gòu)

溶解氣氡數(shù)據(jù)是一種時(shí)間序列數(shù)據(jù)，長短期記憶網(wǎng)絡(luò)是一種特殊的循環(huán)神經(jīng)網(wǎng)絡(luò)(recurrent neural network，RNN)，能夠從歷史時(shí)間序列中發(fā)現(xiàn)規(guī)律，對未來時(shí)間序列進(jìn)行預(yù)測。LSTM基本結(jié)構(gòu)如圖2所示，LSTM的展開結(jié)構(gòu)可以看出，最初時(shí)刻，LSTM處理單元A接受輸入數(shù)據(jù)x0，得到預(yù)測輸出h0，同時(shí)把h0傳遞到下一層，即h0與下一時(shí)刻的輸入數(shù)據(jù)x1一起作為輸入，共同預(yù)測h1，然后由h1和x2作為輸入一起預(yù)測h2，依此類推。

xt為t時(shí)刻的輸入數(shù)據(jù)，ht為t時(shí)刻的輸出數(shù)據(jù)，A為LSTM處理單元

從圖2可以看出，在LSTM神經(jīng)網(wǎng)絡(luò)中核心部分為LSTM處理單元A。LSTM處理單元A詳細(xì)結(jié)構(gòu)如圖3所示。

圖3 LSTM處理單元詳細(xì)結(jié)構(gòu)

從圖3中可以看出，LSTM處理單元包含三個(gè)“門”的結(jié)構(gòu)：遺忘門、輸入門和輸出門，這三個(gè)門結(jié)構(gòu)對應(yīng)的激活向量分別為ft、it和ot。單元狀態(tài)ct貫穿于三個(gè)“門”結(jié)構(gòu)之間。整體上看，除了h在隨時(shí)間流動，細(xì)胞狀態(tài)c也在隨時(shí)間流動，細(xì)胞狀態(tài)c就代表著長期記憶。這種精心設(shè)計(jì)的“門”結(jié)構(gòu)，使得LSTM具有數(shù)據(jù)記憶功能和數(shù)據(jù)遺忘功能。

LSTM處理單元中第一個(gè)“門”結(jié)構(gòu)是遺忘門。遺忘門決定LSTM單元從單元狀態(tài)c中忘記哪些信息，它檢查來自前一個(gè)LSTM處理單元的輸出向量ht-1和當(dāng)前時(shí)刻的輸出向量xt，并輸出0～1的數(shù)，0表示徹底忘記，1表示完全保留。遺忘門的激活向量ft計(jì)算公式如下：

ft=σg(Wfxt+Ufht-1+bf)

(2)

式(2)中：xt為t時(shí)刻輸入到LSTM處理單元中數(shù)據(jù)，W、U和b為參數(shù)，這些參數(shù)將在訓(xùn)練過程中從訓(xùn)練樣本中學(xué)習(xí)出來，σg為sigmoid函數(shù)，其方程如下：

(3)

it=σg(Wixt+Uiht-1+bi)

(4)

(5)

式(5)中：σc為tanh激活函數(shù)，其計(jì)算公式如下：

(6)

(7)

式(7)中，*運(yùn)算符為向量的Hadamard積運(yùn)算。

LSTM處理單元中第三個(gè)門為輸出門，該門根據(jù)狀態(tài)變量的值和輸出門激活函數(shù)的值來計(jì)算整個(gè)處理單元的輸出ht。輸出門的激活函數(shù)值ot如式(8)所示，LSTM處理單元的輸出ht計(jì)算計(jì)算公式如下：

ot=σg(Woxt+Uoht-1+bo)

(8)

ht=ot*σc(ct)

(9)

至此，t時(shí)刻的輸入數(shù)據(jù)xt和前一時(shí)刻的輸出數(shù)據(jù)ht-1，通過遺忘門、輸入門、輸出門結(jié)構(gòu)，得到了t時(shí)刻的輸出數(shù)據(jù)ht。

從整個(gè)LSTM處理單元數(shù)據(jù)流動可以看出，在LSTM結(jié)構(gòu)中，隨著時(shí)間流動的除了輸出數(shù)據(jù)ht之外，還有單元狀態(tài)ct。從圖3中可以看出，LSTM單元的狀態(tài)ct就像是一個(gè)傳送帶，將遺忘門的激活向量ft、輸出門的激活向量it和輸出門的激活向量ot串聯(lián)起來，傳遞向下一個(gè)處理單元，這是LSTM的核心。

1.3.2 長短期記憶網(wǎng)絡(luò)訓(xùn)練

前面介紹了LSTM的基本結(jié)構(gòu)，其中網(wǎng)絡(luò)的參數(shù)W、U和b需要在訓(xùn)練中確定。在訓(xùn)練階段，W、U和b隨機(jī)初始化，輸入數(shù)據(jù)xt，然后根據(jù)1.3.1節(jié)中的公式計(jì)算ft、ot、it、ct和ht,其中ht為t時(shí)刻的預(yù)測值，這個(gè)過程稱為前向計(jì)算過程。假設(shè)訓(xùn)練序列為{x1，x2，…,xn}，通過前向計(jì)算，可以得到預(yù)測值{h1，h2，…,hn}，而對于序列中的數(shù)據(jù)xt來說，其真實(shí)輸出yt，根據(jù)時(shí)間序列的性質(zhì)可知yt=xt+1。則模型在該訓(xùn)練數(shù)據(jù)上的誤差為

(10)

將網(wǎng)絡(luò)在訓(xùn)練集上的誤差L作為損失函數(shù)。進(jìn)而利用誤差反向傳播算法求解使得在訓(xùn)練集上損失函數(shù)L取得最小值的系數(shù)W、U和b為最終訓(xùn)練好的網(wǎng)絡(luò)參數(shù)。

2 實(shí)驗(yàn)與結(jié)果分析

2.1 數(shù)據(jù)和實(shí)驗(yàn)流程介紹

數(shù)據(jù)集1：姑咱水化1 461 d溶解氡連續(xù)觀測數(shù)據(jù)。姑咱泉位于四川省康定縣姑咱鎮(zhèn)，北緯30°07′，東經(jīng)102°09′，海拔1 410 m的康藏高原邊緣，正處于中國南北地震帶中部，是四川省三條主要地震帶即鮮水河斷裂帶、龍門山斷裂帶和安寧河斷裂帶的交匯部位。姑咱泉為多眼上升泉，泉水匯成一個(gè)深3.5 m，直徑約6 m的小池，水溫常年9 ℃左右。

數(shù)據(jù)集2：西昌川32井1 827 d溶解氡連續(xù)觀測數(shù)據(jù)。該井屬于西昌地震中心站地下流體觀測點(diǎn)，位于北緯27°85′，東經(jīng)102°26′，處于則木河斷裂上盤，海拔1 530 m，井深410 m，含水層巖性為礫巖，水質(zhì)類型屬重碳酸鈉型。

數(shù)據(jù)集3：雅安地震臺1 826 d溶解氡連續(xù)觀測數(shù)據(jù)。該臺站位于龍門山斷裂帶南段與滎經(jīng)—馬邊斷裂帶交匯部位，觀測點(diǎn)位于北緯30°17′，東經(jīng)103°03′，海拔700 m。為黏土層下覆泥巖深部層間水受層壓沿孔隙冒出的上升泉。

原始數(shù)據(jù)經(jīng)過缺失數(shù)據(jù)、平穩(wěn)性檢驗(yàn)、差分處理之后，對數(shù)據(jù)進(jìn)行建模并預(yù)測，然后進(jìn)行反差分處理，并統(tǒng)計(jì)誤差。實(shí)驗(yàn)流程如圖4所示。

圖4 基于LSTM的溶解氣氡預(yù)測實(shí)驗(yàn)流程

論文中所選擇的三個(gè)原始數(shù)據(jù)均缺失數(shù)據(jù)較少，采用最近鄰方法進(jìn)行缺失數(shù)據(jù)處理之后，如圖5所示。

圖5 論文中使用的原始觀測數(shù)據(jù)

從原始數(shù)據(jù)中可以看出，三組數(shù)據(jù)的均值隨著時(shí)間的變化有較明顯的周期性變動，為非平穩(wěn)時(shí)間序列，因此在預(yù)測前需要對原始數(shù)據(jù)進(jìn)行差分處理，處理之后的數(shù)據(jù)如圖6所示。

圖6 差分處理后的數(shù)據(jù)

從圖6中可以看出，經(jīng)過差分處理之后，數(shù)據(jù)的均值和方差不隨時(shí)間發(fā)生變化，變成平穩(wěn)的時(shí)間序列數(shù)據(jù)，可以進(jìn)行時(shí)間序列數(shù)據(jù)建模。

2.2 LSTM網(wǎng)絡(luò)結(jié)構(gòu)和參數(shù)

論文共選擇三個(gè)地震臺站的溶解氡觀測數(shù)據(jù)，分別為姑咱水化1 461 d觀測數(shù)據(jù)、西昌川32井1 827 d觀測數(shù)據(jù)和雅安地震臺1 826 d觀測數(shù)據(jù)。三個(gè)數(shù)據(jù)集分別經(jīng)過差分處理變成平穩(wěn)時(shí)間序列數(shù)據(jù)之后，即可采用LSTM網(wǎng)絡(luò)進(jìn)行建模。仿真實(shí)驗(yàn)采用MATLAB2018深度學(xué)習(xí)工具箱來搭建模型。本實(shí)驗(yàn)中的LSTM神經(jīng)網(wǎng)絡(luò)模型包括一個(gè)序列輸入層(sequence input layer)，用于接受輸入數(shù)據(jù)、一個(gè)LSTM層，用于對數(shù)據(jù)進(jìn)行建模、一個(gè)全連接層和一個(gè)回歸層，用于做預(yù)測。LSTM網(wǎng)絡(luò)建模最重要的參數(shù)為節(jié)點(diǎn)數(shù)量，節(jié)點(diǎn)數(shù)量采用實(shí)驗(yàn)的方法確定。LSTM網(wǎng)絡(luò)訓(xùn)練參數(shù)如表1所示。

訓(xùn)練過程中，將差分處理后的平穩(wěn)數(shù)據(jù)集分割為測試集和訓(xùn)練集兩部分，其中前90%的數(shù)據(jù)作為訓(xùn)練數(shù)據(jù)，其余10%作為測試數(shù)據(jù)。模型的評價(jià)指標(biāo)為均方根誤差(root-mean-square error，RMSE)，即選擇在訓(xùn)練集上均方根誤差最小值所對應(yīng)的模型。均方根誤差的計(jì)算公式為

表1 LSTM訓(xùn)練參數(shù)

(11)

LSTM神經(jīng)網(wǎng)絡(luò)的訓(xùn)練過程如圖7所示。圖7中顯示了隨著迭代的進(jìn)行，均方根誤差的變化曲線。網(wǎng)絡(luò)訓(xùn)練初始階段，均方根誤差的數(shù)值較高，隨著迭代的進(jìn)行，均方根誤差不斷減小，初始時(shí)，減小速度較快，隨著迭代的進(jìn)行，均方根誤差減少速度變緩，當(dāng)?shù)螖?shù)達(dá)到60以后，均方根誤差趨于平穩(wěn)，達(dá)到飽和，此時(shí)，神經(jīng)網(wǎng)絡(luò)訓(xùn)練完畢，可以用該網(wǎng)絡(luò)對訓(xùn)練數(shù)據(jù)進(jìn)行預(yù)測。

圖7 LSTM訓(xùn)練過程中RMSE的變化趨勢

神經(jīng)網(wǎng)絡(luò)訓(xùn)練好之后，可以應(yīng)用該網(wǎng)絡(luò)對未知數(shù)據(jù)進(jìn)行預(yù)測。在預(yù)測階段，將10%的測試數(shù)據(jù)集輸入訓(xùn)練好的LSTM網(wǎng)絡(luò)，得到預(yù)測數(shù)據(jù)，然后進(jìn)行反差分處理，反差分處理之后的數(shù)據(jù)為最終的預(yù)測值。

2.3 實(shí)驗(yàn)結(jié)果

應(yīng)用前面介紹的參數(shù)訓(xùn)練網(wǎng)絡(luò)，得到最優(yōu)模型并在測試集上進(jìn)行預(yù)測，并將LSTM方法與傳統(tǒng)的AR方法、ARMA方法進(jìn)行了對比實(shí)驗(yàn)。AR和ARMA的建模過程中，同樣先進(jìn)行缺失數(shù)據(jù)處理，然后進(jìn)行差分運(yùn)算將數(shù)據(jù)變?yōu)槠椒€(wěn)序列，進(jìn)而對數(shù)據(jù)集分割，90%訓(xùn)練，10%測試，然后對測試結(jié)果進(jìn)行反差分處理并統(tǒng)計(jì)RMSE。RMSE計(jì)算公式如式(11)所示。

三種方法在姑咱地震臺上的預(yù)測結(jié)果誤差如如圖8和圖9所示；西昌地震臺預(yù)測結(jié)果和誤差如圖10和圖11所示；雅安地震臺預(yù)測結(jié)果和預(yù)測誤差如圖12和圖13所示。模型的預(yù)測誤差與預(yù)測所消耗的時(shí)間如表2所示。

圖8 姑咱地震臺水化溶解氣氡預(yù)測結(jié)果對比

圖9 姑咱地震臺水化溶解氣氡預(yù)測誤差對比

圖10 西昌地震臺溶解氣氡預(yù)測結(jié)果對比

圖11 西昌地震臺溶解氣氡預(yù)測誤差結(jié)果對比

圖12 雅安地震臺溶解氣氡預(yù)測結(jié)果對比

圖13 雅安地震臺溶解氣氡預(yù)測結(jié)果對比

從圖8～圖13和表2可以看出：①三種模型中，LSTM模型所花費(fèi)的時(shí)間比較多，這是由于LSTM模型是非線性模型，模型的復(fù)雜度和參數(shù)數(shù)量遠(yuǎn)遠(yuǎn)高于AR、ARMA方法。模型的訓(xùn)練過程中需要進(jìn)行大量的迭代，因此消耗時(shí)間較多；②在三個(gè)數(shù)據(jù)集上，對比三種模型實(shí)驗(yàn)結(jié)果，可以看出LSTM模型預(yù)測的誤差(RMSE)最小，AR模型的RMSE最大，ARMA次之。也就是說，LSTM的預(yù)測精度更高，效果最好。

表2 三個(gè)數(shù)據(jù)集上LSTM 與AR、ARMA的性能對比(RMSE)

3 結(jié)論

氡廣泛存在于地殼中，在地震孕育、發(fā)生核構(gòu)造過程中，地殼介質(zhì)受到力的作用而發(fā)生形變，賦存于介質(zhì)中的氡逃逸出來，引起地下水中氡濃度發(fā)生變化，溶解氣氡觀測是目前地震前兆觀測最主要的手段之一。研究溶解氣氡變化規(guī)律，對溶解氣氡進(jìn)行預(yù)測，對快速檢測氡異常、研究震-氡機(jī)制意義深遠(yuǎn)。LSTM作為目前最流行的深度學(xué)習(xí)網(wǎng)絡(luò)模型，在自然語言處理等領(lǐng)域取得了巨大的成就，將LSTM技術(shù)引入溶解氣氡時(shí)間序列預(yù)測，拓展了LSTM的應(yīng)用領(lǐng)域，也為地震前兆數(shù)據(jù)分析引入了新思路；將LSTM與傳統(tǒng)的AR、ARMA方法進(jìn)行了對比實(shí)驗(yàn)。實(shí)驗(yàn)結(jié)果表明，LSTM方法作為一種新興的非線性時(shí)間序列分析方法，雖然耗時(shí)較多，但是與傳統(tǒng)的AR方法、ARMA方法相比，預(yù)測精度明顯優(yōu)于傳統(tǒng)方法。