我國西北地區東部時間序列NDVI數據集重建方法比較研究

王瑋，郭鈮，沙莎，胡蝶，王小平，李耀輝

(中國氣象局蘭州干旱氣象研究所，甘肅省干旱氣候變化與減災重點實驗室，中國氣象局干旱氣候變化與減災重點實驗室，甘肅 蘭州 730020)

?

我國西北地區東部時間序列NDVI數據集重建方法比較研究

王瑋，郭鈮*，沙莎，胡蝶，王小平，李耀輝

(中國氣象局蘭州干旱氣象研究所，甘肅省干旱氣候變化與減災重點實驗室，中國氣象局干旱氣候變化與減災重點實驗室，甘肅 蘭州 730020)

高質量、長時序歸一化植被指數(NDVI)數據集不僅是連續監測陸地表面特征的基礎，也是研究氣候與陸地生態系統變化的重要參數。本研究以生態環境較為脆弱的西北地區東部為例，借助多種時間序列重建方法對LTDR NDVI數據集中的噪聲進行擬合重建，并結合農業氣象資料和高質量NDVI數據，對不同重建方法的擬合結果開展適用性評價分析，結果表明，1)下墊面類型是影響重建方法擬合效果的重要因素。根據不同植被類型或作物生長特點，每種重建方法對其噪聲消除能力有所不同；2)在年均NDVI較高(NDVI≥0.3)，且NDVI曲線具有明顯季節變化的草地、林地以及牧草等作物種植區域內，經過D-L擬合重建的NDVI具有較高的保真能力和適應性；3)在年均NDVI較低(NDVI<0.3)，且植被季節生長變化不明顯或NDVI曲線不呈季節對稱性變化的稀疏植被區，以及以冬小麥為典型作物種植的區域內，經過S-G濾波重建的NDVI數據表現出相對較好的保真能力和適應性。

NDVI；時間序列重建；植被遙感；AVHRR

植被在地球系統中扮演著重要的角色，是氣候和人文因素對環境影響的敏感指標[1-2]。遙感技術的快速發展為人類準確獲取陸地表面大范圍、多尺度的植被信息提供了強有力的手段。利用可見光/近紅外波段形成的歸一化植被指數(NDVI)已被廣泛地應用于區域和全球環境變化、農作物估產、干旱監測以及模擬新的植被指數等方面，是全球應用最廣的一種植被指數[3-6]。

獲取高質量、長時序NDVI數據不僅是長期監測陸地表面特征的基礎，也是研究氣候與陸地生態系統變化的前提。由于衛星傳感器成像本身具有瞬時性和周期性(衛星軌道周期和衛星生命周期)的特點，同時還會受到潛在噪聲的影響(如云、氣溶膠等)，因此在使用衛星傳感器來獲取NDVI時間序列數據時，數據的質量和連續性將會受到極大的挑戰[7]。特別是在大空間尺度的時間序列監測研究中，這一問題尤為嚴重。研究表明，大氣中的水汽、臭氧、氣溶膠和瑞利散射等會影響紅光波段與近紅外波段的反射率，使得衛星傳感器在接收到來自地面目標信號的同時，還會收到部分大氣噪聲信號，導致NDVI數值與實際情況出現較大的偏差[8-10]。盡管研究對象是植被，提取的是NDVI數據，但在植被生長的不同階段，傳感器實際接收到的信息也包括除植被以外的土壤背景。而土壤背景影響的本質是土壤本身對紅光和近紅外波段具有不同的反射率，包括潮濕地面、雪、枯葉、粗糙度、土壤類型等因素[11]。同時，在每次觀測時“太陽-地物-傳感器”三者之間的關系變化都可能直接影響植被指數計算帶來的不確定性[12-13]。因此，建立高質量、長時序的NDVI數據集是進一步開展植被定量遙感監測、氣候變化及陸地生態系統響應等相關研究的關鍵。

雖然目前常用的時間序列NDVI數據集在發布前都經過了較為嚴格的預處理來控制數據的質量，但由于受云、氣溶膠、土壤背景、太陽/觀測角度、數據傳輸錯誤、植被空間結構以及軌道和傳感器信號衰退等因素的影響，使這些數據集中仍然殘留較多噪聲，導致時間序列NDVI曲線出現陡升或者陡降的情況，從而影響了數據集的質量與連續性，造成在宏觀生態系統研究中基礎觀測數據存在不可靠的問題[14]。尤其是AVHRR傳感器，在設計之初并不是以監測地表特征為目的，導致其生成的NDVI時間序列數據更容易受到潛在噪聲的影響[15-17]。作為AVHRR NDVI的完善和升級，自2000年美國航空航天局開始免費分發全球MODIS NDVI產品后，使其快速成為應用研究的熱點[18-19]。雖然MODIS NDVI具有較窄的波段設置和較高的精度，但時間序列長度有限，在年代際尺度的變化趨勢研究中，仍然需要結合AVHRR數據對其進行時間序列的擴展和完善[20]。然而由于AVHRR數據已經積累了35年的時間長度，是目前時間序列跨度最長的光學遙感數據，對研究全球變化具有重要作用。因此，只有不斷訂正AVHRR NDVI數據集中的噪聲，插補缺測資料，發揮該數據長時序的優勢，為相關領域研究提供高質量的NDVI數據集將具有十分重要的現實意義。

時間序列重建是利用多種統計和數值分析方法，實現數據的插補和去噪，并提高數據質量的一種手段，目前已被成功地應用于NDVI數據時序重建工作中[21-22]。迄今為止，國內外學者在NDVI時間序列數據集重建方面開展了大量研究，目前時序重建方法已超過20種，依據計算方式可將其歸并為：閾值法、濾波法、函數擬合法以及綜合方法等[8,11]。雖然重建方法眾多，但每種重建方法對不同性質噪聲點去噪、平滑和保真能力不同。Hird等[23]通過6種重建算法對MODIS 16d NDVI產品進行重建比較分析發現，6種重建方法均有較好的去噪能力，但隨植被類型的變化重建效果具有明顯的差異性。黎治華[24]利用多種重建算法相結合的方案對MODIS的時間序列產品進行重構后發現，不同重建算法各有優劣，相關研究工作還存在很大的不確定性。由此可見，多種重建算法的提出，雖然為長時間序列NDVI數據集的構建提供了良好的技術支持，但對這些重建算法進行系統分析和比較的研究工作還不夠深入。對不同植被覆蓋或作物類型而言，研究者很少關注已有重建方法的適用性，在重建效果定量分析和精度驗證方面還需進一步開展相關研究工作。因此，系統比較和分析常用重建方法的擬合效果，開展以業務應用為目的重建算法適用性評價研究是很有必要的。

2014年6月美國戈達德太空中心發布了V4版本的LTDR(Land Long Term Data Record)產品，其中包括1981-2010年的地表反射率數據，這為建立高質量、長時序的NDVI數據集奠定了較好的基礎。因此，本研究以生態環境相對脆弱的西北地區東部為例，借助時間序列重建軟件(TIMESAT)[25]，通過結合農業氣象臺站資料，在系統比較與分析常用重建算法擬合效果的基礎上，對不同重建方法展開適用性評價，確定研究區內最佳擬合重建方法，以期為建立高質量、長時序的NDVI數據集提供參考。

1　材料與方法

1.1研究區概況

本研究以地處干旱區和半干旱區過渡帶的我國西北地區東部省份作為研究區。該地區介于 32°44′11″ N-44°07′24″ N、 88°47′17″ E-112°43′00″ E 之間，面積為1.48×106km2，約占我國陸地總面積的1/6。我國西北地區東部幅員遼闊，自然條件惡劣、生態環境較為脆弱，以山區降水和冰川融水補給為主，年降水量空間分布極不均勻，從南部450 mm以上過渡到北部沙漠地帶的50 mm以下，耕地灌溉主要依靠自然降水。特殊的地理位置和水熱氣候條件，造成研究區自然植被空間分布差異較大，形成了由西南草地到東北荒漠逐漸傾斜過渡的自然植被空間分布特征(圖1)。尤其是自20世紀70年代中期全球變暖加劇以來，該地區氣候條件復雜多樣，極端水文事件增加，地表植被覆蓋對環境變化敏感。因此，建立高質量、長時序的NDVI數據集，為準確監測該地區的植被生長狀況，了解氣候變化對植被的影響，維護干旱半干旱區生態系統的穩定具有重要指導意義。

圖1　西北地區東部植被分類圖Fig.1　The classification of vegetation in the east of Northwest China　SⅠ：牧草 Pasture；SⅡ：春小麥+玉米Spring wheat+Corn；SⅢ：冬小麥+玉米Winter wheat+Corn；SⅣ：冬小麥Winter wheat.

1.2研究數據與預處理

1.2.1LTDR地表反射率數據通過美國航空航天局戈達德太空中心(http://ltdr.nascom.nasa.gov/cgi-bin/ltdr/ltdrPage.cgi)，下載了1981-2010年的LTDR 每日地表反射率產品(AVH09 Surface Reflectance Product)。該數據是以NOAA07(1981-1985年)、NOAA09(1985-1988年)、NOAA11(1988-1994年)、NOAA14(1995-1999年)、NOAA16(2000-2006年)和NOAA18(2005-2010年)上的AVHRR觀測資料為基礎，通過軌道篩選、輻射定標、云檢測、大氣校正、衛星漂移校正及雙向反射率分布函數(bidirectional reflectance distribution function, BRDF)處理后生成的全球逐日格網數據集。數據分發采用與MODIS相類似的業務流程，以HDF分層格式存儲，包括1～3波段地表反射率、4～5波段大氣上界亮溫、太陽高度角、傳感器高度角、相對方位角以及質量控制文件。其中地表反射率數據是以空間分辨率為0.05°(約5 km)的格網進行存儲，質量控制文件對每個格網制定1～10級的質量標記。為了便于計算，本研究利用ENVI 5.1軟件將1～2波段地表反射數據處理為Geo-TIFF格式，投影方式轉為Albers，空間分辨率定義為5 km。

1.2.2植被類型數據植被類型數據來源于2001年中國科學院中國植被圖編輯委員會發布的中國植被數據集。該數據集是由我國著名植被生態學家侯學煜院士主編出版的《1∶1000000 中國植被圖集》數字矢量化而來，目前在國家自然科學基金委員會資助的“中國西部環境與生態科學數據中心”(http://westdc.westgis.ac.cn)提供下載共享，數據格式為Shapefile格式，投影方式為Albers。中國植被圖是開展自然資源和自然地理特征研究的基本圖件，它全面反映了我國11個植被類型組、54個植被型的796個群系和亞群系植被。在我國西北東部地區形成了以草地或草甸(以下簡稱草地植被)、稀疏植被、林地為主的植被分布格局(圖1)。

1.2.3氣象數據通過中國氣象科學數據共享服務網(http://cdc.nmic.cn)，下載了研究區101個農業氣象站在1992-2010年的逐旬資料。資料主要包括種植作物名稱、發育期名稱、發育期時間、發育程度等觀測指標。經初步統計分析表明，研究區主要種植作物分為4種類型，分別是牧草、冬小麥(Triticumaestivum)、冬小麥+其他種植作物[玉米(Zeamays)、馬鈴薯(Solanumtuberosum)等]、春小麥(Triticumaestivum)+其他種植作物(玉米、一季稻等)。

1.2.4典型作物樣區數據根據研究區作物種植特點，分別選取資料記錄較為完整，且具有典型代表性的青海河南SⅠ(牧草)、甘肅張掖SⅡ(春小麥+玉米)、甘肅西峰SⅢ(冬小麥+玉米)和陜西蒲城SⅣ(冬小麥)4個農業氣象臺站資料作為本研究結果的驗證數據。與此同時，為了盡量避免混合像元等不確定性因素對驗證結果的影響，本研究利用地面實測資料和Landsat8 OLI高分辨率資料，選取了臺站附近地勢較為平坦、耕地分布均勻，種植作物類型與臺站記錄相同、且種植面積大于5 km×5 km的下墊面作為典型樣區(圖1)。

1.3研究方法

雖然AVHRR傳感器均可提供時間分辨率為1 d的數據，但逐日數據易受云、氣溶膠及壞線等因素的影響。因此，本研究首先提取1981-2010年LTDR地表反射率產品中的1～2波段資料生成每日NDVI數據，其次結合數據的質量控制文件，根據最大合成法(max value composite, MVC)[26]，將逐日NDVI數據集中的明顯噪聲進行初步去除，并生成逐旬初始NDVI數據集，然后利用TIMESAT軟件對初始LTDR NDVI數據進行時序擬合重建，最后結合典型樣區資料和同時期高質量NDVI數據，對不同重建方法開展定量分析和適用性評價，確定研究區不同植被或作物類型上的最佳重建擬合方法。具體流程如圖2。

圖2　技術路線圖Fig.2　The flow chart of research method

1.3.1初始NDVI數據合成本研究利用LTDR 1～2波段地表反射率資料，根據NDVI計算公式生成1981-2010年每日的NDVI數據(O-N)，公式如下：

NDVI=(ρNir-ρRed)/(ρNir+ρRed)

(1)

式中，ρNir和ρRed分別為紅光波段和近紅外波段的反射率，LTDR資料相對應為第1波段和第2波段的地表反射率值。

結合數據自帶的質量控制文件(quality assessment, QA)，根據MVC合成法[26]，將每日LTDR NDVI合成初始逐旬NDVI數據，并生成相應時期的QA文件。QA文件是在像元尺度上能夠綜合評價和反映該時段遙感數據質量信息的資料。

1.3.2基于TIMESAT軟件的時間序列重建方法TIMESAT軟件是由J?nsson等[25]共同開發，并用于植被指數時間序列數據集重建及提取植被生長物候信息的程序包。軟件集成了目前已被國內外學者廣泛應用的非對稱高斯函數擬合(asymmetric Gaussian model function)、雙邏輯斯蒂擬合(double Logistic function)以及Savizky-Golay 濾波3種重建算法，并取得較好的效果[27-28]。該軟件有專門網頁介紹相關信息(http://www.nateko.lu.se/timesat/timesat.asp)，可同時處理時間序列數據和影像數據(二維空間陣列)。最初的版為FORTRAN語言編譯，之后由FORTRAN和MATLAB兩種語言編譯成升級兼容版本，在功能上有了較大的提升。主要包括數據預處理、數據處理以及后處理3個大模塊。可結合數據的質量控制文件實現不同遙感時序數據的重建。在數據預處理模塊中可以選擇其中一種或幾種重建算法進行重建，并可在同一窗口中顯示不同重建方法的預覽效果。TIMESAT軟件最大的特點是可根據不同的植被覆蓋類型設置不同的重建算法及參數，提取不同植被特征下的物候參數。

1)非對稱高斯函數擬合法(A-G)。A-G擬合法是一個從半局部擬合到整體擬合的過程，該算法可分為區間提取、局部擬合和整體連接三部分[22]。首先根據時間序列點相對應的值來提取區間，然后使用局部擬合函數即高斯函數對每個區間的時間序列數據進行擬合，最后通過整體擬合實現時間序列的重建，具體步驟如下。

第1步：通過局部擬合公式進行擬合，公式如下:

f(t;c1,c2,a1,……，a5)=c1+c2g(t;a1,……，a5)

(2)

式中，f(t;c1,c2,a1,……，a5)為擬合結果；c1和c2控制曲線的基準和幅度；g(t;a1,……，a5)為高斯函數，計算公式為：

(3)

式中：a1決定最大值和最小值的位置；a2、a3、a4、a5控制曲線左半部分和右半部分的寬度和陡峭度。

第2步：參數確定，通過極小化局部擬合公式進行擬合計算，公式如下：

(4)

式中：σi為第i個數據點的測量不確定性，可認為已知(通過質量控制文件獲取)，Ii為第i點的觀測值。

第3步：調整擬合曲線至植被生長的上包絡線：若擬合曲線的NDVI值低于原始觀測值，則通過調整低值擬合點的σi值進行再次擬合，從而使擬合曲線接近植被生長曲線的上包絡線。

第4步：通過整體擬合實現時間序列重建，將局部擬合的左邊最小值，中間最大值和右邊最小值分別表示為fL(t)，fC(t)和fR(t)，得到一個VI*(t)函數來糾正局部擬合函數，擬合公式為：

(5)

式中：α(t)和β(t)為裁切函數。

2)雙邏輯斯蒂擬合法(D-L)。D-L算法與A-G擬合法類似，也是一種半局部擬合方法。首先將整個時間序列中時間點對應的值按極大或極小值分成多個區間，分別對各區間進行雙邏輯斯蒂函數局部擬合，其局部擬合方式與非對稱高斯擬合方法類似，公式如下：

f(t;c1,c2,a1,……，a4)=c1+c2g(t;a1,……，a4)

式中：g(t;a1,……，a4)根據雙邏輯斯蒂函數擬合得到，計算公式為：

(6)

式中：a1和a3決定曲線左邊和右邊變化的位置，a2和a4決定曲線左邊和右邊的變化率。以上參數的獲取通過極小化x2獲得，如下式：

(7)

最后利用公式(5)進行整體擬合，并實現長時間序列的重建。

3)Savizky-Golay 濾波法(S-G)。S-G濾波是一種權重滑動平均濾波方法，主要是基于最小二乘卷積擬合法來實現噪聲消除的。其權重取決于一個濾波窗口內做多項式最小二乘擬合的次數，如下式：

(8)

式中：Y為原始的NDVI值；Y*為平滑后的NDVI值；Ci為平滑窗口中第i個NDVI值的權重系數；N為卷積數，等于平滑窗口的大小(2m+1)；指數j是原始坐標數據的滑動指數；m為滑動窗口的一半寬度。Ci可根據Madden等[29]提供的等式計算得到。具體步驟如下：

第1步：噪聲NDVI值的線性內插，通過質量控制文件識別標記受殘留噪聲影響的數據，使用線性內插對這些數據進行修正，得到新的時間序列數據Ni0。

第2步：使用S-G濾波擬合得到長時間變化的趨勢曲線Nitr；通過選擇適當的參數值m和j擬合NDVI序列長期變化的趨勢線。

第3步：確定NDVI時序每一個點的權重：根據以上兩步得到的新時間序列數據Ni0和序列的趨勢曲線Nitr來確定新時間序列上每個點的權重Wi，計算公式如下：

(9)

式中：Ni0和Nitr分別為新序列和趨勢線上第i點的NDVI值，di=|Ni0-Nitr|，dmax為di的最大值。

第4步：根據Ni0和Nitr生成新的時間序列Ni1，用于修正Ni0中存在的某些明顯的噪聲點，生成的新序列更真實地反映NDVI的變化趨勢，接近原始NDVI數據的上包絡線。計算公式為：

(10)

第5步：使用S-G濾波對Ni1曲線進行擬合，得到新的擬合曲線Ni2。

第6步：計算擬合效果指數Fk，用于估計擬合過程中擬合的NDVI值接近較高權重值的程度，擬合公式為：

(11)

本研究在利用TIMESAT軟件對初始逐旬LTDR NDVI數據進行A-G、D-L和S-G 3種擬合重建處理時，還需要對一些相關處理參數進行調試，主要包括：NDVI有效值域(range，0～10000)、噪聲去除閾值(spike，2)、滑動窗口大小(w-size，5～6～7)、擬合峰值參數(altitude，5)、迭代次數(fitting steps，3)等。

1.3.33種重建方法擬合效果比較分析為了定量分析以上3種重建方法的擬合效果和保真性能，本研究分別對典型耕地樣區和不同植被類型區的時序重建NDVI數據進行分析與比較。

1)典型耕地樣區擬合前后物候期比較分析。已有研究表明NDVI對植被的物候期非常敏感[30]，當作物從一個物候階段轉向另一個階段時，NDVI曲線曲率的變化能夠被用來確定物候期的轉變[30-32]。Zhang 等[33]將植被返青日期定義為NDVI的曲率導數首次達到局部最大值的日期。在作物生長周期內，通過計算NDVI擬合曲線的曲率，并求其極值，由此能夠確定植被生長階段轉變的時間起點有6個，分別對應返青、生長速率加快、成熟、衰老、衰老速率加快和枯黃死亡期(圖3)。曲率計算公式如下：

(12)

式中：da是沿時間曲線移動單位弧長時切線轉過的角度；ds為單位弧長；z=exp(a+bt)，a和b為時間t的擬合參數；c為NDVI最大值。

因此，根據以上植被物候期特征遙感識別方法[33]，本研究在典型樣區分別提取A-G、D-L、S-G以及初始NDVI時間序列曲線物候特征信息的基礎上，利用蒲城、海宴、河南、西峰和張掖4個樣區農業氣象臺站實際觀測的作物生育期資料，分析比較每個典型區內不同NDVI時間序列重建曲線表征的物候期(如返青期、抽穗和枯黃期)時間節點與農氣臺站實測作物生育期之間的誤差。為了使誤差定量可比，本研究將農氣臺站記錄的作物返青、成熟和枯黃的日期作為實測值，將通過NDVI重建曲線提取到的各相應生育期的時間節點作為模擬值，最后利用均方根誤差(RMSE)和平均絕對誤差(Pd)等相關統計量，定量分析模擬值與實測值之間的誤差精度。

其中，RMSE可以描述重建前后NDVI數值之間的平均差異程度，其值越小，擬合值的代表性越強；Pd能夠深入反映各算法彼此間保真能力的差異。RMSE和Pd公式如下：

(13)

(14)

式中：xi、yi分別表示實測樣本數據及其對應的重建樣本數據；Zai為第i個樣點的實際觀測值，Zei為重建擬合值；n為用于驗證的樣本數目。

2)不同土地類型區誤差定量分析。參照植被覆蓋類型和數字高程模型，在研究區下墊面較為均一、地勢較為平坦的草地、耕地、林地以及稀疏植被區分別選取初始NDVI時序數據中高質量的樣本數據(即QA=1的樣本數據)，并利用這些樣本數據定量分析重建前后NDVI數值之間的相關系數和均方根誤差等相關統計量。

2　結果與分析

2.1典型作物樣區擬合重建結果特征比較

通過分析1981-2010年青海河南(SⅠ)、甘肅張掖(SⅡ)、甘肅西峰(SⅢ)和陜西蒲城(SⅣ)4個典型作物樣區內的時間序列NDVI擬合重建數據，結果表明，3種擬合重建算法在對初始NDVI像元中的殘留噪聲進行消除的同時，也將典型像元中的時序NDVI進行了不同程度的擬合調整。其中經過D-L和A-G方法擬合曲線的結果比較接近，僅在NDVI峰值處兩者之間的擬合結果略有不同，相對誤差為±0.7%(圖4)。S-G濾波法對噪聲比較敏感，其重建的NDVI值與前兩種擬合重建結果有較大的差異，尤其是在生長季開始和結束時段，S-G擬合曲線相對于初始(O-N)、D-L和A-G擬合曲線的峰值存在較大程度的“左偏”和“右移”的情況(圖4)。

植被從返青期到枯黃期階段，經過3種重建算法擬合后的NDVI整體較初始NDVI值提高了11.6%～23.7%，且重建序列曲線“拐點”對應的時段也不一致。在生長季階段(4-10月)，青海河南和甘肅張掖典型樣區的初始NDVI曲線呈較為明顯的年內對稱性，且3種擬合重建后的結果也較為接近。但以冬小麥為主要作物種植區的陜西蒲城和甘肅西峰(圖4)，在生長季階段的NDVI曲線具有明顯的年內不對稱性，波峰波谷更替較為頻繁，且初始NDVI曲線與3種擬合重建后的曲線具有較大差異。由此可見，不同作物的種植結構對擬合重建效果具有較大影響。

2.2典型作物樣區擬合重建結果誤差分析

通過比較分析1992-2010年各典型區初始NDVI曲線(O-N)，以及3種重建NDVI曲線曲率所表征的物候期(返青、成熟和枯黃期)節點與相應臺站記錄作物生育期之間的差異，可以看出：由于受到噪聲的影響，初始NDVI曲線反映的各物候期節點與臺站記錄的相應生育期之間存在較大誤差，且隨著作物不斷生長，即由返青期到枯黃期階段，平均絕對誤差Pd和均方根誤差RMSE也隨之升高，其中Pd和RMSE分別介于2.36～3.68和1.73～5.95之間。

與初始NDVI曲線表征的各物候期相比，經過擬合重建后的3種NDVI曲線所反映的物候期節點與臺站記錄的相應生育期之間的各項誤差均有所降低，其中在青海河南(SⅠ)和甘肅張掖(SⅡ)典型區域內，通過D-L重建后的NDVI擬合曲線反映的各生育期節點與臺站記錄的相應時期最為接近，Pd和RMSE達到最小，分別在1.89和4.21以下；在陜西蒲城(SⅣ)以及甘肅西峰(SⅢ)典型區域內，經過S-G濾波重建的NDVI曲線所反映的各生育期節點與實測值之間具有較小Pd和RMSE，分別介于2.35～2.47和1.57～4.69之間(表1)。

由此可見，3種擬合重建算法均可以對初始NDVI數據中的殘留噪聲進行消除，但根據不同作物生長特點或作物種植結構，每種重建方法對其噪聲消除能力是不同的。其中在生長季階段(4-10月)，作物生長曲線具有明顯年內對稱性的SⅠ和SⅡ典型區域內，經過D-L擬合重建結果在相應生育期內的各項誤差較小，其適應性較好；而在以冬小麥為典型種植區的SⅢ和SⅣ區域內，經過S-G濾波擬合重建的NDVI曲線對各生育期節點的反映與實測值更為接近，具有較好的適應性。

2.3不同植被或作物類型NDVI擬合重建結果誤差分析

根據NDVI質量控制文件，在獲取1981-2010年各植被類型區高質量(QA=1)NDVI樣本數據，以及經過上述3種擬合重建后同期NDVI樣本數據的基礎上，通過比較分析重建前后NDVI數據之間的相關系數和RMSE(圖5)，結果表明，在草地(Ⅰ)和林地(Ⅲ)為主的自然植被區域內，以及牧草(Ⅳ)和春小麥+其他作物(Ⅴ)的耕地作物區域中，經過D-L擬合重建的NDVI數據與重建前同期高質量NDVI數據之間的相關系數達到最高(≥0.680)，且RMSE較低(≤0.026)。經過A-G擬合后的結果與D-L較為接近，而S-G擬合結果與同期高質量NDVI數據之間具有較大差別，其相關系數和RMSE分別介于0.620～0.780和0.019～0.028之間。與此同時，在稀疏植被區(Ⅱ)、冬小麥+其他作物(Ⅵ)和冬小麥(Ⅶ)區域中，經過S-G濾波擬合重建后的NDVI數據與同期重建前高質量NDVI數據之間具有較高的相關系數(≥0.56)和較低的RMSE(≤0.027)。

圖4　1981-2010年典型樣區逐旬NDVI時間序列重建前后對比Fig.4　Comparison of original and reconstructed NDVI of ten-days period from 1981 to 2010 in sample areas

表1　典型作物樣區NDVI曲線物候期誤差分析

圖5　不同植被或作物類型區內高質量NDVI重建前后保真能力分析Fig.5　Comparison of the ability on keeping the main characters by NDVI in different vegetation coverage areasⅠ～Ⅶ分別表示草地、稀疏植被、林地、牧草、春小麥+其他作物、冬小麥+其他作物和冬小麥區Ⅰ to Ⅶ are respectively grassland, sparse vegetation, woodland, pasture, spring wheat+other crops, winter wheat+other crops, and winter wheat.

由此可見，隨著下墊面植被或作物類型的變化，3種擬合重建結果具有較大的差異性。其中，在年均NDVI較高(NDVI≥0.3)的草地、林地區、牧草和春小麥+其他作物區域中，經過D-L擬合重建的NDVI數據具有較高的保真能力和適應性，而在年均NDVI較低(NDVI<0.3)的稀疏植被和冬小麥+其他作物區域中，S-G濾波重建方法表現出相對較好的保真能力和適應性。

3　討論

綜上所述，由于目前已發布的時間序列 LTDR NDVI數據中仍然殘留較多噪聲，會造成NDVI曲線出現陡升或者陡降的情況，由此導致植被生長曲線表征的各物候期時間節點與臺站實際記錄的生育期之間存在較大偏差。因此，在實際應用過程中需要對該數據集中的殘留噪聲進行消除。利用時間序列重建算法對NDVI數據集進行擬合重建，不僅能夠消除數據中的殘留噪聲，同時還可以插補缺測數據，是提高數據集質量，獲取連續時間序列的一種有效手段。

根據本研究可以看出，當植被生長具有明顯的季節變化特征時(如草地或部分作物區域)，利用D-L或A-G擬合方法重建的NDVI時間序列曲線，對返青期、成熟期以及枯黃期等時間節點的反映較為準確，并與重建前高質量NDVI數據具有較小的誤差和較高的保真性。這主要是由于A-G與D-L擬合方法均是一個從半局部擬合到整體擬合的過程，兩者擬合時序數據的出發點是相同的，只是后期分別使用了高斯函數和邏輯斯蒂函數進行擬合。因此，經過以上兩種方法得到的NDVI曲線結果較為接近，僅僅只是對NDVI峰值的擬合存在差異。而S-G濾波重建主要是利用移動加權滑動平均濾波對整個NDVI時序數據的長期變化趨勢進行擬合，并將NDVI值分為“真值”點和“偽真值”點兩類。最后通過局部迭代擬合的方式將“偽真值”點用濾波值進行代替，從而得到新的擬合曲線。與A-G與D-L擬合方法相比較，S-G濾波方法對曲線局部細節的波動更為敏感，因此，在植被季節生長變化特征不明顯的荒漠和部分耕地作物地區，利用S-G濾波方法重建的NDVI曲線具有較好的擬合效果。Hird等[23]在加拿大落基山脈地區，利用MODIS 16d NDVI產品進行NDVI重建比較分析發現，A-G和D-L擬合結果優于S-G。曹云鋒等[27]在長白山自然保護區進行NDVI重建比較研究也得出了以上相似結果。而在本研究中可以看出，只有在植被生長具有明顯季節變化特征的草地、林地、牧草和部分作物類型上，經過A-G和D-L擬合方法才具有較好的擬合效果，在其他植被類型條件下，S-G濾波法均優于A-G和D-L的擬合效果。宋春橋等[28]對我國西藏北部地區的MODIS NDVI數據進行重建分析發現，A-G和D-L方法在草原和灌叢類型上具有較好的重建效果，相對而言，S-G濾波對噪聲比較敏感，在處理局部較強的噪聲時，可能使其“過度”擬合而引入新的噪聲。本研究一方面通過重建前高質量NDVI數據和GIS空間分析等技術手段，從整個空間尺度上對上述觀點進行了定量化的分析與比較，另一方面利用相應時期的農業氣象資料，反映出植被季節生長曲線變化特征對不同重建方法的擬合重建效果具有重要影響。由此可見，在系統比較和分析多種重建方法的基礎上，融合多種重建方法的優勢，構建一種綜合重建方法是下一步需要完善的工作。

由于本研究使用的農氣資料記錄時間有限，而無法對1981-1991年的NDVI數據進行定量分析，同時在自然植被類型區，由于缺乏長期的地面觀測資料，造成精度評價并不完整。因此，收集更長時間序列的相關實測資料，從各個方面定量分析評價不同擬合重建方法的結果，對建立高質量、長時序的NDVI重建數據集具有重要的現實意義。

4　結論

本研究以生態環境相對脆弱的西北東部地區為例，結合農業氣象臺站資料和高質量NDVI數據，在系統比較與分析A-G、D-L和S-G 3種重建算法擬合效果的基礎上，對不同重建方法展開適用性評價，得出如下結論：

1)下墊面類型是影響不同重建方法效果的一個重要因素。經過3種擬合重建算法獲取的NDVI數據，均可以對初始NDVI數據中的殘留噪聲進行消除，并提高數據集的質量，但根據不同植被類型或作物生長特點，每種重建方法對其噪聲消除能力是不同的。

2)在年均NDVI較高(NDVI≥0.3)，且具有明顯季節變化特征的草地、林地、牧草以及春小麥+其他作物區域中，經過D-L擬合重建的NDVI曲線，在植被或作物的返青期、成熟期以及枯黃期等時間節點上與臺站實測生育期較為接近，且與重建前高質量NDVI數據的誤差較小，具有較高的保真能力和適應性。A-G擬合重建結果與D-L結果較為接近，只是根據D-L和A-G擬合算法重建的NDVI峰值結果略有差異。

3)在年均NDVI較低(NDVI<0.3)，且在生長季內植被季節生長變化特征不明顯或曲線不具對稱性的稀疏植被區，或以冬小麥為典型種植作物的區域中，經過S-G濾波重建的NDVI數據與臺站實測資料較為接近，且與重建前高質量NDVI數據的誤差較小，表現出相對較好的保真能力和適應性。

References：

[1]Sudipta S, Menas K. Interannual variability of vegetation over the Indian sub-continent and its relation to the different meteorological parameters. Remote Sensing of Environment, 2004, 90: 268-280.

[2]Li F, Zeng Y, Li X S,etal. Remote sensing based monitoring of interannual variations in vegetation activity in China from 1982 to 2009. Science China: Earth Sciences, 2014, 57: 1800-1806.

[3]Hou M T, Zhao H Y, Wang Z,etal. Vegetation responses to climate change by using the satellite-derived normalized difference vegetation index: A Review. Climatic and Environmental Research, 2013, 18(3): 353-364.

[4]Arnon K, Nurit A, Rachel T P,etal. Use of NDVI and land surface temperature for drought assessment: merits and limitations. Journal of Climate, 2010, 23: 618-633.

[5]Barbosa H A, Huete A R, Baethgen W E. A 20-year study of NDVI variability over the Northeast region of Brazil. Journal of Arid Environments, 2006, 67: 288-307.

[6]Guo N, Wang X P. Advances and developing opportunities in remote sensing of drought. Journal of Arid Meteorology, 2015, 33(1): 1-18.

[7]Piao S L, Wang X H, Ciais P. Changes in satellite-derived vegetation growth trend in temperate and boreal Eurasia from 1982 to 2006. Global Change Biology, 2011, 17(10): 3228-3239.

[8]Gu J, Li X, Huang C L. Research on the reconstructing of Time-series NDVI Data. Remote Sensing Technology and Application, 2006, 21(4): 391-395.

[9]Han P, Yao J, Li T H. Comparison of 3 NDVI datasets and the application at Yanhe basin, China. Journal of Basic Science and Engineering, 2014, 4: 661-674.

[10]Wang Z X, Liu C, Huete A. From AVHRR-NDVI to MODIS-EVI: Advances in vegetation index research. Acta Ecologica Sinica, 2003, 23(5): 979-987.

[11]Geng L Y, Ma M G. Advance in method comparison of reconstructing remote sensing time series data sets. Remote Sensing Technology and Application, 2014, 29(2): 362-368.

[12]Ma M G, Song Y, Wang X F,etal. Development status and application research of the time series remote sensing data products based on AVHRR、VEGETATION and MODIS. Remote Sensing Technology and Application, 2012, 27(5): 663-670.

[13]Du J Q, Shu J M, Wang Y H,etal. Comparison of GIMMS and MODIS normalized vegetation index composite data for Qinghai-Tibet Plateau. Chinese Journal of Applied Ecology, 2014, 25(2): 533-544.

[14]Michishita R, Zhenyu J, Jin C,etal. Empirical comparison of noise reduction techniques for NDVI time-series based on a new measure. ISPRS Journal of Photogrammetry and Remote Sensing, 2014, 91: 17-28.

[15]Beck P S, Goetz S J. Satellite observations of high northern latitude vegetation productivity changes between 1982 and 2008: Ecological variability and regional differences. Environmental Research Letters, 2011, 6: 5501-5511.

[16]Jeremy L W, David S G, Julia E A,etal. Long-term vegetation monitoring with NDVI in a diverse semi-arid setting, central New Mexico, USA. Journal of Arid Environments, 2004, 58(2): 249-272.

[17]Steven M D, Timothy J M, Frédéric B,etal. Intercalibration of vegetation indices from different sensor systems. Remote Sensing of Environment, 2003, 88(4): 412-422.

[18]Chen Y L, Long B J, Pan X B,etal. Differences between MODIS NDVI and AVHRR NDVI in monitoring grasslands change. Journal of Remote Sensing, 2011, 15(4): 831-845.

[19]Fensholt R, Proud S R. Evaluation of earth observation based global long term vegetation trends—Comparing GIMMS and MODIS global NDVI time series. Remote Sensing of Environment, 2012, 119(16): 131-147.

[20]Mao D H, Wang Z M, Luo L,etal． Integrating AVHRR and MODIS data to monitor NDVI changes and their relationships with climatic parameters in Northeast China. International Journal of Applied Earth Observation and Geoinformation, 2012, 18: 528-536.

[21]Geng L, Ma M, Wang X,etal. Comparison of eight techniques for reconstructing multi-satellite sensor time-series NDVI data sets in the heihe river basin, China. Remote Sensing, 2014, 6: 2024-2049.

[22]Li R, Zhang X, Liu B,etal. Review on methods of remote sensing time-series data reconstruction. Journal of Remote Sensing, 2009, 13(2): 335-341.

[23]Hird J N, McDermid G J. Noise reduction of NDVI time series: An empirical comparison of selected techniques. Remote Sensing of Environment, 2009, 113: 248-258.

[24]Li Z H. A Study on the Eco-environment Evolution of Yangtze River Delta Region Based on the Retrieval & Reconstruction of MODIS Time Series Datasets[D].Shanghai: East China Normal University, 2011.

[25]J?nsson P, Eklundh L. TIMESAT-a program for analyzing time-series of satellite sensor data. Computers & Geosciences, 2004, 30(8): 833-845.

[26]Holben B N. Characterististics of maximun-value composite images from temporal AVHRR data. International Journal of Remote Sensing, 1986, 7(11): 1417-1434.

[27]Cao Y F, Wang Z X, Deng F P. Fidelity performance of three filters for high quality NDVI time-series analysis. Remote Sensing Technology and Application, 2010, 1(1): 118-125.

[28]Song C Q, You S C, Ke L H,etal. Analysis on three NDVI time-series reconstruction methods and their applications in North Tibet. Journal of Geo-Information Science, 2011, 1(1): 133-143.

[29]Madden H H, Chem A. Comments on the Savitzky-Golay convolution method for least-squares-fit smoothing and differentiation of digital data. Analytical Chemistry, 1978, (9): 1383-1386.

[30]Moulin S, Kergoat L, Viovy N,etal. Global-scale assessment of vegetation phenology using NOAA/AVHRR satellite measurements. Journal of Climate, 2010, 10(6): 1154-1170.

[31]Liu H. Spring Phenology Model of Grassland Based on Soil Moisture and Air Temperature and Vegetation Reactions to Drought[D]. Beijing: Tsinghua University, 2012.

[32]Wang L, Ding J L. Vegetation index feature change and its influencing factors and spatial-temporal process analysis of desert grassland in the Ebinur Lake Nature Reserve, Xinjiang. Acta Prataculturae Sinica, 2015, 24(5): 4-11.

[33]Zhang X, Friedl M A, Schaaf C B,etal. Monitoring vegetation phenology using MODIS. Remote Sensing of Environment, 2003, 84(3): 471-475.

[3]侯美亭, 趙海燕, 王箏, 等. 基于衛星遙感的植被NDVI對氣候變化響應的研究進展. 氣候與環境研究, 2013, 18(3): 353-364.

[6]郭鈮, 王小平. 遙感干旱應用技術進展及面臨的技術問題與發展機遇. 干旱氣象,2015, 33(1): 1-18.

[8]顧娟, 李新, 黃春林. NDVI 時間序列數據集重建方法述評. 遙感技術與應用, 2006, 21(4): 391-395.

[9]韓鵬, 姚娟, 李天宏. 3種不同數據源NDVI的比較分析及其在延河流域的應用研究. 應用基礎與工程科學學報, 2014, 4: 661-674.

[10]王正興, 劉闖, Huete A. 植被指數研究進展: 從AVHRR-NDVI到MODIS-EVI. 生態學報, 2003, 23(5): 979-987.

[11]耿麗英, 馬明國．長時間序列NDVI數據重建方法比較研究進展．遙感技術與應用, 2014, 29(2): 362-368.

[12]馬明國, 宋怡, 王旭峰, 等. AVHRR、VEGETATION和MODIS時間序列遙感數據產品現狀與應用研究進展. 遙感技術與應用, 2012,27(5): 663-670.

[13]杜加強, 舒儉民, 王躍輝, 等. 青藏高原MODIS NDVI與GIMMS NDVI的對比. 應用生態學報, 2014, 25(2): 533-544.

[22]李儒, 張霞, 劉波, 等. 遙感時間序列數據濾波重建算法發展綜述.遙感學報, 2009, 13(2): 335-341.

[24]黎治華. 基于MODIS反演重構時間序列數據的長江三角洲地區生態環境演變研究[D]. 上海： 華東師范大學, 2011.

[27]曹云鋒, 王正興, 鄧芳萍. 3種濾波算法對NDVI高質量數據保真性研究. 遙感技術與應用, 2010, 1(1): 118-125.

[28]宋春橋, 游松財, 柯靈紅, 等. 藏北地區三種時序NDVI重建方法與應用分析. 地球信息科學學報, 2011, 1(1): 133-143.

[31]劉慧. 基于土壤水分和氣溫的草地返青模型及植被干旱研究[D]. 北京： 清華大學, 2012.

[32]王璐, 丁建麗. 艾比湖保護區荒漠植被時空過程變化及其植被指數影響因素分析. 草業學報, 2015, 24(5): 4-11.

*Comparative studies of reconstruction methods for the long term NDVI dataset in the east of Northwest China

WANG Wei, GUO Ni*, SHA Sha, HU Die, WANG Xiao-Ping, LI Yao-Hui

KeyLaboratoryofAridClimaticChangeandReducingDisasterofGansuProvince,KeyOpenLaboratoryofAridChangeandDisasterReductionofCMA,InstituteofAridMeteorology,ChinaMeteorologicalAdministration,Lanzhou730020,China

A high-level time-series NDVI dataset is not only the basis for continuous monitoring of the land surface, but also an important tool for studying change related to climate and land use factors in terrestrial eco-systems. We reconstructed the noise component of the LTDR NDVI data for the east of Northwestern China where the ecosystem is fragile, using various time-series reconstruction methods. This paper use agrometeorological data and high-level NDVI data to evaluate the accuracy of different reconstruction methods. The results show that: 1) The vegetation or crop land cover is an important factor affecting fitted results of the various reconstruction methods. Each reconstruction method has a different noise reduction ability depending on differences in vegetation or crop growth characters; 2) The D-L reconstruction method has a better noise reduction ability and applicability in those areas of grassland, and woodland for which the annual average NDVI data is higher (NDVI≥0.3) and the NDVI curve has obvious seasonal changes; 3) The S-G reconstruction method has better fidelity ability and applicability in some areas of crop land in winter wheat and in areas of sparse vegetation for which annual average NDVI data are lower (NDVI<0.3) and where the NDVI curve have no obvious seasonal changes.

NDVI; time-series data reconstruction; vegetation remote sensing; AVHRR

10.11686/cyxb2015489http://cyxb.lzu.edu.cn

王瑋, 郭鈮, 沙莎, 胡蝶, 王小平, 李耀輝. 我國西北地區東部時間序列NDVI數據集重建方法比較研究. 草業學報, 2016, 25(8): 1-13.

WANG Wei, GUO Ni, SHA Sha, HU Die, WANG Xiao-Ping, LI Yao-Hui. Comparative studies of reconstruction methods for the long term NDVI dataset in the east of Northwest China. Acta Prataculturae Sinica, 2016, 25(8): 1-13.

2015-10-22；改回日期：2016-01-04

甘肅省氣象局氣象科研項目(GSMAMs2016-10)，中國博士后科學基金項目(2015M582734)和公益性行業(氣象)科研專項(重大專項)(GYHY201506001-5)資助。

王瑋(1985-)，男，甘肅蘭州人，助理研究員。 E-mail: wangwei9969@163.com

Corresponding author. E-mail: guoni0531@126.com