Self-attention for raw optical Satellite Time Series Classification

In plain words

Satellites such as Sentinel-2 photograph the same field every few days, so each field becomes a time series of light readings. Before a classifier sees these readings, experts usually clean them up: they remove cloudy days, correct the atmosphere and fill gaps. This paper asks how much that cleaning matters, and which neural network copes best with the messy raw readings. It compares a Transformer (the self-attention model) with a recurrent LSTM, two convolutional networks and a random forest, on crop types in Bavaria.

Cleaning helps every model, and once the data is clean the choice of model matters much less. On raw data the Transformer and the LSTM do better than the convolutional networks. A gradient-based check shows why: they hardly react to cloudy days, and the Transformer's attention lands on a few useful observations.

For a practitioner, the lesson is that a model can learn to skip noise and so reduce hand-built preprocessing. The paper's own numbers still show that clean data scores higher, so the trade-off is real.

TL;DR

• Crop-type classification on 13-band Sentinel-2 time series, on raw and on industry-preprocessed data, for three Bavarian regions §7 §7.2.
• Preprocessing raised every model's scores. The choice of model mattered less on preprocessed data §8.1 tab. 1.
• On raw data, a Transformer encoder and a bidirectional LSTM beat the two CNNs (MS-ResNet, TempCNN) on kappa and accuracy, author-reported tab. 1.
• Gradient analysis and attention maps suggest the Transformer ignores cloud peaks by attending to a few observations §8.3 §8.4.

Problem

Only a small share of the Sentinel-2 imagery published each year is downloaded and analysed §1. Standard pipelines rely on atmospheric correction, cloud masking and gap filling, which need region-specific expertise and considerable computing effort §1 §7.2.2. Sentinel-2 observes the same spot every two to five days §7.2, so cloud cover and gaps are a constant issue. The question is whether a model can learn to handle raw top-of-atmosphere observations without that step §1. The comparison covers recurrence, 1-D convolution and self-attention, which the authors name as three mechanisms for this task §1 §3.

Method

The paper trains four neural networks and one random forest on the same parcels and compares them with accuracy, Cohen's kappa and class-mean F1 §8.1. Each parcel is one time series: one vector of 13 spectral bands per date, averaged over the parcel's pixels §7.2.1. The network outputs a softmax over the crop classes §4.2 §3.

The Transformer is the focus §4.2. It reuses the encoder of the original Transformer (arXiv:1706.03762) but drops word embedding, because reflectance values are already continuous §4.2. A positional encoding is added, because self-attention on its own ignores order §4.2. Each block contains multi-head self-attention, then dense layers applied to each date, with skip connections and layer normalisation §4.2. After the last block, features are max-pooled over time into one vector, then a dense layer and softmax give class scores §4.2.

Self-attention in the paper is the standard form: queries, keys and values are linear maps of the same input, and the output mixes the values with weights that come from a softmax over scaled query-key products §3.4 eq. 4. The paper writes the scale as the square root of the hidden dimension D_h eq. 4; the code uses the square root of the per-head key size (see Paper vs code).

\mathbf{H} = \mathrm{softmax}\!\left(\frac{W_Q^{\top} X \, (W_K^{\top} X)^{\top}}{\sqrt{D_h}}\right) W_V^{\top} X \qquad [\text{eq. 4}]

Architecture as implemented (derived from code; the paper leaves these details out). The order of operations below follows the repository:

# Transformer classifier for one batch of raw series  [§4.2]
# Input:  series (B, T, C) with C = 13 bands, T = 70 dates
# Output: log-probabilities (B, K) over K classes

def classify(series):
    # Normalise the 13 bands with a LayerNorm over channels
    # [code: src/models/TransformerEncoder.py#L43]
    x = layer_norm(series)                          # (B, T, C)

    # Project the bands to the hidden size with a 1x1 convolution
    # [code: src/models/TransformerEncoder.py#L46]
    x = linear(x)                                   # (B, T, D_h)
    x = layer_norm(x)                               # (B, T, D_h)

    # Add a fixed sinusoidal position code for positions 1..T
    # [§4.2] [code: src/models/TransformerEncoder.py#L50-L51]
    x = x + sinusoidal_table[1:T + 1]               # (B, T, D_h)

    # Stack L encoder blocks
    # [§4.2] [code: src/models/transformer/Models.py#L76-L78]
    for block in range(L):
        x = encoder_block(x)                        # (B, T, D_h)

    x = layer_norm(x)                               # (B, T, D_h)
    # Max over the time axis gives one vector per series
    # [§4.2] [code: src/models/TransformerEncoder.py#L60]
    pooled = max_over_time(x)                       # (B, D_h)
    # Linear layer without bias, then log-softmax
    # [code: src/models/TransformerEncoder.py#L62-L70]
    return log_softmax(linear_no_bias(pooled))      # (B, K)

One encoder block and one attention head, as implemented:

# One encoder block: attention, residual, then per-date dense layers  [§3.4]
# [code: src/models/transformer/Layers.py#L17-L25]
def encoder_block(x):
    # Multi-head self-attention: every head sees the same x  [eq. 4]
    mixed = multi_head_attention(x, x, x)          # (B, T, D_h)
    # Dropout, skip connection, then layer normalisation
    # [code: src/models/transformer/SubLayers.py#L59-L60]
    x = layer_norm(x + dropout(mixed))             # (B, T, D_h)
    # Two dense layers applied to each date separately, ReLU in between  [§4.2]
    # [code: src/models/transformer/SubLayers.py#L74-L81]
    hidden = relu(dense_1(x))                      # (B, T, 4 * D_h)
    x = layer_norm(x + dropout(dense_2(hidden)))   # (B, T, D_h)
    return x


# One attention head  [eq. 4] [§3.4]
# [code: src/models/transformer/Modules.py#L16-L28]
def attention_head(inputs, d_k):
    # Queries, keys and values are linear maps of the same series  [§3.4]
    queries = inputs @ w_query                     # (B, T, d_k)
    keys = inputs @ w_key                          # (B, T, d_k)
    values = inputs @ w_value                      # (B, T, d_v)
    # Scaled dates-by-dates scores; paper writes sqrt(D_h), code uses sqrt(d_k)
    scores = queries @ keys.transpose(1, 2)        # (B, T, T)
    scores = scores / sqrt(d_k)                    # (B, T, T)
    # Each date's weights over all dates sum to one  [code: Modules.py#L14]
    weights = softmax(scores, axis=-1)             # (B, T, T)
    return weights @ values                        # (B, T, d_v)

The learning-rate schedule used for the Transformer. The paper describes warm-up and then exponential decay §5; the code implements the schedule below instead:

# Learning rate at optimiser step `step` (counted from 1)  [§5]
# [code: src/utils/scheduled_optimizer.py#L11-L15] [#L32-L35] [#L40-L41]
def learning_rate(step, d_model, warmup_steps):
    # Base scale is d_model^-0.5, from the model width
    # [code: src/utils/scheduled_optimizer.py#L15]
    base = d_model ** -0.5                          # scalar
    # Linear warm-up, then inverse-square-root decay
    # [code: src/utils/scheduled_optimizer.py#L33-L35]
    scale = min(step ** -0.5, warmup_steps ** -1.5 * step)
    return base * scale                             # scalar

Analysis tools used in the paper:

• Temporal feature importance: the gradient of the largest predicted probability with respect to each input band and date. It needs no labels and works on any differentiable model §8.3 fig. 6 03_FeatureImportance.ipynb.
• Attention maps for the first of three self-attention layers, shown as matrices and as bipartite graphs over the input dates §8.4 fig. 7.
• A comparison of self-attention with the soft-attention used in DuPLO §8.5 fig. 8.
• t-SNE, UMAP and PCA embeddings of features after each layer §8.6 fig. 9.

Hyperparameters

The table lists everything the paper or code states for the Transformer and the shared training setup. "Not stated" means neither the paper nor the code fixes the value for the reported results.

Results

Author-reported. Setup: Sentinel-2 level-1C top-of-atmosphere reflectance, averaged per field parcel, from Google Earth Engine (raw) §7.2.1; an industry-preprocessed series from GAF AG (cloud-free dates, gap filling) for the same parcels §7.2.2. Three Bavarian regions (Hollfeld, Krumbach, Bavarian Forest) with 4500 m blocks split 4:1:1 into training, validation and evaluation blocks, 94, 21 and 29 blocks §7. Two class sets: 23 land-use classes and 12 land-cover classes §7.1. Models were fitted on the training and validation blocks and scored on the held-out evaluation blocks tab. 1 caption.

Source: tab. 1 (the table's subtables give kappa, accuracy and F1 for each class set). Values are mean ± standard deviation over the three hyperparameter sets. DuPLO was not tuned; the authors used its published configuration §4.5 §6 §8.1.

Key statements from the text, all author-reported:

• Preprocessing raised accuracy, F1 and kappa for all models §8.1.
• On raw data, RNN and Transformer were slightly better than MS-ResNet and TempCNN, most clearly on the 12-class set §8.1.
• Random forest reached competitive accuracy on preprocessed data but fell behind on raw data §8.1.
• The gradient analysis showed zero gradient on cloudy dates for RNN and Transformer, and non-zero gradient on some cloudy dates for the convolutional models §8.3.

Paper vs code

Code: marccoru/crop-type-mapping at commit fcc040181de93ee466e525150983ddd6822f86d5, MIT licence. The paper links this repository in a footnote §4.2, and the repository readme names it as the source code of the paper. The MIT copyright line names Marc Rußwurm.

Differences between the paper and the code:

• Learning rate ignored for the Transformer. The paper tunes a base learning rate for every model §5. For the Transformer the code builds an Adam optimiser with that rate, but the wrapper overwrites it before every step scheduled_optimizer.py [#L37-L44]. The tuned learning rate has no effect there tune.py.
• Decay shape. Paper: warm-up, then exponential decay §5. Code: linear warm-up, then inverse-square-root decay scheduled_optimizer.py.
• Adam settings. The paper gives β2 = 0.98 for "the RNN, MS-ResNet and TempCNN" §5. The code uses β = (0.9, 0.999) and ε = 1e-8 for those three train.py. The Transformer gets β = (0.9, 0.98) and ε = 1e-9 train.py.
• Attention scale. Paper: divide by the square root of D_h eq. 4 §3.4. Code: divide by the square root of the per-head key size d_k = D_h // heads Modules.py train.py. With 8 heads and D_h = 256, that is √32 in code and √256 in the formula.
• Search ranges in code. Paper: learning rate in [1e-8, 1e-1] and weight decay in [1e-12, 1e-1] §5. Code: hp.loguniform(-1, -8) and hp.loguniform(-1, -12) tune.py. Hyperopt's loguniform takes bounds in log space, so the code ranges are roughly 3.4e-4 to 0.37 for the learning rate and 6e-6 to 0.37 for the weight decay. This is derived from the hyperopt definition, not stated in the repository.
• Warm-up values. Paper: 10, 100, 1000 §5. Code: 1, 10, 100, 1000 tune.py.
• Tuning schedule. Paper: 60 epochs, four brackets, 10-epoch grace period, 34 random starting points §6. Code: 40 epochs, 2-epoch grace period, 4 starting points tune.py [#L205] [#L24] [#L299] [#L41].
• Input layers. The code normalises the 13 bands with a LayerNorm and projects them with a 1x1 convolution before the encoder, and normalises after the projection TransformerEncoder.py [#L43-L48]. The paper mentions none of this.
• Output layer. The final linear layer has no bias TransformerEncoder.py.
• Positional encoding. The code uses a fixed sinusoidal table, frozen, for positions 1 to T Models.py [#L72-L74] TransformerEncoder.py. The paper says it adds one and cites the original Transformer §4.2.
• Random forest class cap. The 500-sample-per-class cap appears in the RF and SVM scripts only random_forest.py, not in the neural-network training path train.py.
• Data. The preprocessed GAF AG series used in the paper is not public. The repository's GAF loader needs a dataset that download.sh does not fetch download.sh. The BavarianCrops raw data comes from a Zenodo archive download.sh.
• Hyperparameters chosen. The paper does not state the selected hyperparameters for each model. The code's fallback Transformer configuration is hidden 128, 3 heads, 3 layers, warm-up 1000, dropout 0.262, weight decay 0.000413 and learning rate 0.255, which is not the same as the paper's search hyperparameter.py.
• Early stopping uses the test partition. In the final-run configuration the test partition is test experiments.py, and the early-stopping check reads the loss logged in test mode trainer.py. The paper does not say which partition drives the stopping rule §5.

Pitfalls

• The learning-rate search is wasted for the Transformer. Tuning a base learning rate does nothing in the released code. Tune warm-up and width instead, and expect a peak learning rate of about D_h^-0.5 × warm-up^-0.5 from the schedule, for example about 2e-3 at D_h = 256 and warm-up 1000 (derived from the code formula).
• Sampling is random at every read. During training and testing, each item draws 70 dates without replacement and keeps them in time order BavarianCrops_Dataset.py. Test scores therefore depend on the random draw. The paper does not mention this sampling.
• Series must have at least 70 valid dates. NaN dates are dropped first BavarianCrops_Dataset.py. The draw without replacement then fails on shorter series. This follows from the code and is not discussed in the paper.
• Loss is unweighted. The code computes inverse-frequency class weights BavarianCrops_Dataset.py but passes none to the loss trainer.py. Implement a weighted loss only if you want to deviate from the reference.
• Head width does not match D_h when heads do not divide it. With D_h = 16 and 3 heads, the per-head size is 5, so the attention output is 15 wide before the output projection train.py SubLayers.py.
• Dropout placement. Dropout sits on the attention output, on the feed-forward output and on the attention weights at a fixed 0.1 SubLayers.py SubLayers.py Modules.py. The paper's tunable dropout is only one of these §6.
• Max-pooling over a fixed length. The pooling layer is sized to T = 70, so the model needs exactly that many dates TransformerEncoder.py.
• Paper-internal inconsistency on the tuning budget. §6 reports 2400 model evaluations, which equals 300 sets × 4 architectures × 2 datasets. §8.1 says the optimal hyperparameters were also chosen separately for the 23-class and 12-class sets, which would double that number §6 §8.1.
• Random forest is handicapped differently. The authors say the balanced 500-per-class subset favours the random forest §4.6. Compare its numbers with care.
• Gradient importance is interpretive. The gradient of the top probability with respect to inputs is zero where the model's output does not change with small input changes. The paper reads this as the model suppressing cloud dates §8.3. It is one example per class and one model run, so check it on your data before relying on it.
• Attention maps show one layer. Only the first of three attention layers is visualised §8.4. Attention weights are not a complete explanation.
• Cost. Self-attention scores are T × T per head, so memory grows with the square of the number of dates. The paper does not report run time or memory §8.
• Results are only comparable inside this study. The paper compares models on its own splits and class sets. The preprocessed GAF AG series is not public, so those rows cannot be reproduced from the repository. The raw BavarianCrops archive is public download.sh.

Lineage

• Builds on: the Transformer encoder of Vaswani et al. (arXiv:1706.03762, no card), which the paper adapts to time series §3.4 §4.2. Also the authors' earlier LSTM model for crop types, which the RNN follows §4.1.
• Compared against: TempCNN (Pelletier et al., no card) §4.4; MS-ResNet, a multi-scale residual 1-D network (Wang et al., no card) §4.3; DuPLO, a CNN, GRU and soft-attention model (Interdonato et al. 2019, no card) §4.5; and a random forest with spectral indices §4.6.
• Superseded by: not stated in the paper, and no later card is linked here.