彭檢貴，宋亞斌，寧小斌，江騰宇

(國家林業(yè)和草原局中南調(diào)查規(guī)劃設(shè)計院, 長沙 410014)

0 引言

植被的物候變化直接反映了周邊地區(qū)的氣候變化，成為植被動態(tài)模擬的關(guān)鍵[1-2]。同時，植被物候數(shù)據(jù)可以為植被動態(tài)監(jiān)測、區(qū)域農(nóng)業(yè)規(guī)劃、區(qū)域資源配置與決策、區(qū)域氣候變化研究等提供重要參考。

植物物候監(jiān)測一般可分為傳統(tǒng)監(jiān)測方法和遙感監(jiān)測方法。傳統(tǒng)物候監(jiān)測方法需要耗費大量的人力、物力，而且受時間限制，監(jiān)測周期短，覆蓋區(qū)域小[3-4]。而遙感物候監(jiān)測方法僅需要耗費相對較低的費用而獲取大尺度區(qū)域的物候參數(shù)，能及時、有效地對植物物候進(jìn)行監(jiān)測[4-8]。已有大量學(xué)者基于NDVI和EVI植被指數(shù)的各種遙感數(shù)據(jù)源的時間序列提出了大量物候特征參數(shù)遙感提取方法，如White等[3]將傳統(tǒng)物候觀測模型與遙感物候參數(shù)提取結(jié)合，針對特定生物群落生態(tài)系統(tǒng)，提出NDVI比率閾值物候參數(shù)提取方法；Reed等[4]基于NDVI時間序列數(shù)據(jù)，使用滑動平均方法計算提取了美國4種代表性植物的(荒漠灌叢、森林、草地、農(nóng)作物)物候特征；Yu等[9]提出了最大變化斜率方法，通過限制不同閾值的最大變化斜率范圍，計算了中亞東部區(qū)域植物生長季開始期；Zhang等[10]通過對每年NDVI時間序列進(jìn)行分段Logistic函數(shù)擬合，利用擬合后的曲線的曲率變化的特點，得出植物的生長季開始期、結(jié)束期、生長季長度以及生長季峰值等。已有的研究表明，基于遙感數(shù)據(jù)的物候特征提取方法主要有經(jīng)驗公式法、NDVI比率閾值法、滑動平均法、閾值法、最大變化斜率法及Logistic函數(shù)擬合法等。不同的物候特征提取方法適應(yīng)的地區(qū)也不一樣，因此需要分析與對比。

洞庭湖流域在長江流域水系中占十分重要的比重，也是我國商品糧的主要生產(chǎn)基地，棉麻等經(jīng)濟(jì)作物也非常豐富。為此，本文以洞庭湖流域為研究區(qū)域，基于時間序列EVI數(shù)據(jù)，采用滑動平均法、動態(tài)閾值法、最大變化斜率法以及Logistic函數(shù)擬合法4種方法進(jìn)行物候參數(shù)提取，并對四種物候參數(shù)提取方法進(jìn)行對比分析，為洞庭湖流域植被物候提供技術(shù)方法支持。

1 研究區(qū)域與數(shù)據(jù)

洞庭湖流域即洞庭湖水系流經(jīng)的廣大地區(qū)，屬于典型的亞熱帶季風(fēng)氣候，年均氣溫 16.5～17.0 ℃。流域多年平均年降水量1427 mm，多年平均年徑流量為2016億m3，約占長江流域地表水資源的21%，其比重為長江流域各水系之首。洞庭湖流域覆蓋湖南省大部和湖北省、廣西壯族自治區(qū)、貴州省和重慶市部分地區(qū)，流域面積約262800 km2，占長江流域面積的14%。洞庭湖平原區(qū)為湖南農(nóng)業(yè)主產(chǎn)區(qū)，以種植糧食、棉花為主，也是中國主要淡水養(yǎng)殖區(qū)之一，耕地面積約占湖南全省的1/6，居住人口約占湖南全省1/9，水、氣、生物等資源豐富。

本文選用已進(jìn)行畸變校正和定標(biāo)之后的MODIS 13Q1產(chǎn)品的EVI數(shù)據(jù)與MODIS 12Q1遙感數(shù)據(jù)經(jīng)過數(shù)據(jù)提取、影像鑲嵌、投影轉(zhuǎn)換、裁剪等一系列處理后進(jìn)行研究。其中，用TIMESAT軟件對MODIS 13Q1產(chǎn)品中的EVI數(shù)據(jù)進(jìn)行了S-G濾波重構(gòu)處理；MODIS 12Q1為土地覆蓋數(shù)據(jù)，主要包括林地、草地、耕地、水域、城市用地等類型。

2 研究方法

2.1 滑動平均法

滑動平均法主要是將植被指數(shù)時間序列數(shù)據(jù)突然增長時對應(yīng)的時間代表植物開始生長的時期，是在回歸滑動平均模型的基礎(chǔ)上，將植被指數(shù)時間序列曲線與計算后的滑動平均曲線進(jìn)行比較，以反映實際植被指數(shù)時間序列曲線和既定趨勢之間的偏差。滑動平均模型的公式為：

Yt=(Xt+Xt-1+Xt-2+…+Xt-(w-1))/w

(1)

其中，Yt是時間t時滑動平均模型的計算值，Xt表示時間t時的標(biāo)準(zhǔn)植被指數(shù)時間序列值，表示滑動平均模型的窗口大小，即滑動平均時間間隔[11-13]。

2.2 動態(tài)閾值法

動態(tài)閾值法是在固定閾值法的基礎(chǔ)上，由Jonsson等提出，相較于固定閾值法其考慮了植被指數(shù)時間序列值的變化幅度[14-17]，在一定程度上消除了植被類型和土壤背景值的影響。具體研究時，根據(jù)植被指數(shù)時間序列曲線的變化特點，將曲線最初階段的植被指數(shù)值快速增加處定義為植被生長開始期，可具體確定為植被指數(shù)時間序列重構(gòu)曲線的左邊起始值或最小值增長到指定程度(如整體增幅的10%)，其他物候特征說明見表1。

2.3 最大變化斜率法

最大變化斜率法將植物在春季開始迅速生長的時期定義為植被生長季的開始。具體的方法步驟為：第一步，建立標(biāo)準(zhǔn)的植被指數(shù)時間序列數(shù)據(jù)；第二步，利用最大變化斜率模型計算植物的生長季開始點，并設(shè)置條件閾值剔除掉結(jié)果中不合理的計算值。最大變化斜率法的公式為：

Δyt=yt-yt-1，θt=arctg(Δyt)

Δyt+1=yt+1-yt，θt+1=arctg(Δyt+1)

Δθ(t+1)=θt+1-θt

(2)

表1 動態(tài)閾值法植被物候特征定義物候參數(shù)參數(shù)含義生長季開始時間EVI增加至擬合函數(shù)的左半部分EVI振幅的20%的時刻生長季結(jié)束時間EVI降低至擬合函數(shù)的右半部分EVI振幅的20%的時刻生長季長度開始生長時間與停止生長時間之間的時間間隔基準(zhǔn)值擬合函數(shù)左半部分和右半部分EVI最小值的平均值生長季中期時刻EVI增至擬合函數(shù)左半部分EVI 振幅的80%的時刻和EVI降至擬合函數(shù)右半部分EVI振幅的80%的時刻的平均值生長季EVI最大值EVI的最大值生長季EVI振幅EVI 峰值與左半部分最小值和右半部分最小值的均值之間的差值生長速度擬合曲線左邊振幅20%和80%之間的EVI差值與相應(yīng)的時間差值的比值減緩速度擬合曲線右邊振幅20%和80%之間的EVI差值與相應(yīng)的時間差值的比值生長季EVI活躍累積量EVI擬合曲線與基準(zhǔn)值以上的區(qū)域圍成的面積生長季EVI總累積量EVI擬合曲線與作物開始生長到生長結(jié)束時間段內(nèi)的分量

其中，yt-1指t-1時的植被指數(shù)值，yt是指t時的植被指數(shù)值，yt+1是指t+1時植被指數(shù)值，Δθt和Δθt+1指t和t+1時的角度變化值[18-22]。

2.4 Logistic函數(shù)擬合法

該方法通過曲率變化率的極值點反映植物各物候期的轉(zhuǎn)換[23-26]。具體的方法步驟為：

第一步，建立標(biāo)準(zhǔn)的植被指數(shù)時間序列數(shù)據(jù)，并根據(jù)實際情況劃分出植被指數(shù)值持續(xù)增長或下降的階段。

第二步，利用Logistic函數(shù)對劃分出來的各個階段進(jìn)行擬合，生成擬合曲線，擬合函數(shù)為：

(3)

其中，t為具體的時間值，y(t)為隨時間t變化的植被指數(shù)擬合值，d為標(biāo)準(zhǔn)的植被指數(shù)值；a，b，c為擬合函數(shù)的參數(shù)。

第三步，計算Logistic擬合后的植被指數(shù)時間序列曲線的曲率K及變化率，將曲率變化率的極大值定義為植物生長的開始時間。曲率K及變化率K′的計算公式如下：

(4)

其中，z=e(a+bt)，α為曲線在時間t時的單位弧長轉(zhuǎn)過的角度，s為單位弧長。

(5)

本文利用MATLAB軟件，通過編程，實現(xiàn)上述4種方法，并利用4種方法對洞庭湖流域植被的物候信息進(jìn)行了提取與對比分析。

3 結(jié)果與分析

3.1 不同方法下的物候特征計算結(jié)果

利用4種方法提取的2005與2015年洞庭湖流域的部分物候特征如表2，包括生長季開始(start of season， SOS)、生長季結(jié)束(end of season， EOS)以及生長季長度(length of season， LOS)，EVI最大值(max of EVI， MOE)。

表2 不同方法下的物候特征參數(shù)d提取方法SOSEOSLOSMOE2005年2015年2005年2015年2005年2015年2005年2015年動態(tài)閾值法24753213532982790.830.91最大變化斜率法33332252731922400.830.91Logistic函數(shù)擬合法1849324916366700.830.91滑動平均法51423262372761960.830.91 注:像元坐標(biāo)為x=265,y=415

依據(jù)湖南氣象網(wǎng)農(nóng)業(yè)氣象旬、月報上的相關(guān)實地觀測資料，洞庭湖流域主要農(nóng)作物為油菜、水稻、棉花和玉米，其物候分別為：油菜一般為上一年9月底10月初種植，翌年5月中下旬左右成熟收割；水稻當(dāng)年4月上旬左右為返青期，雙季稻7月左右收割第一季，然后種第二季，11月份收割，而單季稻于當(dāng)年10月份左右收割；棉花當(dāng)年4月左右為播種出苗期，9月左右開始收獲，直到11月份；春玉米當(dāng)年四月上旬左右為播種出苗期，6月左右為成熟期?？紤]到研究區(qū)域內(nèi)有一年一熟、一年兩熟的作物，這次實驗提取的為年內(nèi)植被生長的開始期和結(jié)束期。由于實驗所采用的數(shù)據(jù)分辨率為250 m，一個像元內(nèi)的植被可能并不是單一的，綜合對比實驗提取的物候特征參數(shù)和實地觀測資料，可以看出：①動態(tài)閾值法提取2005年像元的SOS值為第24天，早于4種農(nóng)作物觀測數(shù)據(jù)的生長開始時間；提取2015年像元的EOS值為第353天，相比其他三種方法的提取結(jié)果過晚，也晚于4種農(nóng)作物觀測數(shù)據(jù)生長結(jié)束時間。②最大變化斜率法提取的2005年和2015年的SOS普遍過早，且由于最大變化斜率法的模型問題，影像一個時間間隔內(nèi)的斜率變化角相同，使得最終的計算結(jié)果存在誤差。③Logistic函數(shù)擬合法，因研究區(qū)域內(nèi)多種土地覆蓋類型，針對個別土地覆蓋類型如水體等，其植被指數(shù)時間序列曲線呈鋸齒狀重復(fù)出現(xiàn)，若直接采用傳統(tǒng)的Logistic模型方法對原始的植被指數(shù)時間序列影像進(jìn)行分段擬合，MATLAB在求解Logistic函數(shù)時會出現(xiàn)問題，導(dǎo)致出錯，因此利用MOE前后的植被指數(shù)時間序列數(shù)據(jù)，分段用植被指數(shù)時間序列累積曲線來進(jìn)行Logistic函數(shù)擬合。Logistic函數(shù)擬合法提取2005年像元的SOS值為第184天，相比其他三種方法的提取結(jié)果過晚，也晚于4種農(nóng)作物觀測數(shù)據(jù)的生長開始時間；提取2015年像元的EOS值為第163天，相比其他三種方法的提取結(jié)果過早，且也早于4種農(nóng)作物的綜合生長結(jié)束時間。④滑動平均法提取2005年和2015年像元的SOS值分別為第51天和第42天，EOS值分別為第326天和237天，對比其他三種方法，該方法的提取結(jié)果與4種農(nóng)作物觀測數(shù)據(jù)生長開始時間和結(jié)束時間接近，準(zhǔn)確性較高。綜合考慮到上述4種作物物候特征參數(shù)提取方法的模型原理和實驗數(shù)據(jù)物候參數(shù)提取結(jié)果，滑動平均法模型簡單，效果較好，更適用于整個研究區(qū)域的物候特征參數(shù)提取。

3.2 近10年洞庭湖流域植被物候特征參數(shù)變化分析

選擇滑動平均法對2005年與2015年整個洞庭湖流域EVI植被指數(shù)時間序列影像進(jìn)行物候特征參數(shù)提取，提取結(jié)果如圖1。

從圖1可看出，2005年SOS值從0～第120天均有一定區(qū)域的分布，但2015年SOS值大部分分布在0～第80天之間，集中在第25～75天之間；且2005年圖中中部區(qū)域的SOS的值在第80～120天之間的，在2015年時主要分布在第40～80天之間，提前了40 d左右；2005年研究區(qū)東北部和南部區(qū)域的SOS值在第40～50天之間的，在2015年時主要分布在第60～75天之間，推遲了20 d左右；2005年和2015年，洞庭湖流域EOS值整體上都分布在第220～350天之間，但2005年研究區(qū)北部和中東部區(qū)域的EOS值在第300～350天之間的，在2015年時主要分布在第230～320天之間，提前了30～70 d左右；2005年研究區(qū)中西部區(qū)域的EOS值在第220～260天之間的，在2015年時主要分布在第300～350天之間，推遲了85 d左右； 2005年北部和東部區(qū)域的LOS值主要集中在第270～320天之間，在2015年主要分布在第200～270天之間，減短了60 d左右；中西部區(qū)域2005年LOS值集中在第200～260天之間的，在2015年分布在第240～300天之間，延長了50 d左右。

圖1 2005年和2015年物候特征參數(shù)

續(xù)圖1 2005年和2015年物候特征參數(shù)

利用MODIS 12Q1數(shù)據(jù)與洞庭湖流域物候特征參數(shù)數(shù)據(jù)疊加，并進(jìn)行統(tǒng)計分析，得到了洞庭湖流域的林地、草地和耕地這三種主要植被覆蓋類型的物候特征參數(shù)，見表3。

表3 不同植被覆蓋類型的物候特征參數(shù)d植被覆蓋類型SOS均值EOS均值LOS均值2005年2015年2005年2015年2005年2015年林地51.46652.294305.23300.362254.213248.538草地52.16853.205304.646297.831252.93245.098耕地52.85652.6311.253297.647258.847245.522

從表3可以看出，2005年林地、草地、耕地的SOS均值分別為51.466，52.168，52.856；2015年林地、草地、耕地的SOS均值分別為52.294，53.205，52.6；兩年不同植被覆蓋類型的SOS均值相差不遠(yuǎn)。2005年林地、草地、耕地的EOS均值分別為305.23，304.646，311.253；2015年林地、草地、耕地的EOS均值分別為300.362，297.831，297.647，所以2015年林地和草地的EOS均值比2005年提前了6 d左右，2015年耕地的EOS均值比2005年提前了14 d左右。2005年林地、草地、耕地的LOS均值分別為254.213，252.93，258.847；2015年林地、草地、耕地的LOS均值分別為248.538，245.098，245.522。分區(qū)統(tǒng)計結(jié)果顯示，兩年SOS均值變化不大，2015年EOS均值比2005年提前了6～14 d左右，而2015年LOS值比2005年減短了6～13 d左右。

這些變化引起的原因有氣候變化、農(nóng)作物類型發(fā)生變化、耕種方式發(fā)生變化等，但從湖南氣象網(wǎng)了解到，農(nóng)作物類型和耕種方式變化不大，氣候變化是主要原因。2005年年均氣溫18 ℃、年降水量約1 300 mm、年日照時數(shù)1 400 h左右，2015年年均氣溫17.9 ℃、年降水量約1 580.5 mm、年日照時數(shù)1 180 h左右，降水量較常年偏多。且2015年相比2005年，年初氣溫異常偏高，出現(xiàn)強(qiáng)“暖冬”現(xiàn)象，清明節(jié)前后，卻出現(xiàn)“倒春寒”現(xiàn)象，盛夏時氣溫異常偏低，立冬后又出現(xiàn)罕見冬汛，所以這些氣候變化是造成物候特征參數(shù)值空間格局發(fā)生變化的主要原因，且氣候變化對對農(nóng)作物的影響較大，對林地、草地影響較小。

4 結(jié)論

以洞庭湖流域為研究區(qū)域，采用MATLAB GUI將四種物候特征參數(shù)提取方法通過編程實現(xiàn)，然后分析了四種物候特征參數(shù)提取方法的優(yōu)劣性，最后分析了2005年和2015年物候參數(shù)的空間變化及原因。結(jié)論如下：

1)采用滑動平均法計算各類植被生長季開始時間簡單、有效，并且克服了閾值法在不同土壤背景、光照以及植被類型等條件下需要設(shè)置不同閾值的問題，但該方法如果受氣候事件如雪融影響，所得結(jié)果可能會早于植物實際生長季開始時間。

2)動態(tài)閾值法計算生長季開始和結(jié)束時間簡單、耗時短，但需要多次實驗才能確定適合區(qū)域植被覆蓋類型的最佳閾值。

3)最大變化斜率法充分考慮植被指數(shù)時間序列曲線特征，將生長季開始時間評估限制在合理生長期內(nèi)，以減少對稀疏植被地區(qū)不經(jīng)歷主要生長季開始時間事件的錯誤計算，且設(shè)定閾值條件可提高植物生長季開始時間的計算效率與精度，但相關(guān)閾值的設(shè)置須主觀綜合考慮研究區(qū)域、植被類型及背景條件的具體實際情況。閾值條件選擇合理與否，直接影響計算結(jié)果及運行效率。

4)Logistic函數(shù)擬合法監(jiān)測大區(qū)域植被動態(tài)變化靈活方便，結(jié)果理想；逐個處理像素時不需設(shè)置閾值或經(jīng)驗限制條件，因此具有適用性；由于是分段式擬合，能對一年內(nèi)農(nóng)田和半干旱地區(qū)常見的多熟物候行為進(jìn)行監(jiān)測；不受植被物候歷法時期限制，為監(jiān)測接近真實時間的植被物候提供了可能，但實際的植被指數(shù)時間序列曲線并非理想的規(guī)則曲線，突變點較多，擬合精度的高低直接影響植被物候參數(shù)的確定。

5)2005年與2015年洞庭湖流域物候特征參數(shù)的空間格局的變化主要是由氣候變化引起，且氣候變化對林地和草地的影響較小，對耕地的影響較大。

總體而言，從多時相遙感影像中獲取的植被指數(shù)時間序列數(shù)據(jù)，經(jīng)時間序列重構(gòu)后可以較為準(zhǔn)確地提取植物物物候特征參數(shù)，反映植物的生長過程，并可以反映出區(qū)域氣候變化，適合大范圍、快速監(jiān)測植物的物候變化。