崔 颖 蔺宏宏 谢 云 刘素红
北京师范大学地表过程与资源生态国家重点实验室, 北京 100875
水是作物生长的重要因子, 提高作物水分利用效率, 既可降低农作区水资源消耗, 弥补水资源短缺, 并可有效提高作物产量。作物生长模型是研究水分利用效率和作物产量的重要工具, 能在单点尺度上定量描述给定作物在一定环境下整个生长周期内的生长发育和产量形成过程[1-2], 进而建立起作物生长对土壤及气象环境变化的响应系统[3], 为精准农业生产提供切实可行的研究方法。目前应用比较广的作物模型主要有AquaCrop、WOFOST、CERES等[4], AquaCrop 虽然比后2 个模型开发得晚, 但其核心部分的作物需水量计算方法早在20 世纪70 年代就已被广泛应用[5]。这些模型都是以日为步长模拟作物的生长过程与产量形成, 但模拟过程略有差异。AquaCrop 模型重点考虑了水分胁迫, 以水分利用效率为基础模拟产量[6]。WOFOST 和CERES 模型都是先模拟环境适宜条件下的最大产量, 然后考虑热量、水分和氮磷等肥力胁迫。不同的是WOFOST对呼吸作用考虑更加全面, 且与AquaCrop 模型类似,首先模拟生物量, 然后利用收获指数与生物量相乘得到作物产量[7-8]。CERES 模型则是模拟作物的穗数、穗粒数和穗粒重后计算产量[9]。旱地农业受水分胁迫为主, AquaCrop 模型强调水分胁迫机制、模拟过程相对简单、且所需输入参数较少, 使其得到广泛应用[10]。
目前已有多位国内外学者应用AquaCrop 模型在不同地区进行了作物生长动态模拟并开展了模型适应性评估研究: 刘兴冉等[11]将AquaCrop 模型应用于石家庄市栾城区夏玉米水分研究中, 证明AquaCrop 模型能够较好地模拟夏玉米的产量、生物量和冠层发育过程以及表层土壤水储量的动态变化;刘琦等[12]在山西省寿阳县进行了覆膜和露地春玉米种植对比试验, 以验证AquaCrop 模型对土壤含水率、农田蒸散和冠层覆盖度的模拟精度, 认为该模型能较好地模拟晋中地区旱作覆膜春玉米的耗水、生长和产量形成过程; Iqbal 等[13]应用AquaCrop 模型对石家庄市栾城区冬小麦的生长过程进行了模拟,分析了不同程度水分胁迫条件下的模拟精度, 发现在水分严重亏缺时模拟结果与实际产量出入较大,而当参数率定后则能有效提高产量模拟精度。基于此, AquaCrop 模型在应用前需对强敏感性参数进行率定。模型参数敏感性分析是筛选强敏感性参数的主要方法, 其中OAT (one factor at a time)[14]方法是最快速有效的敏感性方法之一, 即每次只改变1 个参数并根据参数改变量与模拟结果变化量的相关程度进行筛选。已有部分研究对其他作物生长模型参数敏感性进行探究, 但关于AquaCrop 模型参数敏感性的研究仍有待进一步补充和完善: 刘刚等[15]使用OAT 方法在山东省禹城市对ALMANAC 模型参数进行了敏感性分析, 为后续应用模型模拟黄淮海平原地区冬小麦和夏玉米的生长提供了科学依据; 宋明丹等[16]采用Morris 和EFAST 的敏感性分析方法,分析了陕西省杨凌区CERES-Wheat 模型中对开花期、成熟期、产量、地上生物量4 个输出变量的敏感性参数, 发现2 种方法筛选出的强敏感性的作物参数基本一致, 其中Morris 法所需的计算量更小,操作性更强, 而EFAST 法精确度更高, 可定量分析各个参数的敏感性, 在进行参数的不确定性分析方面更有优势。可见, 简单有效的敏感性分析方法即可实现强敏感性参数筛选的目的, 从而有针对性的进行参数率定, 这对于提高AquaCrop 模型的模拟精度和区域适用能力而言十分必要。
东北黑土区是世界重要黑土带之一, 土壤肥沃,气候适宜, 主要出产玉米和大豆等作物, 是国家粮食安全战略性保障基地[17-18]。特定的土壤和自然地理条件使得其雨养农业生产模式有别于其他粮食主产区, 大气降水成为该地作物产量的主要影响因子之一[19]。为能应用AquaCrop 模型模拟东北黑土区作物产量的水响应机制, 弥补该模型在东北地区应用研究的缺失, 本研究采用OAT 敏感性分析方法在分析模型输入参数敏感性的基础上, 基于作物试验观测数据对敏感参数进行率定, 然后利用大田观测资料验证参数率定后的AquaCrop 模型在东北黑土区的适用性, 以便通过控制和减小环境变量效应,调整和优化作物种植措施, 为该地区农业高产可持续发展提供科学依据。
研究区位于黑龙江省嫩江县鹤山农场, 地处小兴安岭西南麓、由小兴安岭向松嫩平原过渡的鹤北流域,地理位置(48°43′—49°03′N, 124°56′—126°21′E)。该区海拔310~390 m, 地形起伏较小, 有坡长坡缓的特征[20], 坡度大多在3°~6°之间, 坡长多为2000~4000 m。属寒温带大陆性季风气候, 冬季寒冷干燥, 夏季温热多雨, 年均温约0.4oC; 多年平均年降水量约534 mm, 6 月至8 月降水量约占全年的66.6%[21]。区内土壤以黑土和草甸土为主, 属于东北典型黑土区, 因有机质含量丰富而呈黑色。黑土土壤剖面表层为20~40 cm 的黑土层, 土壤质地比较黏重, 土壤容重为1.39 g cm–3, 田间持水量为0.28 g g–1, 凋萎湿度为0.18 g g–1。向下为过渡层逐渐到达母质层, 包括黄黏土、黄沙土或砾石等。研究区开垦于20 世纪40 年代, 流域内农业用地占总土地面积的62.4%, 采用大型机械化作业和集中农业管理方式, 主要农作物为大豆和玉米每年轮作, 一般于每年5 月初播种至当年9 月底收获。玉米生长期为150 d, 大豆生长期为130 d。
AquaCrop 模型是由联合国粮食及农业组织(Food and Agriculture Organization of the United Nations, FAO)第33 号灌排文件Yield response to water[5]中的作物-水分响应方程演变而来。该方程假设, 作物产量(yield,Y) 的高低是对作物蒸散量(evapotranspiration,ET)大小的响应结果:
式中,xY、aY分别为作物最大产量(kg m–2)和实际产量(kg m–2)。ExT 、EaT 分别为作物潜在蒸散量(mm)和实际蒸散量(mm)。潜在蒸散量指充分供水条件下的蒸散量, 首先根据 FAO 第 56 号文件[22]中的Penman-Monteith 公式计算出参考蒸散量ET0[23], 然后再乘以作物系数得到。实际蒸散量是指实际供水情况下的蒸散量, 会受水分供给不足导致的土壤水分胁迫影响。ky为作物产量对土壤水分亏缺响应的敏感系数, 随作物生育期变化。建立的AquaCrop 作物生长模型以日为模拟步长, 通过输入气候、作物、土壤, 田间管理如灌溉、施肥及除草等数据, 可以模拟作物的光合、呼吸、蒸腾等过程, 最终输出冠层覆盖度、生物量、产量等结果。
作物产量是作物成熟时的生物量(biomass, Bx)与收获指数(harvest index, HI)相乘得到。作物生物量则是由归一化的水分生产力即单位耗水量累积的生物量(normalized water productivity, WP*)乘以作物蒸腾量与参考蒸散量的比值得到:
式中,Bi为日地上生物量(t hm–2); Tri为日蒸腾量(mm), ETo,i为日参考蒸散量(mm)。WP*为归一化水分生产力(gram m–2), 是将作物水分生产力除以标准蒸散条件下的作物蒸散量。实际计算时, 可根据大气实际CO2浓度、C3或C4作物类型、以及作物生长期和产量成熟期进行参数调整。
利用Penman-Monteith 公式计算出参考蒸散量后, 根据冠层覆盖度最大时的作物蒸腾系数和实际生长过程的冠层覆盖度(canopy coverage, CC), 可以将作物生长过程中的作物蒸腾量(transpiration,Tr)从蒸散量中分离出来, 剔除了土壤蒸发的非生产性消耗性用水(evaporation, E)干扰。此外随着冠层扩张,因遮阴和冠层对空气平流的影响都会影响作物蒸腾,需要对CC 进行修订, 用CC*表示:
式中, CC*为考虑了冠层遮阴和空气流动影响后的冠层覆盖度(%); Kcbx为充分供水条件下参考作物蒸散的蒸腾系数。实际作物蒸腾受土壤水分胁迫影响,分别表现为缺水情况下的气孔导度胁迫, 用胁迫系数Kssto表示, 以及多水情况下的土壤通气性胁迫,用胁迫系数Ksaer表示。二者都是0~1 的无量纲参数,表示胁迫程度, 值为1 时无胁迫。
冠层覆盖度随作物生长发生变化, 直接影响作物的光合作用。模型采用指数形式表达CC 的增衰变化[24], 包括指数增长期和平稳增长期(式4), 衰落期(式5)。
式中, CC 为冠层覆盖度(%);t为从出苗开始累积的时间; CC0为初始冠层覆盖度(%), 一般取90%出苗时的平均幼苗覆盖度; CCx为冠层覆盖度达到最大时的值(%); CGC 为冠层覆盖度增加速率(canopy growth coefficient), 表示单位生长度日冠层覆盖度的增加量(%, GDD–1); CDC 为冠层覆盖度衰减速率(canopy decline coefficient), 表示单位生长度日冠层覆盖度的减少量(%, GDD–1)。如果受到土壤水分胁迫, 需要对CGC 和CDC 进行修正, 修正公式如式(6)~(7)。
式中, CGCadj和CDCadj分别表示受水分胁迫影响的CGC 和CDC, Ksexp为水分对冠层覆盖增加的胁迫系数; Kssen为水分对冠层覆盖减少的胁迫系数。此外, 收获指数也会受到水分胁迫的影响, 缺水时利用Kssto调整; 滞水时利用调整。
植物水分依靠根系吸水, 因此模型考虑了根区中的根系生长, 及其对水分胁迫的响应。给定土壤水分上限和下限, 当根区水分消耗达到或超过上限时,会影响冠层覆盖增加, 但不会影响根系加深, 因为水分胁迫对冠层的影响强于对根系的影响。如果水分胁迫低于叶片气孔关闭的阈值, 根系受到胁迫, 不会加深, 否则可以加深。根系增长速率计算公式如下:
式中,Z为t时刻的有效根深(m);Zini为根系初始深度(m);Zx为最大有效根深(m);t0为90%出苗的时间;tx为根系达到最大有效根深的时间;t为根系生长时间。
AquaCrop 模型涉及气候、土壤、作物、管理等因素, 同时考虑了充分供水和水分胁迫对作物生长发育和产量形成的影响。但其中的控制方程和参数都是基于试验观测获得的经验方程和不同地区的参数。考虑到不同地区之间气候、土壤、作物品种、管理等的差异, 将模型应用于某地区时, 应进行参数的率定和模型的验证, 使其适用于当地条件。
研究收集了所在地鹤山农场气象站观测数据,包括日最高/低气温(°C)、日降水量(mm)、太阳辐射量(MJ m–2d–1)、水汽压(kPa), 并利用FAO 发布的ET0Calculator 工具计算ET0(mm d–1), 建立气象数据库。通过测定各试验点土壤剖面的容重、田间持水量和凋萎湿度等土壤属性数据, 建立了东北黑土区土壤特性数据库, 部分参数如表1 所示。
表1 试验点的土壤物理参数表Table 1 Soil physical parameters of the experimental locations
参数率定所用数据来自2016—2017 年位于黑龙江省嫩江县鹤山农场的北京师范大学九三水土保持试验站作物试验。种植的玉米品种为德美亚1 号, 大豆品种为1467。分别以鹤山农场玉米实际平均种植密度9.25 万株 hm–2和当地大田大豆种植密度45 万株 hm–2为标准密度, 按标准密度的60%、80%、100%、120%、140%共5 个密度, 3 个重复, 2 种作物各15 个小区种植, 单个小区面积为5 m × 8 m, 小区间隔0.5 m。试验播种时按当地种植方式每公顷施肥量为0.4 kg, 3 种底肥尿素、磷酸二铵和硫酸钾的配比为1.0∶1.3∶0.3, 并依据管理标准按时中耕、除草以保证作物的正常生长。每2d 观测1 次作物生育及土壤水分动态, 记录各个小区作物出苗期、开花期、成熟期等关键生育期的出现日期。玉米成熟后(约150 d), 采生长状况接近平均水平的植株12 株风干后将所有植株籽粒混合均匀后采样放入85℃烘箱烘干至恒重后称重。大豆成熟后(约130 d), 各小区采2 m × 2 m 的大豆样本地上部分用恒温鼓风烘干箱105℃杀青2 h 后转80℃烘干至恒重后称重。
模型验证使用2011—2018 年鹤山农场选取的16 块农田监测点观测的数据。在生长季内(5 月至9月)每10 d 观测1 次监测点的作物生长动态, 观测指标包括株高、地上生物量和收获指数等。对试验小区和农田实测作物数据进行整理与分析, 建立鹤北小流域作物特征参数数据库。
AquaCrop 模型部分参数不随地理位置、种植时间和管理措施变化而被称为保守参数, 无需调整[25]。为了认识非保守参数中对模拟结果的影响程度, 采用敏感性分析的方法评估模型参数对模拟结果的影响程度[26]。本文采用OAT 法, 参考FAO 作物参数手册[27]及实际作物特性确定拟进行率定的模型参数(表2), 使用2011—2018 年大田数据, 在AquaCrop模型中逐一改变某特征参数为标注范围极值或默认值±10%及±20%的调整值, 保持其余变量为原数据集值不变, 对比作物产量的模型模拟结果, 分析产量对不同变量的敏感性。敏感性分析结果由相对敏感度(relative sensitivity, RS)[28]表示(式9), 其值经过了标准化处理,可以进行参数之间敏感性大小的比较[29]。
式中, RS 为相对敏感度,Y为模拟的作物产量,X为变量, ΔX为变量的调整值。按RS 值由大到小可将变量对产量影响的敏感程度划分为: 极敏感(RS>0.6)、一般敏感(0.6>RS>0.1)和低敏感(RS<0.1)。
本研究采用2016—2017 年不同种植密度下的作物田间试验观测数据率定模型参数, 对 RS 较低的参数使用默认值, 仅对 RS 值较高的参数进行率定。
表2 AquaCrop 模型待敏感性分析参数介绍Table 2 Introduction of parameters to be sensitive analysis in AquaCrop model
基于建立的气候、作物、土壤数据库和率定后的参数, 采用2011—2018 年大田监测点实测的18组玉米和16 组大豆产量对模型进行验证。模型模拟精度评价采用决定系数(coefficient of determination,R2)、均方根误差(root mean square error, RMSE)、标准均方根误差(normalized root mean square error,NRMSE)和模拟效率(model efficiency, ME)作为精度评定指标, 评估预测结果与实测结果的差距, 计算公式如式(10)~(13)所示,R2和ME 越接近1、RMSE和NRMSE 越接近0 说明预测精度越高。
式中,n为总观测数;Oi和Si分别为第i次观测的实测值和预测值, t hm–2;和分别为平均实测值和平均预测值, t hm–2。
表3 AquaCrop 模型部分参数相对敏感度Table 3 Relative sensitivity of parameters in AquaCrop model
应用AquaCrop 模型对东北黑土区主要作物生长过程进行模拟时, 两种作物的参数相对敏感度计算结果如表3 所示(相对敏感度为0 的参数表中已略)。对不同作物而言, 同一参数对产量的敏感性有较大差异。玉米产量对以下参数表现为极强的敏感性: 冠层形成和枯萎前的作物系数(crop coefficient before canopy formation and senescence, KcTr,x)、WP*、参考收获指数(reference harvest index, HI0)、CDC、最大有效根深(maximum effective rooting depth,Zx)。大豆产量的极敏感性参数较少, 包括限制冠层伸展的水分胁迫系数曲线的形状因子(shape factor for water stress coefficient for canopy expansion,Pexshp)、KcTr,x、根区根系伸展曲线的形状因子(shape factor describing root zone expansion, Rexshp)。部分参数的相对敏感度对两种作物而言基本一致, 如限制冠层伸展的土壤水分消耗上限阈值(soil water depletion threshold for canopy expansion-upper threshold, Pexp-up)、因衰老、氮元素亏缺导致的作物系数下降速率(decline of crop coefficient as a result of ageing, nitrogen deficiency, DeKcTr,x)、产量形成期的归一化水分生产力(water productivity normalized for ET0and CO2during yield formation, WP*yf)、CCx、CGC 和KcTr,x。RS 结果描述了参数变化对最终产量值的影响, 在一定程度上忽略了各参数间的相互影响及不同生长发育阶段的影响程度变化, 如CGC 的减小会导致最大冠层覆盖日的延迟到来; CGC 降低或最大冠层覆盖度过大均会使HI0降低等。尽管该方法得到的敏感性分析结果存在一定的不确定性, 但已达到了筛选敏感性参数的目的, 可服务于参数率定。
基于敏感性分析结果, 玉米及大豆的敏感性参数差异较大, 需对两种作物分别进行参数率定。玉米作物主要率定了HI0、Zx、WP*等参数, 率定值与默认值最大差距 230%, 率定后保持不变的有Pexp-up 和WP*yf。大豆作物主要率定了CGC、CDC、Zx、WP*yf 等参数, 率定值与默认值最大差距166%,率定后保持不变的有KcTr,x和Rexshp。表4 展示了AquaCrop 模型模拟东北黑土区玉米及大豆生长的率定后参数。对水分胁迫参数调整较大, 尤其是Pexp-up、限制冠层伸展的土壤水分消耗下限阈值(soil water depletion threshold for canopy expansion-lower threshold, Pexp-lw)、Pexshp 等限制冠层生长的水分胁迫参数在率定后使模拟结果精确度有了很大提高, 更适用于研究区土壤及气候大环境。率定后的AquaCrop 模型参数为精确模拟东北黑土区玉米和大豆产量奠定了基础, 也为后续该地区作物研究提供了参考。
表4 AquaCrop 模型模拟东北黑土区典型作物生长的参数率定结果Table 4 Calibrated parameters for simulating black soil area in Northeast China
利用大田监测点的18 组玉米和16 组大豆的作物实测产量, 对率定参数前后的实测产量与模型预测产量进行对比, 结果见图1。图中可见: (1)参数率定前后, 模型预测精度提升显著。玉米的预测产量回归系数由0.34 提升至0.89, 回归方程截距由7.50 t hm–2降至1.42 t hm–2,R2由0.4797 提升至0.7745,RMSE 由3.271 t hm–2降至1.076 t hm–2, NRMSE 由0.293 降至0.097, ME 由1.333 降至0.747; 大豆的预测产量回归系数由0.80 提升至0.88, 回归方程截距由0.34 t hm–2降至0.29 t hm–2,R2由0.5811 提升至0.7789, RMSE 由0.941 t hm–2降至0.299 t hm–2,NRMSE 由0.561 降至0.178, ME 由1.674 降至0.747。参数率定修正了预测值在趋势上的偏离和系统性误差, 提高了模型的适用能力。(2)玉米和大豆的实测产量与率定后预测产量的R2分别为0.775 和0.779,RMSE 分别为1.076 t hm–2和0.299 t hm–2, NRMSE分别为0.097 和0.178, ME 分别为0.747 和0.730, 均表明率定后的模型具有较好的预测效果。从玉米和大豆的预测结果可见, 针对不同作物的不同产量状态, AquaCrop 模型的预测效果存在差异。例如玉米低产时的模拟精度较低, 回归方程截距较大, 但相较于公顷产量而言误差在可接受范围内; 而大豆无论是高产还是低产时, 其模拟精度基本保持无偏。
相对敏感性分析结果发现, 对玉米和大豆产量比较敏感的参数可分为3 类: 一是与冠层相关的参数, 如Pexshp 和CDC; 二是与根系相关的参数, 如Rexshp 和Zx; 三是与作物遗传特征相关的参数, 如KcTr,x、WP*和HI0。叶片既是光合作用的场所又是植物蒸腾的器官, 植物所需水分来自根系吸水, 作物系数水分生产力决定了作物需水量和水分生产效率,收获指数则决定了生物量中经济产量的高低, 因此上述参数的微小变化会使产量发生较大变化。然而玉米和大豆却对3 类参数的敏感度有所差异。在与冠层相关的参数中, 玉米对冠层衰减系数更为敏感,而大豆则对影响冠层生长的水分胁迫参数更为敏感。从根系参数看, 由于大豆根系一般比玉米根系分布浅, 大豆产量对根系生长的形状参数更为敏感,而玉米产量则对最大有效根深更为敏感。由于玉米需水量大且对水分变化敏感, 其对影响作物需水量大小的作物系数和水分生产力更为敏感, 大豆则表现为一般敏感。两种作物都对影响经济产量的收获指数非常敏感。综上, 敏感性分析不仅反映出模型模拟作物产量的关键影响过程, 还揭示了对不同作物的差异, 具有很好的适用性。实际应用中, 应确保这些敏感参数输入的准确性。
本文采用的OAT 方法属于局部敏感性分析方法的一种, 计算效率较高, 分析的是单个参数变化对模型输出变量的影响, 一定程度上忽略了参数间的相关性。与之相比, 以扩展傅里叶幅度检验法(extended fourier amplitude sensitivity test, EFAST)为代表的全局敏感性分析方法[30-31]考虑了参数间交叉作用和作物在不同生长阶段的参数变化, 更为全面。已有研究采用EFAST 方法对AquaCrop 模型进行了敏感性分析[32], 结果表明在北京地区对冬小麦产量的强敏感性参数包括: Pexp-up、HI0、KcTr,x、CDC、WP*。与本文OAT 方法得出的玉米产量敏感参数具有很高的一致性, 与大豆产量敏感参数也存在部分重叠, 也说明采用 OAT 方法分析的AquaCrop 模型敏感性参数是合理的。AquaCrop 模型在松嫩平原预测春小麦产量时[33], 模拟产量与实测产量的R2在0.887~0.969 之间, RMSE 为0.11 t hm–2, 模拟效果好于本研究, 说明该模型在东北黑土区具有良好的适用性和发展前景, 应选择更多的地区验证玉米和大豆的产量模拟结果。
AquaCrop 模型中许多假设条件均基于理想情况, 诸如土壤层分布均匀、无病虫害、无杂草、无养分限制等, 未考虑自然或人为影响造成的突发情况, 导致模拟结果与实际产量存在偏差。对敏感性高的参数进行率定可以提高模型的模拟精度。作物生长机理模型在单点上已经有了较高的模拟精度, 将其扩展到更大的区域范围进行连续模拟时, 会由于单点数据的随机性、模拟过程的复杂性以及变量的自相关等, 导致模拟结果具有不确定性。如何融合单点观测和多源观测如遥感数据, 是将作物模型扩展到更大空间尺度, 并提高模拟精度的发展方向。
遥感卫星数据[34]能够快速获取大尺度地表信息而在农业监测与作物估产中有巨大应用空间。集合卡尔曼滤波(ensemble kalman filtering, EnKF)[35]等数据同化算法能以实时观测数据对先验的模型数据结果进行再分析, 改进了系统的可靠性, 从而为作物产量预测提供更大参考价值。研究还需结合多源数据和同化方法进一步优化AquaCrop 模型在东北黑土区的作物产量模拟结果, 控制和减小模型模拟不确定性对产量预测造成的误差。
AquaCrop 模型的参数敏感性分析结果表明水分胁迫等环境变量对作物产量有较大影响, 主要集中在以下几个参数: Pexshp、Rexshp、KcTr,x、WP*、HI0、CDC 和Zx。不同作物的参数敏感性差异较大。精准的田间观测可实现模型参数的率定, 有效弥补了模型在客观生物理化意义上的缺失, 降低了作物模型中敏感性参数引起的误差, 可提高模型的预测精度。率定的参数越丰富, 模型模拟效果越好, 率定参数的个数与模拟精度呈正相关。玉米和大豆的实测产量与参数率定后AquaCrop 模型模拟的产量之间差异较小, 其R2分别为0.775 和0.779, RMSE 分别为1.076 t hm–2和0.299 t hm–2, NRMSE 分别为0.097 和0.178, ME 分别为0.747 和0.730。表明参数率定后的AquaCrop 模型可较准确地估算东北黑土区玉米和大豆产量, 适用于作物品种和环境条件相似的作物产量模拟与预测。