Modeling and intercomparison of field and laboratory hyperspectral goniometer measurements with G-LiHT imagery of the Algodones Dunes

Abstract. We compare field hyperspectral bidirectional reflectance distribution function (BRDF) measurements acquired by a hyperspectral goniometer system known as the goniometer of the Rochester Institute of Technology (GRIT) during an experiment in the Algodones Dunes system in March 2015 with NASA Goddard’s light detection and ranging, hyperspectral, and thermal imagery of the site acquired during the experiment. We augment our field spectral data collection with laboratory hyperspectral BRDF measurements of samples brought back from the Algodones Dunes site using GRIT and our second-generation goniometer GRIT-two (GRIT-T). In these laboratory experiments, we vary geophysical parameters such as sediment density and grain size distribution of the sediments that would typically impact observed BRDF with the goal of extending the range of applicability of our resulting BRDF spectral libraries. Geotechnical measurements on site confirm the variability of geophysical parameters such as density and grain size distributions within the dune system, and measurements with GRIT and GRIT-T demonstrate the impact on observed spectral variation. By augmenting field spectral libraries with laboratory BRDF, we show that a greater proportion of the dune system is more faithfully represented in the expanded spectral library. Beyond developing appropriate calibration data for airborne and satellite imagery of the Algodones Dunes, laboratory and field studies also support goals to develop reliable retrieval methods for geophysical quantities such as sediment density directly from spectral imagery. We consider approaches based on the Hapke model. Our approaches use the invariance of the observed functional forms of the single scattering phase function, which must be invariant to differences in the illumination geometry. Fill factor is retrieved and correlates with expected direct measurements of sediment density in a laboratory setting.


Introduction
In March 2015, a joint field experiment and airborne measurement campaign was conducted in the Algodones Dunes, California.The Algodones Dunes has been selected by NASA as a potential site for satellite intercalibration, and the overall goal of the experiment was to better characterize the spectral variability of the dunes.The experiment was a multiinstitutional effort involving NASA, South Dakota State University, the University of Arizona, the University of Lethbridge, and the Rochester Institute of Technology (RIT).These dunes have significant topographic relief.This fact, combined with natural variations anticipated within a large dune system, such as Algodones, made the direct measurement of the bidirectional reflectance distribution function (BRDF), which characterizes the angular variation in scattered light from the surface for a given illumination geometry, of paramount importance.
Within the context of the larger team effort, our goals were to provide representative hyperspectral BRDF measurements at a number of sites within the dune system and to characterize the geophysical properties of the surface that might influence observed variations in BRDF at these locations and throughout the dune system.With a longer-term view of improving representative BRDF of the Algodones Dunes system, we captured sediment samples from multiple sites within the dune system, both where our BRDF field measurements occurred as well as at other sites that represented surface variations in geophysical properties of the dunes where it had not been possible to bring our field goniometer system within the time constraints of the field experiment.Our plan was to bring these samples back to the laboratory and deliberately manipulate their geophysical properties to mimic variations that might be expected within the Algodones Dunes system.Thus, we report here our laboratory BRDF measurements in which we have manipulated individual geophysical variables such as sediment density and grain size of the sediment samples as well as the illumination geometry.We describe the accuracy of using these BRDF libraries from manipulated specimens in comparison with NASA Goddard's light detection and ranging (LiDAR), hyperspectral, and thermal (G-LiHT) 1 imagery.We also report on our efforts to invert radiative transfer models based on the work of Hapke to directly estimate sediment fill factor from the observed BRDF.
Representative BRDF libraries play a critical role in calibration of both airborne and satellite systems.For example, these type of libraries have been used within scene simulation models, such as DIRSIG, 2 which can predict the expected scene top-of-atmosphere radiance.DIRSIG currently incorporates standard BRDF models, such as the Ross-Li 3 model.While models such as the Ross-Li model have been widely used in NASA satellite data products, they have a limited number of available parameters to model variations in expected BRDF of the surface.Another goal of this work, therefore, is to undertake analysis of more complex models, such as those developed by Hapke, 4 for future inclusion in DIRSIG and use in satellite and airborne sensor calibration.
The rest of the paper is organized as follows: Sec. 2 outlines methods used to acquire data in the field and in the laboratory where samples were manipulated.In this section, we also outline approaches to radiative transfer modeling and inversion based on Hapke's model. 4In Sec. 3, we report comparisons between the spectral libraries acquired in the field and developed in the laboratory to airborne G-LiHT imagery acquired over the dune system.In this section, we also examine the results of our approaches to geophysical parameter inversion.These approaches rely on a comparison of BRDF data acquired at two different illumination geometries using the invariance of the single scattering properties of the surface.In this context, we consider the problem of inversion of the Hapke model and a variant of this model to obtain the fill factor, and compare to laboratory measurements of density, which are expected to correlate with the fill factor.Finally, in Sec. 4, we provide a summary and conclusions and describe the next steps needed to eventually transition the approach to retrieval of the fill factor directly from the imagery time series.

Field Measurements of the Angular Dependence of Hyperspectral Reflectance at the Algodones Dunes
During the March 2015 campaign, our team collected BRDF data at multiple sites within the dune system (Fig. 1).Due to time limitations, our measurements were collected from six locations using the goniometer of the Rochester Institute of Technology (GRIT). 5,6The GRIT positioned the fiber optic of an ASD spectrometer with 5 deg fore-optic over a set of hemispherical positions between 0 deg and 65 deg in zenith and 0 deg and 360 deg in azimuth.
Although the aperture size of the sensor is relatively small, measurements such as these are sometimes referred to as hemispherical conical reflectance factor (HCRF) measurements. 7,8he use of the term conical reflects the finite size of the sensing aperture, and "hemispherical" indicates that illumination comes from both direct and diffuse sources in an outdoor setting.
Our ASD FR4 spectrometer measurements onboard the GRIT used a 5-deg fore-optic.Due to the limited number of days available for measurement, we also collected sediment samples at a number of sites within the dune system to capture the typical variation in geophysical parameters.These samples were brought back to the lab where additional testing, described below, was performed to manipulate the geophysical properties of the sample to characterize the potential impact on the HCRF and to compare with radiative transfer models.Relevant parameters include the grain size distribution, sediment density, sediment composition, moisture content, surface roughness, and local topography.Radiative transfer models depend on parameters such as these, and, for example, the fill factor of the sediment (closely related to density) directly parametrizes approximate solutions to the radiative transfer equation (RTE), 4 which we will review in greater detail below.
Field spectral measurements were accompanied by a set of geotechnical measurements, which included sediment density, grain size distribution, and moisture content, as well as mechanical properties of the sediment, such as shear strength.For example, sediment density was measured with a sand cone density apparatus (Fig. 1).Samples collected on the dune were also analyzed to determine moisture content and grain size distributions using, respectively, a drying oven and a mechanical sieve shaker, both brought to the site in our mobile laboratory.To determine grain size distributions, sediments were dried in the oven prior to sorting into size fractions with the mechanical sieve shaker.At sites where GRIT HCRF measurements were made, all geotechnical measurements were completed.Additional samples were collected from other areas to provide examples of sediment variation found across the dune system.These were returned to our laboratory for use in laboratory-based experiments using GRIT (Fig. 1) and later a second-generation hyperspectral goniometer system known as GRIT-two (GRIT-T), shown in Fig. 5. 9 In these experiments, individual variables such as grain size distribution and density were manipulated, along with the illumination geometry, and used in a series of laboratory biconical reflectance factor (BCRF) 7,8,10 measurements.
The primary sites where field HCRF measurements were taken with GRIT in tandem with geotechnical measurements were on the western side of the Algodones Dunes system (Fig. 1).Additional sites where samples were collected for later laboratory use were on the northern edge, central-northern region, and central western edge of the dune system, as shown in Fig. 1.In each case, efforts were made to collect samples from the top, bottom, and middle of individual dunes.Grain size and density variations were observed across sites and even across positions within an individual dune.The next section describes the geotechnical measurements conducted in tandem with the GRIT hyperspectral HCRF measurements.

Geotechnical Measurements in the Algodones Dunes
Geophysical properties vary considerably across the dune system.At the sites where we measured HCRF, as well as at those sites where samples were collected for later laboratory analysis, grain size distributions show some significant differences.The average size factor, X ¼ πD λ , where D represents the average particle diameter and λ is the wavelength of light of an image band, varied by as much as a factor of four across all collected samples, while density varied by as much as 43% (Fig. 2).These density measurements were conducted in the field with a sand cone apparatus according to an American Society for Testing and Materials (ASTM) standard. 11Grain size distributions were determined onsite by drying the sediments in a drying oven over a 24-h period, and then sieving the sediment in a Humboldt mechanical sieve shaker for a 30-min period.Moisture content of the sediment was determined using the weights of the sample before and after drying.

Laboratory Measurements of the Angular Dependence of Hyperspectral
Reflectance: Controlling Geophysical Parameters In a laboratory setting, we conducted a series of BRDF measurements using sediments brought back from the Algodones Dunes site.Since grain size distribution and density both play an important role in the observed BRDF, measurements focused on manipulating these individual parameters.In addition, the illumination geometry was also changed to provide a wide variety of overall measurement scenarios.
To achieve samples of varying density, samples were pluviated using an apparatus in our laboratory following an ASTM standard (Fig. 3). 11At the top of the pluviation apparatus, a canister with a trapdoor releases sediment, which falls through a fixed set of coarse sieves placed at a specific drop height above the sample holder.The process is designed to mimic Aeolian deposition.In many of the experiments, we prepared samples over a wide range of densities obtained by varying the drop height from the coarse sieves to the sample holder in the pluviation device.
For the laboratory studies, BRDF measurements were conducted using two different generations of the goniometer system, GRIT 5,6 and GRIT-T. 9In a laboratory setting, BRDF measurements are sometimes referred to as BCRF 7,8,10 measurements when there is a finite aperture for both the single directional illumination source and the sensor.The term "factor" signifies the referencing of sample radiance data to the radiance of a Lambertian standard.In our experiments, Spectralon ® panels, which approximate a Lambertian surface, 12 were used as the reference standard.These panels were calibrated by LabSphere.Figure 4 shows some examples of BCRF of Algodones Dunes sediment samples for varying density.

Radiative Transfer Models for Retrieval of Sediment Geophysical
Properties: Fill Factor 8][19][20][21][22][23][24][25][26] The isotropic multiple scattering approximation (IMSA) is an approximate solution to an RTE based on the Chandresekhar-Ambartsumian method of invariance. 4 In the Hapke treatment of the RTE, five orders of scattering are considered, and the single scattering contribution through the single scattering phase function is treated in an exact form; however, for multiple scattering contributions, the phase function is assumed to be isotropic. 4The resulting approximate solution (IMSA) has the form E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 0 1 ; 1 1 6 ; 3 7 8 rðθ i ; θ e ; gÞ ¼ K wðλÞ 4 The factors μ i ¼ cos θ i and μ e ¼ cos θ e are the direction cosines for the incident and scattered zenith angles θ i and θ e ; pðg; λÞ is the single particle scattering phase function; g is the phase angle (angle defined by light source, ground, and sensor); K is the "porosity factor;" wðλÞ is the single scattering albedo; and Hðμ i ∕KÞ and Hðμ e ∕KÞ are the incident and view angle Chandrasekhar-Ambartsumian H-functions, which have an approximate solution of the form 4,15 E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 0 2 ; 1 1 6 ; 2 4 4 where 1−wðλÞ p is the diffuse hemispherical reflectance.The term B S0 B S ðg; K; λÞ is a correction to the single-scattering term and is a description of the influence of the shadow hiding opposition effect (SHOE), an increase in brightness observed at smaller phase angles (usually ≲20 deg). 4,14In this expression, B S0 is a constant that must be determined in an inversion procedure.The functional form of the opposition effect B S ðg; K; λÞ is discussed further below.The term B C0 B C ðg; K; λÞ describes the coherent backscatter opposition effect (CBOE), which is observed over only a very narrow range of phase angles, typically ≲2 deg, where coherent addition of scattering pathways in equal and opposite directions can occur. 4,15B C0 is another free constant that must in principle be determined, while B C ðg; K; λÞ describes its functional form, which, similarly to the SHOE, is large at small phase angles and decreases with increasing phase angle.In this paper, we neglect the CBOE since our measurements do not cover phase angles this small.Note that in Eq. ( 1), factors of μ i and 1∕π are removed since our measurements involve reflectance factors (BCRF and HCRF), which are the ratio of radiance values obtained from the sediment sample and a Spectralon ® white reference, an approximately Lambertian standard. 4,7,8,10,12mportantly, a geophysical parameter, the sediment fill factor ϕ (the fraction of volume occupied by particles), explicitly parameterizes the approximate IMSA solution since the nonlinear porosity factor K ¼ KðϕÞ is a function of the fill factor ϕ, and, for equant particles, it has the form 4,16 E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 0 3 ; 1 1 6 ; 6 1 6 Thus, inversion of the IMSA solution provides an explicit estimate of an important geophysical parameter, the sediment fill factor.In this work, our goal was to obtain a relative fill factor through inversion of the radiative transfer model.We wanted to demonstrate that the resulting retrieved fill factor correlates consistently with the directly measured sediment density measured in the laboratory, even if we are not sure of its absolute magnitude.
A second geophysical property that is also implicitly present in the IMSA solution is the sediment grain size distribution.Within the IMSA solution, this plays a role in the term related to the SHOE.Hapke has derived an approximate form for the SHOE, which takes the form 4,14 E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 0 4 ; 1 1 6 ; 4 7 3 For a number of explicit examples, Hapke has shown that the width parameter of the SHOE is proportional to the product of KðϕÞϕ, with the constant of proportionality depending on the shape and extent of the grain size distribution. 4One such example is a unimodal case, having the form re − r∕r , where r is the particle radius and r is the average radius.This grain size distribution leads to a SHOE width parameter, h S ¼ ð 3 8 Þ 3 2 KðϕÞϕ. 4,27Grain size distributions have been the focus of earlier studies, one of which included a retrieval of a weighted average grain size from a desert site. 28he actual shape of distributions found in nature may be more complicated.Figure 2 shows several examples of the grain size distributions derived from our sediment samples acquired in the Algodones Dunes.Several of these show grain size distributions with more than one mode.In many of the examples shown in the figure, the secondary mode is a smaller peak of fine sediments, such as fine sand and/or silt.However, there is at least one example from the sediment samples that has a complicated set of multiple peaks (sample 1202-M-02-B).For sample 1202-M-02-B, the full distribution extends beyond the range shown in Fig. 2, and this is addressed later in our discussion of laboratory experiments.
While a significant portion of the distributions in the samples was unimodal, many were not, and since multimodal distributions are quite common, we have modeled the width of the SHOE peak in this study as being proportional to KðϕÞϕ, with the proportionality constant being a free parameter ϵ that must be optimized in our inversion of the IMSA solution E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 0 5 ; 1 1 6 ; 1 8 9 h S ¼ ϵKðϕÞϕ: (5) As each free parameter is introduced into an inversion procedure, the optimization process grows progressively more complicated.Thus, adding the free parameter ϵ, while logical from a modeling perspective, adds to the complexity of the model residual surface being searched, especially for a model as nonlinear as IMSA.In order to avoid this pitfall, we have taken a different approach, building upon an earlier workflow that we previously developed. 6We have added some important refinements, which allow the approach to be more successful.Another variation of this approach, which we will introduce later, takes advantage of the explicit form of the SHOE width found in Eq. ( 5).
In that initial exploration, 6 we inverted the Hapke IMSA solution solving Eq. ( 1) for the product of the single particle phase function and the SHOE correction E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 0 6 ; 1 1 6 ; 7 This was done for two different BCRF scans of the same sample sediment surface at two different illumination geometries (Fig. 5).The approach then involved a search over the two-parameter space of the single scattering albedo wðλÞ and the fill factor ϕ.This was done in such a manner that the quantity pðg; λÞ½1 þ B S0 B S ðg; K; λÞ was an invariant function.This function describes the underlying single particle characteristics, which was assumed to be identical between the two illumination conditions.This was ensured by minimizing Here, the numerical subscripts refer to the two illumination conditions in which the source zenith angle is varied in the two sets of BCRF measurements.We have found that regularization is necessary to achieve stable solutions numerically, so that in practice we minimize One subtlety of the minimization approach, described in Eqs. ( 7) and ( 8), is that the set of sample phase angles in the two illumination conditions is not guaranteed to be the same, and we elaborate here on some of the technical details required for this to successfully converge to a meaningful solution.An overall flow diagram of the processing steps is shown in Fig. 6.
In the inversion process, the two sets of BCRF scans must be interpolated to a common phase angle grid to facilitate the difference calculation specified in Eq. ( 8).Although we tried a number of different interpolation schemes, we found that Steffen's algorithm 29 provided a stable, wellbehaved interpolation of the BCRF data.Additional refinements, which decrease the overall residual, include the use of a low-pass filter (Sovitsky-Golay filter) 30 during the interpolation onto the common phase angle grid.In the optimization procedure, we typically used a Nelder-Mead simplex method, 31 although we also evaluated the Levenberg-Marquardt algorithm, 32 which did not produce appreciably different results.
A new refinement, which we introduce in this paper, involves an additional constraint that improves the convergence of the model to a consistent solution.This constraint is motivated by the goal that the resulting fill factor retrieved should be consistent across wavelengths.This is, of course, a physical requirement of our approach, but without insisting on this constraint, the processing steps described in Eqs. ( 6) and ( 8) and illustrated in Fig. 5 are essentially wavelength agnostic.The system of equations can be applied independently to each wavelength, and therefore there is a potential to obtain solutions that do not report a consistent fill factor ϕ across all wavelengths.In contrast, we expect to obtain a unique value of wðλÞ for each wavelength λ.
To ensure this constraint, the Nelder-Mead optimization step is set up with an additional constraint.We want to constrain the individual estimates obtained for the fill factor at a particular wavelength ϕ λ so that they do not differ from the average by more than a specified amount δ E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 1 0 ; 1 1 6 ; 3 9 5 To implement this constraint on the optimization procedure, we must first run the optimization over all wavelengths λ [Eqs.( 8) and ( 6)] without the constraint found in Eq. (10).We then repeat the same steps, but this time the Nelder-Mead simplex optimization step enforces the constraint of Eq. (10).We can reiterate this constraint if desired and let the size of δ be given by δ ¼ δðtÞ, where δðtÞ is a decreasing function of the iteration step t.
We also explored another variation of this inversion process based on an earlier proposed modification to the Hapke model, in which we reintroduced directional dependence into the Fig. 6 Workflow for inversion of the IMSA model using the BCRF scans from two different illumination geometries based on Eqs. ( 6), (8), and (10).
multiple scattering term. 5The inclusion of this factor was not derived in a rigorous manner; however, by insisting that the single particle phase function and SHOE correction be inserted in front of the multiple scattering term, this approach can be thought of ensuring that the final scattering event in a multiple scattering pathway has the directional dependence of a single particle, while the other scattering events are assumed to average out and not contribute to the directional dependence.The form used in the earlier work was 5 E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 1 1 ; 1 In the original version of this modified approach, the factor fpðg; λÞ½1 þ B S0 B S ðg; K; λÞg est was first estimated self-consistently as the difference in Eq. ( 6); 5 then the optimization of Eq. ( 8) proceeded with fpðg; λÞ½1 þ B S0 B S ðg; K; λÞg est being used on the right-hand side of Eq. (11).In this paper, however, we use a more refined version of this approach that considers some of the underlying parameters that determine the form of pðg; λÞ½1 þ B S0 B S ðg; K; λÞ explicitly in separate optimization steps.fpðg; λÞ½1 þ B S0 B S ðg; K; λÞg est in Eq. ( 11) is also replaced by ηpðg; λÞ in the present version, which requires us to separate out the functional form of pðg; λÞ as part of our processing chain.The processing steps that we describe are outlined in Fig. 7.
The initial steps of this inversion process follow a similar approach outlined in Eqs. ( 6) and ( 8).This provides us with an initial estimate of the fill factor ϕðt ¼ 0Þ and single scattering albedo wðλ; t ¼ 0Þ.These values are then passed through Eq. ( 6) to obtain a first estimate for both illumination conditions: fpðg; λÞ½1 þ B S0 B S ðg; K; λÞg est:;1;2 .Here the subscripts "1" and "2" refer to the two illumination conditions.These estimates are interpolated onto the common phase angle grid during each step of the initial optimization as before, and then with the initial estimates obtained, ϕðt ¼ 0Þ and wðλ; t ¼ 0Þ, a separate optimization step is used to find an explicit parameterization of fpðg; λÞ½1 þ B S0 B S ðg; K; λÞg est:;1;2 in terms of the free parameters B S0 , ϵ [Eq.( 5)], and the coefficients of the three-parameter Henyey-Greenstein function b 1 ; b 2 , and c 4 E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 1 2 ; 1 1 6 ; 3 6 2 p HG ðg; λÞ ¼ The goal of this step is ultimately to find an estimate of pðg; λÞ as a separate function.In this optimization step, we again use the Nelder-Mead simplex method, where the residual between the forward-propagated value of fp HG ðg; λÞ½1 þ B S0 B S ðg; K; λÞg fwd using the parameters B S0 , Fig. 7 Workflow for the modified model incorporating directional dependence in the multiple scattering term in a multistage optimization procedure.
ϵ, b 1 ; b 2 , and c is compared against the first estimates of the single scattering properties fpðg; λÞ½1 þ B S0 B S ðg; K; λÞg est:;1;2 for the two illumination conditions at the common phase angle grid obtained in the earlier optimization step.From the estimated values of b 1 ; b 2 , and c resulting from this step, we obtain an estimate of pðg; λÞ as a separate factor from the original estimates of fpðg; λÞ½1 þ B S0 B S ðg; K; λÞg est:;1;2 .We insist at this new step that the minimization be over a common set of coefficients B S0 , ϵ, b 1 ; b 2 , and c for both illumination conditions.Therefore, the optimization is ; t e m p : i n t r a l i n k -; e 0 1 3 ; 1 1 6 ; 6 5 1 Having now an estimate of pðg; λÞ from this step, we then return to an equation similar to Eq. ( 11), but modified to include just the phase function and a free scale parameter η; our modified model using our estimate of pðg; λÞ is ; t e m p : i n t r a l i n k -; e 0 1 4 ; 1 1 6 ; 5 5 0 We then calculate Eq. ( 14) for both illumination conditions, interpolating the result to a common grid, followed by optimization of Eq. ( 8) using the Nelder-Mead method.This time, however, the optimization is over the parameters ϕ, wðλÞ, and η.This optimization step is done once first without the constraint of Eq. ( 10) to estimate φ.Subsequent iterations of the estimation procedure then enforce the constraint of Eq. ( 10) during the optimization.

Field HCRF and Laboratory BCRF Measurements: Impact of Geophysical Parameters
Above, we noted that density and grain size distributions can both impact measured BCRF and HCRF data.In previous studies, we have observed strong density-dependent effects, and the trend for these effects has been observed to depend somewhat on the relative optical contrast of the underlying mixture of materials in the granular sediment in combination with the illumination zenith angle, which determines the degree of multiple scatter within the medium, and thus the potential for darker optically contrasting minerals to diminish multiple scatter. 5,6,33Grain size effects also play an important role, since they are not entirely independent from another free variable that affects BCRF and HCRF, namely surface roughness.Although not modeled in the present work, Hapke has also developed a correction to the IMSA model for macroscopic surface roughness. 4,34These corrections are applicable to the majority of the grain sizes found in our distributions, which all fall within the geometric optics region, with the exception of very fine clay particles, which typically have size ≲2 μm.In the laboratory studies in which grain size was explicitly manipulated as an independent variable, this directly impacted the apparent surface roughness.Algodones sediment sample 1102-MS-02 from the north-central region of the dune system (Fig. 1) was prepared in the laboratory, where different subsets of the grain size distributions were separated through mechanical sieve shaking (Fig. 3).Sieved fractions were included or removed, and these new samples with manipulated grain size distributions (Fig. 8) formed the basis of several experiments in which density was also manipulated through pluviation (Fig. 3).Grain size manipulations resulted in the apparent differences in surface roughness (Fig. 8).
In general, altering the grain size distribution to include significant gaps (case 3 of Fig. 8) or measuring from the original distribution that already had significant gaps and therefore surface roughness tended to produce results showing the greatest change in the angular distribution when different illumination zenith angles were chosen.We have observed trends like this in the past.Samples with the most pronounced roughness (case 3 in Fig. 8) tend to exhibit the greatest backscatter when the light source is more oblique (the 65-deg illumination case, Fig. 9, top row, last column), similar to what we have observed in previous studies. 27When illumination angles are close to nadir, the amount of diffuse multiple scatter often significantly increases (second row of Fig. 9), as observed in previous studies. 5ensity of granular sediments also plays an important role in observed BCRF and HCRF. 4,5,33In general, reflectance tends to increase with density within certain conditions and only up to a maximum density.In our past studies, we have observed that the degree of Fig. 9 GRIT BCRF measurements using the three sample preparations of grain size distributions described in Fig. 8   multiple scatter present and the degree of optical contrast of constituents both play a role in the observed trend. 33The trend of increasing reflectance with density tends to be for cases when there is less multiple scatter, as may be expected often at more oblique illumination angles. 5At angles close to nadir, where greater multiple scatter may be present, if the sediment mixture contains fractions that are smaller and darker minerals, then increasing the density may lead to a decrease in overall reflectance 5 since densifying the material will tend to force the smaller, darker particles to more optimally occupy the available pore space, shutting down multiple scattering pathways. 33When illumination is closer to nadir and therefore more likely to produce diffuse multiple scatter, this leads to a decrease in reflectance. 5At oblique angles with less multiple scatter, if the trend does follow the predictions of the Hapke model, then the reflectance should increase up to a critical maximum at which point particles are so close together that they begin to behave like larger single particles, and reflectance tends to decrease. 4We have observed this trend in the past under these types of conditions. 5The overall variance between the spectra may vary with illumination angle and thus degree of multiple scatter.In Fig. 9, we see greater diffuse multiple scatter when the illumination zenith angle is closer to nadir at 25 deg compared to the oblique case of 65 deg, and in Fig. 10, we see that the spectral plots of the full set of HCRF measurements for the 20-deg zenith angle illumination also have greater overall spectral variance compared to the more oblique case at 40 deg.
The relative amount of forward to backward scatter may also change with density, 5 and we see this to some extent in Fig. 10, row 3, where the extent of the backward scatter and relative amount of the backscatter compared to the amount of forward scatter changes with density for the 450-nm wavelength.At the longer SWIR wavelength (fourth row of Fig. 10), the main difference with increasing density is a more gradual change in the extent of the backscatter region (especially along the azimuthal dimension), which grows larger with increasing density.These effects with density may be stronger when the illumination zenith angle is even larger, as we have observed in earlier work. 5

Matching Field and Laboratory HCRF and BCRF to G-LiHT Hyperspectral Imagery
A spectral library was created from all the field hyperspectral HCRF measurements recorded with GRIT during the March 2015 field experiment at Algodones as well as the hyperspectral BCRF (using GRIT and GRIT-T) of samples brought back from the dunes and prepared as described earlier.The total size of the spectral library generated for this study included 13,519 reflectance spectra.Although to date, the laboratory experiments have not exhausted all the collected sediment samples, the laboratory BCRF makes an important contribution to the representation of the spectral characteristics of the dune system (Figs.11 and 12).The figures portray a Euclidean spectral match between the spectral library and the G-LiHT imagery after the library has been resampled to the spectral grid of the G-LiHT imagery.The figures show that manipulating geophysical parameters of field samples brought back to the lab can obtain better matches over significant parts of the imagery acquired with G-LiHT during the March 2015 field campaign.In the sequences of laboratory experiments, we systematically explored varying the sediment density, grain size distribution, and illumination zenith angle of the source relative to the sediment surface.All of these libraries match parts of the G-LiHT imagery better than field measurements in at least some part of the image subsets shown in these figures, and especially in Fig. 12, where roughly half of the scene is better matched by the laboratory-manipulated spectral BCRF libraries.Note that in Fig. 12, the manipulation of grain size to prepare three very different types of distributions followed by pluviation to produce different densities had a highly significant role in the spectral matching shown in Fig. 12.Likewise, density effects are highlighted at the top of the dune ridges in this scene subset.This is consistent with what we observed in our field geotechnical data, which showed differences in density between the bottom and top of a dune, for example, between sites BC01-01 and BC01-02 (Fig. 2).Note that in both cases in Figs.11 and 12, none of the GRIT field measurement locations were within the boundaries of the scene.These results are significant because it means that laboratory manipulations of field samples can be used to successfully augment the representative spectral library for the field site.Since the Algodones Dunes is envisioned as an ongoing vicarious calibration site, this means that scene simulation packages such as DIRSIG will have more reliable spectra as input.

Radiative Transfer Model Inversion Using Laboratory Biconical Reflectance Factor Measurements
We used the optimization procedure for inverting the original Hapke model described in Sec.2.4.
For this particular analysis, a set of five densities was created for the sample from site 0903-T-B (Fig. 1) using the air pluviation technique previously described (Fig. 3).These densities ranged from 1.59 to 1.82 g∕cm 2 .The retrieved fill factor was then compared against the measured density of the material in the sample holder.The regression between retrieved fill factor and measured density is shown in Fig. 13 for each of the two approaches that we described in Sec.2.4.We found that the first iteration of the constraint equation provided the best-correlated estimate with direct laboratory measurements of sediment density.At this stage of processing, the VNIR and SWIR parts of the spectrum give slightly different estimates of the fill factor, and it is clear that we should interpret the results only as a measure of the relative fill factor rather than an absolute measure of this parameter.Correlations of the retrieved fill factor with directly measured density were comparable for the VNIR, with R 2 VNIR ¼ 0.73 for the inversion of the original IMSA model and R 2 VNIR ¼ 0.71 for the modified model, which includes an angular dependence in multiple scattering.In the SWIR, the estimates were less correlated with directly measured density, with R 2 SWIR ¼ 0.50 for inversion of the original IMSA model and R 2 SWIR ¼ 0.60 for the modified model.The weaker correlation in the SWIR is likely due to lower SNR in this spectral region, although that could be improved by recording a larger number of spectral samples at each position.However, it is clear that the modified model was better correlated in this region of the spectrum.In both cases, the estimates of the single scattering albedo (Fig. 14) appear consistent with each other for all wavelengths and across all densities.In only one case, for the density of 1.79 g∕cm 2 , the solution obtained by the original IMSA model exhibited more noise (a spike in the spectral shape ∼700 nm) in the VNIR, which is not observed in the modified model.

Conclusions
We have seen that field-measured HCRF with the GRIT hyperspectral goniometer system is well matched to the hyperspectral imagery collected during the Algodones Dunes campaign.As additional samples brought back to the laboratory have been studied and manipulated for BCRF studies with GRIT and GRIT-T, we have seen that they provide a better match to different parts of the dune system where GRIT measurements were not obtained while on site.It is an important point that laboratory measurements, which manipulate the geophysical parameters, such as density and grain size distribution, can provide important data that augment field BRDF spectral libraries to make the library more representative of the dune system in general.Because it is likely that NASA will use the Algodones Dunes as an ongoing vicarious calibration site for satellite systems, these observations are particularly important.This means that more representative data can be created to augment costly field campaigns and be useful in scene simulation/ prediction packages such as DIRSIG, which are important for calibration and rely on representative BRDF as a foundational input to the overall processing chain.In addition, the data collected during this campaign have also served the dual purpose of allowing us to develop and test radiative transfer models that are more sophisticated than those being used, such as Ross-Li, in packages, such as DIRSIG.
Beyond the significance for calibration, radiative transfer models like those developed by Hapke (and variants described here) are important for geophysical parameter inversion from imagery.In this work, we used the samples brought back from Algodones to develop and test the idea of using BCRF measurements from two illumination geometries to invert an important geophysical parameter of the radiative transfer model, namely the sediment fill factor, which is ubiquitous in the approximate solution to the radiative transfer equation developed by Hapke.The laboratory-based results that were used in the inversion allowed us to demonstrate an inversion of the fill factor that strongly correlates with measured sediment density, an observation that suggests the validity of the inversion process.This does not confirm the absolute scale of the fill factor obtained; however, it does confirm the expected trend of relative fill factor with respect to directly measured density.In the future, a source of potential improvement to our approach would be to include a more sophisticated constraint in the spectral domain, as even in imposing the relatively weak constraint across wavelength used in this paper, convergence of the inversion and consistency across wavelengths was greatly improved.
To adapt this approach that we have taken for geophysical parameter inversion to a practical strategy where this can be accomplished for imagery, such as G-LiHT, there are still some technical hurdles to overcome, but the initial results are promising.The additional hurdles include the fact that the number of relative geometries over which the surface is both illuminated and viewed for a given position on the ground will depend on factors, for example, such as satellite revisit schedule or length of campaign and design of imaging flight patterns and times of day for airborne studies.The success of the approach that we described depended on the ability to interpolate an estimate of the single scattering phase function over the same phase angles with SHOE correction and a comparison of a densely sampled set of views over a hemisphere for two different illumination cases.In an imagery time series, the combination of illumination direction and view direction in a time series may be more complex; however, one point is clear: the better represented the range of phase angles is in the imagery time series, the more likely the interpolation and comparison can be successful.In the G-LiHT imagery flown over Algodones, sets of lines were flown along a particular heading, with successive lines having a significant percentage of overlap.Since the lines were short, this meant that at least several looks at the same point on the ground occurred from multiple view directions with roughly similar illumination conditions.On successive days, other orientations were flown, providing a diversity of potential geometries, and AM and PM lines were flown on each flight day so that the illumination conditions also varied.Thus, comparison of pixels from the same point on the ground among the various flight lines could be used to compare and interpolate the single scattering properties needed in our approach for two or more different illumination geometries.Likewise, in a satellite time series with frequent enough revisit time to ensure sites are sufficiently stable during the acquisition of the time series, a similar strategy might be followed.The essential point of interpolating an estimate of the single scattering properties from observations and modeled multiple scattering appears to be feasible.

Fig. 1
Fig. 1 GRIT during measurements at the Algodones Dunes experiment: (a) sediment HCRF measurement at a dune site with main components of the instrument highlighted (ASD spectrometer, two Vectornav GPS/IMU units, Ximea FoV Camera; (b) measurement of a Spectralon ® white reference during overflight by the NASA G-LiHT.(c) Google Earth view of the Algodones Dunes system showing the location of RIT GRIT HCRF validation sites and ancillary sites where additional samples were collected for laboratory studies.(d) Geotechnical measurements included onsite density measurements, shown here using a sand cone density apparatus.

Fig. 2
Fig. 2 Variations in grain size distribution in the top inch of the sediment at (a) the six field sites where GRIT measurements were undertaken at the Algodones Dunes site and (b) at sample sites in the central portion of the dunes and (c) central western edge of the dunes.(d) Density variations were also considerable at the six field sites where HCRF measurements were undertaken.

Fig. 3
Fig. 3 Sample preparation apparatuses.(a) Air pluviation device: at top, a canister with trapdoor releases sample through coarse sieves with sample holder beneath; height of sieves above the holder determines resulting density of sediment.(b) Humboldt mechanical sieve shaker used to separate grain sizes; (c) Humboldt drying oven used to dry sediments.

Fig. 4
Fig.4BCRF data obtained using GRIT-T for varying density of the same sample: spectral reflectance curves measured in the laboratory for varying density of sediment samples from the central dune system in Algodones.Color code corresponds to the azimuth of the recorded sample.

E 2 þ αwðλÞ 2 þ βKðϕÞ 2 : ( 8 )Fig. 5 (
Fig. 5 (a and b) Paradigm used to collect laboratory BCRF data with GRIT and GRIT-T for inversion of radiative transfer models.(a) Large illumination zenith angle case, and (b) small illumination zenith angle case, which generally will have a greater degree of multiple scatter.Laboratory hyperspectral goniometric measurements using (c) GRIT and (d) second-generation GRIT-T.Both instruments were used to develop the spectral libraries matched to G-LiHT imagery.
Fig.9GRIT BCRF measurements using the three sample preparations of grain size distributions described in Fig.8for an example wavelength in the SWIR (1200 nm) from the onboard ASD spectrometer.Illumination zenith angle of 65 deg (a-c) and 25 deg (d-f).(a and d) Original sample grain size distribution; (b and e) intermediate particle sizes only; (c and f) largest and smallest particles.Light source azimuth was 336 deg.GRIT self-shaded region is removed.

Fig. 8
Fig. 8 Algodones sample 1102-MS02-2 used in BCRF measurements: (a) original sample; (b) grains in central peak only; this peak is the part of the distribution of 1102-MS02-2 shown in Fig. 2. (c) Largest and smallest grains with central peak of distribution removed.Gaps in the original distribution and the case with largest and smallest particles exhibit significantly greater roughness, as expected.

Fig. 10 (
Fig. 10 (a) Spectral BCRF measured by GRIT-T for two different illumination zenith angles: 20 deg (a) and 40 deg (b) for a density of 1.59 g∕cm 2 , and a comparison of the angular distribution of the reflectance at selected wavelengths throughout the spectrum.(c) Comparison of the BCRF for selected wavelengths from the ASD spectrometer at 450 nm as a function of increasing density (from left to right); (d) The same sequence for a SWIR wavelength (2250 nm).(e) Picture of the sample (0903-T-02) at the highest density (1.82 g∕cm 2 ); (f) measured grain size distribution.GRIT-T has a relatively small self-shaded region, and samples in these experiments were every 15 deg in azimuth and 10 deg in zenith.Plotted results are interpolated between these points.Illumination azimuth angle is 7.5 deg.

Fig. 11
Fig. 11 Spectral matching of HCRF and BCRF from field and laboratory settings.(a) Subset of one of the G-LiHT "A"-lines flown on March 13, 2015.(b) Detailed spectral matching with the combined field and laboratory spectral HCRF and BCRF library containing 13,519 spectral reflectance elements color-coded to individual elements of the library.(c) Spectral Euclidean distance between the best-matched library spectrum and the G-LiHT spectral pixel.This region of the dunes contained no GRIT field measurements.(d) Breakdown of library elements by source (field HCRF, lab BCRF with various types of manipulation, such as grain size and density, as well as illumination zenith angle).

Fig. 12
Fig. 12 Spectral matching of HCRF and BCRF from field and laboratory settings for (a) a subset of one of the large area flight lines flown by G-LiHT on March 10, 2015; (b) closest spectral matching library element; (c) spectral Euclidean distance to best-matching library element (114 G-LiHT spectral band space); and (d) breakdown of types of categories providing best spectral match:note that for this image subset, GRIT lab BCRF measurements involving grain size, density, and illumination variations for a site where no GRIT field measurements were taken (yellow) and GRIT lab BCRF involving density and illumination variations for samples taken at GRIT field sites (red) provided the best spectral match in a large fraction of the scene.

Fig. 13
Fig. 13 Average fill factor (average over wavelength) versus directly measured density in the laboratory for VNIR and SWIR regions.(Left) For the original IMSA model using our inversion approach, R 2 VNIR ¼ 0.73 and R 2 SWIR ¼ 0.50.For our modified model, which includes an angular dependence in multiple scattering, R 2 VNIR ¼ 0.71 and R 2 SWIR ¼ 0.60.Results are shown after the first constrained iteration.

Fig. 14
Fig. 14 The corresponding inversion of the single scattering albedo w ðλÞ for our inversion of (a) the original IMSA model and (b) the modified model, which incorporates directional dependence in the multiple scattering term.