CN107544069B - 基于平面近似模型的多基线相位解缠绕方法 - Google Patents

Abstract

本发明属于雷达信号处理领域，公开了一种基于平面近似模型的多基线相位解缠绕方法，包括：获取干涉相位图，计算干涉相位图的缠绕相位梯度，然后进行基于平面近似的模糊数梯度估计，再通过MCF模型求解模糊数，根据模糊数与缠绕相位计算绝对相位，得到相位解缠绕结果；能够解决在噪声干扰下多基线相位解缠绕的噪声鲁棒性差的问题。

Description

基于平面近似模型的多基线相位解缠绕方法

技术领域

本发明属于雷达信号处理技术领域，尤其涉及一种基于平面近似模型的多基线相位解缠绕方法，适用于多基线干涉合成孔径雷达(InSAR)系统在配准、滤波与去平地之后的多基线相位解缠绕问题。

背景技术

作为传统SAR的发展，InSAR系统能够全天候、全天时工作，现已广泛应用于地形测绘、地表形变监测与自然灾害检测等诸多方面，因此，InSAR技术的发展一直受到了高度重视。

与单基线InSAR相比，多基线InSAR能够利用多幅干涉相位图，突破“相位连续性假设”，实现复杂地形(例如有峡谷、陡峭山崖等)的高精度测绘，这一优势使得多基线InSAR在军事及民用领域都具有很广泛的应用前景。

在多基线InSAR系统进行地形测绘中，多基线相位解缠绕是一个关键步骤，在传统的多基线相位解缠绕方法中，因其理论基础——中国余数定理本身对噪声的影响十分敏感，这直接导致了多基线InSAR在实际中的应用受到很大限制，如果能提高多基线相位解缠绕方法的噪声鲁棒性，将极大的推动复杂地形的高精度测绘以及相关领域的发展。

发明内容

针对上述问题，本发明的目的在于提供一种基于平面近似模型的多基线相位解缠绕方法，能够解决在噪声干扰下多基线相位解缠绕的噪声鲁棒性差的问题。

为达到上述目的，本发明采用如下技术方案予以实现。

一种基于平面近似模型的多基线相位解缠绕方法，所述方法包括如下步骤：

步骤1，从多基线InSAR系统中获取两个基线分别对应的干涉相位图，记为第一干涉相位图和第二干涉相位图，两幅干涉相位图的大小相同，且为m×n，其中，m、n分别为大于零的正整数；

步骤2，根据第一干涉相位图中每个像素的缠绕相位，确定第一干涉相位图中水平方向相邻像素间的水平缠绕相位梯度，记为水平缠绕第一相位梯度矩阵，其大小为m×(n-1)，并确定第一干涉相位图中垂直方向相邻像素间的垂直缠绕相位梯度，记为垂直缠绕第一相位梯度矩阵，其大小为(m-1)×n；

根据第二干涉相位图中每个像素的缠绕相位，确定第二干涉相位图中水平方向相邻像素间的水平缠绕相位梯度，记为水平缠绕第二相位梯度矩阵，其大小为m×(n-1)，并确定第二干涉相位图中垂直方向相邻像素间的垂直缠绕相位梯度，记为垂直缠绕第二相位梯度矩阵，其大小为(m-1)×n；

步骤3，设定预设窗口大小，根据预设窗口大小选定给定像素位置对应的估计窗口；

步骤4，给定像素位置(i,j)，1≤i≤m,1≤j≤n-1；

根据给定像素位置(i,j)，以及给定像素位置(i,j)对应的估计窗口，估计得到第一干涉相位图的水平方向模糊数梯度矩阵，其大小为m×(n-1)；以及估计得到第二干涉相位图的水平方向模糊数梯度矩阵，其大小为m×(n-1)；

步骤5，给定像素位置(i,j)，1≤i≤m-1,1≤j≤n；

根据给定像素位置(i,j)，以及给定像素位置(i,j)对应的估计窗口，估计得到第一干涉相位图的垂直方向模糊数梯度矩阵，其大小为(m-1)×n；以及估计得到第二干涉相位图的垂直方向模糊数梯度矩阵，其大小为(m-1)×n；

步骤6，将第一干涉相位图在水平方向的模糊数梯度矩阵中的元素以及干涉相位图1在垂直方向的模糊数梯度矩阵中的元素带入设定的优化模型，得到第一干涉相位图的模糊数矩阵1；

将第二干涉相位图在水平方向的模糊数梯度矩阵的元素以及第二干涉相位图在垂直方向的模糊数梯度矩阵的元素带入所述设定的优化模型，得到第二干涉相位图的模糊数矩阵2；

步骤7，给定像素位置(i，j)，1≤i≤m，1≤j≤n；

根据第一干涉相位图在像素位置(i，j)处对应的缠绕相位和第一干涉相位图在像素位置(i，j)处的模糊数，计算得到第一干涉相位图在像素位置(i，j)处的绝对相位；

根据第二干涉相位图在像素位置(i，j)处对应的缠绕相位和第二干涉相位图在像素位置(i，j)处的模糊数，计算得到第二干涉相位图在像素位置(i，j)处的绝对相位；

步骤8，令像素位置(i，j)遍历集合{(i，j)|1≤i≤m，1≤j≤n}中的所有元素，重复执行步骤7，从而得到第一干涉相位图的绝对相位矩阵，以及第二干涉相位图的绝对相位矩阵。

本发明的有益效果为：(1)本发明所提出的方法能够有效降低噪声对于模糊数梯度估计的影响，从而提高了多基线相位解缠绕方法的噪声鲁棒性；(2)本发明所提出的方法通过有效估计模糊数梯度，从而能够得到相对稀疏的残点分布，为后续最小费用流(MCF)模型的求解提供了良好的铺垫，也可以移植到其它相位解缠绕方法中。

附图说明

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

图1为本发明实施例提供的一种基于平面近似模型的多基线相位解缠绕方法的流程示意图；

图2为本发明方法与传统TSPA方法的处理结果示意图，所有图的横坐标对应相应矩阵的列，所有图的纵坐标对应相应矩阵的行；

图2中，图2(a)为参考绝对相位图；图2(b)为仿真干涉相位图一；图2(c)为仿真干涉相位图二；图2(d)为由TSPA方法求解得到的相位解缠绕结果图；图2(e)为由本发明方法求解得到的相位解缠绕结果图；图2(f)为图2(d)与2(a)的差值图；图2(g)为图2(e)与2(a)的差值图；

图3为本发明方法与传统TSPA方法的处理结果示意图，所有图的横坐标对应相应矩阵的列，所有图的纵坐标对应相应矩阵的行；

图3中，图3(a)为Google Earth图像图3(b)为短基线干涉相位图；图3(c)为长基线干涉相位图；图3(d)为TSPA方法的相位解缠绕结果图；图3(e)为本发明方法的相位解缠绕结果图；图3(f)为TSPA方法所得到的残点分布局部图；图3(g)为本发明方法得到的残点分布局部图。

具体实施方式

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

从多基线InSAR系统可以获得多幅干涉相位图，以下以基线数为2的情况为例说明本发明，对于基线数大于2的情况，可以同理推得。

本发明实施例提供一种基于平面近似模型的多基线相位解缠绕方法，如图1所示，所述方法包括如下步骤：

步骤1，从多基线InSAR系统中获取两个基线分别对应的干涉相位图，记为第一干涉相位图和第二干涉相位图，两幅干涉相位图的大小相同，且为m×n，其中，m、n分别为大于零的正整数。

设多基线InSAR系统中，有两个基线，即为基线1和基线2，基线1的基线长度为B₁，基线2的基线长度为B₂，具有两个基线的多基线InSAR系统可以获得两幅干涉相位图，基线1对应第一干涉相位图，基线2对应第二干涉相位图，两幅干涉相位图的尺寸均为m×n，干涉相位图中的元素称为缠绕相位，第一干涉相位图中位置(i,j)的缠绕相位记为

第二干涉相位图中位置(i,j)的缠绕相位记为

步骤2，根据第一干涉相位图中每个像素的缠绕相位，确定第一干涉相位图中水平方向相邻像素间的水平缠绕相位梯度，记为水平缠绕第一相位梯度矩阵，其大小为m×(n-1)，并确定第一干涉相位图中垂直方向相邻像素间的垂直缠绕相位梯度，记为垂直缠绕第一相位梯度矩阵，其大小为(m-1)×n；

根据第二干涉相位图中每个像素的缠绕相位，确定第二干涉相位图中水平方向相邻像素间的水平缠绕相位梯度，记为水平缠绕第二相位梯度矩阵，其大小为m×(n-1)，并确定第二干涉相位图中垂直方向相邻像素间的垂直缠绕相位梯度，记为垂直缠绕第二相位梯度矩阵，其大小为(m-1)×n。

步骤2具体包括如下子步骤：

(2a)获取第一干涉相位图中像素位置(i,j)处的缠绕相位

1≤i≤m,1≤j<n；获取第二干涉相位图中像素位置(i,j)处的缠绕相位

1≤i≤m,1≤j<n；像素位置(i,j)处的像素值为像素位置(i,j)处的缠绕相位；

(2b)记第一干涉相位图的水平缠绕第一相位梯度矩阵为

中第i行第j列元素

1≤i≤m，1≤j≤n-1；

记第二干涉相位图的水平缠绕第二相位梯度矩阵为

中第i行第j列元素

1≤i≤m，1≤j≤n-1；

(2c)记第一干涉相位图的垂直缠绕第一相位梯度矩阵为

中第i行第j列元素

1≤i≤m-1，1≤j≤n；

记第二干涉相位图的垂直缠绕第二相位梯度矩阵为

中第i行第j列元素

1≤i≤m-1，1≤j≤n。

步骤3，设定预设窗口大小，根据预设窗口大小选定给定像素位置对应的估计窗口。

步骤3中，设定预设窗口大小为(2p+1)×(2q+1)；

根据地形的先验信息确定正整数p和q作为预设窗口尺寸参数。预设窗口大小的选择(即参数p和q的大小)，通常需要根据对地形的先验信息进行选择，预设窗口选取的原则是：保证预设窗口内像素所对应的地形变化相对平缓，可以近似的认为是在同一个平面内。

(3a)对于水平缠绕第一相位梯度矩阵或者水平缠绕第二相位梯度矩阵中给定的像素位置(i,j)，1≤i≤m，1≤j≤n-1，其对应的估计窗口内所包含的水平缠绕相位梯度记为

或

其中，a₀≤i'≤a₁，b₀≤j'≤b₁，a₀，a₁，b₀，b₁的值由下式确定：

(3b)对于垂直缠绕第一相位梯度矩阵或者垂直缠绕第二相位梯度矩阵中给定的像素位置(i,j)1≤i≤m-1，1≤j≤n，其对应的估计窗口内所包含的垂直缠绕相位梯度记为

或

其中，c₀≤i″≤c₁，d₀≤j″≤d₁，c₀，c₁，d₀，d₁的值由下式确定：

步骤4，本步骤的目的是，利用窗内所有水平缠绕相位梯度，对给定位置的水平模糊数梯度进行估计。

给定像素位置(i,j)，1≤i≤m,1≤j≤n-1；

根据给定像素位置(i,j)，以及给定像素位置(i,j)对应的估计窗口，估计得到第一干涉相位图的水平方向模糊数梯度矩阵，其大小为m×(n-1)；以及估计得到第二干涉相位图的水平方向模糊数梯度矩阵，其大小为m×(n-1)。

步骤4具体包括：

(4a)对于给定像素位置(i,j)，第一干涉相位图的水平缠绕相位梯度

和第二干涉相位图的水平缠绕相位梯度

以及尺寸为(2p+1)×(2q+1)的估计窗口，记估计窗口内所有水平缠绕相位梯度位置(i',j')的集合

定义R_pq误差

其中，

其中，符号[·]表示取四舍五入，

为第一干涉相位图在像素位置(i,j)的水平方向模糊数梯度，

为第二干涉相位图在像素位置(i,j)的水平方向模糊数梯度，B₁表示基线1的长度，B₂表示基线2的长度，

和

为水平梯度修正因子，用于修正窗内像素位置(i,j)处的水平绝对相位梯度与窗内其它位置(i',j')处的水平绝对相位梯度因为相位缠绕所带来的梯度不一致；

(4b)求解如下优化模型：

其中，Z表示整数集合，第一干涉相位图在像素位置(i,j)的水平方向模糊数梯度估计值

和第二干涉相位图在像素位置(i,j)的水平方向模糊数梯度估计值

是该优化模型要求解的变量。

在实际应用中，可以根据关于地形的先验知识，确定

和

的取值范围，通过在可能的取值范围内遍历所有

和

的取值组合，计算相应的R_pq误差，经过比较即可获得满足优化模型的最优解。

步骤5，本步骤的目的是，利用窗内所有垂直缠绕相位梯度，对给定位置的垂直模糊数梯度进行估计。

给定像素位置(i,j)，1≤i≤m-1,1≤j≤n；

根据给定像素位置(i,j)，以及给定像素位置(i,j)对应的估计窗口，估计得到第一干涉相位图的垂直方向模糊数梯度矩阵，其大小为(m-1)×n；以及估计得到第二干涉相位图的垂直方向模糊数梯度矩阵，其大小为(m-1)×n。

步骤5具体包括：

(5a)对于给定像素位置(i,j)，第一干涉相位图的垂直缠绕相位梯度

和第二干涉相位图的垂直缠绕相位梯度

以及尺寸为(2p+1)×(2q+1)的估计窗口，记估计窗内所有垂直相位梯度位置(i″,j″)的集合

定义R_pq误差

其中，

其中，符号[·]表示取四舍五入，

为第一干涉相位图在像素位置(i,j)的垂直方向模糊数梯度，

为第二干涉相位图在像素位置(i,j)的垂直方向模糊数梯度，B₁表示基线1的长度，B₂表示基线2的长度，

和

为垂直梯度修正因子，用于修正窗内像素位置(i,j)处的垂直绝对相位梯度与窗内其它位置(i”,j”)处的垂直绝对相位梯度因为相位缠绕所带来的梯度不一致；

(5b)求解如下优化模型：

其中，Z表示整数集合，第一干涉相位图在像素位置(i,j)的垂直方向模糊数梯度估计值

和第二干涉相位图在像素位置(i,j)的垂直方向模糊数梯度估计值

是该优化模型要求解的变量。

在实际应用中，可以根据关于地形的先验知识，确定

和

的取值范围，通过在可能的取值范围内遍历所有

和

的取值组合，计算相应的R_pq误差，经过比较即可获得满足优化模型的最优解。

步骤6，将第一干涉相位图在水平方向的模糊数梯度矩阵中的元素以及干涉相位图1在垂直方向的模糊数梯度矩阵中的元素带入设定的优化模型，得到第一干涉相位图的模糊数矩阵1；

将第二干涉相位图在水平方向的模糊数梯度矩阵的元素以及第二干涉相位图在垂直方向的模糊数梯度矩阵的元素带入所述设定的优化模型，得到第二干涉相位图的模糊数矩阵2。

步骤6具体包括：

(6a)将第一干涉相位图在水平方向的模糊数梯度矩阵中的元素

以及第一干涉相位图在垂直方向的模糊数梯度矩阵中的元素

带入如下设定的优化模型：

其中，Z为整数集合，

和

表示设定的权值(它们可以利用InSAR图像的质量图生成，当不需要加权的时候，取

即可)，

和

为自由变量，k₁(i,j)、k₁(i+1,j)、k₁(i,j+1)分别为第一干涉相位图在像素位置(i,j)、(i+1,j)、(i,j+1)处的模糊数，它们是优化模型的输出。

该优化模型可以利用线性规划优化工具进行求解，求解结果记为模糊数矩阵1，对应于第一干涉相位图，该矩阵的尺寸为m×n。

(6b)将第二干涉相位图在水平方向的模糊数梯度矩阵中的元素

以及第二干涉相位图在垂直方向的模糊数梯度矩阵中的元素

带入如下设定的优化模型：

其中，Z为整数集合，

和

表示设定的权值(它们可以利用InSAR图像的质量图生成，当不需要加权的时候，取

即可)，

和

为自由变量，k₂(i,j)、k₂(i+1,j)、k₂(i,j+1)分别为第二干涉相位图在像素位置(i,j)、(i+1,j)、(i,j+1)处的模糊数它们是优化模型的输出。

该优化模型可以利用线性规划优化工具进行求解，求解结果为模糊数矩阵2，对应于第二干涉相位图，该矩阵的尺寸为m×n。

步骤7，给定像素位置(i,j)，1≤i≤m,1≤j≤n；

根据第一干涉相位图在像素位置(i,j)处对应的缠绕相位和第一干涉相位图在像素位置(i,j)处的模糊数，计算得到第一干涉相位图在像素位置(i,j)处的绝对相位；

根据第二干涉相位图在像素位置(i,j)处对应的缠绕相位和第二干涉相位图在像素位置(i,j)处的模糊数，计算得到第二干涉相位图在像素位置(i,j)处的绝对相位。

步骤7具体包括：

根据第一干涉相位图在像素位置(i,j)处对应的缠绕相位

和第一干涉相位图在像素位置(i,j)处的模糊数k₁(i,j)，计算得到第一干涉相位图在像素位置(i,j)处的绝对相位

根据第二干涉相位图在像素位置(i,j)处对应的缠绕相位

和第二干涉相位图在像素位置(i,j)处的模糊数k₂(i,j)，计算得到第二干涉相位图在像素位置(i,j)处的绝对相位

步骤8，令像素位置(i,j)遍历集合{(i,j)|1≤i≤m,1≤j≤n}中的所有元素，重复执行步骤7，从而得到第一干涉相位图的绝对相位矩阵(相位解缠绕结果)，以及第二干涉相位图的绝对相位矩阵(相位解缠绕结果)。

本发明的有效性可通过以下实验作进一步说明。

(一)InSAR仿真数据的实验

1、数据说明

仿真数据采用Colorado的Isolation Peak地区的数字高程模型DEM，分别采用基线长度112.1米和基线长度389.2米得到短基线和长基线的干涉相位图，并在干涉相位图中分别加入了相干系数为0.7和0.65的随机噪声。

2、仿真内容和结果分析

为了展示本发明方法的有效性，将本发明方法与传统TSPA方法进行了对比实验。图2(a)所示为由DEM得到的参考绝对相位。为了展示本发明方法克服噪声的能力，在进行干涉相位图仿真时加入了随机噪声，图2(b)为由图2(a)得到的干涉相位图，垂直基线(normalbaseline)长度为112.1米，加入了相干系数为0.7的随机噪声，图2(c)为由图2(a)得到的干涉相位图，垂直基线(normal baseline)长度为389.2米，加入了相干系数为0.65的随机噪声。图2(d)是采用TSPA方法求解的解缠绕结果，图2(e)是采用本发明方法求解的解缠绕结果，从这两个解缠绕相位图可以看出，TSPA方法的解缠绕结果中存在一些较为明显的色斑，而本发明方法几乎不存在类似的色斑，解缠绕结果看起来比TSPA方法相对“光滑”。图2(f)是TSPA方法的解缠绕结果与参考绝对相位的差值，差值范围为-9.3558～11.9075，图2(g)是本发明方法解缠绕结果与参考绝对相位的差值，差值范围为-6.8004～8.8865。图2(f)和图2(g)反映了求解结果与理想的参考绝对相位之间的误差大小，根据误差范围可以看出，本发明方法比传统TSPA方法的误差小。

为了能够定量的反映两种方法在精度上提高，本发明采用下式衡量解缠绕结果的精度ξ：

其中，h是由参考绝对相位得到的向量，

由解缠绕结果得到的绝对相位向量。符号||·||表示2范数。经过计算，图2(f)的ξ为0.1003，图2(g)的ξ值为0.0769，这说明了本发明方法比TSPA方法具有更高的相位解缠绕精度。

(二)实测数据实验

1、数据说明

实验数据采用我国甘肃省兰州地区的一组实测数据，基线长度分别为6.54米和40.29米。获取该干涉数据的SAR图像由TerraSAR-X系统在2010年于2月25日(轨道14975)、3月8日(轨道15142)与3月19(轨道15309)分别获得，系统轨道高度为512千米，入射角为28.7°，波长为0.03米。

2、仿真内容和结果分析

图3(a)为该地区的Google Eearth图像(该地区的纬度约为36.0405°，经度约为103.7802°)。该地区有许多陡峭的山崖与山谷，对于相位解缠绕而言，这样的复杂地形是有难度的。图3(b)是短基线干涉相位图，基线长度是6.51米，图3(c)为长基线干涉相位图，基线长度为40.29米。图3(d)为TSPA方法的解缠绕结果，图3(e)为本发明方法的解缠绕结果，为了方便对比，图3(d)与图3(e)采用了相同的颜色条(Colorbar)。通过对比可以看出，在TSPA方法的相位解缠绕图像中存在一些色斑(在图中用白色的矩形方框标识)，这说明在这些局部，解缠绕结果精度较低，而本发明方法在这些区域的图像质量却相对良好，并且通过将两种方法的解缠绕结果与图3(a)所示的谷歌地球(Google Earth)图像进行对比，本发明方法的求解结果更加接近于真实地形。为了进一步说明TSPA方法与本发明方法求解精度不同的原因，本发明根据两种方法得到的梯度估计值分别计算了残点，并在图3(f)和图3(g)中分别给出了白色方框内区域的残点分布图。对比图3(f)和图3(g)，TSPA方法估计的模糊数梯度产生了90078个残点，整个区域几乎布满了残点，本发明方法得到的残点数量为仅为18050，相对集中在干涉相位图中条纹分布变化不规律的位置。在TSPA方法和本发明方法中，在完成了模糊数梯度的估计之后，都采用了MCF模型对模糊数进行求解，然而，该优化模型在残点多为偶极子且分布较为稀疏的情况下求解的精度较高，对于残点密集的区域求解精度将会大大下降，因此，本发明方法通过众值滤波得到了更为合理的模糊数梯度估计值和相对稀疏的残点分布，因而有利于MCF模型提高求解精度。总之，通过本实验，可以说明对于实测的复杂地形，本发明方法相比于传统TSPA方法，相位解缠绕的精度得到了有效的提升。

综上，我们分别采用了仿真数据与实测数据验证了本发明的有效性。

以上所述，仅为本发明的具体实施方式，但本发明的保护范围并不局限于此，任何熟悉本技术领域的技术人员在本发明揭露的技术范围内，可轻易想到变化或替换，都应涵盖在本发明的保护范围之内。因此，本发明的保护范围应以所述权利要求的保护范围为准。

Claims (7)

1.一种基于平面近似模型的多基线相位解缠绕方法，其特征在于，所述方法包括如下步骤：

步骤1，从多基线InSAR系统中获取两个基线分别对应的干涉相位图，记为第一干涉相位图和第二干涉相位图，两幅干涉相位图的大小相同，且为m×n，其中，m、n分别为大于零的正整数；

步骤2，根据第一干涉相位图中每个像素的缠绕相位，确定第一干涉相位图中水平方向相邻像素间的水平缠绕相位梯度，记为水平缠绕第一相位梯度矩阵，其大小为m×(n-1)，并确定第一干涉相位图中垂直方向相邻像素间的垂直缠绕相位梯度，记为垂直缠绕第一相位梯度矩阵，其大小为(m-1)×n；

根据第二干涉相位图中每个像素的缠绕相位，确定第二干涉相位图中水平方向相邻像素间的水平缠绕相位梯度，记为水平缠绕第二相位梯度矩阵，其大小为m×(n-1)，并确定第二干涉相位图中垂直方向相邻像素间的垂直缠绕相位梯度，记为垂直缠绕第二相位梯度矩阵，其大小为(m-1)×n；

步骤3，设定预设窗口大小，根据预设窗口大小选定给定像素位置对应的估计窗口；

步骤4，给定像素位置(i，j)，1≤i≤m，1≤j≤n-1；

根据给定像素位置(i，j)，以及给定像素位置(i，j)对应的估计窗口，估计得到第一干涉相位图的水平方向模糊数梯度矩阵，其大小为m×(n-1)；以及估计得到第二干涉相位图的水平方向模糊数梯度矩阵，其大小为m×(n-1)；

步骤5，给定像素位置(i，j)，1≤i≤m-1，1≤j≤n；

根据给定像素位置(i，j)，以及给定像素位置(i，j)对应的估计窗口，估计得到第一干涉相位图的垂直方向模糊数梯度矩阵，其大小为(m-1)×n；以及估计得到第二干涉相位图的垂直方向模糊数梯度矩阵，其大小为(m-1)×n；

步骤6，将第一干涉相位图在水平方向的模糊数梯度矩阵中的元素以及干涉相位图1在垂直方向的模糊数梯度矩阵中的元素带入设定的优化模型，得到第一干涉相位图的模糊数矩阵1；

将第二干涉相位图在水平方向的模糊数梯度矩阵的元素以及第二干涉相位图在垂直方向的模糊数梯度矩阵的元素带入所述设定的优化模型，得到第二干涉相位图的模糊数矩阵2；

步骤7，给定像素位置(i，j)，1≤i≤m，1≤j≤n；

根据第一干涉相位图在像素位置(i，j)处对应的缠绕相位和第一干涉相位图在像素位置(i，j)处的模糊数，计算得到第一干涉相位图在像素位置(i，j)处的绝对相位；

根据第二干涉相位图在像素位置(i，j)处对应的缠绕相位和第二干涉相位图在像素位置(i，j)处的模糊数，计算得到第二干涉相位图在像素位置(i，j)处的绝对相位；

步骤8，令像素位置(i，j)遍历集合{(i，j)|1≤i≤m，1≤j≤n}中的所有元素，重复执行步骤7，从而得到第一干涉相位图的绝对相位矩阵，以及第二干涉相位图的绝对相位矩阵。

2.根据权利要求1所述的一种基于平面近似模型的多基线相位解缠绕方法，其特征在于，步骤2具体包括如下子步骤：

(2a)获取第一干涉相位图中像素位置(i，j)处的缠绕相位

获取第二干涉相位图中像素位置(i，j)处的缠绕相位

像素位置(i，j)处的像素值为像素位置(i，j)处的缠绕相位；

(2b)记第一干涉相位图的水平缠绕第一相位梯度矩阵为

中第i行第j列元素

记第二干涉相位图的水平缠绕第二相位梯度矩阵为

中第i行第j列元素

(2c)记第一干涉相位图的垂直缠绕第一相位梯度矩阵为

中第i行第j列元素

记第二干涉相位图的垂直缠绕第二相位梯度矩阵为

中第i行第j列元素

3.根据权利要求2所述的一种基于平面近似模型的多基线相位解缠绕方法，其特征在于，步骤3中，设定预设窗口大小为(2p+1)×(2q+1)；

(3a)对于水平缠绕第一相位梯度矩阵或者水平缠绕第二相位梯度矩阵中给定的像素位置(i，j)，1≤i≤m，1≤j≤n-1，其对应的估计窗口内所包含的水平缠绕相位梯度记为

或者

其中，a₀≤i＇≤a₁，b₀≤j＇≤b₁，a₀，a₁，b₀，b₁的值由下式确定：

(3b)对于垂直缠绕第一相位梯度矩阵或者垂直缠绕第二相位梯度矩阵中给定的像素位置(i，j)，1≤i≤m-1，1≤j≤n，其对应的估计窗口内所包含的垂直缠绕相位梯度记为

或者

其中，c₀≤i″≤c₁，d₀≤j″≤d₁，c₀，c₁，d₀，d₁的值由下式确定：

4.根据权利要求3所述的一种基于平面近似模型的多基线相位解缠绕方法，其特征在于，步骤4具体包括：

(4a)对于给定像素位置(i，j)，第一干涉相位图的水平缠绕相位梯度

和第二干涉相位图的水平缠绕相位梯度

以及尺寸为(2p+1)×(2q+1)的估计窗口，记估计窗口内所有水平缠绕相位梯度位置(i＇，j＇)的集合

定义R_pq误差

其中，

其中，符号[·]表示取四舍五入，

为第一干涉相位图在像素位置(i，j)的水平方向模糊数梯度，

为第二干涉相位图在像素位置(i，j)的水平方向模糊数梯度，B₁表示基线1的长度，B₂表示基线2的长度，

和

为水平梯度修正因子，用于修正窗内像素位置(i，j)处的水平绝对相位梯度与窗内其它位置(i＇，j＇)处的水平绝对相位梯度因为相位缠绕所带来的梯度不一致；

(4b)求解如下优化模型：

其中，Z表示整数集合，第一干涉相位图在像素位置(i，j)的水平方向模糊数梯度估计值

和第二干涉相位图在像素位置(i，j)的水平方向模糊数梯度估计值

是该优化模型要求解的变量。

5.根据权利要求3所述的一种基于平面近似模型的多基线相位解缠绕方法，其特征在于，步骤5具体包括：

(5a)对于给定像素位置(i，j)，第一干涉相位图的垂直缠绕相位梯度

和第二干涉相位图的垂直缠绕相位梯度

以及尺寸为(2p+1)×(2q+1)的估计窗口，记估计窗内所有垂直相位梯度位置(i″，j″)的集合

定义R_pq误差

其中，

其中，符号[·]表示取四舍五入，

为第一干涉相位图在像素位置(i，j)的垂直方向模糊数梯度，

为第二干涉相位图在像素位置(i，j)的垂直方向模糊数梯度，B₁表示基线1的长度，B₂表示基线2的长度，

和

为垂直梯度修正因子，用于修正窗内像素位置(i，j)处的垂直绝对相位梯度与窗内其它位置(i″，j″)处的垂直绝对相位梯度因为相位缠绕所带来的梯度不一致；

(5b)求解如下优化模型：

其中，Z表示整数集合，第一干涉相位图在像素位置(i，j)的垂直方向模糊数梯度估计值

和第二干涉相位图在像素位置(i，j)的垂直方向模糊数梯度估计值

是该优化模型要求解的变量。

6.根据权利要求4或5所述的一种基于平面近似模型的多基线相位解缠绕方法，其特征在于，步骤6具体包括：

(6a)将第一干涉相位图在水平方向的模糊数梯度矩阵中的元素

以及第一干涉相位图在垂直方向的模糊数梯度矩阵中的元素

带入如下设定的优化模型：

其中，Z为整数集合，

和

表示设定的权值，

和

为自由变量，k₁(i，j)、k₁(i+1，j)、k₁(i，j+1)分别为第一干涉相位图在像素位置(i，j)、(i+1，j)、(i，j+1)处的模糊数；

(6b)将第二干涉相位图在水平方向的模糊数梯度矩阵中的元素

以及第二干涉相位图在垂直方向的模糊数梯度矩阵中的元素

带入如下设定的优化模型：

其中，Z为整数集合，

和

表示设定的权值，

和

为自由变量，k₂(i，j)、k₂(i+1，j)、k₂(i，j+1)分别为第二干涉相位图在像素位置(i，j)、(i+1，j)、(i，j+1)处的模糊数。

7.根据权利要求6所述的一种基于平面近似模型的多基线相位解缠绕方法，其特征在于，步骤7具体包括：

根据第一干涉相位图在像素位置(i，j)处对应的缠绕相位

和第一干涉相位图在像素位置(i，j)处的模糊数k₁(i，j)，计算得到第一干涉相位图在像素位置(i，j)处的绝对相位

根据第二干涉相位图在像素位置(i，j)处对应的缠绕相位

和第二干涉相位图在像素位置(i，j)处的模糊数k₂(i，j)，计算得到第二干涉相位图在像素位置(i，j)处的绝对相位

Images