2017, Vol. 37
2. 中国科学院大学，北京市玉泉路甲19号，100049
2016-02-06 03：57(北京时间)台湾省高雄市美浓区(22.92°N，120.54°E)发生MW6.4地震。震后，全球多家地震监测机构(CWB、USGS、Harvard大学等)快速给出相关震源机制解[1]，认为此次地震为逆冲型地震，震源深度在17 ~ 26 km之间。主震后4 s，在震源西部较浅地区又发生一次破裂，有两个可能的破裂面，分别为东西走向的低角度北倾断层面274.81°/41.74°/17.02°(strike/dip/rake)，或是南北走向的垂直断层面171.95°/78.76°/130.47°(strike/dip/rake)[2]。GPS数据表明，震中西部主要表现为地表抬升以及向西运动，东部则表现为地表下降。Huang等[3]利用大地测量数据研究了该地震的形变场，通过联合升、降轨Sentinel-1A数据和ALOS-2数据获取震区东西方向和垂直方向的位移场。结果表明，在台南东部大部分地区地表向西运动，在新化和关庙附近局部地区地表表现为向东运动，在台南市和震中之间地表抬升(最大垂直位移可达13 cm)，震中东部表现为地表沉降；利用地震波和大地测量数据联合反演该地震震源机制认为，本次地震为双主震，分别由N61°W/15°NE/45°(strike/dip/rake)、深度为17 km的主断层和N8°E/30°/100°(strike/dip/rake)、深度为5~10 km的浅层断层引起。
InSAR技术具有空间分辨率高、覆盖范围广的优势，能够获得空间连续的地表形变场，国内外许多学者用其研究地震并获得成功[4-5]。本文主要利用一对升轨SAR影像获取美浓地区的空间形变场，将高精度GPS数据与InSAR观测结果进行对比，利用一次趋势面改正InSAR结果，并有效减弱InSAR形变场的解缠点误差和轨道误差等，能保持震区空间形变特征，获得更精确的InSAR形变场, 并快速、准确地获取震区地表变化和宏观震中位置，为理解地震形变过程、发生机理、未来地震演化趋势等提供重要的观测资料。
1 区域概况台湾岛位于欧亚板块和菲律宾海板块的交界处(图 1，白色方框为Sentinel-1A数据覆盖范围)。在台湾岛北部，菲律宾海板块自新生代以来一直朝西北移动，沿琉球海沟向北-西北方潜入欧亚板块；在南部，中国南海(欧亚板块)沿马尼拉海沟向东潜入菲律宾海板块。这两个板块沿着西北-东南方向以约80 mm/a的速度聚合[6]，并在东部海岸山脉、中央山脉、西部麓山带及西部海岸平原区形成一系列南北走向的褶皱断层[7]。这些断层具有很高的活动性，在历史上造成许多灾害性地震。图 1给出上世纪以来此次地震周围100 km范围内大于MW6.0的历史地震，其中包括1906年M7.9嘉义地震、1935年新竹-台中M7.1地震、1941年嘉义M7.1地震、1964年百禾M6.5地震、1999年集集M7.6地震[8]、2010年甲仙MW6.3地震[9]。
|
图 1 台湾地区地形、断层、余震、历史地震及影像覆盖范围 Fig. 1 Distribution of topography, faults, aftershocks, historical earthquakes and InSAR coverage in Taiwan province |
选用欧洲航天局Sentinel-1A卫星影像获取该地区同震形变场。选择震后9 d(2016-02-14)与震前3 d(2016-02-02)的一对升轨SAR数据组成干涉像对，其时间基线与垂直基线分别为12 d和24.26 m，并且震中靠近城市区域如美浓、旗山等(图 1)，保持较高的相干性，有利于获得高质量的干涉图。
选用瑞士GAMMA合成孔径雷达干涉测量软件进行处理[11-12]，所用数据为ESA提供的Level-1干涉宽幅模式(IW)SAR数据，采用D-InSAR技术获取覆盖台湾美浓地区的地表InSAR形变场，通过由粗到精的配准方法以及必要的迭代过程达到超高精度的配准要求。同时，使用90 m分辨率的SRTM DEM[13]消除地形相位。轨道数据采用ESA提供的精密轨道文件。为降低干涉相位噪声，采用依据数据功率谱的非线性滤波方法[14]对干涉图进行滤波，得到比较平滑的干涉条纹图。干涉条纹图中，每个条纹代表 2π弧度的变化量，C波段SAR影像在视线向(LOS)形变量约为2.8 cm[15]。为得到此次地震造成的视线向变化量，采用最小费用流(MCF)技术和不规则三角网(TIN)技术进行相位解缠[16]。由于宏观震中区域相干性较高，而东部山区相干性较低，为提高解缠速度与精度，并保留可靠区域的信息，对低相干区域进行掩膜处理，并且解缠点选择在地形平坦、相干性较高且形变量比较小的地面点。
3 结果与分析选用高精度的GPS观测结果[2]检验和定量分析InSAR形变场。该结果是对392个GPS连续观测站数据, 利用GAMIT/GLOBK 10.6软件和IGS快速星历的单日解，GPS震前测量时间为2016-01-29 00：00：00~ 02-05 19:57:00，震后时间为2016-02-05 19:58:00 ~02-09 23:59:30，采用ITRF2008参考框架。GPS水平速度场精度为1 ~ 3 mm，垂直位移场精度小于10 mm。图 2(蓝色箭头为GPS水平位移场)和图 3(红色和蓝色箭头为GPS垂直位移场，深红色三角形为校正InSAR结果的42个GPS点)给出震中附近120°~121.5°E、22.3°~23.5°N范围共138个GPS站的观测结果。
|
图 2 InSAR形变场和GPS水平位移场 Fig. 2 InSAR displacement field and GPS horizontal displacement field |
|
图 3 InSAR形变场和GPS垂直位移场 Fig. 3 InSAR displacement field and GPS vertical displacement field |
由于一维的InSAR形变结果无法与三维的GPS形变结果进行直接比较，需要把GPS观测结果投影到雷达视线向[17]。根据InSAR卫星与地面的空间几何关系，有：
|
(1) |
式中，α为卫星的方位角，θ为入射角。使用Sentinel-1A/IW数据，每一景影像有3个Sub-Swath，其中IW1的入射角为30.8°~36.6°，IW2为36.4°~41.1°，IW3为41.7°~45.4°。宏观震中在IW1覆盖范围内，选定θ为IW1入射角的平均值(33.7°)，可以得出地表E、N、U方向在雷达视线向的贡献为：
|
(2) |
根据InSAR观测结果的空间分布特征，选择震中附近均匀分布的同震GPS点(共42个)投影到雷达视线向。另外，InSAR是一种相对测量方式，为相对解缠点的变化量，而GPS为绝对测量，利用这两种数据评估形变场时，不能直接用于对比。因此，为便于两种不同数据的比较，以及获得连续分布、更为准确的地表形变信息，本文借助外部高精度GPS数据对InSAR结果进行一次趋势面改正[18-19]：
|
(3) |
式中，dGPS为GPS观测值在雷达视线向的投影，选择2σ以内的GPS观测值参与计算(共39个)；dSAR为InSAR观测值；a0为常数项；a1和a2为一次项；x和y为点位坐标。把GPS观测作为参考值，取GPS位置周围9×9像素的InSAR观测值进行平均，与GPS作差，利用最小二乘法对差值进行一次趋势面改正，得到图 4(红色表示抬升(LOS)，蓝色表示沉降(LOS)，圆点为InSAR和GPS在雷达视线向差值(dSAR-dGPS))所示的InSAR形变图。
|
图 4 改正后的视线向形变InSAR形变场 Fig. 4 Corrected InSAR displacement field |
经过一次趋势面改正后，InSAR形变场的形变量有一定程度的减小，但形变趋势与改正前基本一致，InSAR与GPS的相关系数由0.69提升到0.75。图 4中不同颜色的圆为改正后的误差值，其均值为0.3 cm，除少数误差大的点如GS80、GS76、GS53、LNCH等外，其他大部分点的误差在1 cm以内。总体来看，经过趋势面改正后的InSAR观测值与GPS观测值吻合较好，较好地反映了同震地表形变信息。对于存在较大差异的观测点(＞1 cm)，可能的原因是InSAR观测结果中每个像素大小约150 m×150 m，难以分辨细微的地表形变。
图 5给出跨越旗山断层和小岗山断层以及附近GPS点的3条剖面线, 分别为图 4中的虚线AA′、BB′、CC′。可以看出，GPS点与InSAR形变趋势一致，两者的相关系数分别为0.96、0.87、0.86；LNCH点与剖面线相距较远，可能导致其与InSAR差值较大，进一步表明利用InSAR技术能获得较好的空间形变场。
|
图 5 雷达视线向的剖面线 Fig. 5 Profilesof radar line-of-sight |
从形变图(图 4)可以看出，在宏观震中附近如震中西部和东部，以及剖面线形变最大处两边约5 km处有一些空白区域，由于该部分地区相干性较差，在数据处理过程中被剔除。从InSAR形变结果可以看出，地震造成以宏观震中为中心、方圆30 km的一个主要形变区域，该形变区域在雷达视线向表示为抬升，最大抬升量约12 cm(LOS)；在宏观震中的东部与南部地区，表现为沿雷达视线向远离雷达运动(约3 cm)。从整个形变图来看，总体形变量大约为-4 ~ 12 cm，在30~ 40 km区域内变化不到20 cm。另外，形变的空间连续性较好，说明宏观震中附近地表未发生明显破裂，与实地调查一致。对比分析可知，GPS与InSAR的整体运动趋势一致，InSAR经过校正后，两者的相关系数达到0.75，说明InSAR在监测地表形变方面可靠，但InSAR获得了更高空间分辨率的形变场，有利于整体把握形变场的运动趋势。
另外，使用InSAR技术获得的宏观震中与震中位置并不重合，而是位于震中西部约15 km处。震中为微观震中，是指根据地震波到时反演的震源初始破裂在地面上的投影；而宏观震中则对应着地表破裂能量释放最大点，表现为地表破坏或影响最严重的区域[17]，两者必定存在偏差。另外，宏观震中与微观震中的偏离程度受发震构造类型、震源深度、震级大小、场地效应、居民分布等因素的影响[18]。由于震源区有复杂的破碎地层和密集的断层带，地表形变最大处不一定位于震源的正上方。
从图 5可以看出，3条剖面线的形变趋势基本一致，剖面线呈现一个开口向下的二次曲线形状。以剖面线顶点为界，剖面线左边10 km以外有一段变化平缓的区域，平均变化量约为1.5 cm；剖面线两边各10 km范围内基本对称，并且剖面线右边的形变量要大于左边。在跨越两个断层的剖面线中，形变变化比较平缓，没有明显的突变。
4 结语本文使用欧洲航天局Sentinel-1A卫星，采用D-InSAR技术获取台湾省美浓地震的地表空间分布特征。为获得高精度、连续分布的地表形变场，将高精度GPS观测值投影到雷达视线向，对InSAR形变场进行一次趋势面校正。结果表明，经过趋势面校正后的InSAR观测值与GPS观测值吻合较好，形变主要发生在台南市东南部约17 km处，形变范围约30 km×30 km，主要表现为沿雷达视线向抬升，最大抬升量约12 cm。InSAR技术可以快速确定宏观震中及地表变化范围，但尚不能分辨出地表运动的具体方向。
致谢: 欧洲航天局提供Sentinel-1A数据，台湾“中央”研究院地球科学研究所许雅儒研究员提供GPS数据，在此表示衷心感谢！
| [1] | |
| [2] | |
| [3] |
Huang M H, Tung H, Fielding E J, et al. Multiple Fault Slip Triggered above the 2016 MW6.4 Meinung Earthquake in Taiwan[J]. Geophysical Research Letters, 2016, 43(14): 7 459-7 467 DOI:10.1002/2016GL069351
(  0) |
| [4] |
Massonnet D, Rossi M, Carmona C, et al. The Displacement Field of the Landers Earthquake Mapped by Radar Interferometry[J]. Nature, 1993, 364(6 433): 138-142
( 0) |
| [5] |
Fialko Y, Sandwell D, Simons M, et al. Three-Dimensional Deformation Caused by the Bam, Iran, Earthquake and the Origin of Shallow Slip Deficit[J]. Nature, 2005, 435(7 040): 295-299
( 0) |
| [6] |
Send T. A Model for the Motion of the Philipine Sea Plate Consistent with NUVELWI and Geological Data[J]. Journal of Geophysical Research Atmospheres, 1993, 98(B10): 17 941-17 948 DOI:10.1029/93JB00782
( 0) |
| [7] |
Yu S B, Chen H Y, Kuo L C. Velocity Field of GPS Stations in the Taiwan Area[J]. Tectonophysics, 1997, 274(1): 41-59
( 0) |
| [8] |
Ma K F, Lee C T, Tsai Y B, et al. The Chi-Chi, Taiwan Earthquake: Large Surface Displacements on an Inland Thrust Fault[J]. EOS Transactions American Geophysical Union, 1999, 80(50): 605-611
( 0) |
| [9] |
Huang M H, Dreger D, Bürgmann R, et al. Joint Inversion of Seismic and Geodetic Data for the Source of the 2010 March 4, MW6.3 Jia-Shian, SW Taiwan, Earthquake[J]. Geophysical Journal International, 2013, 193(3): 1 608-1 626 DOI:10.1093/gji/ggt058
( 0) |
| [10] |
Wegmüller U, Werner C, Strozzi T, et al. Sentinel-1 Support in the GAMMA Software [C].Fringe 2015 Conference, ESA ESRIN, Frascati, Italy, 2015 http://adsabs.harvard.edu/abs/2015ESASP.731E..33W
( 0) |
| [11] |
Wegmuller U, Werner C, Strozzi T. SAR Interferometric and Differential Interferometric Processing Chain[C].International Geoscience and Remote Sensing Symposium, 1998
( 0) |
| [12] |
Farr T G, Rosen P A, Caro E, et al. The Shuttle Radar Topography Mission[J]. Reviews of Geophysics, 2007, 45(2): 1-33
( 0) |
| [13] |
Goldstein R M, Werner C L. Radar Interferogram Filtering for Geophysical Applications[J]. Geophysical Research Letters, 1998, 25(21): 4
( 0) |
| [14] |
廖明生, 林珲. 雷达干涉测量:原理与信号处理基础[M]. 北京: 测绘出版社, 2003 (Liao Mingsheng, Lin Hui. Radar Interferometry: Basic Principle and Signal Processing[M]. Beijing: Surveying and Mapping Press, 2003)
( 0) |
| [15] |
Werner C, Wegmuller U, Strozzi T, et al. Interferometric Point Target Analysis for Deformation Mapping[C].Geoscience and Remote Sensing Symposium, 2003
( 0) |
| [16] |
Hanssen R F. Radar Interferometry: Data Interpretation and Error Analysis[J]. Journal of the Graduate School of the Chinese Academy of Sciences, 2001, 2(1): 577-580
( 0) |
| [17] |
中国地震局.GB/T18207.1-2000防震减灾术语第一部分: 基本术语[S].北京: 中国标准出版社, 2000 (CEA.GB/T18207.1-2000 Terminology for Earthquake Prevention and Disaster Reduction(PartⅠ): Basic Terminology[S].Beijing : Standards Press of China, 2000)
( 0) |
| [18] |
兰晓雯, 田家勇, 谢周敏, 等. 发震断层对微、宏观震中偏离的影响分析:以新疆地区为例[J]. 现代地质, 2011, 25(1): 114-121 (Lan Xiaowen, Tian Jiayong, Xie Zhoumin, et al. Analysis on the Impact of Seismogenic Fault on the Micro- and Macro-Epicenter Deviation Taking Xinjiang Area as an Example[J]. Geoscience, 2011, 25(1): 114-121 DOI:10.3969/j.issn.1000-8527.2011.01.015)
( 0) |
| [19] |
Liu C S, Huang I L, Teng L S. Structural Features of Southwestern Taiwan[J]. Marine Geology, 1997, 137(3): 305-319
( 0) |
2. University of Chinese Academy of Science, A 19 Yuquan Road, Beijing 100049, China