基于机载激光雷达点云空间分布特征的树木点云提取方法与流程

本发明属于机载激光雷达点云数据处理技术领域，具体涉及一种基于机载激光雷达点云空间分布特征的树木点云提取方法。

背景技术：

激光雷达(lightdetectionandranging，lidar)是近年来发展十分迅速的主动遥感技术，主要通过测定传感器发出的激光在传感器与目标物体之间传播的距离，分析目标地物表面的反射能量大小以及反射波谱的幅度、频率和相位等信息，呈现目标物精确的三维结构信息。可以直接、快速、精确地获取研究对象的三维空间坐标，具有独特的优势。

森林是地球自然资源中的一个重要组成部分，利用遥感技术获取森林冠层结构参数是当前研究热点。机载激光雷达可以获取树木三维结构信息，为了提高运算效率和降低存储量，发展准确高效的算法提取树木点云非常必要。

国内外学者利用机载激光雷达数据已经开展了树木点云提取工作，但大多集中于森林地区。传统林区植被的测量大多通过现场测量或卫星图像进行，其中还面临着现场控制、控制点选择、图像匹配等困难。随着机载激光雷达的出现，由于其对树木等穿透物体的垂直结构具有很强的探测能力，林木三维结构参数准确反演取得了重大突破。

波佩斯库等人基于机载lidar点云数据，通过不断调整窗口大小进行局部过滤处理，有效地提取林区不同层次的植被信息。blair等人利用大点激光雷达技术对热带茂密森林进行扫描和测量(覆盖率为95％)，证明可以准确提取植被信息。美国宇航局开发了激光植被成像传感器(lvis)系统，可用于森林资源的调查和管理，包括计算植被参数和森林垂直结构，如树高、树冠直径、树密度、植被生长、木材数量、树种等。庞勇等人首先对实验区的激光雷达点云数据进行分类处理，获取了地面点、植被点和高程规范化植被点三大类，然后计算出上四分位数的高度基于高程规范化的植被点，并与测得的数据进行比较，最后进行了实际分析。

但是上述方法都采用逐项分离技术，普适性较差，并且在复杂区域内的树木点云提取效果不佳。目前城区的机载雷达点云数据区分植物和其他地物通常采用人工分离或人机交互等方式。人工分离工程量巨大且繁琐，而人机交互的方式则是通过激光雷达首末次回波的强度差与高程等多种方法相结合来进行粗分离，之后再用人力来剔除一些异常点云数据。人机交互的方法虽然比纯人工效率高了不少，但是对于不同的数据仍需要进行针对性的参数调整，这种方法依赖雷达回波强度差值的准确性以及对城市建筑结构的先验知识，也缺乏普适性。

技术实现要素：

针对上述存在的问题或不足，本发明提供了一种基于机载激光雷达点云空间分布特征的树木点云提取方法，即基于原始激光雷达点云数据，通过目标点与周围点云的分布函数运算得到目标点的离散程度值，再通过每个点的离散程度值来进行树木点的判定。

本发明具体技术方案如下：

步骤一、确定目标点与运算点的加权距离影响因子，其中原始点云数据指的是原始数据中的所有点云数据，目标点指的是当前处理的点，运算点指的是原始点云数据中除了目标点本身并与目标点计算的点；

(1)确定目标点与运算点的正投影欧式距离。

设目标点p0正投影坐标为(x0,y0),运算点pi正投影坐标为(xi,yi),则目标点与运算点的正投影欧式距离distanceoi:

(2)确定目标点与运算点之间的生成距离disoi，其中rresolution为最小距离单位长度；rresolution取决于处理数据的扫描分辨率，rresolution∈[0.05,0.001]；

(3)得出加权距离影响因子d0i

(4)得出目标点p0与运算点pi高程差h0i(h0i≥0)；

步骤二、点云运算区块的确定：

(1)由于激光雷达点云数据量大，如果让目标点与数据库中的所有点云都进行一次运算，会极大的增加运算的时间，并且过多的树木点云会影响最后运算的合理性与准确性。因此在运算前划分好每一个目标点的运算区块很有必要。

设原始点云数据的最小x坐标为xmin,最大x坐标为xmax，最小y坐标为ymin,最大y坐标为ymax。

如若xmax-xmin≥100(m),则对原始点云数据进行分割，将点云数据图像以x坐标为变量，平均分成n份，

在待处理的点云数据(即原始点云数据)已经分块完成的基础上，如若ymax-ymin≥100(m),则对原始点云数据进行分割，将原始点云数据以y坐标为变量，平均分成m份，最终原始点云数据共计分为n×m块运算区块。

(2)设定目标点搜索半径rth：rth取决于处理数据的点云密度，rth∈[1,10](m)；

遍历目标点所在的运算区块内的所有点云，当distanceoi小于rth时进行下一步操作，即步骤二(3)，否则跳过该点云，遍历下一点云直到该运算区块内所有点云被遍历完成。

(3)设定目标点高程滤波阈值hth：hth取决于处理数据的点云密度,th∈[5,10]m；

为了减少点云数据中的高程突变点影响，我们只进行hth高程范围内的点云运算，从而将其他地物剔除出去。当h0i小于hth时则开始树木点因子函数值运算，否则跳过该点云，遍历下一点云直到该运算区块内所有点云被遍历完成。

步骤三、计算目标点树木点因子数值v0：

(1)初始化目标点树木点因子数值v0，即令v0＝0。

(2)通过目标点与运算点进行自定义函数运算，通过自定义函数得到的值与d0i成反相关，与h0i成正相关，并累加到目标点树木点因子数值v0上。自定义函数为算法通过局部范围内的起伏变化来判定树木点云，所以从点云的空间分布角度来说，就是目标点周围越靠近的运算点高程起伏程度越大，那么目标点就越有可能是树木点云。因此我们将h0i除以d0i得到一个数值，我们称之为运算点对目标点贡献的树木点因子数值vi，通过目标点与周围的运算点进行不断的运算，从而得到最终的v0，即这样一来就综合考虑到了目标点所在搜索半径内的起伏变化趋势。通过局部范围内的起伏变化趋势，我们可以将点云进行树木点云和非树木点云的归类。

步骤四、提取出树木点云：

设目标点的树木点因子数值阈值为vth，在结束所有目标点的运算之后，通过遍历每一个点云的目标点树木点因子数值v0，比较vth与v0的大小，将v0小于vth的点云删去，保留下来的点云即为树木点云。

步骤二所涉及的原理：

因为目标数据的点云密度分布并不均匀，并且考虑到点云可能出现重复扫描的状况，所以在处理数据的时候会出现诸多点云正投影后过度聚集，从而影响加权距离因子的合理性，所以在计算加权距离因子时并不能直接采用原始的欧式距离数据。如果一个点云中存在着z方向上高度一样但是并不相邻的点时，就会导致那些点的树木点因子数值v0趋向于无穷大，无法和其他地物点进行区别。其次，面对地表起伏较大的数据，输出结果的误差就会被放大。再次，对于建筑物的边缘交界处，存在较大的高度差，不利于建筑物点云与树木点云的区分，在这里我们引入最小距离单位长度rresolution，通过将distanceoi除以rresolution后向上取整得到disoi，这样可以避免在扫描过程中因为重复扫描对结果带来的不良影响。

本发明首先通过对目标点与周围点云关系的确定，划分出目标点运算的范围；再通过目标点与周围点云起伏的变化计算得到树木点云因子数值的累加，从而确定目标点是否是树木点云。其流程如图1所示。

本发明的技术效果：

本发明是通过利用机载激光雷达获取三维激光点云数据，并根据分析树木的云分布特征与其他地物点云分布特征的不同，根据点与点之间的关系，经历点云间关系的确定、运算区域划分、计算树木点云因子数值的处理，建立基于原始激光点云从空间中提取出树木点云的方法。

综上所述，本发明可以摆脱树木点云提取对先验数据的依赖，克服因点云密度分布不均匀所导致的归类错误的问题，以数值计算的方式代替掉传统的逐步分离法，具有首创性和理论上的优越性，提取效果优异。

附图说明

图1为本算法的流程示意图；

图2为样方1编号为test-1的原始点云与树木点云提取对比图(右为原始点云)；

图3为样方2编号为tset-2的原始点云与树木点云提取对比图(右为原始点云)。

具体实施方式

以下通过一个样方的实施例并结合附图对本发明作进一步详细说明：

开发环境为dev-c++5.11，编程语言为c++。

步骤1，以中国内蒙古呼伦贝尔和额尔古纳之间的上库里农场(120°36'至120°52'e，50°21'至52°24'n)为研究区域，并选取实地测量的六个样方作为单木分割验证区，利用leicalas60扫描获取该区域的机载激光雷达数据，航飞时间为2012年9月,具体系统参数参见表1。以下步骤为对样方1编号为test-1的树木点云提取详细步骤。

表1机载激光雷达系统参数

步骤2，在确定了原始机载激光雷达点云数据的文件位置之后，将文件的地址导入到树木点云提取程序中，并且将树木点云提取程序要输出的文件夹地址导入到程序中，并运行树木点云提取程序。

步骤3，树木点云提取程序输出的文件名为“输出数据”，文件格式为txt格式，其中包含的数据即为提取出的树木点云。

步骤4，将提取出的树木点云数据文件“输出数据”导入到lidar360程序中查看，即可获得其三维可视化图形。

样方编号为test-1的输入数据和输出数据的点云三维可视化效果如图2所示。样方编号为test-1共有39589个点云，其中被提取为树木点云的数量为727，被正确提取的数量为648，真实的树木点云数量为648，其卡帕系数为0.94，提取效果具有较高的精度。