2021, Vol. 41
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地震。
|
红色矩形代表T50降轨范围, 3个红色震源机制解分别由GCMT、AFAD、USGS提供, 蓝色震源机制解表示2020-06-15发生的MW5.5余震 图 1 构造地震背景 Fig. 1 Tectonic and seismic background |
图 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所示。
|
表 1 Sentinel-1A的详细参数 Tab. 1 Detailed parameters of sentinel-1A |
使用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)。
|
图 2 土耳其地震降轨的同震干涉图 Fig. 2 Coseismic interferogram of descending track of Turkey earthquake |
|
图 3 降轨LOS向形变场GACOS大气校正 Fig. 3 GACOS atmospheric correction for LOS deformation field of descending track |
从图 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向形变场的形变中心附近。
|
图 4 LOS向形变场的截面分析 Fig. 4 Cross section analysis of LOS deformation field |
为加快反演进程，采用四叉树方法[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，最终收敛状况都较好，且失配值随着迭代的不断进行而不断地降低。
|
图 5 参数序列图 Fig. 5 Sequence plots for all parameters |
图 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°。
|
红色曲形实线为高斯核密度; 红色虚线为参数的平均值; 90% 置信区间表示为中间的红色阴影区, 阴影区最宽边界代表最小、最大值; 黑色实线表示GCMT参数 图 6 2020年MW5.7土耳其地震断层模型参数的分布 Fig. 6 Parameter distribution of the 2020 MW5.7 Turkey earthquake fault |
|
|
表 2 反演的断层几何形状参数 Tab. 2 Inversion of fault geometry |
在反演断层几何参数时，得到的残差结果如图 7所示。图中圆点表示采样点，矩形表示反演断层所在位置，矩形中的黑点表示震中，断层左上角点经纬度为(40.780 4° E，39.390 8° N)，对应深度为1.308 km。在靠近震中附近，LOS向视向线模型值都能够较好地拟合观测值，降轨数据残差范围是-0.018~0.012 m，残差的RMS为6 mm，总体来看模拟效果较好。
|
图 7 断层几何参数反演获得的视线向形变残差 Fig. 7 LOS deformation residuals of fault geometry inversion |
基于上述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。
|
图 8 利用InSAR数据反演的土耳其地震滑动分布 Fig. 8 Coseismic slip distribution of Turkey earthquake inverted by InSAR |
|
图 9 滑动分布模型的形变残差 Fig. 9 The deformation residuals of the distributed slip model |
根据表 3，InSAR反演的震中与CENC测定的震中最接近，相差1.3 km。结合图 4分析，InSAR得到的震中更符合实际震中特征，更为精确，且更接近最大形变值所处的位置，而与USGS、AFAD测定的震中在平面空间上分别相差约4.1 km、4.4 km，与GCMT测定的震中距离最远，相差约7.7 km，这可能主要受测震台网的分布、测震站台的位置以及地震震级(较小)等影响。近场数据测定的震中明显要比USGS、GCMT更加准确，这表明用近场InSAR形变场数据可以更好地确定中强地震的震中位置。
|
|
表 3 震源参数比较 Tab. 3 Focal parameters comparison |
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的区域，尤其是北安那托利亚断层的部分区域，如果处于应力集中状态且接近临界值，那么未来其地震危险性都将大大提高。
|
绿色五角星表示本文反演的震中 图 10 深度位于4.278 km时的库伦破裂应力变化 Fig. 10 The Coulomb failure stress (CFS) change at a depth of 4.278 km |
本文利用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