基于P波位移增长特征的矩震级时域估算方法

王 江,马 强,陶冬旺,解全才,薛 韬

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

0 引言

破坏性地震发生后,地震预警(earthquake early warning, EEW)系统快速获取震源破裂信息,抢在破坏性地震波到达之前,为震区提供地震警报信息。近30 a来,EEW系统已经在多个地震频发的国家和地区运行和测试,例如墨西哥、日本、美国、中国等[1-6]。震时的震级快速估计是EEW系统中关键且困难的一环。震级表征震源能量的大小,决定了地面强震动强弱程度,是EEW系统发布预警信息的关键内容,震级估计结果的可靠性直接影响地震动预测精度和震中附近预警区范围[7]。实时震级估算基本原理是利用早期破裂辐射的有限波形估算最终震级,其理论依据也是震源理论的热点问题:地震破裂是否存在确定性,即地震的初始破裂是否控制最终震源规模。然而,相关研究并未形成统一认识[8-14],不过近年地震学界逐渐接受破裂弱确定性的观点[15-16]。在EEW系统实际应用中,采用早期破裂有限波形也难以获取理想震级估算结果,一方面离线数据估算结果离散性大;另一方面常出现大震震级严重低估及小震震级高估问题。目前,实时震级估算主流方法是通过分析早期破裂辐射的P波和S波携带震源信息,提取与震级相关的特征参数,采用单一或组合特征参数与震级经验关系实现震级估算。NAKAMURA[17]最早提出利用初始地震动卓越周期参数估算震级的方法。ALLEN等[18],OLSON等[8]和KANAMORI[9]先后提出了利用P波频率成分预测震级的方法,认为地震破裂普遍存在确定性。另一个估算震级的特征参数是P波位移幅值,WU等[19]使用3 s P波的位移幅值(Pd)预测峰值地震动(PGA和PGV),实现地震破坏潜力的快速评估。在此研究基础上,WU等[10]和ZOLLO等[20]通过建立Pd的衰减关系实现震级快速估计。随后考虑地震波频率和能量的特征参数不断涌现,试图得到与震级更准确可靠的经验关系[21-24]。近年来,人工智能方法融合初始破裂辐射地震动的多种特征参数和波形数据,一定程度上提升了实时估算震级的准确度[25-27]。然而,估算或预测震级统计方法采用的特征参数并不表征地震的绝对能量,存在对震源破裂传播过程考虑不足,不涉及破裂尺度,不能直接量化震源破裂能量,且面向地震场景应用仍存在震级估算结果离散,大震震级低估不足等问题,难以保证地震动预测精度,这在一定程度上制约着EEW系统震时防震减灾能效。

本文从远场P波位移参数与地震矩的理论关系出发[28-29],尝试从时域Pd曲线提取特征参数来估计地震破裂能量。在理论上,远场位移波形记录与震源时间函数(source time function, STF)基本一致[30]。中、小地震的震源破裂过程相对简单,可将其STF形状简化为等腰三角形,此时STF三角形面积就是地震的地震矩。由此,本文推导了小震和中震的地震矩与P波位移波形参数(位移峰值Pd和峰值时刻tpd)的理论关系。在观测分析中,一方面统计了日本地震的STF的矩率峰值等参数与地震矩关系;另一方面分析强震动观测记录中不同震级地震P波位移增长曲线(简称Pd曲线)的特征,从观测的Pd曲线中提取表征地震矩特征的构造参数Pdt,并验证其估算矩震级的可行性。本文从STF函数和强震动观测记录两方面分别验证地震矩和观测P波记录中相关特征参数的理论关系,为EEW系统实时矩震级估算提供了一种可靠的新方法。

1 方法

1.1 远场P波位移参数与地震矩的理论关系

远场位移波形正比于震源地震矩释放过程,即与震源时间函数STF的形状一致。因此,P波位移增长曲线(Pd曲线)能够反映出地震震级和破裂尺度等震源参数信息[28]。在均匀弹性半空间介质中,位错点源模型的远场P波位移波形可描述为[31]:

(1)

(2)

中震和小震震源规模小,破裂过程相对简单,其STF一般可简化近似为等腰三角形,因此,远场P波波形也近似为等腰三角形,其积分表达为:

(3)

式中:Pd为P波位移峰值;Td为P波脉冲持时。

(4)

由此可知,当从P波位移波形中测定参数Pd和tpd时,即可确定该次地震的地震矩大小。

此外,地震矩与矩震级的理论关系式为:

logM0=1.5Mw+9.1

(5)

式(4)代入到式(5)中得P波位移波形特征参数与矩震级的表达式:

log(Pd×tpd)=logM0+log (A/R)=1.5Mw+log (A/R)+9.1

(6)

1.2 STF特征参数与矩震级的统计关系

图1 地震STF的特征参数示意图和特征参数与矩震级的统计关系

2 数据筛选和处理

本文采用日本K-net强震动数据库的竖向分量记录。考虑浅源地震破坏严重,设置以下筛选准则:震源深度H≤60 km;震源距小于100 km;P波起始3 s加速度峰值大于3 Gal,记录3 s P波与事前记录的信噪比大于8。选用发生于1997年1月至2011年6月间的地震共计467次,地震震级介于4.0≤Mj≤7.3,数据集记录共计2745条。图2为研究区地震动的震中位置分布,图3为采用强震动数据统计直方图。数据库中地震事件的矩震级目录来自F-net数据库。

图2 研究区地震震中分布图

图3 强震动数据统计直方图

为了观测P波位移增长特征及其传播衰减效应,数据预处理对强震动记录减前5 s均值完成基线归零,P波初至采用STA/LTA粗捡拾和AIC精捡拾组合方法拾取[34],人工检查拾取结果,修正拾取有误的记录。然后通过P波和S波理论到时差明确S波初至时刻(采用Vp=5.333 km/s,Vs=3.2 km/s)。对基线校正后加速度数据进行二次积分获得位移时程,用低频截止频率为0.075 Hz的二阶巴特沃斯高通滤波器去除积分引起的长周期漂移。在时间窗从0.01 s逐渐扩大到10 s(间隔0.01 s)的位移时间序列上,持续测量Pd参数得到Pd曲线。由于地震动传播的几何扩散和衰减效应,不同台站的Pd曲线需校正到统一参考震源距,在此,采用相似数据集统计的Pd衰减关系进行校正[35],式中系数见表1。

表1 Pd衰减关系系数表

logPd=A+BM+ClogR

(7)

由于台站空间分布的几何原因,大震源距记录的数量要远多于近源记录,如图3所示,以往地震初始破裂特征观测中常直接使用一定范围记录的均值以消除震源辐射模式效应,这会导致大震源距记录成为均值结果的主要贡献者,引入地震动长距离传播高频衰减的干扰。为分析地震动传播衰减效应,尤其是对参数Pd和tPd的影响,本节将强震动记录按震源距分到5个区间内([0,20),[20,40),[40,60),[60,80),[80,100] km),每个震源距段内记录校正的参考震源距是其对应震源距段的中间值(如震源距段0～20 km的校正参考距离Dr为10 km),再计算各震源距段内Pd曲线的平均值,计算过程如图4所示,图4中彩点指示采用记录的S波理论初至时刻。图5按震源距范围汇总了所有震级记录的位移增长均值曲线。

图4 Pd曲线计算过程图和典型地震的5个震源距段内Pd均值曲线

图5 5个震源距段内Pd均值曲线和Pd增长模型

3 结果

3.1 Pd曲线特征

图4(b)中6个不同震级地震的Pd均值曲线表明近源(震源距20 km以内)地震动幅值远大于其他区域的地震动,随传播距离增加,Pd曲线幅值幅度快速衰减,引起Pd曲线初始上升斜率的下降。排除S波干扰,P波位移增长达到峰值水平的时间基本不受震源距影响。按震源距范围对比不同震级的P波位移增长过程如图5所示。图5(a)～(e)中不同震级地震的P波的Pd均值曲线分为2个主要阶段:上升段和稳定段。其中上升段Pd均值曲线先陡增后缓赠,陡增段遵循指数上升趋势,随后放缓直至饱和,陡增和缓增阶段间无明显界限;在稳定段,大震P波位移峰值幅值大,位移上升时间长。根据P波位移增长特征,建立了Pd增长模型,包括陡增与缓增的上升阶段和稳定水平段,并定义了曲线的拐角特征参数Pd和tpd,如图5(f)所示。

由图5可知,Pd增长模型中拐角参数(tpd,Pd)显然与震级有关,震级越大,拐角参数越大。而其中参数Pd作为EEW地震实时震级估算的特征参数,Pd曲线表明小震记录中采用很短时间窗便可测量稳定Pd,而大震记录的测量稳定Pd需要时间要远远高于小震。因此,采用统一时间窗记录估算震级的策略不适用所有震级范围的地震事件。

3.2 构造参数Pdt

根据地震动传播衰减特征,将各台站记录的位移时程按式(7)校正至参考震源距10 km。考虑单一记录的Pd曲线对高频噪声信号十分敏感,每次地震Pd均值曲线按下式优化计算:

(8)

图6 构造参数Pdt计算过程

根据P波位移参数与地震矩的理论表达(式(6)),构造参数Pdt为参数Pd和tpd的乘积:

(9)

由此,便可计算随P波时间窗逐渐增大的构造参数Pdt参数时间序列(图6(c))。

按式(10)形式统计全P波窗内Pd均值曲线中特征参数(Pd,tpd)和构造参数Pdt与地震震级的关系,包括日本震级Mj和矩震级Mw,各参数线性回归结果如图7所示。

图7 P波特征参数与地震震级的统计结果

log(Para)=aM+b

(10)

式中,Para为特征参数,分别为Pd,tpd和Pdt;M为震级,统计公式系数a,b,如表2和表3所示。

表2 P波特征参数与Mj统计系数表

3.3 矩震级估算

P波的Pd参数是地震预警实时震级预测的主要参数。本节采用表2和表3拟合系数,对每次地震Pd均值曲线测量不同时窗内的构造参数Pdt与参数Pd预测矩震级,预测结果分别如图8和图9所示。图8中6个时间窗下构造参数Pdt预测结果表明,在本文所涉及震级范围内地震,随时窗增大Pdt参数预测结果逐渐趋向于目录矩震级;而在相同时间窗下,构造参数Pdt的预测震级的误差标准差σ均小于Pd参数预测结果。对5.5级以下的地震,构造参数Pdt在2 s时间窗就可以可靠估算震级(图8(b))。构造参数Pdt在5 s以上时窗的预测结果未出现小震的高估和大震的低估。注意根据式(10),图7中σ是对数下特征参数,即log10(Para)的拟合误差标准差,而图8和图9中σ是对预测矩震级误差标准差,即Mwpre=[log(Para)-b]/a。因此图8和图9的预测矩震级的误差标准差不同于图7中的拟合标准差。

图8 P波构造参数Pdt预测矩震级结果

图9 P波特征参数Pd预测矩震级结果

图10绘制了采用3 s和5 s时窗P波参数Pdt和Pd估算的矩震级误差随目录矩震级的分布图。由图10可知,采用3 s P波的参数Pdt预测结果要优于Pd参数,但是2个参数对6级以上地震均有不同程度矩震级低估,这表明3 s破裂能量还不能表征最终震级大小。相比来看,采用5 s P波,2个特征参数对6级以上地震的震级低估得到有效改善。对比2个时窗的矩震级预测结果可知,参数Pdt表现均优于参数Pd,证实了构造参数Pdt估算矩震级的可靠性。

图10 P波构造参数Pdt和特征参数Pd预测矩震级误差对比

为验证构造参数Pdt在实际地震震例中的性能,本节计算了图6中震例的矩震级增长过程。该次地震矩震级6.6级,最终预测结果为6.3级,稳定震级所需时间约4.3 s,预测误差为0.3级(图6(c))。除此以外,在该节还利用发生于日本宫城县北部的Mj6.2地震开展测试工作,本次地震发震时间是2003年7月26日,矩震级为6.1级,震源深度12 km。经筛选,使用该次地震震源距100 km以内29个台站的高质量数据。由参数Pdt预测矩震级的结果和震例的震中和台网位置如图11所示。由图可知,构造参数Pdt随时间窗增加,预测矩震级逐渐增大,在6.6 s达到最终的6.4级,预测误差也为0.3级。应注意的是,不同震源距记录的全P波波形长度差异,设定小时窗下可用的记录多,而设定大时窗可用的记录较少,因此随P波时窗增加,所用记录个数逐渐减少。考虑单一台站计算误差较大,建议设置P波时窗在少于3条记录时终止计算,以此时的预测矩震级结果作为最终矩震级。这2次地震均考虑使用地震记录较丰富的震例,预警矩震级结果可靠,未出现较大误差。

图11 Mw6.1地震随P波时窗增加构造参数Pdt预测矩震级结果

4 讨论与结论

P波位移增长特征一方面可用于识别不同震级地震初始破裂的差异特征[13,16,35],另一方面可快速计算震源参数,包括估计地震矩、破裂尺度和应力降[29,33]。本文按震源距范围分类,弱化传播衰减对Pd曲线的影响,结果表明不同震级地震强震动记录Pd曲线初始陡增趋势并无明显差异,但随震源距增大,大震记录的陡增段斜率放缓,揭示了地震的初始破裂并无明显差别,相关破裂初始差异的报道可能忽略了传播衰减的影响。

利用Pd曲线参数估计矩震级的方法仅适用于可近似为等腰三角形的小震和中震,因此,对破裂过程复杂的大震的通用性需进一步研究。本文未考虑场地效应对Pd曲线的影响,实际地震动还应考虑场地效应校正,可能会对本方法预测震级结果有一定改善。

在理论上,地震矩可由P波位移波形参数直接计算。本文采用日本地震的震源时间函数进行验证。通过分析近场强震动记录P波位移增长特征,本文提取了Pd均值曲线中位移峰值Pd和峰值时刻tpd,并构建预测矩震级的特征参数Pdt。通过离线数据的进行地震矩震级估计测试,得到如下结论:

1)P波波形参数位移峰值Pd和峰值时刻tpd均由地震震级控制,震级越大,特征参数越大。

2)中、小地震的矩震级可由地震P波位移波形参数Pd和tpd的乘积计算,理论推导和观测记录(包括STF和近场强震动记录)中均得到证实。

3)采用固定时窗(如3 s)的记录估算震级的策略不适用所有震级范围的地震,尤其是大震。

4)完全由经验关系估算地震震级在实时地震预警系统是可行的,但偏离了理论地震矩与矩震级的关系,这可能是导致震级估算误差较大和大震低估的一个原因。

致谢：感谢日本Kik-net数据中心(https://www.kyoshin.bosai.go.jp/)为本文研究提供的强震动数据,感谢SCARDEC数据库(http://scardec.projects.sismo.ipgp.fr/)为本文提供的震源时间函数。