Next Article in Journal
Atmospheric Correction of Airborne Hyperspectral CASI Data Using Polymer, 6S and FLAASH
Next Article in Special Issue
Estimation of Salinity Content in Different Saline-Alkali Zones Based on Machine Learning Model Using FOD Pretreatment Method
Previous Article in Journal
TomoSAR 3D Reconstruction for Buildings Using Very Few Tracks of Observation: A Conditional Generative Adversarial Network Approach
 
 
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

An Overview of Neural Network Methods for Predicting Uncertainty in Atmospheric Remote Sensing

1
Remote Sensing Technology Institute, German Aerospace Center (DLR), 82234 Oberpfaffenhofen, Germany
2
Institute of Mathematics, University of Augsburg, 86159 Augsburg, Germany
*
Author to whom correspondence should be addressed.
Remote Sens. 2021, 13(24), 5061; https://doi.org/10.3390/rs13245061
Submission received: 4 November 2021 / Revised: 8 December 2021 / Accepted: 10 December 2021 / Published: 13 December 2021
(This article belongs to the Special Issue Recent Advances in Neural Network for Remote Sensing)

Abstract

:
In this paper, we present neural network methods for predicting uncertainty in atmospheric remote sensing. These include methods for solving the direct and the inverse problem in a Bayesian framework. In the first case, a method based on a neural network for simulating the radiative transfer model and a Bayesian approach for solving the inverse problem is proposed. In the second case, (i) a neural network, in which the output is the convolution of the output for a noise-free input with the input noise distribution; and (ii) a Bayesian deep learning framework that predicts input aleatoric and model uncertainties, are designed. In addition, a neural network that uses assumed density filtering and interval arithmetic to compute uncertainty is employed for testing purposes. The accuracy and the precision of the methods are analyzed by considering the retrieval of cloud parameters from radiances measured by the Earth Polychromatic Imaging Camera (EPIC) onboard the Deep Space Climate Observatory (DSCOVR).

Graphical Abstract

1. Introduction

In atmospheric remote sensing, the retrieval of atmospheric parameters is an inverse problem that is usually ill posed. Due to its ill posedness, measurement errors can lead to large errors in the retrieved quantities. It is therefore desirable to characterize the retrieved value by an estimate of uncertainty describing a range of values that probably produce a measurement [1].
The retrieval algorithms are mostly based on deterministic or stochastic optimization methods. From the first category, the method of Tikhonov regularization and iterative regularization methods deserve to be mentioned, while from the second category, the Bayesian methods are the most representative. The Bayesian framework [2,3] provides an efficient way to deal with the ill-posedness of the inverse problem and its uncertainties. In this case, the solution of the inverse problem is given by the a posteriori distribution (the conditional distribution of the retrieval quantity given the measurement) that accounts for all assumed retrieval uncertainties. Under the assumptions that (i) the a priori knowledge and the measurement uncertainty are both Gaussian distributions, and (ii) the forward model is moderately nonlinear, the a posteriori distribution is approximately Gaussian with a covariance matrix that can be analytically calculated. However, even neglecting the validity of these assumptions, the method is not efficient for the operational processing of large data volumes. The reason is that the computations of the forward model and its Jacobian are expensive computational processes.
Compared to Bayesian methods, neural networks are powerful tools for the design of efficient retrieval algorithms. Their capability to approximate any continuous function on a compact set to an arbitrary accuracy makes them well suited to approximate the input–output function represented by a radiative transfer model. An additional feature is that the derivatives of a neural network model with respect to its inputs can be analytically computed. Thus, a neural network algorithm can produce, in addition to the approximation of the radiative quantities of interest, also an estimation of their derivatives with respect to the model inputs. While the training of a neural network may require a significant amount of time, a trained neural network may deliver accurate predictions of the forward model and its Jacobian, typically in a fraction of a millisecond.
Neural networks were initially developed for the emulation of radiative transfer models [4,5,6,7,8] and subsequently for atmospheric remote sensing. The latter include:
  • Neural networks that are used to emulate the forward model and are applied in conjunction with a Bayesian approach to solve the inverse problem [9,10,11,12];
  • Neural networks that are developed to directly learn the retrieval mappings from data [13,14,15,16,17,18,19,20]; and
  • Neural networks that are designed to directly invert the atmospheric parameters of interest which are then used as initial values in an optimization algorithm, e.g., the method of Tikhonov regularization [21,22].
In the first case, the total uncertainty sums up the contributions of the uncertainties in the data and in the neural network model, whereby the model uncertainty is computed from the statistics of the errors over the data set. However, the retrieval algorithm is still based on the assumption that the forward model is nearly linear, which is generally not true. In the second case, the probabilistic character of the inverse problem is neglected and uncertainty estimates are provided as mean errors computed over the data set; this is a disadvantage compared to Bayesian methods. Remarkable exceptions to the approaches listed above are the works of Aires [23,24], in which a Bayesian framework in conjunction with the Laplace approximation was used to model the retrieval errors (estimated from the error covariance matrix observed on the data set), and of Pfreundschuh et al. [25], in which a quantile regression neural network was designed. Unfortunately, the Laplace approximation is only suitable for small training sets and simple networks, while quantile regression is suitable for the retrieval of scalar quantities (it is unclear whether a reasonable approximation of the quantile contours of the joint a posteriori distribution can be obtained).
This paper is devoted to an overview of neural networks methods for predicting uncertainty in atmospheric remote sensing. In addition to the method based on a neural network for simulating the radiative transfer model and a Bayesian approach for solving the inverse problem (Case 1 above), methods relying on Bayesian networks are described. These methods, in which the network activations and weights can be modeled by parametric probability distributions, are standard tools for uncertainty predictions in deep neural networks, and can be applied to nonlinear retrieval problems, and to the best knowledge of the authors, have not yet been used in atmospheric remote sensing. The paper, which is merely of pedagogical nature, is organized as follows. In Section 1, we present the theoretical background of this study. This is then used in Section 2 to develop several neural network methods and apply them to a specific cloud parameters retrieval problem. Section 3 contains a few concluding remarks.

2. Theoretical Background

We consider a generic model y = F ( x ) , where x R N x is the input vector, F is some deterministic function and y R N y is the output vector. For an atmosphere characterized by a set of state parameters, the signal measured by an instrument at different wavelengths can be computed by a radiative transfer model R . Specifically, we will refer to the measurement signals as data, and split the vector of state parameters in (i) the vector of atmospheric parameters that are intended to be retrieved, and (ii) the vector of atmospheric parameters that are known to have some accuracy, but are not included in the retrieval (forward model parameters). In this study, we will use neural networks to model the radiative transfer function R , as well as its inverse R 1 . In this regard and in order to simplify the notation, we will consider two cases. In the first case, referred to as the direct problem, the input x is the set of atmospheric parameters, the output y is the set of data and the forward model F coincides with the radiative transfer model R , while in the second case, referred to as the inverse problem, the situation is reversed: the input x includes the sets of data and forward model parameters, while the output y includes the set of atmospheric parameters to be retrieved (in the absence of forward model parameters, the forward model F reproduces the inverse of the radiative transfer model R 1 ).
In machine learning, the task is to approximate F ( x ) by a neural network model f ( x , ω ) characterized by a set of parameters ω [26,27]. For doing this, we consider a set of inputs X = { x ( n ) } n = 1 N and a corresponding set of outputs Y = { y ( n ) } n = 1 N , given by y ( n ) = F ( x ( n ) ) , where N is the number of samples. In a regression problem, D = { ( x ( n ) , y ( n ) ) } n = 1 N forms a data set—or more precisely, a training set—from which the neural network model f ( x , ω ) can be inferred. Traditional neural networks are comprised of units or nodes arranged in an input layer, an output layer, and a number of hidden layers situated between the input and output layers. Let L + 1 be the number of layers, and N l be the number of units in layer l, where l = 0 , , L . The input layer corresponds to l = 0 , the output layer to l = L , so that N x = N 0 and N y = N L . In feed-forward networks, the signals y i , l 1 from units i = 1 , , N l 1 in layer l 1 are multiplied by a set of weights w j i , l , j = 1 , , N l , i = 1 , , N l 1 , and then summed and combined with a bias b j , l , j = 1 , , N l . This calculation forms the pre-activation signal u j , l = i = 1 N l 1 w j i , l y i , l 1 + b j , l which is transformed by the layer activation function g l to form the activation signal y j , l of unit j = 1 , , N l in layer l. Defining the matrix of weights W l R N l × N l 1 and the vector of biases b l R N l by [ W l ] j i = w j i , l and [ b l ] j = b j , l , respectively, and letting ω = { W l , b l } l = 1 L be the set of network parameters, the feed-forward operations can be written in matrix form as
y 0 = x ,
u l = W l y l 1 + b l ,
y l = g l ( u l ) , l = 1 , , L ,
f ( x , ω ) = y L ,
where [ y l ] i = y i , l and [ u l ] j = u j , l . Deep learning is the process of regressing the network parameters ω on the data D . The standard procedure is to compute a point estimate ω ^ as the minimizer of some loss function by using the back-propagation algorithm [28]. In a stochastic framework, the loss function is usually defined as the log likelihood of the data set, with an eventual regularization term to penalize the network parameters. From a statistical point of view, this procedure is equivalent to a maximum a posteriori (MAP) estimation when regularization is used, and maximum likelihood estimation (MLE) when this is not the case.
In this section, we review the basic theory which serves as a basis for the development of different neural network architectures. In particular, we describe (i) the methodology for computing point estimates; (ii) the different types of uncertainty; and (iii) Bayesian networks.

2.1. Point Estimates

A data model with output noise is given by
y = f ( x , ω ) + δ y ,
δ y N ( 0 , C y δ ) ,
where here and in the following, the notation N ( x ; x ¯ , C x ) , or more simply and when no confusion arises, N ( x ¯ , C x ) stands for a Gaussian distribution with the mean x ¯ and covariance matrix C x . When the true input x is hidden (so that it cannot be observed) but samples from a random vector z = x + δ x with δ x N ( 0 , C x δ ) , i.e., p ( z | x ) = N ( z ; x , C x δ ) are available, the pertinent model is the data model with input and output noise, that is:
y = f ( z , ω ) + Δ y ,
Δ y N ( 0 , C y δ ) .
The error Δ y sums up the contributions of the output error and of the input error that propagates through the network in the output space. Specifically, when the noise process in the input space is small and the linearization:
f ( x , ω ) = f ( z , ω ) + K x ( z , ω ) ( x z ) ,
K x ( z , ω ) = f x ( z , ω ) ,
is assumed, we find (cf. Equations (5), (7), and (9)) Δ y = δ y + K x ( z , ω ) ( x z ) , and further:
C y δ ( z , ω ) = C y δ + K x ( z , ω ) C x δ K x T ( z , ω ) .
To design a neural network, we consider a data set D associated to each data model, meaning that:
  • An exact data set D = { ( x ( n ) , y ( n ) ) } n = 1 N , where y ( n ) = F ( x ( n ) ) ;
  • A data set with output noise D = { ( x ( n ) , y ( n ) ) } n = 1 N , where y ( n ) = F ( x ( n ) ) + δ y ; and
  • A data set with input and output noise D = { ( z ( n ) , y ( n ) ) } n = 1 N , where z ( n ) = x ( n ) + δ x and y ( n ) = F ( x ( n ) ) + δ y .
In a stochastic framework, a neural network can be regarded as a probabilistic model p ( y | z , ω ) ; given an observable input z and a set of parameters ω , a neural network assigns a probability to each possible output y . In view of Equations (7) and (8), the a priori confidence in the predictive power of the model is given by
p ( y | z , ω ) = N ( y ; f ( z , ω ) , C y δ ( z , ω ) ) .
The process of learning from the data D can be described by the posterior p ( ω | D ) = p ( ω | Z , Y ) , which represents the Bayes plausibility for the parameters ω given the data D . This can be estimated by using the Bayesian rule:
p ( ω | D ) = p ( D | ω ) p ( ω ) p ( D ) p ( D | ω ) p ( ω ) exp [ E ( ω ) ] ,
where p ( D | ω ) is the likelihood or the probability of the data, p ( ω ) is the prior over the network parameters, p ( D ) = p ( D | ω ) p ( ω ) d ω the evidence, and:
E ( ω ) = E D ( ω ) + E R ( ω )
is the loss function. The first term E D ( ω ) in the expression of the loss function E ( ω ) is the contribution from the likelihood p ( D | ω ) , written as the product (cf. Equation (12)):
p ( D | ω ) = p ( Y | Z , ω ) = n = 1 N p ( y ( n ) | z ( n ) , ω ) exp [ E D ( ω ) ] ,
E D ( ω ) = 1 2 n = 1 N [ y ( n ) f ( z ( n ) , ω ) ] T [ C y δ ( z ( n ) , ω ) ] 1 [ y ( n ) f ( z ( n ) , ω ) ] ,
while the second term E R ( ω ) is the contribution from the prior p ( ω ) , chosen for example, as the Gaussian distribution:
p ( ω ) = N ( ω ; 0 , C ω ) exp [ E R ( ω ) ] ,
E R ( ω ) = 1 2 ω T C ω 1 ω .
In this regard, point estimates with regularization are computed by maximizing the posterior p ( ω | D ) :
ω ^ = ω MAP = arg max ω log p ( ω | D ) = arg min ω E ( ω ) ,
while point estimates without regularization are computed by maximizing the likelihood p ( D | ω ) :
ω ^ = ω MLE = arg max ω log p ( D | ω ) = arg min ω E D ( ω ) .
Some comments are in order.
  • For a data model with output noise, we have z = x ,   C y δ = C y δ , and D = { ( x ( n ) , y ( n ) ) } n = 1 N with y ( n ) = F ( x ( n ) ) + δ y . Moreover, for the covariance matrix C y δ = σ y 2 I , where I is the identity matrix, we find:
    E D ( ω ) = n = 1 N 1 2 σ y 2 | | y ( n ) f ( z ( n ) , ω ) | | 2 ,
    or more precisely:
    E D ( ω ) = n = 1 N 1 2 σ y 2 | | y ( n ) f ( z ( n ) , ω ) | | 2 + N y 2 log σ y 2 ,
    Assuming C ω = σ ω 2 I and using Equation (21), we infer that the point estimate ω ^ = ω MAP is the minimizer of the Tikhonov function:
    E ( ω ) = 1 2 n = 1 N | | y ( n ) f ( x ( n ) , ω ) | | 2 + α | | ω | | 2 ,
    where α = σ y 2 / ( 2 σ ω 2 ) is the regularization parameter.
  • A model with exact data can be handled by considering the data model with output noise and by making σ y 2 0 in the representation of the data error covariance matrix C y δ = σ y 2 I . For σ y 2 0 , the relation y ( n ) = F ( x ( n ) ) + δ y yields y ( n ) F ( x ( n ) ) , while Equation (23) and the relation α = σ y 2 / ( 2 σ ω 2 ) 0 gives:
    ω ^ ω MLE = arg min ω E ( ω ) , E ( ω ) = 1 2 n = 1 N | | y ( n ) f ( x ( n ) , ω ) ] | | 2 .
    Thus, when learning a neural network with the exact data, the maximum likelihood estimate minimizes the sum of square errors. Note that in this case, regularization is not absolutely required, because the output data are exact.
  • For a data model with input and output noise, the computation of the estimate ω ^ is not a trivial task, because the covariance matrix C y δ ( z , ω ) , which enters in the expression of E D ( ω ) , depends on ω . Moreover, in Equation (11), C y δ ( z , ω ) corresponds to a linearization of the neural network function under the assumption that the noise process in the input space is small. This problem can be solved by implicitly learning the covariance matrix C y δ ( z , ω ) from the loss function [29]. Specifically, we assume that C y δ ( z , ω ) is a diagonal matrix with entries σ j 2 ( z , ω ) , that is:
    C y δ ( z , ω ) = diag [ σ j 2 ( z , ω ) ] j = 1 N y ,
    implying:
    log p ( y | z , ω ) j = 1 N y 1 2 σ j 2 ( z , ω ) [ y j μ j ( z , ω ) ] 2 + 1 2 log σ j 2 ( z , ω ) .
    Here, we identified μ f , and set μ j = [ μ ] j and y j = [ y ] j . To learn the variances σ j 2 ( z , ω ) from the loss function, we use a single network with input z , and output [ μ j ( z , ω ) , σ j 2 ( z , ω ) ] R 2 N y ; thus, in the output layer, N y units are used to predict μ j and N y units to predict σ j 2 . In practice, to increase numerical stability, we train the network to predict the log variance ρ j = log σ j 2 , in which case, the likelihood loss function is:
    E D ( ω ) = n = 1 N E D ( n ) ( ω ) ,
    with:
    E D ( n ) ( ω ) = 1 2 j = 1 N y { exp [ ρ j ( z ( n ) , ω ) ] [ y j ( n ) μ j ( z ( n ) , ω ) ] 2 + ρ j ( z ( n ) , ω ) } .
    It should be pointed out that from Equation (25) that it is apparent that the model is stopped from predicting high uncertainty through the log σ j 2 term, but also low uncertainty for points with high residual error, as low σ j 2 will increase the contribution of the residual. On the other hand, it should also be noted that a basic assumption of the model is that the covariance matrix C y δ ( z , ω ) is diagonal, which unfortunately, is contradicted by Equation (11).
  • In order to generate a data set with input and output noise, that is, D = { ( z ( n ) , y ( n ) ) } n = 1 N , where z ( n ) = x ( n ) + δ x , δ x N ( 0 , C x δ ) , and y ( n ) = F ( x ( n ) ) + δ y , we used the jitter approach. According to this approach, at each forward pass through the network, a new random noise δ x is added to the original input vector x ( n ) . In other words, each time a training sample is passed through the network, random noise is added to the input variables, making them different every time it is passed through the network. By this approach, which is a simple form of data augmentation, the dimension of the data set is reduced (actually, the data set is D = { ( x ( n ) , y ( n ) ) } n = 1 N with y ( n ) = F ( x ( n ) ) + δ y ).

2.2. Uncertainties

In our analysis, we are interested in estimating the uncertainty associated with the underlying processes. The quantity which exactly quantifies the model’s uncertainty is the predictive distribution of an unknown output y given an observable input z , defined by
p ( y | z , D ) = p ( y | z , ω ) p ( ω | D ) d ω .
If p ( y | z , D ) is known, the first two moments of the output y can be computed as
E ( y ) = y p ( y | z , D ) d y ,
E ( y y T ) = y y T p ( y | z , D ) d y ,
and the covariance matrix as
Cov ( y ) = E ( y y T ) E ( y ) E ( y ) T .
From Equation (28), we see that the predictive distribution p ( y | z , D ) can be computed if the Bayesian posterior p ( ω | D ) is known. Unfortunately, computing the distribution p ( ω | D ) by means of Equation (13) is usually an intractable problem, because computing the evidence p ( D ) = p ( D | ω ) p ( ω ) d ω is not a trivial task. To address this problem, either:
  • The Laplace approximation, which yields an approximate representation for the posterior p ( ω | D ) ; or
  • A variational inference approach, which learns a variational distribution q θ ( ω ) to approximate the posterior p ( ω | D ) ,
can be used. In the first case, we are dealing with deterministic (point estimate) neural networks, in which a single realization of the network parameters ω is learned, while in the second case, we are dealing with stochastic neural networks, in which a distribution over the network parameters ω is learned.
The Laplace approximation is of theoretical interest because it provides explicit representations for the different types of uncertainty. In Appendix A it is shown that under the Laplace approximation, the predictive distribution is given by [23,24,30]
p ( y | z , D ) exp 1 2 [ y f ( z , ω ^ ) ] T C y 1 ( z , ω ^ ) [ z f ( z , ω ^ ) ] ,
C y ( z , ω ^ ) = C y δ ( z , ω ^ ) + C y e ( z , ω ^ ) ,
C y e ( z , ω ^ ) = K ω ( z , ω ^ ) H 1 ( ω ^ ) K ω T ( z , ω ^ ) ,
where:
K ω ( z , ω ) = f ω ( z , ω )
is the Jacobian of f with respect to ω and:
[ H ( ω ^ ) ] i j = E ω i ω j ( ω ^ ) ,
is the Hessian matrix of the loss function. Equation (32) provides a Gaussian approximation to the predictive distribution p ( y | z , D ) with mean f ( z , ω ^ ) and covariance matrix C y ( z , ω ^ ) . From Equation (33), we see that the covariance matrix C y ( z , ω ^ ) has two components.
  • The first component C y δ ( z , ω ^ ) is the covariance matrix in the distribution over the error in the output y . This term, which is input dependent, describes the so-called aleatoric heteroscedastic uncertainty measured by p ( y | z , ω ) . For a data model with output noise, we have (cf. Equation (11)) C y δ = C y δ , and this term, which is input independent, describes the so-called aleatoric homoscedastic uncertainty measured by p ( y | x , ω ) .
  • The second component C y e ( z , ω ^ ) reflects the uncertainty induced in the weights ω , also called epistemic uncertainty or model uncertainty. The sources of this uncertainty measured by p ( ω | D ) are for example: (i) non-optimal hyperparameters (the number of hidden layers, the number of units per layer, the type of activation function); (ii) non-optimal training parameters (the minimum learning rate at which the training is stopped, the learning rate decay factor, the batch size used for the mini-batch gradient descent); and (iii) a non-optimal optimization algorithm (ADAGRAD, ADADELTA, ADAM, ADAMAX, NADAM).
or model uncertainty, which refers to the fact that we do not know the model that best explains the given data. For NNs, this is uncertainty from not knowing the best values of weights in all the trainable layers. This is often referred to as reducible uncertainty, because we can theoretically reduce this type of uncertainty by acquiring more data.
Some comments can be made here.
  • The heteroscedastic covariance matrix C y δ ( z , ω ) = diag [ σ j 2 ( z , ω ) ] j = 1 N y can be learned from the data by minimizing the loss functions (26) and (27).
  • The computation of the epistemic covariance matrix C y e ( z , ω ^ ) from Equation (34) requires the knowledge of the Hessian matrix H ( ω ^ ) . In general, this problem is practically insoluble because the matrix H is very huge, and so, is very difficult to compute. However, the problem can be evaded by using the diagonal Hessian approximation
    [ H ( ω ^ ) ] i j = δ i j 2 E ω i 2 ( ω ^ ) ,
    where δ i j is the Kronecker delta. The diagonal matrix elements can be very efficiently computed by using a procedure which is similar to the back-propagation algorithm used for computing the first derivatives [31].
  • The covariance matrix C y ( z , ω ^ ) can be approximated by the conditional average covariance matrix E ( C y | D ) of all network errors ε = y f ( z , ω ^ ) over the data set D , meaning that:
    E ( C y | D ) = 1 N n = 1 N [ ε ( n ) E ( ε ) ] [ ε ( n ) E ( ε ) ] T ,
    where ε ( n ) = y ( n ) f ( z ( n ) , ω ^ ) and E ( ε ) = ( 1 / N ) n = 1 N ε ( n ) . Note that this is a very rough approximation, because in contrast to C y ( z , ω ^ ) , E ( C y | D ) does not depend on z .
  • For exact data, we have C y ( x , ω ^ ) = B C y e ( x , ω ^ ) . Thus, when learning a neural network with exact data, the predictive distribution p ( y | x , D ) is Gaussian with mean f ( x , ω ^ ) and (epistemic) covariance matrix C y e ( x , ω ^ ) .

2.3. Bayesian Networks

Stochastic neural networks are a type of neural networks built by introducing stochastic components into the network (activation or weights). A Bayesian neural network can be regarded as a stochastic neural network trained by using Bayesian inference [32]. This type of neural network provides a natural approach to quantify uncertainty in deep learning and allows to distinguish between epistemic and aleatoric uncertainties.
Bayesian neural networks equipped with variational inference bypass the computation of the evidence p ( D ) = p ( D | ω ) p ( ω ) d ω , which determines the Bayesian posterior p ( ω | D ) via Equation (13), by learning a variational distribution q θ ( ω ) to approximate p ( ω | D ) , that is:
q θ ( ω ) p ( ω | D ) ,
where θ are some variational parameters. These parameters are computed by minimizing the Kullback–Leibler (KL) divergence:
KL ( q θ ( ω ) | p ( ω | D ) ) = q θ ( ω ) log q θ ( ω ) p ( ω | D ) d ω ,
with respect to θ . In fact, the KL divergence is a measure of similarity between the approximate distribution q θ ( ω ) and the posterior distribution obtained from the full Gaussian process p ( ω | D ) , and minimizing the KL divergence to be the same as minimizing the variational free energy defined by
F ( θ , D ) = q θ ( ω ) log p ( D | ω ) d ω + KL ( q θ ( ω ) | p ( ω ) ) = q θ ( ω ) log q θ ( ω ) log p ( ω ) log p ( D | ω ) d ω .
Considering the data model with output noise, assuming C y δ = σ y 2 I , and replacing p ( ω | D ) by q θ ( ω ) in Equation (28), which gives an approximate predictive distribution:
p ( y | x , D ) q θ ( y | x ) = p ( y | x , ω ) q θ ( ω ) d ω
which can be approximated at test time by
q θ ( y | x ) 1 T t = 1 T p ( y | x , ω t ) ,
where ω t is sampled from the distribution q θ ( ω ) . As a result, for p ( y | x , ω ) = N ( y , f ( x , ω ) , C y δ = σ y 2 I ) , the predictive mean and the covariance matrix of the output y given the input x , can be approximated, respectively, by (Appendix B):
E ( y ) 1 T t = 1 T f ( x , ω t ) ,
Cov ( y ) σ y 2 I + 1 T t = 1 T f ( x , ω t ) f ( x , ω t ) T E ( y ) E ( y ) T .
The first term σ y 2 I in Equation (43) corresponds to the homoscedastic uncertainty (the amount of noise in the data), while the second part corresponds to the epistemic uncertainty (how much the model is uncertain in its prediction). The predictive mean (42) is known as model averaging, and is equivalent to performing T stochastic passes through the network and averaging the results. Note that for exact data, i.e., σ y 2 0 , the covariance matrix simplifies to:
Cov ( y ) 1 T t = 1 T f ( x , ω t ) f ( x , ω t ) T E ( y ) E ( y ) T
and describes the epistemic uncertainty.
In this section, we present the most relevant algorithm that is used for Bayesian inference (Bayes-by-backprop), as well as two Bayesian approximate methods. For this purpose, we consider a data model with output noise and assume C y δ = σ y 2 I ; in this case, p ( D | ω ) is given by Equation (15) in conjunction with Equations (21) or (22).

2.3.1. Bayes by Backpropagation

Bayes-by-backprop is a practical implementation of variational inference combined with a reparameterization trick [33]. The idea of the parameterization trick is to introduce a random variable ϵ , and to determine ω by a deterministic transformation t ( ϵ , θ ) , such that ω = t ( ϵ , θ ) follows the distributions q θ ( ω ) . If the variational posterior q θ ( ω ) is a Gaussian distribution with a diagonal covariance matrix, i.e., q θ ( ω ) = N ( ω ; μ ω , diag [ σ ω j 2 ] j = 1 W ) , where σ ω j = [ σ ω ] j and W = dim ( ω ) = dim ( μ ω ) = dim ( σ ω ) , a sample of the weight ω can be obtained by sampling a unit Gaussian ϵ N ( 0 , I ) , scaling it by a standard deviation σ ω and by shifting a mean μ ω . Actually, to guarantee that σ ω is always non-negative, the standard deviation can be parametrized pointwise as σ ω = σ ω ( ρ ω ) = log ( 1 + exp ( ρ ω ) ) or as σ ω = σ ω ( ρ ω ) = exp ( ρ ω / 2 ) (i.e., ρ ω j = log σ ω j 2 ). Thus, the variational posterior parameters are θ = ( μ ω , ρ ω ) . By using a Monte Carlo sampling with one sample, we compute the variational free energy (39) by using the relations:
ϵ N ( 0 , I ) , ω = μ ω + σ ω ( ρ ω ) ϵ , F ( θ , D ) = log q θ ( ω ) log p ( ω ) log p ( D | ω ) ,
where ∘ denotes point-wise multiplication. In Ref. [33], the prior over the weights p ( ω ) is chosen as a mixture of two Gaussian with zero mean but differing variances, meaning that:
p ( ω ) = β N ( ω ; 0 , σ 1 2 I ) + ( 1 β ) N ( ω ; 0 , σ 2 2 I ) ,
while in [34], it is shown that, for q θ ( ω ) = N ( ω ; μ ω , diag [ σ ω j 2 ] ) and p ( ω ) = N ( ω ; 0 , I ) , the KL divergence KL ( q θ ( ω ) | p ( ω ) ) can be analytically computed; the result is:
KL ( q θ ( ω ) | p ( ω ) ) = q θ ( ω ) log q θ ( ω ) log p ( ω ) d ω = 1 2 j = 1 W [ μ ω j 2 + σ ω j 2 log ( σ ω j 2 ) 1 ] ,
where μ ω j = [ μ ω ] j , σ ω j = [ σ ω ] j . Thus:
F ( θ , D ) = log p ( D | ω ) + 1 2 j = 1 W [ μ ω j 2 + σ ω j 2 log ( σ ω j 2 ) 1 ] .
After the training stage, i.e., after the variational posterior parameters θ = ( μ ω , ρ ω ) have been learned, we consider the set of samples { ω t } t = 1 T , where ω t = μ ω + σ ω ( ρ ω ) ϵ t and ϵ t is sampled from the Gaussian distribution N ( 0 , I ) , and compute the predictive mean and covariance matrix according to Equations (42) and (43).

2.3.2. Dropout

Dropout, which was initially proposed as a regularization technique during the training [35,36], is equivalent to a Bayesian approximation. This result was proved in [37] by showing that the variational free energy F ( θ , D ) has the standard form representation of the dropout loss function (as the sum of a square loss function and a L 2 regularization term). A simplified proof of this assertion is given in Appendix C. The term “dropout” refers to removing a unit along with all its connections. The choice of a dropped unit is random. In the case of dropout, the feed-forward operation of a standard neural network (2) and (3) becomes:
u l = W l Z l 1 y l 1 + b l ,
y l = g l ( u l ) ,
where:
Z l 1 = diag [ z k , l 1 ] k = 1 N l 1 ,
z k , l 1 Bernoulli ( p ) .
Essentially, the output of the unit k in layer l 1 , y k , l 1 is multiplied by the binary variable z k , l 1 to create the new output z k , l 1 y k , l 1 . The binary variable z k , l 1 takes the value 1 with probability p and the value 0 with probability 1 p ; thus, if z k , l 1 = 0 , the new output is zero and the unit k is dropped out as an input to the next layer l. The same values of the binary variable are used in the backward pass when propagating the derivatives of the loss function. In the test time, the weights W l are scaled as p W l . Thus, we retain units with probability p at training time, and multiply the weights by p at test time. Alternatively, in order to maintain constant outputs at training time, we can scale the weights by 1 / p , and do not modify the weights at test time. The model parameters ω = { W l , b l } l = 1 L are obtained by minimizing the loss function (21) with possibly an L 2 regularization term.
Uncertainty can be obtained from a dropout neural network [37]. For the set of samples { ω t } t = 1 T , where ω t corresponds to a realization from the Bernoulli distribution Z l 1 ( t ) and W l ( t ) = W l Z l 1 ( t ) for all l = 1 , , L , the predictive mean and covariance matrix are computed by means of Equations (42) and (43).

2.3.3. Batch Normalization

Ioffe and Szegedy [38] introduced batch normalization as a technique for training deep neural networks that normalizes (standardizes) the inputs to a layer for each mini-batch. This allows for higher learning rates, regularizes the model, and reduces the number of training epochs. Moreover, Teye et al. [39] showed that a batch normalized network can be regarded as an approximate Bayesian model, and it can thus be used for modeling uncertainty.
Let us split the data set D into N b mini-batches, and let the nth mini-batch B ( n ) contain M pair samples ( x ( n , m ) , y ( n , m ) ) , meaning that:
D = n = 1 N b B ( n ) , B ( n ) = { ( x ( n , m ) , y ( n , m ) ) } m = 1 M .
In batch normalization, the input u j , l ( n , m ) = i = 1 N l 1 w j i , l y i , l 1 ( n , m ) of the activation function g l corresponding to unit j, layer l, sample m, and mini-batch n, which is first normalized:
u ˜ j , l ( n , m ) = u j , l ( n , m ) μ j , l ( n ) v j , l ( n ) + ε ,
and then scaled and shifted:
u ^ j , l ( n , m ) = α j , l u ˜ j , l ( n , m ) + β j , l ,
where:
μ j , l ( n ) = 1 M m = 1 M u j , l ( n , m ) and v j , l ( n ) = 1 M m = 1 M ( u j , l ( n , m ) μ j , l ( n ) ) 2
are the mean and variance of the activation inputs over the M samples (in unit j, layer l and mini-batch n), respectively, α j , l and β j , l are learnable model parameters, and ε is a small number added to the mini-batch variance to prevent division by zero. By normalization (or whitening), u ˜ j , l ( n , m ) becomes a random variable with zero mean and unit variance, while by scaling and shifting, we guarantee that the transformation (53) can represent the identity transform. In a stochastic framework, we interpret the mean μ j , l ( n ) and variance v j , l ( n ) , corresponding to the nth mini-batch B ( n ) , as realizations of the random variables μ j , l and v j , l , corresponding to the data set D . The model parameters include the learnable parameters:
θ = { w j i , l , α j , l , β j , l | j = 1 , , N l , i = 1 , , N l 1 , l = 1 , , L } ,
and the stochastic parameters:
ω = ( μ , v ) = { ( μ j , l , v j , l ) | j = 1 , , N l l = 1 , , L } ,
where for the mini-batch B ( n ) :
ω ( n ) = ( μ ( n ) , v ( n ) ) = { ( μ j , l ( n ) , v j , l ( n ) ) | j = 1 , , N l l = 1 , , L }
is a realization of ω = ( μ , v ) . Optimizing over mini-batches of size M instead on the full training set, the objective function for the mini-batch B ( n ) becomes [39]:
F ( n ) ( θ ) = 1 2 M m = 1 M | | y ( n , m ) f ω ( n ) ( x ( n , m ) , θ ) | | 2 + Ω ( θ ) ,
where Ω ( θ ) = α l = 1 L W l 2 2 with [ W l ] j i = w j i , l is the regularization term and the notation f ω ( n ) ( x , θ ) indicates that the mean and variance ω ( n ) = ( μ ( n ) , v ( n ) ) are used in the normalization step (52). At the end of the training stage, we obtain:
  • The maximum a posteriori parameters θ ^ = θ MAP ;
  • The mean and variance realizations { ω ( n ) = ( μ ( n ) , v ( n ) ) } n = 1 N b of the stochastic parameters ω = ( μ , v ) ; and
  • The moving averages of the mean and variance over the training set ω ¯ = ( μ ¯ = E ( μ ) , v ¯ = E ( v ) ) .
Some comments can be made here.
  • A batch-normalized network samples the stochastic parameters ω once per training step (mini-batch). For a large number of epochs, the ω ( n ) become independent and identical distributed random variables for each training example;
  • The variational distribution q θ ( ω ) corresponds to the joint distribution of the weights induced by the stochastic parameters ω ;
  • The equivalence between a batch-normalized network and a Bayesian approximation was proven in [39] by showing that (cf. Equations (39) and (55)) KL ( q θ ( ω ) | p ( ω ) ) / θ = ( N / σ y 2 ) Ω ( θ ) / θ . The proof relies on the following simplified assumptions: (i) no scale and shift transformations; (ii) batch normalization applied to each layer; (iii) independent input features in each layer; and (iv) large N and M.
In the inference, the output for a given input x is f ω ¯ ( x , θ ^ ) . To estimate the predictive mean and covariance matrix, we proceed as follows. For each t = 1 , T , we sample a mini-batch B ( t ) from the data set D = n = 1 N b B ( n ) , and for the corresponding mean and variance realization ω ( t ) { ω ( n ) } n = 1 N b , compute f ω ( t ) ( x , θ ^ ) and then the predictive mean and covariance matrix from Equations (42) and (43) with f ω ( t ) ( x , θ ^ ) in place of f ( x , ω t ) .

3. Neural Networks for Atmospheric Remote Sensing

In this section, we design several neural networks for atmospheric remote sensing. To test the neural networks, we considered a specific problem, namely the retrieval of cloud optical thickness and cloud top height from radiances measured by the Earth Polychromatic Imaging Camera (EPIC) onboard the Deep Space Climate Observatory (DSCOVR). DSCOVR is placed in a Lissajous orbit about the Lagrange-1 point, and provides a unique angular perspective in an almost backward direction with scattering angles between 168 and 176 . The EPIC instrument has 10 spectral channels ranging from the UV to the near-IR, which include four channels around the oxygen A- and B-bands; two absorption channels are centered at 688 nm and 764 nm with bandwiths of 0.8 nm and 1.0 nm, respectively, while two continuum bands are centered at 680 nm and 780 nm with bandwiths of 2.0 nm. These four channels are used for inferring the cloud parameters. To generate the database, we use a radiative transfer model based on the discrete ordinate method with matrix exponential [40,41] that uses several acceleration techniques, as for example, the telescoping technique, the method of false discrete ordinate [42], the correlated k-distribution method [43], and the principal component analysis [44,45,46]. The atmospheric parameters to be retrieved are the cloud optical thickness τ and the cloud top height H, while the forward model parameters are the solar and viewing zenith angles and the surface albedo. Specifically, we consider the fact that the cloud optical thickness varies in the range 4.0 , , 16.0 , the cloud top height in the range 2.0 , 10.0 km, the solar and viewing zenith angles in the range 0 , , 60 , and the surface albedo in the range 0.02 , , 0.2 (only snow/ice free scenes are considered). The simulations are performed for a water-cloud model with a Gamma size distribution p ( a ) a α exp ( α a / a mod ) with parameters a mod = 8 μ m and α = 6 . The droplet size ranges between 0.02 and 50.0 μ m, the cloud geometrical thickness is 1 km and the relative azimuth angle between the solar and viewing directions is 176 . The O 2 absorption cross-sections are computed using LBL calculations [47] with optimized rational approximations for the Voigt line profile [48]. The wavenumber grid point spacing is chosen as a fraction (e.g., 1 / 4 ) of the minimum half-width of the Voigt lines taken from HITRAN database [49]. The Rayleigh cross-section and depolarization ratios are computed as in [50], while the pressure and temperature profiles correspond to the US standard model atmosphere [51]. The radiances are solar-flux normalized and are computed by means of the delta-M approximation in conjunction with the TMS correction. We generate N = 20,000 samples by employing the smart sampling technique [52]. Among this data set, 18,000 samples were used for training and the other 2000 for prediction. The noisy spectra are generated by using the measurement noise δ mes N ( 0 , diag [ σ mes j 2 ] j = 1 4 ) , where for the jth channel, we use σ mes j = 0.1 I ¯ j with I ¯ j being the average of the simulated radiance over the N samples.
The neural network algorithms are implemented in FORTRAN by using a feed-forward multilayer perceptron architecture. The tool contains a variety of optimization algorithms, activation functions, regularization terms, dynamic learning rates, and stopping rules. For the present application, the number of hidden layers and the number of units in each layer are optimized by using 2000 samples from the training set for validation. To estimate the performances of different hyperparameter configurations, we used the holdout cross-validation together with a grid search over a set of three values for the number of hidden layers, i.e., { 1 , 2 , 3 } , and a set of eight values for the number of units, i.e., {25, 50, 75, 100, 125, 150, 175, 200}. The mini-batch gradient descent in conjunction with Adaptive Moment Estimation (ADAM) [53] is used as optimization tool, a mini batch of 100 samples is chosen, and a ReLU activation function is considered.

3.1. Neural Networks for Solving the Direct Problem

For a direct problem, the input x is the set of atmospheric parameters, the output y is the set of data, and the forward model F reproduces the radiative transfer model R .
We consider a neural network trained with exact data. For the predictive distribution p ( y | x , D ) given by Equation (32), we assume that the epistemic covariance matrix C y e ( x , ω ^ ) ( = C y ( x , ω ^ ) ) is computed from the statistics of ε = y f ( x , ω ^ ) , that is, C y e E ( C y | D ) , where f ( x , ω ^ ) is the network output. Furthermore, let y δ = y + δ y with δ y N ( 0 , C y δ ) , be the noisy data vector. Using the result:
p ( y δ | y ) exp 1 2 ( y δ y ) T ( C y δ ) 1 ( y δ y ) ,
we compute the predictive distribution for the noisy data p ( y δ | x , D ) by marginalization, that is:
p ( y δ | x , D ) = p ( y δ | y ) p ( y | x , D ) d y = exp 1 2 [ y δ f ( x , ω ^ ) + Δ y ] T [ C y e ( ω ^ ) ] 1 [ y δ f ( x , ω ^ ) + Δ y ] × exp 1 2 Δ y T ( C y δ ) 1 Δ y d Δ y exp 1 2 [ y δ f ( x , ω ^ ) ] T C y 1 ( ω ^ ) [ y δ f ( x , ω ^ ) ] ,
where (cf. Equation (33)) C y ( ω ^ ) = C y e ( ω ^ ) + C y δ and Δ y = y y δ .
In the next step, we solve the inverse problem y δ = f ( x , ω ^ ) by means of a Bayesian approach [54,55]. In this case, the posteriori density p ( x | y δ , D ) is given by
p ( x | y δ , D ) = p ( y δ | x , D ) p ( x ) p ( y ) ;
whence, by assuming that the state vector x is a Gaussian random vector with mean x a and the covariance matrix C x , meaning that:
p ( x ) = N ( x ; x a , C x ) exp 1 2 ( x x a ) C x 1 ( x x a ) T ,
we obtain:
log p ( x | y δ , D ) = 1 2 V ( x | y δ ) + C ,
where:
V ( x | y δ ) = [ y δ f ( x , ω ^ ) ] C y 1 ( ω ^ ) [ y δ f ( x , ω ^ ) ] T + ( x x a ) C x 1 ( x x a ) T
is the a posteriori potential and C is constant. The maximum a posteriori estimator, defined by
x ^ = x MAP = arg max x log p ( x | y δ , D ) = arg min x V ( x | y δ ) ,
can be computed by any optimization method, as for example, the Gauss–Newton method. If the problem is almost linear, the a posteriori density p ( x | y δ , D ) is Gaussian with the mean x ^ and covariance matrix [54]:
C x ( x ^ , ω ^ ) = [ K x T ( x ^ ) C y 1 ( ω ^ ) K x ( x ^ ) + C x 1 ] 1 .
It should be pointed out that we can train the neural network for a data set with output noise D = { ( x ( n ) , y δ ( n ) ) } n = 1 N , in which case, the covariance matrix C y ( ω ^ ) can be directly computed from the statistics of ε = y δ f ( x , ω ^ ) .
Essentially, the method involves the following steps:
  • Train a neural network with exact data for simulating the radiative transfer model;
  • Compute the epistemic covariance matrix from the statistics of all network errors over the data set;
  • Solve the inverse problem by a Bayesian approach under the assumption that the a priori knowledge and the measurement uncertainty are both Gaussian;
  • Determine the uncertainty in the retrieval by assuming that the forward model is nearly linear.
In Figure 1, we plot the histograms of the relative error over the prediction set:
ε x = x pred x x ,
where x stands for τ and H, and x pred and x are the predicted and true values, respectively. Also shown here are the plots of x pred and x pred ± 3 σ x versus x. The results demonstrate that the cloud optical thickness is retrieved with better accuracy than the cloud top height, and that for the cloud top height, the predicted uncertainty are especially unrealistic, because the condition:
x pred 3 σ x x x pred + 3 σ x
is not satisfied. The reason for this failure might be that the forward model is not nearly linear.

3.2. Neural Networks for Solving the Inverse Problem

For an inverse problem, the input x includes the sets of data and forward model parameters (solar and viewing zenith angles, and surface albedo), while the output y includes the set of atmospheric parameters to be retrieved (cloud optical thickness and cloud top height).

3.2.1. Method 1

Let f ( x , ω ^ ) be the output of a neural network trained with exact data, and assume that the predictive distribution p ( y | x , D ) given by Equation (32), and in particular, the epistemic covariance matrix C y e ( x , ω ^ ) ( = C y ( x , ω ^ ) ) are available. For the noisy input z = x + δ x with p ( z | x ) = N ( z ; x , C x δ = σ x 2 I ) and under the assumption that the prior p ( x ) is a slowly varying function as compared to p ( z | x ) , the predictive distribution of the network output can be approximated by [56]
p ( y | z , D ) = p ( y | x , D ) p ( x | z ) d x = 1 p ( z ) p ( y | x , D ) p ( z | x ) p ( x ) d x p ( y | x , D ) p ( z | x ) d x .
Thus, if p ( x ) varies much more slowly than p ( z | x ) = p ( z x ) , we assume that p ( y | z , D ) is the convolution of the predictive distribution for an uncorrupted input p ( y | x , D ) with the input noise distribution p ( z x ) . In the noise-free case, that is, if σ x 0 , we have lim σ x 0 p ( z | x ) = δ ( z x ) , yielding p ( y | z , D ) = p ( y | x , D ) . Using Equation (64), we obtain:
E ( y | z ) = y p ( y | z , D ) d y = y p ( y | x , D ) d y p ( z | x ) d x = f ( x , ω ^ ) p ( z | x ) d x
and:
E ( y y T | z ) = y y T p ( y | z , D ) d y = y y T p ( y | x , D ) d y p ( z | x ) d x = [ C y e ( x , ω ^ ) + f ( x , ω ^ ) f ( x , ω ^ ) T ] p ( z | x ) d x .
Equations (65) and (66) show that in general E ( y | z ) f ( x , ω ^ ) , and a noise process blurs or smooths the original mappings.
To compute the predictive mean E ( y | z ) and the covariance matrix Cov ( y | z ) , the integrals in Equations (65) and (66) must be calculated. This can be done by Monte Carlo integration (by sampling the Gaussian distribution p ( z | x ) ) or by a quadrature method. In the latter case, we consider a uniform grid { x i } i = 1 N int in the interval, say [ z 2 σ x , z + 2 σ x ] , define the weights:
v i = exp 1 2 σ x 2 | | z x i | | 2 i = 1 N int exp 1 2 σ x 2 | | z x i | | 2 ,
and use the computational formulas:
E ( y | z ) = i = 1 N int v i f ( x i , ω ^ ) ,
and:
Cov ( y | z ) = i = 1 N int v i C y e ( x i , ω ^ ) + i = 1 N int v i f ( x i , ω ^ ) f ( x i , ω ^ ) T E ( y | z ) E ( y | z ) T .
The neural network trained with exact data can be deterministic or stochastic, as for example, Bayes-by-backprop, dropout, and batch normalization. In this regard, for each noisy input z , we consider a uniform grid { x i } i = 1 N int around z , calculate the weights v i by means of Equation (67), and for each x i , compute C y e ( x i , ω ^ ) as follows:
  • For a deterministic network, we approximate C y e ( x i , ω ^ ) by the conditional average covariance matrix E ( C y e | D ) of all network errors over, that is, C y e E ( C y e | D ) , yielding i = 1 N int v i C y e ( x i , ω ^ ) = E ( C y e | D ) ;
  • For a Bayes-by-backprop and dropout networks, we compute C y e ( x i , ω ^ ) by means of Equation (44) with E ( y ) as in Equation (42);
  • For a batch normalized network, we compute C y e ( x i , ω ^ ) as
    Cov ( y ) = 1 T t = 1 T f ω ( t ) ( x i , θ ^ ) f ω ( t ) ( x i , θ ^ ) T E ( y ) E ( y ) T
    with E ( y ) = ( 1 / T ) t = 1 T f ω ( t ) ( x i , θ ^ ) .
Note that for a Bayes-by-backprop network, the output is f ( x i , ω ^ = μ ^ ω ) , for a dropout network, the output f ( x i , ω ^ ) is computed without dropout, and for a batch normalized network, the output is f ω ¯ ( x i , θ ^ ) .
In summary, the method uses:
  • Deterministic and stochastic networks trained with exact data to compute the epistemic covariance matrix; and
  • The assumption that the predictive distribution of the network output is the convolution of the predictive distribution for an uncorrupted input with the input noise distribution to estimate the covariance matrix.
Under the assumption that the noise process is Gaussian, the convolution integrals are computed by a quadrature method with a uniform grid around the noisy input.
The results for Method 1 with deterministic and stochastic networks are illustrated in Figure 2, Figure 3, Figure 4 and Figure 5.

3.2.2. Method 2

In order to model both heteroscedastic and epistemic uncertainties, we used the approach described in Ref. [57]. More precisely, we considered the data model with input and output noise, and used dropout to learn the heteroscedastic covariance matrix C y δ ( z , ω ) = diag [ σ j 2 ( z , ω ) ] j = 1 N y from the data (see Section 2.1). The network has the output [ μ j ( z , ω ) , σ j 2 ( z , ω ) ] R 2 N y and is trained to predict the log variance ρ j = log σ j 2 , in which case, the likelihood loss function is given by Equations (26) and (27). Considering the set of samples { ω t } t = 1 T , where ω t corresponds to a realization of the Bernoulli distribution Z l 1 ( t ) and W l ( t ) = W l Z l 1 ( t ) for all l = 1 , , L , we compute the predictive mean and covariance matrix as [57]
E ( y ) 1 T t = 1 T μ ( z , ω t ) ,
Cov ( y ) 1 T t = 1 T diag [ σ j 2 ( z , ω t ) ] j = 1 N y + 1 T t = 1 T μ ( z , ω t ) μ ( z , ω t ) T E ( y ) E ( y ) T ,
for each noisy input z . The first term in Equation (71) reproduces the heteroscedastic uncertainty, while the second and third terms reproduce the epistemic uncertainty.
Instead of a dropout network, a Bayes-by-backprop network can also be used to learn the heteroscedastic covariance matrix from the data. In this case, F ( θ , D ) is given by Equation (47) with log p ( D | ω ) = E D ( ω ) = n = 1 N E D ( n ) ( ω ) and E D ( n ) ( ω ) as in Equation (27).
The results for Method 2 with a dropout and a Bayes-by-backprop network are illustrated in Figure 6 and Figure 7.

3.2.3. Method 3

Let ω ^ = { W l , b l } l = 1 L be the parameters of a dropout network trained with exact data. Further, assume that the input data are noisy, i.e., z = x + δ x with p ( z | x ) = N ( z ; x , C x δ = diag [ σ x k 2 ] k = 1 N x ) . In order to compute the heteroscedastic uncertainty, we forward propagate the input noise through the network. This is done by using assumed density filtering and interval arithmetic.
  • Assumed density filtering (ADF). This approach was applied to neural networks by Gast and Roth [58] to replace each network activation by probability distributions. In the following, we provide a simplified justification of this approach, while for a more pertinent analysis, we refer to Appendix D. For a linear layer ( g l ( x ) = x ) , the feed-forward operation (without dropout) is:
    y k , l = j = 1 N l 1 w k j , l y j , l 1 + b k , l , k = 1 , , N l ,
    By straightforward calculation, we find:
    μ k , l = E ( y k , l ) = j = 1 N l 1 w k j , l E ( y j , l 1 ) + b k , l = j = 1 N l 1 w k j , l μ j , l 1 + b k , l ,
    and:
    E ( y k , l ) E ( y k 1 , l ) = j = 1 N l 1 j 1 = 1 N l 1 w k j , l w k 1 j 1 , l μ j , l 1 μ j 1 , l 1 + b k , l μ k 1 , l + b k 1 , l μ k , l b k , l b k 1 , l ,
    yielding:
    E ( y k , l y k 1 , l ) E ( y k , l ) E ( y k 1 , l ) = j = 1 N l 1 j 1 = 1 N l 1 w k j , l w k 1 j 1 , l [ E ( y j , l 1 y j 1 , l 1 ) E ( y j , l 1 ) E ( y j 1 , l 1 ) ] .
    Assuming that the y k , l are independent random variables, in which case, the covariance matrix corresponding to the column vector [ y 1 . l , , y N l l ] T is diagonal, meaning that:
    E ( y k , l y k 1 , l ) E ( y k , l ) E ( y k 1 , l ) = δ k k 1 E ( y k , l 2 ) E 2 ( y k , l ) = δ k k 1 v k , l ,
    we obtain v k , l = j = 1 N l 1 w k j , l 2 v j , l 1 . In summary, the iterative process for a linear layer is:
    μ k , 0 = x k , v k , 0 = σ x k 2 .
    μ k , l = j = 1 N l 1 w k j , l μ j , l 1 + b k , l ,
    v k , l = j = 1 N l 1 w k j , l 2 v j , l 1 , l = 1 , , L .
    For a ReLU activation function ReLU ( x ) = max ( 0 , x ) , it was shown that:
    μ ReLU k , l = v k , l ϕ ( α ) + μ k , l Φ ( α ) ,
    v ReLU k , l = μ k , l v k , l ϕ ( α ) + ( μ k , l 2 + v k , l ) Φ ( α ) μ ReLU k , l 2 ,
    where μ k , l and v k , l are given by Equations (78) and (79), respectively, and:
    α = μ k , l v k , l , ϕ ( x ) = 1 2 π exp 1 2 x 2 , Φ ( x ) = x ϕ ( y ) d y .
  • Interval Arithmetic (IA). Interval arithmetic is based on an extension of the real number system to a system of closed intervals on the real axis [59]. For the intervals X and Y, the elementary arithmetic operations are defined by the rule X Y = { x y | x X , y Y } , where the binary operation ⊕ can stand for addition, subtraction, multiplication, or division. This definition guarantees that x y X Y . Functions of interval arguments are defined in terms of standard set mapping, that is, the image of an interval X under a function f is the set f ( X ) = { f ( x ) | x X } . This is not the same as an interval function obtained from a real function f by replacing the real argument by an interval argument and the real arithmetic operations by the corresponding interval operations. The latter is called an interval extension of the real function f and is denoted by F ( X ) . As a corollary of the fundamental theorem of interval analysis, it can be shown that f ( X ) F ( X ) . Interval analysis provides a simple and accessible way to assess error propagation. The iterative process for error propagation is (compared with Equations (77)–(79)):
    Y k , 0 = [ x k σ x k , x k + σ x k ] ,
    U k , l = i = 1 N l 1 w k j , l Y j , l 1 + b k , l ,
    Y k , l = G l ( U k , l ) , l = 1 , , L ,
    while the output predictions μ k , L and their standard deviations σ k , L = v k , L are computed as
    μ k , L = 1 2 [ Y ̲ k , L + Y ¯ k , L ] ,
    σ k , L = 1 2 [ Y ¯ k , L Y ̲ k , L ] , k = 1 , , N L ,
    where G l ( U ) is the interval extension of the activation function g l , and X ̲ and X ¯ are the left and right endpoints of the interval X, respectively, that is, X = [ X ̲ , X ¯ ] .
By assumed density filtering and interval arithmetic, the forward pass of a neural network generates not only the output predictions μ L = [ μ 1 , L , , μ N L , L ] T but also their variances v L = [ v 1 , L , , v N L , L ] T . Following Ref. [2], we consider now a network with dropout, that is, in Equations (78), (79) and (83), we replace w k j , l by w k j , l z j , l 1 , where z j , l 1 Bernoulli ( p ) and the dropout probability p is the same as that used for training. For the set of samples { ω t } t = 1 T , where ω t corresponds to a realization from the Bernoulli distribution Z l 1 ( t ) and W l ( t ) = W l Z l 1 ( t ) for all l = 1 , , L , we denote by μ L ( x , ω t ) and v L ( x , ω t ) the outputs of the network for an input x corrupted by the noise δ x N ( 0 , C x δ = diag [ σ x k 2 ] k = 1 N x ) , and compute the predictive mean and covariance matrix as (Appendix B)
E ( y | x ) 1 T t = 1 T μ L ( x , ω t ) ,
Cov ( y | x ) 1 T t = 1 T diag [ v k , L ( x , ω t ) ] k = 1 N y + 1 T t = 1 T μ L ( x , ω t ) μ L ( x , ω t ) T E ( y | x ) E ( y | x ) T .
From Equations (87) and (88), it is obvious that the prediction ensemble is not generated by the output of the network f ( x , ω t ) , but by the prediction μ L ( x , ω t ) . In summary, the algorithm involves the following steps:
  • Transform a dropout network into its assumed density filtering or interval arithmetic versions (which does not require retraining);
  • Propagate ( x , v 0 ) through the dropout network and collect T output predictions and variances;
  • Compute the predictive mean and covariance matrix by means of Equations (87) and (88).
The results for Method 3 using assumed density filtering and interval arithmetic are illustrated in Figure 8 and Figure 9, respectively. The main drawback of this method is that it requires the knowledge of the exact input data x . Because in atmospheric remote sensing that is not the case, Method 3 will be used for a comparison with the other two methods.

3.3. Summary of Numerical Analysis

In Table 1, we summarize the results of our numerical simulations by illustrating the average relative error and the standard deviation over the prediction set, E ( ε x ) ± E ( [ ε x E ( ε x ) ] 2 ) and E ( σ x ) , respectively, where x stands for the cloud optical thickness τ and the cloud top height H. The accuracy of a method is reflected by the bias of the error E ( ε x ) and the interval about the mean with length E ( [ ε x E ( ε x ) ] 2 ) , while the precision is reflected by the standard deviation E ( σ x ) (which determines the length of the uncertainty interval). Note that (i) E ( [ ε x E ( ε x ) ] 2 ) reproduces the square root of the diagonal elements of the conditional average covariance matrix E ( C y | D test ) of all network errors over the test set D test ; and (ii) roughly speaking, the epistemic (model) uncertainties are large if there are large variations around the mean.
The results in Figure 1, Figure 2, Figure 3, Figure 4, Figure 5, Figure 6, Figure 7, Figure 8 and Figure 9 and Table 1 can be summarized as follows:
  • For the direct problem, the neural network method used in conjunction with a Bayesian inversion method provides satisfactory accuracy, but does not correctly predict the uncertainty. The reason for this is that the forward model is not nearly linear, which is the main assumption for computing the uncertainty in the retrieval.
  • For the inverse problem, the following features are apparent.
    (a)
    Method 1 using a deterministic and Bayes-by-backprop network yields a low accuracy, while the method using dropout and batch-normalized networks provides high accuracy;
    (b)
    Method 2 using a dropout network has an acceptable accuracy. For cloud top height retrieval, the method using a Bayes-by-backprop network has a similar accuracy, but for cloud optical thickness retrieval, the accuracy is low. Possible reasons for the loss of accuracy of a Bayes-by-backprop network can be (i) a non-optimal training and/or the use of the prior p ( ω ) = N ( ω ; 0 , I ) instead of that given by Equation (45) (recall that for p ( ω ) = N ( ω ; 0 , I ) , the KL divergence KL ( q θ ( ω ) | p ( ω ) ) can be computed analytically).
    (c)
    Method 3 with assumed density filtering and interval arithmetic is comparable with Method 2 using a dropout network, that is, the method has a high accuracy.
    (d)
    From the methods with reasonable accuracy, by using a dropout network, Method 2 predicts similar uncertainties to Method 3 with assumed density filtering and interval arithmetic. In contrast, using dropout and batch normalized networks, Method 1 provides lower uncertainties, dominated by epistemic uncertainties. In general, it seems that Method 1 predicts a too small heteroscedastic uncertainty. Possibly, this is due to the fact that we trained a neural network with exact data, and for this retrieval example, the predictive distribution, represented as a convolution integral over the input noise distribution, does not correctly reproduce the aleatoric uncertainty.
    (e)
    Because E ( [ ε x E ( ε x ) ] 2 ) < E ( σ x ) , we deduce that the conditional average covariance matrix E ( C y | D test ) does not coincide with the predictive covariance matrix Cov ( y ) , which reflects the uncertainty.

4. Conclusions

We presented several neural networks for predicting uncertainty in an atmospheric remote sensing. The neural networks are designed in a Bayesian framework and are devoted to the solution of direct and inverse problems.
  • For solving the direct problem, we considered a neural network for simulating the radiative transfer model, computed of the epistemic covariance matrix from the statistics of all network errors over the data set, solving the inverse problem by a Bayesian approach, and determined the uncertainty in the retrieval by assuming that the forward model is nearly linear.
  • For solving the inverse problem, two neural network methods, relying on different assumptions, were implemented:
    (a)
    The first method uses deterministic and stochastic (Bayes-by-backprop, dropout, and batch normalization) networks to compute the epistemic covariance matrix and under the assumption that the predictive distribution of the network output is the convolution of the predictive distribution for a noise-free input with the input noise distribution, estimates the covariance matrix;
    (b)
    The second method uses dropout and Bayes-by-backprop to learn the heteroscedastic covariance matrix from the data.
In addition, for solving the inverse problem, a third method that uses a dropout network and forward propagates the input noise through the network by using assumed density filtering and interval arithmetic was designed. Because this method requires the knowledge of the exact input data, it was used only for testing purposes.
Our numerical analysis has shown that a dropout network that is used to learn the heteroscedastic covariance matrix from the data is appropriate for predicting the uncertainty associated with the retrieval of cloud parameters from EPIC measurements. In fact, the strengths of a dropout network are (i) its capability to avoid overfitting and (ii) its stochastic character (the method is equivalent to a Bayesian approximation).
All neural network algorithms are implemented in FORTRAN and incorporated in a common tool. In the future, we intend to implement the algorithms in the high-level programming language Python and use the deep learning library PyTorch. This software library has a variety of network architectures that provide auto-differentiation and support GPUs to enable fast and efficient computation. The Python tool will be released through a public repository to make the methods available to the scientific community.

Author Contributions

Conceptualization, A.D. (Adrian Doicu) and A.D. (Alexandru Doicu); software, A.D. (Adrian Doicu), A.D. (Alexandru Doicu) and D.S.E.; formal analysis, D.L. and T.T.; writing—original draft preparation, A.D. (Adrian Doicu) and A.D. (Alexandru Doicu); writing—review and editing, D.S.E., D.L. and T.T. All authors have read and agreed to the published version of the manuscript.

Funding

This research received no external funding.

Conflicts of Interest

The authors declare no conflict of interest.

Appendix A

According to the Laplace approximation, we expand the loss function:
E ( ω ) = 1 2 n = 1 N [ y ( n ) f ( z ( n ) , ω ) ] T [ C y δ ( z ( n ) , ω ) ] 1 [ y ( n ) f ( z ( n ) , ω ) ] + 1 2 ω T C ω 1 ω ,
around the point estimate ω ^ = ω MAP and use the optimality condition E ( ω ^ ) = 0 , to obtain:
E ( ω ) = E ( ω ^ ) + 1 2 Δ ω T H ( ω ^ ) Δ ω ,
Δ ω = ω ω ^ ,
[ H ( ω ^ ) ] i j = E ω i ω j ( ω ^ ) ,
and further (cf. Equation (13)):
p ( ω | D ) exp [ E ( ω ) ] exp 1 2 ( ω ω ^ ) T H ( ω ^ ) ( ω ω ^ ) .
Inserting Equations (12) and (A5) in the expression of the predictive distribution as given by Equation (28) we find:
p ( y | z , D ) = p ( y | z , ω ) p ( ω | D ) d ω exp 1 2 [ y f ( z , ω ) ] T [ C y δ ( z , ω ) ] 1 [ y f ( z , ω ) ] × exp 1 2 ( ω ω ^ ) T H ( ω ^ ) ( ω ω ^ ) d ω .
In Equation (A6), the model function f ( x , ω ) can be approximated by a linear Taylor expansion around ω ^ , that is:
f ( x , ω ) f ( x , ω ^ ) + K ω ( x , ω ^ ) ( ω ω ^ ) ,
where:
K ω ( x , ω ) = f ω ( x , ω )
is the Jacobian of f with respect to ω . Substituting Equation (A7) into Equation (A6), approximating C y δ ( z , ω ) C y δ ( z , ω ^ ) , and computing the integral over ω gives the representation (32) for the predictive distribution p ( y | z , D ) .

Appendix B

For p ( y | x , ω ) = N ( y , f ( x , ω ) , C y δ = σ y 2 I ) , we compute the predictive mean of the output y given the input x , as
E ( y ) = y p ( y | x , D ) d y y q θ ( y | x ) d y = y p ( y | x , ω ) q θ ( ω ) d ω d y = y N ( y ; f ( x , ω ) , σ y 2 I ) d y q θ ( ω ) d ω = f ( x , ω ) q θ ( ω ) d ω 1 T t = 1 T f ( x , ω t ) ,
and by using the result:
E ( y y T ) = y y T p ( y | x , D ) d y y y T q θ ( y | x ) d y = y y T p ( y | x , ω ) q θ ( ω ) d ω d y = y y T N ( y ; f ( x , ω ) , σ y 2 I ) d y q θ ( ω ) d ω = [ σ y 2 I + f ( x , ω ) f ( x , ω ) T ] q θ ( ω ) d ω σ y 2 I + 1 T t = 1 T f ( x , ω t ) f ( x , ω t ) T ,
the covariance matrix as
Cov ( y ) σ y 2 I + 1 T t = 1 T f ( x , ω t ) f ( x , ω t ) T E ( y ) E ( y ) T .
The predictive mean and covariance matrix of a dropout network with assumed density filtering given by Equations (87) and (88), respectively, can be computed in the same manner by taking into account that in this case, the predictive power of the network is given by
p ( y | x , ω ) = N ( y ; μ L ( x , ω ) , diag [ v j , L ( x , ω ) ] j = 1 N y ) ,
where μ L = [ μ 1 . L , , μ N L L ] T and v L = [ v 1 . L , , v N L L ] T are the output predictions and their variances.

Appendix C

In this appendix, which is borrowed from [37], we show that the variational free energy has the standard form representation of the dropout loss function (as the sum of a square loss function and an L 2 regularization term).
For ω = { W l , b l } l = 1 L and W l = [ w k , l ] k = 1 N l 1 R N l × N l 1 , we construct the variational distribution q θ ( ω ) as
q θ ( ω ) = l = 1 L q θ ( W l , b l ) = l = 1 L q ( W l ) l = 1 L q ( b l ) ,
with:
q ( W l ) = k = 1 N l 1 q ( w k , l ) ,
q ( w k , l ) = p N ( m k , l , σ 2 I N l ) + ( 1 p ) N ( 0 , σ 2 I N l ) ,
and:
q ( b l ) = N ( n l , σ 2 I N l ) .
Here, p [ 0 , 1 ] is an activation probability, σ > 0 a scalar, and M l = [ m k , l ] k = 1 N l 1 and n l are a variational parameter to be determined; thus, θ = { M l , n l } l = 1 L . The key point of the derivation is the representation of q ( w k , l ) as a mixture of two Gaussians with the same variance (cf. Equation (A11)). When the standard deviation σ tends towards 0, the Gaussians tend to Dirac delta distributions showing that when sampling from the mixture of the two Gaussians is equivalent to sampling from a Bernoulli distribution that returns either the value 0 with probability 1 p or m k , l with probability p , that is:
q ( w k , l ) = m k , l 0 with probability p with probability 1 p .
As a result, we obtain:
W l = [ w 1 , l , w 2 , l , , w N l 1 , l ] = [ m 1 , l , m 2 , l , , m N l 1 , l ] z 1 , l 1 0 0 0 z 2 , l 1 0 0 0 z N l 1 , l 1 = M l Z l 1 ,
where Z l 1 = diag [ z k , l 1 ] k = 1 N l 1 with z k , l 1 Bernoulli ( p ) . Note that the binary variable z k , l 1 corresponds to the unit k in layer l 1 being dropped out as an input to layer l. For b l , we take into account that q ( b l ) = lim σ 0 N ( n l , σ 2 I N l ) = δ ( b l n l ) ; hence, in the limit σ 0 , b l is approximately deterministic, and we have b l n l .
The variational parameters M l and n l are computed by minimizing the variational free energy (39), that is:
F ( θ , D ) = q θ ( ω ) log p ( D | ω ) d ω + KL ( q θ ( ω ) | p ( ω ) ) .
The two terms in the above equation are computed as follows.
  • For the first term, written as
    q θ ( ω ) log p ( D | ω ) d ω = n = 1 N q θ ( ω ) log p ( y ( n ) | x ( n ) , ω ) d ω ,
    we use the following reparameterization trick. Let ϵ q ( ϵ ) be an auxiliary variable representing the stochasticity during the training, such that ω = t ( ϵ , θ ) for some function t. Assuming that q θ ( ω | ϵ ) = δ ( ω t ( ϵ , θ ) ) , we find:
    q θ ( ω ) = q θ ( ω | ϵ ) q ( ϵ ) d ϵ = δ ( ω t ( ϵ , θ ) ) q ( ϵ ) d ϵ ,
    implying:
    q θ ( ω ) log p ( y ( n ) | x ( n ) , ω ) d ω = δ ( ω t ( ϵ , θ ) ) q ( ϵ ) log p ( y ( n ) | x ( n ) , ω ) d ω d ϵ = q ( ϵ ) log p ( y ( n ) | x ( n ) , t ( ϵ , θ ) ) d ϵ .
    Computing the above integral by a Monte Carlo approach with a single sample ϵ ^ q ( ϵ ) , yields:
    q θ ( ω ) log p ( y ( n ) | x ( n ) , ω ) d ω = log p ( y ( n ) | x ( n ) , ω ^ ) ,
    where ω ^ = t ( ϵ ^ , θ ) . In our case and in view of Equations (A10)–(A12), we reparametrize the integrands by setting:
    W l = ( M l + σ S l ) Z l 1 + σ S l ( I N l 1 Z l 1 ) ,
    b l = n l + σ ϵ l ,
    where:
    S l = [ s k , l ] k = 1 N l 1 , s k , l N ( 0 , I N l ) ,
    Z l 1 = diag [ z k , l 1 ] k = 1 N l 1 , z k , l 1 Bernoulli ( p ) ,
    ϵ l N ( 0 , I N l ) ,
    to obtain:
    q θ ( ω ) log p ( D | ω ) d ω = n = 1 N log p ( y ( n ) | x ( n ) , ω ^ n ) ,
    where:
    ω ^ n = { W ^ l ( n ) , b ^ l ( n ) } l = 1 L ,
    W ^ l ( n ) = ( M l + σ S ^ l ( n ) ) Z ^ l 1 ( n ) + σ S ^ l ( n ) ( I N l 1 Z ^ l 1 ( n ) ) ,
    b ^ l ( n ) = n l + σ ϵ ^ l ( n ) ,
    for the realizations S ^ l ( n ) , Z ^ l 1 ( n ) , and ϵ ^ l ( n ) given by Equations (A20)–(A22). Taking the limit σ 0 , we find that the realizations W ^ l ( n ) and b ^ l ( n ) can be approximated as
    W ^ l ( n ) M l Z ^ l 1 ( n ) , b ^ l ( n ) n l .
  • In the case of W l , the second term in Equation (A13) is the KL divergence between a mixture of Gaussians and a single Gaussian, that is:
    KL ( q ( W l ) | p ( W l ) ) = q ( W l ) log q ( W l ) p ( W l ) d W l ,
    where q ( W l ) is as in Equations (A10) and (A11) and p ( W l ) = k = 1 N l 1 p ( w k , l ) with p ( w k , l ) = N ( 0 , I N l ) . This term can be evaluated by using the following result: For K , N N , let p = ( p 1 , , p K ) be a probability vector, q ( x ) = k = 1 K p k N ( x ; μ k , Σ k ) with x R N a mixture of Gaussians with K components, and p ( x ) = N ( x ; 0 , I N ) . Then, for sufficiently large N, we have the approximation:
    KL ( q ( x ) | p ( x ) ) k = 1 K μ k T μ k + tr ( Σ k ) N ( 1 + log 2 π ) log | Σ k | .
    Consequently, for large numbers of hidden units N l , l = 1 , , L , we find:
    KL ( q ( W l ) | p ( W l ) ) N l N l 1 ( σ 2 log ( σ 2 ) 1 ) + p 2 k = 1 N l 1 m k , l T m k , l + C ,
    where C is a constant. In the case of b l , the KL divergence KL ( q ( b l ) | p ( b l ) ) , where (cf. Equation (A12)) q ( b l ) = N ( n l , σ 2 I N l ) and p ( b l ) = N ( 0 , I N l ) , is a mixture of two single Gaussian and can be analytically computed as
    KL ( q ( b l ) | p ( b l ) ) = 1 2 [ n l T n l + N l ( σ 2 log ( σ 2 ) 1 ) ] + C .
Collecting all the results, we obtain:
F ( θ , D ) = n = 1 N log p ( y ( n ) | x ( n ) , ω ^ n ) + 1 2 l = 1 L p M l 2 2 + n l 2 2 .
Thus, the variational free energy F ( θ , D ) has the standard form representation as the sum of a square loss function and an L 2 regularization term.

Appendix D

In this appendix, we describe the uncertainty propagation based on assumed density filtering by following the analysis given in [58].
The feed-forward operation of a neural network can be written as (cf. Equations (2) and (3)):
y l = f l ( y l 1 ; ω l ) = g l ( W l y l 1 + b l ) ,
where ω l = { W l , b l } . Thus, each layer function f l ( y l 1 ; ω l ) is a nonlinear transformation of the previous activation y l 1 parametrized by ω l . The deep neural network can then be expressed as a succession of nonlinear layers:
f ( x , ω ) = f L ( f L 1 ( f 1 ( y 0 ; ω 1 ) ) .
To formalize the deep probabilistic model, we replace each activation, including input and output, by probability distributions. In particular, we assume that the joint density of all activations is given by
p ( y 0 : L ) = p ( y 0 ) l = 1 L p ( y l | y l 1 ) , p ( y l | y l 1 ) = δ ( y l f l ( y l 1 ; ω l ) ) , p ( y 0 ) = k = 1 N l N ( y k , 0 ; x k , σ x k 2 ) ,
where p ( y 0 : L ) = p ( y 0 , , y L ) and C x δ = diag [ σ x k 2 ] k = 1 N x . Because this distribution is intractable, we apply the assumed density filtering approach to the network activations. The goal of this approach is to find a tractable approximation q ( y 0 : L ) of p ( y 0 : L ) , that is:
p ( y 0 : L ) q ( y 0 : L ) ,
where:
q ( y 0 : L ) = q ( y 0 ) l = 1 L q ( y l ) , q ( y l ) = k = 1 N l N ( y k , l ; μ k , l , v k , l ) , q ( y 0 ) = p ( y 0 ) ,
and μ k , l and v k , l are the mean and variance of the activation of unit k in layer l, respectively. Thus, starting from input activation q ( y 0 ) = p ( y 0 ) , we approximate subsequent layer activations by independent Gaussian distributions q ( y l ) . To compute the approximant q ( y 0 : L ) , we use an iterative process (layer by layer) initialized by q ( y 0 ) = p ( y 0 ) . In particular, for a layer l 1 , we assume that q ( y 0 ) , , q ( y l 1 ) are known, or equivalently, that { ( μ l 1 , v l 1 ) } l 1 = 0 l 1 are known, and aim to compute ( μ l , v l ) , where μ k , l = [ μ l ] k and v k , l = [ v l ] k . For this purpose, we take into account that the layer function f l transforms the activation y l 1 into the distribution:
p ^ ( y 0 : l ) = p ( y l | y l 1 ) q ( y 0 : l 1 ) = p ( y l | y l 1 ) l 1 = 0 l 1 q ( y l 1 ) ,
where p ( y l | y l 1 ) = δ ( y l f l ( y l 1 ; ω l ) ) is the true posterior at layer l and q ( y 0 : l 1 ) = l 1 = 0 l 1 q ( y l 1 ) the previous approximating factor. Furthermore, we compute the first and second-order moments of p ^ ( y 0 : l ) . This will be done in two steps. In the first step, we derive the moments of an activation variable y k that belongs to all layers excluding the last layer, i.e., y k is an element of y 0 : l 1 = { y 0 , , y l 1 } , while in the second step, we assume that y k is an activation variable contained in the last layer y l = { y 1 , l , , y N l , l } . Thus,
  • For y k y 0 : l 1 , we use the relations:
    p ^ ( y 0 : l ) = δ ( y l f l ( y l 1 ; ω l ) ) q ( y k ) q ( y k ¯ , 0 : l 1 ) ,
    and δ ( x x 0 ) d x = 1 (yielding ( y l f l ( y l 1 ; ω l ) ) d y l = 1 ) to obtain:
    E p ^ ( y k ) = p ^ ( y 0 : l ) y k d y = q ( y k ) y k d y k = E q ( y k ) ( y k ) ,
    where q ( y k ¯ , 0 : l 1 ( 0 : l 1 ) ) corresponds to the density of all variables excluding y k , and d y = l 1 = 0 l d y l 1 , while:
  • For y k y l , we use the relations:
    p ^ ( y 0 : l ) = δ ( y k ¯ , l f k ¯ , l ( y l 1 ; ω l ) ) δ ( y k f k , l ( y l 1 ; ω l ) ) q ( y 0 : l 1 ) ,
    and x δ ( x x 0 ) d x = x 0 (yielding δ ( y k f k , l ( y l 1 ; ω l ) ) d y k = f k , l ( y l 1 ; ω l ) ), to obtain:
    E p ^ ( y k ) = p ^ ( y 0 : l ) y k d y = q ( y l 1 ) f k , l ( y l 1 ; ω l ) d y l 1 = E q ( y l 1 ) ( f k , l ( y l 1 ; ω l ) ) .
Replacing y k with y k 2 and repeating the above arguments, we find
E p ^ ( y k 2 ) = E q ( y k ) ( y k 2 ) , y k y 0 : l 1 ,
E p ^ ( y k 2 ) = E q ( y l 1 ) ( f k , l 2 ( y l 1 ; ω l ) ) , y k y l .
From Equations (A31) and (A33), we see that for all layers except for the lth layer, the moments remain unchanged after the update. The moments for the lth layer will be computed by means of Equations (A32) and (A34). For a linear activation function, i.e., f k , l ( y l 1 ; ω l ) = i = 1 N l 1 w k i , l 1 y i , l 1 + b k , l , we find that the expressions of μ k , l and v k , l are given by Equations (80) and (81), respectively, while for a ReLU activation function ReLU ( x ) = max ( 0 , x ) , these are given by Equations (80) and (81), respectively.

References

  1. Tibshirani, R. A Comparison of Some Error Estimates for Neural Network Models. Neural Comput. 1996, 8, 152–163. [Google Scholar] [CrossRef]
  2. Loquercio, A.; Segu, M.; Scaramuzza, D. A General Framework for Uncertainty Estimation in Deep Learning. IEEE Robot. Autom. Lett. 2020, 5, 3153–3160. [Google Scholar] [CrossRef] [Green Version]
  3. Oh, S.; Byun, J. Bayesian Uncertainty Estimation for Deep Learning Inversion of Electromagnetic Data. IEEE Geosci. Remote Sens. Lett. 2021, 1–5. [Google Scholar] [CrossRef]
  4. Chevallier, F.; Chéruy, F.; Scott, N.A.; Chédin, A. A Neural Network Approach for a Fast and Accurate Computation of a Longwave Radiative Budget. J. Appl. Meteorol. 1998, 37, 1385–1397. [Google Scholar] [CrossRef]
  5. Chevallier, F.; Morcrette, J.J.; Chéruy, F.; Scott, N.A. Use of a neural-network-based long-wave radiative-transfer scheme in the ECMWF atmospheric model. Q. J. R. Meteorol. Soc. 2000, 126, 761–776. [Google Scholar] [CrossRef]
  6. Cornford, D.; Nabney, I.T.; Ramage, G. Improved neural network scatterometer forward models. J. Geophys. Res. Ocean. 2001, 106, 22331–22338. [Google Scholar] [CrossRef] [Green Version]
  7. Krasnopolsky, V.M. The Application of Neural Networks in the Earth System Sciences; Springer: Dordrecht, The Netherlands, 2013. [Google Scholar] [CrossRef]
  8. Efremenko, D.S. Discrete Ordinate Radiative Transfer Model With the Neural Network Based Eigenvalue Solver: Proof of Concept. Light Eng. 2021, 1, 56–62. [Google Scholar] [CrossRef]
  9. Fan, Y.; Li, W.; Gatebe, C.K.; Jamet, C.; Zibordi, G.; Schroeder, T.; Stamnes, K. Atmospheric correction over coastal waters using multilayer neural networks. Remote Sens. Environ. 2017, 199, 218–240. [Google Scholar] [CrossRef]
  10. Fan, C.; Fu, G.; Noia, A.D.; Smit, M.; Rietjens, J.H.; Ferrare, R.A.; Burton, S.; Li, Z.; Hasekamp, O.P. Use of A Neural Network-Based Ocean Body Radiative Transfer Model for Aerosol Retrievals from Multi-Angle Polarimetric Measurements. Remote Sens. 2019, 11, 2877. [Google Scholar] [CrossRef] [Green Version]
  11. Gao, M.; Franz, B.A.; Knobelspiesse, K.; Zhai, P.W.; Martins, V.; Burton, S.; Cairns, B.; Ferrare, R.; Gales, J.; Hasekamp, O.; et al. Efficient multi-angle polarimetric inversion of aerosols and ocean color powered by a deep neural network forward model. Atmos. Meas. Tech. 2021, 14, 4083–4110. [Google Scholar] [CrossRef]
  12. Shi, C.; Hashimoto, M.; Shiomi, K.; Nakajima, T. Development of an Algorithm to Retrieve Aerosol Optical Properties Over Water Using an Artificial Neural Network Radiative Transfer Scheme: First Result From GOSAT-2/CAI-2. IEEE Trans. Geosci. Remote Sens. 2021, 59, 9861–9872. [Google Scholar] [CrossRef]
  13. Jiménez, C.; Eriksson, P.; Murtagh, D. Inversion of Odin limb sounding submillimeter observations by a neural network technique. Radio Sci. 2003, 38, 27-1–27-8. [Google Scholar] [CrossRef]
  14. Holl, G.; Eliasson, S.; Mendrok, J.; Buehler, S.A. SPARE-ICE: Synergistic ice water path from passive operational sensors. J. Geophys. Res. Atmos. 2014, 119, 1504–1523. [Google Scholar] [CrossRef]
  15. Strandgren, J.; Bugliaro, L.; Sehnke, F.; Schröder, L. Cirrus cloud retrieval with MSG/SEVIRI using artificial neural networks. Atmos. Meas. Tech. 2017, 10, 3547–3573. [Google Scholar] [CrossRef] [Green Version]
  16. Efremenko, D.S.; Loyola R, D.G.; Hedelt, P.; Spurr, R.J.D. Volcanic SO2 plume height retrieval from UV sensors using a full-physics inverse learning machine algorithm. Int. J. Remote Sens. 2017, 38, 1–27. [Google Scholar] [CrossRef] [Green Version]
  17. Wang, D.; Prigent, C.; Aires, F.; Jimenez, C. A Statistical Retrieval of Cloud Parameters for the Millimeter Wave Ice Cloud Imager on Board MetOp-SG. IEEE Access 2017, 5, 4057–4076. [Google Scholar] [CrossRef]
  18. Brath, M.; Fox, S.; Eriksson, P.; Harlow, R.C.; Burgdorf, M.; Buehler, S.A. Retrieval of an ice water path over the ocean from ISMAR and MARSS millimeter and submillimeter brightness temperatures. Atmos. Meas. Tech. 2018, 11, 611–632. [Google Scholar] [CrossRef] [Green Version]
  19. Håkansson, N.; Adok, C.; Thoss, A.; Scheirer, R.; Hörnquist, S. Neural network cloud top pressure and height for MODIS. Atmos. Meas. Tech. 2018, 11, 3177–3196. [Google Scholar] [CrossRef] [Green Version]
  20. Hedelt, P.; Efremenko, D.S.; Loyola, D.G.; Spurr, R.; Clarisse, L. Sulfur dioxide layer height retrieval from Sentinel-5 Precursor/TROPOMI using FP_ILM. Atmos. Meas. Tech. 2019, 12, 5503–5517. [Google Scholar] [CrossRef] [Green Version]
  21. Noia, A.D.; Hasekamp, O.P.; van Harten, G.; Rietjens, J.H.H.; Smit, J.M.; Snik, F.; Henzing, J.S.; de Boer, J.; Keller, C.U.; Volten, H. Use of neural networks in ground-based aerosol retrievals from multi-angle spectropolarimetric observations. Atmos. Meas. Tech. 2015, 8, 281–299. [Google Scholar] [CrossRef] [Green Version]
  22. Noia, A.D.; Hasekamp, O.P.; Wu, L.; van Diedenhoven, B.; Cairns, B.; Yorks, J.E. Combined neural network/Phillips–Tikhonov approach to aerosol retrievals over land from the NASA Research Scanning Polarimeter. Atmos. Meas. Tech. 2017, 10, 4235–4252. [Google Scholar] [CrossRef] [Green Version]
  23. Aires, F. Neural network uncertainty assessment using Bayesian statistics with application to remote sensing: 1. Network weights. J. Geophys. Res. 2004, 109. [Google Scholar] [CrossRef] [Green Version]
  24. Aires, F. Neural network uncertainty assessment using Bayesian statistics with application to remote sensing: 2. Output errors. J. Geophys. Res. 2004, 109. [Google Scholar] [CrossRef]
  25. Pfreundschuh, S.; Eriksson, P.; Duncan, D.; Rydberg, B.; Håkansson, N.; Thoss, A. A neural network approach to estimating a posteriori distributions of Bayesian retrieval problems. Atmos. Meas. Tech. 2018, 11, 4627–4643. [Google Scholar] [CrossRef] [Green Version]
  26. Arnold, V. On the representation of functions of several variables as a superposition of functions of a smaller number of variables. Math. Teach. Appl. Hist. Matem. Prosv. Ser. 2 1958, 3, 41–61. [Google Scholar]
  27. Cybenko, G. Approximation by superpositions of a sigmoidal function. Math. Control. Signals Syst. 1989, 2, 303–314. [Google Scholar] [CrossRef]
  28. Rumelhart, D.E.; Hinton, G.E.; Williams, R.J. Learning representations by back-propagating errors. Nature 1986, 323, 533–536. [Google Scholar] [CrossRef]
  29. Nix, D.; Weigend, A. Estimating the mean and variance of the target probability distribution. In Proceedings of the 1994 IEEE International Conference on Neural Networks (ICNN–94), Orlando, FL, USA, 28 June–2 July 1994. [Google Scholar] [CrossRef]
  30. Wright, W.; Ramage, G.; Cornford, D.; Nabney, I. Neural Network Modelling with Input Uncertainty: Theory and Application. J. VLSI Signal Process. 2000, 26, 169–188. [Google Scholar] [CrossRef] [Green Version]
  31. LeCun, Y.; Denker, J.; Solla, S. Optimal Brain Damage. In Advances in Neural Information Processing Systems; Touretzky, D., Ed.; Morgan-Kaufmann: Burlington, MA, USA, 1990; Volume 2, pp. 598–605. [Google Scholar]
  32. MacKay, D.J.C. A Practical Bayesian Framework for Backpropagation Networks. Neural Comput. 1992, 4, 448–472. [Google Scholar] [CrossRef]
  33. Blundell, C.; Cornebise, J.; Kavukcuoglu, K.; Wierstra, D. Weight Uncertainty in Neural Networks. In Proceedings of the 32nd International Conference on International Conference on Machine Learning, ICML’15, Lille, France, 6–11 July 2015; pp. 1613–1622. [Google Scholar]
  34. Kingma, D.P.; Welling, M. Auto-Encoding Variational Bayes. arXiv 2014, arXiv:1312.6114. [Google Scholar]
  35. Hinton, G.E.; Srivastava, N.; Krizhevsky, A.; Sutskever, I.; Salakhutdinov, R.R. Improving neural networks by preventing co-adaptation of feature detectors. arXiv 2012, arXiv:1207.0580. [Google Scholar]
  36. Srivastava, N.; Hinton, G.; Krizhevsky, A.; Sutskever, I.; Salakhutdinov, R. Dropout: A Simple Way to Prevent Neural Networks from Overfitting. J. Mach. Learn. Res. 2014, 15, 1929–1958. [Google Scholar]
  37. Gal, Y.; Ghahramani, Z. Dropout as a Bayesian Approximation: Representing Model Uncertainty in Deep Learning. arXiv 2016, arXiv:1506.02142. [Google Scholar]
  38. Ioffe, S.; Szegedy, C. Batch Normalization: Accelerating Deep Network Training by Reducing Internal Covariate Shift. arXiv 2015, arXiv:1502.03167. [Google Scholar]
  39. Teye, M.; Azizpour, H.; Smith, K. Bayesian Uncertainty Estimation for Batch Normalized Deep Networks. arXiv 2018, arXiv:1802.06455. [Google Scholar]
  40. Efremenko, D.S.; Molina García, V.; Gimeno García, S.; Doicu, A. A review of the matrix-exponential formalism in radiative transfer. J. Quant. Spectrosc. Radiat. Transf. 2017, 196, 17–45. [Google Scholar] [CrossRef] [Green Version]
  41. Efremenko, D.; Kokhanovsky, A. Foundations of Atmospheric Remote Sensing; Springer: Cham, Switzerland, 2021. [Google Scholar] [CrossRef]
  42. Efremenko, D.; Doicu, A.; Loyola, D.; Trautmann, T. Acceleration techniques for the discrete ordinate method. J. Quant. Spectrosc. Radiat. Transf. 2013, 114, 73–81. [Google Scholar] [CrossRef]
  43. Goody, R.; West, R.; Chen, L.; Crisp, D. The correlated-k method for radiation calculations in nonhomogeneous atmospheres. J. Quant. Spectrosc. Radiat. Transf. 1989, 42, 539–550. [Google Scholar] [CrossRef]
  44. Efremenko, D.; Doicu, A.; Loyola, D.; Trautmann, T. Optical property dimensionality reduction techniques for accelerated radiative transfer performance: Application to remote sensing total ozone retrievals. J. Quant. Spectrosc. Radiat. Transf. 2014, 133, 128–135. [Google Scholar] [CrossRef]
  45. Molina García, V.; Sasi, S.; Efremenko, D.S.; Doicu, A.; Loyola, D. Radiative transfer models for retrieval of cloud parameters from EPIC/DSCOVR measurements. J. Quant. Spectrosc. Radiat. Transf. 2018, 213, 228–240. [Google Scholar] [CrossRef] [Green Version]
  46. del Águila, A.; Efremenko, D.S.; Trautmann, T. A Review of Dimensionality Reduction Techniques for Processing Hyper-Spectral Optical Signal. Light Eng. 2019, 27, 85–98. [Google Scholar] [CrossRef] [Green Version]
  47. Schreier, F.; Gimeno García, S.; Hedelt, P.; Hess, M.; Mendrok, J.; Vasquez, M.; Xu, J. GARLIC—A general purpose atmospheric radiative transfer line-by-line infrared-microwave code: Implementation and evaluation. J. Quant. Spectrosc. Radiat. Transf. 2014, 137, 29–50. [Google Scholar] [CrossRef] [Green Version]
  48. Schreier, F. Optimized implementations of rational approximations for the Voigt and complex error function. J. Quant. Spectrosc. Radiat. Transf. 2011, 112, 1010–1025. [Google Scholar] [CrossRef]
  49. Gordon, I.; Rothman, L.; Hill, C.; Kochanov, R.; Tan, Y.; Bernath, P.; Birk, M.; Boudon, V.; Campargue, A.; Chance, K.; et al. The HITRAN2016 molecular spectroscopic database. J. Quant. Spectrosc. Radiat. Transf. 2017, 203, 3–69. [Google Scholar] [CrossRef]
  50. Bodhaine, B.A.; Wood, N.B.; Dutton, E.G.; Slusser, J.R. On Rayleigh Optical Depth Calculations. J. Atmos. Ocean. Technol. 1999, 16, 1854–1861. [Google Scholar] [CrossRef]
  51. Anderson, G.; Clough, S.; Kneizys, F.; Chetwynd, J.; Shettle, E. AFGL Atmospheric Constituent Profiles (0-120 km). AFGL-TR-86-0110; Air Force Geophysics Laboratory: Hanscom Air Force Base, MA, USA, 1986. [Google Scholar]
  52. Loyola R, D.G.; Pedergnana, M.; García, S.G. Smart sampling and incremental function learning for very large high dimensional data. Neural Netw. 2016, 78, 75–87. [Google Scholar] [CrossRef] [Green Version]
  53. Kingma, D.P.; Ba, J. Adam: A Method for Stochastic Optimization. arXiv 2017, arXiv:1412.6980. [Google Scholar]
  54. Rodgers, C. Inverse Methods for Atmospheric Sounding: Theory and Practice; Wolrd Scientific Publishing: Singapore, 2000. [Google Scholar]
  55. Doicu, A.; Trautmann, T.; Schreier, F. Numerical Regularization for Atmospheric Inverse Problems; Springer: Berlin, Germany, 2010. [Google Scholar]
  56. Tresp, V.; Ahmad, S.; Neuneier, R. Training Neural Networks with Deficient Data. In Advances in Neural Information Processing Systems; Cowan, J., Tesauro, G., Alspector, J., Eds.; Morgan-Kaufmann: Burlington, MA, USA, 1994; Volume 6, pp. 128–135. [Google Scholar]
  57. Kendall, A.; Gal, Y. What Uncertainties Do We Need in Bayesian Deep Learning for Computer Vision? arXiv 2017, arXiv:1703.04977. [Google Scholar]
  58. Gast, J.; Roth, S. Lightweight Probabilistic Deep Networks. arXiv 2018, arXiv:1805.11327. [Google Scholar]
  59. Jaulin, L.; Kieffer, M.; Didrit, O.; Walter, É. Applied Interval Analysis; Springer: London, UK, 2001. [Google Scholar] [CrossRef]
  60. Kearfott, R.B.; Dawande, M.; Du, K.; Hu, C. Algorithm 737: INTLIB—a portable Fortran 77 interval standard-function library. ACM Trans. Math. Softw. 1994, 20, 447–459. [Google Scholar] [CrossRef]
Figure 1. Retrieval results for the direct problem. The plots in the upper panels show the histograms of the relative error over the prediction set, while the lower plots show the predicted values (red) and the uncertainty intervals (gray) versus the true values.
Figure 1. Retrieval results for the direct problem. The plots in the upper panels show the histograms of the relative error over the prediction set, while the lower plots show the predicted values (red) and the uncertainty intervals (gray) versus the true values.
Remotesensing 13 05061 g001
Figure 2. Retrieval results obtained with Method 1 using a deterministic network.
Figure 2. Retrieval results obtained with Method 1 using a deterministic network.
Remotesensing 13 05061 g002
Figure 3. Retrieval results obtained with Method 1 using a Bayes-by-backprop network.
Figure 3. Retrieval results obtained with Method 1 using a Bayes-by-backprop network.
Remotesensing 13 05061 g003
Figure 4. Retrieval results obtained with Method 1 using a dropout network.
Figure 4. Retrieval results obtained with Method 1 using a dropout network.
Remotesensing 13 05061 g004
Figure 5. Retrieval results obtained with Method 1 using a batch normalized network.
Figure 5. Retrieval results obtained with Method 1 using a batch normalized network.
Remotesensing 13 05061 g005
Figure 6. Retrieval results obtained with Method 2 using a dropout network.
Figure 6. Retrieval results obtained with Method 2 using a dropout network.
Remotesensing 13 05061 g006
Figure 7. Retrieval results obtained with Method 2 using a Bayes-by-backprop network.
Figure 7. Retrieval results obtained with Method 2 using a Bayes-by-backprop network.
Remotesensing 13 05061 g007
Figure 8. Retrieval results obtained with Method 3 using assumed density filtering.
Figure 8. Retrieval results obtained with Method 3 using assumed density filtering.
Remotesensing 13 05061 g008
Figure 9. Retrieval results obtained with Method 3 using interval arithmetic. For a practical implementation of the algorithm, we use the interval arithmetic library INTLIB [60].
Figure 9. Retrieval results obtained with Method 3 using interval arithmetic. For a practical implementation of the algorithm, we use the interval arithmetic library INTLIB [60].
Remotesensing 13 05061 g009
Table 1. Average relative error E ( ε x ) ± E ( [ ε x E ( ε x ) ] 2 ) and standard deviation E ( σ x ) over the prediction set for the methods used to solve direct and inverse problems. The parameter x stands for the cloud optical thickness τ and cloud top height H. In the case of the inverse problem, the results correspond to Method 1 with a deterministic network (1a), Bayes-by-backprop (1b), dropout (1c), and batch normalization (1d); Method 2 with dropout (2a) and Bayes-by-backprop (2b); and Method 3 with assumed density filtering (3a) and interval arithmetic (3b).
Table 1. Average relative error E ( ε x ) ± E ( [ ε x E ( ε x ) ] 2 ) and standard deviation E ( σ x ) over the prediction set for the methods used to solve direct and inverse problems. The parameter x stands for the cloud optical thickness τ and cloud top height H. In the case of the inverse problem, the results correspond to Method 1 with a deterministic network (1a), Bayes-by-backprop (1b), dropout (1c), and batch normalization (1d); Method 2 with dropout (2a) and Bayes-by-backprop (2b); and Method 3 with assumed density filtering (3a) and interval arithmetic (3b).
MethodxErrorStd. Deviation
Direct Problem τ 2.94 × 10 2 ± 6.67 × 10 2 2.04 × 10 1
H 5.87 × 10 2 ± 1.72 × 10 1 3.66 × 10 1
1a τ 2.78 × 10 2 ± 1.30 × 10 1 8.92 × 10 1
H 1.96 × 10 2 ± 1.29 × 10 1 4.36 × 10 1
1b τ 3.39 × 10 2 ± 1.66 × 10 1 9.23 × 10 1
H 1.83 × 10 2 ± 1.33 × 10 1 8.11 × 10 1
1c τ 8.75 × 10 3 ± 2.25 × 10 2 3.23 × 10 1
H 3.45 × 10 3 ± 4.11 × 10 2 2.31 × 10 1
1d τ 1.01 × 10 2 ± 2.41 × 10 2 3.54 × 10 1
H 4.37 × 10 3 ± 2.73 × 10 2 2.29 × 10 1
2a τ 1.16 × 10 2 ± 4.05 × 10 2 8.24 × 10 1
H 1.63 × 10 2 ± 4.21 × 10 2 6.72 × 10 1
2b τ 3.10 × 10 2 ± 1.38 × 10 1 9.63 × 10 1
H 3.31 × 10 2 ± 4.57 × 10 2 4.88 × 10 1
3a τ 6.67 × 10 3 ± 2.48 × 10 2 6.55 × 10 1
H 8.18 × 10 4 ± 3.18 × 10 2 4.76 × 10 1
3b τ 7.08 × 10 3 ± 2.53 × 10 2 7.82 × 10 1
H 1.56 × 10 3 ± 3.33 × 10 2 6.53 × 10 1
Publisher’s Note: MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Share and Cite

MDPI and ACS Style

Doicu, A.; Doicu, A.; Efremenko, D.S.; Loyola, D.; Trautmann, T. An Overview of Neural Network Methods for Predicting Uncertainty in Atmospheric Remote Sensing. Remote Sens. 2021, 13, 5061. https://doi.org/10.3390/rs13245061

AMA Style

Doicu A, Doicu A, Efremenko DS, Loyola D, Trautmann T. An Overview of Neural Network Methods for Predicting Uncertainty in Atmospheric Remote Sensing. Remote Sensing. 2021; 13(24):5061. https://doi.org/10.3390/rs13245061

Chicago/Turabian Style

Doicu, Adrian, Alexandru Doicu, Dmitry S. Efremenko, Diego Loyola, and Thomas Trautmann. 2021. "An Overview of Neural Network Methods for Predicting Uncertainty in Atmospheric Remote Sensing" Remote Sensing 13, no. 24: 5061. https://doi.org/10.3390/rs13245061

Note that from the first issue of 2016, this journal uses article numbers instead of page numbers. See further details here.

Article Metrics

Back to TopTop