机载激光雷达森林资源调查系统的设计与试验

瞿 帅，张晓丽，朱程浩，霍朗宁，刘会玲

(北京林业大学 精准林业北京市重点实验室，北京 100083)

森林是陆地生态系统的主体，森林及其变化对陆地生物圈及其他地表过程有着重要影响[1]。森林资源的可持续开发利用依赖于有效的森林资源调查和经营，而传统的森林资源调查方式需要投入大量的人力和物力，调查周期长，数据更新速度远远不能满足当今信息化时代对森林资源调查数据的实时更新和获取需求[2-4]。

随着遥感技术的发展，森林资源调查拥有了强大的技术支撑[5]。近年来，光学遥感技术已广泛应用于森林参数(如森林蓄积量、叶面积指数、生物量等)的提取及反演[6-11]，但是光学遥感受多种因素的影响，如云、混合像元等[12],而且由于光学遥感不具备三维空间探测的能力，从而限制了光学遥感在森林结构参数提取及反演方面的应用。

激光雷达是一种主动方式的遥感探测技术，具有较强的穿透能力，能够用来探测森林空间结构及林下地面信息，具有光学遥感无可比拟的优势[13],在森林资源调查及经营管理方面具有重要的应用价值[14]。

现阶段利用机载激光雷达进行森林资源调查的工作已经展开，但对于获取的大数据量的激光雷达点云数据面向林业需求的处理工作往往比较滞后，如何有效利用机载激光雷达数据及时获取成果以保障森林资源调查工作的如期完成是一个亟待解决的问题。因此，本研究将结合实际生产需要，利用机载激光雷达数据结合地理信息系统技术进行森林资源调查系统的设计与试验，以期实现自动化、高精度的森林资源调查，满足森林资源数据库及时更新的需求。

1 系统分析与设计

1.1 功能需求分析

林业资源调查获取的机载激光雷达点云数据一般数据量较大，系统需要在足够大的内存下，支持千万级别数量的点云数据载入，并进行流畅的数据渲染显示；能够加载影像；支持不同格式点云数据的转换；可以实现点云去噪并进行高精度的点云滤波处理，生成数字高程模型；从单木尺度实现单木定位、单木树高、冠幅的提取，从林分尺度实现株树密度、林分平均树高的提取。

1.2 系统结构设计

本系统以QT为图形界面框架，结合VTK三维图形引擎库与GDAL影像处理库，在Visual Studio 2013集成开发环境下采用C++语言开发。系统结构设计使用3层结构，分别是应用层、逻辑层、数据层。

应用层也可以称为表现层，主要表现为人机交互，是对用户公开的部分；逻辑层主要实现数据组织管理、各种数据处理算法；数据层主要实现不同数据格式的数据与本系统进行数据交换的接口。本系统所使用的数据以文件形式提供，系统为不同类型及格式的数据设计相应的数据解析接口，并封装成数据操作类库，方便本系统对外部数据进行操作管理及分析，并具备较高的可扩展性。3层结构间的关系具体表现为：应用层调用逻辑层实现的相关数据处理算法完成应用层的功能实现，逻辑层通过数据操作类库完成对数据层相关数据的组织及管理。系统结构示意图如图1所示。

图1 系统结构

机载激光雷达森林资源调查系统主要分为数据输入、地形信息提取、林木参数提取、成果输出4大功能模块。系统主界面如图2所示。

1.3 核心模块实现

1.3.1 地形信息提取模块 机载LiDAR能够穿透植被遮挡，直接获取地面点信息，因而成为数字高程模型(DEM)获取的首选方式，而原始点云数据中还包括建筑、植被等非地面信息，要获取地面点，就必须采取一定的算法将非地面点从原始点云中剔除，这个过程称为点云滤波[15]。常规点云滤波算法主要包括数学形态学滤波[16-17]、逐步加密滤波[18-20]等算法。针对林区地形特征较为复杂的特点，本系统中实现了一种渐近不规则三角网加密滤波算法，算法流程如图3所示。

图2 系统主界面

图3 算法流程

算法基本步骤：

1)对原始点云数据进行中值滤波去噪，去除由于系统误差或者飞鸟激光反射点等噪声点云，避免对后续步骤算法的影响。

2)点云格网化，设定一格网尺寸，根据点云中空间点的X、Y坐标将去噪后点云划分入规则格网内，并将每一规则格网内高程最低点作为初始地面种子点。

3)将初始地面种子点作为构建初始不规则三角网的种子点，在构建初始不规则三角网过程中，引入坡度阈值判断，即在构建三角网时先进行相邻构网节点的坡度计算，如果小于坡度阈值(此处坡度阈值应设定较大值)则该点加入构网，否则标记为非地面点，以后步骤中不再参与运算，最终获得初始不规则三角网。

4)在初始不规则三角网基础上，对未划入地面点云的点进行坡度及高差阈值的判断，如果满足条件则划入地面点，对地面不规则三角网进行迭代加密。

5)重复步骤4，当不再有新的地面点加入三角网时，停止迭代，得到满足条件的三角网点云，即为最终地面点云。

在软件实现方面，本系统中将此算法封装成一个算法类，格网尺寸、坡度阈值、高差阈值等参数可在软件滤波算法界面上进行设定，作为算法传入参数，执行滤波最终得到地面点云。为了根据获得的地面离散点云得到数字高程模型，系统采用GDAL影像处理库中离散三维点反距离加权插值算法将地面点云插值形成DEM影像，以标准GeoTIFF格式输出。插值函数为GDALGridCreate，插值算法设置为GGA_InverseDistanceToAPower。

1.3.2 林木参数提取模块 此模块主要包含单木参数提取及林分参数提取功能，单木参数主要包含单木树高、冠幅以及单木位置，林分参数主要包括株树密度及林分平均高。

1.3.2.1 单木参数提取 机载LiDAR点云数据已成功用于森林垂直结构参数及水平结构参数的估测[21]，可以直接提取单木树高及冠幅。在点云滤波处理的基础上可以得到林区点云数据包含的地面点云和非地面点云，林区非地面点云基本相当于植被点云，对原始植被点云进行高程归一化处理，则输出点云中每个空间点的高程即为该点相对于地面的绝对高度，然后利用归一化点云进行单木点云分割(本系统中实现了一种基于归一化植被点云高度分层的K-Means聚类点云分割算法)，最终依据单木点云分割结果进行树高和冠幅的提取。算法流程如图4所示。

图4 算法流程

算法基本步骤：

1)点云归一化：将点云滤波得到的植被点云数据以水平X、Y坐标为依据，根据生成的DEM影像分辨率尺寸划分到不同格网，每个落在格网内的植被点高程减去对应DEM格网的值，即得到归一化后的植被点云，此时植被点云中点的高度就相当于植被高度。

2)单木点云分割：按照归一化后植被点云高度分布情况，将植被点云分成若干层，在每一层设定一邻域检测窗口尺寸，进行高度局部最大值检测，对每一层的植被点云以局部最大值点为初始聚类中心进行三维空间点K-Means聚类，在每层点云聚类完成之后，从最上层开始，判断相邻上下2层各聚类点云中各聚类中心间水平距离是否小于设定的聚类中心距离阈值，如果小于阈值则合并上下2层中对应的聚类点云，直到所有层比较合并完毕，最终得到单木分割点云。

3)参数提取：对得到的单木分割点云，统计每个单木点云的最高点高度，将其作为提取的树高，并将该点的水平坐标作为单木定位坐标，计算每个单木分割点云投影到水平面上的凸包面积，以该面积为树冠近圆形投影面积，计算圆直径作为单木平均冠幅。

在本系统中将点云归一化、单木点云分割以及树高、冠幅提取分步进行，针对每个子功能封装为一个算法类，并且设置了不同的软件界面窗口，用于设定各子功能参数，最终实现了单木树高及冠幅参数的提取，提取结果以表格形式输出，记录了每棵树的位置信息、树高、平均冠幅。

1.3.2.2 林分参数提取

1)株树密度：在实际林业调查中常使用圆形样地法[22-23]来确定局部区域的林木株树密度，该方法可以提高调查效率。本研究中以方形样地为调查区域，根据样地单木分割结果，可以统计出样地的林木株树，通过式(1)计算出样地的株树密度值。

(1)

式中,N为株树密度(株/hm2)；n是样地内树木株数；S是样地面积(hm2)。

2)林分平均高：本系统中实现的林分平均高提取方法，以庞勇[24]等的研究为依据，样地归一化植被点云数据上四分位数处的高度与实测树高相关性高。

先选定一定数量的样地调查数据作为训练样本，计算样地归一化植被点云上四分位高度，具体计算为：对样地范围内的归一化植被点云根据高度进行排序，然后计算总高度的上四分位处的高度[20]，即为样地归一化植被点云上四分位高度，再将其与实测林分平均高建立线性回归方程，然后可以用预留检验样本进行精度验证。最终根据样地点云及回归方程实现待测样地的林分平均高提取。

其中，实测林分平均树高计算采用断面积加权计算方法，计算公式如下：

(2)

式中，H为林分平均高，hi为第i棵树的树高，gi为第i棵树的胸高断面积，k为林分株数。

2 试验区案例研究

2.1 试验区概况

试验区位于甘肃省张掖市大野口森林区域(97°20′-103°50′E，36°40′-39°50′N)，位于祁连山自然保护区，海拔2 700～3 200 m。样地内树种为天然青海云杉纯林，分布茂密，林龄组成结构主要为成熟林。林下植被种类单一，地表覆盖物主要为苔藓。

2.2 数据采集及预处理

2.2.1 LiDAR数据获取 本次试验区的机载激光雷达数据由RIEGL公司研制的LiteMapper 5600系统进行测量，配置的激光扫描仪为Riegl LMS-Q560，工作波长1 550 nm，光斑平均直径约为0.38 m。飞机的相对航高约700 m，整个测区点云的平均密度为3.43个/m2，地表点的定位精度为水平方向0.1 m、垂直方向0.03 m。LiDAR点云数据采用WGS84坐标系和UTM投影(北半球6度第48带)。

2.2.2 样地调查数据 设立一边长100 m正方形范围的大尺寸样地，为方便测树调查，将其划分成16个子样地，在试验区已知的大地控制点上建立GPS基准站，使用静态差分GPS对大样地四角点及各子样地的中心坐标和四角点位置坐标进行精确定位，以子样地为单位进行每木测树调查，数据主要包括2个部分：一是林木实测数据，包含实测样地内单木胸径、树高以及单木冠幅(南北、东西)等。二是地面定位点数据，采用全站仪对样地内所有挂牌林木的地面位置进行精确测定。利用大样地及子样地角点定位坐标，使用系统点云裁剪功能获取大样地及每个子样地的原始点云数据。大样地点云数据如图5所示，以高程渲染显示。

图5 点云数据

2.3 地形信息精度验证

本次研究对系统的数字高程模型提取值进行了精度验证，以16块子样地中心位置实测高程为真值，将该系统生成的数字高程模型对应位置的高程值与实测数据作了对比(表1)。对比数据分析可得，系统高程提取值的最大偏差为1.102 m，平均偏差为0.737 m。

2.4 林木参数精度验证

为验证软件系统提取的单木树高及冠幅，将该试验区实测数据与软件提取结果进行了对比。根据统计，大样地共调查1 465株树，去除枯立木、死木和缺少单木位置的单木，同时对于融合的树进行拆分，共获得1 445株树，大样地单木分割结果共提取1 223株树(图6)。根据统计输出单木提取位置信息依据最近邻匹配原则，在提取单木中选取距离实测单木位置最近单木作为其匹配单木，将匹配单木实测树高与提取值进行对比，可得单木树高平均相对误差为13.66%；将实测单木的东西、南北向冠幅平均值与提取单木树冠进行对比，可得单木冠幅的平均相对误差为14.18%(表2)。

表1 软件提取高程与实地测量值比较

图6 单木分割结果

对于各子样地中株树密度的精度验证，根据大样地单木分割结果以及16个子样地中各匹配单木，将其与16个子样地实地调查结果进行对比，可得软件提取的株树密度与实地调查株树密度值平均相对误差为5.83%；而在软件提取林分平均高参数时，选取1～10号共10块子样地作为训练样本，构建样地归一化植被点云上四分位高度与实测林分平均高的线性回归方程为Y=0.727 7X+5.659 9，R2为0.752 9，其中Y为实测林分平均高，X为样地归一化植被点云上四分位高度，利用剩余6个子样地点云结合回归方程计算林分平均高并与实测数据进行对比，可得软件提取的林分平均高与实测林分平均高的相对平均误差为9.86%(表3)。

表2 单木参数提取与实地测量值比较

表3 林分参数提取与实地测量值比较

3 结论与讨论

基于机载激光雷达的森林资源调查系统充分利用激光雷达可以探测森林垂直空间结构的能力，体现出了其在林业调查领域的优势。与传统调查方式相比，该系统可以快速、高效的获取调查区域的地形信息及相关林木参数信息。以甘肃省张掖市大野口关滩森林站超级样地试验区为例，进行机载激光雷达森林资源调查系统的验证研究，该系统可以有效获取试验区域地形信息，提取单木树高、冠幅的平均相对误差分别为13.66%、14.18%，提取株树密度、林分平均高的平均相对误差分别为5.83%、9.86%。

针对本系统进行林木参数提取方面，由于样区内点云密度不高，去除样区地面点之后，实际林木点云更加稀疏，而且激光点不能够保证扫描到树顶，因而会造成利用单木分割算法提取计算的单木树高普遍偏低，同样也会造成单木树冠提取结果偏低，而当有高大树木遮挡邻近较矮树木时，单木分割时，可能会将其划分为一棵树木，有可能会造成个别提取的单木冠幅大于实测值。对于林分参数，株树密度的精度与单木分割的结果有直接的关系，而以样地归一化植被点云上四分位高度结合回归方程计算林分平均高的精度对比单木树高提取精度有所提高，与激光点云不一定能够扫描到树顶从而影响单木树高提取精度的分析相印证。

现阶段，本系统尚有诸多不足之处，仅仅针对林分结构单一的针叶纯林进行了林木参数提取及验证，对于不同植被类型以及林分结构较为复杂的情况，系统中设计实现的单木分割算法具有一定的局限性，林木参数提取的精度会有较大程度的降低。因此后期考虑针对不同植被类型，将机载激光雷达系统搭载的光学相机获取的高分辨率光学影像与激光雷达数据有机结合，特别是针对阔叶林，单纯利用激光点云数据难以进行较高精度的林木参数提取，考虑以光学影像进行单木树冠分割，然后以单木树冠分割结果限定单木激光点云的划分，再进一步利用限定后的点云进行林木参数的提取，从而提高林木参数的提取精度，这将是下一步研究的重点方向。