引言

地震发生时，地震预警系统（earthquake early warning system，缩写为EEWs）可以在破坏性地震波到达前向预警目标区发布预警信息，这对于减轻人员伤亡和减少次生灾害至关重要（Allen，Kanomori，2003）。目前，世界上地震较活跃的国家和地区都已经建立或者正在建立地震预警系统，如日本（Kamigaichi et al，2009）、美国（Cua et al，2009）、意大利（Zollo et al，2009）、中国福建（Zhang et al，2016）和中国台湾（Wu，2002）地区等。震级作为地震预警系统中的关键预警信息，其准确性对地震预警系统极为重要。地震预警系统主要依据震级的大小来评估地震对目标区域的影响程度和范围，因此，如何从一个正在发生的地震事件中较为准确地估算出震级，是地震预警系统研发中亟需解决的关键难题之一。

初至地震波（通常为3 s以上）包含了与震级相关的信息（Kanamori，2005；Zollo et al，2006），因此，地震预警中震级估算普遍利用初至地震波的特征参数来建立预测震级的经验公式。这些特征参数在频域和时域内代表初至地震波中与震级相关的信息。在频域内，一般认为震级越大地震波周期越长，据此研究人员提出了卓越周期τ _p （Nakamura，1988）、最大卓越周期${\tau ^{\max }_{\mathrm{p}}} $ （Allen，Kanamori，2003）和平均卓越周期τ_c （Kanamori，2005）等参数；在时域内，相同震源距下震级越大地表变形越大，据此研究人员提出了位移幅值P _d （Wu，Zhao，2006）、速度平方积分（integral of the squared velocity，缩写为IV²）（Festa et al，2008）和累积速度绝对值（cumulative absolute velocity，缩写为 CAV）（Alcik et al，2009）等参数。除了前述早期提出的经典特征参数外，近年来也不断提出新的特征参数，如平均对数周期τ_log （Ziv，2014）和累积位移绝对值（cumulative absolute displacement，缩写为CAD）（马美帅等，2022）等，这些新的特征参数改善了震级预测效果，并促进了震级预测理论的发展。然而，震级估算十分复杂，需要综合考虑震源、传播路径、场地条件、监测仪器等多方面因素的影响（林彬华等，2021），尽管单一特征参数与震级有一定的相关性，但其仅能代表初至地震波中与震级相关的某方面信息，较难反映出地震的全部特征，从而导致震级估算的效果不稳定（Shieh et al，2008；Ziv，2014）。为了改善震级估算效果，一些学者提出采用两个或两个以上特征参数来估算震级，例如将τ_c与P _d相乘来判断地震是否具有破坏性（Wu， Kannmori，2005），或将${\tau ^{\max }_{\mathrm{p}}} $和τ_c各自估算震级的平均值作为最终震级（Hsiao et al，2009），或采用${\tau ^{\max }_{\mathrm{p}}} $和P _d各自估算后取平均震级（Chung et al，2019），以及将频域参数${\tau ^{\max }_{\mathrm{p}}} $和时域参数CAD耦合为新的特征参数S_dτ等（Wang et al，2022）。这些研究均表明了利用更多的特征参数比采用单一特征参数更有利于震级估算。

近年来，随着人工智能的兴起，机器学习在地震预警领域的应用进入了新阶段，在震相捡拾（Perol et al，2018；Wang et al，2021）、地震定位（Perol et al，2018）、事件判别（Li et al，2018）、地震动峰值预测（Hsu et al，2013）等方面都取得了显著的成果。相较人为定义的经验公式，机器学习的优势在于可以自动地建立更为复杂的输入与输出关系，实现更为精准的预测。机器学习的震级估算方法，根据输入的不同可以分为两类：一类是直接将初至地震波作为输入，自动从初至地震波中提取特征，实现端到端的震级估算，以深层卷积神经网络最为典型（Mousavi，Beroza，2020；位栋梁等，2022；Wang et al，2023）；另一类是将五个以上的特征参数（时域参数和频域参数）作为输入，融合初至地震波的多方面信息从而实现震级的估算，主要有浅层人工神经网络（Böse et al，2012；杨黎薇等，2018）、支持向量机（朱景宝等，2021）等方法。尽管当前的机器学习方法显著地改善了震级估算效果，但忽略了地震预警系统估算震级过程中的信息变化。在地震预警系统的实际应用过程中，随着触发台站数量的不断增加，震源距估算是从无到有，从不准确到准确的过程。目前，机器学习方法并不能考虑震源距在地震预警系统的前述变化过程。

本文基于机器学习领域的高斯过程回归（Gaussian processes regression，缩写为GPR）提出了一种多输入的震级估算方法。该方法包括两个版本：无需震源距参数的GPR-M和利用震源距参数以提高估算准确性的GPR-M-R方法。拟利用日本的强震记录对该方法进行训练和测试，并采用智利的强震记录测试GPR方法的泛化能力，之后再通过中国的典型地震事件检验GPR方法的实际应用效果。此外，将GPR方法与广泛使用的${\tau ^{\max }_{\mathrm{p}}} $方法和P _d方法进行对比，以进一步评估其性能，以期为地震预警提供一种更准确、更灵活的震级估算方法。

1.   地震数据

本文选取日本Kiban-Kyoshin Network （KiK-net）和智利SIBER-RISK （Simulation Based Earthquake Risk and Resilience of Interdependent Systems and Networks）数据库的地面加速度记录开展算法的研究和测试，在选取加速度记录时遵循以下原则：① M_W≥4.0；② 因为地震预警系统存在“预警盲区”的问题（郭凯等，2016），为保证初至地震波中至少具有3 s P波用于估算震级（Kanamori，2005；Wu，Zhao，2006；Hsiao et al，2011；金星等，2012；Zhang et al，2016；Cuéllar et al，2018；Kohler et al，2020），且尽可能地包含日本近海地震事件（王延伟等，2020），所以震源距为25—200 km；③ 地震记录的三分向合成加速度的峰值需大于2 cm/s² （Cheng et al，2014），且信噪比不小于10 dB （Perol et al，2018）。对已筛选出的记录进行以下常规处理：① 检查基线并校正；② 采样频率统一为100 Hz；③ 自动捡拾P波到时并人工校验。经以上数据选取和处理后，从日本KiK-net数据库总共筛选出3万3 698条竖向加速度记录，包含了1997年10月21日至2021年12月18日的3 271次M_W4.0—9.0地震事件。将筛选出的KiK-net强震记录数据按照7 ∶ 3的比例随机划分为训练数据集和测试数据集，其中训练数据集含2万3 588条记录，测试数据集含1万零110条记录。从智利SIBER-RISK数据库中总共筛选出5 353条竖向加速度记录，包含了1985年3月3日至2021年7月21日的1 033次M_W4.0—8.8地震事件，用于测试GPR方法的泛化能力。训练、测试和泛化数据集的记录数量的关系如图1所示。

图  1  强震动记录随震级和震源距的分布（每个点代表一条地震动记录）

Figure  1.  Distribution of strong motion records with magnitude and hypocentral distance （each point represents a ground motion record）

下载: 全尺寸图片 幻灯片

2.   高斯过程回归震级估算方法

2.1   高斯过程回归

GPR是一个非参数模型，具有严格的统计学习理论基础，对处理高维、非线性等复杂问题具有很好的适应性，同时GPR还具备避免过拟合、泛化能力强、超参数自适应获取以及输出具有概率意义等优点，这使得GPR成为机器学习领域的重要方法（Seeger，2004；Rasmussen，Williams，2006；Park et al，2011），并在时间序列分析（Rohani et al，2018）、图像处理（Hua et al，2023）和自动控制（Nguyen-Tuong et al，2008）等领域中取得了显著成果。GPR的基本原理是通过协方差核函数的运算，将样本数据映射到高维空间，基于训练数据的训练过程就是对超参数的选择过程。关于GPR的详情请参阅Seeger （2004）与Rasmussen和Williams （2006），此处不再赘述。GPR的网络结构如图2所示。

图  2  GPR的网络结构示意图

x为输入，K（${x^{i}}\text{，} x^{j} $）为协方差函数，y为输出

Figure  2.  Schematic diagram of the GPR network structure

x is the input，K（${x^{i}}\text{，} x^{j} $） is the covariance function，and y is the output

下载: 全尺寸图片 幻灯片

2.2   基于GPR的震级估算

在地震预警系统中，随着触发台站的增加，震源距参数也被不断更新，因此震级估算方法需适应震源距的变化（Wu，Kanomori，2005）。据此，本文在有、无震源距两种情况下估算震级，将特征参数作为输入的震级估算定义为GPR-M方法，将特征参数和震源距作为输入的震级估算定义为GPR-M-R方法，两种方法都将预测的震级作为输出。本文选取频域、时域和时频域三类共十个特征参数作为GPR-M和GPR-M-R方法的输入。

频域类参数估算震级不需要震源距，估算震级的依据是：震级越大，地震波的长周期成分越丰富（Nakamura，1988）。将以下三个典型频域参数作为特征参数输入GPR-M方法：

1） 最大卓越周期$ {\tau ^{\max }_{\mathrm{p}}}$。Allen和Kanamori （2003）提出利用卓越周期τ_p （Nakamura，1988）的最大值估算震级，即${\tau ^{\max }_{\mathrm{p}}} $方法，计算公式为：

式中，$ {X_{{i}}} ＝ \alpha {X_{{{i }- 1}}} ＋ x_{{i}}^{\text{2}} $，$ {D_{{i}}} ＝ \alpha {D_{{{i }－ 1}}} ＋ （ {\text{d}}x/{\text{d}}t ） _{{i}}^{\text{2}} $，τ_p为i时刻测定的卓越周期，x_i和D_i分别为平滑后的竖向速度值及其导数，α为平滑系数，一般取0.99 （Allen，Kanamori，2003）。

2） 平均卓越周期τ_c。Kanamori （2005）提出利用卓越周期τ_p的平均值来估算震级，即τ_c方法，计算公式如下：

式中：u为地震波位移；v为地震波速度；Δt为地震波时长，单位为s。

3） 平均对数周期τ_log。Ziv （2014）提出使用功率谱的平均对数周期τ_log估算震级，其计算公式为

式中，w为功率谱的频率，${P}_{{i}}^{\mathrm{*}} $为0.1—10 Hz频率范围内每0.1个log （频率）的功率谱纵坐标值。

时域类参数估算震级需要已知震源距，估算震级的依据是：在相同震源距情况下，一般震级越大，地表变形越大。以下六个典型时域参数作为特征参数输入GPR-M方法：

1） 位移、速度和加速度的幅值。Wu和Zhao （2006）首先利用初至3 s P波的位移幅值P _d预测震级，并且给出了P _d、震源距和震级的经验公式；速度幅值P _v和加速度幅值P _a与地表变形的最大速度相关（Wu，Kanamori，2005），因此，P _v和P _a与震级存在相关性。

2） 速度平方积分IV²。Festa等（2008）基于地表变形的能量变化提出IV²估算震级，计算公式为：

式中，v为地震波速度，Δt为地震波时长，单位为s。

3） 累积绝对速度值CAV。该参数包含了地表加速度的过程信息，最早被伊斯坦布尔地震预警系统用于判别地震大小（Alcik et al，2009），其计算公式为：

式中a （t）为t时刻的地震波加速度。

4） 累积绝对位移值CAD。该参数包含了地表位移的过程信息，马美帅等（2022）首先提出利用该参数估算震级，其计算公式为：

时频类特征参数在有无震源距时都可以估算震级，估算震级的依据是地震波的时域和频域都包含了与震级相关的信息，综合利于两方面信息更有利于震级估算。Wang等（2022）将CAD和${\tau ^{\max }_{\mathrm{p}}} $相乘得到的时频类特征参数S_dτ用于震级估算，S_dτ的计算公式为

3.   结果及分析

3.1   GPR的训练

由于GPR方法中协方差函数在有限输入点集上要求是正定的，且是一个满足Mercer条件的对称函数，故协方差函数等价于核函数（何志昆等，2013）。常见的协方差函数有径向基核函数、马顿（Matérn）核函数、指数核函数、二次有理核函数等（Seeger，2004；孙斌等，2012）。为选择合适的协方差核函数，将均方误差（mean squared error，缩写为MSE）$ ［ $式（8）$ ］ $作为损失函数。在初至地震波 3 s 时对样本进行训练分析，径向基核函数、马顿核函数、指数核函数和二次有理核函数的MSE分别为0.114 5，0.126 7，0.118 4和0.110 3，其中二次有理核函数的MSE小于其它协方差函数，故本文采用二次有理协方差函数$ ［ $式（9）$ ］ $训练GPR方法。

式中，n为样本数，y_l为估算震级，y_i为实际震级。

式中：M为l ⁻²的对角矩阵，l为关联性测度超参数；$ {\boldsymbol{\sigma}} _{\text{f}}^2 $为核函数的信号方差；α为核函数的形状参数。参数集合$ {\boldsymbol{\theta }}＝ \left\{ {M\text{，} {\boldsymbol{\sigma}} _{\text{f}}^2 \text{，} {\boldsymbol{\sigma}} _{\text{n}}^2} \right\} $为包含所有超参数的向量，通过极大似然法求得。首先建立训练样本条件概率的负对数似然函数$L （ {\boldsymbol{\theta}} ） ＝ － {\log {{p}}} （ {{\boldsymbol{y}} \,|{\boldsymbol{X}} \text{，} {\boldsymbol{\theta}} } ） $，并令其对超参数θ求偏导，然后采用共轭梯度法对偏导数进行最小化以得到超参数的最优解，即：

式中，tr （）表示矩阵的迹。

3.2   确定P _d及$ {{{{\boldsymbol{\tau}} ^{{\mathbf{max}} }_{\mathbf{p}}}}}$的经验公式

利用训练数据集拟合震级与P _d和${\tau ^{\max }_{\mathrm{p}}} $的经验公式，如式（12）所示。选取P _d和${\tau ^{\max }_{\mathrm{p}}} $作对比的原因是：${\tau ^{\max }_{\mathrm{p}}} $方法与震级的经验公式中未使用震源距参数，与GPR-M方法的参数选择具有共同因子，且频域参数震级估算方法的对比结果表明${\tau ^{\max }_{\mathrm{p}}} $方法估算震级的准确性高于其它频域参数方法（Nazeri et al，2017；王延伟等，2020），该方法已经应用于地震预警系统 （Hsiao et al，2009；Chung et al，2019）；P _d方法与震级的经验公式中使用了震源距参数，与GPR-M-R方法的参数选择具有共同因子，且特征参数震级估算方法的对比研究表明P _d方法估算震级的准确性明显优于其它参数方法（Leyton et al，2018；王延伟等，2020），该方法也被广泛应用于地震预警系统 （Hsiao et al，2009）。据此本文选择了以上两种特征参数方法进行对比。

式中：Para为特征参数（P _d或${\tau ^{\max }_{\mathrm{p}}} $）；M为震级；R为震源距；a，b，c为待定的回归系数；${\tau ^{\max }_{\mathrm{p}}} $不需要震源距R项时，c＝0；P _d需要震源距R项时，c ≠ 0。

地震预警系统中的震级估算是随时间持续更新的过程，利用更长的初至地震波可以提高震级估算的准确性（Kanamori，2005；Wu，Zhao，2006；Ziv，2014；Peng et al，2017），初至P波到达3—10 s时，特征参数的经验公式可由训练数据集得到，利用最小二乘法回归得到经验公式（12）的回归系数，如表1所示。

表  1  初至3—10 s地震波的P_d和$ {\tau ^{\max }_{\mathrm{p}}} $方法的拟合系数表

Table  1.  Fitting coefficients for the P_d and $ {\tau ^{\max }_{\mathrm{p}}} $ methods based on the first 3−10 s seismic waves

初至地震波 窗长/s	$P_{{\mathrm{d}}} $方法				$ {\tau ^{\max }_{\mathrm{p}}}$方法
	a	b	c		a	b
3	0.639 4	−3.987 2	−0.840 7		0.309 6	−2.049 4
4	0.681 4	−3.910 9	−0.966 4		0.320 6	−2.048 3
5	0.708 0	−3.740 4	−1.104 4		0.329 8	−2.041 8
6	0.723 6	−3.518 2	−1.248 1		0.339 8	−2.016 5
7	0.738 2	−3.292 8	−1.378 0		0.348 6	−1.975 3
8	0.740 9	−3.070 7	−1.478 3		0.352 6	−1.955 5
9	0.737 8	−2.897 5	−1.538 8		0.358 3	−1.911 7
10	0.733 4	−2.791 4	−1.561 7		0.364 7	−1.875 0

下载: 导出CSV 

| 显示表格

3.3   利用初至3 s P波测试数据集的震级估算结果

在地震预警系统的研究和应用中，普遍认为初至3 s P波估算的震级较好地兼顾了地震预警的时效性和准确性（Kanamori，2005；Wu，Zhao，2006；Peng et al，2017；Wang et al，2022）。因此，改善初至3 s P波估算震级的效果至关重要。图3为GPR-M方法与${\tau ^{\max }_{\mathrm{p}}} $方法估算震级的误差分布及其直方图。从误差的定性分布来看：GPR-M方法估算震级的误差离散性明显小于${\tau ^{\max }_{\mathrm{p}}} $方法且误差更多地分布在$ ［ $−0.5， 0.5$ ］ $之内，GPR-M方法比${\tau ^{\max }_{\mathrm{p}}} $方法估算震级的误差更小、误差更接近于零。随着震级的增加，两者对M_W≥6.5地震事件的震级估算偏小（误差分布下倾），使得M_W4.0—6.4范围内的准确率明显高于M_W6.5—9.0范围内的准确率。这与众多研究中的震级饱和问题是一致的（Wu，Zhao，2006；王延伟等，2020；Wang et al，2022）。从误差的定量结果看：GPR-M方法的误差标准差、误差均值和误差绝对值均值明显低于${\tau ^{\max }_{\mathrm{p}}} $方法（GPR-M方法分别为0.353 9，−0.091 0和0.274 7，${\tau ^{\max }_{\mathrm{p}}} $方法分别为0.745 5，−0.142 8和0.687 4）；GPR-M方法的误差标准差比${\tau ^{\max }_{\mathrm{p}}} $方法降低了约52.53%；在M_W4.0—6.4范围内，GPR-M方法震级估算的准确率约是${\tau ^{\max }_{\mathrm{p}}} $方法的1.59倍（GPR-M方法为88.75%，${\tau ^{\max }_{\mathrm{p}}} $方法为55.94%）；在M_W6.5—9.0范围内，即使两者表现出震级饱和现象，GPR-M方法震级估算的准确率也约是${\tau ^{\max }_{\mathrm{p}}} $方法的1.41倍（GPR-M方法为47.89%，${\tau ^{\max }_{\mathrm{p}}} $方法为33.98%）。

图  3  初至3 s P波时GPR-M （a）和$ {\tau ^{\max }_{\mathrm{p}}}$（b）估算震级的误差分布及其直方图

图中的散点表示每条记录估算震级的误差，蓝色实线和红色虚线表示±0.5，σ为误差标准差，μ为误差均值，ω为误差绝对值均值，准确率定义为误差在$ ［ $−0.5， 0.5$ ］ $的记录数与记录总数的比值，下同

Figure  3.  Distribution of estimated magnitude errors by GPR-M （a） and $ {\tau ^{\max }_{\mathrm{p}}}$ （b） at initial 3 s P-wave

Scatter points show magnitude estimation errors based on per record. Blue solid lines and red dashed lines indicate ±0.5 . σ is error standard deviation，μ is mean error，ω is mean absolute error，accuracy is the ratio of records with errors in $ ［ $−0.5，0.5$ ］ $ to total records，the same below

下载: 全尺寸图片 幻灯片

图4为GPR-M-R方法与P _d方法估算震级的误差分布及其直方图。从误差的定性分布来看：GPR-M-R方法估算震级的误差更多地分布在$ ［ $−0.5， 0.5$ ］ $之内，随着震级的增加，对于M_W≥6.5的地震事件同样表现出了震级低估现象。从误差的定量结果看：GPR-M-R方法的误差标准差、误差均值和误差绝对值均值都优于P _d方法（GPR-M-R方法分别为0.338 4，0.078 1和0.240 7；P _d方法分别为0.570 2，−0.094 6和0.444 3）；GPR-M-R方法的误差标准差比P _d方法降低了约40.65%；在M_W4.0—6.4范围内，GPR-M-R方法震级估算的准确率是P _d方法的1.45倍（GPR-M-R方法为90.54%，P _d方法为62.28%）；在M_W6.5—9.0范围内，GPR-M-R方法震级估算的准确率是P _d方法的1.06倍（GPR-M-R方法为48.67%，P _d方法为45.89%）。

图  4  初至3 s P波GPR-M-R （a）和P _d （b）方法估算震级的误差分布

Figure  4.  Distributions and histograms of estimated magnitude errors for GPR-M-R （a） and P _d （b） for initial 3 s P-wave

下载: 全尺寸图片 幻灯片

3.4   不同初至地震波时长的测试数据集的震级估算结果

将初至地震波时长从3 s增加至10 s，测试GPR-M方法与$ {\tau ^{\max }_{\mathrm{p}}}$方法、GPR-M-R方法与P _d方法持续估算震级的效果。图5a为GPR-M方法与$ {\tau ^{\max }_{\mathrm{p}}}$方法所估算震级的误差标准差随初至地震波时长的变化，可见：两种方法的误差标准差随着初至地震波时长的增加逐步减小，但GPR-M方法的误差标准差始终小于$ {\tau ^{\max }_{\mathrm{p}}}$方法，GPR-M方法的误差标准差比$ {\tau ^{\max }_{\mathrm{p}}}$方法降低了52.53%—57.21%。值得注意的是，GPR-M方法在3 s时的误差标准差（0.353 9）仍小于$ {\tau ^{\max }_{\mathrm{p}}}$方法在10 s时的误差标准差（0.668 6）。图5b为GPR-M-R方法与P _d方法所估算震级的误差标准差随初至地震波时长的变化，可见：两种方法的误差标准差同样随着初至地震波时长的增加逐步减小，GPR-M-R方法的误差标准差始终小于P _d方法，GPR-M-R方法的误差标准差比P _d方法降低了约37.72%—40.65%，并且GPR-M-R方法在3 s时的误差标准差（0.338 4）小于P _d方法在10 s时的误差标准差（0.404 8）。

图  5  初至3—10 s地震波的不同方法估算震级的误差标准差对比

（a） GPR-M与$ {\tau ^{\max }_{\mathrm{p}}}$方法对比；（b） GPR-M与P _d方法对比

Figure  5.  Comparison of the standard deviation of errors in estimating magnitudes for initial 3−10 s seismic waves

（a） Comparison of GPR-M method with $ {\tau ^{\max }_{\mathrm{p}}}$ method ；（b） Comparison of GPR-M method with P _d method

下载: 全尺寸图片 幻灯片

图6和图7给出了GPR-M方法与$ {\tau ^{\max }_{\mathrm{p}}}$方法。GPR-M-R方法与P _d方法估算不同范围震级的准确率随初至地震波时长的变化，可见：在M_W4.0—6.4范围内，GPR-M方法的准确率由88.75%提高至91.68%，$ {\tau ^{\max }_{\mathrm{p}}}$方法的准确率由55.94%提高至61.23%，GPR-M方法的准确率始终约是$ {\tau ^{\max }_{\mathrm{p}}}$的1.50—1.59倍（图6a）；在M_W6.5—9.0范围内，GPR-M方法的准确率由47.89%提高至71.45%，$ {\tau ^{\max }_{\mathrm{p}}}$方法的准确率由33.98%提高至46.81%，GPR-M方法的准确率始终约为$ {\tau ^{\max }_{\mathrm{p}}}$方法的1.41—1.53倍（图6b）；在M_W4.0—6.4范围内，GPR-M-R方法的准确率在由90.54%提高至94.80%，P _d方法的准确率由62.28%提高到79.02%，GPR-M-R方法的准确率始终是P _d方法的1.20—1.45倍 （图7a）；在M_W6.5—9.0范围内，GPR-M-R方法的准确率由48.67%提高至78.12%，P _d方法的准确率由45.89%提高至68.17%，GPR-M-R方法的准确率约为P _d方法的1.06—1.15倍（图7b）。综上表明，GPR-M方法和GPR-M-R方法估算震级准确率的稳定性高于$ {\tau ^{\max }_{\mathrm{p}}}$方法和P _d方法。

图  6  初至3—10 s地震波的GPR-M方法和$ {\tau ^{\max }_{\mathrm{p}}}$方法估算不同范围震级的准确率随时间变化趋势的对比

Figure  6.  Comparison of the trends in estimation accuracy of different magnitude ranges based on the GPR-M and $ {\tau ^{\max }_{\mathrm{p}}}$ methods using initial 3−10 s seismic waves

（a） 4.0≤M_W≤6.4；（b） 6.5≤M_W≤9.0

下载: 全尺寸图片 幻灯片

图  7  初至3—10 s地震波的GPR-M-R方法和P _d方法估算不同范围震级的准确率随时间变化趋势的对比

Figure  7.  Comparison of the trend of magnitude estimation accuracy change over time by the GPR-M-R and P _d methods using initial 3−10 s seismic waves

（a） 4.0≤M_W≤6.4；（b） 6.5≤M_W≤9.0

下载: 全尺寸图片 幻灯片

3.5   泛化数据集的震级估算结果

上述测试中，训练数据集和测试数据集均采用了日本的地震记录，两数据集都包含一些相似的区域信息，例如相似的震源、传播路径、场地条件、监测仪器等，这些区域信息可能会使得GPR-M和GPR-M-R方法仅适用于日本地区，而无法用于其它地区。因此，为了检验日本地震记录训练的GPR-M方法和GPR-M-R在其它地区估算震级的效果，需要测试其泛化能力。泛化数据集选取智利的地震记录，因为智利的地震记录没有用于训练GPR-M方法和GPR-M-R方法，也未被用于拟合$ {\tau ^{\max }_{\mathrm{p}}}$和P _d方法的经验公式。考虑到误差随震级的分布可以最为详细地展示震级估算结果，故给出了泛化数据集的震级估算误差分布图（图8，图9），以P波到达后3 s，5 s，8 s，10 s的结果为例。

图  8  P波到达后3 s，5 s，8 s，10 s GPR-M （a）和$ {\tau ^{\max }_{\mathrm{p}}}$ （b）方法估算泛化数据集震级的误差分布

Figure  8.  Distributions of estimated magnitude errors in the generalized dataset using the GPR-M （a） and $ {\tau ^{\max }_{\mathrm{p}}}$ （b） at the P-wave first arrival times of 3 s，5 s，8 s，and 10 s

下载: 全尺寸图片 幻灯片

图  9  P波到达后3 s，5 s，8 s，10 s GPR-M-R （a）和P _d （b）估算泛化数据集震级的误差分布

Figure  9.  Distributions of estimating magnitude errors in the generalized dataset using the GPR-M （a） and P _d （b） at the P-wave arrival times of 3 s，5 s，8 s，and 10 s

下载: 全尺寸图片 幻灯片

由图8和图9可以看出GPR-M方法和GPR-M-R方法估算震级的误差分布更集中在$ ［ $−0.5， 0.5$ ］ $范围内，两者估算震级的误差标准差、误差均值和误差绝对值均值分别优于$ {\tau ^{\max }_{\mathrm{p}}}$方法和P _d方法，GPR-M方法的误差标准差比$ {\tau ^{\max }_{\mathrm{p}}}$方法降低了53.08%—55.13%，GPR-M-R方法的误差标准差比P _d方法降低了35.88%—36.59%，这表明GPR-M方法和GPR-M-R方法的泛化能力和估算震级的稳定性分别优于$ {\tau ^{\max }_{\mathrm{p}}}$方法和P _d方法。

3.6   震例测试

为了检验GPR震级估算方法对于实际震例的可靠性与合理性，选取了我国的三次典型震例进行验证分析。需要注意的是，我国常用震级为M_S，对于M_W＜8.5的地震，M_W与M_S之间的差异可以忽略不计（Peng et al，2017），因此，本文给出的GPR，$ {\tau ^{\max }_{\mathrm{p}}}$和P _d方法都可以直接用于我国地震的震级估算。选取的三次震例分别为：2019年6月22日四川宜宾M_S5.4地震，震中位置为（28.4°N，104.8°E），震源深度为10 km；2021年5月21日云南漾濞M_S6.4地震，震中位置为（25.7°N，99.9°E），震源深度为8 km；2008年5月12日四川汶川M_S8.0地震，震中位置为（31.0°N，103.4°E），震源深度为14 km。三次地震的震中位置及所用的最先触发的四个台站分布如图10所示。在未知震源距和已知震源距两种情况下，将首个监测台站的P波到时记为0时刻，对每个触发台站P波到达后3 s 至 10 s，以1 s间隔持续估算震级。

图  10  用于检验GPR震级估算效果所选取的我国三次地震的震中及所用台站分布

Figure  10.  Distribution of epicenters and the stations used for testing the effectiveness of GPR magnitude estimation for three earthquakes in China

下载: 全尺寸图片 幻灯片

图11为震源距未知、三次地震最先触发四个台站时，利用$ {\tau ^{\max }_{\mathrm{p}}}$方法和GPR-M方法持续估算震级的结果。可以看出，$ {\tau ^{\max }_{\mathrm{p}}}$方法明显低估了三次地震的震级，并且估算震级偏差明显大于GPR-M方法：对于宜宾M_S5.4地震，采用$ {\tau ^{\max }_{\mathrm{p}}}$方法的最大震级偏差为−1.59，采用GPR-M为−0.54；对于漾濞M_S6.4地震，采用$ {\tau ^{\max }_{\mathrm{p}}}$方法的最大偏差为−2.34，采用GPR-M为−0.39；对于汶川M_S8.0地震，采用$ {\tau ^{\max }_{\mathrm{p}}}$方法的最大偏差为−3.08，采用GPR-M为−1.40。此外，GPR-M估算震级的离散性比$ {\tau ^{\max }_{\mathrm{p}}}$方法更小。

图  11  使用GPR-M方法（左）和$ {\tau ^{\max }_{\mathrm{p}}}$方法（右）持续估算我国2019年宜宾M_S5.4 （a）、2021年漾濞M_S6.4 （b）和2008年汶川M_S8.0 （c）地震震级的结果

Figure  11.  Results of continuous magnitude estimation for the 2019 Yibin M_S5.4 （a），2021 Yangbi M_S6.4 （b），and 2008 Wenchuan M_S8.0 （c） earthquakes using the GPR-M-R （left panels） and $ {\tau ^{\max }_{\mathrm{p}}}$ （right panels） methods

下载: 全尺寸图片 幻灯片

图12为震源距已知、三次地震最先触发四个台站时，用P _d方法和GPR-M-R持续估算震级的结果。可以看出，GPR-M-R估算的震级更接近实际震级，且偏差小于P _d方法：对于宜宾M_S5.4地震，GPR-M-R的最大偏差为＋ 0.47，P _d方法为−0.56；对于漾濞M_S6.4 地震，GPR-M-R的最大偏差为−0.27，P _d方法为−0.57；对于汶川M_S8.0 地震，GPR-M-R的最大偏差与P _d方法相同，均为−1.59，但GPR-M-R估算的震级更接近实际震级，且离散性更小。

图  12  使用GPR-M-R （左）和P _d （右）方法持续估算我国2019年宜宾M_S5.4 （a）、2021年漾濞M_S6.4 （b）和2008年汶川M_S8.0 （c）地震震级的结果

Figure  12.  Results of continuous magnitude estimation for the 2019 Yibin M_S5.4 （a），2021 Yangbi M_S6.4 （b），and 2008 Wenchuan M_S8.0 （c） earthquakes using the GPR-M-R （left panels） and P _d （right panels） methods

下载: 全尺寸图片 幻灯片

4.   讨论与结论

为改善地震预警系统中震级估算的准确性，本文提出了能够同时适用于有震源距和无震源距两种情况下进行震级估算的GPR震级估算方法。该方法以时域、频域以及时频域多个特征参数作为输入（GPR-M），或以特征参数和震源距作为输入（GPR-M-R），以震级作为输出，通过融合初至地震波中多方面信息实现震级估算。GPR-M方法在未使用震源距的情况下提高了震级估算的准确性，GPR-M-R方法在使用震源距的情况下进一步提高震级估算的准确性。在这两种情况下，利用日本的强震记录进行训练和测试，利用智利的强震记录进行泛化能力测试，利用我国的三次震例进行可靠性和合理性检验，并与当前广泛采用的$ {\tau ^{\max }_{\mathrm{p}}}$方法和P _d方法进行了对比，结论如下：

1） 利用日本测试数据集初至3—10 s地震波估算震级时，GPR-M方法和GPR-M-R方法估算震级效果优于$ {\tau ^{\max }_{\mathrm{p}}}$方法和P _d方法，GPR-M方法估算震级的误差标准差比$ {\tau ^{\max }_{\mathrm{p}}}$方法降低了52.53%—61.20%，GPR-M-R方法比P _d方法降低了37.72%—41.21%。

2） 对于震级较大的地震（M_W≥6.5），GPR-M方法和GPR-M-R方法的震级饱和现象无$ {\tau ^{\max }_{\mathrm{p}}}$方法和P_d方法明显，GPR-M方法估算震级的准确率比$ {\tau ^{\max }_{\mathrm{p}}}$方法提高了1.4—1.5倍，GPR-M-R方法估算震级的准确率比P _d方法提高了1.2—1.45倍。

3） 基于智利强震记录的泛化能力测试结果表明，GPR-M方法和GPR-M-R方法震级估算的准确性相较于$ {\tau ^{\max }_{\mathrm{p}}}$方法和P _d方法更高，泛化能力更强，受区域特性影响更小。

4） GPR-M方法和GPR-M-R方法较好地估算了我国三次地震的震级，且准确度高于$ {\tau ^{\max }_{\mathrm{p}}}$方法和P _d方法，表明基于GPR的震级估算方法应用于我国地震预警系统具有较好的可靠性与合理性。

本文所提出的GPR震级估算方法在测试数据集和泛化数据集中表现出了较好的震级估算效果，有效地提高了震级估算的准确性。但对于震级较大的地震，震级饱和现象依然存在。震级饱和问题已经成为了震级估算研究中的学术难题，因为大地震的断层破裂过程有几十秒至上百秒，且破裂过程十分复杂，是否可以使用初始几秒的断层破裂信息来确定最终破裂的规模在理论上仍无统一结论（Olson，Allen，2005；Ide，2019）。此外，本文在进行训练GPR方法时，训练数据集里中小地震记录（4.0≤M_W≤6.4，2万2 207条）的占比远大于大地震（6.5≤M_W≤ 9.0，1 381条），数据的不均匀性是否影响了GPR-M方法的震级估算效果，也需要进一步研究。尽管GPR方法在智利地震和我国震例的测试中表现出良好的震级估算效果，但区域性的差异（如传播介质和场地条件等）对GPR方法的影响程度尚待明确。

日本防灾科学技术研究所、智利SIBER-RISK数据库和中国地震局工程力学研究所强震动观测中心为本研究提供数据支持，审稿专家为本文提出了诸多宝贵的意见，作者在此一并表示感谢。