利用PALSAR数据反演2010年玉树地震断层的同震滑动分布

王永哲，朱建军，李志伟，欧自强

中南大学 地球科学与信息物理学院，湖南 长沙410083

1 引 言

2010年4月14日，在中国青海省南部的玉树县发生了MW6.9级强烈地震，震中位于33.2°N、96.6°E，震源深度17km（USGS），此次地震造成了2220人死亡，70多人失踪（据中华人民共和国民政部网站）。玉树县地处巴颜喀拉地块南边界，该次地震被认为是发生在巴颜喀拉地块南边界的甘孜—玉树断裂带上［1］。甘孜—玉树断裂带是一条北西走向的左旋走滑断裂带，长约700km［2］，自晚新生代以来，构造活动剧烈，引发了一系列的大地震［3］。由于印度板块向北推挤、青藏高原东扩、巴颜喀拉块体向东挤出，使得块体边界应力不断积累并达到一定程度，最终以玉树地震的形式得到了释放［4］，对玉树地表造成了严重的破坏。

利用大地测量形变数据来反演地震发震断层的滑动分布，可以为地震研究提供可靠依据。文献［5—6］曾利用GPS形变观测数据，对汶川地震发震断层的滑动分布进行过研究，并得到了比较可靠的结果。但是在玉树地区，GPS和其他传统形变测量数据严重缺乏。利用震前、震后的SAR数据进行干涉处理，充分发挥InSAR技术高空间分辨率、大面积获取地表信息的特点［7］，可以获得雷达视线方向的同震形变场，为地震地区的形变研究提供详细资料。InSAR技术首先成功应用于美国Landers地震同震形变场的监测［8］，结果很好地反映了地震引起的地表形变情况。随着InSAR技术的成熟，该技术已广泛用于地震引起的地表形变的监测［9－12］。文献［1，13］分别利用InSAR技术获取了玉树地震的同震地表形变场。在这些形变干涉处理结果中，存在着比较明显的误差残余相位信息，其中包括轨道误差和大气延迟误差等，形变场并不准确，且未对引起地表形变的断层滑动情况进行研究和分析。

为深入了解玉树地震断层的滑动分布情况，文献［14］最先通过获得的玉树地震InSAR同震形变场，基于弹性位错理论反演了5段不连续断层的滑动分布，但其InSAR同震形变场中仍存在一定量的轨道残差相位，反演的结果在地表的表现比中国地震局实地调查的结果要小。文献［15］利用InSAR和地震波的数据反演了该次地震的同震断层滑动分布，反演的最大滑动量为1.5m，在地表以下4km深处，尽管反演结果与实测InSAR数据拟合较好，但他们将断层分成3段不连续的分段断层，且断层错动引起的地表滑动量与中国地震局实地调查的结果相差较大。文献［16］利用宽幅和条带SAR数据获取了此次地震的形变场数据，基于均匀弹性半空间位错模型解算了断层的滑动分布，他们在反演过程中，将隆宝滩及玉树的两处断层作为连续断层考虑。根据文献［3］的实地调查，隆宝滩和玉树两处地表破裂带并未连续，且呈错开分布。针对观测数据分析和断层分段中存在的问题，本文首先通过最小二乘多项式拟合的方法去除残差相位，得到质量较好的玉树地震同震形变场，然后根据实地调查和InSAR形变场将断层分为隆宝滩和玉树两部分，再将玉树部分分成4段处理，采用最小二乘方法对分段断层的滑动分布进行了反演，并与实地调查的结果对比，最后对结果的不确定性进行了分析。

2 数据处理及分析

地震发生前后，日本宇航局ALOS卫星L波段PALSAR传感器对地震地区进行了观测。玉树地区尽管海拔较高，但起伏不大且植被稀疏，因此比较有利于进行干涉处理。根据空间基线和时间基线选取的原则，本文选择了一对影像对进行干涉处理，选择的图像对很好地覆盖了地震发生的区域，图像对信息如表1所示。

表1 ALOS／PALSAR数据信息Tab.1 Details of ALOS／PALSAR data

干涉处理所使用的软件是瑞士的GAMMA商用软件，利用软件将影像对进行配准、重采样处理，采用二轨差分的方法获取形变干涉图。使用美国国家航空航天局水平分辨率为90m（3″）的SRTM数据去除地形相位。由于空间基线较长，且PALSAR数据没有精密轨道信息，利用垂直基线模拟的地形和平地相位不准确，因此，经过地形和平地相位去除后，干涉图中仍存在着由于轨道误差引起的残余相位趋势。为了去除差分干涉图中的轨道残差相位，本文采用最小二乘方法拟合轨道残差多项式模型，从而削弱轨道残差的影响。形变主要集中在图像中部区域，离震区较远的地区无地表形变发生，其对应的相位为轨道残差相位，也就是说，利用该模型去除轨道残差影响，不会造成形变场信息的畸变。

为了消弱噪声对结果的影响，改进的Goldstein滤波方法［17］用于对干涉图进行滤波处理。干涉图中的相位为缠绕的相位，其并不能完全反映地表真实的形变情况。本文利用最小费用流方法进行解缠处理。玉树地震发生前后时间段内，该地区气候干燥，尽管大气影响对形变场获取影响较小，但本文也使用简单的最小线性回归的方法削弱了大气延迟的影响，该方法主要是利用大气相位与地形之间的简单线性关系来去除大气的相位。最终得到了玉树地震的同震形变场（图1所示为形变干涉图和形变图，黑色曲线为根据形变陡变确定的断层迹线，五角星为USGS的主震和最大余震的震中位置）。

形变干涉图基本上反映了玉树地震的地表形变分布情况，每个条纹代表雷达视线方向0.118m的形变变化。地震发生时，由于地表发生了破裂，破裂带附近的地表形变较大，地物的反射特性发生了较大的变化，从而引起了干涉失相关，因此，在图1中掩膜了这一部分失相干区域。同震形变干涉图中，在USGS震中地区（实地地质调查的隆宝滩地区）（Ⅰ区）以及震中东南方向距离玉树县西北20km左右的地区（Ⅱ区），干涉条纹比较集中，可以认为，地震引起的形变主要集中在这两个区域，其中Ⅰ区形变相对较小，两个条纹变化。Ⅱ区形变相对较大，3～4个条纹变化，可见在玉树县地区的地表形变较大，这也是造成玉树城区受灾比较严重的原因。

形变图很好地反映了同震地表的变化情况，地震引起的地表形变主要集中在70km×70km范围内。在东北部地区，地表为沿雷达视线方向的抬升；在西南部地区，地表为沿雷达视线方向的下降，地震地表形变最大值为－0.442，位于Ⅱ区。东北和西南的形变陡变带清晰可见，在Ⅱ区地表的陡变带与实地的地质调查所得地表破裂带吻合较好，长度大约为50km。而在Ⅰ区，尽管地表基本没有发生破裂，但可以认为该陡变带为断层在地表的错动迹线。

3 分段断层滑动分布反演

地震发生后，文献［18］利用地震台网的地震波形资料，反演得到了玉树地震的矩张量解，认为该次地震的发震断层为玉树断层，走向119°，倾角83°。在确定断层倾角时，本文作者进行了简单的模拟试验，发现当倾角大小改变时，模拟值与观测值残差的改变量并不大。因此，在进行滑动量反演时，采用与地震波形反演结果相同的83°作为断层的倾角。中国地震局地质研究所的科考工作人员对地表破裂情况进行了实地调查，地表破裂发生在为甘孜—玉树断裂的西北段，同时根据实地调查结果及InSAR形变陡变带分析，本文构建了玉树地震的断层几何，主要包括断层走向角、长度和宽度等，断层由5段组成，共长79.80km，断层几何参数如表2所示。

表2 5段断层模型几何参数Tab.2 Geometric parameters of five segmented fault models

3.1 形变场重采样

由于InSAR数据高空间分辨率的特点，将InSAR结果作为反演的观测值，数据量大得惊人，本文处理得到的地面观测点数据达到106，如果将所有的地面观测点都作为约束进行反演，不但浪费时间也没有必要。因此，需对形变场数据进行合理的重采样处理。在选点时，近场数据对断层参数的约束更敏感，而远场数据对参数的约束同样也有贡献作用。考虑到此次地震引起的地表形变梯度较小，影响范围不大，所以本文采用平均选点法进行重采样，采样间距仅为150m。相比适用于形变梯度大的四叉树选点方法［19］，本文采样方法不但保证了近场有足够密度的采样点，也同样保证了远场数据的采样数量，重采样后用作反演的地面观测点总数为2132个。

3.2 反演方法

首先将分段后的断层分成若干个小的断层单元，每个断层单元具有均匀滑动特征。基于均匀介质的弹性半空间位错模型［20］，介质的泊松比设为0.25，在确定了断层的几何参数之后，InSAR观测数据与地震断层滑动量之间的关系可以表示为

式中，dLOS为观测得到的InSAR视线向形变数据；G为格林函数矩阵；s为小断层单元的滑动分量；ε为观测误差。

在反演过程中，为了避免结果的不稳定，加入了平滑约束，使滑动量的二维二次导数最小，此时构建的方程［19］为

式中，D为拉普拉斯平滑算子；k2为平滑因子，根据函数的残差和粗糙度关系来确定。

在式（2）中，InSAR形变数据与每个小断层单元的滑动量之间呈线性关系，使用最小二乘方法，来反演分段断层的同震滑动分布。

3.3 滑动分布反演

根据确定的断层几何参数，将断层沿走向方向均匀划分成37列，深度方向均匀划分为10行，断层单元总数为370，根据断层单元的几何参数，生成格林函数矩阵G。反演过程中由于失配值和粗糙度之间存在一定的折衷对比关系（关系曲线如图2所示），根据折衷曲线，本文选择的光滑因子为2.3。

图1 同震形变干涉图（a）与同震形变场（b）Fig.1 Coseismic interferogram（a）and coseismic deformation field（b）

图2 失配值和粗糙度之间的折衷关系曲线，选择点是光滑因子为2.3的位置Fig.2 The curve of trade－off between the misfit values and roughness，the selected point is the position of smoothing factor 2.3

经过反演计算，得到了断层滑动的分布，模拟的形变值与反演所用观测值拟合较好，均方根误差为0.9cm，最大的观测值拟合残差为6.3cm，拟合残差直方图如图3。图4为反演结果的三维显示，断层以左旋走滑为主，在断层滑动比较大的地方，存在一定的逆倾滑动，该结果与甘孜—玉树断裂带的性质比较一致。断层的滑动出现了两个峰值区：第1个滑动峰值区发生在第1段断层的沿断层面以下14km深处，最大值为2.084m；第2个峰值区发生在第4段断层的沿断层面以下6km范围深处，最大值为地表处的1.774m。从断层滑动分布图中看出，地表破裂应该主要发生在第1段断层处，尽管断层滑动量比较大，但在地表显现出来的错动位移很小，这与地质调查［3］的隆宝滩破裂位移结果比较吻合。第4段断层处，断层的滑动有两个区域，上部滑动直接出露到地表，并在地表形成了破裂，将出露地表的滑动结果与地质调查的破裂位移结果［4］进行了对比，比较结果列于表3中。本文的结果与地质调查的结果相当，但有一定的差别，该差别主要是由断层位置和反演断层单元大小划分的影响造成的：一方面由于在反演前，将断层在地表划分为5段，在细节上点位位置与实地调查的位置有一定距离，造成数据有偏差；另一方面，在反演时，第4段断层上每个单元沿走向方向上为2.5km大小，所得滑动值为整个单元上滑动量的值，与具体的对应点的实测结果有偏差。正是由于以上两个原因，使得本文结果与地质调查的结果有偏差，且在3、4号点处偏差较大。下部区域有一定的右旋逆冲的滑动量存在。在第2段断层处，地表以下也发生了少量的断层滑动，对地表的影响非常小，这也符合干涉形变图的实测结果。

图3 拟合残差直方图Fig.3 Histogram of the fitting residuals

图4 断层滑动分布图Fig.4 Fault slip distribution

表3 断层错动量与地表实地破裂比较Tab.3 Surface slip comparison between field survey and inversed result

断层平均滑动量为0.669m，从文献［21］中得到的川滇地区的岩石剪切模量μ＝32Gpa，根据滑动量计算的地震矩张量为3.4×1019Nm，对应的矩震级为MW7.0，与利用地震波形反演的结果相当。

3.4 结果分析及讨论

为了估计反演滑动量的不确定性，首先对InSAR形变观测值的误差进行分析。尽管干涉处理中对各种误差进行了削弱和去除，但这些误差的影响在一定程度上仍然存在，尤其以大气延迟误差影响最大。为简化分析，假设整个噪声场中都是大气误差影响，服从各向同性特征［22－24］。因此，各点之间的相关关系只与点之间的距离有关，可以使用经验的地统计学方法来确定观测点之间的协方差，引入协方差函数［25］

式中，C（h）为两点间的协方差；h为两点间距离；a、b为待定系数。

整个形变干涉图中，任意两点之间的协方差都满足以上的关系式，在对观测点进行统计求协方差函数时，忽略震区的观测点，只选择离震区较远地区的数据来拟合函数，认为这些远区数据中不含形变信息而只含有误差信息，使用非线性最小二乘方法拟合函数中的待定系数。求得的模型参数分别为：a＝36.57mm2，b＝23.41km。利用相干性确定的方差［22］以及式（3）确定的协方差函数，求得用于反演的观测值方差－协方差矩阵。根据此方差－协方差矩阵，采用蒙特卡洛方法模拟100组噪声，将噪声分别加到观测数据中，对观测数据进行扰动［26］，然后用新生成的100组数据分别进行反演来确定结果的不确定性。反演结果的不确定性分布如图5所示，误差平均值为0.162m，误差最大值位于第一段断层的14km深处，可以认为在此周围误差相对较大。由于采集的PALSAR图像数据只覆盖了隆宝滩断层影响的一部分区域，而图像则覆盖了玉树断层影响的大部分区域，使得InSAR形变结果对隆宝滩断层的滑动量反演的约束能力要差一些，对玉树断层的约束能力要强，该差别也反映在了滑动量的误差分布图上。

利用反演得到的断层滑动分布，模拟了整个地震区地表PALSAR视线方向的同震形变场以及形变干涉图，形变干涉图如图6（a）所示。从模拟的形变干涉图中可以看出，条纹数和形变区域与实测PALSAR形变场基本一致。将实测值与模拟值做差，得到形变残差干涉图及形变残差图，如图6（b）和图6（c）。从残差图中可以发现，模型对观测值的拟合较好。从形变残差图中可以看出，不存在系统性的残差，整体形变残差的均方根误差为1.4cm，残差较大的观测点主要集中在断层两侧附近位置。在反演时，将断层在地表的迹线简化为几段直线处理，断层附近的点与断层关系发生较大改变，从而引起这些点的拟合残差较大。在相位残差图中可以发现，地震区的大部分区域残差相位较小，整体拟合较好，相位差最大的地方在隆宝滩的第1段断层和第4段断层处的上盘部分，但都小于一个条纹。由以上分析可知，由本文断层滑动模型模拟的形变和相位与观测值拟合的较好，且结果比较可靠。

图5 断层滑动不确定性分布图Fig.5 Distribution of fault slip uncertainties

4 结 论

玉树地震引起的地表形变，主要集中在隆宝滩和玉树两个区域，断层东北部地区沿雷达视线方向抬升，西南部地区沿雷达视线方向下降，在玉树县附近的形变最大。

断层的滑动以左旋走滑为主，滑动量主要集中在两个区域：第1、第4段断层。第1段和第4段断层处有地表的出露；玉树附近地区的第4段断层处，断层滑动引起了地表的严重破裂，地表错动量达到1.774m，该处为地震的宏观震中位置。反演的结果在第2、3、4、5断层处的误差较小，第1断层处较差。根据反演结果计算的地震矩3.4×1019Nm，震级MW7.0。

本文将PALSAR数据作为观测数据进行反演约束，且干涉处理得到是雷达视线方向的形变，反演约束能力并不是非常强，如果能将不同视线方向的InSAR数据和其他大地测量数据进行联合反演，将可以得到更加符合实际的断层滑动量和分布。

图6 模拟的PALSAR同震形变干涉图（a），相位残差图（b）和形变残差图（c）Fig.6 Predicted coseismic PALSAR interferogram（a），phase residuals（b）and deformation residuals（c）

致 谢：研究所用的PALSAR数据由日本宇航局提供，部分图件使用GMT［27］绘制。

［1］ ZHANG Guifang，QU Chunyan，SHAN Xinjian，et al.The Surface Rupture and Coseismic Deformation Characteristics of theMS7.1Earthquake at Qinghai Yushu in 2010［J］.Chinese Journal of Geophysics，2011，54（1）：121－127.（张桂芳，屈春燕，单新建，等.2010年青海玉树MS7.1级地震地表破裂带和形变特征分析［J］.地球物理学报，2011，54（1）：121－127.）

［2］ ZHANG Peizhen，DENG Qidong，ZHANG Guomin，et al.Active Tectonic Blocks and Strong Earthquake in the Continent of China［J］.Science in China：D，2003，33（sup）：12－20.（张培震，邓起东，张国民，等.中国大陆的强震活动与活动地块［J］.中国科学：D，2003，33（sup）：12－20.）

［3］ MA Yinsheng，ZHANG Yongshuang，HU Daogong，et al.The Surface Ruptures and the Macroscopical Epicenter of YushuMS7.1Earthquake［J］.Journal of Geomechanics，2010， 16（2）：115－128.（马寅生，张永双，胡道功，等.玉树地震地表破裂与宏观震中［J］.地质力学学报，2010，16（2）：115－128.）

［4］ CHEN Lichun，WANG Hu，RAN Yongkang，et al.TheMS7.1Yushu Earthquake Surface Rupturesand Historical Earthquakes［J］.Chinese Science Bulletin，2010，50（13）：1200－1205.（陈立春，王虎，冉勇康，等.玉树MS7.1级地震地表破裂与历史大地震［J］.科学通报，2010（13）：1200－1205.）

［5］ XU Caijun，LIU Yang，WEN Yangmao.MW7.9Wenchuan Earthquake Slip Distribution Inversion from GPS Measurements［J］.Acta Geodaetica et Cartographica Sinica，2009，38（3）：195－201.（许才军，刘洋，温扬茂.利用 GPS资料反演汶川MW7.9级地震滑动分布［J］.测绘学报，2009，38（3）：195－201.）

［6］ LI Zhicai，ZHANG Peng，JIN Shuanggen，et al.Wenchuan Earthquake Deformation Fault Inversion and Analysis Based on GPS Observations［J］.Acta Geodaetica et Cartographica Sinica，2009，38（2）：108－113.（李志才，张鹏，金双根，等.基于GPS观测数据的汶川地震断层形变反演分析［J］.测绘学报，2009，38（2）：108－113.）

［7］ LIAO Mingsheng，LIN Hui.Synthetic Aperture Radar Interferometry－principle and Signal Processing ［M ］.Beijing：Surveying and Mapping Press，2003.（廖明生，林珲.雷达干涉测量：原理与信号处理基础［M］.北京：测绘出版社，2003.）

［8］ MASSONNET D，ROSSI M，CARMONA C，et al.The Displacement Field of the Landers Earthquake Mapped by Radar Interferometry［J］.Nature，1993，364（6433）：138－142.

［9］ XU Caijun，LIU Yang，WEN Yangmao，et al.Coseismic Slip Distribution of the 2008MW7.9Wenchuan Earthquake from Joint Inversion of GPS and Insar Data［J］.Bulletin of the Seismological Society of America，2010，100（5B）：2736－2749.

［10］ SUN Jianbao，XU Xiwei，SHEN Zhengkang，et al.Parameter Inversion of the 1997Mani Earthquake InSAR Coseismic Deformation Field Based on Linear Elastic Dislocation Model－I：Uniform Slip Inversion［J］.Chinese Journal of Geophysics，2007，50（4）：1097－1110.（孙建宝，徐锡伟，沈正康，等.基于线弹性位错模型及干涉雷达同震形变场反演1997年玛尼MW7.5级地震参数－I：均匀 滑 动 反 演 ［J］.地球物理学报，2007，50（4）：1097－1110.）

［11］ FUNNING G J，PARSONS B，WRIGHT T J，et al.Surface Displacements and Source Parameters of the 2003BAM（Iran）Earthquake from Envisat Advanced Synthetic Aperture Radar Imagery［J］.Journal of Geophysical Research－Solid Earth，2005，110（B9）：39－49.

［12］ LIU G X，DING X L，LI Z L，et al.Pre－and Co－seismic Ground Deformations of the 1999Chi－Chi，Taiwan Earthquake，Measured with SAR Interferometry［J］.Computers and Geosciences，2004，30（4）：333－343.

［13］ YAO Xin，ZHANG Yongshuang，YANG Nong，et al.D－InSAR Observation of Earth Surface Deformation in theMS7.1Yushu Earthquake［J］.Journal of Geomechanics，2010，16（2）：129－136.（姚鑫，张永双，杨农，等.玉树地震地表变形InSAR观测及初步分析［J］.地质力学学报，2010，16（2）：129－136.）

［14］ ZHA Xianjie，DAI Zhiyang，GE Linlin，et al.Fault Geometry and Slip Distribution of the 2010Yushu Earthquakes Inferred from InSAR Measurement［J］.Bulletin of the Seismological Society of America，2011，101（4）：1951－1958.

［15］ LI Z H，ELLIOTT J，FENG W P，et al.The 2010MW6.8Yushu（Qinghai，China）Earthquake：Constraints Provided by InSAR and Body Wave Seismology［J］.Journal of Geophysical Research，2011，116（B10）：302－316.

［16］ TOBITA M，NISHIMURA T，KOBAYASHI T，et al.Estimation of Coseismic Deformation and a Fault Model of the 2010Yushu Earthquake Using PALSAR Interferometry Data［J］.Earth and Planetary Science Letters，2011，307（3－4）：430－438.

［17］ LI Zhiwei，DING Xiaoli，HUANG C，et al.Improved Filtering Parameter Determination for the Goldstein Radar Interferogram Filter［J］.ISPRS Journal of Photogrammetry and Remote Sensing，2008，63（6）：621－634.

［18］ LIU Chao，XU Lisheng，CHEN Yuntai.Quick Moment Tensor Solution for 14April 2010Yushu，Qinghai，Earthquake［J］.Acta Seismologica Sinica，2010，32（3）：366－368.（刘超，许力生，陈运泰.2010年4月14日青海玉树地震快速矩张量解［J］.地震学报，2010（3）：366－368.）

［19］ JONSSON S，ZEBKER H，SEGALL P，et al.Fault Slip Distribution of the 1999MW7.1Hector Mine，California，Earthquake，Estimated from Satellite Radar and GPS Measurements［J］.Bulletin of the Seismological Society of America，2002，92（4）：1377－1389.

［20］ OKADA Y.Internal Deformation Due to Shear and Tensile Faults in a Half－space［J］.Bulletin of the Seismological Society of America，1992，82（2）：1018－1040.

［21］ WANG Hui，LIU Jie，SHI Yaolin，et al.Dynamic Simulation of Interactions between Major Earthquakes on the Xianshuihe Fault Zone［J］.Science in China：Series D，2008，38（7）：808－818.（王辉，刘杰，石耀霖，等.鲜水河断裂带强震相互作用的动力学模拟研究［J］.中国科学：D辑，2008，38（7）：808－818.）

［22］ HANSSEN R F.Radar Interferometry：Data Interpretation and Error Analysis［M］.New York：Kluwer Academic Publishers，2001.

［23］ BIGGS J，WRIGHT T，LU Z，et al.Multi－interferogram Method for Measuring Interseismic Deformation：Denali Fault，Alaska［J］.Geophysical Journal International，2007，170（3）：1165－1179.

［24］ 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］.Geophysical Journal International，2006，164（1）：202－217.

［25］ SUDHAUS H，JONSSON S.Source Model for the 1997 Zirkuh Earthquake（MW＝7.2）in Iran Derived from JERS and ERS InSAR Observations［J］.Geophysical Journal International，2011，185（2）：676－692.

［26］ FUNNING G.Source Parameters of Large Shallow Earthquake in the Alpine－Himalayan Belt from InSAR and Waveform Modelling［D］.Oxford：University of Oxford，2005.

［27］ WESSEL P，SMITH W.New，Improved Version of Generic Mapping Tools Released［J］.Eos，Transactions American Geophysical Union，1998，79（47）：579－579.