基于物候信息的山東省冬小麥長勢遙感監測

侯學會， 隋學艷， 姚慧敏， 梁守真， 王 猛

(1.山東省農業可持續發展研究所,濟南 250100； 2.農業部華東都市農業重點實驗室,濟南 250100)

0 引言

作物長勢監測是農情監測最核心的內容之一。作物長勢信息客觀描述了作物生長、土壤墑情、肥力及植物營養的綜合狀況，是作物估產的必要前提。遙感技術因其可獲得在空間和時間上的連續信息，已成為大范圍作物長勢監測的有效手段[1]。國內外很多組織建立了政府性的農作物遙感監測系統，為相關部門及時提供農情信息[2]。

植被指數，如歸一化差值植被指數(normalized difference vegetation index，NDVI)、增強型植被指數(enhanced vegetation index，EVI)和歸一化差值紅外指數(normalized difference infrared index，NDII)等，是利用遙感技術研究植被長勢的主要指標[3-4]。國內外學者基于遙感技術進行作物長勢監測已經取得了一系列成果。黃青等[5]基于MODIS NDVI 數據，通過與近5 a 作物長勢的平均狀況對比，發現2013年東北地區不同作物長勢在整個生育期內的時空分布存在較大差異； 鄒文濤等[6]利用MODIS NDVI 序列數據，綜合采用實時監測、過程監測和時間序列聚類監測方法，以2008―2012年間平均長勢為基礎，對印度2013年作物長勢進行監測，優化了Cropwatch作物長勢監測方法； Esquerdo等[7]、Duveiller等[8]和Becker-Reshef等[9]也基于SPOT-VGT，NOAA-AVHRR和MODIS遙感影像開展了大范圍的作物長勢監測。但綜合來看，目前已有的成果一般都是根據作物多年平均長勢狀況對當年相應時間點或時間段的作物長勢進行研究，而在大范圍的作物長勢監測時，不同區域的作物會因所處物候期的不同導致長勢出現差異； 且受氣候因子和人為因素的影響，同一區域作物的物候在時間序列上也會發生變化，這種因為物候變化產生的長勢差異與作物本身長勢狀況的變化混在一起，增加了對作物長勢在時間序列上變化分析的不確定性。因此，開展大區域、長時間序列作物長勢監測研究，首先需要對作物長勢參數進行物候修正，以消除或減少物候差異對作物長勢監測結果的影響。

本文以山東省冬小麥為研究對象，以MODIS 8 d合成反射率數據為主要數據源，結合地面調查數據，在提取2001―2015年間冬小麥抽穗期的基礎上，進行抽穗期冬小麥長勢監測； 并對冬小麥長勢狀況的時空格局進行分析，初步探索如何消除長時間和大區域作物長勢監測中物候差異的影響。

1 數據源與預處理

1.1 遙感數據源

本文使用的遙感數據為8 d最大值合成的MODIS反射率數據集——MOD09A1，空間分辨率為500 m，從MODIS數據免費分發網站(https: //ladsweb.nascom.nasa.gov/search/)獲得，數據覆蓋每年的3―5月，年際范圍為2001―2015年。MOD09A1的原始數據為HDF格式的Sinusoidal投影，利用MRT軟件對MOD09A1進行預處理，提取植被指數構建所需的4個波段(B1―3和B6)的反射率和質量標記文件refl_500m_qc，并轉為tiff格式，投影到Albers坐標系下。根據表1，基于反射率數據計算得到NDVI，EVI，NDII，產品改進-NDVI(product improve-NDVI,PI_NDVI)和產品改進-EVI(product improve-EVI,PI_EVI)等 5種常用的植被指數； 然后基于像元質量文件，利用Savitzky-Golay濾波法對時序植被指數進行像元尺度上的去噪處理[10-11]，進一步去除影像中的噪聲，提高監測精度。為了與驗證數據使用的儒略日一致，將濾波之后的8 d合成植被指數重采樣為1 d的數據。

表1 各種植被指數的定義Tab.1 Definitions of several vegetation indices

表1中，Rnir，Rred，Rblue和Rswir分別為MOD09A1數據中B1,B2,B3和B6波段的反射率值；L為冠層背景的調整因子;C1和C2分別為權重系數，用于減少大氣氣溶膠影響。文中L,C1和C2分別取值為1.0，6.0和7.5[13]。

1.2 冬小麥種植區

本文中的冬小麥種植區依據中科院地理所發布的2005年山東省土地利用數據中的耕地數據進行掩模運算得到(圖1)。

圖1 山東省冬小麥種植區分布及農業臺站位置

因為冬小麥種植面積變化不是本文研究的重點，因此本文未對研究時間段內冬小麥種植區的分布變化對研究精度的影響進行分析。

1.3 地面驗證數據

利用中國氣象科學數據共享網(http: //data.cma.cn/)獲得的山東省11個站點(圖1)的冬小麥2001―2013年間的物候數據，對提取的研究時間段內山東省冬小麥抽穗期進行驗證。冬小麥地面物候觀測數據為“年/月/日”格式，為了與利用遙感數據提取的第N天數據對應，利用Excel中days360()函數將年/月/日轉換成儒略日。

2 研究方法

2.1 冬小麥抽穗期數據提取

基于遙感技術的作物物候提取方法已有很多，本文基于濾波重構后的逐日植被指數時間序列提取植被指數(vegetation index，VI)峰值出現的時間，作為冬小麥抽穗期的時間標識。NDVI是目前基于遙感進行作物物候提取的最常用指標，但在植被生長旺盛時期，NDVI易出現飽和； 而EVI因加入了藍光波段，可降低土壤背景和大氣因素的影響，消除NDVI易飽和的缺陷[16]。 因此，本文基于NDVI和EVI時序數據，分別進行山東省冬小麥抽穗期提取； 然后根據11個農業臺站觀測的冬小麥抽穗期數據，依據相關系數r、平均絕對誤差MAE、偏差bias和納什系數NSE，對2種指數的提取結果進行判定，選擇提取精度較高的指標作為本文冬小麥抽穗期提取的遙感指標。

2.2 冬小麥長勢狀況判定

首先，逐像元計算出2001―2015年間抽穗期的NDVI，EVI，NDII，PI_NDVI和PI_EVI值； 然后，計算所有年份抽穗期中各指數均值； 最后，提取各指數的逐年距平值(即指數值與指數均值的差值)。正距平指數表示植被長勢比多年平均狀況好，負距平指數則表示植被長勢比多年平均狀況差。陳維英等[17]基于平均植被指數(average VI，AVI)對干旱狀況進行了分級，(-0.2,-0.1)表示干旱，(-0.6,-0.3)表示嚴重干旱。以此為參考，結合地面調查資料，本文將山東省冬小麥長勢分為5個等級。由于缺少直接表征冬小麥長勢的地面驗證數據，據已有研究，與其他生長期相比，小麥抽穗期的綜合長勢情況與產量呈顯著相關[18]。因此，為了評價各植被指數在山東省小麥長勢監測中的適應性，首先，計算出研究區每年的各指數和； 然后依據山東省統計局發布的2001―2014年間山東省冬小麥產量數據，評價各指數監測山東省冬小麥長勢精度的適宜性； 從而選擇出能表征山東省小麥長勢狀況的遙感指標，并對近15 a以來的山東省抽穗期冬小麥長勢進行評價。

山東省冬小麥長勢等級劃分如表2所示。

表2 山東省冬小麥長勢劃分等級Tab.2 Grade of growth of winter wheats inShandong Province

3 結果分析與驗證

3.1 2001―2015年山東省冬小麥抽穗期分析

3.1.1 冬小麥抽穗期遙感監測結果驗證

根據經緯度分別提取農業臺站站點的冬小麥抽穗期NDVI和EVI監測結果，并利用站點數據進行驗證，結果如圖2所示。

(a) NDVI指數 (b) EVI指數

圖2冬小麥抽穗期遙感提取結果與地面實測數據對比

Fig.2Comparisonbetweenheadingdatesofwinterwheatextractedusingremotesensingdataandgroundmeasureddata

基于NDVI和EVI提取的冬小麥抽穗期與地面實測數據均呈極顯著相關關系(NDVI：r=0.469，p<0.001； EVI：r=0.559，p<0.001)，說明二者均可作為冬小麥抽穗期遙感判識指標。分析二者提取結果的其他評價指標，基于EVI提取的MAE為9.45，小于基于NDVI提取結果(MAE=10.1)； 且基于EVI提取結果的bias和NSE分別為0.06和-2.39，優于基于NDVI提取結果的0.07和-2.64。說明本文基于EVI指數提取的冬小麥抽穗期優于基于NDVI指數提取的冬小麥抽穗期，可以作為研究區冬小麥抽穗期長勢評價的依據。

3.1.2 近15 a冬小麥抽穗期時空格局

圖3分別示出基于EVI提取的山東省近15 a以來冬小麥抽穗期均值的空間分布(圖3(a))和時間變化趨勢(圖3(b))。根據圖3(a)可知，山東省冬小麥抽穗期主要集中在年積日第100―120天(即4月中旬―下旬)，并從南向北、自西向東逐漸推遲； 從時間變化(圖3(b))來看，近15 a以來山東省冬小麥抽穗期的時間變化趨勢并不明顯，變化在-0.5～0.5 d/a的區域占冬小麥種植區總面積的63%以上，但僅有0.07%區域的冬小麥抽穗期變化達到顯著水平。

(a) 冬小麥抽穗期均值 (b) 冬小麥抽穗期時間變化

圖32001—2015年山東省冬小麥抽穗期時空格局

Fig.3Spatio-temporalpatternofheadingdatesofwinterwheatinShandongProvinceduring2001—2015

3.2 2001―2015年山東省冬小麥長勢評價

3.2.1 冬小麥長勢時空格局

表3示出研究區各植被指數值逐像元累積之和與2001―2014年間山東省冬小麥總產量數據的相關關系。

表3 各植被指數表征的長勢與產量相關關系Tab.3 Correlation between VIs and yield of winter wheat

根據表3，綜合考慮作物綠度和水分指數(PI_NDVI和PI_EVI)與作物產量的相關性達到顯著水平，優于只考慮綠度的指數(NDVI和EVI)和只考慮水分的指數(NDII)，其中以PI_NDVI的相關關系最優，為極顯著相關(r=0.727，p=0.003)。因此，本文以PI_NDVI為指標，對近15 a以來的山東省冬小麥長勢提取結果進行評價。

為更好地反映山東省冬小麥長勢演變，本文取研究區每年PI_NDVI的均值代表該年山東省冬小麥整體長勢，結果如圖4所示。從圖4可以看出，山東省冬小麥抽穗期間長勢在2001―2015年間整體呈上升趨勢，年變化趨勢為0.006(r=0.442)； 年際間波動較大，如2004和2014年全省小麥抽穗期時長勢明顯優于其他年份，2004和2014年PI_NDVI的距平值分別為0.214 5和0.209 3，分別比近15 a平均距平值高218.25%和210.53%，這可能主要是因為2004年2月份[19]和2014年2月份[20]山東全省降水偏多，且氣溫偏高，為冬小麥越冬期之后的生長提供了適宜的生長環境，使小麥長勢較常年好。其余年份長勢狀況差別不明顯，比較接近于多年平均狀況(PI_NDVI=0.067 4)。

圖4 2001—2015年山東省冬小麥抽穗期

圖5為2001—2015年山東省冬小麥抽穗期長勢的空間分布。

(a) 2001年 (b) 2002年 (c) 2003年(d) 2004年

(e) 2005年 (f) 2006年 (g) 2007年

圖5-12001—2015年山東省冬小麥抽穗期長勢空間分布

Fig.5-1Distributionofgrowthconditionofwinterwheatduringheadingdatesfrom2001to2015inShandongProvince

(h) 2008年 (i) 2009年 (j) 2010年(k) 2011年

(l) 2012年 (m) 2013年 (n) 2014年(o) 2015年

圖5-22001—2015年山東省冬小麥抽穗期長勢空間分布

Fig.5-2Distributionofgrowthconditionofwinterwheatduringheadingdatesfrom2001to2015inShandongProvince

從總體來看，近15 a山東省冬小麥長勢很差的區域所占比重很小，2002年占比0.32%、2005年占比0.89%、2008年占比0.51%、2010年占比0.69%、2013年占比1.07%。部分區域的小麥長勢嚴重低于多年均值，主要分布在魯西地區； 但2013年受春季大范圍干旱的影響，小麥長勢很差的區域范圍較廣，魯西、魯中和魯南地區均有出現，長勢較差的區域占當年冬小麥種植面積的20.36%。但從總體上來看，在研究時間段內，絕大部分區域的冬小麥長勢都比較一致，長勢與15 a間平均長勢持平的區域占63%以上，其中以2009年所占比重最大(達 92.94%)，而2013年所占比重僅為63.56%。雖然2015年冬小麥整體長勢不如2014年(圖4)，但2015年冬小麥長勢好于多年平均的區域所占比重最大(達27.85%)。

從空間分布來看，魯西南地區冬小麥長勢在整體上比較一致，除2001年和2002年魯西南部分地區小麥長勢較差外，其余年份的小麥長勢比較一致； 而魯西北地區的小麥長勢總體上呈上升趨勢。這可能是因為魯西地區是山東省主要農業區，政府在該地區對農田水利設施的投入較多，基本上能保障農田需水； 而魯東地區冬小麥長勢狀況在年際間波動較大，這與該地區作物主要靠雨水供養有關，降水較少的年份小麥長勢明顯差于正常年份。

3.2.2 評價結果驗證

冬小麥長勢是降水、墑情/干旱狀況的綜合反映。本文利用已有的相關研究對山東省冬小麥長勢評價結果進行驗證。據楊麗萍等[21]的研究，2008年4月上旬山東全省平均降水量為15.9 mm，比常年偏多148%； 但魯西北大部、半島及魯中地區平均降水量在10 mm以下，墑情較低。而基于遙感數據和地面實測數據反演結果顯示，2008年4月上旬山東省大部分地區墑情適宜，魯北、魯中及半島地區土壤墑情低，出現了一定的干旱，魯北地區旱情尤為嚴重，這與本文結果基本一致。另外，劉暢等[22]基于氣象觀測資料分析指出，2014年春季在山東省冬小麥主產區的魯西北和魯南大部分地區降水較常年偏多，且氣溫較常年較高，日照時數偏多，導致冬小麥生長期進程加快； 段海霞等[23]的研究指出2014年4月上旬山東省出現不同程度的旱情，但到4月下旬旱情已得到緩解，僅魯中和魯西的部分地區出現輕微干旱。而本文研究得出2014年山東省冬小麥92%以上區域長勢好于多年平均或持平，僅魯中和魯西的部分地區呈現長勢較差狀況，這與劉暢等人基于氣候資料分析的結果比較一致。據2012年山東省氣候影響評價，2012年初夏山東省農田干旱總面積約為120萬hm2，其中重旱面積12.4萬hm2。而本文基于MOD09A1數據反演冬小麥長勢，得到2012年抽穗期冬小麥比多年平均長勢差的面積為149.14萬hm2，長勢較差區域面積為7.47萬hm2。盡管因研究空間尺度、數據源和監測時間以及統計對象的差異導致本文研究結論與已有成果在空間特性上略有差異，但本文結論與前人成果仍具有一定的可比性。因此，本文的監測結果基本上反映了研究區冬小麥的長勢狀況，具有一定的可靠性。

4 結論

本文基于MODIS 8 d最大值合成反射率數據，在提取冬小麥抽穗期的前提下，分析了2001―2015年間山東省冬小麥抽穗期時長勢狀況。

1)用加入藍光波段構建的EVI指數提取的山東省冬小麥抽穗期具有一定的可靠性。冬小麥抽穗期主要集中在4月中、下旬，且從南向北、從西向東逐漸推遲。

2)遙感監測的冬小麥長勢是植被綠度和水分綜合信息的反映，基于PI_NDVI指數，近15 a間山東省冬小麥長勢狀況整體呈上升趨勢； 但年際間波動較大，且因受自然氣候因素和人工田間管理的影響，相同年份、不同區域的冬小麥長勢存在顯著差異，但大部分區域與多年狀況持平。

3)與自然植被不同，冬小麥長勢狀況受自然氣候因素和田間管理綜合影響，尤其是近年來山東省加大農田基本設施建設，在平原地區基本上能保證農田需水，因此導致基于自然氣候因素分析農田作物農情的精度不高。如何結合田間管理提高農情遙感監測精度，需要進一步深入研究。

參考文獻(References)：

[1] 蒙繼華,杜 鑫,張 淼,等.物候信息在大范圍作物長勢遙感監測中的應用[J].遙感技術與應用,2014,29(2):278-285.

Meng J H,Du X,Zhang M,et al.Integrating crop phenophase information in large-area crop condition evaluation with remote sensing[J].Remote Sensing Technology and Application,2014,29(2):278-285.

[2] 吳炳方,蒙繼華,李強子.國外農情遙感監測系統現狀與啟示[J].地球科學進展,2010,25(10):1003-1012.

Wu B F,Meng J H,Li Q Z.Review of overseas crop monitoring systems with remote sensing[J].Advances in Earth Science,2010,25(10):1003-1012.

[3] 趙 虎,楊正偉,李 霖,等.作物長勢遙感監測指標的改進與比較分析[J].農業工程學報,2011,27(1):243-249.

Zhao H,Yang Z W,Li L,et al.Improvement and comparative analysis of indices of crop growth condition monitoring by remote sensing[J].Transactions of the Chinese Society of Agricultaual Engineering,2011,27(1):243-249.

[4] Guindin-Garcia N,Gitelson A A,Arkebauer T J,et al.An evaluation of MODIS 8- and 16-day composite products for monitoring maize green leaf area index[J].Agricultural and Forest Meteorology,2012,161:15-25.

[5] 黃 青,唐華俊,周清波,等.東北地區主要作物種植結構遙感提取及長勢監測[J].農業工程學報,2010,26(9):218-223.

Huang Q,Tang H J,Zhou Q B,et al.Remote-sensing based monitoring of planting structure and growth condition of major crops in Northeast China[J].Transactions of the Chinese Society of Agricultaual Engineering,2010,26(9):218-223.

[6] 鄒文濤,吳炳方,張 淼,等.農作物長勢綜合監測——以印度為例[J].遙感學報,2015,19(4):539-549.

Zou W T,Wu B F,Zhang M,et al.Synthetic method for crop condition analysis:A case study in India[J].Journal of Remote Sensing,2015,19(4):539-549.

[7] Esquerdo J C D M,Zullo Júnior J,Antunes J F G.Use of NDVI/AVHRR time-series profiles for soybean crop monitoring in Brazil[J].International Journal of Remote Sensing,2011,32(13):3711-3727.

[8] Duveiller G,López-Lozano R,Baruth B.Enhanced processing of 1-km spatial resolution fAPAR time series for sugarcane yield forecasting and monitoring[J].Remote Sensing,2013,5(3):1091-1116.

[9] Becker-Reshef B,Vermote E,Lindeman M,et al.A generalized regression-based model for forecasting winter wheat yields in Kansas and Ukraine using MODIS data[J].Remote Sensing of Environment,2010,114(6):1312-1323.

[10] Chen J,Jonsson P,Tamura M,et al.A simple method for reconstructing a high-quality NDVI time-series data set based on the Savitzky-Golay filter[J].Remote Sensing of Environment,2004,91(3/4):332-344.

[11] 吳文斌,楊 鵬,唐華俊,等.兩種NDVI時間序列數據擬合方法比較[J].農業工程學報,2009,25(11):183-188.

Wu W B,Yang P,Tang H J,et al.Comparison of two fitting methods of NDVI time series datasets[J].Transactions of the Chinese Society of Agricultaual Engineering,2009,25(11):183-188.

[12] Rouse Jr J W,Haas R H,Schell J A,et al.Monitoring vegetation systems in the Great Plains with ERTS[C]//Proceedings of the 3rd Earth Resources Technology Satellite-1 Symposium. College Station,TX,United States:NASA,1974:3010-3017.

[13] Huete A,Justice C,Liu H.Development of vegetation and soil indices for MODIS-EOS[J].Remote Sensing of Environment,1994,49(3):224-234.

[14] Hardisky M A,Klemas V,Smart R M.The influence of soil salinity,growth form,and leaf moisture on the spectral radiance of Spartina alterniflora canopies[J].Photogrammetric Engineering and Remote Sensing,1983,49(1):77-83.

[15] Gonsamo A,Chen J M.Continuous observation of leaf area index at Fluxnet-Canada sites[J].Agricultural and Forest Meteorology,2014,189-190:168-174.

[16] Liu H Q,Huete A.A feedback based modification of the NDVI to minimize canopy background and atmospheric noise[J].IEEE Transactions on Geoscience and Remote Sensing,1995,33(2):457-465.

[17] 陳維英,肖乾廣,盛永偉.距平植被指數在1992年特大干旱監測中的應用[J].環境遙感,1994,9(2):106-112.

Chen W Y,Xiao Q G,Sheng Y W.Application of the anomaly vegetation index to monitoring heavy drought in 1992[J].Remote Sensing of Environment China,1994,9(2):106-112.

[18] 肖偉中.小麥長勢與產量信息空間變異性及相關性研究[D].鎮江:江蘇大學,2009.

Xiao W Z.Research on Spatial Variability and Relationships of Wheat Growth and Yield[D].Zhenjiang:Jiangsu University,2009.

[19] 耿 勃,張 颯.2003年冬季(2003年12月—2004年2月)山東天氣評述[J].山東氣象,2004,24(1):48-50.

Geng B,Zhang S.Weather review of the winter 2003(Dec.2003-Feb.2004) in Shandong Province[J].Journal of Shandong Meteorology,2004,24(1):48-50.

[20] 韓永清,劉 暢,王 娜,等.2013年冬季(2013年12月—2014年2月)山東天氣評述[J].山東氣象,2014,34(1):77-79.

Han Y Q,Liu C,Wang N,et al.A review on the weather of Shandong in the winter of 2013(December,2013—February,2014)[J].Journal of Shandong Meteorology,2014,34(1):77-79.

[21] 楊麗萍,隋學艷,楊 潔,等.山東省春季土壤墑情遙感監測模型構建[J].山東農業科學,2009(5):17-20.

Yang L P,Sui X Y,Yang J,et al.Construction of remote sensing monitoring model for spring soil moisture in Shandong Province[J].Shandong Agricultural Sciences,2009(5):17-20.

[22] 劉 暢,王 娜,孟祥新.2014年春季(2014年3—5月)山東天氣評述[J].山東氣象,2014,34(2)：62-64．

Liu C,Wang N,Meng X X.A review on the weather of Shandong in the spring of 2014 ( from March to May,2014)[J].Journal of Shandong Meteorology,2014,34(2)：62-64．

[23] 段海霞,王素萍,馮建英.2014年全國干旱狀況及其影響與成因[J].干旱氣象,2015,33(2):349-360.

Duan H X,Wang S P,Feng J Y.Drought events and its influence in 2014 in China[J].Journal of Arid Meteorology,2015,33(2):349-360.