CN112558068A - 多基线InSAR相位估计方法及系统 - Google Patents

Abstract

本发明涉及一种多基线InSAR相位估计方法及系统。该方法包括：计算N个不同基线干涉图的概率密度；根据各概率密度，计算多基线联合概率密度的对数值；根据干涉图中的第一幅干涉图和最大缠绕数，构造流网络；计算流网络中交互边的权值；根据第一幅干涉图和权值，计算先验项的凸函数近似值；基于多基线联合概率密度的对数值和凸函数近似值，确定流网络的权值；基于权值利用最小割算法优化流网络，计算各个像素的缠绕数；更新循环变量，得到各个像素更新后的缠绕数；根据第一幅干涉图和更新后的缠绕数，计算未缠绕的真实相位，完成多基线InSAR相位估计。本发明能够解决基于变分能量项的最大后验估计方法在自然地形中性能欠缺的问题。

Description

多基线InSAR相位估计方法及系统

技术领域

本发明涉及信号处理领域，特别是涉及一种多基线InSAR相位估计方法及系统。

背景技术

多基线InSAR技术相比于单基线InSAR技术的优势在于：多基线InSAR技术可以综合利用长短基线的信息，获得更加精确的高程信息。而信息的融合即在多基线InSAR相位估计步骤中。多基线InSAR相位估计方法需要将缠绕的相位信息映射到地形高程信息。从干涉图中得到的相位是缠绕相位，其取值范围在(-π,π]之间。而我们所需的高程信息可以由真实相位与相关雷达参数计算得到。由缠绕相位求解真实相位的问题是一个病态问题，求解难度较大，往往需要额外的信息辅助。单基线InSAR技术通过相位解缠技术恢复真实相位，并在相位解缠后通过成像几何关系可以计算出逐点的高程。在多基线InSAR技术中，数据的多样性使得相位解缠的处理难度得以下降，甚至，可以绕过相位解缠步骤，直接从干涉图中估计原始相位。

多基线InSAR相位估计是多基线干涉SAR数据处理流程中的最后一步，也是至关重要的一步。相位估计的有效性和准确性对最终DEM的精度形成最为直观的影响。所以，多基线InSAR技术的研究大多集中于相位估计方法的提出与改进。其研究主要分为两个方面：对各种地形的适应能力，以及对噪声的敏感程度。相位估计方法对地形的适应能力决定了其应用场景。而相位估计方法对噪声的敏感程度决定了方法的鲁棒性。提升方法对各类地形的综合重建性能与增强方法的鲁棒性，以便更加准确地估计原始相位一直都是多基线InSAR相位估计方法所追求的目标。

传统的多基线InSAR相位估计方法可以分为两类。一类是基于非统计信号处理的方法，如中国余数定理法(Chinese Remainder Theorem method)、投影法(Projectionmethod)和线性组合方法(linear combination method)等。而另一类方法是基于统计信号处理的方法，如最大后验估计方法(Maximum a posteriori estimation,，MAP)、最大似然估计方法(Maximum likelihood estimation，MLE)等。相比于基于非统计信号处理的方法，基于统计信号处理的方法的性能较优，且效率较高。近年来，对多基线InSAR相位估计方法大多基于统计框架。其中，较为常用的便是基于变分能量项的最大后验估计方法。基于变分能量项的最大后验估计方法将变分能量项作为干涉相位的先验项，用相邻像素间的相位差绝对值来约束优化函数，以获得更加准确的解。方法在平坦区域和不连续区域均能获得较佳的性能，且处理效率较高。但是在自然地形区域(即连续但不平坦区域)，方法的估计性能明确下降，甚至发生错误。因此，对于提升最大后验估计方法在自然地形的重建性能十分必要。

发明内容

本发明的目的是提供一种多基线InSAR相位估计方法及系统，能够解决基于变分能量项的最大后验估计方法在自然地形中性能欠缺的问题。

为实现上述目的，本发明提供了如下方案：

一种多基线InSAR相位估计方法，包括：

S1：获取不同基线的缠绕相位数据，所述缠绕相位数据包括N个不同基线的干涉图、N个不同基线干涉图的相干系数、N个干涉图的基线、N个不同基线干涉图的模糊高度图、平衡常数、邻域参数、最大缠绕数以及最大迭代次数；

S2：根据所述不同基线的干涉图和所述相干系数，计算N个不同基线干涉图的概率密度；

S3：根据各所述概率密度，计算多基线联合概率密度的对数值；

S4：根据所述干涉图中的第一幅干涉图和所述最大缠绕数，构造流网络；

S5：根据所述邻域参数、所述平衡常数和所述第一幅干涉图，计算流网络中交互边的权值；

S6：根据所述第一幅干涉图和所述权值，计算先验项的凸函数近似值；

S7：基于所述多基线联合概率密度的对数值和所述凸函数近似值，确定所述流网络中数据边、约束边和交互边的权值；

S8：基于所述流网络中数据边、约束边和交互边的权值利用最小割算法优化所述流网络，计算各个像素的缠绕数；

S9：更新循环变量，重复S5-S8，直至达到最大迭代次数或所述流网络的最小割的值达到收敛，得到各个像素更新后的缠绕数；

S10：根据所述第一幅干涉图和所述更新后的缠绕数，计算未缠绕的真实相位，完成多基线InSAR相位估计。

可选地，所述根据所述干涉图中的第一幅干涉图和所述最大缠绕数，构造流网络，具体包括：

将所述第一幅干涉图数据矩阵按列转换为一个1×M的数据向量φ_vector，其中，M为干涉图φ₁的像素总个数；

将所述数据向量φ_vector扩展为一个K_max×M的矩阵φ_matrix，矩阵的每一行均为数据向量φ_vector；

将所述数据矩阵φ_matrix映射到流网络的结构中，所述数据矩阵中的数据都作为流网络中的点，完成流网络的构造。

可选地，所述根据所述第一幅干涉图和所述权值，计算先验项的凸函数近似值，具体包括：

根据所述第一幅干涉图和所述权值，计算先验项的凸函数近似值：

其中，E_approx为先验项的凸函数近似值，S为权值,β为平衡常数,w_p,q为邻域参数，P代表当前处理的像素，φ(P)代表像素P的缠绕相位值，i代表像素P的邻域像素，N_P代表像素P的邻域，M为干涉图φ₁的像素总个数。

可选地，所述基于所述多基线联合概率密度的对数值和所述凸函数近似值，确定所述流网络中数据边、约束边和交互边的权值，具体包括：

基于所述多基线联合概率密度的对数值，确定所述流网络中数据边的权值：

其中，K_max为最大缠绕数，F_log为多基线联合概率密度的对数值，s为流网络的源点，t为流网络的终点，下标(p,k)代表流网络中节点的位置，p为列数，k为行数；

确定所述流网络中约束边的权值：

其中，下标(p,i)代表流网络中节点的位置，p为列数，i为行数；

基于所述凸函数近似值，确定所述流网络中交互边的权值：

其中，g(i-j)与先验能量项相对应，g(φ(P)-φ(i))与先验项的凸函数近似值E_approx的关系为：

可选地，所述根据所述第一幅干涉图和所述更新后的缠绕数，计算未缠绕的真实相位，具体包括：

根据所述第一幅干涉图和所述更新后的缠绕数采用公式ψ＝φ₁+2Kπ，计算未缠绕的真实相位；

其中，ψ为未缠绕的真实相位，φ₁为第一幅干涉图，K为更新后的缠绕数。

一种多基线InSAR相位估计系统，包括：

缠绕相位数据获取模块，用于获取不同基线的缠绕相位数据，所述缠绕相位数据包括N个不同基线的干涉图、N个不同基线干涉图的相干系数、N个干涉图的基线、N个不同基线干涉图的模糊高度图、平衡常数、邻域参数、最大缠绕数以及最大迭代次数；

概率密度确定模块，用于根据所述不同基线的干涉图和所述相干系数，计算N个不同基线干涉图的概率密度；

对数值确定模块，用于根据各所述概率密度，计算多基线联合概率密度的对数值；

流网络构造模块，用于根据所述干涉图中的第一幅干涉图和所述最大缠绕数，构造流网络；

交互边权值确定模块，用于根据所述邻域参数、所述平衡常数和所述第一幅干涉图，计算流网络中交互边的权值；

凸函数近似值确定模块，用于根据所述第一幅干涉图和所述权值，计算先验项的凸函数近似值；

流网络有向边权值确定模块，用于基于所述多基线联合概率密度的对数值和所述凸函数近似值，确定所述流网络中数据边、约束边和交互边的权值；

缠绕数计算模块，用于基于所述流网络中数据边、约束边和交互边的权值利用最小割算法优化所述流网络，计算各个像素的缠绕数；

缠绕数更新模块，用于更新循环变量，直至达到最大迭代次数或所述流网络的最小割的值达到收敛，得到各个像素更新后的缠绕数；

未缠绕的真实相位确定模块，用于根据所述第一幅干涉图和所述更新后的缠绕数，计算未缠绕的真实相位，完成多基线InSAR相位估计。

可选地，所述流网络构造模块，具体包括：

转换单元，用于将所述第一幅干涉图数据矩阵按列转换为一个1×M的数据向量φ_vector，其中，M为干涉图φ₁的像素总个数；

扩展单元，用于将所述数据向量φ_vector扩展为一个K_max×M的矩阵φ_matrix，矩阵的每一行均为数据向量φ_vector；

映射单元，用于将所述数据矩阵φ_matrix映射到流网络的结构中，所述数据矩阵中的数据都作为流网络中的点，完成流网络的构造。

可选地，所述凸函数近似值确定模块，具体包括：

凸函数近似值确定单元，用于根据所述第一幅干涉图和所述权值，计算先验项的凸函数近似值：

其中，E_approx为先验项的凸函数近似值，S为权值,β为平衡常数,w_p,q为邻域参数，P代表当前处理的像素，φ(P)代表像素P的缠绕相位值，i代表像素P的邻域像素，N_P代表像素P的邻域，M为干涉图φ₁的像素总个数。

可选地，所述流网络有向边权值确定模块，具体包括：

流网络有向边权值确定单元，用于基于所述多基线联合概率密度的对数值，确定所述流网络中数据边的权值：

其中，K_max为最大缠绕数，F_log为多基线联合概率密度的对数值，s为流网络的源点，t为流网络的终点，下标(p,k)代表流网络中节点的位置，p为列数，k为行数；

确定所述流网络中约束边的权值：

其中，下标(p,i)代表流网络中节点的位置，p为列数，i为行数；

基于所述凸函数近似值，确定所述流网络中交互边的权值：

其中，g(i-j)与先验能量项相对应，g(φ(P)-φ(i))与先验项的凸函数近似值E_approx的关系为：

可选地，所述未缠绕的真实相位确定模块，具体包括：

真实相位计算单元，用于根据所述第一幅干涉图和所述更新后的缠绕数采用公式ψ＝φ₁+2Kπ，计算未缠绕的真实相位；

其中，ψ为未缠绕的真实相位，φ₁为第一幅干涉图，K为更新后的缠绕数。

根据本发明提供的具体实施例，本发明公开了以下技术效果：

(1)本发明提出了一种多基线InSAR相位估计方法，与现有的变分能量项的最大后验估计方法相比，有效提升了方法在各类地形的综合估计性能。

(2)本发明提出了一种多基线InSAR相位估计方法，能够适应不同地形的多基线InSAR相位估计需求。另外，本发明给出了基于非凸能量项的最大后验估计方法的处理框架，后续可以寻找选择更加贴近于干涉相位的先验模型作为能量项，进一步提升方法的综合性能。

附图说明

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

图1为本发明多基线InSAR相位估计方法流程图；

图2为流网络结构示意图；

图3为仿真使用的DEM数据的二维视图；

图4为仿真使用的DEM数据的三维视图；

图5为由原始DEM仿真得到的干涉相位图(左)和含噪声干涉图(右)，仿真所用基线为100m；

图6为由原始DEM仿真得到的干涉相位图(左)和含噪声干涉图(右)，仿真所用基线为143.4m；

图7为对不同基线含噪声干涉图像进行相位估计的结果示意图；

图8为相位估计结果与真实相位之间的残差示意图；

图9为本发明多基线InSAR相位估计系统结构图。

具体实施方式

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

本发明的目的是提供一种多基线InSAR相位估计方法及系统，能够解决基于变分能量项的最大后验估计方法在自然地形中性能欠缺的问题。

为使本发明的上述目的、特征和优点能够更加明显易懂，下面结合附图和具体实施方式对本发明作进一步详细的说明。

图1为本发明多基线InSAR相位估计方法流程图。如图1所示，一种多基线InSAR相位估计方法包括：

S1：获取不同基线的缠绕相位数据，所述缠绕相位数据包括N个不同基线的干涉图φ_i,i＝1,2,…,N，N个不同基线干涉图的相干系数图γ_i,i＝1,2,…,N，N个干涉图的基线B_i,i＝1,2,…,N,N个不同基线干涉图的模糊高度图h_ai,i＝1,2,…,N，平衡常数β，邻域参数w_p,q，最大缠绕数K_max，最大迭代次数Iter_max。

S2：根据所述不同基线的干涉图和所述相干系数，计算N个不同基线干涉图的概率密度。

结合S1获取的不同基线的干涉相位φ_i,i＝1,2,…,N，相干系数γ_i,i＝1,2,…,N，视数L，完成N个不同基线干涉相位的概率密度pdf(φ_i,h),i＝1,2,…,N的计算，具体方法为：

α＝|γ|cos(φ-φ₀) (3)

其中，Γ(·)代表gamma函数，L为多视数，₂F₁(·)为Gauss超几何分布函数，φ为未缠绕的真实相位，φ₀为缠绕相位。

S3：根据各所述概率密度，计算多基线联合概率密度的对数值。

结合S2得到的各个基线的概率密度，计算多基线联合概率密度F(φ₁,φ₂,…,φ_N)的对数值F_log，具体方法为：

S4：根据所述干涉图中的第一幅干涉图和所述最大缠绕数，构造流网络，具体包括：

将所述第一幅干涉图数据矩阵按列转换为一个1×M的数据向量φ_vector，其中，M为干涉图φ₁的像素总个数。

将所述数据向量φ_vector扩展为一个K_max×M的矩阵φ_matrix，矩阵的每一行均为数据向量φ_vector。

将所述数据矩阵φ_matrix映射到流网络的结构中，所述数据矩阵中的数据都作为流网络中的点，完成流网络的构造。

S5：根据所述邻域参数、所述平衡常数和所述第一幅干涉图，计算流网络中交互边的权值。

结合S1获取的邻域参数w_p,q、平衡常数β和第一幅干涉图φ₁，计算权重S。具体方法为：

先验能量项函数为：

结合S1获取的邻域参数w_p,q、平衡常数β和第一幅干涉图φ₁，先验能量项可以具体表示为：

其中，P为图像中的一个像素点，N_P为P点的邻域。

上述先验能量项为凹函数，对其直接优化，易陷入局部最小值。借助“超梯度”将非凸能量项转化为与之相近的凸能量项。将优化问题转化为一迭代优化问题，从而与凹能量项相近似的凸能量项。超梯度的定义为：

结合上述的先验能量项E_prior，权重可以具体表现为：

S6：根据所述第一幅干涉图和所述权值，计算先验项的凸函数近似值，具体包括：

根据所述第一幅干涉图和所述权值，计算先验项的凸函数近似值：

其中，E_approx为先验项的凸函数近似值，S为权值,β为平衡常数,w_p,q为邻域参数，P代表当前处理的像素，φ(P)代表像素P的缠绕相位值，i代表像素P的邻域像素，N_P代表像素P的邻域，M为干涉图φ₁的像素总个数。

S7：基于所述多基线联合概率密度的对数值和所述凸函数近似值，确定所述流网络中数据边、约束边和交互边的权值，具体包括：

基于所述多基线联合概率密度的对数值，确定所述流网络中数据边的权值。

确定所述流网络中约束边的权值。

基于所述凸函数近似值，确定所述流网络中交互边的权值。

图2为流网络结构，其中，s为流网络的源点，t为流网络的汇点。

(a)数据边为数据节点与源点和汇点之间的有向边，其值为：

其中，K_max为最大缠绕数，F_log为多基线联合概率密度的对数值，s为流网络的源点，t为流网络的终点，下标(p,k)代表流网络中节点的位置，p为列数，k为行数。

(b)约束边为流网络相邻行之间数据节点的有向边，其值为：

w(v_p,i+1,v_p,j)＝+∞,i＝0,1,…,K_max-1 (12)

其中，下标(p,i)代表流网络中节点的位置，p为列数，i为行数。

(c)交互边为同一行数据节点之间的有向边，其值为：

其中，g(i-j)与先验能量项相对应。g(φ(P)-φ(i))与先验项的凸函数近似值E_approx的关系为：

S8：基于所述流网络中数据边、约束边和交互边的权值利用最小割算法优化所述流网络，计算各个像素的缠绕数。

在能量函数最小化问题中，最小割算法是一中快速且有效的方法。在本文中，其求解的能量最小化问题为：

最小割算法即是在流网络中上找到具有最小权值(代价)总和的图割cut。图割cut通过切断流网络图的有向边，将图中的所有节点划分为两个不相交的子集S和T，其中源点s在子集S中，汇点t在子集T中。因此，图割cut定义为其所割断边的权值(代价)w(p,q)总和。

S9：更新循环变量，重复S5-S8，直至达到最大迭代次数或所述流网络的最小割的值达到收敛，得到各个像素更新后的缠绕数。

S10：根据所述第一幅干涉图和所述更新后的缠绕数，计算未缠绕的真实相位，完成多基线InSAR相位估计，具体包括：

根据所述第一幅干涉图和所述更新后的缠绕数，计算未缠绕的真实相位：

ψ＝φ₁+2Kπ (16)

其中，ψ为未缠绕的真实相位，φ₁为第一幅干涉图，K为更新后的缠绕数。

实施例1：

本实施例提出了一种多基线InSAR相位估计方法，其仿真计算过程所需的参数如表1所示。

表1 实施例参数

实施例1具体包括以下步骤：

步骤一：获取多个不同基线的缠绕相位数据。

获取仿真的多个不同基线的缠绕相位数据，具体包括：N个不同基线的干涉图φ_i,i＝1,2,…,N，N个不同基线干涉图的相干系数图γ_i,i＝1,2,…,N，N个干涉图的基线B_i,i＝1,2,…,N,N个不同基线干涉图的模糊高度图h_ai,i＝1,2,…,N，视数L，平衡常数β，邻域参数w_p,q，最大缠绕数K_max，最大迭代次数Iter_max。

步骤二：结合不同基线的干涉相位φ_i,i＝1,2,…,N，相干系数γ_i,i＝1,2,…,N，视数L，完成N个不同基线干涉相位的概率密度pdf(φ_i,h),i＝1,2,…,N计算N个不同基线干涉图的概率密度pdf(φ_i,h),i＝1,2,…,N,如式(1)所示。

步骤三：依据式(4)和式(5)计算多基线联合概率密度F(φ₁,φ₂,…,φ_N)的对数值F_log。

步骤四：构造流网络图G，流网络的结构如图2所示。具体流程为：

(a)将第一幅干涉图数据矩阵φ₁按列转换为一个1×M的数据向量φ_vector。其中，M为干涉图φ₁的像素总个数。

(b)将数据向量φ_vector扩展为一个K_max×M的矩阵φ_matrix，矩阵的每一行均为数据向量φ_vector。

(c)将数据矩阵φ_matrix映射到流网络的结构中，数据矩阵中的数据都作为流网络中的点，完成流网络G的构造。

步骤五：结合步骤一读入的邻域参数w_p,q、平衡常数β和第一幅干涉图φ₁，依据式(9)计算权重S。

步骤六：结合步骤一读入的第一幅干涉图φ₁和步骤五得到的权重S，依据式(10)，计算先验项的凸函数的近似值E_approx。

步骤七：确定流网络中有向边的权值，具体方法为：

图2为流网络结构，其中，s为流网络的源点，t为流网络的汇点。

(a)数据边的权值依据式(11)进行计算；

(b)约束边的权值依据式(12)进行计算；

(c)交互边的权值依据式(13)和式(14)进行计算；

步骤八：利用最小割算法优化流网络，计算各个像素的缠绕数K。

步骤九：更新循环变量i，重复步骤五至九，直至达到最大迭代次数Iter_max或网络的最小割的值达到收敛，得到各个像素的缠绕数K。

步骤十：依据式(16)计算未缠绕的真实相位ψ。

经过上述步骤的处理，得到估计的未缠绕的真实相位。图3为仿真使用的DEM数据的二维视图。图4为仿真使用的DEM数据的三维视图。原始DEM中既包含自然地形，又包含不连续区域，能够较好的验证方法对不同地形的重建性能。图5为由原始DEM仿真得到的干涉相位图(左)和含噪声干涉图(右)，仿真所用基线为100m。图6为由原始DEM仿真得到的干涉相位图(左)和含噪声干涉图(右)，仿真所用基线为143.4m。含噪声干涉图的相干系数为0.5，满足大部分星载InSAR图像的要求。图7为对不同基线含噪声干涉图像进行相位估计的结果示意图，图8为相位估计结果与真实相位之间的残差示意图。由图8可以明显看出，本发明所提出方法能够从含噪声地干涉相位图像中正确地恢复未缠绕地真实相位，方法的准确度得到保证。

图9为本发明多基线InSAR相位估计系统结构图。如图9所示，一种多基线InSAR相位估计系统包括：

缠绕相位数据获取模块101，用于获取不同基线的缠绕相位数据，所述缠绕相位数据包括N个不同基线的干涉图、N个不同基线干涉图的相干系数、N个干涉图的基线、N个不同基线干涉图的模糊高度图、平衡常数、邻域参数、最大缠绕数以及最大迭代次数。

概率密度确定模块102，用于根据所述不同基线的干涉图和所述相干系数，计算N个不同基线干涉图的概率密度。

对数值确定模块103，用于根据各所述概率密度，计算多基线联合概率密度的对数值。

流网络构造模块104，用于根据所述干涉图中的第一幅干涉图和所述最大缠绕数，构造流网络。

交互边权值确定模块105，用于根据所述邻域参数、所述平衡常数和所述第一幅干涉图，计算流网络中交互边的权值。

凸函数近似值确定模块106，用于根据所述第一幅干涉图和所述权值，计算先验项的凸函数近似值。

流网络有向边权值确定模块107，用于基于所述多基线联合概率密度的对数值和所述凸函数近似值，确定所述流网络中数据边、约束边和交互边的权值。

缠绕数计算模块108，用于基于所述流网络中数据边、约束边和交互边的权值利用最小割算法优化所述流网络，计算各个像素的缠绕数。

缠绕数更新模块109，用于更新循环变量，直至达到最大迭代次数或所述流网络的最小割的值达到收敛，得到各个像素更新后的缠绕数。

未缠绕的真实相位确定模块110，用于根据所述第一幅干涉图和所述更新后的缠绕数，计算未缠绕的真实相位，完成多基线InSAR相位估计。

所述流网络构造模块104，具体包括：

转换单元，用于将所述第一幅干涉图数据矩阵按列转换为一个1×M的数据向量φ_vector，其中，M为干涉图φ₁的像素总个数。

扩展单元，用于将所述数据向量φ_vector扩展为一个K_max×M的矩阵φ_matrix，矩阵的每一行均为数据向量φ_vector。

映射单元，用于将所述数据矩阵φ_matrix映射到流网络的结构中，所述数据矩阵中的数据都作为流网络中的点，完成流网络的构造。

所述凸函数近似值确定模块106，具体包括：

凸函数近似值确定单元，用于根据所述第一幅干涉图和所述权值，计算先验项的凸函数近似值：

其中，E_approx为先验项的凸函数近似值，S为权值,β为平衡常数,w_p,q为邻域参数，P代表当前处理的像素，φ(P)代表像素P的缠绕相位值，i代表像素P的邻域像素，N_P代表像素P的邻域，M为干涉图φ₁的像素总个数。

所述流网络有向边权值确定模块107，具体包括：

流网络有向边权值确定单元，用于基于所述多基线联合概率密度的对数值，确定所述流网络中数据边的权值：

其中，K_max为最大缠绕数，F_log为多基线联合概率密度的对数值，s为流网络的源点，t为流网络的终点，下标(p,k)代表流网络中节点的位置，p为列数，k为行数。

确定所述流网络中约束边的权值：

其中，下标(p,i)代表流网络中节点的位置，p为列数，i为行数。

基于所述凸函数近似值，确定所述流网络中交互边的权值：

其中，g(i-j)与先验能量项对应。g(φ(P)-φ(i))与先验项的凸函数近似值E_approx的关系为：

所述未缠绕的真实相位确定模块110，具体包括：

真实相位计算单元，用于根据所述第一幅干涉图和所述更新后的缠绕数采用公式ψ＝φ₁+2Kπ，计算未缠绕的真实相位。

其中，ψ为未缠绕的真实相位，φ₁为第一幅干涉图，K为更新后的缠绕数。

本说明书中各个实施例采用递进的方式描述，每个实施例重点说明的都是与其他实施例的不同之处，各个实施例之间相同相似部分互相参见即可。对于实施例公开的系统而言，由于其与实施例公开的方法相对应，所以描述的比较简单，相关之处参见方法部分说明即可。

本文中应用了具体个例对本发明的原理及实施方式进行了阐述，以上实施例的说明只是用于帮助理解本发明的方法及其核心思想；同时，对于本领域的一般技术人员，依据本发明的思想，在具体实施方式及应用范围上均会有改变之处。综上所述，本说明书内容不应理解为对本发明的限制。

Claims (10)

1.一种多基线InSAR相位估计方法，其特征在于，包括：

S1：获取不同基线的缠绕相位数据，所述缠绕相位数据包括N个不同基线的干涉图、N个不同基线干涉图的相干系数、N个干涉图的基线、N个不同基线干涉图的模糊高度图、平衡常数、邻域参数、最大缠绕数以及最大迭代次数；

S2：根据所述不同基线的干涉图和所述相干系数，计算N个不同基线干涉图的概率密度；

S3：根据各所述概率密度，计算多基线联合概率密度的对数值；

S4：根据所述干涉图中的第一幅干涉图和所述最大缠绕数，构造流网络；

S5：根据所述邻域参数、所述平衡常数和所述第一幅干涉图，计算流网络中交互边的权值；

S6：根据所述第一幅干涉图和所述权值，计算先验项的凸函数近似值；

S7：基于所述多基线联合概率密度的对数值和所述凸函数近似值，确定所述流网络中数据边、约束边和交互边的权值；

S8：基于所述流网络中数据边、约束边和交互边的权值利用最小割算法优化所述流网络，计算各个像素的缠绕数；

S9：更新循环变量，重复S5-S8，直至达到最大迭代次数或所述流网络的最小割的值达到收敛，得到各个像素更新后的缠绕数；

S10：根据所述第一幅干涉图和所述更新后的缠绕数，计算未缠绕的真实相位，完成多基线InSAR相位估计。

2.根据权利要求1所述的多基线InSAR相位估计方法，其特征在于，所述根据所述干涉图中的第一幅干涉图和所述最大缠绕数，构造流网络，具体包括：

将所述第一幅干涉图数据矩阵按列转换为一个1×M的数据向量φ_vector，其中，M为干涉图φ₁的像素总个数；

将所述数据向量φ_vector扩展为一个K_max×M的矩阵φ_matrix，矩阵的每一行均为数据向量φ_vector；

将所述数据矩阵φ_matrix映射到流网络的结构中，所述数据矩阵中的数据都作为流网络中的点，完成流网络的构造。

3.根据权利要求1所述的多基线InSAR相位估计方法，其特征在于，所述根据所述第一幅干涉图和所述权值，计算先验项的凸函数近似值，具体包括：

根据所述第一幅干涉图和所述权值，计算先验项的凸函数近似值：

其中，E_approx为先验项的凸函数近似值，S为权值,β为平衡常数,w_p,q为邻域参数，P代表当前处理的像素，φ(P)代表像素P的缠绕相位值，i代表像素P的邻域像素，N_P代表像素P的邻域，M为干涉图φ₁的像素总个数。

4.根据权利要求1所述的多基线InSAR相位估计方法，其特征在于，所述基于所述多基线联合概率密度的对数值和所述凸函数近似值，确定所述流网络中数据边、约束边和交互边的权值，具体包括：

基于所述多基线联合概率密度的对数值，确定所述流网络中数据边的权值：

其中，K_max为最大缠绕数，F_log为多基线联合概率密度的对数值，s为流网络的源点，t为流网络的终点，下标(p,k)代表流网络中节点的位置，p为列数，k为行数；

确定所述流网络中约束边的权值：

w(v_p,i+1,v_p,i)＝+∞,i＝0,1,…,K_max-1；

其中，下标(p,i)代表流网络中节点的位置，p为列数，i为行数；

基于所述凸函数近似值，确定所述流网络中交互边的权值：

其中，g(i-j)与先验能量项对应，g(φ(P)-φ(i))与先验项的凸函数近似值E_approx的关系为：

5.根据权利要求1所述的多基线InSAR相位估计方法，其特征在于，所述根据所述第一幅干涉图和所述更新后的缠绕数，计算未缠绕的真实相位，具体包括：

根据所述第一幅干涉图和所述更新后的缠绕数采用公式ψ＝φ₁+2Kπ，计算未缠绕的真实相位，其中，ψ为未缠绕的真实相位，φ₁为第一幅干涉图，K为更新后的缠绕数。

6.一种多基线InSAR相位估计系统，其特征在于，包括：

缠绕相位数据获取模块，用于获取不同基线的缠绕相位数据，所述缠绕相位数据包括N个不同基线的干涉图、N个不同基线干涉图的相干系数、N个干涉图的基线、N个不同基线干涉图的模糊高度图、平衡常数、邻域参数、最大缠绕数以及最大迭代次数；

概率密度确定模块，用于根据所述不同基线的干涉图和所述相干系数，计算N个不同基线干涉图的概率密度；

对数值确定模块，用于根据各所述概率密度，计算多基线联合概率密度的对数值；

流网络构造模块，用于根据所述干涉图中的第一幅干涉图和所述最大缠绕数，构造流网络；

交互边权值确定模块，用于根据所述邻域参数、所述平衡常数和所述第一幅干涉图，计算流网络中交互边的权值；

凸函数近似值确定模块，用于根据所述第一幅干涉图和所述权值，计算先验项的凸函数近似值；

流网络有向边权值确定模块，用于基于所述多基线联合概率密度的对数值和所述凸函数近似值，确定所述流网络中数据边、约束边和交互边的权值；

缠绕数计算模块，用于基于所述流网络中数据边、约束边和交互边的权值利用最小割算法优化所述流网络，计算各个像素的缠绕数；

缠绕数更新模块，用于更新循环变量，直至达到最大迭代次数或所述流网络的最小割的值达到收敛，得到各个像素更新后的缠绕数；

未缠绕的真实相位确定模块，用于根据所述第一幅干涉图和所述更新后的缠绕数，计算未缠绕的真实相位，完成多基线InSAR相位估计。

7.根据权利要求6所述的多基线InSAR相位估计系统，其特征在于，所述流网络构造模块，具体包括：

转换单元，用于将所述第一幅干涉图数据矩阵按列转换为一个1×M的数据向量φ_vector，其中，M为干涉图φ₁的像素总个数；

扩展单元，用于将所述数据向量φ_vector扩展为一个K_max×M的矩阵φ_matrix，矩阵的每一行均为数据向量φ_vector；

映射单元，用于将所述数据矩阵φ_matrix映射到流网络的结构中，所述数据矩阵中的数据都作为流网络中的点，完成流网络的构造。

8.根据权利要求6所述的多基线InSAR相位估计系统，其特征在于，所述凸函数近似值确定模块，具体包括：

凸函数近似值确定单元，用于根据所述第一幅干涉图和所述权值，计算先验项的凸函数近似值：

其中，E_approx为先验项的凸函数近似值，S为权值,β为平衡常数,w_p,q为邻域参数，P代表当前处理的像素，φ(P)代表像素P的缠绕相位值，i代表像素P的邻域像素，N_P代表像素P的邻域，M为干涉图φ₁的像素总个数。

9.根据权利要求6所述的多基线InSAR相位估计系统，其特征在于，所述流网络有向边权值确定模块，具体包括：

流网络有向边权值确定单元，用于基于所述多基线联合概率密度的对数值，确定所述流网络中数据边的权值：

其中，K_max为最大缠绕数，F_log为多基线联合概率密度的对数值，s为流网络的源点，t为流网络的终点，下标(p,k)代表流网络中节点的位置，p为列数，k为行数；

确定所述流网络中约束边的权值：

w(v_p,i+1,v_p,i)＝+∞,i＝0,1,…,K_max-1；

其中，下标(p,i)代表流网络中节点的位置，p为列数，i为行数；

基于所述凸函数近似值，确定所述流网络中交互边的权值：

其中，g(i-j)与先验能量项对应。g(φ(P)-φ(i))与先验项的凸函数近似值E_approx的关系为：

10.根据权利要求6所述的多基线InSAR相位估计系统，其特征在于，所述未缠绕的真实相位确定模块，具体包括：

真实相位计算单元，用于根据所述第一幅干涉图和所述更新后的缠绕数采用公式ψ＝φ₁+2Kπ，计算未缠绕的真实相位，其中，ψ为未缠绕的真实相位，φ₁为第一幅干涉图，K为更新后的缠绕数。

Images