# Resonant generation of high-order harmonics in nonlinear electrodynamics

###### Abstract

We study the process of resonant generation of high-order harmonics in a
closed cavity in the model of vacuum nonlinear electrodynamics. Concretely, we
study the possibility of resonant generation of the third harmonic induced by
a single electromagnetic mode in a radiofrequency cavity, as well as resonant
generation of a combined frequency mode induced by two pump modes
($\omega_{1}$ and $\omega_{2}$). We explicitly show that the third harmonic as
well as the $2\omega_{1}+\omega_{2}$ combined frequency mode are not
resonantly amplified, while the $2\omega_{1}-\omega_{2}$ signal mode is
amplified for certain cavity geometry. We discuss the process from the point
of view of quantum theory.

a Lomonosov Moscow State University, Faculty of Physics b Institute for
Nuclear Research of RAS

PACS: 44.25.$+$f; 44.90.$+$c

## Introduction

Vacuum nonlinearity is one of phenomena theoretically predicted at the dawn of
quantum electrodynamics[1, 2] but have not experimentally detected yet due to
its extreme smallness. A part of the effects of vacuum nonlinearity related to
electrodynamics mimic to similar effects in nonlinear optical crystals, which
include vacuum birefringence of a photon in external field, and the generation
of high-order harmonics. The latter can be tested both in optical and radio
wave range. The generation of high-order harmonics for radio modes in cavity
may be potentially detected only in case of resonant increase [3], see partial
solutions in [4, 5]. An interesting issue is to study this nonlinear process
for arbitrary set of cavity mode, which can be done analytically for
rectangular shape of cavity. Thus, studying nonlinear effects for single mode
in one-dimensional cavity it turns out that there is no resonant generation of
the third harmonics [6]. We generalise the method to the case of two pump
modes and 3D rectangular cavity. More details are given in [7]; in addition to
that article we discuss the quantum aspects of the obtained result in the
conclusion.

## General theory for resonance

We start from Euler-Heisenberg effective Lagrangian [1, 2],

$\mathcal{L}=-\frac{1}{4}\mathcal{F}+\kappa\,(\mathcal{F}^{2}+\beta\,\mathcal{G}^{2}),\qquad\kappa=\frac{\alpha_{e}^{2}}{90\,m_{e}^{4}},\quad\beta=\frac{7}{4},$
(1)

The electromagnetic field invariants have the standard form,

$\mathcal{F}\equiv
F_{\mu\nu}F^{\mu\nu}=-2\left(\textbf{E}^{2}-\textbf{B}^{2}\right),\qquad\qquad\mathcal{G}\equiv
F_{\mu\nu}{\widetilde{F}}^{\mu\nu}=-4\left(\textbf{E}\cdot\textbf{B}\right).$
(2)

Varying the Lagrangian, one obtains modified Maxwell equations [1]:

$\begin{aligned} \,\text{rot}\textbf{B}&=\frac{\partial\textbf{E}}{\partial
t}+\left[\frac{\partial\textbf{P}}{\partial
t}-\,\text{rot}\textbf{M}\right],\\\
\,\text{rot}\textbf{E}&=-\frac{\partial\textbf{B}}{\partial t},\\\
\end{aligned}\hskip 56.9055pt\begin{aligned} \,\text{div}\textbf{B}&=0,\\\
\,\text{div}\textbf{E}&=-\,\text{div}\textbf{P},\end{aligned}$ (3)

where P and M denote vacuum polarization and magnetization respectively,

$\displaystyle\textbf{P}(\textbf{x},t)$ $\displaystyle\equiv
16\,\kappa\left[\left(\textbf{E}^{2}-\textbf{B}^{2}\right)\textbf{E}+2\beta\left(\textbf{E}\cdot\textbf{B}\right)\textbf{B}\right],$
(4) $\displaystyle\textbf{M}(\textbf{x},t)$ $\displaystyle\equiv
16\,\kappa\left[\left(\textbf{E}^{2}-\textbf{B}^{2}\right)\textbf{B}-2\beta\left(\textbf{E}\cdot\textbf{B}\right)\textbf{E}\right].$

Wave equations both for amplitudes for electric and magnetic fields are
modified as well [1],

$\displaystyle\Box\textbf{E}$ $\displaystyle=\frac{\partial}{\partial
t}\,\text{rot}\textbf{M}+{\,\text{grad}\,\text{div}\textbf{P}}-\frac{\partial^{2}\textbf{P}}{\partial{t}^{2}},$
(5) $\displaystyle\Box\textbf{B}$ $\displaystyle=\frac{\partial}{\partial
t}\,\text{rot}\textbf{P}-{\,\text{grad}\,\text{div}\textbf{M}}+\Delta\textbf{M}.$

These equations are nonlinear; the solving them seems to be a hard issue.
However, for the our goal we may apply the perturbation theory.

We call the mode, initially given in the cavity, as a ‘‘pump mode’’ (electric
field $\textbf{E}^{p}$), and look for the evolution of ‘‘signal mode’’
(electric field $\textbf{E}^{sig}$) which was not initially present in the
cavity. Assuming the hierarchy
$|\textbf{E}^{sig}|\sim\kappa\left(|\textbf{E}^{p}|\right)^{3}\ll|\textbf{E}^{p}|$,
one obtains in the zeroth order homogeneous wave equations for the pump modes
$\Box\textbf{E}^{p}=0,\ \Box\textbf{B}^{p}=0$, and in the first order
inhomogeneous linear wave equations for signal mode amplitudes:

$\displaystyle\Box\textbf{E}^{sig}$ $\displaystyle=\frac{\partial}{\partial
t}\,\text{rot}\textbf{M}(\textbf{E}^{p},\textbf{B}^{p})+{\,\text{grad}\,\text{div}\textbf{P}(\textbf{E}^{p},\textbf{B}^{p})}-\frac{\partial^{2}\textbf{P}(\textbf{E}^{p},\textbf{B}^{p})}{\partial{t}^{2}},$
(6) $\displaystyle\Box\textbf{B}^{sig}$
$\displaystyle=\frac{\partial}{\partial
t}\,\text{rot}\textbf{P}(\textbf{E}^{p},\textbf{B}^{p})-{\,\text{grad}\,\text{div}\textbf{M}(\textbf{E}^{p},\textbf{B}^{p})}+\Delta\textbf{M}(\textbf{E}^{p},\textbf{B}^{p}).$

These inhomogeneus wave equations may have resonantly growing solutions if the
r.h.s. contains terms which coincide with the solution of corresponding
homogeneous equations. This resonance means the linear growing with time (see
1D component for example),

$\Box E_{x}=\cos(\omega_{r}t)\sin(\omega_{r}x)\,+\,...\ \ \to\ \
E_{x}^{(growing)}=\frac{t}{\omega_{r}}\sin(\omega_{r}t)\sin(\omega_{r}x).$

In a real world this linear growth stops when the dissipation effect become
significant. Introducing dissipation coefficient $\Gamma$ ‘‘by hands’’, on
obtains the saturation of the linear growth:

$(\Box+\Gamma\partial_{t})E_{x}=\cos(\omega_{r}t)\sin(\omega_{r}x)\ \ \to\ \
E_{x}^{(steady)}=\frac{1}{\Gamma\omega_{r}}\sin(\omega_{r}t)\sin(\omega_{r}x),$

the amplitude of resonant mode inverse proportianal to $\Gamma$.

We apply this approach to single and two pump modes in 1D and 3D cavities.

## One-dimensional cavity

We call by ‘‘one-dimensional cavilty’’ a rectangular cavity with one spatial
dimension smaller than others, say $L_{x}\ll L_{y},L_{z}$. A pump mode
configuration reads (only nonzero components),

$E^{p}_{y}(x,t)=F_{0}\cdot\sin(k_{n}x)\sin(\omega_{n}t),\qquad
H^{p}_{z}(x,t)=F_{0}\cdot\cos(k_{n}x)\cos(\omega_{n}t),$ (7)

where the eigenfrequencies and wave vectors take discrete number of values:
$\omega_{n}=k_{n}=n\cdot\frac{\pi}{L_{x}},\ n\in\mathbb{N}$. The corresponding
electromagnetic field configuration is shown at Fig. 1.

Figure 1: Field configuration of an electromagnetic mode in one-dimensional
cavity ($n=1$).

To study the evolution of a signal mode initiated by one pump mode in cavity,
we write the linearized wave equations (6) with r.h.s. calculated at the pump
mode configuration (7)111Analytical calculations were made in “wxMaxima
21.02.0” computer algebra system [8] (cf. [6]),

$\displaystyle\left(\Box+\Gamma\partial_{t}\right)E_{y}^{sig}$
$\displaystyle=8\kappa
F_{0}^{3}\omega_{n}^{2}\Bigl{[}2\sin({\omega_{n}}x)\sin({\omega_{n}}t)-3\sin({\omega_{n}}x)\sin({3\omega_{n}}t)+\sin({3\omega_{n}}x)\sin({\omega_{n}}t)\Bigr{]},$
(8) $\displaystyle\left(\Box+\Gamma\partial_{t}\right)H_{z}^{sig}$
$\displaystyle=8\kappa
F_{0}^{3}\omega_{n}^{2}\Bigl{[}2\cos({\omega_{n}}x)\cos({\omega_{n}}t)-\cos({\omega_{n}}x)\cos({3\omega_{n}}t)+3\cos({3\omega_{n}}x)\cos({\omega_{n}}t)\Bigr{]}.$

The first terms in both lines of the r.h.s. of eqs. (8) represent solution of
homogeneous equation so there is a resonant mode with the pump mode frequency
$\omega_{n}$; the other two terms are non-resonant. Surprisingly, there are no
resonant terms with triple frequency in the r.h.s of eqs. (8) so the third
harmonics of the pump mode does not appear.

For the simplicity of presentation of the following results we present the
r.h.s. as a table, see Table 1.

wavenumbers | $n$ | $3n$
---|---|---
eigenfrequencies | $\omega_{n},~{}\omega_{3n}$ | $\omega_{n}$
Table 1: Examination of the resonance criterion for a single pump mode in
1D-cavity.

We see that the only eigenfrequency connected with wave vector is
$\omega_{n}$.

The next step is to consider two pump modes in one-dimensional cavity.
Generally, we should take into account arbitrary angle $\alpha$ between
polarization plane of two modes with wavenumbers $n$ and $p$. Following the
same algorithm as in the the previous case, we calculate the r.h.s. on the
pump mode combination which is shortly presented as Table 2 (symmetric part
with $n\leftrightarrow p$ omitted).

wavenumbers | ${n}$ | $3n$ | ${2n-p}$ | ${2n+p}$
---|---|---|---|---
eigenfrequencies | ${\omega_{n}},~{}\omega_{2p\pm n},~{}\omega_{3n}$ | $\omega_{n}$ | $\omega_{p},~{}{\omega_{2n+p}}$ | $\omega_{p},~{}{\omega_{2n-p}}$
Table 2: Examination of the resonance criterion for two pump modes in
1D-cavity.

We see that the only resonant modes have wavenumber $n$; the mixed signal
modes with wavenumbers $2n\pm p$ do not resonate.

## 3D rectangular cavity.

The next step is to consider one or two modes in 3D rectangular cavity of
dimensions ($L_{x},L_{y},L_{z}$). The cavity modes are classified as
$TE_{npq}$ and $TM_{npq}$ modes, see [9], where $n,p,q$ numerate wavenumbers
for each spatial dimension. The wavevector for each mode is determined as
$\textbf{k}_{npq}=\left(\frac{\pi n}{L_{x}},\,\frac{\pi p}{L_{y}},\,\frac{\pi
q}{L_{z}}\right)$, the same for $TE$ and $TM$ modes; the corresponding
frequency is

$\omega_{npq}=|k_{npq}|=\pi\sqrt{\frac{n^{2}}{L_{x}^{2}}+\frac{p^{2}}{L_{y}^{2}}+\frac{q^{2}}{L_{z}^{2}}}.$
(9)

First, consider single mode in the cavity, and apply aforementioned algorithm
in WxMaxima for resonant signal mode searching. The result is shown in Table 3
.

wavenumbers | $n,p,q$ | $3n,p,q$ | $3n,3p,q$ | $\left(n\leftrightarrow p\leftrightarrow q\right)$ | $3n,3p,3q$
---|---|---|---|---|---
eigenfrequencies | $\omega_{npq},~{}3\omega_{npq}$ | $\omega_{npq}$
Table 3: Examination of the resonance criterion for a single pump mode in
3D-cavity.

The only resonant mode has the same frequency as the pump mode; the third
harmonics does not resonantly amplified.

Continue with the two pump mode configuration. This is the most difficult case
due to complicated analytical calculations even using computer algebra system.
First, due to trigonometric relations the r.h.s. of eqs. (8) contain terms
only of the form:
$A\,h(\omega_{sig}t)\,h(k_{sig,x}x)\,h(k_{sig,y}y)\,h(k_{sig,z}z),$ where the
function $h()$ means sine or cosine, and $\omega_{sig}$ and $k_{sig}$ take the
values of the following sets,

$\begin{matrix}\hskip
5.69054pt\omega_{sig}~{}\in~{}\\{&\omega_{1},\quad&\omega_{2},\quad&2\omega_{1}\pm\omega_{2},\quad&2\omega_{2}\pm\omega_{1},\quad&3\omega_{1},\quad&3\omega_{2}\quad\hskip
2.84526pt\\},\\\ k_{sig,x}~{}\in~{}\\{&k_{1x},\quad&k_{2x},\quad&2k_{1x}\pm
k_{2x},\quad&2k_{2x}\pm k_{1x},\quad&3k_{1x},\quad&3k_{2x}\quad\\},\\\
k_{sig,y}~{}\in~{}\\{&k_{1y},\quad&k_{2y},\quad&2k_{1y}\pm
k_{2y},\quad&2k_{2y}\pm k_{1y},\quad&3k_{1y},\quad&3k_{2y}\quad\\},\\\
k_{sig,z}~{}\in~{}\\{&k_{1z},\quad&k_{2z},\quad&2k_{1z}\pm
k_{2z},\quad&2k_{2z}\pm k_{1z},\quad&3k_{1z},\quad&3k_{2z}\quad\\}.\\\
\end{matrix}$ (10)

The r.h.s. of the linearized wave equations (6) generally contain terms from
different columns of (10), thus the mixed terms like
$\sin(3\omega_{1}t)\cos[(2k_{1x}-k_{2x})x]\cos(k_{1y}y)\sin(3k_{2z}z)$ might
hypothetically appear. However, the corresponding amplitude may vanish, as for
the third harmonics. As previously, we calculate the r.h.s. in WxMaxima
system, and obtain the result which can be presented in Table 4. The terms of
the Table 4 can be grouped into two sectors: the triple wavenumbers are not
associated with combined frequencies $2\omega_{1}\pm\omega 2$ (means the
corresponding amplitude vanishes) while the combined wavenumbers are not
associated with triple frequencies. In addition, it turns out that the
amplitude vanishes for parallel wavevectors,
$\textbf{k}_{1}\parallel\textbf{k}_{2}$.

wavenumbers | $\begin{array}[]{c}n_{1}\\\ p_{1}\\\ q_{1}\end{array}$ | $\begin{array}[]{c}3n_{1}\\\ p_{1}\\\ q_{1}\end{array}$ | $\begin{array}[]{c}n_{1}\\\ 3p_{1}\\\ q_{1}\end{array}$ | $\begin{array}[]{c}\\\ \cdots\\\ \end{array}$ | $\begin{array}[]{c}3n_{1}\\\ 3p_{1}\\\ 3q_{1}\end{array}$
---|---|---|---|---|---
eigenfrequencies | $\omega_{1},~{}3\omega_{1},~{}2\omega_{2}+\omega_{1},~{}2\omega_{2}-\omega_{1}$ | $\omega_{1},~{}3\omega_{1}$ | $\omega_{1}$

wavenumbers | $\begin{array}[]{c}2n_{2}\pm n_{1}\\\ p_{1}\\\ q_{1}\end{array}$ | $\begin{array}[]{c}n_{1}\\\ 2p_{2}\pm p_{1}\\\ q_{1}\end{array}$ | $\begin{array}[]{c}n_{1}\\\ p_{1}\\\ 2q_{2}\pm q_{1}\end{array}$ | $\begin{array}[]{c}\cdots\end{array}$ | $\begin{array}[]{c}2n_{2}\pm n_{1}\\\ 2p_{2}\pm p_{1}\\\ 2q_{2}\pm q_{1}\end{array}$
---|---|---|---|---|---
eigenfrequencies | $\omega_{1},~{}2\omega_{2}+\omega_{1},~{}2\omega_{2}-\omega_{1}$

Table 4: Examination of the resonance criterion for two pump modes in
3D-cavity.

Let us prove that the resonant generation does not appear for signal modes
with frequency $2\omega_{2}+\omega_{1}$. Write the dispersion relation for
this signal mode and apply the triangle inequality,

$2\omega_{2}+\omega_{1}=|2\textbf{k}_{2}|+|\textbf{k}_{1}|\geqslant|2\textbf{k}_{2}+\textbf{k}_{1}|=\sqrt{\sum_{i=1}^{3}(2k_{2,i}+k_{1,i})^{2}}.$

In case of non-parallel wavevectors the triangle inequality holds if at least
for one $i$ $k_{sig,i}=3\times\max\left(k_{1,i},~{}k_{2,i}\right)$ —
explicitly the case for which the amplitude vanishes. Thus, we show that the
signal mode with frequency $2\omega_{2}+\omega_{1}$ does not appear.

## Resonant solution for $2\omega_{1}-\omega_{2}$.

In this section we show explicitly the resonant solution describing the
generation of the combined mode with frequency $2\omega_{1}-\omega_{2}$.
Consider modes $TE011$ and $TM110$ for concreteness. Table 4 shows two
possible options for the signal mode:
$2\omega_{011}-\omega_{110}=\omega_{130}$ and
$2\omega_{011}-\omega_{110}=\omega_{132}$. Checking additionally the
dispersion relation (9) for each mode, we came that the second option cannot
be realised for any cavity dimensions, but the first indeed can. The following
condition on the ratio of cavity dimensions reads,

$\left(\frac{L_{z}}{L_{x}}\right)^{2}\left(\frac{L_{z}}{L_{y}}\right)^{2}+\left(\frac{L_{z}}{L_{x}}\right)^{2}+3\left(\frac{L_{z}}{L_{y}}\right)^{2}=1.$

Assuming additionally $L_{x}=L_{y}$, one obtains the resonant ratio of cavity
dimensions $L_{x}:L_{y}:L_{z}=1:1:r,\quad r=\sqrt{\sqrt{5}-2}\approx 0.486$.
Resonant conditions for other set of modes are expected to be computed
similarly.

The corresponding signal amplitude, even being resonantly amplified, is rather
small but detectable by single-photon detectors, see [5, 7] for details.

## Summary and Discussion.

We summarize the conclusions as follows. First, we have shown explicitly that
the third harmonic is not resonantly amplified both in 1D and 3D rectangular
cavities for arbitrary pump mode. Combined ‘‘plus’’ harmonic (frequency
$2\omega_{1}+\omega_{2}$) is not resonantly amplified in 1D and 3D rectangular
cavities for arbitrary set of pump modes. From the other hand, the ‘‘minus’’
mode ($2\omega_{1}-\omega_{2}$) is resonantly amplified in 3D cavity of
certain resonant ratio of dimensions.

Let us have one more glance on the process
$2\omega_{011}-\omega_{110}\to\omega_{130}$. Note that the wavevector
components combinate independently of each over with $\pm$ signs, the
wavenumbers $(130)$ of the signal mode results as follows, $(1=0+0+1)$,
$(3=1+1+1)$, $(0=1-1+0)$.

From the point of view of quantum theory, the generation of the signal mode is
associated with the process $3\to 1$, the merging of three quanta of cavity
modes into a single one. At first glance it seems to be a contradiction: the
only possible result of such process is the energy of final quanta
$2\omega_{1}+\omega_{2}$ due to the energy conservation, while the classical
approach lead to the energy $2\omega_{1}-\omega_{2}$.

This apparent contradiction can be explained as follows. The classical waves
(and cavity modes) are the coherent states. Each coherent state decomposes to
the linear combination of the states with definite number of particles;
nonzero contribution to the process are given from the amplitudes like

$|N_{\omega_{1}},M_{\omega_{2}}\rangle\to|(N-2)_{\omega_{1}}(M+1)_{\omega_{2}}1_{2\omega_{1}-\omega_{2}}\rangle,$

for arbitrary values $M$ and $N$. The calculation of this amplitudes in a pure
quantum approach will be presented in a following paper.

### Acknowledgments

The Authors thank Maxim Fitkevich, Dmitry Kirpichnikov, Dmitry Levkov, Valery
Rubakov, Alexey Rubtsov and Dmitry Salnikov for helpful discussions. The work
is supported by RSF grant 21-72-10151.

## References

  * [1] _Euler H., Kockel B._ The scattering of light by light in Dirac’s theory // Naturwiss. —- 1935\. —- V. 23, no. 15. —- P. 246–247.
  * [2] _Heisenberg W., Euler H._ Consequences of Dirac’s theory of positrons // Z. Phys. —- 1936\. —- V. 98, no. 11-12. —- P. 714–732. —- arXiv:physics/0605038.
  * [3] _Brodin G., Marklund M., Stenflo L._ Proposal for Detection of QED Vacuum Nonlinearities in Maxwell’s Equations by the Use of Waveguides // Phys. Rev. Lett. —- 2001\. —- V. 87. —- P. 171801. —- arXiv:physics/0108022.
  * [4] _Eriksson D., Brodin G., Marklund M., Stenflo L._ A Possibility to measure elastic photon-photon scattering in vacuum // Phys. Rev. A. —- 2004\. —- V. 70. —- P. 013808. —- arXiv:physics/0411054.
  * [5] _Bogorad Z., Hook A., Kahn Y., Soreq Y._ Probing Axionlike Particles and the Axiverse with Superconducting Radio-Frequency Cavities // Phys. Rev. Lett. —- 2019\. —- V. 123, no. 2. —- P. 021801. —- arXiv:1902.01418.
  * [6] _Shibata K._ Intrinsic resonant enhancement of light by nonlinear vacuum // Eur. Phys. J. D. —- 2020\. —- V. 74, no. 10. —- P. 215.
  * [7] _Kopchinskii I., Satunin P._ Resonant generation of electromagnetic modes in nonlinear electrodynamics: Classical approach // Phys. Rev. A. —- 2022\. —- 01\. —- V. 105. —- P. 013508. —- URL: https://link.aps.org/doi/10.1103/PhysRevA.105.013508.
  * [8] https://github.com/Ilia-Ko/Supplemental-Materials/tree/main/Nonlinear-ED/Part-I.
  * [9] _Hill D._ Electromagnetic Fields in Cavities: Deterministic and Statistical Theories // Antennas and Propagation Magazine, IEEE. —- 2014\. —- 02\. —- V. 56. —- P. 306–306.


# SEN12MS-CR-TS: A Remote Sensing Data Set for Multi-modal Multi-temporal
Cloud Removal

Patrick Ebel, Yajin Xu, Michael Schmitt, and Xiao Xiang Zhu, This work was
partially supported by the Federal Ministry for Economic Affairs and Energy of
Germany in the project “AI4Sentinels– Deep Learning for the Enrichment of
Sentinel Satellite Imagery” (FKZ50EE1910). The work of X. Zhu is jointly
supported by the European Research Council (ERC) under the European Union’s
Horizon 2020 research and innovation programme (grant agreement No.
[ERC-2016-StG-714087], Acronym: So2Sat), by the Helmholtz Association through
the Framework of Helmholtz AI (grant number: ZT-I-PF-5-01) - Local Unit
“Munich Unit @Aeronautics, Space and Transport (MASTr)” and Helmholtz
Excellent Professorship “Data Science in Earth Observation - Big Data Fusion
for Urban Research”(grant number: W2-W3-100), by the German Federal Ministry
of Education and Research (BMBF) in the framework of the international future
AI lab ”AI4EO – Artificial Intelligence for Earth Observation: Reasoning,
Uncertainties, Ethics and Beyond” (grant number: 01DD20001) and by German
Federal Ministry of Economics and Technology in the framework of the ”national
center of excellence ML4Earth” (grant number: 50EE2201C).
(Corresponding Author: Xiao Xiang Zhu)P. Ebel and Y. Xu are with Data Science
in Earth Observation, Technical University of Munich, 80333 Munich, Germany.
(e-mail: patrick.ebel@tum.de).M. Schmitt is with the Remote Sensing Technology
Institute, German Aerospace Center, 82234 Wessling, Germany, and also with the
Chair of Earth Observation, Bundeswehr University Munich, 85577 Neubiberg,
Germany . (e-mail: michael.schmitt@unibw.de).X. X. Zhu is with Data Science in
Earth Observation, Technical University of Munich, 80333 Munich, Germany and
also with the Remote Sensing Technology Institute, German Aerospace Center,
82234 Wessling, Germany. (e-mail: xiaoxiang.zhu@dlr.de).

###### Abstract

About half of all optical observations collected via spaceborne satellites are
affected by haze or clouds. Consequently, cloud coverage affects the remote
sensing practitioner’s capabilities of a continuous and seamless monitoring of
our planet. This work addresses the challenge of optical satellite image
reconstruction and cloud removal by proposing a novel multi-modal and multi-
temporal data set called SEN12MS-CR-TS. We propose two models highlighting the
benefits and use cases of SEN12MS-CR-TS: First, a multi-modal multi-temporal
3D-Convolution Neural Network that predicts a cloud-free image from a sequence
of cloudy optical and radar images. Second, a sequence-to-sequence translation
model that predicts a cloud-free time series from a cloud-covered time series.
Both approaches are evaluated experimentally, with their respective models
trained and tested on SEN12MS-CR-TS. The conducted experiments highlight the
contribution of our data set to the remote sensing community as well as the
benefits of multi-modal and multi-temporal information to reconstruct noisy
information. Our data set is available at
https://patrickTUM.github.io/cloud_removal.

###### Index Terms:

Image Reconstruction, Cloud Removal, SAR-Optical, Data Fusion, Time Series,
Sequence-to-Sequence.

##  I Introduction

The majority of our planet’s land surface is covered by haze or clouds [1].
Such atmospheric distortions impede the capability of spaceborne optical
satellites to reliably and seamlessly record noise-free data of the Earth’s
surface. The presence of clouds is detrimental to typical remote sensing
applications, for instance land cover classification [2], semantic
segmentation [3, 4] and change detection [5, 6].

Figure 1: Example observations and cloud-free predictions. Columns: Samples at
two different time points. Rows: S1 data (in grayscale), cloudy S2 data (in
RGB), predicted cloud-free $\hat{S2}$ data, reference cloud-free S2 data of a
later point in time. The results highlight that our network is able to
integrate multi-modal and multi-temporal information to predict a clear-view
sequence of multi-spectral observations, even in the presence of heavy cloud
coverage.

Hence, the need for cloud-free earth observation gave rise to a rapidly
growing number of haze and cloud removal methods [7, 8, 9, 10, 3, 11, 12, 13,
14]. Most previous methods focus on a multi-modal approach [8, 15, 13, 14] to
reconstruct cloud-covered pixels via information translated from synthetic
aperture radar (SAR) or other sensors more robust to atmospheric disturbances
[16], yet focus on only a single time-point of observations. In comparison,
recent models attempt a temporal reconstruction of cloudy observations by
means of inference across time-series [12, 17, 18], utilizing the circumstance
that the extent of cloud coverage over a particular region is variable over
time and seasons [1].

Figure 2: Spatial distribution of the ROI constituting SEN12MS-CR-TS. Areas
belonging to the training split are denoted in blue, regions of the testing
split are colored in green. The ROI of SEN12MSCR [14], non-overlapping and
compatible with our data set, are depicted in gray. Any graphical overlap of
the semi-transparently plotted dots is rendered in darker tones so close-by
dots can easier be discerned.

The work at hand aims to combine both preceding approaches and thus considers
the challenge of cloud removal in optical satellite imagery by integrating
information across time and within different modalities. For this purpose, we
curate a new data set called SEN12MS-CR-TS, which contains multi-temporal and
multi-modal satellite observations. Specifically, SEN12MS-CR-TS consists of
1-year long time-series of co-registered radar Sentinel-1 (S1) as well as
multi-spectral Sentinel-2 observations (S2) acquired in a paired manner,
covering regions of interest (ROI) from all over the world. We highlight the
benefits of the proposed data set by training and testing two different models
on our data set: First, a multi-modal multi-temporal 3D-Convolution Neural
Network that predicts a cloud-free image from a sequence of cloudy optical and
radar images. Second, a sequence-to-sequence translation model that predicts a
cloud-free time series from a cloud-covered time series. Both approaches are
evaluated experimentally, with their respective models trained and tested on
SEN12MS-CR-TS. The conducted experiments highlight the contribution of our
curated data set to the remote sensing community as well as the benefits of
multi-modal and multi-temporal information to reconstruct noisy information.

Figure 3: Example data, preprocessed as stated in section II-B. Rows: S1 data
(in grayscale), S2 data (in RGB), binary cloud masks (as per s2cloudless
[19]). Columns: Samples of five different time points. The illustrations show
that the observed region is affected by variable atmospheric disturbances and
covered by a dynamic extent of clouds, changing over time. The detected cloud
coverage at the individual time points is 36, 49, 23, 48 percent, with an
average of about 39 percent across all illustrated samples. While some pixels
are clear at least at one point in the series and may thus be reconstructed by
integrating across time, others are cloud-covered throughout the sequence and
require spatial context or cloud-robust sensor information to be
reconstructed.  Figure 4: Statistics of cloud coverage of SEN12MS-CR-TS train
split, computed on full-scene images via the detector of [19]. On average,
approximately $44\%$ $(\pm 42\%)$ of occlusion is observed. The empirical
distribution of cloud coverage is bimodal and ranges from cloud-free views to
total occlusion. Figure 5: Statistics of cloud coverage of SEN12MS-CR-TS test
split, computed on full-scene images via the detector of [19]. On average,
approximately $50\%$ $(\pm 42\%)$ of occlusion is observed. The empirical
distribution of cloud coverage is bimodal and ranges from cloud-free views to
total occlusion.

###  I-A Related Work

As the presence of clouds in optical satellite imagery poses a severe
hindrance for remote sensing applications, there has been plenty of preceding
research on cloud removal methods [7, 8, 9, 20, 10, 3, 12, 13, 14]. The focus
of this overview is on data sets for cloud removal methods. Much of the early
work on cloud removal considered data of simulating cloudy observations [3].
Copying cloudy pixel values from one image to another clear-view one [3]
captures the spectral properties of naturally cloudy observations more
faithfully than synthetic noise (e.g. Perlin noise [21]) [7, 8, 20], but
neither precisely reproduce the statistics of satellite images containing
natural cloud occurrences [14]. Consequently, our data set contain cloud-free
as well as naturally occurring cloud-covered optical satellite recordings. The
SEN12MS-CR data set [14] provides a globally distributed collection of co-
registered mono-temporal Sentinel-1 as well as cloudy and cloud-free
Sentinel-2 observations. Our data set is an extension of SEN12MS-CR in the
sense that we collect repeated measures per ROI and thereby provide a time-
series of co-registered S1 and S2 observations, gathered such that matched
observations of both modalities are no more than two weeks apart. In
comparison to the preceding data set, ours allows integrating information not
solely across different sensors, but also across different points in time
distributed throughout the year. Similarly, the work of [12] allows for time-
series cloud removal by providing a collection of tri-temporal
RGB(NIR)-channel optical data and corresponding models. Our contribution
extends this work by providing true multi-modal data recorded by two distinct
sensors, synthetic aperture radar Sentinel-1 measurements as well as 13-bands
multi-spectral Sentinel-2 observations. Furthermore, the length of each time
series is increased considerably, from 3 to 30 samples. Finally, [12] exclude
observations with greater than $30$ percent cloud coverage from their data
set, which deviates from real conditions. Our approach aims to model the
complete spectrum of cloud coverage, including conditions commonly encountered
by remote sensing practitioners. In sum, our work and its main contribution, a
large-scale multi-modal multi-temporal data set for cloud removal in optical
satellite imagery, build on a history of research and improve upon the current
state of image reconstruction in remote sensing by providing a novel,
carefully curated data set.

##  II Data

This work introduces SEN12MS-CR-TS, a multi-modal and multi-temporal data set
for training and evaluating global and all-season cloud removal methods. The
data set consists of 53 globally distributed regions of interest, curated as
detailed in section II-A. The ROI are over $4000\times 4000$ $px^{2}$ each,
covering about $40\times 40$ $km^{2}$ of land such that the total surface area
covered by the data set is over $80,000$ $km^{2}$. Of all collected regions of
interest, 40 are defined as a training split and 13 as a hold-out split to
evaluate cloud removal approaches on. For every ROI we collect 30 co-
registered and paired S1 and S2 full-scene images evenly spaced in time
throughout the year of 2018. Each acquired image is inspected and quality-
controlled manually. The cloud-free Sentinel-2 (RGB-channel) observations of
four example ROI illustrating the diversity of our data set are illustrated in
Fig. 6. The spatial distribution of all ROI is depicted in Fig. 2 and
highlights the global sampling of our data set. The empirical distribution of
the cloud coverage of all optical observations is computed as detailed in
section II-C and the statistics are presented in Fig. 4 and 5 for the train
and the test splits, respectively. Importantly, the data set is curated
without excluding any interval of cloud coverage such that the collected
observations also reflect conditions of high cloud coverage as commonly
encountered in practice [1]. The data is made available under
https://patrickTUM.github.io/cloud_removal. It is about 2 Tb in size and
compatible with the SEN12MS-CR data set [14]. That is, no train ROI of
SEN12MS-CR are part of our data set’s test ROI and vice versa.

Figure 6: Four different regions contained in SEN12MS-CR-TS, highlighting the
diversity of sampled landcovers. The depicted S2 observations (RGB channels)
are cloud-free samples of their respective time series. The average ROI covers
about $40\times 40$ $km^{2}$ and is split into over 700 patch samples, each
patch of size $256\times 256\>px^{2}$.

###  II-A Data collection

All curated data are recorded via the SAR Sentinel-1 and multi-spectral
Sentinel-2 (level 1-C top-of-atmosphere reflectance product) instruments of
ESA’s Copernicus mission. The recorded observations are acquired via Google
Earth Engine [22] and a custom semi-automatic processing pipeline. We randomly
sample the geospatial locations of $53$ regions of interest from SEN12MS-CR
[14]. To minimize mosaicing, observations of cells covered by a single pass
are collected. The samples are referenced within the World Geodetic System
1984 (WGS84) coordinate system. For every ROI, 30 time intervals are evenly
spaced throughout the year of 2018. For every time interval a co-registered,
geo-referenced and full-scene S1 image as well as a paired full-scene S2 image
(level 1-C) are collected. The acquisition within the same interval window is
such that corresponding multi-modal images are no more than two weeks apart.
Further statistics regarding the pairing of observations are provided in
appendix A.

###  II-B Preprocessing

To prepare the collected raw data and translate it into a format that neural
networks for cloud removal can handle the following preprocessing steps are
taken: Each band of every observation is upsampled to 10m resolution (i.e. to
the native resolution of Sentinel-2’s bands 2,3,4 and 8). Every full-scene
image is sliced into non-overlapping patches of dimensions $256\times
256px^{2}$. The S1 observations are processed via the Sentinel-1 toolbox [23]
(including border and thermal noise removal, radiometric calibration,
orthorectification) and decibel-transformed. An example patch-wise tuple of
paired S1 and S2 data is illustrated in Fig. 3. Input patches to any ResNet
model [24] are preprocessed in line with the pipeline of [13] as follows: the
VV, VH channels of $S1$ observations are value-clipped in the ranges
$[-25;0],[-32.5;0]$ and rescaled to the interval $[0;2]$, while $S2$ patches
are value-clipped to $[0;10000]$ and normalized to the range $[0;5]$. For all
other networks with a different backbone architecture, preprocessing is done
as follows: each patch is value-clipped and then rescaled for every pixel to
take normalized values within the unit range of $[0,1]$. The modalities $S1$
and $S2$ are value-clipped within the intervals of $[-25;0]$ and $[0;10000]$,
respectively. This way, we follow the preprocessing protocol of the preceding
work and avoid any unnecessary adjustments, for the sake of simplicity. For
evaluation, the pixel values of all input patches, target images and
predictions are re-mapped to the unit interval $[0,1]$, where the goodness of
predictions is assessed according to the metrics stated in section IV-A.

###  II-C Cloud detection & mask computation

In order to analyse the statistics of cloud coverage in SEN12MS-CR-TS and to
model the spatio-temporal extent of clouds we compute binary cloud masks $m$.
For each optical image, the masks $m$ are computed on-the-fly via the cloud
detector of s2cloudless [19], which provides a binary mask of pixel-wise
values in $\left\\{0,1\right\\}$ that indicate cloud-free and cloud-covered
pixels, respectively. The cloud mask accuracy of s2cloudless is reported to be
on par with the multi-temporal classifier MAJA [25], but the considered
detector can be applied on mono-temporal satellite observations. Note that,
alternatively to s2cloudless, the masks $m$ may be computed via a dedicated
neural network for cloud detection [26, 27]. However, s2cloudless has proven
to be lightweight and provide sufficient performance at little extra
computational cost in run time or memory, making it an appealing cloud
detector to be applied on a large-scale date set such as SEN12MS-CR-TS.
Example cloud detections are illustrated in Fig. 3.

##  III Methods

We consider two distinctively different methods to highlight the benefits of
our curated data set and the diverse tasks it allows to approach. The first
method is a neural network reconstructing cloud covered pixels in time series
of multi-modal data to predict a single target image acquired at a cloud-free
time point. The second approach introduces a neural network that performs
sequence-to-sequence cloud removal, that is, it predicts a time series of
cloud-free observations the size of the cloudy input sequence.

###  III-A Multi-temporal multi-modal cloud removal

For multi-temporal multi-modal cloud removal we consider a deep neural network
that builds on the generator of [12]. Our model receives a sequence of
$t=1,...,n$ input tuples $\left(S1,S2\right)_{t}$ and predicts a cloud-removed
multispectral image $\hat{S2}$. The architecture of the proposed model uses a
ResNet [24] backbone, with Siamese residual branches processing the individual
time points until their information gets integrated. That is, we replaced the
pairwise concatenation of 2D feature maps in [12] by stacking features in the
temporal domain, followed by 3D convolutions. Moreover, as the first part of
the generator of [12] is effectively a single time point cloud removal sub-
network (as each time point is processed individually up to this point), we
substitute this component by the established ResNet-based [24] cloud removal
network of [13]. Subsequently, the feature maps are stacked in the temporal
dimension and 3D convolutions are applied to integrate information across
time. The output of the network is a single cloud-free image prediction
$\hat{S2}$. A schematic overview of the described architecture is shown in
Fig. 7.

Figure 7: A conceptual illustration of the sequence-to-point cloud removal
architecture $G_{seq2point}$. The network is based on the architecture of [12]
and consists of $n$ siamese ResNet branches [13] doing single time point cloud
removal on $n$ individual time points. Subsequently, the feature maps are
stacked in the temporal dimension and 3D convolutions are applied to integrate
information across time. The output of the network is a single cloud-free
image prediction.

###  III-B Internal Learning for sequence-to-sequence cloud removal

The sequence-to-sequence cloud removal method [28] follows the 3D Encoder-
Decoder architecture of [29], constituted of an encoder as well as a decoder
component. Both components are arranged symmetrically in the style of U-Net
[30] and linked via skip connections between paired layers. The input to the
network is a sequence of multi-temporal $S1$ samples and its output is a
sequence of multi-temporal cloud-removed $S2$ predictions. With regards to its
input-to-output mapping, the proposed architecture resembles earlier SAR-to-
optical translation method [31, 32]. Similar to these earlier domain
translation approaches, our network learns information of the target domain
(i.e. the optical imagery) via the supervision signal. Different from these
approaches, the internal learning framework described below removes clouds and
directly learns to denoise the target image sequence.

The architecture of the network is summarized in Fig. 8. Note, that the key
difference between the given model and the sequence-to-point method of section
IV-B (depicted in Fig. 7) is in the output dimensions: Whereas the sequence-
to-point architecture maps a sequence of $n$ cloudy inputs to a single cloud-
removed prediction, the sequence-to-sequence approach preserves the temporal
information by mapping to a time series of $n$ cloud removed outputs.
Moreover, the point estimator receives tuples of $S1$ and $S2$ inputs, whereas
the network of Fig. 8 is driven solely by $S1$ data (or Gaussian noise, as
proposed in [33, 29]). Finally, the sequence-to-point network of Fig. 7 builds
on the Siamese architecture of [12] with a ResNet backbone [13] plus 3D
convolutions, whereas the sequence-to-sequence approach of Fig. 8 follows a 3D
convolutional variant of U-Net [30], as proposed in [29].

Figure 8: A conceptual illustration of the 3D Encoder-Decoder architecture
$G_{seq2seq}$ employed in the sequence-to-sequence cloud removal model [28].
The network is based on the architecture of [29] and consists of encoder and
decoder parts arranged symmetrically in the style of U-Net [30], with skip
connections between paired layers. Input to the network is a batch of multi-
temporal $S1$ observations. The output is a predicted batch of multi-temporal
multi-spectral $S2$ observations. For the ablation model considered in section
IV-D, Gaussian noise is used as an input as in [33, 29].

The training procedure of the sequence-to-sequence network follows that of
internal learning for image inpainting [33, 29], which is formalized in
Algorithm 1. In this framework, for a given target sequence, a neural network
is trained from scratch directly on the target sequence (without any need for
additional or cloud-free training data) in order to reconstruct its noisy
pixels. The observations exhibit spatio-temporal regularities and patterns
(i.e. signal in the data), which is first modeled and learned by the network.
The irregularities in the sequence (i.e. noise in the target data) are only
internalized after, similar to a conventionally-trained network overfitting to
noise on training data. The internal learning approach exploits this signal-
noise dichotomy and teaches a model to reconstruct cloud-covered pixels in the
target sequence of $S2$ observations, without need for any external or cloud-
free training data. In detail, a neural network is initialized and trained
from scratch directly on the target sequence. At each iteration, the model
receives input driving its activations (e.g. Gaussian noise or $S1$
recordings) and predicts a sequence $\hat{S2}$. The predictions $\hat{S2}$ are
compared against the target sequence $S2$ (e.g. according to a cost function
${L}_{all}$ as in eq. 4) and the network learns to reproduce the cloud-free
pixels. The training stops before the network overfits to internalizing the
cloudy pixels.

Algorithm 1 Internal Learning to Remove Clouds

1:procedure Seq2SeqDeclouding($S1,S2,iterMax$)

2: $G_{S1\rightarrow S2}$ = init. new NeuralNetwork()

3: iterCount = $0$

4: while iterCount $<$ iterMax do

5: $\hat{S2}=G_{S1\rightarrow S2}(S1)$

6: $G_{S1\rightarrow S2}.backpropagate(\mathcal{L}_{all}(S2,\hat{S2}))$

7: iterCount = iterCount + 1

8: Return $\hat{S2}$

With respect to its application and functionality, our sequence-to-sequence
neural network resembles classical low-rank & sparse signal decomposition
methods [34, 35, 36, 37]: First, while neural networks are typically trained
on a dedicated training data set separated from the test observations, numeral
signal decomposition methods can be directly utilized on the data of interest.
Similarly, our model can be directly applied on the test data. Second,
unmixing of signals is very generic and can be applied to matrices as well as
tensors. In comparison, the deep image prior approach applies to single images
as well as time series [33, 29], too. Finally, the decomposition itself is
into a low-rank part and a sparse component. The low-rank part denotes the
data‘s compact representation and regularities. That is, spatial, spectral or
temporal (auto-)correlations such as the land cover mapped by a satellite. The
sparse component consists of the irregular part of the data which has only a
few non-zero entries, such as the appearance of clouds. In comparable terms,
the internal learning technique allows our network to discover the
regularities in the data and generalizing it to cloud-covered samples, before
overfitting to the noise.

##  IV Experiments and Results

This methods details the experimental design and the corresponding results on
the considered cloud removal methods as well as their ablation variants.
Section IV-A specifies the measures of goodness used to assess the quality of
the individual techniques’ predictions. Section IV-B introduces the baselines
compared against the proposed model of IV-B on the sequence-to-point cloud
removal task. Sections IV-C and IV-D detail the experiments and outcomes for
the sequence-to-point and sequence-to-sequence cloud removal tasks,
respectively.

###  IV-A Metrics

We evaluate quantitative performance in terms of normalized root mean squares
error (NRMSE), Peak Signal-to-Noise Ratio (PSNR), structural similarity (SSIM)
[38] and Spectral Angle Mapper (SAM) [39], defined as

$NRMSE(x,y)=\sqrt{\frac{1}{C\cdot H\cdot
W}\sum_{c=h=w=1}^{C,H,W}(x_{c,h,w}-y_{c,h,w})^{2}}$ $PSNR(x,y)=20\cdot
log_{10}\left(\frac{1}{NRMSE(x,y)}\right)$
$SSIM(x,y)=\frac{(2\mu_{x}\mu_{y}+\epsilon_{1})(2\sigma_{xy}+\epsilon_{2})}{(\mu_{x}+\mu_{y}+\epsilon_{1})(\sigma_{x}+\sigma_{y}+\epsilon_{2})}$
$SAM(x,y)=cos^{-1}\left(\frac{\sum^{C,H,W}_{c=h=w=1}x_{c,h,w}\cdot
y_{c,h,w}}{\sqrt{\sum^{C,H,W}_{c=h=w=1}x_{c,h,w}^{2}\cdot\sum^{C,H,W}_{c=h=w=1}y_{c,h,w}^{2}}}\right)$

with images $x,y$ compared via their respective pixel-values
$x_{c,h,w},y_{c,h,w}\in[0,1]$, dimensions $C=3$, $H=W=256$, means
$\mu_{x},\mu_{y}$, standard deviations $\sigma_{x},\sigma_{y}$, covariance
$\sigma_{xy}$ as well as constants $\epsilon_{1}$, $\epsilon_{2}$ to stabilize
the computation. NRMSE belongs to the class of pixel-level metrics and
quantifies the average discrepancy between target and predicted pixels in
units of the measure of interest. PSNR is evaluated on the whole image and
quantifies the signal-to-noise ratio of the prediction as a reconstruction of
the target image. SSIM is another image-wise measure that builds on PSNR and
captures the structural similarity of the prediction to the target in terms of
perceived change, contrast and luminance [38]. The SAM measure is a third
image-level metric that provides the spectral angle between the bands of two
multi-channel images [39]. For further analysis, the pixelwise NRMSE is
evaluated in three manners: 1) over all pixels of the target image (as per
convention), 2) only over cloud-covered pixels (visible in neither of any
input optical sample) to measure reconstruction of noisy information, as well
as 3) only over cloud-free pixels (visible in at least one input optical
patch) quantifying preservation of information. The pixel-wise masking is
performed according to the cloud mask given by the detector of [19].

###  IV-B Baseline Methods

To put the performances of our proposed model and ablations into context, we
consider the following baseline methods. First (”least cloudy”), taking the
least-cloudy input observation and forwarding it without further modification
to be compared against the cloud-free target image. This provides a measure of
how hard the cloud removal task is with respect to the extent of cloud-
coverage present in the data. Second (”mosaicing”), we perform a mosaicing
method that averages the values of pixels across cloud-free time points,
thereby integrating information across time. That is, for any pixel, if there
is a single clear-view time point then its value is copied, for multiple
cloud-free samples the mean is formed and in case no cloud-free time point
exists, then a value of 0.5 is taken as a proxy. This is to avoid any extreme
values, such as cloudy pixels of high-intensity. The mosaicing technique
provides a measure of how much information can be reconstructed across time,
from multi-spectral optical observations exclusively. Third, ResNet refers to
a residual neural network as described and trained in sections IV-B as well as
IV-C. The architecture is based on the model of [13], and serves as a relevant
baseline because parts of this model are used as Siamese residual branches
within our model, as detailed in section IV-B. It provides an estimate of how
well a point-to-point cloud removal model can perform as a baseline. Fourth,
the baseline STGAN denotes the ”Branched ResNet generator (IR)” architecture
of [12]. It is a sequence-to-point cloud removal model, and the architecture
of our own sequence-to-point neural network closely follows its design, as
detailed in section . In sum, the purpose of assessing these baselines is to
analyze whether trivial solutions to the multi-modal multi-temporal sequence-
to-point cloud removal problem exist, and how any more sophisticated deep
learning approach compares against these methods and our proposed model
trained on SEN12MS-CR-TS.

###  IV-C Sequence-to-Point Cloud Removal

This section details the training specifics of the sequence-to-point cloud
removal architecture introduced in section IV-B. As detailed in section IV-B,
up to the temporal concatenation layer, we use a version of the ResNet-based
[24] cloud removal network of [13] and pre-trained it on SEN12MS-CR [14]
according to the training specifics of [13]. All our considered sequence-to-
point cloud removal networks and ablation models share this pre-trained
single-temporal cloud removal network as a starting point for the sake of
comparability and in order to reduce the duration of training. The networks
are trained for a total of 10 epochs on 1 tuple of patches per location for
every ROI in the training split. For training, the input $S2$ patches are
filtered to display within $0$ to $50$ percent of cloud coverage. The target
$S2$ patch is selected to be the sample showing the minimum cloud coverage
over the given time series, i.e. it is not necessarily temporally preceding or
following the input patches. For the first 25000 steps in the training
procedure, the networks are trained with the initial ResNet Siamese components
frozen, exclusively optimizing the subsequent 3D-convolution layers. After the
steps with the pre-trained weights frozen and once the deeper layers have been
calibrated to the initial network’s latent feature maps, the full network is
trained end-to-end for the remainder of the process. During training, the
network minimizes the loss $\mathcal{L}_{\text{all}}$

$\mathcal{L}_{\text{all}}=\lambda_{L1}\mathcal{L}_{L1}+\lambda_{perc}\mathcal{L}_{perc}$
(1)

$\mathcal{L}_{\text{L1}}=||S2-\hat{S2}||_{1}$ (2)

$\mathcal{L}_{perc}=||VGG16(S2),VGG16(\hat{S2})||_{2}$ (3)

with $\lambda_{L1}=100$ according to [12] and $\lambda_{perc}=1$ as
hyperparameters weighting the individual pixel-wise loss
$\mathcal{L}_{\text{L1}}$ and the perceptual loss $\mathcal{L}_{perc}$. The
perceptual loss is computed by means of an auxiliary VGG16 network [40]
resulting in sharper image reconstructions [41]. In comparison to other VGG16
pre-trained on classical computer vision data sets such as ImageNet [42] and
thus limited to RGB channel data, we pre-trained a VGG16 for landcover
classification on the SEN12MS data set [43] according to the training protocol
of [2]. The proposed sequence-to-point cloud removal network (and its ablation
variants) are optimized via ADAM [44], with a learning rate of $0.0002$ and
momentum parameters $\left[0.5,0.999\right]$ as in [12]. A batch size of $1$
tuple of samples per iteration is used for training.

To evaluate performances on the test split, samples containing $S2$
observations from the complete range of cloud coverage (between $0$ and
$100\%$) are considered for input. Table I compared the results of our
proposed model with the baselines detailed in section IV-B. The results show
that the proposed network outperforms the baselines in the majority of
metrics, except for PSNR (where mosaicing comes first) and the NRMSE (clear)
preservation metric (where the ”least cloudy” approach performs best). This
demonstrates that a deep neural network approach can typically outperform
trivial solutions to the multi-modal multi-temporal cloud removal problem.
Exemplary outcomes for the considered baselines on four different samples from
the test split are presented in Fig. 9. The considered cases are cloud-free,
partly-cloudy, cloud-covered with no visibility except for a single time point
and cloud-coverage with no visibility at any time point. The results show that
the considered models typically outperform the simple heuristics. One
exceptional case is least cloudy in the absence of clouds, which manages to
accomplish a faithful prediction in such settings. Moreover, the illustrations
underline that multi-temporal and multi-modal data may benefit image
reconstruction: While most methods perform well in the cloud-free or partly-
cloudy cases, multi-source integration is needed if individual time points
contain dense cloud coverage over wide areas. When all input data is covered
by thick clouds, then this poses a severe challenge for all approaches
considered. To analyze the benefits of including $S1$ SAR data, we perform an
ablation study and compare a multi-sensor model against one only utilizing
multi-spectral $S2$ input. Table II compared the results of the multi-modal
model with an ablation version not using $S1$ SAR data. The comparison
illustrates the benefits of including SAR data when reconstructing cloud-
covered pixels. Next, we conduct an ablation experiment to assess the
additional benefits of utilizing the introduced perceptual loss. Table III
compared the results of our proposed model with an ablation version not using
the perceptual loss (i.e. setting $\lambda_{perc}=0$ in eq 1). The outcomes
imply that the usage of a perceptual loss results in cloud-removed predictions
of a higher quality. Finally, we consider the extension of the proposed model
into networks integrating 4 and 5 time points of input information. Table IV
compared the performance of our model as a function of input time points
($n=3,4,5$). The results indicate that considering longer time series may
provide further improvements in terms of reconstructing cloud-covered
information. In a final experiment on sequence-to-point cloud removal, Table V
reports the performance of our proposed model(n=3, with S1 and perceptual
loss) as a function of cloud coverage. That is, for a given interval of cloud
coverage, all $n=3$ input images are sampled to contain a corresponding extent
of clouds. The outcomes show that image reconstruction performance is highly
dependant on the percentage of cloud coverage. While performance decrease is
not strictly monotonous with an increase in cloud coverage, a strong
association persists.

model | NRMSE (all) | NRMSE (cloudy) | NRMSE (clear) | PSNR | SSIM | SAM
---|---|---|---|---|---|---
least cloudy | 0.079 | 0.082 | 0.031 | — | 0.815 | 0.213
mosaicing | 0.062 | 0.064 | 0.036 | 31.68 | 0.811 | 0.250
ResNet | 0.060 | 0.062 | 0.040 | 26.04 | 0.810 | 0.212
STGAN | 0.057 | 0.059 | 0.050 | 25.42 | 0.818 | 0.219
ours (n=3) | 0.051 | 0.052 | 0.040 | 26.68 | 0.836 | 0.186

Table I: Quantitative evaluation of the proposed sequence-to-point model with
baseline approaches in terms of normalized root mean squared error (NRSME),
peak signal-to-noise ratio (PSNR), structural similarity [38] (SSIM) and the
Spectral Angle Mapper (SAM) [39] metric. Our model performs best in the
majority of metrics, demonstrating that a deep neural network approach yields
additional benefits over trivial solutions to the multi-modal multi-temporal
cloud removal problem.

Figure 9: Exemplary predictions and cloud-free target images for all baselines
reported in Table I. Columns: Four different samples from the test split. The
considered cases are cloud-free, partly-cloudy, cloud-covered with no
visibility except for a single time point and cloud-covered with no visibility
in any time point. Rows: Predictions of least cloudy, mosaicing, ResNet,
STGAN, ours (n=3) as well as the cloud-free reference image. The results show
that the considered models outperform the simple heuristics. Moreover, the
illustrations underline that multi-temporal and multi-modal data may benefit
image reconstruction.

model | NRMSE (all) | NRMSE (cloudy) | NRMSE (clear) | PSNR | SSIM | SAM
---|---|---|---|---|---|---
ours (no S1) | 0.054 | 0.057 | 0.054 | 25.35 | 0.832 | 0.194
ours (with S1) | 0.051 | 0.052 | 0.040 | 26.68 | 0.836 | 0.186

Table II: Comparison of the proposed sequence-to-point model including SAR
observations versus an ablation version without SAR observations in terms of
normalized root mean squared error (NRSME), peak signal-to-noise ratio (PSNR),
structural similarity [38] (SSIM) and the Spectral Angle Mapper (SAM) [39]
metric. The comparison illustrates the benefits of including SAR data when
reconstructing cloud-covered pixels.

model | NRMSE (all) | NRMSE (cloudy) | NRMSE (clear) | PSNR | SSIM | SAM
---|---|---|---|---|---|---
ours (no percept.) | 0.052 | 0.053 | 0.039 | 26.66 | 0.835 | 0.180
ours (with percept.) | 0.051 | 0.052 | 0.040 | 26.68 | 0.836 | 0.186

Table III: Comparison of the proposed sequence-to-point model including
perceptual loss versus an ablation version without perceptual loss in terms of
normalized root mean squared error (NRSME), peak signal-to-noise ratio (PSNR),
structural similarity [38] (SSIM) and the Spectral Angle Mapper (SAM) [39]
metric. The outcomes imply that the usage of a perceptual loss during training
results in cloud-removed predictions of a higher quality at test time.

model | NRMSE (all) | NRMSE (cloudy) | NRMSE (clear) | PSNR | SSIM | SAM
---|---|---|---|---|---|---
ours (n=3) | 0.051 | 0.052 | 0.040 | 26.68 | 0.836 | 0.186
ours (n=4) | 0.049 | 0.050 | 0.041 | 27.10 | 0.845 | 0.172
ours (n=5) | 0.048 | 0.048 | 0.032 | 27.07 | 0.846 | 0.178

Table IV: Quantitative evaluation of the proposed sequence-to-sequence model
with varying numbers of time points ($n=3,4,5$) in terms of normalized root
mean squared error (NRSME), peak signal-to-noise ratio (PSNR), structural
similarity [38] (SSIM) and the Spectral Angle Mapper (SAM) [39] metric. Our
multi-temporal network with SAR guidance outperforms the multi-temporal
ablation model without prior SAR information.

% cloud coverage | NRMSE (all) | NRMSE (cloudy) | NRMSE (clear) | PSNR | SSIM | SAM
---|---|---|---|---|---|---
0-10 % | 0.041 | 0.046 | 0.041 | 28.59 | 0.870 | 0.143
10-20 % | 0.044 | 0.046 | 0.043 | 27.69 | 0.848 | 0.166
20-30 % | 0.046 | 0.047 | 0.044 | 27.25 | 0.841 | 0.169
30-40 % | 0.048 | 0.050 | 0.045 | 26.77 | 0.830 | 0.169
40-50 % | 0.047 | 0.048 | 0.045 | 26.86 | 0.830 | 0.167
50-60 % | 0.049 | 0.494 | 0.048 | 26.55 | 0.825 | 0.185
60-70 % | 0.052 | 0.052 | 0.043 | 26.10 | 0.817 | 0.184
70-80 % | 0.049 | 0.050 | 0.044 | 26.59 | 0.816 | 0.179
80-90 % | 0.050 | 0.050 | 0.044 | 26.54 | 0.820 | 0.175
90-100 % | 0.063 | 0.063 | — | 24.79 | 0.786 | 0.222

Table V: Performance of our sequence-to-point cloud removal method (n=3, with
S1 & with perceptual loss) as a function of cloud coverage. For a given
interval, all n=3 input images are sampled to contain a corresponding extent
of clouds. The outcomes show that image reconstruction performance is highly
dependant on the percentage of cloud coverage. While performance decrease is
not strictly monotonous with an increase in cloud coverage, a strong
association persists.

###  IV-D Sequence-to-Sequence Cloud Removal

A key characteristic of training the sequence-to-sequence cloud removal model
described in section III-B is the model being trained directly on the time
series of images one aims to removes clouds from, without the use of any
external training data as in [33, 29]. More specifically, the training
procedure teaches the network to replicate cloud-free pixels and inpaint
cloud-covered ones in the target sequence $S2$ according to the cost function
$\mathcal{L}_{all}$ formulated in [29] as

$\mathcal{L}_{all}=\lambda_{L2}\mathcal{L}_{L2}+\lambda_{perc}\mathcal{L}_{perc}$
(4)

$\mathcal{L}_{L2}=||S2\cdot(1-m),\hat{S2}\cdot(1-m)||_{2}$ (5)

$\mathcal{L}_{perc}=||VGG16(S2)\cdot(1-m),VGG16(\hat{S2})\cdot(1-m)||_{2}$ (6)

where $\lambda_{L2}=1$ and $\lambda_{perc}=0.01$ refer to hyper-parameters
that linearly combine the terms constituting $\mathcal{L}_{all}$.
$\mathcal{L}_{2}$ is a pixel-wise reconstruction loss evaluated over the
cloud-free pixels via an auxiliary VGG16 network [40] as explained before. The
pseudo-code formalizing the intrinsic learning procedure is given in Algorithm
1 described in section III-B and further justifications are stated in the
original work of [33]. For a given target sequence, the network is trained for
20 passes with batches of $n=5$ samples consisting of temporally adjacent
images, for 100 iterations per pass. The network is optimized via ADAM [44]
with a learning rate of $0.01$ and the hyperparameters of Algorithm 1 set as
stated in [29].

To quantitatively evaluate the considered model on SEN12MS-CR-TS, we propose
the following protocol for a sequence-to-sequence cloud removal task: For a
given target sequence, the least cloud-covered S2 observation is identified
and denoted as a target image $S2_{t}$. The most cloudy S2 sample is observed
and denoted as a source image $S2_{s}$. The cloud-covered pixels of $S2_{s}$
according to a cloud mask $m$ are alpha-blended with the cloud-free pixels of
$S2_{t}$ similar to the approach of [3]. Finally, the cloud-removed prediction
$\hat{S2}_{t}$ is then compared against the originally cloud-free $S2_{t}$ in
order to get a measure of goodness of cloud removal.

Table VI shows the results of the proposed network on the sequence-to-sequence
cloud removal task following the aforementioned protocol. Furthermore, the
considered model is compared against an ablation model, conditioned on random
Gaussian noise as in [33, 29] in place of the meaningful $S1$ input
observations. Example outcomes of sequence-to-sequence cloud removal on a
given ROI are depicted in Fig. 1. Furthermore, Fig. 11 provides a qualitative
comparison between the predictions conditioned on SAR versus no prior
information, underlining the benefits of multi-modal information. The results
highlight that the internal learning approach can learn to reconstruct cloud-
covered pixels on a very limited amount of data. Furthermore, the results
demonstrate that including SAR data results in performance benefits over the
single-sensor baseline.

Figure 10: Illustration of baseline methods for the sequence-to-sequence cloud removal task. The presented results show a cloudy image to be de-clouded, as well as the predictions via Riemannian Robust Principal Component Pursuit (RPCP) [45], Non-negative Matrix Factorization Incremental Subspace Learning (NMFISL) [46], Probabilistic Non-negative Matrix Factorization (PNMF) [47], Manhattan Non-negative Matrix Factorization (MNMF) [48] and Online Stochastic Tensor Decomposition (OSTD) [49]. The results indicate that the presence of large and dense clouds poses a severe challenge for the considered methods. Most baselines de-cloud the image except for some residual artifacts, some techniques display discolorization. For comparison with ours (no S1), ours (with S1) and the cloud-free target image, see Fig. 11 model | NRMSE (all) | PSNR | SSIM | SAM
---|---|---|---|---
RPCP [45] | 0.403 | 7.911 | 0.264 | 30.567
NMFISL [46] | 0.312 | 10.262 | 0.450 | 29.285
PNMF [47] | 0.317 | 10.135 | 0.432 | 29.801
MNMF [48] | 0.361 | 8.945 | 0.361 | 28.685
OSTD [49] | 0.303 | 10.853 | 0.402 | 35.454
seq2seq (no S1) | 0.298 | 11.434 | 0.494 | 28.127
seq2seq (with S1) | 0.274 | 11.590 | 0.512 | 27.733
Table VI: Quantitative evaluation of baseline methods and the proposed
sequence-to-sequence model in terms of root mean squared error (RSME), peak
signal-to-noise ratio (PSNR), structural similarity (SSIM) and the Spectral
Angle Mapper (SAM) [39] metric. Our multi-temporal network with SAR guidance
outperforms the considered baselines as well as the multi-temporal ablation
model without prior SAR information.

Figure 11: Illustrations on the effect of prior guidance via SAR information.
Columns: SAR input to the SAR-conditioned model, cloud-free prediction of the
model conditioned on Gaussian noise, cloud-free prediction of the model
conditioned on SAR information, cloud-free observation as a reference image.
The structural information provided by the SAR input provides a strong prior
to the model, guiding it towards learning to remove clouds in the cloudy input
time series.

##  V Discussion

The main contribution of this work is in curating and providing SEN12MS-CR-TS,
a multi-modal multi-temporal data set for cloud removal in optical satellite
imagery. Our large-scale data set covers a heterogeneous set of ROI sampled
from all over Earth, acquired in different seasons throughout the year. Given
that the contained observations cover clear-view, filmy, as well as non-
transparent dense clouds, the objective of reconstructing cloud-covered
information poses a challenging task for the considered methods and future
approaches. For the sake of demonstrating the usefulness of the presented data
set, we propose a sequence-to-point as well sequence-to-sequence cloud removal
network. The considered methods are evaluated in terms of pixel-wise and
image-wise metrics. We provide evidence that taking time series information
into account is facilitating the reconstruction of cloudy pixels and that
including multi-sensor measurements does further improve the goodness of the
cloud-removed predictions, justifying the design of SEN12MS-CR-TS to include
multi-temporal and multi-modal data. The major difference to the preceding
mono-temporal SEN12MS-CR data set [15] for cloud removal is that SEN12MS-CR-TS
features a time series of 30 samples per ROI. This allows for developing
methods that integrate information across time to more faithfully reconstruct
cloud-obscured measurements. The sensitivity to temporal information may be
particularly valuable for future research investigating the benefits of cloud
removal to time-sensitive applications, such as change detection. On the other
side, there is a trade-off in terms of size, and while SEN12MS-CR-TS is more
than twice as large as its mono-temporal precursor, the latter contains about
two times as many ROI sampled over all continents. However, both data sets are
fully compatible; meaning that holdout ROI of one belong to the test split of
the other data set and vice versa. As there is no geo-spatial overlap across
splits between both data sets, they can be combined for training or validation
purposes. Finally, the two data sets exhibit a comparable extent of cloud
coverage—about 50 and 48 % respectively, both covering the full spectrum from
semitransparent haze to thick and dense clouds. A discrepancy between both
data sets is in SEN12MS-CR having between 25 and 50 % overlap between
neighboring patches (following the design of [43]), whereas SEN12MS-CR-TS has
no intersection between adjacent samples. SEN12MS-CR contains 122,218 patch
triplets of S1, cloudy S2 and cloud-free S2 data, whereas SEN12MS-CR-TS
consists of 30 time samples for each of the 15578 patch-wise observations, for
every S1 and S2 measurement. Due to the differences in pre-processing the two
data sets are not co-registered patch-wise but, importantly, they share a
common definition of ROI as well as train and test splits. This way, they are
compatible with one another such that SEN12MS-CR-TS can be utilized for time-
series cloud removal while SEN12MS-CR can provide further geospatial coverage
of additional ROI on individual time points. Thanks to the different designs
of both data sets, they may prove beneficial facilitating a variety of
downstream tasks, such as semantic segmentation [43], scene classification [2]
or change detection [5], even in the presence of clouds.

Beyond the design of our novel data set, additional contributions of this work
are in introducing the internal learning approach to cloud removal in optical
satellite data, as well as demonstrating that SAR-to-optical cloud removal
performs better than the original noise-to-optical translation framework.
While our data set aims to provide a global distribution of samples, we think
that the internal learning approach to cloud removal may be of particular
interest for remote sensing practitioners focusing on a single a spatially
confined ROI, as no further external data is necessary.

##  VI Conclusion

As a large extent of our planet is covered by haze or clouds at any given
point in time, such atmospheric distortions pose a severe constraint to the
ongoing monitoring of Earth. To approach this challenge, our work presented
SEN12MS-CR-TS, a multi-modal and multi-temporal data set for training and
evaluating global and all-season cloud removal methods. Our data set contains
Sentinel-1 and Sentinel-2 observations from over $80,000$ $km^{2}$ of
landcover, distributed globally and recorded through the year. The globally
distributed ROI each are large-sized, and capture a heterogeneous mass of
landcover. We demonstrated the practicality of SEN12MS-CR by considering two
methods: First, a model for sequence-to-point cloud removal. Second, a network
for sequence-to-sequence cloud removal which, to our knowledge, provides the
first case a model preserving temporal information is proposed in the context
of cloud removal. Both methods benefited from the presence of co-registered
and paired SAR measurements contained in our data set. The conducted
experiments highlight the contribution of our curated data set to the remote
sensing community as well as the benefits of multi-modal and multi-temporal
information to reconstruct noisy information. SEN12MS-CR is made public to
facilitate future research in multi-modal and multi-temporal image
reconstruction.

## Acknowledgment

The authors would like to thank ESA and the Copernicus program for making the
Sentinel observations accessed for this submission publicly available.
Furthermore, we would like to thank Rewanth Ravindran for assisting us in the
data curation process.

##  Appendix A Temporal coincidence of paired observations

Full-scene observations of Sentinel-1 and Sentinel-2 are collected within a
14-days time window in a paired manner, as specified in section II-A. To
further analyze the temporal distance within paired data, Fig. 12 illustrates
the empirically observed coincidences within SEN12MS-CR-TS. The mean time
differences across all paired observations is $2.61$ ($\pm$ $2.41$), which is
considerably smaller than the interval bound and implies a close proximity
between paired samples.

Figure 12: Histogram of temporal differences between paired observations. The
mean time differences across all paired observations is $2.61$ ($\pm$ $2.41$),
indicating a close proximity between paired samples.

## References

  * [1] M. D. King, S. Platnick, W. P. Menzel, S. A. Ackerman, and P. A. Hubanks, “Spatial and temporal distribution of clouds observed by modis onboard the terra and aqua satellites,” _IEEE Transactions on Geoscience and Remote Sensing_ , vol. 51, no. 7, p. 3826–3852, Jul 2013.
  * [2] M. Schmitt and Y.-L. Wu, “Remote sensing image classification with the SEN12MS dataset,” _ISPRS Ann. Photogramm. Remote Sens. Spatial Inf. Sci._ , vol. V-2-2021, pp. 101–106, 2021.
  * [3] M. Usman Rafique, H. Blanton, and N. Jacobs, “Weakly supervised fusion of multiple overhead images,” in _Proceedings of the IEEE Conference on Computer Vision and Pattern Recognition Workshops_ , 2019, pp. 0–0.
  * [4] M. Schmitt, J. Prexl, P. Ebel, L. Liebel, and X. X. Zhu, “Weakly supervised semantic segmentation of satellite images for land cover mapping–challenges and opportunities,” _arXiv preprint arXiv:2002.08254_ , 2020.
  * [5] P. Ebel, S. Saha, and X. X. Zhu, “Fusing multi-modal data for supervised change detection,” _The International Archives of Photogrammetry, Remote Sensing and Spatial Information Sciences_ , vol. 43, pp. 243–249, 2021\.
  * [6] S. Saha, P. Ebel, and X. X. Zhu, “Self-supervised multisensor change detection,” _IEEE Transactions on Geoscience and Remote Sensing_ , pp. 1–10, 2021.
  * [7] K. Enomoto, K. Sakurada, W. Wang, H. Fukui, M. Matsuoka, R. Nakamura, and N. Kawaguchi, “Filmy cloud removal on satellite imagery with multispectral conditional generative adversarial nets,” _2017 IEEE Conference on Computer Vision and Pattern Recognition Workshops (CVPRW)_ , p. 1533–1541, Jul 2017, arXiv: 1710.04835.
  * [8] C. Grohnfeldt, M. Schmitt, and X. Zhu, “A conditional generative adversarial network to fuse SAR and multispectral optical data for cloud removal from Sentinel-2 images,” in _IGARSS 2018 - 2018 IEEE International Geoscience and Remote Sensing Symposium_. IEEE, Jul 2018, p. 1726–1729. [Online]. Available: https://ieeexplore.ieee.org/document/8519215/
  * [9] P. Singh and N. Komodakis, “Cloud-GAN: Cloud removal for Sentinel-2 imagery using a cyclic consistent generative adversarial networks,” in _IGARSS 2018 - 2018 IEEE International Geoscience and Remote Sensing Symposium_. IEEE, Jul 2018, p. 1772–1775. [Online]. Available: https://ieeexplore.ieee.org/document/8519033/
  * [10] J. D. Bermudez, P. N. Happ, R. Q. Feitosa, and D. A. B. Oliveira, “Synthesis of multispectral optical images from SAR/optical multitemporal data using conditional generative adversarial networks,” _IEEE Geoscience and Remote Sensing Letters_ , vol. 16, no. 8, p. 1220–1224, Aug 2019.
  * [11] Z. Gu, Z. Zhan, Q. Yuan, and L. Yan, “Single remote sensing image dehazing using a prior-based dense attentive network,” _Remote Sensing_ , vol. 11, no. 24, p. 3008, 2019.
  * [12] V. Sarukkai, A. Jain, B. Uzkent, and S. Ermon, “Cloud removal from satellite images using spatiotemporal generator networks,” in _The IEEE Winter Conference on Applications of Computer Vision_ , 2020, pp. 1796–1805.
  * [13] A. Meraner, P. Ebel, X. X. Zhu, and M. Schmitt, “Cloud removal in sentinel-2 imagery using a deep residual neural network and SAR-optical data fusion,” _ISPRS Journal of Photogrammetry and Remote Sensing_ , vol. 166, pp. 333–346, 2020.
  * [14] P. Ebel, A. Meraner, M. Schmitt, and X. X. Zhu, “Multisensor data fusion for cloud removal in global and all-season sentinel-2 imagery,” _IEEE Transactions on Geoscience and Remote Sensing_ , 2020.
  * [15] P. Ebel, M. Schmitt, and X. X. Zhu, “Cloud removal in unpaired sentinel-2 imagery using cycle-consistent GAN and SAR-optical data fusion,” in _IGARSS 2020-2020 IEEE International Geoscience and Remote Sensing Symposium_. IEEE, 2020, pp. 2065–2068.
  * [16] R. Bamler, “Principles of synthetic aperture radar,” _Surveys in Geophysics_ , vol. 21, no. 2-3, pp. 147–157, 2000.
  * [17] S. Oehmcke, T.-H. K. Chen, A. V. Prishchepov, and F. Gieseke, “Creating cloud-free satellite imagery from image time series with deep learning,” in _Proceedings of the 9th ACM SIGSPATIAL International Workshop on Analytics for Big Geospatial Data_ , 2020, pp. 1–10.
  * [18] Q. Zhang, Q. Yuan, Z. Li, F. Sun, and L. Zhang, “Combined deep prior with low-rank tensor svd for thick cloud removal in multitemporal images,” _ISPRS Journal of Photogrammetry and Remote Sensing_ , vol. 177, pp. 161–173, 2021.
  * [19] A. Zupanc, “Improving cloud detection with machine learning,” Sentinel-Hub.https://medium.com/sentinel-hub/improving-cloud-detection-with-machine-learning-c09dc5d7cf13, 2017, accessed: 2019-10-10.
  * [20] W. Sintarasirikulchai, T. Kasetkasem, T. Isshiki, T. Chanwimaluang, and P. Rakwatin, “A multi-temporal convolutional autoencoder neural network for cloud removal in remote sensing images,” in _2018 15th International Conference on Electrical Engineering/Electronics, Computer, Telecommunications and Information Technology (ECTI-CON)_. IEEE, 2018, pp. 360–363.
  * [21] K. Perlin, “Improving noise,” in _Proceedings of the 29th annual conference on Computer graphics and interactive techniques_ , 2002, pp. 681–682.
  * [22] N. Gorelick, M. Hancher, M. Dixon, S. Ilyushchenko, D. Thau, and R. Moore, “Google Earth Engine: Planetary-scale geospatial analysis for everyone,” _Remote Sensing of Environment_ , vol. 202, pp. 18–27, 2017.
  * [23] L. Veci, P. Prats-Iraola, R. Scheiber, F. Collard, N. Fomferra, and M. Engdahl, “The Sentinel-1 toolbox,” in _Proceedings of the IEEE International Geoscience and Remote Sensing Symposium (IGARSS)_. IEEE, 2014, pp. 1–3.
  * [24] K. He, X. Zhang, S. Ren, and J. Sun, “Deep residual learning for image recognition,” in _Proceedings of the IEEE conference on computer vision and pattern recognition_ , 2016, pp. 770–778.
  * [25] V. Lonjou, C. Desjardins, O. Hagolle, B. Petrucci, T. Tremas, M. Dejus, A. Makarau, and S. Auer, “Maccs-atcor joint algorithm (maja),” in _Remote Sensing of Clouds and the Atmosphere XXI_ , vol. 10001. International Society for Optics and Photonics, 2016, p. 1000107.
  * [26] J. H. Jeppesen, R. H. Jacobsen, F. Inceoglu, and T. S. Toftegaard, “A cloud detection algorithm for satellite imagery based on deep learning,” _Remote sensing of environment_ , vol. 229, pp. 247–259, 2019.
  * [27] D. López-Puigdollers, G. Mateo-García, and L. Gómez-Chova, “Benchmarking deep learning models for cloud detection in landsat-8 and sentinel-2 images,” _Remote Sensing_ , vol. 13, no. 5, p. 992, 2021.
  * [28] P. Ebel, M. Schmitt, and X. X. Zhu, “Internal learning for sequence-to-sequence cloud removal via synthetic aperture radar prior information,” in _IGARSS 2021-2021 IEEE International Geoscience and Remote Sensing Symposium_. IEEE, 2021, p. In press.
  * [29] H. Zhang, L. Mai, N. Xu, Z. Wang, J. Collomosse, and H. Jin, “An internal learning approach to video inpainting,” in _Proceedings of the IEEE International Conference on Computer Vision_ , 2019, pp. 2720–2729.
  * [30] O. Ronneberger, P. Fischer, and T. Brox, “U-net: Convolutional networks for biomedical image segmentation,” in _International Conference on Medical image computing and computer-assisted intervention_. Springer, 2015, pp. 234–241.
  * [31] M. Fuentes Reyes, S. Auer, N. Merkle, C. Henry, and M. Schmitt, “SAR-to-optical image translation based on conditional generative adversarial networks—optimization, opportunities and limits,” _Remote Sensing_ , vol. 11, no. 17, p. 2067, Sep 2019.
  * [32] L. Wang, X. Xu, Y. Yu, R. Yang, R. Gui, Z. Xu, and F. Pu, “SAR-to-optical image translation using supervised cycle-consistent adversarial networks,” _IEEE Access_ , vol. 7, pp. 129 136–129 149, 2019.
  * [33] D. Ulyanov, A. Vedaldi, and V. Lempitsky, “Deep image prior,” in _Proceedings of the IEEE Conference on Computer Vision and Pattern Recognition_ , 2018, pp. 9446–9454.
  * [34] F. De la Torre and M. J. Black, “Robust principal component analysis for computer vision,” in _Proceedings Eighth IEEE International Conference on Computer Vision. ICCV 2001_ , vol. 1. IEEE, 2001, pp. 362–369.
  * [35] D. L. Donoho, “Compressed sensing,” _IEEE Transactions on information theory_ , vol. 52, no. 4, pp. 1289–1306, 2006.
  * [36] X. X. Zhu and R. Bamler, “Tomographic sar inversion by $l_{1}$-norm regularization—the compressive sensing approach,” _IEEE Transactions on Geoscience and Remote Sensing_ , vol. 48, no. 10, pp. 3839–3846, 2010.
  * [37] ——, “A sparse image fusion algorithm with application to pan-sharpening,” _IEEE Transactions on geoscience and remote sensing_ , vol. 51, no. 5, pp. 2827–2836, 2012.
  * [38] Z. Wang, A. Bovik, H. Sheikh, and E. Simoncelli, “Image quality assessment: From error visibility to structural similarity,” _IEEE Transactions on Image Processing_ , vol. 13, no. 4, p. 600–612, Apr 2004.
  * [39] F. A. Kruse, A. Lefkoff, J. Boardman, K. Heidebrecht, A. Shapiro, P. Barloon, and A. Goetz, “The spectral image processing system (SIPS)-interactive visualization and analysis of imaging spectrometer data,” in _AIP Conference Proceedings_ , vol. 283, no. 1. American Institute of Physics, 1993, pp. 192–201.
  * [40] K. Simonyan and A. Zisserman, “Very deep convolutional networks for large-scale image recognition,” _arXiv preprint arXiv:1409.1556_ , 2014.
  * [41] J. Johnson, A. Alahi, and L. Fei-Fei, “Perceptual losses for real-time style transfer and super-resolution,” in _European conference on computer vision_. Springer, 2016, pp. 694–711.
  * [42] O. Russakovsky, J. Deng, H. Su, J. Krause, S. Satheesh, S. Ma, Z. Huang, A. Karpathy, A. Khosla, M. Bernstein _et al._ , “ImageNet large scale visual recognition challenge,” _International journal of computer vision_ , vol. 115, no. 3, pp. 211–252, 2015.
  * [43] M. Schmitt, L. Hughes, C. Qiu, and X. Zhu, “SEN12MS-a curated dataset of georeferenced multi-spectral Sentinel-1/2 imagery for deep learning and data fusion,” _ISPRS Annals of Photogrammetry, Remote Sensing and Spatial Information Sciences_ , pp. 153–160, 2019.
  * [44] D. P. Kingma and J. Ba, “Adam: A method for stochastic optimization,” in _International Conference on Learning Representations( ICLR)_ , 2015.
  * [45] M. Hintermüller and T. Wu, “Robust principal component pursuit via inexact alternating minimization on matrix manifolds,” _Journal of Mathematical Imaging and Vision_ , vol. 51, no. 3, pp. 361–377, 2015.
  * [46] S. S. Bucak, B. Gunsel, and O. Gursoy, “Incremental nonnegative matrix factorization for background modeling in surveillance video,” in _2007 IEEE 15th Signal Processing and Communications Applications_. IEEE, 2007, pp. 1–4.
  * [47] D. Lee and H. S. Seung, “Algorithms for non-negative matrix factorization,” _Advances in neural information processing systems_ , vol. 13, 2000.
  * [48] N. Guan, D. Tao, Z. Luo, and J. Shawe-Taylor, “MahNMF: Manhattan non-negative matrix factorization,” _arXiv preprint arXiv:1207.3438_ , 2012\.
  * [49] A. Sobral, S. Javed, S. Ki Jung, T. Bouwmans, and E.-h. Zahzah, “Online stochastic tensor decomposition for background subtraction in multispectral video sequences,” in _Proceedings of the IEEE International Conference on Computer Vision Workshops_ , 2015, pp. 106–113.

|  Patrick Ebel<EMAIL_ADDRESS>received the B.Sc. degree in cognitive
science from University of Osnabrück, Germany, in 2015. He received the M.Sc.
degree in Cognitive Neuroscience and the M.Sc degree in Artificial
Intelligence from Radboud University Nijmegen, The Netherlands, in 2018. He is
working as doctoral candidate in the Data Science in Earth Observation lab,
Department of Aerospace and Geodesy of Technical University of Munich (TUM),
Germany. His research interests include machine learning as well as its
applications in computer vision and remote sensing. Specifically, he is
working on multi-modal and multi-temporal data fusion and automated image
reconstruction methods.
---|---
|  Yajin Xu<EMAIL_ADDRESS>received the B.Sc. degree in engineering with
distinctions from Wuhan University, Wuhan, China in 2018 and is pursuing an
M.Sc. double degree in Earth Oriented Space Science and Technology (ESPACE) at
Wuhan University, Wuhan, China and Technical University of Munich (TUM),
Munich, Germany. His study focus is machine learning applied to remote sensing
data. In 2021 he was a research assistant in the Data Science in Earth
Observation (SiPEO) group at TUM and Remote Sensing Technology Institute of
DLR, investigating deep learning-based approaches for cloud removal in optical
satellite data. His interests are in geospatial data analysis.
---|---
|  Michael Schmitt<EMAIL_ADDRESS>received his Dipl.-Ing. (Univ.)
degree in geodesy and geoinformation, his Dr.-Ing. degree in remote sensing,
and his habilitation in data fusion from the Technical University of Munich
(TUM), Germany, in 2009, 2014, and 2018, respectively. Since 2021, he has held
the Chair for Earth Observation at the Department of Aerospace Engineering of
the Bundeswehr University Munich in Neubiberg, Germany. Before that, he was a
professor for applied geodesy and remote sensing at the Munich University of
Applied Sciences, Department of Geoinformatics. From 2015 to 2020, he was a
senior researcher and deputy head at the Professorship for Data Science in
Earth Observation at TUM; in 2019 he was additionally appointed as Adjunct
Teaching Professor at the Department of Aerospace and Geodesy of TUM. In 2016,
he was a guest scientist at the University of Massachusetts, Amherst. His
research focuses on image analysis and machine learning applied to the
extraction of information from multi-modal remote sensing observations. In
particular, he is interested in remote sensing data fusion with a focus on SAR
and optical data. He is a co-chair of the Working Group “SAR and Microwave
Sensing” of the International Society for Photogrammetry and Remote Sensing,
and also of the Working Group “Benchmarking” of the IEEE-GRSS Image Analysis
and Data Fusion Technical Committee. He frequently serves as a reviewer for a
number of renowned international journals and conferences and has received
several Best Reviewer awards. He is a Senior Member of the IEEE and an
associate editor of IEEE Geoscience and Remote Sensing Letters.
---|---
|  Xiao Xiang Zhu (S’10–M’12–SM’14–F’21) received the Master (M.Sc.) degree,
her doctor of engineering (Dr.-Ing.) degree and her “Habilitation” in the
field of signal processing from Technical University of Munich (TUM), Munich,
Germany, in 2008, 2011 and 2013, respectively. She is currently the Professor
for Data Science in Earth Observation (former: Signal Processing in Earth
Observation) at Technical University of Munich (TUM) and the Head of the
Department “EO Data Science” at the Remote Sensing Technology Institute,
German Aerospace Center (DLR). Since 2019, Zhu is a co-coordinator of the
Munich Data Science Research School (www.mu-ds.de). Since 2019 She also heads
the Helmholtz Artificial Intelligence – Research Field “Aeronautics, Space and
Transport”. Since May 2020, she is the director of the international future AI
lab ”AI4EO – Artificial Intelligence for Earth Observation: Reasoning,
Uncertainties, Ethics and Beyond”, Munich, Germany. Since October 2020, she
also serves as a co-director of the Munich Data Science Institute (MDSI), TUM.
Prof. Zhu was a guest scientist or visiting professor at the Italian National
Research Council (CNR-IREA), Naples, Italy, Fudan University, Shanghai, China,
the University of Tokyo, Tokyo, Japan and University of California, Los
Angeles, United States in 2009, 2014, 2015 and 2016, respectively. She is
currently a visiting AI professor at ESA’s Phi-lab. Her main research
interests are remote sensing and Earth observation, signal processing, machine
learning and data science, with a special application focus on global urban
mapping. Dr. Zhu is a member of young academy (Junge Akademie/Junges Kolleg)
at the Berlin-Brandenburg Academy of Sciences and Humanities and the German
National Academy of Sciences Leopoldina and the Bavarian Academy of Sciences
and Humanities. She serves in the scientific advisory board in several
research organizations, among others the German Research Center for
Geosciences (GFZ) and Potsdam Institute for Climate Impact Research (PIK). She
is an associate Editor of IEEE Transactions on Geoscience and Remote Sensing
and serves as the area editor responsible for special issues of IEEE Signal
Processing Magazine. She is a Fellow of IEEE.
---|---