CN109472834A - 一种基于小波变换的卡尔曼滤波相位展开方法 - Google Patents

Abstract

本发明公开了一种基于小波变换的卡尔曼滤波相位展开方法，涉及相位展开技术领域，解决的技术问题是解决卡尔曼滤波相位展开方法中运算成本高以及对计算机内存要求较高问题，该方法包括以下步骤：(1)对干涉图进行N次分解；(2)根据高频相位梯度信息创建掩膜矩阵A；(3)展开最高阶低频干涉相位图s^N；(4)判断相位展开是否完成，检查当前阶是否为最低阶s⁰，是，结束；否则，转入下一步，继续展开较低阶缠绕相位图；(5)展开较低阶缠绕相位图；(6)判断相位展开是否完成，检查当前阶是否为最低阶s⁰，是，结束；否则，N＝N‑1,转第(5)步。本发明在相位展开中运算成本低，在维持高效率的情况下，依然保证了精度，稳健性较强。

Description

一种基于小波变换的卡尔曼滤波相位展开方法

技术领域

本发明涉及相位展开技术领域，尤其涉及一种基于小波变换的卡尔曼滤波相位展开方法。

背景技术

相位展开是指把干涉纹图中的缠绕相位加上2π的整数倍还原出真实相位的过程。传统经典相位展开方法总体分为三类，第一类是基于路径跟踪策略，包含了Goldstein枝切法、质量图引导法和最小不连续相位展开方法；第二类是最小范数方法，代表是最小二乘相位展开法；第三类是基于网络流的各种方法。还有一类方法可以展开相位的同时进行滤波，如卡尔曼滤波类相位展开方法。非线性卡尔曼滤波能在状态估计的同时实现噪声滤波，具有较强噪声抑制能力，将传统质量引导策略与无味卡尔曼滤波方法(UKF)相结合，可以利用相位质量图来引导相位展开路径，使无味卡尔曼滤波方法沿高质量干涉图向低质量干涉图进行相位展开，既可以准确的恢复真实相位，又可以避免直接穿过噪声较大的像元区域导致精度下降。但是卡尔曼滤波方法比传统相位展开方法拥有较高精度的同时，也需要比其他大部分方法更多的运行时间，随着干涉图规模的增大，尤其是针对大规模低信噪比干涉相位图进行相位展开时，质量引导策略的无味卡尔曼滤波相位展开方法其计算成本将急剧增大，故难以直接应用到大规模干涉图相位展开中。

发明内容

针对现有技术的不足，本发明所解决的技术问题是解决卡尔曼滤波相位展开方法中运算成本高以及对计算机内存要求较高问题。

为解决上述技术问题，本发明采用的技术方案是一种基于小波变换的卡尔曼滤波相位展开方法，包括以下步骤：

(1)对干涉图进行N次分解，具体过程如下：

根据干涉图的规模与噪声选取合适的小波分解阶数，然后依次对干涉图的行方向和列方向进行小波分解，高频细节部分为干涉相位梯度信息，之后的分解只需对低频缠绕相位图处理，分解公式如下：

dⁱ⁺¹(n)＝W(sⁱ(2n+1)-sⁱ(2n))

其中，sⁱ⁺¹(n)为第i+1层分解后低频信号，dⁱ⁺¹(n)为第i+1层分解后高频细节，W(·)为缠绕算子，保证每一层分解结果在相位主值区间。

(2)根据高频相位梯度信息创建掩膜矩阵A，具体过程如下：

干涉相位图经过小波分解后得到每一阶的低频缠绕相位与高频梯度信息，低频缠绕相位保存为需要进行相位展开的待展开相位，而由于相位梯度需满足|Δⁱ⁺¹(n)|≤π的条件，所以对于噪声较大的缠绕相位需要在更低阶进行相位展开，计算公式如下：

dⁱ⁺¹(n)＝W(sⁱ(2n+1)-sⁱ(2n))

由公式可以看出：低阶缠绕相位梯度不大于π/2时，高一阶低频缠绕相位梯度可满足约束条件|Δⁱ⁺¹(n+1)|≤π，其对应缠绕相位可在高阶展开，否则需在较低阶进行相位展开，即掩膜矩阵A由高频部分梯度信息是否大于π/2来生成，缠绕相位展开情况与掩膜矩阵元素对应值如下所示：

0：当前相位已经在更高阶被展开；

1：当前相位应在当前阶被展开；

2：当前相位应在更低阶被展开；

低频缠绕相位是否展开取决于对应高频梯度信息，即取决于掩膜矩阵，掩膜矩阵的正确与否直接影响到最终相位展开的精确度与效率。

(3)展开最高阶低频干涉相位图s^N，具体分步骤如下：

1)获取最高阶低频缠绕干涉相位图s^N的相位质量图，指导相位展开方法展开路径；

2)根据所得干涉相位质量图，选取干涉图中相位质量最高的非边界像元作为起始像元，其展开相位为缠绕相位，并设定其估计误差方差，估计误差方差值的范围为(0，1)，根据经验设定；

3)在干涉图中以起始像元为中心，标记与起始像元邻接的上、下、左和右4个缠绕像元为邻接像元，将所有在掩膜矩阵A_N中对应位置值‘1’的邻接像元嵌入待展开像元组；

4)比较所得待展开像元组中所有像元对应的像元质量值，选取最佳待展开像元标记为像元x，然后利用卡尔曼滤波相位展开方法展开此缠绕像元，具体展开步骤如下；

设所述最佳待展开像元x为干涉缠绕相位图中的(m，n)像元，卡尔曼滤波相位展开方法利用干涉图中邻接待展开像元(m，n)的上、下、左、右、上左、下左、上右和下右8个像元中的已展开像元的信息递推估计待展开像元像元 (m，n)的干涉相位；设Ψ为待展开像元(m，n)的8个邻接像元中已展开像元的集合，设像元(k，l)是干涉图中待展开像元像元(m，n)的8个邻接像元中的已展开像元，即(k，l)∈Ψ，像元(k，l)的状态估计及误差方差分别为和

首先，对干涉图中的像元(m，n)干涉相位值按如下卡尔曼滤波方法进行预测估计：

上式中，(m，n)像元为干涉图中的待展开像元，χ_j，(k，l)是(k，l)像元状态估计的Sigma point，为干涉图中(m，n)像元与(k，l)像元之间的相位梯度估计值，SNR_(k，l)是干涉图中(k，l)像元信噪比，为干涉图中(m，n)像元 Sigma point预测值，Q_{(m，n)|(k，l)}为干涉图中的(m，n)像元与(k，l)像元之间的相位梯度估计误差方差，和为相应调节权值系数；

其次，对干涉图中的(m，n)像元状态进行估计：

其中，y_(m，n)和v_(m，n)分别为干涉图中(m，n)像元含噪声的观测矢量及其附加噪声矢量，u_(m，n)和v_(m，n)分别为附加在复干涉信号实部和虚部的噪声，h[x_(m，n)] 为干涉图中(m，n)像元真实观测矢量，表示干涉图中像元(m，n)观测矢量预测值，Π_(m，n)表示干涉图中(m，n)像元增益矩阵，R_(m，n)表示干涉图中的(m，n) 像元观测误差方差，和分别为干涉图中的(m，n)像元状态估计及误差方差；

5)从分步骤3)步获得的待展开像元组中删除像元x，标记s^N中相邻x 的非边界缠绕像元为邻接像元，将所有在掩膜矩阵中对应位置值为‘1’的邻接像元嵌入待展开像元组；

6)检查待展开像元组是否为空，若没有待展开像元，干涉图s^N相位展开完成；否则，转入分步骤4)步。

(4)判断相位展开是否完成

检查当前阶是否为最低阶s⁰，是，结束；否则.转入下一步，继续展开较低阶缠绕相位图；

(5)展开较低阶缠绕相位图，具体分步骤如下：

1)对s^N中每一个展开像元和d^N中对应的相位梯度进行逆小波变换，得到s^N-1中两个对应像元的展开相位，实现对干涉图的s^N-1更新；更新之后的干涉相位图s^N-1含有已在较高阶展开的像元、应在当前阶展开的像元、应在较低阶展开的像元(最低阶s⁰只含有前两者)；

2)标记中s^N-1所有展开像元的非边界相邻缠绕像元为邻接像元组，将所有在A_N-1中对应位置为‘1’的邻接像元嵌入待展开像元组；

3)对待展开像元组中质量最高的像元x利用卡尔曼滤波相位展开方法展开；

4)从第2)步获得的待展开像元组中删除像元x，标记s^N-1中相邻x的非边界缠绕像元为邻接像元，将所有在掩膜矩阵A_N-1中对应位置值为‘1’的邻接像元嵌入待展开像元组；

5)检查待展开像元组是否为空，若没有待展开像元，相位干涉图s^N-1展开完成；否则，转入第3)步。

(6)判断相位展开是否完成，检查当前阶是否为最低阶s⁰，是，结束；否则，N＝N-1，转第(5)步。

与现技术相比本发明有益效果：

(1)通过对较大的原始相位干涉图进行小波分解，把对大尺寸原始干涉图的展开任务主要转移到了对干涉图最高阶的低频缠绕相位图的展开中，利用质量图引导策略的卡尔曼滤波方法对最高阶缠绕相位展开，既大大提高了效率，也维持了较高精度；

(2)与基于质量引导策略的卡尔曼滤波相位展开方法比较，数据处理展开方法略低，运行时间远远小于基于质量引导策略的卡尔曼相位展开方法，表明本发明在维持高效率的情况下，依然保证了精度，稳健性较强。

附图说明

图1为本发明流程示意图；

图2为实施例1所针对的多山地形场景模拟相位图；

图3为对图2模拟相位图进行缠绕得到的缠绕相位图；

图4为实施例1对最高阶低频缠绕相位展开结果；

图5为实施例1相位展开结果；

图6为实施例1相位展开误差直方图；

图7为基于质量引导策略的卡尔曼滤波相位展开方法相位展开结果；

图8为基于质量引导策略的卡尔曼滤波相位展开方法相位展开误差直方图；

图9为实施例2针对的部分Etna火山实测干涉图；

图10为实施例2对部分Etna火山实测数据相位展开结果。

具体实施方式

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

实施例1：

采用多山地形场景干涉图进行相位展开，图2为多山地形场景的模拟真实相位图，其横坐标为距离向像素，纵坐标为方位向像素，规模为256×256，灰度表示干涉相位值，单位为弧度；图3为多山地形场景的含噪声缠绕相位图，横、纵坐标与图2相同，其信噪比为6.1dB。

图1示出了一种基于小波变换的卡尔曼滤波相位展开方法，包括以下步骤：

(1)对干涉图进行分解，具体过程如下：

根据噪声与干涉图规模选取小波变换分解阶数为1阶。对图3二维干涉图进行模小波变换分解。先对干涉图行方向进行小波分解，高频信息为相位梯度，不需进行处理，然后对低频缠绕相位进行列方向小波分解，即完成对初始干涉图的1阶小波分解。

(2)根据高频相位梯度信息创建掩膜矩阵A：

经小波分解后的高频系数，也就是梯度，根据梯度是否大于π/2生成掩膜矩阵。对于地形较为陡峭或者噪声较大的区域，将被放在较低阶干涉相位图进行相位展开。

(3)展开最高阶低频干涉相位图s¹，具体分步骤如下：

1)获取最高阶干涉相位质量图，反映二维最高阶低频干涉相位图s¹像元质量，引导相位展开路径。本实施例采用的干涉相位质量图为相位微分偏差图；

2)根据所得干涉相位质量图，选取干涉图中相位质量最高的非边界像元作为起始像元，其展开相位为缠绕相位，并设定其估计误差方差，本例中设为0.6；

3)在干涉图s¹中以起始像元为中心，标记与起始像元邻接的上、下、左和右4个缠绕像元为邻接像元，将所有在掩膜矩阵A_N中对应位置值为‘1’的邻接像元嵌入待展开像元组；

4)比较所得待展开像元对应相位质量值，选取最佳待展开像元标记为像元x，然后利用卡尔曼滤波相位展开方法(UKF算法)展开此缠绕像元；

5)从分步骤3)步获得的待展开像元组中删除像元x，标记s¹中相邻x 的非边界缠绕像元为邻接像元，将所有在掩膜矩阵中对应位置值为‘1’的邻接像元嵌入待展开像元组；

6)检查待展开像元组是否为空，若没有待展开像元，干涉图s¹相位展开完成；否则，转入分步骤4)步。干涉相位图s¹展开图如图4所示。

(4)判断相位展开是否完成

当前展开的为低频干涉相位图s¹，，转入下一步，继续展开较低阶缠绕相位图s⁰；

(5)展开较低阶缠绕相位图，具体分步骤如下：

1)对s¹中每一个展开像元和d¹中对应的相位梯度进行行方向逆小波变换，得的s⁰中两个对应像元的展开相位；

2)标记中s⁰所有展开像元的非边界相邻缠绕像元为邻接像元组，将所有在A₀中对应位置为‘1’的邻接像元嵌入待展开像元组；

3)对待展开像元组中质量最高的像元x利用卡尔曼滤波相位展开方法展开；

4)从第2)步获得的待展开像元组中删除像元x，标记s⁰中相邻x的非边界缠绕像元为邻接像元，将所有在掩膜矩阵A₀中对应位置值为‘1’的邻接像元嵌入待展开像元组；

5)检查待展开像元组是否为空，若没有待展开像元，低频缠绕干涉相位图展开完成；否则，转入第3)步。

(6)判断相位展开是否完成，检查当前阶是否为最低阶s⁰，是，结束；否则，N＝N-1，转第(5)步。

图5为本实施例基于小波变换的卡尔曼滤波相位展开方法对图3展开的结果示意图，示意图横纵坐标与图2相同；图6为本实施例相位展开结果的误差直方图，其横坐标为相位展开误差，纵坐标为误差统计频数。

图7为质量引导策略的卡尔曼滤波相位展开方法对图3展开的结果，示意图横纵坐标与图5相同；图8为质量引导策略的卡尔曼滤波相位展开结果的误差直方图，示意图横纵坐标与图6相同。

由图5～图8可以看出本发明实施例方法对干涉相位展开效果良好，具有一定的抗噪抑制能力，且精度较高；表1列出了本发明与质量引导策略的无味卡尔曼滤波相位展开方法在同一MATALB计算环境下展开较大规模干涉相位图所消耗时间，可以看出本发明展开大规模干涉相位图的运行时间远远小于质量引导策略的无味卡尔曼滤波相位展开方法。

表1相位展开方法运行时间

实施例2：

针对Etna火山干涉图进行相位展开，图9为预滤波后的部分Etna火山干涉图。

采用与实施例1所述的步骤(1)至步骤(5)步展开该缠绕干涉图，本例相位展开结果见图10，可以看出本发明方法展开后相位连续且平滑，反映本发明相位展开方法的有效性。将图10重缠绕后相位图条纹与原始缠绕相位图条纹一致，说明本发明方法在实测数据处理中获得有效结果。

本发明是质量引导策略卡尔曼滤波方法与小波变换的相结合的结果，在实测与模拟数据处理实验中其有效性明显。

与现技术相比本发明有益效果：

(1)通过对较大的原始相位干涉图进行小波分解，把对大尺寸原始干涉图的展开任务主要转移到了对干涉图最高阶的低频缠绕相位图的展开中，利用质量图引导策略的卡尔曼滤波方法对最高阶缠绕相位展开，既大大提高了效率，也维持了较高精度；

(2)与基于质量引导策略的卡尔曼滤波相位展开方法比较，数据处理展开方法略低，运行时间远远小于基于质量引导策略的卡尔曼相位展开方法，表明本发明在维持高效率的情况下，依然保证了精度，稳健性较强。

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

Claims (5)

1.一种基于小波变换的卡尔曼滤波相位展开方法，其特征在于，包括以下步骤：

(1)对干涉图进行N次分解；

(2)根据高频相位梯度信息创建掩膜矩阵A；

(3)展开最高阶低频干涉相位图s^N；

(4)判断相位展开是否完成，检查当前阶是否为最低阶s⁰，是，结束；否则，转入下一步，继续展开较低阶缠绕相位图；

(5)展开较低阶缠绕相位图；

(6)判断相位展开是否完成，检查当前阶是否为最低阶s⁰，是，结束；否则，N＝N-1,转第(5)步。

2.如权利要求1所述的基于小波变换的卡尔曼滤波相位展开方法，其特征在于，步骤(1)具体过程如下：

根据干涉图的规模与噪声选取合适的小波分解阶数，然后依次对干涉图的行方向和列方向进行小波分解，高频细节部分为干涉相位梯度信息，之后的分解只需对低频缠绕相位图处理，分解公式如下：

dⁱ+¹(n)＝W(sⁱ(2n+1)-sⁱ(2n))

其中，sⁱ⁺¹(n)为第i+1层分解后低频信号，dⁱ⁺¹(n)为第i+1层分解后高频细节，W(·)为缠绕算子，保证每一层分解结果在相位主值区间。

3.如权利要求1所述的基于小波变换的卡尔曼滤波相位展开方法，其特征在于，步骤(2)具体过程如下：

干涉相位图经过小波分解后得到每一阶的低频缠绕相位与高频梯度信息，低频缠绕相位保存为需要进行相位展开的待展开相位，而由于相位梯度需满足|Δⁱ⁺¹(n)|≤π的条件，所以对于噪声较大的缠绕相位需要在更低阶进行相位展开，计算公式如下：

dⁱ⁺¹(n)＝W(sⁱ(2n+1)-sⁱ(2n))

由公式可以看出:低阶缠绕相位梯度不大于π/2时，高一阶低频缠绕相位梯度可满足约束条件|Δⁱ⁺¹(n+1)|≤π，其对应缠绕相位可在高阶展开，否则需在较低阶进行相位展开，即掩膜矩阵A由高频部分梯度信息是否大于π/2来生成，缠绕相位展开情况与掩膜矩阵元素对应值如下所示：

0：当前相位已经在更高阶被展开；

1：当前相位应在当前阶被展开；

2：当前相位应在更低阶被展开；

低频缠绕相位是否展开取决于对应高频梯度信息，即取决于掩膜矩阵，掩膜矩阵的正确与否直接影响到最终相位展开的精确度与效率。

4.如权利要求1所述的基于小波变换的卡尔曼滤波相位展开方法，其特征在于，步骤(3)具体分步骤如下：

1)获取最高阶低频缠绕干涉相位图s^N的相位质量图，指导相位展开方法展开路径；

2)根据所得干涉相位质量图，选取干涉图中相位质量最高的非边界像元作为起始像元，其展开相位为缠绕相位，并设定其估计误差方差，估计误差方差值的范围为(0，1)，根据经验设定；

3)在干涉图中以起始像元为中心，标记与起始像元邻接的上、下、左和右4个缠绕像元为邻接像元，将所有在掩膜矩阵A_N中对应位置值‘1’的邻接像元嵌入待展开像元组；

4)比较所得待展开像元组中所有像元对应的像元质量值，选取最佳待展开像元标记为像元x，然后利用卡尔曼滤波相位展开方法展开此缠绕像元，具体展开步骤如下；

设所述最佳待展开像元x为干涉缠绕相位图中的(m,n)像元，卡尔曼滤波相位展开方法利用干涉图中邻接待展开像元(m,n)的上、下、左、右、上左、下左、上右和下右8个像元中的已展开像元的信息递推估计待展开像元像元(m,n)的干涉相位；设Ψ为待展开像元(m,n)的8个邻接像元中已展开像元的集合，设像元(k,l)是干涉图中待展开像元像元(m,n)的8个邻接像元中的已展开像元，即(k,l)∈Ψ，像元(k,l)的状态估计及误差方差分别为和

首先，对干涉图中的像元(m,n)干涉相位值按如下卡尔曼滤波方法进行预测估计：

上式中，(m,n)像元为干涉图中的待展开像元，χ_j,(k,l)是(k,l)像元状态估计的Sigmapoint，为干涉图中(m,n)像元与(k,l)像元之间的相位梯度估计值，SNR_(k,l)是干涉图中(k,l)像元信噪比，为干涉图中(m,n)像元Sigma point预测值，Q_(m,n)|(k,l)为干涉图中的(m,n)像元与(k,l)像元之间的相位梯度估计误差方差，和为相应调节权值系数；

其次，对干涉图中的(m,n)像元状态进行估计：

其中，y_(m,n)和ν_(m,n)分别为干涉图中(m,n)像元含噪声的观测矢量及其附加噪声矢量，u_(m,n)和v_(m,n)分别为附加在复干涉信号实部和虚部的噪声，h[x_(m,n)]为干涉图中(m,n)像元真实观测矢量，表示干涉图中像元(m,n)观测矢量预测值，Π_(m,n)表示干涉图中(m,n)像元增益矩阵，R_(m,n)表示干涉图中的(m,n)像元观测误差方差，和分别为干涉图中的(m,n)像元状态估计及误差方差；

5)从分步骤3)步获得的待展开像元组中删除像元x，标记s^N中相邻x的非边界缠绕像元为邻接像元，将所有在掩膜矩阵中对应位置值为‘1’的邻接像元嵌入待展开像元组；

6)检查待展开像元组是否为空，若没有待展开像元，干涉图s^N相位展开完成；否则，转入分步骤4)步。

5.如权利要求1所述的基于小波变换的卡尔曼滤波相位展开方法，其特征在于，步骤(5)具体分步骤如下：

1)对s^N中每一个展开像元和d^N中对应的相位梯度进行逆小波变换，得到s^N-1中两个对应像元的展开相位，实现对干涉图的s^N-1更新；更新之后的干涉相位图s^N-1含有已在较高阶展开的像元、应在当前阶展开的像元、应在较低阶展开的像元(最低阶s⁰只含有前两者)；

2)标记中s^N-1所有展开像元的非边界相邻缠绕像元为邻接像元组，将所有在A_N-1中对应位置为‘1’的邻接像元嵌入待展开像元组；

3)对待展开像元组中质量最高的像元x利用卡尔曼滤波相位展开方法展开；

4)从第2)步获得的待展开像元组中删除像元x，标记s^N-1中相邻x的非边界缠绕像元为邻接像元，将所有在掩膜矩阵A_N-1中对应位置值为‘1’的邻接像元嵌入待展开像元组；

5)检查待展开像元组是否为空，若没有待展开像元，相位干涉图s^N-1展开完成；否则，转入第3)步。