Cluster Low-Streams Regression Method for Hyperspectral Radiative Transfer Computations: Cases of O₂ A- and CO₂ Bands

Abstract

Current atmospheric composition sensors provide a large amount of high spectral resolution data. The accurate processing of this data employs time-consuming line-by-line (LBL) radiative transfer models (RTMs). In this paper, we describe a method to accelerate hyperspectral radiative transfer models based on the clustering of the spectral radiances computed with a low-stream RTM and the regression analysis performed for the low-stream and multi-stream RTMs within each cluster. This approach, which we refer to as the Cluster Low-Streams Regression (CLSR) method, is applied for computing the radiance spectra in the O₂ A-band at 760 nm and the CO₂ band at 1610 nm for five atmospheric scenarios. The CLSR method is also compared with the principal component analysis (PCA)-based RTM, showing an improvement in terms of accuracy and computational performance over PCA-based RTMs. As low-stream models, the two-stream and the single-scattering RTMs are considered. We show that the error of this approach is modulated by the optical thickness of the atmosphere. Nevertheless, the CLSR method provides a performance enhancement of almost two orders of magnitude compared to the LBL model, while the error of the technique is below 0.1% for both bands.

1. Introduction

Radiative transfer models (RTM) are a key part of the remote sensing algorithms, which are used to retrieve atmospheric parameters from Earth observation data. The new atmospheric composition sensors with high spatial and spectral resolution require accurate and efficient computational algorithms. An accurate method to simulate spectral radiances in the absorption bands is the line-by-line (LBL) approach [1]. However, due to the high spectral variability of the gas absorption coefficient k in the absorption bands, the LBL-approach is very time-consuming because it requires up to several thousands of monochromatic computations per absorption band. In this regard, specifically designed hyperspectral RTMs are required. They reduce the number of monochromatic computations considerably without compromising accuracy. The basic idea behind the hyperspectral RTMs is to get rid of redundancies in hyperspectral data. This principle can be traced back to Ambartsumian [2], who noted that the transmission within a spectral interval does not depend on the LBL variation of k, but rather on the distribution of values of k within the spectral interval. This idea was used in the k-distribution method [3,4,5,6], in which the wavelengths are grouped into a smaller number of bins with similar values of k, and the LBL calculations are replaced by a smaller set of radiative transfer simulations. In [7,8], the transmission function for a given spectral interval is fitted by a sum of exponentials, while the corresponding fitting coefficients are computed from a reduced number of monochromatic computations. A similar approach is described in Moncet et al. [9], where the fitting weights and the most representative wavelengths are chosen appropriately.

The state-of-the-art hyperspectral RTMs employ dimensionality reduction techniques such as principal component analysis (PCA). In [10,11], PCA is applied to the spectral radiance data to establish a set of empirical orthogonal functions (EOFs), so that an arbitrary spectrum at full spectral resolution can be reconstructed as a weighted sum of EOFs. The weights are found by performing monochromatic simulations at a reduced number of wavelengths. To accelerate the computations in the O₂ A-band, Natraj et al. [12] proposed a fundamentally different PCA-based radiative transfer model, in which the dimensionality of the optical properties data is reduced. A two-stream radiative transfer model was used as an approximate model, and the dependency of the corresponding correction factor on the optical parameters was modeled by a second-order Taylor expansion about the mean value of the optical parameters in the reduced optical data space. This approach was extended to other dimensionality reduction techniques [13] and spectral ranges [14,15,16]; moreover, it was implemented in conjunction with PCA for spectral radiances [17] and with the k-distribution method [18]. The errors of these approaches are usually below 0.1% for the spectral radiances, while the performance enhancement may reach several orders of magnitude depending on the spectral region and the required level of accuracy.

In Efremenko et al. [19] it was shown that, after parallelizing the PCA-based RTM computations, as much as half of the computational time is due to the PCA itself; consequently, according to Amdahl’s Law [20], no further acceleration of a PCA-based RTM is possible. In this regard, it is of great interest to develop acceleration techniques for hyperspectral RTM computations that do not exploit PCA. A possible candidate, which deserves this purpose, is the low streams interpolation method proposed by O’Dell [21]. In this method, (i) each spectral point is assigned to a specific bin according to the values of the gas optical depth, (ii) an approximate two-stream model is used to compute the radiances at all spectral points, and (iii) the two-stream radiances are corrected to the multi-stream radiances by calling a multi-stream model only once for each bin. In this paper, we propose the so-called Cluster Low-Streams Regression (CLSR) method, in which (i) wavelengths are grouped according to the radiance values obtained with the approximate model (rather than the optical properties), and (ii) the single-scattering and the two-stream models are used as approximate models. The CLSR method is tested against the PCA-based RTM for the O₂ A-band at 760 nm and the weak CO₂ band at 1610 nm. These spectral bands are of great interest for retrieving aerosol and cloud parameters [22], as well as CO₂ concentrations [23].

The paper is organized as follows. In Section 2, the methodology is presented, including the descriptions of the PCA-based RTM adopted in our research and the CLSR method. In Section 3, the accuracy and efficiency of the CLSR method for several atmospheric scenarios in the O₂ A- and weak CO₂ absorption bands are analyzed and compared with the PCA-based RTM. In addition, the accuracy of the radiance spectra convolved with the slit functions corresponding to GOME-2 [24], TROPOMI on board Sentinel 5-P [25] and GOSAT [26] instruments is examined. The paper concludes with a summary.

2. Methodology

2.1. Reference Radiative Transfer Model

As a reference RTM, the scalar method of discrete ordinates with matrix exponential (DOME) [27] is used. The model is supplied with the left eigenvector approach [28,29,30] based on the scaling transformation [31,32]. In this method, the number of discrete ordinates (streams) in the polar hemisphere $N_{\mathrm{do}}$ controls the computational time and accuracy. The model is called multi-stream (MS) when $N_{\mathrm{do}}\ge 2$ and low-stream (LS) otherwise. Specifically, the case $N_{\mathrm{do}}=1$ is the two-stream model. The simplest RTM is the single-scattering model, which solves the radiative transfer equation neglecting the integral term [33].

The gaseous absorptions for the O₂ A- and CO₂ bands are computed with the LBL model Py4CAtS [34], while the gas absorption cross-sections are taken from the HITRAN 2016 database [35]. The uncertainties in the spectroscopic parameters are not considered because their role is irrelevant for this study [36]. Continuum (also referred to as collision-induced absorption (CIA) [37,38]) contributions to molecular absorption are not taken into account. We note that the CIA process gives a broad and smooth contribution [39] and hence, does not impose difficulties for the regression techniques considered in this study.

2.2. Study Cases

We consider the reflected spectral radiance at the top of the atmosphere (TOA). In our simulations, the atmosphere is discretized into 35 layers with a step of 1 km between 0 and 25 km, and a step of 2.5 km between 25 km and 50 km. We assume a Lambertian surface with an albedo of 0.3. The solar zenith angle, the viewing zenith angle and the relative azimuth angle are 45°, 35° and 90°, respectively.

In the O₂ A-band the spectral sampling is 0.001 nm in the spectral range 755–775 nm, while in the CO₂ band the spectral sampling is 0.0015 nm in the spectral interval 1590–1620 nm. Thus, on each band, 20,000 spectral points are considered. The Rayleigh cross-sections and depolarization ratios are computed as in [40], while the pressure and temperature profiles correspond to the US standard model atmosphere [41]. The computations are performed for the unit solar irradiance at the TOA.

We define 5 atmospheric scenarios: ‘Clear sky’, ‘Aerosol 1’, ‘Aerosol 2’, ‘Cloud 1’ and ‘Cloud 2’. The ‘Clear sky’ scenario corresponds to an atmosphere without clouds and aerosols. In the ‘Aerosol 1’ and ‘Aerosol 2’ scenarios, the atmosphere contains the clean continental and the polluted continental aerosols taken from the OPAC database [42], respectively. The aerosol optical depths at the middle of the absorption bands are shown in

Table 1. Aerosol optical thickness used in the simulations for O₂ A- and CO₂ bands at the middle of the corresponding absorption band.

	O2 A-Band (760 nm)	CO2 Band (1610 nm)
Aerosol 1	0.2	0.08
Aerosol 2	1.2	0.41

. In the ‘Cloud 1’ and ‘Cloud 2’ scenarios, a continental clean cumulus cloud with a modified Gamma size distribution [42,43]:

$$p\left(a\right)\propto a^{\alpha }exp\left[-\frac{\alpha }{\gamma }{\left(\frac{a}{a_{\mathrm{mod}}}\right)}^{\gamma }\right],$$

is considered. The size distribution parameters are $a_{\mathrm{mod}}=4.8$ µm, $\alpha =5$ and $\gamma =2.16$, the droplet size ranges between 0.02 and 50.0 µm, and the cloud optical depths are $\tau =10$ (‘Cloud 1’) and $\tau =20$ (‘Cloud 2’). In both cases, the cloud-top height is 5 km and the cloud geometrical thickness is 1 km.

The number of streams considered for the multi-stream model is the same for all scenarios. We might need a different number of streams for different scenarios, but in order to simplify and ensure convergence of relative errors for all cases, we use 32 streams.

2.3. PCA-Based RTM

Since a PCA-based RTM can be implemented in several ways (e.g., [12,13,15,16,17,44,45]), we describe below the main features of the PCA-based RTM used in this paper. The idea of the method is to refine the two-stream solution by a correction function estimated in the reduced space of the optical properties. Considering a discretization of the atmosphere in L layers, we define a $2L$-dimensional vector ${\mathbf{x}}_w$ for each wavelength ${\left\{{\lambda }_w\right\}}_{w=1}^W$, by

$${\mathbf{x}}_w^{\mathrm{T}}=\left[ln{\tau }_1\left({\lambda }_w\right),\dots ,ln{\tau }_L\left({\lambda }_w\right),ln{\omega }_1\left({\lambda }_w\right),\dots ,ln{\omega }_L\left({\lambda }_w\right)\right],$$

(1)

where ${\tau }_k$ and ${\omega }_k$ are the optical thickness and the single-scattering albedo of the k-th layer, respectively. Thus, the vector ${\mathbf{x}}_w$ encapsulates the wavelength variability of the optical parameters, which are the input parameters of the radiative transfer code. By applying PCA, a M-dimensional subspace is found which is spanned by a set of linear independent vectors (empirical orthogonal functions) ${\left\{{\mathbf{q}}_k\right\}}_{k=1}^M$ such that the centered (mean-removed) data ${\mathbf{x}}_w-\overline{\mathbf{x}}$ lie mainly on this subspace, i.e.,

$${\mathbf{x}}_w\approx \overline{\mathbf{x}}+\sum _{k=1}^M{y}_{wk}{\mathbf{q}}_k,w=1,\dots ,W,$$

(2)

where $y_{wk}$ are the principal component (PC) scores.

Following Natraj et al. [12], we introduce a correction function

$$f\left({\mathbf{x}}_w\right)=ln\frac{I_{\mathrm{MS}}\left({\mathbf{x}}_w\right)}{I_{\mathrm{LS}}\left({\mathbf{x}}_w\right)},$$

(3)

where $I_{\mathrm{LS}}\left({\mathbf{x}}_w\right)$ is the radiance provided by the low-stream RTM, and $I_{\mathrm{MS}}\left({\mathbf{x}}_w\right)$ is the radiance simulated by the multi-stream model. Setting $\mathsf{\Delta }{\mathbf{x}}_w={\sum }_{k=1}^M{y}_{wk}{\mathbf{q}}_k,$ we consider a second-order Taylor expansion of $f\left({\mathbf{x}}_w\right)$ around $\overline{\mathbf{x}}$ in the direction $\mathsf{\Delta }{\mathbf{x}}_w$. Approximating the directional derivatives with central differences (see [13,17] for mathematical details) we obtain

$$\begin{array}{cc}\hfill f\left({\mathbf{x}}_w\right)& \approx f\left(\overline{\mathbf{x}}\right)+\frac{1}{2}\sum _{k=1}^M[f(\overline{\mathbf{x}}+{\mathbf{q}}_k)-f(\overline{\mathbf{x}}-{\mathbf{q}}_k)]y_{wk}+\frac{1}{2}\sum _{k=1}^M[f(\overline{\mathbf{x}}+{\mathbf{q}}_k)-2f\left(\overline{\mathbf{x}}\right)+f(\overline{\mathbf{x}}-{\mathbf{q}}_k)]y_{wk}^2.\hfill \end{array}$$

(4)

Please note, since $M<2L$ (and in practice, $M\ll 2L$), the estimation of the correction function in the reduced data space is much faster than in the original space. Also note that the total number of calls of the multi-stream RTM is $2M+1$. Once $f\left({\mathbf{x}}_w\right)$ is computed, the low-stream radiances can be converted into the multi-stream radiances by using Equation (3). Finally, the predicted multi-stream radiance ${\tilde{I}}_{\mathrm{MS}}$ can be expressed as follows:

$${\tilde{I}}_{\mathrm{MS}}=I_{\mathrm{LS}}exp\left[f\left({\mathbf{x}}_w\right)\right].$$

(5)

In del Águila et al. [17], it was shown that a higher order expansion of $f\left({\mathbf{x}}_w\right)$ in Equation (4) does not substantially improve the accuracy of the solution. In our simulations, an appropriate value for M is found iteratively, i.e., M is increased until the approximation error is sufficiently small.

2.4. Regression Relationship between Multi-Stream and Low-Stream Models

The LBL computations can be accelerated by finding a regression relationship between the multi-stream radiances and the optical parameters. In O’Dell [21], the correction function of the low-stream radiance as a function of the gas optical thickness ${\tau }_{\mathrm{gas}}$, i.e., the relative error

$$f\left({\tau }_{\mathrm{gas}}\right)=\frac{I_{\mathrm{LS}}\left({\tau }_{\mathrm{gas}}\right)-I_{\mathrm{MS}}\left({\tau }_{\mathrm{gas}}\right)}{I_{\mathrm{MS}}\left({\tau }_{\mathrm{gas}}\right)},$$

was analyzed (Figure 1). The idea is to compute $f\left({\tau }_{\mathrm{gas}}\right)$ at certain values of ${\tau }_{\mathrm{gas}}$, and then to use a regression model for computing $f\left({\tau }_{\mathrm{gas}}\right)$ at the remaining values of ${\tau }_{\mathrm{gas}}.$ However, for the present application, such an approach has two major drawbacks:

1.

the dependence of $f\left({\tau }_{\mathrm{gas}}\right)$ on ${\tau }_{\mathrm{gas}}$ is non-linear and, therefore, the application of the regression model requires a binning of the ${\tau }_{\mathrm{gas}}$ values;

2.

in a mathematical sense, $f\left({\tau }_{\mathrm{gas}}\right)$ is not a function (for a value of ${\tau }_{\mathrm{gas}}$, there are several values of $f\left({\tau }_{\mathrm{gas}}\right)$) and therefore, even for a fine binning, the regression model will be not accurate.

To overcome these drawbacks, O’Dell [21] devised a regression model based on the so-called adjusted gas optical thickness, defined as the gas optical thickness from the TOA down to the layer in the atmosphere where the scattering optical depth equals some critical value.

A similar idea, which we exploit in the following, is to use a regression model for the radiances computed by a low-stream RTM rather than in the ${\tau }_{\mathrm{gas}}$-space. Figure 2 shows a comparison between the multi- and low-stream radiances for the ‘Aerosol 2’ scenario. Although the errors of the low-stream RTM are significant (e.g., they may reach 200%), the low- and multi-stream radiances have a similar spectral behaviour. Moreover, as shown in Figure 3, the dependence between the low- and multi-stream radiances is almost linear. In this respect, it seems reasonable to cluster the spectral points according to the radiance values computed with a low-stream RTM (rather than according to the optical properties) and to establish a regression model between the low- and multi-stream radiances within each cluster. In order to decrease the errors of the regression estimates, an additional wavelength-dependent parameter will be included in the regression model: the direct transmittance T.

2.5. Cluster Low-Streams Regression Method

The Cluster Low-Streams Regression (CLSR) method can be formulated as follows.

1.

Consider a high-resolution spectrum ${\left\{I_{\mathrm{LS}}\left({\lambda }_i\right)\right\}}_{i=1}^N$ computed at N spectral points ${\left\{{\lambda }_i\right\}}_{i=1}^N$ by means of a low-stream RTM.

2.

Sort the radiance set ${\left\{I_{\mathrm{LS}}\left({\lambda }_i\right)\right\}}_{i=1}^N$ in ascending order, and let ${\left\{{\widehat{I}}_{\mathrm{LS},i}\right\}}_{i=1}^N$, with ${\widehat{I}}_{\mathrm{LS},i}\le {\widehat{I}}_{\mathrm{LS},i+1}$, be the sorted radiance set (Figure 4a).

3.

Consider C clusters in ${\left\{{\widehat{I}}_{\mathrm{LS},i}\right\}}_{i=1}^N$ with $N_C=N/C$ radiance points (Figure 4b), and let the c cluster be defined by the radiance set ${\left\{{\widehat{I}}_{\mathrm{LS},i}^c\right\}}_{i=1}^{N_C}$ for $c=1,\dots ,C$.

4.

Select n equidistant radiance points in the c cluster, i.e., ${\left\{{\overline{I}}_{\mathrm{LS},q}^c\right\}}_{q=1}^n$, and for the corresponding wavelengths compute the multi-stream radiances ${\left\{{\overline{I}}_{\mathrm{MS},q}^c\right\}}_{q=1}^n$ (Figure 4c).

5.

Assume that in each cluster c we have the linear relationship

$${\widehat{I}}_{\mathrm{MS},i}^c={\alpha }^c{\widehat{T}}_i^c+{\beta }^c{\widehat{I}}_{\mathrm{LS},i}^c+{\gamma }^c,i=1,\dots ,N_C,$$

(6)

where ${\alpha }^c$, ${\beta }^c$ and ${\gamma }^c$ are the regression coefficients of the c-th cluster, and $\widehat{T}$ is the corresponding direct transmittance.

6.

Compute the regression coefficients ${\alpha }^c$, ${\beta }^c$ and ${\gamma }^c$ as a solution to the least square problem

$$\left({\alpha }^c,{\beta }^c,{\gamma }^c\right)=arg\underset{{\alpha }^c,{\beta }^c,{\gamma }^c}{min}\sum _{q=1}^n{\left[{\overline{I}}_{\mathrm{MS},q}^c-\left({\alpha }^c{\overline{T}}_q^c+{\beta }^c{\overline{I}}_{\mathrm{LS},q}^c+{\gamma }^c\right)\right]}^2.$$

(7)

7.

Use the values of $({\alpha }^c,{\beta }^c,{\gamma }^c)$ found in the previous step to restore the multi-stream radiances ${\left\{{\tilde{I}}_{\mathrm{MS},i}\right\}}_{i=1}^N$ for all the spectral points according to Equation (6) (Figure 4d).

Here, the “hat” notation $\widehat{I}$ refers to the sorted radiances, the “bar” notation $\overline{I}$ refers to the equidistant radiances entering the regression model, and the “tilde” notation $\tilde{I}$ refers to the predicted radiances. Please note that the total number of regression points, and thus the number of calls of the multi-stream model, is $nC$.

2.6. Efficiency and Computational Performance Estimation

To estimate the accuracy of the acceleration techniques, we consider the residual error for the radiance at each spectral point ${\lambda }_i$:

$$\mathsf{\Delta }{I}_{\mathrm{res},i}=\frac{{\tilde{I}}_{\mathrm{MS},i}-I_{\mathrm{MS},i}}{I_{\mathrm{MS}}^{\mathrm{cont}}}·100,$$

(8)

where ${\tilde{I}}_{\mathrm{MS},i}$ is the predicted radiance calculated with either the PCA-based RTM (cf. Equation (5)) or the CLSR method (cf. Equation (6)), while $I_{\mathrm{MS}}^{\mathrm{cont}}$ is the radiance without absorption (i.e., the continuum radiance, which is used to avoid radiance values close to zero in the denominator of Equation (8), when strong gas absorption is present [15]). The mean relative error is computed as follows:

$$\epsilon =\frac{{\sum }_{i=1}^N\left|\mathsf{\Delta }{I}_{\mathrm{res},i}\right|}{N}.$$

(9)

Additionally, the median and interquartile range (IQR) of the residual are computed. Essentially, these parameters serve as robust metrics to evaluate the variability of the residuals in the corresponding spectral band, thereby with a lower sensitivity to single outliers [16].

To estimate the performance enhancement, we define the speedup factor as the ratio between the computational time of a multi-stream LBL calculation to that of a certain acceleration technique. To exclude the hardware-related factors from our analysis, we estimate the speedup factor for the PCA-based RTM as [16]

$$S_{\mathrm{PCA}}=\frac{t_{\mathrm{MS}}\times N}{t_{\mathrm{LS}}\times N+t_{\mathrm{PCA}}+(2M+1)(t_{\mathrm{MS}}+t_{\mathrm{LS}})},$$

(10)

where $t_{\mathrm{MS}}$ and $t_{\mathrm{LS}}$ are the computational times for a single monochromatic calculation corresponding to the multi- and low-stream RTMs, respectively, while $t_{\mathrm{PCA}}$ is the computational time of the PCA. Thus, the numerator in Equation (10) is the time required for the LBL computations. Please note the PCA-based RTM requires $2M+1$ calls to the multi-stream model.

The speedup factor for the CLSR method is defined by

$$S_{\mathrm{CLSR}}=\frac{t_{\mathrm{MS}}\times N}{t_{\mathrm{LS}}\times N+t_{\mathrm{LSM}}\times C+t_{\mathrm{MS}}\times n\times C},$$

(11)

where $t_{\mathrm{LSM}}$ is the computational time needed for the least-squares method according to Equation (7). Recall that the total number of regression points and the corresponding calls of the multi-stream model is $nC$.

3. Results and Discussion

3.1. Performance of the Low- and Multi-Streams Models

In this section, we assess the performance of the low- and multi-stream RTMs before applying the acceleration techniques. The spectral radiances for three atmospheric scenarios (‘Clear sky’, ‘Cloud 1’ and ‘Cloud 2’) are illustrated in Figure 5.

Table 2. Mean relative error $\epsilon$ for the low-stream models (single-scattering and two-stream models) compared with the multi-stream model for the O₂ A-band.

Scenario	Single-Scattering Model	Two-Stream Model
	$\mathit{\epsilon }$ (%)	$\mathit{\epsilon }$ (%)
Clear sky	3.4	0.14
Aerosol 1	36	1.0
Aerosol 2	78	4.6
Cloud 1	93	7.1
Cloud 2	93	1.4

and

Table 3. Mean relative error $\epsilon$ for the low-stream models (single-scattering and two-stream models) compared with the multi-stream model for the CO₂ weak band.

Scenario	Single-Scattering Model	Two-Stream Model
	$\mathit{\epsilon }$ (%)	$\mathit{\epsilon }$ (%)
Clear sky	0.203	0.004
Aerosol 1	14.7	0.40
Aerosol 2	55.5	2.1
Cloud 1	97.56	6.13
Cloud 2	97.91	2.00

provide the mean relative errors of the low-stream RTMs with respect to the multi-stream RTM for the O₂ A- and CO₂ absorption bands, respectively.

Table 4. Computational time in seconds of the radiative transfer solution as a function of the number of streams, and speedup factors with respect to the case $N_{\mathrm{do}}=32$.

${\mathit{N}}_{\mathbf{do}}$	Time (s)	Speedup Factor
0	0.4	5800
1	3.2	725
2	12.4	187.1
4	34.4	67.4
8	110	21
16	550	4.2
32	2320	-

shows the computational performance for a different number of discrete ordinates $N_{\mathrm{do}}$. The results show that (i) the relative errors of the single-scattering model are higher than those of the two-stream model, (ii) the relative error increases as the optical thickness of the atmosphere increases, and (iii) the speedup factors of the single-scattering model are considerably higher than those of the two-stream model.

3.2. Acceleration Techniques: Accuracy Evaluation

3.2.1. Spectral Residuals

The residuals (cf. Equation (8)) are computed for (i) different numbers of principal components in the case of the PCA-based RTM, and (ii) different numbers of regression points per cluster in the case of the CLSR method. Figure 6 shows these residuals for the ‘Cloud 2’ and the ‘Aerosol 2’ scenarios.

From Figure 6 the following conclusions can be drawn.

1.

The residuals decrease when increasing the number of PCs and regression points.

2.

The interquartile range for the CLSR method is substantially reduced when switching from 1–2 to 3 regression points for the CO₂ band, and from 3 to 4 regression points for the O₂ A-band.

3.

In both spectral bands, the interquartile range for the PCA-based RTM decreases gradually with the number of principal components.

4.

The residuals in the O₂ A-band are systematically higher than those in the CO₂ band. This behaviour is more pronounced when the gas optical thickness is large, thus resulting in larger discrepancies for the PCA-based RTM and almost negligible discrepancies for the CLSR approach.

5.

In the O₂ A-band, the median values of the residuals Equation (8) are higher than 2%. The median values remain almost constant with the number of principal components and regression points. However, the median values as well as the interquartile range for the PCA method are generally higher than those of the CLSR method. This trend remains coherent using the single-scattering RTM instead of the two-stream RTM, although the residuals are substantially higher for the PCA-based RTM.

3.2.2. Estimation of the Required Parameters for the Acceleration Techniques

In order to select the optimal number of PCs, we compute the explained variance ratio of the optical data, which accounts for the variance associated with a given number of PCs. The results in Figure 7 show that in all scenarios almost 99% of the optical data variance can be explained within the first four PCs. Moreover, from Figure 6a we can infer that (i) the convergence of the two-stream radiances is reached for 4-5 PCs, while (ii) the convergence of the single-scattering radiances requires 4-7 PCs (Figure 6b).

To select the number of clusters and regression points for the CLSR method, we estimate the mean error (cf. Equation (9)). As an example, we illustrate in Figure 8 the mean errors for the ‘Aerosol 2’ scenario. The results show that 4–5 clusters and 3–5 regression points guarantee a small mean error.

3.2.3. Accuracy of the CLSR and PCA-Based Methods

In this section, we assess the accuracy of the PCA-based and CLSR methods. The residuals (Equation (8)) and interquartile ranges corresponding to the PCA-based and CLSR methods, as well as for the two-stream and single-scattering models, are shown in Figure 9.

The results can be summarized as follows:

1.

For both low-stream models and cloudy scenarios, (i) the residuals of the PCA-based method in the O₂ A-band are substantially larger than those of the CLSR method (on average, for the two-stream model in the O₂ A-band, the residuals are above 1% for the PCA-based method and below 0.01% for the CLSR method), while (ii) the PCA-based and CLSR methods yield comparable accuracies in the CO₂ band (on average, for the two-stream model in the CO₂ band, the residuals are below 0.1% for the PCA-based method and below 0.01% for the CLSR method). Please note similar results were established in previous studies, i.e., the errors of the PCA-based method generally increase with the optical depth [16]. The superiority of the CLSR method is also demonstrated by the interquartile ranges: these are much smaller for the CLSR method.

2.

For both O₂ A- and CO₂ bands, (i) the residuals of the single-scattering model with the PCA-based method are higher than those corresponding to the two-stream model, while (ii) the residuals of the two-stream and single-scattering models with the CLSR method are comparable (on average, for the single-scattering model, the residuals of the CLSR method are below 0.2% in the O₂ A-band and below 0.1% in CO₂ band, while for the two-stream model the corresponding errors are below 0.01% in both spectral bands).

Thus, (i) the two-stream model with the PCA-based and CLSR methods yields accurate results, (ii) the efficiency of the PCA-based method decreases when increasing the optical thickness, and (iii) both the two-stream and the single-scattering models with the CLSR method provide reasonable accuracies.

Please note the considered atmospheric scenarios implicitly involve a set of assumptions which causes differences between RTM calculations and actual (physically realistic) atmospheric conditions. Furthermore, the discrete ordinate method itself is an approximate technique which in principle should be validated against rigorous analytical methods [46]. However, we use the same atmospheric definition inputs for comparing acceleration techniques.

3.3. Acceleration Techniques: Computational Performance

In this section, we evaluate the computational efficiency of the CLSR method and compare it with the PCA-based method. The speedup factors computed by using Equations (10) and (11) for the PCA-based and CLSR methods, respectively, are given in

Table 5. Speedup factor of the PCA-based ($S_{\mathrm{PCA}}$) and CLSR methods with: the two-stream model ($S_{\mathrm{CLSR}}^{\mathrm{TS}}$) and single-scattering model ($S_{\mathrm{CLSR}}^{\mathrm{SS}}$).

${\mathit{S}}_{\mathbf{PCA}}$	${\mathit{S}}_{\mathbf{CLSR}}^{\mathbf{TS}}$	${\mathit{S}}_{\mathbf{CLSR}}^{\mathbf{SS}}$
534	505	1294

. Since the single-scattering model with the PCA-based method provides large errors (see Figure 9b), the table includes the speedup factor $S_{\mathrm{PCA}}$ for the two-steam model with the PCA-based method, as well as for the two-stream and single-scattering models with the CLRS method, $S_{\mathrm{CLSR}}^{\mathrm{TS}}$ and $S_{\mathrm{CLSR}}^{\mathrm{SS}}$, respectively.

According to the optimal values of the PCs and regression points and clusters, the number of calls to the multi-stream model for the PCA-based method is 9, for the two-stream model with the CLSR method 12, and for the single-scattering model with the CLSR method 20. The computational time to perform the PCA calculation $t_{\mathrm{PCA}}$ is 3 orders of magnitude higher than that of the least-squares calculation $t_{\mathrm{LSM}}$. However, since the most time-consuming part is due to the number of calls to the multi-stream model, the efficiencies of the two-stream model with the CLRS and the PCA-based method are comparable. In contrast, the computational performance of the single-scattering model with the CLRS method is much higher than that with the PCA-based method due to the neglect of multiple scattering computations. Comparing the speedup factors of this study with those of other authors (e.g., [21]), we find that our values are of the order of their speedup factors and one order of magnitude higher when considering the single-scattering model for the CLSR technique.

3.4. Computation of Convolved Spectra

In this section, we examine the efficiency of the CLSR method for computing convolved spectra. The high-resolution spectra in the O₂ A-band are convolved with the slit functions corresponding to GOME-2 and TROPOMI instruments, while the radiance spectra in the CO₂ band are convolved with the GOSAT slit function. In this paper, slit functions are modeled with a Gaussian function. The corresponding full widths at half maximum (FWHM) are listed in

Table 6. Spectral ranges and FWHM of the Gaussian slit functions of the instruments used in the study: TROPOMI, GOME-2 and GOSAT.

Instrument	Spectral Range	FWHM
TROPOMI	710–775 nm	0.183 nm
GOME-2	590–790 nm	0.51 nm
GOSAT	1.56–1.69 µm	0.2 cm−1

. The FWHM considered for the O₂ A-band are based on pre-launch calibrations [47] and for the CO₂ band on [48].

The examples of convolved spectra corresponding to the ‘Cloud 1’ scenario are shown in Figure 10. The computations are performed by using the PCA-based and the CLSR methods (in conjunction with the two-stream model) as well as the LBL approach.

Table 7. Mean relative error $\epsilon$ for the convolved spectra compared with the multi-stream spectra for the PCA-based and CLSR methods and for the different atmospheric scenarios considered for the O₂ A-band. In all cases, the low-stream model considered is the two-stream model. The instruments analyzed are GOME-2 and TROPOMI, which are compared with the non-convolved values.

Scenario	O2 A-Band
	GOME-2		TROPOMI		Non-Convolved
	${\mathbf{\epsilon }}_{\mathbf{PCA}}(\%)$	${\mathbf{\epsilon }}_{\mathbf{CLSR}}(\%)$	${\mathbf{\epsilon }}_{\mathbf{PCA}}(\%)$	${\mathbf{\epsilon }}_{\mathbf{CLSR}}(\%)$	${\mathbf{\epsilon }}_{\mathbf{PCA}}(\%)$	${\mathbf{\epsilon }}_{\mathbf{CLSR}}(\%)$
Clear sky	0.021	0.004	0.021	0.004	0.022	0.006
Aerosol 1	0.049	0.007	0.050	0.007	0.061	0.011
Aerosol 2	0.837	0.019	0.837	0.019	0.856	0.026
Cloud 1	1.23	0.011	1.23	0.011	1.24	0.017
Cloud 2	2.92	0.006	2.92	0.006	2.93	0.009

and

Table 8. Mean relative error $\epsilon$ for the convolved spectra compared with the multi-stream spectra for the PCA-based and CLSR methods and for the different atmospheric scenarios considered for the CO₂ band. In all cases, the low-stream model considered is the two-stream model. The instrument analyzed is GOSAT, which is compared with the non-convolved values.

Scenario	CO2 Band
	GOSAT		Non-Convolved
	${\mathbf{\epsilon }}_{\mathbf{PCA}}(\%)$	${\mathbf{\epsilon }}_{\mathbf{CLSR}}(\%)$	${\mathbf{\epsilon }}_{\mathbf{PCA}}(\%)$	${\mathbf{\epsilon }}_{\mathbf{CLSR}}(\%)$
Clear sky	0.0003	0.0003	0.0006	0.0006
Aerosol 1	0.0075	0.0067	0.010	0.012
Aerosol 2	0.035	0.012	0.044	0.017
Cloud 1	0.011	0.007	0.013	0.008
Cloud 2	0.011	0.005	0.016	0.006

show the mean relative errors for the PCA-based (${\epsilon }_{\mathrm{PCA}}$) and CLSR (${\epsilon }_{\mathrm{CLSR}}$) methods for the O₂ A-band and the CO₂ band, respectively. In addition, the residuals for non-convolved spectra are shown for comparison. For the ‘Clear sky’ and aerosol scenarios the accuracies of both methods are comparable, while for cloud scenarios the CLSR method is more accurate. Please note the residuals estimated for the convolved spectra are very close to those for the non-convolved ones and the value of residuals according to Equation (8) are robust.

4. Conclusions

In this study, we developed the Cluster Low-Streams Regression (CLSR) method for fast radiative transfer simulations of the O₂ A- and CO₂ absorption bands. The CLSR method exploits a strong close-to-linear relationship between the radiances computed with a low-stream model (which is either the two-stream or the single-scattering model) and the radiances computed with the multi-stream model. The spectral points are grouped in several clusters according to the values of the low-stream radiances. For each cluster, the regression model is established between the low-stream and multi-stream models, where the corresponding regression coefficients are found by using the least-squares method. This approach can be regarded as a variation of the low-streams interpolation method explained in [21], in which the binning is performed in the space of gas optical depths.

The CLRS method was compared with the PCA-based RTM in several atmospheric scenarios. For the ‘Clear sky’ scenario, the performance of both acceleration methods is comparable. For the cloud scenarios, the CLSR method shows more accurate results than the PCA-based method. However, to improve the performance of the PCA-based method, in [14,15,16] the optical depth binning technique was applied, i.e., the spectral points were grouped into bins according to the optical depth values and the PCA-based method was applied independently for each bin. Such an approach improves the accuracy of the PCA-based method, although at the cost of increasing the number of multi-stream computations. Comparing the results of [14,15,16] with the results obtained in our study, we can conclude that the PCA-based binned approaches and the CLSR method are comparable in terms of accuracy.

Our analysis has shown that, although the CLSR method requires more calls to the multi-stream model than the PCA-based model, the corresponding speedup factors are very similar, with slightly better results for the CLSR method. However, this difference is not significant. In addition, we note that the CLSR method does not require to deal with the eigenvalue problem as the PCA-based RTMs do. The CLSR method can be used in conjunction with either the two-stream or the single-scattering model providing a performance enhancement of almost two orders of magnitude while keeping the maximum error below 0.1% (note that the use of the single-scattering model instead of the two-stream model within the framework of the PCA–based method drastically reduces the accuracy).

As future goals, following [18,49] we plan to extend the CLSR method for computing the Stokes parameters and analyze the possibility of hybrid use of the CLSR method, the PCA-based method and the correlated k-distribution technique. It is also planned to analyze different atmospheric scenarios such as cirrus clouds [50] as well as the efficiency of the CLRS method for improving the accuracy of other types of approximate models. In this case, one can expect reducing the number of regression points if instead of the single scattering model, the double scattering approximation is considered [51]. For optically thick media, as an approximate model, asymptotic radiative transfer models can be used [52].