实验所用的大白菜样本均采购于哈尔滨市某大型商超。 实验农药为40%毒死蜱(市售, 河北八源生物制品有限公司)、 58.00 mg· L^-1乐果(网购, 中国计量科学研究院)、 90%灭多威(市售, 山东华阳农药化工集团有限公司)和4.5%氯氰菊酯(市售, 安徽尚禾沃达生物科技有限公司)。 按照《中华人民共和国食品安全国家标准GB2763— 2019》中大白菜最大农药残留限量要求, 分别将毒死蜱、 乐果、 灭多威和氯氰菊酯分别配置0.10, 1.00, 0.20和2.00 mg· kg^-1的药溶液分别装入4个喷壶内, 喷射距离8~10 cm, 每种药溶喷洒10个样本, 共计40个样本, 均匀喷后平放置于室温(温度20 ℃, 湿度65%RH)自然吸收12 h后采集高光谱影像数据(共50个样本, 其中10个为无农药残留样本)。

高光谱采集设备使用美国Headwall Photonics公司生产的高光谱成像系统, 由面阵CCD、 光栅光谱仪、 高光谱成像镜头、 均匀光源、 一维电动台、 USB1394图像采集卡、 高性能图形工作站及相关采控软件构成。 实验所用高光谱相机分辨率为1 000× 164 Pix, 位深度为24 bit, 线性阵列扫描成像方式; 光栅光谱仪的光谱范围为400~1 000 nm, 共有800个波段, 光源由2个200 W溴钨灯构成, 位于电动台两侧, 入射角呈45° ; 采控软件通过RS232串口控制电动台移动速度及方向, 如图1所示。

图1

Fig.1

	Figure Option ViewDownloadNew Window
	图1 高光谱成像系统结构示意图Fig.1 Schematic diagram of hyperspectral imaging system

高光谱图像采集过程中会受到两方面干扰, 一方面是由于大白菜样本形状各异、 表面高低不平等原因会产生漫反射; 另一方面是由于电网中暗电流和谐波的存在, 导致相机工作状态不稳定, 产生较大的图像噪声, 影响光谱采集效果。 因此, 在采集前必须对高光谱相机进行参数设置和黑白校正, 经多次调试和效果比对, 最终曝光时间设置为0.03 s, 电动台的移动速度参数设置为3.0 mm· s^-1, 镜头垂直向下距电动台450 mm。

后续所使用的数据分析处理软件包括: ENVI 5.2(ITT Visual Information Solutions, Boulder, Co., > USA), OriginPro 8(OriginLab Co., Ma, USA), The Unscrambler X 10.4, Python 3.6, TensorFlow 2.0等软件。

为方便数据采集将大白菜样本裁剪为16 cm× 8 cm长方条, 使用ENVI软件在每个大白菜样本上避开主茎干选取40个感兴趣区域(region of interest, ROI), 选取位置如图2所示。 每个ROI选择的范围是30× 30 Pix, 计算出这900个像素点平均光谱反射值作为一条光谱记录, 共计2 000条(5× 400)光谱记录。

图2

Fig.2

	Figure Option ViewDownloadNew Window
	图2 大白菜样本ROI选取示意图Fig.2 Schematic representation of selection of ROI on Chinese cabbage sample

在高光谱图像采集过程中, 会受到来自样本自身因素干扰(如样本平整度、 色差和水分含量等)和环境因素干扰等, 为了减少上述干扰的影响, 需要对高光谱图像进行光谱预处理以消除干扰影响。 使用Unscrambler软件对高光谱数据进行多元散射校正(multiple scattering correction, MSC)预处理; 5组样本各随机选取200条光谱在预处理前与后的谱形如图3所示。 图3(a)为未经光谱预处理的原始样本曲线, 图3(b)为经MSC预处理后样本曲线。

图3

Fig.3

	Figure Option ViewDownloadNew Window
	图3 光谱数据预处理 (a): 原始光谱; (b): 经MSC处理Fig.3 Spectral data preprocessing by MSC (a): Without preprocessing; (b): With MSC preprocessing

如图3(b)所示, 靠近400和1 000 nm边缘波段光谱反射曲线波动较剧烈, 说明此处受到干扰较大, 数据严重失真, 影响后期建模分类效果; 因此, 需要进行噪声裁剪, 为每个样本掐头去尾剔除噪声较大400~414和912~1 000 nm的边缘波段, 保留415~911 nm(对应原800波段中的20~749, 共计729波段)用于建模分析。

小波变换(wavelet transform, WT)是一种信号变换分析方法, 提供一种随频率改变的“ 时间-频率” 窗口, 通过小波基函数的伸缩平移运算实现信号(函数)多尺度、 多分辨细化, 可同时提取细节分量(高频部分)和近似分量(低频部分), 连续小波变换表达式为

WT(j, k)=< f(λ), Ψj, k(λ)> =1|j|∫-∞+∞f(λ)Ψλ-kjdλ(1)

式(1)中, j为缩放因子, k为平移因子, j和k均为连续变量, 量, Ψ (λ )为小波基函数。

离散小波变换^[6](discrete wavelet transform, DWT)是小波变换的一种, 可实现在离散尺度上进行信号分解, 通过高低通滤波器来实现高低频信号分解, 离散变换公式为

DWT(a, b)=< f(λ), ϕa, b(λ)>(2)

式(2)中, a和b分别为第a层分解和第b个小波系数, ϕ _{a, b}(λ )为离散小波基函数。

将高光谱数据中的波长-反射率与信号的时间-频率相对应, 利用离散小波变换多尺度分解特性实现降维, 既每层变换均略去细节分量(高频部分), 将分解后的近似分量(低频部分)作为建模数据, 以此类推, 每层变换均在上层变换的近似分量基础上进行再变换, 每层变换数据维度减半, 达到数据降维目的。

1.7.1 多层感知机

所使用的多层感知机(multilayer percetron, MLP)^[7]具有4个全连接层和1个输出层, 网络结构如表1所示; 设置学习率learning_rate(LR)为0.001, 迭代次数epochs为500, 批量大小为batch_size(BS)为64, 测试集验证频率validation_freq为1, 优化器使用Adam算法。

表1

Table 1

表1(Table 1)

表1 多层感知机网络结构 Table 1 MLP network structure 层 神经元个数 激活函数 全连接层FC1 512 relu 全连接层FC2 256 relu 全连接层FC3 256 relu 全连接层FC4 128 relu 输出层 5 softmax 注: relu是线性修正单元激活函数 Note: relu is rectified linear unit

表1 多层感知机网络结构 Table 1 MLP network structure

1.7.2 卷积神经网络

卷积神经网络^[8](convolutional neural networks, CNN)是一种在多层感知机基础上建立的深度学习网络, 通过增加若干个卷积层(convolution layer)来增强网络的特征抓取能力。 虽然经典模型如AlexNet(8层)、 VGGNet(11~19层)、 GoogLeNet(Inception V1 22层)、 ResNet(18~152层)分类效果较好, 但模型层数较多、 参数量巨大, 对硬件性能要求较高。 本工作重新设计卷积神经网络结构如图4所示; 高光谱大白菜样本图像经ROI提取、 预处理和噪声裁剪后数据为2 000× 1× 1× 729, 为了便于二维卷积运算, 将光谱数据变形为2 000× 27× 27× 1(all_batch× width× height× channels), 卷积神经网络结构参数设置如表2所示。

图4

Fig.4

	Figure Option ViewDownloadNew Window
	图4 卷积神经网络结构Fig.4 CNN structure

表2

Table 2

表2(Table 2)

表2 卷积神经网络结构参数设置 Table 2 Parameter settings of the CNN structure 层 离散小波变换层数 1~4 5~6 卷积炽C1 3× 3× 16 3× 3× 16 激活层A1 relu relu 舍弃层D1 0.3 0.3 卷积层C2 5× 5× 32 3× 3× 32 正则层L21 0.001 0.001 激活层A2 relu relu 舍弃层D2 0.5 0.5 池化层MaxPooling 2× 2× 2 2× 2× 2 全连接层FC1 256 256 全连接层FC2 256 256 正则层L22 0.001〗〗0.001 输出层 5 5 激活层A3 softmax softmax 注: 卷积层n× n× m表示为m个卷积核大小为n× n的过滤器; 池化层n× n× m表示池化核大小为n× n, 步长为m Note: Convolution layer n× n× m represents m filters with the convolution kernel size of n× n; Pooling layer n× n× m represents the pooling kernel size of n× n and the step size of m

表2 卷积神经网络结构参数设置 Table 2 Parameter settings of the CNN structure

卷积神经超参数的选取直接影响到识别模型分类精度, 为了确定关键超参数Learning_rate, Batch_size和Epochs, 高光谱数据经多次、 多参数训练, 如图5所示; 图5(a)为Batch_size=50, Epochs=1 000时, 4种不同Learning_rate下的模型训练损失值, 从中可以看出当Learning_rate为0.005和0.01时均出现Training loss value(TLA)为不变值1.6, 说明此时Learning_rate过大, 导致模型无法收敛; 当Learning_rate为0.001时模型收敛速度较快。 图5(b)为Learning_rate=0.001, Epochs=1 000时, 6种不同Batch_size下的模型训练损失值, 中可以看出当Batch_size为70, 90, 110和130时Training loss value下降缓慢且振荡剧烈, 说明此时模型不稳定, 泛化能力较弱; 仅当Batch_size为50时Training loss value下降迅速且起伏较小。 图5(c)为Learning_rate=0.001, Batch_size=50时, 6种不同Epochs下的模型训练损失值, 从中可以看到随着Epochs不断增加, 总体精度(overall accuracy, OA)呈明显上升趋势, 当Epochs为1 000时, OA到达顶峰。 综合考量, 最终确定超参数Learning_rate为0.001, Batch_size为50, Epochs为1000, 优化器采用Adam算法。

图5

Fig.5

	Figure Option ViewDownloadNew Window
	图5 卷积神经网络超参数选择 (a): 基于4种Learing rate的模型训练损失值; (b): 基于6种Batch size的模型训练损失值; (c): 基于6种Epochs的模型总体精度Fig.5 CNN hyperparameter selection (a): TLA for different LR; (b): TLA for different BS; (c): OA for different Epochs

为了更全面比较模型性能和精度, 同时还采用支持向量机(support vector machine, SVM)和K最邻近(K-nearest neighbor, KNN)机器学习算法进行建模。

由于样本数量较少, 容易发生模型过拟合, 在建模过程中引入数据增强技术(data augmentation), 在原有的高光谱数据基础上增加高斯噪声, 将样本数据量扩充为原来的2倍; 除此之外, 还采用Dropout和Regularization L2技术提高模型泛化能力和鲁棒性。

为了从光谱曲线角度清楚地区分不同类型的样品之间的光谱信息的差异, 将每类样品的400条光谱曲线取平均值, 得到该平均光谱曲线如图6所示。

图6

Fig.6

	Figure Option ViewDownloadNew Window
	图6 大白菜样本平均光谱曲线Fig.6 Average spectra of chinese cabbage samples

由图6可以看出, 在400~700 nm波段内, 出现3次较为明显反射波谷, 分别出现在410 nm紫光、 490 nm青光和667 nm红光附近; 在667~710 nm波段内反射率陡然上升, 710 nm波段后反射率逐渐平稳, 符合绿色叶片光谱反射率规律。 在530~667 nm波段内5条光谱反射曲线区分度较高, 印证了叶片光谱反射率与农药胁迫存在相关联理论^[9], 为实现通过光谱对大白菜农药残留进行识别提供了科学依据。

2.2.1 离散小波变换降维

为了比较降维效果, 以db1小波基为例, 在平均光谱曲线上进行离散小波变换, 每次变换后均只保留低频分量, 舍弃高频噪声分量, 1~6层离散变换降维后得到的数据维度分别为364, 182, 91, 45, 22和11。

降维后的光谱曲线如图7所示, 可以看到随着变换层数的增加, 曲线形状与原始曲线差异逐渐变大; 图7(a)为经过1层离散小波变换后的低频分量曲线, 与图6相比具有较高的相似性, 能够较好的表征原曲线之间相对位置关系, 但所用数据维度只占原数据的1/2; 随着离散小波变换层数的增加, 数据维度以2-m递减(m为小波变换层数), 图7 (f)为经过6层离散小波变换后的曲线, 虽然数据维度只有11维, 仍能够看出在关键的节点处曲线相对位置关系。

图7

Fig.7

	Figure Option ViewDownloadNew Window
	图7 基于db1小波基函数离散小波变换后的低频部分曲线 (a)— (f)分别经1~6层离散小波变换Fig.7 Low frequency portions of wavelet transform based on db1 function (a)— (f) corresponding to 1~6 layers of DWT, respectively

2.2.2 PCA和CARS降维

采用主成分分析(principal component analysis, PCA)算法对所有裁剪后的大白菜高光谱数据进行降维, 前11, 22, 45, 91, 182和364主成分累计贡献率分别为99.11%, 99.48%, 99.76%, 99.91%, 99.98%和99.99%。

采用竞争性自适应重加权算法(competitive adaptive reweighted sampling, CARS)进行全样本降维, 其中, 采样频率设置为0.8, 筛选阈值为0.8, 蒙特卡洛采样次数为300, 交互验证组为15, 筛选过程如图8。

图8

Fig.8

	Figure Option ViewDownloadNew Window
	图8 CARS筛选特征波长过程Fig.8 Process of selecting characteristic wavelength by CARS

在图8(a)和(b)中可以看到随着运行次数的增加, 特征波长数量逐渐减少, RMSECV呈现先下降后上升的趋势, 当在140次运行时RMSECV到达最小值; 图8(c)表示特征波长变量回归系数的变化趋势, “ * ” 表示RMSECV最小的位置, CARS共筛选出特征波长39个, 占全波段的6.12%, 特征波长分别为430.04, 445.81, 448.06, 448.81, 449.56, 472.09, 472.84, 509.64, 534.42, 538.17, 625.28, 626.78, 627.53, 703.38, 749.19, 749.94, 755.19, 756.70, 757.45, 767.21, 767.96, 797.25, 804.76, 805.51, 807.01, 812.27, 823.53, 830.29, 831.04, 831.79, 832.54, 833.29, 840.05, 843.05, 873.09, 892.62, 909.89, 928.66, 936.17 nm。

2.3.1 基于多种小波基离散变换降维模型评价

不同的小波基函数降维效果不同, 为了获得最佳分类识别效果, 选取Daubechies小波族中的db1, Symlets中的sym2, Coiflets中的coif1, Biorthogonal中的bior2.2和ReverseBior中的rbio1.5分别对数据进行1~6层小波变换后建立识别模型, 建立流程如图9所示。 样本划分采用Kennard-Stone算法, 样本总数的75%作为训练集, 25%作为预测集。

图9

Fig.9

	Figure Option ViewDownloadNew Window
	图9 基于高光谱离散小波变换的卷积神经网络分类流程图Fig.9 Flow chart of hyperspectral image classification based on DWT and deep learning

4种算法的总体精度如图10(a)— (d)所示。 总体来看, 基于CNN和MLP算法的总体精度曲线平稳度均高于基于KNN和SVM算法, 说明CNN和MLP算法学习与预测能力较强, 各算法总体精度随离散变换层数的增加成波浪式上下浮动, 经bior2.2降维后的总体精度高于其他。 图10(a)中除db1和rbio1.5小波基外, 其余小波基降维后基于CNN算法建模总体精度在78.40%~91.20%之间波动; 其中, 经coif1第2层变换时总体精度取得最优值为91.20%; 图10(b)中经coif1和rbio1.5第4层变换后MLP算法总体精度均取得最优值为83.20%; 图10(c)中经sym2第2层变换后KNN算法总体精度取得最优值为66.40%; 图10(d)中经bior2.2第2层变换后SVM算法总体精度取得最优值为90.40%。

图10

Fig.10

	Figure Option ViewDownloadNew Window
	图10 基于DWT, PCA和CARS降维的建模结果 (a)— (d)分别为基于DWT降维的CNN, MLP, KNN和SVM模型总体精度; (e)基于PCA降维4种模型总体精度; (f)基于CARS降维4种模型总体精度Fig.10 Modeling results based on DWT, PCA and CARS (a)— (d) Overall accuracies of CNN, MLP, KNN and SVM models based on DWT; (e) Overall accuracies of four models based on PCA; (f) OAs of four models based on CARS

2.3.2 基于PCA和CARS降维模型评价

为充分验证本研究所提出的离散小波变换降维算法的可靠性, 分别取PCA降维后前364, 182, 91, 45, 22和11个主成分(分别对应离散小波1~6层变换后的数据维度)进行建模, 4种算法的总体精度如图10(e)所示。 随着数据维度的减少, CNN和MLP算法总体精度波动不大, 但KNN和SVM算法的总体精度波动较大, 相比而言CNN和MLP算法学习能力相对较强, 模型稳定性较好。 SVM算法在取前45个主成分时获得最优总体精度为62.80%。 经CARS降维后4种算法的总体精度如图10(f)所示, CNN取得该降维算法最高总体精度82.40%, KNN算法总体精度最低为11.20%, 可见KNN算法不适合该类数据样本。

2.3.3 不同降维和建模算法横向评价

3种降维和4种建模算法总体精度和Kappa系数如表3所示。 从整体来看在总样本未经降维过程的分类模型预测用时最长, 精度不高; 基于CARS降维算法分类总体精度最低; 经DWT降维后的分类模型KNN, SVM, MLP和CNN总体精度分别为66.40%, 90.40%, 83.20%和91.20%。 此外, 经DWT降维后4种分类算法预测用时略高于其他2种降维算法。

表3

Table 3

表3(Table 3)

表3 不同算法的总体精度和Kappa系数 Table 3 Overall accuracies and Kappa coefficients of different algorithms 模型 降维算法 None CARS PCA DWT OA /% Kappa Time /ms OA /% Kappa Time /ms OA /% Kappa Time /ms OA /% Kappa Time /ms KNN 12.00 -0.10 16.98 11.20 -0.11 10.00 56.00 0.49 4.00 66.40 0.58 20.02 SVM 89.60 0.87 74.77 26.40 0.08 4.99 68.00 0.60 4.99 90.40 0.88 14.03 MLP 67.20 0.59 80.78 72.80 0.66 59.99 66.40 0.58 67.01 83.20 0.79 63.23 CNN 82.40 0.78 698.13 26.40 0.08 4.99 64.00 0.55 94.01 91.20 0.89 86.01 注: KNN分别在PCA(11)和DWT(sym2-2)取得最优值, SVM分别在PCA(22)和DWT(bior2.2-2)取得最优值, MLP分别在PCA(45)和DWT(coif1-4)取得最优值, CNN分别在PCA(22)和DWT(coif1-2)取得最优值; 其中, PCA(n)表示基于PCA降维算法前n个主成分建模, DWT(W-g)表示基于W小波基离散第g层变换降维算法建模 Note: KNN reaches optimal values with PCA (11) and DWT (sym2-2); SVM reaches optimal values with PCA (22) and DWT (bior2.2-2); MLP reaches optimal values with PCA (45) and DWT (coif1-4); and CNN reaches optimal values with PCA (22) and DWT (coif1-2); PCA (n) denotes the modeling based on the first n principal components, and DWT (W-g) represents the modeling based on the dimension reduction algorithm using discrete transform of g(th) layer W wavelet basis

表3 不同算法的总体精度和Kappa系数 Table 3 Overall accuracies and Kappa coefficients of different algorithms

不同算法的分类精度如表4所示, 从整体来看灭多威和乐果分类精度较低(分类精度分别为92.00%和96.00%), 其原因可能是两种药物分子式均含有类似的甲基, 导致光谱特征区分度较低; 毒死碑和氯氰菊酯分别在CARS-SVM和DWT(bior2.2-2)-SVM算法取得最佳分类精度100%, 说明SVM识别分类能力较好; 乐果和无残留均在DWT(coif1-2)-CNN获得最高分类精度100%, 说明DWT-CNN算法有较好的降维、 特征提取和分类识别能力。

表4

Table 4

表4(Table 4)

表4 不同算法的分类精度 Table 4 Overall accuracies of the various prediction algorithms 降维 算法 模型 药物种类 毒死蜱 氯氰菊酯 灭多威 乐果 无残留 None KNN 16.00 12.00 8.00 24.00 0.00 CARS 16.00 12.00 4.00 20.00 4.00 PCA 58.82 42.00 63.79 44.44 69.56 DWT 60.00 80.00 48.00 64.00 80.00 None SVM 96.00 92.00 84.00 84.00 92.00 CARS 100.00 8.00 20.00 4.00 0.00 PCA 68.62 66.00 70.69 66.66 67.39 DWT 96.00 100.00 80.00 84.00 92.00 None MLP 88.00 88.00 52.00 12.00 96.00 CARS 56.00 52.00 92.00 72.00 92.00 PCA 60.78 68.00 75.86 46.67 78.26 DWT 88.00 84.00 72.00 76.00 96.00 None CNN 76.00 84.00 76.00 80.00 96.00 CARS 76.00 80.00 76.00 80.00 100.00 PCA 60.78 74.00 63.79 44.44 76.08 DWT 88.00 88.00 84.00 96.00 100.00

表4 不同算法的分类精度 Table 4 Overall accuracies of the various prediction algorithms