Open AccessArticle

A Spectral Acceleration Approach for the Spherical Harmonics Discrete Ordinate Method

Adrian Doicu

^*,

Dmitry S. Efremenko

and

Thomas Trautmann

German Aerospace Center (DLR), Remote Sensing Technology Institute, 82234 Oberpfaffenhofen, Germany

Author to whom correspondence should be addressed.

Remote Sens. 2020, 12(22), 3703; https://doi.org/10.3390/rs12223703

Submission received: 12 October 2020 / Revised: 9 November 2020 / Accepted: 10 November 2020 / Published: 11 November 2020

(This article belongs to the Special Issue Current Advances in Radiative Transfer Modeling for Satellite Optical Remote Sensing Applications)

Download

Browse Figures

Versions Notes

Abstract

A spectral acceleration approach for the spherical harmonics discrete ordinate method (SHDOM) is designed. This approach combines the correlated k-distribution method and some dimensionality reduction techniques applied on the optical parameters of an atmospheric system. The dimensionality reduction techniques used in this study are the linear embedding methods: principal component analysis, locality pursuit embedding, locality preserving projection, and locally embedded analysis. Through a numerical analysis, it is shown that relative to the correlated k-distribution method, PCA in conjunction with a second-order of scattering approximation yields an acceleration factor of 12. This implies that SHDOM equipped with this acceleration approach is efficient enough to perform spectral integration of radiance fields in inhomogeneous multi-dimensional media.

Keywords:

SHDOM; PCA; three dimensional radiative transfer

Graphical Abstract

1. Introduction

The retrieval of trace gas products from UV/VIS spectrometers is strongly affected by the presence of clouds. For the forward-model simulation of satellite measurements from instruments with high spatial resolution, it is important to account for the sub-pixel cloud inhomogeneity, or at the least, to assess this effect on the radiances at the top of the atmosphere, and so, on the retrieval results. For the new generation of atmospheric composition UV–VIS-NIR sensors, such as Sentinel-5P, and Sentinel-4 and −5, fast and accurate models accounting for the cloud inhomogeneities are crucial.

The spherical harmonics discrete ordinate method (SHDOM) developed by Evans [1,2] is one of the most efficient and widely used multi-dimensional deterministic (explicit) methods in the atmospheric sciences. SHDOM adopts (i) both spherical harmonics and discrete ordinates representation of the radiance field during different parts of the solution algorithm, and (ii) a sort of successive order of scattering solution method (Picard iteration [3]). The spherical harmonics are employed for computing the source function including the scattering integral, while the discrete ordinates are used to integrate the radiative transfer equation spatially. Specifically, the radiative transfer equation is solved iteratively by (i) transforming the source function from spherical harmonics to discrete ordinates, (ii) integrating the product of the source function and transmission along discrete ordinates, (iii) transforming the discrete ordinate radiances to spherical harmonics, and (iv) calculating the source function from the radiance in the spherical harmonics space. To achieve higher accuracy with a limited amount of memory, an adaptive grid is implemented. By this technique, regions where the source function is changing more rapidly have a higher density of grid points.

The computation of the radiative transfer through an atmosphere consisting of gas molecules is a demanding task because it requires an accurate modeling of the spectral absorption and scattering signatures of the gases plus additional particulate material under consideration. SHDOM uses a correlated k-distribution approach [4] for the integration over a spectral band. Note that in the version of March 1996, the parameterization of Fu and Liou [5] is the basis of the broadband integration of the gaseous line absorption, while longwave and shortwave Rapid Radiative Transfer Model (RRTM)-based k-distribution programs were added in 2003.

In addition to the k-distribution method, dimensionality reduction techniques applied on the optical properties of an atmospheric system become a standard broadband acceleration tool for one-dimensional radiative transfer [6,7,8]. Natraj et al. [9,10] designed a radiative transfer model based on the principal component analysis, while Efremenko et al. [11] generalized this model to include local linear embedding methods, as for example, locality pursuit embedding, locality preserving projection, and locally embedded analysis.

The aim of this paper was two-fold. First, to extend the applicability range of dimensionality reduction techniques to multi-dimensional radiative transfer, and second, to design a spectral acceleration approach for SHDOM that combines the correlated k-distribution method with some dimensionality reduction techniques. Although in the literature these methods are used separately, we apply them together. In this way, we expect to increase the computation speed of the k-distribution-based SHDOM.

2. The Spectral Acceleration Approach

Let us consider the solar radiative transfer in a rectangular prism of lengths

L_{x}

L_{y}

and

L_{z}

as shown in Figure 1. The top and bottom faces of the prism are denoted by

S_{t}

and

S_{b}

, respectively, while

S_{1 x}

(

x = 0

S_{2 x}

(

x = L_{x}

S_{1 y}

(

y = 0

), and

S_{2 y}

(

y = L_{y}

) correspond to the lateral faces. The boundary-value problem for the total radiance at point

r

in direction

Ω

consists of the inhomogeneous differential equation

\begin{matrix} \frac{d I}{d s} (λ, r, Ω) = - σ_{ext} (λ, r) I (λ, r, Ω) + \frac{σ_{sct} (λ, r)}{4 π} \int_{Ω} P (λ, r, Ω, Ω^{'}) I (λ, r, Ω^{'}) d Ω^{'}, \end{matrix}

(1)

the top-of-atmosphere boundary condition

I (λ, r_{t}, Ω^{-}) = \frac{F_{0} (λ)}{μ_{0}} δ (Ω^{-} - Ω_{0}), r_{t} \in S_{t},

(2)

and the Lambertian surface boundary condition

\begin{matrix} I (λ, r_{b}, Ω^{+}) & = \frac{A (λ)}{π} \int_{Ω^{-}} I (λ, r_{b}, Ω^{-}) | μ^{-} | d Ω^{-}, r_{b} \in S_{b} . \end{matrix}

(3)

At the horizontal boundaries, periodic boundary conditions are assumed, i.e.,

I (λ, r_{1 x}, Ω) = I (λ, r_{2 x}, Ω), I (λ, r_{1 y}, Ω) = I (λ, r_{2 y}, Ω),

(4)

for

r_{i x} \in S_{i x}

and

r_{i y} \in S_{i y}

with

i = 1, 2

. In Equations (1)–(4),

λ

is the wavelength,

σ_{ext} (λ, r)

and

σ_{sct} (λ, r)

the extinction and scattering coefficients, respectively,

P (λ, \cdot)

the phase function,

A (λ)

the surface albedo,

Ω_{0} = (- μ_{0}, φ_{0})

with

μ_{0} > 0

, the solar direction,

F_{0} (λ)

the solar flux,

r_{2 x} = r_{1 x} + L_{x} i

and

r_{2 y} = r_{1 y} + L_{y} j

with

(i, j, k)

being the Cartesian unit vectors. Moreover,

Ω^{+}

and

Ω^{-}

stand for an upward and a downward direction, respectively,

Ω

is the unit sphere, while

Ω^{+}

and

Ω^{-}

denote the upper and lower unit hemispheres, respectively,

An instrument placed at point

r_{m}

measures the radiance at the top of the atmosphere in direction

Ω_{m}^{0} = (μ_{m}^{0}, φ_{m}^{0})

μ_{m}^{0} > 0

. Denoting by

g (λ)

the slit function of the instrument, s the slit width, and

S_{tm}

the footprint of the instrument on the top face

S_{t}

, the spectral signal measured by the instrument at wavelength

λ_{m}

in the spectral interval

[λ_{\min}, λ_{\max}]

is given by

\begin{matrix} I (λ_{m}) & = \int_{λ_{\min} - s / 2}^{λ_{\max} + s / 2} g (λ_{m} - λ^{'}) I_{FOV} (λ^{'}) d λ^{'}, \end{matrix}

(5)

where

I_{FOV} (λ^{'}) = \int_{S_{t}} h (r_{t}) I (λ^{'}, r_{t}, Ω_{m} (r_{t})) d S_{t},

(6)

is the signal integrated over the field of view of the instrument at wavelength

λ^{'}

Ω_{m} (r_{t})

the unit vector along the line connecting the points

r_{t} \in S_{tm}

and

r_{m}

, i.e.,

Ω_{m} (r_{t}) = (r_{m} - r_{t}) / | r_{m} - r_{t} |

(see Figure 1), and

h (r_{t}) = \{\begin{matrix} \frac{1}{A_{tm}}, r_{t} \in S_{tm} \\ 0, otherwise \end{matrix},

(7)

and

A_{tm}

the characteristic function and the area of the instrument footprint, respectively. Assuming that the distance from the top of the atmosphere to the instrument is large, we approximate

I (λ^{'}, r, Ω_{m} (r)) \approx I (λ^{'}, r, Ω_{m}^{0})

in Equation (6).

For an atmosphere consisting of gas molecules and a cloud, we suppose that

the optical coefficients of the gas molecules depend on the altitude level and the wavelength, and
the optical coefficients of the cloud depend on the spatial coordinates but not on the wavelength.

The extinction coefficient is then computed as

σ_{ext} (λ, r) = σ_{ext}^{cloud} (r) + σ_{sct}^{mol} (λ, z) + σ_{abs}^{gas} (λ, z),

(8)

where

σ_{ext}^{cloud} (r)

is the extinction coefficient in the cloud,

σ_{sct}^{mol} (λ, z)

the molecular scattering coefficient due to Rayleigh scattering,

σ_{abs}^{gas} (λ, z)

the gas absorption coefficient, and

(x, y, z)

the Cartesian coordinates of point

r

. For the scattering coefficient, we have a similar representation, namely,

σ_{sct} (λ, r) = ω_{cloud} σ_{ext}^{cloud} (r) + σ_{sct}^{mol} (λ, z),

(9)

where

ω_{cloud}

is the single-scattering albedo of the cloud.

The goal of our analysis is to design an approach for accelerating the computation of the spectral signal

I (λ_{m})

at a set of

W_{m}

spectral points

{λ_{m w}}_{w = 1}^{W_{m}}

characterizing the instrument. For this purpose, we combine the correlated k-distribution method with some dimensionality reduction techniques.

2.1. Correlated k-Distribution Method

The correlated k-distribution method [4,12] is based on grouping spectral intervals according to absorption coefficient strength. In this section we provide a description of the correlated k-distribution method, which will enable us to connect this method with dimensionality reduction techniques.

Let

{{\bar{λ}}_{k}}_{k = 1}^{\bar{W}}

be a discrete set of

\bar{W}

equally spaced wavelengths in the interval

[λ_{\min} - s / 2, λ_{\max} + s / 2]

with

{\bar{λ}}_{1} - \frac{Δ λ}{2} = λ_{\min} - \frac{s}{2} and λ_{\bar{W}} + \frac{Δ λ}{2} = λ_{\max} + \frac{s}{2},

where

Δ λ

the discretization step. The forward-model spectral set

{{\bar{λ}}_{k}}_{k = 1}^{\bar{W}}

is finer than the measurement spectral set

{λ_{m w}}_{w = 1}^{W_{m}}

; in practice,

\bar{W}

is chosen as a multiple of

W_{m}

, i.e.,

\bar{W} = c W_{m}

with say, c being an integer greater than

20 .

The convolution integral giving the expression of the spectral signal measured by the instrument can be approximated by

\begin{matrix} I (λ_{m}) = \int_{λ_{\min} - s / 2}^{λ_{\max} + s / 2} g (λ_{m} - λ^{'}) I_{FOV} (λ^{'}) d λ^{'} \approx \sum_{k = 1}^{\bar{W}} g (λ_{m} - {\bar{λ}}_{k}) \int_{{\bar{λ}}_{k} - Δ λ / 2}^{{\bar{λ}}_{k} + Δ λ / 2} I_{FOV} (λ^{'}) d λ^{'} . \end{matrix}

(10)

As gas absorption has stronger spectral variation than molecular and particulate scattering, we write formally (at a given altitude level z)

I_{FOV} (λ) = I_{FOV} (σ_{abs}^{gas} (λ))

. The accurate technique for computing the integral in Equation (10) is a line-by-line (LBL) calculation of the gas absorption coefficient versus wavelength. Contrarily, the correlated k-distribution method is based on the fact that for a homogeneous atmosphere, the transmission within a spectral interval is independent of the LBL variation of

σ_{abs}^{gas}

with respect to the wavelength, but rather depends only on the distribution of

σ_{abs}^{gas}

within the spectral interval [13]. In this regard, let

F = F (σ_{abs, k}^{gas})

be the cumulative density function of

σ_{abs}^{gas} (λ)

in the spectral interval

[{\bar{λ}}_{k} - Δ λ / 2, {\bar{λ}}_{k} + Δ λ / 2]

, and

σ_{abs, k}^{gas} (F)

the quantile function (i.e., the inverse function of the cumulative distribution function). The spectral signal measured by the instrument can then be computed as

\begin{matrix} I (λ_{m}) = Δ λ \sum_{k = 1}^{\bar{W}} g (λ_{m} - {\bar{λ}}_{k}) \int_{0}^{1} I_{FOV} (σ_{abs, k}^{gas} (F)) d F = Δ λ \sum_{k = 1}^{\bar{W}} g (λ_{m} - {\bar{λ}}_{k}) \sum_{l = 1}^{N_{q}} ϖ_{l} I_{FOV} (σ_{abs, k}^{gas} (F_{l})), \end{matrix}

(11)

where

{F_{l}, ϖ_{l}}_{l = 1}^{N_{q}}

is a set of

N_{q}

quadrature points and weights in the interval

[0, 1]

. The

σ_{abs, k}^{gas} (F_{l})

can be computed by inverting the cumulative density functions of the LBL gas absorption coefficients, or in the case of the “exponential sum fitting of transmittance” method [14], by solving a nonlinear least squares problem.

In the final step, we define the sets of wavelengths and weights

{λ_{w}}_{w = 1}^{W}

and

{ω_{w}}_{w = 1}^{W}

with

W = \bar{W} N_{q}

, respectively, through the relations

λ_{w} = {\bar{λ}}_{k}

and

ω_{w} = Δ λ ϖ_{l}

, where

w = l + (k - 1) N_{q}

for

k = 1, \dots, \bar{W}, l = 1, \dots, N_{q}

. Accordingly, we set

σ_{abs}^{gas} (λ_{w}) = σ_{abs, k}^{gas} (F_{l})

. Note that

{λ_{w}}_{w = 1}^{W}

contains

\bar{W}

groups of

N_{q}

identical wavelengths. By this construction, Equation (11) becomes

I (λ_{m}) = \sum_{w = 1}^{W} ω_{w} g (λ_{m} - λ_{w}) I_{FOV} (σ_{abs}^{gas} (λ_{w})) .

(12)

Note that Equation (12) is a quadrature rule for the convolution integral (10) in the case of the correlated k-distribution method.

The computation of the spectral signal by means of the correlated k-distribution method requires the computation of the absorption coefficients

σ_{abs}^{gas} (λ_{w}, z)

at the set of wavelengths

{λ_{w}}_{w = 1}^{W}

, and so, W monochromatic radiative transfer calculations. For applications with

c \geq 20

and

N_{q} \geq 4

W = \bar{W} N_{q} = c W_{m} N_{q}

can be in the order of hundreds, and a further acceleration is needed. This can be achieved by applying dimensionality reduction techniques on the optical parameters.

2.2. Dimensionality Reduction of Atmospheric Optical Parameters

Following Ref. [11], the dimensionality reduction problem of the optical parameters can be formulated as follows. For each wavelength

λ_{w}

, we define an N-dimensional vector

x_{w}

x_{w}^{T} = [ln σ_{abs, 1}^{gas} (λ_{w}), \dots, ln σ_{abs, L}^{gas} (λ_{w}), ln σ_{sct, 1}^{mol} (λ_{w}), \dots, ln σ_{sct, L}^{mol} (λ_{w}), ln F_{0} (λ_{w})],

(13)

where

σ_{abs, k}^{gas}

and

σ_{sct, k}^{mol}

are the optical coefficients in the kth level,

N = 2 N_{z} + 1

, and

N_{z}

is the number of altitude levels. For simplicity, the surface albedo A was not included in

x_{w}

(we assumed that the variations of A are negligible over the width of a molecular absorption band). Thus, the wavelength variability of the optical parameters is encapsulated in the vector

x_{w} \in R^{N}

, and there is a one-to-one correspondence between the discrete set of wavelengths

{λ_{w}}_{w = 1}^{W}

, and the discrete set of optical parameters

{x_{w}}_{w = 1}^{W}

. For the N-dimensional data set

{x_{w}}_{w = 1}^{W}

\bar{x} = (1 / W) \sum_{w = 1}^{W} x_{w}

is the sample mean of the data. The goal of a dimensionality reduction technique is to find an M-dimensional subspace (

M < N

) spanned by a set of linear independent vectors

{a_{k}}_{k = 1}^{M}

, such that the centered data

x_{w} - \bar{x}

belong to this subspace; that is,

x_{w} \approx \bar{x} + \sum_{k = 1}^{M} y_{w k} a_{k} = \bar{x} + A y_{w}, w = 1, \dots, W .

(14)

In Equation (14), the matrix

A = {[a_{k}]}_{k = 1}^{M} \in R^{N \times M}

, encapsulating the column vectors

a_{k}

, represents the inverse transform from the low-dimensional space to the high-dimensional space, and

y_{w k}

is the kth component of

y_{w} \in R^{M}

. The vector

y_{w}

is given by the forward transform from the high-dimensional space to the low-dimensional space,

y_{w} = A^{†} (x_{w} - \bar{x}),

(15)

where

A^{†} = {(A^{T} A)}^{- 1} A^{T} \in R^{M \times N}

is the pseudoinverse of

A

An approximate model for estimating the integrated signal of the instrument at wavelength

λ

is of the form

ln \frac{I_{FOV} (λ)}{I_{FOV}^{a} (λ)} = f (λ),

(16)

where

I_{FOV}^{a} (λ)

is the integrated signal computed by an approximate radiative transfer model, and

f (λ)

a correction factor. The correction factor

f (λ)

can be calculated efficiently and accurately by means of dimensionality reduction techniques applied on the optical parameters. Let us assume that the scalar function

f (x_{w})

given by Equation (16) is not too nonlinear in

x_{w}

. The computational process is organized as follows.

For (cf. Equation (14))

$Δ x_{w} = x_{w} - \bar{x} = \sum_{k = 1}^{M} y_{w k} a_{k} = A y_{w},$

(17)

$f (x_{w})$ is approximated by a second-order Taylor expansion around $\bar{x}$ ; that is,

$f (x_{w}) \approx f (\bar{x} + Δ x_{w}) \approx f (\bar{x}) + Δ x_{w}^{T} \nabla f (\bar{x}) + \frac{1}{2} Δ x_{w}^{T} \nabla^{2} f (\bar{x}) Δ x_{w},$

(18)

where $\nabla f$ and $\nabla^{2} f$ are the gradient and the Hessian of f, respectively.
For $a_{k} = | | a_{k} | | {\hat{a}}_{k}$ and in view of Equation (17), the second and third terms in Equation (18) are written as

$Δ x_{w}^{T} \nabla f (\bar{x}) = \sum_{k = 1}^{M} y_{w k} | | a_{k} | | \frac{\partial f}{\partial {\hat{a}}_{k}} (\bar{x})$

(19)

and

$Δ x_{w}^{T} \nabla^{2} f (\bar{x}) Δ x_{w} = \sum_{k = 1}^{M} \sum_{l = 1}^{M} y_{w k} y_{w l} | | a_{k} | | | | a_{l} | | \frac{\partial^{2} f}{\partial {\hat{a}}_{k} \partial {\hat{a}}_{l}} (\bar{x}),$

(20)

respectively.
The mixed directional derivatives in Equation (20) are neglected, while the first- and second-order directional derivatives are approximated by means of central differences; that is,

$\frac{\partial f}{\partial {\hat{a}}_{k}} (\bar{x}) \approx \frac{f (\bar{x} + h {\hat{a}}_{k}) - f (\bar{x} - h {\hat{a}}_{k})}{2 h},$

(21)

and

$\frac{\partial^{2} f}{\partial {\hat{a}}_{k}^{2}} (\bar{x}) \approx \frac{f (\bar{x} + h {\hat{a}}_{k}) - 2 f (\bar{x}) + f (\bar{x} - h {\hat{a}}_{k})}{h^{2}},$

(22)

respectively.
For $h = | | a_{k} | |$ , Equations (18)–(22) yield the computational formula:

$\begin{matrix} f (x_{w}) \approx f (\bar{x}) + \frac{1}{2} \sum_{k = 1}^{M} [f (\bar{x} + a_{k}) - f (\bar{x} - a_{k})] y_{w k} + \frac{1}{2} \sum_{k = 1}^{M} [f (\bar{x} + a_{k}) - 2 f (\bar{x}) + f (\bar{x} - a_{k})] y_{w k}^{2} . \end{matrix}$

(23)

From Equation (23), it is apparent that the computation of

f (x_{w})

requires

2 M + 1

calls of the exact and approximate models, while the computation of

I_{FOV} (λ_{w})

requires, in addition, W calls of the approximate model. In this regard and taking into account that

M ≪ W

, we are led to a substantial reduction of the computational time.

The most widely used dimensionality reduction techniques are the linear embedding methods which are summarized below.

The principal component analysis (PCA) [15] performs a dimensionality reduction by projecting the original N-dimensional data on the M-dimensional subspace spanned by the dominant singular vectors of the data covariance matrix.
The locality pursuit embedding (LPE) [16] performs a principal component analysis on local nearest neighbor patches to reveal the tangent space structure on the M-dimensional subspace.
The locality preserving projection (LPP) [17] solves a variational problem that optimally preserves the neighborhood structure of the data set.
The locally embedded analysis (LEA) [18] uses an embedding strategy based on a linear eigenspace analysis to minimize the local reconstruction error.

It should be pointed out that PCA preserves only the global structure of the data, and may fail to preserve the local structure if the data lie on a nonlinear manifold. In contrast, LPE, LPP, and LEA are local linear approaches, which optimally preserve local neighborhood information (the local structure of the data) in a certain sense.

The PCA, LPE, LPP, and LEA methods are illustrated in Algorithms A1–A4 of Appendix A. The outputs of these algorithms are the linear mapping

A

and the vectors of parameters

y_{w}

w = 1, \dots, W

. Note that in the PCA language, the column vectors

a_{k}

k = 1, \dots, M

A

are the so-called empirical orthogonal functions, while the vectors of parameters

y_{w}

w = 1, \dots, W

are the principal components.

A last problem which has to be solved is the choice of an approximate model. Two special features of SHDOM facilitate this choice, namely, the use (i) of spherical harmonics and discrete ordinates to represent the radiance field, and (ii) of a sort of successive order of scattering solution method (Picard iteration). In a first approximate model, we choose 2 zenith discrete angles per hemisphere to represent the radiance field; analogous to the one-dimensional radiative transfer, the resulting method is called the four-stream approximation. In a second approximate model, we consider only 2 steps of the method of Picard iteration; the resulting method is called the second-order of scattering approximation.

In summary, the spectral acceleration approach includes two steps. The first step is a pre-processing step, consisting in the computation of the sets of wavelengths

{λ_{w}}_{w = 1}^{W}

, the absorption coefficients

σ_{abs}^{gas} (λ_{w}, z)

, and the weighting factors

{ω_{w}}_{w = 1}^{W}

in the framework of the correlated k-distribution method. The second step is a computational step, consisting of the calculation of

the data vectors $x_{w}$ , $w = 1, \dots, W$ , the empirical orthogonal functions $a_{k}$ , $k = 1, \dots, M$ , and the principal components $y_{w}$ , $w = 1, \dots, W$ in the framework of dimensionality reduction techniques,
the correction factor $f (x_{w})$ , $w = 1, \dots, W$ , and the integrated signal $I_{FOV} (λ_{w})$ , $w = 1, \dots, W$ according to Equations (23) and (16), respectively, and finally,
the spectral signal measured by the instrument $I (λ_{mw})$ , $w = 1, \dots, W_{m}$ by means of Equation (10).

3. Numerical Analysis

Before presenting some numerical results, we describe an implementation of SHDOM that is devoted to the retrieval of atmospheric trace gas concentrations from space-borne spectral measurements of radiation reflected through the Earth’s atmosphere.

3.1. SHDOM Implementation

For an atmosphere consisting of gas molecules, we have to consider a domain of analysis with a large vertical extent, especially when dealing with trace retrievals. The reason is that gas molecules have a vertical concentration profile up to several tens of kilometers. Unfortunately, for large problems, SHDOM is limited by the amount of memory available for calculation, in particular on smaller computers. To handle this problem, in our implementation of SHDOM, we divide the domain of analysis into a number of subdomains along the vertical direction, and for each subdomain, we store the radiance and source function in separate files. The optical properties and the direct solar beam are calculated at each point on the base grid of the entire domain of analysis, while (i) the radiative transfer equation integration along each discrete ordinate during the solution iterations and (ii) the integration along the viewing directions for the intensity output are done successively moving from one subdomain to the neighboring one. Each subdomain has its own grid obtained by the cell splitting (adaptive grid) procedure. At the boundaries of a subdomain, the radiances are calculated not only at the grid points of the subdomain itself, but also at the grid points corresponding to the boundaries of the neighboring subdomains (Figure 2). This is done by integrating the radiative transfer equation and by applying the standard interpolation approach (the entering and exiting values of the extinction and extinction-source function product are computed using bilinear interpolation of the four gridpoint values of the faces pierced by a ray). To ensure the continuity condition for the radiance field, the boundary values of the radiance are passed between the neighboring subdomains. The data are transferred to and from files.

This implementation of SHDOM does not give exactly the same results as the standard implementation because of (i) the splitting decision for the cells situated at the boundaries of a subdomain, and (ii) the differences between the rays along which the radiative transfer equation is integrated (a ray is traced backward until the transmission falls below some minimum specified value only inside a subdomain). However, as compared to the overall accuracy of the algorithm, these errors are small. On the other hand, this implementation increases the computation time by the times which are needed (i) to write and read data to and from files, and (ii) to compute the radiances at the additional grid points on the boundaries of each subdomain. It should be pointed out that the number of subdomains is set by the user and depends on the problem to be solved.

The code can be reformulated as a parallel code by using the Message Passing Interface (MPI). In this context, it should be pointed out that a parallel implementation of SHDOM, but in which the domain of analysis is divided along the horizontal directions, has been described in Ref. [19]. However, in the case of horizontal splitting, the transfer of data between different neighboring subdomains is more challenging for the implementation.

3.2. Numerical Results

The computational domain is a rectangular prism of lengths

L_{x} = L_{y} = 15

km and

L_{z} = 40

km. The discretization steps along the horizontal directions are

Δ x = Δ y = 0.5

km; thus, the corresponding numbers of base grid points are

N_{x} = N_{y} = 30

. Along the vertical direction, the domain is discretized into 8 subdomains; the corresponding altitude levels and discretization steps are listed in Table 1. The cloud is homogeneous in the vertical direction and is placed between 3 and 4 km (second subdomain). The cloud extinction field is modeled as

σ_{ext}^{cloud} (x, y) = σ_{0} f (x, y),

where

σ_{0} = 6 {km}^{- 1}

and

f (x, y)

is the indicator function taking the values

f (x, y) = 1

inside the cloud and

f (x, y) = 0

inside the clear sky region. The indicator function

f (x, y)

is generated by a two-dimensional broken cloud model [20] with a cloud fraction of about

0.4 .

In order to avoid abrupt changes of the extinction field in the horizontal plane,

σ_{ext}^{cloud} (x, y)

is smoothed at the boundary of a cloudy region (by linear interpolation along two discretization steps in each direction). To simplify the analysis, a Henyey–Greenstein phase function with an asymmetry parameter

g = 0.8

is considered, while the single-scattering albedo is chosen as

ω_{cloud} = 0.99

. The number of discrete zenith and azimuth angles are

N_{μ} = 16

and

N_{φ} = 2 N_{μ}

, respectively, the solar and instrument zenith angles are

θ_{0} = 60^{\circ}

and

θ_{m}^{0} = 30^{\circ}

, respectively, the relative azimuth angle is

Δ φ = 0

, and a Lambertian reflecting surface with the surface albedo

A = 0.2

is assumed. The footprint of the detector is considered to be a square of length

2 a = 10 Δ x

centered at

x_{0} = y_{0} = L_{x} / 2

, and

z_{0} = L_{z}

. The calculations are performed for a mid-latitude summer atmosphere [21] and a wavelength-dependent slit function corresponding to the TROPOspheric Monitoring Instrument (TROPOMI).

SHDOM is run by using an adaptive grid with a splitting accuracy of

10^{- 4}

. The solution accuracy is

10^{- 4}

, and the simulations are performed for periodic boundary conditions by using the delta-M approximation along with the untruncated phase function single-scattering solution (TMS correction of Nakajima and Tanaka [22]). The dimensionality reduction parameter M, which determines the approximation accuracy of the linear embedding methods, is chosen as

M = 4

. The simulations were performed on a computer Intel Core i5-3340M CPU 2.70GHz with 8 GB RAM.

3.2.1. Oxygen A-Band Test Problem

In our first example, we consider, in addition to the scattering and absorption by a cloud, molecular Rayleigh scattering, and the absorption in the Oxygen A-band. The measurement spectral grid contains 107 spectral points between 758 nm and 771 nm, while the number of spectral points for the correlated k-distribution method is

W = 936

. For this test, a two-dimensional slice at

y = 7

km is selected from a three-dimensional broken-cloud field (Figure 3). Figure 4 illustrates the spectral signal computed by using SHDOM with the correlated k-distribution method. Using these results as a reference, we show in Figure 5 the relative differences in the spectral signal corresponding to the linear embedding methods, and in Table 2, the RMS of the relative difference and the computational time in hours, minutes, and seconds. The computation time of SHDOM with the correlated k-distribution method was 36 minutes and 28 s. Several conclusions can be formulated:

For the dimensionality reduction parameter $M = 4$ , the relative differences are smaller than $1.25 \times 10^{- 3}$ over the entire spectral domain, while the RMS values are smaller than $6.4 \times 10^{- 4}$ .
The computation time corresponding to the second-order of scattering approximation is smaller than that corresponding to the four-stream approximation.
For the second-order of scattering approximation, the fastest linear embedding method is PCA and the most accurate is LPP.
Relative to the correlated k-distribution method, the acceleration factor is of about 7–12.

3.2.2. ${NO}_{2}$ -Test Problem

In the second example, the absorption of

{NO}_{2}

, ozone (

O_{3}

), oxygen dimer (

O_{4}

), and water vapor is considered. The measurement spectral grid contains 119 spectral points between 425 nm and 450 nm, while the number of spectral points for the correlated k-distribution method is

W = 268

. Note that in the correlated k-distribution method, the auxiliary gases are treated as weak absorbers meaning that for these gases, the average values of the LBL cross sections are used in the calculation. The test problem is three-dimensional, and the indicator function of the broken cloud is illustrated in Figure 6. In Table 3, we indicate the amount of memory required for storing the radiance and source function for each subdomain. As a result that the domain of analysis is not excessively large, the total amount of disk memory is 3.168 GB. As a result, we infer that SHDOM can be also run without splitting the domain of analysis. Taking these results as a reference, we found that the relative errors due to the domain splitting procedure are smaller than

10^{- 3}

. The spectral signal computed by using SHDOM with the correlated k-distribution method is shown in Figure 7, the relative differences in the spectral signal corresponding to the linear embedding methods are plotted in Figure 8, and finally, the RMS of the relative difference and the computational time in hours, minutes, and seconds are given in Table 4. For this test example, the computation time of SHDOM with the correlated k-distribution method was 12 h, 18 min, and 32 s. As in the first test example, we see that

for the dimensionality reduction parameter $M = 4$ , the relative differences are smaller than $5 \times 10^{- 4}$ over the entire spectral domain, while the RMS values are smaller than $2.2 \times 10^{- 4}$ ;
the second-order of scattering approximation is faster than the four-stream approximation;
for the second-order of scattering approximation, the fastest linear embedding method is PCA and the most accurate is LPE;
the acceleration factor is about 9–12.

4. Conclusions

A spectral acceleration approach for SHDOM has been designed. This approach combines the correlated k-distribution method and the linear embedding methods PCA, LPE, LPP, and LEA.

We have performed a numerical analysis of these techniques for a two-dimensional atmosphere in the Oxygen A-band, and a three-dimensional atmosphere consisting of nitrogen, ozone, oxygen dimer, and water vapor in a spectral interval ranging from 425 nm and 450 nm. The main conclusions of our analysis are as follows.

SHDOM with the correlated k-distribution and linear embedding methods has a sufficiently high accuracy. The relative differences in the spectral signal are smaller than $1.25 \times 10^{- 3}$ (over the entire spectral domain) in the case of a two-dimensional atmosphere, and $5.0 \times 10^{- 4}$ in the case of a three-dimensional atmosphere. In three of the four cases, the most accurate linear embedding method is LPE.
The linear embedding methods based on the second-order of scattering approximation are faster than those based on the four-stream approximation. Specifically, in the case of a two-dimensional atmosphere, the computation time of the linear embedding methods is on average about 4 min and 30 s for the four-stream approximation, and 3 min and 20 s for the second-order of scattering approximation, while in the case of a three-dimensional atmospheres, the corresponding times are 1 h and 21 min for the four-stream approximation, and 1 h and 2 min for the second-order of scattering approximation. The fastest linear embedding methods is PCA followed by LPE. For the test examples considered in our numerical analysis, PCA in conjunction with a second-order of scattering approximation yields an acceleration factor of 12 relative to the correlated k-distribution method.

In conclusion, from the point of view of accuracy and efficiency, it appears that LPE based on the second-order of scattering approximation is the most suitable acceleration method. Note that the accuracy and efficiency of the linear embedding methods are determined by the dimensionality reduction parameter M (as the accuracy and the computation time increase with M, an optimal value for M should be a compromise between accuracy and efficiency).

SHDOM equipped with the spectral acceleration approach can be used to analyze the impact of cloud inhomogeneities on trace gas retrievals [23]. In particular, for a specific cloudy scene, a spectral signal can be simulated by SHDOM and then included in a one-dimensional retrieval algorithm to derive for example, the total column amount of a trace gas. By comparing the retrieved total column with the true value, a bias due to cloud effects can be deduced [24,25], and possible correction strategies can be explored. Going one step further, a multi-dimensional trace gas retrieval algorithm can be designed. In this case, the spectral acceleration approach should be adapted to a linearized version of SHDOM (based on a forward or a forward-adjoint approach as described in Ref. [26]). This is the topic of a companion paper.

Author Contributions

Conceptualization, A.D.; software, A.D. and D.S.E.; formal analysis, D.S.E. and T.T.; writing—original draft preparation, A.D.; writing—review and editing, D.S.E. 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. The funders had no role in the design of the study; in the collection, analyses, or interpretation of data; in the writing of the manuscript, or in the decision to publish the results.

Appendix A

The PCA, LPE, LPP, and LEA methods are described in Ref. [11]. For convenience of the reader, the related algorithms are summarized below.

Algorithm A1: PCA.

Step 1.: Stack all centered data $x_{w} \leftarrow x_{w} - \bar{x}$ , $w = 1, \dots, W$ , into the columns of a matrix $X \in R^{N \times W}$ , i.e., $X = {[x_{w}]}_{w = 1}^{W}$ .
Step 2.: Calculate the covariance matrix $C_{x} = (1 / W) X X^{T} \in R^{N \times N}$ .
Step 3.: Calculate a singular value decomposition of the covariance matrix $C_{x} = U Σ U^{T}$ , where $Σ = diag {[σ_{k}]}_{k = 1}^{N} \in R^{N \times N}$ is the diagonal matrix of the singular values $σ_{1} > σ_{2} > \dots > σ_{N} > 0$ , and $U = {[u_{k}]}_{k = 1}^{N} \in R^{N \times N}$ is the orthogonal matrix of the singular vectors $u_{k}$ .
Step 4.: Take $A = U_{M} = {[u_{k}]}_{k = 1}^{M} \in R^{N \times M}$ yielding $A^{†} = U_{M}^{T} \in R^{M \times N}$ , or alternatively, consider the scaled orthogonal vectors ${\bar{u}}_{k} = \sqrt{σ_{k}} u_{k}$ , and take
$A = {\bar{U}}_{M} = U_{M} Σ_{M}^{1 / 2} \in R^{N \times M}$ yielding $A^{†} = {\bar{U}}_{M}^{†} = Σ_{M}^{- 1 / 2} U_{M}^{T} \in R^{M \times N}$ , with $Σ_{M} = diag {[σ_{k}]}_{k = 1}^{M}$ .
Step 5.: Compute the vectors of parameters $y_{w} = A^{†} x_{w}$ for $w = 1, \dots, W$ .

Algorithm A2: LPE.

Step 1.: Stack all centered data $x_{w} \leftarrow x_{w} - \bar{x}$ , $w = 1, \dots, W$ , into the columns of a matrix
$X \in R^{N \times W}$ , i.e., $X = {[x_{w}]}_{w = 1}^{W}$ .
Step 2.: For each $w = 1, \dots W$ , consider the sets $X (w) = {x_{N (j)}^{(w)}}_{j = 1}^{k}$ containing the k nearest neighbors to $x_{w}$ with respect to the Euclidean distance.
Step 3.: Set $L = 0 \in R^{W \times W}$ .
Step 4.: For $w = 1, \dots W$ , compute the diagonal matrix $D_{w} = diag {[I_{v \in X (w)}]}_{v = 1}^{W} \in R^{W \times W}$ , where the indicator function $I_{v \in X (w)}$ is defined by $I_{v \in X (w)} = 1$ if $x_{v} \in X (w)$ , and $I_{v \in X (w)} = 0$ otherwise.
Step 5.: For all $w = 1, \dots W$ , update the matrix $L$ recursively as $L \leftarrow L + D_{w} - (1 / k) D_{w} 1 1^{T} D_{w}^{T}$ , where $1 \in R^{W}$ is a column vector taking 1 at each entry.
Step 6.: Calculate the symmetric and positive semi-definite matrix $A = X L X^{T} \in R^{N \times N}$ , compute a singular value decomposition of $A$ , $A = U Σ U^{T},$ and as in Algorithm A1, take $A = U_{M} Σ_{M}^{1 / 2} \in R^{N \times M}$ and $A^{†} = Σ_{M}^{- 1 / 2} U_{M} \in R^{M \times N}$ .
Step 7.: Compute the vectors of parameters $y_{w} = A^{†} x_{w}$ for $w = 1, \dots, W$ .

Algorithm A3: LPP.

Step 2.: Stack all centered data $x_{w} \leftarrow x_{w} - \bar{x}$ , $w = 1, \dots, W$ , into the columns of a matrix
$X \in R^{N \times W}$ , i.e., $X = {[x_{w}]}_{w = 1}^{W}$ .
Step 2.: For $w = 1, \dots W$ , consider the sets $X (w) = {x_{N (j)}^{(w)}}_{j = 1}^{k}$ containing the k nearest neighbors to $x_{w}$ with respect to the Euclidean distance.
Step 3.: For $w, v = 1, \dots W$ , compute the elements $W_{w v}$ of the sparse matrix $W \in R^{W \times W}$ as $W_{w v} = \exp (- | | x_{w} - x_{v} {| |}^{2} / t)$ if $x_{v} \in X (w)$ , and $W_{w v} = 0$ otherwise. Alternatively, chose $W_{w v} = 1$ if $x_{v} \in X (w)$ , and $W_{w v} = 0$ otherwise.
Step 4.: Calculate the elements $D_{v}$ of the diagonal matrix $D = diag {[D_{v}]}_{v = 1}^{W} \in R^{W \times W}$ as
$D_{v} = \sum_{w} W_{w v}$ , and the Laplacian matrix $L = D - W \in R^{W \times W}$ .
Step 5.: For $A = X L X^{T} \in R^{N \times N}$ and $B = X D X^{T} \in R^{N \times N}$ , solve the generalized eigenvalue problem $A z_{k} = λ_{k} B z_{k}$ as follows: (i) compute a singular value decomposition of the matrix $B = U Σ U^{T}$ , (ii) calculate the matrix $D = Σ^{- 1 / 2} U^{T} \in R^{N \times N}$ and its inverse $D^{- 1} = U Σ^{1 / 2} \in R^{N \times N}$ , and finally, (iii) compute the matrix $C = D A D^{T} \in R^{N \times N}$ and its singular value decomposition $C = V Λ V^{T}$ . The eigenvalues $λ_{k}$ are the singular values of $C$ , while the eigenvectors $z_{k}$ are given by $z_{k} = D^{T} v_{k}$ , where $v_{k}$ are the singular vectors of $C$ .
Step 6.: Construct the matrix $V_{M} = {[v_{k + N - M}]}_{k = 1}^{M} \in R^{N \times M}$ , where $v_{k}$ is the kth column vector of $V \in R^{N \times N}$ , and take $A = D^{- 1} V_{M} \in R^{N \times M}$ and $A^{†} = V_{M}^{T} D \in R^{M \times N}$ .
Step 7.: Compute the vectors of parameters $y_{w} = A^{†} x_{w}$ for $w = 1, \dots, W$ .

Algorithm A4: LEA.

Step 1.: Stack all centered data $x_{w} \leftarrow x_{w} - \bar{x}$ , $w = 1, \dots, W$ , into the columns of a matrix $X \in R^{N \times W}$ , i.e., $X = {[x_{w}]}_{w = 1}^{W}$ .
Step 2.: For each $w = 1, \dots W$ , consider the sets $X (w) = {x_{N (j)}^{(w)}}_{j = 1}^{k}$ containing the k nearest neighbors to $x_{w}$ with respect to the Euclidean distance.
Step 3.: Calculate the elements of the sparse matrix $W \in R^{W \times W}$ as follows: (i) for $w = 1, \dots W$ , calculate the elements $G_{i j}$ of the matrix $G \in R^{k \times k}$ as $G_{i j} = {(x_{w} - x_{N (i)}^{(w)})}^{T} (x_{w} - x_{N (j)}^{(w)})$ , and (ii) set $W_{w N (j)} = w_{j}$ with $w_{j}$ being the elements of the column vector $w = G^{- 1} 1 / (1^{T} G^{- 1} 1) \in R^{k}$ .
Step 4.: Calculate the elements $D_{v}$ of the diagonal matrix $D = diag {[D_{v}]}_{v = 1}^{W} \in R^{W \times W}$ as $D_{v} = \sum_{w} W_{w v}$ .
Step 5.: As in Algorithm A3, solve the generalized eigenvalue problem $A z_{k} = λ_{k} B z_{k}$ , for $A = X {(D - W)}^{T} (D - W) X^{T} \in R^{N \times N}$ and $B = X (D^{T} D + W^{T} W) X^{T} \in R^{N \times N}$ , and let $C = V Λ V^{T}$ .
Step 6.: Construct the matrix $V_{M} = {[v_{k + N - M}]}_{k = 1}^{M} \in R^{N \times M}$ , where $v_{k}$ is the kth column vector of $V \in R^{N \times N}$ , and take $A = D^{- 1} V_{M} \in R^{N \times M}$ and $A^{†} = V_{M}^{T} D \in R^{M \times N}$ .
Step 7.: Compute the vectors of parameters $y_{w} = A^{†} x_{w}$ for $w = 1, \dots, W$ .

References

Evans, K. The spherical harmonic discrete ordinate method for three-dimensional atmospheric radiative transfer. J. Atmos. Sci. 1998, 55, 429–446. [Google Scholar] [CrossRef]
Spherical Harmonic Discrete Ordinate Method (SHDOM) for Atmospheric Radiative Transfer. Available online: http://coloradolinux.com/shdom/ (accessed on 11 May 2020).
Kuo, K.S.; Weger, R.; Welch, R.; Cox, S. The Picard iterative approximation to the solution of the integral equation of radiative transfer—Part II. Three-dimensional geometry. J. Quant. Spectrosc. Radiat. Transf. 1996, 55, 195–213. [Google Scholar] [CrossRef]
Goody, R.; Yung, Y. Atmospheric Radiation: Theoretical Basis; Oxford University Press: Oxford, UK, 1996. [Google Scholar]
Fu, Q.; Liou, K.N. On the Correlated k-Distribution Method for Radiative Transfer in Nonhomogeneous Atmospheres. J. Atmos. Sci. 1992, 49, 2139–2156. [Google Scholar] [CrossRef] [2.0.CO;2" target='_blank'>Green Version]
Natraj, V. A review of fast radiative transfer techniques. In Light Scattering Reviews 8; Springer: Berlin/Heidelberg, Germany, 2013; pp. 475–504. [Google Scholar] [CrossRef]
Del Águila, A.; Efremenko, D.S.; Trautmann, T. A Review of Dimensionality Reduction Techniques for Processing Hyper-Spectral Optical Signal. Light Eng. 2019, 85–98. [Google Scholar] [CrossRef] [Green Version]
Del Águila, A.; Efremenko, D.; Molina García, V.; Xu, J. Analysis of Two Dimensionality Reduction Techniques for Fast Simulation of the Spectral Radiances in the Hartley-Huggins Band. Atmosphere 2019, 10, 142. [Google Scholar] [CrossRef] [Green Version]
Natraj, V.; Jiang, X.; Shia, R.; Huang, X.; Margolis, J.; Yung, Y. Application of the principal component analysis to high spectral resolution radiative transfer: A case studyof the O2 A-band. J. Quant. Spectrosc. Radiat. Transf. 2005, 95, 539–556. [Google Scholar] [CrossRef]
Natraj, V.; Shia, R.; Yung, Y. On the use of principal component analysis to speed up radiative transfer calculations. J. Quant. Spectrosc. Radiat. Transf. 2010, 111, 810–816. [Google Scholar] [CrossRef]
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]
Goody, R.; West, R.; Chen, L.; Crisp, D. The correlated k-method for radiation calculations in nonhomogeneous atmosphere. J. Quant. Spectrosc. Radiat. Transf. 1989, 42, 539–550. [Google Scholar] [CrossRef]
Ambartzumyan, V. The effect of the absorption lines on the radiative equilibrium of the outer layers of the stars. Publ. Obs. Astron. Univ. Leningr. 1936, 6, 7–18. [Google Scholar]
Wiscombe, W.; Evans, J. Exponential-sum fitting of radiative transmission functions. J. Comput. Phys. 1997, 24, 416–444. [Google Scholar] [CrossRef]
Jolliffe, I. Principal Component Analysis; Springer: Berlin/Heidelberg, Germany, 2002. [Google Scholar] [CrossRef]
Min, W.; Lu, K.; He, X. Locality pursuit embedding. Pattern Recognit. 2004, 37, 781–788. [Google Scholar] [CrossRef]
He, X.; Niyogi, P. Locality Preserving Projections. In NIPS’03: Proceedings of the 16th International Conference on Neural Information Processing Systems; MIT Press: Cambridge, MA, USA, 2003; pp. 153–160. Available online: https://papers.nips.cc/paper/2359-locality-preserving-projections.pdf (accessed on 1 June 2020).
Fu, Y.; Huang, T. Locally Linear Embedded Eigenspace Analysis. Beckman Institute for Advanced Science and Technology, University of Illinois at Urbana-Champaign. 2005. Available online: http://citeseerx.ist.psu.edu/viewdoc/download?doi=10.1.1.93.8003&rep=rep1&type=pdf (accessed on 29 May 2020).
Pincus, R.; Evans, K.F. Computational Cost and Accuracy in Calculating Three-Dimensional Radiative Transfer: Results for New Implementations of Monte Carlo and SHDOM. J. Atmos. Sci. 2009, 66, 3131–3146. [Google Scholar] [CrossRef] [Green Version]
Alexandrov, M.; Marshak, A.; Ackerman, A. Cellular Statistical Models of Broken Cloud Fields. Part I: Theory. J. Atmos. Sci. 2010, 67, 2125–2151. [Google Scholar] [CrossRef]
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]
Nakajima, T.; Tanaka, M. Algorithms for radiative intensity calculations in moderately thick atmos using a truncation approximation. J. Quant. Spectrosc. Radiat. Transf. 1988, 40, 51–69. [Google Scholar] [CrossRef]
Efremenko, D.S.; Doicu, A.; Loyola, D.; Trautmann, T. Fast Stochastic Radiative Transfer Models for Trace Gas and Cloud Property Retrievals Under Cloudy Conditions. In Springer Series in Light Scattering; Kokhanovsky, A., Ed.; Springer: Berlin/Heidelberg, Germany, 2018; pp. 231–277. [Google Scholar] [CrossRef]
Tarasenkov, M.V.; Kirnos, I.V.; Belov, V.V. Observation of the Earth’s surface from the space through a gap in a cloud field. Atmos. Ocean. Opt. 2017, 30, 39–43. [Google Scholar] [CrossRef]
Zhuravleva, T.B.; Nasrtdinov, I.M.; Russkova, T.V. Influence of 3D cloud effects on spatial-angular characteristics of the reflected solar radiation field. Atmos. Ocean. Opt. 2017, 30, 103–110. [Google Scholar] [CrossRef]
Doicu, A.; Efremenko, D. Linearizations of the Spherical Harmonic Discrete Ordinate Method (SHDOM). Atmosphere 2019, 10, 292. [Google Scholar] [CrossRef] [Green Version]

Figure 1. Scheme of the radiative transfer problem.

Figure 2. At the lower boundary of the subdomain

k + 1

, the downward radiances are also computed at the points A, B, C, and D, which are grid points on the upper boundary of the subdomain k. During the downward iteration step, these boundary values are passed to the subdomain k.

Figure 2. At the lower boundary of the subdomain

k + 1

Figure 3. Upper panel: indicator function

f (x_{i}, y_{j})

with

x_{i} = i Δ x

and

y_{j} = j Δ y

for the

O_{2}

A-band test problem. Lower panel: a slice of the indicator function

f (x)

y = 7

km.

Figure 3. Upper panel: indicator function

f (x_{i}, y_{j})

with

x_{i} = i Δ x

and

y_{j} = j Δ y

for the

O_{2}

A-band test problem. Lower panel: a slice of the indicator function

f (x)

y = 7

km.

Figure 4. Spectral signal for the

O_{2}

A-band test problem. The differences between SHDOM with the correlated k-distribution method with and without linear embedding methods are not visible in this plot.

Figure 4. Spectral signal for the

O_{2}

A-band test problem. The differences between SHDOM with the correlated k-distribution method with and without linear embedding methods are not visible in this plot.

Figure 5. Relative differences in the spectral signal for the linear embedding methods and the

O_{2}

A-band test problem. The plots in the left panel correspond to the four-stream approximation, while the plots in the right panel correspond to the second-order of scattering approximation.

Figure 5. Relative differences in the spectral signal for the linear embedding methods and the

O_{2}

A-band test problem. The plots in the left panel correspond to the four-stream approximation, while the plots in the right panel correspond to the second-order of scattering approximation.

Figure 6. Indicator function

f (x_{i}, y_{j})

with

x_{i} = i Δ x

and

y_{j} = j Δ y

for the

{NO}_{2}

-test problem.

Figure 6. Indicator function

f (x_{i}, y_{j})

with

x_{i} = i Δ x

and

y_{j} = j Δ y

for the

{NO}_{2}

-test problem.

Figure 7. Spectral signal for the

{NO}_{2}

-test problem.

Figure 7. Spectral signal for the

{NO}_{2}

-test problem.

Figure 8. Relative differences in the spectral for the linear embedding methods and the

{NO}_{2}

-test problem. The plots in the left panel correspond to the four-stream approximation, while the plots in the right panel correspond to the second-order of scattering approximation.

Figure 8. Relative differences in the spectral for the linear embedding methods and the

{NO}_{2}

-test problem. The plots in the left panel correspond to the four-stream approximation, while the plots in the right panel correspond to the second-order of scattering approximation.

Table 1. Discretization of the domain of analysis along the vertical direction.

Subdomain	1	2	3	4	5	6	7	8
Altitude levels [km]	0–3	3–4	4–7	7–10	10–14	14–22	22–30	30–40
Discretization step [km]	0.5	0.1	0.5	0.5	1.0	2.0	2.0	5

Table 2. Root mean square (RMS) of relative difference and the computational time (CPU) in hours:minutes:seconds corresponding to the linear embedding methods for the

O_{2}

A-band test problem.

Table 2. Root mean square (RMS) of relative difference and the computational time (CPU) in hours:minutes:seconds corresponding to the linear embedding methods for the

O_{2}

A-band test problem.

Method		Four-Stream	Second-Order Scattering
PCA	RMS	$6.37 \times 10^{- 4}$	$5.19 \times 10^{- 5}$
PCA	CPU	0:04:21	0:03:01
LEA	RMS	$5.69 \times 10^{- 4}$	$4.00 \times 10^{- 5}$
LEA	CPU	0:04:45	0:03:17
LPE	RMS	$5.37 \times 10^{- 4}$	$3.29 \times 10^{- 5}$
LPE	CPU	0:04:35	0:03:05
LPP	RMS	$6.00 \times 10^{- 4}$	$3.14 \times 10^{- 5}$
LPP	CPU	0:05:01	0:03:53

Table 3. Amount of memory required for storing the radiance and source function for each subdomain.

Subdomain	1	2	3	4	5	6	7	8
Memory [MB]	549	924	549	486	228	228	228	52

Table 4. RMS of relative difference and the computational time (CPU) in hours:minutes:seconds corresponding to the linear embedding methods for the

{NO}_{2}

-test problem.

Table 4. RMS of relative difference and the computational time (CPU) in hours:minutes:seconds corresponding to the linear embedding methods for the

{NO}_{2}

-test problem.

Method		Four-Stream	Second-Order Scattering
PCA	RMS	$3.84 \times 10^{- 5}$	$1.25 \times 10^{- 4}$
PCA	CPU	1:21:24	1:02:56
LEA	RMS	$7.41 \times 10^{- 5}$	$1.34 \times 10^{- 4}$
LEA	CPU	1:22:18	1:03:48
LPE	RMS	$3.73 \times 10^{- 5}$	$1.12 \times 10^{- 4}$
LPE	CPU	1:21:53	1:03:21
LPP	RMS	$7.22 \times 10^{- 5}$	$1.26 \times 10^{- 4}$
LPP	CPU	1:22:43	1:04:09

Publisher’s Note: MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

© 2020 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).

Share and Cite

MDPI and ACS Style

Doicu, A.; Efremenko, D.S.; Trautmann, T. A Spectral Acceleration Approach for the Spherical Harmonics Discrete Ordinate Method. Remote Sens. 2020, 12, 3703. https://doi.org/10.3390/rs12223703

AMA Style

Doicu A, Efremenko DS, Trautmann T. A Spectral Acceleration Approach for the Spherical Harmonics Discrete Ordinate Method. Remote Sensing. 2020; 12(22):3703. https://doi.org/10.3390/rs12223703

Chicago/Turabian Style

Doicu, Adrian, Dmitry S. Efremenko, and Thomas Trautmann. 2020. "A Spectral Acceleration Approach for the Spherical Harmonics Discrete Ordinate Method" Remote Sensing 12, no. 22: 3703. https://doi.org/10.3390/rs12223703

APA Style

Doicu, A., Efremenko, D. S., & Trautmann, T. (2020). A Spectral Acceleration Approach for the Spherical Harmonics Discrete Ordinate Method. Remote Sensing, 12(22), 3703. https://doi.org/10.3390/rs12223703

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

Article Menu

A Spectral Acceleration Approach for the Spherical Harmonics Discrete Ordinate Method

Abstract

1. Introduction

2. The Spectral Acceleration Approach

2.1. Correlated k-Distribution Method

2.2. Dimensionality Reduction of Atmospheric Optical Parameters

3. Numerical Analysis

3.1. SHDOM Implementation

3.2. Numerical Results

3.2.1. Oxygen A-Band Test Problem

3.2.2. ${NO}_{2}$ -Test Problem

4. Conclusions

Author Contributions

Funding

Conflicts of Interest

Appendix A

References

Share and Cite

Article Metrics

Article Access Statistics

Further Information

Guidelines

MDPI Initiatives

Follow MDPI

Article Menu

A Spectral Acceleration Approach for the Spherical Harmonics Discrete Ordinate Method

Abstract

1. Introduction

2. The Spectral Acceleration Approach

2.1. Correlated k-Distribution Method

2.2. Dimensionality Reduction of Atmospheric Optical Parameters

3. Numerical Analysis

3.1. SHDOM Implementation

3.2. Numerical Results

3.2.1. Oxygen A-Band Test Problem

3.2.2. NO 2 -Test Problem

4. Conclusions

Author Contributions

Funding

Conflicts of Interest

Appendix A

References

Share and Cite

Article Metrics

Article Access Statistics

Further Information

Guidelines

MDPI Initiatives

Follow MDPI

3.2.2. ${NO}_{2}$ -Test Problem