1 引 言

航空激光雷达(airborne LiDAR，或简称LiDAR)系统获取的原始数据也称为点云(point cloud)，其数据结构是离散的，并且以不规则方式分布在地形表面上。对离散的航空激光雷达点云数据进行分类滤波处理是该系统最主要的组成部分，这种处理通常称为数据后处理(post processing)。正常情况下，这种分离处理可以用滤波来表达，简单地说，就是将分布在地形表面上的点(即地面点)与那些非地形表面上的点(即非地面点，如在车辆、植被以及房屋上的点)进行有效而准确的分离。

到目前为止，绝大部分LiDAR点云数据滤波方法实际上是一种直接分离方法，例如“滑动最小值”以及“凸-凹罩的方法”^[1]，基于三角网TIN的方法等^[2]。前一种方法需要通过内插得到一个规则分布的高程栅格网，然后在这个规则分布的高程栅格网中寻找最低点来实现离散点云数据的分类。后一类方法是通过点云构造TIN，然后在每一个三角网内寻找满足一定条件的最低点作为地面点。这类方法中要引入所谓的“坎线长度”或“高度阀值”作为控制参数，如果在寻找的范围内，最低高程值小于这个阀值，则此点被接受为地面点，否则被假设为非地面点。

另外一种直接滤波方法是“最小二乘预测推估法”^[3]，这种方法从所有的点云数据出发，通过一个“平差”过程来达到滤波目的。这个平差过程是逐步迭代进行的，在迭代开始时，所有的原始点云数据假设具有同样的“权”。随着迭代过程的重复，那些较低的点被逐渐赋予较大的权系数，相反，较高的点其权值逐渐减少，每一次迭代后，每一个离散点的权系数被重新计算。与此方法类似的还有移动曲面及类似的方法^[4]，这种方法首先设置一个曲面函数，然后通过迭代的方法完成点云数据的分类滤波。

基于数学形态学理论的LiDAR点云数据滤波同样是一种直接分类方法。其中文献^[5]最早将数学形态学原理用于分类LiDAR点云数据，提出了一维数学形态学算法用于点云数据的处理。其他一些学者对传统数学形态学算法进行了改进，提出了不同的算法用于LiDAR点云数据的分类^{[6, 7, 8, 9]}。

这类直接处理方法是基于“地形表面(DSM)上的非地面点高于地面(DEM)点”和“激光扫描得到的最低点是足够准确的地面点”的前提而衍生的算法^{[10, 11, 12, 13]}，根据不同的判别准则，或直接比较在一定邻域内的点的高程值，或通过直接平差的方法完成对点云数据的分类。还有融合影像与LiDAR点云进行分类的算法等^{[14, 15]}。

“最小二乘预测推估法”用于LiDAR点云数据分类时，是假设所有的原始LiDAR点云数据具有同样的“权”，这种假设在地形变化不剧烈时结果比较可靠，但是，如果地形变化很大，迭代计算不理想。因此，笔者将数学形态学算法与地形逼近思想结合在一起，提出了一种“分步”处理的思想用于LiDAR点云数据的分类。其原理是将整个处理过程分成两步：第1步通过一个相对简单的处理过程将原始的离散点云数据进行“粗分离”，经过粗分离后，原始离散点云数据被分成两组，一组是近似的地面点，即“地面点假设”，另外一组是近似的非地面点，即“非地面点假设”^[16]；在粗分离的基础上，第2步引入顾及因果关系的二维自回归car(causal auto-regressive)模型对地形表面进行逼近和模型化处理，它包括模型化和假设检验。第2步的处理也称为“细分离”或“精处理”。“粗分离”处理过程中采用计算速度较快的数学形态学算法，“细分离”则采用在统计学理论中成熟的car模型。相比单独的数学形态学算法和“最小二乘预测推估法”，这种思想更合理，结果更可靠。

本文主要采用了数学形态学“开”算子。“开”算子是两个基础算子“腐蚀”(erosion)和“膨胀”(dilation)的组合，与“开”算子相对应的是“关”算子。

假设f(x)和g(z)是任意两个集，在一维情况下，集g(z)称为结构元素或核；在二维图像或目标时，集g(z)可以称为结构元素或结构窗口(structuring element or window)。在数字图像处理中，集g(z)表示操作算子点集，f(x)表示图像(或目标)点集，结构元素可以定义不同的形式，例如菱形结构元素、三角形结构元素、8邻居结构元素或者斜形结构元素等。“开”和“关”算子是这样定义的^[17]：

(1) “开”算子(opening,用“”表示)。一幅图像或目标f(x)通过结构窗口g(z)被打开，意味着图像或目标f首先通过结构窗口g被腐蚀，然后在腐蚀后的图像或目标的基础上再进行膨胀操作，在数学上这个过程可以表示为

点击浏览公式

(2) “关”算子(closing,用“·”表示)。一幅图像或目标f(x)通过结构窗口g(z)被关闭，意味着图像或目标f首先通过结构窗口g被膨胀，然后在膨胀后的图像或目标的基础上再进行腐蚀操作，在数学上这个过程可以表示为

点击浏览公式

线性回归估计模型表达如下：

参数估计的过程是在一般的线性模型

点击浏览公式

中寻找参数h_k的过程。式中，x_t－i，i=…,-2,-1,0,1,2,…表示随机变量；e_t~N(0,σ_e²)表示白噪声。对此模型进行扩展，可以引入一体化的自回归处理过程用于一维数据的处理^{[18, 19]}。

笔者对这种一维处理过程进行了扩展，引入更有效的二维估计方法，其中之一是所谓的“顾及因果关系的自回归模型”(causal auto-regressive process，简写为car模型)。一个阶为max(p,q)的，顾及因果关系的二维自回归模型car(p,q)可以写成^[20]

点击浏览公式

式中，a_kl表示系数；x(m,n)表示随机变量；e(m,n)表示白噪声。

“最小二乘预测推估法”用于LiDAR点云数据分类时，假设所有的原始LiDAR点云数据具有同样的“权”。这种假设不是很合理，特别是当地形变化剧烈或存在地形断裂线时，因为没有很好的初值，这种方法需要较多的迭代，计算比较费时。笔者采用的“分步”处理过程简述如下。

首先采用计算相对简单的数学形态学“开”算子对原始点云数据进行“粗分类”，并且对传统“开”算子进行了扩展和改进，即在“开”算子的基础上增加一个“带宽”参数。因此，传统的数学形态学“开”算子则改进为

点击浏览公式

式中，bw_th为“带宽”。由于本文采用的是分步处理方式，经验系数的大小仅影响“地面点假设”权的初值，此值的选择不决定点云是地面点还是非地面点，因此其值的选择可以比较宽松，文献^{[5, 21]}对此就行了详细描述。

“开”处理过程为：① 对每一个离散点首先进行“腐蚀”处理，然后在“腐蚀”处理后得到的结果的基础上进行“膨胀”处理；② 计算“膨胀”操作时当前点的高程与“腐蚀”处理得到的最低高程之差Δh；③ 根据差值Δh判断当前点的类型，如果Δh大于bw_th则为“非地面点假设”，反之，如果高程之差Δh小于bw_th，则判断此点为“地面点假设”；④ 根据高程之差值Δh的大小，将“地面点假设”点赋予不同的权，其权值为

点击浏览公式

式中，Δh为按照形态学“开”算子计算得到的高程之差值，其值总是为正。

按照上述思想对每一个点云点进行处理后，得到了两类假设，“地面点假设”是地面点的可能性高，其权值按照公式(6)计算得到；“非地面点假设”是地面点的可能性相对较低，其权值可以统一给定为0，也可以赋予相对小的权重。经过数学形态学“粗”处理后，每一个点都被赋予了不同的权。

在粗分离的基础上，第2步是采用顾及因果关系的car模型对地形表面进行逼近和模型化处理。在上述“粗分离”后得到的结果基础上，将所有的点按照式(4)进行地面拟合和假设检验。这一步的处理可以按照传统的最小二乘法进行，处理过程通过迭代的方式完成，根据迭代计算得到每一个点上的残差值。如果残差值满足假设检验的条件，则判断该点为地面点，如果超过给定的限差，则判断该点为非地面点。而每一次迭代时，每一个点的权值通过上一次迭代后得到的参数值重新计算，这种方法也称为IRLS方法(iteratively reweighted least square method)。在迭代计算时，权函数可以采用式(7)的形式^[3]，也可以选用文献^[22]中介绍的权函数。

点击浏览公式

式中，gw_i表示当前点的权值；v_i表示当前点前一次解算后的残差；参数a和b则主要用于决定上述权函数的陡峭程度，例如可以分别选择1和4。算法的整体流程如图 1所示。

图 1 分步处理流程图 Fig. 1 Flowchart of the separated filtering method

图选项

根据上述算法原理，本文选择了两块LiDAR试验数据进行了滤波试验。对第1个试验区域的试验结果采用三维显示和断面显示的方式进行演示，第2个试验区域采用三维显示的方式，试验区域基本情况见表 1。图 2则演示了两个试验区域的LiDAR点云高程值图，高程值图的生成方式按照密度分割的原理，图 2中灰度值越低，对应的高程越低，反之越高。

图 2 试验区域 Fig. 2 Test areas

图选项

在两个试验区域内，分别采用car(5,5)、car(7,7)和car(9,9) 3种模型测试不同参数条件下滤波分类的结果，以内部精度、处理时间和分类后地面点/非地面点的数量来考察不同参数的影响，见表 2。从表 2可以看出，在本试验区域情况下，采用不同参数的3种处理模型的内部精度和分类后地面点数量差别不大，而处理时间差别较大，采用car(5,5)模型可以满足要求，因此，下列试验结果都采用car(5,5)模型。

按照本文提出的计算流程，首先对原始点云数据按照数学形态学“开”算子进行处理，处理时改变通常数学形态学处理要首先进行格网化处理的方法，不进行点云数据的内插，直接对离散数据进行处理；“开”算子处理后，所有的点被分成两类点云数据假设，即“地面点假设”和“非地面点假设”，并且具有不同的权；将两类数据假设按照car模型进行处理，并进行假设检验；根据模型计算的结果，得到最终的“地面点”和“非地面点”。两个试验区域试验结果见图 3、图 4和图 5。

图 3 试验区域Ⅰ三维显示图 Fig. 3 3D display view of test area Ⅰ

图选项

图 3(a)表示试验区域Ⅰ原始LiDAR点云数据三维视图，图 3(b)演示了经过car(5,5)模型处理后得到的地面点三维视图。图 4(a)和图 4(b)分别表示了原始数据的一幅高程值断面图和经过处理后的相应的高程值断面图。

图 4 试验区域Ⅰ高程值断面图 Fig. 4 Height profile of test area Ⅰ

图选项

图 5 试验区域Ⅱ三维视图 Fig. 5 3D display view of teat area Ⅱ

图选项

试验区域Ⅱ的分类滤波结果见图 5，试验区域只给出了三维显示的结果。由两个试验地区的图示可以看出，这种滤波方法用于LiDAR点云数据处理取得了令人满意的结果，试验区域内的森林、建筑物、特别是低矮的植被被准确滤除，与此同时，地形所包含的坎线特征则保持得很好。

为了对上述试验结果进行验证，按照文献[22—23]提出的定量分析方法进行分析对比，即计算Ⅰ类误差：×100%；Ⅱ类误差：×100%；总误差Ⅲ：×100%。其中，a为滤波得到的地面点总数；b为非地面点总数；c为被错误判断为非地面点的地面点个数；d为被错误判断为地面点的非地面点个数。

按上述误差计算方法，笔者采用人机交互的方法统计了两个试验区域的参数a、b、c和d，并计算了Ⅰ类误差、Ⅱ类误差和总误差Ⅲ。根据计算结果，两个试验区域Ⅰ类误差小于Ⅱ类误差，准确率基本上高于95%，这一结论与文献相似。

定量分析后进行了滤波对比试验，采取与基于Terrasolid软件滤波结果进行对照的方法，在上述两个区域按照Terrascan软件进行自动滤波，然后进行人工交互编辑，检查建模后的DSM模型和地面点的数量，采用本文的滤波方法分类为地面点的数量略少于采用Terrascan滤波后的数量，二者之间相差不超过1%。

除了与Terrascan软件自动滤波结果进行了对比试验外，与传统“最小二乘预测推估法”(least squares prediction，LSP)滤波结果也进行了相应对比，试验结果见表 3。采用分步处理算法分类为地面点的数量略少于采用“最小二乘预测推估法”滤波后的地面点数量，表明算法更稳健。

建立在“直接处理”基础上的“最小二乘拟合预测法”针对的是整个数据集，分块大小固定，在处理城市区域数据时往往会出现误处理，在处理地面存在断裂线或突变地形时无法完全保持地形特征^[22]。并且基于“最小二乘预测推估法”进行LiDAR点云数据滤波时，假设所有的原始LiDAR点云数据具有同样的“权”，然后将原始点云数据引入到平差模型中进行处理，根据平差后得到的残差的大小判断每一个点属于地面点还是非地面点。

针对上述不足，本文提出的算法作了几方面的改进:

(1) 对原始的LiDAR点云数据不进行格网内插，而使用数学形态学滤波模型对点云数据进行“粗处理”，直接通过数学形态学“开”算子和自回归处理，可以不改变原始数据的结构，避免了内插误差，同时提高了数据处理的效率。

(2) 提出了“地面点假设”和“非地面点假设”的思想，在“开”算子的基础上增加“带宽”参数作为判断“地面点假设”的方法，对传统的数学形态学算法进行了扩展。

(3) 采用了“分步”处理的思想，在“粗处理”的基础上，借助于car(p,q)模型对地形进行模型化处理，对数据进行“假设检验”。因而，相对于“直接处理方法”而言，笔者提出的算法具有更好的稳健性，滤波分类的结果也更可靠。