Artificial neural networks (ANNs) are used to substitute computationally expensive radiative transfer models (RTMs) and inverse operators (IO) for retrieving optical parameters of the medium. However, the direct parametrization 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 results and those published by other authors. In particular, two techniques are proposed. In the first method, the ANN substitutes the eigenvalue solver in the discrete 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 provides robust results even without applying the Tikhonov regularization method.
Machine learning techniques have become of paramount importance in the geosciences,
lighting engineering and remote sensing [
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
computations 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
properties 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 [
ANNs are used to parameterize an inverse operator in retrieval problems. The
advantage 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
procedures. 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 [
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
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
scattering, multiple scattering plays an important role. The one-dimensional radiative transfer
equation for the radiance L reads as follows [
L ( ; ; ') +
2 1 ! ( ) Z Z 4 0 1 p ( ; 0; ; ' '0) L ( ; 0; '0) d 0d'0; (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. [
We consider the cosine azimuthal expansion of the radiance field and the phase function providing an equation for the azimuthal component Lm ( ; ). Then, the radiance field is discretized in the -domain by considering Ndo Gaussian points f ig per hemisphere and corresponding weights fwig, 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 discretization 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 scattering 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): where d d i" ( ) i# ( ) =
i" ( ) Al i# ( ) + bl ( ) ; l
l+1;
i"# ( ) i = Lm ( ;
i) ; i = 1; :::; Ndo;
while A and b are the so called layer matrix and layer source vector, respectively
(expressions for them can be found in [
Z l 0 eAltbl (t) dt: To evaluate eAl l , the eigenvalue decomposition method can be applied: eAlt = Vle ltV 1 l ; 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 il#+1 =
Rl" Tl" " il# # Tl# Rl# il"+1 +
P"
Pl# ;
l
where R"# and T"# are the reflection and transmission matrices, respectively, while
l l
Pl"# is the source term. They can be expressed in terms of Vl and l (ready-to-compute
formulas are given in [
By using the matrix operator method [
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 dependencies 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 radiance 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
[
Two Machine Learning Based Schemes for Solving Direct and... 5
inputs. The particle size distribution is described by the modified gamma distribution
[
f (r) = ar e br ; 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, and thus, f (r) is defined by three independent parameters. The examples of the calculated 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 inputs 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 configuration 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 samples. 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 dependencies 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
[
L ( )
Lref ( ) Lref (315 nm) ; (9) 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
operator 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 [
Two Machine Learning Based Schemes for Solving Direct and... 7
Now we consider the ANN-based approach for accelerating the inverse operator.
Following Tikhonov [
RTM(x)k2 + (x)o ; (10) (11) 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. x = arg min n x ky
ANN(x)k2 +
(x)o :
Then the minimization problem can be solved either by using the local or global
optimization algorithms. In the former case, ANN has to be linearized to compute the
Jacobian matrix [
However, the retrieval problem is ill-posed and the resulting ANN has to be regularized. Essentially, the instability comes from the facts that the inverse operator is highly sensible to perturbations in the measurement data space. As RTMs cannot reproduce 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.
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 radiances (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 retrieved 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 [
Two Machine Learning Based Schemes for Solving Direct and... 9 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
We have analyzed several machine learning based schemes for solving direct and inverse 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 ordinate radiative transfer models by neural networks, we obtain a flexible tool with moderate 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 algorithms with full RTM simulations. Thus, the trained ANN already captures the physics behind the measurement process (expressed by the RTM) as well as the instrumentrelated features. This approach seems to be the most stable from considered techniques. However, a big amount of data is required for the training procedure.
Taking into account the high computational performance of ANNs, our future
research will be focused on applying machine learning based approaches to more
accurate radiative transfer models taking into account stochastic cloud inhomogeneities of
the atmosphere [
Two Machine Learning Based Schemes for Solving Direct and... 11