基于分水岭优化思想的单木信息分割算法

刘方舟，刘浩，云挺

(南京林业大学，南方现代林业协同创新中心，南京 210037)

单木结构参数信息，如单木位置、树高、冠幅和树木种类等，对于树木生长研究、森林资源调查等具有重要意义[1]。机载激光雷达(light detection and ranging，LiDAR)作为一种主动遥感技术，能够大范围、高精度地获取森林三维点云，刻画森林冠层垂直结构信息，可以大幅降低人工测量的难度且能在较高的精度下保证森林探测的空间完整性和时间一致性[2]。机载激光雷达已在森林中得到了广泛的应用，例如生物量估算、树种分类(结合多光谱数据)、树木结构属性的确定和单木层级信息的获取[3-4]。

基于机载激光雷达的单木分割方法主要分为两类：基于冠层高度模型(canopy height model，CHM)和基于点云数据。在基于CHM的单木算法中，大部分方法可分为3类：1)添加曲线轮廓拟合算法来检测树冠边界[5]；2)将局部表面最大值作为种子点，以识别CHM中的树冠边界[6]；3)使用各种尺寸圆形和方形构成的滑动窗口，以提高估算样地水平树高度局部最大值选择的准确性[7]。其次是图像处理中的一些方法，例如使用图像梯度的变化、二维空间小波和图割方法从CHM中获取单棵树木的各种信息。改进的分水岭(局部极大值算法结合形态学方法[8])和标记控制分水岭(局部极大值指导激光雷达数据的分割[9])方法也都是基于CHM分割的重要组成部分。除此之外，人工智能的算法对研究单木分割提供了一些新的方向[10-12]，例如激光点云被压缩成深度图像，并使用R-CNN(Region-CNN，基于区域的卷积神经网络)模型来进行训练，以学习分割单棵植物的能力[13]。

基于点云数据的方法则是直接对三维点进行大量计算，以达到减少单木层次信息丢失的目的[14]。点特征，如点密度和点分布，对于原始激光雷达数据的单木分割非常重要。均值漂移算法[15]( mean shift algorithm，MSA)及其自适应参数调节的优化算法[16]也已应用于同时分割树冠的垂直和水平结构。K-均值聚类算法的分割精度会受到随机选择的种子点影响，但它仍已被应用于从原始的点云中获取单棵树木的各种信息[17]。在聚类算法基础上优化的点云区域增速长算法也被研究并应用于人工松树林[18]和混合针叶林[19]等不同种类的林地。随着点云数量的增加，计算成本也逐渐增加，将点转换为由点分布确定的体素可以有效地减少纯点云方法的计算量[20]。

尽管人们已经提出了许多方法来检测树冠顶点和分割树冠，但对于高密度林分，树冠相交严重，树冠形状和高度不同的情况，使得个体树冠边缘的分割出现一定问题。因此，需要一种具有鲁棒性、通用性和提供树冠自动精确分割能力的树冠分割算法，作为未来精确林业实践的基础。

本研究旨在针对不同树种及立地条件，采集无人机数据和有人机数据，提出一种新的精确分割树冠边界的方法。基于数字表面模型(digital surface model，DSM)，提出一种分层级并构建能量函数优化树冠边界分割的算法，并与经典的分水岭算法进行比较。本研究将该算法和传统的分水岭算法应用于中国南方高峰林场的桉树和浙江遂昌的白皮松，在不同的种植密度和立地条件下进行研究，通过现场测量得到数据，验证该算法的有效性。

1 材料与方法

1.1 研究区域概况

桉树研究区位于广西亚热带人工林南宁高峰林场(108°7′E，22°49′N)，高峰林场位于北回归线南侧，属亚热带湿润季风气候类型。日照不强，雨量充沛，霜冻少，无雪，气候温和，年均气温21.6 ℃，年均降水量1 304 mm，海拔78～468 m。研究区面积52 km2，地形以缓坡为主，坡度20°～47°。研究区主要树种为桉树(Eucalyptusrobusta)。

白皮松研究区域位于中国浙江新安遂昌(118°41′E，28°13′N)。遂昌属于亚热带季风气候区，温暖湿润，四季分明，雨量充沛，山地垂直气候差异明显。年均气温16.8 ℃，无霜期251 d，年均降水量1 510 mm，地形以陡坡为主。研究区域的主要树种为白皮松(PinusbungeanaZucc.)。

1.2 激光雷达数据

南宁高峰林场信息使用Riegl LMS-Q680i激光扫描仪搭载有人机进行数据获取，数据采集时间为2018年2月，有人机飞行高度750 m，飞行速度为180 km/h，旁向重叠为65%。传感器记录了完整的激光脉冲返回波形，时间间隔为3 ns。扫描仪脉冲发射频率为300 kHz，扫描频率为80 Hz，发射波长1 550 nm，扫描角度为±30°，视场角为60°。光束发散角为0.5 mrad，地面光斑大小为37.5 cm。数据平均点间距为0.45 m，样地平均点密度约为9.58个/m2。最终提取的点云以LAS 1.2格式存储(美国摄影测量与遥感学会，美国马里兰州贝塞斯达)。

浙江遂昌信息使用Velodyne32线激光扫描仪搭载无人机进行数据获取，数据采集时间为2018年3月，无人机飞行高度为80 m，传感器发射波长为905 nm，垂直视场角为40°，水平视场角度为360°，角分辨率(水平/方位角)0.1°～0.4°，旋转频率5～20 Hz；距离精度小于2 cm，扫描频率为10 Hz。为达到对比效果，对无人机数据进行抽稀，使数据密度与南宁的数据保持一致。

1.3 样地实测数据

南宁研究区域内的实测数据于2018年2月进行测量：根据树木密度和地形在数据采集点建立3个正方形样地(边长20 m)，样地主要树种为桉树(1～3号样地)。浙江研究区域内的实测数据于2018年3月进行人工测量，并在数据采集点建立3个正方形样地(边长20 m)，样地主要树种为白皮松(4～6号样地)。

两块研究区域均使用Trimble GEOXH6000全球定位系统(GPS)(Trimble，Sunnyvale，CA，USA)对样地中心进行识别，并使用JSCORS接收到的高精度实时差分信号对这些点进行校正，使其精度为亚米级。在本研究中，有效数据为高度超过2 m的树木，并使用顶点高度测量仪(Haglof, Langsele, Sweden)来测量树冠顶点的高度。树冠冠幅是从树冠顶点位置沿东西、南北两个垂直方向测量长度的平均值。上述样地可分为3类：低密度(样地1、4，树木数量0～25)、中密度(样地2、5，树木数量25～35)和高密度(样地3、6，树木的数量超过35)。

1.4 单木分割算法

1.4.1 传统分水岭算法

分水岭算法是将图像看作地理学的拓扑地貌，可以理解为是一个三维的立体表面。图像每一个像素的灰度值代表海拔，图像中存在许多极小值点。在算法中假设每一个极小值点持续有水流出，并开始逐渐浸没附近的区域，形成一个个集水盆。随着集水盆的水位不断升高，直至水面开始汇聚在一起，为了不让水面汇聚，而直接筑起一道堤坝，也就是分水岭算法的分割线。当所有的图像都别浸没或者筑坝之后，算法结束，被分割出来的集水盆就是算法得出的分割区域。

1.4.2 优化分水岭算法

传统的分水岭算法在众多学者的研究下进行了各种优化，并提出了模糊聚类的分水岭、标记分水岭、形态学分水岭等多种优化方式。但是对于分水岭的筑坝，即为了防止水面汇聚筑起来的堤坝，会对分割区域的边界产生一定的误差，因为分割线在图像边界上依旧会占用像素而造成在密集物体(比如密度较高的林地)使用分水岭算法时，树冠边界精度的损失会显著提升，造成分割的树冠边界比实际的树冠小。针对该问题，本研究重构了分水岭算法的实现方式并进行优化，采用局部最大值算法探查树冠顶点，区域增长算法探寻树冠边界，以及构建能量函数对边界的归属进行判断。

1.5 精度评估

使用优化算法的计算结果与现场测量数据进行比较，以验证优化分水岭算法的准确性。通过以下方程对研究区域内探测到的树冠顶点与观察到的树冠顶点的准确度进行评估。

(1)

(2)

(3)

式中：r为树冠探测率；p为树冠准确率；f为被探测出树的总体准确度；T为被正确探测到树的数量；N为未探测到树的数量(漏分错误)；P为由算法识别但样地中不存在的树。

为检测单个树冠轮廓与人工测量树冠的匹配程度，将实测树冠半径和优化分水岭算法检测的树冠半径进行比较，以验证树冠边界轮廓结果的准确性。本研究选择了正确探测到的树木来比较树冠半径，利用R2、均方根误差(RMSE)和相对均方根误差(RRMSE)评估树冠冠幅的准确性。

2 优化算法的构建

优化算法的流程图见图1。

图1 本研究算法流程图Fig. 1 The flow chart of the proposed method

2.1 DSM生成和树冠顶点探测

首先，研究使用渐进的形态滤波分类器将采集的点云数据分为地面点云和地上点云(树木点云数据)[21]。然后，将地上点云使用栅格化的方法生成对应的DSM，每个栅格(像素)的数值为当前栅格内点云的高度最大值(h)，每个栅格的单元大小为d。不包含数据的栅格由邻近的单元进行最近邻插入填充，其中i、j表示生成栅格集C(即DSM)中的第i行和第j列。

采用具有可变尺寸和参数的高斯核平滑滤波器[22]对生成的DSM进行平滑，为后续探测树冠顶点做准备。在DSM平滑之后，本研究使用可变半径(R)的动窗口来探测平滑之后DSM中每个树冠顶点的位置。公式如下：

(4)

式中：R为滑动窗口的半径；h为中心滑动窗口处的高度值。探测得出的局部极大值被认为是探测所得出的树冠顶点，并作为区域增长算法的种子点。

2.2 分层的区域增长算法

为了准确提取树木的各项参数，研究了分层的区域增长算法以达到从DSM中获取树冠轮廓。将样地分成L1至L5共5层，每一个红色点代表探测出来的树冠顶点，也是区域增长算法中的种子点(图2)。算法按照树冠顶点的高度来决定区域增长的顺序，每一次的区域增长只会在当前分的层级中，在当前层数未被搜索完成时，不会进入下一级进行搜索。

注：a1、b1、c1、d1和f1为激光雷达数据侧视图，显示从最高层到最低层同步向每个树冠增长。a2、b2、c2、d2和 f2分别对应于a1、b1、c1、d1和f1所示运行过程中的DSM。对于LiDAR数据的侧视图和DSM， 扩展过程中采用相同的颜色来表示。黑色线段代表算法每一步中填充的树冠边界。图2 区域增长算法提取树冠边缘Fig. 2 Regional growing algorithm for tree crown boundary delineation

2.3 构建能量函数

a. DSM梯度方向计算图；b. 图a中红色矩形的放大显示；c. 树边界像素与两个最近的树冠之间的高度差；d. 使用最近邻原则判断边界 像素两个最近树冠顶点，绿色线段表示欧几里得距离，淡蓝色表示某一次计算时最外圈的像素。图3 树冠边界归属判定准则说明示意图Fig. 3 Schematic representation of the criteria of the proposed method for canopy boundary determination

用计算机程序不能实现各树冠区域增长完全同步。在同一高度间隔内，计算机只能以随机顺序处理每个树冠的边界归属问题，这种不确定性可能导致相邻树木在相同高度上与树冠相交时出现不准确的树冠分割。因此，需要一个标准来控制每棵树冠随着分层区域增长过程中产生的边界分割问题。本研究提出了一种结合DSM中各像素的梯度方向和高度差信息，构建一个能量函数以优化计算过程中各树冠的边界分割问题(图3)，具体描述如下：

(5)

(6)

(7)

利用高度差来实现树冠的最优分割结果是能量函数的第2部分。两棵相邻的树之间必然存在一个鞍点(最低点)到两棵树冠顶点之间直线上点的高度差一定是最大的。因此，随着区域增长的不断进行，当前搜索的树冠边界上的像素与最近的两个相邻树冠顶点之间的高度差也随之增大。在当前计算的像素高度差达到最大值时，表明该像素很大概率是真实的树冠边界像素。

综上所述，在寻找每个树冠真实边界的搜索过程中，每次迭代计算出的树冠轮廓公式(4)的值会保持下降的趋势，直到达到最小值，这表明程序的搜索达到了树冠的真实边界。反之，如果公式(4)的值在运行过程中出现上升，则代表相交树冠区域会出现错误分割，需要进行回滚(roll back，RB)来调整所有树冠的增长顺序。回滚操作代表当判定现在扩张的树冠是错误的时候，会回到上一次扩张的区域并更换搜索像素开始重新计算。树冠边界归属判定准则说明见图3。

图4 桉树和白皮松分割效果Fig. 4 The segmentation results of pure Eucalyptus tree and Pinus bungeana Zucc. plots

2.4 两种算法应用效果比较

本研究将分水岭算法和优化分水岭算法应用于白皮松、桉树共计6个边长为20 m×20 m的样地。为了验证算法的准确性，本研究展示了密度从小到大的桉树地(1～3号样地)和白皮松样地(4～6号样地)(图4)。图4中，桉树样地的密度(棵/m2)从低到高为0.06(1号样地)、0.09(2号样地)、0.17(3号样地)。白皮松样地的密度(棵/m2)，从低到高分别为0.03(4号样地)、0.07(5号样地)、0.14(6号样地)。优化的分水岭算法平均检测率(r)为0.91，并且6个图的最高精度(f)为0.92。在桉树群(1～3样地)中，2号样地的r值最大(0.97)，3号样地的p值最大(0.89)，f值最大(0.92)。在白皮松(4～6样地)中，5号样地的r值最大(0.90)，f值最大0.88，4号样地p最大值(0.91)。总体的f分数优化的分水岭算法胜于传统的分水岭算法。

在优化的分水岭算法中，桉树冠幅的RMSE范围为0.13～0.23 m，平均值为0.18 m，冠幅的RRMSE为10.83%～13.29%，样地均值为12.04%(表1)。3号样地(RMSE=0.13 m，密度=0.17棵/m2，f=0.92)树冠半径的RRMSE最低(10.83%)。白皮松样地冠幅的RMSE为0.15～0.34 m，平均值为0.26 m，树冠冠幅RRMSE为9.51%～14.05%，平均为12.03%。6号样地(RMSE=0.15 m，样地密度=0.14棵/m2，f=0.87)树冠半径的RRMSE最低(9.51%)。在所有样地中，最高RMSE和RRMSE分别小于0.34 m和14.05%。结果表明，优化的分水岭算法在分割树冠信息方面具有很大的优势。

表1 两种算法分割结果比较Table 1 The comparison of segmentation results

3 结论与讨论

本研究分析了研究区域的树冠顶点检测率。根据样地密度将选择的6个样地分为3类，包括低密度(样地1和4)、中密度(样地2和5)和高密度(样地3和6)，在桉树中，当样地密度较高(0.17棵/m2)时，探测树冠顶点(f)的总体准确度最高，这是由于桉树的种植均匀且有规律。传统分水岭算法在桉树的树冠顶点探测中，密度的改变并未呈现出明显的规律，但样地的总体准确度均低于优化分水岭算法得出的结果。在优化分水岭算法中，白皮松样地密度的增加，被探测的树冠准确率维持一个平稳的趋势，密度对于树冠顶点的探测并没有过大的影响。传统分水岭算法在白皮松样地的探测率也均低于优化分水岭算法的结果。在高密度样地下，白皮松的精度低于桉树，这是由于不同的特性和影响这两个树种的人为因素造成的。桉树是一种速生的经济树种，对土壤肥力具有很强的竞争力，种植密度较高，林下几乎没有其他植物。而白皮松生长环境中不存在过多的人为干预，所以在高密度情况下对桉树的探测率略高于白皮松。中低密度则相反，桉树生产的区域为非经济区域，较少存在人为干预，而白皮松却是典型的针叶林，在局部最大值算法中，比较易于探测树冠顶点。由此可以得出：在优化的分水岭算法下，白皮松树冠顶点探测的准确度对样地密度的敏感性较弱；对于桉树树冠顶点探测的准确度会随着样地密度的增加而上升；高密度下，白皮松树冠顶点探测的准确度低于桉树，中低密度则相反。

本研究还分析了树木冠幅结果的准确性。优化分水岭算法检测的树冠冠幅精度随着两种植物样地密度的增加而提高。在桉树和白皮松中，高密度样地和低密度样地之间的RRMSE之差分别为2.46%(1号和3号样地)和4.54%(4号和6号样地)。这一发现表明，在优化的分水岭算法下，样地密度的增加会影响树木冠幅的正确率。这主要是由于高密度情况下出现极多的树冠相邻的情况，而本研究提出的优化分水岭算法利用分级的区域增长算法保证树冠在同一高度下的探测同步，采用能量函数来进一步保证边界分割时的正确率。在2、5号样地中，两个样地的密度相似，桉树冠幅的RRMSE比白皮松低0.52%。结果表明：在优化的分水岭算法下，样地树木冠幅的探测准确率会随着密度的上升而上升；在密度相同的情况下，桉树树冠冠幅的探测准确率会低于白皮松。

在无人机与有人机数据样本中，无人机原始点云密度比有人机高。高密度的点云会使本研究的算法在生成的DSM有更高的分辨率。分辨率的提高会更加精细单个像素所包含的信息，而本研究提出的优化分水岭算法就是以像素为基本单位进行搜索。研究结果表明，在无人机高密度点云的情况下，对于树冠的分割会将误差控制在更小的单位中进而提高单木分割的效率。

本研究提出的优化分水岭算法使用了分层的区域增长算法、由梯度信息和高度差构建的能量函数作为边界的判定标准来完成单木分割。分层的区域增长算法使每一次搜索像素的过程更加精细，防止数据出现从树木的一面全部搜索完成的情况，并且每一层的搜索会呈现同步的状态，保证了树冠在分割过程中的一致性。构建的能量函数作为每一次边界函数的判断标准，极大地提高了高密度样地树冠与树冠贴合紧密状态下边界像素的归属问题，提高了树木冠幅的分割准确率。本研究的优化分水岭算法可应用于大面积野生自然林以及复杂山地的林业资源调查中。