基于迭代高斯模型的干涉DEM-滤波算法

第卷第期遥感技术与应用Vol.No.2012年月REMETESENSINGTECHNOLOGYANDAPPLICATION,2012收稿日期：2012-3-23修订日期：基金项目:“国家863项目”（编号：2007AA120302）支持作者简介：王峰（1988-），男，山东日照，博士研究生。主要从事遥感图像处理研究。Email：[e-mail]。尤红建（1989-），男，江苏如皋，研究员，博士生导师。主要从事遥感图像处理和应用领域的研究。Email:[e-mail]。傅兴玉（1984-），男，山东菏泽，助理研究员。主要从事地理空间信息处理与应用。Email:[e-mail].通讯作者：王峰（1988-），男，山东日照，博士研究生。主要从事遥感图像处理研究。E-mail：[e-mail]。基于迭代高斯模型的干涉DEM滤波算法王峰1,2,3，尤红建1,2，傅兴玉1,2（1.中国科学院空间信息处理与应用系统技术重点实验室，北京1001901；2.中国科学院电子学研究所，北京100190；3.中国科学院研究生院，北京100049）摘要：提出了一种基于迭代高斯模型的干涉DEM滤波算法。该算法首先对含有噪声的干涉SAR高程数据进行噪声检测，对不符合高斯模型分布的数据，标记为噪声点，建立噪声标记矩阵；然后对标记为噪声的点，依据其邻域窗口内数据的统计特性，采用曲面拟合算法估计得出噪声点处的真实高程值。通过实际DEM的实验结果表明，文章中提出的算法在有效消除噪声点的同时，可以较好地保持地面的高程值；与传统的低通滤波、中值滤波、Sigma滤波的实验结果相比，该算法的滤波结果较为理想。关键词：干涉DEM滤波；噪声标记矩阵；迭代高斯模型；曲面拟合中图分类号：P237.3文献标识码：A1引言DEM（DigitalElevationModel，数字高程模型）是地形起伏的数字表达，是地学分析中的基础数据。利用遥感途径获取DEM的手段主要包括：光学遥感立体像对、激光雷达和干涉SAR（SyntheticApertureRadar）。其中，干涉SAR利用两部天线接收到的回波信号的相位差获得地面点的高程值[1]。干涉SAR具有其他DEM获取方式所不具备的优点：它可以全天时全天候工作，同时测绘带较宽，可以进行大面积快速测绘。干涉SAR主要利用相干成像的特点进行干涉测量，从而获得高程数据，但由于SAR数据的一个分辨单元内包含有许多的微散射体，SAR图像上每个像素的信号数据实际上是这些电磁波与微散射体相互之间加强或减弱作用的集成，在影像中以斑点的形式表现出来，会产生相干波成像固有的斑点噪声（Speckle）问题[2]。在干涉SAR处理中，相位信息的保持是获得高精度地形拓扑图的关键因素之一。但是由于系统热噪声和斑点噪声以及相位突变情况的存在，造成某些区域的相位计算存在较大误差，计算得到的地面高程值精度较差，影响DEM数据的质量。斑点噪声是SAR图像所固有的，其成因要从SAR的成像原理和回波信号的统计特征来分析[3][4]。热噪声引起干涉去相关，进而造成高程奇异值斑块，影响干涉SAR获得的DEM的应用[5]。干涉SAR利用同一地区获得的不同SAR复图像的相位信息进行干涉合成，然后通过多视处理、去平地效应、相位解缠，获取地面高程信息[6]。但是对于河流、湖泊等相干性较差的区域，干涉相位求解过程会受到影响，无法解缠得到正确的相位值，导致得到的DEM数据中存在高程畸变，在数据中表现为突兀的斑块。由于解缠过程中的相位误差传播等因素，造成此类斑块的面积较大，不能使用一般的滤波算法进行消除。综上，InSAR（InterferometricSyntheticApertureRadar）处理得到的DEM数据中，一般存在随机分布的斑点噪声以及面积较大的高程畸变区域，这些随机分散的高程奇异值斑块，明显高于或低于周围像素点的高程值，严重影响DEM数据的质量和应用。从SAR的成像原理和回波信号的统计特性来分析[3][4]，同时考虑到地表地形一般的特性，DEM数据一般具有空间连续性和空间自相关性[7]，即地表高程值空间连续变化，并且两个点的平面距离越小，高程值越接近。同时考虑到干涉SAR获取的DEM中奇异值斑块空间分布具有一定的随机性，我们提出首先建立窗口限制空间范围，利用迭代高斯模型遥感技术与应用检测噪点，然后用二次曲面拟合插补干涉SAR获取DEM数据中高程奇异值斑块的算法，可以有效抑制噪声数据的影响。2基于高斯模型滤波算法噪声严重影响着干涉SAR的图像质量，使得目标高度精度降低。为获得高质量的干涉SAR图像，必须对噪声进行有效的抑制，同时尽量保持干涉SAR图像的分辨率。针对干涉SAR图像中出现的随机噪声，抑制的方法大体可以分为两类：频域的多视处理与图像域滤波。图像域的滤波，如低通滤波器，中值滤波器等是以牺牲空间分辨率为代价的。这会影响实际DEM数据的空间相关性和空间连续性，不能直接应用于DEM数据的滤波。为了尽量保持原始数据的正确性，我们在检测噪声的基础上再进行曲面拟合。因此算法步骤主要有两个部分：（1）通过高斯模型迭代检测出噪声点。（2）针对每个噪声点，建立一个以此点为中心的窗口，根据真实地势变化情况，根据未标记为噪声点的真实高程数据，拟合二次或三次曲面，计算出噪声点处的高程值。2.1高斯模型检测噪声点一般情况下，实际地形起伏具有连续性和内部相关性的特点，高程变化比较平缓，很少会出现跳变剧烈的点，从InSAR数据的统计特性分析，结合实际地势起伏的空间相干性，DEM数据在一定窗口内的变化应该是平缓的，高程值大部分分布在均值附近，较大的波动点是很少存在的，换句话说，应该近似符合高斯分布。如图1(a)为一块300×300像素，无斑点噪声的DEM数据，图(b)为对应的高程数据直方图分布，近似满足高斯分布。这就表明以高斯分布估计窗口区域的DEM数值统计特性是合理的。(a)无噪声点DEM数据(b)直方图分布图1无噪声点的DEM数据及其直方图分布Fig.1DataandhistogramofDEMwithoutnoises对于获取DEM数据的InSAR来说，在地势变化平缓，相干性好的区域，产生散斑噪声的可能性较小；而在起伏比较剧烈的山区，相干性较差的地区，产生散斑噪声的可能性明显加大。因而，结合图像不同区域的相干特点，利用相干性测度，如公式（1）所示，合理设置图像散斑噪声检测中高斯模型的系数以及窗口大小，可以充分利用区域特征，尽量保持DEM高程数据的真实性和连贯性。*122212EEEµµγµµ=(1)为了最大限度的使用像素领域内的数值统计特性，应以待检测点为中心建立窗口，分析整个窗口内的高斯分布情况。若待检测点符合窗口内的高斯分布情况，则认为是正确的DEM数据，如果不是，则标记为噪声点。2(,)((,))fxyEfxygασ−=(2)如公式(2)所示，其中(,)xy是窗口中心点坐标，(,)fxy是窗口中心点的灰度值，((,))Efxy是窗口范围内，除中心点以外的灰度的均值，σ是窗口范围内灰度方差，a是高斯分布检测阈值。当窗口中心点数据符合整个窗口范围内的高斯分布时，满足公式(2)。高斯分布检测阈值a，会影响检测的正确率。当a值设置偏大时，检测的虚警较低，但是会有一些噪声点没有检测出来；当a值设置偏小时，会出现虚警，把正确数据误检测成错误数据。而对于不同数据，噪声对数据的影响程度不同，具体阈值的选择应该根据实际数据实验得到。图2高斯函数在不同值时的分布曲线Fig2.Distributingcurveofgaussatdifferenta如图2是一维高斯函数（公式3）在不同a值时的分布曲线，在实际检测中，a的取值范围可以选取1.95~2.58之间。22()212xfxeµσπσ−−=(3)遥感技术与应用对原始数据进行一次噪声点检测，可以检测出部分噪声点。对于数据中，分布稀疏的孤立噪声点可以很好的被填补上，同时算法只对检测出来的噪声点进行插补，保留原来正确数据，原始数据中的正确细节被完整的保存下来。但是对于噪声点分布密集的区域，由于检测窗口内存在较多错误数据的影响，会出现噪声点没有被检测出来的问题。分析原因，对于噪声点密集的区域，当以密集噪声点为中心建立窗口时，窗口内包含的大部分是错误点，而正确数据所占的比例相对较少，从而在建立高斯模型的时候，错误数据所占的权重超过了正确数据的权重。为解决这个问题，需要用迭代算法检测密集区域的噪声点。2.2迭代检测噪声点如上所述，进行一次高斯模型检测，建立噪声标记矩阵iM，假设噪声点的个数为iN。在此基础上，结合噪声标记矩阵iM，只利用窗口内认为正确的DEM数据重新建立高斯模型，对上一步骤中未标记为噪声点的点，再进行一次高斯模型检测，得到新的噪声标记矩阵1iM+，假设现在噪声点的个数为1iN+，当1()/0.05iiiNNN+−时，表明新增加的噪声数据占到前一次检测数据的5%以上，继续迭代。当1()/0.05iiiNNN+−时，表明未标记为噪声数据的点，基本都符合高斯模型，停止检测，迭代过程结束。考虑到实现效率，在每次迭代过程中，只需要对标记为噪声的周围一定范围内的数据，建立检测窗口，进行高斯模型检测即可。3基于曲面拟合的噪声点插补算法检测出一系列的噪声点之后，需要利用噪声点周围数据的统计特性，估计噪声点处的DEM高程值。考虑到在实际地势起伏中，具有连续性和内部相关性的特点，且不会出现剧烈跳变的情况，近似符合曲面分布的，所以可以考虑以二次曲面的方式拟合。针对DEM数据，认为噪声点附近一定区域内，各点高程分布符合二次曲面函数，利用其周围正确的DEM数据进行曲面拟合，计算出拟合系数，进而确定一定区域内DEM高程值与坐标的关系。最后利用这个拟合函数系数，反求噪声点坐标处的高程值。假设在噪声点附近一定区域内，任一点的平面坐标为(,)xy，高程值为(,)Hxy，则它们之间的关系可以用二次曲面方程表达，如公式(4)所示：02212345(,)Hxyaaxayaxyaxay=+++++(4)式中，xy、分别为数据的横坐标、纵坐标，015,,,aaa…为拟合系数，(,)Hxy为点(,)xy处的实际高程值。由公式(4)可知，二次曲面拟合方程有六个系数015,,,aaa…。因此，至少需要6个正确DEM数据点，拟合系数可以通过拟合窗口内正确数据通过最小二乘原理求得。假设参与拟合的点数为n，由式(4)可列误差方程的矩阵形式221111111122222222222233333330124532211.11nnnnnnnnvHxyxyxyvHxyxyxyvHxyxyxyvHxyxyxyααααα=−……………………(5)式中n表示拟合窗口内非噪声点的个数，1,,nvv…表示拟合误差，上述方程的矩阵表达形式如下VXAH=−(6)由公式(4)分析，如果拟合区域内非噪声点个数大于6个，可以用最小二乘原理求得系数矩阵，如公式(7)所示：1()TTAXXXH−=(7)假设待拟合点的坐标为00(,)xy，利用求得的系数矩阵A，代入公式(4)就可以求出待拟合点的高程值00(,)Hxy。在实际运算中为了提高曲面拟合精度和效率，将待拟合点作为拟合窗口的中心点，并以待拟合点为原点建立相对坐标系进行最小二乘曲面拟合。4实验结果与分析如图3所示，为对图3(a)中实际某地区10001000×像素大小的DEM数据块进行迭代高斯滤波实验的结果。实验中采用的参数如下：高斯分布检测阈值2.0a=，检测窗口大小为3131×，拟合窗口大小为2121×。下表为每次迭代检测到的噪声点数目：表1实际DEM数据迭代滤波检测到的噪点数Table1NumberofnoisesdetectedonpracticeDEMdata迭代次数新增检测数共检测数目新增噪声点数比率157135713N/A2584629710.22%324065373.81%实验表明，对实际DEM数据块，通过3到4遥感技术与应用次迭代即可收敛。随着迭代过程