引言

场地地表地震动主要受震源特性、地震波传播路径及局部场地条件耦合影响，自1986年Andrews提出在频域分离震源、路径及场地效应的广义反演思想与求解方法后，国内外学者开展了一系列相关研究，并对该方法进行了不断完善。除反演方法外，地震记录的数量和质量也是决定反演结果稳定可靠的关键。随着地震台网的不断建设与发展，高质量地震观测记录逐渐增加，利用地震观测记录反演地震震源参数、区域介质特性及近地表土层反应已成为众多地震工程学者的主要研究课题。

根据表征衰减的不同形式和计算原理等因素形成了4类反演方法，如假设路径衰减为已知模型的参数化广义反演方法（Andrews，1986），该方法目前应用较广泛，但存在忽略介质不均匀性造成的品质因子在部分频段计算结果异常的现象（周影等，2021）。为解决该问题，Castro等（1990）加入对几何扩散复杂性的考虑，提出利用离散变量表征衰减的非参数化广义反演方法，有效避免了衰减失真对震源谱和场地谱的影响，由于其计算分为两步，也被统称为两步非参数化广义反演方法。Oth等（2011）结合参数和非参数反演方法的优势，提出了优化计算后的单步非参数化广义反演方法。刘杰等（2003）提出结合Atkinson方法和Moya遗传算法的联合迭代反演方法，在所有台站场地反应为1的假设前提下，通过反复迭代使震源的谱残差为极小，从而确定品质因子、场地反应及震源参数，该方法虽有效提高了计算精度，但其假设条件与残差定义均给结果带来了一定不稳定影响（任叶飞，2014）。按数据类型可将我国反演方法研究分为两类，一类是基于测震记录利用Atkinson遗传算法联合反演开展的地区化实际应用，如臧阳等（2021）、周少辉等（2018）、秦敏等（2018）分别以青藏高原东北部、山东长岛及云南盈江为研究区，得到反应局部地质构造特性和构造应力环境的品质因子、震源参数及参与反演的台站场地反应；另一类是以强震动记录为数据基础，通过求解线性方程组的最小二乘解分离品质因子、震源参数及参与反演的台站场地反应（王宏伟等，2021）。

滇西南地处我国怒江、澜沧江地震带，具有地震强度大、频率高、分布密集等地震活动特点，近20年来发生过8次6级以上强震，330余次弱震与中强震。由于亚欧板块和印度洋板块长期挤压碰撞，板块构造变形强烈，使该区域成为扬子地块、保山地块和腾冲地块等多个板块汇聚地带，并沿澜沧江河谷发育了贯穿云南西部的南北走向澜沧江断裂带。因此，利用该区域地震记录开展反演研究不仅对滇西南地区的地震灾害评估和工程相关的抗震设计工作具有实际应用价值，更有利于研究区的地震预测。

本文基于滇西南地区242组强震动记录，采用两步非参数化广义反演方法得到了反映该地区介质特性的品质因子Q、26个台站的场地反应和35个地震事件的相关震源参数及其分布关系。

1.   研究数据

至2021年，云南全省内布设的强震观测台站地理分布情况如图1（a）所示，主要分为3块区域，西北地区台站主要分布于德宏州、保山市、大理市及丽江市，中部地区及部分东北地区台站主要分布于昆明市、邵通市、玉溪市及红河州，西南地区台站主要分布于普洱市及西双版纳州。本文利用西南地区（99°E～102°E，21°N～24°N）37个台站2007—2020年发生的63次3.0≤M_S≤7.6地震事件捕获的492组地震三分量记录作为初始数据库，并按以下要求对数据进行筛选：①震源距主要集中分布于200 km以内，在保证数据量充足的情况下尽可能避免数据过于离散，本文选择了震源距＜150 km的记录，如图2所示；②为避免土层非线性效应的发生，选取PGA＜0.2 g的记录；③选取被不少于3个台站捕获的地震事件及捕获数量不少于3组记录的台站。

图 1  研究区范围与数据情况

Figure 1.  The research area and selected earthquakes and stations

下载: 全尺寸图片 幻灯片

图 2  初始数据库中M_S与震源距分布情况

Figure 2.  Distribution of M_S versus hypocentral distance

下载: 全尺寸图片 幻灯片

最终，在数据库中选取了26个台站、35次地震的242组记录，震级范围为M_S3.2～6.5，震源深度为4～15 km，震源距为9～150 km。台站与震源的地理分布情况如图1（b）所示，台站与震源连线覆盖滇西南约210 km×225 km区域。

反演工作开展前，根据所选记录的信噪比对研究频段进行选择。首先，对原始的三分量记录进行基线校正；然后，对记录中S波与噪声窗进行截取，本文利用Husid函数与累计均方根函数分别拾取S波的到达时刻与结束时刻，并截取P波到时前与信号窗等长的时段为噪声窗口，同时，选用2.5%汉宁窗对信号与噪声窗口的首尾两端进行锥化处理，防止矩形窗截取引起的高频畸变；接着，在0.1～40 Hz频段上按对数等间距取300个频率点计算傅里叶谱，并采用Konno等（1998）算法进行平滑；最终，分别计算两者的水平向傅里叶谱矢量和（$\sqrt {{\rm{N}}{{\rm{S}}^2} + {\rm{E}}{{\rm{W}}^2}} $）作为信号谱与噪声谱，两者谱比即得信噪比。信噪比≥3的记录数量百分比如图3所示，由图3可知，0.1～25.74 Hz频段上信噪比≥3的记录数量超过了记录总量的75%，符合选取期望。此外，本文对记录采用0.25～30 Hz巴特沃斯带通滤波，有效频段定义在高通拐角频率的1.25倍与低通拐角频率的1/1.25倍之间（Abrahamson等，1997），综合考虑信噪比及有效频段，本文最终采用0.5～20 Hz为研究频段。

图 3  信噪比≥3的记录数量百分比

Figure 3.  Percentage of records with SNR≥3

下载: 全尺寸图片 幻灯片

对原始记录进行基线校正、带通滤波、截取S波等处理后，在0.5～20 Hz频段上按对数等间距取300个频率点计算三分量的傅里叶谱，并以两水平向傅里叶谱矢量和表示频域内的水平地面运动，即得到本研究反演所需的观测加速度傅里叶谱。

2.   研究方法

本文采用两步非参数化广义反演方法（Castro等，1990）进行研究，具体步骤如下。

（1）计算不同频率地震波与距离相关的衰减：

式中，${U_{ij}}(f{\text{，}}{r_{ij}})$为i事件在j台站的观测加速度傅里叶谱，$ {r_{ij}} $为i事件到j台站的震源距；$ {\tilde S _i}(f) $为与i事件震级相关的标量，从物理意义上来说，其不仅表示i事件的震源谱，且包含了记录到i事件台站的平均场地反应；$A(f{\text{，}}{r_{ij}})$为与震源距和频率相关的离散变量，用于描述地震波在传播介质中的衰减特性，假设地震波于参考距离$ {r_0} $处开始发生衰减，即假设$A(f{\text{，}}{r_0}) = 1$。

选取反演数据库中最小震源距作为参考距离（$ {r_0} $=9.4 km），并以固定距离10 km为分段划分为N个区间进行计算。此外，引入控制$A(f{\text{，}}{r_0}) = 1$的加权因子$ {\mathrm{\omega }}_{1} $和抑制衰减变量随震源距单调递减、避免复杂场地条件对其影响的加权因子$ {\mathrm{\omega }}_{2} $，加权因子取值均需根据不同反演数据及条件进行试算，本文不再详述。

加入加权因子后可将式（1）写成线性方程组Ax=b，其中，A为与路径、震源相关的系数矩阵，x为目标矩阵，b为观测傅里叶谱组成的已知矩阵。最终利用x=(A^TA)⁻¹A^Tb计算得到不同频率地震波随震源距变化的$A(f{\text{，}}{r_{ij}})$。

（2）利用第一步确定的衰减结果$A(f{\text{，}}{r_{ij}})$对路径项进行衰减修正，并将震源项与场地项分离：

式中，${R_{ij}}(f{\text{，}}{r_{ij}})$为经衰减校正后的傅里叶谱，包含了剔除路径项后的震源项与场地项结果；${S _i}(f)$为i事件的震源谱；$ {Z_j}(f) $为j台站的场地反应。

本文通过确定参考场地以降低自由度，从而解决震源项与场地项之间的权衡问题，并使用最小二乘法求解得到场地谱与震源谱。

3.   研究结果

3.1   场地反应

为限制反演的场地反应结果，在26个台站中选取场地类别为基岩场地的景谷台（53 JGX）作为目标参考场地，并根据水平/垂直谱比法（HVSR）对该台站的地震记录进行计算，以此判断其是否为理想参考场地。HVSR谱比计算结果如图4所示，由图4可知，该台场地反应在研究频段上无明显放大或缩小，曲线近似平坦，无清晰峰值，表明该台站是作为参考台站的良好选择。

图 4  景谷台HVSR谱比计算结果

Figure 4.  HVSR results at the 53JGX

下载: 全尺寸图片 幻灯片

根据式（3）确定了26个台站在0.5～20 Hz频段上的场地反应，并将其与HVSR结果进行对比，如图5所示。由图5可知，非参数化广义反演方法确定的场地反应相较于HVSR法整体偏高，部分台站在1 Hz附近的低频段上由HVSR法确定的场地反应略高于广义反演方法，大于该频段上2种方法确定的场地反应谱形相似。部分台站通过2种方法均可确定峰值频率，且二者的结果较一致，这与任叶飞等（2013）利用以上2种方法确定汶川地区场地反应给出的相关结论相符。而部分台站通过HVSR谱比计算结果无法有效辨识峰值频率，这可能与台站可用记录数量较少有关，但这些台站在非参数化广义反演方法给出的结果中可辨识到明显峰值，对于记录量不足的台站来说，利用非参数化广义反演方法确定场地卓越频率较HVSR法可靠。

图 5  2种方法确定的场地反应结果对比

Figure 5.  Site effect results by using GIT and HVSR method

下载: 全尺寸图片 幻灯片

3.2   品质因子

由式（1）计算得到的衰减变量结果是地震波经一系列复杂现象的综合体现，其中包含了所有导致地震波能量衰减的影响因素，从性质上可分为几何因素和物理因素。其中，几何因素即因波阵面扩展导致的地震波能量衰减，通常用几何扩散函数G（r）描述；物理因素即地震波通过界面折射、反射及散射时造成的衰减与地震波在岩石中传播时发生的内在能量耗散，统称为非弹性衰减，通常用品质因子Q进行估计。两类因素的相互耦合导致用方程表述路径项即地震波衰减时需考虑权衡问题，为此已有学者提出了不同的解决方法，如约束某一震源距范围内几何扩散模型及几何扩散指数n，通常包括线性模型（Oth等，2008；Bindi等，2009）、双线性模型（Castro等，2013；Wang等，2020）、三线性模型（Pacor等，2016），并根据选取的不同衰减模型固定几何扩散指数n；预先约束品质因子Q与频率的关系（Dutta等，2001）；根据频率同时考虑n和Q（Bindi等，2004；Wang等，2019）。本文采用以下模型同时考虑几何衰减与非弹性衰减，并计算与频率相关的品质因子Q：

式中，G(r)用模型$ {({r_0}/r)^n} $表示，v_S为S波速度取3.5 km/s。对式（4）两边取对数后可写成线性方程组，如式（5）所示。

本文采用SVD分解同时求解不同频率地震波的几何扩散系数和品质因子，并用$Q = {Q_0}{f^n}$对每个频率点的反演结果进行回归拟合，最终得到了该地区Q值拟合模型：Q=94.23f ^0.43，如图6所示。

图 6  品质因子Q拟合结果

Figure 6.  Fitting results of quality factor

下载: 全尺寸图片 幻灯片

与许亚吉等（2020）基于遗传算法的Atkinson法拟合的金沙江-红河断裂以西的品质因子Q值模型及苏有锦等（2006）在云南不同地区的品质因子Q值模型结果相比，本文确定的拟合结果偏小。主要原因可能包括：品质因子反映的是一定区域内地震波传播路径上的介质特性，如果开展反演的研究区尺寸不同，地震记录的震源距范围等因素不同，会导致品质因子出现差异；相比其他学者的研究，本文选取的记录地震射线长度短，震源深度小，反演得到的品质因子主要反映了浅层介质对地震波的影响，而由于浅层介质的横向不均匀性大于深层介质，地震波在浅层的衰减更明显，因此浅层介质的品质因子较小。

3.3   震源参数

基于Aki于1967年对震源模型的研究，理论位移震源谱可用地震矩M₀、拐角频率f_c等参数表示为：

式中，C为与震源处介质特性相关的常数；$ {\Re _{\theta \varphi }} $为S波辐射花样系数，设置为0.55；V=$ 1/\sqrt 2 $，为地震波水平分量含剪切波能量的比例；F为自由地表放大系数，设置为2；ρ为震源处介质密度，取云南地区平均介质密度2.9 g/cm³。

本文利用网格搜索法在一定理论范围内确定合适的地震矩M₀、拐角频率f_c及震源模型参数$ \gamma $，使研究频段内观测谱与理论谱良好拟合，即二者相对面积差值最小：

式中，m为研究频段内设置的频率点数；$ {{\text{S}}_{{{{\rm{theory}}}}}}({f_i}) $为式（6）表示的理论位移震源谱，$ {{\text{S}}_{{{{\rm{observe}}}}}}({f_i}) $为反演得到的观测位移震源谱。

搜索每个地震事件震源参数前，以面波震级M_S±1、变化步长0.01限制矩震级范围，并通过式（9）确定地震矩搜索范围。拐角频率f_c搜索范围设定为0.01～5 Hz，变化步长0.01。震源模型参数$ \gamma $搜索范围设定为2±0.5，变化步长0.1。

根据以上步骤对参与反演的35次地震事件进行了理论与观测谱拟合，部分拟合结果如图7所示。

图 7  位移震源谱拟合结果

Figure 7.  Fitting results of partial displacement source spectrum

下载: 全尺寸图片 幻灯片

根据以上搜索原理及方法得到了参与反演的35次地震事件的相关震源参数，其中拐角频率、地震矩、应力降分别 为0.20～1.25 Hz、0.04～32.175×10¹⁶ N·m、0.02～1.5 MPa。在双对数坐标下对地震矩和拐角频率进行最小二乘拟合，得到拟合式为$ \lg {{{f}}_{\text{c}}} = (3.36 \pm 0.47) - (0.23 \pm 0.03){{{M}}_0} $。由图8（a）可知，随着地震矩的增大，拐角频率呈逐渐下降趋势，二者呈现出明显的负相关关系。当数据样本足够时，拐角频率不仅可反映地震矩（地震释放的能量），还可说明一定区域内的应力状态。为利用地震矩和拐角频率得到相应的应力降情况，本文根据Brune圆盘模型（Brune，1970）进行了计算：

图 8  震源参数间的分布关系及拟合结果

Figure 8.  Distribution and fitting results of some source parameters

下载: 全尺寸图片 幻灯片

地震矩与应力降、震源深度的关系如图8（b）所示，由图8（b）可知，随着地震矩的增大，应力降未明显增大或减小，表明应力降与地震矩无明显的相关性；震源深度与应力降的相关性不明显，这与左可桢等（2021）在四川省长宁县开展的震源参数研究结果一致。

对地震矩与面波震级在半对数坐标下进行拟合，结果如图8（c）所示。由图8（c）可知，地震矩与面波震级呈正相关关系，拟合关系为$ \lg {{{M}}_0} = (11.92 \pm 0.41) + (0.78 \pm 0.09){{{M}}_{\text{S}}} $，随着面波震级的增大，地震释放能量增高，地震矩与面波震级的线性关系不仅体现了该地区地下介质状态，也反映了地震活动性的强弱。在不考虑面波震级与里氏震级的转换关系下，本文拟合结果与周少辉等（2020）、杨晶琼等（2010）得到的结果相似。

4.   结果与讨论

本文利用两步非参数化广义反演方法反演得到了滇西南地区S波品质因子Q、26个强震台站的场地放大结果及35个地震事件的加速度震源谱，结合网格搜索法确定了各地震事件的拐角频率、地震矩等震源参数。通过将本文研究结果与其他学者研究结果进行对比，得出以下结论：

（1）在0.5～20 Hz频段内，非弹性衰减系数品质因子Q与频率的拟合关系式为Q=94.23f ^0.43，低Q₀值和回归系数说明了地震波在滇西南地区衰减吸收较快、频率依赖性偏弱的特点，同时也反映了研究区由于地壳运动导致断裂及褶皱发育较多，构造变形活跃，有较强的地震活动性。

（2）本文利用两步非参数化广义反演方法和HVSR法确定了26个台站的场地反应，二者曲线走势相似，2种方法给出的部分台站场地放大频段较一致，但个别台站的对比结果略有不同。总体来说，2种方法均可根据明显的峰值较准确地估计场地卓越频率，但非参数化广义反演方法确定的场地反应偏大，这可能与竖向地震动放大或区域地质构造有关，需结合建台报告等地质资料进行深入研究。

（3）通过网格搜索法估计了各地震事件的地震矩、拐角频率及震源模型参数，并计算了相应的应力降。应力降为0.02～1.5 MPa，平均值为0.30 MPa，处于较低水平，表明该地区无异常应力增强。其中，地震矩与拐角频率呈负相关关系；应力降与地震矩、震源深度均无明显依赖性；地震矩与面波震级呈正相关关系。

（4）总体上来说，在广义反演的理论支撑下，围绕地震波在传播介质中的衰减模型问题，两步非参数化广义反演方法相比传统的参数化反演，不仅给出了包含所有影响地震波传播的综合因素概念，还对传统的衰减表达方式进行了更符合实际情况的函数优化，可在简化反演计算的同时，有效避免反演结果不符合物理意义的情况出现。但是，两步非参数化广义反演方法中第一步从观测傅里叶谱分离地震波随频率变化的衰减函数时，计算中包含了同一地震事件所有不同台站记录的平均意义震源谱和场地反应，这可能会对后续结果造成一定影响，还需通过进一步研究减小误差，提高结果的可靠性。

致谢 感谢中国地震局工程力学研究所及中国地质科学院地质力学研究所张斌助理研究员为本研究提供的数据支持。