本发明涉及相位解缠领域，涉及基于秩信息滤波的相位解缠，尤其是一种基于秩信息滤波的相位解缠方法。
背景技术：
相位解缠是干涉合成孔径雷达测量(insar)、光学干涉测量、合成孔径声呐(insas)以及核磁共振成像等数据处理过程中的重要步骤之一。由于三角函数的周期性，从干涉测量中获取的表征目标参数的干涉相位被限制在其相位主值(-π,π]区间，俗称称为缠绕相位；从缠绕相位中恢复反应目标参数的真实干涉相位，即为所谓的相位解缠。在理想情况下，通过对缠绕相位加上2π的整数倍即可实现相位解缠。然而，干涉测量过程中的相位跳变、相位噪声、雷达阴影等干扰因素均会显著增加相位解缠的难度。
传统相位解缠方法，例如枝切法、最小二乘法、质量引导法等，这些方法在噪声较小且残差点较少的情况下解缠效果较好，一旦噪声较大且残差点分布密集，这些方法解缠结果误差较大。随后，一类基于非线性滤波与状态估计的解缠算法被相继提出，包括扩展卡尔曼滤波相位解缠算法(ekfpu)、无味卡尔曼滤波相位解缠算法(ukfpu)、容积卡尔曼滤波相位解缠算法(ckfpu)、无味信息滤波相位解缠算法(uifpu)等。这些算法在一定程度上弥补了传统相位解缠算法的不足，但高效与精确地解缠噪声干涉图仍然是一个非常具有挑战性的问题。
技术实现要素：
本发明针对上述现有技术的缺陷，提供一种基于秩信息滤波(rif)的相位解缠方法，这种方法能从噪声缠绕干涉图获得较为稳健的结果。
实现本发明目的的技术方案是：
一种基于秩信息滤波的相位解缠方法，包括如下步骤：
1)根据已有的非线性滤波相位解缠模型，把基于修正矩阵束模型(ampm)的局部相位梯度估计技术与秩信息滤波器结合起来，建立基于秩信息滤波的相位解缠模型，把干涉图相位解缠问题转化为秩信息滤波框架下的状态估计问题；
2)根据步骤1)中得到的秩信息滤波相位解缠模型，建立一维秩信息滤波相位解缠算法；
3)根据步骤2)中得到的一维秩信息滤波相位解缠算法，用二维坐标代替一维坐标，建立二维秩信息滤波相位解缠算法。
进一步地，步骤1)中，利用干涉图相邻像元解缠相位之间的关系，干涉图相位展开系统方程可表示为如下：
式中，xk表示干涉图k像元解缠相位，作为待评估的状态变量；μk-1表示干涉图k-1像元真实相位梯度，可利用基于ampm的局部梯度估计技术来获取其估计值
进一步地，步骤2)中，设干涉图k-1像元状态估计值及其估计误差方差分别为
2-1)根据秩采样原理生成秩信息采样点；
2-2)对干涉图k像元进行状态预测，得到干涉图k像元状态变量一步预测值以及一步预测误差方差；
2-3)对干涉图k像元进行状态空间向信息空间转换，得到干涉图k像元状态变量一步预测的信息矩阵yk和信息向量
2-4)对干涉图k像元进行状态更新，得到干涉图k像元解缠相位以及相位估计误差方差。
进一步地，步骤2-1)中，包括根据秩采样原理生成的秩信息采样点：
式中，
进一步地，步骤2-2)中，干涉图k像元状态变量一步预测值：
干涉图k像元状态变量一步预测误差方差：
式中，
进一步地，步骤2-3)中，干涉图k像元状态变量一步预测的信息矩阵yk和信息向量
利用h-∞算子对信息矩阵yk进行优化，则式(5)中的信息矩阵变为：
式中，γs为衰减因子，通常取0.8≤γs≤2，i是与
进一步地，步骤2-4)中，互协方差
zi,k＝h[χi,k](8)
信息状态分布ik和相应的信息矩阵分布ik为：
式中，ηk为干涉图k像元状态变量观测矢量残差，rk为干涉图k像元状态变量观测噪声方差，
更新的信息矩阵yk和信息向量yk分别为：
更新后的状态估计值及其状态估计误差方差分别为：
其中，
进一步地，步骤3)中，用二维坐标(m,n)代替一维k，设干涉图像元(m,n)为待解缠像元，则该像元初始状态预测值
式中，干涉图(a,s)像元为待解缠像元相邻8个像元中的已解缠像元，
干涉图(m,n)像元的秩信息采样点χi,(m,n)，可按如下计算：
状态变量一步预测值
状态变量一步预测误差方差
其中，q[(m,n)|(a,s)]为基于ampm的局部梯度估计技术导致的过程误差方差，利用基于堆排序的快速路径跟踪策略来指导相位解缠路径，引导rif相位解缠程序沿高质量像元到低质量像元的路径完成对干涉图缠绕像元递推估计。
本方法把高效稳健的秩信息滤波器应用于干涉图相位解缠，结合基于ampm的局部相位梯度估计技术与基于堆排序的快速路径跟踪策略，提出一种基于秩信息滤波的相位解缠算法。该方法建立基于秩信息滤波的相位解缠程序，利用基于ampm的局部相位梯度估计技术获取基于秩信息滤波的相位解缠程序的相位梯度信息；利用基于堆排序的快速路径跟踪策略来指导相位解缠路径，保证基于秩信息滤波的相位解缠程序沿高质量像元到低质量像元路径解缠干涉图。模拟数据以及实测数据实验结果表明本算法的有效性，能从噪声缠绕干涉图获得较为稳健的结果。
附图说明
图1a-图1f为模拟干涉图，其中图1a-图1c为三幅干涉图的真实解缠相位，图1d为图1a真实干涉相位的噪声缠绕相位图，图1e为图1b真实干涉相位的噪声缠绕相位图，图1f为图1c真实干涉相位的噪声缠绕相位图；
图2a-图2c为本发明方法解缠图1d的结果，其中图2a表示展开相位，图2b表示相位展开误差，图2c表示相位展开误差直方图；
图3a-图3c为本发明方法解缠图1e的结果，其中图3a表示展开相位，图3b表示相位展开误差，图3c表示相位展开误差直方图；
图4a-图4c为本发明方法解缠图1f的结果，其中图4a表示展开相位，图4b表示相位展开误差，图4c表示相位展开误差直方图；
图5a-图5c为本发明方法解缠实测数据的结果，其中图5a为局部etna火山干涉图，图5b为解缠相位，图5c为解缠相位重缠绕结果。
具体实施方式
下面将结合本发明实施例中的附图，对本发明实施例中的技术方案进行完整地描述，显然，所描述的实施例仅仅是本发明一部分实施例，而不是全部的实施例。基于本发明中的实施例，本领域普通技术人员在没有做出创造性劳动前提下所获得的所有其他实施例，都属于本发明保护的范围。
实施例：
一种基于秩信息滤波的相位解缠方法，包括如下步骤：
1)根据已有的非线性滤波相位解缠模型，把基于ampm的局部相位梯度估计技术与秩信息滤波器结合起来，建立基于秩信息滤波的相位解缠模型，把干涉图相位解缠问题转化为秩信息滤波框架下的状态估计问题；
利用干涉图相邻像元解缠相位之间的关系，干涉图相位展开系统方程可表示为如下：
式中，xk表示干涉图k像元解缠相位，作为待评估的状态变量；μk-1表示干涉图k-1像元真实相位梯度，可利用基于ampm的局部梯度估计技术来获取其估计值
2)根据步骤1)中得到的秩信息滤波相位解缠模型，建立一维秩信息滤波相位解缠算法；
设干涉图k-1像元状态估计值及其估计误差方差分别为
2-1)根据秩采样原理生成的秩信息采样点：
式中，
2-2)对干涉图k像元进行状态预测，得到干涉图k像元状态变量一步预测值以及一步预测误差方差；
干涉图k像元状态变量一步预测值：
干涉图k像元状态变量一步预测误差方差：
式中，
2-3)对干涉图k像元进行状态空间向信息空间转换，得到干涉图k像元状态变量一步预测的信息矩阵yk和信息向量
利用h-∞算子对信息矩阵yk进行优化，则式(5)中的信息矩阵变为：
式中，γs为衰减因子，通常取0.8≤γs≤2，i是与
2-4)对干涉图k像元进行状态更新，得到干涉图k像元解缠相位以及相位估计误差方差；互协方差
zi,k＝h[χi,k](8)
信息状态分布ik和相应的信息矩阵分布ik为：
式中，ηk为干涉图k像元状态变量观测矢量残差，rk为干涉图k像元状态变量观测噪声方差，
更新的信息矩阵yk和信息向量yk分别为：
更新后的状态估计值及其状态估计误差方差分别为：
其中，
3)根据步骤2)中得到的一维秩信息滤波相位解缠算法，用二维坐标代替一维坐标，建立二维秩信息滤波相位解缠算法。
用二维坐标(m,n)代替一维k，设干涉图像元(m,n)为待解缠像元，则该像元初始状态预测值
式中，干涉图(a,s)像元为待解缠像元相邻8个像元中的已解缠像元，
干涉图(m,n)像元的秩信息采样点χi,(m,n)，可按如下计算：
状态变量一步预测值
状态变量一步预测误差方差
其中，q[(m,n)|(a,s)]为基于ampm的局部梯度估计技术导致的过程误差方差，利用基于堆排序的快速路径跟踪策略来指导相位解缠路径，引导rif相位解缠程序沿高质量像元到低质量像元的路径完成对干涉图缠绕像元递推估计。
为了验证各方法的性能，不同算法包括迭代最小二乘法(ils)、质量引导法(qgpu)以及本方法，在同一matlab软件环境(inteli5-8265u@1.60gcpu 8gbram)下，解缠模拟以及实测干涉图，并对各算法解缠结果进行比较分析。
图1为三幅不同模拟干涉图，大小为256×256像素，其中图1a-图1c为真实干涉相位，图1d-图1f分别为其含噪声的缠绕相位图，其信噪比依次为7.44db、2.18db、0.73db，用秩信息滤波相位解缠(rifpu)算法对三幅干涉图进行解缠。
图2a-图2c为本发明方法解缠图1d的结果，其中图2a表示展开相位，图2b表示相位展开误差，图2c表示相位展开误差直方图；
图3a-图3c为本发明方法解缠图1e的结果，其中图3a表示展开相位，图3b表示相位展开误差，图3c表示相位展开误差直方图；
图4a-图4c为本发明方法解缠图1f的结果，其中图4a表示展开相位，图4b表示相位展开误差，图4c表示相位展开误差直方图；
图5a-图5c为本发明方法解缠实测数据的结果，其中图5a为局部etna火山干涉图，图5b为解缠相位，图5c为解缠相位重缠绕结果。
本方法解缠图1d-图1f所示噪声缠绕相位图的结果，见图2a-图4c，可以看出本方法获得解缠相位与真实干涉图一致性较好，其相位解缠误差较小。表一列出迭代最小二乘法(ils)、质量引导法(qgpu)以及本方法解缠不同信噪比干涉图的均方根误差，可以看出本方法的均方根误差远小于qgpu和ils方法。
表一各算法相位解缠误差
实测数据实验，图5a为局部etna火山干涉图，本方法解缠相位如图5b所示，其解缠相位重缠绕结果如图5c所示。可以看出本方法获得的解缠相位是连续一致的，其重缠绕相位图条纹亦与原始干涉图条纹一致，这表明该方法获得了较为有效的解缠结果。
以上公开的本发明的优选实施例，只是帮助阐述本发明，不限制本发明仅为所述的具体实施方式。显然，根据本说明书的内容，可作很多的修改和变化。本说明书选取并具体描述这些实施例，是为了更好地解释本发明的原理和实际应用，从而使所属技术领域技术人员能很好地理解和利用本发明。
1.一种基于秩信息滤波的相位解缠方法，其特征在于，包括如下步骤：
1)根据已有的非线性滤波相位解缠模型，把基于ampm的局部相位梯度估计技术与秩信息滤波器结合起来，建立基于秩信息滤波的相位解缠模型，把干涉图相位解缠问题转化为秩信息滤波框架下的状态估计问题；
2)根据步骤1)中得到的秩信息滤波相位解缠模型，建立一维秩信息滤波相位解缠算法。
3)根据步骤2)中得到的一维秩信息滤波相位解缠算法，用二维坐标代替一维坐标，建立二维秩信息滤波相位解缠算法。
2.根据权利要求1所述的基于秩信息滤波的相位解缠方法，其特征在于，步骤1)中，利用干涉图相邻像元解缠相位之间的关系，干涉图相位展开系统方程可表示为如下：
式中，xk表示干涉图k像元解缠相位，作为待评估的状态变量；μk-1表示干涉图k-1像元真实相位梯度，可利用基于ampm的局部梯度估计技术来获取其估计值
3.根据权利要求1所述的基于秩信息滤波的相位解缠方法，其特征在于，步骤2)中，设干涉图k-1像元状态估计值及其估计误差方差分别为
2-1)根据秩采样原理生成秩信息采样点；
2-2)对干涉图k像元进行状态预测，得到干涉图k像元状态变量一步预测值以及一步预测误差方差；
2-3)对干涉图k像元进行状态空间向信息空间转换，得到干涉图k像元状态变量一步预测的信息矩阵yk和信息向量
2-4)对干涉图k像元进行状态更新，得到干涉图k像元解缠相位以及相位估计误差方差。
4.根据权利要求3所述的基于秩信息滤波的相位解缠方法，其特征在于，步骤2-1)中，包括根据秩采样原理生成的秩信息采样点：
式中，
5.根据权利要求3所述的基于秩信息滤波的相位解缠方法，其特征在于，步骤2-2)中，干涉图k像元状态变量一步预测值：
干涉图k像元状态变量一步预测误差方差：
式中，
6.根据权利要求3所述的基于秩信息滤波的相位解缠方法，其特征在于，步骤2-3)中，干涉图k像元状态变量一步预测的信息矩阵yk和信息向量
利用h-∞算子对信息矩阵yk进行优化，则式(5)中的信息矩阵变为：
式中，γs为衰减因子，通常取0.8≤γs≤2，i是与
7.根据权利要求3所述的基于秩信息滤波的相位解缠方法，其特征在于，步骤2-4)中，互协方差
zi,k＝h[χi,k](8)
信息状态分布ik和相应的信息矩阵分布ik为：
式中，ηk为干涉图k像元状态变量观测矢量残差，rk为干涉图k像元状态变量观测噪声方差，
更新的信息矩阵yk和信息向量yk分别为：
更新后的状态估计值及其状态估计误差方差分别为：
其中，
8.根据权利要求1所述的基于秩信息滤波的相位解缠方法，其特征在于，步骤3)中，用二维坐标(m,n)代替一维k，设干涉图像元(m,n)为待解缠像元，则该像元初始状态预测值
式中，干涉图(a,s)像元为待解缠像元相邻8个像元中的已解缠像元，
干涉图(m,n)像元的秩信息采样点χi,(m,n)，可按如下计算：
状态变量一步预测值
状态变量一步预测误差方差
其中，q[(m,n)|(a,s)]为基于ampm的局部梯度估计技术导致的过程误差方差，利用基于堆排序的快速路径跟踪策略来指导相位解缠路径，引导rif相位解缠程序沿高质量像元到低质量像元的路径完成对干涉图缠绕像元递推估计。
技术总结