to only those acquired within one week of field data collection. Thus, a substantially greater number of observations were available against which the investigated retrieval algorithms could be assessed. The results of leave-one-out cross-validation indicated that the error associated with this interpolation-based approach was small. The largest absolute errors were observed when observations during the onset of senescence were removed, although similar errors occurred when observations during the onset of greenness were left out in the case of LCC. This has important implications for the timing of field campaigns. It is during these periods, when rates of change in LAI and LCC are highest, that the collection of field data is of most importance. Thus, future efforts should place particular focus on these phenological stages.

#### *4.2. Choice of Retrieval Algorithms for Forest-Related Applications*

Previous evaluations of the SNAP L2B retrieval algorithm over crop canopies have demonstrated good performance. For example, Vuolo et al. [36] reported a strong relationship and high retrieval accuracy for SNAP L2B retrievals of LAI over an agricultural site comprised of maize, onion, potato, sugarbeet, and winter wheat (*r*2 = 0.83, RMSE = 0.32, NRMSE = 12%), while Vanino et al. [39] reported comparable results over a tomato crop (*r*2 = 0.69, RMSE = 0.56, NRMSE = 25%). Similarly, an overall RMSE of 0.98 was reported by Djamai et al. [37] in a recent evaluation over alfalfa, black bean, canola, corn, oat, soybean, and wheat. Nevertheless, previous work has demonstrated poor performance when retrieval algorithms based on one-dimensional RTMs such as SAIL are applied over forest environments [23]. The moderate relationships and retrieval accuracies demonstrated by the SNAP L2B retrieval algorithm in our study support these findings, as do the results of our forward modelling experiments.

When compared to our deciduous broadleaf forest site, the crop canopies investigated in previous studies better conform to the turbid medium assumption adopted by SAIL [23,77]. Due to their dense, leafy, and homogeneous nature, a one-dimensional RTM such as SAIL can provide a good approximation of these canopies, enabling high retrieval accuracies to be achieved by SAIL-based retrieval algorithms. In contrast, deciduous broadleaf forest canopies are more heterogeneous, contain a greater number of woody elements, and are subject to increased crown transmission, foliage clumping, and shadowing. As SAIL does not account for these factors, underestimation of LAI and CCC occurs when the SNAP L2B retrieval algorithm is applied over deciduous broadleaf forest. The issue may be compounded by the ANN-based design of the retrieval algorithm, as ANNs are known to perform unpredictably when inputs strongly deviate from their training data [10].

Given that five out of eight bands were considered to be badly modelled by SAIL, it is possible that improved retrieval accuracies could be obtained using a SAIL-based retrieval algorithm trained on a subset of the best modelled bands. As our objective was to evaluate the SNAP L2B retrieval algorithm available to users, which makes use of all considered bands, such an approach was not explicitly investigated in this study. However, methods to determine and select optimal subsets of bands have successfully been applied to improve retrieval accuracies in previous work [23,28,74–76,78]. For example, in a recent study using INFORM to retrieve LCC from MSI data, Darvishzadeh et al. [28] demonstrated increased retrieval accuracies when using a spectral subset of only MSI's red-edge bands as opposed to the full band set. Similarly, Inoue et al. [29] observed that models utilising a large number of spectral bands did not improve the retrieval of CCC from hyperspectral data when compared to those incorporating just two optimally selected bands. It is worth noting that, unlike Darvishzadeh et al. [28], we did not observe large discrepancies between modelled and observed MSI spectra in the near-infrared shoulder. In contrast, B8A (centred at 865 nm) was the band best modelled by both investigated RTMs in our study.

The INFORM-based retrieval algorithm was characterised by increased overall retrieval accuracies, better reflecting previous studies that have successfully applied MSI and its red-edge bands for forest biophysical and biochemical variable retrieval. For example, using a statistical approach to retrieve LAI from MSI data over a boreal forest site, Korhonen et al. [32] reported a comparable relationship and retrieval accuracy (*r*2 = 0.73, RMSE = 0.60, NRMSE = 20%). Similarly, Darvishzadeh et al. [28] obtained an NRMSE of 33% when using LUT inversion of INFORM to retrieve LCC from MSI data

over spruce stands. Beyond MSI, our results reflect those of Schlerf and Atzberger [40], who trained an ANN for retrieval of spruce LAI from airborne hyperspectral data using INFORM, achieving a similar relationship and retrieval accuracy (*r*2 = 0.73, RMSE = 0.58, NRMSE = 18%). Comparable results were reported by Heiskanen et al. [79], who adopted the semi-physical forest reflectance model PARAS to retrieve LAI from High-Resolution Visible and Infrared (HRVIR) and Enhanced Thematic Mapper (ETM+) data over a boreal forest site. Using ANNs, Heiskanen et al. [79] obtained an RMSE of 0.59 (NRMSE = 25%). In a later study, Schlerf and Atzberger [72] applied LUT inversion of INFORM to retrieve LAI from multi-angular Compact High-Resolution Imaging Spectrometer (CHRIS) data, reporting similar results over beech (*r*2 = 0.57, RMSE = 0.94, NRMSE = 26%) and spruce (*r*2 = 0.51, RMSE = 0.74, NRMSE = 18%).

In terms of phenological variations in the relationship between MSI derived retrievals and field data, the results of our study reflect those reported by Heiskanen et al. [80], who used a statistical approach to retrieve LAI from HRVIR and Hyperion data over boreal forest throughout several seasons. As in our study, Heiskanen et al. [80] reported strong relationships with LAI during the onset of greenness (DOY 116 to 153, *r*2 = 0.69 to 0.74). These were followed by weaker relationships during peak greenness and, in particular, the onset of senescence (DOY 273 to 279, *r*2 = 0.29 to 0.40). An explanation for this pattern is that during the onset of greenness, a wide range of LAI values are experienced, whereas during peak greenness, this range is substantially reduced. The weaker relationships during the onset of senescence may be related to changes in leaf biochemistry, to which the MSI derived LAI retrievals are sensitive, but to which the field estimates of LAI (which do not discriminate between green and senescent leaves) are not.

Although improved overall retrieval accuracy was obtained by the INFORM-based retrieval algorithm in the case of CCC, overestimation occurred during the onset of greenness, while some degree of underestimation persisted throughout peak greenness. In terms of the apparent overestimation during the onset of greenness, a small bias was also observed in the INFORM LAI retrievals. Being the product of LAI and LCC, this may have translated into a larger bias in terms of CCC. One explanation for the apparent underestimation observed during peak greenness could be related to the field measurements of LCC themselves, which, due to the height of the forest, could only be performed on leaves from middle or bottom of the tree crown. The sunlit leaves at the top of the crown may have been characterised by reduced LCC when compared to the measured leaves, as LCC commonly decreases with increased light availability [81]. As such, the observed discrepancies may have resulted from sampling biases in the field data. Finally, while our study analysed the performance of the two retrieval algorithms throughout the growing season, it is important to note that it was restricted to a single site. Further work is required to confirm the applicability of our results over additional sites, including those characterised by different forest types.

#### **5. Conclusions**

Although it is not optimised for forest environments, because of its ease of use and integration within the image processing software, it is expected that many users will adopt the SNAP L2B retrieval algorithm for forest LAI and CCC retrieval as a first port of call. Using field data collected throughout the growing season at a deciduous broadleaf forest site in Southern England, we evaluated its performance and that of an alternative retrieval algorithm optimised for forest environments, trained using INFORM. We also developed and successfully applied an interpolation approach to address the temporal mismatch between field and satellite observations. When compared to previously investigated crop canopies, the SNAP L2B retrieval algorithm appears less well-suited to LAI and CCC retrieval over forest environments. Its moderate retrieval accuracies highlight the importance of selecting an RTM that can accurately describe the structure of the canopy of interest, as do the improvements associated with the INFORM-based retrieval algorithm. This finding is corroborated by the results of our forward modelling experiments: over forest environments such as our study site, SAIL appears less able to reproduce observed spectra than INFORM. Based on these results, for forest-related applications using MSI data, we recommend users seek retrieval algorithms optimised for forest environments.

**Author Contributions:** Conceptualization, L.A.B., B.O.O. and J.D.; methodology, L.A.B.; software, L.A.B.; validation, L.A.B., B.O.O. and J.D.; formal analysis, L.A.B.; investigation, L.A.B.; resources, L.A.B. and J.D.; data curation, L.A.B.; writing—original draft preparation, L.A.B.; writing—review and editing, L.A.B., B.O.O. and J.D.; visualization, L.A.B.; supervision, B.O.O. and J.D.; project administration, L.A.B.; funding acquisition, J.D.

**Funding:** This research was supported by the European Space Agency and a University of Southampton Vice Chancellor's Scholarship.

**Acknowledgments:** The authors thank the European Space Agency and the European Commission for making the Sentinel products and associated tools freely available through the Copernicus programme. The authors are also grateful to the Forestry Commission for granting permission to undertake this study, to Tracy Adole, Andrew Maclachlan, Alan Smith and Qiaoyun Xie for their assistance in the collection of field data, and to Clement Atzberger for providing a copy of INFORM.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **Appendix A**

**Table A1.** Theoretical performance of all trained ANNs (assessed using the testing subset of INFORM simulations). The best performing ANNs, indicated by the lowest RMSE, are marked with \*.

#### **References**

© 2019 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).

## *Article* **A Simplified and Robust Surface Reflectance Estimation Method (SREM) for Use over Diverse Land Surfaces Using Multi-Sensor Data**

**Muhammad Bilal 1,**†**, Majid Nazeer 2,3,**†**, Janet E. Nichol 4, Max P. Bleiweiss 5, Zhongfeng Qiu 1,\*, Evelyn Jäkel 6, James R. Campbell 7, Luqman Atique 8, Xiaolan Huang 1 and Simone Lolli 9,10**

Received: 8 May 2019; Accepted: 2 June 2019; Published: 4 June 2019

**Abstract:** Surface reflectance (SR) estimation is the most critical preprocessing step for deriving geophysical parameters in multi-sensor remote sensing. Most state-of-the-art SR estimation methods, such as the vector version of the Second Simulation of the Satellite Signal in the Solar Spectrum (6SV) radiative transfer (RT) model, depend on accurate information on aerosol and atmospheric gases. In this study, a Simplified and Robust Surface Reflectance Estimation Method (SREM) based on the equations from 6SV RT model, without integrating information of aerosol particles and atmospheric gasses, is proposed and tested using Landsat 5 Thematic Mapper (TM), Landsat 7 Enhanced Thematic Mapper plus (ETM+), and Landsat 8 Operational Land Imager (OLI) data from 2000 to 2018. For evaluation purposes, (i) the SREM SR retrievals are validated against in situ SR measurements collected by Analytical Spectral Devices (ASD) from the South Dakota State University (SDSU) site, USA; (ii) cross-comparison between the SREM and Landsat spectral SR products, i.e., Landsat Ecosystem Disturbance Adaptive Processing System (LEDAPS) and Landsat 8 Surface Reflectance Code (LaSRC), are conducted over 11 urban (2013–2018), 13 vegetated (2013–2018), and 11 desert/arid (2000 to 2018) sites located over different climatic zones at a global scale; (iii) the performance of the SREM spectral SR retrievals for low to high aerosol loadings is evaluated; (iv) spatio-temporal cross-comparison is conducted for six Landsat paths/rows located in Asia, Africa, Europe, and the United States of America from 2013 to 2018 to consider a large variety of land surfaces and atmospheric conditions; (v) cross-comparison is also performed for the Normalized Difference Vegetation Index (NDVI), the Enhanced Vegetation Index (EVI), and the Soil Adjusted Vegetation Index (SAVI) calculated from both the SREM and Landsat SR data; (vi) the SREM is also applied to the Sentinel-2A and Moderate Resolution Imaging Spectrometer (MODIS) data to explore its applicability; and (vii) errors in the SR retrievals are reported using the mean bias error (MBE), root mean squared deviation (RMSD), and mean systematic error (MSE). Results depict significant and strong positive

Pearson's correlation (r), small MBE, RMSD, and MSE for each spectral band against in situ ASD data and Landsat (LEDAPS and LaSRC) SR products. Consistency in SREM performance against Sentinel-2A (r = 0.994, MBE = −0.009, and RMSD = 0.014) and MODIS (r = 0.925, MBE = 0.007, and RMSD = 0.014) data suggests that SREM can be applied to other multispectral satellites data. Overall, the findings demonstrate the potential and promise of SREM for use over diverse surfaces and under varying atmospheric conditions using multi-sensor data on a global scale.

**Keywords:** Landsat 8; surface reflectance; LEDAPS; LaSRC; 6SV; SREM; NDVI

#### **1. Introduction**

Due to the cost effectiveness and ready availability of data, satellite remote sensing is now extensively used for deriving various geophysical parameters at a global scale; but, it mostly depends on accurate retrievals of the surface reflectance (SR), i.e., "the fraction of incoming sunlight that the surface reflects", from the remotely sensed data. SR is thus the most basic remotely sensed parameter in the solar reflective spectral bands (visible and infrared) that is used as an essential input parameter to obtain many parameters including vegetation indices [1], leaf area index [2], burned area identification [3], land cover classification [4], aerosol optical depth [5–7], and water quality parameters [8]. It is estimated from the reflectance received by satellites at top of the atmosphere (TOA). However, TOA reflectance is affected by atmospheric constituents that introduce nonnegligible offset and uncertainty in the satellite data. Therefore, before performing any qualitative or quantitative analysis, it is critical to eliminate atmospheric contributions and accurately estimate the SR of a ground target.

Normally, the elimination of the atmospheric contributions and estimation of SR are performed using image-based and physical methods. Image-based methods such as, dark object subtraction (DOS) [9], the empirical line method (ELM) [10], and histogram matching [11] do not employ any physical parameters, e.g., atmospheric direct and diffuse transmissions, water vapor, and ozone, etc., to estimate SR, as they obtain the required ancillary information (solar and sensor viewing geometry) from the image metadata [12]. The most common physical methods are the Atmospheric Correction (ATCOR) [11], the Fast Line-of-sight Atmospheric Analysis of Spectral Hypercubes (FLAASH) [13], the Image Correction (iCOR, previously known as OPERA) [14], the Framework for Operational Radiometric Correction for Environmental monitoring (FORCE) [15], the Landsat Ecosystem Disturbance Adaptive Processing System (LEDAPS) [16], and the Landsat Surface Reflectance Code (LaSRC) [17]. All of these methods estimate surface reflectance based on physical parameters and precalculated comprehensive lookup tables (LUT) which are constructed by radiative transfer (RT) models [18–21]. The ATCOR, FLASH, and iCOR use the Moderate-resolution Atmospheric Transmission (MODTRAN) RT model [18], FORCE is based on the Simulation of Satellite Signal in the Solar Spectrum (5S) [19], and the LEDAPS and LaSRC applies the vector version of the Second Simulation of the Satellite Signal in the Solar Spectrum (6SV) [20,21] to simulate atmospheric conditions. Such methods are complex in nature compared to image-based methods due to the requirement of several ancillary parameters to simulate atmospheric conditions and correct the data degraded by atmospheric constituents. For example, the 6SV, the most commonly used method in the remote sensing community to eliminate atmospheric contributions and estimate accurate surface reflectance, requires information on vertical profiles of air pressure, water vapor concentration, air temperature, ozone concentration, a digital elevation model (DEM), aerosol optical depth, and an aerosol model as well as the solar and sensor viewing geometry to simulate atmospheric conditions and construct the LUT. The attributed accuracy of 6S, of within 1% [22] is nevertheless subject to inherited errors from the ancillary input parameters [12,23], which are obtained from different sources.

Previous studies have evaluated the various atmospheric correction methods over different types of land covers, atmospheric conditions, and geographical locations [12,24–27]. For example, Nazeer et al. (2014) [12] validated the SR from both image-based (DOS and ELM) and physical methods (ATCOR, FLAASH, and 6S) over sand, artificial turf, grass, and water surfaces using in situ measured SR, and found the 6S to be robust and more accurate for SR estimation compared to the other methods. Nguyen et al. (2015) [24] tested the adequacy of the DOS, FLAASH, and 6S for above-ground biomass (AGB) estimations of the Gongju and Sejong regions of South Korea and found that 6S outperforms the other methods. López-Serrano et al. (2016) [25] compared the ATCOR, COST (cosine of the Sun zenith angle), FLAASH, and 6S for estimating AGB in the temperate forest area of northeast Durango, Mexico, and concluded that the 6S method is more efficient and reliable than other methods. These validation exercises suggested that (i) the physical methods performed much better than the image-based methods, and (ii) the 6S is the most reliable physical method.

To date, operational SR satellite products are available from the Moderate Resolution Imaging Spectroradiometer (MODIS) [28], the Visible Infrared Imaging Radiometer Suite (VIIRS) [29], the Landsat 4–7 Thematic Mappers (TM) and Enhanced Thematic Mapper Plus (ETM+) (LEDAPS) [16], the Landsat 8 Operational Land Imager (OLI, LaSRC) [17], and the Sentinel-2 A/B Multispectral Instruments (MSI) [30]. In addition to the operational SR products from the Landsat 8 OLI and Sentinel-2 MSI sensors, a Harmonized Landsat and Sentinel (HLS) SR product is generated to improve the global coverage for every 2 to 3 days at a spatial resolution <30 m [31]. A number of efforts have also been made to validate the LEDAPS, LaSRC, and S2 MSI SR products under different conditions and reference data sets [12,17,32–38].

As evident from all these publications, an accurate estimation of SR based on RT models requires precise retrieval of AOD and vice versa. In other words, these two parameters, i.e., SR and AOD, complement and depend on each other for their inversions. The available state-of-the-art SR methods based on the RT algorithms do not provide a meaningful solution without incorporating information on aerosol particles and atmospheric gases which is a daunting task given the inherent errors in satellite AOD retrievals, which vary spatially and seasonally across the globe [39,40]. Therefore, the prime intent of this study is to provide a user-friendly SR method, which can easily be applied. This study proposes a new Simplified and Robust Surface Reflectance Estimation Method (SREM) based on the equations of the 6SV RT model that can perform SR inversion without using precalculated comprehensive LUT, and information on aerosol particles and atmospheric gases and furnishes results similar to well-known state-of-the-art methods. The outline of the manuscript is as follow: Section 2 is related to the datasets and selection of the validation sites, Section 3 is based on the research methodology for SREM and evaluation process, results and discussions are described in Section 4, and conclusions of this study are summarized in Section 5.

#### **2. Datasets**

In this study, archived datasets are used from the satellite sensors, i.e., Landsat 5 (L5) TM, Landsat 7 (L7) ETM+, and Landsat 8 (L8) OLI, for the development of SREM. For validation purposes, in situ SR data from South Dakota State University (SDSU: grassland site) in South Dakota (Figure 1) are taken by Maiersperger et al. (2013) [35]. For comparison purposes, LEDAPS and LaSRC SR products are obtained for 11 urban, 13 vegetated, and 11 desert (arid) sites (Figure 1) from January 2000 to October 2018.

**Figure 1.** Sites used for validation and comparison of SREM derived surface reflectance (SR) products for Landsat and Sentinel-2A. Please refer to the "S/N" in Table 2 for labeled values, and WRS-2 (Worldwide Reference System) path/row numbers.

#### *2.1. Satellite Data*

#### Landsat TM/ETM+/OLI

The Landsat TM, ETM+, and OLI sensors were launched in March 1984, April 1999, and February 2013, respectively. Landsat ETM+ scenes acquired since May 30, 2003, have data gaps due to the Scan Line Corrector (SLC) failure and the scenes acquired before this date are defined as non-SLC. All sensors have a spatial resolution of 30 m for the multispectral bands and a revisit time of 16 days. In this study, Landsat TM, ETM+, and OLI data for the visible to near-infrared bands (Table 1) are used in the SR retrieval by applying SREM. The new SR data set is compared against the readily available SR products (i.e., LEDAPS and LaSRC) from the respective sensors.

**Table 1.** Spectral bands, band numbers, and central wavelengths (nm) of the L8, L7, and L5 data used in this study.

The LEDAPS algorithm is applied to process the L5 TM and L7 ETM+ Level-1 products to SR (Level-2), where the SR is derived automatically from the calibrated TOA reflectance using the 6S atmospheric correction method and atmospheric parameters [16], similar to the MODIS SR products (i.e., MOD09 for Terra and MYD09 for Aqua sensors). In contrast to the LEDAPS algorithm, the LaSRC SR product [17] for the L8 OLI sensor includes improved estimation of the atmospheric parameters (pressure, water vapor, air temperature, ozone, and AOD) essential as input to the RT-model-based SR estimation. The atmospheric parameters for the LaSRC and LEDAPS are obtained from the MODIS Climate Modeling Gridded (CMG) data products and National Centers for Environmental Prediction (NCEP) gridded products, respectively, both based on the 6SV RT model.

For this study, the Collection-1 Landsat data are taken from the United States Geological Survey (USGS), Earth Resources Observation and Science (EROS), Center's Science Processing Architecture (ESPA) on-demand interface (https://espa.cr.usgs.gov) as Level-1 and Level-2 products for the period of January 2000 to October 2018.

#### *2.2. In Situ Surface Reflectance Data*

The ground-truth SR data were collected by Maiersperger et al. [35] from the South Dakota State University's (SDSU) grassland site in South Dakota (Figure 1) using the Analytical Spectral Devices (ASD) FieldSpec spectrometer. Generally, this site is used by the Committee on Earth Observation Satellites (CEOS) as a reference site for vicarious calibrations. For this site, a total of 10 L5 and L7 scenes were found coincident with the ASD data (Appendix A) for validation.

#### *2.3. Site Selection for Comparison Purpose*

In order to compare the SREM with the standard algorithms, e.g. LEDAPS and LaSRC, three broad land cover types (urban, vegetation, and arid) are selected based on the global mosaics of the standard MODIS land cover type data product (MCD12Q1). This product describes land cover properties derived from yearly MODIS observations. The primary land cover scheme identifies 17 land cover classes defined by the International Geosphere Biosphere Programme (IGBP), which includes 11 natural vegetation classes, 3 developed and mosaicked land classes, and 3 nonvegetated land classes [4]. The spatially aggregated MODIS Collection 6 land cover product at a spatial resolution of 500 m [41] is obtained from the NASA Earthdata Search (https://search.earthdata.nasa.gov/) for the year 2017. A total of 35 AERONET (Aerosol Robotic Network) and CEOS pseudo-invariant and instrumented sites located in urban/built-up (11), vegetated (13) and arid/desert (11) regions are selected based on the MODIS land cover classification (Table 2). For the AERONET sites, only the L8 OLI sensor data were used (from 1 April, 2013 to 15 October, 2018), while for the CEOS sites, the L5 TM, L7 ETM+, and L8 OLI sensors data are used from January 2000 to May 2012, July 1999 to May 2003 (non-SLC affected data), and April 2013 to 15 October, 2018, respectively. For each selected site, the Landsat scenes are ordered through the ESPA on-demand interface (https://espa.cr.usgs.gov).

**Table 2.** Aerosol Robotic Network (AERONET) and Committee on Earth Observation Satellites (CEOS) sites involved in a comparison between the Landsat surface reflectance (Landsat Ecosystem Disturbance Adaptive Processing System (LEDAPS) and Landsat 8 Surface Reflectance Code (LaSRC)) and SREM.

a AERONET sites (L8 OLI sensor's data was used from 1 April, 2013 to 15 October, 2018). b CEOS sites (L5 TM, L7ETM+ and L8 OLI sensors data were used from 1 January, 2000 to 15, October 2018).

#### **3. Methodology**

#### *3.1. Surface Reflectance Inversion*

In general, SR retrievals are derived from the TOA reflectance that can be simulated based on the following 6SV RT model equation (Equation (1)) for the Lambertian uniform target [17,20,42]:

$$\begin{aligned} &\rho\_{TOA}(\lambda\_{\star}\theta\_{\star\star}\theta\_{\overline{\nu}\star}\rho\_{\tau}\tau\_{\mathfrak{a}\prime}\omega\_{\mathfrak{a}\prime}P\_{A\prime}, \mathcal{U}\_{H\_{2}O}, \mathcal{U}\_{O\_{3}})\\ &= T\mathcal{g}\_{\mathcal{O}G}T\mathcal{g}\_{O\downarrow}[\rho\_{\rm atm}(\lambda\_{\star}\theta\_{\star\star}\theta\_{\overline{\nu}\star}\rho\_{\tau}\tau\_{\mathfrak{a}\prime}\omega\_{\mathfrak{a}\prime}P\_{A\prime}, \mathcal{U}\_{H\_{2}O})\\ &+ T\_{s}(\lambda\_{\star}\theta\_{\overline{\nu}\star}\pi\_{\mathfrak{a}\prime}\omega\_{\mathfrak{a}\prime}P\_{A\prime})T\_{v}(\lambda\_{\star}\theta\_{\overline{\nu}\star}\pi\_{\mathfrak{a}\prime}\omega\_{\mathfrak{a}\prime}P\_{A\prime})\frac{\rho\_{\tau}(\lambda)}{1-\mathsf{S}\_{\text{atm}}(\lambda\_{\tau}\pi\_{\mathfrak{a}\prime}\omega\_{\mathfrak{a}\prime}P\_{A\prime})\rho\_{\tau}(\lambda)}\mathcal{T}\mathcal{g}\_{H\_{2}O}] \end{aligned} \tag{1}$$

where

ρ*TOA* = reflectance received by satellite at the top of the atmosphere,

ρ*atm* = atmospheric intrinsic path reflectance,

λ = wavelength

*Ts* = atmospheric transmittance of sun-surface path (downward),

*Tv* = atmospheric transmittance of surface-sensor path (upward),

ρ*s* = surface reflectance to be estimated,

*Satm* = atmospheric backscattering ratio to count multiple reflections between the surface and atmosphere,

θ*s* = solar zenith angle,

θ*v* = sensor zenith angle,

ϕ = relative azimuth angle,

*UH*2*O*, = the integrated water vapor content,

*UO*3, = the integrated ozone content,

τ*a*, ω*a*, *PA* = aerosol optical depth, aerosol single scatter albedo, and aerosol phase function, respectively, and

*TgH*2*O*, *TgO*3 , *TgOG* = gaseous transmission by water vapor, ozone, and other gases, respectively.

The atmospheric intrinsic reflectance can be approximated using Equation (2) [17]:

$$\rho\_{\text{film}}\left(\lambda,\theta\_{\text{s}},\theta\_{\text{l}},\varphi\_{\text{l}},\omega\_{\text{l}s},\mathcal{P}\_{\text{A}},\mathcal{U}\_{\text{H}\_{2}\text{O}}\right) = \rho\_{\text{R}}(\lambda,\theta\_{\text{s}},\theta\_{\text{l}},\varphi) + \left(\rho\_{\text{A}+\mathbb{R}}(\lambda,\theta\_{\text{s}},\theta\_{\text{l}},\varphi) - \rho\_{\text{R}}(\lambda,\theta\_{\text{s}},\theta\_{\text{l}},\varphi)\right)\mathcal{P}\_{\text{H}\_{2}\text{O}}\tag{2}$$

where

ρ*R* = atmospheric reflectance due to Rayleigh scattering and

ρ*A*+*R* = combined atmospheric reflectance due to Rayleigh and aerosols.

The objective of this study is to perform an SR inversion using an equation based on the 6SV RT model without using aerosol information such as τ*a*, ω*a*, and*PA* (i.e., ρ*A* = 0) and other atmospheric parameters such as, *UH*2*O*, *UO*3, and OG (i.e., *TgOG*, *TgO*3 , and*TgH*2*O* = 1). Therefore, these parameters are neglected on the right-hand sides of Equations (1) and (2), such that the TOA reflectance can be approximated as Equation (3):

$$\rho\_{\rm TOA}(\lambda, \theta\_{\rm s}, \theta\_{\rm v}, q, \tau\_{\rm a}, \omega\_{\rm u}, P\_A) = \rho\_{\rm R}(\lambda, \theta\_{\rm s}, \theta\_{\rm v}, q) + T\_s(\lambda)T\_{\rm v}(\lambda) \frac{\rho\_{\rm s}(\lambda)}{1 - S\_{\rm atm}(\lambda)\rho\_{\rm s}(\lambda)} \tag{3}$$

From Equation (3), ρ*s* for the SREM method is approximated as Equation (4), and for simplicity, Equation (4) is expressed as Equation (5):

$$p\_s(\lambda) = \frac{\rho\_{TOA}(\lambda, \theta\_{ss}, \theta\_{\upsilon\prime}, \varphi, \pi\_{\mathfrak{u}}, \omega\_{\mathfrak{u}}, P\_A) - \rho\_R(\lambda, \theta\_{\mathfrak{s}\prime}, \theta\_{\upsilon\prime}, q)}{(\rho\_{TOA}(\lambda, \theta\_{\mathfrak{s}\prime}, \theta\_{\mathfrak{v}\prime}, q, \pi\_{\mathfrak{u}}, \omega\_{\mathfrak{u}}, P\_A) - \rho\_R(\lambda, \theta\_{\mathfrak{s}\prime}, \theta\_{\upsilon\prime}, q))S\_{\mathfrak{u}\mathfrak{m}}(\lambda) + T\_{\mathfrak{s}}(\lambda)T\_{\upsilon}(\lambda)}\tag{4}$$

$$\rho\_s = \frac{\rho\_{TOA} - \rho\_R}{(\rho\_{TOA} - \rho\_R)S\_{atm} + T\_sT\_v} \tag{5}$$

where ρ*s* = SREM estimated surface reflectance.

It should be noted that SREM SR is different than the Rayleigh corrected TOA reflectance, which can be obtained by simple subtraction of Rayleigh reflectance from the TOA reflectance.

In Equation (4), the TOA reflectance and Rayleigh reflectance is computed using Equations (6) and (7) [43], respectively:

$$
\rho\_{TOA}(\lambda, \theta\_s, \theta\_{\upsilon}, \varphi, \pi\_{\mathfrak{a}}, \omega\_{\mathfrak{a}}, P\_A) = \frac{\pi L\_{TOA}(\lambda, \theta\_s, \theta\_{\upsilon}, \varphi) d^2}{ESIN\_{\lambda} \mu\_{\mathfrak{s}}} \tag{6}
$$

where

*LTOA* = radiance received by satellite at the top of the atmosphere, *d* = distance between the Earth and Sun in the astronomical unit, *ESUN* = mean solar exoatmospheric radiation,

μ*s* = cosine of solar zenith angle, and

λ = wavelength.

$$\rho\_R(\lambda, \theta\_s, \theta\_{v\prime} q) = P\_R(\theta\_{s\prime}, \theta\_{v\prime} q) \frac{\left(1 - e^{-M\pi\_r}\right)}{4(\mu\_s + \mu\_v)} \tag{7}$$

where

*M* = air mass calculated using Equation (8) [43], τ*r* = Rayleigh optical depth calculated using Equation (9) [44], *PR* = Rayleigh phase function calculated using Equation (10) [43], and

μ*v* = cosine of sensor zenith angle.

$$M = \frac{1}{\mu\_s} + \frac{1}{\mu\_v} \tag{8}$$

$$
\pi\_r = 0.008569(\lambda)^{-4} \left( 1 + 0.0113(\lambda)^{-2} + 0.0013(\lambda)^{-4} \right) \tag{9}
$$

$$P\_R = \frac{3A}{4+B} \{ 1 + \cos^2 \Theta \}; A = 0.9587256, B = 1 - A \tag{10}$$

where

Θ = scattering angle, and

*A* and *B* are coefficients that account for the molecular asymmetry.

In Equation (4), the atmospheric backscattering ratio and total atmospheric transmission, without integrating aerosol information, is expressed as Equations (11)–(13) [7,45,46], respectively:

$$S\_{atm}(\lambda) = (0.92\tau\_r)e^{-\tau\_r} \tag{11}$$

$$T\_s(\lambda) = e^{(-\tau\_r/\mu\_s)} + e^{(-\tau\_r/\mu\_s)} \left| e^{(0.52\tau\_r/\mu\_s)} - 1 \right. \tag{12}$$

$$T\_{\overline{v}}(\lambda) = e^{\left(-\tau\_r/\mu\_v\right)} + e^{\left(-\tau\_r/\mu\_v\right)} \left| e^{\left(0.52\tau\_r/\mu\_v\right)} - 1 \right>\tag{13}$$

The SREM SR retrievals were estimated for each Landsat band (TM and ETM+: B1-B5 and B7, and OLI: B1-B7) using Equation (4), and for a clear understanding of the SREM, the step-by-step methodology is described in Figure 2.

**Figure 2.** Systematic methodology of the SREM (Simplified and Robust Surface Reflectance Estimation Method).

#### *3.2. Evaluation Process*

The SREM is validated against a range of criteria and features, extensive in itself, in order to test its robustness and explore its potential application. The evaluation process comprises eight steps: (1) The SREM estimated SR and LEDAPS (TM5 and ETM+) SR observations are compared with in situ SR measurements collected by Maiersperger, Scaramuzza, Leigh, Shrestha, Gallo, Jenkerson, and Dwyer [35]. The SREM and LEDAPS SR retrievals are averaged from the spatial window of 3 × 3 pixels if at least 2 out of 9 pixels are available centered on the measurement site. (2) The SREM and Landsat (LEDAPS and LaSRC) SR retrievals are compared for 35 sites located over urban (2013–2018), vegetated (2013–2018), and desert surfaces (2000 to 2018) (Figure 1 and Table 2). To obtain the collocated SREM and Landsat data, (i) retrievals are filtered for the quality flag "66" (clear and low-confidence cloud) for LEDAPS and "322" (clear and low-confidence cloud) for LaSRC, and (ii) matched for the same time and location. (3) In order to evaluate the performance of the SR inversion methods during different aerosol loadings, SR retrievals for each channel were filtered based on the AOD at 550 nm obtained from AERONET sites. (4) A spatio-temporal cross-comparison is conducted for six Landsat paths/rows (122/44, 199/26, 201/32, 170/78, 15/33, and 25/39) located over different regions and climatic zones. For a comprehensive comparison for diverse land surfaces and varying atmospheric conditions, 3000 data points are randomly selected from each image-pair of SREM and LaSRC for each path/row. For this comparison, those Landsat retrievals, that may have values outside the theoretical limits, i.e., 0 < SR < 1, are removed. These unusual retrievals are available due to over-correction for atmosphere and Landsat calibration errors [47–50] or retrievals with SR > 1 might be available for those surfaces that reflect more strongly than Lambertian surfaces [51]. (5) The Normalized Difference Vegetation Index (NDVI, Equation (14)) [1,52], Enhanced Vegetation Index (EVI, Equation (15)) [1], and Soil Adjusted Vegetation Index (SAVI, Equation (16)) [1] are calculated using SREM and Landsat SR data and compared with each other to demonstrate the ability of the SREM to monitor vegetation and crops.

$$NDVI = \frac{(NIR - Red)}{(NIR + Red)}\tag{14}$$

$$EVI = 2.5 \times \left\{ \frac{(NIR - Red)}{\left(NIR + 6 \times Red - 7.5 \times Blue + 1\right)} \right\} \tag{15}$$

$$SAVI = 1.5 \times \left\{ \frac{(NIR - Red)}{(NIR + Red + 0.5)} \right\} \tag{16}$$

(6) To further explore its applicability, the SREM is applied to the Sentinel-2A and MODIS datasets and compared with the Sentinel-2A SR observations estimated by the latest version (2.5.5) of the Sen2Cor and MOD09 level 2 surface reflectance products, respectively. For this purpose, Beijing, a city with mixed bright urban surfaces that mostly remains under frequent haze and dust pollution effects, is selected as a test site. (7) In order to calculate the slope and intercept between the SREM and Landsat retrievals, the reduced major axis (RMA) is used, which can simultaneously account for errors in both dependent and independent variables [53,54]. In RMA, slope (β) and intercept (α) are determined using Equations (17) and (18):

$$
\beta = \frac{\sigma\_y}{\sigma\_x} \tag{17}
$$

$$
\alpha = \overline{Y} - \left(\frac{\sigma\_y}{\sigma\_x}\right)\overline{X} \tag{18}
$$

where

*X* and *Y* = means of X and Y, respectively, and σ*x* and σ*y* = standard deviations of X and Y, respectively.

(8) To report the consistency and errors in the SREM SR product, the Pearson's correlation coefficient (r), mean bias error (MBE, Equation (19)), root-mean-squared difference (RMSD, Equation (20)), and mean systematic error (MSE, Equation (21)) are computed. The MSE is useful to report the difference between the trend of X and Y data; small MSE indicates a good trend.

$$MBE = \frac{1}{n} \sum\_{i=1}^{n} (Y\_i - X\_i) \tag{19}$$

$$RMSD = \sqrt{\frac{1}{n} \sum\_{i=1}^{n} \left(Y\_i - X\_i\right)^2} \tag{20}$$

$$MSE = \frac{1}{n} \sum\_{i=1}^{n} \left(\hat{Y}\_i - X\_i\right)^2 \tag{21}$$

where *Y*ˆ = predicted value based on RMA relationship (Y = βX + α) between X and Y.

#### **4. Results and Discussion**

#### *4.1. Cross-Comparison of ASD, LEDAPS, and SREM SR Data*

The SR data collected by ASD for the SDSU site are available for only 10 days [35] (Appendix A) and are compared with the LEDAPS and SREM SR retrievals (Table 3). Table 3 shows comparable values of Pearson's correlation coefficient (r) for LEDAPS and SREM with ASD data. Overestimation is observed in LEDAPS from B1 to B3, and underestimation in SREM is from B4 to B7. The maximum positive MBE for LEDAPS is 0.006 for B2 and for SREM it is 0.018 for B1. Similarly, the maximum negative MBE is for B5, with −0.009 and −0.035 for LEDAPS and SREM, respectively. The results for SREM are satisfactory with a high value of r and the reason for large negative values of MBE are investigated in the following analysis.

**Table 3.** Validation summary of LEDAPS and SREM SR retrievals against Analytical Spectral Devices (ASD) FieldSpec spectrometer data.

1 StDev = Standard deviation. 2 CV = Coefficient of variations (StDev/average).

For this purpose, LEDAPS and SREM SR retrievals are compared with the TOA reflectance observations obtained for the same dates (Table 4). The hypothesis of this analysis is that the method, LEDAPS or SREM, with the larger negative value of MBE would be considered as superior, as it represents the greater removal of atmospheric effects. On the other hand, a large positive value represents "under-correction"; hence, the respective method is unable to remove atmospheric effects significantly. The results (Table 4) show that LEDAPS has a larger negative MBE for B1 compared to SREM, whereas, SREM has a larger negative MBE for B2 and B3 than LEDAPS, which indicate the better performances of LEDAPS for B1 and SREM for B2 and B3. For B4 to B7, LEDAPS has larger positive values of MBE than SREM, which might be due to its sensitivity to the absorption by atmospheric gases in the infrared spectral region and shows "under correction (lack of atmospheric correction)" of LEDAPS. These results (Tables 3 and 4) suggest that SREM is less sensitive to the absorption by atmospheric gases and performs within the expected range, as the average values of SR retrievals are less than the TOA reflectance, whereas LEDAPS values are even greater than the TOA reflectance; especially, for B4 to B7. Therefore, it can be concluded that SREM, without integrating aerosol information and a comprehensive precalculated LUT in the inversion, can provide SR retrievals that are comparable with both the ASD and LEDAPS observations.

**Table 4.** Summary of cross-comparison between LEDAPS and SREM SR retrievals and top of atmosphere (TOA) reflectance observations.

#### *4.2. Cross-Comparison between SREM and Landsat SR Retrievals*

A cross-comparison between SREM and Landsat (LEDAPS and LaSRC) SR retrievals is conducted over urban as well as vegetated surfaces from 2013 to 2018, and desert (arid) surfaces from 2000 to 2018. Figure 3 shows the scatter (dashed line = 1:1 line) and line plots (black line = Landsat, gray line = SREM) between SREM and Landsat retrievals. The results are summarized in Table 5 which shows that the SREM SR retrievals for coastal aerosol (Figure 3a–c) and blue (Figure 3d–f) spectral bands over desert (arid) sites (Figure 3d,f) are well correlated with the LaSRC and LEDAPS SR retrievals with r ~ 0.990–0.991 (Table 5), and small values of MBE ~ 0.004–0.022, and RMSD ~ 0.009–0.024, compared with urban (Figure 3a,d) and vegetated sites (Figure 3b,e). The values of MSE ~ 0.000–0.002 indicate only minor differences between the trend of SREM and Landsat (LaSRC and LEDAPS) retrievals independent from the surface type. The performance of the SREM for the blue band is better than for the coastal aerosol band with significant small values of MBE and RMSD. The SREM retrievals appear overestimated as indicated by the positive values of MBE, and these values over urban and vegetated sites are high compared to the results from desert sites and are acceptable according to previous studies [31,32,48]. Overestimation in coastal aerosol and blue bands may be due to the enhanced aerosol extinction and Rayleigh contribution in these wavelengths.

The performance of the SREM for green, red, NIR, SWIR1, and SWIR2 bands over the urban and vegetated sites is robust with r ~ 0.951–1.00, MBE ~ −0.01–0.011, RMSD ~ 0.001–0.012, and MSE ~ 0.00. The low values of MBE, RMSD, and MSE represent (i) very small differences between the SREM and Landsat retrievals, (ii) scatter points close to the 1:1 line, and (iii) minimal differences between the trend of SREM and Landsat SR retrievals. The comparison of all these bands over the desert sites was relatively less reliable with larger differences (MBE and RMSD) due to underestimation in the SREM retrievals, especially for SR > 0.50. The underestimation in the SREM retrievals in comparison with Landsat, especially for the desert surfaces, might be due to the "under-correction" of the Landsat atmospheric correction algorithm as observed in Section 4.1, which describes the cross-comparison between SR retrievals and TOA reflectance observations. Therefore, overall performance of the SREM appears robust, as results show a high consistency in the SREM with very high values of r, and small values of MBE and RMSD which suggest that the SREM renders consistent spatial (i.e., from site to site) and temporal (i.e., from 2000 to 2018) variations in SR as generated by LaSRC and LEDAPS products over heterogeneous surfaces.

**Figure 3.** Cross-comparison between coincident SREM and Landsat (LEDAPS and LaSRC) SR retrievals over urban and vegetated sites from 2013–2018, and desert (arid) sites from 2013–2018. Where, coastal aerosol band = (**a**) urban sites, (**b**) vegetated sites, and (**c**) desert sites; blue band = (**d**) urban sites, (**e**) vegetated sites, and (**f**) desert sites; green band = (**g**) urban sites, (**h**) vegetated sites, and (**i**) desert sites; red band = (**j**) urban sites, (**k**) vegetated sites, and (**l**) desert sites; NIR = (**m**) urban sites, (**n**) vegetated sites, and (**o**) desert sites; SWIR1 = (**p**) urban sites, (**q**) vegetated sites, and (**r**) desert sites; SWIR2 = (**s**) urban sites, (**t**) vegetated sites, and (**u**) desert sites; the black line = Landsat retrievals; the grey line = SREM retrievals; and the dashed line = 1:1 line.

β

**Table 5.** Summary of cross-comparison between coincident SREM and Landsat (LEDAPS and LaSRC) SR retrievals over urban and vegetated sites from 2013–2018, and desert (arid) sites from 2000–2018.

1 LC = Land cover; 2 TP = Time Period; 3 n = Total number of observations; 4 β = Slope; 5 α = Intercept; 6 r = Pearson's correlation.

To investigate the underestimation in the SREM retrievals over the desert sites for green to SWIR2 bands, cross-comparisons between the Landsat (LEDAPS and LaSRC) and SREM SR retrievals and TOA reflectance observations are conducted, similar to that in Section 4.1. The results (Table 6) show that the Landsat retrievals have positive MBE of 0.014 for green, 0.017 for red, 0.025 for SWIR1, and 0.036 for SWIR2 bands, compared with the SREM retrievals which are within the expected range of below or equal to TOA. These results represent "under-correction" of data by the Landsat atmospheric correction algorithms which might be due to their sensitivity to the atmospheric scattering and absorption in the visible and infrared spectral regions, respectively. Therefore, these results suggest that the apparent "underestimation" in the SREM retrievals over desert sites when compared to Landsat, is mainly due to the under-correction (positive bias) of the Landsat retrievals.

**Table 6.** Summary of cross-comparison between Landsat (LEDAPS and LaSRC) and SREM SR retrievals and TOA reflectance observations over desert sites.

*4.3. Impact of Aerosol Particles on SR Retrievals*

In order to evaluate the performance of the SREM method during low to high aerosol loadings, the SR retrievals for each band are filtered based on five levels of AOD at 550 nm obtained from the AERONET sites, i.e., (i) 0.0 < AOD < 0.1, (ii) 0.1 < AOD < 0.2, (iii) 0.2 < AOD < 0.3, (iv) 0.3 < AOD <

0.4, and (v) 0.4 < AOD < 1.1. The results are presented in Figure 4, where different colors represent different levels of AOD. The cross-comparison is summarized in Table 7, showing that the number of coincident retrievals decreases with the increase in AOD levels. Figure 4 shows that most of the scatter points for each band are close to the 1:1 line and hence well correlated with each other with a value of r from 0.881 to 1.00. According to the statistical summary (Table 7), the values of MBE increase with the increase in aerosol loadings for the coastal aerosol and blue bands, which suggests that the accuracy of SR retrievals for these bands are affected by the aerosol loadings. However, no direct or linear relationship between MBE and aerosol loadings was observed for the other bands (green to SWIR2). This suggests that the performance of the SREM improves for longer wavelength bands (green to SWIR), independent of the aerosol load. Interestingly, it is observed for these bands that the MBE and RMSD for high aerosol loading (0.4 < AOD < 1.1) is smaller than for low aerosol loadings (0.0 < AOD < 0.1), which suggests that SREM retrievals are less sensitive to the high aerosol load. Overall, results are significant and robust, showing consistency between the SREM and LaSRC retrievals during low to high aerosol loadings, and these justify the application of SREM, without integrating information of aerosol particles and atmospheric gases, to estimate SR similar to the LaSRC product.

**Table 7.** Statistical summary of cross-comparison between coincident Landsat and SREM SR retrievals for low to high AOD levels.

1 AOD = AERONET AOD at 550 nm; 2 n = Total number of observations; 3 β = Slope; 4 α = Intercept; 5 r = Pearson's correlation.

p

**Figure 4.** Cross-comparison between coincident Landsat (LaSRC) and SREM SR retrievals over AERONET sites for low to high AOD levels from 2013 to 2018. The dashed line is the 1:1 line.

#### *4.4. Spatio-Temporal Cross-Comparison between SREM and LaSRC Data*

For performing spatial cross-comparison between SREM and Landsat SR products, six paths/rows located in Asia (122/44), Africa (170/78), Europe (119/26 and 201/32), and the United States of America (15/33 and 25/39) are selected to represent the diversity of land cover types, and climatic as well as air quality conditions. Figure 5 shows the LaSRC (Landsat 8) and SREM SR displayed as RGB false composites of bands 6, 5, 4, and the results indicate that the SREM yields SR images, which are spatially comparable with the LaSRC SR images. A spatial differences map between SREM and Sentinel-2A SR retrievals for the mentioned paths/rows is added as Figure S1 (supplementary data). A careful visual comparison of any land cover feature, from any panel in Figure 5, exhibits strong alikeness and agreement. These results show that the SREM has the ability to remove atmospheric effects without incorporating atmospheric parameters and a precalculated comprehensive LUT based on the RTM.

**Figure 5.** Spatial comparison between LaSRC (left image in each panel) and SREM (right image in each panel) corrected images for different path/rows including (**a**) 122/44, (**b**) 170/78, (**c**) 201/32, (**d**) 199/26, (**e**) 15/33 and (**f**) 25/39. All images are composed using the "natural looking" false color composite of 654 as RGB.  No stretch/contrast is applied to the images.

Temporal analysis is also conducted to consider diverse surface types and atmospheric conditions by selecting 3000 random points from each image of each path/row. This approach found 13 (122/44), 66 (170/78), 21 (199/26), 52 (201/32), 22 (15/33), and 26 (25/39) image-pairs, from 2013 to 2018. Total numbers of coincident points for analysis are 531,462 for blue, 532,179 for green, 533,102 for red, 533,140 for NIR, 534230 for SWIR1, and 533,717 for SWIR2 bands (Figure 6). Overall, results reveal a very good correlation between SREM and LaSRC SR products, with r close to unity (r = 0.993 to 1.00) and small values of MBE ≤ −0.002 for green, red, and NIR bands. A large value of MBE (0.020) is observed for the blue band, which may be due to the enhanced aerosol extinction and Rayleigh contribution. The SREM retrievals for the SWIR1 and SWIR2 are also correlated well with LaSRC retrievals but a slight underestimation is found as indicated by the negative value of MBE (−0.009 to −0.011). This underestimation may be due to under-correction by the Landsat atmospheric correction algorithm as discussed in the previous sections on cross-comparison with TOA reflectance. This investigation considers Landsat data as a "true and standard" data for cross-comparison which in fact has its own uncertainties due to aerosol retrieval algorithm, cloud contamination, and under or over atmospheric correction [47–49,55]. Overall, all these findings demonstrate the robust promise of SREM to retrieve SR for diverse surfaces and under varying atmospheric conditions, without incorporating aerosol and atmospheric parameters, in good agreement with the Landsat SR product.

**Figure 6.** Cross-comparison between SREM and LaSRC data for six paths/rows located in Asia (122/44), Africa (170/78), Europe (119/26 and 201/32), and the United States of America (15/33 and 25/39) from 2013 to 2018 for a large variety of surface types under varying atmospheric conditions. The dashed line is the 1:1 line and the color bar represents the relative frequency of the coincident points.

#### *4.5. Application of SREM to Derive Vegetation Indices*

Vegetation indices such as NDVI, EVI, and SAVI data are computed using SREM and compared with Landsat vegetation indices (Figure 7) for the urban (2013–2018), vegetated (2013–2018), and desert sites (2000–2018), in order to test the suitability of the SREM data for vegetation and crop monitoring. Results reveal high consistency in the SREM computed vegetation indices NDVI (Figure 7a), EVI (Figure 7b), and SAVI (Figure 7c) compared to Landsat, as most of the observations are found close to the 1:1 line (dotted line) with slope from 0.951 to 1.086, intercept from 0.013 and 0.017, Pearson's correlation from 0.995 to 0.997, and MBE from 0.007 to 0.024. This comparison is worthy, as an error in

the surface reflectance can introduce error in the indices and their potential applications. For example, the SREM SR slightly overestimates in the blue band compared to Landsat, which leads to a larger MBE (0.024) and slopes in the SREM EVI (which uses blue, red, and NIR bands) compared to the NDVI and SAVI, which do not incorporate the blue band. These results show that the SREM SR product is faithful and reliable and can be used for vegetation mapping and monitoring on a global scale.

**Figure 7.** Cross-comparison of NDVI (**a**), EVI (**b**), and SAVI (**c**) data based on the SREM and Landsat (LEDAPS and LaSRC) SR products from 2000 to 2018 for the 35 selected urban, vegetated, and desert sites. The dashed line is the 1:1 line; the grey line = SREM retrievals; the black line = Landsat retrievals which are partially hidden by SREM.

#### *4.6. SREM Implementation in Sentinel-2A and MODIS Data*

To further substantiate the applicability of the SREM, preliminary analyses are conducted by SREM using Sentinel-2A and Aqua-MODIS data for Beijing, a city with mixed bright urban surfaces and under effects of severe air pollution episodes. The SREM is applied to cloud-free green band images of Sentinel-2A at 10 m spatial resolution from 8 January to 18 May, 2017, and MODIS at 500 m spatial resolution for the year 2014 (Figure 8). For comparison purposes, the Sentinel-2A SR images are processed using the latest version 2.5.5 of the Sen2Cor atmospheric correction processor, and the Aqua-MODIS Level 2 surface reflectance swath product (MYD09) is used. Figure 8 indicates that most of the scatter points are on or close to the 1:1 line with a high value of r from 0.925 (MODIS) to 0.994 (Sentinel-2A) and small values of MBE from −0.009 (Sentinel-2A) to 0.007 (MODIS) and RMSD of 0.014 for both. The slope between SREM and Sentinel-2A is less than the slope observed for SREM vs. MODIS due to "under-correction "of the Sentinel-2A SR data by the atmospheric correction algorithm, i.e., Sentinel-2A SR values are greater than TOA reflectance over bright surfaces, whereas, SREM SR values are less than TOA reflectance over these surfaces (Figure S2). Overall, these preliminary results suggest that the SREM has the potential to estimate SR also for other multispectral satellite data.

**Figure 8.** Cross-comparison of SREM vs. Sen2Cor 2.5.5, and SREM vs. MYD09 for the Beijing site. Cloud-free images of Sentinel-2A at 10 m spatial resolution from 8 January to 18 May, 2017, and MODIS at 500 m spatial resolution for the year 2014 are used for cross-comparison. The dashed line represents the 1:1 (y = x) line and the color bar represent the relative frequency of the coincident points.

#### **5. Conclusions**

The prime objective of this study was to develop a new Simplified and Robust Surface Reflectance Estimation Method (SREM) based on the Satellite Signal in the Solar Spectrum (6SV) radiative transfer (RT) model equations, without integrating information on aerosol particles and atmospheric gases. The SREM surface reflectance (SR) retrievals were validated against in situ measurements collected by an Analytical Spectral Devices (ASD) spectrometer, and cross-compared with Landsat (LEDAPS and LaSRC) SR products for diverse land surfaces and varying atmospheric conditions, as well as tested on Sentinel2A and MODIS data products. This study concluded that the SREM is capable of accurately estimating spectral surface reflectance (SR) without incorporating information on aerosol particles and atmospheric parameters, and the SR retrievals are comparable with the SR data collected by the ASD spectrometer as well as those provided by Landsat SR products (LEDAPS and LaSRC) which use the 6SV model. Larger positive values of MBE were observed for coastal aerosol band compared to longer wavelengths, which may be related to increase scattering effects at lower wavelengths. Large negative values of MBE were observed in SREM from green to SWIR2 bands when compared to Landsat, which were mainly due to "under-correction (lack of atmospheric correction)" of data by the Landsat atmospheric correction algorithms when compared to TOA reflectance. The preliminary analysis implies that SREM has a strong potential for augmenting vegetation and crop monitoring and it can be implemented with Sentinel-2A and MODIS data or other multispectral satellite data sets.

**Supplementary Materials:** The following are available online at http://www.mdpi.com/2072-4292/11/11/1344/s1, Figure S1: A spatial difference map between SREM and LaSRC corrected images for different path/rows including (**a**) 122/44, (**b**) 170/78, (**c**) 201/32, (**d**) 199/26, (**e**) 15/33, and (**f**) 25/39. The color bar represents the spatial differences between −0.02 and +0.02; Figure S2: Map shows the "under-correction" of Sentinel-2A SR data over bright surfaces compared to the SREM SR data, i.e., Sentinel-2A SR values are greater than the TOA reflectance. No stretch/contrast is applied to the images.

**Author Contributions:** Conceptualization, M.B. and M.N.; methodology, M.B. and M.N.; validation, M.B., M.N., J.E.N., M.P.B., Z.Q., X.H. and E.J.; formal analysis, M.B. and M.N.; investigation, M.B., M.N., J.E.N., M.P.B., and E.J.; resources, M.B., M.N., J.E.N., M.P.B., Z.Q., E.J., J.R.C., S.L., X.H. and L.A.; data curation, M.B. and M.N.; writing—original draft preparation, M.B. and M.N.; writing—review and editing, M.B., M.N., J.E.N., M.P.B., Z.Q., E.J., J.R.C., S.L., and L.A.; funding acquisition, M.B. and Z.Q.

**Funding:** This research is supported by the Special Project of Jiangsu Distinguished Professor<PHONE_NUMBER>003), the Startup Foundation for Introduction Talent of NUIST (2017r107), the National Key Research and Development Program of China by Jiangsu Chair Professorship (2016YFC1400901), and Jiangsu Provincial Programs for Marine Science and Technology Innovation (HY2017-5). Additional support came from the New Mexico State University College of Agriculture Consumer and Environmental Sciences' Agricultural Experiment Station.

**Acknowledgments:** The authors would like to acknowledge USGS for Landsat data, and principal investigators of AERONET sites for aerosol data. The authors also would like to thank Devin White (Oak Ridge National Laboratory) for MODIS Conversion Tool Kit (MCTK), and Vitor S. Martins (Iowa State University, USA) for data sampling procedure, and Yingjie Li (Jiangsu Normal University, China) for Sentinel-2A data processing.

**Conflicts of Interest:** The authors declare no conflict of interest. The funders had no role in the design of the study; in the collection, analyses, or interpretation of data; in the writing of the manuscript, or in the decision to publish the results.

#### **Appendix A**

**Table A1.** List of the Landsat 5 and 7 images used coincident with the ASD spectrometer data for SDSU site obtained from Maiersperger et al. (2013) [35].

#### **References**

© 2019 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).

## *Article* **A Comparison of Hybrid Machine Learning Algorithms for the Retrieval of Wheat Biophysical Variables from Sentinel-2**

**Deepak Upreti 1, Wenjiang Huang 2, Weiping Kong 2,3, Simone Pascucci 4, Stefano Pignatti 4, Xianfeng Zhou 5, Huichun Ye 2 and Raffaele Casa 1,\***

Received: 8 February 2019; Accepted: 21 February 2019; Published: 26 February 2019

**Abstract:** This study focuses on the comparison of hybrid methods of estimation of biophysical variables such as leaf area index (LAI), leaf chlorophyll content (LCC), fraction of absorbed photosynthetically active radiation (FAPAR), fraction of vegetation cover (FVC), and canopy chlorophyll content (CCC) from Sentinel-2 satellite data. Different machine learning algorithms were trained with simulated spectra generated by the physically-based radiative transfer model PROSAIL and subsequently applied to Sentinel-2 reflectance spectra. The algorithms were assessed against a standard operational approach, i.e., the European Space Agency (ESA) Sentinel Application Platform (SNAP) toolbox, based on neural networks. Since kernel-based algorithms have a heavy computational cost when trained with large datasets, an active learning (AL) strategy was explored to try to alleviate this issue. Validation was carried out using ground data from two study sites: one in Shunyi (China) and the other in Maccarese (Italy). In general, the performance of the algorithms was consistent for the two study sites, though a different level of accuracy was found between the two sites, possibly due to slightly different ground sampling protocols and the range and variability of the values of the biophysical variables in the two ground datasets. For LAI estimation, the best ground validation results were obtained for both sites using least squares linear regression (LSLR) and partial least squares regression, with the best performances values of *R*2 of 0.78, rott mean squared error (RMSE) of 0.68 m2 m−2 and a relative RMSE (RRMSE) of 19.48% obtained in the Maccarese site with LSLR. The best results for LCC were obtained using Random Forest Tree Bagger (RFTB) and Bagging Trees (BagT) with the best performances obtained in Maccarese using RFTB (*R*2 = 0.26, RMSE = 8.88 μg cm−2, RRMSE = 17.43%). Gaussian Process Regression (GPR) was the best algorithm for all variables only in the cross-validation phase, but not in the ground validation, where it ranked as the best only for FVC in Maccarese (*R*2 = 0.90, RMSE = 0.08, RRMSE = 9.86%). It was found that the AL strategy was more efficient than the random selection of samples for training the GPR algorithm.

**Keywords:** LAI; LCC; FAPAR; FVC; CCC; PROSAIL; GPR; machine learning; active learning

#### **1. Introduction**

Accurate quantitative estimation of biophysical variables is of crucial importance for different agricultural and ecological applications. Such variables include leaf area index (LAI), leaf chlorophyll content (LCC), fraction of absorbed photosynthetically active radiation (FAPAR), fractional vegetation cover (FVC), and canopy chlorophyll content (CCC). Their knowledge is valuable, for example, for precision agriculture [1], crop traits monitoring [2], and improved yield prediction and reduction of fertilizer usage [3,4]. The retrieval of vegetation biophysical variables from multispectral optical satellite data has been performed by many studies in the last decades [2,5]. In time, there have been significant developments in the algorithms used for such tasks, with a shift from empirical to more physically-based approaches. Starting from simple parametric regressions between vegetation indices and biophysical variables such as LAI and LCC [6–8], the current state of the art relies on hybrid methods exploiting at the same time radiative transfer models (RTM) and non-linear non-parametric regression algorithms [5,9,10]. The methods based on the relationship of, for example, vegetation indices and biophysical variables by the means of fitting functions, typically use the information provided by two or a few spectral bands. This limits the strength of such methods in today's scenarios where tens or even hundreds of spectral bands are available, respectively, in current super-spectral [11] or forthcoming hyperspectral [12,13] spaceborne sensors.

On the other hand, many of the approaches that have been more recently developed for the retrieval of biophysical variables exploit machine learning regression algorithms (MLRAs). These algorithms have gained wide popularity since they allow the use of more complex models, by requiring only one time-consuming training phase, but allowing their fast application any time thereafter for the retrieval [14,15]. Another advantage of these algorithms is the possibility of training them with full spectral information, thus overcoming issues of band selection or transformation, as in the case of parametric regression methods [16]. These algorithms are adaptive and have the ability to cope with the strong non-linearity inherent in remote sensing data [17,18]. Despite their advantages, these algorithms are computationally demanding, so for some MLRAs it is difficult to carry out a training phase with a large number of samples, though there are some that perform well also when trained with small datasets [19]. Some of the MLRAs are also considered black boxes, e.g., neural networks (NNs) [11,20], in which no insight is given about the physical processes linking the spectral reflectance with the biophysical variables. Others, such as partial least squares regression (PLSR) or Gaussian processes regression (GPR), can instead provide some information on what spectral bands are more relevant. In hybrid methods, MLRA can be trained using physically-based modelling approaches that describe the transfer and interaction of radiation inside the canopy based on physical laws, thus providing explicit relations between the biophysical variables and canopy spectral reflectance. When used in the direct mode, biophysical variables are used as input to the physical based RTM, which in turn simulates the top-of-canopy (TOC) spectral reflectance [9]. Using turbid-medium RTMs it is easy to generate during a short time a large simulated database, with realistic ranges of variation of biophysical variables and their corresponding simulated spectra, within the limits of the assumptions of the RTM on canopy properties. These pairs of biophysical variables and reflectance spectra can be used either in numerical optimization or in look up table (LUT) retrieval methods, by minimizing a cost function expressing the differences between observed (e.g., from satellite data) and RTM model's simulated reflectance and looking up the corresponding values of the biophysical variables [21]. For operational applications, due to the pixel-by-pixel calculations, these algorithms are computationally very demanding. These methods are also strongly affected by noise and measurement uncertainties [11]. As inversion methods of physically-based RTMs, they suffer from the limitation of the ill-posed nature of model inversion [22,23], by which different combinations of canopy variables lead to local minima having similar reflectance spectra [24].

In the hybrid methods of retrieval, the database obtained from RTM simulations can be used to train MLRA which are able to establish complex non-linear, non-parametric models linking the biophysical variables and the spectral reflectance [9,11]. Hybrid methods have found widespread application and have reached the operational stage, in particular with the application of NNs [25–27] trained with the RTM PROSAIL [28,29], such as the algorithm implemented in Sentinel Application Platform (SNAP) biophysical processor tool [29], developed by the European Space Agency (ESA). Leaf area index retrieval using NNs trained with PROSAIL is also carried out operationally for Cyclopes and Moderate Resolution Imaging Spectrometer (MODIS) products [30]. However, NN training is a delicate phase and requires the tuning of multiple parameters, which greatly impacts the robustness of the approach [18]. Thus, in recent years, increasing attention has been paid to alternatives to NNs that are simpler to train and have better potential for retrieval accuracy, such as kernel-based methods [31]. These methods solve non-linear regression problems by transferring the data to a higher dimensional space by the means of a kernel function [16]. Kernel-based algorithms have been suggested to offer some advantages in comparison to NNs. Kernel ridge regression (KRR) has proven to be simple for training and to yield competitive accuracy [18]. Gaussian process regression (GPR) is partially transparent compared to the black box nature of NNs. It allows the use of different kernel functions, ranging from simple to very complex, and it also provides uncertainty estimates with the mean value of prediction [10,32]. Gaussian process regression and KRR have been used to estimate successfully LAI and LCC [10,33]. Unfortunately, some kernel-based algorithms, such as GPR and KRR, are computationally very expensive if trained on large sets of simulations [9,10]. In order to have a general-purpose database, including a wide range of vegetation types, generally a huge number of simulations are performed [10]. However, not all pairs of the reflectance and corresponding biophysical variables will be relevant. Few attempts have been made by the researchers to optimize the simulations that are generated by the physically-based models [10,34]. Active learning (AL) has been proposed as a useful strategy to reduce the size of the RTM-generated database, to make the training of the kernel-based algorithms such as GPR more feasible [10]. Active learning is a sub-field of machine learning, also called optimal experimental design in statistics [35]. It initially starts with a small subset of the samples, and then, based on query strategies using either uncertainty or diversity measurement criteria, adds iteratively new samples to the initial training set of samples. In this way, the most informative samples in a dataset are selected, avoiding redundancy [10]. Active learning techniques have a great potential to optimally sample sets generated by RTM.

Very few studies [9,10] have reported the potential of kernel-based methods using AL techniques in relation to the retrieval of biophysical variables such as LAI, LCC, FVC, FAPAR, and CCC. In particular, although the estimation of these variables from hyperspectral data has been demonstrated [28], the suitability of multispectral satellites has yet to be fully proven, especially for biochemical variables such as LCC. The European Space Agency Sentinel-2 mission has been shown to provide data of a high radiometric quality [36] and has a higher revisit frequency than what is planned for hyperspectral satellites in the near future. Because of the availability of a larger number of spectral bands than other multispectral satellites, and of the inclusion of the red edge region of the spectrum, Sentinel-2 has a good potential for the estimation of these variables [9].

This work explores the application of kernel-based GPR using AL for biophysical variables retrieval from Sentinel-2, comparing the potential of this algorithm with other MLRAs and in particular of the version implemented operationally in SNAP.

The main objectives of this study are thus: (1) to compare the performances of different MLRAs, in particular with respect to the algorithm based on NNs implemented in the biophysical processor tool of the ESA SNAP toolbox, by using the same database of PROSAIL simulations [29]; (2) to assess the accuracy of estimation of the biophysical variables LAI, LCC, FAPAR, CCC, and FVC in the wheat crop, for kernel-based (GPR) and non-kernel-based MLRA hybrid methods of retrieval; (3) to explore the feasibility and potential of the use of AL strategies to optimally sample redundant PROSAIL simulations, to minimize computational time and complexity and allow the use of computationally demanding MLRAs (such as GPR) in hybrid methods.

#### **2. Materials and Methods**

The methodology adopted in this work (Figure 1), for the comparison of hybrid methods of biophysical variables retrieval, employed simulations carried out with the model PROSAIL [37] for training MLRA and finally applying them to Sentinel-2 data for assessment with ground measurements collected in two study areas, respectively in Italy and China. The procedure involved the following steps, illustrated in detail in the following sections: (1) generation of simulations with the model PROSAIL; (2) training MLRAs with model simulations and cross-validation; (3) Sentinel-2 data pre-processing; and (4) validation using ground data.

**Figure 1.** Flow-chart of the methodology used for the present study. Retrieval of biophysical variables: Leaf area index (LAI), Leaf chlorophyll content (LCC), Fraction of vegetation cover (FVC), Fraction of absorbed photosynthetically active radiation (FAPAR) and Canopy chlorophyll content (CCC). Radiative Transfer Model (RTM) parameters abbreviations are reported in Table 1. The RTM was used to train Machine Learning Algorithms (MLRA) including Gaussian Process Regression (GPR) using Active Learning (AL). Accuracy was assessed using the root mean squared error (RMSE) and relative RMSE (RRMSE).

#### *2.1. Generation of PROSAIL Simulations*

The model PROSAIL [37] is a widely used RTM for generating realizations of TOC canopy reflectances by considering measured biophysical variables of the crop of interest. It is a combination of the leaf level model PROSPECT [38] and canopy level model SAIL [39]. It has been inverted over a large range of vegetation canopies [40–46] with varying degrees of success. Reflectance and transmittance of a leaf from the PROSPECT model was used as an input for the SAIL model with the soil optical properties and illumination and observation geometry.

The benchmark against which the different hybrid methods were assessed in this study was the NN algorithm implemented in the ESA SNAP Sentinel-2 Toolbox, as biophysical processor. For this reason, PROSAIL simulations were generated using the same input parameter values and configuration as used for SNAP (M.Weiss personal communication), by following the Sentinel-2 Toolbox level 2 products Algorithm Theoretical Basis Document (ATBD) [29], with the exception of the geometry of observation and illumination which was set only for the conditions of our validation dataset (Section 2.3). The PROSAIL parameters were sampled using a fully orthogonal experimental plan [47], using a component of the BV-NNET (Biophysical Variables Neural Network) tool developed by Reference [48]. This procedure consists of identifying classes of values for each variable. Then all the combinations of classes are sampled once. Finally, the actual values of each variable are randomly drawn within the range of variation defined by the corresponding class, according to the distribution law specified for the variable considered. This process allows to take into account all the interactions among parameters, while having the range of variation for each variable densely and near randomly populated [29].

Table 1 reports the parameter ranges, statistical distribution, and number of classes of the PROSAIL model that were used in the present work, which were the same as those employed for training NNs in SNAP (M. Weiss personal communication). A total of 41,472 simulations were generated, representing an extensive range of vegetation properties.

The PROSAIL model top of the canopy (TOC) full spectra (at 1 nm resolution) were then resampled, using Sentinel-2 spectral response functions, to the same eight bands as used in the Sentinel-2 level 2 products ATBD [29], i.e., bands 3 to 7, 8a, 11, and 12, excluding the bands most affected by atmospheric attenuation. The spectral resampling was carried out according to Equation (1):

$$\mathfrak{p} = \sum\_{i}^{N} \left[ \frac{\mathfrak{f}(\lambda)}{\sum\_{i}^{N} \mathfrak{f}(\lambda)} \right] \* \mathfrak{p}(\lambda) \tag{1}$$

where ρ and ρ(λ) are, respectively, the resampled Sentinel-2 reflectance and the PROSAIL reflectance at wavelength λ. β(λ) represents the weight of the band's spectral response function of the Sentinel-2 MSI sensor, whereas i and N are, respectively, the minimum and maximum wavelengths of the Sentinel-2 band. Similarly, a Gaussian noise model was used to better describe actual Sentinel-2 characteristics as well as the ability of the RTM used to represent actual reflectance. The noise was computed as follows [29]:

$$\mathbf{R}^\*(\lambda) = \mathbf{R}(\lambda)(1 + \mathbf{M}\mathbf{D}(\lambda) + \mathbf{M}\mathbf{I})/100) + \mathbf{A}\mathbf{D}(\lambda) + \mathbf{A}\mathbf{I} \tag{2}$$

where R(λ) and R\*(λ) are the raw simulated reflectance and reflectance contaminated with noise, respectively. Multiplicative wavelength dependent noise and multiplicative wavelength independent noise are represented as MD and MI. Similarly, AD and AI are the additive wavelength dependent noise and independent noise, respectively. A value of 0.01 was used for AD and AI and a value of 2% was used for MD and MI for all the bands.

**Table 1.** Input parameter value ranges, number of classes, and statistical distributions used for generating the training database with the PROSAIL model.

#### *2.2. Training Machine Learning Algorithms*

Some of the most widely-used parametric and non-parametric regression algorithms were trained (or calibrated) with the simulations generated by the PROSAIL model, including linear regression, non-linear regression, hierarchical tree-based approaches, and kernel-based algorithms (Table 2).

The non-kernel-based regression methods evaluated in this study were: least squares linear regression (LSLR), partial least square regression (PLSR), NNs, regression trees (RegT), bagging trees (BagT), boosting trees (BooT), random forest fit ensemble (RFFE), and random forest tree bagger (RFTB). We also assessed the kernel-based MLRA GPR, by combining it with an AL strategy. A brief description and the references where these methods are described in detail are reported in Table 2.

#### **Table 2.** Summary of the algorithms compared in this study.

The simulation dataset of 41,472 spectra was used for training the different algorithms. Because of the long processing time required for the training, and of the necessity of having replicates in order to carry out a statistical analysis of the results, 10 subsets of 2500 samples were extracted randomly (with replacement) from the PROSAIL simulated spectra dataset. All the MLRA were trained (and subsequently validated) ten times with the subsets of 2500 simulations each, except for GPR. Since the latter was more computationally demanding, we adopted an AL procedure for selectively further reducing the size of the training set.

Active learning (AL) is a technique which starts from a small set of initial training data and optimally selects samples from a larger data pool, based on different criteria, such as uncertainty or diversity indicators. Samples from the data pool are picked up and added to the initial training dataset, based on diversity criteria. Selecting samples by diversity ensures that added samples are dissimilar from those already accounted for. To identify the best method for selecting the samples from the data pool, we have used three different diversity-based criteria which performed as the best in a previous study [10]. These are angle-based diversity (ABD), Euclidean distance-based diversity (EBD), and cluster-based diversity (CBD). In the implementation used in this work, these criteria were applied also to the label information (i.e., the biophysical variables), since this was available in the data pool (Verreslt J., personal communication). Each process in AL was started with 50 samples as initial training data and the remaining 2450 samples in a data pool. Subsequently, a maximum of 50 samples were added per iteration, with stopping criteria of 100 iterations or a RMSE decrease lower than 50%. Only when performance was improved were the samples added. The whole AL process was repeated 10 times starting from the different subsets of 2500 records each.

In the comparison of the different diversity-based criteria, the full training set of 2500 samples was also used as a reference. Following this assessment, since EBD was the best performing method (see Results section), we only used EBD to further train and validate GPR with AL.

For all the algorithms described above, a cross-validation was performed with a k-fold strategy, where k is the number of folds with a value of 10.

The implementation of all the MLRAs was performed with the Matlab ARTMO Toolbox version 3.24 [58].

#### *2.3. Validation with Sentinel-2 Data and Ground Measurements*

Once the regression algorithms described in the previous section were trained using the PROSAIL simulations, they were applied to Sentinel-2 data acquired over two sites, Maccarese (Central Italy) and Shunyi (China). For both sites, cloud-free Sentinel-2 images were selected for the closest date to ground measurement collection (Table 3). Only the same 8 Sentinel-2 bands that had been used to resample PROSAIL simulations, i.e., bands 3 to 7, 8a, 11, and 12, were used. This band selection was established as optimal for biophysical variables retrieval, by reducing the noise introduced by atmospheric effects [29,59]. The images were initially processed for atmospheric correction with the Sen2cor tool implemented in SNAP (version 6.0) to obtain level 2A TOC reflectance and then resampled to a 10-m pixel size.

Maccarese (41.833◦ lat. N, 12.217◦ long. E, alt. 8 m a.s.l.) is located on the west coast of Central Italy, near Rome. It is a private farm of 3200 ha in a flat area with large fields. Field campaigns to measure biophysical variables on the durum wheat (*Triticum durum* Desf.) crop were carried out for this location from January 2018 to April 2018, for different growing conditions of the crop of interest, at dates close to Sentinel-2 acquisitions. The variables were sampled according to an elementary sampling unit (ESU) scheme, to capture the variability within and among different fields. Each ESU consisted of a quadrat of 20 m by 20 m size, to easily accommodate the Sentinel-2 10-m pixel resolution. A total of 15 ESUs, placed at different locations were employed at different sampling dates. Each of the ESUs contained nine points, where LAI and LCC measurements were collected. Leaf area idex was measured using the LAI-2000 or the LAI-2200C Plant Canopy Analyzer (LI-COR, Lincoln, NE, USA). Since LAI 2000/2200C measurements were not always acquired under diffuse light conditions, data were pre-processed by applying the recommended scattering correction procedure [60]. Four readings were taken above the canopy, with and without diffused cap and six readings below the canopy with 90◦ cap. While processing for scattering corrections five angles of directions 7◦, 23◦, 38◦, 53◦, and 68◦ were considered, i.e., LAI measurements including the extreme rings.

**Table 3.** Sentinel-2 acquisitions and ground measurements dates.

Chlorophyll measurements were obtained using the Force-A Dualex leaf clip reader on the top-most leaves. The calibration of Dualex on durum wheat was previously performed [61] with an estimation accuracy with root mean square error (RMSE) values ranging between 7 and 11 μg cm−2. The following equation [61] was used to recalibrate and compute LCC from the Dualex measured readings

$$\text{LCC} = -3.12 + 1.55 \times (\text{Dualex}) \tag{3}$$

LAI 2000/2200C measurements were also used to derive other biophysical variables, i.e., FVC and FAPAR. FVC was obtained from the LAI values measured on the ground by using the following equation [62]:

$$\text{FVC} = \left(1 - \mathbf{e}^{-0.43 \times \text{LAI}}\right)^{0.52} \tag{4}$$

Previous studies have shown that the instantaneous FAPAR value at 10:00 (or 14:00) solar time is very close to the daily integrated value under clear sky conditions [29,48]. A Matlab script was written to compute FAPAR from the raw LAI-2200 data measured on the ground. The script extracts time dimensions, angles, and gaps fraction data from raw data files, computes day of the year and calculates FAPAR at 10:00 hour solar time. The CCC was estimated by multiplying LAI with chlorophyll content at leaf level.

A dataset of ground measurements corresponding to Sentinel-2 data was also collected in the Shunyi district (40.130 lat. N, 116.654 long. E), Beijing, in China. Ground campaigns were carried out in the period between April 2016 and May 2016 on the winter wheat (*Triticum aestivum* L.) crop. Leaf area index was measured using Licor LAI 2000/2200C and chlorophyll was measured using Force-A Dualex readings and calibrated using laboratory analysis. From April 2016 to May 2016, at four different dates, with 24 observations per date, LAI and LCC were collected in different wheat fields. The FVC was computed as mentioned using Equation (4). Due to the unavailability of raw LAI data, FAPAR was not computed for this test site.

The ground validation dataset was thus composed of 45 samples for Maccarese (each one the average of an ESU) and 96 samples for Shunyi.

#### *2.4. Comparative Assessment of MLRAs*

We used *r*2, root mean square error (RMSE), and relative RMSE (RRMSE %) to assess the accuracy of the retrievals for both the k-fold cross validation and the ground validation. The RRMSE was used to compare the performances across different MLRAs and variables [18]. In general, the performances of the models were estimated by comparing the differences among RRMSE of the estimated and measured values of different biophysical variables. Lower RMSE values corresponded to higher accuracy of the algorithms for the retrieval of biophysical variables from Sentinel-2.

To assess the presence of statistically significant differences among the performances of the models, we applied a Friedman test to the computed metrics [63]. Whenever the Friedman test revealed significant differences (*p* < 0.001) among the metric means, which indicated the presence of significant differences among the models, a Friedman's aligned ranks post-hoc test, followed by Benjamini/Hochberg (non-negative) adjusted for *p*-values, was performed for multiple pairwise comparisons.

#### **3. Results**

The full simulation dataset generated with PROSAIL was used as a source of subsets for training the different MLAs, as described in the methods section, except for GPR for which a data reduction procedure based on AL was employed. Preliminarily, for GPR we compared three different diversity-based criteria used to perform AL, with respect to the full training with 2500 samples. Figure 2 shows the comparison of the different methods of sample selection in terms of *R*2 and RMSE with respect to the number of iterations, for all the variables, applied to the Maccarese dataset. Only when the performance was improved were the samples added. This explains why fewer samples were added than potentially possible, i.e., less than 50 samples were added for each iteration (e.g., in Figure 2a the final sample size for CBD was of 1347 at 75 iterations).

The Euclidean distance-based diversity metric EBD surpassed ABD and CBD, achieving higher *r*2 and lower RMSE values, i.e., higher accuracy and lower error, with a lower number of samples and iterations. The best performing EBD was thus subsequently used in all further applications of AL in the present work. It should be noted that the training of GPR with the full set of 2500 samples, i.e., without AL, provided worse performance than AL in all cases except for FAPAR, despite using a larger dataset for training, revealing some redundancy in the full set.

**Figure 2.** *Cont*.

**Figure 2.** Comparison of three diversity criteria of data pool samples selection in AL, i.e. angle-based diversity (ABD), Euclidean distance-based diversity (EBD), and cluster-based diversity (CBD) for the different biophysical variables for the Maccarese dataset. (**a**,**f**) LAI; (**b**,**g**) LCC; (**c**,**h**) FVC (**d**,**i**) FAPAR; (**e**,**j**) CCC. The values in parentheses in the legend indicate the final size of the training dataset after the Active learning process. The star symbol indicates the final samples after reaching the AL stopping criteria (see text). The horizontal dashed line indicates the performance of the GPR training with the full dataset of 2500 samples without AL.

In the comparison to all the algorithms tested, the cross-validation metrics give an indication of their performance, though the actual accuracy can be more realistically assessed using ground validation results. Tables 4 and 5 report the summary of RMSE estimates for both cross-validation and ground validation for the Maccarese and Shunyi datasets, respectively, for all the algorithms tested.

**Table 4.** Mean RMSE values of the 10 replicates for the estimation of LAI, LCC, CCC, FVC, and FAPAR using different machine learning algorithms (MLRAs) for cross-validation (CV) and ground validation (GV), for the Maccarese test site. Abbreviations for the MLRAs are reported in Table 2. Values sharing the same letters are not significantly different according to Friedman's aligned rank post-hoc test. The best results for each variable are highlighted in bold.

**Table 5.** Mean RMSE values of 10 replicates for the estimation of LAI, LCC, CCC, FVC and FAPAR using different machine learning algorithms (MLRA) for cross-validation (CV) and ground validation (GV), for the Shunyi test site. Abbreviations for the MLRAs are reported in Table 2. Values sharing the same letters are not significantly different according to Friedman's aligned ranks post-hoc test. The best results for each variable are highlighted in bold.

For Maccarese (Table 4), LSLR and PLSR were the best performing algorithms for LAI retrieval when validated with ground data, in terms of RMSE, with values of 0.68 (*R*2 = 0.78, RRMSE = 19.48%) and 0.69 (*R*2 = 0.78, RRMSE = 19.84%) respectively. These were also the fastest computing algorithms among all the MLRAs in terms of time required for training (data not shown). The cross-validation and ground validation results (Figure 3a) indicated that LSLR apparently did not show saturation, even up to LAI values of 5 or 6. The retrieval of LCC was best performed by RFTB (Figure 3b) and BagT, with RMSE values of 8.88 μg cm−2 (*R*2 = 0.26 and RRMSE = 17.43%) and 8.90 μg cm−2 (*R*2 = 0.27 and RRMSE = 17.46%). An overestimation of LCC was observed (Figure 3b). The lowest RMSE of 40.44 g cm−2 (*R*2 = 0.74 RRMSE = 22.7%) (Figure 3e), was obtained for the retrieval of CCC with PLSR. The retrieval accuracy of LSLR was not significantly different from that of PLSR for CCC (RMSE = 40.86 g cm−2, *R*2 = 0.74 and RRMSE = 22.9%). Indeed, PLSR was not significantly different from LSLR in terms of RMSE for all the variables considered. The lowest RMSE for FVC was 0.08 (*R*2 = 0.9 and RRMSE = 9.8%), obtained with the GPR algorithm (Figure 3c). For FVC, RFTB and BagT were also not significantly different from GPR in terms of RMSE. The FAPAR was best retrieved with a RMSE of 0.10 (*R*2 = 0.41 and RRMSE = 12.06%) using the RFTB approach (Figure 3d), although GPR was not significantly different in terms of RMSE, which was only slightly higher. It was observed

that, even in the best performing algorithms, an underestimation of low values of FVC and FAPAR was apparent (Figure 3c,d). It should also be noted that these latter variables show a smaller range of variation in the measured values as compared to the other biophysical variables. As can be seen, GPR was always the best performing algorithm in the cross-validation results, but this was not confirmed by the ground validation tests, with the exception of FVC and FAPAR estimation.

**Figure 3.** *Cont*.

**Figure 3.** Best performing algorithms for the retrieval of (**a**) LAI, (**b**) LCC, (**c**) FVC, (**d**) FAPAR, and (**e**) CCC for the Maccarese test site. The horizontal error bars represent the standard deviation of measurements inside each ESU. Vertical error bars represent uncertainty estimates using GPR for FVC (**c**).

As can be observed from the horizontal error bars of Figure 3, considerable variability occurred among ground measurements inside ESUs, despite the visually apparent homogeneity of the crop canopy in the sampled area of each ESU.

For the Shunyi datasets (Table 5), the ground variability could not be reported for LAI, for which single-point measurements were carried out, but an estimate of the variability for chlorophyll could be made since multiple measurements were available for each point (Figure 4b).

**Figure 4.** *Cont*.

**Figure 4.** Best performing algorithms for (**a**) LAI, (**b**) LCC, (**c**) FVC, and (**d**) CCC retrieval for the dataset from Shunyi, China. The horizontal error bars for LCC represent the standard deviation of single chlorophyll measurements.

The best accuracy for LAI retrieval for the Shunyi site was a RMSE of 1.00 m2 m−2 (*R*2 = 0.73 and RRMSE = 24.12%) obtained with LSLR (Table 5), i.e., the same algorithm providing the best results also for Maccarese. For LAI, PLSR had a similar performance as that of LSLR, with a RMSE = 1.01 m2 m−2 and (*R*2 = 0.72 and RRMSE = 24.50%), i.e., also consistently with the results of LAI retrieval for Maccarese. Similarly, for Maccarese, LCC was best retrieved with an RMSE value of 16.77 μg cm−2, (*R*2 = 0.41 and RRMSE = 40.19%) using RFTB (Table 5). Although RFTB produced the lowest error for retrieving LCC, its retrieval using BagT approach was not significantly different, with RMSE = 17.28 μg cm−2, *R*2 = 0.41 and RRMSE = 41.41%. For CCC retrieval, the lowest RMSE, 56.51 g cm−2 (*R*2 = 0.7, RRMSE = 31.63%), was achieved by the RFTB method. In this case, BagT was not significantly different from RFTB (*R*2 = 0.71, RMSE = 56.72 and RRMSE = 31.75%). For this test site, FVC was retrieved with an RMSE of 0.17 (*R*2 = 0.73 and RRMSE = 23.65%) with NNs (Figure 4c). Also, for the Shunyi test site, GPR was always the best performing algorithm in cross-validation but not in ground validation tests.

A summary of the results in terms of RRMSE for all the biophysical variables tested with different MLRAs for both sites is presented in Figure 5.

**Figure 5.** RRMSE (%) for LAI, LCC, FVC, CCC, and FAPAR for Maccarese, Italy (**left**) and LAI, LCC, FVC, and CCC for Shunyi, China (**right**).

It can be seen from Figure 5, that the k-fold cross-validation (k = 10) provided in some circumstances worse results (higher RRMSE) than the validation with ground data (e.g., Figure 5a). This is particularly evident for LAI and CCC for the Maccarese site. Usually, since it was carried out with subset of the same dataset as used for the training of the MLRAs, the cross-validation error was lower than the ground validation error, but this was not the case here, as can be seen from Figure 5a–e. The different statistical distributions of training and ground validation data might be a possible reason. It appears that the range of variation of the training sets was much larger than that of the ground validation sets, with more extreme values. This could have possibly led to more unreliable

retrievals in the cross-validation when extreme values were sampled, which did not happen in the ground validation.

For the Shunyi test site, RRMSE values were generally higher than for the Maccarese (Figure 5), probably because of the larger differences in the statistical distributions of the ground validation dataset from the training set, as compared for the Maccarese and possibly for the larger error introduced by the sampling protocol.

In general, for LAI, algorithms such as LSLR and PLSR, seemed to outperform more complex MLRAs such as NNs, for example providing better results than those found by applying the biophysical processor implemented in the ESA SNAP toolbox.

For leaf chlorophyll (LCC), much worse results were obtained for Shunyi than for Maccarese, with a clear advantage of RFTB for both sites for ground validation tests. The performance of the NN implementation of SNAP was particularly poor for Shunyi, though it should be noted that this variable was back-calculated from the CCC variable generated by SNAP.

Fractional ground cover (FVC) showed the smallest error among all the biophysical variables tested, alongside FAPAR which was only available for Maccarese. Neural networks provided the best accuracies with the lowest values of relative RMSEs for the latter variable.

Although different algorithms from different families of MLRAs performed well for different biophysical variables, the GPR was considered particularly interesting, as it provided uncertainty estimates with the mean value of prediction (Figure 3c). However, despite its good cross-validation results, GPR ranked as the top algorithm only for FVC and FAPAR for the Maccarese test site (Table 4). The AL procedure only selected the most informative samples of the dataset, thus the training of GPR was performed with a smaller set of PROSAIL simulations compared to the other algorithms, showing that it was quite efficient and robust.

The time taken for training the models ranged from 0.003 seconds for LSLR to 47.5 seconds for NNs and a maximum 548.3 seconds for GPR. Although, GPR took longer to train, it should be noted that it performed AL using EBD and used fewer samples to train and achieved high accuracy in comparison to other models.

When considering single ground sampling dates separately, generally the retrieval performances were improved (Table 6) compared to the bulked data (Table 4). In some cases, e.g., LAI and CCC, the worst results were obtained at the latest date, due to the fact that it was generally more difficult to estimate higher LAI values.

**Table 6.** Ground validation estimation results at each ground sampling date, for the best performing algorithms for the Maccarese test site.

#### **4. Discussion**

There are different methods for the retrieval of biophysical variables from remote sensing data and all the methods have their pros and cons [16]. Currently, hybrid methods have the capability of combining physical and statistical methods and are considered state of the art. In this paper, different MLRAs (kernel-based and non-kernel based) were trained and applied for the retrieval of the crop biophysical variables LAI, LCC, FAPAR, FVC, and CCC with the same configuration and settings of RTM PROSAIL simulated spectra as those implemented in ESA SNAP biophysical variables retrieval toolbox, with the exception of the observation geometry. With the same of configuration of PROSAIL parameters (Section 2.1), a total of 41,472 simulated spectra was generated, which turned out to be rather redundant and inefficient for performing the training steps. Subsets of 2500 randomly extracted simulated spectra were used to perform the training of the algorithms and this procedure was repeated ten times with different subsets. This was done in order to allow a comparison of alternative MLRAs to a well-established operational algorithm [29].

The SNAP algorithm relies on the use of PROSAIL simulations for training NNs, which are the most widely-used tools for operational biophysical variables retrieval [25–28]. The downside of the NNs are that they require a relatively long time for training, tuning of parameters is a difficult task, they are black box in nature, and can be unpredictable if training and validation data deviate from each other even slightly [5]. In this paper, various alternative MLRAs outperformed NNs. For example, LAI retrieval was best performed with LSLR and PLSR, consistently for both tests, although with different accuracies, i.e., RMSE of 0.68 for Maccarese and 1 for Shunyi [5], compared with different retrieval strategies for LAI and reporting the best performing algorithms for each category. In the case of parametric regression, the best performing algorithm was Tian 3-band formulation (RMSE = 0.615 and *R*2 = 0.823), whereas VH-GPR performed best among non-parametric regression algorithms (RMSE = 0.436 and *R*2 = 0.902). However, it should be noted that the accuracies found by these authors, somehow better than those of the present study, could be explained by the fact that they used cross-validation in which the training was carried out using the ground dataset, not independent model simulations such as in the present work.

The LCC was best retrieved by the RFTB method for both sites, revealing a general overestimation, but the error of estimation was higher for the Shunyi test site compared to the Maccarese (Figures 3b and 4b), particularly for SNAP. Also, the error of CCC estimates was very high for the SNAP tool compared to the other algorithms tested in this work for the Shunyi ground validation (Figure 5i). It was previously shown by Reference [64] that Sentinel-2 bands at a 10 m spatial resolution are suitable for estimating LAI, LCC, and CCC. They retrieved LAI (*R*2 = 0.809), LCC (*R*2 = 0.696) and CCC (*R*2 = 0.818) with vegetation indices approach. Multiple vegetation indices were compared for identification of potential vegetation index for the retrieval of LAI, LCC, and CCC in [65], the best correlation for LCC was with the Meris Terrestrial Chlorophyll Index (MTCI) (*R*2 = 0.77) and Sentinel-2 red-edge position index (S2REP) (*R*2 = 0.91), for LAI, inverted red-edge chlorophyll index (IRECI) (*R*2 = 0.77) and NDI45 (*R*2 = 0.62), Normalised Difference Vegetation Index (NDVI) for CCC (*R*2 = 0.70). These results are better than those found in the present work, but again, empirical approaches employing the ground datasets were employed by these authors for the calibration of the models.

The FVC variable was more accurately predicted than LAI or LCC, though with slightly higher error for Shunyi. With hierarchical tree-based methods such as BagT, RFTB, and GPR, the error reached was even lower than 10%, which is in line with Global Monitoring for Environment and Security (GMES) goal accuracy [66] (Figure 5c). In the case of CCC retrieval, it can be noted that the lowest RRMSE was provided by RFTB and BagT, though also GPR performs similarly to RFTB. The FAPAR variable was only estimated for the Maccarese test site, due to the unavailability of the ground measurements for validation, it was not estimated for the Shunyi test site. The results showed that this variable was best retrieved with the RFTB method, with the lowest RMSE (Table 3).

The GPR was proved to be an efficient and powerful regressor for the biophysical variables retrieval in previous reports [18]. A study conducted by Reference [18], retrieved LAI, Chl, and FVC using different methods, such as NNs, kernel ridge regression (KRR), support vector regression, and GPR for different Sentinel-2 and Sentinel-3 configurations. They found that an overall good performance throughout all Sentinel configurations was provided by GPR. A study carried out by Reference [10] also found that GPR using AL techniques is an efficient and robust method for the retrieval of LAI and LCC from Sentinel-3 OLCI spectra. However, this was not the case in our study. The GPR method, coupled with AL procedure, provided the best results in terms of computational time and low error only for FVC and FAPAR ground validation tests (Figure 5c,e,h). On the other hand, GPR generally performed as the best algorithm for the retrieval of all the biophysical variables, only for the cross-validations tests (Tables 4 and 5). Hence, of importance for further analysis to investigate how accurately GPR performs when validated against independent ground data. As Reference [10]

pointed out, GPR is based on non-parametric regression in a Bayesian framework, it provides insights in bands carrying relevant information and also in theoretical uncertainty estimates, thus partially overcoming the black box problem [5]. These uncertainties are a useful tool for the assessment of upscaling capabilities of biophysical variables from airborne or spaceborne platforms and their respective scales [32]. A study conducted by [9] introduced the AL approach with GPR and SVM regressions to deal with the problems of training sample collection for biophysical variables estimation. Their results obtained on simulated MEdium Resolution Imaging Spectrometer (MERIS) and real SeaWiFS Bio-optical Algorithm Mini-Workshop (SeaBAM) datasets were characterized by higher performances in terms of both accuracy and stability with respect to a completely random selection strategy. The present work seems to support these results highlighting the efficiency of the AL procedure, since, when using a smaller set of selected training data, comparable results were obtained than with a random selection of larger size (Figure 2).

#### **5. Conclusions**

Hybrid methods of biophysical variables retrieval rely on the generation of simulated spectra using physically-based RTM for the training of MLRAs, under the assumption of a more general applicability as compared to the training carried out using measured data, since RTMs allow to simulate a wider range of leaf and canopy properties. Operationally, so far, only NNs have been implemented as a hybrid method, e.g., in the ESA SNAP toolbox. Due to complexity of parameter tuning and the black box nature of the algorithm, in recent years, several studies have focused on alternative MLRAs such as kernel-based methods. They offer an interesting alternative to NNs in terms of performance and computational demand, sometimes partially overcoming the black box nature of NNs. In the present study, it was shown, using the same set of PROSAIL simulations, that some MLRAs could provide better results than NN-based algorithms implemented in SNAP. This was particularly evident for chlorophyll and fractional vegetation cover. It resulted that the best performing algorithms varied according to the biophysical variables to be estimated: e.g., LSLR and PLSR worked better for LAI, whereas RFTB and BagT for LCC. This study also showed that these variables' results were consistent for two different datasets respectively collected in Italy and China on wheat, though the retrieval accuracies depended somehow on errors introduced during the ground sampling protocol. The absolute values of the retrieval accuracies are in some cases higher than the best performances reported in previous studies, but most of these previous studies did not train the regression algorithms with completely independent datasets, as was done in the present study. Further studies should be conducted to investigate alternative strategies, e.g., of active learning procedures, using different kernel functions, and validate the retrievals using a wider range of agricultural crops, also assessing their performance for canopies that depart more from the assumptions of the turbid medium used in the PROSAIL RTM.

**Author Contributions:** D.U. implemented the methods, collected and analyzed the data for the Italian test site, and wrote the manuscript; R.C. designed the field experiments and collected the data in Italy, developed the idea and methodology, supervised its implementation and wrote the manuscript; S.P. (Simone Pascucci) and S.P. (Stefano Pignatti) collected and analyzed the data in Italy and edited the manuscript. W.H., X.Z., H.Y., and W.K. collected and analyzed the ground data for the China test site and edited the manuscript. All authors read and approved the final manuscript.

**Funding:** This research was funded by the European Space Agency (ESA) and the Ministry of Science and Technology of China (MOST), under the Dragon 4 Programme, Project ID32275-Topic1; it was also supported by the National Natural Science Foundation of China (41871339)and by the Italian National Research Council within the CNR-CAS Bilateral Agreement 2017–2019.

**Acknowledgments:** The authors gratefully acknowledge MIUR (Italian Ministry for Education, University and Research) for financial support (Law 232/216, Department of excellence) to the DAFNE Department. The authors would like to thank Marie Weiss (INRA Avignon, France) for providing the Matlab code of the BV-Net tool and the information on the SNAP biophysical processor. The authors would also like to thank Davide Tahani and Fabiola Fontana for helping in data collection and the Maccarese farm for their support and for granting access to the fields, and Jochem Verrelst for his support with the ARTMO toolbox.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **References**

© 2019 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).

## *Article* **Physically-Based Retrieval of Canopy Equivalent Water Thickness Using Hyperspectral Data**

#### **Matthias Wocher \*, Katja Berger, Martin Danner, Wolfram Mauser and Tobias Hank**

Department of Geography, Ludwig-Maximilians-University Munich, Luisenstraße 37, 80333 Munich, Germany<EMAIL_ADDRESS>(K.B<EMAIL_ADDRESS>(M.D<EMAIL_ADDRESS>(W.M<EMAIL_ADDRESS>(T.H.)

**\*** Correspondence<EMAIL_ADDRESS>Tel.: +49-892-180-6695

Received: 19 October 2018; Accepted: 27 November 2018; Published: 30 November 2018

**Abstract:** Quantitative equivalent water thickness on canopy level (EWTcanopy) is an important land surface variable and retrieving EWTcanopy from remote sensing has been targeted by many studies. However, the effect of radiative penetration into the canopy has not been fully understood. Therefore, in this study the Beer-Lambert law is applied to inversely determine water content information in the 930 to 1060 nm range of canopy reflectance from measured winter wheat and corn spectra collected in 2015, 2017, and 2018. The spectral model was calibrated using a look-up-table (LUT) of 50,000 PROSPECT spectra. Internal model validation was performed using two leaf optical properties datasets (LOPEX93 and ANGERS). Destructive in-situ measurements of water content were collected separately for leaves, stalks, and fruits. Correlation between measured and modelled water content was most promising for leaves and ears in case of wheat, reaching coefficients of determination (R2) up to 0.72 and relative RMSE (rRMSE) of 26% and in case of corn for the leaf fraction only (R2 = 0.86, rRMSE = 23%). These findings indicate that, depending on the crop type and its structure, different parts of the canopy are observed by optical sensors. The results from the Munich-North-Isar test sites indicated that plant compartment specific EWTcanopy allows us to deduce more information about the physical meaning of model results than from equivalent water thickness on leaf level (EWT) which is upscaled to canopy water content (CWC) by multiplication of the leaf area index (LAI). Therefore, it is suggested to collect EWTcanopy data and corresponding reflectance for different crop types over the entire growing cycle. Nevertheless, the calibrated model proved to be transferable in time and space and thus can be applied for fast and effective retrieval of EWTcanopy in the scope of future hyperspectral satellite missions.

**Keywords:** hyperspectral; spectroscopy; equivalent water thickness; canopy water content; agriculture; EnMAP

#### **1. Introduction**

The quantification of water stored in agricultural plants plays an essential role in understanding the impact of cultivated areas on the earth's water cycle. Due to its close association to biochemical factors, such as vegetation transpiration [1] and net primary production [2], the knowledge of quantities of water contained within agricultural crops is crucial, particularly for the development of environmental process models [3,4]. Moreover, quantifying canopy water content is important in regards to the water use efficiency of plants [5], evaluation of plant physiological status and health [6,7], and crop ripening monitoring [8].

Within the optical spectral domain (400 nm–2500 nm), absorption by vegetation liquid water occurs in the near-infrared (NIR) at 970 nm and 1200 nm and in the shortwave infrared (SWIR) at 1450 nm and 1950 nm [9,10]. Due to a higher absorption coefficient in the SWIR [11] most of the early studies combined those wavelengths with water insensitive wavelengths in the NIR to define empirical narrow-band indices for water content retrieval [12–15]. However, the strong absorption by water may saturate those bands at high water contents in optically thick canopies [16]. Moreover, absorption by atmospheric water vapor at 1450 nm and 1900 nm results in noisy measurements which renders these spectral regions unsuitable for further analysis [9], both for top-of-atmosphere (TOA) and top-of-canopy (TOC) spectroscopy. Some vegetation biophysical variables may disturb the signal of water: for instance Jacquemoud, et al. [17] noted that the leaf area index (LAI) masks the water signal between 1000 nm and 1400 nm and advised caution when using such indices for water retrieval. The comparatively low 970 nm water absorption depth is embedded in an area of generally high vegetation reflectance in the NIR. Due to low absorption it is expected that radiation at 970 nm penetrates deeper into the canopy reflecting a larger portion of its total water content without a tendency to saturation [18–22]. Therefore, Peñuelas, et al. [6,7] developed the 970 nm water index (WI) to retrieve relative plant water concentration (PWC). In the following, other studies also focused on the 970 nm absorption to estimate canopy water content [3,5,23,24].

Methodologically, the definition of a narrow-band spectral index to retrieve vegetation water content information constitutes the parametric regression type of methods. Their simplicity and thus computational feasibility make them highly desirable for large-scale remote sensing applications. However, a fundamental problem of parametric regression methods is their lack of generality and transferability [25]. Since indices are not solely influenced by liquid water, but also affected by leaf internal structure and leaf dry matter [26] or canopy structure, LAI and soil background [15,27,28], the established regression-based relationships and estimated quantities of water stored in a canopy are limited to local conditions [29]. Accordingly, the obtained results are site-, time- and crop-specific [30]. Moreover, as more hyperspectral image data with a continuous spectral coverage become accessible, the limited use of a small number of bands does not correspond to the up to date possibilities in view of the available data information density.

For the implications given, physically based model inversion methods have been introduced as a promising alternative to retrieve biochemical and biophysical vegetation variables. Radiative transfer models (RTM) describe interactions between solar radiation and vegetation constituents using physical laws. Their ability to generate an infinite number of simulated spectra with known input parameters conversely allows their inversion in order to estimate the underlying parameters. For the inversion of RTMs, a variety of strategies have been applied. These include numerical optimization algorithms, look-up table (LUT) approaches, artificial neural networks (ANN) and other machine learning algorithms (for an overview please refer to Verrelst, et al. [25,31]). Although RTM-based inversion methods are considered to be physically sound, the techniques require profound knowledge, are often computationally demanding and are mathematically highly non-linear [31,32]. Another limitation of RTM-inversion is the ill-posed nature or equifinality of model inversion. Many different parameter sets may be equally valid in terms of their ability to reproduce a measured reflectance spectrum (for a discussion of this topic see Atzberger and Richter [33]).

In view of future hyperspectral satellite missions like Italian PRISMA [34], US HyspIRI [35], Israeli-Italian SHALOM [36], European CHIME [37], and German EnMAP [38] fast and efficient retrieval methods for large datasets are required. Mathematically simpler physically-based approaches have been applied before to circumvent the equifinality problem and to reduce the computational effort of RTM-based model inversion. Green, et al. [39,40] originally incorporated the Beer-Lambert law to separate liquid water from atmospheric water vapor and determine both to allow the retrieval of surface reflectance from measured AVIRIS radiance. Thereby, the Beer-Lambert law was applied to directly infer the path length through optically active liquid water, i.e., the equivalent water thickness (EWT), from a measured reflectance spectrum using water absorption coefficients for pure liquid water [41,42]. Since multiple NIR scattering, and the attendant increase in optical path length at both the leaf and the canopy scale are not accounted for in the simple Beer-Lambert law [43,44], absolute quantification of EWT can only be achieved by calibration. Subsequently, validation has to be performed on accurate in-situ measurements. Studies that aimed at separating all three phases of water were not designed to quantify canopy water content in absolute terms and therefore accurate measurements were not carried out [45–48]. On the other hand, studies which derived water content explicitly by applying the Beer-Lambert law often relied on the assumption that upscaling leaf EWT

to canopy water content (CWC) could be done by a simple multiplication with the leaf area index (LAI) (see references [5,23,25,49,50]). In other publications, biomass sampling strategies have not been designed to deduce the single water components of a canopy that an optical sensor can actually detect (e.g., references [21,51–53]). Consequently, validation of these approaches could not approve translation into transferable and generally applicable retrieval tools [42].

Therefore, the objective of the present study was to test the performance of the Beer-Lambert law to retrieve crop water content from spectra with a contiguous spectral coverage around 970 nm and perform validation separately for leaves, stalks, and fruits by means of the two very different crop types: corn and winter wheat.

#### **2. Materials**

#### *2.1. Munich-North-Isar Test Site*

#### 2.1.1. Biomass Sampling and Water Content Determination

Biomass collection was performed in 2015, 2017, and 2018 at three winter wheat fields (*triticum aestivum*) and two corn fields (*zea mays*) of communal farmland 30 km north of Munich (southern Germany) east of the river Isar (Table 1, Figure 1).

**Figure 1.** Munich-North-Isar test sites overview (**left**) and exemplary 2017 simulated 30 × 30 m corn EnMap-pixel with 9 ESUs (**right**).

**Table 1.** Munich-North-Isar winter wheat and corn test sites, locations, periods of sample collection, number of biomass samples, and number of spectral measurements at cloud-free days.

Within the fields, three different sampling points were selected based on long-term biomass distribution pattern observations (TalkingFields Base Map: www.talkingfields.de) representing low, medium, and high persistent relative fertility.

Plant leaf water content is commonly expressed as equivalent water thickness (EWT, Equation (1)) corresponding to a hypothetical thickness of a single layer of water averaged over the whole leaf area [10]:

$$\text{EWT} = \frac{\text{FW} - \text{DW}}{\text{A}} \left[ \text{g } \text{cm}^{-2} \right] \text{ or } [\text{cm}]\_{\prime} \quad \left[ \text{kg } \text{m}^{-2} \right] \text{ or } [\text{mm}] \tag{1}$$

where FW is the fresh sample weight, DW is the oven dry weight and A is the leaf area. While EWT refers to the water content on leaf level, canopy water content (CWC, Equation (2)) is commonly derived through extrapolation by means of the LAI:

$$\text{CWC} = \text{EWT} \ast \text{LAI} \tag{2}$$

Due to the linkage of LAI to the whole canopy, CWC may be biased towards the leaf fraction. Furthermore, CWC does not allow inferring the actual water detectability of plant components in a canopy from total detected water. Consequently, in this study total EWTcanopy (EWTleaf + EWTstalk + EWTfruit) will be defined as the above-ground total equivalent water layer averaged over one square meter of ground surface (Equation (3)).

$$\text{total EWT}\_{\text{canopy}} = \sum (\text{FW}\_{\text{leaves}+\text{stalks}+\text{fruits}} - \text{DW}\_{\text{leaves}+\text{stalks}+\text{fruits}}) \ast \text{A}\_{\text{g}}^{-1} \tag{3}$$

where Ag denotes the ground area. To monitor the development of total amounts of water stored in the canopy throughout the growing season, biomass samples were collected on a weekly basis. In case of wheat, a minimum transect of 50 cm along a sowing track or an area of 0.25 m2 was cut at soil level. For corn, 2–3 plants were cut. In-field plant density was obtained by counting plants and rows per meter. The samples were separated into leaf, stalk and fruit compartments, weighed in fresh state and oven-dried until constant weight for 24 h at 105 ◦C before dry weight was determined. EWTleaf, EWTstalk, EWTfruit (EWTear and EWTcob, respectively) and total EWTcanopy per cm2 (Table 2) were calculated from laboratory results (specific water contents per ground area) and from farm management metadata (plants per meter and row spacing). The phenology was determined according to secondary growth stages of the BBCH-scale [54].

**Table 2.** Statistics (range, mean, standard deviation) for in-situ measured EWTleaf, EWTstalk, EWTfruit, total EWTcanopy and BBCH-range. Values correspond to measurements with available spectral reflectance data.

#### 2.1.2. Spectroscopic Measurements

At each study site, a 30 × 30 m grid of nine 10 × 10 m squares was marked out delineating the elementary sampling units (ESU). This grid layout was designed to trace the geometric properties of one EnMAP pixel; hence, regarding the viewing geometry and grid location, the sensors descending orbit and inclination angle of 97.96◦ was accounted for (Table 1). At each sampling date with clear sky conditions (Table 1) all nine ESUs were revisited and spectral measurements were taken using an Analytical Spectral Devices Inc. (ASD, Boulder, CO, USA) FieldSpec3 Jr. spectroradiometer with an effective spectral resolution of 3 nm in the VIS (≤700 nm) and 10 nm in the NIR and SWIR (≤2500 nm). Five nadir measurements were conducted per ESU at a height of 25 cm above the canopy, the same

height at which the white reference panel (OptoPolymer, Munich, Germany) could be fully observed with the instruments field of view of 25◦. Throughout the measurements, the sensor was slightly moved back and forth manually while maintaining the observation angle to obtain a representative spectral sample of the canopy. The five recorded spectra were averaged per ESU and a spatial mean of the nine ESUs was calculated to provide reflectance of the complete 30 × 30 m grid. Further, post-processing included splice-correction, white reference baseline calibration, and slight smoothing using a Savitzky-Golay-Filter with a frame size of 13 nm.

Note that it was not possible to conduct destructive sampling at exactly the same locations where the continuous spectral measurements were taken. However, due to averaging of spectral sampling points over the 30 × 30 m sampling area, it was possible to capture the in-field variability and therefore to represent average field water conditions.

Altogether, the collected dataset comprised destructively measured, plant compartment specific water content samples with corresponding spectral measurements at 26 dates for wheat and 14 dates for corn over three and two years, respectively (see Tables 1 and 2).

#### *2.2. Leaf Optical Data*

Preliminary tests of the EWTcanopy retrieval model presented in this study were performed on two different leaf optical datasets. The LOPEX93 database was established in 1993 by the Joint Research Centre (JRC, Ispra, Italy). It associates transmittance and reflectance in the range of 400–2500 nm with biophysical and biochemical measurements of 66 leaf samples from 45 species [55]. In total, the database comprises 330 spectra with corresponding measurements of EWT. Secondly, tests were performed on the ANGERS database containing 276 reflectance spectra and EWT measurements of 43 species [56]. While woody species make up the majority of the ANGERS database, both datasets represent a large variety of leaf internal structure and spectra.

#### *2.3. Radiative Transfer Models and Look-Up Tables*

To check consistency between leaf optical data and modelled spectra, large look-up tables (LUT) using the RTMs PROSPECT and PROSAIL were created. PROSAIL [17] is coupling the *Leaf Optical Properties Spectra* model PROSPECT [57] and the turbid medium canopy reflectance model 4SAIL (*Scattering by Arbitrary Inclined Leaves*) [58,59]. The latest recalibrated version PROSPECT-D [60,61] simulates bidirectional-hemispherical reflectance and transmittance in the optical domain as a function of leaf pigments (chlorophyll a+b content Cab, carotenoids Car, and anthocyanins Canth), dry matter Cm, and brown pigments Cbrown as well as a leaf mesophyll structure parameter N, and EWT(Cw). The canopy model SAIL calculates a bidirectional reflectance factor of 1-D turbid medium plant canopies. With regard to leaf optical properties and reflectance of the underlying soil (psoil), it implements canopy structure (LAI), average leaf inclination angle (ALA) or optionally, ellipsoidal leaf inclination distribution (LIDF), and hot spot size parameter (hspot) for a given illumination and viewing geometry (observation zenith angle (OZA), relative azimuth angle (rAA) between sun and sensor, and the solar zenith angle (SZA)).

Considering the impact of plant foliar water on the 970 nm absorption band, during LUT generation, all parameters with sensitivity in the NIR region were uniformly distributed over a wide value range (Table 3). Leaf pigments, having no effect on reflectance in the NIR [26,44], remained constant.

**Table 3.** Parameter ranges for PROSPECT-D and PROSPECT-D + 4SAIL (PROSAIL) LUT. Specified ranges are uniformly distributed, single values are fixed.

#### **3. Methods**

*3.1. The Beer-Lambert Law and Retrieval Method Development*

The Beer-Lambert law is mathematically formulated as Equation (4):

$$
\Phi = \Phi\_0 e^{-\mathfrak{a}(\lambda)d}.\tag{4}
$$