Fast Computations of the Top-of-the-Atmosphere Radiance in a Spectral Range 400–2500 nm Using the PYDOME Tool †

: Accurate computations of the radiance in the gaseous absorption bands typically require fine wavelength steps. In this paper, a fast technique for computing a radiance spectrum in the wavelength region of 400–2500 nm is proposed. Our approach draws inspiration from the established k -correlation distribution model, with k denoting the absorption coefficient. However, our method expands upon this by considering both the direct transmittance and the scattering coefficient as predictors. At selected spectral points, the full radiative transfer simulations are performed, and the mathematical relation between a predictor and the radiance is established. Then, the radiance is restored on a fine wavelength grid. This approach can be used to enhance the accuracy of the convolved spectrum computations based on precomputed approximately monochromatic lookup tables and reduce the size of lookup tables. Numerical analysis demonstrates the method’s applicability to scenarios characterized by aerosol optical thicknesses not exceeding two.


Introduction
Precise computations of the atmospheric radiative processes are crucial for analyzing the remote sensing data obtained by hyperspectral and multispectral air-and space-borne spectrometers.Accurate simulations of the solar radiation spectrum absorbed and scattered by the terrestrial atmosphere require a fine wavelength grid due to the many absorption features in the spectrum.For instance, simulations of narrow absorption bands (like the O 2 A band between 758 and 770 nm) may require as many as 10,000 spectral lines to be considered.Even the simulations of the multi-spectral measurements require, first, an accurate fine wavelength resolution spectrum that will be convolved with a more or less sampled channel spectral response function (SRF).In this regard, fast radiative transfer models (RTMs) are required, which are specifically designed to perform hyperspectral (with the spectral sampling distance (SSD) between 0.1 and 0.4 nm) computations in an efficient manner.
Typically, such techniques involve two steps.In the first step, radiative transfer computations are performed at selected spectral points, while in the second step, some approximations are used to guess the radiance values at 'skipped' wavelength points using those computed in the previous step.The most common approach to accelerate hyperspectral computations is the k-distribution method (where k stands for the absorption coefficient), which is due to Ambartzumian [1].It is based on the assumption that the transmission within a spectral interval is not dependent on the line-by-line (LBL) variation of k with respect to wavelength λ but rather on the distribution of absorption coefficient within the spectral interval.In this regard, the cumulative frequency distribution of k, namely, G(k), is introduced (referred to as the k-distribution function); the convolution of the spectrum with the instrument slit function is performed not in the λ space but rather in the G space.Unlike k(λ), the function G(k) is a smooth function, and so fewer quadrature points for numerical integration are required.This method was extended by Goody et al. [2] to the cases of inhomogeneous atmospheres, assuming that there is a correlation between k-distributions at different pressure levels.The performance enhancement due to the k-distribution method is about an order of magnitude [3,4].
However, the k-distribution methods have some drawbacks.First, the radiance depends not only on the absorption coefficient but also on the scattering coefficient and hence, there is no one-to-one correspondence between the k-value and the radiance value (however, it is implicitly assumed in the k-distribution method).Consequently, the k-distribution method may be inaccurate, as the spectral range increases and is usually applied within a given absorption band.Secondly, due to the multiple scattering impact and screening effect of a turbid medium, the impact of the absorption processes can be reduced, and the radiance dependency on k may be not purely exponential, which leads again to errors in the computations.In fact, as discussed in [5], the k-distribution method is valid only in the context of a single absorption line, periodic lines, and the strong-and weak-line limits.Finally, the result of the k-distribution method is the convolved spectrum on a coarse wavelength grid, rather than the spectrum on an initial fine grid.
In this paper, we adapt the idea of the k-distribution method for computations of lookup tables of the reflected radiation in the wide spectral range.Our goal is to design a method that allows to simulate a radiance spectrum in a reasonable time (e.g., several hours) suitable for lookup table calculations.

Methodology 2.1. Radiative Transfer Model PYDOME
Our simulations are conducted using the PYDOME radiative transfer model.PYDOME is a Python-based model that employs the discrete ordinate method with matrix exponential techniques [6].The left eigenvector approach [7] and neural network-based eigenvalue solver [8] implemented in PYDOME enhance the efficiency of the eigenvalue decomposition in the discretized radiative transfer equation compared to traditional methods.
The molecular scattering is modeled following [9].PYDOME has an interface to the Py4Cats model [10] for absorption coefficient calculations and a database of measured absorption cross sections [11].It accommodates various aerosol models, including OPAC and MODTRAN definitions.PYDOME features a lookup table mode, optimizing computations through parallel processing.In this mode, it generates lookup tables with spectra and atmospheric functions if necessary on a fine resolution grid, enabling subsequent convolution with specific sensor SRF.By default, eight discrete ordinates per hemisphere are used.The viewing zenith angles (VZAs) coincide with the discrete ordinates.In addition, the false discrete ordinate method [7] is used to compute radiances at given VZAs.
For this study, the computations are performed in spectral ranges 400-1000 nm and 1000-2500 nm.To compute a reference spectrum, for each of the two subdomains, N = 50, 000 spectral points are selected equidistantly in the wavelength domain.The simulations are performed for a Lambertian surface albedo of 0.1.The gaseous profiles correspond to the midlatitude summer atmosphere.We consider a clear atmosphere (no aerosol, no clouds) and the rural aerosol [12] with the optical thicknesses 0.3, 0.85 and 2.0 at 550 nm.In addition, one cloud scenario is considered with an optical thickness 10, a cloud geometrical thickness of 1 km, and a cloud top height of 5 km.

Formulation of the Modified k-Distribution Method
The proposed method can be formulated as follows.For a set of N spectral points with the wavelengths λ i , i = 1, . . ., N, and set of layers j = 1, . . ., N l , we compute scattering σ ij and absorption k ij coefficients.The radiances R m = R(λ ms ) are computed at each sth spectral point, where m = si (hereinafter, s is referred to as the sampling rate) at M spectral points in total (M < N).The transmittances T σi and T ki are computed as where h j is the geometrical thickness of the jth layer, while µ 0 and µ are the cosines of the solar zenith angle and viewing zenith angle, respectively.Then, a predictor vector p i can be taken as Note that the third term 1/λ 4 i corresponding to Rayleigh scattering is added empirically to improve the results in the spectral range 400-600 nm.Assuming a linear model between the components of p i and the radiance value R i = R(λ i ), we have: Then, we choose a spectral range containing K consequent spectral points from R m array, for which the radiances are known, and make an over-determined system of Equation ( 5).The unknown coefficients a, b, c and d are solved by least squares minimization.Finally, the values of radiances within a considered spectral range are computed using the found coefficients in Equation (5).As shown in Figure 1, in red points, the full radiative transfer computations are performed, while the predictor computations based on optical parameters are evaluated both at the white and red spectral points.After the coefficients a, b, c and d are found from the 'red' points, the radiances are restored at the white points.

Note that by setting
we derive an analog of the k-distribution model, while by setting with Ri , the radiance is found by using the two-stream RTM, and the low stream regression model [13] is obtained.

Results of Simulations
The described method was applied to test scenarios.The restored spectra and the reference spectra are convolved with the Gaussian SRF.The relative error is computed across the spectrum and averaged over geometries within each scenario type.The values of the sampling rates are listed in Table 1. Figure 2 shows the examples of spectral computations for the clear and aerosol and cloud cases.For the clear sky and aerosol cases, SRF FWHM is 0.04 nm and sampling rate is s = 2048, while for the cloud scenario, the sampling rate is s = 64.The errors of the method for the spectral range 400-1000 nm are summarized in Table 1.The error increases with the sampling rate.However, even for s = 2048, the mean error is below 1%.The maximum errors correspond to some absorption lines in the water vapor region around 950 nm.For the cloud case, the mean error is about 1% at the sampling rate s = 128, while the maximum error is larger than 15%.Also, note that the results for the cloud scenario are somewhat unstable since the maximum error for s = 128 is less than that for s = 16.
The results for the spectral range 1000-2500 nm are summarized in Table 1.Overall, the errors are lower than those for the 400-1000 nm range.It can be explained by the fact that the Rayleigh scattering cross section is negligibly small in the infrared region and so, the impact of the scattering cross section and multiple scattering is small.The maximum error (about 50%) is at 1400 nm, where the absolute value of the radiance is close to zero.For clear and aerosol cases, the sampling rate s = 2048 keeps the mean error around 0.5%.
For the cloud scenario, the sampling rate below s = 128 keeps the mean error below at least 1%.
A practical conclusion from our analysis is that as long as the aerosol loading is relatively low (below 2), the full radiative transfer simulations can be performed on a coarse wavelength grid with a step of about 20 nm in the 400-1000 nm range and 60 nm in 1000-2500 nm.In this case, using the described method, a high-resolution spectrum can be restored with an error than 0.5%.The simulations in the case of optically thick clouds require a finer discretization step.In fact, for this case, the principal component-based methods [14,15] look more promising.Setting s = 2048, the computation time required for simulation of the spectrum in the 400-2500 spectral range is about 20 seconds on a personal computer with 16 GB RAM and processor Intel Core i7-6700 CPU 3.40 GHz.Finally, we test our approach on the Sentinel 2 bands.The reference spectra are computed using a fine resolution spectrum with 0.04 nm step convolved with the Sentinel-2 SRFs.They are compared against the spectrum computed with 0.4 nm step convolved with the SRF and spectra computed with 0.4 nm step and restored at 0.04 nm using the described method.The results for four different geometries are shown in Table 2 .Overall, we can conclude that the proposed method is accurate for small aerosol optical thicknesses.Band 9 corresponding to the water vapor band has the largest error, which is below 0.5%.

Conclusions
We designed a method for fast simulations of the hyperspectral radiances in the spectral range of 400-2500 nm.The method is based on computing a predictor from optical parameters and relating it to the radiance value at a given wavelength point.Even though the same idea is exploited in the k-distribution method, there are several major differences: First, unlike the k-distribution method, the presented technique computes the high-resolution spectrum, which is convolved afterwards.Secondly, the developed method can be applied to the wide spectral range containing several absorption bands as well as regions with a strong influence of Rayleigh scattering.The required least squares method does not introduce a performance bottleneck in the whole algorithm.Therefore, the obtained acceleration factor is very close to the theoretical maximum that would be equal to s.Also, it has been shown that the k-distribution method can be regarded as a particular case of the presented approach.However, the described method looks more convenient for lookup table computations.Indeed, the lookup tables store the approximately monochromatic radiances on a coarse grid and optical parameters on a fine grid.In this way, the high-resolution spectrum can be readily restored and convolved with an arbitrary slit function.

Figure 1 .
Figure 1.Illustration of the method.Here, the sampling rate is s = 5.

Figure 2 .
Figure 2. Radiance spectra for clear sky, aerosol and cloud scenarios convolved with SRF with FWHM = 0.04 nm and the sampling rate s = 2048 (s = 64 for the cloud scenario): blue -reference spectra, orange -spectra computed with the proposed method.

Table 1 .
The efficiency of the method when the SRF with FWHM = 50 nm is applied.