# MT-HCCAR: Multi-Task Deep Learning with Hierarchical Classification and
Attention-based Regression for Cloud Property Retrieval

Xingyan Li1, Andrew M. Sayer23, Ian T. Carroll23, Xin Huang4, and Jianwu Wang1
1 Department of Information Systems, University of Maryland, Baltimore County,
Baltimore, MD, USA  2NASA Goddard Space Flight Center, Greenbelt, MD, USA
3Goddard Earth Sciences Technology and Research (GESTAR) II, UMBC, Baltimore,
MD, USA  4Dept of Computer & Information Sciences, Towson University, Towson,
MD, USA

###### Abstract

In the realm of Earth science, effective cloud property retrieval,
encompassing cloud masking, cloud phase classification, and cloud optical
thickness (COT) prediction, remains pivotal. Traditional methodologies
necessitate distinct models for each sensor instrument due to their unique
spectral characteristics. Recent strides in Earth Science research have
embraced machine learning and deep learning techniques to extract features
from satellite datasets’ spectral observations. However, prevailing approaches
lack novel architectures accounting for hierarchical relationships among
retrieval tasks. Moreover, considering the spectral diversity among existing
sensors, the development of models with robust generalization capabilities
over different sensor datasets is imperative. Surprisingly, there is a dearth
of methodologies addressing the selection of an optimal model for diverse
datasets. In response, this paper introduces MT-HCCAR, an end-to-end deep
learning model employing multi-task learning to simultaneously tackle cloud
masking, cloud phase retrieval (classification tasks), and COT prediction (a
regression task). The MT-HCCAR integrates a hierarchical classification
network (HC) and a classification-assisted attention-based regression network
(CAR), enhancing precision and robustness in cloud labeling and COT
prediction. Additionally, a comprehensive model selection method rooted in
K-fold cross-validation, one standard error rule, and two introduced
performance scores is proposed to select the optimal model over three
simulated satellite datasets OCI, VIIRS, and ABI. The experiments comparing
MT-HCCAR with baseline methods, the ablation studies, and the model selection
affirm the superiority and the generalization capabilities of MT-HCCAR.

###### Index Terms:

cloud property retrieval, remote sensing, deep learning, multi-task learning

##  I Introduction

The retrieval of cloud properties, denoting the estimation of diverse
characteristics of clouds through the analysis of data acquired from remote
sensing instruments, plays an essential role in atmospheric and environmental
studies. This significance arises from the inherent role of clouds as integral
components of the Earth’s atmospheric system, wielding substantial influence
over the planet’s energy dynamics, climate regulation, and hydrological
cycle[1]. The modulation of surface temperatures and the facilitation of
precipitation by clouds underscore their paramount importance in sustaining
ecosystems and shaping broader environmental processes.

Satellites serve as invaluable tools for global monitoring of cloud
properties, offering an extensive data record of Earth’s imagery collected by
various instruments since the latter part of the previous century. A prominent
category of these instruments, commonly referred to as “imagers”, passively
collect measurements of the Earth’s surface across multiple wavelength ranges,
encompassing the ultraviolet (UV), visible (VIS), near-infrared (NIR),
shortwave infrared (SWIR), and thermal infrared (TIR). The data obtained,
primarily consisting of observations of reflected solar and/or emitted thermal
radiation, undergo routine processing by algorithms. This processing
transforms these measurements into geophysical quantities of interest, about
both atmospheric and surface characteristics, in a process commonly termed
“retrieval”. Routine retrieval of geophysical quantities is key to enhancing
the comprehension of Earth’s climate [2]. At the end of such cloud retrieval
process, each observed pixel is equipped with estimated values for many cloud
properties including cloud mask (whether a pixel is clear or cloudy), cloud
thermodynamic phase (liquid droplets or ice crystals), and cloud optical
thickness (COT).

One practical motivation for our work is NASA’s Plankton, Aerosol, Cloud,
ocean Ecosystem (PACE) mission [3], which will launch in 2024 and advance
science in its eponymous disciplines. Heritage cloud masking algorithms are
not directly applicable to PACE’s main sensor, the Ocean Color Instrument
(OCI) [4], due to different spectral bands compared to algorithms developed
for existing sensors. OCI has similarities with some previous spaceborne
imagers used for retrieval of cloud (and other) properties such as Moderate
Resolution Imaging Spectroradiometer (MODIS) [5], Visible Infrared Imaging
Radiometer Suite (VIIRS) [6], and Advanced Baseline Imager (ABI) [7] sensor
types. However, two key differences are that 1) OCI has continuous
hyperspectral coverage in the UV-NIR, plus discrete SWIR bands, while the
others are only multi-spectral (up to a dozen or so discrete relevant bands);
and 2) OCI lacks TIR bands. Some of the most commonly used cloud masking tests
for those sensors are therefore not applicable, and adapting a subset of those
tests would miss out on OCI’s unique abilities, so a new approach is
warranted.

Numerous investigations have delved into the application of machine learning
techniques for the extraction of cloud properties from satellite sensor
observations. These inquiries encompass both singular cloud property
retrievals [8, 9, 10, 11], as well as the concurrent retrieval of multiple
cloud properties[12, 13, 14, 15]. Nevertheless, challenges and unresolved
inquiries persist within the current landscape of machine learning-based cloud
retrieval. First, the role of atmospheric domain knowledge, specifically the
utilization of physical relationships between cloud mask and cloud phase, in
conjunction with advanced machine learning techniques, remains unclear in
enhancing the accuracy of cloud retrieval. Second, despite the deployment or
design of numerous satellite sensors catering to analogous cloud retrieval
tasks (expounded upon in Section III), the feasibility of employing a unified
machine learning architecture to facilitate retrieval across different sensors
and the generalization ability of such a machine learning model remain
uncertain. Third, while there are various metrics to measure machine learning
model performance, it remains ambiguous how an overarching ranking system can
be established to select an optimal model, especially about determining
whether the heightened complexity inherent in newly proposed models is
justified in terms of their capabilities.

To address the above challenges, we propose an end-to-end Multi-Task learning
architecture with Hierarchical Classification and Attention-assisted
Regression (MT-HCCAR) for cloud masking, cloud thermodynamic phase (liquid or
ice) prediction, and COT regression. The model consists of a Hierarchical
Classification network for cloud masking and cloud phase prediction, and a
Classification-assisted Attention-based Regression network (MT-HCCAR) for COT
prediction. The classification task is a hierarchical classification (HC) sub-
network to capture the hierarchical relationship between ice-crystal cloud and
liquid cloud. The COT prediction task is a sub-network integrated with coarse
cloud thickness classification (C) using attention mechanism (A) to further
enhance regression (R) task accuracy. Our implementation is open-sourced
at[16]. Our work’s major contributions are as follows.

  * •

We propose an end-to-end multi-task learning model MT-HCCAR for cloud masking,
cloud phase classification, and COT regression. The model utilizes a
hierarchical classification approach for cloud masking and cloud phase
prediction. Meanwhile, to enhance COT regression accuracy, we use a cross-
attention mechanism to learn similarities across the classification task and
the regression task. Our extensive experiment results show performance
improvement in all tasks using our approach than baseline models.

  * •

We conduct experiments on three satellite sensors (OCI, VIIRS, and ABI,
introduced later) to compare different model performances across different
satellite products. The results show the generalizability of our proposed
model and we analyze our findings from different sensors’ performance.

  * •

To help identify the best model for practical usage, we design a novel model
selection method that ranks the performance of different models across
multiple metrics and multiple datasets.

The subsequent sections of this paper are structured as follows: Section II
presents an overview of related work in the domains of cloud property
retrieval and Multi-Task Learning (MTL) methods for remote sensing. Section
III provides an in-depth exploration of the application background, the
methodology employed for data simulation, and a precise articulation of our
defined problem. Our proposed method, MT-HCCAR, is introduced in Section IV.
The assessment metrics and outcomes of our experiments are elucidated in
Section V. Section VI outlines our strategies for optimal model selection and
presents the corresponding results. Lastly, Section VII encapsulates the paper
with concluding remarks.

##  II Related work

Our study focuses on the application of multi-task learning on pixel-wise
spectral features from satellite images. This section discusses the existing
research on machine learning based cloud property retrieval methods and
related work of multi-task learning applied to remote sensing data.

###  II-A Machine learning based cloud property retrieval

In the field of remote sensing, cloud properties encompass a diverse set of
characteristics and parameters associated with clouds, which can be observed
and measured from a distance using satellite or other remote sensing
technologies. These properties include, among others, cloud mask, cloud
(thermodynamic) phase, cloud vertical distribution, and COT. Machine learning
approaches have demonstrated their efficacy in processing high-dimensional
remote sensing data, automatically learning intrinsic features from spectrum
data.

In recent years, various machine learning techniques, such as Random Forest
(RF), Multi-Layer Perceptron (MLP), and Convolutional Neural Network (CNN),
have been employed for retrieving different cloud properties. For instance,
Liu et al. [8] utilized Artificial Neural Networks (ANN) and RF for cloud
detection using spectral observations and reference labels from the Cloud-
Aerosol Lidar with Orthogonal Polarization (CALIOP). Wang et al. [9] proposed
an RF-based algorithm using passive spectral observations for cloud detection
and phase classification across multiple instruments, evaluating their models
against CALIOP and existing VIIRS cloud products. Giuffrida et al. [10]
introduced CloudScout, a custom CNN deployed on the Hyperscout-2 nanosatellite
payload for onboard cloud detection on hyperspectral images. He et al. [11]
proposed a lightweight network, DABNet, for high-accuracy cloud detection
based on CNN. Mohajerani et al.[17] presented Cloud-Net, a deep learning-based
model using convolutional layers for cloud detection in Landsat 8 imagery.
Shao et al. [18] proposed a multi-scale features CNN method for simultaneous
classification of imagery into thick, thin, and non-cloud regions,
outperforming other cloud detection methods. Jeppesen et al. [19] introduced
RS-Net, a deep learning-based cloud detection algorithm for optical satellite
imagery, achieving high performance even with a single spectral band. Segal-
Rozenhaimer et al. [20] proposed a domain adversarial neural network (DANN) to
enhance the transferability of cloud masks. Huang et al.[14] presented DAMA-
WL, an end-to-end deep domain adaptation approach for cloud phase detection,
incorporating domain alignment of heterogeneous active and passive satellite
remote sensing data. Guo et al.[21] introduced CDnetV2, a CNN-based cloud
detection method for remote sensing imagery with cloud-snow coexistence,
utilizing attention-based modules for adaptive feature fusion. López-
Puigdollers et al.[22] conducted a benchmark of deep learning models for cloud
detection in multispectral images, including both rule-based approaches and
machine learning approaches.

While these machine learning approaches for cloud property retrieval
predominantly leverage spectral features, the main limitations lie in two
aspects: 1) many studies do not consider background knowledge, such as the
task order (e.g., cloud mask prediction, cloud phase prediction, and COT
prediction), and 2) several studies conduct different tasks independently,
lacking parameter sharing between classification and regression tasks. These
limitations underscore the need for more integrated and informed methodologies
in cloud property retrieval studies.

###  II-B Multi-task learning for remote sensing data

Multi-task (MT) learning has emerged as a valuable approach in Earth science
and remote sensing applications, facilitating the retrieval of related
products from shared features. Recent studies have demonstrated the efficacy
of the multi-task learning paradigm in jointly enhancing the performance of
diverse tasks.

Ilteralp et al. [23] proposed an MT Convolutional Neural Network (CNN)
architecture capable of leveraging unlabeled satellite imagery to predict
Chl-a values, incorporating an auxiliary classification task. The proposed MT
learning model was compared against alternative approaches, including its
single-task regression version, support vector regression, Multi-Layer
Perceptron (MLP), morphology-based spectral-spatial approach (EFAL), CNN, and
statistical spectral band characteristics known as the band ratio. Zhang et
al. [24] introduced MTDR-Net, an MT-driven and reconfigurable network designed
to improve cloud detection in regions with coexisting cloud and snow cover in
very high-resolution remote sensing images. This model employs an MT gradient
flow guidance module (MTGFG) to establish a feature-sharing strategy across
different tasks. Chen et al. [25] proposed an MT learning method based on Long
Short-Term Memory (LSTM) with an attention mechanism to classify imbalances in
wind turbine blades caused by ice accretion. Sarafian et al. [26] presented a
meteorology-based model for dust event forecasting, utilizing an MT framework
for predicting regional PM10 fields and local in situ measured PM10 values.
Kim et al. [27] proposed MT CNN and ConvLSTM models for simultaneous
forecasting of pixel-wise sea ice concentration, considering both spatial sea
ice concentration (SIC) images and temporal sea ice extent (SIE) values. Qiu
et al. [28] developed an end-to-end CNN with a shared backbone network for
shared feature learning, attention modules for task-specific feature learning
for human settlement extent (HSE) regression, and local climate zone (LCZ)
classification. Hanna et al. [29] proposed a U-Net-based architecture to
simultaneously predict the pixel area covered by plumes, the type of fired
fuel, and the power generation rate from a single-satellite image of a power
plant. Ahn et al. [30] presented a systematic approach to forecasting regional
wave heights based on different LSTM structures. Zhao et al. [31] introduced a
CNN-based framework that combines self-supervised learning and scene
classification to improve feature extraction and generalization abilities.

There are also studies conducting ML learning specifically for cloud property
retrieval. Yang et al.[12] developed an MLP-based method to retrieve cloud
macrophysical parameters (cloud mask, cloud top temperature, and cloud top
height) using Himawari-8 satellite data. Wang et al.[13] proposed TIR-CNN
based on the U-Net model to retrieve cloud properties including cloud mask,
COT, effective particle radius (CER), and cloud top height (CTH) from thermal
infrared radiometry. The architecture consists of encoding and decoding
layers, convolutional blocks, batch normalization layers, and leaky rectified
linear units. The results of applying the model to thermal infrared spectral
data from MODIS are used to compare model performance for daytime and
nighttime data.

While these above studies offer valuable insights that are beneficial to our
design, our work cannot be compared with most of them directly via
experiments. It is because of the distinctive nature of our model’s output,
which manifests at the pixel level rather than the image level. Consequently,
architectures employing convolutional layers, such as U-net architecture, are
ill-suited for our study’s requirements.

##  III Problem statement and data simulation

In this section, we first present the contextual background of the application
which includes a forward simulation process and an inverse retrieval process.
Then, we delve into a comprehensive exposition of the employed data simulation
method in Section III-A, and define our cloud property retrieval problem based
on the generated simulation data in Section III-B.

The forward simulation and inverse retrieval processes are characterized as
two opposite operations within the context of cloud property analysis. In the
forward simulation, spectral bands and physical variables are generated as
outputs using cloud properties as input parameters. Conversely, cloud property
retrieval involves utilizing spectral features obtained from satellite
observations as input, with the resultant output comprising cloud properties
such as cloud mask, cloud phase, and COT values. To provide a more explicit
formulation, the forward simulation is represented as the mathematical
process:

$\mathbf{X}=f(\mathbf{L}_{cm},\mathbf{L}_{cph},\mathbf{Y}_{COT},...)$ (1)

, where $\mathbf{X}=[\mathbf{x_{1}},\mathbf{x_{2}},\ldots,\mathbf{x_{M}}]$
denotes the spectral and physical features, $f(\cdot)$ represents the
simulation algorithm, and $\mathbf{L}_{cm},\mathbf{L}_{cph},\mathbf{Y}_{COT}$
are the labels corresponding to cloud mask, cloud phase, and COT values,
respectively. Additional properties are also used as inputs (see Section
III-A), but they are not the focus of this work.

Conversely, the inverse retrieval process can be expressed as:

$\mathbf{\hat{L}}_{cm},\mathbf{\hat{L}}_{cph},\mathbf{\hat{Y}}_{COT}=g(\mathbf{X})$

, where
$[\mathbf{\hat{L}}_{cm},\mathbf{\hat{L}}_{cph},\mathbf{\hat{Y}}_{COT}]$
represents the predicted values of cloud mask, cloud phase, and COT, as
determined by the retrieval algorithm $g(\cdot)$ given the input spectral and
physical features $\mathbf{X}$. Figure 1 illustrates how these two processes
operate in our study and the colors of the three sensors show their wavelength
band differences.

Figure 1: Illustration of forward simulation and inverse retrieval. Forward
simulation uses cloud properties as input parameters to generate spectral
bands and physical variables as outputs. Conversely, inverse retrieval entails
using spectral features obtained from satellite observations as inputs to
generate cloud properties such as cloud mask, cloud phase, and COT values.

###  III-A Forward radiative transfer simulation

Given that OCI has not been launched, simulated datasets for OCI, VIIRS, and
ABI are employed for assessing sensors and architectures against a
standardized source. Radiative transfer (RT) models play a crucial role by
capturing optical processes like scattering and absorption. These models map
geophysical parameters, representing surface and atmospheric properties, to
spectral intensities at specific detection points—essentially serving as the
intermediary link between geophysical quantities and observations made by
space-borne sensors.

We cover from ultraviolet (UV) to shortwave infrared (SWIR) and then convolve
these hyperspectral simulations with the solar spectrum from [32] and sensor
relative spectral response functions for OCI (nominal pre-launch) [4], ABI
[7], and VIIRS [6] in order to generate the spectral top of atmosphere
reflectance signal that the instruments would observe. This provides our
simulated observations and reference truth (i.e., cloud classification and
COT), with band centers simulated for each sensor as below:

  * •

OCI: 233 bands in total; 226 from the hyperspectral detectors from 350-890 nm
(with 2.5 nm spacing between bands) plus seven discrete NIR/SWIR bands
centered near 940, 1040, 1250, 1378, 1620, 2130, and 2260 nm.

  * •

VIIRS: 10 bands in total, centered near 412, 445, 488, 555, 672, 865, 1240,
1380, 1610, and 2250 nm.

  * •

ABI: 6 bands in total, centered near 471, 640, 860, 1370, 1600, and 2200 nm.

Each ABI band has a close analog in VIIRS, and likewise, each VIIRS band has a
close analog in OCI, so the trio represent a set of increasing complexity.
Both ABI and VIIRS also have thermal infrared (TIR) bands, but these are
excluded as the primary motivation for this work is the development of a new
methodology for OCI (which does not), and the application to ABI and VIIRS
demonstrates the broader applicability of our new architectures.

Inputs for the RT simulation include cloud (phase, microphysical and optical
properties, and vertical structure), atmospheric profile (aerosol, gas,
Rayleigh scattering, temperature, and pressure), surface (spectral
reflectance), and observation (solar and view angles) properties. An ensemble
of 250,000 simulations, parameterized by independent samples from realistic
distributions of these geophysical properties, form the labeled dataset for
model training and evaluation. Distributions and further specifics of the
simulations are described in Sayer et al. [33], with two deviations: 1) that
study consisted of purely cloud-covered parameterizations; here, $1/3$ of
simulations corresponded to cloud-free conditions, and the rest to single-
layer clouds, and 2) that study considered a narrow spectral range so assumed
a spectrally-flat surface albedo; here, we use 20-class extended International
Geosphere-Biosphere Program (IGBP) albedo system [34]. These 20 classes
correspond to 16 types of land surface spectra (1 desert, 15 non-desert), 1
water surface, and 3 snow/sea ice surfaces. We sampled such that a given
simulation had a 1/34 chance of being drawn from each of the 16 “land”
spectra, a 9/34 chance of using the “water” spectrum, and a 3/34 of being
drawn from each of the “snow/ice” spectra, and for network training classes
collapse the 20 spectra into these four main categories (water; snow/ice;
desert; non-desert land). The simulation does not include multi-layer or
mixed-phase cloud systems, as most satellite retrieval processing algorithms
assume a single-layer single-phase cloud system. However, future work could
extend the network to be trained (or just evaluated) on an ensemble including
such systems.

We emphasize that the software simulates individual point observations (i.e.,
a pixel) which means that (as is standard for most cloud retrievals in Earth
science) the network operates on single pixels at once as opposed to larger
scenes (i.e., no spatial context). It would be computationally prohibitive to
simulate a globally representative ensemble of satellite scenes, particularly
the inclusion of 3D radiative transfer effects. As a result, existing works
involving 3D effects[35, 36] have typically been limited to a few case
studies.

###  III-B Problem definition for inverse cloud property retrieval

Our primary objective is to accurately model the base-10 logarithm of COT
(throughout the remaining text, COT means this logarithm unless explicitly
described as the “original” COT) for pixels labeled as cloudy. We work in log
space for COT as this has a more linear relationship with the brightness seen
by satellite than the original COT. Simultaneously, the model should be able
to accurately classify the cloud phase [cloudy, cloud-free, cloudy-liquid,
cloudy-ice] for each pixel in the dataset to aid COT prediction. Therefore,
our problem consists of two tasks: 1) a classification task to predict cloud
mask and phase, and 2) a regression task to predict COT values. Assuming that
the labels are $\bar{C}=$ ‘cloud-free’, $C=$ ‘cloudy’, $CL=$ ‘cloudy-liquid’,
and $CI=$‘cloudy-ice’, the details of our study are outlined below.

  * •

Model input. Our model takes 2-dimensional (2D) data as input that contains a
set of sample pixels $S={s_{1},s_{2},\ldots,s_{N}}$, where $N$ is the number
of labeled pixels. Each pixel has the following features and target variables.

    * –

Explanatory features: The input features are represented as
$\mathbf{X}=[\mathbf{x_{1}},\mathbf{x_{2}},\ldots,\mathbf{x_{M}}]$, where $M$
is the dimension of the available features. The exact input feature variables
include: 1) surface pressure in millibar (mbar), 2) total column water vapor
content in millimeters (mm), 3) total column ozone content in the Dobson unit
(DU), 4) types of Earth surface with 4 categories (land, snow, desert, and
ocean or water), 5) top of atmosphere reflectance at different wavelengths
collected by each spaceborne sensor, 6) viewing zenith angle, solar zenith
angle, and the relative azimuth angle, in degree.

    * –

Cloud mask/phase label (denoted as $\mathbf{l^{cls}}$): Cloud mask/phase
labels are used to train the classification task of the model. The set of
possible label values is $l=\\{C,\bar{C},CL,CI\\}$.

    * –

True COT value (denoted as $\mathbf{y^{cot}}$): the values are used to train
the regression task of the model.

  * •

Model output

    * –

Predicted cloud mask/phase class $\mathbf{\hat{l}^{cls}}$ with probabilities
$\mathbf{u}$ of each pixel belonging to each of the four classes.

    * –

Predicted COT value $\mathbf{\hat{y}^{cot}}$ for each pixel.

    * –

Model parameters $\mathbf{\beta}$.

Specifically, there are three pieces of background knowledge related to the
targeted variables $\mathbf{y^{cot}}$ and $\mathbf{l^{cls}}$ in the dataset
for model training.

Knowledge 1: Original COT values are not available for pixels with no cloud
cover, so not every pixel has COT value. That is, if ${l^{cls}_{i}=0}$, then
$y^{cot}_{i}=N/A$.

Knowledge 2: There is a hierarchical relationship between the labels, as
pixels covered by liquid cloud or ice cloud are both cloudy pixels. Thus, $CL$
and $CI$ can be two subclasses of label $C=$‘cloudy’.

Knowledge 3: Each pixel is either a cloudy pixel or a cloud-free pixel, and
the cloudy pixel is either classified as a cloudy-liquid pixel or a cloudy-ice
pixel. There is no coexistence between liquid cloud and ice cloud in our data.
That is, each pixel only be labeled as $\bar{C}$, or as $C$ together with one
of $\\{CL,CI\\}$.

With the input features $\mathbf{X}$ and true values of the target variables
$[\mathbf{y^{cot}},\mathbf{l^{cls}}]$, we aim to learn a model
$\mathcal{M}(\beta)$ such that
$[\mathbf{\hat{y}^{cot}},\mathbf{\hat{l}}^{cls}]=\mathcal{M}(\mathbf{X|\beta})$.

##  IV MT-HCCAR Model

Figure 2: MT-HCCAR: an end-to-end multi-task learning model with hierarchical
classification (HC) and cross attention assisted regression (CAR). The HC sub-
network consists of the cloud masking module $C_{Mask}(\cdot)$ and the cloud
phase classification module $C_{Phase}(\cdot)$. The CAR sub-network consists
of the auxiliary coarse classification module $C_{Aux}(\cdot)$, the cross
attention module $A(\cdot)$, and the regression module $R(\cdot)$. On the
right is the structure of $A(\cdot)$.

Our research objective involves the training of deep learning models to
accomplish two tasks: 1) the classification of cloud mask and phase for each
pixel based on its reflectance values, and 2) the subsequent prediction of
Cloud Optical Thickness (COT) values for pixels classified as cloudy. To
address these objectives, we propose an end-to-end Multi-Task Learning (MTL)
model, denoted as MT-HCCAR, which integrates a Hierarchical Classification
network (HC) and a Classification-assisted with Attention-based Regression
network (CAR). Illustrated in Figure 2, the model commences with an encoder,
and three branches including Decoder, HC, and CAR follow. Subsequent sections
will expound on the individual modules within MT-HCCAR.

###  IV-A Encoder-decoder sub-network

An encoder-decoder sub-network, containing an encoder module and a decoder
module, is employed in our model to learn latent feature parameters that can
be shared for other learning tasks. To integrate the COT regression task and
cloud mask/phase prediction task, we adopt the soft-parameter sharing approach
of MTL where the shared parameters are derived by encoder-decoder. Previous
studies [37, 38] have used similar encoder-decoder techniques to transform
original features into more relevant features for better downstream prediction
accuracy. In our model, the decoder reconstructs the input feature
$\mathbf{X}$ into $\mathbf{\hat{X}}$, and a loss is between $\mathbf{X}$ and
$\mathbf{\hat{X}}$ is minimized during training to minimize the distortion of
input feature throughout feature extraction layers.

The shared encoder is formed by the first three dense layers, each of which is
followed by a ReLU activation function. This feature extractor is a wide-to-
narrow structure where the number of kernels at each layer decreases
gradually. The decoder and reconstruction branch form a regularized
autoencoder which encourages the model to learn features in the aid of the
classification task and regression task.

###  IV-B Hierarchical classification (HC) sub-network

Utilizing the shared parameters acquired by the encoder-decoder sub-network,
MT-HCCAR incorporates a hierarchical classification sub-network. This sub-
network comprises a cloud mask classification module and a cloud phase
classification module, facilitating cloud mask and phase predictions. To
enhance the physical interpretability of the model, the architecture is
crafted to mirror human cognitive processes in understanding cloud labels.
Fundamentally, to predict cloud phase and Cloud Optical Thickness (COT) values
from an Earth science perspective, the model must discern between cloudy and
non-cloudy pixels before further classifying liquid-phase and ice-phase pixels
and predicting COT values for cloudy pixels.

Hierarchical classification has the capability to categorize instances
according to label levels, forming a tree-like structure where each label
functions as a node on the tree[39]. The HC network used in MT-HCCAR consists
of two classifiers $C_{Mask}(\cdot)$ and ${C_{Phase}}(\cdot)$, the first of
which distinguishes cloudy pixels from the cloud-free pixels, and the output
is uncertainties for the two labels
$[\mathbf{\hat{u}^{C}},\mathbf{\hat{u}^{\bar{C}}}]$. The output uncertainties
are binarized to $[\mathbf{\hat{l}^{C}},\mathbf{\hat{l}^{\bar{C}}}]$ by a
threshold of 0.5. Before the beginning of the second classifier
${C_{Phase}}(\cdot)$ input feature $\mathbf{X}$ is multiplied by the predicted
cloud label, $\mathbf{\hat{y}^{C}}$, to filter out cloud-free pixels. Then
${C_{Phase}}(\cdot)$ classifies liquid-cloud and ice-cloud and the sigmoid
uncertainties $[\mathbf{\hat{u}^{CL}},\mathbf{\hat{u}^{CI}}]$ are transformed
in the same way to binary cloud-phase labels
$[\mathbf{\hat{l}^{CL}},\mathbf{\hat{l}^{CI}}]$.

###  IV-C Classification assisted regression sub-network based on cross
attention mechanism (CAR)

Taking inputs from the encoded feature set and cloud mask classification
results as input, our regression sub-network is used for further downstream
regression prediction of a cloudy pixel’s COT value. Two novel efforts were
taken to improve COT prediction accuracy. First, instead of a direct
regression module, an auxiliary coarse classification module is added to
predict which sub-range of COT each pixel falls into. Second, a residual-based
cross-attention mechanism inspired by [40] is developed to enable the
regression module and the auxiliary classification module to share relevant
correlations and insights.

As illustrated in Figure 2, there is a connection between internal features
from the auxiliary classifier and those from the regression network,
facilitated by a residual-based cross-attention module. The close alignment in
tasks between the auxiliary classifier and the regression network enables the
cross-attention mechanism to selectively attend to features relevant to one
task while simultaneously executing the other. This contrasts with
relationships involving the regression task and other tasks, such as cloud
phase classification, where associations are forged through joint optimization
of their respective losses during training iterations, rather than through the
shared utilization of internal features.

Our auxiliary coarse classification involves discretizing continuous values
into coarse groups, serving as a preprocessing step before regression to align
pixels with similar characteristics in the feature space. In MT-HCCAR, an
auxiliary classifier is employed to categorize continuous COT values into
three distinct levels: thin cloud, moderate cloud, and thick cloud.
Specifically, we define thin cloud pixels as those with logarithmic COT values
within the range of [-1.5, 0], moderate cloud in [0, 1], and thick cloud in
[1, 2.5]. Given the analogous nature of the tasks in COT coarse classification
and COT regression, the features extracted by the auxiliary classifier are
integrated with those obtained from the regression network.

To learn the joint features more effectively, we utilize cross-attention
layers to enhance the integration of features from both regression layers and
coarse classification layers, facilitating deeper feature learning and more
cohesive feature amalgamation. Based on the cross-attention mechanism, adding
a residual connection guarantees that the module has stable output and related
research such as Wang et al. [41] propose a non-local operation network using
a residual connection to insert blocks to the network. By using the residual
connection in our module $A(\cdot)$, the output is

$A(\vartheta_{1},\vartheta_{2})=W_{z}y_{A}(\vartheta_{1},\vartheta_{2})+\vartheta_{1}$
(2)

, where as shown in Figure 2 $\vartheta_{1}$ is the feature from the
regression network, $\vartheta_{2}$ is the feature from the auxiliary
classification module, $W_{z}$ is the weight matrix to map $y_{A}$ to
$\vartheta_{2}$, and $y_{A}$ denotes the computation within cross-attention
mechanism. The calculation of $y_{A}$ involves three variables: a query ($Q$),
a key ($K$) and a value ($V$)[40]. They are calculated based on corresponding
weights matrices $W_{Q}$, $W_{K}$, $W_{V}$ learned through model training.
$Q=W_{Q}\vartheta_{2}$, $K=W_{K}\vartheta_{1}$, and $V=W_{V}\vartheta_{1}$.
Then the $y_{A}$ in Equation 2 is:

$y_{A}=Softmax(QK^{T})V$ (3)

###  IV-D Model training of MT-HCCAR

The loss function $\mathcal{L}_{MT-HCCAR}$ for training MT-HCCAR model is
formulated as the weighted sum of four components: a hierarchical
classification loss $\mathcal{L}_{HC}$, a regression loss $\mathcal{L}_{CAR}$,
a reconstruction loss $\mathcal{L}_{Rec}$, and a Lasso regularization loss
$\mathcal{L}_{Lasso}$. That is

$\mathcal{L}=\mathcal{L}_{HC}+\mathcal{L}_{CAR}+\mathcal{L}_{Rec}+\mathcal{L}_{Lasso}$
(4)

The four loss components are calculated using different rules. The
reconstruction loss $\mathcal{L}_{Rec}$ describes the difference between the
input features $\mathbf{X}$ and the reconstruction of the input generated
through the encoder and the decoder. It is Mean Square Error (MSE) between
$\mathbf{D}(\mathbf{X})=\mathbf{\hat{X}}$ and $\mathbf{X}$:

$\mathcal{L}_{Rec}=\sum_{i=1}^{n}(\mathbf{X_{i}}-\mathbf{\hat{X_{i}}})^{2}$
(5)

The Lasso regularization loss component $loss_{l1}$ is an additional penalty
to regularize the training process by increasing the sparsity of the model
with $l1$ regularization loss:

$\mathcal{L}_{Lasso}=\lambda\sum_{p=1}^{P}|\beta_{p}|$ (6)

The classification loss and regression loss are a weighted sum of losses from
the supplement modules. The HC loss is the summation of Binary Cross Entropy
loss from cloud masking module $C_{Mask}(\cdot)$ and cloud phase
classification module $C_{Phase}(\cdot)$.

$\mathcal{L}_{HC}=\mathcal{L}_{C_{Mask}}+\mathcal{L}_{C_{Phase}}$ (7)

We have a sigmoid output of cloudy and non-cloud uncertainties
$[\hat{u}^{C},\hat{u}^{\bar{C}}]$ from cloud mask classifier $C_{Mask}(\cdot)$
and sigmoid output of liquid cloud and ice cloud uncertainties
$[\hat{u}^{CL},\hat{u}^{CI}]$ from cloud phase classifier $C_{Phase}(\cdot)$,
then for Equation 7$,\mathcal{L}_{C_{Mask}}$ and $\mathcal{L}_{C_{Phase}}$ are
calculated as:

$\mathcal{L}_{C_{Mask}}=-\frac{1}{N}\sum_{i=1}^{N}l^{C}_{i}\cdot\log(\hat{u}^{C}_{i})+l^{\bar{C}}_{i}\cdot\log(\hat{u}^{\bar{C}}_{i})$
(8)
$\mathcal{L}_{C_{Phase}}=-\frac{1}{N}\sum_{i=1}^{N}\hat{u}^{C}_{i}l^{CL}_{i}\cdot\log(\hat{u}^{C}_{i}\hat{u}^{CL}_{i})+\hat{u}^{C}_{i}l^{CI}_{i}\cdot\log(\hat{u}^{C}_{i}\hat{u}^{CI}_{i})$
(9)

The CAR loss $\mathcal{L}_{CAR}$ is the summation of Cross Entropy loss from
the auxiliary classifier $C_{Aux}(\cdot)$ and $l1$ loss from regression.

$\mathcal{L}_{CAR}=\mathcal{L}_{R}+\mathcal{L}_{C_{Aux}}$

We have predicted COT values and $\hat{y}^{cot}$ and sigmoid output of
thickness group uncertainties $[\hat{u}^{thin},\hat{u}^{mod},\hat{u}^{thick}]$
from auxiliary classifier, then the regression loss can be formulated as:

$\mathcal{L}_{R}=\sum_{i=1}^{N}|{y}^{cot}_{i}-\hat{y}^{cot}_{i}|$ (10)
$\mathcal{L}_{C_{Aux}}=-\sum_{i=1}^{N}\sum_{c\in[thin,mod,thick]}l^{c}_{i}\ln{\hat{u}^{c}_{i}}$
(11)

##  V Experiments

This section presents the comparison with baseline methods and ablation
studies, where models are trained by the OCI simulated dataset. Section V-A
presents the experimental setup for evaluations. Subsequently, Section V-B
introduces the evaluation metrics. Following this, Section V-C and Section V-D
undertake a comparative analysis between our model and two baseline methods,
evaluating performance in both classification and regression tasks. Ablation
experiments are then conducted to discern the significance of the modules
employed in our model from an earth science perspective.

###  V-A Experiment setting

The dataset comprises satellite data simulation of instruments OCI, VIIRS and
ABI. The three datasets encompassing $N=250,000$ instances. We split data into
62.5%, 22.5%, and 10.0% for training, validation, and test sets, respectively.

To evaluate our model’s performance, we compare MT-HCCAR with six other models
on the same dataset OCI. The first two are baseline methods derived from prior
studies. The baseline methods are selected from related work with similar task
objectives and data formats. Because our dataset contains independent pixels,
the baseline methods that can be adapted to our experiments are:

  * •

Chen et al. (2018): A Multi-layer Perceptron (MLP) based method proposed by
Chen et al.[42] for cloud property retrieval. We use 1 hidden layer with 10
nodes, which are the same as the authors used.

  * •

Liu et al. (2022): A random forest (RF) based method implemented by Liu et
al.[8] for cloud detection. We apply one RF to cloud masking and the same RF
to cloud phase classification, respectively. The parameters involved are
$ntrees=100$ and $mdepth=10$, which are the best parameters for RF found by
the authors.

Besides the above two baseline methods, four ablation study models are also
used to evaluate the effectiveness of different modules in our model. In these
four models, specific modules are systematically removed from our proposed
model, MT-HCCAR, allowing us to assess the impact of these modifications on
model performance.

  * •

SEQ: A sequential structure model consists of three separate networks. The
first one classifies [cloudy, cloud-free], the second one classifies [liquid-
cloud, ice-cloud], and the third network makes regressive predictions on COT
for cloudy pixels, respectively.

  * •

MT-CR: MTL method combines the classification network and regression network
in an end-to-end model.

  * •

MT-HCR: Multi-task learning model with the HC network classification task and
a regression network in the multi-task learning approach.

  * •

MT-HCCR: Multi-task learning model with hierarchical classification and
classification-assisted regression. Only the cross-attention module is removed
from MT-HCCAR.

After comparing different models on the same dataset OCI, we further use
10-fold cross-validation [43, 44] on four MTL-based models MT-CR, MT-HCR, MT-
HCCR, and MT-HCCAR to assess the usefulness of modules and compare the
difference in model performance between our models and two existing methods
for cloud property retrieval. Additionally, each of the four models will be
trained by data from three sensors (OCI, VIIRS, and ABI) respectively.

###  V-B Evaluation metrics

We use three types of metrics to evaluate our work: 1) metrics for the
classification task, 2) metrics for the regression task, 3) metrics from an
Earth science perspective.

####  V-B1 Metrics for classification performance

We use Accuracy and Average precision to evaluate classification performance.

Cloud masking accuracy (ACCbi): Accuracy is the fraction of correct
predictions by our model. We have binary accuracies (ACCbi) between the two
big categories [cloudy, cloud-free] and accuracy for each class of [cloud-
free, liquid-phase, ice-phase].

$Accuracy_{i}=(TP_{i}+TN_{i})/(TP_{i}+TN_{i}+FP_{i}+FN_{i})$ (12)

, where $TP_{i}$, $TN_{i}$, $FP_{i}$, and $FN_{i}$ are the number of true
positives, true negatives, false positives, and false negatives for class $i$,
respectively.

Area under precision-recall curve ($AU(\overline{PRC})_{w}$,
$AU(\overline{PRC})_{c}$): Weighted area under the precision-recall curve
($AU(\overline{PRC}_{w}$)) is the weighted mean of precisions at each
threshold $h$ for class $i$, with the weight as substract between recall at
threshold $h$ and the recall at threshold $h-1$. The area under the precision-
recall curve for each label $AU(\overline{PRC})_{c}$ is calculated for each
binary label in $l=\\{C,\bar{C},CL,CI\\}$. For the weighted area under the
precision-recall curve, the recall and the precision at threshold $h$ are:

$R^{h}_{w}=\frac{\sum_{l}TP^{h}_{c}}{\sum_{l}TP^{h}_{c}+\sum_{l}FN^{h}_{c}},P^{h}_{w}=\frac{\sum_{l}TP^{h}_{c}}{\sum_{l}TP^{h}_{c}+\sum_{l}FP^{h}_{c}}$
(13)

And for the area under the precision-recall curve of each label $c$, the
recall and the precision at threshold $h$ are:

$R^{h}_{c}=\frac{TP^{h}_{c}}{TP^{h}_{c}+FN^{h}_{c}},P^{h}_{c}=\frac{TP^{h}_{c}}{TP^{h}_{c}+FP^{h}_{c}}$
(14)

With Equation 13 and Equation 14, $AU(\overline{PRC})_{w}$ and
$AU(\overline{PRC})_{c}$ are:

$AU(\overline{PRC})_{w}=\sum_{h}(R^{h}_{w}-R^{h-1}_{w})P^{h}_{w}$ (15)
$AU(\overline{PRC})_{c}=\sum_{h}(R^{h}_{c}-R^{h-1}_{c})P^{h}_{c}$ (16)

####  V-B2 Metrics for regression performance

We use MSE and $R^{2}$ to evaluate regression performance.

Mean squared error (MSE) is the average squared difference between the true
and predicted COT values. We calculate MSE across all classes and MSE for each
class of [cloud-free, liquid-phase, ice-phase].

$MSE=\sum_{i=0}^{n}(y_{i}^{cot}-\hat{y}_{i}^{cot})$ (17)

Coefficient of determination ($R^{2}$) measures how close the predicted value
is to the true value. $R^{2}$ across all classes and $R^{2}$ for each class.

$R^{2}=1-\frac{\sum_{i=0}^{n}(y_{i}^{cot}-\hat{y}_{i}^{cot})^{2}}{\sum_{i=0}^{n}(y_{i}^{cot}-\overline{y}^{cot})^{2}}$
(18)

, where $\overline{y}^{cot}$ is the mean value of all true COT values
$[y_{1}^{reg},y_{2}^{reg},\ldots,y_{n}^{reg}]$.

####  V-B3 Earth science metrics

Besides the above commonly used data mining metrics, we also evaluate how a
machine learning model could be used for actual satellite missions.

Fraction of pixels meeting PACE goals (FMG) is an evaluation metric defined by
the PACE Validation Plan [45] based on scientific requirements and
expectations. For each pixel, the relative error
$Error_{i}=|\frac{y^{(cot)_{i}}-\hat{y}^{(cot)_{i}}}{y^{(cot)_{i}}}|$. Then
fraction meeting goals represents the percentage of pixels whose relative
error is less than 0.25 (for liquid clouds) or 0.35 (for ice clouds).

The Validation Plan defines a predictive model with satisfactory performance
as one where this goal is met for 65% of cloudy pixels. Note that this mission
goal only applies to cases where the true original COT value is 5 or more
(log10 COT > 0.7).

TABLE I: Model performance Comparison with Baseline Models and Ablation Study - classification results. Models are trained by OCI. Model | ACCbi | $AU(\overline{PRC})_{w}$ | $AU(\overline{PRC})_{c}$
---|---|---|---
|  |  | Cloudy | Cloud-free | Liquid cloud | Ice cloud
Liu et al.,2022[8] | 0.968 | 0.985 | 0.955 | 0.949 | 0.951 | 0.965
Chen et al., 2018[42] | 0.942 | 0.638 | 0.974 | 0.882 | 0.385 | 0.348
SEQ | 0.963 | 0.750 | 0.997 | 0.978 | 0.618 | 0.658
MT-CR | 0.965 | 0.955 | 0.925 | 0.925 | 0.954 | 0.967
MT-HCR | 0.975 | 0.993 | 0.999 | 0.989 | 0.982 | 0.992
MT-HCCR | 0.982 | 0.995 | 0.999 | 0.992 | 0.985 | 0.993
MT-HCCAR | 0.984 | 0.996 | 0.999 | 0.995 | 0.990 | 0.996
TABLE II: Model performance Comparison with Baseline Models and Ablation Study - regression results. Models are trained by OCI. Model | MSE | MSE each label | $R^{2}$ |  $R^{2}$ each label | FMG
---|---|---|---|---|---
|  | Liquid | Ice |  | Liquid | Ice | Liquid | Ice
Liu et al.,2022[8] | - | - | - | - | - | - | - | -
Chen et al., 2018[42] | 0.185 | 0.193 | 0.178 | 0.235 | 0.247 | 0.212 | 43.92% | 60.64%
SEQ | 0.126 | 0.133 | 0.118 | 0.440 | 0.274 | 0.548 | 4.00% | 41.11%
MT-CR | 0.062 | 0.075 | 0.049 | 0.728 | 0.661 | 0.788 | 44.46% | 68.21%
MT-HCR | 0.038 | 0.045 | 0.030 | 0.843 | 0.810 | 0.874 | 65.93% | 82.28%
MT-HCCR | 0.031 | 0.035 | 0.028 | 0.869 | 0.854 | 0.884 | 66.38% | 82.43%
MT-HCCAR | 0.026 | 0.032 | 0.021 | 0.891 | 0.869 | 0.914 | 67.96% | 82.70%

###  V-C Comparison with baseline models

Tables I and II present a comprehensive evaluation of our proposed model MT-
HCCAR against two baseline methods, demonstrating superior performance across
both classification and regression metrics. Specifically, in the assessment of
classification tasks, ACCbi pertains to cloud masking performance, while
$AU(\overline{PRC})_{w}$ and $AU(\overline{PRC})_{c}$ are indicative of cloud
phase classification performance.

Upon close examination of Table I, the Random Forest (RF) classifier, as
implemented by [8], attains superior ACCbi, $AU(\overline{PRC})_{w}$, and
$AU(\overline{PRC})_{c}$ in comparison to the SEQ model and the simplest
Multi-Task Learning (MTL)-based model, MT-CR, within the ablation study.
However, MTL-based models featuring the HC module including MT-HCR, MT-HCCR,
and MT-HCCAR, surpass the performance of the RF classifier in these metrics.

Turning attention to the Cloud Optical Thickness (COT) regression task, the
Multi-Layer Perceptron (MLP)-based baseline method introduced by [42] yields
results comparable to the SEQ model, but much worse than MTL-based models.
This observation underscores the efficiency of MTL-based models for of COT
regression.

###  V-D Ablation study

Figure 3: Scatter plot with PDF distribution of true COT values (y-axis) and
predicted COT values (x-axis).

The outcomes obtained from models SEQ through MT-HCCAR in Table I and Table II
underscore the efficacy of our MTL structure and the incorporated modules,
namely Hierarchical Classification (HC) and Cross-Attention Residual (CAR), in
enhancing both classification and regression performances. A comparative
analysis reveals notable improvements across all metrics when contrasting the
results of models through SEQ to MT-HCCAR. Specifically, when compared to the
MT-CR model, the introduction of hierarchical classification MT-HCR
demonstrates superior accuracy in binary classification ([cloudy, cloud-free])
and further refines the classification of liquid and ice cloud phases, thereby
contributing to enhancements in all regression metrics. The performance of the
model does not see any drop from MT-HCR to MT-HCCAR, even though the models
focus on improving regression by using CAR.

Moreover, the efficacy of the CAR module is validated in the comparison
between MT-HCCAR and MT-HCR, showcasing its utility in improving both
classification and regression outcomes. Furthermore, the adoption of the
cross-attention module is substantiated by the superior performance observed
in MT-HCCR compared to MT-HCCAR, underscoring the module’s effectiveness in
facilitating information exchange between hierarchical classification and
regression tasks.

In order to visually depict the enhancements facilitated by the HC module and
the CAR module, Figure 3 presents a scatter plot illustrating the distribution
of ground truth Cloud Optical Thickness (COT) values and the corresponding
predicted COT values for all instances within the test set. The integration of
the HC module and CAR module in our model yields more accurate predictions
compared to MT-CR and MT-HCR, as evidenced by two notable observations in the
figure. Firstly, the Probability Density Function (PDF) distribution of
predictions closely aligns with the true distribution, indicating an improved
fidelity to the actual COT values. Secondly, the scatter points exhibit a
greater concentration along the diagonal line, further corroborating the
model’s enhanced precision when incorporating the HC and CAR modules.

###  V-E Earth science evaluation

####  V-E1 Applicability for satellite mission

Satellite-derived cloud masks are generally evaluated either through
comparison against a human-labeled image set, or against ground-based,
airborne, or spaceborne observations from sensors that have much greater
sensitivity for cloud detection than the image in question (typically active
sensors such as radar and/or lidar). Accuracies tend to depend on surface type
and illumination conditions; for daytime scenes such as those simulated here,
some reported accuracies for NASA’s widely-used MODIS cloud mask are 0.850
(Ackerman et al.[46]), 0.778 (Wang et al. [47]), and 0.894 (Kotarba [48]).
Accuracy tends to be worse than average in polar regions and better for non-
polar land and ocean [49]. Most commonly it is optically thin clouds (e.g.
subvisual cirrus) that are missed. For instance, Ackerman et al. [46] found
that 90% of misclassified pixels in MODIS over a vegetated ground site in the
USA had a COT under 0.4. Chen et al. [42] developed an NN cloud mask for
MODIS, also based on RT simulations, and found an accuracy of 0.985 on
training data but lower accuracies in the range of 0.739-0.885 on real data
depending on region and time of year. Although real-world conditions are more
complex than a simulation, and each sensor’s measurement characteristics
result in different sensitivity to clouds, the results suggest our models’
accuracies are competitive with current approaches. Note, that the cloud phase
is less commonly validated in this way.

Obtaining a validation-grade COT measurement is exceptionally challenging, and
reference datasets for COT lack the comprehensive coverage available for other
cloud properties like cloud mask and cloud height[50]. Consequently, existing
comparisons primarily focus on assessing agreement levels among different
remotely sensed datasets. As previously mentioned, the PACE mission aims to
achieve an error rate of less than 25% (for liquid clouds) or 35% (for ice
clouds) in determining COT for clouds with true COT > 5, at least 65% of the
time. Our network’s performance indicates its potential to meet these
objectives. Notably, the larger allowable uncertainty for ice clouds in the
PACE mission is rooted in the expectation that accurate determination would be
more challenging due to significant variations in ice crystal shape, size, and
roughness in real-world scenarios. In contrast, our results demonstrate that
the predicted uncertainty for ice COT can be smaller than that for liquid COT.
Conventional approaches often assume specific ice crystal shapes and
roughness, leading to increased uncertainty in retrieved ice COT. Our findings
suggest that satellite measurements inherently contain information related to
these parameters, which the networks are adept at learning.

####  V-E2 Performance of different satellite sensors

The classification tasks exhibit high accuracies across all sensors, with
virtually indistinguishable performance from a scientific perspective. This
uniformity is unsurprising, despite the considerable variance in the number of
bands among sensors. The fundamental nature of cloud masking in Earth science
prompts many sensors to incorporate a common set of bands proven to be
effective for this purpose, with additional bands often designed for diverse
applications. For instance, OCI’s hyperspectral bands support the measurement
of different ocean plankton species [3], revealing subtle spectral differences
not readily discernible in multispectral data. In the regression task, OCI and
VIIRS outperform ABI, indicating the utility of additional bands in predicting
COT. The unexpectedly slightly superior performance of VIIRS over OCI may
suggest that OCI’s substantially larger feature space makes finding an optimal
solution during training more challenging or could be attributed to stochastic
variation. Adding training epochs for OCI may lead to better results.

##  VI Strategies for Selecting the Optimal Model

Based on the evaluation metrics and experiment results expounded in Section V,
our analysis extends to the application of the MT models to three different
sensors simulated datasets, namely OCI, VIIRS, and ABI. Given the diverse
evaluation metrics and the utilization of different satellite datasets, it is
imperative to identify the singular optimal model among our multi-task
learning-based models, namely MT-CR, MT-HCR, and MT-HCCAR, to be implemented
on different datasets.

The process of model selection is based on the four metrics: ACCbi,
$AU(\overline{PRC})_{w}$, MSE, and R2. The process includes 3 steps: 1)
employing K-fold cross-validation, 2) executing one standard error rule (1SE
rule), and 3) computing performance scores. Initially, K-fold cross-validation
is employed to compare the model’s absolute performance using the mean value
of each metric over K folds. This involves aggregating the mean values of all
four metrics across models trained on each dataset over K folds. Subsequently,
the one standard error rule is applied to assess the balance between model
absolute performance and computational complexity.

While the first and second steps prove effective in selecting the optimal
model within the confines of a single metric and a single dataset, extending
this evaluation to encompass all datasets and all metrics requires the
introduction of absolute performance scores ($P_{ab}$ and 1SE performance
scores), denoted as $P_{1SE}$. These proposed scores serve as comprehensive
indicators for selecting a single model that excels across diverse datasets
and evaluation metrics. Table III shows collected metrics and performance
scores for each model and each dataset. The following paragraphs will describe
K-fold cross-validation, one standard error rule, and our proposed performance
score, respectively.

We note our model selection method is also useful for other multi-task
learning applications due to the multi-objective nature of such applications
and our method’s capability of considering multiple metrics and datasets.

TABLE III: Model selection using K-fold cross-validation and 1SE rule. Dataset | Model | $\mu(Metric)$ | $SE(Metric)$ | 1SE Range | $P_{ab}^{\mathcal{M,D}}$ | $P_{1SE}^{\mathcal{M,D}}$
---|---|---|---|---|---|---
| Mean value of ACCbi over 10 folds |
OCI | MT-CR | 0.969 | 7.608E-04 | [0.969, 0.970] | -1.54% | 0
MT-HCR | 0.984 | 6.741E-04 | [0.983, 0.985] | -0.04% | 1
MT-HCCAR | 0.985 | 1.057E-03 | [0.983, 0.986] | 0.00% | 0
VIIRS | MT-CR | 0.976 | 6.866E-04 | [0.975, 0.976] | -1.39% | 0
MT-HCR | 0.989 | 5.027E-04 | [0.988, 0.989] | -0.04% | 1
MT-HCCAR | 0.989 | 4.961E-04 | [0.989, 0.990] | 0.00% | 0
ABI | MT-CR | 0.973 | 1.229E-03 | [0.972, 0.974] | -1.37% | 0
MT-HCR | 0.986 | 4.716E-04 | [0.986, 0.987] | -0.03% | 1
MT-HCCAR | 0.987 | 6.180E-04 | [0.986, 0.987] | 0.00% | 0
| Mean value of $AU(\overline{PRC})_{w}$ over 10 folds |
OCI | MT-CR | 0.953 | 1.60E-03 | [0.951, 0.954] | -4.29% | 0
MT-HCR | 0.995 | 2.83E-04 | [0.995, 0.995] | -0.04% | 0
MT-HCCAR | 0.996 | 2.08E-04 | [0.995, 0.996] | 0.00% | 1
VIIRS | MT-CR | 0.939 | 3.17E-02 | [0.908, 0.971] | -5.78% | 0
MT-HCR | 0.996 | 1.88E-04 | [0.996, 0.997] | -0.06% | 0
MT-HCCAR | 0.997 | 1.37E-04 | [0.997, 0.997] | 0.00% | 1
ABI | MT-CR | 0.969 | 2.19E-03 | [0.967, 0.971] | -2.75% | 0
MT-HCR | 0.996 | 1.62E-04 | [0.996, 0.996] | 0.00% | 1
MT-HCCAR | 0.996 | 1.67E-04 | [0.996, 0.996] | -0.02% | 0
| Mean value of MSE over 10 folds |
OCI | MT-CR | 0.055 | 2.916E-03 | [0.058, 0.052] | -100.00% | 0
MT-HCR | 0.034 | 5.970E-04 | [0.034, 0.035] | -25.27% | 0
MT-HCCAR | 0.027 | 4.910E-04 | [0.028, 0.027] | 0.00% | 1
VIIRS | MT-CR | 0.043 | 2.075E-03 | [0.045, 0.041] | -76.48% | 0
MT-HCR | 0.030 | 8.104E-04 | [0.031, 0.029] | -21.65% | 0
MT-HCCAR | 0.025 | 1.282E-03 | [0.026, 0.023] | 0.00% | 1
ABI | MT-CR | 0.057 | 1.519E-03 | [0.059, 0.056] | -78.86% | 0
MT-HCR | 0.038 | 5.714E-04 | [0.039, 0.038] | -19.48% | 0
MT-HCCAR | 0.032 | 3.121E-04 | [0.032, 0.032] | 0.00% | 1
| Mean value of R2 over 10 folds |
OCI | MT-CR | 0.758 | 1.24E-02 | [0.770, 0.746] | -14.25% | 0
MT-HCR | 0.854 | 2.55E-03 | [0.857, 0.851] | -3.39% | 0
MT-HCCAR | 0.884 | 2.35E-03 | [0.882, 0.886] | 0.00% | 1
VIIRS | MT-CR | 0.810 | 9.62E-03 | [0.800, 0.820] | -9.72% | 0
MT-HCR | 0.874 | 3.59E-03 | [0.870, 0.878] | -2.58% | 0
MT-HCCAR | 0.897 | 5.30E-03 | [0.892, 0.902] | 0.00% | 1
ABI | MT-CR | 0.751 | 7.44E-03 | [0.743, 0.758] | -13.25% | 0
MT-HCR | 0.840 | 2.20E-03 | [0.838, 0.843] | -2.91% | 0
MT-HCCAR | 0.866 | 1.51E-03 | [0.864, 0.867] | 0.00% | 1

###  VI-A K-fold cross-validation and 1SE rule

K-fold cross-validation is used to compare the absolute performance of each
model $\mathcal{M}$ trained by each dataset $\mathcal{D}$ based on mean values
of the four evaluation metrics over K folds. We use the number of folds $K=10$
because this value has been found through experiments to generally result in a
model skill estimate with low bias and modest variance [51, 44]. For each
fold, $(K-1)N/K$ pixels are in the training set, and the left $N/K$ pixels are
in the test set. Mean values $\mu(Metric)$ of four evaluation metrics over 10
folds are calculated for each experiment to compare differences in model
absolute performance. Among the four metrics, Accbi and
$AU(\overline{PRC})_{w}$ access the classification task while MSE and R2
access the regression task. The mean value of a certain metric $Metric$ over
10 folds $\mu(Metric)$ is:
$\mu(Metric)=\frac{1}{K}\sum_{k=1}^{K=10}Metric^{k,\mathcal{M}^{\mathcal{D}}}$,
where $Metric^{k}$ is $Metric$ for fold $k$ of model $\mathcal{M}$ trained by
dataset $\mathcal{D}$.

With the mean values from K-fold cross validation, we then use one standard
error rule (1SE rule) [52] to select the most parsimonious model for each
dataset. The mean of evaluation metrics, computed over 10 folds, provides an
absolute measure of model accuracy. However, it is essential to consider the
trade-off between model complexity and absolute performance, as this balance
is crucial for assessing the overall quality of a model. Therefore, in
addition to the mean value, we also calculate the region defined by one
standard error rule (1SE region) of each metric. With the 1SE range of each
metric, we need to select “1SE region of the best”, which is the 1SE range of
the model with minimum mean error. Then the simplest model whose mean falls
within 1SE region of the best is the most parsimonious model according to 1SE
rule [53, 43]. Under our experiment setting, the standard error for a certain
metric $SE(Metric)$ is calculated as:

$SE(Metric)=\frac{s}{\sqrt{K}}$ (19)

, where $s$ is the square root of variance $s^{2}$ for $Metric$ over 10 folds.
Then the 1SE region of each model is defined as: $[\mu-SE,\mu+SE]$.

In our experimental analyses, we derive the following conclusion regarding the
application of the 1SE rule: if none of the $\mu(Metric)$ values for
alternative models lie within the 1SE region of the optimal model, the model
exhibiting the highest mean metric value is the preferred choice based on this
specific metric.

###  VI-B The absolute score and the 1SE rule score

Despite conducting K-fold cross-validation and 1SE rule, we still need a
method to describe the general model performance across all four evaluation
metrics and all three datasets. Here we introduced the absolute performance
score and the 1SE rule score.

To describe the absolute performance based on the mean value as well as the
model selection results based on the 1SE rule for each model over all datasets
and all metrics, we propose to use absolute performance score $P_{ab}$ and 1SE
score $P_{1SE}$. The absolute performance score $P_{ab}$ describes how worse
other models are than the model with the best absolute performance. The
performance score for each model is calculated by the following formula:

$\begin{split}P_{ab}^{\mathcal{M}}&=\sum_{D}\sum_{Metric}P_{ab}^{\mathcal{M,D}}\\\
&=\sum_{D}\sum_{Metric}\frac{\mu({Metric})^{\mathcal{M}^{\mathcal{D}}}-\mu({Metric})^{\mathcal{M}^{\mathcal{D}\star}_{\mu}}}{\mu({Metric})^{\mathcal{M}^{\mathcal{D}\star}_{\mu}}}\end{split}$
(20)

, where components $P_{ab}^{\mathcal{M,D}}$ is the score for a specific model
$\mathcal{M}$ applied to a specific dataset $\mathcal{D}$,
$\mu({Metric})^{\mathcal{M}^{\mathcal{D}}}$ is the mean value of a specific
metric $Metric$ of model $\mathcal{M}^{\mathcal{D}}$ applied to dataset
$\mathcal{D}$, and $\mathcal{M}^{\mathcal{D}\star}_{\mu}$ is the model with
the best mean value for dataset $\mathcal{D}$.

The 1SE performance score is the weighted sum of counts that the model is
selected according to the 1SE rule across all metrics and all datasets. The
component $P_{1SE}^{\mathcal{M,D}}$ is the 1SE performance score for model
$\mathcal{M}$ applied to dataset $\mathcal{D}$. The calculation of the total
score $P_{1SE}$ for model $M$ is:

$\begin{split}P_{1SE}^{\mathcal{M}}&=\sum_{D}\mathbf{w_{m}}\sum_{Metric}P_{1SE}^{\mathcal{M,D}}\\\
&=\sum_{D}\sum_{Metric}\psi(\mu(Metric),\mu(\mathcal{M},\mathcal{D})\end{split}$
(21)

, where $\psi(\cdot)$ is a rule that if Model $\mathcal{M}$ trained by dataset
$\mathcal{D}$ is selected by 1SE rule based on the mean value of a certain
metric, then $\psi(\cdot)$ outputs 1, and otherwise 0. Weights
$\mathbf{w_{m}}=[w_{Acc_{bi}},w_{AU(\overline{PRC})_{w}},w_{MSE},w_{R^{2}}]$
define the significance of four metrics. We set $\mathbf{w_{m}}=[1,1,1,1]$ in
our experiments because we consider the classification task and regression
task equally important in our application. We will further discuss the model
performance shown in Table III based on 1SE region of the best to show that
MT-HCCAR is optimal based on all four metrics of ACCbi,
$AU(\overline{PRC})_{w}$, MSE, and R2.

###  VI-C Model selection results

This section discusses the superior performance of our proposed MT-HCCAR
model, evident in its highest absolute performance score
$P_{ab}^{\mathcal{M}}$ score and the top 1SE rule score
$P_{1SE}^{\mathcal{M}}$. The last two columns of Table III are components of
$P_{ab}^{\mathcal{M}}$ and $P_{1SE}^{\mathcal{M}}$. Utilizing Equation 20, the
absolute performance scores for each model are calculated as follows:
$P_{ab}^{{MT-CR}}=-4.12\%$, $P_{ab}^{{MT-HCR}}=-0.21\%$, and $P_{ab}^{{MT-
HCCAR}}=-0.02\%$. The highest $P_{ab}^{\mathcal{M}}$ indicates that MT-HCCAR
exhibits the best mean value across all evaluation metrics, including Accbi,
$AU(\overline{PRC})_{w}$, and R2, across all satellite datasets.

Furthermore, MT-HCCAR achieves the most substantial $P_{1SE}^{\mathcal{M}}$.
By using Equation 21, the 1SE performance scores for each model are calculated
as follows: $P_{1SE}^{{MT-CR}}=0$, $P_{1SE}^{{MT-HCR}}=4$, $P_{1SE}^{{MT-
HCCAR}}=8$. Achieving the highest $P_{1SE}^{\mathcal{M}}$ suggests that all
other models exhibit mean values beyond the 1 standard error (1SE) range of
MT-HCCAR. This implies that MT-HCCAR avoids overfitting in comparison to other
models, considering model complexity. This finding further establishes MT-
HCCAR as the optimal model for implementation across various satellite
datasets.

Additionally, MT-HCCAR exhibits higher mean values and a narrower 1SE range
compared to alternative models, indicating superior generalization. In
summary, MT-HCCAR, adhering to the 1SE rule, demonstrates superior performance
in both classification and regression tasks, surpassing other models without
succumbing to overfitting.

##  VII Conclusions

In this study, we present MT-HCCAR, an end-to-end multi-task learning model
tailored for cloud property retrieval on simulated OCI satellite dataset in
the PACE project, tackling tasks including cloud masking, cloud phase
classification, and Cloud Optical Thickness (COT) prediction. The model is
implemented on three sensors’ simulated datasets (OCI, VIIRS, and ABI),
respectively, to examine its generalization. Comparative analyses against two
baseline methods and ablation studies underscore the effectiveness of the
Hierarchical Classification (HC) module and the Classification-Assisted with
Attention-based Regression (CAR) module, enhancing performance in both
classification and regression tasks. Moreover, our model selection approach,
incorporating K-fold cross-validation, the 1SE rule, and two novel performance
scores, establishes MT-HCCAR’s superior performance across different datasets
and multiple evaluation metrics. The positive results affirm our model’s
capability to address real-world challenges in cloud property retrieval and
other multi-task applications. Future research endeavors will involve applying
the model to actual OCI data, co-located in space and time with VIIRS and ABI
post-PACE launch, to assess its performance, and consistency, and enable
detailed comparisons with traditional (non-machine-learning) approaches.

## Acknowledgment

This work is supported by a student fellowship from Goddard Earth Sciences
Technology and Research (GESTAR) II, UMBC, grant OAC–1942714 from the National
Science Foundation (NSF) and grant 80NSSC21M0027 from the National Aeronautics
and Space Administration (NASA).

## References

  * [1] Aiko Voigt, Nicole Albern, Paulo Ceppi, Kevin Grise, Ying Li, and Brian Medeiros.  Clouds, radiation, and atmospheric circulation in the present-day climate and under climate change.  Wiley Interdisciplinary Reviews: Climate Change, 12(2):e694, 2021.
  * [2] R. Hollmann, C. J. Merchant, R. Saunders, C. Downy, M. Buchwitz, A. Cazenave, E. Chuvieco, P. Defourny, G. de Leeuw, R. Forsberg, T. Holzer-Popp, F. Paul, S. Sandven, S. Sathyendranath, M. van Roozendael, and W. Wagner.  The esa climate change initiative: Satellite data records for essential climate variables”.  Bulletin of the American Meteorological Society, 94(10):1541 – 1552, 2013.
  * [3] P. Jeremy Werdell, Michael J. Behrenfeld, Paula S. Bontempi, Emmanuel Boss, Brian Cairns, Gary T. Davis, Bryan A. Franz, Ulrik B. Gliese, Eric T. Gorman, Otto Hasekamp, Kirk D. Knobelspiesse, Antonio Mannino, J. Vanderlei Martins, Charles R. McClain, Gerhard Meister, and Lorraine A. Remer.  The plankton, aerosol, cloud, ocean ecosystem mission: Status, science, advances.  Bulletin of the American Meteorological Society, 100(9):1775 – 1794, 2019.
  * [4] Ocean Color Instrument (OCI) Sensor.  https://oceancolor.gsfc.nasa.gov/data/pace/characterization/.
  * [5] Moderate Resolution Imaging Spec- troradiometer (MODIS) Sensor.  https://modis.gsfc.nasa.gov/.
  * [6] Visible Infrared Imaging Radiometer Suite (VIIRS) Sensor.  https://ncc.nesdis.noaa.gov/VIIRS/VIIRSSpectralResponseFunctions.php.
  * [7] Advanced Baseline Imager (ABI) Sensor.  https://ncc.nesdis.noaa.gov/GOESR/ABI.php.
  * [8] Chao Liu, Shu Yang, Di Di, Yuanjian Yang, Chen Zhou, Xiuqing Hu, and Byung-Ju Sohn.  A machine learning-based cloud detection algorithm for the himawari-8 spectral image.  Advances in Atmospheric Sciences, 39(12):1994–2007, 2022.
  * [9] Chenxi Wang, Steven Platnick, Kerry Meyer, Zhibo Zhang, and Yaping Zhou.  A machine-learning-based cloud detection and thermodynamic-phase classification algorithm using passive spectral observations.  Atmospheric Measurement Techniques, 13(5):2257–2277, 2020.
  * [10] Gianluca Giuffrida, Lorenzo Diana, Francesco de Gioia, Gionata Benelli, Gabriele Meoni, Massimiliano Donati, and Luca Fanucci.  Cloudscout: A deep neural network for on-board cloud detection on hyperspectral images.  Remote Sensing, 12(14):2205, 2020.
  * [11] Qibin He, Xian Sun, Zhiyuan Yan, and Kun Fu.  Dabnet: Deformable contextual and boundary-weighted network for cloud detection in remote sensing images.  IEEE Transactions on Geoscience and Remote Sensing, 60:1–16, 2021.
  * [12] Yikun Yang, Wenxiao Sun, Yulei Chi, Xing Yan, Hao Fan, Xingchuan Yang, Zhanshan Ma, Quan Wang, and Chuanfeng Zhao.  Machine learning-based retrieval of day and night cloud macrophysical parameters over east asia using himawari-8 data.  Remote Sensing of Environment, 273:112971, 2022.
  * [13] Quan Wang, Chen Zhou, Xiaoyong Zhuge, Chao Liu, Fuzhong Weng, and Minghuai Wang.  Retrieval of cloud properties from thermal infrared radiometry using convolutional neural network.  Remote Sensing of Environment, 278:113079, 2022.
  * [14] Xin Huang, Sahara Ali, Chenxi Wang, Zeyu Ning, Sanjay Purushotham, Jianwu Wang, and Zhibo Zhang.  Deep domain adaptation based cloud type detection using active and passive satellite data.  In 2020 IEEE International Conference on Big Data (Big Data), pages 1330–1337, 2020.
  * [15] Xin Huang, Chenxi Wang, Sanjay Purushotham, and Jianwu Wang.  Vdam: Vae based domain adaptation for cloud property retrieval from multi-satellite data.  In Proceedings of the 30th International Conference on Advances in Geographic Information Systems. ACM SIGSPATIAL, 2022.
  * [16] GitHub Repository of Multi-Task Deep Learning with Hierarchical Classification and Attention-based Regression for Cloud Property Retrieval (MT-HCCAR ).  https://github.com/AI-4-atmosphere-remote-sensing/MT-HCCAR.
  * [17] Sorour Mohajerani and Parvaneh Saeedi.  Cloud-net: An end-to-end cloud detection algorithm for landsat 8 imagery.  In IGARSS 2019-2019 IEEE International Geoscience and Remote Sensing Symposium, pages 1029–1032. IEEE, 2019.
  * [18] Zhenfeng Shao, Yin Pan, Chunyuan Diao, and Jiajun Cai.  Cloud detection in remote sensing images based on multiscale features-convolutional neural network.  IEEE Transactions on Geoscience and Remote Sensing, 57(6):4062–4076, 2019.
  * [19] Jacob Høxbroe Jeppesen, Rune Hylsberg Jacobsen, Fadil Inceoglu, and Thomas Skjødeberg Toftegaard.  A cloud detection algorithm for satellite imagery based on deep learning.  Remote sensing of environment, 229:247–259, 2019.
  * [20] Michal Segal-Rozenhaimer, Alan Li, Kamalika Das, and Ved Chirayath.  Cloud detection algorithm for multi-modal satellite imagery using convolutional neural-networks (cnn).  Remote Sensing of Environment, 237:111446, 2020.
  * [21] Jianhua Guo, Jingyu Yang, Huanjing Yue, Hai Tan, Chunping Hou, and Kun Li.  Cdnetv2: Cnn-based cloud detection for remote sensing imagery with cloud-snow coexistence.  IEEE Transactions on Geoscience and Remote Sensing, 59(1):700–713, 2020.
  * [22] Dan López-Puigdollers, Gonzalo Mateo-García, and Luis Gómez-Chova.  Benchmarking deep learning models for cloud detection in landsat-8 and sentinel-2 images.  Remote Sensing, 13(5):992, 2021.
  * [23] Melike Ilteralp, Sema Ariman, and Erchan Aptoula.  A deep multitask semisupervised learning approach for chlorophyll-a retrieval from remote sensing images.  Remote Sensing, 14(1):18, 2021.
  * [24] Guangbin Zhang, Xianjun Gao, Jinhui Yang, Yuanwei Yang, Meilin Tan, Jie Xu, and Yanjun Wang.  A multi-task driven and reconfigurable network for cloud detection in cloud-snow coexistence regions from very-high-resolution remote sensing images.  International Journal of Applied Earth Observation and Geoinformation, 114:103070, 2022.
  * [25] Jianjun Chen, Weihao Hu, Di Cao, Zhenyuan Zhang, and Zhe Chen.  A novel multi-task learning method with attention mechanism for wind turbine blades imbalance fault diagnosis.  In 2022 4th Asia Energy and Electrical Engineering Symposium (AEEES), pages 857–862. IEEE, 2022.
  * [26] Ron Sarafian, Dori Nissenbaum, Shira Raveh-Rubin, Vikhyat Agrawal, and Yinon Rudich.  Deep multi-task learning for early warnings of dust events implemented for the middle east.  Npj climate and atmospheric science, 6(1):23, 2023.
  * [27] Eliot Kim, Peter Kruse, Skylar Lama, Jamal Bourne, Michael Hu, Sahara Ali, Yiyi Huang, and Jianwu Wang.  Multi-task deep learning based spatiotemporal arctic sea ice forecasting.  In 2021 IEEE International Conference on Big Data (Big Data), pages 1847–1857. IEEE, 2021.
  * [28] Chunping Qiu, Lukas Liebel, Lloyd Haydn Hughes, Michael Schmitt, Marco Körner, and Xiao Xiang Zhu.  Multitask learning for human settlement extent regression and local climate zone classification.  IEEE Geoscience and Remote Sensing Letters, 19:1–5, 2020.
  * [29] Joëlle Hanna, Damian Borth, and Michael Mommert.  Physics-guided multitask learning for estimating power generation and co 2 emissions from satellite imagery.  IEEE Transactions on Geoscience and Remote Sensing, 2023.
  * [30] Seongho Ahn, Trung Duc Tran, and Jongho Kim.  Systematization of short-term forecasts of regional wave heights using a machine learning technique and long-term wave hindcast.  Ocean Engineering, 264:112593, 2022.
  * [31] Zhicheng Zhao, Ze Luo, Jian Li, Can Chen, and Yingchao Piao.  When self-supervised learning meets scene classification: Remote sensing scene classification based on a multitask learning framework.  Remote Sensing, 12(20):3276, 2020.
  * [32] Odele M Coddington, Erik C Richard, Dave Harber, P Pilewskie, Thomas N Woods, K Chance, Xiong Liu, and Kang Sun.  The tsis-1 hybrid solar reference spectrum.  Geophysical research letters, 48(12):e2020GL091709, 2021.
  * [33] Andrew M Sayer, Luca Lelli, Brian Cairns, Bastiaan van Diedenhoven, Amir Ibrahim, Kirk D Knobelspiesse, Sergey Korkin, and P Jeremy Werdell.  The chroma cloud-top pressure retrieval algorithm for the plankton, aerosol, cloud, ocean ecosystem (pace) satellite mission.  Atmospheric Measurement Techniques, 16(4):969–996, 2023.
  * [34] Claudia Emde, Robert Buras-Schnell, Arve Kylling, Bernhard Mayer, Josef Gasteiger, Ulrich Hamann, Jonas Kylling, Bettina Richter, Christian Pause, Timothy Dowling, et al.  The libradtran software package for radiative transfer calculations (version 2.0. 1).  Geoscientific Model Development, 9(5):1647–1672, 2016.
  * [35] Zhibo Zhang, Andrew S Ackerman, Graham Feingold, Steven Platnick, Robert Pincus, and Huiwen Xue.  Effects of cloud horizontal inhomogeneity and drizzle on remote sensing of cloud droplet effective radius: Case studies based on large-eddy simulations.  Journal of Geophysical Research: Atmospheres, 117(D19), 2012.
  * [36] Steven Massie, Heather Cronk, Aronne Merrelli, Sebastian Schmidt, and Steffen Mauceri.  Insights into 3d cloud radiative transfer for oco-2.  Atmospheric Measurement Techniques Discussions, 2022:1–40, 2022.
  * [37] Shubham Toshniwal, Hao Tang, Liang Lu, and Karen Livescu.  Multitask learning with low-level auxiliary tasks for encoder-decoder based speech recognition.  arXiv preprint arXiv:1704.01631, 2017.
  * [38] Ryohei Kuga, Asako Kanezaki, Masaki Samejima, Yusuke Sugano, and Yasuyuki Matsushita.  Multi-task learning using multi-modal encoder-decoder networks with shared skip connections.  In Proceedings of the IEEE International Conference on Computer Vision Workshops, pages 403–411, 2017.
  * [39] Fábio M Miranda, Niklas Köhnecke, and Bernhard Y Renard.  Hiclass: a python library for local hierarchical classification compatible with scikit-learn.  Journal of Machine Learning Research, 24(29):1–17, 2023.
  * [40] Ashish Vaswani, Noam Shazeer, Niki Parmar, Jakob Uszkoreit, Llion Jones, Aidan N Gomez, Łukasz Kaiser, and Illia Polosukhin.  Attention is all you need.  Advances in neural information processing systems, 30, 2017.
  * [41] Xiaolong Wang, Ross Girshick, Abhinav Gupta, and Kaiming He.  Non-local neural networks.  In Proceedings of the IEEE conference on computer vision and pattern recognition, pages 7794–7803, 2018.
  * [42] Nan Chen, Wei Li, Charles Gatebe, Tomonori Tanikawa, Masahiro Hori, Rigen Shimada, Teruo Aoki, and Knut Stamnes.  New neural network cloud mask algorithm based on radiative transfer simulations.  Remote Sensing of Environment, 219:62–71, 12 2018.
  * [43] Stuart J Russell and Peter Norvig.  Artificial intelligence a modern approach.  London, 2010.
  * [44] Gareth James, Daniela Witten, Trevor Hastie, Robert Tibshirani, et al.  An introduction to statistical learning, volume 112.  Springer, 2013.
  * [45] NASA PACE Validation Plan, 2020.  https://pace.oceansciences.org/docs/PACE_Validation_Plan_14July2020.pdf.
  * [46] Steven A Ackerman, RE Holz, R Frey, EW Eloranta, BC Maddux, and M McGill.  Cloud detection with modis. part ii: validation.  Journal of Atmospheric and Oceanic Technology, 25(7):1073–1086, 2008.
  * [47] Tao Wang, Eric J Fetzer, Sun Wong, Brian H Kahn, and Qing Yue.  Validation of modis cloud mask and multilayer flag using cloudsat-calipso cloud profiles and a cross-reference of their cloud classifications.  Journal of Geophysical Research: Atmospheres, 121(19):11–620, 2016.
  * [48] A. Z. Kotarba.  Calibration of global MODIS cloud amount using CALIOP cloud profiles.  Atmospheric Measurement Techniques, 13(9):4995–5012, 2020.
  * [49] RE Holz, SA Ackerman, FW Nagle, R Frey, S Dutcher, RE Kuehn, MA Vaughan, and B Baum.  Global moderate resolution imaging spectroradiometer (modis) cloud detection and height evaluation using caliop.  Journal of Geophysical Research: Atmospheres, 113(D8), 2008.
  * [50] Steven Platnick, Kerry Meyer, Galina Wind, Robert E Holz, Nandana Amarasinghe, Paul A Hubanks, Benjamin Marchant, Steven Dutcher, and Paolo Veglio.  The nasa modis-viirs continuity cloud optical properties products.  Remote sensing, 13(1):2, 2020.
  * [51] Max Kuhn, Kjell Johnson, et al.  Applied predictive modeling, volume 26.  Springer, 2013.
  * [52] Leo Breiman.  Classification and regression trees.  Routledge, 1984.
  * [53] Jerome Friedman, Trevor Hastie, and Rob Tibshirani.  Regularization paths for generalized linear models via coordinate descent.  Journal of statistical software, 33(1):1, 2010.