Using Enhanced Gap-Filling and Whittaker Smoothing to Reconstruct High Spatiotemporal Resolution NDVI Time Series Based on Landsat 8, Sentinel-2, and MODIS Imagery

Abstract

Normalized difference vegetation index (NDVI) time series data, derived from optical images, play a crucial role for crop mapping and growth monitoring. Nevertheless, optical images frequently exhibit spatial and temporal discontinuities due to cloudy and rainy weather conditions. Existing algorithms for reconstructing NDVI time series using multi-source remote sensing data still face several challenges. In this study, we proposed a novel method, an enhanced gap-filling and Whittaker smoothing (EGF-WS), to reconstruct NDVI time series (EGF-NDVI) using Google Earth Engine. In EGF-WS, NDVI calculated from MODIS, Landsat-8, and Sentinel-2 satellites were combined to generate high-resolution and continuous NDVI time series data. The MODIS NDVI was employed as reference data to fill missing pixels in the Sentinel–Landsat NDVI (SL-NDVI) using the gap-filling method. Subsequently, the filled NDVI was smoothed using a Whittaker smoothing filter to reduce residual noise in the SL-NDVI time series. With reference to the all-round performance assessment (APA) metrics, the performance of EGF-WS was compared with the conventional gap-filling and Savitzky–Golay filter approach (GF-SG) in Fusui County of Guangxi Zhuang Autonomous Region. The experimental results have demonstrated that the EGF-WS can capture more accurate spatial details compared with GF-SG. Moreover, EGF-NDVI of Fusui County exhibited a low root mean square error (RMSE) and a high coefficient of determination (R²). In conclusion, EGF-WS holds significant promise in providing NDVI time series images with a spatial resolution of 10 m and a temporal resolution of 8 days, thereby benefiting crop mapping, land use change monitoring, and various ecosystems, among other applications.

1. Introduction

The normalized difference vegetation index (NDVI) is one of the most commonly used vegetation indexes. In particular, it is important for various research areas such as crop mapping [1,2], vegetation monitoring [3,4], phenology extracting [5,6], drought monitoring [7,8], and land cover change monitoring [9]. Numerous studies have employed NDVI time series data derived from MODIS [10], Landsat-8 [11], and Sentinel-2 [12] images, however, the calculated NDVI time series may display spatial and temporal discontinuities due to the limitations of clouds and individual sensor constraints [13]. To address this challenge, numerous researchers have endeavored to employ time series reconstruction and fusion techniques involving remote sensing images from multiple sources, thereby providing valuable data for remote sensing investigations.

Over the past decades, many novel methods have been developed to reconstruct spatiotemporally continuous NDVI time series data [14]. These methods can be classified into four types: (1) temporal-based methods [15], (2) frequency-based methods [16], (3) hybrid methods [17], and (4) multi-source fusion methods [18]. This paper focuses on temporal-based methods, due to their widespread popularity and well-established methods in the field. Temporal-based methods can be further subdivided into four categories: (a) temporal interpolation-replacement methods [19], (b) temporal filter methods [20], (c) temporal function-fitting methods [21], and (d) temporal deep learning models [22]. Temporal interpolation-replacement methods are the most commonly used reconstruction approaches among temporal-based methods. Carreiras et al. [19] proposed the maximum value compositing (MVC) method to mitigate the effect of clouds and cloud shadows on optical images, where the actual NDVI would be smaller than the real NDVI due to the presence of clouds and cloud shadows [23]. However, MVC did not fully exploit the temporal information of optical remote sensing images. Among the temporal filter methods, the multi-temporal remote sensing images are decomposed into a one-dimensional time signal, and then the NDVI time series data of each pixel is reconstructed with different filters [20]. For instance, Chen et al. [24] first introduced the Savitzky–Golay (SG) filter into NDVI time series reconstruction, based on least squares fitting to eliminate the effect of noise in the target window [25]. In contrast, the SG filter aims to approximates the envelope of the NDVI time series, rather than simulating the phenological curve of vegetation [14]. Consequently, the SG filter exhibits poor performance in scenarios with excessive cloud contamination or during periods of vegetation growth [26]. Furthermore, a priori knowledge is needed to determine the sliding window size and polynomial order before using the SG filter, which may limit its performance. The Whittaker smoothing (WS) filter is a special case of B-spline smoothing [27], which employs a simple algorithm to fit discrete time series datasets. The WS filter is known for its computational efficiency and ability to balance the fidelity of the raw data with the smoothness of the fitted curve [28], making it a popular choice in numerous studies. The temporal function-fitting method employs a specific function to simulate the vegetation growth patterns for generating smooth time series curves. For instance, the asymmetric Gaussian (AG) function-fitting algorithm is a nonlinear least squares fitting method based on an asymmetric Gaussian function [21]. This method can accurately represent seasonal vegetation growth curves and estimate phenological parameters [19]. However, it is essential to identify a series of maximum and minimum values. When the original NDVI time series contains substantial noise, or the growing season is too short, extracting these two parameters becomes difficult [29], potentially leading to a loss of detail concerning vegetation changes. Deep learning has garnered significant attention in the remote sensing field due to its powerful feature-extraction capability [30]. Zhao et al. [31] have proposed the genetic algorithm–artificial neural network (GA-ANN) algorithm to reconstruct NDVI time series. Nevertheless, the need for a large number of training samples and high computational resources in training the model are disadvantages of temporal deep learning models. Therefore, the generalization ability of deep learning is limited [22]. In general, time-based methods merely utilize the temporal information of remote sensing data without fully leveraging the spatial information [14]. While time-based methods may perform well in reconstructing a small number of missing pixels, their effectiveness in reconstructing long-term gaps is limited.

Numerous novel spatiotemporal fusion algorithms have been developed to take full advantage of the temporal and spatial information in remote sensing images [32,33,34]. To reconstruct NDVI time series data with high temporal and spatial resolution using MODIS and Landsat images [35,36], fine pixels are predicted by calculating the neighborhood information of pixels within the target window. For instance, Gao et al. [37] proposed the spatial and temporal adaptive reflectance fusion model (STARFM), which fuses Landsat and MODIS images by leveraging information from similar neighboring pixels for refined increment estimation [38]. Subsequently, Zhu et al. [39] introduced an improved version, ESTARFM, which incorporates additional data pairs on the benchmark date and adopts a linear hybrid model to enhance the fusion performance in a heterogeneous region [40]. Rao et al. [41] presented a demixing model known as the linear mixed growth model (LMGM), which decomposes coarse increments into fine increments. Adapting the LMGM for heterogeneous landscapes, Zhu et al. [42] proposed the flexible spatiotemporal data fusion (FSDAF) technique, based on unmixing spectral analysis and a thin plate spline interpolator, and utilizing MODIS and Landsat data. Subsequently, Liu et al. [43] introduced an enhanced version of IFSDAF, followed by the development of FSDAF 2.0 by Guo et al. [44]. Nonlinear models have also been employed for reconstructing NDVI time series. Song et al. [45] developed the sparse representation-based spatiotemporal reflectance fusion model (SPSTFM). In contrast to traditional methods, Liu et al. [46] proposed the extreme-learning-machine-based spatiotemporal fusion model (ELM-FM), which focuses on learning a mapping function applied directly to different images rather than utilizing a complex feature representation followed by a simple mapping. Liu et al. [47] introduced the spatiotemporal fusion incorporating spectral autocorrelation (FIRST) model to optimally exploit the valuable spectral autocorrelation present across multiple bands.

However, reconstructing long-term and large-scale NDVI time series data using these algorithms is challenging. Several of the mentioned algorithms rely on MODIS and Landsat image pairs to predict accurate NDVI images. Meanwhile, these algorithms demand substantial computational resources and storage space. To address these issues, Google Earth Engine (GEE) emerged as an online processing platform, providing a convenient technical tool for remote sensing data processing [48]. Chen et al. [49] have developed a practical method to reconstruct high-quality Landsat–MODIS NDVI time series data employing a gap-filling and Savitzky–Golay filter (GF-SG), using GEE. The GF-SG method reconstructed NDVI time series with a spatial resolution of 30 m and a temporal resolution of 8 days by fusing MODIS and Landsat images. Compared with the three previous methods, the results of the study showed that GF-SG performed better than IFSDAF [43], STAIR [50], and fill and fit [51]. Hu et al. [52] conducted a comparative analysis of four spatiotemporal fusion models, namely STARFM, ESTARFM, FSDAF, and GF-SG, in reconstructing NDVI time series. The results show that the GF-SG method outperforms the other models in terms of accuracy, as it effectively generates NDVI images and constructs time series data with high spatial and temporal resolution.

Unfortunately, the NDVI with a spatial resolution of 30 m no longer meets the requirements for precision agriculture applications [53]. The Sentinel-2 satellite provides a wealth of spectral information and a high spatial resolution of 10 m to address this imbalance [54]. Our goal is to utilize image data from diverse satellite sources to reconstruct high-resolution NDVI time series data, yielding invaluable information for crop mapping and phenology extracting. In this study, we propose an enhanced gap-filling and Whittaker smoothing (EGF-WS) method to reconstruct NDVI time series using GEE. High spatial and temporal resolution NDVI time series were recovered using a gap-filling and Whittaker smoothing filter by employing three image collections with different spatial and temporal resolutions, namely MODIS, Landsat-8, and Sentinel-2. With reference to the all-round performance assessment (APA) metrics [55], the performance of the EGF-WS method was compared with the GF-SG approach in Fusui County, Guangxi Zhuang Autonomous Region.

The rest of this paper is organized as follows: Section 2 describes in detail the experimental area and data sources. Section 3 describes the EGF-WS method in detail and its accuracy evaluation method. Quantitative and qualitative perspectives on the accuracy of the EGF-WS reconstruction of NDVI are analyzed in Section 4. The selection of each parameter and the advantages and disadvantages of the EGF-WS are discussed in Section 5. Last, the conclusion is given in Section 6.

2. Study Area and Data

2.1. Study Area

The study domain, as illustrated in Figure 1, Fusui County is located in the southwest of Guangxi Zhuang Autonomous Region. The geographical coordinates are 107°3′–108°6′ E, 22°11′–22°57′ N, covering an area of approximately 2841 km². It has a subtropical monsoon climate, with 130–200 days of rainfall annually. The frequent cloudy and rainy weather in Fusui County poses challenges for obtaining cloud-free optical images. As a result, it becomes necessary to reconstruct high-quality cloud-free normalized difference vegetation index (NDVI) images utilizing multi-source remote sensing data fusion techniques for remote sensing application studies. To assess the performance of the reconstructed NDVI time series, two regions at distinct locations were chosen in Fusui, shown in Figure 1c,d. The land cover data was downloaded from https://viewer.esa-worldcover.org/worldcover (accessed on 6 September 2022) at a spatial resolution of 10 m. The dominant land cover types in the region include open water, trees, built-up areas, cropland, and barren land.

2.2. Data

2.2.1. MODIS Image Collection

To obtain NDVI time series data with high temporal resolution, we collected surface reflectance images from MODIS [56] (“MODIS/061/MOD09A1”) provided by NASA for the entire year of 2021. Each MOD09A1 image features seven spectral bands, all of which have undergone atmospheric correction. These images have a temporal resolution of 8 days and a spatial resolution of 500 m. Based on image coverage, clouds, and aerosols, a value was selected for synthesis from all images synthesized over an 8-day period. The B2 and B1 bands of the MOD09A1 were used to calculate the MODIS NDVI. To generate the cloud-free MODIS NDVI time series data, the “StateQA” band of each image was employed to identify clouds, cloud shadows, and cirrus masks. We note that there exist different MODIS image collections in GEE, such as the MOD09Q1, which have the same temporal resolution (8 days) and a high spatial resolution (250 m). Unfortunately, the MOD09Q1 images have only two spectral bands, but the MOD09A1 images have seven spectral bands. Consequently, more vegetation indices, such as EVI and LSWI, can be calculated using the MOD09A1 images. We aim to utilize EGF-WS not only for reconstructing NDVI time series, but also for reconstructing EVI and LSWI time series. For this reason, the MOD09A1 was selected as the reference image collection for high temporal resolution NDVI in this study.

2.2.2. Landsat-8 Image Collection

For the medium spatial resolution images, we collected the Landsat-8 [57] surface reflectance image collection (“LANDSAT/LC08/C02/T1_L2”) via the GEE for the whole year 2021. These images, having undergone atmospheric corrections, encompass four visible bands, one near-infrared band, and two short-wave infrared bands. With a 16-day revisit cycle and a spatial resolution of 30 m, the Landsat-8 satellite offers reliable and consistent data. To calculate the Landsat NDVI, we used the B5 and B4 bands of the Landsat-8. Furthermore, we employed the “QA_PIXEL” band to mitigate the impact of clouds, cloud shadows, and cirrus clouds.

2.2.3. Sentinel-2 Image Collection

In addition, the Sentinel-2 [58] surface reflectance image collection (“COPERNICUS/S2_SR”) provided by ESA for the whole year of 2021 was adopted as the high spatial resolution image collection. The images in this collection have also been atmospherically corrected with a revisit period of 5 days and a spatial resolution of 10 m. However, the revisit period of the Sentinel-2 satellite is 10 days at low latitudes. The B8 and B4 bands of Sentinel-2 were used to calculate Sentinel NDVI, and the “QA60” band was employed to remove the effects of clouds, cloud shadows, and cirrus clouds.

The number of all available MODIS, Landsat-8, and Sentinel-2 images in Fusui County are shown in Figure 2. The numbers of MODIS, Landsat-8, and Sentinel-2 imagery for the whole year of 2021 are 46, 68, and 438, respectively. The number of cloud-free observations for each remote sensing image: 1 to 37 for MODIS, 2 to 31 for Landsat-8, and 15 to 115 for Sentinel-2. This variance demonstrates the poor quality of optical images, with the number of cloud-free images significantly diminishing in response to inclement weather conditions.

3. Methodology

In this study, we proposed the EGF-WS method to reconstruct high spatiotemporal resolution NDVI time series. This process involves three main steps, that is, firstly, the MODIS image collection was preprocessed to generate cloud-free MODIS NDVI as a reference NDVI time series, providing high temporal resolution NDVI. The Sentinel–Landsat NDVI (SL-NDVI) time series at 16-day intervals were generated using the maximum synthesis method. Secondly, a new NDVI time series, Sentinel–Landsat–MODIS NDVI (SLM-NDVI), was generated using the enhance gap-filling method to fill in the missing values of the SL-NDVI images, thereby providing NDVI with high spatial resolution. Finally, we smoothed the SLM-NDVI time- series data by employing the Whittaker smoothing (WS) method to generate Sentinel–Landsat–MODIS Whittaker smoothing NDVI (SLM-WS-NDVI) time series with both high temporal and high spatial resolution. The schematic diagram illustrating the flow of EGF-WS processing is shown in Figure 3.

3.1. Preprocessing

In this study, we used a preprocessing step to eliminate the presence of cloud contamination in MODIS images. To fill the missing MODIS NDVI pixels, the missing values in the MODIS NDVI images were calculated using a linear interpolation method, as shown in Equation (1). We used all available MODIS NDVI pixels within 60 days to calculate the missing pixel information.

$$NDV{I}_i=NDV{I}_{i-1}+\frac{NDV{I}_{i+1}-NDV{I}_{i-1}}{t_{i+1}-t_{i-1}}\left(t_i-t_{i-1}\right)$$

(1)

where $NDV{I}_i$ represents the NDVI of DOY $i$, and $NDV{I}_{i-1}$ and $NDV{I}_{i+1}$ represent the NDVI at DOY $i-1$ and $i+1$, respectively. $t_i$ represents the DOY $i$, and $t_{i-1}$ and $t_{i+1}$ represent the DOY $i-1$ and DOY $i+1$, respectively.

We replaced missing NDVI pixels with 0 when NDVI was unavailable within 30 days. In this study, the MODIS NDVI was sharpened from 500 m to 10 m using a bicubic interpolation method in order to ensure consistency of resolution between different satellite sensors. The interpolated time series were smoothed using Whittaker smoothing to reduce residuals in MODIS images caused by cloud contamination, and linear interpolation was employed to generate high-quality, cloud-free MODIS NDVI time series, as shown in Figure 4a.

In this study, Landsat-8 NDVI and Sentinel-2 NDVI were combined to generate SL-NDVI, Figure 4b. First, the preprocessed Landsat-8 image was sharpened to 10 m using the bicubic interpolation method. Then, Landsat-8 NDVI and Sentinel-2 NDVI data over 16 days were synthesized using the maximum synthesis method (MVC), as indicated in Equation (2). The MVC was used because the acquisition time between images is different, and the MVC can effectively eliminate low-value noise.

$$\mathrm{SL}-\mathrm{NDVI}=Arcmax(Sentinel\_NDV{I}_{i+16}^i, Landsat\_NDV{I}_{i+16}^i)$$

(2)

where $Sentinel\_NDV{I}_{i+16}^i$ and $Landsat\_NDV{I}_{i+16}^i$ are all NDVI within 16 days for Sentinel-2 images and Landsat-8 images, respectively.

3.2. Enhance Gap-Filling

The MODIS NDVI time series data were used as reference data to fill missing NDVI pixels in SL-NDVI according to the change of time series curve shape. That is, according to the first theorem of geography, everything is related to everything else, but near things are more related to each other. In this study, we searched for neighboring MODIS NDVI pixels (named $M\_serie{s}_{\left(xi,yi\right)}$) within a window of 200 m × 200 m of each SL-NDVI pixel (named $SL\_trage{t}_{\left(xi,yi\right)}$). Then, the linear correlation coefficients ($w_{\left(xi,yi\right)}$) between $M\_serie{s}_{\left(xi,yi\right)}$ and $SL\_trage{t}_{\left(xi,yi\right)}$ were calculated to estimate the reference MODIS NDVI time series ($M\_referenc{e}_{\left(x,y\right)}$). as shown in Figure 4b, $M\_referenc{e}_{\left(x,y\right)}$ was calculated as follows:

$$M\_referenc{e}_{\left(x,y\right)}={{\displaystyle \sum }}_{i=1}^n{w}_{\left(xi,yi\right)}\times M\_serie{s}_{\left(xi,yi\right)}$$

(3)

where $n$ is the number of pixels in the search window and different $M\_serie{s}_{\left(xi,yi\right)}$ within the window have different weights $w_{\left(xi,yi\right)}$, this can be expressed as:

$$w_{\left(xi,yi\right)}=\frac{R_{i·\left(x,y\right)}}{{{\displaystyle \sum }}_{i=1}^n{R}_{i·\left(x,y\right)}}, where R_{i·\left(x,y\right)}=\frac{co{r}_{i·\left(x,y\right)}-co{r}_{min·\left(x,y\right)}}{co{r}_{max·\left(x,y\right)}-co{r}_{min·\left(x,y\right)}}$$

(4)

where $co{r}_{i·\left(x,y\right)} \left(=1,2,···,n\right)$ is the correlation coefficient between the $i_{th}$ pixel and $SL\_trage{t}_{\left(xi,yi\right)}$ in the search window; and $co{r}_{min·\left(x,y\right)}$ and $co{r}_{max·\left(x,y\right)}$ are the minimum and maximum correlation coefficients of all pixels and $SL\_trage{t}_{\left(xi,yi\right)}$ in the search window, respectively. Indeed, for each pixel within the window, the higher the similarity, the higher the correlation coefficient, and the higher its contribution to the generation of the reference time series.

Since the target window was only 200 m × 200 m, it may be difficult to generate $M\_referenc{e}_{\left(x,y\right)}$ for a few $SL\_trage{t}_{\left(xi,yi\right)}$. For these pixels, the $M\_referenc{e}_{\left(x,y\right)}$ was obtained by calculating the time series average of neighboring pixels within a 200 m × 200 m window. Then, the SLM-NDVI time series was obtained by modifying the $M\_referenc{e}_{\left(x,y\right)}$ based on least squares linear transfer, as shown in Figure 4c. SLM-NDVI is calculated as follows:

$$SLM-NDVI=M\_referenc{e}_{\left(x,y\right)}\times a\left(x,y\right)+a_0\left(x,y\right)$$

(5)

where $SLM-NDVI$ represents the NDVI time series fused with Sentinel-2, Landsat 8, and MODIS imagery; and $a\left(x,y\right)$ and $a_0\left(x,y\right)$ are the slope and intercept in the linear transfer function, respectively.

3.3. Whittaker Smoothing

We combined SLM-NDVI and SL-NDVI to generate the integrated NDVI (Figure 4c). The integrated NDVI time series may be unsmooth due to cloud contamination and residual noise after linear interpolation. To cope with this, the Whittaker smoothing algorithm was employed in this study. Then, the smooth SLM-WS-NDVI time series, namely, EGF-NDVI, with a 10 m spatial resolution and an 8-day temporal resolution, was constructed, as shown in Figure 4d. The Whittaker smoothing is a simple and convenient smoothing that performs optimally with evenly spaced data. For $y$ (integrated NDVI time series) with noise, a smooth $z$ (EGF-NDVI) was expected to be generated by Whittaker smoothing. Compared with SG filtering, WS filtering requires only one parameter ($\mathsf{\kappa }$) to smooth the NDVI time series. In this experiment, $\mathsf{\kappa }=1$ was chosen to smooth the integrated NDVI time series. The fidelity of $\mathrm{y}$ and the roughness of $\mathrm{z}$ need to be balanced during the smoothing process, and the combination $Q$ can be used to represent the fitting effect. $Q$ is calculated as follows:

$$Q=S+\kappa R$$

(6)

$$S={\displaystyle \sum }_i{\left(y_i-z_i\right)}^2$$

(7)

$$R={\displaystyle \sum }_i{\left(\left(z_i-z_{i-1}\right)-\left(z_{i-1}-z_{i-2}\right)\right)}^2$$

(8)

where fidelity $S$ is the sum of squares of the difference, $R$ is the roughness of the smoothed data, $\kappa$ controls the smoothness; the larger the $\kappa$ is, the worse the fidelity of the $y$ data and the smoothing $z$ are, and vice versa.

3.4. Accuracy Metrics

To evaluate the performance of EGF-WS in this study, the all-round performance assessment (APA) proposed by Zhu et al. [55] was adopted to evaluate the accuracy of EGF-NDVI. Three metrics, RMSE, AD, and edge were selected to quantitatively evaluate the accuracy of EGF-NDVI. RMSE represents the error of EGF-NDVI, with values ranging from 0 to 1, with smaller RMSE values indicating lower errors in EGF-NDVI, and higher values indicating higher errors. AD represents the average deviation between the reconstructed NDVI and the reference NDVI, with values ranging from −1 to 1. An AD value closer to 0 signifies a smaller deviation in the reconstructed NDVI, while a value between −1 and 0 indicates that EGF-NDVI is underestimated, and a value between 0 and 1 implies overestimation. Edge describes the edge information of the EGF-NDVI, with values ranging from −1 to 1. An edge value closer to 0 indicates a perfect match of edge features in EGF-NDVI, while a value between −1 and 0 suggests over-smoothing of edge features, and a value between 0 and 1 indicates over-sharpening. Three accuracy metrics are calculated as follows:

$$RMSE=\sqrt{\frac{{{\displaystyle \sum }}_{i=1}^N{\left(F_i-R_i\right)}^2}{N}}$$

(9)

$$AD=\frac{1}{N}{{\displaystyle \sum }}_{i=1}^N\left(F_i-R_i\right)$$

(10)

$$Edge=\left|D_{i,j}-D_{i+1,j+1}\right|+\left|D_{i,j+1}-D_{i+1,j}\right|$$

(11)

where $F_i$ is the NDVI of a pixel $i_{th}$ in the EGF-NDVI, $R_i$ is the NDVI of a pixel $i_{th}$ in the reference image, $N$ is the total number of pixels, and $D_{i,j}$ is the NDVI of the $i_{th}$ row and $j_{th}$ column.

4. Results

Yang et al. have previously demonstrated that GF-SG outperforms three commonly used methods (IFSDAF, STAIR, and fill and fit) in reconstructing Landsat NDVI time series. In this study, to verify the accuracy of our proposed method, we selected GF-SG as a benchmark and compared it with our approach. Specifically, NDVI time series were reconstructed using EGF-WS and GF-SG methods, referred to as EGF-NDVI and GF-NDVI, respectively. These reconstructed time series were then compared with Sentinel-2 NDVI in Region A and Region B to assess the performance of both methods. Subsequently, we conducted a comprehensive analysis of the performance of EGF-WS and GF-SG, considering both quantitative and qualitative perspectives.

4.1. Qualitative Assessment

To verify the feasibility of the reconstruction algorithm and the reasonableness of the result accuracy, The EGF-NDVI was verified in this study using Sentinel-2 NDVI. Consequently, Sentinel-2 NDVI served as the true NDVI, and EGF-NDVI and GF-NDVI were used for qualitative evaluation and comparison. The EGF-NDVI and GF-NDVI images were dated 6 December 2021, while the Sentinel-2 NDVI image was dated 10 December 2021.

We evaluated the consistency of the two methods in terms of spatial patterns. Based on a visual comparison, as shown in Figure 5, the spatial distribution of EGF-NDVI and GF-NDVI are consistent with those of Sentinel-2 NDVI in Region A. From Figure 5, it is evident that EGF-NDVI captures the spatial details of the farmland; however, the farmland boundaries appear blurred in GF-NDVI, with the area marked in red on Figure 5d–f. The road within the farmland is clearly visible in both Sentinel-2 NDVI and EGF-NDVI; however, it is not well-represented in GF-NDVI. This discrepancy stems from differences in spatial resolution between EGF-NDVI and GF-NDVI, which result in variations in spatial detail within the NDVI outcomes. The spatial resolution of EFG-NDVI is 10 m, consistent with Sentinel-2 NDVI, and much higher than the 30 m resolution of GF-NDVI. In Region B, EGF-NDVI displays open water boundaries similar to those of Sentinel-2 NDVI, while GF-NDVI presents unclear boundaries, present in Figure 5g–i. Regarding buildings near the open water, building boundaries in EGF-NDVI more closely resemble those in Sentinel-2 NDVI, while GF-NDVI displays blurred boundaries, marked in red in Figure 5g–i. EGF-NDVI captures building boundary lines and corner points, whereas GF-NDVI exhibits curved boundary lines and inaccurate corner points.

EGF-NDVI offers more precise spatial detail than GF-NDVI due to its 10 m spatial resolution, compared to GF-NDVI’s 30 m resolution. This allows EGF-NDVI to capture the boundaries of farmland and buildings more accurately. These findings demonstrate that the NDVI after EGF-WS reconstruction clearly depicts farmland, open-water-body, and building boundaries, which is beneficial for extracting these features using NDVI images.

To assess the reasonableness of the NDVI time series accuracy after EGF-WS reconstruction, considering its practical application, two distinct land cover types were selected in Region A and Region B: evergreen vegetation in Figure 6a, and agricultural land in Figure 6b. Subsequently, time series curves were plotted for these land cover types, and the RMSE of EGF-NDVI and GF-NDVI was calculated in comparison with Sentinel-2 NDVI to assess their accuracies.

The NDVI time series curves of evergreen vegetation using EGF-WS and GF-SG are shown in Figure 6a. As we all know, strong reflection in the near-infrared band and absorption characteristics in the red band yield NDVI values above 0.5. As shown in Figure 6a, the NDVI of evergreen vegetation exceeds 0.5. In the study area, evergreen vegetation remains green year-round, and both EGF-NDVI and GF-NDVI were consistent with the phenological characteristics of evergreen vegetation. During the DOY 200–300 period, the EGF-NDVI generally outperforms GF-NDVI, accurately reflecting high NDVI around DOY 170 and maintaining the retrieved NDVI within a reasonable range. Figure 6b demonstrates the time series curve of sugarcane within agricultural land. Sugarcane, an annual crop in Fusui, has the highest subsurface soil exposure between January and March. As shown in Figure 6b, GF-NDVI overestimates NDVI during the DOY 30–60 period before crop emergence, whereas EGF-NDVI captures this decrease more effectively. Moreover, the EGF-NDVI time series curve closely approximates Sentinel-2 NDVI during the sugarcane maturation period (between November and December), indicating that EGF-WS can effectively reconstruct NDVI time series data and there is superior accuracy with EGF-NDVI.

To illustrate the differences between predicted and actual NDVI in this study, we present the RMSE values obtained from the reconstructed NDVI using two different methods, EGF-NDVI and GF-NDVI, in both evergreen vegetation and agricultural land, as shown in Figure 6. For evergreen vegetation (Figure 6a), the RMSE values of EGF-NDVI and GF-NDVI are identical at 0.11. This result implies that both methods exhibit equally well in estimating the NDVI for evergreen vegetation, demonstrating their suitability and precision for this type of land cover. Conversely, for agricultural land (Figure 6b), the RMSE values display a significant discrepancy between the two methods. The EGF-NDVI produces a significantly lower RMSE value of 0.04, while the GF-NDVI results in a higher RMSE value of 0.10. This observation indicates that the EGF-NDVI surpasses the GF-NDVI in estimating NDVI for agricultural land, particularly sugarcane fields. The reduced RMSE value achieved with the EGF-NDVI indicates that it offers more accurate and dependable NDVI estimations for agricultural land, consequently improving our capacity to monitor crop growth and identify crop types. This finding supports the potential value of the EGF-NDVI for agricultural remote sensing applications.

4.2. Quantitative Assessment

To compare the accuracy of EGF-NDVI and GF-NDVI, we used three metrics, RMSE, AD, and edge, for quantitative assessment, in accordance with the methods described in Zhu et al. The accuracy metrics of EGF-NDVI and GF-NDVI were calculated separately, using Sentinel-2 NDVI as the reference image, and the results are shown in

Table 1. All-round performance assessment (APA) metrics (AD, RMSE, EDGE, and LBP) of EGF-NDVI and GF-NDVI, in regions A and B. The best results are marked in bold.

Test Area	Reconstructed NDVI	AD	RMSE	EDGE
Region A	EGF-NDVI	3.29 × 10−6	1.03 × 10−5	−0.344
	GF-NDVI	5.71 × 10−6	1.34 × 10−5	−0.561
Region B	EGF-NDVI	0.72 × 10−6	0.67 × 10−5	−0.010
	GF-NDVI	3.79 × 10−6	1.445 × 10−5	−0.5817

. The EGF-NDVI and GF-NDVI images were dated 6 December 2021, while the Sentinel-2 NDVI image was dated 10 December 2021.

The APA performance metrics for EGF-NDVI and GF-NDVI predictions are presented in Table 1. It can be seen that both EGF-NDVI and GF-NDVI overestimate Sentinel-2 NDVI to different degrees in Region A and Region B. However, the EGF-WS outperforms GF-SG, and EGF-NDVI has lower AD and RMSE values than GF-NDVI. Moreover, the AD and RMSE of EGF-NDVI are only 0.72 × 10⁻⁶ and 0.67 × 10⁻⁵, respectively, while GF-NDVI exhibits more severe overestimation in Region B (AD = 3.79 × 10⁻⁶, RMSE = 1.445 × 10⁻⁵), which means that the deviation between EGF-NDVI and Sentinel-2 NDVI is smaller and the accuracy of EGF-NDVI is higher. Both methods exhibit a certain smoothing phenomenon in terms of spatial accuracy. The smoothing effect of EGF-NDVI (EDGE = −0.344) is smaller than that of GF-NDVI (EDGE = −0.561) in Region A, and EGF-NDVI shows better fidelity (EDGE = −0.01) in Region B. This indicates that the spatial features of EGF-NDVI are more compatible with those of Sentinel-2 NDVI, and the texture information is more accurate. This may be attributed to the higher spatial resolution (10 m) of the reference image compared to GF-NDVI (30 m), further demonstrating that EGF-NDVI has a higher spatial resolution and richer texture information than GF-NDVI.

To further validate the authenticity of EGF-NDVI, in this paper, the scatter plots of EGF-NDVI and GF-NDVI images with reference image were calculated, as shown in Figure 7.

As can be seen in Figure 7, both EGF-NDVI and GF-NDVI show a reasonable agreement with Sentinel-2 NDVI, with generally improved performance from EGF-NDVI, as indicated by R² values ranging from 0.915 to 0.926. Similar results were observed in Region B. This indicates that the performance of EGF-WS surpasses that of GF-SG, and EGF-NDVI closely resembles Sentinel-2 NDVI, exhibiting a strong correlation in terms of the experimental results.

5. Discussion

Optical imagery is subject to spatial and temporal discontinuities, as remote sensing sensors are influenced by clouds and precipitation. Under such circumstances, Chen et al. [49] have proposed a GF-SG algorithm to reconstruct the NDVI time series. However, the spatial resolution of GF-NDVI is limited to 30 m, which is insufficient for the evolving demands of precision agriculture. In this study, we proposed an EGF-WS algorithm to reconstruct NDVI time series with a spatial resolution of 10 m and a temporal resolution of 8 days by blending MODIS, Landsat-8, and Sentinel-2 images. To determine the optimal parameters for the EGF-WS, we explored the effects of different window sizes and Whittaker smoothing parameter settings on the reconstructed NDVI time series. We also calculated a series of performance indicators and statistical metrics to evaluate the quality of EGF-NDVI. Finally, we analyzed the benefits and challenges associated with the EGF-WS approach.

5.1. Parameter Determination of the Window Size

To evaluate the accuracy of NDVI after reconstruction using different search windows and smoothing windows, the RMSE of EGF-NDVI and Sentinel-2 NDVI were calculated in step 5, as shown in Figure 8.

Figure 8 illustrates the accuracy of the reconstructed NDVI with different window sizes. The larger the search window size, the more pixels are available, which provides more reference data for calculating the NDVI time series of the target pixel but also requires more computational resources. In this paper, we examined the impact of different search window sizes on the accuracy of the reconstructed NDVI. A smaller smoothing window size results in higher fidelity of the smoothed image and less noise in the EGF-NDVI time series, and vice versa. As can be seen in Figure 8, the RMSE of the reconstructed NDVI image decreases as the size of the search window increases. Meanwhile, the size of the smoothing window also influences the accuracy of EGF-NDVI. A smaller smoothing window size leads to a reduced RMSE of the EGF-NDVI image. When the smoothing window is 5, the accuracy of EGF-NDVI is higher. Considering the balance of computation time and reconstruction accuracy, we chose the optimal parameters with a search window size of 20 and a smoothing window size of 5.

5.2. Parameter Determination of the Whittaker Smoothing

The Whitaker smoothing filter minimizes noise errors and generates smooth NDVI time series curves without distortion. The WS filter was chosen to smooth the EGF-NDVI time series in this paper because it is simple, efficient, and easy to implement. It requires only one parameter, $\mathsf{\kappa }$, to smooth and eliminate the noise effects in the original time series data. Larger $\mathsf{\kappa }$ may distort NDVI time series curves, while smaller $\mathsf{\kappa }$ cannot adequately suppress noise. For this reason, we experimented with the effect of different $\mathsf{\kappa }$ on the EGF-NDVI time series to determine the optimal parameters, the results are shown in Figure 9.

We calculated the EGF-NDVI time series curves based on different $\mathsf{\kappa }$ within 0–10. From Figure 9a, differences in $\mathsf{\kappa }$ can decrease the residual noise in the NDVI time series curves. The NDVI of DOY 233 is below normal due to cloud-detection errors. We observed that the impact of noise on NDVI diminishes as $\mathsf{\kappa }$ increases. When $\mathsf{\kappa }$ is set to 5, the NDVI time series curve can be effectively smoothed to remove noise and without overfitting the NDVI time series. Therefore, $\mathsf{\kappa }=5$ was chosen as the best parameter for the Whittaker smoothing filter in this study.

5.3. Performance Metrics and Statistical Metrics

The performance metrics part is very important for the evaluation of application results. For this aim, we compare and analyze the performance of two reconstructed NDVI methods: EGF-NDVI and GF-NDVI in two different geographic regions (Region A and Region B). Meanwhile, we have calculated additional performance evaluation metrics. These include Nash–Sutcliffe efficiency (NSE), where NSE values range from -inf to 1. Higher NSE values, closer to 1, indicate better model performance, while values below 0 indicate poor model performance. In root squared relative error (RSR), the closer the RSR value is to 0, the better the model performance, while a value close to 1 indicates poor model performance. In mean absolute error (MAE), the smaller the value of MAE, the smaller the model prediction error and the better the model performance. In coefficient of determination (R²), R² ranges from 0 to 1, where a higher R² value indicates a stronger explanatory power of the model and a better fit to the observed data. On the other hand, a lower R² value suggests a weaker ability of the model to explain the data and a poorer fit. The results are shown in

Table 2. Performance metrics of EGF-NDVI and GF-NDVI in Regions A and B.

Test Area	Reconstructed NDVI	NSE	RSR	MAE	R2
Region A	EGF-NDVI	0.7277	0.4046	0.0709	0.926
	GF-NDVI	0.5246	0.5285	0.0968	0.897
Region B	EGF-NDVI	0.5805	0.4564	0.0781	0.915
	GF-NDVI	0.233	0.5644	0.1049	0.836

.

First, by observing the data in Table 2, we find that the NSE values of EGF-NDVI in Region A and Region B are 0.7277 and 0.5805, respectively, while the NSE values of GF-NDVI are 0.5246 and 0.233, respectively. This indicates that the EGF-NDVI method has a better fit in both regions than the GF-NDVI method. Particularly in Region B, the NSE value of EGF-NDVI is much higher than that of GF-NDVI, indicating that the predictive performance of the EGF-NDVI method is more outstanding in this region. Second, we observe the RSR data. The RSR values of EGF-NDVI in Region A and Region B are 0.4046 and 0.4564, respectively, while those of GF-NDVI in the two regions are 0.5285 and 0.5644, respectively. A lower RSR value implies smaller errors, so in both regions, the EGF-NDVI method has smaller errors than the GF-NDVI method. Next, we focus on the MAE data. The MAE values of EGF-NDVI in Region A and Region B are 0.0709 and 0.0781, respectively, while those of GF-NDVI are 0.0968 and 0.1049, respectively. This again confirms that the EGF-NDVI method has smaller errors than the GF-NDVI method in both regions. Finally, we compare the R² values of the two methods in both regions. The R² values of EGF-NDVI in Region A and Region B are 0.926 and 0.915, respectively, while those of GF-NDVI are 0.897 and 0.836, respectively. A higher R² value indicates a better fit of the model to the data. Therefore, in both regions, the EGF-NDVI method has a better fit than the GF-NDVI method.

In summary, through a detailed discussion of the table data, we can draw the following conclusions: in both study regions, the EGF-NDVI method has better fitting performance and predictive performance than the GF-NDVI method.

Statistical properties can provide insights into the central tendency, dispersion, and distribution of data. We computed statistical metrics, including measures such as min, max, mean, median, and coefficient of variation (cv), for Sentinel-2 NDVI, EGF-NDVI, and GF-NDVI to characterize the statistical properties of the image. The results are shown in

Table 3. Statistical metrics of Sentinel NDVI, EGF-NDVI, and GF-NDVI in Regions A and B.

Test Area	Reconstructed NDVI	Min	Max	Mean	Median	CV
Region A	Sentinel NDVI	−0.5901	0.9098	0.5240	0.6055	0.4842
	EGF-NDVI	−0.1466	0.9919	0.5571	0.6133	0.3481
	GF-NDVI	−4.0850	1.7273	0.5813	0.6413	0.3197
Region B	Sentinel NDVI	−0.4833	0.9175	0.5352	0.6290	0.4810
	EGF-NDVI	−0.0619	0.998	0.5806	0.6289	0.3025
	GF-NDVI	−0.1918	1.2546	0.5733	0.6174	0.2817

.

Table 3 presents the reconstructed NDVI values and additional statistics for different test areas, namely, Region A and Region B. In Region A, the Sentinel NDVI range from -0.5901 to 0.9098, with a mean of 0.5240 and a median of 0.6055. The CV is relatively high at 0.4842, indicating a moderate level of variability. For EGF-NDVI, the reconstructed NDVI values show a narrower range from −0.1466 to 0.9919, with a higher mean of 0.5571 and median of 0.6133 compared to Sentinel NDVI. The CV is lower at 0.3481, indicating reduced variability. On the other hand, the GF-NDVI yields a wider range of reconstructed NDVI values from −4.0850 to 1.7273, with a mean of 0.5813 and median of 0.6413. However, the CV is the lowest at 0.3197, indicating relatively low variability despite the wider range. In Region B, similar trends are observed with Sentinel NDVI, EGF-NDVI, and GF-NDVI. The EGF-NDVI shows the highest mean (0.5806) and median (0.6289) compared to the other two NDVI, with a relatively low CV of 0.3025. The GF-NDVI method shows the lowest CV at 0.2817, with a mean of 0.5733 and median of 0.6174.

In summary, in both study regions, the statistical metrics of EGF-NDVI, in comparison to GF-NDVI, demonstrate a closer resemblance to the true NDVI (Sentinel NDVI), suggesting that the NDVI reconstructed using EGF-WS exhibits distribution patterns and trend characteristics that are more similar to Sentinel NDVI.

5.4. Advantages and Limitations of the EGF-WS Method

Compared with the GF-SG algorithm, the NDVI reconstructed by EGF-WS demonstrates significant improvement in spatial resolution while maintaining the same temporal resolution. In addition, the EGF-WS method facilitates the convenient reconstruction of NDVI time series, requiring only the study area’s vector to generate high spatial and temporal resolution NDVI time series data, thus eliminating the need for cropping and mosaicking during data processing. Moreover, EGF-WS was implemented in GEE, which was not limited by local computer hardware and required only a browser to execute the code. However, EGF-WS is disadvantageous for long time series studies since GEE does not offer Sentinel 2 images before 2017.

6. Conclusions

In this study, we proposed a method called EGF-WS, which aids in reconstructing NDVI time series data with high spatiotemporal resolution, thus facilitating research in land use change monitoring, agriculture, and various ecosystems. The EGF-WS method leverages MODIS NDVI time series as a reference and integrates Landsat-8 and Sentinel-2 imagery through gap-filling and Whittaker smoothing techniques. We produced an NDVI time series for Fusui County, featuring an 8-day temporal resolution and a 10 m spatial resolution. A qualitative and quantitative comparison was conducted between EGF-NDVI and GF-NDVI. The results demonstrated that EGF-NDVI, with its superior 10 m spatial resolution compared to GF-NDVI’s 30 m resolution, accurately captures the intricacies of land features. EGF-NDVI precisely reflects NDVI variations across different land cover types, demonstrating a strong correlation to the reference NDVI (R² = 0.926, RMSE = 0.67 × 10⁻⁵). However, due to the optical image acquisition method, cloud contamination is inevitable. Sentinel-1 can penetrating clouds and is not affected by weather conditions. In future work, we will explore the utilization of Sentinel-1 imagery for NDVI reconstruction.