InSAR相位解缠最小范数算法的研究

1InSAR相位解缠最小范数算法的研究第一章绪论1.1论文研究的背景1.1.1合成孔径雷达干涉测量概述合成孔径雷达干涉测量(InSAR)是20世纪60年代末发展起来的一项技术，在近20年来受到了世界各国的广泛关注获得了迅猛发展并逐渐趋于成熟。由于合成孔径雷达干涉测量主要是利用主动微波式传感器，它的出现大大地扩展了合成孔径雷达、光学传感器等的应用领域。它不仅能够获取高精度的高程信息，同时还可以全天时、全天候监测陆地表面和冰雪表面地形等的微小变化，监测的时间间隔从几天到几年，监测精度可达毫米级，并且它对某些目标物体还具有一定的穿透能力。其更令人瞩目的是，这项技术还可用于研究过去长时间无法到达的冰川和冰源的变化情况，也可用于一些灾害性地表形变的探测，如地震、火山爆发、等以及地表三维的重建，因而成为了遥感研究的热点川。1.1.2相位解缠研究的现状相位解缠技术最早出现在20世纪60年代末70年代初，当时主要是信号处理的需要，所研究的主要是一维问题。除合成孔径雷达干涉测量中应用外，还在合成孔径声纳、光学干涉、微波干涉、核磁共2振等方面有重要应用。二维相位解缠始于20世纪70年代末。在过去的30多年里，InSAR的相位解缠的方法发展十分迅速，达到了三、四十种，文献(王超，2002)列出了多种算法，但以上基本上可以分为两大类，即路径跟踪法(PathFollowing)和最小二乘法(LeastSquare)，路径跟踪法基于像元到像元的局部运算来解缠，而最小二乘法是通过使解缠后解缠前相位的梯度差整体最小来进行求解的。各种算法都有其自身的优缺点，适用于特定条件的数据，普适性都不是很好，因此算法的选择一般应根据实际情况而定。1.2本文研究内容我国是一个地质灾害频繁的国家，近些年来各种地质灾害接踵而来，如地震、滑坡、地面沉降等，这些地质灾害以地表形变为直接特征，严重影响了人民生命与则一产的安全，因此对地表形变的监测显得尤为重要。合成孔径雷达技术能够利用雷达信号中的相位信息来提取地表的三维信息，精度可达毫米级，己成为目前DEM生产的主要技术手段之一，在地下资源探测以及军事目标探测等方面都具有其独特的优越性和发展潜力。相位解缠作为InSAR技术应用处理中至关重要的一个环节，也因此显得尤为重要。本文主要研究内容包括以下几个方面:1、对相位解缠中最小范数算法的理论进行归纳和研究.2、从对合成孔径雷达干涉测量的常用数据分析入手，在C#编程语言的基础上，结合WPS,GIS等技术和手段，对基于最小范数算法的InSAR相位解缠软件的四种基于最小范数相位解缠算法，包括3DCT、无权重多网格算法、有权重多网格算法及PCG算法进行实现。3、对软件实现的基于最小范数相位解缠算法进行精度分析。4、从数据实例中得到了最小范数相位解缠算法的某些特征，方便用户选择和使用最小范数相位解缠算法。第二章最小范数的相位解缠算法2.1引言一切将相位从主值或相位差值恢复为真实值的过程称为相位解缠。从数学上讲，给定一个二维相位矩阵j,i，为了解缠该矩阵，需要对每个点(i,j)加上或减去2二的整数倍而得到一个连续的函数ji,，即有（2.1-1）其中jiK,为整数。进行相位解缠必须兼顾两个方面:一致性和精确性。一致性是指解缠后的矩阵ji,中任意两点间的相位差与这两点之间的路径无关，精确性则是指解缠后的相位要能忠实地恢复原始相位函数。目前所有相位解缠可分为两个步骤:第一是基于缠绕相位计算解缠相位的相位梯度估计值;第二是进行积分。根据所采用的积分方法，相位解缠方法分为两大类，即路径跟踪法和最小二乘法。现有的相位解缠算法都基于这样一个假设:有可能推导出缠绕相位的离散偏导数，即4邻近像元的相位差，并且这些相位差的绝对值都小于二。通过这些离散的偏导数，可以重建解缠相位。在理想状态下，干涉相位呈周期变化(见图1)，由0渐变到2π，再由2π迅速下降到0，然后又渐变到2π，周而复始，呈现周期性;其变化轮廓分明，层次均匀，突变点为相位周期分界点。因此，在理想的状态下如能提取出垂直向和水平向的相位偏导数，然后沿两个方向积分，就达到相位解缠的目的。实际上，星载或机载InSAR数据普遍存在着由于地形起伏较大而引起的密集的干涉条纹。由于地形起伏引起的顶底位移、雷达阴影、雷达信号低信噪比和各种原因造成的去相干现象，还由于在原始雷达信号处理过程中引入的相干噪声、噪声和伪信号等造成的相位数据的不连续，解缠相位的离散相位梯度估算值将不能保持一致，也就是说，它们不能形成一个“无旋”的向量场。此时虽然干涉相位呈一定的趋势和周期性，但并不明显。特别是当干涉相位从2π二变化到0时界限模糊(见图2)，因而为各周期的分离带来极大的困难，所以很难用简单的积分方法达到恢复真实相位的效果。52.2数学表达对于给定的缠绕相位函数ji,(i=0,1,2..,.，M-1;j=0,1,2,...，N-1)，设其对应的解缠相位函数为ji,，要使解缠相位函数ji,的相位差在最小pL范数意义下满足缠绕相位函数，则ji,，需满足6（2.2-1）其中x,ji、yji,，表示像元(i,j)的缠绕相位差，设全变差用J来表示，则上式有:（2.2-2）其中2,1，的定义如下:（2.2-3）（2.2-4）令（2.2-5）（2.2-6）将（2.2-5），（2.2-6）代入（2.2-3）（2.2-4）中有：7（2.2-7）（2.2-8）令边界满足ji,a=ji,b=0且J0且s}=o，并将i.j互换，得到，（2.2-9）（2.2-10）将上式代入（2.2-2）有（2.2-11）假设扰动，是随机的，则应有（2.2-12）其中0≤-i≤-M-1:0≤j≤-N-1。要求解上式，需要知道j1-，，jm,，1,i，N,i及对应的j1-，jm,，1,iNi,'值，根据边界条件ji,a=0，i=-1M-1;0-j-N-1:以8及jib,=0,r=-1,N一1;0_i_M一1有(2.2-13)(2.2-14)展开式（2.2-12），并对上式加上边界条件，将此问题转换为最小范数问题，如下:（2.2-15）其中（2.2-16）（2.2-17）当p=2时，式(2.2-15)转换为了离散形式的泊松方程：（2.2-18）9其中（2.2-19）其边界条件为（2.2-20）（2.2-21）即当p=2时，基于最小pL范数的相位解缠问题就转化为了求解牛曼边界的泊松方程。2.3无加权最小二乘相位解缠算法最小二乘相位解缠问题，无论是有无权重，都等效于要用有效的数学方法去解决离散的偏微分方程。最小二乘法为线性方程，其算法结果又是有效的，可跟踪的，因此使用非常广泛最小二乘相位解缠方法的基本思想:使缠绕相位的离散偏微分导数和解缠相位的离散偏微分导数差的平方和最小，即，10（2.3-1）将式(2.2-18)写成KxK线性系统下的矩阵形式，有（2.3-2）其中，P是式2.2-18)等式左边的离散拉普拉斯算子的稀疏矩阵。为一维矩阵向量，K=MN是像元数。一维矢量中的和二维之间的关系为:k=i+jM,i=0,1,M-1,j二。,1,...N-1，同样的Pk和其对应的jiP,，也有此关系。2.4加权最小二乘相位解缠算法无加权的最小二乘算法的缺点是穿过奇点区域而不是绕过它们进行解缠，这可以通过引入相位数据的权重而得到克服。对于那些因噪声或其他退化作用而被打断的相位，可将其权重赋予0,使其不影响解缠过程。如果为给定数据赋予权重(0≤jW,i≤1)，那么最小二乘问题就转化为加权最小二乘问题，故式（2.4-1）11可改写为（2.4-2）其中U（i,j）和州V(i.j)为梯度权重:(2.4-3)于是，加权最小二乘解ji,，，由如下方程求解:(2.4-4)其中jiC,为相位拉普拉斯算子:(2.4-5)求解式(2.4-4)的高斯赛德尔松弛算法也可以采用加权的迭代形式:(2.4-6)12与无权重的形式相同，高斯赛德尔松弛算法的收敛速度较慢，因此可采用加权多重网格技术，加快算法的执行速度。2.5无权重的多网格算法多网格法是一类经典的在较大网络基础上解决偏微分方程的快速方法。这类方法是基于在更粗、更小的网格上应用高斯一赛德尔松弛算法的思想。特别地，对于解决M×N网格上离散PDE问题，多网格法与直接傅里叶变换速度一样快。但是，与直接算法不同的是，多网格法可以解决更多的问题，包括非线性的PDE问题。并且从实用角度上来看，多网格还具有一个FFT法和DCT法不具备的优点，即它对网格的大小没有限制，不要求数组的大小为2的幂次。高斯一赛德尔松弛法是一个平滑算子，能很快的去除误差的高频部分，但去除低频部分则非常的慢。因此，多网格法的关键在于将误差的低频成分转换为高频成分，再将这些高频成分用高斯一赛德尔松弛法去除。在实际处理中，则是通过把这个问题转移到更粗的网格中去实现。更粗网格的较低样本率增加了剩余误差的空间频率。如果有一个金字塔形的网格(见图1)，下层网格的分辨率只有它上面一层的一半，分辨率降低了，而空间频率在密网格和粗网格中不断传递，将在粗网格中产生不断增加的较高的空间频率，换句话说，在密网格的全局信息、将变成粗糙网格中的局部信息。把计算结果转移到粗网格的操作数称为约束操作数或密到粗操作数，而把计算结果转移到更密网格的操作数称为延伸操作数或粗到密操作数。13图1多网格模式下的分辨率变化金字塔约束操作数的一般选择为注入操作数，它每隔一个像元进行采样，定义为jiC,=jif2,2，其中jif,是MxN的密网格，jiC,是M/2xN/2的粗网格。但是，在实际应用中常采用一种效果更好的约束操作数，称为“全重”，其定义为:14（2.5-1）实际上，全重操作数是用下面的滤波模板先平滑密网格（2.5-2）然后每隔一个点进行采样。延伸操作数为一个双线性插值，其定义为:（2.5-3）约束操作数把方程的“残差误差”转换到更粗的网格对其作松弛处理，等同于对方程木身作松弛处理。如果把相位解缠问题表示成线性问题，，其中P是矩阵，和P是向量，那么“残差方15程’，可以表示如下:（2.5-4）其中是的近似，是未知误差。用初始的估计值对方程进行松弛运算，等同于用初始估计值=0对剩余方程进行松弛运算，其中:是未知的剩余误差。实际上，近似值必是在一个网格上通过高斯一赛德尔松弛法得到的中间值，初始的估计值=0是在下一个更粗网格的高斯一赛德尔松弛法的起点。在粗网格得到的结果可依次视为一个中间值，它的“剩余”误差被限制在更粗的网格。另外，这个解法可以被延伸(或插值)到更密的网格，并且加上近似值功，在更密网格中得到更好的结果。2.6加权多网格法与无权重算法类似，加权多网格法的核心思想是:将缠绕相位的离散偏导数xji,和y,ji存储在不同的数组中，并且在每一个从密网格到粗网格的循环中对在较粗网格边界的导数进行纠正，这个过程保证了边界条件在更粗网格上收敛到正确结果。对于加权多网格法的挑战是:对更粗网格上加上正确的权重限制。另外，加权算法还需要另外两个数组用于存储权重，一个数组存储用于行导数的权重值，一个数组存储用于列导数的权重值。因此，此算法共需要5个数组，每个数组的大小与输入相位数据的大小相同:(解缠相位)、xji,(缠绕相位行导16数)、y,ji(缠绕相位列导数)、(行导数的权重)和(列导数的权重)。这个算法所需要的内存空间是5*(1+1/4十1/16+...)≈7倍的输入相位数据所需的内存空间。我们把权重数组和:分别称为“权重”和“权重”，这些权重的初始值被定义为输入的权重值。这些权重由于粗糙网格约束操作数的定义而不同，约束操作数表示为Ro。式中的权重系数定义如下:（2.6-1）其中和:分别称为“权重”和“权重”。约束操作数风，(约束粗糙网格导数的残差)在应用全重操作数生成一个新的粗网格值之前必须给密网格加权重。为了约束d/dx残差导数，Rw在用加权的滤波模板对密网格值进行平滑后每隔一个点进行重采样，此加权的滤波模板定义为（2.6-2）其中，jiS,是所有输入矩阵和。如果所有权重都是0，那么就采用上式的