Compositing, smoothing, and gap-filling techniques

C H A P T E R

3 Compositing, smoothing, and gap-filling techniques O U T L I N E 3.2.3 Temporal spatial filter algorithm 3.2.4 Smoothing and gap-filling algorithm based on the wavelet transform 3.2.5 Time series surface reflectance reconstruction 3.2.5.1 Surface reflectance screening 3.2.5.2 Surface reflectance composition 3.2.5.3 NDVI reconstruction 3.2.5.4 Cloud detection of surface reflectance 3.2.5.5 Surface reflectance reconstruction

3.1 Multitemporal compositing techniques 108 3.1.1 Maximum vegetation index composite 108 3.1.2 Minimum band reflectance composite 108 3.1.3 Maximum surface temperature composite 109 3.1.4 Mixing criteria compositing 109 3.1.5 MODIS vegetation index compositing technique 110 3.2 Time series data smoothing and gap filling 3.2.1 Curve fitting method 3.2.1.1 Adaptive SG filtering 3.2.1.2 Asymmetric Gaussian function and double logistic function fitting 3.2.2 Ecosystem-dependent temporal interpolation technique

111 113 113

114

115

118 118 121 121 121 123 123

3.3 Summary

128

References

128

generate land surface parameter products that are temporally and spatially continuous. A multitemporal compositing method is designed to select the optimal remote sensing data as the pixel value in a composite time window, according to certain criteria. This is a common processing method used to obtain cloud-free and spatially continuous images. Section 3.1 of this chapter gives a detailed description of several commonly used compositing algorithms. Here, time series data smoothing and gap-filling methods are applied to reconstruct high-

Abstract Because of the influences of cloud cover, seasonal snow, and many other factors, time series of land surface parameters extracted from remote sensing data often suffer from discontinuities and missing data, which have seriously restricted the application of extracted land surface parameters to research in global change and other fields. This chapter introduces several different techniques for compositing, smoothing, and gap filling of remotely sensed time series data, which Advanced Remote Sensing, Second Edition https://doi.org/10.1016/B978-0-12-815826-5.00003-9

117

107

© 2020 Elsevier Inc. All rights reserved.

108

3. Compositing, smoothing, and gap-filling techniques

quality time series data with temporal continuity and spatial integrity, thereby eliminating the influences of cloud. Section 3.2 introduces several typical smoothing and gap-filling algorithms for time series data.

3.1 Multitemporal compositing techniques Because of the influences of clouds, aerosols, and other factors, single-phase remote sensing images are often subject to missing data, seriously limiting the application of extracted land surface parameters to research in global change and many other fields. Through multitemporal remote sensing data compositing, the influences of these factors can be reduced or eliminated (Holben, 1986; Cihlar et al., 1994). Multitemporal compositing uses set criteria to select the pixel values with the highest quality among multiscene matching data within a specific time range. For different application targets, a variety of data compositing criteria have been proposed, each with the purpose of eliminating the influences of cloud cover, aerosols, and other factors as much as possible (Cihlar et al., 1994; Qi and Kerr, 1997). Several typical data compositing algorithms are briefly introduced in this section.

3.1.1 Maximum vegetation index composite Vegetation index is usually obtained by linear or nonlinear combination operations on remotely sensed red and near-infrared (NIR) reflectance data, which are simple and effective parameters for characterizing vegetation cover and growth status. The maximum vegetation index composite chooses the maximum value of the vegetation index within a given time range as the selection criterion for the remote sensing data. The maximum value composite (MVC), first proposed by Holben (1986), was applied to the maximum normalized difference vegetation index (NDVI) to synthesize multitemporal AVHRR data. Based on the multiscene matching satellite observations within a given time range, the corresponding NDVI images are calculated,

and then a comparison is performed pixel by pixel to choose the remote sensing data with maximum NDVI. This method is simple and easily implemented, which is why it has been widely applied to monitor changes in global vegetation cover (Illera et al., 1996; Kasischke et al., 1993; Peters et al., 2002; Potter and Brooks, 2000). However, MVC tends to select data that are far away from the subsatellite point in practical applications (Cihlar et al., 1994), yielding poor performance when eliminating cloud cover for certain vegetation types (Sousa et al., 2002).

3.1.2 Minimum band reflectance composite To avoid some of the problems associated with MVC, selection based on the minimum band reflectance was proposed by Qi and Kerr (1997). Generally, these methods involve selecting the minimum blue band reflectance (mBlue) (Vermote and Vermeulen, 1999) and minimum red band reflectance (mRed) (Cabral et al., 2003; Chuvieco et al., 2005). In the red band, the reflectance of cloud is significantly higher than that of vegetation. If the minimum reflectance of the red channel is selected, the probability of selecting cloudcontaminated pixel values should be reduced. However, this method still retains a large number of cloud-contaminated pixel values because the reflectance of the red channel in cloud-shadow areas is also very low (Qi and Kerr, 1997). To eliminate the influences of cloud shadows, selection of the third lowest near-infrared reflectance (NIRm3) and the third darkest value (Darkm3) has been suggested by some scholars. Selection of the third lowest values is based on the assumption that there is a low likelihood of a cloud shadow falling on a given pixel more than twice, over a period of 1 month (Cabral et al., 2003). The composite results of daily VEGETATION data obtained by four compositing techniques are shown in Fig. 3.1. All algorithms performed similarly with respect to cloud elimination. However, mRed retained cloud shadows, and

3.1 Multitemporal compositing techniques

109

FIGURE 3.1 Comparison of different compositing techniques: (1) maximum value composite, (2) mRed, (3) NIRm3, and (4) Darkm3 (Cabral et al., 2003).

MVC yielded highly heterogeneous images because of the different acquisition geometries of neighboring pixels (Cabral et al., 2003). Comparatively speaking, the NIRm3 and Darkm3 were able to generate high-quality images which appeared to be visually smoother, while effectively eliminating the influences of clouds and cloud shadows.

3.1.3 Maximum surface temperature composite For burned land mapping, the maximum surface temperature (MaxTs) is an accurate criterion to composite remote sensing data. This is because the temperature in recently burned areas is higher than that of the surrounding unaffected vegetation and is also higher than those of clouds and cloud shadows (Pereira et al., 1999).

In composites of MODIS and AVHRR data, the MaxTs criterion has achieved good results in distinguishing burned areas from the surrounding vegetation. Fig. 3.2 shows the MODIS multitemporal composites generated by MVC and MaxTs over an area affected by a large forest fire. By comparing the results of the two methods, it can be seen that the MaxTs composite can distinguish burned areas from the surrounding vegetation more clearly than the MVC method.

3.1.4 Mixing criteria compositing In multitemporal compositing, the use of mixing criteria to select remote sensing data was proposed by Qi et al. (1993). The above criteria for choosing maximum vegetation index, maximum surface temperature, and maximum view zenith angle can be combined to form a

110

3. Compositing, smoothing, and gap-filling techniques

FIGURE 3.2 Composite images of MODIS multitemporal data over a burned area generated with different criteria: (1) maximum value composite and (2) MaxTs (Chuvieco et al., 2005).

variety of data compositing techniques (Cabral et al., 2003; Carreiras and Pereira, 2005; Carreiras et al., 2003; Cihlar et al., 1994). Fig. 3.3 shows the results of applying several criteria for multitemporal compositing to MODIS multitemporal images collected in September 2003. Visually, all the results appear to be significantly better than the daily observation data, but clouds and cloud shadows were still present in the composite results for MaxNDVI, MinR1, and MinR2. This means that MaxNDVI, MinR1, and MinR2 were unable to properly eliminate the cloud-contaminated data. The MaxTs method yielded the best composite results, with nearly all clouds and cloud shadows eliminated from the daily remote sensing data. The compositing techniques that combine criteria for choosing maximum temperature and minimum view zenith angle or minimum NIR reflectance have been demonstrated to be capable of eliminating all clouds and cloud shadows.

3.1.5 MODIS vegetation index compositing technique In the MVC approach, data far away from the subsatellite point are more likely to be selected than those closer by. Therefore, the selection of remote sensing data with a minimum view zenith angle has been proposed (Chuvieco

et al., 2005). However, before data compositing, strict data screening is required to eliminate cloud-contaminated data. To avoid the selection of data far away from the subsatellite point, a bidirectional reflectance distribution function (BRDF) model correction is adopted to synthesize MODIS vegetation indices (Van Leeuwen et al., 1999; Schaaf et al., 2002). Fig. 3.4 shows the flowchart for the MODIS vegetation index compositing algorithm. Based on cloud identification and other information, high-quality directional reflectance data are selected within a 16-day time window. When the number of high-quality directional reflectance data within the time window is greater than 5, the Walthall BRDF model is used for fitting the reflectance of each band. The Walthall BRDF model can be expressed by Eq. (3.1) (Walthall et al., 1985) as follows: rl ðqv ; fs ; fv Þ ¼ al q2v þ bl qv cosðfv  fs Þ þ cl (3.1) where al , bl , and cl are model parameters determined by least squares fitting. The nadir reflectance for each pixel within a 16-day interval can then be derived. When the number of cloud-free data during a 16-day period is less than 5, the constrained view angle, maximum value composite (CV-MVC) approach is adopted to select the highest NDVI

3.2 Time series data smoothing and gap filling

111

FIGURE 3.3 Composite results for MODIS multitemporal images in September 2003. (1) MinR1, (2) MinR2, (3) MaxTs, (4) MaxNDVI, (5) MinAZMaxTs, (6) MinR2MaxTs, and (7) MaxTsMinR2 (Chuvieco et al., 2005).

based on two cloud-free pixels with their view angles closest to nadir. First, all cloud-free reflectance data are sorted within the time window by view zenith angle, and then two cloud-free reflectance data with view angles closest to the view zenith angle are selected to calculate NDVI. The reflectance data with maximum NDVI values are selected as the composite results within the time window. If only one cloud-free observation exists within the time window, this observation will be automatically selected to represent the optimal value within the composite time. If 16-day observation data are all subject to cloud contamination, then the vegetation indices are calculated for all days, and the pixel value with the maximum NDVI among all the observations is selected.

The MODIS vegetation index compositing algorithm was applied to 1 year of daily AVHRR data to generate 23 consecutive global composites of the NDVI. Fig. 3.5 shows the results of the NDVI composite from August 13 to 28, 1989. It can be seen from Fig. 3.6 that, over the entire time profile, NDVI values obtained by the MODIS vegetation index compositing technique are slightly smaller than those obtained by MVC.

3.2 Time series data smoothing and gap filling If the composite time period is too long, the time compositing technology is unable to reflect

112

3. Compositing, smoothing, and gap-filling techniques

FIGURE 3.4 Flowchart for the MODIS vegetation index compositing algorithm applied to AVHRR test datasets (Van Leeuwen et al., 1999).

FIGURE 3.5 Global normalized difference vegetation index image using the MODIS vegetation index compositing algorithm (Van Leeuwen et al., 1999).

real changes in land surface parameters; if it is too short, the influences of clouds cannot be effectively eliminated, especially in cloudy areas. The smoothing and gap-filling method for time series data aims to reconstruct high-quality time series data and eliminate the influences of clouds. To date, a variety of smoothing and gap-filling methods for time series of remote sensing data have been developed (Viovy et al.,

1992; Hermance, 2007; Julien et al., 2006; Roerink et al., 2000; Moody et al., 2005). Generally speaking, these methods can be divided into two categories: the first is smoothing and gap filling in the time domain, mainly by asymmetric Gaussian function fitting (Jonsson and Eklundh, 2002), weighted least squares linear regression (Sellers et al., 1994), or a SavitzkyeGolay (SG) filter (Chen et al., 2004; Cao et al., 2018); the other is

3.2 Time series data smoothing and gap filling

BRDF/CV-MVC/MVC MODIS (NDVI)

0.3

NDVI

35

MVC (NDVI)

30

MODIS (RD %)

0.25

25

0.2

20

0.15

15

0.1

10

0.05

5

0

1

65 33

129 97

193

161

257

225

321

289

RD (%)

0.35

0

353

DOY

FIGURE 3.6 Comparison between global means of normalized difference vegetation index obtained by MVC and the MODIS vegetation index compositing algorithm (Van Leeuwen et al., 1999).

smoothing and gap filling in the frequency domain, e.g., the fitting method with a Fourier function (Roerink et al., 2000). This section will focus on several common smoothing and gap-filling methods for time series data.

113

Note that the atmospheric effects on surface reflectances generally cause increases in reflectance in the red band rather than in the NIR band, resulting in decreases in the retrieved LAI (Fernandes et al., 2003; Chen et al., 2006). To account for negatively biased noise, the fitting is implemented in two steps (J€ onsson and Eklundh, 2002). The first step minimizes the cost function when all data points have the same weights. It is generally believed that data points with values greater than those simulated by the fitting model in the first step are less affected by noise and are thus of higher quality. Therefore, in the second step, the minimization is repeated with the weight of low-quality data points decreased by a factor of 10. This yields the smooth envelope curve for the NDVI time series, as shown in Fig. 3.7. In Fig. 3.7, thin solid lines indicate the NDVI time series; dotted lines indicate the fitting curves after the first iteration; and thick solid lines indicate the NDVI envelope curve after the two-step iteration. 3.2.1.1 Adaptive SG filtering To suppress noise and to smooth data, a polynomial equation f ðtÞ ¼ c1 þ c2 t þ c3 t2 is used for fitting 2n þ 1 data within a sliding

3.2.1 Curve fitting method Based on least squares fitting, multiple fitting methods for the outer envelope of NDVI time series data were proposed by J€ onsson and Eklundh (2004). These methods mainly include adaptive SG filtering, an asymmetric Gaussian function and double logistic function fitting. Given time series data ðti ; yi Þ;i ¼ 1; 2; .; N, the parameter c in the fitting model f ðt; cÞ can be obtained through minimizing the following cost function: c2 ¼

N X

½ui ðf ðti ; cÞ  yi Þ2

(3.2)

i¼1

where ui is the weight of the ith data point.

FIGURE 3.7 The upper envelope of normalized difference vegetation index time series data obtained by two-step fitting (J€ onsson and Eklundh, 2004).

114

3. Compositing, smoothing, and gap-filling techniques

window at each position i of the time series data ðti ; yi Þ;i ¼ 1; 2; .; N. Data are equally spaced in SG filtering (Savitzky, and Golay, 1964). Then, after fitting the polynomial equation, yi is replaced by the simulated value gi at position i. Considering the negative bias in remote sensing data noise, the above two steps are undertaken to obtain the smooth envelope curve of time series data. The sliding window size n not only determines the smoothness of the results of adaptive SG filtering but also greatly influences the ability to fit rapidly changing time series data by an adaptive SG filter. Fig. 3.8 shows the filtering results obtained with sliding windows of five and three points, respectively. It can be seen that with sliding window of 5, some anomalous peak values are not fitted properly; meanwhile, with the sliding window of 3, the filtering results closely follow the anomalous peak values. 3.2.1.2 Asymmetric Gaussian function and double logistic function fitting

(3.4)

where x1 is the position at which the basis function reaches its maximum or minimum; x2 and x3 are the width and flatness of the right basis function, respectively; and similarly x4 and x5 are the width and flatness of the left basis function, respectively. To ensure the smooth results of model simulation corresponding to time series data, parameters x2 , x3 , x4 , and x5 should lie within specified ranges. In double logistic function fitting, the basis function is as follows: gðt; x1 ; x2 ; x3 ; x4 Þ ¼

In both of the function fitting methods, the fitting model usually has the following form: f ðt; cÞ ¼ c1 þ c2 gðt; xÞ

In asymmetric Gaussian function fitting, the following basis function is selected:

(3.3)

where x is the parameter of the basis function gðt; xÞ.



1   x1  t 1 þ exp x2

1   x3  t 1 þ exp x4

(3.5)

FIGURE 3.8 Results of adaptive SG filtering with different sliding window lengths (J€onsson and Eklundh, 2004). (A) the sliding window of 5 and (B) the sliding window of 3.

3.2 Time series data smoothing and gap filling

where x1 is the position of the left inflection point and x2 is the slope of the left inflection point of the basis function; similarly, x3 is the position of the right inflection point and x4 is the slope of the right inflection point. As in the asymmetric Gaussian function fitting, specified value ranges are also required for these parameters to ensure smooth results in the model simulation. MODIS LAI time series data were used for comparing the above fitting methods by Gao et al. (2008) (Fig. 3.9). The asymmetric Gaussian function fitting was applied to smoothing and gap filling of the MODIS LAI product to generate high-quality LAI products with temporal and spatial continuity.

3.2.2 Ecosystem-dependent temporal interpolation technique Within a small region, pixels of the same ecological category show similar growth and phenological behaviors, which is the basis for such categorization in describing how vegetation parameters of an ecological category change over time. Each ecological category enters consecutive phenological stages (green-up, full maturity, senescence, and dormancy) at approximately the same time. Therefore, the parameters

FIGURE 3.9

115

at pixel level in the region will represent very similar ecological behaviors. However, under different growth conditions, discrepancies can be found between pixels in terms of the magnitude of ecological behavior variations. Based on these inherent properties, Moody et al. (2005) proposed an ecosystem-dependent temporal interpolation technique. This method mainly uses pixel-level and regional ecosystemedependent phenological behavioral curves to obtain high-quality time series data of pixel parameters. Using this algorithm, Moody et al. (2005) interpolated the missing data in the standard MOD43B3 albedo product. Fig. 3.10 shows the MOD43B3 whitesky albedo product over four seasons. Fig. 3.11 shows the albedo product after interpolation. The albedo product obtained using ecosystemdependent temporal interpolation is not only able to maintain the same quality as that of the original data but also yields the statistical properties of the filled data and an estimate of the quality of the processed data. Ding et al. (2017) improved the ecosystemdependent temporal interpolation technique to address the uncertainties in interpolating missing data in satellite-derived LAI time series for heterogeneous landscape, particularly

Comparison of results obtained by several fitting methods (Gao et al., 2008).

116

3. Compositing, smoothing, and gap-filling techniques

FIGURE 3.10

Global MOD43B3 white-sky albedo for four 16-day periods: (A) January 1e16, (B) April 3e18, (C) July 12e27, and (D) September 30eOctober 14, 2002 (Moody et al., 2005).

FIGURE 3.11 Spatially complete white-sky albedo after the temporal interpolation technique was applied to four 16-day periods: (A) January 1e16, (B) April 3e18, (C) July 12e27, and (D) September 30eOctober 14, 2002 (Moody et al., 2005).

117

3.2 Time series data smoothing and gap filling

heterogeneous grasslands. The improved algorithm can more accurately predict missing LAI value for different seasons and proportions of missing data (Ding et al., 2017).

3.2.3 Temporal spatial filter algorithm The MODIS standard LAI product is not continuous in time and space. To fill gaps and improve product quality, Fang et al. (2008) proposed a temporal spatial filter (TSF) algorithm to integrate both the spatial and temporal characteristics of different plant functional types and generate continuous LAI data. The purpose of the TSF algorithm is to process pixels with low-quality and missing data in the original products (quality control identification (QC) > 32). The algorithm flowchart is shown in Fig. 3.12 and comprises three main steps: (1) Calculate the background value. For vegetation types with the same function, the pixel with missing data is first supplemented by a multiyear mean. If the pixel has no value for all years, the pixel does not have a multiyear mean and instead the pixel’s background value is calculated by vegetation continuous fieldeecological curve fitting (VCF-ECF). The new VCF-ECF emphasizes the dependence of time series data for each target pixel on regional VCFdependent phenological behaviors, thereby maintaining the temporal and spatial integrity at pixel level. (2) Calculate the observed value. The original MODIS LAI data are treated as observed values. When 32  QC < 128 (LAI is obtained by an alternate algorithm), the original MODIS LAI value can be considered as an observed value. The same quality identification is adopted to calculate the LAI variance of vegetation types with the same function. When QC > 128, there will be no retrieved value of LAI, and local time filtering is needed to acquire the observed

FIGURE 3.12 Temporal spatial filter flowchart. Sourced from Fang, H., Liang, S., Townshend, J.R., Dickinson, R.E., 2008. Spatially and temporally continuous LAI data sets based on and integrated filtering method: examples from North America. Remote Sens. Environ. 112, 75e93.

value of LAI in an adjacent time. If no desired QC value can be found in an adjacent time, then filtering, such as SG filtering (see Section 3.2.1.1), is used to calculate the observed value from the background values. (3) Time-space filtering. The target value xa is calculated using the background value xb and the observed value xo is obtained by the above two steps as follows:

xa ðri Þ ¼

2 E2 b xb ðri Þ þ Eo wðri ; rj Þfxb ðri Þ þ ½xo ðrj Þ  xb ðrj Þg 2 E2 b þ Eo

(3.6)

118

3. Compositing, smoothing, and gap-filling techniques

where E2b and E2o are the error variances of the background and observed values. Assuming that errors in the background and observed values are homogeneous and not spatially related, E2b and E2o are considered as locally independent. wðri ; rj Þ is the weighting function, which is dependent on the time distance di;j between point ri and point rj : ! R2  d2i;j wðri ; rj Þ ¼ max 0; 2 (3.7) R þ d2i;j where R is a prescribed influence radius. When processing MODIS LAI data, Fang et al. (2008) set the influence radius as 16. Fang et al. (2008) employed the TSF algorithm to process the MODIS LAI product and generated a continuous LAI product for North America. Fig. 3.13A shows the MODIS LAI product for North America on the 225th day of 2000. Because of cloud cover, especially at high latitudes, the MODIS LAI product has a lower quality. Processed LAI data are shown in Fig. 3.13B. It can be seen that the LAI product has improved spatial integrity after time-space filtering. These results confirm that the temporal and spatial continuity of the LAI product after TSF filtering is significantly improved. The TSF algorithm can also be applied when processing other land surface parameters. Fang et al. (2007) used the TSF method when filtering the MODIS albedo product to generate a 1-km resolution albedo product with spatial continuity in North America.

3.2.4 Smoothing and gap-filling algorithm based on the wavelet transform The wavelet transform (WT) can be used to analyze signals in timeefrequency space and reduce noise, while retaining the important components in the original signals. In the past 20 years, WT has become a very effective tool in signal processing. Currently, WT is widely used in numerous remote sensing applications,

such as image fusion (see Chapter 4, Ranchina et al., 2003), the extraction of hyperspectral data characteristics (Pu and Gong, 2004), and phenological detection (Sakamoto et al., 2005). Lu et al. (2007) proposed an algorithm for applying the WT to eliminate noise in time series observations and to fill missing data. First, based on quality identification and blue band data, linear interpolation is performed on the time series data; then, the time series data are decomposed into different time scales, and a set of adjacent time scales with the highest correlation is used to construct the new time series data. Figs. 3.14 and 3.15 show time series data for the MODIS LAI and albedo products before and after WT, respectively. Compared with the original data, the reconstructed LAI and albedo time series data are relatively smooth, have no major fluctuations, and are more consistent with the variations of these parameters in nature.

3.2.5 Time series surface reflectance reconstruction Currently, several datasets based on AVHRR data were produced using different data processing methods (James and Kalluri, 1994; Los et al., 2000; Pedelty et al., 2007; Tucker et al., 2005; Pinzon and Tucker, 2014). These datasets, with various spatial and temporal resolutions, have different temporal extent. Among them, the representative dataset is now being routinely generated by the Land Long-Term Data Record (LTDR) project, the goal of which is to produce a global land surface climate data record (CDR) from the AVHRR, the Moderate Resolution Imaging Spectroradiometer (MODIS), and the Visible Infrared Imaging Radiometer Suite (VIIRS) instruments using a MODIS-like operational production approach (Pedelty et al., 2007). The LTDR project has created daily surface reflectance and NDVI products from AVHRR observations.

3.2 Time series data smoothing and gap filling

119

Comparison of MODIS LAI products before and after temporal spatial filter (TSF) filtering: (A) MODIS LAI (Day 225, 2000) and (B) LAI map after TSF filtering.

FIGURE 3.13

FIGURE 3.14

Comparison of MODIS LAI products before and after wavelet transform (Lu et al., 2007).

120

3. Compositing, smoothing, and gap-filling techniques

FIGURE 3.15

Comparison of MODIS albedo products before and after wavelet transform (Lu et al., 2007).

Although the quality of the LTDR AVHRR products was greatly improved over earlier versions (Nagol et al., 2014; Sobrino and Julien, 2016; Tian et al., 2015), some cloudy or partially cloudy pixels remained in the surface reflectance and NDVI products, and on some occasions, several days of surface reflectance and NDVI values were missing from the time series. Residual cloud and aerosol contamination in the surface reflectance and NDVI products significantly limits their applications of land surface monitoring and results in temporal and spatial inconsistencies in subsequent downstream products. Hence, it has been widely recognized that consistent and gap-filled land surface reflectance and NDVI time series are fundamentally important. As mentioned above, a variety of methods were developed to reconstruct temporal NDVI trajectories (Chen et al., 2004; Hermance, 2007;

Julien et al., 2006; Lu et al., 2007; Viovy et al., 1992; Xiao et al., 2015). However, there are only a few ways to reconstruct surface reflectance time series. Tang et al. (2013) developed a time series cloud detection (TSCD) algorithm to remove cloud contamination in the MODIS surface reflectance product (MOD09A1) and fill in missing pixels. The TSCD algorithm performs very well, particularly when the land surface is stable or changing only slowly (Tang et al., 2013). However, the TSCD algorithm is dependent on surface reflectance in the blue band and other auxiliary information, which makes it difficult to reprocess the AVHRR surface reflectance data. In view of the shortcomings of currently available NDVI and surface reflectance reconstruction algorithms, Xiao et al. (2015) developed a temporally continuous vegetation indicesebased land surface reflectance reconstruction (VIRR) method

3.2 Time series data smoothing and gap filling

to reconstruct time series of surface reflectance. The daily LTDR AVHRR surface reflectance data were first aggregated into 8-day intervals to maintain a temporal resolution consistent with the MODIS and VIIRS surface reflectance data and the aggregated 8-day surface reflectance data were used to calculate NDVI. Compared with surface reflectance, NDVI exhibits obvious seasonal changes that make it relatively easy to reconstruct temporally continuous NDVI. Considering the negative bias of NDVI when surface reflectance values are contaminated by clouds, NDVI upper envelopes were constructed. The cloud-contaminated surface reflectance values were identified using the NDVI time series and its upper envelope. Then the surface reflectance time series was further reconstructed from cloud-free surface reflectance values by incorporating the NDVI upper envelopes as constraints. Fig. 3.16 shows the general process of the VIRR method to reconstruct temporally continuous NDVI and surface reflectance from LTDR AVHRR data. Detailed descriptions of the important parts of the VIRR algorithm are given below. 3.2.5.1 Surface reflectance screening An initial screening process was applied to the daily LTDR AVHRR surface reflectance data to determine whether the surface reflectance values were valid. The characteristic feature of vegetated target spectra is minimum reflectance in the red band due to strong absorption of photosynthetically active radiation by chlorophyll in green leaves and maximum reflectance in the NIR band because of strong reflection by leaf cells. However, the NIR reflectance in cloud-contaminated or snow-covered pixels is lower than reflectance in the red band (Luo et al., 2008). Generally, cloud contamination and snow cover are considered as noise for retrieval of land surface parameters such as LAI and FAPAR. Therefore, for each data point, the surface reflectance values were assumed to be invalid if the reflectance in the red band was

121

higher than the NIR reflectance. At the same time, the two-band enhanced vegetation index (EVI2) should be lower than the NDVI under normal conditions (Jiang et al., 2008; Zhang, 2015). If the EVI2 values were higher than the NDVI values, the surface reflectance values were also assumed to be invalid. In this study, only valid surface reflectance values were used for compositing. 3.2.5.2 Surface reflectance composition The daily surface reflectance data were formed into 8-day composite images from 3 days before through 4 days after each composite time. If the number of data points with valid surface reflectance values over an 8-day compositing period was two or more, the CV-MVC method was used. The view zenith angles were taken for all data points with valid surface reflectance values over the 8-day compositing period, after which these view zenith angles were sorted in ascending order. Then, the NDVI values were computed for the data points with the two smallest view zenith angles, after which the surface reflectance with the higher NDVI value was selected (Van Leeuwen et al., 1999). If only one valid surface reflectance value was available over the composite period, this surface reflectance value was automatically selected to represent the period. If all surface reflectance values during the 8-day period were invalid, a multiyear mean was calculated to fill the gap. 3.2.5.3 NDVI reconstruction Composite surface reflectance data are still susceptible to spurious data points. Therefore, the NDVI values derived from surface reflectance are also temporally discontinuous. Considering the negative bias of NDVI when surface reflectance values are contaminated by clouds, a penalized least squares regression based on a three-dimensional discrete cosine transform (DCT-PLS) developed by Garcia (2010) was used to reconstruct the NDVI upper envelopes.

122

3. Compositing, smoothing, and gap-filling techniques

LTDR AVHRR daily surface reflectance in red and NIR bands

Constrained-view angle,

n=1 n=0

Number of good observations in 8 daily images for each pixel (n)

n>1

maximum value composite Aggregated 8-day surface reflectance in red and NIR bands

Calculating multi-year mean

NDVI reconstruction

Smooth and continuous upper envelopes of NDVI

NDVI calculation

Cloud dection

Reconstructed 8-day surface reflectance in red and NIR bands

High quality surface

Fitting models to

reflectance data

surface reflectance data Y

Simulated NDVI

Cost function

Minimum reached? N

NDVI calculation

Simulated time series surface reflectance data

Adjust control variables

Fitting models to surface reflectance data

FIGURE 3.16 Schematic description of the method used to reconstruct normalized difference vegetation index and surface reflectance time series from LTDR AVHRR.

3.2 Time series data smoothing and gap filling

Let X be a vector that contains a series of equally spaced NDVI values xi ; i ¼ 1; 2; /; n. Let W be the diagonal matrix that contains the weights, wi, corresponding to the NDVI value xi. The DCT-PLS consists in minimizing the following functional to find the best smoothing b of X (Garcia, 2010). estimate, X,     2   b 2 b ¼ b X  (3.8) F X W1=2 X  þ sD X where D is the Laplace operator, and s is a real positive scalar that controls the degree of b can be achieved smoothing. Minimizing on X using an iterative procedure. With type-2 DCT and inverse discrete cosine transform (IDCT), b can be expressed as the X b fkþ1g X

     b fkg þ X b fkg ¼ IDCT G DCT W X  X (3.9)

b fkg refers to X b calculated at the kth iterawhere X tion step, and G is a diagonal matrix with the components  1 Gi;i ¼ 1 þ sð2  2 cosðði  1Þp=nÞÞ2 (3.10) Note that atmospheric effects on surface reflectance generally cause an increase in reflectance in the red band rather than in the NIR band, resulting in decreases in the NDVI. To account for negatively biased noise, we used the weight function 8 1:0 ui > 0:0 > > > > 2  >   > < ui 2 ui < 0:0 1.0 < 1:0  wi ¼ 4:685 4:685 > > > > ui > > : < 1:0 0:0 4:685 (3.11)

123

where ui is the studentized residual, calculated by sffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 11 1 þ 1 þ 16 s ui ¼ ri @1:4826 MADðrÞ 1  pffiffiffipffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi A 2 1 þ 16 s 0

(3.12)

where ri ¼ xi  b x i is the residual of the ith observation, and MAD denotes the median absolute deviation (Rousseeuw and Croux, 1993). This iteratively reweighted procedure leads to a smoothed curve adapted to the upper envelope of the NDVI in the time series. 3.2.5.4 Cloud detection of surface reflectance Atmospheric effects on surface reflectance generally cause negatively biased noise within the vegetation index values. The smoothing method described above was first used to calculate the upper envelopes of NDVI, denoted by NDVI_Env. If NDVI values at each time satisfy the following conditions, the surface reflectance data were deemed to be contaminated by clouds or other factors. Otherwise, they were considered to be cloud free and of high quality. jNDVI  NDVI Envj > a  NDVI Env (3.13) where a is a threshold and was set to 0.4 in this study. The method was applied to check the time series surface reflectance, and any data points assessed as contaminated were filtered out. 3.2.5.5 Surface reflectance reconstruction The surface reflectance time series was reconstructed according to cloud-free surface reflectance values by incorporating the reconstructed NDVI upper envelopes as constraints. Let ðti ; Ii Þ; i ¼ 1; 2; /; m be a series of data points, where ti is time and Ii is the surface

124

3. Compositing, smoothing, and gap-filling techniques

reflectance in the various bands. In this study, only surface reflectance in the red and NIR bands NIR T were reconstructed. Therefore, Ii ¼ rred , i ; ri NIR and r are the surface reflectance where rred i i for the ith position in the red and NIR bands, respectively. For each data point Ii , a quadratic polynomial function, f ðtÞ ¼ at2 þ bt þ c, was fitted to all 2n þ 1 surface reflectance values in each band in a moving time window using the least squares method. An optimization method searched for the trajectories of the quadratic polynomial functions that best fitted the surface reflectance in the given time window, i.e., minimizing n  X

JðXi Þ ¼

j ¼ n

þ

n X j ¼ n

þ

n X

fred ðtiþj Þ  rred iþj

2

 2 fNIR ðtiþj Þ  rNIR iþj ðNDVI Simiþj  NDVI Enviþj Þ2

j ¼ n

(3.14) where fred ðti Þ and fNIR ðti Þ are the quadratic polynomial functions for the red and NIR bands, T respectively; Xi ¼ ðared ; bred ; cred ; aNIR ; bNIR ; cNIR Þ is a control vector consisting of the coefficients of the quadratic polynomial functions; NDVI Simi is a simulated NDVI calculated using the values of the quadratic polynomial functions; and n is the size of the moving time window. In this study, n was set to 2. Then the reconstructed surface reflectance values at the ith position were set to the values of these quadratic polynomial functions according to the optimal values of the control variables. Detailed information about the reconstruction method has been given by Xiao et al. (2015). Xiao et al. (2017) used the VIRR method to generate GLASS AVHRR NDVI and surface reflectance products from 1982 to 2015. To illustrate the performance of the VIRR algorithm, time series of the GLASS AVHRR NDVI and surface reflectance products were compared with

those of the LTDR AVHRR NDVI and surface reflectance products over a sample of sites with different biome types, and the GLASS AVHRR NDVI product was also compared with the reconstructed NDVI from LTDR AVHRR data using the SG filter (denoted by SG AVHRR) over these sites. In addition, the GLASS AVHRR surface reflectance product was compared with the LTDR AVHRR surface reflectance product in space. Fig. 3.17 shows time series of the GLASS AVHRR NDVI and surface reflectance, the LTDR AVHRR NDVI and surface reflectance over a sample of sites with different biome types for 1982e2015. And time series of the SG AVHRR NDVI are also shown in Fig. 3.17 for comparison to evaluate the performance of the refined VIRR algorithm. There were gaps in the LTDR AVHRR surface reflectance in late 1994 and 2000. Fig. 3.17 demonstrates that the refined VIRR algorithm provides reasonable NDVI and surface reflectance values for the missing pixels over these sites, and the SG AVHRR NDVI also shows complete temporal profiles. Over these sites, most of the LTDR AVHRR surface reflectance values in the red band are greater than 0.2, and the corresponding NDVI values are around zero, which may result from serious contamination due to residual clouds. The refined VIRR algorithm successfully removed surface reflectance values contaminated by residual clouds and reconstructed temporally continuous time series of NDVI and land surface reflectance (Fig. 3.17). Excellent agreement was achieved between the time series of GLASS AVHRR NDVI and the upper envelope of the time series NDVI calculated from the LTDR AVHRR surface reflectance product, and the time series of GLASS AVHRR surface reflectance in the red and NIR bands were in good agreement with the lower envelopes of the LTDR AVHRR surface reflectance over these sites. The GLASS AVHRR and SG AVHRR NDVI profiles demonstrated excellent agreement for these years except over the Yucheng site.

3.2 Time series data smoothing and gap filling

FIGURE 3.17 Time series of normalized difference vegetation index (NDVI) and surface reflectance in the red and near-infrared bands from GLASS AVHRR and LTDR AVHRR for the (A) Zhangbei, (B) Yucheng, (C) Puechabon, (D) Larose, (E) Wankama, and (F) Turco sites from 1982 to 2015. For comparison, the reconstructed NDVI profiles using the SG filter over these sites are also shown.

125

cont’d.

3. Compositing, smoothing, and gap-filling techniques

FIGURE 3.17

126

3.2 Time series data smoothing and gap filling

The biome type for Yucheng site is cropland with double annual vegetation seasons. The SG AVHRR NDVI profile over this site was not able to follow sudden decrease of underlying data values between two growing seasons for some years. Examples of the reconstructed surface reflectance are shown in Fig. 3.18. Fig. 3.18A shows an RGB image of the GLASS AVHRR surface reflectance on January 9, 2010, which was composited from the LDTR AVHRR daily surface reflectance from January 6 to 13, 2010, and then reconstructed using the refined VIRR method. For comparison, an RGB image of the LDTR AVHRR surface reflectance on January 9, 2010, is shown in Fig. 3.18B. The images were produced with an RGB color scheme that uses the red band (as red color), the NIR band (as green color), and the red band (as blue color), and the same color enhancement was used. In

(A)

GLASS AVHRR, January 9, 2010

(C)

GLASS AVHRR, July 12, 2010

127

Fig. 3.18B, the LDTR AVHRR surface reflectance is contaminated with residual clouds, and large clouds can in fact be seen over the image. The RGB image in Fig. 3.18A shows that the refined VIRR method effectively removed the residual clouds in the LDTR AVHRR surface reflectance. Fig. 3.18C shows an RGB image of the GLASS AVHRR surface reflectance on July 12, 2010, and an RGB image of the LDTR AVHRR surface reflectance on July 12, 2010, is shown in Fig. 3.18D. The images were produced with the same RGB color scheme and color enhancement as in Fig. 3.18A. In Fig. 3.18D, most areas of the image were covered by residual clouds. The refined VIRR algorithm effectively identified these residual clouds based on the temporally continuous NDVI. The RGB image of the GLASS AVHRR surface reflectance from the refined VIRR method (Fig. 3.18C) shows that almost all contamination due to clouds was removed.

(B)

LTDR AVHRR, January 9, 2010

(D)

LTDR AVHRR, July 12, 2010

FIGURE 3.18 RGB images of the GLASS AVHRR and LTDR AVHRR surface reflectance on January 9 and July 12, 2010. The RGB images were produced with an RGB color scheme that used the red band (as red color), the near-infrared band (as green color), and the red band (as blue color), and the same color enhancement was used.

128

3. Compositing, smoothing, and gap-filling techniques

The VIRR method is a fully automated procedure and provides fast smoothing of data. Therefore, the method is particularly suitable for operational applications with a long time series of global data. Another advantage of the refined VIRR method is that it simultaneously reconstructs time series of surface reflectance in the red and NIR bands, which avoids inconsistent estimation of surface reflectance in the different bands.

3.3 Summary Time series of remotely sensed land surface parameters act as important model inputs when studying global environmental change. However, noise in the data impairs the characterization and prediction performance of models. Through multiphase data synthesis, smoothing, and gap filling, time series of land surface parameters can be reconstructed to improve the quality of the existing land surface parameter products. On the other hand, these methods have their respective shortcomings. For instance, empirical analysis is required to determine the smoothing window size in SG filtering. Current synthesis products have lower atmospheric interference than the originals, but atmospheric influences are not totally eliminated during the synthesis process. Thus, these compositing techniques may not be applicable for improving land surface parameter products in areas with data missing over long time periods (such as tropical rain forests).

References Cabral, A., De Vasconcelos, M.J.P., Pereira, J.M.C., Bartholome, E., Mayaux, P., 2003. Multi-temporal compositing approaches for SPOT-4 VEGETATION. Int. J. Remote Sens. 24, 3343e3350. Cao, R.Y., Chen, Y., Shen, M.G., Chen, J., Zhou, J., Wang, C., Yang, W., 2018. A simple method to improve the quality of NDVI time-series data by integrating spatiotemporal

information with the Savitzky-Golay filter. Remote Sens. Environ. 217, 244e257. Carreiras, J.M.B., Pereira, J.M.C., 2005. SPOT-4 VEGETATION multitemporal compositing for land cover change studies over tropical regions. Int. J. Remote Sens. 26, 1323e1346. Carreiras, J.M.B., Pereira, J.M.C., Shimabukuro, Y.E., Stroppiana, D., 2003. Evaluation of compositing algorithms over the Brazilian Amazon using SPOT-4 VEGETATION data. Int. J. Remote Sens. 24, 3427e3440. Chen, J., Jonsson, P., Tamura, M., Gu, Z., Matsushita, B., Eklundh, L., 2004. “A simple method for reconstructing a high-quality NDVI time-series data set based on the SavitzkyeGolay filter. Remote Sens. Environ. 91, 332e344. Chen, J.M., Deng, F., Chen, M.Z., 2006. Locally adjusted cubic-spline capping for reconstructing seasonal trajectories of a satellite-derived surface parameter. IEEE Trans. Geosci. Remote Sens. 44 (8), 2230e2238. Chuvieco, E., Ventura, G., Martin, M.O., Gomez, I., 2005. Assessment of multitemporal compositing techniques of MODIS and AVHRR images for burned land mapping. Remote Sens. Environ. 94, 450e462. Cihlar, J., Manak, D., D’Iorio, M., 1994. Evaluation of compositing algorithms for AVHRR over land. IEEE Trans. Geosci. Remote Sens. 32, 427e437. Ding, C., Liu, X.N., Huang, F., 2017. Temporal interpolation of satellite-derived leaf area index time series by introducing spatial-temporal constraints for heterogeneous grasslands. Remote Sens. 9, 12. Fang, H., Liang, S., Townshend, J.R., Dickinson, R.E., 2008. Spatially and temporally continuous LAI data sets based on and integrated filtering method: examples from North America. Remote Sens. Environ. 112, 75e93. Fang, H., Liang, S., Kim, H., Townshend, J.R., Schaaf, C.L., Strahler, A.H., Dickinson, R.E., 2007. Developing a spatially continuous 1 km surface albedo data set over North America from Terra MODIS products. J. Geophys. Res. Atmos. 112, D20206. Fernandes, R., Butson, C., Leblanc, S.G., Latifovic, R., 2003. Landsat-5 TM and Landsat-7 ETMþbased accuracy assessment of leaf area index products for Canada derived from SPOT-4 VEGETATION data. Can. J. Remote Sens. 29 (2), 241e258. Gao, F., Morisette, J.T., Wolfe, R.E., Ederer, G., Pedelty, J., Masuoka, E., Myneni, R., Tan, B., Nightingale, J., 2008. An algorithm to produce temporally and spatially continuous MODIs-LAI time series. IEEE Geosci. Remote Sens. Lett. 5 (1), 60e64. Garcia, D., 2010. Robust smoothing of gridded data in one and higher dimensions with missing values. Comput. Stat. Data Anal. 54, 1167e1178.

References

Hermance, J.F., 2007. Stabilizing high-order, non-classical harmonic analysis of NDVI data for average annual models by damping model roughness. Int. J. Remote Sens. 28, 2801e2819. Holben, B.N., 1986. Characteristics of maximum-value composite images for temporal AVHRR data. Int. J. Remote Sens. 7, 1435e1445. Illera, P., Fernandez, A., Delgado, J.A., 1996. Temporal evolution of the NDVI as an indicator of forest fire danger. Int. J. Remote Sens. 17 (6), 1093e1105. James, M.E., Kalluri, S.N.V., 1994. The Pathfinder AVHRR land data set: an improved coarse resolution data set for terrestrial monitoring. Int. J. Remote Sens. 15 (17), 3347e3363. Jiang, Z.Y., Huete, A.R., Didan, K., Miura, T., 2008. Development of a two-band enhanced vegetation index without a blue band. Remote Sens. Environ. 112, 3833e3845. J€ onsson, P., Eklundh, L., 2002. Seasonality extraction by function fitting to timeseries of satellite sensor data. IEEE Trans. Geosci. Remote Sens. 40 (8), 1824e1832. J€ onsson, P., Eklundh, L., 2004. TIMESATda program for analyzing time-series of satellite sensor data. Comput. Geosci. 30, 833e845. Julien, Y., Sobrino, J.A., Verhoef, W., 2006. Changes in land surface temperatures and NDVI values over Europe between 1982 and 1999. Remote Sens. Environ. 103, 43e55. Kasischke, E.S., French, N.H.F., Harrell, P., Christensen, N.L., Ustin, S.L., Barry, D., 1993. Monitoring of wildfires in Boreal Forests using large area AVHRR NDVI composite image data. Remote Sens. Environ. 45, 61e71. Los, S.O., Collatz, G.J., Sellers, P.J., Malmstrom, C.M., Pollack, N.H., DeFries, R.S., et al., 2000. A global 9-yr biophysical land surface dataset from NOAA AVHRR data. J. Hydrometeorol. 1 (2), 183e199. Lu, X., Liu, R., Liu, J., Liang, S., 2007. Removal of noise by wavelet method to generate high quality temporal data of terrestrial MODIS products. Photogramm. Eng. Remote Sens. 73, 1129e1139. Luo, Y., Trishchenko, A.P., Khlopenkov, K.V., 2008. Developing clear-sky, cloud and cloud shadow mask for producing clear-sky composites at 250-meter spatial resolution for the seven MODIS land bands over Canada and North America. Remote Sens. Environ. 112, 4167e4185. Moody, E.G., King, M.D., Platnick, S., Schaaf, C.B., Gao, F., 2005. Spatially complete global spectral surface albedos: value-add datasets derived from Terra MODIS land products. IEEE Trans. Geosci. Remote Sens. 43, 144e157. Nagol, J.R., Vermote, E.F., Prince, S.D., 2014. Quantification of impact of orbital drift on inter-annual trends in AVHRR NDVI data. Remote Sens. 6, 6680e6687. Pedelty, J., Vermote, E., Devadiga, F.S., Roy, D., Schaaf, C., Privette, J., et al., 2007. Generating a long-term land data record from the AVHRR and MODIS instruments.

129

In: IEEE International Geoscience and Remote Sensing Symposium, pp. 1021e1025. Pereira, J.M.C., Sa, A.C.L., Sousa, A.M.O., Martin, M.P., Chuvieco, E., 1999. Regional-scale burnt area mapping in Southern Europe using NOAA-AVHRR 1 km data. In: Chuvieco, E. (Ed.), Remote Sensing of Large Wildfires in the European Mediterranean Basin. Berlin SpringerVerlag, pp. 139e155. Peters, A.J., Walter-Shea, E.A., Ji, L., Vine a, A., Hayes, M., Svodoba, M.D., 2002. Drought monitoring with NDVIbased standardized vegetation index. Photogramm. Eng. Remote Sens. 62 (1), 71e75. Pinzon, J., Tucker, C., 2014. A non-stationary 1981-2012 AVHRR NDVI3g time series. Remote Sens. 6, 6929e6960. Potter, C.S., Brooks, V., 2000. Global analysis of empirical relations between annual climate and seasonality of NDVI. Int. J. Remote Sens. 19 (15), 2921e2948. Pu, R.L., Gong, P., 2004. Wavelet transform applied to EO-1 hyperspectral data for forest LAI and crown closure mapping. Remote Sens. Environ. 91, 212e224. Qi, J., Huete, A.R., Hood, J., Kerr, Y., 1993. Compositing multitemporal remote sensing data sets. Pecora 12, 206e213. Sioux Falls7 American Society of Photogrammetry and Remote Sensing. Qi, J., Kerr, Y., 1997. On current compositing algorithms. Remote Sens. Rev. 15, 235e256. Ranchina, T., Aiazzib, B., Alparonec, L., Barontib, S., Wald, L., 2003. Image fusion-The ARSIS concept and some successful implementation schemes. ISPRS J. Photogrammetry Remote Sens. 58, 4e18. Roerink, G.J., Menenti, M., Verhoef, W., 2000. Reconstructing cloudfree NDVI composites using Fourier analysis of time series. Int. J. Remote Sens. 21, 1911e1917. Rousseeuw, P.J., Croux, C., 1993. Alternatives to the median absolute deviation. J. Am. Stat. Assoc. 88, 1273e1283. Sakamoto, T., Yokozawa, M., Toritani, H., Shibayama, M., Ishitsuka, N., Ohno, H., 2005. A crop phenology detection method using time-series MODIS data. Remote Sens. Environ. 96, 366e374. Savitzky, A., Golay, M.J.E., 1964. Smoothing and differentiation of data by simplified least squares procedures. Anal. Chem. 36, 1627e1639. Schaaf, C.B., Gao, F., Strahler, A.H., Lucht, W., Li, X., Tsang, T., et al., 2002. First operational BRDF, albedo and nadir reflectance products from MODIS. Remote Sens. Environ. 83, 135e148. Sellers, P., Tucker, C., Collatz, G., Los, S., Justice, C., Dazlich, D., Randall, D., 1994. A global 1 by 1 NDVI data set for climate studies. Part 2: the generation of global fields of terrestrial biophysical parameters from the NDVI. Int. J. Remote Sens. 15, 3519e3546. Sobrino, J.A., Julien, Y., 2016. Exploring the validity of the long-term data record V4 database for land surface monitoring. IEEE J. Sel. Top. Appl. Earth Obs. Remote Sens. 9, 3607e3614.

130

3. Compositing, smoothing, and gap-filling techniques

Sousa, A.M.O., Pereira, J.M.C., Silva, J.M.N., 2002. Evaluating the performance of multitemporal image compositing algorithms for burned area analysis. Int. J. Remote Sens. 24, 1219e1236. Tang, H., Yu, K., Hagolle, O., Jiang, K., Geng, X., Zhao, Y., 2013. Cloud detection method based on a time series of MODIS surface reflectance images. Int. J. Digit. Earth 6, 157e171. Tian, F., Fensholt, R., Verbesselt, J., Grogan, K., Horion, S., Wang, Y., 2015. Evaluating temporal consistency of long-term global NDVI datasets for trend analysis. Remote Sens. Environ. 163, 326e340. Tucker, C.J., Pinzon, J.E., Brown, M.E., Slayback, D.A., Pak, E.W., Mahoney, R., et al., 2005. An extended AVHRR 8-km NDVI dataset compatible with MODIS and SPOT vegetation NDVI data. Int. J. Remote Sens. 26, 4485e4498. Van Leeuwen, W.J.D., Huete, A.R., Laing, T.W., 1999. MODIS vegetation index compositing approach: a prototype with AVHRR data. Remote Sens. Environ. 69, 264e280. Vermote, E.F., Vermeulen, A., 1999. Atmospheric Correction Algorithm: Spectral Reflectances (MOD09). MODIS Algorithm Theoretical Basis Document.

Viovy, N., Arino, O., Belward, A.S., 1992. The best index slope extraction (bise)da method for reducing noise in NDVI time-series. Int. J. Remote Sens. 13, 1585e1590. Walthall, C.L., Norman, J.M., Welles, J.M., Campbell, G., Blad, B.L., 1985. Simple equation to approximate the bidirectional reflectance from vegetative canopies and bare soil surfaces. Appl. Opt. 24, 383e387. Xiao, Z., Liang, S., Wang, T., Liu, Q., 2015. Reconstruction of satellite-retrieved land-surface reflectance based on temporally-continuous vegetation indices. Remote Sens. 7, 9844e9864. Xiao, Z., Liang, S., Tian, X., Jia, K., Yao, Y., Jiang, B., 2017. Reconstruction of long-term temporally continuous NDVI and surface reflectance from AVHRR data. IEEE J. Sel. Top. Appl. Earth Obs. Remote Sens. 10 (12), 5551e5568. Zhang, X., 2015. Reconstruction of a complete global time series of daily vegetation index trajectory from long-term AVHRR data. Remote Sens. Environ. 156, 457e472.