# Fusion of Multitemporal Multisensor Velocities Using Temporal Closure of Fractions of Displacements

Laurane Charrier, Yajing Yan, Emmanuel Trouvé, Elise Colin Koeniguer,

Jérémie Mouginot, Romain Millan

► To cite this version:

Laurane Charrier, Yajing Yan, Emmanuel Trouvé, Elise Colin Koeniguer, Jérémie Mouginot, et al. "Fusion of Multitemporal Multisensor Velocities Using Temporal Closure of Fractions of Displacements." IEEE Geoscience and Remote Sensing Letters, 2022, 19, pp.2001805. 10.1109/LGRS.2022.3227413. hal-04030559

HAL Id: hal-04030559

https://hal.science/hal-04030559

Submitted on 15 Mar 2023

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Fusion of multi-temporal multi-sensor velocities using temporal closure of fractions of displacements

Laurane Charrier, *Student Member, IEEE*, Yajing Yan, *Member, IEEE*, Emmanuel Trouvé, *Senior Member, IEEE*
Elise Colin Koeniguer, Jérémie Mouginot, and Romain Millan

***Abstract***—Numerous glacier velocity observations, derived from spaceborne imagery, are available online, but it remains difficult to analyze them because they are measured with different temporal baselines, by various sensors. In this study, we propose a novel formulation of the temporal closure to fuse multi-temporal multi-sensor velocity observations without prior information on the displacement behavior and the data uncertainty. We establish a system of linear equations between combinations of displacement observations and fractions of estimated displacements. The proposed approach provides a velocity time-series with a regular and optimal temporal sampling, the latter representing a compromise between the temporal resolution and the signal-to-noise ratio. The proposed approach is first evaluated on synthetic datasets and second on Sentinel-2 and Venus velocity observations over the Fox glacier in New Zealand. The results show the intra-annual variability of Fox glacier surface velocity with a reduced uncertainty and complete temporal coverage.

***Index Terms***—velocity time series, multi-sensor fusion, temporal inversion, glacier, intra-annual variations

## I. INTRODUCTION

While numerous glacier velocity observations are available online [1]–[4] with a temporal resolution up to 2 days and a spatial sampling up to 50 m [4], [5], intra-annual variations of glacier velocity remain poorly understood at regional scales, especially in mountain areas. This is due to the difficulty to analyze scene-pair velocity observations since they span different temporal baselines [6], are based on images from different sensors, and still contain noise, artifacts, and missing data when outliers have been filtered out. Therefore, a question arises: how to exploit the amount of available multi-sensor and multi-temporal glacier velocities in order to study intra-annual velocity variations with a reasonable uncertainty at a global scale?

A way to fuse multi-temporal velocity observations is to use a Temporal Inversion (TI) based on the temporal closure of the displacement observation networks, inspired from the Small BAseline Subset (SBAS) in the InSAR community. Currently, most of the TI approaches consider mono-sensor datasets [6]–[9] or apply TI separately to each sensor datasets [10], [11]. A few studies imply multi-sensor datasets inside the same system of linear equations [12], [13] whereas it could improve the results by adding more redundancy as stressed

by [11]. In [12], [13] irregular Leap Frog (LF) time series (time series of displacements between consecutive dates) are inverted. This kind of time series are adapted to reconstruct a Cumulative Displacement (CD) time series afterward. While CD time series are useful to study slow-moving landslides or dunes [9], [10], [12], it is more adapted to convert the LF displacement time series into LF velocity time series for fast-moving targets with strong seasonality patterns, such as glaciers [8]. Because velocities correspond to the average of instantaneous velocity over the temporal baseline, the temporal sampling of the velocity series should be regular (all the LF velocities should span the same temporal baseline). Hence, [6] proposed to include combinations of velocity observations in the TI to obtain a Regular LF (RLF) with a given temporal sampling. This approach, which will be referred as Temporal Inversion using Combinations of displacements (TICO), allows obtaining a RLF with a temporal sampling matching a compromise between temporal resolution and signal-over-noise ratio.

This study is an extension of [6] in order to obtain a RLF with an optimal temporal sampling by taking advantage of every multi-sensor observation, even if outliers have been removed beforehand. For multi-sensor datasets, each sensor does not acquire images on the same temporal grid: hence, there is less redundancy and more disconnected subsets of observed displacements. This is even more true when outliers have been removed beforehand. Therefore, we introduce fractions of displacement in the TICO formulation proposed in [6]. This extension is called Temporal Inversion using Combinations of Fractions of displacements (TICOF).

Besides, evaluation of observation uncertainties remain crucial to invert velocity time series [6], [8], [12]. The shape of the similarity function used in image matching appear to be a robust estimation of precision for both off- and on-glacier pixels [6], [12], [14], [15]. However, most of the online datasets do not provide the full similarity function values [1]–[4]. Hence, there is a need to get a robust estimation of uncertainty without any information on surface correlation values.

Therefore, this study has two objectives: 1) to take into account that velocity observations have different uncertainties which are not necessarily provided by the authors 2) to propose a new system to invert a RLF time series with an optimal temporal sampling from multi-sensor datasets, which can present missing data. We apply the proposed TICOF approach to Sentinel-2 and Venus velocity observations over the Fox glacier located in the Southern Alps of New Zealand.

L. Charrier, Y. Yan, Emmanuel Trouvé are with LISTIC, Université Savoie Mont Blanc, Annecy, France. E. Colin Koeniguer and L. Charrier are with DTIS, ONERA, Paris, France. J. Mouginot is with IGE, Université Grenoble Alpes, Grenoble, France. R. Millan is with Department of Geosciences and Natural Resource Management, University of Copenhagen, Copenhagen K, Denmark.
## II. STUDY AREA AND DATA

The study focuses on the Fox glacier, in the Southern Alps of New Zealand. This glacier is a fast-flowing glacier with a strong seasonality [4], [16], situated under a maritime climate (Fig. S1). The considered velocity observations are obtained from Sentinel-2 and Venus optical images. Sentinel-2 has a repeat cycle of 5 days and Venus 2 days. Displacement offsets have been computed by [4] using a modified version of the cross-correlation algorithm AMPCOR from ROI_PAC [4], [5]. Outliers have been removed by rejecting offsets that deviate more than three pixels from the median offset computed over a spatial window of 9 × 9 pixels [4]. Velocity observations span temporal baselines from 5 to 100 days and from 330 to 400 days. The spatial sampling of both Sentinel-2 and Venus velocity maps is 50 m. The central dates of Sentinel-2 velocity observations range from July 21, 2016 to October 16, 2018, whereas those of Venus velocity observations cover the period of January 26, 2018 to May 21, 2019. The mean percentage of remaining values after outliers removal is 26 % for Sentinel-2 and 14 % for Venus as illustrated in Fig. S2.

## III. METHOD

### A. Data quality indicators

In many cases, uncertainties are not provided with velocity observations. In order to quantify the quality of velocity observations without any information given by the authors, two criteria are defined. The first one relies on the spatio-temporal coherence of velocity directions, the second on the spatio-temporal coherence of respectively x- and y- velocity component values [17]:

1) The cosine of the angle between each displacement vector and the spatio-temporal median [14] vector, called Median Angle (MA) here, is defined as:

$$ \cos(\vec{u}_{i,j,t} \cdot \vec{v}_{i,j}^w) = \frac{\vec{u}_{i,j,t} \cdot \vec{v}_{i,j}^w}{\|\vec{u}_{i,j,t}\| \|\vec{v}_{i,j}^w\|} \quad (1) $$

where $\vec{u}_{i,j,t}$ is the displacement vector at the pixel $(i, j)$ and time $t$ and $\vec{v}_{i,j}^w$ is the median vector in time over a spatial window $w$, here set to $3 \times 3$ pixels, centered on the pixel $(i, j)$. If the two vectors are co-linear, as expected, this criterion will be equal to 1. The negative values are set to 0 to penalize vectors which have a direction from 90° to 180° away from the median direction.

2) The modified z-score, a standardized and robust score, measuring the outlier strength [18], is computed respectively for the x and y component of the velocity vector as:

$$ \text{MZscore} = \frac{x_{i,j,t} - \tilde{x}_{i,j}^w}{1.483 \text{ MAD}} \quad (2) $$

where $x_{i,j,t}$ stands for the x or y component of the velocity at the pixel $(i, j)$ and time $t$, $\tilde{x}_{i,j}^w$ is the median of the x or y component of the velocity in time over a spatial window $w$, here set to $3 \times 3$ pixels, centered on the pixel $(i, j)$. MAD is the Median Absolute Deviation defined as: $\text{median}|x_{i,j,t} - \tilde{x}_{i,j}^w|$.

Note that the MZscore considers a robust estimate of the dispersion: 1.483 MAD, unlike the standard deviation which is not robust to outliers [1], [4].

The final confidence value of the two velocity components is the same for a given pixel and a couple of dates. It is computed as the multiplication of the three criteria scaled between 0 and 1: the scaled MZscore for x and y components of the velocity and the MA. This confidence value will be used in the inversion as an initial weight since it indicates data quality.

### B. Temporal Inversion using Combinations of Fractions of displacement (TICOF)

<figure><img src="image_2.png" /><figcaption>Fig. 1. Illustration of a) the traditional TI formulation b) the proposed TICOF formulation. The displacement observations from satellite 1 is in blue while the displacement observations from satellite 2 is in orange.</figcaption></figure>

The classical TI formulation is $AX = Y$, where $Y$ is a vector of dimension $n$ containing the displacement observations $d_{t_i, t_j}$ between dates $t_i$ and $t_j$, $X$ is a vector of dimension $p$ representing the LF displacement estimation $\hat{d}_{\tau_k, \tau_{k+\Delta\tau}}$ between dates $\tau_k$ and $\tau_{k+\Delta\tau}$. $\Delta\tau$ days is the temporal sampling which is identical for every LF displacement. $A$ is the design matrix of dimension $n \times p$ which links the vectors $X$ and $Y$. As shown in [6], the classical TI does not enable the user to obtain a LF time series using all available displacement observations if the temporal sampling $\Delta\tau$ is larger than the minimal time interval between two observed dates. For example, in Fig. 1, the minimal time interval between two observed dates is 1 day. For a temporal sampling of $\Delta\tau = 5$ days, $d_{t_0, t_2}$, $d_{t_2, t_4}$, $d_{t_2, t_{10}}$ and $d_{t_6, t_{10}}$ (in orange) cannot be written as any combinations of $\hat{d}_{\tau_k, \tau_{k+\Delta\tau}}$.

However, it can be established that $\hat{d}_{t_0, t_5} + \hat{d}_{t_5, t_{10}} = d_{t_0, t_2} + d_{t_2, t_{10}}$. Therefore, [6] proposed a novel TI formulation, namely Temporal Inversion using Combinations of displacement observations (TICO), to include combinations of displacement observations. The TICO is formulated as $BY = A'X$ where $B$ is a design matrix of dimension $n' \times n$ and $A'$ is the design matrix of dimension $n' \times p$, with $n'$ the number of displacement combinations as specified in [6].

The TICO formulation is adapted when a high amount of observed dates in $Y$ match the dates $\tau_k$ or $\tau_{k+\Delta\tau}$ because it makes possible to link a combination of $\hat{d}_{\tau_k, \tau_{k+\Delta\tau}}$ with a combination of $d_{t_i, t_j}$. However, when the displacement observations are based on different sensors, the observed dates may not be on the same temporal grid. The dates where acquisitions of both satellites exist depend on the common multiples of the repeat cycles of satellites 1 and 2, here $\Delta t_1$ and $\Delta t_2$. There
are even less common acquisition dates between sensors when missing data exist due to outliers removal.

Therefore, we propose an extension of TICO [6] for datasets containing multi-sensor displacements and missing data. When a displacement observation can not be included in the TICO formulation, a linear relation between $d_{t_i, t_j}$ and a fraction of displacements $\hat{d}_{\tau_k, \tau_{k+\Delta\tau}}$ is established, assuming linear displacement in the considered temporal interval. For example, in Fig. 1, $\frac{2}{5}\hat{d}_{t_0, t_5} = d_{t_2, t_4}$. In a more general form, $d_{t_i, t_j}$ is written as:

$$d_{t_i, t_j} = \sum_k w_k \hat{d}_{\tau_k, \tau_{k+\Delta\tau}} \quad (3)$$

$$\text{with } w_k = \begin{cases} 0, & \text{if } \tau_{k+\Delta\tau} < t_i \\ 0, & \text{if } \tau_k > t_j \\ 1, & \text{if } \tau_k \ge t_i \land \tau_{k+\Delta\tau} \le t_j \\ \frac{\tau_{k+\Delta\tau} - t_i}{\tau_{k+\Delta\tau} - \tau_k}, & \text{if } \tau_k < t_i < \tau_{k+\Delta\tau} \land \tau_{k+\Delta\tau} < t_j \\ \frac{t_j - \tau_k}{\tau_{k+\Delta\tau} - \tau_k}, & \text{if } \tau_k < t_j < \tau_{k+\Delta\tau} \land t_i < \tau_k \\ \frac{t_j - t_i}{\tau_{k+\Delta\tau} - \tau_k}, & \text{if } t_i > \tau_k \land t_j < \tau_{k+\Delta\tau} \end{cases}$$

The TICOF formulation is: $BY = A'X$ where $A'$ is filled, in consequence, with 0, 1 and/or fractions (equation S.1). Since the velocity observations have different but not necessarily known uncertainties, the inversion of $BY = A'X$ is obtained by an Iterative Reweighted Least Square (IRLS) using the Tukey's biweight function and a Tikhonov regularization [6]. The weight at the first iteration corresponds to the confidence indicators presented in section III-A.

The use of fractions of displacements is also found in [19] to adjust the first and last date of descending and ascending displacements. This idea assumes a constant velocity over the relevant temporal sampling. For large temporal sampling and/or non-linear displacement behavior, a bias may be introduced in the final estimation of the LF time series. To our knowledge, this bias has never been discussed in the literature.

## IV. SYNTHETIC DATASET

To assess the bias introduced by the use of fraction in the case of non-linear displacements (e.g. seasonal variation), synthetic simulations are performed. The synthetic instantaneous velocity is taken as: $v(t) = a + b \sin(\frac{2\pi}{T}t) + c \cos(\frac{2\pi}{T}t)$ with $T = 365.25$ as in [20]. $N_1$ displacements based on Sentinel-2 repeat cycle and $N_2$ displacements based on Venus repeat cycle are simulated. To make sure that the coefficients $a, b$ and $c$ represent well the data, instead of an arbitrary choice, these coefficients are estimated by an IRLS inversion by adding a regularization term corresponding to a displacement model with the coefficients $a, b$ and $c$ as parameters as performed in [21]. The corresponding system of equations is given in equation S.2. The system is solved for Sentinel-2 data on the point represented in blue in Fig S1. The coefficients are found to be: $a = -0.49, b = -0.0788$ and $c = 0.018$.

In Fig 2, LF time series obtained from a moving average of $\Delta\tau$ days (Fig. 2b), the traditional TI (Fig. 2c) and the TICOF formulation (Fig. 2d) are presented. The moving average corresponds to the average of all the velocities which

overlap each period of $\Delta\tau$ days. The results are represented for a temporal sampling of 20 days. Because velocities are the average of the instantaneous velocity over the temporal baseline, long temporal baseline velocities tend to pull the moving average results toward the annual averaged velocity, here equal to $a * 365.25 = -178.85$ m/y. In addition, because no displacement observations can be included in the traditional TI when only Venus data are available, the LF displacements in Fig. 2d) are equal to 0 after October 2018.

Moreover, the Root Mean Square Error (RMSE) between the synthesised and the estimated velocities is more than 80 m/y lower with the TICOF formulation than with the TI or TICO formulation, and about 20 m/y with TICOF than with the moving average, for different temporal sampling (Fig. 3). However, the mean residual of the inversion (the mean difference between AX and Y or A'X and BY) characterizes the bias introduced by the assumption of a constant velocity over the considered temporal sampling. For a synthetic dataset without noise, the residual of the inversion should be null. This is the case for the traditional TI or TICO (Fig. S3). On the contrary, the increase of the TICOF residual which ranges from 0 to 3 m/y, for 0 to 120 days temporal sampling, is caused by the bias introduced by the fractions. This bias is around 25 times lower than the benefit of using TICOF rather than TICO, in case of multi-sensor datasets (Fig. 3). The use of TICOF is advantageous: 1) for multi-sensor datasets and/or 2) datasets with a lot of missing data and/or 3) when the temporal sampling is not a multiple of the satellite repeat cycle (Fig. 3 and S4). In the absence of any of these conditions, Fig. S4 illustrates that the RMSE from TICO is up to 5 m/y better than the one from TICOF. This error also corresponds to the bias introduced by the fractions.

## V. REAL DATASET

### A. Uncertainty of the results for different temporal samplings

To analyze the uncertainty of the LF time series according to their temporal sampling, the RMSE is computed over the study area delimited by the orange square in Fig. S1, assuming that the true velocity magnitude is null over the stable areas. RMSE is performed over ice-free areas according to the Randolph Glacier Inventory (RGI) V6.0.

Fig. 4 reveals that the RMSE is lower after inversion than before inversion for a temporal sampling ranging from 5 to 20 days (the RMSE is 68% lower for a 5-days temporal sampling to 13% lower for a 20 days temporal sampling) and remains slightly higher with larger temporal sampling (around 12 m/y higher). According to the RMSE, the performance of the TICOF inversion is degraded with a temporal sampling larger than 20 days. It may be caused by the validity of the assumption beyond the use of fractions, which increases with the temporal sampling. Hence, the use of TICOF inversion appears beneficial for temporal samplings, here, lower than 20 days. From 5 to 20 days, the lowest RMSE is reached for a temporal sampling of 20 days, hence 20 days is selected to match a compromise between uncertainty and temporal resolution.
<figure><img src="image_3.png" /><figcaption>Fig. 2. a) Synthetic velocities spanning a given temporal baseline are represented by blue and orange bars, for velocities based on Sentinel-2 and Venus images respectively. These velocities correspond to the average of the instantaneous velocity (represented by a black curve) over the considered temporal baseline. Hence, small temporal baseline velocities are closer to the instantaneous velocity than long temporal baseline velocities. The dots represent the central date of each velocity. RLF time series obtained using the b) moving average c) traditional TI formulation d) TICOF formulation, with a 20-day temporal sampling. The velocities are about 0 using the TI formulation after October 2018 because the velocity observation based on Venus cannot be included inside the system.</figcaption></figure>

<figure><img src="image_4.png" /><figcaption>Fig. 3. RMSE between the synthetic multi-sensor velocities and the estimated velocities over different temporal baselines using the traditional TI formulation (in red), the TICOF formulation (in blue), the TICO formulation (in green) and the moving average (in gold).</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>Fig. 4. RMSE over stable ground for different temporal samplings before (in gold) and after TICOF inversion (in blue).</figcaption></figure>

## B. Spatio-temporal analysis of the results

The spatio-temporal evolution of the velocity magnitude is computed over a longitudinal profile defined in Fig. S1. The results from the TICOF formulation using a 20-days temporal sampling are compared with the results from a 20-days moving average. It can be seen in Fig. 5 that there are no available velocity observations from April 2017 to November 2017 at

the center part of the glacier tongue.

After inversion, there are visually fewer outliers (e.g. there is no more velocity with a zero value at 3.8 km). Moreover, even if TICOF is pixel-wise, the LF time series seems to be spatially coherent: the velocity field is smooth through the period. A seasonal trend can also be observed: the velocities increase between the two periods September to November 2017 and April to June 2018 with a median increase of 40% along the longitudinal profile. From April/June 2018 to September/November 2018, the velocities decrease with a median decrease of 30%. By using all the available velocity observations, TICOF recovers a seasonal trend when TI or TICO fails because only few velocity observations could have been included in TI or TICO in this case (Fig. S5).

## VI. CONCLUSION

The proposed approach enables the fusion of glacier velocity observations from different temporal baselines and sensors without prior information on the glacier dynamics and the data uncertainties. It relies on a new formulation of the Temporal Inversion (TI) based on the temporal closure of the displacement observations network. This new formulation establishes linear equations between combinations of displacement observations and fractions of LF displacements and is named TICOF (Temporal Inversion using COmbinations and Fractions of displacements). The inversion is solved using an Iterative Reweighted Least Square inversion where the initial weights are based on a combination of two statistical criteria: the Median Angle and the modified z-score. The use of displacement fractions assumes a constant velocity over the considered time interval. The bias caused by the validity of this
<figure><img src="image_6.png" /><figcaption>Fig. 5. Spatio-temporal evolution of velocity magnitude along a longitudinal profile defined in Fig. S1. The left-hand figure corresponds to the results of the TICOF inversion. The right-hand figure corresponds to the 20-days moving average, the results on which are compared to the results of the TICOF inversion.</figcaption></figure>

assumption is evaluated on synthetic data. For a sinusoidal signal with a period of 1 year, the bias introduced is below 0.1 m/y for a temporal sampling of 5 to 20 days. This is compared to the improvement in velocity uncertainty obtained after TICOF. For time sampling less than 20 days, the use of TICOF is still beneficial (the uncertainty is reduced by 13%). The 20-day LF velocity time series makes it possible to highlight the strong seasonality of the Fox glacier with lower uncertainty, fewer outliers, and no more gaps. The fluctuation of Fox surface velocities is on average 35 % between spring and autumn from 2017 to 2018. Since this approach does not require prior information, it can be applied to any kind of available velocity datasets computed by any research team.

## REFERENCES

[1] A. S. Gardner, G. Moholdt, T. Scambos, M. Fahnstock, S. Ligtenberg, M. Van Den Broeke, and J. Nilsson, "Increased west antarctic and unchanged east antarctic ice discharge over the last 7 years," *Cryosphere*, vol. 12, no. 2, pp. 521-547, 2018.

[2] M. Fahnstock, T. Scambos, T. Moon, A. Gardner, T. Haran, and M. Klinger, "Rapid large-area mapping of ice flow using landsat 8," *Remote Sensing of Environment*, vol. 185, pp. 84-94, 2016.

[3] P. Friedl, M. Braun, T. Seehaus, S. Lipl, and P. Malz, "The retreat-project: derivation and distribution of remotely sensed variables for glaciers and ice caps outside of the polar ice sheets." in *Geophysical Research Abstracts*, vol. 21, 2019.

[4] R. Millan, J. Mouginot, A. Rabatel, S. Jeong, D. Cusicanqui, A. Derkacheva, and M. Chekki, "Mapping surface flow velocity of glaciers at regional scale using a multiple sensors approach," *Remote Sensing*, vol. 11, no. 21, p. 2498, 2019.

[5] R. Millan, J. Mouginot, A. Rabatel, and M. Morlighem, "Ice velocity and thickness of the world's glaciers," *Nature Geoscience*, vol. 15, no. 2, pp. 124-129, 2022.

[6] L. Charrier, Y. Yan, E. C. Koeniguer, S. Leins, and E. Trouvé, "Extraction of velocity time series with an optimal temporal sampling from displacement observation networks," *IEEE Transactions on Geoscience and Remote Sensing*, vol. 60, pp. 1-10, 2022.

[7] P. Berardino, G. Fornaro, R. Lanari, and E. Sansosti, "A new algorithm for surface deformation monitoring based on small baseline differential sar interferograms," *IEEE Transactions on geoscience and remote sensing*, vol. 40, no. 11, pp. 2375-2383, 2002.

[8] B. Altena, T. A. Scambos, M. Fahnstock, and A. Kääb, "Extracting recent short-term glacier velocity evolution over southern alaska and the yukon from a large collection of landsat data," *The Cryosphere*, vol. 13, no. 3, pp. 795-814, 2019.

[9] P. Lacroix, G. Araujo, J. Hollingsworth, and E. Taipe, "Self-entrainment motion of a slow-moving landslide inferred from landsat-8 time series," *Journal of Geophysical Research: Earth Surface*, vol. 124, no. 5, pp. 1201-1216, 2019.

[10] C. Ding, L. Zhang, M. Liao, G. Feng, J. Dong, M. Ao, and Y. Yu, "Quantifying the spatio-temporal patterns of dune migration near minqin oasis in northwestern china with time series of landsat-8 and sentinel-2 observations," *Remote Sensing of Environment*, vol. 236, p. 111498, 2020.

[11] A. Dille, F. Kervyn, A. L. Handwerger, N. d'Oreye, D. Derauw, T. M. Bibentyo, S. Samsonov, J.-P. Malet, M. Kervyn, and O. Dewitte, "When image correlation is needed: Unravelling the complex dynamics of a slow-moving landslide in the tropics with dense radar and optical time series," *Remote Sensing of Environment*, vol. 258, p. 112402, 2021.

[12] N. Bontemps, P. Lacroix, and M.-P. Doin, "Inversion of deformation fields time-series from optical images, and application to the long term kinematics of slow-moving landslides in peru," *Remote Sensing of Environment*, vol. 210, pp. 144-158, 2018.

[13] A. Pepe, M. Bonano, Q. Zhao, T. Yang, and H. Wang, "The use of c-/x-band time-gapped sar data and geotechnical models for the study of shanghai's ocean-reclaimed lands through the sbas-dinsar technique," *Remote Sensing*, vol. 8, no. 11, p. 911, 2016.

[14] E. W. Burgess, R. R. Forster, and C. F. Larsen, "Flow velocities of alaskan glaciers," *Nature communications*, vol. 4, no. 1, pp. 1-8, 2013.

[15] B. Altena, A. Kääb, and B. Wouters, "Correlation dispersion as a measure to better estimate uncertainty of remotely sensed glacier displacements," *The Cryosphere Discussions*, pp. 1-23, 2021.

[16] H. Purdie, M. Brook, and I. Fuller, "Seasonal variation in ablation and surface velocity on a temperate maritime glacier: Fox glacier, new zealand," *Arctic, Antarctic, and Alpine Research*, vol. 40, no. 1, pp. 140-147, 2008.

[17] C. Lüttig, N. Neckel, and A. Humbert, "A combined approach for filtering ice surface velocity fields derived from remote sensing methods," *Remote Sensing*, vol. 9, no. 10, p. 1062, 2017.

[18] R. A. Maronna, R. D. Martin, V. J. Yohai, and M. Salibián-Barrera, *Robust statistics: theory and methods (with R)*. John Wiley & Sons, 2019.

[19] S. Samsonov, K. Tiampo, and R. Cassotto, "Measuring the state and temporal evolution of glaciers in alaska and yukon using synthetic-aperture-radar-derived (sar-derived) 3d time series of glacier surface flow," *The Cryosphere*, vol. 15, no. 9, pp. 4221-4239, 2021.

[20] C. A. Greene, A. S. Gardner, and L. C. Andrews, "Detecting seasonal ice dynamics in satellite images," *The Cryosphere*, vol. 14, no. 12, pp. 4365-4378, 2020.

[21] P. López-Quiroz, M.-P. Doin, F. Tupin, P. Briole, and J.-M. Nicolas, "Time series analysis of mexico city subsidence constrained by radar interferometry," *Journal of Applied Geophysics*, vol. 69, no. 1, pp. 1-15, 2009.

**Acknowledgment:** We acknowledge CNES (French Spatial Agency) and ISA (Israel Space Agency) for the VENUS observations acquired through the VENUS Scientific Mission, and ESA (European Space Agency) for the Sentinel-2 observations. JM and RM acknowledge support from the CNES.
<figure><img src="image_7.png" /><figcaption>Fig. S1. Fox glacier in the Southern Alps of New Zealand. The red dots are spaced 1 km apart. The background image is an optical image from the Google Earth collection.</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>Fig. S3. Residual of the inversion on synthetic data with the traditional TI formulation (in red) and the TICOF formulation (in blue).</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>Fig. S2. Temporal baselines of the velocity observations according to their central date for Sentinel-2 velocities in blue and Venus velocities in orange. The colormap which ranges from light to dark blue or orange indicates the percentage of non-masked data for each temporal baseline and central date.</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>Fig. S4. RMSE between the synthetic mono-sensor velocities and the estimated velocities over different temporal samplings using the traditional TI formulation (in red), the TICOF formulation (in blue), the TICO formulation (in green) and the moving average (in gold). The vertical grey lines correspond to the temporal samplings which are a multiple of the satellite repeat cycle, here 2 days. For these temporal sampling, TICO is more advantageous than TICOF. For other temporal sampling, TICOF performs better. Note that, for some temporal sampling, it is not possible to include any observations inside the TI system, therefore there is no corresponding value of RMSE.</figcaption></figure>
The system of equations corresponding to the TICOF is:

$$
\begin{bmatrix} 1 & 0 & 0 & \dots & 0 & 0 \\ 0 & 1 & 1 & \dots & 0 & 0 \\ \vdots & \vdots & \vdots & \ddots & \vdots & \vdots \\ 0 & 0 & 0 & \dots & 1 & 0 \\ 0 & 0 & 0 & \dots & 0 & 1 \end{bmatrix} \begin{bmatrix} d_{\tau_1, \tau_1+2\Delta\tau} \\ \hat{d}_{\tau_1, \tau_1+\frac{3}{2}\Delta\tau} \\ \hat{d}_{\tau_1+\frac{3}{2}\Delta\tau, \tau_1+3\Delta\tau} \\ \vdots \\ \hat{d}_{\tau_1+(n-2)\Delta\tau, \tau_1+(n-\frac{1}{2})\Delta\tau} \\ \hat{d}_{\tau_1+(n-1)\Delta\tau, \tau_1+n\Delta\tau} \end{bmatrix} = \begin{bmatrix} 1 & 1 & 0 & \dots & 0 \\ 1 & 1 & 1 & \dots & 0 \\ \vdots & \vdots & \vdots & \ddots & \vdots \\ 0 & 0 & 0 & \dots & 1 & \frac{1}{2} \\ 0 & 0 & 0 & \dots & 0 & 1 \\ \vdots & \vdots & \vdots & \ddots & \vdots & \vdots \\ 0 & 0 & 0 & \dots & 0 & 0 \\ 0 & 0 & 0 & \dots & \frac{\lambda}{\Delta\tau} & -\frac{\lambda}{\Delta\tau} \\ 0 & 0 & 0 & \dots & -\frac{\lambda}{\Delta\tau} & \frac{\lambda}{\Delta\tau} \\ \vdots & \vdots & \vdots & \ddots & \vdots & \vdots \\ 0 & 0 & 0 & \dots & \frac{\lambda}{\Delta\tau} & -\frac{\lambda}{\Delta\tau} \\ 0 & 0 & 0 & \dots & 0 & 0 \end{bmatrix} \begin{bmatrix} d_{\tau_1, \tau_1+\Delta\tau} \\ d_{\tau_1+\Delta\tau, \tau_1+2\Delta\tau} \\ \vdots \\ d_{\tau_1+(n-2)\Delta\tau, \tau_1+(n-1)\Delta\tau} \\ d_{\tau_1+(n-1)\Delta\tau, \tau_1+n\Delta\tau} \end{bmatrix} \quad (S.1)
$$

where $\lambda$ corresponds to the damping factor of the regularization. $\Delta\tau$ is the temporal sampling of the LF time series $[\hat{d}_{\tau_1, \tau_1+\Delta\tau}, \dots, \hat{d}_{\tau_1+(n-1)\Delta\tau, \tau_1+n\Delta\tau}]$. $[d_{\tau_1, \tau_1+2\Delta\tau}, \dots, d_{\tau_1+(n-1)\Delta\tau, \tau_1+n\Delta\tau}]$ stands for the displacement observations.

The system of equations including a model based regularization used to determine simulation parameters in section IV is:

$$
\begin{bmatrix} 1 & 1 & 0 & \dots & 0 \\ \vdots & \vdots & \vdots & \ddots & \vdots \\ 0 & 0 & 0 & \dots & 1 \\ \lambda & 0 & 0 & \dots & 0 \\ \vdots & \vdots & \vdots & \ddots & \vdots \\ 0 & 0 & 0 & \dots & \lambda \end{bmatrix} \begin{bmatrix} 0 \\ -\frac{T}{2\pi}[\sin(\frac{2\pi}{T}\tau_1+\Delta\tau)] - \sin(\frac{2\pi}{T}\tau_1)] \\ -\frac{T}{2\pi}[\cos(\frac{2\pi}{T}\tau_1+\Delta\tau)] - \cos(\frac{2\pi}{T}\tau_1)] \\ \vdots \\ -\frac{T}{2\pi}[\cos(\frac{2\pi}{T}\tau_1+n\Delta\tau)] - \cos(\frac{2\pi}{T}\tau_1+(n-1)\Delta\tau)] \\ 0 \end{bmatrix} \begin{bmatrix} \hat{d}_{\tau_1, \tau_1+\Delta\tau} \\ \hat{d}_{\tau_1+(n-1)\Delta\tau, \tau_1+n\Delta\tau} \\ a \\ b \\ d \end{bmatrix} = \begin{bmatrix} 0 \\ -\frac{T}{2\pi}[\sin(\frac{2\pi}{T}\tau_1+\Delta\tau)] - \sin(\frac{2\pi}{T}\tau_1)] \\ -\frac{T}{2\pi}[\cos(\frac{2\pi}{T}\tau_1+\Delta\tau)] - \cos(\frac{2\pi}{T}\tau_1)] \\ \vdots \\ -\frac{T}{2\pi}[\cos(\frac{2\pi}{T}\tau_1+n\Delta\tau)] - \cos(\frac{2\pi}{T}\tau_1+(n-1)\Delta\tau)] \\ 0 \end{bmatrix} \quad (S.2)
$$

where $T$ is the period of the sinusoidal signal. $a, b, c, d$ are the coefficients of the model.

<figure><img src="image_11.png" /><figcaption>Fig. S5. Spatio-temporal evolution of velocity magnitude along a longitudinal profile defined in Fig. S1. The results are from the TI (left) and TICO (right) inversion, with a temporal sampling of 20 days. Because the velocity observations are multi-sensor, with missing data due to outlier removal, all the velocity observations cannot be included inside the system. Therefore, the estimations are either not possible (from 2.4 to 5.1 km with TI) or mainly constrained by the Tikhonov regularization. This results in very smooth variations, incoherent with the observations.</figcaption></figure>