因人为或自然因素导致的滑坡、地面沉降、泥石流等地质灾害已经严重影响生命安全和财产安全，特别是因大范围工程施工或地下水过度抽取而引发的地表形变带来的影响已不容忽视^[1]。GPS因其具有高时间分辨率、可获得高精度三维形变等优点已被广泛应用于形变监测等领域^[2]。另一种空间对地观测技术——InSAR，因其几乎不受天气影响继而可获得大范围、高精度厘米级甚至毫米级形变也备受青睐^[3]。然而在实际应用中两种技术均存在先天不足，例如GPS的空间分辨率低、建站成本高，而InSAR仅能反映出视线(LOS)向相对形变量等。能否将GPS和InSAR进行集成和融合，以综合各自优势反演高精度多维形变场是值得关注和研究的问题^[3-4]。

在融合GPS和InSAR监测结果解算多维形变场方面已有众多成果^[5-11]。2002年，文献[5]使用Markov随机场理论及模拟退火法，建立了GPS与InSAR解算地表形变的函数模型，基于能量函数最小原则获得最优解，并较好应用于冰岛Reykjanes Peninsula地区。此方法计算时间长，解算结果并不能保证全局最优。2006年，文献[6-7]简化了能量方程求解算法，通过对方程求导，获得地表形变的最优解，但此方法仍有解算参数不稳定的情况。2008年，文献[8]用直接解法将GPS的水平形变插值，结合降轨InSAR获得的LOS向形变获取垂向形变，此方法的水平形变对GPS插值结果依赖性强。文献[9]于2012年研究了InSAR及GPS两类观测量的方差分量估计确定权比，解算三维形变场。2013年，文献[10]基于目标函数是凸二次函数的特点，利用BFGS方法巧妙地避免了局部最优解。但以上方法仅基于一种InSAR观测数据。

基于上述背景，本文采用一种简洁有效的融合方法，即利用GPS插值南北向形变场及升降轨InSAR两类观测数据获取形变区的三维形变场，并以西安市为例，验证融合方法。西安市位于汾渭盆地，该地区因过度抽取承压水和大面积施工建设而导致的地表形变一直备受关注，而精确获取西安市地表多维形变场对认识其复杂的地质环境尤为重要。以往文献更多关注西安市地面沉降^[12-13]，而鲜有大范围反演多维形变场。考虑到南北向形变对InSAR观测量的不敏感性，本文基于ALOS PALSAR(PALSAR)和EnviSat ASAR(ASAR)两种数据源和GPS观测数据，首先由GPS观测信息，获取研究区域南北向形变场，结合西安地区2009-2010年两种InSAR数据源提取的LOS向年形变速率，解算出该区域东西向和垂向形变速率场(形变场)，并通过同期水准数据对比，其垂向形变精度呈毫米级。

1 形变速率场解算原理 1.1 GPS和InSAR融合策略

InSAR观测量仅反映了实际形变矢量在LOS向的投影大小^[14]，为得到多维形变场，首先分析N、E、U方向形变量与LOS向形变量之间关系。假设d_Nth、d_Est、d_Up为地表形变在地平坐标系下坐标轴方向的分解量，d_LOS为LOS向形变量，根据图 1可建立如下关系式^[15-16]

式中，d_LOS靠近卫星方向为+；d_Nth、d_Est、d_Up指向N、E、U方向为+；α为卫星飞行方向方位角；θ为入射角。

理论上，在没有任何先验信息情况下，解算三维形变量d_Nth、d_Est和d_Up至少需要3种数据源。目前国际上常用数据来源于ALOS PALSAR、EnviSat ASAR和TerraSAR，其传感器参数如表 1所示。

将表 1中参数代入式(1)，得到

式中，d_{LOS_A}、d_{LOS_E}、d_{LOS_T}分别为PALSAR、ASAR、TerraSAR观测的LOS向形变，将上述方程整理成误差方程形式

式中，；=[d_Nth d_Est d_Up]′；L=[d_{LOS_A} d_{LOS_E} d_{LOS_T}]′；V为观测值的残差向量。由上式即可解求三维形变场。

然而，如对于线性方程组Ax=l，A的条件数度量了x对l误差或不确定度的敏感性，其数学定义为：A范数与A逆范数的乘积，如A的条件数较大，l的微小变化就能引起x较大的改变，接近1时为良好条件矩阵，数值稳定性好。式(3)的系数阵条件数为130.2，已呈现一定程度的病态，故直接解算会带来N、E、U形变分量解的不稳定。考虑到InSAR观测量对南北向形变不敏感，为此有学者提出忽略南北向形变^{[12, 14]}，如此可以减少一维参数，改善方程性质。

在实际应用中，一方面大部分研究区域不足3种数据源，另一方面对于南北向具有较大形变的情况，忽略其形变将在一定程度上影响其余两方向的解算精度。顾及到GPS结果平面精度高，如果用GPS点在南北向内插结果d_NthGPS代替InSAR南北方向的形变量d_Nth，然后用ASAR和PALSAR数据解算其余两方向形变，可望能在一定程度上改善解算精度，于是有

此时，系数矩阵条件数1.71与式(3)的系数阵条件数130.2相比降低了两个数量级，方程解较为稳定。

1.2 InSAR形变速率场解算

InSAR解算结果为形变矢量在LOS向的投影，且是与参考像素的相对值。由于不同数据源选择的参考像素一般不同，并且参考像素的形变量不为零值，因此其解算的LOS向形变场与真实形变场之间存在着系统偏差，即参考点的形变量，如此在进行融合解算之前需先确定此偏差常数^[17]。如将形变场中GPS点形变量代替所在像素的形变，存在下列方程

式中，d_{NthGPS_i}、d_{EstGPS_i}、d_{UpGPS_i}为PALSAR形变场GPS点所在像素的绝对形变量；d_{LOS_Ai}为i像素PALSAR形变场的LOS向观测量；d_{NthGPS_j}、d_{EstGPS_j}、d_{UpGPS_j}为ASAR形变场GPS点所在像素j的绝对形变量；d_{LOS_Ej}分别是ASAR在j像素的LOS向形变量；C_N、C_E、C_U和C_N′、C_E′、C_U′分别是PALSAR和ASAR两种数据源的参考点在N、E、U方向的形变量。用D_{LOS_A}、D_{LOS_E}表示PALSAR和ASAR的系统偏差组合

分析式(5)，两种InSAR数据源解算的形变场中均有n个GPS点，可列2n个方程解算2个未知数。为避免GPS点可能含粗差，采用抗差估计方法^[18-19]，即可降低精度低GPS点的权重，从而获得最佳估值和。

由于GPS水平观测值精度较高，可插值GPS的南北向观测量至InSAR相同的时空分辨率，结合PALSAR、ASAR数据源的LOS向观测结果可解算东西和垂向形变场，即

式中，左边为已知量，右边为未知参数。其他参数与式(5)中参数代表含义相同。如此可以解算每个像素的东西向和垂向形变量。

2 形变速率场解算试验

西安位于中国西北部的渭河盆地，地面沉降最早发现于1959年，随后因地面承压水过度抽取导致地面沉降快速发展，在1996年部分地面最大年沉降速率达到191 mm^[13]，随后黑河引水，沉降速率有所减缓，而后因大面积建筑施工，沉降中心转移至南郊和西郊^{[12, 20]}。

2.1 数据准备

西安地区有20个同期观测的GPS点及11个水准点，分布情况如图 2所示。解算出GPS点在3个方向的年形变速率和水准点的垂向形变速率。在众多内插方法中，Kriging插值基于测量点间基于测量点间统计关系，可预测表面模型并精度较高，成为地学统计中常用方法之一。因此，采用选取线性变异函数的Kriging插值，将GPS点的南北和东西向形变量插值到InSAR格网。

关于InSAR数据源，现有17景C波段的ASAR降轨数据，时间跨度2009-01-03-2010-09-25，空间基线阈值为280 m，选择34个干涉对采用小基线和stacking^[21]方法提取LOS向年形变速率，干涉图分布及解算形变结果如图 3(a)、(c)所示。另外，有12景L波段的PALSAR升轨数据，时间跨度2009-02-04-2010-12-26，选择2010-02-07作为参考影像，空间基线阈值为2500 m，选择36个干涉对用小基线和stacking方法解算LOS向年形变速率，结果如图 3(b)、(d)所示。根据形变场及GPS点分布情况选择研究范围，其像素大小为717×1317，像元分辨率为22.56×2.56 m，统一坐标基准及空间分辨单元，并将相干性低于0.3的失相干区域进行Kriging插值，结果如图 3(e)、(f)所示。

2.2 试验部分

为计算西安市多维绝对形变量，需先解算组合系统偏差，然后结合GPS内插结果得出东西向和垂向形变量。以水准观测作为真值，比较解算精度。

2.2.1 系统偏差解算

根据式(5)，结合表 2观测数据，得系统偏差组合估值、分别是29.2 mm和-4.8 mm。

2.2.2 三维形变场提取

获取系统偏差估值后，结合GPS与InSAR解算结果提取东西向和垂向年形变速率。此时，为验证本文所提出方法的有效性，分别采用3种方案进行解算。方案1：将GPS点南北向观测量内插，结合升降轨InSAR的LOS向观测量解算东西和垂向形变量；方案2：忽略南北向形变，即，将方案1中GPS点南北向形变置为零，解算东西和垂向形变量；方案3：将GPS点南北向和东西向观测量内插解算两方向形变场，选用相干性较好的PALSAR解算垂向形变。

由方案1及方案2均可解算东西向和垂向形变场，方案3可解算垂向形变场。最后将同期水准观测量作为真值，分析其在垂向的精度，结果如表 3所示。

分析表 3及诸图可以看出：

(1) 方案1由于充分顾及了南北向形变场，不仅提取了东西向形变场，其垂向形变场精度也在毫米级，误差最大值4 mm，均方差2 mm，明显高于其他方案。方案2忽略南北向形变，虽也解出了东西向形变，但其垂直形变精度降低，水准点上误差最大值达16 mm，均方差9 mm。方案3尽管也顾及了南北向的形变场，但由于东西向形变场由GPS结果内插得到，降低了监测精度，因此垂向形变场解算精度也受到了影响。

(2) 对比分析方案1与方案2发现，图 5中两方案结果，因形变量较大并统一了形变范围，并不能直观看出两者区别。但是与水准结果相比，方案1提取的垂直分量结果明显更接近水准观测量，这也进一步证实了西安市南北向形变不可直接忽略。

(3) 方案1与方案3的结果对比说明了尽管GPS的平面精度较高，但是内插降低精度，其结果并不会与研究区域形变特征严格符合，如此，所提取的E、U两方向的形变量精度不如两系统联合解算得到的结果精度。

因此，考虑到InSAR观测量对南北向不敏感的特点，而研究区域南北形变又不可忽略时，需要借助外界数据或者寻求其他算法改正南北向形变量。

3 形变速率场结果分析

西安市垂向形变分析。2009-2010年均形变速率最大的地区达到了101 mm，呈V字型分布的沉降区，有4个沉降中心：鱼化寨、电子城、三爻村、曲江新区，前3个沉降中心均位于高新技术开发区，开发区自1991年以来，建区累计面积已达到15.20 km²，竣工面积9.10 km²。曲江新区为陕西文化产业示范区，自2005年建成以来，核心区域面积达51.2 km²。因此，大面积的建筑施工造成了地层压密或地表形变，除此之外，无施工建设的区域也有较小范围的沉降。图中也显示形变区基本上位于地裂缝两侧之间。因此形变主要是施工建设的影响，地裂缝的存在破坏了地层形变的连续性。

西安市东西向形变分析。从图 5(a)所示的东西向形变速率场可以清晰识别出高新开发区3个形变中心，即鱼化寨、电子城和三爻村，其中，鱼化寨在东西向形变速率达到120 mm，远远大于其在垂向的60 mm的下沉速率，因此在关注垂直形变的同时，水平形变也需要给予充分的重视，这很可能是施工建设导致地层倾斜在东方向的表现。电子城则有小部分区域向东移动，主要还是以垂直运动为主。

西安市南北向形变分析。图 4(a)南北向形变场由GPS点插值得到，其精度完全取决于GPS点精度及插值模型与实际形变的符合程度。整个研究区域均是向南运动，图中显示XJA1所在的三爻村年形变速率最大，为25 mm，另一个形变中心位于鱼化寨，形变量超过15 mm。

西安形变场整体分析。西安位于渭河盆地，此盆地在整体上表现出一种不连续的逆时针旋转特征，并因东部甘青块体向东挤压、北部鄂尔多斯地块左旋、南部和东部的华北和华南地块的不协调SE向伸展，渭河盆地表现出向东南方向运动^[22]，但局部地层表面受建筑物压力，而表现出非均匀性运动。

4 结 论

三维形变场的提取有助于灾害机理的反演，InSAR技术尽管能提供高精度的形变监测信息，但只能提供LOS向的形变场，难以满足实际需求。考虑到InSAR观测量对南北向形变不敏感，而当南北向形变又不可被忽视的情况下，对现有的GPS南北向监测成果进行插值，进一步与不同数据源的InSAR监测数据进行融合解算，可望能获得整个研究区域的高精度三维形变场。

本文验证了GPS与InSAR融合达到较好反演效果，同时，从条件数的角度验证了函数模型较为稳定，计算出西安市2009-2010年间在三维绝对形变场，并且4个沉降中心主要分布在西安市的东南和西南方向，呈V字型分布，东西向形变集中在高新开发区，形变量之大应引起足够重视。