邰文飛，申文明，蔡明勇，申 振，張新勝

（1. 生態(tài)環(huán)境部衛(wèi)星環(huán)境應(yīng)用中心，北京 100094）

地表植被是陸地自然生態(tài)系統(tǒng)的主要構(gòu)成部分，是生態(tài)環(huán)境變化的敏感指標(biāo)[1]。研究全球或區(qū)域一段時(shí)間內(nèi)的植被覆蓋變化，是氣候與環(huán)境變化監(jiān)測(cè)和環(huán)境質(zhì)量評(píng)價(jià)的主要手段。歸一化植被指數(shù)（NDVI）與植被分布密度呈線性相關(guān)[2]，是目前衛(wèi)星遙感監(jiān)測(cè)植被特征應(yīng)用最廣泛的指數(shù)之一[3-4]。受地面狀況、大氣干擾、傳感器狀態(tài)等因素影響，NDVI 時(shí)間序列數(shù)據(jù)集包含大量噪聲，需要根據(jù)不同研究區(qū)域和植被覆蓋類(lèi)型，選取適用的濾波算法或擬合方法[5]。常用方法主要包括非線性擬合、濾波平滑、閾值去除3類(lèi)。TIME?SAT 軟件集成了非對(duì)稱的高斯函數(shù)擬合法（AG）、雙Logistic 曲線擬合法（D-L） 和Savitzky-Golay 濾波（S-G）3種方法，可針對(duì)不同類(lèi)型的時(shí)序數(shù)據(jù)和研究目標(biāo)設(shè)置最優(yōu)的算法參數(shù)，獲得最好的重建效果。

針對(duì)遙感時(shí)間序列數(shù)據(jù)集重建算法，國(guó)內(nèi)外學(xué)者開(kāi)展了大量研究，如曹云鋒[6]等以遼寧省長(zhǎng)白山自然保護(hù)區(qū)為研究區(qū)，利用AG 方法、D-L 方法和S-G 方法重建了MOD13Q1中的NDVI時(shí)序數(shù)據(jù)集，并對(duì)比了3種方法對(duì)原始高質(zhì)量數(shù)據(jù)的保真性，結(jié)果表明AG方法擬合效果最優(yōu)，S-G方法最差；范德芹[7]等以洞庭湖區(qū)為研究區(qū)，首先利用狄克松檢驗(yàn)法對(duì)MODIS NDVI時(shí)序數(shù)據(jù)進(jìn)行噪聲檢測(cè)，再利用S-G 方法和變權(quán)重濾波法重建噪聲監(jiān)測(cè)后的數(shù)據(jù)，然后結(jié)合植被類(lèi)型圖選取520 個(gè)樣本點(diǎn)，分析和評(píng)價(jià)重建結(jié)果，結(jié)果表明時(shí)序數(shù)據(jù)重建前利用狄克松檢驗(yàn)法進(jìn)行監(jiān)測(cè)預(yù)處理能有效提高NDVI時(shí)序數(shù)據(jù)的重建質(zhì)量。目前，還沒(méi)有一種被大多數(shù)專(zhuān)家和學(xué)者認(rèn)可的實(shí)際應(yīng)用效果較好的普適性重建方法。由于各種重建方法的機(jī)理和算法存在很大差異，重建效果因地區(qū)分布、植被覆蓋特點(diǎn)和其研究?jī)?nèi)容而各異[8-11]。大量研究結(jié)果表明，AG方法、D-L方法和S-G 方法的重建效果較好，被廣泛應(yīng)用于遙感時(shí)間序列數(shù)據(jù)重建[12-16]。本文以弗吉尼亞州韋茲縣、陜西省神木縣和府谷縣、山東省兗州市為研究區(qū)，利用AG方法、D-L方法和S-G方法重建了Landsat TM和MOD13Q1的NDVI時(shí)間序列數(shù)據(jù)集，從像元尺度、影像尺度、保真性和統(tǒng)計(jì)值等方面選取了多個(gè)評(píng)價(jià)指標(biāo)，對(duì)比研究了3 種方法針對(duì)不同區(qū)域和不同植被類(lèi)型的優(yōu)缺點(diǎn)和適用性，為相關(guān)研究中NDVI 時(shí)序數(shù)據(jù)重建方法使用和選擇提供了參考和依據(jù)。

1 研究區(qū)概況與數(shù)據(jù)來(lái)源

1.1 研究區(qū)概況

韋茲縣屬于阿巴拉契亞地區(qū)，地處西南弗吉尼亞煤田的中西部，地表大部分被森林覆蓋，NDVI 值總體較高，該地區(qū)露天開(kāi)采剝離了大量地表植被，對(duì)生態(tài)環(huán)境破壞巨大；神木縣和府谷縣（以下簡(jiǎn)稱神府地區(qū)）隸屬于陜西省榆林市，地處丘陵、森林草原向沙漠、干草原的過(guò)渡地帶，近年來(lái)由于煤礦開(kāi)采，植被變化劇烈且明顯，NDVI 時(shí)序曲線有較大波動(dòng)；兗州市地形以平原為主，耕地面積大、產(chǎn)量高，近年來(lái)森林覆蓋率持續(xù)增高，作為該市最具優(yōu)勢(shì)的礦產(chǎn)資源，煤炭分布面積約占其總面積的37.06%，平原面積占總面積的99%以上。3個(gè)地區(qū)幾乎處于同一緯度（表1），煤炭資源均較豐富，采礦等人類(lèi)活動(dòng)頻繁劇烈，地表植被破壞較嚴(yán)重，NDVI 時(shí)序數(shù)據(jù)變化受人類(lèi)活動(dòng)影響可能呈現(xiàn)一定的規(guī)律或趨勢(shì)；同時(shí)，由于3 個(gè)地區(qū)在氣候、海拔、地形和主要植被等方面存在差異，以這3個(gè)地區(qū)為研究區(qū)有利于評(píng)價(jià)NDVI時(shí)序數(shù)據(jù)重建方法對(duì)不同類(lèi)型研究區(qū)和不同植被覆蓋類(lèi)型的重建效果和適用性。

表1 研究區(qū)概況

1.2 數(shù)據(jù)來(lái)源

1）MODIS 時(shí)序數(shù)據(jù)集。本文采用2014 年山東省兗州市MODIS（MOD13Q1）數(shù)據(jù)產(chǎn)品，利用16 d 最大值合成（MVC），空間分辨率為250 m；利用MRT工具對(duì)影像進(jìn)行空間拼接、投影轉(zhuǎn)換和重采樣等處理，共構(gòu)建23期MODIS-NDVI時(shí)間序列數(shù)據(jù)。

2）Landsat TM時(shí)序數(shù)據(jù)集。本文選取弗吉尼亞州韋茲縣、陜西省神府地區(qū)、山東省兗州市1996—2011年的Landsat TM影像，韋茲縣為13幅影像，神府地區(qū)和兗州市均為15 幅影像，數(shù)據(jù)產(chǎn)品級(jí)別為1 T，空間分辨率為30 m。本文最大程度地獲取了時(shí)相差較小的Landsat TM 數(shù)據(jù)，受數(shù)據(jù)質(zhì)量限制，有幾期數(shù)據(jù)時(shí)間相差較大，但大部分?jǐn)?shù)據(jù)獲取時(shí)間集中在6—8月的植被生長(zhǎng)季。

3）土地利用數(shù)據(jù)。GLC_2010土地覆蓋數(shù)據(jù)集的空間分辨率為30 m，分類(lèi)十分詳細(xì)，總體精度也達(dá)到了80%以上，因此選擇該數(shù)據(jù)集作為本文分類(lèi)結(jié)果驗(yàn)證時(shí)的真實(shí)數(shù)據(jù)[17]。GLC_2010土地覆蓋數(shù)據(jù)集的獲取網(wǎng)址為http://www.globallandcover.com/GLC30Download/index.aspx。

2 研究方法

2.1 時(shí)序數(shù)據(jù)重建方法

AG 方法由PerJo?nsson和Lars Eklundh 在2002 年首次提出，是一種基于不對(duì)稱高斯函數(shù)的非線性最小二乘擬合算法。該方法采用分段擬合的思想，可以確保擬合結(jié)果更加接近符合當(dāng)前時(shí)段的真實(shí)變化規(guī)律，可很好地描述作物生長(zhǎng)和交替過(guò)程中植被指數(shù)與時(shí)間的關(guān)系[18]。李儒[19]等將其擬合過(guò)程分解為區(qū)間提取、局部擬合和整體擬合3個(gè)步驟。

局部擬合函數(shù)的表達(dá)式為：式中，c1、c2為曲線的基準(zhǔn)和振幅；t為時(shí)相；a1為曲線最大值或最小值的位置；a2、a3、a4、a5分別為左右半邊曲線的寬度和峭度。

整體擬合函數(shù)的表達(dá)式為：

式中，[tL,tR] 為時(shí)間序列中尚未擬合部分的NDVI 區(qū)間；fL(t)、fc(t)分別為[tL,tR] 區(qū)間內(nèi)左邊最小值、中間最大值和右邊最小值對(duì)應(yīng)的局部擬合函數(shù)；α(t)、β(t)為介于0和1之間的剪切系數(shù)。

D-L方法的原理與AG方法基本一致，同樣是先進(jìn)行局部擬合再進(jìn)行整體擬合。Beck認(rèn)為D-L方法對(duì)植被生長(zhǎng)季提取的準(zhǔn)確率更高[20]。與AG方法相比，D-L方法的局部擬合函數(shù)為雙邏輯形式，且公式中少一個(gè)參數(shù)[21]。

其局部擬合函數(shù)的表達(dá)式為：

式中，a2、a3、a4、a5分別為曲線左、右部分的拐點(diǎn)位置以及拐點(diǎn)處的變化速率。

其整體擬合函數(shù)參照AG 方法的整體擬合函數(shù)公式。

S-G 方法利用最小二乘卷積算法，通過(guò)計(jì)算一組相鄰值達(dá)到數(shù)據(jù)濾波的目的[22]。其公式為：

式中，Y為濾波前的NDVI 值；Y*為濾波后NDVI值；Ci為該像元第i期NDVI值濾波時(shí)的系數(shù)；N為迭代次數(shù)，是滑動(dòng)窗口的寬度（2m+1）；j為原始NDVI數(shù)組的系數(shù)。

2.2 多維度時(shí)序數(shù)據(jù)重建結(jié)果評(píng)價(jià)

本文在分析3種NDVI時(shí)序數(shù)據(jù)重建方法原理的基礎(chǔ)上，從不同維度出發(fā)，利用3種方法重建NDVI時(shí)序數(shù)據(jù)，并選取多個(gè)評(píng)價(jià)指標(biāo)對(duì)比分析了3種方法的優(yōu)劣和適用性，為今后NDVI時(shí)序數(shù)據(jù)重建方法的選擇提供參考和依據(jù)。第一個(gè)維度，對(duì)比分析3種方法對(duì)3個(gè)地區(qū)Landsat-NDVI 時(shí)序數(shù)據(jù)重建的效果；第二個(gè)維度，對(duì)比分析3 種方法對(duì)同一地區(qū)不同覆被類(lèi)型Land?sat-NDVI時(shí)序數(shù)據(jù)重建的效果。技術(shù)路線如圖1所示。

圖1 技術(shù)路線圖

1）以弗吉尼亞州韋茲縣、陜西省神府地區(qū)、山東省兗州市為研究區(qū)，這3 個(gè)地區(qū)幾乎處于同一緯度帶，但氣候、海拔和地形等因素相差較大，主要植被覆蓋類(lèi)型分別為林地、草地和耕地，這些差異導(dǎo)致3個(gè)地區(qū)的NDVI時(shí)間序列數(shù)據(jù)具有不同特征。本文選用Landsat TM 影像構(gòu)建時(shí)間序列影像，采用3 種方法重建NDVI 時(shí)序數(shù)據(jù)，輸出重建后的影像，并對(duì)比分析3種方法的重建效果。

2）以兗州市2014年23期16 d合成MODIS影像為數(shù)據(jù)源，采用GLC_2010 土地覆蓋數(shù)據(jù)集，將地表覆蓋類(lèi)型分為森林、耕地、草地、濕地、裸地等10 類(lèi)。由兗州市土地利用圖（圖2）可知，兗州市境內(nèi)的主要土地利用類(lèi)型為耕地，林地主要分布在兗州市北部和中部，草地分布極少。因此后續(xù)研究中剔除草地，主要研究耕地和林地。本文利用TIMESAT軟件分別重建耕地和林地的時(shí)序數(shù)據(jù)，并對(duì)比分析3 種方法對(duì)于不同類(lèi)型植被的重建效果。

圖2 兗州市土地利用圖

2.3 基于統(tǒng)計(jì)量的時(shí)序數(shù)據(jù)重建結(jié)果評(píng)價(jià)

本文選取相關(guān)系數(shù)和回歸估計(jì)標(biāo)準(zhǔn)差（RMSE）定量分析3種方法的重建效果。

相關(guān)系數(shù)表征的是兩個(gè)樣本數(shù)組值之間的相關(guān)度。本文對(duì)比了3 種方法擬合前后的相關(guān)系數(shù)，進(jìn)而評(píng)價(jià)3種方法擬合后NDVI 時(shí)序曲線保持原始NDVI曲線特征的能力。相關(guān)系數(shù)越大，表示相關(guān)性越高，保真性越好。其計(jì)算公式為：

本文計(jì)算得到每個(gè)像元的相關(guān)系數(shù)，并統(tǒng)計(jì)所有像元的平均值。

RMSE 表示兩個(gè)樣本數(shù)組值之間的差異程度。本文用于對(duì)比重建前后NDVI 值之間的平均差異程度，其值越小，擬合值的代表性越強(qiáng)。其計(jì)算公式為：

3 研究結(jié)果與分析

式中，NDVIoi、NDVIpi分別為時(shí)間序列中第i期擬合處理前后的NDVI 值；N為像元總數(shù)。

3.1 不同研究區(qū)3種方法重建結(jié)果評(píng)價(jià)

3.1.1 像元水平重建結(jié)果對(duì)比

本文選取具有高度代表性的像元對(duì)重建效果進(jìn)行定性分析，結(jié)果如圖3 所示，可以看出，韋茲縣境內(nèi)大部分被森林覆蓋，重建前整個(gè)時(shí)間序列的NDVI 值在0.8左右波動(dòng)，2000年NDVI值突降至0.5，2007年突降為0，2007年明顯是噪點(diǎn)，在重建時(shí)應(yīng)去除。3種方法重建后，NDVI 時(shí)序曲線更加平滑，2000 年的NDVI值出現(xiàn)不同程度的提升，3種方法重建時(shí)均可去除較明顯的噪點(diǎn)（2007 年）；相較于S-G 方法，AG 方法和D-L 方法重建后的曲線更加光滑，但D-L 方法重建后幾乎變?yōu)橹本€，與原始曲線偏離較大，出現(xiàn)了過(guò)度擬合的現(xiàn)象，這一方面是與其算法有關(guān)，另一方面是由原始時(shí)序曲線的變化幅度較小所致。

圖3 3個(gè)研究區(qū)典型像元的NDVI時(shí)序曲線

相較于韋茲縣，神府地區(qū)原始NDVI 時(shí)序曲線上下波動(dòng)較大，但整體變化趨勢(shì)明顯，大體呈先升后降的變化特征。3 種方法重建后的曲線更加平滑，顯示出明顯的周期，2001 年和2009 年突降點(diǎn)在重建時(shí)被去除。AG方法和D-L方法重建結(jié)果較好地反映了曲線整體周期性變化特征，而S-G 方法對(duì)曲線細(xì)節(jié)的擬合較好，更加接近原始曲線包絡(luò)線。AG方法和D-L方法的基本思想是利用分段高斯函數(shù)來(lái)模擬植被生長(zhǎng)過(guò)程，因此這兩種方法對(duì)于具有明顯生長(zhǎng)周期的NDVI時(shí)序曲線擬合效果更好。兗州市的主要植被覆蓋類(lèi)型為耕地，NDVI 值整體較高，主要分布在0.6～0.8 之間，但2003 年、2004 年出現(xiàn)了連續(xù)異常低值點(diǎn)，這可能與影像獲取時(shí)間有關(guān)。整體上，3 種方法重建結(jié)果有不同程度的“抬高”， NDVI 變化曲線無(wú)明顯規(guī)律，3種方法重建效果差別較大。

綜上所述，從單個(gè)像元時(shí)間序列曲線對(duì)比來(lái)看，AG 方法和D-L 方法能反映出NDVI 時(shí)序曲線的整體變化特征，但D-L 方法出現(xiàn)了過(guò)度擬合的現(xiàn)象；S-G 方法能較好地反映曲線細(xì)節(jié)特征，但曲線整體不夠平滑，難以反映整體時(shí)序變化特征。

3.1.2 基于統(tǒng)計(jì)量的重建結(jié)果對(duì)比

本文計(jì)算得到3 種方法重建后的相關(guān)系數(shù)和RMSE，如圖4 所示，可以看出，AG 方法和D-L 方法的相關(guān)系數(shù)遠(yuǎn)高于S-G方法，AG方法的相關(guān)系數(shù)略高于D-L 方法；AG 方法和D-L 方法的RMSE 小于S-G 方法，AG方法略小于D-L方法。

圖4 3個(gè)地區(qū)3種方法重建后相關(guān)系數(shù)和RMSE

3.2 不同植被類(lèi)型3種方法的重建效果評(píng)價(jià)

3.2.1 像元水平重建效果對(duì)比

通過(guò)3 種方法重建兗州市2014 年23 期MODIS 影像后，利用土地利用圖，結(jié)合Google Earth 中兗州市2014年前后影像，選取各地類(lèi)的純凈樣本像元；在分析大量耕地和林地像元NDVI 時(shí)間序列曲線噪聲特征以及重建效果的基礎(chǔ)上，選取既能反映整體季節(jié)變化特征又包含明顯噪聲的典型像元，分析重建效果。

兗州市主要農(nóng)作物是冬小麥和夏玉米。該地區(qū)氣候適宜，一年兩季，一季小麥、一季玉米，其時(shí)序曲線（圖5）具有明顯的“雙峰”特征，對(duì)應(yīng)一年中兩個(gè)完整的生長(zhǎng)周期。原始曲線在第一個(gè)生長(zhǎng)季開(kāi)始時(shí)，上下起伏較大，包含較多噪聲，第二個(gè)生長(zhǎng)周期曲線比較平滑。3 種方法重建后時(shí)序曲線幾乎重合，對(duì)于包含較多噪聲的第一個(gè)生長(zhǎng)季，3 種方法重建效果較好，去除了大部分噪聲，基本準(zhǔn)確還原了冬小麥生長(zhǎng)季NDVI 時(shí)序變化特征。對(duì)于夏玉米生長(zhǎng)季，原始曲線本身噪聲較少，3種方法較好地保持了原始曲線特征，但重建后生長(zhǎng)季整體向后“偏移”?？傮w而言，3種方法對(duì)耕地像元時(shí)序重建效果比較理想，重建結(jié)果較一致，因此對(duì)于耕地這3種方法重建差別較小。

兗州市境內(nèi)林地分布較少，由于空間分辨率的限制，純凈像元較少，大部分樣本像元受到周?chē)氐母蓴_，因此時(shí)序曲線無(wú)法較好地反映林地真實(shí)的時(shí)間變化特征。如圖5 所示，原始曲線整體呈現(xiàn)一個(gè)完整的生長(zhǎng)周期，符合林地一年一個(gè)生長(zhǎng)季的特點(diǎn)，但受到耕地干擾，第161、177天NDVI值出現(xiàn)“突降”，這應(yīng)該屬于噪聲；重建之后，第161、177 天NDVI 值得到提升，整個(gè)曲線比原始曲線更加平滑。總體而言，相較于耕地，3 種方法對(duì)林地的重建效果差別較大，僅從像元水平對(duì)比發(fā)現(xiàn)，AG方法和S-G方法對(duì)曲線末端兩個(gè)非噪聲點(diǎn)的擬合效果較差，曲線出現(xiàn)明顯的“上揚(yáng)”，而D-L方法不僅很好地反映了曲線整體變化特征，而且較好地反映了曲線的細(xì)節(jié)特征。

圖5 耕地和林地典型像元NDVI時(shí)序曲線

3.2.2 重建質(zhì)量評(píng)價(jià)

本文計(jì)算得到耕地和林地樣本點(diǎn)的相關(guān)系數(shù)和RMSE，進(jìn)而對(duì)比分析3種方法的重建質(zhì)量。耕地重建后，AG方法、D-L方法和S-G方法的相關(guān)系數(shù)分別為0.91、0.92、0.90；林地重建后，AG 方法、D-L 方法和S-G方法的相關(guān)系數(shù)分別為0.94、0.96、0.70。耕地的相關(guān)系數(shù)較高且3種方法相差較小，而林地AG方法和D-L 方法的相關(guān)系數(shù)遠(yuǎn)高于S-G 方法。因此，對(duì)于耕地而言，3 種方法的重建質(zhì)量較高且相差較??；對(duì)于林地而言，AG 方法和D-L 方法的重建質(zhì)量較高，S-G 方法的重建質(zhì)量較差；無(wú)論是耕地還是林地，重建質(zhì)量最高的是D-L方法，其次是AG方法，S-G方法最低。RMSE 可反映保真性高低，其值越小，保真性越高。耕地和林地的RMSE 大小均為D-L 方法＜AG 方法＜S-G 方法，即D-L 方法保真性最高，AG 方法次之，S-G方法保真性最差。

綜上所述，不同植被覆蓋類(lèi)型適用的重建方法不同，3 種方法耕地的重建效果相差較小，而林地D-L方法和AG方法的重建質(zhì)量較好，S-G方法較差。

4 結(jié) 語(yǔ)

本文以弗吉尼亞州韋茲縣、陜西省神府地區(qū)、山東省兗州市為研究區(qū)，利用AG方法、D-L方法和S-G方法重建了NDVI 時(shí)間序列數(shù)據(jù)；并從像元尺度、影像尺度、保真性和統(tǒng)計(jì)值等方面選取多個(gè)評(píng)價(jià)指標(biāo)，對(duì)比分析了不同區(qū)域和不同植被類(lèi)型3 種方法的優(yōu)缺點(diǎn)和適用性，可為相關(guān)研究中NDVI 時(shí)序數(shù)據(jù)重建方法的選擇提供參考和依據(jù)。

1）氣候、海拔、地形和主要植被類(lèi)型等因素將影響數(shù)據(jù)重建結(jié)果。3 種方法均有一定的去噪能力。從像元尺度、影像尺度、保真性和統(tǒng)計(jì)值等方面進(jìn)行對(duì)比發(fā)現(xiàn)，S-G 方法屬于動(dòng)態(tài)濾波，重建過(guò)程中更關(guān)注“細(xì)節(jié)”，但保留了大量噪聲；AG 方法和D-L 方法更加注重整體效果，重建后時(shí)序曲線較平滑，符合植被生長(zhǎng)變化的連續(xù)性規(guī)律，但會(huì)忽略時(shí)序曲線的真實(shí)局部變化特征。

2）NDVI時(shí)序數(shù)據(jù)重建結(jié)果受植被覆蓋類(lèi)型的影響較大，不同植被覆蓋類(lèi)型適用的重建方法不同。3種方法耕地的重建效果均較好且相差較小，AG方法和D-L 方法對(duì)于林地的重建質(zhì)量高于S-G 方法。因此，重建耕地NDVI時(shí)序數(shù)據(jù)可任選3種方法之一，重建林地NDVI時(shí)序數(shù)據(jù)應(yīng)選用AG方法或D-L方法。