集成深度学习模型在地震烈度预测中的应用——以意大利数据为例

Application of an Ensemble Deep Learning Model in Earthquake Intensity Prediction—A Case Study Using Italian Data

摘要: 本研究探讨了一种基于深度学习的集成策略来预测地震烈度的方法。通过采用意大利地震数据集(INSTANCE)的数据，研究结合了卷积神经网络(CNN)、长短期记忆网络(LSTM)和全连接网络(FCNN)，并运用Bagging算法以提高模型的泛化能力和预测精度。实验结果表明，本研究所提出的集成模型能够有效地预测地震烈度，并对不同烈度级别进行了准确区分。

Abstract: This study explores a deep learning-based ensemble strategy for predicting earthquake intensity. Using data from the Italian earthquake dataset (INSTANCE), the research integrates Convolutional Neural Networks (CNN), Long Short-Term Memory networks (LSTM), and Fully Connected Neural Networks (FCNN), while employing the Bagging algorithm to enhance the model’s generalization ability and prediction accuracy. The experimental results demonstrate that the proposed ensemble model effectively predicts earthquake intensity and accurately distinguishes between different intensity levels.

1. 引言

1.1. 研究背景与意义

地震是自然界中最具破坏力的自然灾害之一，其烈度评估对于提高应急响应效率和减轻灾害损失具有重要意义。

传统的地震监测方法虽然提供了关于地震震源、震级和震源深度等震源参数的信息，但对于地震的烈度评估还存在一定的局限性。

随着地震数据采集技术和计算机计算能力的进步，利用机器学习或者深度学习技术对地震烈度进行有效评估提供了新的可能性。近年来深度学习技术在地震数据处理和特征提取方面表现出了突出的优势。尤其是在地震波形数据的自动识别和分析方面，卷积神经网络(CNN)和长短期记忆网络(LSTM)的结合使用，能有效捕捉波形数据中的时空特征，此外，通过集成学习方法，可以提高预测模型的稳定性和准确性。本研究基于深度学习技术，探索对地震烈度的高效预测，为地震预警与灾害减灾提供新的科学依据。

1.2. 研究思路与方法

本研究旨在通过意大利地震数据集(INSTANCE)中的三分量地震动数据，预测地震烈度。首先考虑到数据集的元数据中不包括烈度指标，研究将采用ShakeMap标准来计算每个地震样本的烈度指标。具体来说，将通过地震波的峰值地面加速度(PGA)、峰值地面速度(PGV)和其他相关参数，结合ShakeMap提供的经验公式，估算地震的烈度等级。

随后，为了实现地震烈度的精确预测，本研究将开发一个深度学习框架，该框架结合了卷积神经网络(CNN)、长短期记忆网络(LSTM)和全连接神经网络(FCNN)。CNN将用于从地震动数据中自动提取空间特征，而LSTM的引入能够捕捉这些特征随时间的动态变化，最后通过全连接层进行特征融合和烈度预测。

此外，为了提升模型的泛化能力并减少预测误差，将采用Bagging集成学习算法。通过在不同的训练子集上训练多个深度学习模型并将它们的预测结果进行集成，可以有效地提高预测的稳定性和准确性。具体而言，每个单独的模型都将基于Bootstrap抽样技术生成的随机训练数据集独立地完成训练过程。最终，通过对各个模型输出的结果采取投票或平均的方式进行汇总，以期获得更加稳定可靠的烈度预测值。

采用这一综合性策略，本研究不仅能实现地震烈度的精准预测，还为地震应急管理部门提供了科学依据，有助于显著降低地震带来的潜在损失。

1.3. 文献综述

地震烈度预测作为地震预警系统的核心组成部分，近年来受到广泛关注。通过分析地震动参数与烈度之间的关系来预测烈度是一种有效的方法。传统研究指出，不同地震动参数与烈度之间的相关性存在显著差异。例如，金星等(2013) [1]提出了一种基于多参数综合评估的地震仪器烈度计算方法，旨在提高地震烈度评估的准确性和适用性。马强等(2014) [2]基于中国大陆地区的地震数据，建立了多种地震动参数与烈度的回归公式，结果表明地震峰值速度PGV与烈度之间具有较好的相关性。李亮等(2018) [3]深入分析了美国地质调查局(USGS)和日本气象厅(JMA)的仪器烈度算法，并验证了这一方法在2008年汶川地震中的应用，显示出较高的准确性和可靠性。

深度学习技术的兴起为烈度预测提供了全新思路。多种高质量地震数据集已被开发并广泛用于机器学习研究。具体来说，Mousavi等(2019) [4]收集的STEAD数据集包括120万条三通道地震波形数据，支持全球范围内的地震信号检测和分析；Magrini等(2020) [5]开发的LEN-DB数据集涵盖了全球范围内的局部地震，包含120万条三通道地震波形数据；Yeck和Patton(2020) [6] [7]的NEIC数据集支持全球实时地震检测，Yeck等学者使用三个独立的卷积神经网络模型训练了130万个地震相位到达，以预测到达时间开始、相位类型和距离；Michelini等人(2021) [8]开发的INSTANCE数据集是专为机器学习设计的意大利地震数据集，覆盖意大利超过100万个地震事件的数据，成为本研究的主要数据来源。

基于这些数据集，深度学习在地震学领域的应用成果显著。Mousavi 等(2020) [9]提出了CRED系统结合了卷积层和双向长短期记忆单元(LSTM)，能够在噪声环境下检测微震信号；同年，Mousavi 等(2020) [10]提出了一种称为“Earthquake Transformer”的模型，这种基于注意力机制的方法能够同时进行地震检测和震相识别；此外，Mousavi and Beroza (2020) [11]研究开发了能直接预测地震震级、定位及震中距的深度学习模型，这些模型准确地从原始波形预测地震信息，提高了地震定位的精确性，同年Mousavi and Beroza (2020) [12]利用两个贝叶斯神经网络来估计单个观测站地震定位，该方法能够估算震中距离和P波形成时间，并预测地震的起始时间和深度。Nicolis等人(2021) [13]及Ristea and Radoi (2022) [14]也利用深度学习预测地震的强度和位置，他们通过复杂神经网络实现了震中距、深度和震级的高精度估计。Mohamed S等(2023) [15]基于INSTANCE数据集开发了一种结合极端梯度提升(XGB)算法的模型，用于地震烈度的快速预测。

总之，传统研究与深度学习技术相结合，显著提升了地震烈度预测的精度和可靠性，为预防和减轻地震灾害造成的损失提供有力的理论和技术支撑，随着深度学习技术的不断发展和完善，在地震预警尤其是预测烈度方面的技术将会不断提高。

2. 理论基础与模型构建

2.1. 统计理论与机器学习方法

1. 卷积神经网络

卷积神经网络(CNN)作为深度学习中重要的模型架构，广泛应用于处理图像、时间序列等网格结构数据，其特点在于局部连接与参数共享。CNN的基本组件包括卷积层、池化层和全连接层，其中卷积层通过滑动小规模卷积核提取输入数据的局部特征，参数共享机制有效降低了模型复杂性。池化层用于特征降维，通过最大池化或平均池化减少特征维度的同时，保留关键信息。全连接层则将特征整合为全局表示，用于分类或回归任务。

在地震学领域，卷积神经网络因其能够自动提取特征而日益受到关注，本研究采用图1中的一维卷积架构，以高效处理地震波形数据并优化预测精度，一维卷积通过在时间序列数据上滑动卷积核，提取局部时间特征并生成输出序列，相较于二维卷积，其训练更加轻量化且适应长度可变的地震波形数据。

Figure 1. One-dimensional convolution diagram

图1. 一维卷积示意图

2. 循环神经网络

循环神经网络(RNN)是深度学习中用于处理时间序列数据的常用模型。其通过图2中的内部循环连接结构，在网络状态中积累时间步之间的动态信息，从而捕捉长距离的数据依赖关系。然而，传统RNN在处理长时间序列时易受梯度消失或爆炸问题的影响。

Figure 2. RNN structure

图2. RNN结构

为此，长短期记忆网络(LSTM)作为RNN的变体被提出。LSTM单元通过图3中独有的输入门、遗忘门、记忆和输出门结构，有效地捕捉了时间序列中的长期依赖关系，解决了梯度消失的问题，具体的计算流程如公式(1)所示。

$\begin{array}{c}{I}_t=\sigma \left(X_t{W}_{xi}+H_{t-1}{W}_{hi}+b_i\right)\\ F_t=\sigma \left(X_t{W}_{xf}+H_{t-1}{W}_{hf}+b_f\right)\\ O_t=\sigma \left(X_t{W}_{xo}+H_{t-1}{W}_{ho}+b_o\right)\\ {\tilde{C}}_t=\text{tanh}\left(X_t{W}_{xc}+H_{t-1}{W}_{hc}+b_c\right)\\ C_t=F_t\odot C_{t-1}+I_t\odot {\tilde{C}}_t\\ H_t=O_t\odot \text{tanh}\left(C_t\right)\end{array}$ (1)

其中，$X_t$ 代表时间步t下的输入数据，$H_{t-1}$ 代表上一个时间步的隐状态，$C_{t-1}$ 代表上一个时间步的记忆，$W_{xf}$ ，$W_{xi}$ ，$W_{xo}$ ，$W_{xc}$ ，$W_{xf}$ 是输入数据$X_t$ 到各个门或者候选记忆的权重矩阵，$W_{hi}$ ，$W_{hf}$ ，$W_{ho}$ ，$W_{hc}$ 是上一个隐状态$H_{t-1}$ 到各个门或者候选记忆的权重矩阵，$b_i$ ，$b_f$ ，$b_o$ ，$b_c$ 是各个门或者候选记忆的偏置，$\sigma$ 为sigmoid激活函数，$\text{tanh}$ 为双曲正切激活函数，$\odot$ 为哈达玛积。

Figure 3. LSTM unit structure

图3. LSTM单元结构

本研究通过双向LSTM架构，进一步增强了对地震波形数据前后时序特征的建模能力，为烈度预测提供了更准确的特征提取方式。

3. 集成学习

集成学习通过结合多个学习器的预测结果，有效提升模型的稳定性和准确性。其主要方法包括序列化的Boosting和并行化的Bagging。Bagging通过对同一数据集的不同子集进行抽样训练多个模型，并集成其结果，能够有效降低泛化误差。

本研究采用Bagging集成学习框架，将卷积神经网络(CNN)、长短期记忆网络(LSTM)等模型作为基学习器，通过独立训练不同子集上的多个深度学习模型，进一步提高预测的稳定性与精度。最终模型的预测结果通过简单平均法进行集成，使地震烈度预测在多种场景下具备更强的适应性和鲁棒性。

4. 烈度衡量指标

为了与意大利地震数据集相适应，本研究使用基于ShakeMap标准计算的麦加利地震烈度表作为烈度的衡量标准。麦加利地震烈度被美国、韩国、中国香港等地广泛使用，并以“修订麦加利地震烈度表”为基础延伸出由印度、各前苏联加盟国、以色列等地采用的MSK烈度；继而再衍生欧洲采用的EMS烈度，以及中国大陆的烈度标准。

ShakeMap计算地震烈度的标准结合了Wald [16]的研究，综合水平方向最大值的地震动峰值加速度PGA和水平方向最大值的地震动峰值速度PGV，Wald认为当烈度在Ⅴ度以下时，烈度与PGA的相关性更高，当烈度在Ⅶ度以上时，烈度与PGV的相关性则更高，当烈度介于Ⅴ度到Ⅶ度之间时，烈度与PGA和PGV的相关程度差别不大，具体的运算公式如2~5所示。

$MMI=3.66\lg PGA-1.66$ (2)

$MMI=3.47\lg PGV+2.35$ (3)

$MMI=2.20\lg PGA+1.00$ (4)

$MMI=2.10\lg PGV+3.40$ (5)

首先使用公式2进行初步计算，如果$MMI\le V$ (5度)，则使用公式4计算作为最终的烈度，如果$MMI\ge VII$ (7度)，则使用公式3计算$MMI$ ，进一步地，如果$MMI>V$ (5度)，则该值作为最终的烈度值。若公式3计算得到的$MMI\le V$ (5度)，则使用公式4.4计算得到最终的烈度值，如果公式2的初步计算结果介于5度和7度之间，则使用公式2和3或5得到的$MMI$ 值进行加权平均得到最终的烈度值。

2.2. 模型架构

本研究提出的基于神经网络的地震烈度预测模型如图4所示，网络中包括卷积层、循环层以及全连接层。网络输入的是三分量的以物理地震动为单位的数据，因此输入为[3, 12000]的数据，其中第一个维度代表地震波形的东、北、垂直三个通道，第二个维度代表采样率为100 Hz情况下120秒的时间长度。

Figure 4. Earthquake intensity prediction model based on ground motion unit

图4. 基于地震动单位的地震烈度预测模型

首先，输入的地震动数据通过一系列卷积层进行特征提取，卷积层的通道数从3逐渐增加到16，再增至32，以此增强模型的特征捕捉能力。所有卷积层的核大小统一设置为15，步长设为1，且在每次卷积操作中应用填充(padding)以保持特征维度不变。经过各卷积层处理后，三分量地震动数据的尺寸依次变为[16, 3000]、[32, 750]、[64, 187]。维度的缩减主要通过最大池化层完成，每次池化后都会减少波形的空间尺寸。此外，每个卷积层后均应用ReLU函数作为激活函数，以引入非线性处理能力。

处理完卷积层后，数据流向两个双向长短时记忆(LSTM)层，以适应时间序列分析需求。为此，数据被调整为[187,64]的格式。双向LSTM层能够捕获数据的前向和后向时间依赖性，输出的特征维度因此翻倍，达到每个时间步256维的特征输出。在获取到最终时间步的输出特征后，模型通过两个全连接层进一步提取特征。首先是一个含128个神经元的层，采用ReLU激活函数，随后的输出层将特征维度压缩至单一的烈度预测值。

3. 数据集描述与预处理策略

3.1. 数据来源与结构分析

本研究采用了INSTANCE数据集，该数据集由意大利国家地震网络记录，包含约50,000次地震事件和130,000个噪声记录，数据均为三分量波形记录，采样率为100 Hz，提供120秒的时间迹线，是分析地震波形和开展烈度预测的重要数据来源。由于意大利地质构造复杂，位于欧亚板块和非洲板块的交界处，因此地震活动频繁且多样化，涵盖了从浅源到深源地震，这使得意大利成为研究地震行为和改进地震预测技术的理想地。

此外数据集还包括地震迹线的元数据，包括115个参数，部分参数如表1所示，详细记录了关于来源数据、台站元数据、迹线元数据和路径元数据。

Table 1. Parameters included in metadata (part)

表1. 元数据中包含的参数(部分)

Metadata Parameter Name	Description (Chinese)
Source_id	地震和噪声ID (分别为INGV和UTC时间)
source_origin_time	位置首选起源时间(YYYY-MM-DDTHH: MM.SSZ)
source_latitude_deg	位置首选纬度(度)
source_longitude_deg	位置首选经度(度)
source_depth_km	位置首选深度(千米)
station_network_code	两个字符的FDSN网络代码(例如IV)
station_elevation_m	台站海拔(米)
path_ep_distance_km	震中距离
path_hyp_distance_km	震源距离
path_backazimuth_deg	从台站位置到事件震中的方向(度)
path_travel_time_[P,S]_s	P或S波行进时间(秒)
trace_name	HDF5文件中的波形名称
trace_start_time	波形轨迹UTC开始时间(YYYY-MM-DDTHH: MM.SSZ)
trace_[P,S]_arrival_sample	波形轨迹上P和S波起始的样本编号(整数)
trace_pga_cmps^2	最大水平分量的PGA值(厘米/秒²)
trace_pgv_cmps	最大水平分量的PGV值(厘米/秒)
trace_pga_perc	最大水平分量的PGA值(% g)
trace_deconvolved_units	HDF5容量中轨迹的地面运动单位(例如，mps和mps^2分别为米/秒和米/秒^2)

3.2. 数据描述与预处理

对元数据中关于地震样本中的震级、震中距、震源深度和反向方位角进行直方图绘制，结果如图5所示。

为了消除古登堡–里克特定律¹对数据不平衡造成的影响，数据集的作者尝试去汇编一个没有被大量小震级样本扭曲的数据集，最终数据集包括了所有较大震级的迹线，以及逐步减少了小震级样本的迹线数量。可以直观地发现数据中大多数样本的震级均分布在2到3级之间，这说明尽管数据集作者试图平衡不同震级的样本数量，但样本分布不均的现象仍然无法避免。对于震中距离而言，样本分布在0至600公里以上，但绝大多数地震样本均分布在250公里以内。而在震源深度的分布图中，虽然在100到300公里的深度范围内发生了几千次地震，但大多数的地震都属于浅层地震，在400公里以上的深度下，只有几百条或者更少的样本。在反向方位角的分布图中显示，绝大多数p波和s波的样本都位于西北向东南的方向，即沿亚平宁山脉的西北–东南走向。

Figure 5. Histogram of main parameters of earthquake

图5. 地震主要参数直方图

Figure 6. Visualization of randomly selected seismic waveforms

图6. 随机选择的地震波形可视化

通过python中的h5py和matplotlib等库可以对三分量地震波形进行可视化，其中图6(a)~(c)指的是从INSTANCE数据集中包含的宽带HH通道随机挑选的震级介于2到3级的地震波形示例，其中蓝色和红色的垂直线表示p波和s波的到达时间标志。图6(d)~(f)指的是宽带HH通道中随机挑选的有问题的地震波形示例，数据中可能存在单个样本中存在多个地震事件信号或者信噪比过低从而导致信号难以从噪声中区分开来等问题。图6(g)~(i)指的是原始地震波形通过去除仪器响应后得到的EH通道中地震动加速度百分比大于9.3e−4%的地震动数据，图6(m)~(o)则是随机选择的噪声波形示例。

为了筛选出符合训练要求的数据，首先筛选出地震类型为earthquake′且震级类型为里氏震级的事件，随后筛选出在HDF5文件中以米/秒^2为单位的地震动数据，并对信噪比进行过滤，最终保留141,357条有效记录，作为模型训练与测试的输入数据。

4. 模型估计与统计分析

4.1. 模型配置与估计方法

整个模型架构在PyTorch框架下实现，并部署在RTX4090 GPU上以确保计算效率。模型采用均方误差$MSE=\frac{1}{n}{{\displaystyle {\sum }_{i=1}^n\left(y_i-{\widehat{y}}_i\right)}}^2$ 作为损失函数，并使用Adam优化算法，学习率设为0.0001。为了提高模型

的泛化能力并防止过拟合，考虑时间成本与计算成本，研究采用了含有5个基础模型的Bagging集成学习方法，每个基础模型训练时都从训练集中随机抽取样本和特征，采用不放回的抽样方式。此外，在模型中引入了dropout技术，dropout率设置为0.5，以进一步降低过拟合风险。整个训练过程中，实施了早停机制，在连续20个训练周期内，如果验证集上的损失未见明显下降，则自动停止训练。

通过集成学习的方式，模型能够综合各个基础模型的预测结果，通常采用平均法来确定最终的预测结果。这种方法有效地提升了模型对地震烈度的预测精度和鲁棒性，确保了模型在多种地震场景下的适应性和可靠性。

4.2. 统计分析与模型评估

Figure 7. Loss function convergence

图7. 损失函数收敛情况

模型训练的损失函数收敛情况如图7所示。可以看出，训练损失从一个较高的数值快速下降，在第20个迭代后趋于平稳，表明了模型在训练过程中能对训练数据集进行了有效学习。验证损失与训练损失类似，随着迭代次数的增加而逐渐降低，表明模型在没有见过的样本上有着一个不错的泛化能力。

通过训练后的模型，将测试集输入模型得到的结果与真实值绘制散点图如图8所示，可以直观地发现预测值与真实值主要集中在y = x的直线上，表明预测值和真实值之间存在正相关关系，通过计算得到的平均误差为0.09，平均绝对误差为0.238，标准差为0.4，皮尔逊相关系数约为0.881，决定系数为0.76，这表明预测值与真实值之间有较强的正相关性，尽管模型在中高烈度值的预测上表现良好，但在低烈度范围内的预测表现出一定的不稳定性，特别是在烈度为1的情况下，模型会对数据构成一个误判，导致预测烈度出现过高的情况出现。

Figure 8. Scatter plot of true and predicted intensity values

图8. 烈度真实值与预测值的散点图

此外，使用测试集综合震级以及震源深度绘制三维散点图如图9所示，可以发现随着震级的增大，预测烈度通常呈现上升趋势。这表明更大震级的地震往往预测出更高的烈度，这与地震的能量释放量和潜在破坏力成正比。而震源深度虽然对预测烈度有一点影响，但这种影响不如震级显著。

Figure 9. Three dimensional scatter diagram of earthquake magnitude, focal depth and intensity

图9. 震级、震源深度、烈度三维散点图

5. 结论与建议

5.1. 研究总结

本研究结合了深度学习与集成学习策略，对意大利地震数据集(INSTANCE)中的三分量地震动数据进行地震烈度预测。通过融合卷积神经网络(CNN)、长短期记忆网络(LSTM)和全连接网络(FCNN)，本研究有效地捕捉到了时间序列内隐含的空间和时间特征，并通过Bagging算法显著提高了模型的泛化能力和预测精度。实验结果表明，该模型能够准确预测地震烈度，并对不同烈度级别进行了有效区分。此外，通过集成不同子模型的预测结果，进一步提升了预测的稳定性和可靠性。

此外，本研究在INSTANCE数据集上的成果具备迁移应用的可能性，可应用于中国的地震烈度预测。尽管中国和意大利在地质结构和地震类型上存在差异，但两国均位于地震频发地带，经常面临地震灾害的威胁。中国的地质环境复杂且地震活动多样化，特别是在青藏高原和云贵高原等地震活跃区域。利用中国国家地震科学数据中心(CSNCD)发布的完备数据集，运用本研究所提出的方法对中国境内的地震波形及其烈度进行预测分析，不仅能深化对我国地震现象的认识，亦可为地震预警机制建设和灾害防控策略制定提供关键的科学依据和技术支持。

5.2. 研究限制与未来方向

尽管本研究取得了一定的成就，但仍然面临若干局限。首先，地震数据固有的复杂性和不确定性制约了模型预测精度的进一步提升。其次，即便已采取多种措施以减轻过拟合现象，该模型在应对低强度地震事件时仍显示出一定程度的波动性。另外，虽然所使用的数据集具有广泛的覆盖面，但由于缺乏高强度地震案例，导致数据分布极不平衡，从而削弱了模型的泛化性能。

在后续的工作中可以通过融合来自不同地理位置及各类地震事件的数据集来提升模型的泛化能力和适应性。此外，探索诸如Transformer等前沿的深度学习架构与算法的应用，为地震数据分析带来全新的视角，并提高预测精度。

NOTES

¹古登堡–里克特定律（Gutenberg-Richter law）是地震学中的一个定律，表示震级与某一地区大于等于该震级数量之间的关系，由地震学家宾诺·古登堡与查尔斯·弗朗西斯·里克特与1956年提出。该定律的表达式为：${\log }_{10}N=a-bM$ ，其中N表示震级≥ M的地震数量，a与b为常数。