长短时记忆神经网络在中期地震预报中的探索——以川滇地区为例

李林芳，石耀霖，程术

中国科学院大学地球与行星科学学院，计算地球动力学重点实验室，北京 100049

0 引言

自1970年以来，中国大陆区域6.0级以上震级强震发生了200余次(熊子瑶，2020).这些地震造成严重的人员伤亡，破坏建筑物和基础设施，导致了巨大的经济损失.时间尺度为年度至数年的中长期地震趋势预测，对减轻地震灾害风险具有重要意义，同时也是一个复杂的科学问题(马宗晋，1980；张国民和李宣瑚，1988；Knopoff，1996；Sobolev，2015).

人工智能及大数据的迅猛发展，在许多行业领域里引起了广泛关注.机器学习作为当今世界发展最快的技术之一，是人工智能及数据科学的核心(Jordan and Mitchell，2015).机器学习最基本的做法，是使用算法来解析数据、从中学习，然后对真实世界中的事件做出决策和预测.随着我国地震台站监测和前兆台站监测的增加，地震数据量不断增多，大数据有利于挖掘其中隐藏的规律和趋势，提高地震预测的可靠性(张晁军等，2015).

早在20世纪90年代，机器学习已被用于进行地震趋势预测.M8系列算法(M8，MS，CN)利用地震目录资料，基于模式识别手段对强震发生概率增长的时段进行判定预测(Keilis-Borok and Kossobokov，1990；Tiampo and Shcherbakov，2012).黄德瑜等(1995)运用CN算法对1970年以来的23次强震进行回溯性预测，取得了较好的地震预测效果.近些年来，国内外学者将人工神经网络用于中、长期地震预测，取得了不错的效果(Panakkat and Adeli，2007；Adeli and Panakkat，2009；聂红林等，2012；Asim et al.，2018；Huang et al.，2018；Wang et al.，2020；Li et al.，2020).Mignan和Broccardo(2020)综述了1994年以来的77篇有关运用人工神经网络进行地震预测的文章，认为将人工神经网络应用于地震预测领域的热度正在逐渐增加，且人工神经网络趋于更复杂的长短时记忆神经网络(LSTM)和卷积神经网络(CNN).

川滇地区(21°N—36°N，96°E—108°E)位于青藏高原东南缘，是欧亚板块与印度板块相互作用的边缘地带(图1).作为中国大陆地震活动最强烈的地区之一，川滇地区成为了地震学家的重点关注对象.Li等(2017)通过三维黏弹性有限元模型对川滇地区地壳变形进行了数值模拟，并与GPS观测结果对比，得出下地壳的流动在川滇地区地壳变形中起着重要作用；Jiang和Wu(2010)基于图像信息学(Pattern Informatics，PI)算法，在特定的目标震级和时空窗口下实现了对汶川地震的回溯性预报，但目标震级及网格的增加会导致回溯性预报失败；程佳等(2020)采用地震危险性标准计算方法，建立了考虑超过川滇地区历史最大震级的地震危险性预测模型；Wang 等(2015)综合地震目录、地壳应变速率、地震矩累积速率对川滇地区的地震危险性进行评估，认为多条断裂带可能发生M>7.5地震.

本文选择地震活动频繁的川滇地区作为研究对象，基于地震目录数据建立长短时记忆神经网络模型，运用该模型进行回溯性地震预报，取得理想的效果，并给出2021年川滇地区地震预测参考意见.

1 数据

1.1 地震资料

本文数据来源于中国地震台网中心，选取川滇部分地区(24°N—32°N，98°E—106°E；图1)1970年1月1日至2019年12月31日的地震目录.该段时间内发生地震事件496846个，最大震级为MS8.0.其中，6.0～6.9级地震事件45个，7.0级及以上地震事件7个，如表1所示.

表1 川滇地区1970年以来7级以上地震(标粗为发生在研究区的地震)Table 1 Earthquakes of more than M 7 in the Sichuan-Yunnan region since 1970 (The bold mark is the earthquake that occurred in the study area)

图1 川滇地区区域构造背景及地震活动图(根据Li et al.，2017修改)Fig.1 Regional geologic tectonics and earthquake events in Sichuan-Yunnan region (modified from Li et al.,2017)

1.2 最小完备震级检测

影响地震数据好坏的一个很重要的因素就是地震最小完备震级Mc的选取.地震最小完备震级是指Mc以上的地震均被仪器记录到，无遗漏及不完备现象(苏有锦等，2003).

检测地震完备性的步骤是：计算不同Mc时的最大似然b值，计算不同Mc对应b值的标准差cb，当b的差值Δb小于误差cb可以判断为最小完备震级.

最大似然法具有计算简便、不易受个别较大地震影响、计算结果较为稳定的优点(王光明等，2020)，具体计算公式为：

(1)

(2)

式中，N是计算b值所用的样本量.

以Mc为横轴，b为纵轴做曲线，刚开始b随Mc增大而增大，当b基本变平，再利用b的差值Δb小于误差σb可以判定为最小完备震级.1970—2019年b-Mc变化趋势图如图2所示.

图2 1970—2019年研究区Mc-b趋势图Fig.2 Trend chart of Mc-b in the study area from 1970 to 2019

经计算，1970—2019年研究区的最小完备震级Mc为3.0.研究区在1970—2019年发生了13425次3级以上地震.

2 方法

2.1 长短时记忆神经网络

长短时记忆神经网络(Long Short-Term Memory，LSTM)是一种特殊的循环神经网络(Recurrent Neural Network，RNN)，最初由Hochreiter和Schmidhuber(1997)提出，后来由Gers等(2000)进行了优化.LSTM提升了神经网络的长时间记忆能力，适用于有关时间序列数据的记忆，近些年来已被广泛应用于解决文本识别(Liu and Guo，2019)、动作感知(Liu et al.，2018)、交通预测(Ma et al.，2015；Zhao et al.，2017)、金融市场预测(Fischer and Krauss，2018)、降水预测(Kratzert et al.，2018)、地下水泉涌量预测等问题(Cheng et al.，2021).

作为一种特殊的RNN，LSTM神经网络由一个输入层、一个输出层和一系列记忆单元循环相连而成的隐藏层组成.记忆单元结构如图3所示，它由一个或多个自我循环的记忆细胞和三个门(输入门it、输出门ot、遗忘门ft)组成，这些机制为记忆细胞提供连续的读、写和重置操作.在t时刻，输入给该记忆单元的为t-1时刻的隐藏层状态变量ht-1、记忆单元状态变量Ct-1和t时刻的信息Xt，经过遗忘门、输入门、输出门后获得t时刻的隐藏层状态变量ht和记忆单元状态变量Ct，最终ht会进入输出层生成LSTM在t时刻的计算结果yt，同时与ct一起传入下一个记忆单元进行计算.

图3 LSTM神经网络记忆单元结构图Fig.3 The structure of one memory block in LSTM neural network

2.2 LSTM模型输入输出确定

2.2.1 输入的地震预报因子选择

地震预报因子能够用来评估一个区域潜在发生的地震事件.Panakkat和Adeli(2007)选取基于G-R关系的8个地震参数作为神经网络的输入端，Reyes 等(2013)选取G-R关系中b值的增量为神经网络的输入端，Asencio-Cortés 等(2018)综合了前人对地震参数的选择，提出16个地震预报因子，取得了较为理想的预报结果.因此，地震预报因子的选择将直接影响地震预报的准确性和稳定性.

本文基于地震目录，采用16个地震预报因子来进行地震预报(表2).其中时空窗口内的地震频度N，最大地震震级Mmax，基于Giternberg-Richter关系的最小二乘b值和a值是经验预报中经常用到的参量，描述最小二乘偏差特征的最小二乘G-R拟合时的均方差σG-R和最大震级欠缺ΔM也是前人神经网络模型中使用过的特征因子(Asencio-Cortés et al.，2018).最小二乘b值和最大似然b值，虽然在理论上应该一样，在实际计算中总会有一定差异，这种差异是否对地震预报有意义并不清楚，因此我们将二者都包含在内.我们在预报因子选取中故意保持一定的冗余，宁可后续发现效果重复可以剔除，但在早期的试探中不要漏掉可能有用的因子.

表2 基于地震目录的预报因子Table 2 Seismic indicators calculated from seismic catalog

本文还新引入Latmean、σLat、Lonmean、σLon、K、LatE、LonE共7个表征地震空间位置的特征因子.它们反映了地震活动空间位置、地震成团或成带分布情况，能量释放空间位置等.这些参量在某种程度上与我国经验预报中的地震条带、空区等空间分布特征有关.

2.2.2 数据时空滑动窗口的设置

选定预报因子后，还要选取合适的时空窗口来计算预报因子.例如最小二乘b值，如果一个时空窗口内的地震数目太少，计算的b值将激烈起伏极不可靠.本文使用1970年以来的地震目录，早期的地震台网监测能力有限，仅能做到3级以上的地震大体没有遗漏，这样就需要比较长的时间和比较大的空间窗口.一般认为7级左右和7级以上地震，其“临近”前兆可能的持续时间可达年度尺度或更长时间(“地震预测预报二十年发展设计”工作组，2017)，本文选取两年作为地震数据的时间窗口，来预报未来一年内发生地震的最大震级.但为了适应地震部门往往每个月都需要会商震情的实际工作需要，做到每个月都能够补充最新资料和更新分析结果，因此将时间窗口进行滑动，滑动时间为1个月(图4).

图4 地震数据时间窗口分区Fig.4 Temporal window partitioning of seismic data

此外，还需要确定地震数据的空间窗口.首先将川滇地区划分为64个1°×1°的小单元.虽然采用了两年的时间窗口，但若以1°×1°的单元格为空间窗口，则在一个时空窗口内的地震数目仍不满足b值等预报因子计算的要求.因此根据尝试结果，我们采用6°×6°的矩形区块为空间窗口.与滑动时间窗口类似，为了在空间窗口较大的情况下充分利用资料提高分辨能力，空间窗口也进行滑动.从研究区西北角(32°N，98°E)开始向东滑移一格，到区域边缘后再向南一格开始新一行的滑移，最终将研究区划分为9个空间区块(1西北部，2北部，3东北部，4西部，5中部，6东部，7西南部，8南部，9东南部)，如图5所示.经过时间窗口及空间窗口的设置，单个时空窗口内的地震频度至少为107.

图5 地震数据空间窗口分区Fig.5 Spatio window partitioning of seismic data

2.2.3 模型输入和输出

由于LSTM模型的结构特点，生成最终计算结果yt的隐藏层状态变量ht会传入下一个记忆单元进行计算(图3)，这就要求模型的输入端包含前一预报窗口的结果.因此，LSTM模型的输入包括1～9号区在单位时间窗口(两年)内的16个预报因子和前一个预报窗口的最大震级共17个参数，输出1～9号区各自下一年的最大震级.例如：输入1～9号区各自在1970年1月至1971年12月间的16个预报因子和1971年12月至1972年11月的最大震级，输出1～9号区各自在1972年1月至1972年12月的最大震级.

2.3 LSTM模型建立

2.3.1 参数确定

经过多次模型优化，最终将LSTM结构中隐藏层层数设置为1层，隐藏层内神经元数量设置为64.在训练方面，学习率设定为0.001，最大训练代数设定为100代.

2.3.2 归一化与反归一化

LSTM神经网络的输入数据需要进行归一化处理.本文采用MinMaxScaler函数，将输入的16个地震预报因子进行标准化.输出端则需要反归一化处理，采用Inverse_Transform函数将输出的最大震级反归一化，进而与实际发生震级进行比较.

2.3.3 训练集、测试集设置

为了验证LSTM对汶川地震的预报效果，基于川滇地区1970—2019年地震目录数据，本文在设置训练集∶测试集=8∶2(汶川地震位于训练集)的基础上，另设置训练集∶测试集=7∶3(汶川地震位于测试集).在训练集：测试集=8∶2模型中，训练集预报起始时间窗口为1972年1月至2009年6月，测试集预报起始时间窗口为2009年7月至2018年11月；在训练集：测试集=7∶3模型中，训练集预报起始时间窗口为1972年1月至2004年10月，测试集预报起始时间窗口为2004年11月至2018年11月.

2.3.4 评价指标的选取

为了量化LSTM模型进行中期地震预报的效果，本文采用三种评价指标，分别为均方误差(Mean Square Error，MSE)、平均绝对误差(Mean Absolute Error，MAE)和均方根误差(Root Mean Square Error，RMSE).具体计算公式如下：

(3)

(4)

(5)

3 结果

3.1 中期预报结果

使用LSTM神经网络模型进行中期地震预报的结果如图6所示.通常认为预报震级位于真实震级±0.5级范围内的为准确预报(陈运泰，2009).在训练集∶测试集=8∶2模型中，测试集含1017个时空窗口，即预报地震事件1017次，预报最大震级为MS7.44，准确预报地震事件797次，准确预报率为78.37%；在训练集∶测试集=7∶3模型中，预报地震事件1521次，预报最大震级为MS8.46，准确预报地震事件1068次，准确预报率为70.22%.由表3可知，训练集∶测试集=7∶3模型的误差要大于训练集∶测试集=8∶2模型.

图6 LSTM测试集预报结果(a)训练集∶测试集=8∶2；(b)训练集∶测试集=7∶3.Fig.6 The forecast results of test set in LSTM(a)Training set∶test set=8∶2;(b)Training set∶test set=7∶3.

表3 不同模型的预报结果评价指标统计Table 3 Forecast result evaluation index statistics of different models

将预报震级-真实震级>0.5定义为虚报，预报震级-真实震级<-0.5定义为漏报.地震学界长期存在“宁可虚报也不漏报”和“宁可漏报也不虚报”2种截然相反的主张和争论(吴中海和赵根模，2013).本文LSTM神经网络给出的地震预报结果，其虚报、漏报和准确预报分布情况如图7所示.在训练集∶测试集=8∶2模型中，虚报127次，虚报率为12.49%，漏报93次，漏报率为9.14%；在训练集∶测试集=7∶3模型中，虚报285次，虚报率为18.74%，漏报168次，漏报率为11.05%.无论是虚报率还是漏报率，训练集∶测试集=7∶3模型均差于训练集∶测试集=8∶2模型.这显示随着数据的积累，有更多年份的资料可以用于学习训练，预测的质量有望得到提高.

图7 LSTM测试集预报情况(a)训练集∶测试集=8∶2；(b)训练集∶测试集=7∶3.横轴为先空间窗口、再时间窗口序列.即按照时间窗口依次显示1～9各区的预报情况.Fig.7 The forecast distribution of test set in LSTM(a)Training set∶test set=8∶2;(b)Training set∶test set=7∶3.The horizontal axis is the sequence of space window first,then time window.That is,according to the time window,the forecast distribution of subarea 1～9 is displayed in turn.

地震预报中最关键的是对破坏性地震(6级以上)的预报.LSTM神经网络对破坏性大地震的预报结果如表4所示.

表4 LSTM神经网络对破坏性大地震预报结果Table 4 Destructive earthquake forecast results of LSTM neural network

在训练集∶测试集=8∶2模型中，发布6级以上预报340次，其中准确预报225次，预报成功率为66.18%，虚报66次，虚报率为19.41%，漏报49次，漏报率为14.41%.因此，当该模型发布6级以上预报时，80.59%的可能发生预报震级或更大震级的地震；发布7级以上预报52次，其中准确预报48次，预报成功率为92.31%，虚报4次，虚报率为7.69%，漏报0次.因此，当该模型发布7级以上预报时，92.31%的可能发生预报震级或更大震级的地震.

在训练集∶测试集=7∶3模型中，发布6级以上预报455次，其中准确预报259次，预报成功率为56.92%，虚报157次，虚报率为34.51%，漏报39次，漏报率为8.57%.因此，当该模型发布6级以上预报时，65.49%的可能发生预报震级或更大震级的地震；发布7级以上预报93次，其中准确预报57次，预报成功率为61.29%，虚报26次，虚报率为27.96%，漏报10次，漏报率为10.75%.因此，当该模型发布7级以上预报时，72.04%的可能发生预报震级或更大震级的地震.

3.2 R值评分

表5 R值评分计算参数Table 5 R-value calculation parameters

R=c-b=d-a,

(6)

R值越大，预报效果越好.随机猜测预报时R为0，完全准确预报时R为1(石耀霖等，2000).为了评价LSTM神经网络对6级以上地震的预报效果，选取预报成功率高的训练集∶测试集=8∶2模型进行R评分计算.由于R值计算通常选用1°×1°网格，需将本文设置的6°×6°空间窗口转换为1°×1°空间窗口，1°×1°网格的预报震级计算公式如下：

(7)

表6 LSTM神经网络训练集∶测试集=8∶2模型地震R评分参数Table 6 R-value of training set∶test set=8∶2 model in LSTM neural network

3.3 地震预测回溯性检验

在地震预测研究中，回溯性检验对模型的发展和改进至关重要(Field，2007).前已述及，训练集∶测试集=7∶3模型的设置是为了对汶川地震进行回溯性检验.LSTM神经网络对汶川地震的回溯性预报结果如图8所示.在本研究中汶川地震(31°N，103.4°E)坐落于西北部、北部、东北部、西部、中部、东部共6个空间窗口内.同时，在2007年6月至2008年5月这一中期预报窗口汶川地震首次出现.当预报窗口起始时间由2007年5月滑动至2007年6月时，除南部和东南部区域外，其余区域的地震危险性急剧升高，其中西北部、北部、东北部、西部4个区域在该时间窗口内(2007年6月至2008年5月)预报震级升至7.5～8级之间.当预报窗口起始时间由2007年8月滑动至2007年9月时，整个研究区的地震危险性再次升高，西北部、北部、东北部、西部、中部、西南部6个区域在该时间窗口内(2007年9月至2008年8月)可能发生7.5～8级地震.综上，在2007年6月至2007年9月这段时间内，LSTM神经网络能够发出大震警报，预测2008年5月至2008年8月在研究区的北-中部可能发生8级左右地震.虽然在1970年到2004年的训练数据内并没有发生过8级以上地震，但LSTM方法还是能够成功回溯性预报MS8.0汶川地震.

图8 汶川地震回溯性检验Fig.8 Retrospective test of Wenchuan earthquake

此外，LSTM神经网络能够成功回溯性中期预报研究区2010年以后发生的6.0级以上地震(图9).从2012年5月开始，预测的最大震级显著增大，特别在西北部、北部、东北部三个地区未来一年会有高于7级的地震发生，到2012年9月六个地区预测未来一年内有6.6级左右地震发生，结果在2013年4月四川省芦山县(30.3°N，103.0°E)发生了7.0级地震；芦山地震后各区域预测震级显著降低，但在2013年9月各区域预测震级再次显著升高，预测2013年9月至2014年8月这段时间可能有6.5级以上地震发生，结果在2014年8月云南省鲁甸县(27.1°N，103.3°E)发生了6.5级地震；震后虽然曲线有所下降，但九个区域平均预测仍有6.2级左右地震会在2014年9月至2015年8月期间发生，结果在2014年11月四川省康定县(30.3°N，101.7°E)发生了6.3级地震；此后曲线一直停留在较低水平，直到2017年6月后逐渐升高，特别在2018年7月后，普遍预测在2019年6月前可能有6.1级左右地震发生，结果在2019年6月四川省长宁县(28.34°N，104.90°E)发生了6.0级地震.

图9 研究区2010年后6级以上地震回溯性检验Fig.9 Retrospective test of earthquakes above M 6 in the study area since 2010

3.4 2021年中期地震形势估计

虽然十多年来陆续有不少采用机器学习方法进行预报的尝试，但是大多停留在研究阶段.我国具有每年年终对下一年地震风险进行预测的传统(Shi et al.，2001)，能否把机器学习的研究成果转变为日常实际预测的一种参考手段呢？我们对此进行了尝试.由于训练集∶测试集=8∶2模型的预报效果更好，将该模型应用到最近时段地震目录资料，到本文写作截止地震目录为2020年1月至2021年1月(实用中可以陆续把新的每个月地震资料补充进来)，据此对2022年2月前可能发生的地震进行预报，其结果如图10所示.当预测起始时间滑动至2021年2月时，即对研究区2021年2月至2022年1月可能发生的地震进行预报，西部、中部、东部、西南部可能发生5.1～5.3级地震.

图10 LSTM神经网络2021年中期预报结果Fig.10 2021 intermediate earthquake forecast results from LSTM neural network

4 讨论

我们曾经使用过几种神经网络方法，对研究区域进行预报研究，结果都没能取得理想结果，而LSTM方法取得了有希望的效果.在运用LSTM时本研究设置的两组训练测试模式，在误差评价指标、准确预报率、R值上，训练集∶测试集=8∶2均优于训练集∶测试集=7∶3，这说明随着时间流逝而数据增加，LSTM神经网络的学习训练效果有可能进一步改进.

大地震的非频发性，是地震预报困难的关键因素之一(吴中海和赵根模，2013).例如美国南加州在1950—1990年共492个月中仅有一个月发生了6.0～6.5级地震，Panakkat和Adeli(2007)采用反向传播(Back Propagation，BP)神经网络和循环神经网络均没能成功预测该地震事件，导致R评分为0.本文采用时空滑动窗口，增加了可供学习训练的地震数据量，例如汶川地震通过滑动时空窗口，在6个空间窗口和12个时间窗口中出现，即在训练中出现了72次，极大的提高了大地震的训练次数.

LSTM神经网络所输入的地震预测因子，作为学习训练的对象，也影响着地震预报的结果.本文为了在早期的试探中不漏掉有效果的因子，神经网络的输入端可能包含了冗余的预报因子.因此，预报因子的增选或剔除仍需进一步研究，可以通过比较剔除某个预报因子前后的R值来判断该因子是否为冗余的.此外，在增选的预测因子中加入前兆资料，尤其是近20年来取得极大进展的GPS/GNSS资料，定量化的综合使用地震活动性“以震报震”和前兆手段，应该可以进一步提高地震预报的可靠性.但前兆资料具有明显的时空分布不均衡性，早期前兆资料匮乏，观测点有限，因此在选择和应用时仍有许多困难需要克服.

不同的LSTM神经网络结构和参数，即模型构建的差异也会导致不同的预报结果.对于模型的层数、各层神经元数量，训练时学习率和最大训练代数的设定等，都还需要进一步试验探讨.除了LSTM神经网络，其他人工智能处理大数据的方法也值得进一步探索，例如随机森林(Rouet-Leduc et al.，2017)、支持矢量机(Asim et al.，2018)、模式识别(Asim et al.，2017)、K邻近算法(Last et al.，2016)等.要通过回溯性预报检验，更要通过未来实用化进行预报检验来不断改进方法.

5 结论

本文选择地震活动频繁的川滇部分地区(24°N—32°N，98°E—106°E)作为研究对象，基于地震目录数据，采用了滑动的时空窗口选取16个反映地震时空强度分布特征的地震预报因子，设置了训练集∶测试集=8∶2和训练集∶测试集=7∶3两个不同的训练测试模型，建立了长短时记忆神经网络模型.运用该模型对研究区9个子区块未来一年的最大地震震级进行预测，得到如下结论：

(1)训练集∶测试集=7∶3模型能够利用1970年1月至2004年9月的地震目录进行学习，在2007年6月至2007年9月间发出大震警报，成功回溯性预报2008年汶川地震；

(2)训练集∶测试集=8∶2模型利用1970年1月至2009年5月资料进行训练，回溯性预报6级及以上地震的R评分为0.407，回溯性7级地震预报时准确率高达92.31%，预报效果显著；

(3)加入最近时段地震目录资料后，预报2022年2月前研究区存在发生5.1～5.3级地震的潜在危险性.

(4)输入端的地震预报因子仍需进一步的改进，通过剔除冗余的因子以及引入可靠的前兆资料作为预报因子，来增加LSTM模型的准确率.

致谢感谢刘杰研究员和中国地震台网中心提供地震目录资料.感谢三位匿名评审专家的中肯意见和建议.