基于不同滤波方法的机载激光雷达数据处理及表面模型构建1)

刘 芳 张 琼 范文义 陈 成 杨树文 李 典

(东北林业大学，哈尔滨，150040)

激光雷达(LIDAR)是一种通过位置、距离、角度等观测数据直接获取对象表面点三维坐标、回波和强度信息，实现地表信息提取和三维场景重建的对地观测技术［1］。激光雷达技术在对地探测方面有强大的优势，具有空间与时间分辨率高、动态探测范围大、能够部分穿越树林遮挡、直接获取真实地表的高精度三维信息等特点。这使得激光雷达可以快速高精度地获取地表物体的水平和垂直结构信息。机载激光雷达作为一种新兴的获取信息的技术，因其独特的优势而被广泛应用于各个不同的领域。

激光雷达数据处理是获取信息的前提，而激光雷达点云数据的滤波和分类则是数据后处理的主要任务之一。利用激光雷达点云数据可以获取反映地形起伏变化的数字高程模型(DEM)和目标区域的表面三维信息——数字表面模型(DSM)以及数字城市模型等。从原始点云数据中获取DEM和DSM这两种表面模型最关键的步骤是如何有效地、可靠地滤除噪声点，并分别提取出仅反映地形表面信息的地面点和只反映地物表面信息的地物点，这个过程即滤波分类。笔者综合考虑研究区特有的地形、植被等因素，对凉水国家级自然保护区不同类型的区域采用不同的算法进行滤波试验，取得了较好的效果。

1 研究区概况及数据获取

研究区域为黑龙江省凉水国家级自然保护区，位于北纬 47°6'49″～ 47°16'10″，东经 128°47'8″～128°57'19″。保护区全境为山地，总面积6394 hm2，属小兴安岭南部达里带岭支脉东坡。

LIDAR数据于2009年9月采用LiteMapper 5600激光雷达系统获得。激光扫描仪为Riegl LMSQ560，其波长为1550 nm，回波宽度分辨率0.15 m，垂直精度可达0.15 m。采样点云密度大于2点/m2，飞行高度达800 m，平均对地飞行速度是180 km/h。此外同步获取了空间分辨率为0.2 m的CCD影像共189幅，总覆盖面积约为1.2 km×1.6 km。LIDAR数据保存成 LAS格式，投影方式为UTM，参考椭球为WGS84，每个激光点包含了激光点三维坐标值、强度值、回波类型等信息。

2 数据处理

选择了3种最为普遍有效的滤波方法对研究区进行LIDAR数据滤波处理。

2.1 滤波原理

①不规则三角网算法。Axelsson提出了一种渐进式不规则三角网加密算法。其基本原理是:首先选取区域内的少量低点作为三角网的种子点，然后将点云中的点逐步加入到三角网中，如果要加入的点到三角形的距离(d)以及由于加入该点所形成的新三角形的内角(α、β、γ)均小于所设定的阈值，则该点为地面点。在每轮加入点运算前，都须通过运算重新生成阈值。当没有点可再加入三角网时(即没有点符合阈值标准时)，运算结束，形成最终的三角网［2］。这种方法通过反复迭代不断向模型中添加新点，每个新加的点都使模型更加接近实际地表面。算法如图1所示。

图1 不规则三角网加密算法

图2 回波垂直剖面示意图

②多重回波算法。大部分激光雷达系统除了记录高度信息外也记录了回波及回波的强度信息，打在不同地物上的激光脉冲回波能够反映不同的信息。对于植被，一个激光脉冲可以穿过树叶、枝干到达地面，得到不同的回波信号，即返回不同的三维点信息［3］。激光脉冲经发射后到达地面的过程中被不同高度的植被返回，首次回波一般是植被冠层表面返回的信息，而中间表面则包括第二次和中间次回波，地面表面则往往是由最后一次回波返回的点构成的。多重回波过滤就是利用回波数目及回波次数的特性对点云进行过滤，由最后一次回波可以得到DEM，而由首次和第二次回波能得到DSM。图2是由不同高度的植被冠层反射回来的LIDAR激光脉冲信息［4］。

③迭代线性最小二乘内插法。此法最初由奥地利维也纳大学的Krasu和Pfeifer等人提出［5］。其核心思想是基于地物点的高程比对应区域地形表面激光脚点的高程大，经过线性最小二乘内插后，激光脚点高程拟合残差相对于拟合后的地形参考面不服从正态分布(见图3)，高出地面的地物点高程拟合残差都为正值且数值较大，而地面点的拟合残差较小且可能为负值［6］。该方法首先用所有激光脚点的高程观测值按等权计算出初步曲面模型(该曲面介于真实地面和植被覆盖面之间，其结果使拟合后真实地面脚点的残差为负值的概率变大)。然后依据每一个激光点到该新生成的表面的距离和方向利用稳健估计的权函数关系式(Ⅰ)来计算其权值p:

式中:v为高程值与拟合面的残差;p表示内插时使用的权值;a和b决定内插权函数的陡峭程度;g决定着哪种点的权值为1;地上偏移参数w指决定某一点对中间表面是否有影响的上限值。经过几次迭代后最终的裸露地表点从中间表面中抽取出来，所有高程值满足权函数关系式(1)中第一二种情况的点都被作为裸露地面点。然后用这些点进行内插处理生成最终的DEM。

图3 定权示意图

图4 去除明显高点及孤立点

2.2 滤波试验

为了说明这些方法的可行性，分别选择了研究区内不同区块的点云数据进行处理，这些区域既包含裸露地、水域，也有成片的树林及房屋等地物。为了有效地过滤原始数据，使所获取的DEM更加符合实际区域地形，首先对原始数据进行粗差点剔除，包括剔除明显的低点和孤立点以及空中悬浮点，再进行滤波分类试验。

图5 消除噪声点后按高程显示的数据

①TIN三角网算法滤波效果。图6是用于检验TIN三角网算法滤波效果的原始点云。该区域包含植被、河流湖泊、房屋及道路。首先根据试验区地形选择不同的最陡坡度、最大内插角度和最大内插距离。由于该区域比较平坦，上述参数分别设置为45°、55°、1 m时效果最好，设置最大建筑物为30 m，分离出地面点后便得到图7的结果。对比原始点云及同步获取的影像发现几乎所有的地面点被有效地保留下来。然后利用Terrascan中的其他分类工具分出植被和建筑物点，这些过程都需要根据该区域的特点选择合适的参数，才能得到比较理想的结果。

图6 原始点云

图7 滤波后获得的地面点

②多重回波算法滤波效果。笔者选取试验区中包含植被和人工建筑物的一片区域进行滤波分类。图8为根据不同回波类型记录值分离出来的首次激光脉冲脚点，对比首尾回波可以发现直接或穿透植被打在地面上的点比较少，不利于构建地面模型;图9为尾次激光脉冲脚点;对比两图可以得出如下结论:对于森林地区，大部分打在植被冠层上的激光脉冲束不是直接返回而是向下继续穿透直至抵达地面。但是能够到地面的激光脚点数量还是不多。此外对于人工建筑物特别是屋顶，其上的激光脚点则往往只有一次回波，所以这部分脚点既属于首次回波也属于尾次回波。

图8 首次回波点云

图9 尾次回波点云

③迭代线性最小二乘内插法效果。笔者选取了图10左侧区域(均为植被)进行滤波实验，设置参数 a、b、g、w 分别为0.5、2、0、0.5，迭代 8 次，并设置滤波的格网大小为3 m，得到的结果如图10右侧所示，其中右上部分为植被部分的点云，而右下部分为过滤出的地面点云。

图10 迭代线性最小二乘内插法分类示意图

图10从激光点云角度反映了这一算法在提取大型建筑物及低矮灌木植被的缺陷。借助激光雷达同步获取的影像发现图11对应的试验区有大面积建筑和低矮灌丛，对于该试验区，笔者试了不同的a、b、g和w参数，但是效果一直不理想，总是有残存的植被点和建筑物。

3 模型生成

采用上述不同算法滤波后笔者将各自的滤波结果保存成ASCII文件，再生成相应的数字高程模型与数字表面模型。

图11 迭代线性最小二乘内插算法错分的地面点

3.1 由TIN滤波结果构建

图12即为用TIN滤波所得地面点构建的DEM，与传统的DEM相比，由该方法构建的地面模型更能反映地表的实际情况，微小的地形起伏很好地被展现出来。而图13为DEM渲染及等高线叠加显示图，充分反映了TIN三角网滤波算法的适用范围，此图右半部分地形比较平坦，而左半部分是起伏坡度较大的山，可以看出平坦地区的等高线特别细碎，构建的DEM也很破碎;而用左边陡峭地区滤出的地面点生成的等高线平缓又光滑，证明该部分DEM质量很好。

图12 TIN三角网构建的DEM

图13 平坦与陡峭地区DEM对比

3.2 由多重回波滤波结果构建

图14和图15分别是首回波和尾回波生成的部分DSM，对比可看出首回波蕴含的植被信息。图16为首回波脉冲脚点的高程同尾次回波脉冲脚点的高程作差后形成的正规化数字表面模型(nDSM);在森林地区nDSM被称作冠层高度模型(CHM)，它反映的是冠层的高度信息。可看出这个模型几乎呈平面，没有山体那种大的起伏或坡度，因此正规化的表面模型，消除了传统表面模型中地形起伏变化对地物高程及形状的干扰，能获取更准确的地物形态和高度信息。

图14 首回波生成的DSM

图15 尾回波生成的DSM

图16 首尾两次回波差生成的nDSM

3.3 由迭代线性最小二乘内插结果构建

图17和图18分别为利用迭代线性最小二乘法过滤出的地面点与非地面点生成的DEM及DSM。可知，对于植被生长状况良好的林区，迭代线性预测算法的效果很好，几乎所有地面点都能被分离出来，故生成的DEM能很好地反应微小的地形起伏。图19反映了该算法在提取大型建筑物及低矮灌木植被的不足。

图17 迭代线性法分离出的地面点生成的DEM

图18 迭代线性法分离出的植被点生成的DSM

图19 迭代线性法失效生成的DEM

4 效果评价

激光雷达点云滤波分类之后还需要进行滤波效果评价，评价分为定性评价与定量评价两种。定量分析是统计滤波结果中有多少非地面点被误分为地面点，有多少地面点被误分为非地面点，从而得出TypeⅠ型错误TypeⅡ型错误的百分比［7］，然而对不同地区的海量数据进行统计往往是很困难的。因此本研究仅从定性的角度评价3种算法的滤波质量，具体为:可视化分析(利用分出的地面点在ArcGIS中生成渲染图，通过缩放观测表面模型的构建状况)、等高线分析(根据等高线的疏密、平缓、光滑程度等特性来判断滤波结果好坏)、剖面图分析(直接针对滤波后的点云进行分析，从立面评价滤波效果，判断点是否错分)。

5 结论

采取这些手段评价分析滤波分类效果后，得出或验证了各方法的适用范围及优缺点，结论如下:①对高程突变地物，TIN不规则三角网迭代算法的过滤效果较好。这是由于高大建筑物和植被与其邻近地面点之间形成了明显的高程突变，但在过滤灌丛或低矮的地面物体时，产生过大误差。②利用激光雷达的不同地物对应回波的高程差不同进行滤波的多重回波算法简单有效，根据植被点和与其他点的回波类型情况不同能够很快地把植被点与其他激光脚点分开，特别适用于分类植被激光脚点;但由于穿透率的问题，末回波点也即地面点稀少，因此利用末回波构建的DEM不能很好地反映出精细地形。③迭代线性最小二乘内插算法借助邻近激光脚点间的高程突变来区分地面点与非地面点，能很好地获得地形趋势面。但是该方法在有大型建筑物存在的区域往往不适用，它仅能削去房屋的棱角，而无法完全滤掉建筑物。对地形起伏变化不大的森林地区，其效果则是其他方法无法比拟的。④笔者根据凉水自然保护区特有的地形、植被等因素进行分块，然后对每一分块采用适合其地形特点的滤波算法进行处理，所得效果较佳。对平坦地区设置不同的参数采用迭代线性最小二乘内插法进行处理，对有许多高大植被和建筑物的地区采用不规则三角网算法滤波，对森林地区采用多重回波算法分离出地面点集，再利用各小分块间重叠区域的特征进行拼接得到凉水地区完整的地面点集并生成DEM。

［1］周淑芳，李增元，等.基于机载激光雷达数据的DEM获取及应用［J］.遥感技术与应用，2007，22(3):356－360.

［2］曾齐红.机载激光雷达点云数据处理与建筑物三维重建［D］.上海:上海大学，2009.

［3］张小红.利用机载LIDAR双次回波高程之差分类激光脚点［J］.测绘科学，2006，31(4):48－50.

［4］Akay A E，Oˇguz H，Karas I R，et al.Using LiDAR Technology in Forestry Activities［J］.Environ Monit Assess，2009，151:117－125.

［5］Krasu K，Pfeifer N.Determination of Terrain Models in Wooded Areas with Airborne Laser Scanner Data［J］.ISPRS Journal of Photogrammertry and Remote Sensing，1998，53(4):193－203.

［6］Pfeifer N，Reiter T，Briese C，et al.Interpolation of High Quality Ground Models from Laser Scanner Data in Forested Areas［J］.International Archives of Photogrammertry and Remote Sensing，1999，32(3/W14):31－36.

［7］Sithole G，Vosselman G.Experimental Comparison of Filter Algorithms for Bare－Earth Extraction from Aiborne Laser Scanning Point Clouds［J］.ISPRS Journal of Photogrammertry ＆ Remote Sensing，2004，59:85－101.