相位解缠算法研究

1、、引言合 成 孔 径 雷 达 干 涉 测 量 技 术 ( synthetic aperture radar interferometry, InASR )将合成孔径雷达成像技术与干涉测量技术 成功地进行了结合，利用传感器高度、雷达波长、波束视向及天线基 线距之间的几何关系， 可以精确的测量出图像上每一点的三维位置和 变化信息。合成孔径雷达干涉测量技术是正在发展中的极具潜力的微波遥 感新技术，其诞生至今已近 30 年。起初它主要应用于生成数字高程 模型 (DEM)和制图，后来很快被扩展为差分干涉技术 ( differential InSAR , DInSAR) 并应用于测量微小的地表形变，它已在

2、研究地震形 变、火山运动、 冰川漂移、城市沉降以及山体滑坡等方面表现出极好 的前景。特别， DInSAR 具有高形变敏感度、高空间分辨率、几乎不 受云雨天气制约和空中遥感等突出的技术优势， 它是基于面观测的空 间大地测量新技术， 可补充已有的基于点观测的低空间分辨率大地测 量技术如全球定位系统 (GPS)、甚长基线干涉 (VLBI) 和精密水准等。 尤其 InSAR在地球动力学方面的研究最令人瞩目。二维相位解缠是 InSAR 数据处理流程中重要步骤之一，也是主 要误差来源， 无论是获取数字高程模型还是获取地表形变信息， 其精 确程度都高度依赖于有效的相位解缠。 因此，本人在课程期间对相位 解缠

3、的相关文献进行了阅读。二、 InSAR基本原理用两副雷达天线代替两个光源 S1， S2 ，对地面发射相干信号， 将得到类似的条纹图。 因为雷达信号与光线本质上都是电磁波， 所以 只要保证雷达天线载具运行轨道的稳定， 那么两个信号到达地面上某 一点处的路程差是确定的， 只与该点在地面上的位置有关。 在 InSAR 干涉测量中有两种模式，一种是在载具(卫星或飞机)上搭载一具天 线，而载具两次通过不同轨道航线飞经目标地域上空， 此种称之为单 天线双航过模式； 另一种在载具上搭载两副天线， 只飞经目标地域上 空一次，此种方式称之为双天线单航过模式。 不论是哪种方式都可以 用图 2.2 来模拟并作出几何

4、解释。在测量中两副天线或两次航过接收的数据可以各获得对地面同 一区域的两幅包含幅值与相位信息的二维复数据图像， 分别以 S1 ，S2 表示为j4 r 1S1 |S1 | exp( 1) | S1 |exp( j4 r1)j4 rS2 |S2 |exp( 2) |S2 |exp(2)(2.6)其中 | S1 |和| S2 |表示幅值信息， 1和 2表示相位信息。 将两幅图像 作共轭乘，可得S1 S2* |S1| |S2 |exp( 1 2) |S1| |S2|exp(j4 (r1 r2)(2.7 )j4 (r1 r2) 为两幅图像中相对应的像点的相位差，由路程差决定 的，由余弦定理有r22 r1

5、2 B2 2 Br1 cos()(2.8 )可得2 2 2arccos(r2 r1 B ) ( 2.9 ) 2Br1根据式( 2.7 )的结论，两路雷达波路程差与相位差成正比r r1 r2(2.10 )4式(2.8 )可以进一步得到arccos(2r1 r) r B2)(2.11 )2Br1于是h H r1cos(2.12 )上式中 B 为基线长，由此可以获得地面的高程信息。这里关键 是利用了路程差与相位差成正比这样一个关系， 应该注意的是两天线 接收到的信号的路程差 r 并不很大，但是由于高频的雷达信号的波长 很小，所以 4 r 可以很大，即两个信号的相位差可以比 4大 很多。但是由式( 2

6、.7 )计算相位差时会以 2为模来取值，得到的 相位只会在 ( , 之间，称为相位的主值或缠绕相位， 它与真实相位的关系是相差 2 的整数倍，即有下式的关系2k k=0 ，1，2(2.13 )根据缠绕相位得到真实相位的处理过程就叫做相位解缠，是 InSAR 干涉测量的关键步骤。三、相位解缠基本原理3.1 引言 在上节提到利用相位差能获得精确的路程差进而获得地面的高 程信息，因此获得准确的相位差就是实现测量的关键。 由于复数对其 相位的周期性， InSAR 根据两幅 SAR 复图像获得的干涉相位差值是 被周期折叠后位于 ( , 之间的相位主值， 它与真实的相位差值 之间存在着 2 k 差别。由式

7、 (2.13) 可以表示它们之间的基本关系。 其中 代表解缠相位， 代表缠绕相位。必须对 进行相位解缠，恢 复被模糊掉的相位周期， 获得目标在两次成像中的真实相位差， 才能 得到目标的正确高度信息。相位解缠是 InSAR 三维成像处理中的关 键步骤之一，其准确程度将直接决定数字高程图( DEM)和地表形变 探测的精度。3.2 相位缠绕和解缠理想情况下，图像的采样率满足 Nyquist 采样定理，采样频率 必须大于信号最高频率的两倍， 解缠绕的干涉相位中相邻像素点之间 的相位差值不可能超过半个周期(一个) 。当满足此条件时必然能 由缠绕相位解缠出正确的解缠绕相位，并且可以通过积分进行解缠。记 (

8、m) 为周期缠绕前的真实相位值， (m) 为相应的缠绕相位，定义相位缠绕算子 ，相位缠绕的过程可以用式( 3.1 )表示( (m) (m) (m) 2 k(m)(3.1 )( (m)结果是得到主值属于 ( , 区间的缠绕相位。 定义差分算子 ，根据 Nyquist 采样定理对于解缠相位有(m) (m 1) (m)(m) 对相邻缠绕相位进行差分运算得(m) (m 1) (m) (m) 2 k(m) 对该相位差也使用缠绕算子得(m) (m) 2 k(m) 2 k(m) 根据缠绕算子的定义，其结果必须属于 ( (m) 也必须属于 ( , 区间，所以有k(m) k(m) 0式(3.3 )变为(m) (

9、m) 由式( 3.2 )可得 m1(m) (0) (n) n0(3.2)(3.3)(3.4) 区间，而 (3.5)(3.6)(3.7)由式( 3.7 )可以看出，通过对相邻缠绕相位之差积分可以实现 相位解缠，条件是满足 Nyquist 采样定理。对于一维的情况， 可以简单的使用如下的公式进行解缠计算， 记 (m) 为周期缠绕前的真实相位值， (m) 为相应的缠绕相位， 计算干涉图中一个点到下一个点的相位变化， 即计算相位梯度， 然后从一固 定点开始积分使相位值的变化平稳连续，从而恢复失去的相位周期。 即下式：(1) (1) (m 1) (m) (m)(3.8 )若有如下的一维相位序列0.2 ,

10、 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 。 由于一维序列的积分路径是唯一的， 所以其解也是唯一的。 但由 于是逐个积分， 如果受到相位噪声的影响， 或者碰到地形起伏本来就 不满足相邻缠绕相位差的绝对值小于 的条件，使其中一

11、点的解缠 绕相位发生错误， 则错误会后向传播， 导致之后所有相位的解缠结果 与真实相位相差甚远。为了说明相位缠绕与解缠原理，选取如图 3.1 所示的人工模拟 的简单缠绕相位图进行解释。 在理想状况下， 发生缠绕的干涉相位呈 现周期性变化，由 渐变到 ，然后由 突变为 ，如此反复， 从图像上表现为灰度值由浅渐渐变深， 然后突变为浅色， 再向深色渐 变，形成如图 3.1 (a)所示的条纹图。从图 3.1 ( a)中沿 y 轴方向取一条一维数据， 以像素位置为横坐标， 以灰度强弱代表的相位值 为纵坐标将其表示出来将如图 3.1 (c)所示，其形状如锯齿状，代 表了图 3.1 ( a)中黑白交替变换的

12、条纹。理想情况下的解缠绕只需进行简单的积分将突变消除， 整幅图像的条纹变成了连续的面， 相位 恢复连续变化。如图 3.1 (b )所示。在图 3.1 (b)中也取一条一维数据在坐标图中画出，将如图 3.1 (d )所示四、常用相位解缠算法4.1 常用相位解缠算法概述到目前为止， 针对相位解缠问题已经提出很多解决方案。 主要的解缠算法大致可以分为三类：一类可以称之为路径跟踪解缠算法， 他们的共同特点是采用路径 积分来 实现相位 解缠，以 1988 年 Goldstein 提 出的 枝切 法 ( Branch-Cut )为代表。枝切法通过探测残差点，用枝切线连接残差 点，然后进行路径积分来实现解缠

13、， 在路径积分时以不穿越枝切线为 原则。 Wei Xu 和 Cumming提出的区域生长法( Region-Growing )不 考虑残差点，不布置枝切线， 而是依据额外信息将干涉图划分为高质 量低质量区域， 在各个区域内按照从高质量像元到低质量像元的方向 进行路径积分 。 Flynn 的掩模分割法( Mask- Cut )和最小不连续 法( Minimum-discontinuity )等也属于该类算法。另一类算法着眼于整体， 采用最优化的思想， 寻求最小二乘意义 下的最优解缠结果， 包括用 FFT/DCT方法求解的无加权最小二乘算法， Pritt 的多重网格迭代法求解加权最小二乘相位解缠法

14、， Ghiglia 的 最小范数法 等。这类算法不探测残差点，不布置枝切线，通过建立 一个离散型泊松目标函数， 并用各种数学的方法求解它以实现相位解 缠。第三类方法为最小费用流方法，以 Costantini 的基于网络规划 的解缠方法为代表， 引入图论中的网络模型， 将解缠问题转变为解一 个网络最小费用流的问题， 利用网络规划理论中成熟高效的算法求解。4.2 基于路径的相位解缠算法两幅 SAR图像经过干涉以后，我们可以获得一幅缠绕相位图像，各像元上的值为对应的干涉相位的主值。 根据 Nyquist 定理，当相邻 像元上的相位差小于二时， 可以通过积分的算法来恢复相位的真实值。 基于路径跟踪的相

15、位解缠算法就是通过积分相邻缠绕相位的差分值 来恢复相位的真实值的。假设我们己知在像元 r0 上的相位，那么在其 它像元 r 上的相位可以通过以下公式来获得 :(r ) dr (r0)c(4.1)符号 (r)为像元 r 上的解缠相位， ( r0 )为像元 r0上的已知解缠相 位，C为积分路径，根据积分理论：I F (r)drc(4.2)上式 F(r) 为积分函数， C为积分路径，(4.2)的线性积分不仅依 赖于积分路径 C的起点和终点， 还依赖于积分路径 C本身。要使积分 与路径无关，则要求一下闭合积分成立?F (r)dr 0(4.3)在二维相位解缠中，公式 (4.3) 常用来作为探测积分是否与

16、路径 无关的条件。 InSAR缠绕相位数据中，不是所有的积分路径都满足公 式 (4.3) ，有些像元上的缠绕相位数据由于受到噪声的影响，或者由 于其它的原因，导致通过这些像元的闭合积分不能满足公式 (4.3) ， 这些像元上的相位在 InSAR 中被称为“残差点 (residue) ”，或者“电 荷”( 具有正负性，见随后的讨论 )，在路径跟踪的相位解缠算法中， 关键的问题在于如何判断这些电荷并将它们相连 (称为“分枝”) 以达 到正负抵消，且防止积分路径穿过这些分枝。路径跟踪法的基本策略是将可能的误差传递限制在噪声区内， 通 过选择合适的积分路径，隔绝噪声区，阻止相位误差的全程传递。几 十年

17、来，研究者研究出了许多的相位解缠算法，至今为止，基于路径 跟踪的相位解缠算法有枝切法、区域法、 Mask-cut 算法、像元扩散 法、最小生成树法、条纹检测法、区域生长法，最小不连续算法 ( 简 称 Rynn算法) 等算法。路径跟踪的相位解缠算法一般步骤如下所示 :输入 : 缠绕相位步骤 1:相位连续性 /不连续性检测 :识别残差点，生成枝切线。步骤 2:计算/ 建立相位质量图。步骤 3:相位积分 :在枝切线周围或在质量图的指导下处理。下面对其中比较典型的算法作详细介绍 :4.2.1Goldstein 枝切法Goldstein 枝切法是较经典的路径跟踪法， 是 1988年 Goldstein

18、等人提出的， 它识别正负残差点， 并连接邻近的残差点对或多个残差 点，实现残差点“电荷”平衡，生成最优的枝切线，确定积分路径， 防止误差沿积分路径传递。基本步骤为 :(l) 识别残差点 ;(2) 生成枝切线 ;(3) 绕过枝切线进行积分。具体如下 :首先按一定的顺序寻找残差点，定义一个 2x2 像元的缠绕相位为节点， 将四个像元串接起来， 即为影像中的最小闭合路径(dosedloop) 。沿最小闭合路径将缠绕相位梯度累加起来，如果之和 为零则这四个点是一致的，否则左上角的像元就称为残差点(residual) 。缠绕相位节点图与最小闭合路径图如图所示缠绕相位节点，最小闭合路径图 ( 每个像元数值

19、乘以劫才表示缠绕相位真实值 ) ，i 为各方向相位差 :i, j 1i,ji 1, j 1i, j 1i 1, ji 1, j 1i,ji 1, j计算得1 0.12 0.33 0.04 0.24i0i1说明相位是致的再给出一个例子：缠绕相位节点图与最小闭合路径图如图所示40.30.30.30.1i1i1则左上角的像元为残差点，当找到第一个残差点以后， 从该残差点开始，继续搜索，找到下一个残差点后，用一枝切线将两者连接起 来计算残数和， 如果和为零则完成了该树枝的生长， 继续搜索直至搜 索完全部残差点， 如果和不为零则不断加入残差点每次计算总的残数 和，直至和为零。在 Goldstein 的枝

20、切法中，有两个步骤是极为关键的 :(l) 当搜索 窗口找到新的残差点， 无论该残差点是否与其他的残差点相连， 都将 该残差点与窗口中心的残差点相连 ;(2) 当搜索窗口到达图像的边界， 则将残差点与边界相连，以阻止积分路径。枝切法最大的优点是 : 在实际计算中速度比较快 ; 在噪声比较低、 残差点比较少的情况下，精确度非常高。缺点是 : 当残差点较多且分 布较密集时， 该算法难以准确连接枝切线， 导致无法选择合理的积分 路径，有时会造成错误的阮跳跃，导致误差的传递。但由于该算法的 速度优势，使之成为一种常用的相位解缠算法。4.2.2 质量引导法这种算法不识别残差点， 也不设置枝切线。 而是在进

21、行相位解缠 时，通过相位质量图 (quality map)来定义相位数据的质量，将积分 路径总是沿“高”质量的像元进行，最后解缠“低”质量像元。质量引导法的关键步骤就是在相位质量图的引导下进行像元扩 散，其基本操作过程如下 : 从高质量像元点出发，检测它的四个邻近 像元，对邻近像元进行解缠，将解缠后的像元的邻接像元(未解缠 )存储在“邻接列”中，依据相位质量从“邻接列”移出高质量像元进 行相位解缠，更新“邻接列” ，重复上述步骤直至所有的像元解缠完 毕。质量引导法成功地进行相位解缠的前提是必须有可靠的相位质 量图。相位质量图主要有四种 : 相干图、伪相干图、相位导数变化图、 最大相位梯度图、掩

22、模图 (mask) 。(1) 相干图 最常用的质量图是相干图，相干值的高低表明图像不同 区域的相干性，是最直观的干涉质量评价图。同时相干系数 的变化也表征了在图像获取期间地物的变化情况，所以相干 图也用于地物的分类等。(2) 伪相干图 (Pseudo-correlation)当无法获得 InSAR 图像对的强度值时，常常用伪相干图 来模拟相干图。 这时把 InSAR图像对的强度定义为 1，那么伪 相干定义为：( cos i,j)2 ( sin i,j )2zm,n2kk 为视数。伪相干图的一个最大缺陷是： 对于陡峭地形区， 它 标志为低质量数据区 ( 即使这些相位数据质量很好， 并且没有 噪声

23、) 。这时就需要新的质量图，用来评价相位导数的统计变 化特征。(3) 相位导数变化图 (Phase Derivative Variance) 相位导数变化定义如下 :zm,n( xi, jm,n)2m,n( yi ,jym,n)2m,nk2相位导数变化不同于相干和伪相干。 例如，在倾斜地表，如果相位变化率保持一定， 则相位导数变化为 0，而伪相干不为 0。从严格意义上讲，相位导数变化表征的是相位数据的 “坏” (badness) ，而不是“好” (goodness) ，但我们可以假定 : 如 果相位导数变化是可以忽略的话， 那么相位数据就是好数据。 实验表明，在无法获得相干图的情况下，相位导数

24、变化图是 最可靠的相位质量图。(4) 最大相位梯度图从相位图上可以看出，在噪声相位区往往相位梯度也很 大，所以可以用最大相位梯度来表征相位数据的质量。一般 最大梯度定义为：xymax i, jmax i,j最大相位梯度图也有伪相干的缺陷， 在地形很陡峭 ( 即相 位变化显著，但无噪声 ） 也表征为低质量数据。质量引导法完全依赖于质量图像来指导解缠路径的选择， 因此在缺乏高质量的质量图像时，解缠效果将会很不理想。 另外该算法不识别残差点，因此解缠路径就不可避免的可能 会包围非平衡的残差点， 这样就有可能会产生 2k 二的周期累 加错误。（5）掩膜图（ mask）在解缠过程中，有的区域失相关严重

25、（如水面 ）导致相位 不连续，有的区域地势平坦不需要进行滤波。 在这种情况下， 我们可以制作一个掩模图将这些区域掩盖起来，使解缠和滤 波不涉及这些区域。掩模图的生成方法一般是采用一定的闲 值来进行判断，或者手工确定需要掩盖的区域，通常被掩盖 区域的干涉图像元用 0 来表示，未被掩盖的像元 1 来表示。 第四章对西安数据的解缠过程中均采用了掩模图，即掩盖掉 的像元不参与解缠，以防止误差的传播。质量引导法解缠结果的好坏很大程度上依赖于质量图的 “好”与“坏”，如果有好的质量图引导，能得到优于枝切法 的解缠结果。因为不设置枝切线，积分路径就可能包围残差 点导致错误的 2k的累积错误。4.2.3 掩膜

26、切口法（ mask-cut 法）1966年 Rynn首先对质量图像引导分割路径的算法作了详细的描 述，在该方法中使用逐步增长的像素掩模来连接残差点， 这种像素掩 模称为掩模切口， 与分割路径的作用相似。 掩模切口算法与质量引导 法有些相似， 它在有些地方可以看作是后者的逆向算法， 它并不是从 高质量的区域开始解缠， 而是从残差点开始， 沿低质量的区域逐渐扩 展像素掩模。掩模切口的扩展过程一直到它连接了等量的正残差点和 负残差点或者到达图像边界时为止， 在前一种情况下， 残差点是平衡 的，相位解缠与积分路径无关，而在后一种情况下，掩模将残差点隔 离开来，可以保证不会有解缠路径可能包围这些残差点。

27、该算法实际上是将 Goldstein 枝切法和基于质量图路径积分法 的结合起来，掩模分害算法虽然也识别残差点并产生分割路径，但 总的来讲， 该算法还是在质量图像的指导下实现的， 因此算法的结果 在很大程度上也仍然依赖于高质量的质量图像。 在残留点不仅仅分布 在低质量相位区的情况下， mask-cut 算法与质量指导的路径跟踪算 法都无法正常工作，在这种情况下，枝切法更为有效。其解缠步骤为 :(l) 识别残差点；(2) 生成 mask-cut ；(3) 细化 mask-cut ；(4) 沿 mask-cut 的路径积分。4.3 最小范数法 理想情况下解缠相位梯度等于缠绕相位梯度的假设， 相位解缠

28、可 以看成一个优化问题。 最小范数法将相位解缠问题转化为数学上的最小范数问题， 目前使用最广泛的是最小二乘法， 最小二乘法是一种广 泛使用的优化方法，有无权和加权两种形式。无权形势下，相位解缠 是求取一个平滑的解缠相位。就是求解一个 Neumann件下的 Poisson 方程。可以通过离散余弦变换 DCT、离散傅里叶变换 FFT或无权多级 格网法来有效的解决。 由于干涉图上各像素点相关系数差别较大， 存 在相位的不连续， 因此无权最小二乘虽然获得了平滑的相位解缠曲面， 但造成局部的噪声在最小均方意义下的全局传播而产生与真实相位 值偏差较大的解。 加权最小二乘法可以在一定程度上弥补无权最小二 乘

29、法的这一缺陷。 例如共轭梯度法使用的权系数是经过二值化的质量 图，将干涉图中由于残余点的存在而破坏的区域赋予零权， 阻止它们 对相位解缠的破坏。不带权的最小二乘法相位解缠又分为基于基本迭代法的最小二 乘相位解缠、无权多极格网法、基于 FFT/DCT的最小二乘相位解缠、 基于误差方程的最小二乘相位解缠。加权最小二乘法包括 picard 算 法、 PCG算法、加权多级格网法。4.2.1 无权最小二乘法（1）基本迭代法基本迭代法有 3 种，第一种是。 -Jacobi 迭代法，第二种是Gauss-Seidel 迭代法，第三种是SOR法，其迭代公式分别如下：1) -Jacobik1i, j(1) ik,

30、 j( ik 1,jkkki 1,j i, j 1 i,j 1i,j)/42) Gauss-Seidelk1i, j 1i,j)/4ik,j1(1 ) ik, j( ik1,jki 1, jik,j11ik,j11i,j) /4在每个松弛迭代中，SOR法的计算量与-Jacobi 迭代法与Gauss-seidel 迭代法相当如果选择好最佳松弛因子，该方法的收敛速度比 -Jacobi 代法与Gauss-seidsl 迭代法高一个量级。k 1 k k k1i,j ( i 1,j i 1,j i, j 13) SOR（2）无权多级网络法高斯一赛德尔松弛算法的收敛速度很慢， 为了提高运算效率， 可 采用多级格网技术。 高斯一赛德尔松弛算法是一种典型的局域平滑算 子，它可以迅速的去除信号中的高频成分， 但对低频成分的滤除速度 却非常慢， 因此导致它很难收敛。 多级格网方法的核心思想就是将低 频成分转化为高频成分， 加速高斯一赛德尔松弛算法的处理速度。 这 种由低频到高频的转换是通过格网重采样实现的， 粗格网的低采样率 提高了残差点误差