CN105158761A - 基于枝切法和曲面拟合的雷达合成相位解缠方法 - Google Patents

Abstract

本发明公开了一种基于枝切法和曲面拟合的雷达相位解缠方法，其主要思路是：SAR图像经过图像配准后的图像数据与该SAR图像经过图像配准后的图像数据共轭相乘，得到干涉相位图，进而得到干涉相位图的伪相关系数；设定所述干涉相位图的伪相关系数阈值，得到第一质量区域图和第二质量区域图，利用枝切法对第一质量区域图进行相位解缠，得到第一质量区域图的解缠相位，并采用多项式曲面拟合法对干涉相位图的解缠相位进行拟合，得到干涉相位图的多项式拟合函数的系数列向量和存在误差的干涉相位图的解缠相位曲面，采用遗传算法对存在误差的干涉相位图的解缠相位曲面进行优化，得到优化后的干涉相位图的解缠相位曲面，即干涉相位图的最终解缠相位。

Description

基于枝切法和曲面拟合的雷达合成相位解缠方法

技术领域

本发明属于干涉合成孔径雷达(InSAR)相位解缠领域，进一步涉及一种基于枝切法和曲面拟合的雷达合成相位解缠方法。

背景技术

干涉合成孔径雷达(InterferometricSyntheticApertureRadar，InSAR)雷达技术首先将目标生成两幅合成孔径雷达(SyntheticApertureRadar，SAR)图像，然后经过图像预处理、相位解缠得到高精度的数字高程图。InSAR雷达技术作为一种新型的空间对地观测技术，是SAR雷达技术与射电天文学的完美结合，它不仅能够准确高效地识别目标的数字高程图，而且还能够全天候、全天时工作，上述优点使得InSAR雷达技术成为21世纪新的研究热点。

SAR雷达图像配准后生成的两幅单视复图像(SingleLookComplex，SLC)图像共轭相乘后提取相位信息，即可得到干涉相位图像，对干涉相位图像的相位解缠是InSAR雷达技术中的关键环节，相位解缠结果的精度直接影响数字高程图的准确度。若相位解缠方法选择不当，直接会导致干涉相位图像的相位无法解缠或相位解缠结果错误。因此，对相位解缠算法的研究至关重要。

InSAR相位解缠的算法根据积分方式的不同可以分为局部法和全局法。局部法解缠速度快且精度高，但该算法在干涉相位图像的相位残差点分布较密集区域解缠效果差，甚至会出现无法解缠的情况；局部法以枝切(BranchCut)法为代表，该算法通过最近邻原则连接干涉相位图像的相位留数点，形成枝切线来隔离噪声，然后对干涉相位图像的相邻像元的差分相位进行积分来实现相位解缠。

但由于枝切法只利用了干涉相位图像的相位留数信息，使得使用该算法枝切线会形成闭环回路或枝切线放置位置不当而导致噪声传递的现象。全局法对干涉相位图像的相位残差点分布密集的区域也能进行有效的解缠，但该算法会将其解缠误差传递到非噪声区域，从而造成相位解缠精度的下降；并且，全局法以最小二乘法为代表，该算法利用相位解缠梯度与缠绕相位梯度的最小均方差实现相位解缠，能够提高相位解缠精度，稳定性增强，但是该算法计算量大，实现起来需要大量的时间。

本发明人随着对InSAR雷达相位解缠算法的深入研究发现，将多种算法结合起来对干涉相位图像进行相位解缠的结果要优于采用单一算法得到的相位解缠结果。

发明内容

针对以上现有技术存在的不足，本发明的目的在于提出一种基于枝切法和曲面拟合的雷达相位解缠方法，该方法将枝切法与多项式曲面拟合法相结合，既发挥了枝切法在干涉相位图像中残差点少的区域解缠精度高、解缠速度快的优势，又利用了多项式曲面拟合法平稳的优点，能够克服干涉相位图像中残差点少引起的数据质量差区域解缠精度低或无法解缠的困难。

一种基于枝切法和曲面拟合的雷达相位解缠方法，其特征在于，包括以下步骤：

步骤1，将SAR图像经过图像配准后的图像数据与该SAR图像经过图像配准后的图像数据共轭相乘，得到干涉相位图，并得到所述干涉相位图的伪相关系数；设定所述干涉相位图的伪相关系数的阈值，进而将干涉相位图分为第一质量区域图和第二质量区域图；

步骤2，利用枝切法对第一质量区域图进行相位解缠，得到第一质量区域图的解缠相位φ'_a,b；其中，(a,b)表示第一质量区域图中像素点的坐标；

步骤3，利用第一质量区域图的解缠相位φ'_a,b，采用多项式曲面拟合法对干涉相位图的解缠相位进行拟合，得到干涉相位图的n阶多项式拟合函数f(x,y)的系数列向量A和存在误差的干涉相位图的解缠相位曲面；

步骤4，根据n阶多项式拟合函数f(x,y)的系数列向量A，并采用遗传算法对得到的存在误差的干涉相位图的解缠相位曲面进行优化，得到优化后的干涉相位图的解缠相位曲面，即干涉相位图的最终解缠相位。

本发明相对于现有技术主要优点在于：

第一，本发明的创新点在于将枝切法和曲面拟合的方法相结合，并应用于雷达相位解缠领域中。

第二，本发明在对干涉相位图进行相位解缠时具有较高的相位解缠精度；

第三，由于本发明在残差点少的区域利用枝切法进行相位解缠，故本发明继承了枝切法解缠精度高、解缠速度快的优点；

第四，由于本发明在数据质量差的区域采用的多项式曲面拟合法进行相位解缠，因此本发明也继承了多项式曲面拟合法稳健性的优点。

附图说明

下面结合附图和具体实施方式对本发明作进一步详细说明。

图1是本发明的一种基于枝切法和曲面拟合的雷达相位解缠方法的流程示意图；

图2A是无噪声的真实相位示意图；

图2B是加入噪声的真实相位示意图；

图3是加入噪声的缠绕相位曲面示意图；

图4是伪相关系数图；

图5是统计直方图；

图6是利用伪相关系数生成的质量图；

图7是叠加后的质量图；

图8是高质量区域解缠的相位示意图；

图9是完整解缠的相位示意图；

图10A是Goldstein枝切法解缠的相位图；

图10B是最小二乘法解缠的相位图；

图10C是质量图引导法解缠的相位图；

图10D是本发明解缠的相位图；

图11A是Goldstein枝切法解缠的误差分布图；

图11B是最小二乘法解缠的误差分布图；

图11C是质量图引导法解缠的误差分布图；

图11D是本发明解缠的误差分布图；

具体实施方式

参照图1，为本发明的一种基于枝切法和曲面拟合的雷达相位解缠方法的流程示意图，该种基于枝切法和曲面拟合的雷达相位解缠方法，包括以下步骤：

步骤1，将SAR图像经过图像配准后的图像数据与该SAR图像经过图像配准后的图像数据共轭相乘，得到干涉相位图，并得到所述干涉相位图的伪相关系数；设定所述干涉相位图的伪相关系数的阈值，进而将干涉相位图分为第一质量区域图和第二质量区域图；

具体地，SAR图像经过图像配准后的图像数据与该SAR图像经过图像配准后的图像数据共轭相乘，得到干涉相位图。为了提高相位解缠后的干涉相位图的相位精度，需要对干涉相位图进行质量区域划分。常用质量图主要包括伪相关系数图和相位梯度变化图：

1a)SAR图像经过图像配准后的图像数据与该SAR图像经过图像配准后的图像数据共轭相乘，得到干涉相位图，并得到所述干涉相位图的伪相关系数；设定所述干涉相位图的伪相关系数的阈值，进而将干涉相位图分为第一质量区域图和第二质量区域图。

伪相关系数是干涉相位图的质量划分标准，利用伪相关系数得到伪相关系数图，伪相关系数可定义为|Z_r,s|，其表达式为：

其中，设干涉相位图的中心像素点坐标为(r,s)，r表示干涉相位图中心像素点的横坐标，且r∈[0,M-1]，M表示干涉相位图单行像素点的总个数；s表示干涉相位图中心像素点的纵坐标，且s∈[0,N-1]，N表示干涉相位图单列像素点的总个数，i表示所述干涉相位图中心像素点(r,s)k阶邻域内的任意一像素点的横坐标，j示所述干涉相位图中心像素点(r,s)k阶邻域内的任意一像素点的纵坐标，∑表示求和符号，且该∑的范围是所述干涉相位图中心像素点(r,s)的k阶邻域，表示中心像素点(r,s)的k阶邻域范围内像素点(i,j)的缠绕相位值，cos表示求余弦，sin表示求正弦。

取干涉相位图中心像素点的坐标为(r,s)，并得到该中心像素点(r,s)的k阶邻域范围内所有待解缠相位的正弦和与余弦和的平方值，该平方值除以k²后得到的结果，即为所述伪相关系数图中心像素点(r,s)的伪相关系数。对干涉相位图中每一个像素点进行计算，即可求得干涉相位图中所有像素点的伪相关系数，伪相关系数越大，表示干涉相位图的数据质量越好。

然后利用干涉相位图的伪相关系数和设定的该干涉相位图的伪相关系数的阈值，对干涉相位图进行质量区域划分，得到第一质量区域图和第二质量区域图。

1b)SAR图像经过图像配准后的图像数据与该SAR图像经过图像配准后的图像数据共轭相乘，得到干涉相位图，并得到所述干涉相位图的相位梯度变化系数设定所述干涉相位图的相位梯度变化系数的阈值Q₂，进而将干涉相位图分为第一质量区域图和第二质量区域图。

相位梯度变化系数是干涉相位图进行质量划分的另一个标准，利用相位梯度变化系数得到相位梯度变化图，相位梯度变化系数可定义为其表达式为：

$\left|{\stackrel{\‾}{Z}}_{r,s}\right|=\frac{\sqrt{\Σ{\left({\mathrm{\Δ}}_{i,j}^x-{\stackrel{\‾}{\mathrm{\Δ}}}_{r,s}^x\right)}^2}+\sqrt{\Σ{\left({\mathrm{\Δ}}_{i,j}^y-{\stackrel{\‾}{\mathrm{\Δ}}}_{r,s}^y\right)}^2}}{k^2}---\left(2\right)$

其中，干涉相位图的中心像素点坐标为(r,s)，r表示干涉相位图中心像素点的横坐标，且r∈[0,M-1]，M表示干涉相位图单行像素点的总个数；s表示干涉相位图中心像素点的纵坐标，且s∈[0,N-1]，N表示干涉相位图单列像素点的总个数，i表示所述干涉相位图中心像素点(r,s)k阶邻域内的任意一像素点的横坐标，j表示所述干涉相位图中心像素点(r,s)k阶邻域内的任意一像素点的纵坐标，∑表示求和符号，且该∑的范围是所述干涉相位图中心像素点(r,s)的k阶邻域，表示所述相位梯度变化图横坐标方向的相位缠绕梯度，表示所述相位梯度变化图纵坐标方向的相位缠绕梯度，表示k阶邻域内的平均值，表示k阶邻域内的平均值。

相位梯度变化图与伪相关系数图不同，相位梯度越大，表示干涉相位图的数据质量越差；相反，若相位梯度越小，则表示干涉相位图的数据质量越好。利用相位梯度变化图和设定的相位梯度变化图阈值，对干涉相位图进行质量区域划分，得到第一质量区域图和第二质量区域图。

比较伪相关系数图和相位梯度变化图，发现该两类方法均可以对干涉相位图进行有效的质量区域划分，且效果相似。本发明以伪相关系数图为划分标准。

步骤2，利用枝切法对第一质量区域图进行相位解缠，得到第一质量区域图的解缠相位φ'_a,b；其中，(a,b)表示第一质量区域图中像素点的坐标。

步骤2的具体子步骤为：

2a)根据残差点的计算公式

$\underset{c=1}{\overset{4}{\Σ}}\widehat{\▿}\mathrm{\φ}\left(c\right)=2\mathrm{\π}q$

得到第一质量区域图中的正残差点坐标和负残差点坐标，并根据该正残差点坐标和负残差点坐标分别将正残差点和负残差点的位置在第一质量区域图中标注出来，得到第一质量区域图的正负残差点分布图。

具体的，设第一质量区域图中任意一个像素点的坐标为(a,b)，表示该像素点(a,b)与该像素点(a,b)的第c个相邻像素点之差；那么残差点就可以定义为：当第一质量区域图中像素点(a,b)的相邻四个像素点的中心旋度均不等于0时，就可将该像素点(a,b)为定义为残差点；然后根据式(3)得到q值，当q为1时，该残差点(a,b)表示为正残差点；当q为-1时，该残差点(a,b)表示为负残差点；

当第一质量区域图中像素点(a,b)相邻的四个像素点的中心旋度均等于0时，即q取0时，该像素点(a,b)定义为非残差点。残差点的计算公式为

$\underset{c=1}{\overset{4}{\Σ}}\widehat{\▿}\mathrm{\φ}\left(c\right)=2\mathrm{\π}q---\left(3\right)$

其中，表示像素点(a,b)与所述像素点(a,b)的第c个相邻像素点之差，c表示第一质量区域图中像素点(a,b)的第c个像素点，c∈{1,2,3,4}，q表示判决常数。

根据残差点的计算公式，得到第一质量区域图中的正残差点坐标和负残差点坐标，并根据该正残差点坐标和负残差点坐标分别将正残差点和负残差点的位置在第一质量区域图中标注出来，得到第一质量区域图的正负残差点分布图。

2b)选取第一质量区域图的正负残差点分布图中的第一个正残差点为中心，并在搜索半径为r×r的范围内进行搜索，若搜到负残差点，则将该正残差点与搜到的负残差点定义为近点对；相反，若没有搜到负残差点，则将该正残差点定义为远点；其中，r×r表示设定的初始搜索半径；

若该正残差点是远点，则以该远点为中心，在第一质量区域图的正负残差点分布图中依次在搜索半径为r×r，(r+2)×(r+2)，……的搜索范围内再次进行搜索，直到搜到负残差点，并与该远点组成远点对；若搜索到第一质量区域图的正负残差点分布图边界处仍未搜到负残差点，则将第一质量区域图的正负残差点分布图边界点作为负残差点，并与该远点组成远点对，进而得到第一质量区域图的正负残差点分布图中的第一个远点对；其中，r×r表示设定的初始搜索半径，本发明中初始搜索半径为5×5；

逐个选取第一质量区域图的正负残差点分布图中的每一个正残差点重复此过程，进而得到第一质量区域图的正负残差点分布图中的P个近点对和Q个远点对；其中，P表示第一质量区域图的正负残差点分布图中的近点对个数，Q表示第一质量区域图的正负残差点分布图中的远点对个数；

2c)选取第一质量区域图的正负残差点分布图中的第一个远点对的正残差点为中心、以设置的枝切线长为半径的范围内进行搜索，得到所述范围内的C₁个远点对和D₁个近点对，并将所述范围内的C₁个远点对和D₁个近点对进行重新组合，得到第一质量区域图的正负残差点分布图中的第一个远点对的正残差点对应范围内重组后的E₁个近点对和重组后的F₁个远点对，进而得到第一质量区域图的正负残差点分布图中的第一个远点对的正残差点对应范围内的G₁个优化后的枝切线。

逐个选取第一质量区域图的正负残差点分布图中每一个远点对的正残差点重复此过程，直到得到第一质量区域图的正负残差点分布图中的Q个远点对的正残差点各自对应范围内重组后的E₁～E_Q个近点对和重组后的F₁～F_Q个远点对，进而得到第一质量区域图的正负残差点分布图中的Q个远点对的正残差点各自对应范围内的G₁～G_Q个优化后的枝切线。

具体地，选取第一质量区域图的正负残差点分布图中的第q个远点对的正残差点为中心、以设置的枝切线长为半径的范围并进行搜索，得到该范围内的C_q个远点对和D_q个近点对，并将该范围内的C_q个远点对和D_q个近点对进行重新组合，然后判断根据式(4)计算得到的重新组合后评价标准值是否比重新组合前的评价标准值减小：若重新组合后评价标准值减小，则得到第一质量区域图的正负残差点分布图中的第q个远点对的正残差点对应范围内重组后的E_q个近点对和重组后的F_q个远点对，进而得到第一质量区域图的正负残差点分布图中的第q个远点对的正残差点对应范围内的G_q个优化后的枝切线；若重新组合后评价标准值没有减小，则将该范围内的C_q个远点对和D_q个近点对再次进行重新组合，直到所述重新组合后得到式(4)的评价标准值减小为止。

$fitness=\underset{f_q=1}{\overset{C_q+D_q}{\Σ}}{\[{(x_{f_q}^{+}-x_{f_q}^{-})}^2+{(y_{f_q}^{+}-y_{f_q}^{-})}^2\]}^{1/2}---\left(4\right)$

其中，fitness表示评价标准值，即表示在第一质量区域图的正负残差点分布图中，该远点对在以设置的枝切线长为半径的范围内搜索到的所有远点对的距离和所有近点对的距离之和，表示第一质量区域图的正负残差点分布图中的第q个远点对的正残差点对应范围内的第f个正残差点的横坐标，表示第一质量区域图的正负残差点分布图中的第q个远点对的正残差点对应范围内的第f个负残差点的横坐标，表示第一质量区域图的正负残差点分布图中的第q个远点对的正残差点对应范围内的第f个正残差点的纵坐标、表示第一质量区域图的正负残差点分布图中的第q个远点对的正残差点对应范围内的第f个负残差点的纵坐标，C_q表示第一质量区域图的正负残差点分布图中的第q个远点对的正残差点为中心、以设置的枝切线长为半径的范围并进行搜索，得到该范围内远点对的个数，D_q表示第一质量区域图的正负残差点分布图中的第q个远点对的正残差点为中心、以设置的枝切线长为半径的范围并进行搜索，得到该范围内近点对的个数，q∈{1,2,…,Q}，Q表示第一质量区域图的正负残差点分布图中的远点对个数。

逐个选取第一质量区域图的正负残差点分布图中每一个远点对的正残差点重复此过程，直到得到第一质量区域图的正负残差点分布图中的Q个远点对的正残差点各自对应范围内重组后的E₁～E_Q个近点对和重组后的F₁～F_Q个远点对，进而得到第一质量区域图的正负残差点分布图中的Q个远点对的正残差点各自对应范围内的G₁～G_Q个优化后的枝切线。

随着评价标准值fitness的不断减小，设置的枝切线长度也随之变小。

2d)根据第一质量区域图的正负残差点分布图中的Q个远点对的正残差点各自对应范围内的G₁～G_Q个优化后的枝切线，对第一质量区域图中的缠绕相位梯度进行积分，得到第一质量区域图的解缠相位φ'_a,b；其中，(a,b)表示第一质量区域图中像素点的坐标。

步骤3，利用第一质量区域图的解缠相位φ'_a,b，采用多项式曲面拟合法对干涉相位图的解缠相位进行拟合，得到干涉相位图的n阶多项式拟合函数f(x,y)的系数列向量A和存在误差的相位干涉图的解缠相位曲面。

具体地，多项式曲面拟合法是指利用已知数据选用多项式对待求曲面进行拟合估计。在相位解缠技术中，利用解缠后的第一质量区域图的解缠相位φ'_a,b选用多项式函数对干涉相位图进行相位曲面估计，其过程为：

首先定义干涉相位图的n阶多项式拟合函数f(x,y)，用α₀,α₁,α₂……α_m表示干涉相位图的n阶多项式拟合函数的系数，则干涉相位图的n阶多项式拟合函数f(x,y)可以表示为：

f(x,y)＝α₀+α₁x+α₂y+α₃x²+α₄xy+α₅y²+...+α_myⁿ(5)

干涉相位图的n阶多项式拟合函数f(x,y)的最小均方误差E可以用式(6)表示：

$E=\underset{s_1}{\Σ}{\[{\mathrm{\φ}}_{a,b}^{\′}-f(x_{a,b},y_{a,b})\]}^2---\left(6\right)$

其中，φ'_a,b表示解缠后的第一质量区域图中像素点(a,b)的解缠相位，干涉相位图的n阶多项式拟合函数f(x,y)在第一质量区域图中像素点(a,b)处的值用f(_xa,b,y_a,b)表示，x_a,b表示第一质量区域图中像素点(a,b)在x轴方向的坐标，y_a,b表示第一质量区域图中像素点(a,b)在y轴方向的坐标，s₁表示第一质量区域图中的所有像素点，∑表示求和符号，该的范围是第一质量区域图，所述干涉相位图的n阶多项式拟合函数f(x,y)的第m个系数α_m的下标m和所述干涉相位图的n阶多项式拟合函数f(x,y)的阶数n的关系为：

m>n(7)

由于解缠后的第一质量区域图的解缠相位φ'_a,b是由步骤2得到的，而干涉相位图的n阶多项式拟合函数f(x,y)的系数是未知的，故式(7)可以转换为求解干涉相位图的n阶多项式拟合函数的系数的问题。

为使相位干涉图的n阶多项式拟合函数f(x,y)的最小均方误差E最小，可令其一阶导数为0，推导过程如式(8)所示：

$\begin{array}{c}\frac{dE}{{\mathrm{d\α}}_0}=2\underset{s_1}{\Σ}\[{\mathrm{\φ}}_{a,b}^{\′}-f(x_{a,b},y_{a,b})\]\frac{df}{{\mathrm{d\α}}_0}=0\\ \frac{dE}{{\mathrm{d\α}}_1}=2\underset{s_1}{\Σ}\[{\mathrm{\φ}}_{a,b}^{\′}-f(x_{a,b},y_{a,b})\]\frac{df}{{\mathrm{d\α}}_1}=0\\ \frac{dE}{{\mathrm{d\α}}_2}=2\underset{s_1}{\Σ}\[{\mathrm{\φ}}_{a,b}^{\′}-f(x_{a,b},y_{a,b})\]\frac{df}{{\mathrm{d\α}}_2}=0\\ .\\ .\\ .\\ \frac{dE}{{\mathrm{d\α}}_m}=2\underset{s_1}{\Σ}\[{\mathrm{\φ}}_{a,b}^{\′}-f(x_{a,b},y_{a,b})\]\frac{df}{{\mathrm{d\α}}_m}=0\end{array}---\left(8\right)$

对式(8)进行求解，得相位干涉图的n阶多项式拟合函数f(x,y)的系数列向量A。

根据相位干涉图的n阶多项式拟合函数f(x,y)的系数列向量A，得到存在误差的干涉相位图的解缠相位曲面，但此时得到的干涉相位图的解缠相位曲面误差很大。

步骤4，根据n阶多项式拟合函数f(x,y)的系数列向量A，并采用遗传算法对得到的存在误差的干涉相位图的解缠相位曲面进行优化，得到优化后的干涉相位图的解缠相位曲面，即干涉相位图的最终解缠相位。

具体地，为提高相位解缠精度，本发明采用遗传算法对干涉相位图的n阶多项式拟合函数的系数进行优化，即缩小缠绕方向的梯度与解缠相位梯度差，步骤4的具体子步骤为：

4a)根据干涉相位图的n阶多项式拟合函数f(x,y)的系数列向量A，得到初始种群P(0)。

具体地，利用遗传算法将干涉相位图的n阶多项式拟合函数f(x,y)的系数列向量A，即α₀,α₁,α₂……α_m作为一条染色体；其中，α₀,α₁,α₂……α_m分别对应该条染色体上的0…m号基因。

遗传算法需要一组初始种群P(0)，而该初始种群P(0)种群中的每一条染色体对应一种解，也就是每一条染色体对应的多项式系数代表了相应的相位干涉图的n阶多项式拟合函数f(x,y)。

初始种群P(0)中含有H₀条染色体，设第h条染色体上的0…m号基因分别表示为α_0,h,α_1,h,α_2,h……α_m,h，多项式拟合的系数α₀,α₁,α₂……α_m为第一条染色体，即α_0,1,α_1,1,α_2,1……α_m,1是已知的。初始种群中除了第一条染色体α_0,1,α_1,1,α_2,1……α_m,1以外，其余的染色体将根据式(9)得到：

α_d,h＝α_d,1+SmallNum(9)

其中，SmallNum是一个随机数，具体的大小将根据第一条染色体基因值α_0,1,α_1,1,α_2,1……α_m,1的大小而调整，α_d,h表示第h染色体的第d号基因值大小，且d∈[0,m]，h∈[2,H₀]。

4b)设定最大进化代数为T，并根据评价函数ε²的计算公式

${\mathrm{\ϵ}}^2=\left|\right|\▿\mathrm{\Φ}-\▿\mathrm{\Ψ}|{|}^2,$

和初始种群P(0)，得到下一代种群P(t)；其中，t表示迭代次数，t的初始值为1；其中，表示存在误差的干涉相位图的解缠相位曲面的缠绕相位梯度列向量，表示干涉相位图的n阶多项式拟合函数f(x,y)中解缠相位梯度列向量。

具体地，为了检验染色体的最优性，评价函数ε²如式(10)所示

${\mathrm{\ϵ}}^2=\left|\right|\▿\mathrm{\Φ}-\▿\mathrm{\Ψ}|{|}^2---\left(10\right)$

上式中，表示存在误差的干涉相位图的解缠相位曲面的缠绕相位梯度列向量，表示干涉相位图的n阶多项式拟合函数f(x,y)中解缠相位梯度列向量。根据多项式拟合曲面的均方误差ε²对最优染色体进行保留。

选择策略是选择合适的染色体保留到下一代，是遗传算法中的关键环节。根据H_t-1条染色体的评价函数ε²值的大小，在上一代种群P(t-1)中选择k_t-1条最优的染色体，即对应k_t-1种解，并保留到下一代；

交叉操作就是在k_t-1条最优的染色体随机选择两条染色体，并对该两条染色体相应位置的基因进行交叉操作，计算并比较交叉操作前后的均方误差ε²的值，若均方误差ε²的值减小，则保留交叉操作结果，得到经过交叉操作的最优的染色体；

变异操作是在经过交叉操作的最优的染色体中随机选择一条染色体中的一个基因，使其根据式(9)进行变异，比较变异前后的均方误差，若均方误差减小，则保留变异结果，得到下一代种群P(t)；其中，t表示迭代次数，t的初始值为1。

对初始种群P(0)依次进行选择操作、交叉操作、变异操作，得到下一代种群P(t)；其中，t表示迭代次数，t的初始值为1；

4c)将迭代次数t的值加1，重复子步骤4b)，直到迭代次数t的值等于设置的最大进化代数T，进而得到优化后的相位干涉图的n阶多项式拟合函数f(x,y)的系数。

4d)根据优化后的n阶多项式拟合函数f(x,y)的系数，得到优化后的干涉相位图的解缠相位曲面，即得到干涉相位图的最终解缠相位。

下面结合仿真实验，对本发明的效果做进一步验证。

(一)实验条件

为了进一步说明本发明方法较其它相位解缠方法的优越性，做如下仿真实验；本实验对真实相位曲面为128×128的模拟山峰进行仿真，并加入信噪比为0.0000005dB的噪声

(二)实验内容

本实验对真实相位曲面为128×128的模拟山峰进行仿真，并加入信噪比为0.0000005dB的噪声。

图2A是无噪声的真实相位示意图，对图2A加入信噪比为0.0000005dB的噪声，得到结果如图2B所示，对图2B进行缠绕处理后得到加入噪声的缠绕相位曲面示意图如图3所示，利用伪相关系数划分干涉相位图进行质量区域划分，取邻域阶数k＝3，求得伪相关系数图如图4所示。

参照图5，图5中的黑色小框表示伪相关系数在X:0.2266处的像素点总个数Y最多，为13个，四舍五入后0.2266取为0.23；分析图5可以发现，伪相关系数在0.23处分布最密集，故将干涉相位图的伪相关系数的阈值设为0.23，根据阈值划分干涉相位图，得到通过伪相关系数生成的质量图如图6所示。质量图可以通过伪相关系数低于阀值0.23的区域以及高于干涉相位图的伪相关系数的阈值0.23区域的叠加而得到，如图7所示，其中白色部分表示第一质量区域，黑色部分表示第二质量区域。

利用枝切法对质量图显示为第一质量的区域进行相位解缠的仿真结果如图8所示。利用多项式拟合法对整个相位曲面进行拟合，最终得到的解缠相位的仿真结果如图9所示。

分别利用Goldstein枝切法、最小二乘法、质量图引导法实相位曲面为128×128的模拟山峰进行相位解缠最终得到的解缠相位的仿真结果分别如图10A、图10B、图10C、图10D所示；其中，图10A是Goldstein枝切法解缠的相位图；图10B是最小二乘法解缠的相位图；图10C是质量图引导法解缠的相位图；图10D是本发明解缠的相位图；(三)实验结果分析

对比图9、图2A和图2B可以发现，解缠相位与真实相位相似度高，说明本发明不仅能有效的进行相位解缠，而且相位解缠的准确度高。

对比图10A、图10B、图10C、图10D，即使用本发明方法分别与使用Goldstein枝切法、最小二乘法、质量图引导法对相位解缠进行对比可以发现，使用本发明方法得到的解缠相位与真实相位最为接近，说明了本发明方法的优越性。

图11A是Goldstein枝切法解缠的误差分布图；图11B是最小二乘法解缠的误差分布图；图11C是质量图引导法解缠的误差分布图；图11D是本发明解缠的误差分布图。

对比图11A、图11B、图11C、图11D，即使用本发明方法得到的误差分布分别与使用Goldstein枝切法得到的误差分布、使用最小二乘法得到的误差分布、使用质量图引导法得到的误差分布进行对比可以发现，本发明的误差分布较为均匀，且误差值偏小，说明本发明方法对相位进行解缠的准确度高。

表1是分别使用Goldstein枝切法、最小二乘法、质量图引导法以及本发明方法进行相位解缠后得到的均方误差对比关系表。

表1

由表1可以看出，首先，本发明的雷达相位解缠方法没有未解缠区域，说明本发明方法解缠全面彻底，解缠区域覆盖整个干涉相位图；其次，本发明的雷达相位解缠方法的解缠相位的均方误差较小，说明本发明方法的解缠准确度高，提高了干涉相位图解缠的准确度。

综上所述，仿真实验验证了本发明的正确性，有效性和可靠性。

显然，本领域的技术人员可以对本发明进行各种改动和变型而不脱离本发明的精神和范围；这样，倘若本发明的这些修改和变型属于本发明权利要求及其等同技术的范围之内，则本发明也意图包含这些改动和变型在内。

Claims (6)

1.一种基于枝切法和曲面拟合的雷达相位解缠方法，其特征在于，包括以下步骤：

步骤1，将SAR图像经过图像配准后的图像数据与该SAR图像经过图像配准后的图像数据共轭相乘，得到干涉相位图，并得到所述干涉相位图的伪相关系数；设定所述干涉相位图的伪相关系数的阈值，进而将干涉相位图分为第一质量区域图和第二质量区域图；

步骤2，利用枝切法对第一质量区域图进行相位解缠，得到第一质量区域图的解缠相位φ'_a,b；其中，(a,b)表示第一质量区域图中像素点的坐标；

步骤3，利用第一质量区域图的解缠相位φ'_a,b，采用多项式曲面拟合法对干涉相位图的解缠相位进行拟合，得到干涉相位图的n阶多项式拟合函数f(x,y)的系数列向量A和存在误差的干涉相位图的解缠相位曲面；

步骤4，根据n阶多项式拟合函数f(x,y)的系数列向量A，并采用遗传算法对得到的存在误差的干涉相位图的解缠相位曲面进行优化，得到优化后的干涉相位图的解缠相位曲面，即干涉相位图的最终解缠相位。

2.如权利要求1所述的一种基于枝切法和曲面拟合的雷达相位解缠方法，其特征在于，在步骤2中，所述得到第一质量区域图的解缠相位φ'_a,b，得到第一质量区域图的解缠相位φ'_a,b的具体子步骤为：

2a)根据残差点的计算公式

$\underset{c=1}{\overset{4}{\Σ}}\widehat{\▿}\mathrm{\φ}\left(c\right)=2\mathrm{\π}q$

得到第一质量区域图中的正残差点坐标和负残差点坐标，并根据该正残差点坐标和负残差点坐标分别将正残差点和负残差点的位置在第一质量区域图中标注出来，得到第一质量区域图的正负残差点分布图；其中，表示像素点(a,b)与所述像素点(a,b)的第c个相邻像素点之差，c表示第一质量区域图中像素点(a,b)的第c个像素点，c∈{1,2,3,4}，q表示判决常数。

2b)选取第一质量区域图的正负残差点分布图中的第一个正残差点为中心，并在搜索半径为r×r的范围内进行搜索，若搜到负残差点，则将该正残差点与搜到的负残差点定义为近点对；相反，若没有搜到负残差点，则将该正残差点定义为远点；其中，r×r表示设定的初始搜索半径；

若该正残差点是远点，则以该远点为中心，在第一质量区域图的正负残差点分布图中依次在搜索半径为r×r，(r+2)×(r+2)，……的搜索范围内再次进行搜索，直到搜到负残差点，并与该远点组成远点对；若搜索到第一质量区域图的正负残差点分布图边界处仍未搜到负残差点，则将第一质量区域图的正负残差点分布图边界点作为负残差点，并与该远点组成远点对；其中，r×r表示设定的初始搜索半径；

逐个选取第一质量区域图的正负残差点分布图中的每一个正残差点重复此过程，进而得到第一质量区域图的正负残差点分布图中的P个近点对和Q个远点对；其中，P表示第一质量区域图的正负残差点分布图中的近点对个数，Q表示第一质量区域图的正负残差点分布图中的远点对个数；

2c)选取第一质量区域图的正负残差点分布图中的第一个远点对的正残差点为中心、以设置的枝切线长为半径的范围内进行搜索，得到所述范围内的C₁个远点对和D₁个近点对，并将所述范围内的C₁个远点对和D₁个近点对进行重新组合，得到第一质量区域图的正负残差点分布图中的第一个远点对的正残差点对应范围内重组后的E₁个近点对和重组后的F₁个远点对，进而得到第一质量区域图的正负残差点分布图中的第一个远点对的正残差点对应范围内的G₁个优化后的枝切线。

逐个选取第一质量区域图的正负残差点分布图中每一个远点对的正残差点重复此过程，直到得到第一质量区域图的正负残差点分布图中的Q个远点对的正残差点各自对应范围内重组后的E₁～E_Q个近点对和重组后的F₁～F_Q个远点对，进而得到第一质量区域图的正负残差点分布图中的Q个远点对的正残差点各自对应范围内的G₁～G_Q个优化后的枝切线。

2d)根据第一质量区域图的正负残差点分布图中的Q个远点对的正残差点各自对应范围内的G₁～G_Q个优化后的枝切线，对第一质量区域图中的缠绕相位梯度进行积分，得到第一质量区域图的解缠相位φ'_a,b；其中，(a,b)表示第一质量区域图中像素点的坐标。

3.如权利要求2所述的一种基于枝切法和曲面拟合的雷达相位解缠方法，其特征在于，所述搜索半径为r×r表示设定的初始搜索半径，所述初始搜索半径为5×5。

4.如权利要求1所述的一种基于枝切法和曲面拟合的雷达相位解缠方法，其特征在于，在步骤3中，所述干涉相位图的n阶多项式拟合函数f(x,y)表示为：

f(x,y)＝α₀+α₁x+α₂y+α₃x²+α₄xy+α₅y²+...+α_myⁿ

其中，α₀,α₁,α₂……α_m表示相位干涉图的n阶多项式拟合函数的系数；所述干涉相位图的n阶多项式拟合函数f(x,y)的第m个系数α_m的下标m和所述干涉相位图的n阶多项式拟合函数f(x,y)的阶数n的关系为：m>n。

5.如权利要求1所述的一种基于枝切法和曲面拟合的雷达相位解缠方法，其特征在于，在步骤3中，所述干涉相位图的n阶多项式拟合函数f(x,y)，其最小均方误差E为：

$E=\underset{s_1}{\Σ}\[{\mathrm{\φ}}_{a,b}^{\′}-f(x_{a,b},y_{a,b}){\]}^2$

其中，φ'_a,b表示解缠后的第一质量区域图中像素点(a,b)的解缠相位，相位干涉图的n阶多项式拟合函数f(x,y)在第一质量区域图中像素点(a，b)处的值用f(x_a,b,y_a,b)表示，x_a,b表示第一质量区域图中像素点(a,b)在x轴方向的坐标，y_a,b表示第一质量区域图中像素点(a,b)在y轴方向的坐标，s₁表示第一质量区域图中的所有像素点，∑表示求和符号，该的范围是第一质量区域图。

6.如权利要求1所述的一种基于枝切法和曲面拟合的雷达相位解缠方法，其特征在于，在步骤4中，所述采用遗传算法对得到的存在误差的干涉相位图的解缠相位曲面进行优化，其具体子步骤为：

4a)根据相位干涉图的n阶多项式拟合函数f(x,y)的系数列向量A，得到初始种群P(0)；

4b)设定最大进化代数为T，并根据评价函数ε²的计算公式

${\mathrm{\ϵ}}^2=\left|\right|\▿\mathrm{\Φ}-\▿\mathrm{\Ψ}|{|}^2,$

和初始种群P(0)，得到下一代种群P(t)；其中，t表示迭代次数，t的初始值为1；其中，表示存在误差的干涉相位图的解缠相位曲面的缠绕相位梯度列向量，表示干涉相位图的n阶多项式拟合函数f(x,y)中解缠相位梯度列向量；

4c)将迭代次数t的值加1，重复子步骤4b)，直到迭代次数t的值等于设置的最大进化代数T，进而得到优化后的干涉相位图的n阶多项式拟合函数f(x,y)的系数；

4d)根据优化后的n阶多项式拟合函数f(x,y)的系数，得到优化后的干涉相位图的解缠相位曲面。