Scaling of FAPAR from the Field to the Satellite

The fraction of absorbed photosynthetically active radiation (FAPAR) is a critical biophysical parameter in eco-environmental studies. Scaling of FAPAR from the field observation to the satellite pixel is essential for validating remote sensing FAPAR product and for further modeling applications. However, compared to spatial mismatches, few studies have considered temporal mismatches between in-situ and satellite observations in the scaling. This paper proposed a general methodology for scaling FAPAR from the field to the satellite pixel considering the temporal variation. Firstly, a temporal normalization method was proposed to normalize the in-situ data measured at different times to the time of satellite overpass. The method was derived from the integration of an atmospheric radiative transfer model (6S) and a FAPAR analytical model (FAPAR-P), which can characterize the diurnal variations of FAPAR comprehensively. Secondly, the logistic model, which derives smooth and consistent temporal profile for vegetation growth, was used to interpolate the in-situ data to match the dates of satellite acquisitions. Thirdly, fine-resolution FAPAR products at different dates were estimated from the Advanced Spaceborne Thermal Emission and Reflection Radiometer (ASTER) data using the temporally corrected in-situ data. Finally, fine-resolution FAPAR were taken as reference datasets and aggregated to coarse resolution, which were further compared to the Moderate Resolution Imaging Spectroradiometer (MODIS) FAPAR product. The methodology is validated for scaling FAPAR from the field to the satellite pixel temporally and spatially. The MODIS FAPAR manifested a good consistency with the aggregated FAPAR with R2 of 0.922 and the root mean squared error of 0.054.


Introduction
The Fraction of Absorbed Photosynthetically Active Radiation (FAPAR) describes the efficiency of light absorption by plant and is a critical biophysical parameter in eco-environmental studies [1,2].It is defined as the fraction of Photosynthetically Active Radiation (PAR) absorbed by vegetation, where PAR refers to the incoming solar radiation within the spectral range 0.4-0.7 µm [3].FAPAR provides critical information on vegetation growth condition and reflects the energy-and-carbon exchange between biosphere and atmosphere.As a clearly defined quantitative parameter, it has been widely used as a key input into regional or global climate models [1,4].
In view of the difficulties in ground measurements, sensors onboard satellite now provide the only available FAPAR estimates from local to global scales [5][6][7], such as moderate-resolution imaging spectroradiometer (MODIS) [8], Sea-Viewing Wide Field-of-View Sensor (SeaWiFS) [9], SPOT-VEGETATION [10] and the Medium Resolution Imaging Spectrometer (MERIS) [11].The validity of these satellite FAPAR products now makes FAPAR a key terrestrial remote sensing product in the Global Climate Observing System [12].In particular, in Production Efficiency Models (PEMs) for Gross and Net Primary Productivity (GPP/NPP), FAPAR often acts as the only satellite derived variable [5].Because multiple satellite FAPAR products provide users various choices, it is essential to validate these global datasets before application at local or regional scale [13].Large uncertainties resulting from data quality and retrieval accuracies were found in the currently available FAPAR products [5,14].Apparent inconsistencies among these FAPAR products were revealed in comparative studies, such as those conducted in [5,[14][15][16][17][18].These uncertainties or inconsistencies have made continual validation and comparison with ground measurements imperative for real application or reducing the uncertainties.
However, the validation of remote sensing FAPAR products is greatly limited as the spatially and temporally comparable in-situ FAPAR data are very rare [7,19].Especially FAPAR is highly temporally variable, and its diurnal variations affected by the changing irradiation conditions can be dramatic [20].This hindered the collection of FAPAR ground measurements.Most of the executed or ongoing field campaigns were focused on leaf area index (LAI) rather than FAPAR [3,15,21].To take good advantage of the limited and valuable in-situ FAPAR data, scaling of FAPAR from the field to the satellite is critical to achieve the comparability between in-situ and satellite FAPAR data.
To accomplish the scaling, spatial mismatch naturally poses the major challenge.While in-situ FAPAR are representative of several meters at plot scale, a single pixel on the coarse resolution product extends a distance of 250 m-5 km, as defined in [22].As land surface are commonly heterogeneous, the coarsening of spatial resolution results in the mismatches of spatial representative areas between field data and satellite data [3].To guarantee the comparability between in-situ and satellite data, homogeneous sites are always preferred and direct comparison can be performed after the representativeness of field sites to coarse-resolution pixel was verified, as conducted in [13,15,23,24].However, at heterogeneous sites that are more widely distributed, direct comparison between the pair of in-situ and satellite data is unfeasible.This urges the need to consider the sub-pixel variability at coarse resolutions in scaling in-situ data to coarse-resolution satellite product [16].Morisette et al. [25] proposed a bottom-up methodology for scaling ground measurements from field to elementary sampling units (ESUs), fine-resolution data and then moderate or coarse resolution data sequentially.It has been widely used in scaling field observations to coarse-resolution pixels, as those conducted in [9,13,16,26,27].
Spatial discrepancies can be bridged by the use of fine-resolution data, while temporal mismatch is another critical issue and even a more challenging task for scaling FAPAR from the field to the satellite.Two aspects should be considered here: temporal normalization to match the specific time stamp and temporal interpolation to match the acquisition dates of field measurements to satellite data.
The dramatic diurnal variations of FAPAR lead to fluctuations in ground measurements at varying hours.Ideally FAPAR in-situ data are expected to be measured at the satellite overpass time.However, given the limited labor, time and instrument, it is impractical to measure FAPAR for all field sites at the same time.Multiple sites are always measured one by one at different hours.Temporal normalization is thus required to match the in-situ FAPAR to the time of satellite overpass in the scaling-up process.Fensholt et al. [13] found that daily averages of FAPAR from 9:00 a.m. to 3:00 p.m. approximated the value at 10:30 a.m. from the field measured FAPAR diurnal courses.However, such diurnal FAPAR data cannot be conveniently obtained at multiple sites in field measurements, which limits the use of this normalization method in future applications.The diurnal FAPAR variations are attributed to the changing conditions of incident solar radiations, including the solar zenith angle (SZA) and the diffuse skylight proportions [20].Cao et al. [7] attempted to characterize the effects of SZA and the fraction of skylight on diurnal FAPAR by fitting the in-situ FAPAR at different SZAs.The fitted result represented the diurnal variation of FAPAR simply, which could be the basis for temporal normalization.However, in real application, difficulties in field measurements are also major limits.Huemmrich et al. [28] suggested that the effects of SZA on FAPAR measurements should be directly addressed in future studies.However, to date, few studies have addressed the issue in temporal normalization based on easily available in-situ data.
Despite the fact that coarse-resolution FAPAR products are now available daily, field measurements do not always coincide with the acquisition dates of satellite data [29].Fensholt et al. and Martinez et al. [13,16] interpolated satellite products to the acquisition dates of ground measurements linearly, while Claverie et al.
[29] employed a cubic function to interpolate the in-situ data to the dates of satellite data.The linear interpolation adopts a simple linear assumption for crop growth, yet it can hardly capture the smooth temporal shape at the fast-growing stage of vegetation.The cubic interpolation could derive smooth and consistent temporal shape at high measurement frequency, but might produce unfeasible results when the available measurements are limited or random errors occur in the limited measurements.A more precise representation of crop growth is thus necessary for temporal interpolation.
In brief summary, scaling of FAPAR from the field to the satellite confronts spatial and temporal mismatches.While spatial mismatch has been widely concerned and studied, the temporal mismatches received far less attention.Further investigation is urged to strengthen the spatial and temporal comparability between in-situ and satellite FAPAR in the scaling.In this study, the aim is to develop a general methodology for scaling the in-situ FAPAR to the coarse-resolution satellite product.To achieve the objective, the spatial scaling framework by Morisette et al. [25] and the general scale transfer framework in [30,31] were combined.The spatial and temporal mismatches were addressed in three steps.Firstly, an FAPAR analytical model was integrated with an atmospheric radiative transfer model to simulate FAPAR with different conditions accurately.Then the coefficients for temporal normalization can be derived by fitting those simulated data with a simplified semi-experienced model.Secondly, the logistic model was used to interpolate the ground measurements to match the dates of satellite data.Thirdly, the reference FAPAR values at 15 m were estimated from fine-resolution satellite data to bridge the spatial discrepancies between in-situ data and coarse-resolution satellite product.The reference FAPAR values at coarse resolution were finally generated through spatial aggregation as spatially and temporally comparable to coarse-resolution FAPAR products.
The study is performed in Heihe River Basin, China, where arid and semi-arid climate dictates.It is a part of the Heihe Watershed Allied Telemetry Experiment Research (HiWATER), which was implemented in 2012 and collected a large volume of field-measured biophysical parameters [32].HiWATER provides valuable and reliable field FAPAR measurements for this study.The MODIS FAPAR product MOD15A2 C5 [33] was taken as a study case for coarse-resolution satellite data for its popularity in eco-environmental applications.The definition of FAPAR within MODIS product includes absorption of PAR from both direct solar radiation and diffuse skylight, which is in accordance with the ground measurements and is a particularly typical case for temporal normalization.The ground FAPAR measurements were scaled to the coarse-resolution satellite product following the proposed scaling methodology.The derived coarse-resolution FAPAR reference data were used to validate the MODIS product.Finally the performances of the scaling methodology were evaluated quantitatively.This study provides valuable reference for understanding how spatial and temporal mismatches occur and how to address these issues in scaling FAPAR, a variable with high spatial and temporal variations, from the field to the satellite.

Study Area
The study area is located in Wuxing Village, Zhangye City, China, which is at the midstream of the Heihe River Basin (Figure 1).The Heihe River is the second longest inland river going through distinct landscapes from mountain cryosphere at the upper stream, the artificial oasis at the midstream and the natural oasis at the downstream.The study area is one focus experimental area (FEA) at the midstream that belongs to one key experimental area (KEAs) for intensive and long-term observations in HiWATER [32].The study area features a compound of croplands (72%), residential land (24%) and woodland (4%).Major crops are maize and small patches of orchards and vegetables.observations in HiWATER [32].The study area features a compound of croplands (72%), residential land (24%) and woodland (4%).Major crops are maize and small patches of orchards and vegetables.

Measurement and Data Preparation
The data used in the study are summarized in Table 1, including in-situ measurement, coarseresolution FAPAR products, fine-resolution satellite imagery (ASTER-Advanced Spaceborne Thermal Emission and Reflection Radiometer) and other ancillary data.All data can be openly accessed from the HiWater website at http://westdc.westgis.ac.cn/data/.

In-situ Measurement
To guarantee their representativeness to the study area, the in-situ sites were configured following a mean of surfaces with non-homogeneity (MSN) sampling strategy, which considers the autocorrelation and stratified non-homogeneity of the area using a fine-resolution classification map [34,35].Thirteen cornfield sites were distributed throughout the area almost uniformly (Figure 1) within an extent of 5 × 5 km.From 24 May to 10 July 2012, canopy structural parameters and FAPAR

Measurement and Data Preparation
The data used in the study are summarized in Table 1, including in-situ measurement, coarse-resolution FAPAR products, fine-resolution satellite imagery (ASTER-Advanced Spaceborne Thermal Emission and Reflection Radiometer) and other ancillary data.All data can be openly accessed from the HiWater website at http://westdc.westgis.ac.cn/data/.To guarantee their representativeness to the study area, the in-situ sites were configured following a mean of surfaces with non-homogeneity (MSN) sampling strategy, which considers the autocorrelation and stratified non-homogeneity of the area using a fine-resolution classification map [34,35].Thirteen cornfield sites were distributed throughout the area almost uniformly (Figure 1) within an extent of 5 ˆ5 km.From 24 May to 10 July 2012, canopy structural parameters and FAPAR values were measured repeatedly every five days.The field measurement has been interrupted from 18 June to 3 July by bad weather and irrigation schedule, leaving a data gap during the period.Canopy structural parameters include the average canopy height, leaf numbers, plant density and LAI etc.The canopy FAPAR values were measured by the AccuPAR LP-80 Ceptometer (http://www.decagon.com)and four values were acquired at each measurement: the incident PAR above the canopy (PAR Ó a ), reflected PAR above the canopy (PAR Ò a ), incident PAR below the canopy (PAR Ó b ) and reflected PAR by ground (PAR Ò b ).According to the device instructions [2,36], the canopy FAPAR can be calculated by At each site, a 10 m ˆ10 m sample plot was enclosed and measurements were made 8 to 10 times along the diagonal of the sample plot repeatedly and continuously in a short time.To reduce the measurement errors affected by the spatial representativeness of instrument, operator-dependent errors and changing weather conditions (such as cloud or wind), an averaged FAPAR was calculated from the repeated measurements as the value representative for the plot.Details of the in-situ dataset could be referred to [37].

Satellite Data and Preprocessing
MODIS FAPAR product is among the most popular biophysical products in environmental studies [38].It is generated through a look-up-table algorithm to retrieve FAPAR parameter from a three-dimensional radiative transfer model [8].Numerous efforts have been exerted on validating MODIS FAPAR product (e.g., [13,27]) or comparing it with other global products (e.g., [5,29]).The selected MODIS FAPAR products in Table 1 were extracted, scaled and rectified to UTM projection, WGS-84 coordinate system using MODIS re-projection tool (MRT).To facilitate the multi-scale analysis starting from ASTER data at 15 m resolution, the MODIS data were resampled to 960 m resolution using nearest neighborhood algorithm.
ASTER L1B data at 15 m resolution provides the bridge between field measurements with small support area and coarse-resolution data with large pixels.By extending the field measurements to the ASTER data at 15 m resolution, the ASTER FAPAR can be aggregated to the resolution of 960 m that are spatially comparable to MODIS data.A series of processing were performed to ASTER L1B images sequentially, including radiometric calibration, atmospheric correction and geometric correction.The software of the second simulation of a satellite signal in the solar spectrum (6S) [39] and MODIS aerosol products were used for atmospheric correction, as well as solar zenith angle and azimuth angle data in ASTER metadata file.
Geometric registrations were performed on ASTER data and further MODIS product, to ensure the geometric comparability among the multi-scale satellite data.The SPOT-5 image acted as the base map, while the classification data at 1 m resolution [40] was used to testify the geometric co-registration of the multi-source data.By selecting control points, e.g., crossovers of roads and rivers, the registration accuracy of ASTER data was carefully controlled within 2 m, while that of MODIS products was within 15 m.

Methodology
Spatial and temporal mismatches are the major challenges in scaling FAPAR data from the in-situ measurements to a coarse resolution satellite product.The general scaling methodology in [30,31] provides a conceptual framework to achieve the goal, via a combined use of prior knowledge and fine-resolution data.The validation framework by Morisette et al. [25] holds a similar belief that the use of high resolution data could reduce the spatial scale discrepancy between in-situ and coarse-resolution data.Borrowing the ideas from the two methods, this paper proposed a spatial and temporal scaling methodology for in-situ FAPAR, as shown in Figure 2. The temporal and spatial mismatches were addressed in three steps as temporal normalization, temporal interpolation and spatial aggregation.
After serious control over the spatial and temporal comparability between the pair of in-situ and satellite data, the derived coarse-resolution FAPAR reference data were used to validate the MODIS FAPAR product.
Remote Sens. 2016, 8, 310 6 of 23 pair of in-situ and satellite data, the derived coarse-resolution FAPAR reference data were used to validate the MODIS FAPAR product.

FAPAR Definition
FAPAR is highly temporally variable and is easily affected by illumination.Thus FAPAR has to be defined at specific time step under specific illumination conditions.The use of FAPAR in dynamic vegetation models could run at a time step smaller than a day or at daily time step, which requires instantaneous FAPAR values or daily integrated FAPAR values [3].According to illumination conditions, FAPAR could be defined as the black-sky FAPAR (FAPAR from direct solar radiation) and the white-sky FAPAR (FAPAR from diffuse solar radiation) [3].The instantaneous FAPAR is always a weighted sum of black-sky FAPAR and white-sky FAPAR, where the weight is related to the proportions of skylight.Depending on the inclusion of PAR absorbed by either green (live leaves) or all (including woody materials and dead leaves) part of vegetation, FAPAR could be discriminated as green FAPAR and total FAPAR [41].For MODIS FAPAR product, FAPAR is defined as the instantaneous FAPAR at the time of satellite overpass [3] using direct and diffuse incoming radiation [41].The more detailed overview of product definitions could be referred to [5,41].
In-situ FAPAR values are generally measured by radiation balance monitoring [13,42] or by measuring directional gap fraction with hemispherical photos [43].In the latter case, the measured FAPAR is approximated by the light interception efficiency and only includes absorption of direct solar radiation.
The study attempts to validate the MODIS FAPAR product using ground measurements by AccuPAR device (the radiation balance monitoring approach).Essentially the definition of FAPAR between the satellite product and ground measurements agrees, which lies the basis of the study.

Temporal Normalization
Limited by time, labor and available instruments, mobile measurements were performed in a sequential manner across multiple field sites.The problem naturally arose as how to match the FAPAR values measured at varying hours to the time of satellite overpass.Since the effects of subpixel cloud could be neglected in fine-resolution ASTER data, the time of in-situ measurement

FAPAR Definition
FAPAR is highly temporally variable and is easily affected by illumination.Thus FAPAR has to be defined at specific time step under specific illumination conditions.The use of FAPAR in dynamic vegetation models could run at a time step smaller than a day or at daily time step, which requires instantaneous FAPAR values or daily integrated FAPAR values [3].According to illumination conditions, FAPAR could be defined as the black-sky FAPAR (FAPAR from direct solar radiation) and the white-sky FAPAR (FAPAR from diffuse solar radiation) [3].The instantaneous FAPAR is always a weighted sum of black-sky FAPAR and white-sky FAPAR, where the weight is related to the proportions of skylight.Depending on the inclusion of PAR absorbed by either green (live leaves) or all (including woody materials and dead leaves) part of vegetation, FAPAR could be discriminated as green FAPAR and total FAPAR [41].For MODIS FAPAR product, FAPAR is defined as the instantaneous FAPAR at the time of satellite overpass [3] using direct and diffuse incoming radiation [41].The more detailed overview of product definitions could be referred to [5,41].
In-situ FAPAR values are generally measured by radiation balance monitoring [13,42] or by measuring directional gap fraction with hemispherical photos [43].In the latter case, the measured FAPAR is approximated by the light interception efficiency and only includes absorption of direct solar radiation.
The study attempts to validate the MODIS FAPAR product using ground measurements by AccuPAR device (the radiation balance monitoring approach).Essentially the definition of FAPAR between the satellite product and ground measurements agrees, which lies the basis of the study.

Temporal Normalization
Limited by time, labor and available instruments, mobile measurements were performed in a sequential manner across multiple field sites.The problem naturally arose as how to match the FAPAR values measured at varying hours to the time of satellite overpass.Since the effects of subpixel cloud could be neglected in fine-resolution ASTER data, the time of in-situ measurement should ideally be normalized to 10:30 a.m.local time [13].Thus, the aim of temporal normalization is to normalize the FAPAR measured at different time to the time of satellite overpass.
For the purpose of temporal normalization, a simple equation that can characterize the FAPAR measured at a specific time (at a given SZA) as a function of canopy structures, SZA and illumination conditions is preferred.Based on this, temporal normalization can be achieved by deriving the relationship between the instantaneous FAPAR at different time.The section is organized into: (i) the introduction to the FAPAR-P model (Section 3.2.1);(ii) the estimation of skylight proportion from the 6S model (Section 3.2.2);(iii) the simplified FAPAR model integrating FAPAR-P and 6S models (Section 3.2.3);and (iv) the derivation of the simplified SZA correction model (Section 3.2.4).

FAPAR-P Model
Physical models characterize the effects of SZA and skylight proportions on the diurnal variations of instantaneous FAPAR and are thus the basis for normalizing in-situ FAPAR measured at varying time.
Previous studies [6-8] clearly show the effects of skylight proportion and SZA on the instantaneous FAPAR.However, to correct the effects, a simple and easy-to-retrieve physical model is preferred.
The FAPAR-P model developed by Fan et al. [6] was used in the study as the theoretical basis for temporal normalization.It considers the multiple scattering of incident photons in the canopy by introducing the concept of recollision probability, a parameter only dependent on canopy structure as defined in [44].Based on energy conservation law, the fraction of absorption is computed as the probability of intercepting photons minus the total probability of scattering.By using a complete consideration of canopy structure, incident solar radiation and spectral characteristics, the model is capable of capturing the temporal variations at different time.By using effective LAI instead of LAI, the FAPAR-P model can be extended to row crops (e.g., maize) [6].The relevant input parameters and variables of the FAPAR-P model were summarized in Table 2, as well as the data sources to run the model.The leaf reflectance and transmittance were acquired from LOPEX database corresponding to maize, while the soil reflectance were measured with an ASD spectrometer in the field.The basic assumption of G = 0.5 was inherited from the development of the FAPAR-P model and has been verified by Monte Carlo simulation.More details can be found in [6].According to the study of Fan et al.
[6], the diurnal variations of FAPAR were mainly attributed to LAI, SZA and skylight proportions based on the sensitivity analysis of the FAPAR.The other parameters in Table 2, including the spectra of soil, leaf reflectance and transmittance, as well as p and G can be given as constants under the normal condition of the given scene, which can simplify the problem of temporal normalization of FAPAR.Thus, the use of the FAPAR-P model is to investigate the effects of LAI, θ i and β on FAPAR, while θ i and β mainly account for the temporal discrepancies among the in-situ measurement at varying time.

Estimation of Skylight Proportion
Based on the FAPAR-P model, the most relevant factors accounting for the diurnal variations of FAPAR are SZA and skylight proportions.The proportion of skylight is also highly variable with SZA, as well as aerosol optical depth (AOD), visibility (VIS) and cloud cover.As real-time measurement of skylight proportions are not existed, it is difficult to determine the value of skylight proportion.The major task becomes how to determine the proportion of skylight from relevant environmental factors.
The 6S atmospheric radiative transfer model is employed to model the skylight proportion.Occasional cloud was not considered here, as it is random and hard to measure, which would further complicate the modeling.For a small study area stretching several kilometers, the AOD was regarded as constant spatially and temporally, while the VIS was assumed to be constant in space at the same day.Three VIS cases, including mildly hazy, moderately clear and clear air conditions, were considered, corresponding to VIS = 5 km, 15 km and 30 km in 6S.The atmospheric condition and aerosol model in 6S were set to standard midlatitude summer and continental type, respectively.By running 6S with three VIS cases and varying SZAs, the relationships between the skylight proportions and SZAs were derived and were further fitted using exponential equations (Figure 3), expressed as: which corresponds to VIS = 5 km, 15 km and 30 km, sequentially.The goodness of the fits were evaluated in root mean squared error (RMSE) and the correlation coefficient (R 2 ), where RMSE = 0.007, 0.014, and 0.013 and R 2 = 0.999, 0.996, and 0.996 for Equations ( 2)-( 4), respectively.

Estimation of Skylight Proportion
Based on the FAPAR-P model, the most relevant factors accounting for the diurnal variations of FAPAR are SZA and skylight proportions.The proportion of skylight is also highly variable with SZA, as well as aerosol optical depth (AOD), visibility (VIS) and cloud cover.As real-time measurement of skylight proportions are not existed, it is difficult to determine the value of skylight proportion.The major task becomes how to determine the proportion of skylight from relevant environmental factors.
The 6S atmospheric radiative transfer model is employed to model the skylight proportion.Occasional cloud was not considered here, as it is random and hard to measure, which would further complicate the modeling.For a small study area stretching several kilometers, the AOD was regarded as constant spatially and temporally, while the VIS was assumed to be constant in space at the same day.Three VIS cases, including mildly hazy, moderately clear and clear air conditions, were considered, corresponding to VIS = 5 km, 15 km and 30 km in 6S.The atmospheric condition and aerosol model in 6S were set to standard midlatitude summer and continental type, respectively.By running 6S with three VIS cases and varying SZAs, the relationships between the skylight proportions and SZAs were derived and were further fitted using exponential equations (Figure 3), expressed as: which corresponds to VIS = 5 km, 15 km and 30 km, sequentially.The goodness of the fits were evaluated in root mean squared error (RMSE) and the correlation coefficient (R 2 ), where RMSE = 0.007, 0.014, and 0.013 and R 2 = 0.999, 0.996, and 0.996 for Equations ( 2)-( 4), respectively.

FAPAR Model Integrating FAPAR-P and 6S Models
Using Equations ( 2)-( 4), the skylight proportions estimated from SZA and VIS can be substituted into the FAPAR-P model.Thus the effects of skylight proportions on instantaneous FAPAR were integrated into the effects of SZA.If only the skylight proportion determined by SZA and VIS were considered, the diurnal variations of FAPAR for a given LAI were induced by changing SZAs.

FAPAR Model Integrating FAPAR-P and 6S Models
Using Equations ( 2)-( 4), the skylight proportions estimated from SZA and VIS can be substituted into the FAPAR-P model.Thus the effects of skylight proportions on instantaneous FAPAR were integrated into the effects of SZA.If only the skylight proportion determined by SZA and VIS were considered, the diurnal variations of FAPAR for a given LAI were induced by changing SZAs.
Using the integrated FAPAR model, the effects of SZA and the associated skylight proportions on instantaneous FAPAR were firstly investigated (Figure 4).As shown in Figure 4a, larger SZA increases FAPAR values, especially at low LAI values.At LAI = 1, the increase in FAPAR from 0 ˝to 70 ˝can reach as high as 50%.This effect could be explained by the lengthening of the photons' paths going through the canopy which increase their probabilities of being absorbed or scattered.The higher the value of LAI, the less is the effect of SZA on FAPAR.At extremely high LAI, i.e., LAI = 8, FAPAR is very close to 1 and the effect of SZA is minimal.In Figure 4b, the FAPAR increases with LAI increasing, exhibiting an exponential shape.The increase in FAPAR with large SZA could also be observed at the same LAI.Using the integrated FAPAR model, the effects of SZA and the associated skylight proportions on instantaneous FAPAR were firstly investigated (Figure 4).As shown in Figure 4a, larger SZA increases FAPAR values, especially at low LAI values.At LAI = 1, the increase in FAPAR from 0° to 70° can reach as high as 50%.This effect could be explained by the lengthening of the photons' paths going through the canopy which increase their probabilities of being absorbed or scattered.The higher the value of LAI, the less is the effect of SZA on FAPAR.At extremely high LAI, i.e., LAI = 8, FAPAR is very close to 1 and the effect of SZA is minimal.In Figure 4b, the FAPAR increases with LAI increasing, exhibiting an exponential shape.The increase in FAPAR with large SZA could also be observed at the same LAI.Using the integrated FAPAR model, the effects of SZA and the associated skylight proportions on instantaneous FAPAR were firstly investigated (Figure 4).As shown in Figure 4a, larger SZA increases FAPAR values, especially at low LAI values.At LAI = 1, the increase in FAPAR from 0° to 70° can reach as high as 50%.This effect could be explained by the lengthening of the photons' paths going through the canopy which increase their probabilities of being absorbed or scattered.The higher the value of LAI, the less is the effect of SZA on FAPAR.At extremely high LAI, i.e., LAI = 8, FAPAR is very close to 1 and the effect of SZA is minimal.In Figure 4b, the FAPAR increases with LAI increasing, exhibiting an exponential shape.The increase in FAPAR with large SZA could also be observed at the same LAI.To correct the effects of SZA, the key is the slope of FAPAR varying with independent variable LAI/cos(θ).The steeper the slope is, the more the correction is needed.It could be inferred that such slope is much steeper at lower LAI values; given the same LAI, the slope decreases with the increasing SZAs.Thus the effects of SZA are intertwined with both LAI and VIS.However, for each case with given LAI and VIS, the FAPAR can be fitted by an exponential formula, as: FAPAR " k 1 ´k2 ˆexp p´G ˆLAI{cos pθq ˆk3 q (5) where k 1 , k 2 and k 3 are fitting parameters.The parameters of k 1 and k 2 are associated with the composition of direct solar radiation and diffuse skylight, as well as the multi-scattering within the canopy layer.The parameter k 3 is related to the vegetation structure.It is worth noting that when k 1 , k 2 and k 3 all equal to 1, Equation ( 5) becomes FAPAR " 1 ´e´G¨LAI{cospθq , as shown in the solid purple line in Figure 5.This is in accordance with the form of the FAPAR model based on gap probability [45].
The equation only stands under the assumption of "black leaves" (r l pλq " τ l pλq " 0) that makes p equal to zero and only the direct solar radiation is considered.At SZA lower than 40 ˝, simple SZA correction could be done based on this simpler form of FAPAR model.The precise correction can be achieved by fitting each dashed line in Figure 5 by Equation ( 5), based on which the relationship between FAPAR measured at different time (SZAs) can be derived specifically.In this way a look up table of k 1 , k 2 and k 3 with varying LAI, SZA and VIS values could be established, which is precise but is complex.The RMSE of fitting FAPAR with LAI and SZA by k 1 , k 2 and k 3 were ranging 0.0002~0.0017,which were very low and were thus negligible.
The changes of k 1 , k 2 and k 3 with LAI values followed the specific pattern, such as an exponential form for k 1 .If the relationship of LAI and coefficients k could be modeled further, the correction can be performed easily by deriving k 1 and k 2 from the simultaneously measured LAI values.The fitted models of k 1 and k 2 with LAI followed the pattern as k 1 = a 1 ¨exp(-a 2 ¨LAI)+a 3 and k 2 = a 1 ¨exp(a 2 ¨LAI)+ a 3 ¨exp(a 4 ¨LAI).The fitting of k 1 and k 2 to LAI values were also very accurate, with RMSE ranging from 0.0002 to 0.0100.

Simplified SZA Correction Model
Given the instantaneous FAPAR measured at the time of satellite overpass FAPAR pθ 0 q, the instantaneous FAPAR at any other time is FAPAR pθ i q, then FAPAR pθ i q could be corrected to FAPAR pθ 0 q following FAPAR pθ 0 q " k 1 ´" where the correction is related to k 1 and k 2 , while the parameter k 3 has been eliminated by comparing FAPAR measurements at two SZAs.
To diminish the effects of random measurement errors, two or four instantaneous FAPAR values measured around 10:30 a.m.(θ 0 ) symmetrically were selected and corrected, and the average corrected value was taken as the value of FAPAR pθ 0 q.

Interpolating in-Situ Measurements to ASTER Dates
Previous studies tended to apply linear interpolation to either satellite data or field data to achieve the temporal consistency between them, as conducted by Martinez et al. [16] and Mu et al. [35].This is feasible when the apart dates between the observation of in-situ data or satellite data are short or the parameters (such as fractional vegetation cover) are assumed to be steady in a short time.Claverie et al.
[29] interpolated the in-situ measurements to the dates of satellite observations by a cubic model, which derived a smooth and consistent temporal shape.However, the cubic interpolation is sensitive to the amount of samples and the outliers.Random errors in limited in-situ data might result unreasonable temporal curves.As mentioned above, a data gap from 18 June to 3 July existed in the temporal profiles of the in-situ measurements, overlapping with the fast growing period of vegetation.Either linear or cubic interpolation might not be capable of capturing the temporal curves of crop growth precisely.More precise representation of FAPAR time series is thus required.Zhang et al. [46] proposed a method to represent vegetation phenology as series of piecewise logistic functions of time.Mu et al. [47] applied the logistic model to fit the time series of vegetation index.The logistic representation is advantageous, because the growth speed at different stages can be described accurately in a simple form.Therefore, it is employed in the study for temporal interpolation, expressed as: where b 1 and b 2 are fitting parameters; b 3 `b4 is the maximum FAPAR value and b 4 is the initial background FAPAR value; and FAPAR t is the FAPAR value at the observation date t.The parameter b 4 is set to 0 for bare soil background; the value of b 3 `b4 is set to 1.0 for complete absorption of PAR by vegetation.In implementation, only sites with more than three dates of observations were considered to ensure the reliability of the fitted equations.The time series of thirteen field sites were fitted using the logistic model respectively and then temporal interpolation was performed.

Interpolating MODIS FAPAR to ASTER Dates
The logic behind the 8-day compositing of MODIS FAPAR product is to select the maximum FAPAR across the 8 days for a given pixel [48].For simplicity, it is assumed that the maximum FAPAR value in the 8-day interval corresponds to the end day of the composite period, rather than the start day.In the study, the acquisition dates of ASTER data are very close to the end day of the 8-day intervals of MODIS FAPAR product.Thus, a linear interpolation approach was regarded as acceptable to interpolate the coarse-resolution FAPAR products to the dates of ASTER data, which is suggested by Mu et al. [47] when the intervals were short.

Multi-Scale FAPAR Truth Data
As the support scale of in-situ measurement is limited, the study used fine-resolution ASTER data to expand the "point" scale of field measurement to coarse-resolution "pixel" scale.As the homogeneities of the areas surrounding the in-situ sites were carefully controlled, the pixel at 15 m resolution could be representative of each site.The derived ASTER FAPAR at fine resolution and the aggregated ASTER FAPAR at coarse resolution were thus taken as the truth reference data at two scales.

Estimating Fine-Resolution FAPAR
The normalized difference vegetation index (NDVI) is widely believed to be significantly related to FAPAR [8,49].The linear or approximately linear relationship between NDVI and FAPAR has been testified by satellite data [13,50], radiative transfer models [20,[51][52][53] and field data [13].The NDVI data were calculated from ASTER reflectance.An empirical regression approach was employed to transfer ASTER NDVI to FAPAR, as conducted in [26], in the form as: where c 1 and c 2 are the coefficients converting NDVI to FAPAR.To reduce the effects of outliers, a robust regression using iteratively reweighted least squares (IRLS) [54,55] was applied.It is used to find the optimized estimates of a linear model by an iterative method in which each step involves solving a weighted least squares problem.Compared to the ordinary least square regression, the regression using IRLS is advantageous to mitigate the influence of outliers that derive large residual errors and are assigned lower weights in the iteration.After the fit, the derived regression coefficients were used to estimate FAPAR values from ASTER NDVI.

Coarse-Resolution FAPAR
ASTER satellite data were geometrically registered using a SPOT image at 10 m resolution.Coarse-resolution (MODIS) FAPAR products were then registered to match the ASTER data.
Based on a careful geometric registration between ASTER and coarse-resolution FAPAR products, fine-resolution FAPAR data were aggregated to coarse resolutions to form the truth reference FAPAR data.The MODIS FAPAR product and the coarse-resolution reference FAPAR data were compared pixel by pixel.

Sensitivities of FAPAR to Extra Skylight Proportions
In the proposed temporal normalization and analysis, only the proportion of skylight associated with SZA was However, extra changes in skylight proportion could also be integrated into FAPAR-P model, which is interesting for further considerations on temporal correction issues.
Two cases of SZA = 0 ˝and SZA = 60 ˝were conducted.Thus the minimum of skylight proportion is the proportion solely determined by SZA, while the maximum is 100%, as in Figure 6.It could be inferred that the effects of extra skylight are nearly linear and are more obvious at low LAI values.Besides, the changes of FAPAR with extra skylight show increasing trend at SZA = 0 ˝, and a mild decreasing trend at SZA = 60 ˝.In both cases the FAPAR value at 100% skylight remains the same.The opposite trend could be attributed to the more significant effects of SZA than those of skylight proportion.The findings were in accordance with those revealed in [7], where a higher skylight proportion tends to increase FAPAR when SZA < 21 ˝and to reduce it when SZA > 30 ASTER satellite data were geometrically registered using a SPOT image at 10 m resolution.Coarse-resolution (MODIS) FAPAR products were then registered to match the ASTER data.
Based on a careful geometric registration between ASTER and coarse-resolution FAPAR products, fine-resolution FAPAR data were aggregated to coarse resolutions to form the truth reference FAPAR data.The MODIS FAPAR product and the coarse-resolution reference FAPAR data were compared pixel by pixel.

Sensitivities of FAPAR to Extra Skylight Proportions
In the proposed temporal normalization and analysis, only the proportion of skylight associated with SZA was considered.However, extra changes in skylight proportion could also be integrated into FAPAR-P model, which is interesting for further considerations on temporal correction issues.
Two cases of SZA = 0° and SZA = 60° were conducted.Thus the minimum of skylight proportion is the proportion solely determined by SZA, while the maximum is 100%, as in Figure 6.It could be inferred that the effects of extra skylight are nearly linear and are more obvious at low LAI values.Besides, the changes of FAPAR with extra skylight show increasing trend at SZA = 0°, and a mild decreasing trend at SZA = 60°.In both cases the FAPAR value at 100% skylight remains the same.The opposite trend could be attributed to the more significant effects of SZA than those of skylight proportion.The findings were in accordance with those revealed in [7], where a higher skylight proportion tends to increase FAPAR when SZA < 21° and to reduce it when SZA > 30°.As shown in Figure 6, given specific SZA a change of 20% of skylight proportion affects FAPAR in a subtle way, even at low LAI values.The effect of SZA is proved to be more influential than the extra proportion of skylight.

Validating Temporal Normalization Method
To validate the proposed temporal normalization method, a theoretical analysis of the accuracy of the model and a validation using real in-situ FAPAR were conducted.
From a theoretical perspective, the errors of the temporal normalization were identified to be sourced from: (i) the modeling errors from the FAPAR-P and the 6S models; (ii) the residual errors of fitting the modeled FAPAR in exponential forms and of fitting the correction coefficients to LAI; and As shown in Figure 6, given specific SZA a change of 20% of skylight proportion affects FAPAR in a subtle way, even at low LAI values.The effect of SZA is proved to be more influential than the extra proportion of skylight.

Validating Temporal Normalization Method
To validate the proposed temporal normalization method, a theoretical analysis of the accuracy of the model and a validation using real in-situ FAPAR were conducted.
From a theoretical perspective, the errors of the temporal normalization were identified to be sourced from: (i) the modeling errors from the FAPAR-P and the 6S models; (ii) the residual errors of fitting the modeled FAPAR in exponential forms and of fitting the correction coefficients to LAI; and (iii) the measurement errors of input data (e.g., LAI and FAPAR).Firstly, the modeling errors from the FAPAR-P and the 6S model were found to be both within 1% in [6,56], respectively.As a 1% change in extra skylight proportions has shown negligible effects on FAPAR (Figure 6), the overall modeling errors of integrating 6S and FAPAR-P model was estimated as within 1%.Secondly, the derivation of the temporal normalization method used two-round fittings: the fitting of the FAPAR to k 1 , k 2 and k 3 using Equation ( 5) and the fitting of k 1 and k 2 to LAI under varying VIS cases.The RMSE values for the two-round fittings ranged 0.0002~0.0017and 0.0002~0.0100,respectively, which were testified as very small.To further investigate the effects of the residual errors from the fitting of the integrated FAPAR model, the FAPAR values derived from k 1 and k 2 estimated by LAI (from the two-round fittings) were compared to those fitted FAPAR values (from the first-round fitting) and those directly derived from the integrated model.The comparisons were made for varying LAI values at varying SZAs given VIS = 30 km.RMSEs were calculated for varying SZAs from 0 ˝to 60 ˝for each LAI value, as shown in Table 3.It was found that the errors from the second-round fitting of k 1 and k 2 were less than 0.008.The combined errors of the two-round fittings were also less than 0.008.This validated that the propagation of errors from the fitting of k 1 and k 2 in estimating FAPAR was negligible.The overall errors from the two-round fittings were within 0.01.Thirdly, the effects of measurement errors on temporal normalization were investigated by a sensitivity analysis.The measurement errors were mainly sourced from the measured LAI that determines the correction coefficients k 1 and k 2 and the FAPAR at specific time that is to be corrected.The sensitivity analysis was performed by adding random errors to both LAI and FAPAR measurements and then correcting FAPAR values at varying SZAs to 0 ˝.In implementation, the input LAI and FAPAR values were generated randomly with normal distribution, defined in the form of mean ˘standard deviation (STD) (blue error bars in Figure 7).The STD is an indicator for uncertainty [57].For both LAI and FAPAR, the STD was set to about 20% of their mean values.The mean LAI values were set to 1, 2, 3 and 4, while the mean FAPAR values were set to the true FAPAR values at varying SZAs derived from the integrated model given VIS = 30 km.As shown in Figure 7, the mean values of the corrected FAPAR values (red error bars) were very close to the true FAPAR values at SZA = 0 (grey dashed line), with small values of STD.The correction errors were found to be related to and were less than the errors in the input FAPAR.Significant reduction of uncertainties has been observed in each step.In estimating k 1 and k 2 from input LAI, the uncertainties of the result k 1 and k 2 were less than 10%, while the SZA correction using the resulted k 1 , k 2 and input FAPAR values lead to further reductions of uncertainty, less than 0.02 under all cases.This clearly verified that both errors and uncertainties in the measurement were reduced in the temporal normalization.
Remote Sens. 2016, 8, 310 14 of 23 Figure 8).The LAI was 2.56, from which the diurnal cycle was modeled by the integrated model at VIS = 30 km (blue dotted line in Figure 8).The temporally smoothed in-situ FAPAR were highly consistent with the simulated diurnal cycle.The temporal normalization method was applied to correct the temporally smoothed in-situ data to 10:30 a.m.(red asterisk in Figure 8).The VIS was estimated as clear condition (30 km) by visual inspection.The mean of the corrected FAPAR values was 0.805 with STD = 0.009.The average FAPAR A series of real in-situ data from one day with continuous observations were used to validate the proposed temporal normalization method.The data were obtained at 10 min interval on 5 July 2012.The diurnal variations measured at the fixed position of the field showed dramatic irregular fluctuations (gray square in Figure 8) because of the small representative area of the measurement device and the heterogeneity within canopy.Thus, we smoothed the data by a simple moving-window averaging algorithm with 100 min window to reduce the random turbulence (red line in Figure 8).The LAI was 2.56, from which the diurnal cycle was modeled by the integrated model at VIS = 30 km (blue dotted line in Figure 8).The temporally smoothed in-situ FAPAR were highly consistent with the simulated diurnal cycle.
The temporal normalization method was applied to correct the temporally smoothed in-situ data to 10:30 a.m.(red asterisk in Figure 8).The VIS was estimated as clear condition (30 km) by visual inspection.The mean of the corrected FAPAR values was 0.805 with STD = 0.009.The average FAPAR measured between 10:00 a.m. and 11:00 a.m. was 0.792, taken as the true FAPAR value at 10:30 a.m.The mean error is thus 0.013 with an uncertainty of 0.009.The error and uncertainty from the real in-situ data were consistent with the theoretical analysis of the model accuracy.
The uncertainty analysis based on the simulated and the measured data both validated that the proposed temporal normalization method was highly accurate and reliable.The errors from the modeling and the two-round fittings in all were below 0.02 theoretically.The temporal normalization method was applied to correct the temporally smoothed in-situ data to 10:30 a.m.(red asterisk in Figure 8).The VIS was estimated as clear condition (30 km) by visual inspection.The mean of the corrected FAPAR values was 0.805 with STD = 0.009.The average FAPAR

Comparing Temporal Interpolation Methods
The in-situ data after temporal normalization were fitted using a logistic model respectively, which is further applied to temporal interpolation.The derived parameters of b 1 and b 2 in Equation (7) ranged 12.20~16.52and ´0.097~´0.070,individually.The mean of the residual errors of logistic fitting were 0.003, with STD = 0.065.To illustrate the feasibility of the logistic model, the linear and cubic fittings were also performed for comparison.The R 2 and RMSE were calculated to measure the goodness of fit.
As shown in Figure 9, the logistic fit derived the higher R 2 values and lower RMSE values than the cubic fit does, in 76.9% and 69.2% of the in-situ sites respectively.Compared to the linear model, the logistic model performed significantly better evaluated in both R 2 and RMSE at all sites.For logistic, linear and cubic models, the mean of R 2 values was calculated as 0.978, 0.937 and 0.969 respectively, while that of RMSE was 0.070, 0.090 and 0.075.Although the values of the cubic model were very close to those of the logistic model, the logistic model was much more stable, with R 2 and RMSE ranging 0.951~0.994and 0.040~0.107,while the R 2 and RMSE values of the cubic model ranged 0.914~0.996and 0.027~0.150.In short, the logistic interpolation is preferred as it maintains stable and good performance in fitting existing data series.measured between 10:00 a.m. and 11:00 a.m. was 0.792, taken as the true FAPAR value at 10:30 a.m.The mean error is thus 0.013 with an uncertainty of 0.009.The error and uncertainty from the real insitu data were consistent with the theoretical analysis of the model accuracy.
The uncertainty analysis based on the simulated and the measured data both validated that the proposed temporal normalization method was highly accurate and reliable.The errors from the modeling and the two-round fittings in all were below 0.02 theoretically.

Comparing Temporal Interpolation Methods
The in-situ data after temporal normalization were fitted using a logistic model respectively, which is further applied to temporal interpolation.The derived parameters of and in Equation (7) ranged 12.20~16.52and −0.097~−0.070,individually.The mean of the residual errors of logistic fitting were 0.003, with STD = 0.065.To illustrate the feasibility of the logistic model, the linear and cubic fittings were also performed for comparison.The R 2 and RMSE were calculated to measure the goodness of fit.
As shown in Figure 9, the logistic fit derived the higher R 2 values and lower RMSE values than the cubic fit does, in 76.9% and 69.2% of the in-situ sites respectively.Compared to the linear model, the logistic model performed significantly better evaluated in both R 2 and RMSE at all sites.For logistic, linear and cubic models, the mean of R 2 values was calculated as 0.978, 0.937 and 0.969 respectively, while that of RMSE was 0.070, 0.090 and 0.075.Although the values of the cubic model were very close to those of the logistic model, the logistic model was much more stable, with R 2 and RMSE ranging 0.951~0.994and 0.040~0.107,while the R 2 and RMSE values of the cubic model ranged 0.914~0.996and 0.027~0.150.In short, the logistic interpolation is preferred as it maintains stable and good performance in fitting existing data series.

Estimating Fine-Resolution FAPAR and Spatial Scaling
The linear regression is easily affected by the distribution of data.To identify such effects, the robust regression using IRLS algorithm was applied to the subset of data on each day and to the entire dataset sequentially, as shown in Figure 10.Clearly the fits with each subset and the entire      The histograms of MODIS FAPAR were compared to the upscaled ASTER data, as shown in Figure 12.The distribution of MODIS FAPAR values in general coincided well with the upscaled ASTER FAPAR values.Especially on 15 June, the distribution of data between the pair of reference and MODIS FAPAR were very similar.The only less coincident distribution was observed on 10 July.The upscaled ASTER FAPAR values were more concentrated around 0.75, while the MODIS FAPAR values were more uniformly distributed in the range 0.6~0.8.This might be partially attributed to the tendency of MODIS FAPAR product to derive less high estimates at coarse spatial resolution, which was also observed in [16] that MODIS values rarely exceed FAPAR = 0.6.However, in general, the spatial comparison verified the high spatial agreement between the MODIS FAPAR product and the upscaled ASTER FAPAR data.and MODIS FAPAR were very similar.The only less coincident distribution was observed on 10 July.The upscaled ASTER FAPAR values were more concentrated around 0.75, while the MODIS FAPAR values were more uniformly distributed in the range 0.6~0.8.This might be partially attributed to the tendency of MODIS FAPAR product to derive less high estimates at coarse spatial resolution, which was also observed in [16] that MODIS values rarely exceed FAPAR = 0.6.However, in general, the spatial comparison verified the high spatial agreement between the MODIS FAPAR product and the upscaled ASTER FAPAR data.The statistics of differences between the MODIS and the upscaled ASTER FAPAR also revealed a high degree of agreement (Table 5).Except on 10 July, the mean absolute and the STD of the differences were less than 0.04 and 0.015.The MODIS FAPAR values on 10 July were less consistent and concentrated, with a higher value of 0.072 and 0.034 for the mean absolute and the STD of the differences.Except on 30 May, the MODIS product showed slight underestimates compared to the reference FAPAR as the mean values were all negative from 15 June to 10 July.In general the relative errors were all very small.To further quantify the consistency of MODIS FAPAR with the upscaled ASTER FAPAR, the scatter diagrams were analyzed (Figure 13).The MODIS FAPAR values agreed well with the truth values at 1:1 line and almost all values were concentrated within the 0.1 envelope of the 1:1 line.This evidently confirmed the high consistency of MODIS FAPAR with the truth reference data.The overall RMSE of MODIS was 0.054, which reinforced the good performance of the MODIS FAPAR product.The statistics of differences between the MODIS and the upscaled ASTER FAPAR also revealed a high degree of agreement (Table 5).Except on 10 July, the mean absolute and the STD of the differences were less than 0.04 and 0.015.The MODIS FAPAR values on 10 July were less consistent and concentrated, with a higher value of 0.072 and 0.034 for the mean absolute and the STD of the differences.Except on 30 May, the MODIS product showed slight underestimates compared to the reference FAPAR as the mean values were all negative from 15 June to 10 July.In general the relative errors were all very small.To further quantify the consistency of MODIS FAPAR with the upscaled ASTER FAPAR, the scatter diagrams were analyzed (Figure 13).The MODIS FAPAR values agreed well with the truth values at 1:1 line and almost all values were concentrated within the 0.1 envelope of the 1:1 line.This evidently confirmed the high consistency of MODIS FAPAR with the truth reference data.The overall RMSE of MODIS was 0.054, which reinforced the good performance of the MODIS FAPAR product.

Temporal Validation
The time series of MODIS FAPAR at all in-situ sites were compared to the upscaled ASTER FAPAR values (Figure 14).In the up-scaling process, the STD values of FAPAR within each 64 ˆ64 ASTER pixel block corresponding to each MODIS pixel were also calculated, shown as red bars in Figure 13.The STD represented the degree of sub-pixel heterogeneities within each MODIS pixel.Generally increasing trends of STD with time were also observed.The increased sub-pixel heterogeneities within each MODIS pixel might also be a reason for the less consistency of MODIS FAPAR with the truth reference data on 10 July.On MODIS pixel scale, the No. 4 site has the lowest subpixel heterogeneities while the No. 13 site has the highest subpixel heterogeneities.Although the degrees of subpixel heterogeneities were various, the MODIS FAPAR agreed very well with the upscaled ASTER FAPAR.This is in accordance with the above analysis in spatial comparisons.

Discussion
The scaling methodology proposed in the study directly addressed the spatial and temporal mismatches between in-situ and coarse-resolution satellite data.However, there are still many issues associated with the further developments of the scaling methodology or the validation attempts.
First, as the measurements of skylight proportions were unavailable here, the skylight proportions were decomposed into the part affected by SZA and the part affected by atmosphere.The effects of cloud were not considered here for simplicity.In further studies, when measurements of skylight proportions are available, it could be used in the model directly.
Second, the leaf angle distribution of maize was assumed to be spherical (G = 0.5) and the effects of SZA on the parameter G were assumed to be constant in the development of the temporal normalization method.This assumption was verified as feasible in developing FAPAR-P model.However, the leaf angle distribution of maize is also found to be variable with seed orientation [58] and phenological stages [59].The use of more precise description of leaf angle distribution could thus be considered in other work.Similarly, for other types of vegetation (e.g., grass characterized by erectophile leaf orientation), the assumption for the parameter G should be modified in applying the temporal normalization method.

Temporal Validation
The time series of MODIS FAPAR at all in-situ sites were compared to the upscaled ASTER FAPAR values (Figure 14).In the up-scaling process, the STD values of FAPAR within each 64 × 64 Third, the atmospheric visibilities in temporal normalization were simplified into three typical cases of mildly hazy (5 km), moderately clear (15 km) and clear air conditions (30 km).It is practical because the normalization of FAPAR is not sensitive to the VIS and three VIS cases are distinguishable.However, the simplification could be a potential source of error.In further applications where accurate VIS data are available, more precise classification of VIS can be used.
Fourth, the estimated reference FAPAR data come with uncertainties from measurement errors, temporal correction, geometric errors and the linear fit.The errors inherent in the measurement were identified as a major source of error.The temporal corrections, including both temporal normalization and interpolation, introduced low systematic errors and were found to reduce the uncertainties in measurement significantly.The systematic error of the fine-resolution reference FAPAR was thus very low with RMSE = 0.077, which could be regarded as acceptable.The errors of the coarse-resolution reference FAPAR would be further decreased as the fit errors and geometric errors were decreased in spatial aggregation [35,60].Besides, the NDVI-FAPAR relationship established on maize fields were applied to the entire study area, which would also be a potential source of error.However, the error was limited as the study area was dominated by maize.In heterogeneous areas, site-specific functions of transferring NDVI to FAPAR should be established.From a general perspective, further efforts could be exerted on long lasting in-situ data collection and extending the scaling method to more land cover types.
Last, the suitability of the proposed methods, especially temporal correction methods, should be discussed.On the one hand, the differences in ground measurements of FAPAR leads to the different requirements of temporal normalization.The proposed temporal normalization method is targeted at the measurement by energy balance monitoring.On the other hand, differences in FAPAR definitions between satellite products and ground measurements should be considered.For the satellite products that include both direct and diffuse FAPAR (e.g., MODIS\MISR), the proposed methods are feasible and valid.However, for the black-sky FAPAR products (e.g., GEOV1\CYCLOPES), the light interception measurement is preferred or the absorbed PAR from direct radiation should be separated from the measurement by energy balance monitoring.

Conclusions
To address the spatial and temporal mismatches between in-situ measurements and satellite observations, this study proposed a general methodology for scaling FAPAR from the field to the satellite.The mismatches were addressed in three steps including temporal normalization, interpolation and spatial aggregation.Main conclusions were: (i) The method of temporal normalization was feasible and reliable.By integrating the atmospheric radiative transfer model (6S) and the FAPAR-P model, the effects of multiple influential factors (e.g., VIS, SZA and LAI) on instantaneous FAPAR were characterized accurately.The fitting on the model simulations made the final correction only dependent on the existed FAPAR measurement, SZA and LAI, which is easy to perform.The accuracy of the temporal normalization method was validated using a set of continuous in-situ measurement from one day as 0.013 (RMSE); (ii) The logistic model for temporal interpolation was reliable and accurate.It maintained higher goodness and robustness of fit on in-situ data than either linear or cubic model did.It captured the smooth and consistent temporal shape even with very limited ground measurements during the rapid growth stage.The RMSE and R 2 of the final linear fit between the temporally corrected in-situ data and ASTER NDVI were 0.077 and 0.934 respectively; (iii) By comparing the MODIS and upscaled ASTER FAPAR temporally and spatially, the MODIS FAPAR manifested a good agreement to the aggregated FAPAR with R 2 of 0.922 and the RMSE of 0.054 in the study area.The performances of MODIS FAPAR product evaluated in the study is comparable to previous validation attempts, which have evaluated the MODIS C5 FAPAR product as maintaining an accuracy of about 0.1 [5,15,16,24,29,60,61].A similar validation over croplands by [24,29] found the accuracy of MODIS FAPAR product was about 0.07~0.08,which are very close to this study (0.054).

Figure 1 .
Figure 1.Map of study area and field measured sites.

Figure 1 .
Figure 1.Map of study area and field measured sites.

Figure 2 .
Figure 2. Methodology of scaling Fraction of Absorbed Photosynthetically Active Radiation (FAPAR) from the field to the satellite.

Figure 2 .
Figure 2. Methodology of scaling Fraction of Absorbed Photosynthetically Active Radiation (FAPAR) from the field to the satellite.

Figure 3 .
Figure 3. Modeling skylight proportions from solar zenith angle (SZA) and visibility (VIS) modeled from the second simulation of a satellite signal in the solar spectrum (6S) (blue signals represent 6S model outputs, whereas red lines represent the fitted equations).

Figure 3 .
Figure 3. Modeling skylight proportions from solar zenith angle (SZA) and visibility (VIS) modeled from the second simulation of a satellite signal in the solar spectrum (6S) (blue signals represent 6S model outputs, whereas red lines represent the fitted equations).

Figure 4 .
Figure 4. Effects of SZA and leaf area index (LAI) on FAPAR simulated from FAPAR-P model under VIS = 30 km: (a) effects of SZA on FAPAR with LAI = 1, 2, 4 and 8; and (b) effects of LAI on FAPAR with given SZA.Further, the combined effects of LAI, SZA and associated skylight proportions on instantaneous FAPAR were investigated, as shown in Figure5.It can be seen that the VIS does affect the derived FAPAR values, but in a subtle way at low SZA angles or in dense vegetation canopies with high LAI values (e.g., LAI = 8).At large SZA angles (>50°), the effects of VIS are amplified in sparse vegetation.

Figure 4 .
Figure 4. Effects of SZA and leaf area index (LAI) on FAPAR simulated from FAPAR-P model under VIS = 30 km: (a) effects of SZA on FAPAR with LAI = 1, 2, 4 and 8; and (b) effects of LAI on FAPAR with given SZA.

Figure 4 .
Figure 4. Effects of SZA and leaf area index (LAI) on FAPAR simulated from FAPAR-P model under VIS = 30 km: (a) effects of SZA on FAPAR with LAI = 1, 2, 4 and 8; and (b) effects of LAI on FAPAR with given SZA.Further, the combined effects of LAI, SZA and associated skylight proportions on instantaneous FAPAR were investigated, as shown in Figure5.It can be seen that the VIS does affect the derived FAPAR values, but in a subtle way at low SZA angles or in dense vegetation canopies with high LAI values (e.g., LAI = 8).At large SZA angles (>50°), the effects of VIS are amplified in sparse vegetation.

Figure 8 .
Figure 8. Validation of the temporal normalization method using in-situ FAPAR data.

Figure 8 .
Figure 8. Validation of the temporal normalization method using in-situ FAPAR data.

Figure 8 .
Figure 8. Validation of the temporal normalization method using in-situ FAPAR data.

Figure 9 .
Figure 9.Comparison of the fitting results from logistic, linear and cubic interpolations in: (a) correlation coefficient (R 2 ); and (b) RMSE.

Figure 9 .
Figure 9.Comparison of the fitting results from logistic, linear and cubic interpolations in: (a) correlation coefficient (R 2 ); and (b) RMSE.

Figure 11
Figure11shows the ASTER FAPAR derived from the empirical transfer function at 15 m resolution.Spatial aggregation was applied to the ASTER FAPAR at 15 m to derive the FAPAR truth data at 960 m.The mean value within each 64 × 64 window was calculated as the value for the new 960 m resolution pixel.This process is feasible as FAPAR is a self-consistent physical quantity across scales.Then FAPAR truth data were available at both 15 m and 960 m resolution, based on which direct comparison is now applicable.

Figure 11 .
Figure 11.(a-d) ASTER FAPAR derived from ASTER NDVI and in-situ FAPAR at four dates.

4. 3 .
Validating MODIS Product 4.3.1.Spatial Comparison The histograms of MODIS FAPAR were compared to the upscaled ASTER data, as shown in Figure 12.The distribution of MODIS FAPAR values in general coincided well with the upscaled ASTER FAPAR values.Especially on 15 June, the distribution of data between the pair of reference

Figure 11 .
Figure 11.(a-d) ASTER FAPAR derived from ASTER NDVI and in-situ FAPAR at four dates.

Figure 12 .
Figure 12. (a-d) Histograms of FAPAR values of upscaled ASTER and MODIS product at four dates.

Figure 12 .
Figure 12. (a-d) Histograms of FAPAR values of upscaled ASTER and MODIS product at four dates.

Figure 14 .
Figure 14.(a-m)Time series of MODIS FAPAR compared to the upscaled ASTER FAPAR at 13 in-situ sites (the bars represent one standard deviation (STD) of FAPAR values of ASTER pixels within the 1 km MODIS pixel cell).

Figure 14 .
Figure 14.(a-m)Time series of MODIS FAPAR compared to the upscaled ASTER FAPAR at 13 in-situ sites (the bars represent one standard deviation (STD) of FAPAR values of ASTER pixels within the 1 km MODIS pixel cell).

Figure 14 .
Figure 14.(a-m)Time series of MODIS FAPAR compared to the upscaled ASTER FAPAR at 13 in-situ sites (the bars represent one standard deviation (STD) of FAPAR values of ASTER pixels within the 1 km MODIS pixel cell).

Table 1 .
Data used in the presented study.

Table 1 .
Data used in the presented study.

Table 2 .
Input parameters and variables of FAPAR-P model.

Table 3 .
Root mean squared error (RMSE) of estimating FAPAR from the fitted k 1 and k 2 given VIS = 30 km.

Table 5 .
Statistics of difference between MODIS and upscaled ASTER FAPAR.

Table 5 .
Statistics of difference between MODIS and upscaled ASTER FAPAR.