华北地区夏玉米生产中农田氮淋失的定量预测*

赵晓莹,王诺婷,崔 斌,尹实磊,杨 轩,孟凡乔

(中国农业大学资源与环境学院/农田土壤污染防控与修复北京市重点实验室 北京 100193)

随着人口的不断增长,2050 年全球将面临满足90 亿人口粮食需求的巨大挑战[1],集约化农业生产中增加氮肥投入和提高农田灌溉覆盖率是粮食增产的重要保障。改革开放以来,中国农业集约化程度不断加强,作物产量不断提高。2020 年,中国玉米(Zea mays)总产量为26 066.5 万t,氮肥施用量达1833.9万t,分别为1980 年的4.16 倍和1.96 倍[2]。化肥和灌溉水投入快速增长,一方面可以保证作物高产稳产,另一方面也引发了氮淋失等面源污染问题[3]。中国玉米作物氮肥利用率仅为31.0%[4],大量氮肥会通过NH3挥发、N2O 排放、硝酸盐和有机氮淋洗等途径损失[5],对空气、水体质量以及全球气候产生负面影响。资料显示,中国夏玉米生长季硝态氮平均淋失量为27.6 kg(N)·hm-2,占氮肥施用量的13.3%[6]。农田氮流失或淋失到陆地水体中,会导致地下水硝酸盐含量升高,对地下水及饮用水安全构成威胁[7]。殷乐宜等[8]针对京津冀重点区域-密怀顺平原区地下水的监测表明,该区地下水硝酸盐含量为0.24～178.0 mg·L-1,平均值为47.3 mg·L-1,远超《地下水质量标准GB/T 14848-2017》Ⅲ类标准,其主要原因在于该区域化肥长期高量施用。分析和量化华北地区农田氮淋失及其影响因素,有助于提高肥料利用率、降低氮素损失,对于中国农业绿色发展和生态文明建设具有重要意义。

农田氮素淋失受多个环境条件和生产管理措施的影响,且不同因素对氮淋失的重要程度各异[9]。目前,大多数研究主要从氮肥和水分两个方面,应用单一模型进行氮淋失量定量模拟。Zhou 等[10]对玉米-小麦(Triticum aestivum)种植体系的硝态氮淋失进行了Meta 分析,夏玉米季硝态氮淋失(Y)与施氮量(X)的线性方程为Y=3.83+0.25X,拟合优度为0.51,即氮淋失与施氮量之间存在着正相关关系。随着研究的深入和增加,近年来氮淋失的预测分析也尝试了其他定量关系,解释变量也增加了氮盈余量等指标。Wang 等[11]建立的夏玉米季总氮淋失量(Y)与氮盈余量(X)之间的指数模型为Y=12.14e0.0061X,拟合优度为0.35;而Cui 等[12]构建的指数模型为Y=2.38e0.0041X,拟合优度为0.63。由于水投入直接影响水分渗漏量、进而影响氮淋失量,胡玉婷等[13]构建的旱地硝态氮淋失量(Y)与灌溉水量(X)间的线性方程为Y=4.2936+0.046X,以及总氮淋失量(Y)与降水量(X)之间线性方程Y=4.1463+0.0192X。随着近年来“双减”工作的开展,“十三五”期间我国在冬小麦和夏玉米生产中开展了大量的定量观测试验,为区域尺度氮淋失定量提供了数据支撑。需要指出的是,氮肥施用和水分渗漏对农田氮淋失的影响较大,定量预测时应对这两个因素重点予以考虑。除了施肥和水分,其他自然条件和农田管理措施也会对氮淋失有影响。过去的研究大多局限于单一因素,对不同环境条件和管理措施的影响程度缺乏分析和比较,在定量关系上对于线性和非线性模型也认识各异。

本研究整合了1980-2021 年底在我国华北地区开展的夏玉米氮淋失试验,选取自然条件、农田管理措施等为自变量,基于线性模型、指数模型、多项式模型和多元回归模型进行定量估测与比较,探索不同氮淋失影响因子的重要程度和影响机制,为华北地区农业面源污染防治提供科学依据。

1 材料与方法

1.1 研究区域

本文研究区域为华北地区。传统意义上华北地区包括北京、天津、河北、山西、内蒙古[2],由于气温、降水、日照和栽培方式等因素和措施较相似,将河南、山东、陕西等省份的研究也列入本次研究范围。

1.2 数据来源

以“华北地区” “氮素平衡” “氮淋失” “夏玉米”等以及相应英文为关键词,在“中国知网”和“ISI-Web of Science”等文献平台上进行检索,从中选取满足本研究条件的文献。所筛选的文献应满足以下条件:1)文献的试验年限应为1980-2021 年。2)氮淋失数据应从大田试验直接测定,测定方法包括但不限于淋溶盘法、渗漏池法、陶土头法等。3)文献中试验研究应基于至少一个完整的夏玉米生长季。4)对于在同一地点多个年度的试验和文献,本研究选取最后一个年度;对于同一试验既出现在期刊论文又出现在研究生论文,本研究只选取其中一篇。5)研究区域需设置在前文定义的华北地区,即北京、天津、河北、山东、陕西、内蒙古等地。

通过列表法、作图法与经验公式法等对文献中的重要信息及数据进行提取,获得土壤有机质、全氮、全磷、全钾、质地(黏粒含量、砂粒含量),农田管理措施数据,包括肥料(总施氮量、氮肥基追比、氮盈余)、水分(水分投入量、水分渗漏量、水平衡)和作物(作物秸秆和籽粒产量、地上部吸氮量),以及淋失的总氮及无机氮量等数据。

本研究数据集共包含26 篇期刊文献与硕士和博士论文,共计291 条数据。土壤类型包括潮土、褐土、土、潮褐土及棕壤。预测变量中的总施氮量为无机氮肥和有机肥输入量,氮盈余量等于总施氮量减去地上部吸氮量。水分投入量为降雨量加灌溉量,水平衡等于水分投入量减去水分蒸散量,试验研究中蒸散量主要通过水文学方法、微气象方法、植物生理学方法测定。文中利用氮盈余和水平衡作为解释变量预测氮淋失量时,由于这两类变量部分数据为负值,本文采用统一的数学处理方法,将全部数据转为正值后进行模拟分析(氮盈余数值加150,水平衡数值加50)。

1.3 数据预处理

1.3.1 缺失值处理

对数据分析表明,缺失值主要集中在土壤全磷、全钾、铵态氮和硝态氮淋失等变量。土壤磷和钾含量对土壤氮迁移转化影响较小[14-15],且两类数据缺失量较大,本研究不考虑土壤全磷和全钾。本次定量分析,只对总氮淋失量进行分析。预处理后,可用数据为263 条。考虑到旱地作物生产中淋失氮以硝态氮为主[16],对于仅有硝态氮淋失的研究,按照硝态氮占总氮比例的95%计算总氮淋失量。

总氮淋失量的数据分布直方图如图1 所示。总氮淋失量中最大值为173.7 kg(N)·hm-2,最小值为0.11 kg(N)·hm-2,中位值为12.4 kg(N)·hm-2,平均值为22.8 kg(N)·hm-2,说明氮淋失的数据偏向于低淋失量端。

图1 夏玉米农田总氮淋失量数据分布直方图Fig.1 Data distribution of total N leaching loss in summer maize fields

1.3.2 离散变量编码

1.4 研究方法

本研究采用Origin (OriginLab Origin 2018)和Microsoft Excel (2019)软件进行线性模拟、指数模拟及多项式模拟,获得拟合方程。采用R 语言(R version 4.1.1)中的lm 函数,进行氮淋失量与两预测变量或多预测变量间的多元回归分析,从拟合优度(R2)和显著性水平(P-value)两方面进行评估。

2 结果与分析

2.1 总氮淋失量与预测变量间的相关性分析

本研究数据集中预测变量较多,直接进行模型建立存在一定的主观性和盲目性。因此,首先采用Pearson 相关系数来衡量各变量间的相关性,相关系数的绝对值越接近1,代表相关性越强。

从相关性的强弱来看(图2),总氮淋失量与水分渗漏量、总施氮量、氮盈余、土壤全氮有较强的相关性,与有机质含量、水分投入量、水平衡等的相关性次之,与土壤pH 和秸秆是否还田相关性较弱。从相关性的方向来看,总氮淋失量与氮肥、水分等为正相关,与土层深度、秸秆还田等为负相关。

图2 夏玉米农田总氮淋失量与预测变量的相关关系Fig.2 Correlation between total N leaching loss and predictive variables in summer maize fields

2.2 总氮淋失量与预测变量间的回归分析

2.2.1 线性回归、指数回归及多项式回归

基于对各变量进行的相关性分析可知,在影响总氮淋失量的因素中氮和水所占权重较大,故对总氮淋失量与总施氮量、氮盈余、水分投入量、水平衡、水分渗漏量分别进行模拟预测,探索氮淋失量与相应影响因素间的定量关系。

夏玉米生产中总施氮量范围为0～573 kg·hm-2,中位数为186.3 kg·hm-2,平均值为185.5 kg·hm-2。研究发现(图3),随着氮投入量的增加,由氮淋失造成的肥料氮损失量也呈快速增加态势。采用总施氮量预测氮淋失量时,指数模型(Y=8.88e0.0044X)、线性模型(Y=0.13X-1.07)和多项式模型(Y=0.000 24X2+0.028X+6.35)对因变量总氮淋失的解释程度均不超过30%,但指数模型呈现出优于其他两种模型的预测趋势。采用氮盈余预测氮淋失量时,指数模型、线性模型和多项式模型拟合优度相差不大。

图3 夏玉米农田总氮淋失量与总施氮量、氮盈余间的回归模拟Fig.3 Regression simulation of total N leaching loss with total N rate and N surplus in summer maize fields

总氮淋失量与水分渗漏量表现为正相关关系,氮淋失会随水分渗漏量的增加而增加(图4)。总体上线性方程(Y=0.1935X+4.6084)、指数方程(Y=4.2982e0.0109X)及多项式方程(Y=0.0005X2+0.0868X+7.7906)的拟合优度在同一水平,其中拟合度最高的为指数模型,为0.2937,最低的为线性模型,为0.2648。

图4 夏玉米农田总氮淋失量与水分渗漏量间的回归模拟Fig.4 Regression simulation of total N leaching loss with water percolation in summer maize fields

总氮淋失量与水分投入量、水平衡量间的线性拟合、指数拟合及多项式拟合的决定系数(R2)分别介于0.0128～0.0206 和0.0037～0.0278,回归关系较差(图5)。二次函数(Y=-0.000 17X2+0.1989X-30.228,Y=-0.000 34X2+0.1647X+8.498)表明,随水分投入量或水平衡量的增加,总氮淋失量呈现先增加后降低的趋势。

图5 夏玉米农田总氮淋失量与水分投入量、水平衡间的回归模拟Fig.5 Regression simulation of total N leaching loss with water input and water balance in summer maize fields

2.2.2 多元回归分析

首先,建立总氮淋失量与总施氮量、氮盈余、水分投入量、水平衡、水分渗漏量等预测变量间的多元回归模型,各模型的决定系数(R2)及回归方程如表1 所示。

表1 夏玉米农田总氮淋失量与氮水两预测变量间的多元回归模拟Table 1 Multiple regression simulation between total N leaching loss and pairwise indexes of N and water in summer maize fields

对比不同预测变量组合发现,决定系数R2最高的组合为总施氮量+水分渗漏量(0.3413),最低为氮盈余+水平衡(0.2188)。不管是总施氮量与3 个水分间的组合还是氮盈余与3 个水分间的组合,均是以水分渗漏量获得的拟合优度最高。用总施氮量和水分渗漏量两个因素预测总氮淋失量时所得回归方程为Y总氮淋失量=-5.46+0.06X总施氮量+0.19X水分渗漏量。

对总氮淋失量与土壤条件(有机质、全氮、黏粒、砂粒、pH、深度)以及农田管理措施(秸秆还田、施氮及水分等)建立多预测变量的多元回归模型(表2)。

表2 夏玉米农田总氮淋失量的多元回归模型参数Table 2 Parameters of multiple regression model for total N leaching loss in summer maize fields

多元回归模型的修正拟合优度(R2=0.4150)高于线性模型、指数模型、多项式模型及两预测变量间的回归模型。估算方程为Y总氮淋失量=-18.13+0.67X有机质含量+13.83X土壤全氮+0.40X黏粒含量-0.14X砂粒含量-0.61X土壤pH-0.03X土层深度-4.45X秸秆还田+0.06X总施氮量+0.16X水分渗漏量,表明总氮淋失量与总施氮量、水分渗漏量、有机质含量、土壤全氮和黏粒含量呈正相关,与土层深度、砂粒含量、土壤pH、秸秆还田呈负相关。回归方程建成后,需要通过线性关系和回归系数的显著性检验才可用于预测估算。经检验,多元回归方程F统计量的P值＜0.01,达极显著水平,表明总氮淋失量与预测变量之间存在极显著的线性关系。总施氮量、水分渗漏量、砂粒含量和黏粒含量回归系数的P值小于0.05,回归系数通过了显著性检验,具有统计学意义。实际工作中,水分投入量(灌溉+降雨)的获取比较容易,我们采用水分投入替换水分渗漏后多元回归方程为Y总氮淋失量=-28.56+0.45X有机质含量+5.58X土壤全氮-0.14X黏粒含量-0.31X砂粒含量+2.02X土壤pH-0.03X土层深度-10.22X秸秆还田+0.13X总施氮量+0.04X水分投入量。

由于每个预测变量对总氮淋失量的影响程度不同,选择全部变量构造方程,效果并不好,表3 采用逐步回归选择合适的变量以建立“最优”的回归方程。

表3 夏玉米农田总氮淋失量的逐步回归模型参数Table 3 Parameters of stepwise regression model for total N leaching loss in summer maize fields

逐步回归后保留的预测变量为总施氮量、水分渗漏量、砂粒含量、黏粒含量和有机质含量,模型的修正拟合优度几乎不变(R2=0.4140),但各项系数的显著性均有提升。多元逐步回归方程为Y总氮淋失量=-23.07+1.14X有机质含量+0.34X黏粒含量-0.13X砂粒含量+0.06X总施氮量+0.18X水分渗漏量。经检验,回归方程F统计量的P＜0.01,达极显著水平,各预测变量的回归系数均达到显著水平(P＜0.05),通过了线性关系显著性检验和回归系数显著性检验。基于预测变量的标准化回归系数判断其对总氮淋失量影响程度的大小,有机质含量、黏粒含量、砂粒含量、总施氮量和水分渗漏量的标准化回归系数分别为0.18、0.11、0.07、0.23、0.31,可以发现水分渗漏量的影响作用最重要,总施氮量和有机质含量次之。将水分渗漏量替换为水分投入后,逐步回归方程为Y总氮淋失量=-18.60+0.64X有机质含量-10.27X秸秆还田-0.30X砂粒含量+0.13X总施氮量+0.04X水分投入量。

3 讨论

3.1 施氮对夏玉米氮淋失的影响

氮淋失是指土壤和肥料中氮素随灌溉或降水向下运移至作物根系活动层以下,不能被作物根系吸收所导致的损失[17]。土壤中氮含量过高时氮素会在土壤剖面累积,多余氮素便会向下运移,造成氮淋失[18]。本研究中平均施氮量为185.5 kg(N)·hm-2,氮淋失平均为22.8 kg(N) ·hm-2,大约占施氮量的12.3%,与Ceccon 等[19]的研究结果类似。除施肥量外,施肥种类、施肥期、施肥方式均会影响氮素淋失特征,今后定量研究应予以关注[9]。

本研究发现氮淋失与氮肥数量之间,指数模型的拟合优度高于线性模型和多项式模型,表明氮肥超过一定阈值后,多余的氮淋失风险极高,这可以从近年来的田间试验中不断得到支持[11,20],表明优化氮肥数量是华北地区夏玉米生产应主要关注的问题。本试验构建的指数方程各项参数,与我们团队前期研究整合获得的参数范围较为一致(系数变化范围 为2.38～42.0,指 数变 化范围为0.0041～0.0095)[3]。需要指出的是,近年研究获得的氮淋失估测数量比早期估测数量低,主要表现在指数方程各参数数值减小,如本研究总施氮量与氮淋失量的指数方程Y=8.88e0.0044X的参数低于Wang 等[11]构建的方程Y=25.87e0.0072X。

氮肥施用后,除了作物吸收氮外,理论剩余氮(即氮盈余,氮素输入量与作物吸氮量差值)比施氮量更能准确预测氮淋失[11,21]。本研究却发现施氮量获得的拟合效果优于氮盈余,原因可能是在实际生产中对总施氮量定量较准确,而氮盈余是由总施氮量减去地上部吸氮量(需要分别测定秸秆和籽粒生物量和含氮量)获得,数据误差较大。

3.2 水分对夏玉米氮淋失的影响

水分投入通过影响土壤水分状况和运动特征对氮淋失产生作用[22],一般而言,氮淋失与降水量、灌溉量呈正相关。康爱林等[23]发现漫灌施肥处理导致无机氮在深层土壤累积,增加氮素淋溶风险,滴灌施肥则有利于氮素在表层土壤积累,降低氮淋失。然而,综合前期研究发现,目前对于水分与氮淋失之间关系的定量研究远远少于与氮肥关系的研究。

本研究发现,总氮淋失量与水分投入、水平衡量的回归关系较差,模型解释程度不足3%,我们猜测可能是由于农业生产中水分淋失同时受到水分投入(灌溉、降水)以及水分损失(蒸发、作物蒸散)等因素的共同影响[24]。虽有研究提出直接采用降水量等可以取得较好预测效果[25-26],但从机理上分析,如果不考虑水分通过蒸腾和蒸发造成的水损失,就很难获得准确氮淋失数量,因此基于水分渗漏量预测农田氮淋失效果更好(R2=0.2937)。水氮联合模型(HYDRUS-1D)、脱氮分解模型(DNDC)和水热碳氮模拟模型(WHCNS)等过程模型中虽然考虑了水分对氮淋失的影响,但其更多是通过对水分渗漏量累积进行计算[27-29],其精度与过程模型模拟效果直接相关。对于水分渗漏的定量,除了直接测定,可以采用经验公式[30-31]或田间水量平衡模型进行计算[32],但这些方法要么缺乏严格物理意义,要么在实践中需要花费较长时间和精力,今后研究中需要重点关注。本研究中,考虑到水分渗漏量需要田间测定、操作复杂,还提供了用水分投入量(降雨+灌溉)替代水分渗漏量构建逐步回归方程,能够提高方程的可用程度,当然模拟精度会受到影响。

总体上水分渗漏量与氮肥投入一起决定了氮淋失量。因此,与以往诸多研究直接选取施肥量进行预测不同,本文将氮(总施氮量、氮盈余)、水(水分投入、水平衡、水分渗漏量)两两组合以探究水氮的最优预测组合,进而用于多因素的多元回归模型构建。研究发现,将总施氮量和水分渗漏量组合预测氮淋失的拟合优度最高,优于其他5 个预测变量的组合。

3.3 其他因素对氮淋失的影响

土层深度是定量氮淋失的关键因素。通常认为作物根区氮素具有较好的生物有效性,而根区以外的氮素则很难被作物利用,氮淋失研究深度以根区下边界为界较为合理(水分渗漏的研究深度与氮淋失深度一致)[33]。丁燕等[34]的研究将氮淋失测定深度定义在0～90 cm,林立等[35]将之定义在0～180 cm。也有学者将根际范围定义为0～200 cm[36],刘学军等[37]通过试验证明100～200 cm 土层的氮难以被作物吸收利用。本研究中,氮淋失研究深度最浅为50 cm,最深为200 cm,一方面可以认为这个深度范围的氮很难被作物吸收利用,可以认为是“淋失”的氮,另一方面在模型计算中我们将深度作为变量参与预测氮淋失数量,因此获得模拟方程能够全面说明深度等因素对于氮淋失的影响。

农田氮素淋失同时受到各类因素的影响,除施肥和水分外,土壤条件、秸秆还田、土层深度等也会对氮淋失有影响。土壤有机质含量的高低直接影响硝化作用的强弱[38],硝化作用强容易造成硝态氮积累增加淋溶风险[39]。砂粒含量多的土壤通气透水性能良好,容易发生氮淋失,沙壤土中氮素淋溶损失量可以达到施氮量的16.2%～30.4%[40]。研究发现,将输入变量替换后,多元回归方程中黏粒含量和土壤pH系数的正负(即增加或减少)发生了变化。对于黏粒含量,考虑到土壤质地由3 种粒径组成,粒径大小配比合适的土壤孔隙率低,氮淋溶发生风险较低,但过沙土壤疏松多孔、孔隙率很高,黏重土壤则会产生缝隙进而发生优先流,都会加剧氮素的淋失。土壤pH 系数正负变化,主要是由于本研究文献绝大多数位于我国华北地区,土壤偏中性或碱性,pH 变化范围较小(6.2～9.1),变量替换后其方程计算导致系数正负发生变化。

需要特别强调的是,作物秸秆作为外源有机碳,还田后可以改善土壤理化性质、对氮肥表现为固持[41-42],从而降低农田养分淋失,与本研究结果一致,即秸秆还田对氮淋失具有抑制效果。这与我们团队在山东省桓台县秸秆还田的田间试验有较好的吻合:与2010 年相比,连续4 年作物秸秆还田可以显著降低玉米季氮淋溶量[43]。因此,在华北地区农业生产中推广应用秸秆还田,对于防治氮淋失具有重大意义。

田间试验测定氮淋失费时费力,计算机模型是解决该问题的一个重要途径[35]。目前,农田氮素损失主要以过程机理模型和经验统计模型为主,前者包括DNDC、氮素淋洗经济分析模型(NLEAP)、土壤营养物淋失模型(LEACHM)和水氮管理模型(WNMM)等,这些模型准确度高,但是所需参数较多,过程复杂。经验统计模型方法包括回归分析、相关性分析、线性混合分析等,所需参数较少,计算方法简单[21]。本研究采用逐步回归优化模型,选择赤池信息准则(AIC)值最小的方程作为最优回归方程,逐步回归后模型不存在多重共线性,方程和回归系数均通过了显著性检验,所获得的方程对于华北地区氮素管理具有重要技术支撑。近年来,计算能力强、可以考虑多因素复杂影响的机器学习优势凸显,随机森林算法作为一种机器学习方法,具有准确率高、防止过拟合能力强、能对变量进行筛选等特点[44],可以尝试探索应用。

4 结论与展望

华北地区夏玉米生产中农田氮淋失受多种环境条件和管理措施的影响。在众多土壤条件和管理措施中,氮投入、水分对氮淋失影响较大,其次为土壤有机质和土壤质地,秸秆还田能有效防止氮淋失。指数模型比线性模型更适合分析氮淋失与氮肥、水分间的定量关系,总氮淋失量与水投入的回归关系较差,但与水分渗漏数量回归关系较好,准确定量水分渗漏是今后华北地区预测氮淋失工作中应重点关注的参数。