程晓菲，武 刚

(北京林业大学信息学院，北京 100089)

单木是组成陆地生态系统的最小乔木实体。随着新一代信息技术特别是遥感影像处理技术的发展与普及，高效采集海量单木信息建立单木数据库，是实现精准林业特别是城市林业精准集约化管理的基础，也是新一代智慧林业发展的目标[1-4]。事实上，国内一些个别的林业领域已实现单木集约管理，如古树名木、北京市杨柳雌株精准管理、果树管理等[5-6]，但这些主要通过传统的地面调查实现，效率低且工作量巨大。基于遥感影像的单木提取是实现单木数据库高效建设的关键技术，其中基于遥感影像的单木定位(探测)是基础。

近年来遥感平台和相应的处理技术发展很快[7]，虽然主动遥感技术(如激光雷达LiDAR等)已用于单木检测(包括单木位置、树高等)和反演模型建立等方面研究[8]，但其成本较高且不适用于大规模获取。而基于高分辨率光学遥感影像的单木提取方法，因为其具有技术基础深厚、影像尺度大且易得性高、便于推广应用等特点，依然是研究热点[9-10]。大量关于被动遥感的研究表明：遥感数据源[11]、季节及树木生长状况[12]、定位方法选择[13]等均为影响单木定位效果的重要因素。针对单木定位问题，学者们提出了以局部最大值法[9]为代表的探测单木树尖的直接定位方法，以及单木树冠探测与绘制(ITCD，Individual Tree Crown Detection and Delineation)等间接定位方法[14-15]，如各种改进的分水岭法[16-17]、谷底跟踪法[18]、区域生长法[19]和标记点过程[3]。近年来学者们也尝试各种定位新方法，如爬峰法[20]、Snake模型[21-22]、机器学习方法[23]等。但目前单木定位方法精度难以达到实用要求，还没有一种方法能够适应各种林分[24-25]。对于ITCD而言，其主要问题是单木排列紧密及冠层交错均会导致单个树冠内、重叠树冠间出现过分割或欠分割问题，导致单木定位产生误判和漏判问题[26-27]，因此需要探究更高效且稳定的单木定位新方法。

Chan等[28]提出了一种基于区域活动轮廓模型的图像分割方法，即Chan-Vese模型(简称CV模型)。CV模型利用像元灰度区域的统计信息，能够较好处理模糊边界或者离散状边界的目标[29]，已广泛应用于医学、农业等领域[30-33]。但是，CV模型应用时，初始轮廓的大小、位置及形状，均影响其分割的效率和结果[34]。实际应用中，通常采用人工介入的方法将包裹目标物体的圆或矩形作为初始轮廓，其效率低且在一定程度上忽略了图片的视觉特征[35]。

因此，本研究基于树木的形态学及其影像光谱特征，在自动提取初始轮廓的基础上，通过对多个实验区高分辨率卫星影像进行多方法比较实验分析，探讨基于CV模型的单木定位方法应用效果。

1 材料与方法

1.1 研究区域及其遥感影像

首先选择了北京林业大学及其附近的2个实验区，通过目视和实地测量，进行了初步试验，并与河北文安界围农场、深圳宝安荔枝园、山东栖霞市苹果园、浙江临安青山湖绿道的4个实验区实验影像进行比较，具体的实验影像数据及其参数如表1所示。试验区1为北京市海淀区的东升八家郊野公园，主要树种为针叶圆柏(Sabinachinensis)，该地区树木密度较高，树冠重叠现象较为严重。裁剪的实验影像1如图1A中蓝框区域所示，其中目视定位且实地验证共有397棵参考单木。试验区2为北京林业大学校园，裁剪的实验影像6(如图1F中蓝框区域)，该影像含有篮球场及其周边，目视定位且实地验证了24棵参考单木，其树冠茂盛，重叠严重，树种主要为阔叶的国槐(Sophorajaponica)和馒头柳(Salixmatsudanavar.matsudanaf.umbraculifera)。它们分别代表现实单木存在的林分和非林分(四旁树)两种形式。

后续为了对比考察CV模型单木定位法对不同类型影像源的适用性及其定位效果，选择了试验区2的”北京二号”影像7(图1G)，以及目视定位可靠且代表不同类型(针叶林、阔叶林、经济林)的4幅林分影像2～5(图1B—1E)[19,36-39]。

图1 研究区实验影像Fig.1 Experimental image of study area

1.2 图像预处理

一般单木定位(探测)的图像预处理包括校正、配准、增强、树木区域提取等环节，它们是准确提取单木位置信息的基础。本研究遥感影像已进行过校正和配准，实验只涉及图像增强和树木区域提取2个环节。

为提高图像的层次感，增强树冠和非树冠之间的差异，减少同一树冠内由于小树枝及阴影引起的噪声，采用高斯滤波方法增强图像[14]。

树木区域提取的目的是尽可能地剔除遥感影像中非树木区，如土壤、街道、房屋、操场、草地等，以保证定位处理方法仅作用于树木区域。通常在高分辨率遥感影像中，采用归一化植被指数(NDVI)用于掩膜运算。

采用绿度分割方法[24,40]：

dgreen=(R-G)/(R+G)。

(1)

其中：dgreen为绿度，R和G分别代表红色波段和绿色波段光谱反射值。基于绿度统计直方图，使用最大类间方差法(OTSU)确认阈值[41]，去除小于阈值的像元点，以此将图像中的道路、草地等无关信息剔除，获取树木区域图像I。

1.3 基于CV模型的定位处理

1.3.1 初始轮廓提取

CV模型常用于获取目标物体的精确轮廓，首先沿着物体边缘大致设置一条闭合曲线，使其在约束条件下，朝着内法线的方向移动，最终在目标边界停止。在传统方法中，通常采用基于人工交互的方式将其设置在目标边缘。单木定位中，因树冠数量诸多，手工设置初始轮廓工作繁重且不实用，宜采用自动化方法实现初始轮廓提取。根据树冠顶点受太阳光照强度较大，在影像中通常具有较高的光谱反射值，而树冠边缘反射值较低[42]，本研究通过计算影像中多个局部最大值区域并进行连通域标记，作为初始轮廓。其中，局部最大值区域是指连通且值相同的像元区，其周围像元值均小于区域的像元值。具体提取过程如下：

1)设树木区域图像I的像元(x,y)灰度值为G(x,y)，其8邻域内像元灰度值为Gneigh(x,y)，遍历每个像元，则由公式(2)可得到一个二值图A。

(2)

2)对二值图A中的局部最大值区域采用4邻域方法进行形态学连通域查找并标记[43]，获取初始轮廓标记图C=C1∪C2∪…∪CL。其中，Ci为第i个连通域，Ci(x,y)=i，若记Ci在树木区域图像I中对应的区域为Ωi,i=1,2,…,L，则Ω=Ω1∪Ω2∪…∪ΩL，即为CV模型演化的初始轮廓。

1.3.2 树冠轮廓分割

基于CV模型对初始轮廓进行演化迭代，主要是通过轮廓曲线的能量函数值的不断优化完成。树冠轮廓曲线由水平集函数φ(x,y,t)表达，曲线把图像划分为内外两个部分，c1和c2分别记为曲线内部和外部的图像像元平均值，将图像的能量函数定义为：

(3)

其中：μ≥0、v≥0、λ1≥0、λ2≥0是各个能量项的权重系数；x、y分别为图像横纵轴坐标；当曲线逼近树冠边缘时取值接近0，且c1=c2，此时曲线ci位于目标树冠边界。H(φ)和δ(φ)分别为Heaviside函数和Dirac函数的正则化表示，用于划分演化区域，限定轮廓在演变过程中在零水平集函数附近取值，其表达式为：

(4)

(5)

式中，ε为数字运算选取的一个小正数(阈值)。

依据变分原理和梯度下降原理，对水平集函数φ求式(3)能量泛函极小化，得到其偏微分方程：

(6)

式中：v为能量项的权重系数，φ0(x,y)为初始水平集函数。

式(6)可通过有限差分法求解，即设空间步长为h，时间步长为Δt，则得到迭代方程为[28]：

(7)

树冠轮廓具体迭代计算步骤为：

①n=0，根据待分割图Ω，建立初始水平集函数φ0(x,y)，并根据式(3)构造能量函数；

②根据当前树冠轮廓迭代结果φn，按照式(7)计算φn+1；

③若|φn+1-φn|≤ε(阈值)，则N=n，转步骤④，否则n=n+1,转向步骤②；

④输出最终轮廓(树冠曲线)φN(x,y)。

1.3.3 单木位置计算

对上述分割树冠轮廓二值图φN(x,y)进行连通域查找，基于图像几何矩原理，用公式(8)计算其重心作为单木位置坐标(x,y)[44]，连通域个数即为提取的单木数量。

(8)

其中：φN(x,y)为二值图的像元值；C10和C01为图像一阶矩，分别表示树冠轮廓x坐标值与像元灰度值乘积的累加和以及树冠轮廓y坐标值与灰度值乘积的累加和；C00为图像零阶矩，表示树冠轮廓像元灰度值的累加和。

1.4 定位结果评价指标

影像1的局部影像机器定位与目视定位的对比结果见图2，其中图2B为CV模型定位结果与目视定位的对比，图2C为标记分水岭定位结果与目视定位的对比。

图2 影像机器定位与目视定位对比Fig.2 Comparison between image machine positioning and visual positioning

2 结果与分析

以配有Intel Core(TM)i5、2.60GHz、12GB内存、64位Windows 10操作系统，利用PyCharm Community Edition 2018.3.5、Python编程语言及其OpenCV、MATLAB. Engine等组件开发的算法实验平台，对7幅影像以目视解译数据作为参考，采用CV模型和3种传统的主流单木定位方法(梯度分水岭法、标记分水岭法和局部最大值法)进行了实验比较分析。

2.1 CV模型的单木位置图像精度

影像2的各算法实验影像见图3。其中图3A—3D分别为图像原图、初始轮廓图、CV模型树冠轮廓图以及基于CV模型的单木位置图像。图3E—3G为局部最大值法、梯度分水岭法以及标记分水岭的单木定位图像。可以看出，基于树冠的光谱反射原理，初始轮廓大概标记了单木位置，并在此基础上，进行水平集函数的构建及迭代，使树冠轮廓逐步收敛到目标边缘，能够较好地识别单木。相较于局部最大值的漏判，梯度分水岭多棵单木误识别为1棵，CV模型能够较好地进行单木定位。7组实验的单木位置精度评价见表2，CV模型均取得了最高匹配率。

在影像1排列紧密、大小均匀且排列较为规则的针叶林中，因其树冠辐射值较高且分布较均匀，局部最大值在3种传统方法中匹配率(M)最高，梯度分水岭及标记分水岭等传统方法将多个小树冠误合并识别为1个树冠，因此漏判率(Rom)较高；影像2～3中，树冠大小不一且排列凌乱，局部最大值方法无法准确获取单木位置，导致了一定程度的漏判和误判，树冠间无明显区分，梯度分水岭及标记分水岭法易将树冠误分割为多部分， CV模型的准确率(Rmat)较高，并且误判率(Rcom)较低，能够较好地识别较小单木，处理漏判现象，具有一定的泛化效果。在影像4～7中，阔叶树树冠较大且树冠重叠严重，传统方法将树冠误分割为多个树冠，出现了一定的过分割现象，也因此导致了较高的误判率，CV模型较好地缓解了这一状况，但存在一部分相邻树冠难以分割的情况，而总体匹配率效果较好。

图3 影像2各算法实验影像图Fig.3 Image 2 experimental results of each algorithms

表2 单木位置精度评价

2.2 CV模型影像定位效果

所有实验影像的整体定位效果见表3，局部最大值方法因其受移动窗口因素影响较大，同一树冠(特别是阔叶树)中易产生多个局部最大值点，因此导致较高的误判率(Rcom)，平均为0.434。梯度分水岭和标记分水岭通过在原始图像上分别计算图像梯度函数及标记树冠大概位置，在其基础上进行分水岭分割，容易产生过分割现象，误判率较高。相较于其他3种定位方法，CV模型在7组影像实验中均取得了最高的匹配率(M)，平均高出其他方法约23%，综合定位效果最好。在同一地区不同影像中(影像6～7)，不同类型树木中均具有较好的定位效果。但CV模型在树木排列紧密区域，目视不清时，也容易存在一定的欠分割问题，导致漏判率较高(影像3～7)。影像6的局部区域实验图像见图4。该区域树冠紧密，CV模型将排列紧密树冠误认为同一个(图4B)，仅识别出1棵单木(图4C)。

表3 整体单木定位精度

图4 CV模型局部实验图像Fig.4 Local experimental image of CV model

3 结 论

本研究提出的CV模型单木定位法，首先利用树木的形态学特征及其在影像中的光谱特征，自动设置初始轮廓；其次，进行基于水平集函数构建及迭代，对树冠轮廓进行演变；最终获取单木位置。该方法基于图像全局信息，利用曲线内外的灰度均值而不是梯度信息进行分割，能够在边界模糊或梯度无意义的图像中取得较好的分割效果，能够快速准确地收敛到目标位置。实验表明：该方法比其他3种方法(梯度分水岭法、标记分水岭法、局部最大值法)具有更高的匹配率，平均高出23%，且与影像源无关，显示了良好的泛化应用潜力。

更全面地评价单木信息还需要进行更多方法(如区域增长、模板匹配、机器学习等)和更多影像的对比实验。同时，实验也表明，虽然CV模型的单木定位法一定程度上解决了单木树冠识别的过分割及欠分割问题，综合定位效果有所提高，但在郁闭度较高、目视都难以识别单木的情况下，仍存在一定的欠分割问题；后续可针对郁闭度差异进行单木定位实验研究，探讨其适用性。随着机器学习方法在目标检测方面的突破，新的单木定位方法值得期待。另外，结合多时像(春、夏、秋、冬)高分辨率影像的融合应用技术的研究，或许对解决郁闭度高、单木难定位识别问题有所帮助。