改进分水岭算法在无人机遥感影像树冠分割中的应用

于旭宅, 王瑞瑞, 陈伟杰(1.北京林业大学林学院，北京 10008；2.北京林业大学精准林业北京市重点实验室,北京 10008；.天津市静海公路工程有限公司，天津 01600)

随着电力产业的飞速发展，输电线路在荒野分布日渐广泛.每年树木生长茂盛期，树冠都会对输电线路的安全运行造成很大威胁[1].利用无人机遥感技术可以快速获取输电线路下行树冠的高分辨率遥感影像，通过提取树冠并结合树冠高度模型及输电线高程数据，计算树冠与上行输电线的垂直距离，这对于开展输电线路安全巡检具有重要意义[2].目前无人机遥感技术在输电线路巡检等领域得到广泛应用[2-3]，如何有效、准确地从遥感影像上提取树冠信息已成为国内外众多学者研究的热点[4-7].在遥感影像分割领域，分水岭算法可应用于多光谱遥感影像、高光谱遥感影像和高分辨率遥感影像的分割[8-9]，逐渐成为一种普遍采用的分割技术.张博等[10]提出了一种改进分水岭变换的高分辨率遥感影像多尺度分割方法，在抑制分水岭过分割现象的同时，还能实现对遥感影像的多尺度分割.Sharif et al[11]提出了一种基于梯度修改的改进方法，利用形态学算子对梯度图像进行修正，减少了过分割现象.

由于无人机多光谱遥感影像中输电线路下行林区存在同一冠层内部纹理细节显著且亮度变化不均匀、不同树冠之间重叠遮挡等现象，导致现有方法在提取树冠时存在同一树冠内过分割、不同冠层间难以分割等问题.为此，本文提出了一种基于NDVI植被指数计算的改进分水岭分割方法，通过归一化植被指数模型(normalized difference vegetation index, NDVI)计算树冠显著性区域图像，利用彩色向量空间算法计算其梯度，并采用基于标记控制的分水岭算法，实现了树冠的准确分割.

1 研究区概况

研究区位于安徽省池州市青阳县东北部的锦苏输电线路下行通道，地理位置为东经117°48′—117°54′，北纬30°40′—30°41′，地势南高北低，海拔1.8～112.2 m.该区属于亚热带季风湿润气候区，夏热冬寒，春秋温和，雨量充沛，日照时间长，无霜期短，梅雨期40 d左右.研究区内架设有输电线路，树种主要有泡柳、松树、水杉、毛竹和杂木等，树木生长茂盛.

数据源为无人机遥感影像，影像获取时间是2016年6月16日，空间分辨率是0.05 m.从该无人机遥感影像可以获取输电线路通道下林木6个谱段的光谱反射率信息，第1、2、3波段分别是红、绿、蓝波段，位于可见光区域；第4、5、6波段的波长分别为600、701、750 nm，带宽均为10 nm.

2 研究方法

2.1 显著性区域提取

利用植被的光谱特性和形态特征，对原始无人机多光谱遥感影像进行预处理，得到树冠的显著性区域图像，主要包括树冠信息的去相关性分析、树冠区域提取、形态学滤波等预处理过程.

2.1.1 多光谱影像的相关性分析 研究区的无人机遥感影像有6个波段，若对每一个波段都进行计算，会造成数据量急剧增加、耗时、冗余增大，因此首先对无人机遥感影像进行波段信息量、波段间相关性分析、比较，选择信息量最大且波段间离散度最大、相关性最小的波段组合作为分割图像.利用ENVI软件中的Statistic工具统计出各波段相关系数，结果如表1所示.从表1可以看出无人机数据的第1、6波段相关性最小，第2、6波段相关性较小，但考虑到后续需要提取植被指数NDVI，故选择第1、2、6波段；该3个波段按照RGB的形式合成(标准假彩色合成)，此时通过目视解译提取树冠，表达效果最好.

表1 无人机影像各波段相关系数矩阵Table 1 The correlation coefficient matrix of each band of UAV images

2.1.2 树冠区域提取 植被的光谱特征能够将其在遥感影像上有效地与其它地物相区别.考虑到无人机遥感影像中输电线路等背景区域的影响，选取第1、6波段，利用NDVI将树冠区域与输电线路等背景区域区分开来.通过比值运算进一步扩大树冠区域与输电线路等背景区域的差距，从而增强树冠区域的亮度，抑制输电线路等背景信息.针对比值运算的结果进行直方图分析，通过设定阈值范围对NDVI计算结果做分类处理，阈值范围外用0表示，即输电线路等背景区域；阈值范围内用1表示，即树冠等显著性区域.并生成掩膜图像，对原始影像做掩膜处理，从而分离出输电线路等背景区域和树冠显著性区域.NDVI比值计算如下：

(1)

式中，NIR为近红外波段反射率，R为红光波段反射率.

2.1.3 形态学滤波 无人机遥感影像分辨率特别高，影像上树冠内部的纹理细节特别清晰.为了缩小纹理细节对后续图像的过分割影响, 采用形态学开闭重建运算对显著性区域影像进行进一步预处理.形态学开闭重建运算是建立在腐蚀和膨胀的基础上.对于原始图像f(x,y)和参考图像r(x,y)，其形态学膨胀运算[12]如式(2)所示.

(2)

式中：B表示圆盘状结构元素.圆盘形结构元素满足旋转不变性,因此不会造成图像特征值的畸变.

原始图像f(x,y)的形态学腐蚀运算[12]如式(3)所示.

(3)

OB(f,r)=DB[(fB),r]

(4)

CB(f,r)=EB[(f·B),r]

(5)

形态学滤波计算中结构元素B的大小和形状对滤波效果有较大影响.如果结构元素过小，一些树枝会被当作树冠；如果结构元素过大，一些较小的树冠会被移除.一般来说，结构元素的形状应与所需提取的地物形状近似，尺寸应小于所有待分割地物尺寸的最小值.因此，结构元素采用圆盘形状，尺寸小于最小树冠半径.通过多次反复比较，结果表明结构元素尺寸大小为1时效果最好，树冠内部的纹理信息可被有效抑制.

2.2 影像梯度计算

一般来说，利用梯度图像作为处理对象比直接在原始图像上分割得到目标边界更为准确.由于分水岭算法与影像的梯度有关，而与影像本身的相关性不大，因此梯度图像的构建对后期分水岭的分割效果有着至关重要的影响.考虑到显著性区域图像含有3个波段信息，且在彩色向量空间直接计算梯度比以单独的分量图像为基础计算梯度具有更高的准确度，因此本文采用彩色向量空间梯度算法[13]，直接在RGB向量空间计算梯度.设r、g、b是RGB彩色空间沿R、G、B轴的单位向量，像素点沿水平方向和垂直方向的彩色梯度可用向量u、v来表述：

(6)

(7)

gxx、guu、gxu分别定义为这些向量的点乘，表述如下：

(8)

(9)

(10)

由此可得到单个像素点(x,y)在最大变化率方向θ的彩色梯度公式为：

(11)

其中，最大变化率方向θ为：

(12)

2.3 树冠标记提取

由于显著性区域图像存在较多局部最大值，直接采用现有的分水岭算法对梯度图像进行分割，会导致图像严重过分割，而有效地选取标记是标记分水岭的重要步骤.采用扩展最小变换(E-minima)[14]提取标记，其本质是一个形态学阈值算子，将大多数无关小区域记为0.梯度图像G的高度阈值(h)扩展最小变换运算表示如下：

E=ME(G,h)

(13)

式中，E为输出的二值图像.具体算法过程为：局部极小值与h进行逐一比较，将盆深大于h的局部极小值标记为1，其他标记为0，最终形成二值标记图像.

由于树冠标记中某些标记连接在一起，使得分割结果出现多个树冠而被误分为一个树冠.为使分割效果更理想，结合显著性区域树冠图像的物理特征(树冠中心附近区域辐射最大)，采用形态学开重建和闭重建技术来修正二值标记图像.首先利用形态学腐蚀运算将连接在一起的标记分开;然后利用形态学距离变换确定每个标记的范围，确保每个标记最终分割得到单一树冠;最后运用强制运算将标记加到梯度图像上，得到修正后的梯度图像.

2.4 冠层分割

通过扩展最小变换预先提取树冠顶点作为标记，将标记强制应用为局部最低值，再进行分水岭分割.标记控制分水岭算法从梯度图像的全局极小值点开始，假定当前值为h+1，每个极小值小于或等于h的集水盆会被分配一个唯一标记.对当前值为h+1的像元，若其邻域像元已经标记过，则给它分配相同标记；若邻域像元都没有被标记过，则被认为是一个新的集水盆，并给它分配一个新标记.如此反复标记，直到图像中每个像元都归属到某个集水盆(即属于某个对象区域)为止，分水岭分割则完成.

3 结果与分析

3.1 显著性区域提取

为了缩小输电线路等背景区域的影响，对研究区的无人机多光谱遥感影像进行了波段甄选、NDVI植被指数计算、冠层掩膜等处理，结果如图1所示.根据NDVI计算结果进行分类处理，建立冠层掩膜.从图1D可以看出，输电线路、树冠间空隙与树冠都被有效分离出来.利用冠层掩膜得到去除输电线路等背景后的冠层区域，结果如图1E所示.从图1E可以看出，输电线路等背景被有效去除，为后续分水岭分割冠层分割提供了可靠的影像数据.

3.2 分割结果与精度

为验证本算法的准确性与有效性，采用基于标记控制的分水岭算法，对显著性区域图像进行树冠提取，结果如图2A所示.由于缺少外业调查数据，利用ArcGIS软件通过人工目视解译方式获得了参照图像(图2B).通过改进算法得到的图像与参照图像的差异度来衡量分割结果质量，由于位于边缘的树冠形状不完整，无法应用分水岭算法，因此评价分析不包括边界上的树冠.

为验证算法的准确性，将论文改进算法与目视解译做对比分析.采用改进算法提取142个树冠，其中输电线路下行的树冠26个，如图2A所示；根据树冠顶端高亮度像元与树冠形态特征进行目视解译，结果如图2B所示，提取127个树冠，其中输电线路下行的树冠有23个.从表2可知，从树冠提取总个数来看，通过改进算法提取的个数多于目视解译.这是因为在目视解译过程中忽略了较多树冠特征不明显的区域，但这些区域又存在高亮度像元，所以改进算法将其当作树冠进行提取.从树冠重叠个数(表2)来看，改进算法提取的树冠总个数和目视解译匹配的个数有117个,改进算法提取树冠总个数匹配精度达到92%;改进算法和目视解译在输电线路下行的树冠匹配个数有18个,改进算法提取的输电线路下行树冠匹配率为78.2%.这是因为在去除输电线路背景、提取显著性树冠时，输电线路将部分树冠分为2个树冠，或者将同一树冠分裂成多个树冠，导致改进算法提取个数多于目视解译个数.

A.无人机原始数据；B.无人机数据增强显示；C.NDVI植被指数计算；D.冠层掩膜；E.去除背景后的冠层区域.图1 显著性区域提取结果Fig.1 The extraction results of significant region

图2 算法提取结果(A)与目视解译结果(B)Fig.2 The extraction results of algorithm (A) and visual interpretation (B)

表2 提取树冠个数对比Table 2 The number contrast of extraction canopy

为验证算法的精度，进一步从空间上对改进算法进行单木树冠提取精度分析.随机且均匀地选取20株单木树冠，通过目视解译判断标记(图3中标记的数字)的树冠是否被分割出来，并从单个树冠所占像元个数进行误差分析(表3).在进行精度分析时，所用相对误差公式为：

(14)

式中T计算表示改进算法提取的单个树冠像元总数，T目视表示目视解译提取的单个树冠像元总数.

样本精度计算公式为：

(15)

式中,n为样本数，δi为第i个样本的相对误差.

表3 单个树冠像元总数计算结果对比Table 3 The result contrast of pixel numbers of a single canopy

由表3可知，改进算法提取的单个树冠像元总个数与目视解译相差不大.改进算法提取的单个树冠的平均像元个数为2 437，而目视解译的单个树冠的平均像元个数为2 457，就单个树冠的平均像元个数而言，二者相差甚小.利用样本精度计算可得，改进算法的平均相对误差为11.7%，形状精度为88.3%.总体而言，改进算法能够有效进行输电线路下行的单木树冠分割，不仅简化了数据分析的计算量，也可以达到精度要求，能够满足电力巡检、林业调查等领域的实际测量需求.

4 小结

本研究结果表明，改进算法能够有效去除输电线路等背景区域的影响，与目视解译提取的树冠结果相比较，获得的准确率得到了明显的提高.主要体现在以下3个方面：(1)对于整个研究区，目视解译提取树冠127个，改进算法提取树冠142个，树冠匹配个数为117个，在整体区域树冠个数匹配精度达到92%；(2)对于输电线路下行区域，目视解译提取树冠23个，改进算法提取树冠26个，树冠匹配个数为18个，在输电线路下行区域树冠个数匹配精度为78%；(3)对于单个树冠所占像元个数，改进算法提取的单个树冠的平均像元个数为2 437，目视解译提取的单个树冠的平均像元个数为2 457，改进算法的形状精度达到88.3%.