基于Hyb-F組合濾波算法的向海自然保護區NDVI時間序列重構

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

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

0 引 言

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

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

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

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

1 材料與方法

1.1 研究區地理概況

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

圖1 向海自然保護區Fig.1 Xianghai natural reserve

1.2 數據來源

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

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

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

1.3 研究方法

1.3.1 設定指標閾值剔除噪聲

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

1.3.2 Grubbs異常點檢驗

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

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

1.3.3 Savitzky-Golay濾波算法

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

(1)

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

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

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

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

1.3.4 非對稱高斯擬合算法

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

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

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

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

(2)

式中：高斯函數的具體形式為：

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

(3)

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

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

F(t)=

(4)

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

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

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

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

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

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

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

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

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

(5)

(6)

(7)

(8)

(9)

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

2 結果與討論

2.1 樣本擬合結果直觀視覺比較

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

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

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

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

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

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

2.2 樣本點重構序列指標定量分析

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

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

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

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

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

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

2.3 區域應用

2.3.1 區域應用視覺效果分析

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

2.3.2 區域應用指標定量分析

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

將組合算法在整個向海自然保護區應用，濾波結果與臨近參考影像間的相關系數為0.8350，數據重構效率系數為0.5904，組合濾波算法能夠有效模擬整個區域土地覆被類型的NDVI實際情況。

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

表4 Hyb-F算法區域應用效果統計指標Table 4 Statistic evaluation for the regional application of the hybrid filtering algorithm

3 結 論

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

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

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

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

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

參考文獻：

[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函數方法擬合多時序NDVI數據的改進研究[J]. 遙感技術與應用, 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] 張霞, 李儒, 岳躍民, 等. 諧波改進的植被指數時間序列重建算法[J]. 遙感學報, 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] 劉珺, 田慶久, 黃彥, 等. 黃淮海夏玉米物候期遙感監測研究[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] 卞建民, 林年豐, 湯潔. 吉林西部向海濕地環境退化及驅動機制研究[J]. 吉林大學學報:地學版, 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]. 氣象與環境學報, 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] 唐蘊, 王浩, 嚴登華. 向海自然保護區濕地生態需水研究[J]. 資源科學, 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] 周增光, 唐娉. 基于質量權重的Savitzky-Golay時間序列濾波方法[J]. 遙感技術與應用, 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] 李天祺, 朱秀芳, 潘耀忠, 等.環境星NDVI時間序列重構方法研究[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.