Development of a New BRDF-Resistant Vegetation Index for Improving the Estimation of Leaf Area Index

Su Zhang 1,2, Liangyun Liu 1,*, Xinjie Liu 1 and Zefei Liu 1,3 1 Key Laboratory of Digital Earth Science, Institute of Remote Sensing and Digital Earth, Chinese Academy of Sciences, Beijing 100094, China; zhangsu@radi.ac.cn (S.Z.); liuxj@radi.ac.cn (X.L.); liuzefei1990@163.com (Z.L.) 2 University of Chinese Academy of Sciences, Beijing 100049, China 3 College of Geometrics, Xi’an University of Science and Technology, Xi’an 710054, China * Correspondence: liuly@radi.ac.cn; Tel.: +86-10-8217-8163


Introduction
The leaf area index (LAI) is one of the most important Earth surface parameters used in modeling ecosystems and their interaction with climate [1,2].It provides key information about crop growth and is highly correlated with crop biomass and productivity [3][4][5].
Vegetation with different LAI values has different spectral reflectance characteristics [6,7].Accordingly, numerous vegetation indices have been developed to estimate the LAI.Delegido et al. proposed a red-edge normalized difference index (NDI) that is correlated strongly with the LAI without saturating at larger values [8].Viña et al. tested the sensitivity of different indices to the LAI and found that the chlorophyll indices (CIGreen, CIRed-edge and the MERIS terrestrial chlorophyll index, MTCI) show strong and significant linear relationships with the green LAI [9].Gonsamo et al. developed a new vegetation index, the reduced infrared simple ratio (RISR), which is significantly less affected by the soil background reflectance while maintaining a high sensitivity to a wide range of LAI values [10].Kross et al. showed that the normalized difference vegetation index (NDVI), the red-edge NDVI and the green NDVI are sensitive to the LAI but insensitive to the crop type [11].The enhanced vegetation index (EVI) is also commonly used to estimate the LAI [12,13].
However, most of these vegetation indices are also sensitive to the effect of BRDF because of bi-directional reflectance distribution function (BRDF) effects [14,15].To reduce the BRDF effects, many studies have focused on using semi-empirical kernel-based BRDF models to normalize the bidirectional reflectance and then using the BRDF-normalized reflectance or vegetation indices to estimate the LAI [16,17].Multi-angular observations make it possible to capture the uneven scattering of sunlight by vegetation.Some multi-angular observation data are already available from sensors on-board satellites-this includes CHRIS/PROBA, MISR/Terra and MODIS/MCD43 data.The BRDF signature has different effects on different vegetation indices [18,19].Lacaze et al. investigated a directional index named the hotspot-dark-spot index (HDS), which was defined as the normalized difference between the reflectance at the hot spot and dark spot and was used to discriminate vegetation types using only a portion of the BRDF [20].It has been reported that the HDS can adequately represent the geometric structure of vegetation.Hasegawa et al. proposed the normalized hotspot-signature vegetation index (NHVI), which was defined as the product of the NDVI and the HDS [21].The NHVI was considered to give a better quantitative estimation of the LAI than the NDVI, especially in boreal forest.However, many BRF data have only one view angle, and it is difficult to acquire the multi-angular BRF reflectance at the hotspot and dark spot.
There are other methods for LAI estimation that are based on physical models of the BRDF.For example, the MODIS LAI product is based on the three-dimensional radiative transfer equation with surface spectral bi-directional reflectance factors (BRFs) as input parameters [22].However, these physical models are complicated, and a number of parameters are needed as inputs.As a result, they are time-consuming to use [23].This may lead to unreliable estimates of biophysical parameters.As a result, the vegetation index approach is still widely used for LAI estimation.
By combining the informaiton in some specific bands, the interference to vegetation indices from other factors (such as atmosphere and soil) could be corrected.For example, in the atmospherically resistant vegetation index (ARVI) [24], the red band used in the conventional NDVI was replaced by a combination of the red and blue bands.This produced an index that was able to self-correct for atmospheric effects.In the green atmospherically resistant index (GARI), the difference between the blue band and red band was used to self-correct for the effects of the atmosphere on the green band in the NDVI [25].
Here, we aim to develop a new BRDF-resistant vegetation index (BRVI) that is sensitive to the LAI but insensitive to the effect of BRDF.By including these features, the BRVI is expected to be able to give accurate estimates of the LAI, even without correction for the BRDF.

BRDF Simulation Experiments
In order to model and validate the new vegetation index, we carried out a set of simulation experiments.In this study, the PROSAIL model, PROSPECT (Leaf Optical Properties Spectra) leaf model [26] + SAIL (Scattering from Arbitrarily Inclined Leaves) canopy reflectance model [27], was employed for the simulation of multi-angular reflectance with different values of the LAI.
To acquire multi-angular canopy spectra, the ranges of values for the PROSAIL input parameters, listed in Table 1, were chosen based on the Leaf Optical Properties Experiment 93 (LOPEX'93) data set [28], the results of other studies such as Haboudane et al. 2002;Dennett and Ishag, 1998;Pinheiro et al. 2005 andVile et al. 2005 [29][30][31][32], and field measurements.The purpose of our study was to develop a new BRDF-resistant vegetation index which can be used to retrieve the LAI with a reduced influence by the BRDF effects; accordingly, the other fixed parameters in the PROSAIL model were set to their default values with the exception of the LAI, the view zenith angle (VZA), the solar zenith angle (SZA) and the azimuthal difference between the solar and observation angles (RAA).

Field Experiments
For the purpose of further evaluation of the proposed BRDF-resistant vegetation index (BRVI), a series of field experiments were also carried out.The evaluation of the performance of the BRVI included two aspects: first, the correlation between the BRVI and the LAI and, second, the variation in the BRVI in multi-angular observations.Data from two sets of field experiments were, therefore, used to investigate these two aspects.

Winter Wheat Spectra Experiments
The first set of field experiments was conducted in 2002 at the China National Experimental Station for Precision Agriculture in Xiaotangshan town, Changping District, Beijing (40 • 11 N, 116 • 27 E).Three winter wheat varieties, Jingdong 8, Jing 9428 and Zhongyou 9507, were planted in 48 plots, each with a plot size of 32.4 m × 30.0 m.For each winter wheat variety, there were four irrigation treatments: no irrigation (W0), and irrigation with water at rates of 225 m 3 /ha (W1), 450 m 3 /ha (W2) (typical application) and 675 m 3 /ha (W3) (over-irrigated).There were also four fertilization treatments: no fertilizer (N0), and nitrogen applied at rates of 50 kg/ha (N1), 200 kg/ha (N2) (typical application) and 350 kg/ha (N3) (excessively fertilized) at both the reviving stage and jointing stage.Consequently, the LAI of the wheat was expected to be different in each plot.In this experiment, winter wheat was planted at the end of September.Winter wheat starts to regrow at the reviving stage in March and matures in mid-June.We obtained multiple LAI datasets for four growth stages, as shown in Table 2, before the wheat reached the heading stage.The LAI was measured using a destructive sampling method.All the winter wheat plants within an area of 60 cm × 60 cm were sampled in the laboratory using the specific leaf weight method [33].
The canopy reflectance spectra of the winter wheat in the 48 plots were measured by an ASD spectrometer (Analytical Spectral Devices, Boulder, CO, USA) between 10:00 and 14:00 (Beijing local time) under clear-sky conditions.The spectra acquired by the ASD spectrometer covered the spectral range of 400-2500 nm with a spectral resolution of 3 nm, an interval of 1.4 nm (350-1050 nm) and a field of view of 25 • .A BaSO4 calibration panel (0.4 m × 0.4 m) was used to determine the reflectance spectrum.The background was normal farmland soil.

Multi-Angular Spectral Measurements
To evaluate how the BRVI varied with different effect of BRDF parameters, a set of multi-angular spectral measurements was also carried out.
The winter wheat canopy BRDF spectra were measured by an ASD spectrometer from 8:55 to 16:22 on 9 April 2016 using an automatic multi-angle observation system (MAOS, Figure 1).The MAOS can automatically collect the canopy radiance at different viewing angles together with the reference radiance [34].The system was composed of a two-dimensional automatic goniometer and a spectrometer.It was controlled by a laptop and could be moved to different zenith and view azimuth angles independently.In the diurnal experiments, we set an increment of ten degrees in the view zenith direction for the solar principal plane and the cross-principal plane; around the hotspot position in the solar principal plane, a higher density of measurements-with a sampling interval of two degrees-was made in order to capture the rapid change in the BRDF there.As a result, 15 or 16 measurements of directional spectral radiance were collected within the range −60 • to +60 • for the view zenith angle in the solar principal plane and cross the solar principal plane respectively.Near to the hotspot direction, the MAOS horizon rod, which was used to fix the sensor, was projected into the field of view of the spectrometer.This resulted in self-shadowing that was 2 cm wide.As the horizontal rod was fixed at 2.0 m above the canopy, the shadow of the rod covered an area of about 89 cm 2 in the 25 • field of view (FOV).However, the self-shadowing covered only about 1.44% of the total FOV and so this effect was not taken into account in this study.The LAI was measured using a destructive sampling method.All the winter wheat plants within an area of 60 cm × 60 cm were sampled in the laboratory using the specific leaf weight method [33].
The canopy reflectance spectra of the winter wheat in the 48 plots were measured by an ASD spectrometer (Analytical Spectral Devices, Boulder, CO, USA) between 10:00 and 14:00 (Beijing local time) under clear-sky conditions.The spectra acquired by the ASD spectrometer covered the spectral range of 400-2500 nm with a spectral resolution of 3 nm, an interval of 1.4 nm (350-1050 nm) and a field of view of 25°.A BaSO4 calibration panel (0.4 m × 0.4 m) was used to determine the reflectance spectrum.The background was normal farmland soil.

Multi-Angular Spectral Measurements
To evaluate how the BRVI varied with different effect of BRDF parameters, a set of multi-angular spectral measurements was also carried out.
The winter wheat canopy BRDF spectra were measured by an ASD spectrometer from 8:55 to 16:22 on 9 April 2016 using an automatic multi-angle observation system (MAOS, Figure 1).The MAOS can automatically collect the canopy radiance at different viewing angles together with the reference radiance [34].The system was composed of a two-dimensional automatic goniometer and a spectrometer.It was controlled by a laptop and could be moved to different zenith and view azimuth angles independently.In the diurnal experiments, we set an increment of ten degrees in the view zenith direction for the solar principal plane and the cross-principal plane; around the hotspot position in the solar principal plane, a higher density of measurements-with a sampling interval of two degrees-was made in order to capture the rapid change in the BRDF there.As a result, 15 or 16 measurements of directional spectral radiance were collected within the range −60° to +60° for the view zenith angle in the solar principal plane and cross the solar principal plane respectively.Near to the hotspot direction, the MAOS horizon rod, which was used to fix the sensor, was projected into the field of view of the spectrometer.This resulted in self-shadowing that was 2 cm wide.As the horizontal rod was fixed at 2.0 m above the canopy, the shadow of the rod covered an area of about 89 cm 2 in the 25° field of view (FOV).However, the self-shadowing covered only about 1.44% of the total FOV and so this effect was not taken into account in this study.

BRDF Characteristics of Canopy Spectra at Different Bands
In order to develop a BRDF-resistant vegetation index (BRVI), the BRDF effect at different bands was first investigated in the experiment described above.The bi-directional reflectance of the canopy spectra at different bands for different values of the LAI (0.5, 1, 2, 3, 4 and 6) in the solar principal plane with an SZA of 30 • is illustrated in Figure 2.

BRDF Characteristics of Canopy Spectra at Different Bands
In order to develop a BRDF-resistant vegetation index (BRVI), the BRDF effect at different bands was first investigated in the experiment described above.The bi-directional reflectance of the canopy spectra at different bands for different values of the LAI (0.5, 1, 2, 3, 4 and 6) in the solar principal plane with an SZA of 30° is illustrated in Figure 2. As shown in Figure 2, the shape of the BRDF curve at the green band seems to be similar to that at the NIR band.Meanwhile, the shape of the BRDF curve at the blue band seems to be similar to that at the red band.
To quantitatively evaluate the similarity of the shape of BRDF curve, we employed the Modified Rahaman-Pinty-Verstraete (MRPV) surface reflectance model to express the bidirectional reflectance at blue, green, red and NIR bands for different LAI values.MRPV model uses three functions to express the bidirectional reflectance ( ) of a surface illuminated from a direction ( , ) and observed from a direction( , ). ρ ( , , , ) = × ( , ; ) × ( ; ) × ( ; G), k, bM, and ρc are the three fitting parameters of the MRPV model.k is the parameter describing the bowl (k < 1) or bell (k > 1) shape of the BRDF curve of the vegetation canopy, bM is the parameter describing the BRDF asymmetry, and ρc is the hotspot parameter [35][36][37].As the hotspot parameter (ρc) is calculated directly using the average of the canopy reflectance, we used the k and bM for comparing the BRDF curve shape at different bands.Engelsen et al. [37] listed the fitted parameters of MRPV model using various datasets at visible and NIR bands.The results (statistics on corn, grassland and wheat) show that, for the visible bands, the range of k is 0.811 to 0.998, the range of bM is −0.535 to 0.084; for the NIR band, the range of k is 0.647 to 0.819, the range of bM is −0.266 to 0.078.As shown in Figure 2, the shape of the BRDF curve at the green band seems to be similar to that at the NIR band.Meanwhile, the shape of the BRDF curve at the blue band seems to be similar to that at the red band.
To quantitatively evaluate the similarity of the shape of BRDF curve, we employed the Modified Rahaman-Pinty-Verstraete (MRPV) surface reflectance model to express the bidirectional reflectance at blue, green, red and NIR bands for different LAI values.MRPV model uses three functions to express the bidirectional reflectance (ρ) of a surface illuminated from a direction (θ 1 , ϕ 1 ) and observed from a direction (θ 2 , ϕ 2 ).
k, b M , and ρ c are the three fitting parameters of the MRPV model.k is the parameter describing the bowl (k < 1) or bell (k > 1) shape of the BRDF curve of the vegetation canopy, b M is the parameter describing the BRDF asymmetry, and ρ c is the hotspot parameter [35][36][37].As the hotspot parameter (ρ c ) is calculated directly using the average of the canopy reflectance, we used the k and b M for comparing the BRDF curve shape at different bands.Engelsen et al. [37] listed the fitted parameters of MRPV model using various datasets at visible and NIR bands.The results (statistics on corn, grassland and wheat) show that, for the visible bands, the range of k is 0.811 to 0.998, the range of b M is −0.535 to 0.084; for the NIR band, the range of k is 0.647 to 0.819, the range of b M is −0.266 to 0.078.So, for the visible bands, the variation of asymmetry is relatively greater than that at the NIR band (with larger range of b M ).The BRDF curve shape at the NIR band is dominated by the bowl-shape parameter (k).
The fitted values of the parameters, together with the fitting correlation coefficient (R) and root mean square error (RMSE) are listed in Table 3.According to the results listed in Table 3, when the LAI is relatively low, the asymmetry of BRDF of the four bands are similar (with similar values of b M ), but considering the bowl-shape parameters (k), the green band shows similarity to the NIR band, and the red band shows similarity to the blue band; when the LAI is relatively high, the bowl-shape of the BRDF curve at the four bands turns to be similar (with similar values of k), but asymmetry at the NIR band seems to be similar to that at the green band (with similar values of b M ), and the red band also shows similarity to the blue band.
Reflectance of the surface underneath the vegetation canopy also has an impact on the BRDF shape.At the NIR band, the reflectance of vegetation is much higher than that of the soil, so the influence of soil background to the BRDF shape at NIR band is not significant.For farmland soil, the reflectance is relatively low, so the influence of soil background is relatively small, even for the visible bands.Actually, we have tested the performance of simulated datasets with both dry (high reflectance) and wet (low reflectance) soil background, and the conclusions are similar.So the conclusions about the similarity in BRDF curve shape of NIR-green band and red-blue band are reasonable to be expected valid for most farmland conditions.Limited by the length of the paper, the results of simulation with wet soil background (for which the impact on the canopy BRDF curve shape is relatively small) are not included in the main text (see Appendix A).
We also analyzed the BRDF characteristics of different bands for different SZA values.We selected datasets with an LAI of 4 and solar zenith angles of 10 • to 60 • (with 10 • steps).The bi-directional reflectances within the solar principle plane are plotted in Figure 3.We still used parameters k and b M in the MRPV model to express the bidirectional reflectance at blue, green, red and NIR bands for different SZA values.The two parameters were listed in Table 4. Considering the curve shape, parameter k and b M comprehensively, the shape of the BRDF curve at the green band seems to be similar to that at the NIR band for all the SZA conditions that were tested.In addition, the shape of the BRDF curve at the blue band was similar to that at the red band.The results are also consistent with Figure 2.
Remote Sens. 2016, 8, 947 7 of 18 We also analyzed the BRDF characteristics of different bands for different SZA values.We selected datasets with an LAI of 4 and solar zenith angles of 10° to 60° (with 10° steps).The bidirectional reflectances within the solar principle plane are plotted in Figure 3.We still used parameters k and bM in the MRPV model to express the bidirectional reflectance at blue, green, red and NIR bands for different SZA values.The two parameters were listed in Table 4. Considering the curve shape, parameter k and bM comprehensively, the shape of the BRDF curve at the green band seems to be similar to that at the NIR band for all the SZA conditions that were tested.In addition, the shape of the BRDF curve at the blue band was similar to that at the red band.The results are also consistent with Figure 2.

Development of the BRDF-Resistant Vegetation Index (BRVI)
The consistency in the shape of the BRDF across different bands, as illustrated in Figures 2  and 3, was employed to develop the new BRDF-resistant vegetation index for estimating the LAI.The reflectance ratios of NIR band to green band and blue band to red band were reasonably assumed to be resistant to the BRDF effects.This idea was inspired by the development of other vegetation indices such as ARVI [24] and GARI [25].
We developed a new BRDF-resistant vegetation index, which is a normalized vegetation index sensitive to the LAI but insensitive to the effect of BRDF.The index is shown as Equation (2).For the BRVI to be sensitive to the LAI, Ratio 1 should increase for larger values of the LAI while Ratio 2 should decrease.
The ratio of the NIR to green reflectances could thus possibly be used for Ratio 1 while the ratio of blue to red reflectances could be used for Ratio 2 .
To determine the form of the BRVI, we examined the variation in the reflectance at four bands for different solar zenith angles.We used a simulated dataset with the LAI set to 4, an SZA of 10 • to 60 • (with 10 • steps) and a VZA that varied between −60 • and +60 • in the solar principal plane.The directional ratio (DR) of the maximum to minimum reflectance for each multi-angular spectral measurement in the solar principal plane [38] was employed to quantitatively calculate the variation in the reflectance at four bands for different solar zenith angles.The situation where the view zenith angle exactly matches the solar zenith angle in the solar principal plane rarely occurs for real measurements.We therefore excluded the data in the hotspot direction when calculating the DR for each band.The DR values are shown in Table 5.As the results show, the BRDF has a larger effect at the green band than that at the NIR band; the effect at the red band is also larger than that at the blue band.In other words, the variation in the amplitude of the bi-directional reflectance in the solar principal plane is different for different bands.Therefore, we could not simply combine the NIR to green and blue to red reflectance ratios in order to develop the new vegetation index.
By examining the shape of the BRDF curve for the four bands shown in Figures 2 and 3, we found that the effect of BRDF on the NIR and green bands were opposite to those on the red and blue bands.As a result, we optimized the resistance of the index to the BRDF effects by introducing the red and green band, still using a self-correcting approach.The divisors in the two simple reflectance ratios were improved by combining the reflectance at the red and green bands.The new BRVI was defined as a normalized combination of the two improved reflectance ratios: We allowed the coefficients for the red and blue bands to vary and optimized them so that the change in the BRVI with changes in the VZA and SZA was minimal based on a simulated dataset.After optimization to minimize the BRDF effects, we confirmed the form of the BRVI with a k 1 of 0.1 and a k 2 of 0.5.

Comparison of the Resistance to the BRDF Effects between BRVI and Other Vegetation Indices
In this section, the resistance of the BRVI to the BRDF effects is compared with that of other vegetation indices (listed in Table 6) using both simulated and experimental datasets.We selected simulated datasets with an LAI of 4 and SZA ranging from 10 • to 60 • .The VZA varied between −60 • and +60 • in the solar principal plane.The resistance to BRDF signatures of four commonly used vegetation indices (Table 6) were chosen for comparison with the BRVI.The different vegetation indices were calculated for each VZA and SZA and the values then were normalized by dividing the maximum values at different view zenith angle.In practice, there are rarely observations just at the hotspot direction.In addition, all spectrometers have a field of view.For example, the field of view of the ASD Spectrometer used in this study is 25 degrees.The hotspot effect will be smoothed within the field of view.So the values at hotspots were not included in this part to make the results more representative for practice.The results are plotted in Figure 4.In addition, for a quantitative comparison, we also calculated the DR for each vegetation index in order to evaluate the directional variation in amplitude for each index for each multi-angular spectral measurement in the solar principal plane, as shown in Table 7.
Remote Sens. 2016, 8, 947 9 of 18 After optimization to minimize the BRDF effects, we confirmed the form of the BRVI with a k1 of 0.1 and a k2 of 0.5.

Comparison of the Resistance to the BRDF Effects between BRVI and Other Vegetation Indices
In this section, the resistance of the BRVI to the BRDF effects is compared with that of other vegetation indices (listed in Table 6) using both simulated and experimental datasets.

Validation Using Simulated Dataset
We selected simulated datasets with an LAI of 4 and SZA ranging from 10° to 60°.The VZA varied between −60° and +60° in the solar principal plane.The resistance to BRDF signatures of four commonly used vegetation indices (Table 6) were chosen for comparison with the BRVI.The different vegetation indices were calculated for each VZA and SZA and the values then were normalized by dividing the maximum values at different view zenith angle.In practice, there are rarely observations just at the hotspot direction.In addition, all spectrometers have a field of view.For example, the field of view of the ASD Spectrometer used in this study is 25 degrees.The hotspot effect will be smoothed within the field of view.So the values at hotspots were not included in this part to make the results more representative for practice.The results are plotted in Figure 4.In addition, for a quantitative comparison, we also calculated the DR for each vegetation index in order to evaluate the directional variation in amplitude for each index for each multi-angular spectral measurement in the solar principal plane, as shown in Table 7.As Figure 4 and Table 7 show, the BRVI is more resistant to the BRDF effects than other vegetation indices and has a much smoother BRDF curve and smaller DR values, especially for SZA values less than 30 • .For larger values of the SZA, the BRVI and NDVI have a similar resistance to the BRDF; however, the effect on the BRVI is still smaller than the effect on the four other vegetation indices.

Validation Using Field Dataset
Multi-angular spectral experimental data, obtained between 8:08 and 17:00 on 9 April 2016, were used to test the resistance of five vegetation indices to the BRDF effects.Five vegetation indices were calculated using these datasets cross the solar principal plane and in the solar principal plane and the values then were normalized by dividing the maximum values at different view zenith angle.The results for different view zenith angle are plotted in Figure 5 (cross the solar plane) and Figure 6 (in the solar plane).Meanwhile, the DR values of the five vegetation indices were calculated and are listed in Table 8 (cross the solar principal plane) and  As Figure 4 and Table 7 show, the BRVI is more resistant to the BRDF effects than other vegetation indices and has a much smoother BRDF curve and smaller DR values, especially for SZA values less than 30°.For larger values of the SZA, the BRVI and NDVI have a similar resistance to the BRDF; however, the effect on the BRVI is still smaller than the effect on the four other vegetation indices.

Validation Using Field Dataset
Multi-angular spectral experimental data, obtained between 8:08 and 17:00 on 9 April 2016, were used to test the resistance of five vegetation indices to the BRDF effects.Five vegetation indices were calculated using these datasets cross the solar principal plane and in the solar principal plane and the values then were normalized by dividing the maximum values at different view zenith angle.The results for different view zenith angle are plotted in Figure 5 (cross the solar plane) and Figure 6 (in the solar plane).Meanwhile, the DR values of the five vegetation indices were calculated and are listed in Table 8 (cross the solar principal plane) and Table 9 (in the solar principal plane).As shown in Figures 5 and 6 and Tables 8 and 9, both for the values cross the solar principal plane and in the solar principal plane, the BRVI was least affected by the BRDF effects because the shape of the BRDF curve for the BRVI is smoother than for the other vegetation indices and the DR value is much smaller.As the SZA decreases, the effect of the BRDF on the vegetation indices becomes larger, except in the case of the SAVI.In the solar principal plane, the effect of the BRDF on the vegetation indices is larger than it is cross the solar principal plane due to the influence of the hotspot effect.
The above analysis thus shows that, for both the simulated and field datasets, the BRVI is less sensitive to the effect of the BRDF than the four other vegetation indices that were studied.As shown in Figures 5 and 6 and Tables 8 and 9, both for the values cross the solar principal plane and in the solar principal plane, the BRVI was least affected by the BRDF effects because the shape of the BRDF curve for the BRVI is smoother than for the other vegetation indices and the DR value is much smaller.As the SZA decreases, the effect of the BRDF on the vegetation indices becomes larger, except in the case of the SAVI.In the solar principal plane, the effect of the BRDF on the vegetation indices is larger than it is cross the solar principal plane due to the influence of the hotspot effect.
The above analysis thus shows that, for both the simulated and field datasets, the BRVI is less sensitive to the effect of the BRDF than the four other vegetation indices that were studied.

Performance of BRVI for Estimation of LAI
To test the performance of the BRVI for estimation of the LAI, both simulated and field datasets were used.We compared the relationship between each of the BRVI, NDVI, SAVI, EVI and SR and the LAI.

Validation Using Simulated Dataset
We selected a simulated dataset that had LAI values of 0.2 to 3, an SZA of 30 • , and view zenith angles of −60 • , −50 • , −30 • , 0 • , +60 • +50 • and + 30 • in the solar principal plane to evaluate the sensitivity of the BRVI to the LAI when the BRDF effect was included.The VZA increment was set as 1 • for the PROSAIL simulation and the simulated data were averaged in order to match the 25 • field-of-view of the ASD spectrometer.We analyzed the relationship between the vegetation indices and the LAI with the BRDF included.The resulting scatter diagrams are depicted in Figure 7.

Performance of BRVI for Estimation of LAI
To test the performance of the BRVI for estimation of the LAI, both simulated and field datasets were used.We compared the relationship between each of the BRVI, NDVI, SAVI, EVI and SR and the LAI.

Validation Using Simulated Dataset
We selected a simulated dataset that had LAI values of 0.2 to 3, an SZA of 30°, and view zenith angles of −60°, −50°, −30°, 0°, +60° +50 ° and + 30° in the solar principal plane to evaluate the sensitivity of the BRVI to the LAI when the BRDF effect was included.The VZA increment was set as 1° for the PROSAIL simulation and the simulated data were averaged in order to match the 25° field-of-view of the ASD spectrometer.We analyzed the relationship between the vegetation indices and the LAI with the BRDF included.The resulting scatter diagrams are depicted in Figure 7. Figure 7 shows that the when trying to estimate the LAI using the vegetation indexes, the index that was least affected by the BRDF was the BRVI.It can also be seen that the BRVI performed well for LAI retrieval, with a coefficient of determination (R 2 ) of 0.97; it also had the lowest RMSE of 0.25.It is worth noting that the BRVI performs best at LAI estimation when the LAI is low.For higher values of the LAI, the performance of the NDVI, EVI and SAVI is good too.This may because, due to the saturation phenomenon, the BRDF effects are not obvious for vegetation with a relatively high LAI.

Validation Using Field Dataset
As mentioned in Section 2.2, we obtained LAI and canopy reflectance data at four different growth stages in 2002.Vegetation indices were calculated using the data for all four growth stages together and their relationship with the LAI then evaluated.Plots of the LAI against the different vegetation indexes are depicted in Figure 8.  Figure 7 shows that the when trying to estimate the LAI using the vegetation indexes, the index that was least affected by the BRDF was the BRVI.It can also be seen that the BRVI performed well for LAI retrieval, with a coefficient of determination (R 2 ) of 0.97; it also had the lowest RMSE of 0.25.It is worth noting that the BRVI performs best at LAI estimation when the LAI is low.For higher values of the LAI, the performance of the NDVI, EVI and SAVI is good too.This may because, due to the saturation phenomenon, the BRDF effects are not obvious for vegetation with a relatively high LAI.

Validation Using Field Dataset
As mentioned in Section 2.2, we obtained LAI and canopy reflectance data at four different growth stages in 2002.Vegetation indices were calculated using the data for all four growth stages together and their relationship with the LAI then evaluated.Plots of the LAI against the different vegetation indexes are depicted in Figure 8.
To test the performance of the BRVI for estimation of the LAI, both simulated and field datasets were used.We compared the relationship between each of the BRVI, NDVI, SAVI, EVI and SR and the LAI.

Validation Using Simulated Dataset
We selected a simulated dataset that had LAI values of 0.2 to 3, an SZA of 30°, and view zenith angles of −60°, −50°, −30°, 0°, +60° +50 ° and + 30° in the solar principal plane to evaluate the sensitivity of the BRVI to the LAI when the BRDF effect was included.The VZA increment was set as 1° for the PROSAIL simulation and the simulated data were averaged in order to match the 25° field-of-view of the ASD spectrometer.We analyzed the relationship between the vegetation indices and the LAI with the BRDF included.The resulting scatter diagrams are depicted in Figure 7. Figure 7 shows that the when trying to estimate the LAI using the vegetation indexes, the index that was least affected by the BRDF was the BRVI.It can also be seen that the BRVI performed well for LAI retrieval, with a coefficient of determination (R 2 ) of 0.97; it also had the lowest RMSE of 0.25.It is worth noting that the BRVI performs best at LAI estimation when the LAI is low.For higher values of the LAI, the performance of the NDVI, EVI and SAVI is good too.This may because, due to the saturation phenomenon, the BRDF effects are not obvious for vegetation with a relatively high LAI.

Validation Using Field Dataset
As mentioned in Section 2.2, we obtained LAI and canopy reflectance data at four different growth stages in 2002.Vegetation indices were calculated using the data for all four growth stages together and their relationship with the LAI then evaluated.Plots of the LAI against the different vegetation indexes are depicted in Figure 8.It can be seen that the correlation between the BRVI and LAI was similar to that between the NDVI, SAVI and EVI, and the LAI with the coefficients of determination (R 2 ) all being higher than 0.8.The RMSE for the estimation of the LAI based on the BRVI is 0.83.Because of this significant correlation, the BRVI could, therefore, be used to estimate the LAI.

Discussion
Compared with traditional vegetation indices (NDVI, EVI, SAVI and SR), BRVI is sensitive to LAI but insensitive to the effect of BRDF.It is worth noting that the BRVI performs best at LAI estimation when the LAI is low.For higher LAI values, the performance of the NDVI, EVI and SAVI is also good as, due to the saturation phenomenon, the effects of the BRDF are not obvious for vegetation with a relatively high LAI.In practice, the situation where the LAI is low and the BRDF effect is strong is, in fact, of more interest.Therefore, the new BRVI has good potential for estimation of the LAI.
Although our results showed that the BRVI is less sensitive to BRDF effects, the new index still has some limitations.The reflectance of farmland soil is relatively low.The background in the simulation experiments used in this paper is dry farmland soil with relatively high reflectance.We also tested the performance of wet farmland soil which has relatively low reflectance.For wet farmland soil, the shape of the BRDF curve at the green band is also similar to that at the NIR band.Meanwhile, the shape of the BRDF curve at the blue band is also similar to that at the red band.As a result, BRVI is suitable to retrieve LAI of crops.However, for other soil types which have higher reflectance, soil scattering dominates much more over leaf scattering.Our BRVI may not effectively eliminate the BRDF effect of LAI retrieval.
The leaf chlorophyll content varies with the growing stage, and has an influence on the leaf albedo.However, in the field experiments in 2002, we found that the leaf chlorophyll content of winter wheat does not vary much (39.53 to 49.72 µg/cm 2 ) during the growth stages we studied.So in this paper, the leaf chlorophyll content was fixed.Actually, we tested the BRDF characteristics using a simulated dataset with relatively low and high values of chlorophyll content, and the results show that the similarity of BRDF characteristics reported in this paper is valid for different chlorophyll content levels.
The development of the BRVI was based on simulations and observations, which meant that only multi-angular spectra of crop were analyzed.However, the BRDF characteristics of other vegetation types could be quite different from those of simple homogenous vegetation and so the approach in this paper might not then be valid.In the case of forests, for example, the most important factor effecting the BRDF characteristics is the interaction of photons with the rough surface of the tree crowns and with the soil in between crown openings.This causes the observed variation in the directional reflectance distribution of plant canopies [43].As a result, we need more simulated data and also field data for different vegetation types in order to verify the ability of the BRVI to estimate the LAI when the BRDF effects are smaller.

Conclusions
Most vegetation indices show obvious BRDF characteristics due to the bi-directional reflectance of the canopy.In this paper, a new BRDF-resistant vegetation index (BRVI), which is sensitive to LAI but insensitive to the effect of BRDF was developed and then compared with traditional vegetation indices (NDVI, EVI, SAVI and SR) to test its sensitivity to BRDF effects and the LAI using both simulated data and in-situ measurements.Simulations for farmland soil conditions and with medium chlorophyll content level of the vegetation indices for various values of the SZA and VZA showed that the influence of the BRDF on the BRVI was weaker than its influence on the other indices.The BRVI also had the lowest DR, which ranged from 1.01 to 1.05-this compared with DR values for the other vegetation indices of as high as 3.20.Similarly, using the field data, the DR values for the other vegetation indices were found to be 1.05 to 3 times larger than that for the BRVI.Simulations with various values of the SZA and VZA also showed that, among the indices tested, the BRVI performed best at LAI retrieval when BRDF effects were included, with the BRVI having the highest R 2 of 0.97 and lowest RMSE of 0.25.In addition, the BRVI performed well at LAI estimation, giving an R 2 of 0.84 and an RMSE of 0.83 for field data with a single view angle.In summary, we developed a new BRDF-resistant vegetation index (BRVI), which is sensitive to the LAI but insensitive to the effect of BRDF for farmland soil conditions.However, the results of this study require further validation using additional datasets.Table A1.Fitting parameters of the MRPV model for the bidirectional reflectance for the 5 multiangular spectral simulations at the blue, green, red and NIR bands for different LAI values (with wet soil background).k is a parameter that describes the bowl (k < 1) or bell (k > 1) shape of the BRDF curve of the vegetation canopy; bM is the parameter for describing the BRDF asymmetry; ρc is the hotspot parameter; R and RMSE are the correlation coefficient and root mean square error of the fitted model, respectively.Table A1.Fitting parameters of the MRPV model for the bidirectional reflectance for the 5 multi-angular spectral simulations at the blue, green, red and NIR bands for different LAI values (with wet soil background).k is a parameter that describes the bowl (k < 1) or bell (k > 1) shape of the BRDF curve of the vegetation canopy; b M is the parameter for describing the BRDF asymmetry; ρ c is the hotspot parameter; R and RMSE are the correlation coefficient and root mean square error of the fitted model, respectively.

Figure 1 .
Figure 1.Photo of instrumentation using the multi-angle observation system (MAOS).Figure 1. Photo of instrumentation using the multi-angle observation system (MAOS).

Figure 1 .
Figure 1.Photo of instrumentation using the multi-angle observation system (MAOS).Figure 1. Photo of instrumentation using the multi-angle observation system (MAOS).

Figure 2 .
Figure 2. The simulated bi-directional reflectance in the solar principle plane at the blue, green, red and the near-infrared (NIR) bands for different leaf area index (LAI) values.A view zenith angle of 0° indicates a direction vertically downwards from the observation point.Negative values of the view zenith angle (VZA) indicate the backwards direction with a relative azimuth angle of 0°; positive values of the VZA indicate the forwards direction with a relative azimuth angle of 180°.A peak in reflectance, known as the hotspot, is assumed to occur for the negative VZA value equal to each corresponding solar zenith angle (SZA) value.

Figure 2 .
Figure 2. The simulated bi-directional reflectance in the solar principle plane at the blue, green, red and the near-infrared (NIR) bands for different leaf area index (LAI) values.A view zenith angle of 0 • indicates a direction vertically downwards from the observation point.Negative values of the view zenith angle (VZA) indicate the backwards direction with a relative azimuth angle of 0 • ; positive values of the VZA indicate the forwards direction with a relative azimuth angle of 180 • .A peak in reflectance, known as the hotspot, is assumed to occur for the negative VZA value equal to each corresponding solar zenith angle (SZA) value.

Figure 3 .
Figure 3.The simulated bi-directional reflectances in the solar principle plane at the blue, green, red and NIR bands for different SZA values.

Figure 3 .
Figure 3.The simulated bi-directional reflectances in the solar principle plane at the blue, green, red and NIR bands for different SZA values.

Figure 4 .
Figure 4. Different vegetation indices calculated using simulated bi-directional reflectances in the solar principle plane for different SZA values.

igure 4 .
Different vegetation indices calculated using simulated bi-directional reflectances in the solar principle plane for different SZA values.

Figure 5 .
Figure 5.The values of five different vegetation indices for different SZA values calculated from canopy bi-directional reflectances of winter wheat cross the solar principle plane, acquired from 8:08 to 17:00 on 9 April 2016.

Figure 5 .
Figure 5.The values of five different vegetation indices for different SZA values calculated from canopy bi-directional reflectances of winter wheat cross the solar principle plane, acquired from 8:08 to 17:00 on 9 April 2016.

Figure 6 .
Figure 6.Different vegetation indices after normalizing each vegetation indices (dividing the maximum value in the different view zenith angle) calculated by canopy bi-directional reflectances of winter wheat in the solar principle plane under different SZA values, acquired from 8:08 to 17:00 on 9 April 2016.

Figure 7 .
Figure 7. Relationship between LAI and vegetation indices.LAI predictive function for VZA varying from −60° to 60° and SZA of 30.

Figure 7 .
Figure 7. Relationship between LAI and vegetation indices.LAI predictive function for VZA varying from −60 • to 60 • and SZA of 30.

Figure 7 .
Figure 7. Relationship between LAI and vegetation indices.LAI predictive function for VZA varying from −60° to 60° and SZA of 30.

Figure 8
Figure8illustrates the relationship between the BRVI and four other vegetation indices, and the LAI.It can be seen that the correlation between the BRVI and LAI was similar to that between the NDVI, SAVI and EVI, and the LAI with the coefficients of determination (R 2 ) all being higher than 0.8.The RMSE for the estimation of the LAI based on the BRVI is 0.83.Because of this significant correlation, the BRVI could, therefore, be used to estimate the LAI.

Figure A2 .
Figure A2.The simulated bi-directional reflectance in the solar principle plane at the blue, green, red and NIR bands for different LAI values (with wet soil background).A view zenith angle of 0° indicates a direction vertically downwards from the observation point.Negative values of the VZA indicate the backwards direction with a relative azimuth angle of 0°; positive values of the VZA indicate the forwards direction with a relative azimuth angle of 180°.A peak in reflectance, known as the hotspot, is assumed to occur for the negative VZA value equal to each corresponding SZA value.

Figure A2 .
Figure A2.The simulated bi-directional reflectance in the solar principle plane at the blue, green, red and NIR bands for different LAI values (with wet soil background).A view zenith angle of 0 • indicates a direction vertically downwards from the observation point.Negative values of the VZA indicate the backwards direction with a relative azimuth angle of 0 • ; positive values of the VZA indicate the forwards direction with a relative azimuth angle of 180 • .A peak in reflectance, known as the hotspot, is assumed to occur for the negative VZA value equal to each corresponding SZA value.

Table 1 .
Parameters used in simulating canopy bi-directional reflectance with the PROSAILmodel, PROSPECT (Leaf Optical Properties Spectra) leaf model + SAIL (Scattering from Arbitrarily Inclined Leaves) canopy reflectance model.

Table 2 .
Dates corresponding to the different growth stages of winter wheat for the experiment carried out in 2002.

Table 2 .
Dates corresponding to the different growth stages of winter wheat for the experiment carried out in 2002.

Table 3 .
Fitting parameters of the Modified Rahaman-Pinty-Verstraete (MRPV) model for the bidirectional reflectance for the five multi-angular spectral simulations at the blue, green, red and NIR bands for different LAI values (with dry soil background).k is a parameter that describes the bowl (k < 1) or bell (k > 1) shape of the bi-directional reflectance distribution function (BRDF) curve of the vegetation canopy; b M is the parameter for describing the BRDF asymmetry; ρ c is the hotspot parameter; R and RMSE are the correlation coefficient and root mean square error of the fitted model, respectively.

Table 4 .
Fitting parameters of the MRPV model for the bidirectional reflectance for the five multiangular spectral simulations at the blue, green, red and NIR bands for different SZA values.

Table 4 .
Fitting parameters of the MRPV model for the bidirectional reflectance for the five multi-angular spectral simulations at the blue, green, red and NIR bands for different SZA values.

Table 5 .
DR values for simulated bi-directional reflectances in the solar principle plane at the blue, green, red and NIR bands for different SZA.

Table 6 .
Multispectral indices employed for retrieving LAI.

Table 6 .
Multispectral indices employed for retrieving LAI.

Table 7 .
Directional ratio (DR) values of different vegetation indices calculated from simulated bi-directional reflectances in the solar principle plane for different SZA values.

Table 9 (
in the solar principal plane).

Table 7 .
Directional ratio (DR) values of different vegetation indices calculated from simulated bidirectional reflectances in the solar principle plane for different SZA values.

Table 8 .
DR values of five different vegetation indices for different SZA values calculated from canopy bi-directional reflectances of winter wheat cross the solar principal plane, acquired from 8:08 to 17:00 on 9 April 2016.

Table 8 .
DR values of five different vegetation indices for different SZA values calculated from canopy bi-directional reflectances of winter wheat cross the solar principal plane, acquired from 8:08 to 17:00 on 9 April 2016.

Table 9 .
The DR values of different vegetation indices calculated by canopy bi-directional reflectances of winter wheat under different SZA values in the solar principal plane, acquired from 8:55 to 16:22 on 9 April 2016.

Table 9 .
The DR values of different vegetation indices calculated by canopy bi-directional reflectances of winter wheat under different SZA values in the solar principal plane, acquired from 8:55 to 16:22 on 9 April 2016.