多基线InSAR图割相位解缠算法研究*

张斌，胡庆荣，李爽，韦立登

(1.北京无线电测量研究所，北京 100854;2.中国航天科工集团有限公司 第二研究院，北京 100854)

0 引言

InSAR技术通过利用2幅或多幅SAR图像的相位信息可以获取大范围、高精度的地表三维信息。其中相位解缠是InSAR处理的关键步骤，其处理精度对DEM(digital elevation model)产生重要影响[1-2]。传统的单基线相位解缠，大都假设Itoh条件满足的前提下进行的，其在相位噪声大和陡峭地形情况下很难满足高精度解缠要求，甚至失效，无法准确还原地形[3]。而多基线相位解缠算法能打破Itoh条件限制，获取更多的观测区的信息，融合多幅干涉图，减小干涉相位的模糊区间，降低解缠难度[4]。

近年来，多基线InSAR相位解缠得到了长足发展，Ghiglia[5]等提出多基线相位解缠绕方法的最小二乘法和最大似然函数法，Ferraiuolo[6]等提出基于马尔可夫随机场的最大后验概率估计方法，夏香根[7]等提出一种基于中国余数定理方法相位解缠，于瀚雯[8]提出基于聚类分析的多基线相位解缠绕算法。而在多基线最大后验概率估计方法中，Ferraioli[9]等提出了基于图割的多基线相位解缠，其直接由多幅干涉条纹图直接计算高程值，从而避开了相位解缠过程。但该类方法由于采用Ishikawa图割[10]，其所建立的高度量化值越多，其相位解缠所用的内存和处理所花的时间越大，因此A.shabou[11]等提出了一种基于二分法的多基线相位解缠方法，其所占用的内存减少，但其每一个循环都要访问所有的高度量化值，相应的计算时间却增加。相比直接估计图像的高程值，用Ishikawa方法估计相位解缠所需的缠绕数(2π的整数倍)所花的内存和处理时间无疑更占优势，而为了进一步节省内存和减少处理时间，本文据此提出了基于迭代的多基线图割方法，在不影响全局最优的基础上通过对图像每个点所需的缠绕数进行分段，通过优化迭代，使其最终达到最优的缠绕数。

1 图割模型

干涉复图像可以表示为

φ=Aejφ+n，

(1)

式中：φ为真实干涉相位；n为零均值的复圆周高斯噪声；φ为带有噪声的干涉相位。

根据贝叶斯准则，基于最大后验概率的相位解缠模型可表示为

(2)

式中：P(φ|φ)为在实际干涉相位下真实相位的概率;P(φ)为实际相位的概率，此处为一常数;P(φ|φ)为假设已知真实相位的情况下实际相位的概率，即数据的似然函数;g(φ)为真实干涉相位的先验概率密度函数。

假设所有的像素点之间是相互独立的，可知

(3)

式中：i为相应的通道(基线)，此处共有C个基线;P为每个通道内的元素，假设数据为M×N个元素；γi,p为第i个通道内第p个像素的相干系数值，φi,p为对应的实际相位值。

在一阶马尔可夫(MRF)模型中，相位的先验概率密度函数可表示为

(4)

对式(2)取对数后，基于MRF模型的图割相位解缠算法可表示为

(5)

式中：

2 基于迭代的图割相位解缠算法

2.1 算法流程

基于式(5)的图割相位解缠算法主要分为2种：第1种为Boykov等[12]提出的α-β移动迭代接近最优值，算法处理速度快，但在噪声比较大情况下或者欠采样(地势陡峭地区等)其整体最优的鲁棒性下降，局部解缠性能不佳;第2种方法为全局算法[9，13-14]，该算法将每个点可能所需的缠绕数作为标签集合，然后运用图割算法在标签集合里寻找到整体最优值，从而进行相位解缠，此方法为整体最优算法，鲁棒性较高，但其占用内存大，计算时间也很大，而A.Shabou[11]针对此问题提出了二分法的相位解缠来解决此问题，此算法虽然占用内存变小但其运算时间大大增加。本文针对Ishikawa算法占用内存大与计算时间长的缺点，提出一种基于迭代的多基线图割算法，算法在不影响整体最优的情况下，将每个点可能所需的缠绕数进行初始分段，然后利用图割算法，不断逼近整体最优值。

算法的流程如下：

(1) 估计整幅图的缠绕数，选取合适的分段(分段间隔为2π)。

(2) 在分段内，更新相位值，构建改进的Ishikawa图网络模型。

2.2 图网络的建立

图网络是图割算法的基础，其将每一条边赋予相应权值并进行计算，利用网络图的最大流最小割来寻求当前迭代中的最优值。

图1单向箭头表示数据边或约束边，双向箭头表示交互边，p,q,r表示数据的元素(p,q为相邻的元素)，割集如图中Cut粗体线段所示。

图1中的顶点集合为干涉图中像素点的数目再加上S,T2个节点，即

i=1,2,…,L}∪{s}∪{t}，

(6)

式中：L为算法中所选取的分段的数目;M为干涉图的数；N为干涉图的列。

图1的有向边集可分为4个部分，分别为数据边εD，约束边εC，交互边εI，补充边εA。

数据边如图1中向上的箭头表示，其代表数据的似然函数项，即

(7)

(8)

数据边的权值为

(9)

(10)

约束边如图1中向下箭头表示是为了在进行最大流最小割时，每个像素有且仅有一个最优值，此时有

(11)

(12)

约束边的权值表示为

(13)

交互边如图1中双向箭头表示，体现能量函数的先验信息，此时有

(14)

(15)

在交互边中为了达到公式(5)的平衡，需要加入补充边，此时有

(16)

(17)

在此可以将εA对应的权值cp,q(i)和cp,q(j)加到εD对应边的权值去。

通过构建网络图模型将能量最小问题转化为求解网络图的最小割，如图1中粗体Cut所示，得到图的最大流最小割的同时也就得到了当前迭代中所有像素点的缠绕数，再依据算法流程向下执行。

3 实验校验

本文提出的方法在占用内存和计算时间方面均有一定程度的改进，为了对上述算法进行验证，本节通过仿真中进行了说明。仿真实验在Matlab上运行，其中图割算法中的最大流最小割算法用C++实现[15]。

用三基线进行仿真实验，其中长基线原始高斯面如图2所示,像素大小为64×64,对原始相位图图2加入相干系数为0.9的高斯白噪声。图3为加入高斯白噪声后对应的缠绕图，从缠绕图中可以看出该高斯面出现了较大面积的欠采样情况，即打破了Itoh条件，此时用传统的单基线相位解缠方法难以得到较好的结果，因此适合采用多基线相位解缠方法。相应的短基线缠绕相位图如图4所示，长基线与短基线的比例为7：3，本次仿真实验与文献[6]中的Ishikawa方法进行了对比，实验结果如图5，6所示。

对相位解缠前后的误差(与原始相位之差)进行分析比较，则有error迭代=5.922 8 rad,errorIshikawa=-232.477 9 rad，迭代算法的误差较小，更为接近原始相位，这与迭代算法分段后能更好地兼顾局部最优有关系。从仿真可以看出迭代算法能打破Itoh条件限制，较好地处理相位变化较快的区域，且运算时间相比Ishikawa算法有所减少。

4 结束语

本文提出了一种基于迭代的多基线InSAR图割相位解缠算法。该算法针对Ishikawa图割方法处理相位解缠时占用内存大、计算时间长的问题，通过对缠绕数进行分段，采取迭代的方式逐步逼近最优值。算法通过构建能量函数模型，利用一阶马尔可夫模型的先验信息，构建改进的Ishikawa图网络模型，通过最大流最小割算法及最小化能量函数准则逐步优化相位值，从而完成相位解缠。仿真实验验证该迭代算法能够减少相位解缠所需的内存和处理时间。实验是在低噪情况下进行，对于噪声比较大的区域，其解缠的鲁棒性有待进一步验证。此外，对实测数据的处理也将在下一步展开。

参考文献：

[1] BAMLER R,HARTL P. Synthetic Aperture Radar Interferometry[J].Inverse Problems,1998,14(4):1-54.

[2] 曲长文,侯海平,周强.合成孔径雷达三维成像技术及特点评述[J].现代防御技术,2010,38(4):110-115.

QU Chang-wen,HOU Hai-ping,ZHOU Qiang.An Overview of Three-Dimensional Imaging Technology and Technical Characteristic of SAR[J].Modern Defence Technology,2010,38(4):110-115.

[3] LOTH K.Analysis of the Phase Unwrapping Problem[J].Applied Optics,1982 ,21(14):2470.

[4] FERRAIUOLO G,PASCAZIO V,SCHIRINZI G.Maximum a Posteriori Estimation of Height Profiles in InSAR Imaging[J].IEEE Geoscience and Remote Sensing Letters,2004,1(2):66-70.

[5] GHIGLIA D C,WAHL D E.Interferometric Synthetic Aperture Radar Terrain Elevation Mapping from Multiple Observations[C]∥Proceeding of IEEE Digital Signal Processing Workshop IEEE Press,1994:33-36.

[6] FERRAIUOLO G,MEGLIO F,vito Pascazio,et al.DEM Reconstruction Accuracy in Multichannel SAR Interferometry[J].IEEE Transaction on Geoscience and Remote Sensing,2009,47(1):191-201.

[7] XIA Xiang-gen,WANG Gen-yuan.Phase Unwrapping and a Robust Chinese Remainder Theorem[J].IEEE Signal Processing Letters,2007,14(4):247-250.

[8] YU Han-wen,LI Zhen-fang,BAO Zheng.A Cluster-Analysis-Based Efficient Multibaseline Phase Unwrapping Algorithm[J].IEEE Transaction Geoscience and Remote Sensing,2011,49(1):478-487.

[9] FERRAIOLI G,SHABOU A,TUPIN F,et al.Multichannel Phase Unwrapping with Graph Cuts[J].IEEE Geoscience and Remote Sensing Letters,2009,6(3):562-566.

[10] ISHIKAWA H.Exact Optimization for Markov Random Fields with Convex Priors[J].IEEE Trans Pattern Anal Mach Intell,2003,25(10):1333-1336.

[11] SHABOU A,DARBON J,TUPIN F.Multilabel Partition Moves for MRF Optimization[J].Image and Vision Computing,2013,31(1):14-30.

[12] BOYKOV Y,VEKSLER O,ZABIH R.Fast Approximate Energy Minimization via Graph Cuts[J].IEEE Trans.Pattern Anal Mach Intell,2001,23(11):1222-1239.

[13] JOSÉ M,Dias Bioucas,Valãdao Gonçalo.Phase Image :Unwrapping and Denoising with Diversity and Multi-Resolution[J].IEEE Transaction on Geoscience and Remote Sensing,2001,33(3):579-589.

[14] Jerome Darbon,Marc Sigelle.Image Restoration with Discrete Constrained Total Variation Part II:Levelable Functions,Convex Priors and Non-Convex Cases[J].Journal of Mathematical Imaging and Vision,2006:277-291.

[15] BOYKOV Y,Vladimir Kolmogorov.An Emperiment Comparison of Min-Cut/Max-Flow Algorithms for Energy Minimization in Vision[J].IEEE Trans.Pattern Anal Mach Intell,2004,26(9):1124-1137.