Whittaker-Henderson smoother for long satellite image time series interpolation

Abstract

Whittaker smoother is a widely adopted solution to pre-process satellite image time series. Yet, two key limitations remain: the smoothing parameter must be tuned individually for each pixel, and the standard formulation assumes homoscedastic noise, imposing uniform smoothing across the temporal dimension. This paper addresses both limitations by casting the Whittaker smoother as a differentiable neural layer, in which the smoothing parameter is inferred by a neural network. The framework is further extended to handle heteroscedastic noise through a time-varying regularization, allowing the degree of smoothing to adapt locally along the time series. To enable large-scale processing, a sparse, memory-efficient, and fully differentiable implementation is proposed, exploiting the symmetric banded structure of the underlying linear system via Cholesky factorization. Benchmarks on GPU demonstrate that this implementation substantially outperforms standard dense linear solvers, both in speed and memory consumption. The approach is validated on SITS acquired over the French metropolitan territory between 2016 and 2024. Results confirm the feasibility of large-scale heteroscedastic Whittaker smoothing, though reconstruction differences with the homoscedastic baseline remain limited, suggesting that the transformer architecture used for smoothing parameter estimation may lack the temporal acuity needed to capture abrupt noise variations such as singleday cloud contamination.

Figures (4)

• Figure 1: Synthetic exemple of a time series with heteroscedastic noise. The dashed black curve is the true function, the light gray circles are the noisy observation and the blue/red curves are the reconstruction using the Whittaker smoother. The noise level is higher between time 6 and 10 and the Whittaker smoother with constant regularization value leads to noisy interpolation (blue curve), while the smoother with a time-varying regularization value leads to a smoother reconstruction (red curve). The right plot shows the regularization value a function of the time: for the heteroscedastic model, a piecewise constant function was used.
• Figure 2: Illustration of the band storage for $T=4$ and $(k+1)=1$ for $\boldsymbol{\Omega}$. The left side of the figure illustrates the conventional full matrix storage, while the right side shows the band storage format, which stores only the diagonal and lower diagonal elements. As for triangular matrix, upper and lower format exists. In this work, the lower format was used.
• Figure 3: Location of the 19 tiles over the French territory. Background map © OpenStreetMap contributors. List of tiles: 31TDJ, 31UDR, 31TCH, 31TGL, 31TGK, 31TCL, 31TDN, 31TFL, 30TYP, 31TGJ, 31TDH, 31UFP, 31UEQ, 31TEN, 31TEK, 31UEP, 30UVU, 30TYT, 30UXU.
• Figure 4: Reconstruction of NDVI SITS. In the upper plot, colored dots and triangles represent the raw acquisitions between 2016 and 2024, dot are valid dates and triangle are non valid dates. The black continuous curve represents the reconstruction obtained with the heteroscedastic model, while the shaded area indicates the credibility intervals calculated using equation (2.2) from biessy:hal-04124043. The dashed teal curve is the reconstruction using the homoscedastic model. In the lower plot, red and black curves represent the smoothing parameters for homoscedastic and heteroscedastic models, respectively.