2. 武汉大学地球空间环境与大地测量教育部重点实验室，湖北 武汉 430079；
3. 格拉斯哥大学地理与地球科学学院，英国 格拉斯哥 G12 8QQ；
4. 国家基础地理信息中心大地测量部，北京 100048
2. Key Laboratory of Geospace Environment and Geodesy, Ministry of Education, Wuhan University, Wuhan 430079, China;
3. School of Geographical and Earth Sciences, University of Glasgow, Glasgow G12 8QQ, UK;
4. Department of Geodesy, National Geomatics Center of China, Beijing 100048, China
1 引 言
北京时间2007年5月5日16时51分，在西藏阿里地区日土县与改则县交界处发生了Mw 6.1级地震(本文称为阿里地震)，地震震中位于81.99°E，34.26°N，震源深度为7.0 km(来自美国地质调查局(United States Geological Survey,USGS)。由于此次地震所在地的大部分区域为无人区，没有造成人员伤亡。该地震是自2005年以来发生在青藏高原一系列正断层事件(主要包括2005年仲巴Mw 6.2级地震、2008年改则Mw 6.4级地震与Mw 5.9级地震、2008年于田Mw 7.1级地震、2008年仲巴Mw 6.7级地震和2008年当雄Mw 6.3级地震等)中震级较小的一个，可能进一步佐证了青藏高原板块内部在重力作用下东西向扩张的事实[1]。
2007年阿里地震震区位于班公错—怒江断裂带以北活动构造带上，该构造带主要由发育较好的NE和NW向活动断裂组成，如托和平错断裂、美马错—德克玛洛断裂和鲁玛江冬错断裂等(图 1)。这一区域的隆起带隆升较强烈，呈单斜式隆升，并且断陷带也较发育。其中，第四纪冲、洪积扇分布广泛并有较大湖泊分布[2]。在过去的30年中，2007年阿里地震发生前，该地区的地震活动并不强烈。据USGS地震目录记载，在1975年至2005年的30年间，主震震中1°范围内发生的最大地震为1985年的Mw 5.0级地震。因此，研究这样一个少有的中强地震的同震位移场与断层参数对认识该地区的构造运动特征具有重要意义。
|
| 图 1 2007年阿里地震震中区域地质构造背景 Fig. 1 Tectonic setting of epicenter region of the 2007 Ali earthquake |
此外，由于阿里地震发生在藏北无人区，野外地质调查和地球物理数据的采集都非常困难，对它的研究起初只能依赖于全球地震台网(GSN)记录的宽频带远震体波数字资料，据此给出粗略的矩张量解和最佳双力偶解。而作为新兴的空间大地测量技术的InSAR技术，能够获取高精度、高分辨率的地表形变信息，且不需要提供地面控制点，非常适合偏远地区的地震形变监测[3, 4,5]。需要特别指出的是，InSAR仅能观测到一维形变量，它获取到的是雷达波视线向(line-of-sight，LOS)的形变，是地表三维形变(通常指东西向、南北向和垂向形变)在视线向上的投影(图 2)。在绝大多数情况下，这种一维形变不能完全反映地表的真实形变，有时候甚至会引起偏差。不同成像几何条件下的干涉数据(如具有不同成像几何特征的升、降轨数据)能够提供更为丰富的形变信息，可以有效地降低视线向的模糊问题[6]，从而对地震震源机制和破裂过程的反演等提供更可靠的数据约束。本文利用EnviSat卫星的升、降轨数据来获取2007年阿里地震的高质量同震地表形变场，然后采用线性和非线性反演方法来确定地震的断层参数和滑动分布，反演结果将深化对阿里地震的成因和机理的认识。
|
| 图 2 右视SAR的升、降轨成像几何 Fig. 2 The geometry of the ascending and descending right-looking SAR acquisitions |
地震发生前后，欧洲航天局(European Space Agency，ESA)的EnviSat卫星C波段ASAR传感器对阿里地震进行了观测。为了研究本次地震的几何结构和运动学特征，笔者通过中国科技部和欧洲航天局联合资助的“龙计划”三期国际合作项目收集了覆盖此次地震发生区域的雷达影像，为了尽可能获取高相干数据和减小非同震形变信号的影响，最终选取了两组像对(其中升、降轨各一组)来进行干涉处理(表 1)。
| Track | 主影像 获取时间 | 从影像 获取时间 | 垂直基线 /m | 中误差 /mm | 衰减距离 /km | |
| IP1 | T155A | 2007-04-13 | 2007-06-22 | 365 | 3.8 | 2.3 |
| IP2 | T205D | 2007-02-06 | 2007-07-31 | 94 | 4.3 | 2.6 |
干涉处理使用的平台为瑞士GAMMA软件[7]。基于该软件，采用二通法[8]对获取的卫星SAR影像进行差分干涉处理来获取阿里地震的同震地表形变场。在数据处理过程中，使用ESA提供的DOR精密轨道和NASA提供的90 m分辨率的SRTM DEM数据[9]来去除地形相位的影响。由于干涉像对的空间基线较长，因此尽管使用了较为精确的DOR轨道数据，但是在得到的干涉图中仍然存在由于残余轨道误差所引起的残余相位影响。对于这个影响，本文采用多项式模型来拟合去除。同时为了降低干涉相位的噪声水平，提高干涉图的信号质量，采用基于能量谱的局部自适应滤波对干涉图进行滤波。最终，采用枝切法来解缠得到差分干涉相位。阿里地区地处高寒山区，平均海拔为5000 m，气候干燥，大气影响相对较小。尽管如此，为了进一步削弱干涉相位中的大气相位贡献，还根据地形与大气相位的相关性拟合去除了与地形相关的大气相位[10]。最终得到了地理编码后的升、降轨阿里地震的高精度同震地表形变场(图 3)。
|
| 图 3 2007年阿里地震的升轨T155A(a)和降轨T205D(b)同震地表形变场 Fig. 3 Coseismic interferograms associated with the 2007 Ali earthquake from (a) ascending track T155A and (b) descending track T205D |
从同震形变场(图 3)中可以看出，该地震引起的地表形变总体表现为椭圆状分布，其NNW向长约40 km，NEE向长约20 km。除了部分被冰雪和湖泊所覆盖的区域外，整个干涉形变相位连续，条纹光滑清晰，特征明显。在T155A的形变图(图 3(a))中，视线向(LOS)的最大位移达到5.8 cm；而在T205D(图 3(b))中，LOS的最大位移为6.7 cm。由断层两盘的相对运动模式估计的断层运动模式与震源机制解相一致。
在InSAR同震地表形变场中，主要包含有轨道误差、大气误差、DEM误差、热噪声和数据处理中的配准、重采样误差等多种误差，这些误差对形变相位的影响方式和量级各不相同。为了估计InSAR同震形变场的精度水平，本文采用1D方差-协方差函数[11]：Ci,j=σ2e-r/α来描述形变相位的误差特征，其中Ci,j为相距r的像素i和j之间的协方差函数，σ2是整个形变场的方差，α是误差衰减距离。据此给出的升、降轨InSAR同震形变场的中误差分别为3.8 mm和4.3 mm，方差-协方差衰减距离分别为2.3 km和2.6 km(表 1)。该结果略优于2008年青海大柴旦地震的同震地表形变场的精度水平[4]，接近InSAR自身的精度水平，这可能也意味着在所得到的升、降轨InSAR地表形变场(图 3)中主要包含的是地震同震形变信息。
3 地震反演模型为了提高反演效率，本文采用基于格林函数分辨率的降采样法[12]来对同震形变场(图 3)进行降采样，并按照得到的采样点位置来计算实际的卫星入射角及其轨道方位角，最终得到了1060个观测数据来反演震源参数，其中包括503个升轨InSAR观测数据和557个降轨InSAR观测数据。
为了获取地震的震源参数，本文采用两步法策略[13,14]来进行反演分析。首先是基于矩形位错模型[15]采用非线性反演方法确定断层的几何结构；然后通过线性反演方法来估计断层面上的精细滑动分布。这里特别需要指出的是，对于线性反演而言，基于非线性反演给出的断层几何结构(尤其是倾角)也许并不是最优[16]，它需要在线性反演中重新进行估计。
本文采用PSOKINV程序[17,18]来进行反演，该程序采用多峰值颗粒群(MPSO)优化算法来确定断层几何参数，如断层位置、长度、深度、走向、倾向、滑动角和滑动量等参数，具有收敛效率高、控制参数少和程序简洁易行等特点。最佳拟合的均匀断层模型参数见表 2，表明发震断层是一个以正倾滑为主兼有少量右旋走滑的断层，其走向近NW-SE方向，倾角为39.2°±5.0°。断层的深度为6.0 km，意味着该发震断层并没有出露于地表。反演得到的地震矩为1.05×1018 Nm，小于GCMT给出的1.54×1018 Nm，但大于USGS给出的6.46×1017 Nm。
| 经度 /(°) | 纬度 /(°) | 走向 /(°) | 倾向 /(°) | 滑动角 /(°) | 深度 /km | 长度 /km | 宽度 /km | 滑动量 /m | Mw | |
| USGS-BW | 81.990 | 34.260 | 153 | 46 | -131 | 7.0 | - | - | - | 5.8 |
| GCMT | 81.97 | 34.33 | 130 | 88 | 178 | 23.9 | - | - | - | 6.1 |
| 均匀滑动 | 81.924 | 34.307 | 158.4 | 39.2 | -122.6 | 6.0 | 24.3 | 8.25 | 0.16 | 5.98 |
| ±0.5 km | ±0.7 km | ±2.4 | ±5.0 | ±5.1 | ±0.8 | ±1.5 | ±0.5 | ±0.05 | ||
| 分布式滑动 | 81.924 | 34.307 | 158.4 | 43 | -131.7 | 0 | 30 | 25 | 0.14 | 6.03 |
在线弹性位错模型中，当断层的几何结构确定后，断层面上的滑动量(走滑和倾滑分量)与地表形变之间呈线性关系。在确定断层精细滑动分布模型过程中，本文采用之前均匀滑动模型给出的断层位置和走向等参数，但对其中的断层倾角进行重新估计。将断层长度延长到30 km，断层面宽度延长到25 km，然后将断层面离散成1 km×1 km大小的750个断层片。这些断层片上的滑动量与观测值之间的关系为
式中，G为格林函数；κ2为光滑因子；L为拉普拉斯二阶平滑算子[19]；D为InSAR观测值；S为待求的滑动量。
由于均匀滑动模型给出的断层倾角并非是线性反演模型所对应的最优倾角[16,17]，这里采用综合考虑模型粗糙度和模型拟合残差的方法[18]确定最优断层倾角和光滑因子(κ2)。在该方法中，定义目标函数lg(ε+ψ)，其中ε为模型拟合残差，ψ为滑动分布模型的粗糙度。图 4为断层倾角、光滑因子(κ2)和目标函数lg(ε+ψ)的关系图，文中选取目标函数最小位置处的断层倾角(43°)和光滑因子(κ2=2)作为线性反演模型最优断层倾角和光滑因子。
|
| 图 4 断层倾角、光滑因子与函数lg(ε+ψ)之间的关系(红色五角星为全局最优位置) Fig. 4 Contour map of lg(ε+ψ) with variations of dips and smoothing factors(The red star indicates the point of global minimum) |
线性反演给出的阿里地震同震滑动分布模型见图 5(a)。从图 5(a)中可以看出，同震滑动分布主要集中在7～12 km深度位置(对应10～18 km宽度位置)，其平均滑动量为0.14 m，平均滑动角为-131.7°；最大滑动量位于9 km深度位置，为0.31 m，滑动角为-111.2°。地震释放的能量为1.24×1018 Nm，相当于矩震级Mw 6.03。为了估计线性反演模型所给出滑动分布的精度水平，文中采用蒙特卡洛误差传递方法[20]来估计模型的精度，得到的同震滑动的误差分布见图 5(b)。从图 5(b)中可以看到，线性反演的同震滑动的误差分布比较均匀，其平均误差为2 cm，最大误差为5.6 cm。
|
| 图 5 2007年阿里地震同震滑动分布(a)及其误差(b) Fig. 5 Coseismic slip distribution (a) and its errors (b) of the 2007 Ali earthquake |
图 6显示的是采用线性反演给出的同震滑动分布模型模拟的升、降轨干涉图以及残差分布，拟合结果表明该滑动分布模型能较好地解释阿里地震的InSAR同震地表形变场，模型生成的干涉图非常清晰，残差结果也很小，观测值与拟合值之间的相关性为96.4%。T155A和T205D的残差均方误差分别为2.7 mm和3.0 mm，与InSAR同震地表形变场的精度相当。
|
| 图 6 分布滑动模型拟合的升、降轨形变场((a)和(b))和残差形变场((c)和(d)) Fig. 6 Modeled (upper left, upper right) and residual (bottom left, bottom right) ascending and descending displacement maps |
2007年5月5日阿里Mw 6.1级地震发生在藏北无人区，对人类的生产和生活造成的影响非常小，但是它对理解藏北地区的运动特征有着重要的意义。尽管该地区的恶劣自然条件使野外地质考察难以完全展开，但是本研究同时利用升、降轨道EnviSat卫星SAR数据获取了该地震的高精度同震地表形变场。在此基础上，分别采用MPSO非线性和最小二乘线性方法反演给出了地震的最佳断层几何模型和均匀滑动分布，为研究地震破裂对周边断层的应力状态的影响、震后形变机理分析和灾害评估等研究提供重要数据支撑。
地震发生后，尽管采用全球地震台网(GSN)记录的宽频带远震体波数字资料可以快速地给出粗略的矩张量解和最佳双力偶解，但是由于在亚洲地区地震台站稀少以及青藏高原地区地壳的不均匀性使得基于地震波资料给出的地震震中位置有着较大的不确定性。本文采用升、降轨InSAR资料给出了2007年阿里地震的精确震中位置，发现震中位于USGS和NEIC基于地震波资料给出的震中位置的西北5 km处，这个结果类似于2008年于田Mw 7.1级地震的结果[1]。
为了分析不同类型观测数据对反演结果的影响，本文还分别单独采用升轨和降轨数据来进行同震滑动分布反演，这些不同数据反演得到的残差均方误差以及滑动量见表 3。从表 3中可以看到，单独采用升轨数据进行的反演可以很好地拟合升轨自身的观测值，却不能较好地拟合降轨观测值，而基于降轨数据的反演结果也是如此。但是联合升、降轨数据反演给出的同震滑动分布模型虽然对于单类数据的拟合效果不如采用单类数据的反演结果，但是其总体拟合效果要更好，其残差均方误差比采用单类数据的结果降低了0.2～0.3 mm(6%～10%)。此外，不同数据反演给出的最大滑动量和平均滑动量大小基本一致。
| 数据 | 残差均方误差/mm | 滑动量/m | ||||
| T155A | T205D | T155A+ T205D | 最大值 | 平均值 | ||
| T155A | 2.6 | 3.5 | 3.1 | 0.31 | 0.13 | |
| T205D | 3.7 | 2.6 | 3.2 | 0.32 | 0.13 | |
| T155A+T205D | 2.7 | 3.0 | 2.9 | 0.31 | 0.14 | |
InSAR地表形变观测表明2007年阿里地震形成了40 km×20 km大小的形变区，引起约4 cm的下沉。反演结果表明发震断层是美马错—德克玛洛断裂的次生断层，走向158°，倾角43°，倾向西南向，以正断层为主，兼有少量右旋走滑分量。分布式滑动分布模型显示同震滑动分布主要集中在7～12 km深度位置，而在断层的上部没有显著的滑动分布，存在着明显的类似于走滑型地震的“浅部滑动缺失”现象[6]，这可能意味着该发震断层的上部区域在将来有着较大的地震可能性。
致谢：感谢中国科技部和欧洲航天局(ESA)联合资助项目(龙计划三期)提供的Envisat卫星SAR数据(ID: 10607)。
| [1] | ELLIOTT J R, WALTERS R J, ENGLANDP C, et al. Extension on the Tibetan Plateau: Recent Normal Faulting Measured by InSAR and Body Wave Seismology[J]. Geophysics Journal International, 2010, 183(2): 503-535. |
| [2] | HAN Tonglin. Active Tectonic of Xizang[M]. Beijing: The Geological Publishing House, 1987. (韩同林. 西藏活动构造[M]. 北京: 地质出版社, 1987.) |
| [3] | WEN Yangmao, XU Caijun. Ms 7.9 Mani Earthquake Slip Distribution Inversion by a Sensitivity-based Iterative Fitting Method[J]. Geomatics and Information Science of Wuhan University, 2009, 34(6): 732-735. (温扬茂, 许才军. 基于敏感度的迭代拟合法反演玛尼Ms 7.9级地震滑动分布[J]. 武汉大学学报: 信息科学版, 2009, 34(6): 732-735.) |
| [4] | WEN Yangmao, XU Caijun, LIU Yang, et al. Source Parameters of 2008 Qinghai Dachaidan Mw 6.3 Earthquake from InSAR Inversion and Automated Fault Discretization Method[J]. Geomatics and Information Science of Wuhan University, 2012, 37(4): 458-462. (温扬茂, 许才军, 刘洋, 等. 利用断层自动剖分技术的2008年青海大柴旦Mw 6.3级地震InSAR反演研究[J]. 武汉大学学报: 信息科学版, 2012, 37(4): 458-462.) |
| [5] | WANG Leyang, XU Caijun, WEN Yangmao. Fault Parameters of 2008 Qinghai Dacaidan Mw 6.3 Earthquake from STLN Inversion and InSAR Data[J]. Acta Geodaetica et Cartographica Sinica, 2013, 42(2): 168-176. (王乐洋, 许才军, 温扬茂. 利用STLN和InSAR数据反演2008年青海大柴旦Mw 6.3级地震断层参数[J]. 测绘学报, 2013, 42(2): 168-176.) |
| [6] | 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(7040): 295-299. |
| [7] | WERNER C, WEGMULLER, STROZZI T, et al. GAMMA SAR and Interferometric Processing Software[C]//Proceedings of the ERS-ENVISAT Symposium. Gothenburg:[s.n.], 2000. |
| [8] | MASSONNET D, ROSSI M, CARMONA C, et al. The Displacement Field of the Landers Earthquake Mapped by Radar Interferometry[J]. Nature, 364(6433): 138-142. |
| [9] | FARRT G, ROSENP A, CARO E, et al. The Shuttle Radar Topography Mission[J]. Reviews of Geophysics, 2007, 45(2): 20-24. |
| [10] | WEN Y, XU C, LIU Y, et al. Coseismic Slip in the 2010 Yushu Earthquake (China), Constrained by Wide-swath and Strip-map InSAR[J]. Natural Hazards and Earth System Sciences, 2013, 13: 35-44. |
| [11] | HANSSEN R F. Radar Interferometry: Data Interpretation and Error Analysis[M]. New York: Kluwer Academic Publishers, 2001. |
| [12] | LOHMAN R B, SIMONS M. Some Thoughts on the Use of InSAR Data to Constrain Models of Surface Deformation: Noise Structure and Data Down Sampling[J]. Geochemistry, Geophysics, Geosystems, 2005, 6(1), doi: 10.1029/2004GC000841. |
| [13] | WRIGHT T J, LU Z, WICKS C. Source Model for the Mw 6.7, 23 October 2002, Nenana Mountain Earthquake (Alaska) from InSAR[J]. Geophysical Research Letters, 2003, 30(18). |
| [14] | WEN Yangmao, HE Ping, XU Caijun, et al. Source Parameters of the 2009 L'Aquila Earthquake, Italy from ENVISAT and ALOS Satellite SAR Images[J]. Chinese Journal of Geophysics, 2012, 55(1): 53-65. (温扬茂, 何平, 许才军, 等. 联合ENVISAT和ALOS卫星影像确定 L'Aquila地震震源机制[J]. 地球物理学报, 2012, 55(1): 53-65.) |
| [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. |
| [16] | BVRGMANN R, AYHAN M E, FIELDING E J, et al. Deformation during the 12 November 1999 Düzce, Turkey, Earthquake, from GPS and InSAR Data[J]. Bulletin of the Seismological Society of America, 2002, 92(1): 161-171. |
| [17] | FENG W P, LI Z H, ELLIOTT J R, et al. The 2011 Mw 6.8 Burma Earthquake: Fault Constraints Provided by Multiple SAR Techniques[J]. Geophysics Journal International, 2013, 195(1): 650-660. |
| [18] | FENG Wanpeng, LI Zhenhong. A Novel Hybrid PSO/Simplex Algorithm for Determining Earthquake Source Parameters Using InSAR Data[J]. Progress in Geophysics, 2010, 25(4): 1189-1196. (冯万鹏, 李振洪. InSAR资料约束下震源参数的PSO混合算法反演策略[J]. 地球物理学进展, 2010, 25(4): 1189-1196.) |
| [19] | JONSSON S, ZEBKER H, SEGALL P, et al. 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): 1377-1389. |
| [20] | PARSONS B, WRIGHT T, ROWE P, et al. The 1994 Sefidabeh (Eastern Iran) Earthquakes Revisited: New Evidence from Satellite Radar Interferometry and Carbonate Dating about the Growth of an Active Fold above a Blind Thrust Fault[J]. Geophysics Journal International, 2006, 164(1): 202-217. |