MODIS时序数据在县级尺度作物长势监测分析中的应用

于 堃， 王志明， 孙 玲， 单 捷， 毛良君

(江苏省农业科学院农业经济与信息研究所，江苏 南京 210014)

农业生产是人类社会存在的基础，随着全球气候的不断变化，农业生产面临着越来越多的灾害威胁［1］。为了提高农业生产对各种灾害的应对能力，降低灾害带来的经济损失，需要人们及时掌握农作物的生长状况。作物长势能够及时反映农情信息，通过对农作物长势的实时监测，不仅能够为农业生产的田间管理提供及时、科学的依据，也可为农作物产量估测提供必要的基础数据［2-3］。传统作物长势监测多依赖地面调查，造成农情信息获取成本高、时效性差，且调查结果常受主观人为因素影响。近年来，随着遥感技术的不断发展，尤其是遥感信息的覆盖范围广、数据获取周期短、费用成本低等特点，为作物长势信息的快速、准确、动态获取提供了重要的技术支撑，使传统作物长势监测中所遇到的难题在一定程度上得以解决［4-5］。

早在上世纪70年代，美国便开始应用遥感技术进行作物长势的监测，并形成了应用遥感数据监测全球作物长势的业务化产品［6-8］。之后，加拿大于上世纪80年代构建了基于NOAA/AVHRR以及MODIS等遥感数据的全国作物长势监测系统，并每7～10天向用户提供大尺度的作物长势数据产品［9］。巴西、法国、俄罗斯、日本、阿根廷和印度等国也相继构建了利用遥感数据获取作物长势信息的业务化运行平台，定期向用户发布作物长势数据产品［10］。中国也是较早利用遥感技术开展作物长势监测的国家之一［11］。上世纪80年代起，中国广泛开展了利用气象卫星监测作物长势及估产等方面的研究［12-14］。1998年中国科学院初步建立了国家级的农情监测系统，利用NOAA/AVHRR遥感数据获取全国尺度上1 km分辨率的旬农作物遥感长势分布图，并将作物长势监测作为一个主要的监测指标［15-16］。然而，对快速、准确的作物长势信息的需求不仅存在于全国尺度，也同样存在于中小尺度，如省级、市级以及县级农业相关管理及决策部门。因为大尺度的作物长势监测信息往往空间分辨率较低(1 km左右)，无法满足中小尺度的实际需求。但海量多源卫星遥感数据快速处理能力限制了其在中小尺度作物长势信息获取及发布中的应用，这造成了中国目前对于中小尺度尤其是县级尺度高时间分辨率作物长势的业务化遥感监测体系仍不完备［5-17］。以江苏省64个县市为例，若要利用MODIS数据获取近12年来覆盖每个县市的作物长势信息，需要处理的数据量多达560 640景。

本研究拟选取江苏省泗洪县为研究目标，利用近12年250 m地面分辨率的MODIS数据评价以旬为时间间隔的该县作物长势情况，结合野外实地观测数据对监测精度进行评价，探索利用中等分辨率的MODIS数据解决县级尺度作物长势信息的方法。由于MODIS遥感数据在实效性以及经济性等方面具有其他遥感数据无法比拟的优势，因此，将该方法推广应用到江苏省其他县市，不仅可以指导各县的农业生产、防灾减灾工作，也将推进江苏省县级尺度农情信息的业务化进程。

1 研究区概况与研究方法

1.1 研究区概况

江苏省是中国产粮大省之一，辖区内有64个县市。泗洪县位于苏北平原西部，东临洪泽湖，西接安徽省泗县，南靠江苏省盱眙县，北接江苏省宿豫县(图1)，属东亚季风区，四季分明，气候温和，光照充足，年均气温14.3℃，年均降水量893.9 mm，无霜期213 d，年均风速3.7 m/s。泗洪县辖区总面积为2 731 km2，其中耕地面积1 331.79 km2，总人口数约为7.808×105，该县的粮食总产量多年一直位列江苏省各县市的前茅，被评为中国商品粮基地县以及粮食生产先进县。泗洪县的粮食作物主要以小麦、水稻为主，经济作物主要有大豆、花生、山芋、西瓜等［18-19］。由于泗洪县特殊的地理位置，使其易受洪涝以及干旱等气象灾害的影响，统计数据表明，该地区平均每2.4年就会发生一次旱灾。同时，泗洪县的洪涝灾害也频繁爆发，如2003年、2007年泗洪县各地均不同程度的受到洪涝灾害的影响［20-21］。

图1 泗洪县地理位置图Fig.1 The location of Sihong county

1.2 数据来源

本研究所利用的遥感数据取自极地轨道环境遥感卫星Terra/Aqua上的MODIS(Moderate Resolution Imaging Spectroradiometer)所采集的遥感影像，数据由美国国家航空航天局(NASA)所属的戈达德航天中心(GSFC)提供。MODIS通过x波段将实时观测数据向全球免费发送。MODIS数据波段范围广，包括36个波段，分布在0.4～14.0 μm电磁波谱范围内，且具有很高的信噪比。其数据空间分辨率包括了250 m、500 m和1 000 m三个尺度，幅宽为2 330 km，垂直观测视场±55°。在对地观测过程中，可同时获得来自大气、海洋和陆地表面的信息，每1～2 d即可获取一次全球观测数据，因此被广泛应用于对陆表、生物圈、大气和海洋的长期全球观测［22-23］。

1.3 MODIS数据预处理

目前NASA针对MODIS开发了若干产品，其中被用于陆地研究的多达10种，常用的陆地产品如MOD09(地表反射率产品)、MOD13(地表植被指数产品)，这些MODIS陆地产品的时间分辨率为8-Day、16-Day和monthly合成值，空间分辨率可分为250 m、500 m、1 km 和 25 km［24］。而中国的农业生产往往以旬(10 d)为时间间隔，因此，NASA所提供的陆地产品在时间间隔上不能够满足中国农作物长势监测的需求。同时，NASA所提供的MOD09产品中各个波段不具有统一的最高250 m的地面分辨率，若想使用该产品需要对数据进行重采样［25］。因此，综合上述原因，本研究选用GSFC提供的MODIS原始数据，也被称为 Level-0(L0)级数据。自2000年2月至2012年11月间，覆盖研究区的L0级数据共计约9 000景，通过快视图筛选出667景无云且无太阳耀斑的数据用于本研究，并对所选取的遥感影像做如下处理:(1)辐射校正。通过辐射校正可消除因传感器自身条件、大气因素、太阳角度及其他噪声引起的卫星观测值与实际反射率之间的差异。因此，本研究采用NASA提供的水色遥感软件Sea-DAS 6.1以及MODIS Science Team提供的针对气溶胶吸收和瑞利散射校正软件对MODIS L0数据进行辐射校正，并将所有MODIS数据归一化为地面反射率。(2)几何校正。Terra和Aqua两颗卫星搭载有外部定位系统，每景MODIS L0数据的HDF文件中均包含按照顺序存放的扫描数据帧、时间码和定位信息，利用SeaDAS 6.1软件可以针对MODIS L0数据中包括的定位信息完成MODIS数据的几何校正。应用该软件包将所选研究区的数据校正为等经纬度切圆柱投影，坐标系为WGS-84，几何校正的误差小于0.5个像元［26］。(3)数据重采样。由于针对陆地设计的MODIS数据1～7波段地面分辨率不统一，即1～2波段为250 m地面分辨率，而3～7波段为500 m地面分辨率。因此，需要对3～7波段进行数据重采样，将其分辨率转换成250 m。考虑到双线性内插法计算量适中且精度较高，故本研究选用双线性内插法对500 m分辨率波段进行重采样。(4)图像裁剪。以泗洪县行政边界为基础，结合2011年3月31日10 m分辨率的ALOS遥感数据，去除泗洪县内的水体、河流、居民点及主要道路，生成用于批量裁切的研究区矢量文件并完成研究区MODIS数据的裁剪［27］。

以上数据处理过程均通过批处理模式完成，大大提高了计算效率，生成单景MODIS数据地面反射率产品仅需5 min，使海量遥感数据用于县级尺度作物长势实时监测成为可能。

1.4 作物长势监测方法

目前，作物长势监测的主要方法有过程监测和实时监测。过程监测主要是通过时序植被指数来构建作物生长过程，通过生长过程的年际间(如高产年、平产年、低产年和去年)对比来反映作物生长的状况，寻找出当年与典型年份曲线间的相似和差异，统计生长过程曲线的上升速率、下降速率、累计值等特征参数来反映当年作物生长趋势上的差异［28］。实时监测是利用实时遥感图像判断具体时间段(旬、月、年)作物长势优劣及其空间分布，评价该时间段内作物长势状况的空间分布，主要是将实时植被指数与去年、多年同期平均以及指定年份同期对比，通过对差异值进行分级、统计来反映区域作物实时生长状况［3，29-30］。过程监测和实时监测都要通过植被指数实现。植被指数是指由多光谱数据，经线性或非线性组合构成的各种数值，其对植物的叶绿素含量、健康程度以及植物含水量等特性非常敏感，因此，是评价植被覆盖、生长活力及生物量的简单有效的度量参数［31］。目前，被广泛使用的植被指数有:归一化植被指数NDVI(Normalized difference vegetation index)和增强型植被指数 EVI(Enhanced vegetation index)［24］，各个植被指数的计算公式如下:

NDVI因无法去除大气干扰(如水汽、瑞利散射、气溶胶等)，其值往往较实际值偏低，同时NDVI在高植被覆盖区域红光波段会出现饱和现象，造成其值在高植被覆盖区域无法与植被覆盖度同步增长。而EVI中增加了大气修正参数C1和C2，一定程度上减弱了大气干扰［24，33-34］，因此，本研究采用 EVI来实现对泗洪县近12年来每旬作物长势的实时监测和过程监测。

2 结果与分析

2.1 作物长势过程监测

将经过预处理的641景MODIS地面反射率数据(250 m分辨率)先按照公式(3)计算EVI值，之后，将每景图像中各像元的EVI值进行加和平均以获取每景图像的EVI平均值。根据短期植被生长过程中EVI表现为稳定的上升或者下降这一特性［35］，利用线性插值法生成泗洪县12年每日平均EVI值序列，共有4 538个平均EVI值，最后以旬为时间间隔按照对应年份及月份计算12年来泗洪县各月份逐旬平均EVI值。

统计数据显示，2011年泗洪县粮食作物产量为9.452×105t，为2000～2011年间的极大值。因此，选取2011年各旬EVI值作为典型丰产年参考序列。同时，选取2000～2012年各旬EVI平均值作为多年平均参考序列(图2)。将2012年监测期内各旬平均EVI值序列与典型丰产年及多年平均参考序列进行对比。结果表明，2011～2012年稻麦轮作期内，除4月中旬外其余各旬小麦长势均好于典型丰产年及多年平均。2011～2012年稻麦轮作期内，除6月中旬外其余各旬水稻长势均好于多年平均，但8月中旬～9月下旬以及11月上旬水稻长势不及典型丰产年。野外实地考察及历史资料均证明监测期内泗洪县的主要作物为小麦和水稻。因此，可以预测2012年泗洪县夏粮产量(小麦)将高于典型丰产年和多年平均;秋粮(水稻)产量将高于多年平均，但由于水稻产量形成的关键期内(8月中旬～9月下旬)其长势不及典型丰产年，所以秋粮的最终产量将不及典型丰产年。

图2 泗洪县2011～2012年稻麦轮作期各旬与典型丰产年及多年平均作物长势对比曲线图Fig.2 Comparison of the 10-day crop growth condition(CGC)of Sihong county among the 2011 －2012 wheat-rice rotation periods，the typical harvest year and the mean of the past 12 years

2.2 作物长势实时监测

以上述作物长势过程监测中所生成的667景EVI指数图为基础，生成2011～2012年稻麦轮作期内各旬平均EVI指数图，同时生成典型丰产年以及多年平均(2000～2012年)各旬平均EVI指数图。将2011～2012年稻麦轮作期内各旬的EVI指数图与同期的典型丰产年以及多年平均EVI指数图进行比较，并按照EVI差值≥0.1、0.1＞EVI差值＞ －0.1、EVI差值≤－0.1将长势划分为优于典型丰产年(或多年平均)、与典型丰产年(或多年平均)持平、劣于典型丰产年(或多年平均)3个级别(图3～6)。

通过统计对比，发现2011～2012年稻麦轮作期内各旬小麦长势劣于典型丰产年的面积至多占研究区总面积的13%，劣于多年平均的面积至多占研究区总面积的10%，这说明该时间段内泗洪县小麦长势大多为优于典型丰产年(或多年平均)或持平。而2011～2012年稻麦轮作期内各旬水稻长势劣于典型丰产年的面积至多占研究区总面积的48%，劣于多年平均的面积至多占研究区总面积的35%，这说明该时间段内泗洪县水稻长势较多为优于典型丰产年(或多年平均)或持平。

泗洪县东北部及西南部地势较高，耕地以旱田为主;而泗洪县中部地势较低，耕地以水旱轮作田为主。由于水旱轮作区冬小麦播种较旱田晚，使得该区域的小麦在拔节期(4月上旬)之前长势要劣于旱田区域，这一结论也得到野外实地验证的证明。但水旱轮作区域的水分及土壤肥力条件均要优于旱田区域，保证了该区域的小麦在4月下旬之后的孕穗期、抽穗期以及成熟期中长势基本与旱田区域生长周期更长的小麦长势相当(图3、图4、图7)。随着5月下旬小麦相继收割，泗洪县的小麦长势开始出现退化，但从图7中的5月下旬～6月上旬的监测结果中不难发现，泗洪县的旱田区域小麦收割进程要先于水旱轮作区域。可见，通过对作物长势的实时监测，不仅能够及时反馈各旬作物长势的空间分布情况，同时也能够做到对作物生长及收获进程的全程监控(图7、图8)。

图3 泗洪县2011～2012年稻麦轮作期小麦长势与典型丰产年同期对比分级图Fig.3 Comparison of the 10-day wheat growth condition of Sihong county between the 2011 －2012 wheat-rice rotation periods and the typical harvest year

图4 泗洪县2011～2012年稻麦轮作期小麦长势与多年平均对比分级图Fig.4 Comparison of the 10-day wheat growth condition of Sihong county between the 2011 －2012 wheat-rice rotation periods and the mean of the past 12 years

图7 泗洪县2011～2012年稻麦轮作期内小麦长势分布图Fig.7 The 10-day wheat growth condition of Sihong county during the 2011 －2012 wheat-rice rotation periods

图8 泗洪县2011～2012年稻麦轮作期内水稻长势分布图Fig.8 The 10-day rice growth condition of Sihong county during the 2011 －2012 wheat-rice rotation periods

2.3 作物长势监测精度验证

为了保证作物长势监测结果的准确性，随机选取均匀分布于研究区内的32个验证点对2012年4月上旬的小麦和2012年9月下旬的水稻长势遥感监测结果进行了野外实地验证(图9)，验证结果表明遥感监测结果能够准确、客观地反映研究区内作物长势的差异。

图9 泗洪县作物长势遥感监测结果及野外验证信息表Fig.9 The CGC monitoring results vs.the ground monitoring results

3 结论

泗洪县的作物长势过程监测结果表明，2011～2012年稻麦轮作期内小麦长势优于典型丰产年以及多年平均;水稻长势优于多年平均，但在关键生长期内不及典型丰产年。而作物长势实时监测结果不仅显示了监测周期内作物长势优劣在空间上的分布情况，同时也揭示了作物生长及收获的变化过程。因此，二者结合监测作物长势，可使地方农业管理部门足不出户就能全面掌握辖区内的作物生长及收割情况，以便及时调整田间管理措施，实现防灾减灾、增产增收的目的。

图5 泗洪县2011-2012年稻麦轮作期水稻长势与典型丰产年同期对比分级图Fig.5 Comparison of the 10-day rice growth condition of Sihong county between the 2011 －2012 wheat-rice rotation periods and the typical harvest year

图6 泗洪县2011～2012年稻麦轮作期水稻长势与多年平均对比分级图Fig.6 Comparison of the 10-day rice growth condition of Sihong county between the 2011 －2012 wheat-rice rotation periods and the mean of the past 12 years

将近12年来泗洪县的平均EVI值与对应年份的粮食产量进行了相关性分析，二者的相关系数为0.66，在置信度水平为95%下表现为显著相关。这表明作物长势在一定程度上反映了粮食产量。因此，如果将定期的作物长势监测与生长模型以及气象数据结合，可以对研究区内的作物产量进行估测。

由于受季风影响，江苏省夏初常伴随持续阴雨天气，这大大增加了以旬为时间间隔的作物长势监测难度。因此，在今后的研究中需开发针对海量MODIS数据的有效自动云像元识别算法，并用与之最近日期无云像元替换对应的云像元，从而实现对监测目标的不间断实时监测。

MODIS数据时间分辨率高且可以免费获取，因此，如果能将利用250 m分辨率的MODIS数据实时监测县级尺度作物长势方法推广到江苏省其他县市，既能满足地方农业管理部门对实时作物长势信息的需求，也将有助于当地的农业防灾、减灾工作。

［1］吴炳方，蒙继华，李强子.国外农情遥感监测系统现状与启示［J］.地球科学进展，2010，25(10):1003-1012.

［2］杨邦杰，裴志远，周清波，等.我国农情遥感监测关键技术研究进展［J］.农业工程学报，2002，3(18):191-194.

［3］吴炳方，张 峰，刘成林，等.农作物长势综合遥感监测方法［J］.遥感学报，2004，8(6):498-514.

［4］MACDONALD R B，HALL F G.Global crop forecasting［J］.Science，1980，208:670-679.

［5］蒙继华，吴炳方，李子强，等.集成化的省级农情遥感监测系统［J］.农业工程学报，2011，27(6):169-175.

［6］USDA FAS.The FAS cropexplorer:A web success story［EB/OL］.(2005-06-01)［2013-03-01］.http://www.fas.usda.gov/info/fasworldwide/2005/06-2005/CropExplorer.pdf.

［7］USDA FAS.Production，Supply and distribution online［EB/OL］.(2007-08-12)［2013-03-01］.http://www.fas.usda.gov/psdonline.

［8］USDA NASS.History of remotesensing for crop acreage［EB/OL］.(2009-05-27)［2013-03-01］.http://www.nass.usda.gov/Surveys/Remotely_Sensed_Data_Crop_Acreage.

［9］STATISTICSCANADA.Overview of the crop condition assessment program ［EB/OL］.(2012-05-07)［2013-03-01］.http://www26.statcan.ca/ccap/overview-aper-cu-eng.jsp.

［10］CONAB.GEOSAFRAS［EB/OL］.(2010-06-18)［2013-03-01］.http://www.conab.gov.br/conabweb/index.php?PAG=81.

［11］陈述彭.遥感在农业科学技术中的应用［M］.北京:科学出版社，1990.

［12］肖淑昭，孟宪钺，张桂宗，等.NOAA/AVHRR资料在中小尺度地区进行冬小麦估产的应用研究［J］.环境遥感，1988，3(3):299-307.

［13］肖乾广.用NOAA气象卫星的定量资料计算冬小麦种植面积的两种方法［J］.环境遥感，1989，4(3):191-196.

［14］李郁竹.北方冬小麦卫星遥感监测及估产业务系统［J］.气象，1992，18(11):14-16.

［15］吴炳方.全国农情监测与估产的运行化遥感方法［J］.地理学报，2000，55(1):25-35.

［16］吴炳方.中国农情遥感监测研究［J］.中国科学院院刊，2004，19(3):202-205.

［17］吴素霞，毛任钊，李红军，等.中国农作物长势遥感监测研究综述［J］.中国农学通报，2005，21(3):319-322.

［18］江苏省统计局等.江苏省农村统计年鉴［M］.北京:中国统计出版社，2009.

［19］江苏省统计局等.江苏省农村统计年鉴［M］.北京:中国统计出版社，2010.

［20］陆一忠，刘春山.2003年洪泽湖洪水分析［J］.江苏水利，2004(7):17-18.

［21］楚恩国，卜贤晖.2004年洪泽湖流域干旱原因分析［J］.水文，2006，26(5):80-82.

［22］刘玉洁，杨忠东.MODIS遥感信息处理原理与算法［M］.北京:科学出版社，2001.

［23］于 堃，陆殿梅，熊黑钢.近7年来渤海海区冬季表层海水叶绿素浓度的遥感反演及其变化分析［J］.遥感信息，2009(6):55-62.

［24］HUETE A R，JUSTICE C.MODIS vegetation index(MOD13)algorithm theoretical basis document［R/OL］.(1999-04-30)［2013-03-01］.http://modis.gsfc.nasa.gov/data/atbd/atbd_mod13.pdf.

［25］程 乾.基于MOD09产品的水稻叶面积指数和叶绿素含量的遥感估算模型［J］.应用生态学报，2006，17(8):1453-1458.

［26］WOLFE R E，NISHIHAMA M，FLEIG A J，et al.Achieving subpixel geolocation accuracy in support of MODIS land science［J］.Remote Sensing of Environment，2002，83:31-49.

［27］于 堃.典型平原湿地成因及近10年来植被变化研究［D］.南京:南京大学，2011.

［28］GENOVESE G，VIGNOLLES C，NEGRE T，et al.A methodology for a combined use of normalized difference vegetation index and CORINE land cover data for crop yield monitoring and forecasting-A case study on spain［J］.Agronomie，2001，21(1):91-111.

［29］武建军，杨勤业.干旱区农作物长势综合监测［J］.地理学报，2002，21(5):593-598.

［30］刘爱霞，王长耀，刘正军，等.基于RS与GIS的干旱区棉花信息提取及长势监测［J］.地理与地理信息科学，2003，19(4):101-104.

［31］SIVANPILLAI R，LATCHININSKY A V.Mapping locust habitats in the Amudarya River Delta，Uzbekistan with multi-temporal MODIS imagery［J］.Environmental Management，2007，39:876-886.

［32］VERMOTE E F，TANRE J L，DEUZE M H，et al.Second simulation of the satellite signal in the solar spectrum，6S:an overview［J］.IEEE Transactions Geoscience and Remote Sensing，1997，35:675-686.

［33］HUETE A R，DIDAN K，MIURA T，et al.Overview of the radiometric and biophysical performance of the MODIS vegetation indices［J］.Remote Sensing of Environment，2002，83:95-213.

［34］王正兴，刘 闯，HUETE A R.植被指数研究进展:从AVHRRNDVI到 MODIS-EVI［J］.生态学报，2003，23(5):979-986.

［35］PAN Y，LI L，ZHANG J，et al.Winter wheat area estimation from MODIS-EVI time series data using the Crop Proportion Phenology Index［J］.Remote Sensing of Environment，2012，119:232-242.