基于SfM的针叶林无人机影像树冠分割算法

杨全月 董泽宇 马振宇 吴 悠 崔 琪 卢 昊

(1.北京农学院计算机与信息工程学院， 北京 102206； 2.北京林业大学信息学院， 北京 100083；3.弗莱堡大学环境与自然资源学院，弗莱堡 79106)

0 引言

树冠是树木获取光能并进行能量转换的主要场所，冠幅反映了树木的生长活力及林木空间属性[1]，准确获取树冠信息有助于监测树木长势、估算树木生物量及小班蓄积量、预防树木病虫害[2]。利用传统的光学遥感手段，可进行面向对象的图像分割树冠提取，一般使用光谱和纹理信息[3]、图像梯度[1]或图像分水岭[4-5]等方法。随着无人机(Unmanned aerial vehicle，UAV)技术的迅速发展，无人机航拍影像广泛应用于高精度的森林资源调查等领域[6]。利用无人机影像可快速获取林区高分辨率影像数据，运用传统遥感图像处理手段提取单木冠幅信息[7]。一般认为，激光雷达(Light detection and ranging，LiDAR)点云具有很高的几何精度和细节反映能力，可精确提取单木树高、冠幅、树冠体积、林隙等信息[8-11]，基于机载LiDAR数据的单木分割算法中，具有代表性的分别为基于冠层高度模型(Canopy height model，CHM)[12-13]的分割算法和基于原始点云的分割算法[14]。

随着摄影测量技术和计算机视觉技术的发展，基于双目视觉原理，可对图像构建立体像对、进行前方交会，求得林下树木几何结构信息[15-17]。由于计算性能的大幅度提高，传统摄影测量的空三加密发展为密集匹配算法[18]，出现了快速有效的运动恢复结构(Structure from motion，SfM)技术。通过密集匹配点云提取林分参数已取得了初步进展，有效地解决了高分辨率图像的过分割问题。曾健等[19]运用倾斜摄影测量技术提取了落叶树人工林地形，发现无人机影像匹配点云与机载LiDAR点云具有很高的一致性。LI等[20]对有人机(Charge-coupled device，CCD)航片进行密集匹配，有效反演了东北地区地上生物量和叶面积指数。陈崇成等[21]利用密集匹配点云CHM模型分割了苗圃单木。刘见礼等[22]提出一种点云分层切片方法，实现了单木树冠提取。许子乾等[23]结合LiDAR数字高程模型(Digital elevation model，DEM)研究亚热带森林林分特征变量，发现Lorey’s高等林分指标与UAV点云指标具有较高敏感性。

分析这些方法发现，无人机拍摄一般在人工林[19]、稀疏林地[22]、苗圃[21]等郁闭度较低的林分中才能直接观测地表，得到DEM；在密林[23]、天然林[18, 20]中，影像密集匹配点云往往需要借助LiDAR等主动遥感手段获取的DEM进行高度归一化，再进行林木参数提取。这种依赖性显然制约了无人机影像在大范围天然林区自动提取森林参数的实用性。本文利用无人机快速机动航拍的优势获取地势陡峭、坡度较大的生长季针叶林高重叠度影像，基于SfM技术进行三维表面重建，借鉴激光雷达点云分水岭分割单木的思想，提出一种无需DEM数据支持的自适应邻域分水岭算法，旨在解决在茂密林区高分辨率无人机影像三维重建无法得到地形信息以及树冠过分割的问题，从而提升无人机影像在森林资源调查中的实用性。

1 实验区概况与数据获取

1.1 实验区概况

本研究区位于北京市西部百花山国家级自然保护区，地处北京市门头沟区清水镇境内，是以保护暖温带华北石质山地次生落叶阔叶林生态系统为主的自然保护区。森林植被垂直分异明显，海拔为1 000～1 900 m。选择区内一处北向坡面为实验林地，海拔为1 100～1 200 m，面积约为8 000 m2，位置分布如图1所示。林地坡度约为40°，地势陡峭，生长植被主要为华北落叶松(Larixprincipis-rupprechtiiMayr)。

图1 实验区位置Fig.1 Overage of experiment region

1.2 无人机数据获取

数据获取时间为2018年9月，实验区树木处于生长季末期，部分林冠叶片颜色开始发黄，区分度较大。使用大疆精灵4型(DJI Phantom 4)无人机，搭载FC330型相机，内置Sony Exmor R CMOS传感器，有效像素为1 200万，感光单元尺寸为1.58 μm，焦距为3.61 mm，因此，在相对飞行高度100 m时地面采样分辨率(Ground sample distance，GSD)约4.4 cm。由于目标区域坡度较大，坡面上下高差超50 m，为保证对不同高度处的林木均有近似分辨率成像，依据山势设计了5种不同高度(1 235、1 244、1 265、1 276、1 318 m)航线进行重叠覆盖。无人机作业采用自动巡航模式，使用DJI GS Pro型地面站进行控制。通过指定成像范围和影像重叠度的方式自动计算飞行轨迹和曝光点，以定点悬停拍摄的方式完成区域成像。通过云台控制镜头始终垂直向下以保证成像姿态角的相对稳定。无人机内置全球定位系统(Global positioning system，GPS)模块记录曝光点位置并以(Exchangeable image file format，EXIF)信息的形式存储在图像中。最终5个不同高度共获得高分辨率无人机图像509幅，实验林地坡面区域内的影像曝光点如图1所示。

2 原理与方法

本文方法总体思路是从二维影像构建三维模型，再从三维模型回归林木二维结构信息，技术流程如图2所示。

图2 总体技术流程图Fig.2 Overview technical flowchart

2.1 无人机影像SfM三维重建

2.1.1密集影像SfM处理

为提取林地三维表面结构，首先对原始无人机影像进行SfM三维重建。SfM算法步骤为：密集影像同名特征提取和匹配。利用同名特征点进行影像相对定向并计算外方位元素。空三处理生成稀疏点云及迭代光束法平差。基于稀疏点云进行加密得到密集点云。本研究使用Pix4Dmapper软件处理无人机影像。特征点搜索使用结合特征匹配和像素匹配的Daisy算子，无人机内置GPS提供的曝光时刻机身坐标作为影像外方位元素的坐标初始值进行迭代优化，最后使用面片多视角立体视觉(Patch-based multi-view stereo，PMVS)算法搜索图像中的像素以得到更多的匹配点，加密点云数据。SfM技术进行同名点搜索、匹配和点云加密的方式如图3所示，图3a下方背景为SfM生成彩色点云，空中蓝色圆点为平差后像主点位置，蓝色矩形表示相机投影方向，一组绿色射线表示地表一个同名点被多个影像同时观测，前方交会于地物空间。图3b表示该同名点在一系列不同影像中的搜索位置。初始相机内参数修正误差为2.48%，联合光束法平差的平均误差为0.2像元，符合测绘要求。最终得到实验区点云密度为600～1 000个/m2。同时，对原始无人机影像进行拼接处理和地理编码，平均采样分辨率为2.78 cm，利用点云高程信息进行了正射纠正并生成空间分辨率为0.05 m的彩色影像作为实验区底图(图1)。

图3 SfM技术密集匹配彩色点云示意图Fig.3 Illustrations of densely matched true color point cloud from SfM

2.1.2三维表面模型处理

基于3D点云分割的单木提取方法依赖树冠的点云包络或较完整的点云冠层模型[14,22]，但由于研究区生长季树叶茂密，林地郁闭度较高，无人机影像难以拍摄到林冠内部细节、林窗及林下地表，或无法在多张影像上构成同名点观测。由图3a可知，绝大部分点集中在树冠顶层，导致上述3D点云分割方法适用性较差。但由于冠层表面三维点云密度极高，可利用表面模型处理方法进行树冠提取。针对该特点，使用点云数字表面模型(Digital surface model，DSM)内插方法生成分辨率0.05 m的DSM模型，并进行了表面平滑处理以消除表面噪声。

2.2 kNN自适应邻域分水岭分割

2.2.1单木检测

单木检测是为了找出树木中心位置，作为后续树冠分割的起始位置。为了充分利用无人机影像高分辨率纹理和SfM三维结构信息，设计了一种结合图像信息和高程信息的检测策略。

(1)基于图像的单木检测

首先，进行图像预处理以增强树冠边缘和边界差异。利用Sobel算子和Laplace算子增强图像中树冠的边缘，使树与树之间、树与地面之间的边界更加清晰明确。利用中值滤波可以在一定程度上去除增强时出现的噪声点。在预处理过程中，两种算子和中值滤波交替使用。形态学开运算和闭运算可以使树冠内部紧密、体积缩小，并且树冠之间的间距增大，直方图均衡化使图像的树冠部分和非树冠部分的灰度差距增大。进行上述处理后，使用滑动窗口的大津阈值法将图像分为多个小块，对每一块使用大津阈值法得到二值图像。再对该二值图像使用像素之间的欧氏距离变换计算初始树冠中心[24]，最后使用图像分水岭法分割出树冠。由于图像进行了形态学开运算和闭运算，树冠图像不能准确代表真实树冠，因此仅保留树木位置信息。

(2)基于高程信息的单木检测

由于树木本身有向上生长的特点，在地形高差变化较大的区域，树冠仍然符合树顶高于其他部分的特点。基于高度信息进行树冠检测，一般利用树冠几何形态特征作为先验知识，树顶位置往往在高程上处于局部最大值点。参照文献[22,25]，设定一个最小树冠半径R，对于表面模型栅格中的每一个点，判断与距离小于R的所有点是否全部低于该点，将符合条件的像素视为初始树顶点。分析发现，一般针叶林林分都较为均匀，定值R足够进行检测，自适应调整的影响不大。实际上，文献中一般使用CHM进行检测，但在坡度较大情况下，树木冠幅内部不同像素对应的高程可能不一样，导致CHM树冠形状发生扭曲，直接应用DSM可在很大程度避免该问题。

图4 自适应kNN邻域示意图Fig.4 Illustrations of adaptive kNN neighborhood

(3)单木树冠检测结果合并

树冠最高点一般是较为准确的树顶点，而图像识别的中心点不一定是准确的树顶点，依据高程信息检测的结果往往更加准确，但图像信息可以在三维重建不精确的位置补充识别部分漏检单木。由于三维模型的准确度有限且针叶树木的特殊冠型(树冠凹凸不平滑)，需设定较大的树冠半径R才能较好地避免过分割，但导致有些树木不能被识别。对未被识别的树，使用图像信息作为补充，二者合并的具体策略是：对每个高程检出点，去除其周围半径为R圆周内的图像检出点；再把剩余图像检出的树顶点加入到高程检出的树顶点，形成完整的树顶点检测结果。

2.2.2自适应kNN邻域检测

(1)

其中dii=0(i=1,2,…,m)。待分割的目标树为T，它的kNN邻域定义为

(2)

(3)

2.2.3邻域分水岭分割

图5 单木检测结果Fig.5 Results of individual tree detection

分水岭分割是一种常用于图像分割的数学形态学方法，基本思想是将输入的灰度图像视为拓扑地貌，基于流溢模型进行分水岭变换。对于图像分割而言，通常在梯度图像上进行分水岭变换，而对于表示冠层高度的栅格模型，分水岭可直接应用，不同之处在于树冠高度模型经过“倒置”，一般由最高点开始“浸水”。由于已检测了单木初始位置，本方法是一种标记控制的分水岭算法，目的是基于kNN(T)和Ttree求得Ctree集合。Ttree为树顶检测得到的顶点集合，表示分水岭算法的初始状态。在每个kNN邻域内，从k+1个单木初始点开始进行分水岭处理，设kNN邻域的高程范围是[Hmin,Hmax]，随着水位从Hmin到Hmax不断增加，邻域内的地形会被水漫过，每个被漫过像素将被分配给不同树冠，当这k+1个树冠两两之间出现连通时，在连通像素上标记阻断。最终，邻域内存在一些点属于两个或多个树冠，这些点构成“分水线”，进而确定了树冠像素集合，如图4c所示。由于邻域是围绕目标树构建的，只有包围目标树冠的分水线最接近自然边界，故只取邻域分水岭结果中目标树的分割区域作为该邻域分割结果，其他邻接树冠予以舍弃，如图4d所示。最后将所有单木作为目标树的分割结果合并成完整的林分树冠分割结果。综上可知，邻域分水岭与传统分水岭最大区别在于邻域内高程范围的不同。传统分水岭作用于全局CHM模型，范围一般为从0到最大树高，而邻域分水岭由于不需高程归一化，随着不同目标树邻域在地形表面移动，自适应地获得较小的高程范围，避免高程范围过大提取错误的分水线。

3 结果验证与分析

本文使用于旭宅等[5]和陈崇成等[21]研究的验证方式，利用地理信息系统(Geographic information system，GIS)在高分辨率无人机正射影像上进行目视解译勾绘树冠轮廓作为验证。从图1可看出，实验区坡面西南侧树木冠幅相对较小，以2～4 m居多，东北侧树木冠幅相对较大，分布在4～6 m不等。以此为依据划分了2块验证区域：区域A约为66 m×50 m，包括目视解译的树木246棵；区域B约为54 m×38 m，包括目视解译的树木112棵。区域A可加强表现分割算法的总体适用性，区域B树木可分析局部效果。对GIS勾绘的冠幅多边形计算几何中心作为参考的单木位置，多边形面积作为参考冠幅面积。

3.1 单木检出精度

样地树顶检测结果如图5a所示，算法识别单木位置为图中黑点，综合了图像检出点和高程信息检出点。单木检测精度以检测率表示，检测结果有正检、漏检和过检3种情况。正检定义为在一个参考树冠上，有且仅有一个检测出的树顶点(图5b)；漏检定义为一个参考树冠上没有检测到树顶点(图5c)；过检定义为在一个参考树冠上检测到两个或多个树顶点或不是树冠的位置检测到树顶点(图5d)。图5b～5d中，黄点为目视解译得到树顶点，红圈为算法识别树顶点。正检率D是评价结果的主要指标，用户精度U突出过检(错分)，生产者精度P突出漏检(漏分)，精度定义为

(4)

(5)

(6)

式中Nd——正检数

Nc——树冠内检测出树顶的参考树冠数量

Nr——参考树冠数量

Nall——算法检测出的树冠数量

统计2个区域各项检测指标如表1所示。区域A检出率为91.46%，区域B检出率为92.86%，略高于区域A。可发现总体检出率均较高，对于2～6 m冠幅的树木而言，冠幅变化对检出率的影响不大。但由于算法主要依赖高程信息进行树顶检测，树冠越大，DSM冠型越趋于实际形态，因此检测率略有提高。

表1 实验区单木检测结果Tab.1 Individual tree detection results of experiment regions

3.2 树冠面积提取精度

检出率反映了算法对初始单木位置的识别率，但检出的树顶点在树冠分割过程中未必与参考冠幅有较好的匹配。为了进一步评价分割算法的总体性能，定义分割冠幅面积相对误差为

(7)

式中S——算法提取的树冠面积

Sr——对应位置参考树冠面积

基于相对误差，可使用位置和树冠面积综合判定的策略：在树木中心位置符合正检判定的前提下，冠幅面积相对误差低于判定阈值(δt)时视为正确提取，故正确率定义为

(8)

两个区域的单木树冠提取结果如图6所示，彩色多边形分块表示为本文算法分割树冠结果。可看出本文方法对树冠的提取无论在位置上还是形态上均基本符合实际树冠情况。少数树冠存在“直边”分水线，发现是由于三维重建时DSM模型本身对茂密冠层没有较深的透入，不能准确反映树冠边界轮廓所致，这一现象对面积提取影响不大。

图6 提取树冠结果Fig.6 Results of individual tree crown segmentation

图7 相对误差阈值对提取结果的影响Fig.7 Influence of relative error threshold on extraction results

对人工勾绘树冠面积和算法提取面积进行线性回归并计算决定系数R2，由于相对误差阈值δt对正确率存在影响，计算不同δt下正确提取树冠数量和面积R2的关系如图7所示。由图7可看出，R2随δt增加而减小，因相对误差阈值增加时，保留了更多面积误差较大的树木，从而降低了面积提取精度，同时树木检测数量(右侧纵轴)随之增加。

当δt取50%时，区域A正确率为76.02%，R2=0.834 9；区域B正确率为82.14%，R2=0.817 2。线性回归如图8所示。可以发现，主要为小树的区域A树冠面积精度略高于大树居多的区域B，区域B内部算法提取树冠面积总体偏大，但区域B正确率高于区域A。面积精度的差异主要来自DSM模型：树冠之间存在连通性。对于常规DEM支持的分割算法，一般剔除一定高度阈值以下的CHM像素，以避免低矮植被和地物的影响，从而很大程度消除了树冠之间的连通性。本研究由于不使用高度归一化，故有少数树冠之间的像素误分为树冠边缘(图4a～4c)。通过密集匹配点云提取单木，因难以得到冠层以下及地表观测，SfM三维重建时仅在最顶层以拉伸填充的方式表示冠层，在茂密林区不可忽视。由于区域B树冠较大，株数密度较低，导致树冠边缘像素引起的误差也更明显。对于正确率的差异，由于区域B大树的检出率较高，且面积分布更为均匀，相对误差阈值变化引起的正确率变化较小。

图8 树冠面积回归关系Fig.8 Regression relations of tree crown areas

林地树冠参数提取结果如表2所示。由表2可知，本文结果满足林业资源调查对树冠提取的需求。

表2 单木树冠提取结果与对比Tab.2 Results and comparisons of individual tree crown extraction

4 讨论

4.1 算法比较

利用多种方法进行了树冠分割提取的比较。实现了传统的应用于全图的分水岭分割算法[12](以下称为全局分水岭)和SILVA等[13]提出的针对LiDAR数据CHM的分割算法(以下称为SILVA2016)，正确率和面积R2比较如图9所示。可见，不同算法在区域A的正确率差异不大(76.00%～79.00%)，区域B本文算法正确率为82.14%，略高于全局分水岭(80.36%)，明显高于SILVA2016(75.00%)。因不同算法使用的DSM模型相同，且均主要基于高程信息检测初始单木，同时受到相对误差剔除的影响较小。对于面积精度，在区域A本文算法(R2为0.834 9)明显高于全局分水岭(R2=0.685 8)和SILVA2016(R2=0.509 4)，区域B内本文算法与全局分水岭接近(R2为0.817 2和0.799 0)，远高于SILVA2016(R2=0.485 2)。由比较可知，本文算法避免了在未高程归一化条件下的分割错误问题，在区域A体现最为明显。区域A地形高差变化较大、树冠为中等尺寸，全局分水岭在没有高程归一化的条件下很难准确分割。

图9 算法精度横向比较Fig.9 Comparison of algorithm accuracy

图10为全局分水岭和SILVA2016方法的树冠面积回归情况，可以看出全局分水岭同样存在对大树树冠面积的高估(图10c)。与图8对比可发现，本文算法对面积提取的准确度更高，集中度更好。

图10 算法的面积精度比较Fig.10 Area accuracy compare of algorithms

4.2 精度影响因素分析

本文单木树冠分割算法结合了三维和二维数据处理流程，因此各环节处理效果对最终树冠提取结果均有直接或间接影响。研究发现，SfM三维重建对林木参数提取影响最大，其精度直接决定树冠形态是否真实。SfM重建受到几个关键因素影响：①图像本身拍摄质量，反映为图像的清晰度、畸变和空间分辨率等特征。②图像重叠度。③GPS初始坐标精度。较差的影像质量对同名点搜索匹配质量造成不可逆的降低，很难在算法处理环节恢复。因此，无人机作业应尽量在光照稳定、无风条件下进行，避免光照阴影和过度曝光，防止枝叶的摇晃造成同名点配准误差。尽量提高影像重叠度以增加冗余观测，可有效消除影像质量不佳对SfM总体结果的影响，同时可在不同航高叠加拍摄避免地形高差变化带来航拍漏洞。在SfM处理中，利用区域网平差可有效调整曝光点坐标，因此无人机GPS精度对结果的影响相对可控。

同名点搜索与匹配直接影响点云坐标的精度，错误的冠层点云会使树冠形态发生扭曲和形变，降低冠层的天然连续性，最有可能形成过检(图5d)。由于生长季树叶茂密，无人机无法观测到林冠下层，表面模型在树冠边缘和间隙普遍存在一定的模糊性，无法像激光雷达点云一样准确地描述树梢细节，会导致分水岭分割算法出现“直边”现象，也使分水岭提取树冠面积略高于真实面积。为了解决该问题，需提高无人机成像的清晰度和空间分辨率，避免在树冠间隙进行错误的点云加密。本文算法选取kNN邻元参数为4，可有效进行分割，对更复杂林分可尝试增大k或进行二次分水岭对边缘进行优化。

复杂地形条件对无人机遥感单木识别具有较大的挑战性，本研究区坡面地形条件复杂，坡度陡峭，代表了典型山地环境下的林分地形条件。一方面，复杂地形增大了无人机航线规划难度，为了获取完整坡面数据，使用了多个航高叠加的成像方式。另一方面，复杂地形条件下，地形坡度大、高程变化剧烈，无人机航拍获取林下DEM存在很大的精度问题，采用DSM和DEM求差获取CHM的方法受到了制约，逐像元的差值运算甚至可能扭曲原始冠型。本研究针对性地采用了DSM分割方法解决了这一问题，避免了地形变化对冠层提取的精度造成影响，可达到与传统CHM分割方法相当的精度。

5 结束语

针对目前无人机技术在森林资源调查中的应用现状和实际林区的特点，提出了一种基于SfM三维重建的kNN自适应邻域分水岭分割算法，充分利用密集航片的数据冗余特点进行三维表面重建，利用高程和纹理信息结合的方式实现了较高的树木检出率。自适应邻域分水岭分割方法无需无人机拍摄林下地表，无需第三方DEM数据进行高度归一化，即可准确分割树冠，同时，避免了高分辨率图像的过分割问题，有效提升了无人机影像在森林资源调查中的实用性。