GPS可以获取高精度的三维形变信息，具有很高的时间分辨率，但空间分辨率较低^[1]。而InSAR具有很高的空间分辨率^[2]，但时间分辨率较低。如何有效融合两种对地观测手段以得到监测区域高分辨率的三维形变场，是近年来很多学者关注的问题。Gudmundsson等^[3]基于Markov随机场理论(MRF)和Bayes统计方法建立了获取地表形变的能量方程，但计算复杂，不能保证敛散性。对此，Samsonov等^[4]通过简化Gibbs能量方程，对每一变量求一阶偏导得到最优解析解，但解析优化法系数矩阵条件数很大，容易出现法方程病态，导致数值不稳定^[5]。针对这一问题，罗海滨等^[6]提出以GPS观测值为基准，对InSAR和GPS的精度进行评价，进而对二者定权，即改进的解析优化法。当水平方向权值很大时，称为直接分解法^[7]。但该方法在水平方向过度依赖GPS-Kriging插值精度，对垂直方向的精度提高也有限。

InSAR形变监测结果受到诸如卫星轨道误差、大气延迟误差、相位解缠误差、DEM误差等多种因素影响，从而导致形变监测结果产生一个系统性偏移量，给联合解的结果带来严重影响。由于不同误差源的系统误差相互交织在一起，实践中估计这种误差是有困难的。本文建立附加系统参数的综合形变分析模型，以GPS观测数据为基准，在现有GPS-InSAR数据联合解算模型中添加一个与位置有关的多项式函数作为系统误差整体改正，可望提高融合模型的解算精度。

1 GPS-InSAR数据联合解算模型

InSAR在某点i处的LOS向形变L_InSARⁱ与该点在e、n、u方向的三维形变v_eⁱ、v_nⁱ、v_uⁱ的关系如下：

(1)

式中，为InSAR在i点处的单位投影矢量矩阵，为i点处的待求三维形变量^[8]。

对任一个GPS实测点i，可列出观测方程：

(2)

式中，为GPS在i点处的三维形变量，为设计矩阵。

两类观测方程联立，可得：

(3)

令，则式(3)可写为：

(4)

2 附加系统参数的GPS-InSAR综合形变模型

对于InSAR数据中可能存在的系统误差，由于GPS、InSAR数据获取方式不同，处理流程也不同，如果按照误差的特点来确定附加的系统参数是有困难的。因此，在上述模型基础上，对InSAR观测方程加入一个跟点位置有关的二次多项式函数模型作为系统误差改正项：

(5)

式中，，是一个跟点的位置有关的二次多项式函数，为附加的系统参数。将式(5)与式(2)联立，组成新的综合形变模型：

(6)

式中，。则N个GPS观测点的观测方程为：

(7)

写成矩阵形式：

(8)

根据附加系统参数平差的基本原理，求解得到系统误差参数以后，对于任一点的三维形变都可以代入到观测方程(6)中进行求解。

3 系统参数显著性检验

系统参数函数BS的选择是附加系统参数平差的关键。虽然多项式模型不会出现参数的强相关现象，但拟合模型阶数的选择至关重要。为此，需要对附加参数进行统计假设检验，在一定的显著水平下确定所选BS项的合理性^[9]。

假设平差定权时先验单位权方差为σ²，引入附加系统参数后得到的残差为V，多余观测数为r，计算所得单位权方差估值为。对于引入附加参数的平差模型，原假设与备择假设分别为：

由于，因此对于给定的显著性水平α，若，则接受H₀，表明引入的附加参数合理^[10]。

4 算例分析 4.1 模拟数据算例

采用MOGI点源模型验证新方法的性能。用Matlab软件在400×450栅格网点上模拟地表三维形变场，如图 1所示。原始的水平方向形变形状对称，形变范围-30 ~30 mm；垂直形变在中心区域变化非常大，最大可达70 mm^[5]。

为尽量贴合真实的地表形变情况，在水平方向加入±6 mm的随机误差，较大的随机误差通常出现在GPS点稀少以及预期形变比较大的地方。垂直方向的随机误差略大于水平方向。在400×450网格点上随机选取100个GPS点，并提取每个点的三维形变。模拟SAR卫星的单位投影矢量S=[-0.328 0-0.063 7 0.943 6]，得到其LOS向形变。为了符合InSAR监测的实际情况以及更好地检验附加系统参数平差方法的可行性，加入标准差为0.5 mm的Gauss白噪声以及10 mm的系统误差，结果如图 2，其中黑色点为随机选取的100个GPS点的分布。

首先，利用普通Kriging插值方法将随机选择的100个GPS点内插至与InSAR相同的空间分辨率。从图 3可以看出，InSAR数据在加入10 mm系统误差之后，附加系统参数改正的新方法在E、N、U方向改正的系统误差量级大约为-1.66 mm、-0.31 mm、4.73 mm，对InSAR-LOS向形变的改正大约为-5.02 mm。同时计算的系统参数项BS为10.001 mm，极为接近初始所加入的10 mm系统误差，即能够完全改正系统误差的影响。

为定量分析不同方法的精度，分别采用GPS-Kriging插值法、直接分解法、解析优化法以及附加系统参数的平差方法解算不加入系统误差和加入10 mm系统误差的三维形变场。随机模型方面，观测量的权阵通过其标准差得到，其中InSAR形变观测量的标准差可以直接在SAR影像中开一个规则窗口，利用一个5×5的移动窗口进行估计；而GPS观测的标准差可通过GPS测量标准差和Kriging插值标准差计算得到^[11]。用上述4种方法求解得到的三维形变分量与模拟的无噪声三维形变场进行比较，在400×450个网格点上求得它们的均方根误差(RMSE)，结果列于表 1(单位mm)。

从表 1可以看出，不加入系统误差时，GPS-Kriging插值所得RMSE最大，插值精度较低；直接分解法在水平方向的精度跟GPS-Kriging插值相同，只在垂直方向有一定提高；解析优化法和附加系统参数的平差方法所得精度对只含有随机误差的数据有一定改善，但效果不是很明显。InSAR加入10 mm系统误差后，GPS-Kriging插值结果不受影响；但直接分解法对于加入系统误差的InSAR数据没有抵抗力，垂直方向误差最大；解析优化法是线性方程组整体联合求解，加入系统误差后水平方向精度降低。相比于其他方法，新方法在E、N、U方向均有提升，特别是垂直方向，相比于直接分解法跟解析优化法分别提升了96%和95%，这主要是因为InSAR在垂直向的单位投影矢量约为0.94，对垂直方向的形变最为敏感。从表中可以看出，附加系统参数改正的新方法能够最大程度地减弱系统误差的影响，结果跟不加入系统误差时的精度大体相同，这与100个GPS点处的拟合残差分析结果一致。

4.2 西安地区实测数据算例分析

附加系统参数的GPS-InSAR综合形变模型被应用于西安地区2009年的年均三维形变场研究中。InSAR数据采用ENVISAT_ASAR(降轨)卫星经GAMMA软件差分处理得到的2009年年均LOS向形变，该卫星的单位投影矢量大约为S =[0.3552-0.0754 0.9318]^[12]。同时，还获取了该地区20个GPS台站与InSAR观测量同步的形变监测结果，其外业数据监测采用9台带有扼流圈天线的双频GPS接收机按静态相对定位模式进行同步观测，并采用网连式的构网作业方式进行观测，连续观测时间达到72 h以上；内业采用高精度GPS定位软件GAMIT进行基线解算，以及通过诸如粗差分析、统计检验等特殊的数据处理方法进行GPS监测数据处理，以确保监测数据质量。平差得到各网点的形变监测结果，并经坐标转换得到其在E、N、U方向的形变量。为验证所提方法的有效性，计算时在研究区域均匀选取6个GPS点作为外部检核点且不参与计算，其余14个点参与不同方案的计算分析。图 4给出了以InSAR-LOS向形变为底图的14个参与解算的GPS实测点的点位分布。

将GPS-InSAR联合解算基本模型(LS)与附加系统参数的GPS-InSAR综合形变模型(LS+BS)两种方法在各方向的拟合残差(图 5)作对比。可以看出，LS+BS方法所得拟合残差整体上比LS方法更小，其在E、N、U 3个方向上系统误差改善分别为2.19 mm、0.46 mm和5.54 mm，InSAR-LOS向改正的系统误差约为5.98 mm，验证了InSAR对高程形变最为敏感，对南北向形变感知能力最弱。

为检验新模型的合理性，对其进行统计假设检验。给定原假设与备择假设：

由于，给定显著性水平α=0.05，采用χ²分布进行检验。计算可知，2.180 ，因此接受H₀，引入的附加参数合理。

同样，采用GPS-Kriging插值、直接分解法、解析优化法以及附加系统参数的平差方法进行三维形变场的解算，并用上述方法分别计算6个GPS外部检核点的三维形变，将其结果与实测GPS观测值进行比较，得到6个GPS台站的均方根误差(RMSE)，结果列于表 2(单位mm)。

从计算结果可以看出，由于GPS点较少且分布不均匀，GPS-Kriging插值所得RMSE最大，插值精度较低，而且GPS垂向精度没有水平方向高，西安地区又主要以地面沉降为主^[12]，所以U方向的精度最差；直接分解法在U方向精度提高有限，E、N方向不作改变；而新方法在E、N、U 3个方向上都与检核点GPS测量值最为接近，特别是在U方向，相比于Kriging插值法和解析优化法分别提高了8.6 mm和5.4 mm，这与模拟数据计算所得结果相同。

5 结语

本文结合近代测量数据处理中附加系统参数平差的方法，提出以高精度GPS观测值为基准，采用函数拟合的方法对InSAR观测值中的系统误差进行改正，以期能有效融合GPS与InSAR两种技术，提高地表三维监测的精度。新方法在模拟数据和西安地区的实测数据求解三维形变时都取得了更好的效果，相比于GPS-Kriging插值法、解析优化法、直接分解法，在E、N、U 3个方向的精度均有提高，特别是在U方向更为明显。通过假设检验确定系统参数设置的合理性，进一步验证了该方法的有效性。