《测绘通报》测绘科学前沿技术论坛论文集 1 InSAR图像相位解缠的最小费用流法
及其改进算法研究1
蒋廷臣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－），男，汉族，河南商丘人，副教授，主要从事主要从事精密工程测量和测绘教育研
究。
《测绘通报》测绘科学前沿技术论坛论文集 2 提出了一个统一框架，将经典的相位解缠算法进行了理论上的分析，指明了这些算法内在的数学联系。在这个框架下，相位解缠被认为是一个优化问题，存在着一个目标函数和一个比例因子，解缠的目的就是求解目标函数的最小值。Ghiglia和Romero提出的LP范数目标函数方程如下：
?minimize???i(,xj)??ix,j???ix,j?i,jp?yp???i(,yj)??iy???? （1） ,ji,ji,j?方程（1）对大括号里的函数求最小值，式中0?p?2，??(x)和??(x)分别为x方向的解缠相位梯度和缠绕相位梯度；??(y)和??(y)分别为y方向的解缠相位梯度和缠绕相位梯度；?为对应于每一个梯度差所定义的权；i和j分别代表行数和列数；有的行（i）和列（j）求和。
????为求和算子，表示对所
i,j?指距离向和方位向的离散差分，即说明相位差是以从－?到??的间隔进行缠绕的。最
小费用流法是当p=1时，基于最优估计的解缠方法。 2.2最小费用流法原理
网络流算法最早见于Costantini等提出的基于网络流的相位解缠方法，这种方法是将相位解缠问题转化为最小化问题，通过在全局范围内搜索路径和最短枝切来求得最小化问题的最优解。该方法可以应用于规则网络(网格)，也可以用于不规则网络(三角网)，其原理如下：
在一个M*N大小的方格网内，设?和?分别表示解缠和未解缠的相位，则有：
?(i,j)??(i,j)?2?n （2）
式中，n为整数且?(i,j)????,??，相位解缠过程就是从?(i,j)到?(i,j)的过程。
定义相邻像素点间的差分估计：
??1(i,j)??(i?1,j)??(i,j)?2?n1 （3） ??2(i,j)??(i,j?1)??(i,j)?2?n2 （4）
式中，n1(i,j)为基于先验知识选取，使??l(i,j)????,??（l=1,2）成立的整数值。由于积分路径的不同，??l(i,j)????,??（l=1,2）并不能和相邻点的差分保持一致，因而定义以下差分的残差：
?k1(i,j)???k(i,j)????2?12???(i?1,j)????1(i,j)????(i,j?1)???????(i,j)?? （5） ???2?
k1(i,j)和k2(i,j)是很小的数，可以用如下的最小化问题来估算残差k1(i,j)和k2(i,j)：
?min?c(i,j)k(i,j)?c(i,j)k(i,j)??1? （6） ?122?k1,k2??i,ji,j?根据网络流理论，式（4.20）这个最小化问题可以转化为求解网络中的最小费用流来解决，最小费用流问题的输入为各个结点的度(即各残差的值)和每条流的费用,而该问题的输出为各
《测绘通报》测绘科学前沿技术论坛论文集 3 条流的流量,并且费用和最小。在计算出k1(i,j)和k2(i,j)后，再计算相位梯度，最终可以通过（7）式计算出解缠后的相位，得出最终的解缠结果。
?(i,j)??(0,0)????(p?1,0)??(p,0)???(?(i,q?1)??(i,q)) （7）
p?0q?0i?1j?12.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）基本参数信息
干涉图参数 轨道号 9192 6687 基线 (m) 496 垂直基线 入射角 (m) 475 (o) 22.8 高程模糊度 (m) 17.1 时间跨度 (天) 175 多视处理 2：10 《测绘通报》测绘科学前沿技术论坛论文集 4
图1 MCF解缠中的相关图和解缠图
四、改进的最小费用流解缠方法
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 被称为瓦片，它们都处于由行程线和等质区的边缘 （或图廓） 包围的区域。
《测绘通报》测绘科学前沿技术论坛论文集 5
图2 分块算法原理图
根据分块算法原理，对3节中的算例进行了解缠，其解缠效率见表2，由表2可以看出，
利用分块算法大大提高了解缠速度，降低了对计算机的要求。虽然如此，但是利用分块算法进行解缠时，解缠结果在块与块之间的条纹无法连接，如图1，为了得到更好的解缠结果，需要进行融合处理。
表2 分块解缠效率表
Patchs 时间（m） 图像参数 1 大小：343M，5167ｘ3400 大小：225M 3400ｘ3400 1 2 3 4 5 6 7 9 12 16 Failure Failure Failure 16.59 27.2 8.1 7.5 8.1 8.2 9.3 10.6 13.4 5.4 5.3 6 6.4 7 7.6 8.9 10.8 2
4.2基于Contourlet变换的超小波融合方法
由于小波算法无法识别自然图像中固有的线奇异和面奇异 , 同时其捕获的方向性信息也受到限制。2001年, Do和 Vetterli 构建了针对高维信号处理的多尺度几何分析工具 Contourlet, Contourlet继承了小波的多分辨和时频局部化特性,对图像具有更加稀疏的表达能力。同时，Contourlet变换能比小波算法更好地表示稀疏的二维信号, 具有良好的多分辨率、局部化、方向性和各向异性，更适于处理具有超平面奇异性的图像信号。将 Contourlet变换引入干涉图像融合, 可以利用其优良特性更好地提取解缠图像中的几何特征, 从而获得较好的解缠结果。 Contourlet变换采用双重滤波器组结构，首先采用拉普拉斯塔式分解对相位图进行多尺度分解以捕获点奇异，对分解后的低频子带继续使用 LP变换进行迭代分解，便可以将原始图像分解为一系列不同尺度上的低频和高频子带，随后，对 LP分解所得到的高频子带使用方向滤波器组（DFB）进行方向性分析。DFB的作用是捕获图像的方向性高频信息，将分布在同方向上的奇异点合成为一个系数。在计算时采用 l层的树结构分解，在每层先通过扇型滤