引言

广西龙滩水库横跨广西天峨县和贵州罗甸县，地处云贵高原南缘向桂西山区与丘陵过渡的斜坡地带，由南盘江、蒙江、曹渡河、布柳河和红水河等多条河流组成（图1），年平均发电量为1 870亿kW·h，是国内仅次于长江三峡及金沙江溪洛渡的特大型水电站（周斌等，2014；阎春恒等，2016）。该水库自2006年9月蓄水以来，库区地震活动明显增多，截止到2016年底，已记录到近5 000次地震，其中M4以上地震3次，对库区的生产生活造成了较大影响，引起了社会的广泛关注，也促使众多地震研究人员从地球物理、地震地质等不同角度对龙滩水库进行了研究。陈翰林等（2009a，b）率先从精定位与震源机制方面分析这些水库诱发地震的特征，基于精定位后地震的空间分布表明库水的渗透作用与库区的断层共同控制了这些地震的发生，基于震源机制的分析结果则表明载荷作用所引起的剪应力变化不是龙滩水库诱发地震的主要因素，而蓄水所导致的库区介质孔隙压力的改变和库水的渗透润滑弱化作用则可能是这些地震的成因。华卫等（2009，2012）分析了库区地震的衰减特征及震源力学参数，结果表明库区的地震衰减特征及力学参数均与构造地震存在差异。刘耀炜等（2011）分析了龙滩水库诱发地震的孔隙压力扩散特征，认为孔隙压扩散系数与地震活动呈正相关。Zhou等（2011，2012）通过三维速度结构和Q值衰减成像，讨论了库区水渗透深度的空间差异。詹艳等（2012）利用大地电磁的结果分析了库区岩性与电阻率的关系，并推测了库水的渗透深度。周斌等（2014）基于库区地震精定位、震源机制、地质构造、深部速度结构等资料，开展三维孔隙弹性有限元数值模拟研究，讨论了库区诱发地震的时空演化与库水加卸载、渗透过程的动态响应关系及其可能的成因机制。阎春恒等（2015，2016）采用了更为丰富的数字地震波形资料，通过地震精定位和震源机制分析，得到了与陈翰林等（2009a，b）类似的结果。体波速度和衰减成像（Zhou et al，2011，2012）、电阻率成像（詹艳等，2012）、小震精定位的优势深度（陈翰林等，2009a；阎春恒等，2016）的结果均显示龙滩水库的库水渗透深度达到7 km以上，震源机制（陈翰林等，2009b；阎春恒等，2015）、库区孔隙压力扩散特征（刘耀炜等，2011）及数值模拟（周斌等，2014）则揭示了库区不同区域的地震发生受到断裂构造、流体加载和润滑等的差异性作用。

图  1  龙滩水库地质构造

F₁：罗甸—望谟断裂；F₂：高圩—八茂断裂；F₃：凤亭—下老断裂；F₄：马耳—拉浪断裂；F₅：达恒—达良断裂；F₆：党明—桂花断裂；F₇：望谟—逻西断裂；F₈：长里—八南断裂；F₉：龙凤—八腊断裂，下同

Figure  1.  Structural outlines of Longtan reservoir

F₁：Luodian-Wangmo fault；F₂：Gaoyu-Bamao fault；F₃：Fengting-Xialao fault；F₄：Maer-Lalang fault； F₅：Daheng-Daliang fault；F₆：Dangming-Guihua fault；F₇：Wangmo-Luoxi fault；F₈：Changli-Banan fault；F₉：Longfeng-Bala fault. The same below

下载: 全尺寸图片 幻灯片

根据实验研究，刻画地震波衰减的Q值与岩石内部裂纹密度、分布、构造及孔隙流体的相互作用等微观性质密切相关（陈颙等，2009），因此Q值成像被认为是反映库水渗透最有效的手段（王慧琳等，2012； Zhou et al，2011；吴海波等，2016）。理论上，Q值是空间三维坐标的函数，但三维的Q值成像依赖于研究区域的三维速度结构，而速度结构的误差会传递给Q值，造成Q值成像结果存在较大的误差。为此，很多研究（苏有锦等，2006；汪素云等，2007，2008；周龙泉等，2009；裴顺平等，2010）采用二维成像研究Q值的横向非均匀性。二维Q值成像可视为不同深度的Q值加权平均压缩至平面上，能很好地反映Q值的横向非均匀性。考虑到地壳内地震存在优势深度，二维Q值成像事实上主要反映了优势深度以上地壳介质的平均Q值分布特征。马宏生等（2007）和周龙泉等（2009）进一步假设Q值与频率无关，虽然这样会影响Q值的大小，但是不会影响Q值的分布特征。近震波形资料主要包括P波和S波两部分，Q值成像也对应有Q_P和Q_S成像，由于S波对流体作用更为敏感，且振幅相对于P波更大，信噪比更高，因此通常利用Q_S成像来研究库区的衰减特征与库水渗透、库区地质构造之间的关系。鉴于此，本文拟基于龙滩水库遥测台网2006年至2016年近4 000次M_L≥0地震的波形数据，采用基于波形互相关的双差定位方法（Waldhauser，Ellsworth，2000）确定这些地震的精确位置，然后使用周龙泉等（2009）的二维成像方法反演得到龙滩水库库区Q_S值的二维分布图像。与Zhou等（2011）的三维Q_S成像不同的是，本文试图利用更长时间的观测资料，来勾画龙滩水库库区Q_S的横向非均匀性，以期讨论地质构造、岩层岩性、地震活动性与Q_S分布的具体关系。

1.   地质背景与监测台网

龙滩水库（106.6°E—107.3°E，24.8°N—25.5°N）处于云贵高原南缘向桂西北山区与丘陵过渡的斜坡地带，高山峡谷纵横，属于典型的高山峡谷型水电站。大地构造上，该库区位于广西山字形构造西翼（郭培兰等，2006），褶皱强烈，断裂发育。据野外地质考察，库区内主要断裂有9条（图1），除罗甸—望谟断裂（F₁）、达恒—达良断裂（F₅）和望谟—逻西断裂（F₇）为前第四纪活动断裂外，其余断裂在早—中更新世有过不同程度的活动，断裂破碎带发育，裂隙和岩溶沿断裂带呈线状分布（郭培兰等，2006；阎春恒等，2016）。库区出露的地层主要有泥盆系（D）、石炭系（C）、二叠系（P）及三叠系（T）等（图1），其中三叠系约占70%。库区内碳酸盐岩发育，整体透水性较好，但同一种地层东西部的岩性存在区域差异，同一系内上、中、下统的岩性也存在差异，如对于二叠系上统的岩层，东部主要为碎屑岩、硅质岩，西部为碳酸盐岩，而下统东西部均以灰岩为主（郭培兰等，2006）。

龙滩水库数字遥测台网自2006年3月开始正式运行，初期有10个监测台站（图2），包括1套宽频带地震仪（天峨台）和9套短周期地震仪以及2个信号中继站和1个台网中心。为进一步提高监测能力，2012年增加杠里和向阳两个台，配备FSS-3B型短周期地震仪。关于台网观测系统的设计、仪器型号、监测能力等详见相关文献（姚宏等，2008a，b；杨超英，姚宏，2010）。

图  2  重定位后的地震分布及本文所使用的台站分布

Figure  2.  Distributions of relocated earthquakes and the stations used in this study

下载: 全尺寸图片 幻灯片

2.   数据

本文收集了龙滩水库库区2006年9月至2013年5月底精定位后的3 082次M_L≥0地震的波形资料（阎春恒等，2016）和2013年6月至2016年12底经龙滩水库地震监测台网中心初步定位后约330次M_L≥0地震的波形资料。对2013年6月至2016年12月的地震资料，参照阎春恒等（2016）的参数选择，同样采用结合波形互相关的双差定位方法进行重定位，最终得到300次地震精定位的结果。图2给出了精定位后的3 382次地震的分布，由于结果中阎春恒等（2016）的结果占绝大多数，因此整体分布特征与阎春恒等（2016）相同。可以看出，地震丛集特征明显，总共可分为6个地震丛，即罗妥地震丛（ Ⅰ ）、八茂地震丛（ Ⅱ ）、拉浪地震丛（ Ⅲ ）、坝首地震丛（ Ⅳ ）、天峨地震丛（ Ⅴ ）、布柳河地震丛（ Ⅵ ），其中坝首地震丛（ Ⅳ ）和天峨地震丛（ Ⅴ ）的优势深度略浅，处于6—7 km，其余地震丛的优势深度大致为7—9 km，5 km以内的地震数目多于10 km以上的地震，与中国大陆水库地震的整体统计特征（蒋海昆等，2014）一致。库区地震活动与坝前水位变化存在很好的相关性（图3），下闸蓄水后，库区的地震活动性明显增强，至2008年底地震活动强度达到顶峰，之后自2009年开始逐渐衰减，2012年以后，库区每年发生的地震为100次左右。

图  3  龙滩水库库区地震活动M-t图及水位变化

Figure  3.  M-t map of seismicity and the change of water level in Longtan reservoir

下载: 全尺寸图片 幻灯片

3.   方法

3.1   路径衰减计算

将精定位后的震源参数代替震相报告中的震源参数，挑选出对应的事件波形，并去掉仪器响应；然后基于该震相报告中的P波、S波到时，采用自动截取技术（黄玉龙等，2003；苏有锦等，2006；华卫，2007），取S波窗长t_win＝1.084t_S-P＋0.376 （t_S-P为P，S波到时差），逐一截取南北、东西两个分量的S波段和噪声段；再分别对S波段记录和噪声记录作傅里叶变换，同时除以2π f （ f为频率），得到S波位移谱u（ f ）和噪声位移谱N（ f ），二者均取东西、南北两个分量平方和的算术平方根。定义信噪比为S波振幅谱与噪声振幅谱之比的最小值，若信噪比大于3则保留对应的S波振幅谱；最后，从所选取的S波谱中减去对应的噪声谱，进一步获得“干净”的S波振幅谱（图4）。

图  4  仁顶台记录到的2007年12月10日5时44分贵州罗甸M_L1.0地震的两个水平分量波形（a）及其合成位移谱和噪声谱 （b）. 图（b）中观测谱和拟合谱均未作几何扩散校正

Figure  4.  Two horizontal seismic waveforms of a M_L1.0 earthquake recorded by Rending station （a），and the synthesize displacement spectra and noise spectra （b）

The earthquake occurred in Luodian city，Guizhou Province，at 5:44 on December 10，2007. In Fig. （b） the geometrical spreading corrections were conducted neither for the observed spectrum nor the fitted spectrum of S wave

下载: 全尺寸图片 幻灯片

在Brune （1970）震源模型的假设下，去掉仪器响应、几何扩散影响的S波位移谱具有如下的表达形式：

式中，A（ f ）为未去掉非弹性衰减的震源谱，f_c 为震源谱拐角频率，r为震源距，t为S波走时，Ω₀为零频极限。本文中，几何扩散采用Atkinson和Mereu （1992）的三段扩散模型；Q_S为地震射线路径上的整体衰减，引入衰减算子t^*＝t/Q_S，可得

由于Q_S是频率的函数，因此上式中t*也是频率的函数。已有研究表明，10 Hz以下Q_S强烈依赖于频率，而10 Hz以上，Q_S与频率的相关性很弱（Adams，Abercrombie，1998），由于龙滩库区的地震震级普遍较小，10 Hz以上的频谱成分丰富，因此本文忽略频率的影响。对于一次地震，零频极限Ω₀和拐角频率f_c是唯一的，而t*反映了震源至接收台站间路径的衰减，其个数等于记录到该地震并且波形满足信噪比要求的台站个数。如果该地震被n个台站记录到，则需确定零频极限Ω₀、拐角频率f_c和n个t*共n＋2个参数。假定有m个频点，则建立目标函数$\sum\limits_{j {\text{＝}}1}^m {\sum\limits_{i {\text{＝}}1}^n {\lg \left| {A_i^{\rm cal}({f_j}) {\text{－}} A_i^{\rm obs}({f_j})} \right|} } $，并采用遗传算法来拟合这n＋2个参数。本文所采用的震源谱的频率范围为1—30 Hz，根据陆丽娟等（2015）及叶庆东等（2017）对震级与拐角频率关系的研究，该频率范围可以保证拐角频率前后的频点都得以使用，一定程度上避免了因使用的频点过少或频点分布不合理而对拟合结果产生影响。

3.2   二维Q_S反演成像

利用遗传算法对式（2）进行拟合求得的衰减t*为地震射线路径上的整体衰减，是Q_S与S波速度v_S的泛函，在射线理论下可以表示为

显然，Q_S和v_S两个未知量必须知道其一才能反演另外一个量。二维反演中，假定地壳内S波速度均匀，则射线是直线，v_S由走时-震中距直线拟合的斜率确定。根据图5a，最终拟合得到的v_S为3.46 km/s，则式（3）可变为

图  5  S波走时t、衰减算子t^*与震源距Δ的关系

（a） S波走时t与震源距Δ的关系；（b） 衰减算子t*与震源距Δ的关系；（c） 根据一倍均方根准则筛选出的衰减算子t*与震源距Δ的关系

Figure  5.  The travel times t and attenuation operator t* of S wave versus hypocentral distances Δ

（a） The relationship between the travel times t of S wave and hypocentral distances Δ；（b） The relationship between the attenuation operator t* and hypocentral distances Δ；（c） The selected t* based on one-root-mean-square criterion varied with hypocentral distances Δ

下载: 全尺寸图片 幻灯片

式（4）与走时反演的形式相同，由于本文仅考虑二维反演，因此上式与面波频散二维反演、莫霍面首波速度反演类似。当射线较多时，由式（4）组成的线性方程组所包含的方程数目巨大，求解十分费时。本文采用周龙泉等（2009）编写的Matlab反演程序进行反演，该程序同时利用了Matlab存储稀疏矩阵和矩阵快速求逆的优势，能大大地节约计算时间。

4.   结果

根据以上方法，本文共获得14 815个信噪比较高的t^*。由于Q_S值的测量误差较大，还需进一步对数据进行筛选，具体做法是：对观测值t^*与震源距进行直线拟合（图5b），其斜率为区域的平均Q_S值，然后计算t^*观测值与拟合直线之间残差的一倍均方根值，当观测值与拟合值之差小于一倍均方根值时保留，否则舍弃（图5c）。我们得到的Q_S平均值为240，低于周龙泉等（2009）得到的云南地区Q_S平均值400，也低于中国大陆不同块体的平均值（汪素云等，2007，2008），但高于王惠琳等（2012）得到的紫坪铺水库库区和吴海波等（2016）得到的长江三峡库区的Q_S平均值（前者为90，后者为180）。经过进一步的质量控制之后，共获得9 166个t^*数据用于最终成像。由于坪上寨台离6个地震丛均较远，波形记录噪声较大，最后未得到该台的有效记录。

经过10次迭代后，t*残差均方根由最初的0.021 6降至0.017 4，反演前后残差随震中距的变化如图6所示。射线分布和检测板测试是评价反演结果可靠性的常用方法。图7a给出了射线分布，与地震分布相对应，射线整体上呈现东密西疏的特点。图7b给出了0.05°×0.05°的检测板测试结果，由于射线数量多，射线覆盖区基本能达到0.05°×0.05°的分辨。在射线外侧的边缘区域，由于拖尾效应（He et al，2017），造成了一些虚假分辨，为此我们裁掉了没有射线覆盖的外围区域，最终的成像结果见图7c。可以看出，除库区西南部存在明显的Q_S高值区，Q_S低值区围绕着库区大体呈环形分布，这与王慧琳等（2012）得到的紫坪铺水库库区Q_S的分布特征类似；在Q_S低值区附近，大多为断裂与蓄水区的交汇处，其中Q_S极低值出现在库区西北区罗甸—望谟断裂（F₁）、高圩—八茂断裂（F₂）和凤亭—下老断裂（F₃）所包围的区域内，这与Zhou等（2011）极低值出现在坝区不同，造成这种差异的原因是本文的二维成像对横向非均匀性敏感，但只能反映地震震源优势深度以上地壳介质的Q_S平均值。结合图1可以看出，西北侧大面积的Q_S低值区主要分布在二叠系上统（P₂）碳酸盐岩地区，东侧的Q_S低值区分布在三叠系上统（T₂）与二叠系上统（P₂）的硅质灰岩、硅质页岩夹砂岩、页岩碎屑岩地区。由图7c还可以看出，地震大多分布在Q_S高、低值过渡区域。

图  6  反演前（a）、后（b）t^*残差δt^*随震中距Δ的变化

Figure  6.  The residuals of attenuation δt^* varied with hypocentral distances Δ before （a） and after （b） inversion

下载: 全尺寸图片 幻灯片

图  7  龙滩水库库区Q_S成像结果

（a） 射线分布；（b） 检测板测试结果；（c） Q_S二维分布

Figure  7.  Results of Q_S tomography in Longtan reservoir

（a） The distribution of rays；（b） The result of checkboard test；（c） Two-dimensional distributions of Q_S

下载: 全尺寸图片 幻灯片

5.   讨论与结论

本文反演得到了广西龙滩水库库区的高分辨率二维Q_S图像，结合库区的构造、岩性、地震空间分布及Q_S二维成像方法的特征，得到以下结论：

1） 龙滩水库库区的Q_S值横向不均匀变化明显，Q_S低值区围绕着库区近似呈环形分布，但西南部除外。在Q_S低值分布区附近，大多为河流与断裂带的交汇处，表明Q_S低值分布受水和断裂带的影响。一般认为，Q值非均匀性受固有衰减和散射衰减两种机制控制（Hoshiba et al，2001；何静等，2017）。固有衰减是指地震波在传播过程中，介质微粒运动克服微粒之间的内摩擦力，将地震波能量转化成热能，从而导致地震波能量的衰减；散射衰减是指由于介质的不均匀性，地震波遇到散射体而发生散射，进而导致能量重新分配所产生的衰减。虽然两种机制并存，很难通过观测将二者区分开来，但理论上后者更倾向于解释尾波Q值衰减（Aki，Chouet，1975）。根据库区Q_S低值的空间分布与断裂、河流分布的相关性可以初步推断，库水可能沿着断裂上的岩石破碎带及节理、裂隙发育地区向下渗透，使岩石孔隙中充满流体，内摩擦增大，地震波剧烈衰减，从而使Q_S值大幅下降。在库区西南侧，由于没有大的河流穿过，表现为Q_S高值。

2） 从岩性上来看，研究区域西北部面积最大的Q_S低值区分布在透水性很强的二叠系上统（P₂）碳酸盐岩区，东部的Q_S低值区主要分布在三叠系上统（T₂）砂岩和二叠系上统（P₂）的硅质灰岩、硅质页岩夹砂岩和页岩碎屑岩区。尽管岩石地层学结果表明以硅质灰岩、硅质页岩夹砂岩的透水性相对较弱，但由于库区晚二叠系龙滩煤系中的煤层为封存的古岩溶水提供了炭分，增强了对硅酸盐的溶解能力（潘建雄，1989），使硅质灰岩、硅质页岩夹砂岩的含水性和透水性增强，因此与之对应的东部表现为Q_S低值。布柳河地震丛（Ⅵ）所在的区域虽然有断层穿过，而且岩性以含水性和透水性均较好的二叠纪下统（P₁）灰岩为主，但是，震源机制解结果（阎春恒等，2015）显示该区域以挤压逆冲为主，且倾角较小，因此该区在蓄水初期可能库水的渗透润滑作用明显，随着蓄水时间的增长，库水的重力加载作用反而可能会抑制库水的渗透，因此表现为Q_S高值分布，这与岩石物理实验（陈颙等，2009）和数值模拟（周斌等，2014）的结果是一致的。

3） 除布柳河地震丛（Ⅵ）之外，其余5个地震丛基本都分布在Q_S高、低值的过渡区域。通常，均匀介质内不易积累能量，而在Q_S高、低值过渡区域这种“软”、“硬”介质的交界处，容易积累应变能，进而孕育地震。当Q_S高、低值过渡区域受到水的加载作用和润滑作用（Brodsky，Kanamori，2001）时，发震断层的抗剪强度降低，容易产生滑动，从而诱发地震。

4） 二维Q_S成像除了受观测数据的影响之外，理论上还受两个因素的影响：一是二维成像将不同深度的Q_S值加权平均压缩至平面，在某些区域会因为浅部和深部高低相反的Q_S值分布相互抵消而无法显示细节；二是反演中假定S波速度为恒定值，而实际中上地壳S波速度扰动可达到±10%，由于S波速度与Q_S的耦合，低于平均S波速度的区域可能会被Q_S高值部分补偿，反之亦然。基于以上两个原因，本文的二维Q_S成像结果与Zhou等（2011）的三维成像结果的各个深度剖面均存在一定差异；同样也由于这两方面的局限，本文的结果无法像Zhou等（2011，2012）及詹艳等（2012）的结果那样准确地反映出库水最大渗透深度及其所在区域。考虑到Q_S扰动通常较大（本文中最大可达50%），我们认为本文二维反演所揭示的Q_S横向非均匀性对认识龙滩水库库区地质构造、岩性与地震活动性之间的关系具有一定的价值。

感谢中国地震台网中心周龙泉研究员提供计算程序，感谢审稿人的细致审阅和宝贵意见。