引言

由于水库蓄水或排水而引发的库区附近地震活动性（地震的频次和地震的震级）的明显增强，称为水库诱发地震，或简称水库地震（陈颙，2009）。自二十世纪六十年代全球接连发生M6以上水库诱发强震后，水库诱发地震在全世界范围内引起了广泛关注。广西龙滩水库横跨广西天峨县和贵州罗甸县，地处云贵高原南缘向桂西山区与丘陵过渡的斜坡地带，坝高220 m，正常蓄水位为400 m，总库容为273亿m³，是典型的高山峡谷高坝大库容水电站（周斌等，2014；阎春恒等，2015）。该水库自2006年9月蓄水以来，库区地震活动明显增多，截止到2019年底，已经记录到近5 000余次地震（图1），其中M_L＞4.0地震3次，对库区周边人们的生产生活造成较大影响。库区断裂发育（图1），九条主要断层中除罗甸—望谟断裂（F₁）、达恒—达良断裂（F₅）、望谟—逻西断裂（F₇）外，其余均在早—中更新世有过不同程度的活动（周斌等，2014；阎春恒等，2015）。在库区的出露地层中（图1），泥盆系—二叠系地层绝大部分为碳酸盐岩建造，岩溶、裂隙发育程度较高，渗透性较强（阎春恒等，2015；叶庆东等，2018）。作为国内最重要的水利枢纽工程之一，龙滩水库在带来巨大经济价值的同时，也持续吸引着学者们对库区地震地质构造条件（潘建雄，1989）、危险性评估（郭培兰等，2006）、活动性与孕震环境 （陈翰林等，2009a，b；詹艳等，2012；阎春恒等，2015，2020；叶庆东等，2018）、诱发的物理机制（刘耀炜等，2011；周斌等，2014）等方面的关注并开展相关研究，取得了许多重要的认识。

图  1  龙滩水库地质构造及2006年9月至2019年10月的地震事件分布

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

Figure  1.  Structural outlines and the distribution of seismic events from September of 2006 to October of 2019 in the 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-Ba’nan fault；F₉：Longfeng-Bala fault

下载: 全尺寸图片 幻灯片

水库地震主要受库水载荷与流体作用两方面的影响（蒋海昆等，2014）。水库蓄水或排水受降水量及生产需要的影响一般呈年度和季节性变化，载荷与流体作用响应也呈现相同的变化规律，进而导致地震活动性也呈现相同的变化趋势（Nur，Booker，1972）。基于此，近年来一些学者（Saar，Manga，2003；Christiansen et al，2005；Telesca，2010；Telesca et al，2015，2020）直接对地震活动与水位进行相关分析，获取水库地震活动相对于水位的延迟时间，并通过延迟时间估算库区的孔隙压力扩散系数。该系数是一个非常重要的物理参数，其数值大小直接关系到地震活动的机理，并影响库区地震活动的时空尺度，一直是地震和水文地质工作者关注的重点。然而，这方面的研究在国内并不多见，鉴于此，本文拟对龙滩水库2013年3月1日至2019年10月31日的地震活动与水位变化进行相关分析，计算地震活动对库区蓄水响应的延迟时间，进而根据延迟时间估算龙滩水库库区的孔隙压力扩散系数。

1.   数据

龙滩水库从2006年9月开始蓄水，水位在3个月左右上升100余m，几乎与此同时就触发了大量的地震。这些地震属于对库区蓄水作用的即时响应，受到库区自身的应力状态、蓄水量短时间内大量增加所致的加载和孔隙压力扩散等因素的影响。一般认为，库水加载对诱发地震的作用非常有限甚至是抑制地震发生的（Roeloffs，1988；蒋海昆等，2014）；蓄水量短时间内大量增加引起孔隙压力增加的影响时段通过分析误差余补函数的变化来获得，通常为几个月到几年（Talwani，Acree，1984；Roeloffs，1988；Turcotte，Schubert，2002；蒋海昆等，2014）。本文仅考虑蓄水完成后水位波动与地震活动的关系。为了削弱即时响应的影响，我们收集了龙滩水库完成蓄水6年以后，2013年3月1日至2019年10月31日共六年半的地震目录。该目录共包含516个地震事件，其N-t和M-t分布见图2，这些地震均已由广西地震台网进行了精确定位。采用最大曲率法（Wiemer，Wyss，2000）计算得到这个时段内的最小完整性震级为0.5，为此我们仅保留震级不小于0.5的事件信息。图2b还给出了水位的日均值变化，可以看出水位基本呈周年变化，并且对地震活动具有一定的调控作用，在高水位时期地震活动强度较大，反之在低水位时期地震活动相对较弱。

图  2  龙滩水库2013年3月1日至2019年10月31日地震活动和水位变化

（a） 地震活动N-t图；（b） M-t图和水位变化

Figure  2.  The seismicity and the variation of water level from March 1，2013 to October 31，2019 in the Longtan reservoir

（a） The N-t map of seismicity；（b） The M-t map of seismicity and the change of water level

下载: 全尺寸图片 幻灯片

2.   地震活动与水位变化的互相关分析

Nur和Booker （1972）在研究流体孔隙压力触发余震活动时提出，余震活动水平与孔隙压力变化成正比。对水库诱发地震而言，孔隙压力的变化表现为水库水位的变化，地震活动水平表现为一定时段内的地震数目或这些地震释放的地震矩（Saar，Manga，2003；Christiansen et al，2005），本文选择前者作为地震活动水平的表征。我们先将筛选出的地震目录整理成以天为单位的地震次数时间序列N，然后对每天的地震次数N与水位W进行标准化处理，即减去均值后再除以标准差，分别得到新的时间序列n和w （Telesca et al，2012）。类似于Saar和Manga （2003）、Telesca （2010）、Telesca等（2015）的做法，对序列n与w作互相关后，再采用替代数据来检验。假设地震发生与水位不存在关系，则n与w的互相关系数没有明显的极大值，或者极大值仅仅是因为原始数据包含噪声而产生的假象。为了检验这种假设，通过替代数据生成方法产生一系列与w具有相同功率谱但相位谱随机的时间序列w_R来“重新洗牌”，尽可能破坏掉原序列的原始动力结构（Little et al，2006）。如果地震活动与水位不相关，则n与w和n与w_R的互相关系数统计上处于同一水平，在n与w最大互相关系数对应的延迟时间附近，其数值相对于n与w_R的互相关系数没有明显优势。反之，如果地震活动与水位相关，则在n与w最大互相关系数对应的延迟时间附近，其值显著地大于n与w_R的互相关系数。

本文中采用傅里叶变换法（Saar，Manga，2003；刘鑫，鲍长春，2014）来生成替代数据，总共产生了1 000组替代时间序列，每一组w的替代数据w_R都与n作互相关，得到1 000组互相关系数。由于互相关系数的每一个时间点对应一个延迟时间，1 000组数据意味着每一个延迟时间处都有1 000个互相关系数值与之对应，在这些值中取出不小于序列中99%值的最小值（1 000个数据的序列，从小到大排列取第990个点）作为该点的置信水平。对每一个延迟时间对应的1 000个互相关系数作同样处理，得到99%置信水平的曲线。图3a中灰色曲线给出了n与w的互相关系数及基于替代数据得到的99%的置信水平，可以看出最大互相关系数出现在延迟时间37天处，最大值为0.186，尽管相关性不高，但却高于99%的置信水平。

图  3  标准化地震数目与水位的互相关

（a） 以天为单位的标准化后的地震数目与水位的互相关，其中灰色代表滤波之前的结果，黑色为滤波后的结果；（b） 以月为单位的标准化后的地震数目与水位的互相关

Figure  3.  Cross-correlation between the standardized number of earthquakes and water level

（a） Cross-correlation between standardized number of earthquakes and water level in day，where the gray line represents the results before filtering，and the black line represents the results after filtering；（b） Cross-correlation between standardized number of earthquakes and water level in month

下载: 全尺寸图片 幻灯片

与水位相比，地震活动显然随机性更大，频谱成分的差异可能是导致n与w互相关系数小的直接原因。为此，我们对时间序列n和w进行了功率谱分析，其结果见图4。可以看出：w含有比较明显的年度、半年度和季节性变化；而n除了年度变化比较明显之外，半年度及季节性变化不明显，高频成分较多。考虑到库区地震活动可能主要受水位长周期变化的影响，我们对序列w和n进行低通滤波，对于w保留周期大于一个季度的成分，而对于n仅保留周期大于150天的成分。w和n滤波后的时间序列w_s和n_s见图5a，在箭头1—3处，可以看出w_s与n_s相关性明显，地震活动的变化n_s相对于水位w_s略为滞后。此外，我们以月为单位统计地震数目（n_m）和水位数据（w_m），其中水位数据取每个月的平均值，其结果见图5b。相对于以日为单位的统计，以月为单位的统计也有低通滤波的作用。由于时间间隔的拉大和序列长度的缩短，在图5b中地震数目相对于水位变化的滞后更加明显，其中在箭头1，2，4和5处的地震活动相对于水位滞后1个月，在箭头6处滞后3个月。

图  4  标准化后的水位w （a）和地震数目n （b）的功率谱

Figure  4.  The power spectra of standardized water level w （a） and number n of earthquakes （b）

下载: 全尺寸图片 幻灯片

图  5  （a） 滤波后的标准化地震数目n_s和水位w_s；（b） 以月为单位的标准化地震数目n_m和水位w_m

Figure  5.  （a） Standardized number of earthquakes n_s and water level w_s after filtering；（b） Standardized number of earthquakes n_m and water level w_m in month

下载: 全尺寸图片 幻灯片

同处理n和w一样，我们对n_s与w_s和n_m与w_m也进行互相关计算，并用替代数据进行检验。图3a中黑色曲线为n_s与w_s的互相关系数及其99%置信水平曲线，对比滤波前的灰色曲线，可见互相关系数形态基本相同，但是滤波后的最大相关系数达到了0.67，是滤波前的3.6倍，表明水位的变化对地震活动产生的调控确实主要来自于其长周期波动。滤波前后两者的延迟时间均为37天，说明高频变化仅改变了相关系数的数值大小，掩盖不了水位长周期变化是引起库区地震活动的主要因素这一实质。图3b为n_m与w_m的互相关系数及其99%置信水平曲线，最大相关系数为0.58，对应的最大延迟时间为1个月，与n_s和w_s互相关系数的结果一致，进一步确认了延迟时间为37天是可靠的。

3.   扩散系数估算

孔隙压力的增加降低了岩石的有效正应力和岩石摩擦系数，从而使介质的抗剪强度降低。当岩石处于破裂的临界状态时，孔隙压力微小的增加都可能触发地震，因此孔隙压力扩散被认为是触发水库地震的一个主要机制（Nur，Booker，1972；Talwani，Acree，1984；Roeloffs，1996）。孔隙压力扩散由扩散方程描述，在扩散方程中，扩散系数是衡量扩散快慢的物理参数，与介质的渗透系数、流体的动力学黏度及孔隙弹性模量等参数的大小相关。地表附近的孔隙压力扩散系数可以通过相关试验获取，深部的孔隙压力扩散系数获取非常困难。水库库区的水位变化与地震活动的关系提供了一个探测库区孔隙压力扩散系数的“窗口”，是孔隙压力扩散系数原位测试的试验场（Talwani，Acree，1984；龚钢延，谢原定，1991）。

3.1   基于自由扩散模式的扩散系数估算

在基于扩散方程估算孔隙压力扩散系数的已有研究中，大多忽略孔隙压力扩散与应力的耦合（Talwani，Acree，1984；Saar，Manga，2003；Telesca，2010；Telesca et al，2015），仅考虑孔隙压力的自由扩散。在此情形下仅考虑垂向扩散，一维孔隙压力扩散满足的方程为：

式中，P为孔隙压力，D为孔隙压力扩散系数，h为深度，t为时间。在水库坝基处孔隙压力主要受库水周年性变化调控，则边界条件为P（t，h＝0）＝P₀exp（iωt），其中P₀为坝基处（h=0）由水位变化引起的压力周年性变化的振幅，ω＝2π/T，T为周期，T＝1 a，i为虚数单位。式（1）的解为（Turcotte，Schubert，2002）：

式（2）表明扩散方程在周期性的边界条件下，其解可表示为一个振幅衰减的平面波，衰减系数为$ \sqrt {{\omega / {2D}}}$，波速为$ \sqrt {2\omega D}$。指数项$\exp （ - h\sqrt {{\omega \mathord{\left/ {\vphantom {\omega {2D}}} \right. } {2D}}} ） $描述了孔隙压力振幅随深度h的衰减，而谐波项$\exp ［ {\rm{i}}\omega （ t －{h / {\sqrt {2\omega D} }} ） ］ $描述了孔隙压力作为时间t与深度h的函数周期性变化的规律，两者均与扩散系数D相关。在某一深度岩石已处于破裂的临界应力状态，孔隙压力扰动前锋传播至此触发地震，则应满足$t －h / {\sqrt {2\omega D} ＝ 0}$，整理得

式中，t＝37 d为延迟时间，h取趋肤深度（Saar，Manga，2003）或者地震分布的优势深度（Telesca，2010；Telesca et al，2015），本文取后者。图6给出了库区地震的深度分布图，可见深度主要集中在5—7 km，均值μ为5.94 km，标准偏差σ为1.44 km。将ω，t和h的值代入式（3），可得D＝（8.66±4.11） m²/s，其中，不确定度由地震深度分布的一倍标准偏差确定。

图  6  地震的深度-频次分布

Figure  6.  Depth-frequency distribution of the earthquakes

下载: 全尺寸图片 幻灯片

在自由扩散模式下，基于扩散方程在周期性边界条件下的解和库区地震平均深度得到的龙滩水库孔隙压力扩散系数D＝（8.66±4.11） m²/s，处于0.1—10 m²/s的合理范围内（Talwani，Acree，1984）。然而，该结果与刘耀炜等（2011）基于水压致裂中使用的单点触发模型采用抛物线包络拟合方法所得到的龙滩库区蓄水初期的孔隙压力扩散系数（0.04—0.65 m²/s）差异较大。我们认为差异可能源于以下三方面。其一，根据地质资料（广西壮族自治区地质局，1972；潘建雄，1989），龙滩库区以泥盆系—二叠系碳酸盐岩建造为主，岩溶、裂隙发育程度较高，渗透性较好，受地震定位精度的影响，加之抛物线包络拟合具有一定的主观性，刘耀炜等（2011）估算的值可能偏低。其二，刘耀炜等（2011）是在水位快速上升阶段，假定孔隙压力边界条件为阶跃函数得到的结果；而本文结果是在蓄水完成六年以后，基于水位总体稳定仅有周年扰动变化得到的结果，不同时段库区的介质状态不一样，因而扩散系数也不一样。其三，自由扩散忽略孔隙压力与应力的耦合，模式过于简单。

3.2   基于耦合模式的扩散系数估算

虽然自由扩散模式下得到的孔隙压力扩散系数在合理范围内，但物理上孔隙压力与应力是耦合的，自由扩散模式仍过于简化。根据Biot （1941）的固结理论，Roeloffs （1988）得到了周期性边界条件下一维孔隙压力与应力耦合模式扩散方程：

式中，$\alpha ＝ {{B （ 1 ＋{\nu _{\rm{u}}} ） } / { ［ {3 （ 1 － {\nu _{\rm{u}}} ） } ］ }}$，B为Skempton系数，ν_u为无排水条件下孔隙介质的泊松比。式（4）的解为

式（5）表明深度h处的孔隙压力除了与扩散率D和角频率ω有关，还与系数α有关，α为解耦模式与无排水响应模式的加权平均。当α＝0时，式（5）与式（2）等价，对应于孔隙压力与应力解耦的模式；当α＝1时，对应于无排水响应模式。

式（5）中，对于给定$ \sqrt{\omega /2D}h$，相位仅与α相关，可以通过相位的变化分析时间的超前与滞后，相位正代表时间的超前，相位负对应着时间的滞后，最小值对应着最大延迟时间。图7给出了α＝0，0.05，0.30，0.70，1.00时，孔隙压力及相位与$ \sqrt{\omega /2D}h$ 的关系，可见在相同的$ \sqrt{\omega /2D}h$ 情况下，α越大，孔隙压力越大，延迟时间越小。已知延迟时间为37 d，通过搜索相位的最小值求得对应的α＝0.219，对应的$ \sqrt{\omega /2D}h$＝1.43，ω=2π/T，周期T取1 a，h取地震的平均深度，可得扩散系数D＝（1.72±0.82） m²/s。可见，是否考虑孔隙压力扩散与应力耦合，估算出来的扩散系数差异较大。这种差异可以由延迟时间直观地判断，当延迟时间较大（接近180 d）时，是否考虑孔隙压力扩散与应力的耦合对估算D影响不大；当延迟时间较小（远小于180 d）时，是否考虑孔隙压力与应力的耦合则可能对估算D有一定影响。由$ \sqrt{\omega /2D}h$＝1.43得到对应的$ P/{P}_{0}＝0.307 $，根据图1，龙滩水库水位年变化在330—370 m之间，水位变化Δz＝40 m，则${P}_{0}＝\rho g\Delta {\textit{z}}/2$，其中ρ为水的密度，取1 000 kg/m³，g为重力加速度，取9.8 N/kg，则P＝60 kPa，孔隙压力的年变化为120 kPa。假定岩石强度的降低完全由孔隙压力的变化引起，120 kPa的孔隙压力增加会导致近100 kPa的库仑应力变化，在临界状态下足以触发地震（Hardebeck et al，1998；Stein，1999）。

图  7  α取不同值时孔隙压力幅值（a）及相位（b）随$ \sqrt{\omega /2D}h $ 变化

Figure  7.  Variation of the magnitude （a） and phase （b） of pore pressure as a function of $ \sqrt{\omega /2D}h $ for different α

下载: 全尺寸图片 幻灯片

4.   讨论与结论

通过对龙滩水库的水位与地震数目的互相关分析，得到了地震活动相对于水位的延迟时间为37天。与一些蓄水多年的水库（Talwani，Acree，1984；Ferreira et al，1995；Telesca，2010）相比，37天的延迟时间是一个较小的值，由于孔隙压力扩散范围正比于延迟时间（Talwani，Acree，1984；龚钢延，谢原定，1991），这暗示库水的影响主要局限于库区较小的范围内，这与库区地震精定位的结果（阎春恒等，2015，2020）一致。由于水库地震震源深度较浅，通常M3以上地震均会有震感，引发社会舆论，对人们的生活生产造成影响，地震活动强度相对于水位峰值37天的延迟时间为相关部门制定震情跟踪、舆情应对等预案提供了时间节点信息。

在不考虑孔隙压力扩散与应力耦合的情况下，基于周期性边界条件下孔隙压力扩散方程的解和库区地震的深度分布，利用37天的延迟时间得到扩散系数D＝（8.66±4.11） m²/s；在考虑孔隙压力扩散与应力耦合的情况下，得到扩散系数D＝（1.72±0.82） m²/s，两者均处于合理的范围内。数值上，基于孔隙压力与应力耦合模式得到的扩散系数与刘耀炜等（2011）的结果比较接近，也与近年来部分学者（Saar，Manga，2003；Telesca，2010；Telesca et al，2015）利用地震活动与水位变化的互相关分析得到的其它库区的扩散系数接近，但由于这些研究大多采用了自由扩散的假定，因此无法说明D＝（1.72±0.82） m²/s比D＝（8.66±4.11） m²/s更为合理。从物理上以及近年来对断裂带的孔隙压力扩散系数原位测量的结果（Doan et al，2006；Xue et al，2016）来看支持孔隙压力扩散系数取较小的值，这可能说明对龙滩水库而言，孔隙压力扩散与应力耦合是现阶段诱发地震的主要因素。