基于分水岭算法的阔叶林郁闭度提取方法的研究

严 羽，王永众，杨来邦，楼雄伟

基于分水岭算法的阔叶林郁闭度提取方法的研究

严 羽1，王永众2，杨来邦1，楼雄伟1

（1. 浙江农林大学 信息工程学院，浙江省林业智能监测与信息技术研究重点实验室，浙江 杭州 311300；2. 杭州感知科技有限公司，浙江 杭州 311300）

针对现有的人为调查林区郁闭度耗时，费力且林区条件恶劣等问题，本研究实现了一种使用无人机来测定森林郁闭度的方法。通过无人机在面积为50 m×70 m以内的5个阔叶林为主的地块中，在不同间距和高度拍摄4组图像并制作成正射影像图，通过对图像的灰度化以及滤波等预处理，使用改进的标记控制分水岭的分割算法来提取树冠，并与人工提取做比对，经过实际试验该算法有较高的准确度，所得误差在5%左右，并对提取的树冠使用样线法计算郁闭度，结果表明对0.5 ~ 0.9之间的郁闭度有较高的精度。

无人机；分水岭分割算法；郁闭度；阔叶林

郁闭度指森林中树冠在阳光直射下在地面的总投影面积（冠幅）与此林地（林分）总面积的比，它反映林分的密度。郁闭度是反映森林结构和森林环境的一个重要因子，通过郁闭度能估算树木利用生活空间的程度，林分光能利用程度以及森林的生长状况。同时郁闭度大小也会影响着当地的土壤状况和浅层动物、微生物、植被的生长状况。因此合理、准确的测定郁闭度对森林资源的经营有重要意义。郁闭度的测定方式多种多样，有目测法、树冠投影法、样线法、样点法、遥感图像判读法等[1]。从调查的方式来说，为了提高精确度以及减少人力成本，国内采用各种各样的仪器或者软件来测定。刘菊[2]利用GIS来测定郁闭度，先将林木定位，用皮尺和罗盘仪测定东西方向、南北方向的冠幅值，并记录数据；将数据导入GIS软件中，利用其空间分析功能，将每木位置进行定位，并计算每木投影面积，使传统的树冠投影法得到改进，进而实现了森林郁闭度的简洁快速测算。刘芳等[3]基于鱼眼照片对郁闭度的测定，采用Photoshop、ArcGIS等软件对利用鱼眼照片提取森林郁闭度的方法进行研究，通过照片裁切、灰度处理、影像重分类等处理，求算出森林郁闭度。同时与野外观测值进行比较，两两之间差值均在允许误差范围之内，证明该方法有一定的效果。杜晓明等[4]以遥感数据与森林资源一类清查数据为基础，探讨了用Bootstrap方法筛选最优郁闭度估测变量，用偏最小二乘回归方法建立模型估测森林郁闭度的可行性，结果表明无论是用所有变量构造的模型还是用所选最优变量构造的模型，郁闭度估测的相对偏差都在5%左右。

本文利用无人机的高空拍照优势，拍摄多组高分辨率照片，通过正射影像图处理，采用标记控制分水岭算法[5]同时结合样线法来计算郁闭度，并与实际测试值进行比较来验证该方法的可行性和准确度。

1 数据收集及处理方法

1.1 数据来源

试验数据选用位于浙江省杭州市临安区浙江农林大学校园周边以阔叶林为主的天然林样地，该林内地面植被较少且树木分布较均匀，结构简单。选用不同密度、不同树高的5块林地，选择50 m×70 m大小以内的地块作为数据来源。使用大疆悟Inspire 1无人机从上空拍摄，拍摄前使用第三方app Altizure，该软件能够设置好飞行路径、高度、速度、拍摄间隔等参数。为考察不同高度，拍摄距离对制作正射影像图的影响，拍摄时采用了30 m和40 m高度，水平拍摄间隔设定为2 m和5 m。无人机拍摄时相机镜头朝下，其中镜头的有效像素为2 000万，每块样地无人机共飞行4次，即不同高度，水平距离拍摄4组图像。每次拍完的一组照片通过pix4d软件制作成正射影像图，因航向重叠的关系制作出的图像远超过边界，后期再裁剪成原先规定大小的图像。每一块样地4组样张，供后续处理。

1.2 树冠分割方法设计

制作好的正射影像图，如图1 a所示，该图根据软件Altizure显示实际面积大小为50 m×70 m。树冠分割方法流程图如图2，分割步骤如下：

图1 正射影像图

Figure 1 Digital orthophoto map

图2 树冠分割方法流程图

Figure 2 Image segmentation of crown

（1）源图像预处理即均值滤波去掉噪点减少干扰，得到图1 a，将图1a灰度化后得到图1 b；

（2）判断图1 b灰度直方图是否有双峰并计算梯度图像，如图3 b；

（3）若有双峰则平滑灰度图，若没有则迭代法寻找阈值，二值化后分割为前景e和背景d。如图3 c所示，其中前景e范围为灰白部分，背景d范围为黑色部分（为区分d和e，将背景灰度像素直接用黑色代替了）；

（4）将e计算区域极大值作为内部标记，对d分水岭变换取出相邻区域作为外部标记；

（5）修正梯度图像c，使图像c在内部标记和外部标记处有极小值，并使图像平整；

（6）对修正后的梯度图像c进行分水岭变换；

（7）对分割后的图像在连通区域的进行合并，最终获得分割后的前景树冠图像。

图3 内外标记提取过程

Figure 3 Extraction of internal and external tag

1.3 郁闭度测量方法

采用类比样线法来计算林分郁闭度，按下述公式计算：

郁闭度 = 标准地上的两对角线树冠覆盖的总长度与两对角线的总长度之比。

取上述标记控制分水岭分割出来的前景图像计算对角线长度，如图1b所示。统计该图（分辨率为1 263×663）落在树冠上的像素数量为283对角线长2 852。则上述图郁闭度：283/2 852≈0.1。

2 树冠分割具体实现

2.1 分水岭算法

分水岭算法[6]基本思想是把图像灰度化后，每一点像素的灰度值大小表示该点的海拔高度，每一个局部极小值及其影响区域称为集水盆，算法的过程是模拟浸水现象，每一个极小值及影响域向外扩展，在两个集水盆汇合处构筑大坝，形成分水岭，即将图像的低海拔的汇水盆地分离出来。示意图如图4。由图4可知，图4a是灰度图，图4b是图4a的地形图，图4c和图4d表示水浸没的过程，最终的分割效果如图4e所示。

图4 分水岭算法过程

Figure 4 Watershed algorithm

2.2 标记控制分水岭算法

将三维地形图变为二维地形示意图，如图5。图5中a，b，c，d为局部极大值对应山顶；e，f，g，h为局部极小值对应山底。经过分水岭算法后由图4可知会在a，b，c，d边缘处形成“脊线”即分割线，分割区域为A，B，C，D。而在实际图像中，由于噪点或其他干扰因素，以及很多很小的局部极值点存在，使用该算法存在过度分割现象如图3a所示，反映在图5中即g，h点为非必要极小值点，会形成没必要的d边缘分割线，因此要通过合适的方法提取出与分割主体相关的标记，来引导算法使之避开噪声极值区域，从而避免过分割提高分割算法的准确性。

图5 二维地形图

Figure 5 Two dimensional topographic map

2.3 内外标记提取过程

对源图像开始操作前，要先进行降噪处理，处理方式为K近邻均值滤波[7]，在一个与待处理像素邻近的范围内，寻找出其中像素值与之最接近的K个邻点，将该K个邻点的均值替换原像素值，方法如下：

（1）以待处理像素为中心，作一个×的作用模板；

（2）在模板中，选择个与待处理像素的灰度差为最小的像素；

（3）将这个像素的灰度均值（中值）替换掉原来的像素值。

该处理方式的优点是能够抑制噪声的同时，避免边界的平滑处理，可以大大降低对图像边缘的模糊。

为进行后续的分水岭变换，需要对预处理后的图像计算梯度幅值。因为梯度幅度图像沿着物体边缘有较高的像素值，而其他地方则有较低的像素值。处理方法为使用水平Sobel算子[8]，Sobel= { {－1, 0, 1 }, {－2, 0, 2 }, {－1, 0, 1 } };垂直Sobel算子Sobel= { {－1, －2, －1 }, { 0, 0, 0 }, { 1, 2, 1 } }；该算子包括3×3的矩阵，分别为横向和纵向，将与图像作平面卷积分别得出横向x及纵向y的亮度差分近似值，然后开平方根求得灰度值，如图5b。

为选出前景和后景标记，采用阈值法[9]进行初步分割。阈值的选取采用两种方式，若灰度直方图有双峰则平滑双峰，两个峰相对应的就是前景跟背景像素的极大值，由于灰度直方图中的随机波动，两个峰尖的最大值和他们之间谷底的最小值都不能被很好的确定所以要先平滑处理。否则采用迭代法计算阈值，过程如下：

（1）初始阈值，可以自己设置或者随机生成；

（2）根据阈值对每个灰度数据（，）分为前景灰度数据1与背景灰度数据2（为行，为列）；

（3）1的平均值是1，2的平均值是2；

（4）则新的阈值′=（1+2）/2；

（5）回到第二步，用新的阈值继续分前景像素与背景像素数据，继续2 ~ 4步，直到计算出来的新阈值等于上一次阈值。

阈值确定后将图像分为前景跟背景，其中背景像素用黑色（rgb=0）代替而前景树冠像素用原始灰度值显示如图5c所示。

灰度值越大，则物体颜色越明亮，所以树冠灰度值远远比地面等其他物体来得大所以可以在阈值法分割后得到的前景灰度像素上通过提取树冠的局部最大值来作为标记符。

在提取树冠局部最大值时，因为树冠内部有空隙，阴影，暗斑等因素干扰，所以直接提取的话会产生严重的误差，即可能会提取到伪树冠、枝干、噪声等部分所以提取前应该对图像进行平滑，去除暗色间隙以及非树冠亮细节等处理。

上述公式计算出来的值可能会大于255，则可以将这些灰度直接定义为最大值。如果模板的矩阵为一个常数d，则运算公式简化为

如图6所示运算过程如下：A的矩阵为常量灰度，中间有2×2黑色区域，模板为3×3的常值d模板，因为模板尺寸大于黑色区域，膨胀运算消除黑色区域，而且提升了A矩阵范围的图像灰度使得图像变明亮。

图6 膨胀运算示意图

Figure 6 Inflation method

若模板灰度为常量d，则

腐蚀操作示意图如下：

图7 腐蚀运算示意图

Figure 7 Corrosion operation

则图像形态学开闭操作公式为：

通过开操作即先腐蚀操作时可以去除小的亮细节，但会使得图像整体亮度变暗，然后再进行的膨胀操作会提升图像亮度但不会把之前腐蚀操作去除的细节引入，得到结果是去除微小亮细节区域且保持灰度级不变，如图8a，b所示，其中灰度模板为3×3矩阵[（85，255，85）,（255，255，255）,（85，255，85）]。同理闭操作（先膨胀后腐蚀）可以去除暗细节，比如本文图像树冠中的枝干，阴影等部位。通过上述操作能够维持总体灰度级不变和保持较大的明亮区域，且不会对边界产生偏移。

形态学开闭过程能够修正局部区域的极大值和极小值，减少噪声干扰。开操作过程能在提取局部区域极大值时避免提取到伪树冠部分（噪声极大值），剔除上述干扰后使得提取局部极大值时（树冠内部）定位更准确。

树冠灰度局部极大值是指连通区域且满足该区域领域周围灰度值均小于该连通区域的灰度值。上述操作后将会在每个树冠灰度值内部创建多个区域极大值，可以使用MATLAB自带的imregionalmax函数来定位。该函数的运算过程为：对于输入的×行矩阵，构建一个3×3的滑动窗口矩阵，从输入矩阵的头部开始滑动，每滑动一格比较中心元素值与4领域像素灰度值大小，若滑动到边缘位置则添加0值元素，输出一个×行的矩阵，其中局部极大值位置用1表示，其他位置用0表示。示意结果如下所示。

计算出的局部区域极大值作为内部标记，对图5c形态学重建后提取到的极大值叠加到原图上结果如图8c，其中白色像素为提取到的树冠内部局部极大值范围。外部标记不必要求太接近于要分割对象的边缘位置，所以可以对阈值分割后的背景图像结合分水岭算法取出相邻区域间的分界线构成外部标记符。

分水岭算法针对梯度图像进行，由前面可知梯度图像未改进时会产生过分割现象。所以将得到的内部跟外部标记符通过最小覆盖过程来改进梯度图像，该过程使得梯度图像仅在上述标记过的位置有极小值。对照图3b即标记处（目标物体）落在“盆地”处（极小值），而且剔除掉了不该出现的极小值点，从而能得到理想的分割结果。所以标记控制分水岭算法分割时最大限度的避开了原先的噪声极小值区域，这样避免了噪声极值区域而产生的过分割。标记后分割的结果如图8d，图8e显示的是错误标记的结果可以看出树冠较好的分割出来了但不足的是连带了地面部分原因是内部标记将该部分也标记在内，结果出现偏差。

图8 标记控制分水岭

Figure 8 Tag-control watershed

3 实验验证

3.1 树冠分割精度分析

图像分割的主流算法包括阈值法，边缘法，区域法，聚类分析法等，本文选取较为成熟的图像分割算法作为对比，分别是Otsu算法，基于Kmeans的聚类分割算法，以及本文的方法进行处理，得到的效果图进行误差对比。

图9 不同算法不同地块分割效果图

Figure 9 Image segmentation by different algorithms

图9中d，e副图分别是手动分割的跟本方法分割后为方便对比做的二值化图。从分割后的二值化图可以看出手动分割的与本方法分割的较为相似都是基于区域分割，而Kmeans跟Otsu算法只是比较了各个像素点的大小关系，而没有考虑到空间关系当图像的灰度差异不明显即图像比较复杂时分割效果比较差。本文提出的方法类似区域生长的分割方法，按区域属性特征一致的准则决定每个像素的区域归属，将属性接近的连通像素聚集为同一个区域。但对于复杂图像的处理使用区域生长方法会产生过分割现象，因为分裂的每个小区域需要一定的规则将小区域合并成大区域达到分割图像的目的而各个小区域过于分散会造成图像过度分割，还需要其他方法来将分裂的小区域整理成更好的的区域结构来避免过分割。

表1 不同分割方法平均分割精度对比

注：*标记控制分水岭算法。

3.2 算法性能评价指标

为定量分析该方法的优劣，参考文献[11]的评价体系来对算法图像分割与手动图像分割面积进行对比如表1所示。评测指标公式如下：

从数据可以看出Otsu跟Kmeans方法检准率较低而检全率高，侧面反映出手动分割和本方法分割的面积较上诉算法分割的大且本方法分割的面积较手动分割的小，而从Otsu跟Kmeans分割后的二值化图中也能直观的看出来白色像素部分比手动分割的少。

3.3 郁闭度测量结果验证

图10为不同地块的原图和控制分水岭算法分割后的效果图（顺序相互对应）。

图10 不同地块原图与分割后效果对比

Figure 10 Comparison on result of original photo and that after image segmentation

对5块样地中的每组4张正射影像图按照上诉方法流程计算郁闭度以及实际测量的数据记录如表2。

表2 通过算法计算的郁闭度与实际测定对比

注：误差=（手动分割-×控制分水岭算法分割）/手动分割。

统计结果如表3和图11，其中图11a横坐标为实测郁闭度，纵坐标为树冠分割相对误差均值，图11b横坐标为实测郁闭度，纵坐标为标准差，图11c为计算郁闭度与实际郁闭度对比图，图11（d，e，f）以测量的郁闭度均值为轴，实际郁闭度为轴，分别用线性关系、二次多项式、三次多项式对它们之间的关系进行拟合。

由图11a可以看出，当郁闭度比较大时树冠分割的误差较郁闭度较小时分割的误差小，可以看出样地树木密集程度与分割误差成反比，即越密集算法越容易区分前后背景。

由表3得4号样地计算郁闭度均值较实测郁闭度小且相差较大，其主要原因可以从图10（4）发现，即树木分布不均匀对角线上分布较少导致计算的郁闭度偏小。

表3 不同地块郁闭度均值方差与实测值大小比较

由图11b可知，实际郁闭度越大，计算的郁闭度标准差越小，精度更高。图11（c）曲线中除去4号样地的异常数据（蓝线在红线以下部分）可知计算值与实际值接近或偏大，当林分郁闭度较小时即0.5以下时（样地3，5）标准差较大，计算值比实测值偏大，当郁闭度较大时（样地1，2），标准差较小，计算值比实测值稍大，较为接近。

图11 统计图

Figure 11 Statistical chart

从拟合结果图11d，e，f来看3种不同的拟合方式的相关系数R2都在0.97以上，且三种拟合相关系数都较为接近，可以看出本方法测量的郁闭度与实际郁闭度有显著的相关性，随着拟合方程次数的提高相关系数变大，表明两者存在正相关性。

总体上来看，郁闭度较小时即0.5以下时，控制分水岭算法计算的郁闭度比抬头望法计算的郁闭度大。当郁闭度在0.5 ~ 0.9时，计算的郁闭度比抬头望法偏大或接近。其原因如下：控制分水岭算法是计算对角线上树冠的像素比上总的长度像素，但是实际计算郁闭度时不能忽略树叶之间的空隙即阳光正照射下来的亮斑，按理论是本文计算所得要比实际的偏大。当样地比较稀疏时，使用抬头望法根据目视判断样点是否被树冠挡住来测郁闭度，稀疏的林地更会让测量人员判断被挡住的样点偏少使得抬头望法测得的郁闭度比实际的郁闭度偏小，相反林地稀疏使得树冠与树冠之间的距离远大于树叶空隙，减小了空隙带来的误差，所以郁闭度较小时的结果为：本研究样线法≥实地样线法>抬头望法。当郁闭度较大时，因为树木密集容易判断，抬头望法计算的遮盖样点更准确，同时树冠与树冠之间基本上是挨着的，这样使得空隙的误差不能忽略，但是因为比较密集相对于稀疏的树叶来说空隙较少，所以得出的结论为郁闭度较大时本文样线法≥实地样线法≈抬头望法。

当郁闭度0.9以上时，通过试验此时正射影像图上的对角线基本上被树冠所覆盖，计算后的郁闭度接近于1但实际的郁闭度可能在0.9 ~ 1.0之间，因为该情况下无人机是拍不出地面和非树叶部分的，此时不能忽略树叶之间的空隙，所以此情况下误差将可能变大。

结合上述分析：本研究理想的郁闭度测量范围为0.5 ~ 0.9，在此条件下本研究测量的郁闭度比实际郁闭度偏大或接近，两者差值在0.01以内；当郁闭度<0.5时，本研究测量值比实际值偏大，两者差值在0.05以内；当郁闭度>0.9时，经过实验测定本文测量值接近于1，且两者差值在0.1左右，产生的误差较大。

4 结论与讨论

本研究结果表明，通过无人机拍摄使用图像处理技术来测量郁闭度是可行的，在0.9郁闭度以下具有较高的精度，其计算结果为偏大或接近。其中误差产生原因为：（1）树木的高度会产生误差。无人机上空拍摄时，对最上层的树冠拍摄的最为清楚，若下层树冠被上层树冠直接遮挡或者阴影遮挡，该图层就会被认为地块，会直接当作背景给处理，而实际上是这部分有树冠的。（2）样点的位置以及树木的位置也会带来误差。稀疏的林地分布不均匀，即对角线上分布的跟其他位置分布的差别很大，若稀疏的林地树冠都分布在对角线上则计算比实际偏大，相反则偏小，选样地时尽量避免这种情况，但是对于密集分布的林地则影响较小了。（3）无人机拍摄时最好避免阳光，尽量阴天拍摄，拍出的照片光线要均匀。（4）地面植被旺盛同时树木稀疏的地块，这种样张会很难区分前景跟背景，肉眼虽然能识别但是图像处理时树叶跟植被几乎是一样的绿色，前景跟背景容易混在一起，带来分割难度。而树木密集地块，虽说上层树冠遮挡住了植被，但是树叶空隙被下层绿色植被填满了，也会有误差，相对于稀疏地块来得小。所以根据上面分析得出理想样地的条件为：若为稀疏地块，则地面绿色植被要少，树木分布要均匀；若为树木密集地块，树高差距不能过大，郁闭度过大误差也会增大；拍摄条件光线要均匀，拍摄高度不能过高，最好相对于树高5 ~ 10 m。

[1] 李永宁，张宾兰，秦淑英，等. 郁闭度及其测定方法研究与应用[J]. 世界林业研究，2008（01）：40－46.

[2] 刘菊. 基于GIS的林地郁闭度测定方法研究[J]. 山西林业科技，2015（02）：33－35.

[3] 刘芳，杨广斌. 基于鱼眼照片的森林郁闭度快速提取方法研究[J]. 西南林业大学学报，2013（02）：71－74.

[4] 杜晓明，蔡体久，琚存勇. 采用偏最小二乘回归方法估测森林郁闭度[J]. 应用生态学报，2008（02）：273－277.

[5] 郭昱杉，刘庆生，刘高焕，等. 基于标记控制分水岭分割方法的高分辨率遥感影像单木树冠提取[J]. 地球信息科学学报，2016（09）：1259－1266

[6]王国权，周小红，蔚立磊. 基于分水岭算法的图像分割方法研究. 计算机仿真[J]，2009（05）：265－268.

[7] 杨恒伏，孙光田，祖伟. 基于HVS特性的自适应K近邻均值滤波算法[J]. 计算机工程与应用，2009（10）：179－181.

[8] 陆宗骐，梁诚. 用Sobel算子细化边缘[J]. 中国图象图形学报，2000（06）：5.

[9] 韩思奇，王蕾. 图像分割的阈值法综述[J]. 系统工程与电子技术，2002（06）：91－94，102.

[10] 陈延梅，吴勃英. 基于数学形态学的图像增强方法[J]. 哈尔滨工业大学学报，2006（06）：906－908.

[11] 金茜，汪伟. 基于改进控制标记符分水岭的牙体硬组织分割[J]. 计算机应用研究，2018（12）：8.

Canopy Density of Broad-leaved Forest by Watershed Segmentation Algorithm

YAN Yu1，WANG Yong-zhong2，YANG Lai-bang1，LOU Xiong-wei1

（1. School of Information Engineering, Zhejiang A & F University, Zhejiang Provincial Key Laboratory of Forestry Intelligent Monitoring and Information Technology, Hangzhou 311300, China; 2. Hangzhou Ganzhi Software Company of Zhejiang, Hangzhou 311300 ,China）

Four groups were photographed at different height and distance on 5 sample plots of broad-leaved forest with 50 m×70 m in Linan, Zhejiang province by unmanned aerial vehicle. They were made into orthophotoquad, which were pretreated by gray level and median filter, tree canopy was extracted by improved tag-control watershed segmentation algorithm, and compared with that by manual one. The experiment showed that algorithm had higher accuracy, and the error was about 5%. Canopy density was calculated by line transect method on canopy estrated, and the results demonstrated that the canopy density between 0.5 and 0.9 was accurate.

unmanned aerial vehicle; watershed segmentation algorithm; canopy density；broad-leaved forest

S757

A

1001-3776（2019）06-0053-010

10.3969/j.issn.1001-3776.2019.06.009

2019-03-26 ；

2019-09-30

浙江省科技重点研发计划项目（2018C02013)，浙江省科技计划项目（2017C02044）资助

严羽，研究生，从事计算机在农林业上的应用研究；E-mail:875082356@qq.com。

楼雄伟，副教授，博士，从事模式识别、林业信息化研究；E-mail:123000811@qq.com。