基于地基激光雷达技术林木胸径值提取的方法

本发明属于林业工程技术领域，尤其涉及一种基于地基激光雷达技术林木胸径值提取的方法。

背景技术：

胸径是指距离根茎1.3m处的树干截面直径，是评价林木生长状态的重要参数之一，其准确提取对研究森林生物量估计、林业资源监测及动态管理等具有重要意义。传统的胸径测量方法主要采用胸径尺或检径尺人工实地接触式测量，该方法操作简单，但存在作业时间长、劳动强度大等不足，不能满足现代林业信息化、智能化发展要求。在传感器和国家需求的双重驱动下，高分辨率精细遥感技术得到了快速的发展，涌现了基于无人机摄影测量、机载激光雷达、地基激光雷达等多遥感技术的森林参数反演方法，其中地面激光雷达技术采用非接触式高速激光测量方法，在不损伤树木的前提下直接获取林下空间高精度、高密度的三维空间和反射强度信息—点云数据，具有其他遥感技术无法比拟的优势，为森林资源调查和林业研究提供了全新的数据获取手段，已成为林木参数快速提取研究的热点。

目前，针对地基激光雷达点云数据林木胸径提取，国内外学者开展了大量研究，可概括为3个步骤：①胸径切片提取。即从庞大的样地点云数据集中分离得到林木1.3m处一定厚度点云。②切片点云分割。单木胸径切片点云是胸径值计算最小单元，步骤①提取得到的点云包含多个单木胸径切片共存现象，切片点云分割就是将相同林木胸径点聚为一个集合，不同林木胸径点彼此分开，为后续胸径值计算做准备。③胸径值估算。以圆或者椭圆为胸径截面模型，通过求解模型参数(半径r/或长短半轴/a,b)，从而得到胸径值估算结果。较比于步骤①③，切片点云分割研究相对较为薄弱，主要存在以下两个方面的不足：①存在人为干预，自动化程度有待提高。通常需要林木数等样地调查先验知识，降低了胸径提取过程的自动化程度。②单木胸径点提取精度不高。分割处理后的单木胸径切片点云还存在噪声点，这些噪声点离胸径点距离很近，常规方法很难去除，其将作为有效样本点参与胸径值估算，严重影响胸径提取结果。

由此可见，现有的林木胸径提取方法存在胸径切片分割精度低、自动化程度不高的问题。

技术实现要素：

本发明实施例的目的在于提供一种基于地基激光雷达技术林木胸径值提取的方法，旨在解决现有的林木胸径提取方法存在胸径切片分割精度低、自动化程度不高的问题。

本发明实施例是这样实现的，一种基于地基激光雷达技术林木胸径值提取的方法，包括：

获取地基雷达林地点云数据；

对所述地基雷达林地点云数据进行点云高程化处理，提取胸径切片点；

根据改进k均值聚类对所述胸径切片点进行切割处理，提取单木胸径切片点；

对所述单木胸径切片点进行椭圆拟合处理，获得胸径值提取结果。

本发明实施例的另一目的在于一种基于地基激光雷达技术林木胸径值提取的装置，包括：

林地点云数据获取单元，用于获取地基雷达林地点云数据；

胸径切片点提取单元，用于对所述地基雷达林地点云数据进行点云高程化处理，提取胸径切片点；

单木胸径切片点提取单元，用于根据改进k均值聚类对所述胸径切片点进行切割处理，提取单木胸径切片点；以及

胸径值提取结果获得单元，用于对所述单木胸径切片点进行椭圆拟合处理，获得胸径值提取结果。

本发明实施例的另一目的在于一种计算机设备，包括存储器和处理器，所述存储器中存储有计算机程序，所述计算机程序被所述处理器执行时，使得所述处理器执行所述基于地基激光雷达技术林木胸径值提取的方法的步骤。

本发明实施例的另一目的在于一种计算机可读存储介质，所述计算机可读存储介质上存储有计算机程序，所述计算机程序被处理器执行时，使得所述处理器执行所述基于地基激光雷达技术林木胸径值提取的方法的步骤。

本发明实施例提供的基于地基激光雷达技术林木胸径值提取的方法，通过对地基雷达林地点云数据进行点云高程化处理，提取胸径切片点，根据改进k均值聚类对所述胸径切片点进行切割处理，提取单木胸径切片点，对所述单木胸径切片点进行椭圆拟合处理，获得胸径值提取结果。本发明通过改进后的k均值聚类能快速的实现林木胸径点批量化提取，无需林木数目、样地大小等先验知识，具有自动化程度高，抗噪性强的优势，对研究地基激光雷达在林业资源调查及生产管理应用具有一定的实际应用参考价值。

附图说明

图1为本发明实施例提供的一种基于地基激光雷达技术林木胸径值提取的方法的流程图；

图2为本发明实施例提供的另一种基于地基激光雷达技术林木胸径值提取的方法的流程图；

图3为本发明实施例提供的“拐点”法k值确定仿真实验示意图：(a)模拟实验数据；(b)不同k值与目标函数关系曲线；

图4为本发明实施例提供的又一种基于地基激光雷达技术林木胸径值提取的方法的流程图；

图5为本发明实施例提供的再一种基于地基激光雷达技术林木胸径值提取的方法的流程图；

图6为本发明实施例提供的还一种基于地基激光雷达技术林木胸径值提取的方法的流程图；

图7为本发明实施例提供的不同参数求解方法的圆拟合结果示意图：(a)最小二乘法；(b)随机采用一致性算法；(c)ls-ran法；

图8为本发明实施例提供的基于地基激光雷达技术林木胸径值提取的装置的结构框图；

图9为本发明实施例提供的测站点分布及样地点云数据示意图：(a)样地测站点分布(三角形标识)；(b)样地点云侧视图；

图10为本发明实施例提供的样地实景及实测数据获取示意图：(a)胸径测量方式；(b)样地林木编号；

图11为本发明实施例提供的改进均值聚类单木胸径切片分割示意图：(a)胸径切片点提取结果；(b)切片点分割结果；

图12为本发明实施例提供的编号9单木提取误差测量误差来源示意图；

图13为本发明实施例提供的胸径提取结果与实测值线性回归分析示意图。

具体实施方式

为了使本发明的目的、技术方案及优点更加清楚明白，以下结合附图及实施例，对本发明进行进一步详细说明。应当理解，此处所描述的具体实施例仅仅用以解释本发明，并不用于限定本发明。

在本发明实施例中使用的术语是仅仅出于描述特定实施例的目的，而非旨在限制本发明。在本发明实施例和所附权利要求书中所使用的单数形式的“一种”和“该”也旨在包括多数形式，除非上下文清楚地表示其他含义。还应当理解，本文中使用的术语“和/或”是指并包含一个或多个相关联的列出项目的任何或所有可能组合。

应当理解，尽管在本发明实施例中可能采用术语第一、第二等来描述各种信息，但这些信息不应限于这些术语。这些术语仅用来将同一类型的信息彼此区分开。

林木受争夺空间资源等竞争关系影响，其生长状态存在自疏现象。林木胸径切片在点云空间中表现出典型的空间特性：单木胸径切片点云紧密相连成类圆形并聚集为一簇，不同单木胸径切片点云彼此分开且保持一定的距离。因此，可引入效率高、收敛速度快的k均值聚类算法实现单木胸径切片点云分割。首先根据归一化高程点云数据截取胸径切片点，然后通过改进k均值聚类对胸径切片点进行分割，提取单木胸径切片点，最后对分割得到的离散单木胸径切片点拟合圆并计算单木胸径值。

胸径是表征立地质量和林木生长状态的重要参数之一。本发明实施例为了解决现有胸径提取精度不够高、自动化程度低等不足的问题，基于地基激光雷达胸径切片点云数据，提出了一种改进k均值聚类的基于地基激光雷达技术林木胸径值提取的方法。首先利用约束条件优化初始种子点的选择原则，避免随机种子点选择造成聚类结果陷入局部最优；其次采用拐点法自适应确定聚类目标类别数目，提高单木胸径点分割的自动化程度；最后根据反距离加权计算聚类中心并剔除非目标单木点集，通过圆模型参数求解实现胸径值计算。本发明通过改进后的k均值聚类能快速的实现林木胸径点批量化提取，无需林木数目、样地大小等先验知识，具有自动化程度高，抗噪性强的优势，对研究地基激光雷达在林业资源调查及生产管理应用具有一定的实际应用参考价值。

如图1所示，在一个实施例中，提出了一种基于地基激光雷达技术林木胸径值提取的方法，具体可以包括以下步骤：

步骤s101，获取地基雷达林地点云数据。

步骤s102，对所述地基雷达林地点云数据进行点云高程化处理，提取胸径切片点。

在本发明实施例中，林木通过地面连通为一个整体，受地表起伏的影响，不同林木根茎处的高程点不同，即胸径切片点云高程起算面不同，为了方便胸径切片点云批量化自动提取，采用网格高程归一化处理将根茎点高程转换到统一高程面内。

在本发明实施例中，具体操作如图2所示，所述步骤s102，包括：

步骤s201，将林地划分为边长为固定值的正方形网格，索引网格内最低点视为地面点，确定所述正方形网格内的所有点高程与地面点高程的差值。

在本发明实施例中，将样地划分边长为固定值的正方形格网，索引网格内最低点视为地面点，将网格内的所有点高程与地面点高程作差。

步骤s202，根据所述地基雷达林地点云数据中的任意点、所述正方形网格内的所有点高程与地面点高程的差值以及胸径处点云切片厚度，提取胸径切片点。

在本发明实施例中，胸径切片点集s1.3可描述为高程1.3m处一定厚度范围内的点云如公式(1)所示，其中pi为样地点云数据m中任意点，ε为胸径处点云切片厚度，通常情况下取胸径尺的宽度3.0～5.0cm。

s1.3＝{pi＝(xi,yi,zi)∈m|zi-1.3|≤ε}(1)

步骤s103，根据改进k均值聚类对所述胸径切片点进行切割处理，提取单木胸径切片点。

在本发明实施例中，k均值聚类(k-meansclusteringalgorithm)是一种迭代求解的聚类分析方法，具有效率高，操作方便等优点，已广泛应用于市场调查、数据挖掘等领域。其基本思想为：对于一个未标记的数据集，在给定类别数k的前提下，通过随机选择初始种子点即初始目标类簇中心，以对象-类簇中心的欧式距离作为度量指标确定对象类别归属，并把类簇中所有对象的平均值作为新目标类簇的中心，不断迭代直至类簇中心变化范围趋于稳定或达到预期迭代次数，最终产生k个目标类簇子集。k均值聚类的目标是使得类簇内样本间距离尽可能小，不同样本间的距离尽可能大，即k均值聚类算法的目标函数可表示为簇内样本间距之和如公式(2)所示。其中k为目标类簇数，n为样本点总数，xi为属于第j个目标类簇的样本，uj表示第j个目标类簇的中心。当目标函数j取最小值时，算法收敛此时有最佳聚类结果。

原始k均值聚类算法需要目标类别数目k为先验输入，一定程度上降低了算法的自动化程度，同时受随机选择初始种子点的影响，算法容易陷入局部最优产生错误的聚类结果。结合林木胸径切片点云空间分布特性，本发明通过以下三个方面的改进消除k均值聚类在单木胸径提取中的缺陷，以提高林木胸径提取的精度和自动化程度。

(1)优化初始种子点选择原则。根据单木胸径切片点云空间上相互隔离并保持一定距离，且距离远大于胸径值的特性，以林木间距最小值为约束条件选择初始种子点，即如果初始种子点两两距离最小值小于林间距离最小值，则满足条件继续下一步，否则，重新生成初始种子点。通过优化种子点的选择原则，保证初始种子点来源于不同的单木胸径切片点云集，提高了算法聚类精度和收敛速度。

拐点法自适应k值确定。拐点法就是在不同的目标聚类数下计算目标函数j值，然后通过检测变化率的突变确定“拐点”，则突变点对应的k值即为最佳目标类簇数。其原理可解释为随着目标类簇数目的增大，目标类簇集中样本数量逐渐减少，样本与类簇中心的距离更近，通过公式(2)计算则目标函数j值也随之变小。当目标聚类数大于最佳类别数目k时，类簇子集由多个类簇间合并转变为单个类簇分裂形成，此时目标函数j陡然变小并随之保持变缓，因此目标函数的拐点即为最佳聚类数。模拟产生一组试验数据以验证“拐点”法的可行性和有效性，试验数据由三个二维正态分布离散点叠加而成如图3a所示，每个正态分布点集代表一个目标类簇，目标函数j与不同目标聚类数k的对应关系曲线如图3b所示。通过拐点法确定最佳k值为3，其与模拟数据预先设置类别保持一致。

(3)非目标自动识别。由于林地植被生长速率不同，经高程归一化切割提取得到的胸径切片点云中还包含灌木、胸径小于5.0cm的非目标单木，可改进k均值聚类后目标几何中样本与类中心距离分布进行识别：采用反距离加权方法确定簇类中心位置以减弱噪声点的影响，同时计算类簇中样本与中心点的平均距离vd和方差sd，如果vd小于2.5cm或者sd大于方差阈值t，则标记为非目标类簇。公式(3)为反距离加权公式，其中ci为第i类的聚类中心，pj为第i类中的第j个对象点，sj，smax，smin分别代表第j个点与上次迭代聚类中心的距离、距离最大值和最小值。

改进后的k均值聚类算法能根据样本数据集的特征自适应确定目标聚类数k，同时通过林木间距为约束条件使得初始种子点的选择更为合理，避免了迭代计算陷入局部最优，改进后的k均值聚类胸径切片点云分割步骤如图4所示，即所述步骤s103，包括：

输入：胸径切片点云，方差阈值t

步骤s401，根据k均值聚类算法的目标函数计算不同k取值时目标函数值，并通过目标函数变化率确定拐点以确定最佳目标类别数目k值。

在本发明实施例中，根据公式(2)计算不同k取值时目标函数值，通过目标函数变化率寻找拐点，自适应确定最佳目标类别数目k值。

步骤s402，随机选择一组初始种子点并确定种子点间距最小值。

步骤s403，判断所述种子点间距最小值是否小于林木间距最小值；若否，则返回至所述步骤s402；若是，则进入步骤s404中。

在本发明实施例中，随机选择一组初始种子点并确定种子点间距最小值d，如果d小于林木间距最小值，则进入步骤s404，否则重新生成种子点，即返回至所述随机选择一组初始种子点并确定种子点间距最小值的步骤。

步骤s404，根据所述初始种子点，以距离最小化为分类原则，对所有样本点进行归类。

步骤s405，根据反距离加权公式计算类簇内样本点的重心，并更新为新的种子点集合，判断所述种子点集合是否发生变化；若是，则返回至步骤s404；若否，则进入步骤s406中。

在本发明实施例中，按公式(3)计算类簇内样本点的重心，并更新为新的种子点，如果种子点集合不再变化，则进入步骤s406，否则返回至所述根据所述初始种子点，以距离最小化为分类原则，对所有样本点进行归类的步骤。

步骤s406，分别计算类簇几何内样本点与类中心距离的平均距离和方差，以根据所述平均距离和方差确定单木胸径切片点。

在本发明实施例中，分别计算类簇几何内样本点与类中心距离的平均距离vd和方差sd，如果vd小于2.5cm或sd大于t，则将该目标类簇标记为非目标集合，否则标记为单木胸径切片点集。

输出：单木胸径切片点集和非目标类簇点集

在本发明实施例中，通过拐点法确定目标聚类簇数提高k均值聚类的自动化程度，采用优化初始种子点选择原则避免算法陷入局部最优，利用预估胸径值剔除非目标单木点集，从而提高单木胸径提取精度。

步骤s104，对所述单木胸径切片点进行椭圆拟合处理，获得胸径值提取结果。

在本发明实施例中，如图5所示，所述步骤s104，包括：

步骤s501，根据所述单木胸径切片点的胸径中心点坐标以及胸径切片半径，构建胸径截面形态模型。

在本发明实施例中，胸径截面形态可用二维平面圆方程近似表达，其数学表达式如公式(4)所示。其中(a,b)表示胸径中心点坐标，r表示胸径切片半径。则胸径提取值d＝2r。

(x-a)²+(y-b)²＝r²(4)

步骤s502，联合随机采样一致性算法和最小二乘法求解所述胸径截面形态模型的模型参数，以获得胸径值提取结果。

在本发明实施例中，一般情况下，模型参数求解方法有最小二乘法(leastsquare，ls)和随机采用一致性算法(randomsampleconsensus，ransac)，前者以整体样本的模型拟合误差平方最小化原则为约束条件确定最佳模型参数，该方法稳定性高，但模型参数求解精度容易受噪声点的影响。后者采用随机选择模型求解最小样本数据集迭代产生最佳模型参数，该方法能很好地削弱噪声点的影响，但受样本点选择方式影响模型拟合稳定性较差。为了提高模型参数求解精度，联合随机采样一致性算法和最小二乘原理求解最佳模型参数即ls-ran法：首先采用随机采样一致性算法识别离群噪声点，最后对剔除噪声点后的有效样本集利用最小二乘法求解模型参数:。

在本发明实施例中，如图6所示，所述步骤s502，包括：

步骤s601，根据随机采样一致性算法识别离群噪声点，获得剔除噪声点后的有效样本集。

步骤s602，根据最小二乘法对所述剔除噪声点后的有效样本集进行所述胸径截面形态模型的模型参数的求解，以获得胸径值提取结果。

在本发明实施例中，图7为试验结果，，本发明采用的ls-ran方法很好的识别出噪声点尽大可能保留有效点，建模点与模型曲线很好的贴合且均匀的分布在模型曲线周围。

本发明实施例提供的基于地基激光雷达技术林木胸径值提取的方法，通过对地基雷达林地点云数据进行点云高程化处理，提取胸径切片点，根据改进k均值聚类对所述胸径切片点进行切割处理，提取单木胸径切片点，对所述单木胸径切片点进行椭圆拟合处理，获得胸径值提取结果。本发明通过改进后的k均值聚类能快速的实现林木胸径点批量化提取，无需林木数目、样地大小等先验知识，具有自动化程度高，抗噪性强的优势，对研究地基激光雷达在林业资源调查及生产管理应用具有一定的实际应用参考价值。

如图8所示，在一个实施例中，提供了一种基于地基激光雷达技术林木胸径值提取的装置，具体可以包括林地点云数据获取单元810、胸径切片点提取单元820、单木胸径切片点提取单元830以及胸径值提取结果获得单元840。

林地点云数据获取单元810，用于获取地基雷达林地点云数据。

胸径切片点提取单元820，用于对所述地基雷达林地点云数据进行点云高程化处理，提取胸径切片点。

在本发明实施例中，林木通过地面连通为一个整体，受地表起伏的影响，不同林木根茎处的高程点不同，即胸径切片点云高程起算面不同，为了方便胸径切片点云批量化自动提取，采用网格高程归一化处理将根茎点高程转换到统一高程面内，具体操作为：将样地划分边长为固定值的正方形格网，索引网格内最低点视为地面点，将网格内的所有点高程与地面点高程作差。则胸径切片点集s1.3可描述为高程1.3m处一定厚度范围内的点云如公式(1)所示，其中pi为样地点云数据m中任意点，ε为胸径处点云切片厚度，通常情况下取胸径尺的宽度3.0～5.0cm。

s1.3＝{pi＝(xi,yi,zi)∈m|zi-1.3|≤ε}(1)

单木胸径切片点提取单元830，用于根据改进k均值聚类对所述胸径切片点进行切割处理，提取单木胸径切片点。

在本发明实施例中，k均值聚类(k-meansclusteringalgorithm)是一种迭代求解的聚类分析方法，具有效率高，操作方便等优点，已广泛应用于市场调查、数据挖掘等领域。其基本思想为：对于一个未标记的数据集，在给定类别数k的前提下，通过随机选择初始种子点即初始目标类簇中心，以对象-类簇中心的欧式距离作为度量指标确定对象类别归属，并把类簇中所有对象的平均值作为新目标类簇的中心，不断迭代直至类簇中心变化范围趋于稳定或达到预期迭代次数，最终产生k个目标类簇子集。k均值聚类的目标是使得类簇内样本间距离尽可能小，不同样本间的距离尽可能大，即k均值聚类算法的目标函数可表示为簇内样本间距之和如公式(2)所示。其中k为目标类簇数，n为样本点总数，xi为属于第j个目标类簇的样本，uj表示第j个目标类簇的中心。当目标函数j取最小值时，算法收敛此时有最佳聚类结果。

原始k均值聚类算法需要目标类别数目k为先验输入，一定程度上降低了算法的自动化程度，同时受随机选择初始种子点的影响，算法容易陷入局部最优产生错误的聚类结果。结合林木胸径切片点云空间分布特性，本发明通过以下三个方面的改进消除k均值聚类在单木胸径提取中的缺陷，以提高林木胸径提取的精度和自动化程度。

(1)优化初始种子点选择原则。根据单木胸径切片点云空间上相互隔离并保持一定距离，且距离远大于胸径值的特性，以林木间距最小值为约束条件选择初始种子点，即如果初始种子点两两距离最小值小于林间距离最小值，则满足条件继续下一步，否则，重新生成初始种子点。通过优化种子点的选择原则，保证初始种子点来源于不同的单木胸径切片点云集，提高了算法聚类精度和收敛速度。

拐点法自适应k值确定。拐点法就是在不同的目标聚类数下计算目标函数j值，然后通过检测变化率的突变确定“拐点”，则突变点对应的k值即为最佳目标类簇数。其原理可解释为随着目标类簇数目的增大，目标类簇集中样本数量逐渐减少，样本与类簇中心的距离更近，通过公式(2)计算则目标函数j值也随之变小。当目标聚类数大于最佳类别数目k时，类簇子集由多个类簇间合并转变为单个类簇分裂形成，此时目标函数j陡然变小并随之保持变缓，因此目标函数的拐点即为最佳聚类数。模拟产生一组试验数据以验证“拐点”法的可行性和有效性，试验数据由三个二维正态分布离散点叠加而成如图3a所示，每个正态分布点集代表一个目标类簇，目标函数j与不同目标聚类数k的对应关系曲线如图3b所示。通过拐点法确定最佳k值为3，其与模拟数据预先设置类别保持一致。

(3)非目标自动识别。由于林地植被生长速率不同，经高程归一化切割提取得到的胸径切片点云中还包含灌木、胸径小于5.0cm的非目标单木，可改进k均值聚类后目标几何中样本与类中心距离分布进行识别：采用反距离加权方法确定簇类中心位置以减弱噪声点的影响，同时计算类簇中样本与中心点的平均距离vd和方差sd，如果vd小于2.5cm或者sd大于方差阈值t，则标记为非目标类簇。公式(3)为反距离加权公式，其中ci为第i类的聚类中心，pj为第j个对象点，sj,smax,smin分别代表第j个点与上次迭代聚类中心的距离、距离最大值和最小值。

改进后的k均值聚类算法能根据样本数据集的特征自适应确定目标聚类数k，同时通过林木间距为约束条件使得初始种子点的选择更为合理，避免了迭代计算陷入局部最优，改进后的k均值聚类胸径切片点云分割步骤如下：

____________________________________________

输入：胸径切片点云，方差阈值t

step1：根据公式(2)计算不同k取值时目标函数值j，通过目标函数变化率寻找拐点，自适应确定最佳目标类别数目k值；

step2：随机选择一组初始种子点并确定种子点间距最小值d，如果d小于林木间距最小值，则进入下一步，否则重新生成种子点；

step3：根据初始种子点，以距离最小化为分类原则，对所有样本点进行归类；

step4：按公式(3)计算类簇内样本点的重心ci，并更新为新的种子点，如果种子点集合不再变化，则进入下一步，否则返回step3；

step5：分别计算类簇几何内样本点与类中心距离的平均距离vd和方差sd，如果vd小于2.5cm或sd大于t，则将该目标类簇标记为非目标集合，否则标记为单木胸径切片点集

输出：单木胸径切片点集和非目标类簇点集

___________________________________________

在本发明实施例中，通过拐点法确定目标聚类簇数提高k均值聚类的自动化程度，采用优化初始种子点选择原则避免算法陷入局部最优，利用预估胸径值剔除非目标单木点集，从而提高单木胸径提取精度。

胸径值提取结果获得单元840，用于对所述单木胸径切片点进行椭圆拟合处理，获得胸径值提取结果。

在本发明实施例中，胸径截面形态可用二维平面圆方程近似表达，其数学表达式如公式(4)所示。其中(a，b)表示胸径中心点坐标，r表示胸径切片半径。则胸径提取值d＝2r。

(x-a)²+(y-b)²＝r²(4)

一般情况下，模型参数求解方法有最小二乘法(leastsquare，ls)和随机采用一致性算法(randomsampleconsensus，ransac)，前者以整体样本的模型拟合误差平方最小化原则为约束条件确定最佳模型参数，该方法稳定性高，但模型参数求解精度容易受噪声点的影响。后者采用随机选择模型求解最小样本数据集迭代产生最佳模型参数，该方法能很好地削弱噪声点的影响，但受样本点选择方式影响模型拟合稳定性较差。为了提高模型参数求解精度，联合随机采样一致性算法和最小二乘原理求解最佳模型参数即ls-ran法：首先采用随机采样一致性算法识别离群噪声点，最后对剔除噪声点后的有效样本集利用最小二乘法求解模型参数。本发明采用的ls-ran方法很好的识别出噪声点尽大可能保留有效点，建模点与模型曲线很好的贴合且均匀的分布在模型曲线周围。

本发明实施例提供的林木胸径值提取装置，通过对地基雷达林地点云数据进行点云高程化处理，提取胸径切片点，根据改进k均值聚类对所述胸径切片点进行切割处理，提取单木胸径切片点，对所述单木胸径切片点进行椭圆拟合处理，获得胸径值提取结果。本发明通过改进后的k均值聚类能快速的实现林木胸径点批量化提取，无需林木数目、样地大小等先验知识，具有自动化程度高，抗噪性强的优势，对研究地基激光雷达在林业资源调查及生产管理应用具有一定的实际应用参考价值。

试验实例：

1、研究区概况与数据获取

研究区位于云南省迪庆藏族自治州香格里拉县县城北部格咱乡境内，地理范围为东经99°29′～100°11′，北纬27°53′～28°37′。格咱乡属高寒山区，平均海拔3200米，最高5090米，最低2520米，地形以山地为主，坡度变化幅度较大。气候垂直变化尤为明显，年平均气温6.6℃，年降水量600毫米。研究区内主要乔木树种主要有云杉(piceaasperata)、红杉(sequoia)、冷杉(fir)、高山松(pinusdensata)、红豆杉(taxus)、云南松(pinusyunnanensis)等，灌木和藤本植物主要有麦宾草(elymustangutorum)、鹅冠草(roegneriakamoji)、三毛草(trisetumbifidum)、苔草(carex)、蒿草(armoise)、野百合(lilyrootflour)等。采用leicascanstationp40地基激光雷达获取研究区内天然近熟林林地点云数据，仪器性能主要参数为：测角精度8″；测距精度1.2mm+10ppm，扫描速率100万点/秒。扫描时仪器架设高度约1.60m，扫描范围设置为水平方向360°，垂直方向270°，扫描模式为精细扫描。为获得完整的林地点云数据，根据树木分布和样地环境均匀布设5个扫描测站如图9a所示，同时在各测站点相互通视范围内，选择不同高程面布设3个黑白标靶，以leicascanstationp40自带cyclone点云数据处理商业软件为数据预处理操作平台，主要包括测站点云数据配准、点云去噪、数据裁剪等，选择30m*30m的林地点云为试验数据如图9b所示。点云数据采集时间为2017年7月中旬，天气晴朗无风。

实测数据主要为样地林木胸径值，根据《gb/t26424-2010森林资源规划设计调查技术规范》规定，采用胸径尺对样地内大于5.0cm的林木胸径逐个测量，测量方式如图10a所示，为提高胸径实测精度取三次测量的平均值为最终胸径值，同时对样地内已测林木进行编号并测量其位置，编号如图10b所示，以方便后续胸径提取精度评定。实测数据采集与点云数据获取保持时间同步。

2、结果分析

以林木胸径值、单木数目等样地实测数据为标准参考，采用单木识别率p、均方根误差mse、误差最大值m、决定系数r²等指标综合评定改进后k均值聚类单木胸径提取精度，各评定指标的定义如下：

(1)单木识别率p(precisionofsingletreedetection)。单木数目识别结果包括2类误差：欠识别误差即实测是单木但未正确识别的漏检误差；过识别误差即实测不是单木而被错误识别为单木的过检误差，则单木识别率可用公式(5)描述，其中n1为欠识别误差，n2为过识别误差，n林木总数。

p＝1-(n1+n2)/n*100％(5)

(2)均方根误差rmse(rootmeansquareerror)。均方根误差描述了胸径值提取结果与实测值之间误差的离散程度，反映了胸径提取方法的稳定性和可靠性，其值越小则算法的适用性和可靠性越好。均方根误差的定义用公式(6)表示，其中di表示第i个单木胸径提取结果；di表示第i个单木胸径的实测值，n表示林木总数。

(3)相对误差m(relativeerror)。误差描述了胸径值提取结果与实测值之间的偏差，相对误差消除测量值纲量对结果的影响，更客观的评价了测量的精度，也是林业调查胸径精度评定的主要指标，相对误差越小则表明精度越高，定义如公式(7)所示，参数含义同上。

m＝|di-di|/di*100％(7)

(4)决定系数r²(coefficientofdetermination)。决定系数表示一个随机变量与多个随机变量关系的数字特征，用来反映因变量与自变量之间相关程度的统计指标，其值越接近于1，则表明相关性越强。可定义为可解释的变异与总变异的比值，其定义如公式(8)所示，其中单木胸径提取的平均值，为实测单木胸径的平均值。

图11展示了胸径提取与分割的结果，其中图11a显示了胸径切片点云提取结果，绿色点表示提取得到的胸径切片点。经高程归一化截取得到的胸径切片点云投影在二维平面内，各单木胸径切片点聚为一簇，不同的单木胸径点彼此分开并保持一定的距离，其距离远大于其林木胸径值。图11b显示了改进k均值聚类算法胸径切片点云分割结果，其中黑色点为识别得到的非目标类簇，即胸径小于5.0cm的单木胸径切片，其余颜色点为正确提取的单木切片。经拐点法确定最佳聚类数目k值为33，经过反距离加权聚类中心计算，结合阈值识别非目标单木胸径类簇1个如图中红色圈出点，该单木胸径4.2cm。正确提取得到有效目标类簇32个，与样地实测数据保持一致，均无过识别和欠识别现象，单木识别率p＝100％。以上分析可知，改进后的k均值聚类方法能很好的识别非目标类簇，实现单木胸径自动化提取。

表1显示了胸径提取及精度统计结果。胸径提取绝对误差最大值为0.5cm，相对误差m最大值为2.00％如表中加粗所示，其中胸径小于20cm的单木提取绝对误差最大值为0.2cm。根据《国家森林资源连续清查技术规定》规定：胸高直径等于或大于20cm的树木，胸径测量误差小于1.5％；胸径小于20cm的树木，胸径测量误差小于0.3cm，本发明方法除9号外胸径提取均在该规定范围之内。在均方根误差精度评定指标方面，均方根误差rmse为0.23cm，表明误差在很小的范围内波动，本发明方法胸径提取方法稳定性好。9号树木的胸径提取结果为异常值，其绝对误差m为0.5cm。分析原因为：9号树胸径处树体损伤引起该处树身凸起，采用ts-ran胸径截面圆方程模型重建时，损伤凸起部分被当做噪声点剔除，导致胸径提取值较实测值偏小如图12所示，

表1胸径提取结果及精度评定

注：表1中胸径实测结果为通过实测周长计算得到。

图13展示了胸径提取值与实测值的回归结果。以胸径实测值为自变量x，点云数据胸径提取值为因变量y绘制散点图并对实测值与提取值进行线性回归分析，通过计算得一元线性回归模型方程为：y＝0.97769x+0.37498，决定系数r²＝0.97355。决定系数非常接近于1，说明基于点云数据胸径提取结果与实测值具有显著的线性相关性，在一定程度上表明改进k均值聚类点云林木胸径提取方法可以代替钢卷尺测量方法开展林木胸径参数提取。

综上，本发明以leicascanstationp40地基激光雷达获取了香格里拉格咱乡境内天然近熟林样地点云数据，通过胸径切片点云提取、切片点云分割和胸径值计算等步骤，实现了林木胸径高精度提取，重点探讨了改进k均值聚类的单木胸径切片点云聚类提取方法，主要结论如下：

(1)以林木间距为约束条件优化k均值聚类算法初始种子点选择原则，避免了种子点的选择不当导致算法陷入局部最优，通过附加样本点与类簇中心平均距离阈值，自动剔除了林木胸径小于5cm非目标胸径提取，单木识别率为100％，一定程度上提高了单木提取精度。

(2)采用拐点法自适应确定样地林木数目，克服了传统k均值聚类中k值先验输入，改善了单木胸径提取的自动化程度，不受林地郁闭度的影响，对大范围林木胸径快速批量化提取具有一定的实际应用参考价值。

(3)通过与胸径实测值对比，胸径提取结果均方根误差rmse为0.23cm，相对误差m最大值为2.00％，决定系数r²为0.97355，表明：结合最小二乘法和随机采样一致性算法的椭圆模型参数求解方法反演林木胸径精度较高，在一定范围内可替代传统人工测量方法开展林木胸径提取。

本发明提出的胸径提取方法具有自动化程度高、操作方便等优点，为森林参数反演提供了新的研究思路和方法，对地基激光雷达技术在林业生产经营管理应用具有重要意义。

在一个实施例中，提出了一种计算机设备，所述计算机设备包括存储器、处理器及存储在所述存储器上并可在所述处理器上运行的计算机程序，所述处理器执行所述计算机程序时实现以下步骤：

获取地基雷达林地点云数据；

对所述地基雷达林地点云数据进行点云高程化处理，提取胸径切片点；

根据改进k均值聚类对所述胸径切片点进行切割处理，提取单木胸径切片点；

对所述单木胸径切片点进行椭圆拟合处理，获得胸径值提取结果。

在一个实施例中，提供一种计算机可读存储介质，计算机可读存储介质上存储有计算机程序，计算机程序被处理器执行时，使得处理器执行以下步骤：

获取地基雷达林地点云数据；

对所述地基雷达林地点云数据进行点云高程化处理，提取胸径切片点；

根据改进k均值聚类对所述胸径切片点进行切割处理，提取单木胸径切片点；

对所述单木胸径切片点进行椭圆拟合处理，获得胸径值提取结果。

应该理解的是，虽然本发明各实施例的流程图中的各个步骤按照箭头的指示依次显示，但是这些步骤并不是必然按照箭头指示的顺序依次执行。除非本文中有明确的说明，这些步骤的执行并没有严格的顺序限制，这些步骤可以以其它的顺序执行。而且，各实施例中的至少一部分步骤可以包括多个子步骤或者多个阶段，这些子步骤或者阶段并不必然是在同一时刻执行完成，而是可以在不同的时刻执行，这些子步骤或者阶段的执行顺序也不必然是依次进行，而是可以与其它步骤或者其它步骤的子步骤或者阶段的至少一部分轮流或者交替地执行。

以上所述仅为本发明的较佳实施例而已，并不用以限制本发明，凡在本发明的精神和原则之内所作的任何修改、等同替换和改进等，均应包含在本发明的保护范围之内。