Reconstruction of Vegetation Index Time Series Based on Self-Weighting Function Fitting from Curve Features

Abstract

Vegetation index (VI) time series derived from satellite sensors have been widely used in the estimation of vegetation parameters, but the quality of VI time series is easily affected by clouds and poor atmospheric conditions. The function fitting method is a widely used effective noise reduction technique for VI time series, but it is vulnerable to noise. Thus, ancillary data about VI quality are utilized to alleviate the interference of noise. However, this approach is limited by the availability, accuracy, and application rules of ancillary data. In this paper, we aimed to develop a new reconstruction method that does not require ancillary data. Based on the assumptions that VI time series follow the gradual growth and decline pattern of vegetation dynamics, and that clouds or poor atmospheric conditions usually depress VI values, we proposed a reconstruction method for VI time series based on self-weighting function fitting from curve features (SWCF). SWCF consists of two major procedures: (1) determining a fitting weight for each VI point based on the curve features of the VI time series and (2) implementing the weighted function fitting to reconstruct the VI time series. The double logistic function, double Gaussian function, and polynomial function were tested in SWCF based on a simulated dataset. The results indicate that the weighted function fitting with SWCF outperformed the corresponding unweighted function fitting with the root-mean-square error (RMSE) significantly reduced by 26.82–52.44% (p < 0.05), and it also outperformed the Savitzky–Golay filtering with the RMSE significantly reduced by 13.98–54.04% (p < 0.05) for 270 sample points selected in mid-high latitudes of the Northern Hemisphere. Moreover, SWCF showed excellent robustness and applicability in regional applications.

1. Introduction

Vegetation index time series (e.g., normalized-difference vegetation index, NDVI [1]) derived from satellite sensors have been widely used in the estimation of land surface vegetation parameters (e.g., fractional vegetation coverage, leaf area index, vegetation phenology, and vegetation productivity) and the monitoring of land cover and land use change [2,3,4,5,6]. However, the quality of vegetation index (VI) time series is easily affected by the disturbances caused by cloud contamination and atmospheric conditions [7]. Although the maximum value composition [8] and the cloud detection identification algorithm [9] are often applied for the preprocessing of VI products, significant residual noise still remains [10]. Therefore, a number of noise-reduction techniques have been developed to reduce noise and reconstruct high-quality VI time series.

According to the reconstruction principle, the noise-reduction techniques can be divided into four categories [11]: (1) temporal filtering methods, such as Savitzky–Golay filtering [12], mean-value iterative filtering (MVI) [10], moving average filtering (MA) [13], changing-weight filtering [14], best index slope extraction (BISE) [15], and modified BISE [16]; (2) temporal function fitting methods, such as the double logistic function fitting [17], double Gaussian function fitting [18], and polynomial function fitting [19]; (3) frequency-based methods, such as the harmonic analysis of time series (HANTS) [20], fast Fourier transform (FFT) [21], and wavelet transform (WT) [22]; and (4) spatiotemporal methods, such as the spatial–temporal Savitzky–Golay (STSG) method [23]. Although these methods have been successfully applied, certain drawbacks greatly limit their applications. For the frequency-based methods, they assume that the high frequency components are usually noise. They may reduce some reasonable high values when there are frequent variations, so they will fail to maintain vegetation signals completely. For the spatiotemporal methods, for example, STSG uses spatially adjacent pixels to restore the VI of target pixels, but it will fail when large spatially continuous clouds occur. This case happens often on high spatial resolution images, thus greatly limiting the application scope of STSG. For the temporal filtering methods and some frequency-based methods, there are usually no objective rules used for setting the key parameters such as the width of the smoothing window and the degree of the smoothing polynomial for the Savitzky–Golay filtering, the sliding period and acceptable threshold value for BISE, the size of moving window for MA, the key threshold value for MVI, and the control parameters for HANTS. Researchers need to empirically determine the reasonable values for these key parameters according to their study area, which may result in new uncertainties.

For the temporal function fitting methods, their key parameters can be automatically optimized by iterative nonlinear least square fitting, which has less dependence on manual settings. Moreover, the temporal function fitting methods have the advantage of expressing the VI time series as a function, which helps to interpret vegetation signals. However, the function fitting methods are vulnerable to noisy VI points contaminated by clouds and poor atmospheric conditions. Therefore, ancillary data about VI quality are often used to alleviate the interference of noise before fitting. Ancillary data are used to (1) replace the VI points identified as clouds in ancillary data with the linear interpolation values of their neighboring VI points [3,24] and (2) set a fitting weight for each VI point according to the signal-to-noise ratio indicated in ancillary data [25,26,27]. These approaches can effectively improve the performance of the function fitting methods, but meanwhile, are greatly limited by the availability, accuracy, and application rules of ancillary data. Ancillary data are generated by some noise detection algorithms and have limited accuracy [12,25,28]. In particular, the accuracy of ancillary data is not sufficiently high when the background is ice or snow in identifying clouds. For example, the accuracy for the C code based on the Function of Mask algorithm (CFMASK) in identifying clouds is 64.09% in this case [29]. Moreover, some VI datasets may lack available ancillary data about VI quality, such as the GIMMS NDVI dataset [30]. Furthermore, it varies greatly in the marker rules of VI quality for different ancillary data, for which it is difficult to set a uniform weighting method for different datasets [25,26].

In fact, the VI time-series data itself contains the quality information of each data point. That is, good-quality VI points usually follow the gradual process of vegetation dynamics, while suddenly dropping VI points are more likely to be contaminated by clouds and poor atmospheric conditions [12,31]. This pattern provides a new idea to estimate the signal-to-noise ratio of each VI data point, thus avoiding the limitations of ancillary data. Therefore, this study proposes a reconstruction method based on self-weighting function fitting from curve features (SWCF) to reduce noise in VI time series. SWCF utilizes the curve features of VI time series to generate the weight for each VI data point and then implements the weighted function fitting, which is expected to reconstruct high-quality VI time series without using ancillary data.

2. Materials and Methods

2.1. Data Sources

The MOD09GA (version 6) dataset was used to test the newly proposed SWCF. It included daily 500 m surface reflectance data and a 1 km cloud flag as ancillary data. The red and near-infrared reflectance bands were used to calculate NDVI, and the cloud flags were utilized to indicate whether the pixel was contaminated by clouds. NDVI and cloud flag data were used for the simulation experiment in the pixel-level evaluation (Section 2.3.1) and the region-level evaluation (Section 2.3.2). The MOD09A1 (Version 6) product, the 8 day maximum synthetic product derived from the MOD09GA dataset, was used to obtain the demonstrated NDVI time series curve in Section 2.2.1 for the clear introduction to the new method. The demonstrated NDVI time series curve consists of 46 VI points for one year.

The MCD12Q1 (version 6, 500 m) data product provided the annual global land covers identified as the IGBP (International Geosphere Biosphere Programme) classification. IGBP classes for the period of 2004–2019 were used for the selection of vegetation samples. These two datasets are both acquired from the NASA’s Land Processes Distributed Active Archive Center (LP DAAC, https://lpdaac.usgs.gov, accessed on 1 August 2021).

2.2. Algorithm of SWCF

Ideally, a noise-free VI time series changes gradually and continuously, which is consistent with the vegetation dynamics of growth and decline [12,31]. Therefore, SWCF is based on two assumptions: (1) that a yearly VI time series primarily indicates the seasonal dynamics of vegetation and (2) that clouds and poor atmospheric conditions usually tend to depress VI values, causing sudden drops in VI time series. SWCF utilizes the curve features of VI time series to determine a fitting weight for each VI data point and then implements the weighted function fitting to reconstruct the VI time series. The double logistic function, double Gaussian function, and polynomial function (Appendix B) are used for fitting in SWCF. Apparently, the core point in SWCF is how to obtain an appropriate fitting weight indicating the confidence level for each VI data point.

2.2.1. Self-Weighting from the VI Curve Features

SWCF exploits the signal-to-noise ratio of each data point in the VI time series to generate its weight (Figure 1). In line with the two assumptions mentioned above, during one vegetation growing season, VI points consistent with the vegetation variation pattern usually have higher signal-to-noise ratios, while suddenly dropping VI points incompatible with the gradual process are inclined to be degraded by noise to varying degrees, owning relatively lower signal-to-noise ratios. The weight is introduced as an indicator of the signal-to-noise ratio, with a higher signal-to-noise ratio leading to a higher weight.

(1)

Definition of gradually changing points and suddenly dropping points

Following the growth and decline of vegetation during one growing season, the VI temporal profiles will rise and fall monotonically, while clouds and poor atmospheric conditions will cause drops in VI time series, thereby disturbing the trend. Based on this VI temporal change pattern, VI points rising monotonically from both ends to the peak are defined as gradually changing points, while others are defined as suddenly dropping points (Figure 2). The gradually changing points tend to contain actual vegetation signals and the suddenly dropping points are more likely to be disturbed by noise to varying extents.

The specific algorithm is described as follows. As shown in Figure 2, the VI time series is divided into the rising and falling parts (green curve and brown curve, respectively) by the VI peak point. In the rising part, starting from the first point to the peak point, if the VI value is larger than or equal to the maximum VI value of all previous points, then it will be defined as the gradually changing point; otherwise, it will be defined as the suddenly dropping point. Similarly, in the falling part, the gradually changing points and the suddenly dropping points are defined starting from the last point to the peak point. In Figure 2, the light blue solid points are cloud-contaminated VI marked in ancillary data, and DOY denotes the Julian day of year.

(2)

Calculation of the weights for the VI time series

The weight is set to 1 for the gradually changing point (Equation (1)).

$$W_{\mathrm{GCP}}=1$$

(1)

where W_GCP is the weight of the gradually changing point.

For the suddenly dropping point, the weight is determined by two features: one is the magnitude of the drop (Δh), and the other is the temporal proximity to the VI peak point (P) (Figure 3). The weight of the suddenly dropping point is determined by the following four steps:

(I) Linear stretching of the original VI time series. Since the data value range of the VI time series differs greatly in pixels, the linear stretching is adopted to uniformly transform the original VI time series to the same data value range to ensure that various curves with different amplitudes can be processed automatically (Equation (2)).

$$V{I}_{\mathrm{LS} }=\frac{VI-V{I}_{\min }}{V{I}_{\max }-V{I}_{\min }}\times R$$

(2)

where VI_LS denotes the stretched VI value, VI is the original VI value, R (Range) is the maximum value of VI_LS, and VI_max and VI_min indicate the maximum and minimum values of the original VI time series, respectively. The optimal stretching range R will be determined by experimental tests (see details in Section 2.3.3).

(II) Definition of Δh. Δh is defined as the vertical distance to the straight line formed by two adjacent stretched gradually changing points (Equation (3)). Δh reflects the deviation of a suddenly dropping point from the vegetation seasonal change trend and is used to estimate the signal-to-noise ratio. As shown in Figure 3, A and B are suddenly dropping points (red solid points); Δh₁ and Δh₂ are their vertical distances to the straight line formed by two adjacent stretched gradually changing points, respectively.

$$\mathsf{\Delta }h=\left(V{I}_{\mathrm{LS}\_\mathrm{GCP}1}+\left(DO{Y}_{\mathrm{SDP}}-DO{Y}_{\mathrm{GCP}1}\right)\times \left(\frac{V{I}_{\mathrm{LS}\_\mathrm{GCP}2}-V{I}_{\mathrm{LS}\_\mathrm{GCP}1}}{DO{Y}_{\mathrm{GCP}2}-DO{Y}_{\mathrm{GCP}1}}\right)\right)-V{I}_{\mathrm{LS}\_\mathrm{SDP}}$$

(3)

where Δh is the dropping magnitude of the suddenly dropping point; VI_LS__SDP, VI_LS__GCP1 and VI_LS__GCP2 are the VI_LS values of the suddenly dropping point and its adjacent gradually changing points; and DOY_SDP, DOY_GCP1 and DOY_GCP2 are the DOYs of the suddenly dropping point and its adjacent gradually changing points. For the VI_LS data points that span to the next year, their DOYs should be added with 365. Generally, a larger Δh indicates a lower signal-to-noise ratio.

(III) Definition of P. P shows the temporal proximity to the VI_LS peak, which can be calculated by their DOYs (Equation (4)). As shown in Figure 3, P₁ and P₂ represent the temporal proximity of A and B to M, respectively; M is the peak VI_LS point; DOY_A, DOY_B and DOY_M are the DOYs of A, B, and M, respectively; and DOY_min and DOY_max are the DOYs of the first and last VI points, respectively.

$$P=\{\begin{array}{c}\frac{DO{Y}_{\mathrm{SDP}}-DO{Y}_{\min }}{DO{Y}_{\mathrm{m}}-DO{Y}_{\min }}, \mathrm{when} DO{Y}_{\mathrm{SDP}}<DO{Y}_{\mathrm{m}}\\ \frac{DO{Y}_{\max }-DO{Y}_{\mathrm{SDP}}}{DO{Y}_{\max }-DO{Y}_{\mathrm{m}}}, \mathrm{when} DO{Y}_{\mathrm{SDP}}>DO{Y}_{\mathrm{m}}\end{array}$$

(4)

where P is the temporal proximity of the suddenly dropping point to the peak point; DOY_SDP and DOY_m are the DOYs of the suddenly dropping point and peak VI_LS point, respectively; and DOY_min and DOY_max represent the minimum and maximum values of the DOY series, respectively. For the VI data points that span to the next year, their DOYs should be added with 365. Obviously, P ranges from 0 to 1, and a larger P indicates that the suddenly dropping point is closer to the peak of vegetation growth.

(IV) Calculation of the weight for the suddenly dropping point. For the suddenly dropping point, higher noise levels and closer temporal proximity to the peak of growth lead to relatively lower weights. Since Δh and P are two independent variables, we take the multiplication combination of these two to comprehensively calculate the weights in Equation (5):

$$W_{\mathrm{SDP}}=\{\begin{array}{c}\hfill 1-\mathsf{\Delta }h\times P, \mathsf{\Delta }h\times P<1\\ \hfill 0, \mathsf{\Delta }h\times P\ge 1\end{array}$$

(5)

where W_SDP is the weight of the suddenly dropping point; Δh is the dropping magnitude of suddenly dropping point; and P is the temporal proximity of the suddenly dropping point to the peak VI point. Figure 4 shows the weight series for the VI time series.

2.2.2. Weighted Function Fitting

The double logistic function, double Gaussian function, and polynomial function (Appendix B) are used for fitting in SWCF. The parameters for these functions are optimized with the weighted least squares fitting process. It should be noted that SWCF emphasizes the importance of vegetation signals and lowers the weights of contaminated data points, which distinguishes SWCF from the unweighted function fitting methods.

2.2.3. Software Implementations

In this study, the self-weighting method and the function fitting methods are implemented in Google Earth Engine, and the code is available on the website provided in Appendix C. Other common softwares for remote sensing processing (e.g., Python, Matlab and RStudio) are also convenient for the implementation of SWCF.

2.3. Evaluation of SWCF

2.3.1. Evaluation at the Pixel Level

Similar to the evaluation techniques for noise-reduction methods [32,33,34], a simulation experiment was conducted to compare SWCF with other noise-reduction methods. At each test pixel, a reference daily NDVI time series was constructed from the average of multiyear noise-free daily NDVI data. Then, simulated noise would be added to the reference NDVI time series to test the performance of all methods.

(1)

Preselection of test pixels

According to the MCD12Q1 IGBP land cover types, we randomly selected 50 test pixels with a constant IGBP class during 2004–2019 for each land cover type. These types included evergreen needleleaf forests, deciduous needleleaf forests, deciduous broadleaf forests, mixed forests, shrublands, savannas, and grasslands across mid-high latitudes of the Northern Hemisphere.

(2)

Construction of the reference NDVI curve

A. Division of growing season

Since the annual vegetation cycle may not coincide with the calendar year, that is, the growing season may span two years, the actual start date of a growing season must be determined. To highlight the seasonal changes in vegetation, the 8-day maximum value composition was performed on daily NDVI time series to reduce noise levels, and then the 7-point Gaussian filtering was applied to generate a long-term NDVI time series in 2004–2019. Two adjacent growing seasons can be divided by the DOY corresponding to the lowest NDVI point in every year. As shown in Figure 5a, the green points are the minimum NDVI values between adjacent growing seasons from 2004 to 2007, and Date₁, Date₂, and Date₃ are the corresponding dates, respectively. Here, we took the average multiyear DOYs in 2004–2019 as the average start date of all growing seasons.

B. Construction of reference NDVI time series

The reference daily NDVI time series was constructed from the temporal average of multiyear cloud-free NDVI values on the same date in all growing seasons. Since the cloud flag data cannot include all cases in which NDVI data points are contaminated, the 7-point Gaussian filtering was applied to remove potential noisy data and produce a smoother curve (Figure 5b).

(3)

Classifications of the reference NDVI time series

When the environmental conditions differ a lot, even for the same land cover type, the curve shapes of NDVI time series would be quite different [35]. To comprehensively evaluate the performance of SWCF, here we defined 9 curve types (see

Table 1. Classification of the reference NDVI time-series curves.

Amplitude	Width
	Width ≤ 210 d	210 d < Width ≤ 260 d	Width > 260 d
Amplitude ≤ 0.3	A1W1	A1W2	A1W3
0.3 < Amplitude ≤ 0.5	A2W1	A2W2	A2W3
Amplitude > 0.5	A3W1	A3W2	A3W3

) based on the width and amplitude of a growing season derived from the reference NDVI time series (Appendix A). As shown in Table 1, A₁W₁, A₁W₂, A₁W₃, A₂W₁, A₂W₂, A₂W₃, A₃W₁, A₃W₂ and A₃W₃ are curve types of the reference NDVI time series, classified by the width and amplitude of a growing season. The threshold values of 210 and 260 are the 33rd and 66th percentile values of the widths of growing season derived from the reference NDVI time series for all preselected sample points, respectively, and the threshold values of 0.3 and 0.5 are the 33rd and 66th percentile values of the amplitudes of growing season derived from the reference NDVI time series for all preselected sample points, respectively. These nine curve types included almost all possible curve shapes of any land cover type and explicitly reflected the vegetation growth states. Then 30 samples were randomly selected for each curve type from previous preselected test pixels (Figure 6). These selected samples were used to evaluate the performance of different reconstruction methods.

(4)

Construction of the noise NDVI time series

A noise NDVI time series was generated by introducing random noise to a reference NDVI curve. At each test pixel, we simulated a group of random noise according to the actual occurrence of noise. Then, they were added to reduce the NDVI values of the reference NDVI time series based on the reality that clouds and atmospheric interference tend to degrade NDVI.

The noise simulation experiment consisted of the determinations of the number, location, and intensity of noise. The number depended on the average of annual cloudy days in 2004–2019. In a descending order of cloud probability per day, the corresponding DOYs were selected as the noise introduction locations (Figure 7a). That is, for a day, the higher probability of cloud occurrence made it more likely to be selected. For an NDVI temporal profile of the same day over multiple years, we assumed that its level of noise followed the normal distribution. Therefore, the noise intensity should be within 0–3 times its standard deviation. At each location, a random number within 0–3 times its standard deviation was independently generated as the noise intensity.

At each noise introduction location, we subtracted the noise intensity from the NDVI value of the reference NDVI time series (Equation (6)), thus generating a noise NDVI time-series (Figure 7b).

$$NDVI\_nois{e}_i=NDVI\_re{f}_i-Nois{e}_i$$

(6)

where NDVI_noise and NDVI_ref represent the noise NDVI time series and the reference NDVI time series, respectively. Noise is the simulated noise, and i is the DOY of the noise introduction location.

(5)

Reconstruction of the noise NDVI time series

At each test pixel, the noise NDVI time series was reconstructed by the weighted function fitting methods with SWCF, unweighted function fitting methods and Savitzky–Golay filtering, respectively (Appendix B).

(6)

Evaluation of the reconstructed NDVI time series

At each test pixel, the qualitative and quantitative analyses were carried out to evaluate the reconstructed NDVI time series. First, we gave a visual inspection of the overlap degrees between the reference NDVI time series and the reconstructed NDVI time series. Then, the root-mean-square error (RMSE) was used to quantitatively evaluate the performance of each method. The RMSE reflects the error of the reconstructed time series relative to the reference time series (Equation (7)).

$$\mathrm{RMSE}=\sqrt{\frac{{{\displaystyle \sum }}_{i=1}^{365}{\left(NDVI\_re{f}_i-NDVI\_fi{t}_i\right)}^2}{365}}$$

(7)

where RMSE is the root-mean-square error of the reconstructed NDVI time series, NDVI_ref and NDVI_fit represent the reference NDVI time series and the reconstructed NDVI time series, respectively, and i is the ordinal number of points in the yearly NDVI time series. A smaller RMSE indicates a better capacity to restore the noisy NDVI to its true value. For each curve shape type (Table 1), the paired sample t-test is used to analyze whether the RMSE values for these methods are significantly different.

2.3.2. Evaluation at the Region Level

To test the applicability of SWCF applied in large regions, we selected a test area (the geographic center is at 96.3°E, 33.8°N) on the Tibetan Plateau which is easily affected by clouds but rarely disturbed by human beings. The test area includes 782 × 345 MOD09GA 500 m pixels, and its major land cover type is natural grasslands according to IGBP classifications. The NDVI image on 5 August 2018 was considerably contaminated by clouds. For comparison, the weighted double logistic function fitting with SWCF and the unweighted double logistic function fitting were applied to reconstruct this NDVI image, respectively. Then, we visually inspected their reconstructed NDVI images by referencing the average cloud-free NDVI images on the same day in 2017 and 2019.

2.3.3. Determination of the Optimal Stretching Range for SWCF

In SWCF, the VI time series should be linearly stretched to a uniform range (0–R). The optimal R was determined according to the performance of SWCF at different values of R. Similar to the evaluation at the pixel level, we selected a new batch of samples (Figure 8) and constructed a reference NDVI time series and a noise NDVI time series for each of them. For all test pixels, when R varied from 0–30, the reduction rates of RMSE values (Equation (8)) were calculated for the weighted function fitting with SWCF relative to the unweighted function fitting. Obviously, the maximum reduction rate indicated the optimal R for SWCF.

$$D=\frac{RMS{E}_{\mathrm{un}}-RMS{E}_{\mathrm{w}\_\mathrm{R}}}{RMS{E}_{\mathrm{un}}}\times 100$$

(8)

where D is the reduction rate of RMSE_{w_R}, RMSE_un denotes the RMSE for unweighted function fitting, and RMSE_{w_R} represents the RMSE for weighted function fitting with SWCF.

3. Results

3.1. The Optimal Stretching Range of VI Time Series

Figure 9 shows the reduction rates of RMSE values when R varies from 2 to 28 in steps of 2. When R varied from 2 to 10, the RMSE reduction rates for the weighted double logistic function fitting, double Gaussian function fitting, and polynomial function fitting with SWCF showed obvious increases, which indicated that the reconstruction accuracy improved rapidly. When R varied from 10 to 28, the RMSE reduction rates for the weighted double logistic function fitting and polynomial function fitting with SWCF decreased gradually, while those for the weighted double Gaussian function fitting with SWCF did not change obviously. Overall, the weighted function fitting with SWCF showed the best performance when R was equal to 10. Therefore, the optimal stretching range of VI time series is set to 10 for SWCF.

3.2. The Reconstruction Performance at the Pixel Level

Figure 10 shows the NDVI time series reconstructed by the weighted double logistic function fitting with SWCF. For the 9 curve shape types, the reference NDVI time series and the reconstructed time series by SWCF had the highest degree of overlap. In particular, they were quite close in width and amplitude. SWCF could retain the noise-free NDVI well and restore most noisy NDVIs to their true values. In contrast, the unweighted double logistic function fitting and Savitzky–Golay filtering were greatly affected by noise, and their reconstructed NDVI time series showed a large displacement from the reference NDVI time series. In addition, the Savitzky–Golay filtering may introduce some new noise when reconstructing noise NDVI time series with a high noise level (e.g., the curve shape types A₁W₂ and A₁W₃ in Figure 10).

As we can see from

Table 2. Comparison of Reconstruction Quality of Different Methods.

Curve Shape Types	Double Logistic Function Fitting				Polynomial Function Fitting				Double Gaussian Function Fitting				Savitzky–Golay Filtering
	Unweighted	Weighted	RMSE	RMSE	Unweighted	Weighted	RMSE	RMSE	Unweighted	Weighted	RMSE	RMSE
	RMSE	RMSE	Reduction	Reduction	RMSE	RMSE	Reduction	Reduction	RMSE	RMSE	Reduction	Reduction	RMSE Average
	Average	Average	Rate * (%)	Rate # (%)	Average	Average	Rate * (%)	Rate # (%)	Average	Average	Rate * (%)	Rate # (%)
A1W1	0.0574	0.0273	52.44	54.04	0.0601	0.0290	51.75	51.18	0.0601	0.0304	49.42	48.82	0.0594
A1W2	0.0662	0.0350	47.13	43.55	0.0630	0.0335	46.83	45.97	0.0652	0.0392	39.88	36.77	0.0620
A1W3	0.0667	0.0402	39.73	37.48	0.0629	0.0325	48.33	49.46	0.0779	0.0416	46.60	35.30	0.0643
A2W1	0.0922	0.0496	46.20	44.02	0.0969	0.0610	37.05	31.15	0.1023	0.0678	33.72	23.48	0.0886
A2W2	0.1051	0.0593	43.58	41.05	0.1058	0.0614	41.97	38.97	0.1144	0.0685	40.12	31.91	0.1006
A2W3	0.1157	0.0703	39.24	38.01	0.1136	0.0643	43.40	43.30	0.1177	0.0728	38.15	35.80	0.1134
A3W1	0.1022	0.0675	33.95	26.87	0.1085	0.0794	26.82	13.98	0.1090	0.0767	29.63	16.90	0.0923
A3W2	0.1072	0.0689	35.73	29.91	0.1073	0.0704	34.39	28.38	0.1127	0.0784	30.43	20.24	0.0983
A3W3	0.1196	0.0771	35.54	27.74	0.1122	0.0687	38.77	35.61	0.1168	0.0784	32.88	26.52	0.1067

Thirty sample points were selected for each curve type. For each sample point, the VI time series curve contained 365 VI points. * represents the RMSE reduction rate (%) of the weighted fitting reconstruction results relative to the unweighted fitting ones for the same curve shape type and function. ^# represents the RMSE reduction rate (%) of the weighted fitting reconstruction results relative to the Savitzky–Golay filtering method for the same curve shape type. A₁W₁, A₁W₂, A₁W₃, A₂W₁, A₂W₂, A₂W₃, A₃W₁, A₃W₂, and A₃W₃ refer to different curve shape types of NDVI time series (Table 1).

, for the nine curve shape types, the RMSE values for the weighted double logistic function fitting, polynomial function fitting and double Gaussian function fitting with SWCF were all significantly reduced (p < 0.05) by 33.95–52.44%, 29.63–49.42%, and 26.82–51.75% compared with those for the corresponding unweighted ones, respectively, and also significantly reduced (p < 0.05) by 26.87–54.04%, 16.90–48.82%, and 13.98–51.18% compared with those for the Savitzky–Golay filtering, respectively. On the whole, the weighted double logistic function fitting with SWCF showed the best performance on reconstructing NDVI time series of various shapes.

3.3. The Reconstruction Performance at the Regional Level

In comparison with the unweighted double logistic function fitting, the weighted double logistic function fitting with SWCF reconstructed a NDVI image closer to that of the same day in the adjacent years (Figure 11). Obviously, the weighted double logistic function fitting with SWCF could restore most cloud-contaminated NDVIs to their actual values, showing great robustness and effectiveness in regional applications (Appendix C).

4. Discussion

4.1. The Advantages and Potential Applications of SWCF

SWCF is free from ancillary data. SWCF exploits the signal-to-noise ratio information contained in the VI time series itself to generate weight series for fitting. Distinguished from the existing weighting methods [25,26,27], SWCF is free from auxiliary data, thus avoiding potential uncertainties in ancillary data.

SWCF has a strong self-adaptive capacity. In the pixel-level evaluation, SWCF is applicable for 270 selected samples containing all kinds of vegetation types across the mid-high latitudes of the Northern Hemisphere. In addition, the weighted double logistic function fitting with SWCF shows excellent robustness and applicabilities in a grassland test area on the Tibetan Plateau including 782 × 345,500 m pixels from MOD09GA data.

SWCF can reconstruct high-quality VI time series. Evaluation at the pixel level indicates that the weighted function fitting with SWCF outperforms the corresponding unweighted function fitting with the RMSE significantly reduced by 26.82–52.44% (p < 0.05), and it also outperforms the Savitzky–Golay filtering with the RMSE significantly reduced by 13.98–54.04% (p < 0.05). In the regional evaluation, compared with the unweighted double logistic function fitting, the weighted double logistic function fitting with SWCF shows better performance in preserving the original noise-free NDVIs and restoring most cloud-contaminated NDVIs to their actual values.

SWCF can contribute in VI-based studies. The VI time series are important data sources for the estimation of land surface vegetation parameters (e.g., fractional vegetation cover, leaf area index, vegetation phenology, and vegetation productivity) and the monitoring of land cover and land use change [2,3,4,5,6]. SWCF is able to generate high-quality VI time series, thus having great application potential in VI-based studies mentioned above. For example, the cloud-contaminated VI will result in the underestimation of fractional vegetation cover when using the pixel dichotomy model [36]. SWCF can be used to reconstruct the VI time series and restore the cloud-contaminated VI to the actual value, thus contributing in the accurate estimation of fractional vegetation cover.

4.2. Limitations and Prospects

SWCF may have difficulty in fitting VI time series with an irregular or inconspicuous temporal profile due to the limitations of function models. Moreover, proper initial values for function parameters are necessary in nonlinear function fitting [26,27]. In regional applications, a set of fixed initial parameter values may not be applicable for reconstructing VI time series with different curve shapes.

SWCF is designed based on a single vegetation growth cycle. Therefore, it is necessary to first divide the long-term VI time series into each single vegetation growth cycle before applying SWCF. For example, when reconstructing a whole year VI time series with the vegetation growth cycle spanning two years, the VI time series for this year and the adjacent two years should be used to cover a whole vegetation growth cycle. Then, SWCF can be applied for each complete vegetation growth cycle. Finally, we can extract the VI time series of this year from the reconstructed results. Similarly, this strategy can also be applied for the reconstruction of VI time series with multiple growth cycles in a single year.

SWCF assumes that higher VI values usually have higher signal-to-noise ratios, which may sometimes be invalid. For example, satellite data transmission errors may cause abnormally high NDVI values [15], and excessive atmospheric correction may also lead to spurious EVI high values in cloud-contaminated pixels [37]. Therefore, some priori knowledge is required to preprocess the VI time-series data before applying SWCF. For example, an NDVI value that increases more than 0.4 in 20 days should be removed because such a large increase is unlikely to be caused by vegetation changes [12], and spurious EVI values larger than NDVI should be rejected because EVI should be smaller than NDVI under normal conditions [38]. These approaches can help to identify high-value noisy VIs, thus ensuring the effectiveness of applying SWCF.

5. Conclusions

This study proposed a new reconstruction method based on self-weighting function fitting from curve features (SWCF) to reduce noise in the vegetation index time series. SWCF is free from ancillary data. For the vegetation index time series of various shapes, the weighted function fitting with SWCF reduced the root-mean-square error significantly by 26.82–52.47% (p < 0.05) compared with the corresponding unweighted function fitting methods and reduced the root-mean-square error significantly by 13.98–54.04% (p < 0.05) compared with the Savitzky–Golay filtering method. Moreover, SWCF was strongly robust and effective in regional applications. It showed an excellent denoising and self-adaptive capacity and could be applicable to reconstruct high-quality vegetation index time series.