引言

2017年8月8日21时19分，四川省北部阿坝州九寨沟县发生M_S 7.0地震。根据中国地震台网中心发布的信息，本次地震的震中位于33.20°N、103.82°E，震源深度10km。震中在九寨沟县城以西约39km处，东偏北方向距离陇南市105km，南距成都市285km。本次地震的最大烈度为Ⅸ度，共造成25人死亡，525人受伤，6人失联，受灾人数超过17万(中国地震局，2017)。地震导致7万多间房屋不同程度受损，其中76间房屋倒塌。此次地震获取到强震观测记录的台站数量有限，其中距离震中最近的九寨百河台站(台站编号51JZB)，南北向地震记录的峰值加速度(PGA)达到了0.19g，是本次地震中观测到的最大加速度(图 1)。

图 1  断层及强震台站分布图

Figure 1.  Location of the causative fault and strong-motion stations

下载: 全尺寸图片 幻灯片

大量的震害资料表明，地震动是造成结构破坏、诱发地质灾害(例如场地液化、地表裂缝和滑坡)的主要驱动力。地震动记录对于研究地震危险性、地震引发的地质灾害和结构抗震设计具有重要的意义。在缺少强震观测记录的情况下，可以通过有限的实际观测记录，借鉴科学合理的方法，综合考虑震源、传播路径和场地条件人工合成地震动，其结果可作为结构动力反应分析和灾后抗震设计的依据。

本文采用Wang等(2015)改进的随机有限断层方法合成2017年九寨沟地震的近场地震动。首先选取地震的有限断层模型和随机方法的输入参数，并通过与实际记录对比估计地震的应力降；采用随机有限断层法合成震源附近台站的加速度时程；并将得到的模拟地震动与实际记录和烈度图进行对比，验证结果的有效性。

1.   随机有限断层法

根据有限断层法的思想，断层破裂可以离散为N个子断层，每个子断层可以近似认为是1个点源。破裂过程从震源开始，逐渐向断层边缘扩散，当传播到每个子断层时，子断层发生破裂并产生地震波。将子断层产生的地震波按照到达的先后次序进行叠加，即可得到整个断层在某个观测点产生的地震动。Beresnev等(1997)将有限断层模型和随机点源方法(Boore，1983)结合，提出了随机有限断层法，可以用来模拟断层附近的高频地震动。Motazedian等(2005)进一步发展了这种方法，提出了基于动态拐点频率的随机有限断层法(EXSIM方法)，解决了总能量受到子断层尺寸影响等问题。Wang等(2015)改进了Motazedian等(2005)的方法，采用一种双拐点震源谱，能够反映破裂过程中累积地震矩对震源谱的影响。

本文采用改进的随机有限断层法(Wang等，2015)模拟九寨沟地震动分布，综合考虑震源、传播路径、场地的影响，有限断层模型中第j行、第i列的子断层在场点产生地震动的加速度傅里叶谱可以表示为：

其中，f为频率(Hz)；R_ij为第j行、第i列子断层到观测点的距离(km)；$E({M_{0ij}}, {f_{0ij}}, f) $代表震源谱(cm/s)；$G({R_{ij}}) $为几何衰减函数；$ {H_{ij}}$为低频比例因子；$ Q(f)$为地壳品质因子；$ \beta $为剪切波速(km/s)；$ S(f)$为与频率相关的场地放大因子；$ \kappa $为高频衰减因子。第j行、第i列子断层的震源谱按照式(2)计算：

其中，C为常数，按照Motazedian等(2005)给出的公式计算；${M_{0ij}} $是第j行、第i列子断层的地震矩(dyne-cm)；a和b是控制震源谱形状和幅值的动态参数，可由断层已破裂部分的地震矩确定：

其中，$ {M_0}(t)$表示在时间t断层上已破裂区域的地震矩，可根据地震的总地震矩，以滑动分布为权重分配求得。式(2)中拐点频率f_0ij按照式(4)可表示为：

其中，$ \mathit{\Delta }\sigma $代表地震的应力降(bar)。

除震源谱模型外，式(1)中其他几项函数代表传播路径和场地条件对地震动傅里叶谱的影响，如几何衰减、场地放大因子和高频衰减因子。这些函数在相关文献(Boore，1983；Motazedian等，2005)中有详细的介绍，本文不再赘述。另外，Wang等(2015)采用上述改进的随机有限断层法合成了美国1994年Northridge地震的近场地震动，结果和实际记录较为符合，验证了方法的有效性。

2.   九寨沟强震的地震动模拟

2.1   震源模型及输入参数

基于远场P波和SH波的波形数据，王卫民等(2017)反演了九寨沟M_S 7.0地震的震源破裂过程，得到了震源机制以及发震断层滑动分布的初步结果。结果显示地震矩为6.7×10²⁵dyne-cm，对应的矩震级为M_W 6.5；断层的走向、倾角和滑动角分别为148.5°、68.9°和-3.1°，属于左旋走滑型地震。基于此结果，本文得出的断层滑动分布见图 2，从图中可以看出，断层上的破裂滑动大部分集中在震源附近区域，即深度6—12km范围内，最大滑动量为85cm。此外，在震源东南约16km的位置，也有幅值小于40cm的破裂滑动。

图 2  断层滑动分布

Figure 2.  Slip distribution on the fault plane

下载: 全尺寸图片 幻灯片

表 1列出了随机有限断层方法的输入参数，其中大部分是其他学者根据四川及周边地区的大量地震记录分析得到的，反映了该地区地震动衰减特性。其中，几何衰减采用喻烟(2012)提出的汶川地震区地震动估计的经验模型，高频衰减因子采用Sun等(2013)分析汶川余震记录得到的结果，与距离相关的持时采用Atkinson等(1995)提出的分段模型。基于对九寨沟地震记录S波的傅氏谱反演结果，王宏伟等(2017)确定了这次地震震中区域的地壳品质因子为$ Q(f) = 84.9{f^{0.71}}$，体现了较强的地震波衰减特性，本文的地震动模拟中采用了此结果。对于场地效应，地表到30m深度的平均剪切波速(V_S30)通常被用作场地类别划分的依据，例如美国地震减灾项目(NEHRP)。喻烟(2012)根据我国四川与甘肃地区的钻孔数据，假定钻孔底部到30m深度的波速与孔底一致，计算了台站场地的V_S30。本文根据这些结果，确定了台站对应的NEHRP场地类别，并采用Boore等(1997)针对不同类别场地提出的场地放大因子。对于所有台站，断层破裂的脉冲面积百分比采用50%。本文通过对比模拟和观测的地震动，按照误差最小原则估算本次地震的应力降参数。

表 1  模型输入参数

Table 1.  Input parameters of model

参数	值或表达式
断层走向和倾角	148.5°，68.9°
断层尺寸	54km×24km
子断层尺寸	2km×2km
剪切波速(β)	3.5km/s
介质密度(ρ)	2.8g/cm3
破裂速度(VR)	0.8×β
应力降	4.0MPa
品质因子	$ Q(f) = 84.9{f^{0.71}}$
高频衰减因子(κ)	0.015
路径持时T=dR，参数d	0.0(R≤10km) -0.16(10km＜R≤70km) -0.03(70km＜R≤130km) 0.04(R＞130km)
几何衰减G(R)=Rb，参数b	-1.0(R≤150km) -0.5(R＞150km)
场地放大因子	参见Boore等(1997)
脉冲面积百分比	50%

下载: 导出CSV 

| 显示表格

2.2   应力参数估计

首先，选取震源距小于150km、记录的PGA大于10cm/s²的台站，将这些台站取得的记录作为估计应力降的依据。表 2给出了选取的8个台站信息，包括经纬度、震源距、V_S30和相应的NEHRP场地类别。估计应力参数主要采用“试错法”。考虑1.0—6.0MPa的应力降，增量为0.5MPa，对于每个应力降，采用2.1节中的震源模型和输入参数合成8个台站的地震动加速度时程。考虑随机方法Gauss白噪声幅值的波动性，每个台站模拟30次，计算平均的反应谱。定义总体平均误差如下：

表 2  用于估计应力参数的台站信息

Table 2.  The information of stations selected for the estimation of stress drop

台站名称	台站代码	纬度/°N	经度/°E	震源距/km	VS30/m·s-1	NEHRP场地类别
岷县	62MXT	34.4	104.0	135.1	301	D
迭部	62DIB	34.1	103.2	117.8		C
沙湾	62SHW	33.7	104.5	82.5	384	C
九寨百河	51JZB	33.3	104.1	27.2	328	D
九寨勿角	51JZW	33.0	104.2	39.7	428	C
九寨永丰	51JZY	33.2	104.3	42.3	321	D
茂县叠溪	51MXD	32.0	103.7	133.7	268	D
平武木座	51PWM	32.6	104.5	89.8	373	C
注：迭部台由于缺少钻孔资料和VS30，场地类型假定为NEHRP-C。

下载: 导出CSV 

| 显示表格

其中，$ S{A_{{\rm{obs}}, i}}(f)$，$ S{A_{{\rm{sim}}, i}}(f)$分别表示第i个台站记录和模拟的加速度反应谱(cm/s²)，N表示台站数量。对于每个应力降，根据合成结果计算出总体平均误差。需要说明的是，由于随机方法通常用于模拟高频地震动，不一定能够体现出低频地震动的特征，因此本文考虑的频率范围是1—30Hz。根据上述准则，得到误差最小的应力降为4.0MPa。王宏伟等(2017)通过反演本次地震记录S波的傅里叶谱，确定了地震的应力降为3.854MPa，与本文通过“试错法”得到的结果较为接近。原因之一是本文的模拟采用了同一篇文献(王宏伟等，2017)中基于反演地震记录确定的品质因子$Q(f) $。图 3(a)给出了4.0MPa应力降误差随频率的变化以及误差±1个标准差的范围。从图中可以看出，在1—30Hz频率范围内，误差基本较小，绝对值小于0.2。作为对比，我们也计算了EXSIM方法得到的误差(图 3(b))，其中所有输入参数(包括应力降)和本文采用的方法相同。通过对比可以明显看出，EXSIM方法的误差在频率大于3Hz时明显大于本文的结果；同时在1Hz附近，EXSIM方法的平均误差为负，表示模拟地震动在该频段高估了观测地震动的反应谱幅值。由此可见，相比EXSIM方法，采用Wang等(2015)的方法得到的结果与本次地震记录更相符。

图 3  2种有限断层方法得到的误差随频率分布

Figure 3.  Residual v.s. frequency obtained from two finite-fault methods

下载: 全尺寸图片 幻灯片

2.3   模拟地震动的时程和反应谱

按照上述随机有限断层方法，采用表 1的模型输入参数，分别模拟了6个强震台站的加速度时程。在本文选取的台站中，虽然大部分都有剪切波速的钻孔数据，但多个场地的钻孔未达到基岩，如果进行土层反应分析需要进行一定的假设或者延拓剪切波速剖面，存在较大的不确定性。因此，在模拟中根据台站的NEHRP场地类型，采用Boore等(1997)提出的NEHRP-C、NEHRP-D场地的放大系数，近似考虑场地效应的影响。图 4给出了6个台站模拟和观测地震动的加速度时程对比。为了便于进一步定量比较，表 3列出了模拟和观测的地震动PGA值。需要说明的是，随机有限断层法得到的地震动对应于任意水平分量。

图 4  模拟和观测地震动的加速度时程对比

Figure 4.  Comparison of acceleration time histories between synthetic and recorded ground motions

下载: 全尺寸图片 幻灯片

表 3  模拟和观测地震动的峰值加速度

Table 3.  The PGA from synthetic and recorded ground motions

台站代码	PGA(记录)/cm·s-2		PGA(模拟结果)/cm·s-2
	EW分量	NS分量
51JZB	129.50	185.02	174.78
51JZW	73.84	91.70	134.90
51JZY	45.77	66.71	85.10
51MXD	12.21	23.32	11.13
62SHW	18.57	20.48	19.27
51PWM	18.60	20.87	21.26

下载: 导出CSV 

| 显示表格

通过对比模拟和记录的地震动可以发现，在峰值加速度最大的台站51JZB，模拟地震动的PGA值介于东西向和南北向记录的PGA值之间，与实际记录较为符合；在震源距较大的2个台站62SHW和51PWM，拟地震动的PGA与实际记录也很吻合；另外，51JZW和51JZY模拟地震动的PGA稍大于实际记录，51MXD模拟地震动的PGA小于地震记录，原因可能有以下3点：①采用Boore等(1997)的场地放大系数不能充分考虑场地条件对地震动的影响，如非线性土层反应和不规则地形影响；②本文采用的震源滑动分布模型是通过反演远场波形数据得到的，没有考虑近断层地震动和同震位移场，不一定能够完全反映出震源特征；③本文模拟采用的模型参数大多是基于分析汶川地震数据得到的，这些参数也会存在不同地震、不同地区之间的差异。

从图 4可以看出，大部分台站模拟地震动的持时在不同程度上小于实际记录。需要指出，随机有限断层法模拟的实际上是地震动的S波。一方面，在近场地震记录中P波与S波2种成分相互叠加，会导致记录的持时大于模拟地震动；另一方面，部分台站(如51JZY、51MXD和51PWM)的地震记录中面波成分很丰富，无法通过随机方法模拟出来；除此之外，近地表土层反应和不规则地形也可能造成地震记录的持时较长。解决这些问题需要更复杂的地震动模拟方法。

图 5给出了6个台站模拟和观测地震动的加速度反应谱。通过对比可以看出，在本文考虑的周期(T)范围内(0.03—1.0s)，大部分台站模拟结果的反应谱与实际记录较为符合。51JZW台站模拟地震动的反应谱在短周期(T＜0.1s)大于实际记录的反应谱，51MXD台站模拟地震动低估了T＞0.1s地震记录的反应谱。这些差异说明本文采用的随机有限断层方法和模型输入参数仍有可改进之处。综合考虑随机性方法与确定性方法的混合方法，能够弥补随机性方法在模拟低频地震动和面波等方面的不足，但是需要提供更多的信息，如地壳介质波速模型等。

图 5  模拟和观测地震动的反应谱对比

Figure 5.  Comparison of response spectra between synthetic and recorded ground motions

下载: 全尺寸图片 幻灯片

2.4   模拟的PGA分布图

为了进一步研究模拟地震动加速度峰值的空间分布，我们在发震断层附近以5km为间隔选取了1548个网格点。对于每个点，采用2.1节介绍的断层模型和输入参数模拟地震动，其中场地类型假定为NEHRP-C，并采用Boore等(1997)对此类场地提出的放大系数。考虑随机模拟结果的不确定性，每个点模拟20次，最后取平均的PGA作为结果。根据所有网格点的模拟PGA，画出该区域的PGA分布图(图 6(a))。从图中可以看出，模拟PGA在近断层区域的分布明显受到断层位置以及滑动分布的影响；模拟PGA最大的区域位于震中附近，最大值大约为850cm/s²。为了验证结果的有效性，我们对比了模拟PGA的分布图与中国地震局(2017)发布的本次地震烈度图(图 6(b))，可以看出，模拟PGA最大的区域与本次地震最大烈度区(Ⅸ度)的位置基本一致。除此之外，PGA和烈度的等值线都呈椭圆形，沿断层走向方向较长，垂直断层走向方向较短，说明本文模拟的PGA空间分布基本符合观测到的地震烈度分布情况。

图 6  模拟地震动的PGA空间分布(a)和地震烈度分布图(b)

Figure 6.  The spatial distribution of PGA of synthetic ground motions (a) and the intensity distribution map of the earthquake (b)

下载: 全尺寸图片 幻灯片

3.   结论

本文采用改进的随机有限断层法合成了2017年九寨沟M_S 7.0地震的近场地震动。结果表明：

(1)选取合适的震源模型和输入参数，通过“试错法”，按照误差最小原则估计本次地震的应力降为4.0MPa，与王宏伟等(2017)反演地震记录S波傅氏谱的得到结果3.854MPa基本一致。与EXSIM方法相比，Wang等(2015)的方法得到的结果在频域上与实际地震记录更为符合。

(2)按照美国NEHRP场地划分标准，依据V_S30确定了断层附近台站的场地类型。在模拟地震动过程中，通过场地放大因子近似考虑场地效应的影响。最后，得到了强震台站的模拟加速度时程和反应谱，以及研究区域的模拟PGA分布图。部分台站的模拟地震动和实际记录在峰值加速度和反应谱方面较为符合，模拟的PGA分布图反映出地震烈度图的基本特征，验证了本文结果的有效性。

(3)部分台站的模拟结果在幅值、持时以及反应谱方面与实际记录存在不同程度的差异，说明本文采用的随机有限断层法和模型参数仍有可改进之处。

(4)对九寨沟强震的模拟结果说明，在缺少近断层记录的情况下，综合考虑震源、传播路径和场地条件人工合成地震动，总体上能够有效地重现实际地震记录的峰值和反应谱特征，为该地区的地震危险性分析和抗震设计提供一定依据。