Retrieval of Biophysical Crop Variables from Multi-Angular Canopy Spectroscopy

Abstract

The future German Environmental Mapping and Analysis Program (EnMAP) mission, due to launch in late 2019, will deliver high resolution hyperspectral data from space and will thus contribute to a better monitoring of the dynamic surface of the earth. Exploiting the satellite’s ±30° across-track pointing capabilities will allow for the collection of hyperspectral time-series of homogeneous quality. Various studies have shown the possibility to retrieve geo-biophysical plant variables, like leaf area index (LAI) or leaf chlorophyll content (LCC), from narrowband observations with fixed viewing geometry by inversion of radiative transfer models (RTM). In this study we assess the capability of the well-known PROSPECT 5B + 4SAIL (Scattering by Arbitrarily Inclined Leaves) RTM to estimate these variables from off-nadir observations obtained during a field campaign with respect to EnMAP-like sun–target–sensor-geometries. A novel approach for multiple inquiries of a large look-up-table (LUT) in hierarchical steps is introduced that accounts for the varying instances of all variables of interest. Results show that anisotropic effects are strongest for early growth stages of the winter wheat canopy which influences also the retrieval of the variables. RTM inversions from off-nadir spectra lead to a decreased accuracy for the retrieval of LAI with a relative root mean squared error (rRMSE) of 18% at nadir vs. 25% (backscatter) and 24% (forward scatter) at off-nadir. For LCC estimations, however, off-nadir observations yield improvements, i.e., rRMSE (nadir) = 24% vs. rRMSE (forward scatter) = 20%. It follows that for a variable retrieval through RTM inversion, the final user will benefit from EnMAP time-series for biophysical studies regardless of the acquisition angle and will thus be able to exploit the maximum revisit capability of the mission.

1. Introduction

The retrieval of biophysical plant variables from optical imagery has been playing an important role in remote sensing and ecosystem modelling for more than 30 years. With ongoing technical progress of the sensors, there is also a steady demand for improved extraction of information from the gathered data. Especially in the agricultural context, many studies have pointed out the suitability of multispectral data (e.g., [1,2,3,4,5]), hyperspectral data (e.g., [6,7,8,9,10,11]) and a combination of both (e.g., [12,13,14]) for an assessment of crop characteristics. In order to make these benefits available to modern farming, scientific tools and algorithms need to be directly applicable for a broader user community. Variables like the leaf area index (LAI) or leaf chlorophyll content (LCC) are of prime importance for a proper characterization of the canopy and plant biochemistry [15].

Several approaches are known to successfully retrieve hyperspectral canopy variables from measured spectra. The approaches can be divided into empirical and generic methods. The former build up a statistical relationship between vegetation spectral signatures and in situ measured variables as parametric or non-parametric regressions [4]. Due to this site-specific linkage, empirical methods are not transferable in space or time [16,17]. To become independent of in situ measurements, more generic approaches often make use of radiative transfer models (RTMs). They are the intermediate link between biophysical characteristics of the canopy and its geometry, radiometric interaction and the reflected radiation [18]. Location, intensity and quality of the radiation source, atmosphere, vegetation/canopy, soil as well as position and properties of the sensor are important subsystems for the remote sensing of vegetation [19]. RTMs separate exterior parameters from the influence of the target itself, allowing quantitative analysis and the establishment of distinct relationships between signal and object variables [18]. One of the major improvements in RTMs was the incorporation of arbitrarily inclined leaves instead of a representation by plates. The resulting SAIL model [18] (Scattering by Arbitrarily Inclined Leaves) was later coupled with the leaf optical properties model PROSPECT [20] to form the new fusion model PROSAIL [21].

In the direct or forward mode, PROSAIL simulates synthetic spectra from input variables that describe plant physiology and canopy architecture. In the indirect or inverse mode, these variables are obtained from spectral signatures. Inversion techniques are either based on optimization methods, artificial neural networks (ANN), machine learning algorithms (MLA) or look-up-tables (LUT) (see [22] for overview). Their advantages and drawbacks vary with purpose of use. Optimization methods aim at minimizing deviations between modelled and measured spectra [23]. Such minimization algorithms continuously change the input variables of the RTM until the modelled result matches the observation as closely as possible, leading to comparatively long computation times [24]. ANNs and MLAs on the other hand are quicker in training and execution, but they require a priori information, calibration and lack of mathematical transparency [12,25]. Look-up-tables are databases of modelled spectra and their associated input parameter configurations. LUTs are known to be fast and robust methods producing reasonable results (e.g., [3,6,26,27,28]). In a first step, the LUT is built up in forward mode before it can be browsed in inverse mode. For the compilation of the LUT, the user has the choice of size (number of simulations), artificial noise type and noise level of the spectral model output as well as distribution type and constraints for all input parameters. Inversion of RTMs is impeded by the fact that more than one combination of variables can lead to the same model result. This effect has become known as equifinality or ill-posed problem and is dealt with either by restriction of the input range or by inclusion of the n-best performing results rather than just considering the number one fit [24].

For an ideal analysis, communication between sensor and model must be optimal. Since spectral models have been developed in the laboratory with the help of ground-based spectrometers, they basically are of hyperspectral nature. In order to use the models in combination with multispectral data, their spectral resolution normally is toned down using the spectral response functions of the respective instruments. Using hyperspectral data as input allows for making full use of the quasi-contiguous narrowband output of the RTM in forward mode. The retrieval methods tested in this article thus are intended to be applied on hyperspectral data of the Environmental Mapping and Analysis Program (EnMAP). EnMAP is a German spaceborne imaging spectroscopy scientific mission carrying the EnMAP Hyperspectral Imager (HSI) instrument [29]. Currently under development, EnMAP-HSI will deliver data at high spectral resolution of 6.5 nm in the VNIR and 10 nm in the SWIR domain which together cover the full spectral range of 420 to 2450 nm [30]. Competition for actual data is expected to be intense, since the data take capacity of EnMAP is limited and—as of today—only the Italian hyperspectral mission PRISMA [31] may be going to record comparable data by the estimated time of launch in 2019. Repeat cycles of 23 days in quasi-nadir mode will limit the availability of cloud-free scenes [32]. To mitigate this problem, the satellite platform will be capable of a max. ±30° across-track tilt, allowing side looks upon the target with revisit times of up to 4 days [33] near earth’s equator or even less for latitudes of central Europe, e.g., 2.5 days for Munich, Germany [32]. The effects of this off-nadir pointing for the retrieval of biophysical variables have not been tested in the EnMAP context. Therefore, the objectives of this study are (1) to demonstrate the expected impact of the EnMAP-specific sun–target–sensor-geometry (s–t–s-geometry) on reflectance spectra, (2) to quantify the effect on agriculturally relevant variable retrievals, such as LAI and leaf pigments and (3) to introduce a new hierarchical LUT approach for an optimized retrieval of these parameters.

2. Materials and Methods

2.1. Study Area & Sampling Layout

The study area is located in the North of Munich, Bavaria, in Southern Germany. Two study sites at 48°17′31.25″N, 11°42′21.53″E (field 517) and 48°14′51.46″N, 11°42′24.10″E (field 509) were visited regularly during two field campaigns. Both fields are part of communal farmland belonging to the city of Munich. Each was cultivated with winter wheat (triticum aestivum), representing the dominant cereal crop in the area and situated within 1.5 km distance to the Isar river. The average cloud cover in the Munich-North Isar (MNI) region was a bit higher than usual (5.68 instead of 5.44 okta). This indicates the difficulty of recording spectra on a frequent basis, since adequate clear sky conditions occurred only on few occasions. Information on the site management was provided via personal communication by the farm managers; i.e., dates of seeding, fertilization methods and quantities, harvesting dates, etc.

During the first campaign, data were collected at almost weekly intervals from 17 April to 25 July 2014 (14 sampling dates). The second campaign period already started in autumn. Measurements were conducted from 28 November 2014, to 21 July 2015 (13 sampling dates). In this way, the complete growing cycle of the crop from seeding to harvest could be observed. One elementary sampling unit (ESU) was defined as a 10 m × 10 m pixel size. The measurements were then related to a 3 × 3 ESU raster with equal distances of 10 m. All nine ESUs were marked with sticks and revisited for each sampling date. The row azimuth direction of the winter wheat crops was 170°/350° for 2014 and 150°/330° for the 2014/15 season, with the angular definition of 0° = N.

Table 1. Dates of field visits and corresponding availability of spectral observations nadir/angles as well as crop variable measurements, indicated with check marks. Crop variables were measured weekly for both growing cycles of 2014 (left) and 2014/2015 (right). The frequency of spectral observations was subjected to weather conditions, since reflectance measurements require a cloud-free sky.

Date	Spectral Nadir	Spectral Angles	Crop Variables	Date	Spectral Nadir	Spectral Angles	Crop Variables
17 April 2014	√	√	√	28 November 2014			√
23 April 2014	√	√	√	12 December 2014	√	√	√
30 April 2014			√	19 March 2015	√		√
6 May 2014			√	10 April 2015	√	√	√
14 May 2014			√	22 April 2015	√	√	√
9 May 2014	√	√	√	5 May 2015			√
26 May 2014			√	8 May 2015	√	√	√
2 June 2014	√	√	√	3 June 2015	√	√	√
6 June 2014	√	√	√	12 June 2015	√		√
18 June 2014	√	√	√	1 July 2015	√	√	√
26 June 2014	√	√	√	10 July 2015	√	√	√
3 July 2014	√		√	16 July 2015	√		√
17 July 2014	√		√	21 July 2015	√	√	√
25 July 2014	√	√	√
Total observation number	10	8	14	Total observation number	11	8	1

shows the complete list of sampling dates.

2.2. In Situ Measurements

2.2.1. Spectral Data

Spectral data were collected with an Analytical Spectral Devices Inc. (ASD) (Boulder, CO, USA) FieldSpec 3 Jr. Five separate measurements were carried out per ESU and per observation angle to obtain representative values. Outliers were removed and the remaining spectra were averaged and subjected to further processing. The post-processing included splice correction, radiometric calibration to absolute reflectance values and smoothing with a moving Savitzky-Golay-filter [34]. Apart from nadir measurements, the canopy was also measured under observer zenith angles (OZA) of +30° and −30° regarding the solar plane: a sensor inclination towards the sun is defined as a positive OZA, whereas an inclination away from the sun is described as a negative OZA (Figure 1). Due to the backscatter effect, spectra with positive OZA are noticeably brighter than nadir views or negative zenith angles, as they draw nearer to the spot of increased backscatter, also known as the hot spot [35]. The opposite direction shall accordingly be called cold spot or forward scatter and usually leads to reduced reflectances and darker images. For the angular spectral measurements, a microphone stand was modified to hold the ASD glass fiber optic. The horizontal rod of the stand could be raised or lowered to adjust the viewing angle with help of an attached inclinometer. The observer azimuth angles (OAA) matched up with the row azimuth angle of the canopy stands (170° for field 517 and 150° for field 509). EnMAP will operate on a sun-synchronous orbit with 97.96° satellite inclination angle descending node [36] which corresponds to an OAA of 187.96°. The angular effects measured in the presented campaign therefore are assumed to adequately represent the angular effects expected from future EnMAP data.

Spectral information from each of the 3 × 3 ESUs was compiled to pseudo-images with a ground sampling distance of 10 m. Furthermore, the processed signatures were converted into simulated EnMAP spectra via the EnMAP-end-to-end-Simulator (EeteS) [37]. In this process, the sensor-specific radiometric and spectral properties were adapted. The spatial resolution in this case was retained at 10 m to preserve the data population.

The gap fraction is a measure for the probability of a ray of light to penetrate through the canopy undisturbedly [38]. Accordingly, this parameter decreases with density and/or height of a canopy. Canopy height can also be seen as a path length on which energy can interact with plant traits. Assuming identical canopy height, the path length is shortest for nadir views and increases with OZA > 0°. For a 30° deflection from nadir, the travelled path is longer by factor cos (30°)⁻¹ which is 15.5%. This leads to a weaker influence of soil background and to an apparently higher portion of visible leaf surface.

The anisotropy factor (ANIF) [39] yields useful information about the sensitivity of different wavebands towards changes in illumination geometry. It is simply defined as Equation (1):

$$\mathrm{ANIF}=\frac{{\mathrm{R}}_{\mathrm{off}-\mathrm{nadir}}}{{\mathrm{R}}_{\mathrm{nadir}}}$$

(1)

Since the experimental setup covers two different viewing directions, two ANIFs were obtained: one for forward scatter (ANIF_fs) and one for backscatter (ANIF_bs). Additionally, a third ANIF_fs/bs was calculated as the ratio between reflectances per waveband in forward and backscatter direction.

If spectral information of the same target is available for multiple angles, it was found useful to combine them and thus raise predictability. This has been done with surface-near spectrometers that are handheld [40], mounted on a tower [41], on hemispherical devices [42] or UAVs [11]. A prominent example for multi-angular optical remote sensing from space is CHRIS/Proba which allows to record narrowband spectra in the VNIR-domain from five different viewing angles (e.g., [43,44]). EnMAP will be able to perform an across-track satellite tilt, but will keep up this slanting position for longer time than its view upon the target. If the same target shall be observed under different zenith angles, more than just one acquisition will have to be made with a time gap of at least several days or possibly several weeks or months. Canopy parameters that are strongly influencing the bidirectional reflectance distribution function (BRDF) change diurnally as well as during the seasonal growth cycle. For this reason, in this study we concentrated on single looks only, regardless of the possible improvement of results for a combined multi-angular approach.

2.2.2. Biophysical Variables

Agricultural crop variables were measured at the exact same location where spectral signatures were recorded shortly before. The time offset between the variable sampling and the spectral sampling was 45 min on average and 60 min at maximum.

Average Leaf Inclination Angle (ALIA) was measured with a Suunto PM-5/360 inclinometer held along the leaf petiole to display its slope against the horizontal plane [45]. The measurement was repeated at different positions of the leaf and for different leaves within the canopy. Additionally the Leaf Inclination Distribution Function (LIDF) was noted down for a more detailed description of the canopy architecture [46]. Leaf chlorophyll content (LCC) was measured with a Konica-Minolta SPAD-502 handheld device at different heights with focus on the upper canopy layer. The chlorophyll meter had been individually calibrated in a preceding field campaign against destructive measurements of winter wheat leaf chlorophyll content from different senescence states. Coefficients of [47] were used to derive LCC from the samples. Leaf senescence (C_br) was estimated as the fraction of brown leaf parts in the foliage. This variable varies between zero (no brown spots = 100% fresh vegetation) and one (no green spots = 100% senescent vegetation). For a proper estimation of C_br, the approach of [48] was slightly adapted to incorporate the non-linearity of the vertical distribution of brown leaves. The dissociation factor between upper and lower layer thus created consistent results. This was achieved by applying a cosine function of the brownness in the upper layer to the power of two. C_br can be written as Equation (2):

$${\mathrm{C}}_{\mathrm{br}}=\frac{{\mathrm{br}}_{\mathrm{u}}+{\mathrm{br}}_{\mathrm{l}}}{2}-\frac{\left|{\mathrm{br}}_{\mathrm{u}}-{\mathrm{br}}_{\mathrm{l}}\right|}{2}·\cos ²\left(90°·{\mathrm{br}}_{\mathrm{u}}\right)$$

(2)

with br_u as the fraction of brown leaf parts in the upper and br_l in the lower layer of the canopy. For LAI measurements, a LI-COR Biosciences LAI-2200 instrument was used that had been upgraded with the ClearSky Kit to obtain functionalities of the advanced LAI-2200C. Equipped with a GPS sensor and a white diffuser cap, the device allows for nondestructive measurements of leaf area index under sunlight conditions. To obtain green LAI, the measured LAI value was multiplied with the factor 1 − C_br to exclude the impact of non-photosynthetic vegetation on LAI measurements. Multiplication of leaf variables with the LAI value allows their interpretation on canopy level, e.g., canopy chlorophyll content (CCC).

2.3. Radiative Transfer Modelling

For this study, PRO4SAIL-5B (PROSPECT 5B + 4SAIL) was used which operates based on the input parameters listed in

Table 2. Overview of the PROSAIL parameter ranges for creation of the look-up-table. Biophysical parameters were varied via uniform distribution, with typical min & max values. Parameters that control the sun–target–sensor-geometry were varied in discrete steps, covering all conditions that were observed in the field.

Model	Parameter	Description	Unit	Min	Max
PROSPECT	N	Leaf structure parameter	-	1.0	2.5
	LCC	Leaf Chlorophylla+b content	mg cm−2	0.0	80
	LCarC	Leaf Carotenoids content	μg cm−2	0.0	20
	EWT	Leaf Equivalent Water content	cm	0.001	0.05
	LMA	Leaf Mass per Area	g cm−2	0.001	0.02
	Cbr	Fraction of brown leaves	-	0.0	1.0
SAIL	LAI	Leaf Area Index	m² m−2	0.0	8.0
	ALIA	Average Leaf Inclination Angle	deg	20	90
	Hspot	Hot Spot size parameter	-	0.01	0.5
	Skyl	Ratio of diffuse and total incident radiation	-	0.1	0.1
	γ	Soil Brightness Parameter	-	0.0	1.0
Model	Parameter	Description	Unit	Min	Max	Divisions
SAIL	SZA	Sun Zenith Angle	deg	30	55	6
	OZA	Observer Zenith Angle	deg	−30	30	3
	rAA	relative Azimuth Angle	deg	0	65	14

and described in Section 2.2.2.

Following the suggestion of [1], the Skyl-parameter was kept stable at 0.1. The soil brightness parameter scales the dominance of the bright and dark canopy background in the output signal. By default, standard literature soils are used for this. In this study they were replaced by the brightest and the darkest soil spectrum of the campaign, measured directly at the study fields for each date. The background signal gains more weight in the simulated reflectance for vegetation that is sparse in terms of green LAI. It is important to note that spectral signatures of senescent canopies differ from those of small plants that cover the soil only partially, although situations might result in the same low value for green LAI. For this reason, another background type is introduced that was calculated as the mean senescence signal for ripe wheat crops of both seasons. All other leaf and canopy parameters were randomly drawn from uniform distributions with min and max values adjusted according to Table 2. Input parameters regulating the s–t–s-geometry uniformly covered all field scenarios. For example, the minimum SZA observed in the field was 29.14° and the maximum was 52.63°. As a result, SZAs of 30°, 35°, 40°, 45°, 50° and 55° were used for generating the LUT. Variations in the OZA of −30°, 0° and +30° took account of the three experiments of the simulated EnMAP platform tilt. For winter wheat crops it was suggested setting the leaf structure parameter N to a mean of 2.0 with a SD of 0.34 [7]. Since our study data covered the complete vegetation cycle from seeding to harvest, these values were slightly adapted to a wider range of 1.0 to 2.5.

The size of each LUT (n_lut) can be understood as the number of variations of the parameters (n_para) multiplied by the number of variations of the s–t–s-geometry (n_angles). n_lut has a linear influence on the calculation time for the generation and the inversion of the LUT. On the other hand, larger LUT sizes yield more possible parameter constellations, which may improve the quality of the retrieval. Many authors suggest n_para = 100,000 as the best trade-off between calculation time and inversion accuracy (e.g., [8,12,49]). In each of these studies, however, angles of sun and observer were fixed. As described in Table 2, n_angles here needed to cover 252 different geometric constellations which would result in n_lut = 25,200,000 for each soil and senescence background. A LUT-size of 12,600,000 (n_para = 50,000) turned out to perform equally well (Δ_RMSE < 1%), while allowing a quicker inversion and thus the conduction of more experiments in the same period of time. Accounting for the different potential background signals (soil or senescent material respectively) the LUTs are duplicated and only varied by a different background signal. This method, therefore, shall be called duplex LUT.

Finally, artificial noise can be applied to make simulated spectra more realistic and to improve the inversion accuracy (overview given by [6]). The best performing LUT settings have been varied in noise type (Gaussian additive/Gaussian inverse multiplicative) and noise level (0.0%, 0.1%, 1.0%, 2.0%, 5.0%, 10.0%) respectively.

2.4. Step-Wise Inversion of the LUT

Inverting a LUT means comparing measured spectra with all PROSAIL model results and selecting the parameters that led to the best performing LUT members. Different cost functions can be used to quantify the agreement between measurement and model result. Most authors use the Root Mean Squared Error cost function type (RMSE_cft), defined as Equation (3):

$${\mathrm{RMSE}}_{\mathrm{cft}}=\sqrt{\frac{1}{\mathrm{n}}·{{\displaystyle \sum }}_{\mathrm{i}=1}^{\mathrm{n}}{\left({\mathrm{R}}_{\mathrm{measured}}\left({\mathsf{\lambda }}_i\right)-{\mathrm{R}}_{\mathrm{simulated}}\left({\mathsf{\lambda }}_i\right)\right)}^2}$$

(3)

By squaring the distances before extracting the root, larger deviations gain more influence in this term. Consequently, the RMSE_cft favors results for which both spectral signatures match rather closely for all wavelengths. An alternative cost function tested is the Nash-Sutcliffe-Efficiency (NSE_cft [50]) as defined in Equation (4),

$${\mathrm{NSE}}_{\mathrm{cft}}=1-\frac{{{\displaystyle \sum }}_{\mathrm{i}=1}^{\mathrm{n}}{\left({\mathrm{R}}_{\mathrm{measured}}\left({\mathsf{\lambda }}_i\right)-{\mathrm{R}}_{\mathrm{simulated}}\left({\mathsf{\lambda }}_i\right)\right)}^2}{{{\displaystyle \sum }}_{\mathrm{i}=1}^{\mathrm{n}}{\left({\mathrm{R}}_{\mathrm{measured}}\left({\mathsf{\lambda }}_i\right)-{\overline{\mathrm{R}}}_{\mathrm{measured}}\right)}^2}$$

(4)

Weighing the squared sum of distances between measured and simulated reflectance against the squared sum of distances between measured reflectance and the average measured reflectance of the complete spectrum. The mathematically simplest approach is the mean absolute error (MAE) as defined in Equation (5):

$$\mathrm{MAE}={{\displaystyle \sum }}_{\mathrm{i}=1}^{\mathrm{n}}\left|{\mathrm{R}}_{\mathrm{measured}}\left({\mathsf{\lambda }}_i\right)-{\mathrm{R}}_{\mathrm{simulated}}\left({\mathsf{\lambda }}_i\right)\right|$$

(5)

In all three cases, R_measured(λ_i) is the measured and R_simulated(λ_i) is the modelled reflectance at wavelength λ for the ith spectral sensor band, whereas n corresponds to the total number of bands used for the optimization.

For each sampling date, the s–t–s-geometry must be known. Prior to application of the cost function, the correct sub-LUT must be selected. At first, by analysis of the observed spectrum, the corresponding LUT is inquired, depending on the expected canopy background. Senescent vegetation does not only show distinct absorption features by leaf pigments, but also significant features in the SWIR domain. A new index that has been optimized for the EnMAP spectral configuration, the NPVI_EnMAP (Equation (6)) is introduced. NPVI_EnMAP is used to classify the background of a pixel as either soil (type A) or senescent vegetation (type B) based on a simple threshold.

$${\mathrm{NPVI}}_{\mathrm{EnMAP}}=\frac{{\mathrm{R}}_{2218}}{{\mathrm{R}}_{671}}$$

(6)

If NPVI_EnMAP < 1.4, the spectrum is classified as type B and classified as type A for all other cases. Based on the angular constellation for each pixel a decision is made, which of the remaining 252 sub-LUTs shall finally be used for the inversion.

The ill-posedness can be mitigated by narrowing the parameter constraints for the generation of the LUT. In this case, the user needs to have access to a priori information about the expected data range. These constraints make the approach more empirical and thus inconsistent with the proposed generic conviction of the study. For this reason, the ill-posed problem was dealt with by considering more than just the one best performing LUT member and its according parameter configuration [3]. The final results vary with the number of considered best fits (n_bf). A tradeoff between singular (ill-posed) and multiple (over-balanced) solutions needed to be found for an optimal retrieval setup (n_bf ϵ {1, 20, 50, 100, 200, 500, 1000}). For n_bf > 1, the median is used to get the final parameter value. Figure 2 illustrates the necessity to include an adequate amount of LUT-members for the variable retrieval. Parameter constraints for the creation of the LUT can be narrowed down to further increase inversion performance (e.g., [42]). In doing so, the model is calibrated to site-specific characteristics and might not be able to help retrieve variables for other fields, crop types or phenological stages.

The step-wise hierarchical variable retrieval was achieved by several consecutive complete LUT-inversions. A motivation for this approach is the dominance of some parameters (e.g., LAI) that may suppress the signal of others (e.g., LCC) affecting similar spectral domains. For the first inversion run, all available spectral bands were included except for those influenced by the atmospheric water vapor absorptions (1359 nm–1465 nm and 1731 nm–1998 nm) and the VNIR-bands in the detector overlap of the EnMAP-HSI (911 nm–985 nm). Although all variables were obtained in this first step, LAI was the only one of interest at that time. The average inclination of leaves is an important regulation parameter that describes the visibility of photosynthetically active parts of the vegetation for the sensor. Erectophile canopies reveal larger parts of the underlying soil, especially for low SZAs. Planophile and plagiophile canopies on the other hand cover more of the background and lead to stronger signals just like an increased LAI would. ALIA and LAI therefore counterbalance each other. In an attempt to separate their influences on the measured spectra, another pre-selection is investigated for the first inversion run, selecting only those LUT members with an ALIA close to the one estimated in the field.

For the second inversion run, the LAI values resulting from the first run were fixed. A pre-analysis selected only those LUT members containing the retrieved LAI ± an absolute tolerance of 0.01 (m² m⁻²). If this pre-selection left fewer members than two times the size of n_bf, the valid tolerance was expanded by increments of 0.01 until the minimum condition was met. Only then the second run was started during which LCC was retrieved by separately applying the cost function to the variable-specific sensitive wavelengths (LCC @ 423–705 nm).

According to the authors of [51] the performance of an inversion setting shall ideally be assessed by multiple statistical quality criterions when comparing retrieved model parameters to in situ measured values. Most importantly, the relative Root Mean Squared Error (rRMSE) and the slope of the regression line (m) were considered in this study. The coefficient of determination (R²) played a minor role, since it measures the strength of the correlation according to the linear regression rather than the 1:1 relationship between model parameter and in situ variable. A regression model was calculated nonetheless and its slope served as an indicator of the inversion accuracy. A slope of 1.0 suggests a perfectly outbalanced relationship. Slopes > 1.0 reveal an underestimation for low and overestimation for high values. The reverse relationship applies for slopes <1.0. For all following analyses, model runs with a slope <0.7 or >1.3 were not considered in the final results.

3. Results

3.1. Impact of the Observer Zenith on Reflectance Spectra

Analyzing all 3 × 3 images separately would be time consuming and impractical. For this reason, the spatial dimension was partially sacrificed in favor of an additional temporal dimension by mosaicking all images of one growing season below each other. This principle is explained in Figure 3. The combination of several sub-images in one mosaic allows the application of algorithms for all observations at once and visualizes seasonal changes in reflectance.

A strong dependency of the anisotropy factor (ANIF, Equation (1)) towards wavelength can be observed (see Figure 4). A striking anisotropic behavior of the canopy is evident in the short wave visible range for both illumination settings, with decreasing impact towards longer wavelengths. The experiment was repeated for three different growth stages: (1) Day of year (DOY) 128: stem elongation (flag leaf visible, but still rolled), (2) DOY 182: development of fruit (late milk) & (3) DOY 202: Ripening (fully ripe). Each of these three growth stages is represented by an average of all adjacent canopy spectral signatures of that specific date. For both forward and backscatter the ANIF is highest for phenological stage (1). This indicates that during earlier growth stages angular effects in the observer zenith have the strongest influence on the spectral signal. An ANIF of value 1.0 means identical reflectance for nadir and off-nadir observation. Reflectances of the fully ripe canopy (3) draw closest to this equilibrium line for both viewing directions. ANIF_fs marginally drops below 1.0 for the VIS-domain and from 1500 nm to 2500 nm at growth stage (2), due to lower reflectances in the off-nadir compared to nadir observations. For backscatter, ANIF_bs of (2) is situated mostly between (1) and (3). Lowest impact of angular variations can be assumed for the visible range for forward scatter and for the NIR plateau for backscatter observations. For green vegetation, there is a local minimum in ANIF_bs around 550 nm. On the other hand, maximum anisotropy is observed for shortest (450 nm) and longest (680 nm) visible wavelengths. ANIF_fs shows an anisotropic behavior that is exactly reverse to ANIF_bs.

In the special case of ANIF_fs/bs, differences between negative and positive observation angles can be assessed directly. As expected, backscatter angles generally lead to higher reflectances. The 1.0 line is closest for the NIR plateau of green vegetation signatures. For all three settings, the senescent vegetation had a constant ANIF that was less sensitive to wavelength, but more sensitive to viewing direction.

3.2. Impact of the Observer Zenith on the Retrieval of Crop Parameters

Results for the retrieval of LAI and LCC from different observation angles and two different cost functions are shown in

Table 3. Influence of angular spectral measurements, representing the tilt of the EnMAP satellite platform of ±30° across track (OZA). Model inversion was conducted with relative root mean squared error (RMSE)_cft and mean absolute error (MAE) and different sizes of n_bf. Best results for leaf area index (LAI) and leaf chlorophyll content (LCC) for each angle and cost function type are shown.

OZA	RMSE Cost Function				MAE Cost Function
	LAI		LCC		LAI		LCC
(deg)	Slope	rRMSE	Slope	rRMSE	Slope	rRMSE	Slope	rRMSE
−30	0.81	0.27	0.95	0.20	0.84	0.24	0.82	0.22
0	0.94	0.19	0.87	0.24	0.92	0.18	0.89	0.26
+30	0.83	0.25	0.94	0.27	0.82	0.27	0.79	0.28

. For LAI, the lowest rRMSE of 0.18 was obtained for nadir observation with MAE cost function. For OZA = −30° (forward scatter) rRMSE was 0.24 and 0.25 for OZA = +30° which means an error increase of 30%. The associated scatterplots are shown in Figure 5. Subplot (c) illustrates the instable model inversion for LAI from backscatter spectra. Especially for medium LAI values, PROSAIL suggested a widespread variable range and better results for nadir inversion. LCC retrieval on the other hand improved for negative observation angles (rRMSE = 0.20), as the clusters of high and low chlorophyll content moved closer to the 1:1 line. For the lower reflectances in opposite viewing direction, rRMSE = 0.27 was the best possible result. Slope and intercept of the regression line nearly reached the optimum of f(x) = x, but the scattergram shows a non-linear behavior and a reduced rRMSE.

Figure 6 compares the map of in situ variables as measured at the study fields with the results of the inversion for LAI and LCC. Discrete steps in hue and saturation can be found in both in situ values and model results, indicating that the seasonal dynamics were captured in general. Within the 3 × 3 pixels for each field date, variations could not be fully reproduced by the inversion. For this, model uncertainties would have to be lower than the lateral disparity for the 30 × 30 m sub-plots.

By subtracting parameter estimations from in situ measurements, residuals for the inversion of LAI and LCC are obtained. Figure 7 is an illustration of these residuals as a map. For the first season in 2014, the model overestimated LCC most of the time for all angular settings, especially for the beginning of the season in April and May. In April and for one sampling date in July 2015, higher chlorophyll concentrations were measured than predicted by the model. The latter is especially striking, because estimations in nadir were quite close the week before (ΔLCC_DOY182 = 2.56 μg cm⁻²) and after (ΔLCC_DOY202 = 4.24 μg cm⁻²). Interannual changes in model predictability are represented by the mean standard deviation over all 9 pixels and per season. The results are listed in

Table 4. Mean deviations between estimated and in situ measured LAI and LCC, separated by season and angular setting (OZA).

OZA	Season 2014		Season 2014/2015
	LAI	LCC	LAI	LCC
(deg)	(m² m−2)	(μg cm−2)	(m² m−2)	(μg cm−2)
−30°	0.62	8.43	0.99	7.44
0°	0.47	11.86	0.82	7.17
+30°	0.59	7.38	1.08	11.22

. Both LAI and LCC were more homogeneously predicted in 2014, whereas in 2015 residuals tended to be more variable throughout the growing season. Model inversions from nadir observations were 38% less prone to seasonal effects than those from angular observations. The pattern of the residuals for LAI and LCC do not show any relation and statistics suggest that they are independent of each other (R² = 0.10). Residuals of LAI followed a normal distribution (p > 0.05 for all OZA), but those of the LCC estimation did not (p < 0.01 for all OZA). A shift to negative residuals, i.e., an overestimation of LCC by the model, indicated a slight systematic bias.

Retrieval of CCC proved quite successful (see Figure 8). Concentrations of LCC obtained from the second inversion run were multiplied with estimated LAI values and compared to in situ measured CCC (i.e., measured LAI × measured LCC). PROSAIL tended to overestimate CCC with an intercept of 0.31 g cm⁻² for nadir spectra. The relative rRMSE = 0.37 was higher than for both, LAI and LCC. From forward scatter observations, however, CCC was estimated with an rRMSE of 0.33 meaning an improvement of 12%. Retrieval from backscatter spectra suffered from a weaker estimation of LAI and LCC, leading to a decrease of relative RMSE of 0.40 accordingly.

3.3. Improved Look-Up-Table Inversions

Different techniques were explored to find an optimal setting for the inversion of the PROSAIL RTM. At first, the number of best fits (n_bf) was fixed to 100 and artificial noise was added to the modelled spectra (

Table 5. Experimental results for look-up-table (LUT)-based inversion of LAI and LCC with application of artificial noise (additive and inverse multiplicative) on modelled spectra for n_bf set to 100.

Noise Level	Additive Noise				Inverse Multiplicative Noise
	LAI		LCC		LAI		LCC
σ (%)	Slope	rRMSE	Slope	rRMSE	Slope	rRMSE	Slope	rRMSE
0.0	0.89	0.19	1.50	0.34	0.89	0.19	1.50	0.34
0.1	0.89	0.19	1.49	0.33	0.89	0.19	1.33	0.30
1.0	0.89	0.20	1.40	0.29	0.89	0.19	1.33	0.30
2.0	0.88	0.22	1.30	0.29	0.89	0.19	1.34	0.30
5.0	0.80	0.29	1.24	0.28	0.90	0.19	1.32	0.30
10.0	0.64	0.38	1.68	0.31	0.88	0.20	1.27	0.27

). For the LAI retrieval, performance was better, the lower the level of additive noise was assumed (rRMSE = 0.19 for σ = 0.0%). LCC on the other hand could be better retrieved with noise levels of up to 5% (best rRMSE = 0.27). Inverse multiplicative noise created similar results that were almost independent of the absolute level of σ. Interestingly, in contrast to LCC, for LAI inverse multiplicative noise yielded slightly better results than the additive noise type. A parameter-specific tendency can be assumed. Therefore, by default an inverse multiplicative noise level of σ = 2.0% was used for the first and additive noise of the same level for the second inversion run for all analyses.

In the second experiment, n_bf and the cost function type (cft) were varied (

Table 6. Experimental results for minimization by different cost function types. Deviations between modelled and measured spectra were quantified by either the RMSE_cft or the MAE. Additionally, the number of best fits to be averaged was increased from 1 to 1000 to find the optimal setup. Artificial inverse multiplicative noise is set to a level of 2.0% for all results.

Number of Best Fits	RMSE Cost Function				MAE Cost Function
	LAI		LCC		LAI		LCC
	Slope	rRMSE	Slope	rRMSE	Slope	rRMSE	Slope	rRMSE
1	0.88	0.28	1.61	0.44	0.88	0.28	1.55	0.42
50	0.94	0.19	1.44	0.32	0.93	0.18	1.36	0.32
100	0.89	0.19	1.30	0.29	0.92	0.18	1.18	0.28
200	0.90	0.20	1.26	0.28	0.89	0.18	1.13	0.28
500	0.86	0.21	1.03	0.26	0.87	0.20	0.89	0.26
1000	0.83	0.22	0.87	0.24	0.84	0.20	0.77	0.26

). For LAI, the overall best result was achieved for MAE and 50 ≤ n_bf ≤ 200 (rRMSE = 0.18). The MAE generated 3% lower rRMSE in average. NSE as a cost function did not pass the slope threshold test for any of the experiments and was omitted in the analysis of the results. Also, for LCC, the choice of a merit function seemed to be of minor importance in comparison to the number of best fits taken into account for the retrieval. The more LUT members were considered, the better the performance in terms of rRMSE at the expense of precision at the extreme ranges, leading to a reduced regression slope. RMSE_cft performed slightly better than MAE, again suggesting a parameter-specific behavior of the LUT inversion.

4. Discussion

For the interpretation of directional angular effects on spectral reflectance, the anisotropy factor (ANIF) can be consulted: due to self-shading effects, forward scatter images appear darker than nadir or backscatter images. In the latter case, a greater fraction of incident sunlight is reflected back to the direction of its origin and is consequently missing on the opposite viewing direction. This so-called hot spot effect leads to a spectral saturation and superimposes parts of the signal of leaf constituents. Moreover, ANIF is highly correlated with the magnitude of reflectance itself. If the canopy reflectance is higher, discrepancies increase between nadir and forward scatter but decrease between nadir and backscatter. On the other hand, if more radiation is absorbed or transmitted by the canopy, anisotropy decreases for forward scatter and increases for backscatter. Spectrally, high anisotropy occurs for blue and red portions of the solar spectrum from which leaf chlorophyll mainly absorbs radiation to photosynthesize. This was also found by [52] for both directions, but in our study this could only be confirmed for backscatter mechanisms. This phenomenon may be the reason why it was more difficult for the PROSAIL model to reproduce the measured spectra from this direction, leading to a weaker estimation of LCC from backscatter in comparison to forward scatter spectra.

Our main study objective was to assess the effect of off-nadir observations on the prediction accuracy for leaf and canopy variables, namely LAI and leaf chlorophyll content, as it will have major implications for the user community of future EnMAP data. Generally, for both off-nadir observations, accuracy decreased when estimating LAI: rRMSE = 18% at nadir vs. rRMSE = 25% (backscatter) and rRMSE = 24% (forward scatter). For LCC and CCC, the off-nadir mode with forward scatter yielded highest accuracies with rRMSE = 20% and rRMSE = 33% respectively. Once again, the complex structure of the canopy plays an important role for the output of PROSAIL. Turbid medium assumptions are best met for homogeneous crops with least possible complexity in plant structural traits. Winter wheat is thought to be particularly well suited for a representation through RTMs [19]. Nevertheless, the leaf surface of wheat exhibits anisotropic reflectance that is mathematically described by the BRDF [40]. Backscattering leads to glare effects and thus complicates the retrieval of LAI. In the opposite direction, i.e., forward scatter OZAs, the canopy appears darker which seems slightly better suited for the retrieval of LAI. If EnMAP data is only available for backscatter observations, an inversion will still be successful, but the user will encounter larger uncertainties.

Leaf glint generally leads to a reduced accuracy for the inversion of LCC. Senescent plant material absorbs less of the incident radiation and is more subjected to hot spot effects [53]. This could be confirmed (see Figure 5), as LCC > 40 μg cm⁻² was poorly estimated from backscatter spectral images, but well inverted from forward scatter observations. The findings suggest that LCC is best retrieved, the greater the difference of the zenith angle between sun and observer becomes. CCC acts like a linear combination of LAI and LCC in any statistical analysis. The rRMSE for nadir appears to be comparatively high and the intercept of the linear model is t = 0.20 which is 10.6% of the data average. This constant overestimation is caused by the before mentioned overestimation of LAI and LCC from nadir spectra reducing the models accuracy. An improvement of the retrieval of these two parameters also yields an improvement for CCC, as can be seen for results from the forward scatter observations.

Another focus of this work was to test different LUT-based strategies, while keeping a special focus on the zenithal-angular effects for spectral observations from the future EnMAP sensor. Most commonly, the RMSE_cft is used as a merit function to find the best matching LUT members. In fact, there is only a 3.0% mean difference between RMSE_cft and MAE in the resulting parameter estimation. The n_bf to incorporate in the parameter retrieval has a much stronger impact on the success of the inversion. If we assess only the rRMSE as a statistical measure for the inversion performance, it could be concluded that for LCC there is a steady improvement in predictability for larger n_bf. It should be noted, however, that the slope responds conversely, decreasing for larger n_bf and moving away from the optimal value of 1.0. Lower regression slopes indicate that the range of predicted variables becomes more level, cutting off lowest and highest inversion results. Sehgal et al. [54] found an optimal inversion routine with n_bf = 10% of the LUT-size which, in their case, was 5400 members. For larger n_bf the inversion approaches the expected value of the variable as specified before the creation of the LUT, so RMSE is bound to decrease if field observations served as a reference for the original parameter distribution. For all angular settings, statistics deteriorated for n_bf > 100 or 0.2% of the compared LUT-members. This suggests that the LUT composition was optimally set. Accordingly, n_bf = 100 is considered as the optimal setting in the case of this study.

A comparison of the performance with other studies is generally difficult, due to the exploration of different sensor data, LUT-compositions, inversion techniques, crop types and measurement ranges. However, as example, Atzberger et al. [7] retrieved LAI with an RMSE of 0.83 (m² m⁻²) and CCC with RMSE 0.66 (g m⁻²) from winter wheat spectra by training artificial neural networks on PROSAIL which is roughly in the same accuracy range as our findings for nadir observations.

Different sources of errors and uncertainties in the whole inversion process must be considered as limitations to this study: in situ measurements of biophysical variables, spectral measurements, simulation of EnMAP data, model representation and the inversion scheme. For most variables, in situ errors can be reduced by choosing an adequate sampling scheme with multiple repetitions. The median standard deviation for LAI measurements of two seasons was 0.22 (m² m⁻²) and 3.16 (μg cm⁻²) for LCC. Repetitions of the ALIA estimation in the field revealed a mean error of ±7°. Senescent canopies yielded higher uncertainties for the measurement of most variables. Standard deviation of all EnMAP-end-to-end simulations was σ = 0.013 (Refl.) at the NIR-plateau which is 0.28% of the mean reflectance at this wavelength. In comparison to other error sources, this uncertainty played only a minor role. LAI acts as a scaling factor for the leaf constituents. The reflectance signal is ambiguous for substances of lower concentrations within a dense canopy or substances of higher concentrations in a sparse canopy respectively. The hierarchical approach estimates LAI first, fixates it and then finds the other parameters in consecutive inversion steps. This proved to work well for LCC, but not yet for other parameters. For instance, ALIA could not be estimated despite its high sensitivity throughout the covered spectral range [55]. Early in the development of PROSAIL it was stated that ALIA and LAI are highly correlated and therefore can hardly be separated in the inversion [56]. In fact, if LAI is inverted with an accuracy of 20%, there is no autocorrelation (R² < 0.01) of the ALIA residuals, although the parameter itself could not be retrieved in any acceptable way. For an improvement of leaf pigment estimations, the new version of PROSPECT (Prospect-D [57]) is eagerly awaited for a more detailed representation of the leaf-biochemistry, namely the consideration of anthocyanins.

5. Conclusions

Spectral differences from a change in observation geometry depend on the optical properties of foliage and canopy. Complex architecture does not necessarily lead to an increased anisotropy, but the photochemistry of the leaf does very much so. With a step-wise hierarchical variable retrieval based on the PROSAIL model it was shown that longer path lengths do not lead to a higher accuracy in LAI estimation, but still allow a retrieval of this variable with satisfying accuracy. For LCC, on the other hand, the retrieval accuracy did increase when using off-nadir observations. Overall, accuracies are still in the range of about 20% for LAI and LCC. The look-up-table approach was improved when parameters were inverted hierarchically with educated a priori knowledge about the considered wavelengths. The implementation of a non-photosynthetic vegetation background improved the estimation of biophysical parameters especially for senescent phenology states. We agree with the authors of [58] in the assertion that there is no universal convention about the zenith view angles best suited to retrieve canopy structure from. For the final user, the following findings can be synthesized for winter wheat crops:

• Effects of anisotropy are strongest for early phenological stages and backscatter observations;
• LAI is best estimated from near-nadir observations;
• Optimal results for a retrieval of leaf chlorophyll content is achieved for an observer zenith angle opposite to the sun (forward scatter);
• For both variables (LAI and LCC) feasible results are obtained for all considered EnMAP geometrical constellations.

In summary, the off-nadir capability of the future EnMAP sensor will increase the number of available scenes for the user as well as the probability of achieving continuous time-series acquisitions.