A Fast Smoothing Algorithm for Post-Processing of Surface Reflectance Spectra Retrieved from Airborne Imaging Spectrometer Data

Surface reflectance spectra retrieved from remotely sensed hyperspectral imaging data using radiative transfer models often contain residual atmospheric absorption and scattering effects. The reflectance spectra may also contain minor artifacts due to errors in radiometric and spectral calibrations. We have developed a fast smoothing technique for post-processing of retrieved surface reflectance spectra. In the present spectral smoothing technique, model-derived reflectance spectra are first fit using moving filters derived with a cubic spline smoothing algorithm. A common gain curve, which contains minor artifacts in the model-derived reflectance spectra, is then derived. This gain curve is finally applied to all of the reflectance spectra in a scene to obtain the spectrally smoothed surface reflectance spectra. Results from analysis of hyperspectral imaging data collected with the Airborne Visible/Infrared Imaging Spectrometer (AVIRIS) data are given. Comparisons between the smoothed spectra and those derived with the empirical line method are also presented.


Introduction
Since the mid-1980s, hyperspectral imaging data have been collected with different types of imaging spectrometers from aircraft and satellite platforms. Because solar radiation along the sun-surface-sensor path in the 0.4−2.5 μm visible and near-IR spectral regions is subject to absorption and scattering by atmospheric gases and aerosols, hyperspectral imaging data contain atmospheric effects. In order to use hyperspectral imaging data for quantitative remote sensing of land surfaces and ocean color, the atmospheric effects must be removed. There are now a number of model-based atmospheric correction algorithms for retrieving surface reflectances from hyperspectral imaging data. These include, but are not limited to, the ATmosphere REMoval algorithm (ATREM) [1,2], the High-accuracy ATmospheric Correction for Hyperspectral Data (HATCH) [3], the Atmosphere CORrection Now (ACORN) [4], the Fast Line-of-sight Atmospheric Analysis of Spectral Hypercubes (FLAASH) [5], the Imaging Spectrometer Data Analysis System (ISDAS) [6], and a series of Atmospheric and Topographic Correction (ATCOR) codes [7,8].
The surface reflectance spectra retrieved with radiative transfer models often contain residual atmospheric absorption and scattering effects. The reflectance spectra can also contain artifacts due to errors in radiometric and spectral calibrations. Although models have been improving with time, they are not yet at the level where all artifacts are smaller than sensor noise. Figure 1 shows an example of a reflectance spectrum derived with ATREM from AVIRIS data [9] acquired over Cuprite (NV, USA) in June, 1995. Due to small errors in assumed wavelengths and errors in line parameters compiled on the HITRAN database [10], small spikes (particularly near the centers of the 0.94-and 1.14-µm water vapor bands) are present in this spectrum. These spikes have distracted geologists who are interested in studying surface mineral features. Clark et al. [11] pioneered a hybrid approach for the derivation of laboratory-quality surface reflectance spectra from AVIRIS data. They used a combination of ATREM and field spectral measurements over a single ground calibration site. In this case the use of ATREM allows improved atmospheric corrections at elevations that are different from the calibration site, and the ground calibration removes residual errors commonly associated with sensor artifacts and radiative transfer models. In many situations, researchers do not have any field-measured reflectance spectra for suppressing residual errors. Boardman [12] first developed the Empirical Flat Field Optimal Reflectance Transformation (EFFORT) method to suppress residual errors in ATREM-derived surface reflectance spectra without the need for field-measured reflectance spectra. In this method, the complete spectrum for each pixel in the 0.4-2.5 μm range is fitted with a low order polynomial. It is sometimes difficult to make reasonable matches over the entire (or "global") spectral range. In this article, we describe another technique, which fits spectra "locally" in the spectral domain using moving filters derived with a cubic spline smoothing algorithm, for quick post processing of ATREM-derived reflectance spectra from imaging spectrometer data. Results from analysis of AVIRIS data acquired over the Cuprite mining district in Nevada in June of 1995 and over Ivanpah in California in April of 2010 are given. Comparisons between the smoothed spectra and those derived with the empirical line method are also presented.

Methodology
In order to describe the smoothing technique, we first describe the commonly used cubic spline "fitting" technique, then we describe the cubic spline "smoothing" technique.

Cubic Spline Fitting
The cubic spline fitting technique is a powerful numerical method and has been widely used in engineering and scientific computing. For example, Numerical Recipes [13] provides standard subroutines, using cubic spline fitting method, for interpolating data between points. In order to describe mathematically the cubic spline fitting technique, we consider an interval α ≤ x ≤ b, and subdivide it by a mesh of points corresponding to the location of the data at α = X 0 < X 1 <…< X j−1 < X j …< X J = b. An associated set of the observed data is prescribed by y o , y 1 ,..., y j ,..., y J . We seek an interpolating function h(x), which is defined in the interval a,b [ ]. Its first and second derivatives are continuous on a,b [ ] and it coincides with a cubic polynomial in each subinterval X j −1 ≤ x ≤ X j , and satisfies the relationship h j = h X j ( ) = y j . Figure 2 illustrates the function h(x). As adapted from Ahlberg [14], the function h(x) in the interval X j −1 ≤ x ≤ X j can be expressed as (for convenience, we assume the problem of equally spaced samples with a step size of Δ ): : : : where: {s j }, the spline coefficients, can be interpreted as the normalized second derivatives. (2)

Figure 2. An illustration of the interpolating function h(x).
The polynomials (2) in adjacent segments are continuous at the knots: The first derivative is continuous at the knot provided that: The second derivative is continuous at the knots: The polynomials (2) are determined by specification of {s j }. The selection of these spline coefficients can involve any number of imposed weak constraints that characterize the spline fitting. One of the constraints is the minimization of the second derivative. Because: it follows that the quantity s j −1 ] is to be minimized in any kind of variational selections of {s j }. The simplest quadratic form to minimize is However, this is not enough to guarantee continuity of the derivatives at the knots. A method to incorporate the condition: must be found. This is done by introducing Lagrangian multipliers. The simple spline formulation for the minimization is: where λ j s are the Lagrangian multipliers. These conditions are exactly satisfied upon completion of the minimization so that zeros are in effect added to the quantity to be minimized. The procedure of solving those λ j s, therefore, the spline coefficients {s j }, and the interpolating splines {h j }, is similar to that of spline smoothing to be described in the next section.

Cubic Spline Smoothing
In the spline fitting technique described above, the {h j } are taken to represent errorless data or observations, and the spline passes each point y j . However, there can be circumstances that the observations are contaminated and unwanted noise is present. For example, in our case raw spectra exhibit coherent saw-tooth like "noise". Under these circumstances, the data integrity condition should be relaxed. This can be done by adding a weak constraint term, (8), where y j is the observed data, and only a "best fit" should be sought. The smoothed spline { h j } does not necessarily pass original observed data { y j }, unlike the case of the spline fitting. An appropriate discrepancy sum can be formed as: where τ 2 is an adjustable weighting factor. As it increases, the tension of the spline smoothing increases, i.e., the curve "flattens out". On the other hand, as it decreases, the observed data are reproduced more closely at the expense of increased curvature. The variations on the spline coefficients are tabulated as: : : : : : The variations on the multipliers lead to: Since the spline does not pass the data { y j }, the {h j } are no longer fixed; their variations are listed below: : : : Combining terms in (12), we have: Each of the combinations s j −1 + 4s j + s j+1 [ ] ; j = 1,...,(J − 1) in (10) can be replaced by their equivalents from (11) to obtain the following equations: The (14) can be replaced by the groupings in (13): : : : : : where: The {λ j } are then found as solutions of: where A is the pentadiagonal matrix Here the matrix elements are constants and given by: and: Given the {λ j }, the {h j } can be determined directly from (13), and the spline coefficients {s j } can be found as a tridiagonal solution of (10).

Procedures for Post Processing of ATREM-Derived Reflectance Data
The small spikes, as seen in Figure 1, are systematically present in all spectra in an ATREM output data cube (2-d spatial plus 1-d spectral). We hope to make mild "gain" adjustments to remove these small spikes during the post processing of the ATREM data cube. Specifically, we hope to find a gain function, g( λ ), which contains all the weak spikes and which has a mean value close to 1. The multiplication of g(λ ) to the ATREM output spectra should allow the removal of the systematic small spikes. Several steps are involved in the post processing of an ATREM output data cube. They are: (a) The cubic spline smoothing technique described in Section 2 is applied to each of the spectra in the ATREM data cube. As a result, an intermediate "smoothed" data cube is produced. Because the cubic spline smoothing technique fits a spectrum "locally" in the spectral domain, most of the smoothed spectra at this stage match quite well with the ATREM spectra.
If the spectra were fit with low order Legendre polynomials "globally", only a minor fraction of the smoothed spectra would match well with the ATREM spectra. (b) The average reflectance, ρ avg , is calculated for each of the spectra in the ATREM output data cube.
(c) For each pixel, the standard deviation, σ , between the ATREM spectrum and the "smoothed" spectrum is calculated.
(d) For an AVIRIS scene, a scatter plot of σ / ρ avg vs. ρ avg is made. Pixels with σ / ρ avg values in the lower twenty percentile are identified. (e) For each of the pixels identified in Step d, a ratio spectrum ("smoothed" spectrum/ATREM spectrum) is calculated. The desired gain spectrum, g(λ ), is obtained by averaging all the ratio spectra. Figure 3 shows an example of a gain spectrum, which contains a number of weak spikes in the 0.4-2.5 µm spectral region. (f) Each of the spectra in the ATREM output data cube is multiplied by the gain spectrum to obtain the "final" smoothed data cube. Our algorithm for smoothing the ATREM output data cube is fast. It takes approximately 30 s on a Mac computer with a 2.66 GHz Quad-Core Intel Xeon processor to process one complete data cube with a dimension of 614 samples, 972 lines, and 224 bands.

Sample Results
The cubic spline smoothing algorithm described above was implemented on an AVIRIS server computer at the NASA Jet Propulsion Laboratory for post-processing large volumes of ATREM-derived reflectance data cubes from AVIRIS radiances acquired during major field experiments. Below we present results from one set of AVIRIS data acquired over the Cuprite Mining District in Nevada in June, 1995 and another over Ivanpah playa in California in April, 2010. Figure 4 shows a comparison among an ATREM reflectance spectrum (lower line) over a single pixel within the 1995 AVIRIS Cuprite scene , the smoothed spectrum (middle line), and the reflectance spectrum obtained with the well-known empirical line method (upper curve) [15]. For clarity, the spectra in Figure 4 are vertically displaced. The general shapes of these spectra in the 0.4-1.26 µm, 1.5-1.75 µm, and 2.0-2.5 µm wavelength intervals are very similar. Major mineral features in the 2.0-2.5 µm region are seen in all the spectra. The un-smoothed ATREM spectrum has quite a few weak spikes. These spikes are largely removed in the smoothed spectrum. The spectrum derived with the empirical line method shows weak inverse water vapor features near 0.94 and 1.14 µm. This indicates that the method results in a slight over-correction of atmospheric water vapor absorption effects for this pixel.  Figure 5A shows six ATREM reflectance spectra (vertically displaced for clarity). These spectra have distinct mineral absorption features in the 2.0-2.5 µm spectral region. Weak spikes (for example near 1.14 µm) are systematically present in all the spectra. Figure 5B shows the corresponding smoothed spectra, which look very similar to laboratory-measured reflectance spectra, particularly in the 2.0-2.5 µm spectral region. Weak spikes are all removed. A broad iron feature near 0.9 µm is seen nicely in one spectrum-the 4th spectrum from top. Figure 5C shows six spectra derived from the AVIRIS data with the empirical line method. Mineral features in the 2.0-2.5 µm region are recovered quite well with this method. However, water vapor features in the 0.94 and 1.14 µm regions are either over-or under-corrected. The broad iron feature in the 4th spectrum from the top is not clearly seen due to the over-correction of atmospheric water vapor absorption effects. By comparing Figured 5A-C, it is seen that major mineral features are preserved after the spectral smoothing.   Figure 6B is the ATREM-derived surface reflectance spectrum over a soil pixel located at the center of the red box in Figure 6A. The curve is not spectrally smooth, particularly in the 0.86-1.20 μm range, due to residual errors in the ATREM atmospheric correction process. The solid line (vertically displaced for clarity) in Figure 6B is the smoothed spectrum. The spectrum in the 0.86-1.20 μm range becomes much smoother after the application of the cubic spline smoothing algorithm. By comparing the two curves in Figure 6B, it is also seen that major mineral absorption features centered near 2.20 and 2.34 μm are preserved after spectral smoothing, and no artificial absorption features in the entire 0.4−2.5 μm spectral range are introduced by the smoothing algorithm.   Figure 6, except for a spectrum over green vegetation covered area. The general shapes of the spectrum in the 0.5−2.5 μm range after smoothing (solid line in Figure 7B) are similar to those of green vegetation reflectance spectra measured in laboratories.   Figure 8A is similar to Figure 6A. Figure 8B shows comparisons among ATREM-derived surface reflectance spectrum (dotted line) over an Ivanpah playa pixel located at the center of the red box in Figure 8A, the smoothed spectrum (solid line and vertically displaced by 0.075 in reflectance unit), and a field-measured spectrum (dash-dotted line and vertically displaced by 0.15 in reflectance unit). A major mineral feature centered at 2.2 μm is seen in all the spectra.

Ivanpah, California
In order to quantify improvements in smoothness, first derivatives for all the spectra in the Ivanpah scene before and after application of the smoothing algorithm were calculated. After smoothing, the average absolute spectral derivative for all the bands, excluding those bands centered near 1.38-and 1.88-μm strong water vapor band absorption regions, over the scene decreased by 14%. The magnitude of the decrease is larger for a number of bands, including some bands located within the 0.94-and 1.14-μm water vapor band absorption regions where larger spectral residuals are observed (see Figure 7B). The scene-averaged spectral derivative for a band centered near 1.11 μm decreased by 20%. The decrease in spectral derivatives demonstrates the improvement in smoothness after applications of the cubic spline smoothing algorithm.

Discussion and Summary
During our development of the cubic spline smoothing technique for post-processing of surface reflectance spectra retrieved from AVIRIS data, we also tried another well-known filter, i.e., the Savitzky-Golay (SG) filter [16]. The use of SG filter didn't produce satisfactory results. We observed that the peak positions of spectral features after application of the SG filter were not preserved. Therefore, we decided not to use the SG filter for our spectral smoothing purposes.
In summary, we have described a technique, which fits spectra "locally" in the spectral domain based on cubic spline smoothing, for quick post processing of apparent reflectance spectra derived from AVIRIS data using the ATREM code. Results from analysis of AVIRIS data acquired over Cuprite mining district in Nevada in June of 1995 and over Ivanpah in California in April of 2010 are presented. Very good agreement between our results and those of empirical line method in the 2.0-2.5 µm spectral region is obtained. It is expected that the use of ATREM code for retrieving surface reflectance spectra from AVIRIS data plus the application of additional spectral smoothing should yield high quality surface reflectance spectra comparable with those of reflectance spectra measured in laboratory conditions.