27 September 2017 Modeling and intercomparison of field and laboratory hyperspectral goniometer measurements with G-LiHT imagery of the Algodones Dunes
Author Affiliations +
J. of Applied Remote Sensing, 12(1), 012005 (2017). doi:10.1117/1.JRS.12.012005
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.
Bachmann, Eon, Ambeau, Harms, Badura, and Griffo: Modeling and intercomparison of field and laboratory hyperspectral goniometer measurements with G-LiHT imagery of the Algodones Dunes



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–Li3 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.4 In 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,6 The 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,8 The 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.

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.


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.11 Grain 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.

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.



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).11 At 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.

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.


For the laboratory studies, BRDF measurements were conducted using two different generations of the goniometer system, GRIT5,6 and GRIT-T.9 In a laboratory setting, BRDF measurements are sometimes referred to as BCRF7,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.

Fig. 4

BCRF 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.



Radiative Transfer Models for Retrieval of Sediment Geophysical Properties: Fill Factor

Absorption and scattering from sediment surfaces have been modeled with radiative transfer models in both planetary astronomy and Earth remote sensing applications,4,1314.15.16 and Hapke’s approximate solution to the RTEs has been widely used.1718. 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.4 The resulting approximate solution (IMSA) has the form


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 form4,15


where r0=1γ1+γ=11w(λ)1+1w(λ) is the diffuse hemispherical reflectance. The term BS0BS(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,14 In this expression, BS0 is a constant that must be determined in an inversion procedure. The functional form of the opposition effect BS(g,K,λ) is discussed further below. The term BC0BC(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,15 BC0 is another free constant that must in principle be determined, while BC(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,12

Importantly, 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 form4,16



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 form4,14


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.4 One such example is a unimodal case, having the form rer/r¯, where r is the particle radius and r¯ is the average radius. This grain size distribution leads to a SHOE width parameter, hS  =(38)32K(ϕ)ϕ.4,27 Grain 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.28

The 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


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.6 We 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


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+BS0BS(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



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.


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.

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).


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 algorithm29 provided a stable, well-behaved 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 δ


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 multiple scattering term.5 The 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 was5


In the original version of this modified approach, the factor {p(g,λ)[1+BS0BS(g,K,λ)]}est was first estimated self-consistently as the difference in Eq. (6);5 then the optimization of Eq. (8) proceeded with {p(g,λ)[1+BS0BS(g,K,λ)]}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+BS0BS(g,K,λ)] explicitly in separate optimization steps. {p(g,λ)[1+BS0BS(g,K,λ)]}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.

Fig. 7

Workflow for the modified model incorporating directional dependence in the multiple scattering term in a multistage optimization procedure.


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: {p(g,λ)[1+BS0BS(g,K,λ)]}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 {p(g,λ)[1+BS0BS(g,K,λ)]}est.,1,2 in terms of the free parameters BS0, ε [Eq. (5)], and the coefficients of the three-parameter Henyey–Greenstein function b1,b2, and c4


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 {pHG(g,λ)[1+BS0BS(g,K,λ)]}fwd using the parameters BS0, ε, b1,b2, and c is compared against the first estimates of the single scattering properties {p(g,λ)[1+BS0BS(g,K,λ)]}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 b1,b2, and c resulting from this step, we obtain an estimate of p(g,λ) as a separate factor from the original estimates of {p(g,λ)[1+BS0BS(g,K,λ)]}est.,1,2. We insist at this new step that the minimization be over a common set of coefficients BS0, ε, b1,b2, and c for both illumination conditions. Therefore, the optimization is


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


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,33 Grain 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,34 These 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).

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.


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.27 When 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.5

Fig. 9

GRIT BCRF measurements using the three sample preparations of grain size distributions described in Fig. 8 for 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.


Density of granular sediments also plays an important role in observed BCRF and HCRF.4,5,33 In 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 multiple scatter present and the degree of optical contrast of constituents both play a role in the observed trend.33 The 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.5 At 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 reflectance5 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.33 When illumination is closer to nadir and therefore more likely to produce diffuse multiple scatter, this leads to a decrease in reflectance.5 At 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.4 We have observed this trend in the past under these types of conditions.5 The 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.

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/cm2, 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/cm2); (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.


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.

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

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.


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/cm2. 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.

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, RVNIR2=0.73 and RSWIR2=0.50. For our modified model, which includes an angular dependence in multiple scattering, RVNIR2=0.71 and RSWIR2=0.60. Results are shown after the first constrained iteration.


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 RVNIR2=0.73 for the inversion of the original IMSA model and RVNIR2=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 RSWIR2=0.50 for inversion of the original IMSA model and RSWIR2=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/cm2, 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.

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.




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.


1. B. D. Cook et al., “NASA Goddard’s LiDAR, hyperspectral and thermal (G-LiHT) airborne imager,” Remote Sens. 5(8), 4045–4066 (2013). http://dx.doi.org/10.3390/rs5084045 Google Scholar

2. A. A. Goodenough and S. D. Brown, “DIRSIG 5: core design and implementation,” Proc. SPIE 8390, 83900H (2012).PSISDG0277-786X http://dx.doi.org/10.1117/12.919321 Google Scholar

3. J.-L. Roujean, M. Leroy and P.-Y. Deschamps, “A bidirectional reflectance model of the Earth’s surface for the correction of remote sensing data,” J. Geophys. Res. 97(D18), 20455 (1992).JGREA20148-0227 http://dx.doi.org/10.1029/92JD01411 Google Scholar

4. B. Hapke, Theory of Reflectance and Emittance Spectroscopy, Cambridge University Press, Cambridge (2012). Google Scholar

5. C. M. Bachmann et al., “Improved modeling of multiple scattering in hyperspectral BRDF of coastal sediments observed using the goniometer of the Rochester Institute of Technology (GRIT),” Proc. SPIE 9611, 96110J (2015).PSISDG0277-786X http://dx.doi.org/10.1117/12.2188866 Google Scholar

6. C. M. Bachmann et al., “Modeling geophysical properties of the Algodones Dunes from field and laboratory hyperspectral goniometer measurements using GRIT and comparison with G-LiHT imagery,” Proc. SPIE 9972, 99720K (2016).PSISDG0277-786X http://dx.doi.org/10.1117/12.2238270 Google Scholar

7. F. Nicodemus, J. Richmond and J. Hsia, “Geometrical considerations and nomenclature for reflectance,” Sci. Technol. 60, 1–52 (1977). http://dx.doi.org/10.1109/LPT.2009.2020494 Google Scholar

8. G. Schaepman-Strub et al., “Reflectance quantities in optical remote sensing-definitions and case studies,” Remote Sens. Environ. 103(1), 27–42 (2006). http://dx.doi.org/10.1016/j.rse.2006.03.002 Google Scholar

9. J. D. Harms et al., “A next generation field-portable goniometer system,” Proc. SPIE 9840, 98400J (2016).PSISDG0277-786X http://dx.doi.org/10.1117/12.2223180 Google Scholar

10. E. J. Milton, N. P. Fox and M. E. Schaepman, “Progress in field spectroscopy,” in IEEE Int. Conf. on Geoscience and Remote Sensing Symp. (IGARSS 2006), Vol. 113, pp. 1966–1968 (2006). http://dx.doi.org/10.1109/IGARSS.2006.509 Google Scholar

11. D. Lo Presti, S. Pedroni and V. Crippa, “Maximum dry density of cohesionless soils by pluviation and by ASTM D 4253-83: a comparative study BT—maximum dry density of cohesionless soils by pluviation and by ASTM D 4253-83: a comparative study,” Geotech. Test. J. 15(2), 180–189 (1992).GTJODJ0149-6115 http://dx.doi.org/10.1520/GTJ10239J Google Scholar

12. D. Biliouris et al., “A compact laboratory spectro-goniometer (CLabSpeG) to assess the BRDF of materials. Presentation, calibration and implementation on Fagus sylvatica L. leaves,” Sensors 7(9), 1846–1870 (2007).SNSRES0746-9462 http://dx.doi.org/10.3390/s7091846 Google Scholar

13. B. Hapke and E. Wells, “Bidirectional reflectance spectroscopy: 2. Experiments and observations,” J. Geophys. Res. 86(B4), 3055–3060 (1981).JGREA20148-0227 http://dx.doi.org/10.1029/JB086iB04p03055 Google Scholar

14. B. Hapke, “Bidirectional reflectance spectroscopy 4. The extinction coefficient and the opposition effect,” Icarus 67, 264–280 (1986).ICRSA50019-1035 http://dx.doi.org/10.1016/0019-1035(86)90108-9 Google Scholar

15. B. Hapke, “Bidirectional reflectance spectroscopy 5. The coherent backscatter opposition effect and anisotropic scattering,” Icarus 157(2), 523–534 (2002).ICRSA50019-1035 http://dx.doi.org/10.1006/icar.2002.6853 Google Scholar

16. B. Hapke, “Bidirectional reflectance spectroscopy 6. Effects of porosity,” Icarus 195(2), 918–926 (2008).ICRSA50019-1035 http://dx.doi.org/10.1016/j.icarus.2008.01.003 Google Scholar

17. M. Ciarniello et al., “Hapke modeling of Rhea surface properties through Cassini-VIMS spectra,” Icarus 214(2), 541–555 (2011).ICRSA50019-1035 http://dx.doi.org/10.1016/j.icarus.2011.05.010 Google Scholar

18. M. K. Shepard and P. Helfenstein, “A test of the Hapke photometric model,” J. Geophys. Res. Planets 112(3), 1–17 (2007). http://dx.doi.org/10.1029/2005JE002625 Google Scholar

19. P. Helfenstein and M. K. Shepard, “Testing the Hapke photometric model: improved inversion and the porosity correction,” Icarus 215(1), 83–100 (2011).ICRSA50019-1035 http://dx.doi.org/10.1016/j.icarus.2011.07.002 Google Scholar

20. M. K. Shepard and P. Helfenstein, “A laboratory study of the bidirectional reflectance from particulate samples,” Icarus 215(2), 526–533 (2011).ICRSA50019-1035 http://dx.doi.org/10.1016/j.icarus.2011.07.033 Google Scholar

21. P. Helfenstein and J. Veverka, “Photometric properties of lunar terrains derived from Hapke’s equation,” Icarus 72(2), 342–357 (1987).ICRSA50019-1035 http://dx.doi.org/10.1016/0019-1035(87)90179-5 Google Scholar

22. G. J. Yang et al., “Extension of the Hapke bidirectional reflectance model to retrieve soil water content,” Hydrol. Earth Syst. Sci. 15(7), 2317–2326 (2011). http://dx.doi.org/10.5194/hess-15-2317-2011 Google Scholar

23. F. Schmidt and J. Fernando, “Realistic uncertainties on Hapke model parameters from photometric measurement,” Icarus 260, 73–93 (2015).ICRSA50019-1035 http://dx.doi.org/10.1016/j.icarus.2015.07.002 Google Scholar

24. D. Domingue and B. Hapke, “Fitting theoretical photometric functions to asteroid phase curves,” Icarus 78(2), 330–336 (1989).ICRSA50019-1035 http://dx.doi.org/10.1016/0019-1035(89)90181-4 Google Scholar

25. D. Domingue and A. Verbiscer, “Re-analysis of the solar phase curves of the icy Galilean satellites,” Icarus 128(1), 49–74 (1997).ICRSA50019-1035 http://dx.doi.org/10.1006/icar.1997.5730 Google Scholar

26. B. Hartman and D. Domingue, “Scattering of light by individual particles and the implications for models of planetary surfaces,” Icarus 131(2), 421–448 (1998).ICRSA50019-1035 http://dx.doi.org/10.1006/icar.1997.5861 Google Scholar

27. C. M. Bachmann et al., “A flexible field goniometer system: the goniometer for outdoor portable hyperspectral earth reflectance (GOPHER),” J. Appl. Remote Sens. 10(3), 036012 (2016). http://dx.doi.org/10.1117/1.JRS.10.036012 Google Scholar

28. G. S. Okin and T. H. Painter, “Effect of grain size on remotely sensed spectral reflectance of sandy desert surfaces,” Remote Sens. Environ. 89(3), 272–280 (2004). http://dx.doi.org/10.1016/j.rse.2003.10.008 Google Scholar

29. M. Steffen, “A simple method for monotonic interpolation in on dimension,” Astron. Astrophys. 239, 443–450 (1990).AAEJAF0004-6361 http://dx.doi.org/10.1017/CBO9781107415324.004 Google Scholar

30. A. Savitzky and M. J. E. Golay, “Smoothing and differentiation of data by simplified least squares procedures,” Anal. Chem. 36(8), 1627–1639 (1964).ANCHAM0003-2700 http://dx.doi.org/10.1021/ac60214a047 Google Scholar

31. J. A. Nelder and R. Mead, “A simplex method for function minimization,” Comput. J. 7(4), 308–313 (1965). http://dx.doi.org/10.1093/comjnl/7.4.308 Google Scholar

32. K. Levenberg, “A method for the solution of certain non-linear problems in least,” Q. Appl. Math. 2(278), 164–168 (1944).QAMAAY0033-569X http://dx.doi.org/10.1090/qam/10666 Google Scholar

33. C. M. Bachmann et al., “Phase angle dependence of sand density observable in hyperspectral reflectance,” Remote Sens. Environ. 150, 53–65 (2014). http://dx.doi.org/10.1016/j.rse.2014.03.024 Google Scholar

34. B. Hapke, “Bidirectional reflectance spectroscopy. 3. Correction for macroscopic roughness,” Icarus 59(1), 41–59 (1984).ICRSA50019-1035 http://dx.doi.org/10.1016/0019-1035(84)90054-X Google Scholar


Charles M. Bachmann received his AB from Princeton University (1984) and his ScM (1986) and PhD (1990) from Brown University. He was a research physicist (1990–2013) at the Naval Research Laboratory, serving as head of the coastal science and interpretation section in the remote sensing division (1993–2013). He has been Frederick and Anna B. Wiedman Chair at the RIT Chester F. Carlson Center for Imaging Science since 2013. His research emphasizes coastal hyperspectral remote sensing.

Rehman Eon received his BSc degree in mathematical physics and chemistry from Viterbo University in 2015. He is currently pursuing a PhD in imaging science at the Rochester Institute of Technology. His research interests include the use of optical remote sensing for the assessment of earth sediments and vegetation, radiative transfer modeling, and sensor calibration of Earth-observing systems.

Biographies for the other authors are not available.

© The Authors. Published by SPIE under a Creative Commons Attribution 3.0 Unported License. Distribution or reproduction of this work in whole or in part requires full attribution of the original publication, including its DOI.
Charles M. Bachmann, Rehman S. Eon, Brittany Ambeau, Justin Harms, Gregory Badura, Carrie Griffo, "Modeling and intercomparison of field and laboratory hyperspectral goniometer measurements with G-LiHT imagery of the Algodones Dunes," Journal of Applied Remote Sensing 12(1), 012005 (27 September 2017). http://dx.doi.org/10.1117/1.JRS.12.012005 Submission: Received 4 June 2017; Accepted 25 August 2017
Submission: Received 4 June 2017; Accepted 25 August 2017

Bidirectional reflectance transmission function




Radiative transfer

Multiple scattering

Short wave infrared radiation

Back to Top