CN113009566A - 现地地震动预测模型及其构建方法 - Google Patents

Abstract

本发明涉及一种现地地震动预测模型及其构建方法，所述构建方法包括：获取各强震动观测台站的加速度记录数据，对数据进行预处理，并筛选出标准数据；对筛选出的各强震动观测台站事件记录数据，去掉线性趋势，手动选取P波和S波初至时间的记录数据；对记录数据中的各强震动观测台站的垂直向加速度记录进行一次积分和两次积分以分别获取速度时程和位移时程；对位移时程进行连续的Butterworth带通滤波器滤波处理，以获取P波幅值预警参数P_vall和P_aall；按照仪器烈度计算标准要求计算得峰值速度PGV和峰值加速度PGA；对于P波幅值预警参数P_vall和P_aall，进行相关性拟合，得到P_vall与PGV、P_aall与PGA的相关性模型。该模型提高了预测结果的可靠性。

Description

现地地震动预测模型及其构建方法

技术领域

本发明属于地震预警领域，具体为一种现地地震动预测模型及其构建方法。

背景技术

在利用P波预警参数估算现地地震动时，需要提前构建现地地震动预测模型。传统一般使用P波3秒位移幅值参数P_d与峰值速度PGV的相关性来预估现地地震动，而采用其它类别的P波振幅参数(如速度幅值P_v或加速度幅值P_a)几乎没有。在实际应用中，对于强震仪记录到的加速度记录需要两次积分才能获得位移记录，而积分操作会引起低频漂移，对于经过两次积分获得的位移记录的影响则更大。另一方面，地震烈度仪采用MEMS传感器采集地面加速度，相对于传统力平衡加速度计来说，MEMS记录到的数据质量会差许多，特别是对于动态范围较低(90dB)的简易烈度计。在这种情况下，如果还是使用P_d与PGV的相关性完成现地地震动的估算，结果的可靠性将会大打折扣。

发明内容

本发明的目的在于提出一种现地地震动预测模型及其构建方法，本发明拟构建基于不同P波幅值参数的现地地震动预测模型，提高预测结果的可靠性。

为实现上述目的，本发明所采用的技术方案是一种现地地震动预测模型，所述模型为采用Butterworth带通滤波器处理数据拟合获得的速度幅值P_v或加速度幅值P_a与峰值速度PGV或峰值加速度PGA的相关性模型。

进一步的，所述模型包括采用全P波窗速度幅值P_vall与峰值速度PGV的相关性模型，所述P_vall与PGV的相关性模型采用1阶时的Butterworth带通滤波器对数据进行滤波处理获得，其中所述带滤波器的频带范围为0.1～10Hz。

更进一步的，所述模型还包括采用全P波窗加速度幅值P_aall与峰值加速度PGA的相关性模型，所述P_aall与PGA的相关性模型为采用1阶时的Butterworth带通滤波器对数据进行滤波处理获得，其中所述带滤波器的频带范围为0.1～10Hz。

其中，所述相关性模型的均按照式(1)进行相关性拟合而得，所述方程式为：

logPGx＝AlogPx+B (1)

式中，PGx代表PGV或PGA，Px表示P_v或P_a，而A和B为需要拟合获得的系数。

本发明第二方面公开了一种现地地震动预测模型的构建方法，至少包括下述步骤：

获取各强震动观测台站的加速度记录数据，对数据进行预处理，并筛选出标准数据；

对筛选出的各强震动观测台站事件记录数据，去掉线性趋势，手动选取P波和S波初至时间的记录数据；

对记录数据中的各强震动观测台站的垂直向加速度记录进行一次积分和两次积分以分别获取速度时程和位移时程；

对位移时程进行连续的Butterworth带通滤波器滤波处理，以获取P波幅值预警参数P_vall和P_aall；

按照仪器烈度计算标准要求计算得峰值速度PGV和峰值加速度PGA；

对于P波幅值预警参数P_vall和P_aall，按照式(1)进行相关性拟合，得到P_vall与PGV、P_aall与PGA的相关性模型。

一个优选的技术方案，所述对数据进行预处理，并筛选出标准数据具体包括：

对获取的加速度记录数据按照“发震时刻+震级+台站名”的方式进行重新命名，然后编写程序将所有的数据格式统一转换为文本文件格式；

对获取到的各类数据采用如下条件进行筛选：垂直向>5gal；震源距<200km；震级M≥4；同一地震事件至少有3个台站记录到；同时删除在高频能量到来时产生突跳的加速度记录。

一个优选的技术方案，所述对位移时程进行连续Butterworth带通滤波器滤波处理获取P波幅值预警参数P_vall和P_aall的具体方法包括：

对位移时程进行连续Butterworth带通波器滤波处理以移除由于积分操作带来的低频漂移影响，获取P波到后3s内垂直向数据对应的P_d3和全P波窗垂直向数据对应的P_dall；其中，带通滤波器的频带范围为0.075～3Hz；

按照仪器烈度计算标准针对于PGV和度PGA获取所用的Butterworth带通滤波器频带0.1～10Hz，采用1阶Butterworth滤波器对P波3秒和全P波段垂直向速度和加速度时程进行处理，以获取P_vall和P_aall。

进一步的，所述按照仪器烈度计算标准要求计算得3分向合成PGV和PGA。

第三方面本发明还公开了一种应用上述现地地震动预测模型的构建方法的构建系统，所述系统包括：

数据预处理模块，用于对获取的各强震动观测台站的加速度记录数据进行预处理，并筛选出标准数据；

参数计算模块，用于对筛选出的各强震动观测台站事件记录数据，去掉线性趋势，手动选取P波和S波初至时间的记录数据；对记录数据中的各强震动观测台站的垂直向加速度记录进行一次积分和两次积分以分别获取速度时程和位移时程；以及，对位移时程进行连续的Butterworth带通滤波器滤波处理，以获取P波幅值预警参数P_vall和P_aall，同时按照仪器烈度计算标准要求计算得PGV和PGA；

相关性模型建立模块，用于根据对于各类别P波幅值预警参数P_vall和P_aall，统一按照式(1)进行相关性拟合，得到P_vall与PGV、P_aall与PGA的相关性模型。

由于MEMS记录数据质量较差，在两次积分后容易引入较大的噪声干扰，使得利用P_dall与PGV关系式无法得到可靠的结果。而本发明P_vall与PGV、P_aall与PGA这两种模型更加可靠，这两种关系式不仅相关性优于其它关系式，而且仅需要采用低阶(1阶)滤波器进行处理即可，可大大减少计算量，有助于降低处理器负载，特别是对于设备端嵌入式计算(边缘计算)。

附图说明

为了更清楚地说明本发明实施例中的技术方案，下面将对实施例描述中所需要使用的附图作简要介绍，显而易见地，下面描述中的附图仅仅是本发明的一些实施例，对于本领域的普通技术人员来讲，还可根据这些附图获得其他的附图。

图1为现地地震动预测模型的构建方法；

图2为现地地震动预测模型的构建系统。

具体实施方式

为了使本申请的目的、技术方案及有益效果更加清楚明白，以下结合附图及实施例，对本申请进一步详细说明。应当理解，此处所描述的具体实施例仅仅用以解释本申请，并不用于限定本申请。

地震预警数据处理中，通常使用P波3秒时间窗，对于离震中很近的台站(R<10km)，P波3秒数据已包含S波信息，导致像P_d这样的幅值参数出现较大的变化，从而对现地地震动估测结果造成比较大的偏差。因此，本发明采用了不一样的处理方式，不仅使用P波3秒数据进行关系式拟合，还选用整个P波段数据拟合参数关系式。这里的P波段数据定义为从P波触发至S波到达之间的数据段。

下面参考附图对本发明现地地震动预测模型及其构建方法的实施例进行描述。

本发明所述模型为采用Butterworth带通滤波器处理数据拟合获得的速度幅值P_v或加速度幅值P_a与峰值速度PGV或峰值加速度PGA的相关性模型。特别是，采用全P波窗速度幅值P_vall与峰值速度PGV的相关性模型以及采用全P波窗加速度幅值P_aall与峰值加速度PGA的相关性模型；这两个相关性模型均采用1阶时的Butterworth带通滤波器对数据进行滤波处理获得，其中所述带滤波器的频带范围为0.1～10Hz。

附图1示出了本发明现地地震动预测模型的构建方法，所述方法包括下述步骤：

1)获取各强震动观测台站的加速度记录数据，对数据进行预处理，并筛选出标准数据；

在获取到数据后，由于数据文件的命名和格式等都不一致，首先对数据按照“发震时刻+震级+台站名”的方式进行了重新命名，然后编写程序将所有的数据格式统一转换为文本文件格式，具体为只包括实际记录数据，不包括头信息，并将单位统一转换为cm/s²，也就是gal。在充分考虑拟合所用数据质量、估算结果发散性和震级精度的情况下，同时参考Wu和Zhao(2006)的记录挑选标准，对获取到的各类数据采用如下条件进行筛选：①垂直向>5gal；②震源距<200km；③M≥4；④同一地震事件至少有3个台站记录到。此外，在人工查看数据时，发现个别的加速度记录由传感器的弹簧系统或者电子器件引起其在高频能量到来时会产生突跳，需要将该部分记录排除掉。

2)对筛选出的各强震动观测台站事件记录数据，去掉线性趋势，手动选取P波和S波初至时间的记录数据；

3)对记录数据中的各强震动观测台站的垂直向加速度记录进行一次积分和两次积分以分别获取速度时程和位移时程；

4)对位移时程进行连续的Butterworth带通滤波器滤波处理，以获取P波幅值预警参数P_vall和P_aall；

具体操作为：对位移时程进行连续Butterworth带通波器滤波处理以移除由于积分操作带来的低频漂移影响，获取P波到后3s内垂直向数据对应的P_d3和全P波窗垂直向数据对应的P_dall；其中，带通滤波器的频带范围为0.075～3Hz；

按照仪器烈度计算标准针对于PGV和度PGA获取所用的Butterworth带通滤波器频带0.1～10Hz，采用1阶Butterworth滤波器对P波3秒和全P波段垂直向速度和加速度时程进行处理，以获取P_vall和P_aall。

5)按照仪器烈度计算标准要求计算得峰值速度PGV和峰值加速度PGA，采用3分向合成值计算，具体计算步骤参考国家标准《GB/T17742-2020中国地震烈度表》中的附录A；

6)对于P波幅值预警参数P_vall和P_aall，按照式(1)进行相关性拟合，得到P_vall与PGV、P_aall与PGA的相关性模型。

本发明还提供了一种针对上述方法应用的构建系统，图2所示所述系统包括：

数据预处理模块，用于对获取的各强震动观测台站的加速度记录数据进行预处理，并筛选出标准数据；

参数计算模块，用于对筛选出的各强震动观测台站事件记录数据，去掉线性趋势，手动选取P波和S波初至时间的记录数据；对记录数据中的各强震动观测台站的垂直向加速度记录进行一次积分和两次积分以分别获取速度时程和位移时程；以及，对位移时程进行连续的Butterworth带通滤波器滤波处理，以获取P波幅值预警参数P_vall和P_aall，同时按照仪器烈度计算标准要求计算得PGV和PGA；

相关性模型建立模块，用于根据对于各类别P波幅值预警参数P_vall和P_aall，统一按照式(1)进行相关性拟合，得到P_vall与PGV、P_aall与PGA的相关性模型。

下面就上述方法代入现有地震数据进一步论证P_vall与PGV、P_aall与PGA的相关性模型为最优模型。

1.1所用震例

采用来源于国家强震动台网中心2007年1月至2014年12月共计23365条加速度记录数据。其中，2014年的加速度记录数据只有地震应急数据，包括2014年新疆于田M_S7.3地震、2014年云南鲁甸M_S6.5地震、2014年云南景谷M_S6.6地震、2014年四川康定M_S6.3和M_S5.8地震，以及2014年云南景谷M_S5.8地震。由于国内地震事件记录介于M_S6.5-8.0之间的数据非常匮乏，为了确保地震事件数据的均匀性和提高震级区间事件数据的完备性，利用日本强震动观测台网(K-net和KiK-net)的38次M_J6.5-8.0地震事件数据(共计3642条记录)进行补充。根据Tsuboi(1954)和Bormann et al.(2007)的研究，M_S与M_J两种震级之间的相差很微弱，可以忽略，将所有震级统一为M。

1.2数据预处理

按照上述步骤1)的方法最终筛选出2764个台站加速度记录，对应于285次地震事件。

1.3 P波幅值预警参数计算

参照步骤2)-5)的方法对P波幅值预警参数计算，其中，为了探讨不同阶数Butterworth滤波器对拟合参数关系式的影响，这里同时采用1～4阶四种Butterworth滤波器进行处理，也就是每类预警参数计算获得4个值。同理，参考仪器烈度计算标准针对于PGV和峰值加速度PGA获取所用的Butterworth带通滤波器频带(0.1～10Hz)，采用1～4阶四种Butterworth滤波器对P波3秒和全P波段垂直向速度和加速度时程进行处理，以获取P_v3、P_vall、P_a3和P_aall。

1.4统计关系式拟合

对于各类别P波幅值预警参数，统一按照式(1)进行相关性拟合，具体包括1-4阶Butterworth带通滤波器拟合条件下获取的P_d与PGV、P_d与PGA、P_v与PGV、P_v与PGA、P_a与PGV、P_a与PGA的相关性模型。

拟合过程根据离散数据拟合曲线，拟合关系结果如表1所示。

表1 P_d与PGV、P_d与PGA、P_v与PGV、P_v与PGA、P_a与PGV、P_a与PGA的相关性

针对上表1分析可知各参数关系式的相关性状况，举例说明，在P_v与PGV相关性中，随着滤波器极点数的增加，相关系数逐渐降低，这跟速度记录的能量主要集中在0.1～10Hz的带通滤波频带内有关。随着所用滤波器阶数的增加，频带内的有效能量也被滤掉的更多，从而导致相关系数有所降低。此外，基于全P波窗获得的结果要明显优于基于P波3秒数据获得的结果，两者相差超过0.15，且随着滤波器极点数的增加，相关系数相差逐渐变大。P_vall与PGV相关性的相关系数从1阶时的0.8921，微弱降低至4阶时的0.8888，仅减小了0.004，而P_v3与PGV相关性的相关系数却减小了0.023。这说明滤波器阶数对于全P波窗数据的影响很微弱，而对P_v3的影响更大，损失的有效频段内的能量更多。又如在P_a与PGA相关性中，从表中可以看出，P_a与PGA相关性随着滤波器极点数的增加，相关系数的变化很微弱，基本可忽略。另外，基于全P波窗获得的结果要明显优于基于P波3秒数据获得的结果，相关系数接近0.85，而P波3秒结果的相关性很差，仅有0.66左右，基本不能用。

综上分析后，得出最优拟合结果如下表2所示。

表2最优相关性汇总

最优的结果为由1阶带通滤波器处理数据获得的P_vall与PGV的相关性，其相关系数最高，达到了0.8921，比次优结果的相关系数要高0.04；其次为P_aall与PGV的相关性。传统使用的根据P波3秒数据获得的相关性表现很差，远远不如经过全P波窗统计获得的结果。而在所有基于全P波窗获得的结果中，P_dall与PGV的相关性仅排到了第3名，主要是由于对加速度时程两次积分带来的较大频带外噪声(低频漂移)影响有关。

本发明还提供了现地地震动预测模型的测试。

2.1测试数据——川滇预警示范台网MEMS地震事件记录

以部署在川滇交界预警示范区域的简易烈度计GL-P2B记录到的数据作为测试用例，其动态范围达100dB，可记录到20km范围内低至M_L1.2级的地震事件。自2015年起，GL-P2B已有数百台部署在川滇交界区域，以开展区域台网预警示范应用。目前，该款设备已记录到大量的地震事件，从中选取了一部分M_L>3.5级的地震事件记录作为测试用例(24次)，共计551条GL-P2B台站记录，以对现地地震动预测模型进行验证。

2.2测试结果——川滇预警示范台网MEMS地震事件记录

对于选取的各条GL-P2B台站地震事件记录，采用与获取统计拟合关系一致的处理方式，以0.5秒/包的间隔，分别利用1～4阶的Butterworth带通滤波器进行滤波，获取P波后不同时间点的P_d、P_v和P_a值，以及台站记录的PGA和PGV值。同时以VI(4度)作为现地地震动预警阈值，并按照仪器烈度计算标准和获取到的最优拟合关系分别计算P波后不同时间点的预警烈度。此外，为了进行比较，我们还计算了实际地震动超过预警阈值的时间点。在实际使用MEMS地震事件记录计算P_d预警参数时，由于两次积分的影响，再加上本身MEMS数据质量较低，该值变化幅度太大，无法使用。因此，仅采用P_v和P_a两种预警参数进行现地地震动估算，选取的关系式为P_vall与PGV、P_aall与PGA两种类别分别在不同滤波器阶数下的拟合结果。

为了对测试结果进行评估，我们将结果按照四类进行了区分：

①成功不预警：预估现地地震动没有超过阈值，且实际观测地震动也没有超过阈值；

②成功预警：预估现地地震动超过阈值，且实际观测地震动也超过阈值；

③漏报：预估现地地震动没有超过阈值，而实际观测地震动超过阈值；

④误报：预估现地地震动超过阈值，而实际观测地震动没有超过阈值。

下表为经过1阶Butterworth带通滤波器处理后的结果。

表3测试结果汇总(1阶Butterworth带通滤波器)

类别	GL-P2B台站记录数	百分比(％)
			成功不预警	478	86.75
成功预警	40	7.26
			漏报	14	2.54
误报	19	3.45

其中，对于漏报和误报来说，分别有14条和19条台站记录。对这些记录进行分析发现，其中有9条记录计算获得的现地地震动值正好处于阈值边界，5条略大于VI(4度)，4条略小于VI(4度)。此外，有7条台站记录的数据有问题，不是正常的记录波形。对于漏报的测试用例，可以通过引入现地实际记录的地震动值作为额外的检测条件进行避免。在这种情况下，可成功处理的台站记录达到96.55％。

表4为经过1～4阶Butterworth带通滤波器处理后的测试结果汇总。从中可以看到，采用低阶滤波器处理后的测试结果表现更好，这也与前面获取P_vall与PGV、P_aall与PGA两种关系式得到的相关系数的优良情况完全吻合。

表4测试结果汇总(1-4阶Butterworth带通滤波器)

经过使用实际MEMS设备记录到的地震事件记录对获得的最优现地地震动预测模型进行测试，验证了P_vall与PGV、P_aall与PGA这两种模型的可靠性。

以上所述的实施方式，并不构成对该技术方案保护范围的限定。任何在上述实施方式的精神和原则之内所作的修改、等同替换和改进等，均应包含在该技术方案的保护范围之内。

Claims (9)

1.一种现地地震动预测模型，其特征在于，所述模型为采用Butterworth带通滤波器处理数据拟合获得的速度幅值P_v或加速度幅值P_a与峰值速度PGV或峰值加速度PGA的相关性模型。

2.如权利要求1所述的现地地震动预测模型，其特征在于，所述模型包括采用全P波窗速度幅值P_vall与峰值速度PGV的相关性模型，所述P_vall与PGV的相关性模型采用1阶时的Butterworth带通滤波器对数据进行滤波处理获得，其中所述带滤波器的频带范围为0.1～10Hz。

3.如权利要求2所述的现地地震动预测模型，其特征在于，所述模型还包括采用全P波窗加速度幅值P_aall与峰值加速度PGA的相关性模型，所述P_aall与PGA的相关性模型为采用1阶时的Butterworth带通滤波器对数据进行滤波处理获得，其中所述带滤波器的频带范围为0.1～10Hz。

4.如权利要求1-3中任一项所述的现地地震动预测模型，其特征在于，所述相关性模型的均按照式(1)进行相关性拟合而得，所述方程式为：

logPGx＝AlogPx+B (1)

式中，PGx代表PGV或PGA，Px表示P_v或P_a，而A和B为需要拟合获得的系数。

5.一种现地地震动预测模型的构建方法，至少包括下述步骤：

获取各强震动观测台站的加速度记录数据，对数据进行预处理，并筛选出标准数据；

对筛选出的各强震动观测台站事件记录数据，去掉线性趋势，手动选取P波和S波初至时间的记录数据；

对记录数据中的各强震动观测台站的垂直向加速度记录进行一次积分和两次积分以分别获取速度时程和位移时程；

对位移时程进行连续的Butterworth带通滤波器滤波处理，以获取P波幅值预警参数P_vall和P_aall；

按照仪器烈度计算标准要求计算得峰值速度PGV和峰值加速度PGA；

对于P波幅值预警参数P_vall和P_aall，按照式(1)进行相关性拟合，得到P_vall与PGV、P_aall与PGA的相关性模型。

6.如权利要求5所述的现地地震动预测模型的构建方法，其特征在于，所述对数据进行预处理，并筛选出标准数据的具体步骤包括：

对获取的加速度记录数据按照“发震时刻+震级+台站名”的方式进行重新命名，然后编写程序将所有的数据格式统一转换为文本文件格式；

对获取到的各类数据采用如下条件进行筛选：垂直向>5gal；震源距<200km；震级M≥4；同一地震事件至少有3个台站记录到；同时删除在高频能量到来时产生突跳的加速度记录。

7.如权利要求5所述的现地地震动预测模型的构建方法，其特征在于，所述对位移时程进行连续Butterworth带通滤波器滤波处理获取P波幅值预警参数P_vall和P_aall的具体方法包括：

对位移时程进行连续Butterworth带通波器滤波处理以移除由于积分操作带来的低频漂移影响，获取P波到后3s内垂直向数据对应的P_d3和全P波窗垂直向数据对应的P_dall；其中，带通滤波器的频带范围为0.075～3Hz；

按照仪器烈度计算标准针对于PGV和度PGA获取所用的Butterworth带通滤波器频带0.1～10Hz，采用1阶Butterworth滤波器对P波3秒和全P波段垂直向速度和加速度时程进行处理，以获取P_vall和P_aall。

8.如权利要求5所述的现地地震动预测模型的构建方法，其特征在于，所述按照仪器烈度计算标准要求计算得PGV和PGA中采用3分向合成值计算。

9.一种应用权利要求5-8中任一项所述的现地地震动预测模型的构建方法的构建系统，其特征在于，所述系统包括：

数据预处理模块，用于对获取的各强震动观测台站的加速度记录数据进行预处理，并筛选出标准数据；

参数计算模块，用于对筛选出的各强震动观测台站事件记录数据，去掉线性趋势，手动选取P波和S波初至时间的记录数据；对记录数据中的各强震动观测台站的垂直向加速度记录进行一次积分和两次积分以分别获取速度时程和位移时程；以及，对位移时程进行连续的Butterworth带通滤波器滤波处理，以获取P波幅值预警参数P_vall和P_aall，同时按照仪器烈度计算标准要求计算得PGV和PGA；

相关性模型建立模块，用于根据对于各类别P波幅值预警参数P_vall和P_aall，统一按照式(1)进行相关性拟合，得到P_vall与PGV、P_aall与PGA的相关性模型。

Images