基于最小无穷范数的二维相位解缠绕方法与流程

本发明属于雷达信号处理技术领域，尤其涉及一种基于最小无穷范数的二维相位解缠绕方法，适用于干涉合成孔径雷达(insar)系统的相位解缠绕问题。

背景技术：

作为传统合成孔径雷达(sar)的发展，insar系统能够全天候、全天时工作，这项技术在地形测绘、地表形变监测与自然灾害检测等诸多方面均有应用，因此，insar技术的发展一直受到了高度重视，而相位解缠绕(phaseunwrapping，pu)是其中的关键处理步骤。

基于l^p-norm的pu方法是最具代表性的一类方法，在传统观点下，认为p值越大越不利于相位解缠绕，因为p值越大，优化模型就越倾向于让高质量区域分担来自低质量区域(即残点密集区域)的不连续相位梯度误差，这样所获得的pu结果相对平滑，但会造成低质量区域所产生的误差向高质量区域扩散，导致求解精度的下降。

为了保证高质量区域的求解精度，传统pu观点倾向于保证求解结果与输入条纹一致，然而，尽管这种做法保证了高质量区域的求解精度，但是当输入条纹中存在低质量区域时，这些区域几乎不存在有效信息，此时要求pu结果与输入条纹一致是不合理的。

技术实现要素：

针对上述问题，本发明的目的在于提供一种基于最小无穷范数的二维相位解缠绕方法(本发明实施例中称为min方法)，在进行相位解缠绕的过程中具有自适应滤波功能，能够有效改进传统二维相位解缠绕方法的求解精度。

在实际应用中，需要一种方法，不仅能够保证高质量区域的求解精度(即在高质量区域保证条纹一致性)，还应该能对低质量区域进行合理处理(即放宽条纹一致性的限制)。

为达到上述目的，本发明采用如下技术方案予以实现。

一种基于最小无穷范数的二维相位解缠绕方法，所述方法包括：

步骤1，获取干涉相位图，所述干涉相位图的大小为m×n，其中，m、n分别为大于零的正整数；图中位置(i,j)处的像素值记为干涉相位图中位置(i,j)处的缠绕相位从而得到干涉相位图的缠绕相位矩阵；其中，i＝1,2,...,m,j＝1,2,...,n；

步骤2，根据所述干涉相位图的缠绕相位矩阵，分别计算得到干涉相位图的水平缠绕相位梯度矩阵和干涉相位图的垂直缠绕相位梯度矩阵；

步骤3，根据干涉相位图的水平缠绕相位梯度矩阵，计算干涉相位图的水平绝对相位梯度估计值矩阵；

根据干涉相位图的垂直缠绕相位梯度矩阵，计算干涉相位图的垂直绝对相位梯度估计值矩阵；

步骤4，根据所述干涉相位图的水平绝对相位梯度估计值矩阵、所述干涉相位图的垂直绝对相位梯度估计值矩阵，计算得到干涉相位图的绝对相位矩阵，所述干涉相位图的绝对相位矩阵为干涉相位图解缠绕的结果。

本发明技术方案的特点和进一步的改进为：

(1)步骤2具体包括：

(2a)干涉相位图的水平缠绕相位梯度矩阵中位置(i,j)处的水平缠绕相位梯度记为其中，i＝1,2,...,m,j＝1,2,...,n-1；

且：

其中，表示干涉相位图中位置(i,j+1)处的缠绕相位，表示干涉相位图中位置(i,j)处的缠绕相位；

(2b)干涉相位图的垂直缠绕相位梯度矩阵中位置(i,j)处的垂直缠绕相位梯度记为其中，i＝1,2,...,m-1,j＝1,2,...,n；

且：

其中，表示干涉相位图中位置(i+1,j)处的缠绕相位。

(2)步骤3具体包括：

(3a)干涉相位图的水平绝对相位梯度估计值矩阵中位置(i,j)处的水平绝对相位梯度估计值记为其中，i＝1,2,...,m,j＝1,2,...,n-1；

且：

其中，表示干涉相位图的水平缠绕相位梯度矩阵中位置(i,j)处的水平缠绕相位梯度；

(3b)干涉相位图的垂直绝对相位梯度估计值矩阵中位置(i,j)处的垂直绝对相位梯度估计值记为其中，i＝1,2,...,m-1,j＝1,2,...,n；

且：

其中，表示干涉相位图的垂直缠绕相位梯度矩阵中位置(i,j)处的垂直缠绕相位梯度。

(3)步骤4具体包括：

将干涉相位图的水平绝对相位梯度估计值矩阵中位置(i,j)处的水平绝对相位梯度估计值垂直绝对相位梯度估计值以及位置(i+1,j)处的水平绝对相位梯度估计值位置(i,j+1)处的垂直绝对相位梯度估计值带入如下优化模型：

其中，t(i,j)为自由变量，w(i,j)为t(i,j)的权值，集合s＝{(i,j)|i＝1,2...,m-1,j＝1,2,...,n-1)}，ψ(i,j)、ψ(i,j+1)、ψ(i+1,j)和ψ(i+1,j+1)分别表示干涉相位图的绝对相位矩阵在位置(i,j)、位置(i,j+1)、位置(i+1,j)和位置(i+1,j+1)处的绝对相位，是该优化模型要求解的变量。

本发明为一种具有自适应滤波功能的相位解缠绕方法，在进行相位解缠绕求解时，能够自动识别输入的干涉相位图中的高质量区域和低质量区域，在保证高质量区域的条纹一致性的前提下，还能对低质量区域的噪声进行滤波，提高了相位解缠绕的精度。

附图说明

为了更清楚地说明本发明实施例或现有技术中的技术方案，下面将对实施例或现有技术描述中所需要使用的附图作简单地介绍，显而易见地，下面描述中的附图仅仅是本发明的一些实施例，对于本领域普通技术人员来讲，在不付出创造性劳动的前提下，还可以根据这些附图获得其他的附图。

图1为本发明实施例提供的一种基于最小无穷范数的二维相位解缠绕方法的流程示意图；

图2为针对仿真干涉相位图，分别展示了本发明方法与传统方法的处理结果，所有图的横坐标对应相应矩阵的列，所有图的纵坐标对应相应矩阵的行。图2中，图2(a)为参考解缠绕相位图；图2(b)为参考干涉相位图；图2(c)为仿真的干涉相位图；图2(d)为相关系数图；图2(e)为l^∞-norm方法的pu结果图；图2(f)为l²-norm(2范数)方法的pu结果图；图2(g)为l¹-norm方法的pu结果图；图2(h)为snaphu方法的pu结果图；图2(i)为puma方法的pu结果图；图2(j)为min方法的pu结果图；图2(k)为图2(a)与2(e)之间的差值结果图；图2(l)为图2(a)与2(f)之间的差值结构图；图2(m)为图2(a)与2(g)之间的差值结果图；图2(n)为图2(a)与2(h)之间的差值结果图；图2(o)为图2(a)与2(i)之间的差值结果图；图2(p)为图2(a)与2(j)之间的差值结果图；图2(q)为l^∞-norm方法pu结果的再缠绕相位图；图2(r)为l²-norm方法pu结果的再缠绕相位图；图2(s)为min方法pu结果的再缠绕相位图；

图3为针对实测干涉相位图，分别展示了本发明方法与传统方法的处理结果，所有图的横坐标对应相应矩阵的列，所有图的纵坐标对应相应矩阵的行。图3中，图3(a)为实测的干涉相位图；图3(b)为l^∞-norm方法的pu结果图；图3(c)为l²-norm方法的pu结果图；图3(d)为l¹-norm方法的pu结果图；图3(e)为snaphu方法的pu结果图；图3(f)为puma方法的pu结果图；图3(g)为min方法的pu结果图；图3(h)为l^∞-norm方法pu结果的再缠绕相位图；图3(i)为l²-norm方法pu结果的再缠绕相位图；图3(j)为min方法pu结果的再缠绕相位图；图3(k)为图3(a)的滤波结果图；图3(l)为图3(j)的滤波结果图。

具体实施方式

下面将结合本发明实施例中的附图，对本发明实施例中的技术方案进行清楚、完整地描述，显然，所描述的实施例仅仅是本发明一部分实施例，而不是全部的实施例。基于本发明中的实施例，本领域普通技术人员在没有做出创造性劳动前提下所获得的所有其他实施例，都属于本发明保护的范围。

本发明实施例提供一种基于最小无穷范数的二维相位解缠绕方法，如图1所示，所述方法包括：

步骤1，获取干涉相位图，所述干涉相位图的大小为m×n，干涉相位图中位置(i,j)处的像素值记为干涉相位图中位置(i,j)处的缠绕相位从而得到干涉相位图的缠绕相位矩阵；其中，i＝1,2,...,m,j＝1,2,...,n。

步骤2，根据所述干涉相位图的缠绕相位矩阵，分别计算得到干涉相位图的水平缠绕相位梯度矩阵和干涉相位图的垂直缠绕相位梯度矩阵。

步骤2具体包括：

(2a)干涉相位图的水平缠绕相位梯度矩阵中位置(i,j)处的水平缠绕相位梯度记为其中，i＝1,2,...,m,j＝1,2,...,n-1；

且：

其中，表示干涉相位图中位置(i,j+1)处的缠绕相位，表示干涉相位图中位置(i,j)处的缠绕相位；

(2b)干涉相位图的垂直缠绕相位梯度矩阵中位置(i,j)处的垂直缠绕相位梯度记为其中，i＝1,2,...,m-1,j＝1,2,...,n；

且：

其中，表示干涉相位图中位置(i+1,j)处的缠绕相位。

步骤3，根据干涉相位图的水平缠绕相位梯度矩阵，计算干涉相位图的水平绝对相位梯度估计值矩阵；

根据干涉相位图的垂直缠绕相位梯度矩阵，计算干涉相位图的垂直绝对相位梯度估计值矩阵。

步骤3具体包括：

(3a)干涉相位图的水平绝对相位梯度估计值矩阵中位置(i,j)处的水平绝对相位梯度估计值记为其中，i＝1,2,...,m,j＝1,2,...,n-1；

且：

其中，表示干涉相位图的水平缠绕相位梯度矩阵中位置(i,j)处的水平缠绕相位梯度；

(3b)干涉相位图的垂直绝对相位梯度估计值矩阵中位置(i,j)处的垂直绝对相位梯度估计值记为其中，i＝1,2,...,m-1,j＝1,2,...,n；

且：

其中，表示干涉相位图的垂直缠绕相位梯度矩阵中位置(i,j)处的垂直缠绕相位梯度。

步骤4，根据所述干涉相位图的水平绝对相位梯度估计值矩阵、所述干涉相位图的垂直绝对相位梯度估计值矩阵，计算得到干涉相位图的绝对相位矩阵，所述干涉相位图的绝对相位矩阵为干涉相位图解缠绕的结果。

步骤4具体包括：

将干涉相位图的水平绝对相位梯度估计值矩阵中位置(i,j)处的水平绝对相位梯度估计值垂直绝对相位梯度估计值以及位置(i+1,j)处的水平绝对相位梯度估计值位置(i,j+1)处的垂直绝对相位梯度估计值带入如下优化模型：

其中，t(i,j)为自由变量，w(i,j)为t(i,j)的权值，集合s＝{(i,j)|i＝1,2...,m-1,j＝1,2,...,n-1)}，ψ(i,j)、ψ(i,j+1)、ψ(i+1,j)和ψ(i+1,j+1)分别表示干涉相位图的绝对相位矩阵在位置(i,j)、位置(i,j+1)、位置(i+1,j)和位置(i+1,j+1)处的绝对相位，是该优化模型要求解的变量。

通常，w(i,j)可以由干涉相位图的质量图获得，当不需要加权的时候，令w(i,j)＝1即可。上述优化模型是一个线性规划模型，求解该模型即可获得解缠绕结果。

需要说明的是，该优化模型是在l¹-norm(1范数)的基础上，对每个相邻的四个像素上的绝对相位(即ψ(i,j)、ψ(i,j+1)、ψ(i+1,j)和ψ(i+1,j+1))采用l^∞-norm(无穷范数)的进行约束，因而，l^∞-norm所具有的滤波功能被引入到l¹-norm优化模型中，这样，该优化模型既能够通过l¹-norm在干涉相位图的高质量区域保持条纹一致性(即不进行滤波)，同时也能利用l^∞-norm在干涉相位图的低质量区域起到的滤波的作用，由于不需要任何辅助手段来帮助l^∞-norm识别低质量区域以完成滤波功能，可以说，本发明所提出的相位解缠绕方法具有自适应滤波功能。

本发明的有效性可通过以下实验作进一步说明。

(一)insar仿真数据的实验

1、数据说明

仿真数据采用中国河北省尚义地山区的数字高程模型，在仿真生成干涉相位图时，为了增加相位解缠绕的难度，加入了噪声，在高质量区域加入了相关系数为0.8的噪声，在低质量区域加入了相关系数为0的噪声。

2、仿真内容和结果分析

为了展示本发明方法的有效性，仿真实验将本发明方法与l^∞-norm、l²-norm、l¹-norm、snaphu、puma等传统方法进行了对比实验。

图2(a)所示为由数字高程模型dem得到的参考绝对相位。图2(b)为由图2(a)得到的理想的干涉相位图。图2(c)所示为仿真的干涉相位图，图2(d)为在仿真生成图2(c)时所采用相关系数，为了增加pu的难度，在图2(c)中没有进行滤波处理。此外，从图2(d)中能够看到，在仿真的时候给图2(c)加入了几个噪声区(高质量区域仿真用的相关系数为0.8，噪声区域仿真使用的相关系数为0)。图2(e)、图2(f)和图2(g)分别为l^∞-norm方法、l²-norm方法和l¹-norm方法的pu结果。图2(h)为snaphu方法的pu结果，其统计代价模式设置为topo。图2(i)为puma方法的pu结果，puma方法的cliquepotentialexponent参数设置为0.5。图2(j)是本发明方法的pu结果。所有的方法均在没有加权的情况下进行pu求解。图2(k)-图2(p)分别是图2(e)-图2(j)与图2(a)的差值，对图2(k)-图2(p)的统计信息罗列在表1中。

表1.pu结果的统计信息

从图2中可以看到，l^∞-norm和l²-normpu方法能够得到非常平滑的pu结果，但是它们的pu精度十分低。l¹-norm，puma，snaphu和min方法都能够在高质量区域保证pu精度，所以它们总体上的性能优于l^∞-norm和l²-norm方法。然而，由于min方法对低质量区域具有自适应滤波功能，图2(j)所示的min方法pu结果中的噪声区域相对于图2(g)-(i)更为平滑，也就是说，在低质量区域，min方法的pu精度高于l¹-norm，puma和snaphu方法。通过前述对比，可以得出结论，min方法的滤波功能能够降低pu误差。图2(q)、图2(r)和图2(s)展示了图2(e)、图2(f)和图2(j)的再缠绕相位图，因为l¹-norm，puma和snaphu方法能够严格保证条纹一致性，所以图2(g)、图2(h)和图2(i)的再缠绕条纹图与图2(c)一样。通过将图2(c)与图2(q)、图2(r)、图2(s)进行对比，可以看出l^∞-norm和l²-normpu方法对输入干涉相位图产生了破坏性滤波，但是min方法能够在高质量区域保证条纹一致性，并且在对噪声区域进行了自适应滤波。此外，通过对比图2(b)和图2(s)，可以看出，、由min方法恢复出来的噪声区的相位条纹，是比较接近图2(b)所示的理想条纹的。

(二)实测数据实验

图3(a)是实测干涉相位图。可以在其中清晰的看到“好”区域与“坏”区域。图3(b)、图3(c)和图3(d)分别是为l^∞-norm方法、l²-norm方法和l¹-norm方法的pu结果。图3(e)为snaphu方法的pu结果，其统计代价模式设置为topo。图3(f)为puma方法的pu结果，puma方法的cliquepotentialexponent参数设置为0.5。图3(g)是本发明方法的pu结果。所有的方法均在没有加权的情况下进行pu求解。图3(h)、图3(i)、图3(j)分别是图3(b)、图3(c)、图3(g)的再缠绕相位图，由于图3(d)、图3(e)图3(f)的再缠绕条纹与图3(a)相同，故在此不再展示。能够看到图3(b)、图3(c)、图3(g)的噪声区的pu结果比图3(d)、图3(e)、图3(f)的噪声区光滑。然而，通过将图3(a)与图3(h)、图3(i)相比，很明显，l²-norm方法、l^∞-norm方法对输入干涉相位图进行了过分的滤波，即在高质量区域的输入相位条纹被完全破坏了。在这种情况下，高质量区域的pu求解结果的精度被大大降低了。l¹-norm方法、snaphu方法、puma方法严格保证了条纹一致性，因此高质量区域的pu精度是可信的。如果对比图3(a)和图3(j)，可以看出min方法的滤波功能仅仅对于噪声区域发生了作用。因此，min方法不仅保证了好区域几乎与l¹-norm方法、snaphu方法和puma方法一样，而且对坏区域的输入相位条纹进行了滤波。为了进一步展示min方法的滤波效果，采用了简单的滤波方法(即3×3窗口的均值滤波)对图3(a)与图3(j)的条纹进行再次滤波，滤波的结果如图3(k)与图3(l)所示。因为对图3(a)和图3(g)再次滤波并非insar处理中的合理步骤，所以图3(k)和图3(l)的相位条纹过于平滑。可以看到，图3(a)的低质量区域在图3(k)中仍然是噪声。然而，在图3(l)中，对图3(a)噪声区域，产生了一些合理的新条纹。换言之，图3(l)揭示出min方法的滤波效果的趋势。

综上，我们分别采用了仿真数据与实测数据验证了本发明的有效性。

以上所述，仅为本发明的具体实施方式，但本发明的保护范围并不局限于此，任何熟悉本技术领域的技术人员在本发明揭露的技术范围内，可轻易想到变化或替换，都应涵盖在本发明的保护范围之内。因此，本发明的保护范围应以所述权利要求的保护范围为准。