一种适用于大梯度地表沉降监测的时序差分雷达干涉方法

摘要

本申请公开了一种适用于大梯度地表沉降监测的时序差分雷达干涉方法，该方法通过短时间基线差分干涉图筛选、离散点相位解缠、基于短时间基线差分干涉图的形变分量建模和解算、形变分量可靠性检验、差分干涉图相位梯度修正、对上述过程进行迭代以确保形变分量解算正确，以及基于修正后差分干涉图的形变时间序列建模和解算等过程，并形成整体的技术方案，解决原有时序差分雷达干涉技术中由于形变和相位梯度较大导致形变相位模糊度解算困难或失败的问题，最终达到正确提取大梯度地表形变速率和形变时间序列的目的，并起到降低大梯度形变建模和解算所需合成孔径雷达影像数量的效果，节约时序差分雷达干涉的应用经济成本。

说明书

技术领域

本发明属于地表沉降监测技术领域，具体地说，涉及一种适用于大梯度地表沉降监测的时序差分雷达干涉方法。

背景技术

地表沉降(垂直向地表形变，为便于表述，后续采用“形变”代表“沉降”)是发生范围最广的地质灾害之一，具有持续时间长的特点，且多发于城市及其近郊等经济发达和人口聚集区，对经济发展、城市建设和人民生活均会构成持久危害。我国是世界上地表沉降灾害最为严重的国家之一，累积沉降大于200毫米的面积超过15万平方公里，主要集中在华北平原、长江三角洲和汾渭盆地等经济发达地区，且出现了严重的沉降漏斗，造成了严重的经济损失。对地表沉降(尤其是沉降漏斗)开展大范围精确监测，对沉降防控及避免相应的危害具有十分重要的现实意义。

目前，时序差分合成孔径雷达(Synthetic Aperture Radar,SAR)干涉(Differential SAR Interferometry,DInSAR)技术已在地表形变监测中得到广泛应用，且是微波遥感、卫星大地测量以及地球物理学领域研究和应用的热点之一。时序DInSAR(TimeSeries DInSAR,TS-DInSAR)具有覆盖范围广、空间分辨率高、效率高、精度高且不易受云雾和阴雨天气影响的技术优势，非常适用于开展大范围地表形变监测。TS-DInSAR对缓慢累积性地表形变具有较好的监测效果和精度，但当形变较快和形变梯度较大时易造成解算精度较低或解算失败(如形变欠估计和形变漏斗区的“空洞”现象，即结果缺失)。

发明内容

有鉴于此，本申请针对TS-DInSAR难以满足快速或大梯度地表形变监测需求的问题，提供了一种适用于大梯度地表沉降监测的时序差分雷达干涉方法。

为了解决上述技术问题，本申请公开了一种适用于大梯度地表沉降监测的时序差分雷达干涉方法，具体按照以下步骤实施：

步骤1，对所有筛选后的SAR影像进行任意干涉组合并计算时间基线(形成一个干涉对的两幅SAR影像的获取时间差)和空间基线(两次成像时刻SAR传感器之间的空间距离，一般取该距离垂直于SAR视线方向的分量)；

步骤2，设置时空基线阈值进行干涉(干涉组合)对初始筛选，在保证干涉对数量的前提下限制时间失相干和空间失相干，按照干涉组合进行差分干涉处理得到差分干涉图集；

步骤3，永久散射体(Persistent Scatterer,PS)和相干目标(Coherent Target,CT)探测及点集合并，得到相干散射体(Coherent Scatterer,CS)点集，提取CS点集上的差分干涉相位，并以CS为节点构建不规则三角网(Triangular Irregular Network,TIN)，本发明技术方案中采用Delaunay三角网；

步骤4，基于Delaunay三角网和最小费用流(Minimum Cost Flow,MCF)方法的CS相位解缠；

步骤5，线性形变速率和高程误差建模及解算；

步骤6，基于数据模拟或地面测量数据的线性形变解算结果检验，若结果不可靠则执行步骤7，否则转向步骤8；

步骤7，重新筛选参与计算的差分干涉图，并重新执行步骤5和步骤6，当达到计算终止条件时，执行步骤8；

步骤8，从所有原始差分干涉图中减去线性形变相位分量，对残差相位(从差分干涉图中减去线性形变相位后剩余的相位)进行重新解缠，然后再将线性形变相位分量加回重新解缠后的相位中，重新执行步骤5，得到更新后的线性形变速率和高程误差(这一步主要是让更多干涉对参与计算，提高线性形变速率和高程误差的解算精度)；

步骤9，记录步骤8中所得到的线性形变速率，并从所有的原始差分干涉图中减去更新后的线性形变和高程误差相位分量，重新执行基于离散点的相位解缠，然后将步骤8中的线性形变分量重新加回新的解缠相位中，执行形变时间序列建模和解算过程，最终得到形变时间序列。

与现有技术相比，本申请可以获得包括以下技术效果：

本发明技术方案依据TS-DInSAR时序差分干涉图中形变相位大小以及形变相位梯度大小与时间间隔(差分干涉图的时间基线，即形成差分干涉图的两幅SAR影像获取时间差)成正相关关系这一特点，提出短时间基线差分干涉图筛选、离散点相位解缠、基于短时间基线差分干涉图的形变分量建模和解算、形变分量可靠性检验、差分干涉图相位梯度修正、对上述过程进行迭代以确保形变分量解算正确，以及基于修正后差分干涉图的形变时间序列建模和解算这一整体的技术方案，解决原有TS-DInSAR技术中由于形变和相位梯度较大导致形变相位模糊度解算困难或失败的问题，最终达到正确提取大梯度地表形变速率和形变时间序列的目的。

采用本发明技术方案，只要差分干涉图序列中存在至少6个短时间基线干涉对(由4幅SAR影像构成)即可实现线性形变分量的解算，然后开展长时间基线差分干涉图的形变相位梯度修正，促使更多的差分干涉图得到正确解缠，并参与形变分量的重新估算以及形变时间序列的建模和解算。因此，无需很高的SAR影像使用量便可实现大梯度形变的有效提取。

当然，实施本申请的任一产品并不一定需要同时达到以上所述的所有技术效果。

附图说明

此处所说明的附图用来提供对本申请的进一步理解，构成本申请的一部分，本申请的示意性实施例及其说明用于解释本申请，并不构成对本申请的不当限定。在附图中：

图1是一种适用于大梯度地表沉降监测的时序差分雷达干涉方法的实施流程图；

图2是35幅SAR影像任意组合所形成的干涉对的时间基线和空间基线分布图；

图3是空间基线阈值为30米时干涉对的时空基线分布图；

图4是第一次解算所得CS点上线性形变速率结果图；

图5是使用第一次解算所得的线性形变模拟的差分干涉图与原始差分干涉图的对比；其中，a为原始差分干涉图，b为使用第一次解算所得的线性形变模拟的差分干涉图；

图6是经3次迭代后解算所得CS点上的线性形变速率结果图；

图7是使用第三次解算所得的线性形变模拟的差分干涉图与原始差分干涉图的对比；其中，a为原始差分干涉图，b为使用第三次解算所得的线性形变模拟的差分干涉图；

图8是2009年6月23日至2010年7月2日期间的累积形变量；

图9是图8中所标示三个特征点(A、B和C)的形变时间序列。

具体实施方式

以下将配合附图及实施例来详细说明本申请的实施方式，藉此对本申请如何应用技术手段来解决技术问题并达成技术功效的实现过程能充分理解并据以实施。

一种适用于大梯度地表沉降监测的时序差分雷达干涉方法，如附图1所示，具体按照以下步骤实施：

步骤1，对所有筛选后的SAR影像进行任意干涉组合并计算时间基线和空间基线；

在筛选出合适的SAR影像(排除受雨雪等天气以及积雪影响的SAR影像)后，进行任意干涉组合配对，假设有N+1幅SAR影像，通过任意干涉组合可形成N(N+1)/2个干涉对。每个干涉对由主、从两幅SAR影像构成。然后，根据每个干涉对主、从SAR影像的获取时间计算该干涉对的时间基线(即SAR影像获取的时间差)，根据主、从SAR影像的参数文件所记录的SAR传感器运行位置及相应的时间参数计算该干涉对的空间基线(主、从影像成像时SAR传感器的空间距离，实际中一般取该空间距离在垂直于SAR传感器视线方向上的分量为空间基线，也即垂直基线)。

本实施例中采用覆盖天津市精武镇的35幅SAR影像为数据进行展示。附图2所示为35幅SAR影像进行任意组合所形成的干涉对的时间基线和空间基线分布。图2中以“年年年年月月日日”的方式标注了每幅SAR影像的获取时间，如“20090327”代表2009年3月27日。通过虚线相连的两幅SAR影像为一个干涉对，纵轴代表干涉对的空间基线大小，横轴代表时间，时间的差值即为时间基线。

步骤2，设置时空基线阈值进行干涉对的初始筛选，在保证干涉对数量的前提下限制时间失相干和空间失相干，按照干涉组合进行差分干涉处理得到差分干涉图集；

在进行时序差分干涉处理时，为降低时间失相干和空间失相干的影响，常采用时间基线阈值和空间基线阈值方法排除时间基线和空间基线大于某给定阈值的干涉对(受时间失相干和空间失相干影响较显著的干涉对)。大量研究表明，时间基线不宜超过2年，空间基线不宜超过150米。但是，当SAR影像获取时间整体跨度较小时(如第一幅SAR影像和最后一幅SAR影像的获取时间间隔小于2年)可不考虑对时间基线进行限制。在本实施例中，考虑到SAR影像获取时间的整体跨度较短(1年零9个月)，故未设置时间基线阈值，仅考虑对空间基线进行限制，设置空间基线阈值为30米，排除空间基线大于30米的干涉对。附图3所示为排除空间基线大于30米的干涉对后，得到保留的干涉对的时间基线和空间基线分布。基于此干涉组合，进行差分干涉处理，得到相应的差分干涉图集，处理过程如下：

以其中一个干涉对为例，设组成此干涉对的两幅SAR影像分别为I₁和I₂，分别记为主影像和从影像，对于影像中所记录的一点P(即影像像元P)，其在I₁和I₂中所对应的同名像元的值分别为：

Ψ₁＝A₁exp(jφ₁)(1)

Ψ₂＝A₂exp(jφ₂)(2)

其中，exp(·)代表以e为底的指数函数；Ψ₁为目标P在I₁中的像元值，φ₁为目标P在I₁中的相位值；Ψ₂为目标P在I₂中的像元值，φ₂为目标P在I₂中的相位值；A₁为目标P在I₁中信号的振幅值；A₂为目标P在I₂中信号的振幅值。其中，Ψ₁和Ψ₂为复数。由复数的共轭相乘可得P在由I₁和I₂形成的干涉图中的像元值为：

Ψ_int，P＝Ψ₁Ψ₂^*＝A₁A₂exp(j(φ₁-φ₂))>

其中，“*”代表复数的共轭，即Ψ₂^*为Ψ₂的共轭复数。

由式(3)可得P点所对应的干涉相位(即两次观测的相位差表达形式)为：

φ_int,P＝φ₁-φ₂>

对SAR影像中的所有像元进行相同的处理，即可得到由主、从SAR影像I₁和I₂所形成的干涉图。以P点为例，此时，P点的干涉相位可表达为：

φ_int,P＝φ_ref,P+φ_top,P+φ_def,P+φ_atm,P+φ_add,P>

其中，φ_int,P为P点在由SAR影像I₁和I₂所形成的干涉图中的干涉相位；φ_ref,P为P点的参考椭球面相位，由SAR传感器两次对同一地区观测时所处位置不同以及地球曲面引起；φ_top,P为P点的地形相位，由地表真实高程及其变化引起；φ_def,P为P点所发生的地表形变对应的相位；φ_atm,P为P点上的大气延迟相位，由雷达波在大气层中传播发生折射引起；φ_add,P是P点上的附加相位，主要包括各种随机噪声以及失相干噪声(两次成像时地面目标发生变化，导致雷达回波信号产生差异所致)。

所谓差分干涉处理，即在上述干涉处理基础上，去除地球参考椭球面相位和地形相位，得到以形变相位为主的差分干涉相位，形成差分干涉图。首先，在不考虑大气延迟和噪声相位的情况下，P点上的干涉相位可表达为：

其中，等式右边第一项即为参考椭球面相位的表达式，第二项为地形相位的表达式，第三项为形变相位的表达式；B为干涉对的空间基线；θ为雷达侧视角；α为空间基线与水平方向的夹角；λ为雷达波长；R₁为雷达传感器到观测目标之间的斜距；为空间基线B在垂直于雷达侧视方向上的投影分量，即垂直基线；h_P为P点的高程；Δr为P点所发生的地表形变。需要说明的是，Bsin(θ-α)即为干涉对的平行基线(空间基线B在沿雷达侧视方向上的投影分量)且有：

其中，R₁和R₂分别为雷达传感器两次观测成像时到地面点P的斜距；ΔR＝R₁-R₂为两次成像雷达传感器到P点斜距的差。

差分干涉处理即要从干涉相位中减去公式(6)等号右边的前两项。首先需要计算并扣除参考椭球面相位。由公式(6)等号右边第一项以及公式(7)可知，参考椭球面相位与卫星两次观测至参考椭球面间斜距的差ΔR有关，即：

由公式(7)可知，要获得ΔR的值，需首先已知R₁和R₂，此二者可通过两点间距离公式计算：

其中，d(·)为两点间距离方程；M和S分别为对应于主、从影像中P点成像时刻的卫星三维空间坐标矢量，由SAR影像参数文件提供的参数计算得到；P为P点在参考椭球面上的空间三维坐标矢量，可通过卫星成像参数、斜距方程、多普勒方程和椭球方程计算得到。此时，可通过式(8)和(9)得到P点上的参考椭球面相位。

对于地形相位，通过外部数字高程模型(Digital Elevation Model,DEM)可获取地面点P的高程值h_P(美国地质调查局网站公布了全球免费的DEM，本实施例中即采用该DEM)，R₁为主影像斜距，可由P点在主、从影像上的成像参数计算得到。对于垂直基线有：

其中，空间基线B＝d(M,S)，即主、从影像成像时雷达传感器的空间距离；在得到h_P、R₁、和后，可通过式(6)等号右边第二项计算P点的地形相位。

在得到P点的参考椭球面相位和地形相位后，从式(5)中扣除这两部分相位分量即可得到差分干涉相位。对SAR影像中所有的像元点进行上述差分干涉处理，即可得到相应的差分干涉图。

步骤3，PS和CT探测及点集合并，得到CS点集，提取CS点集上的差分干涉相位，并以CS为节点构建不规则三角网(Triangular Irregular Network,TIN)，本发明技术方案中采用Delaunay三角网；

PS探测采用振幅阈值和振幅离差指数(Amplitude Dispersion Index,ADI)阈值双阈值方法。对于SAR影像中某一特定的像元，其在N+1幅SAR影像中的振幅值可直接从SAR振幅影像中提取，则其时序振幅均值和ADI值为：

其中，i为SAR影像索引号(序号)；a_i和分别为该像元在第i幅SAR影像中的振幅值及在所有影像中的振幅均值；D_a和σ_a分别为该像元的ADI值及时序振幅标准差。当该像元的振幅信息满足如下两个条件时，认为其为PS：

其中，和σ_A分别为N+1幅SAR影像所有像元时序振幅值的均值和标准差；c和l分别为像元的列号和行号，C和L分别为影像列数和行数；A_i和a_cl分别为第i幅SAR影像所有像元振幅值的均值和行、列号分别为c和l的某像元的振幅值。公式(13)中0.25即为ADI阈值。

CT探测采用相干系数阈值方法。假设由N+1幅SAR影像形成了L个干涉对，并通过差分干涉数据处理得到L幅差分干涉图。通过下式计算每个像元在所有干涉图中相应的相干系数：

其中，为某像元在第l幅干涉图中的相干系数值；IM_l和IS_l分别为第l个干涉对的主影像和从影像；Z为相干系数估计窗口内像元数目，z为窗口内像元索引；l为干涉图索引。当某个像元满足以下条件时被识别为CT：

其中，min(·)表示取变量的最小值；γ_crit为判别某个像元是否为CT的相干系数阈值，一般可取0.25至0.3，本实施例中采用0.25。

当探测出所有的PS和CT后，将二者进行合并，并去除重复点，得到CS点集，最后从所有差分干涉图中提取CS点上的差分干涉相位值。

步骤4，基于Delaunay三角网和MCF方法的CS相位解缠；

相位解缠一般是基于规则格网进行，常规的相位解缠方法针对差分干涉图进行处理，易受失相干区域的影响，而基于离散点的相位解缠方法可有效避免这一问题。本发明技术方案中采用基于离散点的相位解缠方法对CS上的差分干涉相位进行解缠处理。对于某个CS点P，其在差分干涉图中的相位和该点上的真实相位的关系为：

其中，为P点上的真实相位(即相位解缠要得到的相位值)；φ_P为P点上的缠绕相位值(即P在差分干涉图中对应的相位值)，位于[-π,π)区间内，只记录了不足整周(2π)的小数部分，存在整周模糊度；n_P为整周模糊数。

首先，根据Delaunay剖分法则，以所有的CS为顶点构造Delaunay三角网，然后以网络边端点对应的两个CS之间的缠绕相位差为观测量，估算两点间的绝对相位差。通过最小费用流(Minimum Cost Flow,MCF)方法对与相位不连续性相关的网络费用流进行估算，并寻找最小费用流对应的积分路径，进行相位积分，完成相位解缠。此过程也即求解公式(16)中n_P的过程。

步骤5，线性形变速率和高程误差建模及解算；

在对CS上的差分干涉相位进行相位解缠后，根据步骤4中构建的Delaunay三角网重新计算相邻CS点间的解缠相位差。此处，相邻的CS点是指Delaunay三角网中三角形边的端点。以任意一边的两端点P和Q为例，二者在第l个差分干涉图中对应于的解缠相位为：

其中，和分别为第l幅差分干涉图的时间基线和垂直基线；λ为雷达波波长；θ_P和θ_Q分别为P和Q点处的雷达波入射角；R_P和R_Q分别为SAR传感器到P和Q之间的斜距；和分别为P和Q点上解缠后的差分干涉相位；v_P和v_Q分别为P和Q点的线性形变速率；δh_P和δh_Q分别为P和Q点上的高程误差；和分别为P和Q点在第l幅差分干涉图中的非线性形变相位；和分别为P和Q点在第l幅差分干涉图中的轨道误差相位；和分别为P和Q点在第l幅差分干涉图中的大气延迟相位；和分别为P和Q点在第l幅差分干涉图中的噪声相位。

线性形变和高程误差建模及解算的目的是对v_P和v_Q及δh_P和δh_Q进行估算。在本发明技术方案中，以P和Q上的解缠相位为观测量进行网络邻域差分建模。对于P和Q，二者之间的网络邻域差分相位增量(解缠相位差)为：

其中，和分别为P和Q点处斜距的平均值和入射角的平均值；为P和Q之间的邻域差分相位增量；Δv_PQ为P和Q之间的线性形变速率增量(即线性形变速率差)；Δδh_PQ为P和Q之间的高程误差增量(即高程误差之差)；为第l幅差分干涉图中P和Q之间的相位残差增量(即相位残差之差)，即P和Q点上非线性形变、轨道误差、大气延迟和噪声相位之间的差值之和。

以L幅差分干涉图为例，对于P和Q所对应的网络边而言，可列出L个与式(19)相同的方程，组成相应的线性方程组。将其表达为矩阵的形式有：

ΔΨ＝AX+W(20)

其中，

A＝[κ,η](24)

其中，

此时，方程组中仅有Δv_PQ和Δδh_PQ两个未知数，观测值数量为L，通过最小二乘可得未知参数的估值为：

对于所有通过网络边连接的CS点，均通过该过程解算相邻两点间的线性形变速率增量和高程误差增量。在得到所有CS点间的线性形变速率增量和高程误差增量后，以此作为观测量，以每个CS点的线性形变速率和高程误差为待估参数，在空间域进行最小二乘网络平差计算，解算所有CS点上的线性形变速率和高程误差。对于线性形变速率有：

其中，和分别为点P和Q的线性形变速率估计值(待估参数)；为点P和Q之间的线性形变速率增量估值；r_PQ为最小二乘残差。对于所有的CS点(设数量为G)和网络边(设数量为H)而言，根据公式(28)得到相应的观测方程，表达为矩阵的形式为：

其中，B为系数矩阵，其元素由弧段所对应的方程式确定，由1、-1和0组成；矩阵中为待估参数，即每个CS点的线性形变速率；中为所有CS点线性形变速率增量的估值；r为残差矩阵。由最小二乘平差可得：

其中，P为线性速率增量的权阵，可通过弧段最小二乘残差进行估算得到。

如附图4所示为使用附图3中所展示的干涉对所解算得出的CS点上的线性形变速率。图4中以分级设色的方式表达线性形变速率的量值，LOS DR表示雷达视线向形变速率，mm/yr为形变速率单位，即毫米/年。在步骤6中，将对此结果进行检验。

步骤6，基于数据模拟或地面测量数据的线性形变解算结果检验，若结果不可靠则执行步骤7，否则转向步骤8；

一般地，对TS-DInSAR形变结果进行验证主要有两种途径：基于外部测量数据的验证方法和基于差分干涉图模拟的验证方法。由于外部数据往往不容易获得，因此本发明技术方案中采用基于差分干涉图模拟的验证方法。

该方法是使用估计所得到的年形变速率或累积形变量和相应的干涉基线参数模拟差分干涉图。即将形变量转换为形变相位，形变到相位的转换公式为：

其中，为线性形变在第l幅差分干涉图中对应的相位；v为通过步骤5中的数据处理所得的线性形变速率。

在大形变梯度区域(如沉降漏斗区)，差分干涉图中的形变相位梯度随时间基线的增大而显著提升，在差分干涉图中表现为密集的条纹变化，若估计所得形变量符合或接近真实情况，使用形变结果模拟的差分干涉图与原始差分干涉图应该相似，尤其是形变条纹数目应相等或差别较小。若估计所得形变量不符合真实情况(如欠估计)，则模拟所得条纹数目与原始差分干涉图会存在明显区别。据此，可简单有效地对形变解算结果的有效性进行验证，且这种验证方法能够有效识别形变欠估计现象。

如附图5所示为由步骤5解算所得线性形变模拟所得的差分干涉图与原始差分干涉图的对比。其中，左侧(图5a)为原始差分干涉图，右侧(图5b)为模拟的差分干涉图。二者条纹数目差别明显，模拟所得差分干涉图中条纹数目低于原始差分干涉图，表明发生了形变欠估计现象。因此，第一次解算的线性形变并不可靠。此时，需要执行步骤7。

步骤7，重新筛选参与计算的差分干涉图，并重新执行步骤5和步骤6，当达到计算终止条件时，执行步骤8；

该步骤是实现正确解算线性形变的关键。经过步骤6的检验，确定所解算结果不可靠时，则通过设置更为严格的时间基线重新筛选出短时间基线差分干涉图，重新执行步骤5和步骤6，并进行检验。具体地，重新筛选短时间基线差分干涉图是指设置更小的时间基线，筛选出比上一次解算所用干涉图时间基线更短的干涉图。这是由于在一般的地表沉降区，时间基线越短则形变越小，且形变梯度越小，差分干涉图中的形变条纹数目就越少，更容易对差分干涉图进行正确的相位解缠，为形变的建模和解算提供正确的观测量。例如，本实施例中第一次解算时未对时间基线进行限制，解算所得到的线性形变不可靠，则设置较短的时间基线值(本实施例中采用180天)，筛选出时间基线小于180天的干涉图重新解算线性形变，然后再次进行结果检验，若不可靠，则进行第三次解算，本实施例中第三次解算时选取90天作为时间基线阈值，即筛选出时间基线小于90天的干涉图重新解算线性形变。时间基线逐渐缩短，即所谓的设置更为严格的时间基线限制。

本实施例中，经过3次迭代计算可获取正确的线性形变结果，如附图6所示。附图7所示为利用第三次解算所得线性形变模拟的差分干涉图与原始差分干涉图的对比，其中左侧(图7a)为原始差分干涉图，右侧(图7b)为模拟的差分干涉图。此时，二者条纹数目几乎一致。得到可靠的线性形变后，可执行步骤8。

需要指出的是，在本实施例中，通过3次迭代计算获取了正确的线性形变速率，但对于不同的研究区域，迭代次数是可变的，如增加或减少，以能够获取正确的线性形变为准。

步骤8，从所有原始差分干涉图中减去线性形变相位分量，对残差相位(从差分干涉图中减去线性形变相位后剩余的相位)进行重新解缠，然后再将线性形变相位分量加回重新解缠后的相位中，重新执行步骤5，得到最终的线性形变速率和高程误差(这一步主要是让更多干涉对参与计算，提高线性形变速率和高程误差的解算精度)；

在通过第7步骤得到正确的线性形变后，从所有的原始差分干涉图中减去线性形变所对应的相位分量。为便于阐述，设此时所得到的正确的线性形变为仍以第l幅差分干涉图为例，其所对应的相位分量为：

此时，通过下式获取扣除线性形变分量后的相位：

其中，为扣除线性形变分量后的缠绕相位；φ_diff,l为原始的差分干涉相位；conj(·)为求复数共轭。对所有L幅差分干涉图执行此步运算，得L幅修正时空相位梯度后的差分干涉图。此时，差分干涉图中的相位梯度将显著降低，条纹数目减少，有利于相位解缠的正确执行。

在得到L幅经过时空相位梯度修正的差分干涉图后，重新对干涉图集执行步骤4中所述的相位解缠过程，得到L幅相应的解缠相位图，新的解缠相位表达为：

其中，为扣除线性形变分量后的解缠相位；为高程误差相位分量；为非线性形变相位分量；为轨道误差相位分量；为大气延迟相位分量；为噪声相位分量。在得到修正后差分干涉相位的解缠值后，将步骤7中所得线性形变分量(即所对应的形变相位)与其相加，得到修正后的解缠相位：

在得到修正的解缠相位后，以其为观测值重新执行步骤5中的解算流程，得到最终的线性形变速率和高程误差。为便于后续阐述，此处设最终的线性形变速率和高程误差分别为和则二者所对应的相位为：

现有研究表明，参与运算的干涉对数目越多，线性形变速率和高程误差的解算精度越高，这也是对相位梯度进行修正使更多干涉对参与运算的重要目的之一。此外，在后续的数据处理中，最终的线性形变和高程误差将用于最终的相位梯度修正。

步骤9，记录步骤8中所得到的线性形变速率，并从所有的原始差分干涉图中减去更新后的线性形变和高程误差相位分量，重新执行基于离散点的相位解缠，然后将步骤8中的线性形变分量重新加回新的解缠相位中，执行形变时间序列建模和解算过程，最终得到形变时间序列。

在得到最终的线性形变速率和高程误差结果后，首先保留并输出线性形变速率结果，然后可根据下式将线性形变和高程误差相位(可通过公式(36)和(37)计算)从原始差分干涉相位中扣除，得到新的差分干涉相位值：

其中，为扣除线性形变和高程误差相位后的差分干涉相位，且有：

在得到最终经线性形变和高程误差修正的相位后，重新执行第4步骤中的相位解缠过程。在得到新的解缠相位后，重新将步骤8中得到的线性形变相位加回，得到扣除高程误差后的解缠相位：

其中，为经高程误差校正后的解缠相位；为由公式(36)计算得到的线性形变相位。

对于某个CS而言，其在第l幅扣除了高程误差的解缠相位图中的相位表达为：

其中，为形变相位(包括线性形变和非线性形变)；为轨道误差相位；为大气延迟相位；为噪声相位。对于L个差分干涉图，经高程误差改正后的相位为：

仍以上述N+1幅SAR影像为例，设其获取时刻为T＝[t₀,t₁,…,t_n,…,t_N]，由其所形成的L个干涉对的主(IM)、从(IS)SAR影像集分别为：

IM＝[IM₁,IM₂,…,IM_l,…,IM_L]；IS＝[IS₁,IS₂,…,IS_l,…,IS_L]>

此时，以t₀时刻为参考时刻，其余任意时刻t_n(n＝1,2,…,N)相对于t₀时刻的相位为待估参数，以扣除高程误差后的解缠相位为观测值，建立观测方程并恢复每个时刻对应的相位(即相位时间序列)。为得到符合实际物理意义的形变估计参数，假设相邻两幅SAR影像获取时间间隔内相位的变化符合线性累积趋势，将对相位时间序列的求解转化为对相位变化速率(相邻时间段内的平均相位变化速率)v_ph(注意与步骤5、6、7和8中形变速率v的区别)的求解，此时，待估参数为：

其中，为t_n时刻相对于t₀时刻的累积相位，且有式(44)实际的物理意义是相邻两幅SAR影像获取时间间隔内的相位变化速率，也称为分段相位速率。相应地，该模型可被称为分段线性模型。此处的“段”指相邻两幅SAR影像之间的时间段。此时可得观测方程为：

将式(45)表达为矩阵的形式有：

其中，B为L×N的矩阵，矩阵元素B[l,n]＝t_n-t_n-1(IS_l+1≤n≤IM_l)，其它元素值为零。由于干涉图集可能是不连续的，即在某个SAR影像获取时刻断开，进而造成系数矩阵B发生奇异(此时方程组为欠定状态，不存在唯一解)。因此，实际处理中首先对B进行奇异值分解(Singular Value Decomposition,SVD)处理，然后估算出最小范数意义下各时间段的平均相位速率，通过在时间域上进行积分即可恢复与SAR影像获取时刻相对应的相位时间序列v_ph。B的奇异值分解为：

B＝USW^T>

其中，U为L×L阶正交矩阵，被称作B的左奇异向量，其前N列是BB^T的特征向量；W是N×L阶正交矩阵，被称作B的右奇异向量。S是半正定L×L阶对角矩阵，其元素为相应于BB^T特征向量的均方根，也即奇异值σ_i，且有，

S＝diag(σ₁,…,σ_N-r+1,0,…,0)(48)

其中，diag(·)代表对角矩阵，除对角线元素为σ_i外，其余元素值均为零；N-r+1为矩阵B的秩，r代表差分干涉图子集的数量(即由于设置了时间基线和空间基线阈值限制，导致由N+1幅SAR影像所形成的干涉图集合被分成了多个子集，造成干涉图集不连续)。最终，可对公式(46)进行解算得：

其中，S⁺＝diag(1/σ₁,…,1/σ_N-r+1,…,0,…,0)。

在解算出相邻时间段内的平均相位变化速率后，根据相位变化速率与时间之积并求和得到对应于SAR影像获取时刻的相位序列，即的值。然后，从相位序列中扣除由步骤8所得到的线性形变分量，并在空间域进行低通滤波去除噪声相位，在时间域进行高通滤波得到大气延迟和轨道误差相位序列。最后，从相位序列中扣除大气延迟和轨道误差相位序列即可得到地表形变所对应的相位序列此时可得形变时间序列为：

其中，D为与每幅SAR影像获取时刻相对应的累积形变量所组成的向量；为与每幅SAR影像获取时刻相对应的累积形变相位。即：

如附图8所示为2009年6月23日至2010年7月2日期间的累积形变量。图8中标记了A、B和C三个特征点，附图9所示为三个特征点的形变时间序列。至此，便完成了大梯度形变区域(沉降漏斗区域)形变速率和形变时间序列的解算。

本发明技术方案带来的有益效果：

(1)解决原有TS-DInSAR技术难以进行大梯度地表形变建模和解算的问题

原有TS-DInSAR技术难以适用于大梯度地表形变的建模和解算，本发明技术方案依据TS-DInSAR时序差分干涉图中形变相位大小以及形变相位梯度大小与时间间隔(差分干涉图的时间基线，即形成差分干涉图的两幅SAR影像获取时间差)成正相关关系这一特点，提出短时间基线差分干涉图筛选、离散点相位解缠、基于短时间基线差分干涉图的形变分量建模和解算、形变分量可靠性检验、差分干涉图相位梯度修正、对上述过程进行迭代以确保形变分量解算正确，以及基于修正后差分干涉图的形变时间序列建模和解算这一整体的技术方案，解决原有TS-DInSAR技术中由于形变和相位梯度较大导致形变相位模糊度解算困难或失败的问题，最终达到正确提取大梯度地表形变速率和形变时间序列的目的。

针对本发明技术方案中具体的思路进行有益效果分析如下：

(a)短时间基线差分干涉图筛选。短时间基线差分干涉图中形变相位梯度较小，较容易实现正确的相位模糊度解算(即相位解缠)，得到正确的解缠相位后，可用于后续形变分量建模和解算。此外，短时间基线差分干涉图筛选是一个可以重复执行的过程。

(b)离散点相位解缠。相位解缠也即解算相位模糊度，离散点相位解缠仅针对高质量的观测目标进行处理，避免低质量目标的影响。解缠后的相位将被用于形变分量建模。相比于原有技术基于非解缠相位进行形变分量解算而言，基于解缠相位进行形变分量建模可避免对SAR影像时间采样率的依赖。

(c)基于短时间基线差分干涉图的形变分量建模和解算，并检验形变分量的可靠性。在(a)和(b)基础上进行形变分量(线性形变速率)建模和解算(方法见技术方案的详细阐述部分)，并采用差分干涉图模拟方法检验形变分量解算结果的可靠性。先保证形变分量解算正确，后续将被应用于对差分干涉图的形变相位梯度进行修正，即从差分干涉图中减去形变分量，起到降低形变相位梯度的作用，进而实现所有差分干涉图的正确解缠。

(d)迭代计算。在对形变分量进行可靠性检验时，若发现结果不可靠，则重新进行(a)、(b)和(c)，确保形变分量解算结果可靠。迭代计算是连接各个步骤的关键。

(e)采用正确的形变分量结果重新对差分干涉图进行相位梯度修正，然后重新进行相位解缠，并基于此对形变时间序列进行建模和解算。

综上所述，以上各个过程紧密结合，共同形成了本发明的整体技术方案，最终实现大梯度形变区域形变分量和形变时间序列的精确解算。

(2)可降低SAR影像使用量，降低TS-DInSAR技术的应用经济成本

在大梯度形变解算方面，现有研究指出SAR时间采样率的提高可有效缓解梯度较大造成的问题，但这往往受现有SAR系统观测周期(对同一地区进行重复观测的周期)以及研究区域可用影像数量的限制，并且数据成本会显著提高(如针对大梯度形变，以往的TS-DInSAR技术往往需要很高的SAR影像时间采样率，即较高的数据量需求)。尤其是高分辨率SAR影像，价格较高，数据成本将十分显著。

采用本发明技术方案，理论上只要差分干涉图序列中存在至少6个短时间基线干涉对(由4幅SAR影像构成)即可实现线性形变分量的解算，然后开展长时间基线差分干涉图的形变相位梯度修正，促使更多的差分干涉图(尤其是长时间基线差分干涉图)得到正确解缠，并参与形变分量的重新估算以及形变时间序列的建模和解算。因此，无需很高的SAR影像时间采样率便可实现大梯度形变的有效提取。

上述说明示出并描述了发明的若干优选实施例，但如前所述，应当理解发明并非局限于本文所披露的形式，不应看作是对其他实施例的排除，而可用于各种其他组合、修改和环境，并能够在本文所述发明构想范围内，通过上述教导或相关领域的技术或知识进行改动。而本领域人员所进行的改动和变化不脱离发明的精神和范围，则都应在发明所附权利要求的保护范围内。