相位分块与拟合法结合的InSAR相位解缠算法

马靓婷，卢小平，余振宝

(河南理工大学 自然资源部矿山时空信息与生态修复重点实验室，河南 焦作 454000)

0 引言

合成孔径雷达干涉技术(interferometric synthetic aperture radar，InSAR)是在主动遥感基础下发展的对地观测技术，其利用同一目标区域的SAR复图像对共轭相乘得到干涉图，根据干涉图的相位值，计算出2次成像过程中信号产生的路程差，从而获取监测区域高精度的三维地形信息和微小的地形形变信息[1-2]。因全天时、全天候、方便快捷等优点，InSAR技术被广泛用来监测地面沉降、地裂缝、火山活动、地震形变等[3]。

相位解缠是InSAR处理的重要环节，其结果的优劣直接影响到地形测量的精度。受限于合成孔径雷达的成像和处理方式，直接利用影像获得的干涉图一般含有较大的噪声，局部相位残差点的增多会形成不可靠数据斑块，使该区域相位解缠出现漏解或错解，导致InSAR图像恢复失败，从而影响到形变提取的精度。因此，提高相位解缠精度是InSAR处理中提高形变精度的重要环节[4]。

相位解缠是通过在解缠路径上进行积分从而还原真实目标信息。当干扰因素少、相位质量高时，能很好地还原真实相位信息；当干扰因素多时，误差会通过积分进行积累与传播，得到的相位数据与真实数据会有较大的差异。现有相位解缠方法主要分为路径跟踪解缠法[5-6]、最小范数解缠法[7]和网格规划解缠法[8]3大类。路径跟踪法是通过设置合适的积分路径，将误差限制在一定区域内，防止相位误差全局传递，其涵盖了经典的Goldstein枝切法、质量图引导法、最小范数法等。枝切法是利用残差点的连接得到枝切线，最后沿着枝切线进行积分得到解缠结果，但枝切法易出现孤岛现象[5]。质量图引导法是在质量图的引导下确定积分路径，这种算法对干涉图质量要求较高[9-10]。最小范数法是将相位解缠转换成数学上的最小范数问题，其常用的是最小二乘法，但这种方法穿过残差点会造成误差的全局传递[11]。网格规划法是将相位解缠问题转化为求解费用流的网络优化问题，主要有最小费用流和统计费用流等，但这种方法噪声会沿着积分路径传递，使得解缠结果不理想[12-13]。

针对当前相位解缠算法在相位图像存在严重噪声时解缠效果较差的问题，本文提出了基于相位分区与拟合法结合的InSAR干涉图相位解缠算法。该算法：首先对缠绕相位进行分块，将相位在相同区间的相邻像素进行合并；块内像素属于同一包裹次数，通过线性拟合求取合适的补偿系数K进行块间解缠；最后利用曲面拟合的方法进行残差块解缠，合并所有解缠块得到最终的解缠结果。

1 基本理论

1.1 相位解缠绕

在进行InSAR数据处理时，将主辅影像共轭相乘并取相位信息即可得到复干涉条纹图，但通过共轭相乘得到的相位差，与真实相位差相差2Kπ，即真实相位与缠绕相位的关系如式(1)所示。

Φ=φ+2Kπ

(1)

式中：Φ表示真实相位；φ表示缠绕相位；K表示补偿系数，其值为整数。

相位解缠是从缠绕相位φ确定补偿系数K值，进而估计真实相位Φ的过程。相位解缠必须兼顾一致性和精确性。一致性是指解缠后的相位数据矩阵中任意两点间的相位差与其路径无关；精确性是指解缠后的相位数据能真实地恢复原始相位信息[14]。

1.2 区域生长算法

区域生长算法的基本思想是将具有相似性质的像素合并到一起。对每一个区域要先指定一个种子点作为生长的起点，然后将种子点周围邻域的像素点和种子点进行对比，将具有相似性质的像素合并起来继续向外生长，直到没有满足条件的像素被包括进来为止[15]。

2 本文算法

InSAR干涉图虽然经过了滤波处理，但干涉图中仍存在噪声，噪声使得缠绕相位出现残差点和低相干区域，使干涉图解缠结果不全面或解缠误差较大。基于分块与合并策略的相位解缠算法是一类高效的方法[16-18]，该类方法把整个缠绕相位图分成若干块并执行块间的相位解缠，合并所有块可得到最终的解缠结果。现有的解缠算法无法将误差和孤岛现象同时避免，因此本文提出针对InSAR干涉图的分区与拟合法结合的相位解缠算法。

2.1 相位分区

本文将相位区域扩张法(prelude)中的相位分区用于InSAR干涉图相位解缠中，当分区间隔较大时相位块内会出现缠绕的现象，分区间隔较小时相位块较为零碎，不利于残差块的识别，影响算法的效率和精度。因此，本文选择以π/3为区间长度，对InSAR干涉图缠绕的相位矩阵实行分区，将相位信息在同一个区域的相邻像素点合并为像素块。缠绕相位的相位区间都在(-π，π]间，根据选定的区间长度，将相位分为(-π，-2π/3]，(-2π/3，-π/3]，(-π/3，0]，(0，π/3]，(π/3，2π/3]，(2π/3，π]6个区间。遍历整个相位矩阵，对其中相位属于各个分区的像素点进行标记，如图1(a)所示，将属于同一分区且相连的像素合并为块，得到若干相位处于同一区间的像素块，如图1(b)所示。合并像素采用四邻域法，即像素周围的4个与其相邻的同区像素可以进行合并，如图1(c)所示，块内像素拥有同一包裹数K。设置阈值为50，像素块大于等于50为正常块，像素块小于50为残差块。

图1 相位解缠过程示意图

2.2 拟合法解缠

在理想状态下干涉图不存在噪声，相位梯度小于π，选择一个点作为起始点，可直接进行积分解缠，但现实中噪声、地形起伏和相干性较低等现象给相位解缠带来了困难。利用线性拟合求取补偿系数K值的方法对噪声不敏感，在区域生长解缠绕中利用临近的已解缠绕像素的相位信息求取K值可快速得到相邻像素的真实相位信息，且能够适应较大的噪声水平和相位的快速变化等情况[19]。

拟合法相位解缠是基于相邻像素真实相位变换不大于π，使用相位分区方法对相位图像进行分块，块内相邻像素的相位变化小于给定相位间隔，即相邻像素之间的相位差小于π/3且块内像素拥有同一包裹数K，将块间相位解缠问题转化为线性拟合法求K值。选取块间相邻的2列像素，如图1(d)所示，其2列像素的相位关系为式(2)所示。

Y=X+2Kπ

(2)

式中：X表示已解缠像素块的相位值；Y表示未解缠像素块的相位值；K为补偿系数，其值为整数。

根据式(2)可知k为

(3)

式中：k为线性拟合系数，认为其最接近的整数值为整数补偿系数K；Xi表示已解缠像素的相位值；Yi表示未解缠像素的相位值。

缠绕相位块之间的相位解缠采用区域生长方式进行。把块作为一个解缠处理的基本单元，块内像素拥有相同的包裹次数K，选择第一个块为起始块，距离起始块中心位置最近的块为生长块。生长块的解缠绕即是找到一个最佳的相位包裹次数K，选择起始块与生长块相临近的2组像素，进行线性拟合，找到最优的补偿系数进行相位解缠，并合并解缠后的2块，之后进行下一个块的解缠，直至所有正常的块完成解缠与合并。

2.3 残差块曲面拟合

当所有的正常块完成相位解缠后，利用曲面拟合的方法对残差块进行解缠。根据解缠后相位是连续曲面的原理，通过残差块周围像素的真实相位对残差块进行解缠，将残差块内像素点10×10窗口内所有经过解缠的像素进行最小二乘曲面拟合，将其结果作为残差像素的真实相位值，选择第一个像素为起点利用区域生长算法将残差块内的像素解缠完毕。

2.4 本文算法描述

本文所提相位解缠算法主要分3个部分，即进行相位分区，将像素点分为若干个块；进行块间的相位解缠；残差像素块解缠。其算法的具体实现有以下4步。

1)将输入的InSAR干涉图缠绕的相位矩阵实行相位分区，将相位信息在同一个区间且相邻的像素点合并为像素块，将像素小于50的像素块标记为残差块。

2)实现相邻块间的相位解缠。利用线性拟合法求取补偿系数K进行相邻块间的相位解缠，根据区域生长算法将正常块解缠完毕。

3)当正常像素块全部合并完成后，对标记的残差像素块进行曲面拟合解缠。

4)合并所有解缠结果，完成整个区域的相位解缠工作。

3 实验结果与分析

本文对模拟数据进行实验并用实测数据验证。针对InSAR相位的特性，从解缠的准确度和算法运行的时间对解缠结果进行评价，除主观视觉评价外，用均方根误差和算法运行时间对实验进行客观评价。

3.1 模拟实验

本实验为MATLAB仿真的地形图，用基于雷达传感器参数和轨道数据的方法模拟InSAR干涉图[20]。首先模拟无噪相位图，如图2(a)所示，加入噪声后进行相位缠绕形成缠绕相位，如图2(c)所示,并将其中50像素×50像素大小的范围加重噪声，以验证算法在高强度噪声下的效果。

图2 模拟数据示意图

将本文算法与Goldstein枝切法、质量图引导的路径追踪法、四向剪切最小二乘法解缠结果进行比较与分析。图3(a)～图3(l)为不同解缠方法结果。从目视效果看，Goldstein枝切法解缠结果在高强度噪声区域出现了孤岛现象，无论是曲面图和俯视图都存在明显的错误。质量图引导的路径追踪法解缠结果在噪声较弱、相位质量好的区域解缠效果较好，在噪声密集区域解缠结果不连续。四向剪切最小二乘法解缠结果相位曲面图较为连续，但据其误差直方图可知解缠结果误差较大。本文算法解缠结果与原曲面图相位较为一致，在噪声密集区域依然有好的效果且误差值较小。

图3 模拟数据实验结果

表1是对模拟数据解缠结果的定量比较。从中可看出，四向剪切最小二乘法运算速度最快，但误差全局传递造成其精度最低。质量图引导法均方根误差小于四向剪切最小二乘法，但其运行时间大于四向剪切最小二乘法。本文算法的解缠结果均方根误差最小，相较于四向剪切最小二乘法减少了44%，相较于质量图引导法减少了38%，且运行时间与质量图引导法相当。

表1 不同方法的模拟干涉图解缠结果定量比较

3.2 实测数据验证

为验证本文算法的有效性，利用山西平朔地区2011年12月17日和2012年2月27日的RadarSat-2实测数据进行验证。对其数据进行配准、干涉等处理得到其真实相位干涉数据，截取其中500像素×500像素大小进行实验处理，如图4(a)所示。

图4(b)～图4(h)显示了不同方法的解缠结果。从目视效果看，Goldstein枝切法出现了孤岛现象，存在明显的错误。质量图引导法在低质量区域解缠效果不佳，其解缠结果俯视图上出现多个解缠效果不连续的区域。四向剪切最小二乘法解缠结果较为连续，但误差较大。本文算法结果较为光滑避免了很多尖峰毛刺现象。

图4 实测数据实验结果

表2是对实测数据解缠结果的定量比较。可以看出，Goldstein枝切法运算时间最长且解缠失败，四向剪切最小二乘法运算时间最短但其均方根误差最大。质量图引导法运算时间较长，但其精度高于四向剪切最小二乘法。本文算法解缠结果精度最高且运算时间适中。

表2 不同方法的实测干涉图解缠结果定量比较

4 结束语

本文提出了一种相位分区与拟合法结合的相位解缠算法，该算法首先将获得的相位图像分为多个块，块内相位都在给定的相位区间内，把像素个数小于给定阈值的块归类为残余像素；然后利用区域生长的拟合方法，依次进行块与块之间相位解缠绕和残余像素相位解缠绕，合并后得到最终的解缠结果。通过对模拟数据进行实验并用实测数据验证，实验结果表明，该算法无论是目视效果还是定量指标分析均优于其他算法，减小了误差传播的范围，提高了相位解缠的精度，但该算法受相位连续性约束面对迭掩区域和复杂地形相位跳变问题时仍存在较大的改进空间。