GP断层破裂模型在青海门源MS6.9地震强地面运动模拟的适用性研究

罗超,曹晓雨,高阳,徐飞,冯怀平,王昊*

1 石家庄铁道大学 省部共建交通工程结构力学行为与系统安全国家重点实验室,石家庄 050043 2 石家庄铁道大学 土木工程学院,石家庄 050043 3 石家庄铁道大学 安全工程与应急管理学院,石家庄 050043

0 引言

北京时间2022年1月8日凌晨1时45分青海省海北州门源县发生面波震级为6.9级的地震,震源深度约为10 km(尹晓菲等,2022).本次地震导致的人员伤亡较少,造成海北州门源等三个县千余户受灾,9人受伤,超过四千间房屋损坏(王祖东等,2022).本次地震引发了工程结构的显著破坏,靠近发震断层的兰新铁路硫磺沟大桥和祁连山隧道破坏严重(韩帅等,2022).震害表明,近场地震动对工程结构的破坏力巨大,当类似的地震发生在人口密集和重大工程附近时后果不堪设想.

地震发生后,许多学者和研究机构都对门源地震区域断层和震源参数开展了研究工作.李振洪等(2022)反演了本次地震的震源破裂过程,结果表明本次地震的类型为走滑型地震,子断层最大滑动量为3.5 m,大约在地下4 km处.潘家伟等(2022)通过详细野外考查分析得出,本次地震以左旋走滑运动为主,最大走滑位移约为3.7 m,形成了约27 km长的破裂带,震源深度较浅.尹欣欣等(2022)通过双差层析成像方法对2009—2022年之间的门源及周边地区地震事件进行了精确定位处理,定位本次门源地震的震源深度为13.3 km,震中位置为(37.76°N,101.25°E).

在震后重建强地面运动场,是确定震后受灾区域和估计未来强震影响场的重要手段.同时随着一系列地震动输入方法的提出和广泛应用(Bielak et al.,2003; Luo et al.,2019,Wang et al.,2021a),对局部场地中密集地震动的获取提出了要求.已有许多研究人员开展了强地面运动模拟方面的研究,其中随机有限断层法被广泛使用.丁阳等(2018)使用改进的有限断层法,对九寨沟MS7.0地震进行了模拟,结果显示,部分模拟结果与强震观测记录存在一定差异.周红等(2021)通过在随机有限断层法中考虑破裂过程,模拟了漾濞6.4级地震强地震动场,对地震动场的时空分布规律进行了分析.李春果等(2021)等考虑了玛多MS7.4地震的多个震源破裂模型,基于随机有限断层法开展了强地震动模拟.通过与预测中位值对比,模拟记录的PGA、PGV结果较为理想,但当与远场观测记录相比时,基于模拟记录的地震烈度值整体上偏低1～2度.强生银等(2021)基于随机有限断层三维地震动模拟方法,给出了漾濞地震中28个强震动台站的三分量加速度时程模拟结果.傅磊等(2023)通过随机有限断层法模拟了2022年芦山MS6.1地震,分析了地震中9个距离断层小于100 km的台站的模拟记录.结果显示,模拟记录与观测记录在加速度时程的S波部分符合良好,在8 s以下周期范围内傅里叶振幅谱和反应谱形状和幅值基本一致.综上所述,随机有限断层法虽然计算效率较高,但是只能从宏观上考虑地形和地质条件的影响,不能完全考虑速度结构、局部场地效应对地震动的影响.

近年来,随着计算能力的提高,研究人员也开展了有关确定性模拟方法的研究.常莹等(2012)在模拟汶川MS8.0地震的强地震动场时,通过对离散断层面分级和对各级断层面破裂上升时间分解的方法,增加了震源的高频成分.结果表明,PGA模拟结果的误差不超过50%,傅里叶谱在4 Hz以上频段比观测值小.赵宏阳和陈晓非(2017)基于辽宁省海城地区地质资料和已发表的海城MS7.3地震震源机制反演结果,模拟了1975年海城地震的强地震动场.计算结果表明,地震动场的分布特征受到包括震源、地形和地质条件等因素的影响.但是受到模型精度的影响,理论和实际等震线图存在一定差异.贺春晖等(2017)提出了震源-介质-坝址峡谷场地的多尺度模拟方法,基于谱元法模拟了坝址区域的空间三维地震动场.Wang 等(2021b)在宽频带地震动模拟中,通过将震源分成多个层叠加模拟断层破裂过程实现高频信号的模拟,并以1994年北岭地震为例进行了验证.尹晓菲等(2022)使用曲线网格有限差分法,基于门源区域的地形、速度结构模型以及反演震源破裂过程的结果,模拟了门源地震的强地震动场.受模型精度和计算资源的制约,模拟的地震动场有效频率不超过0.4 Hz,模拟结果的烈度较强震观测记录偏低.上述文献中,均使用了震后反演得到的断层破裂模型作为输入,来模拟设定地震的强地面运动.但是不同地震的断层破裂模型迥异,能否直接应用在其他有潜在发生地震区域的强地震动模拟中尚存疑问.

Graves和Pitarka提供了另一种思路,他们提出了预测可能发生的断层破裂过程的GP法(Graves and Pitarka,2004,2010,2015,2016),可用于开展宽频带地面运动模拟.Pitarka和Graves(Pitarka et al.,2020)使用GP法和GP&IM混合法建立断层破裂模型,对日本熊本MW7.0地震进行了宽频带地震动模拟;Pitarka等(2022)使用GP法对美国加利福尼亚州RidgecrestMW7.1地震进行了宽频带地震动模拟,两篇文章通过将计算结果与强震观测记录对比,展现出了GP法在美国和日本地区优异的适用性.但是GP法是否适用于中国地区的强地震动模拟,目前仅巴振宁等(2023)使用GP模型对漾濞地震进行了模拟,研究工作还有待进一步开展.

研究表明,中国的城市50年内可能发生6级左右直下型地震数量占城市总数的80.4%,且有较多重要城市发生直下型地震的危险性较高(朱刚等,2007).此外,随着我国重特大工程的建设向深地、深海和中西部地区的转移,土木工程结构不可避免的需要跨越活动断层(谢礼立等,2021).着眼城市建设和重大工程中对抗震韧性的需求,虽然已有多种提高结构韧性的减震构件和设计方法(Wang et al.,2023; Zhao et al.,2023; Hu et al.,2023),但是从地震动输入的角度来看,如何合理地预测和模拟近/跨断层地震动,对提高土木工程结构的抗震设防能力,实现城市和基础设施的抗震韧性具有重要的意义和价值.

考虑到建设韧性城市和韧性交通系统中对潜在强震发生区域近/跨断层地震动场模拟的迫切需求,本文基于Graves和Pitarka提出的GP法(Graves and Pitarka,2004,2010,2015,2016;Pitarka et al.,2020,2022),使用在震前可勘测的活动断层的走向、倾角、滑动角以及可进行参数分析的震级、震源深度、断层破裂长度与宽度等震源参数建立断层破裂模型,然后使用“Seismic Waves,4th order”(简称SW4)有限差分计算程序(Petersson and Sjögreen,2018)进行宽频带地震动模拟,重建门源区域的强地震动场.通过观察分析模拟结果的衰减规律、模拟结果的峰值(PGA、PGV)和烈度分布与实测地震动记录对比分析、不同方位台站的观测记录和模拟记录的波形对比,考察和验证在中国区域内基于GP法的宽频带地震动模拟在震前预测和估计地震动场的能力,为近/跨断层工程结构的规划和设计提供参考和科学依据.

1 地震动模拟方法

地震的观测记录常被用来进行结构地震反应分析,但现有的地震观测记录只对地震动场的一小部分进行了采样,地面运动记录通常需要修正或按比例缩放,以此来适应目标场地的地震动特性,采用震后观测记录对地震进行预测和评估仍具有一定的局限性.

现如今,谱元法、有限差分法等方法被广泛应用于宽频带地震动模拟中,要实现宽频带地震动模拟,受到断层破裂模型的频带范围、速度结构模型的精细程度和网格尺寸的影响.本文使用SW4有限差分法软件和具有宽频特性的GP断层破裂模型,通过控制网格尺寸不超过最小波长的1/8来实现0～10 Hz的宽频带地震动模拟.SW4是一种四阶精度的有限差分计算程序(Sjögreen and Petersson,2012; Petersson and Sjögreen,2017a,b),它具有模拟3D材料属性、各向异性、使用曲线网格近似地表地形(Petersson and Sjögreen,2014)、沿地层深度细化网格(Wang and Petersson,2019)、设置吸收边界条件最小化来自远场边界的人工反射波(Petersson and Sjögreen,2015)的功能.

2 门源地震近场区域地震动模拟参数的确定

根据前两节的介绍,本文研究工作的开展需要确定断层破裂模型、三维地形地质结构等参数,下面逐一介绍.

2.1 断层破裂模型

在门源地震发生后,多家相关机构迅速发布了这次地震的震源参数(GFZ,2022;CENC,2022;USGS,2022).本文选用美国地质调查局(USGS,2022)所提供的震源参数,断层破裂面的走向、倾向和滑动角分别为104°、88°和15°,震中位于(37.828°N,101.290°E),震源深度为13 km.本次地震为走滑型地震,根据Wells统计的5.6

log(RLD)=-2.57+0.62×MW,

(1)

log(RW)=-0.76+0.27×MW.

(2)

经计算求得MW6.6地震断层破裂的地下长度和宽度分别为33 km和11 km.考虑到断层的长宽模拟值和最大滑移量应与USGS发布的值尽可能相近,本文模拟的断层地下破裂长度为28 km,破裂宽度为8 km,最大滑移量为431 cm.综上所述,根据USGS所提供的震源数据,作为GP法生成断层破裂模型的参数,震源模型参数汇总见表1.

表1 门源地震震源模型参数

在GP法中,定义整个断层面的平均上升时间与地震矩的关系如式(3)所示(Graves and Pitarka,2016).

(3)

其中τ为上升时间,αT是与震源机制相关的系数,rc为震源函数上升时间系数,M0为地震矩.在GP法中,默认的rc为1.6,对应的上升时间为0.72 s.

考虑到GP法是由美国西部地区的数据统计得出,而中国地区的地震动加速度反应谱相比于美国西部的地震动加速度反应谱在长周期(T>1.0 s)内幅值小(张齐,2016;张辉,2021),同时USGS提供的门源地震的震源函数上升时间也显著大于0.72 s,所以在进行地震动模拟时需增加断层破裂模型的上升时间.因此,为了研究rc对计算结果的影响,增加上升时间为4.2 s(USGS提供门源地震的震源时间函数的上升时间) 时的震源模型和上升时间为2.0 s时的震源模型.即本文将考虑三组上升时间,分别为0.72 s、2 s和4.2 s,所对应的rc为1.6、4.5和9.0.断层破裂模型是使用南加州地震中心(The Southern California Earthquake Center,SCEC)开发的宽频带平台(Broadband Platform,BBP)生成(Maechling et al.,2015),详细的内容可以查看网站https:∥github.com/SCECcode/bbp/wiki.本文基于BBP的现有代码,在仅改变参数rc的情况下生成断层破裂模型,断层的具体参数设置如表1所示,同时设置震源位于断层破裂面的几何中心,子断层的大小为0.1 km×0.1 km,BBP生成的断层破裂模型给出了每个子断层的错动速率时程,进而可计算断层平面上滑移分布、起始破裂时间与滑移率等,断层破裂模型文件详细的内容和格式可以查看网站Standard Rupture Format-SCECpedia.在其他震源模型参数不变的情况下,改变rc时断层滑移分布、起始破裂时间均不变,断层滑移分布、起始破裂时间与滑移率如图1所示.

图1 断层破裂随机滑移分布图(a); 起始破裂时间(b); 图(c)—(e) rc值依次为1.6、4.5、9.0时断层的峰值滑移率

2.2 地震动模拟研究区域

本文分为区域1、区域2、区域3三个区域进行宽频带地震动模拟,区域范围和尺寸如表2和图2所示.受计算能力所限,区域1中难以开展(0～10 Hz)地震动模拟,主要用于模拟速度场和估计区域尺度的烈度分布,地震动模拟的频带范围为(0～2 Hz).区域2用于模拟PGA和PGV随震中距的衰减分布规律,地震动模拟的频带范围为(0～10 Hz).区域3用于模拟近断层区域宽频带地震动场,对比分析近场区域台站观测记录和模拟记录的波形和峰值,地震动模拟的频带范围为(0～10 Hz).本文使用SW4软件进行地震动模拟时,需要用到输入脚本文件、地形模型、速度结构模型和断层破裂模型.输入脚本文件内的代码在计算时会读取地形模型、速度结构模型和断层破裂模型,同时会定义计算区域的原点位置、计算区域的总尺寸、网格尺寸、模拟的总时间、模拟网格点的输出等关键运行参数,如表2所示.

图2 门源地震震中位置、研究区域地形分布图

表2 SW4关键运行参数

2.3 门源地区的地形模型与速度结构模型

地壳介质的波速、密度和品质因子等特性在空间中存在显著的差异,建立和使用合理的地壳速度结构模型是开展地震动模拟的基础.本次门源地震发生在青海省门源县,在本文研究区域内(36.6°N—39.0°N,100.0°E—103.0°E),使用美国地质调查局发布的具有30 m精度的全球地形数据STRM1(Farr et al.,2007)建立研究区域内的地形模型.使用模型USTClitho2.0(Han et al.,2022)来建立速度结构模型,由于模型USTClitho2.0中没有地下岩层的密度,使用全球地壳模型CRUST1.0(Laske et al.,2013)中的地下岩层密度进行补充.P波与S波的深度切片分布图如图3所示,P波与S波竖向剖面分布图如图4、图5所示.地震波的非弹性衰减采用粗粒度方法,品质因子根据Pitarka和Graves给出的经验关系公式计算(Pitarka and Graves,2016):

图3 P波与S波随深度变化切片分布

图4 P波在垂直方向上变化切片分布

图5 S波在垂直方向上变化切片分布

QS=vS×50,

(4)

QP=2×QS,

(5)

其中,QS为P波品质因子,QP为S波品质因子,vS为S波波速.

3 计算结果与讨论

根据前文介绍的模拟方法、震源参数、地形结构模型与地壳结构模型,按照上升时间0.72 s、2.0 s、4.2 s,将计算工况分为三组,对应的rc分别为1.6、4.5、9.0,对目标区域的地震动场开展强地震动模拟.为了验证本文计算模型的合理性,选取区域2内观测点的模拟记录的东西向、南北向的PGA和PGV,与俞言祥等(2013)编制建立的我国青藏地震区基岩场地水平地震动峰值预测方程进行对比,地震动峰值预测方程的计算公式为

lgY=A+BM+Clg(R+DeEM)±σ,

(6)

其中Y代表地震动加速度和速度峰值,A、B、C、D、E代表回归分析得到的系数,M、R和σ分别代表面波震级、震中距和标准差.预测方程与PGA模拟记录对比时,预测方程的回归系数A、B、C、D、E分别取值为2.457、0.388、-1.854、0.612、0.457,M和σ分别取6.9和0.236;预测方程与PGV模拟记录对比时,预测方程的回归系数A、B、C、D、E分别取值为0.443、0.474、-1.696、0.612、0.457,M和σ分别取6.9和0.271.图6绘制了基于不同rc的模拟记录峰值地面加速度(PGA)和峰值地面速度(PGV)随震中距的衰减规律,其中黑色实线代表预测中位值,黑色虚线表示±1倍标准差.

图6 基于不同rc的模拟记录 PGA 和 PGV 随距离的衰减分布

由图6可知:三组不同rc的模拟记录和预测中位值的PGA与PGV随距离的变化总体衰减速率基本一致.当rc为1.6时,模拟结果PGA和PGV较预测中位值偏大.当rc为4.5和9.0时,模拟结果PGA和PGV总体上与预测中位值较符合.随着rc的增加,在相同震中距的情况下,三组模拟结果的PGA和PGV均逐渐降低.这是因为周期为1～2 s的地面运动受rc的影响显著,即rc的值越大,上升时间越大,峰值滑移率越小,地震动速度和加速度的幅值越小.综上所述,当rc为4.5与9.0时,模拟结果与预测中位值符合较好.

表3给出了距离震中最近的两个台站C0028、C0029的观测记录、模拟记录在10 Hz低通滤波后的PGA和PGV.从表3可以看出,当rc为1.6时,台站C0028模拟记录的PGA和PGV在南北向较观测记录偏大,东西向和竖向与观测记录较符合;台站C0029模拟记录的PGA和PGV在三个方向与观测记录较符合.当rc为4.5和9.0时,台站C0028与台站C0029模拟记录PGA与PGV的三个方向与观测记录均较符合.综上所述,当rc为4.5与9.0时,模拟结果与观测记录符合较好.

表3 台站C0028、C0029观测与模拟记录的PGA、PGV

根据《中国地震烈度表》(国家市场监督管理总局和国家标准化管理委员会,2020),计算研究区域地震烈度.结合模拟所得加速度时程与速度时程,计算研究区域地震烈度.本次门源震后国家预警工程自动产出结果(王士成等,2017)见图7a,rc为1.6、4.5和9.0对应的研究区域的地震烈度空间分布分别由图7b、7c和7d表示.图8为不同rc参数下PGA和PGV地表分布图.

图7 国家预警工程自动产出烈度图与地震动模拟产出烈度图

图8 不同rc参数下PGA和PGV相应模拟量的地表分布图

从图7a中可知,断层附近区域极震区的最大烈度可达9度,本文模拟的三幅图与国家预警工程的最大烈度Ⅸ相一致.对比图7a和图7d可知,图7d高烈度区(Ⅵ、Ⅶ、Ⅷ、Ⅸ度区)分布的经纬度范围与图7a较符合,具体的经纬度范围如表4所示;图7a中高烈度区在平行和垂直于断层方向的烈度较其他方向衰减慢,与图7d中高烈度区呈垂直断层和平行断层分布特征相似;从表5所示的各烈度区面积可以看出图7d中高烈度区与图7a符合度较好.

表4 自动产出结果与地震动模拟结果Ⅵ、Ⅶ、Ⅷ与Ⅸ度区范围

表5 自动产出结果与地震动模拟结果Ⅶ、Ⅷ与Ⅸ度区面积

从表5和图7中可以看出,随着rc的增大,沿着断层走向的高烈度区面积逐渐减小,垂直于断层的高烈度区面积基本不变.当rc为1.6时,模拟结果的Ⅶ度区面积大于自动产出结果,Ⅷ、Ⅸ度区面积远大于自动产出结果;当rc为4.5时,模拟结果Ⅶ度区面积与自动产出结果吻合较好,Ⅷ度区面积远大于自动产出结果;当rc为9.0时,与自动产出结果相比,模拟结果Ⅶ、Ⅷ度区面积略偏大.综上所述,随着rc的增加,模拟的各烈度区域面积整体呈现减小趋势,这是由于rc对1～2 s周期地面运动的影响显著,即rc的值越小,速度幅值越大,烈度区的面积越大.

当rc为9.0时,模拟的各烈度区面积与国家预警工程较符合,但是Ⅶ、Ⅷ度区面积较自动产出记录仍然略有偏大.一般来说,台站多布置在地形较为平坦的区域,而地震动幅值在山地区域会有显著放大,从而导致了模拟结果的Ⅶ、Ⅷ度区面积偏大.图9绘制了台站C0029附近区域的地形图和南北方向速度分布云图.可以看出,图中北部和东北部山地区域出现了显著的地震动放大效应.考虑到模拟区域山区分布广泛,本文模拟的Ⅶ、Ⅷ度区面积较自动产出列度结果偏大.

图9 台站C0029周围地形图(a)与台站C0029附近区域南北方向PGV(b)

当rc=9.0时,台站C0028和C0029在10 Hz低通滤波下的观测记录波形和模拟记录波形如图10所示.其中图10a为台站C0028观测记录和模拟记录的加速度波形对比图;图10b为台站C0028观测记录和模拟记录的速度波形对比图.图10c为台站C0029观测记录和模拟记录的加速度波形对比图;图10d为台站C0029观测记录和模拟记录的速度波形对比图.从图10中可以看出,除台站C0029在南北向模拟记录的速度峰值较小,两台站观测记录和模拟记录的加速度和速度在峰值和波形上均较为符合.

图10 台站C0028和C0029在10 Hz低通滤波下的观测记录波形和模拟记录波形

综合模拟记录峰值与青藏地区的地震动峰值预测方程、强震观测记录和震后国家预警工程自动产出烈度区分布的对比结果,同时对比分析台站C0028和C0029在10 Hz低通滤波下的观测记录波形和模拟记录波形.可以看出,当rc为9.0时,地震动模拟结果与观测记录的符合度较好.

4 结论

本文采用Graves和Pitarka断层破裂模型(GP法)和“Seismic Waves,4th order” (SW4)有限差分计算程序,对2022年门源MS6.9地震开展了强地面地运动模拟,开展了以下工作:(1)对比了模拟记录峰值与青藏地区的地震动峰值预测中位值随震中距的衰减规律;(2)将C0028和C0029两个台站的模拟记录与观测记录的PGA和PGV进行对比分析;(3)根据模拟记录给出烈度分布图,并与国家预警工程自动产出烈度图进行对比分析;(4)对比分析C0028和C0029两个台站的模拟记录与观测记录的波形.可以得到如下结论:

(1) 对于三种不同上升时间系数(rc)的震源模型,rc为4.5和9.0时模拟记录的峰值(PGA、PGV)与青藏地区的地震动预测方程的预测中位值符合程度较好.

(2) 在近断层区域宽频带地震动场的计算结果中,对于台站C0028和C0029,rc为4.5和9.0时模拟记录与强震观测记录的PGA和PGV基本吻合.

(3) 本文模拟的各烈度区面积随着rc的增大而减小.当rc为9.0时烈度区走向、范围和烈度区面积与国家预警工程自动产出结果基本吻合.

(4) 当rc为9.0时,台站C0028和C0029的模拟记录和观测记录的加速度和速度在峰值和波形上均较为符合.

(5) 综上所述,本次强地震动模拟所得结果较好的重现了门源地震,在合理选取rc参数的情况下,GP法较好的重现了门源地震动场,但是影响强地震动模拟的因素还很多,未来还需开展破裂速度、震源时间函数等参数的影响和基于中国大陆强震观测记录的震源破裂模型研究.

致谢感谢中国地震局工程力学研究所“国家强震动台网中心”为本研究提供数据支持.感谢审稿专家和编辑对文稿提出了宝贵的修改意见和建议.同时感谢开源绘图软件GMT(Wessel et al.,2019)为本文提供了方便高效的绘图平台.