CN109541593B - 一种改进的最小费用流InSAR相位解缠方法 - Google Patents

Abstract

本发明涉及一种改进的最小费用流InSAR相位解缠方法，最小费用流相位解缠模型是全局优化，在计算到相位梯度之后，进行路径积分时容易造成误差叠加。当合成孔径雷达获取的目标地形过于陡峭，相邻像素间目标的高程差超出了模糊范围时，会造成真实相位值存在跳变，导致最小费用流算法存在解缠错误传递。地面控制点能够辅助InSAR相位解缠，提高解缠成功率。但是现有的相位解缠算法均是在获得解缠结果后利用控制点进行校正，而没有真正将控制点信息融入到最小费用流解缠算法中。本发明通过构造合理的优化目标函数，减小了最小费用流解缠模型因为路径积分引起的噪声叠加，同时在解缠算法的模型中引入控制点信息，使得解缠的结果更精确。

Description

一种改进的最小费用流InSAR相位解缠方法

技术领域

本发明涉及一种改进的最小费用流InSAR相位解缠方法，属于信号处理领域。

背景技术

在InSAR数据处理过程中，相位解缠是合成孔径雷达干涉测量的关键流程，它的准确性直接影响到InSAR生成数字高程模型的精确性。现有的相位解缠方法均假设各相邻像素的干涉相位差的绝对值小于π。但是，阴影、去相关等因素引起的噪声和突变地形往往造成相位数据不连续，给相位解缠带来了极大的困难。目前大部分算法都无法圆满地解决这些问题，解缠的结果常常会存在较大的误差，由此得到的数字高程模型就会与实际情况存在较大的差别。如何能够从质量较差的数据当中提取有用的信息，而抑制噪声对解缠过程的影响，成为一个亟待解决的问题。

传统的最小费用流解缠算法最早见于Costantini等提出的基于网络流的相位解缠方法，这种方法是将相位解缠问题转化为最小化问题，通过在全局范围内搜索路径和最短枝切来求得最小化问题的最优解。该方法可以应用于规则网络(网格)，也可以用于不规则网络(三角网)。

传统最小费用流解缠算法存在的问题：第一，作为全局优化的模型，当某一个点解缠错误时，路径积分会引发通过该点路径上的其他点相位解缠错误；第二，在路径积分的过程中，由于每个像素点的噪声叠加，从而扩大噪声的影响，可能会引起相位的跳变，造成解缠的错误。第三，图像相邻像素对应的地形高程出现了大于模糊高程的突变时，即相邻像素相位真实值本身存在相位跳变时，最小费用流相位解缠往往出现区域性的解缠错误。

本发明一方面构建噪声最小的优化模型，降低了噪声叠加对解缠结果的影响，另一方面在模型中引入控制点矫正一项，能够有效减少传统最小费用流因为复杂地形而引起的区域性解缠错误，提高解缠成功率。

发明内容

本发明技术解决问题：为了解决传统最小费用流解缠的问题，提供一种改进的最小费用流InSAR相位解缠方法，有效减少传统最小费用流相位解缠因为复杂地形引起的区域性解缠错误，提高解缠成功率。

本发明技术解决方案：一种改进的最小费用流InSAR相位解缠方法，实现为：

第一步：利用最小费用流算法分别获得图像相位矩阵距离向和方位向的离散导数残差矩阵；

第二步：基于第一步所述的离散导数残差矩阵，建立噪声最小的优化模型，降低噪声引起的误差，在噪声最小的优化模型中引入多个控制点相位来进一步校正解缠结果，提高解缠的准确率。

第三步：对包含控制点校正项的最优化函数求导，并利用奇异值分解方法求解导数为零的值，也即原函数的最优解，得出的解为解缠相位。

所述第一步具体实现如下：

(11)利用距离向和方位向偏导数表示离散导数残差值；

设大小为M行N列的缠绕干涉相位矩阵为：

F_i,j，(i,j)∈S,S＝{(i,j)|0≤i≤M-1,0≤j≤N-1} (1.1.1)

最小费用流利用缠绕函数来构建解缠函数，首先利用缠绕相位的离散偏导数作为真实相位离散偏导数的初始估计值，缠绕函数的相位梯度是一个有旋场，解缠函数的相位梯度为无旋场，相位解缠的问题转化为寻找离散偏导数的残差：

S₁＝{(i,j)|0≤i≤M-1,0≤j≤N-2}

S₂＝{(i,j)|0≤i≤M-2,0≤j≤N-1}

k₁(i,j)、k₂(i,j)表示图像(i,j)点距离向和方位向离散倒数残差值，Δ₁φ(i,j)、Δ₂φ(i,j)表示图像(i,j)点距离向和方位向解缠函数相位梯度，ψ₁(i,j)、ψ₂(i,j)表示图像(i,j)点距离向和方位向缠绕函数的相位梯度；

(12)利用最小费用流优化模型求解离散导数残差值并得到离散导数残差矩阵；

求解真实相位的梯度差变成求解k₁(i,j)和k₂(i,j)，即最小化所有的k₁(i,j)和k₂(i,j)之和，

其中c₁(i,j)、c₂(i,j)为加权系数矩阵，同时目标函数满足下列约束条件：

k₁(i,j+1)-k₁(i,j)-k₂(i+1,j)+k₂(i,j)＝[ψ₁(i,j+1)-ψ₁(i,j)-ψ₂(i+1,j)+ψ₂(i,j)]/2π (1.1.4)

k₁,k₂都是整数，通过传统最小费用流算法获得k₁,k₂形成矩阵K1,K2。

所述第二步具体实现如下：

(21)构建噪声最小的二范数最优化函数，

考虑干涉SAR获得的缠绕相位中含有噪声项，有下式：

其中nos₁(i,j)、nos₂(i,j)分别表示图像(i,j)点距离向和方位向的缠绕相位梯度中包含的噪声相位，构建最优化函数使得噪声项达到最小，并将离散倒数残差矩阵代入构建的最优化函数中；

其中B₁、B₂对应为距离向和方位向梯度算子系数矩阵，K₁,K₂为将离散导数残差矩阵K1,K2列向量化得到的列向量，ψ为数据缠绕相位列向量化得到的列向量，φ为与ψ相同维度的待求解的解缠相位列向量，通过求解最优化目标函数得到φ向量的数值后，将列向量φ矩阵化恢复出真实解缠相位矩阵；

(22)在最优化函数中加入控制点校正项，所述控制点校正项为二范数最优化模型，并设置权重λ，

A为控制点提取系数矩阵，

为控制点真实相位矩阵。

所述第三步具体实现如下：

(31)对包含控制点校正项的优化目标函数进行二范数计算化简，求解最优化目标函数也即求解极值点，需要求解函数导数为零时的解，对函数求导，取函数的导数为零建立方程；

(32)利用奇异值分解法(SVD)求解步骤(31)的方程，得到的结果为原优化函数的极值点即最优解。

本发明基于最小费用流与多控制点信息的相位解缠方法的优点在于：

(1)实用性。在原有的最小费用流解缠方法的基础上进行改进，降低了难度。

(2)有效性。能够有效校正最小费用流解缠方法在复杂地形出现的区域性相位解缠错误。现有的许多解缠算法，都是依靠路径积分或者以其中一个高质量点作为起始点，向外部扩散的解缠方式，这些解缠方式，当在路径中遇到一个点解缠错误，后面所有经过该点的相位全部都会出现解缠错误。而且，虽然相邻像素点之间的噪声很小，但是随着路径积分噪声会不断的叠加，有可能会突破π的界限从而发生相位跳变，这种情况同样会引起相位解缠错误。本发明首先构建噪声最小的优化模型降低了噪声在路径积分的叠加对解缠成功率的影响。同时引入控制点，当出现解缠错误之后，通过错误点的路径积分靠近控制点时，控制点的存在会校正错误，不会出现区域性相位解缠错误。本发明以200×200大小的复杂地形SAR干涉相位图进行仿真验证，传统最小费用流解缠成功率为72.16％，本发明解缠成功率为99.05％。

附图说明

图1改进的最小费用流InSAR相位解缠方法步骤；

图2基于实际高程数据利用干涉相位模型仿真出的真实相位；

图3传统最小费用流算法解缠结果；

图4改进的最小费用流InSAR相位解缠方法获得的结果；

图5真实相位和最小费用流算法解缠相位做差的结果；

图6真实相位和改进的最小费用流InSAR相位解缠方法获得相位做差的结果。

具体实施方式

下面将结合附图和实施例对本发明作进一步的详细说明。

如图1所示，本发明一种改进的最小费用流InSAR相位解缠方法，具体实现如下：

步骤一：设大小为M行N列的缠绕干涉相位矩阵为：

F_i,j，(i,j)∈S,S＝{(i,j)|0≤i≤M-1,0≤j≤N-1}

Δ₁F_i,j为距离向离散偏导数，Δ₂F_i,j为方位向离散偏导数，则有：

其中：

S₁＝{(i,j)|0≤i≤M-1,0≤j≤N-2}

S₂＝{(i,j)|0≤i≤M-2,0≤j≤N-1}

如果有：

F₁(i,j+1)-F₁(i,j)＝F₂(i+1,j)-F₂(i,j),0≤i≤M-2,0≤j≤N-2 (1.2.3)则该函数的梯度场为无旋场，也就是说，沿着任意路径积分结果都一致。

对于利用缠绕函数来构建解缠函数，首先利用缠绕相位的离散偏导数作为真实相位离散偏导数的初始估计值，由于相位缠绕，两者相差2π的整数倍。

其中整数n₁(i,j)、n₂(i,j)分别表示图像(i,j)点距离向和方位向的缠绕相位离散偏导数与真实相位离散偏导数相差2π的倍数。由于ψ(i,j)为缠绕函数的相位梯度，Δφ_i,j为解缠函数的相位梯度，前者是一个有旋场，后者为无旋场。因此相位解缠的问题就转化为如何寻找离散偏导数的残差：

k₁(i,j)＝[Δ₁φ(i,j)-ψ₁(i,j)]/2π,(i,j)∈S₁

k₂(i,j)＝[Δ₂φ(i,j)-ψ₂(i,j)]/2π,(i,j)∈S₂

求解真实相位的梯度就变成了求解k₁(i,j)和k₂(i,j)，即最小化所有的k₁(i,j)和k₂(i,j)之和，

同时满足下列约束条件：

k₁(i,j+1)-k₁(i,j)-k₂(i+1,j)+k₂(i,j)＝[ψ₁(i,j+1)-ψ₁(i,j)-ψ₂(i+1,j)+ψ₂(i,j)]/2π

k₁,k₂都是整数。通过传统最小费用流算法获得k₁,k₂形成矩阵K1,K2。

步骤二：利用二范数最优化，求取噪声最小时的最优解。降低噪声引起的误差，并且引入多个控制点相位来矫正原结果，提高解缠的准确率。考虑式(1.1.2)中的噪声项，即：

nos₁(i,j)＝2πk₁(i,j)-[Δφ₁(i,j)-Δψ₁(i,j)]

nos₂(i,j)＝2πk₂(i,j)-[Δφ₂(i,j)-Δψ₂(i,j)]

构建优化函数使得噪声项达到最小。

原始图像的解缠相位数据φ为M×N的矩阵，将其化为N₁×1的列向量，其中N₁＝M×N,K₁、K₂是最小费用流算法解出的距离向和方位向离散倒数残差矩阵K1,K2列向量化得到的列向量，列向量K₁的维度为N₂1,其中N₂为数值，其大小为M(N-1)、列向量K₂的维度为N₃×1，其中N₃为数值，其大小为(M-1)×N，设B₁,B₂分别是是距离向和方位向的梯度算子系数矩阵。

B₁为距离向梯度算子系数矩阵，维度为N₂×N₁。定义如下：

其中，ε为正整数，ε∈1,2...N2，则有：

B₂是方位向梯度算子系数矩阵，维度为N₃×N₁，定义如下：

其中ε∈1,2...N₃

步骤三：设SAR图像中共有H个控制点，则有控制点提取系数矩阵A，维度为H×N₁,

其中i∈1,2...H，j的取值分别对应真实相位数据中H个控制点在被列向量化的φ列向量中的H个索引值。例如原始图像为200×200，坐标为(20,70)的控制点在φ中的索引值为69×200+20。

为控制点真实相位向量，其为维度N₁×1的稀疏列向量，与列向量φ维度相同，

向量中，在H个控制点索引值处的数值为利用对应控制点实际地理位置与卫星位置以及空间几何关系得到的相位非缠绕值，即控制点真实相位，其余索引值上的数值为0。

λ为控制点校正项的权重(取值3到10，本发明实施实例取值为5)。ψ为由干涉SAR获取到的原始图像缠绕相位矩阵列向量化的列向量，即缠绕相位列向量。

最终构造的二范数优化目标函数为下：

其中前两项为噪声最小二范数优化函数，第三项为控制点矫正项经过化简求导及令导函数为零，求得最优化目标函数的解方程为：

利用奇异值分解可以对上述方程进行求解得到结果。

为说明本发明的有效性，进行如下验证。

第一步：利用实际的高程数据，用以下的干涉SAR相位与高程的函数关系来仿真相位数据。

模糊高程：

缠绕相位

为仿真的非缠绕相位，即真实相位，λ为波长，Δb_⊥为干涉SAR垂直基线长度，r为卫星与目标距离，θ为雷达下视角，h为目标高程，

为噪声相位，ψ_matrix为将

进行缠绕的缠绕相位矩阵，即数据相位矩阵。

表1真实相位的仿真参数(基于TerraSAR系统)

第二步：对仿真的真实相位数据进行缠绕。

第三步：分别利用传统的最小费用流算法和本发明提出的改进最小费用流相位解缠方法对缠绕相位进行解缠。

第四步：将真实相位与两种算法的解缠结果做差进行比较。

分析解缠得到的结果，图2为仿真出的真实相位，图3为传统最小费用流解缠结果，而图4为本发明提出的解缠方法的结果，根据这三幅图的对比，可以看出，真实相位下部分发生跳变，但是下部分相位数值比中间区域要大，传统最小费用流解缠的结果下部分同样发生相位跳变，但是相位却比中间区域相位小，而本发明提出的解缠方法在相位跳变边界附近其变化趋势与真实相位基本一致。为了更直观比较两种方法的优劣，图5和图6分别为传统最小费用流解缠结果与真实相位的差值，以及本发明的模型解缠结果与真实相位的差值，当真实相位中出现了跳变时，利用最小费用流解缠的时候，整个下部分全部出现了解缠错误，相位误差在6到7之间，而利用改进算法后可以看出，虽然在跳变的像素附近出现了误差，但是对于全体相位而言，大大减少了解缠错误点的数目，有效提高了解缠的准确率。并且通过计算两者与真实相位差值的均值可得最小费用流为1.7256而最优化模型为0.0966，同时计算两者与真实相位值之差是否在在π之内来确定解缠成功率，传统最小费用流成功率为72.16％，本发明提出的模型成功率为99.05％，由此证实相比于传统最小费用流本发明提出的解缠方法大大降低了解缠的误差，提高了解缠成功率。本发明提出的算法通过增加控制点可以降低解缠的错误率，提高精确度，但是同时会增加对控制点信息的需求，以及计算的复杂度。

提供以上实施例仅仅是为了描述本发明的目的，而并非要限制本发明的范围。本发明的范围由所附权利要求限定。不脱离本发明的精神和原理而做出的各种等同替换和修改，均应涵盖在本发明的范围之内。

Claims (2)

1.一种改进的最小费用流InSAR相位解缠方法，其特征在于，包括以下步骤：

第一步：利用最小费用流算法分别获得图像相位矩阵距离向和方位向的离散导数残差矩阵；

第二步：基于第一步所述的离散导数残差矩阵，建立噪声最小的最优化目标函数，降低噪声引起的误差，在噪声最小的最优化目标函数中引入多个控制点相位来进一步校正解缠结果，提高解缠的准确率；

第三步：对包含控制点校正项的最优化函数求导，并利用奇异值分解方法求解导数为零的值，也即原函数的最优解，得出的解为解缠相位；

所述第二步具体实现如下：

(21)构建噪声最小的二范数最优化函数，

考虑干涉SAR获得的缠绕相位中含有噪声项，有下式：

其中nos₁(i,j)、nos₂(i,j)分别表示图像(i,j)点距离向和方位向的缠绕相位梯度中包含的噪声相位，构建最优化函数使得噪声项达到最小，并将离散导数残差矩阵代入构建的最优化函数中；

其中B₁为距离向梯度算子系数矩阵，B₂为方位向梯度算子系数矩阵，K₁,K₂为将离散导数残差矩阵K1,K2列向量化得到的列向量，ψ为缠绕相位矩阵(干涉SAR获取到的相位数据)列向量化得到的列向量，φ为与ψ相同维度的待求解的解缠相位列向量，B₁矩阵分别与列向量φ和ψ相乘可以对应得到距离向解缠相位梯度向量和距离向缠绕相位梯度向量，同理B₂矩阵分别与列向量φ和ψ相乘可以对应得到方位向解缠相位梯度向量和方位向缠绕相位梯度向量；

(22)在最优化函数中加入控制点校正项，所述控制点校正项为二范数最优化模型，并设置权重λ，

A为控制点提取系数矩阵，

为控制点真实相位矩阵；

所述第三步具体实现如下：

(31)对包含控制点校正项的优化目标函数进行二范数计算化简，求解最优化目标函数也即求解极值点，需要求解函数导数为零时的解，对函数求导，取函数的导数为零建立方程；

(32)利用奇异值分解法(SVD)求解步骤(31)的方程，得到φ向量的数值后，将列向量φ矩阵化恢复出真实解缠相位矩阵。

2.根据权利要求1所述的一种改进的最小费用流InSAR相位解缠方法，其特征在于：所述第一步具体实现如下：

(11)利用距离向和方位向偏导数表示离散导数残差值；

设大小为M行N列的缠绕干涉相位矩阵为：

F_i,j，(i,j)∈S,S＝{(i,j)|0≤i≤M-1,0≤j≤N-1} (1.1.1)

其中i表示行坐标即方位向坐标，j表示列坐标即距离向坐标；

最小费用流利用缠绕函数来构建解缠函数，首先利用缠绕相位的离散偏导数作为真实相位离散偏导数的初始估计值，缠绕函数的相位梯度是一个有旋场，解缠函数的相位梯度为无旋场，相位解缠的问题转化为寻找离散偏导数的残差：

S₁＝{(i,j)|0≤i≤M-1,0≤j≤N-2}

S₂＝{(i,j)|0≤i≤M-2,0≤j≤N-1}

k₁(i,j)、k₂(i,j)表示图像(i,j)点距离向和方位向离散导数残差值，Δ₁φ(i,j)、Δ₂φ(i,j)表示图像(i,j)点距离向和方位向解缠函数相位梯度，ψ₁(i,j)、ψ₂(i,j)表示图像(i,j)点距离向和方位向缠绕函数的相位梯度；

(12)利用最小费用流优化模型求解离散导数残差值并得到离散导数残差矩阵；

求解真实相位的梯度差变成求解k₁(i,j)和k₂(i,j)，即最小化所有的k₁(i,j)和k₂(i,j)之和，

其中c₁(i,j)、c₂(i,j)为加权系数矩阵，同时目标函数满足下列约束条件：

k₁(i,j+1)-k₁(i,j)-k₂(i+1,j)+k₂(i,j)＝[ψ₁(i,j+1)-ψ₁(i,j)-ψ₂(i+1,j)+ψ₂(i,j)]/2π (1.1.4)

k₁,k₂都是整数，通过传统最小费用流的算法获得k₁,k₂形成矩阵K1,K2。

Images