Hadhri Hela
email: hela.hadhri@univ-smb.fr

Vernier Flavien
email: flavien.vernier@univ-smb.fr

Atto Abdourrahmane
email: abdourrahmane.atto@univ-smb.fr

Trouvé Emmanuel
email: emmanuel.trouve@univ-smb.fr

Time-lapse optical flow regularization for geophysical complex phenomena monitoring

Keywords: Computer Vision, Time series, Temporal regularization, Remote/proximal sensing, Natural Outdoor Environment, Glaciers monitoring, Tracking

In this paper, we introduce a framework for the tracking of geophysical complex phenomena via time-lapse images. It includes the regularization of the derived surface motion maps time series. The proposed processing chain addresses five main challenges: undesired camera movement, missing frames, important photometric changes, weak/repetitive texture and model-free dense spatial transformations. In the proposed framework the motion maps time series are obtained via robust pre-processing steps and optical flow computing. The contribution consists of regularizing the resulting velocity and position time series to minimize a temporal closure error in a subsequent stage. This step serves to alleviate the limitations of existing methods in the context of geophysical monitoring. The temporal closure errors are formulated as linear mappings to inverse using signal priors and two formulations are defined along with illustrative cases. Related methods are discussed and extensive experimentation on simulated datasets is carried out to validate the approach and compare between the different proposed formulations and resolution schemes. Experimental results are presented on time series acquired by groundbased cameras used for the monitoring of Alpine glaciers. The algorithm is computationally efficient, even considering the quantity of processed and generated data, and is run in parallel on multiple cores for speed-up.

Introduction

Motion tracking and scene understanding from time-lapse images acquired by automatic ground-based cameras is a task that has been studied in many applications, including medical imaging tracking systems, video surveillance or geophysical monitoring. With a controlled environment scenario, this motion is characterized by the flow of pixels with constant intensity. This assumption of intensity conservation is typically unrealistic for outdoor geophysical monitoring data because illumination difference between images is relatively important, and more photometric invariance must be introduced. There are two global families of computation techniques to achieve this: (1) Differential techniques that compute velocity from spatio-temporal derivatives of pixel intensities while incorporating other constancy assumptions and regularity constraints [START_REF] Lannutti | Optical Flow Applied to Time-Lapse Image Series to Estimate Glacier Motion in the Southern Patagonia Ice Field[END_REF][2]; (2) Region-based offset-tracking techniques, which compute the trajectory of pixels via the identification of local displacements that provide optimal similarity of two image patches and which usually create robustness for intensity variation via normalization. In this paper, optical flow (OF) methods refer to the first type of methods while the second type is referred to as normalized cross-correlation (but other similarity measures exist) [START_REF] Travelletti | Correlation of multi-temporal ground-based optical images for landslide monitoring: Application, potential and limitations[END_REF] [START_REF] Schubert | Glacier surface velocity estimation using repeat terrasar-x images: Wavelet-vs. correlation-based image matching[END_REF]. Moreover, geophysical monitoring data exhibit other particularities, such as geometric non-rigid deformations, repetitive patterns, low texture, occlusion, frame mis-registration and weather artifacts (such as snow, fog, water drops/sun beam reflection on the screen, clouds, shadows, etc.) [START_REF] Travelletti | Correlation of multi-temporal ground-based optical images for landslide monitoring: Application, potential and limitations[END_REF]. The main challenge in these applications is to automatically estimate accurate long term motion time series from pairwise image registration despite the data contamination. To overcome the limitations of existing techniques in such conditions, we propose a data processing methodology which includes several steps: [START_REF] Lannutti | Optical Flow Applied to Time-Lapse Image Series to Estimate Glacier Motion in the Southern Patagonia Ice Field[END_REF] pre-processing steps: to discard some of the wrong measurements; (2) parameterized motion models: to handle frame mis-registration due to camera motion; and (3) robust optical flow: to estimate apparent motion of the object of interest over time. To introduce further robustness, we propose a regularization step that takes the generated incomplete pixel-wise Velocity and Position Time Series (VTS and PTS) and forces temporal closure, which is the main contribution of the paper. This regularization is inspired by the Small Baseline Subset (SBAS) technique for Synthetic Aperture Radar (SAR) interferometry [START_REF] Berardino | A new algorithm for surface deformation monitoring based on small baseline differential sar interferograms[END_REF]. In this domain, time series of complex SAR images are acquired to observe ground surface deformations [START_REF] Massonnet | Radar interferometry and its application to changes in the earth's surface[END_REF]. Redundant displacement fields are obtained by computing phase differences from image pairs selected according to their temporal and spatial baselines and the observed phenomena: urban and land subsidence [START_REF] Lanari | Sentinel-1 results: Sbas-dinsar processing chain developments and land subsidence analysis[END_REF][8] [START_REF] Yan | Mexico city subsidence measured by insar time series: Joint analysis using ps and sbas approaches[END_REF], water level monitoring [START_REF] Hong | Small temporal baseline subset (stbas): A new insar technique for multi-temporal monitoring wetland's water level changes[END_REF], volcanic deformation [START_REF] Lanari | An overview of the small baseline subset algorithm: A dinsar technique for surface deformation analysis[END_REF], slope movements [START_REF] Barboux | Mapping slope movements in alpine environments using terrasar-x interferometric methods[END_REF] and so on. These initial observations suffer from different error sources such as noise, atmospheric artifacts and phase unwrapping errors. The SBAS approach consists of deriving a final displacement time series with reduced uncertainty by combining these observations in a linear system inverted by different techniques [START_REF] Lauknes | Insar deformation time series using an l 1 -norm small-baseline approach[END_REF] [START_REF] Goel | An advanced algorithm for deformation estimation in non-urban areas[END_REF]. An overview of the parameters of the SBAS algorithm can be found in [START_REF] Lanari | An overview of the small baseline subset algorithm: A dinsar technique for surface deformation analysis[END_REF]. Similarly to the SBAS workflow, the proposed regularization is formulated as an inverse problem where a temporal closure error within an appropriate range is minimized to reduce uncertainty via redundancy.

Related methods

Data regularization

Major time series regularization techniques are linked and include modelling, query by content, anomaly correction, motif discovery, prediction, clustering, classification, dimensionality reduction, indexing and segmentation. A survey of these methods can be found in [START_REF] Esling | Time-series data mining[END_REF]. Hence, this regularization is often done by:

• rolling statistics [START_REF]Rolling Analysis of Time Series[END_REF] • filter like methods such as a non-linear median filtering or extended Kalman filters [START_REF] Walker | Kalman Filtering of Time Series Data[END_REF][18]

• other stochastic approaches such as Multiple Imputation and Bootstrapping [START_REF] Efron | Missing data, imputation, and the bootstrap[END_REF] • model fitting methods; such as Self-Organizing Maps [START_REF] Guilherme | Time Series Prediction with the Self-Organizing Map: A Review[END_REF], parametric/non-parametric regression [START_REF] Gemai | Parametric and non-parametric modelling of time series -an empirical study[END_REF], AutoregRessive Moving Average (ARMA) models [START_REF] Xu | Real-time camera tracking for marker-less and unprepared augmented reality environments[END_REF], etc.,

• supervised learning methods, such as RANdom SAmple Consensus (RANSAC), Support Vector Machine (SVM), etc.

• spectral methods such as principal component analysis/empirical orthogonal function decomposition, Discrete Wavelet Transform, etc. [START_REF] Ghil | Advanced spectral methods for climatic time series[END_REF] Geo-statisticians often use linear-weighted averaging methods, with data driven weighting functions, such as the interpolation/prediction kriging and collocation methods and deploy ground means, such as corner reflectors or GPS measurements to constrain the problem by fusing measurements from different sources. Recent works for the monitoring of geophysical phenomena via time-lapses include:

• adaptive Least Squares collocation to detect outliers in VTS and PTS during feature tracking [START_REF] Mathias | Monitoring Rock Glaciers by Combining Photogrammetric and GNSS-Based Methods[END_REF],

• non-linear median filtering by averaging the values of Nearest Neighbors on surface change time series to correct random noise and fill missing observations [START_REF] Eltner | Time lapse structure-from-motion photogrammetry for continuous geomorphic monitoring[END_REF]. Machine learning techniques have also been used to learn the expected behavior of time series in order to regularize observations [START_REF] Ludmila | An Experimental Study on Rotation Forest Ensembles[END_REF][27] [START_REF] Li | A vector forecasting model for fuzzy time series[END_REF] [START_REF] Schulz | Artificial neural networks for modeling time series of beach litter in the southern north sea[END_REF] and it has been shown that fully Connected Networks are more suited for time series tasks as they allow for the input and output to have the same dimension and to process in a semantic way. A survey on Deep Learning methods for time series inquiry can be found in [START_REF] Cristian | Deep learning for time-series analysis[END_REF]. On the other hand, optical measurements are generally regularized via computer vision derived methods such as detection/recognition methods using interest points, features, HoGs (Histogram of oriented Gradients), etc. For example, the current best performing optical flow methods often use prior knowledge of the properties of the motion, like local smoothness assumptions [START_REF] Berthold | Determining optical flow[END_REF], structure and motion adaptive assumptions [START_REF] Wedel | Structure-and motion-adaptive regularization for high accuracy optic flow[END_REF] or the assumption that motion discontinuities are more likely at image edges [START_REF] Revaud | Epicflow: Edge-preserving interpolation of correspondences for optical flow[END_REF] to ensure spatial coherence over time. In [START_REF] Xu Chen | Optical flow for incompressible turbulence motion estimation[END_REF] the authors propose a method for optical flow adaptation to fluid monitoring using the continuity equation of 2D incompressible flows and passive scalar transport constraints. There are several methods to optimally balance data fidelity with regularization of the motion fields obtained by differential optical flow methods. Depending on the strength of the local image gradients, these methods yield a smooth trade-off between matching and interpolation, which can be very useful in the presence of noisy optical observations [START_REF] Vogel | Optical Flow for Glacier Motion Estimation[END_REF].

Inversion

Inversion consists of inferring the manifold of structural models that could yield some actual observations y through a given transformation A induced by a physical model. Theoretically, the goal of a reconstructive inversion is often to estimate a 2D or 3D field x ∈ M: the Model space, from a set of observations or measurements y ∈ D: the Data space, of the form y = Ax + b. The measurement or model operator A can be linear or non-linear, full rank or rank-deficient and has a conditioning that can be quantified by means of the L = 2 condition number K defined as the ratio between the biggest and smallest singular values σ of A [START_REF] Cucker | Chapter 1: Normwise Condition of Linear Equation Solving & Chapter 5: Condition Numbers and Iterative Algorithms[END_REF], such as:

K = σ max σ min ,
which quantifies the sensitivity of the estimation to perturbations in the observations. b is the noise or measure uncertainty which hopefully will not be reconstructed by inversion. Inverting observations is now used extensively in many fields like remote sensing [START_REF] Zhang | Interferometric phase reconstruction using simplified coherence network[END_REF], deconvolution, optical tomography, 3D reconstruction [37][38] or geodesy [START_REF] Tarantola | A strategy for nonlinear elastic inversion of seismic reflection data[END_REF]. A typical photogrammetric example would be recovering the descriptions of objects represented in image data x via optical measurements y corrupted by noise b by inverting the optical model A (back-projection). For this, a model misfit cost between y and Ax minimized. The used cost function is usually convex (unless a special penalizer is needed to create specific robustness) if the direct linear mapping is full-rank. Depending on the definition, sparseness and dimension of the direct problem, on the prior information on the estimated solution and on the choice of the cost function, direct (spectral decomposition) or iterative (optimization via a solver) methods can be considered to stably recover x. In this paper, the physical model, that is encoded in A, and to be inverted, is the consistency of the displacement pattern over time: the temporal closure constraint. The inversion is thus done on the temporal axis of the computed OF maps and densely for every coordinate in order to obtain temporally regularized Optical Flow (OF) maps time series of the size of the images. The inversion is pixelwise and two formulations are proposed to encode the linear temporal closure model A.

Contribution

Despite the existence of discontinuity-preserving smoothness constraints to ensure spatial coherence of optical OF fields, one can notice the lack of long-term temporal constraints in recent models proposed for the Middlebury benchmark [START_REF] Baker | A database and evaluation methodology for optical flow[END_REF] due to computational requirements. Only a few methods make use of more than two frames and none of them incorporate temporal information in the underlying model. For off-line remote sensing of geophysical complex phenomena, we propose to revisit the Differential Synthetic Aperture Radar interferometry SBAS technique to minimize a temporal error from a Network of Optical Flow Fields (NOFF) forming the observations. The workflow to pre-process the raw data and compute relevant pairwise optical flow fields in object coordinates is presented followed by the paper's main contribution: the regularization of the pixelwise Velocity Time Series (VTS) and/or Position Time Series (PTS). To validate the whole pipeline under various conditions, 4 realistic datasets of simulated time lapses are constructed. To validate the temporal regularization step, different processing problems such as outliers, missing data and noise are further simulated on the observations and estimations are compared to ground truth solutions. Finally, results for the monitoring of alpine glaciers are presented.

Proposed pipeline

Overview

In a context such as stated before, the computed observations can be contaminated in different ways: bias can be introduced by the motion of the fixed camera; outliers by the photometric variance distorting the calculations in one of the layers of the automatic processing chain and/or the registration failing to capture the non-rigid spatial transformations; and gaps can be caused by targets being occluded and/or outliers due to weather artifacts being detected and removed. Figure 1 shows the layers of the proposed pipeline to process time-lapses in such a context: (A)Data is first acquired automatically via a programmable micro-controller; (B) a texture score is computed to dismiss images where the information is too corrupted and the remaining images are co-registered with regards to a master image to compensate for camera motion.

Data selection

After the acquisition of the time-lapse imagery, images are selected given their meta-data tags and according to the desired time step of the resulting time series. Textural information is quantified to dismiss corrupted images. This information is retrieved by the construction of a bank of multi-parameters Gabor filters designed manually to recognize the known textures (fringes, grains, tales, stripes, ice foliation, etc.) of the object of interest, if not occluded or highly deformed. This bank is constructed by varying some parameters (frequency, orientation and size) of the 2D Gabor function to create invariance to small differences between the images to be selected for further processing. A similar work can be found in [START_REF] Li | Selection of gabor filters for improved texture feature extraction[END_REF]. For every image i, the sum of all the energies E of the obtained response maps is quantified. The energies of all the images in the time-lapse should theoretically follow a Gaussian distribution and the spurious images are detected by thresholding over median absolute deviation (MAD):

Ẽ = median(E i ), (1)
S 1i = median(|E i -Ẽ|) (2)
S 1 is thus the texture score encoding the object of interest visibility in image I i . This way, and as the known orientation of the gradients is also quantified by the designed filters, images that have a lack of textural information (foggy, occluded, damaged, cut off, etc.) are rejected and are no longer part of the pipeline. This classification gave a success rate of 96% for the rejection of data affected by adverse weather conditions and mis-acquisition.

Pairwise image processing

Co-registration

The stack of remaining images are co-registered by estimating a parametric geometric transformation regarding a chosen master image, the first relevant one for example. This step is done to compensate for the undesired geometric differences of the images in the time-lapse (see figure 1 (A)) and is thus done only on the static zone of the observed scene. This geometric difference is due to camera orientation and image plane instability due to extreme weather conditions. As the observed scene is distant, the images can be considered as the same planar surface in space, and a homography H can encode the translation, rotation, the scale and also the tilt (8 degrees of freedom) of the camera between two images I i and I j with 8 parameters param (2 scaling, 2 rotation, 2 translation and 2 vanishing points parameters). The arguments param defining the homography Φ are encoded in non-singular 3×3 matrix computed in terms of the Enhanced Correlation Coefficient criterion [START_REF] Evangelidis | Parametric Image Alignment Using Enhanced Correlation Coefficient Maximization[END_REF] on the static zones by iteratively minimizing the following normalized Euclidean distance:

E ECC (i,j) = ( I i (X) ∥I i (X)∥ - I j (Φ(X; param)) ∥I j (Φ(X; param))∥ ) 2 , (3)
where I i (x) is the master image and I j (Φ(X; param)) is the image warped by inverse composition of pixels of coordinate X (= (x i , y i , 1) T (in homogeneous coordinates notation) with Φ:

          x j y j 1           = Φ(X; param) = H -1 X = H -1           x i y i 1           (4)

Matching and optical flow

Deep-Matching. For every pair of images in between which an Optical Flow (OF) field is computed, an initial matching is first quantified using a semi-dense matching algorithm called DeepMatching (DM) [START_REF] Revaud | DeepMatching: Hierarchical Deformable Dense Matching[END_REF]. In DM, the correlation is quantified by means of Normalized Cross Correlation (NCC) of HOG descriptors of 4×4 pixel patches called atomic patches. Furthermore, this correlation is computed in a pyramidal scheme to, first, aggregate coarser correspondence by going up in the pyramid, and then, refine the correlation correspondences by optimizing the backtracking to the bottom of the pyramid and retrieve optimized matches of atomic patches.

Besides the scale invariance insured by this multi-scale pyramidal scheme, the used operators for recursive construction of the correlation pyramid make further analysis weakly affected by non-rigidness and, more importantly, shows robustness against repetitive patterns. For speedup, the parameters are tuned so the search area, the rotation and scale factors of a match are restricted to a range that can not be exceeded.

Optical flow computing. These dense matches are then used to initialize a further computation of a dense OF field ⃗ F that minimizes a weighted energy functional based on two data terms E DB and E DG , from the brightness constancy assumption and the gradient constancy assumption respectively and a smoothness term E REG . These assumptions respectively stipulate that the color and the edges between two successive frames are conserved, and the regularizer encourages the spatial smoothness of the flow field. For this, DeepFlow is used [START_REF] Weinzaepfel | DeepFlow: Large displacement optical flow with deep matching[END_REF](DF). DF is similar to the well known Brox and Malik [START_REF] Brox | Large displacement optical flow: Descriptor matching in variational motion estimation[END_REF] algorithm but, in addition, takes a matching term E M to further constrain the minimization, such as:

E OF (i,j) = ∫ ∫ Ω (E DG (i, j) + αE DB i, j + σE reg (i, j) + βE M (i, j) ) dX, (5)
Furthermore, DF have a weighting strategy that has also been shown to be robust to repetitive texture. A robust cost function ψ(s 2 ) = √ S ² + ϵ 2 with ϵ = 0.001 is used for optimization and for all terms. This makes the energy functional to minimize non-convex and non-linear. Let I i and I j be 2 consecutive frames defined on Ω with 3 RGB channels: I i , I j : Ω → R 3 . The estimated optical flow field is then ⃗ F = (F x , F y ) : Ω → R 2 and can be viewed as the change of structured light between the images, or in other terms, as the sub-pixelic displacement from an image I i to an image I j of the same scene at different time indexes t i and t j . This can be seen as the velocity of the observed scene projected onto a subspace orthogonal to the line-of-sight (image plane) if the set of dates are separated by unit intervals ∆t = t i+1t i = 1. OF fields are thus similar, in the physical meaning, to the phase observed in DinSAR images (which is projected onto the line of sight). The weights α, σ and β are tuned so to favor the gradient constancy assumption over the brightness constancy assumption. This makes the optical flow computation sturdier to photometric differences.

Scores and errors

To assess the pairwise OF maps computation steps, the following scores and error measures are computed:

• The mean texture score (from section 3.2): S 1 .

• A correlation score per couple of images equal to the sum of all patch similarities along their DM back-tracking path: S 2 .

• The inverse of the mean residual velocity amplitude 1 ∑ S Z |F (i, j) | observed on the static zones: R 1 .

Scaling and geo-referencing

For the result to be more reliable and to avoid parallax confusion on the computed Optical Flow maps ⃗ F, scaling, from sensor coordinates into object space metric coordinates, and georeferencing is needed. For this purpose, and considering the pinhole camera geometric and optical equations, a ground sample distance map ∆mm (in cm/pixel) can be computed from intrinsic camera parameters and the depth map Z with respect to the image against which the time-lapse has been co-registered:

∆mm = x mm f Z, (6)
where x mm is the physical metric distance between pixels inside the hardware, or pixel pitch, and f is the lens focal length. This depth information Z was obtained via a shape from motion framework, whose seminal work can be found in [START_REF] Tomasi | Shape and motion from image streams under orthography: a factorization method[END_REF], of stereo image pairs. This 3D photogrammetric workflow is stabilized by several in-situ control points that were measured on the static zones seen by the cameras. These in-situ measurements allowed us to constrain the search for the orientation of the cameras (extrinsic), to perform auto-calibration (intrinsic) and to georeference the obtained 3D model of the object of interest during the bundle-adjustment step.

The depth maps from the camera centers are then retrieved by ray intersection from camera centers to the 3D points. This way, the components of ⃗ F could be scaled to quasi-3D with known GPS location. Unfortunately, there are two main disadvantages to this technique. First, the flow of the scene that is projected onto the line of sight is neglected and, hence, the main direction of the sought after flow should be parallel to the image plane for an accurate augmentation to 3D. Secondly, classic calibration algorithms (such as with calibration chessboards / corners) can introduce extra uncertainty in the case of geophysical monitoring due to ground constraints, important distances and the lack of geometrically stable features. Therefore, the in-situ measurements must be well scattered over the entire scene and measured with a high accuracy and precision for both the GPS and the corresponding pixel coordinates.

Regularization 3.4.1. Temporal closure errors

The cost functions for regularizing the obtained velocity and position time series (respectively VTS and PTS) are the norm of the temporal closure errors denoted b LF and b CM for the Leap-Frog and Common-Master formulations respectively. These temporal closure errors can be expressed as follows: Let N be the number of images and a the master image for the regularization step. Its first expression, the Leap-Frog formulation LF, is given by the forward and backward forms such that: ∀ a ∈ {1..N} ,∀ ranges n ∈ {1..N -a} and ∀µ ∈ {1..a -1} :

-F a,a+n + a+n-1 ∑ i=a F i,i+1 = b LF , F a,a-µ + 1 ∑ i=a-1 -F i,i+1 = b LF . (7)
The second expression is the Common-Master formulation CM and is given by the form such that:

∀ a ∈ {1..N} , ∀i ∈ {a..N} , ∀n ∈ {1..N -i} :

F i,i+n + F a,i -F a,i+n = b CM . (8)
F (i, j) denotes either component F x or F y of the OF field between two images I i and I j with date indexes i and j and at a given pixel p = (x, y) T . These closed loop error measures should be equal to zero if the flow computation of the involved fields is flawless. In practice, temporal closure can be violated not only by the uncertainty introduced over the different layers of the pipeline but an abrupt event can also take place in the observed area inducing the flow computation to fail capturing the non-smooth large motion.

Temporal closure mappings

Given a master image a, to consider errors in Equations ( 7) and ( 8) for all temporal loop sizes n, we express them as the following pixelwise system of equations:

A LF × F est LF -F obs = b LF , A CM × F est CM -F obs = b CM . (9)
The processed observation vector ⃗ F obs belongs to the Data space D = R (N-m)(N-m-1) and is the same for the CM and LF formulations. It consists of observed OF fields between all possible date index arrangements after the rejection of m unreliable images, or the NOFF values at a given pixel (all bi-directed orange arrows in the graphs in figure 2). The sought after estimates ⃗ F est LF and ⃗ F est CM should minimize a misfit error to the model spaces M LF , M CM ⊂ R N-1 respectively (black simple directed contour arrows of the graphs in figure 2). The model space M LF is the set of temporally coherent time series: ⃗

F est LF = {F i,i+1 , i ∈ (0, 1, ..., N -1)} and M CM is the set of temporally coherent time series: ⃗ F est CM = {F a,i , i ∈ (a + 1, ..., a + N)} relative to a Common Master image I a .
The output of inverting the Common Master (CM) formulation is thus, theoretically, the Riemann sum over the output estimated via the Leap-Frog (LF) inversion starting from the master index a such as:

⃗ F est CM (t k ) = k-1 ∑ i=a ⃗ F est LF (t i ). (10)
Both temporal closure error formulations in Equations ( 7) and ( 8) can thus be deducted from one another and whereas ⃗ F est LF can be interpreted as the motion velocity between adjacent date indexes (at a given pixel), ⃗ F est CM can be seen as the motion position relative to the Common Master date a. The output of the LF inversion is thus a VTS and the output of inverting the CM formulation is a PTS. Illustration case. To illustrate the temporal closure mappings to inverse, let us consider a timelapse of 18 images. Figure 2a shows the Network of Optical Flow Fields (NOFF) index graph, i.e. the graph of indexes in between which optical flow fields are computed (in orange) and the index of the estimations (in black). Figure 3 show the induced linear mappings and their corresponding Singular Values, respectively for the LF and CM formulations. We can see that both mappings are full-rank and have a low condition number (12.06 and 4.35 respectively for the LF and CM mappings) which suggests that the solution of the inversion process won't be too affected by noise in the observations. Reducing the maximum range. Considering that the registration error increases with the temporal baseline and the important time span for geophysical monitoring, reducing the maximum range of the temporal closure loops allows to reduce mis-registration related errors. Figure 2b shows the obtained index graph and figure 4 the corresponding linear mappings. We can see that The rows and columns are arranged as in equations ( 11) and ( 12) but only 1 on 20 homologous indexes are displayed on the vertical axis.

this operation does not significantly affect the condition numbers of the mappings to inverse (5 and 7, 79 respectively for the Leap-Frog and Common Master mappings).

We can also see that the CM mapping presents smaller singular values for both the full and reduced range cases. As a practical matter, small singular values are usually a sign of ill-conditioning, in the one hand, and it can be difficult to distinguish zero and relatively small Singular Values in a large scale problem, on the other hand. For these reasons, we tend to say that, theoretically, inverting the LF formulation is less sensitive to errors in the observations. Missing observations. Let us now suppose that images with date indexes 2, 3, 14, 16 and 17 have been discarded from pre-processing for lack of texture or high deformation 3.2. The index graph in figure 2b becomes as in figure 5 and the mappings become rank deficient (K is infinite).

The obtained cost to minimize, respectively for the errors in Equations ( 7) and ( 8) and at a given pixel p are:

∥b LF ∥ L = || ⃗ F obs -A LF ⃗ F est || L =                               F obs 0,1 F obs 0,2 . . F obs 0,5 F obs 1,0 F obs 1,2 . . F obs 18,13 F obs 18,15                               -                      
. . -1 -1 -1 -1 -1 0 0 0 . . . 0 0 -1 -1 -1                                          F est 0,1 F est 1,2 F est 2,3 . . . F est 16,17 F est 17,18                    L , ( 11
)
∥b CM ∥ L = || ⃗ F obs -A CM ⃗ F est || L =                              F obs 0,1 F obs 0,2 . . F obs 0,5 F obs 1,0 F obs 1,2 . . F obs 18,13 F obs 18,15                              -                       1 0 0 . . . 0 0 0 0 0 0 1 0 . . . 0 0 0 0 0 . . 0 0 0 . . . 0 0 0 0 0 -1 0 0 . . . 0 0 0 0 0 -1 1 0 . . . 0 0 0 0 0 . . 0 0 0 . . . 1 0 0 0 -1 0 0 0 . . . 0 0 1 0 -1                                              F est 0,1 F est 0,2 F est 0,3 . . . F est 0,16 F est 0,17 F est 0,18                        L . ( 12
)
In this example, the Common Master image is the first valid one, i.e. I 0 and both linear mappings have rank 13. The dimension of the Nullity space of A equals the number of rejected image indexes leading to 5 linearly dependent columns out of 18. The condition-number is effectively infinite, the data space D has shrunk and further regularization or Model space M conditioning is needed to lead to a satisfying solution after inversion. To test the behaviour of inverting both systems on the missing data indexes, a random observation vector was generated by assuming a gaussian process scaled by the temporal step. Figure 6 shows the simulated outputs. The values corresponding to missing observations are filled with a constant in F est LF (VTS) and with zeros in F est CM (PTS). To further illustrate how the solutions might be affected by missing observations, figure 7 illustrates the mappings when facing an extreme rejection case (that can occur in situations where the acquisition of the time-lapse is subject to extreme external conditions) and the obtained solutions ⃗ F est with manually simulated data that respects the temporal closure constraint but lacks several data.

Analytical aspect and optimization

As stated before, the mappings from the model spaces M LF and M CM to the Data space D can be full rank, i.e. rank(A) = N, or rank deficient, depending on the number of missing images. Given that the rows of A are by definition redundant constraints, the row rank can never be full and hence the temporal closure linear mappings are never injective. This means that not every observation ⃗ F obs ∈ D has a temporally coherent solution. Apart from specific situations where most of the data is rejected from preprocessing (see illustration case of figure 7), the problem is over-determined. Generally, for overdetermined problems, an optimal solution to the system of equations is sought by resolving the normal equation: A T A ⃗ F est = A T ⃗ F obs either directly or by SVD decomposition. This is equivalent to solving the problem in Equation [START_REF] Lauknes | Insar deformation time series using an l 1 -norm small-baseline approach[END_REF] where the cost function is the squared L = 2 norm Ψ : x → ||y -Ax|| 2 = ||y|| 2 -2y T Ax + ||Ax|| 2 . The linear mapping Ψ is twice differentiable and its Hessian is the positive semi-definite matrix 2A T A. The Least Squares minimization is thus convex and should lead to unique optimal VTS and PTS solutions, respectively, ⃗ F est LF and ⃗ F est CM .

arg min

⃗ F est ||A × ⃗ F est -⃗ F obs || L , ⃗ F obs ∈ D (13)
However, resolving the normal equation directly breaks down for the under-determined case (severe rejection cases) and the rank deficient case (missing observations). A truncated SVD decomposition is the most commonly used method for this scenario, and theoretically, gives the best unbiased Least Squares estimation. Nonetheless, even using a truncated SVD may lead to numerical difficulties as, for sparse and rank deficient problems, the problem is not strictly convex (presence of local minima) and this can lead to divergent solution when the observations are contaminated by noise and especially outliers. Most approaches for avoiding these direct decomposition methods are based on dedicated incremental decomposition methods or on constrained/unconstrained linear and quadratic programming. Iteratively learning the inverse has the advantage of allowing to add equality/inequality constraints, to easily change the cost to be minimized to suit the prior knowledge on the solution and the errors, to add an initialization for the solution, and to optimize memory consumption by taking full advantage of the sparsity of the direct mappings. For minimizing a functional this way, one can extensively explore the model space, systematically (brute force), randomly (Monte Carlo), or use local properties of the cost function space to find a global minimum. There are many optimization algorithms and missing observations and rank-deficiency is shown to be handled differently according to the inner local properties of the used solver. The solvers that are relevant in this context were thus used in order to make the method tractable and will be cited and tested in Section 4.2 with simulated condition acquisition. More about direct and iterative methods for linear sparse problems can be found in [START_REF] Åke Bjurck | Methods for sparse linear least squares problems[END_REF][48] [START_REF] Scales | Robust methods in inverse theory[END_REF].

Missing data and conditioning. As stated above, if there are no missing frames, the problem is full (column) rank. The convex resolution process theoretically leads to similar minimas with direct and iterative methods, after velocity/position conversion and for both LF and CM formulations. In case an image I i has been rejected from preprocessing because of lack of information in the area of interest, the rows of A and ⃗ F obs involving the index i become obsolete and are deleted. For m deleted images over N, we have 2Nmm(m + 1) constraints less. The decay of the rank depends on m and on the contingency of the missing observations. In this case, for r = rank(A) parameter to estimate the problem is still overdetermined (too many constraints involving the same unknowns) and for the Nr values involving a rejected index, the problem is under-determined (redundant columns for LF system and zeroed columns for the CM system).

Reducing the maximum range and conditioning. It is important to note here that reducing the maximum range (as in figure 4) of the temporal closure loops can also have analytical benefits for the LF formulation. In fact, this operation prevents the matrix A LF from having a dense row and, indeed, if a sparse matrix A have even one dense row, A T A becomes full and iterative and direct methods using the normal equation are affected (much slower).

Error propagation simulation. Observational uncertainties (random errors) are assumed to be normally distributed with zero mean and with prior covariance Σ obs . By "prior" we imply that they are not related to the inversion step but to the data and the prior layers of the pipeline. The uncertainties related to the observations propagate through the inversion process, leading to estimations ⃗ F est LF and ⃗ F est CM with covariances Σ v LF (uncertainties of the LF estimated VTS) and Σ p CM (uncertainties of the CM estimated PTS). The covariance of the estimations can be approximated by the rule of linear functions of random variables as in equation:

Σ est = A † Σ obs A † T , (14)
where † stands for the pseudo-inverse of A. As the mappings from VTS to PTS (and vice-versa) are linear, one can deduct Σ p LF (uncertainties of the LF estimated PTS after conversion) from Σ v LF and Σ v CM (uncertainties of the CM estimated VTS after conversion) from Σ p CM . To ilustrate, let us consider statistically independent observational errors with variance 1, in which case Σ obs = σ 2 obs I. figure 8 shows the obtained covariance matrices of the estimates. The variances can be used to state data confidence bounds for the model parameters, e.g. F est (i,j) ± 2σ F est (i,j) (but can be misleading in the case where the estimations are highly correlated). Figure 8 suggests that VTS are less affected by uncertainty accumulation than PTS, on the one hand, and that the LF minimization is less prone to observational uncertainty propagation on the missing data indexes, on the other hand. Weighting. In order to maximize the efficiency of the regularization and tackle heteroscedasticity, which is the fact that the observations are measured with different uncertainties, one can take into account the proper amount of influence over the observations by using the various residual error estimations and scores accumulated over the pipeline (Section 3.3). This is equivalent to preconditioning the matrix A or to minimizing the weighted norm ∥.∥ W of the temporal closure error b where W represents the weights of the inner-product defined by < a, b > W = b * Wa = b T Wa and thus must be a positive-definite bilinear preferably symmetric form. Practically, this can be implemented as follows: arg min

⃗ F est ||W(A ⃗ F est -⃗ Fobs)|| L . (15)
The weighting matrix W is diagonal and the weights have to be estimated up to a proportionality factor. To focus accuracy, every observation F obs (i, j) is given a weight W (i, j) corresponding to how well-registered the images I i and I j are (scores S 1 and S 2 in Section 3.3.3). To discount imprecision of the observations one would want to use an estimator that captures the information about how volatile F obs i, j is (score R 1 ). In the Small Baseline Subset workflow, weighting is often done by multiplying by the precision matrix Σ -1 obs .

Adding priors about the solutions. As the sought VTS and PTS are known to be smooth (quasi-stationary time series), a Tikhonov regularization term can also be incorporated in the cost function to avoid landing on "almost minimizers" that have a big norm. This term explicitly incorporates the regularity requirement (minimum variance) in the formulation of the problem.

The damped system is such as:

arg min ⃗ F est ∥W(A LF × ⃗ F est -⃗ F obs )∥ L + λ 2 ∥ ⃗ F est ∥ 2 ⃗ F obs ∈ D (16)
In Equation ( 16), the first term measures the goodness of fit to the temporal closure constraints (in terms of the least L-norm cost), i.e. how well the estimations ⃗ F est predict the noisy observations ⃗ F obs and having important partial residuals ∥A LF w × ⃗ F est -⃗ F obs w ∥ L means that the solution does not solve the temporal closure problem. On the other hand, we should not make the residual smaller than the average size of the errors in ⃗ F obs as we do not want to fit the noise in the data. The second term λ 2 ∥ ⃗ F obs w ∥ 2 measures the regularity of the solution based on the knowledge that the naive estimations (without damping) can be dominated by high-frequency components with large amplitudes we want to avoid considering as solutions.

These assumptions and empirical experiments will determine the choice of λ that will best balance the problem, which is beyond the scope of this paper. The final resolution framework (choice of the solver for resolution and of the norm L) will be discussed in the result section (Section 4.1). The "regularization" term in this case (Tikhonov) is linear (L = 2 norm), but could have been any hand designed regularization priors, such as sparseness in the wavelet space and so on. Recent examples of cost function analysis and the Tikhonov regularization term for signal restoration by inversion can be found in [START_REF] Budd | Regularization techniques for ill-posed inverse problems in data assimilation[END_REF] and [START_REF] Mazzieri Gisela | Regularization of inverse ill-posed problems with l2-bv penalizers and applications to signal restoration[END_REF].

Results on simulations

Simulated datasets

In order to assess the proposed pipeline in different controlled conditions, synthetic timelapses are created by resampling a template image with different types of displacement fields. This template image consists of a sub-image retrieved from real data. A simulated frameworks that consists of 4 synthetically generated time-lapses of a textured semi-rigid fluid (a glacier) with different parameters is obtained. These datasets have a full resolution (1100*1500 pixels covering a small sub-area of a glacier) and a simulated temporal range of N = 6. The different simulated displacement behaviors are:

• Spatial rigidity: No spatial deformation (100% rigid) / Affine geometric transform GA (non-rigid, global)/ Piecewise-Affine geometric transform GPA (non-rigid, local).

• Temporal evolution of the PTS: Linear LT (VTS = constant) / Polynomial PT (non-linear PTS). The template image I 0 (figure 9a) is retrieved from real data and is transformed to simulate images with indexes i = (1, 2, 3, 4, 5). For dataset 1, the transformation is a sub-pixelic rigid motion (global translation) and the temporal evolution of the translation coefficient is linear (LT). For datasets 2 and 3, the transformation is a global affine transform (figure 9b) with greater rotation for dataset 3: for every image i in the simulated time-lapses 2 and 3, a rotation angle r i (counter-clockwise direction), a shear factor sh i and translation factors t x and t y are defined with respect to image i = 0, which will set the six degrees of freedom of the affine transform. The temporal evolution of the translation coefficient is linear for dataset 2 (LT) and non-linear for dataset 3 (PT). The shear parameter of the affine transform increases along the vertical and temporal axis. This allows the analysis of the consequences of increasing non-rigid deformation on Normalized-Cross-Correlation (NCC) and Optical Flow (OF) computing as well as the ability of the regularization step to decrease the error on the computed VTS and PTS. These results are reported on figure 10 where the error related to the shear is approximated by summming over the temporal and horizontal axis, such as:

∑ 1450 x=0 ( ∑ 4 t=0 (ξ t ))
. ξ being the distance to the ground truth. Dataset 4 is undergoing a piecewise-affine transform (figure 9c). Whereas affine transform exhibits several geometric properties that describe a non-rigid global motion as seen on the image plane, with a piecewise-affine transform, an image is modeled as a set of piecewise-convex areas (triangles), each of which is undergoing a different affine displacement. This geometric property grossly models the motion of a semi-rigid fluid subject to local compressions and dilatations. As for most remote sensing applications, only the object of interest is subject to motion, nearby zones remaining static. These static zones are retrieved from a real remote sensing dataset at different times to assess the invariance of the framework to environmental and photometric changes on motionless subimages.

Legend and evaluation strategy.

Several experiments are conducted to first evaluate the whole pipeline compared to a NCC approach with relatively clean data (Experiments A) and, later, to test the robustness of the inversion step to various acquisition scenarios and with the different proposed formulations (Experiments B, C and D):

• Experiments A: normal acquisition scenario

• Experiments B: noisy acquisition scenario

• Experiments C: missing observations

• Experiments D: systematic errors Some other experiments and results will also be discussed.

The presented results are evaluated in terms of the temporal Root Mean Squared distance to the ground truth solution, denoted by the tilde, ξ = || F -F est || 2 averaged over a 500 × 1000 pixels central patch. For the LF formulation ξ refers to the error in the estimated Velocity Time Series (ξ VTS ) and for the CM formulation ξ refers to the error in the estimated Position Time Series (ξ PTS ). For every series of experiments, the iterative method is a weighted non-constrained damped least-squares minimization (L = 2) solved with the LSQR (Sparse Equations and Least Squares) algorithm [48][52]. This method proved to be the most appropriate as it is adapted to Sparse Linear systems (fewer vectors of working storage), to ill-conditioning and to both over and under determined problems. It is a conjugate-gradient type method that iteratively uses the Golub-Kahan bidiagonalization process to solve the normal equation. For the series of experiments presented in Subsections 4.5 and 4.6, other minimization routines that have shown an interesting behavior are presented (respectively; L = 1 norm for the cost function to deal with outliers and solvers that have shown a different strategy to deal with missing data in the observations).

Experiments A: normal acquisition scenario

In this series of experiments, the whole pipeline is evaluated and compared to Normalized Cross Correlation (NCC). Results with both the Leap-Frog and Common Master formulations and with direct (truncated SVD) and iterative resolution are presented before and after regularization and with different weights (Non-weighted (NW), weighted with scores: S = S 1 × S 2 (WS), R (WR) and mixed scores S × R (WSR) (see Section 3.3) Conclusions A Results are reported in tables 2, 3 and 4. Unlike Optical Flow (OF), normalized Cross-Correlation (NCC), fails to capture non-rigid motion, especially on Dataset 3, which is characterized by a more distorting spatial transformation, and presents many outliers on Dataset 4, which is undergoing piecewise non-rigid spatial transformations (Table 4). Non-rigid deformation of the the processed images exhibited indeed abrupt changes in the correlation similarity measure. The error decreases by regularization for all datasets and weighting the minimization provides the solutions that give the smallest ξ VTS and ξ PTS . The mixed weighting with WSR gives better results most of the time (Tables 2, 3 and 4) and is therefore the one used for the iterative method and the rest of the experiments. The iterative damped Least-Squares resolution (LSQR) gives slightly better results than the Truncated SVD method and as the solver converges within 2 iterations, LSQR is a better choice for the task. The experiments with this acquisition scenario for minimizing the norm L = 1 of the temporal closure cost function shown results equivalent to those with L = 2 but with a slower convergence rate. The inversion being pixelwise, the spatial transformation does not theoretically affect the amount of extra accuracy given by the temporal regularization, but we can see from figure 10 that when several observations are affected by registration errors (due to non-rigidity), inversion may fail to regularize the observations.

Experiments B: noisy acquisition scenario

In this series of experiments, an additive random white noise ω is injected spatially to all observations to simulate a measure uncertainty such as: F obs (i, j) ← F obs (i, j) + ω (i, j) . The ω (i, j) 's are independent, identically distributed and drawn from zero-mean normal distributions with different standard deviations σ ω . Conclusions B Results are reported in tables 5 and 6. The regularization process accurately denoises the obtained velocity fields from the added white noise and we can see in figure 11 that static zones have been cleared from residual displacement. Numerical instability associated with the SVD resolution can be seen in red in Table 6. The damped least-squares solution given by LSQR gives a "regularized" solution.

WLF Dataset σ ω = 0.1 pixel σ ω = 1 pixel σ ω =

Experiments C: missing observations

In this series of experiments, images have been discarded from the input of the pipeline to imitate image rejection because of lack of texture (see Section 3.3). 1, 2 and 2 adjacent images are discarded over 6 to evaluate how the regularization process fills the gaps in the estimated VTS and PTS due to missing observations. The gap filling property is specific to the local approximations made by the used solver on the under-determined sub-systems induced by missing observations (local non-convexity). When solving with SVD or LSQR, gaps in the estimated VTS (F est LF ) are filled with a constant value interpolated from neighboring data, and with zeros in the estimated PTS (F est CM ) (see illustration case in figure 6). Other solvers were found to have other local properties for gap-filling. The COBYLA (Constrained Optimization BY Linear Approximation) algorithm [START_REF] Powell | A view of algorithms for optimization without derivatives[END_REF][54][52] shows results similar to LSQR where no-data is missing and a gap filling strategy that gives the closest solutions to ground truth F, at the cost of relatively slower execution time (approximately doubled). The iterative method in this subsection refers thus to the solution estimated by minimizing the damped least squares cost of the temporal closure errors with the COBYLA solver. Conclusions C Results are reported in tables 7 and 8. Empirical results show that the COBYLA method outperforms all the other solvers on all datasets and pixel samples through its local approximations. When data is missing, the Leap-Frog formulation of the temporal closure constraint gives better results than the Common-Master formulation (Figure 12) and the solutions from iterative methods are closer, in terms of the RMSE distance, to the ground truth velocity F.

Experiments D: Systematic errors

In this series of experiments, noise have been injected into (a part of) the processed velocity fields that include a common image I i . The aim is to simulate bias induced by data, or systematic observational error, in the involved OF maps. This noise is drawn from a normal distribution with standard deviation σ = 1 and a mean value µ 0 , which is the parameter we will vary in the experiments. In real data, this systematic observational error can be caused by a false detected event in the image I i . Conclusions D Results are reported in tables 9 and 10. The proposed weights rely strongly on the residual error in the static zones of the observed scene. In case the bias were injected in the dynamic zones (Figure 13), the weights might not "detect" the biased map and, as expected, using a cost function with L = 1 for the inversion is more robust against outliers (least absolute deviation). Using the extra regularization term (cost function in equation ( 16)) prevents the solution from diverging (Table 10).

Other experiments and conclusions

To summarize the results on the simulated datasets, the optimal resolution scheme depends on the nature of the suspected anomalies in the processed velocity and position time series. In the presence of outliers, such as when a high tolerance in the image rejection phase induced the presence of occluded zones in the scene, a norm L = 1 should be used to initialize a further regularization with a Least Squares cost function. Using the damped cost function in equation ( 16) has allowed us to avoid divergent solutions, which can happen in the presence of important noise, outliers or missing data. Extensive experiments showed that when combining two acquisition scenarios, the conclusions for both scenarios remain. Initializing with observations makes the convergence faster, in terms of the number of iterations. Introducing weights makes the optimization relatively slow because of the effect of floating point arithmetic on the numerical stability. As the weights must be determined up to a propto this day ity factor, one can consider scaling them up or rounding / truncating them.

Glacier monitoring

Motivation

The motivation to apply our pipeline on alpine glacier monitoring is driven by the progressive awareness to climate changes as it is an indicator of the local effects of global warming. Glaciers are indeed thinning and melting after being subjected to a rise in the surrounding temperature and a spectacular retreat has been observed on most of the alpine glaciers in recent decades. The problem is that, to this day, it is still difficult to perform regular ground-based surveys of alpine glacier flow: only once or twice a year and only on a few spatial points due to their relatively unreachable location and dangerous geomorphology (crevasses). Remote studying of glaciers can thus lead us to a more complete understanding and modelling of global warming (along with meteorological data). On the other hand, studying the melting of glaciers is also useful for economic and security reasons. Economic; because they also are water resources. Security; because the melting can cause hazardous events to the surrounding areas, such as destruction by seracs falls, or the so-called "glacier lake outburst flood" events, which occur when glacier lakes start to drain through a subglacial channel underneath the glacier causing flood in the downstream valleys. In this section we will show that daily surface measurements are possible via optical remote sensing if one can tackle the rather difficult radiometry. These surface measurements can later be used for a better understanding of the glacier's rheological [START_REF] Murray | Dynamics of glacier calving at the ungrounded margin of helheim glacier, southeast greenland[END_REF] and morphological [START_REF] Rabatel | Estimation of glacier thickness from surface mass balance and ice flow velocities: A case study on argentière glacier, france[END_REF] parameters, and predict hazardous events [START_REF] Schwalbe | Recognition of drainage tunnels during glacier lake outburst events from terrestrial image sequences[END_REF].

Instrumentation and Set-up

The Argentiere glacier (45 • 55'N, 6 • 57'E) is a temperate alpine glacier located in the Mont-Blanc massif in France. It is ∼ 10km long with a surface of ∼ 12km 2 , and between altitudes 3400m (background) and 1600m (snout or Lognan seracs fall). Its surface flow has been studied numerous times by glaciologists [START_REF] Vincent | Glacier thickening and decay analysis from 50 years of glaciological observations performed on glacier d'argentière, mont blanc area, france[END_REF] and remote sensors [START_REF] Ponton | Observation of the argentière glacier flow variability from 2009 to 2011 by terrasar-x and gps displacement measurements[END_REF] [START_REF] Benoit | Multi-method monitoring of Glacier d'Argentière dynamics[END_REF]. Given the ground constraints and the expected motion the choice of the acquisition system is made (camera sensor size, lens parameters range, disparity between the stereo views for depth maps computation, etc.). After the parametrization of the external (orientation, tilt and distances from the area of interest) and internal (consensus between focal length, depth of field and light dosage) camera parameters and finally, the instrumentation of the acquisition system (number of automatic images per day giving the object known velocity, electric power supply, etc.), the system gives us a stack of images taken automatically that will be inquired offline. To observe the largest area, without a "fish eye" effect and with the best resolution that we can embed in the installation, we chose a Nikon D810 with 20mm lens (f/2.8D) parametrized to have an infinite depth of field. It gives us 36.3-megapixel images with a high ground resolution, to observe as many details as possible, but that makes code optimization needful. The expected inconveniences are indeed prominent in the raw data (see figure 14). The images also show a complex and repetitive texture pattern due to crevasses and cracked ice on the glacier zone. As stated before (section 3.3.4), the flow direction is mainly parallel to the image plane.The choice of the time step was made regarding the temporal closure range: The time separation between images within the temporal closure range must be long enough to increase the signal (the flow of the glacier) but short enough to preserve the features for tracking. A long time separation can indeed decorrelate the images to co-register and match for 3 main reasons: First, as stated before, photometric change, second; features on the glaciers may change due to melting, snowfall or changes caused by windblown snow. And last, good features for correlation, such as the crevasses and the ice foliation (mainly transverse), appear from time to time at the same places and with similar shapes. For example, crevasses open regularly upstream of a serac fall. For this reason, and as the Argentiere Glacier is a temperate glacier with a relatively high velocity, the time step between images in the timelapse is 2 days (taken at 12AM to avoid drop shadows as much as possible) and the temporal closure range was chosen empirically to be within 10 images. This set-up makes the maximum time step between matched images (observations for the regularization) equal to 20 days. The time span of the presented results is from July 06 to December 06 2017. The computed depth map has an average uncertainty of ± 10 m compared to ground truth measurements. The results obtained on the Argentiere site (Figure 14) are illustrated with figure 15b encoding the orientation and magnitude of the mean OF map in the HSV (Hue, Saturation, Value) colorspace. Images 27, 28, 37, 38, 57, 60, 61, 62, 63, 64, 72, 73 and 74 have been discarded from preprocessing mainly because of fog, clouds or droplets on the camera screen. These results were obtained by minimizing a damped Least-squares LF temporal closure error such as in Equation (16) using the LSQR solver. Data parallelism was performed to cope with the amount of pixels to process. A zoom on the results for the front and back zones of the glacier (respectively yellow and red rectangles in figure 14b are shown in figures 16 and 17. x and y components of the mean OF map are shown separately. We can clearly see the effect of the temporal regularization. On the static-zones, outliers caused by environmental changes, such as occlusion by clouds, have been removed. On the glacier zone, the mean OF maps shows a river-like flow with more velocity at the center of the glacier caused by the friction on the rough edges within the bedrock. Figures 16f and17d show OF maps more consistent with the longitudinal velocity and slope profiles presented in [START_REF] Trouve | Combining airborne photographs and spaceborne sar data to monitor temperate glaciers: Potentials and limits[END_REF] Spatially: we can see from figures 15b, 16, 17 and 19 that the amplitude of the surface velocity increases with the slope of the glacier and decrease with the altitude. This is in accordance with other velocity maps [START_REF] Rabatel | Estimation of glacier thickness from surface mass balance and ice flow velocities: A case study on argentière glacier, france[END_REF][62] [START_REF] Berthier | Surface motion of mountain glaciers derived from satellite optical imagery[END_REF][64] on the glacier d'Argentiere and with a downslope flow of a fluid. This dynamic is valable for glaciers all over the world (example of Greenland glacier in [START_REF] Murray | Dynamics of glacier calving at the ungrounded margin of helheim glacier, southeast greenland[END_REF]). Temporally: we see from figure 18 that the horizontal flow velocity time series don't have a noticeable trend at a given spatial point and during the inquired period of time (06 June to 06 December 2017). This can also be seen for the sliding velocity of the same glacier in [START_REF] Vincent | Sliding velocity fluctuations and subglacial hydrology over the last two decades on argentière glacier, mont blanc area[END_REF] (figure 8.c) during the year 2013 along with further analysis of sudden accelerations peaks. The obtained results show a mean (spatial and temporal) horizontal velocity of 18 cm/day at the back of the glacier, as seen by our cameras, of 35 cm/day near the front of the glacier and up to 66 cm/day in the seracs fall zone. At an altitude of 2400 m (cross section 4 in [START_REF] Vincent | Sliding velocity fluctuations and subglacial hydrology over the last two decades on argentière glacier, mont blanc area[END_REF] and [START_REF] Vincent | Glacier thickening and decay analysis from 50 years of glaciological observations performed on glacier d'argentière, mont blanc area, france[END_REF]) the mean velocity is 22 cm/day. We reported published surface velocities over the Argentière Glacier on this cross-section since 1996 in table 11 and we can notice the temporal variability over time along with a general negative trend and local fluctuations. These temporal fluctuations in the surface flow velocity time series can also be seen in long mono-method time series such as in [START_REF] Vincent | Glacier thickening and decay analysis from 50 years of glaciological observations performed on glacier d'argentière, mont blanc area, france[END_REF] (figure 4). In [START_REF] Vincent | Sliding velocity fluctuations and subglacial hydrology over the last two decades on argentière glacier, mont blanc area[END_REF], the authors state that the negative trend is due to glacier thinning (weight loss) and local accelerations can be explained by a rise in temperature and episodes of precipitation. The inquired period (June -December) presents indeed higher temperature averages and several precipitation episodes (the rainy season began in August 2017). The semiannual seasonality between dry/cool and hot/rainy seasons can be seen in [START_REF] Vincent | Sliding velocity fluctuations and subglacial hydrology over the last two decades on argentière glacier, mont blanc area[END_REF] (figure 3). We can see that years 2003 and 2017 have similar mass balances [START_REF] Vincent | A nonlinear statistical model for extracting a climatic signal from glacier mass balance measurements[END_REF] and velocity profiles [START_REF] Berthier | Surface motion of mountain glaciers derived from satellite optical imagery[END_REF] and are averaged over roughly the same period of the year, with a higher cumulative precipitation for year 2017. None of the referenced works over the Argentière Glacier gave spatially precise, temporally regularized time series of the surface velocity. These results can be used, for example, to estimate the glacier basal sliding or the ice thickness [START_REF] Rabatel | Estimation of glacier thickness from surface mass balance and ice flow velocities: A case study on argentière glacier, france[END_REF] [START_REF] Vincent | A nonlinear statistical model for extracting a climatic signal from glacier mass balance measurements[END_REF], which is of interest for potential hydro-power. These experimental results highlight the potential of remote sensing via time-lapses, on the one hand, and the necessity to further regularize the automatically obtained results, on the other hand. As already highlighted section 3.3.4, numerous, precise and accurate in situ measurements are necessary to decrease the uncertainty related to the scaling to object coordinate and to georeference the observations our depth maps have an uncertainty of ±10m which can be propagated to the presented VTS.

Results on the Argentiere observation site

Conclusion and future work

We have proposed a processing chain for accurate displacement time series generation from time-lapse images for complex geophysical phenomena monitoring. The pipeline includes robust pairwise OF maps computation and the formulation of a temporal closure error from a network of optical flow fields. Temporal closure mappings are formulated as an inverse problem where a weighted Least Squares solution is sought. Due to eventual ill-conditioning of the linear forward problem and the presence of noise in the observations, the inversion requires some extra "regularization" in order to generate physically plausible solutions. For this, a damping term was introduced to the cost function. Several formulations, weighting strategies, norms and solvers have been tested on simulated datasets to demonstrate the efficiency of the proposed approach under controlled and realistic conditions and to establish an automatic resolution framework given the a priori information on the observations. The proposed automatic method has been successfully applied to monitor the Argentiere glacier flow resulting in more complete and accurate OF maps time series all over the observed scene during a time span of 6 months. In future works, 2 extensions of the pipeline will be investigated; (1) replacing the manually constructed Gabor filter bank by learned convolutionnal filters, and (2) a scene flow [START_REF] Vedula | Three-dimensional scene flow[END_REF] processing for deriving 3D motion fields from stereo regularized optical flow fields.

  (C) Optical flow (OF) fields ⃗ F = (F x , F y ) at every pixel are then computed between all the co-registered images within the chosen temporal closure range to form the Network of Optical Fields (NOFF) (D) The NOFF contain redundant information and makes it possible to minimize a temporal closure error to regularize the resulting VTS and PTS and interpolate values where observations are missing. The matrices in (C) and (D) are examples of observations in one pixel before and after regularization. The master image in figure 1 and all the illustrative examples of the paper being the first valid image of the time-lapse; I 0 .

Figure 1 :

 1 Figure 1: Proposed pipeline. (A) data acquisition. (B) data selection via texture detection. (C) matching and optical flow computation. (D) regularization of the obtained time series via the construction of a network of optical flow fields.

  (a) Index graph for a full range problem: 18 vertices (size of F est ) and 342 directed edges (size of F obs ). In blue all observations including image I 18 . (b) Index graph for a problem with a reduced range: 18 vertices and 160 directed edges.

Figure 2 :

 2 Figure 2: Polygonal index graph for a full-rank problem.

Figure 3 :

 3 Figure 3: 3a: Respectively, mappings from the data space D to the model spaces M LF (left) and M LF (right) and thair histogram of singular values in figure 3c and 3d .The associated condition numbers are respectively 12.06 and 4.35.The rows and columns are arranged as in equations (11) and (12) but only 1 on 20 homologous indexes are displayed on the vertical axis.

  (a) matrix associated with the model in Equation (7), size = 160 × 18: A LF (b) matrix associated with the model in Equation (8), size = 160 × 18: A CM (c) Singular values (rounded) of the linear mappings in figure 4a (top) and 4b (down).

Figure 4 :

 4 Figure 4: 4a: Mapping from the data space D to the model space M LF and its histogram of singular values in figure 4c (top). 4b: Mapping from the data space D to the model space M CM and its histogram of singular values in figure 4c (bottom). The associated condition numbers are respectively 5.01 and 7.79.

Figure 5 :

 5 Figure 5: NOFF index graph for a rank-deficient problem: 18 vertices and 94 directed edges

Figure 6 :

 6 Figure 6: Inverting a random observation vector: 6a: the output of inverting the LF temporal closure mapping is a VTS. 6b: the output of inverting the CM temporal closure mapping is a PTS. Images 2, 3, 14, 16 and 17 are discarded. The position of the tracked pixel in I 8 (green dot in figure 6a) is the green area in figure 6b.

Figure 7 :

 7 Figure 7: Illustration of a severe rank deficiency case: only images 0, 4, 10 and 18 are available over the 19 image forming the time lapse. The matrices of the LF (7a) and CM (7b) linear mappings became large (12 × 18) and of rank r = 3. 7c: velocity output of inverting a simulated observation vector F estLF and F estCM converted to a VTS. 7d: position output of inverting a simulated observation vector F estCM and F estLF converted to a PTS. .

  (a) Hypothetical uncertainties related to observations ↓ (b) Uncertainties related to the estimated VTS (top) and PTS (bottom) through the LF (left) and CM (right) inversions.

Figure 8 :

 8 Figure 8: Illustration of uncertainty propagation through inversion for the hypothetical data described in 3.4.2.

Figure 9 :

 9 Figure 9: Example of piecewise affine transform on the dynamic zone. 9a: template image. 9b: affine transform warped image. 9c: piecewise affine transform warped image.

  Figure 10: ξ VTS in terms of the shear with NCC and OF for Dataset 2. The shear increases with the vertical and temporal axes resulting in an contamination of the observations.

  (a) Ground truth OF map (b) OF map with additional white noise (c) Regularized OF map (d) OF VTS: LF and CM solutions at the central pixel of the simulated dataset 4.

Figure 11 :

 11 Figure 11: Iterative LF inversion results after injection of a white noise of standard deviation σ ω = 5 pixel: Top: OF map between Images 0 and 1 of Dataset 4. Bottom: OF velocity pixel temporal tracking.

Figure 12 :

 12 Figure 12: OF VTS: with SVD and COBYLA at the central pixel of the simulated dataset 4. Image 2 has been discarded creating a gap of 2 values on the original VTS. The WCM outputs (PTS) are converted to VTS.

Figure 13 :

 13 Figure 13: Iterative LF inversion results after systematic error injection on the dynamic zone related to image 3. OF map between Images 0 and 3 of Dataset 4.

Figure 14 :

 14 Figure 14: Argentiere site: data and artifact examples. In 14f a zoom is made to show a common artifact due to sunbeam reflections. This artifact does not have constant image coordinates.

Figure 15 :

 15 Figure 15: Result on the Argentiere glacier. 15a: velocity map with and without camera motion compensation. This motion can be up to 20 pixels. 15b: mean velocity map in the HSV colorspace. Flow direction is encoded in the color (Hue) and flow amplitude is encoded in the Value.

Figure 16 :Figure 17 :Figure 18 :

 161718 Figure 16: Results on the Argentiere glacier. Front zoom (in red on figure 14b).

Figure 19 :

 19 Figure 19: Example of regularized velocity map converted to cm. Between 06 and 07 June 2017. The gaps in the maps are due to missing depth data.

Table 1 :

 1 Parameters of the generated datasets.

		Geometric transformation	Temporal Shift
	Dataset	GA	PAG	LT	PT
	1	-	-	x
	2	x		x
	3	x			x
	4		x		x

Table 2 :

 2 . OF errors ξ VTS (2a) and ξ PTS (2b) before and after a weighted LF (WLF) regularization.

	Dataset OF		Leap Frog LF		LF-converted to position
			Direct	Iterative
		NW	WS	WR	WSR	NW	WS	WR	WSR
	1	0.01646 0.00171 0.00130 0.00115 0.00117 0.00289	0.00339 0.00216 0.00135 0.00216
	2	0.08366 0.02005 0.16890 0.01590 0.01583 0.01457	0.06105 0.03312 0.0315	0.02861
	3	0.10486 0.01737 0.01254 0.01097 0.01115 0.01036	0.05623 0.02667 0.02759 0.01846
	4	0.04363 0.02834 0.02970 0.02970 0.02834 0.02115	0.01346 0.01970 0.01750 0.02079
			(a) ξVTS					(b) ξPTS
	Dataset OF	Common Master CM		CM-converted to velocity
			Direct	Iterative
		NW	WS	WR	WSR	NW	WS	WR	WSR
	1	0.05385 0.00339 0.00220 0.00135 0.00214 0.00527	0.00171 0.00131 0.00116 0.00118
	2	0.38395 0.06115 0.03311 0.03149 0.02862 0.03450	0.02008 0.01688 0.01585 0.01585
	3	0.42601 0.05623 0.02675 0.02763 0.02188 0.01850	0.01737 0.01256 0.01099 0.01116
	4	0.09780 0.01349 0.01973 0.01748 0.02058 0.01345	0.02837 0.02834 0.02968 0.02117
			(a) ξPTS					(b) ξVTS

Table 3 :

 3 OF errors ξ PTS (3a) and ξ VTS (3b) before and after a weighted CM (WCM) regularization.

	Dataset	NCC	Regularized NCC	Dataset	NCC	Regularized NCC
	1	0.11597	0.07570	1	0.25562	0.19200
	2	0.63350	0.37497	2	2.07360	0.43204
	3	2.43873	1.19262	3	6.91797	3.83019
	4	1.71092	0.30109	4	3.97034	0.51492
			(a) ξ VTS		(b) ξ PTS

Table 4 :

 4 (NCC) errors ξ VTS (4a) and ξ PTS (4b) time series before and after (iterative -WLF) regularization.

Table 5 :

 5 OF erros ξ VTS before and after a weighted LF regularization with different noise levels.

	WCM

Table 6 :

 6 OF errors ξ PTS before and after a weighted CM regularization with different noise levels.

Table 7 :

 7 OF ξ VTS before and after a weighted LF regularization with different configurations of missing observations.

	WLF

Table 8 :

 8 OF ξ PTS before and after a weighted CM regularization with different configurations of missing observations. (nod = no data)

Table 9 :

 9 OF ξ VTS before and after a weighted LF regularization with different outlier mean amplitude (associated with 1 image).

	WCM

Table 10 :

 10 OF ξ PTS before and after a weighted LF regularization with different outlier mean amplitude (associated with 1 image).

  on the same glacier. Examples of resulting time series for points A, B and the seracs fall zone in figure 15b are shown in figure 18.

		1996	2003	2007	2013	2017
	Velocity (cm.day -1 )	24-41[62]	23-50[63][56]		23(G)[64]
	at altitude				18(O)[64]
	2400m (cross-section 4 in [65][66])	30[65][66]	27[65][66]	18[65][66]	16[65]	22 (us)
					15-30(SX) [64]
	Mass balance
	averaged at altitude between
	2380 -2600m in m.w.e.a -1 [67]	-3.05595	-4.9377857	-2.6223	-3.792	-4.993
	cumulative precipitation in mm.a -1
	(MeteoFrance rainfall archives )	80	40	89	45	57

Table 11 :

 11 Velocities of the Argentiere glacier from different sources. [62]: DinSAR between 10 & 11 March. [63]&[56]: SPOT5 between 23August & 18 Sept. [64]: multi-method, O for optical between 23 September & 7 October, SX for TerrasarX between 23 October & 3 November and G for GPS geocubes measurments (geocube 1018 in the paper) 16 to 21 September. [66][65]: GPS measurments, annual average. Us: optical photogrametry between June & December.

Acknowledgment

This work is supported by the ANR-15-CE23-0012 PHOENIX grant of the French National Agency of Research.