2017, Vol. 37
2. 西南交通大学地球科学与环境工程学院，成都市高新区西部园区，611756
2015-07-03 09:07:46新疆维吾尔自治区皮山县发生MW6.5强震，造成严重的生命财产损失[1-3]。此次地震震中(37.5°N，78.1°E)位于塔里木盆地西南部与青藏高原西北部的交界处，是由于印度板块向北推进和塔里木盆地的俯冲所致，发震断层属于阿尔金断裂带系统[3-4](2008年MW7.1新疆于田地震也位于该断裂带[5])。然而，根据《中国活断层构造地图》[4]，在皮山地震约50 km范围内并未发现活断层的存在。因而，研究皮山地震断层参数，探索隐伏断层位置所在，对于理解塔里木盆地地质构造具有重要意义。
皮山地震发生后，美国地质勘探局(USGS)、全球矩心矩张量项目(global centroid moment tensor, GCMT)、中国地震台网的矩心矩张量中心(centroid moment tensor, CMT)以及Wen等[1-3]，均公布了皮山地震的震源机制解。结果表明，该地震是一次逆冲型为主事件。但是，不同机构由于数据源的不同以及构建模型的差异，所得震源机制、断层滑动分布结果等各有不同，并给出了不同的地学解释。
2014-04-03欧空局Sentinel-1A卫星发射升空[6-7]，其对同一区域的重访周期仅为12 d，极大地提升了连续地表雷达数据的获取能力。特别是该卫星的干涉宽幅成像模式(interferometric wide swath mode, IW)，十分有助于对地震完整地InSAR覆盖[7]。因此，本文基于Sentinel-1A IW SAR数据，采用断层几何参数非线性反演和断层滑动参数线性反演相结合的方法，获得2015皮山地震的震源机制解。
1 Sentinel-1A SAR影像干涉处理塔里木盆地位于天山山脉以南，昆仑山脉以东，阿尔金山脉以北，西邻帕米尔高原。2015皮山地震发生于塔里木盆地西南边界，发震断层属于阿尔金断裂带，震源深度为23 km(USGS)，是中国最大的内陆盆地地震[1-4]。发震的动力学诱因为印度板块向北推进，帕米尔高原、昆仑造山带抬升，塔里木盆地俯冲，致使断裂带应力持续积累并达到临界状态，最终积蓄能量从皮山断层释放并爆发本次地震。
皮山地震引发了较强的地表形变，欧空局的Sentinel-1A卫星震后及时完成了对地观测任务，获取了覆盖震区的多幅SAR影像，并免费公布于科学数据中心网站[8]，为后续研究提供了丰富的数据源。本文获取了升轨飞行模式下的皮山地震震前震后两幅SAR影像，其影像参数基本信息如表 1所示，相应的影像覆盖范围及区域地质背景如图 1所示。图 1中黑色矩形框为Sentinel-1A SIW3(第3个子幅宽)数据覆盖范围，地表形变区域包含于该幅宽内。黄色五角星显示皮山地震震中位置(USGS)，红色圆圈标明震后两个月内发生的余震(MW4.0~5.0)震中位置，黑色断层线为活断层，紫色断层线为早中更新世隐伏断层[1, 4]，本次地震发生于隐伏断层附近。
|
表 1 Sentinel-1A卫星SAR影像参数信息 Tab. 1 The parameters of SAR images from Sentinel-1A satellite |
|
图 1 皮山地震Sentinel-1A卫星SAR影像覆盖范围 Fig. 1 The coverage of SAR images from Sentinel-1A satellite of Pishan earthquake |
基于欧空局公布的S1toolbox软件，使用震前震后一对SAR影像，采用二轨差分干涉方法[9-10]提取皮山地震同震形变场。欧空局提供的S1A精轨数据用于估算轨道状态矢量。采用90 m分辨率的SRTM v4.0数字高程模型[11]来扣除地形相位。在方位向和距离向采用5 :14的多视比例来提高信噪比、降低噪声。最终采用Snaphu软件[12]解缠相位获得沿雷达视线方向(line-of-sight, LOS)形变，如图 2所示，图中黑色实线标明皮山地震断层位置，黑色虚线框为断层地表投影(断层几何参数由后续非线性反演结果获得)。
图 2所示为完整的皮山地震主形变区，相干系数小于0.3的部分(噪声湮没)在干涉图中未予显示。由图 2(b)可见，条纹分布中基本无残留轨道条纹，说明Sentinel-1A卫星具有较高的定轨精度。图 2(a)结果说明，皮山地震沿雷达视线方向(LOS向)最大抬升约为0.11 m。但由于同震形变场远场失相干严重，必然存在解缠基准偏差，致使形变场获取不精确，该偏差将在后续反演中予以估算。
|
图 2 vSentinel-1A提取的皮山地震同震形变场 Fig. 2 The coseismic deformation field extracting from Sentinel-1A of Pishan earthquake |
本文以Sentinel-1A卫星获取的皮山地震同震形变场作为反演约束，对本次地震的断层几何参数和滑动参数开展详细分析与研究。为保证反演精度且降低反演计算量，采用由数据分辨率约束的四叉树方法[13]对同震形变场(图 2(a))进行降采样处理，最终获得形变场近场2 257个数据点作为后续反演的输入数据。
采用两步反演方法进行断层参数估计。首先，综合USGS和Wen等[2-3]已有皮山地震断层参数研究成果，假设断层表面均匀滑动的前提下，采用模拟退火算法[14]对断层参数进行非线性搜优，获得断层位置(经度和纬度)、深度、长度、宽度、走向、倾角7个几何参数(表 2)；其次，基于Okada弹性半空间位错模型[15]反演断层表面滑动分布。
|
|
表 2 皮山地震断层几何参数 Tab. 2 Fault geometry parameters of Pishan earthquake |
为进一步确定皮山地震的断层滑动空间分布，在确定最优断层几何参数的基础上，将断层平面沿走向和倾向进行40 km×38 km的格网离散化，获得1 520个子断层。同时，为避免子断层块间滑动存在过大跳变，通常采用拉普拉斯算子[16]对断层滑动分布进行平滑约束。需要特别指出，由于Sentinel-1A IW模式下SAR影像幅宽较宽，使得入射角产生从34.3°到44.3°的变化，简单使用平均入射角进行反演会引起误差。因此，对上述2 257个数据点，在均采用对应位置入射角的基础上，考虑前文解缠基准偏差的情况，构建反演模型：
|
(1) |
式中，G为格林函数，λ为平滑因子，L为拉普拉斯算子，D为形变观测值。等式左边补充0和1矩阵为估算基准偏差。m包含待求解的走滑和倾滑分量。
基于上述模型在非负最小二乘原则下进行断层滑动量的反演计算，获得皮山地震的断层滑动分布结果如图 3所示。结果表明，本次地震为小倾角逆冲型地震兼具少量走滑分量，最大滑动量为0.76 m，断层位于地下约5 km深度，未破裂到地表，与中国地震局地质研究所调查结果一致；滑动集中区位于地表以下9~14 km，与Wen等[3]基于大地测量数据反演结果较为一致。且该模型下估算的基准偏差约为3.77 cm，累计释放地震矩达5.14×1018，约合矩阵级MW6.47，与USGS[2]和Wen[3]等结果均具有较高的一致性，反演结果较为可靠。
|
图 3 皮山地震断层滑动分布结果 Fig. 3 Fault slip distribution of Pishan earthquake |
基于图 3的断层滑动分布结果，使用Okada模型正演计算了沿Sentinel-1A卫星LOS向模拟形变场，结果如图 4所示。与图 2 Sentinel-1A提取的皮山地震同震形变场对比可以发现，模拟结果与观测结果在形变分布上具有较好的一致性，但两者在形变量上存在略微差异(图 4(a))，且模拟结果具有更加平滑的条纹(图 4(b))。
|
图 4 皮山地震模拟同震形变场 Fig. 4 The simulated coseismic deformation field of Pishan earthquake |
为进一步定量比较观测结果与模拟结果，计算两者的残差，结果如图 5所示。由图 5(a)发现，最大的残差大于8 cm；且从图 5(b)也可以发现，A区域和B区域存在残余条纹，可能是残余的大气相位所致。统计平均残差约为4 cm，与模型所估计的解缠基准偏差3.77 cm较为一致，因此上述反演模型还是较好地给出了皮山地震断层滑动分布。
|
图 5 皮山地震观测结果与模拟结果残差 Fig. 5 Residual results between observation and simulated deformation of Pishan earthquake |
本文基于Sentinel-1A卫星观测的一对震前震后SAR影像，提取了2015皮山地震同震形变场。并在考虑渐变入射角和解缠基准偏差的情况下，构建断层滑动反演模型，获取了本次地震断层滑动分布。结果表明，皮山地震是小倾角逆冲型地震，断层破裂于地下约5 km深度；最大滑动量约在地下11 km处，为0.76 m；高滑区集中于地表以下9~14 km。滑移模型估算的解缠基准偏差约为3.77 cm，与平均统计残差4 cm较为一致，说明InSAR形变场的残差大多由解缠基准偏差所致，剩余残差可能是大气残留相位。皮山地震累计释放地震矩达5.14×1018，约合矩阵级MW6.47。
致谢: 感谢欧空局在科学数据中心网站公布Sentinel-1A卫星雷达数据，感谢Wessel和Smith(1998)提供Generic Mapping Tools(GMT)绘图软件。
| [1] |
CENC.[EB/OL].http://www.cea.gov.cn/publish/dizhenj/468/553/101527/index.html.
(  0) |
| [2] |
USGS.[EB/OL].http://earthquake.usgs.gov/earthquakes/eventpage.
( 0) |
| [3] |
Wen Y M, Xu C J, Liu Y, et al. Deformation and Source Parameters of the 2015 MW6.5 Earthquake in Pishan, Western China, from Sentinel-1A and ALOS-2 data[J]. Remote Sensing, 2016, 134(8): 1-14
( 0) |
| [4] |
Deng Q. Map of Active Tectonics in China[M]. Beijing: Seismological Press, 2007
( 0) |
| [5] |
Zhou Y, Wang W M, Xiong L, et al. Rupture Process of 12 February 2014, Yutian MW6.9 Earthquake and Stress Change on nearby Faults[J]. Chin J Geophys, 2015, 58(1): 184-193
( 0) |
| [6] | |
| [7] |
李永生, 冯万鹏, 张景发, 等. 2014年美国加州纳帕MW6.1地震断层参数的Sentinel-1A InSAR反演[J]. 地球物理学报, 2015, 58(7): 2339-2349 (Li Yongsheng, Feng Wanpeng, Zhang Jingfa, et al. Coseismic Slip of the 2014 MW6.1 Napa, California Earthquake Revealed by Sentinel-1A InSAR[J]. Chinese Journal of Geophysics, 2015, 58(7): 2339-2349)
( 0) |
| [8] |
ESA.[EB/OL].https://scihub.copernicus.eu/s2/#/home.
( 0) |
| [9] |
Massonnet D, Rossi M, Carmona C. The Displacement Filed of the Landers Earthquake Mapped by Radar Interferometry[J]. Nature, 1993, 364(8): 138-142
( 0) |
| [10] |
冯万鹏, 许力生, 李振洪. 2008年10月当雄MW6.3级地震断层参数的InSAR反演及其构造意义[J]. 地球物理学报, 2010, 53(5): 1134-1142 (Feng Wanpeng, Xu Lisheng, Li Zhenhong. Fault Parameters of the October 2008 Damxung MW6.3 Earthquake from InSAR Inversion and Its Tectonic Implication[J]. Chinese Journal of Geophysics, 2010, 53(5): 1134-1142 DOI:10.3969/j.issn.0001-5733.2010.05.014)
( 0) |
| [11] |
Farr M, Rosen A, Caro E. The Shuttle Radar Topography Mission[J]. Reviews of Geophysics, 2007, 45(2): 33
( 0) |
| [12] |
Chen C W, Zebker H A. Two-Dimensional Phase Unwrapping with Use of Statistical Models for Cost Functions in Nonlinear Optimization[J]. Journal of the Optical Society of America A: Optics, Image Science and Vision, 2001, 18(S2): 338-351
( 0) |
| [13] |
Lohman R B, Simons M. Some Thoughts on the Use of InSAR Data to Constrain Models of Surface Deformation: Noise Structure and Data Downsampling[J]. Geochem Geophys Geosyst, 2005, 6(1): 359-361
( 0) |
| [14] |
Clarke P J, Paradissis D, Briole P, et al. Geodetic Investigation of the 13 May 1995 Kozani-Grevena (Greece) Earthquake[J]. Geophysical Research Letters, 1997, 24(6): 707-710 DOI:10.1029/97GL00430
( 0) |
| [15] |
Okada Y. Surface Deformation Due to Shear and Tensile Faults in a Half-Space[J]. Bulletin of the Seismological Society of America, 1985, 75(4): 1135-1154
( 0) |
| [16] |
Harris R A, Segall P. Detection of a Locked Zone at Depth on the Parkfield, California, Segment of the San Andreas Fault[J]. Journal of Geophysical Research: Solid Earth (1978-2012), 1987, 92(B8): 7945-7962 DOI:10.1029/JB092iB08p07945
( 0) |
2. Faculty of Geosciences and Environmental Engineering, Southwest Jiaotong University, West High-Tec Zone, Chengdu 611756, China