弓开元,何亮,邬定荣,吕昌河,李俊,周文彬,杜军,于强,7
(1 西北农林科技大学资源环境学院,陕西杨凌 712100;2 国家气象中心,北京 100081;3 中国气象科学研究院,北京 100081;4 中国科学院地理科学与资源研究所,北京 100101; 5 中国农业科学院作物科学研究所,北京 100081;6 中国气象局成都高原气象研究所,成都 610071; 7 西北农林科技大学黄土高原土壤侵蚀与旱地农业国家重点实验室,陕西杨凌 712100)
【研究意义】气候变化对全球粮食安全产生了深远的影响[1]。气候变化主要体现在降水、温度、太阳辐射等气候因子的变化,对作物的生长发育过程产生影响,最终影响作物产量。对不同作物和不同区域,其影响程度差异较大[2]。因此,明确气候变化对特定地区和特定作物的影响,认识区域作物生产潜力,对于应对未来气候变化、制定适应性措施具有重要意义。青藏高原是世界上地势最高的地区,具有独特而复杂的高原气候,是气候变化的敏感区和启动带,也是全球气候变化的驱动器和放大器[3-4]。在作物生长方面,相比于其他种植区,青藏高原由于海拔原因具有最高的太阳辐射、较低的温度和最低的CO2分压。因此,作物的生长发育特别是光合作用对气候变化的响应及敏感性与低海拔的平原地区有较大差异[5]。近30 年来,青藏高原气候变化显著,总体表现为气温上升,降水增加,整体为由干向湿变化[6-7]。青稞作为青藏高原第一大粮食作物,在2002—2004 年,其播种面积和总产量分别占粮食总产的43%和38%[8]。因此,研究青藏高原青稞的生产潜力、产量差及其对气候变化的响应,对促进青藏高原农业发展和粮食安全具有十分重要的科学价值和实际意义。【前人研究进展】在20 世纪70、80 年代,有关作物生产潜力的研究开始进入快速发展期,DE DATTA[9]首次提出产量差概念。2009 年LOBELL 进一步完善产量差概念,并把产量差定义为特定时空下潜在产量与农户实际产量的差值[10],杨晓光等[11]2014年总结了产量差的定义,并根据潜在产量、可获得产量、农户潜在产量和农户实际产量4 个产量水平将产量差分为3 个层次。不同产量水平中,潜在产量被认为是作物在不受水分和养分限制、病虫害良好控制条件下获得的产量,主要取决于太阳辐射和温度[12]。潜在产量的计算方法较多,1978 年由DE WIT 提出了一套作物生产潜力的计算模型[13],被称之为“FAO 生产力计算模型”,得到广泛使用。如BATTISTI 等[14]用该模型计算巴西大豆的潜在产量和产量差。20 世纪80 年代,通过计算机技术开发的各种作物模型(EPIC、APSIM、CERES 等)也得到广泛应用,并成为现今研究生产潜力和产量差的主要方法。如WU等[15]借助WOFOST 模型计算华北平原夏玉米的潜在产量,ESPE 等[16]借助ORYZA 模型得出了美国水稻的潜在产量和产量差;也有学者通过模型对华北平原冬小麦和东北春玉米等作物在区域上的潜在产量和产量差进行研究[17-18]。国内生产潜力的研究主要集中在东北、黄淮海和南方等主要粮食生产区,虽然有研究表明,青藏高原气候变化对粮食生产产生了巨大影响[3],如青藏高原气候变暖,导致雅鲁藏布江和印度河流域农业灌溉用水减少,粮食安全受到威胁[19]。但对青藏高原本身特别是青稞潜在产量和产量差的研究还鲜见报道。有学者通过经验公式的方法,在站点尺度上结合气象数据对过去40 年西藏自治区和过去10 年青海省的青稞光温生产潜力进行计算,得出西藏自治区青稞光温生产潜力呈上升趋势的结论[20-21];此外,还有学者探究青藏高原过去50 年的气候变化对气候生产潜力的影响情况[22],发现在气候生产潜力升高的基础上中部地区增加更为明显。【本研究切入点】现有的青藏高原生产潜力研究多采用气候因素为主的经验公式法,而该方法对作物生长因素考虑欠缺,且精确度和准确度低,普适性较差[23]。此外,青藏高原地理位置特殊、分布广、海拔差异大,直接以青藏高原或行政区域内的生产潜力计算易造成较大误差。目前,有关青藏高原青稞产量差仍没有系统研究。本研究借助DSSAT-CERES- barley 模型,对青藏高原主要耕作区域站点的光温生产潜力及生育期长度进行模拟。DSSAT 模型(decision support system for agrotechnology transfer)是目前应用最广泛的模型之一,其中的DSSAT-CERES-barley 模型[24]可模拟不同水供应环境下大麦的生长发育情况。【拟解决的关键问题】通过分析青藏高原1977—2017 年间气候因子变化,利用农业气象观测资料校正和检验模型,然后模拟青稞光温生产潜力以及生育期长度的变化趋势和幅度,准确量化不同研究站点青稞产量对气候变化的敏感性。同时结合统计产量进行产量差计算分析,以揭示青藏高原过去40 年青稞产量差的时空分布特征及其受气候变化的影响,量化不同产量水平的差距,明确产量提升空间。
1.1.1 研究站点 青藏高原人口和耕地分布相对集中。选择青稞种植最为广泛的2 个区域进行研究。第一个区域为西藏自治区雅鲁藏布江河谷地区[25];第二个区域为青海省东部环湖、海南州地区和甘肃省西南部甘南藏区[26]。根据2015 年中国土地利用现状遥感监测数据,计算得出研究区域内耕地面积占青藏高原总耕地面积的55.13%,数据来源于中国科学院资源环境科学数据中心(http://www.resdc.cn)[27],并选择这2 个区域中具有青稞农气观测资料的7 个气象站点作为研究站点。这7 个站点的地理位置和气象条件分别见表1。
表1 研究站点基本信息表Table 1 Basic information of research stations
1.1.2 数据源 气象数据来源于国家气象信息中心的中国地面气候资料日值数据集,主要包括最高温度、最低气温、平均气温、日照时数和降水量等;土壤数据为1 km 分辨率的面向陆面过程模型的中国土壤水利参数数据集[28];1987—2017 年的青稞统计产量数据来源于统计年鉴[29-30];7 个研究站点的青稞农业气象观测数据来源于国家气象信息中心的农作物生育状况观测记录年报表,表2 所示分为校准集和验证集,分别进行品种参数的校准和验证;农田管理数据来源于实地考察及文献资料[31-33]。
表2 研究站点品种参数及DSSAT-CERES-barley 模型校准和验证数据集Table 2 Description of values for variety parameters in research stations and dataset of calibration and validation for the DSSAT-CERES-barley model
1.2.1 模型模拟 选用表2的7个研究站点的农业气象观测数据对青稞参数分别进行调参校准。在模型中分别输入7 个站点对应年份的土壤数据、气象数据和农田管理措施数据,管理措施未考虑施肥(氮不受限制),其中站点播种期、出苗期、开花期、成熟期、产量参考农业气象观测数据,水分管理措施(灌溉时期、灌溉量和灌溉方式)根据文献资料及当地实际情况[31-33]。完成参数校准后,输入每个站点的品种参数。山南站只有1997 年农业气象观测数据,由图2、表1可知与拉萨站青稞生长季气候无显著差异且地理位置接近,因此使用拉萨站品种参数进行模拟,根据农业气象数据计算每个站点的平均播期(表1),假设播期和品种不变的条件下,在模型中输入土壤数据和1977—2017 年的气象数据,模拟1977—2017 年站点光温生产潜力及生育期长度。
1.2.2 模型评估 为了验证模型的准确性,选用平均绝对预测误差(MAE)、均方根误差(RMSE)、归一化均方根误差(NRMSE)、绝对系数(R2)4 种指标,计算方法如下:
式中,Obsi为实际观测数据,Simi为模拟数,为实际观测数据的平均值,为模拟数据平均值,n为数值个数。
1.3.1 变化趋势分析 由气象站农气资料可知青藏高原青稞生长季为3—8 月,其平均气温、最高气温、最低气温、降水量和日照时长计算方法通过逐日数据进行平均或求和。青稞生长季≥0℃有效积温通过逐日平均气温数据进行求和。青稞生长季的太阳辐射量通过日照时长计算得到,计算方法通过 ÅngstrÖm- Prescott 公式计算[34]:
式中,Rs为太阳辐射量,n为实际日照持续时间,N为最大可能的日照时数,Ra为天顶太阳辐射,由地理纬度等计算而得,as和bs为参数,参考FAO-56,本研究中as=0.25,bs=0.50。
研究使用Mann- Kendall 非参数秩次相关检验 法[35-36]和 Sen’s 坡度估计法[37-38]对过去40 年时间尺度下的降水、温度和太阳辐射等气象因子趋势进行分析和检验,上述2 种方法常结合使用以揭示时间序列变量变化的幅度和趋势[39-40]。同时对站点光温生产潜力和光温生产潜力下青稞生育期长度模拟进行变化趋势分析。由于Mann- Kendall 非参数秩次相关检验法(M-K 检验)的变量可不完全服从正态分布,因此被广泛应用于水分和气象学的趋势分析。本研究基于M-K 和Sen’s 的计算公式[41],通过R 编制计算程序来分析变化趋势及幅度。
1.3.2 相关性和敏感性分析 本研究使用一阶逐差法,该方法能最大程度减少非气象因素对生育期和光温生产潜力的影响[42],公式为:
式中,ΔXi表示在第i年的X的第一个差异,Xi表示在第t年的时间序列X的值,并且Xi-1是第i-1 年的值。此外,逐步回归使用R“base”包,相关系数计算使用R“Hmisc”包。
1.3.3 产量差计算 本研究对绝对产量差和相对产量差进行了计算,相对产量差的计算方法如下:
式中,PY为光温生产潜力,Y为实际产量,GPYY为光温生产潜力与实际产量之间的产量差。
如图1 所示,不同研究站点之间气候差异显著。其中,雅鲁藏布江河谷地区山南、日喀则和拉萨3 个站点平均海拔在3 500 m 以上,空气稀薄,多为晴天,气温和日照时长都显著高于其他站点,生长季内辐射量可达3 500—4 000 MJ·m-2。林芝站相比于这3个站点海拔略低,在3 000 m 左右,但纬度基本相同,因此平均温度和最高、最低气温均显著高于青海东部和甘肃南部地区站点。由于林芝站生长季内降水量显著高于其他站点,致使生长季内的辐射量最低。而海拔较低的青藏高原东北部的甘南州、门源和贵南之间辐射、积温和降水量都有显著性差异(图1),造成了不同站点之间光温生产潜力和生育期长度的不同。
青藏高原7个站点在1977—2017年的气候倾向率表明(表3),青稞生长季的平均气温、最高气温、最低气温和≥0°C 有效积温显著上升,生长季内平均气温的上升幅度为0.26—0.59 ℃·(10a)-1,最高气温的上升幅度为0.28—0.56 ℃·(10a)-1,生长季平均最低气温的变化幅度最大,尤其是门源和拉萨站的上升幅度超过0.70 ℃·(10a)-1。同时,气候变暖的趋势也导致青稞生长季内有效积温显著上升,上升幅度达到42.52—90.64 ℃·d·(10a)-1。虽然青稞为喜凉作物,≥0℃有效积温1 200—1 500 ℃·d 就可满足其正常发育,但气候变化导致的青藏高原青稞生长季变暖也必然会对青稞生产造成影响。其他气象因子中,生长季内降水量基本呈上升趋势,贵南和拉萨站达到显著性水平,上升幅度为29.18 mm·(10a)-1和31.42 mm·(10a)-1。生长季内太阳辐射量基本呈下降趋势,其中林芝站达到显著水平,下降幅度为72.91 MJ·m-2·(10a)-1。从变化趋势看,青藏高原地区青稞生长季内气象总体呈增温增湿趋势,增温更加明显。
表3 1977—2017 年青藏高原站点生长季内气候变化倾向率Table 3 Trend of climatic factors during the growing season of highland barley on Tibetan Plateau from 1977 to 2017
图1 青藏高原研究站点1977-2017 年青稞生长季内的气候分布Fig. 1 Climatic conditions of highland barley growing season from 1977 to 2017 at research stations on Tibetan Plateau
利用模型的自动(DSSAT 中的GLUE 方法)和手动调参相结合的方法,每个站点用农气资料的校准集对这7 个站点青稞参数进行校准,并与站点验证集的实际观测结果对比。模型模拟与站点观测结果的对比情况如图2 所示,模拟的开花期、成熟期和生育期长度与实际观测相比,R2为0.84、0.80 和0.66,NRMSE为4.58%、3.05%和3.76%。产量模拟R2达到0.60。由此可知,DSSAT-CERES-barley 可对青藏高原青稞进行模拟。
图2 DSSAT-CERES-barley 模拟生育期长度,开花期,成熟期,产量与实测验证集比较Fig. 2 Comparison of simulated growing period, flowering, maturity, and grain yield by DSSAT-CERES-barley and observed data
各站点在无水氮胁迫条件下1977—2017 年青稞的生育期长度变化情况如图3 所示。由地理分布可知,海拔越高地区青稞生育期长度越短,所有站点在表1播种期固定的情况下,过去40 年的生育期长度均显著缩短(P<0.01),下降幅度为0.257—0.843 d·a-1,其中,门源站生育期长度缩短幅度最大,可达0.8431 d·a-1。这说明在播期不变的条件下,过去40 年气候变化导致所有站点的青稞生育期长度显著缩短,但不同海拔之 间下降幅度无显著差异。
由图4 可知,青藏高原1977—2017 年不同站点之间光温生产潜力差异较大。山南站平均光温生产潜力最高,达到11 338.98 kg·hm-2,门源站光温生产潜力最低,只有5 662.44 kg·hm-2,山南站平均光温生产潜力为门源站的2 倍,主要是由于山南站生长季有效积温和太阳辐射量分别高出门源站点11.5%和64.6%。此外,甘南州站平均光温生产潜力为5 875.43 kg·hm-2、贵南站为6 027.44 kg·hm-2、林芝站为5 898.78 kg·hm-2、日喀则站为10 425.59 kg·hm-2,拉萨站为7 266.17 kg·hm-2,分别是山南站光温生产潜力的51.8%、53.2%、52.0%、91.9%和 64.1%。总体趋势表明青藏高原高海拔站点,青稞生长季内温度和辐射量越高,光温生产潜力越强。
图3 青藏高原1977—2017 年青稞生育期长度变化趋势Fig. 3 Trends of growth period of highland barley in Tibetan Plateau from 1977 to 2017
图4 青藏高原1977—2017 年青稞光温生产潜力变化趋势Fig. 4 Trends of photo-temperature potential productivity of highland in Tibetan Plateau from 1977 to 2017
由趋势分析可知,甘南州和门源站的光温生产潜力在1977—2017 年间显著升高,幅度为22.96 kg·hm-2·a-1和22.4 kg·hm-2·a-1。而林芝站光温生产潜力显著下降,幅度为10.25 kg·hm-2·a-1。海拔在3 500 m 左右的站点光温生产潜力较为稳定,低海拔站点由于气候变化导致光温生产潜力波动(图4)。
由图5 可知,青稞生育期长度与温度相关的气象因子均呈负相关,随青藏高原气候变暖,所有站点生育期长度显著缩短。但不同海拔站点之间受气候变化影响程度不同,温度较低的低海拔地区平均温度和有效积温变化与生育期长度变化呈负相关,且相关系数最大。高海拔地区站点平均最高温度变化与生育期长度呈负相关且相关系数最大,这主要是由于这些站点积温充足,但平均最高温度的上升仍会导致生育期长度缩短。由此可知,均温较低的甘南州、门源和贵南站,其生育期缩短主要是由于生长季内平均温度和有效积温升高引起,而在平均温度较高的雅鲁藏布江河谷地区站点主要由于平均最高温度升高导致生育期长度缩短。
图5 气候变化与光温生产潜力和生育期长度变化(P<0.05)的相关系数Fig. 5 Correlation coefficient (P<0.05) between climate change and photo-temperature potential productivity and growth period change
为避免过去几十年产量、物候期和气候因子本身的变化趋势对相关性产生干扰,选择光温生产潜力和气候因子的变化进行相关性和回归分析。由于逐步回归过程中平均气温与其他变量之前存在严重的多重共线性。因此,选择除生长季内平均温度之外的气象因子,对光温生产潜力和生育期长度进行逐步回归分析。由图5 可知,对于海拔较低的贵南、门源和甘南州站点,光温生产潜力与生育期长度呈负相关,与温度相关因子成正相关,这表明低海拔站点青稞生育期内温度未能充分满足青稞的生长发育,虽然气候变暖导致生育期长度缩短,但对光温生产潜力的提升作用更大。而高海拔站点光温生产潜力变化与生育期长度变化无相关性,且与气候变化的显著相关系数较小,说明高海拔站点光温生产潜力稳定,生育期内光温条件可满足青稞生长需要。
由表4 可知,生长季内太阳辐射对光温生产潜力均有显著敏感性,但不同站点的敏感程度不同,每升高1 MJ·m-2,光温生产潜力升高1.60—3.95 kg·hm-2,且海拔越高,太阳辐射强度越大的地区,太阳辐射对光温生产潜力的敏感性更高。
根据过去30 年青稞实际统计产量资料,分析了1987—2017 年青藏高原青稞产量差的时空分布。由图6 可知,不同站点之间绝对产量差差别较大,因此通过研究相对产量差的变化趋势可以更好反映不同站点之间的差别。除林芝站点外,过去30 年里光温生产潜力与实际产量之间相对产量差均有所减少,表明实际产量增加是导致产量差缩小的主要原因。青藏高原青稞在1987—1996 年光温生产潜力与实际产量之间的相对产量差在45.3%—65.6%,站点均值为58.2%;1997—2006 年的相对产量差在25.7%—56.4%,站点均值为41.9%;2007—2017 年相对产量差在15.6%—47.4%,站点均值为34.5%。海拔较高的日喀则站点相对产量差变化幅度最大,2007—2017 年相对产量差只有15.6%,减少41.7%。只有林芝站2007—2017 年相对产量差有所增加,因为光温生产潜力和统计产量相比1997—2006 年均有所减少,原因是日照时数和太阳辐射的显著减少所致,从而造成相对产量差降低4.0%。
表4 不同站点的青稞光温生产潜力变化与气候因子的逐步回归方程Table 4 Stepwise regression equation between photo-temperature potential productivity change of highland barley and climatic factors in difference stations
图6 青藏高原1987-2017 年青稞潜在产量与实际产量对比图Fig. 6 Comparison of photo-temperature potential productivity and statistical yield of highland barley in Tibetan Plateau from 1987 to 2017
青藏高原过去40 年所有站点青稞生长季内最高、最低和平均气温均显著上升,最低气温的升高幅度大于平均和最高气温。大部分站点生长季内降水量上升,生长季内太阳辐射量下降,这与青藏高原整体气候变化趋势相吻合[6-7]。因此在假设播期和品种不变的条件下,通过模型对过去40 年光温生产潜力和生育期长度进行模拟,从而分析得出气候变化对潜在产量及产量差的影响。
甘南州、门源、贵南和林芝站点光温生产潜力平均在6 000 kg·hm-2,拉萨站为7 000 kg·hm-2,日喀则和山南站为11 000 kg·hm-2。通过分析可知,光温生产潜力在青藏高原的分布与海拔变化趋势相同。青藏高原地形特殊,海拔从东北向西南逐渐增加,纬度越低海拔越高的站点光温生产潜力越高。这与华北平原光温生产潜力的地理分布不同[15],说明青藏高原不同站点之间海拔和纬度差异较大,同时品种不同造成海拔较高的站点光温生产潜力显著高于低海拔站点。
在无水氮胁迫且播种期不变的条件下,青藏高原气候变暖导致青稞生长季内的升温,造成过去40 年间所有站点青稞生育期长度显著缩短。海拔较低的甘南州、贵南和门源站点,主要是生长季内有效积温和平均温度升高导致。而高海拔站点主要是由于平均最高气温升高导致,且海拔越高敏感性越强。通过相关分析及回归分析可知,过去40 年青稞生长季气候变化对海拔较高的日喀则、山南和拉萨站点的光温生产潜力无显著影响。但是对低海拔站点影响显著,特别是太阳辐射量的下降会导致作物净光合速率下降,干物质积累减少[43],从而使光温生产潜力下降。但只有林芝站累计太阳辐射量下降幅度最大,达到72.91 MJ·m-2·(10a)-1(表3),且林芝站光温生产潜力只与太阳辐射量呈正相关,这也造成过去40 年林芝站光温生产潜力显著下降。而其他站点由于太阳辐射量下降幅度较小,太阳辐射量对光温生产潜力影响较小。
由于甘南州、贵南和门源站在其生长季内太阳辐射和温度较低,光温生产潜力低而不稳定。特别是青藏高原相比其他地区属于高寒区[44],温度升高不仅抵消生育期长度缩短和太阳辐射下降的影响,反而在一定程度上有利于青稞提高光合作用,增加地上部分生物量并延长灌浆时间,提高籽粒灌浆率,提高产量潜力[45]。因此海拔较低的青藏高原东北部地区过去40 年光温生产潜力总体呈显著上升。这也说明DSSAT-CERES-barley 模型考虑了作物、土壤以及气候的动态变化对青稞生长发育的影响,对光温生产潜力的模拟更加准确。通过西藏自治区高海拔站点光温生产潜力的计算可知西藏自治区光温潜力受气候变化影响较小,且海拔较低的林芝站受辐射变化光温生产潜力下降,只有青藏高原东北部低海拔地区光温生产潜力上升。这与过去计算的西藏自治区光温生产潜力(平均为14 000 kg·hm-2)且受气候变化影响呈上升趋势的结果不同[20]。说明虽然经验公式法计算光温生产潜力的运算过程简单,计算量小,但误差较大,缺乏对影响青稞生长发育因素的考虑,且计算结果不仅高估了青藏高原青稞的光温生产潜力,而且趋势分析结果也有较大误差。本研究是假设相同站点品种和播期不变条件下进行模拟的结果,原因是青藏高原地区农业气象观测数据较少,最早有记录的品种和播期数据为1987 年,相同站点青稞播种期接近。虽然该假设会导致光温生产潜力模拟结果的不确定性,使模拟的光温生产潜力没有准确体现现实结果,却可以有效分离气候变化对光温生产潜力的影响,从而反映青藏高原青稞光温生产潜力的时空分布特征并推测未来变化趋势。
结合青藏高原未来气候变化的研究可知[46-47],未来几十年青藏高原的气候尤其是青海和西藏的青稞种植区域会继续保持变暖和增湿趋势,因此位于东北部海拔较低的甘南州、贵南和门源地区青稞光温生产潜力将继续保持上升趋势;而林芝站由于降水量升高,太阳辐射量降低,未来光温生产潜力可能继续下降;通过分析雅鲁藏布江河谷海拔较高的山南、日喀则和拉萨站,发现其光温生产潜力不易受气候变化的影响,因此未来将继续维持较高水平。
通过产量差分析可知,虽然门源、林芝和甘南州的光温生产潜力具有显著性变化,但幅度较小。总体来说,青藏高原光温生产潜力较实际产量相对稳定,而实际产量在过去30 年里增长趋势显著,导致过去30 年间青藏高原青稞的平均相对产量差由58.2%缩小到目前的34.47%。实际产量升高的主要原因有:(1)育种水平提高,且良种储备充足;(2)科技投入增加与技术的进步;(3)西部大开发战略以及政策扶持也起到关键作用[48-49]。地理位置上,不同站点之间的产量差差异较大。日喀则和拉萨站在2007—2017 年间的产量差达到25%以下,其他站点产量差仍大于30%。从时间上看,前30 年到前20 年产量差缩小16.3%;过去20 年到最近10 年,产量差缩小7.5%,只有过去的50%,表明了进一步缩小产量差的困难。但除日喀则和拉萨之外,其他地区至少有5%—10%的产量差空间可通过水肥管理等措施改善。进一步缩小青藏高原青稞的产量差,可以通过进一步提升灌溉和管理条件、提高水肥利用效率实现。
在过去40 年青稞生长季内气候变化的大背景下,海拔较高的站点光温生产潜力高且稳定,生育期长度的缩短主要由平均最高气温升高引起;低海拔站点光温生产潜力较低,易受气候变化和生育期长度变化影响,而生育期长度缩短主要由平均温度和有效积温升高导致。因此低海拔地区可通过改良品种,改变播期等方式增加生育期长度,这相比于高海拔地区可更有效地增加青稞的光温生产潜力。对比光温产量与实际产量,实际产量增加对产量差变化的贡献更大。相比于其他站点,日喀则和拉萨产量差较低,说明这2 个地区青稞种植的管理水平和水肥利用效率较高,而其他站点则有较高的增产潜力。