基于支持向量机的现地地震预警地震动峰值预测

余 聪，宋晋东，李山有

(1.中国地震局工程力学研究所，哈尔滨 150080;2.中国地震局地震工程与工程振动重点实验室，哈尔滨 150080)

地震是对人类生存安全危害巨大的一种突发性自然灾害。地震预警是在地震发后，利用震源附近地震台站观测到的地震波初期信息，快速估计地震参数并预测地震对周边地区的影响，抢在破坏性地震波到达震中周边地区之前，发布各地地震动强度和到达时间等预警信息，能够给重要的工程或人员密集地区提供几秒到数十秒的预警时间，以减轻地震灾害[1]。地震预警已被证明是实时减灾的有效手段之一，经过几十年的发展，世界上多个地震活跃的国家已运行或正在发展地震预警系统，例如墨西哥[2]、日本[3-5]、美国[6-9]、意大利[10-12]、中国[13-14]。

地震预警可以分为异地预警与现地预警[15]。异地预警在潜在震源区布设高密度的地震台网，利用近震源台(P波或S波)记录信息来确定地震参数，进而估计预警目标区地震动，并对地震潜在破坏区发出警报；现地预警原理是利用初至P波和后至破坏性地震波(S波或者面波)的到时差，Kanamori[16]认为P波携带地震信息，因此可从P波中预测后续的峰值地震动，如PGA与PGV，通常的方法是建立单个P波特征参数与目标地震动参数的线性关系进行后续的地震动预测，例如Wu等[17]、Zollo等[18-19]使用的Pd方法，即采用P波触发后3 s内的垂直方向位移峰值Pd预测后续地震动峰值PGV，这是地震预警最常用的一种现地地震动预测方法。此外，Nakamura[20]采用破坏烈度DI方法，Festa等[21]、Brondi等[22]采用的速度平方积分IV2方法等预测后续地震动。但是，利用单个P波参数预测后续地震动只利用了单一的地震动初期信息，Böse等[23]利用人工神经元网络、基于多频带多参数实现了地震预警PGA、PGV的持续预测方法，为综合利用多种地震波初期信息进行现地地震动预测提供了新的思路。

本文以0.1～10 Hz带通滤波后三分向矢量合成地震动峰值PGA与PGV为预测目标参数，利用日本K-net强震台网P波触发后3 s数据，基于加速度幅值Pa、速度幅值Pv、位移幅值Pd、傅里叶谱幅值AMmax、速度平方积分IV2、破坏烈度DI、累积绝对速度CAV、阿里亚斯烈度Ia这8种特征参数构建支持向量机(Support Vector Machine，SVM)预测模型，旨在利用支持向量机样本需求量小、泛化能力强和抗干扰能力强等优点，探索多参数预测后续地震动的新方法，为现地地震预警仪器地震烈度预测提供服务。

1 数 据

本研究从日本K-net强震观测台网中选取地震数据，并依据以下原则筛选[24]:

(1)时间为自2007年10月至2019年8月；

(2)发生于陆地的地震，因为我国的地震多发生于陆地；

(3)震源深度为10 km以内的地震，因为我国大部分地区发生的是浅源地震，同时，在震级相同的情况下，浅源地震会比深源地震造成更大的破坏性；

(4)震中距取100 km以内的强震数据，对5级以上的地震，补充100～200 km范围的数据。

基于上述原则共筛选出2 068次地震共25 308组(75 924条)强震动数据，竖直向原始加速度记录未滤波幅值范围0.14～665.78 cm/s2，水平向原始加速度记录未滤波幅值范围0.78～1 538.16 cm/s2。随机选取训练数据集与测试数据集，比例为7:3。地震震中及台站的分布，如图1所示。震级、震中距、记录数量分布，如图2所示。

图1 本研究训练数据集与测试数据集的地震震中分布(空心圆直径与震级大小成正比)K-net台站(灰色三角形)分布

图2 本研究训练数据集与测试数据集的震中距及震级分布图(三角形为训练数据，圆圈为测试数据)

对上述数据进行如下处理：

(1)采用马强等[25]的方法自动识别P波震相及到时，即先利用长短时平均法对P波震相初步识别、再利用AIC准则精确捡拾P波到时。为确保P波到时的准确性，对自动识别的地震P波到时进行人工校正；

(2)对加速度积分，积分一次得到速度记录、积分两次得到位移记录。为了消除积分后的低频漂移现象，对未处理加速度记录和积分后的记录作0.075 Hz高通滤波。

(3)对加速度记录与积分后的速度记录做0.1～10 Hz带通滤波，进行三分向矢量合成计算得到地震动峰值PGA与PGV。

2 方 法

2.1 预测目标参数

地震烈度作为标定地震引起的地震动及其影响的强弱程度，具有度量地震破坏原因和后果的双重性，其快速的评定对于震后灾害评估、应急救援有着十分重要的作用。但是限于烈度本身的定义模式，震后难以组织人员对受灾区域进行快速评定。随着地震观测技术的持续发展及地震观测仪器数量的不断增加，根据地震观测仪器获取的地震动记录快速定量地计算仪器地震烈度成为可能。仪器地震烈度的计算是地震烈度速报的基础，可为灾情快速判断、地震应急救援决策和行动、工程抢险修复决策等提供科学依据[26]。

《中国地震烈度表》[27]给出了仪器地震烈度的计算过程：

获取地震动记录所用观测仪器应可为三方向(条件不具备时可为正交两水平方向或单水平方向)，地震动加速度和速度记录的每个分向均应采用数字滤波器进行0.1～10 Hz带通滤波，PGA、PGV分别为滤波后三分向的加速度、速度矢量合成的最大值，那么应用PGA得到的仪器地震烈度计算值IPGA可由式(1)得到，应用PGV得到的仪器地震烈度计算值IPGV可由式(2)得到，进而可由IPGA和IPGV根据式(3)计算出仪器地震烈度II。

IPGA=

(1)

IPGV=

(2)

II=

(3)

2.2 特征参数

本研究采用如下8种特征参数作为支持向量机预测模型的输入参数，且计算时间窗均为P波触发后3 s。

(1)加速度幅值Pa：P波触发后3 s内的垂直方向加速度峰值，该参数表示地震动的加速度幅值特征。

(2)速度幅值Pv：P波触发后3 s内的垂直方向加速度峰值，该参数表示地震动的速度幅值特征。

(3)位移幅值Pd：P波触发后3 s内的垂直方向位移峰值，该参数表示地震动的位移幅值特征。

(4)傅里叶谱幅值AMmax：P波前3 s垂直方向记录傅里叶谱的幅值，该参数表示地震动的频谱特征。

(5)速度平方积分IV2(Integral of the Squared Velocity)，该参数与地震早期辐射能量紧密相关，同时它也可用于预测地震动[21]。IV2的计算公式如式(4)，v为垂直向速度记录，tc为P波到时点，Δt为积分计算的时间窗长。

(4)

(6)破坏烈度DI(Destructive Intensity)：通过地震动竖向加速度a和速度v的内积计算，该参数可估计地震的危险程度。DI的计算方法如式(5)，由式(5)可见，DI值表示的是“功率”的概念，DI取加速度和速度的内积绝对值的对数。

DI=lg|a·v|

(5)

(7)累积绝对速度CAV(Cumulative Absolute Velocity)：由地震动加速度的绝对值对时间间隔积分得到的[28]，该参数考虑了地震动的幅值与持时特征，CAV的由式(6)计算:

(6)

(8)阿里亚斯烈度Ia(Arias Intensity)由Arias[29]提出,该参数同时包含了地震动的幅值、频率及持时特征，计算方法如下：

(7)

2.3 支持向量机

支持向量机是人工智能的机器学习方法中一种基于统计学习理论的监督学习方法，能非常出色解决多变量回归问题，其原理为将输入数据映射到高维特征空间,将非线性问题变化为线性问题，通过变换后的线性问题的方法求解原来的非线性问题。

对于给定训练数据D={(xi,yi)},i=1,2，…,mxi∈Rd,yi∈R，支持向量机找到一个回归函数(式(8))，使得f(x)与y尽可能地接近。支持向量回归模型与传统回归模型不同之处在于其容许f(x)与y之间最多有2ε误差，即f(x)与y之间的误差

f(x)=ωTφ(x)+b

(8)

(9)

式中：C为正则化常数，lε为ε—不敏感损失(ε-insensitive loss)函数

lε=

(10)

s.t.f(xi)-yi≤ε+ξi

(11)

(12)

(13)

上述过程需要满足KKT(Karush-Kuhn-Tucker)条件，用核函数K(xi,xj)替代[φ(xi)·φ(xj)]，即可得到预测函数模型：

(14)

3 结 果

本文在研究中遵循以下流程：

(1)以0.1～10 Hz带通滤波后三分向矢量合成地震动峰值PGA与PGV为预测目标参数，且计算时间窗均为P波触发后3 s；

(2)利用训练数据集，以加速度幅值Pa、速度幅值Pv、位移幅值Pd、傅里叶谱幅值AMmax、速度平方积分IV2、破坏烈度DI、累积绝对速度CAV、阿里亚斯烈度Ia这8种特征参数取其所有组合做为输入参数，训练后建立预测模型，并选择其中误差标准差最小的参数组合做为本文的支持向量机预测模型SVM-PGA与SVM-PGV；

(3)利用训练数据集，训练后建立Pd预测模型Pd-PGA与Pd-PGV，形式分别为如下线性形式：

lg(PGA)=k1·lg(Pd)+b1

(15)

lg(PGV)=k2·lg(Pd)+b2

(16)

式中，k1、b1、k2、b2为拟合系数。

(4)利用测试数据集，对比支持向量机预测模型与Pd预测模型的预测结果；

(5)随机选取两次地震事件(不在训练数据集与测试数据集范围内)，对比支持向量机预测模型与Pd预测模型的预测结果。

3.1 训练数据集与测试数据集

利用本文选取的日本K-net强震观测台网2007年10月至2019年8月的强震数据，按照70%数据为训练集进行模型训练，并在剩余30%数据上进行模型性能测试。

基于python3.5环境搭建支持向量机模型，选用高斯核函数，停止训练的误差精度取0.001，选取正则化参数C和核函数参数gamma并采用网格搜索来选取最优组合，其余参数选择默认值。

由于训练数据集是在对数坐标下得出预测模型，因此将实测值对数与预测值对数之差定义为误差xi，如式(17)，并计算误差标准差σ，如式(18)：

xi=lg(实测值)-lg(预测值)

(17)

(18)

表1 不同参数组合模型在训练数据集上预测PGA的误差标准差

图3 不同参数组合在训练数据集上预测PGA的误差标准差

依据训练数据集建立的Pd预测模型如下：

lg(PGA)=0.567lg(Pd)+2.238

(19)

lg(PGV)=0.740lg(Pd)+1.213

(20)

表2列出了基于训练数据集的支持向量机预测模型SVM-PGA与SVM-PGV、以及Pd预测模型Pd-PGA与Pd-PGV的误差标准差σ计算结果。支持向量机预测模型SVM-PGA与SVM-PGV的误差标准差σ分别为0.240与0.254，Pd预测模型Pd-PGA与Pd-PGV的误差标准差σ分别为0.320与0.309。在训练数据集上，支持向量机预测模型SVM-PGA与SVM-PGV的误差标准差均分别小于Pd预测模型Pd-PGA与Pd-PGV。

表2 支持向量机预测模型与Pd预测模型在训练数据集上误差标准差

图4为支持向量机预测模型SVM-PGA与SVM-PGV的误差概率密度曲线、以及Pd预测模型Pd-PGA与Pd-PGV的误差概率密度曲线。可以看到，这四种预测模型的误差概率密度曲线均呈正态分布。图4中可以得到，支持向量机预测模型SVM-PGA与SVM-PGV的误差概率密度曲线形状比Pd预测模型Pd-PGA与Pd-PGV的误差概率密度曲线形状更高、更窄，这表明支持向量机预测模型SVM-PGA与SVM-PGV比Pd预测模型Pd-PGA与Pd-PGV的误差更小、误差更趋近于零。

(a)PGA

基于训练数据集得到预测模型后，使用支持向量机预测模型SVM-PGA与SVM-PGV、Pd预测模型Pd-PGA与Pd-PGV在测试数据集上进行测试。图5表示了预测值与实测值的分布情况，1∶1的实线表示预测值等于实测值，虚线表示训练模型得到的1倍标准差。从图5(a)和图5(c)可以得到，支持向量机预测模型SVM-PGA与SVM-PGV的误差σ分别为0.240与0.257，与训练数据集得到的结果几乎一致。从图5(b)和图5(d)可以得到，Pd预测模型Pd-PGA与Pd-PGV的误差标准差σ分别为0.313与0.304，小于训练数据集得到的结果。支持向量机预测模型SVM-PGA与SVM-PGV的误差标准差均分别小于Pd预测模型Pd-PGA与Pd-PGV，这一点与训练数据集得到的结果一致。

同时，图5的结果也显示，Pd预测模型在整体震级范围内小数值区域出现高估现象、在高数值区域出现低估现象，即“小值高估、大值低估”，支持向量机预测模型在6.5级以下地震时整体趋近1∶1线性比例关系、6.5级以上地震出现了轻微的低估现象。

(a)SVM-PGA

图6显示了支持向量机预测模型误差以及Pd预测模型误差随震中距变化的分布。

图6(a)与图6(c)分别显示了支持向量机预测模型SVM-PGA与SVM-PGV的预测误差随震中距变化的趋势，可以看到，支持向量机预测模型SVM-PGA与SVM-PGV的预测误差在0～200 km距离区间均匀的分布在零误差(虚线)两侧，这表明支持向量机预测模型的预测结果不受到震中距变化的影响，且整体误差位于(-0.6，0.8)区间，这表明依据本文提出的支持向量机预测模型计算仪器地震烈度，最大误差为0.8烈度单位。

(a)SVM-PGA

图6(b)与图6(d)分别显示了Pd预测模型Pd-PGA与Pd-PGV的预测误差随震中距变化的趋势，可以发现，Pd预测模型Pd-PGA的预测误差明显随震中距逐渐增大而减小，Pd预测模型Pd-PGV的预测误差也有随震中距逐渐增大而减小的趋势，且Pd预测模型Pd-PGA与Pd-PGV的整体误差位于(-1.0，1.0)区间，这表明依据本文提出的支持向量机预测模型计算仪器烈度，最大误差为1.0烈度单位。

3.2 震例分析

为了解支持向量机预测模型在具体地震事件上预测PGA、PGV的性能，随机选取两次不在训练数据集和测试数据集范围内的地震事件进行震例分析，这两次地震的震中位置及台站分布如图7所示，其中，图7(a)为2002年9月16日Mj5.3级地震，共计60组强震动记录，图7(b)为2007年3月25日Mj6.9级地震，共计93组强震动记录。

图7 震例分析的地震震中(五角星)及K-net台站(三角形)分布

图8为支持向量机预测模型和Pd预测模型在2002年9月16日Mj5.3级地震的预测值与实测值的分布，1∶1的实线表示预测值等于实测值，虚线表示基于训练数据集得到的1倍误差标准差。图8(a)与图8(c)分别为SVM-PGA与SVM-PGA预测模型的预测结果，可以看到预测值与实测值趋近1∶1线性比例关系，绝大部分数据分布在1倍误差标准差范围内。图8(b)和图8(d)分别为Pd-PGA、Pd-PGV预测模型的预测结果，可以看到Pd预测模型在小数值区域的预测值出现了明显高估，在大数值区域的预测值出现轻微低估。相比Pd预测模型，支持向量机预测模型在2002年9月16日Mj5.3级地震上没有明显的“小值高估、大值低估”现象，预测值更趋近于实测值。

(a)SVM-PGA

图9为支持向量机预测模型和Pd预测模型在2007年3月25日Mj6.9级地震的预测值与实测值的分布，1∶1的黑色实线表示预测值等于实测值，蓝色虚线表示基于训练数据集得到的1倍误差标准差。图9(a)和图9(c)分别为SVM-PGA与SVM-PGA预测模型的预测结果，可以看到预测值与实测值趋近1∶1线性比例关系，虽然大部分数据在1倍误差标准差范围内，但整体出现了轻微的低估现象。图9(b)和图9(d)分别为Pd-PGA、Pd-PGV预测模型的预测结 果，可以看到Pd预测模型预测PGA的结果在小数值区域出现了明显高估、大数值区域出现了明显低估，且PGV的预测结果整体出现了低估。虽然支持向量机预测模型也出现了整体轻微低估的现象，但是这是由于3 s的预测时间窗与大地震事件较长的震源破裂过程持续时间不匹配所造成的[31]，相比Pd预测模型，支持向量机预测模型在2007年3月25日Mj6.9级地震没有明显的“小值低估”现象，且能够明显改善“大值低估”现象。

(a)SVM-PGA

4 结 论

本文以0.1～10 Hz带通滤波后三分向矢量合成地震动峰值PGA与PGV为预测目标参数，利用日本K-net强震台网2007年10月-2019年8月期间P波触发后3 s数据，基于加速度幅值Pa、速度幅值Pv、位移幅值Pd、傅里叶谱幅值AMmax、速度平方积分IV2、破坏烈度DI、累积绝对速度CAV、阿里亚斯烈度Ia这8种特征参数构建支持向量机预测模型。得到的结论如下：

(1)无论是训练数据集、测试数据集还是两次随机震例数据集，支持向量机预测模型SVM-PGA与SVM-PGV的预测值的整体误差标准差均小于常规的Pd预测模型Pd-PGA与Pd-PGV；

(2)对比Pd预测模型Pd-PGA与Pd-PGV，支持向量机预测模型SVM-PGA与SVM-PGV的预测值与实测值更趋近于1∶1的线性比例关系，且预测值的误差不随震中距的变化而变化；

(3)相比较Pd预测模型Pd-PGA与Pd-PGV整体预测值明显的“小值高估、大值低估”现象，支持向量机预测模型SVM-PGA与SVM-PGV的预测值在6.5级以下地震没有明显的“小值高估”现象，在6.5级以上地震中由于预测时间窗与震源破裂时间不匹配的原因而出现轻微低估。

综上所述，本文构建的支持向量机预测模型可用于现地地震预警地震动峰值、即仪器地震烈度的预测。

致谢

日本防灾科学技术研究所(NIED)为本研究提供了数据支持，所有数据均从K-net网站下载(网址:http://www.kyosh-in.bosai.go.jp/),文中图件使用通用制图工具GMT(Genetic Mapping Tools)绘制。