一种基于太赫兹光谱的农产品中农药残留定量检测方法与流程

1.本发明属于食品安全检测技术领域，具体涉及一种基于太赫兹光谱的农产品中农药残留定量检测方法。


背景技术：

2.近些年，我国农产品生产规模正在逐年扩大，产生了巨大的经济效益和社会效益，但农产品生产面临着很多问题，其中最重要的就是农产品中农药残留超标问题屡禁不止。
3.与传统辐射源相比，太赫兹辐射具有强穿透性、瞬态性、宽带性、相干性、安全性、低能性等诸多特性，太赫兹光谱技术在农产品检测领域展现出巨大的应用潜力，结合相关机器学习算法，可实现农产品中农药残留含量的准确鉴定。目前，已有众多学者开展了基于太赫兹光谱的农产品中农药残留的定量检测研究。但是，目前研究主要具有如下两个问题：
4.(1)、基于太赫兹光谱的农产品中多种农药残留定量检测中，由于残留成分复杂，吸收谱线重叠混乱，导致光谱数据与农药残留浓度线性不可分，针对此特点的数据无法获得有效的定量回归模型，实现基于太赫兹光谱的农药残留含量的准确检测；
5.(2)、已报道的基于太赫兹光谱的农产品中多种农药残留定量检测中，农产品中农药残留检测浓度偏高，缺乏实际应用价值。
6.上述问题限制了太赫兹光谱在农药残留检测中的研究与应用。


技术实现要素：

7.本发明所要解决的技术问题是针对上述现有技术的不足提供一种基于太赫兹光谱的农产品中农药残留定量检测方法，基于太赫兹光谱的农产品中农药残留定量检测方法可用于实现农药残留含量快速、准确、无损检测，具有广阔的应用潜力。
8.为实现上述技术目的，本发明采取的技术方案为：
9.一种基于太赫兹光谱的农产品中农药残留定量检测方法，包括：
10.步骤1：将农产品与已知品种的农药混合，且获取混合后的农产品中农药的残留浓度实际值，制作多组不同残留浓度的样本，采用太赫兹光谱仪对所有样本进行检测，得到太赫兹吸收光谱数据集；
11.步骤2：对步骤1中检测得到的吸收光谱数据进行光谱预处理；
12.步骤3：基于预处理后的吸收光谱数据，进行bp神经网络非线性回归模型的构建，其中模型的输入为步骤2预处理后的吸收光谱数据，模型的输出为农药的残留浓度，构建完成后，得到bpnn模型；
13.步骤4：对步骤3构建的bpnn模型进行参数优化；
14.步骤5：采用pso算法从步骤2预处理后的吸收光谱数据中选取多个频段，将每个频段的吸收光谱数据分别输入到参数优化后的bpnn模型，预测得到每个频段对应的农药残留浓度，选择其中预测精度最高的频段，该频段记为最佳频段；
15.步骤6：选取与步骤1的农产品中具有相同农药品种的待测农产品，采用太赫兹光
谱仪对该待测农产品进行检测，得到太赫兹吸收光谱数据，对该太赫兹吸收光谱数据进行光谱预处理，得到预处理后的吸收光谱数据，从预处理后吸收光谱数据中选取与步骤5中最佳频段相同的频段区间，将该频段区间的吸收光谱数据输入到参数优化后的bpnn模型中进行检测，得到农产品的农药残留浓度检测结果。
16.作为本发明进一步改进的技术方案，所述的光谱预处理具体包括：
17.依次采用光谱去噪法和基线校正法对太赫兹吸收光谱数据进行光谱预处理，所述光谱去噪法为小波软阈值去噪法，所述基线校正法为二阶导基线校正法。
18.作为本发明进一步改进的技术方案，所述步骤3中bpnn模型的构建如下：
19.p＝wi；
20.即:
[0021][0022]
其中，w表示神经网络权重矩阵，i表示光谱输入矩阵，p表示预测浓度矩阵；w
mv
表示第m层网络第v个神经元的权重系数，inputn表示第n个样本的太赫兹吸收光谱数据。
[0023]
作为本发明进一步改进的技术方案，所述步骤4具体为：
[0024]
通过误差反向传递法进行参数优化，其中误差表达式为：
[0025]ek
＝y
k-ok；
[0026]
其中，yk表示样本的检测浓度，ok表示样本的预测浓度，ek表示样本的预测误差。
[0027]
作为本发明进一步改进的技术方案，所述步骤1中的农产品含有的农药品种有1种或多种，所述步骤1中某一农产品包含的农药品种与步骤6中待测农产品包含的农药品种相同。
[0028]
作为本发明进一步改进的技术方案，所述农产品中含有的已知品种的农药包含6-苄氨基嘌呤、吡虫啉和6-二氯苯腈。
[0029]
作为本发明进一步改进的技术方案，所述的步骤5具体为：
[0030]
在0-0.50thz峰宽范围内，以0.02thz为峰宽递增间隔，并选用pso算法计算不同峰宽下频段起始位置，继而确定不同峰宽对应的多个频段，将每个频段的吸收光谱数据分别输入到参数优化后的bpnn模型，预测得到每个频段的吸收光谱数据对应的农药残留浓度，选择其中预测精度最高的频段，该频段记为最佳频段。
[0031]
作为本发明进一步改进的技术方案，若样本中农产品含有的农药品种有1种，计算每个频段的吸收光谱数据对应的农药残留浓度预测值的相关系数或者均方根误差，选择其中相关系数最大或者均方根误差最小的频段，该频段为最佳频段；
[0032]
若样本中农产品含有的农药品种有多种，计算单个频段的吸收光谱数据对应的所
有农药残留浓度预测值的相关系数的平均值或者均方根误差的平均值，选择其中相关系数的平均值最大或者均方根误差的平均值最小的频段，该频段为最佳频段。
[0033]
本发明的有益效果为：
[0034]
(1)、本发明实现了多种农药残留成分的低浓度定量检测任务，通过光谱去噪、基线校正，解决由于农药残留成分复杂，造成的吸收谱线重叠混乱、光谱质量不佳的问题。
[0035]
(2)、本发明通过bpnn的构建，解决光谱数据与农药残留浓度线性不可分的问题，且通过网络的调参，进一步提升网络映射的完整性，提高农药残留含量定量检测的精度。
[0036]
(3)、本发明通过太赫兹光谱的波长选取，遴选出高质量的最佳频段区间，选取对应区间内的吸收光谱数据以定量分析，可获得更理想的检测效果。
附图说明
[0037]
图1为农产品小麦粉中三种农药(6-ba、imidacloprid和2，6-d)的制样浓度示意图。
[0038]
图2为使用太赫兹光谱仪进行样本检测，三种农药的thz光谱检测结果示意图。
[0039]
图3为不同浓度农药残留混合物的太赫兹吸收系数谱图。
[0040]
图4为光谱小波去噪后的农药残留混合物的太赫兹吸收系数谱图。
[0041]
图5为光谱小波去噪及基线校正后的农药残留混合物的太赫兹吸收系数谱图。
[0042]
图6中(a)为使用mw-bpnn得到0.35thz峰宽下的理想吸收频段及对应检测结果示意图。
[0043]
图6中(b)为使用mw-bpnn得到0.30thz峰宽下的理想吸收频段及对应检测结果示意图。
[0044]
图6中(c)为使用mw-bpnn得到0.25thz峰宽下的理想吸收频段及对应检测结果示意图。
[0045]
图6中(d)为使用mw-bpnn得到0.20thz峰宽下的理想吸收频段及对应检测结果示意图。
[0046]
图7中(a)为使用pso对混合物光谱进行波长选取，并给出了不同峰宽对应的6-ba定量检测结果示意图。
[0047]
图7中(b)为使用pso对混合物光谱进行波长选取，并给出了不同峰宽的imidacloprid定量检测结果示意图。
[0048]
图7中(c)为使用pso对混合物光谱进行波长选取，并给出了不同峰宽的2,6-d定量检测结果示意图。
[0049]
图7中(d)为使用pso对混合物光谱进行波长选取，并给出了不同峰宽的三种农药平均定量检测结果示意图。
具体实施方式
[0050]
下面根据附图对本发明的具体实施方式作出进一步说明：
[0051]
本实施例中提供一种基于太赫兹光谱的农产品中农药残留定量检测方法，具体提供了一种基于太赫兹光谱的农产品中三种农药残留物，6-苄氨基嘌呤(6-ba)、吡虫啉(imidacloprid)和2，6-二氯苯腈(2,6-d)的定量检测方法，图1为农产品小麦粉中三种农药
(6-ba、imidacloprid和2，6-d)的制样浓度，共28组不同农药残留含量浓度样本，农药浓度间隔为0.5mg，残留含量共计3mg。检测主要包括以下步骤：
[0052]
步骤(1)：采用太赫兹光谱仪进行样本检测，构建农产品与农药的混合物光谱数据集：
[0053]
测试得到三种农药在0.5-2.5thz频段内的吸收系数谱，如图2所示。6-苄氨基嘌呤在2.08thz处有一明显吸收峰；吡虫啉存在五个吸收峰，分别位于1.20,1.45,1.89,2.05及2.41thz；敌草腈存在3个吸收峰，分别位于1.58,1.98及2.32thz处。
[0054]
在此基础上，测试得到了小麦粉中不同农药残留含量混合物光谱，如图3所示，包含28组样本的吸收光谱数据，由于三种农药均在2.0thz左右处存在一个吸收峰，其中，6-苄氨基嘌呤在2.08thz处有吸收峰，吡虫啉在2.05thz处有吸收峰，敌草腈在1.98thz处有吸收峰，图中可以清楚地观察到混合物光谱在2.0thz处展现了一个明显的吸收峰，其他位置由于样品浓度过低，且浓度间隔较少，无法观察到明显的吸收峰，且谱线出现交叉重叠，不易区分，在后续定量分析中，光谱交叉重叠会严重影响定量分析的精度。
[0055]
步骤(2)：根据光谱特征进行相关光谱预处理(光谱去噪、基线校正)：
[0056]
将步骤(1)中检测得到的吸收光谱数据进行相关光谱预处理，改善光谱信号质量，光谱去噪选取小波软阈值去噪法：
[0057]
f(t)＝s(t) e(t)；
[0058]
wtf(y,x)＝wts(y,x) wte(y,x)；
[0059]
其中，f(t)表示样本信号；s(t)表示有效信号，e(t)表示噪声信号；x，y分别为原始函数中自变量及因变量，对应实例中农药残留浓度及太赫兹光谱吸收系数；wtf(y,x)、wts(y,x)和wte(y,x)分别表示样本信号、有效信号和噪声信号的小波变换函数，并选取自适应阈值进行软阈值去噪。
[0060]
本实例中具体表现为：对输入的含噪光谱信号进行多层小波分解，并采用小波软阈值对各层小波的高频系数进行收缩处理，最后对处理后的小波系数进行小波反变换，得到去噪后的光谱信号。如图4所示为小麦粉中不同农药残留含量混合物光谱去噪后数据，可以看出去噪后光谱质量得到明显改善，特别地，在吸收峰2.0thz处附近，交叉重叠现象较原始光谱有所优化。
[0061]
基线校正选取二阶导基线校正法：
[0062][0063][0064]
其中，y
k,t1
表示x＝k处光谱信号的一阶导数，y
k
△
，yk‑△
表示与k距离为
△
的信号强度；y
k,t2
表示x＝k处光谱信号的二阶导数，y
k
△
，yk‑△
表示与k距离为
△
的信号强度，yk表示k处信号强度。
[0065]
本实例中具体表现为：在频段范围内对各离散光谱点进行二阶导数求解，将二阶导数值为零的点进行记录，并以此点为中心实现基线校正，得到基线校正后的光谱信号。如图5所示为小麦粉中不同农药残留含量混合物光谱去噪并基线校正后数据，与未校正前相比，在吸收峰处，光谱信号质量得到较大改善，不同浓度光谱在吸收峰2.0thz附近区分度更加明显。
[0066]
步骤(3)：基于优化后的农产品与农药的混合物光谱数据，进行非线性回归模型bp神经网络(bpnn)的构建。
[0067]
bpnn的构建如下：
[0068]
p＝wi；
[0069]
即选取步骤(2)预处理后光谱数据为模型输入，农药残留浓度为输出，进行模型构建。
[0070][0071]
其中，w表示神经网络权重矩阵，i表示光谱输入矩阵，p表示预测浓度矩阵；w
mv
表示第m层网络第v个神经元的权重系数，inputn表示第n个样本的太赫兹吸收光谱数据。
[0072]
步骤(4)：对步骤(3)构建的bpnn模型进行参数优化。
[0073]
通过不断优化网络误差，反向传递误差实现参数优化：
[0074]ek
＝y
k-ok；
[0075]
其中，yk表示样本的检测浓度，ok表示样本的预测浓度，ek表示样本的预测误差。
[0076]
步骤(5)：结合移动窗口bpnn(mw-bpnn)及粒子群算法(pso)，对混合物光谱进行波长选取，遴选出高质量的最佳频段光谱区间。
[0077]
所述步骤5中mw-bpnn模型中信号窗口宽度的构建如下：
[0078][0079]
其中，input
mw
表示选定峰宽的光谱输入信号矩阵，x
n,k
表示第n个样本的第k个光谱数据点，p表示信号峰宽。
[0080]
本实例中具体表现为：在1.00-2.30thz频段内，分别以0.20thz、0.25thz、0.30thz及0.35thz为对应峰宽、以0.02thz为移动间隔(如以0.20thz为对应峰宽，依次取1.0-1.20thz、1.02-1.22thz、1.04-1.24thz
…
频段；如以0.25thz为对应峰宽，依次取1.0-1.25thz、1.02-1.27thz、1.04-1.29thz
…
频段；如以0.30thz为对应峰宽，依次取1.0-1.30thz、1.02-1.32thz、1.04-1.34thz
…
频段；如以0.35thz为对应峰宽，依次取1.0-1.35thz、1.02-1.37thz、1.04-1.39thz
…
频段)，选取光谱输入信号矩阵，建立mw-bpnn模型进行农药残留含量定量分析，即将选取的每个频段作为光谱输入信号矩阵输入至参数优化后的bpnn内，得到农药残留成分定量检测结果，选取每个峰宽对应的检测精度高的频段。如图6所示，即四种特定峰宽表征下，农产品中农药残留成分定量检测结果。其中图6中(a)、(b)、(c)、(d)中的灰色斜线填充部分分别为对应峰宽下的检测精度高的频段。结果表明，不
同峰宽表征下，定量检测精度最优区间均为包含吸收峰2.0thz的区间，表明了吸收峰在物质吸收光谱中的特殊性。其中r
p-avg
为灰色斜线填充部分频段对应的农药残留浓度预测值(检测值)的相关系数平均值，rmse
p-avg
为灰色斜线填充部分频段对应的农药残留浓度预测值(检测值)的均方根误差平均值。
[0081]
所述步骤5中pso的构建如下：适应度函数公式如下：
[0082][0083]
其中，f(a)表示适应度函数，dmmax表示输入样本与随机样本的最大聚集距离，dmmean表示平均聚集距离。
[0084]
粒子群中每个粒子具有两个属性：速度和位置，速度代表移动的快慢，位置代表移动的方向。每个粒子在搜索空间中单独的搜寻最优解，粒子的速度及位置更新如下：
[0085]
vi＝wvi c1×
rand()
×
(pbest
i-xi)；
[0086]
xi＝xi vi；
[0087]
其中，i＝1,2,3
…
n，i为当前粒子位置，n为粒子总数，vi为粒子运动记忆项，表示当前粒子运动速度大小和方向，pbesti表示设定的跟踪极值，表征当前粒子指向粒子自身最优点的一个矢量，xi表示粒子当前的位置，两者做差求解粒子下一步运动的速度及方向，其中，rand()为0至1之间的随机数，c1为学习因子,w为惯性权重；
[0088]
惯性权重影响粒子速度及搜索精度，为两者间达到平衡，定义动态惯性权重w如下：
[0089][0090]
其中，wmax为最大惯性权重，wmax＝0.8，wmin为最小惯性权重，wmin＝0.4，run为当前迭代次数，runmax为最大迭代次数；
[0091]
设置最大迭代次数300及优化目标相关系数r》0.97，满足其一，则pso运行结束，获取最佳惯性权重w，并确定频段选取起始位置，否则更新粒子速度及位置，循环计算。其中run每循环一次，更新一次最佳惯性权重w，w替代到上式粒子位置和速度公式中，更新一次粒子位置和粒子速度，对应一个频段选取的起始点(吸收峰宽是确定的)，进而确定频段区间。
[0092]
本实例中具体表现为：在0-0.50thz峰宽范围内，以0.02thz为峰宽递增间隔(即在0-0.50thz峰宽范围内，分别选取0.02thz、0.04thz、0.06thz、0.08thz、0.10thz
…
为峰宽，以上峰宽均以0.02thz为移动间隔)，并选用pso算法计算不同峰宽下多个频段起始位置，继而确定不同峰宽对应的多个频段区间，将每个频段区间内的吸收光谱数据分别输入到参数优化后的bpnn模型，预测得到每个频段区间内的吸收光谱数据对应的农药残留浓度，选择其中预测精度最高的频段，该频段记为最佳频段，从而实现不同吸收谱频段表征下农药残留含量检测精度的讨论。图7给出了不同峰宽对应的每种农药残留浓度预测值的相关系数以及三种农药残留浓度预测值的相关系数平均值，其中，图7中的(d)中圆点表示理想定量检测精度的对应峰宽，可以观察到，当吸收峰宽选取0.20-0.25thz区间时，定量检测精度较高，表明该峰宽可较完整地表征光谱，同时验证了pso算法在峰宽选取中的适用性。通过pso选取一个峰宽对应的多个频段，从多个频段对应的农药残留浓度预测值中选取其中最大相
关系数的值作为该峰宽对应的最优相关系数，然后再通过峰宽和对应的最优相关系数绘制图7中(a)、(b)、(c)，图7中(d)的纵坐标为图7中(a)、(b)、(c)的相关系数的平均。图7中纵坐标的相关系数对应的频段区间均在吸收峰附近。
[0093]
步骤(6)：选用参数优化后bpnn模型对预处理及波长选取后光谱进行检测，并对样本数据采用留一交叉验证法(loocv)进行划分，得到预测集的定量回归系数(r)及均方根误差(rmse)。
[0094]
相关系数r、均方根误差rmse的计算如下：
[0095][0096][0097]
其中，n表示检测样本数，yi表示第i个样本对应的农药残留浓度的实际值，表示所有样本对应的农药残留浓度的平均值，y
′i表示第i个样本对应的农药残留浓度的预测值。
[0098]
下面结合实验数据，对本案的发明效果给出说明。
[0099]
表1为仅进行参数优化后bpnn农药残留物的定量检测准确率。
[0100]
表1、参数优化后bpnn农药残留物的定量检测准确率：
[0101][0102][0103]
表2为进行参数优化及波长选取后农药残留物的定量检测准确率，同时给出了两种传统定量检测算法(偏最小二乘回归pls、支持向量机回归svr)的定量检测精度，并与本案进行对比。
[0104]
表2、参数优化及波长选取后农药残留物的定量检测准确率：
[0105][0106]
从表1和表2可以看出，对光谱进行相关预处理，可有效提升光谱质量；同时，对定量模型进行合理调参、光谱波长选取，可进一步提升检测精度。此外，与pls、svr对比，本案中bpnn具有较强的非线性拟合能力，优化后的网络可以更完整地映射太赫兹光谱与残留浓度间的对应关系。
[0107]
此外，对于未知农药残留浓度的待测农产品(训练集内最低已知残留浓度《农药残留浓度《训练集内最高已知残留浓度)，可参考实例实施以下步骤，实现农药残留含量检测：
[0108]
选取与步骤(1)的农产品中具有相同农药品种的待测农产品，采用太赫兹光谱仪对该待测农产品进行检测，得到太赫兹吸收光谱数据，对该太赫兹吸收光谱数据选用小波阈值去噪及直接差分二阶求导进行光谱预处理，得到预处理后的吸收光谱数据，从预处理后吸收光谱数据中选取与步骤(5)中最佳频段相同的频段区间，将该频段区间的吸收光谱数据输入到参数优化后的bpnn模型中进行检测，得到农产品的农药残留浓度检测结果。
[0109]
以上所述仅是本发明的优选实施方式，应当指出：对于本技术领域的普通技术人员来说，在不脱离本发明原理的前提下，还可以做出若干改进和润饰，这些改进和润饰也应视为本发明的保护范围。