Xihan MU, Tian ZHAO, Gaiyan RUAN, Jinling SONG, Jindi WANG, Guangjian YAN,

Tim R. MCVICAR3, Kai YAN1,2,4, Zhan GAO5, Yaokai LIU6, and Yuanyuan WANG7

1 State Key Laboratory of Remote Sensing Science, Faculty of Geographical Science, Beijing Normal University, Beijing 100875, China

2 Beijing Engineering Research Center for Global Land Remote Sensing Products, Faculty of Geographical Science,Beijing Normal University, Beijing 100875, China

3 CSIRO Land and Water, GPO Box 1700, Canberra, ACT 2601, Australia

4 Department of Earth and Environment, Boston University, Boston, MA 02215, USA

5 Jiangsu Institute of Urban Planning and Design, Nanjing 210036, China

6 Aerospace Information Research Institute, Chinese Academy of Sciences, Beijing 100094, China

7 National Satellite Meteorological Center, China Meteorological Administration, Beijing 100081, China

(Received February 17, 2020; in final form July 31, 2020)

ABSTRACT High spatial resolution and high temporal frequency fractional vegetation cover (FVC) products have been increasingly in demand to monitor and research land surface processes. This paper develops an algorithm to estimate FVC at a 30-m/15-day resolution over China by taking advantage of the spatial and temporal information from different types of sensors: the 30-m resolution sensor on the Chinese environment satellite (HJ-1) and the 1-km Moderate Resolution Imaging Spectroradiometer (MODIS). The algorithm was implemented for each main vegetation class and each land cover type over China. First, the high spatial resolution and high temporal frequency normalized difference vegetation index (NDVI) was acquired by using the continuous correction (CC) data assimilation method. Then, FVC was generated with a nonlinear pixel unmixing model. Model coefficients were obtained by statistical analysis of the MODIS NDVI. The proposed method was evaluated based on in situ FVC measurements and a global FVC product(GEOV1 FVC). Direct validation using in situ measurements at 97 sampling plots per half month in 2010 showed that the annual mean errors (MEs) of forest, cropland, and grassland were ?0.025, 0.133, and 0.160, respectively, indicating that the FVCs derived from the proposed algorithm were consistent with ground measurements [R2 = 0.809,root-mean-square deviation (RMSD) = 0.065]. An intercomparison between the proposed FVC and GEOV1 FVC demonstrated that the two products had good spatial–temporal consistency and similar magnitude (RMSD approximates 0.1). Overall, the approach provides a new operational way to estimate high spatial resolution and high temporal frequency FVC from multiple remote sensing datasets.

Key words: fractional vegetation cover (FVC), high spatial resolution and high temporal frequency, data fusion, nor malized difference vegetation index (NDVI), pixel unmixing model, multiple remote sensing datasets

1. Introduction

Fractional vegetation cover (FVC) is defined as the ratio of the vertically projected area of vegetation to the whole area (Yan et al., 2012). In terms of remote sensing techniques, it is most often the ratio of the green vegetation area in a grid cell to the entire pixel area (Lu et al.,2003; Donohue et al., 2009), though recent attempts include the non-photosynthetic component as well (Guerschman et al., 2015), and occasionally defined in the view direction (Purevdorj et al., 1998; Obata et al.,2012). The FVC is sensitive to the vegetation amount and can characterize the surface vegetation status in a horizontal perspective (Gutman and Ignatov, 1998). Additionally, it provides an alternative to the vegetation index for distinguishing the contributions of the soil and the vegetation canopy (Xiao and Moody, 2005; Verger et al., 2015).

Remote sensing is the only practicable way to generate FVC products at regional or global scales on account of its quick and large-scale data acquisition ability (Liang et al., 2013). There are three main methods for estimating FVC using remote sensing data (Xiao and Moody,2005; Jiapaer et al., 2011; Yan et al., 2012; Jia et al.,2015), including: (i) empirical models, (ii) pixel unmixing models, and (iii) physical models. Among these methods, the pixel unmixing model estimates FVC at a subpixel level by decomposing a pixel into at least two portions: (a) green vegetation and (b) non-green background. Linear unmixing modeling assumes that the spectral reflectance or spectral vegetation index (VI) of a pixel is the linearly weighted combination of these two components (Zeng et al., 2000; Yan et al., 2012; Jia et al., 2015). The VI-based mixture model is the most widely used linear unmixing model in high resolution FVC estimation due to its simple model form and computational efficiency when processing large datasets(Gutman and Ignatov, 1998; Zeng et al., 2000; Lu et al.,2003; Montandon and Small, 2008; Wu et al., 2014).

Considerable progress has been made in the generation of FVC products with remote sensing techniques over regional and global scales (Table 1). However, three main issues limiting the widespread use of these products for practical applications remain. First, some products only contain limited land cover types (DeFries et al.,1999; Guan et al., 2012; Sexton et al., 2013; Broxton et al., 2014). Second, some products suffer from a relatively low spatial and/or temporal resolution (Gutman and Ignatov, 1998; Zeng et al., 2000). The operational FVC products, such as the Carbon Cycle and Change in Land Observational Products from an Ensemble of Satellites(CYCLOPES; Baret et al., 2007), GEOLAND2 (Baret et al., 2013), and Polarization and Directionality of the Earth’s Reflectances (POLDER; Roujean and Lacaze,2002; Lacaze et al., 2003), derived from Systeme Probatoire d’Observation de la Terre (SPOT)-VEGETATION and Advanced Earth Observation Satellite (ADEOS)-POLDER data, provided FVCs with global coverage at 10-day and 1–6-km resolutions. Third, efforts have been made to improve the temporal continuity in FVC estimation over large scales (Wu et al., 2014; Ding et al.,2015; Jia et al., 2015; Verger et al., 2015). However, the spatial resolution still requires substantial improvements to meet many operational requirements. Using the lowresolution (herein 1-km) FVC might introduce a maximum error of about 0.35 due to the scale effect in heterogeneous regions (Zhang et al., 2006). For numerous applications, such as monitoring the dynamics of impervious surface area (Zhang et al., 2013), quantitatively analyzing the process of urban expansion and its impacts(Zhang et al., 2015), monitoring the change of urban greenness (Gan et al., 2014), and land management(Naqvi et al., 2013; Pan and Wen, 2014), FVC products with a spatial resolution of 30 m, or better, are needed.Table 1 lists the existing regional or global FVC datasets;note that there is no FVC product at high spatial resolution (herein defined as 30 m) and high temporal frequency (herein defined as once every 15 days). Filling this niche provides the motivation for this study.

High spatial resolution and temporally dense series of remotely sensed data are difficult to obtain because of the long revisit cycles of the satellites, frequent cloud contamination, and other poor atmospheric conditions (Price,1992; Gao et al., 2006; Zhu et al., 2010; Emelyanova et al., 2013). Many algorithms have been proposed to simultaneously enhance resolution in space and time by combining high spatial resolution images with high temporal frequency data. A spatial and temporal adaptive reflectance fusion model (STARFM; Gao et al., 2006) and its enhancement (ESTARFM; Zhu et al., 2010) are the most widely-used data fusion algorithms for simulating daily Landsat-like surface reflectance (Emelyanova et al.,2013). However, these methods suffer from expensive computation and low accuracies in certain situations(Jarihani et al., 2014). Thus, the fusion methods of normalized difference vegetation index (NDVI) series were introduced in temporal–spatial fusion to generate high spatial resolution NDVI time series (Busetto et al., 2008;Cai et al., 2011; Li et al., 2014).

This paper aims to provide a practical approach to retrieving high spatial resolution and high temporal frequency FVC over a large scale using multiple remote sensing datasets. The continuous correction (CC; Cai et al., 2011) model is used to blend the high and low resolution NDVI data and to generate the basic data for FVC estimation. Then, the FVC product was generated with a nonlinear pixel unmixing model by applying the NDVI to FVC transformation coefficients. This paper is organized as follows. Section 2 provides a general description of the study area and datasets. Section 3 details the approaches for retrieving and validating the FVC dataset,and Section 4 reports our results. Discussion and concluding remarks are presented in Sections 5 and 6, respectively.

Definition Fraction of green vegetation cover Fraction of green vegetation cover Fraction of green vegetation cover Fraction of green vegetation cover Fraction of green vegetation cover Fraction of green vegetation cover Fractions of tree, grass,and bare soil Fraction of green vegetation cover Vegetation Continuous Fields of tree cover(percent tree cover)Fraction of green vegetation cover Fraction of green vegetation cover Fraction of green vegetation cover Table 1. Brief summary of the primary existing global or regional FVC datasets. Studies are ordered chronologically and then alphabetically Algorithm Neural network Algorithm (MSUA)works (GRNNs)Pixel unmixing model Period 1985–1990 Pixel unmixing model 1992–1993 Pixel unmixing model 1996–1997,2003 2005–current A statistical approach 2002–2012 Neural network 1999–2007 Neural network 1999–2008 Mean-Sensitivity Unmixing 1998–current Neural network 2000, 2005 Linear regression 1982–2011 Pixel unmixing model 2000–current General regression neural net-2010 Area Southern Africa, and South America Spatial resolution/Temporal Global Global Global Europe, Northern Africa,Europe Global Southeastern Africa Global Global Global Global China frequency 0.15°, 30 days 1 km, 30 days 6 km, 10 days 3 km, 1 day 0.3 km, 30 days/10 days 1 km, 10 days 10 km, 15 days 1/112°, 10 days 30 m, yearly 0.5 km, 8 days 30 m, 15 days Sensor AVHRR AVHRR POLDER SEVIRI MERIS SPOT-VEGETATION AVHRR GIMMS, QuikSCAT SIR, AMSR-E, TRMM SPOT-VEGETATION MODIS, Landsat AVHRR, GIMMS, MODIS 10 km, 30 days MODIS, Landsat (TM,ETM+)MODIS, HJ-1 Study 1–Gutman and Ignatov, 1998 2–Zeng et al., 2000 3–Lacaze et al., 2003;Roujean and Lacaze, 2002 4–García-Haro et al., 2005a, b 5–Bacour et al., 2006 6–Baret et al., 2007 7–Guan et al., 2012 8–Baret et al., 2013 9–Sexton et al., 2013 10–Wu et al., 2014 11–Jia et al., 2015 12–This study Note: AVHRR: Advanced Very High Resolution Radiometer; SEVIRI: Spinning Enhanced Visible and Infrared Imager; MERIS: Medium Resolution Imaging Spectrometer; GIMMS: Global Inventory Modeling and Mapping Studies; QuikSCAT SIR: Quick Scatterometer images using the scatterometer image reconstruction algorithm; AMSR-E: Advanced Microwave Scanning Radiometer-Earth Observing System; TRMM: Tropical Rainfall Measuring Mission; MODIS: Moderate Resolution Imaging Spectroradiometer; TM: Thematic Mapper; ETM+: Enhanced Thematic Mapper.

2. Study area and data

2.1 Study area

The study region includes Chinese mainland and the island of Taiwan, covering over 9.6 million km2. The climate over this region varies regionally and seasonally due to the East Asian monsoon and complex topography,which affects the spatial distribution of vegetation biomes (Fig. 1). The topography of China varies greatly from highly mountainous regions to inhospitable deserts and flat, fertile plains. It is analogous to a staircase descending from west to east. The main vegetation biomes in eastern China, including regions numbered 1, 2, 3, 4,5, 7, 9, 12, and 14 in Fig. 1, are forest and cropland.Whereas, the dominant land cover types in western China, containing regions numbered 6, 8, 10, 11, and 13 in Fig. 1, are grassland, Gobi desert, and saline–alkali land. Additionally, there are some artificial oases in regions 10 and 11, containing irrigated crops and pastures.The percentages of forest, cropland, grassland, and unused land (desert and bare soil) in China are 24.8%,15.6%, 31.0%, and 23.2% in 2010, respectively (Zhang et al., 2014).

2.2 In situ FVC measurements

Fig. 1. The natural vegetation regionalization map of China; this represents ideal (or climax) vegetation if there was no human modification.The 14 vegetation regions are numbered. The solid red triangle and solid red circles represent the Huailai site and other 22 FVC monitoring sites,respectively.

Time-series characteristics of the in situ FVC observations acquired at 22 FVC monitoring sites and the Huailai experimental site are reported in Table 2 and their locations are provided in Fig. 1. The fortnightly measurements at 22 sites commenced in the first half month in January 2010 and ceased in the last half month in December 2010. The time series in situ measurements were used here to assess the seasonality of the FVC product. At Huailai experimental site, where the area is relatively homogeneous at a large scale, measurements were collected to assess the performance of the FVC products at 1-km resolution in Section 4.2. The characteristics of these FVC measurements are illustrated in Table 2.

2.2.1Time-series FVC measurements in smallwatersheds

The time-series FVC measurements were conducted in 22 small watersheds mainly over eastern China, where denser vegetation is located compared to western China.The measurements served the First National Water Resource Survey of China in 2010 for the assessment of water erosion. In each watershed, five to seven sampling plots were chosen for FVC ground measurements over cropland, improved grassland, natural grassland, forest,and orchard. The number of sampling plots in each watershed depends on the richness of land cover types therein (Fig. 2). The measured FVC of a specific sampling plot would not be used if the surveyed land cover type differed from that delineated in the 30-m resolution land cover map (see Section 2.4.1). After filtering the plots that have mismatched land cover types, 97 sampling plots were available for validation (Table 2).

To measure in situ FVC, digital photos were acquired by using a camera mounted on scaffolding. In each plot,5 photos covering 15 m × 15 m were manually acquired for vegetation lower than 2 m. The FVC extraction method generally introduced an absolute error less than 5%(Yan et al., 2012). FVC of the sampling plot was calculated as the average of FVC derived from the five photos.The measurement frequency was once per 15 days, and hence there were 24 observations within the year at each FVC monitoring plot (Table 1).

For trees in the orchards and forests, we placed the digital camera between the overstory and understory, and then downward-looking photos were acquired from the portable scaffolding to characterize the understory vegetation, while extra upward-looking photos were acquired to capture the tree crown. The FVC was then obtained by using:

where fupand fdownare the FVCs extracted from the photographs captured by the upward and downward photography, respectively.

2.2.2Measurements at a homogeneous region

Huailai experimental site (40°21′N, 115°47′E) is located within a flat plain covering ~280 km2. Summer corn covers more than 80% of the experimental region, with other types of agriculture (e.g., fruit orchard, vineyard)covering the remaining ~20%. Measurements covered a complete growth cycle of corn (Zea mays) at 14 sampling plots that were each 30 m × 30 m (Table 2).FVC measured at a sampling plot can represent the FVC of a large homogeneous area. Ground measurements at this site were also obtained by using the digital photography method as outlined above.

2.3 Remote sensing data

We used the HJ-1 multispectral reflectance, MODIS Nadir Bidirectional Reflectance Distribution Function(BRDF)-Adjusted Reflectance (NBAR) product(MCD43B4) and the corresponding quality description product (MCD43B2; https://lpdaac.usgs.gov/products/mcd43a4v006/), MODIS Land Cover Type product(MCD12Q1), and GEOV1 FVC product. The main dataset characteristics are summarized in Table 3. The HJ-1 multispectral reflectance, Collection 5 MCD43B4,MCD43B2, and MCD12Q1 were used to generate the high spatial resolution and high temporal frequency NDVI dataset and the GEOV1 FVC was used to trainone of the coefficients of a nonlinear pixel unmixing model. Additionally, the GEOV1 FVC product and the Global Land Surface Satellite (GLASS) FVC product(Jia et al., 2015) were compared with the product over the Huailai experimental site.

Table 2. Characteristics of in situ FVC measurements used for validation

Table 3. Characteristics of the remote sensing datasets

2.3.1HJ-1 multispectral reflectance

The HJ-1 satellites (HJ-1AandHJ-1B; similar to Landsat in spectral band and spatial resolution) were launched in 2008. The HJ-1 multispectral surface reflectance data (post radiometric calibration and systematic geometric correction) were provided by China Centre for Resources Satellite Data and Application (http://www.cresda.com/EN/). The spatial resolutions of the red band and the near infrared (NIR) band are 30 m. The whole study region consists of 762 tiles with the size of 100 km~north–south × 150 km ~east–west each. HJ-1 images were acquired three times during 2010 (the specific acquisition dates are not fixed). This dataset would be used to calculate instantaneous HJ-1 NDVI and to generate time-series NDVI for data fusion. The HJ-1 data are intensively selected from day-of-year (DOY) 90 to 310 in 2010 (Fig. 3) when vegetation changes more rapidly than the rest days of the year.

Fig. 2. Schematic diagram of a watershed, sampling plot, and digital photography.

2.3.2MODIS reflectance and land cover

MODIS products we used were downloaded from the Land Processes Distributed Active Archive Center (LPDAAC; https://lpdaac.usgs.gov/), including surface reflectance (MCD43B4 and MCD43B2) from 2005 to 2010 and land cover (MCD12Q1) in 2004. MCD43B4 is generated from a semi-empirical kernel-driven bidirectional reflectance model to rectify the actual view direction to the nadir direction (Lucht et al., 2000). Band 1 (red) and band 2 (NIR) were used to calculate the 1-km MODIS NDVI, which was used to form the background field for data fusion and to compute the coefficients of the nonlinear pixel unmixing model. MCD43B2 describes the overall quality of MCD43B4 using quality assessment (QA)values. Only pixels of highest quality for both of the two bands (QA = 0) were used here. MODIS land classification (MCD12Q1) is derived by a supervised decision-tree classification method (Channan et al., 2014). We merged International Geosphere–Biosphere Programme (IGBP)(layer 1) categories to match the classification system of the 30-m national land use/cover (LC, hereafter) database of China introduced in Section 2.4.1. These two land cover datasets were used to acquire the area percentage of each land cover type within a MODIS pixel.

2.3.3GEOV1 FVC product

The GEOV1 FVC product (distributed at http://land.copernicus.eu/global/products/FCover) is an improved version of the CYCLOPES FVC product, which was trained by an artificial neural network approach (Baret et al., 2013). The GEOV1 FVC product performed well when assessing FVC dynamics (i.e., showing reliable spatial distribution and smooth temporal profiles; Camacho et al., 2013), although there was still relatively high residual uncertainty (about 0.1 for cropland; Mu et al., 2015). In this paper, it was only used to fit the nonlinear trend of the nonlinear pixel unmixing model (Section 3.2.2) and compare with the generated FVC product.

Fig. 3. The day-of-year (DOY) frequency distribution of the collected HJ-1 multispectral reflectance data.

2.3.4GLASS FVC product

The GLASS FVC product was generated from MODIS surface reflectance by training regression neural networks with Landsat TM/ETM+ data. For spatial and temporal continuity, the GLASS FVC product was superior to the GEOV1 FVC product (Jia et al., 2015). We used the GLASS FVC product (distributed at http://www.geodata.cn/datapplication/OrderStepList.html) to validate the proposed product over the Huailai experiment site.

2.4 Auxiliary data

2.4.1The 30-m resolution land cover data

An all-China 30-m land cover dataset was updated in 2010 by visual interpretation based on professional fieldbased ecological knowledge. The input data included Landsat TM images with 30-m spatial resolution and, as alternative, the 20-m China–Brazil Earth Resources Satellite (CBERS) data, as well as the 30-m HJ-1 data.The accuracy of the land cover data is about 95.41%(Zhang et al., 2014). Land cover data were used to provide the land cover type of each 30-m resolution pixel, and to acquire the area percentages of land cover types in MODIS pixels.

2.4.2Vegetation regionalization

Vegetation regionalization is a theoretical integration of vegetation studies, mainly based on principles of geographical distribution of vegetation types (Zhang, 1993).It shows the regional distribution and the zonal differentiation of vegetation, and further indicates the regularity of the distribution of vegetation and its relationship with the environment. Chinese vegetation regions were obtained from Wang and Zuo (2010) and numbered from 1 to 14 (Fig. 1). The 14 vegetation regions (Fig. 1) are the basic units to calculate NDVI background and coefficients that are used to transform NDVI into FVC.

3. Methods

We employed the CC data assimilation method and a nonlinear pixel unmixing model (Cai et al., 2011) to estimate the high spatial resolution and high temporal frequency FVC product using MODIS NBAR product, HJ-1 reflectance data, and GEOV1 FVC product respectively for each land cover type and each vegetation region. Figure 4 shows the general processes of the algorithm: acquisition of the high spatial resolution and high temporal frequency NDVI dataset (30-m NDVI, hereafter), training of the transformation coefficients from NDVI to FVC for different vegetation regions and land cover types, and calculation and validation of FVC, corresponding to the following Sections 3.1, 3.2, and 3.3, respectively.

3.1 Generation of high spatial resolution and high temporal frequency NDVI

The 30-m NDVI dataset was acquired by using the CC data assimilation method. We used HJ-1 NDVI as the high resolution input data, and 8-day multi-year (2005–2010) averaged MODIS NDVI time series as the background data (Cai et al., 2011) that represent vegetation phenology for different land cover types.

3.1.1HJ-1 NDVI and MODIS NDVI time series

NDVI is one of the most frequently used vegetation indices to characterize the vegetation status, defined by Eq. (2). HJ-1 NDVI and MODIS NDVI were respectively derived from HJ-1 multispectral reflectance data in 2010 and the MCD43B4 product during 2005–2010.

where ρRand ρNIRare the reflectance of the red band and NIR band.

The geo-registered MCD12Q1 and 30-m land cover data were processed to output the area percentage image by computing the ratio of the number of high resolution pixels for a specific land cover type in a coarse pixel to the number of all high resolution pixels in this coarse pixel (Cai et al., 2011). Only the pixels with the area ratio greater than 95% (pure pixels, hereafter) were used.The multi-year averaged MODIS NDVI time series were acquired by averaging the 8-day-frequency MODIS NDVI time series over the six years, which were extracted over the pure pixels for each land cover type and vegetation region.

3.1.2High spatial resolution and high temporalfrequency NDVI

The CC method for data assimilation adopts the so called “index-then-blend” (i.e., calculate the index at both high and low spatial resolutions and then perform the blending operation) approach, which produces more accurate results than the alternative “blend-then-index”approach (Jarihani et al., 2014), and can fully use the vegetation information from the two data sources simultaneously (Cai et al., 2011; Meng et al., 2011). The model can be expressed as:

where Xp, Xb, and Xhare the predicted NDVI, background field NDVI, and the high resolution input NDVI,respectively; riand rjare the prediction date and the observation date, respectively;are the errors of the high resolution input and background field data, respectively;nis the number of high resolution inputs; and ?(ri,rj)is the weighting factor calculated by the time distance between the background and high resolution input data.

The fused NDVI series for a 30-m pixel was generated at a 15-day frequency with the inputs of high resolution HJ-1 NDVI acquired three times in 2010, background field NDVI (i.e., the multi-year averaged NDVI time series for the vegetation region and land cover type), and the weighting factors. Water body was masked by using land cover types in FVC retrieval. The pixels with snow were not excluded from the HJ-1 NDVI data as the FVC values of snow pixels are generally calculated to be zero given low values of snow NDVI.

3.2 Transformation coefficients and FVC estimation

3.2.1Nonlinear pixel unmixing model

The pixel unmixing model with two endmembers is the simplest and most extensively used model among the various linear unmixing models (Yan et al., 2012). It assumes that each pixel consists of green vegetation and non-green background. The information (vegetation indices or spectral) of the pixel results from the linear weighted synthesis of the two components. The weight of each component is the area proportion of each component in the pixel. The vegetation proportion is the FVC of the pixel, which can be mathematically expressed as(Gutman and Ignatov, 1998):

where FVC is the proportion of vegetation area in the mixed pixel, namely, the FVC of the mixed pixel; NDVI is the NDVI of the mixed pixel; and NDVIvand NDVIsare the NDVI of fully covered green vegetation and nongreen background, respectively. This pixel unmixing model is a linear transformation between NDVI and FVC.

Lu et al. (2003) showed that NDVI was suitable to model FVC. Some studies defined the scaled NDVI [N*,Eq. (6)] and established an identical quadratic relation between the scaled NDVI and FVC [Eq. (7); Choudhury et al., 1994; Carlson and Ripley, 1997]:

It is documented that the linear relationship or the quadratic relationship could reveal good agreement between NDVI and FVC in some cases (Xiao and Moody, 2005; Jiapaer et al., 2011). Instead of using a simple linear or quadratic expression to calculate FVC,an empirical nonlinear function is used to convert NDVI to FVC, which is appropriate for both the linear and quadratic forms (Mu et al., 2015). This model is defined as the nonlinear pixel unmixing model:

where FVC and NDVI are the FVC and NDVI of the mixed pixel, respectively; NDVIvand NDVIsare the NDVI of fully covered green vegetation and non-green background, respectively; andkis the linearity coefficient of the model.

3.2.2Transformation coefficients of the model

It has long been known that NDVIvand NDVIsare biome-specific (Gutman and Ignatov, 1998; Zeng et al.,2000; Montandon and Small, 2008; Wu et al., 2014). The NDVIvdepends upon the vegetation type, geometric structure, chlorophyll content, and physiology (mesophyll) of vegetation (Price, 1992; Carlson and Ripley,1997). The NDVIsis affected by soil reflectance varying with soil types and soil moisture (O’Neill, 1994; Post et al., 2000; Lobell and Asner, 2002; Montandon and Small, 2008). In addition,kvaries according to surface heterogeneity and vegetation type (Mu et al., 2015).Thus, the three coefficients were independently trained for each vegetation region and land cover type over the pure pixels.

The NDVIvand NDVIscan be estimated as maximum and minimum NDVI (i.e., NDVImaxand NDVImin) in the study areas based on the assumption that fully covered vegetation and non-vegetation background exist in the observed time and space (Gutman and Ignatov, 1998).They were determined through spatial and temporal statistical analysis of the 30-m NDVI data, assuming that similar NDVImaxand NDVImincould be obtained for each vegetation region (see Fig. 1) and each land cover type. We re-grouped the land cover types into three main vegetation types: forest, cropland, and grassland in each region. This was performed in three steps. First, the annual maximum NDVI and minimum NDVI of the 30-m NDVI within the extent of the MODIS pure pixels were calculated. Second, the cumulative histogram of the annual maximum NDVI was acquired for the three main vegetation types in each region. The NDVImaxwas taken as the value at 75% of the cumulative histogram for cropland and grassland, and 90% for forest (Zeng et al.,2000). Third, and finally, the averaged value of the annual minimum NDVI for each vegetation type in each vegetation region was defined as the corresponding NDVImin(Wu et al., 2014; Ding et al., 2016). If these three steps fail to meet the criteria for NDVImax(0.70 < NDVImax<0.95) or NDVImin(0.05 < NDVImin< 0.20), 0.84 and 0.07 are assigned to NDVImaxand NDVImin, respectively(Montandon and Small, 2008). Thekdescribes the degree of nonlinearity and was acquired by fitting Eq. (8)with GEOV1 FVC and MODIS NDVI. GEOV1 FVC can reflect the general trends of vegetation growth (Camacho et al., 2013; Ding et al., 2015) and thereby provides reasonable information for the prediction ofk.

3.2.3High spatial resolution and high temporal frequency FVC generation

The FVC product was generated at 30-m resolution and 15-day frequency for different vegetation regions and land cover types by inputting the 30-m NDVI and the transformation coefficients to the nonlinear pixel unmixing model [Eq. (8)].

3.3 Assessment and validation

Assessment and validation of the high spatial resolution and high temporal frequency FVC product were completed through direct validation and intercomparison.The direct validation was conducted by using the in situ FVC measurements at 97 sampling plots in 22 small watersheds and 14 sampling plots at the Huailai experimental site to quantify the overall performance of the product.The intercomparison was carried out between the estimated FVC by the method proposed in Sections 3.1 and 3.2 and the GEOV1 FVC product to analyze their spatial and temporal consistency.

3.3.1In situ measurements

The direct validation was conducted over 22 small watersheds and the Huailai experimental site (Fig. 1). To reduce the potential geo-location errors between the estimated FVC product and the field FVC measurements, we averaged the FVCs of 3 × 3 pixels (each 30-m resolution)centered on the ground point positioning pixel as suggested by Weiss et al. (2007). The bias of the estimated FVC over a sampling plot was calculated as:

where FVCeis the averaged FVC over the 3 × 3 pixels,and F VCris the in situ FVC.

The GEOV1 FVC and GLASS FVC were also compared against the field FVC measurements over the Huailai experimental site.

3.3.2Intercomparison with GEOV1 FVC

The estimated FVC product was respectively mosaicked and resampled to the 1-km spatial resolution and then together with the GEOV1 FVC was monthly synthesized with the arithmetic mean to ensure they have the same temporal compositing periods for intercomparison.

We evaluated the performance of the two products in four ways, including assessing: (i) the spatial and temporal distribution patterns; (ii) differences between the values of the estimated FVC and GEOV1 FVC; (iii) the correlation between the two products; and (iv) the spatial continuity represented by the percentage of missing values. Five statistical metrics, namely, DIFF, Ratio, mean error (ME), root-mean-square deviation (RMSD), and the Pearson correlation coefficient (R), were computed based on the estimated FVC product and the GEOV1 FVC product per pixel over China for 2010; see Eqs. (10–14):

where FVCe(i)and FVCGEO(i)are respectively the estimated FVC and GEOV1 FVC of a pixel atith temporal phase;represent the averaged estimated FVC and the averaged GEOV1 FVC of the pixel, respectively; andnis the number of valid FVCs for the pixel (the number of valid temporal phases), which means that missing values do not contribute ton.

4. Results

4.1 Spatial distribution of the estimated FVC product

Figures 5a–c show the spatial distributions of forest,cropland, and grassland, respectively. The sub-classes under the three main land cover types were specifically defined by Zhang et al. (2014). Forest, shrub, and sparse woods, respectively, consist of trees higher than 2 m with canopy cover greater than 30%, trees less than 2 m in height with canopy cover greater than 40%, and trees with canopy cover of 10%–30%. Other woods represent tea gardens, orchards, groves, and nurseries. Paddy field represents the cropland that has sufficient water supply for planting paddy rice and lotus, while dry land only provides limited irrigation facilities for rain-fed farming crops. Sparse grass, moderate grass, and dense grass, respectively, are the grassland with canopy cover between 5% and 20%, between 20% and 50%, and greater than 50%.

Forests are mainly distributed in Northeast and Southeast China (Fig. 5a), croplands are mainly distributed in central and eastern China with oasis in Northwest China(Fig. 5b), and grasslands are mainly distributed in northern China (Inner Mongolia) with oasis in Northwest China and southwestern Tibetan Plateau (Fig. 5c). FVC of forest is generally higher than that of cropland, and FVC of cropland is higher than that of grassland. There are large deserts in Northwest China (Fig. 1), which have FVC values of almost zero. In general, China’s FVC ascends from west to east, with values approximately from 0 to 1, while the FVC in the oasis class in Northwest China is higher than that of the surrounding regions. The distribution pattern of the annual maximum FVC over China (Fig. 5d) is in agreement with the actual situation,which is substantially affected by the East Asian monsoon and China’s topography (Hu and Zhang, 2006). The water bodies were masked (e.g., Fig. 5d).

4.2 Direct validation

Figure 6 shows the systematic differences between the in situ FVC measurements and the estimated FVC product during 2010 for all FVC monitoring sites in the 22 small watersheds. Basic statistics of the bias time series throughout the whole year are presented with mean and ± 1 standard deviation in Fig. 6. The bias of the estimated FVC product over cropland (Fig. 6b) presents a temporal fluctuation with a general trend of overestimation, while a slight underestimation is found for July and August. The bias over grassland also shows overestimation (Fig. 6c), and the mean bias of the forest FVC is relatively small, ranging from ?0.081 to 0.075 (Fig. 6a).

Fig. 5. Spatial distributions of the three main vegetation types: (a) forest, (b) cropland, and (c) grassland, and (d) the annual maximum FVC over China in 2010. The Albers equal-area conic projection is used and same for Figs. 8, 9, and 11 in this paper.

Fig. 6. The mean bias of the estimated FVCs over FVC monitoring sites for each temporal phase and vegetation type in 2010. (a)–(c) represent the results over forest, cropland, and grassland, respectively; n is the number of measurements involved in the results. The black curves are the time series of mean bias at each half a month, commencing in January 2010. The upper and lower boundaries of the gray ribbon are determined by adding a positive and negative standard deviation, respectively, to the mean bias.

All in situ FVC measurements over the Huailai experiment site were used to validate the proposed method. A good agreement (R2= 0.809, RMSD = 0.065) is observed between the FVC estimates and the ground measurements (Fig. 7). The performance of the proposed FVC is comparable to that of the GEOV1 FVC (R2= 0.612,RMSD = 0.100) and that of the GLASS FVC (R2=0.713, RMSD = 0.110).

4.3 Intercomparison

Figure 8 shows the seasonal FVC distributions across China from the HJ-1/MODIS-based blending method developed here and the GEOV1 FVC product. January,April, July, and October represent the seasonal FVCs in winter, spring, summer, and autumn, respectively. A good agreement in spatial and temporal distribution between the two products is observed. The seasonal distribution shows reasonable temporal pattern: the high FVC values in summer (Figs. 8c, g) and the relatively low FVC values in winter (Figs. 8a, e). Some differences between the two FVC products exist in spring (Figs. 8b,f), when the change of vegetation is more rapid than in other seasons (Xie and Wilson, 2020) and would cause more uncertainty due to the unequal temporal compositing period (1 month for GEOV1 FVC and half a month for the proposed FVC). The FVC values of snow- or icecovered land and lakes are zero due to their low NDVI values. Figure 9 presents the seasonal variation of the difference between the estimated FVC and the GEOV1 FVC, and the ratio of these two products. It is prominent that there are invalid values in GEOV1 FVC maps (percentages are shown in the legends of Fig. 8), whereas no invalid values exist in all the estimated FVC maps (with the waterbodies masked).

Fig. 7. Scatterplots of the estimated FVC and in situ FVC measurements in Huailai experimental site. The dominant land cover type is cornfield(comprising more than 80% of the experimental region). Cropland also includes the deciduous orchard (covering the remaining 20% of the experimental region).

The monthly time series results of the two FVC products over all pure pixels for the cropland, forest, and grassland (Figs. 10a–c) reveal that both products capture the seasonal phenological variation. The FVC is low and has almost no variation in winter and spring (October to March) when the vegetation stops growing or grows very slowly. Then, the FVC increases gradually in spring(April) when the vegetation starts growing. Afterwards,the FVC reaches the maximum value in summer (July)when the vegetation is flourishing (as limited drought conditions were experienced across China in 2010). Finally, the FVC decreases gradually to the minimum value due to the vegetation decaying and dying.

The two products are relatively consistent in their temporal variation trends; however, there are some important differences. Figure 10a shows that, for forest, in the GEOV1 FVC product, there is a bimodal distribution with a small local maximum being estimated in February,which is not consistent with in situ observations. For cropland (Fig. 10b), the maxima of the estimated FVC developed here occur about a month prior to the maxima in the GEOV1 FVC product, and importantly for crop yield estimation (that is typically related to the growing-season integral), the estimated FVC has higher values for longer time (i.e., it is broader) when compared to the GEOV1 FVC product. Additionally, there exist some differences in magnitude; for example, the estimated FVC is systematically higher than GEOV1 FVC product by up to 0.1 throughout the year for grassland (Fig. 10c).

No missing value was observed in the estimated FVC product because gaps in HJ-1 images were filled with MODIS data in the fusion and generation of NDVI dataset. In contrast, Fig. 8 shows that the GEOV1 FVC (i.e.,the preprocessed SPOT-VEGETATION reflectance data)contained numerous missing values (Camacho et al.,2013). Figure 11a shows the percentage map of the missing values of GEOV1 FVC product in 2010. However,the VEGETATION biogeophysical product version 2(GEOV2) used climatology information from version 1 to apply temporal smoothing and gap filling to improve spatial and temporal continuity and consistency (Verger et al., 2014).

Maps of the calculated statistical metricsR, ME, and RMSD between the time series of estimated FVC and GEOV1 FVC products in 2010 are shown in Figs. 11b–d,respectively. Figure 11b indicates that prominent positive correlations (0.8 ≤R≤ 1) between the two FVC products are observed in most of eastern China. Negative correlations and weak correlations (?0.2 to 0.2) are found in the desert areas in Northwest China and the Tibetan Plateau, where FVC remains low (< 0.1) all year.

A good agreement between the estimated FVC and GEOV1 FVC products was observed in the annual ME,which is mainly distributed from ?0.1 to 0.1 (Fig. 11c).The spatial distribution of RMSD between the estimated FVC and GEOV1 FVC products agrees with that of the annual maximum FVC in Fig. 11d. This means that the regions with low FVC values have low RMSD values and high overall accuracy, and vice versa. The regions with large RMSDs are almost consistent with the regions with the large missing data percentage in the GEOV1 FVC product, suggesting that the RMSD values are partly influenced by the amount and the quality of data involved in GEOV1 FVC product calculation.

Fig. 8. Seasonal spatial distributions of the estimated FVC product developed herein (left column) and the GEOV1 FVC product (right column).Sub-plots (a)–(d) are for January, April, July, and October, respectively, for the estimated FVC product. Sub-plots (e)–(h) are the corresponding monthly GEOV1 FVC estimates. The numbers in parenthesis directly to the right of the word “Invalid” in the legend report the percentage area of China (9.6 × 106 km2) containing invalid (or unavailable) pixels each month. Note there is never any invalid FVC estimate in the product developed herein.

Fig. 9. Seasonal variations of the (a–d) difference and (e–h) ratio between the estimated FVC product and the GEOV1 FVC product for (a, e)January, (b, f) April, (c, g) July, and (d, h) October, respectively. The white areas within each map of China show where GEOV1 FVC estimates are invalid.

Fig. 10. Seasonal variations of mean values of the estimated FVC and GEOV1 FVC products over all pure pixels within 2010 for the three main vegetation types (denoted by the thick solid lines): (a) forest, (b) cropland, and (c) grassland. Statistics are based on the monthly composed FVC products. The width of the ribbon belt represents ± 1 standard deviation from the mean.

4.4 Sensitivity of FVC to NDVImax and NDVImin

The accuracy of the FVC estimate is influenced by the accuracy of NDVImax, NDVImin, andkaccording to Eq.(8). We only analyzed the influences of the NDVImaxand NDVIminon the estimation of FVC (Fig. 12), because the coefficientkis very close to 1, and its variation is negligible. The NDVI varies from 0 to 1 with a sampling interval of 0.2. For each NDVI, the changing ranges of NDVImaxand NDVIminwere set from 0.7 to 0.95 and from 0.05 to 0.2, respectively, with 0.05 as the sampling interval, which were determined according to previous studies (Gutman and Ignatov, 1998; Zeng et al., 2000;Montandon and Small, 2008; Jiménez-Mu?oz et al.,2009; Wu et al., 2014) and the boundary conditions in Section 3.2.2. The NDVImaxand NDVIminreference values, i.e., the values to calculate the reference FVC, were set close to the values used in this study when analyzing the influences of NDVImin(Fig. 12a) and NDVImax(Fig.12b). The deviations of FVC were obtained with the changes of NDVImaxand NDVImin. In Fig. 12a, we could see that FVC decreased with the increase of NDVImin,and the effect of NDVIminon FVC would be greater when the NDVI was smaller. When the NDVI was 0.2,the influence of NDVIminon FVC was the largest among our results, and the deviation of FVC from the reference was up to ?0.2. Similarly, we could see in Fig. 12b that the FVC presented a decreasing trend with the increase of NDVImax.

Fig. 11. Statistical characteristics and assessment between the estimated FVC and GEOV1 FVC products in 2010 over China. Sub-plot (a)shows the percentage of annual missing values of GEOV1 FVC, (b) is the map of R values, ranging from +1.0 to ?1.0, (c) reports the annual ME,and (d) reports the RMSD.

Fig. 12. The deviations of FVC estimate derived by using the method developed herein with the variation of NDVImin and NDVImax. (a) The influence of NDVImin and (b) the influence of NDVImax.

5. Discussion

5.1 Fusion of high resolution and high frequency data

The fusion of high resolution and high frequency data combines the “beneficial” characteristics of each input data types. The output of high frequency temporal changes with high spatial details enables the monitoring of vegetation phenology at fine resolutions (Singh et al.,2011; Bhandari et al., 2012; Walker et al., 2014; Zhang and Wu, 2015). Another strategy to generate high resolution and high frequency data is to fill the gaps in high resolution data by using interpolation techniques (Yang et al., 2017). However, high resolution satellites revisit the same place with longer time intervals. Landsat data are acquired every 16 days, causing large gaps due to the cloudiness, especially in South China (Wilson and Jetz,2016; Cai et al., 2017). Therefore, the data fusion technique, such as that developed herein, generates products with better continuity in space and time (Fig. 8; Weiss et al., 2014). In addition, high resolution (30 m) data are rare, especially before 2014, with the exception of Landsat. The data fusion with HJ-1 data (starting from 2008)is important to generate long time series products, asLandsat-8, Sentinel-2, and Chinese Gaofen satellites were launched after 2013. The method implemented on the fusion of MODIS and HJ-1 data (similar to Landsat in spectral band and spatial resolution) can help to generate a more continuous time-series FVC product.

We implemented the data fusion of NDVI with an “index-then-blend” manner, which is more accurate than a“blend-then-index” manner due to less error propagation(Jarihani et al., 2014). The NDVI blending method developed herein showed a deviation of approximately less than 0.1 from field measurements, better than or at least comparable to the widely used STARFM method (Cai et al., 2011). The uncertainty transferred to FVC will be less than 0.13, based on Eq. (5) with NDVIvand NDVIsof default values (0.84 and 0.07 in Section 3.2.2).

However, the temporal frequency of the fused data primarily depends on the frequency of coarse resolution data that have shorter revisit time (Becker-Reshef et al.,2010). The MODIS NBAR product, composed of 16-day/1-km observations, was used as the background of NDVI in this study, as the product is free of angular effects. With the relatively low frequency of the background data, some discrepancies in the FVC estimates were expected. We found that the FVC had discrepancies up to 0.15 against the in situ FVC in some timephases in Fig. 6. The biases of the estimated FVC product over cropland (Fig. 6b) present relatively large uncertainty, especially from October to December. For grassland, larger biases exist in the start and the end of the year. In Fig. 6, larger biases are also observed in the rapid change period of vegetation, e.g., March to May,and September to October. Low frequency MODIS NBAR product may have missed the temporal information when vegetation changes quickly, and induced errors to the 30-m NDVI dataset. Another issue that affects the performance of data fusion is the limited amount of Landsat-like high resolution data (Roy et al.,2008). For each location, HJ-1 images were only acquired at three temporal phases. The data acquisition dates are non-uniformly distributed in 2010 (Fig. 3).There is a small amount of imagery available for spring and winter when the vegetation changes would not be completely reflected in 30-m FVC, as the primary information is only from the 1-km MODIS NDVI background data. The small amount of HJ-1 high resolution data used for the start and the end of the year could account for the discrepancies to some degree. It is expected that the accuracy of FVC can be promoted with more high spatial resolution observations.

The NDVI background data were selected from 6-yr(2005–2010) MODIS NDVI. The selection involved the QA provided by MODIS product, the correction of BRDF effect (MCD43B4), and the assessment of the spatial homogeneity using 30-m resolution land cover map(Section 2.3.2). We used MODIS land cover product in 2004 (Table 3) and 30-m resolution classification map in 2010 to exclude the influence of land cover changes over the years. Only the pixels with the same land cover type can be the background NDVI candidates.

The introduction of classification data benefits the accuracy of fused high resolution and high frequency data(Fu et al., 2013) with 30-m vegetation classification datasets becoming recently operational at the global scale(Gong et al., 2013; Chen et al., 2015). However, the misclassification of land cover types would degrade the quality of the retrieval of transformation coefficients from NDVI to FVC, particularly between categories with substantially different NDVImaxand NDVImin(e.g.,grassland versus broadleaf forest; Zeng et al., 2000).

5.2 Uncertainty in determining NDVImax, NDVImin, and linearity coefficient

Neural network is popular in generating FVC product at a large scale (Table 1). Pixel unmixing model is also convenient to produce FVC product at regional or global scales whereas the determination of both NDVImaxand NDVIminis important. Figure 12 shows that the values of NDVIminand NDVImaxmostly affect FVC estimates over sparse vegetation (low NDVI) and dense vegetation (high NDVI), respectively. The estimated FVC product overestimates FVC for grassland and cropland (Figs. 6, 10), and it may be caused by the uncertainty of NDVImin. Xiao and Moody (2005) showed that the influence of the soil background would produce an overestimation of FVC in sparsely vegetated area. The overestimation of FVC would be larger than 0.2 if using the invariant NDVImin(about 0.05), compared to using the mean NDVI computed from field measured NDVImin(0.2 ± 0.1) (Montandon and Small, 2008). Similarly, the overestimation of FVC would be 0.07 ± 0.08 when the FVC was estimated by using the invariant NDVImin(about 0.05) in Northeast China, and the largest error occurred at a low NDVI level(Ding et al., 2016). In our study, the NDVI of grassland is usually less than 0.4, which indicates that the uncertainty of NDVIminin this study would substantially affect the estimation of FVC. In addition, the reference FVC, say, in situ measurements and GEOV1 FVC, explains the relative overestimation of FVC for sparse vegetation to some degree: the number of sampling plots is not large over grassland (15 plots) and the underestimation of GEOV1 FVC is found on open grassland and sparse vegetation (Ding et al., 2015). The determination of NDVImaxand NDVIminis expected to be improved by introducing physical models to avoid the uncertainty in statistical methods (Song et al., 2017; Mu et al., 2018).

In the proposed method, the linearity coefficientkis acquired by fitting Eq. (8) with GEOV1 FVC and MODIS NDVI, which is impacted by the quality of GEOV1 FVC and may then influence the estimated FVC.However, the fittedkwas found to be approximately 1 and the linear relationship between NDVI and FVC was widely used in many applications (Gao et al., 2020), indicating that the uncertainty induced bykcould be insignificant.

6. Conclusions

In this study, we proposed a retrieval algorithm for green FVC estimation at high spatial resolution and high temporal frequency by the combination of fine resolution images and high temporal frequency images. We chose the multi-year averaged MODIS NDVI time series as the background field and the HJ-1 NDVI as the highresolution inputs for the CC data assimilation method for each vegetation region and each land cover type. Then,we used the produced high spatial resolution and high temporal frequency NDVI and statistically obtained transformation coefficients from NDVI to FVC to estimate an FVC product at 30-m/15-day resolution within 2010 over China. The comparison was implemented with in situ FVC measurements over 22 small watersheds at Chinese soil and water conservation monitoring stations and a spatially homogeneous experiment site (Huailai).Our product was consistent with the vegetation distribution in China in terms of the temporal and spatial distribution pattern of FVC and was very consistent with the ground-measured FVC over the Huailai site (R2= 0.809,RMSD = 0.065). The FVC estimate does not deviate far from the time series measurements over 22 watersheds,especially for forest (less than 0.1). Furthermore, the FVC maps generated from the proposed method had comparable spatial consistency with the GEOV1 FVC data, and better temporal and spatial continuity.

In summary, the proposed method paves the way to operationally generate high spatial resolution and high temporal frequency FVC products using multisource data. However, there are still some limitations that need rectification for further application. Further work should focus on improving the product quality especially at important phenology periods of vegetation and the determination of coefficients in the pixel unmixing model. The background NDVI time series can be extracted at a smaller scale, e.g., at the pixel scale of MODIS data, to capture more temporal information of vegetation.

Acknowledgments. We thank China Centre for Resources Satellite Data and Application for providing HJ-1 data, Resource and Environment Science and Data Center (RESDC) for providing the land cover data, European Space Agency for providing free GEOV1 data, and the LP-DAAC and MODIS science team for providing free MODIS products. Thanks go to Wenwen Cai, Shuai Huang, and Jingyi Jiang for the data processing. We thank Professor Wenbo Zhang for sharing the FVC on soil erosion in watersheds and Jun Chen for drawing some of the important figures. We especially appreciate the contribution from Academician Xiaowen Li. We thank the two reviewers and Editor for helpful comments that improved an earlier version of this paper.