Atmospheric Correction with the Bayesian Empirical Line

Atmospheric correction of visible/infrared spectra traditionally involves either (1) physics-based methods using Radiative Transfer Models (RTMs), or (2) empirical methods using in situ measurements. Here a more general probabilistic formulation unifies the approaches and enables combined solutions. The technique is simple to implement and provides stable results from one or more reference spectra. This makes empirical corrections practical for large or remote environments where it is difficult to acquire coincident field data. First, we use a physics-based solution to define a prior distribution over reflectances and their correction coefficients. We then incorporate reference measurements via Bayesian inference, leading to a Maximum A Posteriori estimate which is generally more accurate than pure physics-based methods yet more stable than pure empirical methods. Gaussian assumptions enable a closed form solution based on Tikhonov regularization. We demonstrate performance in atmospheric simulations and historical data from the “Classic” Airborne Visible Infrared Imaging Spectrometer (AVIRIS-C) acquired during the HyspIRI mission preparatory campaign. © 2016 Optical Society of America OCIS codes: (010.0280) Remote sensing and sensors; (010.1285) Atmospheric correction; (300.6340) Spectroscopy, infrared; (300.6550) Spectroscopy, visible. References and links 1. G. Schaepman-Strub, M. E. Schaepman, T. H. Painter, S. Dangel, and J. V. Martonchik, “Reflectance quantities in optical remote sensingdefinitions and case studies, ” Remote Sens. Environ. 103(1), 27–42 (2006). 2. C. D. Mobley, “Estimation of the remote-sensing reflectance from above-surface measurements, ” Appl. Opt. 38(36),7442–7455 (1999). 3. D. R. Thompson, B. C. Gao, R. O. Green, D. A. Roberts, P. E. Dennison, and S. Lundeen, “Atmospheric correction for global mapping spectroscopy: ATREM advances for the HyspIRI preparatory campaign, ” Remote Sens. Environ. 167, 64–77 (2015). 4. R. Richter and D. Schläpfer, “Atmospheric/topographic correction for satellite imagery,” DLR report DLR-IB, 565–01 (2005). 5. M. W. Matthew, S. M. Adler-Golden, A. Berk, G. Felde, G. P. Anderson, D. Gorodetzky, S. Paswaters, and M. Shippert, “Atmospheric correction of spectral imagery: evaluation of the FLAASH algorithm with AVIRIS data,” Appl. Imag. Patt. Recog. Workshop, 157–163 (2002). #253097 Received 13 Nov 2015; revised 13 Jan 2016; accepted 15 Jan 2016; published 27 Jan 2016 © 2016 OSA 8 Feb 2016 | Vol. 24, No. 3 | DOI:10.1364/OE.24.002134 | OPTICS EXPRESS 2134 Or shoot me an email," David.r.thompson@jpl.nasa.gov" " "


Introduction
Remote visible/infrared reflectance spectroscopy provides unique insight into Earth's ecosystems and surface composition.Investigators measure radiance at many wavelengths, λ , and correct atmospheric interference to retrieve the apparent surface reflectance ρ(λ ), the Hemispherical Directional Reflectance Function [1], or the related quantity R rs (λ ), the Remote Sensing Reflectance [2].Typical atmospheric correction is based on simulations by Radiative Transfer Models (RTMs), an approach which is effective but can be sensitive to uncertainty in atmospheric gas and aerosol state.Consequently, empirical methods play a complementary role for difficult atmospheric conditions.The "empirical line" is one such method that uses in situ measurements of spectrallyinvariant surfaces to fit a linear relationship between sensor readings and reflectance.This requires visiting and measuring many distinct locations in the scene.Assuming invariant targets can be found, differences in remote and in situ perspectives can still cause discrepancies due to atmospheric and surface heterogeneity and the interaction of solar angle with self shading and non-Lambertian surfaces.Overcoming these discrepancies requires tedious measurements from many diverse locations.Anticipated global mapping applications [3] are particularly challenging due to large and inaccessible wilderness areas.With few coincident in situ reference spectra, the resulting corrections can be unstable and sensitive to confounding effects.
Here we present a more general probabilistic formulation to unify RTM and empirical approaches and enable combined solutions.This significantly improves stability when reference spectra are sparse, poorly conditioned or inaccurate, and makes empirical corrections practical for challenging field environments with few reference spectra.The RTM solution places data in a reflectance representation so that subsequent corrections obey a predictable distribution.This permits a Bayesian prior on empirical correction coefficients and calculation of the posterior mode using closed-form Tikhonov regularization.The article begins by reviewing standard atmospheric correction methods.We describe the proposed approach, and explore its behavior with radiative transfer simulations.Finally, we test its performance on a set of reference measurements acquired by NASA's "Classic" Airborne Visible Infrared Imaging Spectrometer (AVIRIS-C) during the HyspIRI mission airborne preparatory campaign.[3].

Background
Traditional atmospheric correction algorithms like ATCOR [4], FLAASH [5], and ATREM [6,7], exploit a relation between the surface reflectance ρ, the atmospheric path reflectance ρ a , and the top of atmosphere reflectance ρ 0 [8]: Table 1 describes our notation.Here F represents extra-terrestrial solar irradiance [9,10], ψ the solar zenith, T the transmission of gases and aerosols, and S the spherical albedo of the sky.One such relationship exists for each wavelength, though our notation omits this for clarity.Typically T , ρ a and S are calculated using a radiative transfer solver such as 6s [11,12,13] or DISORT [14].Some atmospheric and surface parameters may be retrieved on a per-pixel basis using spectral information [15,3].Any inaccuracy in models of gas absorption or calibration can leave residual errors in the reflectance estimate.These errors have been addressed by a range of different postprocessing methods.One approach, EFFORT [16], "polishes" the spectra using a generalized set of reference spectra to quantify and suppress the high frequency noise component of retrievals.Gao and Liu [17] provide an efficient alternative that reduces residuals relative to a cubic spline fit.The ATCOR atmosperhic correction package can classify pixels into plant, snow, soil or other surfaces, then apply appropriate band-by-band multiplicative factors to reduce the high frequency noise that discriminates the retrieval from the ideal spectrum.Similar methods based on multiplicative coefficients have been used throughout the HyspIRI preparatory campaign [3].These factors can smooth residual roughness in reflectance estimates, but do not typically compensate for broad spectral features from calibration or scattering effects.
In contrast to RTMs, the empirical line method [18] posits a direct relationship ρ ≈ x 1 + x 2 L, with free parameters x 1 and x 2 .If 1 − ρS ≈ 1 this has a natural physical interpretation; Eq. 1 reduces to ρ 0 = ρ a + T ρ, making the free parameters x 1 and x 2 proportional to the path radiance and transmittance respectively.As ρS grows, the parameters lose their physical interpretation, but the linear model can still be a simple and effective predictor.For notational simplicity we combine correction coefficients into a vector x = [x 1 , x 2 ] T and introduce a data matrix A with one row A i = [1, L i ] per reference location.Here L i represents the remote radiance measurement corresponding to the location i, with a corresponding in situ reflectance measurement t i .Thus, the predicted reflectance ρ is given by ρ = Ax.We group the target reflectance values in a vector t = [t 1 ,t 2 . ..]T with one element per location.The correction coefficients, written x EL , should minimize sum squared prediction error: The classical solution does not constrain the magnitude of the correction, which is appropriate because its numerical scaling is variable and scene-dependent.Some proposed hybrids combine RTMs and empirical methods.Teillet and Fedosejevs use an in-scene dark target to estimate the energy scattered into the beam [19].This estimates of ρ a , while an RTM calculates the other terms.Moran starts with the empirical line method, but then uses atmospheric radiative transfer models or dark targets to estimate the additive term x 1 [20].This leaves one free parameter, x 2 , to be fit by reference targets.However, in principle both in situ references and models provide some information about the whole system; a handful of in situ reference targets may be insufficient, but the empirical estimate will eventually outperform as more field measurements are added.An optimal solution would incorporate all information with appropriate weight to empirical and data driven methods based on the certainty of each source.

Approach
Bayesian inference provides a formal approach to combine RTM calculations and in situ measurements.The proposed approach uses an RTM, informed by known observing geometry and estimated atmospheric state, θ θ θ , to calculate an initial guess at each reference location i, ω i , in reflectance units.The resulting data matrix B consists of rows We posit correction factors are distributed about the identity transformation, given by a coefficient vector More generally, we take the linear coefficients x to have a Gaussian prior N (µ µ µ, Q −1 ) with mean µ µ µ, precision matrix Q and covariance matrix Q −1 .We posit measurement errors are distributed according to a zero-mean Gaussian N (0, P −1 ).Ascribing prior standard deviations η g , η o , and η m to gain, offset, and the in situ channelwise measurement noise respectively, we have: In practice errors are not necessarily independent in each channel, nor are they Gaussiandistributed.However, the assumption is convenient and proves adequate in the airborne experiments that follow.This leads to the following posterior reflectance: Here z 1 and z 2 are normalizing constants.This is equivalent to Generalized Tikhonov Regression [21], which seeks the solution x BEL minimizing the following objective function: Here x 2 Q signifies the norm x T Qx, and x 2 P is the norm x T Px.This formulation penalizes correction coefficients that depart from the RTM prior.It has a Bayesian interpretation as a Maximum A Posteriori (MAP) estimator, maximizing the posterior density of the Gaussian prior and Gaussian data likelihood terms.The well-known solution has closed form [22]: For the controlled experiments that follow, we ascribe a regularization factor δ to both gain and offset priors such that η o = η g = δ .The use can set δ with domain knowledge or (as in our experiments) formal cross validation.It acts as a regularization term which balances the influence of the RTM and in situ data.We will find that performance is insensitive to this parameter, so it is not necessary to estimate it perfectly.To derive P −1 , we simply estimate diagonal values η 2 m using the sample variance.This leads to a simple procedure, which is performed separately for each wavelength: 1. Acquire a set of in situ reflectance measurements t t t = [t 1 ,t 2 , . ..]T , one per reference location.
2. Using a RTM-based atmospheric correction, transform each reference location's radiance L i to an estimate of surface reflectance, ω i 3. Associate the remote and in-situ spectra, and construct a data matrix B with one row 4. Construct a 2 × 2 covariance matrix Q −1 comprised of small diagonal coefficients δ 2 ; this represents the certainty of atmospheric correction.
5. Construct a diagonal covariance matrix P −1 representing measurement variance for each remote spectrum.
6. Use Eq. ( 7) to find linear correction coefficients 7. Use the linear correction coefficients to correct all RTM-derived reflectance estimates The following section explores the method's performance and its sensitivity to δ .

Evaluation
We perform two experiments to evaluate the Bayesian approach.The first experiment uses simulations to compare its performance with status quo correction methods, and to explore its sensitivity to the prior.The second experiment demonstrates performance on an actual airborne data acquired by the AVIRIS-C instrument.

Simulations
We first characterize performance in simulation with 20 visible/infrared spectra from the USGS spectral library [23].These spectra span a diverse range of shapes and magnitudes, including urban surfaces such as asphalt and building materials as well as natural surfaces such as vegetation, open water, snow, and mixtures of vegetation and soil.We resample these spectra to the wavelengths of NASA's "Classic" Airborne Visible Infrared Imaging Spectrometer, AVIRIS-C [24], and simulate atmospheric interference using the DISORT radiative transfer model [14].This calculates absorption and scattering by a typical 20 layer midlatitude summer atmosphere, with gas optical cross sections and aerosol parameterizations from the LibRadTran suite [25].We simulated the conditions of an actual AVIRIS-C overflight which took place on 7 May 2014 at 20 km altitude, with a nadir-pointed instrument and solar elevation of 68 degrees.Aerosol properties were based on models of Shettle [26].The optical depth was 0.25, with a rural type aerosol in the boundary layer, background aerosol above 2 km, and spring-summer conditions.The atmosphere was a mid latitude summer model with the H 2 O column scaled to contain 2 cm precipitable water vapor.Our experiment considers a set of independent trials (or scenes) s, each containing many spectra j associated with top of atmosphere radiances L s j .To simulate error-prone observations, we perturb these radiances to yield noisy observations L s j .The perturbations incorporate a scene-wide random gain and offset of φ sg and φ so with standard deviations γ sg and γ so respectively.These represent systematic errors due to uncertainty in atmospheric constituents and instrument calibration that the Bayesian empirical line aims to overcome.They also account for modeling approximations such as neglected coupling between absorption and scattering.We use conservative values of 1%, typical of radiometric calibration accuracy for imaging spectrometers in the visible/infrared range [24].
We also apply a separate gain and offset (with standard deviations γ jg and γ jo respectively) to each spectrum.This signifies measurement noise; we use 1%, which is a conservative choice considering typical instrument SNR [24].The perturbed spectrum can be written: The relative performance ranking of different algorithms is insensitive to the magnitude of this noise.After introducing errors, we simulate the correction process of Eq. ( 1) and applied the Bayesian empirical line adjustment.We select 1 to 5 in situ reference spectra as the training set, and calculate prediction accuracy on the held-out remainder based on the Root Mean Squared Error (RMSE) while ignoring water absorption bands.For each number of reference spectra, the different possible combinations of training and testing samples define a mean and standard deviation for performance.
The regularization δ enables a continuum of possibilities ranging from a completely constrained correction, equivalent to the original RTM solution, to a completely unconstrained correction, equivalent to a pure empirical method.Figure 1 illustrates the effect of changing the regularization parameters to favor RTM or empirical behavior.As the number of reference spectra increase, the position of the global error minimum shifts to the right.In other words, the certainty of the empirical correction improves and larger values of δ became optimal.Additionally, the position of the global error minimum shifts downward, but with diminishing returns for large numbers of in situ reference spectra.We note that there is always significant benefit to incorporating a model prediction: some degree of Bayesian regularization always outperforms the pure extremes.In other words, even single reference target or an error-prone atmospheric model provides some beneficial information.Another important result is that the regularized empirical line is stable over a wide parameter range.Performance degrades appreciably only when values of δ depart from the optimum by a factor of five or more.This suggests that the method can glean benefits even when reliable noise and error predictions are not available.
Figure 2 compares the relative performance of the proposed approach, a pure empirical line, spectral polishing [16], and the hybrid correction strategy of Moran [20].We use 1σ perturbations at 1% of the total signal level.The graph shows 1σ error bars, with arrows signifying intervals that extend outside the plotted area.When all accuracies lie outside the plotted area, we signify this with delta symbols at the top of the plot.The classical empirical line was numerically stable after three training spectra, and approached the optimum after five.The refined empirical line approaches asymptotic error with about three training targets.The Bayesian approach, exploiting both model and in situ information, outperforms all alternatives while reducing the standard deviation for all training set sizes.It provides a non-degenerate solution from a single training spectrum.

Airborne data
Our second evaluation uses airborne data: the ATREM-derived HyspIRI preparatory reflectance product [3].This algorithm typifies the RTM approach.We consider multiple reference targets imaged by AVIRIS-C during flights over California during the HyspIRI mission preparatory campaign [3] (Fig. 3).Eight terrestrial reference targets fall within a single long flightline from the 2013 data collection year, the "Soda Straw" spanning several degrees of latitude across the state of California.The targets include a range of light and dark surfaces at varying altitudes (Fig. 4).The reflectance profiles are measured in situ using a Visible Shortwave Infrared (VSWIR) field spectrometer and a reference spectralon panel, following the procedure outlined  in [3].
Figure 5 shows a typical result applying the Bayesian empirical line strategy to a held-out test spectrum.As expected, the reflectance shape and magnitude better match the ground reference data.We also compare performance against other correction methods by fitting these factors using subsets of 1 to 5 training targets, and then evaluate performance using the remainder.The result appears in Fig. 6 with error bars indicating 1σ performance for all cross validation trials.The accuracy profiles and relative performance rankings of the algorithms support the simulation results.Without exception, all methods' performance improves monotonically a more in situ measurements are available.Overall, the Bayesian method provides the best stability and fastest convergence.

Discussion and conclusions
This article synthesizes two traditional approaches for atmospheric correction of visible/infrared spectra: the empirical line and RTM-based methods.Bayesian inference provides a rigorous formal framework for incorporating both information sources, leading to more stable results.In our experiments, it approaches optimal accuracy with only one or two in situ measurements.This capability can lower the barrier to entry for a broader range of scientific, commercial, and public policy spectroscopy applications.Specifically, it makes accurate empirical line corrections practical for large wilderness scenes where it can be costly or impractical to obtain high quality coincident field data.
Future work could investigate alternative approaches to synthesize multiple data sources.Here we treat the RTM and in situ data collection as two separate processes.But if the in situ references did not exhibit significant systematic errors, each reference might instead be incorporated independently through multiple Bayesian updates.This would emphasize the measurement data and possibly simplify selection of regularization parameters (that would relate more directly to instrument SNR).Multiple RTM solutions from different models might be incorporated in a similar fashion, resulting in an ensemble estimator.Such modeling issues have a long history in data fusion research, and imaging spectroscopy on global scales could benefit from that foundation.
The probabilistic method points to the possibility of treating physical atmospheric parameters as random variables rather than targets of a fitting procedure.In this approach, pursued by researchers like Frouin and Pelletier [27], the atmospheric correction estimates the distribution using conventional or Bayesian inference.The ability to characterize the full posterior distribution over atmospheres, and propagate this information into uncertainties about the atmospheric retrieval, holds great potential to improve the quality and interpretability of atmospheric correction.

Fig. 1 .
Fig.1.Performance as a function of the regularization factor δ and number of reference spectra.The empirical solution is favorable as more reference spectra become available.

Fig. 2 .
Fig. 2. DISORT simulations compare multiple correction methods: the traditional empirical line (EL); the refined empirical line of [20] (REL); spectral polishing (SP); and the Bayesian empirical line correction (BEL).Arrowheads indicate 1-σ error bars and/or performance means that lie outside the chart.

Fig. 3 .Fig. 4 .Fig. 5 .
Fig. 3.The reference spectra in our AVIRIS-C study lie in a large flightlines transecting part of the US state of California.

Table 1 .
Mathematical notation BEL Vector of linear correction coefficients from Bayesian Empirical Line µ µ µ Mean of correction coefficient prior distribution, typically the vector [0,1] P Inverse covariance of in situ data, a 2 × 2 matrix Q Inverse covariance of correction factors, a 2 × 2 matrix γ jg , γ jo Standard deviations of perturbing gain and offset, independent in each spectrum γ sg , γ so Standard deviations of perturbing gain and offset, scene-wide η g , η o Standard deviations of prior gain and offset η m Standard deviation of measurement noise δ Regularization factor ascribed to η g and η o