利用背景噪声研究紫坪铺水库水位加卸载对库区地下介质波速的影响

安艳茹 张晓东

1）中国地震局地震预测研究所，北京市复兴路63号 100036

2）中国地震台网中心，北京市西城区三里河南横街5号 100045

0 引言

中国是世界上水库数量最多的国家，截至2011年底，全国已建成水库97246座，总库容8104.1亿m3，在建水库756座，总库容 1219.0亿 m3（孙振刚等，2013）。中国也是世界上水库地震多发地之一，尤其是近些年在地震多发的川滇地区建设了大量高坝、大库容的梯级电站，水库地震时有发生，大型水库的水库地震问题亦引起了地震学者的高度重视。汶川地震后，由于紫坪铺水库与其特殊的时空关系，更是引发广泛关注（雷兴林等，2008；Shemin et al，2009；马文涛等，2011；Deng et al，2010；卢显等，2010）。雷兴林等（2008）认为，紫坪铺水库在蓄水过程中对其地下的龙门山中央断层和山前断层有明显的作用；Shemin等（2009）认为蓄水后的库仑应力增量足以引发地震；马文涛等（2011）认为汶川地震与紫坪铺水库在时空上均有一定的联系；Deng等（2010）认为紫坪铺水库周边的M-t图像没有增强，库水不能渗透到汶川地震深度14km的震源处；卢显等（2010）认为都江堰震群与紫坪铺水库无关。无论如何，水库的蓄水对地下介质产生影响都是个不争的事实，而利用数字地震资料研究库区介质和应力场变化则是十分有效的方法。目前有关库区水位变化与库区地震的关系、库区介质变化与库区地震的关系等方面的研究较多（周斌等，2010；王惠琳等，2012a、2012b；卢显等，2013），但有关库区水位变化与库区介质变化的研究还相对较少。本文基于紫坪铺水库数字地震台网的连续波形数据，利用噪声互相关技术研究水位的加卸载和渗透过程对库区地下介质波速的影响，以求深化对库区地下结构、状态和动力学过程的认识。

1 紫坪铺水库及地震台网概况

紫坪铺水库位于四川省都江堰市区西北方向9km处的岷江上游麻溪乡，距成都市区约60km，距汶川8.0级大地震微观震中最近处仅6km左右。紫坪铺水库于2001年3月29日正式动工兴建，2006年12月竣工，最大坝高156m，总库容11.12×108m3。本文利用四川省地震局水库地震研究所提供的水位数据，绘制了紫坪铺水库水位变化图（图1）。由图1可见，2005年9月30日下闸蓄水，数天内水位迅速升高至死水位817m以上，蓄水后最高水位为875.18m。2005年9月30日～2007年底，紫坪铺水库共经历了3次大规模蓄水、2次泄水过程。2005年12月5日水位达835.91m，2006年10月14日水位上升到最高值875.18m，2007年12月12日水位再次上升到873.39m。

图1 紫坪铺水库水位变化

紫坪铺水库位于青藏高原东缘的龙门山造山带的中段。库区及邻区被主干断裂分隔，可划分为茂汶韧性剪切带、中央推覆构造带、龙门山前缘拆离带、前陆扩展变形带和川西前陆盆地等5个地质构造单元，它们在物质组成、构造层次和变形样式等方面具有明显的差异。紫坪铺水库库体坐落于龙门山前缘拆离带内，邻近中央推覆构造带和前陆扩展变形带（周斌等，2010）。库区地下水分为碎屑岩类裂隙孔隙水和碳酸盐岩岩溶水等2大类（图2），如图2所示，库区西南主要岩性为碳酸盐岩，东北主要岩性为碎屑岩和部分碳酸盐岩。

紫坪铺水库数字遥测地震台网于2004年8月16日正式采集地震信息，2005年6月27日通过验收（胡先明等，2006）。台网包括7个数字遥测地震台，使用短周期地震仪，频带宽度1～40Hz。此外，该台网共享四川区域数字台网中YZP（油榨坪）台的数据。8个地震台大致均匀地分布在库区周边，平均台间距约为14km，台站分布如图3所示，台站及仪器参数见表1。

图2 库坝区地下水类型（据王云基（2001））

图3 紫坪铺水库台网台站分布

表1 紫坪铺水库台网台站及仪器参数表

2 观测资料

本研究使用紫坪铺台网及YZP台垂直分量的连续波形数据，时间范围2005年1月1日～2008年1月1日。图4为各台站数据连续记录情况，可用数据的天数为1003天，最长连续空缺时间为27天。

图4 各台数据连续情况

3 理论基础与数据处理

利用背景噪声测量地下介质波速连续变化的主要技术路径是：比较经验格林函数（CCFs）与参考格林函数（REF）。其中，CCFs表征的是一段时间内地下介质的状态，而 REF代表的是地下介质的背景状态。计算过程分为4步：①计算台站对间不同时间段的噪声互相关函数，即经验格林函数CCFs；②获取参考格林函数REF，并计算每个CCFs与REF的走时延迟Δt；③对不同台站对间不同时间错动的走时延迟进行平均，并使用一个相对速度均匀一致变化（Δv/v=常数）的简单模型（Lecocq et al，2014）来表征库区地下介质波速的变化。

本文使用MSNoise软件（Lecocq et al，2014）处理数据。该软件已应用于奥克兰火山地区（Kasper，2014）及地震前后（Berkeley Seismological Laboratory，2014）的分析研究中，并得到了稳定可靠的结果。

3.1 噪声互相关提取经验格林函数CCFs

3.1.1 单台数据预处理

本文研究的是库区地下介质微小的波速变化，因此能否得到稳定可靠的CCFs非常关键。首先，我们必须通过数据预处理，获得较为纯净的背景噪声，以保证CCFs的质量。将每个台站的连续波形数据合并成每天一个文件，然后进行预处理，包括：去平均、尖灭处理，高通、低通滤波，降采样，时间域归一化，谱白化处理等。通常，波形数据总会存在一个非零的均值，这会影响到数据的分析；为了使数据标准化，适用于各种标准公式，必须在数据分析前去平均值（郑治真，1979）。另外，在对数据进行谱域操作，如FFT、滤波时，若数据的两端不为零，则会出现谱域假象，因而实际数据经常需要做尖灭处理，以使数据两端在短时间窗内逐渐变成零值。本文研究的是库区浅层介质波速的变化，关注的是相对高频的信息，所以我们选择了0.1～2.0Hz的滤波参数。为了减少计算量和存储量，将数据降采样到10Hz。为了减少地震信号、仪器异常信号和台站附近非稳定噪声源的影响，我们对数据进行了时间域归一化处理（Tukey，1962）。为了使信号在不同频率上的能量更为均衡，我们对数据进行了谱白化处理。

3.1.2 互相关计算提取经验格林函数

利用背景噪声提取地下结构信息设想的探索最初可追溯到20世纪50、60年代，Aki（1957）认为可以用背景噪声提取地下结构中的面波频散性质；Claerbout（1968）提出可以利用背景噪声恢复一维层状介质的反射响应。类似想法的首次成功应用是在太阳地震学研究中，即通过对太阳表面噪声进行互相关计算，成功地提取出了声波时距曲线（Duvall et al，1993）。随后，噪声互相关技术在超声学领域获得重大进展。Weaver等（2001）在声学上给出了随机波场的互相关特性，即通过从发散的或随机的波场中提取格林函数来计算地球的弹性响应。证明该特性的实例就是基于一个弹性体中有发散场的模型

式中，x为点的位置；t为时间；i为虚数单位；an为模态激发函数；un和ωn为地球的本征方程和本征频率。发散场的一个重要性质就是模态振幅是不相关的随机变量

式中，δnm为克罗内克函数；F为频谱能量密度。因式（2）的交叉项在取平均时消失，所以场中x、y点的相互关系简化为

式中，τ为时间错动。比较发现，式（3）与x、y点间真实的格林函数相比，仅差一个振幅因子F。因此，根据足够长时间内得到的场和场的关系，可从发散场中提取2点的格林函数，这正是噪声成像的理论基础。之后，该方法在海洋声学（Roux et al，2004）和地震学（Shapiro et al，2004、2005）等领域得到了迅速的发展。

数学上，频率域与时间域的互相关计算是等价的。若2个序列x（t）、y（t）的傅里叶变换分别是X（f）、Y（f），则频率域的互相关公式为

其中，X*（f）为X（f）的复共轭；而C（f）的傅里叶反变换正是时间域的互相关c（t）。

本文将1天的数据分成48段×30min（50%重叠），在频率域对每个台站对的每段数据做互相关运算，并叠加形成每天的经验格林函数CCF。在测量波速时，会因为与REF相关度较差的CCF而产生强烈的变化；因此我们分别以10、30、50天作为滑动窗，叠加获得 CCFs，以提高相关度，获得更为可靠的波速变化。图5中，以BAJ-GHS台站对为例，展示了互相关提取的经验格林函数以50天为滑动窗叠加后的能量干涉图，并绘出了不同长度的滑动窗获得的CCFs与REF的互相关系数曲线。由图5可见，10天滑动窗的互相关系数曲线抖动较明显，30、50天滑动窗的互相关系数曲线较为平滑，能够突出显著的变化。

图5 BAJ-GHS台站对50天叠加的CCFs的能量干涉图（a）、紫坪铺水库水位变化（b）、CCFs与REF互相关系数的变化（c）

3.2 计算走时延迟Δt

为了研究库区水位对地下不同深度介质波速的影响，分析水位压力及渗透的作用过程，我们分别对 1～2s（0.5～1.0Hz）、2～4s（0.25～0.50Hz）、4～8s（0.125～0.250Hz）等 3个周期进行计算。根据瑞雷面波的特性（图6），上述3个周期的波分别对深度0～2、1～4、1～8km的介质变化敏感。

图6 瑞雷面波波速在不同周期的敏感性（据Liu等（2014））

3.2.1 获取参考格林函数REF

获取参考格林函数REF有2种方法，即绝对叠加法和相对叠加法（Lecocq et al，2014）。由于水库蓄水导致的地下介质波速变化并不剧烈，因此我们选择绝对叠加法，即叠加所有的CCF成为REF。图7为获得的28对台站的参考格林函数，此处进行了反序叠加，旨在展示的方便。

3.2.2 计算走时延迟

本文使用移动窗口互谱方法计算相对走时变化，该方法由 Ratdomopurbo等（1995）提出，最早被Poupinet等（1984）用于研究地震对之间的相对速度变化，其优点是在频率域进行计算，可以明确相关函数中相关信号的带宽。Brenguier等（2008a、2008b）利用该方法先后研究了法国富尔奈斯火山地下介质的相对波速变化及2004年Parkfield地区6.0级地震前后沿圣安德烈斯断层地下介质的相对波速变化。Clarke等（2011）详细阐述了移动窗口互谱方法的原理和步骤，并验证了该方法监测地壳波速变化的分辨率和准确率。

首先，我们将序列CCFs与REF分成许多重叠的窗，去平均、尖灭处理，然后进行傅里叶变换得到Fcur（f）和Fref（f）。互谱X（f）被定义为

式中，f为频率，星号表示复共轭性。在频率域，两序列的相似程度通过能量密度间的相关系数C（f）来评估。由式（5）推导得出

图7 28对台站的参考格林函数（反序叠加）

式中，φ（f）为解缠相位，它与频率f成线性比例关系

对于不同的频率f和相位φ，可通过最小二乘加权线性回归获得m，对应的误差为em，对应的权重为wj

式中，Cj是相关系数；Xj为互谱。这种权重公式的优势在于，能够在相关系数相对为常数而互谱能量变化的情况下给出较为不同的权重。根据式（7）中Δt与m的关系，将m与em分别除以2π，即可获得台站对间对应于不同时间错动的走时延迟Δt及误差eΔt。

3.3 计算平均相对波速变化

假设相对波速变化Δv/v在空间上是均匀的，那么它与相对走时变化Δt/t是相反数，即Δv/v=-Δt/t（Ratdomopurboet al，1995）。考虑到台站间距、信号的能量强弱、CCFs与 REF的相关程度等因素，本文选择了互相关正负错动的8～20s部分（图5），此部分为面波及尾波。我们将28组台站对间不同时间错动的走时延迟Δt及对应的误差eΔt进行平均，得到库区地下介质平均的走时延迟及平均误差对于3个周期，我们分别通过最小二乘加权线性回归获得为了减小误差，拟合时要求直线强制通过原点（Lecocqet al，2014），则

其中，pi为权重

b对应的方差为

最终，根据式（9），我们获得了库区地下介质平均的相对波速变化（图 8）。

4 结论与讨论

本文的研究时间段包含了紫坪铺水库3次大规模蓄水和2次泄水的过程，通过噪声互相关提取了库区28组台站对间的经验格林函数CCFs。图5中以BAJ-GHS台站对为例展示了台站对间的互相关结果，这对台站跨越了西南至东北的整个库区（图3），十分具有代表性。从图5中可以看出，水位加卸载时CCFs与REF的互相关系数明显降低，并在时间上有一定延迟，这定性地表征了水位影响了库区地下介质波速的变化。本文利用互相关结果，通过移动窗口互谱方法，分别对 1～2s（0.5～1.0Hz）、2～4s（0.25～0.50Hz）、4～8s（0.125～0.250Hz）等3个周期，计算了28组台站对间不同时间错动的走时延迟，3个周期分别对应着0～2、1～4、1～8km的深度。最终，使用面波及尾波部分，通过加权线性回归捕捉到了库区地下介质平均的相对波速变化（图8）。可以看出，库区地下介质相对波速变化不超过0.1%。水库水位的加卸载主要通过压力和渗透2种作用影响地下介质。Biot（1956）认为，孔隙流体会增加压缩波的速度并降低剪切波的速度。库区以碳酸盐岩为主，具备较好的渗透条件。图8的3个周期中，1～2s对应的2km内的介质相对波速变化最大，变化幅度超过了0.05%，这是由于在同样的水位压力下，浅层的渗透作用更强烈。

图8 库区地下介质平均相对波速变化

4.1 第1次蓄水及泄水

由图8可见，在2005年9月30日第1次大规模蓄水前，水位与相对波速变化没有明显关系。蓄水后数10天内，水位由760.36m迅速升高至835.91m，3个周期对应的3个深度的介质相对波速随之同时迅速降低，说明此时库水压力作用是导致相对波速变化的主要原因，而渗透作用只影响了2km以内的浅层介质。之后，水库开始泄水，对应的相对波速也逐渐升高。

4.2 第2次蓄水及泄水

第2次蓄水开始于2006年8月15日，2006年10月14日水位上升到最高值875.18m，前述3个深度的介质相对波速随之降低，存在明显的时间延迟。其中，1～2、2～4s对应的相对波速几乎同时降到最小值，而4～8s的相对波速在近1个月后才降到最小值。这说明在承受同样的库水压力下，渗透作用已成为影响相对波速变化的主要因素，且已影响至深度4km左右的断层。在第2次泄水过程中，3个深度对应的相对波速变化略有延迟，但几乎同时升高至最大值，这说明此时渗透作用已影响至深度8km左右的断层。

4.3 第3次蓄水

第3次蓄水于2007年12月12日上升到873.39m，3个周期的相对波速几乎同时降到最小值，与水位变化非常一致，互相关系数分别达-0.81，-0.66，-0.76（图9）。此时库水压力和渗透作用共同影响着介质相对波速的变化，水已渗透至深度为8km以下的断层。

图9 第3次蓄水时的水位与相对波速变化的最小二乘拟合曲线

综上，本文分别对紫坪铺水库3次蓄水及2次泄水过程进行了分析，讨论了不同阶段影响相对波速变化的主要因素，追踪了水在不同深度介质中的渗透过程。结果表明：在紫坪铺水库大规模蓄水前，相对波速变化与水位变化没有明显的相关性；而在之后的水位加卸载过程中，表现出了明显的负相关性，即水位升高，相对波速降低，水位降低，相对波速升高。分析认为，第1次蓄水时，压力是导致介质相对波速变化的主要因素；而后2次蓄水时，渗透作用成为主要因素。在第2次蓄水至最高点时，渗透作用影响至深度4km左右的断层，在第2次泄水至最低点时，渗透作用已影响至深度8km左右的断层。

卢显等（2010）通过研究水库水位与水库区域地震频度的关系发现，在2005年9月、2006年9月和2007年9月这3个水位快速上升的月份，地震频次也相应增加，且地震的震源深度绝大部分集中在5km处。根据本文的研究，这3个月份对应着相对波速由0值左右迅速降低的过程。因此对于地下浅层介质相对波速连续变化的研究，有可能作为水库地震预测的依据，具有十分重要的意义。

致谢：本文使用的大量连续波形及水库资料均由四川省地震局提供。研究过程中，本文得到了马延路博士、杨志高博士、周龙泉研究员、杜方研究员、戴仕贵高工、谢蓉华高工和韩进高工的支持与帮助，在此一并表示衷心的感谢。同时，非常感谢审稿专家提出的宝贵意见。