改进的多基线相位解缠绕CANOPUS算法

刘会涛，张 欢，邢孟道，保 铮

（1.西安电子科技大学雷达信号处理国家重点实验室，陕西西安710071；2.西安外事学院商学院，陕西西安710077）

改进的多基线相位解缠绕CANOPUS算法

刘会涛1，张 欢2，邢孟道1，保 铮1

（1.西安电子科技大学雷达信号处理国家重点实验室，陕西西安710071；2.西安外事学院商学院，陕西西安710077）

中国余数定理（Chinese remainder theorem，CRT）方法是多基线相位解缠绕技术的一个基本方法，但是其较差的噪声鲁棒性问题限制了该方法在多基线相位解缠绕中的应用，然而基于聚类分析的多基线相位解缠绕技术能够克服传统的CRT算法噪声鲁棒性差的问题，本文根据CANOPUS算法中的聚类方法，提出用L∞-norm的距离测度定义两点之间的距离，从而减少特小类的产生和降低噪点数，进而提高聚类分析的精度，并且改进CANOPUS算法的算法流程以提高聚类分析算法的执行效率，进而大幅度降低聚类分析的运算时间。通过用仿真数据和实测数据验证可得，本文提出的改进聚类方法的聚类分析精度和算法执行效率更高，有效性通过实测数据实验得到了验证。

干涉合成孔径雷达；多基线；聚类分析；相位解缠绕

0 引 言

干涉合成孔径雷达（interferometric synthetic aperture radar，InSAR）通过利用两幅SAR图像的绝对相位差可以得到地表的高程信息。但是传统的单基线InSAR技术存在无法有效解决复杂地形相位解缠绕问题和高程层叠效应问题［1－3］。然而通过增加不同频率或者基线的雷达以不同的视角对成像场景进行观测可以有效地克服单基线InSAR技术的技术瓶颈。多基线InSAR技术就是利用多个基线或者频率雷达对同一场景进行观测来实现复杂地形的三维测绘的。与单基线InSAR技术相比，多基线InSAR技术不仅可以有效地解决复杂地形的相位解缠绕问题和高程层叠效应问题，而且还可以提高模糊高度，进而有效地解决相位欠采样问题［4］。然而多基线InSAR技术仍然存在相位解缠绕问题，而且与单基线相位解缠绕问题相比，多基线相位解缠绕问题还面临着数据量大但可用信息较少的问题。

不同于单基线相位解缠绕方法，多基线相位解缠绕算法通过利用不同基线的干涉相位图的多样性实现相位解缠绕。理论上，只要基线长度满足互质条件，缠绕相位可以根据中国余数定理（Chinese remainder theorem，CRT）的理论实现解缠绕［5－6］。但是由于CRT较差的噪声鲁棒性问题，通过传统的CRT方法不能实现相位解缠绕。因此，要使传统的CRT算法能够应用于多基线相位解缠绕问题就必须增加CRT算法的噪声鲁棒性。文献［7］通过深入研究CRT算法，提出了具有闭式解的鲁棒的CRT算法。除此之外，还有一些算法通过研究多基线相位解缠绕所面临的实际问题，将多基线相位解缠绕所独有的特殊性质辅助传统的CRT算法，以提高CRT算法的噪声鲁棒性。文献［4］通过研究多基线数据所特有的聚类特性，将具有相同模糊数向量的像素聚为一个类，然后利用类内所有像素降低相位解缠绕所需数据的噪声并将类内所有像素作为一个整体进行解缠绕聚类分析（cluster-analysis，CA）。再者，将CA算法和具有闭式解的鲁棒的CRT算法相结合，从而在两方面提高传统CRT算法的噪声鲁棒性［8］。然而，由于CA算法中的聚类算法本身噪声鲁棒性不强，导致具有相同模糊数向量的像素未必能够同时包含在一类内，为了改善CA算法的噪声鲁棒性，文献［9］使用多基线数据的密度连通特性提高聚类的精确度，进而提高多基线相位解缠绕算法的噪声鲁棒性（cluster-analysis based noise robust phase-unwrapping algorithm，CANOPUS），同时该算法也是目前国际学术界提出的唯一一个能够快速处理大规模多基线相位解缠绕的方法。除了基于传统的CRT算法的多基线相位解缠绕算法之外，文献［10］采用线性规划的思想解决多基线相位解缠绕问题，提出了利用L1范数的多基线相位解缠绕算法。文献［11］提出了用最大后验估计（maximum a posterior estimation，MAP）的方法实现相位解缠绕，并且该算法在文献［12］中得到了改进和发展。文献［13－14］提出了利用最大似然估计（maximum likelihood，ML）的方法实现多基线相位解缠绕方法，文献［15］利用CA算法中的聚类分析提出了一种改进的基于ML方法。文献［16］发展和改进了MAP和ML估计方法，并给出了ML方法的克拉美罗界和MAP方法的误差下界。文献［17］提出用幅度信息辅助实现多基线相位解缠绕的方法。

本文首先回顾了CANOPUS算法，然后提出了改进的CANOPUS聚类算法，以使聚类分析运行速度更快，聚类结果更有利于多基线相位解缠绕，最后通过仿真数据和实测数据实验验证本文提出的改进的聚类方法的有效性。

1 CANOPUS算法

根据多基线InSAR中相位和基线长度的约束关系，可以得到模糊数、缠绕相位和基线长度之间的关系为（以双基线InSAR系统为例）

式中，k、φ和B分别表示模糊数、缠绕相位和有效基线长度，其下标用于区别不同的干涉相位图。与传统的CRT算法相比，CA算法［4］利用具有相同模糊数向量的像素其对应的式（1）中的截距应该相同这一特性，将具有相同模糊数向量的像素聚为一类，并将类内所有像素视为一个整体对截距进行降噪处理，然后再进行相位解缠绕。然而，由于现实因素的影响，具有相同模糊数向量像素的截距不再一成不变，该截距值的变化将导致具有相同模糊数的像素不能被聚为一类和具有不同模糊数向量的像素可能被聚为一类两方面的劣势。为了克服截距空变的影响，CANOPUS算法［9］提出用截距数据（下文称之为待聚类数据）的密度信息来区分不同的类。为了更好地理解本文提出的改进策略，首先详细介绍CANOPUS算法。

与CA算法相比，CANOPUS算法将更多的待聚类数据中的信息用于聚类分析，不仅利用了截距信息，而且还利用了截距在整个待聚类数据中的相对位置信息。以双基线数据为例，CANOPUS算法将三维的待聚类数据中的一点视为聚类模式，包含三方面信息，即点在待聚类数据中的行和列信息以及其截距值。如此选择聚类模式的好处是既能通过截距值识别具有相同模糊数的类又能克服截距空变带来的影响。假设三维待聚类数据中的一点p，定义以点p为球心ε为半径的球体为点p的ε邻域，并定义点p的ε邻域内包含的待聚类数据中的点的数量为点p的密度Nε（p）。根据待聚类数据中截距值变化的特点，即相邻类之间截距值变化较大（多基线数据截距值固有的属性）而同一类内截距值缓慢变化（受现实因素影响导致的截距值的空变现象），则处于类边界的点的密度理应比类内点的密度小，因此只要选择合适的判断准则就能区分出哪些点是边界点（密度小于最小点数的点即为边界点）。图1为重复航过的多基线InSAR系统得到的待聚类数据及其密度图。从图1（a）可以看出，由于现实因素影响导致的截距的空变特性相当严重，除此之外，不同类之间的某些像素存在截距值的相互干扰，从而导致CA算法失效。然而，仅仅通过肉眼，不同的类还是很容易区分开来的，这是因为相邻类之间的边界点存在着明显的截距值的差。图1（b）为图1（a）对应的密度图，从图1（b）可以看出，相邻类之间的边界点处其密度明显低于类内其他像素。

图1 实测数据的待聚类数据及其密度图

从图1可以看出，不同的类被低密度像素区分开来，因此只要选择合适的最小点数即可得到类的边界点。然而如何确定哪些像素属于同一类呢？CANOPUS算法提供了一种基于密度连通的聚类方法，该方法是由模式识别中经典的聚类算法（density-based spatial clustering of applications with noise，DBSCAN）改进而来［11］。在阐述该方法之前需首先定义密度大于最小点数MinPts的待聚类数据中的点为核心点，相反，称其为非核心点。如果某点p到核心点q的距离小于ε，则称点p是从点q直接密度可达的。如果存在一系列的点p1，p2，…，pn，p1＝p，pn＝q，并且点pi是从点pi＋1（i＝1，2，…，n－1）直接密度可达的，则称点p是从点q密度可达的。如果存在点o使得点p和pi＋1（i＝1，2，…，n－1）点q是从点o密度可达的，则称点p和点q是密度连通的。CANOPUS所提供的聚类方法是将被边界点分开的区域内所有的密度连通的点视为一类。其算法流程如下［9，18］：

将待聚类数据定义为Xun；

①设置分类编号m＝1；

②WHILE Xun≠Φ

③ 选择Xun中任意一点p；

④ IF p是非核心点THEN

⑤ 将p赋予噪点集合N；

⑥ Xun＝Xun－｛p｝；

⑦ ELSE IF p是核心点THEN

⑧ 找到所有与点p是密度可达的点p

⑨ 并将其赋予类Cm；

⑩ Xun＝Xun－Cm；

从上述CANOPUS中的聚类算法流程可以看出，聚类分析开始于类内的任何一个核心点，终止于非核心点，而非核心点是由最小点数MinPts唯一决定的。因此只要选择的参数合适，即可实现对待聚类数据的聚类分析，并且该聚类结果仅与选择的参数值ε和MinPts有关。除此之外，该聚类算法的时间复杂度仅为O（n），其中n为待聚类数据的大小。

2 改进的CANOPUS聚类方法

由于CANOPUS算法中的CA方法并未考虑多基线InSAR数据的特点，因此本文从两个方面改进CANOPUS算法。首先，CANOPUS算法中的聚类算法在计算任意两点之间的距离时采用的是L2-norm距离测度，如此选择的好处是其几何意义更直观，但是存在算法执行时间长和聚类精度不理想的问题。因此本文提出用L∞-norm距离测度定义两点之间的距离。选择L∞-norm距离测度的好处主要表现在两方面：一是易于判断任何两点之间是否直接密度可达；二是对距离越远的点其容许的截距变化相对而言比L2-norm距离测度更大，因此聚类过程产生的噪点的数量会降低。其次，在不改变算法实现功能的条件下，改变算法运行的流程以提高算法的执行速度。改进后的算法流程如下：

① 将待聚类数据定义为Xun；

② FORi＝1：n

③ IF Xun（i）是核心点

④ Xun（i）赋值核心点标志位；

⑤ END IF

⑥ END FOR

⑦ 设置分类编号m＝1；

⑧ WHILE Xun≠Φ

⑨ 选择Xun中任意一核心点p；

⑩ 找到所有与点p是密度可达的点并将其赋予类Cm；

与传统CANOPUS算法流程相比，本文中提出的改进算法首先识别并标识核心点与非核心点，然后再进行聚类分析。其优势表现在两个方面，首先，通过优先确认核心点后避免了原CANOPUS算法第4行执行非核心点的可能，从而提高算法执行效率。更为关键的是在传统的CANOPUS算法中执行第9行确定与核心点p密度可达的点时，需要首先确认与核心点p直接密度可达的点并确认哪些点为核心点，然后再以所确认的核心点为源查找所有与这些核心点直接密度可达的点，以此类推，以至于再也没有其他点与此类内的核心点直接密度可达。再在传统的CANOPUS算法流程中，没有提前确定并表示核心点的前提下，执行第9行需要重复计算核心点邻域内的点是否为核心点，存在重复计算的可能，因此其算法运算效率较低。然而在本文中提出的改进的CANOPUS算法中，核心点的确认与标示被提前完成，因此在执行本文中算法第10行时，不存在重复计算某些点是否为核心点的运算，因此本文提出的对传统CANOPUS算法的改进能够极大地提高算法的执行效率。

3 算法的性能分析

本文提出的改进的CANOPUS算法将通过大小为458像素×157像素的Isolation Peak国家公园仿真数据和大小为7 697像素×8 000像素的TerraSAR-X实测数据验证。Isolation Peak国家公园仿真数据的干涉相位图如图2（a）和图2（b）所示，其基线长度分别为345.27m和281.46m。本文中仿真数据的相干系数约为0.8，因此在图2（a）和图2（b）中的干涉相位图的噪声较大。图2（c）是待聚类数据，图2（d）是用传统的CANOPUS算法中的聚类算法得到的聚类结果，其包含74个特小类（类内聚类点数不超过30个），包含41 021个噪点数。图2（e）是用本文提出的改进的聚类算法得到的聚类结果，其包含62个特小类，包含37 430个噪点数。与传统CANOPUS算法中的聚类算法相比，本文提出的改进措施在减少特小类和减少噪点数两方面都有改进，而且需要指出的是噪声越严重，本文提出的改进方法对聚类算法性能的改善越大。图2（f）是用本文提出的改进的聚类算法实现的相位解缠绕结果。

图2 仿真数据的多基线相位解缠绕实验结果

TerraSAR-X实测数据的主要参数如表1所示。本文实验采用的实测数据为单颗卫星多次对同一地区进行重复观测获得的，即重复航过多基线数据。通过其主要参数可以发现，本次实验所用数据存在较严重的时间去相干，因此数据受噪声影响较严重。相位解缠绕结果如图3所示，图3（a）为短基线的干涉相位图，图3（b）为长基线的干涉相位图，图3（c）为利用长短基线干涉图和基线数据根据式（2）得到的截距图。通过用本文提出的改进聚类方法对图3（c）进行聚类分析得到的聚类结果如图3（d）所示，与其对应的多基线相位解缠绕结果如图3（e）所示，从图3（e）中可以看出，解缠绕结果中仍然存在大量的噪点，而导致噪点产生的主要原因是河流湖泊等相干性差的区域的存在。为了能够清晰地对比本文中改进的聚类算法与CANOPUS算法中的聚类算法的聚类结果的差别，本文将传统的CANOPUS算法中的聚类分析结果呈现在图3（f）中。为了减小图像大小，图2中所有的图像均经过16倍降采样。

表1 实验数据主要参数

本文从3个方面比较了CANOPUS中的聚类算法以及改进的聚类算法的聚类效果。首先从聚类分析结果中包含的聚类数来看，如图3（d）和图3（f）所示，本文中提出的改进算法聚类数降低了约69.1%，说明采用改进的聚类方法后能够有效降低特小类的产生，而特小类的产生实际上是因为聚类算法聚类的过程中将原本属于一类的点分成多个类，如图3（d）和图3（f）中白色矩形框内所示聚类结果。特小类的产生不会导致多基线解缠绕CRT算法失败，但是会降低基于聚类分析的CRT算法的噪声鲁棒性（这是因为同一类内点数越少，根据类内点实现的滤波效果越差）。其次，从聚类结果中的噪点数来看，改进后的聚类算法产生的噪点数降低了约13.7%，说明不可解缠绕点减少了，解缠绕后的相位图像更加精细。从聚类结果中的聚类数和噪点数两方面来看，本文中提出的采用L∞-norm距离测度代替L2-norm距离测度的方法是行之有效的。最后从聚类分析的执行时间来看，本文提出的改进后的聚类方法其运算时间比原聚类算法降低了24.8%（CPU：Intel（R）Core（TM）i5 3.20GHz，RAM：6GB），说明本文中对聚类算法的流程的改进能够极大促进算法的执行效率。表2给出了上述3点聚类算法性能分析的对比结果。

4 结束语

与单基线相位解缠绕技术发展的程度相比，多基线相位解缠绕技术仍然面临较多难题。为了克服传统CRT算法噪声鲁棒性差的问题，相位解缠绕的科技工作者们提出了用聚类分析的方法来增强传统CRT算法的噪声鲁棒性，本文提出的聚类分析方法是对CANOPUS算法中的聚类分析方法的改进，以使其聚类分析结果更加精细，运算时间更短。通过实测数据验证可以发现，本文提出的改进的聚类算法能够减少特小类的产生，降低噪点数和降低聚类分析时的运算时间。

［1］Yu H W，Xing M D，Bao Z.A fast phase unwrapping method for large-scale interferograms［J］.IEEE Trans.on Geoscience and Remote Sensing，2013，51（7）：4240－4248.

［2］Zhang Y，Feng D Z，Qu X N.Hybrid phase unwrapping algorithm combining branch-cut and surface-fitting for InSAR［J］.Journal of Xidian University，2012，39（5）：47－53.（张妍，冯大政，曲小宁.支切法与曲面拟合结合的InSAR相位展开算法［J］.西安电子科技大学学报，2012，39（5）：47－53.）

［3］Yu H W，Li Z F，Bao Z.Residues cluster-based segmentation and outlier-detection method for large-scale phase unwrapping［J］.IEEE Trans.on Image Processing，2011，20（10）：2865－2875.

［4］Yu H W，Li Z F，Bao Z.A cluster-analysis-based efficient multibaseline phase-unwrapping algorithm［J］.IEEE Trans.on Geoscience and Remote Sensing，2011，49（1）：478－487.

［5］Xu W，Chang E C，Kwoh L K，et al.Phase-unwrapping of SAR interferogram with multi-frequency or multi-baseline［C］∥Proc.of the IEEE International Geoscience and Remote Sensing Symposium，1994：730－732.

［6］Jin G W，Zhang H M，Xu Q，et al.Phase unwrapping algorithm with CRT for multi-band InSAR［J］.Journal of Xidian University，2011，38（6）：97－102.（靳国旺，张红敏，徐青，等.多波段InSAR的CRT相位解缠绕方法［J］.西安电子科技大学学报，2011，38（6）：97－102.）

［7］Wang W J，Xia X G.A closed-form robust chinese remainder theorem and its performance analysis［J］.IEEE Trans.on Signal Processing，2010，58（11）：5655－5665.

［8］Yuan Z H，Deng Y K，Li F，et al.Multichannel InSAR DEM reconstruction through improved closed-form robust chinese remainder theorem［J］.IEEE Geoscience and Remote Sensing Letters，2013，10（6）：1314－1318.

［9］Liu H T，Xing M D，Bao Z.A cluster-analysis based noise robust phase-unwrapping algorithm for multi-baseline interferograms［J］.IEEE Trans.on Geoscience and Remote Sensing，2015，53（1）：494－504.

［10］Yu H W，Bao Z.L1-norm method for multi-baseline InSAR phase unwrapping［J］.Journal of Xidian University，2013，40（4）：37－41.（于瀚雯，保铮.利用L1范数的多基线InSAR相位解缠绕技术［J］.西安电子科技大学学报，2013，40（4）：37－41.）

［11］Ferraiuolo G，Pascazio V，Schirinzi G.Maximum a posteriori estimation of height profiles in InSAR imaging［J］.IEEE Geoscience and Remote Sensing Letters，2004，1（2）：66－70.

［12］Ferraioli G，Shabou A，Tupin F，et al.Multichannel phase unwrapping with graph cuts［J］.IEEE Geoscience and Remote Sensing Letters，2009，6（3）：66－70.

［13］Pascazio V，Schirinzi G.Multifrequency InSAR height reconstruction through maximum likelihood estimation of local planes parameters［J］.IEEE Trans.on Image Processing，2002，11（12）：1478－1489.

［14］Fornaro G，Pauciullo A，Sansosti E.Phase difference-based multi-channel phase unwrapping［J］.IEEE Trans.on Image Processing，2005，14（7）：960－972.

［15］Yuan Z H，Deng Y K，Li F，et al.Improved multichannel InSAR height reconstruction method based on maximum likelihood estimation［J］.Journal of Electronics ＆Information Technology，2013，35（9）：2161－2167.（袁志辉，邓云凯，李飞，等.改进的基于最大似然估计的多通道InSAR高程重建方法［J］.电子与信息学报，2013，35（9）：2161－2167.）

［16］Ferraiuolo G，Meglio F，Pascazio V，et al.DEM reconstruction accuracy in multichannel SAR interferometry［J］.IEEE Trans.on Geoscience and Remote Sensing，2009，47（1）：191－201.

［17］Shabou A，Baselice F，Ferraioli G.Urban digital elevation model reconstruction using very high resolution multichannel InSAR data［J］.IEEE Trans.on Geoscience and Remote Sensing，2012，50（11）：4748－4758.

［18］Theodoridis S，Koutroumbas K.Pattern recognition［M］.4th ed.Burlington：Academic Press，2009：764－862.

Improved CANOPUS method for multi-baseline interferograms

LIU Hui-tao1，ZHANG Huan2，XING Meng-dao1，BAO Zheng1
（1.National Lab of Radar Signal Processing，Xidian University，Xi’an 710071，China；2.College of Bussiness，Xi’an International University，Xi’an 710077，China）

As a basic technique against the multi-baseline phase unwrapping problem，Chinese remainder theorem（CRT）restricts itself in the application of multi-baseline phase unwrapping due to its worse noise robustness.However，multi-baseline phase unwrapping algorithms based on the cluster-analysis are able to overcome the drawbacks of the traditional CRT method.Based on the cluster-analysis method in the CANOPUS algorithm，an L∞-norm is employed to define the distance between two elements to decrease the production of very small clusters and the number of noise points so as to improve the cluster-analysis（CA）performance.Besides，the procedure of the CA method in the CANOPUS algorithm is changed to improve the efficiency of the method and decrease the consuming time of the algorithm.According to the experiments on a simulated and repeat-pass real interferometric synthetic aperture radar（InSAR）dataset，the effectiveness and efficiency of the improved cluster-analysis method are tested.

interferometric synthetic aperture radar（InSAR）；multi-baseline；cluster-analysis（CA）；phase unwrapping

TN 957.52

A

10.3969／j.issn.1001-506X.2015.08.08

刘会涛（1986－），男，博士研究生，主要研究方向为InSAR相位解缠绕处理、多基线相位解缠绕。

E-mail：huitaoliuxd＠gmail.com

张 欢（1984－），女，助教，硕士，主要研究方向为应用统计学。

E-mail：huanzhangxd＠126.com

邢孟道（1975－），男，教授，博士，主要研究方向为雷达成像和目标识别。

E-mail：xmd＠xidian.edu.cn

保 铮（1927－），男，院士，主要研究方向为雷达信号处理。

E-mail：piaofei8＠gmail.com

1001－506X201508－1767－06

网址：www.sys－ele.com

2014－10－27；

2015－01－26；网络优先出版日期：2015－03－30。

网络优先出版地址：http：／／www.cnki.net／kcms／detail／11.2422.TN.20150330.1111.015.html

国家自然科学基金优秀青年基金项目（61222108）资助课题