基于InSAR与GPS数据融合的地表三维形变监测方法

摘要

本发明公开了一种基于InSAR与GPS数据融合的地表三维形变监测方法，其步骤为：1.布设GPS点；2.采集GPS信号和收集SAR图像；3.解算GPS点；4.抑制SAR图像斑点噪声；5.建立GPS与SAR数据之间的坐标转换关系；6.利用GPS数据改正SAR卫星轨道误差；7.形成InSAR差分干涉图；8.利用GPS点反演大气延迟；9.改正InSAR大气延迟误差；10.解缠InSAR差分干涉图；11.将解缠相位转换为形变信息；12.InSAR形变信息地理编码；13.融合InSAR和GPS形变信息获取高精度的地表三维形变结果；14.利用卡尔曼滤波估计并获得高时空分辨率的地表三维形变结果。本发明由于采用了InSAR与GPS数据融合的地表三维形变监测技术，既突破单一监测技术的应用局限，又极大地改善了三维形变监测结果的时空分辨率。监测精度高，尤其是在垂直方向上的精度，同时系统的整体可靠性得到加强。

权利要求

1.一种基于InSAR与GPS数据融合的地表三维形变监测方法，其步骤为：步骤一、布设GPS点

布设的GPS点有两类：第一类是监测区内稳定的特征点，用于将GPS数据和InSAR数据统一到同一参考坐标系下、改正SAR卫星轨道误差以及反演大气延迟；第二类是监测区的监测点；

步骤二、采集并记录GPS信号和收集SAR图像采集的GPS信号要与SAR图像获取的时间同步；

步骤三、依据GPS信号解算出两类GPS点的三维坐标信息第一类GPS监测点的三维坐标信息包括：位置信息和高程信息，第二类GPS监测点的三维坐标信息包括：位置信息和形变信息；

步骤四、SAR图像的斑点噪声抑制

步骤五、搜索SAR图像上的特征点，并建立搜索到的特征点与第一类GPS点之间的坐标转换关系步骤六、利用GPS数据改正SAR卫星的轨道误差采用测图方程法和最小二乘估计相结合的方法；

步骤七、形成InSAR差分干涉图

步骤八、利用第一类GPS点反演大气延迟，并将其统一到InSAR的坐标系统中采用随机过程法获取GPS的对流层估计，采用站间时域双差法估计大气延迟改正；

步骤九、利用GPS反演的大气延迟改正InSAR差分干涉图的大气延迟误差采用融合GPS高程信息的改进反距离加权内插法对GPS反演的大气延迟进行加密，然后将加密后的大气延迟从InSAR差分干涉相位图中消除；

步骤十、对大气改正后的InSAR差分干涉图进行解缠采用GPS辅助InSAR Goldstein枝切线解缠算法或引入残差点的基于马尔科夫随机场GPS辅助InSAR解缠算法，对大气改正后的InSAR差分干涉图进行解缠；

步骤十一、将InSAR解缠相位转换为形变信息其中，□Rd为视线向的形变量，φd为解缠相位，λ为雷达波长；

步骤十二、InSAR形变信息地理编码

距离方程：

多普勒方程：

地球模型方程：

式中，R为SAR传感器到地面点间的距离， 分别为卫星和地面点的位置矢量，“·”为矢量点乘，fD为多普勒频移，为卫星速度矢量，为地面点速度矢量，xt，yt和zt为地面点位置矢量 在x，y和z轴上的投影，a为地球长半轴，b为地球短半轴，h为平均正常高；

通过对上述三个方程的求解即可得到InSAR形变信息地理编码；

步骤十三、采用基于改进马尔科夫随机场的GPS和InSAR融合方法，融合InSAR和GPS的形变信息获取高精度的地表三维形变监测结果；

最后利用卡尔曼滤波对格网中的所有点进行估计获得高时空分辨率的地表三维形变结果。

2.根据权利要求1所述的基于InSAR与GPS数据融合的地表三维形变监测方法，其特征在于：步骤十三中所述基于改进马尔科夫随机场的GPS和InSAR融合方法：首先，将InSAR视线向的一维形变结果通过GPS约束分解成三维形变信息；具体算法如下：(1)建立视线向的一维形变结果与三维形变结果之间的关系dlos＝sinθcos(αh-3π/2)dn+sinθsin(αh-3π/2)de-cosθduT

＝[un，ue，uu][dn，de，du]其中，dlos为InSAR视线向形变测量，[un，ue，uu]为视线向单位矢量，[un，ue，uu]＝[0.34，-0.095，0.935]，dn、de和du为北东天方向优化值，θ为雷达视角，其值随像元点至卫星星下点的距离而改变；αh为卫星飞行方向与北方向夹角，其值从SAR数据中读出；

(2)确定总能量函数

其中，

其中， 为InSAR约束权值， 为GPS约束权值， 为InSAR视线向形变测量， 为视线向单位矢量， 和 为北东天方向优化值， 为InSAR观测误差，和 为北东天方向克里格插值， 和 为北东天方向克里格插值标准差；

(3)采用交叉验证定权法确定约束权值，具体如下：对GPS与InSAR观测值进行配准，并将GPS三维观测值换算到InSAR视线方向；

按下式估计GPS观测站对应InSAR像元测量方差：式中，n为GPS观测站个数， 为第i个GPS观测站三维观测值换算到InSAR视线向结果， 为InSAR第i个像元视线向观测值；

InSAR每一像元点的测量方差由下式计算：式中， 和γi分别为InSAR第i个像元点测量方差估值和相干系数，E(γGPS)为GPS观测站对应InSAR像元相干系数均值；

将GPS观测值分成A、B两组，对A组进行普通克里格插值；

求B组GPS三维观测值与其对应插值的方差：式中，n′为B组GPS观测站个数， 和 为北东天方向B组第i个GPS观测站观测值， 和 为北东天方向B组第i个GPS观测站克里格插值；

计算克里格插值方差改正因子：

式中， 和 分别为B组GPS观测站北东天方向克里格插值方差均值；

对克里格方差进行改正：

式中， 和 为第i个像元点北东天方向克里格插值方差；定权：(4)采用直接解析法获取三维形变场的最大后验估计当方程系数矩阵不为零时，求得一组唯一解，该解即为高空间分辨率的三维形变场的最大后验估计；

然后，对高时间频率GPS数据在时间域上对上述三维形变场进行克里格插值和加密，从而将准GPS结果在时间域内加密成准GPS时间序列，获取高时间分辨率的三维形变。

说明书

基于InSAR与GPS数据融合的地表三维形变监测方法

技术领域

[0001] 本发明涉及一种空间对地监测地表三维形变的方法，具体说是一种基于InSAR与GPS数据融合的地表三维形变监测方法，适用于城市环境、滑坡和地震等区域高时空分辨率、高精度的地表三维形变监测。

背景技术

[0002] 近年来，随着卫星定位理论的日臻成熟以及接收设备性能的提高，全球定位系统(GPS)已广泛应用于形变监测。GPS可以提供高时间分辨率、高精度的三维形变测量值，而且不受天气的影响、点位间毋需通视、布站灵活、对地表状况适应能力强。但由于接收机数量和布网阵列限制，GPS空间分辨率较低、覆盖区域小，可能遗漏一些重大安全隐患。具体说来，现有GPS形变监测方法主要存在以下不足：(1)采用布网观测的方法，劳动强度大、监测系统成本高；(2)采用点观测的方法，因而监测数据的空间分辨率低，难以满足大区域变形分析和地质灾害预测的要求；(3)在深山峡谷或城市建筑物密集区，GPS接收机天线受到遮挡，这使得接收到的GPS卫星数减少，导致GPS定位精度大大下降，达不到形变观测精度要求。

[0003] 合成孔径雷达干涉测量(InSAR)是近年来迅速发展起来的地表探测新技术，已有越来越多的国家和地区利用InSAR的差分技术来探测诸如由矿山开采、地震、火山运动等引起的地面形变现象。与GPS技术相比，InSAR技术具有垂直监测精度高、监测范围大、空间近连续性等优点。但InSAR技术时间分辨率低且对大气延迟误差、卫星轨道误差、地表状况以及时态不相干等因素非常敏感，这造成了InSAR技术在地表形变探测应用中的困难。

发明内容

[0004] 本发明的目在于，克服现有技术存在的缺陷，综合利用InSAR与GPS技术，突破单一技术的应用局限，提供一种基于InSAR与GPS数据融合的地表三维形变监测方法。

[0005] 本发明的基本原理为：

[0006] InSAR技术是近年来迅速发展起来的地表探测新技术，其原理是通过处理接收到的不同时相地表回波来监测地表的形变。与GPS技术相比，InSAR技术具有垂直监测精度高、监测范围大、空间近连续性等优点。但InSAR技术时间分辨率低且对大气延迟误差、卫星轨道误差、地表状况以及时态不相干等因素非常敏感，这造成了InSAR技术在地表形变探测应用中的困难。

[0007] 本发明考虑将GPS与InSAR数据融合起来，突破单一技术应用的局限，发挥其各自优势，极大地改善空间域和时间域的分辨能力。但是将GPS与InSAR的数据融合存在以下需要解决关键技术：(1)利用GPS改正InSAR卫星轨道误差；(2)利用GPS改正InSAR数据的大气延迟影响；(3)GPS辅助InSAR相位解缠；(4)InSAR与GPS数据在时空域融合。

[0008] 考虑到要发明一种具有高精度的监测方法，就要对观测数据的误差进行改正，为此根据InSAR数据处理的流程，将GPS数据融合其中，利用GPS改正的InSAR误差，然后采用基于改进马尔科夫随机场的GPS和InSAR融合方法，实现InSAR与GPS数据在时空域的融合，获取高时空分辨率的地表三维形变信息。

[0009] 1、GPS改正InSAR数据的轨道误差

[0010] 在InSAR图像配准、相位到高程的转换和地理编码等步骤都要用到SAR卫星的轨道信息。可见，获取精确的轨道信息对提高InSAR测量的精度至关重要。由于InSAR数据自身带的轨道信息精度不是很高，如，ERS的基本轨道数据的精度在飞行方向2-4m，垂直轨迹方向1-2m，径向5m，不能满足测量精度的要求，为此，本发明引入GPS观测数据，利用GPS精确测得InSAR图像上特征地物(如河流和道路的交叉口，水塔等)或人工角反射器的地面坐标，再利用地理编码的逆过程求得SAR卫星精确轨道坐标，达到改正InSAR数据的轨道误差的目的。

[0011] 2、GPS反演大气延迟改正InSAR大气延迟误差

[0012] SAR传感器发射的电磁波在穿过大气时，传播速度将发生变化，传播路径也将发生弯曲，从而产生延迟现象。同一区域不同时间观测得到的SAR影像的大气延迟不完全相同，因而不可避免的对InSAR干涉相位图产生污染。为此，本发明利用GPS反演大气延迟，首先采用随机过程法获取GPS的对流层估计，其次采用站间时域双差法估计大气延迟改正，考虑到大气延迟与高程具有较强的相关性，根据地理空间相关性理论，提出采用融合GPS高程信息的改进反距离加权内插法对其加密，然后将加密后的大气延迟改正从InSAR干涉相位图中消除，从而改正InSAR大气延迟误差。

[0013] 3、GPS辅助InSAR相位解缠

[0014] 为了将相位转换为高程或形变信息，必须对干涉图中主值区间内的缠绕相位进行恢复得到其绝对相位，即要进行相位的解缠。残差点的存在使得解缠相位不连续甚至出现无法解缠的“孤岛”。为此，本发明引入少量的GPS控制点，以这些GPS控制点作为新的解缠起算点，消除“孤岛”并使整个区域的解缠相位连续。

[0015] 本发明根据解缠算法误差传播的特点，提出了两种不同的GPS辅助InSAR相位解缠的方法，即，GPS辅助InSAR Goldstein枝切线解缠算法和引入残差点的基于马尔科夫随机场GPS辅助InSAR解缠算法。

[0016] 4、基于改进马尔科夫随机场的InSAR和GPS数据融合方法

[0017] 对于像滑坡、地震这样的地质灾害，会在短时间内引起大范围较大的形变，InSAR一个月左右的重复观测周期是不够的，需要在时间域内插值估计。另外，GPS是基于有限点的观测，点的空间密度不能满足监测的需求，需要在空间域内插值估计。本发明采用基于改进马尔科夫随机场的InSAR和GPS数据融合方法，将InSAR视线向的一维形变结果通过GPS约束分解成三维形变信息，获得高空间分辨率的三维形变场；然后对高时间频率GPS数据(其所采集的地面沉降量)在时间域上对上述高空间分辨率的三维形变场进行克里格插值和加密，从而将准GPS结果在时间域内加密成准GPS时间序列，获取高时间分辨率的三维形变场；最后利用卡尔曼滤波对格网中的所有点进行估计获得高时空分辨率的地表三维形变结果，实现地表形变高时空分辨率的监测。

[0018] 本发明基于InSAR与GPS数据融合的地表三维形变监测方法，其步骤是：

[0019] 1、布设GPS点

[0020] 布设的GPS点，既要作为控制点，用于GPS和InSAR数据坐标转化和改正InSAR误差改正，又要作为监测点，对测区的重点部位进行监测。因此，布设的GPS点分为两类：第一类是监测区内稳定的特征点(称为第一类GPS点)(稳定的特征点的稳定是指不受形变影响)，如河流、道路的交叉处，这类点主要用于将GPS数据和InSAR数据统一到同一参考坐标系下、改正SAR卫星轨道误差以及反演大气延迟；第二类是监测区的重点监测部位(称为第二类GPS点)。

[0021] 在布设的第一类GPS点时，要注意点的分布，尽量使其均匀分布在测区内。

[0022] 2、采集并记录GPS信号和收集SAR图像

[0023] 因为要用GPS数据改正InSAR数据的轨道误差和大气延迟误差，采集的GPS信号要与SAR图像获取的时间同步，GPS信息的采集时间段最好能向前和向后顺延10分钟左右。

[0024] 3、依据GPS信号解算出两类GPS点的三维坐标信息

[0025] 两类GPS点分开解算，解算的第一类GPS点的三维坐标信息包括：平面位置信息和高程信息，第二类GPS点的三维坐标信息包括：平面位移和垂直沉降。

[0026] 4、SAR图像的斑点噪声抑制

[0027] 有效地抑制SAR图像的斑点噪声，有利于提高InSAR监测结果精度和确保监测结果的可靠性。本发明采用改进的Lee算法(是现有的算法)。

[0028] 5、搜索SAR图像上的特征点，并建立搜索到的特征点与第一类GPS点之间的坐标转换关系

[0029] 为了将InSAR和GPS的数据进行融合，必须将它们的数据统一到同一参考坐标系下。搜索InSAR图像上特征地物，并用GPS测定其地面坐标。假设特征点目标的影像坐标为(Irow，Icol)，其地面坐标为(Glat，Glan)。在统一投影关系下，两套坐标系统的变换可表达为：

[0030]

[0031] 或

[0032] L＝BX (2)

[0033] 其中，L为特征点目标的GPS坐标矩阵，B为其对应影像坐标矩阵，x是3×2阶的转换矩阵。因此，至少需要3个特征点数据的数据才可以确定转换矩阵中的6个参数。当有更多的特征点数据时，可用最小二乘平差求得6个转换参数。假设观测误差为 ，则观测方程为：

[0034] L＝BX+□ (3)

[0035] 转换矩阵的估计值为：

[0036]

[0037] 其中， ，D□是 的方差矩阵。最后，利用转换矩阵将GPS求得的大气延迟改正、轨道改正等数据平移、旋转转换为InSAR影像坐标系下的数据。

[0038] 6、利用GPS数据改正SAR卫星的轨道误差

[0039] 采用测图方程法和最小二乘估计相结合的方法，具体是：

[0040] 在轨道误差的改正中，噪声抑制后的SAR图像作为输入数据，搜索InSAR图像上稳定的特征地物，并用GPS测定其地面坐标。根据特征地物的雷达坐标(即特征地物到InSAR图像的斜距和视角)与其地面坐标的转换关系，构建地形测图方程，即

[0041]

[0042] 其中，X，Y，Z为特征点的地面坐标，a1，a2，…，c3为SAR卫星姿态角 构成的旋转矩阵元素，λ0是比例因子，r为征地物到SAR图像的斜距，θ为视角，Xs，Ys，Zs为SAR卫星的轨道位置。

[0043] 然后，将SAR卫星的轨道参数和基线参数表示成为时间的线性函数，即：

[0044] Xs＝Xs0+Xs(t-t0)

[0045] . (6)

[0046] B＝B0+B(t-t0)

[0047] 对地形测图方程线性化，可得观测误差方程：

[0048] V＝AX-L (7)

[0049] 式中，T

[0050] V＝[VX，VY，VZ]

[0051]

[0052]

[0053] L＝[LX，LY，LZ]

[0054] 其中，ΔXs0，ΔYs0，ΔZs0和 分别为t0时刻的轨道位置和姿态改正数；和 分别为相应的变化率改正数；Δλ0，ΔB0分别为t0时刻的比例尺
和基线改正数； 分别为相应的变化率改正数；Δφ0为相位常数改正数；a11，a12，…a317为各未知参数的偏导数；VX，VY，VZ分别为X，Y，Z的已知值与其近似值的差。

[0055] 由于误差方程含有17个未知数，因此至少需要6个以上的第一类GPS点解算这些未知参数。本发明采用最小二乘法估计观测误差方程中的17个未知参数，由此得到SAR卫星轨道数据误差的改正值。

[0056] 7、形成InSAR差分干涉图

[0057] 首先对主副SAR图像进行精配准，配准精度要达到1/8像元；其次将副图像重采样到主图像的坐标系中；再次最后将配准后的主副图像作复共轭相乘，生成干涉图；然后消除平地效应相位；最后采用两轨法或三轨法进行差分干涉处理，获得差分干涉图。

[0058] 8、利用第一类GPS点反演大气延迟，并将其统一到InSAR的坐标系统中

[0059] 利用GPS反演大气延迟，首先采用随机过程法获取GPS的对流层估计，其次采用站间时域双差法估计大气延迟改正，并将其统一到InSAR的坐标系统中。其中，站间时域双差法为：

[0060] 假设A点是SAR图像上的一个GPS点，B是同一幅SAR图像上的另外一GPS点，如j j果将GPS解算得到的A、B两点在历元j的对流层延迟分别记为DA、DB，则站间对流层延迟的差值为

[0061]

[0062] 以A为参考站，利用式(8)同理可推得其它GPS点的对流层延迟站间单差。

[0063] 假设两个站点A、B，在两个历元j(主图像获取时刻)、k(副图像获取时刻)由式(8)建立的两个单差分别为：

[0064]

[0065]

[0066] 由式(9)的两个单差可以得到历元双差为：

[0067]

[0068]

[0069] 由式(10)可以看出，估计大气改正的双差有两种形式：一种是站间时域双差，即先站间差然后再历元间差，另一种是时域站间双差，即先历元间差再站间差。因为站间差仅与单幅SAR图像有关，从而可以更自由、灵活的与其它具有站间差的SAR图像组合形成历元间差，所以，本发明采用更有利于实际应用的站间时域双差。

[0070] 9、利用GPS反演的大气延迟改正InSAR差分干涉图大气延迟误差

[0071] 考虑到大气延迟与高程具有较强的相关性，根据地理空间相关性理论，引入高程影响因子，构建空间异相关模型，提出采用融合GPS高程信息的改进反距离加权内插法对其加密。对GPS反演的大气延迟进行加密，然后将加密后的大气延迟从InSAR干涉相位图中消除。其中，融合GPS高程信息的改进反距离加权内插法的基本原理如下：

[0072] 反距离加权内插法是一种常用的空间插值方法，它假设两个事物的相似程度随着彼此间距离的缩短而增加。因此，反距离加权(IDW)内插法以估值点和观测点间的距离为权重进行加权平均，距离估值点越近的观测点赋予的权重越大。反距离加权(IDW)内插法可表示为：

[0073]

[0074] Z(x0，y0，z0)是点(x0，y0，z0)处的估计值，Z(xi，yi，zi)是点(xi，yi，zi)处的观测值，n是用于估计的观测值的个数，λi是点(xi，yi，zi)处观测值对应的权，di0(i＝1…N)是观测点i到估值点的距离，p(p＞0)为距离的幂，且有 实际应用中，通常p＝2。

[0075] 由式(11)可以看出，反距离加权内插法将估值点和观测点间的水平距离和高程差对估值点的影响看成是等同的。但实际并非如此，如大气延迟的插值问题，估值点和观测点间的水平距离和高程差对估值点的影响并不相同。为了体现估值点和观测点间的水平距离和高程差对估值影响的不同，根据地理空间相关性理论，引入高程影响因子α，构建空间异相关模型，提出融合高程信息的改进反距离加权内插法，具体的插值公式如下：

[0076]

[0077] 其中，Z(x0，y0，z0)是点(x0，y0，z0)处的估计值，Z(xi，yi，zi)是点(xi，yi，zi)处的观测值，n是用于估计的观测值的个数，λi是点(xi，yi，zi)处观测值对应的权，hi为点s(xi，yi，zi)处的高程，h0为点(x0，y0，z0)处的高程，b0、b1为待求拟合系数，di0 为观测点ih
到估值点的水平距离，di0 为观测点i与估值点的高程差，k1、k2分别为大气延迟的差值随水平距离和高程差变化关系的线性拟合曲线的斜率，abs(·)代表取绝对值。

[0078] 融合高程信息的改进反距离加权内插法是根据地理空间相关性理论，在空间自相关假设基础上，进一步假设多种地理要素的空间分布具有内在关联性，由一种要素的空间分布推测其它要素空间分布的空间异相关模型。因此，融合高程信息的改进反距离加权内插法更有利于当待估站点高程较周围采样站点高程有较大变化时插值精度的提高。

[0079] 10、对大气改正后的InSAR差分干涉图进行解缠

[0080] 本发明设计了两种GPS辅助InSAR解缠算法：GPS辅助InSAR Goldstein枝切线解缠算法和引入残差点的基于马尔科夫随机场GPS辅助InSAR解缠算法。这两种方法有各自的优缺点，GPS辅助InSAR Goldstein枝切线解缠算法采用的是精确的积分算法，解缠精度高且计算速度快，但只能在局部解缠，且解缠范围和精度受噪声水平影响较大；引入残差点的基于马尔科夫随机场GPS辅助InSAR解缠算法是一种全局解缠算法，抑噪能力强，但由于采用的是模拟退火寻优搜索算法，解缠速度较慢，且不可避免地出现陷入局部最优的可能，从而导致解缠的误差。因此，用户可根据不同的需要选择不同的解缠算法。两种解缠算法的基本原理具体如下：

[0081] (1)GPS辅助InSAR Goldstein枝切线解缠算法

[0082] 在利用Goldstein枝切线法进行相位解缠时，会出现这样一种情况：某一小区域被枝切线包围而与干涉图的其他区域完全隔离，这样一小区域被形象地称为“解缠孤岛”。“解缠孤岛”的形成原因很多，如水中小岛、植被包围的空地等。对于“解缠孤岛”，传统的处理方法有：(a)当孤岛范围较小时，可以在解缠时放弃对其进行解缠而仅对其他区域进行解缠，当求得其他区域的最终产品(如DEM或形变量)后，再利用适当的插值算法对孤岛区域进行插值，从而得到孤岛区域的最终产品。这种方法的缺点是忽略了孤岛上的相位信息，且插值过程会引入较大误差。(b)当孤岛范围较大时，在孤岛区域内重新选择一解缠起算点，对孤岛区域进行独立的相位解缠。这种方法虽然用到了孤岛上的相位信息，但由于采用了新的解缠起算点，使得孤岛区域与干涉图其他区域解缠基准不统一，造成最终产品解译上的困难。针对传统处理方法中存在的问题，由于InSAR形变测量精度通常是厘米级，而GPS形变测量的精度可达到毫米级，本发明假设离散的GPS形变观测量为已知值，采用GPS辅助Goldstein枝切线法进行InSAR相位解缠。

[0083] GPS辅助InSAR Goldstein枝切线解缠算法的原理为：

[0084] 引入离散GPS点，选一GPS点为起算点，求其它GPS点相对于起算点的相对高程或形变，并利用下式将其转化为解缠相位值：

[0085]

[0086] 其中，φtopou和φdefou分别为解缠地形干涉相位和解缠形变干涉相位，Δh和Δr分别为相对于总解缠起算点的相对高程和形变。

[0087] 然后以孤岛点为起算点重新解缠，直到所有孤岛点解算完成。

[0088] 此方法是一种局部解缠算法，具有解缠速度快、精度高、基准统一且能评价解缠精度的优点。

[0089] (2)引入残差点的基于马尔科夫随机场GPS辅助InSAR解缠算法

[0090] 残差点的识别与连接是Goldstein枝切线解缠算法的最大特点，它可以有效地隔离误差点，避免误差的全局传播。但Goldstein枝切线解缠可能产生孤岛，且枝切线上的像元点也不能得到解缠，降低了解缠的有效范围。Gudmundsson提出的基于马尔科夫随机场的GPS辅助解缠算法，没有考虑残差点的影响，会造成误差的全局传播。由此，本发明将残差点识别技术引入基于马尔科夫随机场的GPS辅助解缠算法，提出了引入残差点的基于马尔科夫随机场GPS辅助InSAR解缠算法，实现解缠精度与解缠范围的最佳综合。此方法将InSAR解缠过程转化为求能量函数U(K＝k|Y＝Iw)最小值的过程，即：

[0091]

[0092] 采用反复快速退火策略以较高速率和概率获得全局最优解，即给定一较大温度衰减率cool和终止温度Ts，当T＜Ts时，令T＝T0，反复快速退火，直到达到全局最优解。迭代时所用到的初始整周数矩阵k，初始解缠干涉图Iu，解缠干涉图I′u和比例系数r计算公式如下：

[0093]

[0094] 其中，k为初始整周数矩阵，round(□)为取整算子，IuGPS为由GPS插值结果反演的解缠干涉图，Iw为缠绕干涉图，λ为波长，Iu为初始解缠干涉图，I′u为解缠干涉图，r为比例系数。

[0095] 与Goldstein枝切线解缠算法先进行残差点识别再进行枝切线连接的方法不同，本发明采用的残差点识别和标记策略为：计算2×2大小的封闭路径单元累积值S，若S不为零则将所有4个像元点都标记成残差点，从而完成残差点的识别和标记。这样标记残差点的优势在于：(1)避免了传统残差点识别方法中人为指定残差点位置(左上角点)可能导致的残差点位置标记错误，从而更加有效地抑制误差的传播；(2)不再进行复杂的枝切线连接，从而提高计算效率。

[0096] 11、将InSAR解缠相位转换为形变信息

[0097]

[0098] 其中，□Rd为视线向的形变量，φd为解缠相位。

[0099] 12、InSAR形变信息地理编码

[0100] 这一步是将前面处理得到的雷达坐标系下的结果转化到地理坐标系下的过程，采用距离多普勒(R-D)定位模型进行。在R-D定位模型中，SAR影像上任意像元的位置由三个方程确定：距离方程、多普勒质心平面方程和地球模型方程。

[0101] 距离方程：

[0102] 多普勒方程：

[0103] 地球模型方程：

[0104] 式中，R为SAR传感器到地面点间的距离， 分别为卫星和地面点的位置矢量，“·”为矢量点乘，fD为多普勒频移，λ为雷达波长，为卫星速度矢量，为地面点速度矢量，当采用WGS-84坐标系统时，的数值为零，xt，yt和zt为地面点位置矢量 在x，y和z轴上的投影，a为地球长半轴，a＝6378.137km，b为地球短半轴，b＝6356.752km，h为平均正常高。通过对上述三个方程的求解即可完成测量成果的地理编码过程。

[0105] 13、采用基于改进马尔科夫随机场的GPS和InSAR融合方法，融合InSAR和GPS的形变信息获取高时空分辨率的地表三维形变监测结果

[0106] 首先，将InSAR视线向的一维形变结果通过GPS约束分解成三维形变信息，获得高空间分辨率的三维形变场。采用分解优化法获取最优结果，其中对GPS约束条件的权值确定采用交叉验证定权的方法。具体原理如下：

[0107] 构建InSAR视线向的一维形变结果与三维形变结果之间的关系：

[0108] dlos＝sinθcos(αh-3π/2)dn+sinθsin(αh-3π/2)de-cosθdu

[0109] (17)[0110] ＝[un，ue，uu][dn，de，du]T

[0111] 其中，dlos为InSAR视线向形变测量，[un，ue，uu]为视线向单位矢量，[un，ue，uu]＝[0.34-0.095，0.935]，dn、de和du为北东天方向优化值，θ为雷达视角，其值随像元点至卫星星下点的距离而改变；αh为卫星飞行方向与北方向夹角，其值可以从SAR数据中读出。

[0112] 确定利用马尔科夫随机场进行建模计算的总的能量函数：

[0113]

[0114]

[0115] 其中，

[0116]

[0117]

[0118] 式(18)中，U为总能量，dlosi为InSAR视线向形变测量，[uni，uei，uui]为视线向单i i i i ni ei ui位矢量，dn、de 和du 为北东天方向优化值，σins 为InSAR观测误差，dkrig 、dkrig 和dkrigni ei ui
为北东天方向克里格插值，σkrig 、σkrig 和σkrig 为北东天方向克里格插值标准差。

[0119] 式(18)中包含了4×M个非负项，因此，当所有非负项达到最小时，式(18)便达到i i i了全局最小。于是对dn、de 和du 求偏导并令其为零得，可得解析优化法的计算公式，即：

[0120]

[0121] 采用式(18)中的权值确定方式存在如何评价InSAR测量精度和克里格插值精度两个问题，为有效解决这两个问题，对GPS约束条件的权值，本发明采用交叉验证定权的方法确定权值，具体过程如下：

[0122] ●对GPS与InSAR观测值进行配准，并将GPS三维观测值换算到InSAR视线方向，按式(20)估计GPS观测站对应InSAR像元测量方差：

[0123]los

[0124] 式中，n为GPS观测站个数，dGPS，i 为第i个GPS观测站三维观测值换算到InSARlos视线向结果，dins，i 为InSAR第i个像元视线向观测值。

[0125] ●计算InSAR每一像元点的测量方差：

[0126]

[0127] 式中，σins，i2和γi分别为InSAR第i个像元点测量方差估值和相干系数，E(γGPS)为GPS观测站对应InSAR像元相干系数均值。

[0128] ●将GPS观测值分成A、B两组，对A组进行普通克里格插值，求B组GPS三维观测值与其对应插值的方差：

[0129]

[0130] 式中，n′为B组GPS观测站个数，dGPS，in、dGPS，ie和dGPS，iu为北东天方向B组第i个GPS观测站观测值，dkrig，in、dkrig，ie和dkrig，iu为北东天方向B组第i个GPS观测站克里格插值。

[0131] ●计算克里格插值方差改正因子：

[0132]2 2 2

[0133] 式中，E[(σkrig，n)B]、E[(σkrig，e)B]和E[(σkrig，u)B]分别为B组GPS观测站北东天方向克里格插值方差均值。

[0134] ●对克里格方差进行改正：

[0135]n 2 e 2 u 2

[0136] 式中，(σkrig，i)、(σkrig，i) 和(σkrig，i) 为第i个像元点北东天方向克里格插值方差。

[0137] ●确定权值：

[0138]

[0139]

[0140] 然后，对高时间频率GPS数据(其所采集的地面沉降量)在时间域上对上述高空间分辨率的三维形变场进行克里格插值和加密，从而将准GPS结果在时间域内加密成准GPS时间序列，获取高时间分辨率的三维形变场。

[0141] 最后利用卡尔曼滤波对格网中的所有点进行估计获得高时空分辨率的地表三维形变结果。

[0142] 本发明基于InSAR与GPS数据融合的地表三维形变监测方法，由于综合了InSAR与GPS这两种技术的优点，突破单一技术应用的局限，为地表监测提供了一种高时空分辨率、高精度的新技术，适合在大范围地表形变监测中使用。

附图说明

[0143] 图1：本发明方法的数据处理的流程图；

[0144] 图2：GPS改正InSAR数据轨道误差的流程图；

[0145] 图3：GPS反演大气延迟改正InSAR大气延迟误差的原理图；

[0146] 图4：“解缠孤岛”产生示意图；

[0147] 图5：考虑残差点的基于马尔科夫随机场GPS辅助InSAR相位解缠数据处理流程图；

具体实施方式

[0148] 下面结合实施例，对本发明作进一步的详细说明。

[0149] 实施例：

[0150] 1、布设GPS点

[0151] 本方法布设两类GPS点：第一类是监测区内稳定的特征点，如河流、道路的交叉处，这类点主要用于将GPS数据和InSAR数据统一到同一参考坐标系下、改正SAR卫星轨道误差以及反演大气延迟；第二类是监测区的重点监测部位。

[0152] 在布设的第一类GPS点时，要注意点的分布，尽量使其均匀分布在测区内。

[0153] 2、采集并记录GPS信号和收集SAR图像

[0154] 采集时间同步的GPS信号与SAR图像，采集的GPS信息相对于SAR图像的获取时间，向前和向后顺延10分钟左右。

[0155] 3、依据GPS信号解算出两类GPS点的三维坐标信息

[0156] 两类GPS分开解算，解算的第一类GPS点的三维坐标信息包括：平面位置信息和高程信息，第二类GPS点的三维坐标信息包括：平面位移和垂直沉降。

[0157] 4、SAR图像的斑点噪声抑制

[0158] 本方法采用改进的Lee算法抑制SAR图形啊的斑点噪声。

[0159] 5、搜索SAR图像上的特征点，并建立搜索到的特征点与第一类GPS点之间的坐标转换关系

[0160] 首先，搜索3个以上InSAR图像上特征点目标地物，并用GPS测定其地面坐标。然后，采用下式计算两套坐标系统之间的转换矩阵估计值：

[0161]

[0162] 其中，L为特征点目标的地面坐标矩阵，B为其对应影像坐标矩阵， D□是观测误差 的方差矩阵。

[0163] 最后，利用转换矩阵将GPS求得的大气延迟改正、轨道改正等数据平移、旋转转换为InSAR影像坐标系下的数据。

[0164] 6、利用GPS数据改正SAR卫星的轨道误差

[0165] 在轨道误差的改正中，首先，将噪声抑制后的SAR图像作为输入数据，搜索6个以上InSAR图像上特征点目标地物，并用GPS测定其地面坐标。

[0166] 其次，构建地形测图方程，即

[0167]

[0168]

[0169] 其中，X，Y，Z为特征点目标地物的地面坐标，a1，a2，…，c3为SAR卫星姿态角构成的旋转矩阵元素，λ是雷达波长，c为光速，r为特征点目标地物到SAR图像的斜距，θ为视角，Xs，Ys，Zs为SAR卫星的轨道位置，φ为解缠相位。

[0170] 然后，将SAR卫星的轨道参数和基线参数表示成为时间的线性函数，即：

[0171] Xs＝Xs0+Xs(t-t0)

[0172] B＝B0+B(t-t0)

[0173] 对地形测图方程线性化，可得观测误差方程：

[0174] V＝AX-L

[0175] 式中，

[0176] V＝[VX，VY，VZ]T

[0177]

[0178]

[0179] L＝[LX，LY，LZ]

[0180] 其中，ΔXs0，ΔYs0，ΔZs0和 分别为t0时刻的轨道位置和姿态改正数；和 分别为相应的变化率改正数；Δλ0，ΔB0分别为t0时刻的比例尺
和基线改正数； 分别为相应的变化率改正数；Δφ0为相位常数改正数；a11，a12，…a317为各未知参数的偏导数；Vx，VY，VZ分别为X，Y，Z的已知值与其近似值的差。

[0181] 采用最小二乘法估计观测误差方程中的17个未知参数，由此得到SAR卫星轨道数据误差的改正值。

[0182] 7、形成InSAR差分干涉图

[0183] 首先对主副SAR图像进行精配准，配准精度要达到1/8像元；其次将副图像重采样到主图像的坐标系中；再次最后将配准后的主副图像作复共轭相乘，生成干涉图；然后消除平地效应相位；最后采用两轨法或三轨法进行差分干涉处理，获得差分干涉图。

[0184] 8、利用第一类GPS点反演大气延迟，并将其统一到InSAR的坐标系统中

[0185] 首先采用随机过程法获取GPS的对流层估计，然后采用站间时域双差法估计大气延迟改正，并将其统一到InSAR的坐标系统中。其中，站间时域双差法为：

[0186]

[0187]

[0188] 其中，A点是SAR图像上的一个GPS参考点，B是同一张SAR图像上的另外一GPSj j点，j和k分别为主副图像获取时刻，DA 和AB 分别为将GPS解算得到的A、B两点在历元j的对流层延迟。

[0189] 9、利用GPS反演的大气延迟改正InSAR差分干涉图大气延迟误差

[0190] 首先采用融合GPS高程信息的改进反距离加权内插法对GPS反演的大气延迟进行加密，然后将加密后的大气延迟从InSAR干涉相位图中消除。其中，融合GPS高程信息的改进反距离加权内插法为：

[0191]

[0192] 其中，Z(x0，y0，z0)是点(x0，y0，z0)处的估计值，Z(xi，yi，zi)是点(xi，yi，zi)处的观测值，n是用于估计的观测值的个数，λi是点(xi，yi，zi)处观测值对应的权，hi为点s h(xi，yi，zi)处的高程，b0、b1为待求拟合系数，di0 为观测点i到估值点的水平距离，di0 为观测点i与估值点的高程差，k1、k2分别为大气延迟的差值随水平距离和高程差变化关系的线性拟合曲线的斜率，abs(·)代表取绝对值。

[0193] 10、对大气改正后的InSAR差分干涉图进行解缠

[0194] 本方法设计了两种GPS辅助InSAR解缠算法：GPS辅助InSAR Goldstein枝切线解缠算法和引入残差点的基于马尔科夫随机场GPS辅助InSAR解缠算法。

[0195] 这两种方法有各自的优缺点，GPS辅助InSAR Goldstein枝切线解缠算法一种局部解缠算法，具有解缠速度快、精度高、基准统一且能评价解缠精度的优点，但解缠范围和精度受噪声水平影响较大；引入残差点的基于马尔科夫随机场GPS辅助InSAR解缠算法是一种全局解缠算法，抑噪能力强，但解缠速度较慢。用户可根据不同的需要选择不同的解缠算法。

[0196] 1)GPS辅助InSAR Goldstein枝切线解缠算法

[0197] 在GPS辅助InSAR Goldstein枝切线解缠算法中，数据处理的步骤如下：

[0198] (1)将GPS观测站配准到缠绕干涉图中。

[0199] (2)求缠绕相位梯度场主值、识别残差点并进行枝切线的连接。

[0200] (3)在位于孤岛外的GPS观测站中选定一总解缠起算点，如点1，求其它GPS点相对于点1的相对高程或形变，并利用下式将其转化为解缠相位值：

[0201]

[0202]u u

[0203] 上式中，φtopo 和φdefo 分别为解缠地形干涉相位和解缠形变干涉相位，Δh和Δr分别为相对于总解缠起算点的相对高程和形变。

[0204] (4)以点1为解缠起算点用Goldstein枝切线法对干涉图进行解缠，直到所有可解缠点被解缠。以图4所示情况为例，在这个过程中，孤岛外的所有点将被成功解缠，而孤岛内的点不能被解缠。

[0205] (5)对干涉图进行扫描，找到尚未被解缠的GPS点，即点2，并以该点为新的解缠起算点进行新一轮相位解缠。

[0206] (6)重复步骤5，直到整个干涉图解缠完毕。

[0207] (7)找到已解缠但未做过起算点的GPS点，即点3和点4，对其由步骤3反演得到的解缠相位已知值和由步骤4-5得到的解缠相位估计进行统计分析，从而实现对解缠算法的精度评价。

[0208] (8)将GPS反演的解缠相位已知值付给解缠干涉图上相应的点，整个解缠过程结束。

[0209] 2)引入残差点的基于马尔科夫随机场GPS辅助InSAR解缠算法

[0210] 此方法将InSAR解缠过程转化为求能量函数U(K＝k|Y＝Iw)最小值的过程，即：

[0211]

[0212] 采用反复快速退火策略以较高速率和概率获得全局最优解，即给定一较大温度衰减率cool和终止温度Ts，当T＜Ts时，令T＝T0，反复快速退火，直到达到全局最优解。迭代时所用到的初始整周数矩阵k，初始解缠干涉图Iu，解缠干涉图I′u和比例系数r计算公式如下：

[0213]

[0214] 其中，k为初始整周数矩阵，round(□)为取整算子，IuGPS为由GPS插值结果反演的解缠干涉图，Iw为缠绕干涉图，λ为波长，Iu为初始解缠干涉图，I′u为解缠干涉图，r为比例系数。

[0215] 其中，残差点识别和标记策略为：计算2×2大小的封闭路径单元累积值S，若S不为零则将所有4个像元点都标记成残差点，从而完成残差点的识别和标记。残差点标记完成后利用基于马尔科夫随机场GPS辅助InSAR解缠算法对非残差点进行解缠。该步解缠结束后，将所有已解缠点标记为固定域再进行一次只采用平滑约束的基于马尔科夫随机场GPS辅助InSAR解缠算法，从而完成对整个干涉图的解缠。可见考虑残差点的基于马尔科夫随机场GPS辅助InSAR解缠算法首先对质量好的像元点进行解缠，再以质量好的像元点的解缠值为约束进一步对质量不好的点进行解缠，从而最大限度的抑制误差传播。考虑残差点的基于马尔科夫随机场GPS辅助InSAR解缠算法实现了解缠精度与范围的较好综合。考虑残差点的基于马尔科夫随机场GPS辅助InSAR解缠算法的数据处理流程如图5所示。

[0216] 11、将InSAR解缠相位转换为形变信息

[0217]

[0218] 其中，□Rd为视线向的形变量，φd为解缠相位。

[0219] 12、InSAR形变信息地理编码

[0220] 本方法采用距离多普勒(R-D)定位模型，将前面处理得到的雷达坐标系下的结果转化到地理坐标系下的过程。在R-D定位模型中，SAR影像上任意像元的位置由三个方程确定：距离方程、多普勒质心平面方程和地球模型方程。

[0221] 距离方程：

[0222] 多普勒方程：

[0223] 地球模型方程：

[0224] 式中，R为SAR传感器到地面点间的距离， 分别为卫星和地面点的位置矢量，“·”为矢量点乘，fD为多普勒频移，λ为雷达波长，为卫星速度矢量，为地面点速度矢量，当采用WGS-84坐标系统时，的数值为零，xt，yt和zt为地面点位置矢量 在x，y和z轴上的投影，a为长半轴，b为短半轴，h为平均正常高，a＝6378.137km，b＝6356.752km。通过对上述三个方程的求解即可完成测量成果的地理编码过程。

[0225] 13、采用基于改进马尔科夫随机场的GPS和InSAR融合方法，融合InSAR和GPS的形变信息获取高时空分辨率的地表三维形变监测结果

[0226] 首先，将InSAR视线向的一维形变结果通过GPS约束分解成三维形变信息。具体算法如下：

[0227] (1)建立视线向的一维形变结果与三维形变结果之间的关系

[0228] dlos＝sinθcos(αh-3π/2)dn+sinθsin(αh-3π/2)de-cosθdu

[0229] (28)

[0230] ＝[un，ue，uu][dn，de，du]T

[0231] 其中，dlos为InSAR视线向形变测量，[un，ue，uu]为视线向单位矢量，[un，ue，uu]＝[0.34，-0.095，0.935]，dn、de和du为北东天方向优化值，θ为雷达视角，其值随像元点至卫星星下点的距离而改变；αh为卫星飞行方向与北方向夹角，其值可以从SAR数据中读出。

[0232] (2)确定总能量函数

[0233]

[0234]

[0235] 其中，

[0236]

[0237]i ni ei ui i

[0238] 其中，Wins 为InSAR约束权值，Wkrig ，Wkrig ，Wkrig 为GPS约束权值，dlos 为InSARi i i i i i i视线向形变测量，[un，ue，uu]为视线向单位矢量，dn、de 和du 为北东天方向优化值，σinsni ei ui ni ei ui
为InSAR观测误差，dkrig 、dkrig 和dkrig 为北东天方向克里格插值，σkrig 、σkrig 和σkrig为北东天方向克里格插值标准差。

[0239] (3)采用交叉验证定权法确定约束权值，具体如下：

[0240] ●对GPS与InSAR观测值进行配准，并将GPS三维观测值换算到InSAR视线方向。

[0241] ●按下式估计GPS观测站对应InSAR像元测量方差：

[0242]los

[0243] 式中，n为GPS观测站个数，dGPS，i 为第i个GPS观测站三维观测值换算到InSARlos视线向结果，dins，i 为InSAR第i个像元视线向观测值。

[0244] ●InSAR每一像元点的测量方差可由下式计算：

[0245]

[0246] 式中，σins，i2和γi分别为InSAR第i个像元点测量方差估值和相干系数，E(γGPS)为GPS观测站对应InSAR像元相干系数均值。

[0247] ●将GPS观测值分成A、B两组，对A组进行普通克里格插值。

[0248] ●求B组GPS三维观测值与其对应插值的方差：

[0249]

[0250] 式中，n′为B组GPS观测站个数，dGPS，in、dGPS，ie和dGPS，iu为北东天方向B组第i个GPS观测站观测值，dkrig，in、dkrig，ie和dkrig，iu为北东天方向B组第i个GPS观测站克里格插值。

[0251] ●计算克里格插值方差改正因子：

[0252]

[0253] 式中，E[(σkrig，n2)B]、E[(σkrig，e2)B]和E[(σkrig，u2)B]分别为B组GPS观测站北东天方向克里格插值方差均值。

[0254] ●对克里格方差进行改正：

[0255]

[0256] 式中，(σkrig，in)2、(σkrig，ie)2和(σkrig，iu)2为第i个像元点北东天方向克里格插值方差。

[0257] ●定权：

[0258]

[0259]

[0260] (4)采用直接解析法获取三维形变场的最大后验估计

[0261]

[0262] 然后，对高时间频率GPS数据(其所采集的地面沉降量)在时间域上对上述高空间分辨率的三维形变场进行克里格插值和加密，从而将准GPS结果在时间域内加密成准GPS时间序列，获取高时间分辨率的三维形变场。

[0263] 最后利用卡尔曼滤波对格网中的所有点进行估计获得高时空分辨率的地表三维形变结果。