相位解缠算法研究

一、引言
合成孔径雷达干涉测量技术（synthetic aperture radar interferometry, InASR）将合成孔径雷达成像技术与干涉测量技术成功地进行了结合，利用传感器高度、雷达波长、波束视向及天线基线距之间的几何关系，可以精确的测量出图像上每一点的三维位置和变化信息。

合成孔径雷达干涉测量技术是正在发展中的极具潜力的微波遥感新技术，其诞生至今已近30年。

起初它主要应用于生成数字高程模型(DEM)和制图，后来很快被扩展为差分干涉技术( differential InSAR , DInSAR)并应用于测量微小的地表形变，它已在研究地震形变、火山运动、冰川漂移、城市沉降以及山体滑坡等方面表现出极好的前景。

特别，DInSAR具有高形变敏感度、高空间分辨率、几乎不受云雨天气制约和空中遥感等突出的技术优势，它是基于面观测的空间大地测量新技术，可补充已有的基于点观测的低空间分辨率大地测量技术如全球定位系统(GPS)、甚长基线干涉(VLBI)和精密水准等。

尤其InSAR在地球动力学方面的研究最令人瞩目。

二维相位解缠是InSAR 数据处理流程中重要步骤之一，也是主要误差来源，无论是获取数字高程模型还是获取地表形变信息，其精确程度都高度依赖于有效的相位解缠。

因此，本人在课程期间对相位解缠的相关文献进行了阅读。

二、InSAR基本原理
用两副雷达天线代替两个光源
S，2S，对地面发射相干信号，
1
将得到类似的条纹图。

因为雷达信号与光线本质上都是电磁波，所以只要保证雷达天线载具运行轨道的稳定，那么两个信号到达地面上某一点处的路程差是确定的，只与该点在地面上的位置有关。

在 InSAR 干涉测量中有两种模式，一种是在载具（卫星或飞机）上搭载一具天线，而载具两次通过不同轨道航线飞经目标地域上空，此种称之为单天线双航过模式；另一种在载具上搭载两副天线，只飞经目标地域上空一次，此种方式称之为双天线单航过模式。

不论是哪种方式都可以用图 来模拟并作出几何解释。

在测量中两副天线或两次航过接收的数据可以各获得对地面同一区域的两幅包含幅值与相位信息的二维复数据图像，分别以1S ，2S 表示为
2
22224||exp()||exp()j r S S S πϕλ== （）
其中1||S 和2||S 表示幅值信息，1ϕ和2ϕ表示相位信息。

将两幅图像
作共轭乘，可得
*12121212124()||||exp()||||exp(
)j r r S S S S S S πϕϕλ-⋅=⋅-=⋅ （）
124()j r r πλ-为两幅图像中相对应的像点的相位差，由路程差决定的，由余弦定理有
2222112cos()r r B Br αβ=+++ （） 可得
222
211
arccos()2r r B Br βα--=- （） 根据式（）的结论，两路雷达波路程差与相位差成正比
124r r r φλπ
∆∆=-= （） 式（）可以进一步得到
2
11
(2)arccos()2r r r B Br βα+∆∆-=- （） 于是
1cos h H r β=- （） 上式中 B 为基线长，由此可以获得地面的高程信息。

这里关键是利用了路程差与相位差成正比这样一个关系，应该注意的是两天线接收到的信号的路程差r ∆并不很大，但是由于高频的雷达信号的波长 λ很小，所以4r
πφλ∆∆=可以很大，即两个信号的相位差可以比4π大
很多。

但是由式（）计算相位差时会以2π为模来取值，得到的相位只会在 ( π ,π]之间，称为相位的主值或缠绕相位，它与真实相位的关系是相差 2π 的整数倍，即有下式的关系
2k φϕπ=+ k=0，±1，±2…… （） 根据缠绕相位得到真实相位的处理过程就叫做相位解缠，是 InSAR 干涉测量的关键步骤。

三、相位解缠基本原理
引言
在上节提到利用相位差能获得精确的路程差进而获得地面的高程信息，因此获得准确的相位差就是实现测量的关键。

由于复数对其相位的周期性，InSAR 根据两幅 SAR 复图像获得的干涉相位差值是被周期折叠后位于 ( π ,π]之间的相位主值，它与真实的相位差值之间存在着 2 k π差别。

由式可以表示它们之间的基本关系。

其中φ
代表解缠相位， 代表缠绕相位。

必须对 进行相位解缠，恢复被模糊掉的相位周期，获得目标在两次成像中的真实相位差，才能得到目标的正确高度信息。

相位解缠是 InSAR 三维成像处理中的关键步骤之一，其准确程度将直接决定数字高程图（DEM ）和地表形变探测的精度。

相位缠绕和解缠
理想情况下，图像的采样率满足 Nyquist 采样定理，采样频率必须大于信号最高频率的两倍，解缠绕的干涉相位中相邻像素点之间的相位差值不可能超过半个周期（一个π）。

当满足此条件时必然能由缠绕相位解缠出正确的解缠绕相位，并且可以通过积分进行解缠。

记 φ (m)为周期缠绕前的真实相位值， (m)为相应的缠绕相位，定义相位缠绕算子ϖ ，相位缠绕的过程可以用式（）表示
(())()()2()m m m k m ϖφϕφπ==+ （） 结果是得到主值属于 ( π ,π]区间的缠绕相位。

定义差分算子 Δ，根据 Nyquist 采样定理对于解缠相位有
()(1)()()m m m m φφφπφπ
∆=+--<∆≤ （） 对相邻缠绕相位进行差分运算得
()(1)()()2()m m m m k m ϕϕϕϕπ∆=+-=∆+∆ （） 对该相位差也使用缠绕算子得
[][]'()()2()2()m m k m k m ϖϕϕππ∆=∆+∆+ （） 根据缠绕算子的定义，其结果必须属于 ( π ,π]区间，而 Δφ (m)也必须属于( π ,π]区间，所以有
'()()0k m k m ∆+= （） 式（）变为
[]()()m m φϖϕ∆=∆ （） 由式（）可得
[]1
0()(0)()m n m n φφϖϕ-==+∆∑ （）
由式（）可以看出，通过对相邻缠绕相位之差积分可以实现相位解缠，条件是满足 Nyquist 采样定理。

对于一维的情况，可以简单的使用如下的公式进行解缠计算，记 φ (m)为周期缠绕前的真实相位值， (m)为相应的缠绕相位，计算干涉图中一个点到下一个点的相位变化，即计算相位梯度，然后从一固定点开始积分使相位值的变化平稳连续，从而恢复失去的相位周期。

即下式：
(1)(1)φϕ= (1)()()m m m φφ+=+∆ （） 若有如下的一维相位序列
π , π , π , π , π , π , π
以相邻的 π , π 两个数据为例， 0 .5π π)= π，因为 1. 3π < π所以 Δ ( m )= π +2π=π，将 0. 7π加上前一个解缠结果 0. 8π得到该位置的解缠结果为 1 .5π。

其他照此进行，从左向右解缠后的序列为：
π , π , π , π , π , π , π 。

由于一维序列的积分路径是唯一的，所以其解也是唯一的。

但由于是逐个积分，如果受到相位噪声的影响，或者碰到地形起伏本来就
不满足相邻缠绕相位差的绝对值小于π的条件，使其中一点的解缠绕相位发生错误，则错误会后向传播，导致之后所有相位的解缠结果与真实相位相差甚远。

为了说明相位缠绕与解缠原理，选取如图所示的人工模拟的简单缠绕相位图进行解释。

在理想状况下，发生缠绕的干涉相位呈现周期性变化，由π渐变到π，然后由π突变为π，如此反复，从图像上表现为灰度值由浅渐渐变深，然后突变为浅色，再向深色渐变，形成如图（a）所示的条纹图。

从图（a）中沿y 轴方向取一条一维数据，以像素位置为横坐标，以灰度强弱代表的相位值为纵坐标将其表示出来将如图（c）所示，其形状如锯齿状，代表了图（a）中黑白交替变换的条纹。

理想情况下的解缠绕只需进行简单的积分将突变消除，整幅图像的条纹变成了连续的面，相位恢复连续变化。

如图(b）所示。

在图(b）中也取一条一维数据在坐标图中画出，将如图(d）所示。

四、常用相位解缠算法
常用相位解缠算法概述
到目前为止，针对相位解缠问题已经提出很多解决方案。

主要的解缠算法大致可以分为三类：
一类可以称之为路径跟踪解缠算法，他们的共同特点是采用路径积分来实现相位解缠，以1988年Goldstein提出的枝切法（Branch-Cut）为代表。

枝切法通过探测残差点，用枝切线连接残差点，然后进行路径积分来实现解缠，在路径积分时以不穿越枝切线为原则。

Wei Xu
和Cumming 提出的区域生长法（Region -Growing ）不考虑残差点，不布置枝切线，而是依据额外信息将干涉图划分为高质量低质量区域，在各个区域内按照从高质量像元到低质量像元的方向进行路径积分 。

Flynn 的掩模分割法（Mask - Cut ）和最小不连续法（Minimum -discontinuity ）等也属于该类算法。

另一类算法着眼于整体，采用最优化的思想，寻求最小二乘意义下的最优解缠结果，包括用FFT/DCT 方法求解的无加权最小二乘算法，Pritt 的多重网格迭代法求解加权最小二乘相位解缠法，Ghiglia 的最小范数法[等。

这类算法不探测残差点，不布置枝切线，通过建立一个离散型泊松目标函数，并用各种数学的方法求解它以实现相位解缠。

第三类方法为最小费用流方法，以Costantini 的基于网络规划的解缠方法为代表，引入图论中的网络模型，将解缠问题转变为解一个网络最小费用流的问题，利用网络规划理论中成熟高效的算法求解。

基于路径的相位解缠算法
两幅SAR 图像经过干涉以后，我们可以获得一幅缠绕相位图像，各像元上的值为对应的干涉相位的主值。

根据Nyquist 定理，当相邻像元上的相位差小于二时，可以通过积分的算法来恢复相位的真实值。

基于路径跟踪的相位解缠算法就是通过积分相邻缠绕相位的差分值来恢复相位的真实值的。

假设我们己知在像元0r 上的相位，那么在其
它像元r 上的相位可以通过以下公式来获得:
0()()
c r dr r ϕϕϕ=∇+⎰
（）
符号()r ϕ为像元r 上的解缠相位，0()r ϕ为像元0r 上的已知解缠相
位，C 为积分路径，根据积分理论：
()c I F r dr
=⎰
（）
上式()F r 为积分函数，C 为积分路径，（）的线性积分不仅依赖于积分路径C 的起点和终点，还依赖于积分路径C 本身。

要使积分与路径无关，则要求一下闭合积分成立
()0
F r dr =⎰Ñ （）
在二维相位解缠中，公式常用来作为探测积分是否与路径无关的条件。

InSAR 缠绕相位数据中，不是所有的积分路径都满足公式，有些像元上的缠绕相位数据由于受到噪声的影响，或者由于其它的原因，导致通过这些像元的闭合积分不能满足公式，这些像元上的相位在InSAR 中被称为“残差点(residue)”，或者“电荷”(具有正负性，见随后的讨论)，在路径跟踪的相位解缠算法中，关键的问题在于如何判断这些电荷并将它们相连(称为“分枝”)以达到正负抵消，且防止积分路径穿过这些分枝。

路径跟踪法的基本策略是将可能的误差传递限制在噪声区内，通过选择合适的积分路径，隔绝噪声区，阻止相位误差的全程传递。

几十年来，研究者研究出了许多的相位解缠算法，至今为止，基于路径跟踪的相位解缠算法有枝切法、区域法、Mask -cut 算法、像元扩散法、
最小生成树法、条纹检测法、区域生长法，最小不连续算法(简称Rynn算法)等算法。

路径跟踪的相位解缠算法一般步骤如下所示:
输入:缠绕相位
步骤1:相位连续性/不连续性检测:识别残差点，生成枝切线。

步骤2:计算/建立相位质量图。

步骤3:相位积分:在枝切线周围或在质量图的指导下处理。

下面对其中比较典型的算法作详细介绍:
枝切法
Goldstein枝切法是较经典的路径跟踪法，是1988年Goldstein 等人提出的，它识别正负残差点，并连接邻近的残差点对或多个残差点，实现残差点“电荷”平衡，生成最优的枝切线，确定积分路径，防止误差沿积分路径传递。

基本步骤为:
(l)识别残差点;
(2)生成枝切线;
(3)绕过枝切线进行积分。

具体如下:首先按一定的顺序寻找残差点，定义一个2x2像元的缠绕相位为节点，将四个像元串接起来，即为影像中的最小闭合路径(dosedloop)。

沿最小闭合路径将缠绕相位梯度累加起来，如果之和为零则这四个点是一致的，否则左上角的像元就称为残差点(residual)。

缠绕相位节点图与最小闭合路径图如图所示:
缠绕相位节点，最小闭合路径图(每个像元数值乘以劫才表示缠
绕相位真实值)，
为各方向相位差:
i
计算得
说明相位是一致的。

再给出一个例子：
缠绕相位节点图与最小闭合路径图如图所示:
则左上角的像元为残差点，当找到第一个残差点以后，从该残差点开始，继续搜索，找到下一个残差点后，用一枝切线将两者连接起来计算残数和，如果和为零则完成了该树枝的生长，继续搜索直至搜索完全部残差点，如果和不为零则不断加入残差点每次计算总的残数和，直至和为零。

在Goldstein的枝切法中，有两个步骤是极为关键的:(l)当搜索窗口找到新的残差点，无论该残差点是否与其他的残差点相连，都将该残差点与窗口中心的残差点相连;(2)当搜索窗口到达图像的边界，则将残差点与边界相连，以阻止积分路径。

枝切法最大的优点是:在实际计算中速度比较快;在噪声比较低、残差点比较少的情况下，精确度非常高。

缺点是:当残差点较多且分布较密集时，该算法难以准确连接枝切线，导致无法选择合理的积分路径，有时会造成错误的阮跳跃，导致误差的传递。

但由于该算法的速度优势，使之成为一种常用的相位解缠算法。

质量引导法
这种算法不识别残差点，也不设置枝切线。

而是在进行相位解缠时，通过相位质量图(quality map)来定义相位数据的质量，将积分路径总是沿“高”质量的像元进行，最后解缠“低”质量像元。

质量引导法的关键步骤就是在相位质量图的引导下进行像元扩散，其基本操作过程如下:从高质量像元点出发，检测它的四个邻近像元，对邻近像元进行解缠，将解缠后的像元的邻接像元(未解缠)存储在“邻接列”中，依据相位质量从“邻接列”移出高质量像元进行相位解缠，更新“邻接列”，重复上述步骤直至所有的像元解缠完毕。

质量引导法成功地进行相位解缠的前提是必须有可靠的相位质量图。

相位质量图主要有四种:相干图、伪相干图、相位导数变化图、最大相位梯度图、掩模图(mask)。

(1)相干图
(2)最常用的质量图是相干图，相干值的高低表明图像不同
区域的相干性，是最直观的干涉质量评价图。

同时相干系数
的变化也表征了在图像获取期间地物的变化情况，所以相干
图也用于地物的分类等。

(3)伪相干图(Pseudo-correlation)
(4)当无法获得InSAR图像对的强度值时，常常用伪相干图
来模拟相干图。

这时把InSAR图像对的强度定义为1，那么
伪相干定义为：
(5)
z=
,m n
(6)k为视数。

伪相干图的一个最大缺陷是：对于陡峭地形区，它
标志为低质量数据区(即使这些相位数据质量很好，并且没有
噪声)。

这时就需要新的质量图，用来评价相位导数的统计变
化特征。

(7) 相位导数变化图 (Phase Derivative Variance)
(8) 相位导数变化定义如下:
(9) ,m n z =
(10) 相位导数变化不同于相干和伪相干。

例如，在倾斜地表，如果相位变化率保持一定，则相位导数变化为0，而伪相干不为0。

从严格意义上讲，相位导数变化表征的是相位数据的“坏”(badness)，而不是“好”(goodness)，但我们可以假定:如果相位导数变化是可以忽略的话，那么相位数据就是好数据。

实验表明，在无法获得相干图的情况下，相位导数变化图是最可靠的相位质量图。

(11) 最大相位梯度图
(12) 从相位图上可以看出，在噪声相位区往往相位梯度也很大，所以可以用最大相位梯度来表征相位数据的质量。

一般最大梯度定义为：
(13) {},max x
i j ∆ {},max y
i j ∆
最大相位梯度图也有伪相干的缺陷，在地形很陡峭(即相
位变化显着，但无噪声)也表征为低质量数据。

质量引导法完全依赖于质量图像来指导解缠路径的选择，
因此在缺乏高质量的质量图像时，解缠效果将会很不理想。

另外该算法不识别残差点，因此解缠路径就不可避免的可能会包围非平衡的残差点，这样就有可能会产生2k 二的周期累加错误。

(14)掩膜图（mask）
(15)在解缠过程中，有的区域失相关严重(如水面)导致相
位不连续，有的区域地势平坦不需要进行滤波。

在这种情况
下，我们可以制作一个掩模图将这些区域掩盖起来，使解缠
和滤波不涉及这些区域。

掩模图的生成方法一般是采用一定
的闲值来进行判断，或者手工确定需要掩盖的区域，通常被
掩盖区域的干涉图像元用0来表示，未被掩盖的像元1来表
示。

第四章对西安数据的解缠过程中均采用了掩模图，即掩
盖掉的像元不参与解缠，以防止误差的传播。

(16)质量引导法解缠结果的好坏很大程度上依赖于质量图
的“好”与“坏”，如果有好的质量图引导，能得到优于枝切
法的解缠结果。

因为不设置枝切线，积分路径就可能包围残
差点导致错误的2kπ的累积错误。

掩膜切口法（mask-cut法）
1966年Rynn首先对质量图像引导分割路径的算法作了详细的描述，在该方法中使用逐步增长的像素掩模来连接残差点，这种像素掩模称为掩模切口，与分割路径的作用相似。

掩模切口算法与质量引导法有些相似，它在有些地方可以看作是后者的逆向算法，它并不是从高质量的区域开始解缠，而是从残差点开始，沿低质量的区域逐渐扩展像素掩模。

掩模切口的扩展过程一直到它连接了等量的正残差点和负残差点或者到达图像边界时为止，在前一种情况下，残差点是平衡的，相位解缠与积分路径无关，而在后一种情况下，掩模将残差点隔
离开来，可以保证不会有解缠路径可能包围这些残差点。

该算法实际上是将Goldstein枝切法和基于质量图路径积分法的结合起来，掩模分害」算法虽然也识别残差点并产生分割路径，但总的来讲，该算法还是在质量图像的指导下实现的，因此算法的结果在很大程度上也仍然依赖于高质量的质量图像。

在残留点不仅仅分布在低质量相位区的情况下，mask-cut算法与质量指导的路径跟踪算法都无法正常工作，在这种情况下，枝切法更为有效。

其解缠步骤为:
(l)识别残差点；
(2)生成mask-cut；
(3)细化mask-cut；
(4)沿mask-cut的路径积分。

最小范数法
理想情况下解缠相位梯度等于缠绕相位梯度的假设，相位解缠可以看成一个优化问题。

最小范数法将相位解缠问题转化为数学上的最小范数问题，目前使用最广泛的是最小二乘法，最小二乘法是一种广泛使用的优化方法，有无权和加权两种形式。

无权形势下，相位解缠是求取一个平滑的解缠相位。

就是求解一个Neumann件下的Poisson 方程。

可以通过离散余弦变换DCT、离散傅里叶变换FFT或无权多级格网法来有效的解决。

由于干涉图上各像素点相关系数差别较大，存在相位的不连续，因此无权最小二乘虽然获得了平滑的相位解缠曲面，但造成局部的噪声在最小均方意义下的全局传播而产生与真实相
位值偏差较大的解。

加权最小二乘法可以在一定程度上弥补无权最小二乘法的这一缺陷。

例如共轭梯度法使用的权系数是经过二值化的质量图，将干涉图中由于残余点的存在而破坏的区域赋予零权，阻止它们对相位解缠的破坏。

不带权的最小二乘法相位解缠又分为基于基本迭代法的最小二乘相位解缠、无权多极格网法、基于FFT/DCT的最小二乘相位解缠、基于误差方程的最小二乘相位解缠。

加权最小二乘法包括picard算法、PCG算法、加权多级格网法。

无权最小二乘法
（1）基本迭代法
基本迭代法有3种，第一种是。

ω-Jacobi迭代法，第二种是Gauss-Seidel迭代法，第三种是SOR法，其迭代公式分别如下：1）ω-Jacobi
2）Gauss-Seidel
3）SOR
在每个松弛迭代中，SOR法的计算量与ω-Jacobi迭代法与Gauss-seidel迭代法相当。

如果选择好最佳松弛因子，该方法的收敛速度比ω-Jacobi代法与Gauss-seidsl迭代法高一个量级。

（2）无权多级网络法
高斯一赛德尔松弛算法的收敛速度很慢，为了提高运算效率，可采用多级格网技术。

高斯一赛德尔松弛算法是一种典型的局域平滑算子，它可以迅速的去除信号中的高频成分，但对低频成分的滤除速度
却非常慢，因此导致它很难收敛。

多级格网方法的核心思想就是将低频成分转化为高频成分，加速高斯一赛德尔松弛算法的处理速度。

这种由低频到高频的转换是通过格网重采样实现的，粗格网的低采样率提高了残差点误差的空间频率。

如图的网格金字塔，每一级格网分辨率都是其前一级的二分之一，这样原始格网的低空间频率会随着格网逐渐变粗而迅速提高。

该方法计算速度大大提高，所得结果和高斯一赛德尔松弛算法基本相同。

无权多级格网法解缠速度快，并且不需要质量图。

但是由于解缠时不是绕过残差点，而是穿过残差点进行积分运算，这样解缠相位就可能被这些残差点打断，这种现象在残差点附近最为明显，而且其影响可以传播到很远的地方，这就造成了最小二乘解处过低估计了相位梯度值，使相位梯度差无法完全拟合。

加权最小二乘法
不带权的最小二乘法运算速度快，但由于解缠时穿过了相位的不一致区，所以常常造成解缠面的误差。

这个缺陷可以通过引入权重来弥补。

一般来说，权重数据j的取值范围是[0，1]，权重由相位质量图或其他先验知识所确定。

加权重的最小二乘即是求解下式:
网络流法
相位解缠的前两类算法(路径跟踪法和最小范数法)都致力于克服相位场的不一致性，其不同之处在于每种算法克服该困难的方法不一样，由此也导致了各种算法的速度和精确性各异。

另外，对于噪声比较严重的干涉图，目前的算法都未很好的解决。

网络流法同时兼顾
速度和精确性，很好的解决了这两大方面问题。

基于网络流相位解缠算法的主要思想是将解缠相位的导数与缠绕相位的导数之间的差异最小化。

由于通常最小化问题的计算效率比较低，因此算法将其转化为求解最小费用流的网络优化问题，因而大大降低相位解缠算法的时间和空间复杂度。

最小费用流问题的解决有比较成熟的算法，并且计算效率比较高，提高了整个解缠算法的效率。

该方法还可以将整个误差限制在一个小范围内，防止误差的再传递，解缠结果较精确。

网络流法按照所处理的网络不同，又可分为基于规则网络的最小费用流算法和基于不规则网络的最小费用流算法。

基于不规则网络的最小费用流算法的主要思想是:首先提取出相干系数较高的相位，作为高质量的相位数据集合。

然后根据这些相位点的位置建立一个Delaunay三角网。

然后识别出该三角网中的残差点，应用最小费用流算法连接正负残差点对，建立枝切线。

最后，按照绕过枝切线或者穿过枝切线加减2nπ的方法进行积分从而得到解缠相位值。