通过InSAR与GPS数据融合获取汶川地震同震三维形变场

宋小刚　申　星　姜　宇　万剑华　单新建　屈春燕

1）中国地震局地质研究所，地震动力学国家重点实验室，北京　100029

2）中国石油大学（华东），地球科学与技术学院，青岛　266580

通过InSAR与GPS数据融合获取汶川地震同震三维形变场

宋小刚1）申星1，2）*姜宇2）万剑华2）单新建1）屈春燕1）

1）中国地震局地质研究所，地震动力学国家重点实验室，北京100029

2）中国石油大学（华东），地球科学与技术学院，青岛266580

针对InSAR与GPS两类独立的观测数据，利用广义测量平差理论中方差分量估计法，合理分配权重，有效地融合了这两类观测数据，从而估算出了地表在三维方向上的形变场。以四川汶川地震为例，利用InSAR干涉测量结果和一定数量的GPS观测值，通过该方法获取了地震断层两侧高相干区域上三维形变场，清晰地显示了汶川地震的逆冲和右旋走滑分量的位置分布整体特征。结果表明，在EW、SN和UD方向上地表形变量与GPS结果能够很好地吻合，融合结果在三维方向上的均方根误差均不超过5cm，证明了该方法能够得到精度较高的地震同震三维形变场，也揭示了利用方差分量估计法对相互独立数据之间有效融合的可行性。

汶川地震三维形变方差分量估计InSARGPS

0　引言

合成孔径雷达干涉测量（InSAR），因具有空间分辨率高、监测范围广和费用低廉等优点，已经成为监测地表形变的主要方法之一（Massonnet et al．，1993）。然而，在很多的研究中，往往由于数据的限制，只能利用单一轨道SAR数据进行干涉测量，得到的结果是地表各个方向的形变在雷达视线方向上的投影，因而无法获取研究区真实的三维位移场。GPS可以精确测定地表三维方向上的形变，但是GPS站点受环境条件和成本的限制，往往分布稀疏、空间分辨率有限。如何将高空间分辨率的InSAR观测结果与高时间分辨率的GPS观测结果有效地融合，国内外很多学者开展了研究。Ge等（2000）提出了DIDP方法，在时空域内内插InSAR和GPS数据来获取地表三维形变场；Guemundsson等（2002）结合GPS与InSAR数据，应用马尔可夫随机模型和模拟退火法进行干涉图像相位解缠，初步实现了InSAR与GPS数据的融合，得出了冰岛地区高空间分辨率的三维地壳变形场；Samsonov等（2006，2007）在联合InSAR和GPS数据方面提出了解析法，并于2007年将这种方法进行了改进，得到了南加州地区地表三维形变场。这些方法在联合InSAR和GPS数据方面提供了很好的思路，但是在融合2类观测数据时没有对数据进行合理的定权处理。

本文结合已有的GPS观测资料（Wang，2011）和InSAR观测资料（屈春燕等，2008），对四川汶川MS8.0地震展开研究。首先，利用精度较高的GPS数据对InSAR数据进行校正；然后，利用方差分量估计合理定权的方式来融合2种类型的观测数据，获得了汶川地震三维形变场；最后，对融合结果和GPS实测值进行比较和精度验证。结果表明，该方法能够得到精度较高的同震三维形变场，均方根误差均不超过5cm。

1　数据的选取

汶川地震发生在青藏高原东北缘龙门山断裂带上，断层上、下盘之间地形起伏剧烈、植被覆盖茂密。L波段SAR数据波长较长穿透能力较好，获取的SAR图像可以保持较高相干性（孙建宝等，2006），有效地克服了雷达回波中高噪声的影响，提高了干涉相位的信噪比。因此，在处理过程中选取ALOS／PALSAR雷达数据作为干涉数据源。选取汶川震区471～477共7个条带震前、震后共112景不同时段的Level 1．0级数据，采用瑞士GAMMA软件处理平台，按断层上、下盘分别进行干涉处理，再将南北两盘地震干涉形变场进行拼接。通过配准、滤波、相位解缠等一系列处理步骤，最后得到数值化的形变场图（屈春燕等，2009；单新建等，2009）。地震破裂带附近存在严重失相干现象，InSAR无法在该区域获取有效的形变信息，在数据处理中加以去除。汶川地震前中国地壳运动观测网络和国家重点基础研究发展计划项目组在龙门山断裂带两侧分别布设一定数量的GPS观测点，这些点都是进行长期观测的稳定点并用钢筋混凝土进行加固，安装强制对中装置，在地震发生以前已经有多期观测并有一定数据积累（国家重大科学工程中国地壳运动观测网络项目组，2008）。甘肃省、陕西省等地方测绘局也布设了不同等级的GPS点，这些GPS站点为震后同震形变场的获取和研究提供了极为重要的资料。汶川地震后，中国地震局和国家测绘局以及地方测绘局等单位快速实施了震区GPS点的复测，其中Wang（2011）全面地给出了国家连续GPS站点和流动复测的GPS站点在同震三维方向的位移量。为了计算方便，选择研究区域为近断层两侧附近，如图1（30.5°～32.5°E，103°～105°N）。

InSAR获取的LOS形变量是一相对值，常规差分干涉测量受到时间、空间失相干的严重制约和大气延迟以及DEM误差等因素的影响（班保松，2010），所以首先考虑用精度较高的GPS点来校正InSAR观测值。

2　利用GPS数据对InSAR观测值进行校正

InSAR观测值中包含轨道误差、大气延迟误差、失相干、DEM误差等各项误差（Wright et al．，2004），再加上InSAR这个相对观测量本身存在着系统性偏差，所以在与GPS数据融合之前，首先利用精确的GPS形变观测数据对InSAR形变场进行校正。文中使用的InSAR条带数据共5条，利用每个InSAR条带范围内的GPS站点对该条带进行改正。具体改正步骤为：1）利用InSAR成像几何关系将GPS站点在U、E、N 3个方向的观测量投影至LOS向即可得到LOS向的GPS形变值；2）针对每个条带，利用其所覆盖的GPS点，拟合GPS投影值与相同位置处的InSAR观测值间的多项式关系，拟合公式采用：g＝ah2＋bh＋c（其中a，b，c为拟合系数，h为InSAR观测值，g为GPS投影值）；3）利用拟合得到的多项式对每个条带的InSAR观测值进行校正：y＝ax2＋bx＋c（其中a，b，c为拟合系数，x为InSAR观测值，y为InSAR改正后的值）；4）然后将5个矫正后的条带数据进行拼接，得到改正后的InSAR形变场（图3）。

图1　重采样后InSAR视线向形变场和GPS站点分布Fig．1　Down-sampled InSAR LOS deformation field and the distribution of GPS sites．

为了评价改正的精度，抽取76个GPS点用来进行改进前、后结果的比较，部分结果列于表1中。在剔除掉＞2倍中误差的点后，计算得到改正前GPS与InSAR观测值之差的中误差是15.17cm，改正后降为4.72cm，改正效果明显。图2为表1结果的直观显示，可以看出，改正后InSAR与GPS差值整体变小，更接近于0。

3　利用方差分量估计融合InSAR与GPS估算三维形变量

3.1三维形变矢量和LOS向的关系

InSAR测量结果是地表各个方向的形变在雷达视线方向（LOS）上的投影，由单一的InSAR观测值难以计算得出地表在水平和垂直方向的真实形变量，被称为LOS向模糊问题（洪顺英，2010）。如图4所示UU，UE，UN是在（UD向之U，EW向之E，NS向之N）3个方向的分量，θ为卫星入射角，α为卫星轨道方向的方位角。AZ为雷达图像方位向，ALD为距离向。dLOS是卫星视线向形变量，VLOS是测量误差（轨道误差、大气延迟、失相干和DEM误差等）。则雷达视线方向的形变可以用下式表示：

ALOS／PALSAR的雷达卫星是升轨观测，从图像的头文件中获得中心入射角θ为38.7°，卫星轨道方位角α为349.7°，上式可以转换为矩阵的形式：

表1　18个GPS点上InSAR改正前、后比较Table 1　Corrected results for 18 GPS points

图2　InSAR改正前、后与GPS的差值比较Fig．2　Differences between uncorrected，corrected InSAR and GPS observations．

图3　校正后InSAR形变场Fig．3　GPS-corrected InSAR deformation field．

图4　LOS向形变与三维形变分量几何示意图（a），地表位移三维分量水平投影（b）Fig．4　Geometry for LOS displacement and 3-D components（a）and plane projection（b）．

由上式可以看出，地表UD向对LOS向影响最大，EW向其次，影响最小的是NS向，EW向和UD向的影响大小相反，计算最终的结果是绝对值较大的一个，这也是难以单独利用In-SAR来判断地表形变是抬升还是下沉的原因，特别是对于汶川地震这种逆冲兼走滑的运动形式。

3.2方差分量估计的引入

对于多种不同来源数据的融合，合理的权重分配是其中最关键的环节，因为不同的数据集有各自不同的精度评价体系。因此，如何打破各自体系而进行有效的融合，定性的评价一般无法满足需求，科学的方法是能够从理论上给出合理的解释，而测量平差里的赫尔默特方差分量估计能够很好地解决这一问题。

基本思想（崔希璋等，2004）是：对于彼此之间相互独立的2类观测数据集，为了合理地确定2类不同观测量的权，先对2类观测量进行初定权P1、P2（设为单位权），然后进行预平差，利用平差处理后的改正数加权平方和VTPV来估计观测值的方差，再依据观测值的方差重新定权。如此重复，直到2类观测值的权趋于合理（胡俊等，2010）。这种平差方法被称为验后方差分量估计。

对InSAR观测值而言有如下观测方程：

对于GPS观测值而言有如下方程：

设定2类观测值的权阵P1，P2初始值为单位权，相应的构成2组误差方程式为

计算观测值改正数的加权平方和Wθ和系数矩阵S：

那么根据赫尔默特方差分量估计的严密公式（式（8）），可以估算单位权方差＾θ

图5a为EW向形变量，可以看出断层上盘向东移动，最大位移量为201.11cm，断层下盘向西移动，最大位移量为228.14cm，总体显示出右旋的趋势，而在断层的两端出现了不一致性。图5b为SN向形变，可以看出断层上盘向南移动的区域在主震附近，最大位移量为81.07cm，下盘几乎整体向北移动，最大位移量为68.65cm。主震附近上、下两盘相对形变值较大，说明该区域逆冲分量最大，往NE方向逐渐减弱。东北段显示出下盘向南运动，这正好解释了东北段运动的特点：以右旋走滑为主，而逆冲分量微弱。所以从水平形变分量来看，整体显示出右旋走滑兼逆冲的运动特征。图5c为UD向形变，可以看出，发震断层上、下盘靠近断裂带附近出现多处规模不等的局部隆升和下沉，最大隆升量和最大沉降量分别27.81cm和70.65cm。总之，三维方向的形变特征可以清晰地显示出断层局部运动特征，也反映了发震断层运动的复杂性和非均匀性。

4　GPS实测结果与融合结果的比较分析

为了评价数据融合结果的精度，选取了50个GPS点来进行对比检验。如图6，在EW向GPS观测值与融合结果的偏差为－3.79～8.29cm，其均方根误差为3.47cm；SN向偏差为－4.49～1.19cm，其均方根误差为1.31cm；UD向偏差为－3.63～3.97cm，其均方根误差为1.93cm。其中精度评定公式为（n代表观测点数，Δ为GPS实测值和计算值之间的差值）。

从比较结果中可以看出，GPS和InSAR在水平方向和垂直方向都符合得较好，三维位移分量均方根误差均不超过5cm。三维位移分量中EW向与UD向精度相差不大但明显低于SN向，分析汶川地震断层错动方向与实测GPS和InSAR观测数据特征，总结其中的原因可能包含以下3个方面：

（1）汶川地震发生在以NE向的映秀－北川断裂为核心的龙门山断裂带上，而ALOS卫星飞行的轨道方位角是349.7°近SN方向，SN向形变对LOS向贡献很小，所以InSAR观测量主要是EW向和UD向形变的贡献量。在利用方差分量估计融合时，SN向分量以GPS贡献为主，相应的权重较大，因此最后结果的精度会高于EW向和UD向。

（2）InSAR观测值是某一像元对应空间分辨单元内多点的平均值，而GPS实测的结果是某一点的真实值，再加上大气、轨道及地形误差的影响，InSAR观测值精度较低，而在数据融合时，虽然有GPS观测量的存在，但是InSAR对EW向和UD向分量的贡献较大，从而导致这2个方向整体精度偏低。

（3）InSAR测量结果是主震发生前后对获取的SAR影像进行处理，其中数据选取距离主震发生时刻少则几天多则几个月，而此期间对汶川MS8.0特大地震来说，余震和地表弹性恢复更为频繁，而GPS数据获取时间与InSAR获取时间不完全相同，无疑会导致数据精度的不同，从而给最终的融合结果带来不同程度的偏差。

图5　利用方差分量估计得到的汶川地震三维位移场形变融合结果Fig．5　3－D deformations field estimated by VCE（Variance Component Estimation）based on GPS and InSAR．

图6　GPS实测值与融合结果在三维方向上的比较Fig．6　Comparison between GPS and the results estimated from VCE for 3D deformation components．

5　结论

合成孔径雷达监测的是地表LOS向的一维形变，不能反映地表真实的三维形变信息。对于不同来源的数据，由于精度、空间分辨率和时间分辨率的不同，给计算三维形变场带来了不便。本文利用方差分量估计法，通过合理定权的方式融合InSAR和GPS两种不同类型的观测数据，获取了地表三维形变场。以汶川地震为例，证明了该方法的可行性，并且能够得到精度较高的三维形变场，在EW、SN和UD方向上地表形变量与GPS实测结果能够很好地吻合，揭示了利用方差分量估计法对相互独立数据之间有效融合的可行性。

经过以上分析可以得到如下结论：

（1）方差分量估计法不仅可以将InSAR与GPS数据进行融合，而且对于不同轨道、不同卫星和不同入射角的SAR数据也同样适合。在融合计算过程中还可以进一步将水准观测、气象等多种观测数据通过合理分配权重进行三维形变的估算，进一步提高求解精度。

（2）三维形变场能够给出特定方向上更为精细的形变特征，可以清楚地判断断层局部位置运行的方向和性质，且能够反映出各种分量（逆冲和走滑）的大小，这为准确理解断层运动方式和地震发生机制提供了有价值的数据，也为进一步分析地震断层滑动分布反演提供了详实的资料。

班保松，伍吉仓，陈永奇，等．2010．联合GPS和InSAR观测结果计算汶川地震三维地表形变［J］．大地测量与地球动力学，30（4）：25—35．

BAN Bao-song，WU Ji-cang，CHEN Yong-qi，et al．2010．Calculation of 3D terrain deformation of Wenchuanearthquake with GPS and InSAR data［J］．Journal of Geodesy and Geodynamics，30（4）：25—35（in Chinese）．

崔希璋，於宗俦，陶本藻，等．2010．广义测量平差（第二版）［M］．武汉：武汉大学出版社．

CUI Xi-zhang，YU Zong-chou，TAO Ben-zao，et al．2010．Principle of Generalized Surveying（Second Edition）［M］．Wuhan University Press，Wuhan（in Chinese）．

国家重大科学工程“中国地壳运动观测网络”项目组．2008．GPS测定的2008年汶川MS8.0地震同震位移场［J］．中国科学（D辑），38（10）：1195—1206．

Project Group of the Crust Movement Observation Network of China．2008．The coseismic deformation of the Wenchuan MS8.0 earthquake based on GPS observation［J］．Science in China（Ser D），38（10）：1195—1206（in Chinese）．

洪顺英．2010．基于多视线向DInSAR技术的三维同震形变场解算方法研究及应用［D］．［学位论文］．北京：中国地震局地质研究所．

HONG Shun-ying．2010．The resolving method and application of the 3D coseismic deformation field based on the multi-LOS DInSAR technology［D］．Ph D Thesis．Institute of Geology，China Earthquake Administration，Beijing（in Chinese）．

胡俊，李志伟，朱建军，等．2010．融合升降轨SAR干涉相位和幅度信息揭示地表三维形变场的研究［J］．中国科学（D辑），40（3）：307—318．

HU Jun，LI Zhi-wei，ZHU Jian-jun，et al．2010．Inferring three-dimensional surface displacement field by combining SAR interferometric phase and amplitude information of ascending and descending orbits［J］．Science in China（Ser D），40（3）：307—318（in Chinese）．

屈春燕，宋小刚，张桂芳，等．2008．汶川MS8.0地震InSAR同震形变场特征分析［J］．地震地质，30（4）：1076—1083．

QU Chun-yan，SONG Xiao-gang，ZHANG Gui-fang，et al．2008．Analysis on the characteristics of InSAR coseismic deformation of the MS8.0 Wenchuan earthquake［J］．Seismology and Geology，30（4）：1076—1083（in Chinese）．

屈春燕，单新建，张桂芳，等．2009．四川汶川MS8.0地震同震干涉形变场定量分析［J］．自然科学进展，19（9）：963—974．

QU Chun-yan，SHAN Xin-jian，ZHANG Gui-fang，et al．2009．Quantitative analysis of the coseismic inerferometric deformation field associated with the MS8.0 Sichuan Wenchuan earthquake［J］．Progress in Natural Science，19（9）：963—974（in Chinese）．

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

SHAN Xin-jian，QU Chun-yan，SONG Xiao-gang，et al．2009．Coseismic surface deformation caused by the Wenchuan MS8.0 earthquake from InSAR data analysis［J］．Chinese Journal of Geophysics，52（2）：496—504（in Chinese）．

孙建宝，梁芳，徐锡伟，等．2006．升降轨道ASAR雷达干涉揭示的巴姆地震（MW6.5）3D同震形变场［J］．遥感学报，10（4）：489—496．

SUN Jian-bao，LIANG Fang，XU Xi-wei，et al．2006.3D coseismic deformation field of Bam earthquake（MW6.5）from ascending and descending pass ASAR radar interferometry［J］．Journal of Remote Sensing，10（4）：489—496（in Chinese）．

GE L L，Rizos C，Han S，et al．2000．The double interpolation and double prediction（DIDP）approach for InSAR and GPS integration［J］．Int Arch Photogram Rem Sens（IAPRS），33：205—212．

Gudmundsson S，Sigmundsson F，Carstensen J．2002．Three-dimensional surface motion maps estimated from combinedinterferometric synthetic aperture radar and GPS data［J］．Journal of Geophysical Research，107（B10）：2250．doi：10．1029／2001JB000283．

Massonnet D，Rossi M，Carmona C，et al．1993．The displacement field of the Landers earthquake mapped by radar interferometry［J］．Nature，364（8）：138—142．

Samsonov S，Tiampo K．2006．Analytical optimization of a DInSAR and GPS data set for derivation of three-dimensional surface motion［J］．IEEE Geosci Remote Sens Lett，3（1）：107—111．

Samsonov S，Tiampo K，Rundle J，Li Z H．2007．Application of DInSAR-GPS optimization for derivation of fine-scale surface motion maps of Southern California［J］．IEEE Trans Geosci Rem Sens，45（2）：512—521．

Wang Q，Qiao X，Lan Q，et al．2011．Rupture of deep faults in the 2008 Wenchuan earthquake and uplift of the Longmen Shan［J］．Nature Geoscience．doi：10．1038／NGEO1210．

Wright，T J，Parsons B E，Lu Z．2004．Toward mapping surface deformation in three dimensions using InSAR［J］．Geophys Res Lett，31：L01607．doi：10．1029／2003GL018827．

Abstract

The variance component estimation method（VCEM）in generalized surveying adjustment theory，which realizes optimal weights allocation for different data sources，is applied to jointly invert two independent datasets，InSAR and GPS，for 3-D deformation field acquisition in this paper．Illustrated by the case of the Wenchuan earthquake，3-D deformation field in high coherent area near the fault is achieved by using this method，which shows clearly a whole picture for the locations of right-lateral and thrust components movements generated by the earthquake．The 3-D deformation results are quantitatively consistent with GPS observations with RMS errors less than 5cm in 3-D directions，which demonstrates the feasibility to acquire precise 3-D deformation field by employing VCEM to fuse independent deformation datasets．

COSEISMIC 3D DEFORMATION FIELD ACQUISITION OF THE WENCHUAN EARTHQUAKE BASED ON INSAR AND GPS DATA

SONG Xiao-gang1）SHEN Xing1，2）JIANG Yu2）WAN Jian-hua2）SHAN Xin-jian1）QU Chun-yan1）
1）State Key Laboratory of Earthquake Dynamics，Institute of Geology，China Earthquake Administration，Beijing100029，China
2）School of Geosciences，China University of Petroleum（East China），Qingdao266580，China

Wenchuan earthquake，3-D deformation field，variance component estimation，InSAR，GPS

P315．72＋5

A文献标识码：0253－4967（2015）01－0222－10

10．3969／j．issn．0253－4967．2015．01．017

宋小刚，男，1979年生，副研究员，主要从事雷达干涉测量技术（InSAR），电话：010－62009095，E－mail：sxghohai＠163.com。

2013－04－07收稿，2014－12－10改回。

中国地震局地质研究所基本科研业务专项（IGCEA1010）和国家自然科学基金（41111140386）共同资助。

申星，男，在读硕士研究生，E－mail：shenxing＿cumt＠126．com。