基于DBSCAN聚类和四象限的激光点云林木胸径获取方法与流程

基于dbscan聚类和四象限的激光点云林木胸径获取方法
技术领域
1.本发明涉及城市绿化林木调查测量技术领域，特别是一种基于dbscan聚类和四象限的激光点云林木胸径获取方法。


背景技术：

2.随着社会的发展，人们对生活环境的要求越来越高，而城市绿化的功能也逐步被人们所熟知，这些功能与人们的身心健康、生存环境质量的改善和提高密切相关。城市绿化在城市中已经具有不可替代和估量的作用。而林木在城市绿化中占据着重要地位，为了对城市绿化资源进行维护，必须定期进行城市绿化资源调查，城市绿化资源调查的重点是城市绿化林木资源的调查。对林木资源的调查以测量林木的直径最为关键。在城市绿化工程中，林木直径也直接影响工程造价。
3.随着三维激光扫描技术的发展，三维激光扫描仪被广泛应用到林木胸径的获取工作中，通过对三维激光扫描数据进行数据处理得到林木胸径的激光点云数据，而针对激光点云的林木胸径提取工作目前主要采用人工在计算机及其图形设备(cad)中判读出林木胸径点云所在位置，人工剔除噪声点，然后进行单株计算[基于单站三维激光扫描数据的林木胸径提取方法，柯敏等，浙江林业科技，2019，39(2)，87-91]，这种方法需要大量的人工参与，过滤噪声点，筛选点云数据，效率较低。


技术实现要素：

[0004]
本发明为解决公知技术中存在的技术问题而提供一种基于dbscan聚类和四象限的激光点云林木胸径获取方法，该方法可以自动快速地分离每株林木胸径点云，并分割出单株林木最优胸径点集，然后进行计算，计算精确，智能化程度高，效率高。
[0005]
本发明为解决公知技术中存在的技术问题所采取的技术方案是：一种基于dbscan聚类和四象限的激光点云林木胸径获取方法，包括以下步骤：
[0006]
一)采用站式三维激光扫描仪，采集测区林木激光点云数据；
[0007]
二)对林木激光点云数据进行拼接、裁切、截取，获取测区林木胸径处指定厚度点云，生成林木胸径点云文件；
[0008]
三)采用dbscan聚类方法对不同林木的胸径点云进行分组，获取单株林木胸径处及其支撑棍的点云集合；
[0009]
四)采用dbscan聚类方法进行分类，剔除单株林木周边支撑棍的噪声点，获得单株林木胸径处的点云集合；
[0010]
五)过滤噪声点，获取单株林木胸径最优点云集合；
[0011]
六)以最优点集为点集，以圆方程为基础方程，采用最小二乘计算法，计算该点集的最优圆，该最优圆半径即为林木直径，该最优圆中心即为林木中心。
[0012]
所述步骤四)的方法为，对单株林木及林木支撑棍的点云集合，根据dbscan聚类方法，设定扫描半径为0.04米，最少扫描点数为7点，对单株林木与林木周边支撑棍点云数据
进行分类，当点云分为一类时，则该点集为单株林木点集；当点集分为两类时，则点集点数较多的为单株林木点集；当点集分为三类时，首先判断是否有一类点集的点数是其他两类点集的5倍以上，如果是，则该点集为单株林木点集；如果否，计算三个分类点集中心点，形成三角形，角度最大的点集即为单株林木点集；当点集数量为四个以上时，计算各个类别点集的中心点和所有类别点集的中心点，接近所有类别点集中心点的点集为单株林木点集。
[0013]
所述步骤五)的方法为，对单株林木点云集合求取均值，获取单株林木点云集合均值点，以均值点为中心点，将单株林木点云数据划分至四个象限；计算各象限中与中心点距离最远的点，将各象限离中心点最远的点与中心点按顺序构建三角形；计算各三角形单位面积内的点数，确定单位面积内点数最少的三角形为单株林木胸径圆心所在的三角形；
[0014]
计算单位面积内点数最少的三角形的外接圆，获得外接圆圆心和半径；以外接圆圆心为圆心，以外接圆半径扩大0.03米为新半径，形成外接新圆，计算落入外接新圆的点，然后判断该点集占单株林木胸径点集的比例是否大于80％，如果是，该点集为最优点集；如果否，以中心点为圆心，以外接圆半径为半径，形成中心圆，计算单株林木点集中同时落入中心圆和外接新圆的点，以该点集为新林木点集；
[0015]
以外接圆圆心为初始振动中心，以外接圆圆心和中心点的连线为对角线构建矩形区域，以x/y方向分别移动0.003米，在矩形区域内形成振动中心集合；构建振动中心集合时，扣除距中心点0.01米的点；以振动中心为圆心，以振动中心与中心点的距离为半径构建振动圆，以振动圆半径内外各0.01米构建振动环；
[0016]
在振动中心集合内移动振动中心，构建不同的振动环，计算新林木点集中落入各个振动环内的点，当整个振动中心集合循环一遍后，落入点数最多的振动环为最优振动环，该振动环内的点集为最优点集。
[0017]
所述步骤六)的方法为，计算最优点集的中值，求取所有点坐标与中值的差，然后根据该差值，按最小二乘法，计算出该点集的最优拟合圆，从而获得单株林木的位置和胸径。
[0018]
所述步骤三)的方法为，加载林木胸径点云文件，采用dbscan聚类方法，设定扫描半径为0.4米，最少扫描点数为30点，对整个区域林木数据进行分割，滤除点数少于30点的噪声点，从而分割出每株林木及林木支撑棍的点云数据。
[0019]
综上所述，本发明提出了一种基于dbscan聚类和四象限的激光点云林木胸径获取方法，首先采集林木三维激光点云数据，裁切出林木胸径处点云数据，基于dbscan聚类，对林木数据进行单株分割，对单株林木数据和林木支撑棍数据进行分隔，获得单株林木胸径点云数据，然后基于四象限法，过滤噪声点，获得林木胸径最优点集，最后，利用最小二乘法，计算林木位置及胸径。
[0020]
本发明具有的优点和积极效果是：通过将dbscan聚类方法应用于林木胸径激光点云数据的分类，和单株林木胸径激光点云数据的获取，可以自动快速分离每株林木胸径点云，有助于提高林木胸径计算效率。通过采用四象限法过滤噪声点，获取单株林木最优胸径激光点云，然后进行计算，计算精确，智能化程度高，使用方便，效率高，可以直接得出林木的位置和胸径数值。
附图说明
[0021]
图1为本发明的流程图；
[0022]
图2为三维激光扫描仪扫描的林木激光点云数据示意图；
[0023]
图3为截取的林木胸径激光点云示意图；
[0024]
图4为林木初步分类后不同的单株林木及其支撑棍的点云集合示意图；
[0025]
图5为精细分类后的单株林木点云集合示意图；
[0026]
图6为四象限法去除噪声点后的胸径点集示意图；
[0027]
图7为根据最小二乘法拟合的最优圆示意图。
具体实施方式
[0028]
为能进一步了解本发明的发明内容、特点及功效，兹例举以下实施例，并配合附图详细说明如下：
[0029]
请参阅图1～图7，本发明为解决在cad中人工判读林木胸径点云效率低的问题，而提供一种基于dbscan聚类和四象限的激光点云林木胸径获取方法，包括以下步骤：
[0030]
一)采用站式三维激光扫描仪，采集测区林木激光点云数据
[0031]
采用站式激光扫描仪，采集林木的激光点云，从而获取整个测区林木从地面到树冠的激光点云。
[0032]
使用前，将扫描仪器放置在观测环境中30分钟以上，然后进行连接并通电测试，确认无误后开始作业。扫描作业时，设置测站间的间距为50m，逐次设站扫描获取数据。
[0033]
二)对林木激光点云数据进行拼接、裁切、截取，获取林木胸径处指定厚度点云，生成林木胸径点云文件
[0034]
以外业扫描相邻两站点连线的中垂线为界线，与城市绿化林木两侧边界构成的封闭区域为单站点云数据裁切范围，裁切出矩形单站点云数据，然后将各测站经裁切后的点云数据，利用设站点坐标进行拼接，得到整个测区点云数据；根据整个测区地面激光点云数据，建立测区dtm模型，截取林木胸径处指定厚度点云数据，生成林木胸径点云数据文件。
[0035]
三)采用dbscan聚类方法，对不同林木的胸径点云进行分组，获取单株林木胸径处及其支撑棍的点云集合
[0036]
加载林木胸径点云文件，采用dbscan聚类方法，设定扫描半径为0.4米，最少扫描点数为30点，对整个区域林木数据进行分割，滤除点数少于30点的噪声点，从而分割出每株林木及林木支撑棍的点云数据。
[0037]
四)采用dbscan聚类方法分类，剔除单株林木周边支撑棍的噪声点，获得单株林木胸径处的点云集合
[0038]
对单株林木及林木支撑棍的点云集合，根据dbscan聚类方法，设定扫描半径为0.04米，最少扫描点数为7点，对单株林木与林木周边支撑棍点云数据进行分类，获得单株林木胸径处的点云集合。
[0039]
单株林木及其支撑棍的点云具有以下特点，单株林木胸径的点云一般居中，单株林木胸径的点云较为密集，点数较多；支撑棍的点云居于边缘，点数较少，支撑棍一般为三个，对称分布；
[0040]
根据单株林木及其支撑棍的点云特点，当点云分为一类时，则该点集为单株林木
点集；当点集分为两类时，则点集点数较多的为单株林木点集；当点集分为三类时，首先判断是否有一类点集的点数是其他两类点集的5倍以上，如果是，则该点集为单株林木点集；如果否，计算三个分类点集中心点，形成三角形，角度最大的点集即为单株林木点集；当点集数量为四个以上时，计算各个类别点集的中心点和所有类别点集的中心点，接近所有类别点集中心点的点集为单株林木点集。
[0041]
五)过滤噪声点，获取单株林木胸径最优点云集合
[0042]
对单株林木点云集合求取均值，获取单株林木点云集合均值点，以均值点为中心点，将单株林木点云数据划分至四个象限；计算各象限中与中心点距离最远的点，将各象限离中心点最远的点与中心点按顺序构建三角形；计算各三角形单位面积内的点数，确定单位面积内点数最少的三角形为单株林木胸径圆心所在的三角形；
[0043]
计算单位面积内点数最少的三角形的外接圆，获得外接圆圆心和半径；以外接圆圆心为圆心，以外接圆半径扩大0.03米为新半径，形成外接新圆，计算落入外接新圆的点，然后判断该点集占单株林木胸径点集的比例是否大于80％，如果是，该点集为最优点集；如果否，以中心点为圆心，以外接圆半径为半径，形成中心圆，计算单株林木点集中同时落入中心圆和外接新圆的点，以该点集为新林木点集；
[0044]
以外接圆圆心为初始振动中心，以外接圆圆心和中心点的连线为对角线构建矩形区域，以x/y方向分别移动0.003米，在矩形区域内形成振动中心集合；构建振动中心集合时，扣除距中心点0.01米的点；以振动中心为圆心，以振动中心与中心点的距离为半径构建振动圆，以振动圆半径内外各0.01米构建振动环；
[0045]
在振动中心集合内移动振动中心，构建不同的振动环，计算新林木点集中落入各个振动环内的点，当整个振动中心集合循环一遍后，落入点数最多的振动环为最优振动环，该振动环内的点集为最优点集。
[0046]
采用四象限法，筛选出单株林木圆心所在位置，然后运用占比法、交集法和移动圆心法，过滤噪声点，获得单株林木胸径处的最优点集。
[0047]
六)以最优点集为点集，以圆方程为基础方程，采用最小二乘计算法，计算该点集的最优圆，该最优圆半径即为林木直径，该最优圆中心即为林木中心。
[0048]
为防止计算过程中出现奇异矩阵，计算最优点集的均值，求取所有点坐标与中值的差，然后根据该差值，按最小二乘法，计算出该点集的最优拟合圆，从而获得单株林木的位置和胸径。
[0049]
城市绿化林木调查具有数量多，密度大，通视条件差，调查精度高等特点。站式三维激光扫描仪具有扫描精度高，扫描距离远的特点。采用站式三维激光扫描仪，可以快速获取林木激光点云数据。
[0050]
根据扫描后的激光点云数据，进行相应的裁切、拼接和截取，获取需要的林木激光点云数据。
[0051]
初步分类，这一步骤依赖于计算机程序处理。对林木胸径激光点云数据，根据dbscan聚类方法，设定扫描半径为0.4米，最少扫描点数为30点，将数据分割，分割出单株林木及林木支撑棍的点云集合。
[0052]
精细分类，这一步骤依赖于计算机程序处理。对分割出的单株林木及林木支撑棍的点云集合，基于dbscan聚类方法，设定扫描半径为0.04米，最少扫描点数为7点，对单株林
木与林木周边支撑棍数据分割；并根据单株林木和支撑棍点云特点，筛选出单株林木胸径点云数据。
[0053]
对精细分类后的胸径点云，提出了四象限法过滤掉噪声点，确定最优点集。计算均值作为中心点；根据中心点，将点云分为四象限；然后根据中心点距四象限最远点，按顺序构建三角形，确定点密度最小的三角形为林木中心所在位置，然后按照一定的半径构成外接新圆，然后运用占比判断、交集计算、移动圆心滤除噪声点，获得最优点集。
[0054]
对获得的最优点集，计算均值，对最优点集与均值求差，然后采用最小二乘法求取圆心和直径，圆心坐标与均值相加即得林木位置，直径即为林木胸径。
[0055]
结合数值的实施例：
[0056]
以天津市某道路绿化区域为例，区域位于天津市滨海新区，区域成带状，区域内林木众多、分布密集。其中西侧为乔木，东侧为灌木，本次仅对乔木进行数据采集。乔木间距约为2.5m
×
2.5m，苗木树龄约为3a。
[0057]
步骤1，采用maptek i-site 8820型三维激光扫描仪进行扫描测量，采用trimbler8gps双频接收机采集数据。将扫描仪器放置在观测环境中30分钟以上，然后进行连接并通电测试，确认无误后开始作业。描作业步骤如下：
[0058]
用gps-rtk获取扫描设站点及后视点的三维坐标信息。将三维激光扫描仪架设在设计站位，对仪器进行整平，然后找准三维激光扫描仪的后视点，并将扫描设站点及后视点坐标信息输入三维激光扫描仪控制手簿中，扫描区域角度设置为360
°
，扫描质量设置为中档，开始扫描测量。
[0059]
扫描完成后，检查数据完整性，如有缺失则补充扫描，如没有则换站，直至完成整个测区扫描作业。扫描后测区的林木点云效果如图2所示。
[0060]
步骤2，以外业扫描相邻两站点连线的中垂线为界线，与城市绿化林木两侧边界构成的封闭区域为单站点云数据裁切宽度，裁切长度设置为50米，裁切出矩形单站点云数据，然后将各测站经裁切后的点云数据，利用设站点坐标进行拼接，得到整个区域点云数据；根据测区地面激光点云数据，建立测区dtm模型，截取距地面13米处0.04米厚度点云数据，生成林木胸径点云数据文件。如图3所示，为了清晰美观，示意图只展示4棵林木的点云数据
[0061]
步骤3，初步分类，采用dbscan聚类算法，根据林木之间的距离大于0.4米，设定分类半径为0.4米，最少点数为30点，对整个区域林木数据进行分割，滤除点数少于30点的噪声点，从而分割出每株林木及林木支撑棍的点云数据。如图4所示，展示了初步分类后的4颗林木数据。
[0062]
dbscan(density-based spatial clustering of applications with noise)是一个基于密度的聚类算法。与划分和层次聚类方法不同，它将簇定义为密度相连的点的最大集合，能够把具有足够高密度的区域划分为簇，并可在噪声的空间数据库中发现任意形状的聚类。
[0063]
步骤4，精细分类，设定扫描半径为0.04米，最少扫描点数为7点，对单株林木与林木周边支撑棍数据分类；筛选出每株林木的点云数据。如图5所示，过滤掉了支撑棍的点云数据，获得了单株林木点云数据。
[0064]
单株林木和支撑棍的点云具有以下特点，单株林木胸径的点云一般居中，单株林木胸径的点云较为密集，数量较多，支撑棍的点云居于边缘，数量较少，支撑棍一般为三个
或四个，对称分布；根据单株林木和支撑棍的点云特点，当点云分为一类时，则该点集为单株林木点集；当点集分为两类时，则点集数目较多的为单株林木点集；当点集分为三类时，首先判断是否有一类点集数量在其他两类的5倍以上，是的话则该点集为单株林木点集；否则，计算三个分类点集中心点，形成三角形，角度最大的点集即为单株林木点集；当点集数量四个以上时，计算各个类别点集的中心点和所有类别点集的中心点，接近所有点集中心点的点集为单株林木点集。
[0065]
计算三角形钝角时，需要用到两点间距离公式和余弦定理。
[0066]
两点间距离公式如下：
[0067][0068]
余弦定理：
[0069][0070]
步骤5，对步骤4获得的单株林木点集，求平均值作为中心点(x0,y0)，根据点集中各个点的x/y坐标与中心点(x0,y0)的距离大小关系，将点集分为四象限点集。
[0071]
求取各个象限中距中心点(x0,y0)最远的点，将中心点与各象限距其最远的点按顺序构建三角形。例如中心点与第一象限最远点、第二象限最远点构建三角形，中心点与第二象限最远点、第三象限最远点构建三角形等等。如果某一象限无点，则跳过该象限，例如第二象限无点，则跳过该象限，中心点与第一象限、第三象限最远点构建三角形。
[0072]
根据斜率法，判断三点是否共线，当x2＝x1、x3＝x1，或则三点共线，不构成三角形，否则，则构成三角形。
[0073]
根据海伦公式，计算各个三角形的面积。海伦公式如下：假定三角形边长分别为a、b、c，三角形的面积s可由以下公式求得：
[0074][0075]
根据点与三角形的位置关系，计算单株林木点集落入各个三角形的点的数量。
[0076]
判断点与三角形位置关系的方法如下：假设三角形的三个点按照顺时针(或者逆时针)顺序是a,b,c。对于某一点p，求出三个向量pa、pb、pc,然后计算以下三个叉乘(^表示叉乘符号)：t1＝pa^pb、t2＝pb^pc、t3＝pc^pa，如果t1,t2,t3同号(同正或同负)，那么p在三角形内部，否则在外部。
[0077]
根据各个三角形的面积和落入三角形内的点的数量，计算每个三角形的点密度，点密度最小的三角形为单株林木胸径圆心坐落位置。
[0078]
根据外接圆公式，计算点密度最小的三角形的外接圆。
[0079]
计算外接圆公式如下：设三角形三个顶点为a(x0,y0)、b(x1,y1)、c(x2,y2)，
[0080]
令a1＝x1-x0，b1＝y1-y0，c1＝(a1*a1+b1*b1)/2，
[0081]
a2＝x2-x0，b2＝y2-y0，c2＝(a2*a2+b2*b2)/2，
[0082]
d＝a1*b2-a2*b1
[0083]
则：圆心(xc,yc)坐标为xc＝x0+(c1*b2-c2*b1)/d
[0084]
yc＝y0+(a1*c2-a2*c1)/d
[0085][0086]
以外接圆圆心(xc,yc)为圆心，以外接圆半径r扩大0.03米为新的半径，构建同心圆。
[0087]
根据点与圆心的距离，计算单株林木点集中落入同心圆的点，判断外接新圆内的点的数量占单株林木点集数量的比例，如果占比超过80％，则外接新圆内的点为最优点集。
[0088]
如占比小于80％，则以中心点(x0,y0)为圆心，以外接圆半径r为半径，构建中心圆，计算单株林木点集中同时落入中心圆和外接新圆的点，以该点集为新林木点集。
[0089]
以外接圆圆心(xc,yc)为初始振动中心，以外接圆圆心(xc,yc)和中心点(x0,y0)的连线为对角线构建矩形区域，以x/y方向分别移动0.003米，在矩形区域内形成振动中心集合；构建振动中心集合时，扣除距中心点(x0,y0)0.01米的点。
[0090]
以振动中心为圆心，以振动中心与中心点(x0,y0)的距离为半径构建振动圆，以振动圆半径内外各0.01米构建振动环。
[0091]
在振动中心集合内移动振动中心，构建不同的振动环，计算新林木点集中落入各个振动环内的点。当整个振动中心集合循环一遍后，落入点数最多的振动环为最优振动环，该振动环内的点集为最优点集。
[0092]
步骤6，根据步骤5计算的单株林木的最优点集，计算最优点集的均值，将最优点集中点与均值求差，然后根据最小二乘法计算最优圆，圆心坐标与均值相加即得林木位置，直径即为林木胸径。如图7所示。
[0093]
最小二乘法原理如下。假定圆心坐标为(x0,y0),半径为r，在圆上分布了n个二维坐标点，且坐标点位(xi，yi),i＝1,2,...,n；
[0094]
圆曲线的基础方程为：
[0095]
r2＝(x-x0)2+(y-y0)2ꢀꢀꢀꢀꢀꢀꢀꢀ
(1)
[0096]
展开可得:
[0097][0098]
令
[0099][0100]
即
[0101][0102]
则(2)式变换后的为
[0103]
x2+y
2-ax-by+c＝0
ꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀ
(5)
[0104]
该公式推广到n个点后，写成矩阵形式可得：
[0105]
[0106]
令：
[0107][0108]
则
[0109]
ax＝l
ꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀ
(8)
[0110]
根据最小二乘算法可得：
[0111]
x＝(a
t
a)-1at
l
ꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀ
(9)
[0112]
根据公式(9)，可求得a，b，c，代入到公式(4)中，即可求得圆心坐标和半径。
[0113]
综上所述，本发明涉及一种基于dbscan聚类和四象限的激光点云林木胸径获取方法，采用三维激光扫描仪对城市绿化林木进行扫描，然后进行截取，获得林木胸径激光点云数据，提出了一种改进的林木胸径的快速获取方法，尤其涉及dbscan聚类和四象限过滤噪声点，经过初步分类、精细分类、四象限过滤噪声点、以及最小二乘的计算，具体包括以下步骤：1)采用站式三维激光扫描仪，采集林木激光点云数据；2)对林木激光点云数据进行拼接、裁切、截取，获取林木胸径处指定厚度点云，生成林木胸径点云文件；3)基于dbscan聚类方法，对不同林木的胸径点云进行分组，获取单株林木胸径处及林木支撑棍的点云集合。4)基于dbscan聚类方法，剔除单株林木周边支撑棍的噪声点，获得单株林木的点云集合；5)基于四象限法，过滤噪声点，获取单株林木胸径最优点云集合；6)采用最小二乘计算法，计算林木胸径。
[0114]
尽管上面结合附图对本发明的优选实施例进行了描述，但是本发明并不局限于上述的具体实施方式，上述的具体实施方式仅仅是示意性的，并不是限制性的，本领域的普通技术人员在本发明的启示下，在不脱离本发明宗旨和权利要求所保护的范围的情况下，还可以做出很多形式，这些均属于本发明的保护范围之内。