一种基于UAV-LiDAR的杉木人工林地精细地位指数精准评价方法

本发明属于森林立地条件的评价方法领域，具体涉及一种基于uav-lidar的杉木人工林地精细地位指数精准评价方法。

背景技术：

1、森林立地质量不仅影响森林功能的发挥，还影响生态环境和森林资源的可持续发展。森林立地质量精准评价是是人工林精准培育的重要依据，是实现人工林适地适树原则的基础，为有效提升人工林培育质量和林地生产力提供重要信息支撑。评价立地质量常用方法中，地位指数较地位级更为直观，且优势木平均高受营林措施和密度的影响小，比较稳定，是评价立地的主流方法。目标林分的地位指数通常是由经验丰富的林业技术人员经过全林分实地踏查后，目估选择林分优势木，计算林分优势木平均高，结合林分年龄，在编制的地位指数表中查找对应的地位指数得到。由于各自然因子之间的相互作用，复杂变化的地形因子不仅直接影响立地质量，还通过影响其他因子间接影响立地质量。所以对于地形较平坦，空间异质性较小，立地质量差别不大的平原森林，按传统方法调查的地位指数对整个林地立地质量的代表性相对较高。但是对于复杂山地上而言，通过传统方法的得到的地位指数的精度与其对整个林分的代表性较为有限。且传统方式较为依赖调查人员的经验，调查人员的主观性亦会对标准地代表精度产生较大影响。继续沿用传统方式调查地位指数难以满足精确林业发展的需要。

2、上述限制主要是缺乏一种高效、系统的可以准确获取目标林分精细地位指数的立地质量精准评价技术方法。随着科技的发展，高效率、高精度的无人机遥感技术被广泛运用于林业经营与生产，uav-lidar可以摈除复杂山地的多种干扰因素，快速提取单木尺度的林分生长参数，并进行高精度定位。

3、为解决当前森林立地条件的评价工作中存在的地位指数调查精度有限的问题，本发明充分考虑林分内部的空间异质性对立地质量存在的影响，通过局部莫兰指数量化林分空间异质性，并针对性选择不同的分割网格的尺寸，提取网格尺度的精细地位指数，结合机载lidar的数据获取手段和计算机辅助决策，发明了一种更为高效、细致、准确的杉木人工林地精细地位指数评价方法。

技术实现思路

1、本发明的目的在于提高当前森林立地条件的评价工作中地位指数调查的精度，提供一种基于uav-lidar的杉木人工林地精细地位指数精准评价方法。本发明主要围绕林地立地质量的空间异质性，进行等级量化，选用适当的尺寸，参考林分微地形，进行林分网格化，细化林分地位指数，优化林分地位指数提取的精度。

2、本发明所采用的具体技术方案如下：

3、本发明提供了一种基于uav-lidar的杉木人工林地精细地位指数精准评价方法，具体如下：

4、s1、应用uav-lidar技术对目标杉木人工林进行lidar点云数据采集，并对采集的lidar点云数据进行预处理；

5、s2、对s1中预处理后的lidar点云数据进行分离地面点与非地面点操作；对分离的地面点进行模型构建，得到目标林分的数字高程模型；根据分离出的非地面点，构建出目标林分的数字表面模型；根据目标林分的所述数字高程模型与所述数字表面模型同一坐标点位的高程差值，获得冠层高度模型；根据s1中预处理后的lidar点云数据和所述数字高程模型获得归一化lidar点云数据；

6、s3、对s2中所述数字高程模型进行正地形、负地形、反地形的分析和制作，提取目标林地的山脊与山谷并将其矢量化；

7、s4、对s2中所述冠层高度模型进行基于标记分水岭控制法与局部最大值法的单木分割并提取单木树高，获得目标林分的单木空间经纬度投影坐标与树高；

8、s5、对s4中所述单木空间经纬度投影坐标进行矢量化，建立目标林分的单木空间位置与树高数据库；采用局部莫兰指数法分析目标林分单木树高的空间集聚程度，对目标林分的空间异质性进行分级评价；

9、s6、基于s5获得的空间异质性评级结果，根据目标林分的空间需求情况，兼顾微地形对林木的影响，分别设置目标林分不同精细程度的网格尺寸；

10、s7、将s3中矢量化后的山谷与山脊作为微地形边界区，叠加s6中相应尺寸的精细化网格，参考s2中获得的归一化lidar点云数据，进行目标林分范围山谷与山脊的网格化，并根据山谷与山脊矢量化出山坡区，分别得到并标记目标林分的网格化山谷区、山脊区与山坡区，获得目标林分矢量化网格；

11、s8、根据s5中所述单木空间位置与树高数据库以及s7中所述目标林分矢量化网格，建立各个网格内的单木树高和数量数据库，进一步统计筛选出各网格内的优势木平均高sht；

12、s9、根据s8中各网格优势木平均高sht，对照所在省市地区的杉木人工林地位指数表，得到目标林分的网格尺度精细地位指数；再结合目标林分的遥感影像，进行网格尺度精细地位指数可视化制图。

13、作为优选，所述s1中，lidar点云数据的采集是通过大疆经纬matrice 300rtk无人机搭载激光雷达传感器在目标林分执飞作业实现的。

14、作为优选，所述s1中，lidar点云数据的预处理方法包括三维重建及去噪，具体如下：

15、采用进阶版dji terra软件进行lidar点云数据的三维重建；采用cloudcompare软件进行非目标地物的手动去噪；采用lidar 360软件进行飘移噪声点的去噪。

16、作为优选，所述s2中，分离lidar点云数据的地面点与非地面点采用布料模拟滤波算法，具体如下：

17、(1)倒置s1中经预处理后的lidar点云数据，模拟一块足够柔软的布粘铺于倒置的lidar点云数据上方，将组成模拟布的最小纤维分子看为无数个布料粒子，限定布料粒子仅沿垂直方向运动，则布料粒子的位置和速度由布料粒子间的内力和外力所决定，通过比较布料粒子和地形的高度值来分类lidar点云数据中的地面点；

18、(2)设置格网的规格大小，网格化步骤(1)中模拟的lidar点云布，以格网为基本单元进行点云滤波；

19、(3)基于步骤(2)的结果，将目标区域的所有lidar点云与布料粒子投影至相同的二维水平面，其中lidar点云投影在平面上的粒子称为点云粒子；根据点云粒子与布料粒子之间的平面二维坐标距离，找到该二维平面中每个布料粒子对应的lidar点，定义该lidar点的高度为ih，ih为布料粒子能下落的最低高度；

20、(4)基于步骤(3)的结果，计算布料粒子在外力作用下所发生的位移，因粒子在空中未受碰撞力作用，仅受重力影响，根据如下公式计算布料粒子在重力作用下的移动距离：

21、

22、式中，t代表某一时间节点；δt代表时间步长；p(t)代表布料粒子仅受重力作用下在t时间节点所在的位置；p(t+δt)代表布料粒子仅受重力作用下在δt时间后所在的位置；p(t-δt)代表布料粒子仅受重力作用下在该节点的δt时间前所在的位置；m代表布料粒子的质量，g为常数；

23、(5)基于步骤(3)的结果，计算布料粒子在内力影响下产生的位移距离，公式如下：

24、

25、式中，表示布料粒子的位移向量；b为常数，当布料粒子可移动时值为1，不可移动时值为0；代表布料粒子的当前坐标向量；代表与布料粒子相连的相邻粒子的坐标向量；代表垂直方向的单位向量；

26、(6)根据步骤(4)和(5)的结果，将受外力和内力影响后发生位移的布料粒子的高度ch与ih进行比较；若ch≤ih，则将其放置于ih的高度位置，并设为不可移动粒子；若ch>ih，则作为可移动粒子，重复步骤(4)～(5)，迭代至所有布料粒子的最大变化高度足够小或不可移动或超过设定的最大迭代次数为止，此时所有布料粒子组成的平面便为近似的真实地形，该布料模拟迭代完成；

27、(7)根据步骤(6)的结果，计算布料模拟完成后的布料粒子和对应lidar点的距离mh；

28、(8)根据步骤(7)的结果，设定距离阈值oh；比较每个lidar点的mh与oh大小，当mh≤oh时lidar点标记为地面点，当mh＞oh时lidar点标记为非地面点。

29、作为优选，所述s2中，利用克里金插值法构建数字高程模型和数字表面模型，公式如下：

30、

31、式中，n为可以决定该点属性值的点数量；z′0为待估点z0的估计值，zi代表第i个已知点的属性值；wi代表第i个已知点的权重；

32、

33、式中，γij表示第i，j个领域点之间的半方差；zn代表第n个已知点的属性值；

34、选择合适的函数形式，拟合函数γ＝f(d)；待求点位置已知，可得d10，…，dn0，带入函数得到γ10，…，γn0，计算矩阵乘式：

35、

36、式中，λ是常数；第0个领域点表示待求点本身；

37、最后完成待估点z0的估值。

38、作为优选，所述s3具体如下：

39、(1)利用arcgis软件，以s2所述数字高程模型为基础数据，通过地形领域分析中的焦点统计方法获得统计栅格数据focal_st；计算数字高程模型与focal_st在同一个空间位置的差值并获得栅格数据r1_ca；以0为临界值重分类r1_ca，将0以上的r1_ca作为正地形栅格reclassify1，0以下的r1_ca作为负地形栅格reclassify2；

40、(2)以s2所述数字高程模型为基础数据，进行栅格计算并获得反地形栅格数据reclassify3；

41、(31)以s2所述数字高程模型为基础数据，进行该地形水文的填洼、流向、流量分析及栅格计算，获得流量为0的栅格区域数据r2_co；通过对数字高程模型数据进行表面等高线和山体阴影的分析，确定r2_co的山脊分类阈值；根据分类阈值重分类r2_co并获得reclassify4；

42、(32)将reclassify1与reclassify4相乘，获得栅格数据reclassify；重分类保留值为1的栅格即为山脊栅格区，矢量化该结果即可提取山脊；

43、(41)以reclassify3为基础数据，进行该地形水文的填洼、流向、流量分析及栅格计算，获得流量为0的栅格区域数据r3_co；通过对数字高程模型数据进行表面等高线和山体阴影的分析，确定r3_co的山谷分类阈值；根据分类阈值重分类r3_co并获得reclassify5；

44、(42)将reclassify2与reclassify5相乘，获得栅格数据reclassify6；重分类保留值为1的栅格即为山谷栅格区，矢量化该结果即可提取山谷；

45、(5)将reclassify与reclassify6相乘，获得栅格数据reclassify7；重分类保留值为0的栅格即为山坡栅格区，矢量化该结果即可提取山坡。

46、作为优选，所述s4具体如下：

47、(1)根据单株树冠的高度变化特征，使用局部最大值算法识别单木树梢点；识别的结果被视为树木的x、y坐标，并用其标记前景以优化分水岭分割的效果；

48、(2)用基于冠层高度模型的局部最大值探测法实现单木识别，在迭代过程中动态调整窗口大小以确定最优值，并用具有最优窗口大小的滑动窗口搜索局部最大值，最终输出识别到的单木坐标，并标记在dom上；

49、(3)对冠层高度模型的像素值取反，将局部最大值算法识别的树梢点标记为前景，阈值分割确定的非林冠区域标记为背景，其余不能确定部分标记为0；接着使用分水岭算法，每次灌水，标签会被更新，当两种不同颜色的标签相遇就会构建分水堤坝，最终形成封闭的轮廓；每个局部最小值及其影响范围形成积水盆地，积水盆地的边界便是分水岭；

50、(4)取输入冠层高度模型的投影和仿射变换，给图像添加投影坐标，并保存为tif格式的栅格影像；按照树id将栅格影像的栅格数据转换为矢量面数据，转换后的数据中会存在一些小面积的图斑，将这些小图斑合并到距离最近的大图斑中。

51、作为优选，所述s5具体如下：

52、利用arcgis软件对获得的目标林分单木树高值进行局部莫兰指数分析，得到目标林分单木树高的高值与低值集聚模式，根据非高-高集聚模式与低-低集聚模式的无聚类单木点的数量占全部单木点数量的比例，对目标林分的立地质量空间异质性按五级进行评级。

53、进一步的，所述s6中，自定义与s5中五种评级结果对应的五种不同精细程度的网格尺寸。

54、作为优选，所述s8中，对各网格内的单木进行树高大小排列，选择前20％作为目标树种的优势木，计算该网格的优势木平均高sht。

55、本发明相对于现有技术而言，具有以下有益效果：

56、1)全面科学。基于uav-lidar数据以科学的数据全面呈现出目标林分的全部单木的位置与树高参数，弥补了传统人工踏查法仅能根据视域范围内的林分生长情况得到经验值的不足。

57、2)量化标准。基于全林分单木树高的局部莫兰指数，将目标林分的立地质量空间异质性分级量化为五级，以此为依据选择网格尺寸，解决了立地质量空间异质性难量化，地位指数的细化程度难确定的问题。

58、3)高效节约。采用了机载lidar的数据获取手段和计算机辅助决策，有效缩短了外业调查的时间，节约了人力、物力成本。