基于无网格与频率估计的干涉合成孔径雷达相位解缠方法
【技术领域】
[0001] 本发明属于雷达信号处理技术领域，更进一步是一种基于无网格与频率估计的干 涉合成孔径雷达相位解缠方法。通过对干涉合成孔径雷达干涉相位图的像素点之间受到缠 绕相位应用基于无网格法和局部频率的InSAR相位解缠法，得到目标真实相位，提高解缠 精度，从而使InSAR数据处理精度得到显著改善。
【背景技术】
[0002] 上世纪中叶开始，雷达技术获得迅速发展，许多新颖的雷达技术陆续被提出并使 用。干涉合成孔径雷达（InterferometricSyntheticApertureRadar(以下简称InSAR) 技术便是代表之一，它是射电天文学干涉测量技术与合成孔径雷达技术两者相结合而产生 的全新技术。在继承了合成孔径雷达获得目标二位图像的特点上，又具备了获取目标区域 高程信息的能力。
[0003] InSAR是由CarlWiley在SAR的基础上首次提出的。它的工作原理是通过对同一 目标区域的具有相干性的两幅SAR图像进行干涉处理，得到一幅该区域的干涉相位图像， 分析处理干涉相位图像从而获取平面和高程方向的信息，完成三维测量。
[0004] 相位信息蕴含着距离差的信息，在干涉处理环节中，获取相位信息是整个技术的 关键所在。经过一系列的预处理得到的干涉相位图中的相位信息是受到干扰的，相位信息 中的相位值被限制在了 [-JT，JT]之间，也就是说干涉相位图像中的相位都是存在误差的 缠绕相位。对存在误差的缠绕相位进行解缠并获得真实的相位信息，是解决相位缠绕问题 的必要手段，根据已知的缠绕相位，利用邻近像素点间的缠绕相位梯度获得与实际距离差 成比例的相位的过程就是相位解缠的思想。
[0005] 1951年InSAR技术被发明至今，针对相位解缠的方法不断涌现，目前已经有了几 十种相位解缠算法。按照解缠方式的不同，大体可分为三类，基于路径跟踪的相位解缠方 法、基于最小范数的相位解缠方法和基于网络规划的相位解缠方法。最小二乘法作为基于 最小范数类算法的代表，其思想是以缠绕相位梯度与解缠相位梯度的均方误差最小作为约 束。最小二乘法具备了良好的稳健性，解缠过程中不依赖积分路径，但是缺点是不对残差点 进行处理，最终导致残差点引起的误差扩散到整个图像。
【发明内容】
[0006] 针对现有技术的上述缺点，本发明的目的在于提出一种基于无网格与频率估计的 干涉合成孔径雷达相位解缠方法，能够有效避免误差向整个干涉相位图的扩散，并且更加 准确的反应目标区域实际地形的变化。
[0007] 为实现上述目的，本发明采用如下技术方案予以实现。
[0008] -种基于无网格与频率估计的干涉合成孔径雷达相位解缠方法，包括如下步骤：
[0009] 步骤1，获取干涉合成孔径雷达的干涉相位图；
[0010] 步骤2,将所述干涉相位图划分为多个子干涉相位图，并计算每个子干涉相位图中 相邻像素点之间的缠绕相位梯度值；
[0011] 步骤3,确定所述每个子干涉相位图的形函数；
[0012] 步骤4,根据所述每个子干涉相位图中相邻像素点之间的缠绕相位梯度值和所述 形函数，得到所述每个子干涉相位图的缠绕相位梯度函数；
[0013] 步骤5,根据所述每个子干涉相位图的缠绕相位梯度函数，求解所述每个子干涉相 位图的解缠相位值；
[0014] 步骤6,根据所述每个子干涉相位图的解缠相位值，对所述干涉合成孔径雷达的干 涉相位图进行解缠。
[0015] 本发明的特点和进一步的改进为：
[0016] (1)所述计算每个子干涉相位图中相邻像素点之间的缠绕相位梯度值，具体包 括：
[0017] 子干涉相位图中任意像素点（X，y)在y方向上和x方向上的缠绕相位梯度值为
[0018]
[0019] 其中，^>、.(^')表示子干涉相位图中任意像素点〇[，7)在 7方向上的缠绕相位梯 度值，表示子干涉相位图中任意像素点（X，y)在X方向上的缠绕相位梯度值，所 述子干涉相位图为任意一个子干涉相位图。
[0020] (2)步骤3具体包括如下子步骤：
[0021] (3a)将所述子干涉相位图中任意像素点（x，y)的缠绕相位梯度函数表
示为：
[0022]
[0023] 其中，QT(x，y)为给定的基函数，qi(x，y)是基函数QT(x，y)在空间坐标XT= [x，y] 上给定的单项式，8;是q; (x,y)的待定系数，向量a是由ajl成的列向量；
[0024] (3b)根据所述子干涉相位图中每个像素点的缠绕相位梯度值，计算向量a的值， 有a二其中，是由n个像素点的缠绕相位梯度值 组成的列向量
是一个常数矩阵；
[0025] (3c)确定子干涉相位图的形函数Pe(x,y) =QT(x,y)Wi。
[0026] (3)步骤4具体为：
[0027] 子干涉相位图的解缠相位梯度函数（T(x，y)、缠绕相位梯度函数在x方向上的分 量▽<〇,)〇、缠绕相位梯度函数在y方向上的分量表示为：
[0028] (i>e(x,y)=Pe(x,y) <i>e
[0029]
[0030]
[0031]其中，夢=|辞讀…痛是由子干涉相位图中n个像素点的解缠相位值 ^巾2…t组成的列向量，
k由子干涉相位图中n个像 素点的缠绕相位梯度值在x方向上的分量组成的列向量，
是由子干涉相位图中n个像素点的缠绕相位梯度值在y方向上的 分量V识，…组成的列向量。
[0032] (4)步骤5具体包括如下子步骤：
[0033] (5a)在子干涉相位图中，构造代价函数r((T(x，y))，其中，
[0034]
[0035] 其中，劳/fo,W表示对子干涉相位图中的形函数p6(x，y)进行求导运算，S彥示 子干涉相位图中各个像素所组成的闭合区域；
[0036] (5b)根据局部频率估计的方法得到
[0037]
[0038] 其中，厂(a%)〇表示x方向上的局部频率列向量，//<X)0表示y方向上的局部频 率列向量；
[0039] (5c)求所述子干涉相位图中代价函数L6 (巾6 (X，y))的极小值，确定所述子干涉相 位图的解缠相位。
[0040] 本发明与现有的技术相比具有以下优点：
[0041] (1)由于本发明是基于最小二乘法的改进算法，所以继承了最小二乘法在求解 InSAR解缠相位时的良好稳健性；(2)现有技术利用滤波相位差分法求解缠相位，而此法因 为滤波效果好坏存在很大误差；而本发明采用的局部频率估计法则一步计算邻近像素的高 度差，因而使用局部频率估计法代替滤波相位差分法求在最小二乘准则下的相位解缠值， 不仅能够更加准确而且能够更加快速的反应实际地形的变化；（3)现有技术利用全局最小 二乘法准则求解缠相位时残差点处的错误会向整个干涉相位图的扩散，导致结果严重错 误；而本发明在解缠过程中引入无网格法，划分局部网格，将误差控制在区域内部，减少了 整体误差，提高了解缠精度。
【附图说明】
[0042] 为了更清楚地说明本发明实施例或现有技术中的技术方案，下面将对实施例或现 有技术描述中所需要使用的附图作简单地介绍，显而易见地，下面描述中的附图仅仅是本 发明的一些实施例，对于本领域普通技术人员来讲，在不付出创造性劳动的前提下，还可以 根据这些附图获得其他的附图。
[0043]图1为本发明实施例提供的相位解缠方法的流程示意图。
[0044] 图2为本发明实施例提供的局部网格的示意图。
[0045] 图3(a)为模拟真实相位的示意图。
[0046] 图3(b)为模拟加入噪声后的缠绕相位示意图。
[0047] 图4(a)为采用枝切法对图3(b)的解缠结果示意图。
[0048] 图4(b)为最小二乘法对图3(b)的解缠结果示意图。
[0049] 图4(c)为质量图引导法对图3(b)的解缠结果示意图。
[0050] 图4(d)为本发明方法对图3(b)的解缠结果示意图。
[0051] 图5 (a)为加入标准差0. 6500rad的噪声后采用枝切法对图3 (b)的解缠结果的误 差分布不意图。
[0052] 图5 (b)为加入标准差0. 6500rad的噪声后采用最小二乘法对图3 (b)的解缠结果 的误差分布示意图。
[0053] 图5 (c)为加入标准差0. 6500rad的噪声后采用质量图引导法对图3 (b)的解缠结 果的误差分布示意图。
[0054] 图5 (d)为加入标准差0. 6500rad的噪声后采用本发明方法对图3 (b)的解缠结果 的误差分布示意图。
[0055] 图6 (a)为加入标准差1.lOOOrad的噪声后采用枝切法对图3 (b)的解缠结果的误 差分布不意图。
[0056] 图6 (b)为加入标准差1.lOOOrad的噪声后采用最小二乘法对图3 (b)的解缠结果 的误差分布示意图。
[0057] 图6 (c)为加入标准差1.lOOOrad的噪声后采用质量图引导法对图3 (b)的解缠结 果的误差分布示意图。
[005