曹 颖 钱佳威 黄江培 张国权 付 虹
1)云南省地震局,昆明 650224
2)中国科学技术大学,地球和空间科学学院,合肥 230026
据中国地震台网中心测定,2014年10月7日21点49分,云南省普洱市景谷县(23.39°N,100.46°E)发生了MS6.6地震(以下简称景谷地震),震源深度为5km。之后,同年12月6日又发生了2次震级分别为MS5.8(23.32°N,100.49°E,震源深度为9km)和MS5.9(23.33°N,100.50°E,震源深度为10km)的强余震。景谷地震发生在青藏高原南部,处于思茅-普洱地震带和耿马-澜沧地震带之间,该区域构造复杂,NNW 向和NEE向的共轭断裂体系网络相互交叉,相互制约(图1),主要表现为NNW 向展布的澜沧江断裂和NNE向展布的无量山断裂。在历史上,虽然思茅-普洱地震带和耿马-澜沧地震带均有很多大地震发生,但位于两者之间的地带则地震活动稀少,此次MS6.6地震是该地区有史以来记录到的最大地震(毛泽斌等,2019)。景谷地震发生后,针对该地震开展的研究有很多,其中重定位及震源机制解研究(徐甫坤等,2015;陈浩等,2016;李丹宁等,2017;Wangetal.,2018)认为景谷MS6.6地震为一次右旋走滑地震,余震序列主要呈2个优势方向:主震发生后沿NW 向分布;2014年12月6日MS5.8和MS5.9 2个强余震发生后沿近NNW 向分布。针对发震构造进行的研究(徐锡伟等,2014;常祖峰等,2016;吴坤罡等,2016;毛泽斌等,2019,谢张迪等,2019)则存在不同的认识,这主要是由于该地震所处区域的地质构造资料很少,没有已知断裂能与发震构造特征相匹配,故未得到统一的结果。由于景谷地区历史地震稀少,前人对该区域的关注不多,对该区域地下介质速度结构的专门研究很少,大多为大区域、大范围的讨论。其中,李永华等(2014)发现景谷地震序列下方没有出现显著的壳内低速带,但上地壳表现为S波低速异常,并认为可能与地壳强烈破碎及断层、微裂隙中的流体有关。与景谷地震有关的波速变化研究目前尚无报道,因此对地质构造复杂且发生过强震的景谷地区的地下速度结构以及与景谷地震有关的波速变化开展研究十分必要,同时掌握强震发生前后震源区介质物理性质的变化对于了解地震的孕育、发生及动力学过程也具有重要意义。
图1 研究所用台站、地震、网格节点及活动断裂分布Fig.1 Distribution of stations,earthquakes,grid nodes and faults in this study.a所用台站分布;b地震和网格节点分布。红色圆点表示地震,黑色“×”表示网格节点,黑色实线表示断裂,虚线框为研究区域。F1 澜沧江断裂;F2 窝拖寨断裂;F3 南谷断裂;F4 永平盆地东缘断裂;F5 威远江断裂;F6 益香-赵家沟断裂; F7 景谷-云仙断裂;F8 无量山断裂;F9 普文断裂。地震数据为本文重定位数据,断裂构造数据引自毛泽斌等(2019)
探测地壳介质波速变化的研究大多使用被动源数据,方法主要包括剪切波分裂法(Crampinetal.,1990)、重 复 地 震 法(Poupinetetal.,1984)、尾 波 干 涉 技 术(Sniederetal.,2002)、层析成像技术及近十几年来应用较广泛的基于背景噪声互相关法。这些方法中,地震层析成像方法是能够提供最高空间分辨率的方法之一,而进行波速变化的层析成像研究则一般使用传统的层析成像方法。例如,Foulger等(1997)使用SIMUL方法发现在1991—1994年期间,北加州间歇泉地热储层的VP/VS下降约4%;Gunasekera等(2003)也用相同的方法证实了这一观点,并指出该变化主要由VP降低所引起;Patane等(2006)使用相同的方法发现Etane火山在2002—2003年喷发期间,某些区域的VP/VS约增加4.5%。一般而言,传统层析成像方法通常是独立地对不同时期的地震波到时数据集进行反演,并假设所得到的速度差异代表真实的速度变化,这种假设是危险的,因为观测误差和地震射线分布的差异,即使地壳内的速度结构没有改变,独立反演的结果也会不同(Julianetal.,2010),更好的办法是同时反演多个数据集,这样就可以确定数据真正的变化。Julian等(2010)在此基础上提出了依赖时间的层析成像方法,该方法可以同时反演多个数据集使不同时间段速度模型之间的差异和观察到时与计算到时之间的差异达到最小,但在一定程度上仍将受到射线分布不均匀的影响。Zhang等(2015)提出了一种新的时移地震走时层析成像方法,该方法基于小波变换,可在一定程度上减少不同时间段内不同数据分布带来的误差,但仍受制于数据分布差异误差带来的影响。在以往研究的基础上,Qian等(2018)提出了一种新的时移层析成像方法,该方法基于双差层析成像方法(Zhangetal.,2003,2006),同时反演多个时间段的数据集以确定震源参数和2个时间段之间的波速变化,在反演过程中通过使用复杂的迭代策略解决了射线分布差异所产生的假性变化,为了消除模型差异的影响,将一个时间段的速度模型作为另一时间段的初始模型,因此这种新的时移层析成像可确定相对可靠的波速随时间的变化。
本文中,我们将基于双差层析成像的时移层析成像方法应用于由云南数字地震台网所记录的地震资料,得到景谷地震发生前后高精度的震源区P波速度的时空变化特征,并结合已有的地质构造、震源破裂、余震重定位和震源参数等研究成果,详细探讨引起变化的原因。
本文利用基于双差层析成像的时移层析成像方法(Qianetal.,2018)反演得到景谷地震发生前后的P波速度变化的时空演化。为了求解2个时间段之间的速度变化,将双差层析成像方法(Zhangetal.,2003,2006)进行转换,对于同一个观测台站k,分属于2个时间段的事件i和事件j的到时差之差为
式中,事件i发生在时间段1内,事件j发生在时间段2内,δu1和δu2分别是时间段1和时间段2的慢度参数,Tk是发震时刻,Δτ是事件的发震时刻扰动量,xl(l=1,2,3)是三分量的位置扰动。
对于该方法,首先利用双差层析成像方法使用时间段1的数据得到3D速度模型,然后将所得到的3D速度模型作为初始模型使用式(1)反演时间段2相对于时间段1的速度变化δu2。在这个过程中,仅使用不同时间段的事件对构建到时差,时间段1的慢度u1保持固定,到时差用来反演时间段2相对于时间段1的慢度变化δu2,同时也反演2个时间段事件的位置(Qianetal.,2018)。在该方法中采用规则的三维节点,速度值用线性插值方法进行插值(Thurber,1983)。采用伪弯曲射线追踪方法(Umetal.,1987)发现射线并计算走时,选取2个时间段的地震匹配相对到时,对来自时间段1的事件i进行射线追踪时使用的速度模型是初始速度模型,来自时间段2的事件j进行射线追踪时使用的是每次迭代反演更新后的速度模型。反演采用阻尼最小二乘分解算法(LSQR)进行求解。
本研究收集了2008年1月1日—2017年12月31日由云南省数字地震台网记录到的景谷地震震源区及其附近区域(图1)的地震观测资料。云南省数字地震台网包含固定台站、大地震发生后架设的流动台站以及水库监测台站,台站分布会因建设规划及特殊事件的发生而产生变化。为了保证数据的可靠性,经过多次严格的筛选,剔除走时曲线中离散较大的观测数据,只挑选至少由5个震中距在200km以内的台站所记录到的地震事件,最终挑选出了由39个台站(图1a)记录的ML≥0.0地震事件9 674个(图1b),所涉及台站包含30个固定台站、4个景谷地震应急流动台站以及5个糯扎渡水库台网台站。由P波走时曲线(图2)可见震相走时的离散度小,表明原始震相观测报告具有较高的可靠性。
图2 P波走时曲线Fig.2 Curve of P-wave travel times.
为了得到景谷地震发生前后的P波速度变化,首先根据景谷地震的余震发生情况将景谷地震前后划分为几个时间段,以展现景谷地震前后的波速变化的时空演化。根据前人的研究结果(李丹宁等,2017;Wangetal.,2018)可知景谷地震的余震序列在时间上分为2个阶段:第1个阶段是2014年10月7日MS6.6主震发生后至2014年12月6日MS5.8及MS5.9强余震前,在该阶段内,余震分布主要沿NW 向分布;第2个阶段是2014年12月6日MS5.8及MS5.9强余震发生后,余震的分布方向发生了明显变化,主要分布于NNW 向,且震源深度有加深的趋势。综合考虑余震的分布发展和每个时间段内地震个数的均衡,最后共划分为5个时间段,分别为景谷地震发生前(P1)、同震期(P2)、2014年12月6日MS5.8及MS5.9强余震前(P3)、2014年12月6日MS5.8及MS5.9强余震后(P4)及恢复期(P5),5个时间段的具体时间窗及地震数据如表1所示。这5个时间段内的地震个数相差不大,P波射线基本能覆盖景谷地震余震区及其附近区域(图3)。划分时间段后,用双差层析成像方法对每个时间段的数据进行3D速度结构及地震重定位的联合反演。首先进行地震对的匹配,由于地震分布差异,对时间段1和时间段5采用相同的匹配参数,即选择地震对之间的最大距离为30km,每个地震最多可与20个地震组成地震对。其余时间段则采用另一组匹配参数,即选择地震对之间的最大距离为20km,每个地震最多可与10个地震组成地震对。最终构建了相对到时数据(表1)。
图3 5个时间段的P波二维射线分布图Fig.3 Distribution of 2-D P-wave ray paths for the five periods.a时间段1;b时间段2;c时间段3;d时间段4;e时间段5
表1 5个时间段的数据Table 1 Data of five periods
我们选取MS6.6主震的震中(23.38°N,100.48°E)为坐标原点建立坐标系,逆时针旋转坐标轴,使Y轴与景谷地震序列的主破裂方向平行。综合考虑地震和台站分布后,对5个时间段均划分了相同的网格,并在反演前进行了不同的分辨率测试,以寻求最佳的网格分布。分别测试了3km×3km、5km×5km、10km×10km的网格间隔,最后发现地震分布密集的景谷地震震源区的横向分辨率为5km×5km、震源区外围地震较少的区域横向分辨率为10km×10km的网格分布的分辨率较为可靠,网格分布如图1b所示,垂直向网格位于0km、2km、5km、7km、10km、12km、15km、22km、31km、40km深度。准确的初始模型是获得稳定结果的关键,在建立初始一维模型时,以多个地壳模型(徐甫坤等,2015;陈飞,2017;李丹宁等,2017;Wanget al.,2018)为参考(表2)。由于云南地区的地质结构复杂,3D初始模型相对于1D模型的精度更高。为了获得更加准确的模型,在本研究中首先使用初始一维模型(表2)对所有的地震数据进行双差层析成像联合反演,得到3D速度结构和地震重定位结果,然后采用该3D速度结构模型作为初始模型分别对5个时间段的地震数据进行双差层析成像反演,以此得到较为稳定的P波速度结构。
表2 一维P波速度模型Table 2 1D P wave velocity structure
对所有地震数据进行地震匹配时设置地震对之间的最大距离为30km,每个地震最多可与20个地震组成地震对以构建相对到时数据。反演时,通过对数据方差与模型方差及解的方差进行权衡分析,最后确定阻尼因子和平滑参数的最优值分别为170和40。经过20次迭代后到时残差均方根由0.32降低到0.049,得到了较为稳定的结果。图4a为利用双差层析成像方法对所有地震数据进行联合反演后得到的5km、10km、15km深度处的P波速度结构,并利用棋盘测试方法对分辨率进行了测试(图4b),测试结果显示研究区中心即景谷地震震源区及其附近区域的分辨率较好,数据均可得到很好的恢复。由于分辨能力分布存在差异,故主要针对景谷地震震源区及其附近区域开展讨论。由图4a可看出景谷地震余震序列分布于P波高速异常区及低速异常区的交界,以南谷断裂(F3)和威远江断裂(F5)为界,断裂以北为高速异常区,以南为低速异常区。随着深度的增加,高速异常区范围缩小,低速异常区增大,这与李永华等(2014)采用接收函数和面波联合反演方法所得到的结果一致,即景谷地震及其余震位于上地壳高、低速异常交会地带。也与刘瑞丰等(1993)提出的“云南地区大部分MS≥5地震都发生在速度梯度较大的区域,特别是高、低速过渡地区”这一结论相符。徐甫坤等(2015)、李丹宁等(2017)及Wang等(2018)通过对景谷地震余震序列的重定位研究均发现余震深度分布具有北西浅、南东深的特点,且余震数量衰减具有北西快、南东慢的特点。这可能是由于余震序列的北西端处于高速异常区,坚硬脆性的岩石不利于余震的发展,而南东端处于低速异常区,可能存在破碎程度较高或者富含流体的岩石,这些岩石有利于余震的发展所致。
从图4a也可看出与澜沧江断裂相交或距离很近的南谷断裂(F3)和益香-赵家村断裂(F6)处于低速异常区,而与澜沧江断裂(F1)未相交的断裂(如永平盆地东缘断裂(F4)和景谷-云仙断裂(F7))在上地壳内主要处于高速异常区或高、低速异常的交会地带。这可能是由于与澜沧江断裂(F1)相交形成通道,为这些断裂附近区域的上地壳提供了大量流体,从而导致P波速度的低速异常,而澜沧江的流体对走向近SN的2条断裂没有影响。
在反演准备完成后,将2.1节反演出的整个研究区的三维P波速度结构作为初始模型,以重定位结果为初始地震数据分别对5个时间段进行双差层析成像反演,得到每个时间段的三维P波速度结构。在反演中,同样也通过权衡分析确定5个时间段阻尼因子和平滑参数的最优值。对5个时间段均经过16次迭代,最后的结果均可较好地收敛,反演后5个时间段到时残差的均方根显著降低。如表3所示,5个时间段最后得到的3D速度结构相对于初始结构在稳定性和精度上均有很大提升。
表3 5个时间段反演前后到时差的均方根残差变化Table 3 The RMS residuals between observed and predicted differential travel times based on initial 3D model and final 3D model for the five periods
此外,还使用棋盘测试评估模型的分辨率,分辨率结果如图5所示。P2—P4由于景谷地震的余震序列分布较为集中,故震源区的分辨率较好,但景谷地震发生前(P1)的地震分布较分散、震源区地震较少,导致震源区的分辨率不太理想,大部分区域在5km和15km深度处没有分辨率。恢复期(P5)包含景谷地震后期余震序列及其他分散的地震,故在震源区的分辨率也不太理想,但也能达到部分恢复。图6为5个时间段在5km、10km、15km深度处的P波速度结构,为了进一步量化模型分辨率,引入可与棋盘测试结果互补的DWS值,DWS值表示给定网格节点附近所有射线长度(乘以数据权重)的总和。DWS=100(图6中的白色等值线)包围了棋盘测试恢复较好的区域,所包围的区域能表现出较为可靠的P波速度分布。
图5 5个时间段不同深度处的P波速度棋盘测试结果Fig.5 Checkerboard resolution test of P-wave velocity for the five periods.a—e依次对应时间段P1—P5
图6 5个时间段不同深度处的P波速度结构和地震分布Fig.6 Distribution of P-wave velocity structure and earthquakes at different depths for the five periods.a—e依次对应时间段P1—P5。黑色圆点为地震震中,白色五角星为M S 6.6主震震中,白色圆点为M S 5.8和M S 5.9余震震中,白色实线表示DWS>100的区域。F1 澜沧江断裂;F2 窝拖寨断裂;F3 南谷断裂;F4 永平盆地东缘断裂;F5 威远江 断裂;F6 益香-赵家沟断裂;F7 景谷-云仙断裂;F8 无量山断裂。断裂构造数据引自毛泽斌等(2019)
在使用双差层析成像方法得到5个时间段的三维P波速度结构后,根据时移层析成像方法的技术路线,将2个相邻时间段中前一个时间段的速度结构作为初始模型,重定位结果作为初始定位数据,基于这2个时间段的地震数据所构建的相对到时数据利用时移层析成像方法反演得到P波的速度变化。在构建相对到时数据时,要求2个地震事件分别来自2个相邻的时间段,为了构建尽可能多但又不影响结果稳定性的相对到时数据,要求地震对之间的最大距离为50km,一个时间段内的每个事件最多可以与30个事件相连接,最后得到的相对到时数据如表4所示。景谷地震发生前(P1)—同震期(P2)虽然所使用的地震数目不少,但2个时间段定位所使用的台站差异较大,景谷地震发生前(P1)的地震分布分散,而同震期(P2)的分布更加集中,故最后挑选出的符合要求的相对到时数据较少,这会影响反演结果的精度。为了确保最后得到稳定的结果,所做的4次反演均经过32次迭代,反演后得到的观测到时差与理论到时差之间的均方根残差变化如表5所示。4次反演后残差的下降幅度较大,说明最后得到的P波速度变化与2个时间段内的地震事件较合适,P1—P2受数据影响而残差下降幅度较小。
表4 不同时间段的数据Table 4 Data of different periods
表5 不同时间段反演前后到时差的均方根残差变化Table 5 The RMS residuals between observed and predicted differential travel times based on initial 3D model and final 3D model for the different periods
同双差层析成像方法相同,我们也使用棋盘测试方法解释可接受分辨率的范围,为P波速度变化结果的可靠性分析提供参考。图7为棋盘测试分辨率结果分布图,P1—P2受数据数量影响,棋盘分辨率的恢复不太理想,其余的3个不同时间段的分辨率在震源区恢复得相对较好。
图7 不同时间段之间在不同深度处的棋盘测试结果Fig.7 Recovered checkerboard models for different periods.a时间段2相对于时间段1;b时间段3相对于时间段2;c时间段4相对于时间段3;d时间段5相对于时间段4
为了表现出分辨率较好的区域,采用与棋盘测试结果互补的DWS值,图8中DWS=500的红色实线所勾画的区域基本是分辨率相对较好的区域,我们在下文中也将重点讨论该区域的P波速度变化。2014年景谷MS6.6主震的同震P波速度变化是基于景谷地震发生前(P1)和同震期(P2)2个时间段的数据得到的(图8a)。在浅层深度处主震周边区域的P波速度下降,与景谷地震发生前该区域的P波速度相比下降幅度≤0.2% ;而在深度10km处,主震周边区域的P波速度上升,且上升幅度较大,与地震发生前相比上升了1.7%~2.6%;更深处的P波速度轻微下降。12月6日MS5.8及MS5.9 2个强余震发生前(P3)相较于主震同震期(P2)的P波速度变化是基于这2个时间段的数据得到的(图8b)。由图中可看出,主震周边区域的P波速度变化趋势与同震期近似,但余震序列北西端的P波速度在浅层出现上升现象,幅度为0.8%;在10km深度处下降,幅度≤1.5% ;在更深处出现轻微上升。12月6日的MS5.8及MS5.9 2个强余震发生前后的P波速度变化如图8c所示。主震周边区域在5km深度处的P波速度出现下降,下降幅度较之前的阶段更大,≤2.0% ;在10km深度处P波速度仍出现上升,上升幅度相较于之前的阶段更小;更深处的P波速度出现下降,且下降幅度≤2.2%。由之前的研究结果可知,在2个强余震发生后余震分布方向发生了明显变化,主要分布于NNW 向,且震源深度有增加的趋势,主要分布于9~16km处(李丹宁等,2017)。由图8c可看出,NNW 向余震分布区域的P波速度在余震分布较少的浅层区变化较小,而在余震分布较多的10km深度处P波速度则有所下降,下降幅度达3.8%,在15km深度处P波速度基本没有变化。2个强余震发生后约2a内,与2个强余震发生后不到一个月相比(图8d),震源区开始愈合,整个震源区浅层的P波速度出现上升,且上升幅度较大,远大于之前P波速度下降的幅度,在10km深度处的2个强余震震源区附近的P波速度同样上升幅度较大,而北西端余震区域的P波速度则下降。
图8 不同时间段之间在不同深度处的P波速度变化和地震分布Fig.8 Distribution of temporal P-wave velocity changes and earthquakes at different depths for different periods.a时间段2相对于时间段1;b时间段3相对于时间段2;c时间段4相对于时间段3;d时间段5相对于时间段4。黑色圆点表示地震震中,白色五角星为M S 6.6主震震中,白色圆点为M S 5.8和M S 5.9余震震中,红色实线表示DWS>500的区域,DWS<500的区域被阴影化。F1 澜沧江断裂;F2 窝拖寨断裂;F3 南谷断裂;F4 永平盆地东缘断裂;F5 威远江断裂; F6 益香-赵家沟断裂;F7 景谷-云仙断裂;F8 无量山断裂。断裂构造数据引自毛泽斌等(2019)
总体而言,震源区的P波速度变化与余震分布变化相关。2014年10月7日景谷MS6.6主震发生后至12月6日MS5.8及MS5.9 2个强余震发生前,余震主要分布在NW 向,首先在主震附近出现,然后向NW 端发展,主震周边区域浅层处的P波速度下降,下降幅度≤0.2%。而10km深度处的P波速度上升,且上升幅度较大,约为1.7%~2.6%。更深处(15km)的P波速度则下降,下降幅度与5km深度处近似。
2014年12月6日MS5.8及MS5.9 2个强余震发生后,余震的分布方向发生了明显变化,主要分布于NNW 向,朝SEE向发展,震源深度有加深趋势,主要在9~16km深度范围。在2个强余震发生后不到一个月,主震周边区域在5km和15km深度处的P波速度下降,下降幅度与之前的阶段相比更大,≤2.2% ,主震周边区域在10km深度处的P波速度仍上升,但上升的幅度更小。可以发现,主震周边区域在5km和15km深度处的P波速度变化一致,而在10km深度处却恰恰相反,推测10km深度处的变化可能与主震及余震的发生联系较小,可能还受其他因素的影响。2个强余震发生后,NNW 向余震聚集区在余震分布较少的5km深度处的P波速度出现轻微变化,在余震分布较多的10km深度处,该区域的P波速度下降明显,幅度≤3.8% ,而该区域在15km深度处没有变化,可以看出2个强余震的发生对10km深度处的P波速度影响较大。
MS5.8及MS5.9 2个强余震发生后约2a,主震周边区域开始愈合,各深度的结果均显示P波速度有所上升,但在5km和10km深度处上升的幅度远大于之前下降的幅度,故推测除受到愈合过程的影响外,应该还包含其他因素的影响。NNW 向余震分布区域从浅层到深层也开始愈合。
在余震分布的2个主要方向(NW 和NNW 向)(图1b)做了垂直剖面,图9为这2个垂直剖面的P波速度变化分布图,并将2个时间段中的后一个时间段的地震投影到剖面上,以观察速度变化与地震分布的变化关系。由图9a可以看出,景谷MS6.6主震发生后,余震的投影主要分布在AA′剖面上,AA′剖面主震周边的区域在浅层的P波速度出现轻微下降,而BB′剖面在浅层深度没有明显变化。同震期(P2)至2014年12月6日MS5.8及MS5.9余震发生前(P3),图9b中的AA′剖面的浅层P波速度下降,且从主震周边扩散至A端,A端是主震发生后余震较早到达的区域。在BB′剖面上,余震继续向B′端发展,与主震周边区域相连的B端的P波速度出现较为明显的下降。同时发现,2条剖面在5~15km深度处均存在着宽度不等的P波速度上升的条带区域,这个条带区域使得整个深度的变化不连续,条带下方深度的P波速度变化与浅层基本一致。在2014年12月6日MS5.8及MS5.9余震发生后(P4)(图9c),余震的分布方向发生了改变,主要集中在A′端和NNW 向的BB′剖面上,靠近主震的A端和B端的余震减少,2个5级以上强余震周边的余震增多,且余震分布深度相比之前的阶段有变深的趋势。AA′剖面的浅层深度的P波速度下降范围扩大至A′端,深度延续到10km,并逐渐将原来的P波速度上升条带覆盖,但P波速度下降的幅度不大。而余震分布增多的BB′剖面在0~9km深度范围内的P波速度下降至最大幅度,原来的P波速度上升条带也只剩下2km的宽度,说明MS5.8及MS5.9 2个强余震的发生对震源区影响较大。随着余震活动的减弱,震源区开始愈合。由图9d可以看出AA′和BB′2个剖面的P波速度上升,且上升幅度远大于之前的下降幅度。
图9 不同时间段之间的P波速度变化沿垂直剖面AA′和BB′的分布图(剖面位置见图1b)Fig.9 Distribution of temporal P-wave velocity changes along the cross sections AA′,BB′for different periods.a P2相对于P1;b P3相对于P2;c P4相对于P3;d P5相对于P4。红色实线表示DWS>500的区域,DWS<500的区域被阴影化,绿色五角星为M S 6.6主震震源,绿色圆形为M S 5.8和M S 5.9余震震源,黑色圆点表示地震
与大地震相关的速度变化已被广泛观测到(Poupinetetal.,1984;Pengetal.,2006;赵盼盼等,2012;Obermannetal.,2013;Koulakovetal.,2016;Peietal.,2019;温扬茂等,2019;王俊等,2020),地震发生前后出现波速变化的原因有很多,需要结合地质和地球物理结果进行探讨。在大地震发生后会产生地震波速变化的深度尚不清楚,一般来说,速度变化能反映由强地面运动所造成的岩石损伤,预计主要在顶部浅层区域出现,破裂区地下介质内的结构变化,开放裂缝及流体变化导致应力场的变化可能发生在不同深度(Peietal.,2019)。一般而言,7级以上地震往往能在地表形成明显的地震地表破裂带,但也有少量约6.5级的地震也可形成地表破裂(邓起东等,1992)。景谷地震发生后,云南省地震局组织了野外应急科学考察,常祖峰等(2016)根据考察结果发现此次景谷地震震中区域产生了广泛的砂土液化现象和地裂缝,并提出虽然这些地裂缝单条规模不大,但却集中出现在极震区内。这与处于极震区的景谷主震震中附近的浅层深度的P波速度下降相符,表明该区域的P波速度下降可能是由于岩石破坏所致。景谷主震附近的P波速度在主震发生后并没有到达最大波速降,而是在震后2个月后出现最大波速降,这与庞卫东等(2017)及刘志坤等(2010)提出的大地震发生后最大波速降出现在震后1~4个月内的观点相符。且P波速度下降的幅度与已有的相同量级地震发生前后的速度变化研究结果一致,如Li等(1998)计算得到2004年Park fieldM6地震发生后San An dereas断裂附近台站下方的地震波速下降了2.5%;Nishimura等(2000)计算得到1998年Intra PlateM6.1地震发生后震源区的地震波速下降了0.3%~10%;Schaff等(2004)计算得到1984年Morgan HillMS6.2地震同震期的P波速度下降约2%。在P波速度的下降幅度到达最大后景谷主震震源区开始愈合,说明主震震中附近区域主要受主震的影响。
我们发现2个5级强余震发生前在5~15km深度范围内存在一个宽度不等的P波速度上升条带,该深度范围的P波速度不受主震破裂及余震的影响。该P波速度上升条带与其他深度不一致,导致整个震源区的P波速度变化在深度上不连续。2014年12月6日MS5.8及MS5.9 2个强余震发生后,在余震分布较多的区域内,该P波速度上升条带呈现出P波速度下降的现象,且下降幅度较大,甚至超过了主震附近区域P波速度下降的最大幅度,但也符合Nishimura等(2000)的结果,可能是2个5级以上强余震及其后余震的发生对该深度有影响。王烁帆等(2019)通过对景谷地震震源深度的研究得到了岩石强度模型,其认为主震的起始深度为9.5km,矩心深度为5.0km,主震破裂于深部然后向浅部发展,主震破裂区主要位于浅层,而2个5级强余震的起始深度和矩心深度均约为10km,呈圆盘式破裂,故对深10km处上、下的区域均有影响,推测脆韧性转换带位于10km深度处,主要分布于5~15km深度内。他们还与地壳电性结构结果(程远志等,2016)相结合证实了5~15km深度处的介质为高阻,岩石强度较高,与脆韧性转换带的存在一致,如图10所示。故在此推测由于脆韧性转换带的存在导致了P波速度变化在深度上的不连续性。
图10 景谷地震序列分布图与震源区电性结构及岩石强度模型(改自王烁帆等,2019)Fig.10 The hypocentermap of the Jinggu earthquake sequence and the electrical structure and rock strength model of the source area(after WANG Li fang et al.,2019).a景谷地震余震分布图(Wang et al.,2018);b、c余震在2个剖面上的分布(Wang et al.,2018);d景谷地震震源区电性结构模型(改自程远志等,2016),其中白色实线表示岩石强度随深度的变化曲线(数据源自孙玉军等,2013)
15km以下深度的P波速度变化与浅层深度的P波速度变化基本一致,但在该深度范围内不可能受到主震破裂的影响。由王烁帆等(2019)的研究结果可知,该深度范围属于低电阻率的韧性剪切区,岩石强度较低,可能存在部分熔融和脱水现象,因此可能在地震发生后微破裂及体积应变降低的同时伴随着岩石裂纹的增长,而P波速度在干燥裂隙的环境下随裂隙密度的增加而减小,故较深处的P波速度有所下降。在2个5级强余震发生后,余震震源深度的分布范围变深,促使较深区域的岩石裂纹增加,从而导致P波速度下降。
震源区的P波速度受2个5级强余震影响而降低后,整个震源区开始愈合,P波速度上升,但上升的幅度大于之前下降的幅度,推测该区域除景谷地震的愈合过程外,可能还有应力积累的过程。景谷地震发生后的2018年9月8日,在云南墨江发生了MS5.9地震,中国地震台网中心测定的该地震震中(23.38°N,101.53°E)与景谷地震主震震中相距约117km,均属于滇中区域。云南省地震局形变测量中心2017年年度会商报告中提出,大地震发生前,在震中周边应当存在应变场大幅调整的信号,部分区域应变场变化显著可能是强震发生前的重要特征,且存压性活动增强的区域范围较大,而非小区域的变化(1)云南省地震局形变测量中心,2017,2017年年度会商报告。。在该报告中,利用应变场时序分析法(洪敏等,2014)所计算出的云南区域的应变场变化显示,2018年9月8日墨江MS5.9地震发生前滇西南和滇中的挤压应力增强表明滇中地区的应力场正在增强。据此推测,景谷地震震源区的P波速度大幅上升是震源区愈合及应力积累共同作用的结果。
本文基于云南省地震区域台网记录的2008—2017年的地震观测报告数据,使用双差层析成像方法得到了2014年云南景谷MS6.6地震震源区的三维P波速度结构,再使用基于双差层析成像的时移层析成像方法得到了2014年云南景谷MS6.6地震发生前后景谷地震震源区的P波速度时空变化,并对P波速度的变化特征及其机制进行了深入细致的研究,得到以下几点认识:
(1)景谷地震的余震序列分布于P波高速异常区及低速异常区的交界处,以南谷断裂和威远江断裂为界,断裂以北为高速异常区,以南为低速异常区。随着深度的增加,高速异常区域范围缩小,低速异常区域增大,余震分布特征与P波速度分布特征相关。此外,与澜沧江断裂有所相交或距离很近的南谷断裂和益香-赵家村断裂处于低速异常区,可能是受到断层中流体的影响。
(2)与野外应急考察结果相结合发现处于极震区的景谷主震震中附近的浅层P波速度下降是由于岩石破坏导致的,且P波速度下降的最大降幅出现在景谷主震发生2个月后。
(3)在不受主震破裂影响的5~15km深度处存在一个P波速度上升条带,该条带的存在导致震源区的P波速度变化在深度上不连续,推测该条带存在的原因是该深度范围主要为高强度及高阻介质的脆韧性转换带。在2014年12月6日MS5.8及MS5.9 2个强余震发生后,余震的分布方向发生了改变,主要分布于NNW 向,且震源深度分布有加深的趋势。2个5级强余震及其余震主要分布于脆韧性转换带内,从而影响了5~15km深度范围内的介质,导致其P波速度下降,且下降幅度较大,表明2个5级强余震的发生对脆韧性转换带造成了影响,且余震的发生特征通常与P波速度的变化特征一致。
(4)在震后愈合阶段,震源区P波速度上升的幅度远大于之前下降的幅度,P波速度超过震前水平,表明该区域不仅存在震后愈合过程,还包含其他物理过程。与用GPS观测数据所计算的应变场变化结果相结合推测,2018年9月8日云南墨江MS5.9地震发生前的应力积累过程与震后愈合过程相叠加导致P波速度出现大幅度上升。
与大地震发生有关的P波速度变化的原因有很多,对于深度相对较浅的区域通常可以结合地球物理观测和地质调查结果说明原因,但较深区域的未知信息较多,对速度变化机制的探讨可能存在一定问题。在今后的研究中,将结合其他方法所得到的更加精确的速度变化结果,将15km深度以下的速度变化结果与更加广泛的地球物理观测结果(如岩石孔隙水、应力的变化、地温信息和地球动力学)相结合来进行更加深入的研究。
致谢云南省地震局毛泽斌助理研究员为本文提供了断层数据;云南省地震局王光明助理研究员在画图方面提供了帮助;中国科学院测量与地球物理研究所王烁帆博士为本文提供了图10的图件;审稿专家为本文提供了宝贵的意见和建议。在此一并表示感谢!