-
新方法的技术路线如图1所示,主要包括两个步骤:
(1) 借助区域增长方法将原始点云分割成点云块集和散点集。
(2) 基于多尺度表面插值滤波将块和散点分类为地面点和非地面点。
-
经典的基于平滑性约束的点云分块算法以空间距离、法向量和平面拟合残差为判定条件进行区域生长。但在低矮地物(图2区域A)和茂密植被(图2区域B)区域,基于主成分分析(PCA)方法拟合出的法向量和平面拟合残差不准确,容易将平坦区域的地面点和低矮物体(图2区域A)、陡坡区域的地面点和茂密植被(图2区域B)划分为同一块。
针对这一问题,文中提出了一种改进的区域增长算法,其中点法线的估计步骤如下:
(1) 使用尺寸略大于最大对象尺寸的移动窗口,通过扩展的局部最小值法[7]从点云中选择局部最低点,并将选择出的最低点作为初始地面种子点。
(2) 对于点云中的每个点,
从初始地面种子点中搜索k个最近邻点。
对k个邻近种子点通过PCA计算其法线和平面拟合残差。
改进后的区域增长算法首先选择平面拟合残差最小的点作为分块种子点,然后将点云中与种子点的距离小于距离阈值、法线夹角小于夹角阈值、拟合残差与种子点拟合残差的差小于残差阈值的点划分到当前块Sei。如图3所示,低矮物体(图3中区域A)和斜坡上的植被(图3中区域B)的法线虽和地面点法线相似,但与地面点相比其残差较大,故能够准确地将它们与地面点区分开。之后,将平均残差超过阈值,且点数小于一定数目(如20)的块划分为散点Sci。经过分块之后原始点云被划分为一组点云块Se={Se1, Se2, ···, Sek-1, Sek}(k为点云块的个数)和一组独立的散点Sc={Sc1, Sc2, ···Scm}(m为散点的个数)。
-
对划分好的点云块和散点同时通过多尺度表面插值滤波进行分类。该滤波过程主要包括以下优势:(1)通过改进的多尺度表面插值滤波,解决传统多分辨率层次滤波[6]存在的空间位置误差问题;(2)基于同一种方法同时处理点云块和散点,避免引入多种方法导致的不可控性;(3)结合点云块选择地面点,尽可能保留平滑区域中的不连续地形[11],同时,通过散点对复杂地形进行单独处理,避免基于块的分类在复杂地形的不稳定性,提高算法对地形细节的判断能力。
传统的多分辨率层次滤波分三个层次进行点云滤波,DTM分辨率和残差阈值随层次递增。在每个层次结构中,将已选定的地面点通过TPS插值生成参考DTM,并根据待分类点与DTM间的距离来选择地面点。
然而,在通过插值DTM模拟地表进行滤波的过程中,如果某点不在网格中心,则会存在空间位置误差问题,且该误差大小取决于该点距离网格中心的水平距离和当前区域的地形坡度。如图4所示,空间位置误差(
$\Delta h $ )可以表达为$\Delta h = d \times \tan \alpha $ ,其中α表示地形坡度,d表示待分类点距离网格中心的水平距离。由此可见,地形坡度和水平距离越大,空间位置误差越大。为此,文中在多分辨率层次滤波三层TPS插值基础上,将九个邻近的DTM网格中心点拟合为一个斜面来避免空间位置误差问题。例如,当某个点位于DTM表面的(I, J)单元格中,则它的九个最邻近网格为(I−1, J−1), (I−1, J), (I−1, J+1), (I, J−1), (I, J), (I, J+1), (I+1, J−1), (I+1, J)和(I+1, J+1)。为了避免一些地面种子点中错误分类的非地面点对平面拟合的影响,文中通过加权最小二乘来拟合地面参考面。
在每层迭代过程中,设点Pi对应的拟合平面方程为z=ax+by+c,则Pi距离地面的高度差可估计为:
$${d_i} = {{\textit{z}}_i} - (a{x_i} + b{y_i} + c)$$ (1) 式中:(xi,yi,zi)为Pi的三维坐标。
则点云块Sei中点的高度差可记为Di={de1, de2, ···, deij}(j为点云块Sei中所包含点的个数),点云块集Se的高度差记为DSe={D1, D2, ···, Dk},散点集的高度差记为DSc={dc1, dc2, ···, dcm}。
对每个点而言,高度差越小是地面点的可能性越大。因此,根据每个点的高度差,可以从点云块集和散点集中选择地面块和地面散点。针对点云块集Se,笔者通过以下规则从中识别地面块集GSe:
$$GSe = \left\{ {S{e_i} \in Se|{\rm{Num}}({D_i} < {d_{\rm th}}) \geqslant 50 {\text{%}} \times {\rm{Num}}(S{e_i})} \right\}$$ (2) 式中:dth为随迭代层次变化的高差阈值;Num(Di<dth)为点云块Sei中高度差小于高差阈值dth的点个数;Num(Sei)表示点云块Sei中点的总数。
针对散点集Sc,根据每个散点的高度差单独确定属性,将高度差小于高差阈值的散点添加到地面散点集GSe,即:
$$GSc = \left\{ {S{c_i} \in Sc|d{c_i} < {d_{th}}} \right\}$$ (3) 当前层次下识别的地面点GP即为地面块集和地面散点集中全部的点,即
$GP = \left\{ {Pi|Pi \in GSe \cup GSc} \right\}$ ,并将GP添加到地面种子点集。在下一层次的迭代过程中,对未分类的块和散点进一步识别,直到达到最大迭代层次,最终的种子点集即为滤波后的地面点。 -
文中提出的点云滤波伪代码如算法1所示,其流程可概括如下:
(1) 通过扩展的局部极小值方法从原始点云中选择初始地面种子点。
(2) 通过改进后的区域增长算法将点云划分为点云块集和散点集。
(3) 利用地面种子点通过TPS方法插值出分辨率为h的粗糙DTM。
(4) 利用每个网格及其周围八个邻近DTM值,通过加权最小二乘拟合为地面参考面。
(5) 对点云中的每个未分类点,计算其到所属地面参考面的高度差。
(6) 对每个待分类块
根据公式(2)从中识别地面块。
对每个待分类点
根据公式(3)从中识别地面散点。
(7) 根据地面块集和地面散点更新地面种子点集。
(8) 重复步骤(3) 、(4) ,直到没有发现新的地面点。
(9) 改变DTM分辨率和高差阈值,即:h=h/2和dth=dth+0.1;重复步骤(3)~(8),直到达到最大迭代层次。
算法1:PMMF滤波
Inputs: {P}: Raw point cloud, ws: moving window size, θth: angle threshold, rth: residual threshold, h: DTM resolution, dth: elevation difference threshold, Ω(·): extended local minimum method, Φ(·): improved region growing segmentation
Output: Ground seeds {S}
Initialize: Ground seeds {S}=
$ \emptyset $ , Segments index {SE}=$ \emptyset $ , Scatter index (SC)=$ \emptyset $ , Elevation difference {D}=$ \emptyset $ , Number of newly detected ground points n=11 Select ground seeds: {S}←Ω({P}, ws)
2 Segment {P} into a set of segments and one set of scatter: ({SE},{SC})←Φ({P},{S},θth, rth)
3 for i=1 to 3 do
4 while n!=0 do
5 n=0
6 Interpolate DTM using TPS: {DTM}←TPS({S},h)
7 Fit the plane of each grid using weighted least squares: {(a,b,c)} ←WLS({DTM})
8 for j=1 to Num({P}) do
9 Find the grid to which Pj belongs
10 Calculate the elevation difference between Pj and its grid: dj=zj-(ajxj+bjyj+cj)
11 {D}← dj
12 end for
13 for k=1 to Num({SE}) do
14 if Num(D{SEk}<dth)>50%×Num(SEk) then
15 {S}
$\mathop \leftarrow \limits^{insert} $ P{SEk}16 {SEk}
$\mathop \to \limits^{remove} $ {SE}17 n=n+Num({SEk })
18 end if
19 end for
20 for m=1 to Num({SC}) do
21 if D{SCm}<dth then
22 {S}
$\mathop \leftarrow \limits^{insert} $ P{SCm}23 SCm
$\mathop \to \limits^{remove} $ {SC}24 n=n+1
25 end if
26 end for
27 end while
28 h=h/2; dth=dth+0.1
29 end for
30 Return {S}
-
文中使用两种机载雷达数据集来评估该方法的性能,分别是国际摄影测量和遥感学会(ISPRS)第三委员会/WG3提供的低密度基准数据集和一组高密度的实例数据集。
-
该方法需要设置的五个参数: 用于选择初始种子点的移动窗口大小(ws),用于区域增长的法向量之间夹角的角度阈值(θth)和残余阈值(rth),以及用于判断地面点和非地面点的DTM初始分辨率(h)和初始高差阈值(dth)。对这些参数进行反复实验直到滤波结果的总误差最小。优化后的参数值和滤波精度如表1所示。
表 1 15组数据的最优参数和滤波精度
Table 1. Optimized parameters and filtering accuracy of the proposed method on the 15 samples
Sample Optimized parameters Accuracy measures ws/m θth/rad rth/m h/m dth/m T.I T.II T.E. κ 11 16 0.1 0.3 3 0.7 9.02% 10.69% 9.73% 80.15% 12 16 0.1 0.25 2 0.3 2.01% 3.08% 2.53% 94.93% 21 16 0.1 0.1 2 0.2 0.33% 3.13% 0.95% 97.23% 22 28 0.1 0.1 2 0.4 1.47% 8.16% 3.55% 91.61% 23 20 0.1 0.2 3 0.35 3.05% 5.74% 4.32% 91.33% 24 25 0.1 0.1 2 0.35 3.24% 8.65% 4.73% 88.38% 31 23 0.1 0.2 4 0.1 0.44% 1.20% 0.79% 98.40% 41 26 0.1 0.2 3 0.4 2.32% 2.58% 2.45% 95.11% 42 30 0.1 0.05 3 0.35 0.80% 0.77% 0.78% 98.13% 51 15 0.1 0.1 4 0.3 1.11% 4.11% 1.73% 94.92% 52 19 0.1 0.1 6 0.41 1.11% 17.10% 2.79% 84.63% 53 6 0.25 0.4 3 0.33 1.20% 37.87% 2.68% 63.79% 54 16 0.1 0.25 2 0.4 1.31% 3.36% 2.41% 95.17% 61 6 0.3 0.2 2 0.7 0.77% 7.05% 0.98% 86.19% 71 21 0.1 0.05 3 0.6 0.73% 6.33% 1.36% 93.20% Average 1.93% 7.99% 2.79% 90.21% Median 1.20% 5.74% 2.45% 93.20% STD 2.15% 9.27% 2.28% 9.00% 由表1可见,除samp42外,其余样本的I类误差均小于II类误差;平均而言,I类误差大约是II类误差的1/4。这一结果可能是由以下两个原因造成的:(1)由于样本中地面点的数量远大于非地面点的数量,故少数分类错误的非地面点便会导致较大的II类误差;(2)由于文中方法采用了地面点加密策略,当高差阈值较大时,在增加地面点数量的同时也会添加更多的非地面点。
将文中方法与近五年(2015~2019)新提出的15种滤波方法[14]的总误差进行比较,结果如图5所示。文中方法在15个样本中有11个样本的滤波性能优于目前最先进的滤波方法,而在samp11、samp24、samp41和samp52这四个样本中,新方法精度与最佳滤波方法基本一致。文中方法平均总误差为2.79%,与经典滤波算法相比精度至少提高了20.5%,该结果证明了文中方法具有较好的滤波性能。
图 5 文中方法与其他15种滤波方法的总误差比较
Figure 5. Total errors of the proposed method against and the 15 filtering algorithms
文中方法在samp11、smp21、samp51和samp52的参考DEM、滤波后DEM和误差分布如图6所示。这四组数据分别对应四种不同的地形环境。其中samp11为陡坡地形且混有大量的植被和建筑物,samp21中有一窄桥,samp51的陡坡上包括有大量的低矮植被,samp52为不连续断裂地形。在这四种地形中,现有的各类滤波算法都很难取得满意的滤波结果。相较而言,新算法的滤波结果较好,滤波后的DEM与参考DEM有很高的相似度。该算法几乎完全滤除了samp21中的桥梁部分,可精确滤除samp51斜坡上全部的低矮植被,能完美保持samp52断裂地形上的断裂线。然而,新方法也包含一些分类错误的非地面点(如图6中的红框标记部分),但这些点很容易通过人工编辑的方式剔除。
-
ISPRS数据集的平均点间距为1~3.5 m,远低于现今设备所采集到的点云密度。因此,文中通过三组不同地形(s1,s2和s3)的高密度实例数据对该算法进一步验证,分别为覆盖低矮植被的崎岖山地(s1,图7(a))、覆盖不同高度植被的断裂地形(s2,图7(b))和包含植被和建筑物的斜坡地形(s3,图7(c))。各个区域面积约为23公顷(1公顷=104 m2),平均点间距约为0.5 m。每组参考数据均先通过滤波软件进行初始滤波,然后对滤波结果进行人工检查和编辑,以保证参考数据的质量。
除了文中方法外,还将多分辨率层次滤波(MHC)[6]、布料模拟滤波(CSF)[15]、渐进加密三角网滤波(PTD)[16]和渐进形态学滤波(PMF)[17]应用于实例数据滤波,并与文中方法结果比较。如表2所示,在三种地形中,文中方法的性能均优于其他方法。就平均总误差而言,文中方法较精度最低的CSF算法提高75%,较精度最高的MHC算法提高40%。
表 2 该方法与其他四种方法的滤波精度比较
Table 2. Filtering accuracy of the proposed method against the 4 filtering algorithms
Method CSF PTD PMF MHC Proposed Accuracy measure Total Kappa Total Kappa Total Kappa Total Kappa Total Kappa s1 10.38% 67.67% 7.44% 75.49% 7.43% 74.41% 7.44% 76.29% 3.91% 92.15% s2 14.21% 65.09% 7.30% 82.87% 8.03% 81.04% 5.06% 88.20% 2.95% 93.19% s3 15.95% 68.25% 5.40% 89.20% 3.85% 92.31% 4.19% 91.61% 3.19% 93.61% Average 13.51% 69.21% 6.71% 83.07% 6.44% 85.79% 5.56% 84.77% 3.55% 92.62% 五种滤波方法对s2数据滤波后生成的DEM如图8所示。其中MHC算法得到的DEM最平滑,但同时也损失了一些地形细节。PMF算法滤除了过多的地面点,导致DEM产生过大畸变。PTD和CSF算法II类误差过大,导致DEM表面过于粗糙。相比之下,文中算法滤波后生成的DEM与参考DEM相似度最高。
Interpolation-based filtering with segmentation for airborne LiDAR point clouds
-
摘要: 现有机载激光雷达(LiDAR)点云滤波算法在简单地形下取得了较好的滤波效果,但普遍对陡坡地形适应性较差。为提高在不同地形下的滤波性能,提出了基于分块的多尺度表面插值滤波算法。该算法首先通过改进的区域增长分块算法将原始点云分为点云块集和散点集,然后通过构建的多尺度表面插值算法同时对点云块和散点进行分类。利用国际摄影测量与遥感学会(ISPRS)提供的基准数据验证表明,该方法在15个样本中有11个样本滤波效果优于现有滤波方法,对各类地形均有较强适应性,且该方法平均总误差最小。对三种不同地形特征的高密度数据滤波实验,也验证了该方法的良好性能。Abstract: The classical airborne LiDAR filtering algorithms show good results on most landscapes, but they suffer from low-level adaptability on steep slopes. Thus, to improve the filtering performance under different environments, a surface interpolation-based filtering algorithm with segmentation was proposed. Firstly, the original point clouds were grouped into a set of segments and one set of scattered points by an improved region growing method. Then, the segments and the scattered points were classified simultaneously using a weighted least square algorithm. The benchmark dataset provided by International Society for Photogrammetry and Remote Sensing(ISPRS) was used to validate the performance of the proposed method. Results illustrate that the proposed method outperforms the state-of-the-art filtering methods on 11 out of 15 samples, showing its strong adaptability to different terrain environments. Moreover, the proposed method has the lowest average total error. Filtering three samples of high-density with different terrain features also demonstrates the promising performance of the proposed method.
-
Key words:
- airborne LiDAR point clouds /
- filtering /
- interpolation /
- multi-scale /
- segmentation
-
表 1 15组数据的最优参数和滤波精度
Table 1. Optimized parameters and filtering accuracy of the proposed method on the 15 samples
Sample Optimized parameters Accuracy measures ws/m θth/rad rth/m h/m dth/m T.I T.II T.E. κ 11 16 0.1 0.3 3 0.7 9.02% 10.69% 9.73% 80.15% 12 16 0.1 0.25 2 0.3 2.01% 3.08% 2.53% 94.93% 21 16 0.1 0.1 2 0.2 0.33% 3.13% 0.95% 97.23% 22 28 0.1 0.1 2 0.4 1.47% 8.16% 3.55% 91.61% 23 20 0.1 0.2 3 0.35 3.05% 5.74% 4.32% 91.33% 24 25 0.1 0.1 2 0.35 3.24% 8.65% 4.73% 88.38% 31 23 0.1 0.2 4 0.1 0.44% 1.20% 0.79% 98.40% 41 26 0.1 0.2 3 0.4 2.32% 2.58% 2.45% 95.11% 42 30 0.1 0.05 3 0.35 0.80% 0.77% 0.78% 98.13% 51 15 0.1 0.1 4 0.3 1.11% 4.11% 1.73% 94.92% 52 19 0.1 0.1 6 0.41 1.11% 17.10% 2.79% 84.63% 53 6 0.25 0.4 3 0.33 1.20% 37.87% 2.68% 63.79% 54 16 0.1 0.25 2 0.4 1.31% 3.36% 2.41% 95.17% 61 6 0.3 0.2 2 0.7 0.77% 7.05% 0.98% 86.19% 71 21 0.1 0.05 3 0.6 0.73% 6.33% 1.36% 93.20% Average 1.93% 7.99% 2.79% 90.21% Median 1.20% 5.74% 2.45% 93.20% STD 2.15% 9.27% 2.28% 9.00% 表 2 该方法与其他四种方法的滤波精度比较
Table 2. Filtering accuracy of the proposed method against the 4 filtering algorithms
Method CSF PTD PMF MHC Proposed Accuracy measure Total Kappa Total Kappa Total Kappa Total Kappa Total Kappa s1 10.38% 67.67% 7.44% 75.49% 7.43% 74.41% 7.44% 76.29% 3.91% 92.15% s2 14.21% 65.09% 7.30% 82.87% 8.03% 81.04% 5.06% 88.20% 2.95% 93.19% s3 15.95% 68.25% 5.40% 89.20% 3.85% 92.31% 4.19% 91.61% 3.19% 93.61% Average 13.51% 69.21% 6.71% 83.07% 6.44% 85.79% 5.56% 84.77% 3.55% 92.62% -
[1] Liu Zhiqing, Li Pengcheng, Chen Xiaowei, et al. Classification of airborne LiDAR point cloud data based on information vector machine [J]. Optics and Precision Engineering, 2016, 24(1): 210-219. (in Chinese) doi: 10.3788/OPE.20162401.0210 [2] Wang Yanan, Wang Tingfeng, Tian Yuzhen, et al. Improved local convexity algorithm of segmentation for 3D point cloud [J]. Chinese Optics, 2017, 10(3): 348-354. (in Chinese) doi: 10.3788/co.20171003.0348 [3] Liu Zhiqing, Li Pengcheng, Guo Haitao, et al. Airborne LiDAR point cloud data classification based on relevance vector machine [J]. Infrared and Laser Engineering, 2016, 45(S1): S130006. (in Chinese) doi: 10.3788/IRLA201645s1.130006 [4] Feng Fajie, Ding Yazhou, Li Junping, et al. Airborne LiDAR point cloud filtering using saliency division [J]. Infrared and Laser Engineering, 2020, 49(8): 20190439. (in Chinese) [5] Sun Meiling, Li Yongshu, Chen Qiang, et al. Iterative multi-scale filter based on morphological opening by reconstruction for LiDAR urban data [J]. Infrared and Laser Engineering, 2015, 44(1): 363-369. (in Chinese) doi: 10.3969/j.issn.1007-2276.2015.01.061 [6] Chen C, Li Y, Li W, et al. A multiresolution hierarchical classification algorithm for filtering airborne LiDAR data [J]. ISPRS Journal of Photogrammetry and Remote Sensing, 2013, 82: 1-9. doi: 10.1016/j.isprsjprs.2013.05.001 [7] Zhao Chuan, Zhang Baoming, Yu Donghang, et al. Airborne LiDAR point cloud classification using transfer learning [J]. Optics and Precision Engineering, 2019, 27(7): 1601-1612. (in Chinese) doi: 10.3788/OPE.20192707.1601 [8] Zhang J, Lin X. Filtering airborne LiDAR data by embedding smoothness-constrained segmentation in progressive TIN densification [J]. ISPRS Journal of Photogrammetry and Remote Sensing, 2013, 81: 44-59. doi: 10.1016/j.isprsjprs.2013.04.001 [9] Zhao X, Guo Q, Su Y, et al. Improved progressive TIN densification filtering algorithm for airborne LiDAR data in forested areas [J]. ISPRS Journal of Photogrammetry and Remote Sensing, 2016, 117: 79-91. doi: 10.1016/j.isprsjprs.2016.03.016 [10] Chen C, Li Y, Yan C, et al. An improved multi-resolution hierarchical classification method based on robust segmentation for filtering ALS point clouds [J]. International Journal of Remote Sensing, 2016, 37(4): 950-968. doi: 10.1080/01431161.2016.1142687 [11] Yang B, Huang R, Dong Z, et al. Two-step adaptive extraction method for ground points and breaklines from lidar point clouds [J]. ISPRS Journal of Photogrammetry and Remote Sensing, 2016, 119: 373-389. doi: 10.1016/j.isprsjprs.2016.07.002 [12] Bigdeli B, Amini Amirkolaee H, Pahlavani P. DTM extraction under forest canopy using LiDAR data and a modified invasive weed optimization algorithm [J]. Remote Sensing of Environment, 2018, 216: 289-300. doi: 10.1016/j.rse.2018.06.045 [13] Meng X, Lin Y, Yan L, et al. Airborne LiDAR point cloud filtering by a multilevel adaptive filter based on morphological reconstruction and thin plate spline interpolation [J]. Electronics, 2019, 8(10): 1153. doi: 10.3390/electronics8101153 [14] Chen C, Wang M, Chang B, et al. Multi-level interpolation-based filter for airborne LiDAR point clouds in forested areas [J]. IEEE Access, 2020, 8: 41000-41012. doi: 10.1109/ACCESS.2020.2976848 [15] Zhang W, Qi J, Wan P, et al. An easy-to-use airborne LiDAR data filtering method based on cloth simulation [J]. Remote Sensing, 2016, 8(6): 501. doi: 10.3390/rs8060501 [16] Axelsson P. DEM generation from laser scanner data using adaptive TIN models [J]. International Archives of Photogrammetry and Remote Sensing, 2000, 33(4): 110-117. [17] Pingel T J, Clarke K C, Mcbride W A. An improved simple morphological filter for the terrain classification of airborne LIDAR data [J]. ISPRS Journal of Photogrammetry and Remote Sensing, 2013, 77: 21-30. doi: 10.1016/j.isprsjprs.2012.12.002