CN114219262A - 一种区域农产品超标概率预测及误差分析方法 - Google Patents

Abstract

本申请涉及一种区域农产品超标概率预测及误差分析方法，利用研究区域已有土壤监测信息与农产品监测信息，将研究区域划分为不同子区域，并在子区域内进行迭代精确计算，增进农产品风险预测模型的准确性及适用性，准确率大幅提升的同时，可以大大降低现场采样工作量，为产地环境预警工作提供数据支撑。

Description

一种区域农产品超标概率预测及误差分析方法

技术领域

本发明属于农业环境技术领域，具体涉及一种区域农产品超标概率预测及误差分析方法。

背景技术

我国农产品污染逐年严重，分析主要成因、明确关键环节、提出相应对策对管控农产品重金属超标风险具有重要意义。农产品重金属超标存在风险不确定性高、影响因子复杂、积累生理机制不尽相同等特点，而目前风险预测方法较为单一且结果准确性较低。此外，农产品安全与土壤质量息息相关，但目前土壤及农产品监测工作系统性较差且数据利用性低，农产品与土壤重金属含量之间的响应关系研究较少。

CN113020232A在加密布点检测结果基础上，对目标区域进行初步单元划分，依据各单元污染特征进一步分区，选取对应修复治理措施，并对治理效果进行跟踪评估及动态调整；CN112883137A公开了一种农产品产地安全加密调查及安全预警方法，包括农产品超标概率计算、加密调查、分类预警和预警响应等多阶段多流程；上述技术中均涉及加密布点的方法，该方法并不涉及基于影响主因的检测数据计算变成，进而获得加密点位的点位代表面积，代表研究区域特征的效力较低。

上述技术方法存在的主要问题有：(1)土壤及农产品监测数据综合利用性低；(2)农产品风险预测工作开展较少且技术单一，准确性较低且适用性较差；(3)农产品风险预测工作量的且资金、时间消耗多。

发明内容

为了解决上述技术问题，本发明公开了一种区域农产品超标概率预测及误差分析方法，所述方法包括：

(1)主因筛选

(1.1)建立产地环境数据库

收集研究区域的环境信息、历史土壤监测信息和农产品监测信息，建立产地环境数据库；

(1.2)主因筛选

以农产品监测指标j作为研究指标，采用相对熵Ei方法对研究区域土壤监测指标进行筛选，相对熵Ei计算公式具体如下：

其中：xj为农产品研究指标j的监测值；p(xi)为土壤监测指标归一化后数值构成的矩阵；q(xj)为农产品研究指标归一化后数值构成的矩阵；

归一化计算公式具体如下：

xi＝(Xi-MINi)/(MAXi-MINi)

其中Xi为监测指标监测值，MAXi为监测指标最大监测值，MINi为监测指标最小监测值；

基于各土壤监测指标的相对熵，计算各指标相对熵权重Wi，将各指标相对熵权重Wi降序排列后，选取权重和占比为前Y的监测指标，为影响主因，计算公式具体如下：

所述Y的取值范围为0-100％；优选的，所述Y的取值范围为55-85％；

(2)插值误差分析

(2.1)克里金插值

选取适用克里金插值方法，对各项影响主因进行插值运算，获取各栅格内的指标预测标准误差数据；同时确定研究区域栅格数量；

(2.2)插值误差分析

(2.2.1)确定各项影响主因预测标准误差范围，以单项影响主因是否超过标准误差的S为条件，将研究区域分为单指标符合区域及单指标不符区域两类；各项影响主因的单指标不符区域进行空间合并，得到综合指标不符区域，并提取区域内的历史土壤监测点位；

所述S取值范围为0-100％；优选的，所述S取值范围为50-85％；

(2.2.2)基于提取的历史土壤点位影响主因监测数据，采用GS+、Arcgis软件计算各项指标的变程，并以最小变程为边长单位，对综合指标不符区域进行网格划分；若网格内不存在历史土壤监测点位则获取网格中心点，所述网格中心点即为加密点位，并提取加密点位经纬度信息；

(2.2.3)依据加密点位经纬度信息，现场获取各点位的所有土壤监测信息；现场获取的所有土壤监测信息与历史土壤监测信息内容一致；

(2.3)针对综合指标不符区域，结合区域内历史土壤监测点位及加密点位的各项影响主因监测数据，循环进行插值运算及误差分析，满足研究需求则进入下步；所述研究需求为所有综合指标不符区域面积均小于30*研究区域最大点位代表面积；

(3)农产品含量预测误差分析

基于研究区域所有栅格进行农产品含量预测误差分析

(3.1)误差计算

(3.1.1)栅格内存在历史农产品监测点位

分析研究区域历史农产品研究指标监测数据分布特征，选取回归分析方法，获得农产品研究指标j对应的预测模型，进而计算出各历史农产品监测点位的预测含量值，并采用相对误差计算各点位的预测误差(Dj)，计算公式具体如下：

其中，cj为农产品研究指标j的含量预测值，Xj为农产品研究指标j的监测值；

(3.1.2)栅格内无历史农产品监测点位

采用K-近邻法，选取距离目标栅格最近的K个历史农产品监测点位，获取各点位的农产品预测误差及与目标栅格间的距离Lj，并基于Lj获取各点位对应误差权重(Wj)，进而计算出目标栅格农产品含量预测误差(Cj)，具体计算公式如下：

其中，Lk为第k个点位与目标栅格的距离；

Cj＝W1*C1+W2*C2+W3*C3+....+Wk*Ck

(3.2)误差分析

(3.2.1)确定农产品研究指标j的含量预测误差范围，以是否超过含量预测误差的M为条件，将研究区域分为误差符合区域及误差不符区域两类，提取误差不符区域内的历史土壤监测点位；

所述M取值范围为0-100％；优选的，M取值范围为15％-50％；

(3.2.2)基于提取的历史土壤监测点位影响主因监测数据，进行变程计算并获取协同加密点位，计算流程同(2.2.2)步；

(3.2.3)依据协同加密点位信息，进行土壤、农产品样品的采集及检测，记录监测信息；

所述土壤、农产品监测信息与数据库内监测信息一致；

(3.3)针对误差不符区域，结合区域内历史土壤监测点位及协同加密点位的各项监测指标数据，循环进行主因筛选及误差分析，满足研究需求则进入第下一步，所述研究需求为误差不符区域内农产品预测误差均小于15％；

(4)潜在农产品超标风险区确定

依据农产品含量预测模型，计算各栅格内的农产品含量，依据国家标准GB2762中研究指标j的限量值，将农产品分为超标类、不超标类、潜在超标类三类，输出评价结果分布图。

进一步地，所述环境信息包括但不限于区域面积、周边环境、污染源信息等；历史土壤监测信息包括：土壤监测指标、监测日期、施肥量、投入品、点位代表面积等；农产品监测信息包括：农产品监测指标、监测日期、种植制度、农产品类型、年产量等；

进一步地，所述农产品研究指标j为农产品监测指标中的一种或多种；

进一步地，所述土壤监测指标包括Cd、Hg、As、Pb、Cr、pH、CEC、SOM、Cu、Zn、Ni、Mg、Fe、Mn、Ca、Si、Se、K、N、P中的一种或多种；

进一步地，所述农产品监测指标包括Cd、Hg、As、Pb、Cr中的一种或多种；

进一步地，所述克里金插值方法包括但不限于：普通克里金、简单克里金、泛克里金、经典贝叶斯克里金、指示器克里金、概率克里金；

进一步地，所述克里金插值采用Arcgis软件运算；

进一步地，所述回归方法包括但不限于：线性回归方法、岭回归方法、对率回归方法、多元回归方法、逐步回归方法；

进一步地，所述K设为3-5个；

进一步地，所述超标类为研究指标j的含量值大于1.1倍限量值，不超标类为研究指标j的含量值小于等于0.9倍限量值，潜在超标类为研究指标j的含量值介于0.9倍限量值与1.1倍限量值之间；

进一步地，所述插值方法的所有土壤影响主因栅格数量一致。

本发明的一种区域农产品超标概率预测及误差分析方法，具有以下优点：

1.本发明利用研究区域已有土壤监测信息与农产品监测信息，将研究区域划分为不同子区域，并在子区域内进行迭代精确计算，增进农产品风险预测模型的准确性及适用性，准确率大幅提升的同时，可以大大降低现场采样工作量，为产地环境预警工作提供数据支撑；

2.本发明可选用多种预测方法，且预测结果准确度较为相近，降低了对工作人员的专业素质要求；

3.本发明同时考虑土壤数据形态分布及空间分布，增加各栅格内农产品及土壤数据预测准确性，与正常监测工作相比，节约大量工作时间和资金投入；

4.本发明在原有监测数据的基础上，选取最适宜加密布点方法完成加密布点工作，结合农产品含量预测误差、潜在农产品超标风险区确定等步骤，可以快速确定潜在超标风险区域，在此基础上，为地区规划、生态修复提供指导，大大降低了监测费用的同时，提高了环境预警工作准确性。

附图说明

图1为一种区域农产品超标概率预测及误差分析方法技术路线图；

图2为T县历史监测点位分布图；

图3为土壤监测指标相对熵结果图；

图4为变程统计分布图；

图5加密点位分布图；

图6插值误差分析结果图；

图7水稻Cd含量误差不符区域分布图；

图8水稻Cd含量风险分布图；

具体实施方式

下面通过实施例对本发明做进一步的详细说明，以令本领域技术人员参照说明书文字能够据以实施。

应当理解，本文所使用的诸如“具有”，“包含”以及“包括”术语并不排除一个或多个其它元件或其组合的存在或添加。

实施例1

1.主因筛选

①产地环境数据库建立：选取T县内农田为研究对象(图2)，收集该区域环境信息、历史土壤及农产品(均为水稻)监测数据，其中土壤监测指标包括Cd、Hg、As、Pb、Cr、pH、CEC、SOM、P、K、N、Cu、Zn、Ni、Al、Si、Fe、Mn、Ca、Mg、Ce、La、Mo、Sr、Se、Co、Ti、Na共28项；水稻监测指标包括Cd、Hg、As、Pb、Cr共5项；明确点位代表面积范围为25-150亩/点；

②以水稻Cd作为研究指标，对土壤各项监测指标进行相对熵计算，权重从大到小排序后，选取相对熵信息含量占比前70％的指标为主因，包括Pb、Cd、SOM、N、Hg、Ni、Zn、pH、As、CEC、Cr、Cu、K共13项，见图3；

2.插值误差分析

①采用Arcgis软件，选取普通克里金插值方法，对选取出的13项影响主因分别进行插值，同时确定T县栅格共计67017个；获取各栅格影响主因的预测标准误差范围；

②以各项影响主因预测标准误差范围的75％值为限值，将预测标准误差大于限值的区域划定为单指标不符区域；将13项影响主因的不符区域进行空间叠合，则该区域为综合指标不符区域；

③提取综合指标不符区域历史土壤监测点位187个，利用GS⁺软件分别计算各影响主因的变程，其中Cd的变程最小，为2400米(图4)；

④以2400米为单位边长，采用Arcgis软件对综合指标不符区域进行网格划分，提取网格中心点131个为加密点位(图5)，获取各点位经纬度信息，并采用X-荧光方法现场快速获取各点位的所有土壤监测信息；

⑤结合综合指标不符区域内加密点位及历史土壤监测点位的影响主因数据，重复①-④步，直至T县所有综合指标不符区域均小于4500亩，重复三次误差分析后，获得最终土壤监测点位305个(图6)；

3.农产品含量预测误差分析

①依据所有历史监测点位的土壤影响主因数据及水稻Cd数据，采用SPSS软件进行线性回归分析，获取相关回归模型如下：

Cd_(水稻)＝0.847+0.007Pb+0.320Cd+0.007SOM-0.002N+0.160Hg+0.008Ni-

0.002Zn-0.112pH-0.004As+0.004CEC-0.001Cr-0.001Cu-0.001K

②基于预测模型计算获得历史水稻监测点位对应的Cd含量，并依据公式计算对应的水稻含量预测误差；

③选取最邻近的3个历史水稻点位的含量预测误差，分别计算各栅格内的水稻含量预测误差，获取水稻含量预测误差分布图；

④基于水稻含量预测误差分布图，获取区域内所有土壤监测点位的影响主因监测值，通过SPSS计算各指标的变程(最小值为1600米)，并依此得到协同加密点位45个，同步获取当季水稻及土壤的监测信息(过程同第2步)；

⑤依据上述获取的水稻监测数据及土壤监测数据，重复①-④步，直至T县所有误差不符区域内的水稻Cd含量预测误差均小于15％；重复一次后误差不符区域内的水稻Cd含量预测误差最大值为10.4％(图7)；

4.潜在农产品超标风险区确定

查询最新国家标准GB2762-2017(《食品安全国家标准食品中污染物限量》)，明确水稻Cd的限量值为0.2mg/kg，基于水稻含量预测模型计算出的各栅格的水稻Cd含量以0.18mg/kg和0.22mg/kg为限值，将水稻Cd含量分为三类，绘制分布图(图8)，其中潜在水稻超标风险区即为潜在超标类区域。

尽管本发明的实施方案已公开如上，但其并不仅仅限于说明书和实施方式中所列运用，它完全可以被适用于各种适合本发明的领域，对于熟悉本领域的人员而言，可容易地实现另外的修改，因此在不背离权利要求及等同范围所限定的一般概念下，本发明并不限于特定的细节和这里示出与描述的实施例。

Claims (10)

1.一种区域农产品超标概率预测及误差分析方法，其特征在于，所述方法包括：

所述方法包括：

(1)主因筛选

(1.1)建立产地环境数据库

(1.2)主因筛选

以农产品监测指标j作为研究指标，采用相对熵Ei方法对研究区域土壤监测指标进行筛选，相对熵Ei计算公式具体如下：

其中：xj为农产品研究指标j的监测值；p(xi)为土壤监测指标归一化后数值构成的矩阵；q(xj)为农产品研究指标归一化后数值构成的矩阵；

归一化计算公式具体如下：

xi＝(Xi-MINi)/(MAXi-MINi)

其中Xi为监测指标监测值，MAXi为监测指标最大监测值，MINi为监测指标最小监测值；

基于各土壤监测指标的相对熵，计算各指标相对熵权重Wi，将各指标相对熵权重Wi降序排列后，选取权重和占比为前Y的监测指标，为影响主因，计算公式具体如下：

所述Y的取值范围为0-100％；

(2)插值误差分析

(2.1)克里金插值

选取适用克里金插值方法，对各项影响主因进行插值运算，获取各栅格内的指标预测标准误差数据；同时确定研究区域栅格数量；

(2.2)插值误差分析

(2.2.1)确定各项影响主因预测标准误差范围，以单项影响主因是否超过标准误差的S为条件，将研究区域分为单指标符合区域及单指标不符区域两类；各项影响主因的单指标不符区域进行空间合并，得到综合指标不符区域，并提取区域内的历史土壤监测点位；

所述S取值范围为0-100％；

(2.2.2)基于提取的历史土壤点位影响主因监测数据，采用GS+、Arcgis软件计算各项指标的变程，并以最小变程为边长单位，对综合指标不符区域进行网格划分；若网格内不存在历史土壤监测点位则获取网格中心点，所述网格中心点即为加密点位，并提取加密点位经纬度信息；

(2.2.3)依据加密点位经纬度信息，现场获取各点位的所有土壤监测信息；现场获取的所有土壤监测信息与历史土壤监测信息内容一致；

(2.3)针对综合指标不符区域，结合区域内历史土壤监测点位及加密点位的各项影响主因监测数据，循环进行插值运算及误差分析，满足研究需求则进入下步；所述研究需求为所有综合指标不符区域面积均小于30*研究区域最大点位代表面积；

(3)农产品含量预测误差分析

基于研究区域所有栅格进行农产品含量预测误差分析

(3.1)误差计算

(3.1.1)栅格内存在历史农产品监测点位

分析研究区域历史农产品研究指标监测数据分布特征，选取回归分析方法，获得农产品研究指标j对应的预测模型，进而计算出各历史农产品监测点位的预测含量值，并采用相对误差计算各点位的预测误差(Dj)，计算公式具体如下：

其中，cj为农产品研究指标j的含量预测值，Xj为农产品研究指标j的监测值；

(3.1.2)栅格内无历史农产品监测点位

采用K-近邻法，选取距离目标栅格最近的K个历史农产品监测点位，获取各点位的农产品预测误差及与目标栅格间的距离Lj，并基于Lj获取各点位对应误差权重(Wj)，进而计算出目标栅格农产品含量预测误差(Cj)，具体计算公式如下：

其中，Lk为第k个点位与目标栅格的距离；

Cj＝W1*C1+W2*C2+W3*C3+....+Wk*Ck

(3.2)误差分析

(3.2.1)确定农产品研究指标j的含量预测误差范围，以是否超过含量预测误差的M为条件，将研究区域分为误差符合区域及误差不符区域两类，提取误差不符区域内的历史土壤监测点位；

所述M取值范围为0-100％；

(3.2.2)基于提取的历史土壤监测点位影响主因监测数据，进行变程计算并获取协同加密点位，计算流程同(2.2.2)步；

(3.2.3)依据协同加密点位信息，进行土壤、农产品样品的采集及检测，记录监测信息；

所述土壤、农产品监测信息与数据库内监测信息一致；

(3.3)针对误差不符区域，结合区域内历史土壤监测点位及协同加密点位的各项监测指标数据，循环进行主因筛选及误差分析，满足研究需求则进入第下一步，所述研究需求为误差不符区域内农产品预测误差均小于15％；

(4)潜在农产品超标风险区确定

依据农产品含量预测模型，计算各栅格内的农产品含量，依据国家标准GB2762中研究指标j的限量值，将农产品分为超标类、不超标类、潜在超标类三类，输出评价结果分布图。

2.如权利要求1所述的一种区域农产品超标概率预测及误差分析方法，其特征在于，所述产地环境数据库收集研究区域的环境信息、历史土壤监测信息和农产品监测信息。

3.如权利要求2所述的一种区域农产品超标概率预测及误差分析方法，其特征在于，所述环境信息包括但不限于区域面积、周边环境、污染源信息。

4.如权利要求2所述的一种区域农产品超标概率预测及误差分析方法，其特征在于，历史土壤监测信息包括：土壤监测指标、监测日期、施肥量、投入品、点位代表面积；农产品监测信息包括：农产品监测指标、监测日期、种植制度、农产品类型、年产量。

5.如权利要求4所述的一种区域农产品超标概率预测及误差分析方法，其特征在于，所述土壤监测指标包括Cd、Hg、As、Pb、Cr、pH、CEC、SOM、Cu、Zn、Ni、Mg、Fe、Mn、Ca、Si、Se、K、N、P中的一种或多种。

6.如权利要求4所述的一种区域农产品超标概率预测及误差分析方法，其特征在于，所述农产品监测指标包括Cd、Hg、As、Pb、Cr中的一种或多种。

7.如权利要求1所述的一种区域农产品超标概率预测及误差分析方法，其特征在于，所述克里金插值方法包括但不限于：普通克里金、简单克里金、泛克里金、经典贝叶斯克里金、指示器克里金、概率克里金。

8.如权利要求1所述的一种区域农产品超标概率预测及误差分析方法，其特征在于，所述回归方法包括但不限于：线性回归方法、岭回归方法、对率回归方法、多元回归方法、逐步回归方法。

9.如权利要求1所述的一种区域农产品超标概率预测及误差分析方法，其特征在于，所述K设为3-5个。

10.如权利要求1所述的一种区域农产品超标概率预测及误差分析方法，其特征在于，所述超标类为研究指标j的含量值大于1.1倍限量值，不超标类为研究指标j的含量值小于等于0.9倍限量值，潜在超标类为研究指标j的含量值介于0.9倍限量值与1.1倍限量值之间。

Images