Retrieval of Sediment Fill Factor by Inversion of a Modified Hapke Model Applied to Sampled HCRF from Airborne and Satellite Imagery

Abstract

The physical properties of a medium such as density, grain size and surface roughness all influence the angular dependence of spectral signatures. Radiative transfer models, such as the one developed by Hapke, can relate the angular dependence of the reflectance to these geophysical variables. This paper focuses on extracting geophysical parameters, fill factor (decreasing porosity) and the single scattering albedo (SSA), through the inversion of a modified version of the Hapke model of airborne and space-borne imagery. The inversion methodology was validated through controlled experiments within a laboratory setting, where a good correlation (R${}^2$ = 0.72) between the retrieved fill factor and the measured density was obtained. Using the same approach, we also retrieved the sediment fill factor and SSA from airborne data collected by the NASA G-LiHT system, and space-borne data observed by the NOAA GOES imager. The results from these studies provide a mechanism to understand geophysical characteristics of the terrain and may potentially be used for long-term monitoring of the dynamic dunes system.

1. Introduction

Remote sensing techniques are continuously being developed to extract physical information about the Earth’s surface. Over the years, space-borne and airborne sensors have been used for the characterization of surface sediments [1,2,3]. Spectral observations of sediments can be used to effectively identify the physical characteristics of the surface ranging from its texture, roughness, grain size, to its density [1,3,4,5,6,7].

Researchers have recognized soil compaction as an environmental degradation process which hinders the maintenance of soil quality [8,9,10]. The mapping of soil compaction plays a significant role in agriculture, specifically in improving crop production [8,10] and in other environmental and geophysical applications [9,11,12,13,14]. Soil compaction can be quantified by measuring the bulk density or porosity [15]. Ascertaining bulk density at large scales is a difficult task due to variability of the sediment surface [8]. Remote sensing techniques can provide the necessary means of characterizing sediment surfaces over regional and global scales [8].

Physical properties of a sediment surface such as density, grain size, surface roughness, and moisture content influence angular dependence of spectral signatures, specifically the Bidirectional Reflectance Distribution Function (BRDF) [3,4,5,7,16,17,18,19,20,21,22,23,24,25,26,27]. Models based on radiative transfer equations can relate the angular dependence of the reflectance to these geophysical variables [3,5,7,17,18,19,20,28,29]. Extraction of these parameters can provide a better understanding of the Earth’s surface, and plays a vital role in various environmental modeling processes [30,31]. Among the many geophysical parameters describing Earth sediment, the porosity (bulk density) has a significant influence on its spectral properties [3,4,5,6,28,32,33]. Radiative transfer models, such as the one developed by Hapke, typically predicts a correlation between reflectance and density of the material [5,32]. This phenomena has been demonstrated repeatedly in various laboratory studies [3,4,8,34,35,36,37].

The work detailed in Bachmann et al. [3] focused on the retrieval of the fill factor from the inversion of the Hapke model and variants (increasing fill factor corresponds to decreasing porosity). A strong correlation was reported between the retrieved fill factor and the measured sediment density [3]. Modifications were made to the Hapke isotropic multiple scattering approximation (IMSA) model to account for directional dependence in the multiple scattering term [3,38]. In this paper, we take a practical approach to the retrieval of geophysical properties of Earth sediments. We apply the inversion methodology for the modified Hapke model detailed in Bachmann et al. [3] to imagery data collected from airborne and space-borne platforms.

Located in Southern California, our area of study is the Algodones Dunes, a potentially desirable site for the vicarious calibration of space-borne imaging sensors [39]. A major field campaign was conducted in March 2015 to characterize the Dunes. The NASA Goddard LiDAR, hyperspectral, and thermal (G-LiHT) sensor suite [39,40] collected airborne imagery during the field campaign, and the field team collected sediment samples from various sites across the dune system, capturing the variation in its geophysical properties [3,39]. During the campaign, Hyperspectral hemispherical conical reflectance factor (HCRF) measurements of sediments were collected on site using the Goniometer of the Rochester Institute of Technology (GRIT), and later from samples returned to RIT in a laboratory setting using GRIT and a second generation instrument, the Goniometer of the Rochester Institute of Technology Two (GRIT-T) [3,39,41].

In an earlier study, we retrieved fill factor from GRIT-T BRDF measurements of Algodones sediment samples in the controlled conditions of our laboratory [3]. In this study, we demonstrate that fill factor also can be retrieved from airborne imagery from the NASA G-LiHT system collected during the 2015 campaign and from imagery collected by the the Advanced Baseline Imager (ABI) on the Geostationary Operational Environmental Satellite (GOES) series [42]. Specifically, in this study, we retrieved two geophysical parameters of Earth sediments, fill factor and single scattering albedo, through inversion of a modified version [3] of the radiative transfer model developed by Hapke [5]. The specific objectives of this study were to: (1) re-validate the inversion methodology by performing controlled experiments within a laboratory setting using a field sediment sample drawn from a different part of the Algodones Dunes system than what was used to develop the original laboratory-based demonstration of the retrieval in [3], (2) apply the inversion methodology to airborne data from NASA G-LiHT, and (3) to space-borne data collected by the NOAA GOES system.

2. Study Area

The study area, Algodones Dunes, located in southeastern California (${32}^{°}$52′58.44″N, ${115}^{°}$1′8.4″W), is one of the hottest and driest regions in the United States [43]. The Dunes, which are 64 km in length and 6–10 km wide, have formed and continue to evolve through aeolian deposition and can reach heights up-to 100 m (Figure 1) [43,44]. Prevailing winds and the presence of high mountains located in the west and northwest have led to the formation of larger dunes with coarser sands in the west, and drier dunes with finer sediments in the east [44]. The sediments are typically composed of 70–80% quartz, 10–15% feldspar, 5–15% rock fragments, and an assortment of heavy minerals make up the rest [45].

The Algodones Dunes System has recently been identified as a potential site for the inter-calibration of space-borne sensors. The site has shown promise due to its similarity to the pseudoinvariant calibration site (PICS) Libya-4 [39,46]. The 2015 field campaign provided the information required to perform absolute radiometric correction using the Algodones site as a PICS [39]. Portable spectrometers and hyperspectral goniometers obtained ground truth data which characterized both spatial and temporal variability’s within the dunes [3,39,47]. The imaging platform G-LiHT collected complementary hyperspectral and LiDAR data for the characterization of the terrain [39,40].

3. Theoretical Development

3.1. BRDF: Theoretical Background

The bi-directional reflectance distribution function (BRDF) characterizes the angular distribution of reflectance but is a theoretical definition. In practice, due to the finite nature of the sensor and the source, the terms “conical” or “hemispherical” typically describe reflectance distribution functions [33,48,49,50]. Hemispherical conical reflectance factor (HCRF) measurements refer to directional reflectance measurements conducted in outdoor settings. The “conical” term refers to the finite nature of the sensor aperture, while the “hemispherical” term reflects the fact that illumination conditions stem from both direct and diffuse sources [33,48,49]. In laboratory settings, the term biconical reflectance factor (BCRF) describes directional reflectance factor measurements accounting for the finite nature of both the sensor and directional illumination source [33,48,49,50]. The retrieval method described in this research can be applied to both BCRF and HCRF data. The approach is based on inversion of a modified form [3] of the radiative transfer model developed by Hapke [5].

3.2. Hapke IMSA Model

The radiative transfer equations developed by Hapke have been used widely in planetary astronomy and Earth observation research [5,6,18,21,29]. Hapke’s isotropic multiple scattering approximation (IMSA) model is based on the method of invariance and considers five orders of scattering [5,51]. The IMSA model includes terms for single scattering, multiple scattering, the shadow hiding opposition effect (SHOE), and coherent backscatter opposition effect (CBOE) to describe how light scatters from granular materials [5,51,52,53,54]. The solution takes the following form:

$$\begin{array}{cc}\hfill r({\theta }_i,{\theta }_e,g,\lambda )& =K\frac{w\left(\lambda \right)}{4}\frac{1}{{\mu }_i+{\mu }_e}\left(p\left(g\right)[1+B_{s0}{B}_s(g,K,\lambda )]+\left[H\left(\frac{{\mu }_i}{K}\right)H\left(\frac{{\mu }_e}{K}\right)-1\right]\right)\hfill \\ & \times \left[1+B_{c0}{B}_c(g,K,\lambda )\right],\hfill \end{array}$$

(1)

where ${\mu }_i$ and ${\mu }_e$ are the cosines of the incident (${\theta }_i$) and scattered zenith angles (${\theta }_e$), g is the phase angle (the angle defined by the input direction of the light source to a scattering point on the surface, and the direction of observation), $p\left(g\right)$ is the single particle scattering phase function, $B_s(g,K,\lambda )$ models the SHOE with $B_{S0}$ being an accompanying scaling constant, while $B_c(g,K,\lambda )$ represents the angular dependence of the CBOE with $B_{C0}$ a scaling constant, $H\left(\frac{{\mu }_i}{K}\right)$ and $H\left(\frac{{\mu }_e}{K}\right)$ are the Chandrasekhar-Ambartsumian H-functions describing the multiple scattering, $w\left(\lambda \right)$ is the single scattering albedo as a function of wavelength, and K is the porosity coefficient. The porosity factor (K) is a nonlinear function of the fill factor ($\varphi$) [5,32]:

$$K\approx -\frac{ln(1-1.209{\varphi }^{\frac{2}{3}})}{1.209{\varphi }^{\frac{2}{3}}}.$$

(2)

The CBOE dominates at very small phase angles (<2°), where scattering within and between particles produces coherent amplification of observed reflectance [53,55,56]. The models of our experiments described in this paper ignore the CBOE since our BCRF and HCRF measurements do not include sufficiently small phase angles. GRIT-T can not measure at phase angles this small due to the finite extent of the sensor chassis; GRIT-T self-shading occurs for phase angles ≤5° [41]. The IMSA model in Equation (1) also excludes factors of $1/\pi$ and ${\mu }_i$ as the measurements described are typically BCRF and HCRF.

3.3. Shadow-Hiding Opposition Effect

The surge in brightness observed at small phase angles is represented by the SHOE in Hapke’s radiative transfer equation [5,37,54]. The SHOE is explicitly dependent on the fill factor, but implicitly dependent on grain size distribution [5,54]. Hapke’s model of the SHOE takes the following form:

$$B_{SH}\left(g\right)=\frac{B_{S0}}{1+\frac{1}{h_s}tan\left(\frac{g}{2}\right)}.$$

(3)

$h_s$ is the angular width parameter, which depends on both the grain-size distribution and the porosity of the medium. Hapke evaluates the parameter for various different particle size distributions [5]. In a unimodal case, the distribution takes the form:

$$h_s={\left(\frac{3}{8}\right)}^{\frac{3}{2}}K\left(\varphi \right)\varphi .$$

(4)

However, sediment grain size distributions can be quite complicated, with the distribution often consisting of more than one mode as illustrated in the sample that we describe later in Section 4.2. In earlier work, we modified the width parameter to account for multimodal particle distributions [3]:

$$h_s=ϵk\left(\varphi \right)\varphi ,$$

(5)

where $ϵ$ is a scaling constant, with $0\le ϵ\le 1$. We chose this form to retain the functional dependence on $\varphi$ found in the simple distributions modeled by Hapke [5], which differed only by an overall scaling factor. In [3], we optimized $ϵ$ during the inversion process.

3.4. Modified Hapke IMSA Model

The Hapke IMSA model calculates single-scattering contributions exactly; however, the multiple scattering term uses an isotropic scattering approximation in which the single scattering phase function is uniform in all directions: $p\left(g\right)=1$. In recent work, we modified the Hapke IMSA model to include anisotropic scattering of particles by introducing directional dependence into the multiple scattering term [3]. The modification includes an extra factor which depends on the phase function and a scaling parameter $\eta$:

$$\begin{array}{cc}\hfill r({\theta }_i,{\theta }_e,g,\lambda )& =K\frac{w\left(\lambda \right)}{4}\frac{1}{{\mu }_i+{\mu }_e}\left(p\left(g\right)[1+B_{s0}{B}_s(g,K,\lambda )]+\eta p\left(g\right)\left[H\left(\frac{{\mu }_i}{K}\right)H\left(\frac{{\mu }_e}{K}\right)-1\right]\right)\hfill \\ & \times \left[1+B_{c0}{B}_c(g,K,\lambda )\right].\hfill \end{array}$$

(6)

3.5. Inversion Methodology

Our inversion of the Hapke model for the retrieval of geophysical properties relies on the key observation that the single scattering phase function is invariant to the illumination geometry [3], with a geometric dependence only on the phase angle g. Hence, the Hapke model in Equation (6) can be reorganized in terms of a scaled version of the observed reflectance and the multiple scattering term:

$$p\left(g\right)\left[1+B_{s0}{B}_s(g,K,\lambda )\right]=\frac{4}{Kw\left(\lambda \right)}({\mu }_i+{\mu }_e)r({\theta }_i,{\theta }_e,g,\lambda )-\eta p\left(g\right)\left[H\left(\frac{{\mu }_i}{K}\right)H\left(\frac{{\mu }_e}{K}\right)-1\right].$$

(7)

The approach requires the acquisition of BRDF data at different illumination geometries, and minimizes the residual between the single scattering phase function for the two sets of data:

$$\underset{\varphi ,w\left(\lambda \right)}{min}\left[{\left\{p_1\left(g\right)\left[1+B_{s0}{B}_{s,1}(g,K,\lambda )\right]-p_2\left(g\right)\left[1+B_{s0}{B}_{s,2}(g,K,\lambda )\right]\right\}}^2\right].$$

(8)

The two-parameter search over the single scattering albedo, $w\left(\lambda \right)$, and the fill factor, $\varphi$, uses the gradient-based Nelder–Mead simplex method [57]. Bachmann et al. [3] previously described the work-flow using this inversion method.

The inversion of the Hapke parametric model is a daunting problem due to the numerous parameters involved. Earlier studies had only limited success in relating Hapke model parameters to real physical properties of soil sample or regoliths. These challenges are detailed extensively in [6,37], which incorporate both gradient descent and grid search methodologies. Our inversion method, described in [3], takes a different approach in the optimization procedure to improve the retrieval of fill factor and single scattering albedo of sediments from the inversion of a modified Hapke model. In addition to using a modified form of the Hapke model, which incorporates an isotropic multiple scattering term, our inversion relies on the minimization of residuals between the single scattering phase function from at least two different illumination conditions. This choice leads to an inversion depending only on two parameters, the porosity and single scattering albedo, compared to the original ten free parameters of the Hapke model. Other differences between past approaches and our optimization scheme include the fact that past approaches used least-squares gradient descent and grid search optimization techniques [6,37], while our approach [3] uses a simplex-based search method, the Nelder–Mead Simplex Method [57]. Additional advantages of using the Nelder–Mead optimization include a decreased processing time and more repeatable retrieval process. We also found that the best optimization results relied on using regularization to achieve numerically stable solutions [3]. Our modified Hapke model, which no longer assumes isotropic multiple scattering, also provided consistently better and numerically stable results compared to the original model when using our inversion methodology. This was quite evident when we changed the starting values and step size in our optimization using the Nelder–Mead Simplex method. The use of all these techniques within our optimization to invert the Hapke model and our variant of this model has consistently given us a meaningful retrieval of the sediment fill factor and single scattering albedo from both controlled experiments in a laboratory setting, described both in [3] and in this paper, as well as from air-/space-borne data also described in this paper.

Studies have shown the inherent dependence of reflectance on the fill factor (decreasing porosity) of the medium, when the particles are significantly larger than the wavelength [32]. Hapke’s model predicts that the reflectance of the medium typically increases with increasing fill factor. As the fill factor increases (decreasing porosity) and the medium becomes more opaque to the light passing in between the particles [32], eventually a threshold is reached where the particles are sufficiently closely packed that the scattering is more like that observed from larger single particles, with light reflecting off the surface less efficiently [32]. Hence, at high fill factors, a noticeable decrease in reflectance can be observed [5].

The single scattering albedo (SSA), $w\left(\lambda \right)$, at a given wavelength is the ratio of the total amount of light scattered to the total amount of light scattered and absorbed by the particle [5]:

$$w\left(\lambda \right)=\frac{S}{S+A},$$

(9)

where S and A are the fraction of light scattered and absorbed by the particle respectively. The SSA depends on the optical properties of the medium, determined by factors such as composition, grain-size distribution, shape and structure [37,51,52]. Non-absorbing media normally have high SSA, while media with large particles and high index of refraction generally have smaller SSA [37]. For materials such as sediment, where particles are closely packed and large, SSA increases with decrease in particle size [37].

4. Laboratory Studies

4.1. GRIT-T: Design and Instrumentation

The goniometer of the Rochester Institute of Technology-Two (GRIT-T) is a second-generation system designed to obtain BRDF measurements in both laboratory and field settings [3,41,58]. Figure 2 illustrates the different components of the GRIT-T system, deployed in various settings. The instrument consists of two Analytical Spectral Devices (ASD) FR4 spectro-radiometers to simultaneously provide directional radiance measurements from the target plane while recording the directional downwelling radiance from the sky. The ASD spectrometers provide spectral measurements from the visible through the short-wave infrared (350–2500 nm). The design of GRIT-T ensures positioning of the fore-optic over an azimuthal range of 0° to 180°, and −70° to +70° in zenith. This allows sampling of the target surface over the complete hemisphere. The target-plane tracking of the instrument minimizes any variations in the field of view (FOV). As a result, it has an angular position accuracy of $\pm 0.2$° in both azimuthal and zenith sensor positions. The small-form factor sensor head assembly minimizes self-shading. GRIT-T produces digital elevation models using two different methods: from an onboard laser measurement unit [41] and at higher resolution from an on-board field-of-view camera, which provides stereo views of the measurement location as input to a structure-from-motion algorithm [58].

4.2. Laboratory Measurements

Our laboratory studies in this work included BCRF measurements of sediment samples acquired from the Algodones Dunes during the 2015 field campaign. These studies involved manipulating the density of the samples, and examining how this geophysical variable affected the observed BCRF. These data are the basis of the fill factor and single scattering albedo joint retrievals using our modified Hapke optimization model in the laboratory portion of the present study. The Algodones sediment sample 1306-M-03 acquired from the northern end of the dunes system (${32}^{°}$58′45.00″N, ${115}^{°}$7′47.00″W) was used for the laboratory analysis. The grain size distribution along with a microscopic images of the sediment sample is shown in Figure 3. The smaller peaks of the secondary modes are due to the presence of finer sand and silts within the mixture.

After mechanical sieving, we found that the largest grain size fraction was for the 450 $\mathsf{\mu }$m sieve; particles of this size will result in scattering described by the geometric optics regime. The microscopic image in Figure 3 demonstrates the complexity of the material observed in the grain-size distribution and shape as well as its variety in mineral composition. The complexity of the mixture is also reflected in the multi-modal distribution of the particle size. In [3], we introduced a modification to the width parameter of the SHOE, which depends on the grain size distribution, to allow greater flexibility in representing the possibility of multi-modal distributions [3].

The sample was manipulated via pluviation to achieve a series of different densities [59]. The drop height of the sediment in the apparatus is correlated with the resulting relative density of the sediment, and we use this approach because past analysis indicates that results are more repeatable than an existing American Society for Testing and Materials (ASTM) standard [59]. Other approaches to sample preparation, such as the use of a wind tunnel, have also been studied [26,27]. Such approaches arguably can produce distributions, and in particular gradients, across a larger area that are important in analyzing the details of Aeolian processes; however, our goals here were much simpler, namely the production of a uniform relative density of the sediment that was highly repeatable, both of which result from the use of a pluviation approach [59]. In our laboratory, we used pluviation to achieve ten different densities for the sediment sample 1306-M-03 acquired from the northern end of the Algodones Dunes. The GRIT-T system then measured the BCRF for each sample preparation. The objective of the laboratory measurements was to validate the radiative transfer model and the retrieval process detailed in Section 3 and in [3].

The samples were measured at two different illumination geometries (20° and 60°) for each sediment density preparation. The minimization process described in Equations (7) and (8) uses two data sets acquired at the different illuminations as input. The lab set-up for performing the series of BCRF measurements is shown in Figure 4. Figure 4a shows the set-up with the illumination source at a zenith of 60°. The source is a 300 Watt Fresnel lamp, which provided collimated light onto the target plane. The radiance measurements from the ASD spectrometer were referenced to a Spectralon^® panel [60], which approximates a Lambertian surface, shown in Figure 4b.

4.3. Spectral Analysis and Fill Factor Retrieval

The reflectance spectrum for the sample 1306-M-03 at a density of 1.5252 g/cm${}^3$ is shown in Figure 5. The figure displays the reflectance spectra for sensor positions over the complete hemisphere, with the illumination source at a zenith angle of 60°. The colorbar in the figure corresponds to the sensor azimuth. All of the lab measurements in this study followed a similar trend. Several different factors influence the reflectance spectrum of the sediment: the presence of organic matter, mineral composition, roughness, particle size distribution, and density [61]. The reflectance increases with wavelength in the visible part of the spectrum (400–700 nm). There are weak absorptions bands at wavelengths less than 1000 nm due to the presence of iron oxides within the sample [62]. The strong absorption peaks observed at wavelengths 1450 nm and 1950 nm are due to water and hydroxyl bands [61]. The presence of clay materials influences the absorption peaks observed near 2200 nm [62].

Figure 6 shows the BCRF plots observed at the two different illumination angles (20° and 60°) for five different wavelengths from the visible and near infra-red (VNIR) to the short wave infra-red (SWIR). As the figure shows, the geometry of the source strongly influences the structure of the BCRF. There is a greater amount of diffuse scattering at illumination angles closer to the nadir (20°). When the illumination angle is closer to nadir, there is more multiple scattering within the medium, resulting in more diffuse reflectance. At large zenith angles (60°), there are more prominent backward and forward scattering lobes. A “bowl-shape” at off-nadir positions appears due to volumetric scattering is also evident. The influence of wavelength on the BCRF plots is also quite apparent, especially for BCRF measurements at a zenith angle of 60°.

The inversion of the radiative transfer model detailed in Section 3 used the two sets of measurements conducted at different incident illumination geometries. As described in Section 3, the modified Hapke model was inverted to retrieve the fill factor and the SSA. Figure 7 shows the retrieved fill factor and the SSA for ten different densities ranging from 1.4542–1.6904 g/cm${}^3$. The left part of Figure 7 displays the fill factor (red dots), which for each sample density prepared, was an average over the wavelength range of our instrument (350–2500 nm). The inversion jointly retrieved the fill factor and the SSA, which likewise appears in the right part of Figure 7 with the mean and standard deviation of the SSA being over the optimization runs, one for each density. The SSA provides insight into the optical properties of the of the individual particles in the medium. For granular materials, the SSA typically increases with a decrease in the particle size [37]. The final fill factor was an average of the fill factors obtained at each wavelength by the inversion process, since density is invariant to wavelength. In our inversion methodology, we introduced an additional constraint to ensure the retrieved fill factor remains consistent across all wavelengths [3]. However, the results obtained from the inversion show the fill factor values still vary slightly between the visible near infrared and the short wave infrared. A previous study [3] also showed a similar result. The difference in the fill factor value between the VNIR and SWIR falls within the standard deviation observed within the range of the fill factor across all wavelengths. Over the set of densities explored in this study, we found a moderately good correlation between the retrieved fill factor and the measured density, with a R${}^2$ value of 0.72, as shown in Figure 7. The retrieved fill factor from the original IMSA model (not reported in this paper) varied significantly with wavelength, and also produced less numerically stable results when compared to the modified model. We note that the retrieved geophysical properties, which we have obtained for a sample derived from the northern end of the Algodones Dunes system, achieves a level of accuracy similar to this same approach applied to a sample derived from the western side of the central dune system in our earlier study [3], confirming the validity of the approach.

Having validated the inversion methodology in a set of controlled laboratory experiments, we now turn our attention to the retrieval of the fill factor from data collected by airborne (Section 5) and satellite (Section 6) systems.

5. NASA Goddard’s LiDAR, Hyperspectral, and Thermal (G-LiHT) Airborne Imager

5.1. G-LiHT: Design and Instrumentation

An integral part of the 2015 field campaign was obtaining hyperspectral and LiDAR airborne data using NASA Goddard’s G-LiHT system [39,40]. The hyperspectral imager, designed by Headwall Photonics Inc. (Bolton, MA, USA), is a pushbroom system providing spectral measurements from 400 to 1000 nm at a spectral resolution of 1.5 nm with an FOV of 50° [40]. The LiDAR is a VQ-480 airborne laser scanner (Riegl USA, Orlando, FL, USA). The LiDAR [40] provided a detailed digital elevation model of the terrain [39]. The flight line and the spatial coverage of the G-LiHT system over the Algodones Dunes system is shown in Figure 8. The G-LiHT system provided hyperspectral imagery, a LiDAR-derived digital elevation model, as well as thermal imagery at a spatial resolution of 1 m.

Airborne data was collected using the G-LiHT system from Monday 9 March 2015 to Friday 13 March 2015. There were several flight lines (A-,T-, and E-lines) at various different orientations flown over the course of the field campaign. For sets of flightlines flown on the same azimuthal heading, successive lines in the set have significant overlap with neighboring flight-lines. The combination of all these different flights provided hyperspectral imagery obtained from multiple imaging and illumination geometries over the same region on the ground. For one of the typical areas of maximum overlap used in our analysis, G-LiHT obtained hyperspectral imagery from 16 different view geometries over the same 102 m by 139 m region on the ground. Figure 9 shows a mosaic of the 16 different flight lines along with the region of interest (ROI) for our study. The multiple view-geometries provided the necessary range of phase angles to perform the inversion of our modified Hapke model using the acquired hyperspectral G-LiHT imagery.

5.2. Spectral Analysis and Fill Factor Retrieval from G-LiHT Imagery

Figure 10 illustrates the spectral library and polar plots for a pixel located at ${32}^{°}$54′50.7553″N, ${115}^{°}$6′53.8587″W. The spectral library in Figure 10a plots the reflectance for all 16 view-geometries in the visible and near-infrared (400–1000 nm). The spectral characteristics are consistent with the observations in the laboratory studies, mentioned earlier in Section 4. The spectral library illustrates the angular dependence of the reflectance of the target from airborne platforms, with the colorbar representing the phase angle. It should be noted that the solar/sensor angular calculations from G-LiHT take into account the slope aspect of the dunes. We used the accompanying LiDAR data from the flight to correct for the slope of the dunes in our angle final calculations. The overall intensity of the reflectance is distinctly dependent on both illumination and view geometry, as shown in the plot. The scattering of the target pixel is strongest in the backscatter direction, with the phase angle varying from typically 20° to 80°. This range of phase angles was evident for all the pixels within the scene.

Figure 10b,c shows the polar plots at wavelengths 551 nm and 807 nm, respectively. The plots shown are a composite of the 16 different looks from the G-LiHT overpasses for a single pixel within the ROI. The sensor zenith angles are less than 30° for the different looks, which is typical of measurements taken from airborne and space-borne platforms. The view geometries, and the consequent HCRF, for each pixel within the ROI are fairly different from each other, as illustrated in Figure 11. The figure shows the polar plots at a wavelength of 714 nm for four pixels across the ROI. The reflectance measurements are significantly different for each of the pixels, influenced by the geometries as well as the geophysical properties of the regolith.

In contrast with our lab studies using the goniometer in Section 4, the image-derived HCRF from G-LiHT does not provide us with sensor positions over the complete hemisphere. However, the imagery set does provide a comparable range of phase angles due to both varying heading of the aircraft overflight as well as the change in solar position throughout the day. Thus, we could perform inversion for a similar range of phase angles to that used in our laboratory study in inverting our model of the G-LiHT imagery data, and we were able to get similar values for the fill factor and SSA. As a result, we believe the range of the phase angle for each pixel in the ROI provided satisfactory information about the scene to perform inversion of the Hapke model to retrieve the fill factor. The range of view-geometries for each pixel in the ROI provided satisfactory information about the scene to perform inversion of the Hapke model to retrieve the fill factor. In the inversion procedure, we assumed that neighboring/adjacent pixels had the same geophysical properties. Since the G-LiHT system provided imagery at high spatial resolutions (1 m), it was a reasonable assumption for the inversion process. As a result, the spectra from adjacent pixels served as the two input data sets for inversion using Equation (8). Figure 12 shows the retrieved fill factor for the 102 m by 139 m region. Since the fill factor is proportional to the bulk density of the sample [63,64], we can see, therefore, that the sediment density varies across the terrain. The primary goal of 2015 campaign was focused on acquiring calibration data in support of inter-satellite calibration across the extensive dune system [39]. Although flights occurred over sites where geotechnical data, such as bulk density, was collected, the flight plans were not specifically designed to maximize the number of views of each calibration site, but rather to ensure that there was G-LiHT data collected across the entire dune system, and at least some coverage of the calibration sites. For the geotechnical calibration sites, there were typically only a few overlapping flight lines, and not the higher degree of overlap that occurred circumstantially near but not at the precise calibration locations. The particular spot chosen for our imagery ROI was located in the zone of highest overlap where inversion of the radiative transfer could be most successfully applied. Thus, we do not have any ground bulk density data from within the 109 m by 139 m ROI. However, we did take bulk density measurements from locations close to our ROI (1 km away), and we used the regression of the fill factor vs. density from lab studies in Section 4 (shown in Figure 7) to convert the relative fill factor to bulk density for the ROI. We observed that the range of density (based on the regression) varies from 1.1 to 2.15 g/cm${}^3$, which is comparable to bulk density measurements obtained during the field campaign (

Table 1. Bulk density measurements conducted during the 2015 Algodones field campaign from sites near our ROI.

Sample Name	GPS Location	Bulk Density (g/cm${}^3$)
BC01-01	${32}^{°}$55′8.25″N, ${115}^{°}$7′2.42″W	1.96
BC01-02	${32}^{°}$55′9.47″N, ${115}^{°}$7′0.06″W	2.25
BC02-01	${32}^{°}$54′54.77″N, ${115}^{°}$6′33.03″W	2.18
BC02-02	${32}^{°}$54′55.35″N, ${115}^{°}$6′33.91″W	2.08
BC03-01	${32}^{°}$54′55.07″N, ${115}^{°}$6′34.76″W	2.02
BC03-02	${32}^{°}$54′53.98″N, ${115}^{°}$6′34.26″W	2.19

).

The retrieved single scattering albedo is plotted in Figure 13. The figure shows the average retrieved SSA from all pixels over the ROI, and the corresponding standard deviation. We also report the retrieved fill factor for four locations within our ROI, while the HCRF spectra for these four individual pixels are shown in Figure 11.

6. The Geostationary Operational Environmental Satellite (GOES)

6.1. GOES: Design and Instrumentation

The GOES satellite series are a collaborative venture between the National Aeronautics and Space Administration (NASA) and the National Oceanic and Atmospheric Administration (NOAA) [65]. The 16-channel Advanced Baseline Imager (ABI), the primary imaging system on the GOES-R series, is used for a variety of environmental applications relating to the weather, land, ocean, and the atmosphere [42].

The ABI is a multi-spectral imaging system, with 16 different spectral bands from the visible to the longwave infrared (0.47–13.3 $\mathsf{\mu }$m) [42]. The spatial resolution of the system varies from 0.5 km to 2 km depending on the spectral band [42]. The imaging system has two different scanning modes, providing excellent coverage rate. It is capable of collecting a full disk (Western Hemisphere) image every 15 min, images of the Continental United States (CONUS) every 5 min, and also a selectable 1000 km by 1000 km region every 30 s [42]. A pseudo-RGB image of the CONUS on 1 July 2017 is shown in Figure 14. The ABI system also provides on-orbit calibration for all 16 bands, minimizing errors due to degradation over time.

6.2. Spectral Analysis and Fill Factor Retrieval from GOES Imagery

The analysis of the GOES imagery was performed on four bands from the visible to the near infra-red; band 1 (0.47-$\mathsf{\mu }$m), band 2 (0.64-$\mathsf{\mu }$m), band 3 (0.86-$\mathsf{\mu }$m) and band 5 (1.6-$\mathsf{\mu }$m). Band 2 has a spatial resolution of 0.5-km, while the remaining bands used in this study have a resolution of 1-km. The spatial resolution for bands 1, 3, and 5 was interpolated to 0.5-km. Figure 14 shows ABI imagery of the Algodones Dunes, a 23.5 km by 17.5 km region. Although the ABI instrument has the same view-geometry of the surface, as the solar position changes, the excellent temporal coverage of the system provides the necessary range of phase angles to perform the inversion. Unlike the angular calculations for the G-LiHT study in Section 5, we did not account for the slope aspect for the GOES imagery, as we did not have any accompanying LiDAR data and also the spatial resolution of 0.5-km made it difficult to accurately correct for the slope of the dunes. Figure 15 shows the variation of the phase angle from sunrise to sunset on 1 July 2017. Since the phase angles are similar before and after noon, they act as the two sets of data required to perform the inversion in Equation (8).

Figure 16 shows a map of the retrieved fill factor obtained from the GOES imagery through the inversion of the modified Hapke model. There is not a large variance in the retrieved fill factor across the dunes from the GOES imagery, as can be observed from the colorbar in Figure 16. This is possibly due to the fact that GOES-R ABI imagery has a ground sample distance (GSD) of 0.5 km. The low spatial resolution of the system averages a substantial region on the ground, consisting of different materials, which ultimately affects the retrieved fill factor. However, at this spatial resolution, qualitatively the distribution of the fill factor does highlight certain characteristics of the dunes. The fill factor generally increases with height, and we observed this relationship when retrieved fill factor was compared against LiDAR elevation data collected by the G-LiHT sensor during the 2015 field campaign. Figure 17 shows side-by-side images of the fill factor and DEM of the dunes. The spatial resolution of the fill factor image was interpolated from 0.5-km to 1-m to match the GSD captured by the LiDAR system. The sediment fill factor typically increases (decreasing porosity) with the height of the sediment medium [66]. The surface sediments compress the materials underneath, increasing the density with height [66]. Figure 17 shows that high retrieved fill factors correlate well with the peaks of the higher dunes observed in the middle of the desert. The edges of the desert, where there is more of a basin or depression, have lower fill factor values associated with lower density. The northeastern part of the desert in Figure 17, circled in black, consists mostly of rock formations and dense vegetation, an area for which our sediment retrieval model for the fill factor does not apply, and produced a high fill factor. As expected here, our retrieval model has only a weak correlation with the LiDAR data.

7. Monitoring Large-Scale Changes in the Sand Dunes

The GOES imager is an ideal system for monitoring the Algodones Dunes. Its high temporal resolution can observe changes in the fill factor, and thus the bulk density, in response to meteorological phenomena on short time scales.

Wind-blown sands have formed the Algodones Dunes whose topography changes continuously due to the influence of local wind patterns [45]. Dunes are generally steeper in the windward side and shorter on the lee side [67,68]. The common mechanism for the wind induced sand movement is saltation, in which the wind carries aloft individual grains from the surface and eventually deposits them back on the surface [67,69]. The movement of particles within the dunes generally depends on particle characteristics (grain size, shape and density) and the wind speed [69]. Saltation occurs for particles of size approximately 70–500 $\mathsf{\mu }$m, a range encompassed by typical grain size distributions observed in the Algodones Dunes in previous studies [3,68]. This wind induced sand movement takes place almost entirely within the top 0.5 m of the desert surface, with 90% of the movement happening within 2.5 cm [67]. For fine-grained particles, the minimum velocity required for the movement of sand is approximately 5 m/s at 1 m above the surface [67]. A relationship between wind speed and grain size in which saltation occurs is described in further detail by [70].

The high rate of image acquisition by GOES along with the inversion process detailed in this article can be used to monitor changes in the dunes caused by meteorological influences such as wind and precipitation.

Figure 18 illustrates the fill factor retrieved from GOES ABI cloud-free imagery of the Algodones Dunes during the period from 1 July 2017 to 28 July 2017. The figure also displays the standard deviation over the month of July. We obtained meteorological data from the Global Historical Climatology Network (GHCN) [71]. The database provides historical weather data from numerous stations around the globe. The average daily wind speed (AWND), maximum temperature (TMAX), and minimum temperature (TMIN) for each day is also shown in Figure 18. The meteorological data were recorded by the Marine Corps Air Station (MCAS) in Yuma, AZ, USA. Our retrieved fill factor from the GOES imagery of Algodones varies from day-to-day due to the process of saltation induced by the prevailing wind. It should be noted that wind pattern is not the only factor influencing the retrieved fill factor for a particular day. For example, the wind conditions are similar for 2 and 15 July; however, the fill factor is drastically different for these two days due to changing wind patterns leading up to that day as well as other environmental factors. The largest change in the retrieved fill factor occurs in the center of the dunes as illustrated by the standard deviation in Figure 18. This region generally consists of active dunes with medium to fine sands, and the topography is susceptible to the prevailing wind. Winds have less impact on the northern end of the dune system, and, to a lesser extent, also the southern end. The presence of vegetation and rock formations within these regions of the Dunes stabilizes the sediments, minimizing the influence of winds on the topography.

8. Conclusions

The physical properties of sediments such as density, grain size and surface roughness influence the angular dependence of its spectral signature. Models based on radiative transfer equations, such as the one developed by Hapke, can relate the angular dependence of the reflectance to these geophysical variables. This paper focused on extracting geophysical parameters, fill factor (decreasing porosity) and the single scattering albedo, through the inversion of a modified Hapke model of airborne and space-borne imagery. The area of study was the Algodones Dunes located in Southern California. The Algodones Dunes System is considered a potentially desirable site for the vicarious calibration of space-borne imaging sensors, and the study detailed in this paper provides a better understanding of the physical characteristics of the sediment surface.

We validated the inversion methodology by performing controlled laboratory measurements on sediment samples from the Algodones Dunes. We used GRIT-T to measure hyperspectral BCRF of these samples, and we found a good correlation between the retrieved fill factor and the measured density in the present laboratory experiments, corroborating experiments done in previous work [3]. After validating the inversion methodology, we demonstrated the retrieval of the fill factor from angular dependent reflectance data derived from airborne and satellite imagery time series.

We applied the inversion methodology to airborne hyperspectral imagery collected by the NASA G-LiHT system during the 2015 Algodones Dunes field campaign [3,39]. The significant overlap between the G-LiHT flight lines during the campaign provided imagery with multiple geometries over the same region on the ground. In total, there were 16 different view geometries over the same 102 m by 109 m region used in this study. The retrieved fill factor and SSA provided a better understanding of the variability of the terrain.

We also retrieved the fill factor from NOAA GOES ABI imagery. The inversion used four bands from the visible to the near infra-red. The high temporal coverage of the system provided the necessary range of phase angles to perform the inversion. There was not a large variance in the retrieved fill factor due the low spatial resolution (0.5-km) of the GOES system. The low spatial resolution of the system averages a significant portion of the surface, smoothing variations in the retrieved fill factor. However, qualitatively, the retrieved fill factor still highlighted some characteristics of the dune. The fill factor generally increased with height, and this correlation was clearly visible when compared to LiDAR data from the Dunes. We also used the high temporal resolution of the GOES imager to monitor changes in the dunes associated with meteorological phenomena. In this paper, specifically, we observed changes in the fill factor or density as a function of the intensity of the prevailing winds. The density of the dunes seems to vary day-to-day, due to the process of saltation induced by the winds. Prevailing winds had the greatest impact on the center of the Dunes system; the predominant sediment has medium to fine sands there, which are most susceptible to the winds. The presence of vegetation in the northern end of the Dunes, and, to an extent the southern end, mitigates the effects of the winds, and these trends appear in our retrieved fill factor.

The research in this paper detailed the inversion of the modified Hapke model to retrieve the fill factor and the single scattering albedo originally described in [3]. In this paper, we demonstrated that this method, which was originally developed and tested using laboratory data, could be extended successfully to work with inputs from multi-temporal airborne and satellite imagery. Future work may include a more targeted Algodones-2 campaign to further extend the analysis detailed in this paper, i.e., the flight plan for the airborne sensor would focus on a larger overlap region on the surface with significantly more view-geometries, and also a larger number of samples would be taken at various locations within the overlap zone to provide reference data. We have also recently conducted a field campaign in another field setting where both extensive ground reference data and overlapping hyperspectral imagery from different sensing platforms were collected. These data will be used in a rigorous verification of the initial results that we have obtained in the present study.