地震动非参数化谱反演可靠性分析

王宏伟 ，任叶飞 ，温瑞智 ，周 影

（1.中国地震局工程力学研究所, 黑龙江 哈尔滨 150080；2.中国地震局地震工程与工程振动重点实验室, 黑龙江哈尔滨 150080）

地震动谱反演是频域上同时分离地震动震源、传播路径和场地效应的高效方法，该方法普遍应用于地震震源特征[1-4]、地壳介质衰减特性[5-6]、场地效应[7-8]、地震动不确定性[9]等的研究，这些研究对深入理解地震震源物理、探究地震发震机理、可靠预测地震动及地震危险性有着重要贡献.

根据传播路径衰减表示方式的差异，谱反演方法可分为参数化和非参数化两种.Castro等[10]首次提出的非参数化方法将传播路径衰减约束为距离的单调递减平滑曲线，既体现了传播路径衰减的复杂性，又描述了地震波在传播中的不断衰减，解决了参数化方法利用经验衰减模型会导致品质因子反演结果为极大值甚至负值，违背地震波衰减损耗的物理事实.之后非参数化方法开始广泛应用于地震活跃地区的地震波衰减特征及地震震源特征等的研究，例如日本[2]、新西兰[11]、欧洲[3,9,12]、中国四川地区[13]等，相比参数化方法，近年来非参数化方法的应用更普遍.

尽管非参数化谱反演的应用十分普遍，但反演结果的可靠性需要谨慎处理诸多关键环节来确保.谱反演必须考虑震源与场地项之间的权衡，一般通过选择合适的参考场地来实现，参考场地的选取显然会影响反演结果的可靠性，目前尚无研究针对如何选择合理的参考场地给出合理建议.根据反演步骤的差异，非参数化谱反演具体分为两步法[10]和单步法[2]，两种方法的反演结果是否存在明显差异并不清楚，也未有研究评价应优先选用哪种方法.根据理论震源模型[14]，震源谱以地震矩和拐角频率表示，因此震源谱的反演结果常用于地震震源参数的估计及震源物理特性研究，地震矩和拐角频率之间的权衡可能会影响震源参数估计的可靠性[15]，评估震源参数之间的权衡对其估计的影响十分必要.

针对影响地震动非参数化谱反演可靠性的关键环节，包括参考场地的确定、单步法或两步法的选择、震源参数估计，本文以2016—2017年意大利中部地震序列的地震动非参数化谱反演作为典型实例，欧洲强震动数据中心（engineering strong-motion，ESM）发布了地震序列八千余组人工处理的高质量的地震动记录，本研究限定记录峰值地面加速度（PGA）上限（100 cm/s2）、震源距上限（100 km）及台站（地震记录数下限）10 （10），选用了142个台站收集的78次震级（3.5～6.1）级地震的5 135组观测记录用于谱反演研究，截取记录S波并在其前后增加余弦边瓣窗口消除截断误差，S波到达时刻和结束时刻分别通过人工识别和能量法[3]确定，最后以两水平方向S波傅里叶谱均方根代表观测记录水平向S波傅里叶谱.本研究选用大量强震动记录可确保足够的数据冗余，避免因记录数过少造成反演结果的不可靠，通过系统分析反演结果，提出确保反演可靠性的上述关键环节的处理建议，可靠反演结果有助于地震震源物理的准确认识及地震动和地震危险性的有效预测.

1 非参数化谱反演方法

第i个地震（震级Mi）中第j个台站（震源距Rij）观测到的地震动水平向S波的傅里叶幅值谱（Oij）可表示为

式中：f为频率；Si(f,Mi)为第i个地震的震源谱；Aij(f,Rij)为传播路径衰减；Gj(f)为第j个台站的场地效应.

两步法首先从S波傅里叶谱中分离传播路径衰减项（第1步），即

式中：Mai(f)为与第i次地震震级大小相关的量.

之后从传播路径衰减修正的S波傅里叶谱中同时分离震源和场地项（第2步），即

单步法直接求解式（1），从S波傅里叶谱中同时分离震源、传播路径衰减和场地项，求解时需约束传播路径衰减及震源和场地项的权衡.

2 参考场地选择

参考场地一般为场地放大恒为1的无风化出露坚硬岩石场地[16]或平均场地放大为1的多个岩石场地[10]或至少一个放大效应已知的场地[17].实际应用中找到理想的出露坚硬岩石场地作为参考场地十分困难，多数研究采用场地放大效应恒为1的一个或多个岩石类场地作为参考场地[3,9,11].

一般考虑下列条件选取参考岩石场地：1）场地地质资料，例如：场地类别、场地等效剪切波速（例如地下30 m土层的等效剪切波速Vs30）等；2）水平/垂直（H/V）谱比曲线，岩石场地的谱比曲线在整个频段上平坦且幅值较小.为保证反演结果的稳定性和可靠性，Bindi等[17]规定尽可能选取观测记录足够多且记录到震级和方位角范围足够大的岩石台站.Pacor等[3]则假定所有台站均为岩石场地且平均场地放大恒为1，利用谱反演给出所有台站的场地放大效应，再选取整个频段上场地放大小于2的场地作为参考岩石场地.Oth等[11]首先对参考岩石场地进行初选，通过震源谱和场地效应的谱反演试算结果从初选中确定.鉴于理想岩石场地较难获取，上述这些研究提供了选取参考岩石场地的一系列方法，但这些方法均是不得已而为之，参考场地选取直接影响谱反演的可靠性.

根据Eurocode 8规范，本研究的142个观测台站中A类场地（岩石场地或其他接近岩石地质构造的场地，场地表面软弱层不大于 5 m 且Vs30＞800 m/s）台站共33个，其中7个为流动观测台站，A类场地台站作为参考场地的备选，因流动观测台站建设要求较低，受环境影响较大，7个流动观测台站不纳入备选之列.26个备选台站中，分别取记录数≥ 50组、[40～50）组、[30～40）组、[20～30）组、＜ 20 组的台站分别为7、5、2、4和8个，有详细钻孔数据的3个（LSS、AQP 和 PSC，Vs30分别为 1 091、836、1 000 m/s），只有LSS台站的观测记录充足（66组），AQP和PSC台站观测记录较少（分别为15组和14组）.

进一步缩小备选范围，仅以记录数不少于50组（73、111、155、105、60、58、157）的 7个台站（LSS、TOD、FIAM、GNU、MNF、SPM和 MGAB）作为备选，若最终证明均不合适，可适当放宽对备选台站记录数的要求.分别利用 7 个台站的 PGA≤ 50 cm/s2的观测记录计算0.25～20.00 Hz的H/V谱比曲线（见图1）.FIAM台站缺乏地形地貌描述且其H/V谱比曲线在1 Hz附近出现较为明显的峰值；SPM台站位于陡坡边缘，H/V谱比曲线幅值大且6.00 Hz附近出现明显峰值；GNU和MGAB台站的H/V谱比曲线幅值较大，说明这两个台站存在明显场地放大效应.上述4个台站均不适合作为参考岩石场地.LSS台站的观测记录方位角范围广（大约100°）、记录震源距范围较大（＞30 km）、H/V谱比曲线平坦且幅值小；TOD和MNF台站的H/V谱比曲线平坦且幅值小，此3个台站可作为初选的参考岩石场地.

图1 7 个台站的 H/V 谱比曲线Fig.1 The horizontal-to-vertical spectral ratios at seven stations

本研究以H/V谱比作为3个初选参考场地的场地放大效应，反演给出2016—2017年意大利中部地震序列观测记录的传播路径衰减、震源谱和场地效应，以LSS、TOD和MNF分别作为参考场地得到的6次典型地震的加速度震源谱如图2所示，6次地震（1#～6#）为 20 160 824 013 632_ML6.0（2016 年8月 24日 01：36：32发生的ML6.0级地震，余同）、20 160 825 031 716_ML4.5、 20 160 824 232 205_ML3.8、20 161 026 191 806_ML5.9、 20 161 031 070 545_ML4.1、20 161 103 003 501_ML4.8，震源谱大体上符合ω-2理论震源模型[14]，但不同参考场地得到的震源谱的形状及幅值差异明显.LSS和MNF台站作为参考场地的震源谱形状相似，LSS台站作为参考场地的震源谱的幅值更小；TOD台站作为参考场地的震源谱的形状明显不同，低频段（0.50～2.00 Hz）幅值更大，高频段（＞10.00 Hz）衰减更快.上述结果清楚地展示了参考场地选取对反演结果的显著影响.

图2 6个典型地震的加速度震源谱反演结果Fig.2 The inverted acceleration source spectra of six typical earthquakes

ω-2模型[14]和高频衰减模型[18]表示反演震源谱，即

式中：C为常数，定义见文献[4]；M0为地震矩；fc为拐角频率.

网格搜索方法估计地震的震源参数（M0、fc）及高频衰减参数κ，根据M0计算矩震级Mw，利用M0、fc计算地震应力降Δσ[14].

LSS、MNF和TOD台站分别作为参考场地时，利用反演震源谱得到的Mw估计值与INGV正式发布的地方震级ML的差值如图3（a）所示，整体上Mw估计值依次增大，LSS和MNF台站作为参考场地的Mw估计值接近，明显低于TOD台站作为参考场地的结果.ML＜ 5.3 时，随震级增大，Mw估计值与ML之差整体上呈下降趋势，Mw估计值从明显高于ML逐渐变化为低于或接近ML；ML＞5.3 时，Mw估计值普遍低于ML，且随ML增大，Mw与ML的差异逐渐缩小.图3（a）同时比较了TDMT和MTS提供的Mw与ML之差.TDMT和MTS提供的Mw值普遍低于ML，LSS、MNF和TOD台站作为参考场地的Mw估计与TDMT、MTS提供的Mw存在一定偏差，尤其是小震的矩震级估计偏大，这可能与距离相关的能量法截取的小震记录S波窗口过大有关，但TOD作为参考场地的Mw估计与TDMT、MTS提供的Mw的偏差最显著.

图3（b）、（c）给出了 LSS、MNF 和 TOD 台站作为参考场地的fc和Δσ估计值，参考场地不同时fc和Δσ估计值差异明显，LSS作为参考场地的fc估计值明显低于MNF和TOD作为参考场地的结果，相应地LSS作为参考场地的Δσ估计值也最小.参考场地不同时M0-fc和Δσ-Mw的变化趋势较为一致，Mw＜ 5.00时，Δσ估计值随Mw的增大而增大，LSS、MNF和TOD作为参考场地的Δσ估计值为0.11～2.90、0.32～9.97 和 0.55～13.45 MPa；Mw≥5.00时，Δσ估计值随Mw的变化不明显，LSS、MNF和TOD作为参考场地的Δσ为5.40～12.19、12.96～26.50和28.57～78.84 MPa，明显高于全球板间地震平均应力降（3.31 MPa）[19]及 Bindi等[17]估计的 1997—1998年Umbria-Marche地震序列（本研究选用地震序列的西北部）中15次Mw=4.20～6.00余震的应力降（图3（c））.Pacor等[3]利用非参数化谱反演确定的2009年L’Aquila地震序列（本研究地震序列的南部）应力降为 0.25～25.00 MPa，与本研究选用LSS、MNF台站作为参考场地得到的Δσ估计值的范围较为一致（图3（c））.

LSS、MNF和TOD台站作为参考场地的κ估计值依次增大（图3（d）），κ估计值近似服从正态分布，LSS和MNF作为参考场地的κ估计值较为接近，均值分别为 0.017 s和 0.023 s，明显低于 TOD作为参考场地的κ估计值（均值 0.056 s）.κ估计值反映了参考场地中未予考虑的可能由于近地表岩层风化等引起的高频衰减[20].一般认为岩层风化越严重，其κ值越大.Oth等[2]利用日本井下钻孔记录反演得到的井下钻孔台站的κ均值为0.015 s，LSS和MNF台站作为参考场地的κ均值与日本井下钻孔台站κ均值接近，低于一般岩石场地的κ值（0.020～0.040 s），说明 LSS 和 MNF 台站近地表岩层风化不严重，而TOD作为参考场地的κ均值较高则说明该台站近地表岩层风化可能较为严重.

图3 地震震源参数 （ Mw, fc and Δσ ）及高频衰减参数 （ κ ）估计值Fig.3 Estimated source parameters (Mw, fc and Δσ) and high-frequency decay parameters (κ)

通过分析分别选用LSS、MNF和TOD台站作为参考场地的震源谱、震源参数及高频衰减参数的试算结果，LSS和MNF台站均可作为参考岩石场地的合理选择.我们建议在选取参考场地时，首先综合考虑场地地质资料、观测记录H/V谱比曲线幅值和形状初步选择参考岩石场地，再分别以初选作为参考场地进行谱反演，根据震源谱、震源参数、高频衰减参数的试算结果，最终选取合适的参考场地.

3 单步法与两步法

本研究分别采用两步法和单步法完成2016—2017年意大利中部地震序列地震动记录的非参数化谱反演，均以LSS台站作为参考场地.两种方法分别确定的0.25～20.00 Hz传播路径衰减差异明显（图4（a）），差异大小与震源距、频率显著相关（图4（b）、（c））.近场（＜ 20 km）或低频 S 波（＜ 2.00 Hz）的传播路径衰减较为一致，高频 S 波（＞5.00 Hz）传播路径衰减的差异随震源距增大迅速扩大，高频S波（＞5.00 Hz）在远场（＞60 km）的传播路径衰减的差异显著，两者相差可达5倍.与单步法相比，两步法确定的传播路径衰减更强烈（数值更小）.采用几何扩散和非弹性衰减简单表示传播路径衰减曲线，利用奇异值分解方法给出两种方法确定的传播路径衰减的几何扩散指数n和品质因子Q，两步法和单步法的n值较为一致，分别为 -1.01和 -1.07，品质因子的差异明显，分别为 96.45f0.59和 212.45f0.61，说明两种方法路径衰减差异主要为频率相关的非弹性衰减差异.

图4 基于两步法和单步法反演路径衰减的差异Fig.4 Differences between path attenuations inverted by the two-step and one-step methods

利用单步法和两步法分别确定的6次典型地震的加速度震源谱如图5所示，两步法的震源谱幅值普遍高于单步法，大于2.00 Hz的高频段尤为明显.基于反演震源谱，利用网格搜索法估计地震震源参数（Mw、fc）及高频衰减参数（κ）（见图6（a）～（c）），两种方法Mw估计值相近，多数地震Mw估计值相差不超过0.02；两种方法fc估计值差异较小，单步法fc估计值整体上略小于两步法，多数地震单步法fc估计值约为两步法的88%～98 %，我们还注意到fc估计值的差异与fc大小有关，fc估计值的差异近似随fc增大而增大；两种方法κ估计值差异并不明显，单步法κ估计值普遍高于两步法，多数地震单步法κ估计值约为两步法的1.0倍～1.5倍，只有极少数κ估计值＜ 0.01 s的地震，单步法κ估计值约为两步法的2.0倍～6.0倍.利用Mw和fc估计值计算得到Δσ（如图6（d））和辐射地震能量Es（如图6（e）），两种方法Δσ、Es估计值差异不明显，单步法Δσ和Es估计值均小于两步法，且大多数地震单步法Δσ和Es估计值约为两步法的70%～90%.

图5 单步法和两步法反演给出的6次典型地震的加速度震源谱Fig.5 Acceleration source spectra for six typical events derived from the one-step and two-step methods

图6 基于两步法和单步法的震源参数（Mw、fc、Δσ、Es）及高频衰减参数（κ）估计值的差异Fig.6 Differences in estimated source parameters (Mw, fc, Δσ and Es) and high-frequency decay parameter (κ) by the one-step and two-step methods

非参数化谱反演的单步法和两步法给出的远场高频路径衰减差异显著，高频震源谱有一定差异，Mw、fc、Δσ、Es及κ的差异在可接受范围.两步法将观测记录的震源和场地项假定为地震大小相关的量Ma，包含了该地震所有观测记录的平均震源谱和平均场地效应，单个观测记录的实际Ma与平均Ma存在偏差，每个距离分段的传播路径衰减会吸收该距离分段的所有观测记录的Ma偏差（主要来自场地项），如果该距离分段所有观测记录的Ma偏差均值明显偏离0，此时两步法传播路径衰减项出现偏差，推断这是两步法和单步法传播路径衰减偏差明显的根本原因.由于高频路径衰减差异更显著，推断Ma在高频段的变化性更显著，可能与高频场地效应强变化性有关.为避免两步法路径衰减被场地项干扰，我们建议尽可能采用单步法，但只研究地震震源特征时，由于震源参数影响较小，也可用两步法.

4 震源参数求解中的权衡问题

常用网格搜索估计震源谱的特征参数（Mw、fc、κ），从一系列特征参数构建的理论震源谱中，利用网格搜索寻找最接近反演震源谱的一组特征参数表示的理论震源谱，作为震源谱特征参数的最佳估计.估计Mw、fc和κ时，这些参数之间可能存在的权衡会影响结果估计的可靠性.

为评估这些参数之间的权衡，本研究利用一系列Mw、fc和κ先验值，根据式（4）、（5）构建光滑的先验震源加速度谱（352个），Mw取值 3.75、4.25、5.00 和 6.00，fc取值 0.3、0.5、0.8、1.0、1.5、2.0、3.0、4.0 Hz，κ取值为 0～0.05 s范围内以 0.005 s等间隔递增的一系列数值.利用6次典型地震的理论震源谱与反演震源谱的残差修正先验震源谱，最终构建了2 112个先验震源谱.网格搜索给出先验震源谱的Mw、fc和κ估计值，Mw估计值与先验值之差（m）、κ估计值与先验值之差（e）随fc估计值与先验值之比（l）的变化见图7.

图7 震源谱特征参数之间的权衡Fig.7 Trade-offs among those parameters describing source spectra

fc先验值≤ 2.00 Hz时，fc估计值普遍偏小，而Mw估计值普遍偏大，fc估计值偏小的程度越显著则Mw估计值偏大的程度也越显著，这说明fc估计值越小则Mw估计值越大，两者之间存在一定的权衡.fc先验值为 0.50～2.00 Hz 时，fc估计值约为先验值的80%～100%，Mw估计值与先验值相差小于 0.1，Mw、fc估计值与先验值的差异均较小，说明Mw和fc之间的权衡对Mw和fc估计值的影响有限.fc先验值大于2.00 Hz时，fc估计值偏大而Mw估计值略有偏小，不论fc估计值偏大的程度如何Mw估计值偏小的程度近似不变，说明两者之间几乎不存在权衡.fc先验值≤ 2.00 Hz 时，fc和κ估计值均偏小，fc估计值越小κ估计值也越小；fc先验值＞2.00 Hz，fc和κ估计值均偏大，fc估计值越大κ估计值也越大，说明fc与κ之间存在权衡.由于κ主要控制震源谱的高频段，为避免fc与κ之间权衡的影响，确保Mw、fc的可靠估计，应尽量利用中低频段的震源谱估计Mw和fc.

针对6次典型地震，κ为最佳估计值，Mw与fc的组合构建的理论震源谱与反演震源谱的拟合方差见图8.图中：红色实线表示拟合方差最小值对应的fc值.拟合方差的外包络线均呈开口向上抛物线，fc对应于最小拟合方差值，fc＞0.50 Hz 的 4 次地震，外包络抛物线在fc最佳估计附近谷值明显，容易识别；fc≤ 0.50 Hz的 2 次地震，外包络抛物线在fc最佳估计值附近十分平坦，谷值不明显.结果说明，fc＞0.50 Hz时，Mw和fc之间的权衡对Mw、fc估计影响很小，可不予考虑，fc≤ 0.50 Hz时，需考虑Mw和fc之间的权衡，一般情况下fc较小地震的震级较大，其矩震级可直接采用震源机制解的结果以处理Mw和fc之间的权衡，确保fc的可靠估计.

图8 拟合方差随 fc 变化Fig.8 Fitting variance vs.fc

5 结 论

1）反演结果与参考场地的选取密切相关，建议首先综合考虑台站观测记录数、场地地质资料、H/V谱比曲线等初选参考岩石场地，再以初选作为参考场地，根据震源谱、震源参数、高频衰减参数试算结果，最终选取合适的参考场地.

2）单步法和两步法的远场高频路径衰减差异显著，高频震源谱有一定差异，Mw、fc、Δσ、Es及κ差异较小.为避免两步法路径衰减被场地项干扰，建议采用单步法；如只研究地震震源特征，亦可采用两步法.

3）网格搜索估计震源参数（Mw、fc）和高频衰减参数（κ）时，fc和κ之间存在权衡，为了尽可能避免此影响，建议利用中低频段的震源谱估计Mw和fc确保其可靠估计.fc＜ 2.00 Hz 时Mw和fc之间存在一定权衡，但fc＞0.50 Hz 时，此权衡对Mw、fc估计影响很小，fc较小时地震的震级较大，建议利用Mw处理Mw和fc之间的权衡，提高fc估计的可靠性.

致谢：中国地震局工程力学研究所基本科研业务费专项资助（2018B03）.