高郁闭度人工林无人机激光雷达单木分割方法优化*

朱泊东 罗洪斌 金 京 岳彩荣

(西南林业大学林学院 昆明 650224)

森林覆盖地球表面的1/3，是世界上一半以上陆地生物的家园，也是全球75%淡水的来源，在碳循环、水循环、生物多样性保护和可持续发展中发挥着至关重要的作用(Siepielskietal.， 2017)。单木是森林的基本结构单元，其三维结构、生长状况、空间分布等特性是森林资源调查、监管和生态建模研究等所需的关键因子(李增元等， 2016)，准确的单木分割是获取精细化森林参数的前提条件(Meietal.， 2004； Korpelaetal.， 2007； Zhaoetal.， 2018； Nuijtenetal.， 2019)。

传统森林资源调查多以人工为主，存在投入大、成本高、效率低、时空机动性差等明显缺点。激光雷达(light detection and ranging， LiDAR)是一种快速获取目标地物三维空间信息的主动遥感技术，其通过发射激光脉冲到树木的枝叶等组分上，接收包含三维空间位置等信息的散射回波信号，并可生成点云数据，对精确估算冠层三维结构信息尤其在林分尺度和单木尺度的森林垂直结构探测中具有显著优势(Pedersenetal.， 2012； 刘鲁霞等， 2014)。目前，基于机载LiDAR数据进行单木分割的方法大致可分为2类。一是冠层高度模型(canopy height model， CHM)单木法，利用2D栅格化CHM提取单木特征(Brandtbergetal.， 2003； Brandtberg， 2007； Kochetal.， 2006； Popescu，2007； Yuetal.， 2011； 刘清旺等， 2010)，其基本思路是将点云分为地面点和植被点，采用空间插值方法生成相同空间分辨率的数字表面模型(digital surface model， DSM)和数字高程模型(digital elevation model， DEM)，二者相减得到CHM，通过分析CHM局部最大值判断树冠顶点，再依据树冠顶点提取树冠边界(Popescuetal.， 2003； Johanetal.， 2003； Maltamoetal.， 2004)，或先提取树冠边界，再将树冠内的高度最大值作为树冠顶点(Brandtbergetal.， 2003； Zimbleetal.， 2003)。常见的CHM分割方法包括分水岭分割(watershed segmentation)(Tangetal.， 2007)、标记控制分水岭(marker controlled watershed)(Chenetal.， 2006； 郭昱杉等， 2016)、区域增长法分割(regional growth segmentation)(Hyyppä， 1999)等，该类方法的缺点是仅利用冠层表面和地面信息，对于复杂林分的亚冠层、中下层乔木和灌木存在探测和提取困难问题，且对点云数据的三维信息挖掘深度不够，点云栅格化处理过程中会造成点云空间信息丢失，分割精度也受空间插值方法和栅格分辨率大小影响(霍朗宁等， 2021)。二是归一化点云(normalized point cloud， NPC)单木法，利用NPC数据，根据点云之间的空间位置关系设定判别准则对点云进行聚类，从而实现单木分割并提取单木特征(Alexander， 2009； Reitbergeretal.， 2009； Lietal.， 2012)。常见的NPC分割方法包括k均值聚类(Zimbleetal.， 2003)、均值漂移聚类(Ferrazetal.， 2012)、自适应距离聚类(Leeetal.， 2010)等，该类方法的缺点是在高郁闭度林分条件下对中下层林木提取能力较弱，分割精度受点云密度、郁闭度以及林分中每株树水平分布和垂直结构的影响，表现为林分密度或郁闭度较大时分割精度较低(Lietal.， 2012)，郁闭度较小时分割精度较高(Falkowskietal.， 2008； Kaartinenetal.， 2012)。Forzieri等(2009)研究表明，在株数密度较高的林分内，现有单木分割方法对冠层结构复杂和郁闭度高的林分分割效果均不佳。Vauhkonen等(2011)采用6种单木分割方法对针叶林、阔叶林和混交林进行分割，结果发现不同方法在不同林分条件下的分割效果差异较大。Yang等(2019)探讨植被特征(森林类型、叶面积指数、冠层覆盖度、株数密度和高度变异系数)对单木分割(CHM分割、无坑CHM分割、基于点云分割和基于层堆叠种子点分割)精度的影响，得出植被特征差异会对分割效果产生较大影响，多数植被特征均与分割精度呈负相关，但其机制尚不明确。以往研究大多基于单一分割方法对一种或多种复杂森林场景进行单木分割，视森林参数变化采用不同分割方法的单木分割报道较少； 且多数研究虽有对分割方法进行参数优化和方法改进使之更适用于某一林分类型，但不同分割方法在不同林分条件下的分割效果各有优劣，某一林分条件下采用与之对应的分割方法能够进一步提升分割精度； 同时，据目前报道，对高郁闭度林分条件下亚热带人工针叶林内同一树种在不同龄组下单木分割精度差异的相关研究仍较少。

无人机LiDAR技术是获取精细化森林参数的重要技术手段，利用LiDAR点云数据进行单木分割是提取精细化森林参数的基础，森林参数获取精确与否取决于单木分割是否准确。针对高郁闭度林分条件下基于LiDAR点云数据单木分割林木提取困难、总体精度较低等问题，本研究以云南思茅松(Pinuskesiyavar.langbianensis)人工林为研究对象，以无人机LiDAR点云数据为基础，开展多种点云分割方法试验，提出一种基于冠层起伏率结合分水岭算法和基于点云的局部最大值聚类算法的分层分割法，以期为开展无人机LiDAR技术森林资源调查提供技术参考，为提高单木分割总体精度提供新策略。

1 研究区概况与数据

1.1 研究区概况

万掌山林场位于云南省普洱市思茅区北部，思茅区地处云南省南部、普洱市中南部、澜沧江中下游，北邻无量山，东连江城县，22°27′—23°06′N、100°19′—101°27′E。思茅区属低纬高原南亚热带季风气候区，境内立体气候明显，有北热带、南亚热带、中亚热带和北亚热带4种气候类型，具有低纬、高温、多雨、湿润、静风的特点，年均气温17.8 ℃，年均降雨量1 524.4 mm，无霜期318 天，冬季严寒，夏季酷暑，四季温和(欧光龙等， 2014)。万掌山林场占地面积1.95万 hm2，其中包含生态公益林4 777 hm2、商品林5 866 hm2(潘瑶等， 2016)，以思茅松人工纯林为主。研究区位置及样地分布如图1所示。

1.2 遥感数据

以纵横大鹏CW-30油电混合中型垂直起降固定翼无人机为遥感平台，平台搭载CW-30 LiDAR离散回波形激光雷达传感器和JOUAV CA-102 4200万像素全画幅航摄相机，用于获取LiDAR数据和DOM影像。无人机飞行高度900 m，航宽180 m，旁向重叠度≥60%。激光雷达脉冲发射频率820 kHz，点云平均密度19.21 pts·m-2，DOM影像空间分辨率0.1 m，DEM空间分辨率1 m。机载LiDAR数据获取时间为2020年1月5日。

1.3 目视解译样地数据

基于无人机高分辨率影像，结合激光点云的三维形态并实地踏查辅助视觉验证确定目视解译的单木树冠中心位置，以此作为验证数据。在研究区内设置35块目视解译样地(图1)，包含幼龄林、中龄林和近熟林3个龄组，其中幼龄林样地13块、中龄林样地14块、近熟林样地9块，共1 427株树，样地大小均为20 m×20 m。本研究采用的目视解译结果仅包含上层林木。

2 研究方法

将高郁闭度思茅松林分按龄组划分为幼龄林、中龄林和近熟林，采用分水岭算法、基于点云的局部最大值聚类算法进行单木分割，对比2种方法的分割效果，找到不同林分条件下的最佳分割方法和适用条件； 在此基础上提出了基于冠层起伏率结合分水岭算法和基于点云的局部最大值聚类算法的分层分割法，并分析分水岭算法中CHM空间分辨率和DSM插值方法对分割效果的影响。ULS数据单木分割方法对比及精度分析技术路线如图 2 所示。

图2 ULS数据单木分割方法对比及精度分析技术路线Fig.2 Comparison of individual tree segmentation methods for ULS data and technical route of accuracy analysis

2.1 基于点云数据的CHM生成

在对LiDAR点云数据进行去噪、地面点分类(滤波)的基础上，用地面点插值生成DEM，用非地面点插值生成DSM，试验中CHM由不同空间分辨率下采用不同插值方法(克里金插值、反距离权重插值、不规则三角网插值)生成的DSM与同一DEM作差得到。克里金(Kriging)插值是一种基于统计学的插值方法，其计算优化的协方差，利用高斯过程进行插值，插值的整个过程相当于对未知区域化变量的期望值进行加权滑动平均，并将该值作为栅格单元值。反距离权重(inverse distance weighting， IDW)插值基于相近相似原理，利用附近点计算栅格单元值，并以点距栅格单元中心点距离幂次方的倒数为权重进行加权平均，距栅格单元中心点距离越近的点被赋予的权重越大。不规则三角网(triangulated irregular network， TIN)插值基于地面点数据的三角形不规则网络推导每个三角形的双变量函数，并将其用于估计未采样位置(已知地面点之间)的值，每个生成的三角形平面用于插值。

2.2 单木分割方法

2.2.1 分水岭算法 分水岭算法(watershed algorithm)是一种基于拓扑理论的数学形态学分割方法(Meyeretal.， 1990； Vincentetal.， 1991)，其基本思想是把图像看作测地学上的拓扑地貌，图像中每一点像素的灰度值表示地形起伏高低，对灰度值取反； 把图像中一些局部极小值点看作盆地底部，在每个盆地底部钻一个洞，此时地下涌出的水不断上升并逐渐淹没整个地形； 当2个相邻盆地中的水即将汇合到一起时，在水交汇的边缘上修筑一道堤坝阻止盆地间水的交汇，该堤坝即对应的分割线，当众多堤坝将每盆地都包围时，堤坝的集合就是分水岭，即图像分割边界(韩悬等， 2019)。分水岭算法对微弱边缘具有较强的响应，由于图像质量不断提高，图像纹理特征逐渐增多，图像中的噪声、物体表面细微灰度变化均有可能产生过分割现象，因此栅格分辨率大小对图像分割效果具有较大影响(孙钊等， 2021)。本研究采用分水岭算法先确定每株树的树冠边界范围，再在该范围中寻找局部最大值，将局部最大值看作树冠顶部(Wulderetal.， 2000)。

2.2.2 基于点云的局部最大值聚类算法 基于点云的局部最大值聚类算法假定树木之间存在一定间距，通过分析点的高程以及该点到其他点的距离，自上而下地对点云进行聚类； 将归一化点云中的全局最高点作为最高木顶点，基于间距阈值和最小间距判别准则进行迭代判断，设置一个阈值，将离目标单木顶点水平距离小于阈值的点云作为目标单木的点云集群，反之则移除该集群，最终实现单木分割(Lietal.， 2012)。

2.2.3 基于冠层起伏率结合分水岭算法和基于点云的局部最大值聚类算法的分层分割法 高郁闭度林分条件下，龄组差异对单木分割精度具有一定影响，考虑到实际应用中大范围、精确的龄组信息获取难度较大，本研究采用相关分析方法挖掘反映龄组差异的激光雷达点云特征。首先基于激光点云数据提取59个点云特征，包括高度变量(46个)、点云强度变量(10个)、郁闭度、叶面积指数和间隙率，然后分析点云特征与分割精度的相关性，结果发现，在59个点云特征中，冠层起伏率(canopy relief ratio， CRR)与分割精度显著相关，故本研究以冠层起伏率差异作为分割方法选择依据。冠层起伏率(Parkeretal.， 2004)计算公式如下：

(1)

式中： mean为每一统计单元内所有点的平均高度； min为每一统计单元内所有点的最小高度； max为每一统计单元内所有点的最大高度。

对机载LiDAR点云数据进行去噪、滤波、点云分类和归一化，提取归一化点云的冠层起伏率，采用四分位数法将冠层起伏率由低到高分为4级，每一级应用对应分割方法进行单木分割。

2.3 精度验证

以探测率、准确率和F得分为指标对单木分割精度进行验证和评价(Goutteetal.， 2005； Sokolovaetal.， 2006)，计算公式如下：

(2)

(3)

(4)

式中： TP表示被正确探测到的树的数量(正确分割)； FN表示未被探测到的树的数量(欠分割)； FP表示被算法识别到但实际不存在的树的数量(过分割)；r为单木探测率，p为被探测到单木的正确率，r和p可以近似看作制图精度和用户精度；F得分为综合考虑r与p计算得到的加权调和平均值。

3 结果与分析

3.1 基于龄组的单木分割效果

图3所示为分水岭算法和基于点云的局部最大值聚类算法在不同龄组样地的分割效果及其精度验证，表1所示为3个龄组样地的基本概况及2种分割方法在各龄组对应的单木分割精度。

图3 2种分割方法在不同龄组样地的分割效果及其精度验证Fig.3 Results and accuracy of two methods of individual tree segmentation for the forest with different age groups

由表1可知，分水岭算法的F得分由高到低分别为幼龄林(0.78)、中龄林(0.68)和近熟林(0.67)，基于点云的局部最大值聚类算法的F得分由高到低分别为近熟林(0.77)、中龄林(0.66)和幼龄林(0.64)，分水岭算法和基于点云的局部最大值聚类算法在各龄组对应的单木分割精度不同且存在较大差异，但2种方法在分割精度的变化趋势上刚好相反，具有互补性。

表1 3个龄组样地的基本概况及2种分割方法在各龄组对应的单木分割精度Tab.1 The basic profile of the sample plots of the 3 age groups and the segmentation accuracy of individual trees corresponding to the 2 segmentation algorithms in each age group

不同年龄阶段，思茅松树冠形状不同，冠层表面起伏程度不同，林分内部结构也存在较大差异。在幼龄林中，分水岭算法相较基于点云的局部最大值聚类算法有着更高的单木探测率和准确率，这是因为幼龄林大多为高度差异较小、株数密度较高的单层林，冠层点云几乎处于同一水平面，通过探测局部最大值点云的基于点云的局部最大值聚类算法容易造成欠分割，导致分割精度大大降低，而基于拓扑理论的分水岭2D形态学分割方法通过对CHM进行滤波增强图像纹理有助于提升分割精度，但当林分郁闭度过高时2种分割方法均无法得到较好的分割效果，错分割和欠分割现象非常严重。在中龄林中，分水岭算法与基于点云的局部最大值聚类算法的F得分差异不明显，分水岭算法相较基于点云的局部最大值聚类算法单木探测率更高，但准确率较低，基于点云的局部最大值聚类算法与之相反。在近熟林中，基于点云的局部最大值聚类算法相较分水岭算法对上层乔木有着更高的探测率和准确率，这是因为近熟林大多为高度差异较大的复层林，树冠起伏明显，但树木之间冠层交织严重，且许多思茅松近熟林树枝分叉较多，分水岭2D形态学分割方法利用CHM进行局部最值探测时易将树枝识别为树冠导致过分割，而基于点云的局部最大值聚类算法(3D分割方法)利用点云空间分布对点云进行迭代聚类，能够有效降低将树枝错分为树冠的概率。

3.2 基于冠层起伏率结合分水岭算法和基于点云的局部最大值聚类算法的分层分割法

从图4可以看出，高郁闭度思茅松幼龄林点云大多分布在冠层，导致样地内点云平均高度mean上移，由式(1)可知幼龄林的冠层起伏率最大； 中龄林和近熟林因其垂直结构复杂，点云在垂直方向上分布较为均匀，不会导致mean过大，近熟林的冠层起伏率略高于中龄林。

图4 不同龄组样地点云剖面Fig.4 Point cloud profile of different age groups sample

为更好揭示不同龄组样地间冠层起伏率的分布规律，构建小提琴图(图5)。小提琴图基本上是核密度图以镜像方式在箱线图上的叠加，白点表示中位数，白色盒型范围为下四分位点到上四分位点，细白线表示须，外部形状即核密度估计。图5显示，思茅松幼龄林和近熟林的冠层起伏率近似正态分布，中龄林的冠层起伏率近似双峰偏态分布，对比三者分布情况，幼龄林冠层起伏率的中位数、上四分位数和下四分位数均高于中龄林和近熟林； 核密度估计结果亦显示，幼龄林与中龄林和近熟林存在显著差异； 近熟林冠层起伏率的下四分位数、中位数、上四分位数均高于中龄林，说明随着年龄增加林分内点云的垂直结构分布发生改变。

图5 3个龄组样地冠层起伏率的小提琴图Fig.5 Violin plot of canopy relief ratio of three age groups in sample plot

图4、5分析表明，冠层起伏率与龄组之间存在较好的对应关系，以龄组作为分层分割参数可以有效提高分割精度，但在实际应用中，一些区域大范围精细化的龄组信息获取难度较大，且因立地条件差异常导致相同龄组的林分分割效果也存在较大差异，为此，本研究提出对冠层起伏率进行分级代替龄组作为分层分割的重要依据。根据试验结果，采用四分位数法可将冠层起伏率由低到高分为4级： 1级对应中龄林，2级对应近熟林，3～4级对应幼龄林。冠层起伏率1～2级样地采用基于点云的局部最大值聚类算法进行单木分割，冠层起伏率3～4级样地采用分水岭算法进行单木分割，幼龄林与中龄林、近熟林划分的正确率分别为0.67和0.94，幼龄林划分结果包含12块幼龄林样地、2块中龄林样地和4块近熟林样地； 中龄林、近熟林龄组划分结果包括1块幼龄林样地、11块中龄林样地和5块近熟林样地。

为评价基于冠层起伏率结合分水岭算法和基于点云的局部最大值聚类算法的分层分割法的分割效果，本研究对比不分层单一分割方法(分水岭算法和基于点云的局部最大值聚类算法)的分割精度(表2)，基于冠层起伏率结合分水岭算法和基于点云的局部最大值聚类算法的分层分割法的分割精度(0.75)相较分水岭算法(0.71)或基于点云的局部最大值聚类算法(0.68)有明显提升。

表2 单木分割精度验证Tab.2 Accuracy assessments of individual tree segmentation

3.3 CHM空间分辨率和插值方法对分割精度的影响

分水岭算法基于CHM进行图像分割，为评价不同方法生成的CHM对分割精度的影响，采用反距离权重、克里金和不规则三角网3种常用空间插值方法得到的DSM生成CHM，选取0.2 m×0.2 m、0.5 m×0.5 m、0.8 m×0.8 m和1.0 m×1.0 m 4种CHM空间分辨率进行单木分割，并对比其分割精度(图6)。

在3种常用空间插值方法中，当CHM空间分辨率为0.5 m×0.5 m时，反距离权重插值得到DSM生成的CHM分割精度最高(0.81)； 当CHM空间分辨率为0.2 m×0.2 m时，分割精度由高到低依次为反距离权重(0.44)、不规则三角网(0.42)和克里金(0.40)； 当CHM空间分辨率为0.5 m×0.5 m时，分割精度由高到低依次为反距离权重(0.81)、不规则三角网(0.80)和克里金(0.77)； 而当CHM空间分辨率为0.8 m×0.8 m和1.0 m×1.0 m时，不规则三角网插值得到的分割精度最高，由高到低依次是不规则三角网、反距离权重和克里金，0.8 m×0.8 m时F分别为0.59、0.47和0.47，1.0 m×1.0 m时F分别为0.46、0.37和0.37(表3)。

表3 分水岭算法中CHM空间分辨率和插值方法对分割精度的影响Tab.3 Influence of CHM resolution and interpolation method on segmentation accuracy in watershed algorithm

图6 3种空间插值方法在不同分辨率下的分割精度Fig.6 Segmentation accuracy of three spatial interpolation methods at different resolutions

4 讨论

基于高郁闭度思茅松林分无人机激光雷达数据提取59个点云特征，分析点云特征与分割精度的相关性发现，在59个点云特征中，除冠层起伏率外，高度变异系数和累计高度百分位数等垂直结构特征与龄组也表现出一定相关性，但冠层起伏率对龄组的区分度优于其他特征参数。研究曾对样地冠层起伏率尝试多种分级方式，如相等间隔法、自然断点分级法、分位数法和标准差法等，结果表明采用分位数法将冠层起伏率分为4级与样地龄组的一致性最好。

分水岭算法探测准确率较低，过分割现象较多； 基于点云的局部最大值聚类算法探测率较低，欠分割现象较严重，该结果与李平昊等(2018)的研究结论一致，表明分水岭算法和基于点云的局部最大值聚类算法单木分割适用的林分条件不同，并表现出一定互补性，这也是本研究将2种算法相结合的主要依据。

研究发现，当林分郁闭度为1时，冠层起伏率变化出现2种情况： 一是当点云最小高度为0时，冠层起伏率偏大，统计单元内有少数激光点云透过树冠间间隙打到地面上，大部分点云集中在树冠上部，统计单元内所有点云平均高度偏大，导致点云平均高度与最小高度的差变大，即式(1)分子变大，而统计单元内最大高度和最小高度的差不变，即式(1)分母不变，致使冠层起伏率出现异常增大； 二是当统计单元内没有激光点云打到地面上时，冠层起伏率变化规律与龄组无关，此时统计单元内点云平均高度和最小高度均变大，即式(1)分子和分母同时变大，冠层起伏率可能变大、变小或不变，其变化趋势主要取决于统计单元内最低点高度。此外，当郁闭度为1时，上述2种分割方法无法准确进行单木探测和提取，故不讨论其分层分割精度。

本研究主要对高郁闭度林分上层林木进行探测和分割，未涉及中、下层林木，且受点云数据限制，仅针对思茅松林分开展试验，其分割优化方法是否适用于其他森林类型有待进一步研究。分割精度方面，如果能将机载CCD相机获取的高分辨DOM影像与LiDAR数据相结合，在高密度点云基础上加入高空间分辨率的光谱信息或纹理特征，则有望提升其单木分割精度。此外，应引入更先进的单木分割方法，如深度学习等。

因缺少足够数量具有精确地理坐标的实测单木检验样本，本研究单木分割精度检验主要通过无人机DOM影像目视判读进行，这一定程度上影响了评价结果的可靠性。林内精确定位单木位置坐标一直是森林调查中的瓶颈，如果能解决高郁闭度林分林下单木定位精度低的问题，将对准确评价单木分割精度、改进分类方法具有较好促进作用。

5 结论

针对高郁闭度林分条件下基于LiDAR点云数据单木分割林木提取困难、总体精度较低等问题，本研究比较分水岭算法和基于点云的局部最大值聚类算法的单木分割精度，提出基于冠层起伏率结合分水岭算法和基于点云的局部最大值聚类算法的分层分割法。通过提取样地冠层起伏率，确定分水岭算法和基于点云的局部最大值聚类算法单木分割适用的林分条件，可拓宽单一单木分割方法在不同林分条件下的优势，提升单木分割精度。此外，还通过改变栅格化冠层高度模型(CHM)的空间分辨率和数字表面模型(DSM)的空间插值方法，对分水岭算法单木分割精度敏感性进行分析。结果表明：

1) 在幼龄林中，冠层起伏率较大，分水岭算法对单层林的分割效果优于基于点云的局部最大值聚类算法； 在中龄林和近熟林中，冠层起伏率较小，分水岭算法易将思茅松树枝识别为树冠，基于点云的局部最大值聚类算法的分割效果优于分水岭算法；

2) 基于冠层起伏率结合分水岭算法和基于点云的局部最大值聚类算法的分层分割法充分考虑不同龄组的林分结构差异，精度最高(F= 0.75)，优于分水岭算法(F= 0.71)和基于点云的局部最大值聚类算法(F= 0.68)；

3) 在分水岭算法中，当分辨率为0.5 m×0.5 m时采用反距离权重法(IDW)插值得到的CHM单木分割精度最高(r= 0.70，p= 0.94，F= 0.81)，高于克里金插值(F=0.80)和不规则三角网插值(F=0.77)。