基于改进K-means图像分割算法的细叶作物覆盖度提取

吴焕丽 崔可旺， 张 馨 薛绪掌 郑文刚 王 岩

(1.河北工业大学电子信息工程学院， 天津 300401； 2.北京农业信息技术研究中心， 北京 100097)

0 引言

植被覆盖度是反映作物个体/群体动态变化的重要参考量，是一种衡量作物生长状况的主要指标[1-2]，并与植物的光合有效吸收(APAR)相关，在一定程度上反映了植物的光合和蒸腾作用[3]。WANJURA等[4]认为遥感植被指数受覆盖度的影响大于干/鲜生物量或叶面积指数(LAI)等其他参数。因此，估测农作物的覆盖度在农业生产与科研上有非常重要的意义。

方雨晨等[5]研究了不同覆盖度下小麦农田土壤对归一化植被指数(NDVI)的影响。王春梅等[6]研究了叶面积指数和覆盖度的时空变异性。对于农作物覆盖度的估测，有以下几种方法[7-9]：①有经验的农学家依据肉眼进行估测，缺点是这种方法具有很强的主观性，不同的人估测，结果可能会有较大的差别。②在正午前后太阳角度较高时，利用直尺等仪器进行测量，缺点是易受天气因素影响，并且劳动强度过大。③采用相机或摄像头垂直拍摄图像进行逐点判断，优点是准确度很高，可实现在线或远程监测，缺点是信息采集受环境影响，依赖图像处理算法。近年来，作物表型研究在预估作物发育状况方面成为研究热点，具有动态监测作物性状的特点，其中无损图像分析处理已成为获取植物覆盖度重要手段。

获取植被覆盖度最关键是图像背景分割，准确地从背景中分割图像是减少表型分析误差的关键一步。有一些算法可以对植物叶片进行分类和分割，LIU等[10]提出了一种使用三次分水岭标记的自动叶提取算法。DEEPIKA等[11]提出了一种花卉植物分割方法，分别应用梯度特征的金字塔直方图和金字塔局部二元模式提取花卉植物的形状和纹理特征。国内许多学者也对背景分割做了相关研究。万园洁等[12]提出基于改进的全卷积网络图像语义分割方法，改善小麦图像分割中过分割以及欠分割现象。孙立新等[13]提出一种改进果蝇算法优化模糊均值聚类算法的图像分割算法。黄芬等[14]针对经典图像分割算法精度低、过分割等问题提出了基于HSI空间H分量的K-means算法，对小麦冠层图像进行分割。何建斌等[15]采用K-means算法和数学形态学相结合的方法进行分割，充分利用了小麦植株颜色和背景颜色的差异。上述研究分割对象多为大叶植物或作物种植密度较为稀疏的植物，多是从复杂背景中将作物目标分离出来，很少对拍摄对象进行连续试验。

为克服细叶作物分割精度低、图像分割结果未得到应用的问题，本文提出一种融合小波分析按比例萎缩去噪算法、改进果蝇算法以及K-means算法的图像分割算法，以获取小麦覆盖度，通过多幅图像的试验对其性能进行检验，并通过连续试验探索覆盖度与作物系数之间的关系。

1 图像数据采集

数据采集地点为北京市昌平区小汤山国家精准农业研究基地，选取24套小型蒸渗仪上种植的小麦进行图像采集，如图1所示。试验时间为2017年9月30日—2018年6月5日。采集点共24个小区，编号为1～24，每4个小区为一组循环，每个小区面积为0.7 m×1 m，分别按照参比蒸散的30%、60%、90%以及大水漫灌的方式进行灌溉，选择大水漫灌区域(图1中红色区域)的小麦进行图像数据采集和处理，编号分别为4、8、16、20、24。

图1 数据采集地点示意图Fig.1 Sketch of data collection site1.农大212品种小麦 2.中麦1062品种小麦 3.小麦保护带 4.地下小型蒸渗仪 5.蒸渗仪地下监测室

在2018年3—5月期间，共拍摄10次，使用Canon EOS 7D相机(5 184像素×3 456像素)和DM10型红外热像仪(160像素×120像素)，垂直于地面放置，相机距离地面1.3 m、DM10型红外热像仪距离地面1.5 m处进行拍摄。为降低相机几何畸变给图像带来的影响，在拍摄时尽量避免使用镜头焦距的最广角端或最远摄端，并使用较小的光圈。每隔7 d对5块小麦地分别在09:30、11:00、13:00、14:30、16:00 5个时刻拍摄5组数据，每次重复拍摄3幅图像以消除系统误差，每块区域共采集150幅高像素图像以及150幅低像素图像，均以JPEG格式存储到计算机中。因天气条件的变化，大部分图像是在晴天拍摄，也有部分是阴天与有风天气下拍摄。在Matlab R2017b试验平台实现本文图像分割算法对数据的处理。

2 图像预处理

图像预处理是整个图像处理的关键部分，对图像分割和图像分析都有着至关重要的作用。由于待处理图像的采集环境是自然环境，受光照影响较大，针对这种状况，本文图像预处理的重点为图像光照分布。

选择恰当的彩色空间是进行有效分割的基础[16-17]。而面向设备的RGB颜色空间由于其分量间高度线性相关，使其不适用于直接分割彩色图像。由于在HSV色彩空间中的各分量值固定在[0,1]区间内，有利于目标对象特征参数的提取[18]，可以避免果蝇算法随机搜索到范围之外。为便于果蝇算法的计算，本文将图像由RGB色彩空间转换到HSV色彩空间，实现分割空间的选择，通过Matlab函数库中的rgb2hsv以及hsv2rgb函数实现图像在RGB与HSV色彩空间之间的转换。

3 传统K-means算法

K-means算法的基本思想是以空间中k个点为中心进行聚类，对最靠近他们的对象进行归类。通过迭代法逐次更新各聚类中心的值，直至得到最优的聚类结果[19]。

在K-means算法中，可以人为输入参数k，然后将输入的n个数据对象划分为k个聚类，以便使所获得聚类满足同一聚类中的对象相似度较高；而不同聚类中的对象相似度较低[20-22]。假设把样本集分为c个类别，算法描述如下：

(1)适当选择c个类的初始中心。

(2)在第k次迭代中，对任意一个样本求其到c个聚类中心的距离，将该样本归到距离最短的中心所在的类。

(3)利用均值法更新该类的中心值。

(4)对于所有的c个聚类中心，重复步骤(2)、(3)，如果聚类中心值保持不变，则迭代结束，否则继续迭代。

此算法存在的主要问题是易陷入局部最优，因此对其进行改进。

4 算法改进

使用K-means算法作为基础算法，提出两点改进：①去噪和改进的自适应步长果蝇算法解决局部最优问题，得到最优化的初始聚类中心。②用Matlab函数库中的kmeans函数完成图像的K-means算法分割，并将之与本文算法进行比较。自适应果蝇K均值聚类(IFOA-K-means)算法流程图如图2所示。

图2 算法流程图Fig.2 Flow chart of algorithm

IFOA-K-means算法步骤如下：

(1)首先对图像进行预处理。为降低异常点对后续K-means算法的影响，选择小波比例萎缩去噪方法进行图像去噪处理，并将图像色彩空间由RGB空间转换为HSV空间。

(2)初始化聚类中心选择。利用自适应步长果蝇算法找到最优化的初始聚类中心。

(3)分割图像。利用K-means算法对小麦HSV彩色图像进行聚类分割，当准则函数收敛时，迭代中止，调用Matlab函数库中的kmeans函数实现。

(4)边缘处理。利用形态学的开闭操作对分割图像的灰度图像进行校正，使图像分割效果更好。

图3 分割算法流程图Fig.3 Flow chart of split algorithm

(5)将本文小麦图像分割算法与基于遗传算法的最大类间方差法[23-24]以及传统K-means算法做对比。分割算法流程图如图3所示。

4.1 小波比例萎缩去噪算法

由噪声小波系数[25]特性可知，如果能够根据小波系数的局部特征及时地调整去噪处理方法，使其具有很强的自适应性，去噪效果将会得到很大改进。而局部调整的极限就是对每个小波系数都采用不同的处理方法，小波比例萎缩法的自适应性就很好地做到了这一点。本文所处理的图像环境是自然环境[26]，受光照影响较大，因此本文提出只针对V分量做去噪处理，具体步骤如下：

(1)对加噪后图像y(i，j)(设图像尺寸为m×n)做正交小波变换得到小波系数Y(i，j)，1≤i≤m，1≤j≤n。

(2)选用正方形窗口，把所估计的小波系数Y(i，j)放在窗的中央。窗口尺寸可以选择3×3、5×5、7×7等。2(i，j)的计算公式为

(1)

式中Ω(i，j)——包含Y(i，j)的窗口内小波系数，窗口尺寸为M×M

在对边缘系数估计时，还需要根据窗口尺寸对边缘进行延拓，常用的方法是对称延拓。

(2)

其中

median(·)——中值函数

(3)

图4 去噪算法效果图Fig.4 Effect diagrams of denoising algorithm

对采集到的图像进行V分量去噪和H、S、V3个分量同时去噪，由处理效果图可以看出，单独对V分量去噪比同时对H、S、V3个分量去噪平滑程度大，且没有出现新的噪点，而同时对3个分量去噪导致图像上出现了不规则的黑色噪点，主要集中在红色记号笔圈出区域，小波按比例萎缩去噪算法模型无法同时使3个分量的去噪效果都达到最佳，因而在去噪后仍有部分噪点，影响之后的图像背景分割效果。

4.2 改进自适应步长果蝇算法

果蝇算法是一种基于果蝇觅食行为推演出的寻求全局优化的新方法[27-29]。在基本果蝇算法中，迭代寻优采用固定的步长，不利于算法的收敛和稳定，如果步长取值过小，会降低算法的收敛速率，也会使得收敛的精度降低，但步长过长又会导致果蝇算法跳过最优解，稳定性降低，极易出现振荡等问题，因而本文对果蝇算法的步长进行了改进。

为提高果蝇算法的收敛速率、精度以及稳定性能，本文提出一种根据聚类误差自适应调整步长的果蝇算法，根据聚类误差的大小来调节下一次迭代的步长，调整原则如下：①如果当前迭代得到的聚类中心在图像上的聚类误差(果蝇的味道浓度判定值)越小，则步长应该越小，所以这里使用当前的最小味道浓度(loss)作为下一次迭代的参考量。②如果当前迭代的loss大于上次迭代的loss，则表明上次迭代得到的聚类中心是一个比较优的值(可能是局部最优，也可能是全局最优),单纯从这一点来看，后面的迭代应该减小步长。

原则①可以保证在loss较大时快速进行大范围搜索，便于找到理想的极小值点，因此有助于跳出不理想的局部极值点。原则②可以保证在理想的局部极值点处有更精确的收敛(即使用更小的步长)。

进入迭代寻优，迭代步长根据两个自适应步长迭代原则变换，如果最佳味道浓度低于前一迭代最佳味道浓度[30]，或者当前迭代次数大于最大迭代次数，则终止迭代。本文自适应果蝇K均值算法与基本果蝇算法运行对比，如图5所示。

图5 果蝇算法对比Fig.5 Comparison chart of fruit fly algorithm

由图5可知，本文自适应果蝇K均值算法比基本果蝇算法的收敛速率更大，浓度判定值更小，即聚类误差更小。

5 试验结果与分析

5.1 试验效果

5.1.1分割效果

为验证本文算法的性能，将本文算法与传统的K-means算法以及基于遗传算法的最大类间方差法作对比。

(1)不同时期分割效果

为方便看出分割效果的区别，使用红色记号笔对上述图像相同位置做了标记，图6～8分别为3月29日09:30、4月11日09:30、5月3日09:30拍摄的4区小麦图像。由图6d、7d可以看出，在小麦返青期和拔节期利用最大类间方差法虽然分出了小麦的大致形态，但由于小麦初期有小部分泛黄区域，因而在选择阈值进行分割时，将不属于作物的部分也划分为前景。该算法易受强烈光照的影响，如图8d所示右下角部分小麦未被识别出。在小麦返青期和拔节期利用K-means算法可以在一定程度上分割出小麦图像，但是在小麦较为密集的地方出现模糊现象，而且在小麦抽穗期，分割失真情况也比较严重；而利用本文算法分割小麦图像，在麦底和麦尖处的分割都较为清晰，边缘也得到了较好的处理，避免了K-means算法处理图像时出现的模糊现象。

图6 返青期3种分割算法效果Fig.6 Effect diagrams of three segmentation algorithms for return period of wheat

图7 拔节期3种分割算法效果Fig.7 Effect diagrams of three segmentation algorithms for jointing stage of wheat

图8 抽穗期3种分割算法效果Fig.8 Effect diagrams of three segmentation algorithms for heading stage of wheat

(2)不同天气分割效果

选取晴天(712.6 W/m2)、阴天(18.0 W/m2)、有风天气(7.4 m/s)下小麦图像，图9为4月16日13:00拍摄的4区小麦图像，图10为4月17日16:00拍摄的8区小麦图像，图11为4月10日09:30拍摄的16区小麦图像，分别使用K-means、IFOA-K-means、最大类间方差法3种算法对图像做背景分割。使用最大类间方差法处理时，阴天效果最佳，晴天出现强烈光照或有较明显阴影区域对分割效果有很大的影响，如图9d所示。使用K-means算法处理时，均可在不同天气下将小麦背景去除，但是在很多叶子的细节部分，出现黑色孔洞，而IFOA-K-means算法在不同天气下对细节部分处理的效果均较好，如图9c、10c、11c红色记号笔圈出区域。

图9 晴天3种分割算法效果Fig.9 Effect diagrams of three segmentation algorithms under sunny day

图10 阴天3种分割算法效果Fig.10 Effect diagrams of three segmentation algorithms under cloudy day

图11 有风天气3种分割算法效果Fig.11 Effect diagrams of three segmentation algorithms under windy day

通过上述不同时期、不同天气下的分割效果图可以看出，最大类间方差法在处理图像时易受到天气的影响，并且在小麦生长初期阶段，由于小麦、土壤、灌水带之间的分割阈值不易确定，很容易将不是小麦的部分划分为前景。而K-means算法相对来说处理效果比较稳定，只是在有些细节部分的处理效果图有孔洞或模糊现象；本文算法在处理不同时期、不同天气的图像时都有较好的稳定性以及分割效果，小麦叶子细节部分也较为清晰。

5.1.2不同尺寸图像峰值信噪比数据对比

由于最大类间方差法易受小麦生育期、天气等因素的影响，因而后续分析主要针对本文算法以及K-means算法。峰值信噪比(PSNR)[31]是评价图像的客观指标，值越大，说明分割图像越逼真。本文计算了拔节期两种不同尺寸图像用不同算法分割的运行时间、信噪比，结果如表1所示。

表1 不同算法性能评估结果Tab.1 Performance evaluation result of different algorithms

由表1可以看出，图像尺寸会直接影响算法运行的结果。对于同一算法，低像素图像由于总的像素点少，在运行时间上比高像素图像缩短了70%以上，但是在峰值信噪比上，比高像素图像降低了40%以上。本文算法在处理高像素图像时，运行时间比K-means算法缩短了53%，峰值信噪比提高了21%，在对低像素图像进行处理时，由于像素较低导致了迭代速率减小，所以本文算法运行时间略高于传统算法，但是峰值信噪比提高了20%。

5.2 小麦覆盖度

植被覆盖度指包括乔木、灌木、草和农作物在内的所有植被的冠层、枝叶在生长区域地面的垂直投影面积占研究统计区域面积的百分比[32]。图像采集时为降低图像几何畸变应尽量避免使用镜头焦距的最广角端或最远摄端，并使用较小的光圈。经图像处理后将小麦的图像背景去除，利用分割后的灰度图像计算小麦像素点数占总像素点数的比例，即小麦垂直投影面积占总面积的比例，所求结果便是研究区域的小麦覆盖度。

5.2.1不同方法提取的小麦覆盖度

选取小麦具有代表性的3个生长时期图像数据，通过采用AutoCAD软件对图像中的绿叶进行实际勾绘计算覆盖度[33-34]，将其作为小麦覆盖度真值与IFOA-K-means、K-means算法处理结果做对比，选择农大212品种(16区)和中麦1062品种(4区)各1个小区，每个小区每个生长时期选取正午图像各2幅，对比结果如表2所示。

表2 不同算法处理的小麦覆盖度对比Tab.2 Comparison of wheat coverage with different algorithms %

注：同行数据后不同字母表示处理间差异显著。

由表2可知，IFOA-K-means算法与覆盖度真值无显著性差异，K-means算法与覆盖度真值存在显著性差异，本文算法结果更接近于覆盖度真值，准确率在90%以上。

5.2.2基于IFOA-K-means算法的覆盖度结果

对4区和16区小麦覆盖度结果进一步分析，根据2个小区日覆盖度均值做出覆盖度变化趋势图，如图12所示，详细数据如表3所示。2个小麦品种的增长趋势大致相同，拔节期生长迅速，并在5月中旬前后，覆盖度达到最大值后开始下降。返青期初期，由于4区小麦较16区小麦距遮雨棚更近，生长时光照强度低于16区小麦，导致初期中麦1062品种小麦(4区)覆盖度低于农大212品种小麦(16区)；而后期，4区小麦距水泵更近，可能水量较16区稍多，因而在灌水后4区小麦生长速率超过16区。可见，光照和水分都会对小麦覆盖度的变化产生影响。

图12 日覆盖度变化趋势Fig.12 Changing trend of canopy coverage

区域3月22日3月29日4月3日4月11日4月16日平均值标准差平均值标准差平均值标准差平均值标准差平均值标准差4区25.652.0030.601.7748.890.5663.471.5669.421.4216区31.881.7934.221.7240.812.0262.860.37区域4月25日5月3日5月10日5月18日5月23日平均值标准差平均值标准差平均值标准差平均值标准差平均值标准差4区82.470.9585.741.2388.750.5481.800.8776.350.8516区83.602.1485.932.3687.721.8285.321.1277.214.64

5.2.3日覆盖度、蒸渗仪、气象数据联合分析

根据图像获取的日覆盖度数据、蒸渗仪获取的实际蒸散量数据、气象数据计算得到的参比蒸散量数据可得到日覆盖度与作物系数之间的关系，即由图像数据便可推测作物系数。

根据蒸渗仪数据以及气象数据得到作物系数为

Kc=ETc/ET0

(4)

式中ETc——由蒸渗仪土体电压数据计算得到的作物实际蒸散量

ET0——由气象数据计算得到的参比蒸散量

Kc——作物系数

由3—5月获取的5个小区图像数据计算得到的日覆盖度与对应时间计算得到的作物系数进行拟合，得到两者之间的关系

f(x)=12.97x3-18.17x2+8.537x-0.818 1

(5)

式中x——日覆盖度f(x)——作物系数

决定系数R2为0.953 1，拟合曲线如图13所示。

图13 日覆盖度与作物系数关系曲线Fig.13 Relationship curve of canopy coverage and crop coefficient

根据式(5)便可通过图像数据推测作物系数。日覆盖度在20%～50%时，曲线平缓，作物系数随覆盖度的增长缓慢上升，在50%～90%之间时，作物系数随覆盖度的增长而快速上升，即拔节期之后，小麦生长比较迅速，也是灌水的最佳时期，作物管理者也需在这个时期多加注意小麦覆盖度的变化，及时做出相应处理。

6 结论

(1)从不同天气、不同时期的去除效果、算法运行速率以及覆盖度的计算等方面，对最大类间方差法、K-means算法以及本文提出的IFOA-K-means算法作对比，得到本文算法均优于前2个算法，计算覆盖度准确率在90%以上。

(2)光照和水分都会对小麦覆盖度产生不同程度的影响，在返青期，光照充足会使小麦迅速生长；在拔节后期，水分对小麦覆盖度影响较大，如果缺水，将严重影响小麦的品质和产量。因而，在小麦生长过程中，要兼顾这两个因素。

(3)小麦覆盖度和作物系数具有相似的变化规律，前期由小变大，成熟期时达到最大值后逐渐减小。利用Matlab软件平台建立由覆盖度预测作物系数的模型，拟合方程的决定系数R2为0.953 1。