基于Hyb-F組合濾波算法的向海自然保護區(qū)NDVI時間序列重構(gòu)

劉 舒，姜琦剛，朱 航，李曉東,3

(1.吉林大學(xué) 地球探測科學(xué)與技術(shù)學(xué)院, 長春 130026；2.吉林大學(xué) 機械科學(xué)與工程學(xué)院, 長春 130022； 3.白城師范學(xué)院 旅游與地理科學(xué)學(xué)院, 吉林 白城137000)

0 引 言

歸一化植被指數(shù)(NDVI)時間序列是反映植被與地表氣候條件、人類活動間響應(yīng)關(guān)系的重要工具。然而由衛(wèi)星直接獲取的遙感數(shù)據(jù)得到的NDVI時間序列中存在來自云、冰雪、陰影、氣溶膠雙向反射等多種噪聲，在應(yīng)用前應(yīng)去除噪聲，重構(gòu)時間序列[1-3]。

常見NDVI序列重構(gòu)方式有：基于局部數(shù)據(jù)的閾值法、最大值合成法(MVC)、Savitzky-Golay濾波法等；基于全局或半全局數(shù)據(jù)的非對稱高斯擬合法、雙邏輯函數(shù)擬合法[4，5]和時間序列諧波分析法[6，7]等。不同重構(gòu)方法對不同區(qū)域或不同數(shù)據(jù)源的適用性不同，多位學(xué)者提出了各算法的對比結(jié)果，但各算法優(yōu)劣目前并無定論[2]。

目前已有學(xué)者對一些單一類型數(shù)據(jù)重構(gòu)算法進行改進[3,8]，或同時利用兩類濾波算法構(gòu)建組合式濾波算法[2,9,10]，提高NDVI時間序列重構(gòu)精度和重構(gòu)過程的自動化程度，但相關(guān)研究仍存在以下局限。首先，局部濾波算法較好地保留序列細部特征，但濾波結(jié)果不遵循植被生長規(guī)律，而全局濾波算法所得結(jié)果符合植被生長規(guī)律，但在局部存在過度平滑現(xiàn)象。其次，兩類算法所得結(jié)果都容易受噪聲影響，且對噪聲過濾能力有限。為更形象地描述濾波效果，將噪聲分為兩類，一類是能被局部濾波算法識別的局部噪聲，另一類為能被全局濾波算法識別的不符合植被總體生長趨勢的噪聲，稱為全局噪聲。最后，在目前的組合濾波法中，各組成部分仍相對獨立，去噪過程由局部濾波法完成，而未能被局部算法識別的全局噪聲則被保留，繼續(xù)參與NDVI序列重構(gòu)。以上局限都在一定程度上影響重構(gòu)序列的精度。

本文構(gòu)建了一種新型組合濾波算法Hyb-F (Hybrid filtering method based on Grubbs-introduced S-G and AG algorithms)，將局部S-G濾波和全局AG擬合算法相結(jié)合，并在各自過程中引入Grubbs檢測算法。該組合算法優(yōu)點在于分步剔除序列局部噪聲和全局噪聲，構(gòu)建具有精細細部特征，同時反映植被物候變化的NDVI時間序列，并以向海自然保護區(qū)為研究區(qū)，分析組合算法的性能。

1 材料與方法

1.1 研究區(qū)地理概況

向海自然保護區(qū)位于吉林省通榆縣西北部，地理范圍為122°05′～122°35′E，44°50′～45°19′N(圖1(a)，保護區(qū)范圍由通榆縣地圖矢量化獲得)，是吉林省唯一列入《國際重要濕地名冊》的重要濕地。保護區(qū)內(nèi)多種珍稀野生動植物繁育，多樣生態(tài)景觀共存[11]。向海自然保護區(qū)內(nèi)的地表覆被類型如圖1(b)所示，主要是以蒙古黃榆、大果榆為主的灌木林地，楊樹為主的喬木林地，以蘆葦、香蒲為主的沼澤濕地，以羊草為主的草地，以玉米田為主的耕地，大片水域，出露的重度鹽堿地和零星分布的城鎮(zhèn)工礦建設(shè)用地[12,13]。

圖1 向海自然保護區(qū)Fig.1 Xianghai natural reserve

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

本研究選取Landsat8 陸地成像儀(OLI)獲取的多光譜影像組成基礎(chǔ)數(shù)據(jù)集。影像從USGS網(wǎng)站(https://glovis.usgs.gov/)下載，為L1T級數(shù)據(jù)，已經(jīng)過系統(tǒng)輻射校正和精細幾何校正，空間分辨率為30 m，時間分辨率為16天。覆蓋研究區(qū)影像編號為120/029(帶/行)，為盡量保證影像獲取年份一致，避免數(shù)據(jù)集中連續(xù)多期影像存在云雪“污染”，最終選取2015年1月至2016年12月獲取的數(shù)據(jù)構(gòu)成影像序列集，個別時段連續(xù)多期數(shù)據(jù)云覆蓋度大，由2014年相應(yīng)時段數(shù)據(jù)補充。共選取27期影像，基本信息如表1所示。

將影像在ENVI5.1中進行大氣校正，將投影參數(shù)定義為WGS84橢球參數(shù)下的UTM投影(Zone51)。分別利用各期紅光波段和短波紅外波段計算相應(yīng)影像的NDVI。其中，由2015年5月13日數(shù)據(jù)得到的NDVI不參與擬合過程，僅作為輔助數(shù)據(jù)，驗證數(shù)據(jù)序列重構(gòu)實驗效果。將數(shù)據(jù)按年內(nèi)時間順序組合成一幅具有26個波段的NDVI時間序列影像。

表1 Landsat8 OLI數(shù)據(jù)基本特征(行列號:120/029)Table 1 Main characters of Landsat8 OLI imagery (path/row:120/029)

1.3 研究方法

1.3.1 設(shè)定指標閾值剔除噪聲

直接算得的NDVI序列存在超出理論值域等與實際不符的無效點，因此，首先將移動窗口長度設(shè)為5，通過設(shè)定曲線各點取值、變化量、變化量與變化方向的關(guān)系、雪蓋指數(shù)等指標的閾值，識別移動窗口中超出理論值的噪聲和部分云、雪噪聲。剔除窗口中噪聲，并以其他值的二次多項式最小二乘擬合曲線在噪聲點的取值代替原始值。

1.3.2 Grubbs異常點檢驗

Grubbs檢驗是一種常見的離群異常值檢驗方法，通過檢驗觀測值殘差的Grubbs統(tǒng)計量與某水平下的臨界值之間的關(guān)系，判斷異常點位置。

本研究的噪聲點在剔除過程中，分別將Grubbs檢驗過程與S-G濾波算法和AG擬合算法相結(jié)合，判斷時間序列中存在的局部異常點和全局異常點。假定進行S-G濾波前后或AG擬合過程前后時間序列的差值序列為NDVIC，Grubbs檢驗具體實現(xiàn)方法為：將NDVIC由小到大排列，并計算各點Grubbs統(tǒng)計量，查表獲得顯著性水平α下的上下臨界值，找出NDVIC最大值和最小值點中Grubbs統(tǒng)計量較大的點，判斷該點Grubbs統(tǒng)計量與臨界值的關(guān)系。若統(tǒng)計量超過臨界值的范圍，則該點為異常點[8]。將Grubbs檢測的顯著性水平選為0.05，為防止對異常點的誤判斷，本文通過迭代過程實現(xiàn)異常點剔除，每次迭代僅剔除一個偏差最嚴重的點。在兩步去噪過程中，噪聲點的去噪后取值分別為以該點為中心的S-G移動窗口內(nèi)基于其他點的S-G濾波曲線在該點的取值，或剔除該點后的AG擬合曲線在該點的取值。剔除一個異常點后，得到新的待降噪NDVI序列，并計算新的NDVIC序列，繼續(xù)Grubbs檢驗，直到將所有異常點剔除。

1.3.3 Savitzky-Golay濾波算法

Savitzky-Golay濾波算法(S-G)是一種基于最小二乘原理對滑動窗口內(nèi)數(shù)據(jù)進行卷積運算的數(shù)據(jù)平滑方式。對于以第j點為中心，寬度為N的滑動窗口，S-G算法基本數(shù)學(xué)表達式為：

(1)

式中：qi為原始序列第i個數(shù)據(jù)取值；pj為平滑后序列第j個數(shù)據(jù)取值；N為滑動窗口長度；m為窗口中間至一側(cè)端點的長度，N=2m+1；Ci為第i個數(shù)據(jù)的系數(shù)。

當(dāng)高質(zhì)量數(shù)據(jù)在整個時間序列所占比例較高時，S-G算法能夠過濾序列噪聲，較細致地描述序列局部特征[14]。

在實際運算中，S-G濾波曲線由移動窗口內(nèi)數(shù)據(jù)二次多項式的最小二乘擬合得出。由于時間序列中存在噪聲污染的影像所占比例較大，存在局部集中情況，且無其他質(zhì)量控制數(shù)據(jù)，將研究中涉及的S-G濾波過程窗口長度設(shè)為7。

本研究中兩次運用S-G濾波算法，分別用于局部異常值的剔除與最終時序數(shù)據(jù)的重構(gòu)。

1.3.4 非對稱高斯擬合算法

非對稱高斯擬合算法(Asymmetric Gussian，AG)是一種基于高斯函數(shù)擬合的非線性重構(gòu)算法。該算法根據(jù)數(shù)據(jù)序列極值點所在位置將數(shù)據(jù)分段，對每段數(shù)據(jù)分別進行高斯函數(shù)擬合，并利用區(qū)間端點附近的平滑截斷函數(shù)連接每段曲線，實現(xiàn)全局擬合。整個擬合過程分以下2步完成：

首先，將兩個波峰或波谷間的數(shù)據(jù)劃分為左、中、右3部分，3個峰值點位置由tL、tC和tR表示。利用局部高斯平滑函數(shù)分別擬合每一段曲線，該函數(shù)的數(shù)學(xué)表達式為：

f(t)=f(t;b1,b2,a1,a2,…,a5)=

b1+b2g(t;a1,a2,…,a5)

(2)

式中：高斯函數(shù)的具體形式為：

g(t;a1,a2,…,a5)=

(3)

參數(shù)b1和b2控制高斯函數(shù)g的基線位置和振動幅度，a1控制局部擬合區(qū)域函數(shù)峰值的位置，a2、a3、a4和a5分別決定局部擬合區(qū)域右、左部分函數(shù)的寬度和陡緩程度。

最后，利用全局方程F(t)銜接整個[tL,tR]時段各局部函數(shù)。全局函數(shù)數(shù)學(xué)表達式為：

F(t)=

(4)

式中：fL(t)、fC(t)和fR(t)分別為左、中、右3部分的局部高斯擬合函數(shù)；α(t)和β(t)為由0到1平滑過渡的截斷函數(shù)[4，5]。

在本研究中，非對稱高斯擬合主要與Grubbs檢測算法相結(jié)合，應(yīng)用于全局噪聲的剔除。

Hyb-F組合濾波算法具體流程如圖2所示。

1.3.5 NDVI時間序列重建效果評價方法

為評價Hyb-F組合濾波方式對NDVI時間序列的重建效果，在圖1(b)所示覆被區(qū)選擇樣本點，提取各層NDVI取值，利用Hyb-F算法重構(gòu)樣本NDVI時間序列，并采用AG算法和S-G濾波算法進行對比實驗。AG算法和S-G濾波算法的迭代過程對剔除低值噪聲，提高擬合精度起到較大作用。因此，實際對比了Hyb-F濾波結(jié)果與原始序列(Original)、S-G直接濾波結(jié)果(S-G Direct)、S-G三次迭代濾波結(jié)果(S-G with 3 iterations)、AG直接濾波結(jié)果(AG Direct)和AG三次迭代濾波結(jié)果(AG with 3 iterations)。

本研究中，對不同方法的數(shù)據(jù)重構(gòu)效果分別采用視覺直觀驗證和樣本點重構(gòu)序列統(tǒng)計指標定量分析兩種方式進行對比。

統(tǒng)計指標定量分析的基礎(chǔ)在于構(gòu)建準確的參考時間序列。由于實際條件限制，觀測量會有一定偏差，研究者很難通過其獲取真實的NDVI時間序列。本文對比實驗定量評價基于以下假定：

圖2 時間序列濾波算法流程圖Fig.2 Workflow of the hybrid filtering algorithm

若重構(gòu)算法不能完全剔除原始序列中被云和雪污染的數(shù)據(jù)點，則擬合過程不能成功濾除原有高質(zhì)量數(shù)據(jù)中存在的微小偏差，且擬合結(jié)果還可能降低原有高質(zhì)量數(shù)據(jù)本身的精度。為此，將原序列中的高質(zhì)量數(shù)據(jù)與同時刻的重構(gòu)NDVI數(shù)據(jù)提取出來分別作為參考子序列和重構(gòu)子序列，并計算兩者間相關(guān)系數(shù)(Correlation coefficient，CC)、均方根誤差(Root mean square error，RMSE)、平均絕對誤差(Mean absolute error，MAE)、平均相對誤差(Mean relative error，MRE)和效率系數(shù)(Coefficient of efficiency，CE)，來定量評價算法重構(gòu)效果。評價指標計算公式如下：

(5)

(6)

(7)

(8)

(9)

本文Hyb-F算法基于Matlab軟件平臺編程實現(xiàn)，4種對照方法的濾波結(jié)果通過TIMESAT軟件獲得。對照算法各點初始權(quán)重設(shè)為1，移動窗口寬度設(shè)為7，上包絡(luò)線擬合強度設(shè)為2。由Landsat影像直接算得的NDVI序列存在超出理論值域等與實際不符的無效點，因此，采用對照方法濾波前將移動窗口長度設(shè)為5，通過設(shè)定曲線各點取值、變化量、變化量與變化方向的關(guān)系、雪蓋指數(shù)等指標的閾值，識別并剔除移動窗口中超出理論值的噪聲和部分云、雪噪聲，并以保留下來數(shù)據(jù)的二次多項式擬合曲線在噪聲點的取值代替原始值。

2 結(jié)果與討論

2.1 樣本擬合結(jié)果直觀視覺比較

由于各覆被類型的植被生長狀況有所差異，將各覆被類型分為“有植被生長覆被類型”與“非植被生長覆被類型”兩大類進行分析，從各覆被類型樣本中分別選取一個擬合結(jié)果接近的樣本和一個擬合結(jié)果差異較大樣本進行直觀比較(見圖3和圖4)。

圖3(a)和(b)分別為蘆葦沼澤樣本各算法擬合結(jié)果接近和差異較大的曲線，(c)和(d)分別為旱田樣本各算法擬合結(jié)果接近和差異較大的曲線，(e)和(f)分別為喬木林樣本各算法擬合結(jié)果接近和差異較大的曲線，(g)和(h)分別為灌木林樣本各算法擬合結(jié)果接近和差異較大的曲線，(i)和(j)分別為草地樣本各算法擬合結(jié)果接近和差異較大的曲線。由圖3可知，有植被生長覆被類型，植被均呈現(xiàn)一年一峰的生長規(guī)律，峰值均出現(xiàn)在8月。草地樣本的峰值較低，介于0.4～0.5之間，灌木林地樣本峰值稍高，介于0.5～0.6之間，蘆葦沼澤、喬木林地和旱田樣本峰值最高，達到0.7以上。Hyb-F、S-G和AG算法都能平滑NDVI曲線。Hyb-F算法整體上與AG算法和S-G算法的三次迭代運算結(jié)果更接近。該算法的異常值判斷方案分3步剔除噪聲，植被NDVI曲線峰值保存較好，峰值處噪聲對最終擬合結(jié)果影響較??；剔除異常值后的時間序列各點取值符合植被生長變化規(guī)律，基于此算法得到的重構(gòu)結(jié)果能體現(xiàn)植被的年內(nèi)整體生長變化特征，即全局特征。Hyb-F濾波法最終濾波過程由二次S-G濾波完成，能較好地保留NDVI曲線的細部變化特征。

無迭代過程直接得到的S-G和AG濾波結(jié)果受噪聲影響大，尤其是當(dāng)原始曲線峰值附近存在低值噪聲污染時，擬合曲線的峰值偏低，存在明顯峰值損失。兩種算法的迭代方式增加了擬合結(jié)果與原始序列的比較過程，通過調(diào)整曲線低值位置的擬合權(quán)重來調(diào)整重構(gòu)序列取值。迭代運算能夠優(yōu)化含有低值噪聲序列的擬合結(jié)果，但對高值噪聲不敏感，迭代結(jié)果趨近時間序列的上包絡(luò)線。此外，重設(shè)權(quán)重時雖識別了低值噪聲，但低值仍然參與運算過程，影響擬合曲線質(zhì)量。

圖4(a)和(b)分別為鹽堿地樣本各算法擬合結(jié)果較為接近和差異較大的曲線，(c)和(d)分別為城鎮(zhèn)工礦建設(shè)用地樣本各算法擬合結(jié)果較為接近和差異較大的曲線，(e)和(f)分別水域樣本各算法擬合結(jié)果較為接近和差異較大的曲線。由圖4可知，研究區(qū)內(nèi)重度鹽堿地、工礦建設(shè)用地鮮有植被生長，年內(nèi)NDVI取值接近于0，且無明顯變化趨勢。重構(gòu)后的水域樣本NDVI時間序列取值普遍為負值，整體呈現(xiàn)先降低再升高趨勢。冬季曲線取值偏高，可能東北地區(qū)冬季水體表面凍結(jié)引起的水表反射率變化導(dǎo)致。夏季水體NDVI取值出現(xiàn)低矮波峰，可能由水面植被生長造成[17]。對于重度鹽堿地和城鎮(zhèn)工礦建設(shè)用地的理想擬合結(jié)果應(yīng)為接近于0的恒定值，然而在中等分辨率影像上，城鎮(zhèn)內(nèi)部綠化區(qū)域與建筑物構(gòu)成混合像元，鹽堿地表層生長的少量植被與出露鹽堿地構(gòu)成混合像元，使城鎮(zhèn)工礦用地和鹽堿地在影像分辨率尺度上很難達到絕對純凈，NDVI曲線受植被影響，但總體應(yīng)保持原有趨勢。對于非植被生長覆被類型(圖4)，Hyb-F組合濾波法重構(gòu)的時間序列，最大值降低，最小值提高，在一定程度上削弱了曲線體現(xiàn)出的植被生長規(guī)律，達到平滑曲線的效果，但仍存在細小波動。AG和S-G三次迭代濾波算法所得結(jié)果仍然趨近曲線整體上包絡(luò)線，當(dāng)高值突變噪聲存在時，其對鹽堿地(圖4(a)(b))、城鎮(zhèn)工礦用地(圖4(c)(d))的部分樣本的濾波結(jié)果表現(xiàn)出植被生長的變化規(guī)律，與實際情況不符。對于此類型的NDVI曲線擬合，Hyb-F和AG的直接擬合結(jié)果更接近真實情況，AG算法直接擬合效果更佳。

圖3 植被覆蓋區(qū)域不同重構(gòu)算法視覺效果對比.Fig.3 Visual comparison of fitted curve of different methods for each land cover type with vegetation

圖4 非植被覆蓋區(qū)域不同重構(gòu)算法視覺效果對比Fig.4 Visual comparison of fitted curve of different methods for each land cover type without vegetation

2.2 樣本點重構(gòu)序列指標定量分析

采用不同濾波算法重構(gòu)各類地表覆被時間序列，重構(gòu)結(jié)果的精度評價指標如表2所示。表2中前5類覆被類型為有植被生長的覆被類型，后類為非植被生長覆被類型。AG1和SG1分別表示AG、S-G算法無迭代的直接擬合方式，AG3和SG3分別表示兩種算法三次迭代擬合方式，Hyb-F表示本文提出組合濾波算法。表中粗體數(shù)字表示各地類擬合結(jié)果對應(yīng)統(tǒng)計指標的最佳取值，帶下劃線的數(shù)據(jù)表示各地類擬合結(jié)果對應(yīng)統(tǒng)計指標的最差取值。

Hyb-F算法對灌木林地和草地的重構(gòu)結(jié)果子序列與原子序列的相關(guān)系數(shù)分別為0.9215和0.9158，與其他種擬合結(jié)果相比，RMSE、MAE和MRE的取值最小。Hyb-F算法對灌木林地和草地樣本的重構(gòu)效果最好。AG算法的直接擬合結(jié)果對這兩種覆被類型的擬合結(jié)果精度較差。Hyb-F算法對旱田樣本的重構(gòu)結(jié)果與其他算法相比相關(guān)系數(shù)較小，但也達到0.9108，說明Hyb-F重構(gòu)算法也能構(gòu)建高質(zhì)量的旱田NDVI時間序列。在重構(gòu)旱田數(shù)據(jù)時，AG算法的效果最好。Hyb-F算法對蘆葦沼澤的重構(gòu)效果最差，與原子序列的相關(guān)系數(shù)為0.8036，在幾種算法中相關(guān)性最低，RMSE為0.1732，MRE為0.4597，在幾種算法中重構(gòu)誤差最大。然而，蘆葦沼澤的NDVI取值除受植被類型的影響外，還受到水文條件的制約，相比于旱地，NDVI曲線變化規(guī)律較復(fù)雜，這可能是造成擬合結(jié)果偏差較大的原因。Hyb-F算法對喬木林地的重構(gòu)結(jié)果與原子序列的相關(guān)系數(shù)較低，僅高于SG3方式的重構(gòu)結(jié)果，然而，重構(gòu)子序列和參考子序列間的偏差較小，推斷Hyb-F算法對喬木林地的重構(gòu)子序列能夠還原原序列中高質(zhì)量數(shù)據(jù)，較低的相關(guān)性可能由兩子序列同時間點數(shù)據(jù)差值的方向變化導(dǎo)致。以上結(jié)果表明Hyb-F算法能有效模擬植被覆蓋區(qū)域的NDVI曲線變化趨勢。

對于鮮有植被覆蓋的重度鹽堿地和城鎮(zhèn)工礦建設(shè)用地樣本，組合濾波法的指標值優(yōu)于其他種方式，而對于水域樣本的模擬，Hyb-F 組合濾波法所得結(jié)果精度最差。非植被生長的覆被類型，時序曲線的理想值趨于常數(shù)，不能僅通過相關(guān)系數(shù)、絕對誤差和相對誤差判斷濾波效果，在分析該類地物時序曲線重構(gòu)效果時，還需引入完整時序曲線最大值、最小值、均值指標，各指標取值如表3所示。其中Original為濾波前的原始時間序列。

表2 各覆被類型不同濾波算法結(jié)果統(tǒng)計指標對比Table 2 Statistic evaluation of different methods for each land cover type

由表3可知，所有的濾波結(jié)果都能在一定程度上降低序列最大值，提高序列最小值，使NDVI時序數(shù)據(jù)取值跨度減小，更趨近于常數(shù)，達到平滑曲線效果。Hyb-F濾波算法對水域的平滑效果最好，對重度鹽堿地和工礦用地的平滑效果僅次于AG1方法。采用Hyb-F算法能有效重構(gòu)非植被生長的覆被類型的時序數(shù)據(jù)。

表3 非植被生長覆被類型濾波算法結(jié)果額外統(tǒng)計指標對比Table 3 Additional statistic evaluation of different methods for each land cover type without vegetation

2.3 區(qū)域應(yīng)用

2.3.1 區(qū)域應(yīng)用視覺效果分析

遙感影像云噪聲能通過人眼直接識別，圖5為向海自然保護區(qū)云污染數(shù)據(jù)采用Hyb-F濾波前后對比圖，圖5(a)和(b)為4月29日厚云污染數(shù)據(jù)濾波前后對比圖， (c)和(d)為9月19日絮狀云污染數(shù)據(jù)濾波前后對比圖。從圖5可見，以植被指數(shù)時間序列影像為基礎(chǔ)，采用Hyb-F算法重構(gòu)時序數(shù)據(jù)，能識別云噪聲點，并擬合噪聲點處數(shù)據(jù)取值，濾除影像云層。由圖5 (c)和(d)可知，在無云遮擋區(qū)域，組合濾波算法能有效保持原有數(shù)據(jù)的取值水平。

2.3.2 區(qū)域應(yīng)用指標定量分析

為定量分析Hyb-F算法區(qū)域擬合結(jié)果，將一景噪聲污染影像的重構(gòu)結(jié)果與獲取條件接近的一景高質(zhì)量參考影像進行對比，計算定量評價統(tǒng)計指標。由于植被生長速率有限，對同一樣點，一般假定獲取時間間隔在5天內(nèi)的NDVI值近似相等。Landsat衛(wèi)星傳感器相對穩(wěn)定，且研究區(qū)為自然保護區(qū)，較少人為干擾，相鄰年份同時期土地利用類型相對固定，氣候條件穩(wěn)定。選取2015年5月13日與2016年5月15日拍攝影像評價算法區(qū)域應(yīng)用效果，影像特征見表1。由影像拍攝時間的歷史氣象數(shù)據(jù)得知，兩景影像都是晴天獲取，且獲取影像時段的氣溫接近，符合假設(shè)條件，可近似認為兩景影像獲取的NDVI取值相等。由于水體具有流動性，NDVI取值不確定度大，變化迅速，因此在比較區(qū)域應(yīng)用指標時，僅計算非水體部分的擬合結(jié)果統(tǒng)計指標(見表4)。

將組合算法在整個向海自然保護區(qū)應(yīng)用，濾波結(jié)果與臨近參考影像間的相關(guān)系數(shù)為0.8350，數(shù)據(jù)重構(gòu)效率系數(shù)為0.5904，組合濾波算法能夠有效模擬整個區(qū)域土地覆被類型的NDVI實際情況。

圖5 研究區(qū)云污染數(shù)據(jù)濾波前后對比.Fig.5 Comparison of imagery before and after the regional application of the hybrid filtering algorithm

表4 Hyb-F算法區(qū)域應(yīng)用效果統(tǒng)計指標Table 4 Statistic evaluation for the regional application of the hybrid filtering algorithm

3 結(jié) 論

利用Landsat8影像構(gòu)建NDVI時間序列，并采用Hyb-F組合濾波方案剔除序列噪聲，重構(gòu)NDVI時間序列曲線。與常用的S-G局部濾波算法與AG半全局擬合算法相比，Hyb-F算法有以下特點：

(1)能更有效地全面識別序列全局和局部異常點，提高曲線擬合精度，減少低值異常點引起的峰值損失。

(2)Hyb-F算法兼顧時間序列曲線的局部特征和全局變化趨勢，對于有植被生長的覆被類型，濾波結(jié)果與原始數(shù)據(jù)相關(guān)性較大，樣本的重構(gòu)結(jié)果子序列與純凈數(shù)據(jù)子序列間的相關(guān)度達0.8以上，對草地和灌木林地樣本的重構(gòu)精度最高。由于水文因素影響，蘆葦沼澤NDVI的變化較為復(fù)雜，組合式濾波法的擬合結(jié)果精度較差。

(3)對于非植被生長覆被類型，該算法能有效濾除序列的波動信息，平滑時序曲線。

本研究主要的數(shù)據(jù)源為Landsat8 OLI多光譜影像，在今后的研究中可基于其他類型數(shù)據(jù)繼續(xù)實驗，收集研究區(qū)植被物候觀測數(shù)據(jù)，將物候信息提取能力引入到算法精度評價方案中，并盡量獲取野外實測的NDVI值，得到更精確、更全的面驗證結(jié)果。

參考文獻：

[1] Zhou J, Jia L, Menenti M, et al. On the performance of remote sensing time series reconstruction methods—A spatial comparison[J]. Remote Sensing of Environment, 2016, 187:367-384.

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

[3] 劉亞南, 肖飛, 杜耘. Logistic函數(shù)方法擬合多時序NDVI數(shù)據(jù)的改進研究[J]. 遙感技術(shù)與應(yīng)用, 2015, 30(4):737-743.

Liu Ya-nan, Xiao Fei, Du Yun. Improved logistic model for fitting time-series NDVI data[J]. Remote Sensing Technology and Application, 2015, 30(4):737-743.

[4] Jonsson P, Eklundh L. Seasonality extraction by function fitting to time-series of satellite sensor data[J]. IEEE Transactions on Geoscience & Remote Sensing, 2002, 40(8):1824-1832.

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

[6] Zhou J, Jia L, Menenti M. Reconstruction of global MODIS NDVI time series: Performance of harmonic analysis of time series (HANTS)[J]. Remote Sensing of Environment, 2015, 163:217-228.

[8] 張霞, 李儒, 岳躍民, 等. 諧波改進的植被指數(shù)時間序列重建算法[J]. 遙感學(xué)報, 2010, 14(3): 437-447.

Zhang Xia, Li Ru, Yue Yue-min, et al. Improved algorithm for reconstructing vegetation index image time series based on fourier harmonic analysis[J].Journal of Remote Sensing, 2010, 14(3):437-447.

[9] 劉珺, 田慶久, 黃彥, 等. 黃淮海夏玉米物候期遙感監(jiān)測研究[J]. 遙感信息, 2013, 28(3):85-90.

Liu Jun, Tian Qing-jiu, Huang Yan, et al. Monitoring phenology of summer maize in huang-huai-hai using remote sensing technology[J]. Remote Sensing Information, 2013, 28(3):85-90.

[10] Liang S, Ma W, Sui X, et, al. Extracting the spatiotemporal pattern of cropping systems from NDVI time series using a combination of the spline and HANTS algorithms: a case study for shandong province[J].Canadian Journal of Remote Sensing, 2017,43(1):1-15.

[11] 卞建民, 林年豐, 湯潔. 吉林西部向海濕地環(huán)境退化及驅(qū)動機制研究[J]. 吉林大學(xué)學(xué)報:地學(xué)版, 2004, 34(3):441-444,458.

Bian Jian-min, Lin Nian-feng, Tang Jie. The environment degradation and its driving mechanism of Xianghai wetland in the west of Jilin province[J]. Journal of Jilin University (Earth Science Edition), 2004, 34(3):441-444, 458.

[12] 李榮平, 周廣勝, 郭春明, 等. 1981-2005年中國東北榆樹物候變化特征及模擬研究[J]. 氣象與環(huán)境學(xué)報, 2008, 24(5):20-24.

Li Rong-ping, Zhou Guang-sheng, Guo Chun-ming, et al. Phenology characteristics of ulmus pumila and its simulation from 1981 to 2005 in Northeast China[J]. Journal of Meteorology and Environment, 2008, 24(5):20-24.

[13] 唐蘊, 王浩, 嚴登華. 向海自然保護區(qū)濕地生態(tài)需水研究[J]. 資源科學(xué), 2005, 27(5):101-106.

Tang Yun, Wang Hao, Yan Deng-hua. Wetland ecological water demand in Xianghai natural reserve[J]. Resources Science, 2005, 27(5):101-106.

[14] 周增光, 唐娉. 基于質(zhì)量權(quán)重的Savitzky-Golay時間序列濾波方法[J]. 遙感技術(shù)與應(yīng)用, 2013, 28(2):232-239.

Zhou Zeng-guang, Tang Pin. VI-quality-based savitzky-golay method for filtering time series data[J]. Remote Sensing Technology and Application, 2013, 28(2):232-239.

[15] Wu M, Zhang X, Huang W, et al. Reconstruction of daily 30 m data from HJ CCD, GF-1 WFV, landsat, and MODIS data for crop monitoring[J]. Remote Sensing, 2015, 7(12):16293-16314.

[16] Li X, Shen H, Li H, et al. Patch matching-based multitemporal group sparse representation for the missing information reconstruction of remote-sensing images[J]. IEEE Journal of Selected Topics in Applied Earth Observations & Remote Sensing, 2016, 9(8):3629-3641.

[17] 李天祺, 朱秀芳, 潘耀忠, 等.環(huán)境星NDVI時間序列重構(gòu)方法研究[J]. 遙感信息, 2015(1):58-65.

Li Tian-qi, Zhu Xiu-fang, Pan Yao-zhong, et al. NDVI time-series reconstruction methods of China's HJ satellite imagery[J]. Remote Sensing Information, 2015(1): 58-65.