机载激光雷达(LiDAR)作为一种新型的主动式遥感技术，能直接获取高精度的地面三维坐标，即三维激光点云数据^[1]。对于森林区域，机载LiDAR技术能部分穿透森林植被，全天时、全天候获取地物信息，具有光学遥感影像难以比拟的优势^[2-3]。

随着激光雷达扫描技术日益成熟，国内外学者对LiDAR点云数据的研究逐渐深入，然而如何获取高精度的数字高程模型(DEM)依然是一项难题，其关键在于点云的滤波方法，即区分地面点与非地面点^[4-6]。对于森林区域点云数据来说，大部分数据为植被点云，地面点与林下植被点均处在植被点下方，裸露情况较少，且山地地形与城市区域平坦地形间存在较大差异^[7]，这些都增加了滤波难度。

点云滤波作为LiDAR点云数据应用的基础工作之一，目前已有众多研究，形成了数学形态学方法、基于不规则三角网法、基于坡度算法和移动面拟合法等算法^[8-9]。例如, 基于数学形态学的滤波算法，通过找出结构元中高程值最低点云，并在此基础上对结构元中所有点云进行腐蚀、膨胀运算，从而去掉小于结构元尺寸的突出地物，达到滤除地物点的目的。在该算法研究中，Weidner等^[10]首先提出了一种基于点云数据灰度值的形态学滤波算法，随后Wack等^[6]提出使用不同尺寸的窗口迭代滤波的方法，一定程度上对基于可改变窗口下的滤波算法奠定了基础。此后，Hu等^[12]、张熠斌等^[12]、周晓明^[13]、王芃芃^[14]等的研究基于不同地物环境条件对算法进行了一定改进。基于坡度的滤波算法由Vosselman^[15]提出，它是利用地形中坡度的变化，凭借点云的坡度值大小来分离区域内地面点与地物点的方法，并由Lee等^[16]、Sithole^[17]在森林地区地形应用领域进行了一定改进。

近年来，我国研究者张小红^[18]首次提出了点云数据的移动曲面拟合法滤波算法，该算法在初始种子面的基础上不断判断临近点云是否属于种子面，并不断更新拟合平面，从而实现滤波目的。此外，还有陈磊等^[19]的基于区域生长算法，马鼎等^[20]的利用融合回波强度的分块拟合法等。同时，国际学者更加关注无需参数设置的自动化滤波算法，例如Mongus等^[21]参考点云标准差信息自动设置分层参数，Pingel等^[22]借助图像信息计算最适窗口大小等。

然而，目前已有的方法还存在着一些不足，例如对于传统的数学形态学滤波算法而言，在点云的腐蚀与膨胀过程中，通常改变了原始点云数据的高程值^[23]，导致利用结构元搜寻点云过程中，根据结构元尺寸大小的不同，对不同地形容易造成地形平滑或无法准确识别地物的结果^[24]。再如移动面拟合法，其精度较高，但每加入一点时均需要重新拟合平面，在迭代过程中，运算量极大^{[18, 25]}。因此，本文针对以上问题，采用八叉树结构作为滤波准则。该结构仅用于管理数据，保证了原始点云数据的精度不被更改，使地形特征得以保留，并且可以根据点云位置的不同，将大量点云分割^[26]，控制每个节点中点云数量，在保证结果精度的基础上，减少地面平面拟合时的运算量。

对于森林区域地形起伏较大的现实情况^[17]，利用八叉树结构对点云分割，可以将较大的山包或山洼地形区域分割成为多个子区域，使每一个子区域都近似呈现为山坡地形。由此，可将复杂地形分解并统一转换为大量山坡地形，一定程度上降低了复杂地形对滤波算法的影响，提高了滤波算法的适用性，并简化了滤波运算。

本文利用甘肃省张掖市大野口森林区域的机载LiDAR点云数据，改进传统八叉树模型的分割原理，在分割过程中依据点云数据中心点缩小节点尺寸，提高同类点云的分割准确率，并基于坡度判断的方法，对森林区域不同地形下的点云数据进行滤波处理，为后续森林参数的提取奠定基础。

1.   研究区概况

研究区位于甘肃省张掖市祁连山中段西水林场大野口流域，地势西北高东南低，坡度在20°~40°左右。地貌类型丰富，由多条西北—东南走向的平行山脉和宽谷组成。主要的气候类型为大陆性高寒半湿润山地气候，气温年较差大，冬季寒冷干燥，夏季温凉湿润，全年降水集中在5—9月份，年降水量在400 mm左右。土壤类型具有明显的垂直地带性特征，主要有山地灰钙土、山地灰褐土、山地栗钙土、亚高山灌丛草甸土和高山寒漠土。研究区主要树种有分布在阴坡、半阴坡的青海云杉(Picea crassifolia)、山杨(Populus davidiana)和分布在阳坡、半阳坡的祁连圆柏(Sabina przewalskii)等^[27-28]。

LiDAR数据获取于2008年，来源于黑河计划数据管理中心，飞行区域为大野口森林，飞行航线方位角为北偏东115°，飞行高度距地面约760 m。采用LiteMapper 5600的LiDAR系统，使用双频DGPS和惯性导航系统IMU进行定位。其中激光扫描仪为Riegl LMS-Q560，波长是1 550 nm，激光脉冲长度为3.5 ns，激光脉冲发散角小于等于0.5 mrad，点云密度平均约为3.43个/m²，属高密度点云数据^[29]。

2.   研究方法

2.1   八叉树模型

2.1.1   八叉树基本原理

在三维空间中，八叉树结构将整个数据空间分割为大小不等的立方体以存储各类数据^[30]。该结构利用一个立方体包含所有数据，称为根节点。根据数据的类型，可设定不同的数据分割原则，将节点分割为8个相同大小的立方体，即子节点^[31]。八叉树结构示意图见图 1，其中左图为三维空间中分割示意图，右图为对应的节点示意图。

图  1  八叉树示意图

Figure  1.  Schematic diagram of octree

下载: 全尺寸图片 幻灯片

2.1.2   八叉树分割原理

八叉树的分割由根节点开始，分割出8个子节点。若子节点中不包含任何数据，将该节点称为空节点；根据所设定的分割原则，若子节点中所有数据达到一定要求(如数据类型相似或均满足一定条件等)，则不再继续分割，将该节点称为黑节点；否则，将该节点称为灰节点^[32]，并根据分割原则，将该节点进一步分割为8个子节点。根据上述分割步骤进行迭代，迭代停止的条件有如下两种：(1)所有节点均为空节点或黑节点；(2)达到所设定的迭代次数，即最小节点已符合预先设定节点的最小尺寸。

2.1.3   改进的八叉树分割原理

传统的八叉树依据节点尺寸，将节点分割为大小相等的8个子节点，但在分割过程中，或存在部分冗余的空节点，或某些灰节点内只存在少量数据，并可能将同一类型的点云分割到不同节点中，从而增加不必要的计算量。本研究中根据节点中点云的数量、位置，在分割过程中，以点云数据中心点作为节点中心，对节点尺寸予以调整，即每次分割时，均采用最小立方体包含节点中所有点云数据后再进行分割。该方法在不影响节点父子关系及相邻关系的情况下，一定程度上减少了分割的次数及空节点的数量，同时可以提高将地面点与地物点分至不同节点中的概率。

2.1.4   八叉树模型建立

八叉树在存储三维空间数据、管理数据所在位置上拥有独特的优势，以不破坏原始数据为基础，让用户能迅速找到每一个数据所处的位置。其基本原理是利用数字对每一个节点进行编号。具体规则^[33]如下：

(1) 令根节点编号为0。

(2) 将根节点分割成8个子节点，分别编号为1、2、3、4、5、6、7、8。

(3) 若节点为空节点或黑节点，则不再进行分割编号。

(4) 若节点为灰节点(如图 2中节点3)，则继续分割，其8个子节点分别编号为31、32、33、34、35、36、37、38。

图  2  节点编号示意图

Figure  2.  Schematic diagram of node No.

下载: 全尺寸图片 幻灯片

(5) 若节点35依然为灰节点，则其8个子节点分别编号为351、352、353、354、355、356、357、358，以此类推，直到所有子节编号完毕(图 2)。

由此，区域内所有数据均存在于某一个编号的节点中，通过节点编号，可以迅速找出该节点的父节点与子节点以及相邻节点，以便准确快速地定位每一点所在的位置及显示邻域内存在的其他点云^[34]。

2.2   点云滤波方法

原始机载LiDAR点云数据通常杂乱无章，利用点云滤波算法，可将大量离散点云依据算法分离为地面点与地物点。其分离结果可拟合得到区域内DEM^[35]。

本研究采用改进的八叉树模型，将大量点云数据分割到不同节点中，使每个节点中只包含几个到十几个点云以简化数据处理，在此基础上进行点云滤波。具体滤波步骤如下：

(1) 根据改进的八叉树模型，将所有点云依据其所处位置，分配到不同的节点中，并对每个节点编号。

(2) 根据节点编号找出所有最小节点，对其中高程最小的点所在节点，将其高程值最小点视为地面点，并与一定阈值内的邻近点根据最小二乘法进行平面拟合，得到合成初始种子面aX+bY+cZ+d=0。所包含的点均视为地面点。

(3) 对初始面根据节点大小建立的三维坐标系，引入坡度判断^{[15, 36]}，对于初始面所在节点中未分类的点，分别计算每个点到初始面的距离以及坡度，对距离及坡度与设定阈值作比较，当满足条件时，将该点划入地面点集；否则，划入地物点集。

(4) 根据新加入的地面点，对初始面重新拟合，直至节点内左右点数据分类完毕，得到初始拟合平面。

(5) 重复步骤(2)~(4)，将所有节点中的点进行分类，并拟合出若干平面，根据节点编号，将每一个平面与其邻近平面进行合并。合并规则包括：平面之间角度大小不超过设定阈值；当平面角度小于阈值时，平面之间距离不超过设定阈值。每一次平面合并后均需要重新拟合平面。以此类推，对于无法与其他平面合并的平面，将其中点云数据划入地物点。

(6) 根据新合成的平面集，依据每个平面的平均高程值与相邻平面作比较，若高程差大于设定阈值，则将平均高程较高的平面上的点划入地物点。

当节点中存在点云数据较少以至无法拟合的平面时，需要根据节点编号，将该节点中点云与其相邻节点中点云合并，达到平面拟合所需要的数据量。

在步骤(3)中，引入了坡度判断，其原因在于森林区域地形变化较大，山坡上较低的植被点数据与山坡所在平面的直线距离可能较小，若只根据点到平面的距离作为判断依据，则容易将山坡上的较低的植被点划入地面点，引起误差。其具体规则为：

(1) 若点k到初始面距离大于阈值，则视为地物点。

(2) 若点k到初始面距离小于阈值，则进行坡度判断^[37]。其公式为：

式中：G_k为点k所对应的坡度值，H_k为点k的高程值，H_c为初始面中最新加入的一点c的高程，X_k、X_c、Y_k、Y_c分别为点k与点c的水平坐标。

当G_k值小于坡度阈值时，将点k视为地面点，否则剔除至地物点集。其中，该坡度阈值应根据区域内地形平均坡度设定，阈值取值接近于区域内平均坡度。

3.   结果与分析

3.1   基于八叉树原理的点云分割

本研究使用Matlab编程语言，基于八叉树原理，对研究区域内点云数据进行分割，选取小块区域样本对改进的八叉树模型进行分析，其结果如图 3、图 4所示(图中坐标均为点云的空间位置坐标)。

图  3  传统八叉树模型下的点云分割

图中坐标均为点云的空间位置坐标。下同。

Figure  3.  Point cloud segmentation based on conventional octree model

The coordinates in the figure are all spatial position coordinates of point cloud. The same below.

下载: 全尺寸图片 幻灯片

图  4  改进八叉树模型下的点云分割

Figure  4.  Point cloud segmentation based on improved octree model

下载: 全尺寸图片 幻灯片

图 3为传统八叉树模型下点云数据分割结果。其中，图 3a为三维空间中的分割示意图，图 3b为经过水平投影映射到二维空间后的分割示意图。在传统八叉树模型下，依据点云空间位置，被均匀分割到8个子节点中，相邻节点边界及节点内点云数据被赋予不同颜色以便区分。由图 3b可以看出，均匀分割下，大部分植被点与地面点被分割到不同节点中。为便于阐述，令图 3b中上方4个节点为植被点节点，下方4个节点为地面点节点。在本研究所选样例中，植被点节点中，不存在任何地面点，表明八叉树模型在点云分割中具有一定的可行性；而在地面点节点中，还存在少数的植被点，表明传统的八叉树模型还存在一定的改进空间。

图 4为改进八叉树模型下点云数据的分割结果。改进后的模型在分割点云数据时不再仅仅依据节点尺寸大小进行平均分割，而是根据节点中点云数据的位置，缩小节点尺寸，以点云数据聚集的中心作为基准点加以分割。对比图 3b与图 4b可以看出，两幅图中存在分割结果差异的点云数据已用红色椭圆圈出。在传统模型中被分割到地面点节点的植被点，在新的模型下被正确分割到植被点节点中，一定程度上提高了后续滤波算法的精度，也为地面点平面拟合节省了部分时间。但是，在改进八叉树模型结果中，部分地面点节点中还存在一些植被点，这一点不可避免。虽然限定节点中点云数量可以尽量分离地面点与植被点，但当节点过小，容易造成地形特征支离破碎，反而破坏了原始地形。

总体上，改进的八叉树模型较传统模型对于点云的分割更科学，效果更好，虽然在此例中仅有几个点的分割区别，但对于较大区域内数十万甚至上百万的点云数据，新的模型将会发挥重要作用。

3.2   点云滤波实验

基于八叉树原理，对研究区域点云数据依据所述滤波方法进行点云滤波。实验中以3种不同地形下数据样本为例，对滤波结果进行分析。

数据1：实验结果如图 5所示。图 5a为原始点云数据，可以看出此区域内地形起伏不大，其中山包与山洼较少；有一定坡度但坡度较为平缓，其中坡度最大处约为19.5°，坡度最小处近似平坦地面。

图  5  数据1滤波实验结果

Figure  5.  Results of filtering in data 1

下载: 全尺寸图片 幻灯片

图 5b为滤波后地面点数据内插后所得到的格网地面高程模型(DEM)。从图 5b可以看出：边缘部分存在少许齿状格网，其原因产生于实验区域的截取，但对结果的影响可忽略不计。由图 5b可以看出，在左下方的山坡处，存在少许地物点被错误划分为地面点的情况，致使拟合模型中存在一定突起。根据拟合的DEM结果，结合原始点云数据，分类错误的地物点云主要聚集在两处，具体位置在空间三维坐标系中分别处在坐标点(609 052，4 265 847，2 785)与坐标点(609 068，4 265 856，2 781)以及两点邻域内(图中红色椭圆标记处)。在点A邻域内，共存在11个地物点分类错误；在点B邻域内，存在8个地物点分类错误。

分类错误产生的原因在于，两点坐标邻域内，存在一定面积的低矮植被(苔藓)，这一点通过比较图 5a与图 5b中红色椭圆标记处可以看出：图 5b中点云滤波错误的区域，在图 5a中相应的植被点云也较少。这就导致该区域内所获取的点云数据中，存在较多的苔藓点，致使在八叉树的分割中，可能被划分到同一节点或相邻节点中，而这些点云所拟合的平面与真实地面点平面差异较小，使得最终平面合并中出现误差。

对于其他植被点云较为密集的区域，苔藓点被扫描仪获取的概率较小，所以这部分区域滤波效果更好，具体滤波结果错误率见表 1。

表  1  数据1滤波精度验证

Table  1.  Precision verification of filtering in data 1

点云数量 Number of point cloud	第Ⅰ类错误率 Error rate of type Ⅰ	第Ⅱ类错误率 Error rate of type Ⅱ	总错误率 Total error rate
418	6.40%	2.62%	4.57%
注：第Ⅰ类错误率是指将地面点误分类为地物点的错误率；第Ⅱ类错误率是指将地物点分类为地面点的错误率；以目视解译的分类结果作为验证计算错误率。Notes: error rate of typeⅠ: the misclassification rate of classifying ground points into object points; error rate of typeⅡ: the misclassification rate of classifying object points into ground points; classification result of visual interpretation is used to verify the error rate of calculation.

下载: 导出CSV 

| 显示表格

根据滤波结果精度验证，第Ⅱ类错误率较小，表明地面点中存在的地物点数量较小，所拟合的DEM基本反映了样地区域的真实地形；第Ⅰ类错误率较大，表明地物点中存在一些真实地面点，但这一点对于DEM的拟合结果影响不大^[38]。

数据2：数据2实验结果如图 6所示。从图 6a中原始点云数据可以看出：此区域内存在陡坡及少许山洼，其中坡度最大处约为32.2°，坡度最小处约为7.9°。从图 6b中滤波后所得格网DEM可以看出：在图中偏左的山坡处，拟合地形较为平滑；但在坡度减缓处(图中红色椭圆标记处)，存在部分地物点被错误划分为地面点，使拟合模型中凹凸不平。根据拟合的DEM结果，结合原始点云数据，分类错误的地物点主要聚集在坐标点(608 958，4 265 854，2 790)及其邻域内，共存在23个地物点分类错误。

图  6  数据2滤波实验结果

Figure  6.  Results of filtering in data 2

下载: 全尺寸图片 幻灯片

分类错误产生的原因在于，该点坐标邻域内存在坡度急剧变化的情况。这一点在图 6b中红色椭圆标记处可以看出。该区域内所获取的点云数据中，存在一定数量的林下植被点，而在八叉树的分割结果中，该区域节点内初始拟合平面坡度较大，由于坡度急剧变小，平坦地形处部分较低的林下植被点被误认为山坡地面点，从而导致这些地物点被错误划分至地面点集。

对于图 6a与图 6b中绿色椭圆所标记的区域，从图 6a可以看出：该区域存在部分山洼地形，但没有点云数据存在，从而导致DEM内插过程中山洼地形出现平滑现象而消失，这一点与滤波结果精度无关。具体滤波结果错误率见表 2。根据滤波结果精度验证，其中第Ⅱ类错误率相比数据1有一定增加，表明滤波算法对于坡度存在急剧变化的区域还存在一定错误率。第Ⅰ类错误率相比数据1有所下降，表明在较为陡峭的山坡地形中，地面点被划分为地物点的概率更小。两个数据样本的滤波总错误率较为接近，一定程度上表明该滤波算法可以适用于不同地形区域，结果较为稳定。

表  2  数据2滤波精度验证

Table  2.  Precision verification of filtering in data 2

点云数量 Number of point cloud	第Ⅰ类错误率 Error rate of typeⅠ	第Ⅱ类错误率 Error rate of typeⅡ	总错误率 Total error rate
4 191	5.70%	3.35%	4.75%

下载: 导出CSV 

| 显示表格

数据3：数据3实验结果如图 7所示。从图 7a中原始点云数据可以看出：此区域内存在部分山包与山洼，同时地形坡度与地形起伏较大，其中坡度最大处约为34.3°，坡度最小处约为28.1°。从图 7b中滤波后所得格网DEM可以看出：图中偏左的山坡处，存在两处突起(图中红色椭圆标记处)，但与前两个数据的结果相比存在一定的不同。由图 7a可以看出：两处红色椭圆标记处地形较为相似，即山包地形周围存在山洼地形，山洼地形中存在点云空白区域，但两处区域滤波结果有所不同，其中上方区域中，在保留了山包地形的同时，还存在部分地物点错误划分的情况，分类错误的地物点主要聚集在坐标点(608 929，4 265 859，2 781)及其邻域内，共存在11个地物点分类错误；而下方区域中，并不存在地物点分类错误的情况。在DEM图像中微微突起的原因在于，该区域为山包地形，但由于被空白点云的山洼地形所包围，其周围地形被填平，所以在三维显示中，突起部分较小。

图  7  数据3滤波实验结果

Figure  7.  Results of filtering in data 3

下载: 全尺寸图片 幻灯片

在图 7b中偏右处，同样存在山洼地形平滑现象(图中绿色椭圆标记处)，但与其他山洼处无点云数据的情况不同，图中所圈出的区域存在部分低矮植被点，所以该区域被填平的原因在于部分地物点的分类错误而非DEM内插所造成。分类错误的地物点主要聚集在坐标点(608 909，4 265 757，2 790)及其邻域内，共存在8个地物点分类错误。

分类错误产生的原因在于，八叉树分割过程中设定了节点最小尺寸，当山洼地区面积小于最小节点尺寸时，被划分至同一节点的可能性较大，当点云数据被划分至最小节点后便不再继续分割。由于数据3中地形属于陡坡地形，在节点内点云平面拟合种子点选取时，高程值最低点并非山洼最深处一点，以至于平面拟合时，错误地将林下植被点以及部分低矮植被点拟合至山坡平面中而进入地面点集，造成山洼地形被平滑，但当山洼面积较大时，这种情况较少，原因在于大面积山洼地面点云数据会被划分至不同节点中，分解为多个山坡地形，减少了平面拟合中出现错误的可能性。具体滤波结果错误率见表 3。

表  3  数据3滤波精度验证

Table  3.  Precision verification of filtering in data 3

点云数量 Number of point cloud	第Ⅰ类错误率 Error rate of typeⅠ	第Ⅱ类错误率 Error rate of typeⅡ	总错误率 Total error rate
4 372	5.90%	5.79%	5.83%

下载: 导出CSV 

| 显示表格

根据滤波结果精度验证，其中第Ⅱ类错误率相比前两个数据明显增加，表明滤波算法在地形复杂的区域内还需要进一步改进；第Ⅰ类错误率相比数据2较为接近，进一步肯定了山坡地形中，地面点被划分为地物点的概率更小；样本的滤波总错误率有所提升，一定程度上表明该滤波算法对于复杂区域的多变性还存在一定错误率。

4.   结论与讨论

本研究中样本数据虽属于同一地区，但每一块样本区域环境有所不同。以3组点云数据样本为例，分别在地形缓和、地形陡峭及地形起伏不定3种情况下，验证滤波算法的可行性。根据3次实验结果，平坦地形下(数据1)滤波算法效果最好，总错误率为4.57%；陡坡地形下(数据2)，滤波算法能很好地区分山坡上的地面点与植被点，一定程度上验证了滤波算法中所引入的坡度判断条件的实用性，总错误率为4.75%；地形崎岖时(数据3)，对于部分山包与大面积山洼地区滤波效果较好，在保留地形特征的情况下，较为精确地分离出地面点，但部分山包与小面积山洼处掺杂了一定数量的地物点，总错误率为5.83%。

本文虽未使用国际公共测试数据集进行研究，但依照这些数据集的标准将目标区域样地划分为平坦地形(数据1)、陡坡地形(数据2)、崎岖地形(数据3)，对于不同的地形环境条件进行算法的改进与验证，其结果同样具有可比较性。对于ISPRS的激光雷达参照数据集开展的滤波研究，Axelsson^[3]、Elmqvist等^[39]、Pfeifer等^[40]、Shao^[41]均得到陡坡、崎岖地形误差大于缓坡误差的结果。例如：Axelsson^[3]的算法对于缓坡、陡坡、崎岖地形分别得到2.72%、11.52%、13.09%的错误率，Pfeifer等^[40]得到3.71%、19.64%、12.60%的错误率，Shao^[41]得到3.6%、2.80%、5.27%的错误率。从错误率来看，本文提出的算法对于崎岖陡峭的山区滤波效果较为突出，优于Axelsson^[3]的三角网模型、Pfeifer等^[40]的分层内插模型，精度与国际上表现较优的Shao^[41]的移动爬坡法较为相似。与这些算法相比，本文提出的改进八叉树滤波算法可以充分利用数据结构特征和优势，从而快速高效地实现高精度滤波，节约时间成本和运行成本。

综上所述，滤波算法还需要进一步改进，但若改进八叉树分割最小节点尺寸，势必造成其他问题，如迭代次数增多造成运算缓慢、节点中点云数量达不到平面拟合条件造成运算量增加等。由于时间有限，本研究仅能综合各方面因素，选取较好的滤波参数，达到最优化还需要进一步研究。