Modifying Geometric-Optical Bidirectional Reflectance Model for Direct Inversion of Forest Canopy Leaf Area Index

Forest canopy leaf area index (LAI) inversion based on remote sensing data is an important method to obtain LAI. Currently, the most widely-used model to achieve forest canopy structure parameters is the Li-Strahler geometric-optical bidirectional reflectance model, by considering the effect of crown shape and mutual shadowing, which is referred to as the GOMS model. However, it is difficult to retrieve LAI through the GOMS model directly because LAI is not a fundamental parameter of the model. In this study, a gap probability model was used to obtain the relationship between the canopy structure parameter nR and LAI. Thus, LAI was introduced into the GOMS model as an independent variable by replacing nR The modified GOMS (MGOMS) model was validated by application to Dayekou in the Heihe River Basin of China. The LAI retrieved using the MGOMS model with optical multi-angle remote sensing data, high spatial resolution images and field-measured data was in good agreement with the field-measured LAI, with an R-square (R) of 0.64, and an RMSE of 0.67. The results demonstrate that the MGOMS model obtained by replacing the canopy structure parameter nR of the GOMS model with LAI can be used to invert LAI directly and precisely.


Introduction
Leaf Area Index (LAI) is defined as the one-sided green leaf area per unit of horizontal ground surface area [1]; it is an important botany parameter in the study of farming [2][3][4][5][6][7], forestry [8,9], climate of meteorology [10], and ecology [11].Forest ecosystems are essential parts of the terrestrial ecosystem, and the LAI inversion of forest canopy is one of the main objects of forest ecosystem research [12][13][14][15].Due to complex terrain, variable weather conditions, limited field-measured data, and uneven vegetation types, among other factors, the inversion and verification of forest LAI are usually difficult.
Methods of obtaining LAI include field measurements, photography simulations [16,17], and model inversion using semi-empirical [8,18,19] or physical models; in particular, model inversion using physical models is the most widely-used method to achieve forest canopy LAI.The earliest practical canopy reflectance model was determined by Suits [20,21].Verhoef [22] then extended the Suits model by considering the influence of the leaf-angle distribution, to create the SAIL model.The SAIL model can simplify the description of canopy structures.Reyna and Bhadwar [23] later considered the influence of specular reflection.Among all the canopy LAI inversion models, the canopy geometric-optical (GO) model has attracted a lot of interest.The earliest GO model was produced in the mid-1970s [24][25][26].The GO model fully considers the macroscopic geometric structure of the objects and assumes that the objects are geometries with a specific arrangement and known geometric shapes and optical properties.According to the geometric optical principles, the GO model can analyze the influences of the interception, shielding of the incident light, and surface reflection, and determine the directional reflectance of the vegetation canopy.Based on the "scene synthesis model", the GO model considers four components, including sunlit leaf or canopy, shaded leaf or canopy, sunlit background, and shaded background, and establishes the bidirectional reflectance distribution function (BRDF) model under different illumination and observation conditions.In 1985 and 1986, Li and Strahler applied the GO model to remote sensing data using a simple cone model instead of the tree canopy [27,28].In the past, many attempts to optimize the GO model have been made.For example, Schaaf et al. [29] built a geometric optical model that considered terrain factors.Qi et al. [30] combined the statistical model with the optical model, which gives a more accurate estimation of LAI.In order to apply the model to the dense forest, earlier studies [31][32][33] proposed the geometric-optical mutual shadowing (GOMS) model which considers the effects of mutual shadowing and the radiation transfer process [34], and expanded it to the thermal radiation [35], making the GOMS model the most representative GO model.
GO models considering the object and the atmospheric scattering are suitable for discontinuous vegetation (e.g., shrubbery, sparse forests, coniferous forests, and orchards) and rough surfaces for which the radiation transfer model is difficult to apply.The Li-Strahler GOMS model is a pixel-scale GO model which considers the effects of topographic factors.This model can successfully establish the relationship between forest structure parameters (e.g., average vegetation coverage, average tree height, crown size) and the canopy bidirectional reflection distribution function, giving the relationship between the canopy structure parameters (e.g., clear height, crown radius, forest canopy distribution) and the canopy reflection characteristics [28,31].LAI is not an independent parameter in GOMS.Hence, in previous studies [36,37], crown center height (h), horizontal radius (R) or other parameters were firstly retrieved in the GOMS model.Then, canopy LAI was obtained by the statistical relationship between the field-measured h, R or other parameters and the field-measured LAI.
However, this calculated process is very complicated and the accuracy of the retrieved LAI depends on the precision of the field-measured parameters.To address this problem, the GOMS model was improved in the present study by introducing LAI into the model to enable direct inversion of LAI via the modified GOMS (i.e., MGOMS) model.The GOMS model consists of the canopy structure parameter nR 2 , a parameter that describes the crown coverage condition in the nadir observation.Since LAI and nR 2 are both functions of area and are not dependent on the viewing angle, in this paper, the gap probability model is used to build the relationship between LAI and nR 2 .Then, LAI is introduced into the GOMS model by replacing nR 2 with LAI.The accuracy of the MGOMS model is demonstrated using field-measured and remote sensing datasets.

Study Area
In our study, the study site Dayekou forest (100°E, 38.6°N ) is located in the Qilian Mountain area, Gansu Province of China.The Heihe Basin is the second largest inland river basin of the arid region in Northwest China, the annual precipitation in the middle valley is 140 mm [38] and mainly occurs in the summer.Its ecological system is particularly fragile.Dayekou is located in the middle valley of the Heihe River basin and is mostly covered by forest and upland meadow.The main vegetation types in Dayekou are Picea crassifolia, shrubland, and upland meadow, and the dominant forest type is P. crassifolia.The MODIS land cover product (MCD12Q1) was used to select forest pixels.The forest pixels are indicated in green in Figure 1a.

Data Foundation
All the data used in this article were obtained from the Heihe Watershed Allied Telemetry Experimental Research (HiWATER), the field-measured data in the B sample plots were collected at the Dayekou Guantan forest station (100°15′E, 38°32′N), Heihe river basin, China [39].

Field-Measured Data
Field measurements were performed during June 2008 in the field measurement sample plots.Tracing Radiation and Architecture of Canopies (TRAC, Natural Resources Canada, Canada Centre for Remote Sensing, Saint-Hubert, QC, Canada) [40] and LAI-2000 Plant Canopy Analyzer (LAI-2000, LI-COR Inc, Lincoln, NE, USA) [41] were used to measure LAI of each sample plot [42].LAI-2000 was used to measure the effective LAI (LAI e ), and TRAC was used to measure both the effective LAI (LAI e ) and the foliage clumping index (Ω).The true LAI of each sample plot was obtained through the equation LAI t = LAI e / Ω .The field-measured LAI was used to validate the inversion efficiency of the model.
In addition, the forest canopy structural parameters of each tree were measured in the B sample plots.The measured geometrical structural parameters included tree horizontal radius of the tree crown (R), tree height (H), clear bole height (h), and the diameter at breast height (DBH).The tree height of each tree in the B sample plots was measured by laser altimeter (TruPulse 200, Laser Technology Inc. (LTI), Norristown, PA, USA).The field-measured data was used to provide prior knowledge of the canopy structure parameters in the GOMS and MGOMS models.
The protocols for each instrument used in the sample plots and the sample plot layouts were described in detail in a previous study [37].

MODIS, MISR, and SPOT Reflectance Data
In the inversion process, optical multi-angle remote sensing data and high spatial resolution images were used to do the research.Moderate-resolution Imaging Spectro Radiometer (MODIS) and Multi-angle Imaging SpectroRadiometer (MISR) reflectance products were used to provide the multi-angle bidirectional reflectance datasets in this study.SPOT-5 data was used to provide the spectral information.
Ten images of MODIS 500-m daily bidirectional surface reflectance data (MOD09GA) spanning a MODIS 16-day observation period (Figure 3), were selected in the validation process.MODIS 500-m 8-day bidirectional surface reflectance data (MOD09A1) and MISR 1-km bidirectional surface reflectance data, spanning 1 May-28 July 2008, were selected to build the BRF datasets.The BRF datasets were constructed as follows.First, the MISR BRF and MODIS BRF data projections were transferred into UTM, WGS84.Second, the spatial resolutions of the MISR and MODIS data were normalized by dividing each of the MISR pixels into four 500 m pixels with the same BRF value and viewing direction.Third, the bidirectional reflectance and geometric information, including sun zenith angle and azimuth, viewing zenith angle and azimuth of each pixel, were extracted to build the BRF dataset both in the red and near-infrared (NIR) bands.The BRF datasets were used as the input parameters of the model inversion procedure.The main data processing procedures to construct the BRF dataset were explained in detail in a previous study [36].The SPOT-5 multi-band image obtained on 10 August 2008 that covers the study site was used to obtain information on the four component spectra (G, C, Z, T).The spatial resolution of this SPOT-5 data was 10 m.The FLASSH model [43] was used for atmospheric correction of the SPOT-5 image.As noted by Fu et al. [37], the incident radiation was inconsistent during the measuring process.
Hence, we cannot obtain the available spectra of the components from field measurements.Using the four components spectra obtained from remote sensing data with high spatial resolution images can reduce the uncertainty due to scaling effects.Moreover, after the atmospheric correction of the SPOT-5 remote sensing images, the incident radiation can be considered consistent within a MODIS pixel.Thus, we chose the SPOT-5 data instead of the field-measured radiation data to provide the required spectral information.

Airborne LiDAR data
LiDAR data were used to support the canopy height model (CHM) in our study.The LiDAR data we used were acquired on 28 June 2008, using a laser scanner (Riegl LMS-Q560, Riegl, Horn, Austria) with a maximum scan angle of ±0.5 mrad from nadir.The raw LiDAR data were processed using Terrasolid software (Terrasolid software (Version 007), Terrasolid Ltd., Helsinki, Finland) to classify the data into two types: ground points and non-ground points.The DEM was built from ground points, and the digital surface model (DSM) was produced from non-ground points.Then, the CHM was created with a resolution of 0.5 m by subtracting the DEM from the DSM.The CHM was used to provide b and nR 2 in Section 4.2.3.

Geometric-Optic Mutual Shadowing Model
The GOMS model is constructed based on the Li-Strahler geometric-optical model, which considers the mutual shadowing of crowns and makes the geometric optic model more suitable for highly dense canopy forests [31].
The assumption of the GOMS model is that the reflectance of a pixel can be modeled as a sum of the reflectances of its individual scene components as weighted by their respective areas within the pixel [28].The model also assumes that the vegetation canopy BRDF characteristics at the pixel scale can be explained by geometric optical principle when the discontinuous three-dimensional geometry in a pixel is illuminated and observed in different directions.The reflections received by the sensors are the ground reflection and the crown reflection in the field of view A ("A" is only the assumption that the area of the field of view is A).Considering the three-dimensional canopy structure parameters, the influence of sky light and multiple scattering, the received signal can be defined as the combination of four area weighted components: where S is bidirectional reflectance; K g , K c , K z , and K t are the proportions of sunlit background, sunlit crown, shaded background, and shaded crown the GO modeled, respectively; and G, C, Z, and T are the contributions of the sunlit background, sunlit crown, shaded background, and shaded crown, respectively [27,44].K g , K c , K z , and K t can be expressed by a combination of the forest canopy structural parameters.We assume that the tree crown shape is ellipsoidal (Figure 4a) and that the forest canopy structure parameters include R (the horizontal radius of an ellipsoidal crown), b (the vertical half axis of an ellipsoidal crown), h (the height at which a crown center is located), and n (the number of crowns per unit area).When the ellipsoid is stretched in the vertical direction and made b equal to R (Figure 4b), then K g , K c , K z , and K t can be expressed as: where and where O(  ,   , ∅  − ∅  ) is the shaded area in Figure 4b.∅  and ∅  are the solar azimuth and view azimuth, but   and   are not the solar zenith angle and view zenith angle yet; they have been revised in the stretching process: where   ′ and   ′ are solar zenith angle and view zenith angle.
We can simplify the expression of the GOMS model using the function below: S = f(  , ∅  ,   , ∅  ,   , ∅  ,  2 , /, ℎ/, ∆ℎ/, , , , ) (10) where S is bidirectional reflectance;  2 derived from  *  2 which contains the crown density information (n), that represents the crown coverage condition per unit area in the nadir observation; b/R is the crown shape parameter which significantly affects the crown coverage density in the nonnadir direction; h/b is the parameter which represents the crown height from the ground and mainly affects the outward width of the hot spot; ∆h is the variance of the h distribution in one pixel; and ∆ℎ/ is the parameter which describes the discrete degree of the crown height distribution and affects the bowl-shape of the BRDF.  and ∅  are the local slope and aspect, respectively [36,37].  , ∅  ,   , and ∅  are the solar zenith angle, solar azimuth, view zenith angle, and view azimuth, respectively.
In Equation (10), there are many input parameters in the GOMS model, but LAI is not an independent parameter in the model.We therefore attempt to modify the GOMS model by introducing LAI into Equation (10), that is, S = f(LAI, …), to allow LAI to be retrieved by the modified GOMS model directly.

Gap Probability Model
A gap probability model is a GO model that uses the average transmission theory or the simplified radiation transfer equation to simulate the distribution of the leaves by calculating radiation attenuation and consequently determining the radiation characteristics of the vegetation canopy.The most improved gap probability model is the Li-Strahler Gap probability model [45], in which LAI is an independent parameter.
When considering a unique canopy: where s is the distance traveled by a photon within the canopy; τ = K(θ) * LAI/D mainly depends on the leaf density and the foliage angle distribution (FAD); K (θ) is the coefficient of leaf angle distribution function; and D is the average depth of the canopy.
Considering the gap probability within the canopy and among the canopies, and assuming that the canopies are non-overlapping, the gap probability of the canopies is: where P(n) is the probability that a photon will pass through n canopies, P(0) is the gap probability among the canopies, and P(s|n) represents the s-distribution for n canopies.

Modifying GOMS Model Using the Gap Probability Model
According to the descriptions above, both LAI and nR 2 are closely related to area and, thus, we first propose that LAI and nR 2 are closely related and then attempt to build the relationship function between them.We find that nR 2 is similar to vegetation coverage (n * πR 2 is the vertical vegetation coverage)-the vegetation coverage is a function of area.Moreover, vegetation coverage and gap probability are intimately related such that their sum is equal to one.Furthermore, LAI is a parameter in the gap probability model.Thus, it is possible to calculate the relationship between LAI and nR 2 using the gap probability model.
LAI and nR 2 are the parameters independent of the view angle.Hence, the function between LAI and nR 2 when the view zenith angle is zero can be applied in all directions.
Vertical gap probability [46] is the probability that considers the gap probability both within and among canopies when the view zenith angle is zero.
In field of view A, when the forest canopy is non-overlapping, P gap can be written as: where P(s|1) represents the s-distribution for a single individual, and P(0) represents the gap probability when the crown canopy is completely opaque.
When the view zenith angle is zero, the forest canopy is uniformly distributed, i.e., Assuming that: where q is the proportion of the light striking the sphere that passes through the canopy [31,45].Thus, the gap probability model is: Based on Nilson's [47] work, the directional gap fraction can be written as: where Ω(θ) is determined by the distribution of the leaves.When the space distribution of the forests is random, the clumping index is 1; when the space distribution is uniform, the clumping index is greater than 1; when the space distribution is clumped, the clumping index is less than 1.In practical situations, the space distribution is always clumped.Therefore, Ω(θ) is usually defined as the clumping index [40].G(θ) is the foliage projection coefficient that describes the foliage angular distribution [48].Miller [49] simplified G(θ) as: When the view zenith angle is zero, the vertical gap probability can be written as follows: Then, the relationship between LAI and nR 2 is established as the following: In Equation (20), LAI is the true LAI (LAI t ) and Ω * LAI t is the effective LAI (LAI e ).Thus, Equation ( 20) can be written as follows: Equation ( 21) successfully introduces LAI into the GOMS model by replacing the structure parameter nR 2 , to create the modified GOMS (MGOMS) model.

Input Parameters in the Inversion Process of GOMS and MGOMS Models
There are four canopy structural parameters ( nR 2 , b/R, h/b, ∆h/b ), four spectral parameters (G, C, Z, T), and six angle parameters (θ i , ∅ i , θ v , ∅ v , θ s , ∅ s ) in the GOMS model.In the MGOMS model, the parameter nR 2 is replaced by LAI.In this study, we assume that the reflected intensity of the shadow on the ground and on the canopy are the same (i.e., Z equals T).Thus, the model is simplified with three area-weighting components.
Table 1 presents the initial value of the parameters in the inversion process of GOMS and MGOMS models.The four structural parameters (nR 2 , b/R, h/b, ∆h/b) are obtained from the field-measured data, and the three spectral parameters (G, C, Z) are extracted from the atmosphere-corrected image SPOT-5.In B sample plots, we measured the tree crown (R), tree height (H), and clear bole height (h) of each tree.We also measured the number of trees in each sample plot to calculate n.We then obtained the datasets of nR 2 , b/R, h/b, and ∆h/b.The G, C, and Z datasets were derived from SPOT-5 data.We resampled many pixels, and then built the datasets of these three parameters.We then used the mean values of nR 2 , b/R, h/b, ∆h/b , G, C, and Z as the initial values, and the maximum and minimum numbers as the upper and lower limits.We also calculated the weighted value of each parameter (the reciprocal of the variance).The input data in the inversion process also include the multi-angle reflectance datasets constructed in Section 2.2.2 and the terrain information (slope and aspect of each pixel) extracted from the DEM.

Model Inversion Strategy
The goal of our research is to invert LAI by GOMS model directly.We thus modified the GOMS model to the MGOMS model by introducing LAI into it.The MGOMS model has many parameters.For a given observation period, the MODIS or MISR BRF products only include limited bidirectional reflectance factor data in the model inversion process.In this study, the Multi-stage, Sample-direction Dependent, Target-decisions (MSDT) inversion method [50] was adopted to segment invert the observation data and the parameters.The main objective of the iterative inversion process is to adjust the model parameters to make the model output reflectance as close as possible to the input observation reflectance data.The most sensitive observation data are used to invert the most sensitive parameters, and the previous results are used as the prior knowledge in the next inversion stage.The MSDT inversion method is based on the uncertainty and sensitive matrix (USM), which presents the sensitivity of the model parameters to the observational data in different viewing directions, and then determines the parameter inversion order in the inversion strategy.The USM can be expressed as an objective expression of the prior knowledge.Each element of the USM can be expressed as where ∆BRF(p, q) is the maximum difference of BRF calculated by the model under the circumstance that only the parameter q is changed with its uncertainty and the other parameters are fixed at their expected values; and BRF exp (p) is the BRF computed by the model at the pth geometry of illumination and viewed with all parameters at their expected values.
The inversion process in this study was performed according to a previous study [37] (see Section 4.2.1).The cost function [51] shown below was adopted where y n obs and f n (X) are the observational reflectance value and the corresponding modeled reflectance value, respectively.The variables σ n 2 and σ l 2 are the variances of the observational data and the prior distribution of parameters, respectively.The variables x l and x l prior are the parameter values and the initial values in the model, respectively.N is the number of observation samples, and L is the number of parameters.Sequential quadratic programming (SQP), which is an excellent optimization algorithm for solving nonlinear programming problems [52], was adopted to search for the cost function minimum.Such methods can be interpreted as Newton methods for the problem of finding a constrained saddle-point to the Lagrangian function.These methods solve a non-linear programming problem through the solution of a sequence of approximating QP subproblems.The optimal solution of each subproblem is used to define a search direction in which a line search is made; its solution gives the new iteration point [53].The GOMS model can be used to simulate BRF, and the simulation results can accurately reflect the hotspot and bowl-shape effect of the canopy's BRDF.The MGOMS model derived from the GOMS model must possess the same properties as the GOMS model.When calculating BRF using the MGOMS model, all the input parameters (G, C, Z, LAI, b/R, h/b, ∆h/b, θ i , ∅ i , θ v , ∅ v , n, q) are constant as in Table 1, except for LAI, and the influence of the slope and aspect are neglected.In this part, forest canopy reflectance in the principle plane is simulated with gradually increasing LAI in both the red and NIR bands by MGOMS model (Figure 5) (we regard viewing positions near the principal plane on the hotspot side as the "backward observation" and viewing positions near the principal plane opposite the hotspot as the "forward observation").The results show that BRFs decreased with increasing LAI in the red band and increased with rising LAI in the NIR band.The simulated BRFs reveal the hotspot effect of BRDF in the red and NIR bands and the bowl-shape effect in the NIR band.Additionally, the variation of BRFs with gradually increasing LAI also appropriately reflects the optical characteristics of the vegetation spectral curve.

Comparison of the Simulated BRF and MOD09GA BRF
We compared the BRF simulated by the MGOMS model with the MOD09GA BRF products at the 500-m pixel scale (the compared pixels are the pixels where the sample plots located); the solar position, view position, and DOY of the selected MOD09GA BRF data are shown in Figure 3, and the initial data for the model are shown in Table 1.
Figure 6a presents the results of the comparison of the simulated BRF and MOD09GA BRF in the NIR band with an RMSE of 0.0124 and an R 2 of 0.62. Figure 6c presents the results of the comparison of the simulated BRF and MOD09GA BRF in the red band with an RMSE of 0.0016 and an R 2 of 0.75.Figure 6b,d present the results of the comparison of the simulated BRF and MOD09GA BRF along with the view zenith angle in the NIR and red bands, respectively.The comparison results reveal that the MGOMS model can be used to simulate the BRF precisely (Figure 6a,c).The results presented in Figure 6b,d clearly indicate that the GOMS and MGOMS models have a high consistency in simulating BRF.However, the simulated BRFs were higher than the MOD09GA BRFs in the NIR band and lower than the MOD09GA BRFs in the red band.This is mainly because the simulated BRFs are closely related to the initial values, and the initial values of the model are set according to the field-measured data, which only represent the sample plot of 25 m × 25 m.However, at the 500-m spatial resolution of the MOD09GA BRF products, the unmatched pixel size scale would cause some error; for example, the field-measured LAI of the A sample plots was higher than the LAI of a 500 m × 500 m pixel (the sub-pixel coverage is higher than the pixel coverage of the A sample plot, derived from Section 4.2), which lead to overestimation of the simulated BRFs in the NIR band and underestimation in the red band.

Uncertainty and Sensitivity Analysis of the Parameters in MGOMS Model
BRFs in the NIR band can well reflect the vegetation coverage features.In our research, BRFs in the NIR band were used to invert LAI.Table 2 presents the first USM matrix (USM0) in the NIR band.The solar zenith angle, solar azimuth, view zenith angle, and view azimuth were derived from the selected MOD09GA data shown in Figure 5. On these days, the solar positions were nearly identical and, thus, we set the solar zenith angle equal to 23° and the solar azimuth equal to 143°. Figure 7 shows the first uncertainty and sensitivity analysis results in the NIR band in all directions, the solar zenith and solar azimuths are 23° and 0°, respectively, in this analysis process.The results indicate that NIRC was the most sensitive parameter, and the USM0 data of NIRC were the highest.The results also demonstrate that the uncertainty of NIRC is more sensitive in the backward observation than in the forward observation.Hence, we first invert NIRC to reduce the uncertainty of the parameter NIRC, and then calculate the average NIRC as the initial value for the second uncertainty and sensitivity analysis to determine the second most sensitive parameter, etc.Based on the results of these analyses, LAI was the third most sensitive parameter and was in the third order in the inversion process.When the sensitivities of the parameters are similar, the parameters are always inverted simultaneously.As our main purpose was to invert LAI, we did not assess whether the sensitivities of the other parameters were similar to that of LAI.Table 2.The first uncertainty and sensitivity analyses matrix (USM0) results in the NIR band (the solar zenith angle was 23°, and the solar azimuth was 143°).The multi-angles were derived from the MODIS BRF datasets, the DOY of which are shown in Figure 3. (A higher data value indicates a higher sensitivity of the parameter in the model.)According to the sensitivity analysis results, the steps for inversion and segmentation of the BRF data were as follows: (1) Invert NIRC using the NIR band reflectance datasets with backward-viewing directions and large viewing zenith angles; (2) Invert NIRG using the NIR band reflectance datasets with small viewing zenith angles; (3) Invert LAI (LAIJune) using the reflectance datasets for June with forward-viewing directions and large viewing zenith angles.

LAI Inversion Based on Retrieved Parameters by MGOMS Model
The results of the parameter uncertainty and sensitivity analysis indicated that LAI was the third most sensitive parameter in the NIR band of the MGOMS model.Thus, we first inverted NIRC, then used the inversion results as the prior knowledge to invert NIRG, and finally obtained the inversion results of LAI.Compared with the GOMS model, the MGOMS model can effectively abridge the LAI inversion process.The input parameters were described in Section 3.3.1.
Figure 8 shows the LAI inversion results for June 2008 obtained using the MGOMS model; the pixel scale is 500 m with 300 × 380 pixels.Considering the real surface condition in the Dayekou forest, the density of the trees decreases and the shape of the tree canopy becomes smaller with increasing altitude, indicating that LAI will decrease with increasing altitude.From the DEM map of Dayekou shown in Figure 1b, the altitude increases with decreasing latitude and longitude.Therefore, in Dayekou forest, the LAI should show a decrease with decreasing latitude and longitude.As shown in Figure 8, LAI decreases with decreasing latitude and longitude on the whole.The change of the retrieved LAI with decreasing latitude and longitude is approximately in agreement with the real condition in Dayekou forest.This result shows that the MGOMS model can be used to invert LAI in large areas.

Precision of the Retrieved LAI by MGOMS Model
In order to validate the LAI inversion efficiency of the MGOMS model, the retrieved LAI was compared with the field-measured LAI in situ.As described in Section 2.2.1, the forest LAI was measured using LAI-2000 and TRAC in the A (20 m × 20 m) and B (25 m × 25 m) sample plots.The retrieved LAI must be adjusted to the in situ scale to enable a robust validation.The LAI e at different scales (subpixel and pixel) has the relationship shown below [36] LAI subpixel LAI pixel = (b * nR 2 ) subpixel (b * nR 2 ) pixel (24) where LAI subpixel is the field-measured LAI and LAI pixel is the retrieved LAI.The parameters on the right side can be obtained from the LiDAR image.b can be calculated from H. H can be derived from the LiDAR data (CHM), and the statistical relationship between b and H can be built using field-measured b and H. nR 2 is considered to be approximately equal to the canopy cover factor, which can be calculated as the ratio between the number of canopy pixels and the number of corresponding CHM pixels classified with a zero threshold.Due to the sample plots A5-9 not being forest pixels in the MODIS land-cover product, there were only 24 field-measured LAI e validation data points in this part.We also compared the LAI retrieved by Ma et al. [36], i.e., LAI retrieved by the GOMS model (which cannot be directly inverted by the GOMS model, by using the structure parameter retrieved by the GOMS model and other statistical relationships to calculate LAI), with the LAI retrieved by the MGOMS model.Figure 9 shows the validation results.
Figure 9a-c show the effective LAI measured by LAI-2000 and TRAC versus LAI retrieved by the GOMS model, LAI retrieved by the MGOMS model and MODIS LAI, respectively.Figure 9d shows the results of the comparison of LAI e retrieved by the GOMS model and LAI e retrieved by the MGOMS model.Both the GOMS and MGOMS models can achieve high accuracy in calculating LAI and were in good agreement with the field-measured LAI, with R 2 of 0.71 and 0.64, respectively.However, MODIS LAI produced much lower R 2 values compared with the field-measured LAI (Figure 9d).The RMSE value between the field-measured LAI and the MODIS LAI (3.25) is much higher than that between the field-measured LAI and the GOMS and MGOMS LAI (0.71 and 0.67).Figure 9d

Discussions
The MGOMS model, which is derived from the GOMS model, has the same characteristics in modeling BRF as the GOMS model.In Figure 6a,c, the results suggest that the BRF simulated by MGOMS model is closely related to the MODIS BRF data.In the NIR band, the R 2 and RMSE are 0.62 and 0.0124, respectively, and in the red band, the R 2 and RMSE are 0.75 and 0.0016.The relatively high R 2 and low RMSE indicate that the MGOMS model can be used to simulate BRF.In Figure 6b,d, it is clear that the BRF datasets simulated by the GOMS and MGOMS models are nearly equivalent when the field-measured data are used as the initial value.This result further demonstrates that MGOMS model is consistent with GOMS model.
The GOMS model is the most widely-used physical model to invert LAI because the model builds the canopy structure parameters and BRF together.With the multi-angle BRF datasets, such as the MODIS and MISR datasets, it is easy to build the invert function and retrieve all parameters belonging to the GOMS model.The MGOMS model is derived from the GOMS model and, thus, shares this advantage.Previous studies [36,37] of LAI inversion using the GOMS model have inverted the canopy structure parameters (e.g., b/R, h/b, nR 2 ) first and then used the statistical model to calculate LAI.By contrast, the MGOMS model can retrieve LAI directly.The USM results indicated that LAI is in the third order of the MSDT inversion process, and using the MGOMS model abridges the LAI inversion process.
Compared with the field-measured LAI, the LAI e values retrieved by both the MGOMS and GOMS models are in good agreement with the field-measured LAI, with R 2 of 0.64 and 0.71, respectively.Although the R 2 of LAI e retrieved by the MGOMS model is slightly lower than that of LAI e retrieved by the GOMS model, the RMSE of LAI e retrieved by MGOMS (0.67) is lower than that of LAI e retrieved by GOMS (0.72).These results all demonstrate that the MGOMS model can be used to invert LAI.
However, the GOMS model was constructed based on the Li-Strahler geometric-optical model, which considers the mutual shadowing of crowns and makes the geometric optic model more suitable for highly-dense canopy forests.However, in our modified process, when considering the mutual shadowing of crowns, the P(s|n) function is difficult to build, and the P gap function will be very complex.In order to make the function between LAI and nR 2 laconic, we assumed that the crowns are non-overlapping and that the crowns in the field of view are uniformly distributed.Based on this assumption, the MGOMS model may be less precise than the GOMS model in dense forest, with a smaller applicable scope than the GOMS model.Thus, further research must address the development of the LAI and nR 2 relationship function when the crowns are non-overlapping and the crowns are Poisson-distributed.

Conclusions
In this study, the GOMS model was modified by using LAI to replace nR 2 .The modified GOMS (MGOMS) model enabled the direct inversion of LAI, and LAI is a relatively highly-sensitive parameter in the MGOMS model.By using the MGOMS model, the LAI inversion process was effectively reduced.When comparing the retrieved LAI e with the field-measured LAI e , though the LAI retrieved by the MGOMS model had a slightly lower R 2 (0.64) than the LAI retrieved by the GOMS model (0.71), the LAI retrieved by the MGOMS model was closer to the field-measured LAI (with an RMSE of 0.67) than the LAI retrieved by the GOMS model (with an RMSE of 0.72).All of these results indicate that the MGOMS model can be used to invert LAI effectively and precisely.However, many problems remain to be addressed in the future, such as, in the modeling process, we make the assumption that the canopies are non-overlapping and uniformly distributed.We also assume that the view zenith angle is zero.Moreover, we set the parameter G to a constant value of 0.5.All of these assumptions can affect the accuracy of the model.In addition, the precision of the prior information also influences the calculated results of the MGOMS model.Therefore, in a forthcoming study, we will focus on how to expand the applicability of the MGOMS model and how to solve the uncertainty of the input parameters.

Figure 1 .
Figure 1.(a) Dayekou land-cover type from the MODIS land-cover product.The spatial resolution is 500 m.Pixels in green, pale yellow, pink, and red are forest pixels, grassland pixels, cropland pixels, and urban areas, respectively.Pixels in brown yellow represent land for other use; (b) Digital Elevation Map (DEM) map of Dayekou.The size of the map is the same as that in Figure 1a.The spatial resolution is 90 m.In this map, altitude generally increased with decreasing latitude and longitude.

Figure 2 .
Figure 2. Sample plot locations in our study area in the MODIS land-cover products.Pixels where the A sample plots were located are outlined in red, and the pixels where the B sample plots were located are outlined in black.B sample plots are the part of a super sample plot.

Figure 3 .
Figure 3. Solar position and view position of the selected MOD09GA datasets in uncertainty and sensitivity analysis processes.The datasets span a MODIS 16-day observation period.The numbers 154-177 are the DOY of the MODIS BRF products in 2008.

Figure 4 .
Figure 4. (a) Forest canopy shape in ellipsoid.R is the horizontal radius of an ellipsoidal crown, b is the vertical half axis of an ellipsoidal crown, and h is the height at which a crown center is located; (b) Simplified schematic diagram of forest vegetation canopy in the GOMS model[31].θ i and θ v are the revised solar zenith angle and view zenith angle, respectively.∅ i and ∅ v are the solar azimuth and view azimuth, respectively.∅ i − ∅ v is the azimuthal difference between the illumination and viewing directions.τ i and τ v are the sunlit shadow and viewed shadow, respectively.The shaded area is the mutual shadowing area of the sunlit shadow and viewed shadow.

4. 1 . 1 .
Tendency of the BRF Simulated by MGOMS Model along with the Variation of LAI

Figure 5 .
Figure 5. BRF estimated by MGOMS along with the increased LAI against the view zenith angle in the principal plane with a solar zenith angle of 45° in the (a) NIR band and (b) red band.LAI was gradually increased from 2 to 5 with a data change interval of 0.5.Other parameters in MGOMS were fixed.Negative view zenith angle means the forward observation.

Figure 6 .
Figure 6.BRF comparison results for simulated BRFs and MOD09GA BRFs.Comparison of the simulated BRF by the MGOMS model and MOD09GA BRF in the (a) NIR band and (c) red band; Comparison of the simulated BRF (simulated by the GOMS and MGOMS models) and MOD09GA BRF along with the view zenith angle in the (b) NIR band and (d) red band.(Negative view zenith angles means the foreward observation).

Figure 7 .
Figure 7. Polar coordinate graph of the first uncertainty and sensitivity analyses matrix results in the NIR band in all view directions.The solar zenith and solar azimuth were 23° and 0°, respectively.The data values have no unit.A higher data value indicates a higher sensitivity of the parameter in the model.

Figure 8 .
Figure 8. LAI inversion results of Dayekou forest in June 2008.The retrieved LAI data value is 0 to 6.The size of the map is identical to that in Figure 1a.The pixels in the inversion process are only the forest pixels, and the others are set to zero and colored white in this map.
compares the LAI e retrieved by the MGOMS model and the LAI e retrieved by the GOMS model.The high R 2 (0.96) and low RMSE (0.5243) indicate that the MGOMS model was very consistent with the GOMS model, combined with other statistical relationships, for calculating LAI.The LAI e retrieved by the MGOMS model was higher than the LAI e calculated by the GOMS model, but the LAI e retrieved by the MGOMS model was closer to the field-measured LAI e , with a lower RMSE (0.67).

Figure 9 .
Figure 9.Comparison between the retrieved LAI and field-measured LAI e of the sample plots: (a) LAI e retrieved by the GOMS model compared with field-measured LAI e , with R 2 of 0.71 and RMSE of 0.72; (b) LAI e retrieved by the MGOMS model compared with field-measured LAI e , with R 2 of 0.64 and RMSE of 0.67; (c) MODIS LAI compared with field-measured LAI e , with R 2 of 0.26 and RMSE of 3.25; (d) LAI e retrieved by the MGOMS model compared with LAI e retrieved by the GOMS model, with R 2 of 0.96 and RMSE of 0.52.

Table 1 .
Initial values of the input parameters *.
* RG, RC, and RZ are G, C, and Z in the red band, respectively, and NIRG, NIRC, and NIRZ are G, C, and Z in the near-infrared (NIR) band, respectively.The initial values are the mean values of the parameters; the lower and upper limit values are the thresholds of the parameters.