地震是破坏力极大的自然灾害，会引发严重的人员和财产损失。地震预测，尤其是针对6级以上大震的预测，有利于提前部署城市和小区规划，促进重大基础设施的提前完善，从而有效保护人类的生命安全以及减少经济损失。

自1909年起，人类对地震预测的研究已经持续了一个多世纪^[1]。Allen^[2]指出，地震预测涉及地震发生的时间、地点、震级和地震发生的概率等信息。许多研究者曾尝试预测地震，但由于地震的孕育和发生的物理机制太过复杂，并且与地震预测相关的前兆信息也难以识别，导致目前仍然没有找到有效的地震预测方法。

前人采用了许多方法来预测地震，基本可以分为经验预测、统计预测和物理预测等。地震经验预测长期以来是实际预测中常用的基础方法，利用氡气含量、动物异常行为等经验方法进行短临预报，虽然取得了一些有用的成果，但迄今为止没有突破性的进展，主要问题在于此类经验方法含有许多主观因素。利用统计方法则可以减少主观因素影响，并客观严格地总结地震经验。随着20世纪90年代以来地震观测数字化仪器的发展，地震资料海量增长，如何有效地将大量地震资料应用于地震预测成为国内外研究者的主要研究方向。

如今，机器学习在分析及处理大数据方面表现优异，已经吸引了大量的关注^[3-6]。例如，Negarestani等^[7]使用反向传播神经网络(back propagation neural network, BPNN)检测地震引起的氡浓度异常行为；Liu等^[8]利用径向基函数(radial basis function, RBF)神经网络集成预测了中国各个地方的地震；Mirrashid^[9]选择伊朗部分地区[56°E~58°E，27°N~31°N]研究地震预测模型，使用了模糊逻辑、网格划分、减法聚类和模糊c均值等方法，将地震预测视为一个二元分类问题；Ikram和Qamar^[10]开发了一种基于专家系统的地震预测软件；Narayanakumar和Raja^[11]利用8个地震活动性指标作为输入数据，使用BPNN预测喜马拉雅地震带发生地震的概率，并评估了预测效果，提出地震震级在4~6级之间预测效果更好；Asim等^[12]将神经网络和模糊逻辑相结合，使用模式识别神经网络、循环神经网络(recurrent neural network, RNN)、随机森林和线性规划boost集成分类器，探索地震参数与地震发生的关系；Huang等^[13]提出一种利用历史地震事件进行连续地震预测的卷积神经网络深度学习方法；Panakkat和Adeli^[5]计算了8个地震指标，基于BPNN、RNN和RBF神经网络预测下一个月最大震级；紧接着，Panakkat和Adeli^[14]改进了RNN，根据震中位置的经纬度以及随后地震发生的时间预测了南加州地震发生的时间和位置；Wang等^[15]使用长短时记忆(long-short-term memory，LSTM)神经网络了解不同地点地震之间的时空关系，并基于这种关系进行地震预测；李林芳等^[16]使用LSTM神经网络对2008年汶川M8.0地震进行了回顾性预测；Cheng和Shi^[17]使用LSTM神经网络基于16个地震活动指标预测南加州地震，但在预测过程中容易出现过拟合问题。

根据地震的预测时间长度，地震预测一般可分为长期、中期、短期及临震预测。本文侧重于中期预测，预测未来一年可能发生的地震。选取具有长时间记录的南加州地区的地震目录，以数据集中过去多个时间窗口的数据特征(11个预测地震活动性指标)作为输入，以最大震级为标签，基于改进LSTM神经网络构造输入到输出的映射关系并建立地震预测模型，同时使用当前时间窗口的数据测试模型预测效果。

1 研究区域与数据预处理 1.1 研究区域

美国南加州地区处于北美板块和太平洋板块的板块边界，地震活动由板块边界运动、地壳扰动和断层活动引起。该地区由多个地质构造单元组成，包括圣安德烈斯断层带、圣巴巴拉板块、索尔顿板块等。其中圣安德烈斯断层带是南加州最活跃的地震带，该断层带是由太平洋板块和北美板块相互碰撞形成的，长达800英里，分布于加州和墨西哥边境，从圣迭戈一直延伸到圣弗朗西斯科湾，每年都会发生数百次小型地震，这种位于板块边界活动的地区通常可能发生大型地震。

南加州地震目录记录的时间跨度长达91年，地震记录相对完备，是一个较为理想的地震研究区域。故本研究选择南加州地区作为研究区域。图 1绘制了地理范围在[32.3°N~37°N]和[114°W~122°W]的南加州区域构造以及地震目录中记录地震的震级和发生位置。

地震目录来自南加州地震数据中心。本研究使用自1932年1月1日至2021年9月19日的目录，研究区域内一共有831 906条数据。该时段内目录包含6级以上地震26个，7级以上地震4个。表 1列出研究区域内6级以上的地震，并标出了7级以上地震的地点。由表 1可知，最近30年来7级以上地震频度显著升高。

1.2 地震目录完备性分析

古登堡-利特克(Gutenberg-Richter, G-R)公式描述了地震频度与地震震级的关系，即随着震级逐渐增大，地震频度会出现指数级的下降。公式的数学表达形式如下

$\log N = a - bM.$

(1)

其中: N >0，N代表震级大于等于震级M的地震频度；a、b为常数。

图 2绘制了南加州地区1932年1月到2021年9月地震频度N的对数logN和震级M的关系，用来初步估计地震目录完备情况。

从图 2可以看出，震级在小于3时，实际观测频度开始低于式(1)的拟合直线，说明震级在3.0以下的地震目录存在小震缺失的情况。这是由于地震台网侦测能力有限，小震记录不全。因此，地震目录的完备性震级为3.0。并且图 2中大震端也有起伏，这是由于大地震数目较少的缘故。

由式(1)可知, b值反映了震级和频度关系。b值越小，表示较大震级的地震发生频次越高。一般认为b值与介质均匀程度、应力水平有关^[18]。b值一般可以通过2种方法计算，分别为线性最小二乘法(least square method, LSM)和最大似然估计法(maximum likelihood estimation method, MLEM)。

LSM：

$b = \frac{{\sum {_{i = 1}^m{M_i}\sum {_{i = 1}^m\log {N_i} - m\sum {_{i = 1}^m{M_i}\log {N_i}} } } }}{{m\sum {_{i = 1}^mM_i^2 - {{\left( {\sum {_{i = 1}^m{M_i}} } \right)}^2}} }}.$

(2)

MLEM：

$b = \frac{1}{{\ln 10\left( {\overline M - {M_{\rm{c}}}} \right)}}.$

(3)

其中：M是平均震级，M_c是最小完备性震级^①，m是地震发生总次数。图 3绘制了1932—2021年期间不同起算震级下利用MLE法得到的b值，用来精确估计地震目录完备的阈值震级。纵轴是b值，横轴是计算b值时用到的阈值震级M₀。作为地震活动性的表征，一个区域长期的b值应该是一个常数，之所以在低阈值震级时出现低值，是因为小震漏记。

① 完备性震级 M_c 指的是震级大于 M_c 的所有地震都能被有效检测和记录到。

从图 3可以得出，3.2级以上的地震没有漏记。所以本文后续数据处理选取的最小震级为3.2，3.2级以上的地震共16 389个。

1.3 地震活动性指标

本文目标是基于地震目录预测研究区域内某时间段内的最大震级，但仅依靠地震目录不足以直接预测震级，需要从地震目录挖掘出能反映地震时空强度分布特征的信息，本文将这些信息称作地震活动性指标。原则上，如果有几十年来系统完备的前兆资料，也可以将其作为地震活动性指标之一，但是由于缺乏这样的资料，本文仅采用“以震报震”的方法预测(基于历史地震目录预测未来最大震级)。

前人已经尝试使用多种地震活动性指标来预测地震。例如Zamani等^[19]结合基于PCA的数据归一化和特征选择，将此算法用于预测伊朗地震的一组特定地震活动指标。

表 2描述了11种地震活动性指标，这些地震活动性指标是由地震目录计算得到，用以预测未来时间窗口内的最大地震震级。地震活动性指标不仅包括Panakkat等使用的最大震级M_max、频度frequency、均方根误差η、最大震级欠缺ΔM、地震能量平方根$\sqrt E$和平均震级M_mean ^[5]，还包括纬度的均方根误差RMSE_Lat、经度的均方根误差RMSE_lon、按能量加权计算的震中平均纬度LAT_{$\sqrt E$}、按能量加权计算的震中平均经度Lon_{$\sqrt E$}、用最小二乘法计算的经纬度之间的斜率k，这些地震活动性指标能够体现地震发生的空间位置、地震成团或成带分布等^[16]。

2 模型建立与预测 2.1 LSTM神经网络结构

本文选择使用LSTM神经网络作为预测模型，原因是它在处理时间序列数据时具有更好的记忆能力。这一特点已经在处理与时间相关的序列数据(语音或视频)中得到了验证。图 4是LSTM神经网络的基本单元结构原理图。

图 4中，X_t为t时刻采集的输入特征量；y_t-1为t-1时刻的网络输出值；c_t-1为t-1时刻的细胞单元状态^[20]。LSTM网络基本单元中包含遗忘门、输入门和输出门3种门控结构^[20]。其中，遗忘门的输出值直接影响到细胞单元状态是否保留；输入门包含了2种激活函数对特征量的归一化^[20]。

2.2 具体流程

为满足地震活动性指标的计算条件，将输入时间窗口长度设置为1年，按月滑动，输出时间窗口设置为1年，这样的时间窗口设置可以保证每个窗口包含的地震数目在30个以上(地震数目最少的一个时间窗口是49个地震数目)。

全部资料包含1 053个过去时间窗口，把它们分成训练段和测试段，训练段包含842个时间窗口，测试段包含211个时间窗口。

为预测未来一年的最大震级，选取11个地震活动性指标作为输入，时间窗口长度也设为1年。在处理特征的时候，我们考虑了时间序列问题。不仅使用当前时间窗口的数据进行预测，还使用过去数个时间窗口的数据特征作为模型的输入。经过反复试验，最终将使用的过去时间窗口的个数(lookback)设置为6，并将数据构建成一个三维的数据(datalength，lookback，features)，以输入到LSTM神经网络模型中，这样可以兼顾多个时间点的特征。

图 5是该方法的示意图。该程序从地震事件(1932年以来发生在南加州地区)的目录开始。通过不同的输入和输出(预测)时间窗口，计算生成一组11个地震活动性指标，并使用LSTM神经网络对于时序数据进行训练和测试。最后，通过绘制图形和计算误差及均方误差等指标来评估模型的性能。

2.3 预测模型

用1932年1月开始的地震目录数据建立LSTM神经网络预测模型。为了验证其对地震的中期预测效果，设置了如下模型。

模型1：使用整个区域的数据，地理范围在[32.3°N~35°N]和[114°W~122°W]。训练集时间段和测试集时间段之比是0.8∶0.2，训练集的数据是从1932年1月到2002年3月，验证集的数据是从2002年3月到2021年9月。输入和输出的时间窗口都是12个月。

模型2：在模型1的基础上将输入的时间窗口改为24个月，探索输入更长时间窗口对于模型的训练是否有帮助。

模型3：在模型2的基础上将区域范围扩大到地理范围[32°N~35°N]和[114°W~122°W]，这样囊括了1934年12月的6.4级和6.3级大地震、1980年6月的6.3级大地震和2010年4月的7.2级大地震，可以探索大地震数目增多对于模型的训练是否有正面影响。

LSTM神经网络输出层有1个节点，即未来一年内的最大震级。输入包括整个区域在特定时刻前的6个时间窗口内(每个时间窗口大小为12个月)的11个预测因子。经过多次优化寻找设置最合适的神经网络参数，最终将模型1、2、3的隐藏层层数分别设为2、4、4层, 隐藏层内神经元数量都设为36，并且为防止过拟合，均设置了早停算法。经过学习训练和检验, 获得初步结果。

图 6是模型1、2、3整个数据集上的结果。纵坐标是震级，横坐标是年份。在地震预测研究中, 能否预测到6级以上的大地震对模型的发展和改进至关重要。可以看出，3个模型对地震的趋势预测均较为准确，尤其是对1952和1992年的7级地震拟合效果较好。相较于模型1来说，模型2和3的大震预报准确率更高，但是虚报数量也较多。

可以看出，模型1预测出了2004年的M6.0的地震。模型2中虽然测试集部分存在一定的虚报，但是预报出了2003、2004年的大地震。尤其是模型3预测出了2010年M7.2的大地震，可以看出增加大地震的样本有利于大地震的预测。

图 7中直方图的横轴是预测到的震级减去实际观测到的震级，画出了误差的柱状图。

从直方图中可以看出，大部分误差是落在0.5个震级以内，说明预测结果还是有一定可信度的。并且模型2和模型3的误差要显著小一些，说明增大输入时间窗口的跨度对于震级预测尤其是大地震预测有着正面效果，未来可以朝着这个方向继续努力。

图 7中的散点图是描述预测的震级和观测的震级之间的关系，并且画出了预测震级与实际震级相等的辅助线，以及误差为0.5的辅助线。可以看出，对于6级、7级以上的地震而言，相较于模型1，模型2的表现更好，说明增加时间窗口跨度有利于预测大地震。同时相较于模型2，模型3的表现更好，说明增多大地震样本的输入是有效的。

2.4 结果分析

训练后，必须评估预测记录和实际记录之间的差异。为了衡量我们模型的性能，考虑引用混淆矩阵。混淆矩阵的4个指标可以描述如下：

TP：算法预测发生地震并且实际发生地震的地震事件数量。

TN：算法既没有预测发生地震，实际也没有发生地震的地震事件数量。

FP：算法预测发生地震，但实际并未发生地震的地震事件数量。

FN：算法没有预测发生地震，但实际发生地震的地震事件数量。

由于地震预测关注的是大地震的预测，所以表 3着重列出了训练集和测试集的6级和7级以上地震的混淆矩阵的数值。

通常，在混淆矩阵中，TP和TN的数量越大，FP和FN的数量越小，模型的性能越好。然而，在实际研究中很难通过混淆矩阵来衡量模型的质量。因此对于标准化测量，基于混淆矩阵会额外计算另外几个度量，其中除马修斯相关系数(Matthews correlation coefficient, MCC)外的度量值均介于0~1(表 4)，越接近于1，越是理想状态。MCC是在-1~1之间摆动，系数1表示完美预测，0是平均随机预测，-1表示预测和观察之间完全不一致。

根据表 4中的公式，分别根据训练集和测试集的混淆矩阵的数据，计算得出表 5中的结果。

从表中的计算结果可以看出，3个模型的训练集以及测试集的accuracy值是比较接近1的，说明预测正确的比例比较高。但测试集的准确率比训练集的要略小一些，这也是正常的，可能存在过拟合的情况。但是更大原因是由于样本数据并不足够多，尤其是大地震，同时我们选取的地震活动性指标并不能完全准确地把握地震规律，还需要大量的实践研究和艰苦探索。P₀、P₁和S_p的值是比较接近1的，尤其是模型1预测大于7级地震时候的S_p值等于1，是理想状态，说明虚假预测地震的概率很小。模型1的S_n数值比较小，结果不太理想，但是模型2和模型3的S_n数值比较大，尤其是模型2，说明模型2的灵敏度最高。4个指标的算术平均值avg还是比较接近1的。模型1和模型3的MCC只是略大于0，说明只是比平均随机预测稍微好一些，但是模型2的MCC略大一些，说明模型2的整体预测水平更高。

总体而言，3个模型在测试集上的表现各有长处，模型2综合看来表现更优一些。

为进一步评估LSTM神经网络的准确性和可靠性，本文列出了神经网络2个常用的指标：均方误差(mean-square error, MSE)和均方根误差(root mean square error，RMSE)，以及它们的计算公式。

${\rm{MSE}} = \frac{1}{n}\sum\nolimits_i^n {{{\left( {\hat M_i^{{\rm{future}}} - M_i^{{\rm{future}}}} \right)}^2} + \alpha \left\| \boldsymbol{W} \right\|} _2^2，$

(4)

${\rm{RMSE}} = \sqrt {\frac{1}{N}\sum\nolimits_{i=1}^N {{{\left( {\hat M_i^{{\rm{future}}} - M_i^{{\rm{future}}}} \right)}^2}} } .$

(5)

其中N是训练的样本总数。我们的结果是：模型1训练集上：MSE=0.028 5，RMSE=0.168 8。测试集上: MSE=0.048 6，RMSE=0.220 4；模型2训练集上：MSE=0.030 0，RMSE=0.173 1。测试集上: MSE=0.040 4，RMSE=0.201 0；模型3训练集上：MSE=0.020 7，RMSE=0.144 0。测试集上: MSE=0.049 8，RMSE=0.223 1。可以看出模型1和模型3的2个神经网络常用指标相差不大，在训练集上模型2表现好一些。整体看来测试集的误差比训练集的误差稍微大一些，但是整个神经网络的均方误差还是不算大的。

中国地震局的年度地震预测通过R评分进行评估(完全随机猜测的R得分约为0，完全成功的预测的R得分约为1)。预测效能R评分的计算公式一般为

$R = \frac{{{\rm{TP}}}}{{{\rm{TP + FN}}}}{\rm{ - }}\frac{{{\rm{FP}}}}{{{\rm{FP + TN}}}}.$

(6)

本次结果，模型1和模型3的R评分比较低，在0.20左右，虽然略高于20世纪90年代中国地震局年度预测的R评分0.18^[21]，但与近年中国地震局年度预测水平大体相当。而模型2的训练集对于6级以上地震的R评分是0.69，模型2的测试集对于6级以上地震的R评分是0.62，高于近年中国地震局年度预测水平。

3 讨论和结论

南加州地区处于断层带，分布有两百多条大小不等的断层，是地震活动频繁区，具有1932年以来3级左右以上地震基本完备的地震目录，适合作为地震预测实验研究的数据基础。通过对该地区的地震目录进行数据处理，建立LSTM神经网络预测模型，利用地震活动性指标以及不同的时间输入窗口，进行未来一年的最大震级的地震预测。研究发现，LSTM神经网络在对南加州地震预测应用方面有一定可行性，也有一定的效果。尤其是模型2的测试集对6级以上地震预测的R评分为0.62，高于中国地震局近年年度会商的R评分。

利用LSTM神经网络预测模型在南加州地区进行近20年回溯性预测大地震的检验，在模拟训练中发现，仅仅利用1932—1999年的地震目录，凭借着地震目录进行“以震报震”，能够回溯性地预测2003年的M 6.5地震和2004年的M 6.0地震，特别是模型3预测出了2010年4月的M 7.2大地震。所以依靠人工智能神经网络进行回溯性的预测，对于以后的地震预测研究具有一定意义和启发。但遗憾的是在2018年7月到2019年6月并没有预测出2019年7月有发生7级以上大地震的趋势。这是否是因为2019年的地震发震断层在构造上处于东西走向的加洛克断层北部，与圣安德列斯走滑断层体系存在某种程度的隔离，所以与之前的地震目录并不存在物理上的一些规律，尚不清楚。还需要探索加入更多的地震活动性指标，以及进行参数调节和尝试使用更多更复杂的神经网络。

目前本文使用的LSTM神经网络预测模型，可以在几分钟内完成几十年的地震目录计算和数据训练，并输出预测结果，效率很高。因此，今后有条件可以对更多地震活动区域尝试本文的方法，并探索新方法。探索中可以对不同时段的训练和测试效果进行评估，改变输入时间窗口长度，在评价预测效果时通过对输出大震震级加权以增大大震预测效果等，探索是否能在目前基础上有所提高。

由于神经网络具有“黑匣子”属性，即使地震预测趋势得到较好拟合，但其中具体的影响因素和内在关联也难以解释。因此，未来还需要在以下两方面进行深入探讨: 在神经网络预测研究中，要比较和分析各个因子的贡献，逐步剔除贡献不大的因子，扩展相关关系较大的因子；在地震孕育发生的物理研究中，深层次地挖掘孕育地震背后的物理原因，并利用人工智能搭建框架，添加更多的地质数据、构造信息以及地震前兆观测等新的数据，添加具有物理意义的新预测因子，期待预测成果不断更新提升，并且能够应用于实践。