2021, Vol. 64
2. 中国地震局地震工程与工程振动重点实验室, 哈尔滨 150080
2. Key Laboratory of Earthquake Engineering and Engineering Vibration of China Earthquake Administration, Harbin 150080, China
地震预警，指的是地震发生后对即将到来的破坏性地震动进行预测和警报(中国地震局，2015)，是近年发展起来的防震减灾有效手段之一(宋晋东，2013).世界上多个地震活跃国家或地区已运行或正在发展地震预警系统，如日本(Hoshiba et al., 2008)、墨西哥(Espinosa-Aranda et al., 1995)、美国(Kohler et al., 2017)、意大利(Picozzi et al., 2015a)、中国(Zhang et al., 2016；Peng et al., 2020)、中国台湾(Wu et al., 2013)等.
仪器地震烈度可为灾情快速判断、应急救援决策、工程抢险修复等提供科学依据(马强等，2014).快速准确预测地震发生后的潜在破坏，是地震预警的终极目标，仪器地震烈度的预测结果，也是地震预警系统效能的评判准则(Kodera et al., 2016).地震预警可分为区域预警与现地预警.区域预警基于地震波初期信息快速估计震源参数，再利用地震动衰减关系预测并判断地震可能造成的破坏及其影响范围(马强，2008)，但震源参数估计产生的偏差(Olson and Allen, 2005)、地震动衰减关系的不确定性(Iervolino et al., 2009)、目标区域的场地条件(Hoshiba et al., 2010)、大震级事件的破裂尺度(Hoshiba et al., 2011)等诸多因素都会造成目标场点地震动预测结果的误差.现地预警的基本概念是利用地震波初始阶段信息估计台站端地震波的后续地震动峰值，Kanamori(2015)认为现地预警可以比区域地震预警提供更多的预警时间，是地震预警未来的发展方向，现地地震预警中重要的基础工具就是基于地震波初始阶段参数预测后续地震动的预测方程，后续地震动预测的准确与否取决于预测方程的准确性.
早期的现地地震预警主要采用单个幅值参数预测后续的地震动.Nakamura(1998)基于加速度与速度的内积提出了一种能够表征地震破坏的参数DI，当DI时程的P波段峰值PI超过设定的阈值时认为地震将会产生现地破坏，该方法为现地地震预警利用P波段初始信息预测后续地震动峰值的雏形.Wu和Kanamori(2005a)定义0.075 Hz高通滤波记录在P波触发后3 s内的位移幅值为Pd(Peak displacement)，发现Pd与表征仪器地震烈度的速度峰值PGV(Peak Ground Velosity)呈现出极高的线性相关性，这就是地震预警中常用的现地预测Pd-PGV模型.刘辰等(2019a)研究了Pd-PGV模型的距离分段特征，同时，刘辰等(2019b)也提出了基于不同数据集建立Pd-PGV模型的改进方法.Brondi等(2015)利用地震破裂释放能量参数速度平方积分IV2(Integral of the Squared Velocity)预测PGV与豪斯纳烈度.王子珺和赵伯明(2016)也分别建立了P波触发后谱强度SI(Spectral Intensity)、速度平方积分IV2、位移幅值Pd与后续地震动的预测方程.
基于单参数的现地地震动预测方程只利用了地震触发后波形的单一信息、离散性较大，为提高准确性，基于多参数的现地地震动预测方程逐步引入到地震预警中.Wu和Kanamori(2005b)在Wu和Kanamori(2005a)在Pd-PGV模型的基础上，联合P波触发后3 s特征周期τc(characteristic period)判断现地台站的地震潜在破坏.Zollo等(2010)延续Wu和Kanamori(2005b)的概念，提出了“基于阈值的现地地震预警方法”，该方法可在P波触发后3 s快速估计潜在破坏区(PDZ，Potential Damage Zone)半径，Colombelli等(2012a, 2012b)、Carranza等(2013)、Picozzi等(2015b)、Peng等(2015)与宋晋东等(2017)分别验证了数次破坏性地震应用该方法的可行性.Kanamori(2015)分别建立了P波触发后3 s加速度幅值Pa(Peak acceleration)、速度幅值Pv(Peak velosity)、位移幅值Pd与PGV的线性关系作为预测方程，展现了不同幅值参数预测后续地震动结果的成功、误报、漏报的差异.Wu等(2011)通过Pd以及峰值加速度PGA(Peak Ground Acceleration)阈值联合判定地震的潜在破坏.Colombelli等(2015)还采用了基于Pa、Pv、Pd的等权重因子W(t)预测后续地震动PGV的新概念.可见，融合地震波初期多种特征参数有助于提升现地地震动的预测准确性.
除了准确性，后续地震动预测的连续性也需关注，因为地震是一个动态发展的过程，随着地震波的不断传播，地震预警的相关信息也应随时更新，后续的误报解除、警报更新等都需要预测方程实现连续不断的预测能力.Caruso等(2017)、宋晋东等(2018a, 2018b)拓展了时间窗概念，在P波触发后1 s、2 s、3 s分别建立预测方程实现了现地地震动连续预测.刘辰等(2019a)在P波触发后3~10 s时间窗、以1 s为时间间隔分别建立了Pd-PGV模型.王子珺和赵伯明(2016)也认为对于大地震事件，可以通过延长时间窗对目标场地进行较为准确的地震动预测.因此，需要基于更长时间窗以实现地震动连续预测能力，即也需要建立多时间窗条件下的预测方程.
因此，为提升现地地震动预测得准确性与连续性，本文以计算中国仪器地震烈度的0.1~10 Hz带通滤波后三分向矢量合成速度峰值PGV为预测目标参数，利用日本K-net和Kik-net强震台网P波触发后1~10 s时间窗的数据，基于人工智能领域的机器学习方法-最小二乘支持向量机(LSSVM，Least Squares Support Vector Machine)，选取加速度幅值Pa、速度幅值Pv、位移幅值Pd、速度平方积分IV2、破坏烈度DI、累积绝对速度CAV、阿里亚斯烈度Ia这7种特征参数作为输入，构建PGV连续预测模型LSSVM-PGV(Least Squares Support Vector Machine-based Peak Ground Velosity prediction model)，旨在为“国家地震烈度速报与预警工程”项目建设提供参考.
1 数据与处理选取日本K-net台网的土层台强震数据以及KiK-net台网的井上台强震数据，发震时间范围自2007年10月至2019年8月，发震区域为陆地或陆地沿岸近海，震级范围为气象厅震级Mj3.0~8.0，震源深度10 km以内，由于小震和震源距较大的地震信噪比较低，采用马强(2008)基于震源距的数据筛选方法，如式(1)，其中，R为震源距、Mj为气象厅震级：
|
(1) |
依据上述数据筛选原则，选取2069次地震事件，共15935组、47805条强震数据.其中，单独选取2016年4月的熊本地震序列的主震及5级以上前震与余震进行震例分析，剩余数据15469组、46407条随机划分为训练数据集与测试数据集两组，训练数据集与测试数据集的数据比例为7:3.图 1显示了本文所选取地震的震中及记录台站分布图.图中所用的地震震中、震级等震源参数来自于日本防灾科学技术研究所(NIED，National Research Institute for Earth Science and Disaster Prevention).震级与震源距的统计关系如图 2.
|
图 1 地震震中分布(a)与强震台站分布(b) 图a中空心圆，直径与震级大小成正比. Fig. 1 Locations of the epicenters (a) and stations (b) open circles, size in proportion to magnitude. |
|
图 2 强震数据的震源距及震级分布图 Fig. 2 Distribution of the number of strong ground motion records as a function of the magnitude and hypocentral distance |
对选取的强震数据进行如下处理：
(1) 使用马强等(2013)提出的长短时平均结合AIC(Akaike Information Criterion)准则方法自动识别P波震相以及到时，并对自动识别的地震P波到时进行人工校核；
(2) 加速度记录积分一次得到速度记录、积分两次得到位移记录，对积分后的记录进行4阶0.075 Hz高通滤波消除积分后的低频漂移现象；
(3) 依据《中国地震烈度表》(国家市场监督管理总局，国家标准化管理委员会，2020)，将三分向速度记录分别进行0.1~10 Hz带通滤波，对三分向速度记录矢量合成后计算得到速度峰值PGV.
2 方法 2.1 最小二乘支持向量机最小二乘支持向量机LSSVM(Least Squares Support Vector Machine)是通过在支持向量机SVM(Support Vector Machine)中引入最小二乘损失函数，将不等式约束转化为等式约束，是一种遵循结构风险最小化原则的核函数机器学习方法.由于LSSVM只需求解一组线性方程组，大幅降低了求解复杂度，因此构建模型时学习速度快，具有全局最优、适应性强等优点(余艳芳，2004).
设样本训练数据集{(x1, y1), (x2, y2), …(xi, yi), (xn, yn)}, xi, yi∈R, xi为输入向量，yi为输出向量，使用非线性映射φ(·)将样本输入从原空间映射到高维特征空间，则最小二乘支持向量机的回归函数表示为
|
(2) |
式(2)中：φ(x)为核空间映射函数；ω为权向量；β为常数.
依据结构风险最小化原则，最小二乘支持向量机求解二次规划问题可描述为
|
(3) |
式中, ei为松弛变量；γ为正则化常数.在式(3)中引入拉格朗日乘子α，构造拉格朗日函数为
|
(4) |
对ω, β, ei, αi求偏导为零可得：
|
(5) |
|
(6) |
|
(7) |
|
(8) |
将上述问题转化为求解线性方程组，即
|
(9) |
式中, E=(1, 1, …, 1)T；y=(y1, y2, …, yn)T；α=(α1, α2, …, αn)T；I为n阶单位矩阵；核函数K(xi, xj)=φT(xi)φ(xj).
求解方程(9)，得到参数α和β，则最小二乘支持向量机模型为
|
(10) |
在本文中，核函数采用机器学习中常用的高斯核函数(RBF，Radial Basis Function)，高斯核函数RBF将输入样本点映射成新的特征向量，然后得到该特征向量的点乘，因此，高斯核函数RBF的本质是将每一个样本点映射到一个无穷维特征空间，高斯核函数的数学表达式如式(11)：
|
(11) |
其中，θ为高斯核函数的带宽.
最小二乘支持向量机预测模型需要对计算参数进行设置，即正则化参数γ和核函数带宽θ，这两个参数决定了最小二乘支持向量机的学习和泛化能力.其中，正则化参数γ表示训练误差最小化与估计函数平滑性之间的权衡，核函数带宽θ表示核函数的附加参数.本文采用网格搜索方法，以21为间隔、从2-10到210范围内找出最优的正则化参数γ和核函数带宽θ组合值，同时为了防止过拟合，采取10折交叉验证方法，选取交叉验证根均方误差最小的参数构建最小二乘支持向量机模型.
2.2 预测目标参数《中国地震烈度表》给出了“国家地震烈度速报与预警工程”项目中仪器地震烈度计算的标准流程：获取地震动记录所用观测仪器应可为三方向，条件不具备时可为正交两水平方向或单水平方向；地震动加速度和速度记录的每个分向均应采用数字滤波器进行0.1~10 Hz带通滤波，PGA、PGV分别定义为滤波后加速度与速度矢量合成的最大值；利用PGA由式(12)得到的仪器地震烈度计算值IPGA，利用PGV由式(13)得到的仪器地震烈度计算值IPGV，进而可利用IPGA与IPGV依据式(14)计算出仪器地震烈度II.
|
(12) |
|
(13) |
|
(14) |
由《中国地震烈度表》规定的仪器地震烈度的计算流程可知，当IPGA与IPGV同时大于6.0时，仪器地震烈度II由IPGV计算得到，这表明在高烈度区、即仪器地震烈度大于Ⅵ度时，仪器地震烈度由PGV表征，同时，马强等(2014)的研究也表明，仪器地震烈度在高烈度区与PGV的相关性更高，因此，从预测强烈地震破坏的角度考虑，将本文的预测目标参数确定为：0.1~10 Hz带通滤波后三分向矢量合成速度峰值PGV.
2.3 输入特征参数选取以下7种特征参数作为最小二乘支持向量机预测模型的输入参数.
(1) 加速度幅值Pa：P波触发后计算时间窗内的垂直方向加速度最大值.
(2) 速度幅值Pv：P波触发后计算时间窗内的垂直方向速度最大值.
(3) 位移幅值Pd：P波触发后计算时间窗内的垂直方向位移最大值.
(4) 速度平方积分IV2(Integral of the Squared Velocity)：一个与地震早期辐射能量紧密相关参数，可用于现地地震动预测(Festa et al., 2008)，计算方法如下：
|
(15) |
其中，ν为垂直方向速度记录，tc为P波到时点，Δt为积分时间窗.
(5) 破坏烈度DI(Destructive Intensity)：通过计算垂直方向加速度a和速度v内积绝对值的对数来估计地震的危险程度(Nakamura，1998).DI的计算方法如式(16)：
|
(16) |
由式(16)可知，DI值表示“功率”的概念.
(6) 累积绝对速度CAV(Cumulative Absolute Velocity)：由地震动加速度绝对值对时间间隔积分得到(Reed et al., 1988)，CAV计算方法如下：
|
(17) |
(7) 阿里亚斯烈度Ia(Arias Intensity)：一个包含地震动幅值、频率及持时概念的地震动参数，由Arias(1970)提出，计算方法如式(18)：
|
(18) |
由于训练数据集是在对数坐标下得出预测模型，因此将对数预测值与对数实测值之差定义为误差xi，如式(19)：
|
(19) |
并计算误差标准差σ，如式(20)：
|
(20) |
式中，x为误差均值，N为误差数量.
表 1分别展示了P波触发后1~10 s范围及以1 s为间隔建立计算时间窗时，LSSVM-PGV模型基于训练数据集和测试数据集的对数预测PGV与对数实测PGV误差标准差统计结果.由表 1可知，训练数据集与测试数据集的误差标准差非常接近，在1 s时间窗相差0.008，其余时间窗最大相差0.005，这表明本文构建的LSSVM-PGV模型具备泛化性能，即在新鲜数据样本下模型的适应能力强，对具有同一规律的训练数据集以外的数据，经过训练的网络也能给出合适的输出.同时，随着时间窗的增长，误差标准差也呈现出逐渐减小的趋势.
|
表 1 P波触发后1~10 s时间窗LSSVM-PGV模型预测PGV的误差标准差 Table 1 The error standard deviation of the LSSVM-PGV model predicting the PGV in the 1~10 s time window after P wave arrivals |
图 3的离散点连线表示误差标准差随时间窗增长的变化规律，蓝色离散点表示基于训练数据集的误差标准差，红色离散点表示基于测试数据集的误差标准差.图 3中两条离散点连线的间隔展现了与表 1相同的结果，即训练数据集与测试数据集的误差标准差非常接近.同时可以看到，随着P波触发后时间窗的增长，训练数据集与测试数据集的误差标准差显著减小，6 s以后误差标准差的减小逐渐变缓，呈现逐步收敛的趋势.
|
图 3 P波触发后1~10 s时间窗LSSVM-PGV模型预测PGV的误差标准差随时间的变化 Fig. 3 The variation of the error standard deviation of the LSSVM-PGV model predicting the PGV in the 1~10 s time window after P wave arrivals |
图 4为P波触发后1~10 s范围，以1 s为间隔建立计算时间窗时，LSSVM-PGV模型基于测试数据集的预测PGV与实测PGV的比较，圆点表示每一个台站的预测结果，圆点的不同颜色表示震级的不同，黑色实线表示预测PGV与实测PGV的1:1线性比例关系，两条虚线表示±1倍标准差.
|
图 4 P波触发后1~10 s时间窗LSSVM-PGV模型在测试数据集上的预测结果 黑色实线为预测PGV等于实测PGV，黑色虚线为训练数据集的±1倍标准差. Fig. 4 The prediction results of the LSSVM-PGV model on the test data set in the 1~10 s time window after P wave arrivals The solid black line is the predicted PGV equal to the observed PGV, The black dotted line is ±1 standard deviation of the training data set. |
P波触发后1 s时(图 4a)，LSSVM-PGV模型的预测离散性较大，整体偏离预测PGV与实测PGV的1:1线性比例关系，且“小值高估”与“大值低估”现象明显.P波触发后2 s时(图 4b)，LSSVM-PGV模型预测PGV与实测PGV开始趋向1:1线性比例关系，特别是“大值低估”现象明显缓解.P波触发后3 s时(图 4c)至P波触发后5 s时(图 4e)，LSSVM-PGV模型预测PGV与实测PGV逐步趋向的1:1线性比例关系，“小值高估”与“大值低估”现象明显改善.从P波触发后6 s时(图 4f)开始一直到P波触发后10 s(图 4j)，LSSVM-PGV模型预测PGV与实测PGV整体符合的1:1线性比例关系，且变化细微，这与图 3显示的误差标准差变化趋势一致.
利用本文选取的训练数据集进行Pd-PGV模型建模训练，以考察常规的P波触发后3 s时间窗Pd-PGV模型与本文建立的P波触发后3 s时间窗LSSVM-PGV模型预测PGV的差异，Pd-PGV模型的建模结果如式(21)：
|
(21) |
由式(21)可知Pd-PGV模型预测PGV的误差标准差为0.304，大于P波触发后3 s时间窗下的LSSVM-PGV模型预测PGV误差标准差0.257.
图 5显示了P波触发后3 s时间窗及基于测试数据集时，LSSVM-PGV模型与Pd-PGV模型预测结果的比较.图 5a显示了LSSVM-PGV模型预测PGV与实测PGV的比较、图 5d显示了Pd-PGV模型预测PGV与实测PGV的比较，圆点表示每一个台站的预测结果，圆点的不同颜色表示震级的不同，黑色实线表示预测PGV与实测PGV的1:1线性比例关系，两条虚线表示预测PGV误差的±1倍标准差.可以发现，基于测试数据集时，LSSVM-PGV模型预测PGV的误差标准差为0.260、Pd-PGV模型预测PGV的误差标准差为0.304，LSSVM-PGV模型预测PGV的误差标准差小于Pd-PGV模型预测PGV的误差标准差.观察散点分布与1:1线性比例关系后发现，图 5a LSSVM-PGV模型的“小值高估”现象较图 5d Pd-PGV模型明显改善，同时，图 5d Pd-PGV模型的“大值低估”现象较图 5a LSSVM-PGV模型更为明显，特别是对于7级以上地震，Pd-PGV模型明显低估，而LSSVM-PGV模型只存在一定程度的低估.同时，就散点的整体分布而言，LSSVM-PGV模型较Pd-PGV模型更为贴近1:1线性比例关系.
|
图 5 P波触发后3 s时间窗LSSVM-PGV模型(a—c)与Pd-PGV模型(d—f)在测试数据集上预测PGV与实测PGV的比较 Fig. 5 The comparison of predicted PGV and measured PGV on the test data set between the LSSVM-PGV model and the Pd-PGV model |
图 5b与图 5e分别显示了P波触发后3 s时间窗时、基于测试数据集时，LSSVM-PGV模型与Pd-PGV模型预测PGV误差随震源距的变化，黑色实线表示误差为0，虚线表示±1倍误差标准差，点划线表示±3倍误差标准差.图 5b可以得出，LSSVM-PGV模型预测PGV误差在0~120 km震源距区间均匀分布于黑色实线两侧，整体误差位于±3倍标准差、即(-0.780，0.780)范围内.图 5e可以得出，Pd-PGV模型预测PGV误差在0~120 km震源距区间也均匀分布于黑色实线两侧，整体误差位于±3倍标准差、即(-0.912，0.912)范围内.这表明当P波触发后3 s时，本文构建的LSSVM-PGV模型预测PGV明显小于Pd-PGV模型.
考察预测PGV误差在震源距分布上的频数直方图及其对应概率密度分布后发现，LSSVM-PGV模型预测PGV的误差图 5c与Pd-PGV模型预测PGV的误差图 5f均呈现正态分布，均值都位于0附近，但LSSVM-PGV模型预测PGV的误差正态分布更显瘦高，Pd-PGV模型预测PGV的误差正态分布更显扁平，正态分布的形状特征表明了LSSVM-PGV模型预测PGV的误差明显小于Pd-PGV模型预测PGV的误差，这与图 5a与图 5d的结果一致.
4 熊本地震序列震例分析2016年4月16日1时25分(日本标准时间JST，Japan Standard Time)日本熊本县发生了Mj7.3地震，此次地震为熊本地震序列的主震，此外还发生了一系列破坏性前震及多次余震.熊本地震主震的最大地震烈度达到了日本气象厅烈度表规定的最高等级标准7度，地震还引发了大量建筑结构破坏、滑坡、砂土液化以及地表破裂破坏，造成了重大人员伤亡及财产损失(谢俊举等，2017).本文选取熊本地震序列主震、随机选取7次5级以上前震与余震作为震例数据库，且该数据库不在训练数据集与测试数据集中，以分析本文建立的LSSVM-PGV模型在具体地震事件中预测PGV的准确性.表 2为本文所用熊本地震序列震例数据库的地震事件信息，可以看到震例数据库中包含8次地震事件，最大震级事件是主震Mj7.3地震，还包含1次Mj6.4地震，以及6次Mj5以上地震、震级范围Mj5.0~Mj5.9.图 6显示了熊本地震震中和记录的台站分布.
|
|
表 2 本文选取的熊本地震序列事件 Table 2 Kumamoto earthquake sequence events selected in this study |
|
图 6 本文选取的熊本地震序列事件震中及台站分布 Fig. 6 Distribution of epicenters and stations of Kumamoto earthquake sequence events selected in this study |
将熊本地震序列震例数据库中的8次地震事件分为3组，即Mj7.3主震、Mj6.4地震、Mj6.0以下地震，分别考察LSSVM-PGV模型在这3组地震事件中的PGV预测结果.
图 7显示了P波触发后1~10 s时间窗下LSSVM-PGV模型在熊本地震序列Mj6.0以下地震数据组中预测PGV与实测PGV的比较，圆点表示每一个台站的预测结果，圆点的不同颜色表示震源距的不同，黑色实线表示预测PGV与实测PGV的1:1线性比例关系，两条虚线表示训练数据集下LSSVM-PGV模型预测PGV误差的±1倍标准差.从图 7a中可以看到，P波触发后1 s时，LSSVM-PGV模型的大部分预测PGV、特别是震源距0~25 km的近源台站预测PGV与实测PGV符合1:1线性比例关系，部分震源距25~75 km的近源台站出现低估.P波触发后2 s时，如图 7b，LSSVM-PGV模型预测PGV开始整体趋向1:1线性比例关系，从P波触发后3 s(图 7c)~10 s(图 7j)，LSSVM-PGV模型预测PGV整体符合1:1线性比例关系、且整体分布趋于一致，绝大部分结果分布在±1倍标准差范围内，没有出现“小值高估”与“大值低估”现象.
|
图 7 P波触发后1~10 s时间窗熊本地震序列Mj6.0以下地震数据组的预测PGV与实测PGV Fig. 7 Predicted PGV and observed PGV in the 1~10 s time window after P wave arrivals of seismic data set below Mj6.0 of the Kumamoto earthquake sequence |
图 8显示了P波触发后1~10 s时间窗下LSSVM-PGV模型在熊本地震序列Mj6.4地震数据组中预测PGV与实测PGV的比较，图中的标注与图 7一致.图 8a显示P波触发后1 s，LSSVM-PGV模型的大部分预测PGV与实测PGV符合1:1线性比例关系，少部分震源距0~50 km近源台站预测PGV出现低估.P波触发后2 s，如图 8b，LSSVM-PGV模型预测PGV开始整体趋向1:1线性比例关系，但仍有个别震源距0~25 km近源台站预测PGV出现低估.P波触发后3 s，如图 8c，LSSVM-PGV模型预测PGV开始整体符合1:1线性比例关系，近源台站也不再出现预测PGV低估.从P波触发后3 s(图 8c)~10 s(图 8j)，LSSVM-PGV模型预测PGV整体符合1:1线性比例关系且整体分布趋于一致，绝大部分结果分布在±1倍标准差范围内，没有出现“小值高估”与“大值低估”现象.
|
图 8 P波触发后1~10 s时间窗熊本地震序列Mj6.4级地震的预测PGV与实测PGV Fig. 8 Predicted PGV and observed PGV the 1~10 s time window after P wave arrivals of Mj6.4 of the Kumamoto earthquake sequence |
图 9显示了P波触发后1~10 s时间窗LSSVM-PGV模型在熊本地震序列主震Mj7.3地震数据组中预测PGV与实测PGV的比较，图中的标注与图 7一致.图 9a显示P波触发后1 s，LSSVM-PGV模型在所有震源距范围的预测PGV低于1:1线性比例关系，这表明此时预测PGV整体低估.P波触发后2 s(如图 9b)，预测PGV开始趋近1:1线性比例关系，直到P波触发后5 s(如图 9e)，预测PGV与实测PGV整体符合1:1线性比例关系，但0~50 km近源台站仍然部分低估、即“大值低估”现象.P波触发后6 s(如图 9e)，0~50 km近源台站预测PGV开始趋近1:1线性比例关系，“大值低估”现象明显缓解，且随着时间窗的增长，直到P波触发后10 s，如图 9j，预测PGV的整体分布趋于一致，绝大部分结果分布在±1倍标准差范围内，没有出现“小值高估”与“大值低估”现象.
|
图 9 P波触发后1~10 s时间窗熊本地震序列主震Mj7.3地震的预测PGV与实测PGV Fig. 9 Predicted PGV and observed PGV in the 1~10 s time window after P wave arrivals of mainshock Mj7.3 of the Kumamoto earthquake sequence |
以上熊本地震序列震例分析结果表明，对于6.0级以下地震，本文构建的LSSVM-PGV模型在P波触发后2 s即可以得到预测PGV整体符合实测PGV的结果，对于6.4级地震，也可以在P波触发后3 s即可以得到预测PGV整体符合实测PGV的结果，这与地震预警对于6.5级以下地震利用P波触发后3 s信息即可以确定震级并进行地震动预测的研究结果是一致的(Murphy and Nielsen, 2009).而对于6.5级以上地震，LSSVM-PGV模型至少需要P波触发后5 s时间窗才可以给出预测PGV整体符合实测PGV的结果.熊本地震序列主震破裂过程较为复杂，呈右旋走滑为主兼具部分正断分量的震源机制(Yagi et al., 2016；Asano and Iwata, 2016；USGS，2016).Kubo等(2016)利用强地面运动波形，基于多时窗线性波形反演得到熊本地震主震存在两次主要破裂，其中一个在破裂开始后4 s向东北方向的浅海区传播，并且直到16 s还有大的滑动.另一个在2~6 s从震源到地表破裂，6~10 s沿着近地表向东北方向破裂.郑傲等(2017)利用GNSS和InSAR数据，基于二段式断层模型、并结合活动构造资料作为参考构建熊本地震主震的有限断层模型，研究显示沿Futagawa-Hinagu断层带发生的右旋走滑破裂事件，发震位置分为南北两段，断层主要滑动集中在北段.Murphy和Nielsen(2009)认为对于大震级事件仅利用P波触发后3 s信息进行地震预警时，由于对断层破裂过程欠采样，因此无法反映整个破裂断层的规模，存在饱和震级6.5级.熊本地震序列主震破裂的复杂性决定了本文构建的LSSVM-PGV模型难以在P波触发的初期预测出趋近实测PGV的结果，但随着时间窗的增长，熊本地震序列主震的PGV预测随着时间窗增加而逐渐趋近于实测PGV，P波触发后6 s时间窗的预测PGV与实测PGV整体符合1:1线性比例关系，绝大部分结果分布在±1倍标准差范围内，没有出现“小值高估”与“大值低估”现象，这也与金星等(2012)的研究结论符合，即：如果利用更长的预测时间窗，将能够包含更多的信息，P波参数与地震动参数之间的相关性得到提高.
5 结论与讨论本文以计算中国仪器地震烈度的0.1~10 Hz带通滤波后三分向矢量合成速度峰值PGV为预测目标参数，利用日本K-net和KiK-net强震台网P波触发后1~10 s时间窗的数据，基于人工智能领域的机器学习方法-最小二乘支持向量机，选取加速度幅值Pa、速度幅值Pv、位移幅值Pd、速度平方积分IV2、破坏烈度DI、累积绝对速度CAV、阿里亚斯烈度Ia这7种特征参数作为输入，构建PGV连续预测模型LSSVM-PGV.结论如下：
(1) LSSVM-PGV模型在训练数据集与测试数据集上的预测误差标准差变化趋于一致，这表明本文构建的LSSVM-PGV模型具备泛化性能，即在新鲜数据样本下模型的适应能力强，对具有同一规律的训练数据集以外的数据，经过训练的网络也能给出合适的输出.
(2) LSSVM-PGV模型P波触发后3s的预测PGV与实测PGV即可整体符合1:1关系，随着时间窗的增长，PGV预测的误差标准差显著减小，并在P波触发后6s趋向收敛，这表明本文构建的LSSVM-PGV模型具备连续预测能力.对比常用Pd-PGV模型在P波触发后3 s时间窗的预测结果，LSSVM-PGV模型的PGV预测误差标准差明显减小，整体预测误差范围在0.780烈度单位，“小值高估”与“大值低估”现象明显改善，预测准确性相比较常用的Pd-PGV模型得到提升.
(3) 熊本地震序列的震例分析表明，对于6.5级以下地震，LSSVM-PGV模型最多在P波触发后3 s即可预测出与实测PGV整体符合1:1关系的PGV；对于7.3级主震，由于其破裂过程的复杂性，P波触发后3 s的预测结果出现一定程度的低估，但随着时间窗增长至6 s时，预测PGV与实测PGV符合1:1关系、并直到10 s整体趋势保持一致.
本文构建的最小二乘支持向量机PGV预测模型LSSVM-PGV具备连续预测能力，P波触发后3 s比较常规Pd-PGV模型准确性有所提升，可用于现地地震预警仪器地震烈度的预测.但本文的预测方程还是存在一定的离散性，特别是对于熊本地震7.3级主震这样的大震级事件，LSSVM-PGV预测模型在P波触发后3 s时间窗的预测结果依旧出现一定程度的低估现象，虽然通过增加时间窗可以获得较为准确的PGV预测结果，但是在P波触发后较短的时间窗内产出准确的仪器地震烈度预测结果依旧是未来研究的目标，选取与预测目标参数相关性更高的输入参数、引入深度学习方法以降低预测方程的离散性，或将LSSVM-PGV预测模型引入FinDer方法(Böse et al., 2012)以改善大震级事件的仪器地震烈度低估是今后值得研究方向.同时，虽然现地预警的基本概念是利用地震波初始阶段信息估计台站端地震波的后续地震动峰值，无需快速估算震级、没有利用衰减关系估计目标场地地震动、也没有考虑场地效应，但震源、路径、场地这三方面因素对现地预警方法暨LSSVM-PGV预测模型的预测结果影响也值得关注.
致谢 日本防灾科学技术研究所(NIED)为本研究提供了数据支持，所有数据均从NIED所属的K-net与KiK-net网站下载(网址:http://www.kyoshin.bosai.go.jp/)，文中图件使用通用制图工具GMT(Genetic Mapping Tools)绘制.
Arias A. 1970. A measure of earthquake intensity.//Hansen R J ed. Seismic Design for Nuclear Power Plants. Cambridge, Massachusetts: MIT Press, 438-483.
|
Asano K, Iwata T. 2016. Source rupture processes of the foreshock and mainshock in the 2016 Kumamoto earthquake sequence estimated from the kinematic waveform inversion of strong motion data. Earth, Planets and Space, 68(1): 147. DOI:10.1186/s40623-016-0519-9 |
Böse M, Heaton T H, Hauksson E. 2012. Real-time Finite fault Rupture Detector (FinDer) for large earthquakes. Geophysical Journal International, 191(2): 803-812. DOI:10.1111/j.1365-246X.2012.05657.x |
Brondi P, Picozzi M, Emolo A, et al. 2015. Predicting the macroseismic intensity from early radiated P wave energy for on-site earthquake early warning in Italy. Journal of Geophysical Research:Solid Earth, 120(10): 7174-7189. DOI:10.1002/2015JB012367 |
Carranza M, Buforn E, Colombelli S, et al. 2013. Earthquake early warning for southern Iberia:A P wave threshold-based approach. Geophysical Research Letters, 40(17): 4588-4593. DOI:10.1002/grl.50903 |
Caruso A, Colombelli S, Elia L, et al. 2017. An on-site alert level early warning system for Italy. Journal of Geophysical Research:Solid Earth, 122(3): 2106-2118. DOI:10.1002/2016JB013403 |
China Earthquake Administration. 2015. DB/T 59-2015 Technical requirements of instruments in network for earthquake monitoring-Seismic intensity instrument (in Chinese). Beijing: Seismological Press.
|
Colombelli S, Amoroso O, Zollo A, et al. 2012a. Test of a threshold-based earthquake early-warning method using Japanese data. Bulletin of the Seismological Society of America, 102(3): 1266-1275. DOI:10.1785/0120110149 |
Colombelli S, Zollo A, Festa G, et al. 2012b. Early magnitude and potential damage zone estimates for the great Mw9 Tohoku-Oki earthquake. Geophysical Research Letters, 39(22): L22306. DOI:10.1029/2012GL053923 |
Colombelli S, Caruso A, Zollo A, et al. 2015. A P wave-based, on-site method for earthquake early warning. Geophysical Research Letters, 42(5): 1390-1398. DOI:10.1002/2014GL063002 |
Espinosa-Aranda J M, Jiménez A, Ibarrola G, et al. 1995. Mexico City seismic alert system. Seismological Research Letters, 66(6): 42-53. DOI:10.1785/gssrl.66.6.42 |
Festa G, Zollo A, Lancieri M. 2008. Earthquake magnitude estimation from early radiated energy. Geophysical Research Letters, 35(22): L22307. DOI:10.1029/2008GL035576 |
Hoshiba M, Kamigaichi O, Saito M, et al. 2008. Earthquake early warning starts nationwide in Japan. Eos, Transactions American Geophysical Union, 89(8): 73-74. DOI:10.1029/2008EO080001 |
Hoshiba M, Ohtake K, Iwakiri K, et al. 2010. How precisely can we anticipate seismic intensities? A study of uncertainty of anticipated seismic intensities for the Earthquake Early Warning method in Japan. Earth, Planets and Space, 62(8): 611-620. DOI:10.5047/eps.2010.07.013 |
Hoshiba M, Iwakiri K, Hayashimoto N, et al. 2011. Outline of the 2011 off the Pacific coast of Tohoku Earthquake (Mw9.0)-Earthquake Early Warning and observed seismic intensity. Earth, Planets and Space, 63(7): 547-551. DOI:10.5047/eps.2011.05.031 |
Iervolino I, Giorgio M, Galasso C, et al. 2009. Uncertainty in early warning predictions of engineering ground motion parameters:What really matters. Geophysical Research Letters, 36(5): L00B06. DOI:10.1029/2008GL036644 |
Jin X, Zhang H C, Li J, et al. 2012. Research on earthquake early warning magnitude estimate. Acta Seismologica Sinica (in Chinese), 34(5): 593-610. DOI:10.3969/j.issn.0253-3782.2012.05.002 |
Kanamori H. 2015. Earthquake hazard mitigation and real-time warnings of tsunamis and earthquakes. Pure and Applied Geophysics, 172(9): 2335-2341. DOI:10.1007/s00024-014-0964-y |
Kodera Y, Saitou J, Hayashimoto N, et al. 2016. Earthquake early warning for the 2016 Kumamoto earthquake:Performance evaluation of the current system and the next-generation methods of the Japan Meteorological Agency. Earth, Planets and Space, 68: 202. DOI:10.1186/s40623-016-0567-1 |
Kohler M D, Cochran E S, Given D, et al. 2017. Earthquake early warning ShakeAlert system:West Coast wide production prototype. Seismological Research Letters, 89(1): 99-107. DOI:10.1785/0220170140 |
Kubo H, Suzuki W, Aoi S, et al. 2016. Source rupture processes of the 2016 Kumamoto, Japan, earthquakes estimated from strong-motion waveforms. Earth, Planets and Space, 68(1): 161. DOI:10.1186/s40623-016-0536-8 |
Liu C, Li X J, Jing B B, et al. 2019a. The distance segmentation characters of PGV-Pd relationship parameters for earthquake early warning. Chinese Journal of Geophysics (in Chinese), 62(4): 1413-1426. DOI:10.6038/cjg2019L0456 |
Liu C, Li X J, Jing B B, et al. 2019b. Improved method of estimating PGV using Pd for earthquake early warning. Journal of Seismological Research (in Chinese), 42(4): 608-615. |
Ma Q. 2008. Study and application on earthquake early warning[Ph. D. thesis] (in Chinese). Harbin: Institute of Engineering Mechanics, China Earthquake Administration.
|
Ma Q, Jin X, Li S Y, et al. 2013. Automatic P-arrival detection for earthquake early warning. Chinese Journal of Geophysics (in Chinese), 56(7): 2313-2321. DOI:10.6038/cjg20130718 |
Ma Q, Li S L, Li S Y, et al. 2014. On the correlation of ground motion parameters with seismic intensity. Earthquake Engineering and Engineering Dynamics (in Chinese), 34(4): 83-92. DOI:10.13197/j.eeev.2014.04.83.maq.011 |
Murphy S, Nielsen S. 2009. Estimating earthquake magnitude with early arrivals:a test using dynamic and kinematic models. Bulletin of the Seismological Society of America, 99(1): 1-23. DOI:10.1785/0120070246 |
Nakamura Y. 1998. A new concept for the earthquake vulnerability estimation and its application to the early warning system.//Proceeding of Early Warning Conference 98. Potsdam, Germany.
|
Olson E L, Allen R M. 2005. The deterministic nature of earthquake rupture. Nature, 438(7065): 212-215. DOI:10.1038/nature04214 |
Peng C Y, Yang J S, Chen Y, et al. 2015. Application of a threshold-based earthquake early warning method to the Mw 6.6 Lushan earthquake, Sichuan, China. Seismological Research Letters, 86(3): 841-847. DOI:10.1785/0220140053 |
Peng C Y, Ma Q, Jiang P, et al. 2020. Performance of a hybrid demonstration earthquake early warning system in the Sichuan-Yunnan border region. Seismological Research Letters, 91(2A): 835-846. DOI:10.1785/0220190101 |
Picozzi M, Zollo A, Brondi P, et al. 2015a. Exploring the feasibility of a nationwide earthquake early warning system in Italy. Journal of Geophysical Research:Solid Earth, 120(4): 2446-2465. DOI:10.1002/2014JB011669 |
Picozzi M, Colombelli S, Zollo A, et al. 2015b. A threshold-based earthquake early-warning system for offshore events in southern Iberia. Pure and Applied Geophysics, 172(9): 2467-2480. DOI:10.1007/s00024-014-1009-2 |
Reed J W, Anderson N, Chokshi N C, et al. 1988. A criterion for determining exceedance of the operating basis earthquake. EPRI Report NP-5930. Palo Alto, CA: Electric Power Research Institute, 1-28.
|
Song J D. 2013. Research on seismic ground motion indices for operation control and single station earthquake early warning applied for high-speed railway[Ph. D. thesis] (in Chinese). Harbin: Institute of Engineering Mechanics, China Earthquake Administration.
|
Song J D, Li S Y, Wang Y, et al. 2017. Application of a threshold-based earthquake early warning to Italy Mw6.2 earthquake on August 24, 2016. Earthquake Engineering and Engineering Dynamics (in Chinese), 37(6): 15-22. DOI:10.13197/j.eeev.2017.06.15.songjd.002 |
Song J D, Jiao C C, Li S Y, et al. 2018a. Prediction method of first-level earthquake warning for high speed railway based on two-parameter threshold of seismic P-wave. China Railway Science (in Chinese), 39(1): 138-144. DOI:10.3969/j.issn.1001-4632.2018.01.19 |
Song J D, Jiao C C, Li S Y, et al. 2018b. A predicting method for magnitude 1 earthquake alarm of high-speed railways based on seismic early radiated P-wave energy. Journal of Vibration and Shock (in Chinese), 37(19): 14-22, 38. DOI:10.13465/j.cnki.jvs.2018.19.003 |
State Administration for Market Regulation, Standardization Administration of China. 2020. GB/T 17742-2020 The Chinese seismic intensity scale (in Chinese). Beijing: China Standard Press.
|
USGS. 2016. Preliminary finite fault results for the Apr 15, 2016 Mw7.0 1 km WSW of Kumamoto-shi, Japan earthquake (Version 1).
|
Wang Z J, Zhao B M. 2016. Fast and accurate prediction method for track ground motion in earthquake early warning for high speed railway. China Railway Science (in Chinese), 37(6): 128-134. DOI:10.3969/j.issn.1001-4632.2016.06.17 |
Wu Y M, Kanamori H. 2005a. Experiment on an onsite early warning method for the Taiwan early warning system. Bulletin of the Seismological Society of America, 95(1): 347-353. DOI:10.1785/0120040097 |
Wu Y M, Kanamori H. 2005b. Rapid assessment of damage potential of earthquakes in Taiwan from the beginning of P waves. Bulletin of the Seismological Society of America, 95(3): 1181-1185. DOI:10.1785/0120040193 |
Wu Y M, Lin T L, Chao W A, et al. 2011. Faster short-distance earthquake early warning using continued monitoring of filtered vertical displacement:A case study for the 2010 Jiasian, Taiwan, earthquake. Bulletin of the Seismological Society of America, 101(2): 701-709. DOI:10.1785/0120100153 |
Wu Y M, Hsiao N C, Chin T L, et al. 2013. Earthquake early warning system in Taiwan.//Beer M, Kougioumtzoglou I A, Patelli E, et al. eds. Encyclopedia of Earthquake Engineering. Berlin: Springer, doi: 10.1007/978-3-642-36197-5_99-1.
|
Xie J J, Li X J, Wen Z P. 2017. Long period strong ground motion from the April 15, 2016 Kumamoto MW7.0 earthquake in Japan. Chinese Journal of Geophysics (in Chinese), 60(11): 4431-4446. DOI:10.6038/cjg20171129 |
Yagi Y, Okuwaki R, Enescu B, et al. 2016. Rupture process of the 2016 Kumamoto earthquake in relation to the thermal structure around Aso volcano. Earth, Planets and Space, 68(1): 118. DOI:10.1186/s40623-016-0492-3 |
Yu Y F. 2004. The classification approach based on least squares support vector machine and its applications[Master's thesis] (in Chinese). Shanghai: East China University of Science and Technology.
|
Zhang H C, Jin X, Wei Y X, et al. 2016. An earthquake early warning system in Fujian, China. Bulletin of the Seismological Society of America, 106(2): 755-765. DOI:10.1785/0120150143 |
Zheng A, Wang M F, Zhang W B. 2017. Source rupture process of the 2016 Kumamoto prefecture, Japan, earthquake derived from near-source strong-motion records. Chinese Journal of Geophysics (in Chinese), 60(5): 1713-1724. DOI:10.6038/cjg20170509 |
Zollo A, Amoroso O, Lancieri M, et al. 2010. A threshold-based earthquake early warning using dense accelerometer networks. Geophysical Journal International, 183(2): 963-974. DOI:10.1111/j.1365-246X.2010.04765.x |
金星, 张红才, 李军, 等. 2012. 地震预警震级确定方法研究. 地震学报, 34(5): 593-610. DOI:10.3969/j.issn.0253-3782.2012.05.002 |
国家市场监督管理总局, 国家标准化管理委员会. 2020. GB/T 17742-2020中国地震烈度表. 北京: 中国标准出版社.
|
刘辰, 李小军, 景冰冰, 等. 2019a. 地震预警PGV-Pd关系参数的距离分段特征. 地球物理学报, 62(4): 1413-1426. DOI:10.6038/cjg2019L0456 |
刘辰, 李小军, 景冰冰, 等. 2019b. 地震预警系统中利用Pd估测PGV的改进方法. 地震研究, 42(4): 608-615. |
马强. 2008.地震预警技术研究及应用[博士论文].哈尔滨: 中国地震局工程力学研究所.
|
马强, 金星, 李山有, 等. 2013. 用于地震预警的P波震相到时自动拾取. 地球物理学报, 56(7): 2313-2321. DOI:10.6038/cjg20130718 |
马强, 李水龙, 李山有, 等. 2014. 不同地震动参数与地震烈度的相关性分析. 地震工程与工程振动, 34(4): 83-92. DOI:10.13197/j.eeev.2014.04.83.maq.011 |
宋晋东. 2013.高速铁路运行控制用地震动参数及单台地震预警技术研究[博士论文].哈尔滨: 中国地震局工程力学研究所.
|
宋晋东, 李山有, 汪源, 等. 2017. 基于阈值的地震预警方法在2016年8月24日意大利Mw6.2级地震中的应用. 地震工程与工程振动, 37(6): 15-22. DOI:10.13197/j.eeev.2017.06.15.songjd.002 |
宋晋东, 教聪聪, 李山有, 等. 2018a. 基于地震P波双参数阈值的高速铁路Ⅰ级地震警报预测方法. 中国铁道科学, 39(1): 138-144. DOI:10.3969/j.issn.1001-4632.2018.01.19 |
宋晋东, 教聪聪, 李山有, 等. 2018b. 一种基于地震早期辐射P波能量的高速铁路Ⅰ级地震警报预测方法. 振动与冲击, 37(19): 14-22, 38. DOI:10.13465/j.cnki.jvs.2018.19.003 |
王子珺, 赵伯明. 2016. 高速铁路地震预警线路场地地震动快速准确预测方法. 中国铁道科学, 37(6): 128-134. DOI:10.3969/j.issn.1001-4632.2016.06.17 |
谢俊举, 李小军, 温增平. 2017. 日本熊本Mw7.0地震的长周期地震动. 地球物理学报, 60(11): 4431-4446. DOI:10.6038/cjg20171129 |
余艳芳. 2004.基于最小二乘支持向量机的分类方法及应用[硕士论文].上海: 华东理工大学.
|
郑傲, 王铭锋, 章文波. 2017. 利用近场强震记录反演2016年日本熊本县地震震源破裂过程. 地球物理学报, 60(5): 1713-1724. DOI:10.6038/cjg20170509 |
中国地震局. 2015. DB/T 59-2015地震观测仪器进网技术要求 地震烈度仪. 北京: 地震出版社.
|