基于LSTM网络的单台仪器地震烈度预测模型

李山有,王博睿,卢建旗*,王傲,张海峰,谢志南,陶冬旺

1 中国地震局工程力学研究所地震工程与工程振动重点实验室,哈尔滨 150080 2 地震灾害防治应急管理部重点实验室,哈尔滨 150080

0 引言

大震会造成建筑结构的大量破坏及人员伤亡.为了减轻地震灾害,科学家除了采取地震动参数区划(Zhang and Jin,2008)、建筑抗震设计(Li et al.,2008)以及减隔振(Nepal and Saitoh,2020)等震前地震灾害的预防措施以外,还提出了一种实时减灾手段——地震预警系统(Allen et al.,2009).该系统首先识别破裂初期产生的P波信息(马强等,2013; Lu et al.,2020; Wong et al.,2021),进而利用该信息对此次地地震发生的位置(金星等,2012; Satriano et al.,2008; Horiuchi,2005; Long et al.,2020)、地震的大小(彭朝勇等,2013; Peng et al.,2017; 李鸿杰等,2017; 朱景宝等,2022; Yang et al.,2021),以及地震可能造成破坏的范围进行快速估计,为处在可能遭受地震破坏的区域内居民或重大工程发送预警信息,避免次生灾害及人员伤亡.目前,地震预警系统已经在多个地震多发国家和地区开始了运行或试运行,如美国ShakeAlert预警系统(Cochran et al.,2018),日本JMA预警系统(Hoshiba et al.,2008),意大利PRESTo预警系统(Satriano et al.,2011),以及我国在建的国家地震烈度速报与预警系统等.地震预警系统主要有3种类型(Allen et al.,2009):(1)区域预警,区域预警在可能发生地震的区域内均匀布置观测台站.当区域内发生一次地震时,区域地震预警系统首先利用少数台站接收到的P波信息估算本次地震的震源参数.例如:利用台站接收到P波的到时信息进行地震定位,利用P波的峰值位移等参数进行震级估算.在此基础上,利用地震动衰减关系估计此次地震可能造成建筑物破坏的影响范围,并向该范围内的用户发布预警信息.(2)异地预警,异地预警将台网布置在已知发震断层附近,当地震发生时,利用布置在发震断层附近的台站观测数据估计震源参数,结合地震动衰减关系估计预警目标区域可能遭受的破坏,为目标城市或区域发布预警信息.例如墨西哥地震预警系统就是典型的异地预警(Aranda et al.,1995).然而,人们并非对未来潜在震源的位置了如指掌,大部分可能发生地震的位置都是未知的.随着地震观测台网建设成本的降低以及台网密度的增加,这种预警方式逐渐被区域地震预警所取代.(3)现地预警,现地预警在重大工程厂址布设台站,当台站接收到的地震动超过一定阈值时,重大工程根据阈值采取合适的紧急处置措施(Nakamura,1988; Carranza et al.,2013).显然,这种预警方法操作简单,建设成本低,便于管理和维护.但是,阈值预警能够为重大工程提供的预警时间非常有限.为了克服这一缺点,人们发现P波的峰值位移与地震动峰值之间存在一定的比例关系,在此基础上提出了基于P波峰值位移的现地地震动预测方法(Wu and Kanamori,2005a,b; 刘辰等,2019; 宋晋东等,2021; Caruso et al.,2017).相比阈值预警方式,这种方法能够为重大工程提供更长的紧急处置时间.

不论现地地震预警技术还是区域地震预警技术,均采用P波峰值位移等参数估计地震的最终大小或地震动的最终大小,存在大震低估及精度低的问题,从而可能导致漏报或误报(Kodera et al.,2018).为了克服单一预警方式存在的缺陷,人们提出了融合多种方法的预警手段(Wu et al.,2017; Minson et al.,2017).例如,2011年日本311地震以后,Kodera等针对区域预警系统存在的漏报问题,提出了一种不依赖于震源参数估计的PLUM(The Propagation of Local Undamped Motion)预警方法.Caruso等人提出了融合区域台网及单台现地预测的地震预警原型系统(Caruso et al.,2017).尽管融合方法在一定程度上能够解决漏报或误报的问题,但不能解决预测大震低估及精度低的问题.

大量相关研究表明,断层的破裂是一个复杂的时间过程(Lo et al.,2018; Wang et al.,2019; Liu et al.,2019; Tian et al.,2019),由于基于P波峰值及周期等参数的震级及地震动预测方法仅仅使用了P波的低维信息,没法考虑断层破裂过程蕴含的震源破裂过程信息.尤其P波的峰值位移、周期等参数在大震级段均出现了饱和(Trugman et al.,2019; Leyton et al.,2018),难以解决大震低估的问题.想要得到更好的估计效果,需要一种能够考虑时间变化过程的地震动估计方法.而长短时记忆神经网络(Long short-term memory,LSTM)则是一种能够自动学习蕴含在时间序列中的规律的一种网络结构,是目前解决这一问题最适合的一种手段.

因此,本文采用LSTM网络,以地震动的能量、幅值及周期参数的变化过程作为输入参数,以中国仪器地震烈度(GB/T 17742-2020)作为预测目标,利用日本K-NET台网记录的102次地震学习了能量变化过程与烈度的对应关系.结果表明,该方法能够学到能量变化过程与最终烈度的关系,可以为地震预警中影响场的快速估算提供参考.

1 数据集的建立及参数提取

机器学习是一种数据驱动算法,数据集以及特征参数的选取对机器学习结果的好坏有着非常重要的影响.

1.1 数据集

本文选取了日本K-NET所记录的震级在2～7.4级的191次地震,10457条三分量加速度记录,图1为K-NET数据台站以及震中分布图.选取102次地震作为训练集,89次地震作为测试集,从这两部分中的每次地震中随机抽取部分记录组成验证集,保证各数据集中的震级、震源距分布尽可能接近.最终训练集5103条记录,测试集3781条记录,验证集1573条记录(见图2).

图1 K-NET数据台站及震源分布

图2 震级与震源距直方图

1.2 参数提取

(1) 预测目标

烈度是衡量地震动大小及建筑物破坏的一个综合指标,常用于震害评估.为了能够在一次地震发生以后快速获得建筑物破坏情况以用于应急救援,人们通过地震动参数与宏观烈度之间的统计关系建立了仪器地震烈度,用于地震烈度速报及应急救援.在地震预警系统中,仪器地震烈度也是地震预警信息发布的决策参数.因此,本文选择了中国仪器地震烈度作为预测目标.

(2) 输入特征参数

无论用那种手段进行地震动峰值参数的预测,地震动峰值参数的物理空间分布规律是所有方法首先要考虑的物理背景.人们在地震动衰减规律的研究方面已经取得了大量的研究成果,对地震动空间分布规律也取得了公认的认识.例如,忽略地震动的场地效应、近场饱和效应以及震级饱和效应等一系列影响因素之外,地震动峰值加速度(PGA)、P波峰值位移(Pd)以及烈度(I)等参数随距离(R)的衰减满足最简单的衰减规律:

log10(PGA)=a0+a1M+a2log10(R)±σln(PGA),

(1)

其中,PGA为地震动峰值加速度(cm·s-2),R为台站到震源或震中的距离(km),M为震级,a0、a1、a2为回归系数,σln(PGA)为对数标准差.

log10(Pd)=b0+b1M+b2log10(R)±σln(Pd),

(2)

其中,Pd为P波峰值加速度(cm·s-2),b0、b1、b2为回归系数,σln(Pd)为对数标准差.

I=c0+c1M+c2log10(R)±σI,

(3)

其中,c0、c1、c2为回归系数,σI为标准差.

通过上述公式,我们就可以直接建立地震动峰值加速度或烈度与P波峰值位移之间的关系,进而利用P波峰值位移预测该台站处的峰值加速度或烈度:

log10(PGA)=(a0-b0)+(a1-b1)M+(a2-b2)

×log10(R)+log10(Pd),

(4)

I=(c0-b0)+(c1-b1)M+(c2-b2)log10(R)

+log10(Pd),

(5)

由公式(4)、(5)可知,无论峰值加速度还是烈度与P波峰值位移的关系中,都无法忽略震级与距离的影响.由于不同地震动参数衰减规律的不同,衰减公式中的系数如a0与b0、a1与b1、a2与b2等并不相等,无法消除震级项与距离项,仍然需要考虑震级与距离参数的影响.因此,在机器学习的入参数选择中,我们需要考虑哪些参数间接地蕴含了震级的信息和距离信息.在此基础上,本文提出了以下几个用于机器学习的输入特征参数.

首先,震级是一次地震发震断层破裂释放的能量大小的一种量度,地震的震级越大、破裂面积越大、释放的能量越多、能量的释放速率也越大、断层破裂的时间越长.因此,能量是间接体现地震震级大小的一个参数.本文选取了累积绝对速度(Cumulative absolute velocity,CAV)(Campbell and Bozorgnia,2012)作为一个台站接收到由震源辐射的能量的衡量指标(如公式(6)所示).同时,我们选取了CAV增长率作为衡量断层能量释放速率的指标(如公式(7)).除此之外,本文还选择了表征地震动能量大小的速度平方积分(Squared velocity integral,IV2) (Brondi et al.,2015)作为输入特征参数(如公式(8)).另外,大量研究表明,地震的震级越大,产生地震动的卓越周期越长、P波的峰值位移越大,并建立了大量利用卓越周期或P波峰值位移的震级估算方法(Allen and Kanamori,2003; Olson and Allen,2005; Ziv,2014; Hildyard and Rietbrock,2010; Peng et al.,2017).因此,本文也选择了卓越周期及峰值位移作为间接反映震级信息的输入特征参数,其中,本文卓越周期采用了Hildyard等人提出的卓越周期实时计算方法(Hildyard and Rietbrock,2010).

其次,考虑到现在地震预警系统在快速定位方面的技术已经比较成熟,地震预警台网的密度也在不断增加,利用断层破裂初期3～5个台站接收到P波到时后就可以得到合理的定位结果(卢建旗等,2021),这为单台地震动预测中引入距离信息提供了可能.因此,本文选取了震源距作为衡量能量空间衰减的特征参数.

最后,考虑到中国仪器地震烈度是由地震动加速度峰值和速度峰值计算得到,本文也选取了速度幅值、加速度幅值的绝对值作为输入特征参数,计算方法见公式(9)—(10).

假设x=[x1,x2,…,xn]表示一条加速度记录幅值的时间序列,则该记录的CAV的时间序列CAV=[cav1,cav2,…,cavn]可以表示为

(6)

其中,Δt是加速度记录的采样间隔 (s).

CAV增长率的时间序列CAVr=[cavr1,cavr2,…,cavri]可以表示为

cavri=(cavi+0.3/Δt-cavi)/(0.3),

(7)

其中,cavi表示i时刻的CAV值.

速度平方积分定义为

(8)

速度幅值时间序列Pa的取值方法为

(9)

速度幅值时间序列Pv的取值方法为

(10)

v(t)为竖向加速度记录的速度时程(cm·s-1),即对加速度时程进行一次积分后,再在进行4阶巴特沃斯0.075～25 Hz带通滤波得到速度时程,a(t)为竖向加速度记录经4阶巴特沃斯0.075～25 Hz带通滤波得到的加速度时程;T为时窗长度.

2 方法

LSTM神经网络属于循环神经网络的一个分支,是一种用于处理时间序列特征的网络模型.这种网络能够学习蕴含在时间序列中的特征,可以建立时间序列特征与目标输出之间的映射关系,在机器翻译、情感分析以及交通流量分析等领域得到了广泛的应用.为了解决传统循环神经网络在序列较长时容易出现梯度消失和梯度爆炸的问题,该模型引入了记忆单元和门控单元(Hochreiter and Schmidhuber,1997),可以实现对时间序列前段信息的记忆以及对冗余信息的遗忘,以避免梯度消失和爆炸的问题(Gers et al.,1999; Bengio et al.,1994).本文建立了一个用于仪器烈度预测的神经网络模型(见图3).该模型由3个LSTM层和1个全连接层组成,第一层由8个LSTM神经元组成,用于处理100个样点的时间序列特征.为了加强神经网络对前几秒地震动信息的学习能力,我们采用了不等间隔的时间序列.即,0.2～10 s采用0.2 s的样点间隔,10.5～20 s采用0.5 s的样点间隔,21～50 s采用1 s的采样间隔,共100个样点,序列总长度为50 s.

图3 用于预测烈度的LSTM神经网络结构图

本文LSTM神经元中采用的记忆单元及门控单元函数如公式(11)所示.采用均方误差作为损失函数(公式(12)),利用“NAdam”优化器(Ruder,2017) 进行优化.

(11)

其中,ft,it和ot分别表示遗忘门、输入门和输出门的门限函数;ct表示状态函数;ht表示隐含状态函数;⊙表示逐元相乘算子;W表示系数矩阵,其下标第一项表示当前状态,第二项表示所通过的门限;bf,bi,bo,和bg表示偏置向量;σ(*)=1/(1+e-*).

(12)

(13)

表1 部分时窗长度(WL)下评价指标统计表

(14)

(15)

其中,scorek表示第k个模型的得分.

图4a及附表1展示了所有模型在测试集上的均方损失以及平均预警时间,图4b为各个模型在测试集上的综合评分.通过模型的综合评分,最终确定的神经网络模型参数为:层数L=3,输出神经元个数U=8,批量大小B=100,优化器O为NAdam(Ruder,2017),最终确定的神经网络结构如图3所示.

图4 模型超参数优化图

损失函数是展现模型训练成功与否的重要信息,图5展示了本文模型在训练过程中在测试集的损失函数与验证集的损失函数随训练次数(Epoch)的变化曲线图.图中训练集与验证集的损失随着训练次数的增加而减小,而且训练集与验证集的残差水平一致,说明神经网络学到了数据中的普适性规律.

图5 MSE损失随训练次数变化曲线

3 结果

将通过训练集得到的神经网络模型应用于测试集数据,其结果如图6所示.考虑到烈度本身具有较大的不确定性,以预测误差为1度作为分类依据将测试集预测结果同观测烈度的比较结果分为三个区域,即准确区(绿色),高估(Overpredicted,OP)区(蓝色)以及低估(Underpredicted,UP)区(橙色).图中左上角为当前时间窗以及高估比例,右下角为预测标准差(STD)和低估比例.由图6可知,当序列长度T为3 s时,VI度以下烈度的估计较为准确,大于VI度的烈度出现了部分低估,烈度被高估的比例很低.随着序列长度的增加,高烈度低估的现象逐渐消失.由于序列的长度越长,序列中含有的信息越丰富,预测精度随着时间序列长度的增加而提高,说明神经网络学到了时间序列中蕴含的规律.

考虑到在烈度达到VI度后会造成建筑物的破坏,地震预警系统在预测烈度达到VI度及以上都会发送预警信息.为了评价本文模型提供预警信息的准确性,本文将预测结果定义为4类:

真阳性(TP):Iobs≥VI且Ipred≥VI,表示对大于VI度成功预测;

假阳性(FP):Iobs

真阴性(TN):Iobs

假阴性(FN):Iobs≥VI 且Ipred

由于地震记录中大于VI度的样本数量和小于VI度的数量严重不平衡,小于VI度的记录占比较大,而大于VI度的记录占比较小,简单使用大于VI度或小于VI度被成功识别的比例难以客观评价模型的精度.因此,我们用准确率(Accuracy)、漏报率(Miss alarm)和误报率(False alarm)直观有效地表现模型的预测能力.上述评价指标的定义为

(16)

准确率表示在所有样本中,大于等于VI度和小于VI度被成功预测的样本所占的比例.即,根据预测烈度发布的预警次数中有效预警所占的比例.漏报率表示在观测烈度大于等于VI度的样本中,预测烈度小于VI度的样本所占的比例.误报率表示在观测烈度小于VI度的样本中,预测烈度大于等于VI度的样本所占比例.

为客观评价本文模型的性能,我们没有采用训练集和验证集的数据来计算三个评价指标,而是采用测试集数据计算了三个评价指标随序列长度的变化关系曲线(如图7所示).从图中可见,随着序列长度的增加,准确率值迅速增加,最终趋近于1,漏报率和误报率最终趋近于0,说明模型的准确性较好.

图7 LSTM-I模型评估曲线

表1列出了部分时窗长度下LSTM-I模型的准确率、误报率以及漏报率.可见3 s时间窗时预测准确率为94.21%,误报率为1.25%,漏报率为46.78%,至10 s时间窗准确率上升至97.84%,误报率为1.14%,漏报率下降至17.6%.为了评价模型的泛化能力或是否出现了过拟合的情况,本文分别绘制了训练集、验证集以及测试集的烈度预测标准差随序列长度的变化曲线(如图8所示).由图8可知,训练集、验证集以及测试集的预测烈度标准差非常接近,说明模型没有出现过拟合,具有很好的泛化能力.

图8 标准差随时窗长度变化曲线

为了评价模型的时效性,本文计算了观测烈度到达时间与预测烈度到达时间差(Lead time),到时差越大说明能够为预警用户提供的预警时间越长(如图9所示).由图9可见,到时差随震源距的增加而增加.和理论P-S波到时差相比,二者较为接近.说明对于绝大部分记录,模型在接收到P波信息后就对大于VI度的记录提供了正确的预测.然而,本文的预测模型也出现了一些时差小于0的样本,说明预测时间滞后.

图9 LSTM-I训练模型的时效性

4 对比

为考察本文建立的LSTM-I模型的先进性,利用本文选取的训练集进行回归得到P波触发后3 s时间窗Pd-PGV与Pd-PGA模型,继而通过《中国地震烈度表》(GB/T-17742-2020)给出的中国仪器地震烈度与PGA、PGV的换算关系得到Pd与烈度的对应关系.P波触发后3 s的位移峰值Pd的取值方法如公式(17)所示:

(17)

其中,d(t)为竖向加速度记录的位移时程(cm),即对加速度时程进行二次积分后,再进行4阶巴特沃斯0.075～25 Hz的带通滤波得到位移时程,Pd为P波前3 s的峰值位移(cm).

基于训练集5147条地震动记录,用最小二乘法回归得到的Pd与PGV和PGA的关系(如公式(18)、(19),图10所示),图10中黑色实线为最小二乘法回归结果,黑色虚线为±1倍标准差.

图10 PGV与PGA关于Pd的回归结果

log10(PGV)=0.823×log10(Pd)+1.757±0.427,

(18)

log10(PGA)=0.627×log10(Pd)+2.476±0.416,

(19)

图11为利用测试集数据基于Pd方法计算得到的P波触发后3 s时间窗的预测烈度和LSTM-I模型预测结果与实测烈度对比结果.图11a为LSTM-I模型在P波触发后3 s时间窗的预测结果与实际烈度的比较,图11c为Pd预测结果与实测烈度的比较.对比a与c两图可见,基于Pd方法计算得到的预测结果高估率为20.93%,低估率为7.13%.而LSTM-I模型预测结果的高估率为3.52%,低估率为7.45%.说明LSTM-I模型的 “小值高估,大值低估”现象有明显改善.

图11 P波触发3 s时间窗LSTM-I模型(a,b)与基于Pd模型计算预测(c,d)在测试集上的预测结果与实测烈度的比较

图11b与图11d分别展示了P波触发后3 s时间窗时LSTM-I模型与Pd方法预测的误差随震源距的变化,以及误差在震源距分布上的频数直方图和对应的概率密度分布曲线.左侧散点图内黑色实线表示误差为0,虚线和点划线分别代表±1倍标准差和±3倍标准差,右侧灰色曲线表示概率密度分布.从图11b可看出,LSTM-I模型预测烈度在0～500 km震源距区间内的误差均匀分布于均值-0.1两侧,整体误差位于±3倍标准差(-2.64,2.64)内.图11d可看出,基于Pd方法得到的预测结果在0～500 km震源距区间内的误差均匀分布于均值0.42两侧,整体误差位于±3倍标准差(-3.57,3.57)内.继而考察预测误差在震源距分布上的频数直方图以及对应的概率密度分布曲线,两方法对应的误差均呈正态分布,且均值均在0附近,但LSTM-I模型的标准差小于Pd模型的标准差.

5 结论与讨论

本文提出了一种基于时间序列特征的单台仪器地震烈度持续预测LSTM神经网络模型.该模型以累计绝对速度、累计绝对速度增长率、地震波卓越周期、加速度幅值、速度幅值、速度平方积分以及震源距为输入参数,以单台最大观测仪器烈度为预测目标,利用日本K-NET强震台网记录的102次地震数据训练了LSTM-I模型,用89次地震的数据作为测试集检验了模型的泛化能力.模型在训练数据集与测试数据集上的预测误差标准差变化趋于一致,这表明本文构建LSTM-I模型具备较好的泛化性.

从准确度分析,研究结果表明,P波触发后3 s,预测烈度出现部分高烈度低估,很少出现低烈度高估现象.随着时窗长度的增加,预测烈度的精度不断提高.相较于基于常规Pd-PGV和Pd-PGA模型的计算方法准确性有所提升,可用于现地预警仪器地震烈度预测.

从时效性分析,本文计算了每条记录VI度到达的时间与预测VI度到达的时间差,在不考虑系统延时的情况下可以看做该模型为用户能够提供的有效预警时间.结果表明,本文模型能够提供的预警时间与P-S波到时差接近,说明模型具有较高的时效性.

由于烈度本身存在一定离散性,LSTM-I模型在P波触发后3 s时窗的预测结果仍存在部分“大值低估,小值高估”的现象.虽随着时间窗长度增加,低估和高估的现象大幅度减少,可得到较为准确的预测结果.

致谢使用开源平台 Anaconda、Python、TensorFlow和Keras来设计神经网络.日本防灾科学技术研究所为本研究提供数据支持.

附表1 各个模型的训练损失与平均预警时间