Temporal-Difference Graph-Based Optimization for High-Quality Reconstruction of MODIS NDVI Data

Abstract

The Normalized Difference Vegetation Index (NDVI) is a crucial remote-sensing metric for assessing land surface vegetation greenness, essential for various studies encompassing phenology, ecology, hydrology, etc. However, effective applications of NDVI data are hindered by data noise due to factors such as cloud contamination, posing challenges for accurate observation. In this study, we proposed a novel approach for employing a Temporal-Difference Graph (TDG) method to reconstruct low-quality pixels in NDVI data. Regarding spatio-temporal NDVI data as a time-varying graph signal, the developed method utilized an optimization algorithm to maximize the spatial smoothness of temporal differences while preserving the spatial NDVI pattern. This approach was further evaluated by reconstructing MODIS/Terra Vegetation Indices 16-Day L3 Global 250 m Grid (MOD13Q1) products over Northwest China. Through quantitative comparison with a previous state-of-the-art method, the Savitzky–Golay (SG) filter method, the obtained results demonstrated the superior performance of the TDG method, and highly accurate results were achieved in both the temporal and spatial domains irrespective of noise types (positively-biased, negatively-biased, or linearly-interpolated noise). In addition, the TDG-based optimization approach shows great robustness to noise intensity within spatio-temporal NDVI data, suggesting promising prospects for its application to similar datasets.

1. Introduction

The Normalized Difference Vegetation Index (NDVI) is a fundamental remote-sensing measure of land surface greenness that is widely used to assess ecosystem properties such as vegetation phenology, biomass production, soil moisture, and carbon sequestration [1,2,3]. NDVI also supports the further analyses like classification and regionalization [4], evaluating temporal dynamics and trends [5], and change detection [6]. However, the quality of NDVI products is often compromised by factors such as atmospheric conditions (e.g., clouds, dust, etc.), snow cover, varying sun-sensor-surface viewing geometries, and sensor faults. Noise from these factors must be carefully processed to avoid unreliable analyses [7]. Thus, employing data reconstruction (or noise reduction/restoration) techniques is crucial to enhance NDVI quality and ensure the reliability of its applications.

Over the past several decades, numerous techniques have been developed to reconstruct NDVI data. For instance, the resulting NDVI products must account for the bidirectional reflectance distribution function (BRDF) effect, and BRDF normalization is commonly used to mitigate noise introduced by observation and illumination angular effects [8]. Based on the use of temporal and spatial information, NDVI reconstruction methods can be generally categorized into (1) time-series-based methods, which reconstruct NDVI time series at each pixel individually without considering spatial associations between NDVI values at different locations, and (2) spatio-temporal methods, which simultaneously leverage temporal and spatial interrelations for reconstruction. Time-series-based methods have been well developed and widely applied in NDVI reconstruction, remaining prevalent in the field. These methods typically apply a mathematical filter to NDVI time series, either in the temporal domain or the frequency domain, to reduce noise. Specifically, methods in the temporal domain include the 4253H twice filter [9]; the best index slope extraction algorithm (BISE) [10,11]; asymmetric Gaussian (AG) function fitting [12]; the Savitzky–Golay (SG) filter [13]; the Mean-Value Iteration (MVI) filter [7]; double logistic function (DL) fitting [14]; the Whittaker smoother (WS) [15]; the changing-weight filter method [16]; and the RMMEH filter [17]. Methods in the frequency domain include Fourier analysis [18]; harmonic analysis (HA) [19]; and wavelet transformation [20]. Despite their widespread use, time-series-based methods have significant limitations due to their dependence on incomplete information. If an NDVI time series is dominated by data noise (imprecise information) or data gaps (missing information), particularly continuous noisy entries, reconstruction becomes challenging due to the lack of essential information. Unfortunately, low-quality NDVI time-series data are common in NDVI products. To enhance the effectiveness of time-series-based methods, researchers have recently explored additional information from the spatial domain of NDVI data for reconstruction. These are spatio-temporal NDVI-reconstruction methods. For example, Poggio et al. [21] proposed a hybrid Generalized Additive Model (GAM)–geostatistical space-time model (denoted as GAM-Geo) that combines a multidimensional GAM fitting the spatio-temporal trends of NDVI with a geostatistical approach interpolating the residual components of NDVI over space. Oliveira et al. [22] developed a window regression (WR) method, which replaces noisy NDVI values with new values regressed from the relationship between the NDVI time series at the target pixel and one selected from the nearest eight neighbors. Other methods, such as the correction based on spatial and temporal continuity (CSaTC) method [23] and the temporal-spatial iteration (TSI) method [24], iteratively restore some noisy NDVI values using data from high-quality pixels selected either temporally or spatially close to the noisy pixel.

In theory, spatio-temporal methods are expected to outperform time-series-based methods due to the additional information they utilize; however, these emerging methods are still immature, and each has its own drawbacks. For instance, the GAM-Geo method utilizes a multidimensional GAM model to incorporate spatio-temporal information, but fitting this model could be subjective due to the tradeoff between function smoothness and fitting accuracy. Moreover, neither the smooth function used in the GAM model nor the geostatistical approach (specifically, simple kriging interpolation) can precisely model non-continuous patterns of NDVI variability across space. Both the TSI method and the CSaTC method require supplementary data (such as land cover types or ecological zone maps) in addition to NDVI data. They also use spatial and temporal information separately, meaning that for each noisy pixel, the reconstructed NDVI values are determined using either spatial or temporal information alone. This approach can lead to an imbalance in the ratio between the pixel counts of the two aspects. Additionally, the WR method has been shown to perform even worse than time-series-based methods and is practically incapable of reconstructing all noisy pixels [24].

Despite previous efforts, two major obstacles hinder the effective use of spatial information in NDVI data reconstruction. First, distinct from the temporal dimension, NDVI’s spatial information lies in a two-dimensional space, adding complexity to its exploitation. To address this, some methods (e.g., WR, CsaTC, and TSI) extract spatial information from a very limited number of pixels, typically just one, selected within a small sliding window. Second, NDVI variability often exhibits unsmooth patterns over space, especially in high-resolution NDVI products. Abrupt transitions in land cover types or topographical features lead to uneven textures and irregular structures in NDVI spatial profiles. It is also worth mentioning that the introduction of spatial information into time-series reconstruction depends on the ground sampling distance and the related spatial resolution of the input data. Consequently, many common spatial data-processing techniques, such as smoothing, filtering, and interpolation, may fail in NDVI reconstruction. Due to the lack of effective solutions for such issues, methods like WR, CsaTC, and TSI must carefully select the most suitable pixel that is homogeneous to the noisy pixel based on land cover type, phenological pattern, or other criteria.

Thus, the objective of this study is to develop a novel NDVI reconstruction method that can take full advantage of the spatio-temporal interrelations in NDVI data. Inspired by recent advancements in graph signal processing, we applied a temporal-difference graph (TDG) method [25] to NDVI reconstruction. Originally designed to reconstruct partially-sampled time-varying graph signals, the TDG method is particularly well suited for NDVI reconstruction because it addresses the two major obstacles identified in previous efforts. First, rather than relying on local spatial windows, this method models spatial interrelations among data pixels using a graph structure, where every two pixels are directly or indirectly connected by edges. Second, the TDG method does not rely on the spatial smoothness of the data. Instead, it operates on the assumption that the temporal differences of the data (the first-order derivative of data versus time) should be smooth across space. This assumption is more realistic for NDVI data while NDVI values may not be smooth over space; the changes in NDVI over short periods tend to correlate strongly at adjacent locations. Based on these two principles, the TDG method effectively utilizes both temporal and spatial information in an asymmetric and rational manner.

2. Methodology

This section provides a clear description of the method used for reconstructing NDVI data, and Figure 1 illustrates a flowchart of the research methodology.

2.1. Graph Theory and Graph Signal Processing

Graph theory is a fundamental branch of mathematics that investigates structures consisting of vertices and edges, collectively known as graphs. A graph conceptualizes pairwise relations between objects; each vertex represents an object, and each edge represents the relationship between two linked objects. This generalized structure is versatile and has been used to model various types of relations and processes across a wide range of domains, including computer science, physics, biology, and engineering [26,27].

For clarity, we introduce some formal notation. An undirected and weighted graph is denoted as $G(V,E, \mathit{W})$, where $V$ denotes the vertex set $V=\{v_1,v_2, \dots ,v_{\mathrm{N}}\}$, $E$ denotes the edge set $E=\{e_1,e_2,\dots ,e_{\mathrm{M}}\}$, and $\mathit{W}$ denotes the adjacency matrix of the graph. Each individual vertex and edge are denoted as $v_i(1\le i\le \mathrm{N})$ and $e_j(1\le j\le \mathrm{M})$, respectively, with $\mathrm{N}$ and $\mathrm{M}$ representing the total number of vertices and edges, respectively. Each edge links two vertices, and we may alternatively denote an edge as $e^{(i_1,i_2)}$, indicating the edge that links vertices $v_{i_1}$ and $v_{i_2}$. These vertices are called the end vertices of edge $e^{(i_1,i_2)}$. The adjacency matrix $\mathit{W}$ is an $\mathrm{N}\times \mathrm{N}$ matrix defined as follows. If there exists an edge $e^{(i,j)}$ that connects two vertices $v_i$ and $v_j$, the corresponding entry $W_{ij}$ (the entry in the i^th and the j^th row) equals the weight associated with edge $e^{(i,j)}$; otherwise, $W_{ij}$ equals zero. The setting of edge weights is flexible and depends on the problem at hand. A common strategy involves assigning edge weights based on the similarity between the two vertices connected by the edge, and this approach is prevalent in studies focusing on large-scale graph modeling of real-world systems, such as social networks, computer networks, and biological networks. In addition to the adjacency matrix, another common matrix representation for a graph is the Laplacian matrix $\mathit{L}$, defined as $\mathit{L}=\mathit{D}-\mathit{W}$, where D is a diagonal matrix with the $i$^th diagonal element equal to the degree of vertex $v_i$, i.e., the sum of the weights of all the edges directly linked with vertex $v_i$. Note that $\mathit{L}$ is positive-semidefinite, which follows from its symmetry and diagonal dominance.

Graph theory has been developed for centuries, and even nowadays, with the expansion of its applications, new methods and techniques that merge graph theory with other domains are continuously emerging. Among these subdomains, an emerging field known as graph signal processing (GSP), or signal processing on graphs, has received great attention in recent years [28]. Graph signal processing combines graph theory and signal processing to deal with signal-processing problems such as filtering, compression, and reconstruction, specifically for data whose inherent associations or structures are represented by a graph. These data are collectively referred to as a graph signal, with each data sample attached to a vertex of the graph. Intuitively, the vertices in the graph act as signal sources, while the edges depict the relationships of these signal sources. We used the vector $\mathit{x}\in {\mathbb{R}}^{\mathrm{N}}$ to denote a graph signal built on a graph $G(V, E, \mathit{W})$. The $i$th component of $\mathit{x}$, denoted by $x_i$, is the data sample at the $i$^th vertex $v_i$ in the graph. If the data sample at each vertex is not a single value but a time series, this signal is referred to as a time-varying graph signal. In these cases, the signal is represented as the matrix $\mathit{X}\in {\mathbb{R}}^{\mathrm{N}\times \mathrm{T}}$, where $\mathrm{T}$ is the number of time instances. Each column of this matrix denotes the samples at all vertices in the graph at a single time instance, and each row represents the time-series samples at a single vertex throughout time. It is worth mentioning that we only consider graph signals with synchronous time instances for the time series at each vertex.

2.2. Time-Varying Graph Signal Reconstruction: The Temporal-Difference Graph Method

The realm that inspires our solution to the NDVI reconstruction problem is graph signal reconstruction, one of the core research topics in graph signal processing. A graph signal $\mathit{x}$ is considered partially sampled when some of its entries are not observed, prompting the need for methods to recover the missing entries based on the known entries. A fundamental approach to reconstructing graph signals involves assuming that the signal is smooth on the graph [29], meaning the missing entries are assigned to the values that maximize the overall smoothness of the signal. In most of the literature [25,28,30], the smoothness of a graph signal is typically measured as follows:

$$S\left(\mathit{x}\right)={\displaystyle \frac{1}{2}}{\mathit{x}}^{\mathrm{T}}\mathit{L}\mathit{x}$$

(1)

where $\mathit{L}$ is the Laplacian matrix of the related graph. The smaller $S\left(\mathit{x}\right)$ is, the better the smoothness that the graph signal exhibits. If we unfold the equation, we obtain a more intuitive form, as follows:

$$\begin{array}{c}S\left(\mathit{x}\right)={\displaystyle \frac{1}{2}}{\displaystyle \sum _{1\le i,j\le n}}{W}_{i,j}{\left(x_i-x_j\right)}^2\\ \hspace{1em}\hspace{1em}={\displaystyle \frac{1}{2}}{\displaystyle \sum _{e^{(i,j)}\in E}}{W}_{i,j}{\left(x_i-x_j\right)}^2\end{array}$$

(2)

which means the smoothness of a graph signal is measured as the weighted squared sum of the signal difference between every vertex pair linked by an edge. For a time-varying graph signal $\mathit{X}$, its smoothness can be calculated as the sum of the smoothness values at every time instance, as follows:

$$\left(\mathit{X}\right)={\displaystyle \frac{1}{2}}\sum _{1\le t\le \mathrm{T}}{\mathit{x}}_t^{\mathrm{T}}\mathit{L}{\mathit{x}}_t={\displaystyle \frac{1}{2}}\mathrm{t}\mathrm{r}\left({\mathit{X}}^{\mathrm{T}}\mathit{L}\mathit{X}\right)$$

(3)

where $\mathrm{t}\mathrm{r}(·)$ calculates the trace of a given matrix, and ${\mathit{x}}_t$ denotes the $t$th column of $\mathit{X}$, i.e., the graph signal at time instance $t$.

One way to reconstruct time-varying graph signals is to directly minimize $S\left(\mathit{X}\right)$. However, this method essentially reconstructs the graph signal at each time instance separately, disregarding the relationship between temporally-adjacent signals. Moreover, some signals in nature, including NDVI, exhibit non-smooth behavior over space and, therefore, reconstructing signals based solely on spatial smoothness may lead to significant deviations from reality. Due to the influence of geographic conditions on vegetation growth, NDVI signals would inherit spatial patterns from various factors such as topography, land cover, and soil and water variability, resulting in uneven textures and irregular structures in NDVI spatial profiles. To reconstruct such non-smooth time-varying graph signals, Qiu et al. [25] proposed a new method based on the fact that for many time-varying graph signals in nature, the temporal differences of the signal often exhibit better smoothness than the original signals. We refer to this method as the TDG method, and the follow-up Equations (4) and (5) can also be found in Qiu et al. [25].

The temporal differences of a time-varying graph signal $\mathit{X}$ are expressed as:

$$\begin{array}{c}\Delta \mathit{X}=\left[{\mathit{x}}_2-{\mathit{x}}_1,{\mathit{x}}_3-{\mathit{x}}_2,\dots ,{\mathit{x}}_{\mathrm{T}}-{\mathit{x}}_{\mathrm{T}-1}\right]\\ =\mathit{X}\left[\begin{array}{cccc}-1& & & \\ 1& -1& & \\ & 1& \ddots & \\ & & \ddots & -1\\ & & & 1\end{array}\right]=\mathit{X}\mathit{D}\end{array}$$

(4)

where $\mathbf{D}$ is a $\mathrm{T}\times (\mathrm{T}-1)$ matrix used for compact notation. The TDG method recovers missing entries in $\mathit{X}$ by maximizing the smoothness of the temporal difference signal $S\left(\mathit{X}\mathbf{D}\right)$ rather than the smoothness of the original signal $S\left(\mathit{X}\right)$. Mathematically, the reconstruction problem can be interpreted as a constrained optimization problem; determining the optimal $\mathit{X}$ that minimizes $S\left(\mathit{X}\mathbf{D}\right)$ (maximize smoothness) while keeping the sampled entries in $\mathit{X}$ unchanged. This can be formulated as follows:

$$\begin{array}{cc}\hfill \underset{\mathit{X}}{\min }f\left(\mathit{X}\right)& =S\left(\mathit{X}\mathbf{D}\right)={\displaystyle \frac{1}{2}}tr\left({\left(\mathit{X}\mathbf{D}\right)}^{\mathrm{T}}\mathit{L}\left(\mathit{X}\mathbf{D}\right)\right)\hfill \\ \hfill \mathrm{s}.\mathrm{t}. \mathbf{Y}& =\mathbf{J}\circ \mathit{X}\hfill \end{array}$$

(5)

where $f\left(\mathit{X}\right)$ is the objective function, $\mathit{X}$ is the decision variable, and $\mathbf{J},\mathbf{Y},\mathbf{D},\mathrm{a}\mathrm{n}\mathrm{d} \mathit{L}$ are given constant matrices. The constraint $\mathbf{Y}=\mathbf{J}\circ \mathit{X}$ represents the sampling of $\mathit{X}$, where $\mathbf{J}$ is a sampling operator of the same size as $\mathit{X}$, and $\circ$ denotes the entry-wise product. Entries of $\mathbf{J}$ are equal to 1 if the corresponding entries in $\mathit{X}$ are observed, and equal to 0 if the corresponding entries in $\mathit{X}$ are missing. In other words, $\mathbf{Y}$ represents the observed or reliable part of $\mathit{X}$, which is determined to be fixed. Since $\mathit{L}$ is positive-semidefinite and $f\left(\mathit{X}\right)$ is in a variant quadratic form, it can be inferred that $f\left(\mathit{X}\right)$ is a convex (downward) function. Qiu et al. [25] solved this optimization problem using a gradient projection method, employing a backtracking line search approach to determine the step size. According to convex optimization theory, this algorithm theoretically converges to the global optimal solution after sufficient iterations.

The TDG method has great potential for reconstructing NDVI data because it aligns with a key characteristic of NDVI data—the short-term changes in NDVI between adjacent regions tend to be highly consistent. Even in heterogeneous areas, the temporal differences in NDVI often show greater smoothness than the original data. This phenomenon likely arises because weather conditions—the primary factor influencing vegetation growth—tend to exhibit spatial smoothness and act as a “soft” filter applied to the original NDVI data.

2.3. Applying the TDG Model to NDVI Reconstruction

2.3.1. Grid-to-Graph Conversion

To apply the TDG method to NDVI reconstruction, the NDVI data must first be interpreted as a time-varying graph signal. As the most widely used NDVI products, MODIS data are commonly applied in global evaluations of reconstruction methods [31]. MODIS NDVI data are typically gridded, with each pixel positioned at the cells of a regular grid. Hence, we propose a method to establish graph structures on the grid. In this method, each pixel of the original grid serves as a vertex of the target graph directly. To determine the edges, we used a k-nearest-neighbor rule, linking each vertex only to its k nearest neighboring vertices. Figure 2 shows two examples (with four and eight nearest neighbors). The weights of edges are determined by the distances between connected vertices. When two connected vertices are close, they tend to exhibit high similarity in NDVI variability, so we assign a larger weight to their connecting edge. By default, we set the weight of each edge to be inversely proportional to the distance between the two connected vertices. Particularly, when using the four-nearest-neighbor rule, all edges have the same weight.

Accordingly, a weighted graph $G=(V,E,\mathit{W})$ is established on the original grid, and the NDVI data can be viewed as a time-varying graph signal represented by a matrix $\mathit{X}\in {\mathbb{R}}^{\mathrm{N}\times \mathrm{T}}$. Each column of $\mathit{X}$, denoted by ${\mathit{x}}_t$, represents the spatial profile of NDVI over the region of interest at a specific time instance $t (1\le t\le \mathrm{T})$.

2.3.2. Sampling of the NDVI Signal

The NDVI time-varying graph signal is partially sampled due to the presence of noise and data gaps, which we treat as missing entries in a graph signal. Thus, the NDVI reconstruction problem is framed as a time-varying graph signal reconstruction problem. This is further formulated as a constrained optimization problem in Equation (5) using the TDG method. The optimal solution to this problem yields the final reconstructed NDVI data. In this approach, spatial correlations within NDVI data are modeled by the graph structure, while temporal correlations are captured by the temporal difference signals. The utilization of spatial and temporal information for NDVI reconstruction is global, asymmetric, and reasonable.

The identification of noise in NDVI data is not within the scope of this study. Nowadays, many NDVI products provide pixel-level quality assessment information along with the NDVI data. It is straightforward to distinguish between reliable observations and noise based on such information.

2.3.3. Algorithm to Solve the Problem

We employed the gradient-projection algorithm with backtracking line search, as described by Qiu et al. [25], to solve the problem. Given $\mathbf{J},\mathbf{Y},\mathit{L}$, and an initial signal ${\mathit{X}}_0$, the algorithm updates the signal along the negative gradient direction step by step, gradually approaching the optimal solution after sufficient iterations. For a detailed description of the procedures of this algorithm, readers can refer to Qiu et al. [25]. There is one slight difference in our approach; instead of initializing all noisy entries in the NDVI signal to zero, as this algorithm normally does, we used the preprocessed NDVI values (described in Section 3) as the initial signal ${\mathit{X}}_0$. This benefits solving the optimization problem by providing prior estimations of the true NDVI values.

The graph involved in this study is generally large in scale because graphs built on a NDVI grid are very dense. Specifically, we are dealing with graphs consisting of $128\times 128$ vertices and even more edges, which is much larger in scale than the application examples demonstrated by Qiu et al. [25], by about two orders of magnitude. Therefore, solving the problem could be time consuming. To accelerate the computation, we adopted two strategies in practice: First, we utilized vector/matrix operations instead of entry-wise operations as much as possible in our implementation of the solving algorithm. Such operations can be significantly accelerated using Basic Linear Algebra Subprograms (BLAS, see Blackford et al. [32]). Second, we simplified the formula for calculating the smoothness of temporal difference signals. The counterpart of Equation (2) for temporal difference signals is as follows:

$$S\left(\mathit{x}\right)={\displaystyle \frac{1}{2}}\sum _{1\le t<\mathrm{T}}\sum _{e^{(i,j)}\in E}{W}_{i,j}{\left((x_{t+1,i}-x_{t,i})-(x_{t+1,j}-x_{t,j})\right)}^2$$

(6)

where $x_{t,i}$ denotes the signal sample on vertex $v_i$ at time instance $t$. Note that the smoothness related to reliable entries is always constant during optimizing iterations, whereas what really matters is only the smoothness related to noise. Thus, we only need to calculate the noise-related portion of smoothness, that is

$$S\left(\mathit{x}\right)={\displaystyle \frac{1}{2}}\sum _{1\le t<\mathrm{T}}\sum _{e^{(i,j)}\in E_t^{*}}{W}_{i,j}{\left((x_{t+1,i}-x_{t,i})-(x_{t+1,j}-x_{t,j})\right)}^2$$

(7)

where $E_t^{*}$ is a subset of the entire edge set $E$. Each of the edges in $E_t^{*}$ has at least one end vertex that is diagnosed as noise at time instance $t$ or $t+1$. By doing so, the quantity of computation is reduced substantially.

2.4. Savitzky–Golay Filter Method

For comparison purposes, we implemented a time-series-based NDVI reconstruction method known as the SG filter method [13]. While the basic idea of the SG filter method is to smooth NDVI time series, Chen et al. [13] designed an iterative procedure to replace incompatible depressed values in the NDVI time series with corresponding larger values generated by SG filtering. Since its proposal, the SG filter method has been widely used to preprocess NDVI data in numerous studies [33,34,35].

Recent studies comparing various NDVI reconstruction methods suggest that the SG filter method represents the current state-of-the-art. Michishita et al. [36] compared seven NDVI reconstruction methods and identified the SG filter method and the RMMEH method [17] as the top choices for the Poyang Lake area in China. Zhou et al. [37] revealed that the SG filter outperformed four other methods in tropical and subtropical regions, while the AG method [12] excelled in high latitude boreal regions. Overall, although no single method consistently emerges as superior, the SG filter method is widely recognized as a state-of-art approach. In this study, the SG filter method serves as the baseline model for NDVI reconstruction. We aimed to interpret the performance comparison between the TDG method and the SG filter method as a comparison between an advanced spatio-temporal method and the prevailing time-series-based methods.

2.5. Model Evaluation

To assess the performance of the NDVI reconstruction methods, particularly to compare the TDG method and the SG method, we designed a straightforward evaluation framework. First, some artificial noise was introduced into the original NDVI data. Subsequently, each NDVI reconstruction method was separately used to recover the noise-introduced NDVI data. The results were then evaluated by comparing the recovered NDVI values with their original values, considered as ground-truthing. The greater the proximity between the recovered values and the original ones, the better the method’s performance.

In the noise introduction phase, random entries from the reliable data were selected to be contaminated by artificial noise. These chosen entries were then assigned new random values following a uniform distribution within the valid range. Moreover, these entries were labeled as noise (by changing the Pixel Reliability Index (PRI) value as shown in

Table 1. PRI values for MODIS NDVI Data.

Rank Key	Summary Quality Assurance (QA)	Description
−1	Fill/No Data	Not processed
0	Good Data	Use with confidence
1	Marginal Data	Useful, but check detailed QA for more information
2	Snow/Ice	Target covered with snow/ice
3	Cloudy	Target not visible, covered with cloud

This table is extracted from the MODIS Vegetation Index User’s Guide [38].

), ensuring they are targeted for recovery by the NDVI reconstruction methods. To evaluate the accuracy of the reconstructed results, we employed the root-mean-square error (RMSE) metric. The RMSE between the reconstructed NDVI data matrix ${\mathit{X}}_r$ and the original NDVI $\mathit{X}$ is defined as

$$\mathrm{R}\mathrm{M}\mathrm{S}\mathrm{E}={{\displaystyle \frac{1}{\sqrt{\mathrm{s}}}}‖{\mathbf{S}\circ (\mathit{X}}_r-\mathit{X})‖}_F$$

(8)

where the noise-introduction operator $\mathbf{S}$ is a matrix whose entries equal 1 if the corresponding NDVI entries are contaminated with artificial noise, or, otherwise, equal to 0. The scalar $\mathrm{s}$ is the number of non-zero entries in $\mathbf{S}$, namely, the number of noise-introduced entries. The notation ${‖·‖}_F$ denotes the Frobenius norm of a given matrix.

3. Experiment

3.1. Dataset

We utilized the MODIS/Terra Vegetation Indices 16-Day L3 Global 250 m Grid (MOD13Q1) products (version 006) [38], acquired from NASA’s Land Processes Distributed Active Archive Center (LP DAAC). The MOD13Q1 dataset includes two vegetation index products—the standard NDVI and the Enhanced Vegetation Index (EVI). For our study, we only used the NDVI product. The MOD13Q1 products are globally generated (terrestrial area only) at a spatial resolution of 250 m with 16-day intervals. In this study, the downloaded dataset spans from 18 February 2000, to 17 January 2017, comprising a total of 390 time-series images. These images were mapped in the Sinusoidal Grid projection, an equal-area map projection. We maintained the original projection throughout the whole processing to avoid introducing new distortions into the original data.

The MOD13Q1 products include QA information along with the NDVI data, which assesses the quality of the NDVI data at the pixel level (Table 1). This QA information provides details about various factors affecting data quality, such as aerosol quantity and cloud conditions, enabling users to evaluate the reliability of the data they are using. A summary QA layer, known as PRI (see [39]), is attached to each NDVI data file. Note that the PRI information was used to distinguish between reliable NDVI data and noise. We supposed that data with a PRI value of 0 (indicating “good data”) were reliable, while data labeled as “fill/no data”, “marginal data”, “snow/ice”, and “cloudy” were all regarded as unreliable. This approach allowed us to quickly infer the sampling operator $\mathbf{J}$ and the observed matrix $\mathbf{Y}$ in Equation (5) based on the PRI.

3.2. Study Area

As depicted in Figure 3, the study area is located in the Three-River Headwaters Region, Northwest China, encompassing latitudes ranging from 30°N to 39°N and longitudes ranging from 87°E to 105°E. The land cover information shown in this figure is derived from the GlobeLand30-2010 dataset, a global land cover (GLC) product with a resolution of 30 m [40]. This region exhibits a spatial transition of land cover types from forests in the southeast to bare land in the northwest. Between these two distinct regions, grass dominates the vegetation cover across the vast area. The spatial variation of vegetation in this region accords well with the distribution of precipitation. Annual average precipitation decreases from approximately 600 mm to 100 mm northwestwards. Consequently, the region has an arid or semi-arid climate. Moreover, the region serves as the source of three major rivers in China—the Yangtze River, the Yellow River, and the Mekong River (known as the Lantsang River in China), which highlights the significance of vegetation monitoring and ecological assessment in this region [40].

We selected ten sample patches distributed across the study area for testing the NDVI reconstruction methods. Each NDVI patch, extracted from the MOD13Q1 NDVI product, has a size of 128 × 128 pixels. Given that each pixel corresponds to an actual length of 231.66 m in a 250 m Sinusoidal Grid, each sample patch has an area of about 880 km². Along with the NDVI data, the PRI data for these sample patches were also extracted. These 10 sample patches were chosen to cover the best variability in land cover, geography, and NDVI data quality across the study area.

Table 2. Basic Information for 10 Selected Sample Patches.

Patch Number	Central Longitude	Central Latitude	Average Elevation (m)	Annual Average NDVI	Good-Data Rate
1	92°23′E	37°42′N	2897	0.05	15.53%
2	94°42′E	36°21′N	3185	0.06	33.15%
3	101°46′E	32°25′N	3784	0.45	49.92%
4	95°30′E	32°17′N	4638	0.31	47.46%
5	87°36′E	30°31′N	5223	0.16	57.56%
6	102°21′E	34°47′N	3531	0.40	64.77%
7	95°52′E	34°16′N	4611	0.22	47.26%
8	100°31′E	38°12′N	3389	0.19	67.80%
9	103°29′E	31°35′N	2885	0.53	41.37%
10	105°03′E	36°54′N	1764	0.15	61.55%
Mean			3591	0.25	48.46%

Note: The average elevations were calculated using the Advanced Spaceborne Thermal Emission and Reflection Radiometer (ASTER) Global Digital Elevation Model Version 2 (GDEM V2) dataset [41]. Annual average NDVI was calculated based on the MOD13Q1 NDVI product for the period 2001–2016 without the impact of “Fill/No Data” entries.

provides statistics on the patches with regard to properties including average elevation, annual average NDVI, and the percentage of NDVI entries marked as “good data” among all entries in a patch. From these statistics, the region of interest was situated in a highland characterized by poor vegetation, with the NDVI data in this region suffering from severe noise problems. Geographically, the Qinghai-Tibet Plateau lies southwest of this region, resulting in mostly high altitudes for the selected sample patches, with some exceeding 4000 m. The annual average NDVI values for these patches are predominantly below 0.5, with some extreme cases around 0.1, indicating sparse vegetation cover. Moreover, the quality of NDVI data varies spatially, with an average “good data” rate of 48.46%, which means less than half of the entire NDVI data can be used with full confidence. Some patches, such as Patch 1, exhibit an extremely low “good data” rate, falling below 20%. Such an unfavorable condition amplifies the difficulty of NDVI reconstruction.

3.3. Data Preprocessing

Before executing the NDVI reconstruction using either the TDG or SG filter method, we preprocessed the data to eliminate illogical noisy values through linear interpolation along the time dimension. The PRIs accompanying the NDVI data serve as a confidence criterion for determining which part of the data needs interpolation. While all data excluding “good data” are uniformly considered as noise, some noisy data still contain valuable information for reconstruction, especially those marked as “marginal data”. Therefore, we only interpolated extremely noisy entries marked as “snow/ice”, “cloudy”, or “fill/no data” based on NDVI values at other high-confidence entries marked as “good data” or “marginal data”.

It is worth noting that the first step of the SG filter method typically involves linearly interpolating cloudy NDVI values [13], which is very similar to the aforementioned interpolation procedure. Consequently, in our implementation of the SG filter method, we omitted the first step to avoid redundant operations. Thus, the SG filter method used in this study differs slightly from the original version. Specifically, we interpolated NDVI values marked as “snow/ice”, “cloudy”, and “fill/no data” following the replacement method employed by Chen et al. [13]. Another important difference is that the original SG filter method also interpolates NDVI points that increase greater than 0.4 over 20 days, even though they are not identified as noise. However, we chose to discard this rule in our implementation, opting for simplicity in identifying NDVI noise based solely on provided information (the PRI data) without this additional rule. This approach ensures consistency in data preprocessing between the TDG method and the implemented SG filter method in this study, maintaining fairness in the comparison between them.

3.4. Model Configuration

We applied the TDG and SG filter method to reconstruct the preprocessed NDVI data on the selected 10 sample patches over the study area. For the TDG method, we used the four-nearest-neighbor rule to build edges in the graph. The initial learning rate was set to 0.1, and the iterations were stopped after a maximum of 300 steps or when no further perceptible improvement was observed. This configuration was found to produce good reconstruction results in practice. For the SG filter method, we used the same parameters as recommended by Chen et al. [13].

3.5. Details of Model Evaluation

To evaluate the performance of the two methods, NDVI reconstruction was conducted within the evaluation framework described in Section 2.5. On each sample patch, 100,000 artificial noisy entries were randomly introduced into the original NDVI dataset (approximately 1.57% of the total data). For a fair comparison, both the TDG method and the SG filter model shared an identical selection of noise-introduced entries and identical new values as well. Note that the noise introduction was accomplished before data preprocessing. In addition, considering that NDVI reconstruction methods might respond differently to various types of noise, we explored three different ways to introduce artificial noise, as follows: (1) Positive-Marginal (PM) artificial noise, introduced by replacing the values of selected NDVI entries with random values greater than the original ones and identifying their PRIs as “marginal data”, (2) Negative-Marginal (NM) artificial noise, introduced by replacing the values of selected NDVI entries with random values lower than the original ones and also identifying their PRIs as “marginal data”, and (3) No-Data (ND) artificial noise, introduced by identifying the PRIs of selected NDVI entries as “fill/no data”. The NDVI values do not need to be changed when introducing ND noise because during data preprocessing, the values of these entries will eventually be linearly interpolated. The introduced PM and NM noise tested the performance of the NDVI reconstruction methods in handling positive and negative deviations in NDVI, respectively. The introduced ND noise tested the performance of recovering linearly interpolated NDVI values, which is relatively simpler. Each of the three types of artificial noise was introduced independently, meaning that we carried out three separate rounds of performance evaluation for the NDVI reconstruction methods. In each round, the noise-introduced NDVI data were obtained by contaminating 100,000 randomly selected entries with a single type of artificial noise. As such, this approach allowed us to examine the sensitivity of the NDVI reconstruction methods to different characteristics of noise.

4. Results and Discussion

4.1. Comparison of the TDG Method and the SG Method

The RMSEs of the reconstructed NDVI for each sample patch using the TDG method or the SG filter method are listed in

Table 3. RMSEs of Reconstructed NDVI for Each Sample Patch.

Patch	SG Filter			TDG
	PM	NM	ND	PM	NM	ND
1	0.501	0.035	0.006	0.007	0.004	0.004
2	0.483	0.038	0.021	0.021	0.020	0.020
3	0.279	0.111	0.051	0.029	0.030	0.029
4	0.352	0.089	0.050	0.032	0.032	0.032
5	0.448	0.052	0.026	0.019	0.019	0.018
6	0.320	0.085	0.044	0.027	0.027	0.027
7	0.388	0.079	0.047	0.033	0.033	0.033
8	0.422	0.060	0.038	0.026	0.026	0.027
9	0.232	0.148	0.070	0.047	0.047	0.047
10	0.449	0.045	0.027	0.021	0.021	0.021
Average	0.387	0.074	0.038	0.026	0.026	0.026

, arranged by different types of artificial noise. The results reveal the substantial superiority of the TDG method over the SG filter method. The TDG method not only achieves much lower RMSEs than the SG filter method, in most cases, but also exhibits a remarkable stability in performance when processing various types of noise. However, the SG filter method is highly sensitive to the type of noise introduced. The average RMSE of reconstructed NDVI across the 10 patches using the TDG method is approximately 0.026, regardless of the type of artificial noise introduced. Meanwhile, the average RMSEs for the SG filter method are 0.387, 0.074, and 0.038 for PM, NM, and ND noise, respectively. For a more comprehensive assessment, we also calculated the R² and Mean Absolute Error (MAE), both of which indicated the superiority of the TDG method.

For the three types of artificial noise, Figure 4 sequentially presents one example per type of reconstructed NDVI time series using the two methods. The varying background colors in the set of figures indicate the PRI categories of the corresponding entries. As illustrated in Figure 4a, the SG filter method fails to reconstruct PM noise entirely. The reconstructed NDVI time series based on the SG filter method (green line) retains false spikes at noise-introduced entries, leading to an unacceptable high RMSE (0.387) compared to the original true values (red dotted line). Note that if we use the standard version of the SG filter method that forcibly interpolates unnatural spikes in the first step, the extreme spike occurring in 2014 can be wrapped out, but the two mild spikes in 2013 remain the same because they are below the criterion to interpolate (increasing by more than 0.4 over 20 days). In contrast, the TDG method effectively eliminates these false spikes, and the reconstructed values (blue line) are very close to the true values. Since those sudden spikes seriously disrupt the NDVI smoothness both in time and space, the TDG method forces them to flatten out. We would like to mention that the failure of the SG filter method regarding PM noise is expected because this method is designed based on the assumption that noisy NDVI values are usually lower than true values. This assumption would generally hold because the overwhelming noise introducers like cloud, snow, and low sun zenith angles usually cause negative deviations from true values; however, we cannot reject the existence of positive deviations in NDVI arising from unusual reasons. Indeed, we commonly observed spurious rises in the original NDVI (see Figure 5). Therefore, there is practical meaning in inspecting the behavior of NDVI reconstruction methods when facing positively biased noise.

The noise in real NDVI data is mostly negatively biased, so the NM artificial noise can better simulate the real situation of noise contamination in NDVI than the other two types. From the example shown in Figure 4b, both the TDG method and the SG filter method can recover NM noise, producing reconstructed NDVI time series that have well eliminated artificially introduced sudden drops. However, it is noticeable that the reconstructed time series based on the SG filter method (green line) shows a slight but systematic negative bias with respect to the true values (red dotted line) at noise-introduced entries. In contrast, reconstructed NDVI values based on the TDG method (blue line) are generally closer to the true values. The RMSEs of the two reconstructed time series quantify the performance gap between the two methods. The 10-patch average RMSE for the SG filter method is 0.074 but drops to 0.026 for the TDG method, indicating a dramatic improvement.

Regarding ND noise, Figure 4c shows that the reconstructed NDVI values at noise-introduced entries are very close to the true values for both methods. From the RMSE perspective, the TDG method still performs better than the SG method in recovering ND noise, with much smaller RMSEs at all sample patches except Patch 2.

To assess how well the two methods work in real NDVI reconstruction practice, we executed another round of NDVI reconstruction which was free from the evaluation framework, i.e., using the two methods to reconstruct the original NDVI directly without introducing artificial noise. Figure 5 demonstrates four sets of reconstructed results concerning NDVI time series with various good-data rates (i.e., the rate of “good data” entries along the whole time series). When the NDVI is of good quality, such as the time series in Figure 5a (the good-data rate is 76.9%), both methods effectively eliminate sudden drops in NDVI. As data quality degrades, the superiority of the TDG method becomes evident in the following aspects: (1) The TDG model can fix positively biased noise, as previously discussed. A perfect example is the NDVI entry spiking during the spring of 2015 in Figure 5b. (2) The TDG method better recovers the periodic intra-annual variability of NDVI, as demonstrated in Figure 5c,d. Notably, the TDG model successfully recovers the summer peak of NDVI in 2014 in Figure 5d, where the SG filter method fails because several low-value pixels during this period are categorized as “marginal data” but not “cloudy” and, thus, are not interpolated. (3) If the NDVI time series experiences long-term continuous contamination with bad noise, the TDG model can recover the underlying variability of the NDVI, while the SG model tends to inherit the straight slopes from the linear interpolation in the data preprocessing. Typical cases of this can be found in early 2003 and 2012 in Figure 5b. These advantages of the TDG method are attributed to the utilization of both temporal and spatial interrelations within NDVI data. Note that the TDG model only modifies noisy NDVI entries, so the reconstructed NDVI values at “good data” entries are strictly identical to the original values. Accordingly, this explains why RMSE is lower and constant when introducing more noise-affected observations to be corrected.

Figure 6 also demonstrates the reconstructed NDVI results but from the spatial perspective. Four sets of NDVI spatial profiles with varying good-data rates (i.e., the rate of “good data” entries over the spatial profile) are shown. In Figure 6a, where the case is mostly covered by “good data” entries, the NDVI spatial profiles remain the same after reconstruction. In the other three cases, the quality of NDVI significantly improves after reconstruction. The NDVI reconstruction methods can remove noisy spots that are pervasively scattered over the original NDVI spatial profiles, producing new profiles that are much “cleaner”. Particularly, in Figure 6d, the original NDVI profile has no “good data” entries and presents a very flat and uninformative appearance. Nevertheless, the spatial details and fine textures of this profile can be successfully rebuilt after reconstruction. Comparing the reconstructed NDVI profiles based on the TDG method and those based on the SG filter method highlights the superiority of the TDG method; while the SG filter method significantly improves the quality of NDVI data, it still leaves some coarse spots over the reconstructed profiles. This method reconstructs NDVI time series individually per pixel, which may not guarantee a coherent spatial pattern in the reconstructed NDVI profiles. In contrast, TDG-based reconstruction successfully eliminates these coarse spots and restores the unsmooth structures of the NDVI spatial profile. This is because adjacent NDVI pixels are connected via graph structures, enabling the TDG model to obliterate noisy entries where temporal difference signals conflict with nearby pixels.

4.2. Robustness and Further Applications of the TDG Method

The robustness of the TDG method was further examined by focusing on its performance when noise is extremely prevalent in NDVI data. To achieve this, we incrementally increased the level of introduced noise and checked the RMSE of reconstructed NDVI based on the aforementioned evaluation framework. This process was applied to Patch 6, which had the highest good-data rate (64.77%) among the sample patches. For this test, we only introduced NM artificial noise to generate NDVI data with varying good-data rates, defined as the proportion of “good data” entries within the overall spatio-temporal NDVI data block. Figure 7 illustrates the relationship between the RMSE of reconstructed NDVI and the good-data rate after noise introduction, tested on Patch 6. The results demonstrate that the TDG method is highly robust to noise density. The RMSE of reconstructed NDVI remains largely constant if the good-data rate is no less than 20%. Even when the good-data rate drops to 10%, the TDG method still achieves a relatively low RMSE (0.041). In contrast, the SG filter method’s RMSE increases linearly as the good-data rate decreases, failing at an early stage when the good-data rate becomes too low.

We also explored the impact of different edge-establishment strategies on the performance of the TDG method’s reconstruction. While previous analysis consistently used the default four-nearest-neighbors rule, we tested several other rules (8-, 12-, and 20-nearest-neighbor) to evaluate their effect on NDVI reconstruction. The results indicated that while denser edge-establishment rules enhance the connectivity of the graph, the RMSE of reconstructed NDVI remains almost constant. Specifically, the 10-patch average RMSE decreases only marginally from 0.026 using the four-nearest-neighbor rule to 0.025 when using eight or more nearest neighbors. This marginal improvement is because the four-nearest-neighbor graphs are already fully connected, meaning each pair of vertices in the graph are associated via one or more edges. Adding more edges has a minimal effect since the influence of one vertex on another diminishes with increasing distance between them. Therefore, the simplest four-nearest-neighbor rule is considered the optimal choice for edge establishment in the TDG method, which not only produces good results but also reduces computation.

4.3. Limitations of the TDG Method and Future Directions

To reconstruct high-quality NDVI time series derived from satellites, plenty of algorithms have been proposed to smooth spatio-temporal NDVI data. Hybrid methods, i.e., the temporal spatial filter (TSF), the search-and-fill algorithm with moving offset method (SFA-MOM), the spatio-temporal Savitzky–Golay (STSG) method, and the spatio-temporal tensor completion (ST-Tensor) method, are considered to overcome the limitation of a single-feature-dependent approach by incorporating the information of spatial correlation and temporal continuity [1]. Spatial-based methods and temporal-based methods may have good performance under specific circumstances but not in others [42]. Although hybrid methods can tackle long-term data gaps and preserve valid low values, they may also have poor performance in the case of high-spatial-resolution data compared to low-temporal-resolution data, and uncertainties may be introduced due to the low accuracy of reference data [1,43,44]. Distinct from previous NDVI reconstruction methods, the proposed TDG method benefits from its formalized mathematical approach, which eliminates the need for manually overdesigned procedures that rely on prior knowledge or additional assumptions. While the present results are encouraging, this study’s evaluation of NDVI reconstruction methods is somewhat simplified. The TDG method was tested on a region with relatively uniform vegetation types and environmental conditions, was compared with only one time-series-based method, and was evaluated using a single metric (RMSE). As the number of grids and the length of the time series increase, computational requirements will grow despite efficient optimization strategies. This raises concerns about the operational applicability and the optimal spatial scale for this method. Therefore, a more detailed and systematic evaluation is needed to comprehensively investigate the TDG method’s performance.

There are a few promising directions for further improving the TDG method for NDVI reconstruction: (1) Incorporating measurement noise. A variant of the TDG method considers measurement noise [25], which not only reconstructs the low-quality NDVI values explicitly identified as noise but also rectifies some suspect values, though they are marked as reliable data. This variant provides a flexible and promising alternative for NDVI reconstruction but has not been well studied. (2) Optimizing edge weights. In this study, edge weights were completely determined by the distance of two connected vertices, following the simple assumption that vertices close to each other share high similarity in NDVI variability. In fact, the similarity between vertices can be quantitatively measured using historical data, e.g., by calculating the correlation coefficient between NDVI time series at two vertices. Incorporating this information to determine edge weights would enhance the precision of spatial intercorrelation modeling within NDVI data. (3) Fusing multi-source datasets of different resolutions. Utilization of multiple satellite-based datasets would provide additional information to offset the deficiency of spatial or temporal information. Recent studies have developed several spatio-temporal fusion models to produce improved NDVI reconstructed datasets with both high spatial and temporal resolutions [4,45], but various types of noise associated with distinct datasets still need to be clearly addressed, and large uncertainties may remain. It is worth mentioning that machine learning techniques will gain increasing advantages in the reconstruction of NDVI time series, which are being accumulated, such as deep convolutional neural networks (DCNNs), long short-term memory (LSTM) network, and spatial-temporal-spectral deep convolutional neural network (STSCNN).

5. Conclusions

This study proposed a TDG method for NDVI reconstruction and demonstrated its usefulness in the Three-River Headwaters Region of Northwest China. By converting NDVI data into a time-varying graph signal, the TDG method leverages an optimization algorithm to maximize the spatial smoothness of temporal differences in the original signal. A comparative performance analysis between the TDG method and the SG filter method was conducted, which involved introducing artificial noise into the data and evaluating the RMSE between the recovered NDVI values and the actual values. The results revealed a substantial superiority in the TDG method over the SG filter method for NDVI reconstruction. The average RMSE of reconstructed NDVI using the TDG method was 0.026, regardless of whether the noise was positively biased, negatively biased, or linearly interpolated. In contrast, the SG filter method’s average RMSEs were 0.387, 0.074, and 0.038 for the three types of noise, respectively. Moreover, the TDG method exhibited remarkable robustness to noise intensity, requiring only 10% to 20% of the data to achieve high-accuracy reconstruction. Thus, the TDG method shows excellent performance with no preference for noise types and yields good results from both a time-series perspective and spatial-profile perspective.

The TDG method is a very useful tool for reconstructing noisy NDVI products, primarily due to its comprehensive utilization of both the temporal and spatial information in NDVI data. This is achieved through two key aspects—modeling spatial interrelations between data pixels using a graph structure, and exploiting temporal associations of NDVI values through temporal differences, based on the assumption that the temporal-difference signal tends to be spatially smooth. Moreover, the TDG method exhibits universal applicability in its performance for NDVI reconstruction, regardless of noise characteristics, land cover types, vegetation phenology patterns, and other factors. Thus, the proposed method’s potential for extended application to datasets beyond NDVI is highly prospective, and we hope our efforts could provide new insights into remote-sensing data processing. Note that the principle of the TDG method is also applicable to a wide range of remote-sensing-based products. Nevertheless, many issues still need to be addressed. Considering the spatial smoothness in the TDG method, edge effects due to both the temporal and spatial borders would introduce additional uncertainty. Although the computation efficiency in this study is relatively high due to the use of a high-performance server, more tests and simulations are necessary to demonstrate the operational applicability of the TDG method. Furthermore, various combinations of different types of noise should be investigated, including measurement noise, and, also, the principle, configuration, and algorithm of the proposed approach need to be further improved and validated based on various datasets and ground observations in the future.