Impacts of the data quality of remote sensing vegetation index on gross primary productivity estimation

ABSTRACT

As the most commonly used driven data for gross primary productivity (GPP) estimation, satellite remote sensing vegetation indexes (VI), such as the leaf area index (LAI), often seriously suffer from data quality problems induced by cloud contamination and noise. Although various filtering methods are applied to reconstruct the missing data and eliminate noises in the VI time series, the impacts of these data quality problems on GPP estimation are still not clear. In this study, the accuracy differences of the GPP estimations driven by different VI series are comprehensively analyzed based on two light use efficiency (LUE) models (the big-leaf MOD17 and the two-leaf RTL-LUE). Four VI filtering methods are applied for comparison, and GPP data across 169 eddy covariance (EC) sites are used for validation. The results demonstrate that all the filtering methods can improve the GPP simulation accuracy, and the SeasonL1 filtering method exhibits the best performance both for the MOD17 model (∆R² = 0.06) and the RTL-LUE model (∆R² = 0.07). The reconstruction of the key change points in the temporally continuous gaps may be the primary reason for the different performance of the four methods. Moreover, the effects of filtering processes on GPP estimation vary with latitudes and seasons due to the differences in the primary data quality. More significant improvements can be observed during the growing season and in the regions near the equator, where the data quality is relatively poor with lower primary GPP estimation accuracy. This study can guide the preprocessing of the VI data before GPP estimation.

KEYWORDS:

• Filtering method
• gross primary productivity
• vegetation index
• light use efficiency model
• terrestrial ecosystem

1 Introduction

Gross primary productivity (GPP), a measurement of the amount of organic carbon produced by vegetation photosynthesis, stands for the amount of carbon dioxide in the atmosphere fixed by vegetation per unit of space and time (Beer et al. 2010; Guan et al. 2021; Lin et al. 2022). GPP is a crucial indicator for vegetation growth and photosynthesis activity, which is helpful for advancing the understanding of the global carbon cycle. Therefore, accurate estimation of GPP at regional or global scales is particularly important under more and more severe global climate change (Ballantyne et al. 2012; Peng et al. 2022).

Currently, there are many methods to obtain GPP, which can be broadly classified into two categories: in-situ measurement-based and model-based estimates (Xie and Li 2020; Yuan et al. 2014). As a mainstream in-situ measurement-based method, the eddy covariance (EC) flux system enables direct measuring of the carbon exchange flux between the vegetation canopy and the atmosphere, thus providing accurate GPP estimates within a specific space around the tower (Baldocchi et al. 2001; Pastorello et al. 2020). However, the limitation in spatial footprint makes it challenging for large-scale terrestrial ecosystem monitoring (Pastorello et al. 2020; Yu et al. 2006). Therefore, lots of GPP estimation models have been developed in recent years and act as a feasible approach for regional and global studies. Especially with the development of remote sensing (RS) technology, repeated observation of large-area vegetation canopy can be achieved, and RS-based estimation models are going mainstream. There mainly includes three types: statistic-based models, process-based models, and light use efficiency (LUE) models (Guan et al. 2022; Lin et al. 2021). These models commonly utilize satellite-based vegetation indices (VI) as driven data, and lots of leaf area index (LAI), normalized difference vegetation index (NDVI), enhanced vegetation index (EVI), and other VIs products have been produced and applied. In this condition, the LUE models, which can be easily combined with RS VI data, have also become the most popular used models for GPP estimation at large scales, and several big-leaf and two-leaf models have been developed (Guan et al. 2021; He et al. 2013; Running et al. 2004).

Nevertheless, satellite-derived VI data is prone to be contaminated by clouds, aerosols, and other adverse atmospheric compositions, resulting in data quality problems like data missing and noise. These contaminations can cause sudden decreases in VI time series, introducing spatiotemporal discontinuity in the data and hindering the accurate representation of vegetation status and changes, which greatly impede the application of VIs (Huete et al. 2002; Michishita et al. 2014). Since GPP is estimated based on VIs, data quality problems in VIs should also impact the GPP estimations. Therefore, the reconstruction of contaminated information and acquisition of high-quality VI time series is crucial for the precise simulation of GPP (Shen et al. 2015), and numerous temporal reconstruction methods have been proposed to generate high-quality VI time series in the past decades. These methods can be generally grouped into three types according to the way they handle the data (Chu et al. 2021; Shen et al. 2015). The first type is the sliding window filter method, which filters the data within a moving window and the results are highly impacted by the selected filtering parameters (Julien and Sobrino 2010; Ma and Veroustraete 2006; Savitzky and Golay 1964; Viovy, Arino, and Belward 1992). The second type is the function-based curve fitting methods that use mathematical functions to fit the entire time-series VI curve, primarily used to obtain the outer envelope of VI time-series data (Beck et al. 2006; Jonsson and Eklundh 2002; Julien and Sobrino, 2019). The third type is the frequency domain method which reconstructs data and removes noise within the frequency domain (Ghafarian Malamiri et al. 2020; Lu et al. 2007; Verhoef, Menenti, and Azzali 1996). There are also some methods that cannot be classified into the aforementioned three categories, such as variation-based filtering methods like Whittaker (Atzberger and Eilers 2011) and SeasonL1 (Chu et al. 2022). Among them, SeasonL1 combines the regularization of local smoothness and inter-annual seasonal similarity, which can ensure the full utilization of the prior features in the whole time series, and has been demonstrated to outperform several classic filtering methods in reconstructing VI time series (Chu et al. 2022).

Although these filtering methods provide abundant tools to preprocess VI data before GPP estimation, the impact of VI quality problems on GPP estimation remains unclear. Since quality problems in VI are often manifested as abrupt decreases, which should also result in an underestimation in GPP, because the main spatiotemporal changes of GPP may largely come from the RS VI data (Yuan et al. 2014). If there is a false sudden decrease in VI, there should also be an underestimation of GPP at the same time. Consequently, the intermittent underestimation not only makes it difficult to simulate the continuous changes of GPP over time but also hinders the accurate representation of spatial patterns in large regions. However, although this problem is widely recognized, the precise magnitude of this impact has not been evaluated yet. Data filtering methods have been frequently used in pre-processing to improve the quality of remote sensing data, but there has been little in-depth discussion about the underlying mechanism. The impacts of VI quality improvement on GPP simulation have not been quantified in previous research. Furthermore, most current studies merely select a temporal reconstruction method arbitrarily for processing and application, lacking research on the selection of the specific optimal method for GPP estimation. There also lacks corresponding research on the impact of parameter settings in filtering methods. Therefore, it is necessary to systematically investigate the influence of VI quality issues on GPP estimation and understand the effectiveness of different filtering methods in addressing this problem, providing references for subsequent GPP estimation.

In this context, this study systematically quantifies the impacts of RS index filtering on GPP simulation based on the global EC datasets. The main purposes are: (1) to quantify the impact of VI data quality problems on GPP estimation; (2) to compare the performance of different filtering methods; (3) to compare the sensitivity of the big-leaf and two-leaf models to VI data.

2 Data and methods

2.1 Data

2.1.1 EC data

In this study, EC data from the FLUXNET2015 dataset (www.fluxdata.org) is employed to assess the estimation from different LUE models and different VI inputs (Baldocchi et al. 2001). The FLUXNET2015 dataset follows a standardized format and encompasses refined data gathered from various regional networks, covering a total of 212 sites across the globe (Pastorello et al. 2017, 2020). Consistent with prior studies, GPP partitioned using the daytime partitioning (DT) method (referred to as “GPP_DT_VUT_REF”) served as the reference GPP, which is derived from net ecosystem exchange (NEE) through the application of a hyperbolic light response curve (Reichstein et al. 2005). A total of 169 sites (Table S1) with more than one year of valid half-hour observations during the period from 2001 to 2015 are selected. These sites contain 12 vegetation types according to the International Geosphere-Biosphere Program (IGBP) classification system. The number of sites and site-years for each vegetation type are listed in Table 1, and the spatial distribution of sites is displayed in Figure 1.

Figure 1. The spatial distribution of the 169 EC sites.

Table 1. Number of EC sites and site years for different vegetation types in the IGBP classification system.

Vegetation type^a	CRO	CSH	GRA	WET	MF	DBF	DNF	EBF	ENF	OSH	SAV	WSA
Number of sites	13	1	32	18	8	19	1	41	12	11	7	6
Number of site-years	103	4	214	86	76	147	3	326	70	68	45	49

^aCRO: crop; CSH: closed shrub; GRA: grass; WET: wetlands; MF: mixed forest; DBF: deciduous broadleaf forest; DNF: deciduous needleleaf forest; EBF: evergreen broadleaf forest; ENF: evergreen needleleaf forest; OSH: open.

2.1.2 VI data

The Moderate Resolution Imaging Spectroradiometer (MODIS) VI products provide a consistent and reliable means for monitoring global vegetation conditions (Didan et al. 2015). In this study, the MOD15A2H LAI product, provided in an 8-day composition with a spatial resolution of 500 m, is selected to drive LUE models. Moreover, to verify the consistency of results across different indices and enhance their credibility, NDVI and EVI data from the MOD13Q1 product are also selected as the RS-driven data. These two products are synthesized every 16 days and have a spatial resolution of 250 m. These products also include quality assurance (QA) flags that provide information about the quality of the VI products at the pixel level. These data were accessed free from the Land Processes Distributed Active Archive Center (LP DAAC) from the United States Geological Survey (USGS) website on Earth Explorer tool (https://earthexplorer.usgs.gov).

2.2 VI filtering methods

In this study, four commonly used filtering methods are selected to investigate the impact of filtering methods on the simulation accuracy of GPP using LUE models, including Savitzky-Golay (SG), Whittaker, harmonic analysis of time series (HANTS), and SeasonL1.

2.2.1 SG method

The SG method, a representative sliding window filter method, is a local least squares fitting convolution used for smoothing and computing derivatives of a sequence of consecutive values (Chen et al. 2004). A description of the SG filter is as follows (Shen et al. 2015):

(1) ${\hat{x}}_t=\frac{{\sum }_{i=-l}^l\left(b_i{x}_{t+i}\right)}{2l+1}$(1)

where $x_{t+i}$ is the original value of point x acquired at time $t+i$, ${\hat{x}}_t$ is the reconstructed value, $b_i$ is the weight of $x_{t+i}$, and $\left(2l+1\right)$ is the window length of the convolution that determines not only the smoothness of the results but also the filtering efficiency.

2.2.2 HANTS method

The HANTS method (Verhoef, Menenti, and Azzali 1996) is one of the widely used frequency domain filtering methods, its formula can be expressed as (Zhou, Jia, and Menenti 2015):

(2) $\tilde{x}\left(t_j\right)=a_0+\sum _{i=1}^{nf}\left(a_{i}cos\left(2\pi f_i{t}_j\right)+b_{i}sin\left(2\pi f_i{t}_j\right)\right)$(2)

(3) $x\left(t_j\right)=\tilde{x}\left(t_j\right)+\epsilon \left(t_j\right)$(3)

where $x$ represents the original series, $\tilde{x}$ represents the reconstructed series, and $\epsilon$ represents the error series. $t_j$ corresponds to the time at which x is obtained (observed), where j = 1, 2, …, N, and N is the maximum number of observations (samples) in the time series. The symbol $nf$ is the number of periodic terms in the series, i.e. the number of harmonic components associated with the frequencies $f_i$. Generally speaking, these selected harmonics consist of a base frequency and a series of integer multiples of the base frequency. $a_i$and $b_i$ are the coefficients of the trigonometric components with frequencies $f_i$. $a_0$ can be viewed as the coefficient at zero frequency, which is the average of the series. The harmonic fitting can be achieved by solving the equations above using the linear least square method.

2.2.3 Whittaker smoother

The Whittaker filtering method, a commonly used variation-based filtering method (Atzberger and Eilers 2011; Eilers 2003), employs the second-order difference of the time series to obtain smooth results. In particular, when the second-order difference $x_{t+1}-2x_t+x_{t-1}$ approaches 0, it indicates temporal smoothness. As a result, a regularization term is formulated, leading to the creation of the following filtering model.

(4) $\hat{x}=arg\underset{x}{min}\sum _{t=1}^n{\left(x_t-y_t\right)}^2+{\lambda }_1\sum _{t=2}^{n-1}{\left(x_{t+1}-2x_t+x_{t-1}\right)}^2$(4)

where $x$ represents the reconstructed series and $y$ represents the original series. ${\lambda }_1$ is a regularization parameter that controls the contributions of local smoothness.

2.2.4 SeasonL1 method

The SeasonL1 filtering method is developed based on the Whittaker filtering method, by considering the inter-annual seasonal similarity of VI time series data. This method has been demonstrated to offer significant advantages in terms of recovering temporally contiguous missing values and preventing over-smoothing in the inflection points of vegetation growth (Chu et al. 2022). The SeasonL1 model can be described as:

(5) $\hat{x}=\underset{x}{argmin}\frac{1}{2}\sum _{t=1}^T{w}_t{\left(x_t-\hat{x}\right)}^2+{\lambda }_1\sum _{t=2}^{n-1}\left|x_{t+1}-2x_t+x_{t-1}\right|+{\lambda }_2\sum _{t=1}^{n-T}\left|x_t-x_{t+T}\right|$(5)

where $w_t$ represents the weight of the corresponding VI point determined by data quality, while ${\lambda }_2$ is the regularization parameter that controls the contributions of inter-annual seasonal similarity, $T$ indicates the length of the one-year time series. This formula consists of three terms: the data fidelity term, the temporal smoothness regularization term, and the inter-annual similarity regularization term, which are respectively used to reduce the deviation between x and y, capture the temporal patterns inherent in VI time series data, and the non-local similarity patterns across multiple years in VI time series.

2.3 GPP estimation models

Current LUE models can be typically divided into big-leaf and two-leaf categories. As a result, this study selects the MOD17 model to represent the big-leaf models (Running et al. 2004) and the RTL-LUE model to represent the two-leaf models (Guan et al. 2021), to analyze whether the VI filtering effects are stable among different types of LUE models.

2.3.1 MOD17 model

The MOD17 model was proposed based on the radiation conversion efficiency concept of Monteith (1972). It can be described as (Running et al. 2004):

(6) $GPP={\epsilon }_{BL}\times PAR\times fPAR\times f\left(VPD\right)\times g\left(T_a\right)$(6)

where ${\epsilon }_{BL}$ is the maximum LUE varied with different vegetation types; $fPAR$ is the fraction of photosynthetically active radiation (PAR) absorbed by the canopy; $fVPD$and $g\left(T_a\right)$ are the scalars of vapor pressure deficit (VPD) and minimum temperature.

(7) $fPAR=1-e^{-k\times LAI}$(7)

(8) $f\left(VPD\right)=\left\{\begin{array}{c}0\hspace{0.17em}VPD\ge VP{D}_{max}\\ \frac{VP{D}_{max}-VPD}{VPD-VP{D}_{min}}\hspace{0.17em}VP{D}_{min}<VPD<VP{D}_{max}\\ 1\hspace{0.17em}VPD\le VP{D}_{min}\end{array}\right.$(8)

(9) $g\left(T_a\right)=\left\{\begin{array}{c}1\hspace{0.17em}{T}_a\ge T_{max}\\ \frac{T_a-T_{min}}{T_{max}-T_a}\hspace{0.17em}{T}_{min}<T_a<T_{max}\hspace{0.17em}\\ 0\hspace{0.17em}{T}_a\le T_{min}\end{array}\right.$(9)

where k is the light extinction coefficient and set as 0.5; VPD is the daily average VPD at day- time and T_a is the daily minimum temperature; VPD_max, VPD_min, T_max, and T_min are parameters specific to vegetation types.

2.3.2 RTL-LUE model

Considering the significant differences in the photosynthetically active radiation (APAR) and LUE of sunlit and shaded leaves (Propastin et al. 2012; Wang et al. 2015), two-leaf models are developed to improve the GPP estimations. The TL-LUE model is the first two-leaf model (Chen et al. 2020; He et al. 2018; Liu et al. 1997), and then, the RTL-LUE model is modified taking the impacts of radiation intensity on LUE into account. The GPP of a plant canopy is calculated as follows:

(10) $GPP={\epsilon }_{RTL}\times \left(f\left(PPF{D}_{su}\right)\times APA{R}_{su}+f\left(PPF{D}_{sh}\right)\times APA{R}_{sh}\right)\times f\left(VPD\right)\times g\left(T_a\right)$(10)

where ${\epsilon }_{RTL}$ represents the maximum LUE observed across all the leaves within the canopy. The scalar of radiation on vegetation photosynthesis in sunlit and shaded leaves is denoted by $f\left(PPF{D}_{su}\right)$ and$f\left(PPF{D}_{sh}\right)$, respectively. $f\left(VPD\right)$and $g\left(T_a\right)$ are calculated using the same method as in the MOD17 algorithm. $APA{R}_{su}$ and $APA{R}_{sh}$ indicate the incoming PAR absorbed by sunlight and shaded leaves, and are computed using the following formulas:

(11) $APA{R}_{su}=\left(1-\alpha \right)\times \left[\frac{PA{R}_{dir}\times cos\left(\beta \right)}{cos\left(\theta \right)}+\frac{PA{R}_{dif}-PA{R}_{dif,u}}{LAI}+C\right]\times LA{I}_{su}$(11)

(12) $APA{R}_{sh}=\left(1-\alpha \right)\times \left[\frac{PA{R}_{dif}-PA{R}_{dif,u}}{LAI}+C\right]\times LA{I}_{sh}$(12)

where $\alpha$ represents the canopy albedo that is associated with specific vegetation types. Meanwhile, $\beta$ corresponds to the mean leaf-sun angle and is set to 60 degrees for a canopy with a spherical leaf angle distribution according to previous research and empirical knowledge (Chen 1996, 1999; Guan et al. 2021; He et al. 2013; Zhou et al. 2016). $\theta$ represents the solar zenith angle. Additionally, $PA{R}_{dif}$, $PA{R}_{dir}$, $PA{R}_{dif,u}$, and $C$ denote the diffuse and direct components of incoming PAR, the diffuse PAR under the canopy, and the multiple scattering of direct radiation (Chen 1999). $LA{I}_{su}$ and $LA{I}_{sh}$ indicate the LAI of sunlit and shaded leaves, respectively, which are separated based on the canopy LAI, clumping index (Ω), and solar zenith angle.

2.4 Model evaluation

The parameters VPD_max, VPD_min, T_max, T_min, α, Ω, ε_max and a are empirically set in LUE models according to the previous literature (Guan et al. 2021, 2022; He et al. 2013; Zhou et al. 2016), which are listed in Table S2. To assess the accuracy of the GPP simulations from the two models against the EC GPP, the coefficient of determination (R²) and the root-mean-square error (RMSE) are used for quantitative evaluation:

(13) $R^2={\left(\frac{{\sum }_{i=1}^N\left(P_i-\overline{O}\right)\left(O_i-\overline{O}\right)}{\sqrt{{\sum }_{i=1}^N{\left(P_i-\overline{P}\right)}^2{\sum }_{i=1}^N{\left(O_i-\overline{O}\right)}^2}}\right)}^2$(13)

(14) $RMSE=\sqrt{\frac{{\sum }_{i=1}^N{\left(P_i-O_i\right)}^2}{N}}$(14)

where N is the total simulated experiment data point; $P_i$ and $O_i$ represent the predicted GPP from the MOD17 or RTL-LUE models and observed EC GPP, respectively; and $\bar{O}$ and $\bar{P}$ are the mean values of observations and predictions for all experimental data points.

3 Results

This section firstly evaluates the accuracy of GPP simulations calculated by different VI data at all 169 sites and analyzed after grouping them into different vegetation types. Then the differences across different LUE models and VI filtering methods are compared, and the spatial and temporal heterogeneities are further investigated along latitude and season.

3.1 Validation in all sites

The accuracies of the 8-day GPP simulation result using the raw LAI and the four filtering results are evaluated using EC GPP data from 169 sites worldwide. Figure 2 provides the density scatterplot of all data points with different filtering methods. It can be observed that all four filtering methods can improve the accuracy of GPP simulation from both the MOD17 and the RTL-LUE models. When comparing the accuracy within the same model, the utilization of filtered VIs leads to higher R², lower RMSE, and a fitted slope approaching 1, which suggests that four filtering methods can improve the accuracy of LUE-GPP simulations. Due to the deterministic nature of the model, which is hard to conduct significance testing for the R² changes between different models (Casella and Fienberg 2009), we only test the significance of R² for each model. The results indicate all the results passed the significant test with p < 0.05. Among the four filtering methods, SeasonL1 filtering exhibits the most notable enhancement in simulation accuracy (∆R² is 0.06 in MOD17 and 0.07 in RTL-LUE), indicating it is more effective in mitigating GPP underestimation caused by VI data quality. Furthermore, the slope of the MOD17 model is generally lower than that of the RTL-LUE model, indicating a more severe underestimating tendency in the MOD17 model, which has been demonstrated to be due to inadequate consideration of the contribution of shaded leaves (He et al. 2018; Zhang et al. 2012).

Figure 2. The relation between EC GPP and the simulation results from (a)~(e) the MOD17 model and (f)~(j) the RTL-LUE model. Each model has 5 figures, representing (a) raw LAI, (b) SG filtering, (c) Whittaker filtering, (d) HANTS filtering, and (e) SeasonL1 filtering.

3.2 Validation in different vegetation types

3.2.1 Quantitative accuracy comparison

Table 2 and Table 3 present the R² and RMSE between EC GPP and GPP estimations for different vegetation types. It can be observed in Table 2 that the GPP simulation accuracy improves for all vegetation types after using filtering VIs, with a higher R² and lower RMSE compared to the results using raw LAI. Within the MOD17 model, CSH sites show the most improvement in R², reaching 0.171, while the sites with vegetation type EBF experience the greatest reduction in RMSE, with a decrease of 5.405 g C m⁻² 8day⁻¹. There are seven vegetation types that show an R² improvement exceeding 0.05, and six vegetation types that show an RMSE reduction exceeding 1 g C m⁻² 8day⁻¹. Various filtering methods have different performances in different vegetation types. The SeasonL1 method exhibits the highest R² for all the 12 vegetation types with an average value of 0.597, whereas the R² for SG, Whittaker, and HANTS methods are 0.561, 0.577, and 0.581, respectively. Additionally, the RMSE from the SeasonL1 method is also lower in most vegetation types, with a mean value of 16.440 g C m⁻² 8 day⁻¹, while they are 16.579, 16.746 and 16.462 g C m⁻² 8 day⁻¹ for SG, Whittaker, and HANTS methods, respectively. These results suggest that the SeasonL1 method outperforms the other three classical filtering methods in improving the GPP simulation of the MOD17 model in different vegetation types.

Table 2. R² and RMSE of different vegetation types acquired from MOD17.

	R^2					RMSE
Vegetation types^a	Raw LAI	SG	Whittaker	HANTS	SeasonL1	Raw LAI	SG	Whittaker	HANTS	SeasonL1
CRO	0.353	0.316	0.344	0.359	0.375	20.044	21.857	20.801	20.243	19.832
CSH	0.278	0.428	0.442	0.432	0.450	19.430	17.846	17.748	17.596	17.729
DBF	0.659	0.608	0.643	0.658	0.685	25.070	24.163	24.256	23.869	23.640
DNF	0.730	0.576	0.665	0.659	0.734	7.667	7.806	7.575	7.430	7.197
EBF	0.323	0.407	0.406	0.402	0.410	24.280	18.303	18.791	18.574	18.875
ENF	0.663	0.716	0.724	0.723	0.733	18.487	16.603	17.124	16.866	17.005
GRA	0.529	0.532	0.546	0.554	0.574	13.460	12.495	12.443	12.224	12.082
MF	0.641	0.694	0.716	0.715	0.733	26.862	24.685	25.102	24.821	24.873
OSH	0.455	0.467	0.485	0.495	0.501	5.131	4.827	4.920	4.794	4.805
SAV	0.335	0.377	0.392	0.379	0.403	16.060	15.551	15.816	15.512	15.579
WET	0.581	0.610	0.621	0.632	0.646	16.335	15.569	15.941	15.601	15.590
WSA	0.389	0.420	0.436	0.419	0.448	16.441	15.504	15.852	15.683	15.770
ALL^b	0.541	0.561	0.577	0.581	0.597	17.873	16.579	16.746	16.462	16.440

^aThe abbreviations of different vegetation types are introduced in Section 2.1.

^bALL refers to the mean value of R² and RMSE for all 169 sites.

^cThe unit of RMSE is g C m⁻² 8day⁻¹.

Table 3. R² and RMSE of different vegetation types acquired from RTL-LUE.

	R^2					RMSE
Vegetation types	Raw LAI	SG	Whittaker	HANTS	SeasonL1	Raw LAI	SG	Whittaker	HANTS	SeasonL1
CRO	0.541	0.480	0.523	0.555	0.570	28.159	27.464	27.732	26.919	26.765
CSH	0.325	0.585	0.606	0.599	0.624	11.987	8.944	8.757	8.688	8.650
DBF	0.686	0.666	0.703	0.711	0.739	22.883	21.375	21.352	20.767	20.127
DNF	0.779	0.591	0.766	0.717	0.780	7.606	8.426	7.219	7.490	7.097
EBF	0.309	0.426	0.426	0.421	0.424	23.342	16.203	16.888	16.583	15.969
ENF	0.696	0.750	0.762	0.756	0.767	16.797	13.798	14.457	14.398	14.192
GRA	0.615	0.601	0.623	0.644	0.660	14.334	13.531	13.332	13.081	12.810
MF	0.669	0.731	0.750	0.750	0.761	21.868	18.937	18.631	18.397	18.133
OSH	0.606	0.568	0.595	0.640	0.634	5.946	5.874	5.792	5.547	5.531
SAV	0.428	0.450	0.460	0.454	0.469	18.970	17.726	18.199	17.894	17.837
WET	0.662	0.679	0.687	0.712	0.726	14.908	14.011	14.122	13.685	13.482
WSA	0.486	0.463	0.481	0.485	0.509	14.346	13.118	13.441	13.139	13.100
ALL	0.609	0.624	0.644	0.655	0.669	17.606	15.579	15.782	15.424	15.330

The same conclusion can also be found for RTL-LUE models. As shown in Table 3, the largest R² increase within the RTL-LUE model is 0.299 observed in the CSH type, and the largest reduction in RMSE is 7.373 g C m⁻² 8 day⁻¹ observed in the EBF type, indicating the RTL-LUE model exhibits a greater accuracy improvement compared to the MOD17 model. When comparing the results for different filtering methods, the SeasonL1 method also demonstrates higher R² and lower RMSE in almost all the vegetation types, than that of the other three filtering methods. As a result, it indicates that the GPP estimation can be improved in different vegetation types both for big-leaf and two-leaf models after filtering methods enhance the quality of VI data, and the SeasonaL1 method performs the best in almost all vegetation types.

3.2.2 Qualitative interannual variations

To further investigate how filtering methods improve the accuracy of GPP simulation, seasonal variation curves of GPP and their corresponding LAI time series from representative site years are drawn for qualitative assessment, as shown in Figures 3 and 4. Only the results of the filtering method with the best performance, i.e, SeasonL1, are provided to avoid confusion caused by too many data. It can be observed that there are synchronies between the sudden decrease in the raw LAI time series and the underestimation in the corresponding GPP sequence, which is mainly attributed to cloud contamination in the RS-driven data. However, the underestimations of GPP are effectively suppressed when using the reconstructed LAI time series. The results are consistent both for the MOD17 and the RTL-LUE model. The application of the SeasonL1 method mitigates the underestimation of GPP by reconstructing the LAI time series, resulting in GPP simulations being closer to EC GPP because the unreasonable sudden estimations are well corrected. Although the overall tendency of the GPP estimations is similar when using raw LAI or filtered results, the quantitative indexes can be much better for the latter. Furthermore, Figures 3 and 4 reveal that the interannual variations of the RTL-LUE model results exhibit a closer resemblance to those of EC, in contrast to the MOD17 model. It can be observed that there is also a discrepancy between the EC GPP and the model estimations, either using raw LAI or filtered data. This discrepancy could be caused by various sources of uncertainties in the model structure and vegetation index data obtained from satellites, such as observation errors, sensor degeneration, and spatial heterogeneity (Fensholt and Proud 2012; Ruppert and Linstädter 2014).

Figure 3. Intra-annual variation of EC GPP and GPP simulated by the MOD17 model for different vegetation types. Sites in the figure are selected to represent the 12 vegetation types in the IGBP classification system and the detailed location information is provided in table S1. The title of each subplot is named “vegetation types/site year/site name,” DOY means the day of the year. For each subplot, the upper color figure is the variation of GPP, and the under gray figure is the variation of LAI.

Figure 4. Intra-annual variation of EC GPP and GPP simulated by the RTL-LUE model for different vegetation types. The detailed location information is provided in table S1. The title of each subplot is named “vegetation types/site year/site name,” DOY means the day of the year. For each subplot, the upper color figure is the variation of GPP, and the under gray figure is the variation of LAI.

3.3 Spatial and temporal heterogeneity

3.3.1 Variation along latitude

Since the data quality of VI may vary along space and time, its impacts on GPP simulation should also show spatial and temporal heterogeneities. To further explore the spatial patterns of VI reconstruction effects on GPP estimation, all the EC sites are divided into six groups according to their latitudes (i.e. 20 ~ 40°S, 0 ~ 20°S, 0 ~ 20°N, 20 ~ 40°N, 40 ~ 60°S, >60°S), and the improvements of GPP simulation accuracy (ΔR², ΔRMSE) along different latitude zones are shown in Figure 5. The number in each grid denotes the R² and RMSE of GPP simulation using different LAIs, and the color in the right four grids denotes the improvements of using filtered LAI compared with the raw LAI results. It can be observed that the sites located between 20°N and 20°S show the greatest improvement in GPP simulation accuracy after filtering among the six latitude groups, which suggests that the filtering process is more important in low latitude regions around the equator. The color variations also indicate an uneven accuracy change across different latitude zones, demonstrating the spatial heterogeneity of the filtering impacts on the VI data-driven LUE GPP model. Furthermore, as illustrated by the numbers in Figure 5, the simulation accuracy at low latitudes is lower than that at high latitudes in both the MOD17 and RTL-LUE models, which is mainly due to the higher frequency of overcast and cloudy days in low-latitude regions, leading to reduced reliability of satellite VI data. In addition, due to the neglect of diffuse PAR in in the MOD17 model (Yuan et al. 2014), the RTL-LUE model also shows an obvious increase in R² and a decrease in RMSE along all latitude zones (Guan et al. 2021).

Figure 5. Comparison of R² for simulating GPP at different latitudes using the MOD17 model (a) and the RTL-LUE model (b). ∆R² (or ∆RMSE) means the difference between the R² (or RMSE) after SeasonL1 filtering and the R² (or RMSE) without filtering. The unit of RMSE is g C m⁻² 8day⁻¹.

3.3.2 Accuracy differences in different seasons

The quality of VI data also varies in different seasons due to differences in cloudy and rainy weather. As the vegetation above 20°N is considered to exhibit strong seasonality in productivity, the sites in these regions are selected, and the GPP simulation accuracies are compared in the growing and dormant seasons. In general, the growing season for plants in the Northern Hemisphere is roughly considered to be from May to October, and the other six months are the dormant season. As shown in Figure 6, GPP simulations from the two models both exhibit lower R² and higher RMSE during the growing season compared to the dormant season, which can be attributed to the increased frequency of overcast and cloudy days in this period. As a result, more severe cloud cover effects on the VI data are acquired by remote sensing satellites. Furthermore, the filtering process also improves the GPP simulation accuracy more during the growing season, because it can correct more unreasonable slips. Among them, the SeasonL1 method still performs the best in the RTL-LUE model, with an R² improvement of 0.11 and an RMSE reduction of 3.11 g C m⁻² 8 day⁻¹ during the growing season, and the value is 0.07 and 0.85 g C m⁻² 8 day⁻¹ for the dormant season. This may be attributed to the lower quality of the LAI time series in the growing season, which provides greater potential to improve.

Figure 6. Comparison of R² and RMSE during the growing season and dormant season for the two models. ∆R² (or ∆RMSE) means the difference between the R² (or RMSE) after filtering and the R² (or RMSE) without filtering. These sites are all located above 20°N.

4 Discussion

4.1 Impact of NDVI and EVI data quality on LUE-GPP estimation

Except for LAI, remote sensing spectral vegetation indices such as NDVI and EVI are also widely used to calculate the fPAR and drive the big-leaf LUE models (Prince and Goward 1995; Running et al. 1999; Xiao et al. 2004; Yuan et al. 2007), the related equations are shown in Table S3. To further validate whether the obtained results are consistent among different remote sensing indices, the GPP simulation accuracy differences from the big-leaf LUE model using raw and filtered NDVI and EVI are also compared, to enhance the robustness and reliability of results. Due to the NDVI and EVI products being provided with a time interval of 16 days, the 8-day LAI time series are also synthesized into 16-day intervals using the maximum-value composite (MVC) method (Holben 1986). Table 4 presents the GPP simulation accuracy driven by three different remote sensing VI time series, as well as their results using raw data and filtered results. It can be observed that the accuracy using 16-day LAI is higher than that of 8-day LAI in the MOD17 model (Section 3.1), which may be attributed to the MVC method partially mitigating the contamination of remote sensing data (Holben 1986). The GPP simulation accuracies are all improved after using the filtered data for the three VIs, it indicates the importance of the filtering process. However, huge differences are found among the three VIs, with the greatest improvements being found for EVI, and NDVI demonstrating more stable results. Additionally, although the four filtering methods show various performances in three VIs, the SeasonL1 method still exhibits the best improvements. After the process using SeasonL1, the R² for NDVI-driven, EVI-driven, and LAI-driven simulations are increased by 0.03, 0.08, and 0.05, respectively. The results suggest that, in addition to the LAI time series, filtering methods also contribute to similar improvements in the data quality of NDVI and EVI.

Table 4. Statistic of the 16-day GPP simulation from the MOD17 model driven by NDVI and EVI.

	NDVI			EVI		LAI
		R^2	RMSE^a	R^2	RMSE	R^2	RMSE
Raw VI		0.66	39.08	0.59	50.53	0.62	41.60
SG		0.69	37.96	0.66	34.81	0.62	37.83
Whittaker		0.67	38.45	0.66	35.33	0.64	37.74
HANTS		0.68	37.50	0.66	34.49	0.65	37.93
SeasonL1		0.69	36.87	0.67	34.03	0.67	37.53

^aThe unit of RMSE is g C m⁻² 16day⁻¹

4.2 What makes SeasonL1 the most effective?

The results indicate that the process of data filtering can improve the GPP estimation accuracy to a certain extent (∆R² ~0.07), which is even higher than the improvement by changing the model structure from big leaf to two leaf (∆R² ~0.06). It proves that filtering methods can indeed influence the simulation accuracy of the LUE-GPP model by improving the quality of VI data. Among the four filtering methods, SeasonL1 performs the best in improving the accuracy of LUE-GPP simulations, which warrants further investigation. Previous studies have shown that traditional temporal filtering methods do not consider the prior characteristics of the VI time series and perform poorly in cases of continuous time gaps (Atkinson et al. 2012; Chu et al. 2021). However, the SeasonL1 method innovatively introduces the information from adjacent years combined with assistance from temporal neighbors, and the L1-norm is imposed to establish the two corresponding regularization terms and thus can better characterize their statistical distributions (Chu et al. 2022). Figure 7 displays the comparison of interannual NDVI time series between the SeasonL1 method and the SG method. It can be observed that the SG method fails to reconstruct the deep valley during vegetation senescence, while the SeasonL1 method successfully reconstructs the real NDVI trajectory, which is in good agreement with the original time series. The satisfactory and stable performance of the SeasonL1 method in addressing consecutive gaps is mainly because the newly employed inter-annual similarity term introduces valid information from adjacent years. Compared with the other three typical filters, the SeasonL1 filter shows obvious advantages in recovering the time-continuous missing values of vegetation growth inflection points and preventing over-smoothing, which can extend the ability of the temporal methods in advancing the quality of VI time-series data. Thus, using the filtering results from the SeasonL1 method to simulate GPP can certainly obtain higher accuracy than the other methods, with much better agreements in the temporal continuous missing areas.

Figure 7. Comparison of interannual NDVI time series with different filtering methods. The gray elliptical region is the period during which the two methods show great differences.

4.3 Impact of filtering parameter settings on LUE-GPP estimation

Parameters are important for all the VI filtering methods. Current research lacks a discussion on the parameters of filtering methods, and the impact of filtering parameter settings on their effectiveness is still unclear. To further investigate the impact of the filtering parameter on GPP reconstruction results, we selected the SeasonL1 method that shows the best performances in previous results and analyzed the GPP simulation accuracy (R²) changed with the smoothness parameter using the RTL-LUE model. As shown in Figure 8, R² first increases and then decreases with the increasing smoothness parameter, while RMSE is the opposite, which is a common feature for all the smoothing algorithms. The two lines both reach the optimal value around a smoothness parameter of 0.13. The results suggest that the correlation between filter smoothness and GPP reconstruction accuracy is not always positive, and there exists a range where optimal reconstruction accuracy can be achieved. The accuracy will drop beyond this range, which has important implications for the selection and application of filters in subsequent research on terrestrial ecosystem productivity estimation. Therefore, it is necessary to determine the optimal parameter range for GPP simulation.

Figure 8. The sensitivity of GPP simulations to smoothness parameter in the SeasonL1 filter. The unit of RMSE is g C m⁻² 16 day⁻¹.

5 Conclusion

This study has investigated the impacts of the VI data quality on GPP estimation by comparing the GPP simulation results of two LUE models (MOD17 and RTL-LUE), using raw VI data and the results from four data filtering methods (SG, Whittaker, HANTS, and SeasonL1). The main findings are:

• All four filtering methods can improve the accuracy of GPP estimation both using the big-leaf LUE model (∆R² = 0.06) and the two-leaf LUE model (∆R² = 0.07) for most vegetation types. The improvements by VI data filtering are even more significant than the change of model structure from big-leaf to two-leaf.

• The SeasonL1 method outperforms the other three classical filtering methods in improving the GPP simulations using both the two LUE models, because it can effectively recover the missing values and noise in VI data and preserve the key change points of vegetation growth in temporal continuous missing conditions.

• The improvements in GPP simulation accuracy impacted by VI reconstruction vary across different spatial and temporal contexts. The improvements in GPP estimation accuracy are more pronounced for low latitude regions (20 °S ~ 20 °N) and the vegetation growing seasons.

Although our study provides useful insights for future studies related to modeling and forecasting GPP using LUE models, there are also some limitations. Only 169 sites are applied to investigate the impacts of the VI data quality on GPP simulation, which is insufficient to reveal the distribution of these impacts on a global scale. Therefore, it is important to evaluate the quality and uncertainty of VI data and their impacts on GPP estimation at a global scale using multiple approaches and datasets. Furthermore, even concerning previous research parameters, the uncertainty of LUE model parameters may affect the accuracy of the MOD17 model and RTL-LUE model. In addition to the choice of filtering method, the uncertainty in filter smoothing settings can also influence the performance of filtering and, consequently, affect the simulation accuracy of the LUE models. Future research could address these issues and explore more advanced filtering methods and models for better GPP simulation.

Disclosure statement

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

Data availability statement

The authors acknowledge NASA for free access to MODIS VI products and the FLUXNET community for providing the EC data. The VI data are accessed free from the Earth Explorer tool from the United States Geological Survey (USGS) website (https://earthexplorer.usgs.gov). All the flux and meteorological data are available from the website (https://fluxnet.org/data/fluxnet2015-dataset/), which are acquired and shared by global collaborated regional networks supporting the FLUXNET.

Supplementary material

Supplemental data for this article can be accessed online at https://doi.org/10.1080/15481603.2023.2275421.

Additional information

Funding

This research was supported by the National Natural Science Foundation of China (42001371), and the Open Fund of Hubei Luojia Laboratory (220100041).