CN102608584B - 基于多项式反演模型的时间序列InSAR形变监测方法及装置 - Google Patents

Abstract

本发明提供一种基于多项式反演模型的时间序列InSAR形变监测方法及装置，该方法包括：将某一地区N幅SAR单视复影像组合生成M幅干涉图，并生成M幅差分相位图；计算平均相干系数图，提取高相干点，对相邻的两个高相干点的差分相位进行再次差分建立多项式反演模型，求解相邻点的相对多项式形变与相对高程误差，以某一具有已知形变量和DEM误差的高相干点为参考点，分别集成得到每个高相干点的多项式形变和高程误差；得到高相干点的多项式反演模型相位，从高相干点的差分相位中减去该高相干点的多项式反演模型相位，得到残差相位；从残差相位中提取残余形变，与多项式形变叠加得到高相干点的地表形变信息。本发明为高精度的地表形变监测提供了一种解决方案。

Description

基于多项式反演模型的时间序列InSAR形变监测方法及装置

技术领域

本发明涉及合成孔径雷达技术，尤其涉及一种基于多项式反演模型的时间序列合成孔径雷达干涉测量(InSAR，Synthetic Aperture Radar Interferometry)形变监测方法及装置。

背景技术

合成孔径雷达(SAR，Synthetic Aperture Radar)是20世纪50年代发展起来的最重要的对地观测技术，它通过雷达天线在随载体的运动中以一定的时间间隔发射电磁脉冲信号，在不同位置上接收地面物体反射的回波信号，并记录和存储下来，形成地面的高分辨率图像。与可见光、近红外传统遥感技术相比，SAR具有全天候全天时成像能力，不受天气和时间影响，如微波可以穿透云层和一定程度上穿透雨区，可以不依赖于太阳作为照射源进行夜间成像，这是其他遥感手段所不具备的。

SAR技术只能获取地表目标物的二维信息，缺乏获取目标点高程信息和监测目标微小形变的能力。将干涉测量技术与SAR技术结合而形成的合成孔径雷达干涉测量技术(InSAR，Synthetic Aperture Radar Interferometry)提供了获取地面三维信息的全新方法，它通过两副天线同时观测或通过一副天线两次平行观测，获取地面同一景观的复图像对，根据地面各点在两幅复图像中的相位差，得出各点在两次成像中微波的路程差，从而获得地面目标的高程信息。基于InSAR技术的发展，合成孔径雷达差分干涉测量(DifferentialSynthetic Aperture Radar Interferometry，DInSAR)技术是对两幅以上的干涉图或对一幅干涉图加一幅地面数字高程模型(DEM，Digital Elevation Model)进行再处理的一种技术，它通过去除地形引起的干涉相位，获得关于地表形变的信息，在火山监测、地震位移测量、地面沉降等领域具有极大的应用潜力。

基于重复轨道的DInSAR技术容易受空间失相干、时间失相干和大气干扰等因素的影响，难以进行常态化的实际应用。引起空间失相干的原因包括大的垂直基线和大的形变梯度。时间失相干则是由于两次成像时刻环境的不一致导致同一地面像元中散射体的散射特性发生改变而引起的。易变的大气条件可能会导致在两幅干涉影像上不一致的相位延迟，从而引起形变测量误差。

为了克服传统DInSAR技术的这些限制，自上世纪90年代末，一些新的InSAR处理技术被提出。这些技术的共同特点是：基于时间序列SAR影像进行处理，处理的对象不是整幅影像的全部像元，而是其中具有稳定散射特性从而能在较长时间间隔内保持高相干的像元子集，也就是高相干点。这些技术总体上可以概括为两类：以永久散射体干涉(Permanent Scatterer或者Persistent Scatterer Interferometry，或PS-InSAR)为代表的单一主影像时间序列InSAR技术和以小基线集技术(Small baseline subset interferometry，或SBASInSAR)为代表的多主影像时间序列InSAR技术。为叙述方便，将这两种技术统称为时间序列InSAR技术。时间序列InSAR技术对上述三种制约因素均有良好的免疫力，目前已经取代传统的DInSAR技术在火山、地震、滑坡、地面沉降等领域得到大量应用。由于城市地区拥有密集的天然点目标，如路灯、屋顶等，因此时间序列InSAR技术在城市地面沉降测量方面应用最为广泛。

但是，时间序列InSAR技术中仍有一些问题有待完善，其中之一就是关于形变模型。现有的时间序列InSAR处理中，都将形变模式假设为以线性形变为主，然后在后续处理中再恢复出残余的非线性形变分量。然而，当高相干点的形变呈现高度的非线性且在空间不相关的时候，形变反演结果会有很大误差。本发明针对现有时间序列InSAR处理中的线性形变模型的不足，提出了基于多项式的形变反演模型，及相应的时间序列InSAR形变监测方法。

发明内容

本发明实施例提供一种基于多项式反演模型的时间序列InSAR形变监测方法及装置，以解决如何利用时间序列InSAR进行高精度地表形变监测的问题。

一方面，本发明实施例提供了一种基于多项式反演模型的时间序列InSAR形变监测方法，所述时间序列InSAR形变监测方法包括：将某一地区N幅SAR单视复影像根据小基线原则组合成M个干涉像对，生成M幅干涉图，并去除所述M幅干涉图中由数字高程模型DEM模拟得到的地形相位图，生成M幅差分相位图；通过所述M幅干涉图计算平均相干系数图，通过预设相干系数阈值从该平均相干系数图中提取高相干点，并对所述M幅差分相位图的相邻高相干点的差分相位进行再次差分，得到相邻高相干点的二次差分相位；通过所述相邻高相干点的二次差分相位建立多项式反演模型，求解两个高相干点间的相对多项式形变与相对高程误差，以某一具有已知形变量和DEM误差的高相干点为参考点，分别集成两两高相干点间的相对多项式形变与相对高程误差，得到每个高相干点上的多项式形变和高程误差；通过所述多项式形变和高程误差得到高相干点的多项式反演模型相位，从高相干点的差分相位中减去该高相干点的多项式反演模型相位，得到残差相位；从所述残差相位中提取残余形变，将所述残余形变与所述多项式形变叠加得到所述高相干点的地表形变信息。

可选的，在本发明一实施例中，所述小基线原则为对时间基线和空间基线的限制。

可选的，在本发明一实施例中，所述高相干点的差分相位包括：多项式形变相位、高程误差相位、残余形变相位、大气影响相位、噪声相位。

可选的，在本发明一实施例中，对所述M幅差分相位图的差分相位进行再次差分，得到相邻高相干点的二次差分相位，包括：通过德劳内Delaunay三角网连接各高相干点，对所述德劳内Delaunay三角网上边长小于或等于2公里的相邻点作邻域差分，得到相邻高相干点的二次差分相位；其中，所述相邻高相干点为空间距离小于或等于2公里的两个高相干点。

可选的，在本发明一实施例中，通过所述多项式形变和所述高程误差得到所述高相干点的多项式反演模型相位，该多项式反演模型相位包括如下参数：一次形变速率、二次形变速率、三次形变速率和高程误差。

另一方面，本发明实施例提供了一种基于多项式反演模型的时间序列InSAR形变监测装置，所述时间序列InSAR形变监测装置包括：

差分相位图生成单元，用于将某一地区N幅SAR单视复影像根据小基线原则组合成M个干涉像对，生成M幅干涉图，并去除所述M幅干涉图中由数字高程模型DEM模拟得到的地形相位图，生成M幅差分相位图；

相邻高相干点的二次差分相位获取单元，用于通过所述M幅干涉图计算平均相干系数图，通过预设相干系数阈值从该平均相干系数图中提取高相干点，并对所述M幅差分相位图的相邻高相干点的差分相位进行再次差分，得到相邻高相干点的二次差分相位；

多项式反演模型单元，用于通过所述相邻高相干点的二次差分相位建立多项式反演模型，求解两个高相干点间的相对多项式形变与相对高程误差，以某一具有已知形变量和DEM误差的高相干点为参考点，分别集成两两高相干点间的相对多项式形变与相对高程误差，得到每个高相干点上的多项式形变和高程误差；

地表形变信息获取单元，用于通过所述多项式形变和所述高程误差得到高相干点的多项式反演模型相位，从高相干点的差分相位中减去该高相干点的多项式反演模型相位，得到残差相位；并从所述残差相位中提取残余形变，将所述残余形变与所述多项式形变叠加得到所述高相干点的地表形变信息。

可选的，在本发明一实施例中，所述小基线原则为对时间基线和空间基线的限制。

可选的，在本发明一实施例中，所述高相干点的差分相位包括：多项式形变相位、高程误差相位、残余形变相位、大气影响相位、噪声相位。

可选的，在本发明一实施例中，所述相邻高相干点的二次差分相位获取单元，具体用于通过德劳内Delaunay三角网连接各高相干点，对所述德劳内Delaunay三角网上的边长小于或等于2公里的相邻点作邻域差分，得到相邻高相干点的二次差分相位；其中，所述相邻高相干点为空间距离小于或等于2公里的两个高相干点。

可选的，在本发明一实施例中，所述地表形变信息获取单元，具体用于通过所述多项式形变和所述高程误差得到所述高相干点的多项式反演模型相位，该多项式反演模型相位包括如下参数：一次形变速率、二次形变速率、三次形变速率和高程误差。

上述技术方案具有如下有益效果：与现有基于线性反演模型的时间序列InSAR方法相比，上述本发明实施例为高精度的地表形变监测提供了一种解决方案，能更好地拟合实际地表形变。

附图说明

为了更清楚地说明本发明实施例或现有技术中的技术方案，下面将对实施例或现有技术描述中所需要使用的附图作简单地介绍，显而易见地，下面描述中的附图仅仅是本发明的一些实施例，对于本领域普通技术人员来讲，在不付出创造性劳动性的前提下，还可以根据这些附图获得其他的附图。

图1为本发明实施例基于多项式反演模型的时间序列InSAR形变监测方法流程图；

图2为本发明实施例基于多项式反演模型的时间序列InSAR形变监测装置结构示意图；

图3是本发明应用实例通过预设相干系数阈值提取的高相干点示意图；

图4是本发明应用实例高相干点生成的Delaunay三角网示意图；

图5是本发明应用实例用于精度验证的30个水准点分布图；

图6是本发明应用实例得到的太原地区一次形变速率示意图；

图7是本发明应用实例得到的太原地区二次形变速率示意图；

图8是本发明应用实例得到的太原地区三次形变速率示意图；

图9是本发明应用实例得到的太原地区平均形变速率示意图；

图10是现有技术基于线性反演模型InSAR技术得到的太原地区平均形变速率示意图；

图11是本发明应用实例太原地区2004-2006年多项式反演模型得到的形变量和现有技术线性反演模型得到的形变量与同时段水准数据的精度比较表。

具体实施方式

下面将结合本发明实施例中的附图，对本发明实施例中的技术方案进行清楚、完整地描述，显然，所描述的实施例仅仅是本发明一部分实施例，而不是全部的实施例。基于本发明中的实施例，本领域普通技术人员在没有做出创造性劳动前提下所获得的所有其他实施例，都属于本发明保护的范围。

如图1所示，为本发明实施例提供了基于多项式反演模型的时间序列InSAR形变监测方法流程图，所述时间序列InSAR形变监测方法包括：

101、将某一地区N幅SAR单视复影像(SLC，Single Look Complex)根据小基线原则组合成M个干涉像对，生成M幅干涉图，并去除所述M幅干涉图中由数字高程模型DEM模拟得到的地形相位图，生成M幅差分相位图。

小基线原则可以为对时间基线的限制和空间基线的限制两方面，对于C波段的ERS-1/2SAR或ENVISAT ASAR数据而言，时间基线一般不超过3年，空间基线一般不超过400米。基于N幅SAR SLC影像，可组合形成M个小基线干涉像对，生成M幅干涉图，去除由该地区外部DEM模拟得到的地形相位图，可得到M幅差分相位图，各像素的相位模型可表示为

φ_{dif_i}＝φ_{def_poly}+φ_{residual_def}+φ_{errortopo_i}+φ_{atm_i}+φ_{noise_i}    (1)

式中，i为第i幅干涉图，φ_{def_poly}为多项式形变相位，用于拟合真实地表形变，φ_{residual_def}为残余形变相位，φ_errortopo为DEM引起的高程误差相位，φ_atm为大气影响相位，φ_noise为噪声相位，包括热噪声、时间和基线去相干噪声。其中多项式形变相位φ_{def_poly}和高程误差相位φ_errortopo可分别表示为：

${\mathrm{\φ}}_{\mathrm{def}\_\mathrm{poly}}=\frac{4\mathrm{\π}}{\mathrm{\λ}}\·[v_1\·(t_m-t_s)+v_2\·(t_m^2-t_s^2)+v_3\·(t_m^3-t_s^3)]---\left(2\right)$

${\mathrm{\φ}}_{\mathrm{errortopo}}=\frac{4\mathrm{\π}}{\mathrm{\λ}}\·\frac{b}{r\·\sin \mathrm{\θ}}\·\mathrm{\ξ}---\left(3\right)$

式中，λ为雷达载波波长，r为目标到雷达传感器斜距长，b为垂直基线，θ为雷达入射角，ξ为高程误差，v₁、v₂、v₃分别为一次形变速率、二次形变速率、三次形变速率，t_m、t_s表示干涉像对的主、辅影像获取的时间。

102、通过所述M幅干涉图计算平均相干系数图，通过预设相干系数阈值从该平均相干系数图中提取高相干点，并对所述M幅差分相位图的相邻高相干点的差分相位进行再次差分，得到相邻高相干点的二次差分相位。

其中，所述相邻高相干点为空间距离小于或等于2公里的两个高相干点。

计算每幅干涉图中各像元的相干系数，并得到其平均相干系数，选择适当的预设相干系数阈值，将平均相干系数大于该预设相干系数阈值的像元认为是高相干点，如下式：

$\frac{1}{M}\underset{i}{\mathrm{\Σ}}{r}_i\≥\mathrm{thd}\_r---\left(4\right)$

其中r_i表示单个像素在第i幅干涉图上的相干系数，thd_r是预设相干系数阈值。

通过Delaunay三角网连接各高相干点，对三角网上边长小于或等于2公里的相邻点作邻域差分，其两顶点(x_m，y_m)，(x_n，y_n)间的二次差分相位可表示为：

${\mathrm{\δ\φ}}_{\mathrm{dif}}(x_m,y_m,x_n,y_n,T_i)=\frac{4\mathrm{\π}}{\mathrm{\λ}}\·[(v_1(x_m,y_m)-v_1(x_n,y_n))\·(t_m-t_s)+(v_2(x_m,y_m)-v_2(x_n,y_n))\·(t_m^2-t_s^2)+$

$(v_3(x_m,y_m)-v_3(x_n,y_n))\·(t_m^3-t_s^3)]+\frac{4\mathrm{\π}}{\mathrm{\λ}}\·\frac{b\left(T_i\right)}{r\left(T_i\right)\sin \left({\mathrm{\θ}}_i\right)}\·[\mathrm{\ξ}(x_m,y_m)-\mathrm{\ξ}(x_n,y_n)]---\left(5\right)$

$+[\mathrm{\β}(x_m,y_m)-\mathrm{\β}(x_n,y_n)]+[\mathrm{\α}(x_m,y_m)-\mathrm{\α}(x_n,y_n)]+[n(x_m,y_m)-n(x_n,y_n)]$

式中，(x_m，y_m)、(x_n，y_n)为该条边两顶点的位置坐标，T_i为第i幅干涉图的时间基线，T_i＝t_m-t_s，β为残余形变相位，α为大气影响相位，n为噪声相位。

103、通过所述相邻高相干点的二次差分相位建立多项式反演模型，求解两个高相干点间的相对多项式形变与相对高程误差，以某一具有已知形变量和DEM误差的高相干点为参考点，分别集成两两高相干点间的相对多项式形变与相对高程误差，得到每个高相干点上的多项式形变和高程误差。

考虑到大气影响相位是一个低频信号，在空间上存在一个相关距离，一般约小于或等于2km，比如为0.5至2km，当三角网两点间距在大气影响相关距离内时，可以认为它们的大气影响相位相等，即α(x_m，y_m，T_i)≈α(x_n，y_n，T_i)，由于一次、二次、三次形变速率和高程误差是常数，而残余形变相位和噪声相位都是随机信号，因此，通过相邻高相干点间的相对多项式形变和相对高程误差可以建立所述多项式反演模型，即(5)式可以演变为：

${\mathrm{\δ\φ}}_{\mathrm{mode}l}(x_m,y_m,x_n,y_n,T_i)=\frac{4\mathrm{\π}}{\mathrm{\λ}}\·[{\mathrm{\Δv}}_1(x_m,y_m,x_n,y_n)\·(t_m-t_s)$

$+{\mathrm{\Δv}}_2(x_m,y_m,x_n,y_n)\·(t_m^2-t_s^2)+{\mathrm{\Δv}}_3(x_m,y_m,x_n,y_n)\·(t_m^3-t_s^3)]$

$+\frac{4\mathrm{\π}}{\mathrm{\λ}}\·\frac{b\left(T_i\right)}{r\left(T_i\right)\·\sin \left({\mathrm{\θ}}_i\right)}\·\mathrm{\Δ\ξ}(x_m,y_m,x_n,y_n)---\left(6\right)$

其中，Δv₁、Δv₂、Δv₃和Δξ分别是两高相干点间的相对一次形变速率、相对二次形变速率、相对三次形变速率和相对高程误差。为了求得这四个参数量，可以建立如下的目标函数方程，通过最优化方法求解：

${\mathrm{\γ}}_{\mathrm{mode}l}(x_m,y_m,x_n,y_n)=\frac{1}{M}\·\left|\underset{i=0}{\overset{M}{\mathrm{\Σ}}}\exp \right[j\·({\mathrm{\δ\φ}}_{\mathrm{dif}}(x_m,y_m,x_n,y_n,T_i)-{\mathrm{\δ\φ}}_{\mathrm{mode}l}(x_m,y_m,x_n,y_n,T_i))\left]\right|---\left(7\right)$

式中，γ为整体相位相干系数，其最大值在[0，1]之间，反映了相对形变速率(一次、二次、三次)和相对高程误差对二次差分相位的拟合程度；M为干涉像对个数；δφ_dif为相邻高相干点的二次差分相位；δφ_model为相邻高相干点的多项式反演模型相位；(x_m，y_m)、(x_n，y_n)为两相邻高相干点的位置坐标；T_i为第i幅差分相位图的时间基线，T_i＝t_m-t_s。当对所有的边完成最大化求解后，例如可以选择γ≥0.7的边作为可靠的连接，以某一具有已知形变量和DEM误差的高相干点为参考点，分别集成两两高相干点间的相对多项式形变速率Δv₁、Δv₂、Δv₃与相对高程误差Δξ，得到每个高相干点上的多项式形变v₁、v₂、v₃和高程误差ξ。

104、通过所述多项式形变和高程误差得到高相干点的多项式反演模型相位，从高相干点的差分相位中减去该高相干点的多项式反演模型相位，得到残差相位。

由两高相干点(x_m，y_m)，(x_n，y_n)之间的所述多项式反演模型相位(6)，相应地，可建立所述高相干点上的多项式反演模型，表示为：

${\mathrm{\φ}}_{\mathrm{mode}l}\left(T_i\right)=\frac{4\mathrm{\π}}{\mathrm{\λ}}\·[v_1\·(t_m-t_s)+v_2\·(t_m^2-t_s^2)+v_3\·(t_m^3-t_s^3)]+\frac{4\mathrm{\π}}{\mathrm{\λ}}\·\frac{b\left(T_i\right)}{r\left(T_i\right)\·\sin \left({\mathrm{\θ}}_i\right)}\·\mathrm{\ξ}---\left(8\right)$

通过所述多项式形变v₁、v₂、v₃和所述高程误差ξ计算得到所述高相干点的多项式反演模型相位φ_model，并从差分相位中减去该高相干点的多项式反演模型相位，得到高相干点的残差相位φ_res，它包括大气影响相位φ_atm，残余形变相位φ_non-linear以及噪声相位φ_noise。

105、从所述残差相位中提取残余形变，将所述残余形变与所述多项式形变叠加得到所述高相干点的地表形变信息。

为了得到高相干点的完整的地表形变信息，首先对残差相位进行相位解缠，然后对解缠后的残差相位进行时空频谱特征分析，以提取残余形变相位。残差相位的三个分量中，大气影响相位在时间序列上是不相关的，为高频信号，在空间分布是相关的，为低频信号；残余形变相位在时间序列上是低频信号，在空间上是不相关的，为高频信号；噪声相位则是时间和空间都不相关的随机高频信号。利用这些表现特征，可以将三者分离出来。对所述高相干点，首先在时间序列上做频域低通滤波，提取出低频的残余形变相位，然后利用最小二乘方法及干涉组合关系计算出各时刻的残余形变，将其与多项式形变叠加可得到所述高相干点的完整地表形变信息。

对应于上述方法实施例，如图2所示，为本发明实施例基于多项式反演模型的时间序列InSAR形变监测装置结构示意图，所述时间序列InSAR形变监测装置包括：

差分相位图生成单元21，用于将某一地区N幅SAR单视复影像根据小基线原则组合成M个干涉像对，生成M幅干涉图，并去除所述M幅干涉图中由数字高程模型DEM模拟得到的地形相位图，生成M幅差分相位图；

相邻高相干点的二次差分相位获取单元22，用于通过所述M幅干涉图计算平均相干系数图，通过预设相干系数阈值从该平均相干系数图中提取高相干点，并对所述M幅差分相位图的相邻高相干点的差分相位进行再次差分，得到相邻高相干点的二次差分相位；

多项式反演模型单元23，用于通过所述相邻高相干点的二次差分相位建立多项式反演模型，求解两个高相干点间的相对多项式形变与相对高程误差，以某一具有已知形变量和DEM误差的高相干点为参考点，分别集成两两高相干点间的所述相对多项式形变与所述相对高程误差，得到每个高相干点上的多项式形变和高程误差；

地表形变信息获取单元24，用于通过所述多项式形变和所述高程误差得到高相干点的多项式反演模型相位，从高相干点的差分相位中减去该高相干点的多项式反演模型相位，得到残差相位；并从所述残差相位中提取残余形变，将所述残余形变与所述多项式形变叠加得到所述高相干点的地表形变信息。

可选的，所述小基线原则可以为对时间基线和空间基线的限制。所述高相干点的差分相位可以包括：多项式形变相位、高程误差相位、残余形变相位、大气影响相位、噪声相位。

可选的，所述相邻高相干点的二次差分相位获取单元22，具体可以用于通过德劳内Delaunay三角网连接各高相干点，对所述德劳内Delaunay三角网上的边长为小于或等于2公里的两个高相干点作邻域差分，得到相邻高相干点的二次差分相位；其中，所述相邻高相干点为空间距离小于或等于2公里的两个高相干点。所述地表形变信息获取单元24，具体可以用于通过所述多项式形变和所述高程误差得到所述高相干点的多项式反演模型相位，该多项式反演模型相位包括如下参数：一次形变速率、二次形变速率、三次形变速率和高程误差。

本发明上述方法或装置实施例技术方案具有如下有益效果：与现有基于线性反演模型的InSAR方法相比，上述本发明实施例为高精度的地表形变监测提供了一种解决方案，能更好地拟合实际地表形变。

为了更好地说明本发明实施例技术方案的有效性和优越性，现对应用上述技术方案对本发明实施例与现有线性模型反演技术进行如下对比分析：如图3所示，是本发明应用实例通过预设相干系数阈值提取的高相干点示意图，预设相干系数阈值设为0.3，提取了9940个高相干点。如图4所示，为本发明应用实例高相干点生成的Delaunay三角网示意图。如图5所示，是本发明应用实例用于精度验证的30个水准点分布图，水准测量时段为2004～2006年。如图6所示，为本发明应用实例得到的太原地区一次形变速率示意图，单位mm/a。如图7所示，是本发明应用实例得到的太原地区二次形变速率示意图，单位mm/a²。如图8所示，是本发明应用实例得到的太原地区三次形变速率示意图，单位mm/a³。如图9所示，是本发明应用实例通过一次、二次、三次形变速率计算得到的太原地区平均形变速率示意图，单位mm/a。如图10所示，是现有技术基于线性反演模型InSAR技术得到的太原地区平均形变速率示意图，单位mm/a。为了比较本发明应用实例与现有技术线性反演模型的精度，利用太原地区2004-2006年的30个水准实测数据对二者同时段的形变量进行验证，如图11所示，是本发明应用实例太原地区2004-2006年多项式反演模型得到的形变量和现有技术线性反演模型得到的形变量与同时段水准数据的精度比较表，其数值单位为mm。本发明应用实例计算的误差标准差为3.17mm，现有技术计算的误差标准差则达6.11mm。结果表明，本发明应用实例的形变反演精度明显高于现有技术线性反演模型的结果，能更好地拟合实际地表形变。

本领域技术人员还可以了解到本发明实施例列出的各种说明性逻辑块(illustrativelogical block)，单元，和步骤可以通过电子硬件、电脑软件，或两者的结合进行实现。为清楚展示硬件和软件的可替换性(interchangeability)，上述的各种说明性部件(illustrativecomponents)，单元和步骤已经通用地描述了它们的功能。这样的功能是通过硬件还是软件来实现取决于特定的应用和整个系统的设计要求。本领域技术人员可以对于每种特定的应用，可以使用各种方法实现所述的功能，但这种实现不应被理解为超出本发明实施例保护的范围。

本发明实施例中所描述的各种说明性的逻辑块，或单元都可以通过通用处理器，数字信号处理器，专用集成电路(ASIC)，现场可编程门阵列(FPGA)或其它可编程逻辑装置，离散门或晶体管逻辑，离散硬件部件，或上述任何组合的设计来实现或操作所描述的功能。通用处理器可以为微处理器，可选地，该通用处理器也可以为任何传统的处理器、控制器、微控制器或状态机。处理器也可以通过计算装置的组合来实现，例如数字信号处理器和微处理器，多个微处理器，一个或多个微处理器联合一个数字信号处理器核，或任何其它类似的配置来实现。

本发明实施例中所描述的方法或算法的步骤可以直接嵌入硬件、处理器执行的软件模块、或者这两者的结合。软件模块可以存储于RAM存储器、闪存、ROM存储器、EPROM存储器、EEPROM存储器、寄存器、硬盘、可移动磁盘、CD-ROM或本领域中其它任意形式的存储媒介中。示例性地，存储媒介可以与处理器连接，以使得处理器可以从存储媒介中读取信息，并可以向存储媒介存写信息。可选地，存储媒介还可以集成到处理器中。处理器和存储媒介可以设置于ASIC中，ASIC可以设置于用户终端中。可选地，处理器和存储媒介也可以设置于用户终端中的不同的部件中。

在一个或多个示例性的设计中，本发明实施例所描述的上述功能可以在硬件、软件、固件或这三者的任意组合来实现。如果在软件中实现，这些功能可以存储与电脑可读的媒介上，或以一个或多个指令或代码形式传输于电脑可读的媒介上。电脑可读媒介包括电脑存储媒介和便于使得让电脑程序从一个地方转移到其它地方的通信媒介。存储媒介可以是任何通用或特殊电脑可以接入访问的可用媒体。例如，这样的电脑可读媒体可以包括但不限于RAM、ROM、EEPROM、CD-ROM或其它光盘存储、磁盘存储或其它磁性存储装置，或其它任何可以用于承载或存储以指令或数据结构和其它可被通用或特殊电脑、或通用或特殊处理器读取形式的程序代码的媒介。此外，任何连接都可以被适当地定义为电脑可读媒介，例如，如果软件是从一个网站站点、服务器或其它远程资源通过一个同轴电缆、光纤电脑、双绞线、数字用户线(DSL)或以例如红外、无线和微波等无线方式传输的也被包含在所定义的电脑可读媒介中。所述的碟片(disk)和磁盘(disc)包括压缩磁盘、镭射盘、光盘、DVD、软盘和蓝光光盘，磁盘通常以磁性复制数据，而碟片通常以激光进行光学复制数据。上述的组合也可以包含在电脑可读媒介中。

以上所述的具体实施方式，对本发明的目的、技术方案和有益效果进行了进一步详细说明，所应理解的是，以上所述仅为本发明的具体实施方式而已，并不用于限定本发明的保护范围，凡在本发明的精神和原则之内，所做的任何修改、等同替换、改进等，均应包含在本发明的保护范围之内。

Claims (10)

1.一种基于多项式反演模型的时间序列InSAR形变监测方法，其特征在于，所述时间序列InSAR形变监测方法包括：

将某一地区N幅SAR单视复影像根据小基线原则组合成M个干涉像对，生成M幅干涉图，并去除所述M幅干涉图中由数字高程模型DEM模拟得到的地形相位图，生成M幅差分相位图；

通过所述M幅干涉图计算平均相干系数图，通过预设相干系数阈值从该平均相干系数图中提取高相干点，并对所述M幅差分相位图的相邻高相干点的差分相位进行再次差分，得到相邻高相干点的二次差分相位：

${\mathrm{\δ\φ}}_{\mathrm{dif}}(x_m,y_m,x_n,y_n,T_i)=\frac{4\mathrm{\π}}{\mathrm{\λ}}\·[(v_1(x_m,y_m)-v_1(x_n,y_n))\·(t_m-t_s)+(v_2(x_m,y_m)-v_2(x_n,y_n))\·(t_m^2-t_s^2)+$

$(v_3(x_m,y_m)-v_3(x_n,y_n))\·(t_m^3-t_s^3)]+\frac{4\mathrm{\π}}{\mathrm{\λ}}\·\frac{b\left(T_i\right)}{r\left(T_i\right)\sin \left({\mathrm{\θ}}_i\right)}\·[\mathrm{\ξ}(x_m,y_m)-\mathrm{\ξ}(x_n,y_n)]$

$+[\mathrm{\β}(x_m,y_m)-\mathrm{\β}(x_n,y_n)]+[\mathrm{\α}(x_m,y_m)-\mathrm{\α}(x_n,y_n)]+[n(x_m,y_m)-n(x_n,y_n)]$

式中，(x_m,y_m)、(x_n,y_n)为相邻两顶点的位置坐标，T_i为第i幅干涉图的时间基线，T_i=t_m-t_s，β为残余形变相位，α为大气影响相位，n为噪声相位，λ为雷达载波波长，v₁、v₂、v₃分别为一次形变速率、二次形变速率、三次形变速率，t_m、t_s表示干涉像对的主、辅影像获取的时间，r为目标到雷达传感器斜距长，b为垂直基线，θ为雷达入射角，θ_i为第i幅干涉图的雷达入射角，ξ为高程误差；

通过所述相邻高相干点的二次差分相位建立多项式反演模型，求解两个高相干点间的相对多项式形变与相对高程误差，以某一具有已知形变量和DEM误差的高相干点为参考点，分别集成两两高相干点间的相对多项式形变与相对高程误差，得到每个高相干点上的多项式形变和高程误差；

通过所述多项式形变和高程误差得到高相干点的多项式反演模型相位，从高相干点的差分相位中减去该高相干点的多项式反演模型相位，得到残差相位；

从所述残差相位中提取残余形变，将所述残余形变与所述多项式形变叠加得到所述高相干点的地表形变信息。

2.如权利要求1所述时间序列InSAR形变监测方法，其特征在于，

所述小基线原则为对时间基线和空间基线的限制。

3.如权利要求1所述时间序列InSAR形变监测方法，其特征在于，

所述高相干点的差分相位包括：多项式形变相位、高程误差相位、残余形变相位、大气影响相位、噪声相位。

4.如权利要求1所述时间序列InSAR形变监测方法，其特征在于，对所述M幅差分相位图的差分相位进行再次差分，得到相邻高相干点的二次差分相位，包括：

通过德劳内Delaunay三角网连接各高相干点，对所述德劳内Delaunay三角网上边长小于或等于2公里的相邻点作邻域差分，得到相邻高相干点的二次差分相位；其中，所述相邻高相干点为空间距离小于或等于2公里的两个高相干点。

5.如权利要求1所述时间序列InSAR形变监测方法，其特征在于，通过所述多项式形变和所述高程误差得到所述高相干点的多项式反演模型相位，该多项式反演模型相位包括如下参数：一次形变速率、二次形变速率、三次形变速率和高程误差。

6.一种基于多项式反演模型的时间序列InSAR形变监测装置，其特征在于，所述时间序列InSAR形变监测装置包括：

差分相位图生成单元，用于将某一地区N幅SAR单视复影像根据小基线原则组合成M个干涉像对，生成M幅干涉图，并去除所述M幅干涉图中由数字高程模型DEM模拟得到的地形相位图，生成M幅差分相位图；

相邻高相干点的二次差分相位获取单元，用于通过所述M幅干涉图计算平均相干系数图，通过预设相干系数阈值从该平均相干系数图中提取高相干点，并对所述M幅差分相位图的相邻高相干点的差分相位进行再次差分，得到相邻高相干点的二次差分相位：

${\mathrm{\δ\φ}}_{\mathrm{dif}}(x_m,y_m,x_n,y_n,T_i)=\frac{4\mathrm{\π}}{\mathrm{\λ}}\·[(v_1(x_m,y_m)-v_1(x_n,y_n))\·(t_m-t_s)+(v_2(x_m,y_m)-v_2(x_n,y_n))\·(t_m^2-t_s^2)+$

$(v_3(x_m,y_m)-v_3(x_n,y_n))\·(t_m^3-t_s^3)]+\frac{4\mathrm{\π}}{\mathrm{\λ}}\·\frac{b\left(T_i\right)}{r\left(T_i\right)\sin \left({\mathrm{\θ}}_i\right)}\·[\mathrm{\ξ}(x_m,y_m)-\mathrm{\ξ}(x_n,y_n)]$

$+[\mathrm{\β}(x_m,y_m)-\mathrm{\β}(x_n,y_n)]+[\mathrm{\α}(x_m,y_m)-\mathrm{\α}(x_n,y_n)]+[n(x_m,y_m)-n(x_n,y_n)]$

式中，(x_m,y_m)、(x_n,y_n)为相邻两顶点的位置坐标，T_i为第i幅干涉图的时间基线，T_i=t_m-t_s，β为残余形变相位，α为大气影响相位，n为噪声相位，λ为雷达载波波长，v₁、v₂、v₃分别为一次形变速率、二次形变速率、三次形变速率，t_m、t_s表示干涉像对的主、辅影像获取的时间，r为目标到雷达传感器斜距长，b为垂直基线，θ为雷达入射角，θ_i为第i幅干涉图的雷达入射角，ξ为高程误差；

多项式反演模型单元，用于通过所述相邻高相干点的二次差分相位建立多项式反演模型，求解两个高相干点间的相对多项式形变与相对高程误差，以某一具有已知形变量和DEM误差的高相干点为参考点，分别集成两两高相干点间的相对多项式形变与相对高程误差，得到每个高相干点上的多项式形变和高程误差；

地表形变信息获取单元，用于通过所述多项式形变和所述高程误差得到高相干点的多项式反演模型相位，从高相干点的差分相位中减去该高相干点的多项式反演模型相位，得到残差相位；并从所述残差相位中提取残余形变，将所述残余形变与所述多项式形变叠加得到所述高相干点的地表形变信息。

7.如权利要求6所述时间序列InSAR形变监测装置，其特征在于，

所述小基线原则为对时间基线和空间基线的限制。

8.如权利要求6所述时间序列InSAR形变监测装置，其特征在于，

所述高相干点的差分相位包括：多项式形变相位、高程误差相位、残余形变相位、大气影响相位、噪声相位。

9.如权利要求6所述时间序列InSAR形变监测装置，其特征在于，

所述相邻高相干点的二次差分相位获取单元，具体用于通过德劳内Delaunay三角网连接各高相干点，对所述德劳内Delaunay三角网上的边长小于或等于2公里的相邻点作邻域差分，得到相邻高相干点的二次差分相位；其中，所述相邻高相干点为空间距离小于或等于2公里的两个高相干点。

10.如权利要求6所述时间序列InSAR形变监测装置，其特征在于，

所述地表形变信息获取单元，具体用于通过所述多项式形变和所述高程误差得到所述高相干点的多项式反演模型相位，该多项式反演模型相位包括如下参数：一次形变速率、二次形变速率、三次形变速率和高程误差。

Images