两种应用场景下NDVI时间序列数据拟合方法研究

邰文飞,张新胜,蔡明勇,申振,申文明,史雪威,陈绪慧

(生态环境部卫星环境应用中心，北京 100094)

植被覆盖指数(NDVI)时间序列数据(以下简称“NDVI时序数据”)能够模拟地表植被的生长状况和覆盖信息，已成为生态环境质量监测和评价的重要数据来源，广泛应用于覆被变化动态监测、生态环境变化监测与模拟、物候信息识别与提取、人类活动扰动识别等方面[1-4]。NDVI时序数据集包含大量噪声，受地面状况、大气干扰、传感器状态等因素的影响，需要根据不同的区域类型、覆被类型和应用场景，选取适用的滤波算法或拟合方法[5]。目前，NDVI时序数据拟合常用的方法有多项式最小二乘法、非对称高斯函数拟合法(AG法)、双Logistic曲线拟合法(D-L法)、Savitzky-Golay滤波法(S-G法)、时间序列谐波分析法(HANTS法)、Whittaker平滑法等。

张晗等[6]利用S-G法、Whittaker平滑法、HANTS法这3种方法拟合重建陕西省2000—2012年MOD13Q1 NDVI数据，对比3种方法在物候提取和复种指数提取中的应用。李天祺等[7]选取北京市MOD13Q1 NDVI数据和环境星多光谱数据，使用HANTS法、AG法、D-L法和S-G法4种方法拟合重建后，对比典型地物重建效果，结合该地区农作物物候站点数据，评价各种方法物候信息的一致性。关于时序数据拟合方法的对比，此前的研究主要聚焦于方法机理和参数分析上，对不同覆被类型和应用场景的研究较少，目前还没有一种公认的普适性方法。现选取农作物物候参数提取和扰动识别2种NDVI时序数据集应用场景，对比分析在这2种应用场景下，3种拟合方法对NDVI时序数据集拟合的适用性。

1 研究区概况及数据来源

1.1 研究区概况

兖州市地处山东省中部山地泰沂山区西南部的山前倾斜平原，地势整体上东北高、西南低，境内主要地形为平原，面积646.7 km2，占总面积的99.7%，耕地面积大、产量高，近年来森林覆盖率持续增高，煤炭是该市最具优势的矿产资源，分布面积约占总面积的37.06%。

韦兹县位于美国弗吉尼亚州阿巴拉契亚地区，地处西南弗吉尼亚煤田的中西部，地表大部分被森林覆盖，NDVI值总体较高，该地区因露天开采导致大量地表植被剥离，对生态环境破坏巨大，后期又开展了生态恢复，整个NDVI时序数据包含“较高水平—急剧下降—逐渐恢复”的变化规律，有利于在扰动识别应用场景中对比分析3种拟合方法的重建效果。

1.2 数据来源

1.2.1 MODIS时序数据集

采用兖州市2014年MOD13Q1数据产品，空间分辨率为250 m，16 d最大值合成[8]。采用MRT工具对影像进行投影转换、重采样和NDVI计算等处理，全年时序数据由23期NDVI产品组成。

1.2.2 Landsat TM时序数据集

采用韦兹县1996—2011年Landsat TM影像数据，数据产品级别为1 T，空间分辨率为30 m。最大限度地获取了时相差较小的影像，获取时间集中在6—8月的植被生长季。由于2006和2009年该地区数据时相与其他年份相差较大，因此没有使用这2期数据。

1.2.3 土地利用数据

采用GLC 2010土地覆盖数据集，空间分辨率为30 m，其分类详细，总体分类精度达到80%以上，因此该数据集作为本研究分类结果验证时的真实数据[9]。

1.2.4 农作物物候观测数据

采用农作物地面观测站点数据作为验证和参考数据。由于兖州市内无观测站点，因此参照距离兖州市最近的济宁市任城区地面观测站点数据(116.58°E，35.45°N)。

2 研究方法

基于TIMESAT软件，选取农作物物候参数提取和扰动识别2种NDVI时序数据集应用场景，对比分析AG法、D-L法和S-G法这3种方法的拟合效果和适用性。2种应用场景下NDVI时间序列数据拟合方法研究技术路线见图1。

图1 2种应用场景下NDVI时间序列数据拟合方法研究技术路线

兖州市使用3种方法提取农作物物候参数，包括生长开始时间(SOS)、生长结束时间(EOS)和生长周期(LOS)，对比分析3种拟合方法在物候参数提取中对NDVI时序数据集的应用场景的适用性；韦兹县选取基于样本像元评价和Jacknife评价2种评价方法，对比分析3种拟合方法在扰动识别中对NDVI时序数据集的应用场景的适用性。同时引入地面站点数据作为参照和验证数据。

2.1 拟合重建方法

2.1.1 AG法

AG法是一种基于不对称高斯函数的非线性最小二乘拟合算法，该方法的核心是先局部最优化最后进行全局拟合，其灵活性较好。在时间序列中，选取一个局部拟合区间，即最大或最小值区间，使用高斯拟合函数作为局部拟合模型，拟合这一区间的数据，最后使用全局拟合模型合并拟合局部拟合结果。采用分段拟合的思想，可以确保拟合结果更加接近当前时段的真实变化规律，可以很好地描述作物生长和交替过程中植被指数和时间的关系[10]。李儒等[11]将其拟合过程分解为3个步骤：区间提取、局部拟合和整体拟合。

局部拟合函数见式(1)：

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

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

(1)

式中：c1、c2——曲线的基准和振幅；t——时相；a1——曲线最大值或最小值的位置；a2、a3、a4、a5——分别为左、右半边曲线的宽度和峭度。

整体拟合函数见式(2)：

(2)

式中：tL、tC、tR——时间序列中尚未拟合部分的左边最大值、中间最大值、右边最大值对应的时间节点；fL(t)、fC(t)、fR(t)——分别为左边最大值、中间最大值及右边最大值对应的局部拟合函数；α(t)、β(t)——介于0和1之间的剪切系数。

2.1.2 D-L法

D-L法的原理与AG法基本一致，同样是先进行局部拟合再进行整体拟合。Beck[12]认为D-L法对植被生长季进行提取的准确率更高。与AG法相比，该方法为双逻辑形式，同时公式中少一个参数[13]。因此，运用该方法重建时序数据，结果与AG法基本吻合，只有经过大量的统计对比，该法才会产生微弱的劣势。

局部拟合函数见式(3)：

(3)

式中：a1、a2、a3、a4——分别为左、右半边曲线的拐点位置及拐点处的变化速率。

整体拟合函数参照 AG法的函数公式(2)。

2.1.3 S-G法

S-G法应用最小二乘卷积算法，通过计算一组相邻值达到数据滤波的目的[14]，计算公式见式(4)：

(4)

S-G法重建NDVI时序数据时，首先需要人为设定2个参数，即滤波窗口大小和多项式拟合阶数，合理的参数能够保证NDVI时序数据拟合的准确性。滤波窗口越小，产生越多冗余数据，不易获取数据长期变化趋势；滤波窗口越大，时序曲线就越平滑，容易遗漏一些细节变化信息。多项式拟合阶数通常选取2～4，较低阶数使曲线更加平滑，但会保留异常值；较高阶数可以去掉异常值，但会出现过度拟合，容易出现新噪声。对于不同的研究区域和应用场景，需要根据经验反复尝试，才能得到比较合理的设置参数。

2.2 农作物物候参数提取方法

以兖州市2014年23期16 d合成MODIS影像为数据源，使用TIMESAT软件提取植被物候参数。植被物候参数提取方法主要有拟合法、阈值法、最大斜率法和滑动平均法等，TIMESAT软件采用的是动态阈值法，该方法将NDVI值增长(降低)达到当年NDVI值振幅一定百分比的时刻定义为生长季的开始(结束)时间[15-16]。TIMESAT软件容许使用者自由设定阈值百分比，参考大量相关文献，将该值取为20%[17-19]，即生长季开始和结束时刻设定为NDVI值升高和降低达到振幅的20%。农作物站点数据资料显示，兖州市主要农作物是冬小麦和夏玉米，1年2季，1季小麦，1季玉米。兖州市农作物生长各阶段示意见图2。

由图2可见，整个曲线由冬小麦和夏玉米2个完整的生长周期组成。从1月1号开始，NDVI值从冬小麦返青期前最低值开始逐渐升高，抽穗期达到峰值，成熟期后逐渐降低，收获时降到最低；夏玉米播种后，出苗期之后NDVI值迅速升高，抽穗期达到第2个峰值，成熟期后逐渐降低，到11月第2个生长周期结束。

图2 兖州市农作物生长各阶段示意

2.3 扰动识别方法

2.3.1 基于样本像元评价

根据扰动引起的NDVI时序数据变化特征，结合韦兹县原始TM影像和Google Earth影像，选取采矿扰动点和建筑扰动点，使用3种方法重建后，对比重建前、后的扰动像元曲线特征变化，在分析重建后NDVI时序数据识别扰动可行性的同时，对比3种重建方法。NDVI时间序列数据重建的目的是去除噪声的同时，最大限度地保留时序曲线的真实变化特征[20]。矿区常见的扰动为采矿扰动和建筑扰动。扰动使NDVI值突降，而噪声一般也是突降点，因此，分析扰动引起突降后的时序曲线特征，重建时去除噪声，保留扰动引起的真实变化信息，从而达到识别扰动的目的。

2.3.2 Jacknife法评价

在无法获取真实的实地数据作为验证的情况下，可以创建理想的重建模型环境，实现不同重建方法效果的客观评价。使用Ma等[21]提出的Jacknife法对重建结果进行评价，核心思想是在一个时间序列中，随机加入噪声，噪声可以取原值的10%，20%甚至更大，使用多种重建方法重建此时间序列后，可以直观显示各种重建方法的重建效果。

本研究使用Jacknife法，假想1条光滑的时序曲线，改变某一像元或某些像元的NDVI值，使用上述3种方法重建时序数据，运用目视判断和标准误差法评价3种方法的拟合重建效果。结合Jacknife法，设计1条平滑且变化幅度较小的曲线，根据加入噪点的数量和方式的变化，设计a、b、c、d、e5条曲线。a曲线：设置2个不连续噪声后的曲线；b曲线：设置2个连续噪声后的曲线；c曲线：设置3个连续噪声后的曲线；d曲线：模拟建筑扰动，NDVI值骤降，然后变化幅度较小；e曲线：模拟采矿扰动，NDVI值骤降，然后缓慢升高。

2.4 对拟合精度的评价方法

2.4.1 基于统计量的拟合结果评价

选取回归估计标准差(RMSE)，定量分析3种方法的拟合重建效果。RMSE即均方根误差，表示2个样本数组值间的差异程度。现用于对比重建后NDVI年内时间序列与原始值之间的平均差异程度。RMSE大小与原始NDVI整体值大小有关，但还是能反映重建前后NDVI值的代表性强弱，其值越小，拟合值的代表性越强。计算公式见式(5)：

(5)

式中：NDVIpi、NDVIoi——时间序列中第i期拟合处理前、后的NDVI值；N——像元总数。

2.4.2 扰动识别评价

重建后的曲线对原始NDVI时序曲线扰动特征的保持度，其对扰动识别的准确性有很大影响。可通过构建时间差、振幅差和恢复速率差3个参量来表征NDVI时序曲线的扰动特征。

(1)时间差。即重建前、后扰动开始时间之差，表示扰动时间识别的准确性，计算公式见式(6)：

T=t2-t1

(6)

式中：T——时间差；t1——原始曲线扰动开始时间；t2——重建后扰动开始时间。

(2)振幅差。振幅即扰动开始时的NDVI值与扰动后NDVI最小值之差，振幅差即重建前、后振幅之差，表示扰动强度识别的准确性，计算公式见式(7)：

A=a2-a1

(7)

式中：A——振幅差；a1——原始曲线振幅；a2——重建后曲线振幅。

(3)恢复速率差。恢复速率即开始恢复到稳定水平时每年NDVI值增长量，反映在图中即2点连线的斜率，恢复速率差即重建前、后曲线斜率之差，表示恢复速率识别的准确性，计算公式见式(8)：

(8)

式中：R——恢复速率差；hmaxAF——重建后恢复的最大NDVI值；hminAF——重建后开始恢复时的NDVI值；hmaxRD——重建前恢复的最大NDVI值；hminRD——重建前开始恢复时的NDVI值；tmaxAF——重建后恢复最大NDVI值的时间；tminAF——重建后开始恢复的时间；tmaxRD——重建前恢复最大NDVI值的时间；tminRD——重建前开始恢复的时间。

重建前、后扰动特征变化示意见图3。

图3 重建前、后扰动特征变化示意

3 结果与讨论

3.1 农作物物候参数对比

选取SOS、EOS和LOS这3个农作物物候参数。采用TIMESAT软件提取参数，经过处理后，使用站点数据验证提取结果。3种方法重建效果有所不同，提取的物候参数也略有差异。同时，提取耕地和林地生长周期数，基于耕地和林地样本点，对比分析3种方法生长周期数提取结果。

重建后物候参数与站点数据对比见表1。由表1可见，冬小麦的SOS站点数据是2月20日，AG法和D-L法提取时间为2月20日，S-G法提取时间为2月15日。冬小麦和夏玉米的SOS提取结果显示，AG法和D-L法提取结果比S-G法更接近站点数据，3种方法EOS提取时间相同，3种方法冬小麦LOS提取时间均与站点数据相差较大，夏玉米LOS提取时间与站点数据基本一致。因此，在SOS提取中，AG法和D-L法优于S-G法。

表1 重建后物候参数与站点数据对比

兖州市每年耕地有2个生长周期，林地有1个生长周期，耕地、林地重建后获取的生长周期数占比见表2。

表2 耕地、林地重建后获取的生长周期占比

由表2可见，统计耕地100个样本像元，AG法提取的周期数为2的像元占总像元的95%，周期数为1的像元占总像元的5%，提取结果准确率最高的是AG法，其次是D-L法，S-G法提取误差最大。统计50个林地样本像元，对比林地生长周期数提取结果，AG 法和D-L法提取的周期数为1的像元占总像元的56%，S-G法提取的周期数为1的像元占总像元的64%，S-G法比AG法、D-L法准确率更高。综上，AG法在耕地周期数提取中准确率最高，D-L法次之；S-G法在林地耕地周期数提取中准确率最高。

3.2 扰动识别对比

3.2.1 韦兹县扰动像元拟合结果

采矿扰动识别评价结果显示，AG法、D-L法和S-G法时间差分别为0.6，0.7和0.3，振幅差分别为0.298，0.301和0.216，即S-G法重建后扰动时间和振幅与原始时序曲线差距最小。使用原始恢复速率(0.059)减去重建后恢复速率，得到恢复速率差，S-G法小于AG法和D-L法。3种方法重建后恢复速率分别与原始恢复速率做相关性分析，S-G法相关系数达到0.618，远高于前2种方法。重建前、后采矿扰动恢复速率统计见表3。

表3 重建前、后采矿扰动恢复速率统计

建筑扰动一般为永久性扰动，在NDVI时序曲线上，表现为NDVI值突降之后，随时间起伏变化较小。重建前、后建筑扰动像元时序曲线见图4。由图4可见，NDVI值从0.65突降至2000年前后的0.1，2000—2011年，NDVI值在0.25上下变化。对比建筑像元重建前、后NDVI时序曲线，3种方法在建筑扰动识别中准确性均较高，相比之下，S-G法显示曲线细节特征的优势更加明显。

综上，采矿扰动和建筑扰动像元重建中，对比3个参量(时间差、振幅差、恢复速率差)，S-G法均优于AG法和D-L法。

3.2.2 Jacknife法评价拟合结果

Jacknife法可以更加直观地表示各个拟合方法的重建效果，假想NDVI时间序列曲线较平滑，变化幅度较小，3种方法重建后曲线几乎与原曲线重合，因此，对较平滑的NDVI时序曲线，3种方法均较好地保持了原有曲线的特征。重建前、后像元重建效果见图5(a)—(f)。由图5(b)可见，当曲线中出现2个不连续突降噪声，即第3和10期NDVI值突降至原值的50%，重建后曲线与假想NDVI时序曲线几乎重叠，因此，3种方法均能识别不连续突降点并将其去除。由图5(c)可见，当曲线中出现2个连续突降噪声，即第10和11期NDVI值突降至原值的50%，AG法和D-L法重建后曲线较平滑，更加接近假想曲线，而S-G法重建后第11期NDVI值出现小幅降低，这表明AG法和D-L法无法识别连续2期突降噪声，而S-G法则可以。由图5(d)可见，当曲线中出现3个连续突降噪声，即第10、11和12期NDVI值突降至原值的50%，3种方法重建后曲线在第10、11和12期出现不同程度的降低，即3种方法都能够识别出3个连续噪声，但未加入噪声的曲线部分，AG法和D-L法重建后与原曲线偏离较大，而S-G法几乎与原曲线重合。由图5(e)可见，当模拟建筑扰动时序数据时，3种方法都能识别出建筑扰动，但识别出扰动的时间有不同程度的推迟，AG法和D-L法出现偏离原始曲线的现象。由图5(f)可见，当模拟采矿扰动时序数据时，3种方法都能识别出采矿扰动，但S-G法重建后曲线更加接近原始曲线变化特征。

图5 重建前、后像元重建效果

3.2.3 回归估计标准差评价拟合结果

Jacknife法是对重建结果的定性评价，本文还引入RMSE以定量评价重建结果，该值越小，表示重建效果越好，反之，则越差。6类像元重建后NDVI值的RMSE见表4。由表4可见，AG法和D-L法相差较小，S-G法均小于前2种方法。由图5可见，S-G法能够更好地反映时序曲线的细节变化，在出现突变噪声的情况下，可以保持原始曲线的基本特征，因此，如果NDVI时序数据用于识别建筑、采矿等扰动，应选用S-G法对其进行重建。

表4 6类像元重建后NDVI值的RMSE

4 结论

使用TIMESAT软件中的AG法、D-L法、S-G法这3种方法拟合重建兖州市2014年23期NDVI时间序列数据，对比分析3种方法提取的SOS、EOS、LOS3个农作物物候参数的准确度。使用3种方法重建韦兹县采矿扰动点和建筑扰动点，构建时间差、振幅差和恢复速率差3个参量评价3种方法重建后曲线对原始NDVI时序曲线扰动特征的保持度；使用Jacknife法设计不同特征的NDVI时序变化曲线，更直观地评价3种方法拟合效果。主要结论如下：

(1)在物候参数提取应用场景中，3种方法总体差别较小，3种方法拟合重建后均可提取精度较高的物候参数，提取的SOS、EOS、LOS等物候参数接近于站点数据；相比S-G法，AG法和D-L法保持NDVI时序曲线整体变化特征的能力更强，提取的冬小麦和夏玉米的SOS和EOS更接近于站点数据。3种方法都能比较准确地识别和提取植被生长周期，在识别耕地和林地生长周期数方面差异较小。

(2)人类活动扰动识别应用场景中，3种方法都体现出对原始数据曲线很好的保真性和保持度。特别是采矿扰动识别中，S-G法的时间差、振幅差和恢复速率差3个参量的评价结果好于AG法和D-L法，因此S-G法在滤波时能够最大限度地保留时序曲线细节变化，对扰动发生的时间和强度等信息比较敏感，识别精度优于AG法和D-L法。