基于时序MODIS和气象数据的日喀则地区青稞估产

潘莹，仙巍，陈春容，陈军，徐薇

(1.成都信息工程大学 资源环境学院，成都 610225；2.日喀则市气象局，西藏 日喀则 857000)

0 引言

青稞是藏区人民的主食，也是海拔4 500 m以上唯一可以正常成熟的粮食作物。青藏高原是全球青稞产量和种植面积最大的地区，而日喀则地区作为西藏的第二大城市，也是青稞的主产区。因此，快速准确地监测青稞的产量，对藏区的产业发展与经济发展具重要意义。

当前，对青稞估产的相关研究较少，而估产可分为作物的种植面积提取与估产模型建立两部分。日喀则地区地形复杂，青稞种植区域分布零碎，传统的调查估产手段难以得到全面准确的结果，但遥感技术可以解决其无法全覆盖和成本过高的问题，具有较好的准确度和时空连续性的优势，在种植面积提取和产量预测中起着重要的作用[1]。利用遥感手段实现产量预测方法主要有统计预测法和生长模型预测法两大类，其中，生长模型预测法是数值模拟作物的生长过程，需引入大量的影响作物生长的参数与同化后的遥感数据，故在大范围的作物估产中难以实现[2-3]；而统计预测法只需选入少量的影响产量的相关参数进行估产，如气象数据和遥感指数，对大面积的区域估产相对较易实现。遥感数据估产原理是利用作物的光谱反射特征提取其潜在的产量信息，并且光谱植被指数可以定量化，其中最为常见的归一化植被指数(normalized difference vegetation index，NDVI)光谱组成与作物产量的物理关系在机理上已得到证实[4-6]。刘宝锋等[7]运用遥感分析法，结合Landsat TM数据和历年青稞产量数据，建立了甘孜县青稞估产模型并取得了较好的结果。刘新杰等[8]基于时序卫星的定量NDVI遥感数据实现了冬小麦的长势监测和精确估产。通过遥感技术进行农作物估产主要是建立经验模型，而回归模型是典型的代表，王长耀等[9]、任建强等[10]利用MODIS植被指数产品，采用逐步回归法建立农作物估产模型，取得了较高的估产精度。

气候变暖正影响着脆弱的高原生态环境，青藏高原作为中国乃至世界气候变化的敏感区和启动区，农业生产严重依赖于气候资源[11-12]。日喀则地区的青稞产量对气候变化存在明显的响应，特别是对气温和降水量异常敏感[13]，而目前基于单纯的遥感数据的产量预测模型较多，因此本文结合时间序列的MODIS遥感数据和气象数据，拟建立一种适用于青藏高原日喀则地区的综合的青稞遥感气象估产模型。

1 数据与方法

1.1 研究区概况

本文以我国西藏日喀则地区为研究区，其地处青藏高原东南部，位于82°00′E～90°07′E、27°55′N～31°16′N，面积为18.2×104km2，地形复杂多样，平均海拔在4 000 m以上。日喀则地区空气稀薄，太阳辐射强，日照长，气温低，年较小，日较大，常年农作物播种面积为2.6×105ha。青稞是日喀则地区播种面积最大的农作物，主要集中在雅鲁藏布江和年楚河流域，由拉孜—仁布宽谷和江孜—日喀则平原组成，还包括喜马拉雅山脉北侧、藏南高原上的朋曲河谷平原和其他零星的河谷平原。这些谷地坡度平缓，土层深厚，气候宜人，水源较充足，有利于农作物的种植。

1.2 数据准备

本文所需数据有遥感数据和非遥感数据。其中遥感来源于美国国家航天局(National Aeronautics and Space Administration，NASA)的官方网站(https：//ladsweb.modaps.eosdis.nasa.gov)，获取2015—2019年16 d合成的MOD13A1和8 d合成的MOD09A1的NDVI，分别用于青稞的种植区提取与估产。获取行列号为h25v05与h25v06，其中MOD13A1共230景，MOD09A1共460景。使用MRT输入命令行可对以上数据进行重采样和镶嵌的批处理。非遥感数据包括气象数据和统计数据，分别来源于中国气象数据网(http：//data.cma.cn/)和《日喀则统计年鉴》，包括日喀则地区2016—2019年各市(县)的日平均气温、日降水量、青稞单产产量和该地区的播种总面积。

2 研究方法

2.1 青稞种植区提取

随着农作物的物候期更替，近红外波段因作物的叶面积系数和覆盖率增大而增大，红光反射率则逐渐减小，因此在某一时间断面上，日喀则地区的青稞和冬小麦的状态和结构会因物候期的不同而有明显差别[14-15]。日喀则地区的冬小麦比青稞要早一个月进入生长最茂盛的时期即拔节期，利用这种不同作物间物候差异特征，通过MOD13A1的NDVI数据，经过反复的探究与实验，得到式(1)所示的日喀则地区青稞种植区提取方法。该方法不仅可以提取青稞的种植区，还能有效剔除日喀则地区以冬小麦为主的其他农作物与绿色植物的影响。

(1)

式中：QK是青稞种植范围，当QK为1时，为种植区域，当QK为0时，为非种植区；NDVI拔节为拔节期的NDVI；NDVI播种前为播种前的NDVI。

为提高青稞种植范围的提取精度，将利用式(1)提取的结果与高清的天地图进行人工目视解译对比，剔除个别未在耕地上的错误值，最终得到日喀则市青稞2015—2019年的种植面积。

将青稞提取总播种面积与统计总播种面积相比，2019年的提取精度最高，为99.51%，2017年提取精度最低，为97.46%(表1、表2)。历年来，位于年楚河流域的日喀则市、江孜县和白朗县青稞播种面积较多，约占日喀则地区青稞播种总面积的1/2；亚东县的播种面积是所提取的县中最少的，其中2018年和2019年由于青稞播种面积过小，难以提取。分析统计数据发现，2015年总播种面积最小，为1 293.1 km2，在2018年达最大，为1 395.8 km2，2019年总播种面积为1 391.90 km2。

表1 日喀则地区青稞种植面积提取结果 km2

表2 日喀则地区青稞种植面积提取精度

2.2 青稞产量预测方法

1)时间序列分析法。时间序列分析法是一种基于数理统计与随机过程理论的方法，研究的是连续时间的动态数据所遵循的统计规律[16-19]。青稞生长过程曲线的变化可以反映其长势的变化，而长势会影响最终的产量[20]。利用时间序列植被指数曲线分析方法，提取遥感波段反射率，根据青稞生长的时相变化特点可确定遥感监测的周期，从而选择恰当的遥感数据进行估产(图1)。植被指数时间序列滤波法是去除噪声的一种常用方法，利用S-G(savizky-golay)滤波法可去除云和大气对NDVI的影响，可得到去噪后的时间序列NDVI数据。

图1 青稞时间序列植被指数曲线图

2)产量预测模型。农作物的生产是一个复杂的综合系统，不仅会受到农业生物、技术、经济的影响，也会受到各气象条件等环境因素的影响。分解模型法是一种同时考虑农作物产量的构成因素(农业技术、农业环境和随机因素等)以及这些因素对产量的影响的方法。模型分解法可将农作物产量分解为趋势产量、波动产量与随机产量，如式(2)所示[21]。趋势产量为农业技术进步导致的产量部分，波动产量为农作物的自然环境导致的产量部分，随机产量是由突发自然灾害和其他不确定因素导致的产量部分(随机产量只能在模型外单独处理，本研究不予考虑)。趋势产量可通过区域产量序列拟合得出，表现为时间的正函数。波动产量为实际产量与趋势产量之差，可通过波动产量与气象参数进行相关性分析，以及反映作物生长情况的植被指数进行时间序列分析，并与产量进行相关系数分析，找出影响波动产量状况的遥感气象参数，建立波动产量预测模型，再结合趋势产量，建立单产预测模型。

YE=XT+XF+E

(2)

式中：YE为预测单产；XT为趋势产量；XF为基于遥感气象因子的波动产量；E为随机产量。

3 青稞产量预测模型的建立

3.1 植被指数时序特征分析

利用8 d的NDVI数据(MDO09A1)监测青稞的生长情况，提取日喀则地区各县(市)的青稞种植区所对应的NDVI，获得日喀则地区2015—2019年原始NDVI时间序列数据集，利用S-G滤波法去除云、气溶胶等对时间序列NDVI的影响，得到S-G去噪后的NDVI时间序列数据集并绘制对应的时间序列变化图。以江孜县为例(图2)。从图2可知，S-G去噪后的NDVI图像的噪声明显减少，其时间序列的NDVI值变化更稳定。分析去噪后的NDVI时间序列曲线可知，每年都具有一个明显的大波峰，峰值出现时间与日喀则地区青稞的抽穗期相吻合。在整个时间序列上，日喀则地区的青稞抽穗期NDVI值相对稳定，基本为0.6左右，说明青稞种植生长情况基本较稳定。

图2 日喀则地区江孜县NDVI时间序列变化图

3.2 植被指数的相关系数分析

NDVI反映的是作物生长的总体情况，与作物的产量有直接的关系，因此在进行相关性分析时，应该将NDVI与产量直接进行相关分析，然后将NDVI作为影响因素，引入波动产量的预测模型中，作为NDVI对产量的影响。图3是研究区时序NDVI与单产的相关系数随时间的变化曲线，当青稞种植区NDVI值越高，说明青稞长势越好，所对应的单产也越高。通过时序的NDVI与单产的相关系数曲线，不仅可以看出何时的NDVI与单产的相关系数最高，也能了解日喀则地区青稞整体的生长状况。

图3 NDVI与单产的相关系数时序变化曲线

由图3可知，NDVI与单产的相关系数随时间变化总体呈现抛物线形状，呈先上升后下降的趋势。在时间序列曲线的起始部分，NDVI与单产的相关系数较低，这是因为在3月份青稞还没有开始播种。而相关系数曲线大约在4月底五月初开始出现较大幅度的上升，即在137至161天之间，青稞播种之后开始出苗、分蘖、拔节、抽穗、孕穗，是青稞生长的关键时期，其NDVI值增大，相关系数也增大，对于单产多的地区，其青稞长势自然也更好。而在169至233天之间(6月中旬至8月中旬)，日喀则市整体的NDVI值达到饱和并处于相较稳定的阶段，故其与单产的相关系数也比较稳定。观察该曲线，其最大值出现在209天(7月下旬)，即青稞的抽穗期，而日喀则地区的青稞也在7月下旬时生长达到最旺盛的阶段，所以青稞抽穗期的NDVI是日喀则市青稞估产的最优选择。

3.3 气象因子的确定

日喀则地区的青稞产量对气温和降水具有较好的响应特征。统计日喀则地区各市(县)青稞生长阶段的月积温、月平均气温和累计降水量，并提取青稞生长重要物候期内的积温，将以上气象因子与波动产量进行相关性分析，计算结果如表3所示。选取与波动产量相关性最大的一个气象指标可作为构建农业气象估产模型的主要因子。由表3可知，日喀则地区的7月上旬积温与青稞波动产量的相关性最高，因此7月上旬积温可作为构建农业气象估产模型的主要因子。进一步分析表3发现：波动产量与积温和平均气温的相关系数均为负值，这是因为青稞属于喜凉作物，温度的升高会对青稞产量起到负效应[22]；青稞生长的关键期，累积降水与波动产量的相关系数也为负数，可能是因为降水过多时云量大、太阳辐射也减少，而当地过高的土壤含水量导致土壤蒸散增加会使近地表温度降低，进而降低了青稞的光合作用能力降低，从而对青稞产量起到负效应[23]。

表3 日喀则地区各气象因子与波动产量的相关关系

3.4 产量预测模型

选取最为敏感的遥感气象参数作为构建青稞遥感气象估产模型(式(2))的主要因子，并检验参数的稳定性，从而建立一元线性或多元线性回归方程。波动产量的回归方程可以分为三种(表4)：气象波动产量回归方程、遥感波动产量回归方程和遥感气象波动产量回归方程。其中，XF-m为只考虑气象因子影响的波动产量；XF-rs为只考虑NDVI的遥感因子的遥感波动产量；XF为既考虑气象因子又考虑NDVI遥感因子的遥感气象波动产量；Xm为气象因子中的7月下旬的积温；Xrs为遥感因子的7月下旬的NDVI(MDO09A1的209天的NDVI)。

表4 三种波动产量回归方程

基于以上三种计算波动产量的回归方程，然后根据2019年的遥感气象数据计算2019年日喀则地区各县中青稞的遥感气象波动产量，并基于式(2)结合趋势产量得到青稞单产的预测结果(表5)。

表5 2019年日喀则地区各县青稞单产预测结果 kg·hm-2

3.5 误差分析

将预测的青稞单产与统计数据进行比较和相关性分析，验证模型的显著性和拟合精度。由表6可知：建立的青稞产量预测模型中，遥感气象估产模型的预测结果最为稳定，误差最大的绝对值为16.19%，相对误差控制在-14.34%～16.19%；气象估产模型误差最大的绝对值为19.66%，相对误差范围控制在-9.59%～19.66%；遥感估产模型误差最大的绝对值为19.14%，相对误差控制在15.41%～19.14%；平均估产误差(即各市(县)所有模型的相对误差的绝对值之和的平均)为7.53%。由图4可知，气象、遥感、遥感气象预测单产与统计单产的相关系数分别为0.943 6、0.946 1、0.938 2，这说明三种预测模型得到的预测单产与统计单产都具有良好的线性关系。

表6 2019年日喀则地区各县青稞单产预测误差 %

图4 2019年日喀则地区青稞预测单产与统计单产对比

4 结束语

本文提出了一种利用MOD13A1快速准确地提取青稞的方法，2015—2019年青稞的提取精度均大于97.4%。日喀则地区青稞主要分布在雅鲁藏布江和年楚河流域的河谷平原，西北降水少，海拔高，农田很少，提取结果与高子恒等[24]得到的日喀则地区的农田结果相似，具有一定的可靠性。

日喀则地区的遥感估产的敏感因子为7月下旬(抽穗期)的NDVI(MDO09A1的209天的NDVI)。基于时间序列分析，对日喀则地区的青稞进行长势监测，近五年青稞生长曲线波峰的NDVI都在抽穗期，稳定在0.6左右，说明青稞种植生长情况稳定；青稞NDVI与单产的相关系数在生长期呈抛物线，7月下旬(抽穗期)达到最大值0.709，青稞抽穗期的NDVI是日喀则地区青稞估产的最优选择。

日喀则地区气象估产的敏感因子为7月上旬的积温。日喀则地区青稞产量对气候变化存在明显响应，且对气温的敏感性大于降水量。

本文构建的日喀则地区青稞单产预测模型的平均估产误差为7.53%，平均相关系数为0.942 6，在一定程度上能满足当地的青稞估产需求。利用模型分解法可建立青稞的气象单产预测模型、遥感单产预测模型和遥感气象单产预测模型。其中，遥感气象单产预测模型最稳定，最大的相对误差控制16.19%以内。所有预测结果与统计单产都具有良好的线性关系，相关系数均大于0.938 0。因此，本文所构建的估产模型在日喀则地区青稞生产区具有一定的可操作性和适用性。

由于日喀则地区的范围较广，气象数据受到区域范围的影响，会对利用气象因子预测的单产造成一定的误差；而利用遥感数据提取各年份的青稞总播种面积精度均大于97.4%，但由于混合像元的原因，也会对利用遥感因子预测的单产造成误差。另外，个别县(市)的青稞单产数据年际差别较大，也是造成预测结果误差偏大的原因，如定结县、岗巴县、拉孜县和日喀则市。日喀则地区的遥感和气象估产的敏感因子分别为7月下旬即青稞抽穗期的NDVI(MDO09A1的209天的NDVI)和7月上旬的积温。气象因子中积温对波动产量的影响最大，不同的气温累积会直接影响物质的累积与蛋白质、淀粉的含量，从而影响产量[25]。日喀则地区处于高海拔的青藏高原西南部，在青藏高原强烈的选择压力下，海拔和系统发育综合因子会对植物千粒重产生影响，进而会对青稞产量造成一定的影响。日喀则地区的青稞种植区的海拔落差大约有1 700 m，由于较大的海拔差异，对温度影响显著，主要表现随着海拔的升高，温度降低，作物生长季短，因此会造成一定的物候时间差异，影响产量预测的结果[26]。

本文弥补了当前针对高原地区大范围青稞估产的空白，构建了适应于当前发展趋势的综合遥感气象农作物估产模型，对高原地区的粮食生产与经济发展有重要的指导意义。在模型的实际应用中，需要对不同海拔高度的县进行实地调查，利用抽样统计可对作物单产模型进行修正，从而提高模型精度。