基于对称凸包和平滑轮廓的单木通透度计算*

张卫正 张伟伟 李灿林 万瀚文 张秋闻 刘 岩 金保华

(1.郑州轻工业大学计算机与通信工程学院 郑州 450002；2.郑州大学国际学院 郑州 450052)

树冠通透度指基于树木冠形的、透过树木冠层可见天空区域的系数，是评估森林健康和树木健康的重要指标之一(Svein，2009；Spiecker，1996；宋文龙等，2015)。现阶段，树冠通透度缺乏估测的标准方法，在实际作业中，观察者通常手工绘制树木冠层的二维形态和轮廓，估计透过冠层的可见天空区域占树冠区域的百分比，从而实现等级评估(Bond，2012)；也有研究人员通过参考不同通透度等级的照片(卡片)，人工判读通透度的大概程度或等级(Schomakeretal.，2007)。但这些方法均依赖树冠形态和叶量的人工视觉判断，受观察者视角和主观性影响较大，具有高度的可变性，重复性差；而且，由于每株树木的树冠形态和叶片分布各不相同，难以建立统一的判别标准。

将数字图像处理技术用于定量评估树冠通透度，具有快速准确、自动化程度高等优势，能够克服人工判别可变性大、重复性差及易受主观影响等局限。Zhang等(2007)采用分形维数对树冠的空间格局进行了研究；Córica(2009)分析了阳光透过树木冠层后在外墙上产生光斑的分布模式；Barthélémy等(2007)将树冠孔隙度指标与树冠通透度和树木发育阶段联系起来，但是该指标依赖特定环境，并需要根据树种进行校准；Clark等(2003)探讨树冠通透度估算方法，以“收缩包裹”的视角定义树冠光滑边界，提出了广义的冠层区域和梯度变化概念，但没有详细说明如何设置阈值以及阈值对树木分枝、嫩枝等的影响。

本研究基于树木数字图像进行单木结构特征识别与整合，不仅对树冠所形成内部区域的通透度进行计算，采用K-均值聚类算法将树冠中足够大的可见天空区域设定为大孔(通常是由于主要树枝的缺失造成的)，将较小的可见天空区域设定为小孔(大多是由于水分胁迫、叶子病变或过早退化造成的)，并计算大孔、小孔的分布密度，而且还对树冠外部区域的通透度进行研究，当可见天空区域位于树冠轮廓外部时，探索可见天空区域凹陷进入树冠的深度，对树冠边界轮廓进行平滑处理，建立树冠的对称凸包，计算平滑轮廓距凸包的距离和可见天空区域的凹陷深度，通过对距离设置阈值，计算深度凹陷密度。基于树冠中的大孔和小孔密度、树冠轮廓的深度凹陷密度3个可量化指标计算树冠通透度。

本研究的创新之处在于：1) 将2个原始轮廓定义为参考形状，即用于评估树冠深度凹陷的对称凸包，以便于计算小孔和大孔的树冠平滑轮廓；2) 采用K-均值聚类算法，将计算得到的平滑轮廓上各点到凸包最近距离的集合自动分为2类，计算类间阈值(查找小的类中的最大值与大的类中的最小值，计算二者之和的均值作为阈值)，通过阈值区分轻度和深度凹陷，并统计树冠中的大孔和小孔；3) 根据对树冠通透度的贡献，赋予深度凹陷密度、大孔和小孔密度3个参数不同权重。基于以上创新之处，计算单木通透度，以期为单木的健康状况监测和生长状态分析提供技术支持。

1 图像采集与预处理

图像采集与预处理流程如图1所示。

图1 单木通透度计算流程Fig.1 Flowchart for calculating the single tree permeability

1.1 图像采集 采用微软公司开发的便携式平板电脑Surface Pro 4获取单木图像。该设备触控显示屏尺寸为12.3英寸，分辨率为2 736×1 824像素，采用第6代Intel酷睿处理器，系统内存为16 G，具有800万像素后置自动对焦摄像头，配备1 024级压感触控笔，电池续航时间长达9 h，机身厚度8.45 mm，质量786 g，性能强大又不失轻盈，携带方便，便于移动操作，满足本研究所需图像拍摄和图像处理要求。

以雪松(Cedrusdeodara)为研究对象，该树种为常绿乔木，树冠尖塔形，大枝平展，小枝略下垂，在我国多地栽培引种，除了可用作庭园观赏外，也是一种重要的建筑用材。采用设备自带的高清摄像头获取单木图像，图像中可清晰辨识枯枝、新枝等反映树木生长状态的细节信息。在采集图像时，确保观察者处于合适的距离，以相机可拍摄树木全貌且图像显示清晰为准，同时尽量保持仰视角小于60°、俯视角小于30°，水平距离大于树高，如图2所示。标记观察者获取图像时站立的位置、树木种类、天气情况、拍摄日期等信息，便于后续不同时段进行重复观测，分析和评估单木通透度的变化。

图2 拍摄单木图像示意Fig.2 Schematic diagram of taking an image of single tree

1.2 单株树木提取 考虑到实际采集环境，目标树木两侧和背景中可能有建筑、车辆、相邻粘连交叉树木等干扰，会对单木图像分割和提取产生影响，而常用的图像分割和提取方法难以达到理想效果，故本研究利用Surface Pro 4自带的压感触控笔手工圈存图像中冠层区域，以消除图像中的无关信息，增强单木冠层的可检测性并最大限度简化处理数据的复杂度，从而改进特征提取、图像分割、匹配和识别的可靠性(刘博峰，2016)。树冠大都包含树枝和树叶等诠释树木主要形态特征和生物学特性的组件，本研究针对树冠进行通透度分析，排除了树干。

采用Matlab中的函数get(axesHandle,‘CurrentPoint’)获取鼠标点击处的图像坐标，并将这些坐标组成封闭多边形。为了提取多边形所包围的感兴趣区域(即单木区域)，使用函数inpolygon(m1,n1,im,in)判断图像中各像素点是否在多边形内，其中im、in为构成多边形边界的顶点坐标，m1、n1为图像中的像素点坐标。get函数返回结果为逻辑类型的0或1，如果该点在多边形内则返回1，否则为0。运用imwrite()函数将多边形包围的图像区域存储，区域外的图像背景填充为黑色，可实现手动点击产生多边形，并圈存多边形区域，完成单木区域提取，如图3所示。

图3 手工圈存目标树木Fig.3 Artificial circles target treea.原始图像Original image；b.人工圈存轨迹Manually tracing trajectory；c.存储树木区域格Storage tree area grid.

1.3 感兴趣区域的灰度化和二值化 将圈存的单木区域作为感兴趣区域进行图像灰度化和二值化，以便于后续处理。采用平均值法(冯志新等，2013)将RGB图像转换为灰度图像：

Gray=(R+G+B)/3。

(1)

运用最大类间差法(Ostu算法)确定的阈值进行图像二值化，该算法不受图像亮度和对比度的影响，按图像的灰度特性，将图像分成背景和前景2部分(廉宁等，2014)。设灰度图像灰度级为L，则灰度范围为[0,L-1]，运用Ostu算法计算图像的最佳阈值为：

σ=Max{w0(t)×[u0(t)-u]2+w1(t)×
[u1(t)-u]2)}。

(2)

式中：w0为背景比例；u0为背景均值；w1为前景比例；u1为前景均值；u为整幅图像的均值；σ取最大值时的分割阈值t即为最佳阈值T。

使用阈值T对原图像f(x,y)进行分割，将图像分为2部分，g(x,y)即为通过阈值处理得到的二值化图像：

(3)

通过对圈存区域进行图像灰度化和二值化，得到树木二值化图像，如图4a、b所示。

对图3c中圈存树冠区域的绿色通道进行二值化，得到图4c，并将4b与c合并，最大限度保留树冠形态，得到图4d。

图4 感兴趣区域二值化Fig.4 Region of interest binarizationa.图像灰度化Grayscale image；b.图像二值化Image binarization；c.绿色通道二值化Green channel binarization；d.合并b与c Merge b and c.

将图4d中的圈存区域取反，得到图5a。采用Matlab中的imclose函数对图5a进行闭运算，得到图5b，其中闭运算采用直径为2的圆形结构元素D1。闭运算对图像先膨胀后腐蚀，具有填充物体内细小空洞、连接邻近物体和平滑边界的作用。

采用Matlab中的bwlabel函数对图5b进行连通区域标记，并保留面积最大的连通区域，得到图5c。

图5 单株树木的二值化图像Fig.5 Binarized images of single tree a.取反Negate；b.闭运算Closed operation；c.保留面积最大的连通区域Keep the connected area with the largest area.

1.4 树冠平滑轮廓 首先，采用Matlab中的imfill函数对树冠二值化图像进行孔洞填充，如图6a所示。然后，采用半径为2的圆形结构元素对图6a进行闭运算，实现树冠区域的边界平滑，运用bwboundaries函数提取树冠边界点，得到树冠平滑轮廓，如图6b、c所示。

图6 树冠平滑轮廓Fig.6 Smooth outline of the canopy a.孔洞填充Hole filling；b.平滑轮廓Smooth outline；c.显示平滑轮廓Display smooth outline.

2 图像处理

2.1 对称轴确定 正常生长树木的冠形在竖直方向基本保持平衡和对称，其对称轴通常是树木主干所在的直线。确定树冠在竖直方向的对称轴，建立树冠基于对称轴的镜像，构建对称树冠凸包，从而有助于后续与通透度相关参数的计算。

针对已经手工圈存的单木图像，其大小为m×n(m为行数，n为列数)。设定对称轴所在的列为p，则对称轴左侧的树木像素[1～ (p-1)列的所有树木像素]之和与对称轴右侧的树木像素[(p+1)～n列的所有树木像素]之和的差值最小。在Matlab中计算出p，并在图像中画出对称轴，如图7所示。

图7 树冠的对称轴Fig.7 Symmetry axis of the canopya.树冠的对称轴The symmetry axis of the canopy；b.在全景图像中显示对称轴Show the symmetry axis in the panoramic image.

2.2 树冠对称镜像和对称凸包确立 计算树冠通透度时，考虑树冠中的分枝、嫩枝、新枝、叶片以及可能缺失的枝叶，以对称轴为中心左右两侧的树冠分别做镜像，然后合并成为新的树冠，如图8b所示，白色像素为原始树冠，灰色像素为镜像。

图8 对称树冠与凸包Fig.8 Symmetrical crown and convex hulla.树冠Canopy；b.对称树冠Symmetrical canopy；c.对称树冠的Delaunay三角网Delaunay triangulation of symmetrical canopy； d.对称树冠凸包Convex hull of symmetrical canopy.

将图像中表示对称树冠的每个像素当成一个点，每个点都具有行坐标和列坐标，所有这些点构成集合X。先用DelaunayTri函数得到集合X的三角网，后用Convex hull函数得到树冠凸包。

采用Matlab中的DelaunayTri函数建立约束三角网，将集合X中的所有点进行规则化处理，使得集合中每个点位于三角形的顶点且三角形构成的网络具有空圆特性和最大最小角特性(余杰等，2010；Zhanetal.，2017)。采用Matlab中的Convex hull函数获得构成凸包的边界点(Yangetal.，2014)，依次连接凸包边界点构成树冠凸包，如图8d所示。

2.3 深度凹陷密度计算 为了衡量树冠伸展情况，采用凹陷深度表示树冠填充凸包的紧凑程度。深度凹陷密度指树冠平滑轮廓外部可见天空区域凹陷进入树冠平滑轮廓较深的区域面积与凸包面积的比值。

将树冠平滑轮廓和对称树冠凸包添加到树冠区域，如图9a、b所示。图9c显示了构成树冠平滑轮廓和凸包的点，并将其中一部分(黑色矩形框)放大显示。从平滑轮廓中黑色圆点处开始沿平滑轮廓顺时针方向行走，计算平滑轮廓上各点到凸包的最短距离，并画出距离图，如图10a所示。

图9 树冠的对称凸包及平滑轮廓Fig.9 Symmetrical convex hull and smooth outline of canopya.平滑轮廓和凸包Smooth outline and convex hulls；b.原始图像上显示平滑轮廓和凸包Show smooth outline and convex hull on the original image；c.放大显示局部的平滑轮廓和凸包Zoom in to show local smooth outline and convex hull.

采用K-均值聚类算法(颜佩等，2017)确定深度和轻度凹陷阈值，当凹陷深度超过此阈值即设定为深度凹陷。对距离图进行K-均值聚类(K=2)处理获得阈值Ddepression，如图10a中红色线段所示。

将位于树冠平滑轮廓以外、凸包以内的区域作为凹陷区域，如图10b中的白色区域所示。计算该凹陷区域内像素点到凸包边界的最短距离，如果最短距离大于Ddepression，则该点为深度凹陷点，所有深度凹陷点构成深度凹陷区域。遍历凹陷区域得到所有深度凹陷点，如图10c中灰色区域、图10d中红色区域所示。

图10 深度凹陷区域确定Fig.10 Establishment of the deep depression areaa.平滑轮廓上各点到凸包的最短距离The shortest distance from each point on the smooth outline to the convex hull；b.白色区域为凹陷The white area is the depression；c.灰色区域为深度凹陷The gray area is the depth depression；d.红色区域为深度凹陷The red area is the depth depression.

2.4 大小孔检测及密度计算 将树冠平滑轮廓内的区域(狭义的树冠区域)进行连通区域标记，采用K-均值聚类算法(设定为2类)将标记的连通区域分为大孔和小孔，所有标记为小孔的连通区域像素总数除以树冠平滑轮廓所包围的像素总数得到小孔密度，然后计算大孔密度。

采用Matlab中的连通区域标记函数bwlabel对树冠中所有可见天空区域进行标记，统计其像素数，并对所有连通区域的像素数进行K-均值聚类，其中K为2(树冠中所有可见天空区域分为小孔和大孔，其中小孔对应类别1，大孔对应类别2)。

树冠中树干、树枝和树叶的值设为0，小孔的值设为0.333 3，大孔的值设为0.666 6，树冠轮廓外部像素的值设为1，如图11所示。

图11 树冠的小孔和大孔Fig.11 Small and large holes of canopy

根据所提取的树冠、大孔和小孔，计算得大孔密度Dmacro和小孔密度Dmicro分别为0.160 5和0.179 2。

2.5 树冠通透度计算 参考林业专家经验和树木生长规律，采用深度凹陷密度、大孔和小孔密度3个系数定量评估单木通透度Tc：

(4)

深度凹陷和大孔是因树枝缺失造成的，考虑深度凹陷对冠形和通透度的影响高于大孔，对深度凹陷赋予更大加权。此外，式(4)中也加入了小孔对通透度的贡献。

2.6 多视图聚合 从0°、30°、60°、90°、120°和150°共6个角度获取6个树冠图像，分别计算通透度系数Tc1、Tc2、Tc3、Tc4、Tc5和Tc6，以平均值Tca作为单木通透度系数，减小因视角变化引起的通透度波动，尽可能精确反映单木真实状况。如果难以从多个视角获取树木图像，则根据实际获取数量计算平均值作为通透度系数：

(5)

3 系统实现与分析

3.1 系统实现 由Matlab GUI (graphical user interface，为图形用户界面，又称图形用户接口)设计开发的系统界面如图12所示，依次点击打开图像、提取树木、确立对称轴、凸包与平滑轮廓、深度凹陷、大小孔检测等计算单幅图像通透度，分别显示在右下方的表格中；也可在提取冠层后，直接点击“显示”按钮计算通透度。

图12 系统界面Fig.12 System interface

本研究对其他雪松进行图像采集和通透度分析，如图13所示。

图13 实例Fig.13 Example

3.2 树冠对称镜像和对称凸包确立 为了对本研究所提出方法进行试验验证和精度分析，设计一个单木验证模型，如图14a所示。

图14 单木验证模型Fig.14 Single tree verification model a.本研究所设计的单木验证模型The single-wood verification model designed by this research institute；b.人工判读识别验证模型并标示大孔、小孔及深度凹陷区域Manual interpretation to identify the verification model and mark the large holes,small holes and deep depression areas；c.将验证模型显示在液晶平面显示器 上Display the verification model on the LCD flat panel display.

本研究所提出方法对尺度不敏感(目标树等比例缩放对研究结果无影响)，计算出的树冠通透度为一比值，故采用AutoCAD 2010结合Photoshop CC 2017设计单木验证模型(尺寸无需单位)。组成树冠3个等腰直角三角形的腰长为4，从上到下垂直排列；树冠中大孔是边长为1的正方形，小孔是腰长为1的等腰直角三角形。采用人工方法确定该模型凸包，并人工计算深度凹陷阈值，画出深度凹陷区域，如图14b所示。计算树冠中小孔、大孔、树冠区域、凸包和深度凹陷区域的面积分别为3、3、24、40和6.514 8，理想情况下的小孔密度、大孔密度、深度凹陷密度和通透度分别为0.125 0、0.125 0、0.162 9和0.264 6。

将单木验证模型显示在液晶平面显示器上，屏幕尺寸为23.8英寸，分辨率为1 920×1 080像素，垂直和水平可视角度均为178°。采用Surface Pro 4拍摄该显示器，获得单木模型图像。利用本研究所提出方法计算的小孔密度、大孔密度、深度凹陷密度和通透度分别为0.117 8、0.124 1、0.164 0和0.258 6，由此得到本研究所提出方法的精度为97.73%(0.258 6 / 0.264 6)。

3.3 深度凹陷密度计算 采用K-均值聚类算法自动实现大孔和小孔分类以及深度和轻度凹陷分类。比较2个类的相似度，满足以下条件，则认定2个类可以合并：

Mhigher<1.4×Mlower。

(6)

式中：Mhigher为较大的类均值；Mlower为较小的类均值。

当进行凹陷深度计算时，如果出现2个类合并，则认为树冠只呈现轻度凹陷，此时深度凹陷对通透度系数的贡献为零。同样，大孔和小孔均值满足式(6)，则全部显示为小孔，此时大孔对通透度系数的贡献为零。

如果多角度采集单木图像时背景存在黏连、遮挡等问题，导致难以差距6个角度的单木图像，则式(5)中Tca应以实际采集的图像数量为准。采用手工圈存树冠区域方法，未能做到全自动分析，但可减小树冠周围复杂背景的干扰，且人工工作量小。如果被检测树木的背景为其他树木或暗色干扰物时，在人工圈存树冠区域后，可采用He等(2011)提出的基于暗通道的去雾算法进行图像增强，从而提高对冠层大孔和小孔的识别效果。

4 结论

为解决现有单木通透度评估受观察者视角和主观性影响较大、难以建立统一判别标准的问题，本研究采用Surface Pro 4获取单木图像，进行单木结构特征识别与整合，通过建立树冠对称凸包和树冠平滑轮廓，不仅对树冠所形成的内部区域进行研究，标记内部的可见天空区域，并统计各连通区的像素数，采用K-均值聚类算法将各连通区域分为大孔和小孔，而且还对树冠轮廓外部可见天空区域凹陷进入树冠的深度进行探索，采用K-均值聚类算法判别轻度和深度凹陷。通过自动量化大孔和小孔密度、深度凹陷密度3个指标计算单木通透度系数，并对过程和结果进行详细分析。结果表明，本研究所提出方法的精度高达97.73%，测量速度快、人工工作量小，可为单木健康状况监测和生长状态分析提供技术支持，同时该研究思路和方法也可以推广应用到其他树木和作物的监测分析，具有一定的实用价值。