1 引 言

合成孔径雷达干涉测量(interferometric synthetic aperture radar,InSAR)凭借其独特的地面分辨率高、覆盖范围广、基于面观测等技术特点，在应用于数字高程模型(digital elevation model,DEM)的生成、地表形变的监测等方面已表现出巨大的潜力^{[1, 2, 3]}，InSAR正逐渐发展为一种重要的空间对地观测新技术。

InSAR技术是通过处理和分析覆盖同一地区的两幅或多幅SAR影像，从多个信号混叠的干涉相位数据中，分离和提取出与地形起伏、地表形变相关的信号分量，从而解算出地面高程或地表位移信息^{[4, 5, 6, 7]}。InSAR数据处理主要包括影像配准、重采样、相位干涉、参考相位的去除、相位解缠和地理编码等，每个处理环节的精细程度对于InSAR最终测量成果的精度和可靠性都具有极其重要的影响。

相位解缠是InSAR数据处理中的一个关键环节，也一直是SAR干涉处理中的难点和重点。相位解缠是将干涉相位中丢失的相位整周数采用几何解析或数值分析的方法恢复出其相对整周值，从而完成地形相对高差的解算或形变的检测^{[8, 9]}。相位解缠的准确性将直接影响到InSAR生成数字高程模型或地表形变检测的精度。在InSAR相位解缠中，如果相邻像素的相位差满足连续的条件，即相邻像素相位差绝对值小于π，则相位解缠很容易正确完成。然而，在真实的SAR影像干涉图中，造成相邻相位不连续、相位周跳的因素较多，如雷达阴影、失相关噪声、影像配准等因素引起相位数据不连续，这给相位解缠工作带来极大的困难，InSAR相位解缠仍然是一个十分棘手的问题。

目前相位解缠算法较多，主要归结为基于路径控制的积分法和基于最小二乘的整体求解法^{[10, 11, 12, 13]}。积分法的思路是计算像素行向和列向的一阶相位差分，然后对一阶差分连续积分即可求得解缠相位。由于干涉相位图存在奇异点，积分路径应受到约束以避免局部噪声相位误差的传播，这种算法的关键是按一定的原则对奇异点定位并予以连接作为积分路径的“防火墙”，即积分时不能穿越这些路径。最小二乘算法的思想是解缠后的相位梯度与缠绕相位梯度差异平方和为最小意义下的整体求解，使用带权估计方法可削弱奇异相位对解缠结果的影响。除此两种主要方法外，基于网络流的相位解缠算法、遗传算法、神经网络法等，国内外也有过相应的理论研究与应用尝试^{[14, 15, 16, 17, 18, 19]}。

对于干涉图残差点较少的区域，采用基于路径控制的积分法可以得到较为可靠的结果，而对于残差点较多的干涉图，如大多数非Tandem模式、有一定时间间隔的SAR影像对，由于很难正确设置路径控制线(或枝切线)，积分方法很难得到可靠的解缠结果，而最小二乘法则通过使缠绕相位梯度与真实相位梯度之差的平方和为最小来实现相位解缠，它对于残差点较多的区域仍能有效解缠，但其明显不足之处是，相位解缠的同时将误差扩散到了非噪声区域，且解缠相位与缠绕相位之差不是整周值，解缠相位不能满足一致性条件。

本文在传统最小二乘相位解缠算法基础上，提出了一种基于边界推进收敛加速的最小二乘整周相位解缠(phase unwrapping based on least squares with integer ambiguity,Publesia)方法，它通过对干涉图解缠边界的逐步推进和解缠相位的整数化迭代计算，较好地解决了最小二乘收敛性慢和整周值的求解问题。

最小二乘相位解缠是基于解缠相位梯度与缠绕相位梯度之差的平方和为最小的准则来实现相位解缠，它是数学意义上全局范数最小的求解问题^{[10, 11, 12]}，其目标函数可表示为

点击浏览公式

式中，Δ_i,j^x、Δ_i,j^y表示像元(i,j)在行向和列向的相位梯度值；Φ_i,j表示像元(i,j)的缠绕相位值；φ_i,j表示解缠相位值；W表示取相位主值的函数；上标x、y表示干涉图的距离向(range)和方位向(azimuth)；下标i、j分别表示像素的行向和列向坐标。

由于最小二乘相位解缠的目标是使缠绕相位与解缠相位梯度差的平方和为最小，已有研究表明，这种最小二乘法相位解缠可以归结为求解下式离散泊松方程的过程^[13]

点击浏览公式

令

ρ_i,j=Δ_i,j^x-Δ_i-1,j^x+Δ_i,j^y-Δ_i,j-1^y

则有

φ_i+1,j+φ_i-1,j+φ_i,j+1+φ_i,j-1-4φ_i,j=ρ_i,j

该泊松方程满足如下条件方程

点击浏览公式

式中各符号含义与式(1)相同。上述方程为无加权的最小二乘解法，其缺点在于解缠过程中它穿过残差点区域而并非绕过它们进行解缠，这可借助为相位数据引入权重值来加以克服，即加权的最小二乘相位解缠算法。对于那些因噪声或其他因素跳变的相位，可以将其权重赋为0，使它们不影响全局解缠过程。如果以InSAR中的相位质量图或者相干系数图作为相位数据的权重ω_i,j，则有0≤ω_i,j≤1，那么无加权的最小二乘就转化为加权的最小二乘算法，其改进的加权方程式如下^{[10, 11, 12]}

点击浏览公式

式中，ω_i,j^x和ω_i,j^y分别是像素在x和y方向的加权值，其余符号含义同式(1)。

按照上述的传统最小二乘算法，经过多次高斯赛德尔迭代所求得的最小二乘解具有很好的连续性和平滑性，但该解在任意一点的结果可能有噪声误差扩散，并且解缠相位与缠绕相位的差不是整周值，解缠相位数据不能很好地满足一致性条件。此外，该方法还存在迭代收敛较慢且在某些区域不收敛的问题，算法的收敛性有待提高。

根据最小二乘方法在迭代过程中逐步趋近于最优解的思想，笔者提出探测解缠边界的推进式解缠算法，将解缠区域划分为已解缠区域和待解缠区域，对解缠区域的边界逐步分析判断，对每次迭代计算的解缠相位值附加相位增量以正确提取整周相位值，通过逐步推进解缠边界，完成全局相位的一致性解缠处理。具体计算方法如下。

传统最小二乘算法每次计算返回的相位值可由下式计算得到^{[10, 11, 12]}

点击浏览公式

式中，φ_i,jⁿ表示像元(i,j)在第n次迭代计算后的解缠相位值。将同一像元两次迭代的解缠相位之差定义为Δ=φ_i,jⁿ-φ_i,j^n-1，根据Δ的具体数值判断解缠的趋势，设置一个较小的判别阈值T_C，阈值大小的有效性应充分顾及迭代计算的稳定性与时间效率的平衡。然后采用如下关系式对本次解缠的φⁿ_i,j进行计算

点击浏览公式

式中，dφ是为加快迭代收敛速度而设置的一个相位增量值，其值的大小根据试验结果获得，使解缠相位能更快地趋近于最优值。将前一次迭代计算的结果φ_i,j^n-1作为次优解，则本次计算得到的φ_i,jⁿ作为最优解，再对解算结果φ_i,jⁿ进行如下的相位整周计算。

设φ_i,jⁿ为第n次迭代解缠后的相位值，φ_i,j⁰为原始干涉相位值，则有

点击浏览公式

式中第1项Int为提取解缠相位整周值的函数。上述整周解缠计算能逐步正确解算出每个像素点的缠绕相位整周值，确保全局解缠相位的一致性。

基于上述流程，在每次迭代结束后还需进行边界条件的判断，以此判别已解缠相位和未解缠区域，采用如下判别式

点击浏览公式

式中，φ_i,jⁿ表示第n次迭代计算出的解缠相位值；φ_i,jⁿ⁺表示在第n次迭代的基础上再按式(5)进行计算所得的新相位值。设T_η为解缠判断阈值，若η＞T_η，则判断该像素未解缠，反之则认为已解缠。按此方法可探查已解缠与未解缠的边界像素，并将其进行标记。判断解缠边界后，由于仍无法直接判定边界两侧的已解缠和未解缠区域，但可判定边界上的点仍未得以解缠，然后将边界点的相位值都加上一个相位增量δ。则对于η＞T_η，有

点击浏览公式

式中，Δ为同一像元两次迭代解缠相位之差；相位增量δ的绝对大小可取小于π的一个固定相位数值，其有效性应充分顾及迭代计算的稳定性与时间效率的平衡。

依据上述计算流程，经过少量次数的迭代即可将解缠边界上的像素点逐一相位解缠，解缠边界不断向前推进，直至最终解缠边界像素数为零，迭代终止，完成全局相位解缠。整个算法的基本流程如图 1所示。

图 1 Publesia算法流程图 Fig. 1 Flow chart of Publesia method

图选项

需要提及的是，在本算法中考虑了相邻两次迭代相位差与附加相位增量对相位解缠的加速过程，使得算法的时间效率相当高，并且解缠迭代计算中的整数化过程保证了解缠相位的一致性和整周值的正确求解，较好地解决了传统最小二乘算法非整周性解算和收敛速度较慢的问题。

为验证本文相位解缠算法的可行性和可靠性，笔者基于 MATLAB语言工具开发了Publesia相位解缠程序，并使用模拟数据和真实InSAR干涉图对所开发的算法和程序进行验证。

模拟干涉图的主要干涉参数为：基线水平角α=0°，有效基线长度B_⊥=150 m，干涉相位的模糊度高为59 m，雷达波长5.7 cm，影像区域大小为100像素×100像素，区域内高程分布为0~200 m，干涉图中分布有20%区域的低相干区域(复相关系数≤0.2)。

该模拟区域的地形绝对相位如图 2所示，相位分布范围为0~30 rad，在上述干涉参数条件下得到的缠绕干涉图如图 3。干涉图的四角边界区域视为完全失相干区域，通过设置其权值为零以避免这些像素点对全局相位的解缠干扰。采用本文前述的最小二乘整周相位解缠算法(Publesia)对干涉图(图 3)进行相位解缠，并设置相位差阈值T_C=π/10，相位增量dφ=π/6。迭代计算至第4次、6次、8次、10次、12次和14次得到的解缠相位图如图 4(a)~(f)。

图 2 模拟地形绝对相位图 Fig. 2 Simulated topographic phase

图选项

图 3 模拟干涉相位图 Fig. 3 Simulated interferogram

图选项

从图 4中的(a)~(f)可以看出，随着迭代次数的逐次增加(4~14次)，相位解缠的边界正如前述算法中所描述的过程，逐步向前推进，当迭代计算进行到第14次时，解缠相位图与原真实相位图完全一致，正确解缠率达到100%，验证了该算法的正确性。

图 4 边界逐次推进的相位解缠过程 Fig. 4 The process of phase unwrapping by Publesia

图选项

然后使用含有噪声干扰的仿真数据来验证算法的可行性与可靠性。对上述图 3的干涉相位图分别施加弱噪声和强噪声干扰，噪声方差分别为0.27 rad和0.73 rad，采用本文的Publesia方法对弱噪声干涉图和强噪声干涉图分别进行相位解缠。在弱噪声条件下(噪声方差为0.27 rad)，解缠结果与真实相位仍然保持高度一致性，正确解缠率能达到100%。在强噪声条件下(噪声方差为0.73 rad)，逐次迭代计算得到的干涉相位图如图 5所示，(a)~(f)分别表示迭代第4次、6次、8次、10次、12次和15次的计算结果。

图 5 附加强噪声的相位解缠图 Fig. 5 Unwrapped phase from noisy interferogram by Publesia

图选项

为定量评估相位解缠结果的可靠性及解缠精度，笔者对解缠结果(第15次迭代计算)相位值进行数据统计，将解缠相位与原真实相位值进行差异比较，并进行精度评定。相位差异Δ_Φ的统计数据如表 1所示，并采用如下公式计算解缠相位的均方根误差(RMSE)

点击浏览公式

式中，m_φ为解缠相位的均方根误差；φ_unwrapped为采用本文方法解缠所得到的相位值；φ_truth为相位真值。

从表 1中的统计数据与精度评估可以看出，采用本文Publesia方法进行强噪声干涉图相位解缠，全图中约占99.5%的像素其解缠相位值与真实相位值的绝对差异小于0.1 rad，仅占约0.5%的像素其解缠相位值与相位真值差异超过0.1 rad，表明本文解缠算法对于强噪声(0.73 rad)干涉图仍然能够实现正确解缠，解缠精度达到±0.005 rad。

在采用模拟干涉图验证了算法的可行性和可靠性之后，笔者进一步使用欧洲空间局(ESA)的ERS-1/2卫星在香港和深圳地区获取的两幅单视复数影像(SLC)开展雷达干涉试验，SAR影像的局部覆盖范围如图 6所示。这两幅影像获取时间分别为1996-03-18(主影像)和1996-03-19(从影像)，该干涉影像对的有效空间基线长度为98.6 m，模糊度高约为103.9 m，采用DORIS软件进行干涉数据处理^[20]，得到图 6中方框区域的干涉相位条纹图，如图 7所示。

图 6 SAR影像振幅图 Fig. 6 Amplitude map of SAR image

图选项

采用本文解缠算法对图 7的干涉相位开展相位解缠，经迭代计算后得到的解缠相位如图 8所示。为评价本文解缠算法对该真实干涉图相位解缠的精度与可靠性，进行了解缠结果的外符合精度检查，采用斯坦福大学的Curtis W. Chen所开发的Snaphu解缠程序对图 7的干涉图开展相位计算。Snaphu方法充分考虑相位解缠的地形信息与形变提取物理意义，基于极大验后估计费用函数，利用统计特性寻找全局优化解，Snaphu方法对不同类型数据解缠均具有较高的可靠性^{[14, 21]}。因此，当难以获得解缠真值时，选用Snaphu计算结果作为本文解缠算法结果的评判依据，将本文Publesia计算结果与Snaphu解缠图进行差异比较，相位差异统计数据如图 9所示。

图 7 InSAR干涉相位图 Fig. 7 InSAR interferogram derived from DORIS

图选项

图 8 解缠相位图 Fig. 8 Unwrapped phase by Publesia

图选项

图 9 解缠相位差异统计 Fig. 9 The statistics of unwrapped phase difference

图选项

图 9中，横轴表示解缠相位差异区间(单位: rad)，按相位差异区间统计，在(-1.0,1.0)内统计单元长度为0.1 rad，在差异值超过该区间后，将(-∞,-1.0]内的像素统一归算到-1.0 rad处，[1.0,+∞)区间内像素统一归算到1.0 rad处，纵轴为像素个数的百分比。从图 9中的统计数据可以看出，采用本文Publesia方法对ERS-1/2 SAR干涉图进行相位解缠，图中有93%像素的解缠相位与Snaphu方法解缠相位的绝对差异小于0.5 rad，以Snaphu方法解缠相位作为参考，计算两种方法相位解缠的均方根误差(RMSE)，得到解缠精度为±0.12 rad。为进一步验证本文方法的解缠精度，采用加权的最小二乘解缠方法对图 7进行了相位解缠计算，解缠结果与本文方法计算结果的均方根差(RMSE)达到±0.28 rad，表明本文方法与Snaphu方法的解缠结果具有较高的一致性，而加权最小二乘解缠方法的计算结果偏差相对较大。

此外也需提及的是，图 9中Publesia方法与Snaphu方法的解缠结果也存在相位差超过0.5 rad的少量像素，产生这种差异的原因主要来源于解缠方法的差异，Snaphu方法是基于统计费用流的非线性优化约束算法，而本文方法是在全局优化基础上顾及了整周模糊度参数的一致性解算，有效避免了相位质量差的像元引起的全局误差传播，本文方法对于强噪声干涉图相位解缠具有较好的稳健性。

此外，对于解缠时间效率的比较，将Publesia方法、Snaphu方法和加权的最小二乘解缠方法，分别在相同的计算环境下予以测试，Snaphu方法对图 7干涉图相位解缠所耗时间为5.0 s，本文Publesia方法所耗时间为5.7 s，而加权最小二乘算法所耗时间大于30 s，这表明本文Publesia方法与Snaphu方法解缠效率相当，而相对于加权最小二乘方法，解缠效率有明显提高。

本文针对InSAR最小二乘相位解缠方法存在的收敛性差和整周模糊度解算的问题，提出了将干涉图逐步划分为未解缠区域和已解缠区域，对迭代计算的解缠相位值附加相位增量以正确提取整周相位值，通过逐步探测解缠边界以此推进相位解缠过程，最终完成全局干涉相位的一致性解缠处理。

基于所开发的Publesia相位解缠方法，对模拟干涉图和ERS-1/2卫星SAR干涉图分别开展相位解缠试验，在施加强噪声干扰的模拟干涉图中相位解缠的正确率达到99.5%，解缠的均方根误差为±0.005 rad，试验结果验证了本文解缠方法的可行性和可靠性。在ERS-1/2卫星SAR影像覆盖香港和深圳地区的干涉图相位解缠试验中，有93%的相位解缠数值与Snaphu方法解缠结果保持较高的一致性，且本文方法对于强噪声干涉图相位解缠具有较好的稳健性，解缠效率得以显著改善。