F Kröber

G Fernandez

F Guiotte

F Delerue

T Corpetti

S Lefèvre

LEARNING UAV-BASED ABOVE-GROUND BIOMASS REGRESSION MODELS IN SPARSE TRAINING DATA ENVIRONMENTS

Keywords: vegetation biomass regression, random forest, deep learning, semi-supervised learning, transfer learning, UAV

This study aims at recovering above-ground biomass information from ultra-high resolution UAV RGB-NIR orthophotos. We focus on a realistic scenario where a limited number of training samples for a landscape with heterogeneous herbaceous vegetation is given. Consequently, we explore different machine learning methods explicitly addressing the limitations of small training samples and compare their predictions quantitatively and qualitatively. Our results show that random forest models perform similarly well to deep learning models. While simpler machine learning models may, therefore, still be preferable, our study also points the way to promising architectures and regularisation techniques for deep learning approaches. Beyond vegetation cover, accurate regression of other variables, including vegetation height, volume and biomass remains a difficult task regardless of the model choice.

INTRODUCTION

Precise estimation of the amount of vegetation present in a given ecosystem is vital for understanding its functioning, managing its resources and predicting its evolution under different constraints including global change. Therefore, many ecological studies measured different indicators of vegetation development (like vegetation cover or vegetation volume), the most precise one being vegetation biomass. However, traditional field campaigns for measuring biomass are time-consuming, costly and destructive, making remote sensing an attractive alternative. Unmanned Aerial Vehicles (UAVs) have proven to be suitable platforms for performing data collection at the field scale with high spatial resolutions in an efficient and nondestructive manner. From a modelling perspective, UAV-based estimations of biomass-related variables can be framed as a regression task between orthophotos and a set of field measurements. The often limited availability of the latter may explain the prevailing use of robust models based on traditional machine learning. While Convolutional Neural Networks (CNNs) have been applied to UAV data for various tasks, their use for solving regression problems [START_REF] Osco | A review on deep learning in UAV remote sensing[END_REF] and especially biomass-related ones [START_REF] Wang | Applications of UAS in crop biomass monitoring: A review[END_REF] is still underrepresented. In this paper, we explore three possibilities for recovering biomass information from ultra-high resolution (i.e., millimetric resolution) UAV RGB-NIR orthophotos. We aim to compare different machine and deep learning methods which are explicitly addressing the limitations of small training samples: Random Forest (RF) (Model A), transfer learning (Model B) and data augmentation (model C). In addition to the above-ground biomass (AGB) itself, vegetation cover, height and volume as related variables are also targeted. It should be noted that in this study we focus on a particularly complex case: our data campaign concerns post-mining stressed ecosystems, with small, heterogeneous, herbaceous plant communities, characterised by low values of AGB. This constraint makes a precise estimation of biomass-related variables more challenging, especially from images with a millimetric resolution.

RELATED WORKS

Today, most studies to estimate biomass still focus on classical parametric learning methods such as simple or multiple linear regressions [START_REF] Wang | Applications of UAS in crop biomass monitoring: A review[END_REF][START_REF] Morais | The use of machine learning methods to estimate aboveground biomass of grasslands: A review[END_REF][START_REF] Poley | A systematic review of the factors influencing the estimation of vegetation aboveground biomass using unmanned aerial systems[END_REF]. While the authors of [START_REF] Chave | Ground data are essential for biomass remote sensing missions[END_REF] have pointed out that onsite measurements are crucial for satellite-borne biomass missions, the lack of in situ measurements also hinders the training of more complex deep learning-driven models when the data is collected by UAVs. So far, preprocessing, including the creation of 3D models to derive structural information, is frequently carried out to allow the usage of simpler models based on a few samples. For example, in [START_REF] Bendig | Estimating biomass of barley using crop surface models (CSMs) derived from UAV-based RGB imaging[END_REF], the authors first calculate a crop surface model from UAV data, then derive the plant height, and finally use regression models from plant height to biomass. Also, multi-temporal data is used to improve the quality of estimations [START_REF] Sinde-González | Biomass estimation of pasture plots with multitemporal UAV-based photogrammetric surveys[END_REF], along with other complementary information such as meteorological features [START_REF] Jiang | UAV-based biomass estimation for ricecombining spectral, TIN-based structural and meteorological features[END_REF]. Meanwhile, the problem of estimating biomass from monotemporal UAV orthophotos without extensive preprocessing and/or auxiliary data remains an open problem. Our paper is a contribution in this direction.

MATERIALS AND METHODS

Data Description

UAV and in situ measurements were taken for four sites in different altitude zones (1850 to 2100 m) in the Pyrenees, France, from spring to summer 2021 (Fig. 1). Due to a difficult environment for plant growth in the surveyed areas, corresponding plant communities were herbaceous and particularly small (<20 cm). UAVs were equipped with R, G, B, NIR sensors and the obtained images post-processed to derive orthophotos with a spatial resolution of 3 mm. In situ measurements were taken for four variables: fractional vegetation coverage, vegetation volume, vegetation height and AGB. About 30-35 measurements were taken per site leading to 160 measurements. For each point, all variables were recorded based on a 25x25 cm square representative of the 1x1 m area centered around the point. Two datasets have been generated from these data: a labelled dataset and an unlabelled one. For the former, the UAV images were cropped to the in situ measurement geolocations generating tiles with a 1x1 m extent centered around each point measurement, resulting in 160 couples of tiles and in situ measurements. 124 of them were splitted in a spatially stratified random manner into train (75), validation (24) and test (25) sets. Another 36 data points, all drawn from a single site, are used as a spatially independent control site (CMB test site). For the unlabelled dataset, the UAV images have been tiled at 1x1 m with an overlap of 50%, resulting in 18,085 tiles of 256x256 pixels, splitted in train (70%), validation (15%) and test (15%) sets.

Model A: Random Forest Model

Given our low amount of data points, an RF classifier [START_REF] Breiman | Random Forests[END_REF][START_REF] Belgiu | Random forest in remote sensing: A review of applications and future directions[END_REF] with manually engineered features has been used as a robust baseline architecture. First, all tiles of the labelled dataset were segmented with a Felzenszwalb segmentation [START_REF] Felzenszwalb | Efficient graph-based image segmentation[END_REF] (Figure 2). The derived objects were then evaluated with a simple automated tile cover classification, distinguishing between vegetation, bare ground and shadowed areas. To this aim, segmented tiles with mean NDVI values above a given threshold were classified as vegetation. Shadowed areas classification was based on consistently low reflectance values across the R, G and B channels. Given these classified tiles, two types of features were engineered. First, spectral features and their derivates, i.e. the share of vegetation, bare ground and shadow, were used. The shares of vegetation were calculated for three different NDVI thresh-Fig. 2. Schematic representation of the three models. olds (0.3, 0.4, 0.5) to let the RF algorithm choose the most informative variable. Second, the spatial arrangement of these classes was also considered: the size and shape properties of the segmented objects were used. Specifically, the mean value, standard deviation of the size and shape index of all objects classified as vegetation were taken into account as features. In order to represent multi-scale information, the features were calculated for two segmentation scale levels (20,200). This made a set of 40 features derived from the images fed into the RF for the regression of vegetation parameters.

Model B: Transfer Learning

The underlying hypothesis of model B is that an Autoencoder (AE) [START_REF] Bank | Autoencoders[END_REF], trained on the unlabelled dataset, can generate a feature space that can be easily adapted to the regression task on the labelled dataset. To this end, the encoding part of an AE model was designed to have a latent representation of size 32x42x42. This dimensionality of the latent feature space allowed to create a simple regression CNN with less than 10 K parameters in the second step. The regression CNN starts with a 1D convolutional layer intended to reduce the number of feature maps to the ones needed for the task at hand. Subsequently, three 2D convolutional layers are used again. Average pooling is applied to the output, the features are flattened, fed into a linear layer and sigmoid is applied to force the output to take values between 0 and 1. Rescaling (100x) represents the final step to get the common value range of the variables of interest.

Model C: Data augmentation/Semi-supervised approach

Another possibility to exploit the great number of unlabelled tiles is to consider their image similarity to the labelled ones, following the semi-supervised setting from [START_REF] Castillo-Navarro | Semi-supervised semantic segmentation in earth observation: The minifrance suite, dataset analysis and multi-task network study[END_REF]. Therefore, model C (Figure 2) uses a one-shot pixel-wise pseudo-label generator for unlabelled tiles. With a ResNet-18 [START_REF] He | Deep residual learning for image recognition[END_REF], pre-trained on ImageNet and without the last layer, 512 features were calculated for each UAV tile. The proximity of two vectors in this feature space is assumed to indicate the structural similarity of images and can be exploited to assign the pseudo-labels. To circumvent the issues arising from the high dimensionality of the feature space, vectors were mapped to a 2D space with t-Distributed Stochastic Neighbour Embedding (t-SNE) [START_REF] Hinton | Stochastic neighbor embedding[END_REF], before performing the assignment. In detail, the dataset was divided into small batches of 150 tiles, containing an equal number of labelled training tiles and unlabelled tiles. The mean nearest neighbour (NN) distance between all points was calculated and used as a distance threshold for assigning values to unlabelled tiles. If the distance between an unlabelled point and one or multiple labelled ones was lower than the mean NN distance, a pseudo-value was calculated for the unlabelled one using the weighted average of the values of all the labelled tiles being considered. This procedure was repeated for all unlabelled tiles. To handle t-SNE inherent stochasticity, multiple trials were conducted for each batch. The one with minimum Kullback-Leibler divergence and root mean squared error (RMSE) calculated between labelled point pairs less than the mean NN distance apart was chosen to be the best one used for pseudovalue assignment. Afterwards, the datasets were balanced, ensuring a uniform quantile spacing for each variable using a random oversampling technique. Finally, the number of generated pseudo-values was capped at 5000 to limit the amount of data and to speed up the subsequent training. The regression was then performed by placing a regression block (as described in Section 3.3) in place of the last layer of a pre-trained EfficientNet-b0 [START_REF] Tan | Efficientnet: Rethinking model scaling for convolutional neural networks[END_REF].

Validation strategy

For all models, hyperparameter tuning has been performed by measuring the model performance on the validation set. For model A, the number of features considered when looking for the best split at a given tree node and the minimum number of samples required to be at a leaf node were tuned as hyperparameters. They were optimised using the RMSE. For models B and C, the learning rate, optimiser (Adam vs. SGD) and loss (L1 vs. L2) were considered as relevant hyperparameters and chosen based on the mean absolute error (MAE). The stratified sampling splits as described in Section 3.1 were created ten times to perform cross-validation of the results. Only for the AE training within model B, cross-validation was not necessary due to the large number of unlabelled tiles.

RESULTS AND DISCUSSION

The quantitative comparison between the models using MAE, RMSE and Person's correlation coefficient (r) is presented in Table 1. Across all metrics, the three models have similar accuracy. For AGB, model A performs slightly better, while for all other variables model C produced the most accurate results. For each singular variable, the r coefficient is coherent between models. It varies considerably between variables: while predictions and in situ measurements are strongly correlated for coverage, volume and height exhibit intermediate values. For AGB lowest correlations are observed. This reflects the increased difficulty of estimating more complex structural-morphological variables based on 2D UAV data. Interestingly though, the increased complexity of the regression tasks for these structural variables is not handled better by the deep learning models compared to the RF model. A correlation between the models' capacities, the quality of their predictions and the complexity of the task is not apparent. Importantly, a large uncertainty is attached to all metrics introduced by the small training dataset. To investigate the statistical significance of the results, a Bayesian correlated t-test [START_REF] Benavoli | Time for a change: a tutorial for comparing multiple classifiers through Bayesian analysis[END_REF] was conducted. The advantage of this framework is the ability to calculate a posterior distribution of performance differences between the models, which can be interpreted as an actual probability distribution. This allows to evaluate if the observed differences exceed a level for which we would consider two algorithms to obtain practically equivalent results. The regions of practical equivalence (aka ROPEs) were set to 0.05 times the standard deviation of the respective variable of interest. As Figure 3 shows, most of the model comparisons obtain non-significant differences and the models are more likely to differ only for vegetation volume. To evaluate the results qualitatively, predictions for each model were created for the independent CMB test site using a dense sampling interval of 0.5m. Subsequently, the predicted values were linearly interpolated to obtain spatially continuous representations (Fig. 4). As can be seen in these maps, there is a high correlation between the predictions of the models and the different variables. For all variables and models, the predicted distributions generally match the spatial patterns as emergent from the true and false-color visualisations of the site. Differences between the individual models concern a stronger tendency of model C towards high-frequent noise and less accurate predictions of intermediate values especially for vegetation cover and height. At the same time, the predictions across variables seem to be better correlated for model C. While spots may appear to be densely vegetated based on predicted vegetation cover and height, they may counter-intuitively have low biomass and vegetation volume values at these locations for models A and B. Best models for each variable and metric are highlighted.

CONCLUSIONS

In this paper, we proposed learning-based methods to derive biophysical parameters (AGB, vegetation height, volume and coverage) from ultra-high resolution UAV imaging (i.e., millimetric resolution). We explored approaches that are explicitly addressing the limitation of sparse training data: the first one relied on a robust machine learning algorithm, the second one exploited the large quantity of unlabelled UAV data and the last one additionally considered the statistical distribution of the in situ measurements (pseudo-label gen- eration). Strikingly, the considerable uncertainty in predictions did not allow any model to be identified unequivocally as superior or inferior to another. This is another consequence of the lack of in situ measurements but also due to their distribution and intrinsic measurement precision. We note that the study of the uncertainties (via the t-test in our study) is a crucial step to compare the models in such cases, even if a cross-validation approach is already used. Despite the existing uncertainties, it has to be noted that the two deep learning models were not significantly worse than the RF model. The performances of model C were even promising. Even if the literature on AGB's regression from UAV data is still focused mainly on machine learning approaches, this exploratory study was able to pave the way for novel deep learning methods to study biophysical parameters in the case of sparse training data. A transfer to and test with other data sets such as GrowliFlower [START_REF] Kierdorf | Growliflower: An image time-series dataset for growth analysis of cauliflower[END_REF] remains to be done in the future. Finally, from an ecological perspective, the most interesting parameters to study degradated ecosystems are AGB (also the most challenging with height) and vegetation cover. With this in mind, we can affirm that the good performances on the latter variable puts this study in an interesting position for pratical applications.

Fig. 1 .

 1 Fig. 1. Site localisation, an example of UAV images along with the in situ measurements locations and data description: (a) UAV tile in RGB; (b) UAV tile in false-color; (c) distribution of AGB values (g/m 2 ) in the labelled dataset (the green line is the kernel density estimation of the distribution); (d) vegetation height (cm); (e) vegetation volume (dm 3 ); (f) vegetation coverage (%).

Fig. 3 .

 3 Fig. 3. Posterior distributions for mean absolute differences between models. Red bars enclose the ROPEs.

Fig. 4 .

 4 Fig. 4. Mapped predictions for CMB test site.

Table 1 .

 1 Test accuracy measured for the 25 samples of the test set. Mean values and standard deviations across 10-folds are specified.

	Variable Model	MAE	RMSE	r
	AGB	A	160.83 (± 26.99) 220.99 (± 50.36) 0.58 (± 0.17)
		B	176.94 (± 33.29) 234.54 (± 55.87) 0.55 (± 0.18)
		C	164.21 (± 29.09) 235.43 (± 47.00) 0.61 (± 0.17)
	coverage A	8.60 (± 1.43)	11.07 (± 2.06) 0.93 (± 0.03)
		B	9.36 (± 1.13)	11.47 (± 1.36) 0.93 (± 0.03)
		C	8.07 (± 1.84)	10.45 (± 2.94) 0.94 (± 0.04)
	height	A	5.93 (± 0.97)	8.00 (± 1.40) 0.68 (± 0.07)
		B	5.56 (± 1.08)	7.74 (± 1.36) 0.70 (± 0.09)
		C	5.46 (± 1.00)	7.32 (± 1.34) 0.76 (± 0.09)
	volume	A	42.76 (± 7.63) 65.94 (± 14.84) 0.77 (± 0.08)
		B	37.71 (± 9.84) 62.60 (± 16.23) 0.79 (± 0.09)
		C	36.66 (± 8.12) 58.15 (± 12.16) 0.83 (± 0.06)

This work has been supported by the SIXP project, funded by the French National Research Agency (ANR-19-CE02-0013-01).


Mathilde Letard

Thomas Corpetti

Dimitri Lague

BATHYMETRIC LIDAR WAVEFORM DECOMPOSITION WITH TEMPORAL ATTENTIVE ENCODER-DECODERS

Keywords: full-waveform lidar, bathymetric lidar, signal decomposition, waveform decomposition, attention UTime

This paper is concerned with the decomposition of bathymetric lidar waveforms. Because of the presence of water, processing such data remains a challenge since water impacts their shape and signal-to-noise ratio, depending in particular on the associated turbidity. In this paper, we explore the use of attentive autoencoders to decompose bathymetric waveforms and recover their air/water interface, water column, and water bottom components simultaneously, without relying on assumptions about the impulse or target surface nature. On simulated waveforms, the method achieves lower decomposition error than existing approaches, handling overlapping echoes of very shallow waters and weak returns in deeper water. This opens to attractive strategies to process real bathymetric waveforms.

INTRODUCTION

Though airborne LiDAR data are mainly used as 2D products or 3D point clouds, their origin lies in the processing of backscattered laser signals and their conversion into spatial information. The precise extraction of the various signal components linked to the hit targets is thus crucial to interpreting with precision the observed scene. However, this extraction entirely relies on the ability both to process the noise efficiently and to remain sensitive enough to subtle variations denoting low-reflective or vertically close targets. Research on full-waveform lidar data processing revealed the potential of finer analysis to retrieve additional targets and increase the resulting point clouds' density [START_REF] Mallet | Full-waveform topographic lidar: State-of-the-art[END_REF]- [START_REF] Chauve | Advanced full-waveform lidar data echo detection: Assessing quality of derived terrain and tree height models in an alpine coniferous forest[END_REF]. However, handling large datasets while minimizing the level of false detections -which produce noise in the point clouds and complexify their interpretation -is still a challenge, namely for airborne bathymetric lidar (ALB). ALB uses green lasers to penetrate water and produce bathymetric data. The resulting waveforms have two main peaks: one generated by the water surface and another occurring when, and if, the water bottom is hit [START_REF] Philpot | Airborne Laser Hydrography II[END_REF], [START_REF] Lague | Topo-bathymetric airborne LiDAR for fluvial-geomorphology analysis[END_REF]. ALB signals have to compose with additional sources of noise due to the optical scattering of the laser beams in the water column, which 1) exponentially attenuates the reflected intensity of the water bottom and 2) reflects a portion of the signal towards the receptor without signifying a solid target, resulting in a third waveform component that elongates the surface peak [START_REF] Philpot | Airborne Laser Hydrography II[END_REF], [START_REF] Lague | Topo-bathymetric airborne LiDAR for fluvial-geomorphology analysis[END_REF].

In addition, extracting meaningful information from waveforms is not trivial because of the noise that disrupts the signal and gets mixed up with weak information, which is often lost during processing. In this paper, instead of a deconvolution by inversion of the pulse response function that can be challenging due to the noise in the signal, we explore the use of deep neural networks to tackle these problems and perform bathymetric waveform decomposition. We use a temporal attentive autoencoder to reconstruct the three main components of the bathymetric lidar signal simultaneously: the water surface, the water column, and the bottom. Our method targets three main limitations of existing waveform processing methods: (1) it does not rely on an iterative process, and (2) it simultaneously recovers the three main bathymetric waveform components.

WAVEFORM PROCESSING METHODS

The collected waveform is a sum of sub-signals that each result from a convolution between the emitted pulse and the target surface function [START_REF] Chauve | Processing Full-Waveform Lidar Data: Modelling Raw Signals[END_REF]. Traditionally, waveform processing is thus made either with deconvolution methods [START_REF] Zhou | Gold -A novel deconvolution algorithm with optimization for waveform LiDAR processing[END_REF], [START_REF] Wu | A comparison of signal deconvolution algorithms based on smallfootprint LiDAR waveform simulation[END_REF] or decomposition procedures relying on the fitting of mathematical functions to previously detect peaks in the signal [START_REF] Chauve | Processing Full-Waveform Lidar Data: Modelling Raw Signals[END_REF], [START_REF] Zhou | Comparison Analysis of Five Waveform Decomposition Algorithms for the Airborne LiDAR Echo Signal[END_REF]. These approaches rely on various assumptions about the components' impulse function or nature. In [START_REF] Wagner | Gaussian decomposition and calibration of a novel small-footprint full-waveform digitising airborne laser scanner[END_REF] and [START_REF] Zwally | ICESat's laser measurements of polar ice, atmosphere, ocean, and land[END_REF], the authros show that waveforms can be considered as a sum of Gaussians corresponding to the objects in the illuminated cone. The emitted pulse and the target functions are thus often approximated with Gaussian functions of fixed width or more asymmetric shapes to fit the empirical sensor impulses [START_REF] Chauve | Processing Full-Waveform Lidar Data: Modelling Raw Signals[END_REF], [START_REF] Hofton | Decomposition of laser altimeter waveforms[END_REF], [START_REF] Jutzi | Range determination with waveform recording laser systems using a Wiener Filter[END_REF]. When using the decomposition approach, Gaussian mixture models are often used to estimate iteratively the parameters of these Gaussian or other predefined functions, which are then used as proxies on the nature of the targets [START_REF] Niemeyer | Conditional random fields for urban scene classification with full waveform LiDAR data[END_REF]. Deconvolution procedures include, among others, the Richardson-Lucy deconvolution and the Wiener filter deconvolution. Several comparative studies [START_REF] Wu | A comparison of signal deconvolution algorithms based on smallfootprint LiDAR waveform simulation[END_REF], [START_REF] Wang | A comparison of waveform processing algorithms for singlewavelength LiDAR bathymetry[END_REF], [START_REF] Parrish | Empirical comparison of full-waveform lidar algorithms: Range extraction and discrimination performance[END_REF] show the superiority of Richardson-Lucy deconvolution in retrieving weak returns and recovering the amplitude of the different components but underline its costly iterative process. Several limitations appear with these classical approaches:

(1) Decomposition methods relying on the fitting of mathematical functions are not always adapted to the presence of the water column component, which has an entirely different shape; (2) Target surface functions may not always be similar/gaussian, and returns may be poorly approximated or mixed; (3) The denoising methods or post-processing procedures used with these approaches to avoid false detections may lead to loss of information.

PROPOSED DECOMPOSITION NETWORK

We adapt U-Time, the 1D version of U-Net introduced in [START_REF] Perslev | U-Time: A Fully Convolutional Network for Time Series Segmentation Applied to Sleep Staging[END_REF] and applied to bathymetric lidar waveforms in [START_REF] Letard | Classification of coastal and estuarine ecosystems using full-waveform topo-bathymetric lidar data and artificial intelligence[END_REF]. We use the resulting encoder-decoder structure to output three waveform components, as illustrated in Figure 1.

The encoder is a stack of 4 identical blocks comprising two series of 1D convolution, rectified linear unit activation, and batch normalization. Each block is followed by a max pooling operation of size 2, thereby reducing by half the length of the feature vector passed on to the next step. In addition, at each step, the number of convolution filters, initially set to 16, is doubled to increase the network's abstraction ability.

The decoder is a stack of four blocks with a different structure. Each block starts with a cross-attention mechanism, followed by up-sampling of size 2, a typical U-Net skip connection, and two series of 1D convolution, rectified linear unit activation, and batch normalization. Contrary to the encoder, the number of filters at each block is reduced by half, thus returning to 16 in the last convolution layers.

Additive attention, as introduced in [START_REF] Bahdanau | NEURAL MACHINE TRANSLATION BY JOINTLY LEARNING TO ALIGN AND TRANSLATE[END_REF], is used to better pass information about the waveform across the network and allow the model to focus on features from different representation spaces in the decoder, thereby reconstructing waveform components more effectively. Using this setup, each sample of the decoder layer attends to each sample of the corresponding encoder layer.

Between the encoder and the decoder, two series of 128 convolution filters with ReLu activation and batch normalization are applied to the 16 samples-long embedding of the input waveform. At the end of the decoder, three dense layers with sigmoid activation allow the network to output three sequences corresponding to the three waveform components.

To constrain the network to learn both amplitude and shape relevance when reconstructing the waveform components, we designed a custom loss function, relying on the Mean Absolute Error (MAE) and the Bhattacharyya distance.

𝑙𝑜𝑠𝑠 = % 1 𝑁 ) 𝑎𝑏𝑠(𝑝𝑟𝑒𝑑 -𝑡𝑟𝑢𝑒) × (1 + 𝑡𝑟𝑢𝑒) ! "#$ 7 + 81 -) 9 𝑝𝑟𝑒𝑑 ∫ 𝑝𝑟𝑒𝑑 × 𝑡𝑟𝑢𝑒 ∫ 𝑡𝑟𝑢𝑒 ; ( 1
)
where 𝑁 is the number of samples in the waveform, 𝑝𝑟𝑒𝑑 is the predicted sequence, and 𝑡𝑟𝑢𝑒 is the true sequence. Roughly, the idea of the first term is to ensure global consistency in the estimations, this consistency being balanced by the waveform intensity (more weight is applied for locations with high intensities); the second term rather focuses on shapes of the decomposed waveforms. Training is performed on batches of 100 waveforms using an Adam optimizer with a constant learning rate of 0.0001. The model was trained for 500 epochs. One training epoch took, on average, 32 seconds on an Nvidia Quadro RTX 5000 GPU.

BATHYMETRIC WAVEFORMS SIMULATION

To simplify the training process (that requires pairs of input waveforms and associated decompositions), to evaluate the approach, and to fully control the physical conditions and their output, we use simulated data obtained with a bathymetric waveform simulator [START_REF] Abdallah | Wa-LiD: A new LiDAR simulator for waters[END_REF]. For training, validation and testing, we simulated three sets of 15000, 3700, and 10000 waveforms, respectively. In each dataset, five types of environments were simulated by varying the simulation parameters (the water body's depth and turbidity). Different impulse functions were used to generalize the method to different sensors and a systematic noise is added. The parameters used are summarized in Table 1. Table 1. Parameters used to simulate lidar waveforms. Figure 1. Illustration of the proposed waveform decomposition network architecture.

EVALUATION METRICS AND COMPARISON

To evaluate the performance of our method, we compare each retrieved component to its simulated equivalent. For each of the three components, three metrics are computed. They are defined in equations ( 2) -( 4), where 𝑁 is the waveform length, 𝑝𝑟𝑒𝑑 the predicted sequence, 𝑡𝑟𝑢𝑒 the true sequence, and 𝐾𝐿 the Kullback-Leibler divergence.

The weighted MAE, quantifies the mean element-wise amplitude error:

The Kullback-Leibler divergence measures the relative entropy, or difference in information that two distributions contain. As it is not symmetric, we use a custom metric:

The Bhattacharyya distance quantifies how far two samples of two distributions are from each other in time and amplitude. For consistency purposes, we used the metric B': These three metrics range between 0 and 1, where 0 indicates high similarity and 1 indicates high dissimilarity.

We compare our results to Gaussian Mixture Model (GMM) decompositions with two components -as the water column is not Gaussian, and as in [START_REF] Wang | A comparison of waveform processing algorithms for singlewavelength LiDAR bathymetry[END_REF] -and a Lucy-Richardson (RL) deconvolution. Other methods will be tested in the future. Table 3. Results obtained with state-of-the-art approaches

RESULTS

Results obtained with the best model and with the classical approaches on the test data are presented in Tables 2 and3.

Overall, the autoencoder reaches lower wMAEs than the GMM and RL methods, although it has higher amplitude errors than distribution distances, which shows that it tends to estimate component shapes better than their amplitude.

Analyzing the metrics more deeply shows that there is no relationship between the value of wMAE, mKL and B' and the turbidity of the simulated water body. It also highlights the higher decomposition errors when processing waveforms at depth below 0.5 m. Qualitative results, presented in Figure 2 illustrate these observations: in Figure 2.a), which illustrates very shallow waters, the gap between the prediction and the simulation is broader than in Figures 2.b) and 2.c), and the amplitude error is noticeable, while the peaks appear close to their true locations. Figure 2.c) shows the ability of the network to recover very weak returns in deeper waters. Visual analysis of the corresponding GMM and RL results show that RL misses the bottom return in the cases presented in Figure 2.a) and mistakes noise for returns in the 2.c) one, while GMM with two components retrieves two returns in the first case, and falsely decomposes the surface return into two echoes in the second.

DISCUSSION AND CONCLUSION

Our results outline the potential of temporal attentive convolutional neural networks for lidar waveform Compared to GMM and RL it is possible to reconstruct the three components of bathymetric waveforms simultaneously. Amplitude errors are lower than with these two traditional methods, and overlapping or weak echoes not captured with GMM and RL are detected.

𝑤𝑀𝐴𝐸 = 1 𝑁 ) 𝑎𝑏𝑠(𝑝𝑟𝑒𝑑 -𝑡𝑟𝑢𝑒) × (1 + 𝑡𝑟𝑢𝑒) ! "#$ (2) 𝑚𝐾𝐿 = 𝐾𝐿 , 𝑝𝑟𝑒𝑑 ∫ 𝑝𝑟𝑒𝑑 , 𝑡𝑟𝑢𝑒 ∫ 𝑡𝑟𝑢𝑒 / + 𝐾𝐿( 𝑡𝑟𝑢𝑒 ∫ 𝑡𝑟𝑢𝑒 , 𝑝𝑟𝑒𝑑 ∫ 𝑝𝑟𝑒𝑑 ) 2 (3)
𝐵 ! = 1 -5 6 𝑝𝑟𝑒𝑑 ∫ 𝑝𝑟𝑒𝑑 × 𝑡𝑟𝑢𝑒 ∫ 𝑡𝑟𝑢𝑒 (4)
Considering the correlation between the physical characteristics of the targets and their response functions, decomposing lidar waveforms into a sum of surface functions is of great potential to derive information on the nature or the structure of the objects detected. In the case of bathymetric lidar, it offers the potential to retrieve more information about the bathymetry. In addition, it may pave the way for water column turbidity assessments relying on remotely sensed rather than field measurements. Further developments to add physical constraints in the loss function and better estimate the maximum signal-to-noise ratio supported by the approach are being conducted. Application to real data should also give further information on the usability of attentive autoencoders to process large sets of bathymetric waveforms.

Figure 2 .

 2 Examples of waveform decomposition results in a) very shallow, b) shallow and c) deep water.

Table 2 .

 2 Results obtained with the proposed approach. decomposition. Using deep neural networks has the advantage of not necessitating setting fixed parameters.

	Waveform	Median	Median	Median
	component	wMAE	mKL	B'
	Total	0.00333	0.00056	0.00028
	Surface	0.00420	0.00040	0.00020
	Column	0.00302	0.00056	0.00028
	Bottom	0.00333	0.00074	0.00037
	Waveform component	GMM	RL
	Total		0.235	0.122
	Surface		0.147	/
	Column		/	/
	Bottom		0.098	/

ACKNOWLEDGEMENTS

The authors are grateful to the Saur Group, the Région Bretagne and the Nantes-Rennes Topo-Bathymetric LiDAR platform for their financial support.


Dmytro Kushnir

Céline Ruscher

Eckhard Bartsch

Fabrice Thalmann

Pascal Hebraud

Pascal Hébraud

Stress overshoot, hysteresis, and the Bauschinger effect in sheared dense colloidal suspensions

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.

I. INTRODUCTION

Materials deform or flow under the application of stress. Many exhibit both behaviors and their response depends on the amplitude of the applied stress, and very often it depends also on the past history of the materials. A common feature of crystalline materials and amorphous systems, such as foams, emulsions or concentrated suspensions is the existence of a transition between an elastic response at small deformation amplitudes to a plastic response, where the system flows, at larger amplitudes. In general, under an increasing applied deformation, the stress first linearly increases, and reaches a steady state value at large deformations, but goes through a maximum value before it reaches its steady-state value. This stress overshoot is a common feature of many amorphous materials, such as emulsions [START_REF] Papenhuijzen | The role of particle interactions in the rheology of dispersed systems[END_REF], foam [START_REF] Khan | Foam Rheology: III. Measurement of Shear Flow Properties[END_REF] suspensions [START_REF] Nagase | Heterogeneous Behavior after Yielding of Solid Suspensions[END_REF] and its amplitude depends on the increase rate of the deformation [START_REF] Dzuy | Yield Stress Measurement for Concentrated Suspensions[END_REF].

Numerical simulations of binary Lennard-Jones glasses suggest that the overshoot value, σ over , scales with the applied shear rate γ as σ over ∝ ln γ [START_REF] Varnik | A study of the static yield stress in a binary lennardjones glass[END_REF][START_REF] Rottler | Unified description of aging and rate effects in yield of glassy solids[END_REF]. Conversely, numerical simulations of attractive spherical particules exhibit a power law behavior, σ ∝ γ0.5 [START_REF] Whittle | Stress overshoot in a model particle gel[END_REF] and neutralized carbopol microgel suspensions exhibit a weaker dependency σ over ∝ γ0.13 [START_REF] Divoux | Stress overshoot in a simple yield stress fluid: An extensive study combining rheology and velocimetry[END_REF].

Moreover, experiments have shown that attractive colloidal suspensions exhibit two distinct yield points, due to the breaking of inter-cluster bonds and to collisions between clusters respectively [START_REF] Koumakis | Two step yielding in attractive colloids: Transition from gels to attractive glasses[END_REF]. The second yield point is characterized by a stress overshoot whose value depends on the applied shear rate, as σ over ∝ γ0.5 , the exponent being almost independent of the concentration.

The yielding mechanisms of hard sphere suspensions are a priori different as they do not involve the breaking of attractive bonds between particules, but the escape from cages is the main mechanism of yielding. Yield is related to the fact that the Brownian motion is not sufficient to relax the structure formed under stress and that the suspension is able to store stress before yielding [START_REF] Koumakis | Yielding of Hard-Sphere Glasses during Start-Up Shear[END_REF][START_REF] Zausch | From equilibrium to steady state: the transient dynamics of colloidal liquids under shear[END_REF].

Moreover, at large deformation amplitudes, materials very often exhibit a response that depends on their deformation history. In a typical experiment, a pre-shear is applied, consisting in a constant or oscillatory deformation at high shear rate or high amplitude, and the mechanical property of the system is measured at some time t w after the end of this so-called rejuvenating protocol. For instance, the stress overshoot increases with the time elapsed since the rejuvenation, t w . Past history dependence of the mechanical response of a system is a common feature that is well known in the engineering literature. Thus, metals exhibit a smaller resistance to compression when they have been previously submitted to tension. This effect has been discovered by Bauschinger when studying steels [START_REF] Bauschinger | Uber die veränderung der elasticitätsgrenze und des elasticitätsmoduls verschiedener metalle[END_REF][START_REF] Sowerby | A review of certain aspects of the bauschinger effect in metals[END_REF]. A similar effect has been observed in atomistic simulations of amorphous systems [START_REF] Falk | Dynamics of viscoplastic deformation in amorphous solids[END_REF][START_REF] Karmakar | Statistical physics of elastoplastic steady states in amorphous solids: Finite temperatures and strain rates[END_REF], and discussed in the athermal limit, in terms of the elementary irreversible deformations that occur in such systems, in Shear Transformation Zone( STZ) [START_REF] Falk | Dynamics of viscoplastic deformation in amorphous solids[END_REF]. In that context, Patinet et al. [START_REF] Patinet | Origin of the Bauschinger Effect in Amorphous Solids[END_REF] have shown the importance of the role of the local stress barrier anisotropy to understand the origin of the Bauschinger effect in the deformation of athermal glasses.

More generally, the Bauschinger effect is related to hysteresis behavior of the stress when the applied deformation is cyclic and depends on the nature of the materials. It may, moreover, present some other specificities, such as a drastic reduction of the elastic response after the first cycle in filled elastomers (the so-called Mullins effect). As far as concentrated suspensions are concerned, they exhibit a reduction of their viscosity after the application of a high shear rate, a phenomenon called thixotropy. The behavior may be seen as a hysteresis response if one considers the stress vs shear rate behavior : the flow curve is measured upon increasing the applied shear rate and then upon decreasing the shear rate. These two flow curves do not superimpose and the area enclosed by the hysteresis loops decreases when the applied shear rate decreases. Experiments in carbopol suspensions [START_REF] Divoux | Rheological hysteresis in soft glassy materials[END_REF] have shown that this behavior is also related to the appearance of shear banding.

Here, we study experimentally the nonlinear flow behavior of a concentrated suspension of spheres which mimic hard sphere behavior reasonably well. We consider mainly the two behaviors described above : stress overshoot and non-reversibility of the flow, that manifests itself by a hysteresis behavior and a change of the apparent elastic modulus upon flow reversal. Moreover, we use numerical simulations to compute the anisotropy of the structure of the suspension submitted to the same deformation histories and we study how the response is related to the local microstructure.

II. SYSTEM A. Experimental system

We use spherical particles of cross-linked polymer chains, called microgels, that are swollen in a good solvent. Several such systems have been developed, in particular Nisopropylacrylamide particles cross-linked with N,N'-methylenbisacrylamide, whose swelling ratio in water may be controlled by the temperature [START_REF] Schild | Poly(N-isopropylacrylamide): Experiment, theory and application[END_REF][START_REF] Senff | Temperature sensitive microgel suspensions: Colloidal phase behavior and rheology of soft spheres[END_REF]. In the present study, nonthermosensitive PS microgel particles are used [START_REF] Beyer | From nuclei to micro-structure in colloidal crystallization: Investigating intermediate length scales by small angle laser light scattering[END_REF][START_REF] Eckert | Glass transition dynamics of hard sphere like microgel colloids with short-ranged attractions[END_REF]. Their struct. P150-50 and dynamical properties [START_REF] Eckert | The effect of free polymer on the interactions and the glass transition dynamics of microgel colloids[END_REF][START_REF] Kegel | Space-resolved dynamic light scattering probing inhomogeneous dynamics in soft matter[END_REF] have been studied in detail and closely match that of hard sphere suspensions [START_REF] Schneider | On tuning microgel character and softness of cross-linked polystyrene particles[END_REF]. They have been used to study glass transition behavior [START_REF] Eckert | Re-entrant Glass Transition in a Colloid-Polymer Mixture with Depletion Attractions[END_REF] and crystallization kinetics [START_REF] Palberg | Unusual Crystallization Kinetics in a Hard Sphere Colloid-Polymer Mixture[END_REF][START_REF] Iacopini | Crystallization kinetics of polydisperse hard-spherelike microgel colloids: Ripening dominated crystal growth above melting[END_REF][START_REF] Franke | Solidification of a colloidal hard sphere like model system approaching and crossing the glass transition[END_REF]. We use bidisperse suspensions of particles with crosslink ratio 1 : 50 (i.e. there is on average 1 crosslink per 50 monomer units) and radii of 175 nm and 150 nm with number ratio N = 3.1 (3.1 small particles per 1 large particle). The system exhibits a glass transition at ϕ = 0.58 [START_REF] Burger | Influence of the polymer size on depletion attraction-induced gel and glass transitions of microgel colloids[END_REF].

B. Rheological measurements

An Anton-Paar stress controlled rheometer MCR-301 is used in order to perform rheological measurements. All measurements are done in a cone/plate geometry of radius 12.5 mm and angle 2 • . In a typical experiment, the sample is first rejuvenated by the application of an oscillatory deformation of amplitude γ = 400% and frequency 1 Hz, for 300 s. Then, the system is let at rest for t w = 600 s and the measurements are started. During all the experiments we follow a deformation protocol with a quasi-constant absolute value of the shear rate. In order to avoid inertia effects and to make sure that each measurement is corresponds to an equilibrium stress value, we apply a deformation protocol that follows a series of steps. The strain is incremented or decremented by a constant value, δγ = 5 • 10 -3 , separated by a constant time δt (Fig. 1(a)). An example of stress evolution after the application of strain step is given in Fig. 1(b). The rheometer response reaches equilibrium after time of the order of 3 • 10 -1 s, and stress measurements are performed for times longer than 5 • 10 -1 s. We then define the effective shear rate as γ = δγ/δt whose absolute value is kept constant during an experiment. Given the time of stress equilibration after a strain step, the maximum shear rate available is 10 -2 s -1 . The maximal deformation amplitude γ max is varied between 0.03 and 0.2.

C. Simulations

In order to get the suspension microstructure, we perform molecular dynamic (MD) simulations using LAMMPS. We consider a binary mixture of particles, denoted A and B, that interact through an inverse power law (IPL), potential. Experimenal measurements of the evolution of the elastic modulus with the volume fraction allows to infer the interaction potential between the particles, that is well described by [START_REF] Schneider | On tuning microgel character and softness of cross-linked polystyrene particles[END_REF] :

U αβ (r ij ) = k B T σ αβ r ij 42 (1)
where α, β = A, B and r ij s the distance between two particles. We choose σ AA = 1 and σ BB = 0.876. Assuming that the interactions are additive, the cross characteristic length scale is equal to σ AB = 0.938. Particles are subjected to an overdamped Langevin dynamics at T = 1.0 with a damping parameter ξ = 4.0. The time step is fixed to δt = 10 -3 . The shear rate value is determined so that the experimental and numerical Péclet numbers are identical. We have :

P e = γτ B = γ R 2 D = γ Rξ k B T (2)
The experimental Péclet number is equal to 3 • 10 -3 , and the numerical shear rate is chosen in order to match the Péclet number obtained with Eq. 2, taking into account the largest particles, of radius R AA = 0.5. This leads to the numerical shear rate value : γ = 1.5 • 10 -3 . Simple shear deformation is obtained by the application of a succession of strain increments, δγ during times δt, to the cubic box of size L. At each time step, the positions of the particles are updated according to : x ′ = x + δγy, y ′ = y and z ′ = z. Then, the system is allowed to relax through viscous dissipation during time interval ∆t. The applied shear rate is defined as γ = δγ ∆t . Cyclic shear experiments are performed, during which the simulation box is strained up to a value γ = γ max . Then the shear is reversed, i.e the strain increment is -δγ and the system is deformed up to γ = -γ max . Finally the cycle is completed by the application of shear in the positive direction back to γ = 0.

To compare the numerical simulations with the experiments, we have to consider equivalent packing fractions. To this aim, we use predictions of the mode coupling theory to estimate the packing fraction ϕ c = 0.585 at which the breakdown of ergodicity occurs. We then determine the relevant packing fractions ϕ from the separation parameter ϵ = ϕ-ϕc ϕc . As a result, we numerically consider ϕ ∈ {0.607, 0.617, 0.627}, corresponding to the experimental volume fractions {0.60, 0.61, 0.62}. Throughout this work, we consider for each packing fraction 50 samples of size L = 60.

To characterize the anisotropy of the suspension, we consider the fabric tensor F which has been used successfully for instance in granular media and in silica or metallic glasses to track residual anisotropy during cyclic shearing [START_REF] Rountree | Plasticity-induced structural anisotropy of silica glass[END_REF].

Microstructure based rheological models [START_REF] Phan-Thien | Constitutive equation for concentrated suspensions in Newtonian liquids[END_REF] , that aim at deriving the flow properties of suspensions from their microstructure, use tensorial descriptors of the suspensions microstructure, such as the fabric tensor [START_REF] Ozenda | A new rate-independent tensorial model for suspensions of noncolloidal rigid particles in Newtonian fluids[END_REF]. This approach is well suited to describe the flow behavior of non Brownian suspensions and has also been used to describe the deformation properties of Brownian concentrated suspensions. In particular when the interactions between the colloidal are attractive, their structure is strongly anisotropic and the fabric tensor allows to monitor the structure of the network formed by the attractive colloids submitted to shear [START_REF] Jamali | Microstructural Rearrangements and their Rheological Implications in a Model Thixotropic Elastoviscoplastic Fluid[END_REF][START_REF] Boromand | Structural fingerprints of yielding mechanisms in attractive colloidal gels[END_REF], and to describe the structure of the contacts [START_REF] Olsen | Modeling tensorial conductivity of particle suspension networks[END_REF][START_REF] Olsen | Coupled dynamics of flow, microstructure, and conductivity in sheared suspensions[END_REF] or the structure of the stresses at contacts [START_REF] Martys | Stress propagation in a concentrated colloidal suspension under shear[END_REF] between colloids in concentrated suspensions.

Following [START_REF] Rountree | Plasticity-induced structural anisotropy of silica glass[END_REF], we define the local fabric tensor as :

F (i) = 1 N nn i,j∈n.n r ij r ij ⊗ r ij r ij (3)
where i,j∈n.n denotes the sum taken over the nearest neighbors of particle i and N nn is the

F = 1 N N i=1 F (i) (4)
The global fabric tensor F can be diagonalized and its eigenvalues λ 1 , λ 2 and λ 3 measure the degree of isotropy in the x, y, z directions. For a fully isotropic system one expects λ 1 = λ 2 = λ 3 = 1 3 . We define the scalar quantity α which quantifies the degree of anisotropy [START_REF] Rountree | Plasticity-induced structural anisotropy of silica glass[END_REF] :

α = 3 2 3 k=1 λ k - 1 3 2 (5)
For a fully isotropic (resp. anisotropic) system, α = 0 (resp. α = 1).

D. The Generalized Maxwell Model

We use the Generalized Maxwell Model, proposed by Voigtmann et al [START_REF] Voigtmann | Rheology of colloidal and metallic glass formers[END_REF][START_REF] Frahsa | On the Bauschinger effect in supercooled melts under shear: Results from mode coupling theory and molecular dynamics simulations[END_REF] in order to fit our experimental data. This model is based on the Integration Through Transients [START_REF] Fuchs | Theory of nonlinear rheology and yielding of dense colloidal suspensions[END_REF] description of the flow of glassy systems and allows for an analytic expression of the σ(γ) flow curve. It is first assumed that the stress may be expressed as an integral over the past applied shear rates, that involves a memory function, G(t, t ′ ) :

σ(t) = t -∞ dt ′ γ(t ′ )G(t, t ′ ) (6)
where G(t, t ′ ) depends on the applied deformation history. It is written as :

ν 0 1 - γ(t, t ′ ) γ * 2 e - γ(t,t ′ ) γ * * 2 e -(t-t ′ )| γ| γc (7)
The model thus involves four parameters. ν 0 is the elastic modulus at start-up flow when γ -→ 0, γ c /| γ| is the characteristic shear induced relaxation rate value, and γ * and γ * * capture the distribution of the shear induced relaxation rates over different length scales. In particular, it should be noted that this expression of the memory function allows G(t, t ′ ) to take negative values, depending on the relative values of γ * and γ * * .

This model has been shown to exhibit a well-defined stress overshoot, and yields a hysteresis behavior when the deformation is reversed.

III. STARTUP FLOW AND ELASTO-PLASTIC TRANSITION

Let us consider the start-up flow curve of suspensions of volume fractions 0.60, 0.61 and 0.62 under the application of a constant shear rate (Fig. 1 (a)). At low deformations, the response is linear and the elastic modulus is obtained and plotted in the inset of Fig. 1 (a)

. The elastic modulus G ′ at low deformation may be compared to the modulus obtained for suspensions of PMMA particles. We have, at ϕ = 0.60, G ′ 0 = 9.8 Pa = 11.6 k B T ⟨R⟩ 3 Pa, close to the value 7.5 k B T R 3 obtained for PMMA particules suspensions of identical concentration [START_REF] Amann | Overshoots in stress-strain curves: Colloid experiments and schematic mode coupling theory[END_REF]. Moreover, we have G ′ ∝ ϕ m with m = 14 (Fig. 1 (a) inset), in agreement with measurements of Schneider [START_REF] Schneider | On tuning microgel character and softness of cross-linked polystyrene particles[END_REF] for similar system, where m = 14 -17 was reported.

Then after yielding, characterized by stress overshoot, a plateau value of the stress is reached, corresponding to the flowing state. The numerical results obtained for the corresponding conditions are also given in Fig. 1 (b). The evolution of the maximum stress value, σ over both from experiments and numerical simulations are given in Fig. 2 (a). The amplitude of the stress overshoot increases when the shear rate increases : during a given amount of time, more stress is relaxed under a lower shear rate. The overall set of flow curves may be adjusted by the Generalized Maxwell Model. In particular the memory function exhibits negative values, that allow for the description of stress overshoot. The amplitude and the position of this negative value is controlled by the parameters γ * and γ * * . These parameters are slightly varied between adjustment curves at different shear rates and are respectively equal to 0.074 ± 0.013 and 0.092 ± 0.012.

The evolution of the stress overshoot amplitude as a function of the shear rate is well described by σ over ∝ γα with α = 0.14 for the experimental data and α = 0.16 for the simulations. This value is close to the exponent obtained for attractive carbopol gel particules where shear banding is observed [START_REF] Divoux | Stress overshoot in a simple yield stress fluid: An extensive study combining rheology and velocimetry[END_REF]. Although we do not have access to experimental measures of the local particle velocity, the numerical simulations do not exhibit any signature of shear banding, and the computed shear rate is constant through the gap of the cell.

One may also consider the amount of excess stress stored at yielding, relative to the stress corresponding to the maximal sheared state, under plastic flow, that is (σ over -σ ∞ )/σ ∞ .

Experimental and numerical results are given in Fig. 2 (b). This measures the ability of the system to store stress relative to the stress value of the suspension under a fully developed flow. At low shear rates,(σ over -σ ∞ )/σ ∞ increases with the shear rate, reflecting the fact that the suspension is able to store more stress, but at higher shear rates the suspension cannot deform more than in the fully developed flow and the stress excess saturates. Our data are consistent with the results obtained in [START_REF] Laurati | Transient dynamics in dense colloidal suspensions under shear: shear rate dependence[END_REF], where the evolution of the stress also exhibits a saturation of excess of stress at high shear rates.

IV. FLOW REVERSAL BEHAVIOR

Let us now consider the flow behavior after flow reversal. The suspension is first submitted to a constant shear rate γ = 5 • 10 -3 s -1 up to a deformation value γ max and then to a reverse flow of same shear rate. The evolution of the stress is given in Fig. 3, for ϕ = 0.62, both experimental results (a) and numerical results(b) are given. A strong hysteresis behavior is observed, already reported in athermal systems [START_REF] Fiocco | Oscillatory athermal quasistatic deformation of a model glass[END_REF]. The stress/strain relationship after flow reversal is markedly different from the startup flow curve. One gets a strong reduction of the stress overshoot. Adjustment of the reverse flow curves may be obtained from the GMM model using the parameters determined from the adjustment of the startup flow curve, and reproduce the reversal flow behavior. In order to understand this reduction, we consider the memory function G of the GMM model. As already discussed, during the startup flow, the memory function exhibits negative values responsible for the stress overshoot; upon flow reversal, the memory is erased by the e -γ/γc factor and the overshoot is reduced.

Experimental and numerical results exhibit a softening of the suspensions upon flow reversal (Fig. 4 (a) and (b)). We compute the effective shear modulus G back as the slope of the stress vs strain curve when the total stress gets back to zero (the slope is plotted in Fig. 3 (a)). This effective modulus upon shear reversal is smaller than the elastic modulus at rest and decreases when the maximum applied amplitude increases, reaching a plateau value at large γ max . This trend is well reproduced by the Generalized Maxwell Model.

Nevertheless, at the highest studied concentration, the amplitude of the decay is strongly overestimated by the model. More precisely, the GMM cannot describe quantitatively both the value of the overshoot and the amplitude of the reduction of the apparent elasticity after flow reversal and a slight overestimation is also visible for the other volume fractions. This indicates that the model does not capture the overall relaxation mechanisms that play a role in the non-linear response of the suspension, in particular at the largest volume fractions.

Another characteristic feature of the non-linear behavior of the response is its hysteretic behavior. The area enclosed by the hysteresis curves is represented in Fig. 4 (c) and(d).

At low maximal deformation amplitudes, the area enclosed by the hysteresis curves is close to 0, and the response is close to purely elastic. Then, with increasing value of γ max , the amplitude enclosed by the curves increases. Moreover the experimental evolution of the area as a function of γ max exhibits two behaviors, corresponding to γ max larger or smaller than the deformation at which the overshoot is observed, γ over : at large maximum deformation values, the enclosed area increases slower than at smaller deformations, which may be understood as the fact that the stress reaches a constant value at high γ max . This behavior is not recovered in simulations, where the slope of the area as a function of γ max does not decrease at γ over . This is due a slight difference in the width of the hysteretic response at small deformation amplitudes.

These results show that, even though flow non-linearities are qualitatively explained by a simple model of the relaxation mechanisms and may be reproduced by numerical simulations, the exact dissipation mechanisms require a more precise description of the interaction potential between the particles and a more refined relaxation spectra.

V. MICROSCOPIC INTERPRETATION

To understand the microscopic origin of the Bauschinger effect in the colloidal suspensions, we follow the anisotropy during the deformation and, in particular, we use numerical simulation to characterize the possible differences between states when no constraint is applied on the systems, i.e σ = 0. As we perform simple shear in the xy plane, we first look at the evolution of F xy , the xy component of the fabric tensors with strain. The hysteresis cycles of the stress and of the xy component of the fabric tensor, F xy are plotted in Fig. 5 (a), for a series of shear maximum amplitudes, γ max . F xy follows a similar hysteretic cycle as the stress : at low amplitudes, the anisotropy increases linearly with the deformation and takes the same values when the shear is reversed. When the amplitude increases, stress overshoot develops, and corresponds to an overshoot of the anisotropy of the microstructure. Nevertheless, the evolution of F xy as a function of the deformation is delayed, when compared to the evolution of the stress. In particular, after flow reversal, whereas the overall stress comes back to a null value, the microscopic structure still exhibits an anisotropy characterized by a positive value of F xy . The suspension has not relaxed back to its equilibrium isotropic structure and it might be expected that a different linear response is obtained.

By plotting the evolution of F xy with the stress (Fig. 5 (b)), one may quantify the delay of the microstructure anisotropy relative to the stress. At startup flow, during the initial stress increase, the stress linearly increases with F xy up to the stress overshoot (Fig. 5 (b) inset ). At the overshoot, the stress follows a sublinear increase with F xy , and, when the strain is reversed, it evolves according to an ellipse (Fig. 5 (b) main graph), which may be described as

σ = σ 0 cos ωt (8)
F xy = F 0 xy cos(ωt + φ) (9)
where ω = 2π 4 γ γmax , γ max being the maximum shear applied during a cycle. This defines φ, the phase between the stress and the xy component of the fabric tensor. The ellipse parameters are obtained by singular value decomposition of the data and the results are given in Fig. 5 (c) where φ is plotted as a function of the maximum applied deformation. It increases during the initial startup flow phase, until the deformation corresponding to the stress overshoot and then remains constant as a function of the maximum applied deformation. This corresponds to the observation of the decay of the effective modulus, that decays at low maximum shear amplitudes and reaches a plateau at γ ≈ 0.1. This indicates that the Bauschinger effect is associated to a phase lag between the stress and the deformation that builds up during the initial startup flow and remains constant when under the established plastic flow regime.

VI. CONCLUSION

Concentrated colloidal suspensions exhibit a rich flow behavior, and their mechanical properties depend on the previously applied stress. In particular, after flow reversal, the deformation is not reversible. This leads to a hysteresis curve in the σ(γ) plane whose area is the energy dissipated during a cycle. Moreover, the shape of the hysteresis curve does not possess a central symmetry relative to the non deformed state at zero stress, and shearing in an opposite direction leads to a different response from the deformation in the direction from a rest state : the stress overshoot disappears and the apparent elastic modulus is reduced after flow reversal. These observations are in agreement with numerical simulations [START_REF] Frahsa | On the Bauschinger effect in supercooled melts under shear: Results from mode coupling theory and molecular dynamics simulations[END_REF] and our own simulations show that the anisotropy of the suspension microstructure is retarded relative to the macroscopic stress after flow reversal.

FIG. 1 .

 1 FIG. 1. (a) Schematic illustration of the shear deformation protocol. A series of strain steps is applied to the sample, of constant amplitude δγ = 5 • 10 -3 . The step duration δt is varied between 5 • 10 -1 s and 20 s in order to control the effective shear rate between 10 -2 and 2.5 • 10 -4 s -1 . (b) Stress evolution between 10 -1 and 10 s after the application of a step strain of amplitude 5 • 10 -3 . The horizontal continuous line indicates the stress average over the time interval [5 • 10 -1 , 10] s. The dashed lines indicate the standard deviation of stress over this time interval. (c)-(d) Evolution of the stress as a function of the deformation, during the start-up flow. (c) Experimental data for ϕ = 0.62 and for different shear rate values, increasing from bottom. to top : γ = 2.5 • 10 -4 , 5 • 10 -4 , 10 -3 , 2.5 • 10 -3 , 5 • 10 -3 , 10 -2 s -1 . Continuous curves are adjustment with the Generalized Maxwell Model. d) Numerical data for ϕ = 0.627 and γ = 3 • 10 -4 , 6 • 10 -4 , 1.5 • 10 -3 , 3 • 10 -3 , 6 • 10 -3 , 1.5 • 10 -2 , 3 • 10 -2 , 6 • 10 -2 , 1.5 • 10 -1 and 3 • 10 -1 s -1 . Inset of (c) and (d) Evolution of the elastic modulus as a function of the volume fraction. The dashed lines are power-law adjustments, leading to G 0 ∝ ϕ m , with m = 14 for the experimental data and m = 16 for the numerical data.

FIG. 2 .

 2 FIG. 2. Evolution of the stress maximum at the overshoot (a) and of the overshoot stress relative to the stress at high shear rate (b) as function of the applied shear rate for ϕ = 0.62 (experimental data, •) and for ϕ = 0.627 (numerical data, ■). (a) The lines are power law adjustment of the data, leading to σ over ∝ γα with α = 0.14 (continuous line, experimental data) and 0.16 (dashed line, numerical data).

FIG. 3 .

 3 FIG. 3. Stress as a function of the deformation during the first cycle of deformation : the deformation is increased from 0 to γ max ; then decreased from γ max to -γ max and increased back from -γ max to 0. at a constant absolute value of the shear rate. (a) Experimental data for ϕ = 0.62 and γ max = 0.03, 0.05, 0.07, 0.09, 0.11, 0.13, 0.15, 0.17 and 0.20. The continuous black lines are adjustment of the data with the Generalized Maxwell Model using a single set of parameter values for all cycles. The thick segments are linear adjustment of the σ(γ) curves around zero stress, after flow reversal, from with the effective modulus, G back is computed. (b) Numerical data for ϕ = 0.627 and for the corresponding set of γ max values as experimental data.

FIG. 4 .

 4 FIG. 4. (a) and (b) Evolution of the effective modulus, measured as the slope of the σ vs γ curve when the stress is null, after flow reversal (see Fig. 3, for experimental systems (a) and numerical systems (b), as a function of the maximum deformation applied, γ max . In (a) dashed lines are the evolutions of G back obtained from the adjustment of the experimental data with the GeneralizedMaxwell Model (see also Fig.3). The horizontal thick lines correspond to the elastic modulus at rest, defined as the slope of the startup flow curves when γ -→ 0 (in (a), continuous lines : experimental data, dotted lines : GMM data). (c) and d) Evolution of the area of the hysteresis cycles as a function of their maximum amplitude γ max . Vertical dashed line indicates the γ value at the stress overshoot, γ over . Continuous lines are linear adjustments of the data for γ values larger and smaller than the deformation at stress overshoot. Three volume fractions are studied, increasing from black to grey : ϕ = 0.60, 0.61 and 0.62 for the experimental system and ϕ = 0.607, 0.617 and 0.627 for the numerical system.


Alvarez-Vanhard Emilien

Fernandez-Garcia Guglielmo

Corpetti Thomas

Super-resolution by fusing multi-spectral and terrain models: application to water level mapping

Keywords: Super resolution, fusion networks, Sentinel-2, Digital Terrain Models, flood dynamics, wet grasslands

Recently, deep convolutional networks have made great progress on the task of super resolution, i.e. reconstructing images with finer spatial resolution. However, although the reconstructions are visually impressive, they may lack physical consistency. This aspect is sought in remote sensing, where the resolution of satellite imagery (e.g. Sentinel-2) may be too coarse to characterize the physical structure and dynamics of certain landscapes. Through the study of flooding dynamics in wet grasslands, we propose a super resolution approach that allows deriving fine resolution patterns that are visually realistic and physically exploitable. This approach is based on an architecture, Fusion-UNet, allowing the fusion of multispectral data with a digital terrain model (DTM) associated with a loss function combining content, structure and segmentation losses. Our results show that this model can precisely predict water levels while restituting the fine structure of the landscape. This approach allows to refine the production of hydrological and ecological indicators to define the state of the ecosystem.

I. INTRODUCTION

Monitoring the water levels in wet grasslands is of prime importance to recover important ecosystem functions and to help managers to better manage their resources. Indeed, water dynamics of these areas lead ecological, hydrological, physical, biogeochemical, climatic and agronomic functions [START_REF] Rodríguez-Iturbe | Hydrologic dynamics and ecosystem structure[END_REF]- [START_REF] Damgaard | Plant colonization and survival along a hydrological gradient: demography and niche dynamics[END_REF]. However, due to their spatiotemporal variability, flood measurements (water extent, level and residence time) remain a challenge. To this end, Sentinel-2 (S2) data can be an appealing solution.

The arrival of S2 data has indeed caused significant changes in the monitoring of land use/land cover [START_REF] Phiri | Sentinel-2 data for land cover/use mapping: A review[END_REF]. Thanks to its good compromise between resolution on spatial (from 60 m to 10 m depending on spectral bands), temporal (one revisit every 5 ≈ 6 days) and spectral (12 bands) axes, many new use cases have been explored in a wide range of applications related either to agriculture monitoring [START_REF] Vuolo | How much does multi-temporal sentinel-2 data improve crop type classification?[END_REF], imperviousness estimation [START_REF] Lefebvre | Monitoring urban areas with sentinel-2a data: Application to the update of the copernicus high resolution layer imperviousness degree[END_REF], vegetation [START_REF] Wong | High-resolution calculation of the urban vegetation fraction in the pearl river delta from the sentinel-2 ndvi for urban climate model parameterization[END_REF], or water quality for example [START_REF] Toming | First experiences in mapping lake water quality parameters with sentinel-2 msi imagery[END_REF]. S2 has also revealed to be an appealing complementary source to other satellite (Sentinel-1 radar, higher spatial resolution as WorldView-3 satellite with 1.24m) for data fusion [START_REF] Wang | Fusion of sentinel-2 images[END_REF], [START_REF] Wang | Spatio-temporal fusion for daily sentinel-2 images[END_REF] as for example in [START_REF] Hajj | Synergic use of sentinel-1 and sentinel-2 images for operational soil moisture mapping at high spatial resolution over agricultural areas[END_REF] for soil moisture analysis.

Sentinel-2 features make such data widely used, in particular for the monitoring of dynamics, as floods, through the production of temporal series revealing fine patterns as seasonal and monthly intra-variability. However, in some applications, the 10 m spatial resolution may be limited. This is more generally true in the field of ecology where important variables generally belong to higher resolutions [START_REF] Skidmore | Priority list of biodiversity metrics to observe from space[END_REF] such as the biological effects of irregular flooding due to small scale topographic features (e.g. ditches). The particular example of flood dynamics monitoring in wet grasslands is relevant because these areas are composed of fine spatial patterns (below 10 m spatial resolution) that reflect habitat structure and ecosystem fragmentation [START_REF] Rodríguez-Iturbe | Hydrologic dynamics and ecosystem structure[END_REF], [START_REF] Krause | Historical and recent fragmentation of temperate floodplain grasslands: Do patch size and distance affect the richness of characteristic wet meadow plant species?[END_REF]. Moreover, other physical fine-scaled variables play important roles in the characterization of wetlands, such as water level or volume, that have to be evaluated together with the spatial patterns.

The recent progress of machine learning methods, especially with convolutional neural networks (CNN), associated with the availability of high-resolution data make now possible the reconstruction of high-resolution images using S2 data. This task, called super-resolution (SR), consists in learning an upscaling model to improve the spatial resolution of S2 data. In practice, from pairs of S2 images and corresponding high-resolution, the idea is to learn a model able to improve the initial S2 resolution and authorize an analysis of spatial patterns at a finer scale.

Today, some techniques give visually impressive results using adversarial losses [START_REF] Wang | Esrgan: Enhanced super-resolution generative adversarial networks[END_REF]. However, despite being visually plausible, super-resolved images remain not exploitable in terms of physical consistency. The use of Generative Adversarial Networks (GANs) indeed generates patterns which, although realistic, do not correspond to any reality on the field, preventing practical applications related, for example, to vegetation monitoring, urban growth analysis, estimation of flooded areas, etc [START_REF] Yates | Evaluation of synthetic aerial imagery using unconditional generative adversarial networks[END_REF].

In this paper, we tackle the issue of estimating high resolution water level (WL) images from S2 images. To this end, we exploit the fact that flooded areas are connected to local topology and we propose an original neural network that fuses S2 images with Digital Terrain Models (DTM) to estimate water level images with high spatial resolution. Furthermore, to eliminate the use of adversarial losses likely to generate realistic but physically inconsistent patterns, we suggest relying on structural/content losses to derive visually realistic and physically exploitable patterns.

II. RELATED WORKS

Based on the fact that the topography largely controls the spatial distribution of flooded areas, most of the superresolution of water level mapping are today performed by mixing a coarse resolution water surface information (e.g. satellite observations) with finer resolution data, usually issued from DTM. For example, the studies in [START_REF] Ling | Waterline mapping at the subpixel scale from remote sensing imagery with high-resolution digital elevation models[END_REF]- [START_REF] Zwenzner | Improved estimation of flood parameters by combining space based SAR data with very high resolution digital elevation data[END_REF] propose optimal fittings of the position of flood mask profiles according to the underlying terrain information to extract higher resolution flooded information. Such approaches are straightforward but suffer from some drawbacks : 1) sensitivity to classification errors, particularly due to the coarse resolution ; and 2) difficulty to detect a reliable boundary between flooded and non-flooded areas.

On the computer vision side, with the democratization of CNN architectures in the last decade, several techniques have emerged to design relevant reconstruction of high-resolution data from low-resolution ones (see for example [START_REF] Wang | A comprehensive review on deep learning based remote sensing image super-resolution methods[END_REF], [START_REF] Yue | Image super-resolution: The techniques, applications, and future[END_REF]) or to improve the resolution of some spectral bands [START_REF] Gargiulo | Fast super-resolution of 20 m sentinel-2 bands using convolutional neural networks[END_REF], [START_REF] Lanaras | Super-resolution of sentinel-2 images: Learning a globally applicable deep neural network[END_REF]. By definition, the super-resolution task is subject to two major problems: i) ill-posedness: several high-resolution images can give the same low-resolution one and ii) accuracy analysis: the definition of an objective criterion evaluating the quality of an estimated super-resolution image is tricky (usual image-comparison metrics do not discriminate details and high frequencies). Among the various CNN architectures, GANs have brought interesting solutions to these two issues [START_REF] Ledig | Photorealistic single image super-resolution using a generative adversarial network[END_REF]. In particular, they have partly solved the problem of finding an objective criterion to optimize: usual loss functions based on image comparisons are replaced by an adversarial task where a generative network is trained to generate fake data with realistic features and a discriminative network is trained to recognize whether the input is real or fake data. The first network tries to fool the second network and an increase of the error rate of the latter means that it becomes increasingly difficult to differentiate false data from real data, until the generated data is sufficiently realistic. GAN based super-resolution approaches are the current state of the art for SR tasks by directly learning an end-to-end mapping between low-resolution and high-resolution images.

Although for everyday life image analysis applications, super-resolution images produced by GAN-based networks are of excellent quality, when focusing on details related to the physics of the reconstructed objects, one observes errors or a lack of physical consistency. For example, in the case of water level mapping, GANs do not ensure that topographical features (e.g. ditches, riverbanks, depressions, etc.) and their geographical location are respected. This is due to the adversarial losses that are likely to introduce visually realistic images, but physical consistency is not ensured. This is the reason why, in addition to adversarial losses, some studies have proposed to use perceptual losses based for examples of features from Visual Geometry Group (VGG) neural network [START_REF] Romero | Super-resolution of sentinel-2 imagery using generative adversarial networks[END_REF], [START_REF] Johnson | Perceptual Losses for Real-Time Style Transfer and Super-Resolution[END_REF] or on contours [START_REF] Masoud | Delineation of agricultural field boundaries from sentinel-2 images using a novel super-resolution contour detector based on fully convolutional networks[END_REF].

In this study, we suggest to super-resolve S2 images to recover flooded areas by i) efficiently combining the DTM (to rely as much as possible on the physics of wet grasslands) to S2 and ii) exploiting various loss functions to reconstruct realistic patterns, without the use of GANs.

III. MATERIALS AND METHODS

A. Data Description

We address the super resolution's problem for water-level estimation in the Sougéal marsh (175 ha), part of the river basin of the Couesnon valley (France, 48.52 N, 1.53 W, Brittany region). The flooding dynamics of wet grasslands are characterised by high spatiotemporal variability, due to the heterogeneity of the topography, hydraulic conductivity of soils and hydrological discontinuities [START_REF] Hayashi | Hydrology of Prairie Wetlands: Understanding the Integrated Surface-Water and Groundwater Processes[END_REF]. Therefore marshes are usually divided in hydrogeomorphological units, i.e. portions of space whose hydrological response is homogeneous with respect to external stresses. This work used three piezometers for measuring WL in three different units through the study site.

Atmospherically corrected [START_REF] Hagolle | MAJA Algorithm Theoretical Basis Document[END_REF] cloud-free S2 images were downloaded for 74 dates between 01-01-2017 and 23-04-2021, chosen to take seasonal dynamics into account: in order to avoid an imbalance between flooded and non-flooded areas the acquisitions were mainly selected during the wet period. Only bands with ground resolution between 10 m and 20 m were used.

DTM with 5 cm ground resolution was generated from photogrammetric processing of Unmanned Aerial Vehicle (UAV) images [START_REF] Annis | UAV-DEMs for Small-Scale Flood Hazard Mapping[END_REF], [START_REF] -P. Schumann | Rapid Mapping of Small-Scale River-Floodplain Environments Using UAV SfM Supports Classical Theory[END_REF]. The flight was carried out during the summer of 2020, when the grazing is well advanced in order to have a low vegetation cover. High vegetation patches and features that stand out from the scene surface (e.g. cows) have been removed and elevation has been approximated by Inverse Distance Weighting (IDW) interpolation. The effect of vegetation is not fully corrected as vegetation after grazing is around 10 cm high which corresponds to the vertical error of the DTM : 9 cm (+/-6 cm), estimated via 50 positioning surveys by differential Global Navigation Satellite System with Post-Processed Kinematic (centimetric vertical precision). This vertical error identified in the production of ground truth data represents a decimeter bias in assessing the accuracy of predicted water levels. The DTM was resampled at 1 m to reduce the scale factor with the S2 data. We assume that topographic changes are negligible over the studied period.

The data has been, first, divided into tiles of 112x112 pixels, resulting in 4620 tiles. Then, this dataset was divided into three sets, with the train set comprising data from the years 2017-2019 (75% of total tiles), the validation set 2020 (23%) and the test set 2021 (2%). Finally, data mirroring and random deformation techniques were used for data augmentation.

B. Fusion network for super resolution

Our approach, called Fusion-UNet (Fig. 1), consists of merging two sources of information useful for estimating surface water heights: coarse-resolution optical data (S2) and fine-resolution topographic data (DTM). The data fusion is performed through a semantic segmentation model (U-Net [START_REF] Ronneberger | U-Net: Convolutional Networks for Biomedical Image Segmentation[END_REF]) adapted to super-resolution and regression.

Audebert et al. [START_REF] Audebert | Distance transform regression for spatially-aware deep semantic segmentation[END_REF] have proposed different fusion strategies that go beyond the simple concatenation of input data, as treating each data source in its own encoding branch. Inspired by this work, we applied an early fusion using a third "virtual" branch that allows fusion during encoding in a symmetrical way, i.e. by treating the two data sources independently. This approach, which encodes the two data sources separately and in parallel, avoids imbalances arising from the different nature of the data. The "virtual" branch corresponds to the sum at each level of the two encoding branches, which is then used for explicit connections to the decoding branch (and also at the lowest encoding level). In addition, a regularisation step with batch normalisation was added to the original blocks [START_REF] Ioffe | Batch Normalization: Accelerating Deep Network Training by Reducing Internal Covariate Shift[END_REF].

Due to the task at hand, that is pixel-wise regression from fused data, the structure of the loss function plays a critical role. Given an input y ∈ R d 1 ×d 2 ×c , with d 1 and d 2 its sizes and c the number of channels, and the predicted output ŷ ∈ R d 1 ×d 2 , the loss function L(ŷ, y) can be written as a sum of three parts:

L(ŷ, y) = λ p ℓ ϕ,j percep. (y, ŷ)+λ c ℓ cont. (y, ŷ)+λ s ℓ seg. (y, ŷ) (1)
where ℓ ϕ,j percep. (y, ŷ) is a perception function [START_REF] Johnson | Perceptual Losses for Real-Time Style Transfer and Super-Resolution[END_REF], ℓ cont. (y, ŷ) a content function and ℓ seg. (y, ŷ) a segmentation function. All three terms are preceded by a λ hyperparameter to adjust its importance with respect to the others. For the perception term, we used a pre-trained network (VGG-16 [START_REF] Simonyan | Very Deep Convolutional Networks for Large-Scale Image Recognition[END_REF]) as a feature extractor, with ϕ j the features at the jth convolutional layer :

ℓ ϕ,j percep. (y, ŷ) = 1 c j d 1 j d 2 j |ϕ j (ŷ) -ϕ j (y)| (2)
with c j ,d 1 j and d 2 j depending on the layer. This term encourages the network to generate similar features for ŷ and y. On the other hand, ℓ cont. (y, ŷ) is the simple L 1 distance between input y and output ŷ:

ℓ cont. (y, ŷ) = 1 d 1 d 2 |ŷ -y| (3)
Finally, ℓ seg. (y, ŷ) forces to consider non-wet surfaces as 0 m water heights, with reference data classified as "water" and "non-water". This allows us to frame this part as a semantic segmentation term and to use the IoU (Intersectionover-Union) as a loss function.

The hyperparameters have been chosen empirically and kept fixed for all experiments, except for the λ parameters: an SGD optimiser was used with a momentum of 0.9, gamma of 0.2, step of 15. Due to computing limitations, all experiments ran for 50 epochs with a batch of 10 images. The model algorithm is written with Pytorch and is available here : https://gitlab. com/EmilienAlvarez/spim-nn.

IV. EXPERIMENTS AND DISCUSSION

In this section, we aim to evaluate the contribution of the fusion strategy and the loss function for reconstructing highresolution images of flooded areas. For the first task, we compared Fusion-UNet with a standard UNet with two configuration: with (UNet-MNT-S2) and without (UNet-S2) DTM data as auxiliary data in input. To evaluate the performance of the models on N estimated data {ŷ 1 , ..., ŷN }, the mean absolute error (MAE), which compares the ŷ estimates with the y reference data, is often used :

M AE = 1 N × d 1 × d 2 N k d 1 i d 2 j |ŷ k (i, j) -y k (i, j)| (4)
where, i and j correspond to the pixel's indices. Given the strong disparity in the distribution of water levels in the reference data -non-flooded and lightly flooded areas are over-represented in the dataset -the use of the MAE, alone, may then appear biased. Therefore, we adjust this metric by rebalancing the contribution of WLs in its calculation. The sample is then divided into G ranges of 10 cm WL, i.e.:

M AE adj = 1 G G k=0 M AE k (5)
with k an index running on the total G ranges.

A. Contribution of the fusion strategy

The performance evaluation of each model is presented in Table I. The model UNet-S2 shows the worst performance with an M AE adj of 0.16 m on the test dataset. On the same dataset, the models UNet-MNT-S2 and Fusion-Unet both show better performances with an M AE adj of 0.09 m. However, on the training and validation datasets, Fusion-Unet slightly outperforms UNet-MNT-S2 by +0.01 m. Finally, whatever the model, the results on the validation and test dataset are better than on the training. This point accounts for an imbalance between the learning, validation and test sets due to the interannual variability of the flooding dynamics. Indeed, the high ranges of WL, i.e. >140 cm, are rare and only present in the training set.

Visually, we can see differences between the predictions of the different models (Figure 2). It can be seen that the predictions of the model UNet-S2 are correct for coarse spatial patterns but do not recover fine patterns. The predictions of the models UNet-MNT-S2 and Fusion-UNet, which incorporate the fine resolution DTM, are closer to the reference with very similar spatial patterns. However, the UNet-MNT-S2 model restores in its predictions the coarse resolution pixel raster of the Sentinel-2 data, unlike the Fusion-UNet model.

These results show, visually and via the evaluation metric (M AE adj ), that the Fusion-UNet model performs better than the UNet-S2 and UNet-MNT-S2 models. Therefore, the rest of the performance and sensitivity evaluation focuses on this model and validation dataset.

B. Contribution of the loss function

We first evaluated the contribution of the perception loss in the computation of the compound loss. We evaluated the perception loss computed for various positions of the VGG network from low to high abstraction levels. The first layers (low) indeed concentrate fine-scale structures, while deeper layers (high) extract more global scales. We observed that each layer used separately could not significantly reduce alone the M AE adj value contrary to the use of the mix of each layer (Table IV Then we analyzed the efficiency of content vs perception losses, i.e. ℓ cont. and ℓ percep. proposed in relations (3) and (2). To this end, we train a network with λ c = 1 -λ p and λ s = 0 for a set of values of λ p ∈ [0, 1] (Fig. 3). This experiment shows that the perception loss is better optimized with λ p tending to 1, i.e. higher weighting for this component in the compound loss, contrary to the content loss and the M AE adj metric which are optimal for around 0.5-0.6, i.e. with an equal mix. The segmentation loss shows no trend in the variation of the λ p coefficient and therefore does not seem to be influenced by the mixture of these two components. Finally, we evaluated the contribution of the segmentation loss ℓ seg. with λ c = λ p = 1 and a set of values of λ s ∈ [0, 1] (Fig. 4). We observed that this component in the compound loss has no significant influence on the M AE adj metric but enhances the distinction between water and no-water surface as the segmentation loss is minimized with λ s equal to 0.5, moreover it slightly reduces the perception and content losses.

V. CONCLUSION

The proposed approach exploits the super resolution based on CNNs to map water levels at fine resolution from S2 satellite data. This approach allowed to improve water level mapping in wet grasslands and to derive indicators of ecosystem status for better management.

The Fusion-UNet and UNet-MNT-S2 models make restoring physically exploitable patterns possible, but the former provided more visually realistic results. A sensitivity study on the different components of the compound loss function has allowed to improve the results and to show the interest in crossing content, perception and segmentation losses.

Future research should focus on exploiting SAR data as well, in order to avoid the cloud constraint and to increase the acquisition time series. Domain adaptation methods should be used to jointly process these heterogeneous data through a single model.

Fig. 1 .

 1 Fig. 1. Fusion-UNet workflow.

  [0]; (0, 10]; ... (G -10, G]. The adjusted MAE (M AE adj ) is then the average of the MAEs of each G ranges, i.e. :

Fig. 2 .

 2 Fig. 2. Visual comparison of the water level predictions of the models UNet-S2, UNet-MNT-S2 and Fusion-UNet. a) Sentinel-2 false colour composition (R: B8 -G:B4 -B:B3); b) UAV digital terrain model; c) reference water level; d) UNet prediction (S2); e) UNet prediction (S2 + DTM); and f) Fusion-UNet prediction

Fig. 3 .

 3 Fig. 3. Losses and metrics comparison for various weighting of the content and perception loss.

Fig. 4 .

 4 Fig. 4. Losses and metrics comparison for various weighting of the segmentation loss.

TABLE I PERFORMANCE

 I COMPARISON (M AE adj IN METER) BETWEEN THE MODELS UNET-S2, UNET-MNT-S2 AND FUSION-UNET.

	Models	Training	Validation	Test
	UNet-S2	0.27	0.24	0.16
	UNet-MNT-S2	0.17	0.14	0.09
	Fusion-UNet	0.16	0.13	0.09

TABLE II PERFORMANCE

 II COMPARISON (M AE adj IN METER) BETWEEN THE USED OF DIFFERENT VGG LAYERS FOR THE COMPUTATION OF THE PERCEPTION LOSS.

		Low	Mid	Mid-High	High	Mix
	M AE adj	0.13	0.13	0.13	0.13	0.11

ACKNOWLEDGEMENTS

This work has been partly supported by ANR Projet-ANR-18-CE23-0022 MULTISCALE.