, 习晓环2 1. 东北林业大学 森林作业与环境研究中心, 哈尔滨 150040;
2. 中国科学院遥感与数字地球研究所 数字地球重点实验室, 北京 100094
2018年2月9日 收稿; 2018年4月20日 收修改稿
基金项目: 林业公益性行业科研专项（201504319）和国家自然科学基金（41871264）资助
通信作者: 王成，E-mail:wangcheng@radi.ac.cn
摘要: 针对现有机载激光雷达（light detection and ranging，lidar）数据单木分割算法在密集林区中探测精度较低的问题，结合林木冠层空间结构分层的特点，提出一种从机载点云数据直接分离单木的方法。首先，对原始点云数据进行去噪、滤波、高程归一化；然后基于冠层高度模型（canopy height model，chm）计算局部最大值以确定冠层表面的明显树顶，以此作为单木位置的先验知识，继而采用归一化割（normalized cut，Ncut）方法实现冠层的初始分割；最后，以全局最大值代替局部最大值，并将冠层形状、冠层最小点数作为约束条件，再次利用Ncut方法完成对漏检单木的探测，进而实现单木的精确探测。实验结果表明，针对密集林区的单木分割，本方法有效地减少了漏识单木，整体精度达90%以上，将有助于单木三维结构定量描述及参数反演。
关键词: 激光雷达点云Ncut单木分割
A graph cut-based approach for individual tree detection using airborne LiDAR data
WANG Pu1,2, XING Yanqiu1, WANG Cheng2
, XI Xiaohuan2 1. Center for Forest Operations and Environment Research, Northeast Forest University, Harbin 150040, China;
2. Key Laboratory of Digital Earth Science, Institute of Remote Sensing and Digital Earth, Chinese Academy of Sciences, Beijing 100094, China
Abstract: The accuracy of individual tree segmentation using airborne LiDAR (light detection and ranging) data is generally low for dense forests. A new method, which considers the vertical stratification of forest canopy, is proposed in this study to detect individual tree with high accuracy using airborne LiDAR data. First, several data preprocessing steps are conducted, including noise removal, point cloud filtering, and elevation normalization. Secondly, the initial canopy segmentation is achieved by the normalized cut (Ncut) segmentation with a prior knowledge of individual tree position derived from the local maximum of the canopy height model (CHM). Finally, the Ncut method is used to reduce the leakage rate of individual tree detection by setting the global maximum of CHM as the individual tree position and considering the shape and the minimum point number of the canopy. The results show that the proposed method effectively improves the accuracy of tree detection, which contributes to the quantitative description and parameter inversion in individual tree scale.
Keywords: LiDARpoint cloudnormalized cuttree segmentation
机载激光雷达(light detection and ranging, LiDAR)可以直接、高效地获取高精度、高密度的森林植被三维结构信息，而且其光斑直径较小、全数字化记录回波信息(回波次数、回波强度)[1]，适于单木位置的精准探测与树冠的精确分割，可为森林三维结构定量、精细描述[2](如树高[3]、冠幅、树种等)提供数据源保障，因而也成为当前森林资源调查最有效的手段之一[4]。
国内外****利用机载LiDAR在单木结构参数，特别是垂直结构参数方面提出了多种单木自动探测方法，其中基于冠层高度模型(canopy height model, CHM)进行单木分割是主流方法之一，主要包括局部最大值法、区域增长算法、分水岭算法等。通常情况下可以认为从树顶反射的激光点是单木冠层的最高点，即树顶通常表现为局部最高点。因此，局部最大值算法通常以固定或可变窗口尺寸来检测单木树顶。李响等[5]采用动态窗口局部最大值法对郁闭度较高的针叶林进行单木位置自动提取，精度可达85%。区域增长算法从种子像素开始，以迭代的方式从邻域不断生长，直至满足生长停止条件。甄贞等[6]采用标记控制区域增长方法勾绘树冠边界，但针叶林易出现过度生长情况。分水岭算法是区域增长算法的一种特例，将冠层模型倒置，并将其想象成一个盆地，以局部最小值为入口开始注水，水位上升直至盆地的边界。Chen等[7]提出一种基于冠层极值模型(canopy maxima model, CMM)的树顶探测方法，并基于此应用分水岭算法分割冠层。这些方法存在3点不足之处：1)CHM更多地描述冠层上部，缺少对冠层下部及垂直剖面信息的描述；2)插值方法和格网间距的选择易引入空间误差和不确定性；3)CHM表面平滑程度难以自动控制，可能会引起树高的高估或低估[8-9]。
近年来机载LiDAR系统飞速发展，获取的点云数据密度和精度都极大提高，更多研究趋向于直接利用点云数据进行单木探测[8]，不仅提高了冠层上部的单木探测准确率，而且提升了次生冠层单木探测的可能性。Li等[10]提出基于空间距离的单木分割方法，避免CHM插值误差和不确定性，并在针叶混交林得到较好的应用效果。考虑到树冠的各向异性，Zhang等[9]提出基于滑动膨胀环线圈方法(the donut expanding and sliding algorithm)分割单木，但效率较低且难以探测到下层单木。此外，计算机视觉和图像处理领域的一些经典算法(K-means聚类算法、体素化分割算法等)也被用于激光点云数据的单木分割。Kandare等[11]利用二维和三维K-means聚类算法从点云数据中检测冠层上层、下层树木，但阈值不易控制，易造成过分割。Wang等[12]用体素结构和一种分层的形态学方法，先在同一高度层生成冠层区域，最后从不同高度层的树冠区域中提取单木，但该方法易受激光点密度影响。
前人的方法虽然取得了一定的进展，但在应用于密集林区的重叠冠层分割、中下层单木探测等时效果甚微，而茂密植被区域的单木测量更是林业野外调查的难点。本研究试图从森林冠层结构分层的特点出发，结合归一化割(normalized cut, Ncut)方法提取单木，以提高密集林区的单木探测精度。
1 方法与冠层表面单木不同，密集林区的中下层单木受到抑制，特征不凸显，无法直接对其进行探测。考虑到森林冠层垂直结构这种分层的特点，本研究首先利用局部最大值方法和Ncut方法对冠层激光点云数据探测并分割冠层上部的凸显单木；并在此基础之上，以树冠激光点数阈值和树冠形状指数作为约束条件，全局最大值替代局部最大值，再次利用Ncut算法进行漏检单木探测，识别森林中下层单木；由于缺乏单木位置的地面实测数据，单木探测的精度评价选用目视判读所勾绘的单木位置作为参考，具体处理流程如图 1所示。
Fig. 1
 | Download: JPG larger image |
| 图 1 技术路线流程框图 Fig. 1 Flow chart of single-tree segmentation 图 1 技术路线流程框图 Fig. 1 Flow chart of single-tree segmentation --> | |
1.1 点云数据预处理在LiDAR数据采集过程中，由于各种测量误差和随机误差，原始点云数据中存在噪声点，严重影响数据处理和单木分割效果，因而对点云数据处理前需要去除噪声点。本文采用噪声类型适应性较强的离群点去除方法去除点云噪声[13]。然后，利用具有较强地形适应能力的布料滤波算法[14]进行地面点和非地面点分类，为高程归一化等后续处理做准备工作。最后，为降低地形对单木树高提取精度的影响，采用Lee等[15]方法对植被点云进行归一化处理，每个点的归一化高程值表示该点到地面的高度，因此位于树顶的激光点高程值可认为是此单木树高。
1.2 树顶探测采用局部最大值方法探测冠层表面明显的树顶。首先分别对地面点和非地面点采用反距离加权法(公式(1)和(2))插值生成数字高程模型(digital elevation model, DEM)和数字表面模型(digital surface model, DSM)；然后，基于DSM和DEM生成CHM，并利用文献[16]中方法填充CHM冠层中高程值突变的凹坑；最后，采用窗口为5×5的高斯滤波平滑CHM以减少树顶误判，并采用局部最大值滤波器进行树顶识别。
| $Z = \frac{{\sum\limits_{i = 1}^n {\frac{{{z_i}}}{{{p_i}}}} }}{{\sum\limits_{i = 0}^n {\frac{1}{{{p_i}}}} }}, $ | (1) |
| ${p_i} = {\left( {{{\left( {X - {x_i}} \right)}^2} + {{\left( {Y - {y_i}} \right)}^2}} \right)^{\frac{q}{2}}}.$ | (2) |
在冠层形状规则和中低密集度的林区该方法具有良好的适用性，但应用于树冠形状复杂和高密集度的林区易发生误判。本研究通过设定树顶间的水平距离阈值来降低误判率，即在排除灌木影响的情况下如果两树顶间的水平距离小于设定阈值，则对高程较低的树顶进行剔除。
1.3 Ncut方法分割冠层Ncut算法即基于图论的图像分割方法——归一化割，通过由节点间相似性构建的权值矩阵的特征向量分割图像/点集[17]。Reiberger等[18]和Yao等[19]以树顶、树干的位置和归一化割方法相结合用于激光点云数据单木分割，虽然冠层中下层单木探测率明显提高，但树干识别精度易受激光点云密度的影响。本文直接应用于激光点云，先计算局部最大值以确定冠层表面的明显树顶，以此作为单木位置的先验知识并采用Ncut方法进行冠层初始分割，然后以全局最大值代替局部最大值作为Ncut方法的先验，并作规则约束，完成对漏检单木(主要冠层中下层单木)的进一步探测，进而实现森林单木的精确探测。算法描述如下：
将三维点云映射到无向图G=(V, E)中，V为表示激光点的节点，E为连接每一对节点的边。两节点{i, j}∈V之间的相似性通过权重Wij来描述，即权重Wij表示为两节点{i, j}∈V对应的激光点之间特征相似程度，有
| ${W_{ij}} = \left\{ {\begin{array}{*{20}{l}}{{{\rm{e}}^{ - {{\left( {\frac{{D_{ij}^{XY}}}{{{\sigma _{XY}}}}} \right)}^2}}} \times {{\rm{e}}^{ - {{\left( {\frac{{D_{ij}^Z}}{{{\sigma _Z}}}} \right)}^2}}} \times {{\rm{e}}^{ - {{\left( {\frac{{G_{ij}^{{\rm{max}}}}}{{{\sigma _G}}}} \right)}^2}}}, }&{{\rm{ if }}\left( {D_{ij}^{XY} < {r_{XY}}} \right)}\\{0, }&{{\rm{ otherwise }}}\end{array}} \right.$ | (3) |
Ncut分割通过切除图G=(V, E)中的一些边将图G分割成为两个非连接性点集A、B，且满足A∩B=?和A∪B=V，即类内相似度最大、类间相似度最小。归一化割准则如下
| $\mathit{NCut}(A, B) = \frac{{\mathit{Cut}(A, B)}}{{\mathit{Assoc}(A, V)}} + \frac{{\mathit{Cut}(A, B)}}{{\mathit{Assoc}(B, V)}}, $ | (4) |
求解NCut(A, B)最优解是一个NP完全问题(non-deterministic polynomial complete problem)，本研究通过对矩阵的特征值与特征向量的求解实现近似求解，这也是目前比较常用的方式。
| $(\mathit{\boldsymbol{D}} - \mathit{\boldsymbol{W}})y = \lambda {\mathit{\boldsymbol{D}}_y}, $ | (5) |
1.4 漏检单木的探测本文树顶探测算法可能出现漏检情况(如图 2)，即根据局部最大值算法只能检测到冠层顶部明显的树顶。为降低失误率，提高本文方法的普适性，将全局最大值代替局部最大值作为Ncut方法的先验知识，并进行规则约束，对未分类点云中可能存在的单木或树梢进行迭代检测，方法描述如下：
Fig. 2
 | Download: JPG larger image |
| 图 2 多株单木相邻较近 Fig. 2 Example of multiple adjacent trees 图 2 多株单木相邻较近 Fig. 2 Example of multiple adjacent trees --> | |
输入：Pts-未分类的林区冠层点云数据
输出：TreePts-单木冠层激光点云
1) 定义Pts激光点集为U，做其高程值降序处理；
2) 选取最高点作为种子点，定义为pi(pi∈U)，统计半径为dXY的领域点数N，若N小于阈值TD，则将其归属到其最近单木中；
3) 利用1.3中方法进行冠层分割得到种子点的单木冠层点云Ti，此时，U=U-Ti；
4) 利用文献[20]的方法计算Ti的形状指数(shape index, SI)，若判断Ti为树梢，则其属于其周围单木的一部分，否则，TreePts=TreePts+Ti；
5) 重复2)、3)、4)，直至Ui为空。
2 结果与分析利用美国华盛顿州西部的Capitol State森林蓝岭地区(The Blue Ridge area)的机载点云数据验证本文方法的有效性。研究区以丘陵为主，海拔260~425 m，地面坡度0°~45°；主要树种以针叶林为主，如道格拉斯冷杉、西部铁杉和西部红雪松等。机载点云数据由直升机搭载TopoEyeAB的TopEye系统获取，时间为1999年春季，平均点密度约为4.86 pts/m2。本文选取3块100 m×100 m的样地作为研究样地，标记为Plot 1、Plot 2、Plot 3(图 3)。
Fig. 3
 | Download: JPG larger image |
| 图 3 研究区高分辨率影像 Fig. 3 High-resolution image of the study area 图 3 研究区高分辨率影像 Fig. 3 High-resolution image of the study area --> | |
首先，对原始数据(图 4(a))采用1.1中预处理方法进行去噪、滤波，得到地面点和非地面点(图 4(b))。经滤波之后，分别对非地面激光点云和地面激光点云利用反距离加权插值方法生成DSM(图 5(a))和DEM(图 5(b))。为尽量减少树冠边缘信息的损失，降低树顶漏检率，DSM和DEM分辨率d依据研究区域点云数据的平均点密度ρ进行设置，
Fig. 4
 | Download: JPG larger image |
| 图 4 原始点云与滤波结果 Fig. 4 Original point cloud and filter result 图 4 原始点云与滤波结果 Fig. 4 Original point cloud and filter result --> | |
Fig. 5
 | Download: JPG larger image |
| 图 5 点云生成的DSM(a)和DEM(b) Fig. 5 DSM (a) and DEM (b) generated using LiDAR data 图 5 点云生成的DSM(a)和DEM(b) Fig. 5 DSM (a) and DEM (b) generated using LiDAR data --> | |
然后，采用1.2所描述方法对研究区域滤波后的点云进行树顶初步探测，提取出位于冠层上部的明显树顶。由于植被叶片间的空隙及遮挡性等其他因素影响，CHM影像上易出现凹坑。凹坑的存在，使得CHM不能正确表达冠层表面的结构形态，同时可能导致树高的低估，因此很有必要去除CHM局部区域凹坑。本文利用文献[16]中方法，目视确定最佳阈值，进行凹坑去除。用窗口为5×5的局部最大值方法进行探测树顶，为排除可能存在的干扰值，降低树顶的错误率，引入树顶间的水平欧式距离阈值；该阈值若设置过小，则树顶错误率增加；若设置过大，则树顶漏检率增加。研究区域的实地调查，可作为该阈值调整的依据，本文该阈值设为4 m。图 6分别为3个样方CHM影像和树顶的初步探测结果。
Fig. 6
 | Download: JPG larger image |
| 图 6 局部最大值法探测3个区域的树顶结果 Fig. 6 Treetop detection results using the local maximum method 图 6 局部最大值法探测3个区域的树顶结果 Fig. 6 Treetop detection results using the local maximum method --> | |
将以上树顶探测结果作为1.3中算法运行的先验知识，进行迭代式分割处理，并将探测到的树顶数作为迭代次数上限，每次分割出一棵单木。由于冠层的遮挡性和结构分层的特点，利用1.2中方法只能探测到冠层表面主要单木的树顶，以此为先验知识的1.3中的冠层分割方法也只能分割出冠层表面主要单木的冠层，而位于冠层中下层的单木存在漏检情况。对于漏检单木，在未分类点集中采用1.4中方法以激光点高程全局最大值代替局部高程最大值，作为Ncut分割方法的先验条件，从高到低依次进一步探测，最终冠层分割结果如图 7所示。
Fig. 7
 | Download: JPG larger image |
| 图 7 3个区域的单木冠层分割结果 Fig. 7 Segmentation results of three plots 图 7 3个区域的单木冠层分割结果 Fig. 7 Segmentation results of three plots --> | |
为评价本文算法精度，通过对激光点云和高空间分辨率航空影像目视判读方式得到本研究单木位置精度的验证数据，采用召回率r、正确率p、综合考虑r和p的调和值F [10]来评价。召回率表示有效单木探测株数占真实参考数据的比例，正确率表示有效单木探测株数占整个提取结果的比率。单木分割算法将点云或影像分割成多个部分，每个部分可以看作一棵单木，分割的结果大致分为3种情况:正确分割(TP)、欠分割(FN)和过分割(FP)。因此，r，p，F可描述如下：
| $\begin{array}{l}r = \frac{{{T_{\rm{P}}}}}{{{T_{\rm{P}}} + {F_{\rm{N}}}}}, \\p = \frac{{{T_{\rm{P}}}}}{{{T_{\rm{P}}} + {F_{\rm{P}}}}}, \\F = \frac{{2rp}}{{r + p}}.\end{array}$ | (6) |
Table 1
 表 1 本文方法与分水岭算法的对比Table 1 Comparison between the proposed method and the watershed algorithm
| 表 1 本文方法与分水岭算法的对比Table 1 Comparison between the proposed method and the watershed algorithm | |||||||||||||||||||||||||||||||||||||||||||||||||||||||||
3 结束语本研究结合密集林区冠层结构分层特点提出基于图割的单木分割方法，实现冠层凸显单木的分离，试验结果表明本文方法在高株密度林区具有较高的单木探测精度。尽管如此，但当本文方法应用于高密度点云数据时，会增加Ncut的求解负担，降低算法效率，自适应尺度的体素化高密度点云或许是一种解决办法。另外，本研究只选择了针叶林区域进行算法测试和分析，需要进一步尝试算法对冠层结构更加复杂的阔叶林或针阔混交林的适应性。这也是目前相关研究工作者所面临的一个难题，利用分层分割的混合方法或新型手段(如：点云数据+高光谱影像/多光谱影像等)来解决此问题将是未来的研究方向。
参考文献
| [1] | 李亮, 王成, 李世华, 等. 基于机载LiDAR数据的建筑屋顶点云提取方法[J]. 中国科学院大学学报, 2016, 33(4): 537-541. |
| [2] | 李增元, 刘清旺, 庞勇. 激光雷达森林参数反演研究进展[J]. 遥感学报, 2016, 20(5): 1138-1150. |
| [3] | 邢艳秋, 尤号田, 霍达, 等. 小光斑激光雷达数据估测森林树高研究进展[J]. 世界林业研究, 2014, 27(2): 29-34. |
| [4] | 李旺, 牛铮, 王成, 等. 机载LiDAR数据估算样地和单木尺度森林地上生物量[J]. 遥感学报, 2015, 19(4): 669-679. |
| [5] | 李响, 甄贞, 赵颖慧. 基于局域最大值法单木位置探测的适宜模型研究[J]. 北京林业大学学报, 2015, 37(3): 27-33. |
| [6] | 甄贞, 李响, 修思玉, 等. 基于标记控制区域生长法的单木树冠提取[J]. 东北林业大学学报, 2016, 44(10): 22-29. DOI:10.3969/j.issn.1000-5382.2016.10.005 |
| [7] | Chen Q, Baldocchi D, Gong P, et al. Isolating individual trees in a Savanna Woodland using small footprint LiDAR data[J]. Photogrammetric Engineering & Remote Sensing, 2006, 72(8): 923-932. |
| [8] | Zhen Z, Quackenbush L, Zhang L. Trends in automatic individual tree crown detection and delineation:evolution of LiDAR data[J]. Remote Sensing, 2016, 8(4): 333. DOI:10.3390/rs8040333 |
| [9] | Zhang C, Zhou Y, Qiu F. Individual tree segmentation from LiDAR point clouds for urban forest inventory[J]. Remote Sensing, 2015, 7(6): 7892-7913. DOI:10.3390/rs70607892 |
| [10] | Li W, Guo Q, Jakubowski M K, et al. A new method for segmenting individual trees from the LiDAR point cloud[J]. Photogrammetric Engineering & Remote Sensing, 2012, 78(1): 75-84. |
| [11] | Kandare K, Dalponte M, Gianelle D, et al. A new procedure for identifying single trees in understory layer using discrete LiDAR data[C]//Geoscience and Remote Sensing Symposium. IEEE, 2014: 1357-1360. |
| [12] | Wang Y, Weinacker H, Koch B. A LiDAR point cloud based procedure for vertical canopy structure analysis and 3d single tree modelling in forest[J]. Sensors, 2008, 8(6): 3938-3951. DOI:10.3390/s8063938 |
| [13] | 朱俊锋, 胡翔云, 张祖勋, 等. 多尺度点云噪声检测的密度分析法[J]. 测绘学报, 2015, 44(3): 282-291. |
| [14] | Zhang W, Qi J, Wan P, et al. An easy-to-use airborne LiDAR data filtering method based on cloth simulation[J]. Remote Sensing, 2016, 8(6): 501. DOI:10.3390/rs8060501 |
| [15] | Lee H, Slatton K C, Roth B E, et al. Adaptive clustering of airborne LiDAR data to segment individual tree crowns in managed pine forests[J]. International Journal of Remote Sensing, 2010, 31(1): 117-139. DOI:10.1080/01431160902882561 |
| [16] | Benarie J R, Hay G J, Powers R P, et al. Development of a pit filling algorithm for LiDAR canopy height models[J]. Computers & Geosciences, 2009, 35(9): 1940-1949. |
| [17] | Shi J, Malik J. Normalised cuts and image segmentation[J]. IEEE Transactions Pattern Analysis & Machine Intelligence, 2000, 22(8): 888-905. |
| [18] | Reitberger J, Schn?rr C, Krzystek P, et al. 3D segmentation of single trees exploiting full waveform LiDAR data[J]. ISPRS Journal of Photogrammetry & Remote Sensing, 2009, 64(6): 561-574. |
| [19] | Yao W, Krzystek P, Heurich M. Enhanced detection of 3D individual trees in forested areas using airborne full-waveform LiDAR data by combining normalized cuts with spatial density clustering[J]. Biopolymers, 2013, 349-354. |
| [20] | Guo Q, Kelly M, Gong P, et al. An object-based classification approach in mapping tree mortality using high spatial resolution imagery[J]. Mapping Sciences & Remote Sensing, 2007, 44(1): 24-47. |
