Retrieving leaf area index from multi-angular airborne data

This work is aimed to demonstrate the feasibility of a methodology for retrieving bio-geophysical variables whilst at the same time fully accounting for additional information on directional anisotropy. A model-based approach has been developed to deconvolve the angular reflectance into single landcovers reflectances, attempting to solve the inconsistencies of 1D models and linear mixture approaches. The model combines the geometric optics of large scale canopy structure with principles of radiative transfer for volume scattering within individual crowns. The reliability of the model approach to retrieve LAI has been demonstrated using data from DAISEX99 campaign at Barrax, Spain. Airborne data include POLDER and HyMap data in which various field plots were observed under varying viewing/illumination angles. Nearly simultaneously, a comprehensive field data set was acquired on specific crop plots. The inversions provided accurate LAI values, revealing the model potential to combine spectral and directional information to increase the likely accuracy of the retrievals. In addition, the sensitivity of retrievals with the angular and spectral subset of observations was analysed, showing a high consistency between results. This study has contributed to assess the uncertainties with products derived from satellite data like SEVIRI/MSG. Mailing address: Dr. Francisco Javier García-Haro, Department of Thermodynamics, University of Valencia, C/Dr. Moliner 50, 46100 Burjassot, Valencia, Spain; e-mail: J.Garcia.Haro@uv.es


Introduction
The LSA SAF Project is part of the ground segment for the EUMETSAT missions METEO-SAT Second Generation (MSG) and European Polar System (EPS), developed by the ESA.Our aim is to develop robust and operational algorithms for retrieving vegetation parameters from the synergistic use of SEVIRI/MSG and AVHRR-3/EPS.These instruments will offer innovative angular capabilities for determining vegetation products over Europe and Africa thanks to concomitant multiple viewing and illumination geometries ( Van-Leeuwen and Roujean, 2002).Remotely sensed BRF data are the only way to monitor the vegetation on a global scale.Traditional approaches rely on the exploitation of the reflectance variations in the spatial, temporal and spectral domains, in lack of directional information.However, observations demonstrate that the anisotropic behaviour of the surface is an important property and the absence of angular information inevitably produces a bias in the parameter evaluation (Roujean et al., 1992;Leblanc et al., 1997).Our main concern is to retrieve fractional Leaf Area Index (LAI) corrected from the effects of the sun-target-sensor geometry.
Linear mixture is the basis hypothesis in some simple BRDF (Bidirectional Reflectance Distribution Function) models (e.g., geometricaloptical or kernel-driven models) for simulation of scene reflectance or in remote sensing algorithms, e.g., Spectral Mixture Analysis (SMA).However, this technique is ineffective to address the presence of multiple scattering, surface anisotropy and the scaling processes (Qin and Gerstl, 2000).Most BRDF one-dimensional models imply surface homogeneity within an image pixel and are, therefore, limited to address mixed landcovers pixels, which are common in coarse resolution satellites.Different hybrid models have been developed in recent years to describe the radiation regime in forest canopies (Li and Strahler, 1992;Chen and Leblanc, 1997).These models assume a medium consisting of gaps and regions idealized by a turbid volume with a foliage density of small leaves.In the geometric-optical approach, an overlap calculated from the shape of the crowns allows the estimation of the proportion of shadows cast as a function of view direction relative to the hot spot direction.More recently, the GHOST model (Lacaze and Roujean, 2001b) was developed to address the local scale angular structure of the hot spot, which can be judged relevant for patch and regional scales.Although it was specifically designed to simulate the BRDF of boreal forests, it is also suitable to study simpler vegetation canopies like crops and pasture (Lacaze et al., 2002).This study proposes an operational approach, namely DISMA (DIrectional SMA), to deconvolve the angular reflectance into single landcovers reflectances, attempting to solve the inconsistencies of 1D models and linear mixture approaches.The model formalism and the volume scattering formulation are similar to the GHOST model.However, the between-crown gap probability is formulated in terms of the fractional vegetation cover ( fCover) and a geometric variable (h) associated with the shape of plants.A random spatial distribution of plants is assumed to simplify its computation.The invertibility of the model to retrieve LAI was demonstrated using airborne POLDER and HyMap multi-angular measurements corresponding to cropland.The next section describes the model formulation.The inversion algorithm is presented in Section 3. The model is validated in Section 4. Finally, the conclusions and future prospects are presented in Section 5.

Model formulation
The scene consists of subcanopies idealized by geometric elements superposed on a flat background.The reflectance of an individual pixel is assumed to consist of an area weighted linear combination of the soil and vegetation contributions.The vegetation reflectance is expressed as the sum of single scattering (ss) and multiple scattering (ms) reflectances (Hapke, 1981;Lacaze and Roujean, 2001b) (2.1) where leaf albedo ~is the sum of the leaf reflectance (t) and transmittance (x), n s, v= cosis,v, p is the phase angle (i.e. the angle between the sun and the viewer), and the factors I ss and Ims model the proportion of radiation flux which is single/multiple scattered by foliage elements on the downgoing and outgoing optical pathways as a whole, H(n s, v) is the Chandrasekhar function for multiple scattering (Hapke, 1981) and P(p) is a turbid medium phase function (Ross, 1981;Lacaze and Roujean, 2001b).The GHOST model relies on the coupling between a simple hot spot formula (Roujean, 2000) and the Gfunction that describes the canopy geometry.In our model, G controls only the volume component depending on the within-crown element distribution, whereas the external geometric component depending on the crown shape and dimension is evaluated using an average theory of the gap probability.We assume that the geometrical component and the volume component of the radiation fluxes can be decoupled (2.2) The volume single/multiple scattered component can be described using Beers' Law (2.3) where ΩE is the clumping index of the shoots, which quantifies the level of foliage aggregation within the tree crown and is generally dependent on the view angle (Kucharik et al. 1997); ΩE=1 means random foliage distribution and ΩE<1 means clumped foliage; LAI is the total canopy overstorey leaf area index; gc denotes the fCover and c is a band-specific factor, which was assumed to be dependent on the leaf transmittance (e.g., Bégué, 1993).Assuming plants with similar foliar density in the scene, LAI/gc represents a mean value of the LAI of individual plants.This term is more pertinent to describe crown trees transparency than total LAI.∆ denotes the (bidirectional) normalised extinction coefficient for singly scattered radiance.The model adopts the analytical expression derived by Roujean (2000) to the hot spot effect, i.e. coupling the downgoing and outgoing optical pathways where G is the well-known Ross function (Ross, 1981).For the multiple scattered component the hot spot phenomenon is ignored (Qin and Goel, 1995).
The geometric component of the single scattering I ss, geo is determined by the between-crown light penetration and the visibility of illuminated objects.This component is particularly relevant for discontinuous canopies.The downgoing and outgoing terms of the intercepted flux were formulated in terms of horizontal projection of the crown.Let P0(i) denote the (between-crown) monodirectional gap fraction, which corresponds to the fraction of soil seen in the direction i (Nilson, 1971).For homogeneous Poisson distribution, the probability of observing the ground under the tree crowns in any given pixel approaches to P0(is)=(1−gc) (hs+1)  (Jasinski and Eagleson, 1989), where h is defined as the ratio of ground projected shadow to plant area.It absorbs all the geometric factors that relate canopy area to shadowing area into only one variable.Its analytical expression for the most common geometrical bodies is provided in Jasinski and Eagleson (1990).The probability of having sunlit ground/ crown (P ig, c) and viewed ground/crown (Pvg, c) can be expressed as follows: A functional relationship can also be found between subpixel shaded ground Psg and fractional vegetation cover This equation is applicable at large sampling scales, when imaging stands resolutions greater than the size of the tree crowns, as confirmed in a preliminary analysis using simulated data (García-Haro et al., 2002).The probability of observing sunlit crown when P ic and Pvc are not correlated is simply the product of both probabilities PicPvc, where Pic controls the amount of light intercepted by crowns and Pvc controls the contribution of visible crowns to the scene radiance.However, hot spot kernels are necessary to account for the correlation between the two gap probabilities along sun and view directions (Chen and Leblanc, 1997;Qin and Goel, 1995).We assume that the hot spot has a minor influence on the multiple scattered interception at a crown level, i.e.I ms, geo= PicPvc, but introduce a hot spot kernel Fc to modulate the dependence between the optical paths for the single scattered radiation (2.7) where Fc is usually obtained from the overlap function between viewing and illuminated shadows as projected on the background (Li and Strahler, 1992;Schaff et al., 1994).Finally, the soil contribution Rs is expressed as the product of the background bidirectional reflectance cs and the vegetation transmittance T. The vegetation transmittance is the sum of the probability that a solar ray beam will reach the ground without intercepting any crown Tgeo plus the probability of intercepting a crown without hitting any foliage element (1-Tgeo)Tvol, i.e.T = Tgeo+(1 -Tgeo)Tvol.The following expressions were considered for the transmissions: where F G is a hot spot kernel.It tends to zero as the sun and view directions are far apart (p>>0 o ), i.e. when the viewer sees the sunlit ground through a gap different from that of illumination.The hot spot gradient function also decreases as the view zenith angle approaches a horizontal view perspective (iv→90 o ) and as the leaf density decreases.In this work we have used hot spot kernels FG and Fc similar to the hot spot function proposed by White et al. (2001), which match the hot spot region measured in many BOREAS sites (Leblanc et al., 1997).
Our model neglects side scattered diffuse radiation incident on ground surface after multiple scattering with leaves of neighbouring plants, which can be significant over a bright background.The model is thus less accurate at near-infrared wavelengths at which multiple scattering in plant canopies is the strongest within the solar spectrum.Considering also the presence of skylight irradiance, the reflectance can be expressed as a sum of the direct and isotropic diffuse components, where the component of diffuse reflectance is calculated as the average of the measurements in the principal plane.In addition, the model allows us to derive other parameters like the entire BRDF distribution, albedo or absorptance, including the relative contribution of vegetation and soil.

Model inversion
The inversions are achieved simultaneously for all spectral bands, i.e. by coupling the spectral and directional data available.The variables to be retrieved are LAI and fCover.Without a priori information on these variables, the inverse problem typically consists in determining the optimal set of variables that minimizes the distance between observations and modelled values.Given BRDF values r i (i=1,..., N), representing the conditions of observation (i.e.wavelengths and view and illumination geometries), the retrievals are performed by comparing observed and modelled y i (i = 1,.).We must note, however, that prior information about the most probable solution values, their uncertainty and their correlation may significantly improve the accuracy of the retrievals, especially for ill posed inversion problems (Combal et al., 2002).In order to reduce the computations, we used an automatic method to select the solution bounds.Firstly, spectral mixture analysis was applied to obtain fCover from anisotropically corrected reflectance using as input endmembers representing soil and vegetation.Secondly, an empirical relationship was used to derive LAI from estimated fCover (Lacaze and Roujean, 2001a).The estimated values for LAI and fCover were finally used to define the bidimensional bounds in the solutions domain.
The model requires the optical properties of the underlying soil as c s an input.We employed the RPV model (Rahman et al., 1993) to extend the BRF measured at a few angles to the entire BRDF distribution.A variable configuration using two different soil BRDFs was used to represent more realistically the influence on reflec- tance of soil moisture and roughness (see fig. 1).
Although leaf reflectance and transmittance are not equal in many plants and vary greatly between species, it was assumed for simplicity that they are similar.We have proposed a method to estimate the leaf albedo ~from the data themselves.The method relies on the assumption that pixels presenting a negligible contribution of the soil background (i.e.f Cover =1 and high LAI values) can be found in the scene.In these pixels eq. ( 2.1) is inverted in order to retrieve ~using an iterative method.Results have indicated the adequate convergence of the method, providing reasonable results irrespectively of the data set considered (fig.1).

Validation
The model invertibility was evaluated using data acquired during the 1999 Digital Airborne Spectrometer EXperiment (DAISEX) campaign by the sensors HyMap onboard the DLR Do-228 aircraft and POLDER onboard an ARAT plane.The experiment site selected by ESA for the DAISEX campaigns is a 4 km by 4 km area centred at 39º3lN, 2º5lW, which is located 28 km from Albacete (Spain), (Camachode Coca et al., 2002).The POLDER instrument allows a measurement of surface reflectance directional effects at nine spectral bands (centred at 443, 500, 550, 590, 670, 700, 720, 800 and 864 nm wavelength) in the visible and near infrared.The CCD matrix permits collection of bidimensional images in one shot.The along track and cross track FoV is of ±43º and ±51º respectively.Four flights were undertaken during 3-5 June 1999 at a typical airborne altitude of 3000 m, with a spatial resolution of 20 m.Each flight recorded around 140 spectral images.They provide around 50 values of angular reflectance for every pixel, irregularly distrib- uted in the viewing hemisphere, and different for every pixel.The images were calibrated, geo-coded and corrected for atmospheric effects as it is specified in Leroy et al. (2001).In this work the BRDF was interpolated for the full range of view angles for every pixel.The BRDF was then retrieved considering uniform sites of 3×3 pixels (60 ×60 m 2 ).The model was inverted against a set of bidirectional reflectance factors taken along the principal plane (i.e. the region with the strongest anisotropy) and along the orthogonal plane.Figures 2a,b and 3a,b show a few examples of measured BRF's along with the values predicted by DIS-MA.The results demonstrate the potential of the model to accurately explain the spectral and angular variations of the data.DISMA captures the essential BRDF features such as bowl shape, backscattering increase in reflectance and broad hot spot, as well as the spectral contrast between soil and vegetation.The model works better in the visible regions, when the first order scattering effects predominate but in the NIR region (band 8) the discrepancies are higher.In this region directional effects are less apparent due to the reduction of contrast between canopy components and the prevalence of multiple scattering.Another source of error is attributable to unquantified influences of foliage clumping.
The HyMap instrument (http://www.hyvista.com/)has 128 spectral bands and a high   Vegetation in situ measurements corresponding to major agricultural units were taken during the campaign.The measured properties included LAI, fCover, canopy height, biomass and chlorophyll content.The LICOR LAI-2000 instrument was used to measure LAI.Sampling took place at several times along parallel transects inside the fields (García et al., 2001).The comparison between in situ LAI measurements and the retrieved values using different spectroangular datasets is shown in fig.6.A linear relationship was found in all cases, irrespective of the data set considered.Table I shows the coefficients of the linear fit, revealing a significantly high correlation, with r 2 values higher than 0.92 in all cases.A chi-square (| 2 ) test was performed assuming a LAI uncertainty of 0.4.It produced values below the chi-square critical value of 92 in all cases (for 71 degrees of freedom and a probability of 0.05), confirming thus the predicted linear relationship.Moreover, re-sults indicate a good correspondence between field measurements and retrievals, since the linear relationship is close to the 1:1 line, and the RMSE is relatively low (0.35-0.48).Another   important result is that DISMA is able to produce reasonable results under sub-optimal conditions either in the angular sampling (i.e.taking the orthogonal plane) or with a small number of spectral channels, e.g., SEVIRI wavebands.
We must also note that even though the retrieved LAI reproduces the between crops LAI variability observed in the field measurements very well, it fails to address the within crop variability (e.g., LAI retrievals tend to saturate within each individual crop).One reason for this is that pixel-by-pixel comparison is strongly hampered by the inaccuracies in the georegistration of field measurements and imagery.Another reason is the averaging process performed to the images.One important aspect is the consistency of the retrievals, i.e. its sensitvity with respect to the set of spectro-angular measurements taken in the inversion.The crosscomparison between results from different data sets showed a high consistency between them, with RMSE values typically lower than 0.30 (see table II).Finally, we must indicate that although in situ measurements included only a limited number of fCover sampling, the retrieved values were highly coherent with the available field information.

Conclusions and prospects
This study aimed to develop an operational approach to deconvolve the angular reflectance into single landcovers reflectances, attempting to solve the inconsistencies of 1D models and linear mixture approaches.The model relates the spectral and angular variation with the main optical and structural parameters of discontinuous canopies, like LAI.The inversion revealed the model potential to combine spectral and directional information to increase the likely accuracy of the retrievals.Results also indicate the effectiveness of the algorithm using only SEVIRI channels and information of the diurnal sampling as an angular signature of the surfaces.We preferred to simplify the selection of secondary parameters and the inversion algorithm with specific intent that it not to be canopy dependent.However, for the application of DISMA to complex scenarios as in global studies the stratification of scene is convenient in order to optimise the model inputs to the knowledge of ecosystem characteristics, reducing misidentification and saving computations (García-Haro et al., 2003).One important aspect is the sensitivity of retrieved variables to the information used to parameterise vegetation canopy radiative transfer.Although the inversion algorithm was satisfactory, difficulties still arise.For example the estimation of leaf albedo is known to be impaired by canopy-level variables like LAI.Future research is also needed to test the model on BRDF datasets comprising different vegetation types (shrubland, forest). ,

Fig. 1 .
Fig. 1.Spectral inputs introduced to DISMA, which correspond to leaf albedo of green (solid) and dry vegetation (dotted), and reflectance of two different soils.

Fig. 4a ,
Fig. 4a,b.Measured reflectance (dashed lines+symbols) and simulated (solid lines) at five different spectral channels over sugar beet (a) and alfalfa (b).Configurations 1-3 correspond to the N-S flight and 4-6 to the E-W flight.Within each flight line, consecutive numbers correspond to 8:00, 12:00 and 15:00 UTC, respectively.The figures show several examples of TM-like wavebands derived from real HyMap data.

Fig. 6 .
Fig. 6.Relationship between field-measured LAI and airborne-derived LAI using different data sets.The figure at the top corresponds to POLDER data in the principal (left) and orthogonal (right) plane.Different symbols are used to represent corn ( ), sugar beet ( ) and alfalfa ( * ).

Table I .
Statistics of the correlation between in situ measurements and retrievals of LAI.

Table II .
Statistics of the cross-comparison of LAI as retrieved using different data sets.