Model Selection in Atmospheric Remote Sensing with an Application to Aerosol Retrieval from DSCOVR/EPIC, Part 1: Theory

The limited information provided by satellite measurements does not allow for a complete determination of aerosol and cloud properties and, in particular, their microphysical properties. To deal with this problem, a set of models describing the microphysical properties is typically used. For example, the aerosol models implemented in the Moderate Resolution Imaging Spectroradiometer (MODIS) aerosol algorithm over land [1] consist of three spherical, fine-dominated types and one spheroidal, coarse-dominated type, while typical cloud models consist of water and ice clouds with a predefined particle shape and size. In this regard, a priori information such as the selection of a microphysical model is an essential part of the retrieval process. Standard retrieval algorithms usually ignore model uncertainty, i.e., a model is chosen from a given set of candidate models, and the retrieval is performed as if the selected model represents the true state. This approach is very risky and may result in large errors in the retrieved parameters if the model does not reflect the real scenario. An efficient way to quantify the uncertainty in model selection is the Bayesian approach and, in particular, Bayesian model selection and model averaging [2]. By model selection, we mean the specific problem of choosing the most suitable model from a given set of candidate models. In general, model selection is not a trivial task because for a given measurement, there will be several models that fit the data equally well. These tools were used, among others, in Refs. [3,4] to study uncertainty quantification in remote sensing of aerosols from Ozone Monitoring Instrument (OMI) data.

The key quantity in Bayesian model selection is relative evidence, which is a measure of how well the model fits the measurement. This is expressed in terms of the marginal likelihood, which is the integral over the state vector space of the a priori density times the likelihood density (see below). An accurate computation of the marginal likelihood is not a trivial task because some parameters characterizing the a priori and likelihood densities, e.g., the data error variance and the a priori state variance, are not precisely known and must be estimated. Moreover, the integral over the state vector space cannot be computed analytically because the likelihood density depends on the nonlinear forward model, for which an analytical representation is not available.

In this paper, we aim to eliminate the drawbacks of the Bayesian approach by using:

The paper is organized as follows. In Section 2, we derive a scaled data model with white noise in order to simplify the subsequent analysis. Section 3 is devoted to a review of the Bayesian parameter estimation and model selection in order to describe the problem and clarify the nomenclature. In Section 4, we summarize the iteratively regularized Gauss–Newton method and emphasize some special features of this method, including the estimation of the ratio of the data error variance and the a priori state variance. In Section 5, we assume a linearization of the forward model around the solution and, under this assumption, extend some regularization parameter choice methods for linear problems (the maximum marginal likelihood estimation and the generalized cross-validation) for data error variance estimation and model selection. Section 6 summarizes the computational steps of the proposed retrieval algorithm. In Section 7, we apply the algorithm to the retrieval of aerosol optical thickness and aerosol layer height from the Earth Polychromatic Imaging Camera (EPIC)/Deep Space Climate Observatory (DSCOVR) measurements.

Consider $N_{\mathrm{m}}$ microphysical models, and let ${\mathbf{F}}_m\left(\mathbf{x}\right)\in {\mathbb{R}}^M$ be the forward model corresponding to the model m for $m=1,\dots ,N_{\mathrm{m}}$. More specifically, in our analysis, ${\mathbf{F}}_m\left(\mathbf{x}\right)$ is the vector of the logarithms of the simulated radiances at all wavelengths ${\lambda }_i$, $i=1,\dots ,M$, i.e., ${\left[{\mathbf{F}}_m\left(\mathbf{x}\right)\right]}_i=ln{I}_m({\lambda }_i,\mathbf{x})$, and $\mathbf{x}\in {\mathbb{R}}^N$ is the state vector encapsulating the atmospheric parameters to be retrieved, e.g., aerosol optical thickness and layer height.

The nonlinear data model under examination consists of the equations: and: where ${\mathbf{y}}^{\delta }\in {\mathbb{R}}^M$ is the measurement vector or the noisy data vector, $\mathbf{y}$ the exact data vector, ${\delta }_{\mathrm{mes}}\in {\mathbb{R}}^M$ the measurement error vector, and ${\delta }_{\mathrm{aer}m}\in {\mathbb{R}}^M$ the model error vector, i.e., the error due to an inadequate model. Defining the (total) data error vector by: the data model (1) and (2) becomes:

In a deterministic setting, ${\delta }_m$ is characterized by the noise level ${\Delta }_m$ (defined as an upper bound for ${\delta }_m$, i.e., $\Vert {\delta }_m\Vert \le {\Delta }_m$), $\mathbf{x}$ is a deterministic vector, and we are faced with the solution of the nonlinear equation ${\mathbf{y}}^{\delta }={\mathbf{F}}_m\left(\mathbf{x}\right)$. In a stochastic setting, ${\delta }_m$ and $\mathbf{x}$ are random vectors, and the data model (3) is solved by means of a Bayesian approach.

Using the prewhitening technique, we transform the data model (3) into a model with white noise. For this purpose, we take:

If ${\delta }_{\mathrm{mes}}$ and ${\delta }_{\mathrm{aer}m}$ are independent random vectors, then ${\delta }_m$ is a Gaussian random vector with zero mean and covariance matrix: where: is the variance of the data error vector and: is the weighting factor giving the contribution of ${\mathcal{C}}_{\mathrm{mes}}$ to the covariance matrix ${\mathcal{C}}_{\delta m}$. After these preliminary constructions, we introduce the scaled data model: where: and: is the scaling matrix. Taking into account that: we see that by means of the prewhitening technique, the Gaussian random vector ${\delta }_m\sim \mathsf{N}(\mathbf{0},{\mathbf{C}}_{\delta m}={\sigma }_m^2{\mathcal{C}}_{\delta m})$ is transformed into the white noise ${\overline{\delta }}_m\sim \mathsf{N}(\mathbf{0},{\overline{\mathbf{C}}}_{\delta m}={\sigma }_m^2{\mathbf{I}}_M)$. Here and in the following, the notation $\mathsf{N}({\mathbf{x}}_{\mathrm{mean}},{\mathbf{C}}_{\mathrm{x}})$ stands for a normal distribution with mean ${\mathbf{x}}_{\mathrm{mean}}$ and covariance matrix ${\mathbf{C}}_{\mathrm{x}}$.

The following comments can be taken into consideration:

In this section, we review the basics of Bayesian parameter estimation and model selection by following the analysis given in Refs. [2,3].

In the framework of Tikhonov regularization, a regularized solution ${\mathbf{x}}_{m\alpha }^{\delta }$ to the nonlinear equation ${\overline{\mathbf{y}}}^{\delta }={\overline{\mathbf{F}}}_m\left(\mathbf{x}\right)$ is computed as the minimizer of the Tikhonov function: that is,

In classical regularization theory, the first term on the right-hand side of Equation (18) is the squared residual; the second one is the penalty term; while $\alpha$ and $\mathbf{L}$ are the regularization parameter and the regularization matrix, respectively. From Equations (12) and (18), we see that $V_{\alpha }(\mathbf{x}\mid {\overline{\mathbf{y}}}^{\delta },m)=(1/{\sigma }_m^2){\mathcal{F}}_{m\alpha }\left(\mathbf{x}\right)$, and from Equations (13) and (19), we deduce that ${\widehat{\mathbf{x}}}_{m\alpha }^{\delta }={\mathbf{x}}_{m\alpha }^{\delta }$. Thus, the maximum a posteriori estimate coincides with the Tikhonov solution, and therefore, the Bayesian parameter estimation can be regarded as a stochastic version of the method of Tikhonov regularization with an a priori chosen regularization parameter $\alpha$.

In the framework of Tikhonov regularization, the computation of the optimal regularization parameter is a crucial issue. With too little regularization, reconstructions deviate significantly from the a priori, and the solution is said to be under-regularized. With too much regularization, the reconstructions are too close to the a priori, and the solution is said to be over-regularized. In the Bayesian framework, the optimal regularization parameter is identified as the true ratio of the data error variance and the a priori state variance.

Several regularization parameter choice methods were discussed in [5]. These include methods with constant regularization parameters, e.g., the maximum likelihood estimation, the generalized cross-validation, and the nonlinear L-curve method. Unfortunately, at present, there is no fail-safe regularization parameter choice method that guarantees small solution errors in any circumstance, that is for any errors in the data vector.

An improvement of the problems associated with the regularization parameter selection is achieved in the framework of the so-called iterative regularization methods. These approaches are less sensitive to overestimates of the regularization parameter, but require more iteration steps to achieve convergence. A representative iterative approach is the iteratively regularized Gauss–Newton method.

At iteration step k of the iteratively regularized Gauss–Newton method, the forward model ${\overline{\mathbf{F}}}_m\left(\mathbf{x}\right)$ is linearized around the current iterate ${\mathbf{x}}_{mk}^{\delta }$, where: is the Jacobian matrix at ${\mathbf{x}}_{mk}^{\delta }$, and the nonlinear equation ${\overline{\mathbf{F}}}_m\left(\mathbf{x}\right)={\overline{\mathbf{y}}}^{\delta }$ is replaced by its linearization: where: is the step vector with respect to the a priori and: is the linearized data vector at iteration step k.

Since the nonlinear problem is ill-posed, its linearization is also ill-posed. Therefore, the linearized Equation (20) is solved by means of Tikhonov regularization with the penalty term ${\alpha }_k{\Vert \mathbf{L}\mathbf{p}\Vert }^2$, where the regularization parameters ${\alpha }_k$ are the terms of a decreasing (geometric) sequence, i.e., ${\alpha }_k=q{\alpha }_{k-1}$ with $q<1$. Note that in the method of Tikhonov regularization, the regularization parameter $\alpha$ is kept constant during the iterative process. The Tikhonov function for the linearized equation takes the form: and its minimizer is given by: where: is the regularized generalized inverse at ${\mathbf{x}}_{mk}^{\delta }$. The new iterate is computed as:

For the iterative regularization methods, the number of iteration steps k plays the role of the regularization parameter, and the iterative process is stopped after an appropriate number of steps $k^{☆}$ in order to avoid an uncontrolled expansion of errors in the data. In fact, a mere minimization of the residual $\Vert {\mathbf{r}}_{mk}^{\delta }\Vert$, where: is the residual vector at ${\mathbf{x}}_{mk}^{\delta }$, leads to a semi-convergent behavior of the iterated solution: while the error in the residual decreases as the number of iteration steps increases, the error in the solution starts to increase after an initial decay. A widely used a posteriori choice for the stopping index $k^{☆}$ is the discrepancy principle. According to this principle, the iterative process is terminated after $k^{☆}$ steps such that: with $\eta >1$; hence, the regularized solution is ${\mathbf{x}}_{m{k}^{☆}}^{\delta }$. Since the noise level cannot be estimated a priori for many practical problems arising in atmospheric remote sensing, we adopt a practical approach. This is based on the observation that the squared residual $\Vert {\mathbf{r}}_{mk}^{\delta }{\Vert }^2$ decreases during the iterative process and attains a plateau at approximately ${\Delta }_m^2$. Thus, if the nonlinear residuals $\Vert {\mathbf{r}}_{mk}^{\delta }\Vert$ converge to $\Vert {\mathbf{r}}_{m\infty }^{\delta }\Vert$ within a prescribed tolerance, we use the estimate:

The above heuristic stopping rule does not have any mathematical justification, but works sufficiently well in practice.

Since the amount of regularization is gradually decreased during the iterative processes, the iteratively regularized Gauss–Newton method can handle problems that practically do not require much regularization.

The numerical experiments performed in [5] showed that at the solution ${\mathbf{x}}_{m{k}^{☆}}^{\delta }$, (i) ${\alpha }_{k^{☆}-1}$ is close to the optimal regularization parameter, and (ii) ${\mathbf{x}}_{m{k}^{☆}}^{\delta }$ is close to the Tikhonov solution corresponding to the optimal regularization parameter. In this regard, we make the following assumptions:

Some consequences of the second assumption are the following. The optimality condition $\nabla {\mathcal{F}}_{m\widehat{\alpha }}\left({\mathbf{x}}_{m\widehat{\alpha }}^{\delta }\right)=0$, yields: and further, where ${\overline{\mathbf{K}}}_{m\widehat{\alpha }}={\overline{\mathbf{K}}}_{m{k}^{☆}}$, and:

Consequently, we see that: is the minimizer of the Tikhonov function: and that the residual vector of the linear equation ${\overline{\mathbf{y}}}_{m\widehat{\alpha }}^{\delta }={\overline{\mathbf{K}}}_{m\widehat{\alpha }}\mathbf{p}$ is equal to the nonlinear residual vector at the solution, i.e.,

In summary, in the iteratively regularized Gauss–Newton method:

Taking these results into account, we may conclude that the iteratively regularized Gauss–Newton method gives a solution to Problem 1 of the Bayesian parameter estimation.

In this section, we address Problems 2 and 3 related to Bayesian model selection.

To solve Problem 3, we suppose that the forward model can be linearized around the solution ${\mathbf{x}}_{m\widehat{\alpha }}^{\delta }$ [6]. In other words, in the first-order Taylor expansion: the remainder term ${\overline{\mathbf{R}}}_m(\mathbf{x},{\mathbf{x}}_{m\widehat{\alpha }}^{\delta })$ can be neglected. As a result, the nonlinear data model ${\overline{\mathbf{y}}}^{\delta }={\overline{\mathbf{F}}}_m\left(\mathbf{x}\right)+{\overline{\delta }}_m$ becomes the linear model: where ${\overline{\mathbf{y}}}_{m\widehat{\alpha }}^{\delta }$ is given by Equation (21), $\mathbf{p}=\mathbf{x}-{\mathbf{x}}_{\mathrm{a}}\sim \mathsf{N}(\mathbf{0},{\mathbf{C}}_{\mathrm{x}}={\sigma }_m^2{\left(\widehat{\alpha }{\mathbf{L}}^{\mathrm{T}}\mathbf{L}\right)}^{-1})$ and ${\overline{\delta }}_m\sim \mathsf{N}(\mathbf{0},{\overline{\mathbf{C}}}_{\delta m}={\sigma }_m^2{\mathbf{I}}_M)$. The direct consequences of dealing with the linear model (26) are the following results:

The replacement of the nonlinear data model by a linear data model will also enable us to deal with Problem 2. More specifically, for this purpose, we will adapt some regularization parameter choice methods for linear problems, i.e., maximum marginal likelihood estimation [7,8,9] and generalized cross-validation [10,11], to model selection. Furthermore, estimates for the data error variance ${\sigma }_m^2$ delivered by these methods will be used to compute the marginal likelihood density and, therefore, the relative evidence.

In Appendix B, it is shown that:

Denoting by $p_{\mathrm{L}}({\overline{\mathbf{y}}}^{\delta }\mid m)$ and $p_{\mathrm{G}}({\overline{\mathbf{y}}}^{\delta }\mid m)$ the marginal likelihoods corresponding to ${\widehat{\sigma }}_{\mathrm{L}m}^2$ and ${\widehat{\sigma }}_{\mathrm{G}m}^2$, respectively, we define:

Note that for the data error variance estimate (29), the marginal likelihood density becomes: with:

On the other hand, the statements (30) and (32) are equivalent to (cf. Equation (15)): and: respectively. Deviating from a stochastic interpretation of the relative evidence and regarding this quantity in a deterministic setting merely as a merit function characterizing the solution ${\mathbf{x}}_{m\widehat{\alpha }}^{\delta }$ , we define:

Note that $p_{\mathrm{MMLE}}(m\mid {\mathbf{y}}^{\delta })$ and $p_{\mathrm{GCV}}(m\mid {\mathbf{y}}^{\delta })$ do not depend on the data error variance ${\sigma }_m^2$. Note also that ${\widehat{\sigma }}_{\mathrm{G}m}^2$ is defined in terms of the square residual: while, according to Equation (29), ${\widehat{\sigma }}_{\mathrm{L}m}^2$ is defined in terms of the quantity ${\overline{\mathbf{y}}}_{m\widehat{\alpha }}^{\delta \mathrm{T}}({\mathbf{I}}_M-{\widehat{\mathbf{A}}}_{m\widehat{\alpha }}){\overline{\mathbf{y}}}_{m\widehat{\alpha }}^{\delta }$. The same difference exists between the numerators of the cross-validation function $v\left(m\right)$ and the marginal likelihood function $\lambda \left(m\right)$ (see Equations (32) and (30)).

An algorithmic implementation of the iteratively regularized Gauss–Newton method in connection with model selection is as follows:

The control parameters of the algorithm are: (i) the ratio q of the geometric sequence of regularization parameters, (ii) the minimum acceptable value of the regularization parameter ${\alpha }_{min}$, (iii) the tolerance ${\epsilon }_{\mathrm{R}}$ of the residual test convergence, and (iv) the tolerance $\eta$ of the discrepancy principle stopping rule.

We conclude our theoretical analysis with some comments:

The Deep Space Climate Observatory (DSCOVR) flies in a Lissajous orbit about the first Earth-Sun Lagrange point (L${}_1$), which is 1.5 million kilometers from the Earth towards the Sun. The Earth Polychromatic Imaging Camera (EPIC) is one of the two Earth observing instruments on board DSCOVR. This instrument views the entire sunlit side of the Earth from sunrise to sunset in the backscattering direction (scattering angles between 168.5° and 175.5°) with 10 narrowband filters, ranging from 317.5 to 779.5 nm [12]. This unique Earth observing geometry of EPIC compared to other instruments in Sun-synchronous orbits, that rarely view Earth at such large scattering angles [13], makes it a suitable candidate for several climate science applications including the measurement of cloud reflectivity and aerosol optical thickness and layer height [14].

We apply the algorithm proposed in Section 6 to the retrieval of aerosol optical thickness and aerosol layer height by generating synthetic measurements corresponding to the EPIC instrument. Specifically, Channels 7 and 8 in the Oxygen B-band at 680 and 687.75 nm, respectively, and Channels 9 and 10 in the Oxygen A-band at 764 and 779.5 nm, respectively, are used in the retrieval. The linearized radiative transfer model is based on the discrete ordinates method with matrix exponential and uses the correlated k-distribution method in conjunction with the principal component analysis technique to speed up the computations [15,16].

We consider the aerosol models implemented in the Moderate Resolution Imaging Spectroradiometer (MODIS) aerosol retrieval algorithm over land [1]. There are three spherical, fine-dominated model types (non-absorbing, moderately absorbing, and absorbing) and one spheroid, coarse-dominated model type (dust). These aerosol models, which depend on location and season, are the result of a cluster analysis of the climatology of almucantar retrievals from global AERONET (AErosol RObotic NETwork) measurements [17]. Each model consists of two log-normal modes (accumulated and coarse). A single log-normal mode is described by the number size distribution: where $r_{\mathrm{mod}}$ is the modal or median radius of the number size distribution, $\sigma$ is the standard deviation, and: is the total number of particles per cross-section of the atmospheric column. Table 1 displays the log-normal size parameters and the complex refractive indices $\mathrm{m}$ for the four aerosol models, where, for each log-normal mode, is the median radius of the volume size distribution: and: is the volume of particles per cross-section of the atmospheric column.

In our numerical analysis, the state vector is $\mathbf{x}={[\tau ,H]}^{\mathrm{T}}$, where $\tau$ is the aerosol optical thickness and H the aerosol layer height. The true aerosol optical thicknesses to be retrieved are ${\tau }_{\mathrm{t}}=0.25$, 0.50, 0.75, 1.0, 1.25, and 1.5, while for the true aerosol layer height, we take $H_{\mathrm{t}}=1.0$, 1.5, 2.0, 2.5, and 3.0 km. The a priori values, which coincide with the initial guesses, are ${\tau }_{\mathrm{a}}=2.0$ and $H_{\mathrm{a}}=4$ km. A Lambertian surface is assumed, and if not stated otherwise, the surface albedo is $A=0.06$. We generate synthetic measurements by choosing the moderately absorbing aerosol model (Model 2) as the true or exact model. Specifically, we:

The above scheme yields ${\mathcal{C}}_{\delta m}={\mathbf{I}}_M$ and $\mathbf{P}={\mathbf{I}}_M$, showing that the regularized solution does not depend on the scaling matrix. Coming to the regularization matrix, we first note that the choice ${\sigma }_{\mathrm{x}i}={\epsilon }_{\mathrm{x}}{\left[{\mathbf{x}}_{\mathrm{a}}\right]}_i$ for some ${\epsilon }_{\mathrm{x}}$ gives: