2. 引力与固体潮国家野外观测研究站, 武汉市森林大道196号, 430000;
3. 中国地震局地震大地测量重点实验室, 武汉市洪山侧路40号, 430071;
4. 湖北省地震局, 武汉市洪山侧路48号, 430071;
5. 武汉市测绘研究院, 武汉市万松园路209号, 430022
据中国地震台网测定,北京时间2020-06-14 22:24土耳其(40.75°E, 39.40°N)发生MW5.9地震,震源深度约10 km。根据土耳其灾害和紧急情况管理局(Turkey’s Disaster and Emergency Management Authority, AFAD)发布的相关研究报告(https://deprem.afad.gov.tr/depremdokumanlari/1836),此次地震矩震级为MW5.7,持续时间约15.04 s。AFAD依据T50降轨数据计算出,主要的视线向形变量约7 cm。
如图 1所示,研究区构造较为复杂,主要受安那托利亚板块、阿拉伯板块、欧亚大陆板块三大板块的共同作用,汇聚于卡尔勒奥瓦三联点(KTJ),其结果使得安那托利亚板块向西不断运动,发育的主要断层为北安那托利亚断层与南安那托利亚断层[1]。其中,北安那托利亚断层是以右旋走滑为主的断层,且兼具少量的逆冲特性,是一条1 200 km长的走滑断层,从卡尔勒奥瓦向西延伸;南安那托利亚断层以左旋走滑为主[2-3]。GPS数据解算出北安那托利亚断层右旋走滑量为24 mm/a,南安那托利亚断层左旋走滑量为10 mm/a[4]。1976~2020年MW≥5地震目录(https://www.globalcmt.org/CMTsearch.html)统计结果显示,截至2020-07-01,该区域共发生49次5级以上地震,地震主要沿着这两条大断层分布。其中MW≥6.0地震共5次,最近的一次为2020-01-24发生在南安那托利亚断层附近的MW6.8地震。
图 1中蓝色矩形框(40.55°~41.05°E, 39.25°~39.55°N)内为选取的重点研究区域。本次地震发生在该范围内,位于北安那托利亚断层附近,该区域历史地震活动性较强,聚集了大量5~6级地震,研究本次地震的相关特征,有利于分析所在区域的地震危险性。
本次地震附近GPS站点稀疏且未观测到明显的同震形变。另外,地震波数据一般受台站分布的影响,反演得到的发震断层的走向和倾角精度存在较大误差,且存在中强地震震中定位不精确的现象。如图 1所示,GCMT、AFAD、USGS等3个机构提供的震中位置存在差异,这给后续的同震破裂研究带来诸多不便。近场的大地测量数据,比如InSAR数据,能够直观地获取近场形变数据,可以很好地约束发震断层的几何形状参数。
综上,本文利用Sentinel-1A数据,通过D-InSAR技术提取同震形变场,采取贝叶斯自举优化法研究发震断层几何形状特征[5],用有限断层方法反演断层破裂滑动分布[6],为准确认识此次土耳其MW5.9地震的形变特征以及发震构造提供参考。
1 InSAR数据数据处理和形变分析本文使用欧空局(https://sentinel.esa.int/)升降轨Sentinel-1A干涉宽幅数据,对于升轨T43数据,震前影像拍摄时间是2020-06-02,震后影像拍摄时间是2020-06-14(具体时间是15:18:35.017,距震后约1 h),经过处理后的干涉图没有明显的同震信号,是噪声过大所致,所以舍弃该结果。考虑到2020-06-15 06:51:31(UTC)该区域发生MW5.5余震,如果影像时间跨度包含余震事件,则在后续的同震形变提取中难以剔除,所以选取的影像时间跨度上应尽量不包含MW5.5余震信号。综上,本文采用降轨T50数据,震前影像拍摄时间是2020-06-03,震后影像的结束拍摄的准确时间为2020-06-15 03:17:50.848(UTC),影像相关的具体信息如表 1所示。
使用InSAR scientific computing environment (ISCE) [7]对降轨数据进行D-InSAR处理,采用欧空局发布的精确轨道文件来减小轨道误差的影响,采用30 m分辨率的SRTM数字高程模型消除地形的影响[8],利用Goldstein滤波方法提高信噪比[9],基于Snaphu(statistical-cost, network-flow phase-unwrapping algorithm)[10]算法进行相位解缠。然后,通过地理编码获取土耳其地震降轨的同震形变场(图 2)。利用generic atmospheric correction online service(GACOS)模型[11](http://ceg-research.ncl.ac.uk/v2/gacos/)消除对流层大气延迟所引起的误差,最终得到经大气改正后的形变场(图 3)。
从图 3(a)~(c)可以看出,通过GACOS模型改正,削弱了部分对流层的影响,大气校正的数值范围为-0.95~0.8 cm,同震形变场得到改善,改正后的视线向最大沉降量为-7.75 cm,最大抬升量为8.87 cm。结合干涉图可以发现,此次土耳其MW5.9地震的地表形变影响范围约15 km×10 km,降轨形变场长轴方向近EW向,LOS向形变场呈椭圆状(关于东西向基本对称)。
图 4为降轨形变场横截面示意图,图 4(a)中AB、CD表示降轨LOS向形变场中横截面的位置,图 4(b)为对应横截面的LOS向同震位移变化。横截面AB(黑线表示)靠近CENC震中,从A点至B点,LOS向形变值由约-0.08 m变化至0.08 m,且呈现中心对称,显示出走滑地震特性;横截面CD变化相似,范围在-0.06~0.06 m。LOS向形变量绝对值大于0.02 m的区域范围约10 km2,表明此次地震影响范围较小。此外,图中还显示出不同机构确定的震中,黑色五角星表示CENC确定的震中,白色五角星代表本文利用InSAR数据确定的震中,位于LOS向形变场的形变中心附近。
为加快反演进程,采用四叉树方法[12]对降轨数据进行降采样,考虑到近场形变的重要性,在降采样的过程中,对近场形变进行加密采样,对远场区域进行稀疏采样,最终保留150个降轨数据点。
2.1 断层几何形状反演利用贝叶斯自举优化法(Bayesian bootstrap optimization, BABO)估算发震断层相关的几何形状参数[5]。BABO基于直接搜索,从定义的模型空间中抽取随机模型参数,然后计算综合模型值,并将其与目标的观测数据进行比较。BABO优化的核心是对观测数据和预测数据之间的差异进行逐点计算,求取观测值与模型值最小化差异:
$ \begin{array}{l} {\rm{min}}{\left\| {{d_{{\rm{obs}}}} - {d_{{\rm{synth}}}}} \right\|_2} = \\ {(\sum | {d_{{\rm{obs}},i}} - {d_{{\rm{synth}},i}}{|^2})^{\frac{1}{2}}} \end{array} $ | (1) |
式中,dobs为观测值,dsynth为模拟值,i(i=1, 2, 3, …)为观测点个数。
BABO将断层面作为一个均匀面来处理,优化得到断层面参数(断层长、宽、倾角、走向、滑动角等)。本文的观测数据是经过四叉树降采样处理后的T50降轨Sentinel1-A数据,所有优化参数序列如图 5所示,可用于检查收敛性以及查看模型参数是否推动了给定的边界范围,较高的失配值为冷蓝色,较低的失配值为红色,失配值越低表明参数优化越好。对所有估计样本值进行排序,置信区间的下限取排序后的5个百分位数,置信区间的上限取排序后的95个百分位数,取两者之间的数值,得到各参数的90%置信区间。各参数迭代次数均为50 000,最终收敛状况都较好,且失配值随着迭代的不断进行而不断地降低。
图 6表示模型参数的分布,x轴绘图范围表示给定的初始参数边界,y轴表示概率密度函数(probability density function,PDF),本文取高斯核密度函数。图中显示模型参数分布十分集中,误差较小。最终各参数的具体数值见表 2,以下仅列出最优解:东偏移、北偏移表示UTM坐标系下相对于参考点(40.68°E,39.35°N)(GCMT震中)的偏移量,物理含义上表示反演得到的理论震中位置,转化成经纬度即为(40.754°E,39.389°N)。此外,断层上边界中心点深度为1 308 m,断层长度为5 680 m,断层宽度为5 710 m,滑动量为0.223 m,断层走向、倾角、滑动角分别为257.48°、79.69°、154.2°。
在反演断层几何参数时,得到的残差结果如图 7所示。图中圆点表示采样点,矩形表示反演断层所在位置,矩形中的黑点表示震中,断层左上角点经纬度为(40.780 4° E,39.390 8° N),对应深度为1.308 km。在靠近震中附近,LOS向视向线模型值都能够较好地拟合观测值,降轨数据残差范围是-0.018~0.012 m,残差的RMS为6 mm,总体来看模拟效果较好。
基于上述BABO反演获得的几何形状参数,采用有限断层方法对发震断层滑动分布进行反演。为了获得较精细的地表变形,本文的震源矩形断层采用三角格网[13]构建,断层位错引起的地表变形通过弹性位错公式计算,其本质上属于非负最小二乘方法,即寻求观测值拟合度和滑动分布粗糙度最小化:
$ {\left\| {\mathit{\boldsymbol{W}}(GS{\rm{ }} - d)} \right\|^2} + {\beta ^2}{\left\| {LS} \right\|^2} = {\rm{min}} $ | (2) |
式中,d取InSAR视向线数值;W为观测值的权矩阵;G为格林函数,可以通过基于半无限空间地壳分层的弹性位错模型[14]计算获得;β为子断裂滑动平滑因子;L为拉普拉斯二阶差分算子;S为待求参数,包含子断层滑动量以及InSAR轨道相关的改正数[15]。
断层参数设置:通过观测InSAR干涉图和LOS向形变场的空间特征,参考贝叶斯自举优化法得到的结果(表 2),走向取257.5°,倾角取79.7°,延伸反演断层的长宽,设置断层长25.5 km,宽9.0 km。坐标原点取断层左上点(40.780 4°E,39.390 8°N),对应深度1.308 km。反演采用降轨相关数据。
最终获得的土耳其地震同震滑动分布如图 8所示。该地震主要破裂区域长度约8 km,宽度约6 km,破裂没有出露地表,位于地表以下约0.8 km,最大埋深约8.9 km。主要破裂集中在深度3~6 km的区域,较深区域的子断层滑动角几乎垂直,滑动角随着深度的增大而增大,在接近地表的地方子断层几乎呈水平右旋滑动,主要滑动区域的平均滑动角为146.8°。主要破裂区左侧存在少量的滑动现象,可能是由于本文采用单一断层模型所致。同震破裂主要表现为右旋走滑,且兼具少量逆冲特性。断层面上最大滑动量约0.57 m,位于(40.769°E,39.394°N)处,埋深约4.278 km。假设地壳的剪切模量为3×105 bar,此次地震释放的地震矩约4.54×1017 Nm,对应矩震级MW5.7。滑动分布模型的形变残差如图 9所示,模拟效果较好,降轨拟合残差的范围是为-3.9~4.0 cm,残差的RMS为0.7 cm。
根据表 3,InSAR反演的震中与CENC测定的震中最接近,相差1.3 km。结合图 4分析,InSAR得到的震中更符合实际震中特征,更为精确,且更接近最大形变值所处的位置,而与USGS、AFAD测定的震中在平面空间上分别相差约4.1 km、4.4 km,与GCMT测定的震中距离最远,相差约7.7 km,这可能主要受测震台网的分布、测震站台的位置以及地震震级(较小)等影响。近场数据测定的震中明显要比USGS、GCMT更加准确,这表明用近场InSAR形变场数据可以更好地确定中强地震的震中位置。
Poyraz[16]研究认为,NAF最东端的平均滑移率为21 mm/a,闭锁深度为12.72 km。此次MW5.9地震发生在NAF的东段,造成震中局部地区近东西向的右旋剪切作用。为评估此次地震对土耳其地区带来的影响,假设摩擦系数为0.4,接收断层的走向为257.48°,倾角为79.69°,滑动角为154.2°(即上述反演结果,见表 2),得到最大滑动量所处深度4.278 km对应的库伦破裂应力(CFS)变化[17],ΔCFS>0表示库伦破裂应力增加,ΔCFS < 0表示库伦破裂应力减少,如图 10所示。值得注意的是,当ΔCFS≥0.1 bar时,就有诱发后续地震的可能性[18]。由图可见,南安那托利亚断层所在区域库伦破裂应力略有增加,增量小于0.1 bar,表明此次地震对其影响较小。而对于震中的西边部分区域,库伦破裂应力变化较大,区域库伦破裂应力增量达到0.2 bar,2020-06-15 MW5.5地震就发生在这片区域,进一步表明此次土耳其MW5.9地震直接触发了MW5.5余震。图中其他库仑破裂应力增量超过0.1 bar的区域,尤其是北安那托利亚断层的部分区域,如果处于应力集中状态且接近临界值,那么未来其地震危险性都将大大提高。
本文利用Sentinel-1A宽幅数据获取2020-06-14土耳其地震的同震形变场,并采用BABO和有限断层方法反演断层几何形状参数和断层破裂滑动分布,得到以下几点认识:
1) 土耳其地震降轨LOS向同震形变场显示,断层北侧抬升,南侧沉降。降轨LOS向同震形变场最大沉降量约-7.75 cm,最大抬升形变量约8.87 cm,且关于东西向大致对称。
2) 断层几何参数反演显示,破裂断层走向约257.48°±0.65°,倾角约79.69°±0.98°,滑动角约154.2°±3.8°,震中位置为(40.754°E,39.389°N),通过与其他机构公布的震中比较得到,近场InSAR获取的LOS向形变场数据可以更好地确定中强地震的震中位置。
3) 断层滑动分布反演表明,破裂区域长度约8 km,宽度约6 km,破裂没有出露地表,破裂的最浅埋深约0.8 km,最大埋深约8.9 km,平均滑动角约146.8°。此次土耳其地震是一次以右旋走滑为主兼具少量逆冲特性的地震。最大滑动量约0.57 m,对应深度约4.278 km。假设地壳的剪切模量为3×105 bar,此次地震释放的地震矩约为4.54×1017 Nm,对应矩震级MW5.7,与AFAD提供的震级一致。
4) 对于ΔCFS≥0.1 bar的区域,未来危险性值得更多关注。
5) 此次土耳其地震受近东西向的潜伏断层控制,由于靠近北安那托利亚断层,可能还受到该断裂的共同作用。
致谢: 感谢欧洲空间局(ESA)提供Sentinel-1A卫星升降轨SAR数据。
[1] |
Şentürk S, Çakır Z, Ergintav S, et al. Reactivation of the Adyaman Fault (Turkey) through the MW5.7 2007 Sivrice Earthquake: An Oblique Listric Normal Faulting within the Arabian-Anatolian Plate Boundary Observed by InSAR[J]. Journal of Geodynamics, 2019, 131: 101 654 DOI:10.1016/j.jog.2019.101654
(0) |
[2] |
Taymaz T, Eyidoĝan H, Jackson J. Source Parameters of Large Earthquakes in the East Anatolian Fault Zone (Turkey)[J]. Geophysical Journal International, 1991, 106(3): 537-550 DOI:10.1111/j.1365-246X.1991.tb06328.x
(0) |
[3] |
Şengör A M C, Tüysüz O, Ímren C, et al. The North Anatolian Fault: A New Look[J]. Annual Review of Earth and Planetary Sciences, 2005, 33(1): 37-112 DOI:10.1146/annurev.earth.32.101802.120415
(0) |
[4] |
Reilinger R, McClusky S, Vernant P, et al. GPS Constraints on Continental Deformation in the Africa-Arabia-Eurasia Continental Collision Zone and Implications for the Dynamics of Plate Interactions[J]. Journal of Geophysical Research: Solid Earth, 2006, 111(B5)
(0) |
[5] |
Heimann S, Vasyura-Bathke H, Sudhaus H, et al. A Python Framework for Efficient Use of Pre-Computed Green's Functions in Seismological and other Physical Forward and Inverse Source Problems[J]. Journal of Geophysical Research: Solid Earth, 2019, 10(6): 1 921-1 935
(0) |
[6] |
谭凯, 赵斌, 张彩红, 等. GPS和InSAR同震形变约束的尼泊尔MW7.9和MW7.3地震破裂滑动分布[J]. 地球物理学报, 2016, 59(6): 2 080-2 093 (Tan Kai, Zhao Bin, Zhang Caihong, et al. Rupture Models of the Nepal MW7.9 Earthquake and MW7.3 Aftrershock Constrained by GPS and InSAR Coseismic Deformations[J]. Chinese Journal of Geophysics, 2016, 59(6): 2 080-2 093)
(0) |
[7] |
Rosen P A, Gurrola E, Sacco G F, et al. The InSAR Scientific Computing Environment[C]. 9th European Conference on Synthetic Aperture Radar, VDE, 2012
(0) |
[8] |
Farr T G, Rosen P A, Caro E, et al. The Shuttle Radar Topography Mission[J]. Reviews of Geophysics, 2007, 45(2)
(0) |
[9] |
Goldstein R M, Werner C L. Radar Interferogram Filtering for Geophysical Applications[J]. Geophysical Research Letters, 1998, 25(21): 4 035-4 038 DOI:10.1029/1998GL900033
(0) |
[10] |
Chen C W, Zebker H A. Phase Unwrapping for Large SAR Interferograms: Statistical Segmentation and Generalized Network Models[J]. IEEE Transactions on Geoscience and Remote Sensing, 2002, 40(8): 1 709-1 719 DOI:10.1109/TGRS.2002.802453
(0) |
[11] |
Yu C, Li Z H, Penna N T. Interferometric Synthetic Aperture Radar Atmospheric Correction Using a GPS-Based Iterative Tropospheric Decomposition Model[J]. Remote Sensing of Environment, 2018, 204: 109-121 DOI:10.1016/j.rse.2017.10.038
(0) |
[12] |
Jónsson S. Fault Slip Distribution of the 1999 MW7.1 Hector Mine, California, Earthquake, Estimated from Satellite Radar and GPS Measurements[J]. Bulletin of the Seismological Society of America, 2002, 92(4): 1 377-1 389 DOI:10.1785/0120000922
(0) |
[13] |
Meade B J. Algorithms for the Calculation of Exact Displacements, Strains, and Stresses for Triangular Dislocation Elements in a Uniform Elastic Half Space[J]. Computers and Geosciences, 2007, 33(8): 1 064-1 075 DOI:10.1016/j.cageo.2006.12.003
(0) |
[14] |
Wang R J, Martín F L, Roth F. Computation of Deformation Induced by Earthquakes in a Multi-Layered Elastic Crust: FORTRAN Programs EDGRN/EDCMP[J]. Computers and Geosciences, 2003, 29(2): 195-207 DOI:10.1016/S0098-3004(02)00111-5
(0) |
[15] |
Wang Q, Qiao X J, Lan Q G, et al. Rupture of Deep Faults in the 2008 Wenchuan Earthquake and Uplift of the Longmen Shan[J]. Nature Geoscience, 2011, 4(9): 634-640 DOI:10.1038/ngeo1210
(0) |
[16] |
Poyraz F. Uncertainty in the Determination of Fault Locking Depth and Strike Slip Rates by GNSS Measurements[J]. Tehniĉki Vjesnik, 2016, 23(1): 107-114
(0) |
[17] |
King G C P, Stein R S, Lin U J. Static Stress Changes and the Triggering of Earthquakes[J]. Bulletin of the Seismological Society of America, 1994, 84(3): 935-953
(0) |
[18] |
Ziv A, Rubin A M. Static Stress Transfer and Earthquake Triggering: No Lower Threshold in Sight?[J]. Journal of Geophysical Research: Solid Earth, 2000, 105(B6): 13 631-13 642 DOI:10.1029/2000JB900081
(0) |
2. Gravitation and Earth Tide, National Observation and Research Station, 196 Senlin Road, Wuhan 430000, China;
3. Key Laboratory of Earthquake Geodesy, CEA, 40 Hongshance Road, Wuhan 430071, China;
4. Hubei Earthquake Agency, 48 Hongshance Road, Wuhan 430071, China;
5. Wuhan Geomatics Institute, 209 Wansongyuan Road, Wuhan 430022, China