基于质量引导的InSAR快速相位解缠方法

蒋留兵 ，肖志涛，车 俐，赵纪奎

1.桂林电子科技大学 信息与通信学院，广西 桂林 541004

2.桂林电子科技大学 广西无线宽带通信与信号处理重点实验室，广西 桂林 541004

1 引言

干涉合成孔径雷达（InSAR）是测量地表形变和三维信息的一项技术[1-2]，其关键步骤是要得到由两幅或多幅相干复图像衍生出的真实相位信息。但是在实际运算过程中，由于引入了反三角函数，相位图中每个像素点的相位测量值位于(-π,π)之间，是关于2π模糊的，即得到的只是相位主值。因此，要给每个相位测量值加上2π的整数倍，这种由相位主值恢复出真实相位的过程称为相位解缠[3-5]。相位解缠的精度直接影响最终DEM和地形形变测量的精度，因此它是InSAR干涉信号处理中非常关键的一个步骤。

随着国内外InSAR相位解缠技术的快速发展，目前已经出现了诸多算法，总体来说概括为两大类：局部法和全局法。局部法利用相位数据的局部信息，通过采取不同的策略选择最优的积分路径，将相位解缠过程中产生的误差尽可能地控制在噪声区域内，从而阻止相位误差的传递；最典型的有质量引导法、掩膜分割法[6]、Goldstein枝切法[7-8]等。全局法是将相位解缠问题转化为求解最小范数的极值问题，比较经典的有最小二乘法[9]，它可分为加权和不加权[10]的最小二乘法，该方法能得到一个全局最优解，但其结果比较平滑，误差较大；还有基于网络规划的网络流法[11-12]，它在高相干区域能很好地消除相位噪声所带来的影响。

在众多算法中，很难找到一种综合性能很好的相位解缠算法，这也是它值得去一直研究的原因。其中路径跟踪法是研究最久、应用最广泛的。在此特别想强调的是其中的质量引导法[13-14]，它通过借助辅助信息（质量图）来指导解缠路径，从高质量像素点逐次解缠，直到整幅干涉相位图解缠完毕，相对枝切法而言，它具有更高的解缠精度。然而在质量图指导的过程中需要进行大量的排序工作，造成相位解缠的速度很慢，以及质量图的好坏，会影响最终解缠的精度，因此，需要找寻可靠的质量图和高效的排序算法，来提高相位解缠的效率精度和效率。文献[15]提出基于堆排序法，相对传统方法来说，相位解缠的时间复杂度由O(n)降低为O(lbn)，其大大提高了解缠效率，但是对于大型的堆数据而言，它每次插入一个新元素后调整为最大堆或最小堆过程还是相当费时。因此，本文在此基础上进行改进，提出了一种基于索引分段堆排序的质量引导方法，并且结合传统质量图（相位导数方差图、伪相关系数图）的优点，定义了新的质量图（相位导数方差相关图），其不仅大大地减少了解缠的时间，同时解缠精度也有了一定的提高。

2 质量引导法

2.1 解缠路径的跟踪策略

在已知质量图的指导下，按像素点的质量高低逐次完成解缠，其具体操作流程如下：首先选择一个高质量的点作为已知解缠点，并检测它的四个领域像素点，然后解缠其最高质量像素点，并将该像素点的邻域像素点（未解缠）存放至“邻接列”中，并且将“邻接列”所有的像素点按其质量值进行排序，再将“邻接列”中最高质量的像素点移出，对其进行解缠，更新“邻接列”，重复以上步骤直到邻接列为空时，表示整幅相位图中的像素点都已解缠完毕。图1为其路径跟踪示意图。

图1 解缠路径跟踪示意图

2.2 质量图的确定

传统的相位质量图有四种：最大相位梯度图、相位导数方差图、伪相关系数图、相关系数图。在实际应用中，有些情况下可能只能获取一幅InSAR影像，而相关系数图必须要有两幅SAR图像才行，因此在这里，拿伪相关系数图和相位导数方差图来进行探讨。

伪相关系数图的定义如下：

其中ψi,j是干涉图的缠绕相位值；k为取样窗口的大小，可以视实际情况定义不同大小的窗口。该质量图能够直接从干涉相位图中生成，在实际应用中，噪声较少的情况下，它能很好地反映出真实相位的质量信息，但有时把相位梯度变化大的区域误认为低质量区域，而不能很好地评估相位图的质量。

相位导数方差图定义如下：

本文采用相位导数方差相关质量图（QPDVC），它结合了相位导数方差图和伪相关系数图的优点，在一定程度上，既可以消除相位变化大带来的影响，同时又提高了噪声鲁棒性。其定义如下：

3 索引分段堆排序相位解缠方法

堆是一种基于二叉树的数据结构，除根节点外的其他任一节点都含有两个指针，其中一个指向父节点，另一个指针指向子节点。堆排序是一种稳定的排序算法，它有一个属性就是所有的父子对节点数值在数据结构中都是有序的，因此分为两种：最大堆和最小堆。根据确定的质量图，本文将最小堆作为“邻接列”的数据结构，即所有节点的值都大于其子节点的值，如此，堆中根节点的值就是对应“邻接列”中的最小元素值。其中堆排序最基本操作为：删除和插入。图2为插入新元素的示意图；图3为删除根元素的示意图。

图2 最小堆插入新元素示意图

图3 删除根元素后并调整为最小堆示意图

图4 索引分段堆排序的数据结构示意图

图5 相位质量图（8×8）

图6 质量频率分布直方图

尽管堆排序相对传统方法效率提升不少，但是它在插入新元素和删除根元素后，调整为最小堆的过程依然比较费时。由此，本文提出了一种索引分段堆排序的方法，它相对单独的堆排序方法可以实现更高效的解缠速度。图4为索引分段堆排序的数据结构示意图。

索引分段堆排序的目的：将一个大型相位数据分割成多个L小堆，从而实现快速查找。根据总的质量值范围，分割成L个小范围，即每个小堆中的质量值都在一个固定的小范围内，它们属于一个l级别（索引号）。L等级的取值由人为决定，通常情况下，定义如下：

其中，M和N分别为干涉相位图行和列的大小。

例如图5是一个大小为8×8的质量图，其像素点质量值的范围为[1，32]，由上式知，L=8，它的频率分布直方图如图6，由质量值范围知，将每4个步长的质量值视为一个等级，所以将其合成为8个等级，图7为它的质量等级频率分布直方图。

图7 质量等级频率分布直方图

图8为索引分段堆排序查找示意图，首先，要查找不为空的最高级别l，并标识记录它。每当插入新元素时，确定该元素所属的级别并记为l′。如果l′＜l，则更新记录为l′，否则按堆排序步骤把该元素添加到对应等级的堆末尾处，然后将其调整为最小堆。因为每个最小堆中根元素的质量值最高，所以每次解缠的元素，就是记录不为空的最高级别为l的根元素，从该堆中移出根元素解缠后，也按堆排序步骤将其调整为最小堆。图9为本文方法的快速解缠流程图。

图8 索引分段堆排序查找示意图

图9 本文方法的解缠流程图

根据上述所说，索引分段堆排序将一个大型相位数据分割成L个小堆，相对单独的堆排序方法，其无论是在删除根元素，还是插入新元素后，调整为最小堆所花时间节省不少，计算效率更高。这在实际测量中，对于处理大型InSAR相位数据具有重大意义。

4 实验结果与分析

4.1 仿真实验

仿真的真实相位曲面由Matlab软件中的peaks函数产生，图10（a）为它的俯视图，且大小为256×256，实验中加入了高斯白噪声（信噪比为10 db）。图10（b）是图10（a）对应的缠绕相位图；图10（c）为相位导数方差图；图10（d）为本文质量图（QPDVC），由图（d）与图10（c）对比，可以看出图10（d）能更好地反映相位质量信息的好坏。图10（e）为10（c）对应的解缠结果；图10（f）为本文方法的解缠结果，与图10（e）相比，其解缠结果更接近于真实相位曲面，误差区域明显减少。

图10 仿真实验结果

为了能更好地分析本文方法的性能，通过加入不同信噪比的高斯白噪声，计算出本文方法、传统质量引导法（相位导数方差图）和枝切法解缠相位的均方根误差，如表1所示，从表中可以看出：本文方法的均方根误差最小，即解缠结果更精确且稳定性更好。

表1 本文方法与其他方法解缠精度对比

4.2 实测数据实验

实验所用数据为山西某地区的干涉相位数据，本文只截取了其中一部分来作为实验数据，其大小为256×256，如图11（a）所示。其相应的相位导数方差图如图11（b），本文相位质量图如图11（c）。为了验证本文方法的优越性，分别对枝切法、最小二乘法、网络流法、传统质量法和本文方法的解缠结果进行对比。

图11 实测数据的解缠结果

图11（d）～（h）给出了上述五种算法的解缠结果。表2是本文方法与其他方法解缠效果的对比。所有方法使用同一台计算机（Core i5 2.5 GHz CPU，2 GB内存），利用matlab R2013a来实现。

表2 同台计算机解缠效果对比

从仿真结果图可以看出，枝切法在残差点密集的地方，出现了解缠“孤岛”的现象，如图11（d）黑色区域所示，误差较大；最小二乘法的解缠结果最为平滑，如图11（e）所示，运行时间最短，但误差最大；网络流法有明显的误差传递，如图11（f），且部分边缘信息已丢失；而质量引导法通过借助可靠的质量图把不连续区域尽可能地控制在低质量区域，从而避免了误差的传递。因此，在可靠质量图的前提下，质量引导法的解缠精度比其他方法更高。

将图11（b）和图11（c）进行对比，可知本文质量图能更细微地反映质量图的好坏，如图11（g）所示，传统质量法解缠结果中白色方框标识的相位区域，出现了相位跳变，造成解缠结果不是很理想，而且运行时间最长；而本文方法与其对应的相位区域解缠结果没有出现相位跳变，如图11（h）所示，解缠结果理想，运行时间较短且误差最小。

因此，综合考虑解缠精度和运行时间，本文方法具有一定的优越性。

5 结论

本文提出了一种快速的质量引导解缠方法。质量引导法在可靠质量图的前提下，能更好地控制解缠的路径，其解缠精度比枝切法高。本文将QPDVC作为质量图，它比传统的伪相关质量图、相位导数方差图更加可靠。由于质量引导法在指导路径过程中需要大量的排序工作，通过索引分段堆排序法，减少堆排序过程中插入新元素和重新调整为最小堆所花时间，减少了解缠所需时间。

实验结果表明，本文方法解缠精度优于传统质量引导法和上述其他几种方法，且计算效率比堆排序法更高。该方法兼顾了解缠精度和运行时间两种性能指标，通过仿真实验相关数据对比，充分验证了本文方法的高效性和可行性。

[1]王超，张红，刘智.星载合成孔径雷达干涉测量[M].北京：科学出版社，2002：22-24.

[2]宋芳.干涉合成孔径雷达相位解缠[D].成都：四川大学，2005：1-2.

[3]Rosen P A，Hensley S，Joughin I R，et al.Synthetic aperture radar interferometry[J].Proceedings of the IEEE，2000，88（3）：333-382.

[4]Ghiglia D C，Pritt M D.Two-dimensional phase unwrapping：Theory，algorithms，and software[M].New York，USA：Wiley，1998：122-135.

[5]Zebker H A，Lu Y.Phase unwrapping algorithms for radar interferometry：Residue-cut，least-squares，and synthesis algorithms[J].Journal of the Optical Society of America A，1998，15（3）：586-598.

[6]Flynn T J.Consistent 2D phase unwrapping guided by a quality map[C]//International Geoscience and Remote Sensing Symposium，IGARSS’96，1996：2057-2059.

[7]Goldstein R M，Zebker H A，Werner C L.Satellite radar interferometry：Two-dimensional phase unwrapping[J].Radio Science，1988，23（4）：713-720.

[8]魏志强，金亚秋.基于蚁群算法的InSAR相位解缠算法[J].电子与信息学报，2008（3）：518-523.

[9]Pritt M D，Shipman J S.Least-squares two-dimensional phase unwrapping using FFTs[J].IEEE Transactions on Geoscience and Remote Sensing，1994，32（3）：706-708.

[10]龚国静，罗代升，何小海.无权重多重网格的最小二乘相位解缠研究[J].计算机工程与应用，2010，46（4）：166-169.

[11]Costantini M.A novel phase unwrapping method based on network programming[J].IEEE Transactions on Geoscience and Remote Sensing，1998，36（3）：813-821.

[12]Chen C W，Zebker H A.Phase uwrapping for large SAR interferogram：Statistical segmentation and generalized network models[J].IEEE Transactions on Geoscience and Remote Sensing，2002，40（8）：1709-1719.

[13]Xu W，Cumming I.A region-growing algorithm for InSAR phase unwrapping[J].IEEE Transactions on Geoscience and Remote Sensing，1999，37（1）：124-134.

[14]Fornaro G，Sansosti E.A two-dimensional region growing least squares phase unwrapping algorithm for interferometric SAR processing[J].IEEE Transactions on Geoscience and Remote Sensing，1999，37（5）：2215-2226.

[15]李芳芳，占毅，胡东辉，等.一种基于质量指导的InSAR相位解缠快速实现方法[J].雷达学报，2012（2）：196-202.