辅以多特征聚类的复杂地形区点云层次滤波方法

1.本发明涉及一种辅以多特征聚类的复杂地形区点云层次滤波方法。


背景技术：

2.机载激光雷达(light detection and ranging，lidar)技术作为一种主动式的遥感对地观测手段，在大范围实景三维地理信息获取中发挥着支撑性作用。
3.随着机载lidar系统的快速发展，lidar点云已广泛应用于实景三维建设、城市智能管理、生态环境保护等诸多领域。然而，原始lidar点云混杂有地面信息及植被、建筑物、车辆等非地面信息，因此，点云滤波是上述领域应用的前提和关键步骤。
4.近二十年来，国内外学者提出了多种滤波算法。总结起来可归纳为五类：基于形态学的滤波、基于插值的滤波、基于坡度的滤波、基于分块的滤波和基于机器学习的滤波等。其中，相比其他滤波算法，基于插值的滤波整体适用性较好，能有效的滤除非地面点。
5.然而，由于现实场景地形地物的高复杂性、交错性、多态性，以及部分坐落于丘陵的建筑物与地形存在高程相似性，传统插值滤波算法在复杂地形区仍无法有效区分地面点和地物点，极大影响了滤波质量。此外，传统插值滤波算法仅利用了点云的三维空间坐标信息，缺少对回波强度、颜色、回波次数等物理属性信息的充分利用。因此，传统插值滤波算法在地形断裂等复杂环境下滤波精度低，难以有效区分复杂地形下的地物点和地面点。


技术实现要素：

6.本发明的目的在于提出一种辅以多特征聚类的复杂地形区点云层次滤波方法，以解决现有滤波方法在地形断裂等复杂地形区滤波精度低的技术问题。
7.本发明为了实现上述目的，采用如下技术方案：
8.辅以多特征聚类的复杂地形区点云层次滤波方法，包括如下步骤：
9.步骤1.将原始点云数据格网化，依次对每个独立格网内的点云进行多特征聚类，得到多特征聚类结果，将聚类结果映射到三维空间，得到点云簇的三维坐标分布；
10.步骤2.利用顾及断裂地形的准则从点云簇的三维坐标分布中识别出地面簇；
11.步骤3.利用步骤2捕捉到的地面簇中的地面点作为地面种子点，构建初始地面参考面，并结合多尺度层次滤波方法进一步精化原始点云中的地面点。
12.本发明具有如下优点：
13.如上所述，本发明述及了一种辅以多特征聚类的复杂地形区点云层次滤波方法，首先结合点云几何和物理信息多特征聚类，以提高复杂环境中相邻地物/地面的有效划分；接着采用顾及断裂地形的准则识别地面点簇，避免地形细节损失，提升地面点簇识别率；最后利用点到对应插值点的高差代替点到地面参考面网格的高差，克服了陡坡处点云空间位置误差影响。本发明在地形断裂等复杂地形区滤波精度高，有效地区分了复杂地形下的地物点和地面点。
附图说明
14.图1为本发明实施例中辅以多特征聚类的复杂地形区点云层次滤波方法的流程图；
15.图2为地面与地物相互连通区示意图；
16.图3为采样点与其对应的3
×
3邻近网格关系图；
17.图4为多特征聚类结果图；
18.图5为点云簇的三维坐标分布示意图；
19.图6为顾及断裂地形的地面簇识别示意图；
20.图7为传统mhf示意图；
21.图8为本发明实施例中顾及空间位置误差的mhf示意图；
22.图9为沟壑地形的实验数据影像图；
23.图10为断裂地形的实验数据影像图；
24.图11为陡坡地形的实验数据影像图；
25.图12为断裂陡坡地形的实验数据影像图；
26.图13为基于沟壑地形的参考地面点数据生成的数字高程模型(dem)示意图；
27.图14为基于断裂地形的参考地面点数据生成的数字高程模型(dem)示意图；
28.图15为基于陡坡地形的参考地面点数据生成的数字高程模型(dem)示意图；
29.图16为基于断裂陡坡地形的参考地面点数据生成的数字高程模型(dem)示意图；
30.图17为本发明方法与六种代表性算法总误差的对比示意图；
31.图18为本发明方法与六种代表性算法kappa系数对比示意图；
32.图19为本发明方法与六种代表性算法构建的dem的中误差对比示意图；
33.图20为本发明方法与六种代表性算法各个方法构建的dem的平均误差对比示意图；
34.图21为参考dem示意图；
35.图22为proposed算法对plot4滤波后dems示意图；
36.图23为mhf算法对plot4滤波后dems示意图；
37.图24为smrf算法对plot4滤波后dems示意图；
38.图25为csf算法对plot4滤波后dems示意图；
39.图26为msf算法对plot4滤波后dems示意图；
40.图27为sbf算法对plot4滤波后dems示意图；
41.图28为ptd算法对plot4滤波后dems示意图。
具体实施方式
42.针对现实场景中地物的高复杂性、重叠性、多态性的特点，为充分利用点云的几何和物理信息，本发明构建了一种辅以多特征聚类的层次点云滤波方法。
43.其中，点云聚类可以较好地保留各地物/地面边界，从而提高应对复杂环境能力。
44.如图1所示，该辅以多特征聚类的层次点云滤波方法主要包括数据预处理、多特征点云聚类、顾及断裂地形的地面簇识别以及多尺度层次滤波(mhf)四个步骤。
45.该辅以多特征聚类的层次点云滤波方法的实施过程大致如下：
46.首先根据回波信息剔除非末次回波点；然后将点云格网化，并对每个格网内点云多特征聚类；接着采用一种顾及断裂地形的准则识别地面簇，最后使用mhf进一步识别地面点。
47.下面结合附图以及具体实施方式对本发明作进一步详细说明：
48.步骤1.数据预处理以及多特征点云聚类。
49.激光脉冲从空中向地面传播过程中，如果遇到多个反射表面会形成多次回波。
50.通常地面和建筑屋顶等平整表面仅产生单次回波，而植被、建筑物边缘等垂直结构地物会产生至少两次回波。一般地，除最末次回波外，多次回波中其他回波是地面点的概率较低。
51.在点云滤波之前，本实施例首先根据回波信息从原始点云中首先剔除建筑物边缘等非末次回波的点云，以降低模型计算量，从而提高计算效率。
52.传统插值滤波算法在相对简单的场景下表现较好，但在复杂环境下存在准确性低、地形细节损失等难题，它们对高精度点云滤波提出了严峻挑战。
53.点云聚类能够按照特定判别规则将连通的地形/地物点聚类为同一簇，可有效弥补插值滤波算法在地形断裂等复杂场景中的地形细节损失缺陷。
54.传统聚类方法中，dbscan聚类是一种基于密度的空间聚类方法，其利用簇的密度连通性将点云聚类为任意形状的簇，具有不需要提前指定簇的个数，且能够在聚类的同时发现异常点的优势。基于此，本发明采用dbscan聚类方法对点云聚类。
55.首先将原始点云进行格网化处理。
56.理论上，dbscan在处理高密度、大尺度点云数据时，存在运算效率低、稳健性差、聚类效果不佳等问题，无法满足大范围点云同质点的聚合需求。
57.因此，本实施例首先将原始点云数据格网化(格网大小设为ws)，然后依次对每个独立格网内的点云进行多特征聚类的方式，此种方式的优势在于：
58.由此可避免大范围场景内不同类别点特征相似而导致的聚类混乱问题，进而实现同类特征点云的高效聚合；而且，由于仅对邻域点操作，可极大提升聚类效率。
59.多特征点云聚类。
60.通常，机载lidar点云点密度变化不大，且建筑物和地面会被树木点云连通，如图2所示，而且受地形起伏影响，高程相同的两个点可能分别对应地面和地物点，以及同一类型的地物点在不同位置的高程也不尽相同。
61.因此，仅利用三维坐标信息对复杂环境中点云聚类容易导致错误的聚类结果。
62.基于此，本发明提出结合点云几何和物理信息进行多特征点云聚类。
63.首先利用最小二乘法拟合地面参考面，并计算每个点的拟合残差作为点云的辅助特征。
64.为了准确描述地形表面，以每个格网内最低点作为初始地面点，并利用9个邻近格网的初始地面点，拟合切平面来计算格网内每个点的拟合残差。
65.如图3所示，从dtm格网中取一采样点，该点pi位于dtm表面的(i，j)格网中。
66.其中，(i，j)表示横坐标为i，纵坐标为j的格网。
67.则它的9个最邻近格网为(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)。
68.设点pi对应的切平面方程为z＝ax+by+c，则pi的拟合残差可估计为：
69.ri＝z
i-(axi+byi+c)
ꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀ
(1)
70.其中，(xi，yi，zi)为点pi的三维坐标，a、b、c为平面方程系数。
71.由于稠密的低矮植被与地面具有相似的高程及拟合残差信息，在点云聚类时低矮植被易与地面点混杂为一簇。鉴于同一地物常具有相似或连续的rgb颜色信息，且不同地物间颜色变化较大，本发明还将颜色信息作为辅助特征用于点云聚类。
72.为简化计算且突出地物与地面差异，将rgb颜色信息转化为单一值的可见光波段差异性植被指数(visible-band difference vegetation index,vdvi)。
73.其中，vdvi的计算公式为：
74.式中，r、g、b分别表示点云的红、绿、蓝颜色值。
75.本实施例以矩阵[高程，拟合残差，vdvi]作为点云虚拟坐标(z,h,vdvi)，执行基于欧式距离的dbscan聚类，其中，高程即空间几何信息中的z坐标。
[0076]
本实施例中dbscan聚类的具体步骤如下：
[0077]
步骤1.3.1.将原始点云数据中所有点标记为未分类点，计算未分类点中任意两点间的欧氏距离：其中，聚类过程中任意两点间的欧式距离计算公式如下：
[0078][0079]
其中，d
ij
表示i点与j点之间的欧氏距离。
[0080]
zi、zj分别表示i点与j点的高程，ri、rj分别表示i点与j点的拟合残差值，vdvii、vdvij分别表示i点与j点的vdvi值，α,β,γ是比例因子系数，α,β,γ均大于0。
[0081]
步骤1.3.2.从未分类点中随机选择一点p，同时将该点标记为已分类点，然后搜索该点的ε邻域n，并计算n内点的个数；若集合n中的点个数大于最小点个数p
ts
，则将p标记为核心点，同时创建一个新簇c，并把p点添加到簇c中。
[0082]
其中，ε表示搜索邻域点的半径；n为邻域内点的集合；p
ts
表示最小点个数。
[0083]
步骤1.3.3.遍历n中所有未分类点q的ε邻域m，若m内点个数大于p
ts
，则将m中点添加到n中，同时若点q不属于任何簇，则将q添加到簇c中。
[0084]
步骤1.3.4.重复步骤1.3.3，直至没有新的点添加到簇c，则输出簇c，并返回步骤1.3.2；若无新簇输出，则将剩余点归属为噪声点，直至没有标记为未分类的点，则进入步骤1.3.5。
[0085]
步骤1.3.5.输出结果；簇划分为c
′
＝{c1,c2,c3,c4,
……
,ck}。
[0086]
其中，c
′
表示点云簇的集合，c1,c2,c3,c4,
……
,ck分别表示划分的点云簇。
[0087]
多特征聚类后的结果如图4所示，由图4聚类后的结果可以看出，该多特征聚类充分融合了高程、拟合残差、vdvi等多个特征，能够有效聚合特征相似的点云簇。
[0088]
然而，仅仅依靠多特征聚类结果，却无法判定簇的类别(地面或非地面簇)。因此，需要将聚类结果映射到三维空间如图5所示，得到点云簇的三维坐标分布。
[0089]
然后基于点云簇的三维坐标分布，并根据其几何特征判断点云簇的类别。
[0090]
步骤2.顾及断裂地形的地面簇识别。
[0091]
传统方法一般根据点云簇曲率特征和是否包含格网最低点以判断其类别，即将格
网最低点所属的簇或将点云簇曲率σ小于点云簇阈值σ
th
的簇判定为地面簇。
[0092]
其中，点云簇曲率σ的计算过程如下：
[0093][0094]
式中，λ1,λ2,λ3(λ1＜λ2＜λ3)为pca方法计算的点云簇特征值。
[0095]
然而，上述方法容易丢失地形特征处的地面簇或将非地面簇误判为地面簇。
[0096]
如图6所示，传统方法由于地面簇b与其所在格网最低点g2高差较大，被错分为非地面簇，而非地面簇c因曲率较小被误判为地面簇。
[0097]
为此，除了包含格网最低点的簇直接指定为地面簇外，本实施例将点云簇与8个邻近格网最低点的高差、邻近格网边缘的水平距离以及曲率特征准则识别地面簇。
[0098]
即地面簇gp必须满足：
[0099]
式中，σi为第i个点云簇曲率，σ
th
为点云簇曲率阈值。
[0100]
cpi表示第i个点云簇，为第i个点云簇平均高程，h
gj
表示第j个邻近格网最低点高程，其中，i＝1,2,
…
,j＝1,
…
,8，h
th
表示预设高差阈值。
[0101]
li为第i个点云簇中点到邻近网格边缘距离，l
th
表示预设距离阈值。
[0102]
更为具体地，顾及断裂地形的地面簇识别过程如下：
[0103]
首先将点云簇曲率小于σ
th
的点云簇判断为候选地面簇。
[0104]
然后计算该点云簇平均高程与8个邻近格网最低点的高差。
[0105]
若经过计算得到的8个高差中任意一个高差值小于预设高差阈值h
th
，则进一步判断该点云簇内点到该邻近格网边缘距离的最小值l
min
(例如图6中的lb和lc)。
[0106]
若l
min
小于距离阈值l
th
，则该点云簇判别为地面簇，反之，该点云簇为非地面簇。
[0107]
如图6中，由于点云簇b内点到该邻近格网边缘距离的最小值lb小于距离阈值l
th
，该点云簇b判别为地面簇，而lc大于距离阈值l
th
，因此，点云簇c判别为非地面簇。
[0108]
以上准则既可准确识别地形断裂区的地面簇(图6中的点云簇b)，又能避免将曲率较小的非地面簇(图6中的点云簇c)误判为地面簇，极大提高了地面簇的识别率。
[0109]
步骤3.多尺度层次滤波。
[0110]
尽管顾及地形断裂的地面簇能够准确识别大量地面点，但由于点云场景的复杂性，部分地形细节特征由于无法聚类为点云块而被忽略。
[0111]
因此，为进一步识别细节地面点，本实施例利用捕捉到的地面簇中的点作为地面种子点构建初始地面参考面，并在此基础上执行多尺度层次滤波(mhf)。
[0112]
传统mhf采用待分类点到地面参考面网格垂直距离作为准则识别地面点。但当网格分辨率较粗且该网格处地形坡度较大时，该准则受空间位置误差显著影响如图7所示。
[0113]
为此，本实施例提出了一种顾及空间位置误差的mhf，具体过程如下：
[0114]
步骤3.1.首先将地面簇中的地面点作为地面种子点，构建初始地面参考面。
[0115]
步骤3.2.每一层用薄板样条插值构建分辨率为r的地面参考面，并利用四个邻近地面参考面网格中心点对非地面簇中的点双线性插值计算该点插值高程。
[0116]
将该插值高程与该点高程做差，若差值小于预设高差阈值，则判断为地面点。
[0117]
步骤3.3.将滤波出的地面点用于更新步骤3.2中的地面参考面，直到该层没有地
面点加入，更新后的地面参考面用于下一层滤波，重复上述步骤3.2直至滤波完成。
[0118]
本实施例中优选采用三层滤波，当然也可以不限于三层滤波，需要说明的是，随着本实施例中滤波层次的增加，地面参考面的分辨率以及预设高差阈值也逐渐提高。
[0119]
下面通过实验对本发明方法的有效性作进一步详细说明。
[0120]
目前，大部分算法所用的评估数据为国际摄影测量与遥感协会(isprs)发布的基准数据集，但该数据集仅包含点云几何三维坐标和分类信息。
[0121]
随着点云采集设备的发展，点云数据密度急剧提升，同时将激光雷达扫描的三维点云数据与二维影像数据高精度配准可获得具有物理信息的彩色点云。然而isprs基准数据集点密度远低于当前生产的点云数据，且缺乏颜色、回波次数等其他辅助物理信息。
[0122]
因此，本发明实施例选用四组地形复杂、建筑物和植被混杂的高密度机载lidar点云数据为研究对象，将新方法滤波结果与六种代表性滤波算法比较，包括：
[0123]
多尺度层次滤波(mhf)、改进简单形态学(smrf)、布料滤波(csf)、坡度滤波(sbf)、渐进加密三角网滤波(ptd)和移动曲面拟合滤波(msf)。
[0124]
四组实验区域面积均为500
×
500m2，包含多种地形特征，如沟壑地形、断裂地形、陡坡地形以及断裂陡坡地形，分别如图9至图12所示。
[0125]
为了保证参考点云数据的质量，首先使用terrascan软件自动获取初始滤波结果，再通过人工检核的方式剔除错误分类点。
[0126]
根据参考地面点数据生成的数字高程模型(dem)，如图13至图16所示。表1详细描述了4组实验数据的地形特征、地面点数、非地面点数和平均地物高度等统计信息。
[0127]
表1实例数据集的统计信息
[0128][0129]
为了对滤波结果定量评价，本发明以混淆矩阵(i类误差t.i、ii类误差t.ii和总误差t.e)以及kappa系数k作为精度指标。四个评价指标的计算公式为：
[0130]
t.i＝b/(a+b)
×
100％
ꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀ
(6)
[0131]
t.ii＝c/(c+d)
×
100％
ꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀ
(7)
[0132]
t.e＝(b+c)/e
×
100％
ꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀ
(8)
[0133]
k＝(p
0-pc)/(1-pc)
×
100％
ꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀ
(9)
[0134]
其中，a表示正确分类的地面点数；b表示误分为非地面点的地面点数；c表示误分为地面点的非地面点数；d表示正确分类的非地面点数。
[0135]
e＝a+b+c+d；p0＝(a+b)/e；pc＝((a+b)
×
(a+c)+(c+d)
×
(b+d))/e2。
[0136]
为评价算法对地形的还原能力，还进一步计算了各个滤波方法的dem误差。
[0137]
首先通过地形转栅格法对滤波后的地面点插值生成dem，然后计算各方法dem相对于参考dem(如图13至图16)的中误差(rmse)和平均误差(mae)：
[0138][0139][0140]
式中，ei为参考dem与插值dem在第i个格网的高程差，n是dem格网数。
[0141]
本实施例中新方法涉及五个参数：地面参考面初始分辨率r、最大窗口分辨率ws、初始高差阈值d
th
、识别地面点簇的高差阈值h
th
和水平距离阈值l
th
。
[0142]
为了测试新方法的最佳性能，该实施例中通过联合调试的方法获取四组实验数据对应的最优参数以及滤波结果精度如表2所示。
[0143]
表2四组实例数据的最优参数和滤波精度
[0144][0145]
结果分析表明：新方法在plot1取得了最小的总误差(3.59％)及最高的kappa系数(92.58％)，这是因为plot1地形平缓，植被稀疏，滤波难度较低；而在plot4区域受地形坡度较大及附着低矮植被影响，滤波精度略低于其他区域(总误差为4.36％，kappa系数为87.82％)。平均而言，新方法的总误差和kappa系数分别为4.03％、90.39％，且标准差为0.41和2.26，证明了新方法在处理各类地形时效果稳健。相比plot1和plot2，新方法在plot 3和plot 4的i类误差大于ii类误差，这是由于该区域地面点数远少于非地面点数(表1)，少量错分的地面点容易导致较大的i类误差。然而，新方法的平均ii类误差与i类误差仅相差0.36％，表明新方法在减小i类误差的同时有效抑制了ii类误差。
[0146]
新方法与六种代表性算法总误差和kappa系数对比表明，如图17和图18所示。
[0147]
从对比结果来看，本发明方法的滤波精度及稳定性最优。从总误差来看，新方法在不同场景下滤波总误差均最低，且误差浮动范围最小，csf次之，sbf最差，表明新方法总误差稳定性最好。由kappa系数可知，新方法kappa值最高，smrf次之，sbf最差。平均而言，新方法总误差为4.03％，较精度最差的msf算法提高了66.7％，较精度次高的smrf算法提高了20.3％；kappa系数为90.39％，较其他六种方法至少提高了3.0％。
[0148]
七种方法在不同地形下构建的dem精度表明如图19和图20所示。通过对比发现，新方法在所有场景下生成的dem精度均最高。其中，新方法在plot2中精度明显高于其他方法。相比其他区域，各个方法对plot4构建dem的中误差和平均误差较大，这是由于该区域地面点稀疏，少量漏分的地面点就会导致较大的地形损失。平均而言，新方法生成dem的中误差较精度最差的ptd降低了71.5％，较精度最好的mhf降低了31.5％；平均误差较精度最好的mhf降低了39.3％，较精度最差的sbf降低了78.2％。
[0149]
七种滤波方法对plot4滤波后生成的dem如图21至图28所示。该区域地形复杂，包
含大量低矮植被和建筑物，且西南角存在明显的地形断裂，滤波难度较大。
[0150]
结果表明，本发明方法生成的dem与参考dem最为相近，准确地滤除了复杂地形处建筑物，在断裂地形和局部凸起区域对地面点识别率较高，只有在陡坡上存在较少的异常凸起(如图22中的实线椭圆部分)。这可能是由陡坡上少量的低矮植被与地面点特征相似而被错分为地面点导致。相比之下，mhf生成的dem在陡坡地形表面较为粗糙，同时对断裂地形的保持能力欠佳(如图23虚线椭圆部分)。smrf保留了过多的非地面点，导致dem粗糙(图24实线椭圆部分)，但在断裂地形上损失了部分地形细节(图24虚线椭圆部分)。csf在斜坡上有较大的i类误差，即会滤除过多的地形细节，导致地形平滑(图25虚线椭圆部分)。msf、sbf和ptd(图26-28)均有较大的ii类误差，导致dem表面粗糙。
[0151]
为进一步验证新方法的创新点对滤波结果影响，本发明分别以plot2、plot3和plot4为例，对比加入和不加入相应创新点(顾及地形断裂的地面点簇识别、多特征点云聚类、顾及空间位置误差的点云滤波)时新方法的滤波精度，如表3所示。
[0152]
表3加入与不加入创新点时，新方法滤波精度比较(％)
[0153][0154][0155]
以plot2为例，本发明验证了顾及和未顾及断裂地形的地面点簇识别准则对断裂地形细节的保持效果。结果表明(表3)，与未顾及断裂地形的方法相比，采用顾及断裂地形的判别准则选择地面点簇时i类、ii类和总误差分别减小了74.3％、25.6％、52.5％，kappa系数提高了11.0％。借助plot3对比使用与不使用多特征点云聚类的滤波效果表明(表3)，前者较后者i类、ii类和总误差分别减小了17.9％、8.8％、11.6％，kappa系数提高了1.7％。借助plot4对比新方法考虑和不考虑空间位置误差的滤波效果表明(表3)，与未考虑空间位置误差方法相比，改进后方法的i类误差稍微增加，但ii类和总误差分别减小了23.7％和11.0％，kappa系数提高了1.4％。借助四组实验数据将新方法与六种传统滤波方法结果对比表明，本发明方法在所有实验区滤波性能最优，其平均kappa系数为90.39％，平均总误差为4.03％，比其他六种代表性滤波方法精度至少提升了20.3％；同时，新方法的平均i类与ii类误差仅相差0.36％，体现了该方法处理复杂环境点云的稳健性。
[0156]
综上所述，新方法的平均总误差最小，滤波性能最优，稳定性最高。
[0157]
当然，以上说明仅仅为本发明的较佳实施例，本发明并不限于列举上述实施例，应当说明的是，任何熟悉本领域的技术人员在本说明书的教导下，所做出的所有等同替代、明显变形形式，均落在本说明书的实质范围之内，理应受到本发明的保护。