基于高光谱图像波段融合的猕猴桃软腐病早期分类检测

高宏盛,郭志强*,曾云流,丁 港,王逍遥,李 黎

1. 武汉理工大学信息工程学院,湖北 武汉 430070 2. 华中农业大学园艺植物生物学教育部重点实验室,国家柑橘保鲜技术研发专业中心,湖北 武汉 430070 3. 中国科学院植物种质创新与特色农业重点实验室,中国科学院猕猴桃产业技术工程实验室,中国科学院种子创新研究院,中国科学院武汉植物园,湖北 武汉 430074

引 言

猕猴桃富含丰富的维生素和矿质元素,被誉为“水果之王”,深受消费者的喜爱。 然而由于猕猴桃果实质地柔软,皮薄多汁,极易受到病原微生物的侵染而致病。 其中,软腐病是目前猕猴桃采后最为严重的真菌病害,感病果实在贮藏前期完全无症状,在贮藏后期及销售阶段逐渐发病,发病部位表皮凹陷变软,果皮上形成圆形或椭圆形的病斑,病斑边缘呈现水渍状,严重时会导致整个果实腐烂,并具传染性导致整库腐烂[1],造成极大的经济损失。 由于软腐病发病前期变化不明显,导致人工筛选效率低下。 因此,开发一种快速软腐病的早期无损检测的方法成为果农、 销售者、 消费者等共同关注的一个问题。

近年来,基于高光谱成像技术的无损检测研究已经成为热点。 高光谱成像技术融合了光谱技术和机器视觉技术。 机器视觉技术能够通过红、 绿、 蓝三个波段观测病害部位的大小、 纹理等来确定病害类型,但它无法探知病害的生化参数,而光谱技术可以获取样本的生物和化学参数,但它不能反映样本的全局信息。 使用高光谱成像可以一次获得光谱信息和空间信息,快速完成无损检测工作[2],已广泛应用于水果[3-5]、 蔬菜[6]、 肉类[7]等农副产品的相关无损检测。 目前,桃[8]、 南果梨[9]、 草莓[10]、 库尔勒梨[11]、 蓝莓[12]、 柑橘[13]、 苹果[14]等常见水果的高光谱成像病虫害检测均能达到快速、 无损的早期分类检测效果。 这些研究表明,高光谱成像技术在针对水果病虫害的早期无损检测具有很大的应用潜力。

对于猕猴桃样本,高光谱成像技术目前在可溶性固形物含量检测[15]、 冷害分级[16]、 酸度检测[17]、 成熟度检测[18]等方面均有相关报道。 这些研究表明,高光谱成像技术能够反映猕猴桃的内在成分,比如内在品质参数、 成熟度、 损伤等方面,而利用高光谱成像技术对猕猴桃软腐病进行早期检测鲜见报道。 结合前人相关研究,探索和验证使用高光谱成像技术对猕猴桃软腐病进行早期检测的可行性。

上述使用高光谱成像对水果进行无损检测的大部分文献仅使用了光谱特征作为分类依据,而没有利用高光谱图像“空谱合一”的优点。 对于高光谱图像特征的提取,主要是对其进行纹理特征的提取,Sun等[8]基于主成分分析(principal component analysis,PCA)的PC1图像进行了纹理特征的提取,Yu等[9]基于PC1和PC2图像进行纹理特征的提取,两者研究结果均表明通过结合高光谱图像的光谱特征与图像特征能够进一步提高分类的准确率。 然而,使用主成分图像进行图像特征提取时,容易存在较大的光谱失真,不能充分利用高光谱图像多波段图像的信息。 本工作以“云海一号”猕猴桃为研究对象,首先将猕猴桃的特征波段图像融合得到融合图像,并提取融合图像的纹理特征。 同时,提取高光谱图像的光谱特征,然后将光谱特征与纹理特征融合后使用支持向量机(support vector machine,SVM)进行分类,对猕猴桃感染软腐病后进行了早期检测与识别。

1 实验部分

1.1 材料

实验样本采摘自湖北省武汉市初阳奇异果园,选择成熟度处于商业采收期、 可溶性固形物为8%～9%、 果重在80～110 g的“云海一号”猕猴桃作为样本,并剔除了畸形果和疤痕果,然后使用咪鲜胺和抑酶唑进行处理,等待晾干后放入冷库中保存。 用于样品软腐病接种的是拟茎点霉菌(Phomopsissp.),由华中农业大学园艺林学学院相关研究专业人员对猕猴桃进行接种。

由于猕猴桃的球形形状会导致高光谱图像中光谱反射的不均匀性,为了最大限度的减少这种影响,接种区域选择在猕猴桃的中央赤道部位。 首先在无菌条件下使用灭菌后的移液枪进行接种。 实验分为健康对照组和实验组,健康对照组注射无菌生理盐水来消除接种处伤口对分类识别的影响,实验组使用两种不同剂量的接种方法来模拟不同情况下的猕猴桃软腐病发病情况,其中第一组接种6 μL孢子悬浮液,第二组接种10 μL孢子悬浮液。 接种完成后,对所有猕猴桃进行单果套袋并将接种面朝上,室温20 ℃贮藏。 将接种当天记为第0天,并从第1天开始每天下午1点进行高光谱图像的采集,同时关注软腐病的发展状况,直到猕猴桃软腐病发病严重或出现白色菌丝后停止高光谱图像的采集。

1.2 高光谱图像采集系统

采用比利时IMEC公司开发的SNAPSCAN VNIR高光谱成像仪采集高光谱图像,采集系统如图1所示,该系统主要由相机机身、 镜头、 升降平台、 卤素灯光源及扫描控制器组成。 该系统配有官方的图像采集软件HSISnapscan,HSI(hyperspectral image); 传感器有150多个可用的波段,经系统处理后,在470～900 nm之间有150个可用光谱波段,利用软件可自动黑白校正得到校正后的光谱反射率数据。 根据猕猴桃大小选择的成像分辨率为800像素×800像素,相机积分时间为4 ms,模拟增益为1.6 dB。 图像采集时,整个系统放置于一个密闭的黑箱中,防止外部光源的干扰。

图1 高光谱图像采集系统(a): 示意图; (b): 实物图Fig.1 Hyperspectral image acquisition system(a): Schematic diagram; (b): Actual picture

1.3 方案设计

实验方案由三部分组成: 图像采集及处理、 特征提取及融合、 模型建立及评估。 总体实验流程如图2所示。

图2 总体实验流程Fig.2 Overall experimental process

所使用软件除1.2节中提到的HSI Snapscan外,图像感兴趣区域的选择使用了Labelme软件,除此之外其他数据处理及建模分析均由Python3实现。

1.4 高光谱图像处理及数据处理

1.4.1 猕猴桃图像信息及光谱特征分析

图3是不同类型的猕猴桃样品,其中图3(a)代表健康猕猴桃,图3(b)、 (c)和(d)分别代表接种拟茎点霉菌(6 μL)后第四天、 第五天和第七天的猕猴桃。 由图3可以看出,接种拟茎点霉菌后第四天的猕猴桃与健康猕猴桃并无明显区别,而在第五天,猕猴桃接种部位变软,说明从第五天起肉眼观察可以识别出软腐病发病。 在接种后的第七天猕猴桃软腐病发展严重,有明显的水渍状。

图3 猕猴桃样品(a): 健康; (b): 接种第四天; (c): 接种第五天; (d): 接种第七天Fig.3 Kiwifruit samples(a): Healthy; (b): The fourth day of vaccination; (c): The fifth day of vaccination; (d): The seventh day of vaccination

根据猕猴桃软腐病的发展过程,将所采集的猕猴桃样本分为三类: 第一类为健康对照组,第二类为软腐病早期,第三类为软腐病晚期。 软腐病发病早期是根据病斑肉眼看不出明显扩大,其大小为接种口的直径为0～10 mm,晚期是指病斑肉眼可见的扩大,其直径>10 mm。 实验过程中发现,第二实验组比第一组发病提早一天,因此,用于后续分类研究的样本为: 第四天至第七天样本(实验组一),第三天至第七天样本(实验组二)。 最终用于实验的猕猴桃样本有295个,其中第一类样本60个,第二类样本128个,第三类样本107个。

将所有样本分类标定后,每个样本以接种部位为中心,选择对应的感兴趣区域(region of interest,ROI)用于后续的图像特征提取,其中ROI区域的平均反射率光谱代表样本的原始光谱特征。

图4显示了不同类别猕猴桃的平均反射率光谱曲线,可以看出不同病害程度的猕猴桃平均光谱曲线总体趋势上是相似的,在630 nm附近与670 nm均存在波峰和波谷,这主要是由于猕猴桃果肉和表皮中的叶绿素吸收引起的。 软腐病早期和晚期的猕猴桃光谱曲线比较接近,但健康猕猴桃的光谱反射率比软腐病猕猴桃低,这可能是真菌繁殖过程导致猕猴桃果实中的成分发生了变化。 这说明光谱曲线具备将健康果与软腐病果初步区别的潜力,但是无法将发病早期和晚期果准确区分,需对光谱特征处理进一步分析。

图4 不同类别猕猴桃光谱信息Fig.4 Spectra of different categories of kiwifruit

1.4.2 光谱特征提取

为了减少光谱数据中噪声的影响,通常需要对原始光谱数据进行预处理,常用的光谱预处理方法有: SG卷积平滑(savitzky-golay)、 多元散射校正(multiplicative scatter correction,MSC)和变量标准化(standard normalize variate,SNV)。 MSC算法通过拟合样本与理想光谱的关系来对光谱进行校正,而SNV算法则是通过对每条光谱做标准正态变换来对光谱进行校正。 两者均是在已知样本所属类别的基础上对该类别的平均光谱曲线进行校正,然而在之后的模型测试阶段,样本类别是未知的,故不采用这两种预处理方法。 SG卷积平滑预处理可以减少随机噪声对光谱的影响,而本研究所采用的高光谱仪器在导出光谱数据时已进行过校正,有效减少了随机噪声的影响。 对于输入特征只包括光谱特征时,分类器采用SVM,表1展示了使用SG卷积平滑与不使用预处理方法的实验结果对比,可以看出,使用SG卷积平滑后,最好的准确率与不使用预处理相同,但在SPA或CARS光谱特征提取下,使用SG得到的结果要更差一点,这说明,对于本文实验中采集的高光谱数据而言,已不需要再进行数据的平滑,使用SG平滑反而会造成过度平滑恶化结果。 综合考虑,在获取原始光谱特征后,不进行光谱预处理直接进行光谱特征提取。

表1 光谱数据未处理和经预处理的结果对比Table 1 Results of input data with different preprocessing methods

采集的高光谱图像包含150个波段,使用全波段进行建模时,分析时间长,且各波段间相关性较强,很容易造成信息冗余,因此需要对原始光谱进行特征降维或者特征选择。 采用PCA进行特征降维,选择连续投影算法(successive projections algorithm,SPA)和竞争性自适应重加权算法(competitive adaptive reweighted sampling,CARS)进行特征波段的选取。

PCA降维通过对原始波段特征进行线性加权组合形成新的互不相关的主成分变量,根据主成分贡献率的大小对主成分数量进行选择,可以有效降低特征维度。 对原始光谱特征进行PCA降维后,PC1、 PC2、 PC3贡献率分别为81.98%、 14.7%、 2.1%。 前三个主成分在三维空间中的分布情况如图5所示,从图中可以看出,三种类别的样本之间有所交叉,说明仅根据前三个主成分还无法完全对样本进行分类,后续还需借助分类器进行建模分析。 为了保证后续分类器的分类性能,选择保留累计贡献率达到99.9%的前10个主成分作为原始光谱数据降维后的光谱进行下一步分析。

图5 前三个主成分分布Fig.5 The first three principal components distributions

SPA算法根据最大投影向量不断循环,在筛选特征波段的同时降低了变量组之间的共线性关系。 设置波段选择过程中的最少波段数为2,最多波段数为50,发现在选择波段数为3时,均方根误差波动幅度减小,均方根误差(RMSE)为0.49。

CARS算法通过将蒙特卡洛采样与PLS回归系数相结合,在每次自适应加权采样过程中,保留回归系数绝对值大的波段,剔除绝对值较小的波段,经过多次采样选择交叉验证均方根误差(RMSECV)最小的波段子集。 设置蒙特卡洛采样次数为50,5组交叉验证,发现在采样次数为15时,RMSECV最低,此时筛选的特征波段数时39个。 表2展示了SPA和CARS算法所选择的具体波段,SPA选取的波长处于可见光区,而CARS选取的波长在可见光和近红外区均有分布。

表2 特征波段选择结果Table 2 Feature band selection results

1.4.3 图像特征提取

高光谱图像各波段的图像可以清晰的反映软腐病猕猴桃的图像特征,因此,在提取光谱特征的基础上加入波段图像特征,可以进一步强化猕猴桃的病害信息,增强健康组与不同程度病害组猕猴桃之间的区分度。 由于采集的高光谱图像波段数多达150个,如果对每个波段进行图像特征提取会造成极大的计算量,并且其中存在一些对分类贡献较小的波段,因此,在提取图像特征之前,需要进行特征波段的选择。

选择的特征波段是在SPA过程中产生的一个波段范围子集,共有8个波段,分别是第0、 11、 27、 62、 77、 87、 132、 149波段,对应波长为470、 502、 548、 649、 692、 721、 851和900 nm。 没有使用SPA最终选择的3个波段的原因是: 对于光谱特征而言,这3个波段的平均反射率在光谱信息上具有很好的代表性,但对于图像特征而言,不同波长图像下样本对光进行吸收反射后对图像特征的反映也不同,所筛选的8个波段能够很好的覆盖原始150个波段中不同的波段范围。

采用灰度共生矩阵(gray-level co-occurrence matrix,GLCM)进行图像纹理特征的描述。 设置不同的距离(2、 4、 8、 16)以及不同的角度(0°、 45°、 90°、 135°)进行灰度共生矩阵的提取,并计算灰度共生矩阵的对比度(contrast)、 同质性(homogeneity)、 相关性(correlation)、 能量(ASM)、 熵(entropy)。 将4个距离以及4个角度的平均值作为图像的纹理信息,将每幅图片计算得到的5个纹理参数代表该图像的纹理特征。

对于各子带系数的融合规则,通常对低频子带进行加权平均,对高频子带进行绝对值取最大的方法来求解。 经过实验发现,在使用绝对值最大法对高频子带进行融合时,由于每个像素值只用到了一个波段图像信息,最后所得到的高频子带系数将存在一定的信息丢失。 因此,为了在简化融合规则的同时考虑融合更多的细节信息,对于高频子带系数也采取了加权平均的方法。 式(1)和式(2)分别描述了低频子带和高频子带的融合规则。

(1)

(2)

图6 高频子带图像(a): 第132波段; (b): 绝对值最大融合; (c): 加权平均融合Fig.6 High frequency subband image(a): Band 132; (b): Absolute maximum fusion; (c): Weighted average fusion

1.4.4 模型建立与评估

采用Kennard-Stone算法将采集的实验样本按照训练集∶测试集=7∶3进行划分,共获得训练集207个,测试集88个。 由于样本的光谱特征与图像特征量纲不同,在将两种特征融合之前需要进行归一化处理。 采用最大最小值归一化方法将两组特征归一化至[0,1]之间,计算过程如式(3)

(3)

式(3)中,X表示原始特征,Xmin、Xmax分别代表原始特征中的最小值和最大值,Xnorm为归一化至[0,1]之间的特征,作为后续分类器的输入。

选用三种分类模型作为验证,分别是最近邻算法(K-nearest neighbor,KNN)、 随机森林(random forest,RF)和支持向量机(SVM)。 模型在训练过程中采用了十折交叉验证,并通过测试集的分类准确率(accuracy)和混淆矩阵对模型进行评价,分类准确率的计算过程如式(4)

(4)

式(4)中,TP、TN分别表示模型预测正确的正、 负样本,FP表示模型预测为正类的负样本,FN表示模型预测为负类的正样本。

2 结果与讨论

2.1 识别准确率分析

基于不同处理组合对猕猴桃样本数据进行分析,实验结果如表3所示,其中前六种方式利用的是猕猴桃样本的高光谱图像数据或RGB图像数据,分为仅使用光谱特征、 仅使用纹理特征和使用光谱、 纹理融合特征。 最后一种是将RGB光谱特征结合HSI波段融合图像的纹理特征进行分析,即融合了RGB光谱数据和高光谱图像数据。

表3 不同组合下测试集分类结果Table 3 Classification results of test sets under different combinations

从表3可以看出,(1)KNN是三种分类模型中分类效果最差的,RF与SVM两种模型的测试集分类准确率较为接近,但使用RF分类器时,训练集和验证集的分类准确率大部分为100%,而测试集准确率基本处于80%～88%之间,存在较严重的过拟合现象。 使用SVM分类器时,训练集、 验证集和测试集的分类准确率基本相近,有不错的泛化能力,因此,SVM是最适合识别猕猴桃软腐病的分类器。 分类器模型的参数选择是使用十折交叉验证结合网格寻优进行确定的,表3在SVM-acc列的括号处给出了本实验中SVM的具体参数,核函数均使用RBF。 (2)单独使用光谱特征和纹理特征的分类结果都不够理想,而将两种特征进行融合可以获得更好的分类效果,其中将波段降维后的光谱特征和纹理特征进行融合的分类结果要优于全波段光谱与纹理特征的融合,证明了光谱降维的有效性。 (3)和高光谱图像相比,RGB图像的分类准确率偏低,证明了高光谱图像在猕猴桃软腐病分类中的有效性。 此外,RGB图像光谱值与本文提出的波段融合图像的纹理特征融合后也能得到85.23%的准确率,说明本文通过提取特征波段图像融合后的纹理信息能够有效强化猕猴桃软腐病的病害特征,从而提高了模型的预测能力。

使用光谱特征、 RGB图像光谱特征以及融合特征进行SVM分类的混淆矩阵结果如图7所示。 可以看出,单独使用光谱特征对健康猕猴桃和软腐病早期猕猴桃进行分类的准确率不高,并且不能区分发病早期和发病晚期的样本,这与前文的光谱分析结果一致,发病早期和发病晚期样本光谱曲线比较接近,仅使用光谱特征对两类猕猴桃难以区分。 而使用RGB图像的光谱特征进行分类时,不同类别样本的分类结果均不理想,进一步说明RGB图像不适合分类研究。 当使用融合特征进行分类时,各个类别的误判率最低,大多数的猕猴桃软腐病发病早期样本可以得到正确识别,说明本方法能够实现猕猴桃软腐病的早期无损检测,同时在猕猴桃感染软腐病第三天或第四天病斑还未显现之时,通过提取光谱特征与融合图像的纹理特征的方法能够将软腐病早期猕猴桃与健康猕猴桃分离出来,有效减少了成本损失。

图7 模型测试结果(a): 光谱特征; (b): RGB光谱特征; (c): 融合特征Fig.7 Model test results(a): Spectral features; (b): RGB spectral features; (c): Fused features

本研究还对波段图像融合及纹理特征提取过程中的其他情况进行了对比实验,分别是(1)高频子带系数融合采用绝对值取最大方法,记作max-fusion。 (2)纹理特征提取时,分别计算不同角度的特征,每幅图片提取20个纹理参数代表该图像的纹理特征,记作four-angle。 (3)对不进行融合的8个特征波段提取纹理特征得到40维的向量作为图像特征,记作non-fusion。 (4)使用SPA选择的三个波段进行融合,记作SPAband-fusion。 实验结果如图8所示。 max-fusion实验准确率均偏低,与1.4.3节分析一致,采用绝对值最大法融合后提取的纹理特征,高频子带存在一些图像失真以及信息丢失,该纹理特征无论是和哪种光谱特征进行融合,结果均劣于加权平均法融合后所提取的纹理特征。 在four-angle与non-fusion两个实验中,纹理特征与SPA光谱特征融合的分类结果也偏低,这是因为SPA所得光谱特征仅有3维,而两种方法所得纹理特征分别为20、 40维,纹理特征存在冗余,特征融合后削弱了光谱特征在分类中的作用。 SPAband-fusion实验则表明,本工作所选的8个特征波段包含了更多的图像特征,因此结果要优于SPA三个波段的融合。

图8 对比实验结果Fig.8 Comparative experimental results

2.2 与其他方法对比

为了更全面的说明本工作在波段融合后提取纹理特征的有效性,本节与一些相关文献中的图像特征提取方法进行了对比,使用的数据集均为本工作的猕猴桃软腐病数据集,统一使用SVM分类器进行分类。 结果如表4所示,可以看出本方法优于其他方法。 其中文献[8]分别提取了6维颜色特征,16维GLCM特征以及32维直方图特征; 文献[9]使用光谱特征与主成分图像提取的纹理特征融合后进行分类; 文献[19]根据PC1权重曲线波峰波谷获取特征波段,并将特征波段图像的平均灰度值作为图像特征; 文献[20]在SPA获得的波段中选择了一个最为清晰、 亮度均匀的波段图像,利用此图像提取了共16维GLCM特征。

表4 本方法与其他方法的对比Table 4 Comparison between this method and other methods

在进行图像特征提取时,类比的众多文献大多选择主成分图像或特征波段图像进行提取,本工作与文献[9]虽然使用的特征变量都是8维,但主成分图像由于造成了光谱失真导致主成分图像提取的纹理特征对分类的效果要劣于图像融合后提取的纹理特征。 文献[20]虽然使用的特征变量是19维,但是其中有16维是图像特征,削弱了光谱特征的分类能力,此外,仅使用一个波段图像提取的纹理特征也不够充分,故分类准确率较低。 而其他文献中输入分类器的特征维度都较高,特征中可能存在着冗余特征,因此分类结果较差。

3 结 论

以湖北省武汉市“云海一号”猕猴桃为研究对象,通过对健康猕猴桃及感染软腐病的不同时期猕猴桃进行高光谱图像采集,提出了一种特征波段图像融合的猕猴桃软腐病早期分类检测方法。 首先,对样本原始光谱分别使用PCA,SPA,CARS进行光谱特征提取得到光谱特征,然后对SPA求解过程中所得的8个特征波段图像使用NSCT进行图像融合并提取融合图像的纹理特征,最后通过融合SPA光谱特征与纹理特征建立SVM模型对软腐病早期分类效果达到最优,主要得出以下结论:

(1)健康猕猴桃和软腐病猕猴桃光谱差异明显,可以利用光谱特征将两者区分,但是猕猴桃在软腐病早期和晚期的光谱差异较小,仅利用光谱信息难以区分。

(2)相比于单一光谱特征和图像特征,将光谱特征与图像特征进行融合可以显著提高分类准确率,说明图像的纹理特征可以在一定程度上弥补猕猴桃软腐病在光谱特征上的分类缺陷。 通过图像融合后再提取纹理特征,可以充分利用高光谱图像“空谱合一”的优点,实现在降低纹理特征数据维度的同时,进一步提高分类准确率。

(3)本研究使用高光谱成像技术能够在猕猴桃感染软腐病3～4 d时将染病果与健康果成功区分,实现了猕猴桃软腐病的早期无损检测,为猕猴桃的销售分级提供了一定的指导意义。