一种基于地震P波特征参数的有限破裂模板快速匹配方法与流程

一种基于地震p波特征参数的有限破裂模板快速匹配方法
技术领域
1.本发明涉及地震预警技术领域，具体地说是一种基于地震p波特征参数的有限破裂模板快速匹配方法。


背景技术：

2.地震预警是指在地震发生时技术系统利用通信网络以及速度快但强度弱的纵波和速度慢却强度大的横波之间的时间差，在破坏性地震波到达前发出预警警报，减少人员伤亡和财产损失。
3.现有地震烈度速报与预警中，包含测震、强震、烈度仪等三网融合的密集台网，其中核心关键就是如何实时、快速估算断层破裂特征参数，但目前尚未实现，现有技术中仍旧以传统的点源算法为主，虽然按计划拟引入finder系统，但是从前面的分析结果来看，finder算法在时效性上仍然无法完全满足地震预警的高时效性要求。
4.基于以上原因，本发明设计了一种基于地震p波特征参数的有限破裂模板快速匹配方法，选用整个p波段数据拟合结合有限破裂模板匹配，可以进一步提高获取断层破裂特征参数的时效性，在助于盲区范围的减小和有效预警时间的增加。


技术实现要素：

5.本发明的目的是克服现有技术的不足，提供一种基于地震p波特征参数的有限破裂模板快速匹配方法，选用整个p波段数据拟合结合有限破裂模板匹配，可以进一步提高获取断层破裂特征参数的时效性，在助于盲区范围的减小和有效预警时间的增加。
6.为了达到上述目的，本发明提供一种基于地震p波特征参数的有限破裂模板快速匹配方法，包括数据资料收集与预处理、最优地震动预测模型构建、有限破裂模板库生成、有限破裂模板匹配方式和软件系统检验；
7.数据资料收集与预处理，具体步骤为：
8.s1，对数据按照“发震时刻+震级+台站名”的方式进行重新命名；
9.s2，编写程序将所有的数据格式统一转换为文本文件格式，格式只包括实际记录数据，不包括头信息，并将单位统一转换为cm/s2即gal；
10.s3，对获取到的各类数据进行筛选，筛选的条件为：垂直向》5gal；震源距《200km；m≥4，m为震级；同一地震事件至少有3个台站记录到；
11.s4，排除由于传感器的弹簧系统或者电子器件引起的个别加速度突跳的数据；
12.最优地震动预测模型构建包括初期p波预警参数选取与计算，统计关系式拟合和最优相关性汇总；
13.初期p波预警参数选取与计算为：
14.s10，对筛选出的各台站事件记录数据，去掉线性趋势，手动选取p波和s波到时；
15.s20，对各台站的垂直向加速度记录，即ud，进行一次积分和两次积分以分别获取速度时程和位移时程，并对位移时程进行连续巴特沃斯带滤波器滤波处理，滤波器处理频
带范围为0.075～3hz，以移除由于积分操作带来的低频漂移影响，获取p波到后3s内垂直向数据对应的pd值即p
d3
，和全p波窗垂直向数据对应的pd值，即p
dall
；
16.s30，同时采用1～4阶四种巴特沃斯带滤波器进行处理，即每类预警参数计算获得4个值；
17.s40，采用1～4阶四种巴特沃斯带滤波器对p波3秒和全p波段垂直向速度和加速度时程进行处理，以获取p
v3
、p
vall
、p
a3
和p
aall
；
18.统计关系式拟合为对于各类别p波幅值预警参数按照下式进行相关性拟合：
19.式一：logpgx＝alogpx+b
20.pgx代表pgv或pga，px表示pd、pv或pa，而a和b为需要拟合获得的系数；
21.有限破裂模板库生成为：采用图像识别技术并通过估计地震当前的矩心位置、破裂长度l和走向θ来实时地自动检测断层破裂的地表投影，具体公式为：
22.式二：lg(pga)＝2.206+0.532m-1.954lg(r+2.018e
0.406m
)
23.m为震级；r为震中距，单位km，在实际使用该衰减关系式时，进行了一定的变化：即当m《5级时，r为震中距；当m≥5时，r为至有限破裂模板中假定有限断裂对应的断层距，即r
jb
；
24.断层破裂的长度按照以下公式计算：
25.式三：lg(l)＝(m-4.33)/1.49
26.有限破裂模板匹配方式具体步骤为：
27.s100，对各台站利用p波特征参数估测的峰值加速度pga或近台s波到达后的峰值加速度pga进行联立插值，并依据设定的地震动pga阈值，得到一个二值图像；
28.s200，通过相关性匹配，从生成的破裂模板库中检索最大相关位置，将其作为矩心，相关性匹配公式为：
29.式四：r(x,y|l,θ)＝i(x,y)
★
t(x,y|l,θ)
30.t(x,y|l,θ)为生成的破裂模板集，i(x,y)为当前根据插值数据生成的二值图像；
31.将上述式四转换至傅里叶域计算：
32.式五：
33.式五表示将空间相关性转换成傅里叶变换和在小波域(k
x
,ky)的乘积；
34.s300，得到与二进制图像最佳匹配的模板，将其线源参数作为当前检索得到的震源模型，公式为：
[0035][0036]
x'＝0至w-1，y'＝0至h-1，w
×
h为生成的模板库数量；
[0037]
s400，根据更多的台站被触发和每个台站记录到更多的波形数据，持续迭代地进行模板匹配，以对破裂长度和方向进行修正。
[0038]
同现有技术相比，本发明在震后同一时刻得到的结果相对于finder算法结果在速度上要快3s左右，个别震例结果要快5s，从而更快地趋近于标准参考值；可以采用多报预警
信息和由内向外的发布策略，随着破裂长度的增长而逐渐增加报次，从而尽可能的缩小盲区范围，增加有效预警时间。
附图说明
[0039]
图1为本发明所用震例分布情况的示意图。
[0040]
图2为本发明有限破裂模板匹配流程示意图。
具体实施方式
[0041]
现结合附图对本发明做进一步描述。
[0042]
参见图1～2，本发明提供了一种基于地震p波特征参数的有限破裂模板快速匹配方法，包括数据资料收集与预处理、最优地震动预测模型构建、有限破裂模板库生成、有限破裂模板匹配方式和软件系统检验；
[0043]
数据资料收集与预处理，具体步骤为：
[0044]
s1，对数据按照“发震时刻+震级+台站名”的方式进行重新命名；
[0045]
s2，编写程序将所有的数据格式统一转换为文本文件格式，格式只包括实际记录数据，不包括头信息，并将单位统一转换为cm/s2即gal；
[0046]
s3，对获取到的各类数据进行筛选，筛选的条件为：垂直向》5gal；震源距《200km；m≥4，m为震级；同一地震事件至少有3个台站记录到；
[0047]
s4，排除由于传感器的弹簧系统或者电子器件引起的个别加速度突跳的数据；
[0048]
最优地震动预测模型构建包括初期p波预警参数选取与计算，统计关系式拟合和最优相关性汇总；
[0049]
初期p波预警参数选取与计算为：
[0050]
s10，对筛选出的各台站事件记录数据，去掉线性趋势，手动选取p波和s波到时；
[0051]
s20，对各台站的垂直向加速度记录，即ud，进行一次积分和两次积分以分别获取速度时程和位移时程，并对位移时程进行连续巴特沃斯带滤波器滤波处理，滤波器处理频带范围为0.075～3hz，以移除由于积分操作带来的低频漂移影响，获取p波到后3s内垂直向数据对应的pd值即pd3，和全p波窗垂直向数据对应的pd值，即p
dall
；
[0052]
s30，同时采用1～4阶四种巴特沃斯带滤波器进行处理，即每类预警参数计算获得4个值；
[0053]
s40，采用1～4阶四种巴特沃斯带滤波器对p波3秒和全p波段垂直向速度和加速度时程进行处理，以获取p
v3
、p
vall
、p
a3
和p
aall
；
[0054]
统计关系式拟合为对于各类别p波幅值预警参数按照下式进行相关性拟合：
[0055]
式一：logpgx＝alogpx+b
[0056]
pgx代表pgv或pga，px表示pd、pv或pa，而a和b为需要拟合获得的系数；
[0057]
将各种拟合获得的统计关系式中每项最优结果进行汇总，按照相关系数倒序排列，其结果见表1所示：
[0058]
相关性滤波器极点数系数a系数b标准差std相关系数rp
vall vs.pgv10.94770.88560.27790.8921p
aall vs.pga20.84270.94060.26270.8512
p
aall vs.pga10.84860.89600.26340.8503p
dall vs.pgv40.60381.23550.32590.8481p
aall vs.pgv41.0074-0.50460.33570.8379p
vall vs.pga10.71492.09980.28140.8270p
d3 vs.pga30.50342.55020.34000.7338
[0059]
表1
[0060]
由于p
aall
与pga的相关性在1阶和2阶滤波器的处理结果获得的相关系数差别不太，所以这里将1阶和2阶的结果都列举出来。
[0061]
从表中可以看出，最优的结果为由1阶带通滤波器处理数据获得的p
vall
与pgv的相关性，其相关系数最高，达到了0.8921，比次优结果的相关系数要高0.04。传统使用的根据p波3秒数据获得的相关性表现很差，远远不如经过全p波窗统计获得的结果。而在所有基于全p波窗获得的结果中，p
dall
与pgv的相关性仅排到了第3名，这跟前面分析得到的结果一致，主要还是由于对加速度时程两次积分带来的较大频带外噪声(低频漂移)影响有关。
[0062]
有限破裂模板库生成为：采用图像识别技术并通过估计地震当前的矩心位置、破裂长度l和走向θ来实时地自动检测断层破裂的地表投影，具体公式为：
[0063]
式二：lg(pga)＝2.206+0.532m-1.954lg(r+2.018e
0.406m
)
[0064]
m为震级；r为震中距，单位km，在实际使用该衰减关系式时，进行了一定的变化：即当m《5级时，r为震中距；当m≥5时，r为至有限破裂模板中假定有限断裂对应的断层距，即r
jb
；
[0065]
断层破裂的长度按照以下公式计算：
[0066]
式三：lg(l)＝(m-4.33)/1.49
[0067]
将震级范围扩展至m2.5～m8.0，对应的断层破裂长度为l＝0.06～300km，按照震级0.1级为间隔，对应于走向θ＝0
°
，共生成56个有限破裂模板。在进行模板匹配时，会对模板按照1
°
为间隔旋转至180
°
，以获取震源破裂的方向。
[0068]
有限破裂模板匹配方式具体步骤为：
[0069]
s100，对各台站利用p波特征参数估测的峰值加速度pga或近台s波到达后的峰值加速度pga进行联立插值，并依据设定的地震动pga阈值，得到一个二值图像；
[0070]
s200，通过相关性匹配，从生成的破裂模板库中检索最大相关位置，将其作为矩心，相关性匹配公式为：
[0071]
式四：r(x,y|l,θ)＝i(x,y)
★
t(x,y|l,θ)
[0072]
t(x,y|l,θ)为生成的破裂模板集，i(x,y)为当前根据插值数据生成的二值图像；
[0073]
将上述式四转换至傅里叶域计算：
[0074]
式五：
[0075]
式五表示将空间相关性转换成傅里叶变换和在小波域(k
x
,ky)的乘积；
[0076]
s300，得到与二进制图像最佳匹配的模板，将其线源参数作为当前检索得到的震源模型，公式为：
[0077][0078]
x'＝0至w-1，y'＝0至h-1，w
×
h为生成的模板库数量；
[0079]
s400，根据更多的台站被触发和每个台站记录到更多的波形数据，持续迭代地进行模板匹配，以对破裂长度和方向进行修正。
[0080]
我们利用7次历史地震事件(见表2)记录作为软件系统输入，模拟和分析研发的震源破裂特征快速估测软件系统的性能，以便为后期系统的完善与优化提供基础。这里忽略了由数据传输和发布警报所导致的时间延时。在验证过程中，只要一旦两个及以上台站的估测或实测的pga超过设定阈值时，系统就开始进行处理。
[0081][0082][0083]
表2
[0084]
系统测试结果表明：
[0085]
(1)利用本方法在震后同一时刻得到的结果相对于finder算法结果在速度上要快3s左右，个别震例结果要快5s，从而更快地趋近于标准参考值；
[0086]
(2)在破裂的初期，由于受到地震辐射多样性、场地、传播路径等因素的影响，获得的走向θ会存在较大的变化。随着时间的持续，θ才会逐渐收敛至标准参考值；
[0087]
(3)对于m7.0级以下地震，震后6～10s即可获得较稳定的破裂特征参数结果，而对于m7.0+地震，则需要更长的时间才能获得较稳定的破裂特征参数，尤其是对于像汶川这种特大地震，其较稳定结果一般需要到震后40s以后。当然，这与台网密度密切相关；
[0088]
(4)在具体使用时，可以采用多报预警信息和由内向外的发布策略，第一报可以在
估算的破裂长度l达到10km时就发布，并随着破裂长度的增长而逐渐增加报次，从而尽可能的缩小盲区范围，增加有效预警时间；
[0089]
(5)对于所用的震级与破裂长度的经验关系(wells&coppersmith,1994)，从其结果来看，对于m7.0级左右地震事件，即使得到的破裂长度明显超过标准参考值，获得的震级m也明显比标准参考值小。因此，未来有必要对该关系式进行本地化处理，以获取更加适用于国内的震级与破裂长度的经验关系。
[0090]
以上仅是本发明的优选实施方式，只是用于帮助理解本技术的方法及其核心思想，本发明的保护范围并不仅局限于上述实施例，凡属于本发明思路下的技术方案均属于本发明的保护范围。应当指出，对于本技术领域的普通技术人员来说，在不脱离本发明原理前提下的若干改进和润饰，这些改进和润饰也应视为本发明的保护范围。
[0091]
本发明从整体上解决了现有技术中仍旧以传统的点源算法为主，虽然按计划拟引入finder系统，但是从前面的分析结果来看，finder算法在时效性上仍然无法完全满足地震预警的高时效性要求的问题，通过选用整个p波段数据拟合结合有限破裂模板匹配，可以进一步提高获取断层破裂特征参数的时效性，在助于盲区范围的减小和有效预警时间的增加。

技术特征：
1.一种基于地震p波特征参数的有限破裂模板快速匹配方法，其特征在于，包括数据资料收集与预处理、最优地震动预测模型构建、有限破裂模板库生成、有限破裂模板匹配方式和软件系统检验；所述数据资料收集与预处理，具体步骤为：s1，对数据按照“发震时刻+震级+台站名”的方式进行重新命名；s2，编写程序将所有的数据格式统一转换为文本文件格式，所述格式只包括实际记录数据，不包括头信息，并将单位统一转换为cm/s2即gal；s3，对获取到的各类数据进行筛选，所述筛选的条件为：垂直向>5gal；震源距<200km；m≥4，所述m为震级；同一地震事件至少有3个台站记录到；s4，排除由于传感器的弹簧系统或者电子器件引起的个别加速度突跳的数据；所述最优地震动预测模型构建包括初期p波预警参数选取与计算，统计关系式拟合和最优相关性汇总；所述初期p波预警参数选取与计算为：s10，对筛选出的各台站事件记录数据，去掉线性趋势，手动选取p波和s波到时；s20，对各台站的垂直向加速度记录，即ud，进行一次积分和两次积分以分别获取速度时程和位移时程，并对位移时程进行连续巴特沃斯带滤波器滤波处理，所述滤波器处理频带范围为0.075～3hz，以移除由于积分操作带来的低频漂移影响，获取p波到后3s内垂直向数据对应的p
d
值即pd3，和全p波窗垂直向数据对应的pd值，即pd
all
；s30，同时采用1～4阶四种巴特沃斯带滤波器进行处理，即每类预警参数计算获得4个值；s40，采用1～4阶四种巴特沃斯带滤波器对p波3秒和全p波段垂直向速度和加速度时程进行处理，以获取p
v3
、p
vall
、p
a3
和p
aall
；所述统计关系式拟合为对于各类别p波幅值预警参数按照下式进行相关性拟合：式一：logpgx＝alogpx+b所述pgx代表pgv或pga，px表示pd、p
v
或p
a
，而a和b为需要拟合获得的系数；所述有限破裂模板库生成为：采用图像识别技术并通过估计地震当前的矩心位置、破裂长度l和走向θ来实时地自动检测断层破裂的地表投影，具体公式为：式二：lg(pga)＝2.206+0.532m-1.954lg(r+2.018e
0.406m
)所述m为震级；r为震中距，单位km，在实际使用该衰减关系式时，进行了一定的变化：即当m<5级时，r为震中距；当m≥5时，r为至有限破裂模板中假定有限断裂对应的断层距，即r
jb
；断层破裂的长度按照以下公式计算：式三：lg(l)＝(m-4.33)/1.49所述有限破裂模板匹配方式具体步骤为：s100，对各台站利用p波特征参数估测的峰值加速度pga或近台s波到达后的峰值加速度pga进行联立插值，并依据设定的地震动pga阈值，得到一个二值图像；s200，通过相关性匹配，从生成的破裂模板库中检索最大相关位置，将其作为矩心，相关性匹配公式为：式四：r(x,y|l,θ)＝i(x,y)＊t(x,y|l,θ)
所述t(x,y|l,θ)为生成的破裂模板集，i(x,y)为当前根据插值数据生成的二值图像；将上述式四转换至傅里叶域计算：式五：式五表示将空间相关性转换成傅里叶变换和在小波域(k
x
,k
y
)的乘积；s300，得到与二进制图像最佳匹配的模板，将其线源参数作为当前检索得到的震源模型，公式为：所述x'＝0至w-1，y'＝0至h-1，w
×
h为生成的模板库数量；s400，根据更多的台站被触发和每个台站记录到更多的波形数据，持续迭代地进行模板匹配，以对破裂长度和方向进行修正。

技术总结
本发明涉及地震预警技术领域，具体地说是一种基于地震P波特征参数的有限破裂模板快速匹配方法，包括数据资料收集与预处理、最优地震动预测模型构建、有限破裂模板库生成、有限破裂模板匹配方式和软件系统检验，同现有技术相比，本发明在震后同一时刻得到的结果相对于FinDer算法结果在速度上要快3s左右，个别震例结果要快5s，从而更快地趋近于标准参考值；可以采用多报预警信息和由内向外的发布策略，随着破裂长度的增长而逐渐增加报次，从而尽可能的缩小盲区范围，增加有效预警时间。增加有效预警时间。增加有效预警时间。


技术研发人员：彭朝勇 程振鹏 郑钰 徐志强
受保护的技术使用者：中国地震局地球物理研究所
技术研发日：2023.06.09
技术公布日：2023/9/12