3种时序Ndvi重构方法对比分析

丁灿彧 ，胡希军 ，吴德政 ，3，刘伟乐

（1.中南林业科技大学 风景园林学院，湖南 长沙 410004； 2.厦门城市职业学院，福建 厦门361000；3.集美大学 美术学院，福建 厦门 361000；4.中南大学，湖南 长沙 410004）

植被指数是最常用的监测植被变化状况的遥感数据，如Ndvi（归一化植被指数）时间序列数据，它来源于遥感卫星的传感器，可以准确地反映陆地表面植被年际变化特点、生长状态等情况，Ndvi时间序列数据携带了关于土地表面性质的有价值信息，已经被证明适用于长期监测土地利用情况和覆盖变化情况[1-3]。由于植被生长具有一定的规律性，植被指数也就具备规律性变化的特征，因而时序Ndvi反映植被生长进程的曲线一般是连续平滑的。然而，在这些时间序列数据中几乎总是有来自云覆盖、大气变化、和双向反射等的干扰，以及冰雪、地表水等随机因素的干扰，而且这些干扰显示为不可预测的噪音，使实际的时序Ndvi曲线不能清楚地反映植被时令规律性变动趋势，严重影响了对土地覆盖和陆地生态系统的监测与模拟[4-7]。

为此一些研究者常用基于最大值合成法（MVC）来消除Ndvi数据集中的这些噪声，尽管如此，但Ndvi数据集还是存在一些噪声残差，影响数据集的使用[8-10]。随着相关研究的进一步深入，Ndvi数据集中的这些噪声对研究工作产生了很大的影响，使得研究的普遍性和可靠性指数都无法得到保证。因此，构建合理的Ndvi时间序列数据拟合重建方法，来降低Ndvi时序数据中的噪声水平。当前常见的Ndvi时间序列重建方法包括：阈值法、滤波法、曲线拟合法。阈值去噪法包括Viovy等提出的最佳指数斜率提取法(BISE)和Lovell 等在此基础上提出的改进算法；基于滤波的拟合法主要包括Savitzky-Golay滤波、傅立叶变换等；曲线拟合法以双Logistic 函数拟合法和非对称性高斯函数拟合法为代表[11-14]。

研究者在研究时要根据不同的研究区域、植被类型，选取不同的重建方法。由于各种条件的不同，或者相关条件的变化，很难用得到一种适合各种环境的重建方法，因此对最优的重建方法难以定论。由于不同区域的大气环境及地面条件不同，植被覆盖类型及农作物制度也具有一定的差别，结合国内Ndvi时序重建定量研究的现状，本研究以厦门岛作为研究区，以林地为实验类型，利用TIMESAT软件比较3种主要的拟合重建算法（S-G滤波法、双Logistic曲线拟合和非对称高斯函数拟合）对 2013年MODIS 16 d合成的Ndvi时间序列数据进行处理，对比分析拟合之后的时序Ndvi平滑曲线，以及通过Ndvi序列上包络曲线拟合效果和保真性的对比，得到最适合长厦门岛的时序Ndvi重建方法，为基于Ndvi时间序列研究的降噪处理方法选择提供参考。

1　研究区概况及数据

1.1　研究区概况

本研究以厦门岛为研究区。厦门岛面积约128.14 km2，是福建省第4大岛屿。岛内地势南高北低，南部多为丘陵，北部多为阶地、平原，全岛阶梯状地形明显。厦门岛是典型的亚热带海洋性气候，冬无严寒，夏无酷暑，全年温差小，年均气温为20 ℃，年均降水量在1 200 mm左右。

1.2　数据及处理

研究使用的Ndvi时序数据来自MODDIS传感器的产品之一MOD13Q1v005，该传感器是以EOS/Terra卫星为载体进行工作的。获得的Ndvi数据经过大气校正和几何校正的数据预处理过程。研究数据是从美国NASA LP DAAC工作组网站下载，共有23景影像，时间跨度为2013年1月到12月全年。

Ndvi标准值的范围是-1.0～1.0， 而MODIS13Q1NDVI产品得到的是DN值，其范围是-3 000～10 000，其中由于各种外界条件原因导致某些区域未获得相应的遥感数据，-3 000即为这些数据的填充值。可以利用遥感影像数据处理软件ENVI5.3，将ND值转换成Ndvi数据值，转换公式如下：

MODIS植被指数产品（MODIS13 Q1）为基于16 d最大值合成法的植被指数数据，空间分辨率为250 m。该产品是成品，经过辐射定标、大气校正和几何校正等数据预处理过程，产品中时序Ndvi数据用于本实验研究，之后投影转换及裁剪得到研究区2013年Ndvi时间序列数据。

2　研究方法

（1）Savitzky-Golay滤波法（S-G）

首先根据云状态（像元可信度）对Ndvi数列进行线性插值，然后利用S-G滤波器得到插值后曲线的模拟生长趋势线，再根据上包络线得到新的Ndvi曲线，将上述过程进行数次迭代并设置拟合影响系数作为迭代退出条件，最终得到较为平滑又能反映Ndvi数值变化趋势的时间序列曲线[15-17]。

式（2）中：Y是指Ndvi原始值，Y*是Ndvi拟合值，j是原始Ndvi数组的系数，Ci是第i个Ndvi值的卷积系数，N是滑动数组宽度（2m+1），m是待拟合点左右两端各需的点数。

（2）非对称高斯函数拟合法（A-G）

非对称高斯函数拟合法(A-G)，是一种最小二乘拟合算法，该算法是非线性的，并且基于不同的高斯函数[18-20]。拟合函数的基本形式为：

其中为高斯函数：

式（3）、式（4）中，t表示t时刻的Ndvi值；c1和c2为控制曲线的基准和幅度；a1决定峰值和谷值的位置；a2、a3、a4、a5为控制曲线左、右半部分的宽度和陡峭度。

（3）双Logistic函数拟合法（D-L）

双Logistic函数拟合方法（Double Logistic Function- fitting，DLF）是Beck等在2006年提出的一个新算法，它也是一种半局部拟合算法。首先将整个时间序列Ndvi数据按照极大值和极小值，分成若干个区间，然后在区间内利用双逻辑函数进行局部拟合。函数的基本表达式为：

双Logistic曲线拟合法是一种半局部拟合方法，其局部拟合方式与非对称高斯拟合方法类似（式（3））。采用双Logistic函数（式（5）），基于整体拟合函数（式（4））将各局部拟合函数的特征加以综合，重建新的Ndvi时间序列曲线[21-22]。

在基于植被Ndvi的各种应用研究中，植被全年的生长状况一般是用植被生长季峰期的Ndvi值来表现的。厦门岛七八月份雨水较少，天气状况相对良好，此时的遥感影像比较全面，因此可以选取7、8月份的Ndvi数据，分析3种方法拟合之后的上包络线拟合效果。此外，利用时序Ndvi数据提取植被物候信息，或者研究全年及生长季峰期的Ndvi平均状况，都需要选择一种方法，该方法既能去除原始Ndvi数据的噪点，又能较好的保持一些品质好的Ndvi值的真实性。为了比较D-L函数拟合法、S-G滤波法和A-G函数拟合法3种算法对于厦门岛MODIS 16 d合成的时序Ndvi数据拟合重建的效果，将从以下两个方面进行分析：

（1）比较3种算法拟合之后的Ndvi均值与原始Ndvi数据均值的差别，以及计算拟合前后时序Ndvi曲线的相关系数；

（2）分析全年以及生长季峰期重建后的Ndvi数据与高质量Ndvi原始值的偏离程度。

其中，相关系数反映了拟合后的时序Ndvi曲线保持原始植被生长特征的能力，如下式（6）所示：

式（6）中：NDVIoi表示年内时间序列中第i期拟合处理前的Ndvi值，NDVIpi表示相对应的拟合处理之后的Ndvi值；分别表示拟合处理前后的年内时序Ndvi均值[23]。

回归估计标准差（如式7），表示重建前后Ndvi时间序列的拟合效果，该值表明重建后的时序Ndvi数据与原始值之间的平均差别程度，反映重建前后Ndvi值的代表性强弱。回归估计标准差值越小，说明拟合之后的Ndvi值的代表性越强。

式（7）中，NDVIoi表示年内时间序列中第i期拟合处理前的Ndvi值，NDVIpi表示相对应的拟合处理之后的Ndvi值；N表示高质量（DN=0.1）的Ndvi像元样本数量[24]。

3　结果与分析

3.1　3种方法整体拟合效果

3种方法拟合的效果与原始序列的曲线基本一致，对全年23期的数据均有不同幅度的提升，3种方法提升的幅度大概都为4%，有较好的重建效果。同时3种方法也都检测出了Ndvi时间序列中的突增点和突降点，并做出了较好的纠正，例如16、21期的突降点和22期的突增点。对序列中的异常值（如第3期填充值）也都做出了一定的修正。

数据的拟合情况与实验区域本身情况有关，第3期出现异常值，可能原因是由于年初冰雪覆盖。选取的具体像元区域为厦门岛林地，厦门岛林地植被均为四季常青植物，而从第2期到第8期Ndvi值总体来说较小，这是因为这段时期厦门岛正处于梅雨季节，雨水及大气状况、云覆盖等降低了Ndvi的值，还有一部分原因是受气候影响，该时期植物生长缓慢。9～16期Ndvi时间序列数值较大，这是因为此时期内厦门岛气候炎热，植被生长旺盛，恰好反映了植被的生长规律。第16期和22期的突降值也可能是由于天气状况导致的。

3.2　相关系数分析比较

根据式（6），分别计算全年及生长季内（第7到16期）Ndvi时间序列重建之后的保真度，即基于上述3种方法拟合之后的相关系数（如图1）。

可以看出，A-G函数拟合重建的时间序列与植被原始Ndvi时间曲线的相关性最好，全年23期的Ndvi时序数据拟合的相关系数为0.764 0，生长季内Ndvi拟合的相关系数为0.782 1。S-G滤波法比A-G拟合略低，全年和生长季内的相关系数分别为0.731 9、0.752 0；而D-L函数拟合法由于受到噪声的影响，相关系数分别为0.520 0、0.771 3，全年拟合效果没有前两种方法的好，对整个时间序列曲线整体特征的保持性相对较差，但是D-L函数拟合法在生长季期内，拟合效果较前两种方法差别不大，由图2也可以看出相同的结果。

图1 原始的和3种方法重建之后时间序列比较Fig. 1　Comparison of time series after reconstruction of the original and three methods

图2　基于3种方法的全年及生长季内Ndvi序列重建前后相关系数比较Fig. 2　Comparison of correlation coef ficients before and after reconstruction of Ndvi sequences based on three methods in the whole year and growing season

与全年的拟合效果相比，生长季内的时序Ndvi采取任何一种方法重建前后的相关性都更好，因为研究区域在生长季内的时间内天气状况良好，由云覆盖、大气影响等原因导致的噪声相对较少，重建的可靠性较高。而在生长季以外的时间段，噪声较大，有明显的突增点和突降点，重建的可靠性相对较低。

3.3　回归估计标准差比较

可以用生长季的最大值来反映一个地区植被生长最好的情况，厦门岛时间为5—6月（第10～12期）。而由于受到噪声的污染会使得Ndvi值偏低的这种情况会给相关研究结论带来偏差。因此可用式（7）计算3种方法重建之后的回归估计标准差，来描述Ndvi重建前后之间的平均差异程。图3反映了基于A-G拟合法、D-L函数拟合法、S-G滤波法3种算法重建后的Ndvi时间序列在生长季峰期（第10～12期）相对于原始Ndvi时间序列上包络线的回归标准差。图3反映了基于3中算法拟合重建后的Ndvi时间序列在全年相对于原始Ndvi上包络线的回归标准差。

从图3中不难发现，在生长季峰期时间内，S-G滤波法拟合效果要优于A-G函数和D-L函数拟合算法，后两个方法的拟合效果相近，因此A-G滤波法更能反映植被生长的最佳状况。在全年内3种方法的拟合效果相近，A-G和S-G滤波法拟合效果要稍好于D-L函数拟合法。

4　结论与讨论

图3　基于3种算法重建结果相对原始Ndvi序列上包络线在季峰期与全年的回归标准差Fig. 3　Based on 3 kinds of algorithm reconstruction results relative envelope of the original Ndvi sequence in season and annual Fengqi regression standard deviation

本研究使用厦门岛区域2013年全年23期Ndvi时间序列数据集，3种常用的Ndvi时间序列拟合重建算法（A-G函数拟合法、D-L函数拟合法、S-G滤波法），拟合得到了原始的和3种方法重建之后Ndvi时间序列比较图。由于噪声通常会降低Ndvi理论值，因此本研究选择拟合后的时间序列相对于原始曲线的上包络线拟合效果及保持高质量Ndvi点真实值的性能方面比较分析A-G函数拟合、D-L函数拟合，S-G滤波法拟合3种算法的优劣。

结果表明：

（1）A-G函数拟合法和D-L函数拟合法重建的Ndvi时间序列更能反映植被生长的曲线特征。

整体拟合效果方面，3种Ndvi重建拟合方法均能不同程度的降低噪点，同时拟合之后的所有采样点Ndvi平均值均高于Ndvi原始数值，这也就证实了大部分噪声使得Ndvi理论值降低这一观点。D-L函数拟合法和A-G函数拟合法重建得到的拟合结果相类似，S-G滤波法在重建过程中有“过度拟合”的特点，对噪声比较的敏感，该方法处理的Ndvi均值比前两种方法较小。

（2）S-G滤波法和A-G函数拟合法在保持高质量Ndvi真实值特征方面效果较好。

A-G函数拟合重建时序Ndvi曲线与植被原始Ndvi时间曲线的相关性最好，S-G滤波法比A-G拟合略低。D-L函数拟合法由于受到噪声的影响，在全年内拟合效果没有前两种方法的好，对时间序列曲线整体特征的保持性相对较差，但是D-L函数拟合法在生长季期内，拟合效果和前两种方法差别不大。

（3）A-G滤波法拟合重建的Ndvi更能反映植被生长的最佳状况。

在生长季峰期时间内，S-G滤波法拟合效果要优于A-G函数和D-L函数拟合算法，后两个方法的拟合效果相近。

通过比较3种拟合算法的优点，A-G函数拟合和S-G滤波两种方法在一些方面有比较大的优势，S-G滤波法虽然在拟合季峰期Ndvi时间序列数据时效果很好，更能反映植被的最佳生长曲线，但是它对噪声比较敏感，而且也容易受到人为因素的影响，具有一定缺陷。总体来说，A-G函数拟合法在重建厦门岛林地的Ndvi时间序列方面有一定的优势。

参考文献：

[1]李杭燕,颉耀文,马明国.时序NDVI数据集重建方法评价与实例研究[J].遥感技术与应用,2009,24(5):596-602.

[2]杨 磊,张 梅,罗明良,等.基于MODIS NDVI 的川中丘陵区植被覆盖度景观格局变化[J].生态学杂志,2013,32(1):171-177.

[3]田庆久,阂祥军.植被指数研究进展[J].地球科学进展, 1998,13(4):327-333.

[4]Ackerman SA, Strabala KI, Menzel WP,et al.Discriminating Clear Sky from Clouds with MODIS[J]. Journal of Geophysical Research, 1998, 103(32) : 141-157.

[5]Yu F, Price K P, Ellis J,et al. Satellite Observations of the Seasonal Vegetation Growth in Central Asia: 1982-1990 [J].Photogrammetric Engineering & Remote Sensing, 2004, 70(4):461-469.

[6]Deering D W. Rangeland Re flectance Characteristics measured by aircraft and Spacecraft Sensors[D]. Texas,USA: Texas A&M University, 1978:338.

[7]Cihlar J, Ly H, Li Z Q,et al.Multitemporal, Multichannel AVHRR Data Sets for Land Biosphere Studies-artifacts and Corrections[J].Remote Sensing of Environment, 1997(60):35- 57.

[8]刘伟乐,林 辉,孙 华. 基于GF-1遥感影像湿地变化信息检测算法分析[J]. 中南林业科技大学学报,2015,35(11):16-20.

[9]吴文斌,杨 鹏,唐华俊,等.两种NDVI 时间序列数据拟合方法比较[J].农业工程学报, 2009,25( 11) : 183-188.

[10]ViovyN,ArinoO,Belward AS. The Best Index Slope Extraction(RISE):A Method for Reducing Noise in NDVITime-series[J].International Journal of RemoteSensing,1992,13(8):1585-1590.

[11]Sellers P J,Tucker(J,Collatz U J,et al. A Revised Land Surface Parameterization (SiB2) for Atmospheric VCMs. PartII:The Veneration of Vlobal Fields of Terrestrial Biophysical Parameters from Satellite Data[J]. Journal of Climate, 1996,9(4):706-737.

[12]Jnsson P, Eklundh I. Seasonality Extraction by Function Fiting to Time-series of Satellite Sensor Data[J]. IEEE Trans actions on Ueoscience and Remote Sensing, 2002,40(8)1824-1832.

[13]Pettorelli N, Vik J O, Mysterud A,et al.Using the satellite-derived NDVI to assess ecological responses to environmentalchange[J]. Trends in Ecology & Evolution, 2005,20(9): 503-510.

[14]梅安新.遥感导论[C].北京:高等教育出版社,2008.

[15]Holben B N. Characteristic of maximum value composite images for temporal AVHRR data[J]. International Journal of Remote Sensing, 1936, 7(11): 1417-1434.

[16]郭 妮.植被指数及其研究进展[J].干旱气象,2003,21(4):71-75.

[17]Kaufmann R K,et al.Effect of orbital drift and sensor changes on the time series of AVHRRvegetation index data[J]. IEEE Transactions on Geoscience and Remote Sensing, 2000, 38(6):2231-2597.

[18]陈朝晖,朱 江,徐兴奎.利用归一化植被指数研究植被分类、面积估算和不确定性分析的进展[J].气候与环境研究,2004,9(4): 687-696.

[19]陈述彭,赵英时.遥感地学分析[C].北京:测绘出版社,1990.

[20]覃文汉,项月琴.植被结构及太阳/观测角度对NDVI的影响[J].环境遥感,1996,11(4): 285-290.

[21]Burgess D Lewis P, Muller J PAL. Topographic effects in AVHRR NDVI data[J].Remote Sensing of Environment, 1995,54(3):223-232.

[22]杨俊泉,莫伟华.马尾松毛虫危害区植被指数时序变化特征研究[J].国土资源遥感,1997, 34(4):7-13.

[23]Huete AR, Liu H. An error and sensitivity analysis of the atmospheric and soil-correctingvariants of the NDVI for the MODIS-EOS[J]. IEEE Trans. Geosci. Rem. Sens, 1994,32(4):897-905.

[24]陈朝晖,朱 江,徐兴奎.利用归一化植被指数研究植被分类、面积估算和不确定性分析的进展[J].气候与环境研究,2004,9(4): 687-696.