Monotone Piecewise-Linear Harmonization of Multitemporal Multispectral Surface Reflectance With Unsupervised Deep Multiview Pseudoinvariant Features Learning

Lino Garda Denaro; Chia-Hsiang Lin

eess.IVOct 15, 2026

Monotone Piecewise-Linear Harmonization of Multitemporal Multispectral Surface Reflectance With Unsupervised Deep Multiview Pseudoinvariant Features Learning

Authors: Lino Garda Denaro, Chia-Hsiang Lin

Organizations: Research Center for Geoinformatics, National Research and Innovation Agency (BRIN), Cibinong, Indonesia · Department of Electrical Engineering, National Cheng Kung University (NCKU), Tainan, Taiwan

Abstract

Surface reflectance products derived from optical satellite may exhibit systematic inconsistencies across acquisition dates, because atmospheric correction is generally performed independently for each image. Such residual biases distort the temporal behavior of invariant ground targets and reduce the reliability of multitemporal surface reflectance analysis. This study proposes Deep Multitemporal Generalized Canonical Correlation Analysis (DeepMGCCA), an unsupervised deep multiview pseudoinvariant feature learning framework that integrates Deep Generalized Canonical Correlation Analysis (GCCA) with Multivariate Alteration Detection (MAD) to identify deep-pseudoinvariant features (DIFs) for the removal of surface reflectance biases across multitemporal satellite observations. Unlike conventional MAD-based approaches that rely on linear GCCA, the proposed framework replaces date-specific linear projections with a shared convolutional encoder. The framework learns deep invariant features (DIFs) directly from pixel-weighted canonical variates via iterative reweighting strategy, while a reconstruction regularizer preserves surface reflectance information and prevents degenerate collapse from suppressing genuine temporal variability. Experiments were conducted on four annual Landsat-8/9 and Sentinel-2 acquisitions across four regions of interest (ROIs) to evaluate four DIF-based harmonization transforms. A monotone piecewise-linear transform derived from DIFs quantiles consistently achieved the largest bias reduction, lowering the between-date spectral angle by up to 2.54 and RMSE by 0.015, and raising PSNR by up to 4.04~dB. These findings demonstrate the potential of deep invariant feature learning for multitemporal surface reflectance bias removal while highlighting the importance of preserving physically meaningful temporal and spectral information during nonlinear representation learning.

Figures & tables

Fig. 1: Overview of the proposed DeepMGCCA framework and harmonization pipeline. Top row: each input image $\mathbf{X}_{j}$ , $j=1,\ldots,J$ , is mapped to a learned feature representation $\mathbf{Z}_{j}=h_{\phi_{j}}(f_{\theta}(\mathbf{X}_{j}))$ through the shared encoder $f_{\theta}$ followed by a date-specific head output $h_{\phi_{j}}$ , trained end-to-end by backpropagating the multiview correlation and reconstruction objectives. The consensus-aligned canonical variates $\mathbf{U}_{j}$ are then differenced across the connected date pairs ( $C_{jk}=1$ , $j\neq k$ ) and aggregated into the normalized MAD statistic $\mathrm{NMAD}_{j}$ . Bottom row: a dynamic threshold converts $\mathrm{NMAD}_{j}$ into the unchanged map $\mathbf{X}_{\text{mask}_{j}}$ , from which the DIF of all dates are extracted $\mathbf{X}_{\text{DIF}_{j}}$ . The DIF quantiles determine the knot grid of a monotone piecewise-linear transform $\mathrm{T}(\cdot)$ , one per date and fitted over all pairs $j<k$ , which is applied to every pixel of the scene to produce the harmonized images $\mathbf{X}^{\prime}_{1},\ldots,\mathbf{X}^{\prime}_{J}$ .

Fig. 2: DeepMGCCA architecture. Left: for each date $j=1,\ldots,J$ , the input patch $\mathbf{X}_{j}\in\mathbb{R}^{B\times H\times W}$ passes through a shared nonlinear encoder $f_{\theta}$ with weights tied across all $J$ dates, followed by a date-specific head $h_{\phi_{j}}$ that produces the head outputs $\mathbf{Z}_{j}=h_{\phi_{j}}(f_{\theta}(\mathbf{X}_{j}))\in\mathbb{R}^{r\times H\times W}$ , from which the canonical variates $\mathbf{U}_{j}$ are formed by consensus alignment after training. The weighted unit-variance constraint is enforced by weighted center-scaling, and a date-specific decoder $g_{\psi_{j}}$ reconstructs $\mathbf{X}_{j}$ from $\mathbf{Z}_{j}$ , contributing the reconstruction regularizer. Right: the multiview correlation objective and the reconstruction regularizer are jointly optimized over $\theta$ , $\{\phi_{j}\}_{j=1}^{J}$ , and $\{\psi_{j}\}_{j=1}^{J}$ . The trained canonical variates $\mathbf{U}_{j}$ are fed into the change-detection flow: pairwise difference variates $\widetilde{\mathbf{MAD}}_{jk}=\mathbf{U}_{j}-\mathbf{U}_{k}$ are aggregated into the normalized MAD (NMAD) statistic and thresholded to produce the final deep-pesudoinvariant-feature (DIF) image.

Fig. 3: Piecewise constrained harmonization. Each date $j$ (colored curves) is mapped by a monotone piecewise-linear transform $T_{j,b}$ defined on a knot grid $t_{b,0}<\cdots<t_{b,K}$ (dashed vertical lines) placed at the quantiles of the pooled DIF population. The agreement term $w_{1}$ ( $E_{C}$ , Eq. ( 22 )) drives the transformed quantile functions of all date pairs toward a common level. The anchor term $w_{2}$ ( $E_{T}$ ) pulls each transform toward the identity map (dotted diagonal). The band-smoothness term $w_{3}$ ( $E_{S}$ ) couples the hinge coefficients across adjacent bands, and and the curvature term $w_{4}$ ( $E_{K}$ , dashed ellipses) penalizes kinks between consecutive knots. Monotonicity is imposed as a box constraint on the hinge coefficients (Eq. ( 21 )).

Panel of Table 4 — TABLE I: Regions of interest (ROIs), acquisition dates, sensors, spatial resolutions, and environmental descriptions.

	Landsat 8/9 OLI			Sentinel-2 MSI
Band	No.	Wavelength (nm)		No.	Wavelength (nm)
	No.	Range	Center	No.	Range	Center
B1 (Blue)	B2	450–510	480	B2	458–523	490
B2 (Green)	B3	530–590	560	B3	543–578	560
B3 (Red)	B4	640–670	655	B4	650–680	665
B4 (NIR)	B5	850–880	865	B8	785–900	842

TABLE II: Spectral bands used in this study. Landsat 8/9 OLI and Sentinel-2 MSI band numbers, wavelength ranges, and central wavelengths.

Fig. 4: Qualitative harmonization results for the four ROIs of Table I . For each ROI, the four multitemporal acquisitions are mosaicked into false-/true-color composites: ( top ) the original surface-reflectance images, in which interdate brightness and color seams are clearly visible across acquisition dates, and ( bottom ) the harmonized images produced by the proposed transform, in which the interdate drift is removed while genuine spatial detail is preserved.

Fig. 5: Row-wise mean reflectance profiles over the central column window (highlighted as the translucent band in the composites), computed for the four ROIs of Table I from the four multitemporal acquisitions of each ROI ( left : original and right : harmonized by the proposed transform). For each scene, the per-row mean reflectance (horizontal axis) is plotted against the image row (vertical axis, top to bottom), with one curve per acquisition date and spectral band. Before harmonization, the four date-profiles appear as vertically shifted copies of one another, a direct visualization of the interdate reflectance inconsistency (additive and multiplicative bias) on invariant terrain. After harmonization, the profiles of the four dates seamed, indicating that the interdate drift is removed without suppressing genuine spatial or temporal detail.

ROI	Variant	RMSE $\downarrow$	PSNR $\uparrow$	SSIM $\uparrow$	SAM h $\downarrow$	SAM sc $\downarrow$	$\rho\uparrow$
Taipei	Orig.	0.0359	31.91	0.7568	5.445	—	0.7416
	RA [ 21 ]	0.0217	35.33	0.7894	2.512	4.872	0.7407
	SSA [ 29 ]	0.0222	35.26	0.7882	2.522	4.926	0.7407
	GSC [ 30 ]	0.0238	34.84	0.8065	3.630	2.262	0.7409
	Proposed	0.0210	35.81	0.8184	2.907	2.837	0.7106
Saudi Arabia	Orig.	0.0237	34.17	0.9230	3.066	—	0.9499

TABLE III: Per-ROI quantitative comparison of harmonization variants for the proposed DeepMGCCA. Values are means over all date pairs and spectral bands. RMSE, SAM h , and SAM sc , lower indicates better result. PSNR, SSIM, and Pearson’s

\rho

, higher indicates better result. SAM h is the mean spectral angle between the harmonized dates over the DIF population. SAM sc is the mean spectral angle between each image and its own corrected version, undefined for the unprocessed images (—). Orig. is unprocessed. RA is reference-affine [ 21 ] . SSA is sequential stepwise affine [ 29 ] . GSC is globally spectral constrained [ 30 ] . Proposed is the proposed monotone piecewise-linear transform. Best value per column per ROI in bold. in the SAM sc column both GSC and Proposed are highlighted, as the two operate on a deliberate trade-off between spectral self-preservation and between-date agreement (see text). All variants are scored on the same DIF mask.

Remote Sensing

Remote sensing uses satellite, aerial, and ground-based observations to measure and map Earth and other planetary surfaces.

See all 527 papers

HLS-GPT: A Generative Pretrained Transformer (GPT) for Continental-Scale NASA Harmonized Landsat and Sentinel-2 (HLS) Reflectance Reconstruction Across All Bands on Arbitrary Dates
Jun 16, 2026Junjie Li, Hankui K. Zhang, David P. RoyTransformer Remote Sensing
Geospatial Sciences Center of Excellence, Department of Geography and Geospatial Sciences, South Dakota State University, Brookings, SD 57007, USA · Department of Geography, Environment, and Spatial Sciences, & Center for Global Change and Earth Observations, Michigan State University, MI 48824, USA
Robust RPC Bundle Adjustment for Multi-Date Satellite Imagery with Season-Invariant Correspondences
Jul 29, 2026Roger Marí, Elías Masquil, Xavier Bou +2Remote Sensing Image Understanding Multi-View Geometry
Eurecat, Multimedia Technologies, Barcelona, Spain · IIE, Facultad de Ingeniería, Universidad de la República, Uruguay · AMIAD, France +2
LC-TSRNet: A Land-Cover-Guided Two-Stage Network for Annual MODIS Surface Reflectance Time Series Reconstruction
Oct 18, 2026Wenjuan Zhang, Xuesong Ma, Zhen Li +5Remote Sensing Multiscale Time-Series Modeling
State Key Laboratory of Remote Sensing and Digital Earth, Aerospace Information Research Institute, Chinese Academy of Sciences, Beijing 100101, China · Aerospace Information Research Institute, Chinese Academy of Sciences, Beijing 100094, China · College of Resources and Environment, University of Chinese Academy of Sciences, Beijing 100049, China +3

Jun 16, 2026J/K move · Enter open · S save

HLS-GPT: A Generative Pretrained Transformer (GPT) for Continental-Scale NASA Harmonized Landsat and Sentinel-2 (HLS) Reflectance Reconstruction Across All Bands on Arbitrary Dates

Junjie Li, Hankui K. Zhang, David P. Roy

Geospatial Sciences Center of Excellence, Department of Geography and Geospatial Sciences, South Dakota State University, Brookings, SD 57007, USA · Department of Geography, Environment, and Spatial Sciences, & Center for Global Change and Earth Observations, Michigan State University, MI 48824, USA

Transformer Remote Sensing

Open paper AskPDF arXiv

Recent deep learning methods for Landsat and Sentinel-2 reflectance time series reconstruction remain limited by restricted spectral coverage, limited geographic scalability, or patch-based designs with short temporal contexts. We present HLS-GPT, a large-scale generative pretrained Transformer model for reconstructing NASA Harmonized Landsat Sentinel-2 30 m surface reflectance for all bands, any date, and any pixel location. HLS-GPT uses a hierarchical Transformer architecture to handle the different spectral band configurations of Landsat and Sentinel-2 and operates on single-pixel 12-month time series. To capture geographic and seasonal variability, the model was trained with nine years of HLS time series from more than 0.25 million training pixels across the conterminous United States. A random cropping and masking strategy extracts 12-month periods with varying start dates across epochs, masks 50% of valid observations, and trains the model to reconstruct the masked reflectance values from the remaining observations. Evaluation using more than 62,000 independent test pixels shows robust reconstruction under diverse land surface conditions, including complex crop phenology and sparse, irregular observations. Leave-one-observation-out evaluation achieved reconstruction RMSE below 0.026 for all HLS spectral bands, with relative RMSE below 35% for visible bands and below 13% for other bands. Red-edge band errors were comparable to red and near-infrared errors despite the absence of red-edge bands on Landsat. Sensitivity analyses that randomly masked 10% to 90% of test observations showed only modest degradation when 10% to 50% of observations were masked, with all-band RMSE below 0.028. Image reconstruction over nine independent 109 by 109 km CONUS HLS tiles further demonstrates that HLS-GPT outperforms two conventional methods and the NASA-IBM Prithvi model.

Abstract

Remote Sensing

Explore similar work

HLS-GPT: A Generative Pretrained Transformer (GPT) for Continental-Scale NASA Harmonized Landsat and Sentinel-2 (HLS) Reflectance Reconstruction Across All Bands on Arbitrary Dates

Robust RPC Bundle Adjustment for Multi-Date Satellite Imagery with Season-Invariant Correspondences

LC-TSRNet: A Land-Cover-Guided Two-Stage Network for Annual MODIS Surface Reflectance Time Series Reconstruction