引言

2022年1月8日青海省门源县发生M_S6.9地震，导致兰新高铁运行区间大梁隧道、祁连山隧道、硫磺沟大桥等遭受破坏。迄今为止，门源地区共发生3次强震，此次地震灾害较重，引起了广泛关注。

众多学者已开展了门源地区震源机制、地质构造、波场传播等方面的研究。左可桢等（2018）基于台阵资料，采用双差层析成像法联合反演了该地区的地壳三维速度结构，并详细分析了2016年1月21日门源M_S6.4地震地质构造对地震活动性的影响；韩立波（2022）采用CAP法反演了2022年1月8日门源M_S6.9地震的震源机制，结果显示本次地震与前2次M_S6.4地震震源机制解不同，显示出该区域地质构造的复杂性；李振洪等（2022）基于2016年1月26日门源M_W5.9和2022年1月8日门源M_W6.7地震的地表形变确定震源参数，并分析了门源M_W5.9地震对M_W6.7地震的影响机制，结果显示M_W6.7地震为左旋走滑型地震，破裂带呈现NWW-SEE走向；许英才等（2022）根据地震台网资料，计算了现今应力场体系在2022年1月8日门源M_S6.9地震震源机制2个节面产生的相对剪应力和正应力，结果显示门源地区仍存在应力积累，仍可能发生强震；尹晓菲等（2022）采用有限差分法模拟2022年1月8日门源M_S6.9地震并得到烈度分布，结果显示最大烈度位于震源破裂起始点附近，且模拟得出的烈度与仪器烈度分布较为接近，值得指出的是，相较于低分辨率介质模型模拟结果（徐剑侠等，2015），尽管其模拟最高频率（0.4 Hz）得到了一定的提升，但模拟得到的烈度较实测烈度仍偏低，因此，如何提升模拟最高频率值得进一步研究；Xie等（2023）利用发震断层100 km范围内记录的加速度数据，研究了2022年1月8日青海门源M_S6.9地震中近源地震动的空间分布和衰减特征；朱音杰等（2023）和戴丹青等（2023）分别采用不同的方法反演了2022年1月8日青海M_S6.9地震的震源破裂过程，二者给出的地表破裂长度均与现场调查吻合，但二者给出的最大错（滑）动量分别为1.5 m和3.6 m，相差较大。

考虑到复杂地表起伏对地震波场传播存在重要影响，本文采用谱元法（SEM）模拟含有复杂起伏地形构造的地震波场。该方法最早由Patera（1984）等应用于流体动力学计算，具有收敛快、精度高、灵活性强等特点，近些年已广泛应用于多尺度复杂场地地震模拟及震源机制反演等。Komatitsch等（1999）将SEM拓展于地震学，用于求解地震波传播散射问题，并成功地模拟了洛杉矶盆地全过程地震动。Pilz等（2011）使用该方法模拟了智利圣地亚哥大都市区内的强地震动过程；胡元鑫等（2011）采用谱元法模拟了2008年汶川地震的地形效应，结果显示，与二维地形结果相比，三维地形对地震动的影响更为复杂；Magnoni等（2014）采用谱元法对2009年4月6日意大利中部L'Aquila 6.3级地震产生的复杂波场进行了数值模拟；Abraham等（2016）通过高性能谱元程序对1984年沿Gubbio断层发生的5.7级地震进行模拟，分析了与盆缘效应相关的地震动特征；于彦彦（2016）基于改进的高精度并行谱元方法对我国施甸盆地、四川盆地地震效应进行了详细研究，通过对比模拟结果与地震记录验证了该方法的可行性。刘中宪等（2017）基于有限断层假定的动力学震源模型，采用谱元法详细研究了近断层沉积盆地地表的地震动时程和峰值变化规律；Vijaya等（2020）采用谱元程序SPEED模拟了Kutch盆地三维地震响应；巴振宁等（2021）采用谱元法研究了动力学逆断层地震作用下的三维山体动力响应。

基于以上研究，考虑到复杂山体地形效应，本文采用高性能谱元法模拟2022年1月8日门源M_S6.9地震强地面运动，通过网格过渡使最高频率达到2.912 Hz，同时根据USGS提供的滑移分布赋予每个谱元点不同的位移值，以实现设定频率下地震波传播的精细化模拟，以期为地形复杂山体区域地震预测与防震减灾提供参考。

1.   物理模型

2022年1月8日青海门源地震过程中的断层滑动延伸至地表，据现场考察，此次地表破裂由2条断裂带组成，分别为NW向的冷龙岭断裂西段和EW向的托莱山断裂东段（Liu等，2019；Li等，2022），这2条断裂带位置与2次余震分布带吻合（图1）。

图 1  M_S6.9门源地震发震构造及历史地震活动

Figure 1.  Seismogenic structure and historical seismicity around the M_S6.9 Menyuan earthquake

下载: 全尺寸图片 幻灯片

本文研究区域如图2所示，可以看出，地形起伏变化较大，东北方向地形较为平坦，西南方向山体聚集，地形起伏较大。地表高程的数据来源于NASA9.0 m。物理模型如图3（a）所示，其平均高程为3421.86 m，最大高差为2678.47 m，整个模型尺寸为167.8 km×114.4 km×40 km。由于SWChinaCVM-2.0（水平分辨率为0.2°）暂无本文研究区域的速度结构，故本文采用左可桢等（2018）获得的该地区速度结构（水平分辨率为0.3°），该结构精度高于CRUST1.0（水平分辨率为1°）和USTClitho2.0（水平分辨率为0.5°）。由于缺乏相关介质的衰减资料，假定${Q_{\text{s}}}$、${Q_{\text{p}}}$与频率无关，本文${Q_{\text{s}}}$和${Q_{\text{p}}}$分别取${V_{\text{s}}}$和${V_{\text{p}}}$的1/10（Peyrusse等，2014）。地层介质参数如表1所示。

图 2  研究区域

Figure 2.  Study region

下载: 全尺寸图片 幻灯片

图 3  模型参数

Figure 3.  Model parameters

下载: 全尺寸图片 幻灯片

表 1  地层介质参数

Table 1.  Stratum media parameters

地层序号	层底深度/km	ρ/(kg·m−3)	Vs/(m·s−1)	Vp/(m·s−1)	${Q_{\text{s}}}$	${Q_{\text{p}}}$
1	5	2110	3640	6076	364	608
2	10	2720	3609	6049	361	605
3	15	2790	3609	6040	361	604
4	20	2790	3660	6150	366	615
5	30	2790	3683	6176	368	618
6	40	2790	3766	6199	376	620

下载: 导出CSV 

| 显示表格

断层设置是基于长安大学李振洪等（2022）提供的震源机制解（两段式断层及滑移分布）。断层沿走向方向长55 km，沿倾向方向宽20 km。在图3（b）中，断层截面处不同颜色代表不同的位错大小，通过每个子源所在地层的剪切模量、相对位错等参数计算得到相应地震矩，对每个子源赋予相应的地震矩以模拟整个断层的错动，值得注意的是，谱元点均布在2个地层之间，通过此方法来模拟跨越地层的断层。由图3（c）可知，断层有1个主要破裂过程，即断层平面上的深红色区域，位错最大约为3.5 m，对于门源地震，由于断层的大部分位于地层2之中，故本文假定破裂速度为地层2剪切波速的0.8倍，即2.8 km/s，通过子震源与初始破裂点位置计算得到破裂时间，该破裂时间具有明显的扩散性。参考USGS，震源时间函数如图3（d）所示。

2.   计算方法

本文采用开源的谱元程序SPEED ¹进行门源地震效应模拟。模型的网格划分方式采用连续网格划分，为了兼顾网格数量与计算效率，地表与基岩交界位置采用3倍过渡网格。在模型中，除地表外的5个边界面均设置吸收边界。模型共划分为3219840个单元，3343813个节点。模型网格划分如图4所示。

图 4  模型网格划分

Figure 4.  Model grid division

下载: 全尺寸图片 幻灯片

根据谱元法的计算原则，1个波长需要包含5个GLL点，将多项式近似度设置为$N = 2$，使每个单元包含${(N + 1)^3} = 27$个积分点。基于研究区域的地层剪切波速，将地表网格大小设置为400，模型的最大频率约为2.912 Hz。采用基于Linux系统的多节点并行计算技术，时间步长取0.001 s，模拟时长取80 s。

3.   地面运动模拟

3.1   波场快照

图5给出了计算区域EW方向的速度波场快照，可以看出，在t=4 s时，由断层错动激发的地震波向外扩散，从震中向外传播大致呈椭圆形分布；t=8.2 s时，地震波向外扩散，此时在断层的两侧出现反应较为剧烈的波前；t=10.2 s时，随着波前的增大，x向PGV达到41.4 cm/s，与实测的烈度分布相近；t=12.2 s时，波前进一步增大，此时的x 向PGV高达47.2 cm/s，可以看出，在山体分布密集区域的PGV响应较为剧烈，在位于东南方向的平坦地区，地震响应相对较弱；t=15.8 s时，x向PGV降至27 cm/s，随着地震波的扩散，地震动传递能量逐渐衰减，地震动在一些山体附近的持续时间长于其他位置，这可能与山体区域地震波聚集有关；t=22.8 s时，地震动传播过程基本结束。通过观察波场快照还可发现，沿断层走向方向的波场发展的较快、较大，而垂直于断层方向则较慢、较弱，具有明显的地震方向性效应。

图 5  波场快照

Figure 5.  Wavefield snapshot

下载: 全尺寸图片 幻灯片

3.2   PGV分布

图6给出了计算区域E-W、N-S和U-D分量以及3个分量矢量和的PGV分布，从3个方向的峰值云图可以看出，峰值较大的地方主要集中在断层附近，其中N-S方向的速度值较大；从合速度可以看出，沿断层走向方向的地震响应大于垂直断层走向，具有明显的地震方向性效应，这与实际地震破坏的烈度分布相符。

图 6  PGV云图

Figure 6.  PGV nephogram

下载: 全尺寸图片 幻灯片

3.3   烈度分布

根据模拟结果的PGV分布，结合GB/T 17742—2020《中国地震烈度表》（国家市场监督管理总局等，2020）给出的PGV与烈度之间的关系得到对应的地震烈度，如图7（a）所示，结果表明，垂直断层走向地震动能量衰减要明显大于沿断层走向，在模拟过程中PGV的最大值为53.1 cm/s，位于断层附近。本文模拟得到的烈度分布规律与中国地震局工程力学研究所强震动观测组 ²给出的烈度分布相近，烈度较大的部分主要沿着断层走向分布，南侧烈度大于北侧烈度，可能是由于南侧起伏的山地地貌增加了地震波反射，导致地震持续时间增长，地震破坏加剧。

图 7  模拟结果与实测烈度对比

Figure 7.  Comparison between the simulated result and the measured intensity

下载: 全尺寸图片 幻灯片

3.4   与强震记录对比

在2022年1月8日的门源M_S6.9地震中，中国数字强震动台网获得了台站的强震记录，本文计算区域内唯一的观测台站62HUC距震中约70 km（图1），该台站记录为验证本文模拟结果提供了重要参考。

图8给出了本文模拟结果与台站62HUC记录的对比。可以看出，本文模拟的加速度时程与实测记录在波形和峰值上较为符合；对于加速度反应谱，E-W方向、U-D方向与实测较为相符，N-S方向在1.6 s之前与实测较为相符，而在1.6 s之后，实测记录的长周期部分偏大，模拟结果与实测相差较大。其原因可能是本文模型采用的是恒定的破裂速度、滑动角，且地表以下土层分布平坦，没有起伏，而该区域的破裂带、地质构造以及震源的破裂过程均比较复杂。另外，此次地震仅获得1条强震记录，采样对比数量偏少。因此作者认为模拟结果与实测记录之间的差异是可接受的。

图 8  模拟结果与台站记录对比

Figure 8.  Comparison between the simulated results and the station records

下载: 全尺寸图片 幻灯片

4.   结论

本文基于青海门源地区的地表高程、介质参数以及震源函数等，采用谱元法模拟了该地区的地震传播过程，并与强震记录进行了对比，主要结论如下：

（1）PGV较大的地方主要集中在断层附近，且N-S方向较大。沿断层走向的地震响应要大于垂直断层走向，体现出明显的地震方向性效应。模拟得到的烈度分布体现出这一特点，并与实测烈度分布的规律相同。

（2）受复杂山体起伏和强地面运动方向性效应的影响，地震波在断层两侧出现反应较为剧烈的波前；山体分布密集区域的PGV响应较为剧烈，且山体附近的地震动持续时间也较长，而平坦区域的地震响应相对较弱。

（3）通过与台站记录对比，验证了本文方法的合理性与准确性。加速度时程及反应谱计算结果与实测记录在波形和峰值上较为符合，其差异主要源于实际场地的介质构造、地质条件以及震源破裂过程的复杂性。

此次地震发生在以山体为主的地区，由地震引发的诸如滑坡、泥石流、交通损毁、通讯中断等次生灾害对当地居民的影响较大，所以在今后的研究中有必要加强对地表起伏较大区域的研究，以便开展有针对性的防震措施。

致谢 本研究使用了中国地震局工程力学研究所提供的2022年1月8日青海门源M_S6.9地震近场区域的地震观测记录，以及中国地震局兰州地震研究所提供的门源地区地壳三维速度结构，审稿专家为本文提出了宝贵意见，在此一致表示感谢！