刘见礼，张志玉，倪文俭，汪垚，张大凤

(1.中国科学院遥感与数字地球研究所 遥感科学国家重点实验室，北京 100101；2. 中国科学院大学，北京 100049)

0 引言

森林是陆地生态系统的主体，也是人类赖以生存的最重要的自然资源之一[1-3]。森林的单木信息在森林的生态研究、资源监测等方面具有重要作用[4-5]。

单木信息提取的前提是单木的准确识别，目前单木识别主要利用机载或星载的高分辨率光学影像数据的光谱与空间特征信息，以及LiDAR点云数据的三维几何信息[6-7]。由于直观的三维几何信息更易于单木识别，因而成为目前的研究热点[8-9]。

常用的三维几何信息通常来自LiDAR点云数据或其衍生的数字表面模型(digital surface model，DSM)、冠层高度模型(canopy hight model，CHM)数据[10-11]。单木识别的方法主要为“局部最大值”法[12-13]。由于LiDAR数据获取成本较高，限制了其推广应用，而无人机影像匹配点云数据的获取成本较低，在单木识别中具有较高的应用潜力[14]。

目前缺少针对无人机影像匹配点云数据的单木识别算法，因此我们借鉴基于LiDAR数据的“局部最大值”算法进行单木识别研究。“局部最大值”算法根据搜索窗口的不同分为两种：一种是“固定窗口”的局部最大值；另一种是“变化窗口”的局部最大值。变化窗口的局部最大值通常需要精确的林下地形信息，而无人机影像点云在郁闭林区难以获得精确的林下地形信息，因而难以采用变化窗口的局部最大值法进行单木识别[15]。

因此使用固定窗口大小的局部最大值法进行无人机影像点云数据的单木识别。已有的研究表明，固定窗口的局部最大值单木识别中，当窗口设置较大时，识别出的单木准确性较高，但会发生遗漏；窗口设置较小时，识别单木的数量会增多，但易引入了错识别单木[16]。目前基于局部最大值的单木识别缺乏针对错分单木的判别算法。发展有效的单木真伪判别算法，解决“小窗口”局部最大值算法产生的错识别单木问题，是提高单木的识别精度的关键。

为此，本文提出了一种联合“局部最大值”与“树冠结构分析”的单木识别算法，该算法可解决 “局部最大值”算法的窗口大小设置的困境。

主要思路：(1)基于无人机匹配点云数据，通过小窗口的“局部最大值”，获得尽可能多的候选单木；(2)对每个候选单木树冠进行连续水平切片分析，获得单木树冠结构曲线，通过单木树冠结构曲线对候选单木进行真伪判别，滤去伪单木；(3)根据树冠切片中心对单木初始位置进行改正，对位置改正后相互聚拢的单木进行合并，最终实现单木准确识别。

1 研究区与数据

1.1 研究区

研究区位于大兴安岭北段西坡，内蒙古呼伦贝尔市东北部，属于根河林业局。研究区气候为寒温带湿润型森林气候，全年平均气温-5 ℃左右，年平均降水量294 mm。区域内优势树种为兴安落叶松(larix gmelini)，伴生树种有樟子松(pinus sylvestris var)、白桦(betula platyphylla)、山杨(populus davidiana dode)等。

1.2 无人机匹配点云数据

无人机影像的获取时间为2015年8月17日至24日，采用大疆六旋翼无人机，搭载Sony NEX-5T相机，影像大小4 912×3 264，飞行高度300 m，影像分辨率8 cm，航向、旁向重叠度分别为90%、60%。

无人机影像处理采用Agisoft Photoscan软件，主要步骤如下：

(1)加载包含可交换图像文件格式(exchangeable image file format，EXIF)信息的无人机影像，若影像不包含EXIF信息，则需要进行摄像机定标。

(2)影像外方位元素计算，首先获取影像上的特征点，然后进行影像间特征点匹配，通过影像间的匹配点使用“8点法”计算得到影像间的外方位元素信息。

(3)无人机影像匹配点云的生成，无人机影像匹配点云是利用上一步获得的影像外方位元素，逐像素进行特征点的提取与同名点的匹配，然后使用光束法平差进行误匹配噪声点判别与去除。

(4)控制点的放置，通过放置控制点的方法，对影像匹配点云进行地理编码。

(5)DSM的生成，为后续的单木识别算法评价，将无人机影像匹配点云进行栅格化，通过插值生成栅格大小为0.1 m的DSM数据。

(6)结果的导出，最终得到的平均点密度约为40个/m2的无人机影像匹配点云数据，栅格大小为0.1 m的DSM数据，数据坐标采用WGS84坐标系，投影为UTM投影北半球6°分带的第51带。

1.3 样地实测数据

实测样地数据获取于2016年8月15日至24日。使用2015年获取的无人机正射影像数据作为样地布设参考，按树种组成、株密度、高度等布设了8个直径为30 m的圆形样地。以正射影像作为底图参考，通过手持GPS导航样地位置，利用单木特征(如树种、枯木、倒木等)对正射影像中的单木进行现场解译，将解译的单木进行编号并标注于正射影像图中。实测的单木位置作为参考数据用于单木识别算法识别的精度评价，为了便于对单木识别算法进行评价，将实测单木中的不被遮挡，能够目视识别的单木定义为“可见单木”，表1为各样地实测信息。

表1 各样地实测信息

2 研究方法

2.1 单木识别

单木识别有两种常见错误：漏识别、错识别。漏识别是指单木因为遮挡或靠近大树等原因无法被识别；错识别分为两种，一种是因为噪声点等原因将伪单木误识别为单木；另一种是由于单木有较大分支或多个顶点而被识别为多株。本文提出的基于无人机影像匹配点云的单木识别算法，通过“较小的搜索窗口”获取尽可能多的候选单木，然后对候选单木进行“树冠结构分析”降低错识别率。具体过程如下：

第1步，使用“局部最大值搜索”算法确定候选单木顶点。设定搜索窗口大小为d1，遍历点云，判断每个点在以其为中心的搜索窗口内是否为最高点，若为最高点则标记为候选单木顶点。

第2步，候选单木冠层点云水平分层切片。在原始点云中选取到候选单木顶点的水平距离小于d2的点云作为待分析点云。d2应大于研究区最大的单木冠幅半径。对选取的单木待分析点云以厚度t自上而下进行水平分层切片，直到没有点云为止，总层数记为n。对各层自上而下按1到n的顺序进行编号(图1)。

图1 单木点云分层切片示例

第3步，逐层进行拟合分析，具体如下：

(1)第1层点云分析。在第1层点云中选取到候选单木顶点距离为小于d3的点进行圆形拟合，得到圆心坐标记为(x1，y1)、半径为r1。

(2)第2至n层点云分析。考虑到树冠通常呈现一定的几何体的形状(如圆锥体，椭球体或球体)，树冠点云会随着层号的增加而逐渐外扩，因此点云的水平选取范围参数d3需要根据层号和树冠形状进行调整，调整方法如下：

d3，i=r1+t×c×(i-1)

(1)

式中：i为层号，取值范围为[2，n]；d3，i为d3在第i层的取值，r1为第一层拟合圆的半径，t为层厚，c为冠幅外扩参数(即树冠表面坡度的余切值)；每层进行圆形拟合得到圆心坐标为(xi，yi)，半径为ri。

第4步，伪单木顶点剔除(图2(a))。由于匹配点云不能穿透冠层，而只能刻画冠层表面的特征，所以对于理想的树冠而言，其冠幅会随着层号的增大而逐步增大，直到达到最大值为止。当某一层拟合圆圆心位置较顶点位置发生明显偏离，偏移距离超过拟合圆半径或冠幅突然变小时，则可认定点云出现异常，记录其上一层号为nc。对于伪单木而言，通常难以出现圆心保持基本不变而冠幅单调递增的规律性现象。因而，可以通过冠层点云总厚度(即层数n)和点云异常层号nc两个参数进行伪单木顶点识别，当n

第5步，临近过识别单木顶点合并(图2(b))。某些树种会由于丛生、伴生或多个一级枝的存在等因素而造成过识别现象。因此在对所有候选单木顶点完成第2步至第4步的处理后，仍有可能会存在过识别的问题。但可以想见，对于由丛生、伴生或多个一级枝引起的过识别候选单木顶点，其拟合圆心会随着层数的增加而逐渐聚拢。可以根据nc层拟合圆心对初始单木位置进行改正，对于位置改正后距离小于d4的单木进行合并。d4为设定的阈值。

图2 真伪单木与过识别单木示例

2.2 完整度参数

前一部分介绍了单木识别的方法，其中的关键在于拟合圆的圆心、半径的应用，考虑到单木冠层水平切片得到的点云会围绕着圆形或圆弧形存在一定的变动，且随着层数的增加圆形或圆弧可能出现缺口，为了提高拟合圆的可信度，增加一个对点云分布圆形完整度的描述参数：完整度指数(completeness index，CI)。将平面按10度的间隔划分为36个区间，如果某区间内有点云落入则取值为1，否则取值为0。36个区间的总和作为完整度指标。因此完整度最低为1，最高为36。当完整度低于18时，意味着点云所占据的圆弧长度低于拟合圆周长的一半，此时得到的拟合圆可信度降低。如前一部分第2步和第3步所述，对每木点云从第2层开始自上而下逐层分析，取2至nc层间最后一层完整度不低于18的层数作为最终的nc层。

3 实验结果与分析

3.1 真伪单木树冠结构展示

本研究提出的单木识别算法涉及到的参数设置：“局部最大值搜索”使用的搜索窗口大小d1=0.5 m，小于当地最小的冠幅值；单木点云选取中用到的点云到候选单木顶点的距离阈值d2=5 m大于当地最大的冠幅半径；根据点云密度及当地树种树冠特点，点云水平切片厚度t=0.3 m，第1层点云到候选单木顶点距离阈值d3=1.5 m，冠幅外扩参数c=1.0，正确单木冠层点云厚度阈值n0=3，其单调递增连续层数最小值nc0=3，单木间最小距离阈值d4=1.0 m。

图3和图4分别展示了真、伪候选单木冠层点云切片结果(图3和图4所展示的单木分别为图5(a)中的编号为NO.21和NO.16单木)。图3(a)与图3(b)为该单木第4层和第7层结果，蓝色点为候选单木中d2范围内的点，红色点为通过d3或d3，i截取的单木点云，拟合圆的半径为“R”，蓝色射线为点云圆环的完整度指数(completeness index，CI)，即为每层圆环的完整度。图3(c)为该单木树冠结构曲线图，可以看到随着层数的增加半径先持续单调增加，直至由于点云不完整导致半径发生突变；图3(d)展示的候选单木的冠幅中心偏离初始候选单木顶点的距离很平稳，随着深度增加，点云的不完整导致发生距离突变。

图4(a)至图4(b)分别对应第1层与第3层。由图4(c)可知，伪单木的冠幅没有随着层数的增加而增加，相反呈现减少的趋势，原因在于该伪单木为噪音点，不具有普通单木树冠的基本特征，因而半径变化无规律；同时，由图4(d)可以看到，伪单木的冠幅中心偏离初始单木顶点的距离随着层数的增加而增加，这不符合真实单木的特点；通过分析可以看出真实单木的冠层结构与伪单木的冠层结构存在明显的差异，因而通过对候选单木冠层结构的分析，能够对伪单木与真实单木进行判别。

图3 真单木树冠结构注：黑圈为点云拟合圆，R为拟合圆半径；蓝色点为单木5 m范围内的点，红色点为本株单木的点；蓝色射线为圆环完整度间隔线；CI为完整度。

图4 伪单木树冠结构注：黑圈为点云拟合的圆，R为拟合圆半径；蓝色点为单木5 m范围内的点，红色点为本株单木的点；蓝色射线为圆环完整度间隔线；CI为完整度。

3.2 单木识别分步结果

由上述方法对大兴安岭根河林区的8个样地的无人机影像匹配点云数据进行处理，获得每个样地的单木识别结果。图5所示为P1样地单木识别的分步结果，底图为P1样地正射影像图。图5(a)将搜索窗口d1设为0.5 m，得到的候选单木顶点(紫色点)，可以看到，由于窗口设置较小，引入了较多的伪单木，如NO.9、NO.16、NO.40、NO.50、NO.53等；以及过识别单木，如NO.14与NO.19、NO.58与NO.59以及NO.67与NO.69等，NO.表示单木的编号。图5(b)是根据候选单木树冠结构分析，剔除掉伪单木后的结果(绿色点)，可以看到图5(a)中错识别的伪单木单木，如NO.9、NO.16、NO.40、NO.50、NO.53等被剔除，但也存在NO.18、NO.49这种小冠幅单木由于冠形特征不明显被误删除的现象。图5(c)对剔除掉伪单木后的单木进行单木位置改正后合并的结果(红色点)，看以可到，图5(b)的错识别中的过识别单木，如NO.11与NO.14、NO.21与NO.22、NO.43与NO.44、NO.51与NO.53等，进行了合并。

图5 单木识别分步结果

3.3 各样地的单木识别结果

为了验证本算法的效果，选取窗口大小分别为1.0 m、2.0 m的“局部最大值”算法，对每个样地的DSM数据进行单木识别，其中当地的最小冠幅直径为1.0 m，DSM栅格大小为0.1 m。图6至图13为8个样地的单木识别结果，从左到右，分别是：实测可见单木(黑色点)、基于DSM数据窗口大小为1.0 m的局部最大值单木识别结果(亮绿色点)、基于DSM数据的窗口大小为2.0 m的局部最大值单木识别结果(紫色点)，本文算法得到的单木识别结果(红色点)，其中蓝色大圆为样地边界。

图6 P1样地结果

图7 P2样地结果

图9 P4样地结果

图10 P5样地结果

图11 P6样地结果

图12 P7样地结果

图13 P8样地结果

3.4 单木识别精度评价标准

单木既可以由单木树冠顶点来表示，也可以由单木树冠边界轮廓来表示，因此存在“点”与“多边形”两种不同形式的结果。同样，参考数据也存在点与多边形两种形式：一般通过GPS或全站仪测量单木树干的坐标，形式为点；由于识别结果与参考结果的形式不同，其匹配的方法也不同。识别与参考结果的匹配主要存在4种方式：“点对点”“点对多边形”“多边形对点”“多边形对多边形”。

由于实测的可见单木位于单木树冠的顶点，本文的3种单木识别算法得到的也是单木的顶点。因此本文采用“点对点”的匹配方式对单木的识别结果进行评价。将样地内的实测可见单木顶点作为参考，将3种算法得到的单木顶点与参考单木顶点进行匹配。匹配的标准是：若1个实测可见单木顶点存在1个识别单木顶点与其对应，则为正确识别，存在0个识别单木顶点则为漏识别，存在2个以上顶点或0个实测单木顶点存在1个以上的识别单木与其对应，则为错识别，具体如表2所示。

表2 单木识别结果类别

本文的单木识别结果使用信息检索与统计学中的准确率、召回率和F测度进行评价，其表达式分别如下：

(2)

(3)

(4)

式中：Pd称为准确率；Pr称为召回率(或查全率)；F称为F测度；Nc为正确识别单木总数；Nr为实测可见单木总数；Nd为识别出的单木总数。准确率Pd的含义为：在识别出来的所有单木中，真实单木所占的比例。Pr的含义为：识别出来的真实单木占总真实单木数的比例。F测度是对准确率和查全率的综合描述，当能正确检出全部真实单木时，F=100%；当检出的单木全部为伪单木时，F=0%；F测度越高代表结果越好。

3.5 单木识别结果的精度评价

表3为8个样地的单木识别结果统计对表，为便于叙述，“1 m窗口大小的局部最大值识别算法”称为“A”算法，“2 m窗口大小的局部最大值识别算法”称为“B”算法，“本文算法”称为“C”算法。

从表3可知，8个样地整体上的A、B、C 3种算法的准确率上A低于B接近20%，而B与C结果接近，说明随着窗口的缩小，错识别单木增加，而由于C成功实现伪单木的剔除，使得识别单木的准确性最高98.61%。在查全率上A高于B接近10%，C高于A接近9%，说明随着窗口的缩小，单木漏识别问题得到抑制，这与前人的研究结论相符。A、B、C 3种方法的F测度分别为74.82%、77.35%与90.45%，即本文算法相比1 m、2 m窗口大小的局部最大值算法分别高近15%与13%，说明了本算法相比传统的单木识别算法，能够提高单木的识别精度。

针对稀疏样地P2、P3、P5与P7样地，除P5样地，3种算法的查全率都比较高，均超过85%，原因在于稀疏样地中单木较为独立，树冠相互遮挡的问题较轻，因而查全率较高，P5样地虽整体较为稀疏，然而单木树冠相连较多，使得有些单木树顶点不明显，导致查全率较低；3种算法的准确率差别较大，其中B与C的准确率都较高，而A算法的准确率较低，原因在于稀疏样地存在较多的背景噪音被误识别为伪单木，因而使得小窗口的A算法的精度偏低，而算法C能够剔除伪单木，所以精度受影响较小。3种算法的F测度与准确率相似，由于小窗口的局部最大值算法，即A，引入了大量的伪单木使得其F测度低于B与C的F测度。

对于密集的样地，P1、P4、P6、P8，3种算法的查全率都不高，都不大于80%，说明了在密集林区局部最大值单木识别算法存在明显的局限性，原因在于密集林区的部分单木由于树冠之间的相互遮挡，使得单木顶点不明显；而3种算法的单木准确性都较高，原因在于密集林区的背景噪声少。P6样地比较特殊，是由于其样地内部单木多伴生，树冠相连，使得单木识别存在较严重的漏识别问题P1与P6样地的内均存在大量的白桦(表1)，然而P1样地的查全率与F测度均高于P6样地，主要原因在于单木的空间分布特征，P1样地虽然白桦较多，但都较为独立，树冠顶点较为明显，因而易识别，而P6样地内树冠相连较多，使得单木树冠顶点不明显，因而降低了单木识别的结果。

表3 各样地单木识别结果统计对比

注：A、B方法分别为1 m、2 m窗口大小的局部最大值法，C方法为本文的方法。

4 结束语

本文利用无人机影像匹配点云数据，进行了单木识别研究，针对单木漏识别与错识别的问题，提出了联合“局部最大值”与“单木树冠结构分析”的单木识别算法。得出以下结论：(1)本算法相比于传统的“局部最大值”算法，能够获得更精准的单木位置信息，为其他单木参数的获取提供准确的初始参数；(2)对于郁闭林区(缺乏林下地形信息)，利用影像匹配点云获取准确的单木位置与株数等信息，拓宽了固定窗口的“局部最大值”单木识别算法的应用范围。

本文的研究还存在以下的问题：(1)由于本算法涉及多个阈值，而阈值的确定以实测数据的统计为基础，这是本算法的不足之处，下一步考虑通过对无人机影像匹配点云的分析，实现阈值的自适应；(2)单木漏识别的问题，由于部分单木没有明确的单木树冠顶点导致局部最大值算法也存在饱和的现象，接下来将考虑综合使用无人机影像点云的光谱与几何信息，进行单木识别研究。