A Robust Atmospheric Correction Procedure for Determination of Spectral Reflectance of Terrestrial Surfaces from Satellite Spectral Measurements

Abstract

In this work, we propose simple and robust technique for the retrieval of underlying surface spectral reflectance using spaceborne observations. It can be used to process both multispectral moderate resolution satellite data and also multi-zone high spatial resolution data. The technique can work automatically for different types of land surfaces without using huge databases accumulated in advance. The new cloud discrimination and retrieval of the water vapor content in atmosphere procedures are presented. The key point of the proposed atmospheric correction technique is the suggested single-wavelength method for determining the atmospheric aerosol optical thickness without reference to a specific type of underlying surface spectrum. The retrievals of spectral reflectance for various land surfaces with developed technique, performed using computer simulation and experimental data, have demonstrated a high retrieval accuracy.

1. Introduction

The goal of atmospheric correction is to eliminate the influence of the atmosphere on the brightness and spectral contrasts in satellite images and to retrieve the spectral characteristics of the underlying surface, in particular, the reflectance of the Earth’s surface.

The main factors that distort the satellite signals from the Earth’s surface are aerosols and molecular scattering, as well as absorption by gases and particulate matter suspended in the atmosphere. The problem of gas absorption can be mainly solved by choosing spectral channels outside the absorption bands or by taking into account this absorption for known atmospheric gases. Molecular scattering in the atmosphere is quite stable and can be accounted relatively easily using various latitudinal and regional models.

The most variable atmospheric components are aerosol and water vapor. Therefore, to eliminate the distorting effect of the atmosphere, it is most important to know such parameters as the aerosol optical thickness (AOT) of the atmosphere and the amount of water vapor in the atmospheric column. Another important task solved by atmospheric correction is the formation of the so-called cloud mask, that is, the detection and discrimination of clouds that substantially distort the received information about the atmosphere.

The complexity of atmospheric correction of multispectral satellite data is due to the fact that the spectral characteristics of the Earth’s surface and atmosphere should be retrieved simultaneously from the spectral radiance (SR) recorded by a space instrument. Therefore, it is no coincidence that the most advanced algorithms, such as Bayesian aerosol retrieval (BAR) [1] and Deep Blue (DB) [2,3] for processing MODIS data, as well as DG AComp for processing high-resolution satellite WorldView-2 and WorldView-3 data [4] use huge previously accumulated databases. For example, in the BAR algorithm, the monthly-expected reflectance values and their variations are calculated from the database of daily MODIS measurements of the Earth’s reflectance for the previous year. This database is recalculated monthly. When determining AOT of the atmosphere as the expected AOT value for each pixel, the values from the climatological database MAC-V2 [5], which contains monthly AOT values on a 1° × 1° grid, are used. In the DB algorithm [3], designed to determine the atmospheric AOT in areas with a high reflectance of a surface (desert, semi-desert), it is proposed first to estimate the surface reflectance for each pixel at wavelengths of 412, 470, and 650 nm, using, in particular, the previously calculated database created from MODIS and SeaWiFS measurements. Such a global database of reflectance of bright surfaces in the visible spectrum has a resolution of 0.1° × 0.1°.

The algorithms ART [6] and BAER [7] use another approach. Here, the reflectance spectrum of the underlying surface is sought in the form of a linear combination of two basic spectra. In particular, for dark surfaces, the spectra of green vegetation and soil are accepted as basic spectra. For different types of surfaces, it is necessary to choose and justify a different set of basic spectra. Naturally, the applicability of this method is limited to a certain class of the Earth’s surfaces. Moreover, experience shows that such a choice of basic spectra can be difficult and not always successful.

In this work, we propose and verify a new rather simple and fast method of atmospheric correction Robust Atmospheric Correction Enhancement (RACE). The technique can be used to process both multispectral (MS) moderate resolution satellite data and also multi-zone (MZ) high spatial resolution data. It can work automatically for different types of terrestrial surfaces without using huge databases accumulated in advance. Below, without dwelling on the details of the algorithm, its key, mostly original parts, as well as the results of its testing are described.

This article is divided into three main sections. The next section describes the proposed RACE procedure for fast and reliable atmospheric correction of satellite images. Section 3 presents the equations used to retrieve the reflectance of the underlying surfaces in MS and MZ channels. Section 4 includes the results of the verification and checking the correctness, stability, speed, and accuracy of the proposed RACE algorithm to retrieve the parameters of the atmosphere and the reflectance of the underlying surfaces. Testing was carried out using both computer simulation and available experimental data. Section 5 briefly summarizes the presented results.

2. Atmospheric Correction Procedure

This section presents the main ideas and procedures used in the developed atmospheric correction technique. It includes brief descriptions of used atmospheric model and auxiliary procedures, namely, techniques of cloud and snow pixels detection, and of retrieving the content of water vapor in the atmosphere. The final paragraph of this section presents the idea and rationale of the proposed speedy and robust atmospheric correction procedure.

2.1. Atmospheric Model

The basic atmospheric model is described in detail in [6,8]. The stratified atmosphere is conditionally divided into two parts:

• (1). Part “1” is a layer of the troposphere up to a certain height H (about 2–3 km). The aerosol in this layer can vary in time and space. Therefore, it is allowed that the AOT of this layer may vary for different parts of the image.

• (2). Part “2” (layer above height H) includes the stratosphere, as well as the upper and middle troposphere. It is characterized by vertical stratification of aerosol and gas concentrations, pressure and temperature. It consists of a large number of sublayers with optical characteristics averaged over these sublayers.

The profiles of the vertical distribution of pressure, temperature, and the main absorbing gases (ozone, oxygen, water vapor, carbon dioxide, and nitrogen trioxide) can be specified in accordance with the latitudinal seasonal model of the molecular-gas atmosphere. The parameters of the molecular-gas model are adjusted in accordance with the current value of atmospheric pressure at the Earth’s surface level and its height above sea level.

Our experience shows that a simpler model can be used in practice. In this model the atmosphere aerosol stratification is presented in three layers: 1) the lower troposphere to an altitude of 2 km, 2) the middle troposphere to an altitude of about 6 km, 3) the upper troposphere and stratosphere to an altitude of about 50 km. In each of these layers, a single type of aerosol is adopted. In the stratospheric layer, this is usually sulfuric acid (H₂SO₄), which is the result of volcanic emissions. The aerosol height distributions in these layers are assumed to be uniform. At altitudes above 50 km the presence of aerosols can be practically neglected. It is sufficient to take into account only scattering and absorption by the molecular atmosphere.

Note that the described model of the atmosphere was used in all computer simulations described below. The optical thicknesses of the aerosol layers at a wavelength of 550 nm were: 0.019 for the upper stratospheric layer (H₂SO₄ aerosol) and 0.02 for the middle layer (model “Continental”). In addition, it is assumed that the lower layer also contains aerosol “Continental”, but its AOT may vary.

With a large number of pixels in the image frame and taking into account the vertical structure of the optical parameters, retrieval of the AOT in the lower atmosphere layer requires a large amount of computation. However, as numerous computer simulations proved, the radiation characteristics of the lower layer can be calculated for a homogeneous layer with averaged scattering and absorption characteristics and without taking into account the polarization of radiation. Indeed, calculations show that this approximation leads to a relative error of less than 0.2% when calculating the radiance coefficients at the top-of-atmosphere (TOA) for wavelengths in the visible spectrum. The calculation of the radiation characteristics of atmospheric layers has been carried out using the effective RAY code [9].

2.2. Cloud and Snow Pixels Detection

A cloud mask includes the detection of clouds and snow pixels, the separation of one from the other, as they can be indistinguishable in the visible region. A special attention is given to the detection of high cirrus clouds. Methods to discriminate clouds and snow pixels used in various algorithms varies markedly, but, of course, all of them are based on a high reflectance of clouds and snow. As was shown in [10], when clouds completely cover a pixel, $R_{TOA}\left(412\mathrm{nm}\right)\ge 0.2$. However, to discriminate clouds over bright desert-like surfaces, even the criterion $R_{TOA}\left(412\mathrm{nm}\right)\ge 0.4$ ($R_{TOA}\left(\lambda \right)$ is the reflection function (RF) at the TOA) is unsuitable at some observation geometry because it cuts off not only cloudy, but also cloudless pixels (Figure 1) due to the significant contribution of molecular scattering. Therefore, the detection of optically thick cumulus clouds and snow pixels is proposed to be carried out according to the following criterion:

(1)$R_{TOA}\left(412\mathrm{nm}\right)-R_{mol}\left(412\mathrm{nm}\right)>\eta ,$

i.e., after subtracting the contribution R_mol of the molecular-gas atmosphere at the TOA. The proposed threshold value is $\eta =0.25$.

Optically thick clouds are separated from snow pixels with the Normalized Difference Snow Index (NDSI) [11]:

(2)$NDSI=\frac{R_{TOA}\left(559\mathrm{nm}\right)-R_{TOA}\left(1.605\mathrm{nm}\right)}{R_{TOA}\left(0.559\mathrm{nm}\right)+R_{TOA}\left(1.605\mathrm{nm}\right)}.$

It follows for the snow-covered surfaces: $NDSI>0.8$ [11].

The detection and discrimination of high cirrus clouds is performed using the condition $R_{TOA}>0.02$ in the spectral range 1360–1380 nm, where there is a powerful absorption band of water vapor [11]. Therefore, sun light, reflected from the Earth’s surface in this spectral range, practically does not reach the satellite instrument.

2.3. Determination of Water Vapor Content in the Atmosphere

Water vapor is one of the most important components of the atmosphere, affecting the signals recorded by the satellite sensor. The content of water vapor in the atmosphere, in contrast to the content of most atmospheric gases, varies significantly in time and space. Therefore, determining the total content of water vapor in the atmosphere column is important for atmospheric correction of satellite data.

The most common method determining the absorbing gases content in the atmosphere is the differential method, when signals in the absorption band and outside the absorption band at a small distance from it are compared. Wherein it is believed that the remaining parameters of the atmosphere and surface within the absorption band and outside it are the same. Naturally, failure to fulfill this condition can lead to the corresponding error. Moreover, this approach is applicable only in the approximation of the filter model, when we can assume that all absorption takes place at the upper boundary of the atmosphere.

The authors of [12,13,14], consider the Weighting Function Modified Differential Absorption Spectroscopy (WFM-DOAS) technique, where the water vapor content is retrieved only by using the absorption band. However, wherein the proposed models of the atmosphere, underlying surface, and survey geometry are used, the water vapor content is determined in an iterative process based on the ratio:

(3)$\ln J_{\mathrm{meas}}=\ln J_{\mathrm{mod}}\left(W_n\right)+\frac{d\ln J_{\mathrm{mod}}\left(W\right)}{dW}\Delta W_n,,$

(4)$\Delta W_n=\left(W_{n+1}-W_n\right),$

where W is the total content of water vapor in the atmosphere, $J_{\mathrm{meas}}$ is the spectral radiance (SR) measured at the TOA, $J_{\mathrm{mod}}$ is the SR, computed for given atmosphere model. The derivative in Equation (3) is calculated at the point $W_n$. Note that this method is not limited to the applicability of the filter model. However, since its result is the value of W at $J_{\mathrm{mod}}=J_{\mathrm{meas}}$, it is sensitive to inaccuracies of used atmosphere model, value of surface reflectance and radiometric calibration.

In our opinion, it is preferable to use a slightly different scheme to determine the water vapor content, namely, the iterative process based on the logarithm of the ratio of reflectance values:

(5)$F=\ln \frac{J_1}{J_2}.$

Namely, we use iterations as shown below:

(6)$F_{\mathrm{meas}}=F_{\mathrm{mod}}\left(W_n\right)+\frac{d{F}_{\mathrm{mod}}\left(W\right)}{dW}\Delta W_n,$

(7)$\Delta W_n=\left(W_{n+1}-W_n\right),$

where $J_1$ and $J_2$ are SR at the TOA within the absorption band and outside it correspondingly. In this case, the iterative process for retrieving the water vapor content in the atmosphere is arranged according to the scheme:

(8)$W_{n+1}=W_n+\frac{F_{\mathrm{meas}}-F_{\mathrm{mod}}\left(W_n\right)}{d{F}_{\mathrm{mod}}/dW}.$

The function $F_{\mathrm{meas}}$ is calculated from the values of measured SR, and value $F_{\mathrm{mod}}$‘ is determined for the used atmosphere model. In this scheme, as in the WFM-DOAS method, the conditions necessary for the application of filter models and functions, calculated taking into account the atmospheric models and the reflectance of the underlying surface, are not required.

Let us emphasize that since in this scheme the function $F$ is defined as the logarithm of the ratio of signals in close spectral ranges, it should be more stable to all errors listed above than in the WFM-DOAS method.

To determine the content of water vapor W in the atmosphere, one can use absorption bands of about 910 nm, 940 nm, or other channels in the range of 900–970 nm with a reference signal in the range of 860–880 nm. The surface spherical reflectance $r_{surf}$ at the wavelength of the reference signal, for example $r_{surf}\left(865\mathrm{nm}\right)$, with a known model of the atmosphere, is easily determined based on the atmospheric correction Equation (9) (see below). To estimate the value of $r_{surf}$ within the absorption band, one can use linear extrapolation over wavelengths of 778.5 nm and 865 nm to the region of the absorption band of water vapor.

Results of computer simulation showed the stability and rapid convergence of scheme (8). Moreover, the simulations demonstrate possibility to retrieve the water vapor content even in cases when parameters of the atmosphere and surface within the absorption band and in the reference channel do not coincide on condition that they are known with good accuracy.

2.4. Aerosol Optical Thickness Retrieval

Most of the known methods and corresponding atmospheric correction algorithms have been designed to determine characteristics of aerosol over dark areas of the Earth’s surface such as vegetation. These surfaces have two characteristic features. The first feature is low reflectance values in the shortwave region. In particular, in the wavelength region about $\lambda \approx 400\mathrm{nm}$ surface reflectance values of the order of 0.01 − 0.04 are characteristic for different types of vegetation (Figure 2). The second feature, in fact, is a consequence of the first one. With such a reflectance, the modulus of the derivative dr_surf/dτ_aer, which determines the error of surface reflectance retrieval depending on the error in the AOT, is quite significant, namely $\left|d{r}_{surf}/d{\tau }_{aer}\right|~0.2$ (Figure 3). This means that an error $\left|\Delta r_{surf}\right|~0.01÷0.02$ leads to an error $\left|\Delta {\tau }_{aer}\right|~0.05÷0.1$ in determining the AOT of the atmosphere. This seems to be the main reason for the successful determination of the atmospheric AOT over dark surfaces from multispectral data. Conversely, an error in determining the AOT of the atmosphere $\left|\Delta {\tau }_{aer}\right|~0.05÷0.1$ leads to an error $\left|\Delta r_{surf}\right|~0.01÷0.02$ in surface reflectance retrieving. When moving into the long-wave region, where the surface reflectance usually increases, the error value $\left|d{r}_{surf}/d{\tau }_{aer}\right|$ is noticeably smaller. In the wavelength interval where the derivative $d{r}_{surf}/d{\tau }_{aer}$ is close to zero, even a small error in estimating the surface reflectance leads to a large error in determining the atmospheric AOT. Such a strong sensitivity of the error $\Delta {\tau }_{aer}$ to the magnitude $\Delta r_{surf}$ is due to the fact that in this region the surface reflectance (of the order of $0.15÷0.20$) is close to the reflectance of the aerosol layer of semi-infinite thickness. Therefore, even a significant change in the atmospheric AOT weakly changes the total reflectance of the atmosphere + underlying surface system.

Bright surfaces (deserts, various sands) have an reflectance of the order of $0.15÷0.30$ in the vicinity of $\lambda \approx 400\mathrm{nm}$ (Figure 2). The modulus of derivative $d{r}_{surf}/d{\tau }_{aer}$ is very small, of the order $0.01÷0.02$ (Figure 3). Therefore, even small errors in the estimation of the surface reflectance can lead to large errors in determining the aerosol optical thickness of the atmosphere. This circumstance makes to be very problematic successful retrieval atmospheric AOT over bright surfaces in algorithms based on the processing only spectral data. However, this same circumstance makes it possible to retrieve the spectral reflectance of the earth’s surface with satisfactory accuracy, since the error of the reflectance retrieval is weakly sensitive to the error of the atmospheric AOT. Indeed, in the case of a bright Earth’s surfaces, even an error $\Delta {\tau }_{aer}~0.3÷0.5$ should lead to an error of the order of 0.01 − 0.02 in the retrieved surface reflectance in the short-wave region of the spectrum.

From the above, three important conclusions follow that determine the main feature of the RACE algorithm:

• (1). The minimum sensitivity of error $\left|\Delta {\tau }_{aer}\right|$ to error $\left|\Delta r_{surf}\right|$ occurs for dark surfaces in the short wavelength region of the spectrum. Therefore, it is optimal to determine the AOT above dark surfaces in the short wavelength region of the spectrum.

• (2). The sensitivity of error $\left|\Delta r_{surf}\right|$ of retrieved reflectance of dark surfaces to AOT error is at maximum in the short-wave region of the spectrum. When the surface reflectance increases with the wavelength, the retrieving error decreases even if the AOT of the atmosphere is determined with the same error.

• (3). For bright surfaces, the sensitivity of surface reflectance errors $\left|\Delta r_{surf}\right|$ retrieving to the errors $\left|\Delta {\tau }_{aer}\right|$ is significantly less in almost the entire spectral range.

The above conclusions and the narrow range of dark surfaces reflectance variations in the short-wavelength region of the spectrum, noted above, allow us to propose the following:

• to retrieve the AOT from the SR recorded by a satellite sensor at one wavelength in this region, in particular, at a wavelength of 412 nm.

• the surface reflectance (in the absence of additional information) can be taken to be fixed and, for example, equal to the average value ${\overline{r}}_{surf}\left(412\right)\approx 0.028$ for a sample of dark surfaces in Figure 2.

In what follows the atmospheric AOT retrieval with use of these estimations is called single-wavelength retrieval.

We note one more argument in favor of using the shortest wavelength for the retrieval of the atmospheric AOT. Just in this spectral interval, the atmosphere aerosol manifests itself most strongly since it is here where the atmospheric AOT is greatest, and the surface reflectance is minimal.

The single-wavelength method for determining the AOT has some significant advantages in problems of atmospheric correction. The ability to determine the atmospheric AOT without reference to a specific type of underlying surface spectrum makes it possible to use this method even in cases with different types of underlying surfaces, including a water surface, inhomogeneous surfaces (for example, a city), with a mixture of different types of surfaces, etc. Finally, the undoubted advantage of this method is its efficiency, since it does not require complex iterative procedures with large expenditures of computer time.

Moreover, the use of additional wavelengths in some cases can even reduce the accuracy of determining the atmospheric AOT due to the smallness of the derivative $\left|d{r}_{surf}/d{\tau }_{aer}\right|$ and, accordingly, increasing values $\left|d{\tau }_{aer}/d{r}_{surf}\right|$ with growth of the surface reflectance in the long-wavelength region of the spectrum.

It should be noted that if the real surface reflectance at the minimum wavelength noticeably exceeds the selected fixed value, then the retrieved atmospheric AOT value could noticeably exceed its true value. This situation can occur, in particular, in the case of a bright underlying surface.

Unfortunately, when the single-wavelength retrieval algorithm is operating in an automatic mode, it is impossible to determine the true cause of the retrieved high AOT value in the shortwave region, namely, whether it is due to strong atmosphere turbidity or is associated with a high surface reflectance. Therefore, it is proposed to introduce the maximum permissible value of the AOT equal to 0.5 at the wavelength $\lambda =550\mathrm{nm}$. This means that if for the adopted aerosol model the retrieved AOT value at a wavelength of 550 nm exceeds 0.5, it should be taken equal to 0.5.

The rationale for such a choice of the permissible AOT value is the analysis of the real AOT values of the atmosphere in different regions of the Earth according to MODIS measurements [15]. In most regions of the world, AOT at the wavelength $\lambda =550\mathrm{nm}$ does not exceed 0.5. At the same time, the average AOT(550 nm) over 15,894 measurements around the world is 0.2 [15].

If the real reflectance of the surface at the wavelength of 412 nm is noticeably less than 0.028, then the opposite situation may take place. The value of the AOT retrieved may turn out to be noticeably less than the true value and even less than zero. In this regard, a second artificial limitation is introduced. If the retrieved AOT value at the wavelength of 550 nm turns out to be less than 0.05, then its value is set equal to 0.05. The rationale for this limitation is the fact that, firstly, the actual atmospheric AOT less than 0.05 is very rare, and, secondly, this value corresponds to the minimum root-mean-square error at a wavelength of 550 nm in the BAR algorithm for processing MODIS data [1].

3. Retrieval of the Spectral Reflectance of the Underlying Surface

If a pixel size of the image is large compared to the variance of point spread function (PSF) in the atmosphere, or the averaging is performed over a large area (compared to those of the PSF variance), then the average surface reflectance ${\overline{r}}_{surf}$ is defined from the following well-known equation for atmospheric correction [6]:

(9)$R_{TOA}\left(\lambda ,\mu ,{\mu }_0,\phi \right)=R_{atm}\left(\lambda ,\mu ,{\mu }_0,\phi \right)+\frac{t\left(\lambda ,\mu \right)t\left(\lambda ,{\mu }_0\right)}{1-r_{atm}^{*}\left(\lambda \right){\overline{r}}_{surf}\left(\lambda \right)}{\overline{r}}_{surf}\left(\lambda \right),$

where $R_{atm}\left(\lambda ,\mu ,{\mu }_0,\phi \right)$ is the TOA reflectance of the whole atmosphere above a black underlying surface, ${\overline{r}}_{surf}\left(\lambda \right)$ and ${\overline{R}}_{TOA}\left(\lambda ,\mu ,{\mu }_0,\phi \right)$ are surface reflectance and TOA reflectance, averaged over the area larger than PSF of the atmosphere, $t\left(\lambda ,{\mu }_0\right)$ is the total (direct+diffuse) transmittance of the whole atmosphere, $r_{atm}^{*}\left(\lambda \right)$ is the spherical reflectance of the atmosphere at illumination of atmosphere bottom upwards. We emphasize that in Equation (9) and further the spectral reflectance of the underlying surface is considered to be Lambertian.

For high spatial resolution sensors, the basic equation for atmospheric correction can be written as:

(10)$\begin{array}{l}{R}_{TOA,p}\left(\lambda ,\mu ,{\mu }_0,\phi \right)=R_{atm}\left(\lambda ,\mu ,{\mu }_0,\phi \right)\\ +\frac{t\left(\lambda ,{\mu }_0\right)}{1-r_{atm}^{*}\left(\lambda \right){\overline{r}}_{surf}\left(\lambda \right)}\left[t^s\left(\lambda ,\mu \right){\overline{r}}_{surf}\left(\lambda \right)+t^0\left(\lambda ,\mu \right)r_p\left(\lambda \right)\right],\end{array}$

or:

(11)$R_{TOA,p}\left(\lambda ,\mu ,{\mu }_0,\phi \right)={\overline{R}}_{TOA}\left(\lambda ,\mu ,{\mu }_0,\phi \right)+\frac{t^0\left(\lambda ,\mu \right)t\left(\lambda ,{\mu }_0\right)}{1-r_{atm}^{*}\left(\lambda \right){\overline{r}}_{surf}\left(\lambda \right)}\left[r_p\left(\lambda \right)-{\overline{r}}_{surf}\left(\lambda \right)\right],$

where:

(12)${\overline{R}}_{TOA}\left(\lambda ,\mu ,{\mu }_0,\phi \right)=R_{atm}\left(\lambda ,\mu ,{\mu }_0,\phi \right)+\frac{t\left(\lambda ,\mu \right)t\left(\lambda ,{\mu }_0\right)}{1-r_{atm}^{*}\left(\lambda \right){\overline{r}}_{surf}\left(\lambda \right)}{\overline{r}}_{surf}\left(\lambda \right).$

Here $R_{TOA,p}\left(\lambda ,\mu ,{\mu }_0,\phi \right)$ is the TOA reflectance when observing a certain pixel “p”, $t0\left(\lambda ,\mu \right)$ and $t^s\left(\lambda ,\mu \right)=t\left(\lambda ,\mu \right)-t0\left(\lambda ,\mu \right)$ are direct and diffuse components of the atmospheric transmittance, respectively, $r_p\left(\lambda \right)$ is the reflectance of the pixel “p”.

Equation (10) is derived under assumptions that the Earth’s surface reflects incident light according to Lambert’s law, the area of “p” pixel is much smaller than the PSF in the atmosphere and the contribution of the scattered light to the TOA reflectance is the same for all pixels over which averaging is performed.

It follows from Equation (11):

(13)$r_p\left(\lambda \right)=\frac{R_{TOA,p}\left(\lambda ,\mu ,{\mu }_0,\phi \right)-{\overline{R}}_{TOA}\left(\lambda ,\mu ,{\mu }_0,\phi \right)}{D\left(\lambda ,\mu ,{\mu }_0\right)}+{\overline{r}}_{surf}\left(\lambda \right),$

where:

(14)$D\left(\lambda ,\mu ,{\mu }_0\right)=\frac{t^0\left(\lambda ,\mu \right)t\left(\lambda ,{\mu }_0\right)}{1-r_{atm}^{*}\left(\lambda \right){\overline{r}}_{surf}\left(\lambda \right)}.$

The spectral widths of the high spatial resolution channels in the MZ systems could be quite large (say, 10–30 nm) to provide the high enough values of the signal-to-noise ratio. Therefore, one must account for the spectral behavior of optical properties of atmosphere and underlying surface in the MZ channels during the solution of the inverse problem related to the underlying spectral surface reflectance determination.

Let us integrate Equation (11) with respect to the wavelength $\lambda$ with weight $\left({\mu }_0/\pi d^2\right)S_0\left(\lambda \right)K_n\left(\lambda \right)$. Then it follows:

(15)$J_{TOA,p,n}={\overline{J}}_{TOA,n}+\frac{{\mu }_0}{\pi d^2}{\displaystyle \int S_0\left(\lambda \right)K_n\left(\lambda \right)\frac{t0\left(\lambda ,\mu \right)t\left(\lambda ,{\mu }_0\right)}{1-r_{atm}^{*}\left(\lambda \right){\overline{r}}_{surf}\left({\lambda }_n\right)}\left[r_p\left(\lambda \right)-{\overline{r}}_{surf}\left(\lambda \right)\right]d\lambda },$

where:

(16)$J_{TOA,p,n}=\frac{{\mu }_0}{\pi d^2}{\displaystyle \int S_0\left(\lambda \right)K_n\left(\lambda \right)R_{TOA,p}\left(\lambda ,\mu ,{\mu }_0,\phi \right)d\lambda }$

is the spectral radiance at the TOA in the pixel “p”:

(17)${\overline{J}}_{TOA,n}=\frac{{\mu }_0}{\pi d^2}{\displaystyle \int S_0\left(\lambda \right)K_n\left(\lambda \right){\overline{R}}_{TOA}\left(\lambda \right)d\lambda }$

the spectral radiance in the n-th channel, averaged over the area larger than the PSF of the atmosphere, d is the distance from the Earth to the Sun in the astronomical units, $S_0\left(\lambda \right)$ is the solar light irradiance at the TOA at $d=1$, $K_n\left(\lambda \right)$ is the spectral sensitivity of the n-th channel with the normalization condition:

(18)${\displaystyle \int K_n\left(\lambda \right)d\lambda }=1.$

It follows from Equation (16) that one cannot determine the real reflectance of the test pixel “p” using the MZ satellite data. However, one can derive the so-called effective reflectance averaged with respect to the wavelength taking into account the width of the n-th channel as specified below:

(19)$r_{p,n}=\frac{{\displaystyle \int S_0\left(\lambda \right)K_n\left(\lambda \right)\frac{t0\left(\lambda ,\mu \right)t\left(\lambda ,{\mu }_0\right)}{1-r_{atm}^{*}\left(\lambda \right){\overline{r}}_{surf}\left(\lambda \right)}{r}_p\left(\lambda \right)d\lambda }}{{\displaystyle \int S_0\left(\lambda \right)K_n\left(\lambda \right)\frac{t^0\left(\lambda ,\mu \right)t\left(\lambda ,{\mu }_0\right)}{1-r_{atm}^{*}\left(\lambda \right){\overline{r}}_{surf}\left(\lambda \right)}d\lambda }}.$

One can see that the effective reflectance depends on the spectral sensitivity of the receiver $K_n\left(\lambda \right)$.

Based on Equation (15), it can be shown that the effective reflectance $r_{p,n}$ of the pixel “p” is related to the spectral radiance $J_{TOA,p,n}$ in this channel by a linear coupling:

(20)$r_{p,n}=\frac{J_{TOA,p,n}-{\overline{J}}_{TOA,n}}{D_n}+{\overline{r}}_{surf,n},$

where:

(21)$D_n=\frac{{\mu }_0}{\pi d^2}\frac{{\displaystyle \int S_0\left(\lambda \right)K_n\left(\lambda \right)t0\left(\lambda ,\mu \right)t\left(\lambda ,{\mu }_0\right)d\lambda }}{1-r_{atm,n}^{*}{\overline{r}}_{surf,n}}=\frac{F_n(t0t)}{1-r_{atm,n}^{*}{\overline{r}}_{surf,n}},$

(22)${\overline{r}}_{surf,n}=\frac{{\overline{J}}_{TOA,n}-F_n\left(R_{atm}\right)}{F_n\left(tt\right)+r_{atm,n}^{*}\left[{\overline{J}}_{TOA,n}-F_n\left(R_{atm}\right)\right]}$

(23)$r_{atm,n}^{*}\approx \frac{F_n\left(tt{r}_{atm}^{*}\right)}{F_n\left(tt\right)}$

(24)$F_n(t0t)=\frac{{\mu }_0}{\pi d^2}{\displaystyle \int S_0\left(\lambda \right)K_n\left(\lambda \right)t0\left(\lambda ,\mu \right)t\left(\lambda ,{\mu }_0\right)d\lambda }$

(25)$F_n\left(R_{atm}\right)=\frac{{\mu }_0}{\pi d^2}{\displaystyle \int S_0\left(\lambda \right)K_n\left(\lambda \right)R_{atm}\left(\lambda \right)d\lambda }$

(26)$F_n\left(tt\right)=\frac{{\mu }_0}{\pi d^2}{\displaystyle \int S_0\left(\lambda \right)K_n\left(\lambda \right)t\left(\lambda ,\mu \right)t\left(\lambda ,{\mu }_0\right)d\lambda }$

(27)$F_n\left(tt{r}_{atm}^{*}\right)=\frac{{\mu }_0}{\pi d^2}{\displaystyle \int S_0\left(\lambda \right)K_n\left(\lambda \right)t\left(\lambda ,\mu \right)t\left(\lambda ,{\mu }_0\right)r_{atm}^{*}\left(\lambda \right)d\lambda }.$

Thus, to restore the spectral reflectance of the underlying surface in the pixels of the MZ system of a high spatial resolution, it is necessary to calculate four transfer functions of the atmosphere, namely $F_n\left(R_{atm}\right)$, $F_n\left(tt\right)$, $F_n\left(t^{0}t\right)$ and $F_n\left(tt{r}_{atm}^{*}\right)$.

As seen from relations (24)–(27), with the chosen model of the molecular-aerosol atmosphere, these transfer functions depend on the AOT and content of water vapor in the atmosphere. It is emphasized, that the calculation of these functions is the most laborious. Therefore, the number of such calculations is desirable to minimize. Given that the horizontal scale of changes in the water vapor content in the atmosphere, as a rule, is tens of kilometers, the concentration of water vapor in the atmosphere can be considered the same throughout the frame. In this case, the only atmospheric parameter changing within the image frame is the AOT of the atmosphere.

Therefore, it seems appropriate to calculate once the functions $F_n\left(R_{atm}\right)$, $F_n\left(tt\right)$, $F_n\left(t^{0}t\right)$ and $F_n\left(tt{r}_{atm}^{*}\right)$ for a certain set of AOT values, form them into a Look-up tables, and subsequently, after retrieving the AOT for the local image area, to retrieve the spectral reflectance of the underlying surface in pixels for the corresponding case based on these Look-up tables using interpolation.

The described procedure was deployed for retrieving spectral reflectance of underlying surfaces in the MZ channels in cases presented in Section 4.

4. Testing the RACE Algorithm

4.1. Testing Method

As noted above, in the process of atmospheric correction of multispectral satellite data when retrieving the spectral characteristics of the underlying surface, unknown optical parameters of the atmosphere, and, in particular, atmospheric AOT, must be determined from the same spectral measurements. It is clear that such a task cannot be solved without some a priori assumptions. Therefore, the purposes of performed testing, which description and results are given below, were:

• to verify the correctness of a priori assumptions used in the RACE algorithm,

• to evaluate the stability and speed of the algorithm,

• to estimate the accuracy of retrieving atmosphere parameters and reflectance of the underlying surface with a fairly wide variation of satellite imagery conditions.

The main test results were obtained using computer simulation of signals recorded by a satellite sensor. The software for simulating these signals was developed and carefully tested. The radiative transfer in the system atmosphere—underlying surface has been simulated with the reliable and accurate code RAY [9]. Besides testing based on available experimental data was carried out. When analyzing the RACE algorithm, four main possible sources of error in retrieving the atmospheric optical parameters and the spectral reflectance of the underlying surface can be identified:

• -. error of the single-wavelength method to retrieve AOT due to the difference between the real (unknown in practice) surface reflectance and its fixed value at wavelength of 412 nm used in the RACE algorithm;

• -. error of radiometric calibration of sensors in the MS channels;

• -. error of radiometric calibration of sensors in the MZ channels;

• -. mismatch between the a priori model of atmospheric aerosol used and the “real” atmospheric aerosol.

For testing, a basic set of parameters of the atmosphere, underlying surface, and observation geometry was selected the most typical for satellite imagery of the Earth’s surface. This set includes:

• -. three values of aerosol optical thickness of troposphere equal to 0.1, 0.3 and 0.5 at the wavelength of 550 nm;

• -. 12 types of underlying surface, seven of which relate to dark surfaces such as vegetation (green vegetation (GV), oat, grass, rye, beet, lawn, soil), and five surfaces are bright (desert, sand, loam, dry and wet concrete);

• -. two values of the solar zenith angle: 20° and 60°;

• -. viewing zenith angle 0°;

• -. atmospheric water vapor column $W=2g\cdot {\mathrm{cm}}^{-2}$.

For testing, a fixed model of the molecular-gas atmosphere “MidlatitudeSummer” [16], typical for the middle latitudes of the globe, was chosen. The widely used model “Continental” [17], consisting of three fractions with lognormal distributions of particle sizes (29% Water Soluble, 70% Dust, and 1% Soot), was taken as a model of aerosol in the lower troposphere layer. To determine the water vapor content, we used the spectral channels 910–920 nm (absorption band of water vapor) and 865–875 nm (reference channel outside the absorption band). The reflectance of the underlying surface was retrieved in the narrow spectral channels of the MS at wavelengths of 412 nm, 442 nm, 489 nm, 509 nm, 559 nm, 619 nm, 664 nm, 776 nm and 867 nm, and in the MZ spectral channels. The spectral ranges for various MZ channels, the effective wavelengths in them and the parameters of the a priori atmospheric model used for testing are given in Table 1. Table 2 and Table 3 present ranges of reflectance variation of dark and bright surfaces in the MS and MZ spectral channels.

To clarify the role of the main sources of errors of the RACE, algorithm testing was carried out both for the base case, when the a priori models used fully corresponded to the "real" conditions, and for the cases when the models used and the "real" conditions do not coincide. The influence of all four sources of errors noted above was considered. We emphasize that in all cases the retrieval of parameters of the atmosphere and surface was carried out with the same a priori assumptions. For retrieving AOT, the used earth’s surface reflectance at 412 nm was 0.028. The test results are presented in Figure 4, Figure 5, Figure 6, Figure 7, Figure 8 and Figure 9, in Table 4, Table 5, Table 6, Table 7, Table 8, Table 9, Table 10, Table 11, Table 12 and Table 13 and discussed in what follows.

4.2. The Influence of the Error of the Atmospheric AOT Retrieval with The Single-Wavelength Method

First, let us estimate the errors in retrieving the parameters of the atmosphere and the reflectance of the underlying surface in the basic version, when the a priori models used fully correspond to the “real” conditions. Thereby, we estimate the error in the determination of the atmospheric AOT due to the difference between the real (unknown in practice) reflectance of the surface and its fixed value at a wavelength of 412 nm used in the RACE single-wavelength algorithm.

Table 4 presents the results of retrieval AOT at $\lambda =550\mathrm{nm}$ for a few underlying surface models and the root mean square errors (RMSE) for three values of the real AOT and solar zenith angles 60° and 20°. For comparison, the last column of Table 4 shows the results of estimation the RMSE at a wavelength of 550 nm in the BAR algorithm for processing MODIS data [1].

Note that in the cases of a high surface reflectance at the wavelength 412 nm, retrieved AOT values almost in all cases turned out to be 0.5. Therefore, we do not present the results of a similar test for the case of bright surfaces.

Figure 4 and Figure 5 illustrate the accuracy of retrieving reflectance of dark and bright underlying surfaces when determining the AOT by the proposed method. One can see that for dark surfaces the error weakly depends on the solar zenith angle and does not exceed 0.02. The maximum error occurs at the shortest wavelength of 412 nm and decreases with increasing wavelength. For bright surfaces, for which the surface reflectance is significantly higher, the reflectance retrieval error is also greater, especially at the solar zenith angle 60°. It is obvious that in this case the error in retrieval AOT manifests itself more pronounced. The errors of the surface reflectance retrieval calculated in the MS spectral channels from these data are presented in Table 5. Note that RMSE do not exceed 0.01 in the case of dark surfaces and 0.025 in the case of bright surfaces.

Errors of the spectral reflectance retrieval in the MZ channels are presented in Figure 6. Naturally, in this case as well, the maximum error of surface reflectance retrieval for dark surfaces occurs in the short-wave MZ channel, wherein it does not exceed 0.015. Table 6 presents RMSE for the surfaces’ reflectance retrieval computed with these data. These errors are less than 0.01 for dark surfaces, about 0.02–0.03 for bright surfaces at the solar zenith angle of 60°, and about 0.01–0.02 at the solar zenith angle of 20°.

4.3. Influence of Errors of Sensors Radiometric Calibration in the MS and MZ Spectral Channels

Errors in radiometric calibration of sensors in the MS spectral channels affect the results of atmospheric correction in two ways. First, the calibration error at a wavelength of 412 nm leads to errors in the determination of the atmospheric AOT and, as a consequence, to errors in the retrieval of the surface reflectance both in the MS and the MZ channels. In addition, calibration errors in the MS spectral channels affect the accuracy of the retrieval of the water vapor content in the atmosphere. The errors in the sensor calibration in the MZ spectral channels directly affect the accuracy of retrieving the surface reflectance in these channels. Below, when analyzing the test results, the calibration error of the MS and MZ sensors in all spectral channels is supposed to be the same and equal to 5%.

Figure 7 clearly shows that overestimation of the MS calibration at a wavelength of 412 nm and, as a result, overestimation of retrieved AOT value, leads, as one would expect, to underestimation of the retrieved values of the spectral surface reflectance in the MZ channels for dark surfaces and to overestimation for bright surfaces. Table 7 displays the RMSE in retrieved values of surface reflectance in the MZ channels.

Figure 8 demonstrates errors in the retrieved values of surfaces spectral reflectance due to the error in the calibration of sensors in the MZ spectral channels. In addition, the RMSE of retrieved surface reflectance calculated from these data, are presented in Table 8.

4.4. The Effect of Discrepancy between the Used a Priori and “Real” Models of Atmospheric Aerosol

One of the main sources of possible errors in the RACE algorithm products is the mismatch between the used a priori and the “real” models of atmospheric aerosols. The fact is that atmospheric correction of MS satellite data does not allow to determine parameters of “real” atmospheric aerosols. Therefore, in the automatic express version of atmospheric correction proposed here, it is necessary to use a fixed, a priori selected aerosol model, for which we recommend using the described above model with aerosol “Continental” for lower and middle tropospheric layers. In this regard, the question of the error caused by the mismatch between the a priori model of aerosol used and the “real” atmospheric aerosol in the case seems to be non-trivial and especially important.

In a presented computer test, we examined atmospheric correction with three different models of “real” aerosol in the lower tropospheric layer, namely “Belarus”, “Maritime” and “Oceanic” [17]. The signals recorded by satellite sensors at TOA in these cases were simulated with the RAY code [9]. Let us underline that in all cases, the retrieval of the atmosphere parameters and reflectance of underlying surfaces was carried out with the a priori selected aerosol model “Continental”.

Let us explain this choice of “real” aerosol models.

The “Belarus” model is based on statistical data on the microphysical parameters of atmospheric aerosol for the spring-summer-autumn period for the territories of the Republic of Belarus and Poland. It contains relatively small particles. The “Maritime” and “Oceanic” models, on the contrary, contain large watered salt particles and describe the scattering of light in the troposphere over oceans. Thus, these selected models include two extreme cases, and the a priori model “Continental” corresponds to intermediate situations. Table 9 lists some important optical and microphysical parameters of these aerosol models, namely, the effective particle size of the aerosol particles $r_{eff}$, single scattering albedo (SSA), the phase function value at an angle of 120°, $x\left({120}^0\right)$, and the Angstrom parameter $\alpha$ determining the spectral dependence of the atmospheric AOT, ${\tau }_{aer}\left(\lambda \right)~{\lambda }^{-\alpha }$.

An analysis of the test results showed that the errors in spectral reflectance retrieving in cases of “real” “Maritime” and “Oceanic” aerosol models are approximately the same. Therefore, below we restrict ourselves to reviewing the test results only for the “Belarus” and “Oceanic” models. Table 10 presents the RMSE of the AOT retrieval for different models of the dark underlying surface. The RMSE in retrieving spectral reflectance of underlying surfaces in the MZ channels for cases when “Belarus” and “Oceanic” aerosol models are used as “real” atmospheric aerosols are displayed in Figure 9 and Table 11.

It can be seen from Table 10 that the use of the “Continental” aerosol model instead of “real” aerosol models “Belarus” and “Oceanic” in the atmospheric correction procedure can lead to significant errors in the retrieved values of AOT at a wavelength of 550 nm. Added to this is the difference in the spectral behavior of “real” and retrieved AOT values. However, it is a matter of fact that the comparison of Table 11 and Table 6 shows that this discrepancy does not lead to a significant increase in the error of retrieved value of reflectance of the underlying surface (especially in the cases with dark surfaces), even despite the fundamentally different AOT spectral behavior of the “Continental” and “Oceanic” models.

4.5. Resulting Errors in Retrieval of the Underlying Surface Reflectance in the MZ Channels. The Contribution of Various Factors

Let us assume that the contributions of various sources of error to the resulting error in retrieved value of the surface reflectance are independent. In this case, the RMSE can be calculated by the formula:

(28)$\Delta r_{peз}={\left({\displaystyle \sum _i\Delta r_i}\right)}^{0.5},$

where $\Delta r_i$ are the RMSE of the main sources of error.

Table 12 and Table 13 present resulting RMSE in retrieved surface reflectance in the MZ channels and the contribution from various sources of errors. Note, that the error averaged over models “Belarus” and “Oceanic” was accepted as an error due to mismatch between the a priori model of atmospheric aerosol used and the “real” atmospheric aerosol. It can be seen that the resulting error of the spectral reflectance retrieval in the MZ channels is approximately 0.01–0.02 in the case of a dark surface, while in the case of a bright surface it is of the order of 0.03–0.04. Moreover, in the case of a dark underlying surface, just inaccuracies of radiometric calibration made the main contribution to the resulting error, while in the case of a bright surfaces, errors of single-wavelength method, of the MZ sensors calibration, and of atmospheric model are about the same order. The final rows in Table 12 and Table 13 display the resulting errors in the absence of a calibration error in the MZ channels. As seen, in this case, the error in the retrieving reflectance of dark surfaces becomes noticeably smaller, and that of bright surfaces changes slightly.

Table 14 presents a comparison of the accuracy of the RACE algorithm according to the data of Table 12 and Table 13 with the accuracy of the known algorithms Fast Line-of-sight Atmospheric Analysis of Spectral Hypercubes (FLAASH) [18] of the ENVI package and Digital Globe Atmospheric Compensation (DG AComp) [4] according to [19].

It should be emphasized that the test results of the RACE algorithm presented in the paper were obtained for the case of the algorithm operating in a (blind) automatic mode. It can be expected that when operating in operator mode, when it is possible to use information about the type of underlying surface and the type of aerosol in the survey area, the accuracy of determining the surface reflectance may be higher. This primarily applies to areas with a bright surface, where the reflectance at a wavelength of 412 nm is much higher than the a priori accepted ${\overline{r}}_{surf}\left(412\mathrm{nm}\right)=0.028$.

4.6. Testing the RACE Algorithm Using Experimentally Measured Ground Reflectance Data

Testing was carried out using measurement data on three test sites ((La Crau, France; Railroad Valley, NV, USA; Gobabeb, Namibia) with three different types of underlying surfaces, obtained by the Radiometric Calibration Network (RadCalNet) [20]. In all three cases, ground-based measurements provided spectral surface reflectance, temperature and atmospheric pressure, ozone and water vapor content in the atmosphere, and AOT at the wavelength of 550 nm. In addition, the surface height above sea level and the observation geometry are known [20]. Table 15 presents some of the most important data on the test sites and experimental conditions.

The surface reflectance data were converted to top-of-atmosphere reflectance within RadCalNet and provided through a web portal to allow users to either radiometrically calibrate or verify the calibration of their sensors of interest. These converted data were used for testing the RACE algorithm operating in automatic mode. This means that in all three cases retrieving atmosphere parameters and reflectance of underlying surfaces was carried out with the same a priori basic atmospheric model described above (latitudinal model of the molecular atmosphere, three-layer model of the aerosol atmosphere with aerosol “Continental” in the lower and middle tropospheric layers). Beside while determining the atmospheric AOT, the underlying surface reflectance at the wavelength of 412 nm was assumed equal to 0.028. Additionally, real data on the surface height and observation geometry, as well as on the content of ozone and water vapor were used.

The results of retrieving the spectral reflectance of the surface by the RACE algorithm in comparison with the data of real ground measurements are presented in Figure 10 and in Table 16. Besides, the Table 16 displays the measured and retrieved atmospheric AOT values and errors of surface reflectance retrieval calculated as the difference between the retrieved and measured values. It can be seen that despite considerable errors in retrieving AOT over bright surfaces (Railroad Valley, NV, USA and Gobabeb, Namibia) the surface reflectance is retrieved in all cases with acceptable accuracy.

As for significant errors in retrieving values of small AOTs over bright surfaces, let us note that this feature is characteristic for all existing algorithms. Physical reasons have been described above. However, it may be worth recalling that the AOT retrieval is not the objective point for the RACE algorithm, aimed at the effective retrieval of a surface spectral reflectance in an automatic mode in situations with lack of information about underlying surface.

Table 12, Table 13 and Table 16 show that the results of testing the RACE algorithm according to experimental data and computer simulation data are in a good agreement.

5. Conclusions

In this paper we have presented a simple and robust technique for retrieving surface spectral reflectance using spectral reflectance measurements captured by satellite instruments with high or moderate spatial resolution and with limited amount of available auxiliary information.

The key point of this algorithm is the suggested single-wavelength method for estimating the atmospheric AOT at 400 nm. This method is mainly based on two factors.

First, low spectral ground reflectance values and its low variability (of the order of 0.01 − 0.04) in vicinity $\lambda \approx 400\mathrm{nm}$ are characteristic for a wide class of dark areas (vegetation types) of the Earth’s surface.

With such surface reflectance values, the dependence of the error $\left|\Delta {\tau }_{aer}\right|$ in the AOT estimation on the error $\left|\Delta r_{surf}\right|$ in specifying the surface reflectance has a minimum.

Second, although for bright surfaces, such as deserts and semi-deserts, errors in AOT retrieval may be noticeably greater than in the previous case, the sensitivity of surface reflectance errors $\left|\Delta r_{surf}\right|$ to the errors $\left|\Delta {\tau }_{aer}\right|$ is significantly smaller in almost entire spectral range.

Other important classes of objects are seas and oceans. Considering that the reflectance of water basins in the region $\lambda \approx 400\mathrm{nm}$ is approximately the same as that of vegetation, the proposed technique is promising for satellite monitoring of water basins as well.

The retrievals, performed with computer simulations and use of experimental data, have demonstrated a sufficient retrieval accuracy.

Thus, the proposed technique can be used to retrieve the reflectance of various underlying surfaces in an automatic mode without huge databases accumulated in advance.

Author Contributions

Conceptualization, I.L.K. and A.S.P.; methodology, I.L.K. and E.P.Z., software, A.S.P.; validation, I.L.K. and A.S.P.; formal analysis, I.L.K.; investigation, I.L.K., A.S.P.; resources, I.L.K.; data curation, A.S.P.; writing—original draft preparation, I.L.K., A.S.P., E.P.Z. and A.A.K.; writing—review and editing, I.L.K. and E.P.Z.; visualization, A.S.P.; supervision, I.L.K. All authors have read and agreed to the published version of the manuscript.

Funding

This research received no external funding.

Institutional Review Board Statement

Not applicable.

Informed Consent Statement

Not applicable.

Data Availability Statement

Data is available on request.

Acknowledgments

We are grateful to a Radiometric Calibration Network (RadCalNet) for Earth observing imagers operating in the visible to shortwave infrared spectral range.

Conflicts of Interest

The authors declare no conflict of interest.

Publisher’s Note: MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Figures and Tables

Figure 1. RTOA412 nm and the difference RTOA412 nm−Rmol412 nm under different observation conditions: solar zenith angles from 20 to 80 degrees, viewing angles 0, 15 and 30 degrees, azimuths of observation 45, 90, 135 and 180 degrees. Earth’s surface is sand.

Figure 2. Spectral reflectance of certain types of vegetation (a) and of sand and desert (b).

Figure 3. Reflectance spectra (green line) and derivative drsurf/dτaer for grass (a) and sand (b) at various values of AOT.

Figure 4. Reflectance of dark surfaces ((a) and (c)) and bright surfaces ((b) and (d)): true (solid lines with badges) and retrieved with three values of AOT (0.1-short strokes, 0.3-long strokes and 0.5-dash-dot), as well as the SRs at the top-of-atmosphere at AOT 0.03 (solid lines without icons).

Figure 5. Errors in retrieving surface spectral reflectance in the MS channels for dark surfaces (a) and bright surfaces (b) at solar zenith angles 60° (red circles) and 20° (green triangles).

Figure 6. Errors in retrieving surface spectral reflectance in the MZ channels for dark surfaces (a) and bright surfaces (b) at solar zenith angles 60° (red circles) and 20° (green triangles).

Figure 7. Errors in retrieving surface spectral reflectance in the MZ channels for dark surfaces (a) and bright surfaces (b) at solar zenith angles 60° (red circles) and 20° (green triangles). Radiometric calibration error in the MS channel 412 nm is 5%.

Figure 8. Errors in retrieving surface spectral reflectance in the MZ channels for dark surfaces (a) and bright surfaces (b) at Sun zenith angles 60° (red circles) and 20° (green triangles), an error of radiometric calibration in the MZ HR channels being 5%.

Figure 9. Errors in retrieving spectral reflectance of dark surfaces (a,c) and bright surfaces (b,d) in the MZ channels in the cases with “Belarus” (a,b) and ”Oceanic” (c,d) models as a true “real” atmospheric aerosol. Zenith angles of the solar are 60° (red circles) and 20° (green circles).

Figure 10. Spectral reflectance of the surface at test sites La Crau, France (LCFR), Railroad Valley, NV, USA (RVUS), and Gobabeb, Namibia (GONA) according to ground measurements (solid lines) and data retrieved with RACE algorithm (dashed lines).

Table 1

Characteristics of MZ channels and the molecular ${\tau }_{mol}$ and aerosol ${\tau }_{aer}$ optical thicknesses of the a priori atmosphere model used in testing.

Band	Range of Wavelengths at the Level of 0.5, nm	EffectiveWavelength,nm	${\mathit{\tau }}_{\mathit{m}\mathit{o}\mathit{l}}$	${\mathit{\tau }}_{\mathit{a}\mathit{e}\mathit{r}}$ in the Troposphere Layer at $\mathbf{\lambda }\mathbf{=}\mathbf{550}\mathbf{nm}$
				0.1	0.3	0.5
Band 1	470–510	490.5	0.1593	0.114	0.341	0.568
Band 2	510–580	548.5	0.0927	0.100	0.301	0.502
Band 3	650–690	670.5	0.0427	0.080	0.240	0.401
Band 4	780–910	841.5	0.0233	0.062	0.187	0.311

Table 2

Reflectance ranges of dark and bright surfaces in the MS spectral channels.

SurfaceType	Wavelength, nm
	412.5	442.5	490	510	560	620	665	778.5
Dark	0.01–0.043	0.013–0.05	0.02–0.063	0.025–0.07	0.04–0.09	0.04–0.11	0.03–0.12	0.30–0.60
Bright	0.09–0.25	0.10–0.27	0.12–0.29	0.13–0.30	0.16–0.38	0.19–0.48	0.20–0.52	0.23–0.61

Table 3

The reflectance ranges of dark and bright surfaces in the MZ spectral channels.

SurfaceType	Band
	Band 1	Band 2	Band 3	Band 4
	Effective Wavelength, nm
	490.5	548	670	841
Dark	0.02–0.06	0.04–0.09	0.03–0.09	0.18–0.56
Bright	0.12–0.29	0.16–0.38	0.20–0.53	0.22–0.60

Table 4

Values of AOT at $\lambda =550\mathrm{nm}$, retrieved for models of various underlying surfaces under the assumption $r_{surf}\left(412\mathrm{nm}\right)=0.028$, and RMSE, with three values of the real ${\tau }_{aer}\left(550\mathrm{nm}\right)$ and solar zenith angles 60° (upper numbers) and 20° (lower numbers).

${\tau }_{aer}\left(550\mathrm{nm}\right)$	Model of Underlying Surface							RMSE	BARRMSE
	GV	Oat	Grass	Rye	Beet	Lawn	Soil
0.1	0.050	0.069	0.087	0.096	0.105	0.179	0.179	0.048	0.07
	0.050	0.050	0.069	0.087	0.105	0.233	0.243	0.080
0.3	0.215	0.270	0.288	0.298	0.307	0.371	0.371	0.051	0.11
	0.151	0.252	0.270	0.288	0.307	0.426	0.426	0.090
0.5	0.426	0.472	0.481	0.499	0.499	0.500	0.500	0.031	0.15
	0.371	0.453	0.472	0.490	0.500	0.500	0.500	0.053

Table 5

RMSE of reflectance retrieval in the MS spectral channels for dark (upper numbers) and bright (lower numbers) surfaces.

Solar Zenith Angle	Wavelength, nm
	412	442	489	509	559	619	664	776
60°	0.0089	0.0079	0.0065	0.0059	0.0042	0.0036	0.0033	0.0125
	0.0249	0.0186	0.0134	0.0133	0.0191	0.0234	0.0237	0.0235
20°	0.0086	0.0074	0.0058	0.0052	0.0038	0.0032	0.0028	0.0101
	0.0126	0.0097	0.0098	0.0089	0.0113	0.0123	0.0118	0.0117

Table 6

RMSE of reflectance retrieval for dark (upper numbers) and bright (lower numbers) surfaces in the MZS HR channels.

Solar Zenith Angle	Band
	Band 1	Band 2	Band 3	Band 4
60°	0.0064	0.0045	0.0028	0.0034
	0.0137	0.0197	0.0295	0.0353
20°	0.0059	0.0040	0.0025	0.0033
	0.0092	0.0128	0.0163	0.0190

Table 7

RMSE of reflectance retrieval for dark surfaces (upper numbers) and bright surfaces (lower numbers) in the MZ spectral bands with an error 5% of radiometric calibration in the MS channel 412 nm.

Solar Zenith Angle	Band
	Band 1	Band 2	Band 3	Band 4
60°	0.0132	0.0090	0.0054	0.0091
	0.0146	0.0198	0.0295	0.0354
20°	0.0087	0.0056	0.0033	0.0062
	0.0092	0.0128	0.0163	0.0190

Table 8

RMSE of reflectance retrieval for dark surfaces (upper numbers) and bright surfaces (lower numbers) in the MZ spectral channels with an error 5% of radiometric calibration in these channels.

Solar Zenith Angle	Band
	Band 1	Band 2	Band 3	Band 4
60°	0.0147	0.0114	0.0078	0.0192
	0.0274	0.0363	0.0471	0.0543
20°	0.0106	0.0086	0.0063	0.0186
	0.0230	0.0278	0.0329	0.0373

Table 9

Optical and microphysical parameters of aerosol models.

Parameter	“Continental”	“Belarus”	“Maritime”	“Oceanic”
$r_{eff},\mathsf{\mu }\mathrm{m}$	-	0.330	1.137	2.490
SSA	0.890	0.870	0.986	1.000
$x\left(1200\right)$	0.183	0.146	0.096	0.069
$\alpha$	1.116	1.376	0.238	−0.098

Table 10

RMSE of AOT retrieval with various models of underlying dark surfaces in cases with aerosol models “Belarus” and “Oceanic” considered as true «real» atmospheric aerosols.

Aerosol Model	Solar Zenith Angle	AOT (550 nm)
		0.1	0.3	0.5
“Belarus”	60°	0.044	0.091	0.147
	20°	0.066	0.144	0.217
“Oceanic”	60°	0.044	0.163	0.258
	20°	0.079	0.126	0.174

Table 11

RMSE of reflectance retrieval for dark surfaces (upper numbers) and bright surfaces (lower numbers) in the MZ spectral channels in the cases of “Belarus” and “Oceanic” models as a “real” atmospheric aerosol.

Aerosol Model	Solar Zenith Angle	Band
		Band 1	Band 2	Band 3	Band 4
“Belarus”	60°	0.0072	0.0054	0.0042	0.0119
		0.0222	0.0212	0.0280	0.0353
	20°	0.0066	0.0048	0.0037	0.0110
		0.0149	0.0138	0.0152	0.0186
“Oceanic”	60°	0.0074	0.0082	0.0104	0.0135
		0.0277	0.0407	0.0535	0.0614
	20°	0.0103	0.0148	0.0181	0.0343
		0.0323	0.0405	0.0464	0.0525

Table 12

Resulting RMSE of surface reflectance retrieval in the MZ channels and the contribution of various sources of error. Case with a dark underlying surface.

Source of Error	Solar Zenith Angle	Band
		Band 1	Band 2	Band 3	Band 4
Resulting RMSE	60°	0.0186	0.0144	0.0111	0.0249
	20°	0.0137	0.0127	0.0120	0.0301
Method error (Basic version)	60°	0.0064	0.0045	0.0028	0.0034
	20°	0.0059	0.0040	0.0025	0.0033
MS calibration error (5%)	60°	0.0117	0.0080	0.0055	0.0090
	20°	0.0073	0.0047	0.0035	0.0062
MZ calibration error (5%)	60°	0.0124	0.0101	0.0070	0.0194
	20°	0.0085	0.0074	0.0056	0.0190
Error of atmospheric aerosol model	60°	0.0036	0.0045	0.0059	0.0123
	20°	0.0054	0.0082	0.0097	0.0223
RMSE in the absence of a calibration error	60°	0.0073	0.0064	0.0065	0.0127
	20°	0.0080	0.0092	0.0100	0.0225

Table 13

Resulting RMSE of surface reflectance retrieval in the MZ channels and the contribution of various sources of error. Case with a bright underlying surface.

Source of Error	Solar Zenith Angle	Band
		Band 1	Band 2	Band 3	Band 4
Resulting RMSE	60°	0.0329	0.0361	0.0425	0.0461
	20°	0.0289	0.0311	0.0335	0.0352
Method error (Basic version)	60°	0.0137	0.0197	0.0295	0.0353
	20°	0.0092	0.0128	0.0163	0.0190
MS calibration error (5%)	60°	0.0016	0.0007	0.0003	0.0011
	20°	0.0000	0.0000	0.0000	0.0000
MZ calibration error (5%)	60°	0.0227	0.0212	0.0211	0.0221
	20°	0.0176	0.0176	0.0188	0.0204
Error of atmospheric aerosol model	60°	0.0195	0.0216	0.0221	0.0199
	20°	0.0210	0.0223	0.0224	0.0215
RMSE in the absence of a calibration error	60°	0.0238	0.0292	0.0369	0.0405
	20°	0.0229	0.0257	0.0277	0.0287

Table 14

Comparison of the accuracy of surface reflectance retrieval by different algorithms.

Surface	FLAASH	DG AComp	RACE
	RMSE	RMSE	RMSE in the Absence of a Calibration Error	RMSE with a Calibration Error of 5%
Dark surfaces	0.02–0.04	0.01–0.02	0.01	0.01–0.02
Bright surfaces	0.14–0.19	0.05	0.02–0.04	0.03–0.04

Table 15

Data on test sites and experimental conditions used for testing the RACE algorithm.

	La Crau,France	Railroad Valley, NV, USA	Gobabeb,Namibia
Data	16.05.2018	17.07.2018	16.05.2018
Site	LCFR01	RVUS00	GONA01
Latitude	43.56	38.497	−23.60
Longitude	4.86	−115.69	15.119220
Altitude, m	20	1435	510
Solar zenith angle, deg.	24.87	17.56	42.74
AtmosphereAOT (550 nm)	0.119	0.138	0.036
Angstrom parameter	1.629	1.273	1.045
The water vapor content in the atmosphere, g·cm−2	1.82	1.884	1.95

Table 16

The results of retrieving spectral reflectance of surfaces with the RACE algorithm in comparison with the data of reflectance ground measurements.

Wave-Length, nm	La Crau, France			Railroad Valley, NV, USA			Gobabeb,Namibia
	Measured	Retrieved	Error	Measured	Retrieved	Error	Measured	Retrieved	Error
	AOT
0.55	0.119	0.208	0.089	0.360	0.539	0.179	0.138	0.539	0.401
	Surface reflectance
0.410	0.0375	0.0285	−0.0090	0.2254	0.2384	0.0130	0.0999	0.0540	−0.0459
0.440	0.0515	0.0446	−0.0069	0.2576	0.2784	0.0208	0.1235	0.0946	−0.0289
0.490	0.0668	0.0626	−0.0042	0.3001	0.3271	0.0270	0.1564	0.1469	−0.0095
0.510	0.0764	0.0734	−0.0030	0.3158	0.3457	0.0299	0.1730	0.1709	−0.0021
0.560	0.1123	0.1093	−0.0030	0.3680	0.4004	0.0324	0.2346	0.2463	0.0117
0.620	0.1243	0.1220	−0.0023	0.3913	0.4201	0.0288	0.3053	0.3279	0.0226
0.670	0.1221	0.1227	0.0006	0.3978	0.4230	0.0252	0.3311	0.3591	0.0280
0.780	0.2629	0.2526	−0.0103	0.4308	0.4351	0.0043	0.3780	0.3910	0.0130
0.870	0.2832	0.2835	0.0003	0.4281	0.4439	0.0158	0.3679	0.3930	0.0251