2016年阿克陶MW6.6地震的同震断层滑动分布及库伦应力变化

解 柳,吴绍宇,康明哲,曾庆伟,朱毅诚

(1.西安科技大学,陕西 西安 710054;2.陕西铁路工程职业技术学院,陕西 渭南 714000;3.西安市政设计研究院有限公司,陕西 西安 710068)

0 引言

北京时间2016年11月25日22时24分,在新疆维吾尔自治区阿克陶县境内发生了MW6.6地震。此次地震震中位于阿克陶县木吉乡和库斯拉甫乡之间的高寒山区,发震断层为NWW走向的右旋走滑断裂,破裂长度超过77 km[1]。不同机构和学者对此次地震断层运动进行了研究,美国地质调查局(USGS)发布的破裂起始点坐标为(39.273°N,73.978°E),震源深度为17 km,认为它是具有拉张分量的右旋走滑事件。中国地震台网(http://www.cea.gov.cn/publish/dizhenj/464/515/20161125231007712890916/index.Html)测定的地震震中位置为(39.27°N,74.04°E),震源深度10 km,为右旋走滑兼正断层性质(表1)。张勇(http://www.ceaigp.ac.cn/tpxw/275080.html)研究的震源破裂结果表明该地震为单侧破裂,从震中开始向东传播。由于该地区地质、地形也相对复杂,有学者认为该地震为单一破裂。利用雷达数据提取LOS向形变图,参考各机构的震源机制解的参数,进行同震的滑动分布反演[2]。也有学者认为该地震为双断层破裂[3],并利用升降轨数据获取该地震的LOS向形变,以中国地震台网确定的震源机制解的参数进行滑动分布反演,发现其破裂分别位于震中以东约7 km、约33 km处,主要的滑动量集中在地下深度0～20 km处,最大滑动量为0.84m,并基于弹性半空间形变模型采用两段非均匀断层滑动模型进行精细滑动分布反演,结果显示断层西段的第一次事件以右旋走滑为主,走向103°,倾角76°,断层东段的第二次事件为右旋走滑兼具正断作用,走向109°,倾角略缓约55°[4]。

表1 不同机构给出的阿克陶地震震源机制解Table 1 The focal mechanism solutions ofAketao earthquake given by different organizations

本文在采用贝叶斯算法获取断层几何参数的基础上,反演了单一断层的滑动分布;采用粒子群算法获取双断层的最优倾角值,并进行精细滑动分布反演,并以滑动分布反演结果对该地震的库伦应力变化进行了分析。

1 研究区域及数据

1.1 构造背景

2016年发生的阿克陶MW6.6地震位于印度板块沿北东向与欧亚板块碰撞的西触角地区[5],该区域处于南天山西段和帕米尔高原交会地区、木吉断陷盆地西端地区附近(图1)。不同机构给出的震源机制解如下表1。

从2000年以来距震中100 km范围内共发生了2次6级地震、13次5级地震[6]。自阿克陶主震发生后,截止到2016年11月27日24时,记录到5.0～5.9级余震 1次、4.0～4.9级余震5次,最大余震震级为 5.0级[1]。

注:黑白色沙滩球表示不同机构发布的震源机制结果,红色与绿色圆点为余震分布图1 阿克陶地震震中区域及邻区地形、断层、余震Fig.1 Topography,faults,and aftershocks in the epicenter area and adjacent areas of the Aketao earthquake

1.2 数据说明

1.2.1 数据介绍

本文采用欧空局(European Space Agency,ESA)哥白尼计划中的地球观测卫星分别于2014年4月和2016年6月发布的卫星数据Sentinel-lA升轨单视复数影像和Sentinel-1B降轨单视复数影像。其具体信息见表2。

表2 InSAR数据信息Table 2 InSAR data information

1.2.2 同震形变解算

利用SLC影像数据生成强度图并结合DEM影像进行地理编码、差分干涉处理,生成差分干涉影像。为提高干涉图的质量,对干涉图进行自适应滤波处理、相位解缠、地理编码,最终提取LOS向形变信息(图2)。在图2中,(a)为降轨的LOS向形变场分布,地震断层南侧是木吉盆地,其干涉条纹明显,其西部形变值约为14 cm,东部形变值约为23 cm。而北部山区地形起伏,影像失相干,故干涉条纹模糊,但不影响获取形变范围。(b)为升轨的LOS向形变场分布,东西两部分形变量值相当,约为12 cm。相对于西部条纹分布而言,东部条纹分布范围比较小,条纹紧密。

图2 D-InSAR LOS向形变场Fig.2 Deformation field in LOS direction from D-InSAR

1.2.3 剖面分析

为更清晰地了解升降轨LOS向形变场信息,以两个子破裂为中心取同一经度的剖面,对其数据进行四次多项式拟合(图3)。浅蓝色为A1-A2剖面的形变值,最大沉降值约为25 cm,浅蓝色为A1-A2剖面的形变值,最大抬升值约16 cm,东部子破裂区的沉降量明显大于西部[图3(a)]。由剖面C1-C2,D1-D2可看出升轨的形变主要以抬升为主,最大抬升值约为13 cm,而最大沉降值仅约6 cm[图3(b)]。对比图3,由于升降轨的观测模式不同,造成升降轨LOS向形变场的沉降、抬升趋势相反,升轨的形变量值相对降轨的要小一些,但剖面图的形变趋势与整体形变趋势一致。

图3(a)表示降轨的剖面图,浅蓝色表示以东南部破裂为中心,经度为74.399°,纬度为39.05°至39.39°之间的取值(降轨解缠图白线A1-A2),黑色为其拟合曲线;深蓝色表示以西南部破裂为中心,经度为74.197°,纬度为39.05°至39.39°之间的取值(降轨解缠图白线B1-B2),红色为其拟合曲线;图3(b)表示升轨的剖面图,浅蓝色表示以东南部破裂为中心,经度为74.302°,纬度为39°至39.47°之间的取值(升轨解缠图白线C1-C2),黑色为其拟合曲线;深蓝色以西南部破裂为中心,经度为74°,纬度为39°至39.47°之间的取值(升轨解缠图白线D1-D2),红色为其拟合曲线图3 同震差分剖面图Fig.3 Co-seismic difference section

2 断层几何与滑动分布

本文依据升降轨提取的形变场估计断层的几何参数并反演断层的滑动分布。假设为均匀滑动模型,输入相应的断层几何参数。将断层面划分成若干子断层片,使用线性反演方法来估计断层面的滑动分布。

2.1 断层几何与滑动参数

本文结合矩形位错模型[7]采用贝叶斯算法(GBIS)将非线性反演应用于断层几何[8],对于反演,使用结合了Metropolis-Hastings算法的马尔可夫链蒙特卡罗方法(MCMC)估计最佳参数及其不确定性[9-10]。为了提高反演效率和计算速度,使反演结果得到更好的收敛,采用四叉树方法对D-InSAR数据进行二次采样[11],设置降轨的采样阈值为0.024 m,采样后数据点降为1 266个;升轨的采样阈值为0.014 m,采样后数据点降为1 382个。以降轨和升轨数据同时作为约束条件,设置断层走向的限制范围为-180°至180°,倾角的限制范围设置为0°至90°,迭代次数设置为106次,使得所有模型参数都根据其不确定性得到了很好的约束。反演结果显示深度最佳值为17.35 km,与USGS发出的17 km相近,断层的倾角最佳值为73.8°与USGS发布的76°相差不大,断层的走向最佳值为106.9°,与USGS发布的107°基本吻合。震级与USGS也基本吻合,为MW6.64(表3)。

表3 断层几何与滑动参数Table 3 Fault geometry and sliding parameters

图4(a)为计算断层参数时,模拟的降轨同震形变场,存在5个条纹周期,对应的形变量值为14 cm,与图2(b)中西部的形变量值一致,而与图2(a)中东部的形变量值最大差异为6 cm。图4(b)显示出东部的残差值明显大于西部,其最大残差值为～5.6 cm。图4(c)为计算断层参数时,模拟的升轨同震形变场,也存在5个条纹周期,对应的形变量值为14 cm,图4(d)显示出西南部条纹较密集,具有较大残留信号,残差值为～14 cm。根据表3断层参数的模拟LOS方向的形变场结果来看,降轨模式优于升轨模式。

图4为采用表3中的断层参数计算的降、升轨模式的形变场模拟与残差结果,长箭头表示卫星飞行方向,短箭头表示视线方向,每个条纹周期代表2.8 cm图4 同震干涉图的模拟与残差Fig.4 Simulation and residual of coseismic interferogram

2.2 单一断层模型的断层滑动分布反演

根据表3中的断层几何参数结果建立单一断层模型[12],起点坐标设置为(39.342°N,73.732°E),沿走向和下倾方向分别将断层扩展到70 km和20 km,将断层划分为2 km×2 km的子断层块,共计350个。SDM反演断层面上的滑动分布结果如图5所示,显示断层运动以右旋走滑为主,断层滑动分布主要集中在长度方向上30～70 km 之间,深度方向上5～20 km的区域。其中断层面上的两个滑动峰值对应的滑动量分别为 0.83 m、0.66 m,整个断层面的平均滑动量为0.31 m,对应的地震矩震级约为MW6.65。

图5 单一断层模型的滑动分布Fig.5 Slip distribution of the single fault model

图6(c)可以看出,单一滑动分布仍有较大的残差,其残差最大达到5.03 cm,说明单一断层模型并不能很好地解释观测数据,由此判断发震断层可能为两段断层共同破裂而成,所以将断层分为两段断层模型再次进行精细滑动分布反演。

黑色线为木吉断裂图6 单一断层模型均匀滑动分布模拟的形变场Fig.6 Deformation field simulated by uniform sliding distribution of the single fault model

2.3 双断层模型的断层滑动分布反演

以降轨LOS向形变结果为约束,参考单断层反演结果,将断层一分为二。断层面Ⅰ起点设置为(39.308°N,73.78°E),长度为40 km,宽度为20 km,断层面II起点设置为(39.191°N,74.332°E),长度设置为30 km,宽度为20 km,走向仍为106.9°,倾角采用粒子群算法[13],通过大量的反演实验,确定断层面Ⅰ最优倾角为70.79°,断层面Ⅱ最优倾角为55.33°,图7为双断层滑动量分布的反演结果,图7(a)表示断层面Ⅰ的结果,断层面滑动峰值的深度11.61 km,最大值为0.68 m,(b)表示断层面Ⅱ的结果,断层面滑动峰值的深度为7.46 km,最大值为0.77 m,对应的地震矩震级约为MW6.62。

图7 双断层模型的滑动分布Fig.7 Slip distribution of the double fault model

对比双断层滑动分布的结果计算的形变场[图8(b)]与dLOS观测值[图8(a)]发现两者整体分布形态具有一致性,数量级相同,形变量在断层东部存在一定差异,其最大值为3.09 cm[图8(c)],相比单一断层模型图8(c)而言,残差值减小了1.94 cm。

2.4 库仑应力

地震的同震位错会使周围断裂或者临近区域的应力状态发生改变,从而引起库伦应力的变化[14],并可能影响临近断层状态,从而改变地震发生的概率[3]。本文在单一、双断层模型的基础上,计算了震源深度17.35 km处的最优破裂面库伦应力变化(图9)。图9中的黑色点表示2016年11月26日0点至11月29日24点发生的余震,可以看出地震的余震基本发生在断层的南部,与应力的增强区相对吻合;绿色点表示震后自2017年1月1日到2020年10月28日发生的地震,图9(a)、(b)分别表示为单一断层和双断层的库伦应力结果。对比可以看出,双断层在北部和南部的正库伦应力变化增加了,且在双断层滑动分布的情况下计算的库伦应力变化更符合震后地震的分布。

蓝色区域表示库伦应力的减小区,红色区域表示库伦应力的增强区;深蓝色的线表示该区域附近的断层,其中F1为帕米尔主断裂带,F2为木吉断裂,F3为公格尔拉张系北段;黑色的点表示为余震的分布区域,绿色的点表示从2017年至今发生在该区域的大于3级的地震图9 单一断层、双断层最优破裂面的库伦应力Fig.9 Coulomb stress of optimal fracture surface for single fault and double fault

3 讨论

本次地震发生在新疆西南缘高寒地区,综合前人对帕米尔地质构造研究、InSAR数据解算结果,反演阿克陶地震的滑动分布,在认识余震分布的基础上进行库伦应力变化分析,结合上述结果,展开如下讨论:

双断层能否更好地反演断层滑动分布？根据单一断层反演得到的同震破裂长度为64.2 km,宽度为19 km,走向为105°,倾角为73.8°,与USGS第Ⅱ节面近似,但反演结果残差约5 cm。从反演的单一滑动分布结果来看,两个子破裂中间存在一个较小的滑动,可判断出该地震存在复杂震源破裂,可能在短时间内发生了两次破裂。继而进行双断层反演,采用粒子群算法反演出西段断层长度为48.4 km,宽度为19 km,走向为105.68°,最优倾角为70.79°,深度为5.94～22.94 km,东段破裂长度为14.99 km,宽度为19 km,走向为104.08°,最优倾角为55.36°,深度范围为5.82～15.69 km,最大拟合残差值约3 cm。对比拟合残差结果,可见双断层的模型能更好地模拟断层结构。

库伦应力变化的计算,重点是创建地震的震源模型。本文使用反演的有限元断层模型进行计算。输入350个的小断层块的滑动量去计算库伦应力的变化,后期将余震和后期发生的地震与之对应,可以看出,余震发生在正库伦应力变化区的概率均相对较高,从2017年至今发生的地震在基于双断层反演结果计算的正库伦应力变化区达到57%,可见,库伦应力变化为正值时,很有可能会引发后一次地震事件的发生。阿克陶地震前后有多次地震,可见该区域的活动性是相当活跃,推断该地区未来仍有发生强震的可能。

4 结论

本文通过多次对比实验,综合前人的研究结果,对滑动分布再一次进行了研究,给出了该地震的构造模型,得到了以下几点结论:

(1)采用贝叶斯方法反演了阿克陶地震的单一断层几何,得到其最优走向为106.9°,倾角为73.8°。反演滑动分布结果显示存在两个滑动峰值,其中西部最大滑动量为0.66 m ,东部最大滑动量为0.83 m,对应的矩震级为MW6.65。

(2)采用粒子群算法计算阿克陶地震同震断层西段(断层面Ⅰ)与东段(断层面Ⅱ)的倾角分别为70.79°、55.33°。反演结果显示断层面Ⅰ为走滑断层,滑动峰值的深度11.61 km,最大滑动量为0.68 m;断层面Ⅱ是以走滑为主兼正断性质,滑动峰值的深度为7.46 km,最大滑动量为0.77 m;总的矩震级为MW6.62,整体拟合残差比单一断层模型的结果降低了约2 cm。

(3)最优破裂面的库伦应力变化,其应力增强区与余震分布相对吻合。从震后至今该区域发生的3级以上地震,其发生在应力增强区的概率达到50%。

致谢:感谢中国地震局地球物理研究所房立华研究员提供的余震精定位数据;欧空局提供的Sentinel-1A与Sentinel-1B卫星数据;德国GFZ的汪荣江老师提供的SDM反演程序。部分图由GMT绘制。感谢两位评审专家对本文提出的宝贵意见！