2017-08-08四川省九寨沟县发生M_S7.0地震。关于九寨沟地震强地面运动数值模拟的研究取得了一些成果^[1-4]，但仍有不足之处，如应用复合震源模型进行模拟及其震源参数敏感性分析的研究较少、随机有限断层法计算地震动分量方向单一、地震峰值加速度拟合精度不高等。为此，本文应用复合震源模型对九寨沟M_S7.0地震强地面运动数值模拟展开研究，并探讨不同模型参数对合成地震动的影响规律，为九寨沟地区工程场地地震安全评价提供参考依据。

1 复合震源模型基本理论

复合震源模型是用于计算宽频带强地面运动的有限源模型^[5]。该模型通过描述震源的运动学特征，结合地震波在层状介质中的传播过程，将广义反射和透射系数法^[6]得到的格林函数与震源时间函数进行卷积，计算得到宽频带地震动。模型中，强震震源被认为是大量具有恒定应力降圆形子震的叠加。首先将断层面划分成一系列尺寸相同的正方形子断层网格，之后用一定数量和相应半径大小分布服从Frankel自相似模型^[7]的圆形子震填充该断层面，并允许其部分重叠。圆形子震相当于一系列小震，在震源破裂过程中相继被激发，最终释放的能量在断层面上不断被叠加，达到模拟强震的目的。设圆形子震半径为R、数目为N，则有：

$ \frac{\mathrm{d} N}{\mathrm{~d}(\ln R)}=P R^{-D} $

(1)

式中，分形维数D=2，控制断层面上圆形子震数量的分布及圆形子震的上升时间，P为比例常数：

$ P=\left\{\begin{array}{l} \frac{7 M_0}{16 \Delta \sigma} \frac{3-D}{\left(R_{\max }^{3-D}-R_{\min }^{3-D}\right)}, D \neq 3 \\ \frac{7 M_0}{16 \Delta \sigma} \frac{1}{\ln \left(R_{\max } / R_{\min }\right)}, D=3 \end{array}\right. $

(2)

式中，M₀为地震矩，Δσ为应力降，R_max为最大圆形子震半径，R_min为最小圆形子震半径。对式(1)进行积分，可得半径大于R且小于R_max圆形子震的分布数量N(R)：

$ N(R)=\frac{P}{D}\left(R^{-D}-R_{\max }^{-D}\right) $

(3)

式中，最大圆形子震半径R_max取断层宽度的一半^[8]，最小圆形子震半径R_min取R_max/20，每个圆形子震的地震矩取决于该圆形子震半径R_i和应力降Δσ。第i个圆形子震的地震矩M₀ⁱ为：

$ M_0^i=\frac{16}{7} \Delta \sigma R_i^3 $

(4)

震源时间函数由位于断层面上的一系列圆形子震震源时间函数叠加得到，每个圆形子震的破裂时间由其到震源中心的距离和破裂速度决定。在给定地震矩、应力降、破裂速度及子断层尺寸等震源参数的情况下，断层从震源以恒定的速度v_r破裂，当破裂经过圆形子震中心时，圆形子震开始辐射能量。其中，第i个圆形子震的震源时间函数$ \dot{M}_0^i$为：

$ \dot{M}_0^i=\left(2 {\rm{ \mathsf{ π} }} f_c^i\right)^2 M_0^i H(\tau) \tau \mathrm{e}^{-2 {\rm{ \mathsf{ π} }} f_c^i \tau} $

(5)

式中，τ=t－t_i，t为整个断层的破裂时间，t_i为第i个圆形子震的发震时间，t_i=X/v_r，X为圆形子震圆心到震源的距离，v_r为破裂速度，H(τ)为阶梯函数，f_cⁱ为第i个圆形子震的角频率：

$ f_c^i=\frac{2.34 \beta}{2 {\rm{ \mathsf{ π} }} R_i} $

(6)

式中，β为剪切波的速度。断层面上滑动位移量的空间分布由圆形子震滑动位移量叠加得到，第i个圆形子震的滑动位移量ΔU_i为^[1]：

$ \Delta U_i=C \frac{\Delta \sigma}{\mu} R_i $

(7)

式中，C为依赖断层类型的常数，μ为介质剪切模量。

2 地震数据与震源参数

本文研究区域(103°~105°E，31°~35°N)位于巴颜喀拉块体，该块体近年来地震频发。九寨沟地震后，中国地震台网中心等机构公布了其震源位置、地震矩、震源深度、断层节面的走向、倾角和滑动角等震源参数(表 1)，其中，根据美国地质调查局的研究，断层节面Ⅱ为实际断层面^[1]，因此本文取节面Ⅱ进行研究。

基于孟令媛等^[9]对余震分布情况及断层面上滑动位移的反演结果，建立长30 km、宽20 km的断层模型，并将其划分为150个2 km×2 km的正方形子断层。断层面上圆形子震的分布如图 1所示，其中圆形子震地震矩的总和等于总地震矩。通过式(7)计算得到断层面上滑动位移量的不均匀分布，如图 2所示。九寨沟地震震中距150 km区域范围内台站信息如表 2所示，九寨沟地区介质速度信息(图 3)来自全球地壳模型Crust1.0(https://igppweb.ucsd.edu/~gabi/crust1.html)，台站位置分布如图 4所示。对所有台站观测数据进行基线校准和滤波处理，带通滤波范围为0.1~20 Hz。

3 参数敏感性分析

在运动学震源模型中，断层的破裂面积、长度、宽度、埋置深度、走向、倾角和滑移角等是描述断层特征的参数；地震矩、应力降、平均滑动时间、平均破裂速度、破裂方式和初始破裂点位置等是描述平均破裂过程的参数。这些震源参数控制着地震动强度和分布^[10]，分析这些参数的影响规律可为地震动合成精度提供保证。本文首先以无强震观测记录为前提，研究地震矩、应力降、破裂速度、子断层尺寸、震源深度、断层走向、倾角和滑移角等8个震源参数在所有合理取值范围内对复合震源模型计算结果的影响程度；然后根据九寨沟地震观测数据确定适合的震源参数研究范围，并进行敏感性规律分析，在其他参数保持不变的条件下，仅改变单个参数变量，利用评价指标比较分析计算结果，总结震源参数对地震动计算的影响规律。

3.1 评价参数和评价指标

为表征震源参数对台站NS和EW向峰值加速度(peak ground acceleration，PGA)的影响，用误差值E表示计算的PGA相对于观测的PGA的变化：

$ E=\left|A_{\mathrm{OBS}}-A_{\mathrm{SYN}}\right| / A_{\mathrm{OBS}} \times 100 \% $

(8)

式中，A_OBS为观测的PGA，A_SYN为计算的PGA。

地震动强度主要用PGA、峰值速度或谱加速度等地震动参数表达。本文以PGA为例，根据王德才等^[11]的参数敏感性分析模型，通过计算模型偏差和模型标准差讨论不同震源参数对PGA的影响。

用模型偏差B表征计算值和观测值之间的差异，偏差越大，计算值越偏离观测数据：

$ B=\frac{1}{M} \sum\limits_{j=1}^M \ln \left(\frac{A_{\mathrm{OBS}}}{A_{\mathrm{SYN}}}\right) $

(9)

式中，M为总台站数。

用模型标准差S描述数据离散程度，标准差越大，数据波动性越强：

$ S=\left\{\frac{1}{M} \sum\limits_{j=1}^M\left[\ln \left(\frac{A_{\mathrm{OBS}}}{A_{\mathrm{SYN}}}\right)-B\right]^2\right\}^{1 / 2} $

(10)

3.2 参数敏感性分析结果

首先选取地震矩、应力降、破裂速度、子断层尺寸、震源深度、断层走向、倾角和滑移角等8个震源参数的最大合理取值范围，计算得到不同震源参数数值下台站NS和EW向的PGA，然后以模型偏差B和标准差S为衡量标准，分析8个震源参数对复合震源模型计算结果的影响程度。

地震矩根据Hanks等^[12]提出的地震矩与矩震级的关系，在强震6~8级范围内取值计算，间隔为1×10¹⁸ Nm；应力降参考王振宇等^[13]的取值范围10~200 bar，间隔为5 bar；断层的平均破裂速度约为剪切波速的0.6~0.9倍^[14]，本文选取此区间，并设定区间间隔为0.05；子断层尺寸取值分别为1 km×1 km、2 km×2 km、6 km×6 km、8 km×8 km和10 km×10 km；震源深度范围选定为0~30 km，取值间隔为1 km；断层走向、倾角和滑移角依据理论范围分别取0°~360°、0°~90°和－180°~180°，取值间隔分别为10°、5°和10°。

8个震源参数取值范围和梯度各不相同，本文通过式(9)和式(10)计算出相应PGA的偏差和标准差，然后将每个震源参数下PGA偏差、标准差的最大值及最小值作差，消除A_OBS实测数据的影响，得到每个震源参数下PGA偏差和标准差的变化幅度，结果如图 5所示。由图 5可知，地震矩、应力降、破裂速度、子断层尺寸和震源深度对PGA偏差和标准差的影响明显，其中地震矩和应力降影响更突出；断层走向、倾角和滑移角对PGA偏差和标准差的影响极不明显。根据地震学理论，地震矩直接决定地震释放的能量大小，地震矩越大，地震波携带的能量越多，使得整体加速度计算值增加；应力降表示断层上某点破裂前后的应力强度变化，应力降越大，断层上释放的应力越多；破裂速度决定断层面圆形子震的发震时间，间接影响整个断层面的破裂过程；子断层尺寸的划分会对断层面上滑移量分布产生很大的影响；震级相同的地震，震源越浅，造成的破坏越严重。

基于8个震源参数的影响规律，结合表 1不同研究机构和众多学者^{[2-4, 8]}提供的震源参数，确定九寨沟地震震源参数研究范围，进一步研究8个震源参数对九寨沟11个台站NS和EW向PGA的影响大小和规律，取值间隔如图 6所示。

从图 6看出，地震矩、应力降、破裂速度、子断层尺寸和震源深度对台站PGA的影响大；断层走向、倾角和滑移角对台站PGA的影响小。说明在应用复合震源模型计算九寨沟地震时，即使8个震源参数取值范围有所限制，没有取到理论范围内的所有值，但其对计算结果的影响规律基本与图 5一致。

图 7为地震矩、应力降、破裂速度、子断层尺寸和震源深度等5个震源参数对九寨沟地震地震动模拟的影响规律，其中实心点表示本文取值。由表 3和图 7可知，受地震矩影响，PGA的偏差变化范围超过1，PGA的偏差和标准差随地震矩的增加而减小，其中PGA的偏差下降趋势近似为线性变化，地震矩为7×10¹⁸ Nm时，PGA的偏差基本为0，同时PGA的标准差也较小。应力降越大，台站整体PGA越大，PGA的偏差越小，应力降为40 bar时，PGA的偏差接近于0，这与俞瑞芳等^[15]应用随机有限断层法得到PGA随着应力降增大而逐渐增加的结论类似。破裂速度越大，PGA的偏差越小，当破裂速度等于剪切速度0.8倍(2.96 km/s)时，偏差接近于0，此时标准差也最小，与傅磊等^[16]应用随机有限断层法和李宗超等^[4]应用经验格林函数法在地震模拟过程中选取的破裂速度与剪切速度的比例关系一致，PGA的偏差变化范围为0.9。子断层尺寸越小，PGA的偏差越小，最佳尺寸为2 km×2 km，此时标准差也为最小值，与Dang等^[1]利用随机有限断层法的研究结论相同。震源深度在20 km处PGA的偏差几乎为0，同时标准差最小，随着震源深度的增加，PGA的偏差增长趋势变缓；而且如图 6所示，震源深度为25 km时，大部分台站的峰值加速度很小，与刘启方等^[10]运用震源位错模型得到的研究结果相似。

断层走向、倾角和滑移角中，PGA的偏差最大变化幅度为0.08，标准差最大变化幅度仅为0.1，对计算结果影响不明显，与王德才等^[11]采用混合方法^[17]和随机方法^[18]研究的结果一致。

3.3 地震输入参数组合

由上文分析结果可知，震源参数整体对PGA偏差影响幅度大，对标准差影响幅度相对较小，因此以PGA的偏差为主要衡量标准，标准差为次要标准，得到一组与观测记录拟合效果较好的震源参数。从表 1各机构提供的数据中选取地震矩、破裂速度、子断层尺寸、震源深度、断层走向、倾角和滑移角；应力降参考Sun等^[3]的结果，品质因子Q(f)和高频衰减kappa值参考孟令媛等^[9]的结果，具体数值如表 4所示。

4 地震动模拟结果分析

根据式(8)计算11个台站NS和EW向PGA的相对误差，并分别与Dang等^[1]和李宗超等^[4]的计算结果进行对比(图 8)。可以看出，11个台站NS和EW向PGA相对误差均小于50%。62SHW站NS和EW向PGA相对误差最小，均小于3%；51HSL站NS向、51HSD站NS和EW向PGA误差较大，且计算值均大于观测值，主要原因可能是台站震中距较大，导致传播过程中散射的复杂性增加，同时，台站局部山体地形引起的场地效应也会增加计算误差。11个台站NS和EW向PGA相对误差均在20%以内，说明所有台站PGA与实际数据相比误差小。与Dang等^[1]和李宗超等^[4]的计算结果相比，除李宗超等^[4]51JZY站EW向、62ZHQ站NS向和62MXT站NS向外，本文各台站PGA相对误差均较优。为方便直观地进行地震动持时和波形对比，所有时程均截取45 s时间窗(图 9)，台站加速度记录数据中最小PGA为2.75 cm/s²，为避免地震动持时出现零持时这种不合理状态，选取Bolt持时的阈值为0.5倍的PGA作为计算相对Bolt持时的阈值^[19]。从图 9看出，个别台站误差较大，但大多数台站地震动持时基本保持一致。如图 10所示，合成反应谱整体变化趋势与观测反应谱一致，说明11个台站的峰值加速度、地震动持时和反应谱整体拟合效果较好，与观测数据较为符合。

部分台站地震动三要素计算结果精度偏低的可能原因是：1)没有充分考虑地形因素对地震动的影响，部分台站如51HSL、51HSS和51HSD等位于山体上，山体地形的高度及底部延展距离对断层动力学破裂过程影响较大，进而影响到相应的地面地震动分布^[20]，山体的坡角、地震动入射角度会对山体地震动产生抑制作用，导致台站PGA变小^[21]；2)本文在计算地震动时采用各向同性散射模型，但实际地震波在地壳介质中的传播是各向异性的，而且随着传播路径的增加，散射的影响也会逐渐加强；3)有些震源参数的选取进行了简化，在地震动模拟过程中假定某些参数(如应力降和破裂速度)为常数，忽略了参数与时间和空间的相关性；4)多数台站低频段反应谱模拟结果高于观测值，可能是由于本文模型在破裂过程中没有考虑到圆形子震倾角的可变性，或是在合成格林函数中未纳入小尺度地壳非均匀性的散射效应^[22]。

PGA是表征地震能量的重要指标之一，能够为建筑物抗震设防提供合理的参考依据。综合现有11个台站PGA的观测值和计算值可以得到九寨沟地震PGA空间分布特征，如图 11所示。由图 11可知，随着震中距增大，NS和EW向PGA从震源中心呈椭圆形向四周衰减扩散，长轴方向大致呈SW走向，这可能与断层破裂的方向性有关，且2个方向PGA最大值均出现在51JZB站；在51MXD和51PWM站之间有一个PGA值未连通的区域，可能是由于该区域缺少合适的台站记录所致。总体而言，PGA计算值和观测值的空间分布特征较为一致。

5 结语

1) 在最大合理变化范围内，相对于断层走向、倾角和滑移角，地震矩、应力降、破裂速度、子断层尺寸和震源深度对地震动计算结果影响较大。

2) 震源参数敏感性分析规律为九寨沟地震地震动模拟的参数选取提供了依据，进而得到一组合理的震源参数组合，同时能够获得九寨沟地区内任意一点PGA的变化区间。这一结论可以进一步推广到其他缺乏强震观测记录的场地，计算出该区域强度特征的上限值，为建筑物抗震设防提供依据。

3) 复合震源模型强地震动模拟了特定圆形子震在断层上的破裂过程、断层面上滑动位移的不均匀空间分布以及震源时间函数，计算结果可以再现九寨沟地震的地震动特征，模拟区域台站的地震动PGA、地震动持时、反应谱和PGA分布特征与观测值拟合较好，反应谱在高频段拟合精度更高。