引言

对工程场地开展地震危险性分析，可以有效地评估与区域地震活动性趋势相关的地震作用影响。为了解地震活动的危险性并减小其所造成的危害，全球范围内已广泛实施地震危险性分析（Douglas，Edwards，2016）。地震危险性分析方法分为概率方法和确定性方法两大类。合理计算目标区域的震源参数、传播路径参数和局部场地放大系数等基本参数是实施确定性方法的关键环节，而且这些参数的确定也是工程地震学中建立强震动记录稀少地区地震动预测方程的基础（Campbell，2003；Atkinson，2008）。

龙门山地区是我国地震活动性最强的地区之一，且人口较为密集，2008年5月12日发生的汶川M_S8.0地震和2013年4月20日发生的M_S7.0芦山地震，均造成了严重的人员伤亡和经济损失。因此，研究该地区的地震动预测模型参数，对于该地区的强地震动预测和抗震设防具有重要的实际意义。

许多研究人员采用Atkinson和Mereu （1992）提出的迭代拟合法和Moya等（2000）提出的遗传算法计算了四川不同分区的震源参数、品质因子和测震台站的局部场地放大系数（乔慧珍等，2006；张永久等，2007；华卫等，2009；范军等，2012）。但是，由于测震台站均严格建在基岩场地上，因此并未得到工程关注的软弱覆盖层所引起的局部场地放大效应。汶川地震和芦山地震之后，中国国家强震动台网中心收集了一批宝贵的强震动记录，为研究四川地区地表覆盖土层引起的局部场地放大效应提供了数据基础。基于汶川地震余震的强震动记录，喻畑和李小军（2012）采用四分之一波长法计算了浅硬土层场地和深厚土层场地的平均场地放大系数，同时反演得到了13次中强余震的震源参数、四川山区和盆地地区的品质因子以及不同场地类型的高频衰减参数。Ren等（2013）采用广义线性反演方法（generalized linear inversion technique，简写为GIT）计算了龙门山地区28个强震动台站的局部场地放大系数，并分析了土层的非线性效应。除了汶川地震产出的数据，龙门山地区在过去几年间也积累了较为丰富的强震动记录，采用更多的数据开展该地区的震源谱、传播路径参数和局部场地效应的研究，可以对前人的研究成果进行补充，并可能得到新的认识。

鉴于此，本文拟收集2008—2015年期间汶川地震和芦山地震余震的强震动记录，在对数据进行筛选和预处理之后，通过GIT方法计算这些地震的震源参数、龙门山断裂带西部山区和东部盆地地区的品质因子以及83个强震动台站的局部场地放大系数，以期为龙门山地区及其邻区未来的强地震动模拟研究及相关工程实践提供一套可供参考的地震动预测模型参数。

1.   数据筛选和预处理

中国国家数字强震动观测台网于2008年3月正式运行，在龙门山及其邻区共建有147个强震动观测台站（Li et al，2008 ）。本文共收集到2008—2015年期间这些台站得到的约5 000组地震动三分量加速度记录，按照以下标准对其进行筛选：① 为了尽量减小强震的断层尺度效应、远场面波效应和土层非线性效应对反演结果的影响，选用震级范围为3.0≤M_S≤5.0、震源距小于200 km且峰值加速度小于100 cm/s²的加速度记录；② 考虑到点源近似假设的适用性以及震源深度的不确定性，选用加速度记录的震中距必须大于震源深度；③ 为了尽量避免加速度记录中信号噪声的影响，信噪比必须大于3；④ S波时窗长度选为S波到时至80%能量持时，因为太短的S波时窗可能会改变傅里叶振幅谱的大小和形状，因此规定可靠的持时长度必须大于3 s；⑤ 为了保证反演的稳定性，每个地震事件至少需要被3个台站记录到，每个台站至少需要记录到3个地震事件。

按照以上标准最终得到83个强震动台站记录到的137次汶川地震和芦山地震的余震（3.0≤M_S≤5.0）所产出的2 082组水平向S波时程。首先，使用汉宁（Hanning）窗对这些S波时程的两端2.5%持时进行尖灭处理以防止能量泄露；然后，采用标准快速傅里叶变换将S波时程变换到频域，再用Konno-Ohmachi滤波器对谱进行平滑（Konno，Ohmachi，1998）；最后，对每个记录的东西向和南北向水平分量求矢量和以尽可能消除方位角的影响，这样得到每条记录的平均谱为

式中，S_EW和S_NS分别为东西向和南北向水平分量。

图1为所选强震动台站和地震震中的位置以及每个地震到台站的传播路径。因为以龙门山断裂带为界，西部山区和东部盆地地区的地壳波速结构存在显著差异，导致两个地区的传播路径引起的衰减效应不同（胥颐等，2009；Fu，Li，2016）。所以，本文按照地理位置将台站分为山区台站和盆地台站，在反演过程中分别考虑两个地区的品质因子和几何扩散模型。由于震中位置和强震动台站选址的客观约束，本文选用的大部分台站均位于四川盆地西部边缘。对于这部分台站，我们将高程大于1 000 m的台站视为山区台站，小于1 000 m的台站视为盆地台站。基于《建筑抗震设计规范》（GB50011—2010）中的场地分类标准，本文所采用的83个强震动台站的场地类型均划分为II类场地（中华人民共和国住房和城乡建设部，中华人民共和国国家质量监督检验检疫总局，2016）。为了对这些台站的局部场地放大系数开展更细化的研究，本文根据美国下一代衰减模型（next generation attenuation，简写为NGA）场地分类标准将场地类型进行划分（Chiou et al，2008 ）。本文涉及的83个强震动台站的场地，包括3个A类，32个B类，26个C类和1个D类场地，另外21个台站没有场地资料，故无分类（图1）。

图  1  强震动台站和地震事件的地理位置

Figure  1.  Geographical location of selected strong motion stations and events

下载: 全尺寸图片 幻灯片

2.   反演方法

在频域内，S波的理论傅里叶振幅谱可以看作是相互独立的震源效应、传播路径效应与场地效应的乘积，即

式中：i，j和k分别表示地震、台站和频率的编号，R_ij表示第i次地震到第j个台站的距离；E_i（f_k），D（R_ij，f_k）和S_j（ f_k ）依次表示第i次地震、第j个台站和第k个频率的震源项、路径项和场地项。

对于满足远场近似的地震，震源项用ω⁻²点源模型表示为（Brune，1970）

式中：M_0i和f_0i分别为第i次地震的地震矩和震源拐角频率；C为常系数，即

式中：R_θφ为S波的平均辐射花样，本文取R_θφ=0.55；F为自由表面放大系数，本文取F=2；V为S波水平向能量转换系数，本文取V=0.707；ρ和β分别为震源附近的平均密度和S波速度，本文取ρ=2.8 g/cm³，β=3.5 m/s；R₀为参考距离，本文取R₀=1 km （Boore，2003）。

传播路径项包括几何扩散系数G（R_ij）和非弹性衰减函数，

This page contains the following errors:

Below is a rendering of the page up to the first error.

式中，当R_ij≤R₀₁时，台站接收到的地震波主要是直达波；当R₀₁<R_ij≤R₀₂时，地震波包含直达波和地壳内不连续面和莫霍面的反射波；当R_ij>R₀₂时，地震波包含更为复杂的多重反射波和折射波。本文采用经验的几何扩散模型，即b₁=1.0，b₂=0，b₃=0.5，R₀₁=1.5H和R₀₂=2.5H，其中H表示龙门山断裂带上、下盘的平均莫霍面深度（温瑞智等，2015）。基于龙门山地区的莫霍面深度模型，将山区和盆地的平均莫霍面深度分别取为55 km和45 km （胥颐等，2009）。

传统的GIT方法假设参考基岩场地的场地放大系数在全频段均为1，因此在使用时需要事先找到一个基岩场地作为参考场地；而实际中，这样的场地通常难以找到。即使是地表硬基岩场地，在风化作用的影响下其场地放大系数通常也不为1 （Hassani et al，2011 ）。我们采用一种替代的参考场地选择方式，即预先用四分之一波长法计算区域的平均场地放大系数，用谱衰减法计算每个台站的高频衰减参数κ₀，然后用二者的乘积作为经验的参考场地（喻畑，李小军，2012；傅磊，李小军，2017）。这样，场地项可表示为参考场地条件S_ref（f_k）与残留项S_res（f_k）的乘积，即

将式（5）代入式（2）后进行线性对数化处理，可得

其矩阵形式表示为

式中：A为包含震源项、品质因子和场地项等3个非零项的矩阵；x表示式（8）右边所有未知参数的矢量；b为数据矢量；w为场地条件约束，即假设所有台站的场地放大系数中各个频率的残留项对数和为0，即

最后，利用最小平方QR分解法（least-squares QR-factorization，简写为LSQR）求解式（8）的稀疏矩阵以得到震源谱、品质因子和局部场地放大系数（Paige，Saunders，1982a，b）。

3.   结果

3.1   误差分析

通过把反演得到的震源谱、品质因子与局部场地放大系数相乘，本文得到了所有强震动记录的反演傅里叶振幅谱（简称为反演谱），将其与每条记录的观测傅里叶振幅谱（简称为观测谱）进行比较以分析反演结果的误差。图2为反演结果的残差及其标准差，残差定义为观测谱与反演谱的对数差。在数据处理过程中，每条水平向平均傅里叶振幅谱在0.5—30 Hz频段内取50个对数等间隔采样点，其中盆地地区和山区分别包含16 350和35 700个数据点。可以看到：盆地和山区的残差均符合正态分布，并且随着频率、震源距和震级的增大均无明显的偏差；盆地和山区的标准差分别为0.23和0.30。以上结果表明，所采用的模型与数据符合良好。

图  2  盆地地区（a）和山区（b）的残差及其标准差σ

Figure  2.  The residual and standard derivation σ of basin area （a） and mountainous area （b）

下载: 全尺寸图片 幻灯片

3.2   震源参数

图3为各个震级的平均加速度震源谱，拟合模型采用ω⁻²模型。在较小的震级范围（M_S≤3.9）内，反演谱振幅在小于1 Hz左右的低频段大于拟合谱，表明这一频段内的震源谱衰减略小于ω⁻²。除M_S4.9和M_S5.0外，其余震级的高频段拟合结果均较好。M_S4.9和M_S5.0的部分反演谱中，可以观察到明显的高频衰减现象。注意到反演过程中，参考场地已经考虑了每个台站由场地的κ₀引起的高频衰减现象，因此M_S4.9和M_S5.0的反演谱中所呈现的高频衰减现象可能是由震源效应引起，并且可能与震级的大小有关。

图  3  3.0<M_S≤5.0范围内各震级的加速度震源谱

Figure  3.  Acceleration source spectra of the events with M_S3.0 to 5.0

下载: 全尺寸图片 幻灯片

由式（3）可知，控制震源谱的参数主要为应力降Δσ，M₀和f₀，三者之间的关系为（Boore，2003）

因此，接下来本文将围绕这3个参数来展开研究。尽管对于芦山地震是否为汶川地震的余震至今尚存在争议，但是因为汶川地震和芦山地震均发生在龙门山构造带上，本文假设汶川和芦山地震的余震具有可比性。图4为本文所采用的地震事件的应力降随时间的变化趋势，可以看出：从2008年5月12日（汶川地震主震发生）至2013年4月20日（芦山地震主震发生），余震事件的应力降随着时间逐渐增大；在芦山地震的主震发生之后，应力降明显下降，并在一段时间后开始回升。以芦山地震的主震发震时刻为分界，其前后两个时期的应力降变化范围分别为0.32—13.3 MPa和0.36—27.4 MPa。据此，本文认为对于图4中所反映的应力降变化趋势可以作如下定性解释：在汶川和芦山主震发生之后，震源区累积的大部分应力已经得到释放；随后区域的应力重新加载，震级较小的余震对区域的应力恢复影响较小，因此随着时间的增加，余震的应力降表现出增大的趋势，直到下一次较大震级的地震发生，导致该区域的应力场再次产生较大的改变。这一解释与余震的应力降通常小于主震事件的事实相符。值得注意的是，因为在数据筛选过程中所保留余震事件的时间分布并不均匀（2009—2011年的数据相对较少），因此本文的结果无法准确地表示汶川地震序列和芦山地震序列的应力降恢复周期。

图  4  汶川地震和芦山地震余震的应力降Δσ随时间的变化趋势

Figure  4.  Temporal evolution of stress drop Δσ in Wenchuan and Lushan aftershocks

下载: 全尺寸图片 幻灯片

将M_S4.7—4.9的观测震源谱分为两组，第一组为汶川地震和芦山地震主震发生后1个月内的余震数据，第二组为其它时段的余震数据，分别求出两组数据中各震级观测谱的平均谱，如图5所示。可见，在M_S4.7和M_S4.9中，第一组数据的平均震源谱没有观察到高频衰减现象，而第二组数据存在明显的高频衰减现象。因此，由震源引起的高频衰减现象还可能与能够描述区域内震源时间相关性的参数有关。

图  5  两组数据中M_S4.7，4.8和4.9的平均震源谱

第一组为汶川地震和芦山地震主震发生后1个月之内的数据，第二组为其它时段内的数据

Figure  5.  Average source spectrum of M_S4.7，4.8 and 4.9 obtained from two datasets

Group 1 includes events occurred during the first month after the origin of Wenchuan and Lushan mainshocks，group 2 includes events occurred during other periods

下载: 全尺寸图片 幻灯片

图6给出了两组结果的M₀-f₀相关性比较，可见第二组的平均应力降为3.0 MPa，大于第一组的1.61 MPa。任叶飞（2014）利用2008年5月12日至2008年7月16日期间的汶川地震余震记录反演得到的平均应力降为1.24 MPa，第一组结果与其相似。Beresnev （2001）指出，Brune应力降并不是地震发生前后震源构造区的真实应力差，而是一个控制傅里叶谱高频能量大小的应力参数，该参数反映的是破裂传播的最大速率，即

式中，c为常系数，v_S和ρ分别为震源区附近的剪切波速和密度。综合以上结果本文推测震源产生的高频衰减现象可能与震级有关，并且与破裂传播的最大速率有关；当震源的破裂传播速率足够大时，震源谱中可能会出现明显的高频衰减现象。但是这一推论尚需开展更深入的研究加以论证。

图  6  第一组（a）和第二组（b）的M₀-f₀相关性比较，灰色区域为本文结果±1倍标准差

Figure  6.  Comparison of M₀-f₀ correlation for the data of group 1 （a） and 2 （b），where gray areas represent the results of this study ±1 standard derivation

下载: 全尺寸图片 幻灯片

图7为本文得到的M_W-M_S相关性，即

可以看到，M_W整体上略大于M_S，且随着震级的增大两者逐渐接近。

图  7  本文选用地震事件的M_W-M_S的相关性

Figure  7.  M_W-M_S correlation of the events used in this study

下载: 全尺寸图片 幻灯片

3.3   品质因子Q

图8为反演得到的山区和盆地的品质因子，分别为Q（ f ）=264 f ^0.75和Q（ f ）=223 f ^1.01。华卫等（2009）采用汶川余震序列的测震记录反演得到山区和盆地的品质因子分别为Q（ f ） =274.6f ^0.423和Q（ f ）=206.7 f ^0.836，张永久等（2007）采用测震记录计算得到盆地的品质因子为Q（ f ）=217.8 f ^0.816，本文结果与其基本一致，但是本文获得的山区和盆地的频率相关性均较上述文献中的结果偏大。这可能是因为研究区的强震动台站大部分位于山区与盆地之间的构造边缘，距离地壳活动强烈的龙门山断裂带较近，这些地区的区域构造非均匀性通常较强，可能造成品质因子的频率相关性偏高（Roecker et al，1982 ）；另一方面，任叶飞（2014）采用汶川余震序列的强震动记录得到汶川周边地区的品质因子为Q（f ）=100.6 f ^1.1，温瑞智等（2015）采用芦山余震序列的强震动记录得到芦山周边地区的品质因子为Q（ f ）=31.9 f ^1.04，本文结果小于二者的Q₀，但是本文盆地地区的结果与其频率相关性结果较为相似。

图  8  山区和盆地地区的品质因子Q

Figure  8.  Quality factor Q obtained respectively for mountainous and basin regions

下载: 全尺寸图片 幻灯片

3.4   局部场地效应

本文计算83个强震动台站的局部场地放大系数，并与谱比法（horizontal-to-vertical spectral ratio method，简写为HVSR）所得结果进行比较。HVSR计算中，排除峰值加速度大于100 m/s²的记录以尽量避免局部土层非线性效应的影响。图9展示了随机选择的20个台站的场地放大系数。可以看到，GIT和HVSR两种方法所得到的局部场地放大系数的形状较为相似，二者的结果均能较好地体现场地的卓越频率。在小于场地卓越频率的频段内，GIT与HVSR更为相似；在大于场地卓越频率的频段内，GIT的结果则普遍大于HVSR。Parolai等（2004）采用土耳其31个强震动台站的记录比较了HVSR与GIT的差异，本文结果与其结论相似。

图  9  20个强震动台站的GIT与HVSR得到的结果比较

Figure  9.  Local site amplification coefficients of 20 strong motion station sites obtained from GIT and HVSR

下载: 全尺寸图片 幻灯片

图10为A，B和C这3类场地的平均场地放大系数，相应的场地卓越频率分别为7.87，7.24和4.04 Hz，符合场地越硬、卓越频率越大的一般规律。部分B类和C类场地的局部场地放大系数大于A类场地，但是3类场地的平均最大放大系数较为接近。随着场地趋于软弱，平均放大系数的峰值平台区明显向低频部分发育，表明软弱场地对较长周期地震动有更明显的放大效应。因为根据我国的场地分类标准，本文的所有台站均划分为II类场地，从放大系数的峰值大小而言这种分类是合理的，但是如果考虑到放大系数的峰值平台区宽度，NGA场地分类标准则相对更为细致。需要注意的是：A类场地只有3个，结果可能不具有代表性；D类场地只有62WUD，其场地卓越频率为1.48 Hz，不作统计。

图  10  基于NGA场地分类标准的场地放大系数

（a） A类场地；（b） B类场地；（c） C类场地

Figure  10.  Average site amplification coefficients based on the NGA site classification

（a） Site class A；（b） Site class B；（c） Site class C

下载: 全尺寸图片 幻灯片

4.   讨论与结论

本文采用2008—2015年中国强震动台网中心收集的龙门山地区的强震动记录，筛选出83个强震动台站记录到的2 082组水平向强震动记录，研究了龙门山地区的震源参数、品质因子和局部场地放大系数，得到以下主要结论：

1） 汶川地震和芦山地震的余震序列的应力降随着时间的推移均有逐渐增大的趋势。以芦山地震的主震发震时刻为界，其前后两个时段的余震序列应力降变化范围分别为0.32—13.3 MPa和0.36—27.4 MPa。这一现象定性地反映了本文选用地震发生时段内区域应力降随时间的恢复过程。

2） 部分M_S≥4.7地震事件的观测震源谱中呈现明显的高频衰减现象。因为在反演过程中已经考虑了各台站由局部场地条件引起的高频衰减系数，因此认为观测震源谱中所呈现的高频衰减现象是由震源引起。目前，对于震源效应引起的高频衰减效应的研究较少，Papageorgiou和Aki （1983）提出震源产生的高频衰减与破裂断层前端的软熔带宽度以及震源破裂传播速度等震源特性有关，Wen和Chen （2012）研究了汶川地震中高频衰减频率的特性及其与震源之间的物理关系。根据本文结果推断，这一高频衰减现象可能在震源破裂最大滑动速率达到一定程度时出现，这在一定程度上支持了Papageorgiou和Aki （1983）提出的震源产生高频衰减效应的机理。但是对于这一推论，尚需进一步研究加以讨论。

3） 计算得到山区和盆地地区的品质因子分别为Q（ f ）=264 f ^0.75和Q（ f ）=223 f ^1.01。其中，Q₀与张永久等（2007）和华卫等（2009）基于测震记录的研究结果相似，频率相关性偏大可能是因为龙门山地区的强震动台站主要分布在龙门山断裂带两侧较近的区域，受到区域构造较强的非均匀性影响所致。因此，尽管本文得到的山区和盆地地区的品质因子与前人结果较为相似，但作者建议这一结果只适用于本文台站与地震事件连线所覆盖的区域。

4） 分别采用GIT和HVSR方法计算的20个强震动台站的场地放大系数的比较结果表明：两种方法得到的局部场地放大系数均能较好地反映出场地卓越频率，但是总体上GIT的振幅大于HVSR；在小于场地卓越频率的频段内，GIT与HVSR的结果较为接近；在大于场地卓越频率的频段内，GIT的结果大于HVSR。本文结果显示，HVSR可以较好地反映小于场地卓越频率的低频部分的局部场地放大效应，但这一结果尚需更多的研究加以证实。

5） 根据NGA场地分类标准，将有钻孔资料的62个强震动台站分为3个A类，32个B类，26个C类和1个D类场地，并分别计算了A类，B类和C类场地的平均场地放大系数，其场地卓越频率依次为7.87，7.24和4.04 Hz，符合场地越坚硬、卓越频率越大的一般规律。并且，随着场地变得软弱，场地放大系数的峰值平台区明显向低频部分发育。根据我国场地分类标准，本文所用台站均为II类场地。本文采用NGA场地分类标准的结果显示，II类场地可以较好地反映局部场地放大系数的峰值大小，但是难以体现峰值平台区宽度的影响。因此本文认为，对我国场地分类标准开展进一步的细化工作，对于指导建筑抗震设计具有积极的意义。

感谢中国地震局工程力学研究所“国家强震动台网中心”为本研究提供数据支持，感谢两位审稿专家对文章的不足之处提出中肯且具有建设性的宝贵建议。