基于机载激光雷达的国外松片林蓄积量反演研究

马 骧，姜文龙，闫保银，金 悦，崔 立

(1.常州市金坛区林业技术指导站，江苏 金坛 213200； 2.南京国图信息产业有限公司，江苏 南京 210037)

森林是地球生态系统的重要组成部分，在给人类带来林副产品产生经济效益的同时，还具备涵养水源、调节区域气候条件，实现碳平衡、碳循环等功能，产生生态效益。在IPCC发布的《全球温控1.5 ℃特别报告》中明确指出，“林业是实现1.5 ℃温控目标的重要路径之一，通过林业固碳除碳等技术，可实现温室气体的大幅减少，在促进森林生态系统恢复，提高基于生态系统和社区的适应性以及湿地管理方面潜力巨大。”在此背景下，如何快速准确地获取森林资源信息，实现动态监测，为森林保护、森林经营等提供依据，成为当前森林资源管理的重要任务之一。

激光雷达(LiDAR)的出现为森林资源调查方式的转变提供契机。作为一种主动遥感技术，具有精度高、成本低等优势，早在20世纪80年代已用于林业研究[1]。1984年，Nelson等[2]通过研究激光雷达与树木冠层郁闭度关系，指出两者间相关性显著；MacLean等[3]通过研究树冠垂直剖面与森林蓄积量相关性，指出2者的对数呈线性相关。进入21世纪后，越来越多学者也相继证明激光雷达对于估测树高、胸高断面积、生物量等森林参数准确性较高[4-6]。Morsdorf等[7]采用林木冠幅分割算法对瑞士国家森林公园针叶林进行分割，再生成种子点进行聚类分析，获得树冠直径后按圆面积公式反演了林木胸径；刘清旺等[8-9]基于机载激光雷达技术，使用多元线性回归、随机森林模型等方法，对森林平均高、生物量、郁闭度等参数进行提取分析，其反演结果较光学遥感而言精度更高；Zhao等[10]采用2种尺度不变模型估测生物量，结果表明该2种方法R2范围在0.80—0.95之间，可靠性较高。本文以江苏省常州市金坛区碳汇计量山区监测点中的国外松片林作为研究对象，进行机载激光雷达蓄积量反演试验，评价这一技术在江苏省苏南地区国外松林的适用性，以及在蓄积量反演方面的可能性，为今后森林资源调查新技术应用提供参考。

1 研究区概况与数据获取

1.1 研究区概况

研究区位于江苏省常州市金坛区4 km×4 km碳汇计量山区监测点范围内，国家2000大地坐标系(China Geodetic Coordinate System 2000, CGCS2000)范围坐标：西北至(20 718 118.46，3 517 999.08)，西南至(20 718 118.46，3 513 999.07)，东南至(20 722 118.48，3 513 999.07)，东北至(2 0722 118.48，3 517 999.08)。区内地貌以丘陵为主，为北亚热带季风区，年均气温15.3 ℃四季分明，雨量充沛，年降水量1 063.5 mm。日照充足，日照率46%，无霜期228 d，年平均湿度为78%。

研究区面积为16 km2，根据2018年碳汇计量监测成果，区内森林总蓄积量为48 486.04 m3，分布树种主要有湿地松(PinuselliottiiEngelmann)、杉木(Cunninghamialanceolata)、香樟(Cinnamomumcamphora)、女贞(Ligustrumlucidum)、柳杉(Cryptomeriafortunei)等。本次研究对象位于4 km×4 km碳汇计量山区监测点(见图1)，77，250号国外松小班范围内。其中77号小班用于激光雷达森林参数反演，250号小班用于反演数据验证。

图1 试验区范围

1.2 LiDAR数据获取

分别于2021年4月16日、5月12日，利用激光雷达在4 km×4 km碳汇计量山区监测点，77，250号小班范围内，进行机载激光雷达进行试验数据采集。使用大疆经纬M300 RTK作为飞行平台，飞行航速为6 m/s，总覆盖面积约为10 hm2，飞行时天气晴朗，少量高云分布，对于LiDAR飞行影响较少。

在本次飞行中，激光雷达系统采用大疆禅思L1(DJIL1)，集成Livox激光雷达模块、高精度惯导、测绘相机、三轴云台等模块。飞行时设置航高80 m，旁向重叠度50%。镜头角度-90°，激光扫描方式为三回波重复扫描，类似于传统线扫激光雷达，能获取更均匀、精度更高的扫描结果，也更稳定。

1.3 地面样地数据获取

2021年5月对77，250号小班进行外业调查。

本次地面实测调查分为大样地与小样地。大样地以77号小班范围为界，在其中按照20 m×20 m划分小样地，共计34块。在样地范围内记录树种、胸径、树高等参数，所得数据将以小样地为单位进行林分参数反演。

在250号小班内，按照《森林资源规划设计调查技术规程》[11]中小班测树因子调查的相关要求，深入小班内部，选择具有代表性的调查点，以调查点为中心布设样地并每木检尺，记录树种、胸径、树高、郁闭度等测树因子，按照《江苏省主要树种一元材积表》计算公顷蓄积，所得数据用于林分参数反演的数据可靠性验证。

图2 样地设置

2 研究方法

2.1 样地数据处理

根据样地内每木检尺获得的胸径，对照《江苏省主要树种一元材积表》，查询树种为“马1”一列，以及胸径对应的径阶，填写单株立木的材积，最后再统一求和，获得样地内总蓄积量。计算公式为

M树种=V树种/ ΣAi；

M全小班=∑(M树种×A)。

式中，A：小班面积；Ai:第i块小样地面积；V树种：某树种各样方蓄积之和；M树种：某树种小班每公顷蓄积；M全小班：全小班蓄积。

对于77号小班，其范围内的34块小样地需逐块计算平均胸径，用于后期反演模型的构建，250号小班需求得小班平均胸径、每公顷蓄积量，用于反演结果验证。

2.2 LiDAR数据处理

在林分反演模型构建之前，首先需将点云去噪，其次对点云进行地面点分类，以此分离地面点与地物点，然后分别生成数字地面模型(DTM)、数字表面模型(DSM)和冠层高度模型(CHM)，并基于CHM生成种子点，在点云归一化处理(NPC)，消除地形影响后，最后进行单木分割。本次对于LiDAR数据的处理均在LiDAR 360软件中完成。

地面点分类是点云数据处理的基础操作，本文采用改进的渐进加密三角网滤波算法(IPTD)[12]分类地面点。其原理是：对于不含建筑物的点云，以默认值作为格网大小，取格网内最低点作为起始种子点，构建初始三角网。遍历所有待分类点后，查询各点水平投影所落入的三角形，并计算点到三角形的距离d及点到三角形3个顶点与三角形所在平面所成角度的最大值，分别与迭代距离与迭代角度进行比较，如果小于对应阈值，则将此点判定为地面点，并加入三角网中。重复此过程，直至所有地面点分类完毕。

地面点分类后，需经插值处理生成DSM，DTM。LiDAR 360软件中提供了反距离权重插值(IDW)、克里金插值(Kriging)和不规则三角网插值(TIN)方法。本文采用TIN插值法，分别生成DSM和DTM，两者作差后生成CHM。利用生成的DTM高程对点云作归一化处理，消除地形影响。设置地面高度、最小树高再进行单木分割。

2.3 基于单木分割的反演模型构建

胸径是计算林木蓄积量的重要参数之一，若能准确获取胸径值，则可根据地方相应的一元材积表求算蓄积量。

树高、冠幅等参数可通过机载激光雷达直接测得，而胸径则需通过建立回归方程的方式间接获取[13]。本文拟参考刘清旺[14]、何祺胜等[15]的研究方法，以34块小样地为单位，统计各样地内LiDAR估测的单木平均树高、平均冠幅等参数，作为自变量，实测平均胸径作为因变量，建立实测胸径与LiDAR估测参数的一元回归方程，间接估测胸径。共设置3组：

(1)实测胸径与LiDAR树高

(2)实测胸径与LiDAR冠幅

(3)实测胸径的自然对数与LiDAR树高的自然对数

一元回归方程建立后，根据相关系数，确定最优估测方程。

2.4 反演模型可靠性评价

反演模型采用拟合优度(R2)、均方根误差(RMSE)和相对均方根误差(rRMSE)进行可靠性评价。

R2表示模型的拟合优度，一般R2值越高，表明该模型的拟合度就越高。计算公式如下：

RMSE表示实际值与预测值之间的均方根,误差越小，则方程拟合效果越好，计算公式为

rRMSE为相对均方根误差，其值应当越小越好，计算公式为

3 结果与分析

3.1 样地实测数据分析

77号小班的LiDAR数据经数据处理后共获得单木309株，与实测株数381株相比单木识别率达81%。在小班内按20 m×20 m大小划分小样地34块，以小样地为单位，分别统计LiDAR估测数据(含树高、冠幅、冠幅面积和冠幅体积)和实测胸径的平均数、中位数、最大值、最小值以及标准差，结果如表1所示。

表1 34块小样地林分因子统计

在LiDAR估测数据中，平均树高为18.7 m，其最大值为21.1 m，最小值为14.2 m，标准差达到1.5；平均冠幅为7.1 m，其最大值为12.8 m，最小值为3.9 m，标准差2.0；平均冠幅面积与体积标准差分别为24.6和128.3，数据离散程度较高。

在实测胸径方面，平均胸径为29.4 cm，最大值为33.8 cm，最小值为20 cm，标准差为2.6。

3.2 胸径反演模型构建

建模前需对34块小样地数据进行样本异常值诊断和处理。经分析后发现，第16号小样地的LiDAR估测平均树高较低，其原因是由于该样地范围内林下分布植被，导致单木分割时错将林下植被作为国外松进行分割，降低了样地平均树高值。因此需重新设置最小树高，将林下植被进行过滤，修正平均树高后进行胸径反演。

基于LiDAR估测与样地实测数据结果，按照基于单木分割的反演模型构建中的方法，设置3组反演模型：

(1)实测胸径与LiDAR树高

(2)实测胸径与LiDAR冠幅

(3)实测胸径的自然对数与LiDAR树高的自然对数

反演结果如图3、表2所示。

表2 实测胸径与估测胸径回归方程

在第(1)组模型中，以LiDAR估测树高作为自变量，实测胸径作为因变量建立反演模型。结果显示：获得方程y=26.7 lnx-48.7，拟合优度R2=0.741 1，均方根误差(RMSE)=0.004。

在第(2)组模型中，以LiDAR估测冠幅作为自变量，实测胸径作为因变量建立反演模型。结果显示：获得方程y=- 0.233 9x2+ 3.142 3x+ 19.797，拟合优度R2=0.603 1，均方根误差(RMSE)=0.008。

在第(3)组模型中，首先将LiDAR估测树高与实测胸径分别进行自然对数转换，再以LiDAR估测树高的自然对数作为自变量，实测胸径的自然对数为因变量，建立反演模型。结果显示：获得方程y=- 2.184 9x2+13.569x-17.605，拟合优度R2=0.806 3，均方根误差(RMSE)=0.001。

图3 3组胸径反演模型

综上所述，在3组反演模型中，分别将LiDAR估测树高和实测胸径取自然对数后，再建立反演模型的效果最佳，此时拟合优度(R2)达到最高，均方根误差(RMSE)最小。因此，结合本文样地实测数据分析，当LiDAR估测树高范围在14—22 m时，推荐使用第3组方程y=- 2.184 9x2+13.569x-17.605作为胸径的最优估测方程。

3.3 LiDAR估测蓄积量验证

250号小班面积为3.269 6 hm2，在其内部选择2个具有代表性的点位布设样地，并每木检尺，测得该小班平均胸径为25.8 cm，平均树高12.2 m，公顷蓄积量为183.6 m3/hm2。

250号小班经LiDAR飞行后，获得估测平均树高为12.7 m；将该小班单木估测树高代入胸径反演模型y=- 2.184 9x2+13.569x-17.605进行胸径反算，得LiDAR估测平均胸径为26.6 cm。查《江苏省主要树种一元材积表》中“马1”列，根据LiDAR估测的单木胸径统计对应径阶的材积，得估测公顷蓄积量为157.2 m3/hm2。

根据《森林资源规划设计调查技术规程(GB/T 26424-2010)》“13 调查精度”中允许误差的要求，主要小班调查因子允许误差分经营单位性质、小班森林类别等分为A，B，C 3个等级。250号小班位于金坛区薛埠镇东进村金沙湾高尔夫球场内，属金坛区重点林区，因此小班调查因子允许误差适用等级“A”。

经LiDAR估测数据与实测数据对比(如表3)，平均胸径误差为3.1%，平均树高误差为4.1%，公顷蓄积量误差为14.4%，均在《森林资源规划设计调查技术规程(GB/T 26424-2010)》中允许误差范围内(平均胸径、平均树高、每公顷蓄积量误差需控制在5%，5%，15%范围以内)，符合精度要求。

表3 LiDAR估测与实测数据对比

4 结论与讨论

4.1 LiDAR估测受多因素制约

LiDAR技术较传统的人工调查而言可节约大量人力物力，提升工作效率[16]。然而，LiDAR技术在实际应用时易受多因素制约，影响精度。制约因素主要包括：LiDAR采样密度、航带匹配误差和数据处理方法[17]。

在本研究中，选择的250，77号小班均为国外松片林，周围无建筑物遮挡，从单木结构角度来说确保了树冠上部的回波信号，减少了信号的灭失，同时点密度为100个/m2，较好地反映了试验区现状；在航带匹配方面，对LiDAR飞行小班均进行了航带平差处理，保证了定位精度；在数据处理方面，由于试验区地貌均为丘陵，地形趋于平缓，因此在地面点分离方面采用改进的渐进加密三角网滤波算法(IPTD)，将点云分为地面点与植被点，然后分别生成DTM，DSM，CHM和NPC，对点云进行单木分割。总体而言本次LiDAR估测尽量减少了试验误差，估测精度较高。

4.2 LiDAR估测蓄积量可能性探讨

一般而言，林分蓄积量与生物量随平均胸径的增加而增加[18]。研究表明，LiDAR估测树高与实测树高关系密切，而树高也与胸径、蓄积量有较高的相关性[19-20]。

本研究通过建立3组胸径反演模型，结合拟合优度(R2)挑选出优选方程，反算了250号小班的林木胸径。最后查《江苏省主要树种一元材积表》，得公顷蓄积量为157.2 m3/hm2。虽然其结果符合精度要求，但与实测公顷蓄积量183.6 m3/hm2仍有较大差距。分析原因，笔者认为部分区域树冠较为茂密，同时林下植被较为复杂，单木分割时对于树冠边缘难以区分，最终导致获得的单木数量较实际偏少，估测公顷蓄积量较实测数据偏小。

4.3 LiDAR技术应用综合展望

本次采用LiDAR技术对金坛区国外松片林小班进行了估测，并获得了良好的效果，对于森林资源调查领域新技术应用提供了技术参考。

目前对于机载激光雷达用于森林资源调查更偏向于理论研究，实际应用较少。其主要原因在于LiDAR飞行时客观因素制约较多，影响数据的误差，而且飞行成本偏高，若进行区县范围内的全域飞行，将投入大量物力。此外，飞行效果视林分结构而定，若林分结构、树种结构、林下植被较为简单，则LiDAR估测精度较高，若林分结构复杂区域，其精度则会降低[16]。

综上所述，LiDAR技术是未来发展的大趋势，前景应用广阔。在精度符合要求的情况下，可在样地水平成果的基础上大范围应用，节约人力物力，提高工作效率。