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

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

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

图1

Fig.1

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

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

图2

Fig.2

	Figure Option ViewDownloadNew Window
	图2 总体实验流程Fig.2 Overall experimental process

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

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

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

图3

Fig.3

	Figure Option ViewDownloadNew Window
	图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

	Figure Option ViewDownloadNew Window
	图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

表1(Table 1)

表1 光谱数据未处理和经预处理的结果对比 Table 1 Results of input data with different preprocessing methods 预处理方法 准确率/% whole bands PCA SPA CARS SG卷积平滑 72.72 75 68.18 75.00 无预处理 77.27 75 72.73 77.27

表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

	Figure Option ViewDownloadNew Window
	图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

表2(Table 2)

表2 特征波段选择结果 Table 2 Feature band selection results 波段选择方法 所选特征波段/nm SPA 548、 649、 692 CARS 473、 496、 502、 508、 513、 516、 525、 531、 554、 557、 565、 583、 585、 591、 609、 611、 617、 620、 632、 634、 643、 660、 663、 666、 695、 707、 712、 715、 718、 724、 727、 736、 764、 793、 805、 808、 822、 848、 868

表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个纹理参数代表该图像的纹理特征。

虽然特定波长下图像对于提取的纹理特征有更好的描述, 但8个波段图像提取的纹理特征有40维, 一定程度上会存在冗余。 因此, 采用非下采样轮廓波变换(non subsampled contourlet transform, NSCT)对波段图像进行融合。 NSCT主要过程是: 首先使用非下采样金字塔形滤波器对图像进行多尺度分解得到多个子带系数, 然后使用非下采样方向滤波器组对各子带系数进行方向分解, 最终得到不同尺度、 不同方向的子带。 一幅图像进行J级NSCT变换后, 生成1个低频子带和 ∑j=1J2lj个高频子带, 其中 lj表示在第j级分解上进行方向分解的数量。设 Lb1, Lb2, …, Lb8分别代表8个特征波段图像的低频子带, Hb1j, k, Hb2j, k, …, Hb8j, k代表8个特征波段图像在第j级尺度变换上第 k个方向的高频子带。 由于在分解过程中没有进行下采样, 所有子带都和原始图像大小相同, 因此该变换消除了各子带间的频谱混叠。 同时多尺度、 多方向的特性使图像的高频子带拥有更多的细节信息, 融合后的图像既能去除冗余信息, 又可以将各波段的互补信息合并到单张图像中。

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

LF=∑i=18Lbi8(1)

HFj, k=∑i=18Hbij, k8(2)

式(1)和式(2)中, L_F和 HFj, k分别代表子带系数融合后的低频子带和高频子带。 图6分别展示了第132波段的其中一个高频子带、 采用绝对值最大法融合的高频子带以及采用加权平均融合的高频子带。 可以看出, 采用绝对值最大法融合后的高频子带存在一些图像失真以及信息丢失, 此外, 原始波段图像中还存在由于传感器电路偏置电压差异而形成的条纹噪声, 而采用加权平均对高频子带进行融合后, 可以在一定程度上消除条纹噪声的干扰。

图6

Fig.6

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

对L_F和 HFj, k使用NSCT逆变换得到波段融合后的图像B_F, 再对B_F进行纹理特征的提取, 最后得到一个维度为5的特征向量, 即为猕猴桃样本的图像特征。

1.4.4 模型建立与评估

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

Xnorm=X-XminXmax-Xmin(3)

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

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

accuracy=TP+TNTP+TN+FP+FN(4)

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

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

表3

Table 3

表3(Table 3)

表3 不同组合下测试集分类结果 Table 3 Classification results of test sets under different combinations 模型输入 光谱处理方法 变量数 SVM-acc/% RF-acc/% KNN-acc/% HSI光谱特征 whole bands 150 77.27(c=1 000, g=0.1) 79.77 63.64 PCA 10 75(c=10 000, g=0.01) 77.27 69.32 SPA 3 72.73(c=10 000, g=0.1) 71.14 70.45 CARS 39 77.27(c=1 000, g=1) 76.82 70.45 HSI纹理特征 - 5 72.73(c=10 000, g=0.1) 64.32 67.05 HSI光谱特征+纹理特征 whole bands 155 89.77(c=10 000, g=0.001) 82.73 77.27 PCA 15 89.77(c=10 000, g=0.001) 86.82 77.27 SPA 8 92.05(c=10 000, g=0.01) 82.05 87.5 CARS 44 92.05(c=10 000, g=0.01) 82.73 79.55 RGB光谱特征 - 3 62.5(c=10 000, g=0.1) 47.5 47.73 RGB纹理特征 - 5 77.27(c=10 000, g=0.01) 67.95 64.77 RGB光谱+纹理 - 8 77.27(c=10 000, g=0.01) 73.41 64.77 RGB光谱+HSI纹理 - 8 85.23(c=10 000, g=0.01) 79.55 73.86

表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

Fig.7

	Figure Option ViewDownloadNew Window
	图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

	Figure Option ViewDownloadNew Window
	图8 对比实验结果Fig.8 Comparative experimental results

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

表4

Table 4

表4(Table 4)

表4 本方法与其他方法的对比 Table 4 Comparison between this method and other methods 文献 输入变量数 准确率/% 本方法 8 92.05 [8] 204 88.64 [9] 8 88.64 [19] 155 78.41 [20] 19 85.23

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

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