2022年1月8日青海门源6.9级强地震动模拟研究

师涵博 张元生

摘要：为满足数字化时代应急处置的新要求，基于随机有限断层法的发展，利用青海门源地区地下三维速度结构的模型和vS30数据，通过逐步迭代射线追踪法、格林函数位移解析相位谱和有限断层法的联合计算，以青海门源6.9级地震为例，获得具有地表土层放大效应的強地震动模拟数据，进而绘制出研究区内PGA和烈度分布模拟图；并与实际台站记录的PGA和实地调查烈度结果相比，其烈度区划范围基本一致，同时也验证此联合计算方法可用于未来地震灾害的快速评估，并为灾后应急救援提供参考。

关键词：门源地震; 强地震动模拟; 有限断层法; 逐步迭代射线追踪法

中图分类号： P315.3+3      文献标志码：A   文章编号： 1000-0844（2024）03-0680-12

DOI：10.20000/j.1000-0844.20230920001

Strong ground motion simulation of the Qinghai M6.9 earthquake on January 8， 2022

SHI Hanbo， ZHANG Yuansheng

（Lanzhou Institute of Seismology， CEA， Lanzhou 730000， Gansu， China）

Abstract： To address the modern requirements of emergency response efforts in the digital age， we developed an approach based on the stochastic finite fault method. Specifically， utilizing an underground three-dimensional velocity structure model and vS30 data from the Menyuan region of the Qinghai Province， we applied the step-by-step iterative ray-tracing method， phase spectrum of Green's function displacement analytical solution， and finite fault method to record strong ground motion simulation data incorporating ground surface amplification effects. Furthermore， considering the Menyuan M6.9 earthquake that struck the Qinghai Province on January 8， 2022， we created simulation maps depicting the peak ground acceleration （PGA） and intensity distributions within the study area. Subsequently， comparing the simulated results with actual PGA records from monitoring stations and field survey intensity data revealed consistent intensity zonation scopes. Moreover， these results validated the utility of the proposed method for rapid assessments of future earthquake disasters， offering valuable insights for postdisaster emergency rescue efforts.

Keywords：Menyuan earthquake; strong ground motion simulation; finite fault method; step-by-step iterative ray-tracing method

0 引言

根据地震台网中心测定，北京时间2022年1月8日1时45分27秒，青海省海北藏族自治州门源回族自治县（37.77°N，101.26°E）发生MS6.9地震（下文简称为门源地震），震源深度10 km，这是门源地区自1986年8月26日MS6.4地震和2016年1月21日MS6.4地震之后发生的又一次强震活动［1］。此次门源地震的发震构造为带有少量逆断分量的左旋走滑断层，位于青藏高原东北缘冷龙岭断裂、托莱山断裂和肃南—祁连断裂（俄堡段）的阶区部位，破裂带呈NWW—SEE走向［2］，区域构造较为复杂，如图1所示。其中，冷龙岭断裂长约120 km，具有明显的线性特征，属于全新世活跃的左旋走滑断裂；而托莱山断裂则位于河西走廊南部边缘，走向NNW，全长280 km，是一条晚更新世期间活动的左旋走滑断裂［3］。此外，祁连山断裂带是青藏高原东北部最为活跃的断裂带之一，尤其是海原断裂带，在历史上曾发生过多次MW7～8的破坏性地震。在过去约35年间，门源地区已发生了3次MW>6的地震，而且这些地震的主要发震断层均位于海原断裂带中冷龙岭断裂西段［4］。频繁的地震活动对抗震救灾的时效性提出了更高的要求，因此快速、准确地模拟地震所产生的强地面运动及其地表响应具有重要意义。本文拟以门源MS6.9地震为例，进行强地面运动的模拟及烈度评估，并在此基础上展开相关的讨论和分析。

在进行地震波场模拟时，射线追踪是一个绕不开的话题，相较于波动方程法可能产生的频散以及积分方程法计算的复杂性，其原理简单清晰，无需近似等价模型的建立与调整，对波场衰减的模拟更具有简便性和合理性。本研究首先基于张元生［5］、高尔根等［6］提出的逐步迭代射線追踪法进行射线路径的计算，该方法除了追踪路径结果符合射线追踪要求外，还具有计算速度快、精度高且可以克服射线路径非唯一性的特点;接着将此方法与Xin等［7］提供的中国大陆岩石圈三维速度结构模型相结合，同时运用孙晓丹［8］基于随机有限断层法改进的格林函数位移解析法和王卫民（私人通信）利用远震波形数据反演的断层破裂模型，在频率域考虑震源、传播路径及场地效应等影响，计算出相应的傅里叶幅值谱和格林函数位移解析相位谱，并将两者结合，得到相应的地震时程数据，进而得到研究区域的地震动峰值加速度（Seismic Peak Ground Acceleration，PGA）及仪器烈度（Instrument Intensity，II）分布图。

此方法的合理性和稳定性已在文献［8］中有详细论述。本文在此基础上进行了射线路径计算的优化以及参数完备性的补充与讨论，使此方法的原理更为清晰，方便后续的研究和改进，可为准确模拟地震灾害分布以及防震减灾工作提供一定的理论支持与科学依据。

1 理论与方法

1.1 三维速度结构下的射线追踪法

随着地震学以及各种观测手段的不断发展，人们对地下介质结构及其性质有了更加清晰的认识，使得准确、快速地模拟强地震动成为可能。本文旨在模拟近场范围内直达波产生的地震动分布，因此相较于效率不高的波动方程法［9］，采用逐步迭代射线追踪法［5-6］为基础，在三维不均匀块状速度结构中进行射线追踪的难点主要有：当介质横、纵向都不均匀时，对于射线上某一点的调整较为复杂;除了两种轴交点及两种网格交点的基本情况外，还需要考虑各种特殊情况，如算得的解是否可用、射线垂直或水平入射的参数确定，以及寻找最优解等问题。下文给出相关原理的简单推导。

在设立的三维速度结构的坐标系O-XYZ中，O、A、B、C分别表示三维速度结构与地表水平面的四个交点，Z轴垂直向下。S和R分别为随机选取的震源地表投影点（震中）和接收点，如图2所示。

为了简化计算，上述点的坐标都以经、纬度来表示。还需明确X、Y、Z轴分别被速度所划分的刻度值，均以km为单位，从0开始，到对应的速度结构边界结束。在开始计算之前，需要确保所有数据都在同一单位下，对于经、纬度向距离的转换，本文采用Haversine公式，相较于Great-circle distance公式使用大量余弦函数会导致较大的舍入误差，而其采用正弦函数，即使距离很小，也能保证足够的有效数字：

D=2R·arctan sin2（dlat）+sin2（dlon）·cos（lat1）·cos（lat2） 1-sin2dlat2-sin2dlon2·cos（lat1）·cos（lat2） （1）

式中：D表示S（lon1，lat1）和R（lon2，lat2）之间的距离;R为地球平均半径，单位都为km;dlon和dlat分别代表两点间经、纬度之差的绝对值。

在速度结构确定后，根据Koch［10］、张元生等［11］的研究结果可知：进行三维射线追踪时，如果忽略横向折射大约5%的横向速度扰动，走时结果误差在0.1 s左右。本研究选取的速度结构横向不均匀性较小，因此这一误差是完全可以接受的;据此，在三维速度结构中，本文参考二维速度随机分布逐步迭代射线追踪法［6］，自动截取震源点与接收点所在的二维速度结构平面进行不均匀介质下射线路径的计算，可大大减少计算量和计算时间。此方法的原理基于Snell折射定律，如式（2）所示：

sinθ1v1=sinθ2v2 （2）

当射线交点位于横向轴时，路径的偏移变化如图3所示。

从震源点出发，顺序选取路径上的三个点P1、P2、P3进行迭代，P2为中间点，经过一次迭代后，射线路径由图3中的虚线校正为实线，即中间点由P2移动至P21，计算出偏移量Δx，即可得到校正后的射线路径，如式（3）所示：

v1（x3-x2-Δx）（x3-x2-Δx）2+（z3-z2）2=   v2（x2-x1+Δx）（x2-x1+Δx）2+（z2-z1）2（3）

同理，当交点位于纵向轴时，只需将x、z的坐标调换即可，具体如下：

v1（z2-z3+Δz）（z2-z3+Δz）2+（x3-x2）2=   v2（z1-z2-Δz）（z1-z2-Δz）2+（x2-x1）2（4）

当交点位于两轴交点时，需区分上、下偏移两种情况，分别计算出从点P2到P21和P22所对应的偏移量，如图4所示。

结合Snell折射定律，当式（5）、（6）联立无实数解时，求解式（7）、（8）即可。其公式如下：

v1（Δz）Δx2+Δz2=v2（z1-z2-Δz）（z1-z2-Δz）2+（x2-x1）2（5）

v3（Δx）Δx2+Δz2=v2（x3-x2-Δx）（x3-x2-Δx）2+（z2-z3）2 （6）

v1（Δx）Δx2+Δz2=v21（x2-x1-Δx）（x2-x1-Δx）2+（z1-z2）2（7）

v3（Δz）Δx2+Δz2=v21（z2-z3-Δz）（z2-z3-Δz）2+（x3-x2）2 （8）

在迭代计算过程中，首先计算整条射线路径在实数域内的偏移量Δx和Δz;再对其检查是否存在新的或者可以合并的交点;完成优化后，进行下一轮的迭代。在本研究中，我们规定当整条射线路径迭代15次或者最大偏移量低于5 m 时，认为射线路径已达到要求，不再执行计算。

1.2 基于解析格林函数相位的有限断层法

对射线路径优化后，强地震动模拟的准确性可有一定提升。除此之外，强地震动场模拟的另一重要内容是准确计算傅里叶幅值谱与相位谱。在以往的研究中，很多使用随机合成法得到的平均相位谱，使得各点地震动时程间几乎不存在相关性，地震动场的空间特征无法表达。因此，本研究参考孙晓丹［8］提出的基于解析格林函数相位的有限断层法来进行计算，其主要用格林函数位移解析相位谱代替了之前的平均白噪声随机相位谱，并充分论证了该方法的可行性与一致性。此方法还具有计算时间短、结果稳定等优点。具体实现过程如下：

1.2.1 傅里叶幅值谱

傅里叶幅值谱的计算需在频率域将震源谱S（M0，f）、传播路径衰减P（R，f）、场地效应G（f）、地震动类型因子I（f）=（2πf）n（n取值0、1、2分别表示位移、速度和加速度）进行连乘，获取最终的幅值谱FA（M0，f，R），如式（9）所示：

FA（M0，f，R）=S（M0，f）·P（R，f）·G（f）·I（f）（9）

其中，震源模型采用动拐角频率模型［12］，Sij（M0，f）为子源震源谱，可表示为：

Sij（M0，f）=CMOijHij1+ffcij2 （10）

式中：C为表达辐射方向性差异的常数;f表示频率;MOij、fcij和Hij分别表示第（i，j˙）个子源的地震矩、动拐角频率以及高频标度因子。其中，C的表达式可写为：

C=RθφFV4πρc3 （11）

式中：F为自由表面放大，可由视入射角计算得到;V为地震波水平分量含剪切波能量的比例，通常取为0.707;ρ表示震源处的介质密度;c为P波或S波的速度;Rθφ为点源辐射模式，随方位角呈对称性变化［13］。参考Onishi等［14］的研究，其可由滑动角及断层的倾角确定，P、SV及SH波对应的表达式如下：

RP=415 （12）

RSV=12sin2λ1415+13sin2（2δ）+cos2λ415+23cos2δ（13）

RSH=1223cos2λ（1+sin2δ）+13sin2λ［1+cos2（2δ）］（14）

式中：λ代表子断层滑动角;δ则代表子断层的倾角。

为使每个子断层只触发一次地震并且遵循地震矩守恒原则，子源地震矩应为各个子断层滑动量的加权平均值，如式（15）所示：

M0ij=M0DijD （15）

式中：M0为总地震矩［15］，并M0=10（1.5MW+9.1），单位为N·m;Dij、D分别表示第（i，j）个子源的滑动量以及所有子源滑动量的总和。

关于式（10）中fcij的计算方法较多，在此采用孙晓丹［8］改进的公式：

fcij=f01+M0ijM0ave13 （16）

式中：M0ave表示平均地震矩，M0ave=M0/N，N为子断层总数;f0为根据总地震矩估算的拐角频率（定拐角频率），如式（17）所示：

f0=4.9×106βΔσM013 （17）

同理，断层面上第（i，j）个子源所对应的定拐角频率可表示为：

f0ij=4.9×106βΔσM0ij13 （18）

式（17）、（18）中：Δσ表示应力降，单位为0.1 MPa;β表示震源区介质的剪切波速。

为了减小高频地震动的低估，需引入标度因子Hij进行补偿［8］，如式（19）所示：

Hij=∫f1+ff0ij22df∫f1+ffcij22df （19）

此外，地震波在地壳介质中的传播非常复杂，传播路径衰减的影响Pij（Rij，f）通常考虑两个方面：一是几何扩散GS（Rij）;另一个是非弹性衰减D（Rij，f）。其中的关系可以表示为Pij（Rij，f）=GS（Rij）·D（Rij，f），即几何扩散仅与距离有关，但其有多种表达形式。本文采用Atkinson等［16］研究中给出的三段式几何扩散曲线：

GS（Rij）=R-1ij，RijR02（20）

式中：R01=1.5H，R02=2.5H。H代表研究区内的平均地壳厚度。非弹性衰减D（Rij，f）表示由于地球介质对波能的吸收及耗散产生的非弹性衰减，可以写为：

D（Rij，f）=exp-πfRijQ（f）v （21）

式中：v表示P波或S波的波速;Q（f）为P波或S波的品质因子：Q（f）=Q0fn，Q0（频率为1 Hz时的Q值） 的大小与介质均匀程度相关联。根据周连庆等［17］的统计结果，中国大陆的Q值范围在116～585，n取值范围在0.3～1.0。对于泊松介质，P波的品质因子Qα和S波的品质因子Qβ关系［18］可表示为：Qα=9Qβ/4，将对应的v和Q（f）代入式（21）即可得到P波和S波对应的非弹性衰减。

式（9）中，场地效应项G（f）可表示为近地表幅值放大因子A（f）与高频衰减项K（f）的乘积：

G（f）=A（f）·K（f） （22）

式中：K（f）=exp（-πfk0），f為频率，k0为特定场地与路径无关的高频衰减系数。本研究采用稂子平等［19］拟合的高频衰减系数和地表以下30 m深度内的平均剪切波速度（The Time-Averaged Shear-wave Velocity to 30 m Depth，vS30）之间的关系表达式：k0=-0.042 8×lg（vS30）+0.153 3，由此加入地表土层的影响;A（f）的计算过于复杂，为了简便，本文采用Boore等［20］计算得出的北美基岩放大系数A′（f），得到最终的G（f）=A′（f）K（f）。

应用上述公式，即可求出不同子源和接收点对应的傅里叶幅值谱。为了与实际记录有更好的对应，还需将傅里叶幅值谱按视入射角分解为北南（NS）向、东西（EW）向以及垂直（UD）向进行计算，所需P、SV波的真入射角可由本文1.1节中的射线路径坐标点求出，根据万永革［18］所做的相关研究，真入射角和视入射角之间存在着如下联系：

siniP=αβsinlP2， （P波） （23）

tanlS=2cotiPcot2iS， （SV波） （24）

式中：lP、lS分别表示P、SV波的视入射角;iP、iS分别为入射P、SV波的真入射角，α和β分别代表P波和S波的波速，SH波在自由表面的位移为入射SH波的两倍［18］。结合接收点的方位角（AZ）以及自由表面放大F，P、SV和SH波幅值谱在北东下坐标中的投影可以表示为：

PNS=sin （lP）cos（AZ） （25）

PEW=sin （lP）sin（AZ） （26）

PUD=-cos （lP） （27）

SVNS=cos （lS）cos （AZ+π） （28）

SVEW=cos （lS）sin（AZ+π） （29）

SVUD=-sin （lS） （30）

SHNS=2cosAZ+π2 （31）

SHEW=2sinAZ+π2 （32）

利用式（25）～（32）可以得到不同類型的波在三个方向上的傅里叶幅值谱，将同方向的幅值谱叠加并与相应的相位谱结合，即可得到该方向上的地震波时程数据，下文将围绕相位谱进行相关讨论。

1.2.2 相位谱

根据孙晓丹［8］的研究成果，可使用无限均匀介质中的格林函数位移解析相位谱进行计算，结合位移表示定理，点源引起空间中任意一点的弹性位移可表示为：

Un（x，t）=Mpq（ξ，τ）·Gnp，q（x，t;ξ，τ）（33）

式中：Gnp，q（x，t;ξ，τ）为Green函数，表示在x位置t时刻观测到的一对τ时刻在p取向上方向相反、在q取向上分开且作用在ξ上的力，在n方向上的位移响应;Mpq（ξ，τ）代表地震矩张量，根据角动量守恒，可以用9对（p=1，2，3;q=1，2，3）力偶组合来表示;Un（x，t）为引起的地表位移响应，表示在距离震源x位置t时刻由震源引起n方向的位移量。因此，式（33）可进一步表示为：

Mpq·Gnp，q =15γnγpγq-3γnδpq-3γpδnq-3γqδnp4πρ·1r4∫rβrατMpq（t-τ）dτ+ 6γnγpγq-γnδpq-γpδnq-γqδnp4πρα2·1r2Mpqt-rα- 6γnγpγq-γnδpq-γpδnq-2γqδnp4πρβ2·1r2Mpqt-rβ+ γηγpγq4πρα3·1rpqt-rα-γnγp-δnp4πρβ3·1rrqpqt-rβ（34）

式中：γi表示从震源点指向观测点的矢量在i方向的方向余弦;ρ为密度;r表示射线路径长度;τ为上升时间;α、β分别代表P、S波波速;δij=1，i=j0，i≠j。

式（34）可以分为5项：第一项为近场项，与滑动时间函数的积分有关;第二、三项为P、S波中场项，与滑动时间函数直接相关;第四、五项为P、S波远场项，与滑动速度函数有关。本研究选用钟形函数作为滑动速度函数：

（t）=0，t<01τ-1τcos2πτt，0≤t≤τ0，t>τ（35）

则对应的滑动时间函数为：

S（t）=0，t<0tτ-12πsin2πτt，0≤t≤τ1，t>τ（36）

滑动时间函数的积分可表示为：

∫rβrατMpq（t-τ）dτ=MpqHt，rβ-Ht，rα（37）

函数Ht，rc的表达式需分两种情况讨论：

当rc<τ时：

Ht，rc=0，td<0t3d6τ-τ4πtd+τ28π2sin2πτtd，0≤td

当rc≥τ时：

Ht，rc=0，td<0t3d6τ-τ4πtd+τ28π2sin2πτtd，0≤td<τt3d6τ-（td-τ）36τ+τ24π2，τ≤td

式中：td为破裂从起始点至接收点的时间延迟，td=t-tr，tr为断层子源破裂时间。本文以地震学北东下（XYZ）坐标系为基准，求取子源矩张量Mpq，其可由子源地震矩MOij、走向角φ、倾角δ以及滑动角λ来确定［18］。

综上即可得到含有高频成分的地表地震动时程Un（x，t），进行傅里叶变换后，便可提取出不同子源与接收点之间所对应的相位谱。

2 参数的选取

本文选用开源且兼容性较强的Python［21］编程语言开发代码，借助CPU并行计算实现强地震动的快速模拟。在计算之前，需要输入的参数包括：三维速度结构与地面4个交点（本文1.1节中的O、A、B、C点）的经度和纬度、矩震级、平均应力降（MPa）、需要计算波形的时间范围（起始时间，终止时间，时间间隔）、研究区内S波品质因子Q（f）中的Q0和n、震源处的介质平均密度（g/cm3），以及研究区地壳厚度（km）等参数。

本文选用门源地区的平均地壳厚度为55 km;根据Xu等［22］使用强震数据反演此次门源地震的相关参数获取震源处的平均介质密度为2.72 g/cm3;矩震级MW6.6以及研究区各点的vS30数据均从美国地质调查局（United States Geological Survey，USGS）网站获取;赵燕杰［23］基于台网地震记录，利用Atkinson方法反演得到门源地区S波非弹性衰减为Q（f）=137.0f0.834 2;郑瑞等［24］基于反演模型估算的地震平均应力降约为5.9 MPa;研究区域地下三维速度结构来源于Xin等［7］提供精度为0.5°×0.5°的中国大陆岩石圈速度结构模型USTClitho1.0，进行插值后如图5所示。

图5中不同颜色代表的P波速度值与右下角的色标尺相对应，大致范围在4～6.5 km/s，深度范围在地下0～20 km。可以看出，随着深度的增加，此速度结构的P波速度值逐渐增大，但横向不均匀性相对较小。在本研究中，视地下介质为泊松介质，因此P波波速α与S波波速β的关系应为α=3β。此外，本文断层模型采用王卫民（私人通信）所做的门源地震断层破裂模型，断层面上的静态滑动分布如图6所示。

从图6中可以看出，此模型共有19×9=171个子断层，分辨率为2 km×2 km，深度范围在地下1～17 km，可以适用于图5所示的速度结构范围;此断层结构较为复杂，断层面上不同颜色代表的滑动量大小与左下角的色标尺相对应，其中最大滑移量为470 cm，靠近地表位置的滑移量分布集中且相对较大;子断层倾角均为87°，但走向角的变化较为显著，按断层走向方向可分成5组，依次为：76°、106°、104°、121°和131°，计算地震动时所需的子断层参数和部分数据如表1所列。

将表1中的参数与本文1.2.2节的内容相结合，可计算得到无限均匀空间中格林函数的位移解析解;然后对其采用离散数据求导，可得到相应的加速度时程;运用快速傅里叶变换，分别提取出位移和加速度的相位谱;与本文1.2.1节对应的傅里叶位移幅值谱和加速度幅值谱相结合，经过傅里叶逆变换即可得到位移与加速度的时程数据。对于同一个接收点，为了避免出现频谱泄露，将每个断层子源对其合成的地震动时程在时域进行叠加，最终得到高频段的地震动时程波形数据。需注意的是，此方法合成的地震记录与实际地震记录之间有差异是正常的，因为使用平均随机相位谱或格林函数相位谱法合成结果代表的是某一地震引起地震动的平均水平。从统计意义上来说，一次地震发生时的地面运动仅是一个样本的作用结果，并不能代表多次地震引起的地震动的平均水平［8］。

3 结果分析

根据上述的理论方法和断层破裂模型可计算得到门源地区的PGA的分布。对于II的求取，本文采用王海云等提出的经验关系式［25］：

lgPGV=-0.397 7+0.691 4lgPGA，（R=0.863 0）（40）

II=2.259 5+1.344 9lgPGA+1.800 8lgPGV，（R=0.991 7）（41）

式中：PGV代表地震动峰值速度;R为拟合优度指数，其值越接近1，相关度越高。根据同震效应的影响范围，本研究选取的计算区域为（36.8°N～38.8°N，100°E～102.5°E）、单元网格尺寸为 0.1°×0.1°，共计26×21=546个地表接收点，采用克里金插值中的立方插值法，绘制出PGA和II的分布结果如图7所示。

图7（a）为数值模拟得到的研究区PGA等值线图，右侧色标尺的不同颜色对应左图中不同的地震动峰值加速度值。从中可以看出：PGA高于380 cm/s2的范围较小，且靠近发震断层的区域有较为明显的断层走向分布趋势;随着震源距的增加，在PGA大于40 cm/s2的区域，地震动逐渐衰减且趋于断层四周大致呈橢圆状均匀分布;整体来看，研究区PGA在震中NE向有局部小尖端、ES向有较为明显的区域突出。结合上文的论述可以推断，这种异常分布可能是由断层滑动分布和区域速度结构导致的射线路径变化共同引起的。

图7（b）中不同II值的分布与右侧色标尺的颜色相对应。与图7（a）相比，在靠近发震断层位置的区域，两者分布趋势较为相似;随着震中距的增加，烈度分布更趋于椭圆状，除震中ES向外，基本没有明显的异常凸起或凹陷。模拟结果显示，此次门源MS6.9地震的最大烈度为Ⅸ度，且越靠近震中，烈度梯度的变化越大。

在得到数值模拟结果后，还需要将模拟结果与实际观测结果进行对比，以检验方法的可靠性并分析差异产生的原因。因此，本文采用此次门源地震实际台站记录的PGA，以及实地调查烈度等震线与图7的模拟结果进行比对。

将此次门源地震在甘肃省境内的台站PGA记录，结合部分模拟PGA数据（虚拟台站）使用GMT6［26］内置的surface插值后，呈现出的分布趋势如图8所示。

从图8中可以看出，大部分实际台站位于震中的NE方向，其中H、L分别表示PGA的区域最高值和最低值。由于数据有限，无实际合站的区域使用PGA模拟值进行了填补，因此填补区并不具有对比分析的意义，我们将重点对图7（a）和图8右上角的实际合站密集区进行研究。

通过将图7（a）和图8震中NE方向等值线所示数值比较可得：图8中实际台站位置处的PGA值基本与图7（a）一致，具体数值如表2所列。图7（a）NE向的PGA凸起以及右侧边缘处的异常在图8中对应位置也有所体现，PGA分布趋势和数值基本都有很好的对应。初步分析可得，本文的模拟结果可以与实际记录有相似的趋势分布，且地下速度结构的影响不容忽视。

通过表2数据分析可得，PGA台站值与模拟值基本一致：两者比值最大为1.22，最小为0.75。越靠近震中，模拟的PGA值相对偏小，但从具体比值来看，差异都在可接受范围内;除此之外，根据卢育霞等［27］、史大成等［28］对PGA放大效应与vS30的关系研究可知：地震动峰值加速度的放大程度主要与近地表vS30相关，vS30越大地震动峰值加速度的放大越小。对比分析表2中的后两列数据可知：vS30的变化对台站值与模拟值的比值影响并不大，这说明本文所用方法在不同的地表土层放大作用下，可以得到相对准确的地震动峰值加速度。将上述结论与图7（a）和图8所示的PGA分布结合分析可得，本文模拟得到的研究区PGA结果具备一定的合理性和可信度。

本研究将图7（b）与2022年1月11日应急管理部发布的烈度图（下文简称为实际烈度图）进行了对比分析（图9）。从中可以看出，模拟烈度与实地调查的烈度范围基本一致。但在Ⅸ度和Ⅷ度区，两者在震中NE向略有突出，Ⅶ度和Ⅵ度区相对较为吻合。在震中NW方位，模拟值与实地调查值的差别随着烈度的减小不断增大。

为了进一步分析实际烈度与模拟烈度差异形成的原因，本文根据来源于USGS的数据绘制了研究区内的vS30分布，如图10所示。

从图10可以看出，研究区内vS30的变化较为剧烈，但高低值分布的区域性较为明显。高值区基本呈NW—SE方向分布，与图7的趋势分布较为一致。根据前人的研究结论［27-28］，vS30值越大的地方，烈度模拟值与实际值的差异越小，且随着震中距的增加，这种差异越发明显;图9中Ⅸ度和Ⅷ度在震中NE向的突出，模拟烈度与实地调查烈度在NW向有差异的区域，基本都位于图10中的高值区，结合上述研究结论，两者间不应有如此大的差异，所以应该考虑是否受到了地形因素的影响。

结合图1来看，震中NW向的构造相对更为复杂，而山谷、河谷地形对地震烈度也有一定的放大效应［29］，因此除了vS30之外，本文烈度结果的差异应该还受到了地形和地区房屋结构的共同影响。另外，本文仅计算了直达P、S波在地表的地震动响应，待后续代码算法改进及相关数据完善后，将结合多个震例对此进行深入讨论研究。

将图7（b）和图9结合来看，Ⅷ度以上（PGA≥380 cm/s2）的区域主要沿NWW—SEE向延伸，并与实际烈度范围大致相同，也与断层的主要走向方向较为一致。余下的区域大致呈椭圆状分布，除NW向外，在震源的其余方位均与实地调查烈度区域有很好的对应。在考虑了vS30的影响后，这种差异大概率是由于地形的影响而导致的，除此之外，在实地烈度调查中，人为因素的影响也不容忽视。结合地形、vS30以及PGA和烈度的实际观测数据可知，本文无论是PGA还是烈度的模拟结果，均与实际记录和调查有着较好的对应。这说明本研究所采用的模拟方法可以适用于地震发生后快速估计烈度区划范围，为震后灾害评估和应急救援提供参考，同时也表明地形构造的分布特征对地震PGA和烈度的分布会有较大影响。

4 讨论与结论

4.1 讨论

综上所述，本文通过将逐步迭代射线追踪法与格林函数位移解析相位谱的有限断层法相结合，以2022年1月8日青海门源MS6.9地震为例，从基本的公式和理论出发，详细论述了模拟高频地震动时程的过程以及相关参数的选取等问题，结合青海门源地区地下三维速度结构［7］与远震数据反演得到的断层破裂模型，计算了P、S直达波在该区域地表引起的地震动，并将得到的PGA和II分布与实际数据进行了对比，结果如图7～9所示。

从PGA对比结果来看，模拟结果与实际记录有相似的趋势分布，数值上的差异也在合理范围之内。从烈度对比结果来看，两者在高烈度区较为吻合，主要沿NWW—SEE向延伸;结合vS30的分布来看，低烈度区由于地形的影响在NW向有较大的差异。对于强地震动模拟来说，其结果由地形、断层模型、速度结构以及震源模型共同决定，参数多而复杂，采用不同的模型和参数可能会有不同的结果。随着后续代码算法与房屋及地形等数据的优化和充实，可对此进一步研究以提高此方法的适用性，使模拟结果与实际情况更加吻合。

4.2 结论

本文以门源地震为例，进行了强地震动数值模拟的理论应用研究，主要获得以下结论：

（1） 基于前人的研究成果，使用解析格林函数位移解析相位谱代替随机有限断层法中的随机相位谱，并与逐步迭代射線追踪法相结合，对前人模拟强地震动的方法进行了优化。同时对相关公式进行了详细论述并成功模拟得到研究区PGA和II的分布，证明了所用方法的可行性。

（2） 通过将模拟的PGA与实际记录对比分析，结果表明：本文所用的方法可以很好地模拟出与实际类似的趋势分布;结合表2中的数据，两者间的数值差异也在可接受范围之内，充分证明了此方法获取地震动模拟结果的合理性。

（3） 将图7（b）与图9对比分析可得：模拟烈度与实地调查烈度结果基本一致，存在差异的区域分别位于震中NE和NW向;结合图1和图10分析，推断地形对此有一定的影响。尽管如此，相关结果仍可表明本文所用的模拟方法能够较好地应用于震后烈度范围的快速估算。

致谢：本文断层滑动模型数据来源于王卫民研究员，图件采用Python［21］和GMT6软件［26］绘制而成，审稿专家为本文的完善提供了建设性修改意见，使本文描述更为完善，特此致谢！

参考文献（References）

［1］ 左可桢，罗钧，赵翠萍，等.青海门源地区地震活动时空分布特征和孕震环境研究［J］.地球物理学报，2023，66（4）：1460-1480.

ZUO Kezhen，LUO Jun，ZHAO Cuiping，et al.Spatiotemporal distribution characteristics of seismicity and seismogenic environment in the Menyuan area，Qinghai Province［J］.Chinese Journal of Geophysics，2023，66（4）：1460-1480.

［2］ 万永革，黄少华，王福昌，等.2022年门源地震序列揭示的断层几何形状及滑动特性［J］.地球物理学报，2023，66（7）：2796-2810.

WAN Yongge，HUANG Shaohua，WANG Fuchang，et al.Fault geometry and slip characteristics revealed by the 2022 Menyuan earthquake sequence［J］.Chinese Journal of Geophysics，2023，66（7）：2796-2810.

［3］ 王辽，谢虹，袁道阳，等.结合野外考察的2022年门源MS6.9地震地表破裂带的高分七号影像特征［J］.地震地质，2023，45（2）：401-421.

WANG Liao，XIE Hong，YUAN Daoyang，et al.The surface rupture characteristics based on the GF-7 images interpretation and the field investigation of the 2022 Menyuan MS6.9 earthquake［J］.Seismology and Geology，2023，45（2）：401-421.

［4］ WANG J Y，DING L，HE J K，et al.Research of seismogenic structures of the 2016 and 2022 Menyuan earthquakes，in the northeastern Tibetan Plateau［J］.Remote Sensing，2023，15（3）：742.

［5］ 张元生.三维分块模型速度随机分布逐次迭代射线追踪方法［C］// 中国地球物理学会第十四届学术年会论文集.北京：中国地球物理学会，1998：318.ZHANG Yuansheng.Three-dimensional chunking model velocity random distribution by iterative ray tracing method［C］//Proceedings of the Fourteenth Annual Academic Conference of the Chinese Geophysical Society.Beijing：Chinese Geophysical Society，1998：318.

［6］ 高尔根，徐果明.二维速度随机分布逐步迭代射线追踪方法［J］.地球物理学报，1996，39（增刊1）：302-308.

GAO Ergen，XU Guoming.A new kind of step by step iterative ray-tracing method［J］.Chinese Journal of Geophysics，1996，39（Suppl01）：302-308.

［7］ XIN H L，ZHANG H J，KANG M，et al.High-resolution lithospheric velocity structure of continental China by double-difference seismic travel-time tomography［J］.Seismological Research Letters，2019，90（1）：229-241.

［8］ 孫晓丹.强地震动场估计中若干问题的研究［D］.哈尔滨：哈尔滨工业大学，2010.

SUN Xiaodan.Some issues on estimation of strong ground motion field［D］.Harbin：Harbin Institute of Technology，2010.

［9］ 俞伟哲.基于层状弹性介质的全波地震波场数值模拟［D］.青岛：中国海洋大学，2010.

YU Weizhe.Full wave seismic wave field numerical simulation based on the layered elastic medium［D］.Qingdao：Ocean University of China，2010.

［10］ KOCH M.A numerical study on the determination of the 3-D structure of the lithosphere by linear and non-linear inversion of teleseismic travel times［J］.Geophysical Journal International，1985，80（1）：73-93.

［11］ 张元生，李清河，徐果明.联合利用走时与波形反演技术研究地壳三维速度结构（Ⅰ）：理论与方法［J］.西北地震学报，1998，20（2）：8-15.

ZHANG Yuansheng，LI Qinghe，XU Guoming.Combined inversion technique to study 3-D crustal velocity structure by using seismic wave travel-time and wave form （Ⅰ）：theory and method［J］.China Earthquake Engineering Journal，1998，20（2）：8-15.

［12］ MOTAZEDIAN D.Stochastic finite-fault modeling based on a dynamic corner frequency［J］.The Bulletin of the Seismological Society of America，2005，95（3）：995-1010.

［13］ 温瑞智，任叶飞，王宏伟，等.强震动记录分析与应用：芦山MS7.0地震为例［M］.北京：地震出版社，2017：158.

WEN Ruizhi，REN Yefei，WANG Hongwei，et al.Analysis and application of strong earthquake records：taking Lushan earthquake as an example［M］.Beijing：Seismological Press，2017：158.

［14］ ONISHI Y，HORIKE M.The extended stochastic simulation method for close-fault earthquake motion prediction and comments for its application to the hybrid method［J］.Journal of Structural and Construction Engineering （Transactions of AIJ），2004，69（586）：37-44.

［15］ HANKS T C，KANAMORI H.A moment magnitude scale［J］.Journal of Geophysical Research：Solid Earth，1979，84（B5）：2348-2350.

［16］ ATKINSON G M，MEREU R F.The shape of ground motion attenuation curves in southeastern Canada［J］.Bulletin of the Seismological Society of America，1992，82（5）：2014-2031.

［17］ 周连庆，赵翠萍，修济刚，等.利用天然地震研究地壳Q值的方法和进展［J］.国际地震动态，2008，38（2）：1-11.

ZHOU Lianqing，ZHAO Cuiping，XIU Jigang，et al.Methods and developments of research on crustal Q value by using earthquakes［J］.Recent Developments in World Seismology，2008，38（2）：1-11.

［18］ 万永革.地震学导论［M］.北京：科学出版社，2016.

WAN Yongge.Introduction to seismology［M］.Beijing：Science Press，2016.

［19］ 稂子平，俞瑞芳，肖亮，等.局部场地地震动高频衰减系数估计模型［J］.地震学报，2023，45（5）：919-928.

LANG Ziping，YU Ruifang，XIAO Liang，et al.An estimation model of high frequency attenuation coefficient of ground motion for local site［J］.Acta Seismologica Sinica，2023，45（5）：919-928.

［20］ BOORE D M，JOYNER W B.Site amplifications for generic rock sites［J］.Bulletin of the Seismological Society of America，1997，87（2）：327-341.

［21］ Python Language Reference［EB/OL］.Python Software Foundation，2021.[2023-09-23].Available：https：//docs.python.org/3/.

［22］ XU C Y，ZHANG Y，WANG R J，et al.Application of MEMS data to fast inversion of rupture process：tests with recordings from the IRREEW network［J］.Seismological Research Letters，2023，94（4）：1821-1835.

［23］ 赵燕杰.用Moya方法反演门源地区数字地震台站场地响应［J］.地震研究，2016，39（增刊1）：101-106，134.

ZHAO Yanjie.Site response of digital seismic stations in Menyuan area inversed by Moya method［J］.Journal of Seismological Research，2016，39（Suppl01）：101-106，134.

［24］ 郑瑞，王琪，邹蓉，等.2022年青海门源MW6.6地震InSAR同震形变场与震源特征［J］.地球物理学报，2023，66（8）：3218-3229.

ZHENG Rui，WANG Qi，ZOU Rong，et al.InSAR coseismic deformation monitoring and source characteristics of the 2022 Qinghai Menyuan MW6.6 earthquake［J］.Chinese Journal of Geophysics，2023，66（8）：3218-3229.

［25］ 王海云，李強.震后近断层震动图的快速产出研究：以2022年1月8日青海门源地震为例［J］.世界地震工程，2022，38（2）：1-9.

WANG Haiyun，LI Qiang.Study on rapid generation of near-fault shakemaps after an earthquake：a case of Menyuan earthquake on January 8，2022，Qinhai Province［J］.World Earthquake Engineering，2022，38（2）：1-9.

［26］ WESSEL P，SMITH W H F.New，improved version of generic mapping tools released［J］.Eos，Transactions American Geophysical Union，1998，79（47）：579.

［27］ 卢育霞，王良，魏来，等.利用场地表征参数研究岷漳地区地震动相对变化趋势［J］.防灾减灾工程学报，2018，38（2）：359-366.

LU Yuxia，WANG Liang，WEI Lai，et al.Study on relative variability of ground-motion in Min-Zhang Region by using representative site parameters［J］.Journal of Disaster Prevention and Mitigation Engineering，2018，38（2）：359-366.

［28］ 史大成，温瑞智，杜春清.区域性场地vS30及峰值加速度放大系数估算方法［J］.地震工程与工程振动，2012，32（4）：40-46.

SHI Dacheng，WEN Ruizhi，DU Chunqing.Study on regional site vS30 and PGA amplification factor［J］.Journal of Earthquake Engineering and Engineering Vibration，2012，32（4）：40-46.

［29］ 李平，刘红帅，薄景山，等.汶川MS8.0地震河谷地形对汉源县城高烈度异常的影响［J］.地球物理学报，2016，59（1）：174-184.

LI Ping，LIU Hongshuai，BO Jingshan，et al.Effects of river valley topography on anomalously high intensity in the Hanyuan town during the Wenchuan MS8.0 earthquake［J］.Chinese Journal of Geophysics，2016，59（1）：174-184.

（本文编辑：张向红）

基金项目：甘肃省科技重大专项“甘肃重点地区地震预测预警新方法与新技术”（21ZD4FA011）

第一作者简介：师涵博（1999-），男，硕士研究生，研究方向：固体地球物理学。E-mail：1598139042@qq.com。

通信作者：张元生（1965-），男，研究员，硕士生导师，主要从事地球内部结构及动力学环境研究。E-mail：zhangys@gsdzj.gov.cn。

師涵博，张元生.2022年1月8日青海门源6.9级强地震动模拟研究［J］.地震工程学报，2024，46（3）：680-691.DOI：10.20000/j.1000-0844.20230920001

SHI Hanbo，ZHANG Yuansheng.Strong ground motion simulation of the Qinghai M6.9 earthquake on January 8， 2022［J］.China Earthquake Engineering Journal，2024，46（3）：680-691.DOI：10.20000/j.1000-0844.20230920001