2022年3月16日日本近海7.4级地震地震动场三维有限差分模拟

李春果，王宏伟，温瑞智，强生银，任叶飞

（1.中国地震局工程力学研究所地震工程与工程振动重点实验室，黑龙江 哈尔滨 150080；2.地震灾害防治应急管理部重点实验室，黑龙江 哈尔滨 150080）

引言

据日本气象厅（JMA）测定，2022年3月16日23：36：29（UTC+9）日本本州东岸近海发生7.4级地震，震中位于141.623°E、37.697°N，震源深度为57 km，是2011年3月11日9.1级东日本大地震后日本东部地区最大震级地震。据日本防灾科学技术研究所（NIED）震度分布，地震最大震度达到6强，造成了强烈的地面运动、严重的工程结构破坏及人员伤亡。此次地震的发震区位于欧亚板块、太平洋板块、北美板块与菲律宾海板块的交界，该区域板块间相互作用复杂，太平洋板块和菲律宾板块向西俯冲形成狭长的板间地震带，欧亚板块沿边缘东边界向东移动，形成双向对冲式汇聚，拥有活跃的浅源和深源强震活动［1－4］（见图1）。地震发生后，密集的地震动观测台网（K－NET和KiK－net）观测到强烈的地震动及丰富的长周期面波，其中K－NET台网的MYG002台站获得了峰值地面加速度（PGA）最大的记录，其东西、南北和垂直 向 的PGA分 别 为800.353、609.915、414.211 cm/s2。

图1 日本浪江町7.4级地震震中及其3.0≤MJMA≤6.0余震（截至2022年3月21日）Fig.1 Epicenters of the M7.4 Namie,Japan earthquake and its aftershocks with 3.0≤MJMA≤6.0(up to March 21,2022)

日本东部海域地区强震活动十分活跃，研究人员开展了大量的地震波动数值模拟工作，Furumura等［5］、Takemura等［6］分别模拟了不同地震中日本关东盆地产生的长周期地震动，前者分析了盆地内地震波聚焦效应引起的显著的盆地放大，后者则重点讨论了震源参数对长周期地震动模拟的影响；Nakamura［7］等探讨了表面地形、海水层对地震动数值模拟的影响；Petukhin等［8］分析了包含海水层的三维速度结构模型对大阪盆地模拟地震动的影响，解释了地震面波在海洋沉积层中的产生过程，强调了地震波动模拟中考虑海水层的必要性。

文中利用日本7.4级地震的震源运动学破裂模型初步反演结果以及考虑地表地形、海水、板块边界等的三维速度结构精细化模型，基于交错网格有限差分法，模拟大尺度区域内地震波的产生和传播过程，对比强震动台站的观测记录，对数值模拟结果进行检验，重点讨论了复杂地形、海水层对地震动传播的影响。

1 数值模拟方法

有限差分法是一种求解微分方程数值解的近似方法，其主要原理是对微分方程中的微分项进行直接差分近似，从而将微分方程转化为代数方程组求解。交错网格有限差分法是运用较为广泛的方法之一［9－10］，该方法能够处理包括起伏地表和海水平面在内的复杂地形问题，采用局部差分算子增加运算速度，同时减少数值频散现象，提高计算精度［11－12］，结合高阶有限差分能够适当加大网格间距，增加运算效率［13－14］。文中采用交错网格有限差分法的开源程序OpenSWPC［15］（https：//tktmyd.github.io/OpenSWPC/）模拟地震波的产生和传播过程，在空间上选择四阶精度，时间上为二阶精度［10］。该开源程序求解了以速度－应力分量表示的三维连续介质力学运动方程，

式中：vi是第i个复合体中弹性运动的粒子速度：ρ是密度；σij是应力张量的i，j个分量；fi是 体力的第i个分量，采用广义齐纳体（Generalized Zener Body）模型的本构方程表示宽频范围的恒定Q值，将内存变量纳入本构方程以表示粘弹性介质，采用一阶Crank-Nicolson法［16］显式求解原始隐式方程。

采用有限差分法的OpenSWPC程序相较于有限元法、谱元法等程序求解形式简单，占用内存小，计算效率高，且整合了模拟过程必要的前处理和后处理程序，并包含了日本全部地区精细化三维速度结构模型。输入相应的地震震源或有限断层参数，程序将自动分配并行计算机内存生成非均匀介质模型的交错网格，进行地震波传播的数值模拟。

2 震源破裂模型及三维速度结构模型

2.1 震源破裂模型

震后多个机构给出了本次地震的震源机制解和破裂过程反演结果，文中采用了中国地震局地球物理研究所张旭等（https：//www.cea－igp.ac.cn/kydt/278892.html）利用远场体波数据初步反演的此次地震的震源运动学破裂模型。震中为37.702°N、141.587°E，震源深度63.1 km，矩震级Mw7.4，地震矩M0为1.696×1020Nm，破裂面走向和倾角分别为184°和40°，破裂面沿走向的长度和沿倾向的宽度为110 km和70 km，子断层尺寸为10 km×10 km。断层面的破裂滑动主要集中于起始破裂点附近，破裂持续时间约15 s左右，破裂速度VR=3.0 km/s，其SSW侧的破裂滑动更明显，最大滑动量约3 m，该破裂模型充分描述了滑动在破裂面上的不均匀分布，滑动主要集中于破裂起始点附近区域，该区域可看作一个凹凸体，文中也根据Somerville等［12］的方法对破裂面进行了裁剪，裁剪了边缘基本不存在破裂滑动区域（图2虚线所示）。震源时间函数为单周期的Kupper小波函数［18］，根据Ekström等［19］给出的震源上升时间与地震矩经验关系确定震源上升时间TR，根据上升时间得到震源辐射频带为1/2TR～2/TR。

图2 运动学破裂模型的震源滑动分布Fig.2 Source-slip distribution of kinematic rupture model

2.2 三维速度结构模型

地震波传播特性受介质性质的影响，特别是对地震动低频成分具有显著影响，因此精细的三维速度结构模型是准确模拟研究区域内的地震波传播过程的基础。文中模拟了400 km×500 km×100 km范围内的地震波传播过程，使用横向墨卡托投影将模型区域转换为笛卡尔坐标系统以进行球形测量［20］，中心坐标为140.5°E、37.5°N。地表和海底地形起伏数据来自于美国国家海洋和大气管理局发布的全球地表地形ETOPO1模型［21］，陆地自由表面以上的网格为空气层，波速为0，密度为1 kg/m3，为考虑海水层对地震波传播的影响，在海域将海底至海平面的范围设定为弹性海水层，海水S波波速为零，P波波速VP=1.5 km/s，密度为103kg/m3，刚度为0。三维地壳速度结构模型采用日本综合速度结构模型［22］，该模型使用23个等速层表示三维非均质地下结构以及沉积盆地的复杂基底形状，包括不规则基底地形和海底测深、低速沉积层分布、康拉德不连续面和莫霍面深度、太平洋板块边界、洋－壳与洋－幔边界等速度不连续界面等，能够用于评价日本大地震长周期强地震动。研究区域的S波速度模型如图3所示，图中可见，速度结构横向不均匀、板块交界处的显著差异等特点，由于板块交界，计算区域内西南部速度结构较为复杂，存在非均匀的壳内低速体，即分割不同板块的局部低速层。文中建立了包含复杂地表地形、海水层、板块边界等的交错网格，模型边界选择完美匹配层PML厚度为5 km（20个网格）［23－24］作为吸收边界条件。

图3 三维S波速度结构模型透视图Fig.3 Perspective view of model with three-dimensional S-wave velocity structure

为保证数值模拟精度，根据有限差分法的稳定性条件，三维模型的空间离散网格单元尺寸设定为250 m×250 m×250 m，计算区域内共划分了14.4亿个单元，时间步长设定为0.015 s，模拟了250 s内的地震波传播过程。沉积层最小剪切波速为500 m/s，根据网格尺寸、最小波速等，模拟的最高频率fmax为0.8 Hz。本研究利用了国家超级计算天津中心天河二号超算，通过MPI并行来提高计算效率，实现大尺度区域内的地震波传播的数值模拟。

3 讨论

3.1 速度波场快照

图4显示了不同时刻的速度波场快照。20 s时，地震波沿断层双向辐射，P波波峰到达断层北部的沿海地区，地震波传播相对较为规则，主要由震源破裂模型控制，破裂滑动向东侧滑冲，未表现出明显的方向性效应；35～50 s左右，强烈的P波及S波抵达福岛东部及宫城县地区，产生强烈地面运动，在震中东南部海域出现更大的速度峰值。60～80 s，地震波继续传播至关东、新潟地区，受地表地形的影响，波场呈边缘不规则的圆形，并存在局部起伏地形引起的幅值放大，体波抵达后在关东、新潟地区产生了长周期面波，东南部海底出现了破碎的散射波，关东地区低速沉积层引起波形滞后现象，且增大了地震波的幅值，倾斜入射角使得体波转换为面波，体波与面波的相互干涉在沉积边缘产生局部放大效应；100 s时，关东、新潟地区的面波波长变短、传播速度下降，同时观测到关东、新潟地区速度场的持续放大，在模型边界吸收层未发现显著的人工边界反射波。模拟波场快照较好地体现了地震波的传播特征。

图4 模拟速度波场快照Fig.4 Snapshots of the simulated velocity wavefiled

东南部海底散射波的存在表明海底高速介质层间反射和低速软流层对入射地震波的吸收折射作用［25－27］；Noguchi等［28－29］分别研究了瑞利波在海洋区域的传播及转换，瑞利波对海水层的厚度和海底地形起伏等结构非均质性较为敏感，说明了日本海沟及海水层对海底瑞利波传播的影响；关东、新潟等地区的滞后性放大现象表明深厚沉积场地对地震波传播的影响，长周期面波在盆地内低速沉积物中持续传播，造成局部速度波场的振幅放大和持续时间明显增长［30－31］。

3.2 速度波形对比

图5对比了沿海地区固定轨迹AA′（大致平行于断层迹线）上部分台站的模拟与观测记录的三分量速度时程，其中对观测记录依次进行零线校正、记录波形首位加余弦窗并补零、巴特沃斯非因果带通滤波（0.01～0.8 Hz）处理。观测记录速度波形中存在一个主要波峰，这与断层上发生了一次主破裂相对应，模拟记录速度波形也显示了一个主要波峰。震中距较大的台站的模拟与观测速度波形具有较为一致的振幅与震相，但近断层处的台站（例如MYG017、FKS004、FKS001等）模拟记录幅值略偏小。近断层地震动高频成分丰富，并受发震断层的规模、破裂特点等因素影响，表现为近断层强震动的集中性等特征［32－33］，从而使近断层观测速度波形幅值增加；震中距较远的关东及千叶地区台站速度波形出现与地形相关的较为复杂的面波持时和幅值的放大，图4中同样存在局部放大现象，表明该地区地震波在深厚沉积场地发生复杂的散射和波形转换。

图5 （续）Fig.5（Continued）

图5 观测与模拟记录速度波形对比Fig.5 Comparisons of the observed and the simulated velocity waveforms

3.3 峰值参数对比

图6对比了K－NET台网观测与模拟的峰值地面速度PGV，二者具有相似的随断层距的衰减趋势，在近断层区域模拟值普遍偏小，如东西向近断层台站（RJB<100 km）模拟记录PGV约2～8 cm/s，而观测记录PGV约6～20 cm/s，图7对比了观测和模拟记录PGV的空间分布，模拟记录PGV空间分布与观测记录较为一致，在关东、新潟地区均存在明显的地震动放大。

图6 观测与模拟记录的PGV随距离衰减Fig.6 Attenuation of PGVs versus distance for the observed and simulated ground motions

图7 观测记录和模拟记录PGV空间分布Fig.7 Spatial distributions in terms of PGVs for the observed and simulated ground motions

3.4 傅里叶幅值谱对比

图8为部分近场及远场台站速度记录的傅里叶幅值谱，其中近断层台站（MYG015、FKS004）模拟值与观测值具有一致的峰值频带（分别约为0.4、0.3 Hz）且三分量谱值较为相近，而远场台站（CHB007、CHB017）拥有较宽峰值频带（约0.15～0.7 Hz）的同时，水平向模拟谱值略小于观测结果。无论是近场或远场台站，大于0.1 Hz的范围内，速度记录与模拟结果的谱值具有较好的一致性，部分台站在小于0.1 Hz的超低频域内，模拟值与观测值存在一定的偏差，这可能与记录的有效频带有关。

图8 部分台站观测和模拟的速度记录的傅里叶幅值谱对比Fig.8 Comparisons between the Fourier velocity spectra of the observed and simulated recordings in selected stations

模拟与观测记录总体具有较好的一致性，而针对模拟记录偏小的情况，我们分析可能有以下几点原因：

（1）模拟中采用的数值网格尺寸有限，为保证满足运算的波长条件（波长内包含5～10个离散网格），运算中截止速度Vcut=500 m/s，即陆地中波速小于500 m/s的松软表层场地均按500 m/s的坚硬场地进行替代，因此并未考虑近地表更软弱的土层对地震动的影响。而在海底未进行此项替换，因而较好的模拟了由沉积物产生的丰富的面波。

（2）数值模拟选择完美匹配层吸收边界条件，假设吸收边界区域介质为完全弹性，此方法在介质模型中波速比较大的界面处极偶尔会出现不稳定性［18，30］。为保证结果的稳定，运算中忽略速度结构中的波速比过大界面处的低速层，从而使PGV的空间分布中模拟结果偏小。

（3）低估体波速度振幅可能与考虑了海水层的三维速度模型有关，海水层相较于刚性基岩，作为地下结构中的低速体能够吸收一部分地震波，从而降低陆地上的地震波振幅。Maeda［34］曾模拟了2011年日本东北311地震的三维地震动响应，并对比了模型中是否含有海水层的模拟结果，发现含海水层的结果较之不含海水层的结果偏小。

4 结论

文中利用2022年3月16日日本7.4级地震的震源运动学破裂模型初步反演结果，采用综合考虑了地形起伏、海水层、板块边界等的三维速度结构精细化模型，基于交错网格有限差分方法，模拟了地震波的传播过程，得到以下主要结论：

（1）计算区域内的地震动显示出明显的区域变化特征。模拟速度波场快照中存在局部地形引起的放大效应；同时在含有沉积层的关东、新潟等地区呈现显著的盆地效应，出现较为复杂的长周期面波及放大的地震波振幅；海底持时较长的散射波的存在则显示了海底高速介质与地表介质的层间反射和低速软流层对入射地震波的吸收折射作用，体现了计算区域内的地表地形起伏、深厚沉积层和海底结构等复杂地形对地震波传播的影响。

（2）数值模拟结果与K－NET台站的实际观测记录具有较为一致的幅值与震相。通过对比模拟记录与观测记录PGV的空间分布，模拟记录具有相似的空间分布特征但峰值速度略小，可归因于计算模型中的海水层一定程度上降低了模拟的地震波振幅，突出了海水层作为速度模型中的低速层对地震波的吸收作用。

在地震波三维数值模拟过程中，地下速度结构的精度及网格划分的疏密程度对波的传播过程影响较为显著，文中所采用的三维速度模型已相对比较精细，但模拟结果相对观测记录仍偏小。可知除速度模型以外的影响因素，包括断层分布等震源参数的设置和由于有限差分法数值运算的局限性而对三维速度模型进行的处理等，都将影响数值模拟结果的准确度，需进一步进行分析，以得到更加符合实际的模拟结果。

致谢：感谢东京大学地震研究所Takashi Furumura教授对数值模拟过程的指导。中国地震局地球物理研究所张旭提供日本本州近海M7.4地震的震源破裂过程反演结果。K－NET台站的强震动观测数据来自日本科学技术研究所（https：//www.kyoshin.bosai.go.jp/kyoshin/search/index_en.html，NIED）。