机载lidar点云数据的滤波

机载lidar点云数据的滤波

1基于lida的数据及热容量测量高精度数字模型（dem）在评价城市发展、土地三维可视化、道路测量等方面发挥着极其重要的作用，而高精度dem的生成是一项复杂和延迟的工作。传统的生成方法周期长、代价大,LIDAR技术的出现有效地克服了这些不足。LIDAR技术是融合了激光雷达全球定位系统(GPS)和惯导系统(INS)的技术。与传统摄影测量不同,LIDAR技术直接获取飞行区域内所有地形和地物点的水平(x,y)和垂直(z)坐标,避免了传统摄影测量中所必须的定向及影像匹配等步骤,节省了大量人力物力,为DEM快速准确生成开辟了新的有效途径。机载激光雷达在城市规划与管理、道路设计、水利工程和地形监测与分析中都有着广阔的应用前景。LIDAR点云数据是应用LIDAR技术在同一空间参考系下表达目标空间分布和目标表面特性的海量离散点的集合。它包含了建筑物、真实地面、植被等的表面三维坐标。LIDAR数据既包含地面点,又包含非地面点。在生成DEM前,必须剔除非地面点,这个处理过程称为滤波,是LIDAR数据处理中的一个重要环节。滤波是一件有挑战性的工作,机载LIDAR点云的数据量相当庞大,每平方公里可达上百万个点。随着机载LIDAR系统软硬件的更新,数据量还有继续增加的趋势,因此处理如此大的数据量需要消耗大量的存储空间和计算时间。再加上地形的多样性和地物的复杂性,滤波一直是LIDAR数据智能处理研究的难点和热点之一。2基于三种方法的滤波算法目前已有许多学者进行过相关的研究,提出了不同的滤波算法。概括来讲,主要包括基于数学形态学的滤波算法、基于坡度的滤波法、基于TIN的LIDAR点云过滤算法、基于三角网滤波算法、伪扫描线、小波分层等几种算法。2.1传统的数学形态学Zhang等提出了一种基于渐进窗口尺度的形态学滤波的方法,该方法首先将LIDAR点云数据内插为规则格网数据,再使用假定的初始结构窗口和阈值进行开运算,并计算每个点经过膨胀后的高程与原始高程之差,如果该高差小于阈值,就为地面点,否则,就为非地面点。在此基础上,逐渐增大结构窗口的尺寸,调整相应的阈值,重复上述运算,直到结构窗口的尺寸大于最大的建筑物。数学形态学理论用于LIDAR点云数据滤波时一般采用两个复合算子,即“开”算子和“关”算子。为了提取地面点,对传统的数学形态学“开”算子进行了改进,并按照下列流程进行处理:①离散点腐蚀处理。遍历LIDAR点云数据,以任意一点为中心开w×w大小的窗口,比较窗口内各点的高程,取窗口内最小高程值为腐蚀后的高程;②离散点膨胀处理。再次遍历LIDAR点云数据,对经过腐蚀后的数据用同样大小的结构窗口做膨胀。即以任意一点为中心开w×w大小的窗口,此时,用腐蚀后的高程值代替原始高程值,比较窗口内各点的高程,取窗口内最大高程值为膨胀后的高程;③地面点提取。设Zp是p点的原始高程,t为阈值,在每点膨胀操作结束时,对该点是否是地面点作出判断。如果p点膨胀后的高程值和其原始高程值Zp之差的绝对值小于或等于阈值t,则认为p点为地面点,否则为非地面点。该算法有两种滤波方式:一种是按离散点进行滤波,一种是按格网滤波。(1)根据离散点滤波是对每个激光点进行腐蚀和膨胀操作各一次,结构窗口内数据的选取按距离来量度。(2)基于格网的数学滤波方法指将每个格网看成一个“像素”,按照数字图像处理中取邻域的方法来开取结构窗口。腐蚀时,格网的“像素值”即为w×w邻域所包含格网的最小高程值;膨胀时,格网的“像素值”即为w×w邻域所包含格网的最大高程值。数学形态学滤波方法以规则格网组织数据,提出了逐离散点和逐格网滤波的方式来快速提取地面点。按离散点的方式进行滤波,该算法在格网化数据的同时,先统计出各格网的最小高程值,在逐点腐蚀(或膨胀)时,依次处理结构窗口边界所跨格网的最小(大)高程值,如果结构窗口所在边界格网的最小(大)高程值所对应的点不在结构窗口内,那么重新统计边界格网内被结构窗口所包含的数据点的最小(大)高程值,然后比较结构窗口所跨格网的最小(大)高程值,取其最小(大)值作为该点腐蚀(膨胀)后的高程值。按格网滤波的方式执行效率最高,适宜于大数据量的滤波,该算法以格网为“开”运算的对象,取其w×w邻域作为结构窗口,且各格网只腐蚀和膨胀一次,大大提高了LIDAR数据的滤波速度。按离散点滤波的方式减少了手工编辑的工作量,提高了滤波自动化处理的程度,保证了局部复杂区域滤波的可靠性。因此基于格网的数学形态学快速滤波算法,执行效率高,算法稳定,具有良好的可靠性与实用性。缺点是算法计算速度慢、耗时长,易受到滤波窗口大小的影响,且不能很好地处理地形变化剧烈的区域;若数据中存在小建筑物时,给定较小的滤波窗口即可剔除建筑物点,而过大的滤波窗口在消除大建筑物点时,也会错误地将局部地形特征剔除,从而造成误分类。总体上,数学形态学算法存在的主要问题是坡度阈值的人工选取和细节地形的方块效应。如果阈值设定太大,可能保留一些低矮的地物目标,设定太小,则可能削平地形特征。现在各种阈值的选取一般根据研究者的经验设定,或者根据地形特征设定的,没有考虑全局的特征因素,不具有普适性。解决这些问题的方法是根据地形的起伏大小和高程变化自适应的进行滤波窗口调整。但此方法在大范围地区及地形变化强烈山区的有效性还有待进一步研究。从应用上,Lindenberger将数字形态学方法引人到机载激光雷达数据滤波中,首先采用水平结构单元对机载激光测高数据进行开运算,过滤剖面式激光扫描数据,然后利用自回归过程改善了开运算结果。2.2选取偏差阀值的点和参数值基于坡度变化的滤波算法是根据地形坡度变化确定最优滤波函数,对于给定的高差值,随着两点间距离的减小,高程值大的激光脚点属于地面点的可能性就越小。假设A为原始数据集,DEM为地面点集,d是点间距离,那么满足下列滤波函数的点就是DEM的元素:DEM={pi∈A|∀pj∈A∶hpi-hpj≤Δhmax(d(pi,pj))}(1)如果对于给定点Pi,找不到临近点Pj使它们满足关系式(2),那么Pi划分为地面点。hpi-hpj>Δhmax(d(pi,pj))(2)该滤波方法主要是通过比较两点间的高差值的大小,来判断拒绝还是接收所选择的点。两点间高差的阀值定义为两点间距离的函数Δhmax(d)即所谓的滤波核函数。通常该函数是非递减函数,确定该函数的方法主要有合成函数,假设地形坡度不超过a%,且观测值没有误差,则滤波函数定义为:Δhmax(d)=ad(3)通常观测值是有误差的,所以再增加一个置信区间,并假定允许的具有标准偏差的地面点被拒绝,滤波函数就为:Δhmax(d)=ad+1.652√δ(4)Δhmax(d)=ad+1.652δ(4)在绝大多数情况下,很难用一些参数指定具体的滤波函数,因而需要根据具体的地形训练数据子集推求同地形变化特性相符的滤波核函数。这需要选择一个合适的区域作为训练数据子集用这些数据点推求Δhmax(d)。基于坡度的滤波算法具有计算简单、适应性强等特点,但是需要预先知道地形坡度和确定所开窗口的大小,所选点必须同其它所有点进行比较,以确定该点是否为地面点,也需要在整个数据集中,对每一个点进行坡度计算,这样势必造成计算量的增大,速度变慢。同时,高差阀值的选择是整个算法的关键,这些过滤阀值的设置取决于测区的实际地形状况,对于平坦地区,丘陵地区和山区,应该根据不同坡度设置不同的过滤参数值。而上述方法仅根据坡度设置统一的阈值,很可能会滤掉一些真实的地形信息,造成分类误差。要克服这些缺点可以把分块处理的思想引入,将原始点云数据按地形统计特性进行分块,然后每一个分块再按照基于坡度变化的滤波算法进行处理得到各块数据地面点集,最后根据重叠区域特征点将各块拼接,得到完整地面点集。这样不同的分块就得到不同的过滤阈值,避免了阈值的单一性,减少了分类误差。Vosselman使用Delaunay三角网组织数据,根据坡度过滤地物点的方法,通过计算该点与邻域内所有点的坡度值,如果最大坡度值在阈值内,则该点分类为地面点。Sithole在Vosselman的基础之上修改了该算法,采用变化的斜率阈值来提取地面点以适应陡坡地形的算法,不同的地区使用不同的坡度阈值来得到更好的分类结果。2.3基于点云过滤的算法描述TIN是一种重要的表示数字高程的模型,经常用来存储空间离散点之间的邻近关系。对于不规则分布的高程点,可以形式化描述为平面的一个无序点集,将该点集转换成TIN的常用方法就是构建点集的Delaunay三角网,在Delaunay三角网中的邻近点在TIN模型中也邻近。该滤波算法就是利用TIN模型中的地物临近点云高程突变关系,研究利用高差临界值条件和满足该条件的临近点数量等参数来过滤地物点。假设给定一个非空点云Pt_Cloud,并依据区域地形、建筑物、植被等分布及高程变化情况给定高差(threshold_h)和临近点数量(threshold_vn)两个域值条件,并定义Filtered和Unfiltered两个数组分别记录被过滤点和未被过滤点。首先构建Pt_Cloud的Delauney三角网,则点云过滤算法的描述如下:TIN_Filtering(Pt_Cloud,Threshold_h,Threshold_vn,Filtered,UnFiltered)如果Pt_Cloud为空,则结束;初始化Vicinity数组(邻近点集合),vn(邻近点数)和h(高差值)等参数;取出Pt_Cloud的第一个点,记为obj_pt(目标点);在Delauney三角网中检索obj_pt所有邻近点,并将邻近点逐个加入到Vicinity数组;取出Vicinity的第一个点,记为v_pt(邻近点);计算obj_pt与V,Pt的高差值并赋给h;如果h≥Threshold_h,,则vn++;如果Vicinity非空,则返回循环;如果vn≥Threshold_vn,则把obj_pt加入到Filtered,否则,obj_pt加入到UnFiltered;当进行多次(Method)循环逐步过滤点云时,算法描述如下:TIN_Method_Filtering(Pt_Cloud,Threshold_h,Threshold_vn,Method,Filtered,UnFiltered)如果Method<1,则结束;置空UnFiltered;构建Pt_Cloud的二维Delauney三角网;TIN_Filtering(Pt_Cloud,Threshold_h,Threshold_vn,Filtered,UnFiltered)Pt_Cloud=UnFiltered;Method;基于不规则三角网(TIN)的方法,是基于二维邻域搜索的方法,其计算量和算法复杂度相对较大。一般而言,由于高大建筑物和植被与其邻近地面点之间形成明显的高程突变,所以对高程突变地物,算法的过滤效果较好,但在过滤灌丛或低矮的地面物体时,产生过大误差。毛建华等对位于起伏山区的实验数据(数据包括建筑物、道路、植被、旱地和水库等覆盖类型)采用基于TIN滤波算法进行了滤波处理,结果表明植被、建筑物等高程突变比较明显地物的过滤效果较好,在所得到的DEM中,地形的起伏变化基本保持,并可以清晰地显示道路的基本形状,但在旱地平坦区域误差较大。解决上述问题的办法需要在实际的点云过滤操作中,根据具体的地形和地物及其空间分布特征选择合适的约束条件参数组合,并且采取渐进式过滤的方法才可能达到比较满意的结果,在相当多的情况下,还需要根据地面起伏的规律分割连续的地面为若干区块,并根据区块特点选择各自的参数组合。2.4局部邻域二维滤波伪扫描线是指将水平面上二维离散分布的激光点重新组织成一维线状连续分布点序列的一种数据结构,称其为伪扫描线。邬建耀等基于伪扫描线的一维邻域搜索的方法提出基于伪扫描线的滤波方法。基于伪扫描线的滤波是利用高程突变信息来区分地面点和非地面点。其基本思想是:两点之间的高度差是由自然地形的起伏和地物的高度共同引起的。若两个邻近点之间的高度差越大,那么这个高度差是由自然地形引起的可能性就越小,更为可能的是较高点位于地物上而较低点位于地面上,即:假设有两个邻近的激光脚点p1和p2,p1是地面点,p2是它的邻近点。如果它们的高度值h1和h2满足条件:h2-h1≤Δhmax×d(Δhmax是高差的容差,d是他们之间的水平距离),那么就认为p2也是地面点,否则就认为p2是非地面点。基于伪扫描线滤波算法把二维滤波问题简化为一维滤波问题,算法构造简单,有效地减少了滤波的计算量并且保证了准确性,同时该算法只需两个滤波参数,较容易实现自动化。但由于局部邻域二维滤波器大多假设邻域内高程最低点为地面点,当地面点较少的时候,这类滤波方法往往失效。而基于伪扫描线的滤波算法,总能保证每个滤波窗口中都包含有地面点,能得到比较小的一类误差和总的误差,准确地提取出地形点。在平坦地区,伪扫描线滤波效果非常好,在地形比较陡峭地区,它的误差也控制在较小范围内。但是在陡峭的斜坡和高程变化比较剧烈的区域或过滤大型物体时,为了获得可靠的结果,通常要减小高程的域值和滤波窗口的大小;在城市区域,为了全部滤除大型建筑物,则要适当增大滤波窗口,使滤波窗口的大小不小于建筑物的最大尺寸。目前,这两个参数的选取还不能做到完全的自动化,该方法还有待进一步改进。2003年国际摄影测量与遥感协会第三工作组公布了一份关于滤波方法比较的测评报告,尽管伪扫描线方法的结果并不是所有方法中最优的,但是它能得到比较小的一类误差和总的误差。2.5基于多分辨率的lida算法LIDAR激光脚点位于不同地物的点云会表现出不同的高程差异,可以借助邻近激光脚点间的高程突变来区分地面点与非地面点。激光脚点的高程突变通常采用线性预测的方法进行,假设原始LIDAR点云属于集合S,地面点集合ST可以表示为:ST={pi∈S|∀pj∈S∶hpi-hpj≤Δhmax(d(pi,pj))}(5)即在某一给定的距离范围内,若点pi与其它任意一点pj之间的高差均在该距离条件下的最大高差限差范围之内,则该点为地面点。有时,高程突变有可能是由地形本身造成的,此时,单就高程突变信息进行滤波容易产生错误。因此,可以借助方向预测法来进行滤波。方向预测法的思想:对于某一距离范围,若当前点pi与所有方向预测值的差值均大于该距离条件下的最大高差限差,则该点为地物点,否则为地面点。在以上理论支持下,为了顾及大尺度下地形的整体起伏以及小尺度下真实的地形细节,实现不同尺度LIDAR点云的滤波采用类似影像金字塔的方式构建不同分辨率的数据集,以分辨率由低到高的次序依次进行平滑处理。在局部邻域中,利用方向预测法对格网数据集进行平滑处理。如对于某一格网C(i,j),局部邻域大小取3×3,C(i,j)邻域方向的示意图如图1所示。若当前格网C(i,j)与所有方向预测值之差都小于阈值,则格网值取格网中所有点的高程平均值,否则以最小方向预测值作为当前格网值。图1C(i,j)邻域方向示意图Fig.1DirectionsofC(i,j)neighborhood对高一级分辨率的数据集进行平滑时,利用初始分辨率数据集的平滑结果来辅助新数据集的生成,然后继续采用同样的策略对该数据集进行平滑处理。基于多分辨率激光雷达滤波算法结合规则格网滤波和离散点云滤波的优点,结合使用格网数据集以及离散的LIDAR点云进行滤波,不仅能发挥规则格网快速处理的优势,而且保证了原始点云的精度。该方法在滤除地物目标点的同时,多路效应等原因引起的粗差和局外点(outliers)也得到有效的滤除,因此,无需剔除粗差的预处理。基于多分辨率方向预测的点云滤波算法可以实现复杂尺寸的地物目标的剔除处理,由于数据集在数据量上的减少,地面点提取的效率有很大提高,这种优势在处理密集的LIDAR点云时会体现得更明显。当然,该算法仍有一些地方需要进一步改进,比如对于附属在斜坡上的建筑物会出现地物滤除不完全的情况,需要结合航片等辅助数据源来提高滤波精度。多分辨率平滑处理可以得到格网数据集,但格网之间存在缝隙,降低了精度,需要利用点云数据和格网数据集进一步做点云滤波处理。万幼川等利用此滤波方法,对青岛地区的城市激光点云进行了分类实验。数据由LeicaALS50仪器采集,点云的平均密度为1.047点/m2,区点云包含多种地形特征和复杂地物如空旷裸地、植被、尺寸不一的建筑物、道路等。实验结果表明,本算法在正确滤除各种复杂的地物对象的同时,较好地保留了地面上的不连续特征。2.6小波分层滤波小波变换是介于函数的空间域表示和频率域表示之间的一种表示方法。它在空间域上和频率域上同时具有良好的局部化性质,对高频成份采用逐步精细的空间域取样步长,可以“聚焦”到对象的任意细节,从而被誉为“数学显微镜”。小波变换是信号处理中一个经典的理论,它的基本思想是对原始信号建立数据金字塔,从而获得不同尺度上的信号描述。在最大尺度的信号描述(即金字塔的顶层)中获取最初的兴趣域,然后把这个粗略的兴趣域传到下一层中作为当前层兴趣域的初始值,从而减少计算的时间,提高处理结果的精度。把小波变换引入到机载LIDAR数据滤波中来,实现了基于小波变换的数据滤波。在一般情况下,激光雷达的点云数据的地面点表现为局部区域内的高程最低点。因此,可以用特定大小的窗口分割原始数据,然后在每个窗口中选择一个高程最低的点,组成一个新的数据描述。对这些地面点进行组网,从而形成一个粗略的地形表面。然后利用这个粗略的地形表面作为参考面,在下一层进行滤波,获取更多的地面点。必须保证每一个分割窗口中至少有一地面点,它需要分割窗口足够大。采用比目标区域内最大的建筑面积稍微大一点的窗口作为最上面一层数据描述的尺度。接下来确定第一层数据描述的窗口尺度。以目标区域内最小的人造建筑的面积为第一层数据描述的窗口尺度。这样金字塔的层数可以确定下来了,金字塔的层数可以用下面的公式来描述:n=[InSw/Ink]+1(6)Sw为最上层数据描述的窗口尺度;k为第一层数据描述的窗口尺度;n为金字塔的层数。描述点与参考面的相对关系最简单的尺度就是高程差异,当高程差异大于某个阈值时,就认为该点不属于地面点;而当高程差异小于某个阈值时,就认为该点属于地面点。小波分层滤波算法要先进行分层,然后把这个粗略的兴趣域传到下一层中作为当前层兴趣域的初始值,从而减少了计算的时间,提高处理结果的精度,但是每一层的判断结果受到了上一层的影响,如果上一层次的处理出现了错误,这种错误会导致下一层次的数据点类型判断出现错误。另外分割窗口的尺度选择也很重要,最上面一层数据描述的尺度要选择比目标区域内最大的建筑面积稍微大一点的尺度,第一层数据描述的窗口尺度要选择目标区域内最小的人造建筑的面积的大小。还有小波分层滤波算法还需要在数据初值选择和判别规则方面考虑更加细致,剔除数据中的粗差。王刃等用英国Birmingham地区的机载LIDAR数据进行基于小波方法分层的滤波验证,结果几乎所有的非地面点都被有效地过滤掉了,而所有的地面点被有效地保留下来,从而为获取高精度的DEM奠定了基础。地面点、非地面点的误判率低,正确率优于94%,精度高于±0.71m。此外,机载激光雷达点云数据滤波算法还有激光雷达点云平面拟合过滤算法,基于距离限制的激光数据滤波算法,基于区域生长的LIDAR点云数据滤波算法等。3基于多分辨率激光雷达滤波算法的分析与展望(1)在平坦的地区各种