汶川地震同震形变场的GPS和InSAR邻轨平滑校正与断层滑移精化反演

杨莹辉，陈强＊，刘国祥，程海琴，2，刘丽瑶，胡植庆

1西南交通大学遥感信息工程系，成都 610031

2华东交通大学土木建筑学院，南昌 330013

3台湾大学地质科学系，台湾台北

1 引言

2008年5月12日，以汶川县映秀镇为震中发生了MS8.0级特大地震，地震发生机制为龙门山中央断裂的右旋走滑逆冲，其动力诱因为印度板块持续向北运动挤压欧亚大陆板块，使得青藏高原东缘地块向东运动，而与之相邻的四川盆地地块较为坚固（徐锡伟等，2005；单新建等，2009），导致龙门山断层应力持续积累并趋于临界状态，最终积蓄能量从映秀—北川主断裂带释放导致汶川地震的发生.

汶川地震引发龙门山地区强烈的地表形变，同震形变场的准确获取对于分析龙门山地区地表破裂的空间分布特征、震源参数计算以及地球物理反演具有重要的科学意义.汶川地震后，中国地壳运动观测网络项目组及时对龙门山地区的GPS观测站进行了复测，并公布了观测结果（张培震，2008）.GPS形变观测显示，以映秀—北川断裂带为中心，两侧断层上下盘存在显著的相向运动和水平缩短现象，同时断层下盘以下降为主，上盘在地表破裂线周边呈现上升趋势，但在远离破裂线后很快又转为下降运动，GPS形变的这些特征揭示了断层运动以逆冲为主的特性（张培震，2008）.但由于GPS稀疏分布的特性，难以获得发震区完整的高分辨率地表形变场，不利于断层面破裂的精细反演计算.汶川地震前后，ALOS卫星搭载的PALSAR传感器获取了该区域大范围SAR影像数据，InSAR技术可为龙门山地区提供大范围、高分辨率的地表形变观测资料（单新建等，2009）.

汶川地震InSAR同震形变场的观测研究表明，由于雷达卫星数据受轨道误差、大气信号延迟等负面因素的影响，从InSAR干涉图中提取精细化的地表形变场干扰较大，尤其是汶川地震波及范围广，邻轨PALSAR干涉图之间存在严重的形变场不连续问题（Zhang etal.，2008；Hashimoto etal.，2009；Shen etal.，2009；Furuya etal.，2010；Liu etal.，2010；Tong etal.，2010），积极探索邻轨InSAR同震位移的修正方法，有助于准确获取龙门山地区精细的同震形变场与断层滑移参数的精化反演.

为解决InSAR轨道误差残留相位与邻轨形变条纹的不连续问题，本文提出联合GPS形变观测与邻轨相位修正的数据处理策略，利用稀疏分布的GPS位移观测去除PALSAR轨道误差引入的残留平地相位，采用邻轨约束平滑消除干涉图之间的不连续问题，从而获得龙门山地区大范围精细的同震形变场，在此基础上，采用弹性半空间位错模型反演地壳发震断层的滑移模型，并由此计算汶川强震的地震矩数值.

2 ALOS／PALSAR影像干涉处理

使用ALOS卫星PALSAR传感器获取的地震前后雷达影像开展干涉处理，以其中连续6个升轨条带共46幅影像作为同震形变观测数据，各轨道影像的基本信息如表1所示.

雷达影像各条带的场景覆盖范围如图1红色矩形方框所示，图中的绿线为野外调查获得的汶川地震地表破裂线，主要包括中央主断裂（映秀—北川）、小鱼洞断裂和前山断裂（灌县—彭县），小鱼洞断裂带位于中央断裂带和前山断裂带之间，前山断裂带长60余公里，在地下约10km处并入中央主断裂带（李勇等，2009）.图1中的蓝色与棕色箭头表示龙门山地区的GPS同震形变观测向量（张培震，2008），其中棕色箭头代表单位长度形变为40cm，蓝色箭头单位长度形变为10cm.图中的红色小圆为MS＞4.0级的余震震中位置（Zhang etal.，2008；Hashimoto etal.，2009；Shen etal.，2009；Furuya etal.，2010；Liu etal.，2010；Tong etal.，2010）.

表1 汶川地震PALSAR影像的干涉对信息Table 1 PALSAR interferometric pairs used in the study of Wenchuan earthquake

对上述PALSAR影像进行干涉处理，使用ASTER DEM作为外部DEM以去除地形相位分量，获得雷达差分干涉相位图（如图2）.从图2可以看出，相邻轨道间的干涉条纹存在显著的不连续性以及轨道残留相位的干扰，在475—476轨道邻接边界（左上角的黑色方框）存在明显的干涉条纹不连续，而图幅中间黑色方框所标识的473条带则存在显著的轨道残留相位，导致远场存在明显的干涉条纹.此外，在右下角的471—472条带连接处，则存在严重的大气相位，其形变条纹与473—474存在显著的不一致性，造成形变场的不连续.为克服上述干涉条纹形变场的不连续问题，本文提出采用GPS形变位移对PALSAR轨道误差残留平地相位进行建模修正，利用邻轨平滑约束条件校正邻轨干涉条纹.

3 轨道残留相位建模与校正

在InSAR数据处理过程中，平地相位（也称为参考相位）的计算是基于SAR轨道状态矢量（包括位置矢量与速度矢量），而ALOS卫星轨道数据的不确定性导致干涉图中存在残留平地相位，在图2中473条带的干涉相位表现尤为显著，在473条带的黑色矩形框内，地震远场呈现出显著的周期性条纹，从形变场空间分布特征表明，该信号为残留的平地相位.

图1 汶川地震PALSAR影像场景覆盖范围Fig.1 Ground coverage of PALSAR images for the Wenchuan earthquake

图2 InSAR同震形变干涉相位Fig.2 InSAR derived interferogram of coseismic deformation field

图3 轨道残留相位校正前后的干涉相位图（a）校正前的干涉相位；（b）校正后的干涉相位.Fig.3 Local interferogram before（a）and after（b）correction from residual orbit phase modeling

为消除轨道数据误差引入的残留参考相位，本文提出以GPS精密形变测量值作为约束，对同一区域的InSAR与GPS形变观测进行差分建模，从而消除或减弱由于卫星轨道误差引入的残留相位值.

首先将GPS观测得到的三维形变转换至雷达视线方向上的LOS形变（陈强等，2012a），定量的数学转换模型如下：

式中，gLos代表GPS三维形变投影到卫星LOS视线方向的形变，dUp、dN、dE分别为GPS测量得到的竖直向、南北向和东西向形变量，θ代表卫星侧视角，α是卫星飞行方向与北方向的夹角，θ和α可依据卫星参数予以计算.

从图2中未校正的473条带可以看出，轨道残留相位在空间尺度上呈现出显著的线性分布（张磊等，2007），因此构建如下的轨道残留误差修正函数模型：

式中，gLos为GPS形变转换到LOS方向的形变，可通过式（1）计算得到，dLos为InSAR观测得到的形变，l、p分别表示GPS测站在雷达干涉图中的方位向（Azimuth）和距离向（Range）坐标，a0、a1、a2为待求模型系数，采用3个或3个以上的GPS形变观测量，即可求解出上述系数，从而建立轨道残留相位修正模型.

从图2中的三维GPS形变测站分布来看，PALSAR条带的472、473、474、475四个条带中，GPS点数量均大于3个，可用于轨道残留相位的建模，而472条带中的H010号GPS点位于余震分布密集的区域，其InSAR相干性较差，形变观测数据质量不高，不用于轨道误差的校正.

采用上述模型和GPS观测数据分别对473、474和475条带中的残留相位予以建模，获得校正后的差分干涉图，如图3所示（对应于图2中的矩形框区域），对比校正前后干涉相位，可以看出由于轨道不确定性所导致的残留相位在校正后相位图中已得到较好的消除，在473—474条带的边界已经没有明显的多余条纹.

4 邻轨平滑校正

对于PALSAR条带内无足够数量GPS校正点的情况下，如475—476的条带连接处，本文提出了基于邻轨间干涉形变差异建模的相对校正方法.由于主震震级造成的地表形变远大于余震和其他因素引起的形变（如震后余滑和孔隙回弹效应在短时间内造成的形变量较小，与主震形变量级相比不显著），因此可以认为在上述PALSAR影像的时间基线跨度内，相邻条带重叠数据具有平滑一致的形变量.以GPS校正后的干涉图为基准，对邻轨干涉图与基准干涉图之间的形变差异进行建模，从两幅邻轨重叠影像中提取高相关性的目标点，解算形变差异模型系数，进而实现对未校正PALSAR条带形变场与干涉相位的修正.该方法充分利用邻轨重叠区域已校正的干涉相位数据，对无足够可用数量GPS形变的干涉条带进行校正，可进一步减弱未校正条带内的轨道相位残差和大气残留相位.

图4为邻轨条带的振幅影像，蓝色虚线方框为相邻条带间的重叠部分，顾及条带重叠区的同震形变应保持一致性的原则，以其中一幅经GPS校正后的干涉数据作为基准，对另一幅未校正干涉条带相位值进行建模修正.例如，475条带已经过GPS校正，且与474和473形变条纹保持较好的一致性，则以475为基准对邻轨的476干涉条纹进行建模修正，具体算法如下.

首先在邻轨PALSAR影像的重叠区进行同名点识别，基于相关系数最大值匹配方法开展重叠区同名点搜索，结果如图4所示，图中的蓝色和白色点代表所识别的同名点.为保证后续建模的可靠性，提取相关系数大于0.4的同名点（图中白色点）用于数学建模，计算高相关性同名点在两个条带中的形变差，依据其在未校正条带中的方位向和距离向坐标，构建如下式的形变校正模型：式中，Δd表示两条带同名点形变量之差，l、p分别为同名点在被校正条带中的方位向和距离向坐标，bi（i＝0，1，2，3，4，5）为待求模型参数.当高质量的同名点数量大于6，使用最小二乘方法解算出模型系数，然后使用该模型对邻轨条纹相位进行数据校正，获得校正后的差分干涉图，图5为条带475与476经邻轨校正后的干涉相位图，从图5中可以看出，条带之间的不连续性问题已得到较好解决，其邻轨干涉条纹连续性较好.

对上述校正后的差分干涉图进行相邻条带相位数据拼接，获得完整的地表形变干涉相位图，如图6所示，使用Snaphu软件对拼接后的干涉相位进行全局相位解缠和绝对形变解算（Chen and Zebker，2001；Chen and Zebker，2002；陈强等，2012b），获得地表沿卫星雷达视线方向（LOS）的形变场，顾及图2中的JYAN GPS基站所在的InSAR像元具有较高的相干性，且远离断层破裂线，选取该GPS观测点作为相位解缠的形变基准点，最终获得PALSAR影像区域完整的同震形变场，如图7所示.

图4 邻轨PALSAR振幅影像重叠区与联系点Fig.4 Spatial overlapping of adjacent PALSAR amplitude maps and detected tie points

图5 邻轨约束校正前后的干涉相位（a）校正前的干涉图；（b）校正后的干涉图.Fig.5 Local interferogram before（a）and after（b）correction using adjacent geometric constraints

图6 —7所示的校正后InSAR同震形变场与已有研究成果相比（单新建等，2009；Shen etal.，2009；Furuya etal.，2010；Liu etal.，2010），整体趋势基本一致，例如，观测到最大的地表形变场位于北川区域，LOS向形变量级在1.2m左右，汶川—映秀段断层上盘存在量级0.6～0.8m的显著形变区，成都平原存在厘米量级的地表形变（Hashimoto etal.，2009；Liu etal.，2010；Shen etal.，2009；Furuya etal.，2010；Tong etal.，2010）.

然而，校正后的InSAR形变场在局部细节区域呈现出显著的不同，首先，已有研究中条带473存在的严重轨道误差在本文中通过GPS校正方法得以很好去除，并在此基础上，使用473条带的观测数据对条带472和471进行了邻轨平滑校正，显著减弱了条带471—472中存在的大气、轨道等系统误差.此外，依据邻轨校正方法，条带476也基于475进行了误差校正.校正后的干涉相位较好地解决了已有研究成果中显著存在的邻轨拼接不连续现象，使得最终提取的形变场具有较好的邻轨平滑性，有利于后续断层滑动分布模型的精确反演.

图6 校正后的InSAR同震形变干涉条纹Fig.6 Corrected InSAR interferogram related to coseismic deformation

图7 校正后的汶川地震InSAR同震形变场Fig.7 Corrected InSAR surface deformation field of Wenchuan earthquake

5 同震形变场精度评定

5.1 绝对精度分析

为定量分析上述校正后InSAR同震形变场的精度与可靠性，以高精度GPS观测值作为参考，对InSAR形变与同点位GPS形变进行差异比较，统计分析校正后InSAR同震形变场的实际绝对精度.

将PALSAR干涉图中的GPS点位分为两类，干涉相干性在0.3以上的点位划分为高相干点，低于该相干阈值则为低相干点，将GPS观测得到形变位移投影转换到雷达LOS方向上的位移，定量数据比较如表2所示.

从表2中可以看出，经轨道残留相位校正和邻轨平滑校正后的同震形变场，高相干点的形变量残差由8.95cm减小至3.59cm，精度提高了60%，低相干点平均残差由25.21cm减小至15.00cm，精度提高达40%，表明本文提出的轨道残留相位建模和邻轨平滑方法具有较高的可靠性，同震形变场的精度得以显著提高.

5.2 邻轨形变的相对精度分析

对于地震断层错动引发的地表形变，在断层上下盘错动的边界处通常存在显著的形变突变，如逆断层上盘呈抬升运动，而下盘则表现为下降，对于其余方位的形变场则较为平滑.因此，在InSAR同一个条带内的同震形变场，除了断层上下盘分界处的地表破裂带区域形变不连续之外，其余地方则呈现出较好的平滑度和连续性.然而，在相邻两个不同PALSAR条带之间的边界区域，则可能出现前述的由于轨道误差、大气相位等因素引起的干涉相位或形变场不连续现象，因此采用条带边界形变的平滑程度指标，可对同震形变场的可靠性与平滑度进行相对精度评价.

本文引入拉普拉斯平滑参数作为同震形变场连续性的定量评价指标，以中心像元S（i，j）为计算单元，通过邻近像元的差分计算，获得干涉形变场的拉普拉斯二阶导数，如下式所示：（4）式中，i，j为待计算像素位置，S（i，j）为像元位置处的形变量，Δl1和 Δl2为SAR分辨单元尺寸，ρi，j为平滑因子，｜ρi，j｜的数值越小，表示形变场的平滑度越高（Jonsson etal.，2002）.

表2 修正前后的InSAR形变与GPS形变的差异比较Table 2 Displacement comparison between InSAR and GPS derived deformation field

对校正前后InSAR形变场各相邻条带边界处的相位，使用上述公式（4）计算其平滑度指标，得到如图8所示的邻轨边界平滑指数分布.

从图8可以看出，经轨道误差和邻轨约束校正后，条带边界干涉形变场的非零｜ρi，j｜像元数量显著减少，边界平滑度得到明显改善，对所有邻轨边界像元平滑因子ρi，j进行数据统计，直方图如图9所示，统计数据表明，校正后形变场平滑因子多趋于零，标准偏差减小了33%.对平滑因子进行分区间统计，在平滑因子绝对值较小的（-40，40）区间内，校正后形变场平滑因子比校正前提高了7.7%，对于平滑因子绝对值较大的（-200，-40）、（40，200）区间，校正后形变场的平滑因子减少了6.8%.综合以上分析，经过轨道误差校正和邻轨约束校正后的形变场，邻轨边界具有更好的平滑度，形变条纹的不连续性得以较好修正.

图8 同震形变场校正前后的邻轨边界平滑度（a）校正前的边界平滑度；（b）校正后的边界平滑度.Fig.8 The smoothness of coseismic interferogram respectively before（a）and after（b）correction

图9 形变场的邻轨边界平滑因子统计（a）校正前邻轨边界平滑因子；（b）校正后邻轨边界平滑因子.Fig.9 Statistics of adjacent track smoothing factor of derived deformation field

6 断层滑移模型反演

在精确获得同震地表形变场基础上，采用Okada弹性半空间位错模型反演发震断层表面的滑移分布（Okada，1985）.通过对野外地质调查数据分析和综合已有研究成果（Zhang etal.，2008；Hashimoto etal.，2009；Liu etal.，2009；Shen etal.，2009；Furuya etal.，2010；Liu etal.，2010；Tong etal.，2010），本文构建了如图10中所示的F1—F6六段断层分段模型，首先假设断层表面为均匀滑动，使用非线性遗传算法对断层参数进行优化反演，为保证反演结果的可靠性，我们赋予反演模型不同的初始值，进行多次反演搜优计算，并对反演参数集进行统计分析，最终获得了各子断层的几何参数及其对应精度，如表3所示.

表3 汶川地震断层几何Table 3 Fault geometry of Wenchuan earthquake

对于中央断裂带（映秀—北川）长大断层来说，断层面的实际滑动往往呈现为非均匀分布，为更精细获得局部断层滑动的空间特征，对表3中各子断层进行4km×4km的格网离散化处理，共获得1051个子断层，以其作为基本单元精细反演断层滑动的空间分布，采用拉普拉斯平滑约束断层滑动模型，如式（5）所示：

式中，Gi为根据离散子断层单位走滑量和逆冲量计算得到的InSAR数据反演格林函数，Hi为GPS数据反演格林函数，α、β分别对应其权重因子，Los、def为InSAR和GPS观测得到的地表形变量，Li为拉普拉斯约束格林函数，κ代表平滑因子分别代表断层的走滑量和逆冲量.

由于GPS形变观测在垂直方向的精度较弱，且垂直向观测资料较少，因此在反演过程中仅使用了GPS平面观测数据.InSAR形变观测数据的条带471和472在校正后仍可能存在电离层及大气负面影响，为避免其对反演结果造成不利影响，471—472条带中形变场连续性较差的数据不参与反演计算.此外，为尽可能减弱不可靠的形变观测结果对反演的影响，使用InSAR相干性大于0.3的形变观测数据，并对其进行四叉树降采样工作，该采样方法在高形变梯度区域采样量大，在低形变梯度区域采样数据量小，能最大限度地保留形变场的空间特征.降采样过程一方面可避免地震远场InSAR非构造形变噪声对反演结果的影响，另一方面也可有效避免因观测数据过多引起的计算负荷过重和计算效率低下的问题，最终的降采样形变场如图10所示.

图10 四叉树降采样InSAR形变场Fig.10 Distribution of InSAR sampled points based on quad-tree method

对于InSAR和GPS形变观测数据用于断层滑移反演的权重分配，本文将依据其观测精度进行确定.根据前述的形变场精度评定可知，高质量InSAR观测数据的残差为0.036m，在实际反演中使用相干性较高的观测点，此处适当放宽其精度取为0.04m，对于GPS形变观测数据的精度，根据相关研究可知其平面精度约为2～5mm，取3mm作为精度指标（张培震，2008）.因此，InSAR和GPS点的误差之比为13.3∶1.0，反比例计算可得其相应权重分别为α＝0.075、β＝1.0.

在确定了数据权重之后，赋予反演模型连续变化的平滑因子κ，并依据反演结果计算获得相应断层面粗糙度与模型残差的对应关系，如图11所示，从图中可看出，随着粗糙度的增大，模型残差急剧减小，顾及模型残差和断层粗糙度的折中关系，我们选取了如图11中箭头标明处的粗糙度为7.5cm／km处的平滑因子为最优的κ取值.

图11 模型残差与断层面粗糙度关系曲线Fig.11 The trade-off curve with misfit function of InSAR and GPS plotted as function of average solution roughness

图12 断层反演滑动及其精度分布（a）断层滑动三维分布；（b）断层滑动误差分布；（c）断层滑动平面分布.Fig.12 Fault slip distribution and its error（a）Three-dimensional distribution of coseismic slip；（b）Errors of coseismic slip；（c）Coseismic slip distribution on fault surface.

确定了各项模型参数后，在非负最小二乘原则下进行断层滑动量的反演计算，获得汶川地震断层滑动的空间分布，如图12所示，其中图12a为反演断层滑动的三维分布.为了验证滑动反演结果的可靠性，我们基于蒙特卡洛算法生成100组带随机扰动的观测数据集，并基于这些数据反演相应的断层滑动分布，对所有计算得到的断层滑动分布集进行统计分析，从而估算出反演模型的精度分布，如图12b所示.统计表明反演模型的平均精度为0.014m，表明图12a的反演结果是可靠的.同时需要指出的是，最大误差分布位于青川地区16km深度处，分析此现象的可能原因，一方面是该区域的InSAR观测数据由于存在大气影响而未用于参与滑动反演，导致反演模型对此部分断层滑动的约束较弱；另一方面是汶川主震后存在密集余震，InSAR观测的滞后性使得相邻条带形变场中包含了不同余震的形变信息，对反演精度造成了一定影响.

为清晰显示断层面的滑动分布状态，我们绘制了如图12c所示的断层滑动平面分布图，从中可看出，断层滑动整体表现为逆冲为主兼具右旋走滑特性，平面滑动量主要分布于映秀、北川和青川区域，断层最大滑动量为9.0m，位于北川地表附近，导致北川地区地表破裂最为严重.在地壳深度方向，地表以下0～16km范围内孕育了较大的滑动成分，而16km以下滑动量分布并不显著.

由于断层滑动反演所使用的地表InSAR形变场约束以及断层几何模型等的不同，本文的反演结果与已有研究成果存在一定差异.首先，断层北川段的滑动量在沿断层走向分布更广，但在深度方向上相对较浅，而映秀—汶川段的平均滑动量级较已有研究成果偏小（Hashimoto etal.，2009；Shen etal.，2009；Furuya etal.，2010；Liu etal.，2010；Tong etal.，2010）；此外，青川段滑动量分布更为深入地下，这可能与青川段震后密集的余震分布相关（朱艾斓等，2008）.

根据上述反演出的断层滑移模型与断层几何参数，正演计算汶川地震的InSAR地表形变场，如图13中的左下插图（a），与图6中校正后InSAR同震形变场相比，两者具有较高的一致性，尤其是在映秀和北川等形变严重的区域，正演数据与InSAR实际观测值吻合的较好.但需要说明的是，在主断层北端的青川区域，由于InSAR观测数据存在严重大气干扰负面影响，形变观测值可靠性相对较弱，未参与到滑动分布模型的反演，使该部分数据与正演形变场存在显著差异.

进一步分析同震形变场的模型反演残差与精度，图13为GPS的反演残差向量（各矢量箭头）与InSAR的反演残差数值分布.图（b）与（c）为GPS观测值与模型计算值的对比图，GPS在北东方向的反演残差分别为4.3cm和6.2cm.此外，统计表明约有75%的GPS平面残差小于2cm，且仅在近断层附近存在残差较大的现象.GPS的整体形变残差为5.5cm，与同类已有研究基本相当.然而，本文InSAR反演形变残差有显著改善，从图13中的InSAR残差分布中可看出，大部分区域残差小于10cm，形变残差统计值达到9.2cm，与已有同类研究成果（平均残差约15cm）对比（Zhang etal.，2008；Hashimoto etal.，2009；Liu etal.，2009；Shen etal.，2009；Furuya etal.，2010；Tong etal.，2010），InSAR形变反演精度提高约30%，验证了InSAR同震形变场的高精度特性与可靠性，同时也表明前述InSAR形变场校正方法的有效性与准确性.

需要指出的是，断层滑动模型的反演残差一方面来源于地表形变场约束的观测误差（包括InSAR形变场中残留的大气、轨道误差等），另一方面，主要受限于不确定的复杂断层几何、位错理论模型的系统误差等.本文断层滑动反演的统计残差为9.2cm，仍大于InSAR高质量点3.6cm的观测校正残差，这主要是由于断层几何参数的精度、以及Okada位错模型未考虑地壳不均匀特性所带来的系统误差.其中，断层几何对反演结果影响最为显著，已有研究均假定断层几何在地壳不同深度具有相同的断层倾角，这与实际断层地下延伸模式可能存在差异.因而，即使在获得高精度地表形变场的情况下，对于龙门山断裂带复杂的断层几何来说，仍可能出现相对较大的反演残差.但相对于已有研究平均15cm的反演残差，本文反演精度（9.2cm）有较大提高，这也正是InSAR形变场校正对于改善滑动模型反演精度的重要贡献，所获得的断层滑动精细分布有助于进一步推进对汶川地震断层破裂机制的认识和深入理解.

依据上述反演所获得的断层几何参数和滑动分布模型，计算得到本次地震的地震矩为8.0469×1020N·m，对应的矩震级Mw7.9，与已有反演结果基本相当（Tong etal.，2010），但本文计算得到的地震矩稍显偏大，一方面是因为本文使用了经邻轨校正后的精细化地表同震形变场，另一方面因为本文所用的断层模型几何倾角与已有方法存在一定差异.

图13 InSAR和GPS形变的联合反演残差（a）模拟InSAR干涉图；（b）GPS北方向形变与模型正演对比；（c）GPS东方向形变与模型正演对比.Fig.13 Residuals of InSAR and GPS joint inversion（a）Simulated InSAR interferogram；（b）Model prediction versus GPS observation in north direction；（c）Model prediction versus GPS observation in east direction.

7 结论

本文采用轨道误差校正方法与邻轨平滑技术对汶川地震形变场进行修正处理，获得汶川地震精细的同震形变场.GPS形变观测与InSAR同震形变的残差分析表明，校正后形变场高质量数据残差为3.6cm，精度提高约60%，低质量数据精度提高约40%；邻轨边界形变场平滑度分析表明，校正后形变场边界平滑因子的标准差减小33%，验证了本文提出的同震形变场校正方法的有效性和可靠性.

基于Okada弹性半空间位错模型的断层滑动反演结果表明，断层滑动量在沿断层走向上主要分布于映秀、北川和青川区域，在深度方向上主要集中在地壳浅部的0～16km范围内，最大滑动量为9.0m，位于北川县城地表附近，InSAR反演后的模型残差达9.2cm，反演精度约有30%的提高，由此计算出的地震矩数值为8.0469×1020N·m.

Chen C W，Zebker H A.2001.Two-dimensional phase unwrapping with use of statistical models for cost functions in nonlinear optimization.Journal of the Optical Society of America A，18（2）：338-351，doi：10.1364／JOSAA.18.000338.

Chen C W，Zebker H A.2002.Phase unwrapping for large SAR interferograms：statistical segmentation and generalized network models.IEEE Transactions on Geoscience and Remote Sensing，40（8）：1709-1719，doi：10.1109／TGRS.2002.802453.

Chen Q，Liu G X，Hu Z Q，etal.2012a.Mapping ground 3-D displacement with GPS and PS-InSAR networking in the Pingtung area，southwestern Taiwan，China.Chinese Journal of Geophysics（in Chinese），55（10）：3248-3258，doi：10.6038／j.issn.0001-5733.2012.10.007.

Chen Q，Yang Y H，Liu G X，etal.2012b.InSAR phase unwrapping using least squares method with integer ambiguity resolution and edge detection.Acta Geodaetica et Cartographica Sinica（in Chinese），41（3）：441-448.

Furuya M，Kobayashi T，Takada Y，etal.2010.Fault source modeling of the 2008Wenchuan Earthquake based on ALOS／PALSAR data.Bulletin of the Seismological Society of America，100（5B）：2750-2766，doi：10.1785／0120090242.

Hashimoto M，Enomoto M，Fukushima Y，etal.2009.Coseismic deformation from the 2008Wenchuan，China，Earthquake derived from ALOS／PALSAR images.Tectonophysics，491（1-4）：59-71，doi：10.1016／j.tecto.2009.08.034.

Jonsson S，Zebker H A，Segall P，etal.2002.Fault slip distribution of the 1999 Mw7.1Hector Mine，California，Earthquake，estimated from satellite radar and GPS measurements.Bulletin of the Seismological Society of America，92（4）：1377-1389，doi：10.1785／0120000922.

Li Y，Huang R Q，Densmore A L，etal.2009.Basic features and research progresses of Wenchuan Ms8.0earthquake.Journal of Sichuan University（in Chinese），41（3）：7-25.

Liu G X，Li J，Xu Z，etal.2010.Surface deformation associated with the 2008 Ms8.0Wenchuan earthquake from ALOS L-band SAR interferometry.International Journal of Applied Earth Observation and Geoinformation，12（6）：496-505，doi：10.1016／j.jag.2010.05.005.

Liu Z J，Zhang Z，Wen L，etal.2009.Co-seismic ruptures of the 12May 2008，Ms8.0Wenchuan earthquake，Sichuan：Eastwest crustal shortening on oblique，parallel thrusts along the eastern edge of Tibet.Earth Planet.Sci.Lett.，286（3-4）：355-370，doi：10.1016／j.epsl.2009.07.017.

Okada Y.1985.Surface deformation due to shear and tensile faults in a half-space.Bull.Seismol.Soc.Am，75（4）：1135-1154.

Shan X J，Qu C Y，Song X G.2009.Coseismic surface deformation caused by the Wenchuan Ms8.0earthquake from InSAR data analysis.Chinese Journal of Geophysics（in Chinese），52（2）：496-504.

Shen Z K，Sun J B，Zhang P Z，etal.2009.Slip maxima at fault junctions and rupturing of barriers during the 2008Wenchuan earthquake.Nature，10：718-724，doi：10.1038／ngeo636.

Tong X P，Sandwell D T，Fialko Y，etal.2010.Coseismic slip model of the 2008Wenchuan earthquake derived from joint inversion of interferometric synthetic aperture radar，GPS，and field data.Journal of Geophysical Research，115（B4）：1-19，doi：10.1029／2009JB006625.

Xu X W，Zhang P Z，Wen X Z，etal.2005.Features of active tectonics and recurrence behaviors of strong earthquakes in the western Sichuan province and its adjacent regions.Seismology and Geology（in Chinese），27（3）：446-461.

Zhang L，Wu J C，Chen Y L，etal.2007.Determination of coseismic displacement field of Taiwan Chi-Chi earthquake based on Doris.Journal of Geodesy and Geodynamics（in Chinese），27（3）：18-24.

Zhang P Z.2008.Co-seismic deformation field of 2008 Ms8.0 Wenchuan earthquake revealed by GPS data.Science in China Press（in Chinese），38（10）：1195-1206.

Zhang Y H，Gong W Y，Zhang J X，etal.2008.Measuring coseismic deformation of the sichuan earthquake by satellite differential INSAR.International Conference on Earth Observation Data Processing and Analysis，7285：4F-1：4F-8，doi：10.1117／12.816004.

Zhu A L，Xu X W，Diao G L，etal.2008.Relocation of the Ms8.0 Wenchuan earthquake sequence in part：preliminary seismotectonic analysis.Seismology and Geology（in Chinese），30（3）：759-767，doi：10.3969／j.issn.0253-4967.2008.03.014.

附中文参考文献

朱艾斓，徐锡伟，刁桂苓等.2008.汶川Ms8.0地震部分余震重新定位及地震构造初步分析.地震地质，30（3）：759-767，doi：10.3969／j.issn.0253-4967.2008.03.014.

张培震.2008.GPS测定的2008年汶川Ms8.0级地震的同震位移场.中国科学，38（10）：1195-1206.

张磊，伍吉仓，陈艳玲等.2007.基于Doris的台湾集集地震同震位移场计算.大地测量与地球动力学，27（3）：18-24.

徐锡伟，张培震，闻学泽等.2005.川西及其邻近地区活动构造基本特征与强震复发模型.地震地质，27（3）：446-461.

单新建，屈春燕，宋小刚等.2009.汶川Ms8.0级地震InSAR同震形变场观测与研究.地球物理学报，52（2）：496-504.

李勇，黄润秋，Densmore A L等.2009.汶川8.0级地震的基本特征及其研究进展.四川大学学报.41（3）：7-25.

陈强，杨莹辉，刘国祥等.2012（b）.基于边界探测的InSAR最小二乘整周相位解缠方法.测绘学报，41（3）：441-448.

陈强，刘国祥，胡植庆等.2012（a）.GPS与PS-InSAR联网监测的台湾屏东地区三维地表形变场.地球物理学报，55（10）：3248-3258，doi：10.6038／j.issn.0001-5733.2012.10.007.