随机有限断层方法考虑了震源、传播路径和场地的全过程，使用简单高效、输入参数物理意义清晰，是地震工程中主流的强地震动模拟方法。其模拟的地震动记录在频谱、幅值以及持时上与观测记录对应良好^[1]。通过模拟已发生地震的地震动可以验证相对稳定的路径和场地参数，对于未来可能发生的地震利用现有合理输入参数结合其震源信息，使用模拟方法快速得到地震动影响场，进而为后续的地震烈度速报、震后情景构建和灾情评估等震后态势感知系统提供地震动输入。

该方法被广泛应用于世界各地的地震动模拟中^[2-5]，但这些地震多是发生在内陆或震源深度不足25 km的俯冲带浅层地壳区域，目前国内针对浅层地壳以下俯冲带地震引起的地面运动模拟研究较少。文献[6]对俯冲带板间地震的地震动进行了模拟，虽然浅层地壳地震、俯冲带板间地震和板内地震存在一定的关联，但俯冲带板内(subduction slab)事件的震源深度更深(通常超过50 km)，其地震动大于同条件俯冲带板间事件的地震动，地震动峰值衰减慢于浅层地壳地震，短周期谱相比于另外两者会更高^[7-10]。考虑到俯冲带板内和板间地震的地震动差异较大^[11]，随机有限断层法在俯冲带板内地震动模拟中的适用性还有必要进一步验证。作为地震多发国家，中国海域分布在大陆板块与海洋板块共同作用的俯冲带地区，南海北部及华南沿海历史上曾发生7级以上地震，俯冲带地震活动性较强^[12]。然而中国对俯冲带地震的相关研究起步较晚，海域地震动记录稀缺，这对地震风险评估造成严重制约，由于中国台湾地区俯冲带与日本俯冲带构造相似^[13]，日本地区俯冲带地震研究对中国俯冲带地震的地震动模拟以及地震工程研究有重要的参考价值。

北京时间2021年10月7日21时41分，日本千叶县(35.6°N，140.1°E)发生M_j 6.1地震，日本气象厅及日本总务省消防厅发布的地震震源深度为60 km。根据美国地质调查局(USGS)公开的震源机制解，此次地震断层机制为逆断层。本文利用Zhao等^[14]的日本地区地震分类方法，进一步将其确定为俯冲带板内地震，基于随机有限断层法对此次俯冲带板内地震引起的地面运动进行了模拟。

1 随机有限断层法

文献[15]在随机点源法^[16-17]和静拐角频率的有限断层方法^[18-19]基础上，提出了动力学拐角频率模型的随机有限断层方法，摆脱了子断层尺寸划分对地震动模拟结果的影响，解决了同一个子断层在地震中多次破裂的异常问题。此后为了保证随机点源法和有限断层法在远场模拟结果相匹配，文献[20]进一步做出重要改进，将子断层拐角频率的倒数作为震源上升持时，在低频段使用滤波器函数以保证低频傅里叶谱的一致性，这使得随机有限断层法进一步扩大了适用震级范围。

根据随机有限断层模型，将断层面沿走向和倾向划分为N个子断层^[21]，频域上可将第i行j列子断层在距离为R_ij的目标场点加速度傅里叶谱描述为

$ \begin{aligned} A_{i j}(f)= & \left\{\frac{C M_{0 i j} H_{i j}(2 {\rm{ \mathsf{ π} }} f)^2}{1+\left(f / f_{0 i j}\right)^2}\right\} \cdot\left\{\exp \left(-\frac{{\rm{ \mathsf{ π} }} f R_{i j}}{Q \beta}\right) G\left(R_{i j}\right)\right\} \cdot \\ & \{\exp (-{\rm{ \mathsf{ π} }} \kappa f) A(f)\} \end{aligned} $

(1)

式中：第1项为震源部分; M_0ij、f_0ij和R_0ij分别为第i行j列子断层的地震矩、拐角频率和到目标场点的震源距; $为常数，$为辐射图型系数，通常取剪切波的平均值0.55；F为自由表面放大，通常取2.0；V为2个水平分量的分配系数，一般取0.71；ρ和β分别为震源附近的介质密度和剪切波速。第i行j列子断层的拐角频率表示为

$ f_{0 i j}=4.9 \times 10^6 \beta\left(\frac{\Delta \sigma N}{N(t) M_0}\right)^{1 / 3} $

(2)

式中：Δσ为应力降；M₀为整个断层的地震矩；N(t)为t时刻已经破裂的子断层的数量之和。

为了保证高频辐射能守恒，文献[15]引入高频标定因子H_ij对低估的子断层远场辐射能进行补偿，表达式如下：

$ H_{i j}(f)=\sqrt{\frac{N \sum\left\{f^2 /\left[1+\left(f / f_0\right)^2\right]\right\}^2}{\sum\left\{f^2 /\left[1+\left(f / f_{\omega, i j}\right)^2\right]\right\}^2}} $

(3)

式(1)的第2项为路径传播部分，包括非弹性衰减$\exp \left(-\frac{{\rm{ \mathsf{ π} }} f R_{i j}}{Q \beta}\right)$和几何扩散G(R_ij)两部分，其中Q为传播介质的品质因子，R_ij为第i行j列子断层到目标台站的距离。第3项为场地影响，exp(－πκf)表示震源谱高频成分的快速衰减，A(f)为近地表场地的放大效应，包括地壳放大和局部土层放大两部分。

将式(1)表示的所有子断层在目标场点产生的加速度谱经过傅里叶逆变换，考虑断层内的破裂和地震波传播造成的时间上的迟滞，在时域内叠加即可得到该场点最终的地震动加速度时程，可表示为

$ a(t)=\sum\limits_{i=1}^{n_l} \sum\limits_{j=1}^{n_w} a_{i j}\left(t+\Delta t_{i j}\right) $

(4)

式中：n_l和n_w分别为沿断层走向和倾向划分的子断层数目之和；a_ij为第i行j列子断层的加速度时程，即式(1)中A_ij的傅里叶逆变换值；Δt_ij为破裂传播到第i行j列子断层，并从该子断层传播到目标场点的时间延迟。

2 数据处理和参数获取 2.1 地震动记录筛选及处理

本次千叶县M_j6.1俯冲带板内地震引起的地面运动模拟所用井下和井上加速度数据均来自日本防灾技术科学研究所的KiK-net强地震动台网。选取震中距≤100 km的台站，利用文献[22]基于信噪比的评估方法对记录进行筛选以确定记录质量，最终有25个台站的记录通过了阈值为5的筛选标准，可被视为高质量记录，参考台站及震中位置见图 1。

根据KiK-net台网发布的强震动记录信息和台站钻孔资料，将所用台站的震源距、井上水平向峰值加速度(A_pg)和井上台场地分类整理见表 1，其中场地分类按照文献[23]基于场地周期的方法进行划分，对于给出完整钻孔资料的台站，使用土层厚度与对应土层等效剪切波速之比计算出的特征周期(即T=4H/V)确定场地类别；未给出钻孔资料的台站，使用H/V平均反应谱比的峰值特征大致判断场地类别。

由表 1可以看出，KNGH10井上台站同时记录到了本次地震EW向和NS向地表A_pg值，分别为140.6 cm·s^-2和107.8 cm·s^-2，地表水平向A_pg>60 cm·s^-2的台站有4个，分别为CHBH10、CHBH12、CHBH16和KNGH10台站；井上台场地情况复杂，涵盖场地分类标准的全部类型。所有地震动数据经过基线校正后，使用频带范围为0.1~25 Hz的四阶巴通沃斯非因果滤波器进行带通滤波处理。

2.2 输入参数的计算和选取

本文将USGS提供的震源机制信息作为本次模拟的震源参数，矩震级M_w=5.9，震源深度取60.5 km，断层走向、倾角和滑动角取178°、33°和94°。在断层几何信息不明确的情况下，断层尺寸使用文献[24]的俯冲带板内地震震源尺寸与矩震级经验统计关系估算，将其确定为10 km×8 km。结合上述断层尺寸估算出的上缘埋深为58.5 km。应力降利用文献[25]给出的破裂面积、应力降与地震矩的经验关系估计出初始值为3 MPa，之后通过“试错法”最佳匹配拟反应谱(A_ps)残差平均值以确定本次模拟使用的应力降^[26]，增量设置为0.5 MPa，最终将应力降取为6 MPa。本文使用随机滑动分布的方式对子断层滑动权重进行分配。震源附近介质密度和S波波速由Crust 1.0全球地壳模型确定，分别为3.17 g·cm^-3和4.29 km·s^-1。

对于地震动持时，文献[4]使用90%能量持时作为地震动持时模型，这使得模拟加速度记录在持时上更加接近观测记录，本文认为70%能量持时更能体现S波主导的强地震动持续时间。因此，取地震动能量处于5%和75%之间的持续时间(即70%能量持时)作为本次模拟的地震动持时模型，以距离为自变量，回归得到二者的经验关系式

$ T=a+b \cdot R $

(5)

式中：T可被视为包含了震源持时和路径持时的地震动持时模型(以下记为70%能量模型)，a、b为模型系数，R为震源距。由于场地对于地震动持时的影响尚不明确，以上模型可能包含场地项对于持时的影响。本文计算了以上所选台站记录的能量持时，去掉其中持时异常过大和过小的台站(CHBH20和KNGN20)记录，对筛选后的记录使用最小二乘法进行回归，结果见图 2。圆点表示观测记录的能量持时，点划线表示地震动持时模型，对应井下和井上地震动持时模型函数分别为T_down=0.129 7R+0.797 95和T_up=0.090 8R+4.930 64。

路径效应部分包括由几何扩散引起的弹性衰减以及在介质中传播损耗引起的非弹性衰减。几何扩散一般可用分段连续函数来简化表示，本文针对研究区域参考了文献[2]对日本3.11地震的模拟结果将几何扩散取为R_hypo^－1；本文的研究区域位于日本火山弧前，传播路径项中非弹性衰减部分按照文献[27]对火山弧前弧后的研究结果，采用其拟合出的弧前区域品质因子用于本次地震动模拟，形式为$Q(f)= \begin{cases}300(f / 0.1)^0 & f \leqslant 0.64 \mathrm{~Hz} \\ 150(f / 1.0)^{1.30} & f \geqslant 3.66 \mathrm{~Hz}\end{cases}$的模型，当频率0.64 Hz＜f＜3.66 Hz时，Q(f)由分段点处对应值连线构成的一次函数表示。

井下台的场地效应考虑了包括地壳放大和表示局部场地高频快速衰减的κ₀值，在此基础上，井上台还考虑了局部场地效应，其中包含局部土层放大效应以及对井上台κ₀值的修正。

近地表地壳放大通过改进后的四分之一波长法^[28]计算各个台站所在地区的波速模型获得，KiK-net提供了部分台站对应井下土层结构的数据资料，但不足以使用四分之一波长法计算全部台站的近地表地壳放大。对于没有提供井下波速结构的台站，本文结合文献[29]提出的日本综合速度模型作为补充，选取台站所在划分区块的地下结构速度模型，计算出每个场点对应的近地表地壳放大因子，全部结果见图 3。

与场地效应相关的高频衰减部分，本文使用文献[30]提出的κ₀值作为加速度幅值谱在高频部分的衰减，计算κ₀所用地震动数据为2018—2021年内目标区域及邻近台站记录到的90组井下水平向地震动加速度记录，所用记录均来自KiK-net台网，且经过基线校正和滤波处理。按照文献[31-32]的处理流程，取每条记录加速度时程S波到时到能量持时80%的部分进行傅里叶变换。根据文献[33-34]的观点，将傅里叶幅值谱的衰减起始频率f_E控制在5 Hz≤f_E＜10 Hz，衰减终止频率f_max控制在10 Hz≤f_max≤20 Hz，且起始频率与终止频率的间隔>10 Hz。本文所用起始频率为5 Hz，终止频率为20 Hz，对每组傅里叶谱中间的高频衰减部分使用最小二乘拟合，得到EW分量和NS分量的κ关于震源距R的关系式κ=0.047 5+0.000 033 7R和κ=0.048 3+0.000 016 9R(见图 4，圆点表示根据观测记录计算出的κ值)。当R=0 km时κ记为κ₀，表示目标研究区域的局部场地高频衰减效应，本文使用2个水平分量κ₀的算数平均值0.048作为此次模拟的局部场地高频衰减输入参数。

由于KiK-net具有丰富的井上井下观测记录，本文采用物理意义明确的井上井下谱比法对各个台站的局部场地效应进行估计。首先，对所有基线校正和滤波后的井上井下水平向地震动记录进行傅里叶变换，并使用Konno-Ohmachi窗对其幅值谱进行平滑处理。对平滑后的两分量傅里叶幅值谱求几何平均值，再用同一个台站的井上和井下幅值谱均值做谱比，最终得到25个台站的井上井下谱比，将其作为本次模拟的局部场地效应参数，全部结果见图 5。

本次随机有限断层法模拟所用震源、路径和场地部分的输入参数见表 2。

3 模拟结果及分析

考虑到模拟结果的稳定性和随机方法的不确定性，使用有限断层法对全部井上井下台站记录各进行了30次模拟，得到相应阻尼比5%的A_ps和加速度时程，并将每组30次模拟A_ps的几何平均值作为目标台站的A_ps模拟结果。分别选取地表水平向A_pg＞60 cm·s^-2的井上台观测记录和其对应井下台记录，将其A_ps值与模拟值进行比较，并在井上部分的对比中加入地震动预测方程(GMPEs)A_ps预测平均值(以下简称预测值)，其中GMPEs为文献[7]提出的日本地区俯冲带板内水平向地震动衰减关系(以下简称Zhao16)，Zhao16所用场地见表 1，比较结果见图 6。其中，图 6(a)、(b)分别为井下和井上的比较结果。整体来看，井上和井下台模拟记录与观测记录得到的A_ps谱形高度相似，低频部分模拟稍微高估。由图 6(b)中可以看出，Zhao16也可以较好地描述代表地面运动的井上台记录频谱特征，但Zhao16的A_ps预测值在低频上略低于观测值和模拟值，模拟得到的A_ps值更加接近观测值，由于以上台站处于关东盆地，盆地效应及厚覆盖层可能对实际地震动的低频造成一定程度的放大。为验证这种放大与盆地相关，本文在山地区台站中选取海拔>300 m的台站(KNGH21、KNGH18、TKYH13、KNGH20)加入比较，见图 6(c)。可以看出，对于山地区台站，Zhao16的A_ps值并没有出现低频段低估的现象。

对模拟A_ps的评价可用观测记录与模拟记录的A_ps值对数残差表示，在0.1~10 Hz的频带范围内，全部井上井下台站的模拟记录与水平向观测记录A_ps对数残差结果见图 7，图 7(a)为井下结果，图 7(b)为井上结果，实线为对数残差平均值，越接近0表示模拟结果越接近观测值，阴影区域为对数残差的一倍标准差范围。可以看出在所选频段内，所有台站模拟记录的A_ps对数残差平均值均在±0.4以内，其中0.6~10 Hz频段基本在±0.2以内；A_ps对数残差的一倍标准差主要在±0.5范围内。由此可以看出，模拟记录能很好地反映本次地震动的频域特征。

检验模拟与观测记录加速度时程的匹配程度，同样选取地表水平向A_pg＞60 cm·s^-2的台站，从模拟30次得到的水平向加速度时程中挑出A_pg与观测记录A_pg几何平均值相近的井上记录进行对比，同时为了检验70%能量模型的表现，本文将广泛使用的由文献[35]开发的地震动持时模型(以下简称AB98)加入对比，见图 8。可以看出相比于AB98，70%能量模型得到的模拟加速度时程与观测记录在强震动段更相似，其持时也更接近地震动观测记录，模拟地震动记录在时域上的表现良好。

峰值地面加速度(A_pg)是计算仪器烈度的常用参数，与烈度分布有很强的相关性。为研究地表A_pg特征及其随距离的衰减特征，本文将25组井上台站水平向A_pg观测值、A_pg模拟值与Zhao16的A_pg预测值进行比较。地表A_pg观测值使用井上水平向A_pg的几何平均值，将模拟30次井上记录获得的A_pg值进行几何平均作为地表A_pg模拟值，Zhao16预测值所用场地见表 1。比较结果见图 9。可以看出，观测值和模拟值基本落在Zhao16一倍标准差范围内，表明模拟记录与Zhao16均能很好地描述本次俯冲带板内地震的A_pg衰减规律。此外，Zhao16的A_pg预测值普遍高于地表A_pg观测值，地表模拟值更接近观测值，大多数模拟值与观测值匹配良好。由于A_pg主要受地震动高频成分的影响，关东盆地的厚软土覆盖层可能会抑制地震动的高频成分，从而导致Zhao16在A_pg预测上的高估。相比于Zhao16模型，随机有限断层法可以更灵活地考虑研究区域的局部场地效应。

进一步研究本次地震动的A_pg空间分布，基于井上水平向观测记录和模拟记录，利用连续曲率样条插值法方法^[36]绘制地表A_pg等值线，其中地表A_pg模拟值为30次井上模拟值的几何平均值，结果见图 10，可以看出观测与模拟的A_pg空间分布相似，模拟记录得到的A_pg在东京地区略有低估，整体来看具有与观测记录较为一致的空间分布特征。

4 结论

1) 井上井下台模拟记录与观测记录的A_ps幅值整体对应较好，虽然在0.1~0.6 Hz的低频段内，模拟A_ps普遍存在高估现象，但在中高频段上(0.6~10 Hz)，A_ps模拟值与观测值高度吻合。此外，Zhao16也能很好地描述本次地震动的频谱特征，但盆地区台站的预测值在低频出现了一定程度的低估，而山地区台站并未出现上述情况，这可能是由于盆地效应及厚覆盖层对实际地震动造成了低频的放大。

2) 模拟加速度时程与观测记录时程波形相似，基于70%能量持时回归得到的地震动持时模型可以有效改善模拟记录与观测记录在强震动段不匹配的情况，模拟加速度记录的持时也更接近真实记录。

3) 地表A_pg模拟值与观测值匹配良好，随机有限断层法和Zhao16均可以很好地描述本次俯冲带板内地震动的衰减特性。其中，Zhao16的A_pg预测值相比于模拟值和观测值普遍出现高估，这可能是目标研究区域软土层对地震动高频抑制引起的。

4) 井上模拟记录与观测记录获得的A_pg等值线相似度较高，通过随机有限断层法获得的模拟地震动在A_pg空间分布上与此次观测到的俯冲带板内地震引起的地震动有很好的适配性。