自20世纪90年代，合成孔径雷达干涉测量(InSAR)因其高分辨率、全天候、无需地面辅助等优点，已经成功应用于火山喷发、冰川漂移、滑坡、泥石流、地震、矿山沉陷等地质灾害的地表形变监测^[1-4]和全球或者区域的DEM获取^[5]。常规的差分干涉图中，强度反映地表散射特征的变化，相位记录地表与传感器之间距离变化量中不足整周相位的部分，且包含有形变、轨道误差、大气效应、失相干等多种成分，恢复真实形变是提取形变相位的基础。一个周期的解缠误差对形变的贡献量为λ/2(λ为传感器微波波长)，对于C波段而言，反映在干涉图上即是2.8 cm的粗差。因此，解缠误差作为InSAR技术中一个关键问题已不容忽略。

InSAR相位解缠是由二维平面上的缠绕相位恢复到真实相位的过程。相位梯度大于λ/2时，解缠误差会传播到后续的积分路径中，因此，解缠误差常呈片状分布。目前，国内外学者已提出多种解缠方法^[6-10]，根据解缠求取的方式可分为3类：路径跟踪法^[6-9]、全局最小范数法^[10-11]、网络流算法^[12]。路径跟踪法是根据任意闭合路径的相位积分值为零，确定并绕过残差点求解所有像素点的相位，路径的设置方式包括避开残差点^[6]、相位质量图引导路径^[7]、相干图指导枝切线^[8-9]等；全局最小范数法是以缠绕相位和解缠后相位差值最小达到干涉图的全局最优、确定像素点相位；网络流算法是根据相干性确定权重，以网络最优化为原则获取解缠相位^[13-14]。然而，实际解算中，由于InSAR侧视成像的特点，会有叠掩、阴影、收缩现象，研究区域复杂多样的变化特征产生的时空失相干及随机噪声会导致解缠误差的存在。

如何在最大化地保留干涉图的正确解缠相位信息的基础上探测干涉图的解缠误差，许多学者已有研究成果，主要分为两类：1)数据解算中处理。常规小基线技术中，构建包含有地表形变及DEM误差的函数模型，在参数解算过程中，含有较大残差相位的像元被认为含有解缠误差，从而作为粗差滤除。永久散射体技术中，选择裸露岩石等长时间内保持强度一致的相干点组成三角网，因弧段间距离较短，大气误差可忽略，如果函数模型残差值大于设定的阈值，则认为含有解缠误差，并将含有解缠误差的弧段作为粗差去除，同时解算弧段间的形变及DEM误差，通过网平差来最终获取全局相干点的形变量和DEM误差^[15]。2)干涉图预处理，探测出解缠误差。根据干涉图解算其闭合环残差相位，如闭合环残差较大且与周围残差存在2π的整数倍跳变，则认为含有解缠误差，手动识别解缠误差所在的干涉对^[16]，改正解缠误差。将解缠误差仅作为粗差从构建的函数模型中识别时，没有考虑大气误差在干涉图的空间分布，易造成解缠误差的不准确判断。仅用单个闭合环残差检查解缠误差的预处理方法并不能确定误差所在的干涉图，手动选择增加了工作量，而且受部分低相干区域的影响，解缠误差并不总是呈现明显的块状分布，或者干涉对达到上百时，逐个手动检查来定位解缠误差的位置和范围并不现实。

通过以上分析可知，在预处理阶段有效定位解缠误差并最大化地保留干涉图中未受解缠误差影响的相位信息显得尤为重要，而现有方法并不能满足此要求。受闭合环残差思想的启发，本文统计多个闭合环残差，可有效地定位解缠误差。现今多个载有InSAR传感器的卫星的成功发射及运行，提供了大量的SLC数据，更有利于本方法的广泛应用。

1 理论部分

单视复数影像(single look complex，SLC)经配准、干涉、差分、滤波后可获取缠绕的差分干涉图，经过解缠操作，干涉图中的相位表示卫星传感器与地面目标距离的变化量，其构成如下：

$ {\varphi _{ij}} = \varphi _{ij}^{{\rm{disp}}} + \varphi _{ij}^{{\rm{res, orb}}} + \varphi _{ij}^{{\rm{res, dem}}} + \varphi _{ij}^{{\rm{res, unw}}} + \varphi _{ij}^{{\rm{atm}}} + \varphi _{ij}^{{\rm{noise}}} $

(1)

式中，φ_ij指在时间点i和j获取的SLC处理得到干涉图，φ_ij^disp表示地表在视线向的形变量，φ_ij^{res, orb}表示因传感器位置不精确引入的残余轨道误差，φ_ij^{res, dem}为用于地形改正的外部DEM误差，φ_ij^{res, unw}为干涉图上的解缠误差，φ_ij^atm为干涉图上的大气误差，φ_ij^noise表示时间失相干和传感器热噪声等随机误差的综合影响。

假设有φ_ij、φ_ja、φ_ai3个构成闭合环的干涉对，其组成的闭合环残差为：

$ \begin{array}{l} \varphi _{ija}^{{\rm{raw}}} = {\varphi _{ij}} + {\varphi _{ja}} + {\varphi _{ai}} = (\varphi _{ija}^{{\rm{disp}}} + \varphi _{ija}^{{\rm{atm}}}) + {\rm{ }}\\ \;\;\;\;\;\;\;\;\;\;\varphi _{ija}^{{\rm{res, orb}}} + \varphi _{ija}^{{\rm{res, dem}}} + \varphi _{ija}^{{\rm{res, unw}}} + \varphi _{ija}^{{\rm{noise}}} \end{array} $

(2)

式中，φ_ija^disp为闭合环残差中的形变相位分量，φ_ija^disp=φ_ij^disp+φ_ai^disp+φ_ja^disp≡0；φ_ija^atm为大气相位残差，φ_ija^atm=φ_ij^atm+φ_ai^atm+φ_ja^atm≡0；φ_ija^{res, orb}为残余轨道误差，φ_ija^{res, orb}=φ_ij^{res, orb}+φ_ai^{res, orb}+φ_ja^{res, orb}；φ_ija^{res, dem}为DEM相位分量，φ_ija^{res, dem}=φ_ij^{res, dem}+φ_ai^{res, dem}+φ_ja^{res, dem}；φ_ija^{res, unw}为解缠误差分量，φ_ija^{res, unw}=φ_ij^{res, unw}+φ_ai^{res, unw}+φ_ja^{res, unw}；φ_ija^noise为随机噪声，φ_ija^noise=φ_ij^noise+φ_ai^noise+φ_ja^noise。

实际上，在干涉对组成的闭合环残差部分，大气相位和形变误差仅与SLC获取时刻相关，因此在闭合环中恒为零，残差结果仅含有后4项。φ_ija^{res, dem}与φ_ija^noise均表现出随机性，φ_ija^{res, orb}在干涉图上趋势性分布，φ_ija^{res, unw}在干涉图闭合差空间范围内表现为跳变。通过多项式拟合去除趋势面φ_ija^{res, orb}，则可得到φ_ija^{res, unw}所在区域：

$ \varphi _{ija}^{{\rm{corr}}} = \varphi _{ija}^{{\rm{res, unw}}} + n{\rm{ }} $

(3)

式中，φ_ija^corr为φ_ija^raw去除轨道误差的结果；n是随机误差，含有φ_ija^{res, dem}和φ_ija^noise。

将闭合环残差通过多项式估计去除轨道误差。采用多项式p=b₀+b₁x+b₂y+b₃xy，通过最小二乘求解多项式参数，并改正闭合环残差：

$ {\rm{res}}_{ija}^{{\rm{corr}}} = {\rm{res}}_{ija}^{{\rm{raw}}}-({\widetilde b_0} + {\widetilde b_1}x + {\widetilde b_2}y + {\widetilde b_3}xy) $

式中，$ {\widetilde b_i} $为最小二乘估计的多项式系数。

假设闭合环φ_ija^corr不存在解缠误差像素，在组成闭合环的干涉图φ_ij、φ_ja、φ_ai中也不存在解缠误差，即如果|φ_ija^corr| < ε，满足条件的像素在干涉对φ_ij、φ_ja、φ_ai中均不存在解缠误差，仅含有随机性误差。ε为事先设定的阈值，一般取ε小于或者等于0.5×2π。

通过求解与待检测干涉图φ_ij组成的k个闭合环，提取所有闭合环中不含有解缠误差的像素。根据假设，φ_ij的这些像素均不含解缠误差，取所有“干净”像素的并集，可识别出干涉对φ_ij中的解缠误差所在位置：

$ \left\{ \begin{array}{l} \varphi _{ijA}^{{\rm{corr}}} = \varphi _{ij{a_1}}^{{\rm{corr}}} \cup \ldots \cup \varphi _{ij{a_k}}^{{\rm{corr}}}\\ \varphi _{ij}^{{\rm{corr}}} = \varphi _{ijA}^{{\rm{corr}}} \cap \varphi _{ij}^{{\rm{raw}}} \end{array} \right. $

(4)

式中，φ_ijai^corr为不含有解缠误差的闭合环；φ_ij^raw为原始的干涉图；φ_ijA^corr为以φ_ijai^corr组成的模板，仅保留不含有解缠误差的相位。获取的φ_ij^corr既去除了解缠误差，又保留了有效消息。

2 实验部分

欧空局的Sentinel-1目前由两个C波段星座构成：2014年发射的Sentinel-1A和2016年发射的Sentinel-1B。其延续了先前的ERS和ENVISAT SAR^[17]，运行周期6 d，为哥白尼计划提供了大量数据，50 m的重复轨道间隔设计可降低常规差分干涉图的空间基线、提高干涉图的时空相干性^[18]。在地表形变监测中，最常使用的是通过IW TOPS(interferometric wide swath mode terrain observation by progressive scans)宽幅模式获取^[19]的数据，也是本文使用的数据类型。

墨西哥城位于墨西哥盆地南部，为横贯墨西哥火山带上构造火山的内部单元，平均高程在3 000 m以上。因地下水抽取过度及地层疏松，墨西哥城区是世界上地表沉降量较大的区域之一，最大年沉降速率超过35 cm。由于断层的存在，地表的空间形变呈现非均匀性变化，干涉图上解缠时易产生粗差。图 1显示了墨西哥研究区域及影像分布(红色框选区域为解缠误差研究区域)。

在欧空局网站获取墨西哥城区的Sentinel-1A数据，时间跨度为2015-04-13~2016-09-10，总共34景SLC影像，组成560个干涉图(图 1)。基于GAMMA平台，在距离31和方位向6的多视处理获得分辨率150 m×150 m的像素，经自适应滤波、由外部DEM模拟地形相位、去除地形相位的贡献量、最小费用流解缠，获得干涉对。

选择20150413-20150730作为检验的目标干涉对，与其他干涉对组成闭合环的情况如图 2(蓝色线段表示要监测的目标干涉对)所示，共组成31个闭合环。此处仅随机选择一个闭合环，干涉对构成为[20150413-20150730, 20150413-20150531, 20150531-20150730]。因干涉对闭合环残差中具有明显的趋势项误差，需对其进行多项式拟合去除轨道误差，图 3(b)和(d)为去除轨道误差前后相位周数的直方图统计，去除轨道误差前后闭合环残差的均值分别为0.176周和0.00周，标准偏差分别为1.292周和0.527周。为了更方便地可视化解缠误差，图 3中将残余的相位除以2π，以相位的周数(周)为单位。

为了验证检验结果的正确性，将31个闭合环分为2组，分别解算每组干涉对闭合环相位不含有解缠误差像素的并集; 以并集为模板，获取目标干涉对未被解缠误差影响的像素集合，结果如图 4(单位rad)所示。用同样的方法，图 5(单位rad)呈现了20150531-20160302干涉图解缠误差的解缠结果，显示片状分布的解缠误差已经被检测到并滤除。

3 结果分析

1) 闭合环残差中趋势项误差仍然存在。趋势项主要是轨道误差在干涉图上的体现，由图 3中闭合环残差的直方图统计结果可知，干涉对上轨道误差并不仅仅是单个SLC轨道误差模型的线性组合。经多项式拟合改正后，闭合环残差的统计结果基本趋于正态分布。

2) 多个闭合环残差的统计结果对解缠误差的定位。图 4和图 5显示研究区域内含有两个形变场，并展示了20150413-20150730和20150531-20160302两个干涉对的改正结果对比。图 4的干涉图时空基线分别为108 d和90 m，形变场的区域和大小与历史研究结果相符，因时间间隔较短，并没有出现明显的解缠误差。图 5的干涉图时空基线分别为276 d和-118 m，相对于图 4而言，时空基线增加，相干性降低，机场区域监测出有明显突变的片状相位解缠误差，并且两组干涉对闭合环均探测到此误差的存在，表明了此方法的有效性。同时，此方法是逐像素应用，不考虑干涉对像素之间的相关性，片状区均监测出存在2π的整数倍跳变，进一步验证了此方法的正确性。

3) 墨西哥区域形变梯度较大与此方法的适用性。此研究区域形变量较大，并在空间上存在不连续形变，在常规解缠误差探测中，这种形变特点易与解缠误差的探测混淆。此方法不受研究区域大气误差和形变的影响，可有效地探测出解缠误差。

4 结语

1) 本文提出的闭合环残差统计方法可有效地探测出解缠误差，不受研究区域形变特征及大气误差的影响。研究区域存在有大量Sentinel-1数据，可形成较多的干涉对闭合环，在多种形变情况下适用性较强。干涉对闭合环个数越多，解缠误差的探测结果越准确。

2) 闭合环残差中趋势项的存在表明，轨道误差并不是在SLC函数模型之间的线性组合，其直方图统计结果也验证了去除残余轨道误差的重要性。

3) 闭合环相位探测解缠误差中，两组实验结果的一致性及片状解缠误差的分布均验证了此方法的有效性。

4) 解缠误差的探测并去除，极大地简化了后续形变的计算，可提高结果的准确性。此方法适用于解缠误差的探测，今后可进一步研究定位解缠误差后进行解缠误差的改正。