采用的广陈皮出产于广东省江门市的新宝堂陈皮有限公司(网购), 样本均为茶枝柑的干燥成熟果皮, 样品陈化年份分别为5年、 10年、 15年、 20年, 不同年份的广陈皮在外观上很难进行人工区分。

采用的高光谱采集设备是北京卓立汉光仪器有限公司的GaiaField便携式高光谱成像系统, 如图1所示。 该系统主要由面阵CCD探测器, 高光谱镜头(分辨率696× 256 Pix), 均匀光源(4个100 W卤素灯, 位于采样台四周), 光谱仪(波长范围385~1 014 nm, 共256个波段)等组成。

图1

Fig.1

	Figure Option ViewDownloadNew Window
	图1 高光谱成像系统结构Fig.1 Hyperspectral imaging system

由于广陈皮样本形状各异、 表面凹凸不平, 采集的高光谱图像会有较大的差异; 同时光谱采集过程中相机存在暗电流等噪声影响, 为了使采集效果达到最佳, 必须调整相机的镜头参数。 设置高光谱成像系统的相机曝光时间为27 ms, 镜头到样本的距离为42 cm, 平台移动速度为1.38 mm· s-1。^[4]

高光谱需要在黑暗密闭的环境中进行采集, 使用100 W卤素灯模拟太阳光, 调整卤素灯方位角为45° 且距离样本区垂直高度为40 cm。 采样前需对广陈皮样本擦拭去除表面的灰尘。 在相同条件下, 样本内侧朝向相机镜头放置于反射率接近于0的黑色绒布上, 光谱仪镜头尽量对准样本的中心位置。 共获得341个广陈皮样本的高光谱数据, 包括86个5年样本数据、 87个10年样本数据、 82个15年样本数据和86个20年样本数据。

数据分析处理软件包括: SpecVIEW、 ENVI5.3、 The Unscrambler X10.4、 Matlab 2016a、 Python3.10等。

使用SpecVIEW软件对广陈皮的高光谱样本进行镜头校准和反射率校准, 利用自带的镜头校准修复失真和透视; 利用标准反射板的反射率校正数据, 得到更准确的样本高光谱反射图像^[8]。 高光谱图像反射率校正如式(1)所示。

R=(R0-RB)/(RW-RB)(1)

式(1)中, R为校正后高光谱图像, R₀为原始图像, R_B为遮挡镜头得到的黑板图像, R_W为白板反射图像。

使用ENVI5.3软件对采集到的广陈皮高光谱数据进行处理。 在每个样本的四角位置和中间位置分别取点, 总计选取5个感兴趣区域(region of interest, ROI), 每个ROI选择的范围是20× 20 Pix, 因此每个样本可得到2 000(400× 5)条光谱数据, 将数据导入到Matlab 2016a中进行均值化处理。

图2

Fig.2

	Figure Option ViewDownloadNew Window
	图2 5年广陈皮样本ROI取样Fig.2 ROI sampling of 5 years old Cantonese tangerine peel samples

高光谱成像数据受到来自广陈皮自身因素干扰(样本表面平整度、 水分含量等)和其他外界因素干扰(环境、 暗电流等), 因此采集的样本数据中通常包含着无关信息和噪声^[9]。 为了降低以上干扰对数据采集的影响, 需要进行光谱预处理。

1.4.1 多项式平滑算法

多项式平滑算法(savitzky-golay smoothing, SG)是基于最小二乘原理的平滑算法, 其原理是利用中心点以及周边若干点进行最小二乘拟合, 尽可能减少噪音对有用信息的影响, 同时提高光谱数据的平滑性。

1.4.2 多元散射校正算法

多元散射校正算法(multiple scattering correction, MSC)通过理想光谱修正数据的基线偏移现象, 可以有效的消除不同散射水平的光谱差异, 从而增强光谱与数据之间的相关性, 分离散射光谱中的物理光散射信息以及化学光吸收信息, 校正样品中的光谱信息。

1.4.3 去趋势算法

去趋势算法(detrend)将光谱图象分解成静态稳定分量和非线性趋势项, 并消减非线性趋势项, 从而得到静态稳定分量, 其优势是减少相关性分析中数据畸变的概率, 防止低频成分上翘掩盖主频成分, 影响模型鉴别精度。

基于各类预处理算法的广陈皮高光谱数据如图3所示, 其中基于SG_detrend算法处理后的光谱曲线有相对明显的波峰和波谷, 因此更适合于处理广陈皮高光谱数据。

图3

Fig.3

	Figure Option ViewDownloadNew Window
	图3 基于各类预处理算法的广陈皮高光谱数据 (a): 原始光谱数据; (b): 基于SG算法处理后的光谱数据; (c): 基于SG_MSC算法处理后的光谱数据; (d): 基于SG_detrend算法处理后的光谱数据Fig.3 Hyper-spectral data of Cantonese tangerine peel using various preprocessing algorithms (a): Raw spectral data; (b): Spectral data using SG algorithm; (c): Spectral data using SG-MSC algorithm; (d): Spectral data using SG-detrend algorithm

高光谱具有波长多、 谱带重叠等特征, 在处理高光谱信息时可能出现运算资源的浪费和运算时间过长等问题。 同时无关的波段对建立模型还会产生负面影响, 所以需要通过光谱降维去除冗余波段, 这样既提高模型的鉴别准确率, 同时也降低模型过拟合风险和存储运算压力。

1.5.1 连续投影算法

连续投影算法(successive projections algorithm, SPA)是一种使矢量空间共线性最小化的前向变量选择算法, 在降低共线信息和获取有效信息等方面有较好的效果^[10]。 在全波段的光谱数据中选择出提供主要贡献的波段, 减少参与鉴别的光谱维度, 同时降低噪声和其他冗余信息的影响。

将基于SG_detrend算法处理后的广陈皮高光谱数据作为输入样本, 经SPA算法筛选特征波段后, 得到18个特征波段, 如图4所示。

图4

Fig.4

	Figure Option ViewDownloadNew Window
	图4 SPA算法筛选特征变量过程Fig.4 Process of selecting characteristic wavelength by SPA

1.5.2 竞争性自适应重加权算法

竞争性自适应重加权算法(competitive adaptive reweighted sampling, CARS)以波段变量回归系数的绝对值当参考指标, 选取数值大的变量, 筛除数值小的变量, 其中变量个数由衰减指数法(exponentially decreasing function, EDF)确定; 在未被选取的波段变量中, 采用自适应加权算法(adaptive reweighted sampling, ARS)再次选取回归系数绝对值相对较大的波段^[11]。

将基于SG_detrend预处理算法的广陈皮高光谱数据作为输入样本, 设置蒙特卡洛迭代次数为50; 结果表明, 当迭代次数为8时, 筛选效果最好, 经CARS算法筛选后得到50个特征波段, 如图5所示。

图5

Fig.5

	Figure Option ViewDownloadNew Window
	图5 CARS筛选特征变量过程Fig.5 Process of selecting characteristic wavelength by CARS

1.5.3 逐步回归算法

逐步回归算法(stepwise regression, SR)通过逐个引入贡献相对大的显著自变量来增加回归方程的准确性, 同时排除非显著变量, 减少冗余信息干扰, 最终得到最优的回归模型。 由于CARS算法筛选出的特征波段数量较大, 需要进一步降维, 采用SR算法对CARS算法筛选后的广陈皮高光谱特征波段进一步筛选, 最终得到22个特征波段, 如图6所示。

图6

Fig.6

	Figure Option ViewDownloadNew Window
	图6 SR算法筛选特征变量Fig.6 Characteristic variable selection by SR

分类鉴别模型是模式识别与机器学习等领域的重要组成部分, 其表征了模型输入与输出之间的特定映射关系, 能够通过对已有样本数据的学习与训练来实现对广陈皮高光谱数据的精准分类鉴别。

1.6.1 偏最小二乘判别分析模型

偏最小二乘判别分析(partial least-regression, PLS)模型通过将因变量和自变量投影到新的多维空间中, 形成基于最大化自变量相干性的线性回归模型。 PLS将原始光谱中包含的信息放置在排序在前面的若干主成分中, 所选择的主成分之间互不干扰, 这样既解决了高光谱数据中的共线性难题, 同时减少了噪声影响。

1.6.2 支持向量机模型

支持向量机(support vector machine, SVM)模型采用线性方程组得到支持向量, 基于最小化经验误差和最大化几何边缘区的方式来降低结构分类风险和过拟合风险, 进而增加分类模型的稳定性。 但是SVM模型的错误惩罚系数c和核系数g大多基于经验数据获取, 通用性较差, 因此本研究结合样本数据特征进一步优化SVM模型参数。

1.6.3 粒子群算法优化SVM模型

粒子群算法(particle swarm optimization, PSO)的灵感来自于生物体的社会行为, 利用一群粒子通过学习不断调整位置和速度来寻找最佳的个体和全局最优解。 设置PSO算法的种群规模为20, 迭代次数为100, 独立运行次数为20; 将SVM模型系数作为PSO算法的可行解域, 基于最优的适应度值更新粒子位置和速度, 采用当前的全局最优解训练SVM模型, 直至达到最大迭代次数。

1.6.4 蝗虫算法优化SVM模型

蝗虫算法(grasshopper optimization algorithm, GOA) 是一种元启发式仿生优化算法^[12], 模拟蝗虫在自然界中的种群迁移和觅食行为, 具有较高的搜索效率和较快的收敛速度, 并且算法本身特殊的自适应机制能够很好地平衡全局和局部搜索过程, 具有较好的效果。 设置GOA算法的种群规模为20, 迭代次数为100, 独立运行次数为20; 基于GOA算法优化SVM模型的过程与1.6.3节一致。

1.6.5 黑寡妇优化算法优化SVM模型

黑寡妇蜘蛛对能够表现个体状态的性信息素较为敏感, 高性信息素的配偶能够提供足够的生育能力, 低性信息素的个体因倾向于吞噬同类而被种群集体排斥和遗弃。 黑寡妇优化算法(black widow optimization algorithm, BWO)是基于黑寡妇蜘蛛独特行为启发的一种新型元启发式算法^[13], 具有收敛速度快、 寻优效率高等优势。 在BWO算法中, 黑寡妇蜘蛛能够根据上一次迭代过程中最优解所对应的位置和本次迭代过程中其他黑寡妇蜘蛛所在的位置来调整自身的下一个寻优位置。 基于BWO算法优化SVM模型的计算流程如图7所示, 具体步骤如下所述。

图7

Fig.7

	Figure Option ViewDownloadNew Window
	图7 BWO优化SVM流程Fig.7 Flow chart of SVM parameters optimized by BWO

(1)BWO算法参数: 种群规模为20, 迭代次数为100, 独立运行次数为20, 将作为黑寡妇蜘蛛的性信息素警戒值设为0.3, 个体螺旋式移动决定阈值设为0.7; SVM模型参数: 错误惩罚系数c和核系数g的取值范围设为[0.000 1, 1 000]。

(2)采用均方根误差(root mean squared error, RMSE)作为适应度函数fitness, 如式(2)所示。

fitness=argmin(RMSESVMpredict)(2)

对所有的黑寡妇蜘蛛对应的适应度值进行排序, 得到当前最优的SVM模型参数。

(3)计算并更新当前迭代过程中所有黑寡妇蜘蛛对应的性信息素值, 如式(3)所示。 其中fitness(i)是当前黑寡妇蜘纂的适应度, pheromone(i)是当前黑寡妇蜘蛛的性信息素值。

phermone(i)=fitnessmax-fitness(i)fitnessmax-fitnessmin(3)

(4)黑寡妇蜘蛛的移动策略分为线性移动和螺旋式移动, 如式(4)所示, 在[0, 1]范围内随机生成的变量c, 并通过判断其是否小于决定阈值来确定本次迭代过程中的移动策略。 其中 X→i(t)是当前的黑寡妇蜘蛛位置, X→i(t+1)是黑寡妇蜘蛛移动到的新位置, X→*(t)表示前一次迭代中搜索到的最佳黑寡妇蜘蛛, 变量m是位于区间[0.4, 0.9]中随机产生的浮点数, r₁是从1到当前黑寡妇蜘蛛序号之间随机产生的整数, X→r1(t)是在本次迭代过程中第r₁只黑寡妇蜘蛛, β 是定义在[-1, 1]区间内随机生成的浮点数。

X→i(t+1)=X→* (t)-mX→r1(t)X→* (t)-cos(2πβ)X→i(t)(4)

(5)低性信息素的黑寡妇蜘蛛往往是非常饥饿需要吞噬同类的雌性黑寡妇蜘蛛, 这类黑寡妇蜘蛛会受到种群集体排斥, 会被种群遗弃。 在迭代过程中, 低性信息素的黑寡妇个体被遗弃后, 需要及时补充种群数量, 如式(5)所示。 其中r₁、 r₂是从1到当前搜索代理蜘蛛范围内随机产生的不同整数, X→r1、 X→r2是第r₁、 r₂个黑寡妇蜘蛛, σ 是[0, 1]区间内随机生成的浮点数。

X→i(t)=X→* (t)+12[X→r1(t)-(-1)σX→r2(t)](5)

(6)计算并更新黑寡妇群体适应度值, 若达到算法停止条件则输出最优解, 否则重复步骤2— 步骤6。

采用PLS模型检验不同的高光谱数据预处理算法性能, 指标包括决定系数(determination coefficients, R²)、 RMSE和PLS模型鉴别准确率^[14]。 其中, R²值越大表征回归效果越好, RMSE越趋近于0则表示鉴别结果越靠近实际结果。 预处理算法性能结果如表1所示。 从表中可知, 各类高光谱数据预处理算法均对PLS模型的鉴别结果有一定程度的提升作用, 其中SG_detrend预处理算法的性能最佳。

表1

Table 1

表1(Table 1)

表1 预处理算法性能分析 Table 1 Performance analysis of preprocessing algorithms 建模方法 R2 RMSE 准确率/% PLS 0.962 96 0.196 45 98.83 SG-PLS 0.971 34 0.189 94 99.12 SG_MSC-PLS 0.969 99 0.194 35 99.41 SG_detrend-PLS 0.971 89 0.188 10 99.41

表1 预处理算法性能分析 Table 1 Performance analysis of preprocessing algorithms

针对基于SG_detrend预处理算法的广陈皮高光谱数据, 分别采用SPA和CARS_SR两类降维算法来筛选特征波段, 如表2所示, 其中SPA算法筛选得到18个特征波段, CARS_SR筛选得到22个特征波段。

表2

Table 2

表2(Table 2)

表2 广陈皮高光谱特征波段筛选分析 Table 2 Screening and analysis of hyperspectral characteristic bands of cantonese tangerine peel 降维算法 数量 特征波段/nm SPA 18 385.679 993 394.769 989 403.880 005 406.149 994 408.450 012 410.730 011 472.970 001 781.049 988 834.429 993 844.650 024 875.549 988 956.349 976 985.400 024 990.700 012 996.000 000 998.669 983 1 003.969 971 1 011.969 971 CARS_SR 22 403.880 005 410.730 011 422.179 993 489.299 988 512.750 000 538.700 012 560.099 976 584.099 976 634.799 988 659.200 012 681.299 988 703.530 029 718.400 024 748.349 976 781.049 988 834.429 993 875.549 988 898.849 976 911.849 976 959.000 000 980.099 976 956.349 976 998.669 983

表2 广陈皮高光谱特征波段筛选分析 Table 2 Screening and analysis of hyperspectral characteristic bands of cantonese tangerine peel

使用蒙特卡洛算法将各年份广陈皮高光谱数据按比例3:1随机分为训练集和测试集, 共得到训练集258个, 测试集83个。 错误惩罚系数c和核系数g是决定SVM模型性能的重要参数, 本研究分别采用了PSO算法、 GOA算法和BWO算法来优化SVM模型参数。 其中, 参数c和g的取值范围为[0.000 1, 1 000]; 各优化算法的种群规模为20, 迭代次数为100, 独立运行次数为20。 经各算法优化得到的SVM模型最优参数分别为: (PSO-SVM)c=59.931 6, g=0.010 3; (GOA-SVM)c=93.975 9, g=0.506 6; (BWO-SVM)c=80.270 1, g=0.159 5。

为综合分析不同分类鉴别模型性能, 本研究设计了包含PLS、 PSO-SVM、 GOA-SVM和BWO-SVM四种算法的广陈皮陈化年份鉴别实验。 其中, 输入样本数据均采用了SG_detrend预处理算法, 特征波段分别采用了未经降维算法处理的全波段(FULL)、 SPA特征波段集合和CARS_SR特征波段集合。 不同分类鉴别模型的实验结果如表3所示。 首先, 从降维算法的角度看, 与全波段相比, SPA算法和CARS_SR算法均能够在一定程度上剔除冗余信息, 进而提升鉴别模型的各项性能; 其次, 从分类鉴别模型的角度看, 与PLS模型相比, SVM模型能够得到更好的鉴别结果, 尤其是基于全波段的FULL-BWO-SVM模型对广陈皮年份鉴别的准确率(95.18%)高于基于CARS_SR降维算法处理后的SPA-PLS模型(91.57%); 最后, 从优化算法的角度看, 与PSO算法和GOA算法相比, BWO算法优化SVM模型的实验结果具有明显优势, 也验证了BWO算法寻优效率高、 不易于陷入局部最优等特点。 从表中可知, SG_detrend-CARS_SR-BWO-SVM模型性能最好, 对广陈皮陈化年份的鉴别准确率达到97.59%, RMSE为0.060 2, R²为0.952 9, 该误差主要来自高陈化年份的广陈皮样本数据, 验证结果如图8所示; 该模型比次优模型SG_detrend-SPA-BWO-SVM提升了1.20个百分点, 比对照模型CARS_SR-PSO-SVM提升了4.82个百分点。

表3

Table 3

表3(Table 3)

表3 不同分类鉴别模型的实验结果 Table 3 Experimental results of different classification discrimination models 建模方法 准确率/% RMSE R2 FULL-PLS 87.95 0.119 4 0.855 8 SPA-PLS 91.57 0.090 7 0.901 4 CARS_SR-PLS 92.77 0.083 2 0.920 3 FULL-PSO-SVM 89.16 0.216 9 0.837 5 SPA-PSO-SVM 91.57 0.120 5 0.909 3 CARS_SR-PSO-SVM 92.77 0.072 3 0.943 1 FULL-GOA-SVM 89.16 0.216 9 0.838 6 SPA-GOA-SVM 93.98 0.096 4 0.926 8 CARS_SR-GOA-SVM 93.98 0.093 7 0.928 5 FULL-BWO-SVM 95.18 0.084 3 0.935 1 SPA-BWO-SVM 96.39 0.072 3 0.943 5 CARS_SR-BWO-SVM 97.59 0.060 2 0.952 9

表3 不同分类鉴别模型的实验结果 Table 3 Experimental results of different classification discrimination models

图8

Fig.8

	Figure Option ViewDownloadNew Window
	图8 SG_detrend-CARS_SR-BWO-SVM模型验证结果Fig.8 SG_detrend-CARS_SR-BWO-SVM model validation