Probabilistic aftershock hazard assessment for Jiuzhaigou MS7.0 earthquake in 2017
-
摘要: 本文介绍了余震危险性概率分析的概念，完整列出其相关公式的解析表达式并阐述了余震危险性概率分析的具体方法。以2017年九寨沟MS7.0地震为例，首先，求得该地震余震序列的相关参数，结果显示：九寨沟MS7.0地震余震序列理论最大余震震级约为ML5.3；b值约为0.784 1，明显低于中国西南地区同类型的其它地震，表明九寨沟地震余震区应力水平相对较高；p值约为1.109 7，明显高于中国西南地区同类型的其它地震，表明九寨沟地震余震序列随时间衰减较快。其次，计算此次地震余震与主震释放能量的比例关系，结果显示：在九寨沟MS7.0地震事件中，截止到10月22日约99.69%的能量为主震所释放，0.31%的能量为余震所释放。最后，利用九寨沟MS7.0余震序列参数结果并结合衰减关系，分别计算了主震后0—1天、1—10天、10—30天和90—100天内不同断层距内不同水平的峰值加速度和峰值速度的超越概率，结果显示：随着断层距增加，在相同时间间隔内超越概率的值呈下降趋势；在震后第一天内余震危险性最高，随着震后时间增加，相应的超越概率值呈现出明显下降趋势，表明九寨沟MS7.0地震余震危险性主要来自早期余震。本文的相关成果可以为震后地震危险性分析提供参考，并为短时间内应急救援及灾后重建提供辅助决策意见。Abstract: Probabilistic aftershock hazard assessment is important to assess the hazard of aftershocks. In this paper, we introduce the concept of probabilistic aftershock hazard assessment, and give the formulas related to the method of it.In response to the Jiuzhaigou MS7.0 event in Sichuan Province in 2017, firstly, we obtained the related parameters of the aftershock. The results show that the maximum aftershock magnitude of the given aftershock sequence is about ML5.3. The b value of Jiuzhaigou event is about 0.784 1, which is lower than the same type of earthquake in southwest China, indicating that the stress level of the aftershock area is relatively high. The p value is about 1.109 7, obviously higher than events of the same type, which shows that this aftershock sequence decays faster with time. Secondly, according to the parameters estimation, 99.69% of the energy was released by the mainshock and the rest by the aftershocks. Finally, combining the aftershock sequence parameter and the attenuation relationship, we calculated the probability of exceedance certain level of peak ground acceleration and peak ground velocity in 0−1 day, 1−10 days, 10−30 days and 90−100 days intervals within different fault distances after the mainshock. The results show: ① The value of exceedance probability decreases with fault distance in the same time interval; ② The risk of aftershocks is highest in the first day after the mainshock; the exceedance probability value exhibits a downward trend with the increase of time, indicating that the aftershock hazard of the JiuzhaigouMS7.0 event mainly comes from early aftershocks. These results provide reference for aftershock risk analysis, and auxiliary decision-making opinions for short time emergency rescue and post disaster reconstruction.
-
引言
我国属于地震频发国家，如果能够有效预测大地震的时间、位置和震级，便可及时采取应对措施，保护人民的生命财产安全，降低房屋、桥梁和道路等基础设施的破坏程度和经济损失，从而维护社会稳定和发展（陈运泰，2008）。然而，由于地震发生概率低，孕震机制复杂，地球深部观测技术受限，地震预测仍待长期持续的研究（陈运泰，2007）。
随着观测技术的不断进步，大量震例研究表明，大地震的孕育过程通常伴有地表异常，如2008年汶川MS8.0地震（江在森等，2009；Zhu et al，2019）、2011年日本东北MW9.0地震（Hirose，2011）、2017年九寨沟MS7.0地震（高曙德，2020）等。GPS、应力和应变观测资料作为地壳形变指标，在多次大地震前显示出断层闭锁或变幅异常（张希等，2019；顾国华，王武星，2020；牛安福等，2021；Yu et al，2021a，b）。近年来，国内外学者研究发现，地震孕育过程中，岩石破裂还会产生压电效应、动电效应及压磁效应等一系列地震电磁效应（De Santis et al，2017；Parrot，2017；Pulinets et al，2018）。以芦山地震为例，姑咱台钻孔应变观测在主震前半年就记录到脉冲异常（池顺良等，2013），除此之外，附近台站还观测到视电阻率、地电场异常（高曙德，2016）、地磁场异常（廖晓峰等，2018）等。由于电磁波的趋肤深度与震源深度相当，低频电磁辐射也被认为是非常有力的地震前兆（Ouyang et al，2018）。
应变数据可反映地下应力变化，对于地震研究具有重要意义。在青海玉树、日本东北和四川芦山三次大地震中，青海地区的6个钻孔应变台站均记录到了与震级和震源方位相符的应变阶跃（张敏等，2014）；针对汶川地震，邱泽华等（2009，2012）在剔除钻孔应变数据的周期信号后，采用超限率分析的方法提取了高频信号中的幅度异常，发现异常在震前开始增加、震后逐步减少，认为可能与汶川地震有关；对于芦山地震，邱泽华等（2015）通过分析姑咱台在地震数天前记录到的钻孔应变波形数据，排除了远震干扰和实地施工等人为干扰，认为剩余的幅度异常变化与芦山地震呈因果关系；池顺良等（2013）同样通过观察姑咱台的钻孔应变波形数据的幅度变化，以及计算观察到的应变异常方向与地震应力的主方向一致，证明了芦山地震前几天和震前4至6个月记录的异常现象可能反映了地震发震前地层的蠕滑；刘琦等（2014）利用S变换方法也对芦山地震前后的钻孔应变数据进行了时频分析，发现两簇高能量异常表现出与芦山地震较好的相关性；Chi等（2019）采用变分模态分解方法对钻孔应变观测数据进行分解，通过状态空间模型，采用频率对比方式去掉气压和钻孔水位的影响后，发现了芦山地震的震前异常；Zhu等（2019）提出了钻孔应变数据的震前负熵异常提取方法，汶川和芦山地震的负熵结果分析显示，各台站负熵异常的累计值在震前均出现以幂律形式的加速增长；之后，Yu等（2020）又采用多通道奇异谱分析联合分解了多台站钻孔应变数据，构造了多台站的空间拓扑网络，联合提取了川滇地区13次地震前的应变异常。这些研究都表明，大地震前钻孔应变数据能够反映地震孕育过程。
然而，对于钻孔应变数据的震前异常提取研究，目前还停留在震例分析与小样本统计分析阶段。Yu等（2020）分析了2010—2017年多台站网络度的异常，也只是研究了29次能量较大的地震，对于M3.0—5.0中等地震的震前异常目前尚未涉及。另外，大多数研究只关注了应变数据异常与地震时间上的关联，对异常的空间分布以及震源方位的分析研究甚少。
本文拟依据大量钻孔应变观测数据，发挥数据驱动方法的优势，采用卷积神经网络框架，设计地震震级与方位的短临预测模型；根据四分量钻孔应变仪原理，计算分量应变、面应变、剪应变间的相关系数和自洽系数，构造应变特征数据集以作为模型的输入；模型输出分别为震级以及方位，按震级大小将震级划分为三类，根据钻孔应变仪方位角将方位划分为五类；通过混淆矩阵计算准确率、召回率以及F1分数评价与修正模型。最后，通过对我国西南地区的姑咱、永胜、昭通和腾冲台的钻孔应变数据进行训练与测试，探索应变数据与地震发生大小、方位的内在联系。
1. 数据简介
1.1 四分量钻孔应变观测原理
四分量钻孔应变仪一般放置于地下数十米断裂带周围的基岩中，通过内部传感器的双环模型测量孔径变化量。孔径变化量与原孔径的比值为Sθ，表示沿孔径变化方向θ即传感器方向直接观测到的孔径变化。根据该模型，在相对远处有均匀水平应变ε1和ε2时，有（李进武，邱泽华，2014）
$$ S_{\theta} ＝ A （ \varepsilon_{ 1}－\varepsilon_{ 2} ） ＋B （ \varepsilon_{ 1}＋\varepsilon _{ 2} ） \mathrm{cos} ［ 2 （ \theta －\phi ） ］ \text{，} \text{ } $$ (1) 式中：$\phi $为1分量的方位角；A，B为耦合系数，其大小受周围水泥和套筒的尺寸及材质的共同影响。
四分量钻孔应变仪的四个传感器按相邻夹角45°布置，根据式（1），有四分量观测值：
$$\left\{ {\begin{array}{*{20}{l}} {S_{1} ＝ S_{\theta} ＝ A （ \varepsilon_{1} － \varepsilon _{2} ） ＋ B （ \varepsilon_{1} ＋ \varepsilon _{2} ） \cos ［ 2 （ \theta － \phi ） } ］ \text{，} \\ {S_{2} ＝ S_{\theta ＋ {\pi /4 }} ＝ A （ \varepsilon_{1} － \varepsilon _{2} ） － B （ \varepsilon_{1} ＋ \varepsilon _{2} ） {\rm{sin}} ［ 2 （ \theta － \phi ） } ］ \text{，} \\ {S_{3} ＝ S_{\theta ＋ {\pi / 2} }＝ A （ \varepsilon_{1} － \varepsilon _{2} ） － B （ \varepsilon_{1} ＋ \varepsilon _{2} ） \cos ［ 2 （ \theta － \phi ） } ］ \text{，} \\ {S_{4} ＝ S_{\theta ＋ {3\pi / 4}} ＝ A （ \varepsilon_{1} － \varepsilon _{2} ） ＋ B （ \varepsilon_{1} ＋ \varepsilon _{2} ） {\rm{sin}} ［ 2 （ \theta － \phi ） ］ } \text{，} \end{array}} \right.{\rm{ }} $$ (2) 在式（2）的基础上可作变换，令$ S_{13} ＝ S_{1} － S_{3} $和$ S_{24} ＝ S_{2} － S_{4} $，S13和S24称为剪应变。
四分量钻孔应变观测还有一个重要的特性就是观测自洽，它可以判断元件的性能，也可评价观测数据的准确性（于紫凝，2022）。自洽方程为
$$ {S}_{1}＋{S}_{3}＝k （ {S}_{2}+{S}_{4} ） \text{ } \text{，} $$ (3) 式中，k为自洽系数，若k≥0.95，即认为数据满足自洽。实践证明，若应变仪探头布设于完整岩石中时，其观测曲线应大体上符合观测自洽。自洽系数k也常作为震前数据异常的重要依据，在汶川、鲁甸和康定地震前，姑咱台和昭通台均观测到了失恰现象（池顺良等，2014）。
1.2 钻孔应变特征数据集设计
四分量钻孔应变仪每隔45° 布置一路应变传感器，四路应变分别为$ S_{1} $，$ S_{2} $，$ S_{3} $和$ S_{4} $，采样间隔为1 min。为了挖掘到更多的应变信息以及多路应变之间的关联信息，本文又构造剪应变$ S_{1} － S_{3} $，$ S_{2} － S_{4} $及自洽方程两端面应变$ S_{1} ＋ S_{3} $和$S_{2} ＋ S_{4} $，共得到8维应变数据。池顺良等（2014）多次对面应变、剪应变之间的相关性以及四路应变的自洽系数进行分析，发现汶川、鲁甸、康定地震前，钻孔应变可能反映了地震成核的过程。所以，本文将这些应变间的皮尔逊相关系数与自洽系数作为特征数据集，共24维，其物理含义如表1所示。
表 1 24维应变数据特征Table 1. Features of 24-dimensional strain data序号 特征名称 物理含义 1 corr01_02 S1与S2间的皮尔逊相关系数 2 corr01_03 S1与S3间的皮尔逊相关系数 3 corr01_04 S1与S4间的皮尔逊相关系数 4 corr02_03 S2与S3间的皮尔逊相关系数 5 corr02_04 S2与S4间的皮尔逊相关系数 6 corr03_04 S3与S4间的皮尔逊相关系数 7—10 corr01_13—corr04_13 S1−S4与剪应变S1−S3间的皮尔逊相关系数 11—15 corr01_24—corr04_24 S1−S4与剪应变S2−S4间的皮尔逊相关系数 16—19 corr01+13—corr04+13 S1−S4与面应变S1＋S3间的皮尔逊相关系数 20—23 corr01+24—corr04+24 S1−S4与面应变S2＋S4间的皮尔逊相关系数 24 corr13+24 自洽系数 注：S1，S2，S3和S4为钻孔应变四分量观测值。 2. 基于钻孔应变的一维卷积神经网络地震活动性预测模型
2.1 一维卷积神经网络模型原理
为防止模型参数过多导致训练过拟合，本文采用一维卷积神经网络（convolutional neural network，缩写为CNN），提取震前应变数据特征训练模型，实现地震震级及方位的预测。卷积神经网络利用其卷积层参数共享机制的优势，可以有效捕捉到数据中的局部特征（卢宏涛，张秦川，2016）。如图1所示，模型由输入层、隐藏层、输出层三部分组成。模型输入为单输入，输入层24个神经元节点x1， x2， ···， x24分别对应24种特征。模型输出为多输出多分类：多输出是指模型输出具有多个预测值$ f{ （ x_{1}\text{，} x_{2}\text{，} \cdots \text{，} x_{24} ） } ＝ （ y_{1}\text{，} y_{2} ） $，输出层1和输出层2分别对应震级标签$ y_{1} $和方位标签$y_{2} $两个输出预测值；多分类是指每个预测值可以取多个类别，$ y_{1} \in \{ {y_{11}} {\text{，}} {y_{12}} \text{，} {y_{13}} \} $分别对应无震、3.0≤MS<5.0及MS≥5.0三个类别，$y_{2} \in \{ {y_{21}}\text{，} y_{22}\text{，} y_{23}\text{，} y_{24}\text{，} y_{25} \} $分别对应无地震及划分地震方位的四个方向，共五个类别。隐藏层1和隐藏层2分别为卷积区和全连接区。
卷积区由多个卷积层、批量标准化层、最大池化层组成，其结构如图2所示。输入为335×24的应变特征矩阵，第一个卷积层使用64个卷积核提取出329×64的特征图。特征图在批量正则化层进行归一化处理，并输入最大池化层进行降采样，得到输出为164×64的特征图，并将其作为第二个卷积层输入，进行第二轮卷积、归一化、池化操作。两个卷积层采用不同大小的卷积核，具有参数共享和局部连接的特性，可以减少模型的参数量，增强模型的表征能力（周飞燕等，2017）。在两轮卷积、归一化、池化操作后得到51×64的特征图，并将其输入全连接区。
全连接区由全局平均池化层、全连接层和丢弃层组成（图3）。卷积区输入51×64的特征图，通过全局平均池化层对每个通道特征图进行平均运算，得到该通道特征图的平均强度，64个通道特征图取得了长度为64的特征向量并输入全连接层。两个全连接层分别采用32个和16个神经元，用于提取特征并进行分类。每个全连接层后添加丢弃层，防止模型过拟合。
输出层使用Softmax激活函数，可以将输出层1和输出层2的特征向量映射为（0， 1）区间内的预测值，表示震级和方位的类别概率分布，并将概率最大类别作为震级和方位类别的预测结果。
整体网络设置中，ReLU激活函数应用于每个卷积层和全连接层，并在每个卷积层和全连接层采用L2正则化方法。L2正则化通过在损失函数中添加与参数平方和成比例的惩罚项来约束模型复杂度，在卷积层采用$ \lambda ＝ 0.01 $的权重，在全连接层采用$ \lambda ＝ 0.001 $的权重。最小化损失函数和L2正则化项的加权和可以降低模型的复杂性，以避免过拟合。计算方法为
$$ J （ w ） ＝\underset{w}{\mathrm{min}} ［ L （ w ） ＋ \lambda {\Vert w\Vert }_{2}^{2} ］ \text{ } \text{，} $$ (4) 式中，J为加入正则化项的网络代价函数，w为神经网络权重参数，L为网络的损失函数。
2.2 模型编译与参数优化
模型编译使用多元交叉熵损失函数和Adam优化器，并采用学习率衰减策略动态调整学习率。学习率衰减策略根据验证集损失值判断模型训练状态，当验证集损失值在一定轮次（10轮次）内未下降时，学习率则乘以衰减因子（0.1）以使学习率缩小为原来的十分之一，直至学习率达到最小值（1×10−8），从而避免由学习率过小或过大导致的收敛缓慢、梯度爆炸等问题，并可提高模型的收敛速度和准确性。
为了优化模型超参数并提高模型性能，本文使用随机搜索算法寻找最优超参数，构造包含卷积核大小及数量、池化窗口、全连接层神经元个数、L2正则化系数、丢弃层丢弃率、样本批次等参数的参数空间。在参数空间随机选择一组超参数组合，计算此时模型的验证集准确率，并向着模型最优的方向不断动态调整超参数。依次持续迭代，直到模型的性能达到稳定且最优水平为止。选择该组超参数作为模型超参数。
$${ {{x}}^{*}}＝\underset{{\boldsymbol{x}}\in X}{{\mathrm{argmax}}}f （ x ） \text{ }\text{，} $$ (5) 式中，x为超参数集合，$ X $为超参数空间，$ f （ x ） $为模型验证集的性能指标函数，x*为最优的超参数向量。本文构建的1D_CNN模型共两层，超参数数量适中，随机搜索过程基本可以覆盖全部参数空间。
2.3 模型预测结果评价指标
混淆矩阵是用于评估分类模型性能的一种矩阵（表2）。根据模型的预测结果与真实标签的一致性进行统计得到真正例TP、假正例FP、真反例TN和假反例FN四个指标。在混淆矩阵的基础上可以得到准确率Pacc、精确率Ppcs、召回率Prec等常用的模型评估指标。各指标的意义及计算方法为：准确率，预测正确样本数占总样本数的比例，
表 2 混淆矩阵Table 2. Confusion matrix实际正例 实际负例 预测正例 真正例（TP） 假正例（FP） 预测负例 假反例（FN） 真反例（TN） $$ P_{{\mathrm{acc}}}＝\frac{TP＋TN}{TP＋TN＋FP＋FN} \text{；} $$ (6) 精确率，预测为正例中正确预测的样本比例，
$$ P_{{\mathrm{pcs}}}＝\frac{TP}{TP＋FP}\text{；} $$ (7) 召回率，正例中被正确预测的样本比例，
$$ P_{{\mathrm{rec}}}＝\frac{TP}{TP＋FN}\text{；} $$ (8) F1分数（F1-score），精确率和召回率的调和平均数，综合评价分类器的性能，
$$ F_{1}\text＝\frac{{2} P_{{\mathrm{pcs}}} P_{{\mathrm{rec}}}}{P_{{\mathrm{pcs}}} ＋P_{{\mathrm{rec}}} } .$$ (9) 3. 基于钻孔应变数据的地震活动性预测结果及分析
3.1 地震震级与方位数据集的构造
本文选择我国西部永胜、昭通、姑咱和腾冲四个台站2010—2017年的应变数据构造震前应变特征数据集。根据震级大小将地震分为无震、3.0≤MS<5.0及MS≥5.0三类，再根据震源相对于钻孔应变仪的方向，将方位分为五类，制作样本标签。以永胜台为例，台站及地震筛选分布以及震源预测的方位示意图如图4所示。震源方位是以钻孔应变仪1分量为基准，顺时针每隔45°划分为一类，与无震时的样本一起将震源方位划分为五类。
![]()
图 4 钻孔应变台站位置和地震分布右下角为震源方位划分示意图Figure 4. Location of borehole strain stations and earthquakes distribution mapThe inset at the lower right is the schematic diagram of epicenter orientation delineation以四个台站为中心，筛选地震样本基本情况（表3），包括样本数量、最大震级、最小距离、最大距离等。可以看到M5.0以上地震较少，大约是M3.0—5.0地震的1/7，而大地震的影响范围大，孕震时间又长。为了增加大地震的样本数量，对于M5.0以上地震的选择，本文首先将空间半径扩大到500 km。M3.0—5.0地震的选择范围则是距离中心台站200 km以内。其次，对于M5.0以上地震，本文使用窗长7天、步长1天的滑动时间窗口来预测之后7天的地震大小和方位。也就是说，大地震作为预测窗口的第一天和第七天，可构造出7个样本。这样既增加了样本的数量，同时也可保证样本的真实性和有效性。另外，由于无震时的样本数目众多，本文从中均随机抽取300个作为该类别的应变特征数据集。
表 3 永胜、昭通、姑咱和腾冲四个台站筛选的地震基本情况Table 3. Basic information about the earthquake samples from the stations Yongsheng，Zhaotong，Guzan and Tengchong台站 MS3.0—5.0地震数量 MS≥5.0地震数量 最大震级MS 最小距离/km 最大距离/km 永胜 99 47 7.0 5.96 491.86 昭通 138 31 7.0 31.09 497.47 腾冲 172 52 7.6 6.44 498.75 姑咱 287 36 7.0 13.89 431.65 3.2 模型训练与预测结果分析
本文以震前一周的钻孔应变特征预测之后一周内的地震震级和方位。以永胜台为例，图5给出了震级预测和方位预测的模型损失率曲线。模型训练集和验证集的损失率曲线均在150轮训练后平稳收敛。昭通、腾冲和姑咱三个台站的模型训练结果与永胜台均一致，四个台站验证集损失率则稳定在0.1—0.3之间。
![]()
图 5 永胜台预测模型的损失率曲线（a） 震级预测损失曲线；（b） 方位预测损失曲线Figure 5. Loss curves of the prediction model at Yongsheng station（a） Loss curves for the magnitude prediction；（b） Loss curves for the orientation prediction模型训练完成后，本文在测试数据集上分析模型的预测效果。永胜、昭通、姑咱和腾冲四个台站的数据集测试结果列于表4和表5。总的来说，在震级和方位的预测上，模型均表现出了较高的准确率。尤其对于0类标签预测，即非地震标签，有着非常高的准确率和召回率，F1值均接近或达到0.85，这说明正常情况下钻孔应变数据特征是比较容易捕捉的。实际上，钻孔应变仪在没有外界干扰时，其波形数据中记录着的光滑且有周期性的固体潮是数据中最显著的特征（邱泽华等，2015）。
表 4 永胜、昭通、姑咱和腾冲台震级预测结果表Table 4. Magnitude prediction results for the stations Yongsheng，Zhaotong，Guzan and Tengchong标签 精确率 召回率 准确率 F1值 样本数 标签 精确率 召回率 准确率 F1值 样本数 永胜台 0 0.903 0.823 0.861 147 昭通台 0 0.813 0.875 0.843 144 1 0.829 0.907 0.866 107 1 0.836 0.825 0.830 154 2 0.720 0.766 0.742 47 2 0.810 0.567 0.677 30 0.844 301（总） 0.823 328（总） 腾冲台 0 0.837 0.860 0.848 143 姑咱台 0 0.788 0.748 0.768 139 1 0.888 0.864 0.876 184 1 0.868 0.909 0.888 318 2 0.880 0.772 0.822 57 2 0.840 0.636 0.724 33 0.852 384（总） 0.845 490（总） 表 5 永胜、昭通、姑咱和腾冲台震源方位预测结果表Table 5. Orientation prediction results for the stations Yongsheng，Zhaotong，Guzan and Tengchong标签 精确率 召回率 准确率 F1值 样本数 标签 精确率 召回率 准确率 F1值 样本数 永胜台 0 0.879 0.837 0.857 147 昭通台 0 0.801 0.868 0.833 144 1 0.857 0.686 0.762 35 1 0.879 0.886 0.883 140 2 0.754 0.897 0.819 58 2 1.000 0.467 0.636 15 3 0.781 0.862 0.820 58 3 0.882 0.682 0.769 22 4 0.000 0.000 0.000 3 4 0.571 0.571 0.571 7 0.827 301（总） 0.838 328（总） 腾冲台 0 0.810 0.867 0.838 143 姑咱台 0 0.728 0.770 0.748 139 1 0.750 0.632 0.686 19 1 0.840 0.750 0.792 28 2 0.743 0.728 0.735 103 2 0.593 0.574 0.583 61 3 0.824 0.792 0.808 106 3 0.775 0.811 0.793 106 4 0.833 0.769 0.800 13 4 0.818 0.776 0.796 156 0.794 384（总） 0.755 490（总） 对永胜台来说，在震级预测上，模型对1类样本（3.0≤MS<5.0）的识别有较高的准确率、召回率和F1值，均可达到0.82以上。2类样本（MS≥5.0）的F1值达到0.74以上，结果也处于相对高的水平，但较低于1类样本。可能原因有二：一是根据古登堡—里克特定律，随着震级的增大，地震数量呈指数减少，2类样本也就是MS≥5.0地震的样本数量较少，模型没有足够的样本特征去学习，从而导致了一些误判；二是大地震的发生本身就具有一定的复杂性，地震类型、传播途径等因素都会导致震前特征复杂多变，这也增加了模型的预测难度。
对于震源方位预测，与震级预测类似，0类预测的评价指标都很高，可达到0.87以上。除此之外，其它类别预测结果的准确率、召回率和F1指标都与样本数量有直接关系，如永胜台4类结果，该类样本数很少，因此模型未得到充分训练，导致模型评价指标相对较低。
昭通、姑咱和腾冲三个台站的预测结果也与永胜台类似，其震级预测的准确率分别可达到0.823，0.845和0.852，震源方位预测的准确率分别可达0.838，0.755和0.794。这说明，三个台站钻孔应变数据的震前特征也都能够由本文提出的1D-CNN模型挖掘得到，模型在地震发生前一周的应变数据中学习到了十分有效的震前特征。
针对我国五次典型的大地震，本文进一步单独给出了实际预测结果（表6）。从表6可以看出，5次地震震级和方位的预测标签均与实际相符，说明本文模型捕捉到了有效的震前应变信息，有一定的预测能力。
表 6 五次典型强震的实际预测效果Table 6. Prediction results of five typical strong earthquakes序号 发震日期 发震地点 MS 预测震级标签 实际震级标签 预测方位标签 实际方位标签 1 2 010−04−14 青海玉树 7.3 2 2 1 1 2 2 013−07−22 甘肃岷县 6.7 2 2 2 2 3 2 014−08−03 云南鲁甸 6.6 2 2 3 3 4 2 016−10−17 青海杂多 6.3 2 2 1 1 5 2 017−08−08 四川九寨沟 7.0 2 2 2 2 此外，为了观察不同震级类别所对应的最优参数范围，本文还分析了不同震级类别下，不同时间窗口及预测窗口对模型预测效率的影响。由于四个台站的结果相差不大，本文依旧以永胜台站为例，结果如图6所示。
![]()
图 6 各震级标签下不同时间窗口及预测窗口下的预测精确率（a） 0类标签；（b） 1类标签；（c） 2类标签Figure 6. Prediction accuracies under different time windows and prediction windows for each class（a） Class 0；（b） Class 1；（c） Class 2整体上，三个震级的预测精确率都随着时间窗口的增大而上升。这是由于时间窗口越长，含有的应变特征信息就越多，准确率自然就会上升。首先，0类预测基本上在时间窗口为5、预测窗口为3后，精确率就基本上稳定在0.80以上。其次，1类和2类的准确率规律性略差，这可能是因为震前特征丰富，增加了随机性。其中，1类是在时间窗口为7、预测窗口为3后，精确率可达到0.80以上。其余方案准确率都很低。而2类整体上比1类高一些，随机性也更大。这可能反映了钻孔应变数据对大地震具有一定的短临预测能力，同时也说明了不同大地震间的震前特征存在差异，可能有些特征显著、有些则难以提取，有些异常特征持续长、有些则短。因此，如何设计网络提取样本量小、差异化又大的大地震震前异常特征，也需要试验和讨论。最后，由于本文使用了滑动窗口来构造样本，因此在两个连续样本中可能会存在部分重复特征的情况，而时间窗口越长，重复就越多，因此图6右下角的精确率可能会比实际情况偏高一些。
4. 讨论与结论
面对地震短临预测这一科学难题，本文突破传统方法的局限，提出了采用深度学习手段进行地震震级和方位预测的策略。通过钻孔应变仪四分量数据构建震前应变特征数据集，采用1D-CNN搭建了地震震级和方位的预测模型，其中，地震震级划分为无震、3.0≤MS<5.0以及MS≥5.0三类，震源方位标签根据四分量钻孔应变仪1分量方位角方向划分为四类。结果显示，震前应变特征在震级和方位预测上的预测准确性可高达80%以上。这表明，对于不同方向放置的钻孔应变仪，本文搭建的1D-CNN模型都能利用该台站数据预测地震震级与大致方位，也说明了卷积神经网络结构有能力从钻孔应变数据中挖掘到有效的震前异常特征。这为地震前兆观测研究提供了新的思路和方法。
然而，本文将复杂的地震预测问题抽象成深度学习中的分类任务，这与实际地震预测业务还有一定差距。实际预测中，从业人员要根据全新的观测数据给出地震预测意见，而三要素的预测报准率、虚报率、漏报率都会在决策中起到重要影响。因此，地震预测仍是一个科学难题。若要使预测任务更加贴近实际应用，本文仍有一些待改进之处。例如：由于震级及空间位置的双方面限制会导致大地震的样本非常少，本文将震级与方位视为两个独立的标签，想要实现震级与方位的同时预测还存在一定挑战；本文网络模型在四川云南地区的地震预测中表现良好，然而在其它地区是否具有普适性还有待验证。考虑到实际应用价值，未来将设计更完善的网络，将两种标签组合起来，解决样本不平衡问题，实现地震三要素预测；或采用多个应变台站联合，加入震中距、角度等更多的地震位置信息作为预测标签，通过划分分辨率较高的空间网格来提高方位预测的精度；对弱震区或其它多震区的应变观测台站进行试验验证，构建稳健、适用强的短临地震预测模型。
-
![]()
图 1 九寨沟MS7.0地震震中位置以及截至2017年10月22日的余震序列分布
Figure 1. Location of Jiuzhaigou MS7.0 earthquake and its aftershocks （ML≥0） distribution from the occurrence of the mainshock to October 22，2017
![]()
图 2 九寨沟MS7.0地震余震震级-频度图
Figure 2. Magnitude-frequency distribution of aftershocks of Jiuzhaigou MS7.0 earthquake
![]()
图 3 九寨沟地震余震累积数目随时间的变化曲线
Figure 3. Curve changes of cumulative aftershocks of Jiuzhaigou earthquake with time
![]()
图 5 不同时段内九寨沟MS7.0地震余震PGV超越概率变化示意图
Figure 5. Probabilities of exceeding PGV for aftershocks of Jiuzhaigou MS7.0 event during different periods
-
蒋长胜,庄建仓,吴忠良,毕金孟. 2017. 两种短期概率预测模型在2017年九寨沟7.0级地震中的应用和比较研究[J]. 地球物理学报,60(10):4132–4144 doi: 10.6038/cjg20171038 Jiang C S,Zhuang J C,Wu Z L,Bi J M. 2017. Application and comparison of two short-term probabilistic forecasting models for the 2017 Jiuzhaigou,Sichuan,MS7.0 earthquake[J]. Chinese Journal of Geophysics,60(10):4132–4144 (in Chinese) doi: 10.6038/cjg20171038
蒋海昆,郑建常,吴琼,曲延军,李永莉. 2007. 传染型余震序列模型震后早期参数特征及其地震学意义[J]. 地球物理学报,50(6):1778–1786 Jiang H K,Zheng J C,Wu Q,Qu Y J,Li Y L. 2007. Earlier statistical features of ETAS model parameters and their seismological meanings[J]. Chinese Journal of Geophysics,50(6):1778–1786 (in Chinese)
蒋海昆. 2010. 5·12汶川8.0级地震序列震后早期趋势判定及有关问题讨论[J]. 地球物理学进展,25(5):1528–1538 Jiang H K. 2010. Review of tendency judgement of the 5·12 Wenchuan M8 earthquake and discussion on some problems[J]. Progress in Geophysics,25(5):1528–1538 (in Chinese)
苏有锦,赵小艳. 2008. 全球8级地震序列特征研究[J]. 地震研究,31(4):308–316 Su Y J,Zhao X Y. 2008. Characteristics of global earthquake sequences with MW≥8.0[J]. Journal of Seismological Research,31(4):308–316 (in Chinese)
中国地震台网中心. 2017. 全国统一快报目录[EB/OL]. [2017-10-22]. http://www.csi.ac.cn/publish/main/813/5/index.html. China Earthquake Networks Center. 2017. National bulletin[EB/OL]. [2017-10-22].http://www.csi.ac.cn/publish/main/813/5/index.html （in Chinese）.
Båth M. 1965. Lateral inhomogeneities of the upper mantle[J]. Tectonophysics,2(6):483–514 doi: 10.1016/0040-1951(65)90003-X
Boore D M,Atkinson G M. 2008. Ground-motion prediction equations for the average horizontal component of PGA,PGV,and 5%-damped PSA at spectral periods between 0.01 s and 10.0 s[J]. Earthquake Spectra,24(1):99–138 doi: 10.1193/1.2830434
Campbell K W,Bozorgnia Y. 2008. NGA ground-motion model for the geometric mean horizontal component of PGA,PGV,PGD and 5% damped linear elastic response spectra for periods ranging from 0.01 to 10 s[J]. Earthquake Spectra,24(1):139–171 doi: 10.1193/1.2857546
Chiou B J,Youngs R R. 2008. An NGA model for the average horizontal component of peak ground motion and response spectra[J]. Earthquake Spectra,24(1):173–215
Console R,Lombardi A M,Murru M,Rhoades D. 2003. Båth’s law and the self-similarity of earthquakes[J]. J Geophys Res,108(B2):2128
Cornell C A. 1968. Engineering seismic risk analysis[J]. Bull Seismol Soc Am,58(5):1583–1606
Gallovič F,Brokešová J. 2008. Probabilistic aftershock hazard assessment I:Numerical testing of methodological features[J]. J Seismol,12(1):53–64 doi: 10.1007/s10950-007-9072-0
Gutenberg B,Richter C F. 1942. Earthquake magnitude,intensity,energy and acceleration[J]. Bull Seismol Soc Am,32(3):163–191
Hainzl S,Marsan D. 2008. Dependence of the Omori-Utsu law parameters on main shock magnitude:Observations and modeling[J]. J Geophys Res,113(B10):B10309 doi: 10.1029/2007JB005492
Hamdache M,Peláez J A,Kijko A,Smit A. 2017. Energetic and spatial characterization of seismicity in the Algeria-Morocco region[J]. Nat Hazards,86(S2):273–293 doi: 10.1007/s11069-016-2514-7
Helmstetter A,Sornette D. 2003. Båth’s law derived from the Gutenberg-Richter law and from aftershock properties[J]. Geophys Res Lett,30(20):2069
Kisslinger C,Jones L M. 1991. Properties of aftershocksequences in southern California[J]. J Geophys Res,96(B7):11947–11958 doi: 10.1029/91JB01200
Omori F. 1894. On after-shocks of earthquakes[J]. J Coll Sci Imp Univ Tokyo,7:111–200
Rodriguez-Marek A,Montalva G A,Cotton F,Bonilla F. 2011. Analysis of single-station standard deviation using the KiK-net data[J]. Bull Seismol Soc Am,101(3):1242–1258 doi: 10.1785/0120100252
Rodríguez-Pérez Q,Zúñiga F R. 2016. Båth’s law and its relation to the tectonic environment:A case study for earthquakes in Mexico[J]. Tectonophysics,687:66–77 doi: 10.1016/j.tecto.2016.09.007
Shcherbakov R,Turcotte D L. 2004. A modified form of Båth’s law[J]. Bull Seismol Soc Am,94(5):1968–1975 doi: 10.1785/012003162
Shcherbakov R,Turcotte D L,Rundle J B. 2004. Ageneralized Omori’s law for earthquake after shockdecay[J]. Geophys Res Lett,31(11):L11613 doi: 10.1029/2004GL019808
Shcherbakov R,Turcotte D L,Rundle J B. 2005. Aftershock statistics[J]. Pure Appl Geophys,162(6/7):1051–1076
Shcherbakov R,Goda K,Ivanian A,Atkinson G M. 2013. Aftershock statistics of major subduction earthquakes[J]. Bull Seismol Soc Am,103(6):3222–3234 doi: 10.1785/0120120337
Utsu T. 1961. A statistical study on the occurrence of aftershocks[J]. Geophys Mag,30:521–605
Utsu T,Ogata Y,Matsu'ura R S. 1995. The centenary of the Omori formula for a decay law of aftershock activity[J]. J Phys Earth,43(1):1–33 doi: 10.4294/jpe1952.43.1
Wiemer S,Katsumata K. 1999. Spatial variability of seismicity parameters in aftershock zones[J]. J Geophys Res,104(B6):13135–13151 doi: 10.1029/1999JB900032
Wiemer S. 2000. Introducing probabilistic aftershock hazard mapping[J]. Geophys Res Lett,27(20):3405–3408 doi: 10.1029/2000GL011479
Žalohar J. 2014. Explaining the physical origin of Båth's law[J]. J Struct Geol,60:30–45
-
期刊类型引用(1)
1. 李宁，刘纪陆，赵崛. 钢框架结构火灾下抗震性能损伤研究. 建筑钢结构进展. 2022(11): 72-81 . 
百度学术
其他类型引用(1)
下载:
计量
- 文章访问数: 1275
- HTML全文浏览量: 614
- PDF下载量: 63
- 被引次数: 2
