Renaud Hostache

Xijun Lai

Jérôme Monnier
email: jerome.monnier@insa-toulouse.fr

Christian Puech

Assimilation of spatially distributed water levels into a shallow-water flood model. Part II: use of a remote sensing image of Mosel river

Keywords: Hydraulic modelling, roughness parameters, satellite SAR images, Digital Elevation Model, variational data assimilation, hydraulic coherence

With rapid flood extent mapping capabilities, Synthetic Aperture Radar (SAR) images of river inundation prove to be very relevant to operational flood management.

In this context, a recently developed method provides distributed water levels from SAR images. Furthermore, in view of improving numerical flood prediction, a variational data assimilation method (4D-var) using such distributed water level has been developed in Part I of this study. This method combines an optimal sense remote sensing data (distributed water levels extracted from spatial images) and a 2D shallow water model. In the present article (Part II of the study), we also derive water levels with a ± 40 cm average vertical uncertainty from a RADARSAT-1 im-

Introduction

Over the last decades, important efforts have been set up to better understand and manage floods, which is one of the most important natural hazard in the world. To identify flood prone areas and to manage flood events, satellite images could be very useful. Especially because of their all weather image acquisition capability, Synthetic Aperture Radar (SAR) sensors are relevant tools for the spatial characterization of floods. In this context, radar images of floods are nowadays mostly used for instantaneous flood extent extraction.

Nevertheless, as mentioned by [START_REF] Smith | Satellite remote sensing of river inundation area, stageand discharge: A review[END_REF], there is no doubt that earth observation images contain information that goes beyond simple flood extents. In this context, [START_REF] Schumann | High-resolution 3D flood information from radar imagery for flood hazard management[END_REF] shown that high spatial resolution (25 m) SAR spaceborne images allow the estimation of distributed water levels in floodplains with reasonable uncertainty by merging SAR derived flood extent limits with a high-resolution high-accuracy Digital Elevation Model (DEM). Furthermore, hydraulic modelling is of paramount importance in most flood forecasting and management systems. Due to huge stakes in flood management, the reliability of these flood inundation models is of primary concern.

To be reliable, hydraulic models have to be constrained by using various observed data sets. The model calibration consists in forcing the model to provide outputs as close as possible to observed data by searching optimal values of its parameters. Generally, the main parameters to be calibrated in a hydraulic model are the roughness parameters (Manning coefficients) since these are difficult to determine a priori. In an operational context, a "hand" calibration is often done through trial tests with the use of point observations, such as recorded hydrographs at stream gauges. Nevertheless these data are often insufficient to make the calibration reliable [START_REF] Aronica | Uncertainty and equifinality in calibrating distributed roughness coefficients in a flood propagation model with limited data[END_REF] as no reference data is available in the regions in-between those point measurements and a "hand" calibration is not accurate. As a matter of fact, taking into account complementary observations in an advanced calibration could make a better constraining of the model and in return a better reliability. Spatially distributed data of a flood event, e.g.: water levels, would certainly help to enhance calibration reliability, but are very difficult to gauge in the field, especially due to the risk of people being injured. A rather simple alternative way to obtain such distributed data is satellite imagery. Recent studies, [START_REF] Schumann | Deriving distributed roughness values from satellite radar data for flood inundation modelling[END_REF][START_REF] Hostache | Water level estimation and reduction of hydraulic model calibration uncertainties using satellite sar images of floods[END_REF] have shown that SAR derived water levels allow the identification of Manning friction coefficients of one dimensional hydraulic models with a GLUE approach [START_REF] Beven | The future of distributed models -model calibration and uncertainty prediction[END_REF]. Furthermore, the recent improvements of mathematical tools provide more and more efficient means to calibrate numerical models. The variational data assimilation method based on the optimal control theory of partial differential equations (also called 4D-var method) offers a special powerful tool to fuse in an optimal sense measurements (observations) and the mathematical model. In river hydraulics, variational data assimilation methods have been used successfully for shallow water models, see e.g. [START_REF] Ding | Identification of Manning's roughness coefficients in shallow water flows[END_REF][START_REF] Belanger | Data assimilation (4D-VAR) to forecast flood in shallow-waters with sediment erosion[END_REF][START_REF] Roux | Parameter identification using optimization techniques in open-channel inverse problems[END_REF][START_REF] Castaings | Automatic differenciation: a tool for variational data assimilation and adjoint sensitivity analysis for flood modeling[END_REF][START_REF] Honnorat | Dassflow v1.0: a variational data assimilation software for river flows[END_REF][START_REF] Gejadze | On a 2d zoom for 1d shallow-water model: coupling and data assimilation[END_REF][START_REF] Honnorat | Lagrangian data assimilation for river hydraulics simulations[END_REF][START_REF] Marin | Superposition of local zoom model and simultaneous calibration for 1d-2d shallow-water flows[END_REF][START_REF] Honnorat | Identification of equivalent topography in an open channel flow using lagrangian data assimilation[END_REF].

The present study aims at investigating whether the variational assimilation of water levels derived from a flood SAR image could help to enhance the calibration of flood inundation models. This would help getting benefits of both variational assimilation and recently developed remote sensing methods, improving flood model calibration.

Based on the study of [START_REF] Raclot | Remote sensing of water levels on floodplains: a spatial approach guided by hydraulic functioning[END_REF], which provides water levels having an average uncertainty of ± 18 cm by using aerial photographs, the water level estimation method employed here is composed of three main steps [START_REF] Hostache | Estimation de niveaux d'eau en plaine inondée à partir d'images satellites radar et de données topographiques fines[END_REF]: i) extraction of the flood extent limits that are relevant to water level estimation, ii) primary estimation of water levels by merging the relevant limits and a high-resolution and high-accuracy Digital Elevation Model (DEM), iii) minimization of water level estimation uncertainties by hydraulic coherence constraints. This method is preferred in this study to the one developed by [START_REF] Schumann | High-resolution 3D flood information from radar imagery for flood hazard management[END_REF] since it can evaluate the observation uncertainty and does not need a one dimensional flow hypothesis.

Regarding the variational data assimilation, we apply the method detailled to Part I, see [START_REF] Lai | Assimilation of spatial distributed water levels into a shallow-water flood model. part i: mathematical method and test case[END_REF].

The paper is organized as follows.

In Section 2, we give a flood description and present available data, namely ground recorded hydrometric and remote sensing data. In Section 3, we present the methodology used to derive spatially distributed water levels from a SAR image. In Section 4, we recall very briefly the full mathematical and numerical model, including the variational data assimilation process. In Section 5, we show the interest of assimilating SAR derived water levels in order to calibrate the Manning coefficients. This is done in two steps. First, we perform a twin experiment (i.e. observations are generated from direct model). This aims at evaluating the a priori potential efficiency of the assimilation process. Second, the identification of the hydraulic model friction parameters (Manning coefficients) is addressed by assimilating in situ and SAR derived observation. This is done by using various (and classical) land use based spatial distributions of the Manning coefficients. In Section 6, we perform the sensitivity analysis tests in order to evaluate if such a priori land-use based spatial distributions of the friction parameters leads to the most accurate numerical model.

2 Flood description and data extraction.

Study area

The area of interest includes a 28 km reach of the Mosel River between Uckange (France) and Perl (Germany) -see Fig. 1. In this area, the Mosel River meanders in a flat plain having an average width of 3 km and a mean slope of 0.05 %. Most of the villages in the plain are located at an elevation slightly higher than the area covered by the February 1997 flood event. It is worth noting the presence of a narrow valley between Berg/Mosel and Perl cities (Fig. 1). The latter behave as a bottleneck during flood events causing upstream water retention area. In Perl city (downstream boundary), the Mosel basin has a size around 35500 km 2 . The propagation velocity of the flood peak in the study area is low, around 2 km.h -1 . The peak discharge recorded at Uckange city stream gauge (upstream boundary) was around 1450 m 3 s -1 (see Fig. 2), corresponding to a 4-5 year time return period [START_REF]Banque hydro[END_REF]. Because of its representativeness of river inundations in floodplains, this area has already been investigated within the framework of the NOAH European project [START_REF] Noah | New oportunities of altimetry in hydrology[END_REF] and the PACTES French project [START_REF] Puech | Flood forecasting and prevention using remote sensing data (PACTES) final report of cemagref contribution[END_REF] which aimed at analyzing the potential of remote sensing data -especially SPOT HRV images -for hydrological and hydraulic flood modelling. This shows the interest in analyzing this kind of river for generalizable methods.

Ground recorded hydrometric data

As hydrometric data, discharge hydrographs were available at three stream gauges located in the study area (see Fig. 1). These hydrographs are shown in Fig. 2 River during this flood event). Indeed, the water volume and the peak discharge recorded at Perl (downstream gauge station) during the flood events are higher (around 10 %) than those recorded in Uckange (upstream) and EDF gauge stations (middle), although there is no lateral inflow in the study area and the variation of the basin area is negligible. As a matter of fact, these recordings are in opposition to the hydraulic principles of mass conservation and peak discharge attenuation without inflow. Recorded discharge hydrographs araise from calculation using observed water stage hydrographs and rating curves: relationship between discharge and water level computed using discharge in situ measurements. Considering that higher magnitude discharges have been in situ measured in Uckange than in Perl for the the rating curve computation, the Uckange hydrograph has been assumed more reliable.

As a consequence, only the discharge hydrographs in Uckange and EDF stream gauges are used as ground truth information in this study.

Consequently, available ground observations are fairly limited. The time series data of discharges at the EDF gauge station are only available at the beginning and the end of the flood period; the measurements are lacking during the high flood stage because of sensor disability (see Fig. 3).

Remote sensing data

The image used in this study has been acquired at 6:00 AM, during the February 28th 1997 flood event, by the Synthetic Aperture Radar (SAR) sensor of the RADARSAT-1 satellite (descending orbit, standard beam S2, 5.6 cm wavelength, C band, Horizontal-Horizontal polarization). On C band radar images, the flood extent delineation is rather straightforward because the backscatter of smooth open water is very low due to specular reflection. However, wind may induce wavelets on water surface that may roughen it and increase its backscatter. This may occur for the SAR sensor of RADARSAT-1 especially if these wavelets have vertical dimensions higher than 2.8 cm, according to Rayleigh's criterion for a 5.6 cm C Band [START_REF] Ulaby | Handbook of radar scattering statistics for terrain[END_REF]. At the image acquisition time, the wind speed was moderate (7 m.s -1 in Metz city airport, French Meteorological Center record), thus the wind effects on open water surface roughness have been assumed negligible with respect to flooded area detection.

The SAR image, amplitude coded, has a pixel spacing of 12.5 m, resulting from the sampling of a complex image of 25 m spatial resolution. It was acquired a few hours after the flow peak, at the beginning of the recession, as shown in Fig. 2. Furthermore, five air photographs at 1:15 000 scale, flown by the French National Institute of Geography (IGN) in 1999 out of flooding period, have been acquired. They have been used to digitize land use information : buildings, urban areas and sparse habitat, and high vegetation, forests, sparse trees, hedges etc. Indeed trees and buidings perturb the backscattering signal, so we choiced to mask these areas before processing : the whole process runs using only reliable areas.

The topographic and bathymetric raw data have been provided respectively by the North-Eastern French Navigation Services (SNNE) and the DIREN Lorraine as 3D points and 3D lines -calculated by photogrammetry using air photographs at 1:8 000 scale and SONAR sounding. The average altimetric uncertainty on the raw data is about 25 cm (DIREN information) for the floodplain and 1 cm (SNNE information) for the channel. Using a linear interpolation between points and lines, a Triangular Irregular Network (TIN) Digital Elevation Model (DEM) has been generated and then converted to a RASTER DEM, easily superimposable with image data [START_REF] Hostache | Analyse d'images satellitaires d'inondations pour la car-actérisation tridimensionnelle de l'aléa et l'aide à la modélisation hydraulique[END_REF]. Hereafter, DEM acronym will be assigned to the topographic RASTER data. These interpolation choices are sufficient here, because main errors are not in the interpolation process : indeed, the initial data set (xyz points) is very dense, and sufficient to derive a good raster by linear interpolation ; moreover steep banks, potentially badly represented, are eliminated from the final process.

Water level estimation using SAR image

Based on the method developed by [START_REF] Raclot | Remote sensing of water levels on floodplains: a spatial approach guided by hydraulic functioning[END_REF] providing water level estimates with a ± 18 cm mean uncertainty using flood aerial photographs, the water level estimation method used in this study is detailed in [START_REF] Hostache | Estimation de niveaux d'eau en plaine inondée à partir d'images satellites radar et de données topographiques fines[END_REF]. The current section presents the main steps and the general philosophy of this method. Using the RADARSAT image, this approach provides spatially distributed water level estimates within a ± 40 cm mean uncertainty. It can seem strange to reach such accuracy with large pixels (25m). Nevertheless, more than to pixel size, the accuracy in this process relates to DTM accuracy and the fact that the vertical estimates concern only places with very low slopes and no perturbing items (trees, houses ..).

Flood extent extraction

Flood extent mapping using SAR images is now a quite common issue [START_REF] Smith | Satellite remote sensing of river inundation area, stageand discharge: A review[END_REF] because water appears with very low backscatter compared to other objects.

In this study, to discriminate water from no-water on SAR images, radiometric thresholding has been used because it is a robust and reliable way [START_REF] Henry | Systèmes d'information spatiaux pour la gestion du risque d'inondation de plaine[END_REF]. [START_REF] Hostache | Estimation de niveaux d'eau en plaine inondée à partir d'images satellites radar et de données topographiques fines[END_REF] shown that the use of a unique threshold value does not allow a satisfying flood extent extraction since flooded areas may sometimes have the same backscattering values as non-flooded one in a SAR image. To deal with this radiometric uncertainty, the discrimination of flooded and non-flooded pixels is done using two threshold values. The first threshold value, T min , aims at detecting only pixels that correspond to water bodies. T min has been determined as the minimum radiometric value of non-flooded pixels inside grassland areas (outside the floodplain and outside the permanent water surfaces) [START_REF] Hostache | Estimation de niveaux d'eau en plaine inondée à partir d'images satellites radar et de données topographiques fines[END_REF].

The second threshold value, T max , aims at detecting all flooded areas, at the risk of detecting in addition non-flooded areas that have a similar radiometric value to the flooded one. T max has been determined as the maximum radiometric value of water bodies outside the flooded area, using the SAR image pixels located inside the Mirgenbach lake (see Fig. 1). The thresholding of the SAR image using T min and T max provide a flood extent map with fuzzy limits, coded as follows (see [START_REF] Hostache | Estimation de niveaux d'eau en plaine inondée à partir d'images satellites radar et de données topographiques fines[END_REF]), depending on the intensity I of the SAR image pixels: 0 = non-flooded (I > T max ), 1 = flooded (I < T min ), 2 = fuzzy limit

(≈ potentially flooded) (T min ≤ I ≤ T max ).
The innovative point of the SAR image processing is the analysis of the relevance of the remote sensing-derived flood extent limits for hydraulic purpose, and especially for water level estimation. To estimate water levels, the flood extent limits are merged with the underlying DEM. Using such a merging, any erroneous flood extent limit will lead to errors in water level estimation.

Consequently, flood extent limits prone to error have to be identified before the merging. Errors in the flood extent limits are mainly due to emerging objects such as building and high vegetation [START_REF] Hostache | Estimation de niveaux d'eau en plaine inondée à partir d'images satellites radar et de données topographiques fines[END_REF] that may mask water. To treat this potential errors, it has been chosen to remove all SAR derived flood extent limits located in habitat or vegetation areas. The remaining limits will be called relevant limits hereafter. Furthermore, since the flood extent limits are SAR image derived, they are prone to uncertainties arising mainly from the image spatial resolution and its georeferencing precision. To take account of these spatial uncertainties, the relevant limits have been enlarged (on both sides) using a distance equal to the sum of the image spatial resolution (25 m for the RADARSAT image) and the precision of the georeferencing (10 m for the RADARSAT image), as proposed in [START_REF] Hostache | Estimation de niveaux d'eau en plaine inondée à partir d'images satellites radar et de données topographiques fines[END_REF]. Considering that errors have been removed and that radiometric and spatial uncertainties have been taken into account previously, these "enlarged" relevant limits are then assumed to include the real flood extent limits (Hyp. 1 ). These are shaped as small patches which are sparsely distributed along the floodplain (Fig. 4).

Preliminary water level estimation

The second part of the process estimates one range of possible water levels IW LE =[ WL min ; WL max ] for each relevant patch. To do so, the maximum and the minimum elevation values are first extracted inside each relevant patch using the DEM Z values. Next, the DEM altimetric uncertainty (uncert DEM ) is taken into account by being respectively added/subtracted to the maximum/minimum values extracted previously for each relevant patch:

IW LE =[min (Z patch ) -uncert DEM ; max (Z patch )+uncert DEM ].
Since the DEM altimetric uncertainty is taken into account, Hyp. 1 allows to assume that each range of water level estimation -IW LE -includes the real water level (Hyp. 2 ).

Final water level estimation

The last part of the process uses hydraulic rules to constrain the water level estimates and reduce their uncertainty, as proposed firstly by [START_REF] Puech | Using geographical information system and aerial photographs to determine water levels during floods[END_REF][START_REF] Raclot | Remote sensing of water levels on floodplains: a spatial approach guided by hydraulic functioning[END_REF]. In a floodplain, hydraulic laws manage the flow, so that water levels must follow an hydraulic logic : hydraulic energy decreases from upstream to downstream.

With low flow velocity, like in the Mosel floodplain, this hydraulic rule can be simplified into a decrease of water level in the flow direction (Hyp. 3 ). To apply Hyp. 3 on the IW LE intervals, flow directions between patches (locations of the water levels) have to be determined. As proposed in [START_REF] Hostache | Estimation de niveaux d'eau en plaine inondée à partir d'images satellites radar et de données topographiques fines[END_REF], some flow directions between patches have been determined using the shape of the SAR derived flood extent and the lines perpendicular to the elevation contour lines, oriented from the highest to the lowest elevation, called steeper lines hereafter. As a matter of fact, the knowledge of a up-/downstream relationship between two patches is conditioned by the connectivity of those on the SAR derived flood extent. If no hole in the flood extent is present between these two patches, an up-/downstream relationship can be determine. Furthermore, since the image has been acquired at the beginning of the flood recession, the main flow directions are assumed to be convergent toward the river channel, following the steeper lines. The steeper lines around relevant patches have been determined using the contour lines derived from the DEM.

As a matter of fact, the sense of an up-/downstream relationship is given by the orientation of the steeper line. Using this criterion, some up/-downstream relationship between patches have been determined. At the floodplain scale, these relationships constitute an hydraulic hierarchy of the relevant patches.

Consequently, according to Hyp. 3, the water level must decrease from the patch A to the patch B if A is upstream of B. Due to Hyp. 2 this induce the following constraints:

N max (B) ≤ N max (A) (constraint on the maxima), N min (A) ≥ N min (B) (constraint on the minima).
To apply these constraints, the algorithm that has been developed, is flow oriented and impose a decrease on the maxima from upstream to downstream, and vice et versa, an increase on the minima from downstream to upstream. This finally provides intervals of constrained water level estimation IW LE =[WL min ; WL max ], with a half mean range of about ± 40 cm.

As a consequence, the method allow the definition of a set of distributed water levels across the floodplain at the satellite overpass time. These remote sensing-derived water levels are then used as new observation for the hydraulic modelling.

As a matter of fact, remote sensing derived data are assumed to be an observation of water depth along the floodplain. The single image is extracted from the satellite image at 6:00 AM,Feb 28th,1997; The locations where the water level values have been extracted are shown in Fig. 4.

Mathematical model with variational data assimilation (4D-var)

This section describes briefly the mathematical model with the variational data assimilation (4D-var) method. For more details, we refer to Part I of the present study, [START_REF] Lai | Assimilation of spatial distributed water levels into a shallow-water flood model. part i: mathematical method and test case[END_REF], where the method has been validated on a test case.

Numerical computations are performed by using the DassFlow software [START_REF] Honnorat | Dassflow v1.0: a variational data assimilation software for river flows[END_REF].

Mathematical model

The flood flow is modeled by the two-dimensional shallow water equations (2D-SWEs):

∂U ∂t + ∂F(U) ∂x + ∂G(U) ∂y = B(U) (1)
where x,y are the coordinates and t is the time, U is the state vector, F and G are the x-and y-directional flux vector respectively, B is the source term vector. These vectors are defined as follows:

U =(h, hu, hv) T =(h, q x ,q y ) T (2)
F =(hu, hu 2 + 1 2 gh 2 , huv) T (3)
G =(hv, huv, hv 2 + 1 2 gh 2 ) T (4)
B = (0, gh(S 0x -S fx ), gh(S 0y -S fy )) T

where h is the water depth; u and v are the x-and y-directional velocity components; q x = hu and q y = hv are the unit discharge in the x-and ydirections. S 0x and S fx are the bed and friction slopes, respectively along the x axis and similarly S 0y and S fy for the y axis. The friction slopes are evaluated using the Manning formula.

The full inverse model (with variational data assimilation) includes the forward model (2D SWEs), the adjoint model and the minimization algorithm. The 2D-SWEs is solved using a finite volume method on unstructured meshes.

The adjoint model is directly derived from the source codes of forward model by automatic differention [START_REF] Hascoet | Tapenade 2.1 user's guide[END_REF]. The cost function J that has to be minimized is defined in detail later (see [START_REF] Ding | Identification of Manning's roughness coefficients in shallow water flows[END_REF] in the next section). As a matter of fact, the optimization problem to be solved is: min p J(p), where p is the control variable. The latter can be the initial condition, the Manning coefficient and/or inflow discharge.

Preliminary forward run

As a preliminary study, the Manning coefficients are set to an a priori "reason- The computational time step is set to two seconds (this respects the stability condition of the numerical scheme). The Manning coefficients are set to 0.025 all over the computational domain (this value is classically used to model the present flow).

The observed and computed water stage hydrographs in the middle gauge station (EDF) are plotted in Fig. 5. Table 1 presents some basic statistics of the differences between the computed and the measured water levels in the floodplain (extracted from the image).

The computed water level at EDF gauge station (Fig. 5) is very close to measurements; while in floodplain at image time, the computed water level fits much less with available data, see Table 1. Thus, this preliminary forward run shows that the numerical model can reproduce the flow; and the current solution will be in next section an acceptable hand made first guess.

Cost function

As described in Part I of the present study, [START_REF] Lai | Assimilation of spatial distributed water levels into a shallow-water flood model. part i: mathematical method and test case[END_REF], the cost function J to be minimized (see section 4.1) contains three terms:

J obs (p)=J obs + αJ f lux + βJ reg (6)
In ( 6) α and β are weight coefficients to be set and :

(1) J obs corresponds to the discrepancy between the observations and the computed flow state,

(2) J f lux corresponds to the discrepancy of net mass flux,

(3) J reg is a regularization term for smoothing time-dependent control variables p.

The weight coefficients α, β are setted by hand "at best" in order to equilibrate the three terms after convergence. The cost function J depends on the control variable p since the state variable U depends on p.

The term J f lux computes a "mix" net mass flux since it is based on the observed elevation h obs and the computed velocity (u, v) T , see Part I [START_REF] Lai | Assimilation of spatial distributed water levels into a shallow-water flood model. part i: mathematical method and test case[END_REF]. This term enhances significantly the assimilation of the spatial distributed observations in the floodplain (see Part I).

The term (J obs + αJ f lux ) is defined as follows:

J obs + αJ f lux = 1 2σ 2 z m z T imag (x m ,y m ) -z obs T imag (x m ,y m ) 2 [1 + α (u, v)(x m ,y m ) 2 ] + γ 1 2σ 2 Q [ T 1 0 Q(t) -Q obs (t) 2 dt + T T -T 2 Q(t) -Q obs (t) 2 dt] (7)
where σ z and σ Q are the standard deviations of observations, T imag is the image time, (x m ,y m ) is the position of water level measurements, γ 1 is a weighting coefficient, T 1 and T 2 are the time periods of the discharge measurement available at EDF station (see Fig. 3), T is the assimilation time period. The regularization term is:

J reg (p)= 1 2 ∂Q in ∂t 2 (8)
5 Calibration of spatially distributed Manning friction coefficients (land-cover based spatial distribution)

The Manning friction coefficients n, that represents the resistance to the flow in channels and floodplains, are empirical and most of the time cannot be measured directly. It can be locally very fluctuant due to vegetation growing or erosion deposition. In fact we must point out that (n) does not really exists, since it is strongly scale dependent, because it integrates all the friction processes at all scale. This section aims at showing the capability of SAR derived spatially distributed water levels to enhance the Manning coefficient calibration, in comparison to a "hand" calibration using trial-error tests. Thus, the data assimilation method presented here uses "orthogonal data": water heights derived from water limits versus classical water depths in gauge station, see Fig. 3.

In this section, the spatial distribution of Manning friction coefficient is based on land-cover classes as classically done. Various cases have been considered

depending on the total number of classes taken into account (between 1 and 10).

In the first approach, in order to investigate the reliability of our variational data assimilation process, while being free of measurement errors, twin experiments are performed: "observations" are first generated using the numerical model (Manning coefficient being set a priori) and our variational data assimilation process is next used to retrieve the Manning coefficients used to generate the observed data.

In the second approach, the data assimilation is performed with the real observed data described in Section 2.

In all numerical experiments, the assimilation period is 66h from 12:00, Feb 25th (flood event starting time) to 6:00 AM, Feb 28th,1997 (SAR image acquisition time).

From a mathematical point of view the deal is different : here we want to verify that the method is able to retreive hidden values. Within the set of parameters, when (n) is the hidden value the method permits to retreive the initial value, which is in this case unique and well defined.

Land-cover classification

The domain is decomposed by using land-cover classes, one Manning coeffi- We are expecting close n values within banks, and strong differences with n main channel value).

cient
c) five classes consisting of the main channel, the bridges (same n value for all bridges), the small lakes, and the left and right floodplains.

d) ten classes according to the land-cover classification presented in Fig. 6: main channel, bridges, small lakes, the left and right floodplains, vegetation, gravel, grassland, urban areas, downstream of the main channel.

Identification of Manning coefficients using synthetic data (twin experiments)

The identification process in the first approach is based on synthetic data. This means that given some Manning coefficients values, observations are generated by running the model. Then, the aim is to retrieve these Manning coefficients which are used to generate observations. Such experiments are called as twin experiments. They can test the reliability of the identification process, while being free of measurement errors, and then to evaluate the potential efficiency of the variational data assimilation process for identifying correctly the Manning coefficients.

It is worth noting that the synthetic data are generated exactly at the same locations as those of the real observations. The difference is that they is generated by the exact forward model.

For the sake of clarity, we consider here case b) only (see above) i.e. three land-cover classes. Fig. 7 shows the decreases of the normal cost function and its normal gradient vs iterations (convergence of the optimization process).

Table 2 presents the results obtained when considering three land-cover classes (case b). We start from one initial guess value (i.e. a single value for the three classes), and we compute the three optimal values corresponding to the three classes. Table 2 shows that the variational data assimilation process is capable of retrieving the three initial Manning coefficient values ("true" values): in the main channel, in the right and left floodplains.

Let us point out that various numerical experiments have been performed to corroborate this result. Indeed, we started the optimization process with various initial guesses (e.g. n 0 =0.020, 0.033 , 0.050), and all of them converged to the same optimal values, equal to the guesses.

The present identification experiment based on synthetic data, shows that the SAR derived spatially distributed water levels and the stream gauge in-situ measurements can identify properly the spatial distributed Manning coefficients, at least if three land-cover classes only are considered and if the math-ematical model is perfect (i.e. observations corresponds to the modelled flow) and no error is introduced in measurements. These last two features are very important. In next section, we use the real data and we will notice that the (real) identification problem is much more difficult to solve numerically, even with three land-cover classes only.

Identification of Manning coefficients using real data

In the second approach, the same identification process as in the previous section has been performed, but under the real observations. The upstream and downstream boundary conditions are respectively the observed discharge hydrograph in Uckange and the observed water stage hydrograph in Perl (see section 2). The initial condition is calculated as a steady state simulation of the model, using as boundary conditions the discharge and the water level observed respectively in Uckange and Perl at the simulation starting time. Then, given initial guesses for Manning coefficient values (see n 0 values in Table 3), the optimal values, i.e. minimizing the cost function J, are sought. For each case (lumped n, three, five or ten land-cover classes), the same minimum value of J is obtained.

Table 3 presents the initial and optimal computed n values for each case (depending on the land-cover number). In all cases, the computed value of the main channel is the same whatever the initial guess n 0 (e.g. 0.033, 0.020, 0.025 or 0.050) and the total class number. On the contrary, n in other land-cover classes are almost not modified during the optimization process. This result is consistent with the forthcoming sensitivity analysis (see next section), runs which shows that the total cost function is almost insensitive to the n values outside the main channel.

In twin experiments performed previously, data were perfect since they were created by the model. In that case, the cost function was much less insensitive to the two n values outside the main channel. That is why numerical identification was accurate and robust, Table 2.

Using real data, sensitivities can become different, and the identification problems become much more difficult. In the present case, the few sensitivity of the model to the Manning coefficient outside the channel is mainly due to the rather low time return period of the investigated flood event (4-5 year).

Indeed, for the 1997 flood event the most important part of the flow occurred in the channel. As a result, since the Manning equation is used to compute the friction slope, even a drastic change in the floodplain Manning coefficient value implies a low variation of water level when the discharge is low in the floodplain. This consideration shows that the water levels in the floodplain is more influenced by the water level inside the channel (and thus by the channel Manning coefficient) instead of the floodplain Manning coefficient.

In the case of a lumped n, we present in Table 1 some statistics on the difference between the simulated and the SAR derived water levels (where image information is available). Compared to the results obtained with the preliminary forward run, the present results show a significant improvement since the error indicators are drastically lowered. In Fig. 8 and Fig. 9, we plot the difference between the computed and the observed water level where image information is available and when considering one, three and ten land-cover classes. Fig. 10 shows the comparison between simulated and observed discharge hydrographs at the EDF middle gauge station (in function of time).

These results show that after the identification of an "optimal" set of Manning coefficient values, the computed flow state is much closer to the observations than that computed using the forward hand-calibrated model (i.e. with n values resulting of many trial-error tests).

Sensitivity analysis and Manning decomposition

The previous section shows that our variational data assimilation process is capable of enhancing significantly the calibration of the distributed Manning friction coefficients, and consequently the accuracy of the numerical model.

Nevertheless, all computations have been done with an a priori: the Manning coefficients have been defined spatially distributed according to given land-use classes (between 1 and 10). This a priori is rather traditional. Nevertheless, does this necessary lead to the most accurate numerical model ?

In the current section, we aim at answering to this question. To meet this aim, a sensitivity analysis is performed without any a priori for the Manning spatial distribution (i.e. any land-use is defined).

One run of the forward model plus the adjoint model, without running the full optimization process, gives the gradient of cost function J with respect to the Manning coefficients. These gradient values represent the local sensitivity of the cost function with respect to the Manning coefficients. They help one to understand which Manning area is the most important to calibrate.

The ten land-use case

First, we perform sensitivity analysis in the case the Manning is decomposed by the ten land-uses defined previously, see Fig. 6. We show the ten spatially distributed gradient values in Table 4 andFig This sensitivity analysis suggests that the most important Manning value to focus on is in the main channel, than much less important are the vegetation area, the bridge (in main channel), gravel area and grassland area. Others land-use values are negligible.

This sensitivity analysis result is consistent with the calibration process presented in previous section: the optimization algorithm calibrate essentially the main channel value.

No a priori land-use case

Next, a sensitivity analysis with no a priori on the Manning decomposition has been performed: we define one Manning value for each mesh cell, and we compute the gradient value for each cell. In particular in the main channel, there are as many potential values as cells. The computed spatially distributed gradient values are shown in Fig. 12.

The present sensitivity analysis suggests defining few Manning areas inside the main channel (at least four in the present case). More generally, this suggests to define the Manning areas not upon the land-use (or land-cover) but upon "sensitive areas". Furthermore, the combination of the satellite image and runs of sensitivity analysis indicate the most important Manning areas we need to calibrate (including inside the main channel).

Conclusions

Part I of this study [START_REF] Lai | Assimilation of spatial distributed water levels into a shallow-water flood model. part i: mathematical method and test case[END_REF], presents a method to assimilate spatially distributed water levels into a shallow water model. The present article (Part II) deals with an application of the assimilation method into a real test case, in order to investigate the value of remote sensing-derived water levels for improving model reliability, variational methods are used in assimilation. By using a variational data assimilation approach (4D-var), we have investigated potential contributions of SAR derived spatially distributed water levels for the identification of time-independent parameters (Manning coefficients) in a shallow-water flood model.

The spatially distributed water levels have been derived from a SAR image by employing the method developed in [START_REF] Hostache | Estimation de niveaux d'eau en plaine inondée à partir d'images satellites radar et de données topographiques fines[END_REF]. They have been obtained with a ± 40 cm mean uncertainty, using a RADARSAT-1 image of the 1997 flood event of the Mosel river. This has been possible by using both an analysis of the relevance of SAR derived flood extent limits for hydraulic purposes, and a merge between the relevant limits and a highly resolution DEM, under hydraulic coherence constraints inspired from [START_REF] Puech | Using geographical information system and aerial photographs to determine water levels during floods[END_REF][START_REF] Raclot | Remote sensing of water levels on floodplains: a spatial approach guided by hydraulic functioning[END_REF]. Such a water level estimation provides spatially distributed information at the time of a satellite overpass while classical in situ measurements are punctually available.

Furthermore, numerical experiments conducted in this study indicate that a rather dense information in space is of great benefit for the identification of unknown parameters (Manning friction coefficient). Indeed, the assimilation of the SAR derived water levels, in addition to an incomplete discharge hydrograph, proves to be capable of identifying Manning friction coefficients, while the ground data alone does not allow such an identification. Furthermore, a sensitivity analysis conducted by using the SAR derived water levels,

shows that a spatial distribution of the friction coefficient based on land-cover may not necessarily lead to the better model results. Indeed, these water levels, used as a guide in the sensitivity analysis, can define areas of Manning friction homogeneity, without apparent link with land-cover. Such sensitivity analysis may finally base the Manning parameter spatial distribution more on the model hydraulic functioning, than on the land-cover.

The most important issue in flood forecast is to predict the flood peak. We know that remote sensing images cannot be useful in predicting directly flood peaks except in regions with very slow flood events. Our process leads to a better model, more physically constrained, which helps to reduce equifinality, and can be used for future predictions.

This study contributes to the basic understanding of the assimilation of spatial distributed water levels into a shallow-water flood model.

  . These have been provided respectively by the Environment Management Direction of Lorraine (DIREN Lorraine), the French Company for Electricity (EDF) and the German Water Management and Navigation Direction (Wasser und Schifarts Direktion SüdWest). The discharge hydrograph recorded in Perl hydrometric station looks contradict those recorded in EDF and Uckange hydrometric stations (no mass balance of water in the Mosel

  able" constant value, and the available data at inflow and outflow boundaries are used to perform a run of the forward model. Next, the simulation results are compared with the available measurements, namely the water stage hydrograph at the middle gauge station (EDF) and the SAR image derived water levels in the floodplain at the satellite passover time. The computational domain includes the main channel and the floodplain from Uckange to Perl. It is meshed using triangles and quadrilaterals elements and contains 2340 cells and 2430 nodes. The simulation time window begins at 12:00, Feb.25, 1997 and ends at 12:00, Mar. 2, 1997. Initial conditions are computed as the steady state with constant inflow discharge at start time. At upstream boundary, the observed discharge hydrograph is imposed; at downstream, observed water stage hydrograph is imposed.

  value n being set by land-cover class. Four various decompositions have been investigated, depending on the total class number: a) a lumped n value: one constant value all over the domain. b) three classes (distributed n) consisting of the main channel and the left and right floodplains (left and right banks present few differences in landuse.

. 11 .

 11 These gradient values are computed for a constant Manning value over the whole domain (n =0 .033); sensitivity results are similar if computed from different Manning coefficient values (n =0.05 or n =0.01 for example).

Fig. 2 .

 2 Fig. 2. Recorded hydrographs during the flood event at 3 different locations

Fig. 3 .Fig. 5 .

 35 Fig. 3. Observations available: discharge at low water level at EDF gauge station and one (partial) image

Fig. 6 .

 6 Fig. 6. Land cover: 10 classes (Case d) in the text).

Fig. 7 .

 7 Fig. 7. Convergence of the minimization process. Plotted: the normal cost function J (solid line), the normal gradient dJ (dashed line). Case considered: Manning calibration with synthetic data for 3 classes (case b).

Fig. 8 .Fig. 9 .

 89 Fig. 8. Water levels (m) at image blocks i.e. where image information is available (see Fig. 4). Vertical bars corresponds to the measures with estimated uncertainties. The four curves corresponds to the preliminary forward run and the three calibrated models responses (depending on the land classes number considered; VDA-nXX = 1, 3 and 10 land classes). (Recall: VDA=Variational Data Assimilation = 4D-VAR algorithm; direct = forward model without calibration.)

Fig. 10 .

 10 Fig. 10. Measured discharge hydrographs (m 3 /s) at the middle gauge station EDF, and computed ones using the calibrated model (Manning is calibrated in 3 land classes, case b)). "Z": computed values with elevation imposed at downstream; "QZ-org": computed values with water stage hydrograph imposed at downstream.

Table 1

 1 Difference between the computed and the measured water depth (m) in floodplain (data extracted from the image).

	It demonstrates

Table 2

 2 Twin experiments. Identified Manning coefficients n. Example based on three landcover classes: case b).

	Land class True n	Land use	Computed n
	1	0.025	Main channel	0.0251
	2	0.04	Right floodplain	0.0406
	3	0.05	Left floodplain	0.0495

Table 4

 4 Sensitivity analysis. Gradient values for case d) 10 land-uses (at Manning value constant everywhere n =0.033)

	Channel	Vegetation	Bridge	Gravel	Grassland
	881.54 e+6 43.75 e+6	17.76 e+6	3.64 e+6	-1.06 e+6
	Left	Right	Snaked land	Urban	Downstream channel
	0.07 e+6	0	0	0	0

Fig. 11. Sensitivity analysis. Gradient value: cost function with respect to the Manning coefficient if 10 land blocks are considered (case d).

Acknowledgment

This research has been supported by Région Rhône-Alpes, France ( "numerical prevention for floods" project [2003][2004][2005][2006].

The authors would like to acknowledge André Paquier (Cemagref Lyon, France) for his help in the definition of the preliminary forward run, Joel Marin (IN-RIA Rhône-Alpes) for his help in computational issues. The authors would like to thank the reviewers for their detailed comments that has helped to improve the paper. English writing has been improved thanks to National Natural Science Foundation of China (Grant No. 50709034).

Preliminary forward run. Comparison of discharge (m 3 /s) at the middle gauge station EDF: measured, computed with elevation imposed at downstream ("Z"), computed with water stage hydrograph imposed at downstream ("QZ-org").


Don Joven Agravante

Andrea Cherubini

Abderrahmane Kheddar

Using vision and haptic sensing for human-humanoid joint actions

Keywords: Physical Human-Robot Interaction, Human and humanoid skills/cognition/interaction, Human-Robot Collaboration

Human-humanoid haptic joint actions are collaborative tasks requiring a sustained haptic interaction between both parties. As such, most research in this field has concentrated on how to use solely the robot's haptic sensing to extract the human partners' intentions. With this information, interaction controllers are designed. In this paper, the addition of visual sensing is investigated and a suitable framework is developed to accomplish this. This is then tested on examples of haptic joint actions namely collaboratively carrying a table. Additionally a visual task is implemented on top of this. In one case, the aim is to keep the table level taking into account gravity. In another case, a freely moving ball is balanced to keep it from falling off the table. The results of the experiments show that the framework is able to utilize both information sources properly to accomplish the task.

I. INTRODUCTION

Humanoid robots provide many advantages when working together with humans to perform various tasks. This is because humans have an extensive experience in physically collaborating with each other. Humanoids take advantage of this with their human-like range of motion and sensing capabilities. This reduces a human collaborator's need to learn how to interact with the robot. However, many challenges are still present in the various research areas that study physical human-robot collaboration. The particular area of interest in this paper is the integration of vision and force information to enable "human-robot haptic joint actions". These are collaborative tasks requiring both parties to physically interact with each other. It has two important aspects:

1) both robot and human are doing jointly the same task, 2) a haptic interaction exists. An example of such a task is joint carrying/transportation of large objects [START_REF] Yokoyama | Cooperative works by a human and a humanoid robot[END_REF]- [START_REF] Bussy | Human-humanoid haptic joint object transportation case study[END_REF]. This is illustrated in Fig. 1 along with three important information sources: prior task knowledge, vision and force. In this example, the haptic interaction exists through the object -the haptic channel. This means that a force/torque applied on one end of the object is felt by the partner on the other end. Because of this, previous research has focused primarily on regulating interaction forces for safety. The use of vision for the robot has not 1

CNRS-UM2

LIRMM UMR 5506, Interactive Digital Human group, 161 Rue Ada, 34392 Montpellier, France {firstnames.lastname}@lirmm.fr.

2 CNRS-AIST, JRL (Joint Robotics Laboratory), UMI 3218/CRT, Intelligent Systems Research Institute, AIST Central 2, Umezono 1-1-1, Tsukuba, 305-8568, Japan. been investigated before in this context and this is the main focus of this paper. Finally, the prior task knowledge can be used as a guideline on how vision and/or force information should be used for the task. Aside from just being able to do the task, the desire is to make the robot proactive in helping the human [START_REF] Bussy | Human-humanoid haptic joint object transportation case study[END_REF]- [START_REF] Gribovskaya | Motion learning and adaptive impedance for robot control during physical interaction with humans[END_REF]. A proactive behavior aims to lighten the load on the human, as opposed to being a passive follower -only reacting to human motion without any anticipation. A variety of ways are being researched to program proactive behaviors in joint actions. For example in [START_REF] Maeda | Human-robot cooperative manipulation with motion estimation[END_REF] a minimum jerk model of human motion is used to have a good guess of the human's intention. In [START_REF] Gribovskaya | Motion learning and adaptive impedance for robot control during physical interaction with humans[END_REF], machine learning methods are used to obtain a proactive behavior. Another method is studying humans doing the task. An example of this is presented in [START_REF] Bussy | Human-humanoid haptic joint object transportation case study[END_REF], where it was found that a human dyad moves in constant velocity phases when doing joint transportation of large objects (as in Fig. 1). These findings were then applied to a controller to make the humanoid robot more proactive [START_REF] Bussy | Human-humanoid haptic joint object transportation case study[END_REF].

As previously mentioned, most works in this field use data from force/torque sensors only to do the task and make the robot pro-active. Although it is impressive what can be achieved using force data alone (i.e. a "blind" robot), vision is clearly required for some tasks and could possibly help to make the robot proactive. For example in collaboratively carrying a table, force data alone cannot be used to see if an object on top of the table is about to fall down. It is clear that visual information is largely complementary to force information (analogous to the human senses of sight and touch). Combining these might enable a humanoid to perform more complicated tasks, similar to a human [START_REF] Ernst | Merging the senses into a robust percept[END_REF]. Although the benefits are great, there are not many established methods integrating vision and force control.

To start detailing this paper's contribution, the concept of human-robot haptic joint actions is first elaborated in Section II. Following this, a short review on the state of the art in combining vision and force is presented in Section III. The proposed framework to be used is then detailed in Section IV. This is then applied to two case studies in Section V. The results of these are shown in Section VI. Finally Section VII ends with a brief summary and plans for future work.

II. HUMAN-ROBOT HAPTIC JOINT ACTIONS

Haptic joint actions are mainly characterized by a collaborative task and a haptic interaction throughout this task. This applies to several collaborative tasks that humans do together. For example: helping someone to their feet, the process of handing off objects, supporting someone while walking, etc. Another common example is carrying a large object together. This is illustrated in Fig. 1. This work investigates the addition of visual information to such tasks (shown in Fig. 1 by the field-of-view). To better show our contribution with respect to previous works, a few examples of these are shown in Fig. 2. The most commonly studied scenario of human-robot collaboration is making the robot a pure follower as in Fig. 2(a). This illustrates that the task is only known to the human leader a priori. Through the haptic channel, the human communicates his intentions to the robot. With a properly designed controller that takes into account interaction forces, the robot can follow the human's lead [START_REF] Bussy | Human-humanoid haptic joint object transportation case study[END_REF]. This is represented by the bigger arrow in the haptic channel from human→robot. Furthermore, the human is able to obtain visual and haptic data as to how good/bad the robot follows the task.

Another situation is the exact opposite of the robot being a pure follower: it can be the leader of the task. This is illustrated in Fig. 2(b). An example of such a scenario is presented in [START_REF] Bussy | Proactive behavior of a humanoid robot in a haptic transportation task with a human partner[END_REF] where a joystick is used to give the robot direct information on the task. Although a 2nd human provides this task information via the joystick, in the context of collaborative carrying, the robot is the leader. Without knowing the task of the robot, the human partner then tries to help the robot carry the object to a desired location. It is fair to assume that the human uses both the sense of sight and touch to achieve this task. Apart from the two illustrated examples, other possibilities exist. For example, a combined leader-follower strategy has also been developed [START_REF] Evrard | Homotopy switching model for dyad haptic interaction in physical collaborative tasks[END_REF], [START_REF] Kheddar | Human-robot haptic joint actions is an equal controlsharing approach possible?[END_REF]. The concept for this strategy is that a sharing of roles can be possible where one is both leader and follower (with varying degrees).

Finally, it can be noticed that in these previous works the robot is solely reliant on haptic information. No visual data is used in the task as illustrated in Fig. 2 by the lack of the robot field-of-view. The aim of the work here is to move towards the general case of Fig. 1, particularly adding vision as another information channel. In order to do this, the main question is how does one combine vision and force information in the context of haptic joint actions. The context presents significant challenges specifically by having a human in the loop.

III. STATE OF THE ART ON COMBINING VISION AND FORCE FOR MANIPULATION

To our knowledge, there has been no previous work on combining vision and force in the context of human-robot haptic joint actions. However, previous works on combining vision and force for object manipulation tasks can be found in the literature. A brief review of these methods is given here as context for the choice of the control framework detailed in the next section.

The different approaches to combining vision and force information for manipulation tasks can be classified into three general categories [START_REF] Nelson | Improved force control through visual servoing[END_REF]. These are: traded, hybrid and shared. The simplest strategy is termed as "traded" since the final control output is traded between a purely forcebased controller and a purely vision-based controller. The switching between controllers depends on the task execution. A common example of this strategy starts by doing a "guarded move" visual servoing phase. A contact event causes the guarded move phase to stop. After this, the task is completed by force control.

The second category is termed as "hybrid" methods. These are hybrid in the sense that the vision-based controller and force-based controller act at the same time but in different spaces. This requires prior specification of a "taskframe" [START_REF] Mason | Compliance and force control for computer controlled manipulators[END_REF], [START_REF] Baeten | Integrated vision/force robotic servoing in the task frame formalism[END_REF]. The task frame is a Cartesian reference frame that can be divided into vision-controlled and forcecontrolled directions. Doing this decouples vision and force into orthogonal spaces. After this, the controllers can be designed separately and work independently in their own predefined space.

The third and final category is termed as "shared" methods. These are so-called since there is no separation in time (the case of traded control) or in space (the case of hybrid control). The control responsibility is shared by both vision and force throughout the operation. By doing this, all available information can be used [START_REF] Nelson | Improved force control through visual servoing[END_REF]. An example of a shared method is described in [START_REF] Prats | Vision force control in task-oriented grasping and manipulation[END_REF], where force feedback is used to correct the visual servo control trajectory.

In the context of human-robot haptic joint actions, a shared method is needed. With the human in the loop, safety is always a top priority. Therefore, the control must always be compliant in all degrees of freedom (DOFs). This can be achieved by always making use of the force information in all DOFs. Therefore, in order to make use of visual data, a shared control strategy must be used. A good candidate for this is the impedance control framework [START_REF] Hogan | Impedance control -An approach to manipulation. I -Theory. II -Implementation. III -Applications[END_REF]. Impedance control allows a manipulator to be compliant by defining a virtual mechanical impedance. Vision can be used in this framework to provide a reference trajectory that is tracked in the absence of external forces [START_REF] Santis | Human-robot interaction control using force and vision[END_REF], [START_REF] Morel | Impedance based combination of visual and force control[END_REF]. This approach fits the criteria that the robot is always compliant even when moving according to the visual information.

IV. OVERVIEW OF THE ROBOT CONTROL FRAMEWORK

As mentioned in the preceding section, a vision-based controller coupled to an impedance controller is the core concept that is used in this work. However, there are a few more modules that need to be defined to make everything work on a humanoid robot. These are briefly detailed in this section for the sake of completion. As a guide, Fig. 3 shows the overall workflow in a scenario where the goal of the human and robot is to transport the plate while preventing the ball from falling off. This same scenario was tested for the framework here and is discussed more in the next sections as one of the case studies.

A. Visual Tracking Algorithm

To start the whole control framework, visual information needs to be processed. In the example depicted in Fig. 3, this algorithm gives data about the ball on top of the table. The raw data is a combination of an RGB image and a depth map, obtained from an ASUS Xtion sensor mounted as the "eyes" of the humanoid. This raw data needs to be processed into an estimate of the pose or position to be controlled by the vision-based controller. This is done by a tracking a salient visual feature throughout the image sequence and extracting the needed pose or position information from this. However, there is no generic visual tracking algorithm that can work for any and all cases. Furthermore, some knowledge of the task is needed to know what is the important visual information that is needed. This is the reason why the case studies presented here use trivial objects -a cube of known size and color and a ball of known color. However visual tracking is a fairly well-developed field in computer vision and a wide range of algorithms have been developed as reported in this extensive survey [START_REF] Lepetit | Monocular model-based 3D tracking of rigid objects[END_REF]. A future plan is to improve in this area by making use of more state-of-the art algorithms but this is not the focus of this paper. But it should be noted that the main constraints for the vision algorithm are robustness and speed since it needs to be used with a controller. Therefore, visual trackers such as [START_REF] Comport | Realtime markerless tracking for augmented reality: the virtual visual servoing framework[END_REF] that are designed with this in mind are preferred.

B. Vision-Based Control

After the pertinent information has been extracted from the images, this data needs to be used in a controller that does the desired task. Visual servoing is a well-studied research field that uses the visual information directly into a control law [START_REF] Chaumette | Visual servo control. i. basic approaches[END_REF]. However, as with the visual tracking algorithm, a visual servoing controller needs to be specifically designed depending on the task. Hence a prior knowledge of the task is also needed. In the two case studies to be shown, the tasks are different and as such different vision-based controllers are created. The details of these are presented in the next section.

C. Impedance Control

Impedance control is a general framework that regulates the contact interaction of the robot and its environment according to a defined virtual mechanical impedance [START_REF] Hogan | Impedance control -An approach to manipulation. I -Theory. II -Implementation. III -Applications[END_REF].

In the context of human-robot collaboration, this means that interaction forces between the human and robot are regulated at all times. A general formulation of impedance control is:

f = M( Ẍd -Ẍ) + B( Ẋd -Ẋ) + K(X d -X) (1)
where f is the wrench (force-torque screw) composed of the force and torque vectors. The vectors X d , Ẋd and Ẍd are a desired pose and its first and second derivative. Correspondingly, vectors X, Ẋ and Ẍ represent an actual pose and its first and second derivative. Finally, matrices M, B and K are the inertia, damping and stiffness parameters that define the desired virtual mass-spring-damper system [START_REF] Hogan | Impedance control -An approach to manipulation. I -Theory. II -Implementation. III -Applications[END_REF].

The general formulation in Eq. ( 1) can account for "different" control methods. Firstly, it can either be "impedance controlled" by controlling forces/torques (f ) or "admittance controlled" by controlling pose (X). We use the latter, since the force-torque sensors on the HRP-2 wrists provide the signal to be used for f and it is inherently positioncontrolled. Secondly, a variety of specific methods can be defined, depending on the values of M, B, and K. Here, all three matrices M, B, K are utilized. These are tuned from experiments to achieve the desired compliant behavior.

To completely define the impedance controller, the desired pose and its derivatives X d , Ẋd and Ẍd are needed. These are just the output of the vision-based controller as shown in Fig. 3. Therefore in the absence of external forces, the vision-based controller achieves its control objectives.

D. Stack-of-Tasks

The Stack-of-Tasks is a generalized inverse kinematics abstraction layer [START_REF] Mansard | A versatile generalized inverted kinematics implementation for collaborative working humanoid robots: The stack of tasks[END_REF]. As its name implies, the main advantage it gives is the hierarchical organization of different tasks to be executed. This allows efficient integration and abstraction of the different humanoid robot tasks. In the context of collaborative object carrying depicted in Fig. 3, the visionbased control and impedance control are used to servo the humanoid hands. While this is done, other tasks are also in the stack, namely: collision and joint limit avoidance, a task to maintain a good posture and control the robot's center of mass (COM) and finally a task to control the legs for walking [START_REF] Herdt | Walking without thinking about it[END_REF]. The Stack-of-Tasks framework is able to directly control all the humanoid robot's joints considering all the tasks and their priority. Since critical tasks have a higher priority, their execution is ensured [START_REF] Mansard | A versatile generalized inverted kinematics implementation for collaborative working humanoid robots: The stack of tasks[END_REF].

V. CASE STUDIES

To test the framework described in Section IV, two simple cases of human-robot haptic joint actions are created that clearly benefit from visual information. Because the task of collaborative "table-carrying" has been well-studied within the group, this is used as the base task. An object is then placed on top of the table and the additional task is concerned with this object and the table tilt angles (φ x , φ y ). A simplified side-view of the task in Fig. 4 shows φ y and its relation to the height difference z r and the table length l t . Furthermore, three important reference frames are drawn in this image to facilitate the explanations that follow -the control frame {cf }, a local reference frame on the robot {l} and the table frame {t}. The control for the robot can be done just by defining the pose l T cf . This if justified by assuming a rigid grasp during the whole task. This means that the pose of the right and left hands: cf T rh and cf T lh are constant throughout the task and generating the 2-handed control is just a change of frame.

A. Stationary Object -keeping the plane level

In this scenario, a green cube is placed on top of the table as a representative "stationary object". To visually track this object and get its pose, the model based tracker of the ViSP software package [START_REF] Comport | Realtime markerless tracking for augmented reality: the virtual visual servoing framework[END_REF] is used. The additional goal here is to keep this object upright with respect to the gravity field. This implies that the table plane must be kept level. More formally, the objective can be thought of as φ x = 0, φ y = 0. The case of φ x is trivial in this regard (for as long as no conflicting haptic intention is applied by the human). The case of φ y is more interesting as it depends directly on the human as seen in Fig. 4. From the figure, it can be seen that there are 2 ways to affect φ y . First is applying a rotation and hence a torque on the controlled end. Since the table is rigid, the human needs to react to this appropriately by changing his z. This can be uncomfortable for the human, hence the other option is used here which is the exact opposite -the robot changes z r . It should also be noted that a rotational compliance in both partners is needed for such a control strategy. With this, the model of the task is then:

l t sin φ y = z r . (2)
The final control law is obtained by using the method to derive visual servoing controllers described in [START_REF] Chaumette | Visual servo control. i. basic approaches[END_REF]. By differentiating the model and setting an exponential decrease of the error ė = φy = -λφ y the control law of Eq. ( 3) is obtained where λ is a gain parameter.

żr = -λ φy cos φy

To use this in the impedance control framework, numerical integration is done such that z r (t) = z r (tdt) + żr (t)dt where dt is the time step of the controller. A piece-wise constant velocity is assumed making the acceleration null.

B. Moving Object -keeping a ball from falling off

In this scenario, a ball is placed on top of the table. Any disturbance will tend to make the ball fall off the table. In this case, a specific visual tracking method was developed for the ball using well-established computer vision algorithms. Using both rgb and depth information, the ball's position relative to a pre-defined goal location on the table is tracked visually. With this information, the controller can then be defined.

The scenario described here is the classical "ball-on-plate" example problem in control theory. This is well-known to be marginally stable with just a proportional control law but it is also known that a proportional-derivative (PD) control can be enough which is what is tested in the work presented here. Furthermore, the same argument about controlling φ y from the stationary object case can be made here. Hence the control is done in z r and φ x . With these considerations, the control laws are shown as follows:

z r = K p,x (x d -x) -K d,x ẋ φ x = K p,y (y d -y) -K d,y ẏ, (4)
In Eq. ( 4), (x, y) pertain to the ball position, (x d , y d ) its desired location and ( ẋ, ẏ) the error derivatives considering a constant desired location. K represent the control gains which are experimentally tuned.

To use this result in the impedance control framework, numerical differentiation is done to obtain the velocity terms. However, the acceleration is again considered null to prevent noise from differentiation.

VI. RESULTS

The case studies of Section V were tested experimentally on the HRP-2 humanoid robot. These were also integrated with previous work in the group [START_REF] Bussy | Proactive behavior of a humanoid robot in a haptic transportation task with a human partner[END_REF], [START_REF] Bussy | Human-humanoid haptic joint object transportation case study[END_REF] that does the collaborative transportation task. Snapshots of the experiments for the stationary object case is shown in Fig. 5 and the experiments for the moving ball are shown in Fig. 6. In these experiments, the walking gait of the humanoid produces an external disturbance. But even with this, the designed task is achieved showing that the approach presented works.

A. Stationary Object -keeping the plane level

To verify the vision controller design, a step response of the error (φ y ) was obtained in experiments. This is shown in Fig. 7. The top plot shows that the error decreases exponentially to 0. This is the expected behavior from setting ( ė = -λe). This result proves that the implemented system works (although with some noise and noticeable latency of the vision algorithm).

B. Moving Object -keeping a ball from falling off

To show the performance for this case, the resulting ball trajectory while walking is plotted in Fig. 8. This also corresponds to the snapshots of Fig. 6 (bottom). This data is from the visual tracker used and as such is an estimate of the actual result. However, it should be noted that in this experiment, the ball did not fall while walking with the human even when changing directions (which is the difficulty for the faster ball). The results here are from the slower yellow ball, and the desired ball position is set to (0.15, 0), which is 15cm closer to the human than the table centroid. The red border signifies a rough approximate of the table boundaries. This result show that although the ball moves a lot, it doesn't fall off the table during the transportation task.

VII. CONCLUSION AND FUTURE WORK

In this paper, a framework to combine vision and haptic information in human-robot haptic joint actions is presented. The core idea of the framework is using visionbased controllers to define a desired virtual trajectory for the impedance controller which regulates interaction forces. Two case studies are presented here with different visionbased controllers and both are shown to be able to achieve the task. These experiments verify the chosen approach. This is also just preliminary steps for the planned integration of visual information into collaborative tasks. A next area of research is in increasing the proactive behavior of the system by detecting visual cues from the partner.

Fig. 1 .

 1 Fig.1. A diagram of human-robot haptic joint actions. In this general case, both human and robot have a complete knowledge of the task (represented by the blue arrows). Furthermore, each uses both vision (green) and haptic (red) information to achieve this task.

Fig. 2 .

 2 Fig. 2. Particular cases of Human-Robot Haptic Joint Actions

Fig. 3 .

 3 Fig.3. An illustration of the overall workflow of the algorithms running in the robot to enable it to collaborate with the human using both vision and force information.

Fig. 4 .

 4 Fig. 4. A simplified "thin beam" model used to control the table tilt through the height

Fig. 7 .

 7 Fig. 7. Step response of the vision controller. Top: Plot of the error (φy). Bottom: Plot of l z cf

Fig. 5 .Fig. 6 .

 56 Fig. 5. Snapshots of an experiment for the stationary object (green cube)

Fig. 8 .

 8 Fig. 8. Controlled ball trajectory (blue) during the experiment, computed from the visual data. The red border is a rough estimate of the table edges

VIII. ACKNOWLEDGEMENT

This work is supported in part by the FP7 IP Robo-How.Cog project (www.robohow.eu). FP7-ICT-2011-7 Contract No 288533.