InSAR图像相位解缠的最小费用流法及其改进算法研究

InSAR图像相位解缠的最小费用流法
及其改进算法研究
蒋廷臣1，2，焦明连1，史建青1，王秀萍1
（1．淮海工学院测绘工程学院,江苏连云港222001；
2．武汉大学卫星导航定位技术研究中心，武汉430079）
摘要：最小费用流法是基于网络流的相位解缠方法，解决了许多解缠方法无法消除相位噪声对高相干区域影响的问题，在此基础上，本文针对该方法解缠时速度较慢和对计算机性能要求较高的缺点而提出改进算法，即将干涉图像分为若干子区域分别进行处
理，再利用基于Contourlet变换的超小波方法进行融合处理，最后用算例进行了验证，结果表明最小费用流法及其改进算法是一个
较好的解缠方法。

关键词：干涉测量相位解缠最小费用流法分块算法小波融合
一、前言
随着测绘新技术新理论的发展，现代大地测量范畴得到了较大拓宽，现在，合成孔径雷达干涉测量(Interferometry Synthetic Aperture Radar—InSAR)已成为其分支学科。

合成孔径雷达干涉测量 ( InSAR)利用合成孔径雷达数据的相位信息提取地面三维信息，主要用于测量地面的高程和监测其变形。

随着COSMOS和terraSAR卫星的发射成功，该技术日益受到各国政府部门以及科学工作者的重视。

在InSAR数据处理过程中，相位解缠是合成孔径雷达干涉测量的关键流程，它的准确性直接影响到 InSAR生成的数字高程模型的精确性。

现在所有的解缠方法都是基于这样的假设，即
φ差的绝对值小于π。

解缠后的真实相位是平滑且变化缓慢，同时图像各相邻像素的干涉相位
但是，雷达阴影、去相关等因素引起的噪声和伪信号往往造成相位数据不连续，给相位解缠带来极大的困难，目前大部分算法都无法圆满地解决这些问题 ,解缠的结果常常会有较大的误差，由此得到的数字高程模型就会与实际情况存在较大的差别。

如何能够从质量较差的数据当中提取有用的信息，而忽略噪声对解缠过程的影响，成为一个急待解决的问题。

基于上述，本文根据统一的解缠数学模型和网络优化原理，阐述了最小费用流法法的相位解缠方法，并针对该方法解缠时速度较慢而提出分块算法，将整幅图像分为若干子区域分别进行处理 ,再利用超小波方法进行融合处理，从而得到较理想的解缠效果，同时利用算例进行了比较分析，较好地解决了上述问题。

二、最小费用流法解缠原理
2.1统一解缠模型
经过多年对相位解缠方法的研究，现在已有很多的解缠方法。

在1996年，Ghiglia和Romero
第一作者简介：蒋廷臣（1975-），男，汉族，四川蓬安人，武汉大学测绘学院博士生，主要从事GPS与宽幅SAR融合的相关理论与方法研究。

第二作者简介：焦明连（1963－），男，汉族，河南商丘人，副教授，主要从事主要从事精密工程测量和测绘教育研究。

提出了一个统一框架，将经典的相位解缠算法进行了理论上的分析，指明了这些算法内在的数学联系。

在这个框架下，相位解缠被认为是一个优化问题，存在着一个目标函数和一个比例因子，解缠的目的就是求解目标函数的最小值。

Ghiglia 和Romero 提出的L P 范数目标函数方程如下：
⎭
⎬⎫⎩⎨⎧∆-∆+∆-∆∑∑j i j i p y j
i y j i y j i p x j
i x j i x j i imize ,,,,)(,,,)(,min φϕϖφϕϖ （1） 方程（1）对大括号里的函数求最小值，式中20≤≤p ，)(x ϕ∆和)(x φ∆分别为x 方向的解缠相位梯度和缠绕相位梯度；)(y ϕ∆和)(y φ∆分别为y 方向的解缠相位梯度和缠绕相位梯度；ω为对应于每一个梯度差所定义的权；i 和j 分别代表行数和列数；
(),i j
•∑为求和算子，表示对所有的行（i ）和列（j ）求和。

∆指距离向和方位向的离散差分，即说明相位差是以从－π到π+的间隔进行缠绕的。

最小费用流法是当p=1时，基于最优估计的解缠方法。

2.2最小费用流法原理
网络流算法最早见于Costantini 等提出的基于网络流的相位解缠方法，这种方法是将相位解缠问题转化为最小化问题，通过在全局范围内搜索路径和最短枝切来求得最小化问题的最优解。

该方法可以应用于规则网络(网格)，也可以用于不规则网络(三角网)，其原理如下：
在一个M*N 大小的方格网内，设ϕ和φ分别表示解缠和未解缠的相位，则有：
n j i j i πϕφ2),(),(+= （2）
式中，n 为整数且[]ππϕ,),(-∈j i ，相位解缠过程就是从),(j i φ到),(j i ϕ的过程。

定义相邻像素点间的差分估计：
112),(),1(),(n j i j i j i πφφφ+-+=∆ （3）
222),()1,(),(n j i j i j i πφφφ+-+=∆ （4）
式中，),(1j i n 为基于先验知识选取，使[)ππφ,),(-∈∆j i l （l=1,2）成立的整数值。

由于积分路径的不同，[)ππφ,),(-∈∆j i l （l=1,2）并不能和相邻点的差分保持一致，因而定义以下差分的残差：
⎪⎪⎭
⎫ ⎝⎛∆∆-⎪⎪⎭⎫ ⎝⎛++=⎪⎪⎭⎫ ⎝⎛),(),()1,(),1(),(),(2121
21j i j i j i j i j i k j i k ϕφϕϕπ （5）
),(1j i k 和),(2j i k 是很小的数，可以用如下的最小化问题来估算残差),(1j i k 和),(2j i k ： {}⎭
⎬⎫⎩⎨⎧+∑∑),(),(),(),(,min 2,21,121j i k j i c j i k j i c k k j i j i （6） 根据网络流理论，式（4.20）这个最小化问题可以转化为求解网络中的最小费用流来解决，最小费用流问题的输入为各个结点的度(即各残差的值)和每条流的费用,而该问题的输出为各
条流的流量,并且费用和最小。

在计算出),(1j i k 和),(2j i k 后，再计算相位梯度，最终可以通过
（7）式计算出解缠后的相位，得出最终的解缠结果。

()∑∑-=-=ψ-+ψ+ψ-+ψ+Φ=Φ101
0)),()1,(()0,()0,1()0,0(),(i p j q q i q i p p j i （7）
2.3最小费用流法解缠步骤
通过上面对最小费用流解缠方法的原理描述，现给出最小费用流法的解缠步骤：
(1)为相干图选取一个阈值，以该阈值作为评价相位质量的标准，提取出所有相干系数高于该阈值的相位，并将这些相位添加到一个集合。

(2)在相位集中建立一个Delaunay 三角网，根据三角网构造其对偶图，在原三角网中求得残差，并将这些残差映射到其对偶网络中。

(3)在对偶网络中，应用最小费用流方法连接各正负残差点对，计算出最小费用流集合。

(4)根据流的大小和方向对相位矩阵进行积分，得到解缠结果，再从高质量区域向低质量区域进行积分，得到全局结果。

从本质来讲，此种方法仍属于路径跟踪法，但是又有所不同，由于能够有效区分高质量数据和低质量数据，因而在解缠过程中就可以避开质量不好的区域，提高解缠精度。

三、最小费用流法解缠算例
2003年12月26日世界时01:56:52，伊朗东南部的Bam 地区发生了强度为6.6级的剧烈地震，震中位置位于北纬29.01°，东经58.26°，造成了巨大的人员伤亡和财产损失，蔚为壮观的Bam 古城遗址也毁于一旦。

此次地震的发生是由于阿拉伯板块北移与欧亚板块挤压造成的应力造成。

我们选择轨道号为6687和9192两幅ASAR 图像，其干涉图参数信息见表。

在经过轨道基线估算，影像配准和产生干涉图后，利用最小费用流法对以6687和9192的干涉图进行了相位解缠，该算法将干涉相干作为评价相位质量的标准 ,除去数据质量差的相位，见图1，其中图1左边是相位解缠时所用到的有效屏蔽图，图1右是最小费用流法的解缠结果图。

由前述可知，最小费用流法相比于有许多优点，解决了现在一些常用解缠方法不能解决的问题，如噪声对解缠结果的影响，但是，最小费用流法也有一个缺点，即解缠时占用内存大，耗时多。

故，为了快速解缠，本文提出改进的最小费用流法，其主要由分块算法和图像融合两部分组成。

表1干涉图（9192-6687）基本参数信息
图1MCF解缠中的相关图和解缠图
四、改进的最小费用流解缠方法
4.1分块算法
瓦片分割是指将干涉图分割成为不重叠的矩形块，这些矩形块可以进行单独进行相位解缠，就像他们是完全独立的图像。

分割的目的是为了降低相位解缠过程中所需的内存资源，同时能在一定程度上加快解缠速度。

如图2 所示，图中有l1和l2两根直线，i、j、k、m、n共5个区域。

现在的问题是在图像上寻找l1和l2，进而得到i、j、k、m、n等区域。

在此，直线l1和l2被称为行程线，是指平行于x轴（或y轴）的线段，它至少应相切某一等质区的局部最大顶点或最小顶点（通常指沿x轴或y轴方向），其起点可以是图廓或等质区的边缘。

区域i、j、k、m、n 被称为瓦片，它们都处于由行程线和等质区的边缘（或图廓）包围的区域。

图2分块算法原理图
根据分块算法原理，对3节中的算例进行了解缠，其解缠效率见表2，由表2可以看出，利用分块算法大大提高了解缠速度，降低了对计算机的要求。

虽然如此，但是利用分块算法进行解缠时，解缠结果在块与块之间的条纹无法连接，如图1，为了得到更好的解缠结果，需要进行融合处理。

表2分块解缠效率表
Patchs
时间（m）
图像参数
1 2 3 4 5 6 7 9 12 16
1 大小：343M，
5167ｘ3400
Failure Failure 27.28.17.58.18.29.310.613.4
2 大小：225M
3400ｘ3400
Failure 16.59 5.4 5.36 6.477.68.910.8
4.2基于Contourlet变换的超小波融合方法
由于小波算法无法识别自然图像中固有的线奇异和面奇异, 同时其捕获的方向性信息也受到限制。

2001年, Do和Vetterli 构建了针对高维信号处理的多尺度几何分析工具Contourlet, Contourlet继承了小波的多分辨和时频局部化特性,对图像具有更加稀疏的表达能力。

同时，Contourlet变换能比小波算法更好地表示稀疏的二维信号, 具有良好的多分辨率、局部化、方向性和各向异性，更适于处理具有超平面奇异性的图像信号。

将Contourlet变换引入干涉图像融合, 可以利用其优良特性更好地提取解缠图像中的几何特征, 从而获得较好的解缠结果。

Contourlet变换采用双重滤波器组结构，首先采用拉普拉斯塔式分解对相位图进行多尺度分解以捕获点奇异，对分解后的低频子带继续使用LP变换进行迭代分解，便可以将原始图像分解为一系列不同尺度上的低频和高频子带，随后，对LP分解所得到的高频子带使用方向滤波器组（DFB）进行方向性分析。

DFB的作用是捕获图像的方向性高频信息，将分布在同方向上的奇异点合成为一个系数。

在计算时采用l层的树结构分解，在每层先通过扇型滤
波器组 (QFB) 进行扇型方向上的频率切分，随后与旋转重采样操作适当组合以实现图像高频信息方向性分析，捕获图像中的线、面奇异性，Contourlet 变换的原理框图如图 3所示。

图3 Contourlet 原理变换图
Fig 3 Principle of contourlet transform 本文基于区域特性的Contourle 变换，将其应用于最小费用流法解缠后的相位图融合, 算法的实现步骤如下：
（1）Contourlet 分解；将精确配准的相位图像进行Contourlet 变换,得到相应的Contourlet 系数集合，在融合中一般采用 3层或 4层分解。

设两相位图像为 I A , I B , 使用m 级 LP 分解 , 每个尺度上的方向数分别为 l( j) , 则相位图像的 Contourlet 分解过程为
{}B
A n n b n b n b n c
B d B c B c B c B c j i I A d A c A c A c A c j i I j l m m m m m m m B m m m A ,,)(,),(),()())(),(),(,),(),((),())
(),(),(,),(),((),()(,2,1,121121==→→-- （8）
式中m d 是低频子带；m c 为尺度m 上的方向子带集合；k m b ,为尺度m 第k 个方向的高频子带。

(2) 图像融合；对于分解后的低频子带和所有高频子带，使用基于区域特性的融合规则进行判别和融合处理，得到各尺度上融合后的 Contourlet 系数。

(3) Contourlet 重构；重构是分解的逆过程，对融合后的 Contourlet 系数进行Contourlet 逆变换，得到重构的融合图像，该图像包含原有多幅图像中的信息.，则为融合后的结果，此过程可以表示为
),(),,,,,(121j i I d c c c c AB F m F m F m F F →- （9）
根据上面的图像融合步骤，将分块算法得到的不同解缠结果图（本文使用的是块数6和
LP
12的相位结果图）进行融合处理，经过处理后得到了融合相位图，见图4，从图4可以看出，融合结果显得光滑，较仅利用分块算法的结果有了很大改善。

通过这样的处理，不仅未影响解缠的效果，而且降低对计算机性能的要求，提高了解缠的速度。

图4融合后的相位效果图
五、结论
相位解缠是合成孔径雷达干涉测量的重要步骤，但相位解缠中，噪声容易对高相干区域的解缠结果造成影响，从而影响最终的解缠结果，也是现在很多解缠方法所无法解决的问题。

为此，本文针对最小费用流法的原理，提出了最小费用流法的改进算法，然后利用InSAR数据对其进行了验证，结果表明最小费用流法算法可以忽略含有大量噪声的区域,而直接处理数据质量好的区域,避免了解缠过程中噪声区域的误差对高质量区域的影响，而其改进算法则解决了最小费用流法的运算速度慢，对计算机性能要求较高的问题。

参考文献：
[1]Alessandro Ferretti, Andrea Monti-Guarnieri, Claudio Prati, Fabio Rocca.2007.InSAR Principles: Guidelines for SAR Interferometry Processing and Interpretation,The Netherlands : ESA Publications [2]AnaBertran Ortiz.2007. SCANSAR-TO-STRIPMAP INTERFEROMETRIC OBSERVATIONS OF HAWAII[D],USA:STANFORD UNIVERSITY
[3]Chen C.W. 2001.Statistical-cost Network-flow approaches to Two-demensional phase unwrapping for radar interferometry[D]. USA: Stanford university
[4]Chipman J. W. 2001. Topographic mapping of forested areas in the western Great Lakes region using spaceborne synthetic aperture radar interferometry[D]. Madison:The University of Wisconsin
[5]Davis J. L., Herring T A, Shapiro II et al,1985. Geodesy by radio interferometry: Effects of atmospheric modeling errors on estimates of baseline length. Radio Science, Vol.20,1593~1607
[6]Dennis C. Ghiglia, Mark D. Pritt , 1998, Two-dimensional Phase Unwrapping: Theory, Algorithms, and Software . New York: Published by John Wiley and Sons
[7]YAN Jing-wen, ZHUO L in,QU Xiao-bo，A New Research on Image Fusion, Journal of Xiamen University of Technology, Vol. 15 No. 4 2007:44－50
[8]SHAO Yong - she, LI J ing.Image fusion method of r estr aining speckle noise of SAR image .Computer Engineering and Applications, 2008, 44( 6) : 82- 84.
[9]YU Yong ,WANG Chao ,ZHANG Hong ,LIU Zhi ,GAO Xin. A Phase Unwrapping Method Based on Network Flow Algorithm in Irregular Network. JOURNAL OF REMOTE SENSING, 2003, Vol. 7 , No. 6:472-478
[10]CHEN J ia-feng , CHEN Hai-qing , REN HaiOxia. Least squares phase unwrapping algorithm based on the wavelet transform. OPTICAL TECHNIQUE,2007, Vol. 33 No. 4:613-618。