Q值描述了介质吸收特性的强弱，与地壳介质的温度、破碎程度和流体活动等密切相关，能很好地反映构造活动强度。近年来，为研究Q值的精细分布特征，层析成像技术被广泛应用于衰减结构成像并取得了大量成果，地壳衰减结构存在不均匀性基本达成共识^[1-5]。此外，研究认为地震带的活动断层区域地震波衰减很强烈^[6]，介质裂隙分布越广泛、含流体饱和度越高其衰减程度越大^[6-8]，火山裂谷带的地壳介质具有低Q值特征^[9]，高大地热流点基本分布在高衰减区域^[10-11]。因此，研究地壳介质的衰减结构成像可为活动断层分布、介质裂隙分布和流体含量、介质的热力学状态等研究提供有效信息，有利于进一步分析区域地震危险性。

太行山隆起区位于华北地区中部，1900年以来该区地震活动较弱，没有6级以上地震发生。但其作为山西隆起区与渤海湾盆地的边界，周边地区地震活动活跃，北面为张渤地震带，西面为山西地震带，东面为河北平原地震带，这些区域1990年以来发生多次6级以上地震，是地震学家重点关注区域。目前关于华北地区开展的Q值研究基本以省为区域研究对象，研究多地震多台联合反演平均Q值^[12-13]，尚未开展Q值精细成像特征的研究。本文搜集太行山隆起及邻区2009—2018年M_L2.0以上地震的波形资料，利用S波衰减层析成像方法反演该区地壳S波Q值，分析该区地壳Q值分布特征，结合地形、地震活动、大地热流等分析，探讨其与地壳衰减特征的关系。

在频率域，第j个台站上观测到的第i个事件的振幅谱A_ij（f）为

式中：f为频率；S_i（f）为震源谱；I_j（f）为仪器响应；R_j（f）为场地响应；G_ij（r）为沿路径r的几何扩散因子；B_ij（f）为衰减谱。

Brune^[14]将震源谱用长周期振幅谱Ω₀和拐角频率f_c来表达，即

采用Atkinson等^[15]提出的互相衔接的3段几何衰减函数表示几何扩散因子，即

式中：系数b₁、b₂和b₃均与频率无关；当r≤R₀₁时，对应直达波的几何衰减；当R₀₁<r≤R₀₂时，对应过渡区，在该震源距范围内，直达波中加入了在地壳内间断面和莫霍面上的反射波；当r>R₀₂时，对应多次折射反射波的衰减。

沿整个射线路径的衰减谱可以用下式表示^[16-17]:

式中：t_ij为沿射线路径的走时；Q_ij为无量纲的品质因子；$t_{ij}^*$为沿整个射线路径的衰减算子。

仪器响应可以根据台站参数直接扣除，此时振幅谱为$A_{ij}^{'}(f) = {A_{ij}}(f)/{I_j}(f)$。此外，由于大多数区域台站都是建在基岩上，其场地响应可以假定为接近1的常数。因此方程（1）可以写成

其中，$t_{ij}^*$可以写成$1/(Q(s)v(s))$沿震源i至台站j的射线路径的投影^[18-19]，即

式中：v为横波速度，ds为射线路径单元。

方程式（5）中存在3个未知变量Ω₀、f_c和$t_{ij}^*$，对某个地震事件来说，Ω₀和f_c只与震源有关，对同一事件采用多台观测谱联合反演，可以得到震源到各个台站的路径衰减$t_{ij}^*$。然后根据式（6）可采用与走时层析成像相同的方法进行衰减层析成像确定Q。

本文所用事件为太行山隆起及邻区2009—2018年记录到的可定位的震级M_L2.0以上的地震事件，其中2009年1月1日至2017年10月19日的事件波形来自国家数字测震台网数据备份中心，2017年10月20日至2018年12月31日的事件波形由喜马拉雅三期台阵提供。

S波Q值成像首先需用S波位移谱反演求取$t_{ij}^*$。对台站记录的两个水平分量波形进行带通滤波和水平校正处理，取“S波窗”和“噪声窗”。从S波开始到包括S波总能量90%的时间段内的S波即为“S波窗”，“S波窗”在近源距离只包含直达S波，在较远距离包含了从地壳内间断面和莫霍面产生的反射波，在更远的距离还包含Sn和Lg震相。在P波初至之前的波形记录上取512个点（采样频率为100 Hz）或256个点（采样频率为50 Hz）作为“噪声窗”。采用平移窗谱法进行快速傅里叶变换，将“S波窗”内的信号数据和“噪声窗”的噪声数据分别转换成观测振幅谱和噪声谱。由于所有地震台的地震计都是速度计，应将速度振幅谱转换为位移振幅谱。分别处理两个水平分量波形得到其振幅谱，即可合成S波振幅谱^[20]。为了数据的可靠性，选择频率范围1～15 Hz，信噪比至少为3。图1为山西灵丘台（SXLNQ）记录的2017年6月7日山西浑源M_L2.5地震两个水平分量波形及其合成S波位移谱和噪声位移谱，该事件满足信噪比大于3的台站记录共10个。

图  1  2017年6月7日山西浑源M_L2.5地震

据式（5），采用多台观测谱联合反演^[10]，假定一个地震事件被N个台站记录到，那么未知变量为N+2个，即Ω₀、f_c、$t_1^*$、$t_2^*$…，$t_N^*$。假定有M个频率点，第i个台在第j个频率点的观测振幅谱为$A_i^{{\rm{obs}}}({f_j})$，理论振幅谱为$A_i^{{\rm{cal}}}({f_j})$，寻找一组未知变量值，使残差$\displaystyle\sum\limits_{i = 1}^N {\displaystyle\sum\limits_{j = 1}^M {\lg \left| {A_i^{{\rm{obs}}}({f_j}) - A_i^{{\rm{cal}}}({f_j})} \right|} }$最小，这组未知变量值就是多台观测谱联合反演结果。采用遗传算法反演N+2个未知变量，使得N个台观测谱和理论谱在M个频率点上的所有残差之和最小的解即为最优解^[21]。经过噪声和几何衰减校正后，对同一震源多台观测谱进行联合反演，图2为2017年6月7日山西浑源M_L2.5地震多台观测谱联合反演结果，共得10个台站的t^*算子，不同台站位移谱的高频衰减主要受不同路径的Q值影响。从图中可以看出，有3个台站的拟合位移谱误差较大，分别为SXDAX（t^*=0.117 3）、SXSHZ（t^*=0.121 8）和SXZCH（t^*=0.104 2），这3个台站的t^*值远大于其他台站，说明有些台站的t^*算子反演误差较大，在后续处理中应当按一定规则删除（大于1倍均方差）。

图  2  2017年6月7日山西浑源M_L2.5地震多台观测谱联合反演结果

根据式（6），求得衰减算子t^*后，采用与速度层析成像相同的方法就可得到S波Q值层析成像结果。首先，选出满足误差要求的t^*衰减算子，据网格平均射线数选择合适大小的网格进行检测板分辨率实验，根据检测板分辨率实验结果调整网格直到选出合适的分辨率。之后，计算地壳平均Q_S作为初始模型，采用射线追踪伪弯曲法迭代反演Q_S值。为减弱震源深度影响，在用t^*反演Q_S值的过程中用震源距代替震中距。

根据式（5），利用遗传算法反演t^*数据，为降低解的非唯一性，要求每个事件至少能反演出5条t^*数据，共得到1 434个事件的20 985条t^*数据。由于某些台站t^*数据反演误差较大（图2），剔除大于1倍均方差的数据，进一步筛选出13 627条t^*数据，相应的台站、事件及射线分布见图3。将研究区地壳（35°～42°N，112°～116°N）在平面上划分为0.5°×0.5°的均匀网格，利用筛选出的t^*数据对Q_S反演成像，平均网格射线数为221。

图  3  研究区台站（红色三角）、事件（+）及射线分布

反演Q_S前计算研究区S波的平均速度，得到v_S=3.52 km/s。假定该区初始Q值为均匀值，据式（6）可知t^*与震中距成线性关系，利用线性最小二乘法对t^*与震中距进行线性拟合（图4），求得研究区地壳的平均Q_S值为340。将v_S和平均Q_S作为反演的初始模型。经过10次迭代反演后即可得到研究区S波Q值成像结果，其中t^*的均方根（RMS）残差由0.027 2降低到0.025 3。

图  4  t^*数据与震中距之间的线性拟合

采用检测板分辨率实验验证0.5°×0.5°网格反演结果的可靠性（图5），结果表明，研究区内（35°～42°N，112°～116°N）检测板恢复较好，反演结果可靠。由于研究区外部边缘没有地震，只有记录到研究区内事件波形的台站，部分区域检测板恢复不好，这些区域基本不在研究区域内，对反演结果影响不大。

图  5  网格检测板分辨率实验结果

经过10次迭代反演求得太行山隆起及邻区S波Q_S值（图6）。结果显示，研究区的Q_S值范围为180～520，平均值为350。由图6可以看出，太行山隆起及邻区的地壳Q_S值横向变化显著，Q_S值大小代表地壳中地震波能量的衰减程度，与沉积层厚度、构造活动及地壳介质属性密切相关^[22]。总体上看，Q_S值分布与研究区的地形和构造活动明显相关，地形越低、构造活动越强烈的地区Q_S值越低。以太行山山前断裂北段（F₁₁）和晋获断裂带（F₂）为界，西部和北部山区盆地地区Q_S值总体高于东南部平原地区，但在西部和北部总体Q_S值偏高的背景下，盆地区域大多数存在低Q_S值区，如延怀盆地、大同盆地、忻定盆地、太原盆地和长治盆地均存在低Q_S值；在东南部总体Q_S偏低的基础上，河北邢台和山东菏泽等构造活动强烈地区则分布了大片低Q_S值。

图  6  研究区S波Q值分布

从研究区1900年以来6级以上强震分布来看（图6），大多数6级以上地震位于低Q_S值（高衰减）区域边界，尤其是历史强震多发区域。如河北邢台和山东菏泽地区，且研究区2个7级以上地震（1966年邢台7.2级地震和1937年菏泽7.0级地震）均发生在这两个历史强震多发的低Q_S值区域边界，这些区域的低Q_S值可能是强震（尤其是7级以上地震）发生导致大量裂隙，且裂隙中充满流体的结果^[23]。同时，也注意到1976年和林格尔6.2级地震和1998年张北6.2级地震发生在高Q_S值的边缘，这可能是由于1900年以来这两个震源区均只发生了1次6.2级地震，不如强震多发区域地壳介质松散、破碎，Q_S值相对较高，但震源区周边有微弱的低Q_S值区，尤其是1998年张北6.2级地震震源区东面的低Q_S值区更为明显。

Q_S值与温度的关系主要表现为区域Q_S值与区域大地热流分布的联系。大地热流是地球内部热状态和热结构在地表的最直接显示。根据中国大陆地区大地热流数据汇编^[24-26]，本文给出了研究区范围内的热流点数据（热流值大于等于80 mW/m²）。由图6可以看出，研究区反演得到的介质品质因子Q_S值分布与大地热流分布有一定关系，如区域内的高热流点，大多分布在Q_S值相对较低的区域，即热流值与衰减值成负对应关系。

值得注意的是，本文假设Q值与频率无关，理论上反演得到的Q值比频率相关的反演结果略高，但Q值的分布不会因此改变，因此可以采用该结果来分析介质衰减特征^[10]。

采用类似速度层析成像的方法对太行山隆起及邻区进行Q_S值成像，分析讨论其成像特征及其与地震活动、大地热流的关系，得到如下结论：

1）Q_S值分布与研究区的地形和构造活动基本相关，总体上地形越低，构造活动越强烈的地区Q_S值越低；

2）1900年以来6级以上地震大部分发生在低Q_S值区域边界，如河北邢台和山东菏泽2个强震多发且各发生1次7级以上地震的震源区位于低Q_S值区域边界，而1976年和林格尔6.2级地震和1998年张北6.2级地震发生在高Q_S值的边缘，这可能是由于单个6级地震震源区地壳介质不如强震多发（尤其是7级以上地震）的震源区介质松散、破碎；

3）研究区内热流值与Q_S值呈负对应关系，热流值大于等于80 mW/m²的热流点大多位于Q_S值相对较低的区域。