一种无损卡尔曼滤波与粒子滤波相结合的相位展开算法的制作方法

本发明涉及干涉相位图展开技术领域，尤其涉及一种无损卡尔曼滤波(ukf)与粒子滤波(pf)相结合的相位展开算法。

背景技术：

干涉合成孔径雷达(interferometricsyntheticapertureradar,insar)是以sar复图像对的干涉相位为信息源，通过一系列数据处理环节，进而反演得到地表的三维信息和变化信息。insar是一个交叉性很强的领域，涉及微波遥感，电磁波传播，信号处理，图像处理等学科。insar最初的发展目的是为了进行地形测绘，应用最为成功的就是获得高精度dem数据。随着技术的发展，insar应用领域在不断扩展，主要包括：地表形变探测，动目标检测，海洋测绘，森林制图，洪涝监测，交通监测和冰川研究等。干涉相位展开是insar技术应用中尤为关键的环节，其相位展开精度直接影响着insar系统高程测量的精度，故干涉图展开问题受到越来越多的关注。

传统相位展开算法，包括枝切法、质量图引导法、网络流法、最小二乘法等，要在相位展开之前使用前置预滤波器抑制干涉图中的相位噪声，以提高相位展开精度，但在尽可能地除去相位噪声的同时却易模糊掉干涉条纹的边缘特性。

随后，一类基于数据融合的相位展开算法陆续被提出，包括扩展卡尔曼滤波相位展开算法(ekfpu)、不敏卡尔曼滤波相位展开算法(ukfpu)、粒子滤波相位展开算法等。这类算法利用非线性滤波器自身具有的噪声抑制能力，可以在展开缠绕像元的同时抑制相位噪声，在一定程度上弥补了传统相位展开算法的不足,但还是存在着误差大，精度低，稳定性差的缺点。

技术实现要素：

针对现有技术的不足，本发明解决的技术问题是如何解决相位展开技术中容易出现误差大，精度低，稳定性差的问题。

为解决上述技术问题，本发明采用的技术方案是一种无损卡尔曼滤波与粒子滤波相结合的相位展开算法，把无损卡尔曼滤波算法与粒子滤波算法相结合应用于缠绕相位图像的展开中，利用相位质量图沿干涉图高质量到低质量的路径展开缠绕相位，并结合基于修正矩阵束模型的局部相位梯度估计算法,包括以下步骤：

(1)获取解缠缠绕相位图的有关数据，其具体步骤如下：

1)获得干涉相位质量图，对干涉图进行处理，获得反映干涉像元质量的二维数据图，即干涉相位质量图，采用的是相位微分偏差图，其定义如下：

其中分别是干涉图距离向和方位向的相位偏导数，是在k×k局部窗口中的相位偏导数的平均值；

2)利用ampm局部相位梯度估计算法，获取干涉图局部相位梯度信息；

3)利用干涉相位质量图引导算法，沿干涉相位图从高质量像元到低质量像元的路径展开缠绕相位实现相位解缠；

(2)对每个缠绕干涉图像元用无损卡尔曼滤波(ukf)算法进行初步解缠，其具体步骤如下：

1)利用干涉图相邻像元干涉相位之间的关系，把归一化的复干涉的同向分量和正交分量分别作为干涉相位的两个观测值，沿某一路径可得无损卡尔曼滤波(ukf)系统方程：

上式中，x(m,n)表示干涉图(m,n)像元真实干涉相位，和w(m,n)|(k,l)分别为干涉图像(m,n)元和(k,l)像元之间的相位梯度估计值及估计误差，y(m,n)和v(m,n)分别为干涉图(m,n)像元的观测值及其附加噪声；

2)针对上述系统方程，使用无损卡尔曼滤波(ukf)算法预测估计，按如下进行：

上式中，干涉图(m,n)像元为待展开像元，(k,l)像元是干涉图(m,n)像元八个邻接像元中的已展开像元，其状态估计及误差方差分别为和χj,(k,l)是(k,l)像元状态估计的sigmapoint，ψ为干涉图(m,n)像元的八个邻接像元的集合，是干涉图(k,l)像元信噪比，为干涉图(m,n)像元的sigmapoint预测值，q(m,n)|(k,l)为干涉图(m,n)像元和(k,l)像元之间的相位梯度估计误差方差，和为相应调节权值系数；

3)对干涉图(m,n)像元进行状态估计，按如下进行：

其中，y(m,n)和分别表示干涉图(m,n)像元观测值和预测值，w(m,n)表示干涉图(m,n)像元增益矩阵，r(m,n)表示干涉图(m,n)像元观测误差方差，和为干涉图(m,n)像元初步状态估计及误差方差；

4)对经过无损卡尔曼滤波(ukf)算法后的有关参数进行修正；

修正经过无损卡尔曼滤波(ukf)算法后的(m,n)像元误差方差相位梯度估计误差方差qk以及量测误差方差rk，其表达式如下：

qk＝b*qk

rk＝c*rk

其中，a，b，c分别是(m,n)像元误差方差相位梯度估计误差方差qk以及量测误差方差rk的调整参数；

(3)对经过初步解缠过的干涉相位图像元，用粒子滤波(pf)算法进一步解缠，其具体步骤如下：

1)抽取粒子：从均值和方差分别为和的正态分布中抽取粒子

2)对各粒子权值进行计算并归一化处理：

其中，r为量测噪声方差；δz1为相位观测值实部与对应的粒子预测值复相位值实部的差，δz2为相位观测值虚部与对应的粒子预测值复相位值虚部的差；q1(n)和q2(n)分别为第n个粒子所对应的实部与虚部的概率值；为第n个粒子的概率值，它是该粒子实部概率值的平方与虚部概率值的平方和的开根号的值；q(n)为n个粒子概率加权归一化所得到的第n个粒子的权值；

3)进行相位状态估计：

即为所求展开相位。

采用本发明的技术方案具有较高的精度和稳健性，利用无损卡尔曼滤波与粒子滤波本身具有一定的噪声抑制能力，可以降低前置预滤波器的难度与复杂度，甚至可以在处理含噪声干涉图时免去前置预滤波处理步骤，精确和高效地展开缠绕像元。

附图说明

图1为本发明的流程图；

图2为模拟复杂多山地形干涉图，其中：

图2(a)模拟多山地形三维效果图；

图2(b)为模拟多山地形场景真实干涉相位图；

图2(c)为多山地形场景的含有噪声缠绕相位图；

图2(d)为相位残差点图；

图3为模拟数据实验中本发明算法展开结果图，其中：

图3(a)为三峡实测地形干涉相位图；

图3(b)为图3(a)均值滤波器(3×3)滤波后所得干涉图；

图3(c)为相位残差点图；

图4为三峡实测地形干涉图，其中：

图4(a)为用本发明方法对图3(b)的相位展开结果；

图4(b)为重缠绕相位图；

图4(c)为相位残差点图；

图5为实测数据实验中本发明算法展开结果图，其中：

图5(a)为用本发明方法对图4(b)的相位展开结果图；

图5(b)为重缠绕相位图。

具体实施方式

下面结合附图和实施例对本发明作进一步的说明，但不是对本发明的限定。

图1示出了本发明流程，为无损卡尔曼滤波(ukf)与粒子滤波(pf)相结合的相位展开过程。

实施例1

在无预滤波情况下，在模拟数据实验中，针对模拟复杂多山地形场景干涉相位图进行相位展开。

无损卡尔曼滤波(ukf)与粒子滤波(pf)相结合的相位展开方法具体步骤如下：

本发明利用无损卡尔曼滤波(ukf)与粒子滤波(pf)相结合的相位展开方法对干涉相位图进行相位解缠。具体为对缠绕干涉相位图分两步进行解缠，即对于每个缠绕干涉图像元首先用无损卡尔曼滤波(ukf)算法进行初步解缠，再用粒子滤波(pf)算法对经过初步解缠过的干涉相位图像元进一步解缠，从而达到较好的解缠效果。选择起始解缠像元时，在干涉图质量较好的非边界区域选取一像元作为起始像元，并默认它为已展开像元，即把其缠绕相位作为展开相位。

获取解缠缠绕相位图的有关数据，具体步骤如下：

获得干涉相位质量图，对干涉图进行处理获得反映干涉像元质量的二维数据图，即干涉相位质量图，本文中采用的是相位微分偏差图。

其中分别是干涉图距离向和方位向的相位偏导数，是在k×k局部窗口中的相位偏导数的平均值。

本发明利用ampm局部相位梯度估计算法，获取干涉图局部相位梯度信息。

利用干涉相位质量图，引导本发明相位展开方法沿干涉相位图从高质量像元到低质量像元的路径展开缠绕相位，效率及精度高，稳健性强，能够快速实现相位解缠。

用无损卡尔曼滤波(ukf)算法对缠绕干涉相位图像元进行初步解缠的步骤如下：

利用干涉图相邻像元干涉相位之间的关系，把归一化的复干涉的同向分量和正交分量分别作为干涉相位的两个观测值，沿某一路径可得无损卡尔曼滤波(ukf)系统方程：

上式中，x(m,n)表示干涉图(m,n)像元真实干涉相位，和w(m,n)|(k,l)分别为干涉图像(m,n)元和(k,l)像元之间的相位梯度估计值及估计误差，y(m,n)和v(m,n)分别为干涉图(m,n)像元的观测值及其附加噪声。

针对上述系统方程，无损卡尔曼滤波(ukf)算法预测估计按如下进行：

上式中，干涉图(m,n)像元为待展开像元，(k,l)像元是干涉图(m,n)像元八个邻接像元中的已展开像元，其状态估计及误差方差分别为和χj,(k,l)是(k,l)像元状态估计的sigmapoint，ψ为干涉图(m,n)像元的八个邻接像元的集合，是干涉图(k,l)像元信噪比，为干涉图(m,n)像元的sigmapoint预测值，q(m,n)|(k,l)为干涉图(m,n)像元和(k,l)像元之间的相位梯度估计误差方和为相应调节权值系数。

干涉图(m,n)像元状态估计按如下进行：

其中，y(m,n)和分别表示干涉图(m,n)像元观测值和预测值，w(m,n)表示干涉图(m,n)像元增益矩阵，r(m,n)表示干涉图(m,n)像元观测误差方差，和为干涉图(m,n)像元初步状态估计及误差方差。

修正经过无损卡尔曼滤波(ukf)算法后的有关参数。

修正经过卡尔曼滤波后的(m,n)像元误差方差相位梯度估计误差方差qk以及k像元量测误差方差rk，其表达式如下：

qk＝b*qk

rk＝c*rk

其中a，b，c分别是(m,n)像元误差方差相位梯度估计误差方差qk以及量测误差方差rk的调整参数。

用粒子滤波(pf)算法对经过无损卡尔曼滤波(ukf)算法初步解缠过的缠绕干涉图像元进一步解缠的步骤如下：

抽取粒子：从均值和方差分别为和的正态分布中抽取粒子

计算各粒子权值并归一化：

其中，r为量测噪声方差；δz1为相位观测值实部与对应的粒子预测值复相位值实部的差，δz2为相位观测值虚部与对应的粒子预测值复相位值虚部的差；q1(n)和q2(n)分别为第n个粒子所对应的实部与虚部的概率值；为第n个粒子的概率值，它是该粒子实部概率值的平方与虚部概率值的平方和的开根号的值；q(n)为n个粒子概率加权归一化所得到的第n个粒子的权值。

相位状态估计：

即为所求展开相位。

需注意的是，在本发明实验中，ampm局部相位梯度估计算法所使用的局部窗口大小为7×7，起始像元的估计误差方差设为0.6，本发明所使用的路径引导图均为微分偏差图，粒子滤波算法中粒子数ns＝40。

图2(a)为本发明相位展开方法实施例1所针对的模拟多山地形场景的三维效果图；图2(b)为模拟多山地形场景真实干涉相位图，其横坐标为距离向像素，纵坐标为方位向像素，灰度表示干涉相位值，单位为弧度；图2(c)为多山地形场景的含有噪声缠绕相位图，纵横坐标与图2(b)相同，其信噪比为0.73db；图2(d)为相位残差点图。

图3(a)为用本算法解缠缠绕干涉相位展开图；3(b)干涉相位展开误差统计直方图；图3(c)为干涉相位展开误差图。

实施例2

在实测数据实验中，针对实测数据为三峡干涉图的一部分，验证本发明算法的性能。

图4(a)为三峡实测地形干涉相位图，图4(b)为图4(a)均值滤波器(3×3)滤波后所得干涉图，图4(c)为相位残差点图。

图5(a)为用本发明方法对图4(b)的相位展开结果，图5(b)为重缠绕相位图，重缠绕图条纹与原始缠绕相位图条纹的相似程度可以反映本发明方法的有效性，可以看出本发明方法其重缠绕相位图条纹与原始缠绕相位图4(b)条纹一致，可以看出本发明方法在实测数据处理中得到较好的效果。

本发明算法基于改进的无损卡尔曼滤波(ukf)与粒子滤波(pf)相结合的相位展开方法，在模拟和实测数据处理实验中对干涉相位图缠绕相位的解缠均能得到很好的效果，表明本发明算法的有效性。

以上结合附图和实施例对本发明的实施方式做出了详细说明，但本发明不局限于所描述的实施方式。对于本领域技术人员而言，在不脱离本发明的原理和精神的情况下，对这些实施方式进行各种变化、修改、替换和变型仍落入本发明的保护范围内。