基于EFAST的CERES-Wheat模型土壤参数敏感性分析

崔金涛 丁继辉 YESILEKIN Nebi 邓 升 邵光成

(1.河海大学水文水资源与水利工程科学国家重点实验室， 南京 210098； 2.河海大学水文水资源学院， 南京 210098；3.河海大学农业科学与工程学院， 南京 210098； 4.佛罗里达大学农业生物工程学院， 盖恩斯维尔 32608；5.江西省水利科学研究院农村水利研究所， 南昌 330029)

0 引言

随着农业信息技术的兴起，作物模型逐渐成为数字农业和现代农业研究的热点[1]。自提出植被冠层光能截获作物生长模型[2]后，经过50多年的发展已经涌现了大量的作物模型，其中较为成熟的作物生长模型主要有DSSAT(Decision support system for agrotechnology transfer)[3-4]、APSIM (Agricultural production system simulator)[5-6]、AquaCrop[7]、WOFOST(Worldfood studies)[8]和STICS (Simulateur multidisciplinaire pour Ies cultures standard)[9]等。与传统田间试验相比，作物生长模型借助一系列数学公式，通过对农田生态系统中天气变化、养分循环、作物生长等过程进行模拟，为作物栽培、水肥管理、生态环境调控等方面提供决策，具有省时、省力、易重现的优点。因此，作物生长模型的应用越来越广泛。

农业生态系统运行过程十分复杂，涉及众多的物理过程及生化反应，因此作物模型在拥有众多优点的同时，也面临着参数过多、难以获取的问题[10-11]。以DSSAT模型为例，仅要求输入的土壤参数就有10多个，如果再考虑土壤异质性进行分层，则需要输入的土壤参数数量会成倍增加。虽然土壤参数具有明确物理含义，可以根据田间实测数据确定，但是在实际应用中，由于田间土壤的异质性，一些土壤参数的测量值可能变化范围较大[12]，同时一些参数的测量极其复杂繁琐，很难对模型所涉及的全部参数进行测量和标定[13]。通常只测量和标定对模型模拟结果影响较大的参数，而固定或简单处理对模型模拟结果影响较小的参数[11,13-15]。因此，筛选出对模拟结果敏感的关键控制参数，从而减少田间实测参数数量、降低模型输入参数获取难度，对模型参数率定及后续模型本地化、区域化应用至关重要。

敏感性分析可以确定参数对模拟结果的影响，从而筛选出对模拟敏感度较高的参数。全局敏感性分析不仅可以同时检验多个参数变化对模拟结果的影响，还评估了参数对模拟结果的直接和间接影响，因此广泛应用于地学、农学等领域的模型参数敏感性分析。近年来，不少学者已将全局敏感性分析法成功应用到作物生长模型的参数敏感性分析中[10,13,16-20]。然而，已有研究大多集中于品种参数和田间管理对模拟结果的敏感性分析，研究目的多服务于模型调参率定，而对减少模型土壤参数实测数量及土壤参数简单处理的研究尚不多见。近年来，在促进农业生产的同时，开始注重对肥料的高效利用和对农业面源污染的防治[21-22]。氮肥作为一种在农业生产中最为常见、施用最多的肥料[21]，常通过氨挥发、淋溶、径流、硝化-反硝化作用等途径损失。土壤状况是影响氮素运移、转化过程及不同形态氮素分布的重要因素之一[23]，研究氮素在农业生产中的运移、转化过程及土壤对该过程的影响机理，有助于指导人们科学合理施肥[21]，模型氮素分布敏感性研究可为模型在肥料高效利用及面源污染防治方面的应用提供理论参考。然而，已有研究大多侧重于模拟作物生长及产量方面，鲜见涉及氮素运移、转化过程的研究。

本文以江苏省涟水县水利科学研究站为研究区，根据冬小麦田间实测数据，应用EFAST法分析CERES-Wheat模型中的土壤参数对冬小麦生长及氮素分布等模拟的敏感性，旨在降低模型土壤参数输入要求，为模型本地化应用提供技术支撑。

1 材料与方法

1.1 CERES-Wheat模型及参数

CERES-Wheat模型是DSSAT-CERES模型中的一个子模块，以天为步长动态模拟冬小麦在不同的气象、土壤及灌溉、施肥等田间管理下的生长发育[24]。模型通过计算作物冠层截获有效光合辐射来驱动，涉及水分吸收、光合过程、干物质分配和植株生长以及衰老等基本生理生态过程，模拟作物物候、叶面积变化、生长发育及生物产量形成，并在模拟过程中考虑生长过程中温度、水分和氮素等胁迫因素影响[24-25]。模型运行需要输入气象数据、土壤数据、作物参数和田间管理信息4部分[3,24]。气象数据至少包括模拟期内每日最低和最高气温、每日降雨量和太阳辐射量；土壤数据包括土壤容重、孔隙率、凋萎系数、田间持水率等土壤特性参数，并采用分层形式表达每一层土壤物理、水力、形态结构等特征；作物参数包括生态型、物种型和品种型3种参数，描述作物各器官光合作用、呼吸作用等生理过程，一般在模型率定调参时仅考虑品种型参数，不建议修改生态型和物种型参数[26]；田间管理信息包括种植日期、种植密度、行距、灌溉、施肥类型及施肥量等信息[27]。

本研究所需的试验数据来源于江苏省涟水县水利科学研究站(119.27°E，33.78° N，海拔7.5 m)。该地区属于暖温带季风性半湿润气候，多年平均气温14.8℃，日照时数2 280 h，降雨量时间分布不均，多集中于5—9月，多年平均降雨量为979 mm，多年平均蒸发量为1 056 mm。模型所需的气象数据来源于中国气象数据网(http:∥data.cma.cn)淮安站。气象数据主要包括日最高气温、最低气温、降雨量及日照时数等。模型所需的太阳辐射量根据获得的日照时数由埃斯曲郎经验公式估算得到[28]。模型所需的土壤数据主要来源于田间试验的实测数据及模型中相同土壤类型的默认土壤数据。供试区耕地层土壤为砂壤土，0～30 cm土壤容重为1.4 g/cm3，田间持水率为28.7%(质量含水率)。模型所需输入的土壤参数及取值范围如表1所示，土壤参数的上下限为试验区不同地块实测土壤参数最大值和最小值。

表1 选取的DSSAT模型中土壤参数Tab.1 Definition and intervals of soil parameters in DSSAT model

模型所需的作物参数来源于文献[29]，该文献调参并验证了江苏省涟水县水利科学研究站种植的淮麦11的品种参数值，分别为春化天数(P1V)35.57 d，光周期影响因子(P1D)94.81%，灌浆期积温(P5)707.2℃，开花期冬小麦单位冠层生物量籽粒数(G1)28.86粒/g，标准籽粒质量(G2)23.52 mg，成熟期标准条件下单蘖质量(G3)1.228 g，叶热间隔(PHINT)60.00℃。模型所需的田间管理数据来源于试验的观测记录和调查。冬小麦于2006年10月15日播种，行距15 cm，播深5 cm，播种密度为400粒/m2。田间管理措施如灌溉、施肥、除草、病虫害防治等按当地正常田间管理进行。

本研究主要考虑模型的2类输出，即冬小麦生产指标和氮分布指标。其中生产指标包含地上生物量(CWAM)、产量(HWAM)、收获指数(HIAM)；氮分布指标包括表征植物各器官中氮素含量的4个输出量，即地上生物量中氮素含量(CNAM)、籽粒中氮素含量(GNAM)、根中氮素含量(RNAD)及叶中氮素含量(LNAD)，表征土壤中氮素含量变化的4个输出量，即作物吸收氮素量(NUCM)、土壤淋失氮素量(NLCM)、土壤中硝态氮含量(NITD)及铵态氮含量(NHTD)，表征土壤中氮素转化的4个输出量，即氮矿化量(NMNC)、硝化量(NITC)、反硝化量(NDNC)及氨挥发量(AMLC)。

1.2 EFAST方法

扩展傅里叶幅度检验(EFAST)法是一种基于方差分解的全局敏感性分析方法，结合了FAST和Sobol法的优点改进而来。其算法简单介绍如下：

设有模型y=f(x1,x2,…,xm)，通过傅里叶转换可将其转换为y=f(s)，转换函数为

(1)

(2)

(3)

(4)

式中ωi——参数xi的振荡频率，i=1,2,…,m

φi——每个参数xi的随机初相位，取[0,2π]

s——标量变量，取[-π,π]

p——傅里叶变换参数

Ap、Bp——傅里叶振幅

参数xi变化引起的模型输出方差Vi为

(5)

(6)

p∈Z={-∞,…,-1,1,…,+∞}

模型总方差为

(7)

对标量s在[-π,π]中等间隔取样，傅里叶振幅Ap和Bp近似计算式为

(8)

(9)

式中Ns——取样总数

sk——标量s的第k个取样值

模型总方差可分解为

(10)

式中Vi——参数xi自身变化引起的方差

Vij——参数xi通过参数xj作用贡献的方差

V12…m——参数x1通过其余m-1个参数相互作用贡献的方差

通过归一化处理后，参数xi的一阶敏感性指数Si可表示其对模型输出总方差的直接贡献，即

(11)

总敏感性指数为

(12)

式中V-i——不包括参数xi的其他所有参数方差之和

1.3 模型参数敏感性分析方案

全局敏感性分析试验采用SimLab(Version 2.3)进行分析，CERES-Wheat模型运行借助于RStudio进行模拟运行，具体方案如下：

(1)在SimLab中定义模型输入参数的取值范围及分布形式。本研究中土壤参数取值范围如表1所示，并假设参数在范围内服从均匀分布。

(2)在参数范围内进行随机取样，生成多维参数集。EFAST法认为采样次数大于参数个数65倍的分析结果才有效，本研究共产生1 455组参数。

(3)将生成的参数集写入对应的CERES-Wheat模型文件中，运行模型并整理模拟结果。

(4)将模拟生成的数据整理成SimLab可识别的文本格式，通过SimLab进行分析并得到最终的敏感性分析结果。

2 结果与分析

2.1 冬小麦生产指标敏感性分析

由土壤参数对冬小麦生产指标的一阶敏感性指数和全局敏感性指数分析结果(图1)可知，对于地上生物量平均全局敏感指数最大的前3个参数为SDUL、SLHB和SLRO，分别可解释地上生物量变化方差的70%、57%和48% (图1a)。一阶和全局敏感性分析结果表明，SLHB对CWAM的影响主要通过与其他参数的交互作用实现，而SDUL、SLRO对CWAM的影响主要是其直接效应引起的，与其他参数的交互作用对CWAM的影响较小。在参数设定范围内，CWAM随SDUL、SLHB、SLRO参数值增加呈现增加趋势。除SDUL、SLRO、SLHB、SLNI 4个参数外，其余参数全局敏感性指数均小于0.10，表明对CWAM的影响较小，甚至可忽略不计。与CWAM敏感性结果类似，HWAM和HIAM对SDUL、SLRO、SLHB、SLNI较为敏感(图1b、1c)，在参数值获取时需要重点关注，尤其是在目标土壤存在高异质性时。

图1 土壤参数对冬小麦生产指标的敏感性分析Fig.1 Results of global sensitivity analysis of soil parameters for production indexes of winter wheat

2.2 冬小麦氮素分布指标敏感性分析

土壤参数对冬小麦中氮素分布指标的敏感性结果(图2)表明，SLHB是影响地上生物量中氮素含量最敏感的参数，可解释CNAM变化方差的75%，且主要通过与其他参数的交互作用产生影响(图2a)。在参数设定范围内，CNAM随SLHB参数值的增加而增加。SDUL则是对籽粒中氮素含量最敏感的参数，可解释GNAM变化方差的70%，且其直接效应占主导(图2b)。在参数设定范围内，GNAM与SDUL呈正相关，GNAM随SDUL增加而增加。对地上生物量和籽粒中氮素含量平均一阶敏感性指数大于0.10的参数有SDUL、SLRO、SLHB和SLNI，表明这些参数是影响CNAM和GNAM较为敏感的土壤参数，与其余土壤参数相比，在野外获取这4个参数时需充分考虑土壤异质性的影响。对于根中氮素含量而言，仅有SLNI是引起冬小麦根氮素含量的最关键土壤参数，且RNAD与SLNI呈正相关关系；其余参数全局敏感性指数均小于0.10，表明对RNAD变化影响小于10%(图2c)。对于叶中氮素含量，平均一阶和全局敏感性指数大于0.10的参数均为SLNI和SDUL，且其直接效应占主导，而与其他参数的间接效应对LNAD变化的影响较小(图2d)。在参数设定范围内，LNAD与SLNI和SDUL均呈正相关。

2.3 土壤氮素形态指标敏感性分析

由土壤参数对模型模拟的作物吸收氮量、氮素淋失量及土壤中硝态氮、铵态氮含量的敏感性分析结果(图3)可知，影响作物吸收氮量较敏感的参数有SLRO、SDUL、SLNI和SLHB，且SLRO、SDUL和SLNI的直接效应占主导，与其他参数的间接效应影响较小(图3a)。在参数设定范围内，NUCM与SLRO、SDUL、SLNI均呈正相关关系。对氮素淋失影响最敏感的参数为SLDR，可解释其变化方差的65%，且其直接效应占主导(图3b)。对NLCM，平均一阶敏感性指数大于0.10的参数还有SDUL、SBDM和SLNI。虽然SLCL对NLCM的直接影响较小，仅可解释NLCM变化方差的2%，但其通过与其他参数交互作用所引起的间接影响可解释NLCM变化方差的26%，表明在研究模拟氮素淋失时，SLCL也是需要重点测定的土壤参数之一。在CERES模型中，NLCM与SLDR、SBDM、SLNI均呈现正相关关系，而与SDUL呈负相关关系。对于土壤硝态氮含量，仅对SLNI和SBDM较为敏感，分别能解释NITD变化方差的76%和29%(图3c)。在参数设定范围内，NITD分别随SLNI和SBDM增加而增加。土壤铵态氮含量则对SLNI和SLHB变化较为敏感，分别能解释NHTD变化方差的69%和17%，且均是直接效应占主导(图3d)。在参数设定范围内，NHTD与SLNI呈正相关关系，而与SLHB呈负相关关系。在研究模拟氮素分布时，需要重点考虑土壤异质性对SLNI、SLHB和SBDM的影响。

图3 土壤参数对作物吸收氮量、氮淋失量及土壤中硝态氮、铵态氮含量的敏感性分析Fig.3 Results of global sensitivity analysis of soil parameters for nitrogen uptake and leached, content of total soil nitrate and ammonium nitrogen

2.4 土壤氮素转化指标敏感性分析

由土壤参数对模型模拟涉及土壤氮素转化的氮矿化、硝化、反硝化及氨挥发的一阶敏感性指数和全局敏感性指数(图4)可知，NMNC和NITC均对SLNI最敏感，分别解释其方差变化的89%和79%，而且参数间的交互效应较小，参数的直接影响占主导。其次NMNC和NITC对SBDM也较敏感，而对其余参数的敏感性较差，甚至可忽略不计(图4a、4b)。在参数设定范围内，NMNC与SLNI和SBDM均呈正相关关系，与NMNC类似，NITC随SLNI和SBDM增加而增加。对于氮反硝化量，与NITC的分析结果不同，仅有SLDR对其影响较大，可解释其方差变化的92%，其余参数对其影响很小，全局敏感性指数均不足0.10(图4c)。在参数设定范围内，NDNC与SLDR呈负相关关系，NDNC随SLDR增大而减小。对于氨挥发量直接影响较大的参数有SLNI、SLHB和SLDR，且与SLNI、SLHB呈正相关关系，而与SLDR呈负相关关系。虽然SBDM、SSKS和SRGF对AMLC的直接影响很小，但通过与其他参数的交互作用可分别解释AMLC方差变化的12%、11%和10%(图4d)。

图4 土壤参数对氮素矿化量、硝化量、反硝化量及氨挥发量的敏感性分析Fig.4 Results of global sensitivity analysis of soil parameters for nitrogen mineralization, nitrification, denitrification and ammonia volatilization

3 讨论

扩展傅里叶幅度检验法不仅考虑了单个参数对模型输出的直接影响，而且考虑了参数间的交互作用，解决了局部敏感性分析法仅考虑参数直接效应而忽视参数间间接效应的不足[1]。从全局敏感性分析来看，冬小麦生产指标对土壤水分参数SDUL和SLRO较为敏感。田间持水率反映了土壤持水能力，在生产实践中，常用来作为灌溉上限和计算灌水定额的指标。模型模拟土壤水分对作物生长的影响主要通过影响土壤蒸发蒸腾量来影响水分在土层中的分布，进而影响作物根系吸水及作物生长[12]。SLRO用来估算降雨径流量，反映了冬小麦对降雨的利用效率，其取值与土地利用、坡度、土壤类型以及土壤湿度有关。模型输出的3个生产指标敏感性分析结果存在类似趋势，主要是因为这3个指标存在较强的相关关系。

氮素是制约作物生长和产量形成的主要因素[22]。在以模拟作物产量等生产能力指标为研究重点时，需要重点确定SDUL和SLRO取值的同时，也不能忽视SLNI和SLHB的影响。对于地上生物量和籽粒中的氮素而言，最敏感的参数是田间持水率，而对根和叶中的氮素而言，最敏感的参数则是土壤中总氮含量，这可能是因为冬小麦根、叶生长主要受氮素含量限制，而进行干物质分配时主要受土壤水分含量的影响。氮素施用对作物生长及干物质形成具有良好的作用，可增加株高、千粒质量及干物质量[21,30]。同时庄恒扬等[22]指出，开花期后作物叶、茎中的氮素开始向籽粒中转移，而根系的吸氮能力也随生育期进程而不断下降。此外，在CERES模型中，模型假设叶片中的氮素优先转移、补充营养器官[22]，这可能是造成籽粒中氮素含量对田间持水率和土壤总氮含量均敏感而根和叶中氮素含量仅对土壤总氮含量敏感的原因之一。作物吸收氮量敏感性分析结果与作物地上生物量及籽粒中含氮量的结果类似，这也反映了作物吸收的氮素大部分转移到了籽粒中[21-22,30]。当研究侧重于氮素在作物中的分布情况时，应重点考虑土壤异质性对SDUL、SLRO、SLNI和SLHB的影响，对这些参数的取值需更加慎重。

氮素淋失过程实际上是土壤氮素与降雨、径流相互作用的过程[23,31]，氮素淋失量与土壤导流能力及土壤氮素含量有密切联系，因此在测定土壤参数时应考虑土壤异质性对SLDR、SLRO、SDUL、SBDM和SLNI的影响。在模型模拟过程中，土壤氮素淋失量随SLDR、SLRO、SBDM和SLNI的增加而增加，但随SDUL的增大而减少，这可能是田间持水率增加可减少排水量，从而导致氮素淋失量降低。土壤中的含氮物质在各种微生物及一系列物理、化学过程的作用下会发生形态及状态的转化，氮素矿化是其中的一种生物化学过程，氮素矿化过程会受土壤质地、含氮量、温度和湿度等诸多因素的影响[23]。在本研究中，CERES-Wheat模型土壤氮素矿化主要受土壤氮素含量及土壤容重的影响，其他土壤参数的影响效应很小。

硝态氮和铵态氮是土壤中无机氮的两种存在形式[23,30]，土壤总氮含量是影响硝态氮和铵态氮含量最敏感因素，此外，硝态氮含量还受土壤容重的影响，而铵态氮含量还受土壤酸碱度的影响，而且随土壤碱性增加而减少。氮素硝化反应是影响土壤硝态氮含量的重要因素，对氮素硝化量的敏感性分析结果与硝态氮含量一致，土壤总氮含量和土壤容重是影响其变化的敏感参数。而氮素反硝化量主要受排水比率的影响，这主要由于反硝化作用是反硝化细菌在缺氧条件下还原硝酸盐的过程[32]，而模型中排水比率决定了土壤中水分排出速率及氧气状况。对于氨挥发量而言，除土壤总氮含量对其敏感外，与土壤含水率有关的参数和土壤酸碱度对其也比较敏感，这也与氨气易融于水及其溶液呈碱性易与酸性溶液中和的特性相一致。在研究土壤氮素转化及不同形态氮素分布时，需要同时考虑土壤异质性对SLDR、SLRO、SDUL、SLLL、SRGF、SSKS、SBDM、SLNI和SLHB的影响，而其他参数对土壤氮素含量影响较小，可适度降低土壤参数获取难度。

4 结论

(1)对于冬小麦生产指标，田间持水率、径流曲线数、土壤酸碱度及土壤总氮含量是影响冬小麦地上生物量和产量较为敏感的土壤参数，需要考虑土壤异质性的影响，而其他参数的影响较小，可将原来需测参数数量降低至4个，且冬小麦生产指标与这4个参数均呈正相关关系，降低了参数获取难度。

(2)对于冬小麦氮素分布指标，田间持水率和径流曲线数是影响地上生物量及冬小麦籽粒中氮素含量较为敏感的土壤参数，而土壤总氮含量则是影响冬小麦根和叶中氮素含量最敏感的参数。在研究作物氮素分布模拟时，需重点测量土壤田间持水率、径流曲线数、土壤酸碱度及土壤总氮含量4个参数。

(3)对于氮素转化与分布指标，田间持水率和径流曲线数是影响作物从土壤中吸收氮素量较为敏感的参数，而排水比率是影响氮素淋失量、反硝化量最敏感的参数，土壤总氮含量则是影响土壤中硝态氮、铵态氮、氮矿化、氮硝化量最敏感的参数，土壤总氮含量、土壤酸碱度及排水比率则是影响氮素氨挥发量较为敏感的土壤参数。在模拟氮素转化与分布时，需重点考虑土壤异质性对排水比率、径流曲线数、凋萎含水率、饱和含水率、饱和导水率、根系影响因子、土壤容重、总氮含量及土壤酸碱度的影响，原需测量的15个土壤参数可最大简化为9个重点测量参数。