# Self-Attention for Raw Optical Satellite Time Series Classification

Marc Rußwurm and Marco Körner  Chair of Remote Sensing Technology,
Department of Aerospace and Geodesy,
Technical University of Munich
Arcisstraße 21, 80333 Munich, Germany

###### Abstract

The amount of available Earth observation data has increased dramatically in
the recent years. Efficiently making use of the entire body information is a
current challenge in remote sensing and demands for light-weight problem-
agnostic models that do not require region- or problem-specific expert
knowledge. End-to-end trained deep learning models can make use of raw sensory
data by learning feature extraction and classification in one step solely from
data. Still, many methods proposed in remote sensing research require implicit
feature extraction through data preprocessing or explicit design of features.

In this work, we compare recent deep learning models on crop type
classification on raw and preprocessed Sentinel 2 data. We concentrate on the
common neural network architectures for time series, i.e., 1D-convolutions,
recurrence, a shallow random forest baseline, and focus on the novel self-
attention architecture. Our central findings are that data preprocessing still
increased the overall classification performance for all models while the
choice of model was less crucial. Self-attention and recurrent neural
networks, by their architecture, outperformed convolutional neural networks on
raw satellite time series. We explore this by a feature importance analysis
based on gradient back-propagation that exploits the differentiable nature of
deep learning models. Further, we qualitatively show how self-attention scores
focus selectively on few classification-relevant observations.

###### keywords:

Self-Attention , Transformer , Time Series Classification, Multitemporal Earth
Observation , Crop Type Mapping , Vegetation Monitoring , Deep Learning

††journal: Journal of Photogrammetry and Remote Sensing

##  1 Introduction

Today, medium resolution Earth observation data is often provided free-of-
charge since the Landsat program, jointly operated by the United States
Geological Survey (USGS) and National Aeronautics and Space Administration
(NASA), set a precedence by adopting an open data policy in 2008 [1; 2]. In
recent years, the Copernicus program of the European Space Agency (ESA) joined
this initiative by a growing family of Sentinel satellites that capture data
from the Earth’s surface at increasing temporal intervals. For instance, the
Sentinel 2 multispectral satellites acquire data at up to
$10\text{\,}\frac{\mathrm{m}}{}$ resolution in $13$ spectral bands every two
to five days. Acquiring surface reflectance data at a high temporal frequency
is a crucial tool for vegetation-related remote sensing applications, such as
crop type mapping [3; 4], or forest monitoring [5; 6; 7]. This data-abundance,
however, leads to computational and organizational challenges and demands for
methods that are problem- and region-agnostic enough to utilize all available
temporal data on large spatial scales. That these _challenges of Big Geo Data_
are not solved yet can be seen on the low exploitation ratio of Sentinel 2
data in 2018. A large body of
$7.76\text{\,}\frac{\mathrm{TiB}}{}$111$$1\text{\,}\frac{\mathrm{TiB}}{}$=2^{40}$\frac{\mathrm{B}}{}$$
Sentinel 2 data was published on a daily average. Still, only
$7.6\text{\,}\frac{\mathrm{\char 37\relax}}{}$ of all published images in
201822262,317,329 published, 4,737,253 downloaded Sentinel 2 products in 2018
have been actually downloaded [8]. Accordingly, $12$ out of $13$ published
images remained unused. Similar figures can be drawn for the Sentinel 1, 3,
and 5 missions. This is a clear sign that, despite claims made by academia and
industry, methods and principles of _big data analytics_ are hardly applied to
their full potential in the field of Earth observation. The reasons for this
low exploitation ratio are manifold. First, visual inspection of satellite
images is still often the first step of data acquisition. Intermediate results
are often required to be visually interpretable to allow for review by domain
experts. Preprocessing steps, like atmospheric correction or manual or
automatic cloud filtering, are common in remote sensing. All of these
processes require computational resources or visual inspection that scales
poorly when applying methods at dense temporal or global spatial scales. For
instance, the commercial provider that preprocessed the datasets used in this
work required multiple months to acquire, atmospherically correct, and
organize the image time series for a comparatively small region of interest.
In contrast, the raw Sentinel 2 time series were downloaded without active
supervision within few days from Google Earth Engine. This demonstrates the
demand for methods that utilize the entire body of available satellite data
and, thus, requires minimal supervision by experts with region-specific expert
knowledge. In principle, deep learning mechanisms seem well-suited to
approximate preprocessing-like mechanisms by jointly learning feature
extraction and classification within one neural network topology using
gradient backpropagation. Developing models that do not strictly require
extensive data preprocessing are likely a key contribution to utilize all
available published satellite data accordingly in the future. In the scope of
this objective, we assess three deep learning mechanisms on four deep learning
models on selectively available preprocessed and readily available raw
satellite data.

Summarizing, the contributions of this work are two-fold:

  1. i)

We first present a large-scale evaluation of deep learning models on
preprocessed and raw satellite data using three mechanisms—i.e., convolution,
recurrence, and self-attention–for crop type identification.

  2. ii)

We further provide quantitative and qualitative analysis of self-attention in
the context commonly used deep learning models for satellite time series
classification.

##  2 Related Work

Since satellite data is inherently multi-temporal, this temporal dimension has
been a topical focus for vegetation-related applications for a long time [4;
3]. The limited availability of multi-temporal satellite data and the
excessive computational demand for temporal stacks of large-scale satellite
images impeded broad exploitation of the temporal dimension. Methods relied on
separate feature extraction, i.e., calculation of vegetation indices [9; 10;
11; 12; 13] with temporal statistics [9; 12; 14], and classification
mechanisms, e.g., with Random Forest Classifiers [12; 15; 11] or Support
Vector Machines [16; 17; 18], as summarized by [19; 20]. Similarly, the
TimeSat software [21; 22] fits non-symmetric Gaussian curves to satellite time
series. The parameters of these curves, i.e., steepest ascent and descent, and
their times of occurrence are assumed to coincide with key phenological
characteristics, such as the onset of greenness or the date of senescence.
This allows for detailed phenological analyses [23; 24] and can be used as
distinctive feature for further classification [25; 26]. The _Continuous
Change Detection and Classification (CCDC)_ [27] first identifies “clear”
observations by a prior cloud filtering with FMask [28; 29] on atmospherically
corrected images. The FMask [28] algorithm itself is a functionally designed
decision tree that exploits the physical properties of Landsat bands to
identify cloud coverage and was later extended to more satellites [29]. The
CCDC algorithm then explicitly models inter- and intra-annual seasonality by a
sum of periodic functions for each pixel fitted with _Robust Iteratively
Reweighted Least Squares (RIRLS)_ [30; 31] to training data. The parameters of
the periodic functions can be used as features for a subsequent
classification. Similarly, the periodic functions can be used to predict
future reluctance values and compared to measured data for change detection.

The _Breaks for Additive Seasonal and Trend (BFAST)_ algorithm [32] similarly
decomposes a satellite time series into piece-wise seasonal and linear
components where the remainder can be used to detect anomalies and the
periodicity can be used as feature for classification. The LandTrendr [33]
aims at removing high-frequency components in a Landsat time series by finding
a combination of successively simpler linear models, defined by vertices, that
still explain the data sufficiently well. A variety of parameters can be
adjusted to tune the model to specific regions and data through the number of
allowed spikes (robustness to clouds), or number of segments (model
complexity). These methods have been designed to require data preprocessing,
i.e., atmospheric correction of cloud filtering, some degree of direct
supervision, i.e., through the choice of parameters, or are taylored towards
specific types of data, e.g., Landsat for CCDC or MODIS for BFAST.

In contrast, Artificial Neural Networks [34; 35] aim at learning feature
extraction and classification with a single dynamic model solely from the
provided data with minimal supervision beyond the provided data. First success
was achieved for the problem of the recognition of handwritten digits [36]
using multi-layer neural networks trained with gradient descent. However,
comparatively recent advances in parallel processing and data availability
[37; 38] were necessary to outperform common feature extraction and
classification pipelines in computer vision [39; 40; 41; 42; 43] and natural
language processing [44; 45; 46] on a large scale. Following these
developments, mono-temporal Earth observation approaches have adapted 2D
convolutional neural networks from the domain of computer vision with great
success [47; 48; 49; 50]. Rotation-invariant feature extraction [51] for scene
classification [52] and metric learning [53] has been similarly tested for
remote sensing image classification. Methods from multi-temporal Earth
observation [54; 55; 56; 57; 58; 59] adapted models from natural language
processing. Here, recurrent neural networks [60], such as Long Short-Term
Memory (LSTM) [61] or Gated Recurrent Units (GRU) [62], were commonly used in
encode-decoder architectures[63] for generative prediction of words. In Earth
observation, the encoder model was utilized for change detection [56; 57] and
land cover [25], as well as for crop type identification [54; 59; 64; 65]. To
utilize both spatial and temporal features from time series, combinations of
convolutional layers and recurrent layers [58; 57] or convolutional-recurrent
networks [66] have been explored and comprehensively compared [64]. These
recurrent neural network encoders can be augmented by soft-attention
mechanisms, originally developed for machine translation [67], as tested in
[59]. Recently, the light-weight time-convolutional feed-forward neural
networks have shown to be a powerful alternative [68] and have achieved
promising results in crop type identification [69]. While the state-of-the-art
in remote sensing remains focused on convolutional and recurrent
architectures, _self-attention_ [44], in combination with pre-training [45;
46], has started dominating the state-of-the-art in natural language
processing. In the remote sensing context, self-attention has hardly been
explored in the canon of existing work [70] on a comparative analysis with
recurrence and convolutional mechanisms.

##  3 Method

In this section, we introduce the notation throughout this work and provide
background on convolution, recurrence and address self-attention before
employing these mechanisms in four neural network topologies in Section 4.

A deep learning model $f_{\boldsymbol{\Theta}}:\mathcal{X}\mapsto\mathcal{Y}$
approximates a mapping from an input domain $\mathcal{X}$ to a target domain
$\mathcal{Y}$. These deep models are implemented by cascaded neural network
layers that form a problem-agnostic, non-linear differentiable function. The
parameters ${\boldsymbol{\Theta}}$ are determined by minimizing an objective
function
$\mathcal{L}(\boldsymbol{y},f_{{\boldsymbol{\Theta}}}(\boldsymbol{X}))$ that
quantifies the dissimilarity of ground truth labels
$\boldsymbol{y}\in\mathcal{Y}$ and the model predictions
$\hat{\boldsymbol{y}}=f_{{\boldsymbol{\Theta}}}(\boldsymbol{X})$. This
objective function is minimized iteratively via mini-batch stochastic gradient
descent

$\displaystyle{\boldsymbol{\Theta}}_{t+1}$
$\displaystyle\leftarrow{\boldsymbol{\Theta}}_{t}-\mu\boldsymbol{s}_{t},$
$\displaystyle\text{ with }\boldsymbol{s}_{t}$
$\displaystyle=\frac{1}{N}\sum_{i=1}^{N}\frac{\partial}{\partial{\boldsymbol{\Theta}}}\mathcal{L}(\boldsymbol{y}_{i},f_{{\boldsymbol{\Theta}}_{t}}(\boldsymbol{x}_{i}))\quad.$
(1)

The gradients $\frac{\partial}{\partial{\boldsymbol{\Theta}}}\mathcal{L}$ are
averaged over a batch of size $N$, while the learning rate $\mu$ determines
the step size. This optimization scheme guarantees that the chosen parameters
${\boldsymbol{\Theta}}$ are optimal for the observed dataset and objective
function.

Deep learning for time series classification aims at learning the mapping
$\hat{\boldsymbol{y}}=f_{\boldsymbol{\Theta}}(\boldsymbol{X})$ from an input
time series $\boldsymbol{X}\in\mathcal{X}^{T}=\mathbb{R}^{T\times
D}=(\boldsymbol{x}_{0},\boldsymbol{x}_{1},\dots,\boldsymbol{x}_{T})$ of
individual measurements $\boldsymbol{x}_{t}\in\mathbb{R}^{D}$ of $D$ features
to one of $C$ classes. These classes are represented by the _one-hot_ target
vector $\boldsymbol{y}\in\\{0,1\\}^{C}$ where $y_{i}\in\\{0,1\\}$ indicates
the boolean class membership. We denote the approximation of the target vector
as $\hat{\boldsymbol{y}}\in[0,1]^{C}$. A neural network is a non-linear and
differentiable function
$f_{{\boldsymbol{\Theta}}}=f^{(1)}_{{\boldsymbol{\Theta}}_{1}}\circ
f^{(2)}_{{\boldsymbol{\Theta}}_{2}}\circ\dots\circ
f^{(L)}_{{\boldsymbol{\Theta}}_{L}}$ of $L$ cascaded layers
$f^{(l)}_{{\boldsymbol{\Theta}}_{l}}$. Each layer
$f^{(l)}_{{\boldsymbol{\Theta}}_{l}}:\boldsymbol{H}^{(l)}\xmapsto{{\boldsymbol{\Theta}}_{l}}\boldsymbol{H}^{(l+1)}$
consists of a linear transformation ${\boldsymbol{\Theta}}_{l}$ mapping an
input representation $\boldsymbol{H}^{l}\in\mathbb{R}^{T\times D_{h}}$ to
corresponding hidden features $\boldsymbol{H}^{l+1}\in\mathbb{R}^{T\times
D_{h}}$ that encode increasingly higher-level features. In time series
classification, one of the following mechanisms are used to implement this
non-linear mapping.

###  3.1 Fully Connected Layers

A dense or fully connected layer
$\boldsymbol{h}_{t}=f^{\text{fc}}_{{\boldsymbol{\Theta}}}(\boldsymbol{x}_{t})\in\mathbb{R}^{D_{h}}$
is applied independently at time instance $t$ of a series. Each vector
$\boldsymbol{x}_{t}\in\mathbb{R}^{D}$ is transformed into
$\boldsymbol{h}_{t}=f^{\text{fc}}_{\\{{\boldsymbol{\Theta}}_{\text{l}},{\boldsymbol{\Theta}}_{\text{b}}\\}}(\boldsymbol{x}_{t})=\phi\left({\boldsymbol{\Theta}}_{\text{l}}^{\\!\top\\!}\boldsymbol{x}_{t}+{\boldsymbol{\theta}}_{\text{b}}\right)$
by a linear transformation
${\boldsymbol{\Theta}}_{\text{l}}\in\mathbb{R}^{D\times D_{h}}$ with optional
bias translation ${\boldsymbol{\theta}}_{\text{b}}\in\mathbb{R}^{D_{h}}$
followed by a non-linear activation function
$\phi\left(\cdot\right)\in\\{\tanh,\sigma,\text{ReLU},\ldots\\}$.

###  3.2 Recurrent Layers

Recurrent Neural Networks [71] extend fully connected layers by contextual
information from previous time steps $t-1$. This results in a transformation
$\boldsymbol{h}_{t}=f^{\text{rnn}}_{{\boldsymbol{\Theta}}}(\boldsymbol{x}_{t},\boldsymbol{h}_{t-1})=\phi\left({\boldsymbol{\Theta}}_{x}^{\\!\top\\!}\boldsymbol{x}_{t}+{\boldsymbol{\Theta}}_{h}^{\\!\top\\!}\boldsymbol{h}_{t-1}+{\boldsymbol{\theta}}_{\text{b}}\right)$.
This formulation is prone to _vanishing_ and _exploding_ gradients through
time [61; 72; 60; 73] that inhibited the extraction of features from long-term
temporal contexts. While the effect of exploding gradients could be controlled
through gradient clipping, vanishing gradients have been addressed by the
introduction of additional gates. This lead to _Long Short-Term Memory (LSTM)_
recurrent networks [61] that introduced four internal gates
$\boldsymbol{f}_{t}=f^{\text{rnn}}_{{\boldsymbol{\Theta}}_{f}}(\boldsymbol{x}_{t},\boldsymbol{h}_{t-1})$,
$\boldsymbol{i}_{t}=f^{\text{rnn}}_{{\boldsymbol{\Theta}}_{i}}(\boldsymbol{x}_{t},\boldsymbol{h}_{t-1})$,
$\boldsymbol{g}_{t}=f^{\text{rnn}}_{{\boldsymbol{\Theta}}_{g}}(\boldsymbol{x}_{t},\boldsymbol{h}_{t-1})$,
$\boldsymbol{o}_{t}=f^{\text{rnn}}_{{\boldsymbol{\Theta}}_{o}}(\boldsymbol{x}_{t},\boldsymbol{h}_{t-1})$.
An additional _cell state_ matrix $\boldsymbol{c}_{t}$ was introduced as a
container for the long-time temporal context. In each iteration, this cell
state
$\boldsymbol{c}_{t}=\boldsymbol{f}_{t}\odot\boldsymbol{c}_{t-1}+\boldsymbol{i}_{t}\odot\boldsymbol{g}_{t}$
is updated by element-wise multiplications $\odot$ with the values of the
forget, input, and modulation gates $\boldsymbol{f}_{t}$,
$\boldsymbol{i}_{t}$. and $\boldsymbol{g}_{t}$, respectively. If the forget
gate evaluates values close to zero, features in the prior cell state
$\boldsymbol{c}_{t-1}$ will be suppressed (forgotten). The modulation and
input gates control the writing of new features from the input to the cell
state $\boldsymbol{c}_{t}$. The cell output
$\boldsymbol{h}_{t}=\boldsymbol{o}_{t}\odot\tanh(\boldsymbol{c}_{t})$ is
calculated using the result from the output gate $\boldsymbol{o}_{t}$ and the
new cell state $\boldsymbol{c}_{t}$. Overall, these transformations yield the
long short-term memory update
$f^{\text{lstm}}_{{\boldsymbol{\Theta}}}:\boldsymbol{h}_{t},\boldsymbol{c}_{t}\leftarrow\boldsymbol{x}_{t},\boldsymbol{h}_{t-1},\boldsymbol{c}_{t-1}$.
An alternative to LSTMs are _Gated Recurrent Units (GRUs)_ [62] that follow
the same principle of additional gates, but are parameterized with fewer gates
and less weights. Despite computational benefits, empirical evaluations did
not show a significant difference in performance between these
parameterizations [72].

###  3.3 1D Convolutional Layers

Convolutional layers extract features by correlating the input signal
$\boldsymbol{X}\in\mathbb{R}^{T\times D}$ with a set of $D_{\text{h}}$
convolutional kernels
${\boldsymbol{\Theta}}=({\boldsymbol{\Theta}}_{0},{\boldsymbol{\Theta}}_{1},\dots,{\boldsymbol{\Theta}}_{D_{\text{h}}})$
with ${\boldsymbol{\Theta}}_{d}\in\mathbb{R}^{K\times D}$ employing a
convolutional operation
$\boldsymbol{H}=\phi(\boldsymbol{X}\ast{\boldsymbol{\Theta}})$ followed by an
element-wise nonlinear activation function $\phi$. The size of the
convolutional kernel, $K$, determines the receptive field of each layer. In
contrast to recurrent layers, convolutional layers extract features from a
fixed temporal neighborhood. The receptive field increases through the number
of layers. Convolutional layers perform well in the field of computer vision
where features from a local neighborhood of an image are extracted with small
local kernels. In a temporal context, larger kernel sizes are used. There, the
learned kernel resembles a correlation between learned patterns in the time
domain [74; 75] with recent applications in the field of remote sensing [69].

###  3.4 Self-Attention Layers

${\boldsymbol{A}^{\\!\top\\!}}$$\boldsymbol{H}$$\boldsymbol{V}$$\boldsymbol{K}$$\boldsymbol{Q}^{\\!\top\\!}$$\phi\left({\boldsymbol{\Theta}}_{V}^{\\!\top\\!}\boldsymbol{X}\right)$$\phi\left({\boldsymbol{\Theta}}_{K}^{\\!\top\\!}\boldsymbol{X}\right)$$\phi\left({\boldsymbol{\Theta}}_{Q}^{\\!\top\\!}\boldsymbol{X}\right)$
(a) visualization of the matrix multiplications involved in the calculation of
attention scores $\boldsymbol{A}$

$\boldsymbol{V}:$$\boldsymbol{H}:$$\boldsymbol{A}:$ (b) relationship between
values and outputs drawn as bipartite graph

Figure 1: Schematic self-attention operation using
${T_{\text{in}}}={T_{\text{out}}}=4$, ${D_{\text{in}}}={D_{\text{out}}}=2$ and
${D_{\text{h}}}=4$. The hidden output $\boldsymbol{H}\in\mathbb{R}^{4\times
2}$ results from matrix multiplication with attention scores
$\boldsymbol{A}\in[0,1]^{4\times 4}$ with a time series $\boldsymbol{V}\in
R^{4\times 2}$. These attention scores themselves are calculated by a matrix
multiplication between key $\boldsymbol{K}\in\mathbb{R}^{4\times 4}$ and query
matrices $\boldsymbol{Q}\in\mathbb{R}^{4\times 4}$. In self-attention, all
key, query and value matrices originate from an input time series
$\boldsymbol{X}$ transformed via a linear mapping
$\boldsymbol{X}{\boldsymbol{\Theta}}$ with a non-linear activation function
$\phi(\cdot)$

Attention mechanisms [67] allow a neural network to extract features from
observations at specific time steps of an input time series of _values_
$\boldsymbol{V}=(\boldsymbol{v}_{0},\boldsymbol{v}_{1},\dots,\boldsymbol{v}_{T_{\text{in}}})\in\mathbb{R}^{{T_{\text{in}}}\times{D_{\text{in}}}}$.
This is realized by a weighted sum
$\boldsymbol{h}_{i}=\sum_{t=0}^{T_{\text{in}}}\alpha_{t}\boldsymbol{v}_{t}=\boldsymbol{\alpha}^{\\!\top\\!}\boldsymbol{V}$
of attention scores
$\boldsymbol{\alpha}\in[0,1]^{T_{\text{in}}},\sum_{i}\alpha_{i}=1$. Extending
this to a matrix multiplication

$\displaystyle\boldsymbol{H}=\boldsymbol{A}^{\\!\top\\!}\boldsymbol{V}$ (2)

with attention matrix
$\boldsymbol{A}=(\boldsymbol{a}_{0},\dots,\boldsymbol{a}_{T_{\text{out}}})$
and
$\boldsymbol{H}=(\boldsymbol{h}_{0},\dots,\boldsymbol{h}_{T_{\text{out}}})\in\mathbb{R}^{{T_{\text{out}}}\times{D_{\text{out}}}}$,
generates ${T_{\text{out}}}$ output vectors in parallel. To illustrate these
operations and to prepare for numerical results later in Section 8.4, we show
a visual representation of these operations in Fig. 1. Attention scores
themselves are the result of a second matrix multiplication
$\boldsymbol{A}=\boldsymbol{Q}\boldsymbol{K}^{\\!\top\\!}$ in Eq. 3 of _key_
and _query_ matrices
$\boldsymbol{K}\in\mathbb{R}^{{T_{\text{in}}}\times{D_{\text{h}}}}$ and
$\boldsymbol{Q}\in\mathbb{R}^{{T_{\text{out}}}\times{D_{\text{h}}}}$,
respectively. This results in the generic formulation of attention

$\displaystyle\boldsymbol{H}=\underbrace{\text{softmax}\left({\boldsymbol{Q}\boldsymbol{K}^{\\!\top\\!}}\right)^{\\!\top\\!}}_{\boldsymbol{A}^{\\!\top\\!}}\boldsymbol{V}$
(3)

that involves a softmax function to normalize $\boldsymbol{A}$. In self-
attention [44]

$\displaystyle\boldsymbol{H}=f_{\\{{\boldsymbol{\Theta}}_{K},{\boldsymbol{\Theta}}_{Q},{\boldsymbol{\Theta}}_{V}\\}}(\boldsymbol{X})=\text{softmax}\left(\frac{{\phi\left({\boldsymbol{\Theta}}_{Q}^{\\!\top\\!}\boldsymbol{X}\right)}{\phi\left({\boldsymbol{\Theta}}_{K}^{\\!\top\\!}\boldsymbol{X}\right)^{\\!\top\\!}}}{\sqrt{{D_{\text{h}}}}}\right){\phi\left({\boldsymbol{\Theta}}_{V}^{\\!\top\\!}\boldsymbol{X}\right)},$
(4)

this formulation is adopted with unified dimensions
${D_{\text{in}}}={D_{\text{out}}}={D_{\text{h}}}$ and
${T_{\text{in}}}={T_{\text{out}}}$ and where the keys
$\boldsymbol{K}=\phi\left({\boldsymbol{\Theta}}_{K}^{\\!\top\\!}\boldsymbol{X}\right)$,
queries
$\boldsymbol{Q}=\phi\left({\boldsymbol{\Theta}}_{Q}^{\\!\top\\!}\boldsymbol{X}\right)$,
and values
$\boldsymbol{V}=\phi\left({\boldsymbol{\Theta}}_{V}^{\\!\top\\!}\boldsymbol{X}\right)$
matrices originate from the same input matrix $\boldsymbol{X}$ transformed by
three linear transformations with non-linear activation function
$\phi\left(\cdot\right)$. Following the original work [44], the attention
scores are scaled by $\sqrt{{D_{\text{h}}}}$ for better gradient
backpropagation. Today, such self-attention mechanisms are the core of state-
of-the-art natural language processing models [45; 46] and will be
experimentally analyzed later in Sections 8.4 and 8.6.

###  3.5 Soft-Attention

Soft-attention [76], as implemented in the DuPLO model [59], is a special case
of Eq. 3 where values $\boldsymbol{V}=\boldsymbol{X}$ and queries
$\boldsymbol{Q}=\tan\left(\boldsymbol{\Theta}_{\alpha}\boldsymbol{X}\right)$
are obtained from a linear transformation from the input tensor
$\boldsymbol{X}$. In contrast to self-attention, the keys
$\boldsymbol{K}=\boldsymbol{\Theta}_{K}$ are realized by a fixed weight matrix
which is learned by gradient descent. Following the definitions of [59] and
modified from Eq. 4, soft-attention yields

$\displaystyle\boldsymbol{h}=f_{\\{{\boldsymbol{\Theta}}_{K},{\boldsymbol{\Theta}}_{A}\\}}(\boldsymbol{X})=\underbrace{\text{softmax}\left(\tan\left({\boldsymbol{\Theta}}_{A}^{\\!\top\\!}\boldsymbol{X}\right){\boldsymbol{\Theta}}_{K}\right)^{\\!\top\\!}}_{\boldsymbol{\alpha}^{\\!\top\\!}}\boldsymbol{X}$
(5)

where the non-linearity $\phi$ is implemented as tangens $\tan$ without
scaling factor $D_{h}$. We compare self-attention with the simplified soft-
attention in a dedicated experiment in Section 8.5.

##  4 Models

The previous section provided an overview on neural network layers used for
temporal feature extraction and introduced temporal convolution, recurrence,
and self-attention. In this section, we describe four neural network
topologies, each of them using one of these layer mechanisms, which will be
evaluated experimentally in the following333The models are available at
https://github.com/marccoru/crop-type-mapping/tree/master/src/models.

###  4.1 Recurrent Long Short-Term Memory Neural Network

$\boldsymbol{X}$$(\vec{\boldsymbol{c}},\vec{\boldsymbol{h}})$$(\reflectbox{$\vec{\reflectbox{$\boldsymbol{c}$}}$},\reflectbox{$\vec{\reflectbox{$\boldsymbol{h}$}}$})$$(\vec{\boldsymbol{c}},\vec{\boldsymbol{h}})$$(\reflectbox{$\vec{\reflectbox{$\boldsymbol{c}$}}$},\reflectbox{$\vec{\reflectbox{$\boldsymbol{h}$}}$})$$\boldsymbol{c}=\text{concat}\left(\vec{\boldsymbol{c}},\reflectbox{$\vec{\reflectbox{$\boldsymbol{c}$}}$}\right)$softmax$({\boldsymbol{\Theta}}^{\\!\top\\!}\boldsymbol{c})$$\boldsymbol{y}$
(a) LSTM-RNN (recurrence)

$\boldsymbol{X}$ self-attention
feed forward self-attention
feed forward max-
poolingsoftmax$({\boldsymbol{\Theta}}^{\\!\top\\!}\boldsymbol{H})$$\boldsymbol{y}$
(b) Transformer (self-attention)

$\boldsymbol{X}$BN+ReLU$+$BN+ReLU$+$BN+ReLU$+$BN+ReLU$+$BN+ReLU$+$BN+ReLU$+$BN+ReLU$+$BN+ReLU$+$BN+ReLU$+$avg-
poolingsoftmax$({\boldsymbol{\Theta}}^{\\!\top\\!}\boldsymbol{H})$$\boldsymbol{y}$
(c) MS-ResNet (convolution)

$\boldsymbol{X}$BN+ReLUBN+ReLUBN+ReLUflattensoftmax$({\boldsymbol{\Theta}}^{\\!\top\\!}\boldsymbol{H})$$\boldsymbol{y}$
(d) TempCNN (convolution)

Figure 2: A schematic overview over the four network topologies used in this
study.

We designed a multi-layer _Recurrent Neural Network (RNN)_ [71] with _Long
Short-Term Memory (LSTM)_[61] cell architecture, similar to our previous work
[54]. This network encodes a satellite time series to increasingly higher-
level $D$-dimensional representations through $N$ cascaded bidirectional LSTM
layers. The hidden states are initialized with zero-valued vectors. The last
hidden states from forward and backward passes, $\vec{\boldsymbol{c}}$ and
$\vec{\reflectbox{$\boldsymbol{c}$}}$ , are concatenated and produce final
classification scores via a dense output layer with softmax activation
function. This network architecture is illustrated in Fig. 2(a) We tuned the
conditionality of hidden states
$D_{h}\in\left\\{2^{4},2^{5},\dots,2^{8}\right\\}$, the number of recurrent
layers $L\in\left\\{1,2,\dots,7\right\\}$, and applied dropout
$p_{\text{drop}}\in[0,1]$ between the recurrent layers at training time.

###  4.2 Transformer

For the self-attention model, we adopted the encoder architecture of the self-
attention Transformer network [44], as illustrated in Fig. 2(b). We discarded
the step of word embedding, as the satellite time series data lives in a
continuous space of spectral reflectance values. Following Vaswani ifnextchar.
et al ifnextchar, et al. et al. [44], we added a positional encoding to the
time series, since self-attention can not utilize the sequential correlation
of time series. Subsequently, the time series with positional encoding is
transformed into higher-level $D$-dimensional feature representations through
$L$ transformer blocks. Each block encodes features through a multi-head self-
attention mechanism followed by multiple dense layers applied to each time
instance independently. Skip connections are introduced between the layers and
layer normalization is used throughout the model. We refer to [44] for a
detailed description of the layer topology. The representation of the last
layer $\boldsymbol{H}\in\mathbb{R}^{H\times T}$ is then reduced to
$\boldsymbol{h}\in\mathbb{R}^{D}$ by a global maximum pooling through the time
dimension. This reduced representation is then projected to scores for each
class by a final fully-connected layer with a softmax activation function.

For this architecture, we tuned the dimensionality of the hidden states
$D_{h}\in\\{2^{4},2^{5},\dots,2^{8}\\}$, the number of self-attention layers
$L\in\\{1,2,\dots,8\\}$, and the number self attention heads
$H\in\\{1,2,\dots,8\\}$ that determine the number of self-attention mechanisms
applied in parallel.

###  4.3 Multi-Scale Residual Networks

_Residual Networks (ResNets)_ [41] are a common backbone architecture in
computer vision. They consist of residual blocks
$f^{\text{cbn}}_{\boldsymbol{\Theta}}=f^{\text{c}}_{\boldsymbol{\Theta}}\circ
f^{\text{b}}\circ f^{\text{r}}$ of convolution layers
$f_{\boldsymbol{\Theta}}^{\text{c}}$ and batch normalization layers
$f^{\text{b}}$ [77], followed by a _Rectified Linear Unit (ReLU)_ activation
function $f^{\text{r}}$. Characteristically, in each block
$\boldsymbol{H}=f^{\text{cbn}}_{\boldsymbol{\Theta}}(\boldsymbol{X})+\boldsymbol{X}$,
the encoded feature representation
$f^{\text{cbn}}_{\boldsymbol{\Theta}}(\boldsymbol{X})$ is added to the
identity mapping of the input. This forms residual skip connections where each
residual layer numerically adds higher-level features to the forward
propagated representation. These skip connections aid gradient backpropagation
through the network and allow the training of very deep models.

These convolutional residual connections have been adapted to time series
classification [78] where 2D convolutions through the spatial dimensions are
replaced by 1D convolutions through time. Since convolutional layers, by
design, extract features from the local (temporal) neighborhood, architectures
that process time series at multiple scales [79] have shown good results on
time series benchmark datasets [80; 81]. In this work, we utilize the
architecture proposed by Wang ifnextchar. et al ifnextchar, et al. et al.
[78]444Original implementation available at https://github.com/geekfeiw/Multi-
Scale-1D-ResNet where an input time series
$\boldsymbol{X}\in\mathbb{R}^{(T=512)\times D_{\text{in}}}$ is processed in
three separate streams of each three residual blocks with increasing
convolutional kernel sizes $K^{\text{c}}\in\\{3,5,7\\}$. Hence, each stream
extracts features at a different temporal scale, as drawn in Fig. 2(c). The
representation after each stream is then average-pooled with pooling kernel
sizes $K^{\text{p}}\in\\{16,11,6\\}$, yielding three feature representations
$\boldsymbol{h}^{K^{c}=3},\boldsymbol{h}^{K^{c}=5},\boldsymbol{h}^{K^{c}=7}\in\mathbb{R}^{D=256}$,
and concatenated to a common representation
$\boldsymbol{H}\in\mathbb{R}^{D=768}$. A final fully connected layer with
softmax activation function reduces this vector to individual activation
scores per class. Note that the design of this network requires a fixed
sequence length of $T=512$. To obtain this sequence length, we interpolated
the satellite time series by the nearest neighbor method. Since the
architecture is precisely-defined in the original implementation, we solely
tune the number of convolutional kernels
$D_{h}\in\\{2^{4},2^{5},\dots,2^{8}\\}$ for this architecture.

###  4.4 TempCNN

Recently, the temporal convolutional neural network (TempCNN) [69] has been
proposed and evaluated specifically for the purpose of crop type mapping. It
is a comparatively light-weight architecture of three sequential 1D
convolutional layers followed by batch normalization, a ReLU activation
function, and dropout. The encoded features are flattened and passed to a
final fully connected layer with batch normalization, ReLU activation
function, and dropout, as shown in Fig. 2(d). These features are then
projected to scores per class using a final fully connected layer and softmax.

The tuneable hyperparameters of this model are the number of convolutional
kernels $D_{h}\in\left\\{2^{4},2^{5},\dots,2^{8}\right\\}$ that determine the
dimensionality of hidden states, the respective kernel sizes
$K\in\\{3,5,7\\}$, and the dropout probability $p_{\text{drop}}\in[0,1]$.

###  4.5 DuPLO

The _DUal view Point deep Learning architecture for time series classificatiOn
(DuPLO)_ [59] model is a complex deep learning model designed for crop type
classification from sequences of small satellite images of five by five
pixels. Since we do not use sequences of images in this work, we treated the
time series as a single-pixel images. It consists of two streams. A three-
layer CNN stream uses 2d convolutions to aggregated spatial features
independently of time. The second stream implements a 2D-CNN encoder and mono-
directional RNN layer implemented by a Gated Linear Unit (GRU) [62] for
temporal characteristics. The recurrent neural network features
$\boldsymbol{H}\in\mathbb{R}^{H\times
T}\mapsto\boldsymbol{y}\in\mathbb{R}^{H}$ are then aggregated via a simplified
soft-attention mechanism [76], described in Section 3.5. According to the
authors, this model requires a special training procedure through additional
weighted auxiliary losses from additional dense layers for each stream to
simplify gradient propagation through.

###  4.6 Random Forests

As shallow learning baseline, we tested a _Random Forest (RF)_ model using the
scikit-learn framework [82]. We augmented the raw reflectance input features
by additional spectral indices, since feature extraction is usually required
for a good performance. Hence, we added the features _Normalized Difference
Vegetation Index (NDVI)_ , _Normalized Difference Water Index (NDWI)_ ,
_Bridghtness Index (BI)_ , _Inverted Red-Edge Chlorophyll Index (IRECI)_ , and
_Enhanced Vegetation Index (EVI)_ to the time series data. Overall, the random
forest model receives a flattened vector of all spectral values for each time.
This results in a raw feature vector
$\boldsymbol{x}_{\text{raw}}\in\mathbb{R}^{1260}$ (70 times each 13 spectral
bands and 5 spectral indices) and
$\boldsymbol{x}_{\text{pre}}\in\mathbb{R}^{345}$ (23 times each 10 spectral
bands and 5 indices). Also note, that we did not obtain satisfactory results
on the random forest using the entire body of the dataset likely due to class
imbalance of the labels. Hence, we chose to sample a more uniform subset of up
to 500 time series per class to train and test the random forest model. We
argue that this would rather advantage random forest in the classification
results (especially for the macro accuracies).

We tuned the number of trees $\left\\{200,400,\dots,2000\right\\}$, the number
of features to be considered at every split
$\left\\{\text{auto},\text{sqrt}\right\\}$, the maximum depth the the trees
$\left\\{10,20,\dots,110\right\\}$, the minimum number of samples required to
split a node $\left\\{2,3,\dots,10\right\\}$, the minimum number of samples
required at each leaf node $\left\\{1,2,3,4\right\\}$, and whether or not to
use bootstrapping $\left\\{\text{true},\text{false}\right\\}$.

##  5 Training Details

We trained the LSTM-RNN , MS-ResNet , and TempCNN models via mini-batch
stochastic gradient descent with gradients scaled using the Adam optimizer
[83] with momentum parameters $\beta_{1}=0.9$, and $\beta_{2}=0.98$. The
learning rate and weight decay parameters were sampled from log-uniform
distributions over $[${10}^{-8}$,${10}^{-1}$]$ and
$[${10}^{-12}$,${10}^{-1}$]$, respectively, and determined via hyperparameter
tuning, as described in Section 6. Following Vaswani ifnextchar. et al
ifnextchar, et al. et al. [44], we employed a learning rate scheduler for the
Transformer model where the learning rate is first increased linearly for
$N_{\text{warmup}}\in\\{${10}^{1}$,${10}^{2}$,${10}^{3}$\\}$ gradient steps
and then decayed exponentially. We stopped training early when the loss did
not decrease over an average of the last 10 epochs. This condition must have
been true for five epochs in a row to stop the training process.

##  6 Model Selection

The model architectures and training procedures depend on a variety of
hyperparameters that may vary based on the objective and the dataset. We
described the tuneable hyperparameters for each model in the previous
sections. These were, for the deep learning models, the dimensionality of the
hidden vectors, the number of layers, the number of self-attention heads, the
convolutional kernel sizes, the dropout probability, the number of transformer
warmup steps, the learning rate, and the weight decay rate. To determine the
optimal set of parameters, we sampled candidate hyperparameters from a
hyperparameter space, trained on a subset of the training partition, and
evaluated the performance on a validation set. The partitioning scheme and
spatial separation of the datasets are described in Section 7. Here, we first
sampled hyperparameters from the search space at random for the first $34$
models and recorded the validation performance measured in terms of the kappa
metric [84]. Based on these initial points, a kernel density estimate scaled
by the validation performance of the model determined a probability
distribution that was used to sample the next best set of hyperparameters,
following the protocol suggested by Bergstra ifnextchar. et al ifnextchar, et
al. et al. [85] and implemented using their HyperOpt framework. We combined
this with an _asynchronous successive halving_ strategy [86]. This algorithm
splits the $60$ training epochs into four brackets with a grace period of $10$
epochs. After each bracket, only the best-performing half of the models
continued training, while the other runs were stopped at the end of each
bracket. Both of these were implemented using the Ray-Tune framework [87] that
allowes hyperparameter training of four models in parallel per GPU on a
Nvidia-DGX1. For each model architecture, we evaluated on $300$ sets of
hyperparameters for both preprocessed and raw datasets having resulted in
$2400$ model evaluations.

The hyperparameters of the RF models were determined through a random search
on $300$ runs with three-fold cross-validation. To reduce tuning time and the
class imbalance in the data, we used only up to $500$ samples per class. The
architecture of the DuPLO model was precisely-defined in [59] with detailed
specifications for hidden dimensions, learning procedure and learning rate.
Hence, we employed the model directly as specified.

##  7 Data

(a) The three areas of interest in Bavaria, Hollfeld, Krumbach and Bavarian
Forest are located in different regional environments indicated by the
elevation map. We show the field parcels colored by crop type, along with the
random train test split for each of the regions.

(b) Class frequencies per region on a 23 class partition

(c) Class frequencies per region on a 12 class partition

Figure 3: Test regions for labels and satellite time series

We evaluated the mechanisms described in Section 3 that were implemented in
the models of Section 4 on the task of crop type identification. Crop type
identification is a sub task of land cover and land use classification where
the model has to extract classification-relevant features and learn a
discriminative decision function to separate the classes of vegetation. We
analyze the extracted features later in Section 8.6. Vegetation life cycle
events, known as phenology, provides a distinctive temporal signal to identify
types of vegetation using a limited set of spectral channels. This makes the
temporal signal a key source of relevant features when learning a
discriminative model to differentiate various types of vegetation and thus
well-suited to test the mechanisms and models of this work.

We focus on three spatially separate regions in Bavaria, as shown in Fig. 3.
These regions are _Hollfeld_ in Upper Franconia, _Krumbach_ in Swabia, and a
northeastern portion of the _Bavarian Forest_. These regions are located in
approximately $100\text{\,}\frac{\mathrm{km}}{}$ distance from each other.
While the climate conditions are comparatively similar, different elevations
and soil conditions favor differences in the distribution of cultivated crops,
as can be observed from the class distribution histograms in Figs. 3(b) and
3(c).

The label data for this study originates from a joint cooperstion project with
the _Bavarian State Ministry of Food, Agriculture and Forestry (StMELF)_ and
the German remote sensing company _GAF AG_. This enabled us to obtain two
Sentinel 2 time series datasets from the same field parcels. These are one
dataset with raw top-of-atmosphere Sentinel 2 observations acquired with
minimal effort from Google Earth Engine, as described in Section 7.2.1, and
one preprocessed time series dataset provided by GAF AG which can be
considered prototypical for the industry-standard, as outlined in Section
7.2.2.

To obtain spatially separate partitions [88] for _training_ of model
parameters using gradient descent, _validation_ of hyperaparameter sets and
final _evaluation_ of the model, we divided the three regions further into
rectangular blocks of $4500\text{\times}4500\text{\,}\frac{\mathrm{m}}{}$ with
a $500\text{\,}\frac{\mathrm{m}}{}$ margin between blocks, as shown in the
train-test split in Fig. 3(a). These blocks were randomly assigned to
training, validation, and evaluation partitions obeying a 4:1:1 ratio which
resulted in 94, 21, and 29 blocks respectively. We decided for such a block-
wise spatial separation in order to enforce independence of the dataset
partitions without implicit overfitting, as experimentally evaluated and
observed in previous work [54] and as recommended for geospatial data by
further studies [89; 88] focusing the implicit bias of spatial auto-
correlation.

###  7.1 Crop Type Labels

The _Common Agricultural Policy_ of the European Union subsidizes farmers
based on the type of crops they cultivated. Each member country is required to
gather geographical information of the geometry of the parcel and the type of
crop. This information is provided by obligatory surveys, as part of the
subsidy application process, directly from the farmers. National agencies
monitor the correctness either by gathering control samples in-situ or by
utilizing remote sensing and Earth observation technology.

The crop label categories provided for this study were provided by StMELF.
They follow a _long-tailed class distribution_ with more than $269$ distinct
categories. Here, the most common $15$, $26$, and $62$ categories cover
$90\text{\,}\frac{\mathrm{\char 37\relax}}{}$, $95\text{\,}\frac{\mathrm{\char
37\relax}}{}$, and $99\text{\,}\frac{\mathrm{\char 37\relax}}{}$ of the field
parcels, respectively. In cooperation with StMELF and GAF, a set of land-use
and land-cover categories was aggregated and selected with respect to the aims
and objectives of the ministry. From this aggregation, we selected two labeled
datasets: The first one contains $23$ classes, as shown in Fig. 3(b), and
resembles the land-use of the parcels. These categories cover, for instance,
multiple types of grassland. This dataset is a challenging to classify, since
multiple categories (e.g., grassland for cattle-use and for machining) share
similar surface reflectance features measured by the satellite. We also
aggregated these categories further into a second dataset that focuses on $12$
land cover categories, as shown in Fig. 3(c). By evaluating models on two
23-class land-use categorization and 12-class land cover categorization
datasets, we aimed at reporting model accuracies from two objectives of
different difficulty.

###  7.2 Satellite Data

(a) raw time series of a meadow parcel

(b) preprocessed time series of a meadow parcel

Figure 4: An illustration and comparison of a raw and a preprocessed Sentinel
2 time series of the same meadow field parcel. Preprocessing allows for a
visual interpretation. The onset of growth after timestep $t=5$ is clearly
visible. Also, several cutting events can be observed over the time of the
vegetation period. The preprocessed time series, however, contains repeated
values due to temporal interpolation and cloud removal. In the raw time
series, most information from the measured signal is retained. Noise caused
by, e.g., atmospheric effects and clouds, obscures the phenological events.

We utilized data from the optical Sentinel 2 satellite constellation that
consists of two satellites orbiting the Earth in a sun-synchronous orbit on
opposite tracks. These satellites observe the same spot on the Earth’s surface
every 2 to 5 days, depending on the latitude. The data is acquired by a line-
scanner capturing 13 spectral bands, ranging from ultra-violet wavelengths,
for observing atmospheric water vapor, over optical and near-infrared
wavelengths, sensitive to chlorophyll and photosynthesis, up to short-wave
infrared wavelengths, that are sensitive to soil moisture. In regions where
the sensor strips overlap, we observed approximately $140$ measurements of the
same point during the entire year of 2018, while on the stripe centers we only
recorded $70$ observations within the same time range.

####  7.2.1 Raw Dataset

For the raw Sentinel 2 dataset, we utilized the _top-of-atmosphere_
reflectance values at the processing level L1C. This data was acquired from
_Google Earth Engine_ [90] and queried for each field parcel individually.
Pixels located within the boundaries of a field parcel were mean-aggregated
into a single feature vector of $13$ spectral bands at each time. We show
examples of a meadow parcel of the raw Sentinel 2 time series in Fig. 4(a).
Note that cloud coverage, visible as positive peaks within the reflectance
profiles, dominates the signal and makes a visual interpretation of this time
series difficult. This time series dataset is challenging to classify but can
be acquired at minimal effort. Hence, we benchmarked our models on this harder
objective.

####  7.2.2 Preprocessed Dataset

For further evaluation, we had access to a second dataset that originates from
the same publicly available top-of-atmosphere data products. In contrast to
the raw data set, this data was processed by GAF through their tested and
operational preprocessing engine. This process includes common preprocessing
techniques, such as, e.g., atmospheric correction, temporal selection of
cloud-free observations, a focus on observations of the vegetative period, and
cloud masking. We show an example of the preprocessed dataset in Fig. 4(b)
compared to the identical parcel in Fig. 4(a). Here, the cloudy observations
have been identified and filtered by a separate cloud classification model.
This makes this dataset easier to classify by shallow models, as distinctive
phenological features, i.e., onset of growth and cutting of the meadows, can
be visually distinguished. Overall, this preprocessing pipeline can be
considered prototypical for an industry standard but requires significant
computational and design effort to generate.

##  8 Experiments and Results

In this section, we experimentally compare the transformer model, that is
based on self-attention, to a recurrent neural network and two convolutional
neural networks, as described in Section 4. This experimental section is
structured in three parts: First in Section 8.1, we show quantitative results
on all evaluated neural network architectures. Here, we present results on
preprocessed and raw Sentinel 2 data described in Section 7.2 and on two sets
of categories. One 23-class categorization evaluates model performance on land
use classification while the other 12-class categorization focuses on land
cover. Next in Section 8.3, we analyze the ability of the models to suppress
noise, e.g., induced by clouds, in raw time series data by a feature
importance analysis based on gradient backpropagation. Finally in Sections 8.4
and 8.6, we focus specifically on the self-attention mechanism and analyze
activation scores and internal states of the transformer model in detail.

###  8.1 Quantitative Model Evaluation

Table 1: Comparison of models on preprocessed (pre) and raw dataset on the
23-class land use and the 12-class land cover categorization. The values
reported are the mean and standard deviation of three models with the best,
second and third-best hyperparameter sets trained on the training and
validation partitions and tested on the evaluation partition.

$\kappa$ | RF | LSTM-RNN | Transformer | DuPLO | MS-ResNet | TempCNN
---|---|---|---|---|---|---
pre | 0.76 |  0.78±0.01 |  0.79±0.03 | 0.81 |  0.76±0.03 |  0.80±0.00
raw | 0.53 | 0.71±0.01 | 0.71±0.03 | 0.68 |  0.69±0.03 |  0.67±0.02
(a) Kappa metric 23-class dataset

$\kappa$ | RF | LSTM-RNN | Transformer | DuPLO | MS-ResNet | TempCNN
---|---|---|---|---|---|---
pre | 0.86 |  0.87±0.01 |  0.87±0.04 | 0.89 |  0.85±0.01 |  0.88±0.03
raw | 0.62 | 0.83±0.01 |  0.82±0.02 | 0.78 |  0.78±0.02 |  0.71±0.07
(b) Kappa metric 12-class dataset

acc. | RF | LSTM-RNN | Transformer | DuPLO | MS-ResNet | TempCNN
---|---|---|---|---|---|---
pre | 0.83 |  0.85±0.01 |  0.85±0.02 | 0.86 |  0.83±0.02 | 0.86±0.00
raw | 0.71 | 0.81±0.01 |  0.80±0.02 | 0.79 |  0.79±0.03 |  0.79±0.00
(c) Overall accuracy metric 23-class dataset

acc. | RF | LSTM-RNN | Transformer | DuPLO | MS-ResNet | TempCNN
---|---|---|---|---|---|---
pre | 0.91 | 0.92±0.01 | 0.92±0.03 | 0.92 |  0.91±0.01 | 0.92±0.02
raw | 0.80 | 0.90±0.00 |  0.89±0.01 | 0.87 |  0.87±0.01 |  0.83±0.04
(d) Overall accuracy metric 12-class dataset

f1 | RF | LSTM-RNN | Transformer | DuPLO | MS-ResNet | TempCNN
---|---|---|---|---|---|---
pre | 0.38 |  0.47±0.02 |  0.50±0.06 | 0.55 |  0.49±0.03 |  0.50±0.02
raw | 0.18 |  0.43±0.01 | 0.45±0.06 | 0.34 |  0.44±0.01 |  0.36±0.01
(e) Class-mean $f_{1}$ score 23-class dataset

f1 | RF | LSTM-RNN | Transformer | DuPLO | MS-ResNet | TempCNN
---|---|---|---|---|---|---
pre | 0.55 |  0.60±0.02 | 0.66±0.07 | 0.66 |  0.64±0.02 |  0.60±0.06
raw | 0.34 |  0.63±0.01 | 0.64±0.07 | 0.51 |  0.55±0.04 |  0.41±0.07
(f) Class-mean $f_{1}$ score 12-class dataset

We compared the performance of the deep learning models LSTM-RNN [61],
Transformer [44], MS-ResNet [78], the DuPLO model [59], and TempCNN [69], as
well as a Random Forest (RF) classifier as shallow baseline. We determined the
optimal hyper-parameters separately for preprocessed and raw Sentinel 2 time
series datasets, and for the 23 and 12 class categorizations, as described in
Section 6. For each experiment, we trained and evaluated three different
models with the best, second-best, and third-best hyperparameter configuration
and random seeds for parameter initialization and composition of the training
batches. In Table 1, we present the mean and standard deviation of the
accuracy metrics from these three results. All models were trained and
evaluated on block partitions for training and evaluation (cf. Fig. 3) in all
three regions Hollfeld, Krumbach and Bavarian forest. We evaluated these
models on raw and preprocessed satellite time series on identical field
parcels.

It is noteworthy that the overall accuracy measure, as presented in Tables
1(c) and 1(d), over-represents frequent classes. Since the datasets used for
this study show a heavily imbalanced class distribution, the accuracy of
frequent classes dominates this metric. Nevertheless, it is an intuitive
measure and a good representation of the example-wise accuracy, representative
of the visual impression from observing a spatial map classification. To
account for the less frequent classes, we also report Cohen’s kappa metric
[84] in Tables 1(a) and 1(b). This is a correlation score that is frequently
used in remote sensing and normalizes the classification scores by the
probability of a random chance prediction based on empirical class
frequencies. As a further measure of performance, we report the $f_{1}$ score,
i.e., the harmonic mean of precision and recall, for each class and average
over all classes in Tables 1(e) and 1(f). By doing so, all classes get equally
weighted disregarding the the number of samples per category. The class-mean
$f_{1}$ scores are generally lower than accuracy and kappa, as we chose a
large set of classes where some classes semantically overlap or only have few
examples which make it difficult for a data-driven neural network to learn
feature extraction and decision function. Overall, we aimed at comparing
different properties of the classification models. Hence, the overall accuracy
reflects the classification accuracy per field parcel, while the $f_{1}$ score
measures the accurate classification of all classes.

Table 1 reveals that data preprocessing had a positive effect on the accuracy,
$f_{1}$ score, and kappa of all models. It seems that the manual supervision
during preprocessing improved the classification performance. This better
performance on pre-processed time series is, to a certain degree, expected,
since this preprocessing makes a visual identification of classification
relevant events possible, as could be seen in Fig. 4(b). Comparing the visual
examples at Fig. 4 suggests that the classification of the preprocessed
datasets is an easier task than approximating preprocessing-like mechanism
jointly together with classification in a holistic end-to-end model. Overall,
the region- and domain-specific expertise needed to design such preprocessing
pipelines increases the classification performance compared to classification
of the raw datasets were the deep learning models have to learn preprocessing-
like mechanisms solely based on the provided examples without access to model
knowledge. This holds especially true if the number of samples is limited,
some examples are falsely labeled, or the semantic representation of classes
do overlap. Considering this, it is notable that the LSTM-RNN and Transformer
models performed competitively well on raw data compared to preprocessed data.
Especially in the 12-class land cover categorization setting, the difference
in performance was rather minor, speaking of $5\text{\,}\frac{\mathrm{\char
37\relax}}{}$ in accuracy and $0.03$ in kappa score. In terms of the $f_{1}$
score, the LSTM-RNN variant even achieved better performance on raw data
compared to preprocessed data, while is ranging behind the Transformer .
Throughout both the 23-class dataset (cf. Tables 1(c), 1(a) and 1(e)) and
12-class dataset (cf. Tables 1(d), 1(b) and 1(f)) variants, all evaluated deep
learning models performed similarly well on preprocessed data. Interestingly,
for the raw dataset partitions, the LSTM-RNN and Transformer models achieved
slightly better accuracy, kappa, and $f_{1}$ scores values compared to the MS-
ResNet and TempCNN variants. In the 23-class setting, the difference is rather
small with $12\text{\,}\frac{\mathrm{\char 37\relax}}{}$ in overall accuracy
and $0.01$ in kappa score which seems not significant considering their
reported variances. The 12-class case, however, confirms this observation with
a more pronounced difference of $0.030.07$ in kappa metric and
$5\text{\,}\frac{\mathrm{\char 37\relax}}{}11\text{\,}\frac{\mathrm{\char
37\relax}}{}$ in accuracy. In the next Section 8.3, we will investigate this
further by a feature importance analysis monitoring the backpropagation
process.

When comparing the $f_{1}$ scores in Tables 1(e) and 1(f) with the accuracy
scores in Tables 1(c) and 1(d) on raw data, we see that the convolutional
TempCNN model showed similar performance compared to the MS-ResNet model. The
$f_{1}$ scores, however, show a lower score for TempCNN . From these metrics,
we can derive that the TempCNN did classify the majority of field parcels
accurately but achieved lower accuracies on some of the more infrequent class
categories compared to MS-ResNet . This may be attributed to the shallower
network topology. The DuPLO achieves very good classification accuracy on
preprocessed satellite data, as it is on par with the other models. Since the
model architecture is precisely defined by the authors, we could not run the
model with different hyper parameter settings. Still it outperforms the other
models in some tasks on preprocessed data, see Tables 1(b), 1(a) and 1(e). On
the raw satellite data, it is behind Transformer and LSTM-RNN models.

Interestingly, the random forest baseline achieved competitive results to the
deep learning models on preprocessed time series data. It appears that
preprocessing, as a form of feature extraction, helped the random forest
classifier to improve its performance significantly. Without data
preprocessing on the raw dataset, the random forest baseline fell behind the
deep learning models. The kappa metric—that ranged $0.15$ worse in the
23-class land use variant and $0.1$ to $0.2$ worse in the 12-class land use
variant—indicates that infrequent classes were classified less accurately by
the random forest classifier. This got further confirmed by the poor $f_{1}$
score throughout both dataset variants where the random forest achieved worse
performance on raw datasets compared to preprocessed datasets. Overall, this
stresses the necessity of feature extraction for shallow machine learning
methods while demonstrating that—with sufficient manual effort in
preprocessing and, thus, feature extraction—random forests can achieve
competitive accuracies.

It may be concluded that the region- and domain-specific expert knowledge in
data preprocessing helped all models to achieve better accuracies compared to
the raw sentinel 2 time series dataset. All models achieved similar accuracies
for preprocessed datasets. Even the random forest classifier showed a
competitive performance to the deep learning models with similar scores in
overall accuracy, slightly worse kappa and $f_{1}$ scores. Both, the LSTM-RNN
and Transformer models, that rely on recurrence and self-attention, achieved
better accuracies on raw time series data compared to the convolutional models
MS-ResNet and TempCNN . This effect was minor in the land-use categorization
with 23-classes and more pronounced in the land-cover categorization, but
overall consistent throughout this evaluation. We will investigate the
difference of recurrence and self-attention compared to convolution further in
the next section.

###  8.2 Qualitative Model Evaluation

LabelsLSTM-RNNTransformerMSResNetTempCNNfallowgrasslandwinter wheatcornsummer
wheatwinter speltwinter ryewinter barleysummer barleysummer oatwinter
triticalerapeseedcorrectwrong Figure 5: Qualitative comparison of the
classification models from region in 4430000 E/W, 5528000 N/S in EPSG:31468

We conduct a qualitative analysis in addition to the quantitative model
comparison in model analysis in Section 8.1 and Table 1. Here, we use the deep
learning models to predict the 12 land cover classes on the raw time series
dataset. Each time series, as in Fig. 4(a), was obtained form the reflectance
values aggregated over one field geometry, as shown in the left-most figure.
The parcels were colored by crop type and were obtained from 4500m by 4500m
block (4430000 E/W, 5528000 N/S in EPSG:31468) from the test partition. The
figures towards the right show the fields that were correctly and wrongly
classified by the respective models. As indicated by Table 1, the
classification accuracies do not vary significantly throughout the models so
that no clear advantage between models is visible.

###  8.3 Temporal Feature Importance Assessment by Gradient Backpropagation

(a) Example of a _summer barley_ time series

(b) Example of a _corn_ time series

Figure 6: Input feature importance analysis by gradients of the input time
series $\boldsymbol{X}$

(a) Self-Attention Matrix of Head 1

cloudsground (b) Head 1 as bipartite graph visualization in context of the
input time series

(c) Self-Attention Matrix of Head 2

(d) Head 2 as bipartite graph visualization in context of the input time
series

Figure 7: Visual illustration of the self-attention scores from a pre-trained
transformer self-attention model

In the previous section, we observed that models based on recurrence and self-
attention achieved better accuracy metrics on noisy raw data compared to
convolutional models. Here, we investigate this further through a simple, yet
effective, feature importance analysis based on gradient backpropagation.

Deep learning models are differentiable functions $f_{\boldsymbol{\Theta}}$
that approximate a mapping from an input tensor $\boldsymbol{X}$ to a ground
truth tensor $\boldsymbol{y}$, as detailed in Section 3. In the training
process, gradients are backpropagated using Eq. 1 in order to adjust the model
weights ${\boldsymbol{\Theta}}$. Since the model $f_{\boldsymbol{\Theta}}$ is
fully differentiable, we can use the same mechanism to backpropagate gradients
further up to the input tensor $\boldsymbol{X}$. When we start from the
highest predicted score $y^{\ast}=\max(\hat{\boldsymbol{y}})$, we can
propagate gradients through the entire network back to the individual elements
$\boldsymbol{x}_{t}=(x_{t}^{\text{B01}},x_{t}^{\text{B02}},x_{t}^{\text{B03}},\dots)^{\\!\top\\!}$
of the input tensor
$\boldsymbol{X}=(\boldsymbol{x}_{0},\dots,\boldsymbol{x}_{T})$. These
gradients $\frac{\partial y^{\ast}}{\partial\boldsymbol{X}}$ reflect the
influence of each input data element on the current class prediction.
Vanishing gradients indicate that any change in the input has no effect on the
prediction. Positive or negative gradients would suggest an increase or
decrease of the predicted score if the input changed at these times. We would
like to emphasize that this analysis is model-agnostic, can be implemented in
few lines of code in current deep learning frameworks that utilize automatic
differentiation, and does not require ground truth labels.555Our
implementation is available at https://github.com/marccoru/crop-type-
mapping/blob/master/notebooks/FeatureImportance.ipynb.

For this experiment, we estimated the influence of each input time step on the
classification prediction for each of the evaluated networks, i.e., LSTM-RNN
(recurrence) and Transformer (self-attention), as well as MS-ResNet and
TempCNN (convolution). Figure 6 illustrates this by means of two separate
examples of two corn- and summer barley parcels. The top figures each show the
input time series $\boldsymbol{X}$ as a sequence of raw Sentinel 2
reflectances over the year of 2018. In the raw time series, we can identify
atmospheric noise and clouds as positive peaks in the data. These are caused
by the high reflectance values of clouds throughout all spectral bands. The
presence of a cloud at a given point in time does not provide any additional
information about the covered surface and should be, thus, considered
irrelevant for the classification. The following rows display the gradients
$\frac{\partial y^{\ast}}{\partial\boldsymbol{X}}$ for each of the respective
models. These plots indicate the influence of the measurement at the
particular input time in the top row on the classification prediction of the
respective model.

The results of this experiment shows that the gradients through the LSTM-RNN
and Transformer networks were only non-zero for comparatively few
observations. This indicates that only a few key time-points are were
necessary to extract the classification relevant information from the raw time
series. Only time points where the surface was not obscured by clouds had non-
zero gradients. Time instances with positive peaks in reflectance, likely
induced by clouds, did not have influence the classification prediction, as
indicated by vanishing gradients. This experimentally shows that the LSTM-RNN
and Transformer models have been optimized to automatically identify and
suppress time instances that are not classification relevant, e.g., though
cloud coverage, learned purely from data.

The convolution-based model architectures MS-ResNet and TempCNN show, in
general, a similar behavior. However, these models seem to have some non-zero
gradients on cloudy observations. This is especially true for the
comparatively shallow TempCNN architecture where the cloudy observations at
$t=30$ in Fig. 6(a) and $t=38$ in Fig. 6(b) still have influence in the
classification prediction. The MS-ResNet model also seems to to be able to
extract features from the entire temporal range. Since the data ranges from
January to December, it is expected that the particular crop class is only
visible during the summer periods. The MS-ResNet model seems to still include
features from the winter periods.

Summarizing the above, the recurrence and self-attention based models LSTM-RNN
and Transformer seem to extract selective temporal features, while the
convolution-based models utilize the entire time series. Also, the LSTM-RNN
and Transformer models were able to transition from a time observation with
high gradients to zero gradients without any time delay. It appears that the
convolution-based models TempCNN and MS-ResNet require at least two subsequent
observations to reduce the influence of the particular time instance. This can
be observed from Fig. 6(a) at $t=30$ and $t=34$: The MS-ResNet model could
reduce the influence of the observations at $t=30$ (indicated by zero
gradients) when three subsequent observations were cloudy. However, when only
one observation was cloudy, as in $t=34$, MS-ResNet was unable to suppress
this observation completely. In contrast, the models LSTM-RNN and Transformer
were able to suppress these cloudy observations even if they appeared in a
single time instance.

We attribute this difference in the influence of input observations to the
respective mechanisms for feature extraction. Recurrent networks utilize
internal gates that can control the influence of the particular time instance
to a hidden memory state. Our previous work [66], which visualized internal
LSTM states, supports this hypothesis for recurrence. Similarly, self-
attention enables a model to select specific observations by assigning a large
attention score to them. In contrast, convolutions always extract features
from a local neighborhood. Hence, it seems to be more difficult for
convolutional architectures to ignore sudden appearences of irrelevant
observations within a temporal sequence, as indicated by the non-zero
gradients at cloudy observations.

###  8.4 Qualitative Analysis of Self-Attention Scores

In the previous section, we experimentally observed that the Transformer model
was able to suppress the influence of an observation for classification-
irrelevant observations, e.g., clouds. Here, we concentrate on the self-
attention mechanism (cf. Section 3.4) that was implemented in the Transformer
model and analyse the attention scores on the example of the cornfield parcels
from the previous experiment in shown in Fig. 6(b).

The Transformer models realize multi-headed self-attention as multiple
attention mechanisms in parallel. To recall, every self-attention mechanism
calculates an attention matrix $\boldsymbol{A}\in[0,1]^{T\times T}$ using the
softmax operation following Eq. 3. These attention scores define the influence
of an input time feature on a higher-level output time feature. We visualize
the values of this matrix in Figs. 7(a) and 7(c) for three attention heads
from the first self-attention layer (of three). This matrix can be seen as an
adjacency matrix between input nodes and output nodes. In Figs. 7(b) and 7(d),
we alternatively show the same matrix as bipartite graphs. The former matrix
elements are now shown as weighted directed edges between input and output
nodes. The strength of the attention is indicated by the opacity of the edge.
Above the attention graph, we show the input time series of the example as a
reference to allow for visual interpretation of the time instances. Here, the
attention scores are drawn as vertical lines corresponding to the values in
$\boldsymbol{A}$.

From these figures, we can directly observe that the self-attention scores
focus on distinct events with each head. While the first head appears to focus
on the early-middle part of the time series, the second head extracts features
from the early and late time series. Also, the attention scores seem to re-
distribute features over time. Hence, even though hidden feature tensors
throughout the network still maintain a temporal dimension, the temporal
consistency with the input time series loses context. Consistent with the
previous experiment in Section 8.3, we observed that the attention scores did
not focus on the strong positive peaks in the time series which indicate a
cloudy observation. Hence, we can conclude that self-attention mechanisms are
a key tool in suppressing the non-classification-relevant cloudy observations
which explain the zero-gradients of the previous experiment.

###  8.5 Comparison of Self-Attention with Soft Attention as implemented in
DuPLO

(a)

(b)

Figure 8: Comparison of self-attention scores, as implemented in the
transformer model, with soft-attention scores used in the DuPLO model

In the previous experiment, we explored the self-attention mechanism, as
implemented in cascaded layers of the transformer model. The DuPLO model[59],
as described in Section 4.5, proposes a neural network architecture that
combines convolutional layers, recurrent layers and a soft-attention
mechnanism and has been employed for crop type mapping in remote sensing. In
Table 1, observed that the duplo architecture achieved good accuracies on
preprocessed data, but was outperformed by a purely recurrent model and the
self-attention model on the raw time series data. This yields the question of
the effectiveness of soft-attention (Section 3.5) compared to self-attention
(Section 3.4) which we address in this experiment.

In Fig. 8, we draw the attention scores from two test samples shown in the top
row. In the second row, we see the average self-attention scores
$\boldsymbol{\alpha}=\frac{1}{T_{out}}\sum_{i=0}^{T_{out}}\boldsymbol{A}_{i}$
from the three attention heads of the first transformer layer per input time.
This corresponds to the average over the columns in Figs. 7(c) and 7(a). The
third row shows the soft-attention scores, as obtained from the duplo model.

Comparing the attention scores, we observe that, consistent with the previous
experiment, the self-attention scores focus on few single time instance that
obtain classification-relevant features. In contrast, the comparatively
simpler soft-attention is not able to focus on single time instances so that
multiple neighboring time instances have attention scores greater than zero.
Note that the soft-attention scores are calculated within the recurrent branch
of the Duplo model that uses a complex combination of convolutional layers
combined with a single mono-directional recurrent layer. The limitations of
these attention scores may be also caused by the preceding layers. For
instance, the mono-directional encoding of the gated recurrent unit of the
duplo model likely leads to higher-level features in the latter half of the
time series. It is likely that this effect causes the higher attention scores
focusing on the later time instances where more expressive recurrent features
have been extracted from the single mono-directional GRU layer.

Overall, we conclude that the self-attention mechanism can be, by its design
with three dense layers, more expressive than the simpler soft-attention that
requires only one dense layer. However, the architecture of the DuPLO model,
that comprises a mixture of convolutional, recurrent, and attention layers,
makes it difficult to isolate the effects of the individual layers in
comparison to the other model architectures that exclusively rely on one type
of layer architecture.

###  8.6 Feature Separability Analyes through Manifold Embeddings

Deep learning models, in general, extract features of increasingly complexity
throughout the cascaded layer architectures [91]. In this experiment, we
analyzed this property by visualizing the hidden features at varying deeper
layers of the transformer architecture, as shown in Fig. 9. For this anaysis,
we extract $D=128$-dimensional feature vectors after encoding time series
samples from the test-dataset. We project them into a two-dimensional
embedding space via _t-distributed stochastic neighbor embedding (t-SNE)_
[92], _Uniform manifold approximation and projection (UMAP)_ [93], as used in
[94], and _Principal Component Analysis (PCA)_. We used a trained transformer
network with three self-attention layers and $D=128$ hidden dimensions. To
obtain the features, we passed the entire test dataset over all regions
through a pre-trained transformer network and recorded the hidden features
after each self-attention layer and after the last dense layer before max
pooling and softmax activation, as visualized in Fig. 9. The original features
$\boldsymbol{H}\in\mathbb{R}^{(T=70)\times(D=128)}$ were then averaged along
the time dimension to obtain $\boldsymbol{h}^{\prime}\in\mathbb{R}^{D=128}$
and mapped into the final embedding space $\mathbb{R}^{D=2}$. The illustration
shows the representations of the embedded 200 samples for each of the 23 land
use classes. We augmented the final plots with a key indicating the classes
identities and a schematic picture of the network topology. Since the class
definitions are rather broad with, for instance, multiple categories of
grassland, we grouped them by color but assigned varying marker shapes to each
class. Note, that the embeddings rely on the extracted features without taking
the actual ground truth class label into account. One can easily observe that
the transformed features become more and more separable with each deeper
network layer throughout all embedding techniques. This becomes apparent by
the formation clusters throughout the layers. Looking at the last layer, which
is closest to the final prediction, we can identify distinct separable
manifolds for corn () and rapeseed (). This suggests that these classes are
well-separable in the feature space. This observation is clearly observable
for both the non-linear T-SNE and UMAP projections while the separability in
the first two principle components in the linear PCA projection is less
pronounced.

From expert knowledge provided by the _Bavarian Ministry of Agriculture_ , we
account for a certain degree of wrong labels in the dataset. The
identification of these outliers is a key motivation for crop type monitoring.
Hence, these few field parcels may have been assigned wrong labels by the
respective farmer.

In summary, the feature extraction capabilities learned by this self-attention
network is consistent with assumptions based on domain-knowledge of the
particular classes. Methodologically, we observed that features from deeper
cascaded network layers get increasingly distinctive and field parcels of the
same label are mapped in similar embedded regions. We would like to stress
that this type of analysis is not unique to self-attention networks and can be
reproduced for many deep learning architectures that extracts features in
cascaded layers.

(a)(b)(c)UMAPT-SNEPCAUMAPT-SNEPCAUMAPT-SNEPCA

$\boldsymbol{X}$densepooling$\hat{\boldsymbol{y}}$(a)(b)(c)

Figure 9: Feature embedding at varying depths of the Transformer Network. Each
self-attention layer increases separability of the extracted features.

##  9 Conclusion

In this work, we quantitatively and qualitatively analyzed self-attention for
the application in multi-temporal Earth observation. We performed a large-
scale quantitative comparison in Section 8.1 where we evaluated multiple model
architectures that rely on self-attention, recurrence and convolution, and a
random forest baseline. We compared these models from multiple angles by
reporting their performance on preprocessed and raw Sentinel 2 time series as
well as land-use and a land-cover oriented set of classes. Here, we observed
that all models performed equally well on preprocessed data. Even the random
forest baseline achieved competitive overall accuracy. This leads to the
conclusion that the choice of model architecture is less critical when
extensive data preprocessing is utilized as a form of feature extraction aided
by region-specific expert knowledge. For raw unprocessed Sentinel 2 time
series data, the Transformer , and LSTM-RNN architectures were able to achieve
better accuracies compared to the convolutional models. We investigated this
further by a feature importance analysis in Section 8.3 using gradients where
we observed that the mechanisms self-attention and recurrence helped to
suppress information not relevant for classification in the time series. For
the Transformer , this was realized by learning weights in the attention
mechanism which enables the model to specifically focus on some observations,
as could be observed in Fig. 1 in Section 8.4. Finally, we looked at the
larger Transformer network topology and observed how deeper neural network
layers were able to learn increasingly separable representations of the
classes in Section 8.6. Since we chose a challenging set of classes, class
overlaps in the semantic representation of classes with common properties were
present and followed common intuition. This is a typical challenge for land
cover and land use classification with long-tailed class distributions.

To summarize, we saw that self-attention is a promising technique that allows
neural networks to extract features from specific time instances with raw
optical satellite time series. It suits-well in the canon of time series
classification models that utilize recurrent or convolutional layers. We did
not find any mechanism among the compared baseline methods that systematically
achieved better accuracies, but observed that self-attention and recurrence
were, by design of the feature extraction, more robust to noise in the data
and could better suppress cloudy observations in raw time series.

In the future, a robust classification of raw time series data without region-
specific expert knowledge will be key to quantitatively exploit the satellite
data that is published daily. We hope to have contributed a step towards this
direction by this comprehensive study that evaluated model accuracies on
preprocessed and raw data for a variety of mechanisms for time series
classification.

##  10 Acknowledgements

We would like to thank Sabine Erbe from the Bavarian Ministry of Agriculture
for the continuous and persistent support and valuable domain knowledge as
well as Markus Ziegler and Christoph Schmidt from GAF-AG for their
collaboration in the joint crop type classification project which provided
objective, method and preprocessed time series data for this study. The work
of Marc Rußwurm was funded by the _German Federal Ministry for Economic
Affairs and Energy (BMWi)_ under reference _50EE1908_.

## References

  * [1] C. E. Woodcock, R. Allen, M. Anderson, A. Belward, R. Bindschadler, W. Cohen, F. Gao, S. N. Goward, D. Helder, E. Helmer, et al., Free access to landsat imagery, Science 320 (5879) (2008) 1011–1011.
  * [2] M. A. Wulder, T. R. Loveland, D. P. Roy, C. J. Crawford, J. G. Masek, C. E. Woodcock, R. G. Allen, M. C. Anderson, A. S. Belward, W. B. Cohen, et al., Current status of landsat program, science, and applications, Remote sensing of environment 225 (2019) 127–147.
  * [3] B. C. Reed, J. F. Brown, D. VanderZee, T. R. Loveland, J. W. Merchant, D. O. Ohlen, Measuring Phenological Variability from Satellite Imagery, Journal of Vegetation Science 5 (5) (1994) 703–714.  doi:10.2307/3235884.
  * [4] J. B. Odenweller, K. I. Johnson, Crop identification using Landsat temporal-spectral profiles, Remote Sensing of Environment 14 (1-3) (1984) 39–5.  doi:10.1016/0034-4257(84)90006-3.
  * [5] M. A. Wulder, W. A. Kurz, M. Gillis, et al., National level forest monitoring and modeling in canada, Progress in Planning 61 (4) (2004) 365–381.
  * [6] M. Herold, S. Carter, V. Avitabile, A. B. Espejo, I. Jonckheere, R. Lucas, R. E. McRoberts, E. Næsset, J. Nightingale, R. Petersen, et al., The role and need for space-based forest biomass-related measurements in environmental management and policy, Surveys in Geophysics 40 (4) (2019) 757–778.
  * [7] J. Reiche, R. Lucas, A. L. Mitchell, J. Verbesselt, D. H. Hoekman, J. Haarpaintner, J. M. Kellndorfer, A. Rosenqvist, E. A. Lehmann, C. E. Woodcock, et al., Combining satellite data for better tropical forest monitoring, Nature Climate Change 6 (2016) 120–122.
  * [8] Sentinel data access annual report 2019, date 06/05/19, COPE-SERCO-RP-19-0389 (2019).
  * [9] S. Foerster, K. Kaden, M. Foerster, S. Itzerott, Crop type mapping using spectral-temporal profiles and phenological information, Computers and Electronics in Agriculture 89 (2012) 30–40.  doi:10.1016/j.compag.2012.07.015.
  * [10] C. Conrad, S. Fritsch, J. Zeidler, G. Rücker, S. Dech, Per-Field Irrigated Crop Classification in Arid Central Asia Using SPOT and ASTER Data, Remote Sensing 2 (4) (2010) 1035–1056.  doi:10.3390/rs2041035.
  * [11] C. Conrad, S. Dech, O. Dubovyk, S. Fritsch, D. Klein, F. Löw, G. Schorcht, J. Zeidler, Derivation of temporal windows for accurate crop discrimination in heterogeneous croplands of Uzbekistan using multitemporal RapidEye images, Computers and Electronics in Agriculture 103 (2014) 63–74.  doi:10.1016/j.compag.2014.02.003.
  * [12] P. Hao, Y. Zhan, L. Wang, Z. Niu, M. Shakir, Feature Selection of Time Series MODIS Data for Early Crop Classification Using Random Forest: A Case Study in Kansas, USA, Remote Sensing 7 (5) (2015) 5347–5369.  doi:10.3390/rs70505347.
  * [13] J. M. Peña-Barragán, M. K. Ngugi, R. E. Plant, J. Six, Object-based crop identification using multiple vegetation indices, textural features and crop phenology, Remote Sensing of Environment 115 (6) (2011) 1301–1316.  doi:10.1016/j.rse.2011.01.009.
  * [14] L. Zhong, L. Hu, H. Zhou, Deep learning based multi-temporal crop classification, Remote sensing of environment 221 (2019) 430–443.
  * [15] S. Valero, D. Morin, J. Inglada, G. Sepulcre, M. Arias, O. Hagolle, G. Dedieu, S. Bontemps, P. Defourny, B. Koetz, Production of a dynamic cropland mask by processing remote sensing image series at high temporal and spatial resolutions, Remote Sensing 8 (1) (2016) 1–21.  doi:10.3390/rs8010055.
  * [16] P. Kumar, D. K. Gupta, V. N. Mishra, R. Prasad, Comparison of support vector machine, artificial neural network, and spectral angle mapper algorithms for crop classification using LISS IV data, International Journal of Remote Sensing 36 (6) (2015) 1604–1617.  doi:10.1080/2150704X.2015.1019015.
  * [17] B. Zheng, S. W. Myint, P. S. Thenkabail, R. M. Aggarwal, A support vector machine to identify irrigated crop types using time-series landsat ndvi data, International Journal of Applied Earth Observation and Geoinformation 34 (2015) 103–112.
  * [18] R. Devadas, R. Denham, M. Pringle, Support vector machine classification of object-based data for crop mapping, using multi-temporal landsat imagery, Int. Arch. Photogramm. Remote Sens. Spat. Inf. Sci 39 (2012) 185–190.
  * [19] C. Ünsalan, K. L. Boyer, Review on Land Use Classification, in: Multispectral Satellite Image Understanding: From Land Classification to Building and Road Detection, Springer, 2011, pp. 49–64.  doi:10.1007/978-0-85729-667-2_5.
  * [20] Y. Shao, R. S. Lunetta, Comparison of support vector machine, neural network, and cart algorithms for the land-cover classification using limited training data points, ISPRS Journal of Photogrammetry and Remote Sensing 70 (2012) 78–87.
  * [21] P. Jönsson, L. Eklundh, Timesat—a program for analyzing time-series of satellite sensor data, Computers & Geosciences 30 (8) (2004) 833–845.
  * [22] L. Eklundh, P. Jönsson, Timesat for processing time-series data from satellite sensors for land surface monitoring, in: Multitemporal Remote Sensing, Springer, 2016, pp. 177–194.
  * [23] M. A. White, K. M. de Beurs, K. Didan, D. W. Inouye, A. D. Richardson, O. P. Jensen, J. O’KEEFE, G. Zhang, R. R. Nemani, W. J. van Leeuwen, et al., Intercomparison, interpretation, and assessment of spring phenology in north america estimated from remote sensing for 1982–2006, Global Change Biology 15 (10) (2009) 2335–2359.
  * [24] L. Olsson, L. Eklundh, J. Ardö, A recent greening of the sahel—trends, patterns and potential causes, Journal of Arid Environments 63 (3) (2005) 556–566.
  * [25] K. Jia, S. Liang, X. Wei, Y. Yao, Y. Su, B. Jiang, X. Wang, Land cover classification of landsat data with phenological features extracted from time series modis ndvi data, Remote sensing 6 (11) (2014) 11518–11532.
  * [26] M. Singha, B. Wu, M. Zhang, An object-based paddy rice classification using multi-spectral data and crop phenology in assam, northeast india, Remote Sensing 8 (6) (2016) 479.
  * [27] Z. Zhu, C. E. Woodcock, Continuous change detection and classification of land cover using all available landsat data, Remote sensing of Environment 144 (2014) 152–171.
  * [28] Z. Zhu, C. E. Woodcock, Object-based cloud and cloud shadow detection in landsat imagery, Remote sensing of environment 118 (2012) 83–94.
  * [29] Z. Zhu, S. Wang, C. E. Woodcock, Improvement and expansion of the fmask algorithm: Cloud, cloud shadow, and snow detection for landsats 4–7, 8, and sentinel 2 images, Remote Sensing of Environment 159 (2015) 269–277.
  * [30] J. O. Street, R. J. Carroll, D. Ruppert, A note on computing robust regression estimates via iteratively reweighted least squares, The American Statistician 42 (2) (1988) 152–154.
  * [31] W. Dumouchel, F. O’Brien, et al., Integrating a robust option into a multiple regression computing environment, in: Computer science and statistics: Proceedings of the 21st symposium on the interface, American Statistical Association Alexandria, VA, 1989, pp. 297–302.
  * [32] J. Verbesselt, R. Hyndman, G. Newnham, D. Culvenor, Detecting trend and seasonal changes in satellite image time series, Remote sensing of Environment 114 (1) (2010) 106–115.
  * [33] R. E. Kennedy, Z. Yang, W. B. Cohen, Detecting trends in forest disturbance and recovery using yearly landsat time series: 1. landtrendr—temporal segmentation algorithms, Remote Sensing of Environment 114 (12) (2010) 2897–2910.
  * [34] W. S. McCulloch, W. Pitts, A logical calculus of the ideas immanent in nervous activity, The bulletin of mathematical biophysics 5 (4) (1943) 115–133.
  * [35] J. D. Cowan, Neural networks: the early days, in: Advances in neural information processing systems, 1990, pp. 828–842.
  * [36] Y. LeCun, L. Bottou, Y. Bengio, P. Haffner, et al., Gradient-based learning applied to document recognition, Proceedings of the IEEE 86 (11) (1998) 2278–2324.
  * [37] J. Deng, W. Dong, R. Socher, L.-J. Li, K. Li, L. Fei-Fei, Imagenet: A large-scale hierarchical image database, in: 2009 IEEE conference on computer vision and pattern recognition, IEEE, 2009, pp. 248–255.
  * [38] 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 115 (3) (2015) 211–252.
  * [39] K. Simonyan, A. Zisserman, Very deep convolutional networks for large-scale image recognition, arXiv preprint arXiv:1409.1556V6.
  * [40] A. Krizhevsky, I. Sutskever, G. E. Hinton, Imagenet classification with deep convolutional neural networks, in: Advances in neural information processing systems, 2012, pp. 1097–1105.
  * [41] K. He, X. Zhang, S. Ren, J. Sun, Deep residual learning for image recognition, in: Proceedings of the IEEE conference on computer vision and pattern recognition, 2016, pp. 770–778.
  * [42] J. Long, E. Shelhamer, T. Darrell, Fully convolutional networks for semantic segmentation, in: Proceedings of the IEEE conference on computer vision and pattern recognition, 2015, pp. 3431–3440.
  * [43] C. Szegedy, S. Ioffe, V. Vanhoucke, A. A. Alemi, Inception-v4, inception-resnet and the impact of residual connections on learning, in: Thirty-First AAAI Conference on Artificial Intelligence, 2017.
  * [44] A. Vaswani, N. Shazeer, N. Parmar, J. Uszkoreit, L. Jones, A. N. Gomez, Ł. Kaiser, I. Polosukhin, Attention is all you need, in: Advances in neural information processing systems, 2017, pp. 5998–6008.
  * [45] J. Devlin, M.-W. Chang, K. Lee, K. Toutanova, Bert: Pre-training of deep bidirectional transformers for language understanding, arXiv preprint arXiv:1810.04805V2.
  * [46] A. Radford, J. Wu, R. Child, D. Luan, D. Amodei, I. Sutskever, Language models are unsupervised multitask learners, OpenAI Blog 1 (8).
  * [47] D. Marmanis, M. Datcu, T. Esch, U. Stilla, Deep learning earth observation classification using imagenet pretrained networks, IEEE Geoscience and Remote Sensing Letters 13 (1) (2015) 105–109.
  * [48] M. Volpi, D. Tuia, Dense semantic labeling of subdecimeter resolution images with convolutional neural networks, IEEE Transactions on Geoscience and Remote Sensing 55 (2) (2016) 881–893.
  * [49] J. Sherrah, Fully convolutional networks for dense semantic labelling of high-resolution aerial imagery, arXiv preprint arXiv:1606.02585V1.
  * [50] N. Audebert, B. Le Saux, S. Lefèvre, Semantic segmentation of earth observation data using multimodal and multi-scale deep networks, in: Asian conference on computer vision, Springer, 2016, pp. 180–196, v1.
  * [51] G. Cheng, P. Zhou, J. Han, Learning rotation-invariant convolutional neural networks for object detection in vhr optical remote sensing images, IEEE Transactions on Geoscience and Remote Sensing 54 (12) (2016) 7405–7415.
  * [52] G. Cheng, J. Han, X. Lu, Remote sensing image scene classification: Benchmark and state of the art, Proceedings of the IEEE 105 (10) (2017) 1865–1883.
  * [53] G. Cheng, C. Yang, X. Yao, L. Guo, J. Han, When deep learning meets metric learning: Remote sensing image scene classification via learning discriminative cnns, IEEE Transactions on Geoscience and Remote Sensing 56 (5) (2018) 2811–2821.
  * [54] M. Rußwurm, M. Korner, Temporal vegetation modelling using long short-term memory networks for crop identification from medium-resolution multi-spectral satellite images, in: Proceedings of the IEEE Conference on Computer Vision and Pattern Recognition Workshops, 2017, pp. 11–19.
  * [55] X. Jia, A. Khandelwal, G. Nayak, J. Gerber, K. Carlson, P. West, V. Kumar, Incremental Dual-memory LSTM in Land Cover Prediction, in: 23rd ACM SIGKDD International Conference on Knowledge Discovery and Data Mining, 2017, pp. 867–876.  doi:10.1145/3097983.3098112.
  * [56] H. Lyu, H. Lu, L. Mou, Learning a Transferable Change Rule from a Recurrent Neural Network for Land Cover Change Detection, Remote Sensing 8 (12) (2016) 506\.  doi:10.3390/rs8060506.
  * [57] L. Mou, L. Bruzzone, X. X. Zhu, Learning spectral-spatial-temporal features via a recurrent convolutional neural network for change detection in multispectral imagery, IEEE Transactions on Geoscience and Remote Sensing 57 (2) (2018) 924–935.
  * [58] P. Benedetti, D. Ienco, R. Gaetano, K. Ose, R. G. Pensa, S. Dupuy, M3fusion: A deep learning architecture for multiscale multimodal multitemporal satellite data fusion, IEEE Journal of Selected Topics in Applied Earth Observations and Remote Sensing 11 (12) (2018) 4939–4949.
  * [59] R. Interdonato, D. Ienco, R. Gaetano, K. Ose, DuPLO: A DUal view Point deep Learning architecture for time series classificatiOn, ISPRS Journal of Photogrammetry and Remote Sensing 149 (2019) 91–104.  doi:10.1016/j.isprsjprs.2019.01.011.
  * [60] P. J. Werbos, Backpropagation through time: what it does and how to do it, Proceedings of the IEEE 78 (10) (1990) 1550–1560.
  * [61] S. Hochreiter, J. Schmidhuber, Long Short-Term Memory, Neural Computation 9 (8) (1997) 1735–1780.  doi:10.1162/neco.1997.9.8.1735.
  * [62] J. Chung, C. Gulcehre, K. Cho, Y. Bengio, Empirical evaluation of gated recurrent neural networks on sequence modeling, arXiv preprint arXiv:1412.3555V1.
  * [63] I. Sutskever, O. Vinyals, Q. V. Le, Sequence to sequence learning with neural networks, in: Advances in neural information processing systems, 2014, pp. 3104–3112.
  * [64] V. S. F. Garnot, L. Landrieu, S. Giordano, N. Chehata, Time-Space Tradeoff in Deep Learning Models for Crop Classification on Satellite Multi-Spectral Image Time Series, in: IEEE International Geoscience and Remote Sensing Symposium (IGARSS), 2019, pp. 6247–6250.  doi:10.1109/IGARSS.2019.8900517.
  * [65] A. Sharma, X. Liu, X. Yang, Land cover classification from multi-temporal, multi-spectral remotely sensed imagery using patch-based recurrent neural networks, Neural Networks 105 (2018) 346–355.
  * [66] M. Rußwurm, M. Körner, Multi-temporal land cover classification with sequential recurrent encoders, ISPRS International Journal of Geo-Information 7 (4) (2018) 129.
  * [67] D. Bahdanau, J. Chorowski, D. Serdyuk, P. Brakel, Y. Bengio, End-to-end attention-based large vocabulary speech recognition, in: 2016 IEEE international conference on acoustics, speech and signal processing (ICASSP), IEEE, 2016, pp. 4945–4949.
  * [68] H. I. Fawaz, G. Forestier, J. Weber, L. Idoumghar, P.-A. Muller, Deep learning for time series classification: a review, Data Mining and Knowledge Discovery 33 (4) (2019) 917–963.
  * [69] C. Pelletier, G. I. Webb, F. Petitjean, Temporal convolutional neural network for the classification of satellite image time series, Remote Sensing 11 (5) (2019) 523.
  * [70] V. S. F. Garnot, L. Landrieu, S. Giordano, N. Chehata, Satellite Image Time Series Classification with Pixel-Set Encoders and Temporal Self-Attention, arXiv e-printsV1.  arXiv:1911.07757.
  * [71] D. E. Rumelhart, G. E. Hinton, R. J. Williams, Learning representations by back-propagating errors, Nature 323 (6088) (1986) 533–536.  doi:10.1038/323533a0.
  * [72] R. Jozefowicz, W. Zaremba, I. Sutskever, An Empirical Exploration of Recurrent Network Architectures, in: Proceedings of the 32nd International Conference on Machine Learning (ICML), Proceedings of Machine Learning Research, 2015, pp. 2342–2350.
  * [73] B. Yoshua, S. Patrice, F. Paolo, Learning long-term dependencies with gradient descent is difficult, IEEE transactions on neural networks 5 (2) (1994) 157–166.
  * [74] H. Ismail Fawaz, B. Lucas, G. Forestier, C. Pelletier, D. F. Schmidt, J. Weber, G. I. Webb, L. Idoumghar, P.-A. Muller, F. Petitjean, Inceptiontime: Finding alexnet for time series classification, ArXivV1.
  * [75] N. Hatami, Y. Gavet, J. Debayle, Classification of time-series images using deep convolutional neural networks, in: Tenth International Conference on Machine Vision (ICMV 2017), Vol. 10696, International Society for Optics and Photonics, 2018, p. 106960Y.
  * [76] D. Britz, M. Y. Guan, M.-T. Luong, Efficient attention using a fixed-size memory representation, arXiv preprint arXiv:1707.00110V1.
  * [77] S. Ioffe, C. Szegedy, Batch normalization: Accelerating deep network training by reducing internal covariate shiftV3.  arXiv:1502.03167v3.
  * [78] Z. Wang, W. Yan, T. Oates, Time series classification from scratch with deep neural networks: A strong baseline, in: 2017 international joint conference on neural networks (IJCNN), IEEE, 2017, pp. 1578–1585.
  * [79] Z. Cui, W. Chen, Y. Chen, Multi-scale convolutional neural networks for time series classification, arXiv preprint arXiv:1603.06995V4.
  * [80] H. A. Dau, A. Bagnall, K. Kamgar, C.-C. M. Yeh, Y. Zhu, S. Gharghabi, C. A. Ratanamahatana, E. Keogh, The ucr time series archive, arXiv preprint arXiv:1810.07758V2.
  * [81] A. Bagnall, H. A. Dau, J. Lines, M. Flynn, J. Large, A. Bostrom, P. Southam, E. Keogh, The uea multivariate time series classification archive, 2018, arXiv preprint arXiv:1811.00075V1.
  * [82] F. Pedregosa, G. Varoquaux, A. Gramfort, V. Michel, B. Thirion, O. Grisel, M. Blondel, P. Prettenhofer, R. Weiss, V. Dubourg, et al., Scikit-learn: Machine learning in python, Journal of machine learning research 12 (Oct) (2011) 2825–2830.
  * [83] D. P. Kingma, J. Ba, Adam: A method for stochastic optimization, arXiv preprint arXiv:1412.6980V9.
  * [84] J. Cohen, A coefficient of agreeement for nominal scales, Educational and Psychological Measurement 20 (1960) 37–46.
  * [85] J. Bergstra, D. Yamins, D. D. Cox, Hyperopt: A python library for optimizing the hyperparameters of machine learning algorithms, in: Proceedings of the 12th Python in science conference, Citeseer, 2013, pp. 13–20.
  * [86] L. Li, K. Jamieson, A. Rostamizadeh, E. Gonina, M. Hardt, B. Recht, A. Talwalkar, Massively parallel hyperparameter tuning, arXiv preprint arXiv:1810.05934V5.
  * [87] R. Liaw, E. Liang, R. Nishihara, P. Moritz, J. E. Gonzalez, I. Stoica, Tune: A research platform for distributed model selection and training, arXiv preprint arXiv:1807.05118V1.
  * [88] P. Schratz, J. Muenchow, E. Iturritxa, J. Richter, A. Brenning, Performance evaluation and hyperparameter tuning of statistical and machine-learning models using spatial data, arXiv preprint arXiv:1803.11266V1.
  * [89] A. Brenning, Spatial cross-validation and bootstrap for the assessment of prediction rules in remote sensing: The r package sperrorest, in: 2012 IEEE international geoscience and remote sensing symposium, IEEE, 2012, pp. 5372–5375.
  * [90] N. Gorelick, M. Hancher, M. Dixon, S. Ilyushchenko, D. Thau, R. Moore, Google earth engine: Planetary-scale geospatial analysis for everyone, Remote Sensing of Environment 202 (2017) 18–27.
  * [91] I. Goodfellow, Y. Bengio, A. Courville, Deep learning, MIT press, 2016.
  * [92] L. v. d. Maaten, G. Hinton, Visualizing data using t-sne, Journal of machine learning research 9 (Nov) (2008) 2579–2605.
  * [93] L. McInnes, J. Healy, J. Melville, Umap: Uniform manifold approximation and projection for dimension reduction, arXiv preprint arXiv:1802.03426V2.
  * [94] F. Mohammadimanesh, B. Salehi, M. Mahdianpari, E. Gill, M. Molinier, A new fully convolutional neural network for semantic segmentation of polarimetric sar imagery in complex land cover ecosystem, ISPRS journal of photogrammetry and remote sensing 151 (2019) 223–236.