基金项目:国家自然科学基金——基于竖向台阵记录数据的强震动作用下土非线性动力特征的实证研究(3092)、中国铁道科学研究院科研项目(2017YJ137)、云南省科技计划项目(2017 FD088)以及昆明学院校级科研项目(XJL15007)联合资助.
(1.云南省地震局,云南 昆明 650224; 2.昆明学院,云南 昆明 650214; 3.中国铁道科学研究院铁道科学技术研究发展中心,北京 100081; 4.中国地震局地球物理研究所,北京 100081)
(1.Yunnan Earthquake Agency,Kunming 650224,Yunnan,China )(2. Kunming Vniversity,Kunming 650214,Yunnan,China)(3. Railway Science and Technology Research and Development Centre,China Academy of Railway Sciences,Beijing 10081,China)(4. Insititute of Geophysics,China Earthquake Administration,Beijing 10081,China)
earthquake early warning; magnitude estimation; artificial neural networks; characteristic parameters; linear fitting
备注
基金项目:国家自然科学基金——基于竖向台阵记录数据的强震动作用下土非线性动力特征的实证研究(3092)、中国铁道科学研究院科研项目(2017YJ137)、云南省科技计划项目(2017 FD088)以及昆明学院校级科研项目(XJL15007)联合资助.
预警震级测定是地震预警的关键技术环节之一。在满足地震预警系统时效要求的前提下,以国内现有的人工神经元网络构架为基础,考虑采用更多的特征参数,对实时持续计算确定预警地震震级的方法进行研究。通过对日本部分实际强震数据进行持续估算预警震级与实际震级间的偏差情况,对预警震级和实际震级进行线性拟合,提出对预警震级结果的修正公式,进一步完善本方法快速估算预警震级的准确程度。
The magnitude determination is the critical technology of earthquake early warning(EEW). With satisfying the timeliness requirements of the EEW system,this paper focuses on the research of magnitude real-time continuous determination considering adopting more characteristic parameters based on the current Artificial Neural Networks(ANN). We continuously calculate the deviation between estimation magnitude and the real magnitude by using partial strong motion records from Japan,and then suggest the revised model for estimation magnitude of EEW after linear fitting of those two magnitudes to further improve the accuracy of this fast magnitude estimation method.
引言
自然灾害本身不可避免,但如果能预先采取一些合理的防御措施就可以有效地减少这些灾害造成的损失。为了达到减轻地震灾害的目的,除加强城市工程结构抗震设计外,人们最先想到的就是地震预报,但现有科技水平还无法彻底攻克这一难题,地震的预测预报必将长期处于探索和研究阶段。地震预警是目前世界上公认的能够有效减轻地震灾害的新手段之一(金星等,2012; 何少林,2017)。
实时震级计算是地震预警系统中最重要功能模块之一,也是整个地震预警系统中最复杂、最困难的部分(张红才,2013; 杨黎薇等,2017)。预警震级对时效性要求很高,主要利用布设在潜在震源区周围的实时传输地震观测台站,在破坏性地震发生后极短时限内,根据距离震中较近的若干个触发台站信息,迅速判断地震规模,并采用这若干个触发台站数据估算震级,随着触发台站数目不断增多,以不停变更的信息量对震级最初测定结果进行修正,在规定时限内得到最终测定结果。
1 预警震级估算相关研究
目前,国际上也发展形成了一些实用的实时震级测定方法,所采用的参数虽然各不相同,但基本突破点主要建立在有效利用P波段携带地震信息这一基础条件上。Nakamura(1988)最早提取P波段初始数秒内地震信息去评估地震震级大小。在此思路的影响下,地震学家通过对P波段以及S波段记录的深入分析,以实测地震记录为基础,得到了一些比较成熟的地震预警震级测定方法。最常见的计算方法大致可分为周期(频率)参数算法、幅值参数算法、能量参数算法以及其他算法。
1.1 周期(频率)参数算法大多数的地震预警系统使用的地震预警震级估算主要是以周期(频率)参数为依据推演得到的。其代表方法有τmaxp方法及τc方法。
τmaxp方法主要利用实时速度记录去计算地震动的卓越周期,τmaxp值是从台站触发开始的若干时间内(通常为3 s)计算得到的卓越周期最大值。多年来,各国专家学者不断对τmaxp参数法进行改进,其中Shieh 等(2008)所定义的τmaxp公式得到广泛认可:
τmaxp=max(2π((Xi)/(Di))1/2)(1)
Xi=αXi-1+x1i(2)
Di=αDi-1+(dx/dt)2i(3)
式中:α表示平滑参数,该值决定了平滑过程的速度,一般取值0.999; xi表示记录中地面运动的速度时程; Xi表示平滑后地面运动的速度导数平方值; Di表示平滑后地面运动的加速度导数平方值。
Kanamori(2005)对τmaxp改进后提出τc方法。其后,Shieh 等(2008)对τc参数的计算方式进行详细推导:
τc=(2π)/(([∫ t00u2(t)dt〖JB<2/〗∫ t00u2(t)dt])1/2)(4)
式中:积分区间[0,t0]表示记录中P波触发后3 s内的时间窗; u(t)是位移时程。τc 方法的基本思路与τmaxp方法是一脉相承的。τmaxp通过步步积分获取固定时间窗内周期参数; τc通过区间积分获取固定时间窗内周期参数。
周期参数与震级之间并不是简单的线性函数关系。τmaxp参数是幅值和频率的非线性函数,其准确性和稳定性受采样率影响,且与记录的预处理过程密切相关,采用不同滤波器或者不同长度时间窗,计算出的预警震级有明显差异; 而改进后的τc 参数直接将时间窗长度设置为3 s,对于6.5级以下的地震,基本可以根据τc 参数算法做出准确的预警震级估算。从目前的研究现状来看,使用周期参数估算预警震级对于中小震而言,估算结果是较为理想的。
1.2 幅值参数算法幅值参数算法是为了有效利用地震P波初始数秒内的波段信息而引入了记录波形峰值位移的幅值参数,其对位移幅值Pd 参数的定义是初始P波3 s时间窗内的垂直分量峰值位移。Wu等(2007)采用2阶高通巴特沃斯滤波器(低频截止频率为0.075 Hz)进行滤波,利用美国南加州地震记录,选用捡拾到P波后3 s时间窗内的位移幅值Pd的衰减关系去预测震级,以此为基础,张红才(2013)利用基本的震源理论推导出初始P、S波段的位移幅值Pd与最终震级间的关系,具体公式如下:
u(t)=const-1/RM·(t-R/c)=const1/RΔu·∑
=const1/RΔu·CL2(5)
式中:u(t)表示震中距为R处的P、S波位移场; const表示常数; R表示震中距; M·表示地震矩速率; c表示地震波速; Δu·表示断层平均滑动速率; Σ表示断层中初始阶段的滑动断层面积; C表示一阶几何参数; L表示断层线性尺度。
多数学者对幅值参数算法的稳定性及可靠性表示肯定,这种算法利用单个台站触发3 s时间内的信息记录就可得到相关参数,对台网密度的要求有了极大的降低,张红才(2013),金星等(2012)主要推荐将Pd 参数法作为国内地震预警系统中优先采用的方法。这种算法相对较简便,节约时间,但对于M≥6.5的地震,其预警结果误差较大。
1.3 能量参数算法能量参数算法是从能量角度去考虑预警震级的测定方法,它充分利用了累计绝对速度(CAV)作为强地面运动的快速捡拾量,最初是用于伊斯坦布尔的早期地震预警系统,用来决定是否有破坏性地震正在发生。其定义式为:
CAV=∫0∫tmax〖JB<1*|〗a(t)〖JB>1*|〗dt(6)
式中:积分下限0表示从台站触发开始计算,积分上限tmax可自行设定,即累计绝对速度(CAV)是通过对加速度积分计算得到的。该算法的基本思路是设定阈值,当给定台站的CAV超过设定阈值时,第一触发就会发生,有3个台站超过设定阈值触发时初始警报就会产生并对外宣布地震警报; 初始警报过后,系统将自动更新设置一个更高的阈值,同样,当3个台站超过新设置的阈值时第二次警报继续发布。能量参数主要是用于甄别地震是否具有破坏性,并不直接用于震级估算。
1.4 其他算法Odaka等(2003)提出从单个地震记录快速估算震中距与震级的新方法。为了定量分析不同的地震波形,采用Bt·exp(-At)这种简单的函数形式,时间t从P波到达开始,通过对波形包络线最初的部分进行最小二乘拟合,确定出拟合系数A,B,进而将A,B值作为快速估算震级大小的重要参数。用Δ表示震中距时,lgB与lgΔ成反比关系,这种关系适用于不同的地震震级,同时不受高频噪声的影响。
在震级、震源深度、震中距等因素影响下,地震波会在特定来源和观测环境下形成各不相同的包络波形,故本文考虑用1个直观的形式去展现这些包络波形。与后来P波与S波的最大振幅相比,初始P波最开始的振幅通常是非常小的,本文构建1个对数波形,先基线校正消除零点漂移或直流电(DC)组件的影响,再添加1个比标准噪声偏差更小的振幅以免除零振幅现象。Odaka等(2003)使用Kyoshin-Net(K-Net)网上下载的强震动记录,通过对数波形拟合,探寻初始地震波某部分(纵波到达后数秒内)不同的包络系统在形式上对应的地震震级和震中距,进而确定出拟合系数A与B。具体操作如下:(1)绘出竖向加速度记录包络波形,确定P波到时;(2)取P波到时后3 s内的包络波形,拟合Bt·exp(-At)中A、B值,通常采用最小二乘法。
Odaka拟合系数法主要运用于日本UrEDAS系统,根据该系统的运行经验,采用该方法也能够快速准确获取地震震级。可是参数B的取值需从大量的地震记录中统计获取,且参数B具有强烈的区域性,因而在我国实际应用并不广泛。
综上,实际的地震记录是非常复杂的,地震震级的确定涉及到震源过程、传播介质、场地条件、仪器性能等多个方面,并不是单一的周期参数、幅值参数或能量参数等就能准确稳定估算出来,这些参数只能在一定程度上反映地震的规模,不同参数估计得到的预警震级结果也可能存在着一些差异(林华伟等,2016)。因此,以国内较成熟的构建为基本,合理发展多种特征参数,也是非常值得研究的问题。本文以国内现有人工神经元网络构架为基础,除了选用常见的地震震级指示参数,考虑增加τmaxp、τc、Pd、CAV以及记录前半段拟合系数A、B等特征参数,对实时持续计算预警震级进行研究。通过实测拟定相应的预警震级修正公式,进一步完善快速估算预警震级的准确程度。
2 基于人工神经元网络选取特征参数
本文的预警震级测定是在实测地震记录基础上,将地震P波段或者S波段的前几秒记录获取特征参数与实际震级大小相联系,得到预警震级测定的经验拟合关系。由于不同特征参数之间属于非线性关系,采用人工神经元网络来确定预警震级。
2.1 选取人工神经元网络人脑神经元既有局部的计算和存储功能,又可通过联结构成统一体系。人工神经元网络(简称ANN)是采用物理可实现的系统去模仿人脑神经细胞结构与功能的一种信息处理系统。ANN由大量简单处理单元构成,具有巨量并行性、存储分布性、高度非线性、结构变化性及自组协调性等特点(丛爽,2003)。它最大的特点是仅仅借助样本数据,无需建立系统的数学模型,就可对系统实现由Rn空间(n为输入节点数)到Rm空间(m为输入节点数)的高度非线性映射。故而在结构分析中,可以直接使用人工神经网络模型实现结构系统输入参数与输出参数之间的非线性映射,无需建立系统的数学模型(毛健等,2011)。图1是最常见的人工神经元网络结构模型,对输入参数选取控制做了多项尝试,并对一些在线计算非线性控制算法结果进行对比研究,其利用人工神经元网络作为非线性过程模型,结合控制思想来组成控制器,在线寻找最优的ANN过程输入参数,经过多种选择比较,选择出适合于控制思想的输入量。
2.2 多特征参数选取国内常用震级指示参数主要以周期和幅值两方面的参数为主。从周期参数与震级的关系来看,P波的有效位移首脉冲累积宽度、P波位移脉冲有效上升时间以及峰值比Vmax/Amax与震级之间的拟合 关系较好; 而在幅值参数中,P波在某时间段内不同频率的加速度峰值、速度峰值、位移峰值则直接影响着震级的结果(马强,2008)。故而,本文在震级指示参数的基础上,综合考虑增加τmaxp,τc,Pd,CAV,拟合系数A与B等与震级大小关系密切的特征参数,以此提高预警震级的准确程度。
图3显示了具体参数的选取,在综合应用初始P波3 s内的多个震级指示参数、有效特征参数,并采用人工神经元网络中的ANN算法模型进行震级预测时,该网络模型的实例输入多参数分量,可连续估算出预警震级结果。
3 数据选取
以日本KiK-Net网(http://www.kik.basai. go.jp/kik/)下载的强震动记录为主,选取2011—2015年112个M5.0~7.0地震事件(记录数据3 500余条),震级统一使用M表示,挑选出波形记录完整、不同震中距的台站数据,将数据以震中距50 km内、100 km内以及大于100 km分类,具体如图4所示。
4 预警震级统计结果
6 讨论与结论
地震预警关键是首报的准确度,首报越准,越能保证地震预警的实际效益。本文在总结借鉴前人预警震级方法的基础上,引入国际上成熟的特征参数τmaxp,τc,Pd,CAV以及地震记录前半段拟合系数A,B,利用2011—2015年日本KiK-Net网112余组M5.0~7.0地震事件(记录数据3 500余条),提出1套基于人工神经网络和多种特征参数的预警震级估算法。通过持续估算预警震级与实际震级间的偏差,对预警震级和实际震级进行线性拟合,拟定对预警震级结果的修正公式,进一步完善快速估算预警震级的准确程度。
通过对研究结果的统计分析,得到如下结论:
(1)实际的地震记录是非常复杂的,单一的周期参数或者幅值参数估算预警震级时,其离散程度与估算结果准确性息息相关。采用人工神经网络方法持续对震级进行预测,应先检验选取的特征参数与地震震级间相关性是否紧密,即能否解决实质问题。本方法估算出的预警震级结果,其准确程度总体依赖前期工作中测试样本的可信性,故前期的训练极其关键。
(2)内陆地震震级多数为M7.0以内,以内陆地震为主的地震事件,预警震级的估算值通常比实际震级发震值大,可采用公式(8)进行相应的修正。海底地震震级多数在M7.0以上,以海底地震为主的地震事件,预警震级的估算值通常比实际震级发震值小,可采用公式(10)进行相应的修正。对所有M5.0~7.0破坏性地震而言,预警震级的估算值均可采用公式(12)进行相应修正。
(3)最后,本文所得到的结论与统计结果均是建立在日本台网地震事件的基础上,是否适用于中国地区还需更多国内的地震事件进行检验,有针对性完善并改进本文所提出的方法。
综上所述,本文所提出的预警震级估算法,以人工神经网络为基础,综合多个特征参数快速计算预警震级,并拟定了不同破坏性地震发生时可对应采用的修正公式,提高了预警系统震级信息的准确性,为地震预警系统建设提供一定的技术支持。
连续估算预警震级需要设置合理时限,最早估算出的预警震级称为最初预警震级,最后估算出的预警震级称为最终预警震级。以图4统计震例为主,选取主要特征参数进行人工神经网络强震预警震级估算,通过预警震级与实际震级间的偏差统计,寻求预警震级与实际震级间的拟合关系。
4.1 内陆地震中预警震级与实际震级间的差异为了探寻内陆地震与海底地震间预警震级的差异,本文将分别进行统计分析,将2者结果展开讨论。
图5 M5.0~6.0地震事件最初预警震级(a)、最终预警震级(b)与实际震级间的偏差
Fig.5 Deviation of original warning magnitude(a),final warning magnitude(b)and the actual magnitudes of M5.0~6.0 earthquakes图5主要展现了预警震级与实际震级间偏差变化趋势,在该震级区间内,2者偏差以正数为主,即预警震级估算值比实际震级发震值大。图5a中每1个数据点均代表了1个地震事件最初预警震级偏离实际震级的结果,震级相同但发震不同的地震,预警结果也各不相同。随着实际震级逐渐增大,最初预警震级与实际震级间偏离程度逐渐减小。图5b中每1个数据点均代表了1个地震事件最终预警震级偏离实际震级的结果,随着实际震级逐渐增大,最终预警震级与实际震级的偏离程度也逐渐减小。可以看出图5a中数据的离散程度大于图5b。为了寻求两者实质区别,本文计算了M5.0~6.0地震事件的预警震级与实际震级的平均偏差,如表1所示。
由表1可见,实际震级为M5.0时,预警震级与实际震级的平均偏差最大,基本接近0.8; 随着实际震级逐渐增长至M6.0,预警震级与实际震级的平均偏差越来越小,基本接近0.1。仅从表1可看出,最初预警震级与最终预警震级间差别非常小。
表1 M5.0~6.0地震事件预警震级与实际震级的平均偏差
Tab.1 Average deviation of estimated magnitude and actual magnitude of M5.0~6.0 events表2主要计算了 M5.0~6.0地震的预警震级与实际震级的平均标准偏差。表2的结果显示,最终预警震级与实际震级的平均标准差明显小于最初预警震级与实际震级的平均标准差。
表2 M5.0~6.0地震事件的预警震级与实际震级的平均标准偏差
Tab.2 Average standard deviations between estimate magnitude and the actual magnitude of the M5.0~6.0 events综上,最终预警震级的偏差相对更小,故而选择最终预警震级作为最后的预警震级结果,对其结果进行线性拟合,得到:
Mwarn-Mpra=-0.81649Mpra+4.91948(7)
Mwarn=0.18351Mpra+4.91948(8)
式中:Mwarn表示预警震级; Mpra表示实际震级,下同。
由此,对于内陆地区的强震预警,我们可尝试先以主要特征参数进行人工神经网络强震预警震级估算,再结合公式(8)对预警震级进行修正。
4.2 海底地震中预警震级与实际震级间的差异图6为M>6.0地震事件的预警震级与实际震级间的偏差变化,与图5的变化趋势基本一致,但由图6可看出,在该震级区间内,预警震级与实际震级间的偏差以负数为主,即预警震级估算值比实际震级发震值小。比较图6a和6b,两者的离散程度相差不大,偏差程度随着震级增大而逐渐增大(偏差值为负数)。
表3的计算结果显示,实际震级为 M6.1时,预警震级与实际震级的平均偏差为0,属于最佳理想结果; 随着实际震级逐渐增大,预警震级与实际震级的平均偏差越来越大,最大偏差值超过了-2.8。仅从表3可看出,最初预警震级与最终预警震级间差别接近于0。
图6 M>6.0地震事件预警震级与实际震级间的偏差
Fig.6 Deviation between estimated magnitude and the actual magnitudes of all the M>6.0 events表3 M>6.0地震事件的预警震级与实际震级的平均偏差
Tab.3 Average difference of the actual magnitude and testimate magnitude from M>6.0 earthquake events表4主要计算了M>6.0地震的预警震级与实际震级的平均标准偏差。表4结果中有一个明显划分,M6.0~7.0地震的最终预警震级与实际震级的平均标准差小于最初预警震级与实际震级的平均标准差,对比单独使用τc法和Pd法有明显改善; M>7.0地震的,最终预警震级与实际震级的平均标准差大于最初预警震级与实际震级的平均标准差。
表4 M>6.0地震事件的预警震级与实际震级的平均标准差
Tab.4 Average standard deviation of the estimate magnitude and actual magnitude from M>6.0 earthquake events综上,当实际震级 M>7.0地震,预警震级基本无效。本文对M6.0~7.0地震的最终预警震级结果进行线性拟合,得到:
Mwarn-Mpra=-0.91899Mpra+5.5943(9)
即:Mwarn=0.081011Mpra+5.5943(10)
由此,对于海底地区的强震预警,我们可尝试先以主要特征参数进行人工神经网络强震预警震级估算,再结合公式(10)对预警震级进行修正。
5.3 不同震级下预警震级与实际震级间的差异将所有M5.0~7.0地震事件的预警结果与实际结果进行对比,如图7所示,可看出M6.0地震事件预警震级与实际震级基本吻合。M<6.0地震事件预警震级大于实际震级; M>6.0地震事件预警震级小于实际震级(图7a)。取每一个震级所有事件预警结果的平均偏差值,其平均偏差随着实际震级增大呈线性变化(图7b)。
图7 预警震级与实际震级间偏差变化趋势
Fig.7 Difference variation tendency between the warning magnitude and actual magnitudes对所有破坏性地震的预警结果进行线性拟合,在M5.0~7.0地震的预警震级与实际震级间的关系为:
Mwarn-Mpra=-0.84268Mpra+4.95758(11)
即:Mwarn=0.17532Mpra+4.95758(12)
公式(12)可用于所有破坏性地震的预警震级修正,其准确程度还需以国内大量地震数据进行下一步的验证。
- 丛爽.2003.面向MATLAB工具箱的神经网络理论与应用[M].北京:中国科学技术出版社.
- 何少林.2017.地震烈度速度与预警台站选址相关问题探讨[J].地震研究,40(1):15-21.
- 金星,张红才,李军,等.2012.地震预警震级确定方法研究[J].地震学报,34(5):593-610.
- 林华伟,邱勇,刘超,等.2016.联合卓越周期的震级预测新方法[J].价值工程,29:226-227.
- 马强.2008.地震预警技术研究及应用[D].哈尔滨:中国地震局工程力学研究所,80-122.
- 毛健,赵红东,姚倩倩.2011.人工神经网络的发展及应用[J].电子设计工程,19(24):62-65.
- 吴逸民.2006.如何利用地震初达波从事地震预警[J].自然科学简讯,18(1):8-11.
- 杨黎薇,邱志刚,林国良.2017.强震预警中P波到时STA/LTA和贝叶斯BIC双步骤检拾研究[J].地震研究,40(4):629-637.
- 张红才.2013.地震预警系统关键技术研究[D].哈尔滨:中国地震局工程力学研究所,43-70.
- Kanamori H..2005.Real-Time Seismology and Earthquake Damage Mitigation[J].Annual Review of Earth and Planetary Sciences,33:195-214.
- Nakamura Y.1988.On the urgent earthquake detection and system[C].Proceedings of Ninth World Conference on Earthquake Engineering.Ⅶ:673-678.
- Odaka T,Ashiya K,Tsukada S,et al.2003.A New Method of Quickly Estimating Epicentral Distance and Magnitude from a Single Seismic Record[J].Bulletin of the Seismological Society of America,93(1):526-532
- Shieh J T,Wu Y M,Allen R M.2008.A comparison ofτc andτρ max for magnitude estimation in earthquake early warning[J].Geophysical Research Letters,35(20):1-5.
- Wu Y M,Kanamori H,Richard M A,et al 2007.Determination of earthquake early warning parameters,τc and Pd,for southern California[J].Geophys J Int,170:711-717.