A Comprehensive BRF Model for SpectralonⓇ and Application to Hyperspectral Field Imagery

We present a comprehensive hyperspectral bi-directional reflectance model for a pristine optical-grade white Spectralon panel with full azimuthal coverage, a wide range of view zenith angles, and incident illumination zenith angles, throughout ultraviolet (UV), visible near-infrared (VNIR), and shortwave infrared (SWIR) wavelengths (350–2500 nm). Measurements were acquired using the Goniometer of the Rochester Institute of Technology-Two (GRIT-T), which incorporates Analytical Spectral Device (ASD) FieldSpec Full-Range 4 (FR4) spectroradiometers. Residual plots show that the empirical model is accurate to within 1%–2% reflectance within most of the observing hemisphere. We demonstrate an application of the Spectralon panel model to hyperspectral imagery (HSI) acquired during a field experiment at the RIT Tait Preserve where a pair of field Spectralon panels was imaged from various view geometries from both drone and mast-mounted hyperspectral imaging systems.

Reflectance data of various Spectralon products presented in the literature are often limited to the principal scattering plane, a limited number of illumination geometries, and a small range of wavelengths [1], [2], [3], [4], [5], [6], [10], [11], [12], [13], [15], [16], [17].Although these analyses provide sufficient characterization for use in laboratory settings, the extensive possible observational and illumination geometries in Earth remote sensing field applications necessitate a more comprehensive reflectance model of Spectralon.Acquiring such measurements in both observer azimuth and zenith angles requires, at a minimum, a two-axis goniometric system for measuring the sample and an additional angular degree of freedom for varying the illumination source geometry.Accurate spectral measurements over a large range of wavelengths also require a broadband source and a calibrated spectroradiometer.
The necessary infrastructure for more comprehensive directional measurements of samples just described, in general, has inhibited the acquisition and dissemination of such reflectance products.To support the Earth remote sensing community for the processing and calibration of accurate reflectance products, we present an easily accessible comprehensive UV-VNIR-SWIR reflectance model for a pristine optical-grade white Spectralon panel.We then also demonstrate the model in realistic field settings for hyperspectral imaging applications.Measurements and the source code for our model can be found on GitHub [19].To achieve our modeling goals, we expand upon previous measurements and modeling, primarily done within the principal plane in the literature, by including full azimuthal coverage.Our modeling also includes sensor view zenith angles up to 70 • and incident illumination zenith angles between 10 • and 70 • .For angles beyond these ranges, the model can be extrapolated.

A. Geometric Considerations for Reflectance
The bi-directional reflectance factor (BRF) quantifies light scattered from a surface or object of interest at a specific incident illumination geometry to another specific observation direction, referenced to the scattering of light from a perfectly diffuse (Lambertian) surface illuminated and observed in the same geometry [20], [21].Specifically, it is defined as the ratio between the radiance reflected by a sample surface and the radiance reflected by a perfectly reflecting Lambertian surface that is irradiated and observed under the same conditions as the sample.The BRF can be written as BRF(θ i , φ i ; θ e , φ e ) = L r (θ i , φ i ; θ e , φ e ) L id (1) for incident azimuth and zenith angles φ i and θ i , view azimuth and zenith angles φ e and θ e , reflected sample radiance L r (θ i , φ i ; θ e , φ e ), and the ideal Lambertian surface reflected radiance L id .For our analysis, we reduce the azimuthal variables to a relative azimuth such that φ = φ i − φ e .Due to the uniform composition of the Spectralon panel, we assume that the BRF is independent of the panel's azimuthal orientation under any lighting geometry; that is, the azimuthal dependence to characterize is the view azimuth relative to the incident azimuth.Our model also assumes symmetry in reflectance across the principal scattering plane ( φ : 0 • → 180 • ).This simplifies the panel model, making it easier to apply in both field and laboratory settings.We note that the terminology "BRF" is a theoretical, idealized quantity.Our measurements, in general, are more precisely described by the terminology bi-conical reflectance factor (BCRF), which explicitly represents the finite nature of the source and the sensor aperture in its definition [20], [21].However, in our discussion, we will use the generic terminology BRF to describe our practical measurements and modeling.

II. GONIOMETER OF THE ROCHESTER INSTITUTE OF TECHNOLOGY-TWO (GRIT-T)
GRIT-T [22], [23] is a second-generation field and laboratory goniometric system designed and developed at RIT. GRIT-T incorporates two Analytical Spectral Device (ASD) FieldSpec Full-Range 4 (FR4) spectroradiometers and has been used in a wide variety of remote sensing applications [23], [24], [25], [26], [27], [28], [29], [30], [31], [32], [33], [34], [35], all of which also incorporate Spectralon measurements for converting radiance data into reflectance products.In field settings, the two spectrometer fields-of-view (FOV) are co-aligned along a common axis.One spectrometer records downwelling radiance, and the other records scattered radiance from the surface or sample of interest.In laboratory settings, normally just the downward-looking spectrometer is used since the light source is directional and there is no diffuse light present.Following a user-defined scan pattern in its control software, GRIT-T's motorized carriage, arm, and head orient the FOV of the fiber optics of the ASD FR4 spectrometers.
For the downward-looking spectrometer, the system tracks a common measurement point on the surface using a rotating head on the end of the rotating arm to ensure a common measurement area regardless of surface configuration.This is enabled by a laser distance measurement from the nadir direction at the start of the measurement cycle, which allows the system to minimize parallax errors as the rotating arm, combined with the head rotation, determines the zenith angle of the fore-optic.The carriage translates the arm and head assembly in azimuth to allow full coverage of the observation hemisphere.GRIT-T appears in our laboratory configuration in Fig. 1 with the pristine Spectralon panel (model: SRT-99-180) used in our study.A 3 • fore-optic with an optical scrambler was used for all Spectralon panel measurements.The plot in Fig. 2 shows the 1σ uncertainty in the BRF measurements as percentages.The uncertainty was calculated by computing the standard deviations of the repeated nadir measurements of the Spectralon at each incidence angle.Although there is increased uncertainty near 1000 nm and for longer wavelengths in the SWIR, the uncertainty is <1% across all wavelengths.For oblique viewing geometries, the measurement circle elongates into an ellipse.At the largest view and incident zenith angles, both at 70 • , the difference in the irradiance onto the panel is the largest between the front and back ends of the measurement spot, which are −2.3 and +2.3 cm from the measurement origin, respectively.The relative difference in nadir-viewing radiance measurements between the front and back of the elliptical measurement spot with respect to the origin is −4% and +4%, respectively.Within the ellipse, we measured a linear relationship between the radiance as a function of measurement position on the panel; therefore, the effects of differences in irradiance due to oblique illumination and viewing geometries are largely averaged out.As described previously, the second ASD FR4 spectrometer can be used for sky radiance measurements in the field but was not used in the laboratory setting.The spectral range of the ASDs is 350-2500 nm in 1 nm steps with 3 nm spectral resolution in the VNIR and 8 nm in the SWIR.
The illumination source is an unpolarized broadband QTH lamp attached to an external motorized arm.The QTH lamp is also connected to an RA Series SCR preregulated linear power supply from Mid-Eastern Industries that provides stable dc power to our lamp.As a standard procedure, we use a 1 h warm-up period for our illumination source and ASD spectrometers onboard GRIT-T to ensure instrument stability.The external arm is affixed to the same optical table as GRIT-T and rotates to provide any desired incident zenith angle along the principal scattering plane.The orientation of the lamp illumination is aligned with the center of the sample area of GRIT-T's FOV.III.BRF DATA OF SPECTRALON Fig. 3 shows an array of polar colormaps of a spectral subset of the Spectralon panel BRF data measured by GRIT-T.The color axis denotes the BRF values of the panel in the observing hemisphere.Over the observing hemisphere, the data scan pattern provided measurements in 10 • steps in azimuth and 5 • steps in view zenith from 0 • to 70 • for a total of 504 spectral radiance measurements.The same scan pattern was used in each of the BRF measurement scans for seven incident zenith angles ranging from 10 • to 70 • in 10 • steps for a total of 3528 spectral radiance measurements, each covering the spectral range from 350-2500 nm in 1 nm increments.A scan pattern of equal degree increments was used instead of points with equal increments in cosine space; by default, the control software designates scan patterns in equal degree increments, which is what was initially used.The repeated nadir measurements along different lines in azimuth are also used to correct for any radiometric drift during each scan.Points in the backscattering region that were self-shaded by the goniometer head (the black circles in Fig. 3) were not used in our analysis.The measured spectra were initially recorded in raw digital counts, which were calibrated to radiance using standard procedures in ViewSpecPro with the most recent ASD system calibration.At each incident zenith angle, all point measurements were divided by the nadir measurement and converted to BRF values using the procedure described in Section IV.

IV. BRF MODEL OF SPECTRALON
We adopted the general form of the empirical BRF model developed by Lévesque and Dissanksa [13], [14], which consists of four physically motivated scattering terms, the Spectralon reflectance spectrum, and a normalization function.We expand on their framework by including the azimuthal dimension, φ, in the model BRF( φ, θ e , θ i , λ) = D( φ, θ e , θ i , λ) + F( φ, θ e , θ i , λ) where  In the original Lévesque-Dissanska model, the diffuse component depends on low order powers of the view and incident zenith angles θ e and θ i (Table I).During the initial stages of development for our panel model, we attempted to construct a natural extension of the Lévesque-Dissanska model; however, their principal-plane model is defined from negative to positive view zenith angles.Because our BRF data are defined using only positive view zenith angle coordinates around the hemisphere, the different angular coordinate systems could not be merged without discontinuities in the BRFs, particularly for the diffuse scattering term.In our version of the model, therefore, the diffuse component was developed independently and is a sum of two terms: diffuse-power D P and diffuseforward D F .These have the form where Authorized licensed use limited to the terms of the applicable license agreement with IEEE.Restrictions apply.
with two free parameters α D1 and γ D1 associated with D p , five free parameters α D2 , β D2 , γ D2 , α D3 , β D3 , and γ D3 associated with D F , and three free parameters α R1 , β R1 , and γ R1 associated with the wavelength-dependent reddening factor R(λ), which differs in form from the one originally developed by Lévesque and Dissanska, as described later.The diffuse-power term is a consistent inverted bowl-shaped function, and the diffuse-forward term is a Gaussian with increasing scale factor, depending on the zenith angle and obeying a power law.The width of the Gaussian also increases according to a power law as a function of the input zenith angle.The reddening factor R(λ) in this expression accounts for more specular reflection and forward scattering at longer wavelengths [13].
We modeled forward scattering as a product of two unnormalized t-distributions and the reddening factor from ( 8).The t-distributions are functions of the relative azimuth angle and are centered on the forward part of the principal plane.For consistency in form, the heights and widths of the t-distributions also vary according to power laws as a function of the incident or view zenith angle, similar to the Gaussians in our diffuse-forward component.Our forward scattering term takes the form where +1 2 (10) with five free parameters α F1 , β F1 , α F2 , β F2 , and γ F2 .Although the exponential function used in the model of Lèvesque and Dissanska (Table I) is simpler in form than our model, we found that numerical convergence was difficult to obtain when applying an exponential function to data that includes an azimuthal dimension.The alternative Gaussian function that they recommend for data obtained at view zenith angles exceeding 70 • showed reasonable fits, but we found that a peaked function with a higher kurtosis was more appropriate.Because our diffuse model component differs substantially from that of Lévesque and Dissanska, this necessitated a different functional form.Therefore, in our model, we incorporated unnormalized t-distributions as the functional form to better model the full distribution of observations in azimuth and zenith.
While the specular component in our model incorporates a Gaussian function dependent on the relative zenith angle like that in Lévesque and Dissanska's model, our model also incorporates a second Gaussian distribution in the azimuthal dimension.The amplitude of our specular component varies according to a power law as a function of the incident zenith angle rather than the cosine dependence found in Lévesque and Dissanska's model; the Gaussian widths, however, are optimized constants for both models.We also include our reddening factor within our specular term.Our specular term takes the form with four free fitting parameters α S1 , β S1 , α S2 , and α S3 for our model.Like our forward-scattering term, we center the specular component on the forward part of the principal plane.
Finally, the backscatter component has a functional form that is similar to our specular component but without the reddening factor where α B1 , β B1 , γ B1 , α B2 , and α B3 are free parameters.

A. BRF Model Comparison With the Lévesque-Dissanska Model
In Table I, we compare the explicit forms of D, F, S, B, and R used in our BRF model and that of Lévesque and Dissanska, and we also provide the optimized numerical parameters of our model and their units.For the scattering functions developed by Lévesque and Dissanska [13], [14] we also have provided their numerical parameters within the equations.We note, however, that the units of the input angles for the functions in the Lévesque and Dissanska model are in degrees, while our model uses radians, so the units of their model parameters are also different.We have also revised the nomenclature developed for several scattering components by Lévesque and Dissanska.We rename their "deep" scattering to be "diffuse" scattering for our analysis, a term that semantically incorporates both single and multiple scattering as part of the diffuse reflectance.We also rename their "forward" scattering to be "specular" scattering and their "subsurface" scattering to be "forward" scattering for our analysis.
The diffuse and forward-scattering terms are functionally similar between the two models, although we note differences in behavior for the specular and backscattering components, which likely stem from differences in the development of our models.Specifically, for our model, the specular and backscattering components both increase in amplitude with respect to the incident zenith angle.In contrast, the specular and backscattering components from Lévesque and Dissanska decrease and remain constant, respectively.We also note that our reddening function is a multiplicative factor in forward-scattering terms, whereas their reddening function is an additive function with an explicit dependence on the polarization angle.

B. Fitting the BRF Model
Because the Spectralon panel serves as both the sample and its reflectance conversion target in our analysis, additional steps were required to convert the panel measurements into radiometrically accurate BRFs.At each incident zenith angle, all panel radiance measurements in the observing hemisphere were divided by each scan's nadir measurement to obtain a nadir-normalized reflectance r 0,data for panel radiance L r ( φ, θ e , θ i , λ) at each position on the observation hemisphere.The BRF is calculated using the nadir-normalized reflectance according to the following relationship where the normalization coefficient A(θ i , λ) is the same one appearing in (3).However, A(θ i , λ) is not known a priori because the reflectance model is required to estimate the normalization coefficient.Therefore, we fit our reflectance model comprised of D, F, S, B, R, and C, but without the normalization coefficient A, to r 0,data by minimizing the chi-squared χ 2 over all four variables φ, θ e , θ i , and λ where σ = σ ( φ, θ e , θ i , λ) is the estimated radiometric measurement uncertainty, that varies with all four variables at each point in the scan, and r 0,model is the nadir-normalized reflectance model that is mathematically equivalent to BRF × A(θ i , λ).The χ 2 was minimized by using the adaptive Nelder-Mead simplex procedure [36] implemented in Python's SciPy library to vary the 27 free parameters used in the BRF model.We provide the values of the optimized values of the free parameters in Table I.We determined initial guesses for the free parameters in the optimization step by manually varying the free parameters and approximately matching the form of the initial model to the data through visual inspection of successive residuals in juxtaposed plots.

C. Estimating the Normalization Coefficient A(θ i , λ)
Once we had determined the best-fit free parameters for the nadir-normalized reflectance model r 0,model , we then used the model to estimate the normalization coefficient A(θ i , λ) to convert both r 0,data and r 0,model into true BRF Authorized licensed use limited to the terms of the applicable license agreement with IEEE.Restrictions apply.
values, shown in Figs. 3 and 4, according to (14) using the relationship where the differential solid angle is d = sin (θ e ) d φdθ e .
In other words, the normalization coefficient A(θ i , λ) is the hemispherical average of r 0,model .However, r 0,model is an interpolative model, defined only within the range of view zenith angles measured and used to develop the model (θ e ∈ [0 • , 70 • ]).Because ( 16) requires view zenith angles up to 90 • , the BRF model, therefore, was extrapolated linearly for the larger view zenith angles from 70 • to 90 • instead of using values from the BRF model itself.Various orders of polynomials were fit to subsets of the data as a function of the view zenith angle, and reasonable fits produced values of A(θ i , λ) all within 0.5% of the coefficient values found in the linear extrapolation.Given these small differences, we chose the simpler linear extrapolation to avoid overfitting in estimating A(θ i , λ).In other words, to calculate the hemispherical average of the reflectance model, the model itself was used for view zenith angles from 0 • to 70 • , while the linear extrapolation of the model was used for view zenith angles from 70 • to 90 • .

D. Considerations for Model Fitting
Developing an empirical model for the panel BRF enables a user to calibrate the model to their specific panel using fewer measurements.To ensure the accuracy of the model fitting with fewer points, we recommend concentrating more measurements in areas where the gradient of the BRF, i.e., change with respect to viewing angles, is higher for a given sample and illumination geometry.For our analysis, several measures were taken to ensure numerical convergence of the χ 2 minimization to a physically reasonable model.We found that the Nelder-Mead simplex algorithm implemented by SciPy performed more reliably for our model when the range of adjustable free parameters was of a similar order of magnitude, as the procedure initially constructs the smallest possible simplex.To achieve this, we converted all input angles to radians and wavelengths to microns, and we imposed bounds on the free parameters to ensure that the model returned realistic, physical values.We developed various iterations of our model, and, during this process, we eliminated some free parameters and variables to simplify our model.To ensure the convergence of our procedure to a global minimum, we perturbed the initial guesses for our parameters and ran the optimization procedure numerous times; all instances converged to effectively produce the same form of the model.We also note that the first specular intercept term α S1 converged to a lower bound of 0.0, suggesting that this parameter could be eliminated, although we still present it in Table I to indicate that it was varied as a free parameter.For ease of use, on GitHub [19] we have provided our Python implementation of our model that also includes a lookup table of the optimized normalization coefficients that we obtained.

V. RESULTS
The BRF panel model, a subset of data for comparison, and the various scattering components of the model appear in Fig. 4 for an example wavelength (750 nm) represented in an array of polar colormaps.The first four columns show the optimized model components, which are described in Section IV and summarized in Table I, obtained at each of the illumination angles.The fifth column shows the resulting BRF model, and the sixth column illustrates the measured BRF data for which the model was optimized.

A. BRF Residuals
Fig. 5 shows the residuals of the BRF model in an array of polar colormaps in the same layout as Fig. 3.The residuals were calculated by subtracting the panel model BRFs from the measured BRF panel data.Distances between the dashed contour lines in Fig. 5 denote a change in reflectance of 0.01.Within most of the observing hemisphere for all incident zenith angles and wavelengths, the residuals show that the model is accurate to within 1%-2% reflectance, although we note residuals larger in magnitude at θ i = 70 • in the forward-most scattering regions.There, the model over-predicts the forward scattering on the left hemisphere at shorter wavelengths by about 5%, while the model under-predicts the forward scattering on the right hemisphere at more intermediate wavelengths by about 3%.At these large incident and view zenith angles, the panel BRF changes more rapidly with respect to view geometry compared to changes at smaller incident angles.Therefore, radiance uncertainty or pointing error has a more pronounced effect on model accuracy in these regions.

B. Errors From Adjacency Effects
We identify two potential sources of adjacency effects in our data that represent some small sources of error and uncertainty as revealed in the residual plots: 1) GRIT-T's moving carriage, arm, and head as it scans the sample; and 2) the backboard on GRIT-T that shields the ASD FR4 units from the scanning area, shown in Fig. 1.Although GRIT-T and the laboratory room were designed to be minimally reflective, the Spectralon panel itself intensifies all adjacency effects due to its high reflectivity.GRIT-T's movement during the scanning process might cause asymmetry in the panel BRFs, which may then manifest as small asymmetry in the residuals.This asymmetry, for example, can be seen at θ i = 20 • and θ i = 30 • in Fig. 5. Small regions of higher residuals of about 1% reflectance are seen on the left side ( φ ∼ 270 • ) of the hemisphere, which are not present on the right side ( φ ∼ 90 • ).To compensate for the potential of these types of adjacency effects from GRIT-T's movement, we scale all nadir measurements at each scan in azimuth to the same radiance level, along with the other measurements within the same azimuthal scan as the nadir measurement with the same scale factor.In addition, GRIT-T's backboard is coated with 3% reflectance Avian Black to minimize secondary scatter; however, a small amount of stray light may still be reflected back onto the sample area, primarily into the backscattering direction.The solid angle subtended by the backboard for a sample at a typical measurement height is approximately 0.19 sr, which is 0.19 sr/2π sr ≈ 3% of the solid angle of the observing hemisphere.Therefore, at most, only about 3% × 3% = 0.09% of the source irradiance could be reflected back onto the measurement area for a Lambertian sample, which would be a systematic but small source of secondary scatter.

C. Retrieving the Filling Factor of Spectralon
The microscopic structure of Spectralon is similar to that of densely compact sediment; therefore, in addition to developing our empirical model, we applied Hapke's radiative transfer model for semi-infinite granular media [37] to our BRF data.Using Hapke's model, we retrieved the filling factor (ϕ) [38] of our Spectralon panel.The filling factor is a parameter in the radiative transfer solutions that defines the fractional amount of volume occupied by a particle in any given volume of the medium.Hapke [37], [38] examined the effects of porosity on the extinction characteristics of a medium, showing that a nonlinear porosity function is a natural scale variable that appears in closed-form solutions to the radiative transfer equation.He found that this porosity function, K (ϕ), has the approximate form [37], [38] We used a modified procedure that inverted Hapke's solutions for our panel at two selected incident angles, following an approach previously used by us to retrieve porosity from multiangular remote hyperspectral data and imagery [24], [39], [40].The procedure imposes similarity conditions, among many, such as constraining the retrieval to converge to the same phase function for both incident angles.Using a subset of our BRF data at θ i = 30 • and θ i = 60 • , we retrieved a filling factor of ϕ ≈ 0.6147 for our Spectralon panel.In comparison, sediment samples from the northern ends of the Algodones Dunes, CA, USA, with a range of densities from ∼1.45 -1.69 g cm −3 , had retrieved filling factors of ∼ 0.49 -0.57[39].Datasheets from Labsphere state densities of 1.25 -1.5 g cm −3 for Spectralon [41].

VI. FIELD APPLICATION OF THE SPECTRALON MODEL
For notational brevity in this section, we represent a given state of the three angles of φ, θ e , and θ i with e,i herein.
Similarly, e represents the two angles of φ and θ e , while i represents the two angles of φ and θ i .If the Spectralon reference target is approximated as a Lambertian reflector, the reflectance spectrum r mn of pixel (m, n) in a hyperspectral image (HSI) can be calculated as for pixel radiance spectrum L mn (λ), panel radiance spectrum L p (λ), and the 8 where BRF( e,i , λ) is the panel BRF model from Section IV, DDRF(λ) is the diffuse-directional reflectance factor of the panel, and a and b are their respective scaling coefficients.We approximate the DDRF as Lambertian due to the diffuse illumination; that is, the DDRF is approximated to be independent of the view angle.Just as the BRF is a theoretical, idealized quantity, the HDRF is likewise an idealized quantity; HDRF measurements are more precisely described by the hemispherical-conical reflectance factor (HCRF) [21].We, however, use the generic terminology HDRF, applying the same convention to the DDRF.Equation (20) approximately models the skylight as isotropic illumination that is dependent only on the wavelength λ and assumes that the DDRF of the reference panel has been measured using an integrating sphere.The scaling coefficients a and b are field-dependent quantities for the same panel radiance spectrum L p ( e,i , λ) from ( 18) and ( 19) and a shadowed panel measurement L sh ( e , λ).
The scaling coefficient a represents the reflected radiance contribution from the direct sunlight component relative to the full hemispherical illumination, where a = 1 − b.Therefore, in addition to the panel BRF model and panel radiance spectrum L p ( e,i , λ), the panel DDRF(λ) and a shadowed panel measurement L sh ( e , λ) are required to calculate the HDRF in this manner.Ideally, shaded and unshaded reference panels are placed within the scene of interest and present in the imagery to acquire radiance measurements of L p and L sh ; this ensures simultaneous and similar measurement geometries for targets and reference panels.A single panel could also be shaded halfway, so both L p and L sh can be acquired from one panel, provided that the imagery has sufficient spatial resolution.If L p and L sh are measured from a different geometry than that of the imagery, such as with a separate spectrometer, the panel BRF model can also be applied to L p and L sh .We note that the HDRF, while a reflectance factor, depends on not only the reflectance properties of a target but also the illumination conditions of a given site; the HDRF is a field-dependent quantity.
If the angular distribution of the sky radiance L sky ( i , λ) is known, we can replace the integrating sphere measurement of DDRF with a more accurate estimate of the skylight component modeled with respect to view angles; measurements of L sky ( i , λ) can be acquired, for example, by GRIT-T using the second onboard spectrometer [22], although for the dataset described in this work, these measurements were not taken.In this case, also incorporating the panel BRF, the DDRF can be calculated explicitly with respect to sensor view angles as BRF( e,i ′ , λ) dω ′ (22) where

The primed variable ′
i represents φ ′ and θ ′ i , whereas e,i ′ represents φ ′ and θ ′ i , but θ e .Since the panel BRF model is limited to incident zenith angles between 10 • and 70 • , the BRF for incident angles between 0 • -10 • and 70 • -90 • can be estimated in a manner similar to the approach that we used to extrapolate the normalization coefficient A(θ i , λ) for our panel model.
In addition to the HDRF, the BRF of a specific HSI pixel at (m, n) can be computed as BRF mn ( e,i , λ) = L mn ( e,i , λ) − L mn,sh ( e , λ) L p ( e,i λ) − L sh ( e , λ) where BRF p is the panel BRF, and L mn,sh ( e , λ) is a shadowed radiance measurement of the target pixel or that of a shadowed pixel of the same composition and illumination as the target.The numerator of ( 23) can be expressed as where b mn = L mn,sh ( e , λ)/L mn ( e,i , λ).We note the similarity in the form of b mn with that of b from (21); they are both ratios of shadowed and unshadowed measurements of the target and panel, respectively.If the shadowed target radiance L mn,sh ( e , λ) is not acquired, b mn can be estimated as where b is from (21), and DDRF s (λ) is the diffuse-directional reflectance factor of a sample specimen that accurately represents the diffuse scattering within a target pixel.The DDRF s (λ) can be measured with an integrating sphere in a laboratory, for example, and used under the further approximation of isotropic sky illumination.Thus, we can estimate the shadowed target radiance as Authorized licensed use limited to the terms of the applicable license agreement with IEEE.Restrictions apply.

A. Measuring the BRFs for Field-Worn Panels
The reflectances of Spectralon panels used in the field degrade over time.Fig. 6 shows the BRF data of our pristine panel and two of our field panels for four wavelengths measured at θ i = 50 • , which is close to the solar zenith angle at the time of our field experiment.The differences in the BRFs between the pristine panel and the field panels are also computed and shown.Distances between contours in the BRFs denote a change in reflectance of 2%, whereas distances between dotted contours in the difference maps denote a change in reflectance of 1%.These field panels have been used in numerous field experiments over many years [22], [23], [24], [25], [26], [27], [28], [29], [30], [31], [32], altering their reflective properties.We therefore sought to use panel models applied separately for the field panels by replicating the procedures used for the pristine panel.The residuals computed from the field panel BRF data and models are also shown in Fig. 6 in the same color axes as Fig. 5.Although the majority Fig. 7. RGB images derived from the HSI imagery of a mudflat at the RIT Tait Preserve (43 • 8 ′ 29.53 ′′ N, 77 • 30 ′ 8.50 ′′ W) during the field experiment using three mast heights from the mast-mounted hyperspectral imaging system [42] and the georectified drone hyperspectral imagery (HSI).The latitude and longitude of the drone frame corners are provided, and panels are identified with their geometries labeled.The colored labels of the panels correspond to the colored spectra in Fig. 8. of the BRFs within the observing hemisphere are similar, we measured a significantly larger forward-scattering peak that is common to both field panels.Particulates in various field settings inevitably accumulate onto the panels over time, reducing the diffuseness of their reflectances.In Section VI-B, for the application of our model to the HSI acquired during our field experiment, we used the panel models derived for the field panels.

B. Applying Panel Models to HSI in a Field Setting
We acquired HSI under a clear sky from the mast-and drone-mounted imaging spectrometers on a mudflat during a field experiment at the Rochester Institute of Technology (RIT) Tait Preserve, Penfield, NY, USA, on October 5, 2022.The mast-based hyperspectral images were collected with a Headwall Micro-High Efficiency VNIR hyperspectral imaging spectrometer mounted on a General Dynamics pantilt unit [42].This imaging spectrometer was affixed to a BlueSky Mast system, which can elevate the spectrometer platform to various heights [40], [42], [43].Drone images were acquired with a Headwall Nano VNIR hyperspectral imaging spectrometer within a multisensor imaging payload mounted on a DJI Matrice 600 Pro drone platform [30], [32], [44], [45].Both unshadowed and shadowed hyperspectral measurements of Spectralon panels were acquired onsite with an ASD FR4 field spectrometer, while DDRF measurements of the various field panels were acquired in our laboratory with a small Labsphere integrating sphere configured with an ASD spectrometer.Fig. 7 shows RGB frames derived from VNIR hyperspectral images that were used to evaluate our application of the Spectralon BRF panel model to field imagery.Before imagery was acquired, each panel was leveled using two-axis levels and cleaned using dust blowers.These RGB frames were derived from a subset of HSI imagery consisting of images collected from three different mast heights with the mast-mounted hyperspectral imaging system [42] and an airborne drone hyperspectral image acquired from the RIT drone imaging platforms [30], [44].
The plots in Fig. 8 compare the panel models P( e,i , λ) and their respective panel calibration coefficients C(λ) from 400 to 900 nm.The effects of applying the different models to the measured panel radiances are also shown.The results for panel 1 appear in the left column plots, while the results for panel 2 appear in the right column.The first row shows All quantities are spectral.L p and δL p are the measured panel radiances and their radiometric uncertainties, respectively; Authorized licensed use limited to the terms of the applicable license agreement with IEEE.Restrictions apply.
P and δ P are the panel model from (20) and its radiometric uncertainty.The spectra acquired from the drone HSI data were saturated in the wavelengths not presented in the figure because the exposure time for the 12-bit Headwall Nano spectrometer had been optimized to a 50% reflectance target.In contrast, the Headwall Micro HE mast-based imaging spectrometer is a 16-bit system and acquires spectra of both highly and minimally reflective targets.As a result, data from the mast-mounted Micro HE imaging spectrometer are depicted as full VNIR spectra in red, green, and blue lines in the third and fourth rows of Fig. 8.In the same plots, orange lines indicate partial VNIR spectra from the drone Nano imaging spectrometer where the system was not saturated over the white Spectralon panel.Fig. 8 shows that the percentage differences between P and C range from −5% to 0% for geometries associated with panel 1, while differences range oppositely from 0% to almost 5% for geometries associated with panel 2. Equations ( 18), (19), and (23) show that the panel scaling coefficients and the panel models being applied are directly proportional to the reflectance quantity being calculated.In other words, errors would directly propagate into and have a systematic impact on reflectance calculations if the more complete panel reflectance models were not applied.Due to the proportionality, the magnitudes of these systematic errors would be scaled by the ratio of radiances between target and panel pixels.Brighter target pixels would be subject to these larger systematic errors, while dimmer panel pixels would cause larger uncertainty in the reflectance calculations.The bottom row of plots in Fig. 8 shows that the differences between scaled radiances using Labsphere's calibration coefficients and our panel model are greater than the radiometric uncertainty.
To better understand the sources of error, we note that (20) approximates the sky radiance as isotropic illumination.Due to scattering in the atmosphere, the sky radiance distribution is not only anisotropic but also polarized and significantly brighter in blue wavelengths [46], [47], [48].Incorporating better models for the sky radiance distribution and the polarimetry of the illumination and scene would likely improve harmonization of the various scaled radiances L p ; however, this also would require a polarized BRF model of Spectralon with full view coverage in the observing hemisphere.More accurate measurements of the diffuse sky radiance would also reduce sources of error in our experiment.For example, a rotating shadowband radiometer would more precisely measure L sh as opposed to our less precise manually staged shadowband measurements.In addition, the relative polarization biases of the Headwall Micro HE and Headwall Nano systems were not accounted for in this experiment, and importantly, the view geometries of the two systems were significantly different.Skylight radiance is not only polarized but increases in blue wavelengths due to Rayleigh scattering in the atmosphere.The discrepancy in scaled radiance measurements between the mast and drone imaging spectrometers is highest in the shorter wavelengths in Fig. 8.This is likely due to the differences in polarization biases affecting measurements at these wavelengths.] , covering the UV-VNIR-SWIR λ ∈ [350 , 2500 nm] .We measured hyperspectral BRF data with the GRIT-T hyperspectral goniometer system [22], [23], which incorporates ASD FR4 spectrometers, and in this experiment, with a 3 • fore-optic and an optical scrambler.We developed an improved empirical model consisting of various scattering terms and scaling functions, inspired by the Lévesque and Dissanska principal-plane model [13] and fit the model to measured hyperspectral BRF data of the pristine Spectralon panel.
Our improved model extends the earlier model [13] by explicitly characterizing the azimuthal dependence of the panel BRF and incorporating more sophisticated forms for the diffuse term D, forward-scattering term F, specular reflection term S, backward-scattering term B, a reddening factor R, and a normalization coefficient A that is dependent on the illumination geometry.We summarized the explicit forms of the terms and the free parameters of our model in Table I.We have also provided a Python implementation of the BRF model and a lookup table of the coefficients A(i, λ) on GitHub [19].
Residual plots show that our model replicates our measurements to within 1%-2% reflectance within most of the observing hemisphere.We outlined cases for how the Spectralon panel BRF model can be used to estimate the HDRF or BRF of a target pixel within HSI, and we demonstrated an application of the panel BRF model to HSI acquired during a field experiment at the RIT Tait Preserve.We also measured the BRFs of our field-worn panels and compared the differences in BRFs between our pristine panel and our field panels, using the BRFs of the field panels in our application to the field HSI.
The field panels were observed to have stronger forwardscattering peaks.Mast-and drone-based HSI data acquired from four different view geometries for each of the two field Spectralon panels were used to evaluate the effectiveness of our panel model in a field setting.The differences in scaled radiances ( L p ) for the Spectralon panel when normalizing by the panel model P versus the diffuse calibration coefficients C were significantly greater than the radiometric uncertainty (δL p ) for our measurement geometries in the field setting.Although these effects were significant, additional modeling that incorporates polarization effects [2] and the skylight distribution is required for more accurate modeling of the panel radiometry in future work.

Fig. 1 .
Fig. 1. (Left) GRIT-T viewing the pristine Spectralon panel (model: SRT-99-180) at nadir in the dark laboratory room with the external quartz-tungsten-halogen (QTH) light source illuminating the measurement area.(Top right) Close-up view of our stage lamp housing the QTH light source used to provide directional illumination onto the measurement area.(Bottom right) Context photograph showing GRIT-T, the illumination source, and one ASD spectrometer, with an olivine sand sample.

Fig. 2 .
Fig. 2. 1σ uncertainty derived from repeated nadir BRF measurements of the Spectralon panel for each of the seven incidence angles.

Fig. 3 .
Fig. 3. Unpolarized BRF data of the Spectralon panel shown in an array of polar colormaps.Within each circle, the polar angle with respect to the north denotes the relative azimuthal angle φ from 0 • to 360 • .The radius denotes the sensor view zenith angle θ e from 0 • to 70 • where each concentric light gray circle is a 20 • step.Distances between contour lines denote a change in reflectance of 2%, and the black circle masks the self-shaded backscattering region.The color axis ranges from 0.80 to 1.20, although the actual BRF values exceed the upper range for forward-scattering regions at larger incident zenith angles.(Rows: top to bottom) BRF data ordered by incident zenith angle θ i .(Columns: left to right) BRF data ordered by example wavelength λ.

Fig. 4 .
Fig. 4. (Columns: left to right) Spectralon panel BRF model components D P , D F , F, and S + B, the combined model, and the corresponding subset of measured BRF data used to optimize the model at 750 nm.The reddening factor is not applied to the two forward-scattering components, F and D F , in this figure.(Rows: top to bottom) The BRF model components and data are ordered by incident zenith angle θ i as in Fig. 3.The scales of the BRF model components are depicted by their respective color axes at the bottom of each column.The full BRF model and data share the long color axes on the right side, which is identical to the one in Fig.3.

Fig. 5 .
Fig. 5. Polar colormaps of the BRF residuals in the same layout as Fig. 3. (Rows: top to bottom) BRF residuals ordered by incident zenith angle θ i to be consistent.(Columns: left to right) BRF residuals are ordered by example wavelength λ.Residuals were calculated by subtracting the BRF model from measured BRF data.The color axis ranges between ±0.05, although some residual regions where θ i = 70 • exceed these values.Distances between dashed contour lines denote a change in residual reflectance of 0.01.

Fig. 6 .
Fig. 6.BRF data of the pristine panel and two field panels, the BRF differences between pristine and field panels, and BRF residuals computed from fitting our Spectralon model to the field panel BRFs for θ i = 50 • .(Top left columns) Polar BRF plots of the pristine Spectralon panel and the two field panels imaged in Fig. 7 within the indicated columns.Distances between contours denote a change of 2% in the BRF.(Top right columns) Percentage differences in the BRFs between the pristine panel and the two field panels.Distances between dotted contours denote a change of 1% in the BRF differences.(Bottom plots) BRF residuals from fitting our BRF model to the field panel BRFs in the same color axis as that of Fig. 5.

Fig. 8 .
Fig. 8. Spectra for the Spectralon panel models and their corresponding measured radiances for the three heights of the mast-mounted hyperspectral imaging system and the drone hyperspectral image.(First row) Labsphere's provided calibration coefficients of the Spectralon panels C(λ) and complete panel models P( φ, θ e , θ i , λ) for comparison.(Second row) Comparison of measured radiances of the two field Spectralon panels scaled by their pristine coefficients C (dashed lines) and P (solid lines) for each panel.(Third row) Error budget analysis comparing the difference of scaled panel radiances L p = (L p /C)−(L p /P) with the uncertainty δL p given by (27) across the spectrum from 400 to 900 nm.The gray area indicates where L p /δL p < 3.
VII. CONCLUSIONWe presented an empirical unpolarized comprehensive BRF model of a pristine optical-grade Spectralon panel with full azimuthal coverage φ ∈ [0 • , 360 • ] , a wide range of sensor view zenith angles θ e ∈ [0 • , 70 • ] , and a wide range of illumination incident zenith angles θ i ∈ [10 • , 70 • Our model and theirs still share dependencies on the view zenith angle θ e , the incident zenith angle θ i , and the wavelength of light λ.

TABLE I SPECTRALON
BRF MODEL COMPONENTS AND FREE PARAMETERS • /h spectrum of the panel C(λ).Alternatively, a more complete panel model P( e,i , λ) can be substituted for C(λ) to more accurately calculate the hemispherical-directional reflectance factor HDRF mn ( e,i , λ).