Nicolas Baghdadi
email: nicolas.baghdadi@teledetection.fr

Nicolas Lemarquand
email: nicolas.lemarquand@teledetection.fr

Hani Abdallah
email: hani.abdallah@teledetection.fr

Jean Sté Phane Bailly
email: jean-stephane.bailly@teledetection.fr

The Relevance of GLAS/ICESat Elevation Data for the Monitoring of River Networks

Keywords: lidar, ICESat, rivers, lakes, elevation

The Ice, Cloud and Land Elevation Satellite (ICESat) laser altimetry mission from 2003 to 2008 provided an important dataset for elevation measurements. The quality of GLAS/ICESat (Geoscience Laser Altimeter System) data was investigated for Lake Leman in Switzerland and France by comparing laser data to hydrological gauge water levels. The correction of GLAS/ICESat waveform saturation successfully improved the quality of water elevation data. First, the ICESat elevations and waveforms corresponding to water footprints across the transition from the land to water were analyzed. Water elevations (2 to 10 measurements) following the land-water transition are often of lesser quality. The computed accuracy for the ICESat elevation measurements is approximately 5 cm, excluding transitions footprints, and 15 cm, including these footprints. Second, the accuracy of ICESat elevation was studied using data acquired on French rivers with a width greater than the size of the ICESat footprint. The obtained root mean square error (RMSE) for ICESat elevations in regard to French rivers was 1.14 m (bias = 0.07 m; standard deviation = 1.15 m), which indicates that small rivers could not be monitored using ICESat with acceptable accuracy due to land-water transition sensor inertia.

Introduction

Freshwater provisioning will be a major problem in future decades. Currently, 1 billion people depend on lakes for domestic water consumption, which will expand to 5.5 billion in 30 years [START_REF] Calmant | Continental surface waters from satellite altimetry[END_REF]. Monitoring of continental water resource, particularly lakes and rivers, is necessary because many global regions lack sufficient water for domestic consumption, agriculture, and other uses. Hydrologic networks that are generally organized on a national basis provide monitoring of water resource. These networks measure the temporal variations in the water level of rivers, lakes and reservoirs. However, the spatial distribution of hydrological gauges is insufficient in some parts of the world due to freshwater needs. Moreover, a reduction in the number of stations has been observed in conjunction with a decline in the quality of measurements. A better understanding of the global water cycle and the impact of climate change requires a thorough knowledge of water resources on the continental surface. Therefore, there is a real demand for a global, homogeneous, perennial monitoring system of continental water.

Optical and Synthetic Aperture Radar (SAR) imageries (e.g., LANDSAT, NOAA, MODIS, SPOT, ASTER, ERS, ASAR, RADARSAT, and PALSAR/ALOS) are often used to map the extent of water areas [START_REF] Hui | Modelling spatial-temporal change of Poyang Lake using multitemporal Landsat imagery[END_REF][START_REF] Andreoli | Inland Lake Monitoring Using Low and Medium Resolution ENVISAT ASAR and Optical Data: Case Study of Poyang Lake (Jiangxi, PR China)[END_REF]. Determination of changes in the associated water level is estimated using radar altimeters and lidar systems [START_REF] Koblinsky | Measurement of river level variations with satellite altimetry[END_REF][START_REF] Birkett | Surface water dynamics in the Amazon Basin: Application of satellite radar altimetry[END_REF][START_REF] Coe | Calculation of river discharge and prediction of lake height from satellite radar altimetry: Example for the Lake Chad basin[END_REF]. The radar altimetry missions demonstrated the potential of continental water observation with three principal advantages suitable to hydrological observations: all-weather operability, global data coverage, and temporal repetivity (10 days for Jason and Topex and 35 days for ENVISAT and ERS). The radar altimeters used for Topex/Poseidon, ERS-2, ENVISAT, and Jason-1/2 have a coarse spatial resolution "footprint size" (a few hundred meters to few kilometers), which means they can monitor only large lakes, rivers, and reservoirs. Typically, altimetric measurements can range in accuracy from a few centimeters (e.g., Great Lakes, USA) [START_REF] Morris | Evaluation of the TOPEX/Poseidon altimeter system over the Great lakes[END_REF][START_REF] Birkett | Synergistic remote sensing of Lake Chad: Variability of basin inundation[END_REF] to tens of centimeters (e.g., Lake Chad, Africa) [START_REF] Birkett | Synergistic remote sensing of Lake Chad: Variability of basin inundation[END_REF]. Accuracies of 5 cm RMS for Lake Issykkul (Kirgyzstan) and 10 cm RMS for Lake Chardarya (Kazakhstan) were previously determined by Cré teaux and Birkett [START_REF] Cré Taux | Lake studies from satellite radar altimetry[END_REF]. In all of these studies, the obtained accuracies are directly related to the size of the body of water because the satellite measurement compared to ground reality is usually obtained by averaging individual shot measurement along tracks. Therefore, rivers appear more difficult to study because of their width. The first studies relating to rivers were performed in the Amazonian basin, which contains very large rivers. In studies by Birkett et al. [START_REF] Birkett | Surface water dynamics in the Amazon Basin: Application of satellite radar altimetry[END_REF], RMSE (Root Mean Square Error) values between hydrologic gauges and Topex/Poseidon data for a selection of rivers and floodplains in the Amazon Basin exhibited high variability, with an overall mean of approximately 1.1 m RMS and optimal values of approximately 40 cm. The accuracy for rivers is highly variable, being that these measurements are dependent on the target width.

Satellite laser altimetry, such as the Geoscience Laser Altimeter System (GLAS) on the Ice, Cloud, and Land Elevation Satellite (ICESat), provides elevation data that are very useful to many applications and offer the potential to monitor small water areas due to small associated footprints (50 to 90 m). NASA launched the ICESat satellite in January 2003 [START_REF] Schutz | Overview of the ICESat mission[END_REF][START_REF] Zwally | ICESat's laser measurements of Polar ice, atmosphere, ocean and land[END_REF]. ICESat includes two laser altimeters operating at visible (green, 532 nm) and near-infrared (1,064 nm) wavelengths, and each near-infrared laser footprint along the satellite track is 172 m apart, with 40 points (shots) per second (40 Hz). The horizontal geolocation accuracy of the ICESat footprints is 3.7 m. The waveform of each shot is sampled in 544 or 1,000 bins over land area at a temporal resolution of 1 ns. The vertical resolution of land waveforms is 15 cm. The nominal pointing angle of ICESat is approximately 0.3° off nadir, with a potential of up to 5°. The high pointing angles reduce the specular reflection on waters, and therefore, the waveform saturation will increase the elevation error [START_REF] Urban | A survey of ICESat coastal altimetry applications: Continental coast, open ocean island, and inland river[END_REF]. However, GLAS laser altimeter cannot be used for routine monitoring because only some tracks are possible with ICESat due to power safety and meteorological conditions, which is in contrast to radar altimeters where study sites are revisited every 10 days.

Chipman and Lillesand [START_REF] Chipman | Satellite-based assessment of the dynamics of new lakes in southern Egypt[END_REF] assessed the accuracy of the ICESat GLAS laser altimeter on the Toshka lakes west of Lake Nasser in southern Egypt. The standard deviation of GLAS elevation measurements along a single transect was on the order of 3-8 cm. Bhang et al. [START_REF] Bhang | Verification of the vertical error in C-band SRTM DEM using ICESat and Landsat-7, Otter Tail County, MN[END_REF] have demonstrated a difference between GLAS elevation and gauging stations between 2 and 35 cm for lakes in Otter Tail County, Minnesota (USA). Other studies have demonstrated that the accuracy of ICESat is greater than 10 cm, compared to gauge data [START_REF] Urban | A survey of ICESat coastal altimetry applications: Continental coast, open ocean island, and inland river[END_REF][START_REF] Martin | ICESat range and mounting bias estimation over precisely-surveyed terrain[END_REF][START_REF] Braun | ICESat Laser Altimetry in the Great Lakes[END_REF]. However, the quality of ICESat data for rivers has not been thoroughly investigated. Analysis of ICESat measurements (transect of 50 km with 300 footprints/shots) along the Tapajos River, Brazil, provides a RMSE of approximately 3 cm for clear meteorological conditions, 8 to 15 cm for partly cloudy skies, and 25 cm for heavy clouds [START_REF] Urban | A survey of ICESat coastal altimetry applications: Continental coast, open ocean island, and inland river[END_REF]. These studies are indicative of the performance of ICESat for large-size targets, with high spatial coverage. The accuracy of lidar measurements should be degraded for small targets, with only some measurements by track, as in Birkett et al. [START_REF] Birkett | Surface water dynamics in the Amazon Basin: Application of satellite radar altimetry[END_REF], where a factor of ten was observed with radar altimetry between RMSE obtained from rivers and floodplains. Could the elevation of rivers and lakes that span several hundreds of meters be measured with good accuracy (approximately 10 cm) using ICESat (footprints of 50-90 m)?

Moreover, the number of measurements excluded from each ICESat transect corresponding to footprints with possible disturbances related to the land-water transition is never specified. However, the value of this margin will limit the size of targets that could be analyzed by laser systems, such as ICESat.

The response of GLAS/ICESat to the land-water transition was first analyzed. The accuracy of ICESat for Lake Leman (a flat surface) was then evaluated. Finally, the potential use of ICESat measurements for rivers in France (small targets) was studied.

Lake Leman

Description of Study Site and Dataset

Study Site

The accuracy of ICESat measurements was evaluated across Lake Leman located in Switzerland and France (Lat. 46°26′N and Long. 6°33′E). This is one the largest lakes in Western Europe with a surface of 582 km 2 . The average level of water is 372 m, which is controlled by the Seujet Dam near Geneva. The maximum length and width of the lake are 73 and 14 km, respectively. Lake Leman is bordered on the south by the Haute-Savoie massif. The land-water transition in the ascending acquisition mode of ICESat involves a difference in elevation of 550 m for a distance of 5 km (south to north). In the area north of Lausanne, the land-water transition (north to south) in descending mode is softer with a difference in elevation of 260 m for a distance of 5 km.

Dataset

A total of 21 GLAS (Geoscience Laser Altimeter System) transects across Leman Lake were available for comparison with the lake gauge data (Figure 1). ICESat data were acquired between 2 March 2004 and 13 March 2008. The length of transects varied between 8.5 and 10.5 km, with a total number of 774 shots measured. The data used are comprised of GLA01 (waveform) and GLA14 (Land Surface Altimetry) products in the Release 31. The elevation of GLAS laser data was computed using tools developed by the US National Snow and Ice Data Center (NSIDC). Six transects corresponding to cloudy episodes were next excluded from the dataset because the corresponding elevations were strongly disturbed (heights of approximately 3,500 m). Two hydrological stations, Chillon and Saint Prex (Figure 1), were used to evaluate the accuracy of ICESat elevations. As the difference in water level between the two stations is relatively weak (~1 cm), an average elevation was used. The Swiss Federal Office for the Environment (FOEN) [START_REF]Federal Office for the Environment (FOEN)[END_REF] measures water levels at permanent gauging stations for Swiss bodies of water. The water levels of hydrological stations are given in the Swiss height measurement reference system. The reference for all height measurements in Switzerland is the "Repère Pierre du Niton" in the harbor of Geneva (stone). The height of this stone was evaluated in 1902 to be 373.6 m over sea level.

All heights refer to the same vertical datum to conduct a consistent comparison among the two height data sets. In this study, the data comparison was performed in terms of orthometric heights with respect to the WGS84 reference system and the EGM96 geoid model. ICESat elevations available to users correspond to ellipsoidal heights with reference to the Topex ellipsoid. The conversion of the Topex ellipsoid to the WGS84 ellipsoid was approximated by h WGS84 = h Topex -70.7 cm. Next, ellipsoidal elevations (h WGS84 ) were converted into orthometric elevations (H) by applying the EGM96 geoid value (N): H = h WGS84 -N.

Shots of water were then extracted individually track by track, and the mean elevation was calculated for each track. This mean was used for comparison with reference elevations (hydrological gauges).

Results

Analysis of ICESat waveforms demonstrated a saturation phenomenon, which must be corrected. Saturation occurs when the returned energy by a number of 1 ns bins is greater than the threshold function of gain [START_REF] Urban | A survey of ICESat coastal altimetry applications: Continental coast, open ocean island, and inland river[END_REF][START_REF] Abshire | Geoscience Laser Altimeter System (GLAS) on the ICESat mission: On-orbit measurement performance[END_REF]. Therefore, elevations calculated from saturated waveforms are incorrect. A saturation correction is available for the produced data and is recommended for use on calm waters because of the associated specular reflection. The correction value varies between 0 and 1.5 m and is added to the GLAS elevation. However, the correction is not available for some shots and the value "-999.00" is given instead (i.e., elevation should be corrected, but the method used to calculate the correction coefficient cannot provide a correct value). Figure 2 shows some GLAS waveforms obtained for Lake Leman in the case of water only. For the saturated waveforms (calm water), the high level of reflected energy involves a clipping of the maximum peak, which shifts this peak several bins to the left. The derived GLAS elevation is thus underestimated.

Figure 3 shows the ICESat elevations across Lake Leman for the 3 June 2005 transect. This example shows the importance of the saturation effect (errors reach 1.41 m) and demonstrates the need to include the saturation correction.

With correction Without correction Reference elevation Invalid correction

Analysis of ICESat transects across the lake demonstrates that the shots from the land to water could be divided into four categories:

 Category 1: Shots corresponding to land. The saturation correction is valid (between 0 and 1.50 m).

 Category 2: Three to eight shots correspond to "land", "land+water" or "water" with invalid correction.

 Category 3: Two to ten shots correspond to "water" with valid correction, but not adapted (called transition shots).

 Category 4: ICESat shots corresponding to "water" with estimated elevations close to reference elevations. The saturation correction is valid and well adapted.

On the passage of ICESat from land to water, the first laser waveforms corresponding to water exhibit a strong saturation. The elevation of theses shots systematically exhibit a significant variation compared to the reference elevation (Figure 4). Indeed, GLAS elevations of water surfaces without saturation correction are lower than reference elevations from -40 to -90 cm (maximum values). With correction of the saturation effect, elevations are higher than reference elevations from +20 to +50 cm (maximum). The highly saturated elevations without available correction (-999.00) are lower than the reference elevation of approximately one meter (Figure 4(b)). This saturation phenomenon which affects 2 to 10 spots is more important in ascending mode than in descending mode. The instrumental parameters are not different in ascending and descending modes. Shots following the first saturated waveforms provide elevations very close to the reference elevation with weak fluctuations. However, the analysis of the terrain geometry shows that the transition land-lake in ascending mode is very different from the transition land-lake in descending mode. In the South, when the Haute-Savoie massif is thrown into the lake, the transition land-lake in ascending mode is brutal with uneven of 550 m on a distance of 5 km. Conversely, in the Lausanne area in North, the transition land-lake in descending mode is softer with uneven of 260 m on 5 km. The difference in the observations between ascending and descending transects could be related to the laser beam target reflectivity.

Saturation is due to specular reflection of the water that saturates the detector, depending on laser beam target reflectivity. ICESat profiles demonstrate that a certain time is necessary for the gain control to compensate for electronic saturation by a reduction of gain. This could explain the progressive return to normal that can last for 10 measurements (maximum duration of 0.25 seconds and a maximum distance about 1.5 km). Moreover, the time necessary for desaturation of the detector depends on the ICESat pointing angle. Pointing angles larger than nominal will reduce the probability of waveform saturation but will increase the elevation error. The relationship between the pointing angle () and the elevation error () is given by K.  [START_REF] Urban | A survey of ICESat coastal altimetry applications: Continental coast, open ocean island, and inland river[END_REF], where K is a constant = 5cm/°/arcsecond and  is the pointing knowledge error in arcseconds. For our ICESat dataset, the pointing angle was close to 0.3° (nominal value). For  for a fully-calibrated laser campaign is = 1.5 arcseconds (L3a laser campaign), the elevation error  = 2.25 cm.

This slowness in adaptation of gain after the land-water transition phase seriously limits the ICESat for limnimetric monitoring of rivers. With the exception of great rivers, the detector has not enough time after each land-water transition to correctly adjust the saturation correction coefficient. With or without correction, the first estimates of ICESat elevations would have an altimetric error of several tens of centimeters.

The 19 March 2006, transect demonstrates successively the first shots of the transition slightly saturated (valid correction), oversaturated (not valid correction = -999.00), and then slightly saturated (valid correction) (Figure 4(b)). This alternation of "valid correction" and "not valid correction" is surprising. There can be a restitution error for the spot position on the ground related to pointing error. Alternatively, the laser direction can oscillate to approach nadir in certain cases resulting in saturation, whereas the pointing angle is sufficient to avoid saturation in other cases. This second scenario would explain the nonsystematic nature of the phenomenon. Table 1 summarizes the results of the comparison between GLAS elevations (GLA14 product) and lake level data recorded by hydrological stations. The mean lake height for each ICESat transect is determined from averages of all valid height measurements. The mean value difference between the ICESat and reference elevations varies from 0 to 1.09 m for GLAS data without the correction of saturation effect and from 0 to 22 cm for corrected data. The RMS error ranges from 5 cm to 1.19 m without correction and from 5 cm to 30 cm with correction (Table 1). For nine ICESat tracks, the RMSE is greater than 20 cm when the correction is not applied against 2 tracks (Table 1). The mean RMSE for all tracks is approximately 15 cm with the saturation correction and 5 cm excluding transition footprints (first saturated waveforms). The errors are higher when the saturation correction is not applied (33 cm). Transects acquired on 13 March 2008, and 6 April 2007, had high standard deviations in comparison to other transects, which could result from a lack of sufficient correction for atmospheric effects. For the 13 March 2008 transect, only 22 data were used in the accuracy calculation. This result from implementing the saturation correction for continental waters corroborates the previous observations of Urban et al. [START_REF] Urban | A survey of ICESat coastal altimetry applications: Continental coast, open ocean island, and inland river[END_REF].

Rivers of Metropolitan France

The accuracy of GLAS elevations was then studied using data acquired for French rivers wider than the size of ICESat footprint (Figure 5(a)). This comparison allows for analysis of the accuracy of GLAS measurements for rivers and to conclude on the use of ICESat for these small bodies of water. Only GLAS data from laser number 3 (campaign L3) were used because this area had the smallest major axis footprint size (approximately 55 m compared with 90 m for L2 and 149 m for L1), which is better adapted to the width of the French river. 10 ICESat data collection campaigns were used (L3A to L3J) (Table 2). All selected ICESat footprints had the following characteristics: (1) intersecting rivers with a width greater than 50 m, (2) radius of 1 km for gauge stations with a well referred limnimetric scale, and 3) distance far enough from hydraulic systems that can disrupt the water level within a small distance. Moreover, optical aerial photos were used to validate the selection of ICESat footprints. For this comparison, 26 ICESat tracks and 46 footprints were obtained (Figure 5(b)). This weak number for ICESat measurements illustrates the difficulty of monitoring small rivers with the ICESat configuration.

Assuming that measurements at gauge stations were of sufficient accuracy regardless of distance in time and space between GLAS shot and gauge station measurement, the elevation of each GLAS shot was then compared with the one from the closest hydrological station at the closest measured time. For this latter point, the obtained mean deviation between the time of GLAS acquisitions and the stations was 1.8 hours with a standard deviation (std) of 2.25 hours. The RMSE for GLAS elevations obtained on the French rivers is 1.14 m (bias = 0.07 m; std = 1.15 m; 46 footprints) (Figure 6). This corresponds to approximately 10 times the RMSE for wider bodies of water, such as Lake Leman. For concerned ICESat waveforms, 61% were saturated and the correction of this saturation was not possible (correction = -999,000). Moreover, the correction proposed for the first water shots following the passage of ICESat from land to water is not appropriate (data before correction are underestimations, while data after correction are overestimations). The difference between the errors observed for Lake Leman in comparison to French rivers is close to the results that Birkett obtained with altimetry radar for Lake Chad [START_REF] Birkett | Synergistic remote sensing of Lake Chad: Variability of basin inundation[END_REF] and rivers in the Amazon basin [START_REF] Birkett | Surface water dynamics in the Amazon Basin: Application of satellite radar altimetry[END_REF].

Conclusions

GLAS elevation precision for Lake Leman was similar to previously published values in the literature (RMSE between 5 and 30 cm with a mean value of approximately 15 cm). The saturation correction for GLAS elevation was necessary for consequent improvement of error (bias, RMSE). A precision greater than 15 cm for GLAS elevation was possible for bodies of water greater than 1.5 km. Indeed, the first ten footprints for only water after passage of laser from land to water are often saturated with an inappropriate saturation correction. These results demonstrated that several GLAS measurements are necessary to ensure that the GLAS system provides reliable saturation correction. GLAS temporal profiles reveal a slow progressive adaptation of the GLAS sensor before correct elevations can be proposed. This adaptation can reach 0.25 s (10 shots). Therefore, elevation measurements performed for the first hundred meters, which follow the land-water transition, are biased. For example, this adaptation distance is approximately 1.5 km for ten transition shots.

The GLAS precision was then evaluated for data acquired from rivers in France. The determined RMSE was 1.14 m (bias ± std = 0.07 m ± 1.15 m). This result demonstrates the difficulty of using GLAS data for rivers that are relatively smaller in width (inappropriate saturation correction).Therefore, a factor of 10 was determined for the GLAS precision between large lakes and small rivers.

The computed accuracy for ICESat water elevation data can make these data inconsistent depending on their hydrological use. When using water levels data to monitor the dynamic of water balance for large reservoirs, assimilation of water levels data with a metric accuracy can reduce errors of estimated water balance by a factor of two [START_REF] Pereira-Cardenal | Real-time remote sensing driven river basin modeling using radar altimetry[END_REF]. But according to studies on fluvial hydraulic simulation model performances, the computed accuracy for ICESat water elevation is of poor interest [START_REF] Schumann | Progress in integration of remote sensing-derived flood extent and stage data and hydraulic models[END_REF]: for flood modeling in Europe occurring typically in low slopes fluvial environments, a decimeter accuracy of observed water levels, consistent with hydraulic output accuracies and fluvial slopes is required [START_REF] Straatsma | The Effect of Uncertain Floodplain Roughness in Hydrodynamic Flood Simulations[END_REF]. Assimilation of data with only decimeter accuracy adds meaningful value to hydraulic modeling, especially if this data is not frequent in space [START_REF] Alsdorf | Measuring surface water from space[END_REF], which is the case of ICESat data.

Moreover, as the accurate ICESat measurements require significant water body width (hundreds of meters), the ICESat water level measurements are only adapted to major rivers and large catchment scales, or flood plains during extreme events.

Figure 1 .

 1 Figure 1. The Ice, Cloud and Land Elevation Satellite (ICESat) tracks through Lake Leman. (a) Ascending mode and (b) descending mode.

Figure 2 .Figure 3 .

 23 Figure 2. Typical Geoscience Laser Altimeter System (GLAS) waveforms on Lake Leman (water only). (a) Unsaturated waveform, (b) waveform slightly saturated, and (c, d) waveforms with high saturation. These examples are extracted from the 03 June 2005 transect. The reference elevation is 372.18 m. Ordinate = energy in volts and abscissa = bin number (bin spacing = 1 ns).

Figure 4 .

 4 Figure 4. Some examples of ICESat transects in comparison to reference elevation (gauges). Abscissa = shot number (spacing = 0.025 s) and ordinate = elevation (ICESat with/without saturation correction, and hydrological gauges).

Figure 5 .

 5 Figure 5. (a) ICESat tracks (in gray) and French hydrographic network (in blue), (b) GLAS data used in this study (rivers width >50 m, far from any dam or other hydraulic plant, with distance between ICESat shot and hydrological stations <1,000 m).

Figure 6 .

 6 Figure 6. Difference between GLAS and hydrological stations elevations (reference) for French rivers. Each point corresponds to one footprint (for a total of 46 footprints).

Table 1 .

 1 Comparison of Geoscience Laser Altimeter System (GLAS) and reference elevations (hydrological gauges) for Lake Leman. The eight dates marked by a star (*) have weak or null saturation corrections.

	Date	Hydrological gauges	ICESAT Dataset		ICESat elevation (m) Without correction		ICESat elevation (m) With correction
	dd/mm/yy	Reference elevation (m)	Shot number	mean	std	RMSE ICESat-Ref	mean	std	RMSE	ICESat-Ref
	02/03/04 *	371.72	60	371.72 0.05	0.05	0.00	371.72 0.05	0.05	0.00
	12/06/04 *	372.23	47	372.19 0.11	0.11	-0.04	372.24 0.12	0.12	0.01
	04/03/05 *	371.74	61	371.78 0.11	0.12	0.04	371.71 0.09	0.09	-0.03
	16/03/05	371.68	42	370.59 0.50	1.19	-1.09	371.77 0.16	0.18	0.09
	03/06/05	372.18	61	371.49 0.59	0.90	-0.69	372.14 0.09	0.10	-0.04
	15/11/05 *	372.25	49	371.99 0.10	0.27	-0.25	372.05 0.13	0.23	-0.19
	19/03/06	371.82	43	371.53 0.22	0.36	-0.29	371.74 0.15	0.17	-0,08
	07/06/06 *	372.17	63	372.24 0.07	0.10	0.07	372.26 0.08	0.12	0.09
	18/06/06	372.26	45	372.19 0.21	0.22	-0.07	372.42 0.11	0.19	0.16
	07/11/06 *	372.22	63	372.23 0.13	0.13	0.01	372.29 0.08	0.11	0.07
	19/11/06	372.22	50	372.00 0.34	0.40	-0.22	372.18 0.12	0.13	-0.04
	25/03/07 *	371.74	61	371.81 0.10	0.12	0.07	371.87 0.07	0.15	0.13
	06/04/07	371.66	48	371.45 0.28	0.35	-0.21	371.69 0.19	0.19	0.03
	16/10/07	372.21	59	372.02 0.24	0.30	-0.19	372.20 0.07	0.07	-0.01
	13/03/08 *	371.53	22	371.31 0.21	0.30	-0.22	371.31 0.21	0.30	-0.22

Table 2 .

 2 Start date and End date of ICESat campaigns L3A to L3J.

	Laser

campaign Start date (dd/mm/yyyy) End date (dd/mm/yyyy) L3A

		10/03/2004	11/08/2004
	L3B	02/17/2005	03/24/2005
	L3C	05/20/2005	06/23/2005
	L3D	10/21/2005	11/24/2005
	L3E	02/22/2006	03/28/2006
	L3F	05/24/2006	06/26/2006
	L3G	10/25/2006	11/27/2006
	L3H	03/12/2007	04/14/2007
	L3I	10/02/2007	11/05/2007
	L3J	02/17/2008	03/21/2008

Acknowledgements

The authors would like to thank the ICESat team and the National Snow and Ice Data Center for providing ICESat data and expertise. The authors would also like to thank the French National Water Observation Service "EauFrance" for providing access to the national gauge database.


Fabien Racapé

Simon Lefort

Dominique Thoreau

Marie Babel

Olivier Déforges

Fabien Racape
email: fabien.racape@technicolor.com

Olivier Deforges

Characterization and adaptive texture synthesis-based compression scheme

INTRODUCTION

In state-of-the-art compression schemes, pixel-wise redundancy is reduced by using predictions and transformed domain operations. However, classical spatio-temporal approaches, exploiting redundancy based on the mean squared error (MSE) criterion, are not able to take visual redundancy into account. Detailed textures may seem nearly stationary for the Human Visual System, but totally irregular according to MSE criterion. At the same time, texture synthesis algorithms [START_REF] Wei | Fast texture synthesis using treestructured vector quantization[END_REF][START_REF] Ashikhmin | Synthesizing natural textures[END_REF][START_REF] Kwatra | Graphcut textures: image and video texture synthesis using graph cuts[END_REF][START_REF] Kwatra | Texture Optimization for Example-based Synthesis[END_REF] have shown promising results. The purpose of new coding schemes is to detect regions where exact positions of texture patterns are irrelevant for the human eye. The whole regions are not encoded since a few patterns are sufficient enough to synthesize satisfactory regions.

This paper presents a framework that can be integrated in standard coding schemes. One the main novelties of this scheme comes from the fact that texture algorithms have to be parameterized in order to give their best results. Thus, the scheme includes a texture characterization step based on DCT-domain descriptors that outputs the feature size of texture patterns to be synthesized. The remainder of the paper is organized as follows. In section 2 existing such schemes are presented. In section 3 is introduced the proposed framework. In sections 4 and 5 are described encoder and decoder designs in detail. Section 6 finally presents some experimental results.

RELATED WORK

One of the first synthesis-based compression scheme was presented in [START_REF] Dumitras | A texture replacement method at the encoder for bit-rate reduction of compressed video[END_REF]. This approach includes an analysis in order to detect replaceable textures which are finally synthesized by a dedicated texture synthesis algorithm at decoder side.

An interesting framework has been designed in [START_REF] Zhu | Video coding with spatio-temporal texture synthesis[END_REF] where some 8x8 blocks are removed at the encoder side and synthesized by the decoder. The segmentation first classifies blocks into structures, corresponding to boundaries of objects, and textures. The first category is classically encoded and the latter is removed and synthesized at decoder. The segmentation is based on simple edge detector thresholding. To avoid temporal inconsistencies, motion estimation is considered in selecting patches which will be used for texture synthesis. The latter is performed using the algorithm presented in [START_REF] Kwatra | Graphcut textures: image and video texture synthesis using graph cuts[END_REF], which appeared hard to exploit in our tests, since 8x8 blocks offer poor overlap for the Graphcut technique. The work presented in [START_REF] Ndjiki-Nya | Generic and robust video coding with texture analysis and synthesis[END_REF] proposes a closed-loop analysis-synthesis approach. As in [START_REF] Zhu | Video coding with spatio-temporal texture synthesis[END_REF], groups of pictures (GOP) are considered for a spatio-temporal scheme. Each potential region from a GOP is both analyzed and synthesized, using side information from texture analyzer. A first synthesizer is designed for Rigid textures with global motion, which has a great similarity to global motion compensation (GMC). Another synthesizer, inspired by the patch-based approach developed in [START_REF] Kwatra | Graphcut textures: image and video texture synthesis using graph cuts[END_REF], processes Non-rigid textures with local and global deformations. The scheme has recently been upgraded [START_REF] Ndjiki-Nya | Optimization of video synthesis by means of cost-guided multimodal photometric correction[END_REF] with photometric corrections, using Poisson editing and covariant cloning. Another framework for image compression is proposed in [START_REF] Liu | Image compression with edge-based inpainting[END_REF] where both textural and structural blocks can be removed at encoder side. The textures are synthesized at decoder using [START_REF] Kwatra | Graphcut textures: image and video texture synthesis using graph cuts[END_REF], while a binary edge information helps an inpainting method to retrieve structural ones. To keep color coherency, a draught board of blocks is preserved, limiting the bit-rate saving. The problem of these gradated regions is tackled in [START_REF] Xiong | Block-based image compression with parameter-assistant inpainting[END_REF] with parameter assistant inpainting. However, several pixel-based synthesis algorithms from the literature [1, 2, 4] also provide promising visual performance for a large range of textures, until they get the right parameters depending on the patterns to synthesize. Thus this work proposes to explore their efficiency in the compression context, by adding a coherent segmentation and a preprocessing step for texture synthesis that provides such parameters.

FRAMEWORK OVERVIEW

The proposed framework is depicted in figure 1. First, input images are analyzed to detect textured regions that can be synthesized by the algorithm used at decoder side. Thus this part of the framework is divided into two main steps: a segmentation step and a characterization step. A resulting texture map is sent to the encoder and side information, describing texture patterns, can potentially been sent. Then, pointed out regions are partially removed, whereas structural regions are classically encoded. Removed regions are synthesized from neighboring small surfaces of texture, also classically encoded. In the following, the patches denote these samples of texture which serve as inputs for synthesis. At the decoder side, the bit-stream and potential side information to locate encoded regions are both used to build and locate structural parts, whereas the labeled removed regions are synthesized using a new adaptive pixel-based algorithm. Like encoder side analysis, the synthesis contains a characterization step, which outputs required parameters for texture algorithm.

The next section describes more precisely the analysis done at the encoder side.

TEXTURE ANALYSIS

The goal is here to provide a complete and coherent texture analysis process, containing an adapted block-based segmentation and texture characterization.

LAR segmentation

Since the proposed framework requires a robust segmentation process that can be adapted to our context, the contentbased codec LAR (Locally Adaptive Resolution) presented in [START_REF] Deforges | Color LAR codec: a color image representation and compression scheme based in local resolution adjustment and self-extracting region representation[END_REF] has been chosen for its region handling functionality. A spatial coder partitions images into variable-size blocks, relying on a homogeneity criterion based on a gradient thresholding. A quadtree topology is adopted in order to get a non-overlapping distribution of blocks. Then, a segmentation step merges adjacent blocks. This step uses a hierarchical approach and non-symmetrical distances taking into account the surface areas. The distance between two adjacent regions is the average between:

• Mean weighted distance, which is computed from mean Y values, • Gradient distance, which corresponds to the local gradients along shared borders.

Two parameters steer the segmentation. First, the described distance is hierarchically thresholded. Then, a chromatic control involving the mean distance of chroma components between adjacent regions supervises the merge operation. Figure 2 gives, on the left side, an example on the Wool frame where regions are represented with their mean value. The following describes the adaptation of the LAR-region output map to the framework constraints.

LAR Regions Removed texture regions

Figure 2: Adapted regions, black regions are removed to be synthesized at decoder side.

Adapting labels to texture synthesis compliant maps

In order to be compliant with block-based standard compression schemes, segmented regions are aligned on a blockbased grid, i.e. 8x8 or 16x16 blocks. Each block containing several regions is considered as structure like in [START_REF] Zhu | Video coding with spatio-temporal texture synthesis[END_REF]. The others are considered as texture and will be removed by the encoder. Since the grid follows 16x16 macroblocks (MB), no additional information is required when using h.264 encoder. Binary side information may be sent in other configurations to locate structure blocks. Figure 2 shows such regions in the Wool frame with a 16x16 grid. Large textured regions require anchor preserved blocks in order to prevent artefacts. After experiments, it has been decided to encode some anchor macroblocks at strategic locations depicted in figure 2 since the region size exceeds a fix number of MBs in width or height. Coupled with a confidence-based synthesis order described in next section 5.2, they prevent visible seams at synthesis junctions. These segmented textures are then characterized to get input information about texture patterns for synthesis.

Texture characterization

This part of the scheme is used at both encoder and decoder sides. It outputs the parameters required by the texture synthesizer in order to give its best results. Thus, using characterization enables the encoder to decide whether the texture can be synthesized or not. Indeed, if the characterized texture requires impossible parameters, for example too large patterns, it is finally classically encoded. The main parameter for our pixel based algorithm corresponds to a size of neighborhood window. Thus the characterization step computes descriptors from different sizes. The descriptors are derived from those described in [START_REF] Smach | Generalized Fourier descriptors with applications to objects recognition in SVM context[END_REF] which are computed from the Fourier transform. Since the characterization requires little sizes of blocks, i.e. inferior to 32 pixels wide, descriptors from DCT domain have been chosen. Indeed, the Fourier transform outputs descriptors with a lower resolution in frequency, which is a problem to describe a signal from a small window. Like in [START_REF] Smach | Generalized Fourier descriptors with applications to objects recognition in SVM context[END_REF], descriptors are computed from concentric circles represented in figure 3, following:

D DCT (λ ) = π/2 θ =0 |C (λ , θ ))| 2 (1)
where C(λ , θ ) corresponds to the DCT coefficient at the location (λ cosθ , λ sinθ ), except they are integrated on a quarter of circle. The DC coefficient is not taken into account in order to be invariant to changes of average luminance, so the descriptors vector has size of block minus one coefficient.

A large set of blocks inside the texture region are randomly chosen to compute descriptors at different sizes centered at the same position.

Since the scheme uses discrete transform, figure 4 depicts the integral computation on a particular coefficient at position (i, j). The DCT value is weighted with the length L of the arc of a circle, crossing position (i, j) from a to b, given by L = R * θ where

θ = arcsin 1 - b 2 R 2 -arcsin a R . (2)
The average descriptors are then computed for each size and analyzed. According to [START_REF] Wei | Fast texture synthesis using treestructured vector quantization[END_REF] and our experiments, the size of the compared neighborhood has to be greater than the elementary pattern of the texture, to produce a visually good result. The first coefficient of each block corresponds to a single variation of luminance over the block used for DCT computation. Experiments on a large set of texture patches show that if the first coefficient is greater than any others in the descriptors vector, the block size is smaller than an elementary pattern. Indeed, the first coefficient of each block corresponds to a single variation of luminance over the block used for DCT computation. The output size parameter corresponds to the minimum size NxN of DCT block that provides a non-monotone decreasing vector, which is given by

N out = min N B N×N DCT |D(1) = max λ =[2;N-1] {D(λ )} ( 3
)
where B is the DCT block of size N. Thus, computing descriptors at random locations will serve to approximate the feature size of texture patterns. Typically, descriptors of sizes from 4x4 up to 16x16 pixels are computed for each location.

The next section describes the synthesis of the segmented regions, using parameters from texture characterization.

TEXTURE SYNTHESIS

The proposed synthesizer is based on the work presented in [START_REF] Wei | Fast texture synthesis using treestructured vector quantization[END_REF]. This algorithm is pixel-based since the output surface is processed one pixel at a time. It relies on the matching between current pixel neighborhood and those in a texture patch. The criterion used for matching is the Sum of Square Error (SSE). Although this algorithm provides pleasing visual results for texture synthesis only, it has to be welladapted to compression. First, it has been designed to rely on a given texture patch only, whereas our synthesis aims at exploiting all the available data (decoded blocks). Then, algorithm parameters such as synthesis order and neighborhood size and shape have to be reconsidered.

Adaptive neighborhood size using texture characterization

Texture characterization enables the synthesizer to get an approximate neighborhood. In order to refine the neighborhood sizes for matching, a set of neighborhood sizes are tested, for example 7x7 and 9x9 if the best descriptor is 8x8 large. The chosen distance minimizes the norm distance

D m,n (p, c) = ∑ N m,n k=0 (p k -c k ) 2 N m,n (4)
where N m,n is the number of pixels contained in the neighborhood of current pixel, m and n respectively representing its height and width, c k and p k denote the kth pixel's luminance in the current neighborhood and the considered one in the patch respectively. Considering various neighborhood sizes enables to better catch texture patterns size and shape. This process is illustrated in figure 5 where three neighborhood sizes {(3; 5)(4; 7)(5; 9)} are competing.

Experiments show that a synthesis with a lot of changes in neighborhood sizes does not give visually good results. In order to avoid this kind of issue, we propose to favor the neighborhood sizes that have been chosen by a majority of previously synthesized pixels, by mean of a global weight.

Synthesis scan order

The goal is here to exploit available data which are, at this point, the only confident data. The latter corresponds to surrounding pixels, which are whether previously decoded or synthesized. Thus, raster scan order is clearly not adapted. As in [START_REF] Liu | Image compression with edge-based inpainting[END_REF], a scan order depending on a synthesis-coherent confidence map is adopted. The map building law is depicted in figure 6, where initial computing from previously decoded pixels is represented on the left side and the map update during synthesis on the right side. First, previously decoded pixels are initialized with a unique confidence value. In order to be coherent with the synthesis, the neighborhood size used for texture synthesis serves for computing the confidence of available pixels. Thus, this confidence value ω is computed ω(i + k, j + l).

1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 0 0 0 0 0 0 0 0 1 1 1 1 1 1 1 1 1 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 1 1 1 1 1 1 1 1 1 1 1 1 0 0 0 0 0 0 0 0
(

) 5
After having synthesized the most confident pixel, the map is updated by assigning the previously computed ω(i, j), which is shown in figure 6 b).

Texture patch design

Texture synthesizers from the literature process with a rectangular patch. Its relevance is essential to ensure an efficient synthesis in our context. In the literature, patch choices are usually either not described or correspond to cropped proximate blocks [START_REF] Zhu | Video coding with spatio-temporal texture synthesis[END_REF]. However, when removing large regions, one can note that cropping a patch closed to the texture region to be synthesized can lead to inconsistencies between candidate pixels and previously reconstructed borders. For instance, figure 8 a) shows a selected region in which there is a variation of luminance, leading to a failed synthesis. That is, the patch contains the texture around the region to be synthesized, which is depicted in figure 8 b). The patch has to be large enough, according to the size of neighborhood. A compromise is done at encoder side, if texture region are too small and patterns too large, then the region is not classified as texture. Reversely, large texture regions require anchors in order to prevent pixel-based synthesis to fail, so we decide to keep some MBs since the size exceeds a fix width or height. Figure 7 depicts the use of such anchors, following the confidence-based scan order.

Neighborhood optimization

In order to improve the research process, we propose to use the exact neighborhood matching described in [START_REF] Sabha | Texture synthesis using exact neighborhood matching[END_REF]. A shortlist of candidate neighborhoods in the patch is first computed for the current location. It corresponds to the set of neighborhoods in the patch containing at least one direct neighbor which exactly matches respective current pixel neighbor. Then the SSE distance is computed for the all neighborhoods in the shortlist and the best candidate is chosen like in [START_REF] Wei | Fast texture synthesis using treestructured vector quantization[END_REF]. The only change lies in the scan order which defines the causal direct neighbors. This technique enables a great adap- tation to our context in which high frequency discontinuities must be avoided. Some results are provided in the next section, in particular visual comparisons between the persented pixel-based synthesis and a patch-based method.

EXPERIMENTAL RESULTS

The framework has been implemented within the JM joint model [START_REF]JM reference software version 14[END_REF]. Presented results are given for intra frames sequences Wool, Container, Coastguard and Soccer. A set of four Quantization Parameters (QP){15,20,25,30} is considered. The segmentation is aligned on a 16x16 grid. The retained patch has been set to 16 pixels of width and three sizes of neighborhoods {(4,7); [START_REF] Dumitras | A texture replacement method at the encoder for bit-rate reduction of compressed video[END_REF][START_REF] Liu | Image compression with edge-based inpainting[END_REF]; [START_REF] Zhu | Video coding with spatio-temporal texture synthesis[END_REF][START_REF] Deforges | Color LAR codec: a color image representation and compression scheme based in local resolution adjustment and self-extracting region representation[END_REF]} are competing. Table 1 shows the bit-rate saving at similar visual quality for the different sequences and QP. One can see the promising performance, i.e. the scheme achieves up to 20% bit-rate saving on Soccer and 14% on Container. One notes the other good results on Coastguard containing a large part of water texture. Figure 9 shows the visual results for the Coastguard sequences encoded with QP = 20, where 22.2% of the frame has been removed and synthesized. One can note comparable results for most of textures, for example with the trodden earth in the Morocco Trial sequence. However, results on water in the Coastguard sequence points out that the pixel-based technique avoids edge artefacts for smooth textures.

CONCLUSION AND FUTURE WORK

A framework based on a new adaptive texture synthesis has been presented. The synthesis is pixel-based since it relies on the comparison between current pixel neighborhoods and those in an atypically shaped patch. In order to better catch texture patterns, a texture characterization is carried out to parameterize the synthesizer, in order to choose the best candidate among different neighborhood sizes and shapes.

Figure 1 :

 1 Figure 1: Framework overview.

Figure 3 :Figure 4 :

 34 Figure 3: Descriptor computed on a 4x4 DCT block.

Figure 5 :

 5 Figure 5: Finding best matching neighborhood.

Figure 6 :

 6 Figure 6: Confidence map order.

Figure 7 :

 7 Figure 7: Pixel-based synthesis order.

Figure 9 :

 9 Figure 8: Texture patch design.

  Figure 10 aims at comparing this pixel-based scheme and the patch-based from[START_REF] Kwatra | Graphcut textures: image and video texture synthesis using graph cuts[END_REF] on the coastguard and Morocco Trial sequences.

Figure 10 :

 10 Figure 10: Comparison between pixel-based and patch-based synthesis.

Table 1 :

 1 The framework has been validated within an H.264/AVC video codec. Experimental results show significant bit-rate saving compared to H.264/AVC. Ongoing researches focus on the characterization of detected textures to improve synthesizer adaptation as well as combining patch-based and pixel-based approaches. Coastguard CIF 35.1 33.3 31.1 27.0 Soccer 4CIF 20.8 20.9 20.8 20.6 Resulting bit-rate saving.

	QP	15	20	25	30
	Sequences		Bit-rate saving (%)
	Wool SD	11.3 9.6	7.9	6.2
	Container CIF	14.7 11.7 8.4	3.4