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.