# Iterative-Bayesian unfolding of isotopic cosmic-ray fluxes measured by
AMS-02

E. F. Bueno<EMAIL_ADDRESS>F. Barão<EMAIL_ADDRESS>M. Vecchi
<EMAIL_ADDRESS>Kapteyn Astronomical Institute, University of
Groningen,Landleven 12, 9747 AD, Groningen, the Netherlands  Laboratório de
Instrumentação e Física Experimental de Partículas (LIP), 1649-003 Lisboa,
Portugal

###### Abstract

The measurement of the isotopic composition of cosmic rays (CRs) provides
essential insights the understanding of the origin and propagation of these
particles, namely the CR source spectra, the propagation processes and the
galactic halo size. The Alpha Magnetic Spectrometer (AMS-02), a CR detector
operating aboard the International Space Station since May 2011, has the
capability of performing these measurements due to its precise determination
of the velocity provided by its Time of Flight (TOF) and Ring Imaging
Cherenkov (RICH) detector. The correct interpretation of the data requires the
measurements to be deconvoluted from the instrumental effects. The unique
design of AMS-02, with more than one subdetector being used to measure the
same flux, requires a novel approach to unfold the measured fluxes. In this
work, we describe an iterative-Bayesian unfolding method applied in the
context of isotopic flux measurements in AMS-02. The accuracy of the method is
assessed using a simulated flux based on previous measurements and a full
detector response function. We introduce a non-parametric regularization
method for the detector response functions, as well as a single, smooth prior
flux covering the full range of measurements from both detectors, TOF and
RICH. In addition, the estimation of the errors and a discussion about the
performance of the method are also shown, demonstrating that the method is
fast and reliable, allowing for the recovery of the true fluxes in the full
energy range.

###### keywords:

AMS-02 , unfolding , bayesian , cosmic-rays , isotopes

††journal: Nuclear Inst. and Methods in Physics Research, A

##  1 Introduction

The elemental abundance of light isotopes such as 2H, 6Li and 10Be in the
Solar System is low if compared to their abundance in CRs, which indicates
these are particles from secondary origin: they are produced during the
interaction of primary CRs, those produced and accelerated at the sources,
with the insterstellar medium (ISM) [1]. Therefore, the measurement of the
flux of such isotopes provides essential information about the CR source
composition and the propagation processes in our Galaxy [2]. Moreover,
radioactive species such as 10Be also provide important constraints on the
galactic halo size [3]. Magnetic spectrometers such as AMS-01[4], PAMELA[5]
and IMAX[6] have measured the fluxes of some of these isotopes. In such
detectors, the mass identification is performed by combining the measurements
of the rigidity (momentum per unit charge, $R=pc/Ze$) and the velocity,
measured either directly or through energy deposition, $\beta$. Both
measurements can be combined in a single observable, the mass, which is given
by the equation

$m=\frac{RZ}{\beta\gamma}$ (1)

where $Z$ is the magnitude of the charge and $\gamma$ is the Lorentz factor.
The measurements before AMS-02, reaching at most a few GeV/n, are limited by
detector velocity and rigidity resolutions. In AMS-02, isotopes are identified
either in the velocity and rigidity phase space [7], or through the use of
mass templates in events binned in velocity [8].

Additionally to particle identification, one of the main challenges in
particle physics experiments is the reconstruction of the true properties of
the particles. The finite resolution of the detector, energy losses and other
instrumental effects lead to biasing and smearing of quantities such as
velocity and rigidity. In the case of the measurement of isotopic spectra,
where they are generally obtained as a function of the kinetic energy per
nucleon or rigidity, one must face this problem to recover the spectra of the
particles before the particles interact with the detector. Generally speaking,
the true value of a variable in a given bin, $x_{t}^{j}$, is related to its
measured value, $x_{m}^{i}$, through the response function of the detector,
$R(x_{m}^{i}|x_{t}^{j})$:

$x_{m}^{i}=\sum_{j}R(x_{m}^{i}|x_{t}^{j})x_{t}^{j}$ (2)

The problem is then to invert the previous equation, obtaining the values of
$x_{t}^{j}$. Naively inverting this expression leads to numerical instability;
hence other methods have been developed to obtain estimators of $x_{t}$. This
estimation procedure is frequently called unfolding and is widely used in the
field of astroparticle physics, especially in direct measurements of CR
fluxes, where several unfolding techniques have been applied in different
measurements [5][6][7]. In this paper, we will focus on the application of an
iterative approach, based on Bayes’ theorem, similar to the works of
D’Agostini [9] and Mülthei and Schorr [10], in the context of the measurement
of isotopic fluxes with AMS-02, where the events are counted as a function of
the velocity measured by two different subdetectors, TOF and RICH.

This paper is organized as follows: in section 2 we given and overview of
AMS-02, describing the TOF and RICH in more detail. In section 3 we recall the
ingredients for the calculation of the flux, followed by the development of
the flux unfolding expressions based on Bayes’ theorem. In section 4 we
describe in detail the application of the proposed unfolding technique in mock
fluxes in AMS-02, from the construction of the response function to the error
estimation. We conclude in section 5.

##  2 The Alpha Magnetic Spectrometer

The Alpha Magnetic Spectrometer (AMS-02) is a particle physics detector
operating aboard the International Space Station since May 2011, having
collected over 200 billion events so far. AMS-02 consists of different sub-
detectors[11]: the nine-layer silicon tracker, combined with the 0.15 T
permanent magnet, measures the magnitude and sign of the charge, and the
rigidity ($R=pc/Ze$) of the particle; the Transition Radiator Detector, used
to separate leptons from hadrons; the Time of Flight (TOF), responsible for
measuring the velocity of the particle, also acting as the main trigger of
AMS-02; the Ring Imaging Cherenkov Detector (RICH), which measures the
magnitude of the charge, as well as the velocity; the Anti-Coincidence Counter
(ACC), that rejects particles with high-incidence angle; and the
Electromagnetic Calorimeter (ECAL), located right below the RICH, used to
measure the energy of the particle, as well as to separate leptons from
hadrons.

AMS-02 can measure the fluxes of isotopes in CRs in a wide energy range
because of the accuracy and complementarity of the velocity measurements
provided by the TOF and the RICH. The TOF consists of two pairs of layers of
scintillator planes, one above and the other below the magnet. The time taken
for a particle to travel between the pairs of scintillators while traversing
the magnetic field, $\Delta t$, is registered. The length of the path taken by
the particle, $L$, is obtained with the track reconstruction made by the
silicon tracker. In possession of these quantities, the velocity $\beta$ is
simply $\beta=\Delta t/cL$, and has a resolution $\Delta\beta/\beta^{2}\approx
4\%$ for single-charged particles and $\beta\approx 1$ [12]. The RICH is a
proximity-focusing detector which measures $\beta$ using the Cherenkov effect.
It includes a radiator plane in which particles travelling with a velocity
above the speed of light in the medium emit electromagnetic radiation peaked
in the UV region [13]. The radiation emitted is detected below in a plane of
$4\times 4$ pixelized photomultipliers. The PMTs of the detection plane are
arranged to be outside the ECAL acceptance, thus reducing the amount of matter
traversed by the particles before they reach the calorimeter. Considering this
arrangement, the radiator plane is made of two materials; its center consists
of sodium fluoride (RICH-NaF), with a refraction index $n=1.33$ and threshold
of emission $\beta=0.75$. The NaF produces broad Cherenkov rings which avoid
the ECAL hole, increasing the acceptance of the detector. Surrounding this
material there is the aerogel (RICH-AGL) with a refraction index $n=1.05$ and
threshold of emission $\beta=0.96$. the velocity resolution in the NaF and in
the AGL are $\Delta\beta/\beta\approx 0.35\%$ and $\Delta\beta/\beta\approx
0.12\%$, respectively, for particles with $Z=1$ and $\beta\approx 1$ [11].

Monte Carlo (MC) simulations of protons were also used in this study. These
simulations were produced by the AMS collaboration through a dedicated
software based on the GEANT4 package [14]. The software simulates the
interactions of particles with AMS material. It produces detector responses,
which are then used to reconstruct the desired event properties in the same
way as in data.

##  3 Flux unfolding

The flux of a particle can be expressed as

$\Phi_{m}(x)=\frac{N(x)}{A(x)\varepsilon(x)T(x)\Delta x}$ (3)

where $N(x)$ is the number of events in the accumulated number of events in a
finite interval of $x$, $\varepsilon(x)$ the selection efficiency, $A(x)$ the
geometric acceptance of the detector, $T(x)$ is the exposure time for
particles above the geomagnetic cutoff; and $\Delta x$ the interval width. In
the flux expression, $N(x)$ is an observable and therefore is a function of
the measured variable, $x_{m}$; $A$ is a constant factor; $\varepsilon(x)$ is
estimated from Monte Carlo simulations. Hence it is obtained as a function of
the true variable, $x_{t}$; $T(x)$ is, strictly speaking, dependent on the
measured variable, but due to the use of a detector resolution safety factor,
as done in [15], this dependence vanishes.

Since $N(x)$ includes the instrumental effects, the obtention of the true
number of events requires the application of an unfolding method. The
following expression relates the measured number of events in a given bin,
$N(x_{m}^{i})$, to its true value, $N(x_{t}^{i})$

$N(x_{m}^{i})=\sum_{j}R(x_{m}^{i}|x_{t}^{j})\varepsilon(x_{t}^{j})\eta_{0}(x_{t}^{j})N(x_{t}^{j})$
(4)

where the factor $R(x_{m}^{i}|x_{t}^{j})$ is the response function of the
detector, and $\eta_{0}(x_{t}^{i})$ is an efficiency translating the
probability of a true value being observed inside the detection range, that
is, $\eta_{0}(x_{t}^{j})=\sum_{i}R(x_{m}^{i}|x_{t}^{j})$. The efficiency
$\varepsilon(x_{t}^{i})$, although estimated from simulated events as a
function of the true variable, $x_{t}$, needs to be corrected for eventual
discrepancies between data and simulations on the efficiencies of selection
criteria, that is

$\varepsilon(x_{t})\rightarrow\varepsilon(x_{t})\underbrace{\prod_{k}\frac{\varepsilon_{Data}(x_{m}^{k})}{\varepsilon_{MC}(x_{m}^{k})}}_{\delta(x_{m})}$
(5)

where $\varepsilon_{MC}$ and $\varepsilon_{Data}$ are, respectively,
efficiencies as estimated from simulations and experimental data for the set
($k$) of data selection criteria. This ensures that the obtained efficiency
values agree with those from the actual experiment. Hence, equation 4 becomes

$\frac{N(x_{m}^{i})}{\delta(x_{m}^{i})}=\sum_{j}R(x_{m}^{i}|x_{t}^{j})\varepsilon(x_{t}^{j})\eta_{0}(x_{t}^{j})N(x_{t}^{j})$
(6)

The goal of the unfolding procedure is to obtain the true flux,
$\Phi_{t}(x_{t})$, which corresponds to the incident flux purged from
instrumental effects. To this end, the true number of events, $N(x_{t})$ must
be obtained. The inversion of equation 6 is intrinsically unstable and
therefore cannot be done directly. Using Bayes’ theorem, one can estimate the
probability of the true value being $x_{t}$ in the bin $i$, when $x_{m}$ was
measured

$P(x_{t}^{i}|x_{m}^{j})=\frac{P(x_{m}^{j}|x_{t}^{i})P(x_{t}^{i})}{\sum_{k}P(x_{m}^{j}|x_{t}^{k})P(x_{t}^{k})}$
(7)

where he probability $P(x_{m}^{j}|x_{t}^{i})$ represents the likelihood that
$x_{m}$ in bin $j$ was produced by $x_{t}$ in bin $i$, while $P(x_{t}^{i})$
represents the prior knowledge of $x_{t}$. As a result, the probability
$P(x_{t}^{i}|x_{m}^{i})$ is the posterior, or the updated knowledge of
$x_{t}$. In the context of this problem, the likelihood is the response
function of the detector $R(x_{m}|x_{t})$, and the prior is the previous
knowledge about $N(x_{t})$. With an estimate of the response function, and
with some prior knowledge of $N(x_{t})$, it is possible to use equation 7 to
estimate the true value from the measured one. This is done in the following
equation, where the true number of events in the bin $i$ is derived from the
accumulated number of events, the efficiency corrections and the posterior
probability.

$\hat{N}_{t}(x_{t}^{i})=\frac{1}{\varepsilon(x_{t}^{i})\eta_{0}(x_{t}^{i})}\sum_{j}\frac{N_{m}(x_{m}^{j})}{\delta(x_{m}^{j})}P(x_{t}^{i}|x_{m}^{j})$
(8)

$\hat{N}_{t}(x_{t}^{i})$ represents the estimator of $N(x_{t}^{i})$, and is
called the unfolded number of events. The terms outside the sum represent the
efficiency corrections. Since the posterior is naturally dependent on the
prior, a convergent iterative process has to be established, with
$\hat{N}_{t}(x_{t}^{i})$ being used to define a prior for the next iteration.
In each step, the prior is regularized by means of a fit to the data. The
convergence of the method is evaluated through the use of the $\xi$ estimator:

$\xi_{k+1,k}=\sum_{i}\left(\frac{\hat{N}_{t}^{k+1}(x_{t}^{i})-\hat{N}_{t}^{k}(x_{t}^{i})}{\hat{N}_{t}^{k}(x_{t}^{i})}\right)^{2}$
(9)

where $k$ is the number of the iteration.

In the next section, the method will be applied in the context of the AMS-02
experiment using a simulated flux that mimics the actual conditions of the
experiment.

##  4 Application to mock AMS-02 fluxes

In order to show the method at work in the context of single-charged isotopes
flux measurements in AMS-02, it has been applied to mock fluxes. Considering
the precise measurements of the velocity provided by the TOF and the RICH
detectors, the fluxes are measured as a function of the velocity in three
overlapping ranges: $0.5<\beta<0.9$ in the TOF; $0.78<\beta<0.99$ in the RICH-
NaF; and $0.96<\beta<0.996$ in the RICH-AGL. The events were generated
following the AMS-02 proton spectrum [15] parametrization given in [16]. They
were then smeared and biased according to the three detector response
functions. In order to reproduce the effect of data selection efficiencies,
generated events were weighted according to the selection efficiencies in
every range. Figure 1 shows the measured and generated fluxes as a function of
kinetic energy per nucleon, $E_{k}=(\gamma-1)m_{p}$, where $\gamma$ is the
Lorentz factor and $m_{p}$ is the proton mass. The noticeable differences
between the true and measured fluxes emphasize the necessity of using a method
to correct the measurements for the detector response.

Figure 1: True (empty) and measured (solid) fluxes versus the kinetic energy
per nucleon for the TOF (circles), RICH-NaF (squares) and RICH-AGL
(triangles).

In the following subsections, we describe the construction of the posterior
probability for the unfolding procedure.

###  4.1 Detector response function regularization

The detector response functions, often referred to as migration matrices, are
obtained from the Monte Carlo simulations of the detector. As in the isotopic
flux measurements, the events are counted in bins of velocity (or,
equivalently, kinetic energy per nucleon). These matrices correspond to binned
distributions of $\beta_{t}$, the velocity generated in the simulation, versus
$\beta_{m}$, the velocity as reconstructed in the detector. Figure 2 shows
these response functions for the TOF (a), RICH-NaF (b) and RICH-AGL (c). The
TOF has a larger velocity bias at low energies due to the energy loss
undergone by the particles. Both RICH-NaF and RICH-AGL are less penalized by
this effect.

Figure 2: Migration matrices of the TOF (a), RICH-NaF (b) and RICH-AGL (c).
See text for discussion.

Given that the simulations have a finite number of events, regularising the
matrices is a crucial step to avoiding fluctuations caused by the statistical
uncertainties. To this end, the one-dimensional projections used for
calculating the likelihood are regularized to obtain smooth distributions.
Figure 3 shows examples of these projections for the TOF, RICH-NaF, and RICH-
AGL in panels (a), (b), and (c), respectively, for three velocity intervals.
The $\beta_{m}$ distribution in the TOF is composed of a Gaussian core and a
tail corresponding to events whose velocity reconstruction was heavily
affected by the energy loss. For the two RICH regions, tails on the measured
velocity are visible for two main reasons: the low number of Cherenkov photons
in the ring produced by single-charged particles and the presence of
additional hits in the RICH PMT plane caused by secondary particles from
protons interacting with the AMS-02 material. In order to have a smooth
description of the detector response functions, we proceeded to their
regularization by using the kernel-density-estimation technique (KDE)[17].
This non-parametric method eases the smoothing of features that would need
complex parametrizations to be described, albeit at a higher computational
cost.

Figure 3: The resolution functions for the TOF (a), RICH-NaF (b) and RICH-AGL
(c) in different true velocity bins. See text for discussion.

The obtained parametrizations, $f(\beta_{m}|\beta_{t})$, were then used to
calculate the probabilities $P(\beta_{m}^{j}|\beta_{t}^{i})$ by integrating
them over the bins of measured velocity, $\beta_{m}^{j}$, for a given bin of
true velocity, $\beta_{t}^{i}$:

$P(\beta_{m}^{j}|\beta_{t}^{i})=\int_{\beta_{m}^{j}}f(\beta_{m}|\beta_{t})d\beta_{m}$
(10)

This process is repeated for every bin, allowing for the construction of a
smooth migration matrix.

###  4.2 Prior construction

The prior, representing the degree of belief in the true variable, is a
critical component in calculating the posterior defined in equation 7. Due to
the iterative nature of this method, the choice first prior cannot affect the
result. Hence, a natural choice for the initial prior is a flat shape. After
every iteration, the unfolded number of events, $\hat{N}(E_{k}(\beta)_{t})$,
is calculated using equation 8 for every detector range, allowing for the
calculation of the corresponding flux according to the equation

$\Phi(E_{k})_{t}=\frac{\hat{N}(E_{k})_{t}}{A(E_{k})_{t}\ T(E_{k})_{t}\ \Delta
E_{k}}$ (11)

These fluxes are regularized through a fit of the model described in [16],
performed on the three ranges simultaneously as shown in figure 4. This
parametrization represents an updated degree of belief in the flux as a
function of the true variable and is used to obtain the prior for the next
iteration through equation 11: the updated flux is used together with the
acceptance, exposure time and bin width to compute the new prior probability
of $N(E_{k})_{t}$. The advantage of such an approach is that it describes the
data well, providing a smooth description of the flux, even in transition
zones between the ranges, where resolution and threshold effects may cause the
points to fluctuate, as seen in figure 4. This guarantees that the prior for
the next iteration is a smooth, continuous function. In addition, it ensures
that the unfolded fluxes match in the transition zones between the detectors.
It prevents the fluctuations that may come from the unfolding procedure from
propagating to the next iteration, ensuring the method converges.

Figure 4: Model (black dashed line) adjusted on the results of the second
iteration of the TOF (pink circles), RICH-NaF (blue squares) and RICH-AGL
(green triangles).

###  4.3 Error estimation

The calculation of the unfolded number of events depends directly on the
detector response function, as well as on the number of measured events. These
quantities have errors associated with them, which are propagated to the
result of the unfolding. In order to estimate the contributions from each of
these quantities, a bootstrapping method was implemented. This method consists
of generating $n$ randomized sets of data based on the original one and
performing the unfolding procedure in all sets, resulting in a distribution of
unfolding results with $n$ entries per bin. This allows for calculating the
covariance matrices, which contain the total uncertainties in each bin in its
diagonal. Considering iterative nature of the procedure, where the results of
one iteration are used as the prior for the next, the propagation of the
uncertainty is done by sampling the prior.

Assuming that the uncertainties on the number of measured events and on the
migration matrices are uncorrelated, they are treated separately. For the
contribution from the number of measured events, $10^{3}$ sets of measurements
were generated by sampling new values $N(E(\beta)_{m}^{i})$ in every bin
according to a Gaussian distribution. That is, the number of events in each
bin $i$ of each new set, $N^{\prime}(E(\beta)_{m}^{i})$, is obtained by
sampling a Gaussian distribution $G$ centered at $N(E(\beta)_{m}^{i})$, with
RMS given by
$\sqrt{N(E(\beta)_{m}^{i})}$, that is

$N^{\prime}(E(\beta)_{m}^{i})=G(N(E(\beta)_{m}^{i});\sigma_{N})$ (12)

The unfolding procedure was applied to the complete set of measurements,
yielding a distribution of unfolded counts. Using the distribution of the
results, the covariance matrix of the unfolded number of events was obtained.

After every iteration, the prior is updated with the previously obtained
unfolded counts. We associate with the prior the uncertainty obtained in the
previous iteration and propagate it by sampling the unfolded counts according
to their uncertainty, analogously to what is stated in equation 12.

A similar exercise was performed to obtain the errors coming from the
uncertainties on the migration matrices. This time, the measured number of
events was kept fixed, and $10^{3}$ migration matrices were generated
according to

$P(E(\beta)_{m}^{j}|E(\beta)_{t}^{i})=G(P(E(\beta)_{m}^{j}|E(\beta)_{t}^{i});\sigma_{P})$
(13)

where G is a Gaussian distribution centered at the original
$P(E(\beta)_{m}^{j}|E(\beta)_{t}^{i})$, and width given by $\sigma_{P}$, that
represents the uncertainty on each migration matrix bin. Analogously to the
previous case, the unfolding is performed with the $10^{3}$ different
matrices, resulting in a covariance matrix unfolded counts per velocity range,
from where the uncertainties were obtained.

The relative errors for both contributions, coming from the diagonal of the
covariance matrices, are shown for the different ranges and iterations in
figure 5. Panel (a) shows the uncertainty coming from the migration matrices.
In all cases, the errors are larger in the extremes of each range. The
explanation for the increase at the end is the same for the three
subdetectors: at the upper edges of each velocity range, the velocity
resolution degrades, hence increasing the migration of events. Therefore,
fluctuations in the migration matrices lead to fluctuations in the result. The
increase at the beginning of each range is due to different physical
phenomena: in the TOF, it is due to the energy loss, which produces a strong
bias in the response function at the lowest energies. Hence any response
fluctuations in that region of the matrix produce significant changes in the
unfolded results. The errors on the RICH follow a similar trend, but in this
case this behavior is explained by the Cherenkov emission kinetic energy per
nucleon threshold. Events close to the threshold have a smaller probability of
being reconstructed. Hence the migration matrix is less populated in this
region. As a result, the dispersion of the unfolding results increases,
increasing the error. The same behavior is seen in panel (b), where the
uncertainty coming from the number of events is shown.

The errors from both contributions increase progressively with the number of
iterations, as expected from the error propagation done via the prior
sampling. The error increases faster from iteration 0 to 1 than from 1 to 2.
This is due to the change in the prior: while iteration 0 has a prior which is
flat in the entire energy range, iterations 1 and 2 use a model which is
fitted to the results, leading to larger uncertainties.

Once both contributions are estimated, the total error of the unfolding
procedure is given by the sum in quadrature of the uncertainties from the
migration matrix and the number of events.

Figure 5: Relative error on unfolded flux for the migration matrix
contribution (top) and on the number of measured events (bottom). Pink lines
represent the TOF, blue the NAF, and green the AGL. Solid lines are the errors
from iteration 0 (flat prior), while dotted and dot-dashed correspond to the
errors for iterations 1 and 2, respectively.

###  4.4 Number of iterations

Considering the iterative nature of this method, the number of iterations is
an important parameter. Too few iterations lead to a result far from the true
value, while a large number of iterations leads to a larger error due to the
error propagation. Therefore, one must stop iterating as soon as the
differences between two consecutive iterations are small enough. Figure 6
shows the residuals between consecutive iterations, defined as

$R^{k+1,k}_{i}=\frac{N^{k+1}_{i}-N^{k}_{i}}{N^{k+1}_{i}}$ (14)

where $k$ is the number of the iteration and $i$ is the bin index. In
addition, the value of $\xi$ defined in equation 9 is shown for every range in
each of the panels. The largest changes are between iteration 0 and 1 due to
the use of a flat prior in iteration 0. Iterations 1 and 2 already show signs
of convergence, with the residuals displaying small fluctuations in the last
points of the AGL range. The bottom panel shows that iterations 3 and 2 yield
essentially the same results, with $\xi_{3,2}\approx 10^{-8}$ meaning that it
is enough to stop at iteration 2.

The similarity between iterations 3 and 2 highlights the importance of
regularizing the prior in each iteration. The spline fit removes eventual
fluctuation that appear in the calculations and ensures convergence, whereas
the methods with no regularization have oscillations due to positive feedback
of these fluctuations [9].

Figure 6: Residuals and $\xi$ calculated for consecutive iterations in the
three velocity ranges. The top panel shows the residuals between iterations 0
and 1, the mid panel shows between 1 and 2, and the bottom shows between 3 and
2. The different ranges are represented by different markers and colors: pink
circles (TOF), blue squares (NAF) and green triangles (AGL). See text for
discussion.

###  4.5 Comparison with true flux

The top panel of figure 7 shows the unfolded flux after three iterations,
calculated with equation 11, for the different ranges, while in the bottom,
the residuals comparing the results with the true fluxes are shown. Generally,
the agreement between the unfolded and true values is very good, being below
$1\%$ in the entire measurement range. Minor fluctuations are seen at the
beginning of the RICH-NaF and RICH-AGL ranges. These are related to steep
drops in the number of measured events due to the efficiency of the detector,
which makes the true distribution harder to obtain in these regions. In the
case of the RICH the efficiencies drop rapidly below 0.8 and 3 GeV/n,
respectively, due to the emission threshold of Cherenkov radiation. Still, the
method can recover the true distribution quite well in these regions. There is
also a good agreement between the detectors in the overlapping regions.

Figure 7: True and unfolded flux after three iterations (top) and the
normalized residuals (bottom). See text for discussion.

##  5 Conclusion

The study of the isotopic composition of CRs will provide important
information about the propagation of CRs in our own Galaxy, as well as about
the CR source composition and the Galactic halo size. These measurements are
intrinsically challenging due to several factors: large backgrounds, as in the
case of single-charged isotopes, and small mass difference, as is the case of
Be isotopes. As a consequence, the available measurements are limited to the
kinetic energy per nucleon range from sub-GeV/n to a few GeV/n. The AMS-02
experiment, with its precise measurements of the rigidity and the velocity,
especially provided by the Cherenkov detector, has the capability of measuring
the fluxes of isotopes reaching energies up to 10 GeV/n. The correct
interpretation of these data depends on correcting the results for the
detector response. The unique design of AMS-02, with three different
measurements of the velocity, requires special treatment. In this work, we
applied an iterative unfolding method based on Bayes’ theorem in the context
of the isotopic flux measurements in AMS-02. The general mathematical aspects
of the method were discussed, emphasizing the importance of data/simulation
and efficiency corrections. In order to demonstrate the applicability and
performance of the method, a mock flux was generated according to the measured
proton spectrum, which was then smeared and biased according to the response
function of the different sub-detectors included in the analysis.
Regularization was also introduced in two stages of the method. The response
function was regularized with KDEs, which could readily describe the non-
Gaussian features of the measurements. In addition, the unfolding of the flux
in the three different ranges was connected via a common prior, obtained
through the fit of a model using all the available points. This fit leads to
the smoothing of features in the flux, such as the transition between ranges,
and prevents fluctuations from leading the method to calculation to diverge.
The error, estimated through a bootstrapping procedure, has shown strong
dependence on the migration matrix, growing rapidly in regions with high bias
or low statistics. Finally, the results have shown that only a few iterations
of this procedure are enough to recover the true spectrum within the
uncertainties, enabling the true flux to be estimated up to 10 GeV/n. The fast
convergence of this method also makes it suitable for studying the time-
variability of isotopic fluxes, where a large number of flux unfolding
calculations are necessary.

## CRediT authorship contribution statement

E. F. Bueno: Conceptualization, Methodology, Software, Data Curation, Formal
analysis, Visualization, Writing - Original draft preparation. F. Barao:
Conceptualization, Methodology, Software, Writing - Review $\&$ Editing. M.
Vecchi: Writing - Review $\&$ Editing, Funding acquisition, Visualization,
Supervision.

## Declaration of competing interest

The authors declare that they have no known competing financial interests or
personal relationships that could have appeared to influence the work reported
in this paper.

## Acknowledgements

This publication is part of the project ”Statistical methods applied to cosmic
ray anti-deuteron searches with the AMS-02 experiment” with project number
040.11.723 of the research programme Bezoekersbeurs 2019 VW, which is financed
by the Dutch Research Council (NWO). This publication is also part of the
project ”Title Cosmic ray antideuterons as a probe for new physics” with
project number OCENW.KLEIN.387/11680 of the research programme Grant Open
Competition Domain Science - M OC ENW KLEIN which is (partly) financed by the
Dutch Research Council (NWO). We are grateful to Portuguese FCT for the
financial support through grant CERN/FIS-PAR/0007/2021.

## References

  * [1] T. K. Gaisser, R. Engel, E. Resconi, Cosmic Rays and Particle Physics, 2nd Edition, Cambridge University Press, 2016.  doi:10.1017/CBO9781139192194.
  * [2] B. Coste, et al., Constraining galactic cosmic-ray parameters with $z\leq 2$ nuclei, A&A 539 (2012) A88.  doi:10.1051/0004-6361/201117927.
URL https://doi.org/10.1051/0004-6361/201117927

  * [3] S. Hayakawa, K. Ito, Y. Terashima, Origin of Cosmic Rays, Progress of Theoretical Physics Supplement 6 (1958) 1–92.  arXiv:https://academic.oup.com/ptps/article-pdf/doi/10.1143/PTPS.6.1/5269796/6-1.pdf, doi:10.1143/PTPS.6.1.
URL https://doi.org/10.1143/PTPS.6.1

  * [4] M. Aguilar, et al., The alpha magnetic spectrometer (ams) on the international space station: Part i – results from the test flight on the space shuttle, Physics Reports 366 (6) (2002) 331–405.  doi:https://doi.org/10.1016/S0370-1573(02)00013-3.
URL https://www.sciencedirect.com/science/article/pii/S0370157302000133

  * [5] O. Adriani, et al., Measurements of cosmic-ray hydrogen and helium isotopes with the PAMELA experiment 818 (1) (2016) 68.  doi:10.3847/0004-637x/818/1/68.
URL https://doi.org/10.3847/0004-637x/818/1/68

  * [6] G. A. de Nolfo, L. M. Barbier, E. R. Christian, A. J. Davis, R. L. Golden, M. Hof, K. E. Krombel, A. W. Labrador, W. Menn, R. A. Mewaldt, J. W. Mitchell, J. F. Ormes, I. L. Rasmussen, O. Reimer, S. M. Schindler, M. Simon, S. J. Stochaj, R. E. Streitmatter, W. R. Webber, A measurement of cosmic ray deuterium from 0.5–2.9 gev/nucleon, AIP Conference Proceedings 528 (1) (2000) 425–428.  arXiv:https://aip.scitation.org/doi/pdf/10.1063/1.1324352, doi:10.1063/1.1324352.
URL https://aip.scitation.org/doi/abs/10.1063/1.1324352

  * [7] M. Aguilar, et al., Properties of cosmic helium isotopes measured by the alpha magnetic spectrometer, Phys. Rev. Lett. 123 (2019) 181102.  doi:10.1103/PhysRevLett.123.181102.
URL https://link.aps.org/doi/10.1103/PhysRevLett.123.181102

  * [8] E. Bueno, F. Barão, M. Vecchi, A parametric approach for the identification of single-charged isotopes with ams-02, Nuclear Instruments and Methods in Physics Research Section A: Accelerators, Spectrometers, Detectors and Associated Equipment 1031 (2022) 166564\.  doi:https://doi.org/10.1016/j.nima.2022.166564.
URL https://www.sciencedirect.com/science/article/pii/S016890022200167X

  * [9] G. D’Agostini, Improved iterative Bayesian unfolding, in: Alliance Workshop on Unfolding and Data Correction, 2010.  arXiv:1010.0632.
  * [10] H. Mülthei, B. Schorr, On an iterative method for the unfolding of spectra, Nuclear Instruments and Methods in Physics Research Section A: Accelerators, Spectrometers, Detectors and Associated Equipment 257 (2) (1987) 371–377.  doi:https://doi.org/10.1016/0168-9002(87)90759-5.
URL https://www.sciencedirect.com/science/article/pii/0168900287907595

  * [11] M. Aguilar, et al., The alpha magnetic spectrometer (ams) on the international space station: Part ii — results from the first seven years, Physics Reports 894 (2021) 1–116, the Alpha Magnetic Spectrometer (AMS) on the International Space Station: Part II - Results from the First Seven Years.  doi:https://doi.org/10.1016/j.physrep.2020.09.003.
URL https://www.sciencedirect.com/science/article/pii/S0370157320303434

  * [12] V. Bindi, et al., Calibration and performance of the ams-02 time of flight detector in space, Nuclear Instruments and Methods in Physics Research Section A: Accelerators, Spectrometers, Detectors and Associated Equipment 743 (2014) 22–29.  doi:https://doi.org/10.1016/j.nima.2014.01.002.
URL https://www.sciencedirect.com/science/article/pii/S0168900214000102

  * [13] F. Giovacchini, J. Casaus, A. Oliva, The AMS-02 RICH detector: Status and physics results, Nuclear Instruments and Methods in Physics Research Section A: Accelerators, Spectrometers, Detectors and Associated Equipment 952 (2020) 161797.  doi:10.1016/j.nima.2019.01.024.
URL https://doi.org/10.1016/j.nima.2019.01.024

  * [14] J. Allison, et al., Recent developments in geant4, Nuclear Instruments and Methods in Physics Research Section A: Accelerators, Spectrometers, Detectors and Associated Equipment 835 (2016) 186–225.  doi:https://doi.org/10.1016/j.nima.2016.06.125.
URL https://www.sciencedirect.com/science/article/pii/S0168900216306957

  * [15] M. Aguilar, et al., Precision measurement of the proton flux in primary cosmic rays from rigidity 1 gv to 1.8 tv with the alpha magnetic spectrometer on the international space station, Phys. Rev. Lett. 114 (2015) 171103.  doi:10.1103/PhysRevLett.114.171103.
URL https://link.aps.org/doi/10.1103/PhysRevLett.114.171103

  * [16] Boudaud, M., Bueno, E. F., Caroff, S., Genolini, Y., Poulin, V., Poireau, V., Putze, A., Rosier, S., Salati, P., Vecchi, M., The pinching method for galactic cosmic ray positrons: Implications in the light of precision measurements, A&A 605 (2017) A17.  doi:10.1051/0004-6361/201630321.
URL https://doi.org/10.1051/0004-6361/201630321

  * [17] K. Cranmer, Kernel estimation in high-energy physics, Computer Physics Communications 136 (3) (2001) 198–207.  doi:https://doi.org/10.1016/S0010-4655(00)00243-5.
URL https://www.sciencedirect.com/science/article/pii/S0010465500002435


# PhilEO Bench: Evaluating Geo-Spatial Foundation Models

###### Abstract

Massive amounts of unlabelled data are captured by Earth Observation (EO)
satellites, with the Sentinel-2 constellation generating 1.6TB of data daily.
This makes Remote Sensing a data-rich domain well suited to Machine Learning
(ML) solutions. However, a bottleneck in applying ML models to EO is the lack
of annotated data as annotation is a labour-intensive and costly process. As a
result, research in this domain has focused on Self-Supervised Learning and
Foundation Model approaches. This paper addresses the need to evaluate
different Foundation Models on a fair and uniform benchmark by introducing the
PhilEO Bench, a novel evaluation framework for EO Foundation Models. The
framework comprises of a testbed and a novel 400GB Sentinel-2 dataset
containing labels for three downstream tasks, building density estimation,
road segmentation, and land cover classification. We present experiments using
our framework evaluating different Foundation Models, including Prithvi and
SatMAE, at multiple $n$-shots and convergence rates.

Index Terms—  Foundation Models, Earth Observation

##  1 Introduction

Problem. Label-efficient learning is challenging as the performance of deep
learning models depends highly on the size and the quality of the training
data. However, for several real-world applications, labelling large datasets
is laborious, expensive, and time-consuming. This holds for Earth Observation
(EO) and remote sensing [1], where the data usually require domain expertise:
the biologist and the cartographer see different things. Also, the aspect of
the Earth changes over time because of nature (seasonality), man-induced
changes, and natural hazards, including volcanic eruptions and floods. For a
specific geographical region, this can induce a change of label or not. Data
avalanche in EO is another challenge, as many satellites orbit the Earth, such
as the Sentinel-2 constellation [2], and produce a lot of data. Labelling them
would require massive human interactions and frequent updates.

For real-world applications [3, 4], including land cover mapping, training
deep learning models with less labeled data is crucial. This can be achieved
using Self-Supervised Learning (SSL) [1, 5], that is learning supervised on a
pretext task attached to the structure of the data itself.

Fig. 1: The PhilEO Suite, where complex neural networks are trained on massive
data. On the right, the PhilEO Bench evaluates such foundation models on
diverse downstream tasks. Table 1:  Current Foundation Models, trained by
self-supervised learning and evaluated on disparate downstream tasks.

Name |  Dataset |  Dataset features |  Architecture |  Downstream tasks |  Benchmark
---|---|---|---|---|---
Prithvi |  HLS V2 L30 |  Coverage: USA, 30m res, 6 bands, Size: 1TB, Time |  ViT, MAE, 3D, 100M params, Geolocation |  Flood detection, Burn scars, Multi-temp crop |  Sen1Floods11: mIoU 88.7%. HLS Burn Scars: IoU 73% burn class
SatMAE |  fMoW |  Size: 3.5 TB, Global, 10m res, 8 bands |  ViT, MAE, With geolocation, With time |  LULC, Multi-label, Building segmentation |  NAIP: Acc: 72%, BEN: mAP: 82%, SpaceNet: mIoU: 78%
SeCo |  SeCo |  Size: 36 GB, 1M samples, Global, 10m res, RGB |  ResNet, MoCo-v2, 23M params, Geolocation |  LULC, Change detection, Fine-tuning and LP |  BEN: MAP 87%. EuroSAT: Acc: 93%. OSCD: F1: 46.9%
Satlas |  SatlasPre- train, 137 classes |  Size: 30 TB, Global, 10m res & 1m, 3 bands (RGB) and 9, With time |  Swin Transf, Shifted window self-attention, Labeled data, Geolocation |  LULC, Segmentation of roads, buildings, ships. Multi-scale features |  UCM: Acc: 99%, RESISC: 98%, AID: 88%, FMoW: F1-score: 44%, Roads: 87%, Buildings: 88%
|  |  |  |  PhilEO (Ours): Building density estimation, Road extraction, LULC |  Evaluation framework using annotated Sentinel-2 data that are global, 10m res, 10 bands

HLS = Harmonised Landsat Sentinel-2; LULC = Land Use Land Cover; LP = Linear
Probing; BEN = BigEarthNet dataset, OSCD = Onera Satellite Change Detection

Our contributions. The main aim of Foundation Models (FM) is to develop
effective SSL methods to find good feature representations [5, 6]. Thus, they
could be easily adapted to a variety of downstream tasks with good
performances, even without much labels in the target domain. The PhilEO
Bench111Phil in “PhilEO” is from the Greek philia for love, and EO refers to
Earth Observation (EO). Phi is also a wink to ESA $\Phi$-lab ;). is our
proposed evaluation framework for benchmarking geo-spatial FMs on downstream
tasks. The EO Data, the same for all downstream tasks we examine, are global,
come from Sentinel-2, and include all bands. The downstream tasks (with
corresponding labels) are building density estimation, road segmentation, and
land cover classification. For the sake of fair comparison, the only variable
module is the FM, while the following decoder architecture remains unaltered.
The PhilEO evaluation of a FM is shown in Fig. 1. In the PhilEO Bench, we
evaluate different EO FMs, including the state-of-the-art Prithvi or SatMAE
[7, 8].

##  2 Related work

Self-supervised learning (SSL) [9, 10] is a machine learning paradigm that
uses unlabeled data to train models to recognize patterns and make decisions
that can assist with downstream tasks [1, 11]. This is done by training the
model to solve a task that requires understanding of the data (pretext task)
[5]. The potential of SSL has been verified also in EO, e.g. Seasonal Contrast
(SeCo) [12, 3]. Promising directions for applications have been identified in
[1, 13, 14] and in [15, 16]. Using SSL has shown to be effective in few-shot
learning settings on EO tasks [3, 1]. Thanks to the emergence of large
datasets [1, 17], such as SSL4EO-S12 [5] which contains 500 GB of S2 data,
models have been trained and validated for downstream tasks [18], for example
CORSA [19] or the model USat [20].

Foundation Models for Earth science applications are currently developed [7].
ClimaX [21], which is for weather and climate, can be tailored to various
tasks with different data, coverage, and physical backgrounds. SatMAE [8] is a
Vision Transformer (ViT) based on Masked Auto-Encoder (MAE) [22] using
temporal and multi-spectral information, trained on Very-High Resolution (VHR)
and S2 data at 10m resolution. In Table 1, we list characteristics of existing
EO FMs. Prithvi [7] is a geospatial FM for EO. Trained on US data only, it is
not a global/ planetary-scale model and uses six bands at 30m resolution.
Prithvi and Scale-MAE [23] training is based on SatMAE. PRESTO [24] is a pre-
trained Transformer focusing on time-series data. Similarly, CaCo [25]
improves upon SeCo [3] by removing seasonality-invariant training features
unsuitable for change detection tasks.

Geo-Bench [26] is a platform for evaluating FMs (with various backbones or
training) on downstream tasks. Different from [26] that compiles already
existing datasets, the PhilEO Bench introduces a new global dataset apt to
test universality of FMs.

##  3 The PhilEO Evaluation Framework

The aim of Foundation Models is to improve the performance on downstream tasks
[27, 28, 29], while also being generalizable and flexible [30]. However, these
models are often evaluated on a range of datasets with different
characteristics (size, resolution, locations, satellite sources, and capture
dates). There is also a focus on evaluating classification downstream tasks,
while omitting image-to-image (such as segmentation) downstream tasks.
Therefore, it is challenging to fairly compare the performance of these
burgeoning FMs and draw meaningful conclusions. To evaluate FMs, we propose
the PhilEO Bench, an evaluation framework with the aim of providing a
flexible, consistent, and fair benchmark for EO Sentinel-2 FMs. Our main
contributions are: (1) A novel global dataset of S2 images and downstream task
labels designed to assess the generalizablity of models across different
geographic regions and tasks. (2) A new flexible evaluation framework focused
on generating comparable, fair, and reproducible results.

###  3.1 PhilEO Downstream Dataset

The PhilEO dataset is a 400 GB global dataset of S2 images and has labels for
roads, buildings, and land cover, where these are the three downstream tasks.
The data is sampled from geographically diverse regions around the globe
including: Denmark, East Africa, Egypt, Guinea, Europe, Ghana, Israel, Japan,
Nigeria, North America, Senegal, South America, Tanzania, and Uganda. Each
region has up to 200 tiles of varying sizes. Some locations have been
revisited up to $3$ times. The data contain $11$ bands at 10m resolution in
the following order: 0-SCL, 1-B02, 2-B03, 3-B04, 4-B08, 5-B05, 6-B06, 7-B07,
8-B8A, 9-B11, and 10-B12 where SCL is the Scene Classification Layer222For
more information about the Sentinel-2 spectral bands, please see:
http://docs.sentinel-hub.com/api/latest/data/sentinel-2-l2a. As shown in Fig.
2, each S2 tile in PhilEO has a label for each of the downstream tasks:

  * •

ROADS: The labels are expressed as a number of squared meters of roads in a
given pixel. The values are between $0$ and $100$, and for a resolution of
10m, this reflect the percentage of coverage.

  * •

BUILDINGS: The labels are expressed as squared meters of buildings. The values
are between $0$ and $100$. For 10m resolution, this reflect the coverage (in
%).

  * •

LANDCOVER: Land cover labels are taken from ESA World
Cover333http://worldcover2020.esa.int/data/docs/WorldCover$\\_$PUM$\\_$V1.1.pdf:
11 classes, e.g. tree cover and built-up.

All the three downstream tasks can be treated both as segmentation and as
classification tasks.

###  3.2 Evaluation framework

One of the biggest challenges in evaluating FMs is disentangling the
performance impact of various factors such as: model architectures, pre-
training tasks, and downstream task training data. The effectiveness of a FM
can be measured by the quality of its latent representations, and how the key
features learnt through the process of pre-training can boost downstream task
performance [1, 31]. Hence, to provide a fair comparison between different FMs
within our evaluation framework, we minimize the impact of confounding
variables by providing: (1) consistent and repeatable training and evaluation
datasets, and (2) a common downstream task head for all the pre-trained
models.

  * •

Dataset creation: To minimize the impact of variability in the downstream task
datasets, a common hold out Test Set for each downstream task was created
providing comparable results across all evaluated models. A data partitioning
script is also provided that allows for the creation of smaller subsets
($n$-shot or percentage split) from the full PhilEO Downstream Task training
dataset enabling us to evaluate the impact of training set size on model
performance. The indexes of the training samples used are automatically saved.
Hence, different models can be trained on the identical sub-datasets.

  * •

One head to rule them all: To minimize the impact of different decoder heads
on the downstream task performance, we propose the use of a common decoder
head. The latent representations of each of the FMs [7, 32] are fed to a
similar decoder. Hence, the comparable performance of a model is a consequence
of the effectiveness of the pre-training task and the representational
strength of its latent space. For segmentation downstream tasks, we use a
multi-convolution decoder based on the U-Net design. For some models, it is
required to up-sample or down-sample the features to achieve the desired
output image size. For classification downstream tasks, a linear decoder is
used 444For ViT models, the cls token is used as the input to the linear
decoder..

Fig. 2: From left: 1) RGB visualisation of S2 patch, 2) Land cover labels, 3)
Building Density labels, and 4) Road Density.

Our framework supports two training configurations: (1) Fine-tuning, which
allows for updating of all downstream task model weights including the FM
encoder, and (2) Linear probing, where only the decoder head weights are
updated, freezing the FM encoder parameters. PhilEO also contains U-Net [33],
Mixer [34], and ViT [35] architectures. PhilEO supports pre-trained models
such as Masked Auto-Encoder (MAE) [22] ViT, and Pre-trained U-Nets [33, 36],
as well as the models Prithvi [7], SatMAE [8], and SeCo [3]. In addition, the
testbed should be flexible and easy to use. Hence, an Object Oriented
Programming approach is used with an emphasis on modularity, allowing for the
easy addition of other downstream tasks, architectures, and pre-trained
models.

##  4 Evaluating Foundation Models

We use the framework of Sec. 3.2 to evaluate the performances of our own pre-
trained models, including MAE ViTs555Based on the SatMAE implementation and
Pre-trained U-Nets. We also evaluate the state-of-the-art Prithvi, SatMAE, and
SeCo. In the next sections, we present the experiment, implementation details,
and the results.

###  4.1 Experiment for the three downstream tasks

A $n$-shot transfer learning experiment for the three downstream tasks
presented in Sec. 3.1 is set up. The experiment aims to examine the impact of
different downstream tasks and training set sizes on model performance. We
also compare the linear probing and fine-tuning transfer learning paradigms.

###  4.2 Implementation details for the experiment

All models are trained until convergence with an initial learning rate of
$0.0001$666Using a AdamW optimizer and reduce on plateau and a batch size of
$64$. The in-house models we developed are trained on 10-band 10m resolution
128-by-128 patches. The SatMAE models are trained on 10-band 10m resolution
96-by-96 patches. Prithvi is trained on 8-band 30m resolution 224-by-224
patches, while SeCo and ResNet are trained on RGB 10m resolution 224-by-224
patches.

Fig. 3: Land cover classification accuracy (acc) using PhilEO Bench for $n$
samples per region with different architectures and transfer learning
paradigms: (1) Linear probing (lp), and (2) Fine-tuning (ft). The legend is
for both this plot and Fig. 4.

###  4.3 Results and discussion

From the PhilEO Bench, we pick out two downstream tasks and we report the
results we obtain:

Image-to-Image Downstream Tasks:

  * •

Land Cover: In Fig. 3, we compare the land cover segmentation accuracy of
various benchmark and pre-trained models versus the training dataset size. The
pre-trained models outperform their randomly initialized fully supervised
counterparts, with the pre-trained (fine-tuned) Geo-Aware U-Net consistently
giving the best overall performance. The ViTCNN models seem to underperform
when compared to the basic U-Net.

  * •

Building Density: In Fig. 4, we compare the building density estimation MSE of
the benchmark and pre-trained models versus the training dataset size. We
observe that the Prithvi-based models result in the lowest MSE scores.
However, on visual inspection of the results, it is clear that the outputs
generated by the U-Net models are better at picking out finer details than the
Prithvi models, as shown in Fig. 5. The discrepancy between the visual results
and the MSE scores is because the Prithvi models are trained and evaluated on
30m resolution data, skewing the results in their favour. Fairly evaluating
models across different resolutions is a point of future work for the
evaluation framework.

The results show that state of the art Sentinel-2 FMs can be outperformed by a
simple U-Net architecture when applied to image-to-image downstream tasks and
$n$-shot transfer learning experiments777Our code can be found in the GitHub
repository http://github.com/ESA-PhiLab/PhilEO-Bench. These results are not
all that surprising when one considers the encoder-decoder structure of these
pre-trained models. The amount of information that can be stored in the latent
space (bottleneck) is limited. Therefore, only salient high-level features are
compressed into the latent space. The lack of low-level information in the
latent space makes it challenging to reconstruct fine-grained details in the
output image, such as building locations, resulting in poor image-to-image
downstream task performance. The U-Net skip connections address this by
bypassing the bottleneck and providing low-level features to the decoder,
allowing for the reconstruction of fine-grained details in the output image.

Fig. 4: MSE for building density regression using the PhilEO Bench for $n$
samples per region, evaluating different models.

Fig. 5: From left: (1) RGB visualisation of Sentinel-2 patch, (2) Building
Density labels, and (3) Model predictions. Top row is a U-Net model, while
bottom row is the Prithvi model.

##  5 Conclusion

We have proposed the PhilEO evaluation framework which comprises both a new
testbed for benchmarking Sentinel-2 FMs and a novel global dataset. With this
evaluation framework, we hope to promote the advancement of Sentinel-2 FMs by
providing a general, fair, and easy to use evaluation tool. The experiments we
carried out on the PhilEO Bench have exposed a weakness in the current state
of the art EO FMs. By mainly focusing on classification tasks, most existing
models are based on the Masked Auto-Encoder (MAE) Vision Transformer (ViT)
architecture. However, this model architecture does not lend itself well to
image-to-image downstream tasks, with even a basic U-Net outperforming models
such as SatMAE and Prithvi. We believe that this case study underlines the
utility of our evaluation testbed and the importance of having a fair
comparison between different models.

##  6 Acknowledgement

We thank for fruitful discussions Leonardo Bagaglini, Giacomo D. Cascarano and
Giorgio Pasquali (Leonardo Labs).

## References

  * [1] Yi Wang et al.,  “Self-supervised learning in remote sensing: A review,”  IEEE Geosc & Rem Sen Mag, 2022.
  * [2] European Space Agency (ESA),  “Sentinel-2 Operations,”  Enabling and Support, Operations, http://www.esa.int/Enabling_Support/Operations/Sentinel-2_operations.
  * [3] Oscar Manas et al.,  “Seasonal contrast: Unsupervised pre-training from uncurated remote sensing data,”  in Proc. ICCV, 2021.
  * [4] Kumar Ayush, Burak Uzkent, et al.,  “Geography-aware self-supervised learning,” arxiv:2011.09980, 2021.
  * [5] Yi Wang et al.,  “SSL4EO-S12: A large-scale multi-modal, multi-temporal dataset for self-supervised learning in Earth Observation,” arxiv:2211.07044, 2022.
  * [6] Philipe Dias, Abhishek Potnis, Sreelekha Guggilam, et al.,  “An agenda for multimodal Foundation Models for Earth Observation,” In Proc. IGARSS, 2023.
  * [7] Johannes Jakubik, Sujit Roy, C. E. Phillips, Paolo Fraccaro, et al.,  “Foundation Models for Generalist Geospatial Artificial Intelligence,” arxiv:2310.18660, 2023.
  * [8] Yezhen Cong, Samar Khanna, et al.,  “SatMAE: Pre-training Transformers for temporal and multi-spectral satellite imagery,” arxiv:2207.08051, 2023.
  * [9] A. Dosovitskiy et al.,  “Discriminative unsupervised feature learning with exemplar CNNs,”  Trans. PAMI, 2016.
  * [10] Carl Doersch and Andrew Zisserman,  “Multi-task self-supervised visual learning,”  in Proc. ICCV, 2017.
  * [11] Michail Tarasiou and Stefanos Zafeiriou,  “Embedding Earth: Self-supervised contrastive pre-training for dense land cover classification,” arxiv:2203.06041, 2022.
  * [12] Yang et al. Yu,  “An unsupervised convolutional feature fusion network for deep representation of remote sensing images,”  IEEE Geoscience Rem Sens Letters, 2018.
  * [13] Paul Berg et al.,  “Self-supervised learning for scene classification in remote sensing: Current state of the art and perspectives,”  Remote Sensing, 2022.
  * [14] D. Tuia et al.,  “Artificial Intelligence to advance Earth observation: A perspective,” arxiv:2305.08413, 2023.
  * [15] Favyen Bastani et al.,  “Satlas: A large-scale, multi-task dataset for RS image understanding,” In ICCV, 2023.
  * [16] Matias Mendieta et al.,  “Towards geospatial foundation models via continual pretraining,” In Proc. ICCV, 2023.
  * [17] Qian Shi et al.,  “Globe230k: A benchmark dense-pixel annotation dataset for global land cover mapping,”  J Rem Sen AAAS, 2023.
  * [18] Gencer Sumbul et al.,  “BigEarthNet-MM: A large scale multi-modal multi-label benchmark archive for RS image classification,”  IEEE Geosc Rem Sens Mag, 2021.
  * [19] Xenia Ivashkovych et al.,  “CORSA deep Earth Observation semantic compression applied to flood detection,”  ESA OBPDC workshop, 2022.
  * [20] J. Irvin et al.,  “USat: A unified self-supervised encoder for multi-sensor satellite,”  arXiv:2312.02199, 2023.
  * [21] Tung Nguyen et al.,  “ClimaX: A foundation model for weather and climate,” arxiv:2301.10343, 2023.
  * [22] Kaiming He et al.,  “Masked autoencoders are scalable vision learners,” In Proc. CVPR, 2022.
  * [23] Colorado J. Reed, Ritwik Gupta, et al.,  “Scale-MAE: A scale-aware masked autoencoder for multiscale geospatial representation learning,” arxiv:2212.14532, 2022.
  * [24] G. Tseng et al.,  “Lightweight, pre-trained Transformers for remote sensing timeseries,” arxiv:2304.14065, 2023.
  * [25] Utkarsh Mall et al.,  “Change-aware sampling and contrastive learning for satellite images,”  in Proc. CVPR, 2023.
  * [26] A. Lacoste et al.,  “Geo-Bench: Toward Foundation Models for Earth monitoring,” arxiv:2306.03831, 2023.
  * [27] Wenyuan Li et al.,  “Geographical knowledge-driven representation learning for remote sensing images,”  IEEE Trans Geoscience and Rem Sens, 60:1-16, 2022.
  * [28] Zhenmei Shi, Jiefeng Chen, et al.,  “The trade-off between universality and label efficiency of representations from contrastive learning,”  in Proc. ICLR, 2023.
  * [29] K. Klemmer et al.,  “SatCLIP: Global, general-purpose location embeddings,”  arXiv:2311.17179, 2023.
  * [30] Isaac Corley et al.,  “Revisiting pre-trained remote sensing model benchmarks: Resizing and normalization matters,” arXiv:2305.13456, 2023.
  * [31] Danfeng Hong et al.,  “SpectralGPT: Spectral foundation model,”  arXiv:2311.07113, 2023.
  * [32] Randall Balestriero et al.,  “A cookbook of self-supervised learning,”  arXiv:2304.12210, 2023.
  * [33] O. Ronneberger et al.,  “U-Net: Convolutional networks for biomedical image segmentation,”  In MICCAI, 2015.
  * [34] Ilya O Tolstikhin et al.,  “MLP-Mixer: An all-MLP architecture for vision,”  In Proc. NeurIPS, 2021.
  * [35] A. Dosovitskiy et al.,  “An image is worth 16x16 words: Transformers for recognition at scale,”  In ICLR, 2021.
  * [36] M. Kottek et al.,  “World map of the Köppen-Geiger climate classification,”  Meteor Z, v. 15, p. 259-263, 2006.