基于区域识别和区域扩展的相位解缠算法

毕海霞 魏志强

（1.西安电子工程研究所，陕西 西安710100；2.西安交通大学信息与系统科学研究所，陕西 西安710149）

引 言

干涉合成孔径雷达（Interferometric Synthetic Aperture Radar，InSAR）技术是当前空间对地观测领域的研究热点，它通过两部天线或者重复轨道观测获得同一地区的影像.多景图像经过滤波、配准、干涉、去平地效应、解缠和地理编码等处理后，可反演出地面的三维数字高程图（Digital Elevation Models，DEM）或地表形变［1］，其中，相位解缠是InSAR处理的关键.

目前，相位解缠算法主要有路径跟踪算法和最小范数法［2-14］.枝切法［3］是典型的路径跟踪法，通过识别残差点并在它们之间建立分割路径，避免局部的误差传播到整个图像中去.该方法运算速度快，但对低相干区域易造成误差传播.质量指导法属于路径跟踪算法［4-5］，它以相干系数为指导，选择像素作为生长种子和待解缠像素，由种子的邻域像素逐步扩展到低相干的像素，对相干性过低的区域不做解缠.该方法具有解缠精度较高的优点，但它不判断种子像素的分布、不回避相位非连续像素，在低相干区域会引入累积的相位误差.此外，其他解缠算法也均有一定的针对性或局限性.最小范数法［6-7］以缠绕相位的梯度与解缠相位的梯度之间差最小为标准，建立相应的数学模型，以使得解缠后的相位接近于真实的相位.根据实际解缠的对象，选取不同的范数和缠绕相位梯度的权重等效于各种最小二乘相位解缠方法.该方法解缠精度高，整体相位误差小，但由于每个像素的解缠相位都是拟合值，不是真实值，对相干性差异较大的图像，容易造成局部区域的相位整体偏移.

实际的SAR图像中往往存在着多个相干性不同的特征区域，兴趣区域是整幅图像中的局部.例如，在被低相干的水域或植被环绕的混合区域，它们通过狭长地带（如桥梁、堤坝等）与其他区域相连.相位解缠中低相干区域的误差容易在整个图像中传播.为此，本文介绍了一种基于区域识别和区域扩展的解缠算法，该方法能够有效限制低相干区域的误差传播，提高相位解缠的精度.实验数据为日本对地观测卫星装载的相控阵L波段合成孔径雷达（Phased Array L-band Synthetic Aperture Radar，PALSAR）图像，选择香港大屿山附近区域为实验区域.首先，选取相同实验区域的光学图像，根据其光谱特征进行分类，以分离出低相干区域的光学图像作为分类模版对干涉图像进行掩模.然后，针对掩模低相干区域后的干涉图像，以高相干区域的像素作为种子向低相干区域进行邻域扩展解缠.最后，由已解缠的邻域像素确定低相干像素的模糊数，迭代检验，最终确定其解缠相位.为了进一步证明区域识别与扩展方法相对于质量指导区域生长法的优越性，选择伊朗Bam地区为实验区域，使用欧空局环境卫星装载的高级合成孔径雷达（Advanced Synthetic Aperture Radar，ASAR）数据进行比较，实验结果与比较证明了本文方法在解缠精度上的优越性.

1 相位解缠算法

1.1 区域识别方法

实验PALSAR图像包括香港大屿山岛屿和附近的海区.其中陆地区域是目标区域，海区是低相干区域.为避免海区的误差对研究区域的干扰，在干涉处理前先从InSAR图像中识别出兴趣区域.区域识别方法融合可见光图像（来自Google Earth数据库），对可见光图像进行监督分类，从中分离出低相干的海区，再以可见光图像对InSAR图像掩模，分离出InSAR图像中的海区.

监督分类方法以实际地物类别的先验知识为基础，根据每一类别的图像特征进行分类处理及判别函数计算.采用最大似然监督分类算法对可见光图像进行分类［15］，步骤如下：

1）在可见光图像上截取与InSAR图像大致相同的区域.以InSAR图像为主图像，可见光图像为辅图像，对可见光图像进行插值，配准到SAR图像.

2）分别在海区、山区和城区中选定多个典型区域作为监督分类的先验样本，如图1（a）所示.所选区域的大小和数量由特征区域的分布和光谱复杂性决定.Google Earth的光学图像由多幅不同光学性质的图像拼合而成，且每幅图像中的海区、植被覆盖区和城区等区域的光学性质也不同，所以，需选择多个具有代表特征的区域作为先验样本，并求其先验概率.

3）以光学图像中每个像素作为训练样本进行计算，进行统计意义上的归类.以像元间的光谱距离作为分类的依据，采用欧氏距离作为判别函数，对于每个训练样本进行以下计算：

①计算各先验样本R、G、B三个通道的中位数（Rm，Gm，Bm）作为类中心.

②计算各先验样本区域内像素（Ri，Gi，Bi）到其各自类中心三个通道的欧氏距离Di

③计算样本区域像素欧氏距离的标准差σ

4）计算整个光学图像像素到各个样本区域中位数的欧氏距离di，若其标准差小于预先设定基于标准差的阈值，则被认为属于这一允许归类区.

式中：A为阈值和标准差的倍数；置信概率σ可由归一化概率分布函数的反函数计算得到.训练过程中，由小到大逐步调整，将满足条件的像素初次归类于各自的允许分类区.

5）对少量剩余“模糊像素”，采用Bayes判断准则进行统计意义上的二次归类，即

式中：P（i）为类别的先验样本概率；P（X／i）为训练样本的似然概率，表示像素X在第i类中出现的概率.如果式（4）成立，则像素X属于类别i.设可见光图像像素为正态分布，根据Bayes判别规则将可见光图像上各个像素归属概率最大的类.

按上述步骤，反复迭代运算可分类出光学图像中的海区、植被覆盖区和城区等区域.针对研究对象，从图像中分离出海区，保留陆地区域进行下一步处理，如图1（b）所示.

1.2 区域扩展相位解缠算法

图像中的区域扩展处理是将成组的像素或区域发展成更大区域的过程.从种子点集合开始，按照一定的扩展规则，将符合条件的相邻像素合并到该区域，直到处理完所有相邻像素.InSAR图像中的相位解缠可采用区域扩展的策略［4-5］，干涉图像的相位解缠过程可视为由“缠绕相位”的像素逐步扩展为“解缠相位”的像素过程.

区域扩展相位解缠方法的关键步骤是种子的选取、待解缠像素的解缠次序和相位估算.本解缠方法选取一系列高相干的像素作为生长“种子”，按照相干性和邻域已解缠像素的“影响力”的原则选择和估算待解缠像素的相位.由已解缠的像素向其邻域未解缠的像素逐步增长，最终实现干涉图像相位的整体解缠.算法步骤如下：

1）选取种子.设定阈值，选择均匀分布的相干系数较高的多个像素作为生长种子.

图1 融合光学影像识别兴趣区域

2）解缠次序.解缠像素顺序不同会引起解缠结果的不同，通常由可靠性高的像素向可靠性低的方向扩展.与种子像素相邻的像素为待解缠像素，从中选择相干系数最高的和受到周围像素“影响力”最大的像素进行解缠.逐次降低阈值（选取迭代步长为0.05），相邻像素中可能会有多个像素符合条件，因此，从中选择与受到已解缠相邻像素“影响力”最大的像素作为待解缠像素.然后将已解缠的像素归为种子像素集合，更新待解缠像素集合.从每个种子同时进行，解缠区域由其相邻像素逐步向外扩展，直到边界.

待解缠像素由其24邻域像素中已解缠的像素的“影响力”决定，其模型要素如下：

①距离近的已解缠像素“影响力”强.8邻域已解缠像素（图中网格单元）的“影响力”大于其16邻域已解缠像素（图中斜线单元）的“影响力”，如图2（a）所示，图中间的黑色像素表示待解缠像素.

②相位的线性分布强于非线性分布.在已解缠的24邻域像素中，与待解缠像素在同一方向上的像素组（如图2（b）中箭头方向上两像素）对待解缠像素的“影响力”大于其他方向上的像素组.绝大多数相邻的已解缠像素的相位在统计意义上是连续分布的，即已解缠像素的相位梯度场是连续的.

图2 已解缠像素对待解缠像素的“影响力”模型

③“影响力”估算.利用矩阵运算求得待解缠像素的“影响力”，将待解缠像素和其24邻域像素看作5×5的矩阵P5×5或8邻域像素视为3×3的矩阵P3×3，P中已解缠元素的值为“1”，未解缠像素的值为“0”.设定“影响力”权重系数，以距离矩阵R5×5、R3×3和相位梯度矩阵G3×3表示.

计算图2（c）所示的已解缠像素对待解缠像素的“影响力”，黑色像素和其他白色空白像素为未解缠像素，网格像素和斜线像素均为已解缠像素.分别设24邻域像素中的16邻域像素、8邻域像素和梯度场的个数分别为N1，N2和N3，归一化的影响权重分别为ω1，ω2，ω3，且有ω1＜ω2＜ω3.待解缠像素受到24邻域中已解缠像素的“影响力”如下计算

式中，｜［］·［］｜表示两矩阵的分素乘积的模.影响力取决于ω1，ω2，ω3的取值和已解缠像素的数量.设定其数值分别为1／6，1／3，1／2.图中待解缠的像素受到邻域像素的“影响力”为19／6，高于图中其他已解缠像素的邻域像素，所以，图2（c）中黑色像素为待解缠像素.

3）待解缠像素的相位估算.相位解缠的数学过程实际是估算像素相位的整数倍周期数值，即模糊数m.定义未解缠像素的缠绕相位φwrapped和其解缠像素的解缠相位φunwrapped的模糊数为

式中，round（x）为取整函数，即求与x最邻近的整数.

待解缠相位的预测值由24邻域中已解缠像素相位的“影响力”决定.定义待解缠像素24邻域和8邻域中已解缠像素的相位矩阵分别为Hp5×5和HP3×3，矩阵每个元素为已解缠像素的相位值，未解缠像素的相位设为零.

初步估算的模糊数为

所估算的解缠相位为

4）对估算的解缠相位进行迭代检验.在单视复数InSAR图像中，在地形渐变的高相干区域，其相邻像素间的相位差变化较小.而在地形起伏大或随机散射严重的区域，其相邻像素间的相位差变化较大，例如，山区、海洋、植被覆盖区和密集城区等区域的SAR图像，其相位容易产生突变.因此，对解缠像素进行迭代检测时，应遵循以下规则：

①解缠后的相位是连续分布的.解缠过程中，优先检验相干系数高、相邻像素间相位差小的已解缠像素，再逐步向相邻像素位相差较大的、相干系数低的像素扩展.已解缠像素的相位与每个8邻域像素之间的相位差小于一个周期，所估算的模糊数与其8邻域的差不超过“1”.

②计算所估算的解缠相位与解缠的预测值间的差值为

差值小于指定阈值的解缠结果为有效的解缠结果，检验阈值由小到大逐步增加，逐渐放松阈值，迭代检验.直至检测完所有的像素.本文分别取π／4、π／2和π等.③对低相干区域不符合条件的像素，待邻近像素检验通过后，重新采用24邻域像素估算方法估计该像素的解缠相位值，直到所有像素被检验通过.

缠绕相位像素在以每个种子为中心进行扩展的过程中，随着每个种子所在区域的扩大，两个区域可能会重合在一起.当两个种子的增长区域超过3个像素重合时，将两个区域融合一起成为一个新的整体.融合的标准是使得两个种子的增长区域的模糊数差最小.

2 实验结果与比较

2.1 香港大屿山的PALSAR数据实验

本文的实验数据为香港大屿山岛附近的PALSAR数据.大屿山岛通过桥梁与马湾岛、青衣岛、新界相连，如图1（a）所示.该实验区陆地部分由东西走向的山体和香港国际机场组成，其干涉相位如图3（a）所示.下面分别使用本文方法和其他几种解缠方法对干涉图像进行解缠，反演试验区域的DEM，如图3所示.

图3 几种解缠方法的PALSAR数据实验比较

采用枝切法的解缠结果如图3（c）所示.海区噪声过大，残差点过于密集，造成了累积误差的突变.在其他低相干区域，出现了多个突变高程噪声值，尤其在机场中部和左下角，在机场与大屿山之间海湾和沿海的青屿干线附近海区被反演为与陆地相同的高程值.而右上角由桥梁连接的马湾岛区域（图3（c）中圆圈标识区域）被反演为低于海平面的高程值，桥面区域也被扩大.可见，没有经过区域识别产生的边界限制，噪声不仅影响低相干区域，也影响到了邻近的高相干区域.尤其是通过狭小区域（如桥梁）连接的岛屿，在沿桥面向岛屿积分解缠过程中，可能会造成整个区域的误差.

最小范数法反演的DEM如图3（d）所示.最小范数法采用整体拟合的方法虽然有很好的解缠视图效果，但其牺牲了图像的细节，使整体的DEM变得模糊.由于实验区域存在大面积海区和植被覆盖的山区等低相干区域，采用整体拟合的解缠方法易造成部分区域的DEM被低估，而另一部分区域被高估，如图3（d）中被低估的左上角海区和被高估的右下角区域，同时，也造成了邻近的迪斯尼公园区域（图3（d）中右下角圆圈标识区域）的DEM被高估.

质量指导法依赖于干涉图像的相干性，相干图像如图3（b）所示.图中低相干的区域（图中圆圈标识区域）与海区将左下角区域与其他区域隔离开，高相干系数的像素离散分布于机场建筑物、沿海道路、城区等.质量指导法自动选取相干系数最高的像素作为生长种子，对相干系数过低的区域不进行解缠，因此，设置解缠阈值后，大部分的海区无法得到解缠.另外，因相干性的差异过大，整个图像被分割成多个独立解缠的区域，反演后造成大量的高低错位倒置，图像无法进行整体处理.

本文方法选取相干系数大于0.7的像素作为扩展种子，从每个种子开始逐步按照相干系数和“影响力”规则向外扩展解缠，解缠结果如图3（e）所示.由于边界限制，海区部分没有参与解缠，图中海区的高程设为“零”.其他区域解缠的相位，按照PALSAR参数反演DEM.本文方法采用限制区域和由高相干区域向低相干区域扩展的策略，最大限度地控制了反演误差的传播，获得了较高的反演精度.

将上述实验结果与实验区域的DEM进行比较.该实验区域真实的DEM是由香港规划地政局通过航空摄影立体测量获取，分辨率为20m，如图3（f）所示.为了最大限度地体现各方法在反演陆地区域的差异，比较区域选为图3（b）中方框区域.将各个解缠算法反演的陆地区域的DEM与实验室区域的真实DEM逐点进行比较，其相对误差如表1所示.可以看出，区域识别与扩展解缠方法的反演精度高于其他方法.另外，如果选择整个图像中的大屿山陆地区域比较，各方法反演的DEM与真实的DEM相对误差将更大.

表1 各种解缠算法反演的相对准确率

2.2 Bam地区的ASAR数据实验

本文方法与质量指导法相近，均选取高相干像素作为种子开始生长，解缠像素由高相干区域向低相干区域生长.但这两种方法在待扩展邻域像素的选取、解缠相位的估算和初步解缠像素的检验均有所不同.质量指导方法在低相干区域易发生相位跃变和误差传播.另外，由于低相干区域的反演误差较大，质量指导方法对相干过低的区域不进行处理，而本文方法采用多次迭代策略，利用24邻域中已解缠像素的“影响力”估算剩余的低相干像素.

图4 伊朗Bam地区ASAR数据实验比较

为了进一步比较本文方法与质量指导法的优劣，本文采用另一组ASAR数据进行验证，数据区域为伊朗Bam城区.Bam为沙漠中的古城，图像中间和左上角的低相干区域为Bam城区，周边区域为沙漠和沙丘.其干涉图像如图4（a）所示.其相干系数图像如图4（b）所示，城区和山脊处的相干系数较低.分别采用质量指导法和本文方法对其进行解缠.质量指导法反演的DEM如图4（c）所示，本文方法反演的DEM如图4（d）所示.

两种方法在高相干区域的解缠结果是相近的，但在低相干区域，质量指导法出现了相位跃变和误差累计.其原因是：在待解缠像素的选择上，质量指导法只是根据相干系数的大小选择待解缠像素；而本文方法在选择待解缠像素时，不仅考虑了相干系数，而且考虑了其邻域已解缠像素的“影响力”.在估算解缠结果的检验上，质量指导法没有考虑相邻像素的整体连续性.针对单纯的相干系数指导的不足，本文方法在种子选取、待解缠像素的选择和估算等方面做了修正，从而提高了反演精度.

3 结 论

在InSAR图像解缠中，研究区域容易受到其他低相干区域误差的影响，而常用的相位解缠方法在限制低相干区域误差传播方面存在不足，为此，本文提出了一种基于“低相干区域隔离”和“高相干区域解缠像素扩展”的解缠算法.在区域识别中，由于光学图像分辨率高，且容易获取，将光学图像与SAR图像融合.根据光学图像上不同区域光谱特征的差异，识别出研究区域，并去除无关的低相干区域，再进行InSAR处理.为减少误差传播，解缠算法采用了区域扩展方法，由高相干区域向低相干区域扩展，保障了高相干区域反演结果的可靠性，提高了反演精度.在生长种子的选择上，本算法考虑了相干系数的大小和种子在图像中的分布均匀性.在待解缠像素的选取上，同时考虑待解缠像素的相干系数和其24邻域像素中已解缠像素的“影响力”.在待解缠像素的相位估算和检验上，充分考虑解缠像素的相位连续性和方向性，并对个别低相干像素进行反复迭代解缠.

由于解缠前对研究区域进行了区域识别，使得一些关键的研究区域或边界不被邻近的噪声等影响，例如，由桥梁或堤坝等狭长区域连接孤立的岛屿等场景.传统相位解缠方法从陆地沿着桥面向岛屿方向解缠，由于桥梁的宽度较窄，即使桥面的相位准确，解缠过程中也容易受其两侧低相干水体的误差影响.而采用本文方法，先排除水体区域，这样就能避免传统方法的不足.本文提出的区域识别和区域扩展解缠算法运用了多种遥感图像处理的方法，并对质量指导解缠方法进行了修正，提高了InSAR技术对地观测的反演精度.将该方法应用于PALSAR和ASAR的数据，并与其他几种典型方法作了比较，实验结果证明了本文方法在解缠精度上优于其他几种解缠方法.

本文方法也存在一些不足，区域识别和计算已解缠像素的“影响力”等操作，增加了计算时间和资源占用.另外，本文方法适用范围存在一定的限制，它依赖于可见光图像中不同区域清晰的光谱特征和实验数据的相干信息.

致谢：衷心感谢香港中文大学林珲教授、陈富龙和江利明研究员，在本人博士后期间给予的指导和有益讨论，并提供了航空摄影实测的DEM数据，使本文算法得以验证.

［1］ZEBKER H A，GOLDSTEIN R M.Topographic mapping from interferometric SAR observation［J］.J Geophys Res，1986，91（B5）：4993-5001.

［2］GHIGLIA D C，PRITT M D.Two-Dimensional Phase Unwrapping：Theory，Algorithms，and Software［M］.New York：John Weley ＆Sons Inc，1998.

［3］GOLDSTEIN R M，ZEBKER H A，WERNER C L.Satellite radar interferometry：two-dimensional phase unwrapping［J］.Radio Science，1988，23（4）：713-720.

［4］XU W，CUMMING I.A region-growing algorithm for InSAR phase unwrapping［J］.IEEE Trans Geosci Remote Sensing，1999，37（1）：124-134.

［5］FORNARO G，SANSOSTI E.A two-dimensional region growing least squares phase unwrapping algorithm for interferometric SAR processing［J］.IEEE Trans Geosci Remote Sensing，1999，37（5）：2215-2226.

［6］FLYNN T J.Two-dimensional phase unwrapping with minimum weighted discontinuity［J］.Journal of the Optical Society of America A，1997，14（10）：2692-2701.

［7］GHIGLIA D C，ROMERO L A.Minimum Lp-norm two dimensional phase unwrapping［J］.Journal of the Optical Society of America A，1996，13（10）：1999-2013.

［8］WEI Zhiqiang，XU Feng，JIN Yaqiu.Phase unwrapping for SAR interf erometry based on ant colony optimization algorithm［J］.International Journal of Remote Sensing，2008，29（3）：711-725.

［9］WEI Zhiqiang，JIN Yaqiu.The multiple patches and center expansion phase unwrapping algorithms for In-SAR image［J］.International Journal of Remote Sensing，2010，31（10）：2757-2765.

［10］WEI Zhiqiang，LIN Hui.Two InSAR procedures for wide and noise interferograms［C］／／IEEE International Conference on Audio，Language and Image Processing.Shanghai，China，2010：1552-1556.

［11］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.

［12］张 妍，冯大政，曲晓宁.基于改进粒子群算法的二维相位解缠方法［J］.电波科学学报，2012，27（6）：1116-1123.ZHANG Yan，FENG Dazheng，QU Xiaoning.Twodimensional phase unwrapping method based on improved particle swarm optimization［J］.Chinese Journal of Radio Science，2012，27（6）：1116-1123.（in Chinese）

［13］黄海风，王青松，张永胜.一种多基线相位解缠频域快速算法［J］.电波科学学报，2011，26（5）：831-836.HUANG Haifeng，WANG Qingsong，ZHANG Yongsheng.A fast phase unwrapping algorithm for multibaseline interferometry in frequency domain［J］.Chinese Journal of Radio Science，2011，26（5）：831-836.（in Chinese）［14］毛志杰，廖桂生，李 海.自适应图像准多基线In-SAR干涉相位展开方法［J］.电波科学学报，2008，

23（5）：877-882.MAO Zhijie，LIAO Guisheng，LI Hai.Multibaseline InSAR with image auto-coregistration for phase unwrapping［J］.Chinese Journal of Radio Science，2008，23（5）：877-882.（in Chinese）

［15］RICHARDS J A，JIA X.Remote Sensing Digital Image Analysis［M］.4th ed.Springer，2006.