CN105005046A - 基于无网格与频率估计的干涉合成孔径雷达相位解缠方法 - Google Patents

Abstract

本发明提供一种基于无网格与频率估计的干涉合成孔径雷达相位解缠方法，能够提高干涉合成孔径雷达干涉相位图的像素点之间解缠精度。该相位解缠方法包括如下步骤：获取干涉合成孔径雷达的干涉相位图；将所述干涉相位图划分为多个子干涉相位图，并计算每个子干涉相位图中相邻像素点之间的缠绕相位梯度值；确定所述每个子干涉相位图的形函数；根据所述每个子干涉相位图中相邻像素点之间的缠绕相位梯度值和所述形函数，得到所述每个子干涉相位图的缠绕相位梯度函数；根据所述每个子干涉相位图的缠绕相位梯度函数，求解所述每个子干涉相位图的解缠相位值；根据所述每个子干涉相位图的解缠相位值，对所述干涉合成孔径雷达的干涉相位图进行解缠。

Description

基于无网格与频率估计的干涉合成孔径雷达相位解缠方法

技术领域

本发明属于雷达信号处理技术领域，更进一步是一种基于无网格与频率估计的干涉合成孔径雷达相位解缠方法。通过对干涉合成孔径雷达干涉相位图的像素点之间受到缠绕相位应用基于无网格法和局部频率的InSAR相位解缠法，得到目标真实相位，提高解缠精度，从而使InSAR数据处理精度得到显著改善。

背景技术

上世纪中叶开始，雷达技术获得迅速发展，许多新颖的雷达技术陆续被提出并使用。干涉合成孔径雷达(Interferometric Synthetic Aperture Radar(以下简称InSAR)技术便是代表之一，它是射电天文学干涉测量技术与合成孔径雷达技术两者相结合而产生的全新技术。在继承了合成孔径雷达获得目标二位图像的特点上，又具备了获取目标区域高程信息的能力。

InSAR是由Carl Wiley在SAR的基础上首次提出的。它的工作原理是通过对同一目标区域的具有相干性的两幅SAR图像进行干涉处理，得到一幅该区域的干涉相位图像，分析处理干涉相位图像从而获取平面和高程方向的信息，完成三维测量。

相位信息蕴含着距离差的信息，在干涉处理环节中，获取相位信息是整个技术的关键所在。经过一系列的预处理得到的干涉相位图中的相位信息是受到干扰的，相位信息中的相位值被限制在了[-π,π]之间，也就是说干涉相位图像中的相位都是存在误差的缠绕相位。对存在误差的缠绕相位进行解缠并获得真实的相位信息，是解决相位缠绕问题的必要手段，根据已知的缠绕相位，利用邻近像素点间的缠绕相位梯度获得与实际距离差成比例的相位的过程就是相位解缠的思想。

1951年InSAR技术被发明至今，针对相位解缠的方法不断涌现，目前已经有了几十种相位解缠算法。按照解缠方式的不同，大体可分为三类，基于路径跟踪的相位解缠方法、基于最小范数的相位解缠方法和基于网络规划的相位解缠方法。最小二乘法作为基于最小范数类算法的代表，其思想是以缠绕相位梯度与解缠相位梯度的均方误差最小作为约束。最小二乘法具备了良好的稳健性，解缠过程中不依赖积分路径，但是缺点是不对残差点进行处理，最终导致残差点引起的误差扩散到整个图像。

发明内容

针对现有技术的上述缺点，本发明的目的在于提出一种基于无网格与频率估计的干涉合成孔径雷达相位解缠方法，能够有效避免误差向整个干涉相位图的扩散，并且更加准确的反应目标区域实际地形的变化。

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

一种基于无网格与频率估计的干涉合成孔径雷达相位解缠方法，包括如下步骤：

步骤1，获取干涉合成孔径雷达的干涉相位图；

步骤2，将所述干涉相位图划分为多个子干涉相位图，并计算每个子干涉相位图中相邻像素点之间的缠绕相位梯度值；

步骤3，确定所述每个子干涉相位图的形函数；

步骤4，根据所述每个子干涉相位图中相邻像素点之间的缠绕相位梯度值和所述形函数，得到所述每个子干涉相位图的缠绕相位梯度函数；

步骤5，根据所述每个子干涉相位图的缠绕相位梯度函数，求解所述每个子干涉相位图的解缠相位值；

步骤6，根据所述每个子干涉相位图的解缠相位值，对所述干涉合成孔径雷达的干涉相位图进行解缠。

本发明的特点和进一步的改进为：

(1)所述计算每个子干涉相位图中相邻像素点之间的缠绕相位梯度值，具体包括：

子干涉相位图中任意像素点(x,y)在y方向上和x方向上的缠绕相位梯度值为

其中，表示子干涉相位图中任意像素点(x,y)在y方向上的缠绕相位梯度值，表示子干涉相位图中任意像素点(x,y)在x方向上的缠绕相位梯度值，所述子干涉相位图为任意一个子干涉相位图。

(2)步骤3具体包括如下子步骤：

(3a)将所述子干涉相位图中任意像素点(x,y)的缠绕相位梯度函数表示为：

其中，Q^T(x,y)为给定的基函数，q_i(x,y)是基函数Q^T(x,y)在空间坐标X^T＝[x,y]上给定的单项式，a_i是q_i(x,y)的待定系数，向量a是由a_i组成的列向量；

(3b)根据所述子干涉相位图中每个像素点的缠绕相位梯度值，计算向量a的值，有其中，是由n个像素点的缠绕相位梯度值组成的列向量， $W=\left[\begin{array}{ccccc}1& x_1& y_1& ...& q_n(x_1,y_1)\\ 1& x_2& y_2& ...& q_n(x_2,y_2)\\ & & & ...& \\ 1& x_n& y_n& ...& q_n(x_n,y_n)\end{array}\right],$ W是一个常数矩阵；

(3c)确定子干涉相位图的形函数P^e(x,y)＝Q^T(x,y)W^-1。

(3)步骤4具体为：

子干涉相位图的解缠相位梯度函数φ^e(x,y)、缠绕相位梯度函数在x方向上的分量缠绕相位梯度函数在y方向上的分量表示为：

φ^e(x,y)＝P^e(x,y)φ^e

其中，是由子干涉相位图中n个像素点的解缠相位值φ₁ φ₂ … φ_n组成的列向量，是由子干涉相位图中n个像素点的缠绕相位梯度值在x方向上的分量组成的列向量，是由子干涉相位图中n个像素点的缠绕相位梯度值在y方向上的分量组成的列向量。

(4)步骤5具体包括如下子步骤：

(5a)在子干涉相位图中，构造代价函数L^e(φ^e(x,y))，其中，

$\begin{array}{c}{L}^e\left({\mathrm{\φ}}^e\right(x,y\left)\right)={\mathrm{\φ}}^e(x,y)(\underset{S_e}{\∫\∫}\▿p^e(x,y)\·\▿p^e(x,y\left)dxdy\right){\mathrm{\φ}}^e(x,y)\\ +2{\mathrm{\φ}}^e(x,y){\left(\underset{S_e}{\∫\∫}{p}^e\right(x,y)\·\▿p^e(x,y\left)dxdy\right)}_x\▿{\mathrm{\φ}}_x^e(x,y)\\ +2{\mathrm{\φ}}^e(x,y){\left(\underset{S_e}{\∫\∫}{p}^e\right(x,y)\·\▿p^e(x,y\left)dxdy\right)}_y\▿{\mathrm{\φ}}_y^e(x,y)\end{array}$

其中，表示对子干涉相位图中的形函数p^e(x,y)进行求导运算，S_e表示子干涉相位图中各个像素所组成的闭合区域；

(5b)根据局部频率估计的方法得到

$\begin{array}{c}{L}^e\left({\mathrm{\φ}}^e\right(x,y\left)\right)={\mathrm{\φ}}^e(x,y)(\underset{S_e}{\∫\∫}\▿p^e(x,y)\·\▿p^e(x,y\left)dxdy\right){\mathrm{\φ}}^e(x,y)\\ +2{\mathrm{\φ}}^e(x,y){\left(\underset{S_e}{\∫\∫}{p}^e\right(x,y)\·\▿p^e(x,y\left)dxdy\right)}_x*2\mathrm{\π}*f_x^e(x,y)\\ +2{\mathrm{\φ}}^e(x,y){\left(\underset{S_e}{\∫\∫}{p}^e\right(x,y)\·\▿p^e(x,y\left)dxdy\right)}_y*2\mathrm{\π}*f_y^e(x,y)\end{array}$

其中，表示x方向上的局部频率列向量，表示y方向上的局部频率列向量；

(5c)求所述子干涉相位图中代价函数L^e(φ^e(x,y))的极小值，确定所述子干涉相位图的解缠相位。

本发明与现有的技术相比具有以下优点:

(1)由于本发明是基于最小二乘法的改进算法，所以继承了最小二乘法在求解InSAR解缠相位时的良好稳健性；(2)现有技术利用滤波相位差分法求解缠相位，而此法因为滤波效果好坏存在很大误差；而本发明采用的局部频率估计法则一步计算邻近像素的高度差，因而使用局部频率估计法代替滤波相位差分法求在最小二乘准则下的相位解缠值，不仅能够更加准确而且能够更加快速的反应实际地形的变化；(3)现有技术利用全局最小二乘法准则求解缠相位时残差点处的错误会向整个干涉相位图的扩散，导致结果严重错误；而本发明在解缠过程中引入无网格法，划分局部网格，将误差控制在区域内部，减少了整体误差，提高了解缠精度。

附图说明

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

图1为本发明实施例提供的相位解缠方法的流程示意图。

图2为本发明实施例提供的局部网格的示意图。

图3(a)为模拟真实相位的示意图。

图3(b)为模拟加入噪声后的缠绕相位示意图。

图4(a)为采用枝切法对图3(b)的解缠结果示意图。

图4(b)为最小二乘法对图3(b)的解缠结果示意图。

图4(c)为质量图引导法对图3(b)的解缠结果示意图。

图4(d)为本发明方法对图3(b)的解缠结果示意图。

图5(a)为加入标准差0.6500rad的噪声后采用枝切法对图3(b)的解缠结果的误差分布示意图。

图5(b)为加入标准差0.6500rad的噪声后采用最小二乘法对图3(b)的解缠结果的误差分布示意图。

图5(c)为加入标准差0.6500rad的噪声后采用质量图引导法对图3(b)的解缠结果的误差分布示意图。

图5(d)为加入标准差0.6500rad的噪声后采用本发明方法对图3(b)的解缠结果的误差分布示意图。

图6(a)为加入标准差1.1000rad的噪声后采用枝切法对图3(b)的解缠结果的误差分布示意图。

图6(b)为加入标准差1.1000rad的噪声后采用最小二乘法对图3(b)的解缠结果的误差分布示意图。

图6(c)为加入标准差1.1000rad的噪声后采用质量图引导法对图3(b)的解缠结果的误差分布示意图。

图6(d)为加入标准差1.1000rad的噪声后采用本发明方法对图3(b)的解缠结果的误差分布示意图。

具体实施方式

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

下面结合附图对本发明做进一步的描述。

参照图1，本发明的具体实施步骤如下：

步骤1，获取干涉合成孔径雷达的干涉相位图。

由两幅相干SAR单视图得到一幅具有M*N个像素点的干涉相位图，每个像素点的值代表该像素点的缠绕相位的值。缠绕相位的值与解缠相位的值之间具有如下关系：(x,y)∈S x＝0,1,...,M-1；y＝0,1,...,N-1(1.1)

式(1.1)中，S表示具有M*N个像素点的干涉相位图，(x,y)为干涉相位图中的任意一个像素点，φ(x,y)表示任意一个像素点(x,y)的解缠相位值，表示任意一个像素点(x,y)的缠绕相位值。n(x,y)是整数。从缠绕相位值中求得2π周期的整数倍n(x,y)，从而得到解缠相位值φ(x,y)是整个相位解缠的思路。

步骤2，将所述干涉相位图划分为多个子干涉相位图，并计算每个子干涉相位图中相邻像素点之间的缠绕相位梯度值。

其中，根据预设规则，将所述干涉相位图划分为多个子干涉相位图的方法有多种。示例性的，将干涉相位图划分为多个正方形区域的子干涉相位图，所述正方形区域的子干涉相位图由相邻的四个像素点组成；或者，将干涉相位图划分为多个六边形区域的子干涉相位图；或者，将干涉相位图划分为多个三角形区域的子干涉相位图等等，本发明对此不作限制。

示例性的，本发明实施例以将干涉相位图划分为多个正方形区域的子干涉相位图为例进行说明，所述正方形区域的子干涉相位图由相邻的四个像素点组成的。下面将对任意一个子干涉相位图的解缠过程进行描述，其他子干涉相位图的解缠过程可参考该子干涉相位图的解缠过程。

具体的，将干涉相位图划分成多个如图2所示的局部网格(即子干涉相位图)，计算局部网格中相邻像素点之间的相位差(即缠绕相位梯度值)。

相位梯度值表示相邻像素点的相位的差值。干涉相位图的局部网格中任意一个像素点(x,y)在纵向上和横向上的缠绕相位梯度值可表示如下：

表示干涉相位图的局部网格中的任意一个像素点(x,y)在纵向上的缠绕相位梯度值，表示干涉相位图中的局部网格的任意一个像素点(x,y)在横向上的缠绕相位梯度值(需要说明的是，横向即为x方向，纵向即为y方向)。

步骤3，确定所述每个子干涉相位图的形函数。

需要说明的是，以下所述局部网格即为子干涉相位图。

利用无网格法中的插值法，求得局部网格的形函数，根据求得的形函数将局部网格的缠绕相位梯度函数用相邻像素点之间的缠绕相位梯度值表示，此步骤将局部网格中四个独立的缠绕相位梯度值表示成一个连续的缠绕相位梯度函数。

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

(3a)将子干涉相位图中任意像素点的缠绕相位梯度函数表示为：

其中，Q^T(x,y)为给定的基函数，q_i(x,y)是基函数Q^T(x,y)在空间坐标X^T＝[x,y]上给定的单项式，a_i是q_i(x,y)的待定系数，向量a是由a_i组成的列向量，所述子干涉相位图包括n个像素点。

使用点插值法生成形函数，首先在参考图2划分的子干涉相位图(包含4个像素点)内，假设任意点的缠绕相位梯度函数可以近似表示为：

Q^T(x,y)为给定的基函数，基函数是缠绕相位梯度函数的基底函数，由于缠绕相位梯度函数往往复杂难以表示，而基函数则是易于表达的，因此选取适宜的基函数近似替代缠绕相位梯度函数。基函数的物理意义是提供了一种将缠绕相位梯度函数在各个方向上进行分解的方法，类比于二维空间上，一个矢量可分解成为x方向与y方向上的分量。q_i(x,y)是基函数Q^T(x,y)在空间坐标X^T＝[x,y]上给定的单项式。一般情况下，2维空间q_i(x,y)可由下列形式表示：

Q^T(x,y)＝{1 x y x² xy y² … xⁿ yⁿ}     (3.2)

针对如图2的局部网格，Q^T(x,y)取(3.2)中的前4项，即q₁(x,y)＝1，q₂(x,y)＝x，q₃(x,y)＝y，q₄(x,y)＝x²。a_i是q_i(x,y)的待定系数，4个像素点确定后，待定系数a_i不再发生变化。

(3b)根据所述子干涉相位图中的像素点的缠绕相位值，计算向量a的值，有其中，是由n个像素点的缠绕相位梯度值组成的列向量， $W=\left[\begin{array}{ccccc}1& x_1& y_1& ...& q_n(x_1,y_1)\\ 1& x_2& y_2& ...& q_n(x_2,y_2)\\ & & & ...& \\ 1& x_n& y_n& ...& q_n(x_n,y_n)\end{array}\right]$

W是一个常数矩阵。

示例性的，系数a_i的计算，需要用到子干涉相位图内的4个像素点，系数a_i可由4个像素点缠绕相位梯度值满足式(3.1)而得到，对于4个像素点分别有：

(3.3)写成矩阵形式有：

式子中，是4个像素点缠绕相位梯度值的向量。

$W=\left[\begin{array}{ccccc}1& x_1& y_1& ...& q_4(x_1,y_1)\\ 1& x_2& y_2& ...& q_4(x_2,y_2)\\ 1& x_3& y_3& ...& q_4(x_3,y_3)\\ 1& x_4& y_4& ...& q_4(x_4,y_4)\end{array}\right],$ 当4个像素点确定后，它们的坐标便是已知的，所以W是一个常数矩阵。因此可得：

因为在步骤2已经得到的常数列向量，而W是一个常数矩阵，所以a就是一个常数列向量。

(3c)确定子干涉相位图的形函数P^e(x,y)＝Q^T(x,y)W^-1。

将(3.5)带回(3.1)，得到：

P^e(x,y)＝Q^T(x,y)W^-1     (3.7)

我们称式P^e(x,y)就是局部网格的形函数，它将局部网格的缠绕相位梯度函数用相邻像素点之间的缠绕相位梯度值进行了表示。因为基函数Q^T(x,y)与矩阵W^-1现都已知，故可以顺利的求出局部网格的形函数P^e(x,y)。

步骤4，根据所述每个子干涉相位图中相邻像素点之间的缠绕相位梯度值和所述形函数，得到每个子干涉相位图的缠绕相位梯度函数。

求出形函数P(x,y)后，相似地，局部网格内任意像素点横向的缠绕梯度函数是通该像素点所属的局部网格中的4个像素点的横向缠绕梯度来近似表达的，同理局部网格内任意像素点纵向的缠绕梯度函数是通过该像素点所属的局部网格中的4个像素点的纵向缠绕梯度来近似表达的；而解缠相位函数φ^e(x,y)也可以通过形函数P^e(x,y)与4个像素点的解缠相位值来表示。

由公式(3.6)，可类比的将局部网格中解缠相位梯度函数φ^e(x,y)、缠绕相位梯度函数在x方向上的分量缠绕相位梯度函数在y方向上的分量表示为：

φ^e(x,y)＝P^e(x,y)φ^e       (3.8)

其中：

形函数P^e(x,y)在上述步骤中已经求得，为缠绕相位梯度值在x方向的分量，为缠绕相位梯度值在y方向的分量。

步骤5，根据所述每个子干涉相位图的缠绕相位梯度函数，求解所述每个子干涉相位图的解缠相位值。

其中，子干涉相位图的解缠相位值能够使得所述缠绕相位梯度函数和解缠相位梯度函数的均方误差最小。

根据缠绕相位梯度函数，构造代价函数L(φ(x,y))，

其中，

根据最小二乘法的思路，要求解缠相位梯度函数和缠绕相位梯度函数之差的平方和最小时的解缠相位需要构建代价函数L(φ(x,y))，使得此代价函数L(φ(x,y))最小的解缠相位梯度值即是所需的解缠相位。

是从干涉相位图中可直接观察得到的缠绕相位，因此公式(4.1)的就是已知项，因此只需要求得(4.1)式等号右侧的最小即可，所以把代价函数L(φ(x,y))改写成如下形式。

(5a)在每个子干涉相位图中，代价函数L(φ(x,y))表示为L^e(φ^e(x,y))，具体的，若在局部网格中用L^e(φ^e(x,y))表示代价函数L(φ(x,y))的子泛函，那么根据公式(3.8)，(3.9)和(3.10)可将式(4.2)改写为如下的形式：

$\begin{array}{c}{L}^e\left({\mathrm{\φ}}^e\right(x,y\left)\right)={\mathrm{\φ}}^e(x,y)(\underset{S_e}{\∫\∫}\▿P^e(x,y)\·\▿P^e(x,y\left)dxdy\right){\mathrm{\φ}}^e(x,y)\\ +2{\mathrm{\φ}}^e(x,y){\left(\underset{S_e}{\∫\∫}{P}^e\right(x,y)\·\▿P^e(x,y\left)dxdy\right)}_x\▿{\mathrm{\φ}}_x^e(x,y)\\ +2{\mathrm{\φ}}^e(x,y){\left(\underset{S_e}{\∫\∫}{P}^e\right(x,y)\·\▿P^e(x,y\left)dxdy\right)}_y\▿{\mathrm{\φ}}_y^e(x,y)\end{array}---\left(4.3\right)$

其中，表示对子干涉相位图中的形函数进行求导运算。

(5b)因为得到干涉相位图的回波信号是平稳的，所以根据局部频率估计的知识，可以利用最大似然频率估计的方法获得局部频率值，并且以该频率值的2π倍替代缠绕相位梯度值，则式(4.3)可改写为：

$\begin{array}{c}{L}^e\left({\mathrm{\φ}}^e\right(x,y\left)\right)={\mathrm{\φ}}^e(x,y)(\underset{S_e}{\∫\∫}\▿p^e(x,y)\·\▿p^e(x,y\left)dxdy\right){\mathrm{\φ}}^e(x,y)\\ +2{\mathrm{\φ}}^e(x,y){\left(\underset{S_e}{\∫\∫}{p}^e\right(x,y)\·\▿p^e(x,y\left)dxdy\right)}_x*2\mathrm{\π}*f_x^e(x,y)\\ +2{\mathrm{\φ}}^e(x,y){\left(\underset{S_e}{\∫\∫}{p}^e\right(x,y)\·\▿p^e(x,y\left)dxdy\right)}_y*2\mathrm{\π}*f_y^e(x,y)\end{array}---\left(4.4\right)$

式(4.4)中f_x ^e(x,y)表示横向上的局部频率列向量，f_y ^e(x,y)表示纵向上的局部频率列向量。

(5c)求子干涉相位图中代价函数L^e(φ^e(x,y))的极小值，确定子干涉相位图的解缠相位。

对φ^e求导并且要求即求得(4.4)代价函数L^e的极小值。要使就必须令局部网格的代价函数对四个像素点的导数都为0，即令是局部网格中的四个像素点组成的列向量。

对(4.4)式求导：

$\begin{array}{c}\frac{\∂L^e\left({\mathrm{\φ}}^e\right)}{\∂{\mathrm{\φ}}^e}=\frac{\∂}{\∂{\mathrm{\φ}}^e}{\mathrm{\φ}}^e(\underset{S_e}{\∫\∫}\▿P^e\·\▿P^{e}dxdy){\mathrm{\φ}}^e+\frac{\∂}{\∂{\mathrm{\φ}}^e}2{\mathrm{\φ}}^e{(\underset{S_e}{\∫\∫}{P}^e\·\▿P^{e}dxdy)}_x*2\mathrm{\π}*f_x^e\\ +\frac{\∂}{\∂{\mathrm{\φ}}^e}2{\mathrm{\φ}}^e{(\underset{S_e}{\∫\∫}{P}^e\·\▿P^{e}dxdy)}_y*2\mathrm{\π}*f_y^e\end{array}---\left(4.5\right)$

将上式简化

$\frac{\∂L^e\left({\mathrm{\φ}}^e\right)}{\∂{{\mathrm{\φ}}_1}^e}=\[+4{{\mathrm{\φ}}_1}^e-1{{\mathrm{\φ}}_2}^e-2{{\mathrm{\φ}}_3}^e-1{{\mathrm{\φ}}_4}^e-2\mathrm{\π}*(2f_{x(1,2)}+1f_{x(4,3)}+3f_{y(1,4)}+1f_{y(2,3)})\]$

$\frac{\∂L^e\left({\mathrm{\φ}}^e\right)}{\∂{{\mathrm{\φ}}_1}^e}=\[-1{{\mathrm{\φ}}_1}^e+4{{\mathrm{\φ}}_2}^e-1{{\mathrm{\φ}}_3}^e-2{{\mathrm{\φ}}_4}^e-2\mathrm{\π}*(-2f_{x(1,2)}-1f_{x(4,3)}+1f_{y(1,4)}+2f_{y(2,3)})\]---\left(4.6\right)$

$\frac{\∂L^e\left({\mathrm{\φ}}^e\right)}{\∂{{\mathrm{\φ}}_1}^e}=\[-1{{\mathrm{\φ}}_1}^e-2{{\mathrm{\φ}}_2}^e+4{{\mathrm{\φ}}_3}^e-1{{\mathrm{\φ}}_4}^e-2\mathrm{\π}*(-1f_{x(1,2)}+-2f_{x(4,3)}-1f_{y(1,4)}-2f_{y(2,3)})\]$

$\frac{\∂L^e\left({\mathrm{\φ}}^e\right)}{\∂{{\mathrm{\φ}}_1}^e}=\[-1{{\mathrm{\φ}}_1}^e-2{{\mathrm{\φ}}_2}^e-1{{\mathrm{\φ}}_3}^e+4{{\mathrm{\φ}}_4}^e-2\mathrm{\π}*(1f_{x(1,2)}+2f_{x(4,3)}-2f_{y(1,4)}-1f_{y(2,3)})\]$

因为令所(4.6)继续简写成

$\left|\begin{array}{cccc}4& -1& -2& -1\\ -1& 4& -1& -2\\ -2& -1& 4& -1\\ -1& -2& -1& 4\end{array}\right|\·{\left[\begin{array}{cccc}{\mathrm{\φ}}_1^e& {\mathrm{\φ}}_2^e& {\mathrm{\φ}}_3^e& {\mathrm{\φ}}_4^e\end{array}\right]}^T=2\mathrm{\π}*\left|\begin{array}{cccc}2& 1& 2& -1\\ -2& -1& 1& 2\\ -1& -2& -1& -2\\ 1& 2& -2& -1\end{array}\right|\left[\begin{array}{cccc}{f}_{x(1,2)}& f_{x(4,3)}& f_{x(1,4)}& f_{x(2,3)}\end{array}\right]---\left(4.7\right)$

将(4.7)改换为矩阵的表达式为：

K^eφ^e＝b^e      (4.8)

式(4.6)为对局部网格的代价函数L^e求导的矩阵表达式，K^e为系数矩阵f_x(i,j)((i,j)＝{(1,2),(4,3)}表示局部频率在行向的估计值，f_y(i,j)((i,j)＝{(1,4),(2,3)}表示局部频率在列向上的估计值。(4.8)中，K^e与b^e已知，所以可以顺利求出解缠相位φ^e。

本发明实施例还可以将所有子干涉相位图的代价函数联合求解。

在相位解缠过程中，用各个子泛函的和来表示相位干涉图中所有像素点构成的待解缠域的整体泛函L(φ)：

$L\left(\mathrm{\φ}\right)=\underset{e=1}{\overset{(N-1)(M-1)}{\Σ}}{L}^e\left({\mathrm{\φ}}^e\right)---\left(4.9\right)$

利用驻点条件，对上式进行处理便能得如下方程：

$\frac{\∂L}{\∂\mathrm{\φ}}=\underset{e=1}{\overset{N}{\Σ}}\frac{\∂L^e}{\∂{\mathrm{\φ}}^e}=0---\left(4.10\right)$

分别把单元矩阵K^e与单元向量φ^e、b^e组合成系统(整体区域)矩阵K与系统向量φ，如下表示：

$\frac{\∂L}{\∂\mathrm{\φ}}=\underset{e=1}{\overset{N}{\Σ}}\frac{\∂L^e}{\∂{\mathrm{\φ}}^2}=\underset{e=1}{\overset{N}{\Σ}}(K^e{\mathrm{\φ}}^e-b^e)=0---\left(4.11\right)$

Kφ＝b      (4.12)

式(4.12)是(4.11)的矩阵表达式，是对整体泛函L(φ)求导的矩阵表达式，b为已知条件；K为分块对角阵，用迭代法即可解出整体解缠相位φ。

本发明的效果可以通过下述仿真实验得到验证：

首先调用Matlab中的peaks函数生成128×128相位值，以此作为假设的回波信号的真实相位图。如图3(a)所示，接着，将相干系数为0.9分布区间在[-b,b]，b＝1.1258rad的均匀噪声加在真实相位图3(a)中，得到如图3(b)的干涉相位图(即缠绕相位图)。

图4(a)为采用枝切法对图3(b)的解缠结果、图4(b)为最小二乘法对图3(b)的解缠结果、图4(c)为质量图引导法对图3(b)的解缠结果、图4(d)为本发明方法对图3(b)的解缠结果。比较后可以看出本发明方法具有较好的解缠结果。下一步定量分析枝切法、最小二乘法、质量图引导法与本方法的解缠相位误差的大小以及分布的状况。

图4系列为加入标准差为0.6500rad的噪声后的误差分布图；图5系列为加入标准差为1.1000rad的噪声后的误差分布图。

由以上两组误差的分布图可以准确得到如下的结论：在加入标准差为0.6500rad的噪声后本方法相位解缠的误差分布范围比较广，但是总体的误差值较小；质量图引导图法的误差分布范围最广，枝切法误差分布在残差点密集的区域。在加入标准差为1.1000rad的噪声后最小二乘法和本发明结果一样表现为误差范围分布比较广，且相对本方法最小二乘法在相位变化大的区域的误差值会更大。可以看出，本发明方法与现有的方法相比具有最优的解缠效果。

将图4与图5系列的相位误差进行比较，得到如下表1。

表1

通过表1可以得到结论，在加入不同方差的噪声后本发明方法与现有的最小二乘法、质量图引导法、枝切法相比，具有最小的均方误差，也就证明本发明具有更优的相位解缠效果。分析原因，由于划分支持域限制的了误差的传播，局部频率估计更准确的反映了地形变化情况，有效提高了解缠精度。

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

Claims (5)

1.一种基于无网格与频率估计的干涉合成孔径雷达相位解缠方法，其特征在于，包括如下步骤：

步骤1，获取干涉合成孔径雷达的干涉相位图；

步骤2，将所述干涉相位图划分为多个子干涉相位图，并计算每个子干涉相位图中相邻像素点之间的缠绕相位梯度值；

步骤3，确定所述每个子干涉相位图的形函数；

步骤4，根据所述每个子干涉相位图中相邻像素点之间的缠绕相位梯度值和所述形函数，得到所述每个子干涉相位图的缠绕相位梯度函数；

步骤5，根据所述每个子干涉相位图的缠绕相位梯度函数，求解所述每个子干涉相位图的解缠相位值；

步骤6，根据所述每个子干涉相位图的解缠相位值，对所述干涉合成孔径雷达的干涉相位图进行解缠。

2.根据权利要求1所述的基于无网格与频率估计的干涉合成孔径雷达相位解缠方法，其特征在于，所述计算每个子干涉相位图中相邻像素点之间的缠绕相位梯度值，具体包括：

子干涉相位图中任意像素点(x,y)在y方向上和x方向上的缠绕相位梯度值为

其中，表示子干涉相位图中任意像素点(x,y)在y方向上的缠绕相位梯度值，表示子干涉相位图中任意像素点(x,y)在x方向上的缠绕相位梯度值，所述子干涉相位图为任意一个子干涉相位图。

3.根据权利要求1所述的基于无网格与频率估计的干涉合成孔径雷达相位解缠方法，其特征在于，步骤3具体包括如下子步骤：

(3a)将所述子干涉相位图中任意像素点(x,y)的缠绕相位梯度函数表示为：

其中，Q^T(x,y)为给定的基函数，q_i(x,y)是基函数Q^T(x,y)在空间坐标X^T＝[x,y]上给定的单项式，a_i是q_i(x,y)的待定系数，向量a是由a_i组成的列向量，所述子干涉相位图中包括n个像素点；

(3b)根据所述子干涉相位图中每个像素点的缠绕相位梯度值，计算向量a的值，有其中，是由n个像素点的缠绕相位梯度值组成的列向量， $W=\left[\begin{array}{ccccc}1& x_1& y_1& ...& q_n(x_1,y_1)\\ 1& x_2& y_2& ...& q_n(x_2,y_2)\\ & & & ...& \\ 1& x_n& y_n& ...& q_n(x_n,y_n)\end{array}\right],$ W是一个常数矩阵；

(3c)确定子干涉相位图的形函数P^e(x,y)＝Q^T(x,y)W^-1。

4.根据权利要求1所述的基于无网格与频率估计的干涉合成孔径雷达相位解缠方法，其特征在于，步骤4具体为：

子干涉相位图的解缠相位梯度函数φ^e(x,y)、缠绕相位梯度函数在x方向上的分量缠绕相位梯度函数在y方向上的分量表示为：

φ^e(x,y)＝P^e(x,y)φ^e

其中，φ^e＝(φ₁ φ₂ … φ_n)^T是由子干涉相位图中n个像素点的解缠相位值φ₁φ₂…φ_n组成的列向量，是由子干涉相位图中n个像素点的缠绕相位梯度值在x方向上的分量组成的列向量，是由子干涉相位图中n个像素点的缠绕相位梯度值在y方向上的分量组成的列向量。

5.根据权利要求1所述的基于无网格与频率估计的干涉合成孔径雷达相位解缠方法，其特征在于，步骤5具体包括如下子步骤：

(5a)在子干涉相位图中，构造代价函数L^e(φ^e(x,y))，其中，

$\begin{array}{c}{L}^e\left({\mathrm{\φ}}^e\right(x,y\left)\right)={\mathrm{\φ}}^e(x,y)(\underset{S_e}{\∫\∫}\▿p^e(x,y)\·\▿p^e(x,y\left)dxdy\right){\mathrm{\φ}}^e(x,y)\\ +2{\mathrm{\φ}}^e(x,y){\left(\underset{S_e}{\∫\∫}{p}^e\right(x,y)\·\▿p^e(x,y\left)dxdy\right)}_x\▿{\mathrm{\φ}}_x^e(x,y)\\ +2{\mathrm{\φ}}^e(x,y){\left(\underset{S_e}{\∫\∫}{p}^e\right(x,y)\·\▿p^e(x,y\left)dxdy\right)}_y\▿{\mathrm{\φ}}_y^e(x,y)\end{array}$

其中，表示对子干涉相位图中的形函数p^e(x,y)进行求导运算，S_e表示子干涉相位图中各个像素所组成的闭合区域；

(5b)根据局部频率估计的方法得到

$\begin{array}{c}{L}^e\left({\mathrm{\φ}}^e\right(x,y\left)\right)={\mathrm{\φ}}^e(x,y)(\underset{S_e}{\∫\∫}\▿p^e(x,y)\·\▿p^e(x,y\left)dxdy\right){\mathrm{\φ}}^e(x,y)\\ +2{\mathrm{\φ}}^e(x,y){\left(\underset{S_e}{\∫\∫}{p}^e\right(x,y)\·\▿p^e(x,y\left)dxdy\right)}_x*2\mathrm{\π}*f_x^e(x,y)\\ +2{\mathrm{\φ}}^e(x,y){\left(\underset{S_e}{\∫\∫}{p}^e\right(x,y)\·\▿p^e(x,y\left)dxdy\right)}_y*2\mathrm{\π}*f_y^e(x,y)\end{array}$

其中，f_x ^e(x,y)表示x方向上的局部频率列向量，f_y ^e(x,y)表示y方向上的局部频率列向量；

(5c)求所述子干涉相位图中代价函数L^e(φ^e(x,y))的极小值，确定所述子干涉相位图的解缠相位。