本发明属于森林资源调查，森林立地质量评价和森林生产力估测研究等领域，具体涉及一种基于无人机lidar数据的银杏人工林郁闭度提取方法。
背景技术
精确的人工林郁闭度提取，对于森林资源调查，森林立地质量评价和森林生产力估测有重要意义。同时这些信息也可以用于掌握森林冠层空间分布规律，对并森林可持续经营、生态环境恢复与重建及维持碳平衡提供数据支持。传统的人工林郁闭度提取主要依赖于直接法和仪器测量，这些地面实测耗时费力，效率较低，有些直接法往往会对森林冠层具有一定破坏性，且只能获得样地尺度的信息，较难获得大范围的连续郁闭度分布。
无人机激光雷达技术是通过发射激光扫描脉冲到森林表面并分析其返回信号的一种主动遥感技术。通过获取无人机激光雷达点云并提取冠层特征变量，并结合地面实测数据，可以高效率且高精度地获取特定范围的“面”上连续分布的郁闭度。近年来，基于机载激光雷达技术的林分郁闭度反演的研究为：穆喜云等2015年在《东北林业大学学报》第9卷上发表的“基于机载lidar数据的林分平均高及郁闭度反演”，该研究使用机载激光雷达数据，提取高度百分位数特征变量，结合地面实测有郁闭度数据，构建线性回归模型，估测我国内蒙古大兴安岭地区林分郁闭度。moeser等2014年在《agricultural&forestmeteorology》第197卷上发表的“laiandradiationtransferfromairbornelidarsyntheticimages”，该研究使用记载激光雷达合成的图像数据结合半球影像提取的郁闭度去估测瑞士达沃斯地区云杉森林的郁闭度。然而以上方法都使用低密度的有人机激光雷达数据运用单一的特征变量去估测郁闭度，未见高密度的无人机激光雷达点云数据应用于郁闭度反演，同时，更未见全面深入计算无人机激光雷达点云特征及郁闭度提取的方法。
技术实现要素：
发明目的：针对现有技术存在的不足，本发明的目的是提供一种基于无人机lidar数据的银杏人工林郁闭度提取方法。
技术方案：为了实现上述发明目的，本发明采用的技术方案如下：
一种基于无人机lidar数据的银杏人工林郁闭度提取方法，包括对激光雷达离散点云数据滤波，插值生成数字地形模型，使用生成的数字地形模型对激光雷达点云数据归一化处理；然后分别基于归一化点云数据提取冠层结构特征变量；最后，结合地面实测郁闭度和提取的特征变量，在优选变量的基础上进行逐步回归构建预测模型对郁闭度进行估测。
一种基于无人机lidar数据的银杏人工林郁闭度提取方法具体步骤如下：
1)在地面设置样地，借助多旋翼无人机搭载的激光雷达传感器进行激光雷达数据采集；lidar即激光雷达测距传感器(lightdetectionandranging-lidar)；并在样地中记录树种并计数，同时测量每木的胸径和树高，采用点云阈值距离分割法底地面实测郁闭度；
2)对遥感数据进行处理：首先去除无人机激光雷达原始点云数据的噪音点，采用滤波算法去除非地面点，然后通过计算每个像元内激光点高度的平均值，生成数字地形模型，数字地形模型的空间分辨率为0.5m，通过生成的数字地形模型对点云数据进行归一化处理，得到归一化后的无人机激光雷达点云数据；
3)基于归一化点云数据提取冠层结构特征变量；
4)通过相关性分析筛选特征变量，首先筛选特征变量之间相关性低于0.6的特征变量，然后进一步筛选特征变量与郁闭度相关性高于0.6的特征变量；
5)将地面实测郁闭度作为因变量，无人机激光雷达点云特征变量作为自变量，建立多元回归模型。
进一步地，步骤1)中的地面实测郁闭度：首先假定归一化点云数据中最高点为树顶点，然后从该点进行区域增长，分割出一棵树，最后通过迭代算法，直到林分中的单木都被分割出来，对分割出来的样地内单木进行冠幅叠加，计算样地内冠幅总面积(对于重叠部分不累加计算)m，将计算的冠幅总面积m除以样地面积得到地面实测森林郁闭度cc(canopyclosure)，计算公式如下：
cc＝m/(15×15)×π(1)
进一步地，步骤3)中，冠层结构特征变量包括百分位高度变量、各层覆盖度和冠层体积与剖面特征变量。
进一步地，所述百分位高度变量冠层高度分布百分位数(h25，h50，h75，h95)，冠层点云分布平均高度以上的覆盖度(ccmean)，冠层点云分布的变异系数(hcv)。
进一步地，所述各层覆盖度为点云数量在各百分数高度(30th，50th，70th，90th，即d3，d5，d7，d9)以上的点占所有点云的百分比。
进一步地，冠层体积与剖面特征变量：剖面特征变量为weibull函数对冠层高度分布剖面进行拟合得到2个剖面特征量α，β(即weibullα和weibullβ)；冠层各结构类别体积占比为开放层，透光层，低光层和封闭层四个冠层结构类别，每个冠层结构类别的体积所占百分比(即opengap，oligophotic，euphotic，closedgap)。
进一步地，步骤5)中，运用逐步进入法将与实测森林郁闭度的激光雷达特征变量按照相关性由高到低依次引入方程，筛选出相关性高的变量，剔除相关性低的变量，最终得到最优模型。
进一步地，步骤5)中，通过主成分分析计算相关关系矩阵得到控制因子，即最大特征根的平方根与最小特征根的比值k，k小于30则模型进一步被选中。
进一步地，步骤5)中，采用决定系数(r2)、均方根误差(rmse)和相对均方根误差(rrmse)评价回归模型拟合的效果及估测精度：
式中xi为某林分郁闭度实测值；
有益效果：与现有技术相比，本发明具有以下优点：
1)以往的方法都是基于有人机激光雷达数据，数据点云密度低，而本发明则是采用无人机激光雷达数据，然后基于高密度的点云数据获取冠层结构特征；
2)本发明全面深入地提取了多组人工林冠层无人机激光雷达点云特征，并进行了特征变量优选，从而高质量地提取了人工林林分郁闭度。同时，该发明不仅利于特征变量的机理解释，也易于进行方法移植(即在天然林和次生林中也可以进行应用)。
本发明通过获取无人机激光雷达点云并提取冠层特征变量，并结合地面实测数据，可以高效率且高精度地获取是特定范围的“面”上连续分布的郁闭度。而本发明则是采用无人机激光雷达数据，然后基于高密度的无人机激光雷达点云，全面深入地提取了多组人工林冠层无人机激光雷达点云特征，从而高质量地提取了人工林林分郁闭度。验证结果表明，通过本发明对人工林郁闭度进行提取，与使用其他相近遥感方法进行林分郁闭度反演相比，其相对均方根误差降低1.7％。
附图说明
图1是本发明样地单木冠幅提取郁闭度方法图，a为样地正射影像及单木冠幅；b为样地范围内单木提取的冠幅；c为样地范围内冠幅总面积；
图2是基于高度特征变量的模型交叉验证结果；
图3是b为基于高度特征变量与覆盖度特征变量组合的模型交叉验证结果；
图4是c为基于高度特征变量与冠层容积变量组合的模型交叉验证结果。
具体实施方式
下面结合具体实施例，进一步阐明本发明，实施例在以本发明技术方案为前提下进行实施，应理解这些实施例仅用于说明本发明而不用于限制本发明的范围。
近年来，机载激光雷达技术逐步应用于林分郁闭度的获取中。无人机激光雷达技术是通过发射激光扫描脉冲到森林表面并分析其返回信号的一种主动遥感技术。通过获取无人机激光雷达点云并提取冠层特征变量，并结合地面实测数据，可以高效率且高精度地获取特定范围的“面”上连续分布的郁闭度。本发明对无人机激光雷达离散点云数据滤波并插值生成数字地形模型，使用生成的数字地形模型对无人机激光雷达点云数据进行归一化处理从而获得归一化点云；分别基于归一化无人机激光雷达点云提取多组人工林冠层结构特征变量，在优选特征变量的基础上，结合地面实测数据进行逐步回归预测模型反演人工林郁闭度。
通过以下实例对本发明作进一步描述：
研究区位于江苏省北部地区邳州市铁富镇，其地理位置北纬34°33′49″-34°34′23″，东经118°05′1″-118°06′06″。属于半湿润温带季风气候，年降雨量约为903mm，最大降雨量集中在7、8月份梅雨季节，年平均温度约为13.9℃，无霜期211天，主要土壤类型为黑粘土壤，呈酸性。本地区银杏大规模种植始于1993年，总面积约5400hm2。
详细过程如下：
(1)借助多旋翼无人机搭载的激光雷达传感器进行激光雷达数据采集。根据历史森林资源调查数据和前期获取的卫星遥感影像数据在银杏人工林核心分布区内选取了5块1×1km的正方形大样地，然后在5个样地内按照典型选样的方法设置9块半径为15m的圆形小样地，样地的中心位置通过trimblegeoxh6000handhelds手持gps定位(结合jscros广域差分系统)进行定位，精度优于0.5m。并在样地中记录树种并计数，同时测量每木的胸径和树高，地面实测郁闭度采用点云阈值距离分割法所分割的单木进行单木冠幅提取：首先假定了归一化点云数据中最高点为树顶点，然后从该点进行区域增长，分割出一棵树，最后通过迭代算法，直到林分中的单木都被分割出来。对样地内单木进行冠幅叠加，计算样地内冠幅总面积(对于重叠部分不累加计算)m，将计算的冠幅总面积m除以样地面积得到地面实测森林郁闭度cc(canopyclosure)，计算公式如下：
cc＝m/(15×15)×π(1)
(2)遥感数据处理：首先去除无人机激光雷达原始点云数据的噪音点，基于滤波算法去除非地面点，然后通过计算每个像元内激光点高度的平均值，生成数字地形模型(空间分辨率为0.5m)。通过生成的数字地形模型对点云进行归一化处理，得到归一化后的无人机激光雷达点云数据。
(3)冠层结构特征变量提取。提取三组特征变量，即百分位高度、各层覆盖度和冠层体积与剖面特征变量。
)1)百分位高度变量包括：冠层高度分布百分位数(h25，h50，h75，h95)，冠层点云分布平均高度以上的覆盖度(ccmean)，冠层点云分布的变异系数(hcv)；2)各层覆盖度：点云数量在各百分数高度(30th，50th，70th，90th，即d3，d5，d7，d9)以上的点占所有点云的百分比；3)冠层体积与剖面特征变量：采用weibull函数对冠层高度分布剖面进行拟合得到2个剖面特征量α，β(即weibullα和weibullβ)；冠层各结构类别体积占比，包括开放层，透光层，低光层和封闭层四个冠层结构类别，每个冠层结构类别的体积所占百分比(即opengap，oligophotic，euphotic，closedgap)。
(4)通过相关性分析筛选特征变量。即首先筛选特征变量之间相关性低于0.6的特征变量，然后进一步筛选特征变量与郁闭度相关性高于0.6的特征变量。
(5)将地面实测郁闭度作为因变量，无人机激光雷达点云特征变量作为自变量，建立多元回归模型。运用逐步进入法将与实测森林郁闭度的激光雷达特征变量按照相关性由高到低依次引入方程，筛选出相关性高的变量，剔除相关性低的变量，最终得到最优模型。在回归建模过程中，每一步引入变量检验值f都要达到显著水平(p值<0.05)。为了降低自变量之间的相关性，本发明通过主成分分析计算相关关系矩阵得到控制因子k(即最大特征根的平方根与最小特征根的比值)，k小于30则模型进一步被选中。
本发明采用决定系数(r2)、均方根误差(rmse)和相对均方根误差(rrmse)评价回归模型拟合的效果及估测精度：
式中xi为某林分郁闭度实测值；
林分郁闭度模型及模型预测精度见表1，基于无人机lidar特征变量反演郁闭度与样地实测郁闭度交叉验证结果图见图2。
表1不同特征变量组合作为预测变量的森林郁闭度估测模型比较及精度评价
注：h50，h75，为冠层50％，75％的高度分布百分位数；d5，d7为点云数量在百分数高度50th和70th上的点占所有点云的百分比；α，β为weibull函数对冠层高度分布剖面进行拟合得到的2个剖面特征量；open，euphotic冠层开放层各和低光层结构类别体积占比。r2为决定系数；rmse为均方根误差；rrmse为相对均方根误差。