3 種NDVI 時序數據重建方法對比研究

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

（1. 生態環境部衛星環境應用中心，北京 100094）

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

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

1 研究區概況與數據來源

1.1 研究區概況

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

表1 研究區概況

1.2 數據來源

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

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

3）土地利用數據。GLC_2010土地覆蓋數據集的空間分辨率為30 m，分類十分詳細，總體精度也達到了80%以上，因此選擇該數據集作為本文分類結果驗證時的真實數據[17]。GLC_2010土地覆蓋數據集的獲取網址為http://www.globallandcover.com/GLC30Download/index.aspx。

2 研究方法

2.1 時序數據重建方法

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

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

整體擬合函數的表達式為：

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

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

其局部擬合函數的表達式為：

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

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

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

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

2.2 多維度時序數據重建結果評價

本文在分析3種NDVI時序數據重建方法原理的基礎上，從不同維度出發，利用3種方法重建NDVI時序數據，并選取多個評價指標對比分析了3種方法的優劣和適用性，為今后NDVI時序數據重建方法的選擇提供參考和依據。第一個維度，對比分析3種方法對3個地區Landsat-NDVI 時序數據重建的效果；第二個維度，對比分析3 種方法對同一地區不同覆被類型Land?sat-NDVI時序數據重建的效果。技術路線如圖1所示。

圖1 技術路線圖

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

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

圖2 兗州市土地利用圖

2.3 基于統計量的時序數據重建結果評價

本文選取相關系數和回歸估計標準差（RMSE）定量分析3種方法的重建效果。

相關系數表征的是兩個樣本數組值之間的相關度。本文對比了3 種方法擬合前后的相關系數，進而評價3種方法擬合后NDVI 時序曲線保持原始NDVI曲線特征的能力。相關系數越大，表示相關性越高，保真性越好。其計算公式為：

本文計算得到每個像元的相關系數，并統計所有像元的平均值。

RMSE 表示兩個樣本數組值之間的差異程度。本文用于對比重建前后NDVI 值之間的平均差異程度，其值越小，擬合值的代表性越強。其計算公式為：

3 研究結果與分析

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

3.1 不同研究區3種方法重建結果評價

3.1.1 像元水平重建結果對比

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

圖3 3個研究區典型像元的NDVI時序曲線

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

綜上所述，從單個像元時間序列曲線對比來看，AG 方法和D-L 方法能反映出NDVI 時序曲線的整體變化特征，但D-L 方法出現了過度擬合的現象；S-G 方法能較好地反映曲線細節特征，但曲線整體不夠平滑，難以反映整體時序變化特征。

3.1.2 基于統計量的重建結果對比

本文計算得到3 種方法重建后的相關系數和RMSE，如圖4 所示，可以看出，AG 方法和D-L 方法的相關系數遠高于S-G方法，AG方法的相關系數略高于D-L 方法；AG 方法和D-L 方法的RMSE 小于S-G 方法，AG方法略小于D-L方法。

圖4 3個地區3種方法重建后相關系數和RMSE

3.2 不同植被類型3種方法的重建效果評價

3.2.1 像元水平重建效果對比

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

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

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

圖5 耕地和林地典型像元NDVI時序曲線

3.2.2 重建質量評價

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

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

4 結 語

本文以弗吉尼亞州韋茲縣、陜西省神府地區、山東省兗州市為研究區，利用AG方法、D-L方法和S-G方法重建了NDVI 時間序列數據；并從像元尺度、影像尺度、保真性和統計值等方面選取多個評價指標，對比分析了不同區域和不同植被類型3 種方法的優缺點和適用性，可為相關研究中NDVI 時序數據重建方法的選擇提供參考和依據。

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

2）NDVI時序數據重建結果受植被覆蓋類型的影響較大，不同植被覆蓋類型適用的重建方法不同。3種方法耕地的重建效果均較好且相差較小，AG方法和D-L 方法對于林地的重建質量高于S-G 方法。因此，重建耕地NDVI時序數據可任選3種方法之一，重建林地NDVI時序數據應選用AG方法或D-L方法。