植被遙感時間序列數據重建方法簡述及示例分析

李偉光， 侯美亭

(1.海南省氣候中心，海口 570203； 2.海南省南海氣象防災減災重點實驗室，海口 570203； 3.中國氣象局氣象干部培訓學院，北京 100081)

0 引言

遙感生態監測起步于20世紀70年代，從早期利用陸地衛星(Landsat)、氣象衛星(如AVHRR)監測地表植被格局[1-2]，到近年來新興的使用衛星遙感太陽誘導葉綠素熒光(solar-induced chlorophyll fluorescence，SIF)監測植物光合作用和評估生產力[3]、利用無人機等新型平臺監測地表生態狀況等，基于遙感監測地表生態已經發展成為地球生態系統監測的重要組成部分。遙感具有基于現場原位生態觀測不可比擬的空間覆蓋優勢，遙感方法是在局地和全球尺度上對生態系統的短期和長期狀態、壓力、干擾和資源限制進行及時、經濟有效、可重復記錄和評估的唯一方法[4]。

衛星時間序列數據已被證明對闡明生態系統動態變化和驅動因素特別有價值[5-6]。Cord等[7]匯總了衛星遙感應用于生態系統監測的服務類別，指出食物供應和氣候服務(溫室氣體調節)是衛星遙感應用最為廣泛的2大生態服務領域，土地利用/土地覆蓋(land use and land cover，LULC)和歸一化植被指數(normalized difference vegetation index,NDVI)是最為常用的2類衛星遙感產品。NDVI與綠葉密度、凈初級生產力和CO2通量高度相關，是地表生物量在景觀水平上的一個代表[8]。盡管隨遙感技術的發展，空間、光譜和時間分辨率不斷增強，對植被空間異質性和連續性的理解不斷加深，但基于遙感對離散的植被單元進行分類一直是LULC分類方法的重要基礎[9]。遙感反演的植被產品是遙感應用于生態監測領域重要的數據之一。

伴隨著遙感數據的不斷積累，植被遙感產品形成了完善的時間序列數據，利用這些數據分析生態系統變化，例如評估生態系統過程中的擾動、檢測植物物候學的時空變化、分析氣候變化對植被的影響、監測LULC變化、植被時空異質性評價等，基本上構成了近幾十年來最為活躍的生態遙感研究主題。植被覆蓋和生產力的年際變化被認為是理解全球氣候變化對生態系統動力學影響的關鍵[10]。目前，研究者已經開發出許多軟件或程序來處理遙感時間序列數據，例如TIMESAT[11]，BFAST[12]，Timestats[13]，SPIRITS[14]，DBEST[15]，phenor[16]，FORCE[17]，BEAST[18]和DATimeS[19]等。

然而，云、氣溶膠和儀器誤差等因素往往造成植被遙感時間序列數據存在不同程度的缺失或質量變差等情況。無論利用哪種軟件分析植被遙感時間序列數據，確保所使用遙感時間序列數據完整、高質量是至關重要的，高質量的完整時間序列數據才能保證結果的可信度。對存在數據缺失的時間序列數據進行時空重建是準確提取時間序列變化特征的重要前提。盡管國內外有許多研究對比了不同方法重建時間序列數據的優劣[20-25]，但新的研究成果不斷涌現，本研究試圖對以往研究使用的時間序列重建方法進行綜述，并側重于近年來的最新研究成果，通過一些示例，介紹有關新方法的適用性和潛在的應用價值。

1 植被遙感時間序列數據重建方法

衛星數據時間序列可以看作是一種3-D數組，對于這種時間序列，最常見的問題就是數據缺失，也可稱為數據異常。數據缺失問題是遙感數據預處理中面臨的最主要問題之一。缺失數據對定量遙感研究的影響會很嚴重，導致參數估計有偏差、信息丟失、可信度下降等。導致數據缺失/異常的原因可以歸納為氣象條件(云、雪、灰塵和氣溶膠)、儀器誤差/傳感器退化、數據傳輸過程中圖像數據丟失、時間分辨率低等。其中，傳感器退化引起的數據異常往往導致時間序列存在某種趨勢，對于這種因素往往無法通過常見的數據重建方法進行校正，需要數據生產方對傳感器進行評估和校準，比如MODIS C6版本數據(較C5 版本)對Terra/MODIS退化的校正[26]。本研究中的數據重建方法主要針對傳感器退化因素以外引起的數據缺失和異常。

數據缺失/異常產生的噪聲包含了不真實的時間變異，對于這類噪聲，以往研究大致有以下2類處理方法： ①利用平滑方法，對包含噪聲的原始時間序列直接進行平滑處理。例如，Mondal等[27]使用基于離散傅里葉變換(discrete Fourier transform，DFT)的平滑方法來處理原始的MOD13Q1數據，減少NDVI時間序列中的噪聲，相比其他平滑技術，DFT需要的模型參數更少。②只使用原始時間序列中的高質量數據，通過數據插值和平滑方法，重建時間序列。原始數據的質量可基于質量標記產品來判定，或者基于相應數據點的標準偏差σ，一般認為位于平均值±3σ以外的數據點是低質量的[28-30]。例如，Jeganathan等[29]將任一像元1 a內特定時段的GIMMS NDVI值與其歷年相應時段的多年平均值進行比較，然后移除NDVI值在平均值±3σ以外的那些數據點，經過這些處理后，剩余的數據點可認為是“高質量”的，再進一步利用插值方法來填補上述低質量數據消除過程中產生的數據空白。與第一類方法相比，第二類方法在數據重建時首先直接去除了異常值，相對更為可靠。本研究將著重概述第二類方法有關的研究進展。對于該類遙感時間序列數據重構方法，簡而言之，包括插值和平滑2步，即通過利用遙感數據的時空相關性對已有的存在數據缺失的高質量數據序列進行插值、平滑濾波，以達到重建時間序列的目的。

1.1 數據插值

插值是指對存在數據缺失的時間序列進行填補。一些平滑方法對于不等間距的序列可以自動插值，如Whittaker方法[31]、傅里葉方法[19]，這些方法兼顧了插值和平滑的功能。另一些平滑方法(如Savitzky-Golay)的前提就是時間序列必須連續、完整，因此在使用這些平滑方法前需進行數據缺失的填補。插值方法大致可分為基于時間的插值方法和基于空間的插值方法2種類型。

基于時間的插值方法主要是考慮到所插值的數據是與生物過程時間動態直接聯系在一起的[32]，應用尤為廣泛。比較簡單的基于時間的數據填補方法，可以由數據缺失像元的歷史同期平均或最近的有效數據(特別是對于受到積雪影響的植被數據)替代[33]，或者使用2個相鄰有效值的移動平均值替換可能由云和大氣氣溶膠等因素引起的局部不規則變化和短間隙[34-35]。更為常用的時間插值方法為線性插值[36]和樣條插值[37]，其中，樣條法也可以用于空間插值[38]和數據序列平滑[39]。與線性插值方法相比，時間序列諧波分析(harmonic analysis of time series，HANTS)、奇異譜分析(singular spectrum analysis，SSA)是更為復雜的時間插值方法。HANTS算法是快速傅里葉變換(fast Fourier transform，FFT)的一種改進，采用基于諧波分量的最小二乘曲線擬合方法，在植被遙感時間序列插值中展示了良好的性能[40-41]。利用SSA進行插值是一種近年來發展的數據插值方法，其原理在于使用迭代方法對缺失數據進行估計，然后使用這些估計來計算滯后協方差矩陣及其經驗正交函數分解分量，再使用交叉驗證來優化窗口寬度和SSA主模態的數量，以填補數據空白。盡管SSA插值方法在植被遙感數據中應用不多[42]，但在海溫、降水和土壤呼吸等多種類型的地學數據插值中顯示了其價值[43-44]。然而，類似SSA和HANTS這樣的方法是基于數據序列周期性的，如果數據缺失持續的時間長，就不能成功識別周期及振蕩，從而導致利用這些方法得到的結果缺乏可信度[45]。

基于空間的插值方法主要是利用存在數據缺失的像元周邊的(不存在數據異常的)像元信息進行插值，除樣條法以外，常用的方法還包括最鄰近插值、反距離加權插值、地統計方法(如克里格法)等。需要注意的是，數據缺失像元通常存在空間聚集的情況，在進行空間插值時，有必要評價空間自相關效應[46]。而且，由于混合像素的存在，植被指數在短距離內變化可能很大，導致基于空間的插值可能無法代表真實景觀的復雜性。

盡管目前有很多種時間和空間插值方法，但這2類方法都不能在所有情況下超越對方，而時間和空間信息相結合的插值方法正在顯現出更多的優勢。Borak等[47]以MODIS LAI時間序列為例，對幾種時間和空間插值技術進行了評估，發現插值技術的精度取決于下墊面覆蓋，時間插值在非森林覆蓋類型上執行得更有效，空間插值在森林區域執行得更好，而空間和時間方法的結合提供了比任何單一的時間和空間方法更好的插值能力。Pede等[38]發現，對于受到云覆蓋影響的MODIS 8 d LST數據，空間插值方法、空間和時間相結合的插值方法比時間插值方法顯示出了更大的優勢，而且空間和時間相結合的方法受地形、氣候等環境因素的影響更小。近年來更多研究使用來自鄰近像素的歷史數據和時間曲線來結合空間和時間信息進行插值，也顯現出了比僅用時間或空間插值更大的優勢[48-49]。Yan等[50]提出了一種應用于高空間分辨率、低時間分辨率的Landsat數據進行間隙填補的時空相似替代算法，通過合并空間和時間上相鄰的相似像素來分割圖像時間序列，然后將這些片段分組到不同的類中，需要填充的像素位置與屬于同一類的其他像素進行比較，以識別替代像素，不過其算法不考慮地表變化趨勢，因此較適用于在地表沒有發生變化期間獲取的圖像。

最近還涌現了一些基于機器學習的插值算法，如神經網絡、隨機森林和決策樹等，這些方法不同于傳統的參數方法(如線性方法)，信息來自于訓練數據，不需要假設它們的統計分布或變量之間的相互關系，但這些方法的效果如何，還有待深入評價[19]。總體上，插值的效果不僅依賴于插值算法，更依賴于原始數據的質量，當原始數據序列具有較大的高質量數據比例時(比如大于85%)，普通的線性插值一般就可以滿足需要[51]。而且，對于常用的GIMMS NDVI數據，Burrell等[52]認為當原始時間序列的有效值大于85%時，質量即滿足需求，可以直接使用原始序列。另外，在一些情形下(如物候特征提取)，當某像元缺值數據超過一定數量時，比如連續缺值超過1個月，Pastor-Guzman等[53]選擇不使用相應位置的數據，而通過插值方法獲得相應缺失的數據。

1.2 時間序列平滑

由于插值算法本身局限性以及原始時間序列數據缺失程度的影響，插值后的時間序列數據可能仍然包含噪聲，這時候需要進行平滑，平滑的主要目的是降噪，以保證重建時間序列的質量。

簡單的平滑方法包括最大值合成(maximum value composite，MVC)方法、通過使用預定義的閾值移除異常值的最佳指數斜率提取算法(best index slope extraction，BISE)[54]、迭代插值方法(iterative interpolation for data reconstruction，IDR)[32]。Holben[55]是最早提出使用MVC在時間序列中抑制NDVI噪聲的研究之一，而且MVC產品還減少了NDVI時間序列中云污染像素造成的數據間隙。然而，即使在MVC過程之后，特別是當云層覆蓋持續時間長于合成期(如在雨季)時，一些殘留的噪聲可能仍然會出現在數據中。這時就需要使用更加穩健的平滑方法對序列進行降噪。

曲線擬合方法是目前使用較多的平滑方法[36]，如SG(Savitzky-Golay)、TSGF(temporal smoothing and gaplling)、非對稱高斯函數(asymmetric Gaussian，AG)、雙邏輯函數(double logistic，DL)和傅里葉變換等應用于局部時間窗口的曲線擬合方法[56]，以及可應用于整個時間序列的Whittaker濾波方法(一種樣條平滑方法)等。其中，AG，DL和SG方法被集成到時間序列分析軟件(TIMESAT)中，使得這3種方法更易于實現[57]。曲線擬合方法除了應用于平滑降噪外，也被用于提取植被物候變化參數，因為擬合最主要的依據就是植被時間變化特征。特別是在局部窗口內擬合衛星數據時間序列，多采用高階諧波方法或更多參數，能夠更好地捕捉植被動態，尤其適于不規則、不對稱的植被指數曲線，不過它們對局部波動和數據噪聲也更為敏感，某些情況下，可以對經過局部擬合的且仍然包含有不規則數據的序列，再應用局部中值濾波去除局部不規則值[34]。相比基于整個時間序列或低階諧波的方法往往得到更平滑的曲線，但也可能偏離實際植被生長軌跡，導致過度校正，抑制了實際植被波動。

每種時間序列平滑算法都有其獨特的優點和缺點，在選擇平滑算法時，不僅應考慮時間序列的應用目的、數據中噪聲的性質、時間序列曲線的典型形狀和使用方法的方便性等[56]，還應考慮原始時間序列的數據質量[58]。如果數據質量足夠高，可以選擇保持季節內變化的濾波方法，如SG濾波[22, 59]和Whittaker濾波[31, 60]； 在高質量數據較少的情況下，比如觀測頻率相對較低的Landsat和Sentinel-2數據，函數擬合(如DL)可能是實現穩健性的唯一選擇[61]。另外，好的插值算法也能間接提高平滑方法的質量。比如，Padhee等[62]基于局地(尤其是同一土地覆蓋類型)物候時空相似性的考慮，使用參考物候曲線來預測NDVI時間序列中的缺失值的方法對MODIS NDVI時間序列進行預插值，改善了云高覆蓋地區時間序列中存在較大數據間隙而導致的HANTS方法的不足。

盡管目前沒有公認的標準來比較不同的時間序列重建方法[63]，但是可以通過評估重建后的時間序列對原始無云數據的保真度和識別受云污染值的能力來判斷方法合適與否。由于地面真實數據通常不可用，Julien等[32]以GIMMS NDVI為例，提出了以下判斷標準： 重建后的值與原始觀測值的距離、以及與上包絡線的接近度。與原始觀測值的距離旨在量化對未污染數據的保真度。對于上包絡近似準則，它基于大氣污染有降低NDVI值的趨勢這一考慮，這一標準在整個時間序列中被定義為NDVI低于原始觀測值的重建觀測值的頻率(范圍在0～1)。對于該準則，較低的值對應于更接近原始時間序列的上包絡線。這2個標準對于評價時間序列重建方法具有較好的效果，特別是當大氣污染或云對遙感反演的生物物理參數有負面影響時。

2 常用的時間序列數據重建方法比較示例

這里選取線性插值、SSA，Whittaker和HANTS等4類常用的不同類型的時間序列數據重建方法，通過在模擬的NDVI時間序列、原始NDVI完整時間序列中隨機制造不同比例的數據缺失，使用上述4類方法進行數據重建，然后計算重建時間序列與原始時間序列的均方根誤差(root mean square error，RMSE)，以衡量重建序列同原始序列之間的偏差，從而判斷不同方法的優劣。

RMSE亦稱標準誤差，值越小，擬合效果越好。計算方法為：

(1)

式中：y為原始序列；z為擬合序列；N為序列長度。

2.1 模擬NDVI序列

模擬NDVI序列包含短周期的波動和長期趨勢[64]，按照下式生成為：

(2)

對生成的模擬NDVI序列隨機刪除10%，20%，30%和40%的數據，以制造不同比例的數據缺失，然后分別使用前述的4種方法進行數據重建，計算重建后的完整時間序列與原始完整序列之間的RMSE(用RMSEall表示)以及數據缺失位置處的重建序列與相應位置原始序列之間的RMSE(用RMSEgap表示)。 4種方法的效果及相應的RMSE見圖1和表1。

(a) 缺失10%(b) 缺失20%

(c) 缺失30%(d) 缺失40%

圖1 模擬NDVI序列數據重建效果對比

表1 基于模擬NDVI序列的不同數據重建方法的RMSE對比

分析可以發現，線性插值對已有數據不進行改動，僅對數據缺失位置進行線性內插。當數據缺失較少或者數據缺失不在峰、谷端點時，線性插值具有良好的準確性(表1)。SSA方法作為基于識別序列周期性的方法，插值的序列有時具有相對較大的振幅(圖1)，具體原因可能和模擬序列的生成方法有關。Whittaker方法對原始序列的保真及數據缺失位置的插值效果最好，而HANTS方法次之。

2.2 原始NDVI時間序列

本序列來自GIMMS NDVI數據集，首先對GIMMS NDVI數據進行月最大值合成，然后篩選1982年1月—2018年12月NDVI月值質量標示全部為0(flag=0)的格點得到，所選取的序列位于經緯度分別為86.416°E，47.167°N的位置。在該序列中隨即剔除10%，20%，30%，40%數據制造數據缺失(圖2和圖3)，評估插值后數據與真實值之間的RMSE(表2)。

(a) 缺失10%

(b) 缺失20%

(c) 缺失30%

(d) 缺失40%

(a) 缺失10%(b) 缺失20%

(c) 缺失30%(d) 缺失40%

圖3 NDVI序列數據重建方法效果局部放大對比

表2 基于GIMMS NDVI時間序列的不同數據重建方法的RMSE

圖2和表2顯示了針對GIMMS NDVI序列不同數據缺失比例下的4種數據重建方法的性能，圖3截取了其中一段展示其局部效果。可以發現，線性方法保持了模擬序列時的表現，隨著數據缺失比例增加，RMSEgap逐漸上升。而處理模擬NDVI序列時表現最好的Whittaker方法依然保持了良好的性能，而且其只需要通過設置唯一參數即可平衡擬合的保真度和粗糙度，在保證擬合精度的同時，也提高了曲線的平滑度，不過，該方法的RMSEgap有時稍高于SSA。相比，在該示例中，SSA成為最好的數據重建方法，而且隨缺失數據的增多，誤差變化不大。HANTS方法對原始序列的保真和斷點位置的插值基本穩定，但該方法參數設置較為復雜，需要多次試驗才能確定合適的參數，提高平滑度的同時也犧牲了擬合的精度，而且其RMSE在不同缺失比例下都要高于SSA，這也與Ghafarian Malamiri等[42]的研究相一致。

3 結論

完整的高質量植被遙感時間序列是利用遙感準確提取植被生態系統變化特征的重要前提。然而，云覆蓋等因素嚴重制約著植被遙感產品的觀測質量，造成時間序列中存在不同程度的數據缺失或低質量數據。本研究回顧了以往常用的時間序列重建方法，這些方法可以概括為2大類型： ①利用平滑方法對包含噪聲的原始時間序列直接進行平滑處理； ②僅使用原始時間序列中的高質量數據，通過數據插值和平滑方法，重建完整時間序列。后者是目前數據重建的主流方法，且已經發展了相應的多種多樣的數據插值和平滑方法。

在回顧以往數據重建方法的基礎上，本研究以在模擬的NDVI時間序列、原始NDVI完整時間序列中隨機制造不同比例的數據缺失為例，具體對比了線性插值、SSA方法、Whittaker方法和HANTS方法等4類方法的數據重建性能，研究發現，這些方法在NDVI數據缺失插補中均取得了較好的精度。然而，它們有各自的優點和缺點，在不同的情況下填補數據空缺時需要加以考慮。其中，線性插值是一種廣泛使用的方法，然而當缺失數據過多時，這種方法并不適用； 相比Whittaker 方法、SSA方法和HANTS方法顯示出了更好的性能，盡管一些情形下，SSA及HANTS對應的RMSE高于Whittaker方法，但這3種方法的整體差異并不大； 不過，HANTS方法需要設置的參數相對較多，SSA方法在某些情形下的重建序列具有較大的振幅，而Whittaker方法在不同的數據缺失情形下都顯示出了較好的效果，且需要的參數較少，是一種值得推廣的插值方法。當然，除NDVI以外，EVI，EVI2以及近年來的SIF等也是較為常用的植被遙感產品。由于不同植被遙感時間序列最重要的共同特征之一是具有顯著的季節動態，而常見的植被遙感時間序列重建方法一般具有很好的普適性(比如Whittaker 方法)，也多基于序列周期性(比如HANTS方法)，因此它們也適用于NDVI以外的其他植被遙感時間序列。

總體上，本研究主要側重于回顧近年來的植被遙感時間序列重建方法，但僅采用了模擬數據和一個樣點的遙感數據進行時間序列重建方法比較。而隨區域的不同，插值方法的表現可能有所差異，故研究結果還有待進一步推廣試驗。