基于随机森林算法与农户视角的村落级农业面源污染风险分析

朱康文，张 晟*，陈玉成，熊海灵，雷 波，黄昌前，段秋宴

（1.重庆市生态环境科学研究院，重庆 401147；2.西南大学 资源环境学院，重庆 400716；3.西南大学 计算机与信息科学学院，重庆 400716；4.重庆市綦江区生态环境局，重庆 400803；5.重庆市渝北区生态环境监测站，重庆 401120）

农业面源污染是在农村生产与生活过程中，因作物种植、畜禽养殖、村民生活带来的有机或无机物质，通过农田的地表径流和农田渗漏等进入水体引起的污染[1-2]。2020年6月，由生态环境部、国家统计局和农业农村部联合发布了《第二次全国污染源普查公报》，显示2007—2017年中国农业源减排效果明显，但是农业源排放量占总排放量的比重仍然居高不下，其中化学需氧量（chemical oxygen demand，COD）、总氮（total nitrogen，TN）、总磷（total phosphorus，TP）排放量占比分别由2007年的约44%、57%、67%变为2017年的约50%、47%、67%[3]，因此，农业面源污染仍是当前亟须解决的问题。

针对农业面源污染问题在各空间尺度均开展了一系列研究工作。大尺度评估方法以统计数据分析及模型构建分析为主，部分研究借用GIS方法进行空间展示[4-5]；中尺度主要以GIS结合区域地形、土地利用类型、污染源、驱动因子、输出系数法等方法进行分析[6-7]；小尺度主要以GIS结合“源-汇”分析、室内模拟、水文模型等方法进行分析[8]；微尺度主要以农户调研分析等方法为主[9-10]。其中微尺度是农业面源污染政策落地的层级，同时也是反馈政策实施效果的层级。当前针对微尺度研究主要以农户调研为主，能较好地体现农户角度的农业面源污染物输移情况，但是缺乏空间的直观分析，很难体现不同地块之间的农业面源污染风险差异，导致微尺度的农业面源污染风险调控效果较差。为解决此问题，本文认为针对村级范围内泛地块尺度的风险测算研究，可以较好地融入农户角度的调研成果，为微尺度研究提供一种新的研究思路。村落和泛地块分别作为最小行政管理单元和农户行为单元，对于农业面源污染防控具有很好的实际意义，对于农户行为下的风险防控更具指导意义。当前村落范围内泛地块尺度的风险测算需要解决高精度低空遥感数据获取、微尺度泛地块范围划分两个问题，是微尺度地物识别和统计测算的重要基础数据。解决泛地块范围划分问题需要寻找获取高精度地物分类数据的技术手段。当前地物识别基本采用卫星遥感技术，因其具有数据易获取、低成本等优势，在农业领域的应用已较为常见。Nuarsa等[11]利用MODIS数据在印度尼西亚巴厘省结合NDVI、RVI和SAVI 3种植被指数提取水稻种植面积。但同时也存在分辨率较低等问题，尤其难以满足村落范围微尺度的研究。随着无人机低空遥感技术的发展，近年应用无人机多光谱技术开展作物精细分类的研究逐步增多[12]。无人机低空遥感技术在动植物调查、地质调查、精细农业等方面应用呈现爆发式增长[13-15]。高精度数据获取后的地物分类手段采用随机森林算法，它具有人为影响小、多决策优化等优点，目前在国内外地物分类中应用较多[16]，这些已有的研究成果可以为本论文提供很好的借鉴。

因此，本文选取三峡库区典型村落重庆市涪陵区南沱镇睦和村为研究对象，采用多光谱无人机技术获取多光谱遥感数据，采用随机森林算法进行精细化地物分类[17]，并通过区域作物种植类型实地调研、核查等方式校准、验证分类精度，在此基础上开展泛地块范围划分及TN、TP施用水平测算，并对农业面源污染物消减工程的成效进行分析，为村落范围内农业面源污染精细化风险防控提供依据。

1 研究区概况与数据来源

1.1 研究区概况

睦和村位于重庆市涪陵区南沱镇西北部，与长江相邻，区域内存在较大面积的库湾型区域，属于三峡库区重点移民村。经过近年来的产业结构调整，形成了规模化的果树种植模式，部分果林套作青菜头、桑树等，基本属于无粮村种植模式。由于区域生态地位较高且区域农业发展较好，因此，针对南沱镇睦和村的农业面源污染问题研究较多，已有较多研究关注统计数据分析和示范工程建设，较少关注污染风险空间差异分析。2020年11月，研究人员赴南沱镇睦和村现场调研，主要调研农业面源污染负荷排放情况。

（1）人口情况。全村人口493户，分三类情况：T1类，实际有人在家且有种植、养殖的农户，315户；T2类，集中安置的农户，144户（集中安置无养殖）；T3类，长期无人在家的农户，34户。所有生活污水排放均进入沼气池，生活垃圾收集运输至垃圾站。根据调研信息，T1和T2类型中以老年人为主，主要分4种类型：（a）老年人在家，165户；（b）中年夫妇和老年人在家，49户；（c）中年人一人在家且老年人在家，225户；（d）中年人在家，20户。基本每户都有孙子、孙女在家。不计算孙子辈，60岁以下男性121人、女性150人，60岁以上男性353人、女性278人。男女比例52.55∶47.45。

（2）养殖业情况。有养殖的农户315户，其中养猪100户，平均2头/户；养鸡300户，平均14只/户；养鹅或鸭90户，平均7只/户。粪污基本全部进入沼气池后就近还田施用。

（3）种植业情况。睦和村研究区域面积为291.04 hm2，其中种植面积185.21 hm2，户均耕地面积约0.59 hm2。其中桂圆、荔枝、枇杷、脐橙、桑树、柚子、玉米（与青菜头轮作）、花生（与青菜头轮作）、耕地水田的种植面积分别约为32.57hm2、4.40hm2、11.27hm2、98.18hm2、12.05hm2、1.85hm2、20.67hm2、2.06hm2、 2.16 hm2， 化 肥 主 要施用氮磷钾养分含量各15%的复合肥，桂圆、荔枝、枇杷、脐橙、桑树、柚子、玉米、青菜头、花生、耕地水田的化肥施用水平分别为1 012.50 kg/hm2、1 012.50 kg/hm2、421.88 kg/hm2、 506.25 kg/hm2、700.50kg/hm2、0kg/hm2、187.50 kg/hm2、468.75kg/hm2、0 kg/hm2、135.00 kg/hm2（折纯）。

图1 睦和村区位分布图Fig.1 Distribution map of Muhe Village location

1.2 研究路径与数据来源

考虑当前无人机较多存在飞行时间短、飞机过大过重、对飞行环境要求较高等问题，结合研究需求的无人机产品性能情况，选取大疆精灵4多光谱版无人机作为数据获取手段，此无人机具有飞行任务可连续进行、获取数据量丰富、续航时间较长等优势。主要研究路径如下：（1）基于获取的无人机遥感数据，采用成熟的ENVI软件扩展工具中基于决策树的随机森林算法进行地物分类；（2）开展区域作物种植类型实地调研、核查等以校准、验证分类精度；（3）结合区域地形、种植类型、路网、水系网络等进行泛地块网格划分，以利于分析各泛地块网格化肥施用、畜禽养殖等情况，测算各网格TN、TP施用水平，从污染输入角度厘清区域各泛地块网格农业面源污染风险的主要原因；（4）初步开展区域典型点位的水质与风险测算结果之间的一致性分析，以验证测算结果；（5）初步测算区域已有面源污染消减工程的水质改善情况，为区域泛地块尺度的农业面源污染风险调控提出防控策略，为精细化风险防控提供支撑。

数据来源：区域范围边界数据来自规划部门；遥感数据、NDVI（归一化植被指数）等指数数据由无人机获取，数据处理后得到NDRE（归一化差异红边植被指数）、LCI（叶片叶绿素指数）、OSAVI（优化土壤调节植被指数）、GNDVI（绿色归一化植被指数）[18]；地块作物及施肥施药信息、农户信息、养殖信息等来自实地农户调研数据；地形数据来自资源环境数据云平台①http://www.resdc.cn/；河流、道路数据通过无人机数据矢量化得到。数据采用2000国家大地坐标系统。

2 研究方法

2.1 基于随机森林算法的地物精细化分类

随机森林算法是2001年提出的一种集成学习技术的方法[19]，在多源数据进行地物分类方面已开展了大量研究[20]，在多光谱影像识别地物分类方面也进行了一定研究[21]，已有研究者以Landsat ETM+多光谱影像数据作为研究对象，对随机森林、Boosting、Bagging、CART分类方法的精度进行研究，得出随机森林的分类精度较优[16]。随机森林算法进行地物分类包含大量的决策树、分类树和回归树运算，它可以同时处理分类和回归分析问题，通常包括训练样本选取、决策树模型构建、随机森林分类、精度验证等步骤，详见参考文献[22]。

2.2 泛地块网格划分

泛地块网格划分主要依据随机森林算法提取的地物类型数据与地形数据、道路数据、河流数据等进行叠加分析，进行泛地块网格划分，分为一级网格和二级网格，其中一级网格为二级网格划分的过渡数据。本研究将睦和村划分为66个泛地块级尺度一级网格，2 710个泛地块级尺度二级网格。

2.3 污染风险测算方法

农业面源污染风险测算以污染负荷测算为基础，污染负荷测算主要包括种植源、生活源、养殖源，其中种植源的TN、TP污染负荷根据化肥施用量进行测算，生活源、养殖源的TN、TP污染负荷根据沼气池排放量进行测算。研究参考前期研究成果[23]，污染风险按照TN、TP施用强度＞280 kg/hm2、＞250～280 kg/hm2、＞225～250 kg/hm2、＞200～225 kg/hm2、≤200 kg/hm2划分为极高风险、高风险、中风险、低风险、无风险5个等级。

式中：PL为区域TN或TP污染负荷，kg；PS、LS、BS分别为种植源、生活源、养殖源的污染负荷，kg；RS为轮作作物消减的污染负荷，kg；PPi为第i种作物的化肥施用量，kg；Si为第i种作物对应化肥的折纯系数，无量纲。

LS和BS本研究中按照生活源、养殖源产生的污染物经过沼气池后排放量进行测算，根据西南大学污染修复与清洁生产评估中心在《重庆市农村户用沼气系统的环境效应与绩效评估》中的单个户用沼气池对农业面源污染的贡献结果，单个8～10 m3的沼气池每年TN、TP排放量分别相当于施用有机肥29.5 kg和12.6 kg，对睦和村户用沼气池的TN、TP排放量进行测算，同时假定所有的有机肥都均匀还田到所有地块。RS根据实际按照轮作作物消减能力进行测算。由于睦和村存在花生、青菜头轮作的情况，花生种植期间未施用化肥，因此花生实际上在轮作过程中起到了一定的污染物消减的作用。因此，综合已有结论，结合当地实际地理环境、地形坡度、土壤质地等情况[24]，当地常用复合肥养分含量为45%，其中氮、磷、钾比例为1∶1∶1。根据孙虎[25]和万书波等[26]的研究成果，花生对化肥氮素的当季吸收利用率为31.28%～42.15%，但是化肥氮素当季吸收率根据品种差异会存在一定变化。因此结合戴良香[27]、李海东[28]等人研究认为的氮肥施用最佳量分别为77.19 kg/hm2、75 kg/hm2，研究假定花生轮作可以消纳75 kg/hm2的氮素养分（折纯）。

3 结果与分析

3.1 随机森林算法识别地物结果

3.1.1 无人机多光谱数据获取结果

本研究在飞行时为避免飞行过程中辐射校正的误差，选择上午10：00至下午2：00之间飞行，飞行时间为2020年5月22日，飞行航线覆盖整个睦和村研究区域（不含江面）。获取的多光谱数据采用大疆智图软件进行影像重建等预处理工作[29-30]，得到 NDVI、NDRE、LCI、OSAVI、GNDVI 5 种指数数据（图2）。随机选取居民用房、道路、林地、草地、桑树、荔枝、桂圆、脐橙、玉米（与青菜头轮作）、花生（与青菜头轮作）、水域、裸地、其他用地等地类点位，拾取对应区域的指数结果，获取各地物各指数的纹理概略特征，为利用随机森林算法进行地物判别提供特征信息，比如桑树存在“NDVI＞GNDVI＞LCI＞OSAVI＞NDRE＞0”、草地存在“NDVI＞GNDVI＞LCI＞OSAVI＞NDRE＞0”的特征。

图2 睦和村遥感影像和多光谱指数结果分布图Fig.2 Distribution map of RS image and multispectral index in Muhe Village

3.1.2 随机森林算法识别地物与精度验证

随机森林算法识别地物结果显示，研究区域总面积为291.04 hm2，人工表面、道路、林地、草地、水域、其他用地（含裸地）比例分别为7.24%、3.08%、14.06%、7.20%、3.40%、1.38%；耕作用地总面积比例高达63.64%，其中脐橙田块的面积最高，为33.74%，其次为桂圆、玉米、桑树。从图3中的空间分布来看，林地和草地主要分布在水域周围，以1号、8号、13号、35号等一级网格分布最多。

为验证地物分类精度，通过室内判断、实地调研等形式选取181个点位进行精度校验（图3），包括玉米（与青菜头轮作）、花生（与青菜头轮作）、柚子、桑树、脐橙、枇杷、荔枝、桂圆、房屋、道路、林地（竹林、樟树、灌丛、桂花树等）、水面、耕地（水田）、草地、裸地等类型。校验结果发现181个点位分类结果中18个点位与实地核查结果不相符，分类精度为90.05%，符合地物分类精度要求，误判的主要为林地与草地、柚子与林地、玉米与草地。

图3 睦和村泛地块级尺度网格下精细化地物分类结果Fig.3 Results of fine feature classification of Muhe Village in block scale grid

3.2 一级网格的污染负荷排放与风险测算结果分析

根据一级网格的污染负荷测算结果（图4），睦和村TN、TP的总负荷量分别是70.62 t和64.73 t，其中30号网格污染负荷强度均值最小，TN负荷强度均值为89.97 kg/hm2，TP负荷强度均值为81.93 kg/hm2。47号、50号、65号网格的污染负荷强度均值最大，TN负荷强度均值分别为452.24 kg/hm2、459.83 kg/hm2、476.68 kg/hm2，TP负荷强度均值分别为420.07 kg/hm2、427.65 kg/hm2、444.51 kg/hm2。通过对66个一级网格的种植地块TN、TP负荷进行测算发现，仅有1个网格（30号）的TN、TP施用强度低于200 kg/hm2，且其他网格均高于280 kg/hm2。因此，根据风险分级标准，仅有1个网格（30号）属于无风险等级，其他网格均属于高风险等级。

图4 泛地块一级网格TN、TP施用水平分布图Fig.4 Distribution map of TN,TP application level based on first level of block grid

3.3 二级网格的污染负荷排放与风险测算结果分析

根据二级网格的污染负荷测算结果来看（图5），整体上涉及种植的图斑中TN、TP施用强度分别有124个和409个网格低于200 kg/hm2，处于无风险状态。TN施用强度来看，分别有285个、63个、1 309个网格处于低风险、高风险和极高风险状态。TP施用强度来看，分别有205个、1 167个网格处于中风险和极高风险状态。空间上存在3个风险值极高（TN、TP施用强度均高于400 kg/hm2）且聚集分布的区域：东部的54号、55号、63号等一级网格中的54-51、54-56、55-41、55-45、63-83、63-87等二级网格；中部33号、34号等一级网格中的33-18、33-25、34-43、34-48等二级网格；西部的9号、10号等一级网格中的9-118、9-127、10-50、10-59等二级网格。桑树且存在玉米和青菜头轮作的网格的TN、TP施用水平最高，主要是3-2号、12-4号、14-1号、15-37号、15-38号、15-39号、18-3号、27-2号、36-35号、36-36号、58-2号、59-2号、62-18号、62-19号、63-73号、63-74号、63-75号17个二级网格，其次是脐橙且存在玉米和青菜头轮作的网格。

图5 二级网格下TN、TP施用水平与风险等级分布图Fig.5 Distribution map of TN,TP application level and risk level based on second level of block grid

3.4 污染风险强度与区域水质之间的一致性分析

为初步探明研究测算得到的风险测度成果与实地水质之间的一致性，同时避免消落区水位变化对分析的干扰，在常年有水且不易被水淹没的龙王沟及田间沟壑等区域选取了5个点位进行水质分析，结果显示所有点位水质均较差（图6；表1、表2）。MH1点位为容量较大的堰塘，TN超标较严重、COD略高；MH2点位为田间沟壑、MH3和MH4号点位为流水溪沟，TP、TN均超标较严重；MH5点位为田间沟壑，水质最差，氨氮、TN、TP严重超标。根据各采样点汇水范围内的一级网格的土地利用类型和TN、TP施用强度统计结果来看，所有采样点TN、TP施用强度均较高，林草类型面积占比为11.58%～26.74%，占比较低。MH5号点位超标严重与林草类型面积占比低及TN、TP施用强度高密切相关。MH4号点位，虽然同样具有高的TN、TP施用总量，但由于区域内林草类型面积占比相对较高，因此水质相较MH5号点位改善较多。结果表明睦和村由于化肥施用强度高、林草比例低、林草缓冲带少等原因导致区域内水质情况较差。整体上来看，各区域污染负荷（风险强度）与点位水质质量呈现明显的一致性，表明污染负荷（风险强度）测算结果与实地相符。

图6 睦和村水质采样点位的空间分布图Fig.6 Spatial distribution of sampling points of water quality in Muhe Village

表1 睦和村水质测试结果Table 1 Water quality test results of sampling points in Muhe Village

表2 采样点汇水涉及一级网格土地利用、TN和TP施用情况Table 2 Land use,TN and TP application in the first level grid of catchment area of sampling points

4 讨论

4.1 低空遥感与随机森林算法结合适用于地物精细化分类

目前低空遥感数据主要获取源无人机具有很好的机动性、分辨率高、高频次等优点，在小尺度研究中使用得越来越多。低空遥感数据结合随机森林算法的研究逐渐增多，且都取得了较好的研究效果。耿仁方等[31]在广西桂林会仙喀斯特国家湿地公园通过无人机遥感影像和随机森林算法开展岩溶湿地植被识别研究取得了较好的进展。任泽茜[32]在南方平原地区利用无人机和随机森林算法进行农作物种植面积监测。这些已有的研究均表明低空遥感数据和随机森林算法结合可以很好地支撑小区域地物精细化识别、作物种植面积、作物状态等方面的研究工作。本文针对地物分类结果的现场核实情况也反映了随机森林算法在精细化地物识别方面的优势。

4.2 人工湿地构建有利于降低污染风险

结合睦和村已经开展的人工湿地去污植物筛选与群落优化配置工作（图7），对人工湿地入水口、出水口水质进行分析，发现出水口水质明显优于入水口。人工湿地的污染物消减程度非常高，悬浮物、TP、TN、铵态氮、硝态氮、生化需氧量分别消减了近86.67%、54.66%、81.11%、10.67%、83.85%、59.42%。表明人工湿地对于削减农业面源污染的排放强度、风险程度是极有利的方式。同时结合研究前期开展的植被调查等基础工作，认为水位高程是消落带植被恢复最重要的因素。三峡库区蓄水初期自然再生植物狗牙根（Cynodon dactylon）、香附子（Cyperus rotundus）对于阻控、减缓区域农业面源污染具有很好的效益，而且植物篱这种生态修复方式对于农业面源污染有较好的防控效果，同时乡土植物桑树（Morus alba）、黄荆（Vitex negundo）等植物也有一定的耐水性和较好的污染物吸附能力[33]。这与潘杰[34]、余顺慧[35]、莫福孝[36]等学者的研究结果基本一致。

图7 人工湿地构建区域低空遥感图Fig.7 Remote sensing map of the constructed wetland area

5 结论

（1）无人机多光谱数据获取手段结合随机森林算法可有效用于地物精细化分类，睦和村地物分类精度高达90.05%。

（2）基于农户行为调研角度与泛地块网格划分框架下的污染负荷及风险测算，有利于识别村落尺度下的地块级污染负荷及风险差异。睦和村TN、TP的总负荷量分别是70.62 t和64.73 t，一级网格测算发现仅有1个网格（30号）属于无风险等级，其他网格均属于高风险等级；二级网格测算发现涉及种植的用地中，分别有124个、285个、63个、1 309个网格的TN处于无风险、低风险、高风险和极高风险状态，分别有409个、205个、1 167个网格的TP处于无风险、中风险和极高风险状态。研究通过实地采样发现区域污染负荷（风险）强度与区域水质质量存在较好的一致性，侧面验证了测算结果的真实性。

（3）人工湿地构建在大幅消减污染物的同时可产生较好的经济效益，区域内的悬浮物、TP、TN、铵态氮、硝态氮、COD在人工湿地的作用下分别消减了近86.67%、54.66%、81.11%、10.67%、83.85%、59.42%。

研究认为劳动力下降、人口老龄化等问题导致作物复种指数高、高肥高药行为突出，是引起农业面源污染风险较高的主要原因。污染负荷或风险较高的地块在未来农户行为管理中应重点关注并对其进行指导。研究结果可为村落级管理人员在指导农户开展农业生产活动时提供参考，为相关人员制定区域农业面源污染风险防控策略提供依据。未来研究将深入分析轮作最佳模式、人工湿地削减污染的微观机制等，以更好地提升成果的实用价值。