技术领域
[0001] 本
发明涉及森林资源管理与保护技术领域，具体涉及一种基于全波形LiDAR冠层剖面模型的树种识别方法。
背景技术
[0002] 精确的树种分类对于林业调查、
生物多样性研究以及模拟
指定树种的单木生长有重要意义。同时，这些信息也可以用于参数化森林生长模型和生态过程模型从而指导和优化森林资源管理。常规的森林树种调查方法主要依赖于野外调查及利用大比例尺航片判读等，其
精度往往不高，且难于在大区域上实用化推广。
激光雷达(LiDAR，Light Detection and Ranging)是通过发射
激光束击打物体表面并分析其返回
信号的一种主动遥感技术。通过LiDAR可获得高精度的地球表面及地表实体的高度信息，如地形和植被等可达到亚米级的垂直精度。现有研究表明，LiDAR可穿透森林冠层获得其三维结构特征，特别适合植被
覆盖度高且森林结构复杂地区森林信息提取。
[0003] 近年来基于小光斑LiDAR数据进行树种分类研究为： 等2009年在《Remote Sensing of Environment》第113卷上发表的“Classifying species of individual trees by intensity and structure features derived from airborne laser scanner data”，该研究在已有的单木
位置和冠幅信息
基础上，提取了单木的多个高度(如最大高，平均高和高度百分位数等)及统计(如高度分布的峰度和偏度等)特征变量，并据此对挪威寒带森林中的针叶和阔叶树种进行了分类；Li等2013年在《Agricultural and Forest Meteorology》第171卷上发表的“Classification of tree species based on structural features derived from high density LiDAR data”，该研究首先基于LiDAR点
云提取单木冠幅，然后在冠幅内汇总和计算三维质地，
叶片聚集尺度以及孔隙分布等
水平和垂直信息，并依据这些特征信息对北美温带森林进行树种分类。然而，以上分类方法都是基于点云特征变量的，并不包含完整的森林垂直结构信息(即存在信息缺失)。同时，点云数据也存在“盲区”(即两次反射回波需要有一定的垂直间距才能被系统区分开)，会影响林下低矮植被的提取和分类。
发明内容
[0004] 发明目的：针对
现有技术中存在的不足，本发明提出一种基于全波形LiDAR冠层剖面模型的树种识别方法，有效提高树种分类的总体精度，易于推广应用等特点。
[0005] 技术方案：为了实现上述发明目的，本发明采用的技术方案为：
[0006] 一种基于全波形LiDAR冠层剖面模型的树种识别方法，包括以下步骤：
[0007] 1)借助机载小光斑全波形LiDAR
传感器进行
数据采集；传感器记录每束激光脉冲返回的完整波形信息；
[0008] 2)LiDAR波形数据预处理：
[0009] A)噪声水平估计和数据平滑：首先把原始数据转换到
频率域，再将频率较高的低值部分作为噪声水平的判断标准；然后选用高斯
滤波器进行平滑。
[0010] B)高斯拟合分解及波形数据点云化：基于回
波数据是多个高斯函数的累加这一假设，对波形数据采用非线性最小二乘法进行拟合；然后通过局部最大峰值检测滤波
算法从处理后的波形数据上提取离散点云，每个离散点中记录了返回信号的
能量和振幅信息；
[0011] C)生成数字地形：首先对从波形数据中提取出离散点云进行分类，然后对末次回波进行Kraus滤波处理用以去除非地面点，最后使用滤波后的末次回波数据并借助自然邻近法插值生成数字地形模型DTM；
[0013] A)对地面以上点云进行中值滤波，然后将点云中的高度信息
栅格化生成数字表面模型DSM；将DSM减去数字地形模型，得到归一化植被高度CHM；
[0014] B)通过局部最大值法确定单木树顶所在位置；
[0015] C)单木冠幅的确定：首先以树顶为中心拟合16个半径方向上的冠幅剖面，然后计算到局部最小值的水平距离，最后将这些距离值进行水平方向上的平均从而得到冠幅半径；
[0016] 4)结合Weibull分布和三次样条函数模型，拟合单木垂直结构和返回能量信息的冠层剖面模型，并提取模型参数作为特征变量；
[0017] 5)使用
随机森林分类器进行树种分类：
[0018] A)随机森林分类是由很多
决策树分类模型{h(X,Θk),k＝1,2,...}组成的组合分类模型，且参数集{Θk}是独立同分布的随机向量，在给定自变量X下，每个决策树分类模型都由一票投票权来选择最优的分类结果；通过k轮训练，得到一个分类模型序列{h1(X),h2(X),…hk(X)}，再用它们构成一个多分类模型系统，其最终的分类决策函数为：
[0019]
[0020] 式中，H(x)表示组合分类模型，hi是单个决策树分类模型，Y表示输出变量，I为示性函数；
[0021] B)使用特征变量1-4和随机森林分类器进行分类，并借助总体分类精度和卡帕系数这两个定量化指标对分类精度进行评价；卡帕系数计算公式：
[0022]
[0023] 式中，r为总的类别数，xii为对
角线上的像元数，xi+和x+i是列和行的总像元，N是总像元数。
[0024] 步骤1)中，所述的LiDAR传感器为奥地利的Riegl LMS-Q680i。
[0025] 步骤1)中，所使用的遥感平台为运-5
固定翼飞机。
[0026] 步骤1)中，所述传感器的
采样间隔为1ns。
[0027] 步骤3)的B)中，通过圆形的搜索框在CHM上“
滑行”遍历栅格图层，通过逐个比较搜索框内的高度信息来确定最大值点；搜索框的半径通过高度信息和参数β0及β1来确定；β0及β1则通过地面实测数据拟合来获取，公式为：
[0028] CW(m)＝β0+β1×h2 (1)
[0029] 式中，CW为地面实测冠幅半径，h为树高(m)，β0和β1为模型参数。
[0030] 步骤3)的C)中，单木冠幅的确定借助冠幅半径来进行描述，方法为首先以树顶为中心拟合16个半径方向上的冠幅剖面，然后计算到局部最小值的水平距离，最后将这些距离值进行水平方向上的平均从而得到冠幅半径。
[0031] 步骤4)的具体步骤为：
[0032] A)对非地面点云数据进行垂直分布信息汇总：该点云提取自波形数据，且高度信息为相对高度；汇总的垂直间隔为0.3m，每个间隔内汇总点云的数量；
[0033] B)通过Weibull分布拟合点云的累积垂直分布，Weibull的
密度函数为：
[0034]
[0035] 式中，α和β为模型拟合参数，H为最大冠幅高度；其中，参数α为表达了冠层剖面分布的垂直尺度及位置的“特征变量1”，参数β为描述了分布宽度的增减的“特征变量2”；
[0036] C)对非地面点云数据的强度信息汇总：汇总的垂直间隔同样为0.3m，每个间隔内汇总点云中包含的强度信息，通过三次样条函数拟合强度分布信息，并提取峰值所在位置作为“特征变量3”，函数内包含的总面积作为“特征变量4”。
[0037] A)中的相对高度为绝对高度减去地形高度。
[0038] B)中的最大冠幅高度为最高点返回信息。
[0039] 步骤5)中的B)中的总体分类精度为对角线像元数/总体像元数×100％。
[0040] 有益效果：与现有技术相比，本发明的特色和创新点如下：1)所提取的特征变量融合了单木尺度上的冠层垂直结构和返回能量信息，从而“增强”了树种的识别效
力；2)由于提取的冠层剖面特征变量是对单木冠幅内木质结构和叶片分布的整体描述，该发明克服了传统特征变量(如与高度相关的多个因子)彼此相关度高，以至于难以筛选和判别重要性的
缺陷。同时，该发明也利于特征变量的机理解释和方法移植(如在多个森林类型中的应用)。实验验证结果表明，通过本发明对北亚热带天然次生混交林的主要树种进行树种分类，与其他使用遥感方法(注：对比条件为：单一遥感数据源；且在类似复杂森林类型内)进行树种分类的方法相比，在四个主要树种的分类级上总体精度提升了9％左右；Kappa系数提升了0.1左右。
附图说明
[0041] 图1是研究区的6个主要树种照片及其球面镜影像图；
[0042] 图2是单个样地内的树顶和冠幅提取图；
[0043] 图3是部分运算方法和结果的
可视化图；其中，a是原始波形数据；b是基于波形提取的点云数据；c是以树顶为中心拟合16个半径方向上的冠幅剖面；d是Weibull分布拟合点云的累积垂直分布；
[0044] 图4是汇总样木点云垂直分布信息的各树种冠层剖面图。
具体实施方式
[0045] 下面结合具体
实施例对本发明作进一步的说明。
[0046] 实施例1
[0047] 一种基于全波形LiDAR冠层剖面模型的树种识别方法，以针对一个北亚热带天然次生混交林为主要森林类型的林区内树种分类为例。林区海拔20-261m，面积约1100公顷。主要树种为针叶的
马尾松(Pinus massoniana)、杉木(Cunninghamia lanceolata)和湿地松(Pinus elliottii)，以及阔叶的麻栎(Quercus acutissima)、枫香(Liquidambar formosana)和冬青(Ilex chinensis)(见图1)。林区内根据树种组成、林龄和立地状况等布设了12个方形样地(30×30m)，每个样地内人工判别了单木树种，并实测了胸径、树高和冠幅等森林参数；样地的中心通过差分GPS进行定位，每木的相对位置(即距离样地中心的水平距离和方向角)通过皮尺和森林罗盘仪测定(最后换算为每木的绝对位置坐标)。具体过程如下：
[0048] 1)LiDAR数据获取
[0049] 借助奥地利Riegl LMS-Q680i的机载小光斑全波形LiDAR传感器进行数据采集。遥感平台为：运-12飞机(飞行高度900m，速度为：55m/s)。获得的LiDAR波形数据参数为：脉冲发射频率400kHz，扫描频率114lines/sec(扫描角±30°)，激光脉冲的光斑半径为45cm。传感器记录了每束激光脉冲返回的完整波形信息，采样间隔为1ns。获得数据的
2
脉冲间距为0.48m，脉冲点密度2.1pulse/m(航向重叠区域的脉冲点密度更高)。
[0050] 2)LiDAR波形数据预处理
[0051] a)噪声水平估计和数据平滑。首先把原始数据转换到频率域，再将频率较高的低值部分作为噪声水平的判断标准。然后选用高斯滤波器进行平滑(
内核设置：FWHM＝2.35×σ)，这是由于高斯滤波器在有效平滑数据的同时，还可以最大限度地保持原有曲线的趋势。
[0052] b)高斯拟合(分解)及波形数据点云化。基于回波数据是多个高斯函数的累加这一假设，对波形数据采用非线性最小二乘法进行拟合。然后通过局部最大峰值检测滤波算法(即通过移动窗体逐个判断某点周边的最近4个点的脉冲强度信息)从处理后的波形数据上提取离散点云，每个离散点中记录了返回信号的能量和振幅信息(单木内点云提取效果示意见图3.b)。
[0053] c)生成数字地形。LiDAR数据高度归一化的目的是为了得到去除了地形影响的“真实”植被高度，通常采用原始LiDAR数据高度信息减去地形高度得到。因此，精确生成数字地形模型(DTM)是计算归一化植被高度的重要前提。首先对从波形数据中提取出离散点云进行分类，然后对末次回波进行Kraus滤波处理用以去除非地面点(在此基础上进一步采用中值滤波进行平滑，窗体大小3×3)，最后使用滤波后的末次回波数据并借助自然邻近法插值生成数字地形模型(输出空间
分辨率为1m)。
[0054] 3)单木定位和冠幅提取
[0055] a)对地面以上点云进行中值滤波(窗体大小为3×3)，然后将点云中的高度信息栅格化生成数字表面模型(DSM)。将DSM减去数字地形模型(DTM)从而得到归一化植被高度(CHM)。
[0056] b)通过局部最大值法确定单木树顶(即最高点)所在位置。即通过圆形的搜索框在CHM上“滑行”遍历栅格图层，通过逐个比较搜索框内的高度信息来确定最大值点。搜索框的半径通过高度信息和参数β0及β1来确定。β0及β1则通过地面实测数据拟合来获取。其公式为：
[0057] CW(m)＝β0+β1×h2 (1)
[0058] 其中CW为地面实测冠幅半径，h为树高(m)，β0和β1为模型参数。
[0059] c)单木冠幅则是借助冠幅半径来进行描述。其方法为首先以树顶为中心拟合16个半径方向上的冠幅剖面(见图3.c)，然后计算到局部最小值的水平距离，最后将这些距离值进行水平方向上的平均从而得到冠幅半径。根据所测森林结构参数，选取样地中的典型样本：马尾松(75棵)、杉木(69棵)、湿地松(65棵)、麻栎(65棵)、枫香(57棵)和冬2
青(45棵)，拟合模型得到β0为1.781，β1为0.029(模型的R 为0.63)。单木冠幅的确定借助冠幅半径来进行描述，方法为首先以树顶为中心拟合16个半径方向上的冠幅剖面，然后计算到局部最小值的水平距离，最后将这些距离值进行水平方向上的平均从而得到冠幅半径。将LiDAR提取单木与地面实测数据进行空间位置对应“链接”(即LiDAR估算和地面实测的树冠中心水平位置在2m范围内则认定为匹配)，得到单木提取正确率为76％(提取效果见图2)。单木树高的RMSE(均方根误差)为0.63m，冠幅的RMSE为0.55m。共提取(430棵)有效树:马尾松(85棵)、杉木(81棵)、湿地松(70棵)、麻栎(72棵)、枫香(73棵)和冬青(50棵)。
[0060] 4)结合Weibull分布和三次样条函数模型，拟合单木垂直结构和返回能量信息的冠层剖面模型，并提取模型参数作为特征变量。
[0061] a)对非地面点云数据进行垂直分布信息汇总。该点云提取自波形数据，且高度信息为相对高度(即绝对高度减去地形高度)。汇总的垂直间隔为0.3m,每个间隔内汇总点云的数量。
[0062] b)通过Weibull分布拟合点云的累积垂直分布(见图3.d)，Weibull的密度函数为：
[0063]
[0064] 其中α和β为模型拟合参数，H为最大冠幅高度(即最高点返回信息)。该密度函数的导数即为单木的冠层剖面分布模型。其中参数α即为“特征变量1”，表达了冠层剖面分布的垂直尺度及位置；而参数β即为“特征变量2”，描述了分布宽度的增减。
[0065] c)对非地面点云数据的强度信息汇总。汇总的垂直间隔同样为0.3m，每个间隔内汇总点云中包含的强度信息。通过三次样条函数拟合强度分布信息，并提取峰值所在位置(作为“特征变量3”)，函数内包含的总面积(“特征变量4”)。
[0066] 5)使用随机森林分类器进行树种分类。
[0067] a)随机森林分类是由很多决策树分类模型{h(X,Θk),k＝1,2,...}组成的组合分类模型，且参数集{Θk}是独立同分布的随机向量，在给定自变量X下，每个决策树分类模型都由一票投票权来选择最优的分类结果。通过k轮训练，得到一个分类模型序列{h1(X),h2(X),…hk(X)}，再用它们构成一个多分类模型系统，其最终的分类决策函数为：
[0068]
[0069] 其中，H(x)表示组合分类模型，hi是单个决策树分类模型，Y表示输出变量，I为示性函数。
[0070] 使用特征变量(1-4)和随机森林分类器进行分类，并借助总体分类精度(即对角线像元数/总体像元数×100％)和卡帕系数这两个定量化指标对分类精度进行评价。卡帕系数计算公式：
[0071]
[0072] 式中，r为总的类别数，xii为对角线上的像元数，xi+和x+i是列和行的总像元，N是总像元数。
[0073] 将以上4个特征变量作为输入参数代入随机森林分类模型(决策树数量：100；每个分叉设置4个变量)，并借助总体分类精度和卡帕系数等定量指标对分类精度进行评价。3个分类等级的混淆矩阵见表1-3。6个树种分类结果(表1)(总体精度＝64.93％；Kappa系数＝0.579)，4个树种分类结果优于6个树种(表2)(总体精度＝72.35％；Kappa系数＝0.631)，森林类型(针阔叶)分类结果比树种分类结果最好(表3)(总体精度＝83.2％；
Kappa系数＝0.664)。
[0074] 表1 6个树种分类结果验证混淆矩阵
[0075]
[0076] 注：像元数已转换为百分比。
[0077] 表2 4个主要树种分类结果验证混淆矩阵
[0078]
[0079] 注：像元数已转换为百分比。
[0080] 表3 森林类型(针阔叶)分类结果验证混淆矩阵
[0081]
[0082]
[0083] 注：像元数已转换为百分比。
