1.1.1 材料

本实验对象为采购于菜市场的新鲜小白菜, 将其以45 ℃温水洗净后作为施以农药的原材料。 选取高效氯氟氰菊酯为实验农药, 使用喷雾方式确保农药喷洒均匀。 将实验原材料分为四组, 第一组不喷洒农药, 作为对照实验组; 第二、 三、 四组的农药配比分别为1:500, 1:100, 1:20。 阴干后将白菜分解成小块, 便于进行中红外光谱数据的收集。

1.1.2 光谱仪器与分析软件

仪器选用安捷伦Cary 630 FTIR光谱仪, 调至ATR衰减全反射光谱采集模式。 扫描背景与样本各64次, 分辨率调整为8 cm^-1。 配合Microlab PC、 Resolutions Pro软件, 每组样品采集40组数据, 共计160组, 每组数据有971维, 波数变化范围为4 300~590 cm^-1。 白菜傅里叶中红外光谱(Fourier transform mid-infrared spectroscopy, FT-MIR)如图1所示, 采用Matlab R2016b对该光谱数据进行定性分析。

图1

Fig.1

	Figure Option ViewDownloadNew Window
	图1 原始傅里叶变换中红外光谱数据图Fig.1 Raw FT-MIR spectral data

(1)参数初始化

设置权重指数m(m> 1), 测试样本的维数d, 循环迭代次数的初始值r=0, 最大允许迭代次数r_max=100, 迭代的最大误差参数ε =0.000 01; 运行FCM算法并将最终模糊隶属度与类中心作为GKIPCM算法的初始模糊隶属度 UGKIPCM(0)与初始聚类中心 VGKIPCM(0)。

(2)计算第r次迭代时的模糊散射矩阵S_{fi, r}

Sfi, r=∑k=1nuik, r-1m(xk-vi, r-1)(xk-vi, r-1)T∑k=1nuik, r-1m(1)

式(1)中, u_{ik, r-1}为第r-1次迭代后第k个测试样本对第i类的模糊隶属度; n为样本总量; x_k为第k个白菜中红外光谱测试样本; v_{i, r-1}为第r-1次迭代后第i类的类中心(i=1, 2, 3, 4); S_{fi, r}是第r次迭代后第i类样本的模糊散射矩阵。

(3)计算第r次迭代后第i个聚类中心的范数矩阵A_{i, r}

Ai, r=‖Sfi, r‖1d(Sfi, r-1)(2)

(4)由式(3)计算第r次迭代后测试样本x_k到类中心v_{i, r-1}的距离范数D_{ik, r}

Dik, r=(xk-vi, r-1)TAi, r(xk-vi, r-1), ∀i, k(3)

(5)计算第r次迭代后第k个测试样本隶属于第i类的典型值t_{ik, r}

tik, r=exp-Dik, r2‖Sfi, r‖1d(4)

(6)计算第r次迭代后的模糊隶属度矩阵U_r, 其各元素值计算公式为:

uik, r=∑j=1c‖Sfi, r‖1d(1-tik, r)‖Sfj, r‖1d(1-tjk, r)1m-1-1(5)

式(5)中, c为样本类别数。

(7)计算第r次迭代后第i类的类中心v_{i, r}

vi, r=∑k=1nuik, rmtik, rxk∑k=1nuik, rmtik, r(6)

(8)计算相邻两次迭代后两聚类中心的距离范数D_r

Dr=‖vi, r-vi, r-1‖(7)

(9)对已计算数据进行如下判断:

若D_r< ε 或r≥ r_max, 则停止运行, 记录迭代结束时模糊聚类中心V_end, 模糊隶属度矩阵U_end以及典型值矩阵T_end; 否则令v_{i, r-1}=v_{i, r}, u_{ik, r-1}=u_{ik, r}, 回到步骤(2)继续迭代计算。

(10)根据模糊隶属度矩阵与典型值矩阵对样本进行分类:

对于模糊隶属度, 若u_ij为u_j中的最大值, 则判断第j个样本属于第i类; 对于典型值, 若t_ij为t_j中的最大值, 则认为第j个样本有更大的可能性属于第i类。

将样本数据分为训练集与测试集。 其中训练集共4类, 每类22个样本; 测试集共4类, 每类18个样本。 在光谱数据采集的过程中, 散射水平的差异往往造成白菜表面实际农药的残留量与光谱波长的数据相关性的下降。 为了减小此类差异且提高信噪比, 首先使用多元散射矫正(multiplicative scattering correction, MSC)对光谱数据进行预处理。 图2为使用多元散射矫正后傅里叶中红外光谱数据图。

图2

Fig.2

	Figure Option ViewDownloadNew Window
	图2 多元散射矫正后傅里叶中红外光谱数据图Fig.2 FT-MIR spectral data processed by MSC

为观察不同农药残留程度的光谱数据的吸光度有所差异, 计算了每一类光谱数据的平均值并截取了波数介于800与1 500之间的部分如图3所示。

图3

Fig.3

	Figure Option ViewDownloadNew Window
	图3 多元散射矫正后白菜平均中红外光谱图Fig.3 The average FT-MIR spectral data processed by MSC

不同农药浓度样品的中红外光谱吸光度有着明显差值, 直接保证了不同类光谱数据的可分性。

尽管MSC将基线平移并且消除了偏移量, 但数据依然达到了971维。 为减小运算量, 采用主成分分析(principal component analysis, PCA)对数据实现降维并提高分类准确率。 通过训练数据集的协方差矩阵得到前23个主成分, 累计贡献率高达99.60%。 故通过该23个主成分将光谱数据缩小至23维。 为了提取样本数据的相关特征信息并进一步降维, 使用线性判别分析(linear discriminant analysis, LDA)处理训练样本。 此处求取了3个判别向量与特征值。 其中特征值的计算结果为: λ ₁=7.371 3, λ ₂=3.932 0, λ ₃=1.931 7。 将PCA处理后的训练数据投影到判别向量构成的特征空间, 得到线性判别分析得分图, 见图4。

图4

Fig.4

	Figure Option ViewDownloadNew Window
	图4 线性判别分析得分图Fig.4 Scores plot of LDA algorithm

图中, “ · ” , “ * ” , “ ○” , “ × ” 分别表示农药配比为无农药, 1:500, 1:100, 1:20的光谱数据。 无农药的样本与配比为1:20的样本区分度最高, 相比之下农药配比为1:100和1:500的两类样本区分度稍弱。 尽管如此, 四类样本整体的区分度较强, 说明线性判别分析的特征提取结果是可取的。 于是将PCA的23组线性变换(主成分)与LDA的3个鉴别向量应用于测试集, 为模糊聚类做准备。

2.2.1 模糊聚类相关参数的初始化

设置权重指数m=2.0(m≥ 1), 测试样本的维数d=3; 由于FCM, GK, IPCM, GKICPM算法都属于迭代算法, 设置迭代次数初始值r=1, 最大迭代次数r_max=100; 当迭代误差小于ε =0.000 01或r≥ r_max时停止迭代。 FCM的初始聚类中心如式(8)

VFCM(0)=vFCM, 1(0)vFCM, 2(0)vFCM, 3(0)vFCM, 4(0)=-0.01280.0036-0.0024-0.0123-0.00160.0039-0.00620.00360.0053-0.00790.00050.0022(8)

对于GK, IPCM, GKIPCM, 将FCM的最终聚类中心与模糊隶属度作为算法的初始值。

2.2.2 模糊聚类算法的迭代时长

运行GK, IPCM, GKIPCM, 对比迭代次数与收敛时间: GK迭代了100次, 运行时间为0.093 8 s; IPCM迭代了13次, 运行时间为0.062 5 s; GKIPCM迭代了87次, 运行时间为0.218 8 s。 处理器: Inter(R) Core(TM) i5-8300H CPU @2.30GHz(8 CPUs); 内存: 8192MB RAM。

2.2.3 模糊隶属度分类结果

参照2.2.1将各参数初始化后, 分别运行GK, IPCM与GKIPCM。 得到GK, IPCM, GKIPCM对应的模糊隶属度分别为u_{ik, GK}, u_{ik, IPCM}, u_{ik, GKIPCM}, 表示第k个测试样本对第i类的模糊隶属度; 同时, IPCM与GKIPCM还具有典型值t_{ik, IPCM}, t_{ik, GKIPCM}, 表示第k个测试样本属于第i类的可能性, 其中GKIPCM算法的典型值分布如图5所示。 典型值取消了测试样本x_k对各类的隶属度之和为1的约束条件, 减弱了聚类过程中噪声对分类的影响。

图5

Fig.5

	Figure Option ViewDownloadNew Window
	图5 GKIPCM典型值Fig.5 Typical value of GKIPCM

对于模糊隶属度u_ik与典型值t_ik, 若u_ik为u_k中的最大值或t_ik为t_k中的最大值则认为第k个测试样本可以被划分到第i个类别当中。 最终根据模糊隶属度与典型值得到GK, IPCM, GKIPCM的分类精度分别为: 63.89%, 91.67%, 97.22%。

2.3.1 权重指数m

在模糊聚类算法中, 权重指数m对算法的准确程度有着重要影响。 对于GKIPCM, 当m≥ 2时, 模糊隶属度矩阵在迭代过程中有着较好的收敛性; 但由于同一样本对不同类的典型值没有和为1这一约束条件, m过大时式(5)的计算结果会受到干扰, 甚至使下一次迭代时式(1)计算的模糊散射矩阵S_{fi, r}接近奇异值, 影响分类准确度。 此处固定主成分个数为23, 研究权重指数m在区间[2, 6]上时, 各算法的聚类准确率, 结果如表1所示。

表1

Table 1

表1(Table 1)

表1 不同参数m的GK, GKIPCM, IPCM算法聚类准确率 Table 1 Clustering accuracies of GK, GKIPCM, IPCM with different values of m m GK聚类 GKIPCM聚类 IPCM聚类 2.0 63.89% 97.22% 91.67% 2.4 90.28% 97.22% 91.67% 2.8 94.44% 95.83% 91.67% 3.2 94.44% 95.83% 94.44% 3.6 94.44% 94.44% 94.44% 4.0 94.44% 95.83% 94.44% 4.4 94.44% 95.83% 95.83% 4.8 94.44% 97.22% 87.50% 5.2 95.83% 97.22% 86.11% 5.6 95.83% 97.22% 87.50% 6.0 97.22% 93.06% 90.28%

表1 不同参数m的GK, GKIPCM, IPCM算法聚类准确率 Table 1 Clustering accuracies of GK, GKIPCM, IPCM with different values of m

由表1可见, m在有限范围内, GKIPCM聚类准确率高于GK聚类与IPCM聚类。

2.3.2 主成分个数

在PCA降维过程中, 通常使得主成分的累积贡献率达到一个较高的百分比(80%~90%), 否则样本数据维度降低的同时也会损失大量的特征鉴别信息。 此处固定m=2, 通过调整主成分的个数从而调整累计贡献率, 观察对聚类准确率的影响, 结果如表2所示。

表2

Table 2

表2(Table 2)

表2 不同主成分的个数下GK, GKIPCM, IPCM算法的准确率 Table 2 Clustering accuracies of GK, GKIPCM, IPCM with different numbers of principal components 主成分个数 GK聚类 GKIPCM聚类 IPCM聚类 11 56.94% 51.39% 59.72% 13 56.94% 87.50% 54.17% 15 88.89% 87.50% 83.33% 17 86.11% 87.50% 86.11% 19 86.11% 94.44% 88.89% 21 47.22% 95.83% 91.67% 23 63.89% 97.22% 91.67% 25 88.89% 94.44% 91.67% 27 93.06% 94.44% 94.44% 29 86.11% 95.83% 95.83% 31 94.44% 95.83% 95.83% 33 93.06% 95.83% 91.67% 35 81.94% 93.06% 91.67% 37 93.06% 94.44% 90.28% 39 94.44% 93.06% 90.28%

表2 不同主成分的个数下GK, GKIPCM, IPCM算法的准确率 Table 2 Clustering accuracies of GK, GKIPCM, IPCM with different numbers of principal components

由表2知: 随着主成分个数的增加, GK, GKIPCM, IPCM的准确率都呈现上升趋势, GK的准确率相较其他算法对参数敏感一些, GKIPCM聚类的准确率在主成分个数为23时高达到97.22%。 综合表1、 表2, 该GKIPCM算法适合处理农药残留程度的定性分析问题。