In-Situ and Aircraft Reflectance Measurement Effectiveness for CAL/VAL Activities: A Study over Railroad Valley

Abstract

This paper aims to assess the relationship between the surface reflectance derived from ground based and aircraft measurements. The parameters of the Rahman–Pinty–Verstraete (RPV) and Ross Thick-LiSparse (RTLS) kernel based bi-directional reflectance distribution functions (BRDF), have been derived using actual measurements of the hemispherical-directional reflectance factor (HDRF), collected during different campaigns over the Railroad Valley Playa. The effect of the atmosphere, including that of the diffuse radiation on bi-directional reflectance factor (BRF) parameter retrievals, assessed using 6S model simulations, was negligible for the low turbidity conditions of the site under investigation (${\tau }_{550}\le 0.05$). It was also shown that the effects of the diffuse radiation on RPV spectral parameters retrieval is linear for the isotropic parameter ${\rho }_0$ and the scattering parameter $\mathrm{\Theta }$, and can be described with a second order polynomial for the k-Minnaert parameter. In order to overcome the lack of temporal collocations between aircraft and in-situ measurements, Monte Carlo 3-D radiative transfer simulations mimicking in-situ and remote sensing techniques were performed on a synthetic parametric meshed scene defined by merging Landsat and Multianglhe Imaging Spectroradiometer (MISR) remote sensing reflectance data. We simulated directional reflectance measurements made at different heights for PARABOLA and CAR, and analyzed them according to practices adopted for real measurements, consisting of the inversion of BRF functions and the calculation of the bi-hemispherical reflectance (BHR). The difference of retrievals against the known benchmarks of kernel parameters and BHR is presented. We associated an uncertainty of up to 2% with the retrieval of area averaged BHR, independently of flight altitudes and the BRF model used for the inversion. As expected, the local nature of PARABOLA data is revealed by the difference of the anisotropic kernel parameters with the corresponding parameters retrieved from aircraft loops. The uncertainty of the resultant BHR fell within $\pm 3\%$.

1. Introduction

Surface albedo is a crucial parameter in order to understand the behavior of the earth climate system, because it regulates the fraction of solar energy reflected back to space by the surface-atmosphere coupled system [1,2,3]. Myhre et al. [4] reported a wide range of estimates of albedo change due to anthropogenic land use change and its negative radiative forcing related to the higher albedo of bare soil/crops with respect to forests, or positive feedback in polar areas where snow melting produces an opposite effect. As such, the albedo was defined as an essential climate variable (ECV) by the Global Climate Observing System (GCOS) of the World Meteorological Organization (WMO) [3,5]. The GCOS implementation plan (GCOS-200 [5]) has specified a set of global target requirements that in many cases may meet local and regional needs, with a spatial resolution of 200/500 m, resolved on a daily basis, with an accuracy lower than the maximum value of 5% or 0.0025, and a stability error lower than the maximum value of 1% and 0.001.

Satellite Earth observation provides the best temporal and spatial coverage to monitor the evolution of the lithosphere, cryosphere and biosphere [6]. To produce fit-for-purpose outputs, accurate instrument calibration and algorithm assessment are necessary. Adequate knowledge of the limits of retrieval algorithms, throughout the process that transforms the measured physical quantities (typically spectral radiances) into derived land or atmospheric products, needs to be guaranteed. A common way to perform calibration and validation activity is based on specific and intensive field measurement, during which in-situ data are collected to define the surface status. Simultaneously, atmospheric optical properties have to be defined to account for the atmospheric contamination of the signal measured by the satellite [7].

Locations characterized by high signal stability, high signal-to-noise ratio, low impact of the atmosphere and possibly weak and stable surface reflectance anisotropy, such as a desert or dry lake bed, are typically chosen as calibration sites [6,8]. They are employed to assess sensor calibration trends that can affect long term missions by introducing biases and inconsistencies in long term climate data records [9,10].

While assessing the quality of retrieval algorithms, the surface complexity should be as low as possible to minimize problems related directly to its heterogeneity, or adjacency effects, and avoiding issues related to data manipulation such as re-gridding and re-projection. Ground targets should be chosen to be the most homogeneous possible, with stable optical properties, and easily accessible to host intensive measurement campaigns.

In situ radiometric measurements should be performed routinely to assess the variability of the surface reflectance over time and its dependence on the solar position, atmospheric turbidity, and cloud coverage. These measurements should be taken following standardized and traceable methods, whenever available, and adequate spatial sampling approaches. Hence, it can be argued that field campaigns require substantial planning efforts, and a large amount of human and logistic resources [6].

Koukal et al. [11] inverted the Rahman–Pinty–Verstraete (RPV [12]) and the Ross Thick and Li-Sparse (RTLS [13]) BRDF models over digital hemispherical camera reflectances with similar retrieval performances. The study indicated a minimum of ten viewing directions needed to exploit a forest classification task. Hill et al. [14] analyzed the relationship between RPV and RTLS parameters obtained respectively from Multiangle Image Spectro Radiometer (MISR) and MODerate resolution Imaging Spectroradiometer (MODIS) data collected over an Australian tropical savanna, showing that the asymmetry parameter in particular is suitable to differentiate between different canopy structures.

Chen et al. [15] compared the MISR Bi-hemispherical Reflectance (BHR) and directional hemispherical reflectance (DHR) broadband products, showing good agreements for land and overestimation of the spectral BHR for snow/ice. The root mean squared error (RSME) and the correlation coefficient ($r^2$) metrics were used to assess the performances of MISR BHR with respect to ground based measurements of the albedo, highlighting the need for the scaling-up of point measurements.

Román et al. [16] assessed the spatial representativeness of broadband albedo measured from towers over forested landscapes using a hybrid approach based on MODIS and ETM+ data. Cescatti et al. Cescatti et al. [17] investigated the representativeness of satellite albedo products using MODIS and FLuxNET data, discussing the importance of spatial representativeness of point data such as tower measurements, linking comparison performance with the landscape heterogeneity.

In order to be effectively used for calibration and validation of remote sensing products, in situ data should be scaled up to a common resolution [17,18,19]. In particular, Bruegge et al. [19] described the vicarious calibration [20] activities implemented for TOA radiances and performed for different platforms including MISR, MODIS, Landsat, AirMISR over the Railroad Valley Playa (RRV). More recently, Kharbouche et al. [21] investigated the CAR to MISR upscaling process.

The number of measurement campaigns allowing a direct comparison between in-situ and aircraft data is often limited, because for logistical reasons they are often not collocated in space and time, and they can be affected by variability in the surface or sky conditions, preventing a consistent comparison.

Here we propose to show how three-dimensional radiative transfer modeling (3D-RTM) can be employed to help understand the impacts on albedo retrievals, and the corresponding equivalent pixel size, produced by (i) the difference in height between in-situ and aircraft measurements and (ii) the number and sampling of values for parametric model inversion to validate the spectral surface albedo.

2. Methodology

This work aims to assess the consistency of the surface reflectance and albedo, as evaluated through ground based and aircraft measurements at different spatial resolutions over the Railroad Valley Playa vicarious calibration site (RRV, Nevada, USA, 38.50${}^{\circ }$N, 115.69${}^{\circ }$W, 1435 m ASL https://calval.cr.usgs.gov/apps/radiometry-image-gallery) [22,23].

We decided to describe surface optical properties in terms of bidirectional reflectance distribution function (BRDF) parametrization, which easily allows the calculation of the spectral albedo by means of hemispherical integrations [24,25].

The bidirectional reflectance factor (BRF) is a convenient way to express the anisotropic properties of the surface reflectance, and is defined by the ratio between the exiting radiance measured in a certain direction identified by the zenith and azimuth angles (${\theta }_v,{\varphi }_v$), and the one that would be reflected by an ideal lambertian surface in the same direction [24,25]. The hemispherical-directional reflectance factor (HDRF) differs from the BRF because it is measured in ambient conditions, where the incoming radiation is a combination of a direct and a diffuse component, the latter affecting the BRF measurement.

We adopted two different BRF models based on the combination of specific kernel based functions in the ${\sum }_{i=1}^3{p}_i{K}_i$ or ${\prod }_{i=1}^3{K}_i\left(p_i\right)$ forms, with $p_i$ representing the function parameters, and $k_i$ the corresponding kernels, which incorporate the angles describing the illumination and viewing conditions. They were the RTLS model [13] and the RPV model [12], respectively adopted by the Moderate Resolution Imaging Spectroradiometer (MODIS) albedo products’ team and the Advanced Very High Resolution Radiometer (AVHRR) team. Each kernel is meant to catch a particular feature of the reflectance anisotropy. Typically, an isotropic parameter is used in combination with a $K_{iso}=1$ constant kernel, to define a baseline Lambertian reflectance. Then, for RTLS, volumetric ($K_{vol}$) and geometrical ($K_{geo}$) kernels, combined with their multiplying parameters, describe the anisotropy as originating from the radiative transfer occurring within the media volume, and the geometrical features due to surface roughness or related to surface elements’ mutual shadowing. Then, the bi-directional reflectance is parametrically described through three numbers $f_{iso}$, $f_{vol}$ and $f_{geo}$, and the corresponding kernels, which depend on illumination and viewing angles only. Similarly, RPV describes the BRDF by means of an isotropic parameter (${\rho }_0$), a k-Minneart parameter responsible for catching the ’U’ or bell shaped form of BRDF with respect to the viewing zenith angle ${\theta }_v$, and an asymmetry parameter ($\mathrm{\Theta }$) which determines, through a Henyey–Greenstein kernel function, the re-partition of the reflected radiance in forward or backward semi-hemisphere, with respect to the solar position. A modified version of the RPV function, a different asymmetry kernel and parameter ($b_M$) formulation, are adopted to optimize the production of operational Multiangle Imaging Spectroradiometer (MISR) surface albedo.

This study is based on actual measurements, and is supported by simulations produced by 3-D Monte-Carlo radiative transfer modeling [26] to account for the lack of concurrent surface and aircraft based campaigns. The measurements are used principally to define the intrinsic uncertainties of the instruments as operated in the RRV area. The reference instruments adopted here are the Portable Apparatus for Rapid Acquisition of Bidirectional Observations of Land and Atmosphere (PARABOLA) [22,27,28], and the Cloud Absorption Radiometer (CAR) [21,29] for ground and aircraft based measurements, respectively. They both measure the radiance at different angles in the spectral bands listed in

Table 1. Central wavelengths (in nm) relevant to PARABOLA for Cloud Absorption Radiometer (CAR) (1999–2011), Multiangle Imaging Spectroradiometer (MISR), and Thematic Mapper (TM) instruments. The weighed spectral reflectances $S{R}^{\prime }$ in the red and NIR bands are given in parenthesis.

PARABOLA	CAR	MISR	TM
551	-	557.5	561
650 (0.423)	682 (0.438)	671.7 (0.431)	661 (0.429)
860 (0.466)	870 (0.466)	866.4 (0.466)	832 (0.465)

. Figure 1 shows the relative spectral responses of the red and near-infrared (NIR) bands of PARABOLA, CAR, Multiangle Imaging Spectroradiometer (MISR) and the Landsat Thematic Mapper (TM) instruments. Atmospherically corrected reflectances obtained from TM and MISR were used to define a virtual model for the 3-D simulations (Section 5). We computed the the spectrally weighted reflectances in red and NIR bands for each instrument listed in Table 1 adopting the convolution defined as follows: $S{R}^{\prime }={\int }_0^{\infty }{S}_{0\lambda }S{R}_{\lambda }R{R}_{\lambda }d\lambda /{\int }_0^{\infty }{S}_{0\lambda }R{R}_{\lambda }d\lambda$, where $R{R}_{\lambda }$ is the instruments’ relative spectral responses, $S{R}_{\lambda }$ is the surface spectral reflectance taken with a spectroradiometer over the RRV, adapted from [30] (Figure 1), and $S_{0\lambda }$ is the extraterrestrial solar spectrum [31]. For the particular conditions of RRV, SR’ varies between 0.423 (PARABOLA) and 0.438 (CAR) with a relative error of 3% in the red band. Due to the flat reflectance across the wide spectral range from 750 up to 1000 nm, the four instruments provide the same near-infrared $S{R}^{\prime }$ within a 0.2% relative error.

Whilst both are based on the mechanical orientation of the optical part, the sampling strategy of PARABOLA differs from that of CAR. The ground footprint (gFOV) is defined by the projection of the instrument field-of-view (FOV) on the ground surface, and its dimensions can be expressed in terms of the viewing zenith and azimuth angles (${\theta }_v$, ${\varphi }_v$) with respect to the local normal to the surface, the height of the instrument with respect to the surface H, and the FOV itself. In this work, ${\theta }_v$ is assumed to vary between 0${}^{\circ }$ in the nadir direction and 90${}^{\circ }$ at the horizon, while ${\varphi }_v$ is assumed to be 0${}^{\circ }$ at North and increases clockwise.

PARABOLA measures the signal reflected by the ground surface $D{N}_{surf}^{\uparrow }({\theta }_i,{\varphi }_i,{\theta }_v,{\varphi }_v)$ in different viewing direction and solar positions (${\theta }_i$, ${\varphi }_i$), with an angular resolution of 5${}^{\circ }$, from a nominal height of 3 m, with a nominal optical FOV of 5${}^{\circ }$. Viewing angles vary from 0${}^{\circ }$ to 85${}^{\circ }$ in zenith, and from 0${}^{\circ }$ to 355${}^{\circ }$ in azimuth. The area covered by the instrument depends on the height of installation and the footprint varies with the viewing zenith angle. It typically performs a complete scan in about 3 min, a temporal window that guarantees sufficient stability of solar angles and incoming radiation.

The reference used for the BRF calculation is provided by the signal reflected by a Lambertian Spectralon panel with a flat spectral 100% reflectance (e.g., ideal reflectance) positioned in the nadir viewing sector ($D{N}_{lamb}^{\uparrow }$). The BRF is obtained according to its definition [24] as,

$$BRF({\theta }_i,{\varphi }_i,{\theta }_v,{\varphi }_v)=\frac{D{N}_{surf}^{\uparrow }({\theta }_i,{\varphi }_i,{\theta }_v,{\varphi }_v)}{D{N}_{lamb}^{\uparrow }}.$$

(1)

The mechanical rotation of CAR is limited to ${\theta }_v$, and it can sample from nadir to zenith (0${}^{\circ }$ to 180${}^{\circ }$) to allow both ground and sky sampling, or be set to perform a complete across-tracking scan of the surface [32]. Its angular range includes $\pm 5^{\circ }$ offsets allowing the removal of the aircraft roll effect. The sampling resolution is 0.5${}^{\circ }$, FOV is 1${}^{\circ }$ and a complete individual scan is almost instantaneous as it is performed in a few seconds (100 scan lines/minute). The azimuth angle is determined by the aircraft heading, as the instrument optical axes are set perpendicularly with respect to the aircraft’s axes. A typical sampling approach to scan the whole azimuth range is to perform, in 2–3 min, circular paths of a few kilometers in diameter [29,33].

The CAR high sampling rate guarantees a continuous mapping that depends on the aircraft speed and altitude. The resolution of the ground footprint at nadir view ($d=tanFOV\sim FOV$) for a small $FOV$ is 3.5, 8, 11 and 17.5 m for flight altitudes of 200, 500, 645 and 1000 m, respectively.

If we consider an aircraft speed of ∼80 m/s, a complete scan line is obtained every ∼50 m. The acquisition of each of the 382 viewing positions is nominally performed in about 1.5 ms with the aircraft translating by a distance of 12–13 cm. This translation can be considered negligible with respect to the ground footprint dimension of CAR listed above, even at the lowest altitude of ∼200 m at nadir this is equal to ∼350 cm. At higher observation angles and altitudes, the footprint increases while the translation remains constant. Hence, the influences of the aircraft movement can be neglected, and the pixel can be considered a quasi-static image of the ground with no averaging effects. It is important to simulate the measurements in the most realistic way compatible with the implementation model limits.

With respect to PARABOLA, the CAR performs only an absolute measurement of radiance ($L^{\uparrow }$) and the $BH{R}_h$ is computed at the aircraft flight level h, with respect to the TOA spectral irradiance ($E_s^{\downarrow }$) as

$$BH{R}_h({\theta }_i,{\varphi }_i,{\theta }_v,{\varphi }_v,h)=\frac{\pi L^{\uparrow }({\theta }_i,{\varphi }_i,{\theta }_v,{\varphi }_v,h)}{E_s^{\downarrow }\left({\theta }_i\right)cos{\theta }_i}.$$

(2)

The quantity defined in Equation (2) is also known as apparent reflectance (${\rho }^{*}$). This formulation requires an atmospheric correction to obtain accurate information on the surface [32].

In both the equations above, the term BRF is used to describe a physical quantity that is only approximated by real measurements, as these are constrained by the finite value of solid angles, and the smoothing effect provided by the diffuse component of incident radiation.

Given this, another difference of the field measurements with respect to the theoretical formulation of the bidirectional reflectance factor, is that the incident flux consists of a direct ($E_D$) and a diffuse component ($E_d\left({\theta }_i\right)=2\int L_i^{\downarrow }cos{\theta }_{i}sin{\theta }_{i}d{\theta }_i$), which both contribute to the exiting radiance. Following [34], the contributions to the exiting radiance $L^{\uparrow }({\theta }_v,{\varphi }_v,\varphi )$ in a specific viewing direction is given by,

$$L^{\uparrow }({\theta }_i,{\theta }_v,\varphi )=\frac{1}{\pi }BRF({\theta }_i,{\theta }_v,\varphi )E_D\left({\theta }_i\right)cos{\theta }_i+2\int L_i^{\downarrow }({\theta }_i,{\varphi }_i)BRF({\theta }_i,{\theta }_v,\varphi )cos{\theta }_{i}sin{\theta }_{i}d{\theta }_i,$$

(3)

where the BRF is expressed in terms of an azimuth difference $\varphi ={\varphi }_v-{\varphi }_i$, as it is commonly accepted to adopt a cylindrical symmetry of the BRF for natural surfaces. In Equation (3), the first term describes the contribution of the direct solar radiation component to the reflected radiance, which decreases with atmospheric turbidity, while the second term describes the contribution of the diffuse component $L_i^{\downarrow }({\theta }_i,{\varphi }_i)$, which in contrast increases with turbidity.

The effect of the atmosphere was evaluated by performing simulations of reflected radiance in different directions for a known BRF and with increasing optical depths, followed by an optimization of the BRF functions so as to evaluate the kernel parameter dependence on turbidity. These results are presented in Section 3.

Remote sensing techniques to define surface reflectance and albedo are based on the optimization of a BRF function over a limited set of measurements. Accepting an uncertainty cost related to their capability to represent all the natural anisotropic features, this is an approach that has the advantage of being able to determine the directional reflectance in any other direction, as well as its integrals, normally expressed in terms of the DHR and BHR [25], whether purely direct or diffuse illumination are considered. Presented in Section 4 are the PARABOLA and CAR measurements performed over the RRV and the assessment of the performance of two independent BRDF functions, to stress their capability and robustness in representing the reflectance anisotropy at various sub-samples of the full set of measurements. The RPV and RTLS functions were fitted through a non-linear-least-square algorithm [35,36] to obtain the model parameters and the BHR for selected wavelengths of PARABOLA and CAR.

To support findings derived from field measurements which are not collocated in time, we created a virtual scene using a combination of high (decametric) and medium (hectometric) resolution satellite derived products. Then, we simulated the BRF for PARABOLA and CAR from different viewpoints and from different altitudes using the Raytran 3-D radiative transfer code [26,37]. The reflectance arises directly as a ratio between the reflected number of photons exiting the scene in the viewing direction ($N^{\uparrow }$), and those entering the scene ($N^{\downarrow }$), i.e., $BRF({\theta }_i,{\theta }_v,\varphi )=\pi N^{\uparrow }({\theta }_i,{\theta }_v,\varphi )/N^{\downarrow }\left({\theta }_i\right)$.

By excluding explicitly the atmosphere in the 3-D model approach, attention was focused on the scaling aspects of the problem and the sampling approach. The setup of the virtual scene, virtual measurement simulations, and the discussion of the results obtained for both PARABOLA and CAR, are described in Section 5.

3. Effects of the Ambient Diffuse Radiation on BRDF Parametrization

In this section we discuss the atmospheric effects on the surface bidirectional reflectance factor (BRF) as evaluated in ambient light conditions. The reflected radiance at the surface and aircraft levels were calculated using simulations performed with the Second Simulation of the Satellite Signal in the Solar Spectrum (6S), a code based on the successive orders of the scattering method to solve the radiative transfer equation through the atmosphere in the solar spectrum, allowing for a detailed description of the surface anisotropy and aerosol optical properties [11,38,39].

As already mentioned, in situ measurements performed with either the PARABOLA or CAR represent a hemispherical-directional reflectance factor (HDRF, [24]), because reflectance is determined under ambient illumination, which includes a diffuse component that varies with the atmospheric turbidity. Martonchik [40] defined an accurate iterative procedure to obtain BRF values from HDRF measurements. Rahman et al. [12] and Tanre et al. [41] considered that the reflectance at the surface under ambient illumination is a weighted average between the BRF and the spherical reflectance (bi-hemispherical reflectance) due to the diffuse radiation effect.

In this work we provide an estimation of the error through an empirical forward modeling approach. We simulated the reflectance in the principal plane for different parameterizations of the RPV function [12] and increasing turbidity using the 6S code. Then, we performed the inversion of the RPV function on simulated reflectances and compared the optimized parameters with their expected values.

Figure 2 shows the reflectance at 550 and 870 nm in the principal plane, for two solar zenith angles and for various atmospheric turbidities as calculated at the ground level adopting a formulation similar to that presented by Lewis and Barnsley [42] for the blue-sky albedo, expressed in terms of a weighted sum of black-sky and white-sky albedo. By removing the hemispherical integration across viewing angles, we obtained a similar expression for the HDRF as

$$HDRF({\theta }_i,{\theta }_v,\varphi )\approx (1-d)·BRF({\theta }_i,{\theta }_v,\varphi )+d·\frac{1}{\pi }{\int }_0^{\pi /2}{\int }_0^{2\pi }BRF({\theta }_i^{\prime },{\theta }_v,{\varphi }^{\prime })cos{\theta }_i^{\prime }sin{\theta }_i^{\prime }d{\theta }_i^{\prime }d{\varphi }^{\prime }$$

(4)

where $d=E_d/(E_d+E_D)$ is the spectral diffuse ratio at the instrument level. Equation (4) shows that for $d=0$, the HDRF equals the BRF.

In addition to the pristine case in which the simulations represent the BRF, three examples are provided with increasing ${\tau }_{550}$, assuming values of 0.05 for an extremely clear air condition (visibility v∼ 150 km), compatible with the characteristics of the Railroad Valley Playa [23], of 0.3 for intermediate conditions (v∼ 15 km), and of 0.7 for a turbid condition (v∼ 5 km). The US62 standard atmosphere and a continental aerosol model are assumed. The HDRF defined in Equation (4) has been modeled for ${\theta }_i$ of ${20}^{\circ }$ and ${50}^{\circ }$ in the principal plane, using the 6S implementation of the RPV model. The RPV parameters were set to (${\rho }_0$ = 0.17, k = 0.75, $\mathrm{\Theta }$ = −0.12) for 0.55 $\mathsf{\mu }$m and (0.23, 0.76, −0.09) in the NIR. They were obtained from PARABOLA data acquired on 24 June 2009.

As shown in Figure 2, the parameters are representative of a pronounced backward reflection and an appreciable hot-spot, occurring in the principal plane for ${\theta }_v={\theta }_i$, which is more evident at higher sun zenith angles. Increasing turbidity smooths the reflectance simulations with respect to the pristine conditions, in particular around the hot-spot and at higher viewing zenith angles. For nadir viewing, the difference with respect to the black-sky conditions (labeled as noatm in the figure) appears almost negligible. For ${\tau }_{550}=0.05$ and ${\theta }_i={20}^{\circ }$ the differences in the HDRF with respect to the BRF (solid ines) in the principal plane are appreciable only in the green band. At ${\theta }_i={50}^{\circ }$ differences of up to 0.02 (over 0.5, hence ∼4%) are observed in the hot-spot for the green band.

These forward simulations were used as input to an optimization scheme of the RPV function, to assess the errors made on retrieved parameters by neglecting the diffuse radiation effects. The retrieved parameters are shown in Figure 3 as a function of the spectral diffuse ratio (d). For equal aerosol loads (identified by different symbols), the diffuse ratio is greater in the green channel due to both the Rayleigh (gas) scattering effects, and the lower scattering effects of continental aerosols at higher wavelengths. At ${\theta }_i={20}^{\circ }$, it amounts to 0.12 (green) and 0.04 (NIR) for ${\tau }_{550}=0.05$, and increases at higher sun angles (dashed lines) due to the greater air mass. We observed a linear increase of ${\rho }_0$ and $\mathrm{\Theta }$ estimates with turbidity, with a slope in the green band of $+0.06/d$ for ${\rho }_0$, and of $+0.1/d$ for $\mathrm{\Theta }$. Then, assuming a value of $d=0.12$ as a typical diffuse ratio in the green band over RRV, the overestimation of ${\rho }_0$ was found to be around 0.007 (+4% ), and promptly corrected using its linear dependence in d. The ${\rho }_0$ and $\mathrm{\Theta }$ overestimations do not depend considerably on ${\theta }_i$. The k parameter presents a non-linear dependence with respect to the diffuse ratio. The bias on k is of the opposite sign to the two ${\theta }_i$, but it is below 5% even for the most turbid conditions, and below 2% for the typical transparency of RRV.

How the atmosphere affects the BRF measurements performed on board an aircraft depends on the flight altitude and the atmospheric gas and aerosol profiles. 6S can simulate the spectral reflectance (or radiance) at any altitude between 0 and 100km within the atmosphere. The diffuse ratio is provided by default only at the ground level, and the apparent reflectance is normalized with respect to the extraterrestrial spectral radiation ($E_s$), as shown in Equation (2). The simulations of absolute calibrated radiances (with errors between ±1% and ±3% for all spectral channels) are normalized with respect to the TOA solar radiation [32] to obtain an apparent reflectance (${\rho }^{*}$) as the reference is relevant to a different height than the aircraft altitude.

We provide here only a summary of the results obtained for the same geometric conditions and atmospheric properties discussed above, and for altitudes relevant to a specific campaign performed over Railroad Valley (RRV) in 2008 (e.g., 180, 645, 1480 and 3400 m above ground). As can be expected, the atmospheric and environmental contributions to the reflected radiance increases with d and altitude. Figure 4 quantifies this effect for the green and near-infrared bands at different altitudes.

A flat attenuation of the pixel reflectance at all viewing angles is observed. In the green channel gases and Rayleigh scattering cause ${\rho }^{*}$ to decrease by ∼5% at all ${\theta }_v$ in the principal plane and 180 m above ground. In general, the apparent reflectance decreases by up to 10% in the green channel with increasing turbidity (Figure 4).

This effect is due to the atmospheric attenuation of the radiation reaching the surface, and then reflected to the sensor. In cases with a very low surface reflectance the atmosphere produces an increase of the apparent reflectance because of its scattering contribution to the signal.

For the aircraft flying height considered here (180 m) the signal reduction is not compensated by the intrinsic atmospheric reflectance ${\rho }_a$ and the apparent reflectance ${\rho }^{*}$ underestimates the surface reflectance. As the observation altitude h and ${\theta }_v$ increase, the contribution of the atmosphere becomes more evident, in particular in the forward hemisphere. For ${\tau }_{550}$ of 0.3, the apparent reflectance is equal to the surface reflectance at ${\theta }_v\sim {60}^{\circ }$ at 1480m, or ${\theta }_v\le {50}^{\circ }$ at 3400 m.

Additional investigations have shown that the relative contribution of the target with respect to the whole signal, including environment and atmosphere, is maximum in the hot-spot direction and decreases at higher viewing angles. In the green band, with ${\theta }_i={28}^{\circ }$ and ${\tau }_{550}$ = 0.3 (∼15 km visibility), two-thirds of the signal is given by the target for a ${\theta }_v=0^{\circ }$, but increases to 85% for a ${\tau }_{550}$ = 0.05, the typical condition of observation over RRV. This partitioning remains almost constant with the observing altitude h.

Contrarily to RPV, the RTLS model was not originally implemented in 6S. In the future, it would be worth performing similar analysis importing the RTLS BRDF functions into 6S, where we could expect similar behavior for the isotropic parameter $f_{iso}$, at least, due to very similar performances of RPV and RTLS in anisotropy representation.

4. In-Situ and Aircraft Measurements

4.1. Bidirectional Reflectance Data from In-Situ Campaigns (PARABOLA)

Surface based hemispherical-directional reflectance measurements were performed using PARABOLA over RRV during different field campaigns [33,43] and in different site locations shown by P1-P8 pins in Figure 5. We focused on the North-West (N-W) and the South-East (S-E) areas of the RRV, indicatively defined by the two main groups of aircraft loops described in Section 4.2.

Selected HDRFs, measured over P7 (S-E) taken on 2014-06-24, P2 (N-W) taken on 2014-06-25, and a dataset collected on 24 June 2009, were used in this part of the research. Even if unflagged, in terms of exact position, the latter dataset is the most complete in terms of daily coverage and permitted us to explore the full dependence on solar angles and to exploit the grade of symmetry of the reflectance with respect to the local noon.

The collection relative to 24 June 2009 contains 251 sets of HDRF measurements with solar zenith angle ${\theta }_i$ ranging from 15${}^{\circ }$ to 86${}^{\circ }$. In this work we refer to a “set” as each individual PARABOLA acquisition related to a single ${\theta }_i$ and the full combination of viewing angles. Bands 1 (444 nm), 2 (551 nm), 6 (944 nm) and 8 (1650 nm) show a bowl shaped behavior of the reflectance, symmetric with respect to the local noon and increasing with solar zenith angles ${\theta }_i$. In the other bands, reflectance was characterized by an asymmetric shape and with intraday oscillations inconsistent with previous bands and expectations such as symmetry and smoothness. These were removed from further analysis. Spectral HDRF within 0.2 and 0.6 has been observed in the blue and green channels, while values varying from 0.3 to 0.8 are observed in the near-infrared channels.

For the acquisitions performed in 2014, band 4 (650 nm) and band 7 (1028 nm) were still affected by some artifacts and were removed from the BRDF optimization analysis. The spanned range of sun zenith angles was smaller than in the 2009 dataset, as confined below 45${}^{\circ }$, pertaining to the central part of the day, around local noon. The spectral behavior was similar to that observed in 2009, with an increase of the average reflectance from the green (∼0.30) to infrared channels (∼0.40).

4.2. Bidirectional Reflectance Data from Aircraft (CAR)

During a transect flight from Chico California to NASA Wallops Flight Facility (Virginia) taken on 16 May 2008 (flight number 2006, https://car.gsfc.nasa.gov/flights/data-arctas-2008-flight-2006), a series of reflectance measurements were collected from different altitudes over RRV in coincidence with the Terra satellite overpasses. Two characterizations of the surface from 200 m (in S-E area) and 650 m (N-W) above ground level, and an additional exiting characterization of the western area from higher elevations (up to 3400 m), have been performed before noon, between 18:30 and 19:30 UTC time, for ${\theta }_i$ (${\varphi }_i$) ranging between 24${}^{\circ }$ (138${}^{\circ }$) and 20.5${}^{\circ }$ (156${}^{\circ }$), respectively [21].

The CAR data were processed to identify individual loops representing a complete set of viewing angles (i.e., 0–360${}^{\circ }$ in the azimuth and 0–90${}^{\circ }$ in zenith). This step provided 23 loops over RRV, and 11 of them were selected to perform BRF model inversions over the N-W and S-E areas. The scanning regime provided an angular resolution of 0.5${}^{\circ }$ in viewing zenith and ∼1.5${}^{\circ }$ in azimuth. All data including radiance, irradiance, aircraft positions, times, geometries and other metadata are provided within the CAR product [29,32]. Figure 5 shows the aircraft trajectories as well as the estimated center of the individual quasi-circular loops considered in this study.

As shown in Section 3 (Figure 4), the apparent reflectance ${\rho }^{*}$ (Equation (2)) underestimates the surface reflectance due to the effect of the atmosphere. For the clean conditions expected over RRV (${\tau }_{550}<0.1$) an underestimation of the order of 5–6% in the visible and well below 5% in the NIR is reasonably expected.

4.3. Comparison of PARABOLA and CAR BRF at Similar Zenith Angle Ranges

Table 2. Average and standard deviation of the hemispherical-directional reflectance factor (HDRF) as measured for solar zenith angle ${\theta }_i$ < 25${}^{\circ }$ and viewing zenith angle ${40}^{\circ }\le {\theta }_v\le {50}^{\circ }$, with PARABOLA collected on 24 June 2009, 24 June 2014 (over P7) and 25 June 2014 (over P2), and CAR (16 May 2008). CAR loops were aggregated into two groups. Six loops were performed at an altitude of 645 m above the N-W area (L1-L6), and five loops at 200 m over the S-E area (L7-L11). PARABOLA band numbers are flagged with an * for suspicious data on 24 June 2009, ** for suspicious data in the 2014 set. For these bands, despite the misleading daily evolution that prevented BRF function inversions, the indicative average is reported here for reference.

$\mathit{\lambda }\left(\mathbf{band}\right)$	PARABOLA			$\mathit{\lambda }\left(\mathbf{band}\right)$	CAR
nm	2009 (UND)	2014 (P7, S-E)	2014 (P2, N-W)	nm	S-E Loops	N-W Loops
444 (1)	0.25 ± 0.02	0.33 ± 0.04	0.32 ± 0.03	472 (1)	0.30 ± 0.07	0.30 ± 0.05
551 (2)	0.29 ± 0.03	0.34 ± 0.05	0.33 ± 0.02
PAR (3)	0.44 ± 0.12	0.36 ± 0.03	0.35 ± 0.06
650 (4 **)	0.46 ± 0.07	0.48 ± 0.07	0.49 ± 0.03	682 (2)	0.45 ± 0.05	0.43 ± 0.04
860 (5 *)	0.45 ± 0.10	0.40 ± 0.02	0.41 ± 0.02	870 (4)	0.47 ± 0.05	0.45 ± 0.04
944 (6 *)	0.37 ± 0.05	0.39 ± 0.02	0.40 ± 0.02
1028 (7 **)	0.59 ± 0.14	0.72 ± 0.09	0.71 ± 0.06	1036 (5)	0.48 ± 0.05	0.46 ± 0.04
1650 (8)	0.42 ± 0.04	0.45 ± 0.02	0.46 ± 0.02	1656 (9)

shows the average and standard deviation values of the HDRF values measured by PARABOLA and CAR at similar geometries, with ${\theta }_i<{25}^{\circ }$, and ${40}^{\circ }<{\theta }_v<{50}^{\circ }$. CAR loops were grouped into S-E and N-W sets, and the directional reflectance averaged.

Despite the difference in acquisition dates, the agreement, in particular within the PARABOLA dataset, is reasonable. The spectral measurements taken over two different ground points in 2014 (P2/MODIS-N and P7/M20) and belonging to the N-W and the S-E area, respectively, are cross compatible within the standard deviation, suggesting a certain uniformity of the sampled area. The CAR measurements, collected in 2008 show lower reflectances in the red band and higher in the NIR band with respect to the corresponding values collected by PARABOLA, suggesting a different surface moisture status between the two years.

It is also worth mentioning that by filtering CAR data for nadir viewing (${\theta }_v<5^{\circ }$), we obtained average values very similar to those of Table 2, indicating that the surface is actually quite Lambertian. In fact, for the S-E loops, we obtained 0.32 ± 0.10 (472 nm), 0.46 ± 0.05 (682 nm), 0.48 ± 0.05 (870 nm), and 0.48 ± 0.05 (1036 nm). A similar behavior was observed for the N-W area.

4.4. Optimization of RTLS and RPV parameters using PARABOLA measurements

The inversions of RPV and RTLS functions, were performed for PARABOLA data by limiting the range of angles to within 0${}^{\circ }$ and 70${}^{\circ }$ for ${\theta }_i$, and to within 20${}^{\circ }$ and 70${}^{\circ }$ for ${\theta }_v$. The area relevant to ${\theta }_v<{20}^{\circ }$ contains the reference panel. Its relative surface area to the whole hemispherical footprint is small and, as it is in a confined area below the mast, it might also be affected by mounted hardware and shadows. At higher ${\theta }_i$ the uncertainties increase because of the reduced intensity of the illumination and the deviation from the cosine response of the reference panel (in total sky illumination conditions it is evaluated to be <1% for all the bands above the green band in clear-sky (visibility v∼ 100 km), with peaks of 2.5% for the blue band in hazy conditions (v∼ 20 km), according to [33]). Hence, excluding these geometries from the analysis should reduce artifacts and uncertainties.

By using a subset of the full data in the fitting procedure, we stressed the capability of each function to replicate the parameters as obtained from all data, which were presumably internally correlated. Given the restriction on zenith angles provided above, the optimizations were performed by using:

• All the measurements collected every 3 min over the entire day (full);
• Subsetting measurements by balancing the contribution of measurements with different ${\theta }_i$, i.e., avoiding duplicates of reflectance collected around local noon with an almost constant ${\theta }_i$ (balanced);
• Only the principal plane of reflection, i.e., ${\varphi }_i\sim {\varphi }_0$) (pp);
• Only ${\theta }_v\sim {45}^{\circ }$ as in Table 2 (45);
• A limited number of random measurements from 4 to 128;
• Set-by-set inversion and computing the statistics (set-by-set average).

The fitting process consisted of two steps, the first using angular filtered data, the second excluding the outliers identified by points with residuals greater than ±1.5 for the inter-quartile range IQR [44].

The results of these fits for the green band reflectance collected on 24 June 2009 are reported in

Table 3. Bi-directional reflectance factor (BRF) parameter retrieval for the green band (551 nm) for Ross Thick-LiSparse (RLTS) and Rahman–Pinty–Verstraete (RPV) functions using different subsetting of the RRV PARABOLA datasets, as described in the text. The uncertainty reported for sampling methods (5-6) associated with the median value, represents half the interquartile range (IRQ/2). The last column values (BHR) are obtained by integrating the bi-directional reflectance distribution functions (BRDF) functions with the inverted parameters. Sampling methods (Smp) are 1. full, 2. balanced, 3. pp, 4. ${\theta }_v\sim$45${}^{\circ }$, 5. random and 6. set-by-set for the 24 June 2009. For the 2014 dataset only the set-by-set inversion statistics are reported.

Function	Smp	${\mathit{f}}_{\mathbf{iso}}$ or ${\mathit{\rho }}_0$	${\mathit{f}}_{\mathbf{vol}}$ or k	${\mathit{f}}_{\mathbf{geo}}$ or $\mathbf{\Theta }$	RMSD	N$·{10}^3$	BHR
Date: 24 June 2009
RTLS	1	0.372	0.149	0.062	0.015	142	0.315
	2	0.375	0.139	0.063	0.016	65	0.315
	3	0.364	0.153	0.058	0.012	4	0.313
	4	0.387	0.121	0.070	0.014	12	0.314
	5	0.365 (12)	0.181 (11)	0.055 (10)	–	0.016	0.324 (30)
	6	0.372 (25)	0.141 (32)	0.062 (13)	–	–	0.313
RPV	1	0.170	0.750	−0.121	0.018	136	0.319
	2	0.172	0.760	−0.119	0.018	67	0.319
	3	0.170	0.722	−0.114	0.018	4	0.324
	4	0.181	0.812	−0.109	0.017	12	0.317
	5	0.164 (8)	0.741 (39)	−0.132 (15)	–	0.016	0.315 (40)
	6	0.179 (3)	0.812 (11)	−0.110 (18)	–	–	0.317
Date: 24 June 2009 (P7)
RTLS	6	0.317 (5)	−0.16 (7)	−0.01 (1)	<0.02	–	0.300
RPV	6	0.36 (1)	1.04 (3)	0.17 (1)	<0.02	–	0.337
Date: 25 June 2009 (P2)
RTLS	6	0.33 (1)	−0.22 (12)	0.00 (3)	0.02	–	0.288
RPV	6	0.36 (1)	1.08 (2)	0.16 (1)	<0.02	–	0.333

. Method 1 has been assumed to provide the benchmark set of retrieved parameters and BHR since it used the whole dataset. Methods 2 to 4 produced compatible parameter inversions even using fewer BRF measurements. The major difference observed in the BHR was given by method 3 in the case of the RPV model, and by method 5 (random) for the RTLS model. The root mean squared difference, defined as $RMSD=\sqrt{{\sum }_{i=1}^N{(BR{F}_{mod}-BR{F}_{mea})}^2}$/N, represents the capability of the model to fit the measurements. It ranged between 0.012 and 0.018, establishing the typical value for the uncertainty associated with the fitting step.

In order to calculate the performance statistics, for the random cases only, ten random extractions of samples with increasing dimensions of 4, 16, 32, 64 and 128, have been performed. The optimizations were applied for each case, and the average and standard deviation of the parameters provided. Table 3 shows only the results for N = 16 (sampling method number 5). The dispersion $\sigma$ of the retrieved parameters decreases with the number of samples N. Isotropic parameters ${\rho }_0$ and $f_{iso}$ converge to values with relative errors (${\sigma }_x/\overline{x}$) lower than 10% as is apparent from N = 8 onwards, but to obtain relative errors below 5%, a minimum of ∼30 measurements are required. The RTLS’s volumetric ($f_{vol}$) and geometric ($f_{geo}$) parameters, and the RPV’s $\mathrm{\Theta }$ (RPV) parameter, presented relative errors between 10% and 20% until the maximum N considered in this study (N = 128). The k-Minneart parameter of RPV was retrieved with an error below 10% only for N≥ 16.

The BRDF optimization has been performed even for a set-by-set basis (case 6) to assess the representativeness of each single acquisition with respect to the daily inversions (Figure 6). Despite a pronounced daily variation of the RTLS parameters with respect to the RPV ones, the BHR presents similar behavior with increasing reflectance at higher sun zenith angles (${\theta }_i>{60}^{\circ }$) on 24 June 2009. The BHR (white-sky albedo) should not depend on ${\theta }_i$ because it is theoretically removed by the bi-hemispherical integration over viewing and illumination angles. However, the dependence that was found, particularly evident for RPV, might be related to the effect of diffuse radiation that cannot be neglected with increasing ${\theta }_i$ and should be considered during the optimization of the BRF function. The median values of the parameters are reported in Table 3 (case 6) and appear compatible with other approaches if ${\theta }_i<{60}^{\circ }$ is considered. Unfortunately, the red and NIR values cannot be discussed for 24 June 2009 because of the artifacts observed in the time series for this particular date. For the other valid bands (444, 944, 1650 nm), the BRF reflectance residuals are normally distributed with biases around zero and with $\sigma \le 0.02$ (results for case 3). Two additional datasets were collected on June 2014 at stations P2 and P7. Measurements were taken around local noon with a maximum ${\theta }_i\sim {40}^{\circ }$. Despite similar BHR values, the anisotropy was less marked as evidenced by a $k\sim 1$ for the RPV function or $f_{geo}\sim 0$ for the RTLS function.

4.5. Optimization of RTLS and RPV Functions Using CAR Measurements

Since the green channel is missing in CAR, a discussion fully consistent with PARABOLA results was not possible. Nevertheless, we decided to perform the same optimization described in the previous section using CAR data in the red band. Figure 7 shows the retrieved values for each loop of RTLS and RPV parameters, assuming a flat and homogeneous surface [29]. The measurements were taken for a ${\theta }_i$ that varied between 19${}^{\circ }$ and 25${}^{\circ }$, and inversions were performed limiting the range of ${\theta }_v$ below 60${}^{\circ }$.

Figure 7 reveals slight differences between parameters retrieved from individual loops, with a lower variability observed for the isotropic parameters $f_{iso}$ and ${\rho }_0$. Loops 1-6 are relevant to the S-E sector scanned from a lower altitude (∼200 m) with respect to the N-W area, which was scanned from an altitude of ∼650 m. This is likely to be the explanation of the higher BHR variability in the S-E loops, which scan the surface from a lower altitude and thus resolve some smaller scale variability.

For the RPV function the resulting BHR in the red band varies from $0.465\pm 0.007$ to $0.48\pm 0.02$ from higher (N–W) to lower (S–E) observation altitudes. These albedo values are comparable with the HDRF statistics reported in Table 2 that were computed by limiting the observing zenith angles within the 40${}^{\circ }$ to 50${}^{\circ }$ range. The RTLS parameter optimization shows a similar behavior for the two areas and provides consistent values of the BHR.

5. Evaluating the Reflectance Consistencies through a 3-D Model Based Approach

In order to overcome the difficulties encountered in the collection of data of which quality is not always traceable, a modeling approach based on a validated 3-D Monte Carlo radiative transfer model was used. We have chosen to simulate reflectance measurements using the Raytran program [26,37]. Among several virtual measurements, it implements a directional radiance/reflectance model by simply defining the position of the sensor with respect to the surface, its field-of-view, and the angles of observation. In its current release the code is not suitable to simulate atmospheric radiative transfer and we performed the simulations in black-sky conditions only.

5.1. Virtual Scene Setup

The definition of a land-surface virtual scene can be considered the most difficult and time consuming activity in the framework of 3-D radiative transfer modeling. It requires knowledge of the structural and optical properties of the world under investigation, and their adequate abstraction, in order to fit the limits set by calculation time and storage capabilities. Because of the lack of data regarding structural features such as the surface roughness over RRV, we implemented an approach based on medium and high resolution satellite imagery to create a patchy surface of squares, each with its own bidirectional reflectance defined by a triplet of the modified RPV function parameters [45]. This analytical approach is known as parametric meshing [46].

We used the MISR land surface product (MIL2ASLS.002) for 15 July 2009 to map the RRV area with the medium resolution information (at 1.1 km) for the k and $b_M$ anisotropic parameters, the latter replacing $\mathrm{\Theta }$ in the modified version of the RPV function adopted by MISR. In order to refine the spatial resolution we adjusted the value of ${\rho }_0$ to match the spectral surface reflectance $R^{\prime }$ obtained by Landsat-5 Thematic Mapper (TM) for the nearest clear-sky day (14 July 2009). During this process, k and $b_M$ were simply down-scaled (with a bi-linear approach) from 1.1km resolution to 30 m TM resolution, to give $k^{\prime }$, $b_M^{\prime }$. Then, ${\rho }_0^{\prime }$ at 30 m was obtained by inverting the $mRPV$ BRF formulation for each Landsat pixel as

$$R^{\prime }=mRPV({\rho }_0^{\prime },k^{\prime },b_M^{\prime })={\rho }_0^{\prime }\frac{{\left(cos{\theta }_{i}cos{\theta }_v\right)}^{k^{\prime }-1}}{{(cos{\theta }_i+cos{\theta }_v)}^{1-k^{\prime }}}{e}^{-b_M^{\prime }cosg}H\left(G\right)$$

(5)

where, $cosg=cos{\theta }_{i}cos{\theta }_v+sin{\theta }_{i}sin{\theta }_{v}cos({\varphi }_i-{\varphi }_v)$ is the phase angle between the incoming and the outgoing direction, and the hot–spot function is defined by $H\left(G\right)=1+\frac{1-{\rho }_0}{1+G}$, with $G^2={tan}^2{\theta }_i+{tan}^2{\theta }_v-2tan{\theta }_{i}tan{\theta }_{v}cos({\theta }_i-{\theta }_v)$ [12,45]. A complete discussion of the algorithm and of the BRF function, with respect to the original RPV model, is given in the MISR Land Product ATBD [47].

To simplify the inversion, we assumed the Landsat viewing direction (${\theta }_v$) equal to 0${}^{\circ }$, an approximation based on viewing angles ranging between ±7.5${}^{\circ }$, with a nadir-adjusted reflectance that likely varies by less than 5% with respect to the ETM+ measurements [48]. Nevertheless, it should be noted here that the intention was not to compare model results with actual measurements, but simulations performed at different scales based on a reasonable virtual surface model. With this approximation, being $cos{\theta }_v\sim 1$, and $tan{\theta }_v\sim 0$ (e.g., $G\sim tan{\theta }_i$), Equation (5) simplifies to a second order equation to be solved with respect to ${\rho }_0^{\prime }$

$$\left(1+\frac{1-{\rho }_0^{\prime }}{1+tan{\theta }_i}\right){\rho }_0^{\prime }\sim \frac{R^{\prime }}{\frac{{cos}^{k^{\prime }-1}{\theta }_i}{{(cos{\theta }_i+1)}^{1-k^{\prime }}}{e}^{-b_M^{\prime }cos{\theta }_i}.}$$

(6)

Using the root of Equation (6) falling within the range [0:1], the RRV area of interest was meshed with a grid of 30 × 30 m squares with spectral reflectance properties defined in terms of $({\rho }_0^{\prime },k^{\prime },b_M^{\prime })$ triplets, that can be easily ingested by the Raytran code. This process implies 100 squares to represent a 300 × 300 m area, and ${10}^4$ squares to cover a 3 × 3 km area, that is approximately the maximum area of our interest to compute CAR simulations taken from a height of 500 m, and a maximum viewing angle of ∼70${}^{\circ }$. Figure 8 shows the surface reflectance $R^{\prime }$ as measured from Landsat-5 TM and the isotropic parameter ${\rho }_0$ as obtained by the above mentioned methodology. The median values over the plotted area are ${\overline{R}}^{\prime }=0.38\pm 0.06$, ${\overline{\rho }}_0=0.20\pm 0.05$, $\overline{k}=0.86\pm 0.05$ and ${\overline{b}}_m=-0.18\pm 0.10$.

This virtual scene can be considered as an independent scenario that can be used to perform simulations relevant to PARABOLA, CAR or satellite sampling. It is not necessarily compliant with the original reflectance value as it is obtained by combining two products with different resolutions, under certain assumptions and bi-linear smoothing.

The true value of reflectance is known a-priori because it is analytically defined for each elementary square (at 30 m resolution), thus so is the performance of each sampling approach.

5.2. Simulations of PARABOLA and CAR Reflectances

As the implementation of the RRV here is based on a parametric meshing at 30 m resolution, it is expected that the variability associated with a specific sampling strategy depends only on the number of squares scanned by the elliptical instrument footprint. It is a widely accepted assumption, to consider the linear dimension of the area sampled by an instrument with a hemispherical field of view, to be nearly 10 times the height of the instrument installation [49,50]. Cescatti et al. [17], discussing the FLUXNET implementation for albedo measurements, observed that the actual footprint of an albedometer depends on its height above the canopy. For a typical installation of 5–10 m above the canopy top, 80% of the signal originates within 10–20 m from the tower location. Adams et al. [51], using a 3-D modeling approach, have shown how having multiple albedo observations increases the possibility of meeting GCOS requirements over canopies.

In the PARABOLA case, as the useful range of directional observation was limited to ${\theta }_v<72.5^{\circ }$, for an installation height of 3m, the radius of the circular area $d_M$ sampled by all measurements is about 9 m (

Table 4. Distance $d_M$, between the optical axes and local vertical intercepts with the ground surface. Values are given for increasing viewing zenith angles ${\theta }_v$ and installation heights h (in bold). The second part of the table gives the gFOV elliptical area for the same observing conditions. For h of 3, 10 and 30 m the instrument field-of-view f is assumed to be 5${}^{\circ }$ to mimic PARABOLA, while for $h>=200$, f is set to 1${}^{\circ }$ to mimic CAR.

${\mathit{\theta }}_{\mathit{v}}$	$\mathbf{tan}{\mathit{\theta }}_{\mathit{v}}$	${\mathit{d}}_{\mathit{M}}$ (m)						Individual Footprint gFOV = $\mathit{\pi }\mathbf{ab}$ (m${}^2$)
-	-	3 m	10 m	30 m	200 m	500 m	1 km	3 m	10 m	30 m	200 m	500 m	1 km
1	0.017	0.1	0.2	0.5	3.5	8.7	17.5	0.2	2.4	21.6	38.3	239	957
5	0.087	0.3	0.9	2.6	17.5	43.7	87	0.2	2.4	21.8	38.7	242	968
10	0.176	0.5	1.8	5.3	35.3	88.1	176	0.2	2.5	22.5	40.1	250	1002
30	0.577	1.7	5.8	17.3	115	289	577	0.33	3.7	33.2	58.9	368	1473
45	0.707	3	10	30	200	500	1000	0.6	608	61.1	108	676	2707
60	1.732	5.2	17.3	52	346	866	1730	1.73	19.3	174	306	1914	7658
72.5	3.17	9.5	31.7	95.1	634	1587	3170	8.1	90	808	1409	8805	35,223

).

Therefore, in these conditions, PARABOLA samples an area with linear dimensions of ∼30 m, corresponding to the spatial resolution of our surface model. Hence, using our virtual model, the optimization of RTLS and RPV parameters are always performed over a set of reflectances originating from 1, 2 or 4 pixels only, depending on the position of the instrument with respect to the model grid.

For that reason, it was decided to associate with each individual simulation an additional random and normally distributed error based on the spectral values of the root mean square error (RMSE) obtained by fitting the RPV function to real field measurements taken in the corresponding solar angular range ($\sigma =0.025$). The radius $d_M(72.5^{\circ })$ increases to 32 and 95 m, with an area ($\pi d_M^2$) of about 3.5 and 31 equivalent Landsat pixels (900 sqm), for installations of 10 and 30 m, respectively.

On the other hand, the individual ground footprint area $gFOV$, varies consistently with installation height (PARABOLA) or altitude flight (CAR) and the viewing zenith angle. Considering the major a and minor b axes of the projected footprint as $a=h(tan{\theta }_v+FOV/2)-tan({\theta }_v-FOV/2))$ and $b=2htan(FOV/2)sec{\theta }_v$[21], expressed in terms of Landsat resolution, the ground footprints vary from sub-pixel to multiple pixel dimensions while either ${\theta }_v$ or h increase (Table 4). Nominally, at ${\theta }_v=72.5^{\circ }$, gFOV for PARABOLA is 1/100, 1/10 or 1 equivalent Landsat pixel, for 3, 10 and 30 m, respectively. For CAR (which has a narrower FOV), the values are approximately 1 and 10 Landsat pixels at flight elevations of 200 m and 500 m, respectively.

Therefore, the optimization process of a BRF function, performed by using measurements with different surface sampling areas, can lead to considerable errors in heterogeneous ground cover, as re-gridded data are influenced by information from surrounding pixels [52].

The reflectance for PARABOLA was simulated over the eight ground points marked in Figure 5, for 15 July 2009 from 7:00 to 17:00 local time. Within this period the sun zenith ranged between 17${}^{\circ }$ (at noon) to 67${}^{\circ }$. Simulations were limited to the red band, for which the virtual scene was implemented.

To perform PARABOLA simulations from 3 m, the full virtual scene was subsetted over a squared area of 200 × 200 m${}^2$ centered on the instrument position to maintain a reasonable computation time. One million photons equally distributed over the area were generated to produce an average of 25 photons/m${}^2$. Figure 9 (left panel), shows a representation of the BRF simulated for PARABOLA centered on P6 in Figure 5 positioned at 10 m.

The simulation of CAR’s BRF were performed for eleven counterclockwise loops with 1 km radius, centered over L1-L11 (Figure 5). Figure 9 (right panel) illustrates a CAR sampling from an altitude of 645 m. The BRFs have been rendered in the image as ellipse-shaped footprints. Simulations were performed up to a ${\theta }_v$ of 60${}^{\circ }$ from a maximum flight height of ∼1 km. We performed CAR BRF simulations for altitudes from 200 m to 1 km. The same number of incoming photons (${10}^6$) produced a density of 0.25 photons/m${}^2$. Each Landsat equivalent pixel was then sampled by 225 photons.

5.3. Inversion of BRF Using the PARABOLA and CAR Simulations

The RTLS and mRPV retrievals were run using PARABOLA and CAR measurement simulations performed for ${\theta }_v$ between 20${}^{\circ }$ and 60${}^{\circ }$, to be compliant with the inversions performed using real measurements described in Section 3.

The simulations were performed at three different heights for PARABOLA (3, 10 and 30 m) and four heights for CAR flights (200, 500, 645, and 1000 m), for each of the P1-P8 sites (Figure 5), and L1-L11 sites, respectively. The maximum distance D sampled from the local vertical ground point for a ${\theta }_v$ of 60${}^{\circ }$, for example, the limit used for the kernel function inversions, varied from 5 to 50 m for PARABOLA, and from 350 to 1730 m for CAR simulations ($D=h\sqrt{3}$). Figure 10 shows the retrieved parameters for each PARABOLA point and CAR loop as simulated at the different levels. The RMSE of the RTLS and RPV models in representing the BRF, and the corresponding BHR, are also reported.

Excluding results over P1, the values of $f_{iso}$ and BHR relevant to the N-W area (P2–P6, and L1–L6) are internally consistent within ±0.03, for each level the measurement was taken.

For CAR simulations only, the albedo of the S-E area (white symbols), is evidently greater with respect to the N-W one. This behavior agrees with real data shown in Figure 7 relevant to 2008 data. PARABOLA measurements over sites P7 and P8 present quite different reflectance properties and albedo. These measurements are not compliant with the evaluation made by CAR (L7–L11) over the same area, which, excluding L7, are internally compatible (L8–L11) in terms of $f_{iso}$, $f_{geo}$ and BHR. The volumetric parameter $f_{vol}$ retrieved by CAR measurements shows the maximum variability (within 0 and 0.25), and it drives the RTLS fitting performance, as the RMSE mimics its patterns in the plot, and the coefficient of determination $r^2$ between RMSE and $f_{vol}$ varies from 60% to 86% depending on the instrument altitude. On the other hand, $f_{geo}$ characterizes, very markedly, the differences between the N-W and S-E areas with a clear step in each inversion group above 500 m.

Table 5. Coefficient of determination $r^2$ of the RTLS parameters with respect to the RMSE (first set) and the BHR (second set).

h (m)→	3	10	30	200	500	645	1 km
	$f_i$ vs. RMSE ($r^2$)
$f_{iso}$	0.43	0.33	0.01	0.22	0.51	0.60	0.69
$f_{vol}$	0.08	0.00	0.84	0.59	0.75	0.74	0.86
$f_{geo}$	0.00	0.01	0.15	0.34	0.80	0.84	0.86
	$f_i$ vs. BHR ($r^2$)
$f_{iso}$	0.75	0.71	0.67	0.76	0.90	0.77	0.91
$f_{vol}$	0.01	0.03	0.01	0.21	0.47	0.34	0.90
$f_{geo}$	0.01	0.03	0.11	0.65	0.79	0.90	0.91

shows $r^2$ as calculated between the RMSE (and BHR) with respect to each RTLS kernel parameter. As we notice in

Table 6. Statistics of the RLTS and RPV parameters as obtained by optimization of the simulated PARABOLA (h of 3, 10, 30 m) and CAR (200, 500, 1000 m) measurements. Values obtained for PARABOLA over P7 and P8 (S–E) were not averaged and are reported independently separated by the “/” symbol. Points P2–P6 were used to compute the statistics over the N–W area for PARABOLA. Loops L1–L6 and L8–L11 were used to compute the statistics for CAR over the N–W and S–E sectors, respectively. The numbers in parenthesis indicate the standard deviation and refer to the last digit of the average value.

	RTLS Function		RPV Function
$\mathbf{h}$(m)	N–W	S–E	N–W	S–E
	$f_{iso}$		${\rho }_0$
03	0.382 (6)	0.376/0.351	0.206 (4)	0.212/0.202
10	0.381 (8)	0.376/0.347	0.204 (6)	0.215/0.199
30	0.383 (5)	0.376/0.346	0.205 (6)	0.214/0.198
200	0.372 (3)	0.401 (5)	0.194 (12)	0.255 (8)
500	0.374 (3)	0.406 (3)	0.184 (8)	0.257 (11)
1 km	0.368 (3)	0.402 (3)	0.181 (7)	0.256 (6)
	$f_{vol}$		k
03	0.113 (8)	0.115/0.120	0.839 (5)	0.840/0.842
10	0.130 (10)	0.118/0.120	0.830 (7)	0.837/0.816
30	0.127 (14)	0.131/0.135	0.836 (13)	0.821/0.807
200	0.144 (50)	0.066 (32)	0.826 (34)	0.882 (16)
500	0.145 (34)	0.051 (33)	0.819 (26)	0.878 (31)
1 km	0.183 (19)	0.057 (12)	0.776 (22)	0.888 (7)
	$f_{geo}$		$\mathrm{\Theta }$
03	0.052 (1)	0.044/0.038	−0.064 (3)	−0.053/−0.045
10	0.052 (4)	0.043/0.039	−0.067 (4)	−0.047/−0.048
30	0.054 (2)	0.044/0.037	−0.065 (6)	−0.051/−0.046
200	0.053 (5)	0.039 (3)	−0.076 (24)	−0.013 (11)
500	0.060 (3)	0.043 (2)	−0.092 (17)	−0.010 (17)
1 km	0.057 (2)	0.041 (4)	−0.099 (13)	−0.009 (8)
	BHR		BHR
03	0.332 (5)	0.337/0.321	0.330 (6)	0.335/0.317
10	0.334 (6)	0.339/0.316	0.332 (7)	0.338/0.320
30	0.332 (7)	0.340/0.321	0.331 (11)	0.343/0.320
200	0.326 (6)	0.360 (5)	0.321 (1)	0.365 (3)
500	0.319 (4)	0.357 (4)	0.314 (4)	0.368 (3)
1 km	0.324 (2)	0.356 (4)	0.322 (5)	0.363 (4)
Ref BHR	0.322 (15)	0.348 (12)	0.322 (15)	0.348 (12)

, the various quantities are correlated, and $r^2$ can be interpreted as how much the parameter $f_i$ (or its corresponding parameter for the RPV function) influences the albedo or the uncertainty in its retrieval (RMSE).

For PARABOLA, the RMSE variability was given by the error artificially introduced in the forward simulations because of the limited area sampled ($D\sim$ 50 m at ${\theta }_v$ of 60${}^{\circ }$). For CAR the results are more significant, with an increase of 20% to 70% of the RMSE variability explained by the $f_{iso}$ for simulations performed at 200 and 1000 m, respectively. Nevertheless, both $f_{vol}$ and $f_{geo}$ may explain even better the RMSE fitting performance with a greater and similar strength (75 to 85% above 200 m). The isotropic parameter explains 70% of the BHR variability for PARABOLA, and up to 90% of its variability for CAR measurements taken from 1 km altitude.

The volume and geometric kernel parameters do not influence the albedo appreciably for in-situ measurements (<3%). Only for PARABOLA placed at 30 m, where 11% of the BHR variability can be attributed to $f_{geo}$ ($r^2=0.11$).

For CAR, $f_{iso}$ and $f_{geo}$ affect the value of the BHR by 21–90% and 65–90%, respectively, for increasing altitudes. It is evident that, though this is a qualitative interpretation, the co-variance between different parameters is not null, and they are cross correlated, hence not independent from each other. This is a property that is reinforced with observation height, or the area sampled by the instrument.

The RMSE accounts for the capability of the RTLS function to represent the simulated bidirectional reflectance features. The RMSE for PARABOLA measurements is around $5\times {10}^{-3}$ as it is merely determined because an artificial random error was added to the forward simulations to account for the limited sampling footprint of the PARABOLA with respect to the spatial variability of the virtual scene, as previously discussed. The RMSE are constant even for observations performed at 30 m (footprint radius is 52 m at 60${}^{\circ }$ or ∼$3\times 3$ pixels) indicating that locally the virtual surface reflectance varies quite smoothly.

The RMSE increases for CAR measurements (0.010 to 0.025, see the right panel of Figure 10) because of the sampling strategy and the variability of the surface sampled over a loop. For the particular case of RRV the altitude of the aircraft does not affect appreciably the values of the retrieved parameters and RMSE, as the same patterns within each single L1–L11 group in Figure 10 are constantly met above h = 500 m, with small variations only for the group of inversions related to measurements taken from an altitude h of 200 m.

Table 6 shows the statistics for each block of inversions presented in Figure 10. The isotropic parameter statistics obtained at different heights and pertaining to different ground sectors were observed to be compatible over the N-W sites (P2–P6, L1–L6), where it ranges between 0.37 and 0.38 with a standard deviation of ∼0.03, for each installation height. Over the S-E area the PARABOLA measurements give similar parameters ranging between 0.35–0.37 (P7, P8) and 0.40 (CAR).

The volumetric kernel is constrained between 0.10 and 0.15 for the PARABOLA inversions, while it varies for each CAR loop, decreasing toward the S-E sector from 0.20 to 0.05, with an average value of $0.14\pm 0.05$ for measurements taken below 1000 m. While values are slightly different, the geometric parameter presents a similar shape evolution from P1–P8 and L1–L11 with higher values over the N–W sector (>0.05) with respect to the S–E one (0.37–0.44). Similar results were obtained from the optimization of the RPV function parameters.

The reference values for the bi-hemispherical reflectance (Ref BHR) were computed at the pixel level over the various sub-samples used for the loop simulations, and using the original mRPV model used to define the scene. The reference BHR was $0.322\pm 0.015$ for L1–L6, and $0.361\pm 0.012$ over the L7–L11 areas (Figure 10). Table 6 shows that the BHRs retrieved by means of RTLS or RPV parameterizations are always compatible within $1\sigma$ over the N-W sector (for CAR and PARABOLA), and within $2\sigma$ over the S-E one for CAR only. The BHRs obtained from the PARABOLA inversions over P7 and P8 are compatible within 1%, but presented lower values with respect to the CAR retrievals, and evidently are not representative of the wider area sampled by CAR.

6. Conclusions

In-situ PARABOLA and CAR measurements collected during different campaigns and different times, was acquired and investigated. The lack of temporal collocation prevented a direct comparison of the surface reflectance properties. They were expressed in terms of the RTLS [13] and RPV [12] BRF kernel models. We evaluated the effect of the different spectral response functions (SRF) of the red and NIR bands of PARABOLA and CAR, by showing that, over the Railroad Valley the CAR likely overestimates PARABOLA reflectance by up to 3% in the red band, while no differences arise from the SRF in the NIR band (0.2%), due to the flat spectral reflectance.

For the peculiar low-turbidity conditions of RRV (${\tau }_{550}=0.05$) we estimated that, for a sun zenith angle ${\theta }_i={50}^{\circ }$, the differences between HDRF and BRF ranges up to 0.02 (over 0.5, hence ∼4%) in the hot-spot region, for the green band only. In the NIR the differences between the HDRF and BRF (${\tau }_{.55}=0.05$) forward simulations are negligible and we can neglect the atmosphere. Then, by directly using the measurements (HDRF) in the process of inversion of the BRF models, we showed that the ${\rho }_0$ isotropic RPV parameter could be retrieved with an error as high as 4% in the green band at ${\theta }_i={20}^{\circ }$. Moreover, we highlighted a linear relation between the retrieved ${\rho }_0$ and the spectral diffuse ration d, which allowed us to define an easy method to correct the parameters retrieved in ambient illumination conditions (blue-sky).

We also investigated the effect of different angular filtering applied to the full set of PARABOLA measurements, with the aim to reduce the dimension of the dataset, and evaluate the effect on the parameters’ optimizations (Table 3). It has been shown that, using different sampling schemes the BHR in the green band of PARABOLA (∼0.31) can be retrieved with an error as lower as 2–3%, using all data taken during a full day. Performing retrieval of parameters on a set-by-set case, without a sufficient sampling of sun angles (case 6) resulted in an erroneous behavior of the retrieved BHR, which was observed to vary with ${\theta }_i$, while it should remain constant during the day (Figure 6). This affected in particular the RTLS model inverted from PARABOLA observations collected at sun angles ${\theta }_i>{50}^{\circ }$.

Subsequently, we presented a model based approach to investigate the level of compatibility between ground-based and aircraft-based sampling. Because of the lack of detailed field measurements such as the micro-topographic soil structure, digital-elevation model and extended hyper-spectral measurements, we created a virtual scene for Railroad Valley by combining collocated hectometric (MISR) and decametric (Landsat-5 TM) atmospherically corrected land surface reflectance products in the red band.

Such a method can be easily replicated over other areas to assess the uncertainty that would affect similar field campaigns. The method implements a simple down-scaling of BRF function parameters from the medium resolution scale of MISR to the higher resolution of the thematic mapper (TM), by assuming a bi-linear variation of the geometric ($b_M$ ) and volumetric (k) parameters, and calculating the isotropic parameter (${\rho }_0$) by just inverting the modified $mRPV$ equation (adopted in MISR products) with respect to it, and using the TM reflectance as a constraint. This process allowed us to set a parametric meshing of the whole RRV area of interest at 10 m resolution.

We then performed directional reflectance measurements to mimic PARABOLA and CAR sampling approaches, for different installation heights (3–30 m and 200–1000 m for the two systems, respectively) and a fixed sun zenith angle (${\theta }_i={28}^{\circ }$). In this part of the work, we got rid of the spectral differences between the equivalent instruments’ band. The spectral properties of the surface were intrinsically defined by the surface parametric meshing and assumed to be common to the different instruments. Therefore, the residual errors were only related to the difference between the various sampling approaches and to the area averaging process.

The assessment of in-situ and aircraft retrieval of BHR was made with respect to the area averaged values as obtained by calculating BHR over all pixels by means of Equation (5) and averaging them up. We showed in Figure 10 and in Table 6 that retrieved BHR matches the expected area average value at all flight altitudes over the N-W area, within 1$\sigma$ (up to 2%) and that PARABOLA is, as expected, representative of local properties only, as it overestimates and underestimates the area averaged values over the N-W and S-E areas respectively, by ∼3%. The results for BHR do not depend on the usage of either RTLS or RPV functions, confirming their robustness in representing anisotropic features of Railroad Valley Playa. In this study we focused on ground points and aircraft CAR loops as originating from real experiments, highlighting the limitation of non-overlapping samplings. Further studies should be envisaged to account for virtual simulations that guarantee overlapping sampling and direct comparisons of surface based and aircraft directional reflectance measurements.