Q值是描述地球介质的衰减特性、反应构造活动的重要参数.例如在有岩浆活动的火山区，构造活动区的岩石破碎带，充填流体的空隙和裂隙等区域，当地震波穿过时，能量将发生强烈衰减而具有低Q值，而在一些构造稳定的克拉通地区则衰减较少而具有高Q值.动态Q值变化则能够反应构造活动的变化趋势，在构造活动强烈，地震多发的地区，每次地震的发生都是应力不断积累到释放的过程，在该过程中裂隙也会随着应力的变化而发生闭合或开裂变化，从而导致Q值发生变化.因此研究Q值动态变化能够有效提供应力、裂隙、地下流体等的变化信息，为进一步理解地震的孕育发生提供证据.

振幅层析成像方法先后在华北^[1]、川滇^[2]和整个中国大陆地区^{[3, 4]}，以及日本群岛^[5]都取得较好的效果，但上述研究均隐含假定Q值不随时间变化，用所有的资料仅获得一个Q值横向变化的静态图像.实际上，特别是在构造活动强烈的地区，Q值应该是变化的，因为地震在不停地发生，其主要控制因素-应力环境在不断变化.川滇地区地处青藏高原东南缘，南北地震带的南端，地震频繁，构造活动强烈，受欧亚板块和印度板块强烈碰撞的影响，青藏高原物质东移并受四川盆地阻挡，在此发生强烈变形和旋转运动.强烈的地震活动和构造变形，可能使得地壳Q值随时间的变化幅度比较大而被观测到.因此川滇地区是研究Q值动态变化的理想区域之一.同时川滇地区台网密集，地震观测数据丰富，近二十多年来，四川省和云南省积累了完整的地震观测报告，为Q值的动态变化研究提供了数据基础.

本研究尝试利用川滇地区不同时段的数据，分别进行Q值层析成像，从而获得Q值的动态变化趋势.首先仔细分析了单次成像结果的误差和可靠性，也就是确认在数据量较少的情况下，振幅层析成像能否给出误差较小，质量较好的成像结果.然后才有基础分析Q值的动态变化.因此如何评价反演结果的误差和分辨率就显得非常重要.本文在详细论述利用各种技术方法分析振幅层析成像结果的误差和分辨率的基础之上，通过分时段成像获得Q值的动态变化信息.

振幅层析成像方法已经在参考文献[1~5]中作了详细论述，本文为了下文描述方便，现将方法简单陈述.该方法重新推导了振幅谱公式，使其适应振幅成像的需要.振幅谱公式通常可以写为^[6]：

(1)

其中A_ij（f）是地震波从震源j到台站i之间的观测振幅谱，f是该振幅对应的频率，O_j（f）是震源处的振幅，S_i（f）是观测场地和仪器的响应，G_ij（R）是路径距离为R时的几何扩散项，其值为R^k，k是几何扩散因子，B_ij（f，R）是地震波传播时的能量衰减项，它主要由介质的品质因子Q决定，是频率的函数，可用幂函数描述为^[6]：

(2)

其中c=πf^1-η/v，Q=Q₀f^η，v是S波的速度，Q₀和η分别是S波1Hz时的品质因子Q和频率的幂指数.将公式（2）带入公式（1），两边取对数，并整理后得：

(3)

式中台站项a_i=lnS_i（f），震源项b_j=lnO_j（f）.公式（3）与Pn的走时方程非常相似^[7]，因此，可以通过类似Pn成像技术求出Q₀的横向变化.这里假定地壳S波速度均一，射线为直线.通过线性拟合-cR即可以得到平均Q₀值.拟合直线即是区域的平均衰减模型，将方程（1）减去平均模型，并将模型离散化成二维网格，即可获得离散化的扰动方程.通过引入了Laplace光滑算子^[8]的经典LSQR方法^[9]即可求得扰动方程的最小二乘解，即Q₀值的横向变化和台站项及震源项.更详细的公式推导和方法叙述请参见文献[1~5].

本研究取范围为20°N~35°N，95°E~110°E，作为研究重点的川滇地区处于该区域的中部.所用水平分量最大振幅和周期的观测数据主要选自中国地震年报、四川省地震台网观测报告及云南省地震台网观测报告.将这些数据分成6个时间段，每个时间段大致为2年，有大约一年的重复.每个时段所选的数据都必须满足相同的条件：①每个事件至少被4个地震台站记录到；②每个地震台站至少记录到4个事件；③直达S波的卓越周期在0.4~2.0s之间；④震源深度小于15km.最后将所选的数据统计情况列于表 1中，表 1中的标准差是上文方程（3）中的减掉平均模型后的残差的统计结果.每个时间段的月数稍有差异，这是为了使得每个时间段的有效射线数基本相同，射线的分布也尽可能一致，从而可以在反演中使用相同的反演参数，特别是采用相同的阻尼系数.

表 1(Table 1)

表 1 1999~2006年6个时间段的数据统计表 Table 1 The data statistics for 6 time segments from 1999 to 2006

表 1 1999~2006年6个时间段的数据统计表 Table 1 The data statistics for 6 time segments from 1999 to 2006

图 1给出了6个时段的射线分布和网格射线数图像.从图中可以看出，不同时段的射线分布非常相似，网格射线数也很接近.数据残差均在0.78~0.86之间，相互之间差别很小.这些均为下一步反演采用相同的反演参数提供了基础.

图 1

Fig. 1

图 1 各时间段的射线分布（左）和网格射线数（右）图像 (a)2k01;(b)2kl2;(c)2kl3;(d)2k34;(e)2k45, ⑴2k56.左图为射线分布，其中+表示地震，红色三角形表示台站，蓝色细线表示地震射线；右图为网格射线数图像，其中颜色对应于穿过大小为10'×10'的网格内的射线条数，蓝色表示射线密集，红色表示射线稀疏. Fig. 1 Ray paths (left) and its density image (right) for each time segment

为了详细分析单时段成像结果的误差和可靠性，我们选取了2k56时段数据，即2005年1月~2006年12月共两年的数据进行成像和结果误差分析.也就是要确认在大约14000条射线数量的情况下，振幅层析成像能否给出误差较小，质量较好的成像结果.以下分别从分辨率测试、观测误差、参数取值、Bootstrap抽样等方面进行成像可靠性的论证.由于各时段射线数量和分布都非常接近，分辨率测试和误差估计也非常接近，因此仅选取2k56时段进行陈述.在实际反演中，方程（3）中的各参数分别取值为k=1.0，η=0.5，v=3.6km/s，阻尼系数取2000，网格大小取10′×10′，对结果的评估也是基于这些参数来计算的.

分辨率测试的目的是用来检验数据的数量和分布是否能够反演出分辨尺度大小的异常体.它是结果误差分析的基础，也就是说在能够分辨的区域内讨论结果误差才是有意义的.

由于LSQR算法不能直接给出解的分辨矩阵，通常采用检测板方法对结果的分辨率进行检验，即分别用不同尺度的检测板，施加幅值为200，正负相间的函数作为Q₀值扰动模型，根据该模型计算每条射线的残差作为合成观测数据，然后用与实际数据反演相同的参数进行合成数据反演，看能否恢复给定的Q₀扰动模型.如果能恢复即说明解的分辨率能够达到合成数据的检测板尺度.在以往的研究中，阻尼系数也采用与实际数据反演同样的数值，这样就容易引起实际分辨率被低估.因为在实际数据中都存在观测误差、公式近似造成的误差等，并且观测误差可能是非正态分布的，而在检测板测试中的合成数据是精确的，添加的观测误差也基本都是正态分布.在LSQR算法中，阻尼系数的选取一般与数据的误差有关，误差越大则阻尼系数也越大，实际数据和合成数据具有不同的误差水平，因此需要取不同的阻尼系数.本文在反演合成数据时是这样选取阻尼系数的：就是使得合成数据的反演结果的最大值与初始扰动模型的最大值基本一致.

图 2分别给出30′，20′的分辨率图，其反演过程中对合成数据都添加了10%的随机误差，阻尼系数分别取80和30，可以发现该值比实际资料反演的阻尼系数小很多（实际资料反演中阻尼系数取2000）.从图中可见，在研究重点地区分辨率能够达到30′，对比图 1射线数分布图可以发现，射线数超过50条的区域分辨基本可以达到20′.边缘地区射线覆盖较稀且射线交叉性弱，分辨率较低，这也是将实际计算反演区域比重点研究区的范围取更大的原因.

图 2

Fig. 2

图 2 检测板分别为30'(a)，20'(b)的分辨图像 Fig. 2 Resolution from checkerboard test with size 30'(a) and 20'(b)

本文采用在原始数据上加噪音的方法测试反演结果的稳定性，即保守地假定原始观测数据是精确的，然后再加噪音来评估反演结果.具体做法是，分别在原始数据上增加20%和40%的均值为0方差为1的高斯噪音，即

gasdev为产生均值为0方差为1高斯分布的函数，A为原始振幅数据，A_noise为加噪音后的振幅数据，α为噪音水平的值，这里分别取20%和40%.然后采用与实际数据反演相同的参数，对增加噪音后的数据进行反演.进行100次加噪音合成数据-反演操作，即可以获得统计平均和方差，方差可看作是噪音引起的结果误差.这种估计方法是一种保守的估计，因为原始数据已经含有误差.这里没有考虑震源辐射花样，主要是因为，对于高频S波而言，其振幅并不满足双力偶辐射花样，而带有很大的随机性^[10].本文采用的大多是300km以内的小震，频率相对较高.另外，由于采用数量巨大的地震数据，辐射花样的影响在一定程度上能够被光滑掉.此处采用相对误差来模拟观测误差主要是因为，仪器的放大倍数每年只校定一次，放大倍数可能会因温度、湿度等外界参数的变化而变化；另外就是采用相对误差在取对数时很容易处理成加减法.

图 3分别给出了添加20%和40%噪音时的Q₀值误差分布.增加20%的噪音时，结果误差不超过8；增加40%的噪音时，结果误差不超过18（平均Q₀值的6%）.之所以最大加40%的噪音，是考虑到方法中没有考虑辐射花样的影响，实际的测量误差应该是比较小的.从试验结果看，40%的噪音对结果的影响仍然是不大的.主要原因是在反演中需要对振幅取对数，这样即使噪音水平较高，取对数后仍然是一个比较小的值.同时由于大量数据叠加的效果，使得正态分布的噪音的影响很小.

图 3

Fig. 3

图 3 噪音水平分别为20%(a)和40%(b)时反演的结果误差 Fig. 3 Inversion errors for noise level 20% (a) and 40% (b)

S波的速度v和频率的幂指数η均会影响到最终反演的Q值，LSQR算法中的阻尼系数也会影响到Q值变化的幅度，进而影响地震的定位精度.

在川滇地区，台站比较密集，地震的定位精度大部分都在I类精度，即误差小于5km.震源位置的误差主要反映在公式（3）中的震中距R上，从Q₀^-1cR项可以看出，对于一条射线而言，震中距R或c值百分之几的变化会导致Q₀百分之几的变化.川滇地区震中距统计表明R=107±60km，5km的震源误差会产生平均约5%的震中距误差.频率的幂指数η一般在0.3~0.9之间，频率范围在0.4~2Hz之间，从单条射线来看，相对于η=0.5，最大可导致c值约30%的变化，取其一半15%作为平均误差.如果将震中距R和c值合并考虑，则平均约产生20%的误差，对于最后反演的残差方程（4）而言，单条射线时ΔQ₀的误差平均应为20%，即为200×20%=40，但对于大量射线叠加而言，其误差应该小的多.

我们对上述分析进行模拟试验，将c值加上20%的噪音（与cR乘积加上20%噪音是等效的），即取c值为c*（1+20%*gasdev）.进行100次反演计算，结果显示Q₀最大误差不超过5.所以定位误差和幂指数η对反演结果的影响在数据量较大时是可以忽略的.

S波的速度v与Q₀在公式（3）是乘积的关系，速度的增加会导致Q₀值的降低.从一些层析成像的研究结果来看，上地壳的波速变化最大大约5%~6%，因此能造成Q₀值的最大误差也仅约为5%~6%.如果我们假设S波速度是不随时间变化的，那么波速不会影响Q₀值的动态变化趋势，只会稍微影响Q₀的绝对值.

文中几何扩散因子取k=1，只是将理论的几何扩散去掉，并没有将波散去掉，因此本文得到的衰减包括内摩擦造成的能量损失和波散，将二者完全分开还有一定困难，但我们尝试过取不同的k值，得到的结果显示Q值的相对高低基本不变，只是平均Q值有明显的变化.

借鉴Lg波在1Hz附近的频率衰减指数，文中取η=0.5，因为ML振幅数据的频率集中在0.4~2Hz，S波最大振幅一般在直达S波到时之后，已经在Lg波的范围内，所以借鉴Lg波的η值，但在称呼上仍沿用了横波（或S波）的叫法.另外，因为原始振幅数据的频率范围非常窄，所以η的变化对最终Q值分布的影响也是不大的.这些参数虽然对Q值有一些影响，但我们讨论Q值动态变化时都是基于相同的参数取值，甚至数据量和数据分布也力求做到相同，如此才具有可比性.

另外一种检验反演结果的方法就是bootstrap方法.由于LSQR算法不能直接给出解的方差矩阵，解的误差可以通过bootstrap方法^{[11, 12]}求得.Bootstrap方法是由Efron首先提出的一种统计估计方法，其基本思想是在样本组中再抽样，就是从n组原始数据中随机再抽取n组数据进行反演，得到一个结果x_i，重复进行B次试验就可以得到B个计算结果，进行统计计算就可以求得试验的标准差即解的误差.Efron & Tibshirani认为对于标准差的估计，B大约取100是一个比较合理的数值^[13].本文的求解过程中也取B=100.

图 4给出了用Bootstrap方法计算获得的Q₀值误差.其最大误差不超过24（相当于平均Q₀值的8%）.由于每次Bootstrap随机在母本中进行重新抽样，所抽取的样本反映不同的数据分布，因此Bootstrap方法获得的结果误差，不仅包含了原始数据观测误差的影响，而且反映了不同数据分布造成的结果误差.Bootstrap方法获得的结果也显示成像结果的误差比较小.图 3和图 4都是采用和实际反演中相同参数得到的，并且图 3中用到的射线与实际反演用到的数据完全相同，图 4中用到的射线因是随机抽样，与实际反演用到的数据分布大致相同，因此估计误差的分辨率与实际数据反演的分辨率基本相同.在射线较少的边缘地区，误差估计没有分辨，是不可信的.

图 4

Fig. 4

图 4 Boostap方法获得的反演结果误差 Fig. 4 Inversion errors from Bootstrap method

综合分析观测数据误差、地震定位误差、参数取值误差可以发现，在某均值附近这些误差是正态分布的话，即使误差较大，由于大量数据叠加的作用，使得最终结果的误差仍然较小.但参数的整体偏离，会使得平均Q₀值也发生偏离.这要求我们在动态变化反演中各参数尽可能相同，这样才具有可比性.上述分析结果表明，各种误差造成的结果变化基本都在10%以内，因此动态变化反演中Q值随时间的变化明显大于10%的区域，基本可以认为是由于介质变化造成的.

在计算过程中，我们取几何扩散系数k=1.0，因为直达S波是典型球面扩散的体波；另外取Q随频率指数衰减系数η=0.5^{[14, 15]}，由直达S波的到时计算得到地壳平均S波速度v=3.6km/s.将地壳横向离散成10′×10′的网格，通过层析成像的方法反演ML振幅残差即可获得S波Q₀值的横向变化，相对的台站增益以及震源尺度校正.反演过程中阻尼系数起到比较关键的作用，阻尼系数的选取除了满足绝对Q值大于0的物理约束外，还要求在相对均一的高Q值的四川盆地地区，反演获得的Q值水平与该地区直线拟合获得的平均Q值水平相当.通过一系列不同阻尼系数的测试，最后选取阻尼系数为2000.在各个时段的反演计算过程中，均取相同的参数.通过拟合几何扩散校正后的ML振幅观测值与折合震中距之间的线性关系，获得各个时段的平均S波Q₀值，列于表 2中.

表 2(Table 2)

表 2 各时段地震个数统计表1) Table 2 Earthquake number and average Q0 value in each time segment 1)

表 2 各时段地震个数统计表1) Table 2 Earthquake number and average Q0 value in each time segment 1)

1)四川省和云南省地震观测报告.

图 5给出了各个时段的Q₀横向变化图像，图中用颜色的变化表示Q₀的变化，红色表示比平均Q₀低的地区（强衰减），蓝色表示高Q₀地区（弱衰减）.为了便于比较，各个时段的Q₀值均归算到以290为均值的横向变化上来.从图中可以看出，各个时段所得到的整体结果基本相似，均是高Q₀值存在于四川盆地和云南保山-大姚地区，而低Q₀主要存在于川滇菱形块体的西半部和东南边界带.这也说明反演方法和成像结果是比较稳定的.成像结果在整体上与采用全部数据反演的结果^[2]是一致的，与利用波形数据获得的Lg波Q₀值^[16]的基本格局也是相同的.反演获得的Q值主要反映上地壳脆性层S波的衰减特征，因为所用数据震中距较小，统计表明R=107±60km，大部分振幅数据都是S波振幅，同时选用的都是震源深度小于15km的地震.因此Q值横向变化图像基本上是整个脆性层S波衰减的平均效果.

图 5

Fig. 5

图 5 6个时间段的Q0值横向变化图像.红色为低Q0异常，蓝色高Q0.灰色线条为活动断裂[17]，黑色线条为河流.各个时段的Q0值均归算到以290为均值的横向变化上来.色标是相对于平均Q0的扰动值，对应绝对Q0值范围为90~490. Fig. 5 Lateral Q0 variation image for 6 time segments

反演中同时也获得了相对的台站增益以及震源尺度校正，各时段的结果也基本相似，鉴于篇幅有限，此处暂未详列.反演后的残差均有大幅度降低，6个时段分别为0.566，0.567，0.572，0.559，0.555，0.567.平均降幅31%.

通常认为是应力不断积累造成的，当其差应力超过岩石或断裂的强度就会发生地震.应力积累的外部表现就是地震活动性增强，地震个数增加.研究发现区域平均Q值与地震活动性有非常强的正相关关系.表 2给出了6个时段不同震级的地震个数与整个区域平均Q值的对应关系.图 6更直观显示了二者的密切关系，即地震个数越多，区域平均Q值越大.震级越大这种关系越明显，在2003~2004年期间，地震个数最多，平均Q值也最大.

图 6

Fig. 6

图 6 平均Q0值的动态变化及其与地震个数的关系 Fig. 6 The correlation between average Q0 value and earthquake number

地壳中的地震波衰减，一般认为是由于岩石中存在微裂缝，并且裂纹内充填流体是造成波衰减的重要原因^{[18, 19]}.裂纹表面和晶粒边界产生滑动摩擦使得部分地震波能量由于克服摩擦产生滑动而转化为热能耗散掉^{[20, 21]}，而裂纹内流体的存在使这种滑动更容易产生，并且流体的流动将消耗地震波能量以克服液体的黏滞阻抗，造成地震波的衰减^{[18, 22]}.实验室研究表明，不论是干岩石，还是饱水或部分饱水岩石，随着应力水平的增加，Q值都会升高^{[23, 24]}.其内部物理机制一般认为，应力的增加会导致裂隙的闭合，裂纹密度降低，闭合的裂隙是不会产生摩擦的，因此由摩擦而产生的地震波的衰减也随之减少，即Q值升高^[21].本文的研究结果也似乎符合这一认识，地震个数增加，地震活动性增强是区域应力积累的效果，而高Q值与此有非常好的对应关系.

由于每个时间段为2年，时间较长，所获得的平均Q值或者Q值变化是2年的平均结果，并不能直接反映是地震前还是地震后造成的.但也存在另外一种可能，即Q值的变化也可能是地震引起的结果.当岩石临近破裂时，随着应力的增加岩石的体积会不断增大，即扩容^[25]，岩石内的微裂隙会不断产生、增大、扩展，直到贯通破裂.但地震发生之后的过程与扩容阶段刚好相反，地震释放能量后应力急剧下降，相应岩石体积减少，岩石内原本张开的裂隙也将缩小或闭合，从而导致衰减减少，Q值升高.图 6的观测结果似乎也符合这种解释.

区别这两种解释的根本之处在于判断川滇地区的岩石变形是处于扩容前的线性变形阶段，还是处于扩容后的非线性变形阶段.川滇地区是南北地震活动带的中段，中强地震较多，因此Q值的变化似乎用扩容理论来解释更合理一点.这也是下一步工作很值得进一步研究的地方.

本文获得了6个时段的Q₀值横向变化，很容易计算出不同时段间Q₀值的变化量来，图 7（a~e）给出了相邻时段的Q₀变化图像，图 7（f，g）给出了2张相间时段的Q₀变化图像.图中红色表示Q值降低，蓝色表示Q值升高；图 7（a~e）中的红圆圈为结束时段最后一年的震级大于4级的地震分布，粉色圆圈为开始时段第一年的震级大于4级的地震分布，图 7（f、g）中的红圈为结束时段两年的震级大于4级的地震，粉色圆圈为开始时段两年的震级大于4级的地震分布，圈的大小与震级成正比，圈越大表示震级也越大.

图 7

Fig. 7

图 7 不同时间段间Q0值的变化图像 图中红色表示Q0值降低，蓝色表示Q0值升髙异常.灰色线条为活动断裂[17]，黑色线条为河流.黄色和蓝色线条分别为一级和二级活动地块边界[26].图（a~e)中的红圆圈为结束时段最后一年的震级大于4级的地震分布，粉色圆圈为开始时段第一年的震级大于4级的地震分布，图G、g)中的红圈为结束时段两年的震级大于4级的地震，粉色圆圈为开始时段两年的震级大于4级的地震分布，圈的大小与震级成正比，圈越大表示震级也越大. Fig. 7 The changes of Q0 value between different time segments The changes of Q0 value correspond with color, red represents drop and blue denotes increse of Q0.The gray and dark lines are the active faults[17] and rivers, The yellow and blue lines are the boundries of first and second active blocks[26].The red and pink circles represent earthquakes with magnitude larger than 4.0 in former and later time segment, respectively, in figure (f) and (g), but in figure (a~e), only earthquakes happened in firt yearfor former time segment, and last year for later time segment, have been ploted.The circle size is proportional to the earthquake magnitude, the biger circle, the larger earthquake.

从图 7中可以看出，不论是起始时段的地震还是结束时段的地震，地震的分布一般并不集中在Q值增加或降低的中心区域，而是大部分集中在增加和降低峰值区域之间的过渡区域.即存在于红色和蓝色区域的中强地震较少，而存在于其他过渡颜色区域的地震较多.在图 7（f、g）中，由于时段间隔较长，Q值变化幅度较大，这种分布关系显得似乎更加清楚.虽然二者之间的内在机制尚不清楚，但这种分布关系有助于地震危险区的划分.

本文在深入分析探讨了观测误差、参数取值等对反演结果的影响，利用检测板方法测试M_L振幅成像结果的分辨率，并利用原始数据加噪音的方法和Bootstrap方法进行误差分析，在此基础之上，将川滇地区的数据分成6个时段分别进行成像，获得上地壳Q值的动态变化图像.主要表现在以下几个方面：（1）原始数据的观测误差和忽略震源辐射花样对成像结果的影响是比较小的，40%的噪音所造成的结果误差最大不超过6%；（2）用Bootstrap方法获得的结果误差不超过平均Q₀值的8%，说明振幅层析成像方法在具有大量数据的川滇地区所获得的结果是稳定可靠的，误差是比较小的；（3）分辨率测试发现，川滇地区在射线超过50条的区域，分辨率能够达到20′；（4）分别对每个时间段进行成像，从而获得Q值变化的动态图像.每个时段的区域平均Q值与中强地震的个数有明显的正相关关系，即地震个数越多，则区域平均Q值越大，这可能是由于区域整体应力变化的结果.通过分析地震的分布和Q值动态变化的关系还发现，大多数中强地震不是分布在Q值变化最大的区域，而是分布在Q增加和降低最大区域之间的过渡区域，这可能是由于差应力变化在过渡区最大，更容易触发地震造成的.本文对Q值动态变化研究尚处于探索阶段，有些结论可能还需要更多的资料，更多的方法来互相验证，甚至需要室内岩石破裂实验来揭示其内在机理，这也是本文将来的研究方向.

非常感谢匿名评审人对本文提出的建设性意见和建议.感谢中国地震局台网中心、四川地震局和云南地震局提供宝贵数据.