-
近年来，机载激光探测与测距（LiDAR）技术凭借快速、高效、精确捕获地表三维信息的优势，广泛应用于数字地面模型（digital terrain models，DTM)生产[1]、数字城市构建[2]、林业资源调查与管理[3]等领域。原始LiDAR点云既包括目标区域内的地面信息，也包括建筑物、桥梁、植被等地物信息。因此，在进行上述应用之前需要通过点云滤波对地面和地物进行准确分离[4]。为此，研究人员提出了大量的点云滤波方法，根据工作原理可大致分为五类：基于坡度的滤波、基于形态学的滤波[5]、基于插值的滤波[6]、基于分块的滤波和基于机器学习的滤波[7]。
理论上，不同种类滤波方法利用的点云信息也各不相同。大量研究表明，滤波算法利用的信息越多滤波精度越高，表现越稳健。为此，为了利用更多的点云信息，研究人员提出了多种不同点云滤波组合算法。部分学者通过扩充插值滤波的种子点数量，来提高算法对地面点的识别能力。例如，Zhang和Lin[8]在进行迭代加密三角网滤波的过程中，通过对地面种子点的区域增长来加密初始三角网。Zhao等[9]先通过形态学开运算和平移平面拟合法精确获取地面种子点，然后通过种子点构建三角网并迭代加密。Chen等[10]通过区域增长分块获取地面种子点，然后通过基于薄板样条（TPS）插值的分层滤波来识别地面点。另一部分学者通过组合不同的方法以利用不同的几何信息。例如，Yang等[11]先将网格点云划分为一组点云块和一组散点，然后分别通过基于块的滤波和改进的多尺度形态学滤波从点云块和散点中识别地面点。Bigdeli等[12]先通过基于形态学重建的滤波滤除较大地物，然后通过基于坡度的滤波滤除低矮植被等稀疏非地面点。Meng等[13]先通过形态学重建滤波初步确定数字表面模型中的地面区域，然后再对地面区域使用坡度自适应阈值来选择地面点。大量实验表明，现有滤波方法在陡坡和地形断裂线区域的滤波精度仍有待提高。
因此，为了提高滤波性能，更好地保持地形特征，文中提出了一种同时以点和块为基本处理单元的插值滤波方法，主要步骤如下：
(1) 提出的改进区域增长算法将原始点云划分为点云块集和散点集，避免传统区域增长中的过分割和欠分割问题，为基于块的滤波奠定基础。
(2) 两种基元（点和块）由一种基于插值的滤波方法同时处理。相比由不同的滤波方法依次处理，文中方法操作简单，结果可控，而且可以提高滤波方法在不同地形区域的适应性。
(3) 该方法用点到稳健拟合平面的距离代替点到DTM网格的距离。当陡坡区域的待分类点不在相应的DTM网格中心时，该方法可显著降低空间位置误差影响。
-
新方法的技术路线如图1所示，主要包括两个步骤：
Figure 1. Flow chart of the proposed method
(1) 借助区域增长方法将原始点云分割成点云块集和散点集。
(2) 基于多尺度表面插值滤波将块和散点分类为地面点和非地面点。
-
经典的基于平滑性约束的点云分块算法以空间距离、法向量和平面拟合残差为判定条件进行区域生长。但在低矮地物（图2区域A）和茂密植被（图2区域B）区域，基于主成分分析（PCA）方法拟合出的法向量和平面拟合残差不准确，容易将平坦区域的地面点和低矮物体（图2区域A）、陡坡区域的地面点和茂密植被（图2区域B）划分为同一块。
Figure 2. Classical region growing method
针对这一问题，文中提出了一种改进的区域增长算法，其中点法线的估计步骤如下：
(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为散点的个数)。
Figure 3. Improved region growing method
-
对划分好的点云块和散点同时通过多尺度表面插值滤波进行分类。该滤波过程主要包括以下优势：(1)通过改进的多尺度表面插值滤波，解决传统多分辨率层次滤波[6]存在的空间位置误差问题；(2)基于同一种方法同时处理点云块和散点，避免引入多种方法导致的不可控性；(3)结合点云块选择地面点，尽可能保留平滑区域中的不连续地形[11]，同时，通过散点对复杂地形进行单独处理，避免基于块的分类在复杂地形的不稳定性，提高算法对地形细节的判断能力。
传统的多分辨率层次滤波分三个层次进行点云滤波，DTM分辨率和残差阈值随层次递增。在每个层次结构中，将已选定的地面点通过TPS插值生成参考DTM，并根据待分类点与DTM间的距离来选择地面点。
然而，在通过插值DTM模拟地表进行滤波的过程中，如果某点不在网格中心，则会存在空间位置误差问题，且该误差大小取决于该点距离网格中心的水平距离和当前区域的地形坡度。如图4所示，空间位置误差（
$\Delta h $ ）可以表达为$\Delta h = d \times \tan \alpha $ ，其中α表示地形坡度，d表示待分类点距离网格中心的水平距离。由此可见，地形坡度和水平距离越大，空间位置误差越大。
Figure 4. Spatial position error caused by terrain slope
为此，文中在多分辨率层次滤波三层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距离地面的高度差可估计为：
式中：(xi，yi，zi)为Pi的三维坐标。
则点云块Sei中点的高度差可记为Di={de1, de2, ···, deij}（j为点云块Sei中所包含点的个数），点云块集Se的高度差记为DSe={D1, D2, ···, Dk}，散点集的高度差记为DSc={dc1, dc2, ···, dcm}。
对每个点而言，高度差越小是地面点的可能性越大。因此，根据每个点的高度差，可以从点云块集和散点集中选择地面块和地面散点。针对点云块集Se，笔者通过以下规则从中识别地面块集GSe：
式中：dth为随迭代层次变化的高差阈值；Num(Di<dth)为点云块Sei中高度差小于高差阈值dth的点个数；Num(Sei)表示点云块Sei中点的总数。
针对散点集Sc，根据每个散点的高度差单独确定属性，将高度差小于高差阈值的散点添加到地面散点集GSe，即：
当前层次下识别的地面点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所示。
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% Table 1. Optimized parameters and filtering accuracy of the proposed method on the 15 samples
由表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%，该结果证明了文中方法具有较好的滤波性能。
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中的红框标记部分)，但这些点很容易通过人工编辑的方式剔除。
Figure 6. Filtering results of the proposed method on samp11, samp21, samp51 and samp52
-
ISPRS数据集的平均点间距为1~3.5 m，远低于现今设备所采集到的点云密度。因此，文中通过三组不同地形（s1，s2和s3）的高密度实例数据对该算法进一步验证，分别为覆盖低矮植被的崎岖山地（s1，图7(a)）、覆盖不同高度植被的断裂地形（s2，图7(b)）和包含植被和建筑物的斜坡地形（s3，图7(c)）。各个区域面积约为23公顷（1公顷=104 m2），平均点间距约为0.5 m。每组参考数据均先通过滤波软件进行初始滤波，然后对滤波结果进行人工检查和编辑，以保证参考数据的质量。
Figure 7. Digital surface models of the three samples
除了文中方法外，还将多分辨率层次滤波（MHC）[6]、布料模拟滤波（CSF）[15]、渐进加密三角网滤波（PTD）[16]和渐进形态学滤波（PMF）[17]应用于实例数据滤波，并与文中方法结果比较。如表2所示，在三种地形中，文中方法的性能均优于其他方法。就平均总误差而言，文中方法较精度最低的CSF算法提高75%，较精度最高的MHC算法提高40%。
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% Table 2. Filtering accuracy of the proposed method against the 4 filtering algorithms
五种滤波方法对s2数据滤波后生成的DEM如图8所示。其中MHC算法得到的DEM最平滑，但同时也损失了一些地形细节。PMF算法滤除了过多的地面点，导致DEM产生过大畸变。PTD和CSF算法II类误差过大，导致DEM表面过于粗糙。相比之下，文中算法滤波后生成的DEM与参考DEM相似度最高。
Figure 8. (a) Reference data and interpolated DEMs of (b) proposed method, (c) MHC, (d) PMF, (e) PTD and (f) CSF on s2
-
为了对不同地形的机载LiDAR点云精确滤波，文中提出了一种分块插值滤波方法。该方法利用改进的区域增长算法将点云分为点云块集和散点集，并通过加权最小二乘模拟的地形表面来估计点到地表的距离，以解决点不在网格中心时的空间位置误差问题。对ISPRS基准数据集的实验表明，该方法在15个样本中有11个样本的滤波性能优于现有的滤波算法。在三组高密度实例数据实验中，该方法总误差和Kappa亦优于其他四种经典方法。然而，该方法也存在部分缺陷，例如它需要花费较长时间来执行点云分割。此外，该方法需要五个参数，较难确定其在不同地形上的最优值。因此，笔者后续将致力于提高算法在不同地形场景下的计算效率和稳健性。
Interpolation-based filtering with segmentation for airborne LiDAR point clouds
doi: 10.3788/IRLA20200369
- Received Date: 2020-12-20
- Rev Recd Date: 2021-02-25
- Publish Date: 2021-09-23
-
Key words:
- airborne LiDAR point clouds /
- filtering /
- interpolation /
- multi-scale /
- segmentation
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.
DownLoad: