A COMPARISON OF MACHINE / STATISTICAL LEARNING STRATEGIES
FOR MODELING THE RADIATIVE PROPERTIES OF NON-UNIFORM
GASEOUS ATMOSPHERES

X. Aihemaiti¹, T. Ren¹, L. Guilmard², F. André³*

¹China-UK Low Carbon College, Shanghai Jiao Tong University, Shanghai, China, 201306

²Univ. Lyon, CNRS, INSA Lyon, UMR 5008 – CETHIL, Lyon, France

³Univ. Lille, CNRS, UMR 8518 – LOA, Lille, France

ABSTRACT. Fast though accurate methods of gaseous absorption are key components of the forward models used in atmospheric remote sensing studies. In the present work, two techniques founded on Statistical / Machine Learning strategies are considered for the modeling of transmissivities of non-uniform atmospheric paths. One of them is purely data-driven, whereas the second one uses a double mathematical induction strategy. Both approaches are shown to provide accurate estimates of non-uniform path transmissivities at a small fraction of the CPU cost of high-resolution calculations.

# 1. INTRODUCTION

Gas radiation methods with an accuracy higher than 1.0 % (0.5 % in some application cases) are frequently requested in satellite remote sensing studies. This level of precision is hard to ensure with most approximate gas radiation methods (SNB models [1] for instance). Moreover, CPU cost needs to be as small as possible to allow the analysis of the large amount of data generated by recent space instruments. This excludes the use of LBL calculations in many cases. Recently, methods based on machine / statistical learning (ML/SL) and relying on the Universal Approximation Theorem [2,3] have been applied successfully in many fields (finance, biology, etc). Such methods are considered in the present work to model gaseous absorption of non-uniform atmospheric paths for remote sensing applications.

In general, in order to build a ML/SL approximation, a set of training data is provided and the model parameters are constructed by minimizing the Mean Squared Error (MSE) between the training data and the model (the so-called Data Loss Function, DLF) using high dimensional regression techniques. Purely data-driven models mostly interpolate between the input data. Consequently, they may approximate the set of training data very well but may also provide unphysical predictions: generalization performance of this kind of model is indeed rarely guaranteed. Moreover, purely data-driven models usually suffer from a lack of interpretability as the relationship between the inputs and the outputs of the network is in general non-intuitive and strongly non-linear owing to the use of combinations of squashing functions. To address these limitations, various methods have been developed to incorporate physical constraints into ML/SL approximators.

Incorporating physics in a machine / statistical learning framework can be done by three different ways [4]. The first one involves increasing the number of training data to improve the numerical model. This is the weakest though the most intuitive form of embedding of physics in ML. This

<footnote>* Corresponding Author: frederic.andre@univ-lille.fr.</footnote>
approach is called “observational bias” in Ref. [4]. The second method consists of the addition of physical constraints at the learning stage. These constraints usually take the form of a Physical Loss Function (PLF), based on the governing laws of the physical process studied, which is added to the DLF. Adding the physical loss forces the model to follow, as much as possible, the governing equations of the physical problem considered. A typical example of such strategy is found in so-called Physics Informed Neural Networks [5] in which systems of ordinary or partial differential equations together with their corresponding boundary and initial conditions are included in the training process using a PLF. This second method is called “learning bias” in Ref. [4]. The last approach, called “inductive bias” [4], consists of incorporating physical constraints directly inside the structure of the model. A recent example of mathematical induction through the use of Physical Activation Functions is described in Ref. [6]. Mathematical induction is the most principled of the three approaches, but may not be applied in all possible scenarios.

The present paper aims to compare two techniques based on machine / statistical learning strategies for the approximation of non-uniform atmospheric path transmissivities. One method applies a purely data-driven machine learning strategy, whereas the second one uses a double mathematical induction (the structure and choice of activation functions are both founded on physical arguments). Details about the two methods are provided in the paper. The two models are assessed against reference line-by-line (LBL) calculations in realistic atmospheric configurations and compared to each other in terms of computational cost (CPU time) and accuracy when applied for radiative transfer calculations.

## 2. DOUBLE MATHEMATICAL INDUCTION ML/SL STRATEGY

The double mathematical induction strategy considered here finds its roots in the recently proposed $\ell$-distribution framework [7]. This method is presented briefly. Several new developments, specific to the present work, are also described.

The quasi-scaling assumption is the most important element in the construction of the $\ell$-distribution method. The terminology “quasi-scaling” means that for two thermophysical states of the gas, with index 1 and 2 respectively, the band averaged (exponent $\Delta\nu$) transmission functions in the two states are assumed to be related to each other by the relationship:

$$ \tau_2^{\Delta\nu}(L) = \int_0^1 \tau_1^{\Delta\nu}[u(\xi)L] d\xi \quad (1) $$

where the variable $u(\xi)$ plays the role of a non-constant scaling coefficient between the two spectra. Following this definition, two states for which the spectral variables $\kappa_{v,1}$ and $u_v = \kappa_{v,2}/\kappa_{v,1}$ are statistically independent are rigorously quasi-scaled. In this case: 1/ $u(\xi)$ is the inverse of the distribution function of the spectral variable $u_v = \kappa_{v,2}/\kappa_{v,1}$, and, 2/ the equality set by Eq. (1) is exact. It was shown recently [8] that if gas spectra are quasi-scaled and if the spectral variable $u_v$ takes values inside a small interval $u_v \in [u_{min}, u_{max}]$, $|u_{max} - u_{min}| \ll u_{min}$, then the composite function $\ell_1 \circ \tau_2^{\Delta\nu}$ constructed as the combination (symbol $\circ$) of the inverse of the band averaged transmissivity in state 1 written $\ell_1$ s.t. $\ell_1 \circ \tau_1^{\Delta\nu}(L) = L$ and of the band averaged transmissivity in state 2, $\tau_2^{\Delta\nu}$, is of the Bernstein type [9]. It then admits the following Lévy-Khintchine (LK) representation:

$$ \ell_1 \circ \tau_2^{\Delta\nu}(L) = u_{min} \cdot L + \int_0^1 \frac{1 - \exp[-s(0) \cdot v(\xi) \cdot L]}{s(0)} d\xi \quad (2) $$
where $s(0)$ is the rate of change of $1/k_{p,1} \cdot \partial \ln \tau_1^{\Delta\nu} / \partial L$ in $L=0$ and $v(\xi) = u(\xi) - u_{\min}$.

From Eq. (2), two important properties can be shown. The first one (P1) stems from a rewriting of the integral at the RHS in Eq. (2) using a numerical quadrature that shows that the composite function $\ell_1 \circ \tau_2^{\Delta\nu}$ can be approximated as a weighted sum of Exponential Linear Units, ELUs [10]. This property was already noticed and validated thoroughly in Refs. [8,11]. The second important property (P2) is that the transmissivity of a $n$ layers system can be evaluated using a recurrent scheme:

$$ \left\{ \begin{array}{l} \tau_{1..n}^{\Delta\nu}(L_1, \dots, L_n) = \frac{1}{\Delta\nu} \int_{\Delta\nu} \exp(-\kappa_{v,1} L_1 - \dots - \kappa_{v,n} L_n) d\nu = \tau_1^{\Delta\nu}(L_{1..n}) \\ L_{nn} = L_n; \quad L_{i..n} = L_i + \ell_i \circ \tau_{i+1}^{\Delta\nu}(L_{i+1n}) \end{array} \right. \qquad (3) $$

in which all functions $\ell_i \circ \tau_{i+1}^{\Delta\nu}$ admit a Lévy-Khintchine representation of the form of Eq. (2). The recurrent structure of the model Eq. (3) is intimately connected with the functional form set by Eq. (2). The proof of this statement is rather lengthy in the general case and cannot be given here due to limited space. Astonishingly, the iterative process set by Eq. (3) is equivalent to a method proposed intuitively by Godson in the 50s [1,12], as soon as the conditions of applicability of Eq. (2) are fulfilled.

One of the key “problems” with the formulation provided by Eq. (3) is that state 1 (which is rather arbitrary) plays a particular role in the model. Indeed, this is the only state along the non-uniform path for which the length of the sub-path appears linearly in the formulation of the non-uniform equivalent gas pathlength $L_{1..n}$. This particularity can be sidestepped by replacing the first line of Eq. (3) (the second line is not modified) by:

$$ \left\{ \begin{array}{l} \tau_{1..n}^{\Delta\nu}(L_1, \dots, L_n) = \int_0^1 \tau_1^{\Delta\nu}[u_1^+(\xi) \cdot L_{1..n}] d\xi = \tau_1^{\Delta\nu}[\ell_1 \circ \tau_1^+(L_{1..n})] \\ \tau_1^+(L) = \int_0^1 \tau_1^{\Delta\nu}[u_1^+(\xi) \cdot L] d\xi \end{array} \right. \qquad (4) $$

where $u_1^+$ is a bounded random variable with mean value 1 statistically independent of the absorption coefficient in state 1.

In Eq. (4), the composite function $\ell_1 \circ \tau_1^+(L_{1..n})$ is of the Bernstein type. $u_1^+$ plays a role similar to a scaling coefficient, but has no physical meaning in contrast to all the other scaling coefficients involved in the model: $u_1^+$ mostly plays the role of a correction factor. The main effect of the variable $u_1^+$ is to modulate the spectrum in the first layer in such a way that $L_1$ appears non-linearly in the network. Notice that Eq. (3) is a particular case of Eq. (4) for which $u_1^+(\xi) = u_{1,\min}^+ = \text{constant} = 1$. Property P2 remains valid if Eq. (4) is used instead of Eq. (3), but one extra layer (related to $u_1^+$) is added in the network. Properties P1 and P2 are useful to recast the model set by Eq. (4) inside a ML/SL framework. Indeed, they allow a double mathematical induction to construct a network model for the equivalent absorption length $\ell_1 \circ \tau_1^+(L_{1..n})$ of Eq. (4): the use of ELUs relies on P1, whereas the recurrent structure is motivated by P2.
### 3. PRACTICAL IMPLEMENTATION

The training of the Augmented *ℓ*-distribution (ALD) model constructed with the help of the set of Eqs. (3,4), described in Section 2, mostly follows the same steps as in Ref. [8] where interested readers will find extensive implementation details. Due to limited space, we only provide here a few elements that are specific to the model used in the present paper:

1. Usually, *ℓ*-distribution models rely on look-up tables of mapping functions [7]. This parametrization strategy was found to be highly CPU efficient, but also memory expensive. Consequently, in the present work, all integrals that appear in Eq. (2) use a four-point Gauss Legendre quadrature and Eq. (2) in discrete form is used directly inside the system of Eqs. (3,4). Only the mapping function in the first layer is needed. It is tabulated at very high resolution here (100,000 values). The total size of the database for one atmospheric profile (49 distinct thermophysical conditions) and one optical filter is about 2 Mbytes (a full *ℓ*-distribution database based on look-up tables for the same case is roughly 100 times bigger).

2. As noticed in the introduction, most ML/SL regression processes are founded on the minimization of the MSE between the training data and the model. In the present work, MSRE (Mean Squared Relative Error) was used in all cases (C1, C2 and C4 in the next section) for which the minimum non-uniform path transmissivities were found to remain above 0.1 over the full range of training data. MSE was restricted to strongly absorbing scenarios (Case C3). In all calculations, the ADAM optimizer combined with explicit Jacobians was used for the regression. ADAM's hyperparameters are those given in Ref. [13]. The learning rate was set to 0.05 in all calculations.

3. The training set contains 1,440 values of non-uniform path transmissivities (one value every 500 meters starting from the top of the atmosphere down to the ground for 6 values of relative air mass factor e.g. 1, 2, 4, 8, 16 and 24). The same selection of training set was applied in Ref. [8].

For the data-driven machine learning (DDML) model, the absorption length is also divided into 49 homogeneous sub-layers, each layer representing the length traveled by radiation through a specific thermophysical state. These lengths serve as the 49 input features for the machine learning model, as depicted in Figure 1.

<figure><img src="image_1.png" /><figcaption>Figure 1. Schematic of a representative multilayer perceptron (MLP) neural network architecture.</figcaption></figure>

To capture large gradient near the surface, the first sub-layer (closest to the ground) is subdivided into intervals of 10 meters. For the remaining 48 sub-layers, each layer uses increments of 100 meters to represent the distance traveled through higher atmospheric regions. Each input sample in the data set consists of these 49 distances as features, with a single output representing the corresponding transmissivity $τ$. The transmissivity along the entire atmospheric path is modeled as in Eq. (5) where *L* denotes the total path length and $κ_v(u)$ is the spectral absorption coefficient at height *u*:
$$ \tau = \frac{1}{\Delta\nu} \cdot \int_{\Delta\nu} \exp\left(-\int_{0}^{L} \kappa_{\nu}(u)du\right) d\nu \equiv \tau_{1..n}^{\Delta\nu}(L_{1}, \dots, L_{n}), \quad L = L_{1} + \dots + L_{n} \quad (5) $$

The machine learning model aims to learn the relationship between the 49 input features (layer heights) and the output transmissivity. The neural network is implemented using PyTorch. Data manipulation is handled with the Pandas library, while scikit-learn provides additional tools for model evaluation and preprocessing. Training is conducted by minimizing the MSE loss function:

$$ \text{MSE} = \frac{1}{N} \sum_{i=1}^{N} (\hat{\tau}_{i} - \tau_{i})^{2} \quad (6) $$

where $\hat{\tau}_{i}$ denotes the predicted transmissivity for the i-th sample (for which a set of layers' heights $L_{1i}, \dots, L_{ni}$ is given as input), $\tau_{i}$ represents the actual LBL transmissivity for the i-th sample evaluated using Eq. (5), and $N$ is the total number of samples. The dataset is divided into train, validation, and test splits to ensure robust evaluation. It comprises a total of 29,860 samples, with 80% allocated for training, 10% for validation, and 10% for testing. The training set is used to train the model, the validation set monitors performance during training to prevent overfitting, and the test set provides an unbiased assessment of the final model. The model is trained by minimizing the MSE loss function, using the ADAM optimizer to update model parameters. During each epoch, the model processes the 49 input features and compares the predicted transmissivity with the actual LBL transmissivity values.

During training, it was observed that implementing early stopping, combined with a low learning rate, weight decay, and dropout techniques to prevent overfitting, contributed significantly to reducing the loss. Conversely, using an excessively large learning rate may result in gradient explosion, causing the loss to diverge and eventually produce Not a Number (NaN) values, highlighting the importance of careful hyperparameter selection for stable and effective training. To optimize model performance, Bayesian optimization is employed for hyperparameter tuning. By iteratively applying Bayesian optimization on the training set, the optimal combination of parameters—such as learning rate, regularization parameter, and network architecture—can be identified.

Throughout training, the validation dataset is used to observe model performance, reducing the risk of overfitting. After training, the model is evaluated using the relative error metric. Testing on an independent dataset confirms the model's robustness and generalization capabilities, demonstrating reliable transmissivity predictions.

## 4. RESULTS

The four test cases considered in this section are all related to the METImage (Meteorological Imager) instrument on board the Meteorological Operational satellite - Second Generation (MetOp-SG). The space mission is developed by the European Space Agency (ESA) and the European Organization for the Exploitation of Meteorological Satellites (EUMETSAT). METImage was developed by the German Space Agency (DLR). Its main objective is to provide information on cloud, cloud cover, land surface properties, sea, ice and land surface temperatures. Additional information on MetOp-SG mission can be found in MetOp-SG (MetOp-Second Generation Program) - eoPortal. Four channels and associated gaseous molecules are considered for the evaluation of the two methods: the A-band of O₂ in Channel 5 (case C1 in the following), H₂O in Channels 7 (case C2) and 9 (case C3), CH₄ in Channel 11 (case C4). Filter response functions were taken from: nwp-saf.eumetsat.int/downloads/rtcoef_rttov13/ir_srf/rtcoef_metopsg_1 metimage srf.html. Influence of the filter response functions was accounted for in the model by simply replacing in Eqs. (3,5) $d\nu$
by $\phi(\nu)$ $d\nu/\Phi$, where $\Phi$ is the integral of the filter response function. The LBL data used to generate both the training and evaluation sets were calculated at high spectral resolution (0.01 cm$^{-1}$), using HITRAN2020 [14] spectroscopic database. Spectral line profiles were assumed to follow a Voight shape. The Mid Latitude Summer (MLS) atmospheric profile was selected.

Results are plotted in Figures 2 to 5. They are presented in two panels: the right panel shows outcomes derived using a fully data-driven machine learning approach, while the left panel displays results from the Augmented $\ell$-distribution. Although the data-driven machine learning method exhibits larger variability in relative error compared to the Augmented $\ell$-distribution, it always meets the established performance criterion, maintaining error levels within 0.5%. In the case of H$_2$O (Case C3), an increase in relative error is observed near the surface, where transmissivity is exceptionally low (the calculation of relative errors thus involves a division by a very small quantity - this amplifies the relative effect of small deviations).

<figure><img src="image_2.png" /><figcaption>Figure 2. O$_2$ in METImage Channel 5 (case C1)</figcaption></figure>

<figure><img src="image_3.png" /><figcaption>Figure 3. H$_2$O in METImage Channel 7 (case C2)</figcaption></figure>

<figure><img src="image_4.png" /><figcaption>Figure 4. H$_2$O in METImage Channel 9 (case C3)</figcaption></figure>
<figure><img src="image_5.png" /><figcaption>Figure 5. CH₄ in METImage Channel 11 (case C4)</figcaption></figure>

Table 1 provides a comparison of the two models (results for a CKD model with 32 gray gases are also provided for completeness) in terms of CPU cost when applied for the calculation of 1,200 non-uniform path transmissivities.

<table><caption>Table 1. CPU costs in application mode (for each case, the CPU cost is for 1,200 non-uniform path transmissivity calculations, each path is discretized into 49 uniform sub-layers – CPU cost of a CKD model with 32 gray gases is provided for comparison – all calculations were made on a single core of a Xeon Silver 4214R at 2.4 GHz).</caption><thead><tr><th>Method</th><th>CKD32</th><th>DDML</th><th>ALD</th></tr></thead><tbody><tr><td>CPU cost (ms)</td><td>1.8</td><td>3.5</td><td>1.5</td></tr></tbody></table>

From Figures 2-5 and table 1, the following comments can be drawn:

1/ In all the considered test cases, the two models allow achieving a very high level of accuracy, below 0.3% over the full test sets. These test cases are representative of real remote sensing applications. In the strongly absorbing scenario (C3), higher orders of magnitude of relative errors are obtained. They are due to the division by small values of LBL transmissivities in the calculation of relative errors. Once plotted in terms of absolute difference, Figure 3 is similar to the other figures, with absolute differences below $7.10^{-4}$ for the Augmented $\ell$-distribution method.

2/ Both methods are highly efficient compared to LBL calculations. For instance, in the O₂ case (C1), the LBL dataset contains 42,974 spectral values over the filter response function which corresponds to about 2.4 seconds for the calculation of 1,200 values of non-uniform path transmissivities on the same computer as used in Table 1. Both the data-driven and Augmented $\ell$-distribution approaches provide the same set of non-uniform path transmissivities three orders of magnitude faster than LBL, without any significant loss of accuracy.

3/ The higher CPU cost of the data-driven model compared to Augmented $\ell$-distribution is due to the number of neurons required by the model (6 hidden layers are needed, with 256 neurons in layers 1, 2, 5 and 6, and 512 in layers 3 and 4). Even if the activation functions used in this model (linear units) are less CPU expensive than the ELUs used in Augmented $\ell$-distribution, this later model only requires 5 of them for each layer which explains its rather small CPU cost.

4/ Eventually, it should be noticed that the choice of the model requires selecting specific training set and hyperparameters for the regression process. But in both cases model parameters can be found that provide a consistent generalization. For the Augmented $\ell$-distribution model, only 5 coefficients (4 for the LK integrals) per layer are sufficient to obtain a very high accuracy.
## 5. CONCLUSION

In this work, a comparison between two machine / statistical learning strategies for the modeling of radiative properties of non-uniform gaseous atmospheres is described. The first method uses a purely data-driven approach whereas the second one applies a double mathematical induction strategy. If the main advantage of the second method is its interpretability, the set of tests cases considered here have shown that the two approaches yield roughly the same accuracy, higher than what is needed in most application cases. None of the two approaches can thus be considered superior to the other in terms of precision, as they both achieve LBL accuracy. Additionally, the CPU costs of the two models are roughly the same in the realistic scenarios considered in this work, much smaller than high resolution LBL calculations. The Augmented *t*-distribution method is however slightly faster (1.5 ms for 1,200 values of transmissivities) than the data-driven multilayer neural network (3.5 ms for 1,200 values of transmissivities). Acceleration techniques based on path recycling can be used to minimize further the CPU cost of the Augmented *t*-distribution method. Similar acceleration techniques may also exist for the data-driven neural network version.

## ACKNOWLEDGEMENT

This work was supported by the French National Research Agency under the grant ANR-20-CE04-0005 ASGARD.

## REFERENCES

[1] S. J. Young, *Band Model Theory of Radiation Transport*, The Aerospace Press, 2013.

[2] K. Hornik, M. Stinchcombe, H. White, "Multilayer Feedforward Networks are Universal Approximators", *Neural Networks*, vol. 2, pp. 359-366, 1989.

[3] A. M. Schäfer, H.-G. Zimmermann, "Recurrent Neural Networks are Universal Approximators", *Int. J. of Neural Systems*, vol. 17, no. 4, pp. 253-263, 2007.

[4] G. E. Karniadakis, I. G. Kevrekidis, L. Lu, P. Perdikaris, S. Wang, L. Yang, "Physics-informed machine learning", *Nature Review Physics*, vol. 3, Issue 6, pp. 422-440, 2021.

[5] M. Raissi, P. Perdikaris, G. E. Karniadakis, "Physics-informed neural networks: a deep learning framework for solving forward and inverse problems involving nonlinear partial differential equations", *J. Comput. Phys.*, Vol. 378, pp. 686-707, 2019.

[6] J. Abbasi, P. O. Andersen, "Physical activation functions (PAFs): an approach for more efficient induction of physics into physics-informed neural networks (PINNs)", *Neurocomputing*, Vol. 608, 2024.

[7] F. Andre, "The *t*-distribution method for modeling non-gray absorption in uniform and non-uniform gaseous media", *J. Quant, Spectrosc. Radiat. Transfer*, vol. 179, pp. 19-32, 2016.

[8] F. Andre, C. Delage, L. Guilmard, M. Galtier, C. Cornet, "Bridging physics and statistical learning methodologies for the accurate modeling of the radiative properties of non-uniform atmospheric paths", *J. Quant, Spectrosc. Radiat. Transfer*, vol. 320, 2024, 108961.

[9] R. L. Schilling, R. Song, Z. Vondraček, *Bernstein Functions: Theory and Applications*, Second Edition, De Gruyter, 2012.

[10] D.-A. Clevert, T. Unterthiner, S. Hochreiter, "Fast and accurate deep network learning by exponential linear units (ELUs)", In: *Proceedings of the International Conference on Learning Representation ICLR'16*, 2016.

[11] F. Andre, C. Cornet, C. Delage, Ph. Dubuisson, M. Galtier, "On the use of recurrent neural networks for fast and accurate non-uniform gas radiation modeling", *J. Quant, Spectrosc. Radiat. Transfer*, vol. 293, 2022, 108371.

[12] W. L. Godson, "The evaluation of infra-red radiative fluxes due to atmospheric water vapour", *Q. Journal Royal Meteorological Society*, vol. 79, pp. 367-379, 1953.

[13] D. P. Kingma, J. L. Ba, "Adam: a method for stochastic optimization", In: *Proceedings of the 3rd International Conference on Learning Representation ICLR'15*, 2015.

[14] I. E. Gordon, L. S. Rothman, R. J. Hargreaves, R. Hashemi, E. V. Karlovets, F. M. Skinner, E. K. Conway, C. Hill, R. V. Kochanov, Y. Tan, P. Wcisło, A. A. Finenko, K. Nelson, P. F. Bernath, M. Birk, V. Boudon, A. Campargue, K. V. Chance, A. Coustenis, B. J. Drouin...S. N. Yurchenko, "The HITRAN2020 molecular spectroscopic database", *J. Quant, Spectrosc. Radiat. Transfer*, vol. 277, 2022, 107949.