基于FK方法和GP14.3震源模型的2023年土耳其Mw7.8级地震宽频地震动合成

巴振宁,韩书娟,赵靖轩,刘 悦,芦 燕,陈三红

(1. 中国地震局地震工程综合模拟与城乡抗震韧性重点实验室(天津大学),天津 300350; 2. 天津大学 土木工程系,天津 300350;3. 中国地震灾害防御中心,北京 100029)

0 引言

北京时间2023年2月6日9时17分土耳其发生Mw7.8地震,震源深度17.9 km,最高烈度为Ⅸ度,震中(北纬37.17°、东经37.03°)位于东南部的东安纳托利亚断裂带(east anatolian fault zone,EAFZ)[1]。东安纳托利亚断裂带位于安纳托利亚板块(亚欧板块内部的次级板块)和阿拉伯板块的交界处如图1所示,是土耳其境内的一条主干左旋走滑断裂带,起于西南的哈塔伊(Hatay),终止于东北的卡利若瓦(Karliova)三联结点位置[2]。东安纳托利亚断裂带地震活动频繁,20世纪以来多次发生6级以上强震[3]。同日下午18时24分,再次发生Mw7.5地震,震源深度10 km,最高烈度亦为Ⅸ度,震中位于首次地震震中东北部96 km处[1]。这两次地震构成了典型且罕见的“双震”地震序列[4]。

图1 土耳其地震的震中及其相关的构造环境Fig. 1 Epicenter of the Turkey earthquake and its associated tectonic environment

此次地震震级大、震源浅和地面运动强度大,造成了惨重的人员伤亡和严重的建筑结构破坏:41156人死亡、105505人受伤、超过8.4万栋建筑倒塌和严重受损或急需拆除[5]。其中:位于土耳其东南部的人口密集的工业城市-加济安泰普[6],距离此次Mw7.8地震震中仅37 km,是此次地震中受影响最为严重的城市之一。针对该Mw7.8地震,本文采用作者提出的FK方法开展0～10Hz的近断层宽频地震动合成[7],并重点分析加济安泰普地区的地震动时空场分布特征,以期为该地区的震害原因分析和地震区划等工作提供参考。

1 地下结构与震源模型建立

确定性物理宽频地震动模拟依赖于能激发宽频地震波的合理震源模型和能模拟宽频地震波传播的精确高效方法[8]。震源模型采用的是GP14.3运动学混合震源模型[9],该模型由Graves和Pitarka(2010)宽频带地面运动模拟方法(GP2010)改进后得到(GP14.3为在当前南加州地震中心宽带平台上的更新版本名称)[10],其能够同时考虑破裂面上低波数的确定性错动和高波数的随机性错动,具备激发宽频地震波场的能力。地震波传播模拟方面采用作者提出的FK法,该方法是一种具有严密理论基础的地震动合成方法,其适用于跨尺度一维速度结构宽频带地震波的传播模拟。本节详细介绍一维波速结构和运动学混合震源的构建。

1.1 加济安泰普地区地下一维速度结构模型的构建

考虑到后续分析主要针对加济安泰普市,地下一维速度结构基于该市进行构建,相应参数表1,对应的示意图如图2所示。

表1 加济安泰普地区的地下一维速度结构Table 1 1D velocity structure in Gaziantep region

图2 加济安泰普地区的一维速度结构模型示意图

该一维地下结构模型包括3部分:1)依据USGS网站提供的Vs30数据[11],确定0～0.03 km深度范围内的一维速度结构参数,其中剪切波速Vs取为加济安泰普市界范围内Vs30的平均值。2)依据INFANTINO等[12]给出的伊斯坦布尔地区波速模型确定方法,确定0.03～2 km及2～5 km深度范围内的一维速度结构参数:在0.03～2 km之间,Vs随深度z的增加以(z)0.5的速率增加;在2～5 km之间,假定Vs呈线性增加,通过对更密集的深度点进行线性插值、分层取均值得到每层土体的Vs。3)依据全球地壳模型Crust1.0,确定大于5 km深度范围内的一维速度结构参数:包括上地壳(5～17.64 km)、中地壳(17.64～26.54 km)、下地壳(26.54～36.18 km)以及地幔层(>36.18 km)的Vp、Vs和密度。在整个地下一维速度结构模型中,除Crust1.0直接给出的大于5 km深度范围的Vp和密度之外,Vp、密度和品质因子的取值方式如下:Vp的值为1.6Vs[12],土层密度的取值参考文献[13],S波和 P波的品质因子分别由Qs=Vs/10和Qp=Vs/5(Vs的单位为m/s)确定[14]。

1.2 GP14.3运动学混合震源模型的构建

本文采用能同时考虑破裂面上低波数确定性错动和高波数随机性错动的GP14.3运动学混合震源模型,该模型的建立过程如下:

1)基于USGS发布的震源机制(见表2)和断层滑动量分布结果(见图3)[15],建立低波数确定性的凹凸体震源模型。震源参数包括全局震源参数、局部震源参数和其他震源参数[16],其中:全局震源参数表示断层的形状、大小和位置等,如断层走向、倾角和滑动角,断层面沿走向长度、沿倾向宽度,断层顶面深度,断层面平均滑动量,断层面划分的矩形子源尺寸等;局部震源参数表示震源内部的滑动不均匀性,如凹凸体的面积、凹凸体的滑动量等,这些参数由断层的滑动量和分布结果来确定;其他震源参数表示破裂的形成和终止,如破裂起始点和破裂传播方式,这些参数由破裂前端分布结果来确定,其中:破裂传播方式为由破裂起始点向四周传播。以上参数详见表3。

表2 USGS发布的震源机制参数Table 2 Focal mechanism parameters given by USGS

表3 土耳其Mw7.8地震的震源参数Table 3 Local source parameters of the Turkey Mw7.8 earthquake

图3 USGS所给的断层破裂面上的滑动量及破裂前端分布(色谱表示滑动量幅值,五角星表示破裂起始点,灰色箭头表示滑动矢量,等值虚线表示破裂起始时间(以秒为单位))Fig. 3 Distribution of the slip and rupture-front on the fault rupture surface released by USGS (Slip amplitude is shown in color. The point of rupture initiation is denoted by a star. The gray arrow represents the slip vector. The equivalent dotted lines show the rupture initiation time in seconds)

2)在低波数确定性的凹凸体震源模型中引入高波数随机性错动,确保合成宽频地面运动的带宽有效性,进而构建GP14.3运动学混合震源模型[17]。①首先将已有的凹凸体断层面上滑动量的空间分布经二维傅里叶变换至波数域,得到断层面上确定性的滑动波数谱。②其次采用波数衰减满足von Karman自相关函数的波数谱并引入随机数φ=[-π, π]表达断层破裂面上小尺度的随机变量。③然后在波数域中结合确定性的低波数谱和随机高波数谱,并利用二维逆傅里叶变换至空间域中。

以上引入高波数随机性错动的过程中所涉及的Hurst指数H、控制结合的锐度N、沿断层走向和倾向的拐角波数kcx与kcd均参考文献[9]进行取值,进而确定断层面上的滑动量分布,如图4(a)所示。关于破裂前端的分布,首先确定破裂速度的初始分布,根据子源位置、初始破裂速度分布以及各子源滑动量对破裂起始时刻的时间扰动,得到破裂起始时刻的空间分布(图4(a)等值线表示间隔5 s的破裂前端分布)。上升时间表征子源完成破裂过程所需时间,本文采用Hartzell矩率函数,得到各子源的上升时间后根据断层面上的平均上升时间按比例调整各子源上升时间,可得其在断层面上的空间分布如图4(b)所示。

图4 运动学混合震源Fig. 4 Hybrid kinematic source model

2 模型的适用性验证分析

本节通过将地震动合成结果与8个台站的观测记录对比,有效检验了宽频地震动模拟的可靠性和模型参数(包括地下一维速度结构模型和混合震源模型)的适用性,且FK法的精度及高频稳定性本文作者在文献[18-19]中已广泛验证。首先,基于建立的GP14.3运动学混合震源模型,采用 FK 法计算了距离此次Mw7.8地震的震中50 km范围内地震动较强的8个台站的加速度时程(台站位置信息见图5和表4);其次,将地震动合成结果和AFAD提供的强震记录及其反应谱进行对比[20]。该次计算的时间步距取为 0.01 s,时间步的数量为8192。这里需要说明的是:本文所选取的台站记录来源于AFAD网站截止至2月21日的强震记录数据。

表4 所验证的8个台站的信息

图6对比了本文计算的8个台站的三分量加速度时程和强震记录,黑线代表强震记录,红线代表本文合成结果,每条加速度时程的时长均为80 s,PGA均标注于曲线的右上方,单位为cm/s2。强震记录及本文计算结果均经4阶Butterworth 低通滤波至 10 Hz。从地震动时程比较可知:前6个为记录完整的台站,合成地震动在4616 台站偏小,在其它5个台站符合较好;NAR和4615台站记录虽然存在部分缺失,但因其记录到了主要波形成分而仍具有一定的参考意义,本文合成的两个水平方向的地震动与观测记录前段呈现的幅值和波形特征基本上一致,但竖向的地震动比观测记录明显偏小;除4630台站的EW方向、4616台站的NS方向、NAR和4615台站的UD方向的上的计算结果和强震记录的PGA相差较大之外,其他的计算结果和强震记录符合较好;2718、4629和4632台站处的计算结果和强震记录在强震段持时符合较好,但其他台站处有一定差别。总体上,本文计算的加速度时程与强震记录有一定的一致性。

图7对比了本文合成结果与强震记录的5%阻尼比加速度反应谱。为便于比较和观察,将所有EW分量缩小5倍以及UD分量缩小25倍。由图7可知:NAR和4615台站处的竖向合成结果的短周期段小于观测记录,NS和EW方向的合成结果在0.1～10 s内均与观测记录符合较好;除NAR和4615台站外的其他6个台站的三个方向的合成结果在短周期段内均与观测记录符合较好;而2703、4629、4630和4632这4个台站的三个方向的合成结果在长周期段内高于观测记录。综上所述,从平均意义上来看:计算的地震动的短周期段与观测记录符合较好,而长周期段有一定程度的高估现象。根据上述对地震动时程、峰值和反应谱的比较分析表明:本文所构建的一维速度结构模型和混合震源模型适用于2023年土耳其Mw7.8地震的宽频地震动模拟,可进一步用于近场特性的分析。

图7 模拟地震动反应谱与强震记录对比Fig. 7 Response-spectrum comparison between the simulations and the records

3 加济安泰普及附近区域地面运动场合成及分析

本节对土耳其Mw7.8地震进行了三分量宽频地震动合成。本节设定研究目标区域为加济安泰普及附近区域,并在上一节研究基础上合成了该区域内的地面加速度峰值、速度峰值以及烈度分布。首先,设定目标区域:如图8所示,设定东西向180 km×南北向150 km的目标矩形区域(红色矩形,面积为27000 km2),以3 km间隔均匀设置61×51,共计3111个计算点;其次,采用本文方法及模型合成了目标区域地面运动场;最后,给出了目标区域的三分量加速度和速度峰值分布,分析了其空间分布特征及近断层效应;给出了目标区域的烈度分布图,并与USGS给出的土耳其Mw7.8地震的ShakeMap进行对照;给出了部分位置的时程结果,分析此次地震造成的地面永久位移及近断层速度脉冲特征。

图8 模拟区域布置示意图(黑色矩形表示断层投影面、星号表示破裂起始点)Fig. 8 Sketch of the simulation area (The black rectangle indicates the projection plane of the fault and the star indicates the rupture initiation point)

3.1 地面运动峰值分布

图9给出了目标区域地面加速度峰值(PGA)分布和速度峰值(PGV)分布图,其中从左到右依次是NS、EW和竖向(UD)分量,从上到下依次是PGA和PGV。图中矩形和星号分别表示断层面的地表投影范围和震中位置,蓝色虚线区域为加济安泰普市,绿色虚线区域分别为加济安泰普市区区域和伊斯拉希耶镇区域,黑色虚线表示主要道路。

图9 地震动场PGA和PGV分布Fig. 9 Distribute of PGA and PGV

从图9中可以看:加济安泰普的市区范围内的NS向PGA接近180 cm/s2,PGV接近100 cm/s,说明加济安泰普市区范围内的地震动强度较小。分析其原因为:走滑型断层破裂释放的能量主要集中在破裂前方区域,垂直断层走向两侧的区域聚集能量较小,且远离断层投影区域后地震动幅值衰减较快,该市区虽距震中十分近,但其位置在这条西南-东北走向的东安纳托利亚断裂带的东南方向,不属于破裂前方区域,故该市区的地震动强度较小,这也与报道中提及的“加济安泰普市区极少建筑出现崩塌或明显损毁,大部分楼房从外部都看不出明显损伤”的现象相符合[21]。由图9还可知:位于加济安泰普的西南部、直接处于地震断层带上的伊斯拉希耶镇(37.03°N,36.63°E)地区内的EW向PGA高达560 cm/s2,PGV高达150 cm/s,这与附近2718台站的观测记录(EW向PGA=584.2 cm/s2,EW向PGV=113.16 cm/s)较为接近,说明该地区的地震动强度较大。此外,从图9的空间分布样式来看:峰值较大区域集中在断层投影面附近,这与震源模型的凹凸体位置相对应,体现了凹凸体对地面运动峰值的主导控制作用以及近断层地震动的集中性效应。

3.2 烈度分布

在地震发生的初期,通常可利用烈度与强地面运动参数之间的转换关系来快速获取地震烈度分布特征,以快速可靠地估计烈度分布。本节基于强地面运动合成结果,应用WALD等[22]给出的加速峰值(PGA)、速度峰值(PGV)和烈度(IMM)的经验关系式(1),给出震区的烈度分布图。

(1)

式中:PGA单位为cm/s2,PGV的单位为cm/s,IMM为场地的烈度。基于上式即可将得到的地表水平PGA/PGV转换为相应的地震烈度值。

图10(a)和图10(b)分别为USGS给出的土耳其Mw7.8地震的ShakeMap和本文基于强地面运动合成结果得到的地震烈度分布图[23]。将图10(a)范围缩小至本文目标区域的大小,以便于本文结果与其对比。从图10(b)中给出的烈度分布结果可看出:加济安泰普市区范围内的地震烈度为VII度,与图10(a)中加济安泰普市区范围的烈度一致;处于西南部的伊斯拉希耶镇地区的地震烈度高达IX度,与图10(a)中伊斯拉希耶镇地区的烈度一致,且与图9所反映的伊斯拉希耶镇地区的地震动强度较大的结论相符。

图10 烈度分布Fig. 10 The intensity distribution

这里需要说明的是:本文所构建的GP14.3运动学混合震源模型没有考虑弯折,为简单连续的平面断层,而USGS最终更新的震源模型是包含三段的弯折断层,这也是造成本文一些地震动合成结果与台站观测记录存在一定差距的主要原因。

3.3 时程结果分析

据报道,加济安泰普一西部城镇-努尔达伊附近出现几百米长的地表断裂,水平位移约3～4 m[24]。本节首先计算努尔达伊镇中心(37.18°N,36.74°E)地表上一观测点的时程结果,如图11所示,其中从左到右依次是平行走向、垂直走向和竖向分量,从上到下依次是加速度、速度和位移时程。从图11所合成的时程结果中可发现:在水平方向上有比较明显的永久位移:平行于断层走向的地表永久位移达到了105 cm,垂直于断层走向的地表永久位移达到了100 cm,而竖向的地表永久位移约为10 cm。合成结果所反映的现象说明近断层效应较为明显,且这种现象可能导致建筑发生更加明显的破坏。

图11 Mw7.8土耳其地震中努尔达伊地区的三方向的加速度、速度和位移时程Fig. 11 Time histories of acceleration, velocity and displacement in three directions of Nurdagi region in Mw7.8 Turkey earthquake

图12分别给出了在地表上沿断层西南-东北走向以10 km的间隔均匀设置的15个观测点的位置以及观测点平行断层走向、垂直断层走向两个方向的速度时程。如图12(a)所示,R1观测点位于破裂起始点在地表的投影处,R1～R15分别对应震中距为0、10、20、30、40、50、60、70、80、90、100、110、120、130和140 km的各个观测点。图12(b)为本文所合成的平行于断层走向的速度时程结果,在速度时程中可看到单向速度脉冲(特别是R6～R9这四个观测点);图12(c)为本文所合成的垂直于断层走向的速度时程结果,可看到明显的双向速度脉冲,这是因为此次地震的断层为倾角较大的走滑断层,方向性速度脉冲主要体现在垂直断层走向的分量上[25],特别是位于破裂前端的R6～R8这三个观测点,满足方向性效应的条件,双向速度脉冲现象更为明显。由图 12(b)和图12(c)可以看出:在破裂起始点于地表的投影处(R1观测点)速度脉冲现象不明显,在距离断层破裂起始点一定的距离后开始看到较为明显的速度脉冲,比如距断层破裂起始点30 km处的R4观测点,速度脉冲已经很明显。

图12 震中距分别为0、10、20、30、40、50、60、70、80、90、100、110、120、130和140 km的各观测点Fig. 12 Observation points with epicenter distance of 0、10、20、30、40、50、60、70、80、90、100、110、120、130 and 140 km

4 结论

针对2023年2月6日土耳其南部Mw7.8地震开展了近断层宽频地震动合成研究。综合考虑Vs30数据及全球地壳模型Crust1.0,构建了土耳其东南部城市-加济安泰普及附近地区的地下一维速度结构模型。依据USGS发布的信息建立了GP14.3运动学混合震源模型,进而采用作者提出的FK方法合成了8个台站的三分量加速度时程,通过与台站的时程记录及其反应谱对比,验证了本文方法的可靠性及模型的适用性。最后计算并分析了加济安泰普附近27000 km2目标区域内的地震动场,得到如下结论:

1)本文合成的加济安泰普市区范围内的NS向PGA接近180 cm/s2,PGV接近100 cm/s;本文合成的位于加济安泰普的西南部和直接处于地震断层带上的伊斯拉希耶镇地区的EW向PGA高达560 cm/s2,PGV高达150 cm/s,与附近2718台站的观测记录一致。

2)此次土耳其Mw7.8地震呈现明显的近断层地震动集中性效应、地面永久位移以及破裂方向性效应。

3)基于强地面运动合成结果得到的地震烈度分布显示:加济安泰普市区范围内的地震烈度为VII度;处于西南部伊斯拉希耶镇地区的地震烈度高达IX度。

4)从努尔达伊地区的地震动合成结果中可发现在水平方向上有比较明显的永久位移。

由于本文震源模型未考虑实际弯折,同时采用的方法也不能考虑起伏地表的影响和地下结构的横向非均匀性,这些因素可能是造成本文模拟结果与强震记录在某些台站差异还较大的主要原因。因此,在接下来的地震动模拟研究中,有必要对震源进行深入研究并考虑复杂地形的影响,以建立更为可靠的模型,进行更精细化的复杂场地地震动模拟,为地震区划和抗震设防等工作提供参考。