本发明属于水文气象及农业交叉学科数据分析领域，尤其涉及一种作物产量对气候变化响应的诊断方法。
背景技术：
气象要素(水，热，温，湿等)在作物生长生产过程中起着重要作用，决定着作物产量的高低。近年来，水文、气象、农业工作者基于变化环境，对变化环境下水-能源-作物产量相互作用机理开展了大量研究(ahmedetal.,2015；lietal.,2014；lobelletal.,2011；longetal.,2006；monteith&moss,1977；ottmanetal.,2012；porter&semenov,2005；schlenker&roberts,2009；zhaoetal.,2017)。气候变化对作物产量的影响被大量观测资料所证实。研究表明全球变化引起的气温变化和极端气候事件使得全球约三分之一的粮食产量受到影响(rayetal.,2015；vogeletal.,2019)。然而，由于缺乏长期的大范围的观测数据资料，加上外界条件如环境、农业种植施肥灌溉等技术的发展，使得现有的分析方法很难将这些外界条件对作物产量的影响从气候变化对作物产量的影响中剥离出来。
目前，统计分析方法，比如方差分析方法(nicholls,1997；porter&semenov,2005)、线性相关和线性回归方法(challinoretal.,2014；lobell&field,2007；lobelletal.,2011)在作物产量-气候变化相关研究中广泛应用。然而，现有研究结果在空间分布上存在着较大的差异性和地域性。基于非参数弹性分析方法在作物产量-气候变化响应中的应用目前暂无相关研究。
技术实现要素：
发明目的：为解决现有技术的不足，本发明首次将非参数弹性分析方法应用在作物产量对气候变化响应的诊断分析中，提出一种新型作物产量对气候变化响应的诊断方法。
技术方案：为实现本发明的目的，本发明所采用的技术方案是：一种新型作物产量对气候变化响应的诊断方法，包括如下步骤：
步骤1，选择研究区域，收集区域气象要素时间序列和农作物产量时间序列；
步骤2，根据搜集到的气象要素和农作物产量时间序列选取计算时段；
步骤3，对选定时段内的气象要素和农作物产量数据序列去趋势分析；
步骤4，计算去趋势分析后气象要素时间序列均值、偏差和偏差率；
步骤5，计算去趋势分析后农作物产量时间序列均值、偏差和偏差率；
步骤6，根据去趋势气象要素和农作物产量偏差率时间序列，计算敏感指数时间序列，求中位数，获得选定时段内敏感指数εa；
步骤7，根据敏感指数εa判断农作物产量对气候变化的响应。
进一步的，所述步骤1，所选取的研究区域不限，可以是地区、国家或全球尺度的农作物种植区域。所选取的气象要素不限，可以是气温、降水、蒸散发。所选取的农作物种类不限，可以是小麦、水稻、玉米。
进一步的，所述步骤2，所选取的研究时段不限，可以是十年、二十年或更长时段，但选取的气象要素和农作物产量时间序列长度需一致。
进一步的，所述步骤3，采用tfpw(trendfreepre-whitening)去趋势方法，对选定时段内的气象要素和农作物产量数据序列去趋势分析，公式如下：
x′i＝xi-βi(2)
y′i＝x′i-r1x′i-1(3)
yi＝y′i βi(4)
其中，xi和xj表示选定的初始要素x在时刻i和时刻j的量，x是选定的气象要素或作物产量数据序列，β是初始要素的增减趋势，median为取中位数；x′i是初始要素x在第i时刻去趋势后的量；r1是间隔一个时段的去趋势后的要素序列x′i的自相关系数；y′i为去趋势分析后的x′i移除自相关系数后的要素序列；yi为去趋势分析后的要素序列。
进一步的，所述步骤6，获得选定时段内敏感指数的方法如下：
其中，εa为敏感指数，ri为去趋势分析后第i时段的农作物产量，ai为去趋势分析后第i时段的气象要素，
进一步的，所述步骤7，根据敏感指数εa判断农作物产量对气候变化的响应，方法如下：
当εa>0时，表明农作物产量对气候变化的响应为正反馈；当εa<0时，表明农作物产量对气候变化的响应为负反馈；当εa＝0时，表明农作物产量对气候变化不响应。
敏感指数εa的量值大小表明农作物产量对气候变化的响应强度，εa量值愈大，响应强度愈大；当|εa|≥1时，表明响应显著；当|εa|<1时，表明响应不显著。
有益效果：与现有技术相比，本发明的技术方案具有以下有益的技术效果：
本发明通过特定时段气象要素和农作物产量时间序列，计算敏感指数，确定农作物产量对气象要素变化的响应特征，由此获得农作物产量对气候变化的响应，可为农业生产过程中适应气候变化，结合自然条件，为最优化生产灌溉技术提高作物产量提供理论依据，具有重要的社会生产意义。
附图说明
图1是本发明的方法流程图；
图2是研究区域小麦丰收区域空间分布图(单位：公顷)；
图3是研究区域小麦种植期空间分布图；
图4是研究区域小麦收获期空间分布图；
图5是研究区域小麦生长期月平均气温空间分布图(单位：℃)；
图6是研究区域小麦对气温的敏感指数空间分布图。
具体实施方式
下面结合附图和实施例对本发明的技术方案作进一步的说明。以下实施例仅用于更加清楚地说明本发明的技术方案，而不能以此来限制本发明的保护范围。
现有0.5°*0.5°空间分辨率的1961-2014年间的亚洲年平均小麦产量数据(数据来源于faostatisticaldatabase(http://faostat.fao.org))和cru(climateresearchunit)ts3.1数据(http://crudata.uea.ac.uk/cru/data/)，研究区域经纬度为[30e150e,-10n50n]。依据本发明方法，流程如图1所示，以气象要素气温为例，计算研究区域内1961-2014年间小麦产量对气温变化的响应特征。过程如下：
(1)根据小麦的丰收区域(图2)，种植期(图3)和收获期(图4)，在研究区域内单个格点上，以种植期和收获期之间的月平均气温，计算小麦生长期的平均气温。若在某个格点上的小麦的种植期为4月，收获期为同年11月，则该小麦生长期的平均气温为该年4-11月的月平均气温。若在某个格点上小麦的种植期为10月，收获期为次年5月，则该格点上小麦的生长期的平均气温为该年10-12月和次年1-5月的月平均气温。因为，小麦种植期和收获期在不同格点上具有差异性，因而生长期在空间分布上具有差异性。根据该方法计算所得的研究区域内1961-2014年间小麦生长期的月平均气温的空间分布图，如图5所示。
(2)根据研究区域内各格点上小麦生长期的平均气温和小麦年产量时间序列，采用yueetal.(2002)提出的tfpw方法去趋势分析，获取1961-2014年间去趋势分析后的小麦生长期的平均气温和小麦年产量时间序列。
(3)计算研究区域内1961-2014年间各格点的去趋势分析后的小麦生长期的平均气温和小麦年产量的均值；
(4)计算研究区域内1961-2014年间各格点的去趋势分析后小麦生长期的年平均气温和小麦产量的偏差和偏差率时间序列；
(5)计算研究区域内1961-2014年间各格点的去趋势分析后小麦产量对生长期气温的敏感指数时间序列；
(6)计算研究区域内1961-2014年间各格点的去趋势分析后小麦产量对生长期气温的敏感指数时间序列的中位数，即为该格点的敏感指数，绘制研究区域内小麦产量对气温的敏感指数空间分布图，如图6所示；
(7)根据研究区域内各格点的敏感指数空间分布图，如图6所示，判断小麦产量对气温变化的响应。当εa>0时，表明小麦产量对气温变化的响应为正反馈；当εa<0时，表明小麦产量对气温变化的响应为负反馈；当εa＝0时，表明小麦产量对气温变化不响应。εa的量值愈大，表明小麦产量对气温变化的响应强度越大。当|εa|≥1时，表明响应显著；当|εa|<1时，表明响应不显著。
本发明通过特定时段气象要素和农作物产量时间序列，计算敏感指数，确定农作物产量对气象要素变化的响应特征，由此获得农作物产量对气候变化的响应，可为农业生产过程中适应气候变化，结合自然条件，为最优化生产灌溉技术提高作物产量提供理论依据，具有重要的社会生产意义。
以上所述是本发明的优选实施方式，应当指出，对于本技术领域的普通技术人员来说，在不脱离本发明技术原理的前提下，还可以做出若干改进和变形，这些改进和变形也应视为本发明的保护范围。
1.一种新型作物产量对气候变化响应的诊断方法，其特征在于：该方法包括如下步骤：
步骤1，选择研究区域，收集区域气象要素时间序列和农作物产量时间序列；
步骤2，根据搜集到的气象要素和农作物产量时间序列选取计算时段；
步骤3，对选定时段内的气象要素和农作物产量数据序列去趋势分析；
步骤4，计算去趋势分析后气象要素时间序列均值、偏差和偏差率；
步骤5，计算去趋势分析后农作物产量时间序列均值、偏差和偏差率；
步骤6，根据去趋势气象要素和农作物产量偏差率时间序列，计算敏感指数时间序列，求中位数，获得选定时段内敏感指数εa；
步骤7，根据敏感指数εa判断农作物产量对气候变化的响应。
2.根据权利要求1所述的一种新型作物产量对气候变化响应的诊断方法，其特征在于：所述步骤1，所选取的研究区域不限，包括地区、国家或全球尺度的农作物种植区域；所选取的气象要素不限，包括气温、降水、蒸散发；所选取的农作物种类不限，包括小麦、水稻、玉米。
3.根据权利要求2所述的一种新型作物产量对气候变化响应的诊断方法，其特征在于：所述步骤2，所选取的研究时段不限，选取的气象要素和农作物产量时间序列长度一致。
4.根据权利要求3所述的一种新型作物产量对气候变化响应的诊断方法，其特征在于：所述步骤3，采用tfpw去趋势方法，对选定时段内的气象要素和农作物产量数据序列去趋势分析，公式如下：
x′i＝xi-βi(2)
yi′＝x′i-r1x′i-1(3)
yi＝yi′ βi(4)
其中，xi和xj表示选定的初始要素x在时刻i和时刻j的量，x是选定的气象要素或作物产量数据序列，β是初始要素的增减趋势，median为取中位数；x′i是初始要素x在第i时刻去趋势后的量；r1是间隔一个时段的去趋势后的要素序列x′i的自相关系数；yi′为去趋势分析后的x′i移除自相关系数后的要素序列；yi为去趋势分析后的要素序列。
5.根据权利要求3所述的一种新型作物产量对气候变化响应的诊断方法，其特征在于：所述步骤6，获得选定时段内敏感指数的方法如下：
其中，εa为敏感指数，ri为去趋势分析后第i时段的农作物产量，ai为去趋势分析后第i时段的气象要素，
6.根据权利要求3-5任一所述的一种新型作物产量对气候变化响应的诊断方法，其特征在于：所述步骤7，根据敏感指数εa判断农作物产量对气候变化的响应，方法如下：
当εa＞0时，表明农作物产量对气候变化的响应为正反馈；当εa＜0时，表明农作物产量对气候变化的响应为负反馈；当εa＝0时，表明农作物产量对气候变化不响应；
敏感指数εa的量值大小表明农作物产量对气候变化的响应强度，εa量值愈大，响应强度愈大；当|εa|≥1时，表明响应显著；当|εa|＜1时，表明响应不显著。
技术总结