CN104867180A

CN104867180A - 一种集成uav与lidar的森林林分特征反演方法

Info

Publication number: CN104867180A
Application number: CN201510284456.XA
Authority: CN
Inventors: 许子乾; 曹林; 阮宏华; 李卫正; 蒋圣
Original assignee: Nanjing Forestry University
Current assignee: Nanjing Forestry University
Priority date: 2015-05-28
Filing date: 2015-05-28
Publication date: 2015-08-26
Anticipated expiration: 2035-05-28
Also published as: CN104867180B

Abstract

本发明公开了一种集成UAV与LIDAR的森林林分特征反演方法，通过结合UAV与LiDAR数据提取典型林分特征，利用特征点匹配算法和空中三角生成研究区三维点云作为数字表面模型结合LiDAR提取的数字地形模型计算树冠冠层模型，再结合地面实测林分调查数据与联合提取的点云变量构建多元回归估算模型并验证其精度。一方面冠层点云采用UAV数据，相对于LiDAR，具有灵活性高，廉价等特点，数据获取方便，具有长期定时重复观测的可能性。另一方面利用LiDAR的穿透性获取地形数据，而地形相对稳定，可以作为长期观测的数据基础。二者结合，各取所长，是一条廉价、灵活的森林状态监测方式，相关分析表明联合提取的点云变量与森林垂直信息（树高）具有很高的敏感性。

Description

一种集成UAV与LIDAR的森林林分特征反演方法

技术领域

本发明涉及森林资源管理与保护技术领域，具体涉及一种集成UAV与LIDAR的森林林分特征反演方法。

背景技术

林分平均高、平均胸径、林分密度、蓄积量等林分特征变量的精确测定对于森林管理与生态监测具有重要意义。传统的测定方式为人工实地测量，耗时耗力且只能获取小面积数据，难于获取区域尺度或更大尺度的信息。遥感方法能快速的获取大面积森林状态，降低人力和物力的消耗，且在数据的完整性和时间性上具有较高的精度。但光学遥感数据主要记录的是森林的水平结构信息，其遥感信号穿透性较低，难以有效反映森林的垂直结构；且在茂密林地区域，就算是森林具有很高的复杂度和变异度，表现在遥感图像上的光谱差异也很小，导致林分特征信息呈现区域性光谱饱和。

激光雷达(LiDAR)是近年来发展迅速的主动遥感技术，激光雷达通过激光直接照射，返回三维点云和波形数据，可用于森林的生物物理特性和冠层结构的精确估算。但激光雷达数据成本高，较难在林业生产实践中广泛应用。近年来，无人机(UAV)以其在云下低空飞行的能力，以及其低成本和灵活性，弥补了卫星光学遥感和普通航空摄影易受云层遮挡影响的同时，也降低了地物三维信息(如点云)的获取成本。

发明内容

发明目的：针对现有技术中存在的不足，本发明提出一种集成UAV与LIDAR的森林林分特征反演方法，具有廉价、灵活等特点。

技术方案：为了实现上述发明目的，本发明采用的技术方案为：

一种集成UAV与LIDAR的森林林分特征反演方法，包括以下步骤：

1)LiDAR数据获取：LiDAR数据获取装置为Optech Airborne Laser TerrainMapper 3100激光传感器斑离散点云数据；每个激光返回点包含x,y,z三维坐标值、强度值和返回类型信息；

2)UAV数据获取：采用的无人机搭载平台为固定翼飞机，配备惯性测量单元获取飞机实时的三轴姿态角\角速率以及加速度；无人机通过地面遥控系统进行飞行的控制和通讯；遥感数据来源于佳能EOS 5D Mark II高清相机获取的高分辨率、真彩色数码影像；

3)地面实测数据获取：设计地面实测采用圆形样地方式，多个随机分布于整个待研究区；在各个样地中，记录胸径≥5cm的每棵树木的树种，并量测胸径、树高；胸径采用围尺量测，树高通过Vertex IV超声波测高器获取；在样地内进行每木测量，选取样地Lorey’s树高、林分密度和蓄积量为反演的林分特征，计算方法如下：

Lorey’s树高计算公式如下：

H = \frac{Σ_{i - 1}^{k} h_{i} g_{i}}{Σ_{i - 1}^{k} g_{i}}

式中：h_i为第i株树高，g_i为第i株树木胸高断面积，k为样地内株数；

林分密度通过为单位面积上的树木数，以株·hm^-2为单位，通过样地实测林分密度转换每公顷林分密度得到；

蓄积量依据“江苏省立木一元材积式”计算，通过单木胸径汇总到样地级蓄积量，再换算为单位面积蓄积量；

4)利用LiDAR数据生成DEM；

5)通过UAV平台拍摄数字影像，通过IMU实时记录每幅影像的经纬度、高程、翻滚度、俯仰和航向五个参数生成Exif图像文件；

6)借助Pix4D软件获取影像像对的特征点，通过光束法平差进行影像外方位元素的解析，完成图像标定；

7)利用已有三维点、影像定标和外方位元素解析生成三维点云，在剔除明显的异常点后，通过平滑处理实现对“孔洞”的填充，再利用反距离加权法插值为与Lidar数据相适配的1m×1m的栅格生成DSM以便点云的归一化；同时通过影像拼接生成研究区UAV全幅影像，作为进行Lidar点云与UAV点云精度匹配的辅助验证底图；

8)点云的归一化指对点云高度处理，即距离地表的相对高度的获得，方式由DSM与DEM差值而来，归一化的结果为数字树冠高层模型；

9)变量提取：将DSM与DEM进行差值运算，获得DCHM，完成点云数据的归一化；建模数据变量从归一化后的点云中进行提取，选取高度分位数变量组：h10,h25,h30,h40,h60,h75,h85,h90；密度变量组：d10,d25,d30,d40,d60,d75,d85,d90；平均高度，最大高度，最低高度；同时去除了高度2m以下点云；

10)通过获取归一化点云提取的林分特征变量与实测样地4个关键变量之间线性关系进行精度的验证；验证的方式采取多元线性回归方法；建立回归模型运用逐步进入法和决定系数R²的变化情况来选择进入模型的合适变量，所有样地均参与模型构建；

11)精度验证

验证的指标选择：(1)决定系数R²：

R^{2} = 1 - Σ_{i = 1}^{n} {(x_{i} - {\hat{x}}_{i})}^{2} / Σ_{i = 1}^{n} {(x_{i} - {\overset{&OverBar;}{x}}_{i})}^{2}

式中，x_i为实测的某林分特征、为模型估算值、为实测特征值算术平均、n为样地数量；

(2)均方根误差：表达的是观测值与真值偏差的平方和观测次数n比值的平方根；

RMSE = \sqrt{\frac{1}{n} Σ_{i = 1}^{n} {(x_{i} - {\hat{x}}_{i})}^{2}}

(3)相对均方根误差：定义为RMSE与实测结果算术平均值的比值；

rRMSE = \frac{RMSE}{{\overset{&OverBar;}{x}}_{i}} \times 100 % - - - (4)

12)联合提取估算模型构建及精度分析

通过多元回归分析，构建联合提取估算模型，并进行精度分析，结果如下：

林分特征变量为Lorey’s树高/(m)，联合提取估算模型为ln H＝0.36+0.651·lnhmean+0.655·ln d90-5.1215·ln d30，R²为0.95，RMSE为0.06，rRMSE为1.9％；

林分特征变量为林分密度/(株·hm^-2)，联合提取估算模型为ln N＝5.607+1.854·lnh60-1.195·ln h10，R²为0.36，RMSE为0.29，rRMSE为13.74％；

林分特征变量为胸高断面积/(m2·hm^-2)，联合提取估算模型为ln G＝2.752·lnh60-1.841·ln h10-1.126，R²为0.53，RMSE为0.28，rRMSE为16.38％；

林分特征变量为蓄积量/(m3·hm^-2)，联合提取估算模型为ln V＝2.499+1.429·lnH90+0.7·ln d90，R²为0.59，RMSE为0.40，rRMSE为6.93％；

其中，H为Lorey’s树高；N为林分密度；G为胸高断面积；V为蓄积量。

步骤2)中，航向和旁向重叠分别为30％和70％。

步骤2)中，IMU型号为ADIS16355，集成3轴加速度计及3轴陀螺，技术参数为：1)可控最高空速200m·s^–1；2)可控最高地速350m·s^–1；3)可控飞行高度20–5500m；4)导航精度(偏航距)≤±3m；5)定高控制精度≤±2m。

步骤4)中，数据处理的过程包括以下三步：(1)数据预处理，包括异常点剔除和重叠数据处理，首先去除高程噪声点，噪声点分两类，一类为高程负值点，一类为高程突变剧烈点；其次采用以大于激光点云平均间距的平滑滤波进行多次平滑消除重叠带；(2)点云数据分类，将点云数据分为裸地与非裸地两类，激光点云的分类采用滤波和分类的方法提取裸地高程信息，将生成的裸地点云以点状矢量文件输出，再与UAV影像叠加进行目视纠正；(3)DEM生成，将分类后的裸地点采用反距离加权法进行1m×1m的插值，生成DEM。

有益效果：与现有技术相比，本发明的特色和创新点如下：本发明通过结合UAV与LiDAR数据提取典型林分特征，利用特征点匹配算法和空中三角生成研究区三维点云作为数字表面模型结合LiDAR提取的数字地形模型计算树冠冠层模型，再结合地面实测林分调查数据与联合提取的点云变量构建多元回归估算模型并验证其精度。一方面冠层点云采用UAV数据，相对于LiDAR，具有灵活性高，廉价等特点，数据获取方便，具有长期定时重复观测的可能性。另一方面利用LiDAR的穿透性获取地形数据，而地形相对稳定，可以作为长期观测的数据基础。二者结合，各取所长，是一条廉价、灵活的森林状态监测方式，相关分析表明联合提取的点云变量与森林垂直信息(树高)具有很高的敏感性。

附图说明

图1是反演方法的流程图；

图2是牛山无人机影像及样地空间分布图；

图3是UAV与LiDAR联合提取特征变量与实测林分特征之间的相关分析结果图；

图4是各林分特征值与模型估算值对比图；图中，虚线为1:1验证线。

具体实施方式

下面结合具体实施例对本发明作进一步的说明。

实施例1

一种集成UAV与LIDAR的森林林分特征反演方法，基于LiDAR点云获取的林下高精度数字高程模型(DEM)和UAV像对序列构建的植被冠层上层三维点云(DSM)，提取高度和冠层点云密度相关的特征变量并构建林分特征变量预测模型，技术流程如图1所示，具体过程如下：

研究区牛山位于江苏省南京市东郊紫金山余脉，中心地理坐标(118°48'E，32°04'N)，总面积13.90hm²，高程范围16～40m。为北亚热带季风气候，四季分明，年均降雨量1005.9mm。主要森林类型为北亚热带常绿落叶阔叶混交林，主要树种为榉树(Zelkova schnideriana)、构树(Broussonetia papyrifera)、女贞(Ligustrun lucidumAit)、青檀(Ligustrun lucidum Ait)等阔叶树，另有少量雪松(Cedrus deodara)、圆柏(Sabina chinensis)、柳杉(Cryptomeria fortunei)等针叶树。

1)LiDAR数据获取：LiDAR数据获取装置为Optech Airborne Laser TerrainMapper(ALTM)3100激光传感器(加拿大Optech公司，产地加拿大多伦多)斑离散点云数据。每个激光返回点包含三维坐标值(x,y,z)、强度值(intensity)和返回类型(return type)等信息。

2)UAV数据获取：研究采用的无人机搭载平台为固定翼飞机，配备惯性测量单元(IMU，inertial measurement unit)(美国ADI公司，产地上海)获取飞机实时的三轴姿态角\角速率以及加速度。IMU型号为ADIS16355，集成了3轴加速度计及3轴陀螺，主要技术参数为：1)可控最高空速200m·s^–1；2)可控最高地速350m·s^–1；3)可控飞行高度20–5500m；4)导航精度(偏航距)≤±3m；5)定高控制精度≤±2m。无人机主要技术参数：1)翼展长度2710mm；2)机身长1956mm；3)最大推力8kg；4)总质量6kg；5)有效载荷4kg；6)飞行高度500m；7)飞行速度100km·h–1；8)抗风能力6级。无人机通过地面遥控系统进行飞行的控制和通讯。遥感数据来源于佳能EOS 5D Mark II高清相机(日本佳能公司，产地日本九州)获取的高分辨率、真彩色数码影像。获取时间为2012年8月12日。相机的主要参数为：1)250万像素CMOS传感器(1920×1280像素)；2)镜头：适合多种不同镜头，标准8.5mm镜头；3)图像大小与采集速度：每图0.7–0.9MB空间；每图2–5s；4)曝光时间：1/2500s；5)焦距：35mm；6)尺寸与质量：152mm×113.5mm×75mm，810g。航向和旁向重叠分别为30％和70％。

3)地面实测数据获取：结合研究区实际情况并通过前期对研究区林木林分特征变异性的调查分析，设计地面实测采用圆形样地方式，直径10m，30个随机分布于整个牛山研究区(图2)。样地调查时间为2014年5月。样地的圆心点通过TrimbleGeoXH6000Handhelds(Trimble Navigation Limited，2013，美国天宝公司，产地上海)手持GPS定位(集成GPS和GLONASS系统，并配备了双频GNSS天线)。通过接收JSCORS高精度实时差分信号，提升其水平定位精度到0.5m以内。样地的方向和倾角通过罗盘仪测定，边界半径通过皮尺量测。在各个样地中，记录胸径≥5cm的每棵树木的树种，并量测胸径、树高。胸径采用围尺量测，树高通过Vertex IV超声波测高器(瑞典Haglof公司，产地北京)获取。

样地选择参考了研究区内典型的树种组成、立地条件以及树木年龄(参考历史调查数据以及前期考察结果)。在样地内进行每木测量，选取样地Lorey’s树高、林分密度和蓄积量为反演的林分特征，计算方法如下：

Lorey’s树高(即断面积加权树高)计算公式如下：

H = \frac{Σ_{i - 1}^{k} h_{i} g_{i}}{Σ_{i - 1}^{k} g_{i}}

式中：h_i为第i株树高，g_i为第i株树木胸高断面积，k为样地内株数。

林分密度通过为单位面积上的树木数，以株·hm^-2为单位，通过样地实测林分密度转换每公顷林分密度得到。

蓄积量依据“江苏省立木一元材积式”(《江苏省森林资源规划设计调查操作细则(2007)》)计算，通过单木胸径(实测值)汇总到样地级蓄积量，再换算为单位面积(每公顷)蓄积量。结果如表1所示。

表1样地林分特征信息汇总

4)利用LiDAR数据生成DEM，DEM是获取相对树高的前提。数据处理的过程包括以下三步：(1)数据预处理，包括异常点剔除和重叠数据处理，首先去除高程噪声点，噪声点分两类，一类为高程负值点，一类为高程突变剧烈点。第二类点需参考研究区实际情况进行判断。其次去除LiDAR重叠航带的数据，由于LiDAR数据是分航带进行的，需要去除跨带的数据点，采用以大于激光点云平均间距的平滑滤波进行多次平滑消除重叠带。(2)点云数据分类，将点云数据分为裸地与非裸地两类，激光点云的分类采用滤波和分类的方法提取裸地(即地形)高程信息，将生成的裸地点云以点状矢量文件输出，再与UAV影像叠加进行目视纠正。(3)DEM生成，将分类后的裸地点采用反距离加权法进行1m×1m的插值，生成DEM。

5)通过UAV平台拍摄研究区共132幅数字影像，通过IMU实时记录每幅影像的经纬度、高程、翻滚度、俯仰和航向五个参数生成Exif图像文件，每幅影像分辨率为1920×1080。

6)借助Pix4D软件获取影像像对的特征点，平均每幅影像获取特征点为8760个，平均每幅图像建立4992匹配点对，通过光束法平差进行影像外方位元素的解析，完成图像标定，参与平差的二维点共673471个，生成三维点221175，平均投影误差为0.23像素。

7)利用已有三维点、影像定标和外方位元素解析生成研究区三维点云，在剔除明显的异常点后，通过平滑处理实现对“孔洞”的填充，再利用反距离加权法插值为与Lidar数据相适配的1m×1m的栅格生成DSM以便点云的归一化；同时通过影像拼接生成研究区UAV全幅影像，作为进行Lidar点云与UAV点云精度匹配的辅助验证底图。

8)点云的归一化指对点云高度处理，即距离地表的相对高度的获得，方式由DSM与DEM差值而来，归一化的结果为数字树冠高层模型(DCHM,Digtal CanopyHeight Model)。由于DSM和DEM二者数据平台不一致，数据差异性较大，主要有两点：(1)地理坐标不一致：由于传感平台不一致，本实施例中UAV影像数据采用wgs84坐标系统，而LiDAR采用西安80坐标系统。文章采用实地获取的17个典型地物的地面验证点坐标，通过基于布尔莎七参数法的三维坐标差转换模型，将DEM与DSM坐标统一，以便叠合运算。(2)点云数量有差距：研究区Lidar点云数量为22730，约0.17点·m^–2，UAV影像点云数量为772916，约5.56点·m^–2，两者数量差异较大，通过统一插值为1m×1m栅格数据满足空间运算前提。

9)变量提取：将DSM与DEM进行差值运算，获得DCHM，完成点云数据的归一化。建模数据变量从归一化后的点云中进行提取，选取高度分位数变量组(h10,h25,h30,h40,h60,h75,h85,h90)，即森林冠层高度垂直分布的分位数(第10,25,30,40,60,75,85,90分位数)；密度变量组(d10,d25,d30,d40,d60,d75,d85,d90)，即在相应百分位高度等级以上的点云在所有点云中所占的百分比；平均高度(havg)，即归一化高度的平均值；最大高度(hmax)，即归一化高度的最大值；最低高度(hmin)即归一化高度的最小值；同时，为了避免下层植被对冠层分析的干扰，去除了高度2m以下点云。

10)通过获取归一化点云提取的19个林分特征变量与实测样地4个关键变量之间线性关系进行精度的验证(表2)。验证的方式采取多元线性回归方法。建立回归模型运用逐步进入法(stepwise)和观察决定系数R²的变化情况来选择进入模型的合适变量，因样地数据较少，所有样地均参与模型构建。

表2多元回归分析变量表

11)精度验证

验证的指标选择：(1)决定系数(R²)：R²值越大，则因变量与自变量之间相关性越强,这是对回归直线拟合优度的检验。

R^{2} = 1 - Σ_{i = 1}^{n} {(x_{i} - {\hat{x}}_{i})}^{2} / Σ_{i = 1}^{n} {(x_{i} - {\overset{&OverBar;}{x}}_{i})}^{2} - - - (2)

式中，x_i为实测的某林分特征、为模型估算值、为实测特征值算术平均、n为样地数量。

(2)均方根误差(RMSE)：表达的是观测值与真值偏差的平方和观测次数n比值的平方根。

RMSE = \sqrt{\frac{1}{n} Σ_{i = 1}^{n} {(x_{i} - {\hat{x}}_{i})}^{2}} - - - (3)

(3)相对均方根误差(rRMSE)：定义为RMSE与实测结果算术平均值的比值，rRMSE越小，表明模型预测的效果越好。

rRMSE = \frac{RMSE}{{\overset{&OverBar;}{x}}_{i}} \times 100 % - - - (4)

通过相关性分析可见(图3)，归一化点云提取的特征变量与Lorey’s树高的相关度最高(R²＝0.58–0.95)，蓄积量次之(R²＝0.29–0.7)，胸高断面积与林分密度最低。联合提取变量中，均高(h_avg)、高分位高度分位数(h₇₅–h₉₀)和中高分位点云密度(d₆₀,d₉₀)与各林分特征之间的相关性也较高，而高度最值(h_min)与中低分位高度分位数及点云密度与各林分特征之间相关性较低。

12)联合提取估算模型构建及精度分析

通过多元回归分析，构建联合提取估算模型，模型演算变量与实测变量较为拟合(R²为0.36～0.96)(表3)。从模型参数来看，各林分回归方程参数在4个以内，自变量参数包括高度百分位数、点云密度和高度均值等，均与相应的林分特征呈较高的相关性。林分特征变量拟合中，Lorey’s树高拟合度最高，蓄积量次之，胸高断面积和林分密度较低。拟合结果与相关性分析结果一致。从相对均方根误差看，精度从高到低依次为Lorey’s树高(1.9％)、蓄积量(6.93％)、林分密度(13.74％)、胸高断面积(16.38％)。

表3联合提取估算模型及精度评价

注：H为Lorey’s树高；N为林分密度；G为胸高断面积；V为蓄积量；各模型自变量见表2。

由各林分特征散点图(图4)来看，Lorey’s树高识别精度最高(R²＝0.94)，拟合的方程与1:1验证线基本吻合，相关性显著；林分密度和胸高断面积精度相对较低(R²＝0.45和R²＝0.38)，蓄积量精度居中(R²＝0.55)，较为相关。

本申请通过结合UAV与LiDAR数据提取典型林分特征，利用特征点匹配算法和空中三角生成研究区三维点云作为数字表面模型结合LiDAR提取的数字地形模型计算树冠冠层模型，再结合地面实测林分调查数据与联合提取的点云变量构建多元回归估算模型并验证其精度。一方面冠层点云采用UAV数据，相对于LiDAR，具有灵活性高，廉价等特点，数据获取方便，具有长期定时重复观测的可能性。另一方面利用LiDAR的穿透性获取地形数据，而地形相对稳定，可以作为长期观测的数据基础。二者结合，各取所长，是一条廉价、灵活的森林状态监测方式，相关分析表明联合提取的点云变量与森林垂直信息(树高)具有很高的敏感性。

Claims

1.一种集成UAV与LIDAR的森林林分特征反演方法，其特征在于，包括以下步骤：

Lorey’s树高计算公式如下：

H = \frac{Σ_{i - 1}^{k} h_{i} g_{i}}{Σ_{i - 1}^{k} g_{i}}

4)利用LiDAR数据生成DEM；

11)精度验证

验证的指标选择：(1)决定系数R²：

R^{2} = 1 - Σ_{i = 1}^{n} {(x_{i} - {\hat{x}}_{i})}^{2} / Σ_{i = 1}^{n} {(x_{i} - {\overset{&OverBar;}{x}}_{i})}^{2}

RMSE = \sqrt{\frac{1}{n} Σ_{i = 1}^{n} {(x_{i} - {\hat{x}}_{i})}^{2}}

rRMSE = \frac{RMSE}{{\overset{&OverBar;}{x}}_{i}} \times 100 % - - - (4)

12)联合提取估算模型构建及精度分析

2.根据权利要求1所述的集成UAV与LIDAR的森林林分特征反演方法，其特征在于：步骤2)中，航向和旁向重叠分别为30％和70％。

3.根据权利要求1所述的集成UAV与LIDAR的森林林分特征反演方法，其特征在于：步骤2)中，IMU型号为ADIS16355，集成3轴加速度计及3轴陀螺，技术参数为：1)可控最高空速200m·s^–1；2)可控最高地速350m·s^–1；3)可控飞行高度20–5500m；4)导航精度(偏航距)≤±3m；5)定高控制精度≤±2m。

4.根据权利要求1所述的集成UAV与LIDAR的森林林分特征反演方法，其特征在于：步骤4)中，数据处理的过程包括以下三步：

(1)数据预处理，包括异常点剔除和重叠数据处理，首先去除高程噪声点，噪声点分两类，一类为高程负值点，一类为高程突变剧烈点；其次采用以大于激光点云平均间距的平滑滤波进行多次平滑消除重叠带；

(2)点云数据分类，将点云数据分为裸地与非裸地两类，激光点云的分类采用滤波和分类的方法提取裸地高程信息，将生成的裸地点云以点状矢量文件输出，再与UAV影像叠加进行目视纠正；

(3)DEM生成，将分类后的裸地点采用反距离加权法进行1m×1m的插值，生成DEM。