Comparison of Eight Techniques for Reconstructing Multi-Satellite Sensor Time-Series NDVI Data Sets in the Heihe River Basin, China

Abstract

More than 20 techniques have been developed to de-noise time-series vegetation index data from different satellite sensors to reconstruct long time-series data sets. Although many studies have compared Normalized Difference Vegetation Index (NDVI) noise-reduction techniques, few studies have compared these techniques systematically and comprehensively. This study tested eight techniques for smoothing different vegetation types using different types of multi-temporal NDVI data (Advanced Very High Resolution Radiometer (AVHRR) (Global Inventory Modeling and Map Studies (GIMMS) and Pathfinder AVHRR Land (PAL), Satellite Pour l’ Observation de la Terre (SPOT) VEGETATION (VGT), and Moderate Resolution Imaging Spectroradiometer (MODIS) (Terra)) with the ultimate purpose of determining the best reconstruction technique for each type of vegetation captured with four satellite sensors. These techniques include the modified best index slope extraction (M-BISE) technique, the Savitzky-Golay (S-G) technique, the mean value iteration filter (MVI) technique, the asymmetric Gaussian (A-G) technique, the double logistic (D-L) technique, the changing-weight filter (CW) technique, the interpolation for data reconstruction (IDR) technique, and the Whittaker smoother (WS) technique. These techniques were evaluated by calculating the root mean square error (RMSE), the Akaike Information Criterion (AIC), and the Bayesian Information Criterion (BIC). The results indicate that the S-G, CW, and WS techniques perform better than the other tested techniques, while the IDR, M-BISE, and MVI techniques performed worse than the other techniques. The best de-noise technique varies with different vegetation types and NDVI data sources. The S-G performs best in most situations. In addition, the CW and WS are effective techniques that were exceeded only by the S-G technique. The assessment results are consistent in terms of the three evaluation indexes for GIMMS, PAL, and SPOT data in the study area, but not for the MODIS data. The study will be very helpful for choosing reconstruction techniques for long time-series data sets.

1. Introduction

Analysis of Normal Difference Vegetation Index (NDVI) time-series data is becoming increasingly important for ecological research on environmental dynamics and climate change [1–3], vegetation dynamics [4–8], land cover change [9–11], and animal species distribution [12]. Analyzing NDVI time-series data has been a useful tool for studying climate, vegetation, and animal distribution, and performance at large spatial and temporal scales [13]. Although NDVI data sets are pre-processed to reduce noise from sensor resolution and calibration, digital quantization errors, ground and atmospheric conditions, and orbital and sensor degradation, some noise is still present in the downloadable data sets, including noise that results from cloud cover, poor atmospheric conditions, and unfavorable sun-sensor-surface viewing geometries [13–19]. Though the standard maximum value compositing technique (MVC) [14] was used for the main data sets, (e.g., the vegetation index products from Global Inventory Modeling and Map Studies (GIMMS), Advanced Very High Resolution Radiometer (AVHRR) data, Pathfinder AVHRR Land (PAL) data, Satellite Pour l’ Observation de la Terre Vegetation (SPOT) VEGETATION (VGT) data, and TERRA or AQUA Moderate Resolution Imaging Spectroradiometer (MODIS) data), these data sets still include some residual noise. Therefore, NDVI time-series data need to be smoothed before being used.

Several methods have already been developed to reconstruct NDVI time series. They can be grouped into five types: (1) threshold methods, including the best index slope extraction technique (BISE) [20] and the modified BISE technique (M-BISE) [21]; (2) filter-based methods, including running medians (4253H) [22], the ARMD3-ARMA5 filter technique [23], the Savitzky-Golay filter technique (S-G) [24], the mean value iteration filter (MVI) technique [25], the changing-weight filter (CW) technique [26], and phenology-preserving filtering (PP) [27]; (3) function fitting methods, such as the fast Fourier transform (FFT) [28], the temporal window operation (TWO) [29], the harmonic analysis of time series (HANTS) [30], double logistic function fitting (D-L) [31–34], and asymmetric Gaussian function fitting (A-G) [35]; (4) comprehensive methods [36–38]; and (5) other methods, such as wavelet transform (WT) [39], data assimilation (D-A) [40], the Whittaker smoother (WS) technique [41–43], and iterative interpolation for data reconstruction (IDR) [44].

Comparisons of these techniques have shown that each has its own advantages and drawbacks [15,20,24,25,39,44,45]. Viovy et al. [20] compared the BISE and MVC techniques and concluded that BISE is superior to MVC in terms of de-noising. Jönsson and Eklundh [35] demonstrated that A-G is superior to BISE and the Fourier-based technique. Chen et al. [24] demonstrated that S-G, BISE, and the Fourier-based transformation (FT) are effective techniques for constructing high-quality NDVI time-series data sets, and they showed that S-G is the best of those three techniques. However, Jönsson and Eklundh [46] explained that the A-G technique could outperform an S-G filter and an alternative harmonic analysis. Later, Ma and Veroustraete [25] reported that the MVI performed better than the M-BISE and FT. Beck et al. [33] showed that the new version technique of D-L is better than both the FT and A-G because it can handle outliers effectively and estimate phenological event parameters. Lu et al. [39] developed the WT technique and compared it with the BISE, FT, and S-G techniques. The results indicated that the WT is the best for removing noise. Later, Hird and McDermid [15] demonstrated that the general superiority of the D-L and A-G function-fitting techniques over four alternative filtering techniques (S-G, 4523H, MVI, and ARMD3-ARMA5) by performing an empirical comparison. Julien and Sobrino [44] presented the IDR technique and showed that it can provide the best profile reconstruction for most land cover classes (compared with HANTS and D-L). More recently, Atkinson et al. [45] found that it was necessary to tune the A-G, D-L, FT, and WS techniques according to the number of annual growing seasons to produce reliable fits. Zhu et al. [26] presented the CW technique and showed that it was more effective than A-G, D-L, and S-G in preserving the curve shape, as well as the timing and amplitude of the local maxima/minima in the time series, for a broad range of phenologies. Jiang et al. [27] developed the PP technique based on CW, and showed that it performs much better than the CW filter for different levels of noise.

Although many studies have focused on NDVI noise reduction techniques comparison, there are still some aspects that need to be improved and completed. First, previous work has been restricted by a small selection of the available noise reduction techniques, and the majority of the literature compares two or three techniques, typically, single novel technique is compared to one or two widely-known techniques [15]. Only a few studies compare more than three techniques, such as the studies by Atkinson et al. [45], Hird and McDermid [15], and Jiang et al. [27]. Hird and McDermid [15] compared six NDVI de-noising techniques, two of which (the ARMD3-ARMA5 filter and the 4253H) were rarely used in the literature. Moreover, approximately ten new techniques have been developed in the last decade, and some of these new techniques have not been rigorously compared against other techniques, IDR has not been compared to A-G, and WS has not been compared to S-G, and CW. Second, previous comparisons have been focused on the pixels or regional scale, and few studies have paid attention to vegetation types, which could be important for selection optimum noise-reduction techniques [33,45]. Third, almost all noise-reduction technique comparisons have been based on one sensor. Few published studies have tested the differences between sensors. Finally, previous evaluations have principally relied on qualitative rather than quantitative assessments, and they generally have not addressed factors that could affect the performance of noise-reduction techniques [15]. Currently, few studies refer to the effects of different evaluation methods on assessment of outcome.

The main purpose of the present study was to support other basic studies by comparing eight techniques for representing time-series NDVI data to support other basal studies. These techniques include the A-G, M-BISE, CW, D-L, IDR, MVI, S-G, and WS techniques. Several studies have indicated that the A-G and D-L are superior to other fitting techniques and filters [15,33,35,46]. The S-G has been widely used in the literature [47–50] due to the advantages cited by Chen et al. [24]. Additionally the A-G, D-L, and S-G have been made available as software that can be easily downloaded freely from the website [51]. The M-BISE and the MVI techniques are simple and effective de-noising techniques that can be implemented readily with a few lines of code. The CW, IDR, and WS are innovative techniques that have been used recently for time-series reconstruction using remote sensing data. The potential advantages of these techniques could lead to widely application in the future.

Four widely used long time-series NDVI data sets were used in this study (

Table 1. Long time-series Normalized Difference Vegetation Index (NDVI) data sets used in the study.

Products	Time Period	Temporal/Spatial Resolution	Sensors	Data Source
GIMMS	January 1982–December 2006	15 d/8 km	AVHRR	Cold and Arid Regions Science Data Center [52]
Pathfinder	July 1981–July 1994 January 1995–December 2000	10 d/8 km	AVHRR	Cold and Arid Regions Science Data Center [53]
SPOT VEGETATION (S10)	April 1998–April 2013	10 d/1 km	SPOT VGT	VITO Earth observation [54]
MOD13A2	February 2000–February 2013	16 d/1 km	Terra/MODIS	The Land Processes Distributed Active Archive Center [55]

[52–55]). All of these data sets were reconstructed using the eight techniques mentioned above combined with information about the vegetation types in the study area. Three statistical indexes were used to evaluate the reconstruction results. The aims of this research were to identify the optimal de-noising technique for each vegetation type and data set in the study area; then to analyze the differences between the optimal techniques for each vegetation type and the data set; and finally to explore the effects of different evaluation methods on the results.

2. Method

2.1. Study Area

The Heihe River Basin, in the middle of the Hexi Corridor of Gansu Province, is the second largest inland river basin in China. It is located between 97°1′–102°0′E and 37°7′–42°7′N, with an area of approximately 143,000 km². The elevation in the Heihe River Basin ranges from 5500 m in the south to 1000 m in the north [56]. It consists of three regions: the upper mountainous area (the source of the Heihe River), the middle oasis area (including towns, such as Zhangye and Jiuquan), and the lower terminal arid area around Ejina [57], from an upstream to downstream, the environment changed from glacier, to frozen soil, alpine meadow, forest, irrigated farmland, a riparian ecosystem, and a desert (the Gobi) [58]. Figure 1 shows the location and vegetation types of the study area.

2.2. Data and Processing

Four long time-series NDVI data sets of remote sensing were used in this study, and detailed information about the data sets is shown in Table 1. All of these NDVI data sets were compiled for multi-day data using the MVC technique. To ensure data quality, a series of processes have been performed for these data sets. GIMMS data sets had been corrected for calibration, view geometry, and volcanic aerosols and were verified using a stable desert control point [59,60]; PAL data sets pre-processing steps include navigation, inter-satellite calibration and partial correction for Rayleigh scattering [61]; The MODIS data (MOD13A2) have been performed atmospherically corrected bidirectional surface reflectance that had been masked for water, clouds, heavy aerosols, and cloud shadows, and the accuracy had been assessed over a widely distributed set of locations and time periods via several ground-truth and validation efforts [62]. The processes of SPOT VGT NDVI data sets include atmospheric correction, radiometric correction, and geometric correction [63]. A subset was extracted for SPOT VGT NDVI using VGTExtract software (version 2.1.0), a free vegetation extraction tool produced by VITO. The status maps were used to identify useful NDVI. For MOD13A2, the MODIS Reprojection Tool was used to extract the desired bands and to re-project them as WGS 1984. Quality assurance products from Terra MODIS were used to eliminate the obvious error noised data. Only pixel reliability values of 0 and 1 were accepted as reliable pixel. The continental GIMMS and pathfinder NDVI data sets were provided by the Cold and Arid Regions Science Data Center [52,53]. Pixels with obvious error NDVI data, such as values less than 0 or null values, were replaces by the mean value of the same time of multiyear for all of the four NDVI data sets. After the obvious errors were removed, the data were fitted with each of the techniques.

The 2001 Vegetation Map of China [64], with a spatial resolution of 1 km, was used to identify the main vegetation types in the study area. To thoroughly analyze the reconstruction results for the main types of vegetation, pixels with 80% to 100% cover of the same vegetation type were extracted as representative homogeneous pixels in the study area. All of the homogeneous pixels for grassland, meadow, crops, desert, and shrub areas were selected for further analysis.

2.3. Candidate Time Series Reconstruction Techniques and Parameters

Eight de-noising techniques were selected for this research, and the main objective of each technique is shown in

Table 2. Summary of the NDVI time-series reconstruction techniques selected for comparison.

Candidate Techniques (ab.)	Description	References
Modified-best index slope extraction (M-BISE)	Compares the current term value with the previous and the next term within a predefined sliding window, and replaces these values with the mean value of the previous and the next values if the percentage difference is greater than a predefined threshold.	[20,21]
Asymmetric Gaussian function-fitting (A-G)	Fits local, nonlinear functions at intervals around the local maxima and minima, then merges these into a global function describing the full NDVI time series.	[35]
Double logistic function fitting (D-L)	Uses six parameters to model the NDVI time series with a double logistic function. These parameters are the winter NDVI (wNDVI), maximum NDVI (mNDVI), two inflection points, one as the curve rises (S) and one as it drops (A), and the rate of increase or decrease (mS and mA) of the curve at the inflection points.	[33]
Savitzky-Golay filtering (S-G)	Based on a simplified least-squares-fit convolution for smoothing and computing derivatives of a set of consecutive values (a spectrum). The convolution can be understood as a weighted moving average filter with weighting given as a polynomial of a certain degree. The weight coefficients, when applied to a signal, perform a polynomial least-squares fit within the filter window. This polynomial is designed to preserve higher moments within the data and to reduce the bias introduced by the filter.	[24]
Mean value iteration filtering (MVI)	Iteratively compares each date with the average of the dates before and after it, replacing the date with this average if the difference is above a certain threshold. The maximum difference date value will be removed in an iteration process. Iteration will stop when all differences are less than the threshold.	[25]
Whittaker smoother (WS)	Based on penalized least squares, fits a discrete series to discrete data and penalizes the roughness of the smooth curve. In this way, it balances the reliability of the data and roughness of the fitted data.	[41,43,45]
Iterative interpolation for data reconstruction (IDR)	Creates an alternative NDVI time series by computing the mean of the immediately preceding and following observations and comparing it to that of the original time series, replacing the original data with the alternative time-series data if the maximum difference between the alternative and original time series is above a certain threshold.	[44]
Changing-weight filtering (CW)	Filters an NDVI time series with a three-point changing-weight filter and replaces the local maximum/minimum points in a growth cycle.	[26]

. All of these techniques have been used or have the potential to be used in real applications. Detailed descriptions of the techniques can be found in the literature (Table 2). All of the techniques run by compiling Interactive Data Language code (IDL), except for A-G, D-L, and WS, the A-G and D-L are run in TIMESAT [46] and the WS is run in Matlab.

The parameters used for the techniques greatly affect the reconstruction results. As described by Atkinson et al. [45], for any of the tested techniques, finding a single set of technique parameters that is appropriate for the study area can be challenging for a landscape that is diverse and complex. To obtain the optimum parameters for each of the eight techniques, different fitting criteria were set to reconstruct the NDVI time-series for the four sensors over the five main vegetation types in the study area.

By comparison, we used 0.1 as the threshold for the M-BISE technique and a 20% multi-year average for the MVI filter threshold. For S-G, m (the half-width of the smoothing window) was 4, d (the degree of the polynomial) was 6, and the IDR threshold was 0.02. These values are similar to their authors’ recommended. For WS, λ was 15 because the entire research area contains single-season vegetation types [45]. There were three parameters for the CW technique: the window width was 7 for MODIS and GIMMS data and 9 for the SPOT and PAL data. The threshold for removing a false local maximum and minimum point and a noise point were the same for the four data sets, 0.1 and 0.05, respectively. The D-L and A-G techniques were adjusted interactively in the TIMESAT software to arrive at close fitting results. In TIMESAT, the median filter option with a parameter value of 2 was chosen to remove spikes and noise.

2.4. Technique Evaluation

Three statistical indicators were used to evaluate the performance of each noise reduction technique: (1) the root mean square error (RMSE); (2) Akaike’s Information Criterion (AIC); and (3) Bayesian Information Criterion (BIC).

RMSE is a well-accepted absolute goodness-of-fit indicator for continuous response variables that describes the difference between the observed and predicted values in the appropriate units [65]. Here, the RMSE indicates the difference between the mean NDVI time series obtained from the eight noise- reduction techniques (assumed to be accurate) and the corresponding experimental time series to which noise reduction has been applied [15]. It is calculated using Equation (1).

$$\text{RMSE}=\sqrt{\frac{{\Sigma }_{t=1}^N{\left(VI*(t)-VI(t)\right)}^2}{N}}$$

(1)

where VI* is the resultant NDVI value, VI is the mean NDVI value obtained from the eight noise-reduction techniques (mean NDVI), and N is the number of time points.

The Akaike’s Information Criterion (AIC) [66] and the Bayesian Information Criterion (BIC) [67] are criteria used for selecting technique from a finite set of techniques. Both criteria are closely related and can measure the efficiency of the parameterized technique in terms of predicting the data, but the BIC penalizes the free parameters more severely than does the AIC. The AIC and BIC evaluation index have used to evaluate the performance of de-noised techniques by Atkinson et al. [45]. They are calculated using Equations (2) and (3), respectively.

$$\text{AIC}=2k+n\left[\text{ln}(\mathit{RSS})\right]$$

(2)

$$\text{BIC}\hspace{0.17em}=n\left[\text{ln}(\mathit{RSS})\right]+k\text{ln}(n)$$

(3)

where, k is the number of free parameters in the technique, n is the number of input data points, and RSS is the residual sum of squares between the mean NDVI data and the fitted technique. A lower value for AIC and BIC would indicate a preferable technique.

For M-BISE, MVI, and IDR required only one free parameter. For S-G and CW, two and three free parameters were needed, respectively. For the A-G and D-L techniques, seven and six parameters were needed, respectively. For WS, the free parameter is 9.37 when λ is 15 according to Atkinson et al. [45].

3. Results

3.1. Visual Assessment of the Fitted Curves

Figures 2 and 5 present the sample pixel results for the five vegetation types in the four NDVI data sets for the study area. The sample positions are shown in Figure 1.

Due to space limitations, only the results for one year (2000) are shown in the figure. It was observed that most de-noising techniques were effective for reconstructing high-quality NDVI time-series data sets. The results closely resemble the profiles of vegetation growth, especially for crops, grass, shrubs, and meadows. However, some obvious differences were found in the reconstructed NDVI time-series created with the eight techniques. For example, the WS technique produced the smoothest fitted curve, but it showed a lower maximum NDVI than other techniques (Figure 2a–d; Figure 4a,b; and Figure 5a). The S-G technique obtained most close the upper envelope for the input NDVI values in almost all cases. It is hard to judge the results for some techniques because the fitting curves results approximately match the other curves. Therefore, we quantify the differences between the eight techniques by using the RMSE, the AIC and the BIC in the following section.

3.2. Quantitative Evaluation and Regional Application

The RMSE, the AIC, and the BIC were calculated at the pixel level for the eight techniques for all of the data sets after reconstruction. All of the homogeneous pixel samples for each vegetation type were used to calculate the mean value of the three indexes to further analyze the fit of the techniques. To test the performance of the eight techniques at the regional scale, the RMSE was used as an example to display the differences between the techniques in the study area.

3.2.1. GIMMS NDVI

Table 3. RMSE, AIC, and BIC values for the eight techniques of Global Inventory Modeling and Map Studies (GIMMS) NDVI time-series data from January 1982 to December 2006 (best-fitting technique shown in bold).

Evaluation Indexa	Vegetation Type	Techniqueb
		A-G	M-BISE	CW	D-L	IDR	MVI	S-G	WS
RMSE	shrub	0.0071	0.0081	0.0055	0.0067	0.0104	0.0082	0.0055	0.0061
	crop	0.0123	0.0174	0.0093	0.0114	0.0214	0.0153	0.0087	0.0122
	desert	0.004	0.0045	0.0033	0.0038	0.0058	0.0047	0.0033	0.0031
	grassland	0.0128	0.0166	0.009	0.0118	0.0199	0.0152	0.0088	0.0126
	meadow	0.0162	0.0192	0.0101	0.0152	0.0238	0.0175	0.0102	0.0165
AIC	shrub	−2330	−2206	−2559	−1421	−2407	−2170	−2564	−2600
	crop	−1549	−1192	−1833	−802	−1653	−1262	−1908	−1619
	desert	−2889	−2791	−3111	−1968	−2929	−2728	−3080	−3216
	grassland	−1441	−1184	−1851	−759	−1545	−1247	−1873	−1493
	meadow	−1161	−996	−1711	−573	−1241	−1074	−1696	−1161
BIC	shrub	−6137	−6040	−6384	−6219	−5717	−6004	−6393	−6429
	crop	−5357	−5026	−5658	−5465	−4780	−5096	−5737	−5448
	desert	−6697	−6624	−6936	−6740	−6304	−6562	−6909	−7045
	grassland	−5249	−5018	−5676	−5357	−4800	−5081	−5702	−5322
	meadow	−4968	−4829	−5536	−5053	−4581	−4908	−5525	−4990

^aRoot mean square error (RMSE), the Akaike’s Information Criterion (AIC), and the Bayesian Information Criterion (BIC). The smaller the value, the better the fit;

^bTechniques: Asymmetric Gaussian (A-G) technique, the modified-best index slope extraction (M-BISE) technique, the changing-weight filter (CW) technique, the double logistic function (D-L), the iterative interpolation for data reconstruction (IDR) technique, the mean value iteration (MVI) technique, the Savitzky-Golay filter (S-G) technique, and the Whittaker smoother (WS) technique.

shows the RMSE, the AIC, and the BIC for each of the eight techniques for the five vegetation types of GIMMS NDVI time-series data.

The technique performance differs based on the three evaluation indexes. For crop, desert, grassland and meadow, the optimal techniques are S-G, WS, S-G, and CW, respectively, in terms of the three evaluation indexes. However, the best technique is different for shrub in terms of the three evaluation indexes, the optimal technique is S-G for RMSE and AIC, while it is WS for BIC. At the same time, the worst technique is also different in terms of the three evaluation indexes, for the five vegetation types both RMSE and BIC indicate that IDR perform worst, while it is D-L for AIC. For all the eight techniques, the accuracy differs for different vegetation types. For example, the WS technique performs better than the D-L technique for shrub and desert areas, but the situation is reversed for other three vegetation areas. Overall, the S-G, CW and WS techniques are effective for the GIMMS NDVI time series data. The IDR, M-BISE, and MVI techniques perform worse than other techniques for most vegetation types. Evaluation indexes can affect the assessment results for some vegetation types.

Figure 6 presents the RMSE values calculated for all pixels in the GIMMS NDVI data sets using the eight techniques. All of the techniques show good fitting results (i.e., low RMSE values). The performances of the eight techniques in Figure 6 are similar to the performances described in Table 3. Overall, the S-G (Figure 6g) and CW techniques (Figure 6c) performed best; most of the study areas are green in the results (representing lower RMSE values). The WS (Figure 6h) and D-L (Figure 6d) techniques were the next best, showing more yellow areas in Figure 6d,h. The M-BISE (Figure 6b) and IDR (Figure 6e) techniques performed worst, producing more red areas than the other techniques (representing higher RMSE values).

3.2.2. PAL NDVI

There are a few differences between the results for the GIMMS NDVI data and the PAL NDVI data (

Table 4. RMSE, AIC, and BIC values for the eight techniques for the Pathfinder AVHRR Land (PAL) NDVI time-series data from July 1981 to December 2000 (except August 1994–December 1994) (best-fitting technique shown in bold).

Evaluation index a	Vegetation type	Technique b
		A-G	M-BISE	CW	D-L	IDR	MVI	S-G	WS
RMSE	shrub	0.0143	0.019	0.0137	0.0131	0.0251	0.0197	0.0093	0.0103
	crop	0.0198	0.027	0.0152	0.0181	0.033	0.024	0.0119	0.0148
	desert	0.0066	0.0099	0.0073	0.0062	0.0139	0.0101	0.0052	0.0048
	grassland	0.0196	0.0269	0.0173	0.0182	0.0334	0.0249	0.0119	0.0149
	meadow	0.0228	0.0277	0.0178	0.0209	0.0331	0.026	0.0122	0.0174
AIC	shrub	−1528	−1128	−1545	−1656	−788	−1087	−2036	−2037
	crop	−950	−580	−1310	−1075	−316	−698	−1636	−1377
	desert	−2514	−2038	−2393	−2587	−1675	−2053	−2828	−3059
	grassland	−985	−574	−1147	−1087	−289	−662	−1648	−1392
	meadow	−777	−536	−1102	−900	−309	−606	−1605	−1174
BIC	shrub	−5962	−5589	−5997	−6094	−5249	−5547	−6492	−6493
	crop	−5383	−5041	−5761	−5513	−4777	−5159	−6092	−5833
	desert	−6948	−6499	−6844	−7025	−6136	−6514	−7284	−7515
	grassland	−5419	−5035	−5599	−5525	−4750	−5123	−6104	−5848
	meadow	−5211	−4997	−5554	−5338	−4770	−5067	−6061	−5630

^aRoot mean square error (RMSE), the Akaike’s Information Criterion (AIC), and the Bayesian Information Criterion (BIC). The smaller the value, the better the fit;

^bTechniques: Asymmetric Gaussian (A-G) technique, the modified-best index slope extraction (M-BISE) technique, the changing-weight filter (CW) technique, the double logistic function (D-L), the iterative interpolation for data reconstruction (IDR) technique, the mean value iteration (MVI) technique, the Savitzky-Golay filter (S-G) technique, and the Whittaker smoother (WS) technique.

). The three indexes show that the S-G technique performs best for crop, grassland and meadow, while WS technique performs best for desert. However, the best technique is different for shrub in terms of the three evaluation indexes, the optimal technique is S-G for RMSE, while it is WS for AIC and BIC. For all vegetation types, the S-G and WS techniques performs better than the A-G and the D-L techniques in terms of the three evaluation indexes. The IDR perform worse than the other techniques for all vegetation types. The eight techniques performed consistently for most situations according to the three evaluation indexes, but there were subtle differences between shrub and desert.

Figure 7 shows the RMSE calculated for each pixel of the PAL NDVI data sets using the eight techniques. Overall, the S-G (Figure 7g) and WS techniques (Figure 7h) perform best. The next best are the A-G (Figure 7a), the CW (Figure 7c), the D-L (Figure 7d), and the MVI (Figure 7f) techniques, and it is hard to determine which is better because their performances vary with different areas. Again, the IDR (Figure 7e) perform worst for this data set, and next to it is the M-BISE (Figure 7b).

3.2.3. SPOT VGT NDVI

For the SPOT VGT NDVI data, the S-G technique is also the best reconstruction technique for shrub, crop and meadow, and the WS technique outperforms the other six techniques for desert (

Table 5. RMSE, AIC, and BIC values for the eight techniques for SPOT VGT NDVI time-series data from April 1998 to April 2013 (best-fitting technique shown in bold).

Evaluation Indexes a	Vegetation Types	Technique b
		A-G	M-BISE	CW	D-L	IDR	MVI	S-G	WS
RMSE	shrub	0.0111	0.0129	0.0089	0.0111	0.0131	0.0121	0.007	0.0078
	crop	0.0133	0.0142	0.0091	0.0135	0.0144	0.0145	0.0084	0.0099
	desert	0.0043	0.005	0.0039	0.004	0.0057	0.0049	0.0033	0.003
	grassland	0.0128	0.0146	0.0101	0.0125	0.0147	0.0144	0.0083	0.0084
	meadow	0.014	0.0157	0.0108	0.0138	0.0159	0.0146	0.0086	0.0095
AIC	shrub	−1716	−1563	−1897	−1745	−1518	−1605	−2136	−2121
	crop	−1369	−1271	−1719	−1354	−1270	−1248	−1830	−1698
	desert	−2504	−2361	−2607	−2590	−2211	−2377	−2801	−2922
	grassland	−1376	−1238	−1605	−1408	−1235	−1253	−1828	−1839
	meadow	−1292	−1165	−1541	−1321	−1144	−1229	−1795	−1725
BIC	shrub	−5083	−4957	−5281	−5116	−4911	−4998	−5525	−5510
	crop	−4737	−4664	−5103	−4725	−4663	−4641	−5219	−5087
	desert	−5871	−5754	−5992	−5962	−5604	−5770	−6190	−6310
	grassland	−4744	−4631	−4989	−4779	−4628	−4646	−5217	−5228
	meadow	−4659	−4558	−4926	−4692	−4537	−4622	−5184	−5113

^aRoot mean square error (RMSE), the Akaike’s Information Criterion (AIC), and the Bayesian Information Criterion (BIC). The smaller the value, the better the fit;

^bTechniques: Asymmetric Gaussian (A-G) technique, the modified-best index slope extraction (M-BISE) technique, the changing-weight filter (CW) technique, the double logistic function (D-L), the iterative interpolation for data reconstruction (IDR) technique, the mean value iteration (MVI) technique, the Savitzky-Golay filter (S-G) technique, and the Whittaker smoother (WS) technique.

) in terms of three evaluation indexes. However, the best technique is different for grassland in terms of the three evaluation indexes, the optimal technique is S-G for RMSE, while it is WS for AIC and BIC. Overall, the CW and the WS techniques are effective for all the vegetation types, but it is difficult to determine which is better because their performances vary with different vegetation types. The three evaluation indexes are highly consistent for all vegetation types except subtle differences in grassland.

Figure 8 displays the RMSE calculated for each pixel of the SPOT VGT NDVI data set using the eight techniques. All of the techniques show good fitting results (i.e., low RMSE values). The S-G technique (Figure 8g) performs best, and the next best are the WS (Figure 8h) and CW techniques (Figure 8c). The performances of the A-G (Figure 8a) and the D-L techniques (Figure 8d) are next to the WS (Figure 8h) and CW techniques (Figure 8c). The performances of the M-BISE (Figure 8b), the IDR (Figure 8e), and the MVI techniques (Figure 8f) are relatively poor in the study area.

3.2.4. MODIS NDVI

As shown in

Table 6. RMSE, AIC, and BIC values for the eight techniques for the NDVI time-series data of Terra MOD13A2 from February 2000 to February 2013 (best-fitting technique shown in bold).

Evaluation Index a	Vegetation Type	Technique b
		A-G	M-BISE	CW	D-L	IDR	MVI	S-G	WS
RMSE	shrub	0.0138	0.0143	0.0119	0.0141	0.019	0.016	0.0089	0.0122
	crop	0.0182	0.0217	0.0182	0.0194	0.0294	0.0243	0.0125	0.0181
	desert	0.0043	0.0041	0.0032	0.0035	0.0044	0.0041	0.0028	0.0028
	grassland	0.0142	0.0181	0.0147	0.0147	0.0219	0.0195	0.0103	0.0136
	meadow	0.017	0.0175	0.0142	0.0169	0.0234	0.02	0.011	0.0156
AIC	shrub	−1454	−1466	−1600	−1538	−1394	−1455	−1687	−1654
	crop	−1648	−1730	−1884	−1786	−1694	−1737	−1923	−1950
	desert	−1301	−1309	−1460	−1385	−1241	−1293	−1558	−1507
	grassland	−1701	−1758	−1903	−1824	−1715	−1764	−1950	−1982
	meadow	−1701	−1769	−1910	−1826	−1722	−1773	−1952	−1978
BIC	shrub	−2365	−2221	−2255	−2326	−2202	−2270	−2406	−2283
	crop	−2095	−1990	−1994	−2070	−1885	−2072	−2254	−1999
	desert	−3091	−2755	−2804	−3008	−2801	−2779	−2861	−2941
	grassland	−2107	−2016	−2047	−2096	−1969	−2072	−2218	−2061
	meadow	−2111	−1939	−1943	−2085	−1873	−1972	−2095	−2020

^aRoot mean square error (RMSE), the Akaike’s Information Criterion (AIC), and the Bayesian Information Criterion (BIC). The smaller the value, the better the fit;

^bTechniques: Asymmetric Gaussian (A-G) technique, the modified-best index slope extraction (M-BISE) technique, the changing-weight filter (CW) technique, the double logistic function (D-L), the iterative interpolation for data reconstruction (IDR) technique, the mean value iteration (MVI) technique, the Savitzky-Golay filter (S-G) technique, and the Whittaker smoother (WS) technique.

, the S-G technique performs best for all the five vegetation types in terms of RMSE, while the best technique is different for the five vegetation types in terms of the AIC and BIC (shrub excepted). The optimal techniques for crop, desert, grassland and meadow are the WS, S-G, WS and WS in terms of AIC, respectively, but it is the S-G, A-G, S-G, and A-G in terms of BIC, respectively. Overall, both the S-G and the WS techniques are effective techniques for most vegetation types using the MODIS NDVI time-series data. The IDR technique performs worst for most vegetation types and the M-BISE and the MVI techniques perform second worst. The three evaluation indexes are inconsistent for most vegetation types.

Figure 9 shows the RMSE of each pixel of the MODIS NDVI data set using the eight techniques of the study area. All of the techniques show good fitting results (i.e., low RMSE values). The S-G (Figure 9g) technique performs best. The IDR (Figure 9e) and the MVI (Figure 9f) techniques perform worst in the study area. It is difficult to distinguish which one is better for the rest five techniques because the techniques are complementary to each other in different areas.

4. Discussion

4.1. Performance of the Reconstruction Techniques

Overall, the S-G, CW, and WS techniques show better reconstructed effects than the other techniques, and the IDR technique shows generally poor performance in terms of RMSE, AIC, and BIC for most vegetation types in the Heihe River Basin (

Table 7. Performance orders of the eight techniques for different data sources in terms of RMSE, AIC, and BIC.

Data Source	Vegetation Type	The Eight Techniques in Order of Fit a
GIMMS NDVI	shrub	WS	S-G	CW	D-L	A-G	M-BISE	MVI	IDR
	Crop	S-G	CW	D-L	WS	A-G	MVI	M-BISE	IDR
	desert	WS	CW	S-G	D-L	A-G	M-BISE	MVI	IDR
	grassland	S-G	CW	D-L	WS	A-G	MVI	M-BISE	IDR
	meadow	CW	S-G	D-L	WS	A-G	MVI	M-BISE	IDR
PAL NDVI	shrub	WS	S-G	D-L	CW	A-G	M-BISE	MVI	IDR
	Crop	S-G	WS	CW	D-L	A-G	MVI	M-BISE	IDR
	desert	WS	S-G	D-L	A-G	CW	MVI	M-BISE	IDR
	grassland	S-G	WS	CW	D-L	A-G	MVI	M-BISE	IDR
	meadow	S-G	WS	CW	D-L	A-G	MVI	M-BISE	IDR
SPOT VGT NDVI	shrub	S-G	WS	CW	D-L	A-G	MVI	M-BISE	IDR
	Crop	S-G	CW	WS	A-G	D-L	M-BISE	IDR	MVI
	desert	WS	S-G	CW	D-L	A-G	MVI	M-BISE	IDR
	grassland	WS	S-G	CW	D-L	A-G	MVI	M-BISE	IDR
	meadow	S-G	WS	CW	D-L	A-G	MVI	M-BISE	IDR
MODIS NDVI	shrub	S-G	WS	CW	A-G	D-L	M-BISE	MVI	IDR
	Crop	S-G	WS	CW	A-G	D-L	M-BISE	MVI	IDR
	desert	WS	S-G	CW	D-L	A-G	M-BISE	MVI	IDR
	grassland	S-G	WS	A-G	CW	D-L	M-BISE	MVI	IDR
	meadow	S-G	CW	WS	D-L	A-G	M-BISE	MVI	IDR

^aTechniques: Asymmetric Gaussian (A-G) technique, the modified-best index slope extraction (M-BISE) technique, the changing-weight filter (CW) technique, the double logistic function (D-L), the iterative interpolation for data reconstruction (IDR) technique, the mean value iteration (MVI) technique, the Savitzky-Golay filter (S-G) technique, and the Whittaker smoother (WS) technique. The techniques shown at the beginning of the row perform better than those at the end.

). These findings are inconsistent with these of Hird and McDermid [15] who stated that A-G and D-L techniques perform better than S-G in terms of RMSE. In our study, the A-G technique performs better than S-G only for some vegetation types (desert and meadow) of MODIS NDVI in terms of BIC (Table 6). However, the findings are supported by Zhu et al. [26] and Jiang et al. [27], and both of the results indicated that the S-G technique performs better than the A-G and D-L techniques in terms of RMSE under different noise levels. Though Julien and Sobrino [44] revealed that the IDR technique performs better than the HANTS and the D-L techniques in terms of the distance to the raw data and the proximity to its upper envelope, IDR performed worst among the eight techniques in most of our test cases.

Except for the S-G and the IDR techniques, the performance of other techniques is relatively unstable, meaning that they change for different vegetation types and data sources. For example, for SPOT VGT NDVI data sets, the CW and the WS are effective techniques for all vegetation types, performing worse than the S-G technique (desert and grassland excepted) and better than the rest five techniques (Tables 4 and 7). But for GIMMS NDVI data sets these techniques’ performance changed with vegetation types. The WS technique performed better than the D-L for shrub and desert vegetation, while the reverse is true for the other three vegetation types of MODIS NDVI. The D-L and the A-G techniques performed similarly to each other. The D-L technique outperforms the A-G technique for all vegetation types of GIMMS and PAL NDVI, but for MODIS and SPOT VGT NDVI, the situation is changed with vegetation types. Thus, it is difficult to judge which is better. This result is supported by Jönsson & Eklundh [35,46,68] who argued that the two techniques complement each other and that they may be suitable in different areas depending on the behavior of the NDVI signal. Beck et al. [33], and Hird and McDermid [15] also obtained similar findings. In general, the D-L and the A-G techniques outperform the MVI and the M-BISE techniques. The MVI and the M-BISE techniques show pool performance, which only outperform the IDR. It is hard to determine which is better because their performances vary with different vegetation types and data sources. However, for PAL NDVI data sets, the results are consistency with those of Ma and Veroustraete [25] who indicated that the MVI performed better than the M-BISE.

4.2. Effect of Evaluation Index on the Final Reconstruction Results

As ground reference measurements are challenging to obtain due to the medium/coarse resolution of the imagery, the problem of developing a robust, accurate and fast filter is amplified by the difficulty of obtaining reference measurements to use for validation [43]. Thus, choosing an appropriate index is very important for the final evaluation of de-noising techniques. To reduce the uncertainty associated with relying on one evaluation index, we used three different statistical indexes in this study.

Generally, the results of the three indexes are consistent in most cases, especially for worst-performing technique. For this technique, the assessment results for all vegetation types for all four NDVI data sets are in complete agreement (Tables 3 and 6). However, some differences were observed in the other techniques. For the GIMMS NDVI data, for example, the WS technique performance for shrub vegetation is third according to the RMSE, but, according to the AIC and BIC values, it is best. A similar situation was found for desert vegetation of SPOT VGT NDVI data. Different evaluation indexes have great effect on evaluation the performance of each technique for MODIS NDVI. The performance of the A-G and D-L techniques is different for all vegetation types in terms of RMSE, AIC, and BIC. A similar result was found for the CW and WS techniques. Thus, in certain situations, the evaluation results may be arbitrary if only one index is used. Comprehensive consideration of several indexes is very important.

4.3. Factors Influencing Performance

The performance is different between the eight techniques even under the same vegetation type and data source according to the above results (Tables 3 and 6). Thus, one of the main factors affecting the techniques’ performance is the de-noise principle of each techniques. The S-G applies an iterative weighted moving average filter to the NDVI time series, with the weight defined by a polynomial of a particular degree. This polynomial is designed to preserve higher moments within the data and to reduce the bias introduced by the filter [24]. This approach can replace the noise data as well as keep the fidelity of the time series. The results of this study also proved that the S-G technique is an effective de-noise method. The CW technique designed to replace the noise data by a three-point changing-weight filter while preserving the curve shape, as well as the timing and the amplitude of the local maxima/minima in the NDVI time series for a broad range of phenologies [26]. In addition, the WS is based on penalized least squares to fits a discrete series to discrete data and penalizes the roughness of the smooth curve. It balances reliability of the data and roughness of the fitted data. Both the CW and WS techniques also can replace the noise data as well as keep the fidelity of the original data. The two techniques were also proved effective techniques which only next to the S-G in this study. The A-G and D-L techniques are using a series of parameters to model the NDVI time series [33,35]. The difference is the A-G is based on asymmetric Gaussian functions while the D-L is based on a double logistic function to determine the parameters. That is may be the reason of the performance of the two techniques is similar and difficult to judge which is better. The M-BISE, MVI and IDR techniques show pool performance in this study. The M-BISE is to look for a spike (i.e., an increase immediately followed by a decrease or a decrease immediately followed by an increase) [21]. Both the increased range and decreased range greater than a threshold are treated as noise values. Therefore, the key to this approach is to set a suitable threshold. However, the noise data are changed with time and space. It is difficult to replace all noise data by a fixed threshold. The same problem was existed for the MVI and the IDR techniques. The MVI remove the noise data by replacing the maximum difference date value by the average of the dates before and after it. Iteration will stop when all differences are less than a threshold. The IDR is somewhat similar to the MVI, with the difference that the threshold of IDR technique is carried out from the data itself, and not from a comparison to an average of different years [44].

Vegetation type is also a factor affecting the techniques’ performance. For the five vegetation types in our study, the performances of the techniques for each type are unstable for most of the NDVI data. Similar findings were found by the research of Atkinson et al. [45] who compared the performances of A-G, D-L, FT, and WS for four vegetation types using the Terrestrial Chlorophyll Index (MTCI), the RMSE results indicated that the best technique changed with different vegetation types. Hird and McDermid [15] also showed that six techniques perform differently for six different vegetation regions. Another factor that had an influence on performance was the temporal and spatial resolution of the data. Since different data sources, as well as different vegetation types, the characteristic of noise varied. Therefore, Hird and McDermid [15] recommend that the strength and character of the noise present in an NDVI data set be considered when selecting an approach for time-series noise reduction.

As this study was focus on the Heihe River Basin, the vegetation types were relatively simple; most of the study area was covered by desert. Our findings and conclusions may not be as applicable to areas with multiple growing seasons (e.g., subtropical zones). However, our comparison of the performances of eight reconstruction techniques was conducted systematically and comprehensively using four NDVI data sets covering five different vegetation types, and we determined obtained the optimum technique for GIMMS, PAL, SPOT, and MODIS data for the different vegetation types. These findings will be of great reference and actual using values for choosing de-noise techniques and their parameter values.

5. Conclusions

In this study, the performances of eight de-noised techniques were compared for different vegetation types represented in four NDVI data sets in the Heihe River Basin, and the following results can be observed: the S-G, CW, and WS techniques perform better than other techniques for almost all vegetation types according to the RMSE, the AIC, and the BIC. The IDR, M-BISE, and MVI techniques performed worse than the other techniques for most vegetation types using the four sensor data sets. The best technique varies with vegetation types and NDVI data sources. However, the S-G performs best in most situations, the CW and WS techniques are next to it. The assessment results are consistent among the three evaluation indexes for most situations, but subtle differences exist for some vegetation types; and comprehensive consideration of several indexes is very helpful to decide the best technique for certain situations.