環(huán)境星NDVI時(shí)間序列重構(gòu)方法研究

環(huán)境星NDVI時(shí)間序列重構(gòu)方法研究

李天祺1，2，朱秀芳1，2，潘耀忠1，2，劉憲鋒1，2

(1.北京師范大學(xué) 地表過(guò)程與資源生態(tài)國(guó)家重點(diǎn)實(shí)驗(yàn)室，北京 100875；2.北京師范大學(xué) 資源學(xué)院，北京 100875)

摘要：利用目前時(shí)間序列曲線重構(gòu)中較為常用的非對(duì)稱(chēng)高斯函數(shù)擬合、Double-Logistic曲線擬合、S-G濾波和時(shí)間序列諧波分析法對(duì)環(huán)境星NDVI時(shí)間序列進(jìn)行重構(gòu)處理。分析了上述4種植被指數(shù)時(shí)間序列重構(gòu)方法對(duì)環(huán)境星數(shù)據(jù)的適用性。實(shí)驗(yàn)結(jié)果表明，對(duì)于環(huán)境星數(shù)據(jù)，在4種方法中非對(duì)稱(chēng)高斯函數(shù)擬合、Double-Logistic曲線擬合法更適用于對(duì)植被地物的時(shí)間序列進(jìn)行重構(gòu)，對(duì)照參考數(shù)據(jù)，其重構(gòu)曲線對(duì)植被物候的表達(dá)有較高的一致性；而時(shí)間序列諧波分析法對(duì)原始數(shù)據(jù)的擾動(dòng)最小，適用于非植被地物的時(shí)間序列重構(gòu)；S-G濾波在4種方法中的重構(gòu)效果最差。

關(guān)鍵詞：時(shí)間序列；環(huán)境星；非對(duì)稱(chēng)高斯函數(shù)擬合；Double-Logistic曲線擬合；S-G濾波；時(shí)間序列諧波分析法

doi:10.3969/j.issn.1000-3177.2015.01.010

中圖分類(lèi)號(hào)：TP79文獻(xiàn)標(biāo)識(shí)碼：A

收稿日期：2014-02-21修訂日期：2014-04-11

基金項(xiàng)目：國(guó)土資源部公益性行業(yè)科研

作者簡(jiǎn)介：楊進(jìn)生(1959～)，男，高級(jí)工程師，主要從事水文地質(zhì)環(huán)境地質(zhì)遙感技術(shù)方法研究。

NDVI Time-series Reconstruction Methods of China’s

HJ Satellite Imagery

LI Tian-qi1,2,ZHU Xiu-fang1,2,PAN Yao-zhong1,2, LIU Xian-feng1,2

(1.StateKeyLaboratoryofEarthSurfaceProcessesandResourceEcology，BeijingNormalUniversity，Beijing100875；

2.CollegeofResourcesScienceandTechnology，BeijingNormalUniversity，Beijing100875)

Abstract：In this paper，we compared the performance of four NDVI time-series reconstruction methods，including asymmetric Guassian function fitting (A-G)，double logistic function fitting (D-L)，Savitzky-Golay filtering method (S-G) and harmonic analysis of time series (Hants)，on processing the imagery of China’s HJ satellites. Results show that A-G and D-L are more suitable for reconstruction of vegetation features，comparing with the reference data. The reconstruction curves of A-G and D-L have high consistency in expression vegetation phenology. Hants brings minimal disturbance of the original data and fits for reconstruction of architectural features and water features. S-G performs worst.

Key words：time-series；China’s HJ satellite；asymmetric Guassian function fitting；double logistic function fitting；Savitzky-Golay filtering；Hants

1引言

歸一化植被指數(shù)(Normalized Difference Vegetation Index，NDVI)是反映植被生長(zhǎng)狀態(tài)及植被覆蓋度的最佳指示因子，也是季節(jié)變化和人類(lèi)活動(dòng)影響的重要指示器[1]。在遙感研究中，隨著高時(shí)間分辨率數(shù)據(jù)(NOAA/AVHRR、TERRA/MODIS、SPOT/VEGETATION)的應(yīng)用，歸一化植被指數(shù)時(shí)間序列作為反映地表狀況的連續(xù)信號(hào)，已經(jīng)廣泛地應(yīng)用于環(huán)境動(dòng)態(tài)變化監(jiān)測(cè)[2-3]、植被物候信息提取[4-5]、土地覆蓋/利用分類(lèi)及監(jiān)測(cè)[6-7]等諸多領(lǐng)域，是生產(chǎn)研究的重要數(shù)據(jù)源之一[8]。

然而，由于傳感器角度變化、云或霾的干擾、數(shù)據(jù)傳輸誤差、二向性反射或地面冰雪覆蓋的影響，使得NDVI時(shí)間序列曲線有明顯的突升或突降，導(dǎo)致植被物候信息難以辨別[9-10]。NDVI時(shí)間序列中噪聲的存在限制了該數(shù)據(jù)的深入應(yīng)用，因此在應(yīng)用前必須有效地去除噪聲、重構(gòu)時(shí)間序列[11]。對(duì)此，許多關(guān)于NDVI時(shí)間序列降噪、重構(gòu)的方法被提出、評(píng)價(jià)以及用于實(shí)際研究[12]。目前，主要的NDVI時(shí)間序列重構(gòu)方法可分為：基于閾值去噪法(如最佳指數(shù)斜率提取法(BISE))、基于濾波的平滑方法(如Savitzky-Golay濾波法、時(shí)間序列諧波分析法(Hants)、傅里葉變換)和非線性擬合法(如Double-Logistic曲線擬合法和非對(duì)稱(chēng)高斯函數(shù)擬合法)，這些方法已廣泛應(yīng)用于全球不同區(qū)域的研究中[13]。

不同的重構(gòu)方法有其各自?xún)?yōu)點(diǎn)，對(duì)不同區(qū)域或不同數(shù)據(jù)源的適用性也不同。因此，國(guó)內(nèi)外學(xué)者在方法重構(gòu)效果方面做了大量的對(duì)比研究。例如Beck以芬蘭、瑞典、挪威的交界區(qū)域?yàn)檠芯繀^(qū)，利用Double-Logistic曲線擬合法、非對(duì)稱(chēng)性高斯函數(shù)擬合法和傅立葉變換法對(duì)MODIS NDVI產(chǎn)品進(jìn)行重構(gòu)比較研究，認(rèn)為Double-Logistic曲線擬合和非對(duì)稱(chēng)性高斯函數(shù)擬合的重構(gòu)效果相似，且均優(yōu)于傅立葉變換的方法[14]；Jennifer選擇加拿大洛磯山脈中的一個(gè)典型區(qū)域作為研究區(qū)，利用MODIS NDVI產(chǎn)品對(duì)比6種NDVI時(shí)間序列重構(gòu)方法，認(rèn)為Double-Logistic曲線擬合與非對(duì)稱(chēng)高斯函數(shù)擬合在總體上優(yōu)于其他4種濾波方法[15]；Chen利用SPOT Vegetation產(chǎn)品以中國(guó)為研究區(qū)，對(duì)多類(lèi)地物的NDVI時(shí)間序列進(jìn)行重構(gòu)研究，得出通過(guò)改進(jìn)的Savitzky-Golay濾波法在NDVI曲線重構(gòu)上優(yōu)于最佳指數(shù)斜率提取法(BISE)與傅里葉變換(FT)的結(jié)論[16]；Julien利用GIMMS/NDVI數(shù)據(jù)集對(duì)比IDR(Iterative Interpolation for Data Reconstruction)方法與HANTS、Double-Logistic曲線擬合法的重構(gòu)效果，認(rèn)為IDR方法在3種方法中的重構(gòu)效果最好，但在強(qiáng)噪聲中對(duì)NDVI低值存在過(guò)度擬合[17]。

環(huán)境星全稱(chēng)為環(huán)境與災(zāi)害監(jiān)測(cè)預(yù)報(bào)小衛(wèi)星星座，A、B星(HJ-1A/1B星)，于2008年9月投入使用。重訪周期為2天，所搭載的CCD相機(jī)具有可見(jiàn)光及近紅外4個(gè)波段，空間分辨率為30m，具有寬幅、中高分辨率的特點(diǎn)，適用于區(qū)域的大范圍中尺度覆蓋監(jiān)測(cè)。環(huán)境星對(duì)我國(guó)生態(tài)、環(huán)境變化和災(zāi)害監(jiān)測(cè)，反映生態(tài)環(huán)境發(fā)展過(guò)程，以及生態(tài)環(huán)境發(fā)展變化趨勢(shì)的預(yù)測(cè)均有著重要意義[18]。環(huán)境星反演的NDVI時(shí)序序列數(shù)據(jù)已經(jīng)被廣泛的應(yīng)用于植被變化監(jiān)測(cè)、農(nóng)作物面積測(cè)量及長(zhǎng)勢(shì)監(jiān)測(cè)、災(zāi)害監(jiān)測(cè)以及土地覆蓋/利用的劃分中。然而，目前尚無(wú)有關(guān)各類(lèi)重構(gòu)方法在環(huán)境星NDVI數(shù)據(jù)重構(gòu)研究中的適用性的報(bào)道。為此，本文綜合目前國(guó)內(nèi)外研究成果，選取當(dāng)前4類(lèi)主要的重構(gòu)方法，分析適用于環(huán)境星數(shù)據(jù)的植被指數(shù)時(shí)間序列重構(gòu)方法，以改善環(huán)境星在植被動(dòng)態(tài)變化監(jiān)測(cè)、土地覆蓋/利用分類(lèi)等方面的應(yīng)用效果。

2研究區(qū)及數(shù)據(jù)

2.1研究區(qū)

本研究以北京為研究區(qū)。北京位于華北平原西北邊緣，地處115°20′E～117°30′E，39°28′N(xiāo)～41°05′N(xiāo)，地勢(shì)西北高、東南低。西部、北部和東北部三面環(huán)山，東南部是向東傾斜的平原。北京市總面積16807km2，其中山區(qū)約占總面積的62%，平原約占38%。研究區(qū)屬北溫帶大陸性季風(fēng)氣候，具有分明的四季劃分。

2.2實(shí)驗(yàn)數(shù)據(jù)

研究中所用數(shù)據(jù)包括環(huán)境星多光譜數(shù)據(jù)和MODIS產(chǎn)品(MOD13Q1)。其中環(huán)境星數(shù)據(jù)為2012年2月到11月北京地區(qū)環(huán)境星星座多光譜影像(近紅外波段：760nm～903nm，紅波段：630nm～693nm，綠波段：520nm～603nm，藍(lán)波段：430nm～523nm)，影像大小為5926×5987像元，空間分辨率為30m，一共為28景，平均時(shí)間間隔為10天；MODIS植被指數(shù)產(chǎn)品(MOD13Q1)為基于16天最大值合成法的植被指數(shù)數(shù)據(jù)，空間分辨率為250 m，該產(chǎn)品經(jīng)過(guò)幾何糾正和大氣校正，產(chǎn)品中NDVI數(shù)據(jù)和QA(Quality Assessment)用于本實(shí)驗(yàn)研究。

圖1　實(shí)驗(yàn)數(shù)據(jù)

2.3環(huán)境星數(shù)據(jù)預(yù)處理

環(huán)境星數(shù)據(jù)預(yù)處理包括輻射定標(biāo)、大氣校正、幾何糾正和植被指數(shù)時(shí)間序列生成。首先，對(duì)環(huán)境星數(shù)據(jù)的輻射定標(biāo)采用2011年中國(guó)資源衛(wèi)星應(yīng)用中心公布的HJ-A/B星絕對(duì)輻射定標(biāo)系數(shù)，將傳感器DN值標(biāo)定為反射率；其次，大氣校正采用FLAASH模型進(jìn)行校正，該模型在環(huán)境星數(shù)據(jù)的大氣校正中能夠很好地消除大氣的影響[19]；幾何糾正以帶有精確地理坐標(biāo)的TM影像為參考對(duì)環(huán)境星影像進(jìn)行配準(zhǔn)，其中TM影像的觀測(cè)時(shí)間為2011年7月26日，配準(zhǔn)誤差控制在1個(gè)像元以?xún)?nèi)；最后，對(duì)28景影像逐像元計(jì)算NDVI，并進(jìn)行波段融合，但由于大面積云污染等天氣原因及傳感器時(shí)間分辨率的限制，現(xiàn)有的高質(zhì)量環(huán)境星NDVI數(shù)據(jù)并非等時(shí)間間隔，考慮植被在短時(shí)間(10天～16天)的生長(zhǎng)特性，通過(guò)線性插值的方法[12]，將現(xiàn)有28期NDVI數(shù)據(jù)插值為等時(shí)間間隔的27個(gè)波段，時(shí)間間隔為11天。

3研究方法

3.1典型地物原始時(shí)序曲線

研究中選取北京地區(qū)的6種典型地物進(jìn)行分析，分別為林地、草地、單季作物、雙季作物、水體、建筑，其中單季作物以春小麥為主，雙季作物以冬小麥-夏玉米為主。在本研究中，參考北京地區(qū)高分辨率遙感影像及2012年農(nóng)作物目視解譯分類(lèi)結(jié)果，每類(lèi)地物選取50個(gè)～80個(gè)樣本像元，共選取390個(gè)像元，并對(duì)每類(lèi)樣點(diǎn)求其N(xiāo)DVI均值，構(gòu)建該類(lèi)別的原始時(shí)序曲線并對(duì)其進(jìn)行重構(gòu)，原始曲線如圖2所示。

圖2　環(huán)境星原始NDVI時(shí)序曲線

3.2典型地物時(shí)序曲線的重構(gòu)

3.2.1非對(duì)稱(chēng)高斯函數(shù)擬合法(Asymmetric Gaussian model，A-G)

非對(duì)稱(chēng)高斯函數(shù)擬合法是基于高斯函數(shù)的分段最小二乘擬合算法，由J?nsson等在2002年提出[4]，該方法的擬合函數(shù)形式如下：

f(t) =f(t;c1,c2,a1,a2,a3,a4,a5)=

c1+c2g(t;a1,a2,a3,a4,a5)

(1)

其中，g(t;a1,a2,a3,a4,a5)為高斯函數(shù)：

g(t;a1,a2,a3,a4,a5)=

(2)

其中，c1,c2決定曲線的基準(zhǔn)和幅度；a1為曲線最大值或最小值在時(shí)間軸上的位置；,a2,a3,a4,a5分別決定左右半曲線的寬度和陡峭度。非對(duì)稱(chēng)高斯函數(shù)擬合法的處理過(guò)程分為區(qū)間提取、局部擬合和整體連接3個(gè)步驟[20]。

3.2.2Double-Logistic曲線擬合法(D-L)

Double-Logistic曲線擬合法是由Beck等[20]于2006年提出的一種非線性最小二乘擬合算法，與非對(duì)稱(chēng)高斯函數(shù)擬合法類(lèi)似，也屬于分段擬合。其函數(shù)形式如下：

NDVI=wNDVI+(mNDVI-wNDVI)×

(3)

其中，wNDVI和mNDVI分別為冬季NDVI和全年NDVI的最大值，S、A代表曲線上升拐點(diǎn)和下降拐點(diǎn)，mA和mS分別代表這兩點(diǎn)處NDVI的上升和下降速率。在處理過(guò)程上，Double-Logistic曲線擬合法與非對(duì)稱(chēng)高斯函數(shù)法類(lèi)似，即利用曲線極值將時(shí)間序列分成多個(gè)區(qū)間，對(duì)每個(gè)區(qū)間進(jìn)行局部擬合，最后對(duì)擬合出的各段曲線進(jìn)行整體連接。

3.2.3Savitzky-Golay濾波法(S-G)

Savitzky-Golay濾波法最初由Savitzky和Golay于1964年提出，后由Chen等[16]對(duì)其進(jìn)行改進(jìn)，是一種基于最小二乘卷積的移動(dòng)平均算法，其權(quán)重取決于在一個(gè)濾波窗口范圍內(nèi)做多項(xiàng)式最小二乘擬合的多項(xiàng)式次數(shù)[13]。其基本公式如下：

(4)

其中，Y為原始NDVI值，Y*為平滑后的NDVI值，Ci是從過(guò)濾窗口首部開(kāi)始的第i個(gè)NDVI值的權(quán)重，N是窗口寬度，且N=2m+1。其中窗口寬度與多項(xiàng)式擬合的階數(shù)由人為設(shè)定。Savitzky-Golay濾波法的處理過(guò)程如下[16]：

①對(duì)時(shí)間序列的長(zhǎng)期變化趨勢(shì)進(jìn)行擬合，將原始NDVI值分為“真值”和“假值”兩類(lèi)；

②通過(guò)局部循環(huán)迭代，使得原始值中的“假值”點(diǎn)被Savitzky-Golay重構(gòu)值替代，與“真值”點(diǎn)構(gòu)成新的NDVI曲線；

③重復(fù)上述過(guò)程，直至擬合結(jié)果足夠接近NDVI時(shí)序曲線的上包絡(luò)線。

3.2.4時(shí)間序列諧波分析法(Hants)

時(shí)間序列諧波分析法(Harmonic analysis of time series，Hants)是基于傅里葉變換改進(jìn)的一種時(shí)序曲線重構(gòu)方法[21]，通過(guò)傅里葉變換將時(shí)序曲線表示為不同相位、幅度和頻率的正弦函數(shù)組合，其表達(dá)式為：

(5)

其中，A0為諧波的余項(xiàng)，等于時(shí)序曲線的平均值；Aj為各諧波的振幅；ωj=2jπ/N為各諧波的頻率，N為時(shí)間序列的長(zhǎng)度；θj為各諧波的初始相位；m=N-1為諧波個(gè)數(shù)。

原始時(shí)序曲線經(jīng)傅里葉變換后，Hants方法通過(guò)對(duì)新序列中幾個(gè)最為顯著頻率的傅里葉序列進(jìn)行最小二乘擬合對(duì)原始曲線重構(gòu)，該過(guò)程經(jīng)多次迭代直至重構(gòu)曲線中沒(méi)有被污染點(diǎn)或達(dá)到某一迭代閾值。擬合過(guò)程中，原始序列中的低頻部分被保留，高頻部分作為噪聲被去除。Hants算法在對(duì)時(shí)序曲線處理時(shí)，需要設(shè)置顯著頻率、誤差閾值、最大刪除點(diǎn)數(shù)及有效數(shù)據(jù)范圍。

3.3重構(gòu)效果評(píng)價(jià)方法

本研究選取北京地區(qū)的6種典型地物進(jìn)行分析，分別為林地、草地、單季作物、雙季作物、水體、建筑，其中單季作物以春小麥為主，雙季作物以冬小麥-夏玉米為主。針對(duì)植被地物與非植被地物分別采取物候信息一致性與曲線相關(guān)性?xún)煞N評(píng)價(jià)方法對(duì)重構(gòu)效果進(jìn)行評(píng)價(jià)。

3.3.1物候信息一致性評(píng)價(jià)

MODIS植被指數(shù)產(chǎn)品經(jīng)MVC方法的初步處理可以有效地去除云污染和擬合缺損[8]，且MODIS植被指數(shù)產(chǎn)品在植被變化監(jiān)測(cè)及物候信息提取等研究中，其物候信息反映的真實(shí)性在國(guó)內(nèi)外有著廣泛的認(rèn)可。在本實(shí)驗(yàn)中林地、草地樣本的分布具有較大區(qū)域的一致性，因此可排除混合像元對(duì)NDVI值造成的影響。在實(shí)驗(yàn)中選取高質(zhì)量MODIS NDVI數(shù)據(jù)，可較為真實(shí)地反映地表林地與草地的NDVI變化；之后，利用同種重構(gòu)方法對(duì)高質(zhì)量MODIS NDVI時(shí)序曲線與相同地點(diǎn)的環(huán)境星NDVI時(shí)序曲線進(jìn)行重構(gòu)，并分別提取重構(gòu)結(jié)果的物候信息，對(duì)二者進(jìn)行比較，來(lái)評(píng)價(jià)各類(lèi)重構(gòu)方法在環(huán)境星時(shí)序曲線重構(gòu)中的效果。在對(duì)植被生長(zhǎng)的物候信息提取中采用J?nsson和Eklundh改進(jìn)后的動(dòng)態(tài)閾值法對(duì)物候信息進(jìn)行提取，該方法在時(shí)間域和空間域上都具有很好的適用性[20]；其中，閾值設(shè)定為J?nsson在提出動(dòng)態(tài)閾值方法時(shí)建議的生長(zhǎng)季開(kāi)始和結(jié)束的NDVI值是年振幅的20%[20]。

對(duì)于農(nóng)作物(春玉米、冬小麥-夏玉米)，提取重構(gòu)后時(shí)序曲線的作物生長(zhǎng)信息，通過(guò)對(duì)比北京地區(qū)農(nóng)作物物候歷，對(duì)各類(lèi)重構(gòu)方法在環(huán)境星數(shù)據(jù)適用性進(jìn)行評(píng)價(jià)；其中，返青期的提取采用與林地與草地相同的動(dòng)態(tài)閾值法和動(dòng)態(tài)閾值20%，由于農(nóng)作物收割不屬于植被的自然生長(zhǎng)，且收割會(huì)使農(nóng)作物NDVI值在短時(shí)間內(nèi)大幅度下降，因此農(nóng)作物收割時(shí)間采用最大斜率法進(jìn)行提取。

最后，通過(guò)物候信息一致性的比較，以時(shí)間偏差對(duì)各重構(gòu)方法的重構(gòu)效果進(jìn)行定量評(píng)價(jià)。

3.3.2曲線相關(guān)性評(píng)價(jià)

在曲線重構(gòu)中，重構(gòu)方法對(duì)原始數(shù)據(jù)會(huì)產(chǎn)生擾動(dòng)，主要來(lái)自在重構(gòu)過(guò)程中數(shù)據(jù)有效性判斷、異常值剔除以及重構(gòu)函數(shù)擬合這3個(gè)方面，且這3種擾動(dòng)是不可避免的，所以在時(shí)序曲線重構(gòu)中盡可能地縮小擬合值與原始值的差異，以減小重構(gòu)方法對(duì)原始序列數(shù)據(jù)過(guò)度擾動(dòng)[22]。非植被地物的NDVI時(shí)序曲線沒(méi)有明顯的生長(zhǎng)峰值，因此，通過(guò)比較各重構(gòu)方法對(duì)原始時(shí)序數(shù)據(jù)的擾動(dòng)程度，來(lái)進(jìn)行重構(gòu)方法的評(píng)價(jià)。其中，擾動(dòng)程度由原始曲線與重構(gòu)曲線的相關(guān)統(tǒng)計(jì)量比較得出。

4實(shí)驗(yàn)結(jié)果與重構(gòu)效果評(píng)價(jià)

4.1重構(gòu)方法的實(shí)現(xiàn)

基于上述4種時(shí)間序列重構(gòu)方法，研究中采用TIMESAT與IDL實(shí)現(xiàn)時(shí)間序列的重構(gòu)。其中，A-G、D-L、S-G 3種重構(gòu)方法基于TIMESAT 3.1.1平臺(tái)，通過(guò)對(duì)比不同參數(shù)對(duì)重構(gòu)結(jié)果的影響，以及參考國(guó)內(nèi)外相關(guān)研究[16，20，23]。最終將參數(shù)設(shè)定如下：NDVI有效值域(Data Range)為-1～1，噪聲去除閾值(Spike Method)為2，擬合峰值參數(shù)(Adaptation Strength)為3，迭代次數(shù)(Envelope Iterations)為2，滑動(dòng)窗口大小(Window Size)為6或7；Hants方法的實(shí)現(xiàn)基于IDL-HANTS v1.3，其參數(shù)設(shè)置包括：最大下降偏移量(FET)為0.2，NDVI有效值域(Range)為-1～1，異常點(diǎn)最大數(shù)(TAT)為5，迭代次數(shù)(iMAX)為10，以及諧波頻數(shù)。其中，諧波頻數(shù)設(shè)置過(guò)小會(huì)導(dǎo)致擬合曲線過(guò)于平滑而缺失細(xì)節(jié)信息，此外諧波頻數(shù)與植物生長(zhǎng)季及參與擬合點(diǎn)個(gè)數(shù)相關(guān)，在本研究區(qū)內(nèi)植物生長(zhǎng)峰值數(shù)不超過(guò)3，且可滿足參與擬合點(diǎn)個(gè)數(shù)大于最大頻數(shù)的2倍減1。因此，在本研究中，諧波頻數(shù)最小設(shè)置為2，最大設(shè)置為3。

4.2重構(gòu)結(jié)果分析

本文為分析上述4種NDVI時(shí)間序列重構(gòu)方法對(duì)環(huán)境星數(shù)據(jù)的適用性，對(duì)6類(lèi)典型地物(林地、草地、春玉米、冬小麥-夏玉米、建筑及水體)的原始數(shù)據(jù)分別進(jìn)行了時(shí)間序列的重構(gòu)，重構(gòu)結(jié)果如圖3所示。

圖3　環(huán)境星數(shù)據(jù)NDVI時(shí)間序列重構(gòu)結(jié)果

圖3中，經(jīng)過(guò)4種重構(gòu)方法處理后的NDVI時(shí)序曲線較原始曲線去除了由噪聲帶來(lái)的影響，整條曲線更為平滑。其中，植被時(shí)序曲線經(jīng)過(guò)重構(gòu)后更符合植物的生長(zhǎng)特點(diǎn)，林地、草地、春玉米表現(xiàn)出一年單峰的生長(zhǎng)特點(diǎn)，且均在8月達(dá)到峰值，冬小麥-夏玉米表現(xiàn)出一年兩峰的生長(zhǎng)特點(diǎn)，其波峰與波谷的位置也與北京地區(qū)冬小麥-夏玉米的播種—生長(zhǎng)—收割時(shí)間較一致；非植被地物中，建筑經(jīng)重構(gòu)后的NDVI時(shí)序曲線變化不大，近似于某一常數(shù)，介于0～0.1之間，符合該類(lèi)別的光譜反射特征，水體經(jīng)重構(gòu)后的NDVI普遍在0值以下，而夏季部分高值的出現(xiàn)為該時(shí)段水體表面的植物生長(zhǎng)所致。因此，上述4種方法均可有效地對(duì)HJ數(shù)據(jù)進(jìn)行NDVI時(shí)間序列重構(gòu)。

4.3重構(gòu)效果評(píng)價(jià)

為進(jìn)一步區(qū)分各重構(gòu)方法的差異，本研究從以下兩方面對(duì)環(huán)境星數(shù)據(jù)的重構(gòu)效果進(jìn)行定量評(píng)價(jià)：①植被地物重構(gòu)后的環(huán)境星時(shí)序曲線對(duì)物候信息反映的一致性評(píng)價(jià)；②非植被地物重構(gòu)曲線與原始曲線的相關(guān)性評(píng)價(jià)。

4.3.1物候信息一致性評(píng)價(jià)結(jié)果

對(duì)林地、草地MODIS數(shù)據(jù)的重構(gòu)采用與環(huán)境星數(shù)據(jù)相同的實(shí)驗(yàn)方法及重構(gòu)參數(shù)。其中，在S-G濾波中，考慮MODIS數(shù)據(jù)的時(shí)間間隔為16天，而環(huán)境星數(shù)據(jù)為11天，因此滑動(dòng)窗口大小設(shè)置為4或5(環(huán)境星數(shù)據(jù)為6或7)。MODIS重構(gòu)結(jié)果如圖4所示。

對(duì)重構(gòu)后的環(huán)境星與MODIS時(shí)序曲線用同一動(dòng)態(tài)閾值20%，分別提取林地與草地的返青期與枯黃期，并得到不同方法下返青期與枯黃期的時(shí)間偏差，如圖5所示。從圖5中可見(jiàn)，在4種方法中，D-L對(duì)環(huán)境星時(shí)序數(shù)據(jù)的重構(gòu)效果最好，時(shí)間偏差均在5天以?xún)?nèi)，對(duì)林地及草地物候期的提取最為精確；S-G在兩類(lèi)植被地物的物候提取上偏差最大，其平均偏差在7天以上，其原因可能是環(huán)境星原始數(shù)據(jù)中的噪聲較多，滑動(dòng)窗口平均濾波的方法屬于局部擬合，不能完全消除噪聲的影響，而相比其他3種方法對(duì)數(shù)據(jù)整體進(jìn)行半局部和全局?jǐn)M合，能夠更好地描述時(shí)序曲線的變化趨勢(shì)和全局特征。A-G在林地返青期和枯黃期的提取精度均上要優(yōu)于Hants，且Hants在重構(gòu)中對(duì)原始數(shù)據(jù)在生長(zhǎng)季外出現(xiàn)的假高值不能很好地處理，如圖4(a)與圖4(f)所示。在本實(shí)驗(yàn)林地?cái)?shù)據(jù)中，出現(xiàn)在返青期前的假高值，Hants重構(gòu)后并未將其消除，而在原始數(shù)據(jù)中沒(méi)有假高值影響的草地類(lèi)型重構(gòu)中，利用Hants重構(gòu)后的環(huán)境星時(shí)序曲線對(duì)返青期的提取結(jié)果要更為精確。

圖4　MODIS數(shù)據(jù)的4種方法NDVI時(shí)間序列重構(gòu)結(jié)果

對(duì)于農(nóng)作物(春玉米、冬小麥-夏玉米)，提取重構(gòu)后時(shí)序曲線的作物生長(zhǎng)信息，對(duì)比北京地區(qū)農(nóng)作物物候歷，如表1所示。對(duì)圖3中4種方法的重構(gòu)曲線分別提取作物生長(zhǎng)信息，其中，對(duì)春玉米提取返青(出苗)期和收獲期；對(duì)冬小麥-夏玉米提取冬小麥的返青期、生長(zhǎng)峰值期(抽穗)、夏玉米的生長(zhǎng)峰值期(吐絲后期)和夏玉米的收獲期；并將所提取的物候時(shí)間與表1中所對(duì)應(yīng)時(shí)段的中間時(shí)刻進(jìn)行比較，計(jì)算時(shí)間偏差，其結(jié)果如圖6所示。

表1　北京地區(qū)農(nóng)作物物候歷

圖5　不同方法重構(gòu)后環(huán)境星與MODIS物候期差異

圖6　環(huán)境星農(nóng)作物在不同方法重構(gòu)后的差異比較

從圖6可以看出4種方法對(duì)環(huán)境星數(shù)據(jù)的重構(gòu)結(jié)果中，A-G對(duì)單季作物春玉米的物候提取，偏差最小；在雙季作物冬小麥-夏玉米物候的提取中，D-L重構(gòu)效果略?xún)?yōu)于A-G與Hants的提取結(jié)果。相比A-G與Hants兩種方法，在冬小麥返青期的提取中Hants偏差較小，在夏玉米生長(zhǎng)峰值提取中A-G更為準(zhǔn)確；綜合兩種作物的時(shí)間偏差，Hants的偏差最小，但由圖3(d)中可以看出，Hants重構(gòu)后的農(nóng)作物曲線過(guò)于平滑，且峰值兩邊嚴(yán)格對(duì)稱(chēng)，這對(duì)存在人為影響的農(nóng)作物，時(shí)序曲線中的非對(duì)稱(chēng)信息很難提取[24]。同上述林地與草地物候的提取，S-G對(duì)兩類(lèi)作物物候期的提取偏差最大，其原因?yàn)榄h(huán)境星數(shù)據(jù)噪聲的影響與局部擬合方法有一定局限性。

4.3.2曲線相關(guān)性評(píng)價(jià)結(jié)果

選取原始曲線與重構(gòu)曲線的最大值、最小值、均值、均方差及二者的相關(guān)系數(shù)進(jìn)行重構(gòu)方法的擾動(dòng)程度比較，其比較結(jié)果如表2所示。

表2　原始曲線與4種方法處理結(jié)果對(duì)照表

從表2可以看出，經(jīng)4種方法的重構(gòu)，兩種地物的原始曲線均達(dá)到了平滑的效果，其效果表現(xiàn)在曲線最大值的降低、最小值的升高與均方差的減小上；由于植被指數(shù)時(shí)間序列的噪聲主要來(lái)自云、大氣的影響，而造成值的下降，4種濾波方法中均有對(duì)此類(lèi)異常值的判斷與處理。可以看到在A-G、D-L和S-G 3種方法重構(gòu)后，曲線均值相比原始曲線均有提高，而Hants重構(gòu)曲線的均值保持不變，其原因是Hants方法是以傅里葉變換為原理的全局?jǐn)M合。從數(shù)據(jù)整體減小了重構(gòu)過(guò)程對(duì)原始數(shù)據(jù)的擾動(dòng)，這一結(jié)論與表中重構(gòu)曲線與原始曲線的相關(guān)系數(shù)一項(xiàng)的結(jié)果一致。綜合上述結(jié)果，對(duì)于非植被地物，4種方法中Hants在重構(gòu)過(guò)程中對(duì)原始數(shù)據(jù)的擾動(dòng)最小，S-G其次，A-G與D-L由于為特定函數(shù)的半局部擬合，對(duì)原始數(shù)據(jù)的擾動(dòng)均較大。

5結(jié)束語(yǔ)

本文以北京為研究區(qū)，比較分析了非對(duì)稱(chēng)高斯函數(shù)擬合、Double Logistic函數(shù)擬合、S-G濾波和諧波分析法這4種時(shí)間序列重構(gòu)方法在環(huán)境星數(shù)據(jù)中的重構(gòu)效果。通過(guò)植被物候信息比較和重構(gòu)過(guò)程對(duì)原始數(shù)據(jù)的擾動(dòng)分析，定量評(píng)價(jià)了這4種方法在植被地物的時(shí)序曲線重構(gòu)中的效果，總體來(lái)看4種方法重構(gòu)后的曲線都達(dá)到了對(duì)噪聲的去除與對(duì)整體曲線的平滑。在植被物候信息的準(zhǔn)確性上，對(duì)于林地與草地，Double Logistic函數(shù)擬合在環(huán)境星數(shù)據(jù)與高質(zhì)量MODIS數(shù)據(jù)的重構(gòu)中，植被物候信息的一致性最高，非對(duì)稱(chēng)高斯函數(shù)擬合與Hants效果次之。其中，Hants對(duì)非生長(zhǎng)季內(nèi)的假高值不能很好的處理，S-G偏差較大；在農(nóng)作物時(shí)序曲線的重構(gòu)中，對(duì)比農(nóng)作物物候歷，S-G在兩類(lèi)作物(單季與雙季)中偏差均較大，A-G，D-L，Hants的精確性較高，但相比A-G與D-L的重構(gòu)結(jié)果，Hants重構(gòu)曲線無(wú)法反映收割對(duì)NDVI值造成的突降。對(duì)于水體與建筑這兩類(lèi)非植被地物，分析4種重構(gòu)方法對(duì)原始數(shù)據(jù)的擾動(dòng)程度，發(fā)現(xiàn)Hants重構(gòu)后的曲線與原始曲線最為接近，相關(guān)系數(shù)較高；A-G，D-L對(duì)原始曲線的改動(dòng)較大。

綜上，對(duì)于環(huán)境星數(shù)據(jù)，A-G與D-L更適用于重構(gòu)植被地物的時(shí)序曲線，可以很好地用于植被物候提取、植被變化監(jiān)測(cè)以及農(nóng)作物分類(lèi)等研究；對(duì)非植被地物，Hants對(duì)原始時(shí)序數(shù)據(jù)的擾動(dòng)最小，重構(gòu)曲線更接近原始時(shí)間序列。

參考文獻(xiàn)：

[1]趙英時(shí).遙感應(yīng)用分析原理與方法[M].北京：科學(xué)出版社，2003.

[2]LYON J G，YUAN D，LUNETTA R S，et al.A change detection experiment using vegetation indices[J].Photogrammetric Engineering and Remote Sensing，1998，(64)：143-150.

[3]PETTORELLI N，VIK J O，MYSTRUD A，et al.Using the satellite-derived ndvi to assess ecological responses to environmental change[J].Trends in Ecology & Evolution，2005，(20)：503-510.

[4]JONSSON P，EKLUNDH L.Seasonality extraction by function fitting to time-series of satellite sensor data[J].Geoscience and Remote Sensing，IEEE Transactions on，2002，(40)：1824-1832.

[5]REED B C，BROWN J F，VANDER Z D，et al.Measuring phenological variability from satellite imagery[J].Journal of Vegetation Science，1994，(5)：703-714.

[6]DEFRIES R，TOWNSHEND J.Ndvi-derived land cover classifications at a global scale[J].International Journal of Remote Sensing，1994，(15)：3567-3586.

[7]LUNETTA R S，KNIGHT J F，EDIRIWICKREMA J，et al.Land-cover change detection using multi-temporal modisndvi data[J].Remote sensing of environment，2006，(105)：142-154.

[8]李儒，張霞，劉波，等.遙感時(shí)間序列數(shù)據(jù)濾波重建算法發(fā)展綜述[J].遙感學(xué)報(bào)，2009，13(2)：335-341.

[9]MA M，VEROUSTRAETE F.Reconstructing pathfinder avhrr land ndvi time-series data for the northwest of China[J].Advances in Space Research，2006，(37：835-840.

[10]BRADLEY B A，JACOB R W，HERMANCE J F，et al.A curve fitting procedure to derive inter-annual phenologies from time series of noisy satellite ndvidata[J].Remote Sensing of Environment，2007，(106)：137-145.

[11]林忠輝，莫興國(guó).NDVI時(shí)間序列諧波分析與地表物候信息獲取[J].農(nóng)業(yè)工程學(xué)報(bào)，2006，22(12)：138-144.

[12]ZHU W，PAN Y，HE H，et al.A changing-weight filter method for reconstructing a high-quality ndvi time series to preserve the integrity of vegetation phenology[J].Geoscience and Remote Sensing，IEEE Transactions on，2012，(50)：1085-1094.

[13]吳文斌，楊鵬，唐華俊，等.兩種NDVI時(shí)間序列數(shù)據(jù)擬合方法比較[J].農(nóng)業(yè)工程學(xué)報(bào)，2009，25(11)：183-188.

[14]BECK P S，ATZBERGER C，H?GDA K A，et al.Improved monitoring of vegetation dynamics at very high latitudes：A new method using modisndvi[J].Remote sensing of Environment，2006，(100)：321-334.

[15]HIRD J N，MCDERMID G J.Noise reduction of ndvi time series：An empirical comparison of selected techniques[J].Remote Sensing of Environment，2009，(113)：248-258.

[16]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-golayFilter[J].Remote sensing of Environment，2004，(91)：332-344.

[17]JULIEN Y，SOBRINO J A.Comparison of cloud-reconstruction methods for time series of composite ndvidata[J].Remote Sensing of Environment，2010，(114)：618-625.

[18]陳雪洋，蒙繼華，杜鑫，等.基于環(huán)境星CCD數(shù)據(jù)的冬小麥葉面積指數(shù)遙感監(jiān)測(cè)模型研究[J].國(guó)土資源遙感，2010，22(2)：55-58.

[19]孫志群，劉志輝，HJ-1 CCD影像大氣校正模型比較研究[J].International Conference on Remote Sensing，2010，(3).

[20]J?NSSON P，EKLUNDH L.Timesat—a program for analyzing time-series of satellite sensor data[J].Computers & Geosciences，2004，(30)：833-845.

[21]ROERINK G，MENENTI M，VERHOEF W.Reconstructing cloudfreendvi composites using fourier analysis of time series[J].International Journal of Remote Sensing，2000，(21)：1911-1917.

[22]李儒.遙感植被指數(shù)時(shí)間序列數(shù)據(jù)濾波重建算法研究[D].北京：中國(guó)科學(xué)院遙感應(yīng)用研究所，2008.

[23]宋春橋，柯靈紅，游松財(cái)，等.基于Timesat 的3種時(shí)序NDVI擬合方法比較研究—以藏北草地為例[J].遙感技術(shù)與應(yīng)用，2011，26(2)：147-155.

[24]OLSSON L，EKLUNDH L.Fourier series for analysis of temporal sequences of satellite sensor imagery[J].International Journal of Remote Sensing，1994，(15)：3735-3741.

E-mail：yjs010@163.com