网络流与曲面拟合结合的相位解缠方法

， ，

( 1.广西无线宽带通信与信号处理重点实验室， 广西桂林 541004； 2.桂林电子科技大学计算机与信息安全学院， 广西桂林 541004； 3.桂林电子科技大学信息与通信学院， 广西桂林 541004)

0 引言

相位解缠是干涉合成孔径雷达(Interference Synthetic Aperture Radar，InSAR)干涉测量的关键流程[1-3]，它的准确程度直接影响到InSAR生成的数字高程模型的精确性[3]。干涉图的相位展开方法有多种，主要分为4类：第一类是路径跟踪法[4-5]，第二类是最小范数法[6-7]，第三类是网络流法[8]，第四类是其他方法。

路径跟踪法最常用的是枝切法。利用枝切法进行相位解缠是Goldstein等首先提出的[4]。枝切法在相干性差、残差点较多的区域难以选择准确的枝切线，容易形成一个个被枝切线包围的孤立区间，即“解缠孤岛”，积分路径都无法到达从而无法解缠。“解缠孤岛”由于没有明显的联系,单独解缠时会造成很大误差。最小范数法最典型的就是最小二乘法，它是一种全局算法。解缠结果中不存在不能进行相位展开的区域，但由于此算法穿过而非绕过残差点，很容易使相位误差传递到整个区域。这就造成了整体上结果为连续的，但实际上没有一点相位的值是准确的[9]。网络流法是通过全局范围搜索最优路径和最短枝切线来求得最小化问题的最优解。在网络流成熟的条件下网络流算法可大大提高运算的效率[10]，与此同时保证一定的精度，在误差的限制方面有着良好的效果，可以防止误差传递到整个区域内，对于噪声较为严重的干涉图也可以进行解缠并取得较为精确的结果[11-13]。但是在噪声较大的环境下网络流法的解缠结果会将噪声不可避免地沿着积分路径传递，造成解缠结果的严重“尖峰毛刺”现象，使得解缠结果不理想。

结合残差点的思想、曲面拟合法以及网络流法提出了一种新的相位解缠算法。新算法根据相关度找出解缠结果的噪声点，然后运用最小曲面拟合法重构值代替噪声点的数值，最后沿着积分路径更正有误差传递的结果。逐一处理所有噪声点，最终的解缠结果可以看作是由多个最小曲面多次贴合的结果。本文首先介绍相关理论，然后介绍新算法的具体步骤和推导，最后通过实验验证可行性。

1 新算法描述

网络流算法最早见于Costantini等提出的基于网络流的相位解缠方法[8],从根本上来讲此种方法仍属于路径跟踪法。这种方法是将相位解缠问题转化为最小化问题，通过在全局范围内搜索路径和最短枝切线来求得最小化问题的最优解。在一个M×N大小的方格网内,设φ和φ分别表示解缠和未解缠的相位,则有

φ(i,j) =φ(i,j) + 2πn

(1)

式中，n为整数，φ(i,j)∈[-π，π]，相位解缠就是从φ(i,j)到φ(i,j)的过程。根据二维行为解缠原理，定义相邻像素点间的差分估计[14]：

Δφ1(i,j)=φ(i+1,j)-φ(i,j)+2πn1

(2)

Δφ2(i,j)=φ(i,j+1)-φ(i,j)+2πn2

(3)

式中，n1(i,j)，n2(i,j) 为基于先验知识选取，使Δφl(i,j)∈[-π,π)(l=1,2) 成立的整数值。由于积分路径的不同，Δφl(i,j)∈[-π,π)(l=1,2) 并不能和相邻点的差分保持一致，因而定义以下差分残差：

由于k1(i,j)，k2(i,j)都是非常小的数，可以用如下的最小化问题来估计残差：

如图1所示，根据网络流理论，式(5)这个最小化问题可以转化为最小费用流来解决。最小费用流问题的输入为各个节点的度，即各个残差的值和每条流的费用，式中c1(i,j)，c2(i,j)可以根据像素点的相关度得出，而该问题的输出为各条流的流量。式(5)的最小化问题可转化为如下公式：

图1 网络图

(6)

min(0,k1(i,j))

(7)

min(0,k2(i,j))

(8)

在计算出k1(i,j)和k2(i,j)之后，再计算相位梯度，最终有

(9)

由后面的实验结果可以看出，通过网络流方法的解缠结果和枝切法、最小二乘法相比有很大的优势，没有“解缠孤岛”的现象，并且解缠精度要比其他两种方法高。但是在噪声较大的环境下网络流法的解缠结果会将噪声不可避免地沿着积分路径传递，造成解缠结果的严重“尖峰毛刺”现象，造成解缠结果的失真。很明显这些随机出现的“毛刺”点是由于噪声的干扰引起的。可以设计方法识别出这些有误差的点，然后根据这个点与周围其余点的联系重新构造出更准确的值。通过误差点识别方法找出这些由于噪声出现的“毛刺”点的集合，假设这些点(i,j)∈S1(i=0,1,…，n-1，j=0,1,…，m-1)，设集合S1的残差值k1(i,j)，k2(i,j)：

i=0,1,…,n-1,j=0,1,…,m-1

(10)

i=0,1,…,n-1,j=0,1,…,m-1

(11)

误差点的重构值是根据其与周围点相关程度和曲面趋势进行的重构，经过重构后的相位值在一定程度上降低了噪声的干扰，减少了噪声的叠加,所以可以得出

从而得出经过重构后的网络图的目标函数：

结合式(5)、式(12) 和式(13)可以得出

目标函数最小化是相位解缠的本质，以上的推导可以证明通过对网络流算法解缠结果中的噪声点进行精准识别、定位与重构后，可以有效改善“毛刺”现象，并且能够进一步最小化目标函数，从而提高解缠结果的准确度。考虑到解缠算法的时间成本以及内存要求，对于误差点的识别和利用曲面拟合重构数值分别设计了以下方法实现。

1.1 误差点处理

受到枝切法中识别“残差点”思想的启发，按照如下方法对误差点进行识别并处理。根据二维解缠原理连续相位的假设，首先设置π为阈值，再依次检验所有解缠后的像素点。在像素点周围3×3的区域内进行检测，如果像素点与其所有周围的8个点的差值的绝对值都大于阈值，那么识别此像素点为误差点。如图2(a)所示，3×3区域的像素点中心点明显比其周围的8个像素点的值大很多。因为经过网络流处理后的数据在大概率上可以代替真实相位值，那么基于连续相位的解缠假设，如果出现上述的像素点就可以认定此点由于噪声的干扰造成了数据的严重失真，可以将此点识别记录，认定为误差点，称之为E1误差点。

图2 E1,E2误差点三维示意图

式中，DI是像素点的相位值，而CI是根据像素点的相关度得出的权值。最终得到的D5就是利用最小曲面拟合构造的新像素点的值。

如图3(b)所示的12个像点，其中N5和N8为3×4区域识别的误差点，而其余的10个点均是认定的可靠数据。在这种情况下，将这个区域分割成两个3×3的区域，如图4所示。

图3 E1,E2误差点二维示意图

图4 E2误差点分割示意图

(16)

图5 E3误差点俯视图

1.2 新算法步骤

根据以上描述的方法，实现的算法步骤如下：

步骤1：根据缠绕相位信息计算出像素点间相位的相关系数，利用原始相位信息以及各个像素点相关系数等信息，结合图论中图的概念构造网络图；

步骤2：对初步构造的图，以图论的最大流最小割定理(max flow/min cut theory)为基础的增广路径算法求得最小费用路径；

步骤3：根据求得的路径对相位进行积分求得初步的解缠结果，并保存积分路径；

步骤4：对初步解缠结果进行第1轮误差点识别，识别出所有E1点，使用最小曲面拟合方法进行数据重构，并利用重构值沿着步骤3的积分路径修正积分路径的值；

步骤5：进行第2轮误差点识别，识别出所有E2点，使用最小曲面拟合方法进行数据重构，并利用重构值沿着步骤3的积分路径修正积分路径的值；

步骤6：进行第3轮误差点识别，识别出所有E3点，使用最小曲面拟合方法进行数据重构，并利用重构值沿着步骤3的积分路径修正积分路径的值；

步骤7：重复步骤6至少两次，识别出由E4点转而来的E3点，使用最小曲面拟合方法进行数据重构，并利用重构值沿着步骤3的积分路径修正积分路径的值。

2 实验结果

为了验证算法的有效性和正确性，分别用枝切法、等权最小二乘法、网络流法以及改进算法进行实验，并将结果进行对比。实验分为仿真数据实验和实测数据实验。模拟数据中加入不同程度的随机噪声和椒盐噪声，构成两组不同程度的信噪比进行仿真实验，用以验证不同信噪比情况下方法的正确性和精确度[15]。评定的参数用运行时间和均方根误差(REMS)两个参数来衡量，均方根误差计算公式如式(17)所示。实验环境为Matlab R2004a，计算机运算条件为Intel(R) Xeon(R) CPU E3-1241 v3+Windows 7专业版64 bit+8 G内存。

2.1 模拟数据仿真

如图6所示，模拟数据是Matlab软件中peaks函数产生模拟相位曲面，大小为512像素×512像素。对模拟曲面加入随机噪声和椒盐噪声实现信噪比(SNR)为6.11 dB和-0.51 dB的模拟数据。图7(a)是信噪比为6.11 dB仿真相位的曲面图，图7(b)是缠绕相位的俯视图。将数据进行缠绕后分别运用枝切法、等权最小二乘法、网络流法和新方法进行解缠仿真。解缠结果如图8所示，运行结果如表1所示。

图6 模拟数据示意图

图7 SNR=6.11 dB缠绕相位图

算法运行时间/s均方根误差/rad枝切法25.891.002最小二乘法2.82361.4995网络流法24.590.795新算法30.190.5733

图9(a)是信噪比为-0.51 dB模拟数据的缠绕图，图9(b)是图9(a)的俯视图。将数据分别运用枝切法、最小二乘法、网络流法和新算法进行解缠仿真。解缠结果如图10所示，仿真运行结果汇总如表2所示。

图8 SNR=6.11 dB模拟数据实验结果

图9 SNR=-0.51 dB缠绕相位图

算法运行时间/s均方根误差/rad枝切法42.27632.4258最小二乘法2.5721.6207网络流法28.47341.3214新算法35.03620.5861

图10 SNR=-0.51 dB模拟数据实验结果

2.2 真实数据实验

图11 实测数据

算法运行时间/sσ枝切法127.6868—最小二乘法2.57403.1157网络流法48.92193.0344新算法59.05882.6157

图12 实测数据实验结果

3 结束语

本文的方法借鉴枝切法思想，采用网络流法和曲面拟合的方法，对网络流算法增加了最小曲面拟合的步骤。以最小的时间代价，利用网络流解缠算法过程中的中间参数和初步解缠结果，对噪声引起的错误进行了可靠重构，并沿着解缠积分路径更正了积分过程中误差的传递。无论从模拟数据实验还是真实数据实验来看，新方法都大大地增加了相位解缠精度，并证明了方法的可靠性，相比其他方法具有一定的优越性。