强地震动衰减是工程地震学的一个重要研究方向.在地震危险性评估中，不管是概率方法还是确定性方法，都需要借助地震动衰减关系估算场地的地震动期望值.大多数衰减关系是经验性的，研究基础是强地震动的观测数据.在美国西部、日本等观测数据丰富的地区，相关研究一直处于前列^[1-3].近几年，美国新一代地震动衰减关系项目 (NGA) 取得的进展影响很大^[4]，最大的贡献是发布了一个当前最丰富的地震动数据库，不仅包括美国和加拿大西部的数据，也包括中国台湾、土耳其、中国大陆、伊朗、意大利、墨西哥、新西兰以及日本等国家和地区的强震记录^[5].此举曾使人产生错觉，以为在全球数据的基础上，有可能建立全球适用的衰减关系.虽然NGA考虑得更加全面，参数增加很多，成果也鼓励世界上其他地区试用，但项目的成果以NGA-West命名，说明这类关系肯定是区域性的，美国中东部仍然需要另外研究衰减关系，NGA-East就是最有力的证明.从道理上说，区域的地壳结构、地壳介质动力性质及震源机制的不同都会影响地震动衰减规律，世界上大多数国家或地区尚缺少足够的强地震动观测数据，难以直接建立经验性的衰减关系.在中国，对观测数据丰富地区统计结果进一步修正的作法得到广泛应用^[6-7]，亦受到一些质疑^[8-9].与此同时，其他国家的地震学者发展了基于随机振动理论的简化方法^{[1, 10-12]}.十多年来，中国研究人员陆续有相关论文发表^[13-16].本文是在这些相关研究进展基础上的一个总结推进.

1 从基于随机振动理论到基于地震学的强地震动衰减关系

基于随机振动理论的方法，用震源谱表达震源对地震动的影响，用几何衰减和非弹性衰减项表达传播途径对地震动的影响，用场地放大和高频截止项表达局部场地条件对地震动的影响，比经验性衰减关系中只用震级一个参数、一个综合衰减项和一个场地分类项表达这三方面的影响，更全面，更突出强调物理意义^{[9, 15, 17]}.

点源引起地表一点的地震动傅里叶幅值谱可表达为^{[2, 18]}

$ \begin{array}{l} {F_{\rm{A}}}\left( {{M_0},f,R} \right) = C \cdot S\left( {{M_0},f} \right) \cdot G\left( R \right) \cdot D\left( {R,f} \right) \cdot \\ \;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;A\left( f \right) \cdot P\left( f \right) \cdot I\left( f \right). \end{array} $

(1)

式中：M₀为地震矩，f为频率，R为震源距，C为比例系数，可参照震源附近地壳介质密度和剪切波速取值，S(M₀, f) 为震源谱，描述从震源辐射的总能量大小和能量在频域的分布，G(R) 为几何衰减项，描述几何扩散导致的地震动衰减，D(R, f) 为非弹性衰减项，描述地震波传播过程中地壳介质耗散能量引起的地震动衰减，A(f) 为近地表幅值放大因子，描述近地表地壳的速度梯度变化对地震动的影响，P(f) 为近地表高频截止项，描述近地表地震动中不含有很高频率成分的现象，I(f)=(2πf)^z中分别取z=0、1或2，表达地震动位移和速度、加速度的关系.

应用最广泛的震源谱是Brune^[19-20]提出的ω²模型.为避免合成大地震的地震动中子源尺寸对地震动的影响，作者课题组借鉴Masuda^[21]的谱形提出如下改进^{[13, 22]}：

$ S\left( {{M_0},f} \right) = \frac{{{M_0}}}{{{{\left[ {1 + {{\left( {\frac{f}{{{f_0}}}} \right)}^a}} \right]}^b}}}. $

(2)

式中：f为频率，f₀为拐角频率，a、b是同震级有关的系数，a=3.05-0.33M_w，b=2.0/a.

Brune^[19-20]推导了圆盘破裂模式的拐角频率：

$ {f_0} = 4.9 \times {10^6}{\beta _{\rm{s}}}{\left( {\Delta \sigma /{M_0}} \right)^{1/3}}. $

(3)

式中Δσ为应力降，描述地壳应力在地震后相比地震前的降低，不能直接测量，可借助远场记录反演估计.一般认为，Δσ与地震破裂面上的错动成正比，同破裂面的大小成反比.大多数研究者认为大地震的应力降是基本稳定的，对小地震的应力降尚有不同的观点.Shearer^[23]研究论证应力降呈区域性.

几何衰减项G(R) 可以分不同的距离段定义，分别反映近距离处地震动成分以体波 (剪切波) 为主、远距离处以面波为主以及两者之间的一段距离上直达波与来自莫霍面临界反 (折) 射波叠加形成的地震动衰减缓慢甚至基本不衰减 (称为地壳波导效应) 的不同特征.几何衰减项与区域地壳速度结构有密切关系，通常用三段函数表示^[12].

$ G\left( R \right) = \left\{ \begin{array}{l} \frac{1}{R},\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;R \le {R_1};\\ \frac{1}{{{R_1}}},\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;{R_1} \le R \le {R_2};\\ \frac{1}{{{R_1}}}\sqrt {\frac{{{R_2}}}{R}} ,\;\;\;\;\;\;\;\;\;\;\;{R_2} \le R. \end{array} \right. $

(4)

式中R₁和R₂为距离分段点.

非弹性衰减项D(R, f) 可表达为^[18]

$ D\left( {R,f} \right) = {{\rm{e}}^{ - \frac{{\pi fR}}{{Q\beta }}}}, $

(5)

式中Q为品质因子，其值越大，地震波的衰减越慢，可表达为

$ Q = {Q_0}{f^\eta }, $

(6)

式中Q₀为f=1 Hz的品质因子，η为常数.

近地表幅值放大因子A(f) 可以根据区域地壳速度结构按四分之一波长法近似确定^[24].一般认为，近地表高频截止项P(f) 描述的是传播路径和场地效应，或是这两种效应与震源效应的综合.普遍使用的有f_max滤波器和kappa滤波器，分别表达为^[25-26]：

$ P\left( f \right) = {\left[ {1 + {{\left( {\frac{f}{{{f_{\max }}}}} \right)}^8}} \right]^{ - 1/2}}, $

(7)

$ P\left( f \right) = {{\rm{e}}^{ - \pi f\kappa }}. $

(8)

式中f_max为地震动高频截止频率，κ称为“kappa”系数，一般认为与距离和场地条件密切相关.

大地震近场地震动可以借助有限断层震源模型合成，此处不详述.

由上述可知，地震动傅里叶幅值谱的估计中，震源谱中的应力降Δσ、品质因子Q中的Q₀和η以及几何衰减项的两个分段点距离R₁和R₂等5个区域性参数，分别控制了近场地震动的强度和随距离衰减的速率，通常多不能直接量测到.不同研究成果采用的基本数据，甚至公式的形式往往不同，结果差别不小，直接采用随机性过大，对地震学研究基础相对薄弱的地区更不可靠.为克服这一瓶颈，与地震研究者对其中某些参数反演研究的思路相类似^[27-29]，本文作者提出、发展了根据宽频带数字地震监测台网的小震记录，借助遗传算法对这5个参数联合反演估值的方法^[30-31].

遗传算法是一种启发式的、有指导的蒙特卡洛反演方法，全局搜索能力强，不需要计算目标函数的导数，可以解决复杂的大尺度、多变量非线性反演问题.作者采用的微遗传算法是传统遗传算法的一种改进，便于在较小的种群中得到较好的优化结果，应用很广泛.先从一组随机数中进行初始化操作产生个体，计算每个个体的适应度；再通过锦标赛选择、精英策略和交叉操作产生新个体；通过判断种群收敛性确定是否重复以上步骤；在所有的代数中，寻求适应度最大的数值组合个体作为最终结果^[32].

反演中，一般选用震级在3.5~4.5之间的数十次区域浅源构造地震中几十个数字测震台站记录到的上百或几百条速度时程.以相应傅里叶幅值谱或其包络线与按式 (1) 中I(f) 取2πf、根据震级推算的地震矩及选定的一组5个区域参数计算得到的速度傅里叶幅值谱贴近为目标，构造目标函数：

$ {\varphi _j} = \sum\limits_{k = 1}^N {\sum\limits_{i = 1}^{4\;096} {{{\left[ {{F_{{A_0}}}\left( {i,k} \right) - {F_{{A_j}}}\left( {i,k} \right)} \right]}^2}} } . $

(9)

式中：N为小震记录的总数，F_A0(i, k) 为第k个记录的傅里叶幅值谱 (离散成4 096个点) 的第i个幅值，F_Aj(i, k) 为据第k个记录对应的震级、距离和第j代5个区域参数值，按式 (1) 计算得到的傅里叶幅值谱的第i个幅值.采用式 (10)，可以将适应度压缩至 (0, 1].

$ {F_j} = \exp \left( { - \beta {\varphi _j}} \right), $

(10)

式中β为适应度变换系数.进一步设计反演策略，就可以在给定的相当宽泛的取值域中搜索寻优，得到合适的参数值.如此，既扩展了丰富的数据源，又与地震学密切结合，体现了区域特征.

给定一组M₀和R，将上述区域参数分别代入式 (3)~(5)，借助式 (1) 估计得到地震动的加速度傅里叶幅值谱后，可以据其与功率谱的关系，以及功率谱零阶谱矩m₀的定义，按式 (11) 计算其值.

$ {m_0} = 2\int\limits_0^\infty {{{\left| {{F_A}\left( f \right)} \right|}^2}{\rm{d}}f} . $

(11)

由Parseval定理，加速度均方根值A_rms可表示为

$ {A_{{\rm{rms}}}} = {\left( {\frac{{{m_0}}}{{{T_{\rm{d}}}}}} \right)^{1/2}}, $

(12)

式中T_d为强烈震动段的持续时间.综合考虑震源辐射模式和台站方位效应、自由地表放大效应及地震动能量的水平分量后，还要再乘上一个系数0.85^{[3, 33]}.

均方根值A_rms乘以峰值因子γ_m^[34]，可得到相应加速度峰值PGA.

$ {\gamma _{\rm{m}}} = 2\int\limits_0^\infty {\left\{ {1 - {{\left[ {1 - \xi \exp \left( { - {z^2}} \right)} \right]}^{{N_{\rm{e}}}}}} \right\}{\rm{d}}z} . $

(13)

式中：ξ为带宽因子，描述随机变量的分布形式，$\xi = \frac{{{N_{\rm{z}}}}}{{{N_{\rm{e}}}}} = \frac{{{m_2}}}{{{{\left( {{m_0}{m_4}} \right)}^{1/2}}}}, {N_{\rm{z}}}$和N_e分别为T_d内零交次数的和达到极值次数的期望，N_z=2f _zT_d，N_e=2f_eT_d，f_z和f_e分别为零交和极值出现的频率，${f_{\rm{z}}} = \frac{1}{{2\pi }}{\left( {{m_2}/{m_0}} \right)^{1/2, }}{f_{\rm{e}}} = \frac{1}{{2\pi }}{\left( {{m_4}/{m_2}} \right)^{1/2}}$，其中m₂和m₄分别为二阶和四阶谱矩.

按上述方法得到地震动衰减关系是基于随机振动理论.推导中，式 (13)(称为峰值因子) 的推导，均是以平稳随机过程理论为基础的，用于强烈非平稳的地震动时程，即便是只考虑相对平稳的强烈震动段亦颇受质疑.作者课题组随后的研究表明，采用这种简化会导致误差增加10%~15%^[35]，没有必要为节约不很多的计算时间而付出准确性的代价.可以利用地震动的傅里叶幅值谱和一个随机相位谱结合直接生成时程，从中读取峰值，取不同相位谱多次合成时程的峰值平均值作为期望值，故改称为基于地震学的地震动衰减关系^[36].

2 建立基于地震学的强地震动峰值衰减关系的3个实例

上述思路先在数据丰富的地区研究验证，是方法发展与推广应用不可缺少的一个环节，本文介绍的第一个实例就选在日本东北地区.随后，在中国积累有一些强地震动记录的川滇地区试用，进一步与区域强地震动观测的数据比较，即便是不够充分的数据，亦可在一定程度上说明本文方法的适用性.

2.1 例1

区域范围取东经138°~143°、北纬36°~40°.小震数据取自宽频带数字地震台网F-net (www.fnet.bosai.go.jp) 的12个台站，1996年1月至2010年10月间，震源深度不大于30 km的632次小地震 (M_w=3.5~4.5) 的850条记录.检验用的强地震动数据取自K-NET (www.kyoshin.bosai.go.jp) 的89个基岩台站，在同一时间段187次地震 (M_w≥4.5) 的4 323条观测记录.

参考区域地震学研究成果，选定Δσ、Q_o、η、R₁、R₂等5个参数的可能取值范围，用微遗传算法反演得到满意的解，列于表 1^[35].将反演结果代入式 (3)~(5)，对一系列M和R，按式 (1) 计算傅里叶幅值谱，再进一步借助式 (11) 和式 (12) 计算地震动均方根加速度值，与峰值因子相乘，得到地震动峰值加速度衰减的结果.

为了比较、检验，图 1中分别表示了5、6、7级地震的PGA衰减的结果 (实心点) 和相应震级±0.5范围内的强地震动观测数据 (小圈)，以及日本学者发表的3组用强地震动观测数据统计回归建立的区域经验性衰减关系^[37-39].

从图 1可见，本文完全根据小地震的速度记录反演区域参数，一条强地震动数据都没有用，得到的结果与强地震动观测数据符合得很好，与日本学者统计得到的经验性衰减关系大体相同.

2.2 例2

区域范围选择四川省.小震数据取自区域地震观测台网29个数字台站，2001年初至2007年底之间，震源深度不大于30 km的82次小地震 (M_w=3.5~4.5) 的147条记录.检验用的强地震动数据取自中国地震局工程力学研究所“国家强震动台网中心”，选择同一时空范围内强震台网69个台站记录到的118次地震 (M_w≥4.5) 的1237条数据.参考区域地震学研究成果，选取的参数反演搜索范围和得到的区域参数值，列于表 2^[40].

据此建立的衰减关系与区域有限的强地震动数据 (主要是汶川地震余震的记录)，见图 2.图中分别表示了5、6、7级地震的PGA结果和相应震级±0.5范围内的强地震动观测数据，以及5组区域通常使用的衰减关系曲线^[41-45].为便于比较，凡分别长、短轴方向给出衰减关系的，距离均折算为等效圆半径.

从图 2可见，对应5、6级地震的衰减曲线从数据点丛中间穿过，上下两侧数据点的分布比较均匀，“盲测”成功.对于7级地震，本文方法的结果显得衰减慢些，在100 km以上的距离处高于观测值，是观测数据不足，还是其他问题有待进一步研究.

2.3 例3

区域范围取云南省.小震数据取自区域地震观测台网26个数字台站，2001年初至2007年底之间，震源深度不大于30 km的154次小地震 (M_w=3.5~4.5) 的863条记录.检验用的强地震动数据取自国家强震动台网36个台站，同一时空范围内27次地震 (M_w≥4.5) 的78条记录.参数反演搜索范围和得到的区域参数值，列于表 3^[40].

采用同样方法建立区域衰减关系与区域有限的强地震动数据，见图 3.分别表示了5、6、7级地震的PGA结果和相应震级±0.5范围内的强地震动观测数据，以及6组区域通常使用的衰减关系曲线^{[41-42, 44-47]}.

从图 3可见，结果与观测数据亦很匹配，观测数据点较为分散，曲线从中穿过.与区域通常使用的衰减关系相比，在100 km以内没有明显差别，在更远处本文的结果显得衰减偏慢，是否合理还有待观测数据积累来检验.

3 结语

世界上大多数国家或地区都缺少足够的强地震动观测数据，难以直接统计建立经验性地震动衰减关系.本文在作者十几年研究基于地震学的强地震动衰减关系的基础上，讨论了从基于随机振动理论的简化方法向基于地震学的地震动衰减的发展，详细介绍了借助地震监测台网的小震数据反演区域震源谱和衰减参数的技术途径.通过3个研究实例说明据此估计的地震动与区域 (有限的) 强地震动观测的数据比较，说明了方法的适用性.对于5、6级地震，两者相符；对于7级地震，在强地震动数据丰富的日本东北地区符合得很好，对强地震动数据有限的中国川滇地区在距离大于100 km的距离处符合得不够好，是观测数据不足，还是其他方面的问题，有待进一步研究.