引言

2022年1月8日01时45分27秒青海省海北藏族自治州门源回族自治县发生M_S6.9地震，震中位置为（37.77°N，101.26°E），震源深度为10 km。本次地震造成了17 069人受灾，由于震中距离人口稠密地区较远，尚无人员伤亡的相关报道（青海日报，2022），兰新高铁浩门至山丹军马场区间隧道群由于此次地震而发生局部塌方（央视网，2022）。门源M_S6.9地震发生在青藏高原东北缘冷龙岭断裂、托莱山断裂和肃南—祁连断裂的阶区部位。现场工作队在冷龙岭断裂带西段探测到长约22 km的地表破裂带（青海省地震局，2022）。托莱山断裂和冷龙岭断裂现为左旋兼挤压断裂（李强等，2013），历史上无M7以上地震记载（姜文亮，2018），最近一次强震为2016年1月21日门源M_S6.4地震，该地震的发生反映了青藏高原地块向NE向不断推挤生长的过程（胡朝忠等，2016）。门源M_S6.9地震发生后，许英才等（2022）对门源地震早期序列（2022年1月8日至12日）进行了重定位和震源机制研究，其结果表明目前门源地区还存在一定的应力积累且应力尚未得到充分释放，该地区仍有发生强震的危险。由于地震灾害主要是地震在地表产生的强地面运动造成的，因此为了获得地震波场传播过程及其引起的地表响应，对门源M_S6.9地震开展强地面运动初步模拟及烈度估计具有重要意义。

地震波场的正演模拟需要考虑复杂地表。曲线网格有限差分方法（Zhang，Chen，2006；Zhang et al，2012）适用于含起伏地形的地震波场模拟，该方法中地表的形状用任意曲线网格近似，并采用牵引力镜像法处理自由表面条件。曲线网格有限差分的强地面运动模拟方法在地震后的灾害评估中得到了非常广泛的应用。Zhang等（2008）结合三维介质模型、震源破裂模型和地表地形数据模拟了2008年汶川M_S8.0地震的强地面运动，研究表明断层破裂方式和盆地构造控制了地表峰值速度（peak ground velocity，缩写为PGV）的分布，地表起伏剧烈地区往往对应较大的PGV数值，应关注其震害问题。张振国等（2014a，b）对2014年2月12日新疆于田M_S7.3地震和2014年8月3日云南鲁甸M_S6.5地震引起的强地面运动作了初步模拟和烈度预测，研究显示地震动在山峰、山脊处具有较大幅值，该结果可指导震后的灾区重建工作。赵宏阳和陈晓非（2017）利用1975年2月4日辽宁海城M_S7.3地震的地震地质资料，模拟计算了海城地震的波场传播过程，分析了强地面运动的方向性效应、盆地效应和近断层效应，得出理论烈度分布同震后调查烈度分布基本一致，验证了研究所用的震源模型和速度结构的合理性。

本文拟根据张勇^①提供的门源地震震源破裂过程的初步结果，利用曲线网格有限差分方法模拟门源M_S6.9地震的强地面运动，再结合地表峰值速度和烈度间的关系，计算地震烈度，并在此基础上评估地震灾害分布特征，以期为门源地区的防震减灾提供科学依据。

1.   研究区概况及计算方法和模型

研究区范围如图1所示，可以看出，该区域的地貌特征变化大，地势西南高东北低，高程介于1.3—5.0 km之间。门源M_S6.9地震发生后至2022年1月20日，该地区已发生M_S≥4.0余震23次，其中M_S≥5.0余震2次。

图  1  研究区及周边的地质概况和2022年1月8—20日震源区M_S≥3.0 地震分布

Figure  1.  Geological overview of the study area and its surrounding regions，spatial distribution of M_S≥3.0 earthquakes in source zone during the period of 8 to 20 January，2022

下载: 全尺寸图片 幻灯片

采用曲线网格有限差分方法（Zhang，Chen，2006；Zhang et al，2012）对门源M_S6.9地震的强地面运动进行模拟。计算中需要设置描述地形起伏的网格模型、反映地下物质属性的介质模型，以及表示地震破裂过程的震源模型。

地形选取GTOPO30地形数据，其水平分辨率大约为1 km。整个计算区域尺度为350 km×220 km，深度为60 km。将研究区域离散成700×440×120个网格，垂直方向采用等间距排列，单位网格为边长500 m的立方体。

考虑到面波对纵波速度的灵敏度和起伏地形的影响，Han等（2022）提出一种改进的体波和面波数据联合反演方法获得中国大陆地壳和上地幔水平分辨率为0.5°的速度结构（USTClitho2.0模型）。该模型对于青藏高原地块内体波射线覆盖相对稀疏的区域，加入面波频散数据可更好地提高纵波速度和横波速度的准确度，因此本文以USTClitho2.0模型为基础建立研究区域的网格化横波速度和纵波速度模型。每个网格的东向、北向和垂直方向的长度分别为43.75 ，55 和5 km，同时依据Brocher （2005）提出的密度和纵波速度的经验关系确定各网格的密度。模拟采用的横波速度结构如图2所示，可以看出，2—12 km深度存在高速层，沿北东方向地壳厚度逐渐减薄。虽然青海湖处于研究区域内（图1），但是由于其位于模拟计算区域的边缘，且水深较浅，因此在建立介质模型时忽略其影响。

图  2  横波速度结构

Figure  2.  S-wave velocity structure

下载: 全尺寸图片 幻灯片

张勇^①反演了该地震的震源破裂过程，结果显示断层滑动出露地表，走向为103°，倾角为88°，滑动角为−6°，共有31×11＝341个子断层，子断层的空间分辨率为2 km×2 km，采样点为80个，每间隔0.25 s给出一个滑动速率，破裂过程持续时间约20 s，最大滑动量为1.8 m。基于这一结果，本文设置了相应的震源模型，如图3所示。

图  3  本文所用的门源M_S6.9地震震源模型

Figure  3.  Source model of Menyuan M_S6.9 earthquake used in this study

下载: 全尺寸图片 幻灯片

2.   强地面运动模拟结果

强地面运动模拟计算的时间步长为0.01 s，总步数为1万步，模拟时长共100 s，使用640个计算核心进行计算。门源地震x分量的模拟速度场快照如图4所示，可见：地震的主要能量由位于震中附近下方的位错产生，且以水平走向错动为主，因此在断层垂向上产生了能量较强的S波；大约第2.36 s时，地震波到达地表，在初 始破裂时刻，速度的最大值集中在断层破裂的前锋上；随着远离发震断层，地震波场能量逐渐减小。

图  4  门源M_S6.9地震x分量粒子速度波场快照

Figure  4.  Wavefield snapshots of the x component particle velocity of Menyuan M_S6.9 earthquake

下载: 全尺寸图片 幻灯片

通过强地面运动模拟得到各网格点x，y，z三个方向不同时刻的速度，然后对三个方向的速度分量求矢量和获得各网格点的运动速度随时间的变化，根据各网格点的速度最大值确定研究区域的地表峰值速度（PGV）分布，并结合国家市场监督管理总局和国家标准化管理委员会（2021）给出的PGV与烈度之间的关系得到对应的地震烈度，如图5所示，结果表明：沿平行断层走向方向的地震动衰减明显小于垂直断层走向方向；地震的最大烈度为Ⅷ度（PGV＝37.63 cm/s），位于震源破裂起始点附近区域；Ⅷ度区主要涉及门源县、祁连县和肃南县的部分区域；Ⅶ度区主要涉及大通县、永昌县、民乐县等部分区域。中国地震局（2022）野外调查的烈度分布显示等震线长轴呈WNW走向，同理论模拟显示的高烈度区主要沿WNW方向延伸的结果一致。理论烈度分布与中国地震局工程力学研究所强震动观测组（2022）给出的仪器烈度分布（图6）较为接近，即高烈度区主要沿断层走向展布，且断层南侧的地震烈度大于北侧。受模型分辨率和计算成本的限制，本文模拟的地震波场的最高有效频率为0.4 Hz （2.5 s），缺少高频成分，因此获得的烈度较强震仪和烈度仪观测的结果偏低，但相较于使用低分辨率介质模型的模拟结果（徐剑侠等，2015），地震动模拟的最高频率得到了一定的提升。

图  5  门源M_S6.9地震烈度分布图

Figure  5.  The simulated intensity of Menyuan M_S6.9 earthquake

下载: 全尺寸图片 幻灯片

图  6  门源M_S6.9地震仪器观测的烈度分布（引自中国地震局工程力学研究所强震动观测组，2022）

Figure  6.  Instrumental seismic intensity distribution of Menyuan M_S6.9 earthquake （after Strong Motion Observation Group，Institute of Engineering Mechanics，China Earthquake Administration，2022）

下载: 全尺寸图片 幻灯片

地震灾害区域集中在发震断层附近，并向WNW方向和ESE方向延伸，这主要由两个原因造成：一方面，该地震为高倾角的走滑型地震，地震波能量主要沿走向方向传播，表现出强地面运动的方向性效应；另一方面，起伏地表对地震波传播具有重要影响（Zhang et al，2008），发震断层的WNW方向和ESE方向分布有托莱山、大通山、达坂山和冷龙岭等山脉（图1），属于山地地貌，因此该区域的地震动高值可能与山脊地区地震波多次反射有关。

3.   讨论与结论

强地面运动模拟结果由地形数据、介质模型和震源模型共同决定，因此数据选择对结果的可靠性有较大影响，为使结果更为可信，本文参考美国加州综合地震破裂预测模型（Field et al，2014）中提出的“使用可获得最优解”原则对介质模型进行选择。常用的介质模型包括CRUST1.0模型（Laske et al，2012）和CRUST2.0模型（Bassin et al，2000），二者均基于水平层状介质的假设，水平分辨率分别为1.0°和2.0°，而Han等（2022）提出的USTClitho2.0模型的水平分辨率更高，因此研究中选择USTClitho2.0模型作为介质模型。本文基于门源地震震源区附近的地形数据、介质模型和震源模型，使用曲线网格有限差分方法计算了该地震的近场地震波传播过程，得到了理论的速度场快照和地震烈度分布。研究结果表明：总体上垂直断层走向方向的地震动衰减大于平行断层走向方向的地震动衰减，震中最大烈度为Ⅷ度；模拟得到的理论烈度同野外调查的地震烈度分布基本一致。受强地面运动方向性效应和起伏地表的影响，地震灾害主要沿断层的WNW方向和ESE方向分布。由于地震发生区域主要以山地地貌为主，该地震的发生造成了局部边坡崩塌、滚石和冻土开裂以及兰新高铁大桥桥面受损等次生灾害（颉满斌，2022）。考虑到研究区域仍有发生强震的危险（许英才等，2022），今后的防震减灾工作中有必要加强地表起伏剧烈区域特别是山脊地貌的震害防御工作。

北京大学张勇教授为本文提供了震源破裂过程的初步结果，中国科学技术大学张海江教授为本文提供了三维速度结构，中国地震局工程力学研究所马强研究员为本文提供了仪器观测的地震烈度分布图，中国地震局地震预测研究所徐岳仁研究员与作者就青藏高原东北缘的地质概况进行了讨论，审稿专家为本文提出了宝贵意见，作者在此一并表示感谢。