2021, Vol. 64
2. 中国科学技术大学地球和空间科学学院, 合肥 230026
2. School of Earth and Space Sciences, University of Science and Technology of China, Hefei 230026, China
据中国地震台网测定，2014年8月3日16 ∶ 30在云南省昭通市鲁甸县发生了MS6.5地震(震中103.3°E，27.1°N，震源深度12 km)(CENC，http://www.cenc.ac.cn).此次地震震级虽不大，但地震极震区最高烈度可达Ⅸ度，且由于震源浅、震区人口密度高，是云南省近年来典型的“小震”致大灾事件(李西等，2018)，引起国内外的高度关注.鲁甸地震的研究结果，让我们对鲁甸地震的震源特性、运动学特征及发震构造等有了深刻的认识.重定位结果显示，鲁甸地震余震序列的震中分布主要为近EW向和SES向的共轭型分布，震源深度主要分布于13~15 km，主震位于破裂底部(张广伟等，2014；房立华等，2014；王光明等，2018)；结合震源机制解和震源破裂过程研究(张广伟等，2014；徐锡伟等，2014；张勇等, 2014, 2015；许力生等，2014)，可知此次地震为一次高倾角左旋走滑型地震，是在NW向主压应力与NE向主张应力相互作用下发生的以NNW向破裂为主的共轭破裂事件.此次地震周边地质构造复杂，大多数研究(徐锡伟等，2014；李西等, 2014, 2018；房立华等，2014；张勇等, 2014, 2015；程佳等, 2016, 付芮等，2015)在结合了野外考察、余震序列分布、震源机制解、地震烈度分布、库伦应力变化和震源破裂过程等研究结果后认为发震断裂为NNW向的包谷垴—小河断裂.以往也有对鲁甸地震发生前后地下速度结构的研究，Riaz等(2017)基于双差层析成像方法给出了鲁甸地震周边区域的地壳层析成像结果，发现鲁甸地震震源区低P波速度异常被高P波速度区域所包围，大多数地震发生在低P波速度和低泊松比区域，受主震破裂的限制.赵小艳和孙楠(2014)利用震源位置和速度结构的联合反演得出了鲁甸地震余震主要分布在P波速度较高的区域.戴宗辉(2016)利用地震背景噪声互相关方法得到了鲁甸地震前后的相对速度变化，发现较大幅度的同震波速变化主要集中于则木河断裂—小江断裂及莲峰断裂区域.庞卫东等(2017)利用背景噪声自相关方法得到了鲁甸地震前后震源区周边5个台站的相对速度变化，发现鲁甸地震具有明显的同震速度变化，推测可能是青藏高原东南缘与华南地块相互作用而造成了研究区地壳构造应力变化.虽然前人已经做了很多与地壳速度结构变化有关的研究工作，但是均缺乏高精度空间分辨率的结果，所以提供鲁甸地震发生前后高精度的地壳波速变化的空间分布变化，有助于我们更加深入的了解地震孕育、发生和恢复的过程.
技术的进步，尤其是地震台网空间密度的增加，使监测波速变化的灵敏度增加了一个数量级以上.监测波速变化有使用固定的人工震源和被动震源的方法.但是使用人工震源时，对所有危险区域进行监测的费用太高且在某些地区禁止这样做，因此大多研究采用被动震源的方法及数据.而使用被动震源数据的一种方法是使用剪切波分裂方法，这是通过各向异性介质传播的结果(Crampin et al., 1990；高原和吴晶，2008；常利军等，2010).一种方法是比较具有相似位置和机制的重复地震的尾波中的散射波(Poupinet et al., 1984；Schaff and Beroza, 2004；Taira et al., 2008；Zhao and Peng, 2009).还有近10年发展迅速的基于背景噪声的方法，该方法依赖于台站的分布，已经成功应用于研究不同区域的波速变化(Brenguier et al., 2008；Obermann et al., 2014；刘志坤和黄金莉，2010). 上述方法均不能得到速度变化的高精度空间分布，而能得到空间分布的一种方法是将不同时间段的速度模型相减(例如：Chiarabba et al., 2009；Patanè et al., 2006；Westman et al., 2012)，然而Julian和Foulger(2010)指出这种把不同时间段层析成像反演的速度模型相减的结果代表真实的速度变化的假设是值得怀疑的，因为由于地震位置的自然变化引起的地震射线分布变化，即使速度结构没有变化，层析成像结果也会不同，即使地震位置没有变化和地震台站分布是固定的，但由于随机观测误差，在反演出的模型中也会出现误差.更好的方法是同时反演多个数据集，所以Julian和Foulger(2010)提出了一种依赖时间的成像方法，该方法通过最小化模型的差异和不同时间段的到时残差来同时反演多个数据集，该方法对正演数据进行了测试，与传统的直接相减的方法相比，Julian和Foulger(2010)的方法能够得到更可靠的速度随时间的变化，但在一定程度上可能还会受到不同时间段数据分布不均匀的影响.Qian等(2018)提出了一种新的时移层析成像方法，该方法基于双差层析成像方法(Zhang and Thurber, 2003)，通过使用不同时间段的事件对构建的到时差，可以克服部分的观测误差的影响和不同时间段的地震射线分布差异的影响.在反演过程中通过使用复杂的迭代策略从而解决了射线分布差异所产生的假性变化，为了消除不需要的模型差异，一个时间段的速度模型被用作另一个时间段的初始模型，因此这种新的时移层析成像可以确定可靠的波速随时间的变化.
在本研究中，我们利用基于双差层析成像的时移层析成像方法，使用2008年1月至2017年12月的地震数据进行时移层析成像反演来得到2014年云南鲁甸MS6.5地震发生前后震源区及其周边地区的地下介质P波速度变化，并详细分析了P波速度随时间和空间的演化规律，进而探讨引起这种变化的原因.
1 方法和数据 1.1 基于双差层析成像的时移层析成像方法双差层析成像方法通过利用绝对到时和相对到时同时得到地震速度模型和地震事件定位，与常规的地震层析成像方法相比，双差层析成像方法能得到更加精确的源区域的速度结构(Zhang and Thurber, 2003，2006)，双差层析成像方法能表示为
|
(1) |
|
(2) |
其中，(1)式是对于绝对到时数据，(2)式是对于相对到时数据.rki是事件i和台站k的到时差，Δτi是事件i的发震时刻扰动量，xli(l=1, 2, 3)是三分量的位置扰动，δu是慢度扰动.对于双差层析成像方法，事件对是由一定距离范围内的事件构建的，并且要求至少有一定数量的相同台站记录到这两个事件.
在双差层析成像方法中，用来构建事件对的事件均来自同一时间段，即(2)式中的事件i和事件j来自同一时间段，为了求解两个时间段之间的速度变化，则要求事件i和事件j属于两个不同的时间段，即事件i发生在时间段1内，事件j发生在时间段2内，我们将(2)式改变为
|
(3) |
(2) 式和(3)式的不同点在于慢度参数δu1和δu2分别是时间段1和时间段2的慢度参数，并且事件i和事件j分属于不同的时间段，这里为属于不同时间段的事件构建的到时差.
对于基于双差层析成像方法的时移层析成像方法，首先利用双差层析成像方法使用时间段1的数据来得到3D的速度模型，然后将所得到的3D速度模型作为初始模型使用(3)式反演时间段2相对于时间段1的速度变化δu2.在这个过程中，仅使用不同时间段的事件对构建到时差，时间段1的慢度u1保持固定，到时差用来反演时间段2相对于时间段1的慢度变化δu2，同时也反演两个时间段的事件位置.在计算程序中采用伪弯曲射线追踪方法(Um and Thurber, 1987)发现射线并计算走时，选取两个时间段的地震匹配相对到时，对来自时间段1中的事件i进行射线追踪时用的速度模型是初始速度模型，来自时间段2中的事件j进行射线追踪时用的是每次迭代反演更新后的速度模型.模型表示为固定的三维节点，速度值用线性插值方法来进行插值(Thurber，1983).
新的时移层析成像方法不仅如传统层析成像方法一样能提供高空间分辨率，同时能够改进传统层析成像方法的一些缺陷，即能克服部分的观察误差的影响和不同时期地震射线分布的差异的影响，得出较为可靠的不同时间段之间的速度变化.因此这种新方法可以确定可靠的速度随时间变化的空间分布.
1.2 数据为了更好的监测余震的发生，保证余震序列数据的完整性，云南省地震局在鲁甸地震震源区附近架设了3个流动台，并加入了昭通市小孔径测震台网的2个台站以及中国地震局地球物理研究所架设的巧家地震台阵的5个台站，加上原有的云南省区域地震台网以及四川和贵州区域台网的共享台站，使得鲁甸地震震源区及其附近的地震监控能力大大提升.在鲁甸地震发生前，该区域的地震主要由云南区域台网和四川、贵州区域台网的共享台站所记录，台站的分布会随着台网规划和大地震的发生而产生变化.
本研究所使用的数据是由云南区域数字地震台网所记录到的地震事件的观测报告数据，由于鲁甸地震的余震分布集中在一个小区域内，射线分布少且不均匀，为此我们增加了余震区以外区域的地震来提高地震的射线覆盖度.选取的研究区域为鲁甸地震余震区及其附近(26.6—27.6°N，103—104°E)，时间段为2008年1月1日至2017年12月31日，我们选取了由42个台站(图 1)所记录到的5623个地震事件(图 2)，每个事件由震中距小于200 km的5个以上台站所记录.一般来说，5个以上台站所记录到的地震事件的震相都较为清晰，拾取精度小于0.1 s.
|
图 1 2014鲁甸MS6.5地震震源区附近的台站分布黑框表示研究区. Fig. 1 The distribution of stations near Ludian 2014 MS6.5 earthquake source region Black frame denote research area. |
为了确定鲁甸地震震源区及其附近的震后速度变化，我们根据数据集划分了五个地震数目大致相同的时间段(表 1)，每个时段的P波射线数目大致相同.其中P1是鲁甸地震发生前，P2是同震期，P3是震后1—4个月，P4和P5分别是震后约1年和约3年的两个时间段，这两个时间段的地震数目大致相同.为了检查观测报告中震相数据的可靠性，绘制了五个时间段的P波震相走时曲线，并且使用所选初始一维速度模型(表 2)绘制了理论走时曲线，如图 3所示.从图中可看出P波走时在所选取的距离范围内呈现出明确的线性趋势，并且震相的离散度较小，表明观测震相数据的可靠性较高.观测震相走时与理论走时曲线吻合较好，说明初始一维速度模型精度较高.图 4为五个时段内的P波二维平面射线分布图，可看出五个时段内的鲁甸地震震源区的射线密度分布均较密集.
|
表 1 五个时段的数据 Table 1 Data of the five periods |
|
|
表 2 一维P波速度模型表 Table 2 1-D P wave velocity model |
|
图 3 五个时段的P波震相走时 Fig. 3 P phase travel-time curve of the five periods |
|
图 4 五个时段P波二维射线分布 红色点表示地震，黑色三角形表示台站. Fig. 4 Distribution of 2-D P wave ray paths for the five periods Red dots denote earthquakes, black triangles denote stations. |
由时移层析成像方法的技术路线，在划分完时间段后，首先使用双差层析成像方法对每个时段的数据进行反演，以得到每个时段的三维P波速度结构.首先进行地震对匹配，选择地震对之间的最大距离为20 km，每个地震最多可与10个地震组成地震对，最终构建了相对到时数据，如表 1所示.
在离散化的模型空间中进行反演，选取(27.1°N，103.4°E)为坐标原点，每个时段均划分了相同的网格.综合考虑每个时段内的地震分布划分网格，在反演之前进行了大量的分辨率测试，以寻求最佳的网格分布.分别测试了3 km×3 km、5 km×5 km、10 km×10 km的网格间隔.由于地震分布不均匀，在地震分布密集的鲁甸地震震源区的横向分辨率达3 km×3 km，震源区外围区域地震较少，横向分辨率可达5 km×5 km，研究区边缘地震分布很少，分辨率为10 km×10 km.如图 2所示，垂直向网格位于0、2、5、7、10、16、24、30、46 km深度.准确的参考地壳模型是获得稳定结果的关键，所以在本研究中结合多个地壳模型(陈飞，2017；王光明等，2018；房立华等，2014；Riaz et al., 2017)建立了初始一维速度模型，其中陈飞(2017)是通过使用重力数据与面波数据联合反演得到的全国三维速度结构，王光明等(2018)和房立华等(2014)均是对鲁甸地震及其余震进行了重定位，Riaz等(2017)则是对鲁甸地震震源区进行了双差层析成像.在最后得到的模型中，Moho面的深度范围在30~46 km，这与郑勇等(2013)通过背景噪声成像得到的结果相似，每个时间段均采用相同的初始一维速度模型，如表 2所示.
|
图 2 研究区域活动断裂、地震分布及网格划分 红点表示地震，“叉”字表示网格节点，黑线表示断裂. Fig. 2 Geological tectonic settings, seismic stations used in this study and grid nodes division in the studied area Red dots denote earthquakes, crisscrosses denote grid nodes, black lines denote faults. |
在反演中，为了反演出较小尺度的异常并能找到地质上可行特征的反演模型而系统地测试了不同组合的平滑参数(η)和阻尼参数(ε)，最后设置了一个较小的平滑参数(η=6)，然后测试了不同的ε，并选择一个可以确保条件数在150至200之间的ε，保持条件数的数值范围较低但不太窄是很重要的，因为真实的条件数与使用LSQR所计算的条件数不同，可能更小(Zhang and Thurber, 2007).所以最后5个时间段的ε选取会有一定差别.对5个时段的数据均进行了16次迭代以确保最终的结果能够达到令人满意的收敛，前6次迭代，重点放在拟合绝对数据上，包括3个速度与地震位置的联合反演及3个重定位的反演.而在后10次迭代中，重点放在到时差数据上，包括7个速度与地震位置的联合反演及3个重定位的反演.反演后5个时段到时残差的均方根都有显著的降低，如表 3所示，表明了所获得的最终模型的稳定性.
|
|
表 3 五个时段反演前后观测到时差与理论到时差之间的均方根残差变化 Table 3 The RMS residuals between observed and predicted differential travel times based on 1D model and 3D model for the five periods |
同时在本文中对5个时段的反演结果的分辨率进行棋盘测试，在初始一维速度模型上分别增加了±5%的正负相间速度扰动.如图 5所示，在此选取了5、10、16 km深度剖面的分辨率测试结果，可以看出P2~P5这4个时段内的三个深度的鲁甸地震震源区的P波速度均恢复较好，而P1时间段由于射线数目较少，在鲁甸地震震源区恢复的相对较差.为了进一步量化模型分辨率，我们使用了衍生加权总和(DWS，Thurber and Eberhart-Philips, 1999)对模型分辨率进行补充.DWS值代表给定网格节点附近所有射线长度(乘以数据权重)的总和，一般来说，DWS值是解释可接受分辨率区域的有效工具.总的来说，DWS=100(图 6中的白色等值线)包围的基本是棋盘测试恢复良好的区域(图 5).图 6为最后反演出的5个时间段的5、10、16 km深度的P波速度分布图，可看出在5个时段的不同深度处均能较好的表现出P波速度.
|
图 5 五个时段不同深度处的P波速度棋盘测试结果 (a) 时段1；(b) 时段2；(c) 时段3；(d) 时段4；(e) 时段5. Fig. 5 Checkerboard resolution test of P wave velocity for the five periods (a) P1; (b) P2; (c) P3; (d) P4; (e) P5. |
|
图 6 五个时段不同深度处的P波速度和地震分布 黑点表示地震，白色五角星表示主震，白色实线表示DWS值大于100的区域. (a) 时段1；(b) 时段2；(c) 时段3；(d) 时段4；(e) 时段5. Fig. 6 Distribution of P-wave velocity and earthquakes at different depths for the five periods Black dots denote earthquakes, white five-pointed stars denote the main shock, white lines denote the areas that their DWS value are greater than 100. (a) P1; (b) P2; (c) P3; (d) P4; (e) P5. |
在使用双差层析成像方法计算出5个时段的P波速度模型后，接着用时移层析成像方法得到相邻时间段的P波速度变化.在这里首先构建两个时段的事件对的到时差，为了尽量多的构建分属于两个时段的两个事件的事件对，要求地震对之间的最大距离为50 km，一个时间段内的每个事件最多可以与30个事件相连接.从P1到P2，构建了平均空间距离为7 km的616007个到时差，挑选出了70470个分属于两个时段的事件对的到时差；从P2到P3，构建了平均空间距离为4.2 km的667502个到时差，最后挑选出250680个符合要求的到时差；从P3到P4，我们构建了平均空间距离为4.8 km的462766个到时差，最后挑选出符合要求的到时差164765个；从P4到P5，我们构建了平均空间距离为7.3 km的665468个到时差，最后挑选出了206155个符合要求的到时差.
对于时移层析成像方法，使用由双差层析成像方法所得到的第一个时间段的速度模型作为初始模型，重定位结果作为初始位置，然后使用两个相邻时间段所构建的到时差来得到速度随时间的变化.为了确保最后的结果能很好的收敛，反演过程经过了32次迭代，反演后到时残差的均方根变化如表 4所示，4个到时残差的均方根均有显著的下降，说明最后得到的速度变化与两个相邻时间段内的事件非常合适.
|
|
表 4 不同时段利用时移层析成像方法反演前后观测到时差与理论到时差之间的均方根残差变化 Table 4 The RMS residuals between observed and predicted differential travel times based on initial 3D model and final 3D model for the different periods |
同时还进行了棋盘测试，与3D层析成像的棋盘测试相类似，在第一个时间段的初始3D模型中加入了±5%的正负相间速度扰动，图 7显示了4个相邻时段棋盘测试的恢复结果.P1与P2，地震分布和所使用台站差异较大，虽然尽量选取多的地震对，但是仍然过于分散，所以分辨率恢复的并不好，这在一定程度上会影响反演结果的精度.其余时段由于地震分布和所使用台站差异不大，所以分辨率相对较好.但是棋盘测试恢复出来的模型分辨率与反演实际数据获得的速度模型之间没有完全的相互映射，这是因为这两个模型具有不同的速度(Share et al., 2019)，所以同时也使用了DWS值(图 7中的白色虚线)来解释可接受分辨率的区域，如图 7所示.由于在鲁甸地震震源区有相对多数量的射线覆盖，该区域棋盘恢复的相对较好，图中也显示了DWS值大于300的区域基本覆盖了震源区，所以接下来主要讨论鲁甸地震震源区的P波速度变化.
|
图 7 不同时段之间利用时移层析成像方法得到的不同深度处的棋盘测试结果 白色虚线表示DWS值大于300的区域. (a) 时段2相对于时段1；(b) 时段3相对于时段2；(c) 时段4相对于时段3；(d) 时段5相对于时段4. Fig. 7 Recovered checkerboard models for different periods by using time-lapse tomography White dotted lines denote the area that their DWS values are greater than 300. (a) P2-P1; (b) P3-P2; (c) P4-P3; (d) P5-P4. |
根据分辨率和DWS值的分布情况主要讨论鲁甸地震震源区P波速度变化，所以在图 8中主要显示了分辨率较好的鲁甸地震震源区在5、10、16 km深度处的不同时间段之间的P波速度变化.由图 8可以清楚的看出鲁甸地震发生前后的P波速度的降低及恢复.其中图 8a是由P1和P2两个时段的数据得到的鲁甸地震同震速度变化，可以看出每个深度处的震源区的速度均下降，下降的幅度不大，约在0.2%以下，大多为百分之零点零几，其中16 km深度处的速度下降幅度小于5 km和10 km深度处.图 8b是P2和P3两个时段的数据得到的震后1—4个月相对同震期的速度变化，可看出震源区的速度大幅下降，下降幅度大多在0.2%~3%内，这个范围与前人对相同量级地震发生前后速度变化研究结果相一致.例如Li等(1998)计算的2004年Parkfield M6地震发生后速度下降了2.5%；Nishimura等(1998)计算的1998年Intra-Plate M6.1地震发生后速度下降了0.3%~10%.相比P2-P1，P3-P2的速度下降的幅度有很大的变化，达到了最大的速度降幅，说明鲁甸地震发生后速度并没有立即下降到最大，而是在震后1—4个月期间速度下降到了最大.这与庞卫东等(2017)利用背景噪声自相关方法得到的鲁甸地震发生后最大波速降出现在震后1个月的时间相一致，也与刘志坤和黄金莉(2010)计算的汶川地震、王俊等(2020)计算的芦山地震发生后的最大波速降出现在震后1—4个月内相符合.
|
图 8 不同深度处的不同时段之间的P波速度变化和地震分布 (a) 时段2相对于时段1；(b) 时段3相对于时段2；(c) 时段4相对于时段3；(d) 时段5相对于时段4. Fig. 8 Distribution of temporal P-wave velocity changes and earthquakes at different depths for different periods (a) P2-P1; (b) P3-P2; (c) P4-P3; (d) P5-P4. |
在震后速度下降达到最大值后，震源区开始愈合，速度开始上升，由P4-P3和P5-P4的速度变化可以清楚的看出震源区在时空上的愈合过程.愈合从浅层开始，首先是5 km深度处的震源区的速度开始上升，上升幅度大约为0.1%以下，而10 km和16 km深度处的速度仍然是在下降，说明浅层区域开始愈合.随着时间发展，愈合延伸至更深的深度，5 km和10 km深度处的震源区的速度也上升，16 km深度处的震源区的速度虽然在下降，但是有一部分区域的速度在上升，说明愈合在16 km深度处也在进行，上升的幅度约在0.2%~2%，说明基本在向震前的水平进行愈合.
总的来看，2014年鲁甸地震发生后，震源区的速度变化特征为，同震期速度下降，但没有下降至最大，而是在震后1—4个月内速度下降至最大，然后愈合过程开始，愈合从浅层逐渐发展至深层.值得注意的一点是，由于我们的观测同时对应于空间和时间的平均值，因此不应过度解释所有变化幅度的数值大小.
由于鲁甸地震余震分布呈共轭型分布，所以我们还做了如图 2所示的两个剖面的速度变化示意图(图 9)，其中AA′剖面主要沿着地震序列的SES方向穿过震源区，BB′剖面主要沿着地震序列的近EW向穿过震源区，这样我们就能看到两个余震优势分布方向上的速度变化的时空分布.在同震期，AA′和BB′两个剖面下方均为速度下降，且靠近AA′剖面的A′端即SES端的区域和靠近BB′剖面的B′端即近E向端的区域的速度下降幅度略大，这可能是因为这两个区域是鲁甸地震发生后余震最初集中发生的区域.在震后1—4月，两个剖面下方速度下降剧烈，其中靠近A′端的区域和靠近B端的区域的速度下降最为剧烈，这是因为震后余震分布发生了变化，不再局限于主震附近，近EW向分支的余震继续往W向发展，SES向分支的余震则仍然往SES向发展(王光明等，2018)，对于余震发展的方向的解释，Riaz等(2017)和赵小艳和孙楠(2014)均认为是由于震源区的北部存在一个高速异常，对余震有阻碍作用.两个剖面的震后愈合恢复过程均是从浅层开始，AA′剖面余震分布较为集中的SES端的区域及BB′剖面余震最开始发生的靠近近E端的区域的下方至5 km深度的速度上升，随着时间发展，AA′剖面的SES端区域下方速度上升幅度增大，延伸至约8 km深度，同时BB′剖面下方速度上升的区域延伸至近W端，即后来余震发展的方向，且该区域下方速度上升幅度较大，同时速度上升的区域扩大至约15 km深度.
|
图 9 不同时段之间的P波速度变化及地震事件沿垂直剖面AA′和BB′的分布图(剖面位置见图 2) 红色虚线表示DWS大于300的区域. (a) 时段2相对于时段1；(b) 时段3相对于时段2；(c) 时段4相对于时段3；(d) 时段5相对于时段4. Fig. 9 Distribution of temporal P-wave velocity changes along the cross sections AA′, BB′ for different periods(section location is in figure 2) Red dotted lines denote the areas that their DWS values are greater than 300. (a) P2-P1; (b) P3-P2; (c) P4-P3; (d) P5-P4. |
两个剖面的速度变化的发展过程同时也反映了余震分布的发展过程.结合房立华等(2014)和王光明等(2018)对鲁甸地震余震序列的重定位研究结果，可以得出：鲁甸地震发生后余震主要分布在主震附近及其东南侧，这些区域下方的速度开始下降，下降幅度较小，此后主震西侧开始出现余震，而东南侧的余震分布范围基本稳定，余震呈现出不对称的共轭分布特征，主震西侧及东南侧这两个区域下方的速度下降剧烈，速度下降至最大后愈合过程开始，余震最初出现的主震附近及其东南侧区域下方的浅层的速度上升，随着时间的发展，愈合延伸至余震后来发展的区域，即主震西侧，且该区域下方速度上升幅度增大，整个近EW向分支区域下方速度上升的深度扩展至约15 km深度处，而余震分布稳定的主震东南侧区域下方速度上升幅度增大，范围扩大至约8 km.
4 讨论一般认为，引起地震发生前后波速变化的物理机制有以下几种：地震断层区的破坏与愈合(Li et al., 1998；Vidale and Li, 2003)，强地面运动引起的地表破坏(Peng and Ben-Zion, 2006)，断层区及周边地壳介质的应力变化(Wegler et al., 2009；Brenguier et al., 2008；Niu et al., 2008；Pei et al., 2019)，地下水位变化(Sens-Schönfelder and Wegler, 2006)等(刘志坤和黄金莉，2010).而要准确地判定引起地震波速变化的物理机制是非常困难的.鲁甸地震发生后，对鲁甸地震震源运动学、余震分布特征、应力变化以及地质构造等方面的研究有很多，可以将本研究所得到的结果与已有的地质与地球物理研究成果相结合，从而进一步分析造成速度变化的可能机制.
一般来说，速度变化能反映强地面运动所造成的岩石损伤，预计主要发生在顶部几百米(Pei et al., 2019)，而本研究在中上地壳均发现了P波速度的明显变化，因此大地震的强地面运动引起的浅部地壳介质破坏对本研究得到的波速变化结果影响不大，不是主要因素.
鲁甸地震发生后，在由中国地震局组织的野外应急科学考察中并没有发现明显的地震地表破裂现象(徐锡伟等，2014)，但是李西等(2018)基于详细的野外地质地貌调查发现了自龙头山镇以南谢家营盘—光明村—王家坡一段长约8 km断续展布的地表破裂带，呈左行右阶排列，伴有30~35 cm走滑分量，局部段落伴有15~25 cm倾滑分量.这段地表破裂带主要位于鲁甸地震余震带SES向分支上，而主震附近及近EW向分支上并没有观察到地表破裂，而在本研究中整个震源区均有明显的波速变化，表明本研究得到的波速变化对整个震源区的介质性质变化都有所反映.并且地表破裂带应该仅对地表浅层区域有较大影响，对深层区域的介质影响不大，而本研究得到的波速变化从浅层至深层均有明显变化，所以判断地震断层区的破坏并不是影响本研究所得到的波速变化的主要原因，但是对浅层区域及主震东南方向区域介质变化具有一定的影响.
破裂区深度内的结构变化和开放裂缝导致应力场的变化可能发生在不同的深度(Pei et al., 2019)，通过观察图 8发现同震期的5 km和10 km深度处及震后1—4个月的10 km和16 km深度处均呈现出按象限分布的现象，且在以前的很多研究中(Wegler et al., 2009；Cheng et al., 2010；刘志坤和黄金莉，2010；王俊等，2020)均发现了这一现象.为了探究波速变化呈象限分布是否是由地震的应力变化引起的，查找了相关的研究结果.缪淼和朱守彪(2016)计算了鲁甸地震静态库伦应力变化，发现两个方向的余震分布都主要集中在主震产生的库伦应力增大区，付芮等(2015)也发现余震主要位于沿破裂面及主破裂以西的库仑应力增强区.对比图 8发现5 km和10 km深度处及震后1—4个月的10 km及16 km深度处P波速度下降的区域基本为余震主要分布区，说明P波速度下降区主要对应于应力增大区，表明鲁甸地震主震产生的应力变化与P波速度变化之间有紧密的联系，这也从另一角度反映出，地震引起的应力变化是造成鲁甸地震速度变化的一个重要原因.并且在第3节内容中提到鲁甸地震发生后最大速度下降发生在震后1—4个月，并没有发生在同震期，而速度的下降及上升在空间上的分布变化均与余震发生的空间变化有关，这表明地震波速变化与余震的时空分布有关，可能是余震的动态和静态应力造成震源区介质物理属性的改变从而影响了地震波速的变化.
鲁甸地震发生在南北地震带南段，青藏高原东南缘鲜水河—小江断裂系以东走向NE向昭通—莲峰断裂带西段，北临NNW至近NS向大凉山断裂和马边—盐津断裂，属青藏高原东南部川滇菱形块体与华南地块之间I级块体边界构造带(徐锡伟等，2014).庞卫东等(2017)认为鲁甸地震是青藏高原东南缘与华南地块相互作用而造成地壳内部构造应力变化的结果；程佳等(2016)认为鲁甸地震受到了主要来自于则木河断裂带上的应力作用并积累地震能量；洪敏利用应变场时序分析法(洪敏等，2014)计算出的鲁甸地震发生前的云南地区应变场分布发现滇东北地区出现明显的挤压活动增强1).所以区域应力场的变化对鲁甸地震的发生有着很大的影响，而鲁甸地震对后续余震具有触发作用，从而对震源区介质的物理性质变化有着一定的影响.
1) 云南省地震局形变测量中心2018年年度会商报告.
5 结论运用基于双差层析成像的时移层析成像方法能够得到高精度的地震波速变化的时空分布.在本研究中我们利用云南区域数字地震台网所记录到的地震资料，使用该方法来得到了2014年鲁甸MS6.5地震发生前后的P波速度变化的时空分布，结果表明：2014年鲁甸地震发生后，同震期的震源区速度下降，但没有下降至最大，而是在震后1—4个月内速度下降至最大，接着P波速度开始上升，开始愈合，愈合过程从浅层逐渐发展至深层，并逐渐恢复至震前水平.由于鲁甸地震余震区呈共轭形态，所以我们做了分别穿过两个余震分支的垂直剖面，发现P波速度变化与余震分布变化相一致，余震最初分布的主震附近及其东南侧区域下方的速度最先下降，然后主震西侧开始出现余震，则主震西侧及东南侧这两个区域下方的速度下降剧烈，愈合过程则是余震最初出现的区域下方开始，之后延伸至余震后来发展的区域.
进一步将所得到的结果和已有的地质与地球物理观测相结合来分析引起这种变化的可能机制，认为强地面运动引起的地表破坏及地震断层区的破坏并不是影响本研究所得到的波速变化的主要原因，但是对浅层区域介质变化具有一定的影响.由于在空间上鲁甸地震应力变化与波速变化之间有较紧密的联系，且波速变化与余震分布变化相一致，所以我们认为地震引起的应力变化是造成鲁甸地震波速变化的一个重要原因，余震的动态和静态应力会造成震源区介质物理属性的改变从而影响了地震波速的变化.并且由以往的研究也可得出区域构造应力场变化对鲁甸地震的发生有很大作用，从而对震源区地下物理性质的改变造成影响.
致谢 感谢中国科技大学地球和空间科学学院胡景博士在本研究过程中提供的帮助，感谢云南省地震局的王光明助理研究员和洪敏高工对文章的热心帮助，感谢三位匿名审稿专家提出的宝贵意见和建议.
Brenguier F, Campillo M, Hadziioannou C, et al. 2008. Postseismic relaxation along the San Andreas Fault at Parkfield from continuous seismological observations. Science, 321(5895): 1478-1481. DOI:10.1126/science.1160943 |
Chang L J, Ding Z F, Wang C Y. 2010. Variations of shear wave splitting in the 2010 Yushu MS7.1 earthquake region. Chinese Journal of Geophysics (in Chinese), 53(11): 2613-2619. DOI:10.3969/j.issn.0001-5733.2010.11.009 |
Chen F. 2017. Lithospheric shear wave tomography of continental China by joint inversion of surface-wave and satellite gravity data[Master's thesis] (in Chinese). Hefei: University of Science and Technology of China, 30-38.
|
Cheng J, Xu X W, Liu J. 2016. Cause and rupture characteristics of the 2014 Ludian MS6.5 mainshock and its aftershock distribution using the coulomb stress changes. Chinese Journal of Geophysics (in Chinese), 59(2): 556-567. DOI:10.6038/cjg20160214 |
Cheng X, Niu F L, Wang B S. 2010. Coseismic velocity change in the rupture zone of the 2008 MW7.9 Wenchuan Earthquake observed from ambient seismic noise. Bulletin of the Seismological Society of America, 100(5B): 2539-2550. DOI:10.1785/0120090329 |
Chiarabba C, De Gori P, Boschi E. 2009. Pore-pressure migration along a normal-fault system resolved by time-repeated seismic tomography. Geology, 37(1): 67-70. DOI:10.1130/G25220A.1 |
Crampin S, Booth D C, Evans R, et al. 1990. Changes in shear-wave splitting at Anza near the time of the North Palm Springs earthquake. Journal of geophysical Research: Solid Earth, 95(B7): 11197-11212. DOI:10.1029/JB095iB07p11197 |
Dai Z H. 2016. Study on changes of seismic velocity associated with strong earthquakes in Sichuan and Yunnan province from seismic ambient noise[Master's thesis] (in Chinese). Beijing: Institute of Earthquake Science, China Earthquake Administration, 51-56.
|
Fang L H, Wu J P, Wang W L, et al. 2014. Relocation of the aftershock sequence of the MS6.5 Ludian earthquake and its seismogenic structure. Seismology and Geology (in Chinese), 36(4): 1173-1185. |
Fu R, Shan B, Xiong X, et al. 2015. Co-seismic coulomb stress changes and its influences on aftershock distribution and surrounding faults caused by 2014 Ludian earthquake, Yunnan. Seismology and Geology (in Chinese), 37(4): 1084-1095. |
Gao Y, Wu J. 2008. Compressive stress field in the crust deduced from shear-wave anisotropy: an example in capital area of China. Chinese Science Bulletin, 53(18): 2840-2848. |
Hong M, Shao D S, Wang L L, et al. 2014. Using spatial correlation between continuous GNSS time series to achieve regional deformation anomalies. Journal of Geodesy and Geodynamics (in Chinese), 34(2): 175-178, 182. |
Julian B R, Foulger G R. 2010. Time-dependent seismic tomography. Geophysical Journal International, 182(3): 1327-1338. DOI:10.1111/j.1365-246X.2010.04668.x |
Li X, Zhang J G, Xie Y Q, et al. 2014. Ludian MS6.5 earthquake surface damage and its relationship with structure. Seismology and Geology (in Chinese), 36(4): 1280-1291. |
Li X, Xu X W, Zhang J G, et al. 2018. Surface rupture characteristics of the seismogenic structure of the Ludian MS6.5 earthquake and identification and dating of related paleoearthquakes. Earth Science Frontiers (in Chinese), 25(1): 227-239. |
Li Y G, Vidale J E, Aki K, et al. 1998. Evidence of shallow fault zone strengthening after the 1992 M7.5 Landers, California, earthquake. Science, 279(5348): 217-219. DOI:10.1126/science.279.5348.217 |
Liu Z K, Huang J L. 2010. Temporal changes of seismic velocity around the Wenchuan earthquake fault zone from ambient seismic noise correlation. Chinese Journal of Geophysics (in Chinese), 53(4): 853-863. DOI:10.3969/j.issn.0001-5733.2010.04.010 |
Miao M, Zhu S B. 2016. The static coulomb stress change of the 2014 Ludian earthquake and its influence on the aftershocks and surrounding faults. Seismology and Geology (in Chinese), 38(1): 169-181. |
Nishimura T, Uchida N, Sato H, et al. 2000. Temporal changes of the crustal structure associated with the M6.1 earthquake on September 3, 1998, and the volcanic activity of Mount Iwate, Japan. Geophysical Research Letters, 27(2): 269-272. DOI:10.1029/1999GL005439 |
Niu F L, Silver P G, Daley T M, et al. 2008. Preseismic velocity changes observed from active source monitoring at the Parkfield SAFOD drill site. Nature, 454(7201): 204-208. DOI:10.1038/nature07111 |
Obermann A, Froment B, Campillo M, et al. 2014. Seismic noise correlations to image structural and mechanical changes associated with the MW7.9 2008 Wenchuan earthquake. Journal of Geophysical Research: Solid Earth, 119(4): 3155-3168. DOI:10.1002/2013JB010932 |
Pang W D, Yang R H, Xu Y J, et al. 2017. Study on the crustal wave velocity variation around the 2014 Ludian MS6.5 earthquake area by background noise autocorrelation function. Journal of Seismological Research (in Chinese), 40(3): 491-501. |
Patanè D, Barberi G, Cocina O, et al. 2006. Time-resolved seismic tomography detects magma intrusions at Mount Etna. Science, 313(5788): 821-823. DOI:10.1126/science.1127724 |
Pei S P, Niu F L, Ben-Zion Y, et al. 2019. Seismic velocity reduction and accelerated recovery due to earthquakes on the Longmenshan fault. Nature Geoscience, 12(5): 387-392. DOI:10.1038/s41561-019-0347-1 |
Peng Z G, Ben-Zion Y. 2006. Temporal changes of shallow seismic velocity around the Karadere-Duzce branch of the north Anatolian fault and strong ground motion. Pure and Applied Geophysics, 163(2-3): 567-600. DOI:10.1007/s00024-005-0034-6 |
Poupinet G, Ellsworth W L, Frechet J. 1984. Monitoring velocity variations in the crust using earthquake doublets: An application to the Calaveras Fault, California. Journal of Geophysical Research: Solid Earth, 89(B7): 5719-5731. DOI:10.1029/JB089iB07p05719 |
Qian J W, Zhang H J, Westman E. 2018. New time-lapse seismic tomographic scheme based on double-difference tomography and its application in monitoring temporal velocity variations caused by underground coal mining. Geophysical Journal International, 215(3): 2093-2104. DOI:10.1093/gji/ggy404 |
Riaz M S, Zheng Y, Xiong X, et al. 2017. Refined 3D seismic-velocity structures and seismogenic environment of the MS6.5 Ludian Earthquake. Bulletin of the Seismological Society of America, 107(6): 3023-3036. DOI:10.1785/0120170072 |
Schaff D P, Beroza G C. 2004. Coseismic and postseismic velocity changes measured by repeating earthquakes. Journal of Geophysical Research: Solid Earth, 109(B10): B10302. DOI:10.1029/2004JB003011 |
Sens-Schönfelder C, Wegler U. 2006. Passive image interferometry and seasonal variations of seismic velocities at Merapi Volcano, Indonesia. Geophysical Research Letters, 33(21): L21302. DOI:10.1029/2006GL027797 |
Share P E, Guo H, Thurber C H, et al. 2019. Seismic imaging of the Southern California plate boundary around the South-Central transverse ranges using double-difference tomography. Pure and Applied Geophysics, 176(3): 1117-1143. DOI:10.1007/s00024-018-2042-3 |
Taira T, Silver P G, Niu F L, et al. 2008. Detecting seismogenic stress evolution and constraining fault zone rheology in the San Andreas Fault following the 2004 Parkfield earthquake. Journal of Geophysical Research: Solid Earth, 113(B3): B03303. DOI:10.1029/2007JB005151 |
Thurber C H. 1983. Earthquake locations and three-dimensional crustal structure in the Coyote Lake area, central California. Journal of Geophysical Research: Solid Earth, 88(B10): 8226-8236. DOI:10.1029/JB088iB10p08226 |
Thurber C H, Eberhart-Phillips D. 1999. Local earthquake tomography with flexible gridding. Computers & Geosciences, 25(7): 809-818. |
Um J, Thurber C H. 1987. A fast algorithm for two-point seismic ray tracing. Bulletin of the Seismological Society of America, 77(3): 972-986. |
Vidale J E, Li Y G. 2003. Damage to the shallow Landers fault from the nearby Hector Mine earthquake. Nature, 421(6922): 524-526. DOI:10.1038/nature01354 |
Wang G M, Zhao X Y, Fu H, et al. 2018. Relocation of the 2014 Ludian MS6.5 earthquake and its large aftershock sequence. Acta Seismologica Sinica (in Chinese), 40(5): 582-594. DOI:10.11939/jass.20170180 |
Wang J, Zheng D C, Zhang J C, et al. 2020. Seismic velocity changes in the epicentral region of the 2013 Lushan earthquake measured from ambient seismic noise. Chinese Journal of Geophysics (in Chinese), 63(2): 517-531. DOI:10.6038/cjg2020N0275 |
Wegler U, Nakahara H, Sens-Schönfelder C, et al. 2009. Sudden drop of seismic velocity after the 2004 MW6.6 Mid-Niigata Earthquake, Japan, observed with passive image interferometry. Journal of Geophysical Research: Solid Earth, 114(B6): B06305. DOI:10.1029/2008JB005869 |
Westman E C, Luxbacher K D, Schafrik S J, et al. 2012. Time-lapse passive seismic velocity tomography of longwall coal mines: a comparison of methods. //82nd Ann. Internat Mtg., Soc. Expi. Geophys. . Expanded Abstracts, 602-607.
|
Xu L S, Zhang X, Yang C, et al. 2014. Analysis of the Love waves for the source complexity of the Ludian MS6.5 earthquake. Chinese Journal of Geophysics (in Chinese), 57(9): 3006-3017. DOI:10.6038/cjg20140925 |
Xu X W, Jiang G Y, Yu G H, et al. 2014. Discussion on seismogenic fault of the Ludian MS6.5 earthquake and its tectonic attribution. Chinese Journal of Geophysics (in Chinese), 57(9): 3060-3068. DOI:10.6038/cjg20140931 |
Zhang G W, Lei J S, Liang S S, et al. 2014. Relocations and mechanism solutions of the 3 August 2014 Ludian, Yunnan MS6.5 earthquake sequence. Chinese Journal of Geophysics (in Chinese), 57(9): 3018-3027. DOI:10.6038/cjg20140926 |
Zhang H, Thurber C H. 2007. Estimating the model resolution matrix for large seismic tomography problems based on Lanczos bidiagonalization with partial reorthogonalization. Geophysical Journal International, 170: 337-345. DOI:10.1111/j.1365-246X.2007.03418.x |
Zhang H J, Thurber C H. 2003. Double-difference tomography: The method and its application to the Hayward fault, California. Bulletin of the Seismological Society of America, 93(5): 1875-1889. DOI:10.1785/0120020190 |
Zhang H J, Thurber C H. 2006. Development and applications of double-difference seismic tomography. Pure and Applied Geophysics, 163(2-3): 373-403. DOI:10.1007/s00024-005-0021-y |
Zhang Y, Xu L S, Chen Y T, et al. 2014. Rupture process of the 3 August 2014 Ludian, Yunnan, MW6.1 (MS6.5) earthquake. Chinese Journal of Geophysics (in Chinese), 57(9): 3052-3059. DOI:10.6038/cjg20140930 |
Zhang Y, Chen Y T, Xu L S, et al. 2015. The 2014 MW6.1 Ludian, Yunnan, earthquake: A complex conjugated ruptured earthquake. Chinese Journal of Geophysics (in Chinese), 58(1): 153-162. DOI:10.6038/cjg20150113 |
Zhao P, Peng Z G. 2009. Depth extent of damage zones around the central Calaveras fault from waveform analysis of repeating earthquakes. Geophysical Journal International, 179(3): 1817-1830. DOI:10.1111/j.1365-246X.2009.04385.x |
Zhao X Y, Sun N. 2014. Simultaneous inversion for focal location of Yunnan Ludian MS6.5 earthquake sequence in 2014 and velocity structure in the source region. Journal of Seismological Research (in Chinese), 37(4): 523-531. |
Zheng Y, Ge C, Xie Z J, et al. 2013. Crustal and upper mantle structure and the deep seismogenic environment in the source regions of the Lushan earthquake and the Wenchuan earthquake. Science China Earth Sciences, 56(7): 1158-1168. DOI:10.1007/s11430-013-4641-2 |
常利军, 丁志峰, 王椿墉. 2010. 2010年玉树7.1级地震震源区横波分裂的变化特征. 地球物理学报, 53(11): 2613-2619. DOI:10.3969/j.issn.0001-5733.2010.11.009 |
陈飞. 2017. 利用面波与重力的联合反演确定中国大陆三维岩石圈速度结构[硕士论文]. 合肥: 中国科学技术大学, 30-38.
|
程佳, 徐锡伟, 刘杰. 2016. 2014年鲁甸6.5级地震成因、破裂特征和余震分布特征的库仑应力作用. 地球物理学报, 59(2): 556-567. DOI:10.6038/cjg20160214 |
戴宗辉. 2016. 利用地震背景噪声研究川滇地区强震前后波速变化特征[硕士论文]. 北京: 中国地震局地震预测研究所, 51-56.
|
房立华, 吴建平, 王未来, 等. 2014. 云南鲁甸MS6.5地震余震重定位及其发震构造. 地震地质, 36(4): 1173-1185. DOI:10.3969/j.issn.0253-4967.2014.04.019 |
付芮, 单斌, 熊熊, 等. 2015. 2014年云南鲁甸地震同震库仑应力对余震分布及周边断层的影响. 地震地质, 37(4): 1084-1095. DOI:10.3969/j.issn.0253-4967.2015.04.012 |
高原, 吴晶. 2008. 利用剪切波各向异性推断地壳主压应力场: 以首都圈地区为例. 科学通报, 55(23): 2933-2939. DOI:10.3321/j.issn:0023-074X.2008.23.015 |
洪敏, 邵德盛, 王伶俐, 等. 2014. 利用GNSS连续站时间序列空间相关性提取区域形变异常. 大地测量与地球动力学, 34(2): 175-178, 182. |
李西, 张建国, 谢英情, 等. 2014. 鲁甸MS6.5地震地表破坏及其与构造的关系. 地震地质, 36(4): 1280-1291. DOI:10.3969/j.issn.0253-4967.2014.04.026 |
李西, 徐锡伟, 张建国, 等. 2018. 鲁甸MS6.5地震发震断层地表破裂特征、相关古地震的发现和年代测定. 地学前缘, 25(1): 227-239. |
刘志坤, 黄金莉. 2010. 利用背景噪声互相关研究汶川地震震源区地震波速度变化. 地球物理学报, 53(4): 853-863. DOI:10.3969/j.issn.0001-5733.2010.04.010 |
缪淼, 朱守彪. 2016. 2014年鲁甸地震(MS=6.5)静态库仑应力变化及其影响. 地震地质, 38(1): 169-181. DOI:10.3969/j.issn.0253-4967.2016.01.013 |
庞卫东, 杨润海, 许亚吉, 等. 2017. 利用背景噪声自相关研究2014年鲁甸MS6.5地震震区周边地壳波速变化. 地震研究, 40(3): 491-501. DOI:10.3969/j.issn.1000-0666.2017.03.027 |
王光明, 赵小艳, 付虹, 等. 2018. 2014年鲁甸MS6.5地震及其强余震序列重定位. 地震学报, 40(5): 582-594. DOI:10.11939/jass.20170180 |
王俊, 郑定昌, 张金川, 等. 2020. 2013年芦山地震震源区地壳介质地震波速变化的特征分析. 地球物理学报, 63(2): 517-531. DOI:10.6038/cjg2020N0275 |
许力生, 张旭, 严川, 等. 2014. 基于勒夫波的鲁甸MS6.5地震震源复杂性分析. 地球物理学报, 57(9): 3006-3017. DOI:10.6038/cjg20140925 |
徐锡伟, 江国焰, 于贵华, 等. 2014. 鲁甸6.5级地震发震断层判定及其构造属性讨论. 地球物理学报, 57(9): 3060-3068. DOI:10.6038/cjg20140931 |
张广伟, 雷建设, 梁姗姗, 等. 2014. 2014年8月3日云南鲁甸MS6.5地震序列重定位与震源机制研究. 地球物理学报, 57(9): 3018-3027. DOI:10.6038/cjg20140926 |
张勇, 许力生, 陈运泰, 等. 2014. 2014年8月3日云南鲁甸MW6.1(MS6.5)地震破裂过程. 地球物理学报, 57(9): 3052-3059. DOI:10.6038/cjg20140930 |
张勇, 陈运泰, 许力生, 等. 2015. 2014年云南鲁甸MW6.1地震: 一次共轭破裂地震. 地球物理学报, 58(1): 153-162. DOI:10.6038/cjg20150113 |
赵小艳, 孙楠. 2014. 2014年云南鲁甸6.5级地震震源位置及震源区速度结构联合反演. 地震研究, 37(4): 523-531. DOI:10.3969/j.issn.1000-0666.2014.04.006 |
郑勇, 葛粲, 谢祖军, 等. 2013. 芦山与汶川地震震区地壳上地幔结构及深部孕震环境. 中国科学: 地球科学, 43(6): 1027-1037. |