Next Article in Journal
Detection, Localization and Classification of Multiple Mechanized Ocean Vessels over Continental-Shelf Scale Regions with Passive Ocean Acoustic Waveguide Remote Sensing
Next Article in Special Issue
Analyzing Performances of Different Atmospheric Correction Techniques for Landsat 8: Application for Coastal Remote Sensing
Previous Article in Journal
An Improved Algorithm for Discriminating Soil Freezing and Thawing Using AMSR-E and AMSR2 Soil Moisture Products
Previous Article in Special Issue
Inland Water Atmospheric Correction Based on Turbidity Classification Using OLCI and SLSTR Synergistic Observations
 
 
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

Fast Atmospheric Correction Method for Hyperspectral Data

by
Leonid V. Katkovsky
1,*,
Anton O. Martinov
1,
Volha A. Siliuk
1,
Dimitry A. Ivanov
1 and
Alexander A. Kokhanovsky
2
1
A. N. Sevchenko Research Institute of Applied Physical Problems, Belarus State University, BY-220045 Minsk, Belarus
2
VITROCISET, Bratustrasse 7, D-64293 Darmstadt, Germany
*
Author to whom correspondence should be addressed.
Remote Sens. 2018, 10(11), 1698; https://doi.org/10.3390/rs10111698
Submission received: 6 September 2018 / Revised: 23 October 2018 / Accepted: 25 October 2018 / Published: 28 October 2018
(This article belongs to the Special Issue Atmospheric Correction of Remote Sensing Data)

Abstract

:
Atmospheric correction is a necessary step in processing data recorded by spaceborne sensors for cloudless atmosphere, primarily in the visible and near-IR spectral range. In this paper we present a fast and sufficiently accurate method of atmospheric correction based on the analytical solutions of radiative transfer equation (RTE). The proposed analytical equations can be used to calculate the spectrum of outgoing radiation at the top boundary of the cloudless atmosphere. The solution of the inverse problem for finding unknown parameters of the model is carried out by the method of non-linear least squares (Levenberg-Marquardt algorithm) for an individual selected pixel of the image, taking into account the adjacency effects. Using the found parameters of the atmosphere and the average surface reflectance, and also assuming homogeneity of the atmosphere within a certain area of the hyperspectral image (or within the whole frame), the spectral reflectance at the Earth’s surface is calculated for all other pixels. It is essential that the procedure of the numerical simulation using non-linear least squares is based on the analytical solution of the direct transfer problem. This enables fast solution of the inverse problem in a very short calculation time. Testing of the method has been performed using the synthetic outgoing radiation spectra at the top of atmosphere, obtained from the LibRadTran code. In addition, we have used the spectra measured by the Hyperion. A comparison with the results of atmospheric correction in module FLAASH of ENVI package has been performed. Finally, to validate data obtained by our method, a comparative analysis with ground-based measurements of the Radiometric Calibration Network (RadCalNet) was carried out.

1. Introduction

The main task of the atmospheric correction is to find a relation between the surface spectral reflectance and the spectral radiance at the top of atmosphere (TOA), measured by satellite sensors. Atmospheric correction methods can be divided into empirical and based on radiative transfer models [1]. The first class of methods is based on the use of a-priori information about the object with which one can perform atmospheric correction without detailed radiation transfer modeling. The second class of methods is more accurate, but, in turn, requires more complex and time-consuming calculations.
There are a number of radiative transfer codes for the atmospheric correction: ACORN (Atmospheric CORrection Now) [2], based on MODTRAN-4; ATCOR (ATmospheric CORrection) [3], ERDAS Imagine; ATREM (ATmospheric REMoval), [4]; FLAASH (Fast Line-of-sight Atmospheric Analysis of Spectral Hypercubes) [5], RSI ENVI; HATCH (High-accuracy ATmospheric Correction for Hyperspectral data) [6,7], improved ATREM code; Tafkaa [8], based on ATREM. Most of these codes are designed for the specific satellite imaging systems, a certain spectral range, a set of spectral bands, spectral and spatial resolutions. The disadvantages of many techniques is the use of complex algorithms of the RTE solution [9] or pre-calculated look-up-tables, followed by interpolation [5], that requires a considerable time, or is unsatisfactory in accuracy. Generally, at the same time, they require a certain a-priori data about the parameters of the atmosphere. In [10] a fast method for reconstructing the aerosol optical thickness of the atmosphere is proposed, based on a combination of analytical and numerical solutions for various scattering layers positioned in atmosphere, which can also be successfully used to obtain the reflectance spectra of the underlying surface.
A more detailed description of the existing methods and software is given by Katkovsky [11]. A fast method of atmospheric correction presented in this paper does not require complex calculations and usage of a-priori data.

2. Atmospheric Model

Our atmospheric model includes the following processes:
  • molecular (Rayleigh) scattering,
  • aerosol absorption and scattering,
  • absorption by water vapor, oxygen and ozone.
The development of the optical model of the atmosphere is a necessary step in solving any inverse problem. Atmospheric models used in atmospheric radiative transfer studies usually account for the latitudinal zones, different seasons and vertical profiles of the scattering layers [12]. However, our analysis of the numerous spectral radiance calculations for various solar/observation angles and atmospheric models has shown that the spectrum of outgoing radiation [13] is weakly dependent on the vertical stratification of the atmosphere. A similar conclusion has been reached in [10], in relation to the stratification of the lower layer of the troposphere. Essential parameters used by us for the calculation of the spectrum of outgoing radiation are:
  • measurement geometry,
  • the spectral land surface reflectance,
  • vertical optical thickness of molecular scattering,
  • aerosol optical thickness,
  • scattering phase function parameter (average cosine of scattering angle),
  • single scattering albedo,
  • integrated content of water vapor (in the column of the atmosphere), oxygen, and ozone.
The spectra of outgoing radiation (outside the absorption bands of atmospheric gases) are most sensitive to variation of (in descending order of influence) surface reflectance, total vertical optical thickness and the average cosine of the scattering angle.
In particular, we have used the parameters of the optical model of the cloudless atmosphere discussed below.
(1)
Vertical optical thickness of the atmosphere at the wavelength λ (excluding optical thickness for gaseous absorption) is defined as
τ λ = τ sca , m + τ sca , a + τ abs , a ,
where τ sca , m is optical thickness of molecular (Rayleigh) scattering; τ sca , a is optical thickness of aerosol scattering; τ abs , a is optical thickness of aerosol absorption.
(2)
Single scattering albedo (quantum survival probability) ω λ is calculated as follows:
ω λ = τ sca , m + τ sca , a / τ sca , m + τ sca , a + τ abs , a .
(3)
The vertical molecular optical thickness is determined by the atmospheric model (season and location) [14]
τ sca , m = F P s , T s λ ( B + C λ + D / λ ) T s T 0 P 0 P s .
Constants in Equation (3) are the following: the pressure P s and the temperature T s at the Earth’s surface are taken in the selected atmospheric model (for example, Midlatitude Summer); P 0 and T 0 are real values of pressure and temperature, accordingly, at the Earth’s surface at the measurement time; constants B, C, D are independent on atmospheric models. We have used six types of atmospheric models, values of the coefficients for which are given in Table 1.
The spectral dependence of the aerosol optical thicknesses introduced in the atmospheric model is approximated by the power law function:
τ sca , a = τ sca , a 0 ( λ 0 / λ ) β ,
where τ sca , a 0 is the corresponding optical thickness at the reference wavelength λ 0 , β is Ångstroem exponent. The value of τ abs , a is supposed to be not dependent on the wavelength (it is usually small as compared to τ sca , a ).
(4)
The spectral reflectance of the underlying surface ρ λ is considered to be Lambertian (isotropic).
(5)
The total scattering phase function is given as the weighted-average function of the Rayleigh x m ( γ ) and x a ( γ ) the aerosol scattering phase functions as follows:
x ( γ ) = x m ( γ ) τ sca , m / τ sca , m + τ sca , a + x a ( γ ) τ sca , a / τ sca , m + τ sca , a ,
where aerosol scattering phase function is approximated by Henyey-Greenstein function. In particular, it follows:
x m ( γ ) = 3 4 ( 1 + γ 2 ) x a ( γ ) = 1 g a 2 / 1 + g a 2 2 g a γ 3 / 2 γ = μ μ 0 + ( 1 μ 2 ) ( 1 μ 0 2 ) cos φ ,
with an average cosine of scattering angle, g a , which is assumed to be independent on the wavelength. However, by means of Equation (5), the spectral dependence of the total scattering phase function on the wavelength is introduced. In Equation (6), γ is cosine of the scattering angle, the Sun zenith angle Θ 0 = arccos μ 0 , μ is cosine of the zenith angle of observation Θ , φ is azimuth angle of radiation propagation direction with respect to the solar vertical plane. The scattering angle γ (and x ( γ ) ) is used later only in the term for atmospheric haze radiation. According to [10], the choice of the scattering phase function is not important in determining the surface reflectance, because the equation includes the product of aerosol optical thickness and the phase function, each of which may be incorrectly defined, but the product and the value of the surface reflectance is correctly determined.
Accounting for the absorption bands of the major gaseous components of the atmosphere (water vapor, ozone and oxygen) is performed by filter method, that is, we use the general expression for spectral reflectance at TOA recorded excluding gas absorption components, multiplied by the transmissions of the three gas each components:
T g ( λ ) = T H 2 O λ T O 2 λ T O 3 λ .
To account for the transmission in the bands of absorption of these gases, atmospheric transmissions were calculated for standard absorbing mass values for each of these components with a spectral resolution of 2 nm, at the zenith Sun’s location, nadir viewing, twice the radiation passage (up- and downwelling) and average surface reflectance 0.2, as shown in Figure 1 in [11]. The standard conditions used to calculate the transmissions correspond to the following values of various parameters:
  • standard surface temperature 293 K;
  • standard surface pressure 101.3 kPa;
  • standard total integrated precipitable water: 4.20 g/cm 2 ;
  • standard integrated ozone amount 0.330 atm · cm.
For transmission with the absorbing mass, which differs from the standard for other Sun’s and observation angles, it is necessary to use the following expressions:
T H 2 O λ = T H 2 O λ 0 m 1 , T O 2 λ = T O 2 λ 0 m 2 , T O 3 λ = T O 3 λ 0 m 3 ,
where T H 2 O λ 0 , T O 2 λ 0 , T O 3 λ 0 are standard transmission functions of respective gaseous components for the conditions given above.
The following assumptions were used for the concentration profiles of atmospheric gases:
  • The fixed oxygen content is assumed, so the parameters m 2 depend only on the air mass.
  • In contrast, due to considerable variability, the transmission function of water vapor depends on the unknown concentration. The parameter m 1 is adjustable, including both the concentration of water vapor and the effective path of radiation (effective air mass).
  • Although variations in the ozone concentration are usually minor, the public datasets of the European Center for Medium-Range Weather Forecasts (ECMWF) were used to estimate the actual ozone concentration in the selected area.
Thus, in order to reduce unknown parameters in the solution of the inverse problem, the parameters m 2 , m 3 are assumed to be known (calculated) when solving the first stage of the inverse problem. In the next step, aimed at finding the surface reflectance of each pixel, parameters are re-estimated.
It should be noted that the use of a relatively simple model of the aerosol component of the atmosphere (see Equations (4)–(6)) also agrees with the conclusions of the paper [10] that, on the one hand, the spectra of outgoing radiation on TOA in the range 400–900 nm do not contain sufficient information to uniquely select aerosol type (four unknown parameters are entered for describing aerosol ( τ sca , a 0 , τ abs , a , β , g a ), but, on the other hand, retrieved spectra of surface reflectance are sufficiently stable with respect to the aerosol model of the atmosphere. This means that the spectrum of outgoing radiation depends on the aerosol model, but this dependence is not so strong as on the parameters mentioned above.

3. Approximate Analytical Equation for Spectral Radiance at TOA

For the at-sensor (satellite) spectral reflectance, we use the following expression:
R λ ( μ , μ 0 , φ ) = R atm λ ( μ , μ 0 , φ ) + E λ ( μ 0 , ρ e λ ) T λ dir ( μ ) ρ λ + ρ e λ T λ dif ( μ ) T g λ ,
where R atm λ ( μ , μ 0 , φ ) = B atm λ ( μ , μ 0 , φ ) / μ 0 S λ , B atm λ ( μ , μ 0 , φ ) is spectral radiance of atmospheric haze, S λ solar brightness constant, R atm λ ( μ , μ 0 , φ ) is the spectral atmospheric reflectance of haze (corresponds to zero surface reflectance), ρ e λ is the average surface reflectance around the current pixel with the reflectance ρ λ , E λ ( μ 0 , ρ e λ ) is the spectral illuminance of the Earth’s surface, normalized to that on TOA one, E 0 λ = π S λ μ 0 , depending on the average surface reflectance ρ e λ and zenith solar angle; T λ dir ( μ ) , T λ dif ( μ ) are direct and diffuse atmospheric transmission, respectively, from the surface to satellites sensor excluding gas transmission (which is included in the term T g λ ). Total transmission (without gaseous absorbing bands) is the sum of direct and diffuse transmissions:
T λ ( μ ) = T λ dir ( μ ) + T λ dif ( μ ) = exp ( τ λ / μ ) + T λ dif ( μ ) .
It should be noted that Equation (9) differs in form, but is completely equivalent to the generally written equation given in [15], which, in addition to the quantities described above, also includes the downwelling transmission (from TOA to the surface) T λ ( μ ) and spherical albedo of the atmosphere s ( λ ) . We obtain an expression for spectral radiance in this form by substituting in (9) the expression for normalized on the Earth’s surface illuminance E 0 λ = π S λ μ 0 , (see in [16,17]):
E λ ( μ 0 , ρ e λ ) = T ( λ ) 1 s ( λ ) ρ e λ .
The basis of analytical approximation for the outgoing radiation (9) consists of expressions for three functions: spectral illuminance E λ ( μ 0 , ρ e λ ) , total transmission of the atmosphere T λ ( μ ) and the spectral reflectance of atmospheric haze R atm λ ( μ , μ 0 , φ ) . The approximation error for each of these functions and Equation (9) as a whole does not exceed 3–5%.
We propose the following approximation for the illuminance (9) with some modification of the Eddington approximation [18,19]:
E λ ( μ 0 , ρ λ ) = ω λ E Ed λ ( μ 0 , ρ λ ) + ( 1 ω λ ) π S λ μ 0 exp ( τ λ / μ 0 )
E Ed λ ( μ 0 , ρ λ ) = 4 4 + 3 ( 1 g λ ) ( 1 ρ λ ) τ λ 1 2 + 3 4 μ 0 + 1 2 3 4 μ 0 exp τ λ μ 0 .
The average cosine of the scattering angle in Equation (13) must be calculated in accordance with the definition of the total phase function of the elementary volume (5). Since an average cosine of the scattering angle for molecular one is zero, it follows:
g λ = g a τ sca , a / τ sca , m + τ sca , a .
The error of approximation Equations (12)–(13) to calculate solar radiation fluxes (illuminance) is 1–2% [18], which was also confirmed by the comparisons we carried out with precise calculations using COART program and Equations (10)–(11) in the spectral range 0.4–1.1 μ m [11].
For a total transmittance function (excluding absorption gaseous bands) T λ ( μ ) , an analytic approximation from [20] was used. This approximation is obtained for the values of parameters 0 < g λ < 0.9 ; 0.2 < μ < 1.0 ; 0 < τ λ < 2 with a maximum error of about 8% for the g λ [ 0 0.9 ] , τ λ [ 0 2 ] , μ [ 0.2 1.0 ] , and an error of less than 4% for τ λ 1.6 , g λ 0.8 and μ [ 0.2 1.0 ] .
The advantage of this approximation of the function for total atmospheric upwelling transmittance T λ ( μ ) is a sufficiently high accuracy. Moreover, it does not contain new unknown parameters and depends only on total vertical optical thickness of the atmosphere τ λ , the cosine of the observation angle μ and the average cosine of the scattering angle g λ .
For the calculation of the spectral reflectance of atmospheric haze R atm λ ( μ , μ 0 , φ ) in Equation (9) the following approximation equation was proposed:
R atm λ ( μ , μ 0 , φ ) = R atm λ S ( μ , μ 0 , φ ) 1 + q ( ω λ τ λ ) p , p = 1.25 R atm λ S ( μ , μ 0 , φ ) = ω λ 4 x ( γ ) μ + μ 0 1 exp τ λ 1 μ 0 + 1 μ ,
in which the contribution of multiple scattering to the radiation of atmospheric haze is taken into account in the form of a quasi-linear amendment to the atmospheric haze reflectance in the single scattering approximation, R atm λ S . Approximations similar to Equation (15) have been proposed and investigated in [21]. We are using this approximation only for atmospheric haze contributions (for the case of ρ λ = 0 ), where the constant q is an unknown adjustable model parameter. We include it to the list of unknown values of the inverse problem. Accounting for the contribution of atmospheric haze in a single scattering approach, as the calculations showed, is not satisfactory, even for very clear atmosphere, while the representation (15) provides a sufficiently high accuracy [11].
The set of Equations (9) and (12)–(15) together with the formulas for the total transmission [20], gives an analytic representation of the spectrum of outgoing radiation, depending on the following seven parameters of the atmospheric optical model and surface:
τ abs , a , τ sca , a 0 , β , g a , q , ρ λ , m 1 .
Note that we put here the index λ for reflectance ρ λ , emphasizing its dependence on the wavelength. Minimizing the variation of the objective function (9) with respect to the measured value, only spectral independent parameters can be found. Therefore, in the first stage of the algorithm for finding the parameters of the atmosphere, one-parameter function ρ λ is introduced.
The test results of analytical approximation of spectral reflectance (9), (12)–(15) are given in [11]. The maximum error of the analytic representation in the spectral range of 0.4–0.65 microns is less than 4% and it can reach 10% in the range 0.35–1.1 microns. However, errors greater than 4% are only observed in the absorption bands of gases, in particular water vapor, suggesting the need for a more accurate account of profiles of gaseous components of the atmosphere.
We found that even better accuracy can be achieved if instead of one parameter m 1 , we use two parameters for the water vapor content ( m 11 and m 12 ): one in the term for the reflection by the haze and other one for the surface. Finally, Equation (9) could be rewritten (to shorten the reference, hereinafter everywhere we have omitted angular variables):
R λ = R atm λ ( T H 2 O λ ) m 11 + E λ ( ρ e λ ) ( T λ dir ρ λ + ρ e λ T λ dif ) ( T H 2 O λ ) m 12 T O 2 λ m 2 T O 3 λ m 3 .

4. Atmospheric Correction of Hyperspectral Imagery

The input data for the algorithm are: the reflection spectra of the hyperspectral image R λ ; the geometric parameters (the Sun’s zenith angle θ 0 , determined by the geographical coordinates and the acquiring time, the recording zenith angle θ (Earth-sensor direction), the azimuth angle ϕ ), the TOA Sun’s brightness function S λ [ m W / m 2 / n m / s r ] , the spectral transmitting functions of oxygen T O 2 λ , ozone T O 3 λ , and water vapor T H 2 O λ (dimensionless).
The algorithm of atmospheric correction see (Figure 1) consists of the following basic steps:
  • The selection of original pixel and area in the image under study for which the atmospheric parameters are determined.
    (a)
    We choose a “dark” pixel, if possible, on the hyperspectral image (low reflectance value, it is important to note that it is “dark” in the “blue-green” part of the spectrum, where there is the largest contribution of atmospheric haze). The user has the option of selecting a “dark pixel” either interactively or automatically based on a high correlation of a hyperspectral pixel with the library of various objects such as water, dark soil, asphalt, coniferous forest, etc.
    (b)
    In the case of absent suitable “dark pixel” in the image, a pixel with an approximately identifiable underlying surface (for example: vegetation, water, soil, sand, etc., or a mixture thereof) is selected, which is determined visually from RGB image.
  • Setup of the initial (zero) algorithm’s iteration for the reflectance of the selected pixel.
    (a)
    In case of a dark pixel surface (step 1a), we use
    ρ λ = ρ λ ( 0 ) = c ,
    in Equation (17). Since in this case the contribution of atmospheric haze to the detected signal exceeds the contribution of the reflected radiation from the surface (low reflectance value), then the spectral dependence of ρ λ can be neglected. Therefore, the reflectance can be assumed constant over the spectrum.
    (b)
    In the case of a homogeneous pixel (step 1b) (“pure”—vegetation, water, soil, sand, etc), we assume that the spectral surface reflectance function ρ λ ( 0 ) of the selected pixel can be presented in the form
    ρ λ = ρ λ ( 0 ) = c ρ λ l i b ,
    where c is the unknown weighted positive parameter, ρ λ l i b is the spectral surface reflectance library function of the identified surface type of the selected pixel.
    If the selected surface for the pixel is an inhomogeneous (consisting of a mixture of several types of surfaces) one, then the reflectance is assumed to be a linear combination of two dominant surface types with unknown parameter c (which here can vary in the range [0, 1]):
    ρ λ = ρ λ ( 0 ) = c ρ λ l i b 1 + 1 c ρ λ l i b 2 ,
    where ρ λ l i b i is the spectral surface reflectance library. By selecting a pixel over land, it is recommended to choose a linear combination of typical surface reflectance of vegetation and bare soil ([10]).
  • The choice of the neighborhood around the selected pixel (as an arbitrary polygon).
    The neighborhood should be (if possible) of the same type (with close reflectance values to the selected source pixel) as the selected pixel itself. The given zero-th order surface reflectance functions for the original pixel are considered as such for the whole selected neighborhood, i.e., are given by one of the formulas (18)–(20).
    The neighborhood area selected here is used to account for an adjacency effect in a non-traditional way: firstly, the initial iteration for atmospheric parameters is made with the average reflection signal via selected area around original (central) pixel, and then the new atmospheric parameters and reflectance of the original (central) pixel are found.
  • Finding the first iteration of the optical atmospheric parameters and average reflectance of the neighborhood area.
    We consider for this case ρ λ = ρ e λ = ρ λ ( 0 ) in (17), where in accordance with the choice of the original pixel, ρ λ ( 0 ) is determined by one of the formulas (18)–(20), and one can write (17) in the form
    R ¯ λ = R atm λ ( T H 2 O λ ) m 11 + ρ e λ ( 0 ) E λ ( ρ e λ ( 0 ) ) T λ dir + T λ dif ( T H 2 O λ ) m 12 T O 2 λ m 2 T O 3 λ m 3 .
    Here R ¯ λ is the average spectral reflectance of the hyperspectral image around the original pixel in the user-selected area (neighborhood, step 2). The next step is non-linear fit of the average spectral reflectance of selected area R ¯ λ by given analytical formula (right-hand side of Equation (21)) with the non-linear least squares Levenberg–Marquardt (LM) algorithm finding a set of unknown optical parameters of the atmosphere (see 16) and average surface reflectance: τ abs , a , τ sca , a 0 , β , g a , q , m 11 , m 12 , c . In (16), the water vapor exponent m 1 is replaced by two ones, m 11 and m 12 , in accordance with (17), and instead of ρ λ the weighted parameter c stands here for the mean reflectance of the chosen neighborhood (from (18) or (20)).
  • Additional smoothing of the spectral curve ρ λ of the current pixel in the water absorption spectral bands by re-fitting, where we vary just m 11 and m 12 , corresponding to the water vapor.
  • Refinement (next iteration).
    The atmospheric parameters τ abs , a , τ sca , a 0 , β , g a , q , m 11 , m 12 and reflectance of the original (central) pixel ρ λ ( 1 ) , which may differ from the neighborhood reflectance and which are determined by one of Equations (18)–(20) with a new unknown value of the parameter c = c 1 are found. In this case, reflectance of the neighborhood pixels remains the same as in the previous iteration (step 4).
    The (LM) non-linear least squares method is re-started using the following equation, which takes into account an adjacency effect (instead of (21)):
    R λ = R atm λ ( T H 2 O λ ) m 11 + E λ ( ρ e λ ( 0 ) ) T λ dir ρ λ ( 1 ) + ρ e λ ( 0 ) T λ dif ( T H 2 O λ ) m 12 T O 2 λ m 2 T O 3 λ m 3 .
    Therefore, the fitting of measured reflection spectrum of the original pixel R λ is performed by the objective function (22).
  • Application of smoothing filter to the current pixel by re-fitting using exponents m 2 and m 3 .
  • Calculation of the first approximation for the surface reflectance ρ λ .
    The surface reflectance is found for all other pixels of the hyperspectral image without the adjacency effect from the following quadratic equation (it is obtained by solving the Equation (17) for the variable ρ λ with putting there ρ e λ = ρ λ and taking into account the expressions (12,13) for E λ ρ λ ). Atmospheric parameters, which found in step 6 are substituted into Equation (17) (the atmosphere is assumed horizontally homogeneous, identical over all pixels). Then it follows:
    a λ ρ λ 2 b λ ρ λ + c λ = 0 ,
    where
    a λ = 3 τ λ ( 1 g a ) ( 1 ω λ ) exp τ λ / μ 0 ,
    b λ = 3 τ λ ( 1 g a ) R 1 λ + 4 ω λ 1 2 + 3 4 μ 0 + 1 2 3 4 μ 0 exp τ λ μ 0 + + 4 + 3 τ λ ( 1 g a ) 1 ω λ exp τ λ / μ 0 ,
    c λ = 4 + 3 τ λ ( 1 g a ) R 1 λ ,
    R 1 λ = R λ / T O 2 λ m 2 T O 3 λ m 3 R atm λ ( T H 2 O λ ) m 11 · T λ ( T H 2 O λ ) m 12 1
    Here R λ is the spectral reflectance of the processed pixel. Atmospheric parameters found above are substituted into Equations (24)–(27).
  • The adjacency effect area specification by pixel-wise.
    At this step we specify a fixed neighborhood of each pixel and calculate the average reflectance to account for the adjacency effect. In this case, the spectral surface reflectance of pixels from the previous stage is used. The contribution to the total spectral surface reflectance from neighboring pixels is accounted for with an exponentially decay weight function depending on distance from the central (current) pixel:
    ρ λ ¯ = 1 N i , j A exp a i 2 + j 2 d ρ λ , i , j ,
    where i , j are pixel numbers in the neighborhood, d is the half-width of the specified window in pixels, a , A are the distribution parameters, N is the number of neighboring pixels (d is selected using iterations).
  • The estimation of final reflectance value.
    The surface reflectance is found from the following formula, which follows from Equation (17) with ρ e λ = ρ λ ¯ :
    ρ λ = R λ T O 2 λ m 2 T O 3 λ m 3 R atm λ T H 2 O λ m 11 ρ ¯ λ E λ ( ρ ¯ λ ) T λ dif T H 2 O λ m 12 · E λ ( ρ ¯ λ ) T λ dir T H 2 O λ m 12 1
    where one needs to substitute the following values: the spectral reflectance of the processed (current) pixel R λ , the average spectral surface reflectance estimated for each current pixel using Equation (28) (calculated based on step 8) and atmospheric parameters found in step 6.

5. Validation

5.1. Validation Using Synthetic Spectra

The developed method Satellite Hypercube Atmosheric Rapid Correction (SHARC) was first tested using the spectra of outgoing radiation, obtained by numerical calculation of radiation transfer equation using the well-known open access programming code LibRadTran [22].
The validation using radiative transfer modeling, say, with LibRadTran, is an important part of the atmospheric correction, since it represents the most “clean” task version (without influencing some negative factors of real measurements), and therefore, it allows to verify the adequacy and accuracy of the atmospheric model and the analytical expressions used. The purpose of validation of the atmospheric correction method using the numerical calculations of the outgoing radiation spectrum is to check the method of atmospheric correction in the form that most corresponds to the objective function used for the solution of the inverse problem, that is the analytical spectral radiance function. What is this correspondence? A numerical calculation of the spectral radiance is carried out for a separate observation beam (as well as Equation (17)), i.e., the problem is a one-dimensional one. The surface and atmosphere are horizontally uniform, whereas the experimental value of the spectral radiance from the satellite sensor is averaged over the finite instantaneous field of view corresponding to a single pixel. In the numerical calculation of spectral radiance, the adjacency effect is not taken into account. It is also omitted in the analytical equations in this validation part (one should use ρ e λ = ρ λ in Equation (17)), but it is taken into account for real satellite data (there is a problem of the neighborhood influence selection, e.g., selection of the adjacency effect radius). Furthermore, usually instrumental errors are present while working with satellite imagery, whereas for numerically calculated spectral radiance, such errors are absent.
As shown by numerical calculations, the accuracy of the emission spectrum at the top of atmosphere by analytical equations is very high, which is confirmed by Figure 2. In Figure 2 we show two curves of the brightness (instead of the reflection coefficient) at TOA direct problem: obtained by numerical calculation of LibRadTran with the underlying surface "grass" and fitting curve of non-linear least squares LM algorithm with the analytical Equation (21).
The reconstruction accuracy of the spectral reflectance function of the surface (inverse task) is shown in Figure 3a–d, where we present the surface reflectance spectra for four different types of underlying surfaces (water, soil, grass, snow) obtained by SHARC technique using numerically calculated (LibRadTran) TOA brightness spectra, in comparison with the true reflectance data given input for LibRadTran. Figure 3 reflects the typical retrieved reflectance spectra errors obtained for different types of underlying surfaces.

5.2. Validation Using Hyperion Measurements

In addition, atmospheric correction of hyperspectral satellite images obtained using the Hyperion sensor on the Earth Observer-1 (EO-1) was performed by us. Three different images with different types of surface were selected, including various types of underlying surfaces such as water, sand, green vegetation, rocks, etc., see Figure 4.
Results of the SHARC method were compared with atmospheric correction data of these images in the well-known FLAASH module, implemented in ENVI software package. Typical comparison results for images of Hyperion sensor are shown in Figure 5, Figure 6 and Figure 7. In the left hand side of the figures, the pixel brightness spectra in absolute units (after radiometric correction) registered by the sensor are shown. The right hand side contains the surface reflectance spectra retrieved by SHARC method and FLAASH module.
It should be noted that the retrieved profiles of SHARC and FLAASH reflectances are in satisfactory agreement. In this case, the SHARC curve is smoother in comparison with FLAASH and shows somewhat better (smoother) behavior in the region of the absorption bands of water vapor around 940 nm. We see the “smoothness” of the spectral reflectance reconstructed by the SHARC method (both in the calculated spectra of LibRadTran and in Hyperion satellite spectrum images) in the absorption of water vapor bands. In addition, one observes that FLAASH often gives negative reflectance values at some spectral ranges (Figure 6 and Figure 7).
In some cases, the Hyperion sensor gives insufficient quality, as shown, in Figure 7 for the underlying water surface, where the original spectrum contains “emission bands” around 900 nm, which are most probably artifact.

5.3. Validation Using Ground Data

Moreover, to validate our data obtained by the SHARC method, a comparative analysis with ground-based measurements of Radiometric Calibration Network (RadCalNet) was carried out. RadCalNet is an initiative of the Working Group on Calibration and Validation of the Committee on Earth Observation Satellites (available online: https://www.radcalnet.org/). This portal provides an open access bottom-of-atmosphere (BOA) and TOA reflectance data measured derived over the network of sites. The data continuously updated, have a 10 nm spectral sampling interval in the spectral range from 380 to 2500 nm.
The Hyperion image that covers the La Crau site with location in the south of France (43.55N 4.85W) was used. This site area has a thin, pebbly soil with sparse vegetation cover. Figure 8 shows the Hyperion image footprint (red rectangle) with location near La Crau (yellow mark). The Hyperion image was acquired on 29 October 2016, 7:48:20 UTC (sun elevation angle: 14.6 deg). RadCalNet data were measured on 29 October 2016, 9:00 UTC (sun elevation angle: 24.11 deg). The results of in-situ measured surface reflectance, that corresponds to the pixel number (1719, 476) of Hyperion image mentioned above, were compared with surface reflectance retrieved by SHARC method for this pixel, as shown in Figure 9. A good correspondence of both spectra is found (except in the vicinity of 450 nm, where there is a slight difference in spectra).

6. Conclusions

A new algorithm for atmospheric correction of hyperspectral images, SHARC, is presented. This algorithm enables the determination of reflectance of the underlying surface “on the fly” without using a priori information about the atmosphere and the surface. In addition, we do not use LUTs or time-consuming numerical solutions of the radiative transfer equation. Instead, the algorithm uses a parameterized analytical solution of the radiative transfer equation for TOA outgoing radiation for a wide range of atmospheric parameters. The type of atmospheric aerosol is not predetermined depending on the season and geographical location, but it is characterized by four unknown parameters accounting for almost all possible types of aerosols.
The advantage of the proposed method is a high accuracy combined with a very short time of performing the procedure of atmospheric correction, even for very large images. It takes, for example, for the hypercube with dimension 1546 over 592 pixels and 68 spectral channels, about 15 s, on a modern personal computer of standard configuration, with a good accuracy of the spectral reflectance retrieval. The limitation of this method is the above-mentioned range of atmospheric parameters, for which approximations of the functions have been obtained: g λ 0 0.9 , τ λ 0 2 , μ 0.2 1.0 , which, however, covers most realizations of the cloudless atmosphere. Another limitation is the assumption of a horizontal homogeneity of the atmosphere within the image, since the parameters of the atmosphere are assumed to be the same for all pixels. If this assumption is not fulfilled (for example, for large images), the proposed correction method can be applied successively to image subsets (for example, using a sliding rectangular window) and within which the independent parameters of the atmosphere can be obtained. We plan to further develop the method, taking into account the possible horizontal inhomogeneity of the atmosphere, as well as the non-flat terrain.

Author Contributions

L.V.K. contributed to the conception of the paper and the development of the method. A.A.K. contributed to the to the development of some functions approximation, conception and design of paper and English translation. D.A.I. and A.O.M. implemented the method in the SHARC algorithm. They were responsible for implementation of this method to software code. L.V.K., A.O.M., V.A.S. and D.A.I. contributed to the analysis, processed the data and discussion. V.A.S. provided the ground-based validation of surface reflectances and LibRadTran data. A.O.M. and V.A.S. prepared figures. D.A.I. contributed to reviewing the manuscript before submission, English spelling and formatted it to the MDPI template. L.V.K., D.A.I. and A.A.K. finalized the manuscript. All authors contributed to the writing of the manuscript and improving the review.

Funding

This research received no external funding.

Conflicts of Interest

The authors declare no conflicts of interest

References

  1. Eismann, M. Hyperspectral Remote Sensing; SPIE Press Monograph; SPIE: Washington, WA, USA, 2012; Volume PM210, p. 725. ISBN 9780819487872. [Google Scholar]
  2. Miller, C.J. Performance assessment of ACORN atmospheric correction algorithm. In Algorithms and Technologies for Multispectral, Hyperspectral, and Ultraspectral Imagery VIII; SPIE: Orlando, FL, USA, 2002; Volume 4725, pp. 438–449. [Google Scholar]
  3. Richter, R.; Schlapfer, D. Geo-atmospheric processing of airborne imaging spectrometry data. Part 2: Atmospheric/topographic correction. Int. J. Remote Sens. 2002, 23, 2631–2649. [Google Scholar] [CrossRef]
  4. Gao, B.-C.; Heidebrecht, K.B.; Goetz, A.F.H. Derivation of scaled surface reflectances from AVIRIS data. Remote Sens. Environ. 1993, 44, 165–178. [Google Scholar] [CrossRef]
  5. Adler-Golden, S.M.; Matthew, M.W.; Bernstein, L.S.; Levine, R.Y.; Berk, A.; Richtsmeier, S.C.; Acharya, P.K.; Anderson, G.P.; Feldeb, G.; Gardner, J.; et al. Atmospheric correction for short-wave spectral imagery based on MODTRAN4. In Imaging Spectrometry V; SPIE: Denver, CO, USA, 1999; Volume 3753, pp. 61–69. [Google Scholar]
  6. Goetz, A.F.H.; Kindel, B.C.; Ferri, M.; Qu, Z. HATCH: Results from simulated radiance, AVIRIS and Hyperion. JIEEE Trans. Geosci. Remote Sens. 2003, 41, 1215–1222. [Google Scholar] [CrossRef]
  7. Qu, Z.; Kindel, B.C.; Goetz, A.F.H. The high accuracy atmospheric correction for hyperspectral data (HATCH) model. IEEE Trans. Geosci. Remote Sens. 2003, 41, 1223–1231. [Google Scholar]
  8. Montes, M.J.; Gao, B.-C.; Davis, C.O. Tafkaa atmospheric correction of hyperspectral data. In Imaging Spectrometry IX; SPIE: San Diego, CA, USA, 2003; Volume 5159, pp. 188–197. [Google Scholar]
  9. Leprieur, C.; Carrere, V.; Gu, X.F. Atmospheric corrections and ground reflectance recovery for Airborne Visible/Infrared Imaging Spectrometer (AVIRIS) data: MAC Europe’91. Photogramm. Eng. Remote Sens. 1995, 61, 1233–1238. [Google Scholar]
  10. Katsev, I.L.; Prikhach, A.S.; Zege, E.P.; Grudo, J.O.; Kokhanovsky, A.A. Speeding up the aerosol optical thickness retrieval using analytical solutions of radiative transfer theory. Atmos. Meas. Tech. 2010, 3, 1403–1422. [Google Scholar] [CrossRef] [Green Version]
  11. Katkovsky, L.V. The parameterization of the outgoing radiation for rapid atmospheric correction of hyperspectral images. Optika Atmosfery i Okeana 2016, 29, 778–784. [Google Scholar]
  12. Ginzburg, A.S.; Romanov, S.V.; Fomin, B.A. The radiation-convective model use to estimate the temperature of greenhouse gas potential. Izvestiya RAN Fizika Atmosfery i Okeana 2008, 44, 324–331. [Google Scholar]
  13. Belyaev, B.I.; Belyaev, M.Y.; Desinov, L.V.; Katkovsky, L.V.; Sarmin, E.E. Spectral and Images Processing from Photospectral System in Space Experiment “HURRICANE” on the ISS. Izvestiya, Atmospheric and Oceanic Physics 2014, 6, 54–65. (In Russian) [Google Scholar]
  14. Bucholtz, A. Rayleigh-scattering calculations for the terrestrial atmosphere. Appl. Opt. 1995, 34, 2765–2773. [Google Scholar] [CrossRef] [PubMed]
  15. Bassani, C.; Cavalli, R.M.; Antonelli, P. Influence of aerosol and surface reflectance variability on hyperspectral observed radiance. Atmos. Meas. Tech. 2012, 5, 1193–1203. [Google Scholar] [CrossRef] [Green Version]
  16. Middleton, W.E.K. Vision Through the Atmosphere; University of Toronto Press: Toronto, ON, Canada, 1952. [Google Scholar]
  17. Schlapfer, D.; Borel, C.C.; Keller, J.; Itten, K.I. Atmospheric Precorrected Differential Absorption Technique to Retrieve Columnar Water Vapor. Remote Sens. Environ. 1998, 65, 353–366. [Google Scholar] [CrossRef]
  18. Vasil’ev, A.V.; Kuznetsov, A.D.; Mel’nikova, I.N. Remote Sensing of the Environment from Space: Practice; Baltic State Technical University: St. Petersburg, Russia, 2008; p. 33. [Google Scholar]
  19. Minin, I.N. Approximate equations for short-wave radiation absorption calculations in cloudless atmosphere. Izv. Akad. Nauk SSSR Fiz. Atm. Okeana 1984, 20, 999–1001. (In Russian) [Google Scholar]
  20. Kokhanovsky, A.A.; Mayer, B.; Rozanov, V.V. A parameterization of the diffuse transmittance and reflectance for aerosol remote sensing problems. Atmos. Res. 2005, 73, 37–43. [Google Scholar] [CrossRef] [Green Version]
  21. Vasil’ev, A.V.; Kuznetsov, A.D.; Mel’nikova, I.N. Approximation of multiply scattered solar radiation in the assumption of single scattering. In Proceedings of the International Symposium “Atmospheric Radiation and Dynamics” (ISARD 2015), Saint-Petersburg, Russia, 23–26 June 2015; p. 131. [Google Scholar]
  22. Emde, C.; Buras-Schnell, R.; Kylling, A.; Mayer, B.; Gasteiger, J.; Hamann, U.; Kylling, J.; Richter, B.; Pause, C.; Dowling, T.; et al. The libradtran software package for radiative transfer calculations (version 2.0.1). Geosci. Model Dev. 2016, 9, 1647–1672. [Google Scholar] [CrossRef] [Green Version]
Figure 1. Pipeline of atmospheric correction of hyperspectral image data.
Figure 1. Pipeline of atmospheric correction of hyperspectral image data.
Remotesensing 10 01698 g001
Figure 2. Spectral radiance on the top of atmosphere (TOA) according to LibRadTran calculations and fitting an analytic function (21).
Figure 2. Spectral radiance on the top of atmosphere (TOA) according to LibRadTran calculations and fitting an analytic function (21).
Remotesensing 10 01698 g002
Figure 3. Surface reflectances for four different surface types (water (a), bare soil (b), grass (c) and snow (d) retrieved by SHARC technique using the numerically calculated (LibRadTran) outgoing TOA radiance spectra, compared to the surface reflectance spectra, initially given in direct problem.
Figure 3. Surface reflectances for four different surface types (water (a), bare soil (b), grass (c) and snow (d) retrieved by SHARC technique using the numerically calculated (LibRadTran) outgoing TOA radiance spectra, compared to the surface reflectance spectra, initially given in direct problem.
Remotesensing 10 01698 g003aRemotesensing 10 01698 g003b
Figure 4. RGB images of the Hyperion based on retrieved spectral surface reflectance functions using the SHARC method.
Figure 4. RGB images of the Hyperion based on retrieved spectral surface reflectance functions using the SHARC method.
Remotesensing 10 01698 g004
Figure 5. Initial radiance spectra of Hyperion (left column) and reflectance spectra for the same pixel retrieved by the SHARC atmospheric correction technique (right column) and FLAASH module (ENVI) (right column) for three types of underlying surfaces (coastal water (a), green vegetation (b), sand (c)), labeled with pins in Figure 4a.
Figure 5. Initial radiance spectra of Hyperion (left column) and reflectance spectra for the same pixel retrieved by the SHARC atmospheric correction technique (right column) and FLAASH module (ENVI) (right column) for three types of underlying surfaces (coastal water (a), green vegetation (b), sand (c)), labeled with pins in Figure 4a.
Remotesensing 10 01698 g005
Figure 6. Initial radiance spectra of Hyperion (left column) and reflectance spectra for the same pixel retrieved by the SHARC atmospheric correction technique (right column) and FLAASH module (ENVI) (right column) for three types of underlying surfaces (green vegetation (a), rocks (b), coastal water (c)), labeled with pins in Figure 4b.
Figure 6. Initial radiance spectra of Hyperion (left column) and reflectance spectra for the same pixel retrieved by the SHARC atmospheric correction technique (right column) and FLAASH module (ENVI) (right column) for three types of underlying surfaces (green vegetation (a), rocks (b), coastal water (c)), labeled with pins in Figure 4b.
Remotesensing 10 01698 g006
Figure 7. Initial radiance spectra of Hyperion (left column) and reflectance spectra for the same pixel retrieved by the SHARC atmospheric correction technique (right column) and FLAASH module (ENVI) (right column) for three types of underlying surfaces (soil (a), green vegetation (b), water (c)), labeled with pins in Figure 4c.
Figure 7. Initial radiance spectra of Hyperion (left column) and reflectance spectra for the same pixel retrieved by the SHARC atmospheric correction technique (right column) and FLAASH module (ENVI) (right column) for three types of underlying surfaces (soil (a), green vegetation (b), water (c)), labeled with pins in Figure 4c.
Remotesensing 10 01698 g007
Figure 8. La Crau site with Hyperion image footprint.
Figure 8. La Crau site with Hyperion image footprint.
Remotesensing 10 01698 g008
Figure 9. Hyperion radiance spectrum (left side) and measured reflectance of RadCalNet network in comparison with calculated reflectance by SHARC method (right side) for the pixel shown in Figure 8.
Figure 9. Hyperion radiance spectrum (left side) and measured reflectance of RadCalNet network in comparison with calculated reflectance by SHARC method (right side) for the pixel shown in Figure 8.
Remotesensing 10 01698 g009
Table 1. Values of coefficients for calculating the Rayleigh optical thickness of the total atmosphere from Equation (3).
Table 1. Values of coefficients for calculating the Rayleigh optical thickness of the total atmosphere from Equation (3).
Atmospheric Model λ = 0.2–0.5 μ m λ > 0.5 μ mPressure P s , Mbar at the Earth’s SurfaceTemperature T s , K at the Earth’s Surface
BFor all models3.552123.99668
C1.355790.00110298
D0.115630.0271393
F P s , T s Tropical0.0065258410.0086800891013300
Midlatitude Summer0.0065155470.0086659971013294
Midlatitude Winter0.0065318960.0086884021018272.2
Subarctic Summer0.0064775390.0086161751010287
Subarctic Winter0.0064958230.0086417421013257.1
1962 US Standard0.0064995950.0086452611013288.1

Share and Cite

MDPI and ACS Style

Katkovsky, L.V.; Martinov, A.O.; Siliuk, V.A.; Ivanov, D.A.; Kokhanovsky, A.A. Fast Atmospheric Correction Method for Hyperspectral Data. Remote Sens. 2018, 10, 1698. https://doi.org/10.3390/rs10111698

AMA Style

Katkovsky LV, Martinov AO, Siliuk VA, Ivanov DA, Kokhanovsky AA. Fast Atmospheric Correction Method for Hyperspectral Data. Remote Sensing. 2018; 10(11):1698. https://doi.org/10.3390/rs10111698

Chicago/Turabian Style

Katkovsky, Leonid V., Anton O. Martinov, Volha A. Siliuk, Dimitry A. Ivanov, and Alexander A. Kokhanovsky. 2018. "Fast Atmospheric Correction Method for Hyperspectral Data" Remote Sensing 10, no. 11: 1698. https://doi.org/10.3390/rs10111698

APA Style

Katkovsky, L. V., Martinov, A. O., Siliuk, V. A., Ivanov, D. A., & Kokhanovsky, A. A. (2018). Fast Atmospheric Correction Method for Hyperspectral Data. Remote Sensing, 10(11), 1698. https://doi.org/10.3390/rs10111698

Note that from the first issue of 2016, this journal uses article numbers instead of page numbers. See further details here.

Article Metrics

Back to TopTop