Two Machine Learning Based Schemes for Solving Direct and Inverse Problems of Radiative Transfer Theory Dmitry S. Efremenko[0000−0002−7449−5072] , Himani Jain and Jian Xu[0000−0003−2348−125X] Remote Sensing Technology Institute, German Aerospace Center (DLR), 82234 Oberpfaffenhofen, Germany {dmitry.efremenko, himani.jain, jian.xu}@dlr.de Abstract. Artificial neural networks (ANNs) are used to substitute computation- ally expensive radiative transfer models (RTMs) and inverse operators (IO) for retrieving optical parameters of the medium. However, the direct parametriza- tion of RTMs and IOs by means of ANNs has certain drawbacks, such as loss of generality, computations of huge training datasets, robustness issues etc. This paper provides an analysis of different ANN-related methods, based on our re- sults and those published by other authors. In particular, two techniques are pro- posed. In the first method, the ANN substitutes the eigenvalue solver in the dis- crete ordinate RTM, thereby reducing the computational time. Unlike classical RTM parametrization schemes based on ANN, in this method the resulting ANN can be used for arbitrary geometry and layer optical thicknesses. In the second method, the IO is trained by using the real measurements (preprocessed Level-2 TROPOMI data) to improve the stability of the inverse operator. This method pro- vides robust results even without applying the Tikhonov regularization method. Keywords: Radiative Transfer, Machine Learning, Trace gas retrieval 1 Introduction Machine learning techniques have become of paramount importance in the geosciences, lighting engineering and remote sensing [22]. Atmospheric parameter retrieval (in par- ticular, trace gas retrieval) is an important application of remote sensing [20]. At- mospheric composition sensors on board satellites provide a huge amount of data of high spatial resolution about atmospheric constituents [30]. The sensor measures the spectral radiances reflected by the terrestrial atmosphere. From these measurements, information can be retrieved about ozone (including tropospheric ozone), nitrogen, formaldehyde, carbon monoxide as well as clouds and aerosols. Radiative transfer mod- els (RTMs) are key components of the algorithms designed for the retrieval of atmo- spheric constituents from remote sensing data. The RTMs encompass our understand- ing of the physics behind the measurement process and relate the optical parameters Copyright © 2020 for this paper by its authors. Use permitted under Creative Commons Li- cense Attribution 4.0 International (CC BY 4.0). 2 D.S. Efremenko et al. of the medium with the signal measured by the sensor. To process a big amount of data comping from the state-of-the-art atmospheric composition sensors, acceleration techniques for RTMs are required. Numerous approaches based on analytical proper- ties of dimensionality reduction techniques have been developed (see e.g. [3, 14, 26] for general review). To boost further the performance of RTMs, the machine learning based approaches are widely used. In particular, artificial neural networks (ANNs) are used to reproduce the results of RTMs, but at low computational costs. As input, ANN takes the optical parameters of the atmosphere, while as output it provides the radiance field. Despite its high efficiency, this approach has a certain drawbacks. For instance, extensive computa- tions of the training datasets are required. Usually, the training is restricted to a chosen model for the single scattering phase function (e.g. Henyey-Greenstein), so no arbitrary phase functions are allowed. Besides, ANN does not necessary capture implicit proper- ties of the radiance field, such as the symmetry relation with respect to the incident and viewing polar angles. Finally, the use of ANNs for the forward RTM parameterization relies on the universal approximation theorem [9, 19] which states that a feed-forward network with a single hidden layer and with a finite number of nodes can approximate any continuous function from a compact interval. However in practice, the ANN-based approach is not fail-save, i.e. the error can be large in some cases which have not been captured by the training procedure. ANNs are used to parameterize an inverse operator in retrieval problems. The ad- vantage of this approach is that the explicit inversion of the RTM is avoided, while the resulting ANN operator is fast. However, this approach requires two regularization pro- cedures. One regularization is required to deal with ill-posedness of the inverse problem and is usually done by means of the Tikhonov regularization in the L2 metric [29]. The second regularization is required to obtain robust weights of the ANN itself and can be done in the L1 metric. In this paper we consider two alternative approaches. In the first method, the ANN is used to substitute a part of RTM, namely, the eigenvalue solver, thereby reducing the computational time. In the second method, the ANN inverse operator is trained on real measurements (namely, preprocessed TROPOMI data) to improve the robustness of retrievals. 2 Overview of the Radiative Transfer Model We consider computations of the spectral radiances in the Huggins band (280-335 nm), which is used for ozone profile retrievals. Due to strong influence of Rayleigh scatter- ing, multiple scattering plays an important role. The one-dimensional radiative transfer equation for the radiance L reads as follows [8]: Z2π Z1 dL (τ, µ, ϕ) ω (τ ) µ = −L (τ, µ, ϕ) + p (τ, µ0 , µ, ϕ − ϕ0 ) L (τ, µ0 , ϕ0 ) dµ0 dϕ0 , dτ 4π 0 −1 (1) where τ is the optical depth, µ and ϕ are the cosine of the polar angle and the azimuth angle, respectively, ω is the single scattering albedo, while p is the single scattering Two Machine Learning Based Schemes for Solving Direct and... 3 phase function. The numerical solution of Eq. (1) is described in numerous textbooks and papers (e.g. [2, 6, 28]). Here we outline main steps to put our considerations in the proper context. We consider the cosine azimuthal expansion of the radiance field and the phase function providing an equation for the azimuthal component Lm (τ, µ). Then, the radi- ance field is discretized in the µ -domain by considering Ndo Gaussian points {µi } per hemisphere and corresponding weights {wi }, where i = 1, ..., Ndo (the positive values of µ correspond to the downwelling radiance, while the negative values of µ correspond to the upward radiance). For an inhomogeneous atmosphere, we consider a spatial dis- cretization using N layers: τ1 < τ2 < ... < τN +1 , where τ1 = 0 and τN +1 = τs . A layer l is bounded above by the level τl and below by the level τl+1 . The single scatter- ing albedo and the single scattering phase function are assumed to be constant within each layer. For layer l with optical thickness τ̄l = τl+1 − τl Eq. (1) is rewritten in a discrete space as the following linear system of differential equations (for simplicity we omit azimuth index m): d i↑ (τ )   ↑  i (τ ) = −A l + bl (τ ) , τl ≤ τ ≤ τl+1 , (2) dτ i↓ (τ ) i↓ (τ ) where  ↑↓  i (τ ) i = Lm (τ, ∓µi ) , i = 1, ..., Ndo , (3) while A and b are the so called layer matrix and layer source vector, respectively (ex- pressions for them can be found in [12]). The integral of Eq. (2) is " # Z  ↑ ↑ τ̄l il Al τ̄l il+1 − ↓ +e ↓ = eAl t bl (t) dt. (4) il il+1 0 To evaluate eAl τ̄l , the eigenvalue decomposition method can be applied: eAl t = Vl eΛl t Vl−1 , (5) where Vl is the eigenvectors matrix and Λl is the eigenvalues matrix for the lth layer. Then, substituting Eq. (5) in Eq. (2) and rearranging the entries, we obtain an equation which relates the radiances at the layer boundary. It has the following form: " #   " ↓ #  ↑ i↑l R↑l T↑l il Pl = + , (6) i↓l+1 T↓l R↓l i↑l+1 P↓l where R↑↓ ↑↓ l and Tl are the reflection and transmission matrices, respectively, while ↑↓ Pl is the source term. They can be expressed in terms of Vl and Λl (ready-to-compute formulas are given in [17]). By using the matrix operator method [27], two layers can be merged into one layer. Applying this approach iteratively for each atmospheric layer, the equation similar to Eq. (6) is derived for the whole atmosphere. Finally, applying the boundary conditions, i↑0 corresponding to the upward (reflected) radiance at the top of the atmosphere can be found. 4 D.S. Efremenko et al. Table 1 shows the computational times for the eigenvalue decomposition part and the matrix operator method. The eigenvalue part consumes about 70 % of the total com- putational time. Also note, that the computational time of the whole discrete ordinate 2.5 RTM scales as ∼ Ndo due to eigenvalue problem and matrix multiplications in the framework of the matrix operator method. From Table 1 one can see that the eigenvalue problem (5) imposes a performance bottleneck in the whole discrete ordinate RTM. Several approaches were developed to accelerate or to avoid computations of the eigen- values. In particular, by setting Ndo = 1, we obtain the two-stream model, in which the eigenvalues are expressed analytically, and thus, Eq. (2) has a closed form solution. Table 1. Computational times of the eigenvalue decomposition part and the matrix operator method as a function of number of discrete ordinates. Computational time Number of discrete ordinates Eigenvalue solver Matrix operator method 2 0.53 0.2 4 1.78 0.58 8 5.68 1.92 16 26.0 10.3 32 100.1 49.7 3 Parameterization of the Eigenvalue Problem by Means of Artificial Neural Network The classical approach to parametrize the RTM is illustrated in Fig. 1 with a blue line. Here the RTM is used as a black box, and the ANN is capturing the input-output de- pendencies of RTM. In this case, the ANN is trained for a fixed number of layers. Besides varying the geometry (the viewing zenith angle, the relative azimuth angle and the incident angle), for each layer, the following parameters should be given: the optical thickness, the single scattering albedo, and the Legendre expansion coefficients of the phase function. Note that it is quite challenging to capture the dependencies of the radi- ance with respect to high-order expansion coefficients. This is why, usually, the phase function is restricted to a certain model. An alternative approach proposed in this paper is shown with a green line in Fig. 1 and consists in applying the ANN for computations of Λl . We recall again, that Λl depend neither on the layer optical thickness nor the geometry. Hence, as inputs, we require only the single scattering albedo and the phase function. Once Λl are known, Vl can be readily computed. We consider two types of layers. The first type of layers has only Rayleigh scattering and gaseous absorption. The phase function is smooth and assumed to be constant (the weak dependence on the depolarization factor related to the King factor is neglected [4]). The second type of layers has a cloud. The phase function in this case is computed by applying Mie theory [5] using the particle size parameter and the refractive index as Two Machine Learning Based Schemes for Solving Direct and... 5 Fig. 1. Scheme used for radiative transfer model parametrization by means of artificial neural networks. inputs. The particle size distribution is described by the modified gamma distribution [10]: γ f (r) = arα e−br , (7) where r is the particle size, while a, b, α and γ are real positive constants controlling the shape of the distribution. Note that f (r) satisfies the normalization condition, Z∞ f (r)dr = 1, (8) 0 and thus, f (r) is defined by three independent parameters. The examples of the cal- culated phase functions are shown in Fig. 2. The ANN for the Rayleigh layers uses as inputs the single scattering albedo only, while the ANN for the cloudy layers uses as in- puts the single scattering albedo and 3 parameters of the modified gamma distribution, i.e. in total 4 parameters. The outputs of the ANNs are the eigenvalues Λl . Note that we use Ndo = 16 and there are Ndo = 16 independent eigenvalues and eigenvectors for each azimuthal harmonic. The training is performed using Keras API. For the cloudy layers, the ANN config- uration with 1 hidden layer with 25 neurons is used. The training data set consists of 20000 samples, while the testing phase of the ANN training is performed on 5000 sam- ples. For Rayleigh scattering case, the ANN configuration with 1 hidden layer with 10 neurons is used. The training data set consists of 2000 samples, while the testing phase of the ANN training is performed on 500 samples. Note that by excluding dependen- cies on geometries and optical thickness, the dimension of the datasets is significantly reduced. The output functions (i.e.) depend smoothly on chosen inputs and the training procedure is fast (e.g. computations of the training dataset and the training procedure take less than two hours on Intel Core i5-3210). To test this method, a 15 layer atmosphere with the Lambertian boundary at the bottom is considered. One layer contains a cloud. The ozone profile is modeled as in [31]. The computations are performed in the spectral range 280-335 nm. The results are 6 D.S. Efremenko et al. Fig. 2. Examples of calculated phase functions, which are used for training. The blue curve corre- sponds to the Rayleigh phase function, while the black curve is related to the Mie computations. compared against a reference solution, in which full multi-stream RTM with Ndo = 16 is used. The examples of spectra are shown in Fig. 3. The error is computed for each spectrum in each spectral point. As the radiance values around λ =290 nm are very small comparing to radiances at longer wavelengths, the relative error e is computed as follows: L (λ) − Lref (λ) e (λ) = , (9) Lref (315 nm) where Lref stands for the reference model. The histogram of errors is shown in Fig. 3. As we can see, the error of this approach is within 0.1 %. The overall computational time is reduced by approximately 5 times as compared to the reference RTM. Note that as the ANN part is much faster than the matrix opera- tor part, the latter becomes a performance bottleneck. This scheme appears to be quite stable, since no cases with high errors have been captured. We should note that the performance enhancement is not significant (up to 5 times) in comparison to 1000-fold performance enhancement, e.g. in [7, 23]. However, this approach is flexible and the training procedure is fast. In addition, one should note that by replacing the eigenvalue solver with the ANN, the computational time of the RTM is mostly due to the time re- quired to the matrix operator method. That makes it comparable to closed-form RTMs, in which the eigenvalues are found analytically (e.g. the two-stream RTMs). Two Machine Learning Based Schemes for Solving Direct and... 7 Fig. 3. (left) Examples of simulated spectra in the Hartley-Huggins band; (right) the histogram of errors of the method, in which the eigenvalue part of the RTM is substituted with the ANN. 4 Inverse Model Parametrization Using Real Measurements Now we consider the ANN-based approach for accelerating the inverse operator. Fol- lowing Tikhonov [29], the retrieval algorithm finds the state vector x of medium pa- rameters that minimizes the discrepancy between the simulated and measured spectra in the following sense: n o 2 x = arg min ky − RTM(x)k + αΩ (x) , (10) x where y is the vector of measurements (in the case of atmospheric retrievals, spectral radiances at the top of the atmosphere), α is the regularization parameter and Ω is the stabilizing functional. The expression in curly brackets is referred to as the Tikhonov function. In principle, RTM in Eq. (10) can be substituted with the ANN, i.e. n o 2 x = arg min ky − ANN(x)k + αΩ (x) . (11) x Then the minimization problem can be solved either by using the local or global op- timization algorithms. In the former case, ANN has to be linearized to compute the Jacobian matrix [25]. Alternative approach is to train ANN in the backward direction, as shown in Fig. 4. The main advantage of the backward trained ANNs over the classi- cal optimization approach (Eq. (10)) is that the resulting inverse operator can be readily applied to real measurements to retrieve the parameter of interest. The concept of back- ward ANNs was applied to several problems of remote sensing, including CO2 [21] and SO2 plume height [16, 18] and surface albedo [24] and ozone profile [31] retrieval. However, the retrieval problem is ill-posed and the resulting ANN has to be reg- ularized. Essentially, the instability comes from the facts that the inverse operator is highly sensible to perturbations in the measurement data space. As RTMs cannot re- produce accurately the measurements (due to noise in real data, instrumental artifacts and imperfection of the physical model (RTM)), there is always a discrepancy between the simulated and real data. This discrepancy propagates (and grows) into the space of parameters to be retrieved. Consequently, the inverse operator has to be regularized. 8 D.S. Efremenko et al. Fig. 4. Scheme of the ANN training to parameterize the inverse operator. Fig. 5. Training of ANN for retrieval of atmospheric constituents using real measured spectral radiances and parameters retrieved by using conventional retrieval algorithm. The alternative approach is to use real measurements in training, in particular, the data already processed by using traditional retrieval methods. This concept is illustrated in Fig. 5. The ANN is trained in the reverse mode, i.e. the sun-normalized spectral ra- diances (Level-1 data) are fed to the ANN as input, while the retrieved atmospheric parameters (Level-2 data) are considered as outputs. Since Level-2 data has been re- trieved using the conventional RTM and regularization procedures (according to Eq. (10)), the RTM is captured by the ANN implicitly, and thus, the predictions of the ANN are still based on physical models comprising instrument-specific features (e.g. noise, offset, wavelength calibration, detector degradation issues, etc.). To examine this approach, we consider the TROPOMI spectra (date 04.04.18) in the wavelength range 325-335 nm, available at [1]. The spectra together along with the geometry parameters (i.e. solar zenith angle, viewing zenith angle, and relative azimuth angle) and the surface albedo are used to train and test the ANN. The output of the ANN is the ozone total column value. The data is being shuffled and divided into 70% train Two Machine Learning Based Schemes for Solving Direct and... 9 Fig. 6. The histogram of differences between retrieved ozone total column values by using ANNs and those retrieved by using the conventional approach based on the Tikhonov regularization. and 30% test data sets. The ANN configuration with 3 hidden layers with 50, 15 and 5 neurons respectively is trained with 185000 samples and tested on 79067 samples. The histogram of differences between retrieved ozone total column values by using ANNs and those retrieved by using the conventional approach is shown in Fig. 6. The mean relative difference is of about 0.6% (i.e. below 3 DU). Note that the error increases if the ANN trained with the data from one day is applied to another date. Therefore, there are two possible applications of this method. First, to train ANN by using a part of the available data per day, and use the trained ANN to process the rest of the data. Secondly, ANN can be used for finding a first estimation in the iterative retrieval procedure. 5 Conclusions We have analyzed several machine learning based schemes for solving direct and in- verse problems of radiative transfer. The conventional RTMs are the most flexible tools for simulating and processing the optical data; however, the flexibility comes at high computational costs. By substituting the eigenvalue computations in the discrete ordi- nate radiative transfer models by neural networks, we obtain a flexible tool with mod- erate performance enhancement. However, the training datasets are relatively small and the training procedure is not time consuming. To parameterize the inverse RTM operator, the ANNs can be trained in the backward direction. Although the ANN is computationally very fast, some sort of regularization is still required. We have considered the ANN trained using the real radiance spectra instead of synthetic data. The key point here is that the training is performed by using real measurements, which have been already processed via conventional retrieval algo- rithms with full RTM simulations. Thus, the trained ANN already captures the physics behind the measurement process (expressed by the RTM) as well as the instrument- related features. This approach seems to be the most stable from considered techniques. However, a big amount of data is required for the training procedure. 10 D.S. Efremenko et al. Taking into account the high computational performance of ANNs, our future re- search will be focused on applying machine learning based approaches to more accu- rate radiative transfer models taking into account stochastic cloud inhomogeneities of the atmosphere [13, 15] as well as three-dimensional effects [11]. References 1. Sentinel-5P Pre-Operationals Data Hub, https://s5phub.copernicus.eu/dhus, [accessed 08.02.2020] 2. Afanas’ev, V., Basov, A.Y., Budak, V., Efremenko, D., Kokhanovsky, A.: Analysis of the dis- crete theory of radiative transfer in the coupled “ocean–atmosphere” system: Current status, problems and development prospects. Journal of Marine Science and Engineering 8(3), 202 (2020). https://doi.org/10.3390/jmse8030202 3. del Águila, A., Efremenko, D.S., Trautmann, T.: A review of dimensionality reduction techniques for processing hyper-spectral optical signal. Light & Engineering 27(3), 85–98 (2019). https://doi.org/10.33383/2019-017 4. Bodhaine, B., Wood, N., Dutton, E., Slusser, J.: On Rayleigh optical depth calcu- lations. Journal of Atmospheric and Oceanic Technology 16(11), 1854–1861 (1999). https://doi.org/10.1175/1520-0426(1999)016¡1854:orodc¿2.0.co;2 5. Bohren, C., Huffman, D.: Absorption and Scattering of Light by Small Particles. Wiley (1998). https://doi.org/10.1002/9783527618156 6. Budak, V., Klyuykov, D., Korkin, S.: Complete matrix solution of radiative transfer equation for pile of horizontally homogeneous slabs. J Quant Spectrosc Radiat Transfer 112(7), 1141– 1148 (2011). https://doi.org/10.1016/j.jqsrt.2010.08.028 7. Bue, B.D., Thompson, D.R., Deshpande, S., Eastwood, M., Green, R.O., Natraj, V., Mullen, T., Parente, M.: Neural network radiative transfer for imaging spectroscopy. Atmospheric Measurement Techniques 12(4), 2567–2578 (2019). https://doi.org/10.5194/amt-12-2567- 2019 8. Chandrasekhar, S.: Radiative trasnfer. Dover publications, inc. New York (1950) 9. Cybenko, G.: Approximation by superpositions of a sigmoidal function. Mathematics of Control, Signals, and Systems 2(4), 303–314 (1989). https://doi.org/10.1007/bf02551274 10. Deirmendjian, D.: Electromagnetic Scattering on Spherical Polydispersions. Elsevier (1969) 11. Doicu, A., Efremenko, D.S.: Linearizations of the spherical harmonic discrete ordinate method (SHDOM). Atmosphere 10(6), 292 (2019). https://doi.org/10.3390/atmos10060292 12. Doicu, A., Trautmann, T.: Discrete-ordinate method with matrix exponential for a pseudo- spherical atmosphere: Scalar case. Journal of Quantitative Spectroscopy and Radiative Trans- fer 110(1-2), 146–158 (2009). https://doi.org/10.1016/j.jqsrt.2008.09.014 13. Doicu, A., Efremenko, D.S., Loyola, D., Trautmann, T.: Discrete ordinate method with ma- trix exponential for stochastic radiative transfer in broken clouds. Journal of Quantitative Spectroscopy and Radiative Transfer 138, 1–16 (2014) 14. Efremenko, D., Doicu, A., Loyola, D., Trautmann, T.: Acceleration techniques for the dis- crete ordinate method. Journal of Quantitative Spectroscopy and Radiative Transfer 114, 73–81 (2013). https://doi.org/10.1016/j.jqsrt.2012.08.014 15. Efremenko, D.S., Doicu, A., Loyola, D., Trautmann, T.: Fast stochastic radiative trans- fer models for trace gas and cloud property retrievals under cloudy conditions. In: Springer Series in Light Scattering, pp. 231–277. Springer International Publishing (2018). https://doi.org/10.1007/978-3-319-70796-9 3 Two Machine Learning Based Schemes for Solving Direct and... 11 16. Efremenko, D.S., Loyola, D., Hedelt, P., Spurr, R.J.D.: Volcanic SO2 plume height retrieval from UV sensors using a full-physics inverse learning ma- chine algorithm. International Journal of Remote Sensing 38(sup1), 1–27 (2017). https://doi.org/10.1080/01431161.2017.1348644 17. Efremenko, D., Molina Garcı́a, V., Gimeno Garcı́a, S., Doicu, A.: A review of the matrix- exponential formalism in radiative transfer. Journal of Quantitative Spectroscopy and Radia- tive Transfer 196, 17–45 (2017). https://doi.org/10.1016/j.jqsrt.2017.02.015 18. Hedelt, P., Efremenko, D., Loyola, D., Spurr, R., Clarisse, L.: Sulfur dioxide layer height retrieval from Sentinel-5 Precursor/TROPOMI using FP ILM. Atmospheric Measurement Techniques 12(10), 5503–5517 (2019). https://doi.org/10.5194/amt-12-5503-2019 19. Hornik, K.: Approximation capabilities of multilayer feedforward networks. Neural Net- works 4(2), 251–257 (1991). https://doi.org/10.1016/0893-6080(91)90009-t 20. Kataev, M., Lukyanov, A.: Simulation of reflected solar radiation for atmosphere gas com- position evaluation for optical remote sensing from space. Light & Engineering 26(3), 14–21 (2018) 21. Kataev, M., Lukyanov, A., Bekerov, A.: Modification of the empirical orthogonal functions method for solving the inverse task of retrieving of the CO2 total content from satellite data. Journal of Siberian Federal University. Engineering & Technologies 11(1), 77–85 (2018). https://doi.org/10.17516/1999-494x-0011 22. Levit, G.S., Krumbein, W.E., and, R.G.: Space and time in the works of V. I. Vernadsky. Environmental Ethics 22(4), 377–396 (2000). https://doi.org/10.5840/enviroethics20002244 23. Loyola, D.G., Gimeno Garcı́a, S., Lutz, R., Argyrouli, A., Romahn, F., Spurr, R.J.D., Ped- ergnana, M., Doicu, A., Molina Garcı́a, V., Schüssler, O.: The operational cloud retrieval algorithms from TROPOMI on board Sentinel-5 Precursor. Atmospheric Measurement Tech- niques 11(1), 409–427 (2018). https://doi.org/10.5194/amt-11-409-2018 24. Loyola, D., Xu, J., Heue, K.P., Zimmer, W.: Applying FP ILM to the retrieval of geometry- dependent effective Lambertian equivalent reflectivity (GE LER) daily maps from UVN satellite measurements. Atmospheric Measurement Techniques 13(2), 985–999 (2020). https://doi.org/10.5194/amt-13-985-2020 25. Molina Garcı́a, V., Sasi, S., Efremenko, D.S., Doicu, A., Loyola, D.: Linearized radia- tive transfer models for retrieval of cloud parameters from EPIC/DSCOVR measure- ments. Journal of Quantitative Spectroscopy and Radiative Transfer 213, 241–251 (2018). https://doi.org/10.1016/j.jqsrt.2018.03.008 26. Natraj, V.: A review of fast radiative transfer techniques. In: Kokhanovsky, A. (ed.) Light scattering reviews, vol. 8, pp. 475–504. Springer Berlin Heidelberg (2013). https://doi.org/10.1007/978-3-642-32106-1 10 27. Plass, G., Kattawar, G., Catchings, F.: Matrix operator theory of radiative transfer 1: Rayleigh scattering. Applied Optics 12(2), 314 (1973). https://doi.org/10.1364/ao.12.000314 28. Stamnes, K., Tsay, S., Wiscombe, W., Jayaweera, K.: Numerically stable algorithm for discrete-ordinate-method radiative transfer in multiple scattering and emitting layered media. Appl Opt 12, 2502–2509 (1988). https://doi.org/10.1364/AO.27.002502 29. Tikhonov, A., Arsenin, V.: Solutions of ill-posed problems. Scripta series in mathematics, Winston & Sons, Washington, D.C. (1977) 30. Veefkind, J., Aben, I., McMullan, K., Förster, H., de Vries, J., Otter, G., Claas, J., Eskes, H., de Haan, J., Kleipool, Q., van Weele, M., Hasekamp, O., Hoogeveen, R., Landgraf, J., Snel, R., Tol, P., Ingmann, P., Voors, R., Kruizinga, B., Vink, R., Visser, H., Levelt, P.: TROPOMI on the ESA Sentinel-5 Precursor: A GMES mission for global observations of the atmo- spheric composition for climate, air quality and ozone layer applications. Remote Sensing of Environment 120, 70–83 (2012). https://doi.org/10.1016/j.rse.2011.09.027 12 D.S. Efremenko et al. 31. Xu, J., Schussler, O., Loyola, D.G., Romahn, F., A.Doicu: A novel ozone profile shape retrieval using full-physics inverse learning machine (FP-ILM). IEEE Journal of Selected Topics in Applied Earth Observations and Remote Sensing 10(12), 5442–5457 (2017). https://doi.org/10.1109/jstars.2017.2740168