We consider an algorithm of atmospheric correction of satellite images for reconstruction the reflection coefficients of the ground surface in the visible and near-IR range. The algorithm accounts for the effect of inhomogeneity of the ground surface, adjacency effect, additional irradiance by reflected radiation, and polarization of radiation. The capabilities of this algorithm and MOD09 algorithm are compared using areas in Tomsk, Moscow, and Irkutsk regions of the Russian Federation as examples.
the authors of the works [12-14] for solving the direct problem of the solar radiative transfer in the atmospheresurface system. Their estimates indicate that polarization may introduce an error up to 10% into the intensity of received radiation. In our works [15-16] it is shown that, in reconstruction the reflection coefficients of certain surfaces, the neglect of polarization may lead to absolute errors exceeding the value of reconstructed reflection coefficient. Therefore, the effect of polarization of radiation should be taken into consideration for areas covered, e.g., by weakly reflecting vegetation. Below, we will consider a modified algorithm of reconstruction the reflection coefficients, taking into account the polarization of radiation, and will verify it against MODIS images, as an example.
The algorithm of correction was developed assuming that: the atmosphere is spherical and is divided into 32 homogeneous layers; the atmosphere is a scattering and absorbing medium; the ground surface is Lambertian; relief is disregarded; and the atmosphere partially polarizes radiation upon scattering.
This algorithm had been a basis for a software package, the block-diagram of which is presented in Figure 1.
The step-by-step procedure of atmospheric correction is as follows. 1) The direct transmission coefficient Ti of the path from the viewed pixel to receiving system is determined. 2) The Monte Carlo method is used to calculate the intensity of radiation Isun, not having interacted with the ground surface, taking into account the polarization of radiation for 30 nodal directions. Based on these results, the approximate Isun,I values are determined from the formula:
Ai = {
C11μ2d,i + C21 (√1 − μ2d,icosφi) + C22μd,i√1 − μ2d,icosφi − (√1 − μ2d,isinφi) , φi ≤ 900 C11μ2d,i + C31 (√1 − μ2d,icosφi) + C32μd,i√1 − μ2d,icosφi − (√1 − μ2d,isinφi) , φi > 900 2 2 Isun,i = −
Bi+√Bi2−4AiC13
2Aiμd,i 2 2 Bi =
C12μd,i + C23√1 − μ2d,icosφi , φi ≤ 900 { C12μd,i + C33√1 − μ2d,icosφi , φi > 900 (1) (2) (3) where , is the cosine of the viewing zenith angle in observation of the ith pixel; is the azimuth between the directions toward Sun and toward receiver in the ith pixel; and 11, 12, 13, 21, 22, 23, 31, 32, and 33 are the approximation constants, determined by the least-squares method (LSM) from nodal values for a fixed .
Algorithm of the Isun calculation with accounting for polarization and its testing were described in [15,16]. 3) The Monte Carlo method is used to calculate the nodal values of the integral of the point spread function (PSF) of the channel of formation of the adjacency effect . These results are used to determine the boundaries of isoplanar zones θ1,k, (regions on the ground surface within which the same PSF can be used with a specified error ) according to the criterion: {
I
(μk) = I μk+1 = 1 − [1 (I
A (1) − A(1 − μk)N (1) − Isurf(μk))]1⁄N
1+δ δ ≡
I (μk)−I
(μk+1)
Isurf(μk+1) I (μ) ≡ ( ) + ( ) = ( ) + ∫ h(μ, rw, φw)d
S used in the calculations); and A, N are approximation constants determined using LSM.
where =
1, specifies the boundary between the kth and k+1st isoplanar zones; ( ) is the direct transmission on the path from a point on the ground surface to receiving system for the cosine of angle of deviation from direction of nadir ; ℎ( , , ) is PSF of the channel of formation of adjacency effect; r w is the surface distance from the center of the viewed pixel on the ground surface to a point on the ground surface; φw is the azimuth angle on the ground surface between the direction toward the projection of receiving system onto the ground surface and the direction toward a given point away from the pixel viewed; S is the entire area of the ground surface; δ is the maximum admissible error level in using PSF corresponding to μk, instead of PSF corresponding to μk+1 (δ=0.05 was 4) The radius of the region of adjacency effect Rk is determined, outside of which the adjacency effect can be considered to be zero with a specified error δ1. It can be shown that, in order for the condition:
1 ≥ mi n ̃ ≥ 1 adjacency effect (δ1 = 0.95 inourcalculation)s, fulfillment of the condition is sufficient:
to be satisfied for an arbitrary inhomogeneous surface, where Qi ≡ rsurf,iEsum,i is a certain exact value of the surface emissivity in the ith pixel; Q̃i is an approximate value of the emissivity, obtained assuming that no adjacency effect exists outside Rk; and δ
1 is a quantity, characterizing the maximal error due to the use of the radius of the f1(Rk) ≡ ∬S(Rk) h(μk,rw,φw)ds ∬S h(μk,rw,φw)ds ≥ δ1 + (δ1 − 1) 1 πTk ∬S h(μk,rw,φw)ds Isum,i = Isun,i + ∑jN=i1 Ai,jQj + Aout,i̅Q̅̅ii = ̅1̅,̅̅N̅
1 Ai,j ≈ {π Ti + h(μki, 0,0)Si , i = j
h(μki, rw,j, φw,j)Sj , i ≠ j
Aout,i = Idif(μki) − ∑jN=i1 Ai,j
̅Q̅̅i = 1Isum,i−Isun,i
πTi+Idif(μd,i)
(
where Rk is the radius of the adjacency effect for the kth PSF, calculated for the angles (
For MODIS channels considered below, for 0.1≤AOD0.55≤5, and for an arbitrary unknown distribution of
reflection coefficients over the ground surface and different situations, in [6] upper estimates of Rk were obtained, for
which the condition (
For each of k isoplanar zones, the Monte Carlo method is used to calculate PSF of the channel of formation of 5) The Seidel method is used to solve the system of linear algebraic equations (SLAE) for determining the where Isum,i is the intensity of total radiation received by satellite system; N is the number of pixels in the area under consideration; Ni is the number of pixels within the radius Rk around the ith pixel; μki is the boundary of the kth isoplanar zone into which the ith pixel falls; ̅Q̅̅i is the quantity that estimates approximately the surface emissivity outside the region under consideration; and Si is the area of the ith pixel.
Solving the system of equations (
6) The Monte Carlo method is used to calculate the irradiance of the ground surface without accounting for rereflections E0.
7) Radius is calculated for the region of formation of additional irradiance by singly reflected radiation R, outside of which the additional irradiance can be assumed to be zero. It can be shown that, in order for the condition: 1 ≥ mi n
, ≥ 1
̃ ,
to be satisfied for an arbitrary inhomogeneous surface, where rsurf,i is a certain reflection coefficient of the ith
pixel; r̃surf,i is an approximate value of the reflection coefficient, obtained in using the radius of the region of
formation of additional irradiance; δ2 is the quantity characterizing the maximal error due to the use of the radius of
formation of additional irradiance (δ2 = 0.95 inourcalculation)s, fulfillment of the condition is sufficient:
(
γ1 γ1 = ∬S h1(rw)ds Qi = rsurf,i (1 + ∑jM=i1 Ci,jrsurf,j + Cout,i̅r̅s̅u̅̅r̅f̅,i + E0 (̅r̅s̅̅u̅̅r̅f̅,̅iγ1)2) 1−̅r̅s̅̅u̅̅r̅f̅,̅iγ1 Ci,j ≈ h1(rw,ij)Sj Cout,i = γ1 − ∑jM=i1 Ci,j
̅Q̅̅i⁄E0 ̅r̅s̅u̅̅r̅f̅,i = 1+γ1̅Q̅̅i⁄E0 where h1(rw) is the value of PSF of channel of formation of additional irradiance.
For MODIS channels considered below, for 0.1≤AOD0.55≤5, and for different situations, in [6] we estimated R
values, for which the condition (
The Monte Carlo method is used to calculate PSF of channel where additional irradiance of the ground surface is formed by radiation reflected in the atmosphere-surface system h1(rw) within the radius R.
8) The Newton method, with auxiliary SLAE solved iteratively by the Seidel method, is used to solve the nonlinear system of equations in rsurf of the form:
where rsurf,i is the reflection coefficient of the ith pixel of the image; ̅r̅s̅u̅̅r̅f̅,i is the reflection coefficient of the ith pixel, obtained in the approximation of a homogeneous ground surface; and Mi is the number of pixels within the radius R around the ith pixel.
Solving the nonlinear system (
Analysis of convergence conditions for the systems of equations (
To test the performance and to estimate the error of the algorithm of reconstruction the reflection coefficients with accounting for the polarization effect, we considered MODIS images for 5 channels: channel 1 centered at λ=0.649 m, channel 2 at λ=0.860 m, channel 3 at λ=0.469 m, channel 4 at λ=0.555 m, and channel 8 at λ=0.412 m, with the spatial resolution of 1000 m. We considered three test regions 1) area in the south of Tomsk region (55.950 – 56.850N and 84.050 – 84.950E), 7 images from June 17, 2012 to June 23, 2012, 2) area in Moscow region (55.720-55.950N and 37.560-38.100E), 5 images from May 6, 2017 to May 7, 2017, and 3) area in Irkutsk region (51.420 -52.670N and 103.640 -105.470E), 4 images from June 20, 2017 to June 21, 2017. Results from our algorithm with and without accounting for polarization were compared with those from MOD09 algorithm and those, obtained without atmospheric correction.
To estimate the errors of the algorithms, we considered test points at the centers of coniferous forest massifs. The errors of the algorithms were estimated as the difference between the reconstructed coefficients and measurements [17], presented in Table 1.
As analysis showed, the reconstructed reflection coefficients for points at the centers of coniferous forest massifs in MODIS channels centered at λ=0.412, 0.469, 0.555, and 0.649 m differ little from ground-based measurements presented in work [17] when aerosol content is small (AOD0.55≤0.1). Therefore, for these channels and areas, the data from [17] can be used as reference values; and the differences from the data in [17] estimate approximately the errors of the algorithms. The reflection coefficients in the MODIS channel centered at λ= 0.860 m differ from measurements in [17] markedly stronger because the reflection in this channel depends appreciably on the state (productivity) of vegetation. Therefore, data in [17] cannot be used as a reference for this channel in a number of situations.
AERONET data [18] on aerosol optical depth (AOD), particle size distribution, and complex refractive index were used as initial data for specifying the atmospheric model. Profiles of temperature and pressure from MODIS measurements [19] were additionally used. The AOD0.55 values were in the range from 0.1 to 1.52 for the area in Tomsk region, from 0.04 to 0.07 for the area in Moscow region, and from 0.04 to 0.06 for the area in Irkutsk region. An example of reconstructed distributions of reflection coefficients over the ground surface is presented in Figure 2. Figure 3 presents an example of comparison of reflection coefficients of the ground surface, obtained using MOD09 algorithm and our algorithm with accounting for the polarization. From Figure 3 it can be seen that the results from the algorithms well agree for images with low atmospheric turbidity. 0,08 f,renw0,06 rsu 0,04 0,02 0,00 0,00 0,06 rsurf,MOD09 0,02 0,04 0,08 0,10 0,12
The errors of these algorithms were estimated for 3 points (one point for each area): 1) a point on the territory of Tomsk State Nature Reserve (56.20N, 84.30E); 2) a point in Losiny Ostrov National Park (58.850N, 37.830E); and 3) a point in Krasny Yar State Nature Reserve (52.520N, 105.060E).
For each channel and each point, we determined the difference from reference value Δrsurf, averaged over these images. The average Δrsurf values thus obtained are presented in Table 2.
Comparison of Δrsurf of our algorithm with Δrsurf of MOD09 algorithm for the point in Tomsk region shows that our algorithm with accounting for polarization for these test images shows markedly less differences from reference values for channels centered at λ=0.649, 0.469, and 0.555 m and almost identical differences at λ=0.860 m as compared to MOD09 NASA algorithm. However, the MOD09 values are closer to data in [17] at the wavelength of 0.412 m.
For the test point in Moscow region, the reflection coefficients of the ground surface, reconstructed using algorithms considered here (except algorithm without atmospheric correction), deviate from the reference values by almost the same amount. The results for this point at λ=0.860 m differ much stronger from data in [17] than those for the first point. This is probably because results in [17] strongly diverge from actual reflection coefficient in this channel and cannot be used as a reference for the situation, considered here.
Comparison of results for the test point in Irkutsk region shows that our algorithm with accounting for polarization for these images shows somewhat smaller differences from reference values at wavelengths 0.469, 0.555, and 0.412 m than the MOD09 algorithm and almost identical differences at λ=0.860 and 0.649 m. 4
Comparison with MOD09 algorithm shows that our algorithm gives much smaller rsurf reconstruction errors than the MOD09 NASA algorithm at λ= 0.469, 0.555, and 0.649 m, and gives an error of the same order of magnitude at λ= 0.860 m. At λ= 0.412 m, MOD09 algorithm reconstructs rsurf with a smaller error in some cases, our algorithm is more preferable in some other cases, and algorithms give errors of the same order of magnitude in the other cases.