相位解缠算法研究.doc

一、引言合成孔径雷达干涉测量技术（synthetic aperture radar interferometry, InASR）将合成孔径雷达成像技术与干涉测量技术成功地进行了结合，利用传感器高度、雷达波长、波束视向及天线基线距之间的几何关系，可以精确的测量出图像上每一点的三维位置和变化信息。 合成孔径雷达干涉测量技术是正在发展中的极具潜力的微波遥感新技术，其诞生至今已近30年。起初它主要应用于生成数字高程模型(DEM)和制图，后来很快被扩展为差分干涉技术 ( differential InSAR , DInSAR)并应用于测量微小的地表形变，它已在研究地震形变、火山运动、冰川漂移、城市沉降以及山体滑坡等方面表现出极好的前景。特别，DInSAR具有高形变敏感度、高空间分辨率、几乎不受云雨天气制约和空中遥感等突出的技术优势，它是基于面观测的空间大地测量新技术，可补充已有的基于点观测的低空间分辨率大地测量技术如全球定位系统(GPS)、甚长基线干涉 (VLBI)和精密水准等。尤其InSAR在地球动力学方面的研究最令人瞩目。 二维相位解缠是 InSAR 数据处理流程中重要步骤之一，也是主要误差来源，无论是获取数字高程模型还是获取地表形变信息，其精确程度都高度依赖于有效的相位解缠。因此，本人在课程期间对相位解缠的相关文献进行了阅读。二、InSAR基本原理用两副雷达天线代替两个光源， ，对地面发射相干信号，将得到类似的条纹图。因为雷达信号与光线本质上都是电磁波，所以只要保证雷达天线载具运行轨道的稳定，那么两个信号到达地面上某一点处的路程差是确定的，只与该点在地面上的位置有关。在 InSAR 干涉测量中有两种模式，一种是在载具（卫星或飞机）上搭载一具天线，而载具两次通过不同轨道航线飞经目标地域上空，此种称之为单天线双航过模式；另一种在载具上搭载两副天线，只飞经目标地域上空一次，此种方式称之为双天线单航过模式。不论是哪种方式都可以用图 2.2 来模拟并作出几何解释。在测量中两副天线或两次航过接收的数据可以各获得对地面同一区域的两幅包含幅值与相位信息的二维复数据图像，分别以， 表示为 （2.6）其中和表示幅值信息，和表示相位信息。将两幅图像作共轭乘，可得 （2.7）为两幅图像中相对应的像点的相位差，由路程差决定的，由余弦定理有 （2.8）可得 （2.9）根据式（2.7）的结论，两路雷达波路程差与相位差成正比 （2.10）式（2.8）可以进一步得到 （2.11）于是 （2.12）上式中 B 为基线长，由此可以获得地面的高程信息。这里关键是利用了路程差与相位差成正比这样一个关系，应该注意的是两天线接收到的信号的路程差并不很大，但是由于高频的雷达信号的波长 很小，所以可以很大，即两个信号的相位差可以比4大很多。但是由式（2.7）计算相位差时会以2为模来取值，得到的相位只会在 ( ,之间，称为相位的主值或缠绕相位，它与真实相位的关系是相差 2 的整数倍，即有下式的关系 k=0，1，2 （2.13）根据缠绕相位得到真实相位的处理过程就叫做相位解缠，是 InSAR 干涉测量的关键步骤。三、相位解缠基本原理3.1引言在上节提到利用相位差能获得精确的路程差进而获得地面的高程信息，因此获得准确的相位差就是实现测量的关键。由于复数对其相位的周期性，InSAR根据两幅 SAR 复图像获得的干涉相位差值是被周期折叠后位于 ( ,之间的相位主值，它与真实的相位差值之间存在着 2 k差别。由式(2.13)可以表示它们之间的基本关系。其中 代表解缠相位， 代表缠绕相位。必须对 进行相位解缠，恢复被模糊掉的相位周期，获得目标在两次成像中的真实相位差，才能得到目标的正确高度信息。相位解缠是 InSAR 三维成像处理中的关键步骤之一，其准确程度将直接决定数字高程图（DEM）和地表形变探测的精度。3.2相位缠绕和解缠理想情况下，图像的采样率满足 Nyquist 采样定理，采样频率必须大于信号最高频率的两倍，解缠绕的干涉相位中相邻像素点之间的相位差值不可能超过半个周期（一个）。当满足此条件时必然能由缠绕相位解缠出正确的解缠绕相位，并且可以通过积分进行解缠。记 (m)为周期缠绕前的真实相位值， (m)为相应的缠绕相位，定义相位缠绕算子 ，相位缠绕的过程可以用式（3.1）表示 （3.1）结果是得到主值属于 ( ,区间的缠绕相位。定义差分算子 ，根据 Nyquist采样定理对于解缠相位有 （3.2）对相邻缠绕相位进行差分运算得 （3.3）对该相位差也使用缠绕算子得 （3.4）根据缠绕算子的定义，其结果必须属于 ( ,区间，而 (m)也必须属于( ,区间，所以有 （3.5）式（3.3）变为 （3.6）由式（3.2）可得 （3.7）由式（3.7）可以看出，通过对相邻缠绕相位之差积分可以实现相位解缠，条件是满足 Nyquist 采样定理。对于一维的情况，可以简单的使用如下的公式进行解缠计算，记 (m)为周期缠绕前的真实相位值， (m)为相应的缠绕相位，计算干涉图中一个点到下一个点的相位变化，即计算相位梯度，然后从一固定点开始积分使相位值的变化平稳连续，从而恢复失去的相位周期。即下式： （3.8）若有如下的一维相位序列0.2 , 0.5 , 0.6 , 0.8 , -0.5 , -0.4 , -0.2以相邻的 0.8 , -0.5 两个数据为例， 0 .5 (0.8)= 1.3，因为 1. 3 所以 ( m )= 1.3 +2=0.7，将 0. 7加上前一个解缠结果 0. 8得到该位置的解缠结果为 1 .5。其他照此进行，从左向右解缠后的序列为：0.2 , 0.5 , 0.6 , 0.8 , 1.5 , 1.6 , 1.8 。由于一维序列的积分路径是唯一的，所以其解也是唯一的。但由于是逐个积分，如果受到相位噪声的影响，或者碰到地形起伏本来就不满足相邻缠绕相位差的绝对值小于 的条件，使其中一点的解缠绕相位发生错误，则错误会后向传播，导致之后所有相位的解缠结果与真实相位相差甚远。为了说明相位缠绕与解缠原理，选取如图 3.1 所示的人工模拟的简单缠绕相位图进行解释。在理想状况下，发生缠绕的干涉相位呈现周期性变化，由 渐变到 ，然后由 突变为 ，如此反复，从图像上表现为灰度值由浅渐渐变深，然后突变为浅色，再向深色渐变，形成如图 3.1（a）所示的条纹图。从图 3.1（a）中沿 y 轴方向取一条一维数据，以像素位置为横坐标，以灰度强弱代表的相位值为纵坐标将其表示出来将如图 3.1 （c）所示，其形状如锯齿状，代表了图 3.1（a）中黑白交替变换的条纹。理想情况下的解缠绕只需进行简单的积分将突变消除，整幅图像的条纹变成了连续的面，相位恢复连续变化。如图 3.1 (b）所示。在图 3.1 (b）中也取一条一维数据在坐标图中画出，将如图 3.1 (d）所示。四、常用相位解缠算法4.1常用相位解缠算法概述到目前为止，针对相位解缠问题已经提出很多解决方案。主要的解缠算法大致可以分为三类：一类可以称之为路径跟踪解缠算法，他们的共同特点是采用路径积分来实现相位解缠，以1988年Goldstein提出的枝切法（Branch-Cut）为代表。枝切法通过探测残差点，用枝切线连接残差点，然后进行路径积分来实现解缠，在路径积分时以不穿越枝切线为原则。Wei Xu和Cumming提出的区域生长法（Region-Growing）不考虑残差点，不布置枝切线，而是依据额外信息将干涉图划分为高质量低质量区域，在各个区域内按照从高质量像元到低质量像元的方向进行路径积分 。 Flynn的掩模分割法（Mask- Cut）和最小不连续法（Minimum-discontinuity）等也属于该类算法。另一类算法着眼于整体，采用最优化的思想，寻求最小二乘意义下的最优解缠结果，包括用FFT/DCT方法求解的无加权最小二乘算法，Pritt的多重网格迭代法求解加权最小二乘相位解缠法，Ghiglia的最小范数法等。这类算法不探测残差点，不布置枝切线，通过建立一个离散型泊松目标函数，并用各种数学的方法求解它以实现相位解缠。第三类方法为最小费用流方法，以Costantini的基于网络规划的解缠方法为代表，引入图论中的网络模型，将解缠问题转变为解一个网络最小费用流的问题，利用网络规划理论中成熟高效的算法求解。4.2基于路径的相位解缠算法两幅SAR图像经过干涉以后，我们可以获得一幅缠绕相位图像，各像元上的值为对应的干涉相位的主值。根据Nyquist定理，当相邻像元上的相位差小于二时，可以通过积分的算法来恢复相位的真实值。基于路径跟踪的相位解缠算法就是通过积分相邻缠绕相位的差分值来恢复相位的真实值的。假设我们己知在像元上的相位，那么在其它像元r上的相位可以通过以下公式来获得: （4.1）符号为像元r上的解缠相位，为像元上的已知解缠相位，C为积分路径，根据积分理论： （4.2）上式为积分函数，C为积分路径，（4.2）的线性积分不仅依赖于积分路径C的起点和终点，还依赖于积分路径C本身。要使积分与路径无关，则要求一下闭合积分成立 （4.3）在二维相位解缠中，公式(4.3)常用来作为探测积分是否与路径无关的条件。InSAR缠绕相位数据中，不是所有的积分路径都满足公式(4.3)，有些像元上的缠绕相位数据由于受到噪声的影响，或者由于其它的原因，导致通过这些像元的闭合积分不能满足公式(4.3)，这些像元上的相位在InSAR中被称为“残差点(residue)”，或者“电荷”(具有正负性，见随后的讨论)，在路径跟踪的相位解缠算法中，关键的问题在于如何判断这些电荷并将它们相连(称为“分枝”)以达到正负抵消，且防止积分路径穿过这些分枝。路径跟踪法的基本策略是将可能的误差传递限制在噪声区内，通过选择合适的积分路径，隔绝噪声区，阻止相位误差的全程传递。几十年来，研究者研究出了许多的相位解缠算法，至今为止，基于路径跟踪的相位解缠算法有枝切法、区域法、Mask-cut算法、像元扩散法、最小生成树法、条纹检测法、区域生长法，最小不连续算法 (简称Rynn算法)等算法。路径跟踪的相位解缠算法一般步骤如下所示:输入:缠绕相位步骤1:相位连续性/不连续性检测:识别残差点，生成枝切线。步骤2:计算/建立相位质量图。步骤3:相位积分:在枝切线周围或在质量图的指导下处理。下面对其中比较典型的算法作详细介绍:4.2.1Goldstein枝切法Goldstein枝切法是较经典的路径跟踪法，是1988年Goldstein等人提出的，它识别正负残差点，并连接邻近的残差点对或多个残差点，实现残差点“电荷”平衡，生成最优的枝切线，确定积分路径，防止误差沿积分路径传递。基本步骤为:(l)识别残差点;(2)生成枝切线;(3)绕过枝切线进行积分。具体如下:首先按一定的顺序寻找残差点，定义一个2x2像元的缠绕相位为节点，将四个像元串接起来，即为影像中的最小闭合路径 (dosedloop)。沿最小闭合路径将缠绕相位梯度累加起来，如果之和为零则这四个点是一致的，否则左上角的像元就称为残差点(residual)。缠绕相位节点图与最小闭合路径图如图所示:缠绕相位节点，最小闭合路径图(每个像元数值乘以劫才表示缠绕相位真实值)，为各方向相位差:计算得 说明相位是一致的。再给出一个例子：缠绕相位节点图与最小闭合路径图如图所示: 则左上角的像元为残差点，当找到第一个残差点以后，从该残差点开始，继续搜索，找到下一个残差点后，用一枝切线将两者连接起来计算残数和，如果和为零则完成了该树枝的生长，继续搜索直至搜索完全部残差点，如果和不为零则不断加入残差点每次计算总的残数和，直至和为零。在Goldstein的枝切法中，有两个步骤是极为关键的:(l)当搜索窗口找到新的残差点，无论该残差点是否与其他的残差点相连，都将该残差点与窗口中心的残差点相连;(2)当搜索窗口到达图像的边界，则将残差点与边界相连，以阻止积分路径。枝切法最大的优点是:在实际计算中速度比较快;在噪声比较低、残差点比较少的情况下，精确度非常高。缺点是:当残差点较多且分布较密集时，该算法难以准确连接枝切线，导致无法选择合理的积分路径，有时会造成错误的阮跳跃，导致误差的传递。但由于该算法的速度优势，使之成为一种常用的相位解缠算法。4.2.2质量引导法这种算法不识别残差点，也不设置枝切线。而是在进行相位解缠时，通过相位质量图 (quality map)来定义相位数据的质量，将积分路径总是沿“高”质量的像元进行，最后解缠“低”质量像元。质量引导法的关键步骤就是在相位质量图的引导下进行像元扩散，其基本操作过程如下:从高质量像元点出发，检测它的四个邻近像元，对邻近像元进行解缠，将解缠后的像元的邻接像元(未解缠)存储在“邻接列”中，依据相位质量从“邻接列”移出高质量像元进行相位解缠，更新“邻接列”，重复上述步骤直至所有的像元解缠完毕。质量引导法成功地进行相位解缠的前提是必须有可靠的相位质量图。相位质量图主要有四种:相干图、伪相干图、相位导数变化图、最大相位梯度图、掩模图(mask)。(1) 相干图 最常用的质量图是相干图，相干值的高低表明图像不同区域的相干性，是最直观的干涉质量评价图。同时相干系数的变化也表征了在图像获取期间地物的变化情况，所以相干图也用于地物的分类等。(2) 伪相干图(Pseudo-correlation) 当无法获得InSAR图像对的强度值时，常常用伪相干图来模拟相干图。这时把InSAR图像对的强度定义为1，那么伪相干定义为：k为视数。伪相干图的一个最大缺陷是：对于陡峭地形区，它标志为低质量数据区(即使这些相位数据质量很好，并且没有噪声)。这时就需要新的质量图，用来评价相位导数的统计变化特征。(3) 相位导数变化图 (Phase Derivative Variance) 相位导数变化定义如下: 相位导数变化不同于相干和伪相干。例如，在倾斜地表，如果相位变化率保持一定，则相位导数变化为0，而伪相干不为0。从严格意义上讲，相位导数变化表征的是相位数据的“坏”(badness)，而不是“好”(goodness)，但我们可以假定:如果相位导数变化是可以忽略的话，那么相位数据就是好数据。实验表明，在无法获得相干图的情况下，相位导数变化图是最可靠的相位质量图。(4) 最大相位梯度图 从相位图上可以看出，在噪声相位区往往相位梯度也很大，所以可以用最大相位梯度来表征相位数据的质量。一般最大梯度定义为： 最大相位梯度图也有伪相干的缺陷，在地形很陡峭(即相位变化显著，但无噪声)也表征为低质量数据。质量引导法完全依赖于质量图像来指导解缠路径的选择，因此在缺乏高质量的质量图像时，解缠效果将会很不理想。另外该算法不识别残差点，因此解缠路径就不可避免的可能会包围非平衡的残差点，这样就有可能会产生2k二的周期累加错误。(5) 掩膜图（mask） 在解缠过程中，有的区域失相关严重(如水面)导致相位不连续，有的区域地势平坦不需要进行滤波。在这种情况下，我们可以制作一个掩模图将这些区域掩盖起来，使解缠和滤波不涉及这些区域。掩模图的生成方法一般是采用一定的闲值来进行判断，或者手工确定需要掩盖的区域，通常被掩盖区域的干涉图像元用0来表示，未被掩盖的像元1来表示。第四章对西安数据的解缠过程中均采用了掩模图，即掩盖掉的像元不参与解缠，以防止误差的传播。 质量引导法解缠结果的好坏很大程度上依赖于质量图的“好”与“坏”，如果有好的质量图引导，能得到优于枝切法的解缠结果。因为不设置枝切线，积分路径就可能包围残差点导致错误的2k的累积错误。4.2.3掩膜切口法（mask-cut法）1966年Rynn首先对质量图像引导分割路径的算法作了详细的描述，在该方法中使用逐步增长的像素掩模来连接残差点，这种像素掩模称为掩模切口，与分割路径的作用相似。掩模切口算法与质量引导法有些相似，它在有些地方可以看作是后者的逆向算法，它并不是从高质量的区域开始解缠，而是从残差点开始，沿低质量的区域逐渐扩展像素掩模。掩模切口的扩展过程一直到它连接了等量的正残差点和负残差点或者到达图像边界时为止，在前一种情况下，残差点是平衡的，相位解缠与积分路径无关，而在后一种情况下，掩模将残差点隔离开来，可以保证不会有解缠路径可能包围这些残差点。该算法实际上是将Goldstein枝切法和基于质量图路径积分法的结合起来，掩模分害算法虽然也识别残差点并产生分割路径，但总的来讲，该算法还是在质量图像的指导下实现的，因此算法的结果在很大程度上也仍然依赖于高质量的质量图像。在残留点不仅仅分布在低质量相位区的情况下，mask-cut算法与质量指导的路径跟踪算法都无法正常工作，在这种情况下，枝切法更为有效。其解缠步骤为:(l)识别残差点；(2)生成mask-cut；(3)细化mask-cut；(4)沿mask-cut的路径积分。4.3最小范数法理想情况下解缠相位梯度等于缠绕相位梯度的假设，相位解缠可以看成一个优化问题。最小范数法将相位解缠问题转化为数学上的最小范数问题，目前使用最广泛的是最小二乘法，最小二乘法是一种广泛使用的优化方法，有无权和加权两种形式。无权形势下，相位解缠是求取一个平滑的解缠相位。就是求解一个Neumann件下的Poisson方程。可以通过离散余弦变换DCT、离散傅里叶变换FFT或无权多级格网法来有效的解决。由于干涉图上各像素点相关系数差别较大，存在相位的不连续，因此无权最小二乘虽然获得了平滑的相位解缠曲面，但造成局部的噪声在最小均方意义下的全局传播而产生与真实相位值偏差较大的解。加权最小二乘法可以在一定程度上弥补无权最小二乘法的这一缺陷。例如共轭梯度法使用的权系数是经过二值化的质量图，将干涉图中由于残余点的存在而破坏的区域赋予零权，阻止它们对相位解缠的破坏。不带权的最小二乘法相位解缠又分为基于基本迭代法的最小二乘相位解缠、无权多极格网法、基于FFT/DCT的最小二乘相位解缠、基于误差方程的最小二乘相位解缠。加权最小二乘法包括picard算法、PCG算法、加权多级格网法。4.2.1无权最小二乘法（1）基本迭代法基本迭代法有3种，第一种是。-Jacobi迭代法，第二种是Gauss-Seidel迭代法，第三种是SOR法，其迭代公式分别如下：1）-Jacobi2）Gauss-Seidel3）SOR在每个松弛迭代中，SOR法的计算量与-Jacobi迭代法与Gauss-seidel迭代法相当。如果选择好最佳松弛因子，该方法的收敛速度比-Jacobi代法与Gauss-seidsl迭代法高一个量级。（2）无权多级网络法高斯一赛德尔松弛算法的收敛速度很慢，为了提高运算效率，可采用多级格网技术。高斯一赛德尔松弛算法是一种典型的局域平滑算子，它可以迅速的去除信号中的高频成分，但对低频成分的滤除速度却非常慢，因此导致它很难收敛。多级格网方法的核心思想就是将低频成分转化为高频成分，加速高斯一赛德尔松弛算法的处理速度。这种由低频到高频的转换是通过格网重采样实现的，粗格网的低采样率提高了残差点误差的空间频率。如图的网格金字塔，每一级