机载LiDAR数据的林分胸高断面积反演研究

龙 飞,岳彩荣,金 京,李春干,罗洪斌,徐婉婷

(1.西南林业大学 林学院,昆明 650224；2.广西大学林学院,南宁 530004)

森林生态系统在水平和垂直维度上都体现出较强的异质性,不仅在维持全球生物多样性上有重要作用,同时还参与组成全球生物碳循环系统[1-2]。森林生物量则是作为评价森林生产力的重要评价指标,故实现森林生物量快速高效地监测,对全球生物碳循环有着重要的现实意义[3-4]。传统估测森林生物量的主要手段是通过大规模的森林资源调查来实现,但此类调查方式弊端明显,耗财耗力的同时,甚至需要破坏原有的森林生态系统来完成取样,不仅人为产生的误差相对较大,还会因为环境条件限制部分地方无法实地测量[5-6]。机载激光雷达(Airborne Light Laser Detection and Ranging,LiDAR)的出现很好地弥补和解决了这一问题,给未来的森林资源调查提供一种新的选择。

机载LiDAR相较于光学遥感具有较强的穿透性和主动性,精度上也相对较高,可以获取森林的三维结构信息,此类信息在林业方面得到广泛应用[7-8]。胸高断面积(Basal Area,BA)是指树木距地面1.3m处树干的横切面面积,一般是通过测量树干横切面的直径得到,是预测林分材积和森林生长的良好指标[9-11],也是林业工作者在进行森林资源调查时必不可少的林木参数之一[12]。前期,有关机载LiDAR数据在森林特征变量提取方面的研究主要是围绕估测树高、冠幅、林分生物量、材积和碳储量等方面[13-15],用机载LiDAR数据反演林分胸高断面积的研究相对较少。对林木胸高断面积的估测,主要集中在运用不同类型的激光雷达进行多元线性回归模型[16]、点云密度[17]、k近邻算法[18]、均方平均冠层高度[19]等方面,但采用机载LiDAR数据考虑滤波以及随机森林(Random Forest,RF)和迭代的决策树算法(Gradient Boost Regression Tree,GBRT)估测林木胸高断面积的研究较少。

本研究以广西壮族自治区的高峰林场为实验区,以VUX1LR激光雷达系统获取的LiDAR数据和实测的样地数据为基础,比较渐进不规则三角网(Progressive Triangulated Irregular Network,PTIN)、布料模拟滤波算法(Cloth Simulation filter,CSF)、渐进形态学滤波算法(Progressive Morphological Filter,PMF)和基于插值的滤波方法(Interpolation-Based Filtering,IBF)对点云滤波效果的影响,同时采用随机森林和迭代的决策树进行建模,选取模型精度高、适应性好的模型完成森林胸高断面积制图。通过本文的研究,旨在探索一种快速高效获取林木胸高断面积的方法,以期为今后进一步研究提供参考。

1 研究区与数据

1.1 研究区概况

研究区位于广西壮族自治区国营高峰林场,地理位置为22°49′～23°15′N,108°08′～108°53′E,是目前广西最大的国有林场,整体是一个近矩形区域(长11.2km,宽4.2km,面积约为47.04km2)(图1)。实验样地主要以高峰林场界碑、东升、延河3个分场为主,此外,还包括武鸣区、兴宁区、西乡塘区的少量林地[20]。高峰林场:属低山地貌,高程范围在88.76～462.38m,坡度多数在20～35°,海拔由西南向东北逐渐升高；属南亚热带季风气候区,该气候条件下植被多为南亚热带常绿阔叶林；森林覆盖率较好,保持在90%以上；天然林占比较低,约为5%,其余为人工林；按照树木生长周期可将其划分为短周期经营的桉树人工林和长周期经营杉木人工林。由2016年森林资源调查数据可知[20],研究区主要树种为尾叶桉(EucalyptusurohpyllaS.T.Blake)、巨尾桉(Eucalyptusgran-dis×E.urohpylla)、马尾松(PinusmassonianaLamb)、湿地松(PinuselliottiiEngelmann)、杉木(Cunninghamialanceolata(Lamb.)Hook)、八角(llliciumverumHook.f.)等。样地情况如图1所示。

图1 研究区样地分布情况示意图Fig.1 Schematic diagram of the sample plots distribution in the study area

1.2 样地数据获取

样地数据获取于2016年5月至12月。考虑到森林类型、林龄和树种组成等因素对实验准确性的影响,采用手持双频差分GPS接收器记录样地四个角点的地理坐标,用激光测距仪和森林罗盘仪设置30m×30m的方形样地共计105块。对样地内胸径(DBH)≥5cm的树木进行每木检尺并记录树种、胸径、树高、枝下高、冠幅等实验数据,优势树种主要有杉木、马尾松、巨尾桉和一般阔叶树。

1.3 机载LiDAR数据的获取

2016年9月,采用R4直升机搭载集激光测距、全球定位系统和惯性导航系统为一体的VUX1LR激光雷达系统来获取样地点云数据,其中,激光器波长为1 550nm,激光发射角0.5mrad,发射频率820KHz,点云密度约为2.9point/m2。采用CCD相机来获取研究区航空影像,其中,像素为3 000万,空间分辨率为0.2m。R4直升机获取样地时的实际飞行高度约为1 000m,覆盖范围为23°0′28″～22°56′3″N,108°20′9″～108°27′33″E,飞行3次共计15h,覆盖面积约为55km2。

2 研究方法

2.1 样地数据处理

按照样地植物叶的形态特征将森林类型划分为针叶林和阔叶林,并计算林分平均高、样地蓄积和林分胸高断面积。通过样地胸高断面积和平均树高,查阅形高表计算样地的蓄积量。样地数据统计情况如表1所示。

表1 样地数据统计表Tab.1 Plot data statistics table

2.2 机载LiDAR数据的处理

通过去噪、滤波和分类对机载LiDAR数据进行预处理。为提高地面点分类效果,滤波时就渐进不规则三角网(PTIN)、渐进形态学滤波算法(PMF)、布料模拟滤波算法(CSF)和基于插值的滤波方法(IBF)在样地中的滤波效果进行对比,同时结合实验区环境选择可适性较高的方法完成点云滤波。点云滤波方法如下:

1)渐进不规则三角网算法(PTIN)在森林点云数据滤波研究中常被使用。基本思路:首先,根据一定策略选择初始地面点完成初始地形的建立；其次,采用相应准则不断迭代增加地面点；最后,将相关数据叠加至所有的点完成归类。经典的三角网加密方法在坡度陡峭的区域滤波效果不佳,故此本研究使用郭庆华团队提出改进的渐进三角网方法进行滤波[21]。

2)渐进形态学滤波算法(PMF)的理论基础是形态学。该方法主要是利用具有特定形状的结构元素提取图像中的拓扑结构信息,从而用于图像分类,其基本操作为膨胀和腐蚀[22]。

3)布料模拟滤波算法(CSF)是基于“lidR”包实现,主要将激光LiDAR数据进行反转,并用硬布覆盖反转表面。基本流程:首先,通过分析分布节点和相应激光LiDAR数据点之间的相互作用,从而确定分布节点来生成近似地面点；其次,通过对比原始激光雷达点与生成的表面之间的关系来提取地面点；最后,布料模拟本身是由质量和相互连接的粒子组成的网格,这些粒子共同确定布料三维位置和形状[23]。

4)基于插值的滤波方法(IBF)是通过ENVI的BCAL-tool插件来实现。一般流程:首先,依靠一定的规则形成一个初始地表面；其次,采用特定规则对未分类的数据进行判定；最后,不断迭代加入新的地面点,迭代到符合停止加密的规则为止。

点云数据滤波效果评价主要采用Ⅰ类误差(地面点误分为非地面点)、Ⅱ类误差(非地面点误分为地面点)和kp系数为主要指标(表2)[2],利用LiDAR360软件进行人工检查和修改。受到森林结构特征影响,在复杂地形条件或高郁闭度条件情况时,滤波算法很难准确完成地面点分类,因此在滤波处理后,结合目视解译才能得到研究区准确地面点。将机载LiDAR数据分为地面点和非地面点,针对提取的地面点,采用插值的方法生成数字高程模型(Digital Elevation Model,DEM),其分辨率为0.5m。使用DEM对点云进行归一化处理,提取特征变量。

表2 Ⅰ类误差 Ⅱ类误差 总误差和kp系数的计算Tab.2 Calculation of typeⅠ error,typeⅡ error,total error and kp coefficient

2.3 特征变量提取

利用样地矢量数据对预处理后的点云进行裁剪,从裁剪后的样地点云数据中提取34个特征变量(23个高度变量,10个密度变量和1个郁闭度变量)。各点云特征变量描述如表3所示。

表3 机载激光雷达数据中提取的特征变量指标Tab.3 Feature variable indexes extracted from ALS data

2.4 建模方法

通过Python语言将机载LiDAR数据获取的34个变量,利用随机森林算法(RF)和迭代的决策树算法(GBRT)结合样地实测数据分别建模[25-26]。随机森林算法(RF)由Breiman最早提出,是一种集成学习方法。其主要思想是为了让机器学习具有更强的泛适性,从而在将无数的弱机器学习组合在一起形成的合集[27]。迭代的决策树算法(GBRT)是一种由梯度提升(Gradient Boosting)、缩减(Shrinkage)、和回归决策树(Regression Decision Tree)等3个概念组成的新型迭代决策树算法。该算法的最终结果是多棵决策树的结果之和,最早由Jerome Friedman在1999年提出[28]。

通过上述两种算法完成建模后采用独立样本检验的方式来进行模型校验,随机将样本数据的70%用来建模,剩下的30%则作为检验样本。通过十折交叉验证法和网格寻优法对建模参数进行选择,采用决定系数(R2)、均方根误差(RMSE)和相对均方根误差(rRMSE)来对模型进行检验,从而比较随机森林算法和新型迭代决策树算法在机载LiDAR数据反演林分胸高断面积时的适应能力,以便选择最优样本来参与建模。R2,RMSE和rRMSE公式:

(1)

(2)

(3)

3 结果分析

3.1 滤波方法对比

对点云数据进行滤波处理是激光LiDAR数据应用的基础,旨在正确有效地区分地面点和非地面点。滤波效果的好坏将对后期提取的点云特征变量和森林参数估测产生影响。本研究随机选择4块样地,样地信息如表4所示。

表4 样地信息Tab.4 Sample plot information

通过比较PTIN,PMF,CSF和IBF滤波方法在4块样地中的滤波效果,来选择适合具有普适性的方法对研究区点云进行滤波。滤波效果如表5所示,当林地坡度较大、林分郁闭度较低时,IBF滤波效果相较于其他算法较好；当地形平缓,郁闭度较低时,4种滤波方法的精度有一定差异。森林的郁闭度、坡度和点云密度对滤波效果的好坏都会产生影响,选择滤波算法时,应考虑研究区环境特征。研究区大部分坡度在25～33°之间,相对平缓,郁闭度在0.5～0.7之间,雷达获取数据时的点云密度约为2.9 point/m2,样地2与研究区的情况基本一致。结合表5可知,样地2中,CSF滤波算法对该样地滤波效果最好。综上,选用CSF滤波算法对研究区点云数据进行滤波处理可得到最好效果。

表5 4种滤波算法的精度评价Tab.5 Accuracy evaluation of 4 filtering algorithms

3.2 特征选择

从LiDAR数据中提取34个特征变量用于构建RF和GBRT模型。研究表明,采用RF进行变量筛选方法被广泛应用,且效果较好[2,20,30],故采用RF对LiDAR特征变量的重要性进行排序筛选,胸高断面积反演时,选择累计贡献值达到80%的特征变量参与胸高断面积反演。从图2看出,高度变量占所选变量的比例较大,其中Hmax和H99重要性都达到0.17,累计贡献到34%,说明Hmax和H99与胸高断面积具有较强的相关性。

图2 特征变量重要性Fig.2 Importance of characteristic variables

3.3 模型效果分析

构建反演模型时,采用网格寻优法和十折交叉验证法对RF反演参数进行寻优。结果显示:当构造决策树的数量为100,树节点分裂的最小特征数为2,最小叶子节点为1时,模型效果最优。胸高断面积模型的建模精度和验证精度如表6所示。

表6 胸高断面积反演模型建模和检验精度Tab.6 Modeling and test accuracy of inversion model of BA

从表6得出,建模结果中RF和GBRT模型的R2分别为0.77和0.83,RMSE为3.99m2/hm2和3.42m2/hm2,rRMSE为17.76%和15.22%；但在RF和GBRT模型进行独立样本检验的结果中，R2分别为0.66和0.64,RMSE为3.27m2/hm2和3.37m2/hm2,rRMSE为14.73%和15.18%。综合考虑模型建模和检验精度的R2,RMSE和rRMSE之间值的关系,RF模型建模精度和检验精度相对稳定,表明该模型的泛化能力优于GBRT模型,故选用RF模型完成胸高断面积估测。

3.4 胸高断面积反演

综合考虑模型的建模精度、检验精度和模型的泛化能力等条件,最终决定采用RF模型实现研究区林木胸高断面积的反演。胸高断面积反演结果如图3所示。

从图3中看出:胸高断面积在9～12m2/hm2区间的林分呈离散分布；胸高断面积在12～18m2/hm2和18～23m2/hm2区间的林分面积最大,且两者交错分布；胸高断面积在23～30m2/hm2的林分总体面积不大且呈零星块状分布；胸高断面积在>30(m2/hm2)区间的林分最少。上述结果与林地实际情况基本吻合。

图3 胸高断面积反演制图Fig.3 Inversion mapping of BA

4 结论与讨论

基于机载LiDAR数据进行广西高峰林场林分胸高断面积估测研究,在点云数据的滤波处理中,对比了渐进不规则三角网(PTIN)、渐进形态学滤波算法(PMF)、布料模拟滤波算法(CSF)和基于插值的滤波方法(IBF),结果表明,CSF滤波算法效果最好。在滤波基础上提取了34个点云特征变量,利用随机森林算法(RF)和迭代的决策树算法(GBRT)进行林分胸高断面积估测,研究表明,RF模型估测效果优于GBRT模型。在RF模型反演中,Hmax和H99两个特征变量贡献值最大,表明点云最大及99%分位高度变量与林分胸高断面积呈显著相关。这一结论与Silva等[16]利用机载LiDAR数据对火炬松人工林胸高断面积进行估测研究结论有一致性。

相关学者针对激光雷达估测森林胸高断面积的研究表明,胸高断面积不仅受到高度变量影响,还与林龄结构、林分密度、平均树高和平均冠幅等变量有关[29-30]。同时,从测树学角度来看,胸高断面积的大小同林分密度有显著相关性,这与本研究结论有差异,究其原因可以归纳为:1)高峰林场相较于其他学者的研究区在地理环境上存在差异,因此林龄结构、林分密度、平均树高等森林结构因子有一定差异,导致所提取的主要特征变量存在差异；2)由于研究区多数为人工速生丰产林(如桉树林),林分密度无较大差异,故在采用随机森林对胸高断面积进行反演时体现林分密度的变量贡献值较低,而Hmax,H95和H99等3个高度变量贡献率较大。

本研究利用机器学习的方法进行胸高断面积反演,反演方法在模型机理表达方面欠佳,且研究区95%的林木为人工林,天然林范围小,致使本研究在天然林以及大面积混交林上是否适用有待研究,期望在未来的研究中探索更多森林类型以及建模方法(如线性回归和联立方程组等)在胸高断面积反演研究中的应用潜力。

致谢:广西壮族自治区林业勘测设计院,为本文的实验提供了机载 LiDAR 点云数据及野外调查地面数据。在此,表示衷心的感谢！