干涉合成孔径雷达技术(InSAR)是当前空间对地观测领域的为获得高程信息的热门技术。它通过两部天线对同一地区的多重的单景成像,然后通过对多重的单景成像进行相应的干涉处理,然后通过对缠绕的相位进行解缠的处理之后,三维数字高程图(digital elevation models,DEM)就可以通过解缠后的相位结合干涉合成孔径雷达的模型从而更进一步得到地表形变等物理信息,测量原理来源于杨氏双缝干涉实验,其中,如何从缠绕的相位信息中获得与高程直接对应的解缠相位是InSAR处理的关键一步[1-7]。
相位解缠算法主要有路径跟踪算法[1]和最小范数法[4]。路径跟踪算法中的经典算法是枝切法,它通过在残差点之间建立路径来分割极性不为零的残差点,从而避免解缠时的局部误差传播到别的待解缠区域。质量图法也属于路径跟踪算法,但是不同于相位解缠的经典枝切算法,它是直接利用能反应缠绕相位质量的质量图,然后从质量图的高质量部分开始解缠缠绕相位,依次向低质量缠绕相位区域进行解缠。质量图法的明显缺陷是没有像经典枝切法一样考虑相邻残差点引入的误差,从而在低质量区域,解缠效果会非常不理想。此外,其他解缠算法在相位解缠过程中也都有一定的局限性。最小范数法[7]以全局误差作为算法的衡量指标,容易造成局部区域误差的偏移,从而影响其他高质量区域的解缠效果。
本文提出一种枝切法和质量图法相结合的一种合成算法。该算法不同于以往的解缠算法,不仅利用了枝切法中的残差点的概念,以残差点的分布为依据,优化质量图中阀值的确定方法[8],还以此指导质量图法的相位顺利解缠,同时为了进一步的提高算法的运行效率,在质量图法中引入了堆排序的概念,使质量图算法的效率得到改进和提高。
由于枝切法和质量图法的解缠都可以归属为路径积分的一个类别,所以这2种方法都很好地满足相位解缠的一致性,即展开相位重新缠绕后与原始缠绕相位有很好的一致。所以提出的新的基于枝切法和质量图法的融合算法就保证了相位解缠的一致性问题。通常,枝切法得不到相位解缠的完整解,而由于质量图法解缠过程没有考虑到“残差点”的概念,所以很容易在解缠过程中引入由残差点造成的相位误差。所以针对这2种算法各自的优劣,提出一种枝切法和质量图法相结合的新型算法,仿真和实测数据处理结果表明该方法很好地克服了枝切法和质量图法在相位解缠过程中的缺点,同时保留了两者的优点,相比于其他方法[9-12]具有更高的精度。
1 堆排序算法先简单的介绍下计算机的对象结构,堆数据结构。它是一棵完全二叉树[13]。由于每一次的存取分别对应堆的根节点和叶子节点,所以这棵完全二叉树有时候并不总是满的,对于最后一层的叶子节点,往往跟数据的插入操作有关,因此它的大小总是保持着动态的变化。现在讨论堆{r1,r2,…,rn}应该满足的性质,并且做一些传统堆排序算法的简单介绍和相关性质的整理:
(1) |
式中i=1,2,…,n/2。该定义是2种常见的堆结构,分别是最小堆和最大堆。本文使用的堆结构要满足“邻接列”的性质,即取出其中最大的质量对应的相位的位置,然后将该点邻接的相位的质量系数放入里面,进行大小比较,然后放到对应的位置上。通过需求描述的对比,可以得知我们需要的堆结构应该是最大堆,它的根节点的元素的值就是所有节点的最大的那个,取出的操作也非常简单。
2 一种基于堆排序的质量法快速实现本节通过仿真意大利某火山的干涉数据分析验证基于堆排序的相位解缠算法的整体效率,如图 1所示。
图 1(a)为做过去平地处理和相位滤波后的缠绕相位,该干涉数据的大小为426 pixel×432 pixel,图 1(b)为通过所选用的相位导数变化图生成的相干系数,然后通过堆结构的排序算法,得到的解缠结果如图 1(c)。
通过比较不同尺寸数据的解缠时间进一步的衡量加入堆排序的质量图法的效率问题,如表 1。
图像尺寸/pixel | 本文算法/s | 传统质量图法/s | 加速比 |
128×128 | 6.17 | 14.38 | 2.33 |
256×256 | 16.25 | 126.94 | 4.51 |
426×432 | 145.76 | 1138.38 | 7.81 |
表 1中,把传统的基于质量图的路径跟踪算法与新提出的基于堆排序的质量图算法做了对比和验证,通过对不同大小的干涉图的解缠绕,得出本文算法的运行时间,可以很明显地看到所提出算法的优势。本次实验的环境为基于联想Y460(英特尔 Core i3 M 330 @2.13 GHz,6 GB内存)和MATLAB R2014a来进行的。
3 合成相位展开算法本文提出的方法主要是基于路径积分的思想。该方法首先通过使用传统的枝切法,解缠出大范围的缠绕相位,然后针对枝切法无法解缠的孤立区域以及残差点处的缠绕相位,通过质量图法中的伪相干系数图为参考,对剩余未解缠相位进行进一步的解缠。具体内容和步骤介绍如下。
3.1 利用枝切法展开较大范围相位由复变函数的相关理论可知闭合回路的积分满足如下等式:
(2) |
式中Nr为闭合积分回路里面的残差点极性的代数和。当式(2)里面的闭合路径的残差点的极性为零,或者等效极性为零时,解缠的相位正确性与解缠路径无关。
但是这部分的解缠过程只针对的是残差点比较少的部分。但是对于残差点过于密集的区域,枝切线无法完全放置下去,并且容易形成孤立区域,使相位解缠的过程无法进行,所以这部分的缠绕相位是枝切法无能为力的,这部分枝切法无法解缠的相位可通过质量图法进一步解缠。
3.2 利用质量图法展开剩余区域利用枝切法解缠大范围缠绕区域后,通过质量图法解缠剩余的未解缠区域。图 2为本文算法示意图,其中图 2(a)中黑色表示枝切法解缠出来的相位,灰色为邻接相位。假如箭头指向为质量图法的基准相位,则解缠过程通过该基准相位向外蔓延。通过判断待解缠相位及其四周的解缠状况进行相应处理。如果是已解缠相位,则直接跳过该点并将其标记为已解缠相位,如果是未解缠相位,则将其四周各点放入邻接矩阵,并扫描质量最高的相位是否为已解缠相位,是则将该相位视为参考相位对待解缠点进行解缠,否则重复将该点四周相位放入邻接矩阵并重复操作。从图 2(a)到图 2(b)的过度是假设图 2(b)的箭头所指的相位对应的质量数在图 2(a)的灰色相位中为最高,因此该点先得到解缠,并且将该点及其周围未接触相位标记为灰色。
基本操作步骤如下:
1) 利用3.1中描述的枝切线法,布置枝切线并通过像元扩散法对缠绕相位进行积分来获得解缠相位,此时的解缠未完成,因为对于大规模缠绕相位,很容易形成闭合的枝切线而无法解缠。
2) 利用伪相干系数图的定义得到相位质量图,然后标记其中的已解缠相位,并且把这些相位对应的质量图中的质量数修改为0,因为对于伪相干系数图来说,数值越小,则表示其相应的相位质量越高。
3) 以步骤1)中基准相位为开始,对其周围的4个相邻点分别进行遍历,对其中质量最高的相位进行解缠,若其已在步骤2)中标记为解缠,则跳过该点,否则标记其为已解缠相位。然后将该点周围的4个相邻点中未被处理的点标记到“邻接矩阵”中。
4) 从“邻接矩阵”中找出质量值最小,也即质量最高的点,对其重复步骤3)的处理,这一过程是通过堆排序实现的,不在详述。
5) 重复步骤4),直到“邻接矩阵”为空,也即所有像素点都得到处理。
4 实验结果与分析本节实验拟通过从仿真数据和实测数据两个维度反应出本文算法在相位解缠过程中,对解缠过程指标的提高。这其中提高的指标反应在成像上的对比和精确度上的对比,从而验证本文算法的性能。
4.1 仿真实验这里的质量图选的是相位倒数变化图,然后数据是采用大小为像素点128×128的仿真相位,它由MATLAB软件中peaks实现,效果如图 3(a)。在仿真过程中,为了模拟现实数据的噪声的干扰,手动加入相干系数为0.9的均匀噪声。该相干系数对应的噪声方差为0.422 5 rad,它的分布区间假设为[-b,b),则可得b=1.125 8 rad[14]。对加入噪声的待处理相位进行缠绕处理后如图 3(b)。通过枝切法得到的残差点如图 3(c)所示。较高质量区域内缠绕相位的展开如图 3(d)。引入质量图进行相位展开后如图 3(e),从而得到最终解缠结果如图 3(f),其与图 3(a)的未缠绕相位几乎一致。
为了对本文算法进行定性地分析和比较,我们首先做出的是各类算法解缠后对应的解缠误差的二维分布比较。其中解缠误差的定义是该点展开相位与未加均匀噪声前的待处理的真实相位上对应点相位值之差的绝对值。与常用的解缠算法对比,效果如图 4。可以看到本文方法误差没有很集中,整体的分布范围可以通过色调偏暗得出。并且可以看出误差的幅度也不会太大;最小二乘法的误差分布因为算法本身的问题,其误差范围最大,并且在相位变化较大缠绕相位处的解缠误差明显较大。同本文算法主要是误差的范围较大,这也是最小二乘法的局限;网络流法的误差也主要分布于相位变化较大的区域,虽然误差范围较小,但是误差的幅度非常明显,这也是跟本文相比的不足之处。
最后,为了进一步的检验本文算法相较于最小二乘法和网络流算法的稳健性,我们通过改变加入均匀噪声的方差的幅值,然后计算出各类算法展开后所得到的均方根误差,见表 2。从表 2中可以看到:本文方法的误差值要小于其他的2种算法,而这充分说明了本文算法比最小二乘法和网络流算法有更好的稳健性。
噪声方差 | 对应的相干系数 | 本文算法误差 | 最小二乘法误差 | 网络流法误差 |
0.422 5 | 0.9 | 2.498 9 | 8.523 5 | 3.562 4 |
0.810 0 | 0.8 | 4.084 4 | 9.945 6 | 7.432 3 |
1.210 0 | 0.7 | 5.972 7 | 11.790 8 | 10.659 9 |
4.2 实测数据实验
实测数据采用的是文献[2]提供的某山地干涉相位数据,该数据为512 pixel× 512 pixel的缠绕相位,如图 5(a)所示;其对应的残差点分布如图 5(b);部分未孤立的相位通过枝切法展开得到的枝切线放置如图 5(c);部分展开相位如图 5(d);剩余的没有展开的区域留待进一步通过质量图法进行解缠,在图 5(e)中用黑色表示;利用相位倒数变化法得到的质量图如图 5(f);通过质量图法展开图 5(e)中未展开的区域;图 5(g)为最后的展开相位。为了直观地比较和分析本文算法对缠绕相位解缠的有效性和一致性,把最终得到的展开相位重新缠绕如图 5(h),该缠绕相位与原始缠绕相位非常吻合,充分说明本文方法能够恢复出完整且较为平滑的相位曲面,具有一定的现实可推广性。
5 结束语为了充分利用经典枝切法和质量图法的一致性和有效性的两大优点,本文提出了一种即结合了枝切法又兼顾到质量图法的新算法,该方法对于下面的问题有针对性的诠释。
1) 根据枝切法的算法特点,有效地避免了质量图法误差传递的缺陷,并利用质量图法的特点有效地解决了枝切法解缠不完全以及容易形成无法解缠的孤立区域的缺陷。
2) 结合堆排序算法进一步提高算法整体的运行效率。
3) 能够恢复出完整且较为平滑的三维曲面。
4) 提高了整体的缠绕相位的解缠精度。
[1] | GOLDSTEIN R M, ZEBKER H A, WERNER C L. Satellite radar interferometry-two-dimensional phase unwrapping[J]. Radio science , 1988, 23 (4) : 713-720 DOI:10.1029/RS023i004p00713 |
[2] | GHIGLIA D C, PRITT M D. Two-dimensional phase unwrapping: theory, algorithms, and software[M]. New York: Wiley, 1998 . |
[3] | WANG C, ZHANG H, LIU Z. Spaceborne synthetic aperture radar interferometry[M]. Beijing: Science Press, 2002 . |
[4] | PRITT M D, SHIPMAN J S. Least-squares two-dimensional phase unwrapping using FFT's[J]. IEEE transactions on geoscience and remote sensing , 1994, 32 (3) : 706-708 DOI:10.1109/36.297989 |
[5] | COSTANTINI M. A novel phase unwrapping method based on network programming[J]. IEEE transactions on geoscience and remote sensing , 1998, 36 (3) : 813-821 DOI:10.1109/36.673674 |
[6] | 刘双亚, 李世强, 冯锦. 改进的基于频谱分割的InSAR绝对相位确定方法[J]. 国外电子测量技术 , 2016 (7) : 27-33 |
[7] | 张金芝, 黄海军, 毕海波, 等. SBAS时序分析技术监测现代黄河三角洲地面沉降[J]. 武汉大学学报: 信息科学版 , 2016, 42 (2) : 242-248 |
[8] | 余应淮, 谢仕义, 梅其祥. 基于核回归修正的上采样相位相关精确运动估计算法[J]. 计算机应用 , 2016, 36 (8) : 2316-2321 |
[9] | 陈强, 杨莹辉, 刘国祥, 等. 基于边界探测的InSAR最小二乘整周相位解缠方法[J]. 测绘学报 , 2012, 41 (3) : 410-416 |
[10] | 钟何平, 唐劲松, 张森. 一种基于质量引导和最小不连续合成的InSAR相位解缠算法[J]. 电子与信息学报 , 2011, 33 (2) : 369-374 |
[11] | 李芳芳, 占毅, 胡东辉, 等. 一种基于质量指导的 InSAR 相位解缠快速实现方法[J]. 雷达学报 , 2012, 1 (2) : 196-202 |
[12] | 张小京, 刘冰. 数据共轭复用策略InSAR干涉相位估计[J]. 测绘科学 , 2016, 41 (8) : 28-32 |
[13] | 肖枫, 伍吉仓, 刘朝功, 等. 不同质量图在相位解缠算法中的比较分析[J]. 大地测量与地球动力学 , 2010, 30 (2) : 80-85 |
[14] | LEE J S, HOPPEL K W, MANGO S, et al. Intensity and phase statistics of multilook polarimetric and interferometric SAR imagery[J]. IEEE transactions on geoscience and remote sensing , 1994, 32 (5) : 1017-1028 DOI:10.1109/36.312890 |