云南小湾水库蓄水后P波速度结构的双差地震层析成像研究

曹颖 黄江培 付虹

云南省地震局，昆明市北辰大道148号 650224

0 引言

目前，对水库诱发地震的形成有多种认识和解释，但都一致认为，水库所在区域地质和地球物理环境、库区从地表到深部的岩石组成和物质结构状态、区域断裂构造分布、产状、力学性质、发育程度、现今活动性等，都是发生水库诱发地震的重要因素（任金卫，1993）。近年来，随着数字地震观测技术的发展，深部地球物理探测研究是获得区域深部构造环境和孕震环境的重要手段。层析成像方法不仅能用于一般区域，对于水库地区也同样适用。利用层析成像方法探测地下速度结构，并结合地震震源分布揭示地下介质的特点，为研究水库诱发地震的形成提供了新的途径。王亮等（2015）利用速度结构和小震精定位结果得出紫坪铺水库水的渗透作用对汶川地震的发生没有直接影响的结论；钟羽云等（2010）利用速度结构和震源参数结果研究认为，由于水库蓄水后库水下渗，地震大多发生在低速异常区内，水库诱发地震之初的几年中震源深度有一个逐渐变大的过程；李强等（2009）对三峡水库坝址及邻区中上地壳P波速度结构研究发现，库水渗透作用对地壳浅层速度结构产生影响。

tomoDD双差地震层析成像方法由Zhang等（2003）首先提出，该方法避免了复杂速度结构引起的射线路径异常对定位结果的影响，同时，利用观测走时、走时残差进行地震定位和速度结构联合反演，不仅能够获得较常规地震层析成像更为精确的震源密集区附近的速度结构，还能够获得地震发生的绝对位置，其地震定位的精度与双差地震定位的精度相当。在国外双差层析成像方法被广泛应用于地震区精细速度结构研究，并已获得了从区域尺度到几百米尺度的高分辨率成像和地震定位结果（Shelly et al，2006；Zhang et al，2004；Okada et al，2005；Pei et al，2010），在我国也有很多应用（于湘伟等，2010；李海鸥等，2011；邓文泽等，2014；王小娜等，2015；吕子强等，2016）。

云南小湾水库是滇西澜沧江中下游河段梯级电站的龙头水库，水库设计最大坝高292m，设计正常蓄水位高程1242m，水库总库容约149×108m3，水库区跨越大理、临沧、保山3个地区，由西支干流澜沧江和东支支流黑惠江组成，库长分别为178、124km，水库靠近库坝区的最高海拔约1160m，靠近库尾及库区正北的黑惠江流域最高海拔约1210m，总体呈东南低、西北高的地形特点，水库于2008年12月16日开始蓄水。小湾水库位于构造活动较强烈的地区，在唐古拉-昌都-兰坪-思茅褶皱系和冈底斯-念青唐古拉褶皱系的接合部位（毛玉平等，2004），库区出露的地层主要是变质岩和沉积岩，局部亦有岩浆岩发育（任金卫，1993）。库区及其周边主要断裂有怒江断裂、保山-施甸断裂、澜沧江断裂、兰坪-云龙断裂、维西-乔后断裂、红河断裂、程海-宾川断裂、无量山断裂、南汀河断裂、昌宁断裂、柯街断裂等（图1（b））。历史上该区域地震多发，距坝址仅75km处有7级强震的活动记录。对于小湾水库的研究已有很多，邬成栋等（2010）利用小湾电站水库诱发地震监测台网记录的数字波形资料，使用遗传算法反演了小湾水库近场328个中小地震的地震矩、应力降和拐角频率等震源参数；李永莉等（2012）计算了小湾水库蓄水前后的水库地震波速比变化；姜金钟等（2016）利用地震精定位方法分析了小湾水库的地震活动性；柯乃琛等（2016）使用 Velest程序计算得到小湾库区蓄水前的最小一维速度结构，然后以该最小一维速度结构作为初始模型，使用Simul2000程序计算了小湾水库蓄水前后库区上地壳介质三维速度结构和震源位置。在上述研究的基础上，本文使用双差地震层析成像方法（tomoDD）联合绝对到时、相对到时计算小湾水库库区及周边区域水库蓄水后不同时间段内的地震重定位结果和三维P波速度结构，分析由蓄水所导致的该区域地下介质波速的变化，以期深化对小湾水库地震发震构造、发震机理等的认识。

1 资料与方法

1.1 地震资料

图1 研究区域活动断裂、地震分布、网格划分和所用台站分布

为监测小湾水库及其周边的地震活动，2005年6月1日建成并运行小湾水库地震台网12个一期台站，台站包围并均匀分布于库坝区，2009年9月1日开始运行二期台站4个，主要分布于库首区，2010年4月小湾一、二期台站合并运行。2012年8月1日漫湾电站水库地震台网建成，并与小湾水库台网联合运行。但从2014年开始水库台网运行率下降，基本无法记录到库坝区的地震，所以在本文中加入了距研究区内地震约250km范围内的云南区域地震台网的49个台站，包括30个固定台站、3个腾冲火山台网台站、3个下关小孔径台网台站以及13个流动台站（图1（a））。本文选取的研究范围为小湾水库库区及其周边地区（24.5°～25.5°N，99.0°～100.5°E），数据为至少有 6个台站记录到的地震事件，时间段为小湾水库蓄水后的 2008年12月16日～2016年12月31日，经过筛选共挑选出 11249个ML≥－0.5地震（图1（b）），再去除走时曲线中离散较大的震相，最后得到77306个P波绝对到时资料（图2）。

图2 所选取的走时-震中距曲线

由小湾水库水位资料（图3）可知，小湾水库开始蓄水后水位快速上升，至2009年7月，水位由蓄水前的996m上升至1100m，2010年10月水位上升至1210m，由图3可见，2次水位快速上升期间库区内部分区域的地震活动频次明显增多。随后，水位以年为周期在1200～1250m间变化（姜金钟等，2016）。根据水位变化，选取2008年12月16日～2011年6月30日及2011年7月1日～2016年12月31日2个时间段来分析研究区地下速度结构，其中，第1个时间段是蓄水后水位快速上升至平稳变化的时间段，在这个时间段内地震活动明显增多，第2个时间段是每年水位周期变化的时间段。由柯乃琛（2016）的研究可知，从小湾水库开始蓄水至2011年6月地震频次大量增加，此后一直保持在较高的水平。第1个时间段内符合选取条件的地震事件有4474个，并产生31921条绝对P波到时，最大震级地震为2010年6月1日施甸ML4.8地震；第2个时间段内有6775个地震，产生45385条绝对P波到时，最大震级地震是2015年10月30日昌宁MS5.1地震，其震中距水库大坝仅77km，距澜沧江仅11km。在此基础上进行地震对匹配，选择每次地震与2个地震之间的最大距离为20km，每个地震最多可与30个地震组成地震对，最终构建了相关到时数据。其中，2008年12月16日～2011年6月30日有270189条P波相对到时，2011年7月1日～2016年12月31日有381176条P波相对到时。

图3 小湾水库水位变化与地震月频次

1.2 双差地震层析成像方法

反演中选取的坐标原点为（25°N，99.75°E）。在划分网格时，较小的网格分辨率较高，但若网格太小又会影响地震层析成像图像的质量，所以，根据所选地震事件和台站的分布情况，在反演之前进行了大量的分辨率测试，以寻求最佳的网格分布，我们分别测试了5km×5km、10km×10km、15km×15km、20km×20km的网格间隔。虽然2个时间段使用的台站不一样，但差别并不大，且2个时间段的地震分布密集区基本相同，均为黑惠江和昌宁-施甸一带，2个时间段的横向、纵向分辨率均可达到10km，研究区东侧和北侧外围地震分布较稀疏，故使用 20km×20km的间隔（图1（b））。2个时间段垂直向网格相同，均为 0、3、5、10、15、20、25、30km。选择一维模型作为初始模型，初始一维速度模型为使用Kissling方法得到的小湾水库最小一维速度模型（表1）（柯乃琛等，2016），2个时间段均使用该初始模型。空间任意点的速度利用线性插值求得（Thurber，1983）。

由于阻尼参数和平滑权重的大小对反演结果的稳定性有较大影响，因此，对不同平滑权重和阻尼参数进行了权衡分析，构建不同平滑权重、阻尼参数的解的方差与数据方差之间的均衡曲线，图4给出了2个时间段的均衡曲线，2个时间段的平滑权重和阻尼参数的搜索范围一致，分别为0～600、10～2000。图4中所示的解的方差和数据方差仅有相对意义，解的方差包含了事件定位和速度模型参数的影响。由图4可见，对于不同的平滑权重，最佳阻尼参数约为750；由阻尼参数为750时的一系列平滑权重的数据方差与速度模型方差的均衡曲线得到的最优平滑权重约为150。

表1 一维P波速度和波速比模型表

2 层析成像结果

2个时间段的数据经过11次迭代后，第1个时间段2008年12月16日～2011年6月30日的数据在最后一次反演后拾取的到时残差（RMSCT）从0.36s降到0.15s；第2个时间段2011年7月1日～2016年12月30日的从0.28s降到0.06s。图5为2个时间段定位模型的残差变化图，以0.1s为间隔统计所有地震的均方根残差。由图5可见，初始均方根残差分布较为分散，最后均方根残差则紧密地分布于0s附近。由于初始定位结果采用简单的一维速度模型并且用于定位的数据类型较少，仅使用了绝对到时数据，因此地震定位精度较低；而使用地震波绝对、相对到时数据联合反演震源参数以及速度结构会使得反演后的地震定位及速度模型精度都有了显著的改进。

2.1 重定位结果

图4 平滑权重和阻尼参数均衡曲线

图5 蓄水后2个时间段的残差变化

从研究区域（图1（b））2个时间段地震重定位结果来看（图6），研究区内主要有 A、B、C、D、E等地震聚集区域，A区为黑惠江，B区为小湾水库回水澜沧江段，C区为澜沧江昌宁段，D区为施甸一带，E区为距澜沧江仅11km的2015年10月30日昌宁MS5.1地震余震序列分布区。第1个时间段的地震主要分布在 A、B、C、D区，由姜金钟等（2016）及柯乃琛（2016）的研究结果可知，A、B、C三个区域蓄水前地震活动频次不高，蓄水后频次开始增加，并随水位涨落而起伏。D区地震大多为2010年6月1日施甸ML4.8地震余震序列，且D区的地震一直以来都较多，与水库蓄水前后水位变化间的关系不明显。由图6（c）可见，A、B、C三个区域大部分地震的震源深度不超过10km；D区地震震源深度为0～30km，与姜金钟等（2016）的结果相符。

图6 重定位后地震分布

第2个时间段的地震分布与第1个时间段大致相同，分别是A、B、C、D、E区，D区地震一直较多，其中包含了2012年9月11日施甸MS4.7地震余震序列，E区主要为距澜沧江仅11km的2015年10月30日昌宁MS5.1地震序列，该地震为水库蓄水后发生的震级最大的地震。由图6（b）可见，A、B、C区地震的震源深度大多小于10km，与第1个时间段差不多，但A区发生的2016年1月4日ML4.0地震的震源深度约为15km。C区昌宁MS5.1地震的余震序列震源深度大多小于10km。D区施甸MS4.7地震余震序列的震源深度为0～30km。2个时间段内小湾水库库坝区地震活动均未有明显增强，表明小湾水库大坝位置地质构造较为稳定。

由姜金钟等（2016）的研究结果可知，A、B、C三个区域地震的震源深度平均值均小于10km，D区地震的震源深度平均值大于10km，小于15km，姜金钟等（2016）将此结果与天然构造地震的深度结果（张国民等，2002；杨智娴等，2004；Yang et al，2005）进行对比，得出 A、B、C三个区域存在明显的水库触发地震，而D区域的地震则属于构造地震。

2.2 P波三维速度结构

利用棋盘测试方法评价成像结果的分辨率。在初始一维速度模型上分别增加±6%的正负相间速度扰动，从而获得初始三维速度模型，根据实际地震数据正演计算合成走时数据集，使用合成走时数据集进行双差地震层析成像，通过对反演结果中每一个节点的速度与理论值进行对比，了解成像结果的分辨率。对2个时间段都进行了棋盘测试，由于20、25、30km深度处的反演结果基本为空白或失真，所以分辨率测试只选取了0、3、5、10、15km剖面。图7为2个时间段的棋盘测试结果。由图7可见，由于不同深度处射线的分布不同，因此，0km深度处恢复不好；3km深度处的中心区域恢复较好；5、10、15km深度处恢复都较好。

选择0、3、5、10、15km深度处水平层面的速度分布进行讨论，对于2.1节所讨论的地震集中分布的区域将着重分析。由0km深度水平层面P波速度结构（图8）可知，2个时间段研究区的速度异常与地表地形起伏之间有密切关系，高速区基本沿着澜沧江分布，而澜沧江两岸均为高山。

由图8可见，2个时间段3km深度处地震均集中分布在 A、B、D区，D区第1个时间段地震大多为2010年6月1日施甸ML4.8地震的余震序列，第2个时间段地震分布较为离散，其中包括2012年9月11日施甸MS4.7地震的余震序列。所有这些地震集中分布的区域均为高速区，第1个时间段该高速区的 P波速度约为4.5km／s，小湾水库库区的 P波速度为3.8～4.0km／s；第2个时间段地震分布集中区域的P波速度为4.2～4.5km／s，水库库区P波速度为 4.2～4.5km／s。

2个时间段5km深度处地震仍然主要分布于A、B、D区，D区2个时间段的地震仍然主要是2次施甸地震的余震序列，第2个时间段E区地震主要为距澜沧江仅11km的2015年10月30日昌宁MS5.1地震的余震序列。A区、B区和库坝区第1个时间段P波速度为5.8～6.0km／s，第2个时间段 A区和库坝区为5.6～5.8km／s，B区为5.9～6.0km／s；C区第1个时间段速度为6.0～6.1km／s，第 2个时间段为 5.8～6.0km／s；D区第 1个时间段速度为 5.9～6.0km／s，第 2个时间段为 5.8～5.9km／s。

图7 不同深度处的P波速度棋盘测试结果

2个时间段10km深度处A区基本为低速区，P波速度均为5.8～6.0km／s；库坝区及B区为高速区，B区第1个时间段P波速度为6.2～6.4km／s，第2个时间段为6.1～6.4km／s，库坝区第1个时间段P波速度约为6.6km／s，第2个时间段为6.2～6.4km／s；C区第1个时间段 P波速度约为5.8km／s，第2时间段约为5.7km／s；D区第1个时间段 P波速度为6.1～6.3km／s，第2个时间段为5.8～6.2km／s；E区2个时间段P波速度均约为6.0km／s。

图8 不同深度的P波速度、P波速度差和地震震中分布

至15km深度处，D区2个时间段地震分布仍相对较多，而0～10km深度地震分布集中的A、B、C区在15km深度处却很少有地震发生。该深度处第1个时间段A区P波速度为5.9～6.0km／s，第 2个时间段为 6.0～6.2km／s；B区第 1个时间段 P波速度约为5.8km／s，第 2个时间段约为6.4km／s；库坝区2个时间段P波速度为6.1～6.3km／s；C区第1个时间段P波速度约为5.8km／s，第2个时间段约为6.2km／s；D区第1个时间段 P波速度为6.0～6.2km／s，第2个时间段为 6.2～6.4km／s；E区第 1个时间段 P波速度约为5.6km／s，第2个时间段为5.6～5.8km／s。

为了对比2个时间段研究区内P波速度的变化，将第2个时间段的速度减去第1个时间段的速度得到图8（c）所示结果。由图8（c）可见，A区在0、3、5、10km深度处 P波速度均为下降，本文认为这可能是由于黑惠江是一条年轻的支流，蓄水后库水沿着该区域岩石破碎形成的通道向四周进行渗透，使得周边介质空隙含水率增加，强烈地改变了介质的性质，造成了此处P波速度呈现低速异常变化。对于B区，0～10km深度内，在地震分布的区域，P波速度随深度增加而降低，至5km深度处，速度降低幅度很小，至10km深度处，基本为P波速度上升，可能是由于该区域地下岩体主要为二长花岗岩，库水不易渗透，最多只能渗透至5km深度。C区在0、3、5、10km深度处P波速度均为下降，C区位于澜沧江断裂与柯街断裂的交界处，由柯乃琛（2016）对小湾水库蓄水水位与地震活动性间的相关性研究可知，蓄水回水位沿着澜沧江超过了C区后，C区的地震大量增加，说明该区域P波速度降低的原因是库水沿澜沧江断裂与柯街断裂的岩石破碎带向下渗透从而改变了地下介质的性质。D区在0、3、5、10km深度处P波速度均为下降，但是D区远离澜沧江，蓄水前后地震活动都较活跃，故本文认为该区域构造活动较为活跃，与水库蓄水无关。E区在0、3、10km深度处P波速度降低，5km深度处地震主要分布于速度降ΔvP上升的区域，E区位于柯街断裂上，E区发生的地震与蓄水回水无关，初步判定是构造地震。水库库坝区的P波速度变化没有规律，表明了该区域地质构造的稳定性。

图9为 A-A′、B-B′、C-C′、D-D′、E-E′剖面的 P波速度结构分布。由图9可见，2个时间段内地表低速层基本不变，A-A′剖面上的小湾水库回水澜沧江段的地震集中区域在第1个时间段震源深度基本在5km以内，第2个时间段的地震震源深度大于5km，小于10km，P波速度未有太大的变化。B-B′和C-C′剖面均横穿过黑惠江上的地震集中区域，2个剖面上地震聚集区的震源深度在第1个时间段内均小于5km，而在第2个时间段内大于 5km，小于10km，在地震聚集区的P波速度稍微变小。D-D′剖面上红色箭头所标注地即为澜沧江保山段上地震聚集区域，该地震聚集区的震源深度在2个时间段内均很浅，大多并未超过10km，第2个时间段震源聚集区域低速区变厚。E-E′剖面上红色箭头所标注地为2015年10月30日昌宁MS5.1地震余震区，在第1个时间段内该区域并没有地震，余震深度小于10km，P波速度基本无变化，第1个时间段内E-E′剖面西端的地震为2012年9月11日施甸MS4.7地震的一部分余震，该区域波速相对较低。

本文认为由不同地震聚集区在2个时间段的震源深度及P波速度变化可得到一些认识，即黑惠江上的地震集中区（A区）和小湾水库回水澜沧江段的地震集中区（B区），随着蓄水时间的延长，库水沿着破碎带向下渗透，地震震源深度随之增加。澜沧江保山段上地震聚集区（C区）的震源深度并未随蓄水时间延长而增加，2个时间段的震源深度分布相差不大，表明该区域的地下通道和岩性均更利于库水渗透，蓄水回水至该段后，库水渗透至较深处，而在10km深度以下地质构造较稳定，故库水难以向下渗透。

3 讨论

图9 P波速度沿垂直剖面深度分布

已有研究认为，只有当水体和库区的地质构造条件、水文地质条件、岩性条件和应力条件等因素适当地结合后才有可能触发较大的地震（丁原章等，1989；李安然等，1992；Simpson et al，1988）。水库水体的加载和渗透主要对库底岩石及断层产生3种作用（陈翰林等，2009），即弹性效应、孔隙压变化及断层弱化等，其中，水库蓄水载荷会导致断层面上正应力和剪应力增大，正应力增大会使得深部断层增强，而剪应力增加，断层变弱还是变强则取决于断层走向与区域应力场方向之间的关系。水体载荷的影响深度非常浅，载荷和渗透引起的孔隙压变化则首先增加浅部的孔隙压，继而使得孔隙压向深部扩散，从而改变介质的强度和断层的摩擦阻力。所以，结合地震的震源深度和速度变化可初步判断不同的地震聚集区在小湾水库蓄水后地震增加的原因。在本文中，基于已有的认识可知：

（1）A、B区的地震在第1个时间段内震源深度小于5km，在第2个时间段内大于5km，但并未超过10km，第2个时间段的P波速度相比第1个时间段下降，说明这2个区域的地震与水库蓄水有关，是由于水体渗透导致孔隙压变化，并随着时间的变化孔隙压变化朝着更深的部位扩散，从而导致介质变化，P波速度降低。

（2）C区的地震震源深度小于10km，大多在5km以内，并未随时间而向深部发展，第2个时间段P波速度相比第1个时间段下降，在蓄水回水位到达该区域后地震开始增多，说明C区地震的发生是由于库水的渗透作用，该区域地下10km深度以上渗透通道有利于库水的渗透，而10km深度以下地质构造较为稳定。

（3）小湾水库蓄水后，在水库大坝附近并没有明显的地震增多现象，2个时间段内的P波速度也无规律性的变化，表明了小湾水库大坝位置地质构造的稳定性。

（4）D区历来地震多发，由震源深度以及速度结构变化可知，震源深度为0～30km，P波速度变化无规律，说明蓄水后施甸一带发生的地震，如2010年6月1日施甸 ML4.8、2012年9月11日施甸MS4.7等地震均属于构造地震，与水库蓄水关系不大。

（5）发生在E区的2015年10月30日昌宁MS5.1地震及其余震震中位于柯街断裂与澜沧江断裂的交汇处，震源深度为0～10km，大多在5km以内，由主震震源机制解结果①云南省地震局监测中心，2014，大震应急数据产品产出判断，NE向柯街断裂为发震断裂，具正断性质。在第1个时间段内该区域基本无地震，2个时间段内该区域的P波速度上升。初步判定该地震是与小湾水库蓄水有关的构造地震，但需进一步对该震群单独作计算分析才能得到更准确的结果。

致谢：中国科技大学的张海江教授为本研究提供了tomoDD程序，两位审稿专家为本文提出了宝贵的修改意见，在此一并表示感谢。