衛(wèi) 煒，吳文斌，李正國，楊鵬，胡瓊，周清波※

(1．中國農(nóng)業(yè)科學(xué)院農(nóng)業(yè)資源與農(nóng)業(yè)區(qū)劃研究所，北京 100081;2．農(nóng)業(yè)部農(nóng)業(yè)信息技術(shù)重點實驗室，北京 100081)

植被指數(shù)是多光譜遙感數(shù)據(jù)經(jīng)過空間變換或者不同波段之間的線性、非線性組合構(gòu)成，其基本原理是利用植被的光譜特征，通過強化可見光 (主要是紅光)與近紅外波段反射率之間的差異來反映植被的生長狀況［1］。時間序列植被指數(shù)能夠較好反映植被的季節(jié)性生長發(fā)育過程，是監(jiān)測不同尺度植被覆蓋及其動態(tài)變化特征，提取植被生長季開始、生長最盛和結(jié)束等物候信息，描述全球氣候變化及陸地生態(tài)系統(tǒng)響應(yīng)的重要工具［2］。然而，星載傳感器在數(shù)據(jù)采集和處理過程中受到各種因素 (如太陽高度角、傳感器觀測角、云、水汽和氣溶膠等)干擾，使得植被指數(shù)時序曲線波動大并出現(xiàn)許多噪聲［3］。因此，在進(jìn)行各種趨勢分析和信息提取之前，有必要先對時間序列植被指數(shù)進(jìn)行去噪和平滑處理，即時序植被指數(shù)重構(gòu)［4］。

盡管目前大多數(shù)植被指數(shù)產(chǎn)品都采用了最大值合成法來部分消除一些干擾因素的影響［5］，但是得到的時序植被指數(shù)仍然存在大量的噪聲。為此，研究者們提出了許多時序植被指數(shù)重構(gòu)方法，總體上可以歸為濾波方法和函數(shù)擬合方法兩大類。濾波方法分為時域濾波和頻域濾波，時域濾波是在一個給定大小的濾波窗口內(nèi)，采用某種方法對數(shù)據(jù)中的噪聲進(jìn)行處理，濾波窗口的大小與最終的濾波效果有關(guān)，是一種局部處理的方法;頻域濾波則是通過數(shù)學(xué)變換將數(shù)據(jù)從時間域轉(zhuǎn)換到頻率域，噪聲通常存在于高頻部分，利用一個低通濾波器即可實現(xiàn)去噪，是一種整體處理的方法。函數(shù)擬合方法則是通過某種形式的函數(shù)對時序數(shù)據(jù)曲線進(jìn)行最小二乘擬合，用擬合后的平滑曲線代替原來的時序數(shù)據(jù)曲線來實現(xiàn)平滑去噪。同樣，函數(shù)擬合方法也存在局部處理和整體處理之分，其中經(jīng)過局部擬合后的曲線還需要一個整體函數(shù)將它們連接起來以構(gòu)成最終的擬合曲線。常用的濾波方法包括Savitzky-Golay濾波法、最佳指數(shù)斜率提取法、滑動中值/均值濾波法和基于傅里葉/小波變換的頻域低通濾波法等;而常用的函數(shù)擬合方法包括非對稱高斯函數(shù)擬合法、雙Logistic函數(shù)擬合法和傅里葉函數(shù)擬合法等。這些方法有其各自的特點和區(qū)域適宜性，因此根據(jù)應(yīng)用目的和研究區(qū)域?qū)Σ煌椒ㄟM(jìn)行對比和評價具有十分重要的意義。國內(nèi)外一些學(xué)者已經(jīng)對部分重構(gòu)方法進(jìn)行了比較研究［6-9］，然而上述研究多是針對森林、草地等自然植被植被指數(shù)進(jìn)行，且僅針對不同方法的重構(gòu)效果進(jìn)行了對比，很少考慮原始數(shù)據(jù)質(zhì)量的影響和重構(gòu)結(jié)果的應(yīng)用情況。

耕地資源是保障國家糧食安全的基礎(chǔ)性資源，耕地植被指數(shù)是進(jìn)行農(nóng)作物長勢監(jiān)測與估產(chǎn)以及農(nóng)業(yè)生態(tài)系統(tǒng)對氣候變化響應(yīng)等研究的重要工具［10］。東北三省是我國重要的商品糧生產(chǎn)基地，在我國糧食安全保障體系和農(nóng)業(yè)生產(chǎn)中占有重要地位［11］，同時該地區(qū)也是受氣候變化影響最為顯著的區(qū)域之一［12］。該研究基于耕地植被指數(shù)，在東北三省利用MODIS數(shù)據(jù)產(chǎn)品對非對稱高斯函數(shù)擬合、雙Logistic函數(shù)擬合和Savitzky-Golay濾波3種常見的時序植被指數(shù)重構(gòu)方法進(jìn)行比對研究。從原始數(shù)據(jù)質(zhì)量對時間序列植被指數(shù)重構(gòu)的影響、不同重構(gòu)方法對原始數(shù)據(jù)的保真性和基于不同重構(gòu)方法的耕地生長季長度空間格局等方面進(jìn)行分析，以此來評價不同方法在該地區(qū)耕地時序植被指數(shù)重構(gòu)中的適宜性。

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

1.1 研究區(qū)域

東北三省包括黑龍江省、吉林省和遼寧省，區(qū)域總面積79．18萬km2，其中耕地面積為26．44萬km2，約占全國耕地總面積的16．5%［13］(圖1)。該地區(qū)屬寒溫帶、中溫帶濕潤 (半濕潤)氣候，大部分地區(qū)≥0℃積溫為2 000～4 200℃，≥10℃積溫為1 600～3 600℃，夏季平均氣溫為20～25℃，無霜期為140～170天，年降雨量500～800mm，其中60%集中在7～9月。受熱量條件限制，該地區(qū)作物種植制度為1年1熟制，主要作物類型包括大豆、春玉米、春小麥和一季稻［14］。

圖1 研究區(qū)域耕地范圍及分布

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

該研究采用的遙感數(shù)據(jù)為搭載在EOS/Terra衛(wèi)星上的MODIS傳感器MOD13A2數(shù)據(jù)產(chǎn)品，版本為v005。涉及到的數(shù)據(jù)內(nèi)容包括基于最大值合成 (MVC)方法合成的16天1km空間分辨率的NDVI及其質(zhì)量控制數(shù)據(jù)，該數(shù)據(jù)產(chǎn)品已經(jīng)完成了幾何校正和大氣校正等處理。數(shù)據(jù)從美國NASA LP DAAC(Land Processes Distributed Active Archive Center)工作組官方網(wǎng)站下載 (https://lpdaac．usgs．gov/products/modis_products_table)，時間范圍為2010年1月～2012年12月，3年時間共69期影像。數(shù)據(jù)產(chǎn)品中NDVI的DN值范圍為-3 000～10 000，其中-3 000為填充值，將DN值轉(zhuǎn)換成標(biāo)準(zhǔn)植被指數(shù)值的公式為:VI=DN*0．000 1。該研究使用的質(zhì)量控制數(shù)據(jù)為MOD13A2中的Data Pixel Reliability數(shù)據(jù)集，其在時間和空間分辨率上與NDVI數(shù)據(jù)集相匹配。

該研究使用的土地利用數(shù)據(jù)來自中國1:25萬土地覆蓋遙感調(diào)查與監(jiān)測數(shù)據(jù)庫的2005年土地覆蓋數(shù)據(jù)，該數(shù)據(jù)是在全國1:10萬土地利用數(shù)據(jù)和遙感分類數(shù)據(jù)的基礎(chǔ)上，并結(jié)合輔助資料、野外考察記錄、遙感圖像分析等通過全數(shù)字制圖完成的。數(shù)據(jù)內(nèi)容包括森林、草地、農(nóng)田、聚落、濕地與水體、荒漠等6個一級類型，空間分辨率為100m。該研究將東北三省的農(nóng)田數(shù)據(jù)提取出來作為耕地圖層，并將空間分辨率重采樣到1km使之與MOD13A2的NDVI數(shù)據(jù)相匹配，然后基于該圖層提取東北三省耕地的植被指數(shù)進(jìn)行分析。

2 研究方法

2.1 時序植被指數(shù)重構(gòu)方法

星載傳感器在獲取地表信息時會受到諸多因素的影響，即使在經(jīng)過了一定處理的遙感數(shù)據(jù)產(chǎn)品中仍然殘留了很多噪聲，導(dǎo)致時間序列植被指數(shù)曲線的周期性變化趨勢不明顯［15］。因此，對時序植被指數(shù)進(jìn)行去噪和平滑等預(yù)處理，重構(gòu)其時間序列曲線至關(guān)重要［16］。該研究主要對3種常見的時序植被指數(shù)重構(gòu)方法:非對稱高斯函數(shù)擬合法、雙Logistic函數(shù)擬合法和Savitzky-Golay濾波法進(jìn)行對比和分析。

2．1．1 非對稱高斯函數(shù)擬合法

非對稱高斯函數(shù) (Asymmetric Gaussian，AG)擬合法是一種由局部擬合到整體擬合的方法，使用高斯函數(shù)來分段模擬植被的生長過程，最后通過平滑連接各段高斯擬合曲線實現(xiàn)時間序列的重構(gòu)［17］。其主要過程大致可以分為區(qū)間提取、局部擬合和整體連接3個步驟。首先提取原始時序數(shù)據(jù)曲線中的谷值和峰值，采用高斯函數(shù)分別擬合曲線的左右部分。針對曲線突出部分?jǐn)M合效果欠佳的問題，將曲線劃分成左邊谷值區(qū)、中部峰值區(qū)與右邊谷值區(qū)，分別用不同的局部擬合函數(shù)進(jìn)行描述，最后，再利用各局部擬合函數(shù)構(gòu)建整體擬合函數(shù)。局部擬合函數(shù)為:

其中g(shù)(t;a1，…，a5)為高斯函數(shù)，式中c1和c2控制曲線的基準(zhǔn)和幅度，a1決定峰值和谷值的位置;a4、a5和a2、a3分別控制曲線左、右部分的寬度和陡峭度。整體擬合函數(shù)為:

其中 ［tL，tR］是時序數(shù)據(jù)中待擬合部分的變化區(qū)間，fL(t)、fC(t)和fR(t)分別代表 ［tL，tR］區(qū)間內(nèi)左邊谷值、中間峰值及右邊谷值對應(yīng)的局部擬合函數(shù)，α(t)和β(t)為介于0和1之間的剪切系數(shù)。通過整體擬合函數(shù)將局部擬合函數(shù)連接起來是該方法的關(guān)鍵之一，這種從局部到整體的擬合策略避免了整體數(shù)據(jù)對局部擬合的干擾，擬合后的曲線更加接近真實情況。

2．1．2 雙Logistic函數(shù)擬合法

雙Logistic函數(shù) (Double Logistic，DL)擬合法同樣是基于從局部到整體的擬合思想，處理過程與AG擬合方法相類似，區(qū)別之處在于局部擬合函數(shù)為雙Logistic形式且比AG函數(shù)少一個參數(shù)［18］。DL擬合的局部擬合函數(shù)為:

其中，g(t;a1，…，a4)為雙Logistic函數(shù)，式中c1和c2控制曲線的基準(zhǔn)和幅度;a1、a2和a3、a4分別控制曲線左、右部分的拐點位置及拐點處的變化速率，整體擬合函數(shù)與AG擬合相同。

2．1．3 Savitzky-Golay 濾波法

Savitzky-Golay(S-G)濾波法是Savizky和Golay于1964年提出的，它是一種通過局部多項式回歸模型來平滑時序數(shù)據(jù)的時域低通濾波方法［19］。S-G濾波的基本思想是基于多項式，在濾波窗口內(nèi)利用最小二乘法對數(shù)據(jù)進(jìn)行最佳擬合，可以簡單的理解為一種加權(quán)平均算法，權(quán)重由在一個濾波窗口范圍內(nèi)做最小二乘擬合的多項式次數(shù)來決定［20］。S-G濾波公式為:

其中，Y為原始數(shù)據(jù)，Y'為擬合值，Ci為第i個點的權(quán)重，N=2n+1為濾波窗口的大小。S-G濾波法對濾波窗口的大小非常敏感，濾波窗口的寬度設(shè)置偏小容易產(chǎn)生大量冗余數(shù)據(jù);反之則可能遺漏一些細(xì)節(jié)信息［21］。另外，擬合多項式的次數(shù)也會對平滑效果產(chǎn)生影響，次數(shù)較低時結(jié)果較為平滑，當(dāng)次數(shù)較高時則會導(dǎo)致過度擬合。

2.2 不同重構(gòu)方法對比策略

該研究主要從3個方面對AG擬合、DL擬合和SG濾波3種時間序列植被指數(shù)重構(gòu)方法進(jìn)行比較。首先對NDVI及其質(zhì)量控制數(shù)據(jù)進(jìn)行分析，探討原始數(shù)據(jù)質(zhì)量對時序植被指數(shù)重構(gòu)的影響;然后選擇若干隨機抽樣點，從定性和定量兩個角度來評價不同方法對時序植被指數(shù)的重構(gòu)效果;最后借助研究區(qū)的耕地物候特征，間接對不同的時序植被指數(shù)重構(gòu)方法進(jìn)行對比說明。

2．2．1 原始數(shù)據(jù)質(zhì)量的影響分析

根據(jù)Data Pixel Reliability數(shù)據(jù)集的DN值可將所有像元的數(shù)據(jù)質(zhì)量劃分為5個等級 (表1)［22］。由表1可知，MOD13A2數(shù)據(jù)產(chǎn)品中植被指數(shù)的數(shù)據(jù)質(zhì)量主要受到冰雪覆蓋和云掩的影響。在時序植被指數(shù)進(jìn)行重構(gòu)過程中，這些數(shù)據(jù)質(zhì)量受到影響的像元點可能會使真實的植被指數(shù)曲線發(fā)生扭曲或偏移。本研究以Data Pixel Reliability數(shù)據(jù)集的DN值為依據(jù)對各個數(shù)據(jù)點賦予不同的權(quán)重，使得重構(gòu)后的植被指數(shù)時序曲線最大限度的接近真實數(shù)據(jù)。具體的賦權(quán)規(guī)則如下:DN值為0的點權(quán)重設(shè)為1，DN值為1的點權(quán)重設(shè)為0．8，而DN值為2和3的點權(quán)重設(shè)為0．2［22］。對某一像元的NDVI數(shù)據(jù)分別利用AG擬合法、DL擬合法和SG濾波法進(jìn)行時序植被指數(shù)重構(gòu)，對比未引入質(zhì)量控制數(shù)據(jù)和引入質(zhì)量控制數(shù)據(jù)的擬合曲線，以此來評價原始數(shù)據(jù)質(zhì)量對不同方法重構(gòu)時序植被指數(shù)的影響。

表1 MOD13A2 v005產(chǎn)品Data Pixel Reliability數(shù)據(jù)集DN值說明［22］

2．2．2 抽樣點數(shù)據(jù)重構(gòu)效果對比

在耕地范圍內(nèi)隨機選擇100個抽樣點并獲取其原始的時間序列植被指數(shù)數(shù)據(jù)，分別用AG擬合、DL擬合和SG濾波3種方法對植被指數(shù)時序曲線進(jìn)行重構(gòu)，將不同方法的重構(gòu)結(jié)果以及原始數(shù)據(jù)進(jìn)行定性對比分析。在由時序植被指數(shù)提取物候信息和應(yīng)用全年或生長季植被指數(shù)平均狀況的研究中，要求重構(gòu)算法在去除噪聲的同時最大限度的保留原始高質(zhì)量數(shù)據(jù)的真實性。因此，該研究從重構(gòu)后的植被指數(shù)時序曲線保持原始植被指數(shù)曲線特征的能力及其與原始數(shù)據(jù)之間的平均差異程度兩個方面進(jìn)行分析。其中高質(zhì)量數(shù)據(jù)點是進(jìn)行時序植被指數(shù)重構(gòu)的重要依據(jù)，它們代表了真實的植被指數(shù)時序曲線的波形特征。為了定量評價不同方法的重構(gòu)效果，計算了重構(gòu)前后的植被指數(shù)時序曲線之間的回歸系數(shù) (R2)和均方根誤差(RMSE)，其中R2反映了二者之間的相關(guān)程度，而RMSE則反映了二者之間的平均差異程度。具體做法是:分別計算每個樣本像元利用不同方法重構(gòu)的植被指數(shù)時序曲線與原始植被指數(shù)曲線之間的R2和RMSE，對所有像元統(tǒng)計這兩個指標(biāo)的平均值，當(dāng)R2值越高或RMSE值越低時表明重構(gòu)的植被指數(shù)時序曲線對原始數(shù)據(jù)的保真性越好。考慮到原始數(shù)據(jù)的質(zhì)量問題，在定量評價時分全年數(shù)據(jù)和高質(zhì)量數(shù)據(jù)進(jìn)行比較，其中高質(zhì)量數(shù)據(jù)是指在Data Pixel Reliability數(shù)據(jù)集中DN值為0和1的像元比重較高時段的數(shù)據(jù)。

2．2．3 研究區(qū)耕地物候特征對比

一般而言，耕地的植被指數(shù)在以年為周期的變化中會經(jīng)歷一個從升高到峰值再到降低的動態(tài)過程［23］。這種總體的變化特征可以描述耕地作物的生長過程。因此，時序植被指數(shù)能夠用來提取耕地物候信息。對時序植被指數(shù)進(jìn)行去噪和平滑等預(yù)處理之后，即得到重構(gòu)的植被指數(shù)時序曲線，在此基礎(chǔ)上就可以開展物候信息提取的工作［24］。該研究基于AG擬合、DL擬合和SG濾波3種方法重構(gòu)的植被指數(shù)時序曲線，分別提取研究區(qū)的耕地物候參數(shù)及其空間格局，通過耕地物候特征間接對比不同的時序植被指數(shù)重構(gòu)方法，作為前面抽樣點對比的補充。提取的耕地物候參數(shù)包括耕地的生長季開始期和結(jié)束期并以此來計算生長季長度，物候參數(shù)提取方法采用比例閾值法［17］。比例閾值法能夠消除不同土壤背景和植被類型的影響，但是在比例參數(shù)的選擇上具有一定的經(jīng)驗性［25］。針對農(nóng)作物等植被生長季NDVI曲線不對稱的問題，可以將生長季開始和結(jié)束的比例參數(shù)設(shè)置為不同值。對于東北三省而言，理論上耕地生長季開始期可以對應(yīng)作物的出苗期，而生長季結(jié)束期則對應(yīng)作物的成熟期［26］。該研究結(jié)合地面觀測數(shù)據(jù)將耕地生長季開始和結(jié)束期的閾值分別設(shè)置為10%和50%。

3 結(jié)果與分析

3.1 原始數(shù)據(jù)質(zhì)量對植被指數(shù)時序曲線重構(gòu)的影響

圖2是對東北三省2011年全年23期Data Pixel Reliability數(shù)據(jù)的DN值進(jìn)行統(tǒng)計的結(jié)果，其中DN值等于0和1的點被認(rèn)為是高質(zhì)量數(shù)據(jù)，而DN值等于2和3的點被認(rèn)為是低質(zhì)量數(shù)據(jù)。從中可以發(fā)現(xiàn)，前5期和后3期影像中除第21期低質(zhì)量數(shù)據(jù)比重略小于高質(zhì)量數(shù)據(jù)之外，其它各期影像低質(zhì)量數(shù)據(jù)的比重都大于高質(zhì)量數(shù)據(jù)，特別是前4期影像中低質(zhì)量數(shù)據(jù)比重遠(yuǎn)遠(yuǎn)大于高質(zhì)量數(shù)據(jù)。即這些時段的影像整體上屬于低質(zhì)量數(shù)據(jù)，可能會對時間序列植被指數(shù)重構(gòu)產(chǎn)生影響。而從第6期影像起高質(zhì)量數(shù)據(jù)的比重開始遠(yuǎn)遠(yuǎn)大于低質(zhì)量數(shù)據(jù)，從第7期一直到第20期影像高質(zhì)量數(shù)據(jù)的比重始終保持在一個較高的水平，同時低質(zhì)量數(shù)據(jù)比重維持在一個非常低的水平。因此，這些時段的影像整體上屬于高質(zhì)量數(shù)據(jù)，是進(jìn)行時間序列植被指數(shù)重構(gòu)的重要依據(jù)，也是該研究中對3種方法的重構(gòu)效果進(jìn)行比較時需要重點分析的數(shù)據(jù)。

圖2 東北三省MOD13A2數(shù)據(jù)產(chǎn)品2011年各期影像數(shù)據(jù)質(zhì)量統(tǒng)計

圖3為分別利用AG擬合、DL擬合和SG濾波3種方法在無質(zhì)量控制和引入質(zhì)量控制的情況下得到的單個像元在2011年的時序NDVI重構(gòu)結(jié)果以及原始數(shù)據(jù)的對比。從圖中可以看出，前5期和后3期的原始植被指數(shù)明顯偏低，與前面對Data Pixel Reliability數(shù)據(jù)DN值的統(tǒng)計結(jié)果分析相對應(yīng)。3種方法引入質(zhì)量控制前后的重構(gòu)曲線之間差異明顯，主要表現(xiàn)為:AG擬合和DL擬合的重構(gòu)結(jié)果在低質(zhì)量數(shù)據(jù)時段更加偏離原始曲線，而在高質(zhì)量數(shù)據(jù)時段更加符合原始曲線，且二者的重構(gòu)曲線除后3期外在引入質(zhì)量控制后幾乎趨于一致;SG濾波的重構(gòu)結(jié)果在低質(zhì)量數(shù)據(jù)時段的情況與前兩種方法類似，而在高質(zhì)量數(shù)據(jù)時段引入質(zhì)量控制前后變化不大。原始數(shù)據(jù)質(zhì)量的影響使得無質(zhì)量控制的重構(gòu)曲線在上升和下降階段的最小值明顯低于引入質(zhì)量控制的重構(gòu)曲線，這樣就可能會給下一步的趨勢分析和信息提取帶來誤差。考慮到原始數(shù)據(jù)質(zhì)量影響的問題，該研究在后續(xù)利用3種方法對NDVI時序數(shù)據(jù)進(jìn)行重構(gòu)的過程中統(tǒng)一引入質(zhì)量控制數(shù)據(jù)。

圖3 單個像元引入質(zhì)量控制前后3種方法重構(gòu)的時序植被指數(shù)以及原始數(shù)據(jù)對比

3.2 3種重構(gòu)方法對原始數(shù)據(jù)的保真性對比

圖4中4條曲線分別為100個抽樣點2011年的原始植被指數(shù)和經(jīng)過3種方法重構(gòu)后的植被指數(shù)取平均值得到的時間序列曲線。從整體上看，AG擬合和DL擬合得到的結(jié)果非常接近，而SG濾波得到的結(jié)果則與上述兩種方法存在差異。在低質(zhì)量數(shù)據(jù)時段，即影像的前5期和后3期時段，3種方法得到的重構(gòu)曲線相對于原始的植被指數(shù)曲線都有所提升，但SG濾波結(jié)果的提升幅度小于前兩種方法。在高質(zhì)量數(shù)據(jù)時段，AG擬合和DL重構(gòu)曲線的上升和下降階段與原始曲線吻合度較高，而在曲線峰值附近擬合曲線較原始曲線有所提升;SG濾波得到的結(jié)果在重構(gòu)曲線的上升和下降階段相對原始曲線分別出現(xiàn)左偏和右移，而在曲線峰值附近較原始曲線有所降低。其中SG濾波重構(gòu)曲線的左偏和右移將會導(dǎo)致利用比例閾值法提取的生長季開始期提前和結(jié)束期推遲，從而使得計算的生長季長度延長，另外擬合峰值的降低則進(jìn)一步擴大了這種誤差的范圍。

圖4 3種方法重構(gòu)得到的植被指數(shù)均值時序曲線與原始植被指數(shù)時序曲線對比

圖5為分別用AG擬合、DL擬合和SG濾波3種方法對100個抽樣點2011年全年23期以及代表高質(zhì)量數(shù)據(jù)的第6～20期的植被指數(shù)進(jìn)行重構(gòu)前后的時序曲線之間的回歸系數(shù)和均方根誤差對比。從全年數(shù)據(jù)來看，SG濾波的R2值最高，AG擬合次之，DL擬合最低;RMSE值SG濾波最低，AG擬合居中，DL擬合最高。說明SG濾波法對原始數(shù)據(jù)的保真性最好，AG擬合次之，DL擬合最差，但值得注意的是這3種方法的R2和RMSE值之間差距并不大。對高質(zhì)量數(shù)據(jù)而言，整體上AG擬合和DL擬合的R2和RMSE值非常接近，而SG濾波與上述兩種方法之間差距明顯。其中AG擬合和DL擬合的R2值明顯高于SG濾波，且AG擬合比DL擬合略高;而AG擬合和DL擬合的RMSE值則明顯低于SG濾波，且AG擬合比DL擬合略低。由此說明，利用SG濾波重構(gòu)的植被指數(shù)時序曲線對原始高質(zhì)量數(shù)據(jù)的保真性最差，AG擬合與DL擬合重構(gòu)的植被指數(shù)時序曲線對原始高質(zhì)量數(shù)據(jù)的保真性基本上相接近且AG擬合略占優(yōu)勢。

圖5 3種方法重構(gòu)前后植被指數(shù)時序曲線高質(zhì)量數(shù)據(jù)的回歸系數(shù)和均方根誤差對比

3.3 基于不同重構(gòu)方法的耕地生長季長度空間格局

該研究分別利用AG擬合、DL擬合和SG濾波3種方法對整個東北三省2010～2012年共69期的NDVI時序數(shù)據(jù)進(jìn)行重構(gòu)，然后由比例閾值法提取了2011年該地區(qū)耕地的生長季開始期和結(jié)束期并以此來計算生長季長度。

3．3．1 耕地生長季長度的空間分布特征比較

圖6為分別基于3種時序植被指數(shù)重構(gòu)方法得到的東北三省2011年耕地生長季長度的空間分布特征?？梢钥闯?，基于AG擬合和DL擬合提取的耕地生長季長度空間分布基本一致，而基于SG濾波的提取結(jié)果與前兩種方法差異明顯。在基于AG擬合和DL擬合的提取結(jié)果中，吉林省西部、黑龍江省松嫩平原西北部以及三江平原的部分地區(qū)耕地生長季長度一般為90～120天，三大平原其它大部分地區(qū)的耕地生長季長度為120～135天，另外黑龍江省東部和南部的零星區(qū)域、吉林省中南部以及遼寧省東部和西部的部分地區(qū)耕地生長季長度為135～150天，該結(jié)果與李正國等［26］的研究結(jié)果基本一致;而基于SG濾波的提取結(jié)果中耕地生長季長度較前兩種方法則明顯延長，大部分區(qū)域的耕地生長季長度在150天以上，只有少部分區(qū)域為120～150天。這與之前定性分析3種方法對原始高質(zhì)量數(shù)據(jù)的保真性時的發(fā)現(xiàn)相對應(yīng)，即SG濾波重構(gòu)曲線在上升和下降階段相對于原始曲線分別出現(xiàn)左偏和右移以及在峰值附近有所降低，由此導(dǎo)致利用比例閾值法提取的生長季開始期提前和結(jié)束期推遲，進(jìn)而使得計算的生長季長度延長。

圖6 東北三省2011年耕地生長季長度空間分布特征(從左到右依次為基于AG擬合、DL擬合和SG濾波得到的結(jié)果)

3．3．2 不同生長季長度區(qū)域的重構(gòu)曲線比較

根據(jù)東北三省2011年耕地生長季長度的空間分布情況，從中分別選擇生長季短、中和長3種類型的隨機抽樣點，將原始NDVI曲線和由3種方法得到的重構(gòu)曲線進(jìn)行比較分析，如圖7?？梢钥闯?，在耕地生長季長度較短的區(qū)域3種方法得到的重構(gòu)曲線在低質(zhì)量數(shù)據(jù)時段都較原始曲線有所抬升，其中AG擬合和DL擬合在前5期略有差異但后3期幾乎一致，SG濾波與上述兩種方法則差異明顯;在高質(zhì)量數(shù)據(jù)時段AG擬合和DL擬合結(jié)果與原始曲線比較吻合，但SG濾波重構(gòu)曲線在上升和下降階段較原始曲線分別出現(xiàn)左偏和右移。在耕地生長季長度居中的區(qū)域總體上與短生長季的區(qū)域相類似，不同的是低質(zhì)量數(shù)據(jù)時段AG擬合和DL擬合在前5期幾乎一致而在后3期略有差異，同時SG濾波與上述兩種方法的差異增大;高質(zhì)量數(shù)據(jù)時段AG擬合和DL擬合的結(jié)果更加符合原始曲線，SG濾波重構(gòu)曲線在上升階段的左偏程度減小，擬合峰值也有所降低。在耕地生長季長度較長的區(qū)域原始數(shù)據(jù)質(zhì)量影響不大，3種方法得到的重構(gòu)曲線非常接近且與原始曲線的吻合程度較高。

圖7 不同生長季長度的隨機抽樣點3種方法重構(gòu)的時序植被指數(shù)及原始數(shù)據(jù)對比

4 結(jié)論與討論

4.1 結(jié)論

該研究以東北三省為研究區(qū)域，利用MODIS植被指數(shù)產(chǎn)品對比了AG擬合、DL擬合和SG濾波3種常見的時序植被指數(shù)重構(gòu)方法。主要從原始數(shù)據(jù)質(zhì)量對時間序列植被指數(shù)重構(gòu)的影響、不同重構(gòu)方法對原始數(shù)據(jù)的保真性以及基于不同重構(gòu)方法的耕地生長季長度空間格局3個方面進(jìn)行分析。研究結(jié)果表明:

(1)原始數(shù)據(jù)質(zhì)量對利用3種方法進(jìn)行時序植被指數(shù)重構(gòu)都存在影響，并且對AG擬合和DL擬合重構(gòu)結(jié)果的影響大于SG濾波。主要表現(xiàn)為冰雪覆蓋和云掩造成的植被指數(shù)值降低，導(dǎo)致無質(zhì)量控制的重構(gòu)曲線在上升和下降階段的最小值明顯低于引入質(zhì)量控制的重構(gòu)曲線，這樣就可能給下一步的趨勢分析和信息提取帶來誤差。

(2)從整體上看，AG擬合和DL擬合得到的結(jié)果非常接近，SG濾波得到的結(jié)果則與上述兩種方法存在差異。從全部原始數(shù)據(jù)來看，3種方法中SG濾波重構(gòu)的時序植被指數(shù)的保真性最好，AG擬合次之，DL擬合最差;但是對高質(zhì)量數(shù)據(jù)而言SG濾波的保真性則最差，AG擬合與DL擬合的保真性相接近且AG擬合略占優(yōu)勢。

(3)基于AG擬合和DL擬合的耕地生長季長度在空間分布特征基本上一致，并且與前人在該地區(qū)利用VGT植被指數(shù)得到的研究結(jié)果相符合?；赟G濾波的耕地生長季長度空間格局與上述兩種方法差異顯著，大部分區(qū)域的耕地生長季長度明顯延長。在生長季較短和居中的區(qū)域SG濾波重構(gòu)曲線在上升和下降階段相對原始曲線分別出現(xiàn)左偏和右移，而在生長季長度較長的區(qū)域與原始曲線的吻合程度較高。

總體上看，AG擬合和DL擬合重構(gòu)結(jié)果在各個方面都表現(xiàn)出非常相似的特點，且基于這兩種方法得到的耕地物候特征更加符合實際情況，而SG濾波則與前兩者明顯不同。通過上述研究發(fā)現(xiàn)，在東北三省利用AG擬合、DL擬合和SG濾波進(jìn)行耕地時序植被指數(shù)重構(gòu)時需要引入質(zhì)量控制數(shù)據(jù)，AG擬合和DL擬合在該地區(qū)耕地時序植被指數(shù)重構(gòu)時表現(xiàn)相似且比SG濾波更加適宜。

4.2 分析和討論

造成AG擬合、DL擬合和SG濾波3種方法重構(gòu)結(jié)果差異性的原因主要與這些方法的自身特點密切相關(guān)。AG擬合和DL擬合屬于函數(shù)擬合方法，是從局部到整體的處理方法，受到擬合函數(shù)形式的影響較大且得到的重構(gòu)曲線更加平滑;SG濾波屬于時域濾波方法，是一種局部處理的方法，重構(gòu)曲線的平滑程度與濾波窗口的大小有關(guān)并且受原始數(shù)據(jù)的影響較大。AG擬合和DL擬合在原理上基本相同，只是局部擬合函數(shù)形式有所差異，因此這兩種方法的重構(gòu)曲線比較接近，由此得到的耕地物候空間格局也非常相似。SG濾波的原理與前兩種方法截然不同，因此導(dǎo)致了重構(gòu)曲線以及耕地物候空間格局的明顯差異，其中SG濾波重構(gòu)曲線在上升和下降階段的左偏和右移是所提取的耕地生長季長度明顯延長的直接原因。除此之外數(shù)據(jù)情況也與3種方法重構(gòu)結(jié)果的差異性有關(guān)，該研究使用的植被指數(shù)來自Terra衛(wèi)星上的MODIS傳感器，雖然經(jīng)過了MVC方法合成，但殘余的云掩和冰雪覆蓋影響仍然存在，另外1km的空間分辨率和16天的時間分辨率也造成了部分細(xì)節(jié)信息的損失。

AG擬合、DL擬合和SG濾波3種時序植被指數(shù)重構(gòu)方法各具優(yōu)勢，但是在不同地區(qū)和不同研究目的時仍然需要進(jìn)行適宜性分析。理論上以強調(diào)局部細(xì)節(jié)為特點的SG濾波法對原始數(shù)據(jù)的保真性應(yīng)該最高，然而，在東北三省由于受到冰雪覆蓋及云掩導(dǎo)致的原始數(shù)據(jù)質(zhì)量下降影響，使其對時序植被指數(shù)重構(gòu)的效果大打折扣，進(jìn)而影響到物候信息的提取，以整體擬合為特點的AG擬合和DL擬合則能夠在一定程度上避免這種影響。另外在用比例閾值法提取耕地物候信息時比例參數(shù)的選擇非常關(guān)鍵，20%是目前應(yīng)用比較普遍的比例參數(shù)，然而大部分研究都是以自然植被為對象的。盡管該研究通過結(jié)合地面觀測數(shù)據(jù)將耕地生長季開始期和結(jié)束期的比例參數(shù)分別設(shè)置為10%和50%，但在缺乏地面數(shù)據(jù)的情況下如何為受人為因素影響較大的耕地選擇一個能夠準(zhǔn)確獲取耕地物候信息的比例參數(shù)仍然值得商榷。

［1］王正興，劉闖，Alfredo Huete．植被指數(shù)研究進(jìn)展:從AVHRR-NDVI到MODIS-EVI．生態(tài)學(xué)報，2003，23(5):979～987

［2］Zhang X，F(xiàn)riedl M A，Schaaf C B，et al．Monitoring vegetation phenology using MODIS．Remote Sensing of Environment，2003，84(3):471～475

［3］李儒，張霞，劉波，等．遙感時間序列數(shù)據(jù)濾波重建算法發(fā)展綜述．遙感學(xué)報，2009，13(2):335～341

［4］吳文斌，楊鵬，唐華俊，等．兩種NDVI時間序列數(shù)據(jù)擬合方法比較．農(nóng)業(yè)工程學(xué)報，2009，25(11):183～188

［5］Taddei R．Maximum Value Interpolated(MVI):A Maximum Value Composite method improvement in vegetation index profiles analysis．International Journal of Remote Sensing，1997，18(11):2365～2370

［6］Beck P S A，Atzberger C，H?gda K A，et al．Improved monitoring of vegetation dynamics at very high latitudes:A new method using MODIS NDVI．Remote Sensing of Environment，2006，100(3):321 ～334

［7］Hird J N，McDermid G J．Noise reduction of NDVI time series:An empirical comparison of selected techniques．Remote Sensing of Environment，2009，113(1):248～258

［8］曹云鋒，王正興，鄧芳萍．3種濾波算法對NDVI高質(zhì)量數(shù)據(jù)保真性研究．遙感技術(shù)與應(yīng)用，2010，25(1):118～125

［9］宋春橋，游松財，柯靈紅，等．藏北地區(qū)三種時序NDVI重建方法與應(yīng)用分析．地球信息科學(xué)學(xué)報，2011，13(1):133～143

［10］李正國，唐華俊，楊鵬，等．植被物候特征的遙感提取與農(nóng)業(yè)應(yīng)用綜述．中國農(nóng)業(yè)資源與區(qū)劃，2012，33(5):20～28

［11］任國玉，初子瑩，周雅清，等．中國氣溫變化研究最新進(jìn)展．氣候與環(huán)境研究，2005，(4):701～716

［12］劉志娟，楊曉光，王文峰，等．氣候變化背景下我國東北三省農(nóng)業(yè)氣候資源變化特征．應(yīng)用生態(tài)學(xué)報，2009(09):2199～2206

［13］李正國，楊鵬，唐華俊，等．氣候變化背景下東北三省主要作物典型物候期變化趨勢分析．中國農(nóng)業(yè)科學(xué)，2011，44(20):4180～4189

［14］李正國，唐華俊，楊鵬，等．東北三省耕地物候期對熱量資源變化的響應(yīng)．地理學(xué)報，2011，66(7):928～939

［15］李杭燕，頡耀文，馬明國．時序NDVI數(shù)據(jù)集重建方法評價與實例研究．遙感技術(shù)與應(yīng)用，2009，24(5):596～602

［16］孫華生，徐愛功，林卉，等．基于不同算法的時間序列植被指數(shù)去噪效果分析．江蘇農(nóng)業(yè)科學(xué)，2012，40(5):375～379

［17］J?nsson P，Eklundh L．Seasonality extraction by function fitting to time-series of satellite sensor data．IEEE Transactions on Geoscience and Remote Sensing，2002，40(8):1824～1832

［18］J?nsson P，Eklundh L．TIMESAT—a program for analyzing time-series of satellite sensor data．Computers ＆ Geosciences，2004，30(8):833～845

［19］Chen J，J?nsson P，Tamura M，et al．A simple method for reconstructing a high-quality NDVI time-series data set based on the Savitzky-Golay filter．Remote Sensing of Environment，2004，91:332 ～344

［20］蔡天凈，唐瀚．Savitzky-Golay平滑濾波器的最小二乘擬合原理綜述．?dāng)?shù)字通信，2011，(1):63～68

［21］顧娟，李新，黃春林．NDVI時間序列數(shù)據(jù)集重建方法述評．遙感技術(shù)與應(yīng)用，2006，21(4):391～395

［22］宋春橋，柯靈紅，游松財，等．基于TIMESAT的3種時序NDVI擬合方法比較研究——以藏北草地為例．遙感技術(shù)與應(yīng)用，2011，26(2):147～155

［23］吳文斌，楊鵬，唐華俊，等．過去20年中國耕地生長季起始期的時空變化．生態(tài)學(xué)報，2009，29(4):1777～1786

［24］胡硯霞，黃進(jìn)良，王立輝．基于MODIS NDVI時序數(shù)據(jù)的丹江口庫區(qū)物候格局研究．中國農(nóng)業(yè)資源與區(qū)劃，2012，33(4):25～30

［25］吳文斌，楊鵬，唐華俊，等．基于NDVI數(shù)據(jù)的華北地區(qū)耕地物候空間格局．中國農(nóng)業(yè)科學(xué)，2009，42(2):552～560

［26］李正國，唐華俊，楊鵬，等．基于時序植被指數(shù)的東北地區(qū)耕地生長季特征識別與應(yīng)用研究．北京大學(xué)學(xué)報 (自然科學(xué)版)，2011，47(5):882～892