植被遥感时间序列数据重建方法简述及示例分析

李伟光， 侯美亭

(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以外的其他植被遥感时间序列。

总体上，本研究主要侧重于回顾近年来的植被遥感时间序列重建方法，但仅采用了模拟数据和一个样点的遥感数据进行时间序列重建方法比较。而随区域的不同，插值方法的表现可能有所差异，故研究结果还有待进一步推广试验。