基于MODIS NDVI數(shù)據(jù)的鄭州市植物物候監(jiān)測研究

楊 陽 何瑞珍 田國行 任曉娟 李 巖 鄭景飚

(1. 河南農(nóng)業(yè)大學(xué)林學(xué)院，河南 鄭州 450002；2. 上海工程技術(shù)大學(xué)電子電氣工程學(xué)院，上海 201620；3. 信息工程大學(xué)地理空間信息學(xué)院，河南 鄭州 450001)

植物由于受氣候等環(huán)境因子影響形成的生長發(fā)育節(jié)律稱為植物物候，因其能直接響應(yīng)自然氣候變化的季節(jié)、年際規(guī)律，被作為重要的生態(tài)指示器廣泛應(yīng)用于農(nóng)業(yè)生產(chǎn)、全球氣候變化、生態(tài)學(xué)實踐等[1]。以往的物候監(jiān)測一般采用目視觀察法，直接在野外定點觀察和測量生物物候的年內(nèi)、年際變化[2]。該觀測方法具有準(zhǔn)確性和客觀性的特點，但難以進(jìn)行長周期、大尺度的物候時空變化監(jiān)測和分析，同時會耗費大量的人力和物力。而日益發(fā)展的遙感技術(shù)則為研究物候提供了新的手段。

光譜曲線反應(yīng)了物體對波的吸收和反射狀況，其形狀取決于物體內(nèi)在的組織結(jié)構(gòu)特征。因此，不同物體各具不同的光譜曲線，而同類物體由于內(nèi)部構(gòu)造相近形成了相似的光譜曲線形狀，如植被特有的光譜曲線。遙感技術(shù)正是基于此原理實現(xiàn)對植物物候的監(jiān)測。遙感監(jiān)測覆蓋范圍廣、時間序列較長，同時滿足空間連續(xù)和多時相監(jiān)測，對于指示植被動態(tài)變化與全球氣候變化的關(guān)系具有重要意義[3]。

遙感植被監(jiān)測是利用遙感植被指數(shù)與研究區(qū)植被參數(shù)的相互關(guān)系，以植被指數(shù)的變化來反映植被狀況[4]。NDVI是基于MODIS數(shù)據(jù)的歸一化植被指數(shù)，其利用綠色植物的波譜特性能在大范圍覆蓋區(qū)域內(nèi)反映植被覆蓋度、光合有效輻射等植被參數(shù)。Zhang等[5]采用MODIS 時間序列數(shù)據(jù)提取美國東北部植被的物候期，并分析了不同植被物候期隨緯度變化的情況。劉玲玲等[3]利用2005年AVHRR和MODIS的NDVI數(shù)據(jù)研究了全球植被物候參數(shù)，并通過比較分析，證明2種傳感器提取的結(jié)果在空間變化趨勢上具有一致性。這些研究均表明NDVI能有效的評估植被生長狀態(tài)[6]，可應(yīng)用于監(jiān)測植被、分析物候等。

本研究結(jié)合鄭州市域植被類型的分布特點，基于MODIS NDVI時間序列數(shù)據(jù)進(jìn)行分析，采用滑動平均法獲取鄭州市域植被關(guān)鍵物候期，以期實現(xiàn)鄭州主要植被類型物候的實時監(jiān)測，達(dá)到快速獲取區(qū)域植被信息的目的，對評估中原地區(qū)的氣候變化及其對植物的影響，指導(dǎo)當(dāng)?shù)氐霓r(nóng)業(yè)生產(chǎn)，監(jiān)測生態(tài)環(huán)境變化等具有重要的理論和現(xiàn)實意義。

1 研究區(qū)概況

鄭州市位于東經(jīng)112°42′～114°14′，北緯34°16′～34°58′，地處中原，北臨黃河，西依嵩山，東南為廣闊的黃淮平原。全市總面積7 446.2 km2，市區(qū)面積1 010.3 km2，農(nóng)用地、園林地和草地面積4 562.74 km2。鄭州屬于北溫帶大陸性季風(fēng)氣候，春季時間為每年的3月下旬至5月下旬；夏季為5月下旬至9月初；秋季為9月初至11月初；冬季為11月初至次年3月下旬。夏季炎熱、降雨集中，秋高氣爽、日照充足，冬季寒冷、降雪少。年平均氣溫15.6 ℃，年平均降雨量為542 mm，全年日照時間約1 869.7 h。

本研究根據(jù)鄭州市域的植物分布情況，以林地、草地、農(nóng)田3種代表性植被類型為對象進(jìn)行研究。由于農(nóng)田受氣候、土壤等自然因素影響，不同地區(qū)的農(nóng)田會進(jìn)行各種作物輪種，出現(xiàn)了一年一熟、一年兩熟，甚至于一年多熟的情況[7]。鄭州的農(nóng)田基本上為一年兩熟，以旱地作物為主，包括玉米 (Zeamays)、小麥 (Triticumaestivum)、棉花 (Gossypiumherbaceum)、大豆 (Glycinemax) 等。本研究對于農(nóng)田的復(fù)雜情況，以大面積種植的冬小麥和夏玉米為主要研究對象，進(jìn)行植物物候的監(jiān)測。

2 材料與方法

2.1 數(shù)據(jù)來源

本研究采用的MODIS陸地產(chǎn)品是搭載于EOS/Aqua衛(wèi)星上的MODIS獲取的數(shù)據(jù)，由NASA的EOS數(shù)據(jù)中心提供的16 d合成數(shù)據(jù) (MYD13A1)。MYD13A1是陸地的標(biāo)準(zhǔn)數(shù)據(jù)產(chǎn)品，內(nèi)容主要為NDVI，空間分辨率為500 m，時間范圍為2011年全年。GlobCover為2009年全球陸地覆蓋數(shù)據(jù)，分辨率為300 m。選取MERIS傳感器在2009年1月1日至12月31日期間所接收的較高質(zhì)量的影像數(shù)據(jù)，來進(jìn)行圖像合成，其采用 “美國食品和農(nóng)業(yè)組織的地表覆蓋分類系統(tǒng)” 作為圖例生成標(biāo)準(zhǔn)。

2.2 時間序列諧波分析法

云和大氣能夠降低NDVI數(shù)值，是影響遙感影像最大的噪聲[8]。本研究采用時間序列諧波分析法 (HANTS) 對NDVI數(shù)據(jù)進(jìn)行濾波處理，其核心算法是傅里葉變換和最小二乘法擬合，即把時間序列波譜數(shù)據(jù)的復(fù)雜函數(shù)分解成若干正弦、余弦函數(shù)，通過疊加其中具有時間序列特征的多條曲線，對時序數(shù)據(jù)進(jìn)行平滑。該方法在保持原時間序列周期性的同時，又能使異常像元恢復(fù)正常[9-11]。在NDVI時間序列數(shù)據(jù)中，受大氣干擾及處理誤差影響的諧波頻率一般較高，而反映植被相關(guān)信息的諧波頻率一般較低。因此平滑重建較低頻率的諧波即可減弱噪聲影響，重構(gòu)NDVI時間序列曲線[12-13]，并真實反映曲線的周期性變化規(guī)律。

2.3 物候劃分

將鄭州植物物候分為3個關(guān)鍵期：開始期、成熟期 (抽穗期)、結(jié)束期。開始期指植物開始生長或綠度開始增加的時間，為NDVI值開始增大，NDVI時間序列曲線上升速率迅速增大處；成熟期為植物綠度和光合作用最大化的時間，NDVI值達(dá)到最大；結(jié)束期為植物綠度和光合作用迅速下降的時間，NDVI時間序列曲線下降速率變化最大處；生長季的長度為結(jié)束期和開始期之間的差值。對于農(nóng)作物一年兩熟的情況，根據(jù)對農(nóng)作物的辨別和輪作的規(guī)律，分時段進(jìn)行研究。

2.4 滑動平均法

采用滑動平均法進(jìn)行植物物候關(guān)鍵期的提取。其基本原理是利用滑動平均后的曲線與實際NDVI時間序列曲線相交來獲取植物物候。把實際NDVI曲線第1次超過滑動平均曲線的時間作為開始期，第1次低于滑動平均曲線的時間作為結(jié)束期[14]，實際曲線最大值處對應(yīng)的日期即為成熟期。該方法對NDVI時間序列數(shù)據(jù)的計算更為穩(wěn)定、可靠，但其準(zhǔn)確性與滑動平均時間間隔的選取密切相關(guān)，而滑動平均時間間隔又與不同植被類型特征有關(guān)，需依研究植被實際狀況而定。

3 結(jié)果與分析

3.1 鄭州市域土地分類結(jié)果與驗證

對GlobCover數(shù)據(jù)進(jìn)行裁剪得到鄭州市域地表分類數(shù)據(jù)，通過地表覆蓋分類系統(tǒng)中的唯一值解譯出農(nóng)田、林地、草地3種主要的土地類型，可得到鄭州市域主要植被類型的分布圖 (圖1)。由圖1可知，鄭州市域內(nèi)土地利用以農(nóng)田為主，分布最廣。林地次之，與地形相結(jié)合，集中在西部、南部山區(qū)。草地較少，主要在北部黃河沿岸分布，同時在西部山區(qū)中草地與林地交錯分布。利用Google地球高分辨率影像數(shù)據(jù)對其進(jìn)行抽樣驗證，結(jié)果表明其能滿足鄭州市域尺度的分類精度。另外，采用2009年全球陸地覆蓋數(shù)據(jù)與2011年MODIS數(shù)據(jù)進(jìn)行綜合分析，由于兩者的獲取時間不一致，在采樣中同時借助Google地球高分辨率影像數(shù)據(jù)進(jìn)行數(shù)據(jù)比對調(diào)整，以避免隨著城市的發(fā)展，草地、農(nóng)田和林地的分布變化對研究結(jié)果造成影響。

圖1鄭州主要植被類型分布
Fig.1 Distribution of main plant types in Zhengzhou

3.2 濾波效果比較

頻率個數(shù)、曲線匹配閾值 (FET)、剩余點個數(shù) (DOD) 等參數(shù)是影響曲線平滑效果的重要因子，其設(shè)置沒有客觀標(biāo)準(zhǔn)。本研究經(jīng)過反復(fù)試驗，取頻率個數(shù)為2，F(xiàn)ET為1 000，DOD為8時，擬合效果最好。鄭州林地采樣點NDVI濾波前后的對比結(jié)果見圖2。由圖2可以看出，濾波前NDVI時間序列曲線呈無規(guī)則劇烈波動，難以進(jìn)行分析。濾波重構(gòu)的NDVI時間序列曲線排除了干擾值同時保持了曲線特征。經(jīng)HANTS濾波的NDVI時間序列曲線可直接用于分析植被物候期。

圖2HANTS濾波前后NDVI曲線對比
Fig.2 NDVI curve comparison before and after HANTS filtering

3.3 物候關(guān)鍵期的提取

基于研究對象特點及分布，采用樣點統(tǒng)計分析方法。將鄭州市域主要植被類型分布圖與處理后的研究區(qū)NDVI影像疊加，在草地、農(nóng)田分布區(qū)域內(nèi)分別隨機選取20個采樣點。草地采樣點植物主要包括：白頭翁 (Anemonechinensis)、東風(fēng)菜 (Doellingeriascaber)、卷丹 (Liliumlancifolium)、小糠草 (Agrostisgigantea)、馬唐 (Digitariasanguinalis)、胡枝子 (Lespedezabicolor)、復(fù)盆子 (Rubusidaeus) 等。考慮到林地在城市附近因人為的引種、選育、栽培各種林木而物候期差別甚大。為滿足本研究要求，故在以鄉(xiāng)土樹種為主的鄭州西部山區(qū)天然林地內(nèi)隨機選取非常綠植被采樣點20個。林地采樣點喬木主要包括：楓楊 (Pterocaryastenoptera)、刺槐 (Robiniapseudoacacia)、棗樹 (Zizyphusjujuba)、黃櫨 (Cotinuscoggygria)、櫟屬 (Quercus)、硬闊、軟闊等。逐一分析各類型植物的物候關(guān)鍵期。

經(jīng)過反復(fù)試驗，取林地、草地的滑動平均時間間隔為5個16 d合成周期 (圖3)。因小麥、玉米的生長季長度較短，農(nóng)田滑動平均時間間隔取3個16 d合成周期時效果最好 (圖4)。根據(jù)各植被類型的最適滑動平均時間，利用滑動平滑法對各植被的物候關(guān)鍵期進(jìn)行提取，提取結(jié)果見圖5。

圖3滑動平均前后林地NDVI曲線對比及物候參數(shù)提取
Fig.3 Comparison of NDVI curves and extraction of phenological parameters of forest land before and after the moving average

根據(jù)各植被物候期分布，可知鄭州西部山區(qū)林地大部分在第85～110天開始生長，在第215～240天生長達(dá)到最大程度，在第305～325天停止生長，林地生長季長度為200～220 d，生長時間較長 (圖5a)。大部分草地在第85～110天開始生長，在第220～240天生長旺盛，在第300～325天停止生長、枯萎，草地生長季長度為190～220 d (圖5b)。大部分冬小麥開始期在第40～60天，處于拔節(jié)期。在第100～120天生長達(dá)到極限，在第135～150天逐漸停止生長，進(jìn)入收獲期，小麥生長季長度集中在90～110 d (圖5c)。玉米在小麥?zhǔn)崭詈筮M(jìn)行播種，與小麥采用輪作模式。大部分玉米生長期集中在第165～180天，在第220～235天進(jìn)入抽穗期，在第260～270天生長結(jié)束，進(jìn)入收獲期，生長季長度最短，玉米生長季長度為90～105 d (圖5d)。

圖4滑動平均前后農(nóng)田NDVI曲線對比及物候參數(shù)提取
Fig.4 Comparison of cropland NDVI Curves and extraction of phenological parameters before and after the moving average

在鄭州市域內(nèi)不同地點，同類型植被的物候關(guān)鍵期略有波動，但總體保持一致。天然林地、草地形成了與鄭州自然氣候變化相適宜的物候期。農(nóng)作物小麥、玉米在人為作用下依自然季節(jié)的更替進(jìn)行有規(guī)律的輪作，在各氣候區(qū)均形成了相對穩(wěn)定的物候期。

3.4 物候關(guān)鍵期的驗證

驗證是基于遙感數(shù)據(jù)提取大面積物候信息的重要工作。柳晶等[15]的研究中，春季物候主要為前1～3個月的影響，平均溫度每升高1 ℃，物候期提前2.8～7.9 d、2.7～9.0 d；對于生長季，年平均溫度每升高1 ℃，物候期延長9.5～18.6 d，結(jié)束期推遲。

把劉玲玲等[3]對該區(qū)域的物候期研究作為真實值與本研究得出的物候期進(jìn)行對比分析，結(jié)果見表1。經(jīng)查詢鄭州2011年前3個月的均溫比2005年高1.17 ℃，對應(yīng)開始期提前12 d左右。2011年的年平均氣溫相較2005年升高0.2 ℃，對應(yīng)結(jié)束期推遲20 d左右，生長季長度延長35 d左右。可以發(fā)現(xiàn)，本研究中各項溫度均有所升高，物候期提前，結(jié)束期推遲，生長季長度延長，與上述研究結(jié)論基本是一致的。但由于環(huán)境條件改變及存在的可能誤差，物候期的延長天數(shù)相較柳晶的研究有所增加。另外，根據(jù)對玉米和小麥的實地調(diào)查發(fā)現(xiàn)，小麥返青期基本從2月開始,收割期為5月底。玉米播種基本上為6月初小麥?zhǔn)崭詈螅衩资斋@通常在9月底、10月初進(jìn)行。通過對比以上資料，看出基于MODIS NDVI數(shù)據(jù)獲得的植物物候期與研究和調(diào)查資料具有一定的可比性，表明可利用遙感來監(jiān)測植物關(guān)鍵物候期。

圖5各植被物候期分布
Fig.5 Distribution of vegetation phenology

表1 2005年與2011年物候期對比Table 1 Contrast of phenology between 2005 and 2011

4 結(jié)論與討論

本研究通過時間序列諧波分析法 (HANTS) 平滑重構(gòu)NDVI時間序列曲線，有效減弱了噪聲的干擾，能直接用于提取植物的物候信息。同時結(jié)合全球陸地覆蓋數(shù)據(jù)，解譯出鄭州市域主要植被類型分布圖，并針對林地、草地和農(nóng)田的不同特征，運用滑動平均法分別進(jìn)行各類型植物采樣點的滑動平均分析，實現(xiàn)快速、有效監(jiān)測鄭州市域植物關(guān)鍵物候期，包括開始期、成熟期、結(jié)束期以及生長季長度。通過與前人相關(guān)研究及實際調(diào)查的物候資料比較驗證，結(jié)果基本吻合，表明該方法獲得的物候期具有一定科學(xué)性。

對比于信芳等[1]與劉玲玲等[3]的研究，均表明通過HANTS濾波處理的MODIS NDVI時間序列數(shù)據(jù)可直接用于植物關(guān)鍵物候期的監(jiān)測。而與本研究采用滑動平均法不同，兩者均采用動態(tài)閾值法進(jìn)行物候關(guān)鍵期的提取，該方法通過設(shè)定閾值條件，將植物物候期限制在合理生育期內(nèi)，但閾值本身的確定易受主觀影響。針對本研究，滑動平均法對物候的計算則更為穩(wěn)定、可靠，但容易受春季雪融的影響，所得結(jié)果可能早于植物實際返青期。研究中已查詢鄭州地區(qū)2011年歷史天氣，發(fā)現(xiàn)不存在大量積雪狀況，可排除NDVI數(shù)據(jù)可能受春季雪融的影響。另外，兩者均通過逐個像元運算得到區(qū)域內(nèi)的關(guān)鍵物候期分布格局，相較本研究的隨機采樣統(tǒng)計分析，更加直觀的反應(yīng)了整個區(qū)域關(guān)鍵物候期的分布特點，可為本研究的深入提供借鑒。

本研究通過將NDVI時間序列數(shù)據(jù)與全球陸地覆蓋遙感數(shù)據(jù)相結(jié)合，實現(xiàn)了不同植被類型物候的分類研究、比較。但受遙感影像時間分辨率、空間分辨率的限制，僅研究了植物的3個物候關(guān)鍵期，下一步可利用高精度影像詳細(xì)研究植物各階段的物候特征。另外，隨機采樣可能會對提取的植被物候精確度產(chǎn)生影響，今后也可通過提取多年的物候平均值來做相關(guān)性分析將會更加科學(xué)。

[參 考 文 獻(xiàn)]

[1] 于信芳, 莊大方. 基于MODIS NDVI數(shù)據(jù)的東北森林物候期監(jiān)測[J]. 資源科學(xué), 2006, 28(4): 111-117.

[2] Wang Q, Tenhunen J D. Vegetation mapping with multitemporal NDVI in North Eastern China Transect (NECT) [J]. International Journal of Applied Earth Observation and Geoinformation, 2004, 6(1): 17-31.

[3] 劉玲玲, 劉良云, 胡勇. 基于AVHRR和MODIS數(shù)據(jù)的全球植被物候比較分析[J]. 遙感技術(shù)與應(yīng)用, 2012, 27(5): 754-762.

[4] 強建華, 趙鵬祥, 陳國領(lǐng). 基于NDVI的油松天然林生長狀況的遙感監(jiān)測研究[J]. 西北林學(xué)院學(xué)報, 2007, 22(1): 149-151, 172.

[5] Zhang X Y, Friedl M A, Schaaf C B, et al. Monitoring vegetation phenology using MODIS[J]. Remote Sensing of Environment, 2003, 84(3): 471-475.

[6] 朱源, 彭光雄, 王志, 等. 西藏林芝地區(qū)近30 a來的NDVI變化趨勢研究[J]. 西北林學(xué)院學(xué)報, 2011, 26(4): 69-74.

[7] 張峰, 吳炳方, 劉成林, 等. 利用時序植被指數(shù)監(jiān)測作物物候的方法研究[J]. 農(nóng)業(yè)工程學(xué)報, 2004, 20(1): 155-159.

[8] 侯學(xué)會, 牛錚, 高帥, 等. 基于SPOT-VGT NDVI時間序列的農(nóng)牧交錯帶植被物候監(jiān)測[J]. 農(nóng)業(yè)工程學(xué)報, 2013, 29(1): 142-150, 294.

[9] Julien Y, Sobrino J A. Comparison of cloud-reconstruction methods for time series of composite NDVI data[J]. Remote Sensing of Environment,2010,114(3): 618-625.

[10] Menenti M, Azzali S, Verhoef W, et al. Mapping agroecological zones and time lag in vegetation growth by means of fourier analysis of time series of NDVI images[J]. Advances in Space Research, 1993, 13(5): 233-237.

[11] Roerink G J, Su Z, Menenti M. S-SEBI: a simple remote sensing algorithm to estimate the surface energy balance[J]. Physics and Chemistry of the Earth, Part B: Hydrology, Oceans and Atmosphere, 2000, 25(2): 147-157.

[12] 李正國, 楊鵬, 周清波, 等. 基于時序植被指數(shù)的華北地區(qū)作物物候期/種植制度的時空格局特征[J]. 生態(tài)學(xué)報, 2009, 29(11): 6216-6226.

[13] Julien Y, Sobrino J A, Verhoef W. Changes in land surface temperatures and NDVI values over Europe between 1982 and 1999[J]. Remote Sensing of Environment, 2006, 103(1): 43-55.

[14] 崔凱. 基于遙感技術(shù)的作物物候監(jiān)測方法及動態(tài)變化分析研究[D]. 長沙: 中南大學(xué), 2012.

[15] 柳晶, 鄭有飛, 趙國強, 等. 鄭州植物物候?qū)夂蜃兓捻憫?yīng)[J]. 生態(tài)學(xué)報, 2007, 27(4): 1471-1479.