机载激光雷达(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所示。

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)

式中：(X, Y, Z)为插值点坐标，(x_i, y_i, z_i)为搜索样本点坐标，p为权重，q为幂次，i为样本点坐标编号，n为搜索区域内样本点总数。

在冠层形状规则和中低密集度的林区该方法具有良好的适用性，但应用于树冠形状复杂和高密集度的林区易发生误判。本研究通过设定树顶间的水平距离阈值来降低误判率，即在排除灌木影响的情况下如果两树顶间的水平距离小于设定阈值，则对高程较低的树顶进行剔除。

1.3 Ncut方法分割冠层

Ncut算法即基于图论的图像分割方法——归一化割，通过由节点间相似性构建的权值矩阵的特征向量分割图像/点集^[17]。Reiberger等^[18]和Yao等^[19]以树顶、树干的位置和归一化割方法相结合用于激光点云数据单木分割，虽然冠层中下层单木探测率明显提高，但树干识别精度易受激光点云密度的影响。本文直接应用于激光点云，先计算局部最大值以确定冠层表面的明显树顶，以此作为单木位置的先验知识并采用Ncut方法进行冠层初始分割，然后以全局最大值代替局部最大值作为Ncut方法的先验，并作规则约束，完成对漏检单木(主要冠层中下层单木)的进一步探测，进而实现森林单木的精确探测。算法描述如下：

将三维点云映射到无向图G=(V, E)中，V为表示激光点的节点，E为连接每一对节点的边。两节点{i, j}∈V之间的相似性通过权重W_ij来描述，即权重W_ij表示为两节点{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)

式中: D_ij^XY表示节点i与节点j间的水平距离，D_ij^Z表示它们的垂直距离，G_ij^max=max(D^XY(i, treeTop), D^XY(j, treeTop))，D^XY(i, treeTop)、D^XY(j, treeTop))分别表示节点i、j与最近树顶之间的水平距离；σ_XY、σ_Z、σ_G分别表示D_ij^XY、D_ij^Z、G_ij^max的标准差；r_XY表示节点间相似性的水平方向作用距离，水平距离超过该距离的节点之间无相似性。

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)

式中：$\mathit{Cut}(A, B) = \sum\limits_{i \in A, j \in B} {{W_{ij}}} $，表示A中节点和B中节点之间的权值之和，即A和B之间的相似性；$ Assoc(A, V) = \sum\limits_{i \in A, j \in V} {{W_{ij}}} $，表示A和所有节点之间的权值之和；$ Assoc(B, V) = \sum\limits_{i \in B, j \in V} {{W_{ij}}} $，表示B和所有节点之间的权值之和。

求解NCut(A, B)最优解是一个NP完全问题(non-deterministic polynomial complete problem)，本研究通过对矩阵的特征值与特征向量的求解实现近似求解，这也是目前比较常用的方式。

$ (\mathit{\boldsymbol{D}} - \mathit{\boldsymbol{W}})y = \lambda {\mathit{\boldsymbol{D}}_y}, $

(5)

式中:令n=|V|，权重矩阵W为n×n对称矩阵，W(i, j)=W_ij，{i, j}∈V；D为n×n对角矩阵，且$ D(i, i) = \sum\nolimits_j {W(i, j)} $。由Rayleigh商的基本性质可知，求解NCut(A, B)的最小值问题转化为求解该特征系统的第二最小特征向量，即可利用该特征向量对图进行划分。但由于第二最小特征向量的元素一般呈现为连续实数值，因此需要引入一个分离点来二分，通常用0或特征向量元素的中值作为分离点。为使NCut(A, B)最小，即图G得到最优化分，采用试探法寻找最优分离点^[18]。

1.4 漏检单木的探测

本文树顶探测算法可能出现漏检情况(如图 2)，即根据局部最大值算法只能检测到冠层顶部明显的树顶。为降低失误率，提高本文方法的普适性，将全局最大值代替局部最大值作为Ncut方法的先验知识，并进行规则约束，对未分类点云中可能存在的单木或树梢进行迭代检测，方法描述如下：

输入：Pts-未分类的林区冠层点云数据

输出：TreePts-单木冠层激光点云

1) 定义Pts激光点集为U，做其高程值降序处理；

2) 选取最高点作为种子点，定义为p_i(p_i∈U)，统计半径为d_XY的领域点数N，若N小于阈值TD，则将其归属到其最近单木中；

3) 利用1.3中方法进行冠层分割得到种子点的单木冠层点云T_i，此时，U=U-T_i；

4) 利用文献[20]的方法计算T_i的形状指数(shape index, SI)，若判断T_i为树梢，则其属于其周围单木的一部分，否则，TreePts=TreePts+T_i；

5) 重复2)、3)、4)，直至U_i为空。

2 结果与分析

利用美国华盛顿州西部的Capitol State森林蓝岭地区(The Blue Ridge area)的机载点云数据验证本文方法的有效性。研究区以丘陵为主，海拔260~425 m，地面坡度0°~45°；主要树种以针叶林为主，如道格拉斯冷杉、西部铁杉和西部红雪松等。机载点云数据由直升机搭载TopoEyeAB的TopEye系统获取，时间为1999年春季，平均点密度约为4.86 pts/m²。本文选取3块100 m×100 m的样地作为研究样地，标记为Plot 1、Plot 2、Plot 3(图 3)。

首先，对原始数据(图 4(a))采用1.1中预处理方法进行去噪、滤波，得到地面点和非地面点(图 4(b))。经滤波之后，分别对非地面激光点云和地面激光点云利用反距离加权插值方法生成DSM(图 5(a))和DEM(图 5(b))。为尽量减少树冠边缘信息的损失，降低树顶漏检率，DSM和DEM分辨率d依据研究区域点云数据的平均点密度ρ进行设置，$ d=\sqrt{1 / \rho} $，本研究d近似值0.5 m。

然后，采用1.2所描述方法对研究区域滤波后的点云进行树顶初步探测，提取出位于冠层上部的明显树顶。由于植被叶片间的空隙及遮挡性等其他因素影响，CHM影像上易出现凹坑。凹坑的存在，使得CHM不能正确表达冠层表面的结构形态，同时可能导致树高的低估，因此很有必要去除CHM局部区域凹坑。本文利用文献[16]中方法，目视确定最佳阈值，进行凹坑去除。用窗口为5×5的局部最大值方法进行探测树顶，为排除可能存在的干扰值，降低树顶的错误率，引入树顶间的水平欧式距离阈值；该阈值若设置过小，则树顶错误率增加；若设置过大，则树顶漏检率增加。研究区域的实地调查，可作为该阈值调整的依据，本文该阈值设为4 m。图 6分别为3个样方CHM影像和树顶的初步探测结果。

将以上树顶探测结果作为1.3中算法运行的先验知识，进行迭代式分割处理，并将探测到的树顶数作为迭代次数上限，每次分割出一棵单木。由于冠层的遮挡性和结构分层的特点，利用1.2中方法只能探测到冠层表面主要单木的树顶，以此为先验知识的1.3中的冠层分割方法也只能分割出冠层表面主要单木的冠层，而位于冠层中下层的单木存在漏检情况。对于漏检单木，在未分类点集中采用1.4中方法以激光点高程全局最大值代替局部高程最大值，作为Ncut分割方法的先验条件，从高到低依次进一步探测，最终冠层分割结果如图 7所示。

为评价本文算法精度，通过对激光点云和高空间分辨率航空影像目视判读方式得到本研究单木位置精度的验证数据，采用召回率r、正确率p、综合考虑r和p的调和值F ^[10]来评价。召回率表示有效单木探测株数占真实参考数据的比例，正确率表示有效单木探测株数占整个提取结果的比率。单木分割算法将点云或影像分割成多个部分，每个部分可以看作一棵单木，分割的结果大致分为3种情况:正确分割(T_P)、欠分割(F_N)和过分割(F_P)。因此，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)

利用上述评价标准，本研究方法对漏检单木处理前和处理后的分割结果分别进行评价，并与经典的分水岭方法比较，如表 1所示。因本文测试数据Plot 3的株密度较小，不参与算法结果的定量比较，但仍然可以定性地说明本文方法在该株密度条件下具有较好的适用性。已知Plot 1的株密度大于Plot 2，对比可知，基于CHM的分水岭算法在林区环境比较简单的Plot 2研究区的结果与本文方法相当，但在株密度较高的Plot 1研究区域中的分割精度远低于本文方法，有近27%的单木未探测到。继续分析发现，本文研究方法在对漏检单木处理前的分割结果与基于CHM的分水岭算法相当，且在株密度较高的Plot 1研究区域中分割精度较差，原因是本文树顶探测方法是基于CHM影像，缺少冠层垂直结构信息，因此未能有效地探测到冠层中下层单木。根据表 1，本文考虑到林区冠层结构分层特点的单木探测方法，在茂密林区取得非常好的单木检测率，并能有效地降低单木漏检率。

3 结束语

本研究结合密集林区冠层结构分层特点提出基于图割的单木分割方法，实现冠层凸显单木的分离，试验结果表明本文方法在高株密度林区具有较高的单木探测精度。尽管如此，但当本文方法应用于高密度点云数据时，会增加Ncut的求解负担，降低算法效率，自适应尺度的体素化高密度点云或许是一种解决办法。另外，本研究只选择了针叶林区域进行算法测试和分析，需要进一步尝试算法对冠层结构更加复杂的阔叶林或针阔混交林的适应性。这也是目前相关研究工作者所面临的一个难题，利用分层分割的混合方法或新型手段(如：点云数据+高光谱影像/多光谱影像等)来解决此问题将是未来的研究方向。