基于点的多尺度形态学重建滤波方法

常兵涛， 陈传法， 郭娇娇， 武慧明， 贝祎轩， 李琳叶

1. 山东科技大学 测绘与空间信息学院, 青岛 266590;

2. 山东省基础地理信息与数字化技术重点实验室, 青岛 266590

1 引 言

近年来，机载激光探测与测距技术（LiDAR）发展迅速，逐渐成为地理信息科学（GIS）数据获取的有效工具（Chen 等，2013），并广泛应用于数字地面模型（DTM）生产（胡翰 等，2019）、滑坡监测（马洪超 等，2008）、林业调查与管理等（刘浩 等，2018）。由于原始LiDAR点云包括地面点和非地面点，因此，在进行上述应用之前需要对原始点云滤波以区分地面点和地物点。其中，滤波结果直接影响后续应用精度。目前，国内外研究人员提出了大量机载LiDAR 点云滤波算法，根据其原理可分为5 大类：基于插值的滤波（Axelsson，2000；Chen 等，2013 和2020；Evans和Hudak，2007；Kraus 和Pfeifer，1998；苏伟 等，2009；詹总谦 等，2020）、基于机器学习的滤波（Gevaert 等，2018；Jahromi 等，2011；Lu 等，2009；Luo 等，2017）、基于分块的滤波（Lin 和Zhang，2014；Tóvári 和Pfeifer，2005；Yang 等，2016）、基于坡度的滤波（Liu，2008；Meng 等，2010；Shao 和Chen，2008；Sithole 和Vosselman，2001；Susaki，2012；Vosselman，2000；Wang 和Tseng，2010） 以及基于形态学的滤波（Bigdeli等，2018；Chen 等，2007；Hui 等，2016；Kilian等，1996；Li 等，2014 和2017；Meng 等，2019；Mongus 等，2014；Pingel 等，2013；Zhang 等，2003）。

相比其他滤波算法，形态学滤波以精度高、计算速度快等优势备受关注（Pingel 等，2013）。其中，该方法滤波过程中结构元素的选择显著影响运算精度（Kilian 等，1996）。例如，Kilian 等（1996）通过特定结构元素的简单开运算检测非地面点，但结构元素尺寸严重影响最终滤波精度；为 此，Zhang 等（2003）、Chen 等（2007） 及Pingel 等（2013）通过一组尺寸线性增大的结构元素进行形态学开运算，但单独的形态学开运算难以保持地形细节信息；Mongus 等（2014）、Li 等（2014）及Hui 等（2016）通过形态学顶帽变换增强以提升地形细节保持水平，但这些方法依旧依赖窗口的设置；Li 等（2017）、Bigdeli 等（2018）及Meng 等（2019）通过形态学测地变换对点云滤波，具有收敛速度快、不需人为选择窗口尺寸等优点。然而，形态学测地变换滤波仍然存在一些问题，例如，该方法运算过程中需要对点云栅格化，由此造成点云细节损失（Bigdeli 等，2018；Meng等，2019），而且在地面点选择时仅考虑了点与点间的相对高差，容易导致地形断裂点错分（Li等，2017）。

针对上述问题，本文提出了一种直接基于点的多尺度形态学重建点云滤波方法PMMF（Pointbased Multi-scale Morphological reconstruction Filter）。该方法主要创新点包括：（1）以原始点云的k邻域作为结构元素，避免点云栅格化造成的信息损失；（2）在基于点的测地膨胀过程中，构建了自适应膨胀和高程缓冲区，降低陡坡区地面点错分和漏分问题；（3）发展了自适应坡度方法，用以进一步剔除潜在地面点中的非地面点，降低斜坡区Ⅱ类误差（地物点误分为地面点的比例）影响，提升滤波结果精度。

2 原理与步骤

传统形态学重建中的测地线膨胀算法需要输入两幅图像，即标记图像和掩膜图像（Soille，2003）。测地膨胀即对标记图像重复膨胀，并通过掩模图像限制膨胀的传播，具体表达式为：

式中，G为标记图像，F为掩模图像，B为结构元素，n为迭代次数，⊕为膨胀运算，Δ为取最小值运算。

为避免空间插值导致的精度损失，PMMF将原始点云作为操作对象，即以标记点云和掩膜点云分别替代标记图像G和掩模图像F，以k邻域点集作为结构元素B。PMMF分3层滤波，分别是Ⅰ层、Ⅱ层和Ⅲ层（图1），从底层到高层网格分辨率逐层提高。在每层中，首先对原始点云覆盖一定尺寸网格，并获取网格内最低点用以赋值标记点云和掩膜点云；然后在掩膜点云约束下，借助结构元素和高程缓冲区迭代膨胀标记点云，直至标记点高程收敛；最后利用膨胀后标记点云获取潜在地面点，并通过自适应坡度方法剔除潜在地面点中的地物点。算法流程如图1所示。

图1 PMMF滤波流程图Fig. 1 Flowchart of the proposed method

2.1 基于点的数据组织

首先定义掩模数据为原始点云，标记数据为一组标记点云。其中，标记点云平面坐标与掩模点云相同，但高程z不同。标记点云高程需满足以下条件：（1）不能过大，即地物点在标记点云中的高程应低于其真实高程，否则会导致仅能检测出建筑物脊线和植被冠层；（2）不能过小，即标记点云中地面点高程应尽可能接近其真实高程，以保证不会滤除凸起地形。为此，首先将原始点云用分辨率为r的网格覆盖，然后搜索网格内最低点，最后以最低点高程作为网格内标记点的高程（图2）。

图2 标记点云与掩模点云Fig. 2 Marker point cloud and mask point cloud

传统测地线膨胀需要定义结构元素（通常为n×n邻域网格），并将其在栅格图像上顺次移动实现膨胀运算。鉴于原始点云非结构化特性，传统结构元素无法对原始点云形态学运算（Jesús 等，2020）。因此，本文以k个二维邻近点（k邻域）作为PMMF 算法的结构元素，并按点云在内存中存储顺序依次对其测地膨胀。相比传统的n×n网格结构元素，k邻域结构元素具有各向同性以及对运算结果无偏差等优势。除此之外，由于点云的置换不变性，k邻域结构元素的移动次序不会对膨胀结果产生影响。

2.2 自适应测地膨胀

传统的形态学膨胀以结构元素内最高点作为当前点膨胀高程，但陡坡区最高点通常高于当前点高程，以该高程进行膨胀容易导致此区域低矮灌木等地物点错分为地面点。为此，PMMF构建了坡度自适应形态学膨胀方法，即以结构元素内最高点减去其空间位置误差∆h作为当前点膨胀高程（Zdilated=Zmax‒∆h）。其中，空间位置误差表示最高点和待膨胀点因平面位置不重合和地形起伏引起的高程差（图3），即∆h=d× tanα，式中d表示最高点和膨胀点水平距离，α表示待膨胀点处地形坡度。

图3 标记点云膨胀Fig. 3 Marker point dilation

测地膨胀过程中，受点云采集设备误差和微地形凸起影响，仅以掩膜点云作为约束容易导致陡坡区地面点错分为地物点。为此，PMMF构建了附有高程缓冲区掩膜点云（图4），即如果膨胀后的标记点高程Zmarker高于掩膜点云高程Zmask减去缓冲区宽度eth，则令Zmarker=Zmask。重复此迭代过程，直到标记点云高程不再改变。最后，选择标记点云与掩膜点云高程相同的点作为该层的潜在地面点。

图4 有和无高程缓冲区滤波结果对比Fig. 4 Accuracy comparison between dilations without and with a buffer

2.3 潜在地面点误差剔除

由于受地形场景复杂性影响，自适应测地膨胀滤波后的地面点中不可避免的混有少量地物点。因此，为保证结果滤波精度，PMMF 算法构建了一种自适应坡度滤波方法进一步剔除潜在地面点中的地物点。具体而言，首先平面拟合待求点Pi的m邻域局部最低点Lj（j=1，2，…，m）；然后借助式（2）计算Pi与Lj（j=1，2，…，m）连线与拟合平面的夹角θij（j=1，2，…，m）（图5），并取它们的平均值θi（式3）；最后将θi与坡度阈值δth比较，如果θi>δth，将该点从潜在地面点中剔除。

图5 坡度滤波示意图Fig. 5 Slope-based filtering

在地形断裂等变化剧烈区域（图6），仅以固定坡度阈值作为判断标准，会导致错误剔除地形特征点。为此，本文在坡度阈值中引入地形复杂度因子，以提高算法的地形细节保持能力。坡度阈值δth计算如下：

图6 自适应坡度阈值原理图Fig. 6 Adaptive slope-based threshold

2.4 算法步骤

综上所述，基于点的形态学重建滤波算法整体流程如下：

（1）将所有原始点云初始化为待选地面点PCandidate，初始化外迭代次数iter=1；设置初始网格分辨率r、筛选地面点的固定高差阈值es、坡度阈值δs、尺度缩放因子f。

（2）搜索PCandidate网格最低点，并存储在最低点集Pl中。

（3）根据网格最低点赋值标记点云PMarker，利用PCandidate赋值掩模点云PMask。

（4）搜索每个点Pi在PMarker中的k个二维邻近点（如k=12）并根据邻近点计算局部坡度αi。

（5）测地膨胀计算：根据αi计算Pi的膨胀高程Zdilated，并比较Zdilated与PMask中对应点高程Zmask，若Zdilated

（6）获取潜在地面点：将标记点云PMarker和掩膜点云PMask高程相同的点标记为潜在地面点PPotential。

（7）剔除潜在地面点残差：搜索PPotential中每个点在最低点集Pl中的m个邻近点Ni（如m=6），并拟合地面参考面；计算每个点对应的θi（式（3））和δth（公式4）。如果θi>δth，将该点从潜在地面点集PPotential删除，用PPotential更新PCandidate。

（8） 若iter≤3，则r=r/2，eiter=es- 0.1 ×iter，δiter=δiter- 0.02 ×iter，iter=iter+1，重复步骤2—7；否则，结束计算，并输出PCandidate。

3 实验结果与分析

本文使用两组机载雷达数据集来评估该方法的性能：第一组是国际摄影测量和遥感学会（ISPRS） 第三委员会/WG3 提供的低密度机载LiDAR 点云基准数据集；第二组是高密度机载LiDAR点云实例数据集。

3.1 ISPRS基准数据

该数据集包括多种不同复杂地形：陡坡、断裂地形以及植被和建筑物的混合区等。表1描述了15 组数据对应的地形特征、地面点（BE）和地物点（OBJ）的数量、点密度以及高程标准差（STD）。

表1 15组基准数据特征Table 1 Characteristics of the 15 benchmark samples

新算法须设置4 个参数：网格分辨率r，固定缓冲区宽度es，固定坡度阈值δs以及尺度缩放因子f；对15 组数据的参数最优值和对应滤波精度如表2所示。

表2 15组数据的最优参数和滤波精度Table 2 Optimized parameters and filtering accuracy of the proposed method on the 15 samples

结果分析表明，15 组数据中有12 组数据的Kappa 系数高于90%，平均Kappa 为91.08%，且除sample11 以外，其余数据的总误差均低于4%。平均而言，Ⅰ类误差（地面点误分为地物点的比例）大约是Ⅱ类误差（地物点误分为地面点的比例）的30%，这主要归因为：（1）PMMF构建的缓冲区和自适应坡度阈值使得滤波算法具有很好的地面点保持能力；（2）样本中的地面点数量远大于地物点数量，少数分类错误的非地面点便会导致较大的Ⅱ类误差。整体而言，Ⅰ类误差表现稳定，表明PMMF 对各类地形均有较好的地面点保持能力；Ⅱ类误差表现略差，但残留的非地面点易通过后期人工编辑去除。

将PMMF 与近5 年（2016 年—2020 年）新提出的15 种滤波方法（Chen 等，2021）总误差比较表明（图7），PMMF 在15 个样本中有8 个样本的滤波性能最优，而在其余7 个样本中，PMMF 精度与最佳滤波方法基本一致。比较而言，PMMF平均总误差为2.71%，精度比其他方法至少提高了12.6%。

图7 该方法与其他15种滤波方法的总误差比较Fig. 7 Total errors of the proposed method against the 15 filtering algorithms

samp12、smp21、samp51 和samp53 的 参 考DEM、滤波后DEM和误差分布如图8所示。这4组数据分别对应4 种不同的地形。其中samp12 为典型城区地形，存在大量建筑物和低矮植被，且其东北角坡度较大，地形复杂；samp21 中有一窄桥和高矮不等的大量植被；samp51 的陡坡上存在大量低矮植被；samp53 为不连续断裂地形。结果表明，滤波后DEM 与参考DEM 基本相同；对各类复杂地形能很好保留地形细节（samp12），可完全滤除桥梁信息（samp21），准确滤除斜坡上绝大部分的低矮植被（samp51），较好保持了地断裂线和陡坡地形细节（samp53）。综上所述，本文算法在各种场景中都有较好的地形特征保持能力。然而，PMMF 也残留了少量错分的地物点。例如，如图8 samp12和samp53的滤波后DEM 中椭圆部分，这可能是由该区域较高的地形复杂度导致坡度阈值较大引起。

图8 新算法在samp12、samp21、samp51和samp53结果Fig. 8 Results of the proposed method on samp12， samp21，samp51 and samp53

3.2 实例数据

受当时硬件设备限制，ISPRS 发布的数据点密度较低，与当前生产工作所使用的点云密度相差较大。为进一步验证新算法对高密度点云处理能力，本文对4组不同的高密度实例数据处理（图9），并将计算结果与经典滤波方法比较，包括：渐进形态学滤波（SMRF）（Pingel 等，2013）、布料模拟滤波（CSF）（Zhang 等，2016）、渐进加密三角网滤波（PTD）（Axelsson，2000）和多分辨率层次滤波（MHF）（Chen 等，2013）。所选取的4 组实例数据的面积均为500 m ×500 m，包括两组城区地形（数据1、数据2）和两组森林地形（数据3、数据4）。其中，数据1 为混合有建筑物和植被的沟壑地形，数据2中存在不规则建筑物和覆盖茂密植被的断裂地形，数据3为覆盖茂密植被的斜坡地形，数据4为覆盖不同密度植被的陡坡地形。数据采集设备为Optech Orion H300‒325 系统，其扫描角度为15°，脉冲频率为250 kHz。参考数据通过标准工业流程制作：先将原始点云通过商业软件TerraSolid 自动化处理，然后人工检核并编辑分类错误的点。4组数据参考DEM如图10所示。

图9 4组数据对应区域的正射影像图Fig. 9 Digital orthophoto maps of four private data

图10 4组数据的参考DEMsFig. 10 Reference DEMS of the four private data

由图11 可见，PMMF 对4 组数据滤波后的DEM 与参考DEM（图10）非常接近，表面平滑，只有数据4（图11（d））中有少数凸起地形细节丢失。新算法与4 种经典滤波算法（SMRF、CSF、PTD、MHF）比较表明（表3），各个方法对城市地形滤波精度优于林区地形；林区地形Ⅰ类误差大于Ⅱ类误差，这可能是由于林区地形植被密度较高，地面点稀少，少量分类错误的地面点便会导致较大Ⅰ类误差。对平均总误差而言，新算法为3.24%较SMRF、CSF、PTD 和MHF 分别减小了约12.0%、59.1%、70.1%和53.2%，kappa 系数为92.23%分别提高了约1.1%、13.5%、22.8%和9.7%。

表3 新算法与其他4种经典算法对私有数据精度对比Table 3 Accuracy comparison between the proposed method and the classical filters on the private dataset/%

图11 新算法对4组数据滤波后的DEMsFig.11 DEMs filtered by proposed method

PMMF 及4 种经典滤波算法对数据4 滤波后DEM如图12所示。结果表明，SMRF、PTD、MHC斜坡表面较为粗糙，同时PTD 水平表面也非常粗糙。CSF 斜坡上I 类误差过大，滤除了过多的地形细节。整体而言，PMMF 滤波后DEM 精度最高，只在陡坡上存在较少的异常凸起和滤除了部分地形细节。

图12 参考DEM及5种方法对数据4滤波后DEMFig. 12 Reference DEM and filtered DEMs of data 4 by five filters

4 讨 论

PMMF 较高滤波精度主要得益于以下3 点创新：（1）基于点的滤波方法；（2）在滤波过程中构建了自适应膨胀和高程缓冲区；（3）通过自适应坡度剔除潜在地面点中的地物点。为进一步验证这3 个创新点效果，本文以ISPRS 基准数据为例，比较了新算法分别在加入和不加入相应创新点时精度（表4）。

表4 在ISPRS数据上加入创新点与不加入的滤波精度比较Table 4 Accuracy of the proposed method with and without the contributions on the ISPRS data/%

首先，以samp23 为数据源，对比了基于点的PMMF算法和基于栅格的形态学滤波算法。结果表明，前者较后者Ⅰ类误差减小了24.8%，总误差减小了21.4%。如图13 所示，基于点的PMMF 算法在复杂地形边缘保留了更多的地面点，这是因为基于点的PMMF 算法避免了点云栅格化过程中的细节信息损失。

图13 基于点和基于栅格的形态学重建滤波误差分布Fig. 13 Errors of PMMF and grid-based morphological reconstruction filtering on samp23

其次，以samp52 数据，验证了在形态学重建过程中构建自适应膨胀和缓冲区的影响。结果表明，带缓冲区的自适应膨胀比传统膨胀方法总误差减小了72.3%，其中，Ⅰ类误差减小了84.5%；Ⅱ类误差略有增大可能是由平坦区域高程扩展造成的。由图14 可见，无缓冲区的形态学重建滤波在地形凸起区域滤除了大量地面点，这是由于标记点云的最高点低于凸起地形高程，致使标记点云无法膨胀至凸起地形。而构建缓冲区后，标记点云可通过缓冲区逐层向上膨胀，直至贴合真实地表；同时，自适应膨胀在斜坡地形上滤除了更多的地物点。

图14 构建缓冲区的自适应膨胀和无缓冲区的传统膨胀滤波误差分布Fig. 14 Errors of the proposed method with adaptive dilation and elevation buffer and classical dilation without elevation buffer on samp52

最后，借助samp24 数据验证了自适应坡度剔除潜在地物点对滤波结果的影响。结果表明，加入自适应坡度方法比加入固定坡度剔除和不加入坡度剔除方法滤波总误差分别减小了23.8%和33.3%，其中Ⅱ误差分别减小了36.9%和42.2%。如图15 所示，新方法在斜坡上有效减少了Ⅱ类误差影响。

图15 自适应阈值的坡度滤波、固定阈值的坡度滤波和无坡度滤波的误差分布Fig. 15 Errors of the proposed method with adaptive threshold，with fixed threshold and without removing non-ground points based on slope filter on samp24

5 结 论

为了提升各种场景机载LiDAR 点云滤波精度，本文提出了一种基于点的多尺度形态学重建滤波算法。该算法全流程以原始点云为操作对象，首先通过构建带缓冲区测地膨胀方法筛选潜在地面点，然后通过改进的自适应坡度方法剔除潜在地面点中的地物点。ISPRS基准数据实验表明，新方法在15个样本中有8个样本的滤波性能优于近5年发展的滤波算法，其平均总误差为2.71%，平均Kappa 为91.08%，均高于现有其它滤波方法。处理4组高密度机载LiDAR点云表明，新方法平均总误差为3.24%，比其它经典滤波方法（SMRF、CSF、PTD 和MHF）精度提升至少12.0%，平均Kappa 为92.23%，亦高于其它4 种滤波方法。创新性分析实验表明，PMMF较栅格化形态学滤波、不构建自适应膨胀和缓冲区及不加入自适应坡度滤波的总误差分别减小了21.4%、72.3%和33.3%，证明了创新点的有效性和必要性。然而，新方法需要人为设置4个参数，影响了算法的易用性。因此后续拟研究新方法最优参数的设置策略，以提高点云滤波的自动化和智能化水平。