InSAR作为当前重要的空间对地形变观测手段，能够获取高精度区域地表形变及大范围地壳运动特征等信息，较常规的地壳形变监测方法能更清晰地反映观测区的整体形变情况。但因存在视线(LOS)向模糊的问题，InSAR只能获取雷达脉冲入射方向的一维地表形变量，无法满足地表三维形变的定量描述。GPS技术具有高时间分辨率、高精度的地表三维监测优势，近年来利用InSAR-GPS技术融合解算三维形变场一直是InSAR技术监测地表形变研究的热点。已有学者在不同对地观测手段融合监测三维形变场的方面进行了研究^[1-3]，为InSAR和GPS两种类型的数据融合提供了很好的思路。

本文在利用InSAR技术监测地表三维形变场时，将研究区GPS监测结果作为先验信息，建立附加随机模型约束的形变模型，以解决InSAR技术LOS向观测量三维形变信息不足的问题。同时考虑到SAR技术对监测点南北向形变信息不敏感的问题，以GPS监测的南北信息作为南北向形变的强约束条件，以提高南北向形变量的解算精度。理论上，附有双约束的参数解要优于采用单一约束获得的参数解，因此建立随机模型与函数模型双约束的模型进行GPS与InSAR数据融合，以期提高计算结果精度与准确性。

1 函数模型与随机模型双约束的GPS-InSAR融合模型

InSAR在某一点P处的LOS向形变观测量为L_InSAR，根据雷达成像的几何关系建立观测方程:

$ {\mathit{\boldsymbol{V}}_{{\rm{InSAR}}}} = {\mathit{\boldsymbol{A}}_{{\rm{InSAR}}}}\mathit{\boldsymbol{\hat X}} - {\mathit{\boldsymbol{L}}_{{\rm{InSAR}}}} $

(1)

式中，A_InSAR=[S_x S_y S_z]为方程设计矩阵，由InSAR在P点处的单位投影矢量组成；$\mathit{\boldsymbol{\hat X}} = \left[ {{v_E}\;{v_N}\;{v_U}} \right]$为P点处的三维形变量，V_InSAR为观测值残差。若有3个以上从成像几何差异明显的SAR干涉对中获取的LOS向形变观测值，即可解算P点三维形变。

若研究区同时有其他来源的观测资料，以GPS观测为例，则可将GPS观测获取的三维形变观测向量X=[L_GPS^E L_GPS^N L_GPS^U]^T作为先验形变信息，构成参数向量X的随机模型约束。若以GPS的三维形变量X作为伪观测向量，其协方差矩阵Σ_X作为先验协方差矩阵，P_X为X的权矩阵，结合式(1)可以列出附加随机模型约束的误差方程：

$ \left\{ \begin{array}{l} {\mathit{\boldsymbol{V}}_{\mathit{\boldsymbol{\bar X}}}} = \mathit{\boldsymbol{\hat X}} - \mathit{\boldsymbol{\bar X}}\\ {\mathit{\boldsymbol{V}}_{{\rm{InSAR}}}} = {\mathit{\boldsymbol{A}}_{{\rm{InSAR}}}}\mathit{\boldsymbol{\hat X}} - {\mathit{\boldsymbol{L}}_{{\rm{InSAR}}}} \end{array} \right. $

(2)

则依据经典最小二乘贝叶斯估计原理可得：

$ \mathit{\boldsymbol{\hat X}} = \mathit{\boldsymbol{N}}_{bb}^{ - 1}\mathit{\boldsymbol{G}} $

(3)

$\mathit{\boldsymbol{\widehat X}}$的协因数矩阵为：

$ {\mathit{\boldsymbol{Q}}_{{{\hat X}_S}}} = \mathit{\boldsymbol{N}}_{bb}^{ - 1} $

(4)

式中，N_bb= P_X+ A_InSAR^TP _LInSARA_InSAR，G = P_XX+A_InSAR^T P _LInSARL_InSAR，P_LInSAR为InSAR在LOS向观测值对应的权阵。

考虑到目前SAR卫星的极轨运行方式导致LOS向观测几乎不能捕捉测点南北向的形变信息，顾及GPS结果平面精度较高，考虑用GPS点在南北向的形变信息补偿InSAR南北向形变信息的不足。结合式(1)，可以列出附加函数模型约束的误差方程：

$ \left\{ {\begin{array}{*{20}{l}} {{\mathit{\boldsymbol{V}}_{{\rm{InSAR}}}} = {\mathit{\boldsymbol{A}}_{{\rm{InSAR}}}}\mathit{\boldsymbol{\hat X}} - {\mathit{\boldsymbol{L}}_{{\rm{InSAR}}}}}\\ {\mathit{\boldsymbol{C\hat X}} - \mathit{\boldsymbol{W}} = 0} \end{array}} \right. $

(5)

式中，C为约束条件方程的系数矩阵，W为约束条件方程的常数项。依据经典最小二乘估计原理, 可得：

$ \begin{array}{l} \mathit{\boldsymbol{\hat X}} = \left( {\mathit{\boldsymbol{N}}_b^{ - 1} - \mathit{\boldsymbol{N}}_b^{ - 1}{\mathit{\boldsymbol{C}}^{\rm{T}}}\mathit{\boldsymbol{N}}_c^{ - 1}\mathit{\boldsymbol{CN}}_b^{ - 1}} \right) \cdot \\ \left( {\mathit{\boldsymbol{A}}_{{\rm{InSAR}}}^{\rm{T}}{\mathit{\boldsymbol{P}}_{{\rm{InSAR}}}}{\mathit{\boldsymbol{L}}_{{\rm{InSAR}}}}} \right) - \mathit{\boldsymbol{N}}_b^{ - 1}{\mathit{\boldsymbol{C}}^{\rm{T}}}\mathit{\boldsymbol{N}}_c^{ - 1}\mathit{\boldsymbol{W}} \end{array} $

(6)

$ {\mathit{\boldsymbol{K}}_S} = \mathit{\boldsymbol{N}}_c^{ - 1}\left( {\mathit{\boldsymbol{CN}}_b^{ - 1}\mathit{\boldsymbol{A}}_{{\rm{InSAR}}}^{\rm{T}}{\mathit{\boldsymbol{P}}_{{\rm{InSAR}}}}{\mathit{\boldsymbol{L}}_{{\rm{InSAR}}}} + \mathit{\boldsymbol{W}}} \right) $

(7)

$\mathit{\boldsymbol{\widehat X}}$的协因数矩阵为：

$ {\mathit{\boldsymbol{Q}}_{{{\hat X}_F}}} = \mathit{\boldsymbol{N}}_b^{ - 1} - \mathit{\boldsymbol{N}}_b^{ - 1}{\mathit{\boldsymbol{C}}^T}\mathit{\boldsymbol{N}}_c^{ - 1}\mathit{\boldsymbol{CN}}_b^{ - 1} $

(8)

其中，N_b=A_InSAR^TP_InSARA_InSAR，N_c=CN_b^－1C^T，K_S为联系数向量，其他符号的定义与前面一致。

当采用随机模型和函数模型进行双约束时，组成如下的误差方程:

$ \left\{ {\begin{array}{*{20}{l}} {{\mathit{\boldsymbol{V}}_{{\rm{InSAR}}}} = {\mathit{\boldsymbol{A}}_{{\rm{InSAR}}}}\mathit{\boldsymbol{\hat X}} - {\mathit{\boldsymbol{L}}_{{\rm{InSAR}}}}}\\ {{\mathit{\boldsymbol{V}}_{\bar X}} = \mathit{\boldsymbol{\hat X}} - \mathit{\boldsymbol{\bar X}}}\\ {\mathit{\boldsymbol{C\hat X}} - \mathit{\boldsymbol{W}} = 0} \end{array}} \right. $

(9)

依据最小二乘原理可得：

$ \begin{array}{*{20}{l}} {\mathit{\boldsymbol{\hat X}} = \mathit{\boldsymbol{N}}_{bb}^{ - 1}\left( {\mathit{\boldsymbol{G}} - {\mathit{\boldsymbol{C}}^{\rm{T}}}{\mathit{\boldsymbol{K}}_S}} \right) = }\\ {\mathit{\boldsymbol{N}}_{bb}^{ - 1}\mathit{\boldsymbol{G}} - \mathit{\boldsymbol{N}}_{bb}^{ - 1}{\mathit{\boldsymbol{C}}^{\rm{T}}}{\mathit{\boldsymbol{K}}_S} = {{\mathit{\boldsymbol{\hat X}}}_D}} \end{array} $

(10)

$ {\mathit{\boldsymbol{K}}_S} = \mathit{\boldsymbol{N}}_\alpha ^{ - 1}\left( {\mathit{\boldsymbol{CN}}_{bb}^{ - 1}\mathit{\boldsymbol{G}} - \mathit{\boldsymbol{W}}} \right) $

(11)

$\mathit{\boldsymbol{\widehat X}}$的协因数矩阵为：

$ {\mathit{\boldsymbol{Q}}_{{{\hat X}_D}}} = \mathit{\boldsymbol{N}}_{bb}^{ - 1} - \mathit{\boldsymbol{N}}_{bb}^{ - 1}{\mathit{\boldsymbol{C}}^{\rm{T}}}{\left( {\mathit{\boldsymbol{CN}}_{cc}^{ - 1}{\mathit{\boldsymbol{C}}^{\rm{T}}}} \right)^{ - 1}}\mathit{\boldsymbol{CN}}_{bb}^{ - 1} $

(12)

其中，N_bb=P_X+A_InSAR^TP_InSARA_InSAR，N_cc=CN_bb^－1C^T，G=P_XX+A_InSAR^TP_InSARL_InSAR。

通过上面的推导可得：

$ {\mathit{\boldsymbol{Q}}_{{{\hat X}_D}}} = {\mathit{\boldsymbol{Q}}_{{{\hat X}_S}}} - \mathit{\boldsymbol{N}}_{bb}^{ - 1}{\mathit{\boldsymbol{C}}^{\rm{T}}}{\left( {\mathit{\boldsymbol{CN}}_{bb}^{ - 1}{\mathit{\boldsymbol{C}}^{\rm{T}}}} \right)^{ - 1}}\mathit{\boldsymbol{CN}}_{bb}^{ - 1} $

(13)

则：

$ {\rm{tr}}\left( {{\mathit{\boldsymbol{Q}}_{{{\hat X}_D}}}} \right) \le {\rm{tr}}\left( {{\mathit{\boldsymbol{Q}}_{{{\hat X}_S}}}} \right) $

(14)

同理：

$ {\rm{tr}}\left( {{\mathit{\boldsymbol{Q}}_{{{\hat X}_D}}}} \right) \le {\rm{tr}}\left( {{\mathit{\boldsymbol{Q}}_{{{\hat X}_F}}}} \right) $

(15)

即理论上，双模型约束的参数解优于仅有随机模型约束或仅有函数模型约束的参数解。根据上述基本原理，使用随机模型约束、函数模型约束与双模型约束进一步求解地表的三维形变场。

2 算例分析 2.1 模拟数据进行算例分析

通过公式模拟100×100格网的三维形变速率：

$ \left[ {\begin{array}{*{20}{c}} {{V_x}}\\ {{V_y}}\\ {{V_z}} \end{array}} \right] = \left[ {\begin{array}{*{20}{c}} {V_x^0x}\\ {V_y^0y}\\ {\left[ {a + 2{\rm{ \mathit{ π} }}b\cos \left( {2{\rm{ \mathit{ π} }}t} \right)} \right]{{\rm{e}}^{ - \frac{{{x^2} + {y^2}}}{w}}}} \end{array}} \right] $

(16)

式中，V_x⁰、V_y⁰为水平方向的初始形变速率，a、b为常数，w为尺度参数，t为形变发生的时间间隔。本次研究计算的是年平均速率大小，因此取t=1模拟原始三维形变场。为尽量贴合真实的地表形变情况，首先设地表在东西向和南北向发生了匀速形变，而垂直向则作随时间推移的周期性下沉和抬升^[4]。同时，需要对模拟的数据进行加噪处理，在模拟的三维形变场的水平方向和垂直方向分别加入N(0, 0.25²)cm和N(0, 0.5²)cm的Gauss白噪声，结果见图 1，其中图 1(a)~1(c)为不加入噪声的模拟三维形变场，图 1(d)~1(f)为加入噪声的模拟三维形变场。

模拟研究区的D-InSAR升降轨投影矢量分别为[S_x^A S_y^A S_z^A]=(0.34, -0.095, 0.935)和[S_x^D S_y^D S_z^D]=(－0.34, 0.095, 0.935)^[5]，根据InSAR计算形变的几何原理即可得到升降轨D-InSAR视线向的形变，同时对升降轨InSAR数据加入N(0, 0.25²)cm的Gauss白噪声^[6]，结果如图 2所示。

在100×100的格网中随机选取100个GPS点位，提取点位的三维形变速率作为GPS实测数据，利用克里金插值法将GPS点内插至与模拟InSAR数据相同的空间分辨率中, 如图 3所示。克里金插值的精度主要取决于变异函数的选择，本文变异函数采用球形变异函数模型^[7]。

对模拟的GPS与InSAR数据采用3种方案进行融合计算三维形变场：

1) 方案1，以GPS的先验观测值为参数的随机模型条件，以InSAR的LOS向形变值为观测量，进行随机模型约束的GPS与InSAR数据融合(STMD)。

2) 方案2，以三维形变的南北向观测值作为条件约束，进行函数模型约束的GPS与InSAR数据融合(FNMD)。

3) 方案3，以GPS的先验观测值为参数的随机模型条件，同时以三维形变的南北向观测值作为条件约束，进行函数模型与随机模型双约束的GPS与InSAR数据融合(DCMD)。

观测量的权阵通过标准差计算得到，其中InSAR形变观测量的标准差是对模拟的SAR数据利用一个5×5的移动窗口进行估计，GPS观测的标准差是利用克里金插值的标准差计算^[3]得到。

利用克里金插值法、直接分解法及以上3种方法计算三维形变场，并将结果与原始模拟的三维形变场进行比较，求取100×100格网处的均方根误差(RMSE)，结果见表 1(单位cm)。

从表 1可以看出，直接分解法在E和N方向的精度与克里金插值法相同，在U方向的精度则优于克里金插值法；采用函数模型约束时，E方向的检核精度相较于其他几种方法明显较低，主要是因为在利用N方向进行强制约束时，E和U方向的形变信息主要来源于InSAR数据源，故平面精度低，而U方向精度相较于克里金插值法提升较多，进一步证明InSAR对于高程信息较为敏感；附加随机模型约束的三维形变解的外部检核精度均优于克里金插值法和直接分解法；双模型约束不仅实现了地表三维形变的解算，还使U方向上的形变解外部检核精度较常规的直接分解法提高近20%。

2.2 实测数据算例分析

利用随机模型约束、函数模型约束和双模型约束的方法对西安地区三维形变场进行计算，InSAR数据为ENVISAT卫星获取的10景C波段降轨ASAR数据(时间跨度为2009-12-18~2010-01-03)和ALOS卫星获取的8景L波段升轨PALSAR数据(时间跨度为2009-12-23~2010-12-20)。ALOS_PALSAR(升轨)和ENVISAT_ASAR(降轨)卫星经过GAMMA软件处理得到年均LOS向形变，升轨和降轨卫星的单位投影矢量分别为S_A=[-0.618 1-0.110 8 0.778 3]和S_D=[0.355 2-0.075 4 0.9318]。同时，获取了20个高精度的GPS观测点，采用GAMIT进行基线解算得到三维坐标，进一步转化得到观测点的E、N、U方向形变量。图 4为以升轨LOS向形变为底图的点位分布，为验证各方法的适用性，选取6个点位进行外部检核，点位分布见图 4(b)；剩余的14个点位参与融合模型的解算，点位分布见图 4(a)。

采用随机模型约束融合解算(STMD)、函数模型约束融合解算(FNMD)和双模型约束融合解算(DCMD)进行计算，解算过程中观测量的权阵通过标准差计算得到，其中InSAR形变观测量的标准差是对实测SAR数据利用一个5×5的移动窗口进行估计，GPS观测的标准差是利用GPS测量的标准差和克里金插值的标准差计算得到。3种方案的E、N、U方向拟合精度统计见表 2(单位mm)，E和U方向的拟合残差分布见图 5，3种方案计算参数协因数的迹见表 3。从表 2可以看出：1)双模型约束的解算结果在E和U方向的拟合精度优于单一的随机模型约束和函数模型约束；2)附加函数模型约束的解算结果在N和U方向的拟合精度较差，主要原因是该模型利用N方向的GPS观测结果对模型进行约束，而E和U方向的形变信息仅来源于InSAR数据源，导致解算结果平面精度较低，拟合残差较大。从表 3可以看出，双模型约束协因数的迹小于单一模型约束协因数的迹，与理论分析结果一致。从图 5可以看出，双模型约束在E和U方向的拟合残差小于函数模型约束和随机模型约束在E和U方向的拟合残差。

分别采用克里金插值法、直接分解法、单一随机模型约束、单一函数模型约束和双模型约束5种方法对6个外部检核点的三维形变速率进行解算，并将结果与实测GPS观测值进行比较。3个方向的外部检核精度统计结果见表 4(单位mm)。利用随机模型约束、函数模型约束和双模型约束求解研究区的三维形变场，图 6(a)~6(c)为随机模型约束计算的E、N、U方向三维形变速率，图 6(d)~6(f)为函数模型约束在3个方向的计算结果，图 6(g)~6(i)为双模型约束的计算结果。

从表 4和图 6可以看出：

1) 直接分解法在E和N方向的精度与克里金插值法相同，在U方向的精度则优于克里金插值法。

2) 在E和U方向，双模型约束与单一的随机模型约束的精度优于常规的直接分解法，其中U方向的精度提高最为显著。

3) 函数模型约束在E方向相较于其他几种方法精度较低，主要原因是在利用N方向进行强制约束时，E和U方向的形变信息主要来源于InSAR数据源，故平面精度较低。U方向精度相较于克里金插值法提升较多，进一步说明InSAR对于高程信息较为敏感。E方向精度较差，说明对于东西方向的形变不能单靠InSAR信息获得，而是要充分地利用GPS数据。

4) 双模型约束的精度优于单一函数模型约束和单一随机模型约束的精度，即双模型约束的精度优于单一模型约束的精度。但双模型约束较单一模型约束的精度提高有限，因为研究区GPS点数量非常有限且分布不均匀，对GPS的内插精度影响很大，且没有考虑InSAR观测数据中系统误差的影响。

5) 随机模型约束和双模型约束计算的三维形变场运动趋势基本保持一致，而函数模型约束与随机模型约束、双模型约束在E和U方向存在明显差异。研究区有着向东南方向运动的趋势，与渭河盆地向东南方向运动的结论一致^[8]。从图 6(c)和6(h)可以看出，研究区的西南方向主要发生沉降。

3 结语

本文针对InSAR技术研究地表三维形变时监测信息不足的问题，以GPS监测信息作为先验信息，建立附有随机模型约束的地表三维形变模型。考虑到InSAR卫星近南北向的极轨运行方式使观测量对南北向形变不敏感的问题，可能导致三维矢量模型系数矩阵呈现一定程度的病态，直接解算容易使N、E、U方向形变分量的结果不稳定，但研究区南北向形变又不可忽略，而GPS平面监测结果的精度较高，因此用GPS点在南北向的内插结果来代替InSAR不敏感的南北向形变分量，以提高三维形变解的精度。模拟数据和实测数据算例表明，附加函数模型约束能在一定程度上提高三维模型解算的精度，U方向尤为明显；双模型约束不仅实现了地表三维形变的解算，还使U方向上的形变解外部检核精度较常规的直接分解法有了显著提高，模拟数据体现更为明显；附加随机模型与函数模型双约束的精度优于单一函数模型约束或单一随机模型约束。