基于改进k均值聚类方法的林木冠层孔隙度提取

沈德胜，朱良宽，宋佳音，李克新

(东北林业大学 机电工程学院，哈尔滨 150040)



基于改进k均值聚类方法的林木冠层孔隙度提取

沈德胜，朱良宽*，宋佳音，李克新

(东北林业大学 机电工程学院，哈尔滨 150040)

摘要：在利用林木冠层半球图像获取冠层孔隙度的过程中，针对因植被反光现象导致的图像误分割或分割不精确，以及传统k均值聚类图像分割方法易陷入局部最优值的问题，提出一种新的林木孔隙度提取策略。首先，利用数字图像处理方法消除了冠层半球图像中植被的反光，降低了因植被反光引起孔隙度计算出现较大误差的风险；然后，将初始聚类中心优化的k均值聚类方法与数字图像的数据特点相结合，确定了较为合理的数据密度搜索半径，提出改进k均值聚类方法对冠层图像进行两次聚类分割，计算出冠层孔隙度。最后，将结果与传统k均值方法、Otsu法和WinSCANOPY2006a软件的处理结果进行对比分析，实验结果表明所提的方法是可行且有效的。

关键词：图像反光；k均值聚类；冠层半球图像；孔隙度

0引言

林木冠层孔隙度是植被结构参数中重要的组成部分，也是获取其他冠层结构参数(如叶面积指数、植被郁闭度、叶面积等)的重要依据。冠层孔隙度是冠层图像中植被部分所占的比例，近年来，由鱼眼镜头获得地面植被半球图像计算植被冠层孔隙度的方法深受国内外学者的青睐，如何探寻一种计算准确、精度更高、计算速度更快的孔隙度获取方法具有重要意义[1-4]。

目前，应用数字图像处理技术处理冠层图像获得冠层孔隙度是一种高效实用的方法，并取得了一系列成果。文献[2]通过对植被半球图像的蓝色分量进行处理，可以有效抑制部分光线对冠层孔隙度提取的影响，但是当图像中植被部分有较强反光时孔隙度的计算结果并不理想。文献[3-5]通过对低矮植被图像HIS(色度，饱和度，亮度)颜色空间H分量图像进行图像分割虽然效果较好，但由于树木冠层结构复杂，H分量的植被区域与实际的植被区域差别较大，其对H分量图像的分割结果并不理想。

针对上述问题，本文提出一种基于改进k均值聚类的冠层孔隙度提取方法。首先，通过分析图像RGB(红，蓝，绿)和HSI颜色空间的色彩结构，并结合图像形态学对图像进行预处理，旨在消除因树木、枝叶等反光现象对图像分割的影响；然后，将初始聚类中心优化的思想与冠层数字图像的数据特点相结合，对传统k均值聚类方法进行了改进，利用改进后的k均值聚类方法对去除反光影响的图像进行两次聚类分割，并将两次聚类分割得到的植被区域合并；最后，将两次聚类分割得到的图像进行二值化处理，利用公式计算出冠层孔隙度，并将得到的结果与常规方法处理结果进行了对比验证。

1图像数据来源

实验用冠层半球图像来源于东北林业大学帽儿山实验林场阔叶林实验样地，采用带有鱼眼镜头的WinSCANOPY2006a设备，在样地内由下向上采集阔叶林冠层半球图像，试验用图像以WinSCANOPY2006a软件测定的冠层孔隙度作为真实值，用于验证不同方法计算得到孔隙度的准确性。

2冠层半球图像反光区预处理

2.1RGB与HSI颜色空间色彩结构分析

RGB和HSI空间的数字图像都是三维矩阵，通过相机获得的图像几乎都是RGB颜色空间的图像，RGB颜色空间数字化程度高，三个分量数据信息模式一致，使得数据计算方便快捷。但是，由于RGB三个颜色分量的相关度较高，任一分量的改变都会改变整个像素的颜色信息，在图像处理中容易相互影响，而H分量与I分量相关度很低，因此若将其转换成HSI颜色空间，则可将亮度分量与彩色图像中的彩色信息分开，从而实现只改变亮度信息而不改变颜色信息。

2.2反光区处理

反光区处理的基本流程如图1所示。

图1　反光区处理流程图Fig.1 Flow chart of reflective processing

RGB与HSI颜色空间转换公式如下：

色度分量为：

(1)

其中

(2)

饱和度分量为：

(3)

亮度分量为：

(4)

应用转换公式(1)～公式(4)，将原冠层图像由RGB空间转换为HSI空间，应用形态学的方法确定天空灰度为0的区域，将该区域灰度值修正为正常值，得到新的H分量图像H_NEW数据，再用I-λ×H_NEW得到差值图像(其中λ取值根据反光程度设定，反光越强λ越大，反之λ越小)。通过差值图像可以得到明显的反光区，再进一步确定反光区域的位置集合P，对I中相应区域进行减弱，I_NEW=m×I(P)(其中m为亮度调节参数)得到新的I分量数据I_NEW，将I_NEW与原来的H分量数据和S分量数据合成新的图像HSI_NEW，再将HSI_NEW图像转化为新的图像RGB_NEW。

图2　原图HSI空间H分量 Fig.2 The H component inHSI space for the original image

图3　形态学处理后H分量Fig.3 The H component aftermorphological processing

图2、图4为原图HSI空间中的H分量图和I分量图，图3为形态学处理修正后的H_NEW图，图5为I-λ×H_NEW得到的差值图，图中的亮区即树叶和树干的反光区，图6、图7为原RGB图及去除反光后的RGB_NEW图像。图8为去反光处理的对比结果，其中(a1)，(b1)，(c1)为原图，(a2)，(b2)，(c2)为去反光后结果。可以发现，应用本文方法可有效去除反光，并未改变原始图像非反光区的颜色信息和反光区色度信息，保留了图像的原始数据。

图4　原图HSI空间I分量Fig.4 The I component in HSI space for the original image

图5　差值图像Fig.5 Differencial image

图6　原RGB图像Fig.6 The original RGB image

图7　RGB_NEW图像Fig.7 RGB_NEW image

(a1)

(a2)

(b1)

(b2)

(c1)

(c2)

表1给出了图6中树干和树叶各3处反光区域、同类无反光区和去反光后的平均灰度对比结果，可以看出，本文去除反光的方法计算量较小，易于实现且效果较好，是一种林木冠层图像反光区预处理较实用的方法。

表1　反光区处理结果对比

本文中的图像数据是以鱼眼镜头在样地内由下向上采集阔叶林冠层半球图像，图像只有植被和天空两部分，文中图像分割的目的是获取植被和天空区域的面积即像素数量，用以计算孔隙度，但由于部分树叶和树干有较强的反光现象，直接导致实际的枝干、树叶占空比降低，因此需要对其进行上述(如流程图1)的预处理，以提高孔隙度的计算精度。

3改进k均值聚类法图像分割

3.1k均值聚类算法

设对有N个数据的集合A分为k类，也就是k个族，x表示Ci族中的成员，ci表示Ci族中心值[5-8]一般使用公式(5)：

(5)

为目标收敛函数，dis为集合A的族内成员对族内中心点的误差平方和。k均值聚类迭代过程可参见文献[5-8]，在此不赘述。

3.2改进k均值聚类方法

k均值聚类算法的不足之处在于：算法随机选取聚类中心，如果初始中心为孤立点，会导致迭代过程陷入局部最优值，得不到最优解；迭代的次数与初始聚类中心的选取密切相关，随机选取聚类中心也会使每次聚类时间不同，甚至最后的聚类结果也不同，增加了算法的不稳定性。针对冠层图像的数据特点和k均值聚类算法的不足，本文结合数据密度和最大化最小距离的方法，对k均值聚类方法进行了改进。

首先确立统计数据的区域半径，本文的数据区域半径 为全体数据中任意两个数据的距离之和的均值。计算公式如下：

(6)

以r为区域半径，统计任一区域内数据的数量，如果数量大于Min则表示以该数据x为中心r为半径的区域为高密度区，x为高密度点，统计所有的高密度点和相应的密度。

最大化最小值方法描述如下：[9-10]

(1)假设G={x1,x2,x3,…,xn}为已选好的高密度点数据集，找出其中密度最大的点J1。

(2)在集合G中找到距离J1最远的点为J2。

(3)统计G中其余对象与J1和J2的距离di1和di2，计算其中较小的值为min(di1，di2)。

(4)计算min(di1，di2)中的最大值，即最大的最小距离max(min(di1，di2))，以此xi为J3。

(5)以此类推J4=max(min(di1，di2，di3))最终确定k个初始聚类中心点。

需要注意的是：由高密度点和最大化最小距离方法选取的初始聚类中心点，避免了聚类结果的不确定性，避免了聚类局部最优的风险，可减少迭代次数，由此得到的初始聚类中心方法切实可行。

结合最大化最小距离方法的改进k均值聚类过程描述如下：

Step1：由公式(6)确定数据集A的区域半径r，对A中每个样本计算距离在r内的数据量，若大于Min则为高密度点，统计得到高密度点集G={x1,x2,x3,…,xn}。

Step2：在G中由最大化最小距离的方法确定出k个初始聚类中心。

Step3：由Step2中得到的k个初始聚类中心进行k均值聚类。

3.3冠层图像分割实验

本文采用两次聚类分割的方式：第一次从整个半球图像整体差异上区分天空与植被，第二次以第一次分割出来的含有细小枝叶的天空为整体，此时的待分割对象虽然很相近，但类别很少，易于聚类分割。将分割出来的细小枝叶与第一次分割出来的植被区域合并，剩下的则为天空区域。由于在可见光波段中，绿色植物对蓝光的吸收较强，在蓝色分量下对植物的绿色叶片与天空的区分非常明显，得到的分割结果更为准确，因此选取RGB_NEW中的蓝色分量(如图9所示)进行图像分割。

冠层半球图像像素大小为2211×2211，圆形区域的像素数量为3 821 941个。首先，结合冠层图像的特点，分别取第一次聚类分割的类别数k1=3、4、5、6、7，对去反光图像进行本文方法聚类。综合考虑分析处理结果和耗时长短，选择k1=6对图像进行聚类，结果如图10所示，按灰度值大小分为6类，同一类具有相同的灰度值，对分类后的图像化分为两类，黑色为植被区域，白色为天空区域，结果如图11所示。然后，考虑到还有细小枝干未被分割出来，对未分割完全的天空部分再进行第二次聚类分割，取第二次聚类分割的类别数为k2=2，结果如图12所示。可见利用本文的方法几乎可以将植被部分完全分割出来，取得了较好的冠层图像分割效果。

图9　RGB_NEW蓝色分量图Fig.9 Blue componentdiagram of RGB_NEW

图10　第一次k1=6聚类Fig.10 The first time k1=6clustering result

图11　第一次k1=6聚类效果图及细节Fig.11 The first time k1=6clustering result

图12　第二次k2=2聚类效果图及细节Fig.12 The second time k2=2clustering result

4孔隙度计算

孔隙度为天空区域所占的面积与植被区域和天空区域面积和的比值，计算公式如下[11-12]：

(7)

式中：T为孔隙度；Ns为天空像素总数；Nl为植被部分像素总数，由RGB_NEW的蓝色分量分割结果得出天空像素总数Ns=409 600，植被部分像素总数Nl=3 412 341，经计算得T=10.71%。

表2为针对冠层图像数据利用传统k均值聚类与本文算法选取不同k值时所花费的时间和聚类次数的对比，可以看出，本文方法可在一定程度上减少聚类次数，降低了聚类时间。

表2　传统k均值与本文方法聚类时间对比

目前，关于孔隙度的计算有直接测量法和间接测量法，直接测量法中有落叶收集法、斜点样方法等，测量准确但费时费力实现性差；间接测量法中有基于辐射量的方法和光学仪器法，测量仪器有HemiView、LAI-2000、CI-110、TRAC、DHP等，其中，WinSCANOPY2006a冠层图像分析仪作为一种通过间接测量法获得林木冠层参数的商业化仪器，测量精度较高，操作方便，在同类产品中测量效果较好，在业内受到广泛的认可，因此本文以WinSCANOPY2006a冠层图像分析仪的测量结果作为孔隙度真实值，用以验证不同方法孔隙度计算结果的有效性。

表3为针对采集的10幅冠层半球图像分别利用Otsu法、传统k均值方法及本文方法进行孔隙度提取，再将WinSCANOPY2006a软件的计算结果与不同方法的计算结果相对比。图13为根据10幅图像所得三种方法的孔隙度数据和WinSCANOPY2006a软件的处理结果，结果表明，相比传统k均值方法和Otsu方法，本文方法最接近WinSCANOPY2006a软件的处理结果，所得到的冠层孔隙度更准确。

图13　几种方法计算结果对比图Fig.13 The comparison charts of results with different methods

图 号孔隙度/%Otsu算法传统k均值算法本文方法WinSCANOPY2006a113.2812.4810.3210.67215.2811.8710.169.39317.3513.4511.4910.97414.1711.149.649.34511.8611.869.649.15613.9013.2610.8110.21717.4817.2414.7214.31810.8010.588.108.58912.1611.9310.149.95109.789.457.837.74

表4给出了三种方法与WinSCANOPY2006a处理结果的平均误差对比。由表4可以看出，本文方法平均误差小于5%，经统计分析得到与WinSCANOPY2006a的结果相关系数R=0.978 6，决定系数R2=0.926 1二者具有很高的相关性，说明了本文方法精度较高，结果较Otsu法和传统k均值法更为精确。

表4　三种方法与WinSCANOPY2006a

5结论与讨论

本文提出了一种新的林木冠层孔隙度提取方法，该方法考虑了植被反光现象对图像分割的影响，通过对图像进行预处理有效地去除了半球图像的反光，并利用改进k均值聚类方法对预处理后的图像进行两次聚类分割，计算出冠层孔隙度。通过实验可以看出，本文提出的方法去除冠层半球图像反光效果较好，而且在聚类分割的效率和精确度方面具有明显的优势，相比传统k均值方法和Otsu方法，所得到的冠层孔隙度更准确，为森林生长状况分析提供了有效的依据，为冠层孔隙度的获取提供了一种新的方法和途径。

【参考文献】

[1]潘桂颖.基于半球摄影法的水稻冠层结构参数测量与反演[D].杭州：浙江大学，2011.

[2]凌飞龙，李增元，陈尔学，等.青海云杉林叶面积指数半球摄影测量方法研究[J].地球科学进展，2009，24(7)：803-809.

[3]李明，张长利，房俊龙.基于图像处理技术的小麦叶面积指数的提取[J].农业工程学报，2010，26(1)：205-209.

[4]Han J，Kamber M.Data mining concepts and technique(2nd Ed)[M].Beijing：China Machine Press，2006.

[5]黄芬，于琪，姚霞，等.小麦冠层图像H分量的K均值聚类分割[J].计算机工程与应用，2014，(3)：129-134.

[6]张文君，顾行发，陈良富，等.基于均值-标准差的K均值初始聚类中心选取算法[J].遥感学报，2006，10(5)：715-721.

[7]MacQueen J B.Some methods for clustering and analysis of multivariate observations[M].Berkeley：University of California Press，1967.

[8]Ester M，Kriegel H P，Sander J，et al.A density based algorithm for discovering clusters in large spatial databases with noise[M].Portland：AAAI Press，1996.

[9]袁方，周志勇，宋鑫.初始聚类中心优化的k-means算法[J].计算机工程，2007，33(3)：65-66.

[10]张健沛，杨悦，杨静，等.基于最优划分的K-Means初始聚类中心选取算法[J].系统仿真学报，2009，(9)：2586-2590.

[11]杜春雨，范文义.有效叶面积指数与真实叶面积指数的模型转换[J].东北林业大学学报，2010，38(7)：126-128.

[12]李根柱，王贺新，朱教君.辽东山区长白落叶松叶面积指数和林冠开阔度的月动态[J].东北林业大学学报，2009，37(7)：20-22.

Extraction of Vegetation Hemisphere Image Porosity Based onImproved k-means Clustering Method

Shen Desheng，Zhu Liangkuan*，Song Jiayin，Li Kexin

(College of Mechanical and Electrical Engineering，Northeast Forestry University，Harbin 150040)

Abstract：Forest canopy hemisphere images can be used to capture the canopy porosity.However，the reflective phenomenon of vegetation can easily lead to image segmentation by mistake or inaccuracy，and the traditional k-means clustering image segmentation method is easy to fall into local optimal value，therefore a novel porosity extraction method was proposed in this paper.Firstly，the digital image processing technology was applied to eliminate the reflection of vegetation in the canopy hemisphere images，which reduced the risk of a larger porosity calculation error caused by vegetation reflection.Then，the method of k-means clustering algorithm was combined with the data characteristics of digital image to determine a reasonable searching radius of data density，and the image after being removed the reflective was segmented twice by the improved k-means clustering method.Finally，the canopy porosity was computed and the result was compared with the traditional k-means method，Otsu method，and WinSCANOPY2006a software processing results.The experimental results showed that the proposed method was feasible and effective.

Keywords：reflective image；k-means clustering；canopy hemisphere image；porosity

中图分类号：S718.42

文献标识码：A

文章编号：1001-005X(2016)02-0014-05

作者简介：第一沈德胜，硕士研究生。研究方向：图像处理技术林业应用。*通信作者：朱良宽，博士，副教授。研究方向：林业信息智能检测。Email：zhulk@126.com

基金项目：国家自然科学基金项目(31370710)；高等学校博士学科点专项科研基金资助(20110062110002)

收稿日期：2015-11-17

引文格式：沈德胜，朱良宽，宋佳音，等.基于改进k均值聚类方法的林木冠层孔隙度提取[J].森林工程，2016，32(2)：14-18.