Spectral Invariant Provides a Practical Modeling Approach for Future Biophysical Variable Estimations

This paper presents a simple radiative transfer model based on spectral invariant properties (SIP). The canopy structure parameters, including the leaf angle distribution and multi-angular clumping index, are explicitly described in the SIP model. The SIP model has been evaluated on its bidirectional reflectance factor (BRF) in the angular space at the radiation transfer model intercomparison platform, and in the spectrum space by the PROSPECT+SAIL (PROSAIL) model. The simulations of BRF by SIP agreed well with the reference values in both the angular space and spectrum space, with a root-mean-square-error (RMSE) of 0.006. When compared with the widely-used Soil-Canopy Observation of Photochemistry and Energy fluxes (SCOPE) model on fPAR, the RMSE was 0.006 and the R2 was 0.99, which shows a high accuracy. This study also suggests the newly proposed vegetation index, the near-infrared (NIR) reflectance of vegetation (NIRv), was a good linear approximation of the canopy structure parameter, the directional area scattering factor (DASF), with an R2 of 0.99. NIRv was not influenced much by the soil background contribution, but was sensitive to the leaf inclination angle. The sensitivity of NIRv to canopy structure and the robustness of NIRv to the soil background suggest NIRv is a promising index in future biophysical variable estimations with the support of the SIP model, especially for the Deep Space Climate Observatory (DSCOVR) Earth Polychromatic Imaging Camera (EPIC) observations near the hot spot directions.


Introduction
Canopy reflectance models have been widely used as the forward model for biophysical variable estimations, e.g., leaf area index (LAI), fractional vegetation cover (FVC) and the fraction of absorbed photosynthetically active radiation (fPAR) [1][2][3][4][5].However, there are a large variety of models to be selected, depending on the complexity of the canopy structure [6][7][8].For horizontally homogeneous canopies, e.g., cropland and grassland, there are 1D radiative transfer (RT) models, such as the PROSPECT+SAIL (PROSAIL) model and the two-layer Markov chain canopy reflectance (ACRM) model [9,10]; for discontinuous canopies, e.g., sparse forest, there are geometric-optical (GO)-based or hybrid models, such as the five-scale model, FRT model and SRTE model [11][12][13].There is a lack of a simple unified RT model that is suitable for both continuous and discontinuous canopies with different heterogeneities [14,15].The newly launched Deep Space Climate Observatory (DSCOVR) with the sensor Earth Polychromatic Imaging Camera (EPIC) provides unique near hot spot observations with a sub-hourly temporal frequency, which requires accurate RT forwarding modeling near the hot spot region for biophysical variable estimations [16].
Spectral invariant theory simplifies the radiative transfer process by decoupling the canopy structure and leaf optics [17].There are three key spectral invariant parameters to describe the canopy structure: the canopy interceptance i 0 , recollision probability p and escape probability ρ [18][19][20].The canopy interceptance i 0 is the fraction of photons that are intercepted by leaves in the canopy; the recollision probability p is defined as the probability that a photon scattered from a leaf in a canopy will interact within the canopy again; the escape probability ρ is defined as the probability that a scattered photon will escape the canopy.By the three wavelength-independent parameters, it is no longer necessary to use the successive orders of scattering approximation (SOSA) method for the multiple scattering calculation, which is usually iteration-based and time-consuming [21].The spectral invariant has recently been used in the remotely sensed albedo [22], fPAR [23] and sunlit LAI estimations [22][23][24].However, there are still several aspects of the spectral invariant that need to be improved for future biophysical variable estimations.
The first aspect is how to explicitly describe the spectral invariant by canopy structure parameters, including the leaf angle distribution (LAD) and clumping index (CI) [25,26].Only by linking the canopy structure parameters with the abstract spectral invariant can we estimate the canopy structure parameters from remotely sensed observations.Stenberg [27] found a simple analytical formula for the relationship between LAI and recollision probability, while LAD and CI are still important to explicitly account for.In fact, LAD is one of the key parameters that determine the angular distribution of the reflected photons.It is common to assume the LAD as spherical distribution or a given type for each biome [17], while canopy-level studies have shown that such an assumption can result in as much as a 47% underestimation in LAI if the true LAD distribution is better characterized as planophile [28].Besides this, accounting for the clumping effects is the key to making the 1D RT model applicable to discontinuous canopies.CI measurements are primarily recorded as hemispherically integrated CI [29], despite the fact that CI strongly varies with the view zenith angle [30][31][32].The angular variation effect of CI is still important to describe in the RT process.
The second challenge is how to explicitly account for the hot spot effect in the RT process by the spectral invariant, which is due to the finite size of the leaf element [33].The hot spot effect is that the reflectance reaches the local maximum when the view and solar direction coincide, because there is no shadow in view.Thus, the single scattering contribution by the first-order escape probability needs to be treated separately.The simple semi-physical parameterization (PARAS) model, which is forest spectral invariant-based, does not include the hot spot effect yet [34].Recently Yang, et al. [24] began to incorporate the hot spot phenomenon by separating the escape probability of sunlit leaves and shaded leaves.Up to now, as far as we know, the spectral invariant-based bidirectional reflectance factor (BRF) models still have not been validated at the radiation transfer model intercomparison (RAMI) online model checker (ROMC) platform [35].
The last aspect is how to describe the multiple scattering between vegetation and soil and reduce the soil contamination in the biophysical variable estimations.The multiple scattering between vegetation and soil can still make a considerable contribution to the BRF when the LAI is low and the soil is bright, e.g., sparse forest with snow background, while this effect still has not been accounted for in the PARAS model [34].Besides this, the direct extraction of recollision probability and directional area scattering factor (DASF) by the BRF-fitting approach assumes soil background effects are negligible (black soil problem) [18,36].If not, a full RT model by stochastic radiative transfer equation (SRTE), would be used to remove the soil contribution (soil problem), which reduces the simplicity of the spectral invariant approach [13,24].In fact, DASF is a canopy structure parameter and is equal to the canopy BRF when there is no soil background contribution and the leaves have no absorption (leaf albedo w = 1).Thus, it is important to find a simple index that can be directly extracted from BRF observations, and at the same time is robust to soil brightness in the biophysical variable estimations, e.g., the newly developed vegetation index, the near-infrared (NIR) reflectance of vegetation (NIRv) [37].NIRv is the product of the total NIR reflectance and the normalized difference vegetation index (NDVI), and has been shown to be insensitive to soil background contribution [37].By NIRv, we do not need to account for the single scattering contribution of the soil and the photon interactions between vegetation and soil, and this is therefore promising for future biophysical estimations.
To address the above three challenges, we develop a simple RT model based on spectral invariant properties (SIP) which explicitly describes the leaf angle distribution (LAD) and clumping index (CI) and treats the single scattering contribution separately.Then, we evaluate the SIP model of BRF by the RAMI platform in the angular space [35], compare the SIP model with the PROSAIL model in the spectrum space [10], and evaluate the SIP model of fPAR by the Soil-Canopy Observation of Photochemistry and Energy fluxes (SCOPE) model [38], a widely accepted model that can directly simulate fPAR as the output.The input parameters are varied at different canopy structures, solar angles, and soil backgrounds.Then, the impact of the soil effect on DASF calculation has been evaluated by the proposed SIP model.Finally, the impact of LAD on BRF and NIRv has been evaluated.

Model Description
Accounting for the fact that the first-order escape probability has a stronger anisotropic effect than high-order ones [24], the SIP model proposed in this study adopts the method of calculating the single scattering contribution from the PROSAIL model [10], but meanwhile explicitly describes the clumping effects and the leaf angle distribution.The multiple scattering within the canopy, and the interactions between the canopy and the soil, are calculated in the framework of the spectral invariant, which has the advantage of simplicity compared with the SOSA method [21].The framework of the SIP model is in Figure 1, and detailed descriptions of each item are in Section 2 as follows.
Remote Sens. 2018, 10, x FOR PEER REVIEW 3 of 17 are negligible (black soil problem) [18,36].If not, a full RT model by stochastic radiative transfer equation (SRTE), would be used to remove the soil contribution (soil problem), which reduces the simplicity of the spectral invariant approach [13,24].In fact, DASF is a canopy structure parameter and is equal to the canopy BRF when there is no soil background contribution and the leaves have no absorption (leaf albedo w = 1).Thus, it is important to find a simple index that can be directly extracted from BRF observations, and at the same time is robust to soil brightness in the biophysical variable estimations, e.g., the newly developed vegetation index, the near-infrared (NIR) reflectance of vegetation (NIRv) [37].NIRv is the product of the total NIR reflectance and the normalized difference vegetation index (NDVI), and has been shown to be insensitive to soil background contribution [37].By NIRv, we do not need to account for the single scattering contribution of the soil and the photon interactions between vegetation and soil, and this is therefore promising for future biophysical estimations.To address the above three challenges, we develop a simple RT model based on spectral invariant properties (SIP) which explicitly describes the leaf angle distribution (LAD) and clumping index (CI) and treats the single scattering contribution separately.Then, we evaluate the SIP model of BRF by the RAMI platform in the angular space [35], compare the SIP model with the PROSAIL model in the spectrum space [10], and evaluate the SIP model of fPAR by the Soil-Canopy Observation of Photochemistry and Energy fluxes (SCOPE) model [38], a widely accepted model that can directly simulate fPAR as the output.The input parameters are varied at different canopy structures, solar angles, and soil backgrounds.Then, the impact of the soil effect on DASF calculation has been evaluated by the proposed SIP model.Finally, the impact of LAD on BRF and NIRv has been evaluated.

Model Description
Accounting for the fact that the first-order escape probability has a stronger anisotropic effect than high-order ones [24], the SIP model proposed in this study adopts the method of calculating the single scattering contribution from the PROSAIL model [10], but meanwhile explicitly describes the clumping effects and the leaf angle distribution.The multiple scattering within the canopy, and the interactions between the canopy and the soil, are calculated in the framework of the spectral invariant, which has the advantage of simplicity compared with the SOSA method [21].The framework of the SIP model is in Figure 1, and detailed descriptions of each item are in Section 2 as follows.

BRF of the Vegetation-Soil System
The BRF of the vegetation-soil system at the wavelength of λ and at the solar-sensor geometry Ω contains three components: where BRF V is the contribution of the photons from vegetation with no interaction with the soil background, BRF S is the single scattering contribution of the photons from soil background with no interaction with the vegetation, and BRF M is the multiple scattering contribution of the photons between the vegetation and soil.The solar-sensor geometry Ω contains the solar-target geometry and target-sensor geometry, and both of them can be described by the zenith angle and azimuth angle.
(1) Photons which only interact with the vegetation (BRF V ) BRF V comprises the single scattering of the vegetation (BRF V1 ) and multiple scattering within the vegetation with no interaction with the soil (BRF V M ): where BRF V1 can be derived as w so is the average bidirectional scattering coefficient in the PROSAIL model integrated over all the leaf inclination angles from 0 • to 90 • [10], which is a function of leaf reflectance/transmittance at the wavelength of λ, leaf inclination angle θ l , and the solar-sensor geometry Ω. u L is the leaf area density (m 2 /m 3 ), H is the canopy height, and P so (Ω, z) is the bidirectional gap fraction at height z, which is a function of the directional gap fraction at height z with the solar zenith angle θ s (P s ) and the view zenith angle θ o (P o .), and the hot spot correction function at height z with a phase angle ξ of the solar-sensor direction [33].Here τ and τ are the extinction coefficients in the solar and view directions, respectively; s and s are the optical lengths intersected the vegetated component in the solar and view directions, respectively; and ξ is the scattering phase angle between the solar and view directions.k l = 1/s L , and s L is the characteristic linear dimension of the leaves [33].S AB is the distance between the two points A and B where the solar and the observation rays penetrate into the canopy [15], and S AB can be computed as S AB = s 2 + s 2 − 2ss cosξ, where ξ can be calculated from the solar/view zenith angle θ s and θ o , respectively, and the solar/view azimuth angle ϕ s and ϕ o , respectively, as cosξ = cosθ s cosθ o + sinθ s sinθ o cos(ϕ s − ϕ o ) [15,33].
The directional gap fraction P s (θ s , z) at height z with the solar zenith angle θ s can be described as where G(θ s ) is the leaf projection function at the solar direction and is a function of the leaf inclination angle θ l , µ s = cosθ s , LAI(z) is the culmulative LAI at height z from the top, LAI(H) is the canopy LAI, and CI(θ s ) is the clumping index at the solar zenith angle θ s .When we replace the solar zenith angle θ s in Equation ( 5) by the view zenith angle θ o , we can get the directional gap fraction P o (θ o , z) at height z with the view zenith angle θ o .By Equations ( 3) and ( 5), the leaf inclination angle and multi-angular clumping index have been explicitly accounted for in the SIP model.
The multiple scattering within the vegetation BRF V M can be derived by the spectral invariant as where ω(λ) is the leaf albedo, and the method for linking the canopy structure parameters with the canopy directional interceptance i 0 , the directional escape probability ρ(Ω) and recollision probability p will be described in Section 2.4.
(2) Photons which only interact with the soil (BRF S ) In fact, BRF S is the contribution of the sunlit soil, and can be described as where ρ G (λ) is the soil reflectance, and K g is the fraction of visible sunlit soil in view, which is equal to the bidirectional gap fraction at the bottom of the canopy P so (Ω, H) in Equation ( 4).
(3) Photons which interact with the both the vegetation and soil (BRF M ) BRF M is the multiple scattering between vegetation and soil, and can be described as where T dn and T up are the canopy downwelling and upward transmittance which will be described in Section 2.3, R dn (λ) is the albedo in the downward direction under isotropic diffuse illumination from the bottom of the canopy, t 0 = 1 − i 0 , and t di f 0 is the canopy directional transmittance under diffuse radiation illumination.

Absorption of the Vegetation and fPAR
The total photons absorbed by the vegetation comprise two components: A V (λ) represents the absorbed photons that only interact with the vegetation [17]: A M (λ) represents the photons that interact with both the vegetation and the soil, but are finally absorbed by the vegetation: where R dn (λ) is the albedo in the downward direction under isotropic diffuse illumination from the bottom of the canopy, and A up (λ) is the corresponding absorption by vegetation.The two terms will be described in Section 2.3.Finally, fPAR is the integration of the vegetation absorption A(λ) and the incoming solar irradiance R in (λ) across the range 400 nm-700 nm: Remote Sens. 2018, 10, 1508 6 of 17

Canopy Radiation Transfer Terms by Spectral Invariant
The canopy level radiation transfer terms in Equations ( 8) and (11) will be described by the spectral invariant in this section.The canopy downwelling transmittance with black soil background in Equation ( 8) can be described as (13) where τ hemi is the hemispherical escape probability in the downward direction, that will be described in Section 2.4.The canopy upward transmittance under isotropic illumination from the bottom of the canopy, but with black soil background in Equation ( 8), can be described as 1−pω(λ) (14) where t di f 0 is the canopy directional transmittance under diffuse radiation illumination, and i D is the hemispherical canopy interceptance under diffuse radiation illumination, which will be described in Section 2.4.
The albedo in the downward direction under isotropic diffuse illumination from the bottom of the canopy, but with black soil background in Equation (11), can be described as (15) where ρ hemi is the hemispherical escape probability in the upward direction that will be described in Section 2.4.Similar to Equation (10), the absorption by vegetation under isotropic diffuse illumination from the bottom of the canopy, but with black soil background in Equation (11), can be described as 2.4.Analytical Formula of Spectral Invariant by Canopy Structure (1) Canopy interceptance The directional canopy interceptance at solar direction , where θ s is the solar zenith angle, µ s = cosθ s , and CI(θ s ) is the clumping index at the solar direction.
The hemispherical canopy interceptance under diffuse radiation illumination which can be achieved by Gaussian integral.

2) Escape probability
The directional escape probability for diffuse radiation, or high order (>1) directional escape probability for direct radiation, can be described by directional interceptance i(θ): Similarly, the hemispherical escape probability for diffuse radiation, or high order (>1) hemispherical escape probability for direct radiation, can be described by hemispherical interceptance i D Here, the high order (>1) upward and downward escape probability have the derivation as ρ(Ω) = τ(Ω), and ρ hemi = τ hemi .
(3) Recollision probability The recollision probability for both direct and diffuse radiation can be derived by hemispherical interceptance i D [27]

Materials and Methods
First, we evaluate the SIP model by the RAMI platform (http://romc.jrc.ec.europa.eu/) in the angular space [35].RAMI has been developed since the 1990s and is the most comprehensive benchmark for the validation of radiation transfer models so far [39].The RAMI online model checker provides a web-based platform for model developers to autonomously evaluate the performance of the canopy radiation transfer models, and the "surrogate truth" for the evaluation is generated by the means of credible 3-D models, including Discrete Anisotropic Radiative Transfer (DART), librat, raytran and rayspread [35].The "validate" mode of the canopy scenario in RAMI platform is in Table 1, and the leaves with finite size are randomly distributed in space.As the RAMI was designed for the validation of the multi-angular and multispectral simulations, we compare the SIP model with the PROSAIL model in the spectrum space, because the PROSAIL model has been widely used to simulate the hyperspectral observations at the canopy scale [10].The fPAR estimation by SIP at different canopy structures, solar angles, and soil background has been evaluated by the SCOPE model [38].For canopy structure, we varied the FVC, LAI and LAD.Different soil brightnesses have also been accounted for in the fPAR evaluation.The list of variables and their ranges are listed in Table 2.Then, the impact of the soil effect on DASF calculation has been evaluated.DASF was calculated by the standard method from the fitting of the hyperspectral BRF simulations (BRF/ω vs. BRF, where ω is the leaf albedo) within the wavelength range of 710-790 nm [18].Finally, the impact of LAD on BRF and NIRv has also been evaluated [37].

Evaluation in the RAMI Platform in Angular Space
The evaluation of the SIP model the by the RAMI platform is shown in Figure 2. The simulations of the SIP model at red and NIR bands agreed well with the reference "surrogate truth" by the RAMI On-line Model Checker (ROMC) reference dataset (ROMCREF).The simulation error of SIP at red and NIR bands can be within the ±5.0%bias region, and most of the simulation errors were within the ±2.5% bias region.The scatter plot is shown in Figure 3, which suggests the simulations of the SIP model closely followed the 1:1 line.The root-mean-square (RMS) error of the SIP model was 0.0057, which was less than 0.01.The signal-to-noise ratio (S/R) of the SIP model was 49.0, which was significantly larger than 1.There was a slight underestimation of BRF in the red band and a slight overestimation in the NIR band away from the hot spot direction, while this was more accurate at the hot spot direction, which might be due to the effect of the hot spot function.This can be affected by the leaf size, leaf shape and numerical singularities in the calculation.The hot spot effect has been considered in the first order scattering in both the red and NIR bands.In total, the evaluation results by the RAMI platform show that the high accuracy of the SIP model can meet the accuracy requirement of the BRF simulations in the angular space, especially at the hot spot direction.

Evaluation in the RAMI Platform in Angular Space
The evaluation of the SIP model the by the RAMI platform is shown in Figure 2. The simulations of the SIP model at red and NIR bands agreed well with the reference "surrogate truth" by the RAMI On-line Model Checker (ROMC) reference dataset (ROMCREF).The simulation error of SIP at red and NIR bands can be within the ±5.0%bias region, and most of the simulation errors were within the ±2.5% bias region.The scatter plot is shown in Figure 3, which suggests the simulations of the SIP model closely followed the 1:1 line.The root-mean-square (RMS) error of the SIP model was 0.0057, which was less than 0.01.The signal-to-noise ratio (S/R) of the SIP model was 49.0, which was significantly larger than 1.There was a slight underestimation of BRF in the red band and a slight overestimation in the NIR band away from the hot spot direction, while this was more accurate at the hot spot direction, which might be due to the effect of the hot spot function.This can be affected by the leaf size, leaf shape and numerical singularities in the calculation.The hot spot effect has been considered in the first order scattering in both the red and NIR bands.In total, the evaluation results by the RAMI platform show that the high accuracy of the SIP model can meet the accuracy requirement of the BRF simulations in the angular space, especially at the hot spot direction.

Evaluation by the PROSAIL Model in Hyperspectral Space
The evaluation of the SIP model by the PROSAIL model in hyperspectral space is shown in Figure 4.The simulations by SIP agreed well with that by PROSAIL across all the wavelengths from 400 nm-2500 nm.The valley in the red band at the canopy scale (Figure 4b) was not as significant as the leaf scale (Figure 4a), which was due to the soil background contribution at low LAI (LAI = 1).The canopy scale simulations by SIP were more accurate at low leaf reflectance (red) than at high leaf reflectance (NIR), which means there was smaller uncertainties in the single scattering calculations by SIP, compared with the multiple scattering calculations.In total, the evaluation result by the PROSAIL model shows the high accuracy of the SIP model when the LAI was 1, 3 and 5, respectively.

Evaluation by the PROSAIL Model in Hyperspectral Space
The evaluation of the SIP model by the PROSAIL model in hyperspectral space is shown in Figure 4.The simulations by SIP agreed well with that by PROSAIL across all the wavelengths from 400 nm-2500 nm.The valley in the red band at the canopy scale (Figure 4b) was not as significant as the leaf scale (Figure 4a), which was due to the soil background contribution at low LAI (LAI = 1).The canopy scale simulations by SIP were more accurate at low leaf reflectance (red) than at high leaf reflectance (NIR), which means there was smaller uncertainties in the single scattering calculations by SIP, compared with the multiple scattering calculations.In total, the evaluation result by the PROSAIL model shows the high accuracy of the SIP model when the LAI was 1, 3 and 5, respectively.400 nm-2500 nm.The valley in the red band at the canopy scale (Figure 4b) was not as significant as the leaf scale (Figure 4a), which was due to the soil background contribution at low LAI (LAI = 1).The canopy scale simulations by SIP were more accurate at low leaf reflectance (red) than at high leaf reflectance (NIR), which means there was smaller uncertainties in the single scattering calculations by SIP, compared with the multiple scattering calculations.In total, the evaluation result by the PROSAIL model shows the high accuracy of the SIP model when the LAI was 1, 3 and 5, respectively.1.

Evaluation the fPAR by the SCOPE Model
The fPAR estimation by SIP at different canopy structures (e.g., FVC, LAI and LAD), solar angles, and soil backgrounds was evaluated by the widely-used SCOPE model (Figure 5).The fPAR estimated by SIP agreed well with the estimations by SCOPE, and they were almost on the 1:1 line.The R 2 was 0.99, and the RMSE was 0.006, which show the high accuracy of the fPAR estimation by SIP.It is also obvious that there was no systematic error on the fPAR estimation by SIP.This was in accordance with the high accuracy of the BRF simulation in red band in Figure 2.  2.  1.

Evaluation the fPAR by the SCOPE Model
The fPAR estimation by SIP at different canopy structures (e.g., FVC, LAI and LAD), solar angles, and soil backgrounds was evaluated by the widely-used SCOPE model (Figure 5).The fPAR estimated by SIP agreed well with the estimations by SCOPE, and they were almost on the 1:1 line.The R 2 was 0.99, and the RMSE was 0.006, which show the high accuracy of the fPAR estimation by SIP.It is also obvious that there was no systematic error on the fPAR estimation by SIP.This was in accordance with the high accuracy of the BRF simulation in red band in Figure 2.

Impact of Soil Effect on DASF Calculation
Figure 6 shows the polar map of NIRv and DASF calculated directly from the total BRF of the soil-canopy system (DASF t ), and DASF calculated after the removal of the soil contribution as the black soil condition (DASF bs ).The corresponding scatter plots are shown in Figure 7. Figure 6 suggests the hot spot directions were more easily influenced by the soil due to the larger sunlit soil contribution.In the angular space, NIRv had a similar shape to DASF bs , while DASF t shows a different trend, especially the obvious high values around the hot spot directions.In Figure 7, NIRv was more linearly related with DASF bs , with an R 2 of 0.99 and RMSE of 0.006, while DASF t had a worse relationship with DASF bs , with an R 2 of 0.82 and RMSE of 0.023.It is obvious to see the outlier for DASF t at the hot spot direction in Figure 7.This suggests the simple index NIRv was promising as a good linear approximation of the important canopy structure parameter DASF bs , without the use of a full radiative transfer model to eliminate the soil contribution.

Evaluation the fPAR by the SCOPE Model
The fPAR estimation by SIP at different canopy structures (e.g., FVC, LAI and LAD), solar angles, and soil backgrounds was evaluated by the widely-used SCOPE model (Figure 5).The fPAR estimated by SIP agreed well with the estimations by SCOPE, and they were almost on the 1:1 line.The R 2 was 0.99, and the RMSE was 0.006, which show the high accuracy of the fPAR estimation by SIP.It is also obvious that there was no systematic error on the fPAR estimation by SIP.This was in accordance with the high accuracy of the BRF simulation in red band in Figure 2.  2.

Impact of Soil Effect on DASF Calculation
Figure 6 shows the polar map of NIRv and DASF calculated directly from the total BRF of the soil-canopy system (DASFt), and DASF calculated after the removal of the soil contribution as the black soil condition (DASFbs).The corresponding scatter plots are shown in Figure 7. Figure 6 suggests the hot spot directions were more easily influenced by the soil due to the larger sunlit soil contribution.In the angular space, NIRv had a similar shape to DASFbs, while DASFt shows a different trend, especially the obvious high values around the hot spot directions.In Figure 7, NIRv was more linearly related with DASFbs, with an R 2 of 0.99 and RMSE of 0.006, while DASFt had a worse relationship with DASFbs, with an R 2 of 0.82 and RMSE of 0.023.It is obvious to see the outlier for The LAI is 1, LAD is erectophile, and the solar zenith angle is 30°.The other input parameters are the same as Table 1.
Figure 7.The scatter plot between NIRv or DASF calculated from total scene BRF without the removal of the soil (DASFt), with DASF calculated at the black soil condition (DASFbs).The LAI is 1, and the solar zenith angle is 30°.The other input parameters are the same as Table 1 or Figure 5.

Impact of the Mean Leaf Inclination Angle on BRF and NIRv
The impact of LAD on BRF has been evaluated by SIP at different mean leaf inclination angles (MLA) by the ellipsoidal distribution function adopted in the PROSAIL model [10].Figure 8 suggests that, even when the other parameters were the same as, e.g., LAI, only the difference in LAD can lead to significant BRF variations in both the red band and NIR band.This also suggests it was important The LAI is 1, LAD is erectophile, and the solar zenith angle is 30°.The other input parameters are the same as Table 1.
Figure 7.The scatter plot between NIRv or DASF calculated from total scene BRF without the removal of the soil (DASFt), with DASF calculated at the black soil condition (DASFbs).The LAI is 1, and the solar zenith angle is 30°.The other input parameters are the same as Table 1 or Figure 5.

Impact of the Mean Leaf Inclination Angle on BRF and NIRv
The impact of LAD on BRF has been evaluated by SIP at different mean leaf inclination angles (MLA) by the ellipsoidal distribution function adopted in the PROSAIL model [10].Figure 8 suggests that, even when the other parameters were the same as, e.g., LAI, only the difference in LAD can lead to significant BRF variations in both the red band and NIR band.This also suggests it was important Figure 7.The scatter plot between NIRv or DASF calculated from total scene BRF without the removal of the soil (DASF t ), with DASF calculated at the black soil condition (DASF bs ).The LAI is 1, and the solar zenith angle is 30 • .The other input parameters are the same as Table 1 or Figure 5.

Impact of the Mean Leaf Inclination Angle on BRF and NIRv
The impact of LAD on BRF has been evaluated by SIP at different mean leaf inclination angles (MLA) by the ellipsoidal distribution function adopted in the PROSAIL model [10].Figure 8 suggests that, even when the other parameters were the same as, e.g., LAI, only the difference in LAD can lead to significant BRF variations in both the red band and NIR band.This also suggests it was important to have the LAD information in biophysical parameter estimations, or a large inversion error may be introduced.1.

Discussion
The SIP model proposed in this study separates the first-order escape probability by the bidirectional gap fraction, which keeps the simplicity of the spectral invariant, and meanwhile explicitly describes the clumping effects and the leaf angle distribution.The SIP model shows high accuracy of BRF in the angular space and spectral space.The simplicity and high accuracy show that it is promising for application in future reflection-based biophysical variable estimations.
The SIP model is especially accurate at the hotspot direction, which has potential for biophysical variable estimation with the newly launched DSCOVR EPIC [16].At the hot spot direction, the bidirectional gap fraction in Equation ( 4) would be reduced to the directional gap fraction or because the • or • is equal to 1.This greatly simplifies the radiative transfer process, and would largely reduce the uncertainties in the hot spot correction function in Equation ( 4).The angular distribution of NIRv especially in the hot spot directions shows the sensitivity of NIRv to the mean leaf inclination angle (Figure 9), which suggests that it is promising for estimating the mean leaf inclination angle by NIRv from multi-angular EPIC observations and the proposed SIP model.Usually, the effects of the leaf inclination angle and multi-angular clumping index are coupled,  1.
Figure 9 shows the NIRv angular variation in the principal plane when the solar zenith angle was 50 • and in the hot spot direction.This suggests a different LAD can lead to different NIRv angular variation, especially for the hot spot NIRv.The more erect the leaves are, the larger the angular variation of hot spot NIRv.It is interesting to find that all the NIRv curves converged when the view zenith angle was about 57 • (Figure 9b), and this is because at 57 • the leaf projection G-function gives almost the same values regardless of the leaf inclination angle.This also suggests theoretically that, by using only two angular NIRv measurements (one was at 57 • view zenith angle), it would be possible to extract the mean leaf inclination angle.Of course, more angular samplings of NIRv at the hot spot direction would be possible to achieve more robust mean leaf inclination angle estimations.1.
(a) (b) Figure 9.The NIRv simulated by the SIP model when the solar zenith angle is 50° (a), and when the solar direction coincides with the view direction for hot spot observation (b).MLA represents the mean leaf inclination angle.The LAI is 1, and other input parameters are the same as Table 1.

Discussion
The SIP model proposed in this study separates the first-order escape probability by the bidirectional gap fraction, which keeps the simplicity of the spectral invariant, and meanwhile explicitly describes the clumping effects and the leaf angle distribution.The SIP model shows high  (a), and when the solar direction coincides with the view direction for hot spot observation (b).MLA represents the mean leaf inclination angle.The LAI is 1, and other input parameters are the same as Table 1.

Discussion
The SIP model proposed in this study separates the first-order escape probability by the bidirectional gap fraction, which keeps the simplicity of the spectral invariant, and meanwhile explicitly describes the clumping effects and the leaf angle distribution.The SIP model shows high accuracy of BRF in the angular space and spectral space.The simplicity and high accuracy show that it is promising for application in future reflection-based biophysical variable estimations.
The SIP model is especially accurate at the hotspot direction, which has potential for biophysical variable estimation with the newly launched DSCOVR EPIC [16].At the hot spot direction, the bidirectional gap fraction P so in Equation ( 4) would be reduced to the directional gap fraction P s or P o because the P o •C HS or P s •C HS is equal to 1.This greatly simplifies the radiative transfer process, and would largely reduce the uncertainties in the hot spot correction function C HS in Equation ( 4).The angular distribution of NIRv especially in the hot spot directions shows the sensitivity of NIRv to the mean leaf inclination angle (Figure 9), which suggests that it is promising for estimating the mean leaf inclination angle by NIRv from multi-angular EPIC observations and the proposed SIP model.Usually, the effects of the leaf inclination angle and multi-angular clumping index are coupled, while their retrievals can still be constrained by their angular features, e.g., the NIRv would converge when the view zenith angle was about 57 • for any leaf inclination angle (Figure 9), and the results of Kucharik et al. [32] show that the clumping index with a zenith angle between 20 • -60 • obeys an exponential or even linear function.Due to the sub-hourly high temporal frequency, the EPIC can almost always acquire two observations with a solar-sensor zenith angle of about 57 • each day, which shows the feasibility of retrievals.
Besides this, the estimation of biophysical variables requires the elimination of the soil background contribution.The simple vegetation index, NIRv, was shown to give a good linear approximation of the canopy structure parameter DASF (Figure 7) and was not influenced much by the soil background effect (Figure 6).Thus, it is worth trying NIRv for the cost function of biophysical variable retrievals, in contrast to the independent use of the red and NIR band in the look-up-table [17].By using NIRv, we do not need to account for the single scattering contribution of the soil (BRF S in Equation ( 7)) and the photon interactions between vegetation and soil (BRF M in Equation ( 8)), and thus we only need to calculate the BRF V in Equation (1), which would further simplify the radiative transfer process.The high accuracy of the proposed SIP model across the whole wavelength (400 nm-2500 nm) in Figure 4 also shows it would be promising to retrieve the canopy structure and leaf optics at the same time by the joint use of hyperspectral and multi-angular observations.To improve the retrieval accuracy, the recent progress in the field of RT forwarding modeling, inversion and validation should be accounted for, such as the mixed pixel effect [40,41], the variation of the vegetation density [42][43][44], the topography effect [45][46][47][48], and the spatial-temporal constraints [49][50][51].

Conclusions
We develop a simple radiative transfer model based on spectral invariant properties (SIP).The SIP model adopts the method of calculate the single scattering contribution from the PROSAIL model, but explicitly describes the clumping effects and the leaf angle distribution.The multiple scattering within the canopy, and the interactions between the canopy and the soil, are calculated in the framework of the spectral invariant, which has the advantage of simplicity compared with the SOSA method [21].The SIP model agreed well with the reference values on BRF and fPAR from the RAMI platform, the PROSAIL model and the SCOPE model, with an R 2 of 0.99 and an RMSE of 0.006.We also found that NIRv gave a good linear approximation of the canopy structure parameter, DASF.NIRv was robust to the soil background contribution, but was sensitive to the leaf inclination angle.Results suggest NIRv is a promising index for future biophysical variable estimations by the proposed SIP model, especially for the joint use of the hyperspectral and multi-angular observations.

Figure 1 .
Figure 1.The framework of the spectral invariant properties (SIP) model on how to describe the radiative transfer process from spectral invariants, which are derived bfromcanopy structure parameters.Detailed descriptions of each item are in Section 2.

Figure 1 .
Figure 1.The framework of the spectral invariant properties (SIP) model on how to describe the radiative transfer process from spectral invariants, which are derived bfromcanopy structure parameters.Detailed descriptions of each item are in Section 2.

Figure 2 . 17 Figure 3 .
Figure 2. The bidirectional reflectance factor (BRF) of the SIP model in the red (a) and near infra-red (NIR) bands with the solar zenith angle at 0 • (b), 30 • (c), 60 • (d) in the orthogonal plane evaluated by the RAMI On-line Model Checker (ROMC) reference dataset (ROMCREF)."HOM11" is the number of the standard RAMI scene, "DIS" represents the discrete medium of the leaves, "ERE" represents the erectophile leaf angle distribution (LAD) type, and "OP" means the orthogonal plane.Remote Sens. 2018, 10, x FOR PEER REVIEW 9 of 17

Figure 3 .
Figure 3.The scatter plots between the BRF of the SIP simulations and the RAMI reference values (ROMCREF) in the red and NIR bands in the orthogonal plane at the RAMI platform in Figure 2. The root mean square (RMS) error and signal-to-noise (S/N) were directly calculated with the RAMI platform (http://romc.jrc.ec.europa.eu/).

Figure 4 .
Figure 4.The input leaf/soil spectrum (a) and canopy-scale spectrum simulations by SIP and PROSAIL (b-d).VZA represents the view zenith angle.The leaf area index (LAI) is 1 (b), 3 (c) or 5 (d), respectively, LAD is erectophile, and the solar zenith angle is 0°.The other input parameters are the same as Table1.

Figure 5 .
Figure 5.The evaluation of the SIP model by the widely-used SCOPE model at different canopy structures, solar angles, and soil backgrounds.The list of variables and their ranges are in Table2.

Figure 4 .
Figure 4.The input leaf/soil spectrum (a) and canopy-scale spectrum simulations by SIP and PROSAIL (b-d).VZA represents the view zenith angle.The leaf area index (LAI) is 1 (b), 3 (c) or 5 (d), respectively, LAD is erectophile, and the solar zenith angle is 0 • .The other input parameters are the same as Table1.

Figure 5 .
Figure 5.The evaluation of the SIP model by the widely-used SCOPE model at different canopy structures, solar angles, and soil backgrounds.The list of variables and their ranges are in Table2.

Figure 5 .Figure 6 .
Figure 5.The evaluation of the SIP model by the widely-used SCOPE model at different canopy structures, solar angles, and soil backgrounds.The list of variables and their ranges are in Table2.

Figure 6 .Figure 6 .
Figure 6.The polar map of NIRv, (a) the directional area scattering factor (DASF) calculated without the removal of the soil (b), and DASF calculated after the removal of the soil (c).The LAI is 1, LAD is erectophile, and the solar zenith angle is 30 • .The other input parameters are the same as Table1.

Figure 8 .Figure 9 .
Figure 8.The reflectance simulated by the SIP model in the red band (a) and NIR band (b) at different mean leaf inclination angles (MLA) in the principal plane.The LAI is 1, and the solar zenith angle is 30°.The other input parameters are the same as Table1.

Figure 8 .
Figure 8.The reflectance simulated by the SIP model in the red band (a) and NIR band (b) at different mean leaf inclination angles (MLA) in the principal plane.The LAI is 1, and the solar zenith angle is 30 • .The other input parameters are the same as Table1.

Figure 8 .
Figure 8.The reflectance simulated by the SIP model in the red band (a) and NIR band (b) at different mean leaf inclination angles (MLA) in the principal plane.The LAI is 1, and the solar zenith angle is 30°.The other input parameters are the same as Table1.

Figure 9 .
Figure 9.The NIRv simulated by the SIP model when the solar zenith angle is 50 • (a), and when the solar direction coincides with the view direction for hot spot observation (b).MLA represents the mean leaf inclination angle.The LAI is 1, and other input parameters are the same as Table1.

Table 2 .
Parameters varied in SIP and Soil-Canopy Observation of Photochemistry and Energy fluxes (SCOPE) simulations for the fraction of absorbed photosynthetically active radiation (fPAR) estimation.
BRFbidirectional reflectance factor BRF V BRF contribution by vegetation with no interaction with the soil background BRF V1single scattering of the vegetation in BRF BRF V M multiple scattering contribution of the vegetation in BRF BRF S single scattering contribution of the soil background in BRF BRF M (λ,Ω)multiple scattering contribution between the vegetation and soil in BRF canopy zero-order transmission at the direction of Ω under diffuse radiation t D canopy zero-order hemispherical transmission under diffuse radiation A up absorption by vegetation under isotropic diffuse illumination at bottom R dn albedo in the downward direction under isotropic diffuse illumination at bottom T dn canopy downwelling transmittance with black soil background T up canopy upward transmittance with isotropic diffuse illumination at bottom