基于M- K聚类法的果树上下冠层体积比测算

祁力钧 程一帆 程浈浈 杨知伦 吴亚垒 葛鲁振

(中国农业大学工学院， 北京 100083)

0 引言

果树冠层体积是决定施用农药量的重要指标，在果园精细化作业中具有十分重要的地位[1-3]。随着科技的发展，基于高新传感器的测量技术成为冠层信息的主要获取手段[4-8]。

近年来，由于无人机自身所具有的不受空间限制、自由度高与经济成本低等优点，借助机载LiDAR回传数据，根据ABA(区域分析法)和ITC分析法(单株分析法)[9]对树木冠层进行测量成为农林测绘领域的研究热点。对于ITC法而言，由于无人机的独特视角导致来自机载LiDAR较大比例的激光脉冲被拦截在冠层的顶部，而中部和下部的反馈较少[10]，单纯依靠机载LiDAR回传数据进行冠层测算忽视了林木下冠层体积，致使结果不精确[11]。针对上述问题WEZYK[12]提出了将地面激光扫描(TLS)回传数据和机载激光扫描(ALS)点云转换为同一个坐标系的方法，并对林木特征参数进行半自动提取；MURGOITIO等[13]尝试将ALS和TLS传感器的回传数据整合一起，用于可视化树冠被遮挡部分，其结果表明TLS对ALS数据进行补充可显著完善森林结构模型；LOVELL等[14]将地基激光扫描仪与无人机机载系统的测量数据相结合，以减轻树木上冠层的遮挡对体积测量结果产生的影响。然而，用TLS的方法来获取果树下冠层体积，成本很高、耗时较长且需多点位布控的复杂测量方式不适于果园种植紧密、行列间空隙较少的特殊情况。

为解决上述问题，本文参照LEFSKY等[15]的冠层分割思想，拟将果树分为上下2个冠层，在地面采集果园果树侧视图图像后运用基于M- K聚类法的图像处理技术对冠层进行提取，并在图中定位最大遮挡面[16]，以此为基准实现果树上下冠层的分割；针对密闭型和疏散型两类果树，分别运用散点积分与旋转积分的方法，从而得到果树上下冠层的体积比。将此参数引入到机载LiDAR系统中，对冠层信息缺失的部分进行体积预测，以实现多传感器信息融合。

1 系统和材料

1.1 检测系统组成

本实验检测系统由硬件与软件两部分组成，其中硬件由微型单反相机、云台、三脚架、计算机等组成，如图1所示。相机为2 430万有效像素的Sony α7微型单反相机，搭载15～75 mm变焦镜头，计算机选用Accer aspire V5，Intel core i5处理器，Windows 10操作系统。软件为基于Matlab语言的自编程序。

图1 检测系统的硬件组成Fig.1 Hardware composition of inspection system

1.2 果树图像获取

分别于2017年5月5日10:00在北京市中农富通现代示范果园选择20棵自然纺锤形樱桃树和2017年9月20日10:00在山东省果树研究所示范果园选择23棵高纺锤形苹果树进行实验。

在自然曝光的环境中，运用Sony α7微型单反相机的AUTO模式进行图像采集，采集时相机装载在百诺IT15云台上。

2 冠层体积比的测算

为求得果树冠层上下体积比，通过图像处理法提取冠层图像相关信息，并根据提取的数据进行体积比计算，最后为增加结果的准确性与鲁棒性提出修正方案。

2.1 果树图像分割

图像分割是根据一定特征将图像中目标区域与背景区域进行分离、提取的技术过程。对目标物分割的最终效果将会直接影响到后续的图像边缘提取与最大遮挡面的定位。果树冠层相对于背景区域具有明显的颜色差异，根据该特点，提出运用M- K聚类法进行图像分割。

2.1.1马氏距离

马氏距离是由MAHALANOBIS于1936年提出的一种基于变量间相关性的距离度量。相对于欧氏距离，马氏距离的优势是考虑了数据之间的相关性，且在进行距离计算时具有尺度无关性，它是计算2个未知样本集相似程度的一种有效方法[17]。

本文运用LI等[18]的方法，计算所采集果树冠层图像中每一像素点与经人工分割的标准果树图像之间的马氏距离，并将输出矩阵MDi作为二者相似度的判定标准。其计算公式为

(1)

(2)

式中xi——采集图像中某一像素点包含的颜色信息向量

n——图像的总像素数

X——图像所包含的颜色特征矩阵

图2为经上述公式计算后生成MDi矩阵的三维显示图，其中Z轴为马氏距离，X、Y轴为像素点所在行、列值，顶部平面图像为马氏距离越大，该点对应亮度越高，距离越小亮度越低。由图2可知，图像中部灰色凹陷区域与果树标准图像相似度较高，马氏距离较小，为果树树冠主体，其余部分为地面、天空等非目标区域，与果树标准图像相似度较低，马氏距离较大。

图2 MDi矩阵的三维显示图Fig.2 Three-dimensional display of MDi matrix

2.1.2色度空间的变换

由于园内光照环境复杂，导致所得冠层图像通常存在亮度不均匀的现象。在进行图像处理时，目标区域的亮度差异成为干扰因子，影响了分割的准确性。为排除亮度对图像分割的影响，本文将果树图像由RGB色彩空间转换为Lab色彩空间，转换公式为

(3)

(4)

(5)

(6)

图3 K-means法与M- K法分割效果对比Fig.3 Comparisons of K-means and M- K segmentation effects

(7)

式中r、g、b——图像红、绿、蓝通道灰度

L、a、b——转换后的色彩空间三通道灰度

在Lab色彩空间中L表示亮度，a表示从绿色到红色的色彩范围，b表示从蓝色至黄色的色彩范围[19]，该色彩空间将图像颜色因子与亮度因子进行分离，能够有效解决图像光照不均匀的问题，从而保留图像原本的色彩信息。

2.1.3M- K聚类法

常用的聚类算法有K-means法、模糊c均值法、meanshift法和gmm法等，其中K-means聚类算法运算速度快、结果准确。它的基本思想是：首先从数据对象中随机选择k个对象作为初始聚类中心,然后将剩余的每个对象根据与这些聚类中心的距离,分别赋予与其距离最近的聚类。再重新计算每个新聚类的聚类中心,不断重复,直到聚类中心不再变化[20]。然而对于现代果园内环境复杂、干扰因素众多等特殊情况，K-means聚类法存在错分漏分现象严重、聚类数k难以确定的缺陷。

针对上述问题，本文提出结合马氏距离和K-means算法的M- K聚类法，并运用该算法对果树冠层图像进行分割。首先将前文所得马氏距离MDi矩阵与Lab色度空间中的a、b矩阵进行归一化处理(归一到0～255)，将此3类特征数据结合成新的三维矩阵，以新的MAB矩阵替代原有的RGB矩阵进行K-means聚类。MAB图像中M代表该像素点与目标图像整体的相似程度，A、B包含原有的图像色彩信息。新的三维矩阵排除了干扰因子亮度L的影响，引入了相似度M作为聚类特征，在未损失图像原本的颜色信息A、B的前提下，使得聚类分割的准确度得以提高。如图3所示，对同一幅果树图像分别进行基于K-menas聚类及基于M- K算法的冠层分割，从结果来看，运用M- K法进行分割，冠层完整度较高，错分现象不明显。

图4 马氏距离的空间散点图Fig.4 Space scatter plots of Mahalay distance

由于背景区域(地面、天空)与果树冠层有一定的差异(图4b)，M往往较大，而目标区域与干扰区域(杂草、其他品种果树)可通过颜色信息A、B进行区分(图4c)，故在进行M- K聚类时，统一将聚类数设置为k=3，将图像分为目标区域、背景区域和干扰区域3部分，从而解决K-means算法无法确定初始聚类数的问题。

2.1.4图像形态学处理及边缘轮廓提取

对果树图像进行M- K聚类后，取M最小的聚类中心所在类为目标类，进行二值化。由于树冠中存在空隙，分割所得二值图像存在孔洞现象，且仍有少数干扰物未能完全去除。因此运用形态学孔洞填充和开闭运算进行处理，以达到消除孔洞与残留干扰物的效果，结果如图5c所示。

传统的边缘提取算法有Roberts算法、Sobel算法和Canny算法[21]，其中以Canny算法的检测效果最为优良。对经过形态学处理的二值图像进行基于Canny算法的边缘提取，获得果树冠层边缘轮廓，如图5d所示。

图5 果树冠层边缘轮廓提取Fig.5 Extractions of canopy edge contours of fruit trees

2.2 体积比计算

2.2.1冠层最大遮挡面的定位与图像切割

进行上下冠层体积比的计算时，首先需要确定其最大遮挡面在树冠中所处的位置。在果树图像处理中，一般认为图中最大横向像素距离为果树最大冠幅[22]，从而近似认为其所在位置即为最大遮挡面所处位置。然而，现实中果树枝叶不齐，某些较长的侧枝外伸，破坏了树木的轮廓走势，因而简单的将最大冠幅等同于最大遮挡面容易使测量结果产生较大的偏差。

为解决上述问题，本文首先对冠层图像进行遍历，寻找出果树的顶点，过此顶点将图像竖直切割成两份；然后对左右两部分果树的边缘轮廓进行基于最小二乘法的多项式曲线拟合，如图6a所示，并根据留一交叉验证法进行循环检验，自动选取最佳拟合函数，输出拟合曲线方程，分别解出两曲线极大值所在的位置，定为果树遮挡面，再次进行图像切割，将左右两图切成左上、左下、右上、右下4份，如图6b所示。

图6 最大遮挡面的定位与图像切割Fig.6 Positioning of the largest occlusion surface and cutting of image

以较为平滑的拟合曲线代替果树不规则的轮廓线，排除了过长侧枝的干扰，提高了冠层最大遮挡面定位与分割的稳定性。

2.2.2上下冠层体积比

对于边缘疏散形果树，由于树叶稀疏，分支明显，中空面积较大，运用散点积分的方法求体积易产生误差，所以采用树冠轮廓拟合曲线旋转体积代替原果树体积，具体方法如下：首先运用骨架化法对果树冠层进行迭代骨架化运算，对经过骨架化的图像进行八连通邻域检测，得出树枝的端点，对所得端点进行曲线拟合，如图7所示。再对曲线进行旋转积分。其计算公式为

(8)

式中f(x)——所得拟合曲线

x——拟合曲线横坐标值，像素

l——图像的列像素数

图7 疏散型果树的曲线拟合Fig.7 Curve fitting of evacuation fruit trees

对于边缘轮廓规则的密闭型果树，直接运用散点积分法求果树各部分冠层的像素体积，其计算公式为

(9)

式中r——树冠边缘轮廓与分割中心的距离，像素

最后，根据求得的各部分像素体积进行上下冠层体积比P的计算，公式为

(10)

式中V11——左上部果树的像素体积

V12——右上部果树的像素体积

患者右骶棘肌紧张，腰椎MRI示L4 ～ S1区域增粗变宽的右骶棘肌无信号异常，排除损伤可能。MRI异常影像结合临床症状，分析疼痛原因如下。

V21——左下部果树的像素体积

V22——右下部果树的像素体积

2.3 修正方案

为提升测量结果的精确性与鲁棒性，选择在室内对两棵特征参数已知的果树(樱桃树、苹果树)模型进行探究实验，确定对测量产生影响的外界因素(图像拍摄距离、高度和方向)，并针对产生的影响提出相应的修正方案。

2.3.1影响因素

实验于实验室内进行，实验对象为一棵苹果树模型与一棵樱桃树模型。将相机架设在三角架上，利用水平仪进行校正，使相机镜头主光轴与地面保持平行。如图8所示，4次实验中，分别控制相机对地高度为180、150、120、90 cm；镜头与树干之间的距离为260、310、360、410 cm；每次实验拍摄2个方向，拍摄方向依次为0°、45°、90°、135°、180°、225°、270°、315°(以正东方向为0°)。

图8 实验示意图Fig.8 Experimental diagrams

以拍摄距离和高度为变量，对两棵模型树的补偿体积系数进行测量，结果如表1所示。对结果分别进行方差分析，如表2所示。经F检验得出，在置信度为95%的前提下，改变拍摄距离、拍摄高度对实验结果的影响不明显。

表1 补偿体积系数测量结果Tab.1 Measurement result of compensationvolume coefficient

表2 方差分析Tab.2 Variance analysis result

以拍摄方向为变量进行实验，结果如图9所示，对所得结果进行方差分析，发现置信度为90%的情况下拍摄方向对测量结果产生显著影响。

图9 误差折线图Fig.9 Error lines graph

2.3.2修正方案

在改变拍摄方向的实验中，误差如图9所示，樱桃树所得补偿体积系数与人工测量的最大误差为22.2%，苹果树与人工测量值的最大误差为25.3%，可见最大偏差均在20%以上，结果不够精确。原因是式(8)、(9)中将修形后的果树当作规则的旋转体进行体积计算，然而在实际情况中，由于各种偶然因素的影响，果树很难满足理想条件，当各部分生长差异过大时，单侧果树图像不足以反映整棵果树的形体特征。为了使测量结果足够精确，需采集多幅不同方向的果树图像，以确保获取足够全面的果树信息。

分别对苹果树与樱桃树不同侧面的体积测量值进行取平均，所得结果与人工测量值误差较小(苹果树10.7%，樱桃树3.9%)，因此在采集冠层图像时应根据实际情况采集2幅及以上不同方向的果树图像，以平均值作为最终结果，达到降低测量误差的目的。

3 实验结果与分析

选取20棵樱桃树、23棵苹果树作为研究对象。对每棵果树分别采集东北侧和西北侧的果树图像进行上下冠层体积比计算。在确定冠层最大遮挡平面后，参考WHEATON等[23]、王佳等[24]提出的计算方法，运用卷尺、标杆和手持激光测距仪测量果树冠层各截面所在高度与半径(将各截面近似为圆形)，运用圆台累加法对果树上下冠层的真实体积进行人工测算，以人工测量所得上下冠层体积比作为标准值，与本文所述方法进行对比分析。

分别将两棵果树的东北侧、西北侧摄影测量值与人工测量值进行线性相关分析，如表3所示，得出的决定系数R2最低为0.661，最高为0.818，可见运用单幅果树图像测算上下冠层体积比其结果与人工测量值具有一致性，但不同树种、不同侧面的测量值与人工测量值的线性相关性有较大偏差(23.7%)，方法稳定性较差。对两棵果树东北侧与西北侧的测量值进行算术平均后与人工测量值进行相关分析，结果显示相对于未进行平均之前二者的相关关系具有显著增强，其决定系数分别为0.775和0.832；再进行显著性t检验，P值为0.389 8、0.613 9，均大于0.05，可见本文方法测量值与人工测量值未见显著差异，且2次回归的决定系数差异缩小(7.4%)，说明经修正后稳定性增强。

表3 果树冠层体积的人工测量值与图像处理值Tab.3 Artificial measurements and image processing values of canopy volume

对比图10a和图10b，可见图10a中与回归线偏差较大的奇异点分布较为均匀且数量多，而图10b中奇异点多出现于上下冠层体积比较大的区间且数量较少。这主要是由于秋天的苹果树普遍叶密度较小且树形不够规整，在进行图像处理时存在孔隙过大、枝干明显等问题，难以形成闭合的树冠，导致处理结果存在一定的偏差；而春天的樱桃树枝叶茂密，但幼龄果树未完全发育成型，其下冠层所占体积比例低，上下冠层的界限不够明显，难以寻找合适的最大遮挡面，从而导致测量产生偏差。

图10 线性回归Fig.10 Linear regressions

4 结论

(1)在冠层分割中，对于任意图像M- K算法均将聚类值确定为3，且引入了MDi矩阵作为图像特征值，相对于单一的K-means分割法，M- K算法解决了聚类数k无法确定的缺陷，提高了分割的准确性。

(2)实验结果表明，改变拍摄方向会对测量结果产生显著影响。对原方案经过修正后，误差由最高的25.3%降至10.7%。以23棵苹果树、20棵樱桃树为实验对象，对本文方法测量值与人工测量值进行线性分析和显著性检验，结果表明二者之间有较强的线性相关关系(R2为0.775、0.832)和较小的差异(P为0.389 8、0.613 9)。

1 王万章，洪添胜，李捷，等. 果树农药精确喷雾技术[J]. 农业工程学报,2004,20(6):98-101.

WANG Wanzhang, HONG Tiansheng, LI Jie, et al. Review of the pesticide precision orchard spraying technologies[J]. Transactions of the CSAE, 2004, 20(6): 98-101. (in Chinese)

3 吴亚垒，祁力钧，张亚，等.基于驻波率原理的农药雾滴沉积量检测系统设计与试验[J].农业工程学报，2017，33(15)：64-71.

WU Yalei, QI Lijun, ZHANG Ya, et al. Design and experiment of pesticide droplet deposition detection system based on principle of standing wave ratio[J]. Transactions of the CSAE,2017, 33(15): 64-71.(in Chinese)

4 丁为民,赵思琪,赵三琴,等.基于机器视觉的果树树冠体积测量方法研究[J/OL].农业机械学报,2016,47(6):1-10,20.http:∥www.j-csam.org/jcsam/ch/reader/view_abstract.aspx?flag=1&file_no=20160601&journal_id=jcsam. DOI：10.6041/j.issn.1000-1298.2016.06.001.

DING Weimin,ZHAO Siqi,ZHAO Sanqin, et al.Measurement methods of fruit tree canopy volume based on machine vision[J/OL].Transactions of the Chinese Society for Agricultural Machinery,2016,47(6):1-10,20. (in Chinese)

5 何诚，冯仲科，袁进军，等.基于数字高程模型的树木三维体积测量[J].农业工程学报，2012,28(8):195-199.

HE Cheng,FENG Zhongke,YUAN Jinjun,et al. Three-dimensional volume measurement of trees based on digital elevation model[J]. Transactions of the CSAE，2012，28(8):195-199．(in Chinese)

6 刘峰，龚健雅.基于机载激光雷达技术的茂密林地单株木识别[J]. 农业机械学报, 2011, 42(7):200-203，209.

LIU Feng，GONG Jianya．Individual trees recognition in dense forest based on airborne LiDAR[J]． Transactions of the Chinese Society for Agricultural Machinery，2011，42(7): 200-203，209．(in Chinese)

7 SCHUMANN A W,ZAMAN Q U.Software development for real-time mapping of ultrasonic tree canopy size[J]. Computers and Electronics in Agriculture,2005,47(1):25-40.

8 王万章，洪添胜，陆永超，等.基于超声波优越感器和DGPS的果树冠径检测[J].农业工程学报，2006，22(8)：158-161.

WANG Wanzhang, HONG Tiansheng, LU Yongchao, et al.Performance of tree canopy diameter measurement based on ultrasonic sensor and DGPS[J].Transactions of the CSAE,2006,22(8):158-161.(in Chinese)

9 BREIDENBACH J, N☞SSET E, LIEN V, et al. Prediction of species specific forest inventory attributes using a nonparametric semi-individual tree crown approach based on fused airborne laser scanning and multispectral data[J]. Remote Sensing of Environment, 2010, 114(4):911-924.

10 CHASMER L, HOPKINSON C, TREITZ P. Assessing the three-dimensional frequency distribution of airborne and ground-based lidar data for red pine and mixed deciduous forest plots[J]. International Archives Photogramm Remote Sens Spat Information: Science, 2004, 36(8): W2.

11 HILKER T,LEEUWEN M V,COOPS N C,et al. Comparing canopy metrics derived from terrestrial and airborne laser scanning in a Douglas-fir dominated forest stand[J]. Trees-Structure and Function,2010,24(5):819-832.

12 WEZYK P. The integration of the terrestrial and airborne laser scanning technologies in the semi-automated process of retrieving selected trees and forest stand parameters[J]. Ambiência, 2012, 8(Especial):67-80.

13 MURGOITIO J, SHRESTHA R, GLENN N, et al. Airborne LiDAR and terrestrial laser scanning derived vegetation obstruction factors for visibility models[J]. Transactions in GIS, 2014, 18(1):147-160.

14 LOVELL J L, JUPP D L B, CULVENOR D S, et al. Using airborne and ground-based ranging lidar to measure canopy structure in Australian Forests[J]. Canadian Journal of Remote Sensing, 2003, 29(5):607-622.

15 LEFSKY M A, COHEN W B, ACKER S A, et al. Lidar remote sensing of the canopy structure and biophysical properties of Douglas-Fir Western Hemlock Forests[J]. Remote Sensing of Environment, 1999, 70(3):339-361.

16 JENNINGS S B, BROWN N D, SHEIL D. Assessing forest canopies and understorey illumination: canopy closure, canopy cover and other measures[J]. Forestry, 1999,72(1):59-74.

17 骆志高，李旭东，赵俊丽，等.利用马氏距离判别法准确实现对裂纹的识别[J]. 振动与冲击,2013,32(21):186-188.

LUO Zhigao,LI Xudong,ZHAO Junli,et al.Crack identification with mahalanobis distance discrimination method[J].Journal of Vibration and Shock,2013, 32(21):186-188.(in Chinese)

18 LI N, GRIFT T E, YUAN T,et al. Image processing for crop/weed discrimination in fields with high weed pressure[C]∥2016 ASABE Annual International Meeting, ASABE Papers 162460475,2016.

19 GANESAN P, RAJINI V, RAJKUMAR R I. Segmentation and edge detection of color images using CIELAB color space and edge detectors[C]∥2010 International Conference on Emerging Trends in Robotics and Communication Technologies, 2010: 393-397.

20 李飞，薛彬，黄亚楼.初始中心优化的K-means聚类算法[J].计算机科学, 2002, 29(7):94-96.

LI Fei,XUE Bin,HUANG Yalou.K-means clustering algotithm with refined center[J].Computer Science，2002,29(7)：94-96.(in Chinese)

21 CANNY J．A computational approach to edge detection[J]．IEEE Transactions on Pattern Analysis and Machine Intelligence，1986，8(6):679-698．

22 张富贵，洪添胜，肖磊，等．果树冠幅的检测机理研究[J]．农业工程学报，2008，24(4)：25-29．

ZHANG Fugui, HONG Tiansheng, XIAO Lei, et al. Detection mechanism for fruiter crown diameter[J]. Transactions of the CSAE, 2008, 24(4):25-29.(in Chinese)

23 WHEATON T A, TUMBO S D, WHITNEY J D, et al. Investigation of laser and ultrasonic ranging sensors for measurements of citrus canopy volume[J]. Applied Engineering in Agriculture, 2002, 18(3):367-372.

24 王佳, 杨慧乔, 冯仲科.基于三维激光扫描的树木三维绿量测定[J/OL]. 农业机械学报, 2013, 44(8):229-233.http:∥www.j-csam.org/jcsam/ch/reader/view_abstract.aspx?flag=1&file_no=20130839&journal_id=jcsam.DOI：10.6041/j. issn.1000-1298.2013.08.039.

WANG Jia，YANG Huiqiao，FENG Zhongke．Tridimensional green biomass measurement for trees using 3D laser scanning[J/OL]．Transactions of the Chinese Society for Agricultural Machinery，2013，44(8):229-233．(in Chinese)