CN109784769B - 一种农业面源污染风险识别方法 - Google Patents

Abstract

本发明公开了一种农业面源污染风险识别方法。它包括如下步骤：1)测定对象地区的风险识别指标体系中各指标的数值；2)计算风险识别指标体系中的各指标数值的权重，得到权重，记为W；对风险识别指标体系中各指标进行赋值，得到指标赋值，记为I；3)将权重W和指标赋值I代入下式(1)中，得到农业面源污染风险指数；NPSPRI＝∑W_i×I_i；上式中，NPSPRI为农业面源污染风险指数；W_i为各风险指标在指标体系中的权重值，其值范围为0～1；I_i为指标赋值，i为指标体系中的指标；根据计算得到的农业面源污染风险指数数值，即能识别农业面源污染风险。本发明能够简单、快速、低成本、适应性强、不依赖GIS和遥感技术等高精尖设备和技术，且能非常精确地识别农业面源污染风险。

Description

一种农业面源污染风险识别方法

技术领域

本发明涉及一种农业面源污染风险识别方法，属于农业污染潜在风险预警领域。

背景技术

目前，农业面源污染面广量大，污染源分散且隐蔽，污染发生的时间和空间具有随机性和不确定性，监测、控制难度大。识别农业面源污染的关键源区和高风险区，将有限的资源投入到对水体危害可能性最大而范围相对较小的区域进行重点治理，可大大降低治理难度和提高治理成效。因此建立一种农业面源污染风险识别方法是农业面源污染管理和控制的当务之急。

目前识别农业面源污染风险的方法有：输出系数法、面源污染定量模型法、指标体系法。输出系数法结构简单，所需资料较少，可直接评估和预测农业面源总氮和总磷的污染负荷量，但其在区域尺度上的应用需要大量的实地监测资料。面源污染定量模型包括Sparrow模型、GWLF模型、SWAT模型和HSPF模型等，需要参数较多，而目前农业管理中的数据积累还不够丰富，下垫面情况更复杂，区域性差异大，更是增加了地面基础信息的获取难度，因此国外模型直接移植的难度较大。指标体系法，可综合分析影响农业面源污染物流失的主要因子，能够为农业面源污染风险提供一个更为合理的评价框架，灵活性较强。常用的方法有非点源污染潜力指数法(agricultural pollution potential index，APPI)和磷指数法(phosphorus index，PI)。但常规农业面源污染风险指标体系法中存在考虑的污染来源分类少、指标选择不全面、研究单元太粗等问题。

发明内容

本发明的目的是提供一种农业面源污染风险识别方法，本发明能够简单、快速、低成本、适应性强、不依赖GIS和遥感技术等高精尖设备和技术，且能非常精确地识别农业面源污染风险。

本发明提供的一种农业面源污染风险识别方法，包括如下步骤：1)测定对象地区的风险识别指标体系中各指标的数值；

所述风险识别指标体系包括水文气象指标、土壤地形植被指标以及经济指标；

2)计算所述风险识别指标体系中的各指标数值的权重，得到权重，记为W；

对所述风险识别指标体系中各指标进行赋值，得到指标赋值，记为I；

3)将所述权重W和所述指标赋值I代入下式(1)中，得到农业面源污染风险指数；

NPSPRI＝∑W_i×I_i (1)

式(1)中，NPSPRI为农业面源污染风险指数；W_i为各风险指标在指标体系中的权重值，其值范围为0～1；I_i为指标赋值，i为指标体系中的指标；

根据计算得到的所述农业面源污染风险指数数值，即能识别农业面源污染风险。

上述的方法中，所述水文气象指标包括年降水量、溶解态面源污染物入河系数以及吸附态面源污染物入河系数；

所述土壤地形植被指标包括年植被覆盖度、坡度以及土壤可侵蚀性因子；

所述经济指标包括农田氮平衡量以及农田磷平衡量。

上述的方法中，所述年降水量根据测定所述对象地区各监测站的月降水量数据计算得出；

根据所述年降水量和测定所述对象地区各监测站的月径流量数据，通过下式(2)计算所述溶解态面源污染物入河系数：

式(2)中，CR为溶解态面源污染物入河系数，Prec和Runoff分别为年降水量和年径流量；

根据测定的泥沙含量数据和通过式(4)～(11)计算得到的土壤侵蚀量，通过下式(3)计算所述吸附态面源污染物入河系数；

式(3)中，SDR为吸附态面源污染物入河系数，Sed和Sel分别为年泥沙含量和年土壤侵蚀量；

Sel＝R×K×L×S×C×Ρ (4)

式(4)中，R为降雨侵蚀力因子，K为土壤可蚀性因子，L、S分别为坡长因子和坡度因子，无量纲；C为生物措施因子，无量纲；P为工程措施措施因子，无量纲。

式(5)中，P_i为月降雨量，P_year为年降雨量，i是月份；

根据(6)计算得到土壤可侵蚀性因子；

K＝0.1317×K_China (6)

K_China＝-0.01383+0.51575K_EPIC (7)

式(6)-(8)中，K_EPIC为EPIC模型计算得到的土壤可蚀性因子，K_china是中国土壤可蚀性因子，K为土壤可蚀性因子，0.1317为美国制和国际制的单位转换系数，S_a为砂粒含量，S_i为粉粒含量，C_i为粘粒含量；C是土壤有机碳含量，S_n＝1-S_a/100；

式(9)中，θ是坡度；

式(10)中，λ是坡长；

式(11)中，FVC为年植被覆盖度。

本发明中，对于土壤侵蚀计算中工程措施因子P由土地利用类型决定，具体取值为：当土地利用类型为水田时，P取值为0.15；当土地利用类型为旱地时，P取值为0.5；当土地利用类型为林地时，P取值为1；当土地利用类型为草地时，P取值为0.2；当土地利用类型为灌木时，P取值为0.2；当土地利用类型为水体时，P取值为0；当土地利用类型为滩地时，P取值为0.2；当土地利用类型为建设用地时，P取值为1.0。

本发明中，所述年植被覆盖度数据可从相关的地面数据平台购买，具体如地理国情监测云平台。

本发明中，所述坡度数据可通过相关的地面数据平台购买，具体如地理国情监测云平台。

上述的方法中，根据国家县级农业统计指标数据采用土壤系统收支平衡方法通过式(12)～(15)计算所述农田氮平衡量和所述农田磷平衡量；

Q_bal＝Balance/area×1000 (12)

其中，Balance＝Input-Output (13)

Input＝Ftlz+Mnr+Irg+Seed+Dpzt+Bnf (14)

Output＝Hvst+Gas+water (15)

其中，Q_bal为所述农田氮平衡量或所述农田磷平衡量；area为耕地面积和园地面积之和；1000为单位转换系数；Balance为养分平衡量；Input为养分输入量；Output为养分输出量；Ftlz为化肥养分输入量；Mnr为有机肥养分输入量；Irg为灌溉养分输入量；Seed为种子养分输入量；Bnf为生物固氮氮输入量；Dpzt为干湿沉降养分输入量；Hvst为作物带走养分输出量；Gas为氨气挥发和反硝化氮输出量；Water为淋失和径流养分输出量。

上述的方法中，当所述年降水量的取值分别为≤400mm、400～500mm、500～700mm、>700mm时，其赋值分别为1、2、3、4；

当所述溶解态面源污染物入河系数分别为≤0.018、0.018～0.055、0.055～0.130、>0.130时，其赋值分别为1、2、3、4；

当所述吸附态面源污染物入河系数分别为≤0.018、0.018～0.13、0.13～0.28、>0.28时，其赋值分别为1、2、3、4；

当所述年植被覆盖度的取值分别为>60％、45～60％、30～45％、≤30％时，其赋值分别为1、2、3、4；

当坡度分别为≤8、8～15°、15～25°、>25°时，其赋值分别为1、2、3、4；

当土壤可侵蚀性因子分别为≤0.010t·h·(MJ·mm)^-1、0.010～0.023t·h·(MJ·mm)^-1、0.023～0.027t·h·(MJ·mm)^-1、>0.027t·h·(MJ·mm)^-1时，其赋值分别为1、2、3、4；

当农田氮平衡量分别为≤0t·km^-2、0～15t·km^-2、15～40t·km^-2、>406t·km^-2时，其赋值分别为1、2、3、4；

当农田磷平衡量分别为≤0t·km^-2、0～5t·km^-2、5～15t·km^-2、>15t·km^-2时，其赋值分别为1、2、3、4。

本发明中，所述风险识别指标体系中各指标进行赋值不依赖于GIS技术和遥感技术。

本发明中，对所述农业面源污染风险指数NPSPRI进行自然间断点分级法(Jenks)分级为四个风险等级，详细分级标准为：无风险(0，2.208]、低风险(2.208，2.704]、中风险(2.704，3]和高风险(3，4)，即得到农业面源污染四级风险识别图(图9)；并将各个级别数值范围赋予1～4分，即：无风险为1分，低风险为2分，中风险为3分，高风险为4分。

本发明具有以下优点：

本发明方法能够简单、快速、低成本、适应性强、不依赖GIS和遥感技术等高精尖设备和技术且非常精确对农业面源污染风险进行识别。通过采用本发明的农业面源污染风险识别方法，能够简单、快速、低成本、适应性强、不依赖GIS和遥感技术等高精尖设备和技术且非常精确对农业面源污染风险进行识别。

所谓“简单”、“快速”，是指本发明的风险识别方法中采用了大量易于获得的指标。例如，所采用的指标中存在许多地方监测站常年监测的指标(如年降水量、坡度))，可利用水文气象监测站监测数据、统计数据、中国土壤数据库以及文献数据等非常容易计算的指标(如土壤侵蚀量、溶解态面源污染物入河系数、吸附态面源污染物入河系数、土壤氮平衡量以及土壤磷平衡量)。由于这些指标的存在，极大地减少了风险识别的难度和工作量，使得风险识别工作变得非常简单，数值常年固定的指标可以直接使用，不必再行测定，而常年监测的指标和容易计算的指标也可以快速获得，因此，极大地加快了风险识别工作的速度。

所谓“低成本”，是指本发明的风险识别方法采用了大量的常见指标。很多指标例如：年植被覆盖度、坡度、土壤侵蚀量、土壤氮平衡量以及土壤磷平衡量等的计算中用到的坡长数据、植被归一化指数(NDVI)、化肥数据、畜禽数据、农作物数据、人口数据可以由根据本领域常规方法监测或测定而得到，也可以由各环境监测、统计部门的公开数据获得，即使是需要购买的数据，也由于是常见数据而非常便宜，使得风险识别工作的成本大大降低。

所谓“适应性强”，是指本发明的风险识别方法可以广泛地适用于北方地区各种农田，例如，不同的农作物类型：水稻、小麦、蔬菜、经济作物、果树等；不同的地形：平原或山地等；不同的经济程度：农业经济发达地区或农业经济欠发达地区。由于在本发明的风险识别方法所采用的指标体系中，广泛地考虑了以上不同的情形，使得本发明的风险识别方法具有非常强的适应性。

所谓“不依赖GIS和遥感技术等高精尖设备和技术”，这是本发明的风险识别方法的重要技术效果之一。由于依赖于GIS技术和遥感技术的数据并非随时随地能够获得，需要经过特殊的手续，而且存在无法及时获得的可能性，成为了导致风险识别工作的瓶颈，另外，GIS技术和遥感技术需要极强的专业知识，而往往存在农业面源污染风险的地区的地方环保监测/评价单位不具备这样的专业知识，导致即使有了数据也无法充分地利用的问题。而本发明的风险识别方法中所采用的指标均不依赖GIS技术和遥感技术，一来克服了上述数据不易获得的“瓶颈”，二来降低了风险识别工作的技术难度，使得本发明的指标体系适用性更为广泛，且更具推广性和易操作性。

所谓“非常精确”，这也是本发明的风险识别方法的一个重要技术效果。本发明的发明人采用本发明的风险识别方法进行了实地的风险识别，并将识别结果与遥感分布式DPeRS模型的识别结果进行了比对，结果发现，二者具有非常高的一致性，这说明本发明的风险识别方法实现了非常精确的风险识别。

另外，本发明中更为突出的特点在于，本发明的发明人针对本发明的风险识别方法特意设计了一套指标体系，该指标体系中的各项指标的数值可以简单、快速、低成本地获得，且一级指标和二级指标都不需要依赖GIS和遥感技术。再者，本发明的发明人通过对上述指标体系中的各指标赋予特定的权重从而形成一套互为补充、扬长避短的有机的指标体系，更加显著地实现了以上的技术效果，特别是，通过特定的指标与特定的权重的组合，实现了极为精确的风险识别，能够实现与面源污染定量模型相媲美甚至更优的识别精度。

附图说明

图1是年降水量的分级图；

图2是溶解态面源污染物入河系数分级图；

图3是吸附态面源污染物入河系数分级图；

图4是年植被覆盖度分级图；

图5是坡度分级图；

图6是土壤可侵蚀性因子K分级图；

图7是农田氮平衡量分级图；

图8是农田磷平衡量分级图；

图9是海河流域农业面源污染综合指标风险等级图；

图10是海河流域DPeRS模型农业面源污染风险等级图；

图11是海河流域农业面源污染综合指标风险等级结果与DPeRS模型结果的偏差分析图。

具体实施方式

下述实施例中所使用的实验方法如无特殊说明，均为常规方法。

下述实施例中所用的材料、试剂等，如无特殊说明，均可从商业途径得到。

下面结合实施例，对本发明的具体实施方式作进一步详细描述。以下实施例用于说明本发明，但不用来限制本发明的范围。

实施例、一种农业面源污染风险识别方法

(一)它包括如下步骤：

步骤1：识别目标区域的基础数据的收集

在进行农业面源风险识别之前，需要对识别目标区域的基本情况进行摸底和核实。具体工作包括：

(1)收集以下自然及经济基础数据资料：风险识别目标区域的具体位置(经纬度等)、行政区划、地理边界、面积、地质地形地貌、水系分布、人口数量及分布、经济类型、支柱产业等；

(2)收集以下污染基础数据资料：风险识别目标区域的主要污染源、其强度及其分布、主要污染物质、污染源存在的时间；

(3)收集以下周边基础数据资料：风险识别目标区域的周边区域的污染源状况、周边区域污染对风险识别区域的影响等。

上述基础数据的收集方法包括：

实地调查、文献资料收集、历史档案资料收集、监测资料收集等。

步骤2：风险识别指标体系的确立；

在本发明的风险识别方法中，确立了以下的风险识别指标体系。

上述风险识别指标体系由以下一级指标构成：水文气象指标、土壤地形植被指标以及经济指标。

步骤3：确定上述风险识别指标体系中的各指标的权重，得到权重W；

在本发明中，为确定风险识别方法所采用的指标的权重，采用了层次分析法。步骤4：对上述风险识别指标体系中各指标进行赋值，得到指标赋值I；

由于各个指标的数值范围、单位等均千差万别，因此，计算风险指数之前，为了能够方便地对各个指标进行加和，对各个指标进行了分级赋值。在本发明中，将各个指标根据其不同的数值范围分为四个级别，并将这四个级别分别赋值1、2、3、4四个数值，从而完成各个指标的赋值I。

步骤5：将上述权重W和指标赋值I代入下式(1)中，获得农业面源污染风险指数NPSPRI；

NPSPRI＝∑W_i×I_i (1)

式(1)中，NPSPRI为农业面源污染风险指数；W_i为各风险指标在指标体系中的权重值(通过层次分析法获得)，其值范围为0-1；I_i为指标赋值，i为指标体系中的指标。

步骤6：对上述农业面源污染风险指数NPSPRI数值进行自然间断点分级法(Jenks)分级为四个风险等级，分别是：无风险、低风险、中风险以及高风险；并将各个级别数值范围赋予1～4分，即无风险为1分，低风险为2分，中风险为3分，高风险为4分。

上述自然间断点分级法(Jenks)是美国制图学家George F.Jenks提出的一种基于数据中固有自然分组的地图分级算法，即认为数据分布本身有断点，可利用数据这一特点分级；其算法原理是一个小聚类，聚类结束条件是组间方差最大、组内方差最小。

上述分级法的详细内容可以参考《ArcGIS地理信息系统空间分析实验教程》(第二版)(汤国安等著，科学出版社，2012)。

(二)本发明的风险识别方法中所采用的指标体系

(1)指标体系的内容

本发明涉及一种农业面源污染风险识别方法，其采用指标体系对农业面源污染风险进行识别，其中，所述指标体系由以下一级指标构成：水文气象指标、土壤地形植被指标以及经济指标。

本发明选择的指标的赋值均不依赖于GIS技术和遥感技术。需要说明的是，本发明中采用的部分指标(如坡度、年植被覆盖度)当然是可以通过GIS或遥感技术得到的，但是，这些指标也可以通过其他方法获得，这样的指标也被认为是不依赖于GIS技术和遥感技术的。另外，本发明的风险识别中所采用的相关地图当然是可以采用GIS或遥感技术获得的，但是这些地图并不属于本发明的指标体系的内容。

本发明通过层次分析法对上述水文气象指标、土壤地形植被指标以及经济指标的权重进行了确定，并具体赋值为：0.5396、0.1634以及0.2970。采用层次分析法确定权重时，通过层次分析软件进行一致性检验，最终确定权重。这样确定的权重全面考虑了各指标在风险识别中的影响，与特定的指标组合使用，能够更为精确地进行风险识别。

在上述一级指标之下，进一步设立二级指标。具体地，

水文气象指标由以下二级指标构成：年降水量、溶解态面源污染物入河系数以及吸附态面源污染物入河系数。

土壤地形植被指标由以下二级指标构成：年植被覆盖度、坡度以及土壤可侵蚀性因子K。

经济指标由以下二级指标构成：农田氮平衡量以及农田磷平衡量。

上述二级指标也通过层次分析法确定了指标权重，具体方法为：根据专家打分和文献查阅整理，对各指标相对于其他指标的重要程度进行排序打分，然后通过层次分析法，以通过一致性检验为标准，确定指标权重。

具体确定的指标权重为：

年降水量：0.2662；溶解态面源污染物入河系数：0.1677；吸附态面源污染物入河系数：0.1057；年植被覆盖度：0.0267；坡度：0.0485；土壤可侵蚀性因子K：0.0882；农田氮平衡量：0.1980；农田磷平衡量：0.0990。

在本发明中，通过上述特定的二级指标的组合和特定的权重，形成了一套有机的指标体系，该指标体系中的指标涵盖面广、统计全面，并且通过特定的权重使得各指标之间可以轻重有序、重点突出、相互补充、扬长避短。采用这样指标体系时，能够以简单的方法快速、精确识别出农业面源污染的高风险区。

上述指标体系以及权重的组合，是实现本发明的技术效果的最经济的指标体系，如果进一步增加指标，对于识别精度的提高有限，但是却可能付出非常高的经济成本、技术成本、人力成本或时间成本。

(2)二级指标数值的获取或计算方式

下面对以上三个一级指标所含的二级指标的获得方式或计算方式进行说明。

1.水文气象指标

1-1.年降水量

从目前已经公开的水文年鉴、气象共享资料中获得研究对象地区各监测站的月降水量数据，进而计算出年降水量。

1-2.溶解态面源污染物入河系数

根据1-1中的年降水量和从水文年鉴获得的月径流量数据通过下式(1)计算溶解态污染物入河系数：

其中，CR为溶解态污染物入河系数，Prec和Runoff分别为年降水量和年径流量。

1-3.吸附态面源污染物入河系数

根据《中国河流泥沙公报》以及水文监测站获得的泥沙含量数据和通过式(3)～(11)计算得到的土壤侵蚀量，再利用通过下式(2)计算吸附态污染物入河系数。

其中，SDR为吸附态污染物入河系数，Sed和Sel分别为年泥沙含量和年土壤侵蚀量。

Sel＝R×K×L×S×C×Ρ (4)

其中，R为降雨侵蚀力因子，K为土壤可蚀性因子，L、S分别为坡长因子和坡度因子(无量纲)；C为生物措施因子(无量纲)；P为工程措施措施因子(无量纲)。

其中，P_i为月降雨量，P_year为年降雨量，i是月份。

K＝0.1317×K_China (6)

K_China＝-0.01383+0.51575K_EPIC (7)

其中，K_EPIC为EPIC模型计算得到的土壤可蚀性因子，K_china是中国土壤可蚀性因子，K为土壤侵蚀方程中的土壤可蚀性因子，0.1317为美国制和国际制的单位转换系数，S_a为砂粒含量，S_i为粉粒含量，C_i为粘粒含量；C是土壤有机碳含量，S_n＝1-S_a/100。

其中，θ是坡度。

其中，λ是坡长。

其中，FVC为植被覆盖度。

对于土壤侵蚀计算中工程措施因子P按照现有研究成果，由土壤类型决定。例如，在本申请涉及的农业用地中，可以具体取值为：水田为0.15，旱地为0.5，林地为1，草地为0.2，灌木为0.2，水体为0，滩地为0.2，建设用地为1.0。

2.土壤地形植被指标

2-1.年植被覆盖度

年植被覆盖度数据可从相关的地面数据平台购买(如地理国情监测云平台等)。

2-2.坡度

坡度数据可通过相关的地面数据平台购买(如地理国情监测云平台等)。

2-3.土壤可侵蚀性因子K

土壤可侵蚀性因子K的计算见式(6)～(8)。

3.经济指标

依据国家县级农业统计数据指标采用土壤系统收支平衡方法核算农田氮磷平衡量空间数据(式(12)～(15))，并构建农田氮磷平衡量空间数据库。其农田氮磷平衡量空间数据具体算法为：

Q_bal＝Balance/area×1000 (12)

其中，Balance＝Input-Output (13)

Input＝Ftlz+Mnr+Irg+Seed+Dpzt+Bnf (14)

Output＝Hvst+Gas+water (15)

其中，Q_bal为农田氮磷平衡量空间数据；area为耕地面积和园地面积之和；1000为单位转换系数；Balance为养分平衡量；Input为养分输入量；Output为养分输出量；Ftlz为化肥养分输入量；Mnr为有机肥养分输入量；Irg为灌溉养分输入量；Seed为种子养分输入量；Bnf为生物固氮氮输入量；Dpzt为干湿沉降养分输入量；Hvst为作物带走养分输出量；Gas为氨气挥发和反硝化氮输出量；Water为淋失和径流养分输出量。

本发明的农业面源污染风险识别指标体系可以适用于各种农田，例如，不同的农田类型：旱田、水田或梯田；不同的农作物类型：水稻、小麦等粮食作物、蔬菜、经济作物、果树等；不同的地形：平原或山地等；不同的经济程度：农业经济发达地区或农业经济欠发达地区。具体而言，上述指标中的年降水量以及污染物入河系数等指标即为考虑了以水分划分的农田类型；坡度等指标即为考虑了平原或山地等不同地形；降水量、农田氮平衡量、农田磷平衡量等指标即为考虑不同的农作物类型以及经济程度。

另外，基于本发明的发明人对以上不同情况的农田进行风险识别的结果来看，本发明的风险识别指标体系更加适用于我国北方农田。本发明的风险识别指标体系对于我国北方农田的面源污染风险的识别准确度完全可以与模型法得到的结果相媲美。特别地，本发明的风险识别指标体系更优选适用于我国北方、特别是东北、华北、西北地区的粮食作物、蔬菜、经济作物等农田的面源污染风险识别，特别优选适用于我国北方平原地区农田的农业面源污染风险识别。

(三)具体的一个实施例

为进一步显示本发明的风险识别方法在农业面源污染风险识别中的识别精准度，申请人示出以下的实施例。应当理解的是，该实施例并不是对本发明的限定。

步骤1：识别目标区域的基础数据的收集

首先，根据《中国水资源公报》(2014、2015)，海河区水质为劣，Ⅲ类水及以上水占海河流域评价河长的64.6％，且其主要污染来源为农业面源污染和生活污水污染，因此选取海河流域进行农业面源污染风险识别的实例应用。海河流域总面积31.69万km²，地跨北京、天津、河北、山西、山东、河南、内蒙古和辽宁等8个省(自治区、直辖市)。流域数据分辨率为1km。

步骤2：风险识别指标体系的确立

建立表1所示的指标体系。

表1农业面源污染风险识别指标体系

步骤3：确定上述风险识别指标体系中的各指标的权重

根据上述的方法，采用层次分析法，对以上各指标赋予权重如表2所示。

表2农业面源污染风险识别指标体系及其权重

步骤4：对上述风险识别指标体系中各指标进行赋值

按照下述表3对8个指标进行分级，形成各指标的分级图，即图1～8，并将对应级别的数值范围赋予1～4分，即得到指标赋值I。

表38个指标分级标准

步骤5：将上述权重W和指标赋值I代入下式(1)中，获得农业面源污染风险指数NPSPRI；

NPSPRI＝∑W_i×I_i (1)

式(1)中，NPSPRI为农业面源污染风险指数；Wi为各风险指标在指标体系中的权重值，其值范围为0-1；Ii为指标赋值。

具体而言，对分级后的8个指标数据按各指标以及权重值代入上述式(1)中进行加权，得到农业面源污染风险指数NPSPRI。

步骤6：对上述农业面源污染风险指数NPSPRI进行自然间断点分级法(Jenks)分级为四个风险等级，详细分级标准为：无风险(0，2.208]、低风险(2.208，2.704]、中风险(2.704，3]和高风险(3，4)，即得到农业面源污染四级风险识别图(图9)；并将各个级别数值范围赋予1～4分，即：无风险为1分，低风险为2分，中风险为3分，高风险为4分。

为显示本发明的风险识别方法的风险识别精准度，下面，将综合指标风险等级图(图9)与遥感分布式DPeRS(Diffuse Pollution estimation with Remote Sensing)模型风险等级图(图10)的识别结果进行比较。DPeRS模型是一种基于二元结构的半经验半机理过程的模型，既考虑降水、植被覆盖、地形地貌等等自然因素，同时也考虑了施肥利用效率、人口、牲畜和家禽等社会经济因素，并耦合遥感技术，可以对流域尺度面源污染负荷的时空动态进行精确的定量评估。

为实现两者的可比性，对遥感分布式DPeRS模型的总负荷风险图也进行分级，这里，与图9类似地也分为四个风险等级(即无风险、低风险、中风险以及高风险)，将四个等级数值范围赋予1～4分，得到模型风险等级图，即图10。

然后，为显示采用本发明的风险识别方法进行风险识别而得到的综合指标风险等级图(图9)与DPeRS模型风险等级图(图10)之间的高度一致性，对二者的一致性进行验证。

1、验证方法说明

首先，将赋值后的两幅等级风险图(图9和图10)做差值处理，差值结果图的偏差数值情况及具体描述如下表4～5；其次，根据风险图偏差像素数与总像素个数做比值运算得出偏差百分比并输出偏差图(即图11)，公式计算为：

偏差百分比＝(偏差像素数/总像素数)×100％；

表4偏差图数值情况表

级别(数值)	一级(1)	二级(2)	三级(3)	四级(4)
					一级(1)	0	1	2	3
二级(2)	-1	0	1	2
					三级(3)	-2	-1	0	1
四级(4)	-3	-2	-1	0

表5偏差图数值定义描述表

2.偏差分析

对于海河流域农业面源污染风险采用本发明的指标体系进行识别的结果与模型识别结果的偏差分析结果如下表(表6和表7)。

表6图9与图10的风险等级面积及百分比情况表

注：1)DPeRS模型风险等级图(图10)用A表示；综合指标风险等级图(图9)用B表示；研究区用C表示；

2)数据分辨率为1000m，换算面积公式为：

面积＝(pixel*1000*1000)/1000000，

即像素数pixel为面积值，单位为km²。

表7图10与图9偏差情况表

从表6和表7可以明确看出，图10与图9中，同一风险等级面积相差不超10％，尤其是高风险级别面积相差仅为0.24％；94.46％以上的面积均为偏差小或无偏差。这表明，本发明的风险识别方法具有与遥感分布式DPeRS定量模型法同等水平的风险识别精准度。而本发明的风险识别方法不需要模型法所必须的复杂的指标以及数据积累，也不需要考虑复杂的下垫面情况。因此，本发明的风险识别方法可以实现简单、快速、低成本、适应性强、不依赖GIS和遥感技术等高精尖设备和技术且非常精确地农业面源污染风险识别。

Claims (3)

1.一种农业面源污染风险识别方法，包括如下步骤：1)测定对象地区的风险识别指标体系中各指标的数值；

所述风险识别指标体系为水文气象指标、土壤地形植被指标以及经济指标；

所述水文气象指标为年降水量、溶解态面源污染物入河系数以及吸附态面源污染物入河系数；

所述土壤地形植被指标为年植被覆盖度、坡度以及土壤可侵蚀性因子；

所述经济指标为农田氮平衡量以及农田磷平衡量；

2)计算所述风险识别指标体系中的各指标数值的权重，得到权重，记为W；

对所述风险识别指标体系中各指标进行赋值，得到指标赋值，记为I；

3)将所述权重W和所述指标赋值I代入下式(1)中，得到农业面源污染风险指数；

NPSPRI＝∑W_i×I_i (1)

式(1)中，NPSPRI为农业面源污染风险指数；W_i为各风险指标在指标体系中的权重值，其值范围为0～1；I_i为指标赋值，i为指标体系中的指标；

根据计算得到的所述农业面源污染风险指数数值，即能识别农业面源污染风险；

当所述年降水量的取值分别为≤400mm、400～500mm、500～700mm、>700mm时，其赋值分别为1、2、3、4；

当所述溶解态面源污染物入河系数分别为≤0.018、0.018～0.055、0.055～0.130、>0.130时，其赋值分别为1、2、3、4；

当所述吸附态面源污染物入河系数分别为≤0.018、0.018～0.13、0.13～0.28、>0.28时，其赋值分别为1、2、3、4；

当所述年植被覆盖度的取值分别为>60％、45～60％、30～45％、≤30％时，其赋值分别为1、2、3、4；

当坡度分别为≤8、8～15°、15～25°、>25°时，其赋值分别为1、2、3、4；

当土壤可侵蚀性因子分别为≤0.010t·h·(MJ·mm)^-1、0.010～0.023t·h·(MJ·mm)^-1、0.023～0.027t·h·(MJ·mm)^-1、>0.027t·h·(MJ·mm)^-1时，其赋值分别为1、2、3、4；

当农田氮平衡量分别为≤0t·km^-2、0～15t·km^-2、15～40t·km^-2、>40t·km^-2时，其赋值分别为1、2、3、4；

当农田磷平衡量分别为≤0t·km^-2、0～5t·km^-2、5～15t·km^-2、>15t·km^-2时，其赋值分别为1、2、3、4。

2.根据权利要求1所述的方法，其特征在于：所述年降水量根据测定所述对象地区各监测站的月降水量数据计算得出；

根据所述年降水量和测定所述对象地区各监测站的月径流量数据，通过下式(2)计算所述溶解态面源污染物入河系数：

式(2)中，CR为溶解态面源污染物入河系数，Prec和Runoff分别为年降水量和年径流量；

根据测定的泥沙含量数据和通过式(4)～(11)计算得到的土壤侵蚀量，通过下式(3)计算所述吸附态面源污染物入河系数；

式(3)中，SDR为吸附态面源污染物入河系数，Sed和Sel分别为年泥沙含量和年土壤侵蚀量；

Sel＝R×K×L×S×C×Ρ (4)

式(4)中，R为降雨侵蚀力因子，K为土壤可侵蚀性因子，L、S分别为坡长因子和坡度因子，无量纲；C为生物措施因子，无量纲；P为工程措施因子，无量纲；

式(5)中，P_i为月降雨量，P_year为年降雨量，i是月份；

根据(6)计算得到土壤可侵蚀性因子；

K＝0.1317×K_China (6)

K_China＝-0.01383+0.51575K_EPIC (7)

式(6)-(8)中，K_EPIC为EPIC模型计算得到的土壤可侵蚀性因子，K_china是中国土壤可侵蚀性因子，K为土壤可侵蚀性因子，0.1317为美国制和国际制的单位转换系数，S_a为砂粒含量，S_i为粉粒含量，C_i为粘粒含量；C是土壤有机碳含量，S_n＝1-S_a/100；

式(9)中，θ是坡度；

式(10)中，λ是坡长；

式(11)中，FVC为年植被覆盖度。

3.根据权利要求1或2所述的方法，其特征在于：根据国家县级农业统计指标数据采用土壤系统收支平衡方法通过式(12)～(15)计算所述农田氮平衡量和所述农田磷平衡量；

Q_bal＝Balance/area×1000 (12)

其中，Balance＝Input-Output (13)

Input＝Ftlz+Mnr+Irg+Seed+Dpzt+Bnf (14)

Output＝Hvst+Gas+water (15)

其中，Q_bal为所述农田氮平衡量或所述农田磷平衡量；area为耕地面积和园地面积之和；1000为单位转换系数；Balance为养分平衡量；Input为养分输入量；Output为养分输出量；Ftlz为化肥养分输入量；Mnr为有机肥养分输入量；Irg为灌溉养分输入量；Seed为种子养分输入量；Bnf为生物固氮氮输入量；Dpzt为干湿沉降养分输入量；Hvst为作物带走养分输出量；Gas为氨气挥发和反硝化氮输出量；Water为淋失和径流养分输出量。

Images