基于溶解氡數據和長短期記憶網絡的地震預報

劉海軍, 單維鋒, 耿貴珍

(1.防災科技學院應急管理學院，三河 065201；2.防災科技學院經濟管理學院，三河 065201)

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

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

1 模型與算法

1.1 圖問題數學描述

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

1.2 數據預處理方法

溶解氣氡數據預處理分為兩部分：缺失數據處理和平穩性處理。

(1)缺失數據處理

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

(2)平穩性處理

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

Δxt=xt-xt-1

(1)

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

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

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

1.3 長短期記憶網絡模型

1.3.1 長短期記憶網絡基本結構

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

xt為t時刻的輸入數據，ht為t時刻的輸出數據，A為LSTM處理單元

從圖2可以看出，在LSTM神經網絡中核心部分為LSTM處理單元A。LSTM處理單元A詳細結構如圖3所示。

圖3 LSTM處理單元詳細結構

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

LSTM處理單元中第一個“門”結構是遺忘門。遺忘門決定LSTM單元從單元狀態c中忘記哪些信息，它檢查來自前一個LSTM處理單元的輸出向量ht-1和當前時刻的輸出向量xt，并輸出0～1的數，0表示徹底忘記，1表示完全保留。遺忘門的激活向量ft計算公式如下：

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

(2)

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

(3)

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

(4)

(5)

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

(6)

(7)

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

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

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

(8)

ht=ot*σc(ct)

(9)

至此，t時刻的輸入數據xt和前一時刻的輸出數據ht-1，通過遺忘門、輸入門、輸出門結構，得到了t時刻的輸出數據ht。

從整個LSTM處理單元數據流動可以看出，在LSTM結構中，隨著時間流動的除了輸出數據ht之外，還有單元狀態ct。從圖3中可以看出，LSTM單元的狀態ct就像是一個傳送帶，將遺忘門的激活向量ft、輸出門的激活向量it和輸出門的激活向量ot串聯起來，傳遞向下一個處理單元，這是LSTM的核心。

1.3.2 長短期記憶網絡訓練

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

(10)

將網絡在訓練集上的誤差L作為損失函數。進而利用誤差反向傳播算法求解使得在訓練集上損失函數L取得最小值的系數W、U和b為最終訓練好的網絡參數。

2 實驗與結果分析

2.1 數據和實驗流程介紹

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

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

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

原始數據經過缺失數據、平穩性檢驗、差分處理之后，對數據進行建模并預測，然后進行反差分處理，并統計誤差。實驗流程如圖4所示。

圖4 基于LSTM的溶解氣氡預測實驗流程

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

圖5 論文中使用的原始觀測數據

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

圖6 差分處理后的數據

從圖6中可以看出，經過差分處理之后，數據的均值和方差不隨時間發生變化，變成平穩的時間序列數據，可以進行時間序列數據建模。

2.2 LSTM網絡結構和參數

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

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

表1 LSTM訓練參數

(11)

LSTM神經網絡的訓練過程如圖7所示。圖7中顯示了隨著迭代的進行，均方根誤差的變化曲線。網絡訓練初始階段，均方根誤差的數值較高，隨著迭代的進行，均方根誤差不斷減小，初始時，減小速度較快，隨著迭代的進行，均方根誤差減少速度變緩，當迭代次數達到60以后，均方根誤差趨于平穩，達到飽和，此時，神經網絡訓練完畢，可以用該網絡對訓練數據進行預測。

圖7 LSTM訓練過程中RMSE的變化趨勢

神經網絡訓練好之后，可以應用該網絡對未知數據進行預測。在預測階段，將10%的測試數據集輸入訓練好的LSTM網絡，得到預測數據，然后進行反差分處理，反差分處理之后的數據為最終的預測值。

2.3 實驗結果

應用前面介紹的參數訓練網絡，得到最優模型并在測試集上進行預測，并將LSTM方法與傳統的AR方法、ARMA方法進行了對比實驗。AR和ARMA的建模過程中，同樣先進行缺失數據處理，然后進行差分運算將數據變為平穩序列，進而對數據集分割，90%訓練，10%測試，然后對測試結果進行反差分處理并統計RMSE。RMSE計算公式如式(11)所示。

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

圖8 姑咱地震臺水化溶解氣氡預測結果對比

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

圖10 西昌地震臺溶解氣氡預測結果對比

圖11 西昌地震臺溶解氣氡預測誤差結果對比

圖12 雅安地震臺溶解氣氡預測結果對比

圖13 雅安地震臺溶解氣氡預測結果對比

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

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

3 結論

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