大岗山水库库区蓄水前后三维速度成像研究

庄园旭，丁 勇，杜 瑶，冯 静，康 宏，王余伟

（四川省地震局，四川 成都 610041）

大岗山水电站是大渡河干流的大型水电工程之一，库区位于青藏高原东南缘，地属川西南高原区。库区内山顶海拔一般为2 000 m以上，库区右岸最高峰贡嘎山海拔为7 556 m，左岸马鞍山海拔为4 021 m，坝区枯水期河水位952～957 m，水面宽40～80 m。电站正常蓄水位1 130 m，最大坝高210 m，总库容7.42亿m3，是典型的高山峡谷型高坝大库容水库。大岗山水库于2014年12月导流洞下闸开始第一阶段蓄水，2015年5月导流底孔下闸开始第二阶段蓄水。

库区位于鲜水河断裂带、龙门山断裂带和安宁河断裂带交汇部位，构造活动强烈，地震地质背景复杂，是高烈度地震危险区，穿过水库库区的主要活动断裂包括大渡河断裂、磨西断裂，其中磨西断裂历史上1786年曾发生73/4级地震，震中距大坝仅54 km。大岗山水库具有高坝大库，且位于地震危险性极高的磨西8级地震潜在震源区内的特点，对大岗山水库地震活动监测和相关研究工作十分必要。

1 方法与数据选取

根据射线理论，地震i到地震台k的体波到时T可以表示：

式中，τi为地震的发震时刻，u为慢度场，ds为射线的路径微元。

对于近震层析成像，震源坐标（x1，x2，x3）、发震时刻、射线路径、慢度矢量都是未知的。地震到时和地震定位的关系是非线性的，为了方便计算，采用Taylor级数展开的方法使其线性化，同时用三维网格速度模型，就可以把上式方程写成一个线性的等式：

对两事件i和j到同一台站k的方程作差，有：

假设这两个事件彼此靠近，可以认为事件到同一台站的射线路径几乎相同并且速度结构是已知的，则等式（3）可以简化为：

方程（5）既可以用震相目录的绝对走时，也可以用互相关数据的相对走时差。

双差层析成像方法[1]是在双差定位方法[2]和传统走时层析成像的基础上提出的一种走时层析成像方法。与传统的层析成像方法相比，该方法利用绝对到时以及相对到时实现了三维速度结构的反演和地震重定位，因为震源位置和速度结构存在耦合效应，并且真实的速度结构未知，利用绝对走时数据确定震源区域以外区域的速度结构，利用到时差数据确定震源区的精细速度结构，双差层析成像可以通过适当调整速度结构来提高地震定位的精度，地震定位的相对精确可以降低由于射线路径的不同造成的速度结构反演的误差，从而提高速度结构反演的精度。因此相比于传统的层析成像方法（仅采用绝对到时）双差层析成像方法可以得到更精确的速度结构和地震定位结果。自该方法提出以来，被广泛应用在区域成像中。邹振轩等[3]利用浙江数字地震台网大量地震记录，通过地震波双差走时层析成像，研究浙江省及周边地区上地壳精细速度结构和局部各向异性，并探讨地壳速度结构与活动断裂及地质构造的关系。邓文泽等[4]利用川西流动地震台阵、汶川地震震后应急台网记录到的P波到时资料，对2008年5月至2008年10月期间发生的汶川地震余震序列应用双差层析成像方法进行了地震震源和三维P波速度结构的联合反演，确认了汶川地震的发震构造特征。Dixit等[5]利用双差层析成像的方法，获得了印度Koyna-Warna水库区域高分辨率三维速度结构，并讨论了库水渗透、地震分布和速度异常之间的关系。罗佳宏[6]利用三峡加密观测台网记录的2009年3月到2010年12月高精度地震数据，使用双差层析成像方法，联合反演了三峡库区上地壳的震源位置参数和三维上地壳P波和S波速度结构，讨论了速度结构、微震分布和断裂的关系。

大岗山水库地震台网于2013年3月建成投测。大岗山水库地震台网由8个测震台站组成，布局孔径NS 20 km×EW 40 km，台站沿库区均匀展布，并全部包围了可能诱发地震的重点监视区段，区内监测震级下限为ML0.5级。大岗山水库台网在台站地震观测系统上全部采用反馈式短周期地震计和24位数据采集器，观测频带2s～40Hz。同时，大渡河流域还分布有瀑布沟水库台网、泸定水库台网、黄金坪水库台网等，如图1所示，黑色小方框为研究区（101.9°～102.5°E，29.2°～29.9°N）。此外，四川区域数字地震台网在距离大岗山大坝100 km范围内布设有6个宽频带地震仪，记录了大量观测资料均可作为本研究的有效补充。

图1 研究区域位置及地震台站分布

本文利用大岗山水库数字地震台网及其他台网2014年2月1日至2017年12月31日在研究区内记录到的ML≥0级地震事件12 764次，采用双差层析成像的方法进行反演。在速度结构反演之前，先要在研究区建立直角坐标系并设置反演网格点，研究区坐标系X轴与Y轴分别沿EW向和NS向，没有进行旋转，坐标中心点为大岗山水库大坝（102.2191°E、29.4497°N）。

由于研究区地震分布和台站分布不均匀，根据初始数据显示的研究区地震分布特点和射线覆盖路径，如图2所示，本文采用非均匀网格点设置，X轴网格点为-40 km、-35 km、-30 km、-27 km、-24 km、-22 km、-20 km、-18 km、-16 km、-13 km、-10 km、-5 km、0 km、5 km、10 km、20 km、30 km、40 km，共18个节点；Y轴网格点为-40 km、-30 km、-20 km、-10 km、-5 km、0 km、5 km、10 km、13 km、16 km、18 km、20 km、22 km、24 km、27 km、30 km、35 km、40 km、50 km、60 km，共20个节点，Z轴网格节点为0 km、2 km、5 km、10 km、18 km、28 km、40 km。网格点原点为大岗山水库大坝位置。一维速度模型采用赵珠等[7]给出的四川西部地区P波速度模型，见表1，波速比Vp/Vs取1.73。

表1 研究区域地壳P波速度模型

图2 重定位前射线路径

2 结果及分析

2.1 大岗山库区地震活动性分析

从2014年2月1日至2017年12月31日，研究区内台网共记录到12 764次ML≥0的地震事件，其中ML0～0.9地震9 317次，ML1.0～1.9地震3 097次，ML2.0～2.9地震315次，ML3.0～3.9地震32次，ML≥4.0地震3次分别为2016年3月18日4.2级、2016年3月24日4.0级、2017年1月4日4.0级。大岗山水库分两个阶段蓄水，在第二阶段蓄水位达到正常蓄水位后，水位保持相对恒定。在2015年4月、2016年3-4月、2017年3月出现了地震活动增强的现象，如图3、图4所示。

图3 大岗山库区蓄水位与地震日频次关系图

图4 大岗山库区N-T图

2.2 反演参数选择

本文分为四个阶段对库区三维速度结构和地震分布进行精确定位。第一阶段为蓄水前2014年2月1日至2014年12月31日，第二阶段为蓄水过程2015年1月1日至2015年12月31日，第三阶段为蓄水后2016年1月1日至2016年12月31日，第四阶段为2017年1月1日至2017年12月31日。

四个阶段反演均采用相同的初始速度模型和相同的网格节点，反演中共设置了6组迭代，每组迭代2次，迭代的截断值设为9倍标准差，并同时反演速度结构和地震位置。另外，由于数据量较大，反演解算时使用带阻尼LSQR算法，该算法中平滑系数smooth和阻尼系数damp约束着地震位置和慢度的变化量，对收敛速度和结果平滑程度影响较大。为此，本文采用L-curve方法进行测试以找到合适的系数值，通过归一化模型与走时残差关系曲线（图5）分析显示，当固定平滑系数为100模型较平滑，阻尼系数分别为80、80、140、60时走时残差也相对较小，因此在反演成像中采用这两个值，在合成分辨率测试中smooth取10，damp取值和反演过程取值相同。

图5 阻尼系数选择曲线

2.3 棋盘测试和反演结果

为确保成像结果的可靠性，在实际数据反演之前，采用棋盘格测试[8]的方法对解的质量进行分析，本文棋盘测试中网格设置与反演时网格设置相同，速度扰动取正常值的±5%，图中黑色表示正速度异常，白色表示负速度异常，棋盘测试结果如图6、图10、图14、图18所示，由于地震分布集中在鲜水河断裂西侧并且台站沿大渡河两岸分布，造成无地震和台站区域的检测板测试恢复较差，本文中检测板测试结果恢复较好的范围为X轴-30～10 km，Y轴-10～30 km区域，在深度上，2～10 km深度范围具有较好的模型分辨率，本文只使用恢复较好区域成像结果进行解释。

图6 第一阶段不同深度V p（左）V s（右）棋盘测试

图10 第二阶段不同深度V p（左）V s（右）棋盘测试

图14 第三阶段不同深度V p（左）V s（右）棋盘测试

图18 第四阶段不同深度V p（左）V s（右）棋盘测试

第一阶段为蓄水前，此阶段参与反演地震数2 766条，重定位后获得2 406次地震的精确定位结果，均方根残差从0.415 1 s降低到0.1 s，X、Y、Z三个分向震源位置估算误差降低到336 m、229.5 m、459 m。不同深度震源位置和速度结构切片如图7所示，从图中可以看出，大岗山库区的速度结构存在明显的横向不均匀性。地震较为集中的X=-30 km到X=-10 km，Y=5 km到Y=25 km的磨西断裂西侧区域（以下称为A区），在5 km以上的浅层区域P波速度和S波速度均表现为高速异常，5～10 km深度范围高Vp范围减小，低Vp范围增加。X=-10 km到X=0 km，Y=0 km到Y=10 km的水库蓄水区（以下称为B区），5 km以上P波、S波速度都较低，5～10 km范围表现出P波速度与A区相当，S波依然表现为低速。

图7 第一阶段不同深度V p（左）V s（右）速度结构

第一阶段地震主要分布在A区，呈丛集状分布，定位结果与杜瑶等[9]使用双差定位法得到的结果基本一致。震源深度在10 km以内，呈丛集状分布，且位于高Vp、高Vs区域，如图8、图9所示。从地震分布上看，与磨西断裂关系不大，推测存在小型隐伏构造。

图8 大岗山库区第一阶东西向P波速度和S波速度剖面

图9 大岗山库区第一阶南北向P波速度和S波速度剖面

第二阶段为开始蓄水至正常蓄水位，此阶段参与反演地震2 700个重定位后获得2 497次地震的精确定位结果，均方根残差从0.384 4 s降低到0.08 s，X、Y、Z三个分向震源位置估算误差降低到227.7 m、206.4 m、333.1 m。不同深度震源位置和速度结构切片如图11所示，A区在5 km以上的浅层区域P波速度和S波速度均表现为高速异常，5～10 km深度范围高Vp范围减小，低Vp范围增加。B区5 km以上P波、S波速度都较低，5～10 km范围P波、S波均表现为低速。A区地震分布与第一阶段相似，呈带状分布，震源深度主要在7 km以内，且处于P波、S波高速区，根据图12、图13震源深度特征依然能看到与第一阶段类似的隐伏断裂特征。B区出现明显的地震丛集，震源深度在2～7 km范围内，如图22、图23所示。

图11 第二阶段不同深度V p（左）V s（右）速度结构

图12 大岗山库区第二阶段东西向P波速度和S波速度剖面

图13 大岗山库区第二阶段南北向P波速度和S波速度剖面

第三阶段参与反演地震4 980个重定位后获得4 678次地震的精确定位结果，均方根残差从0.374 s降低到0.071 1 s，X、Y、Z三个分向震源位置估算误差降低到144.9 m、136.2 m、201.2 m。不同深度震源位置和速度结构切片如图15所示，速度分布与第二阶段类似，A区在5 km以上的浅层区域P波速度和S波速度均表现为高速异常，5～10 km深度范围高Vp范围减小，低Vp范围增加。B区5 km以上P波、S波速度都较低，5～10 km范围P波、S波均表现为低速。A区地震分布与前两个阶段相似，呈带状分布，如图16、图17所示。此阶段B区地震数目增加，丛集性更明显，震源深度主要分布在2～8 km范围内，如图22、图23所示。

图15 第三阶段不同深度V p（左）V s（右）速度结构

图16 大岗山库区第三阶段东西向P波速度和S波速度剖面

图17 大岗山库区第三阶段南北向P波速度和S波速度剖面

第四阶段参与反演地震2 311个重定位后获得2 128次地震的精确定位结果，均方根残差从0.365 s降低到0.065 9 s，X、Y、Z三个分向震源位置估算误差降低到208.3 m、207.6 m、282.0 m。不同深度震源位置和速度结构切片如图19所示，A区、B区速度分布和地震分布与第二、三阶段相似，在此不再赘述。

图19 第四阶段不同深度V p（左）V s（右）速度结构

针对蓄水后B区在蓄水后地震活动增强的现象，我们做了Y=5 km、X=-5 km两个剖面，如图22、图23所示，第一阶段地震分布散乱，从第二阶段开始，B区出现了地震丛集，震源深度主要集中在2～8 km，B区在0～5 km范围P波速度和S波速度均没有明显变化，而在5～10 km出现的S波高速异常，是反演过程中地震丛集过于集中该区域射线密度大而其他区域射线分布少造成的。

图22 四个阶段库区东西向P波速度和S波速度剖面

图23 四个阶段库区南北向P波速度和S波速度剖面

图20 大岗山库区第四阶段东西向P波速度和S波速度剖面

图21 大岗山库区第四阶段南北向P波速度和S波速度剖面

3 讨论与结论

综上所述，本文利用水库地震台网2014年2月1日至2017年12月31日在研究区（101.9°～102.5°E，29.2°～29.9°N）内记录到的ML≥0级地震事件12 764次，采用双差层析成像的方法分别对蓄水前、蓄水过程、蓄水后等四个阶段进行反演，获得了水库蓄水区及其临区三维速度结构和精确的地震定位结果。棋盘测试结果显示，各个阶段磨西断裂西侧的A区（X=-30 km到X=-10 km，Y=5 km到Y=25 km）和大岗山水库蓄水去的B区（X=-10 km到X=0 km，Y=0 km到Y=10 km）均有较高分辨率。

从蓄水各阶段成像结果来看，研究区速度结构存在明显的横向不均匀性，各个阶段A区P波、S波速度均高于B区，结合地质资料，B区分布有3条断裂，且大渡河断裂穿过大渡河，河水可能沿断裂渗透，造成该区域P波、S波速度较低。B区各个阶段在0～5 km的浅层P波、S波速度均为未发生较大变化，这与常见的水库蓄水后速度结构发生较大变化不同[5-6，10-12]，分析原因，因大渡河断裂穿过大渡河，在蓄水前，河水可能就已经渗透到地下，孔隙水已达到饱和，而在蓄水后，浅层的水含量未发生较大变化，因此蓄水前后速度结构未发生较大变化。

从蓄水各阶段地震精确定位结果来看，A区地震分布集中且呈带状分布，震中分布方向与磨西断裂走向一致，从各阶段震源深度分布可以清晰分辨出A区存在一条走向与磨西断裂一致倾向WE的隐伏断裂，该区域地震主要受本底构造活动影响。B区蓄水前地震稀少，蓄水后即出现了地震活动增强的现象，震中分布呈丛集状，震源深度主要在2～8 km范围且各阶段震源深度未发生明显变化，P波、S波速度也未发生明显变化，因此B区出现的地震丛集可能是由于库水荷载引起局部应力变化造成的[13]。