Reconstruction of seamless harmonized Landsat Sentinel-2 (HLS) time series via self-supervised learning

The Harmonized Landsat Sentinel-2 (HLS) data, harmonizing Landsat-8/9 and Sentinel-2 imagery, offers frequent 30 m resolution multispectral observations but is often contaminated by clouds, shadows, and snow that reduce the availability of good-quality surface observations. Traditional techniques for reconstructing HLS time series, such as polynomial, logistic, or harmonic functions that model seasonal reflectance changes struggle with complex changes violating the function fitting assumptions. We propose a data-driven time series reconstruction framework based on Transformer, termed self-supervised learning for interpolation (SSLI) with a smoothing constraint to model seasonal reflectance change patterns without any annotated labels. In this study, a year of HLS 30 m data were processed into 3-day surface reflectance composites (i.e., 122 composites). SSLI was trained by using randomly selected 70% of the good-quality 3-day composites in a time series to reconstruct the reflectance for the remaining 30%. The random masking was undertaken independently for each HLS pixel time series and for each training iteration so that the selected composite periods cover all the good-quality periods evenly. The methodology was tested on five diverse regions in the Conterminous United States (CONUS) each using three HLS tiles. Two versions of SSLI were evaluated. SSLI i was trained using time series samples from one tile per region to impose data independence, and SSLI ii was trained using samples from all the tiles to simulate data availability in real-world applications. The results were compared with those of three state-of-the practice approaches for gap-filling reflectance and Normalized difference vegetation index (NDVI) time series, i.e., the fill-and-fit (FF), the dynamic temporal smoothing (DTS), and the double logistic (DL) algorithms. The superior performance of SSLI, reflected in lower RMSE (SSLI i: 0.0192, SSLI ii: 0.0164) and higher R2 (SSLI i: 0.9150, SSLI ii: 0.9349) compared to the other algorithms, demonstrates much higher accuracy in reconstructing complex phenological changes with multiple greenness peaks (e.g., in cropland) and robustness to temporal cadence variations and time series noise. The potential of adapting SSLI for land cover mapping and global-scale time series reconstruction is discussed. The developed codes and training data were made publicly available.