刘海军, 单维锋, 耿贵珍

(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方法相比，预测精度明显优于传统方法。