Abstract: Objectives: Multipath effects are a primary unmodeled error source constraining the precision of global navigation satellite system (GNSS) data processing. Sidereal filtering (SF) relies on strict daily orbital repeatability, and cumulative orbit period drift causes temporal misalignment between historical templates and current observations, leading to performance degradation or negative corrections. Multipath hemispherical maps (MHM) apply spatial grid averaging indexed by elevation and azimuth, but spatial discretization irreversibly smooths high-frequency temporal variations, producing step-like correction artifacts. Neither method addresses dynamic multipath from moving reflectors, which lacks daily repeatability. Methods: This paper proposes a robust two-stage machine learning framework, spatio-temporal XGBoost (ST-XGB), which explicitly captures the spatio-temporal heterogeneity of multipath errors in short-baseline environments, and decomposes multipath errors into a geometric trend component and a local time-varying component. In the first stage, it uses six spatial geometric features including elevation, azimuth, and their four trigonometric transformations, to capture long-period systematic geometric biases. In the second stage, it constructs a residual feature set from the residuals to model short-term dynamic fluctuations. All temporal features use only historical data prior to the current epoch to prevent data leakage, and a 5-day sliding window training strategy balances sample diversity and environmental stationarity. Ablation experiments across five progressively designed configurations and Shapley additive explanations (SHAP) analysis reveal the feature masking effect in single-stage models that lag features dominate learning and suppress spatial feature contributions. The two-stage strategy resolves it by providing independent learning spaces for each feature type. The framework is evaluated on 30-day 1 Hz GPS data from five short baselines at Curtin University (DOY 244–273, 2021) and validated on a TU Delft dataset with different hardware. Results: Ablation experiments show the two-stage configuration incorporating complete residual time-varying features achieves an average residual reduction rate of 28.04%, outperforming all single-stage configurations. Across five baselines and 25 test days, ST-XGB reduces P1, P2, L1, and L2 residual RMS by 26.09%, 38.90%, 20.63%, and 20.49%, respectively, consistently outperforming SF and MHM. Under a fixed model without daily updates, ST-XGB maintains 23.19% overall average residual reduction over 25 days, while SF degrades to 8.41% with negative corrections appearing from day 10 due to orbit drift accumulation. Simulated kinematic positioning yields RMS values of 1.5 mm, 1.8 mm, and 4.4 mm in East, North, and Up, representing improvements of 25.0%, 21.7%, and 20.0% over the uncorrected solution. Independent validation at TU Delft confirms hardware robustness, with positioning improvements of 47.6%, 45.7%, and 50.0%. Dynamic multipath simulation tests show ST-XGB maintains 26.82% average residual reduction under severe interference, versus 13.10% for SF and 16.16% for MHM. Conclusions: ST-XGB resolves feature coupling interference in singlestage models through two-stage decoupling, enabling geometric trend modeling and time-varying feature representation to be fully exploited in independent learning spaces. The method demonstrates consistent robustness across hardware configurations, observation environments, and dynamic interference conditions. Current experiments are limited to GPS single-system short-baseline scenarios. Future work will extend the framework to multi-GNSS (GPS/BDS/Galileo/GLONASS) positioning and explore deep learning architectures (LSTM, Transformer) for enhanced spatio-temporal representation of multipath error sequences.