J M Roger

B Palagos

D Bertrand

E Fernandez-Ahumada

CovSel : Variable selection for highly multivariate and multi-response calibration

I theory

1. algorithm 1. Searching index I 1 corresponding to the predictor closest to the responses, by :

I 1 = ArgMax i (x T i YY T x i ) (1)
2. Collinear information is removed from X and Y by orthogonal projection :

X ← P ⊥ x I 1 X (2)
Y ← P ⊥ x I 1 Y (3)
This process is then repeated for I 2 , I 3 , • • •, I k .

properties

Equation 1 can be written as :

I 1 = ArgMax (diag (X T YY T X)) (4)
Furthermore it can be demonstrated that this equation is equivalent to :

I 1 = ArgMax i Max v,v 2 =1 cov (x i , Yv) 2 (5)
CovSel is a PLS like selection process Test results on the four responses

III results on corn dataset

λ (nm) moisture oil protein starch assignement 1100 × × × × baseline 1190 × × oil 1306 × × 1428 × × × starch 1500 × NH 1592 × × × × 1718 × × × × oil 1886 × × × × 1940 × × × × water 2106 × × × starch 2204 × × × × 2250 × × × starch 2306 × × × × oil 2388 × × 2498 × × × × baseline
Spectral assignment of the selected variables CovSel yields parcimonious, meaningful and little correlated selections A calibration and a validation sets were randomly drawn in the proportion of 2/3 and 1/3, respectively. CovSel was applied on the calibration set, with a predefined number of variables k = 15. Four models were then optimized individually on each response and applied to the validation set.

-Set Wine Grapes : -X = 250 Vis/VNIR spectra of wine grain samples × 256 wavelengths (310 to 1100 nm).

-Y = 250 membership degrees to 3 varieties : carignan (crg), grenache blanc (grb) and grenache noir (grn) ; e.g. y = [0, 1, 0] for a sample belonging to class 2.

A calibration and a validation sets were randomly and equally drawn. The variables selected by CovSel were used as input of a Linear Discriminant Analysis. The observation of the leaveone-out cross-validation results allowed the determination of the optimal number of selected variables. The discriminant model calibrated on this subset was applied on the test set.

IV results on wine discrimination CovSel is efficient on discrimination problems

  Illustration of the iterative erosion of X T YY T X

-

  CovSel is a variable selection method well suited to highly multivariate and multi-response calibration -CovSel acts as the PLS, splitting the covariance between X and Y -Covsel can be used on discrimination problems -CovSel is particularly adapted to the design of multispectral devices II material and methods -Set Corn (http ://software.eigenvector.com/Data/Corn) :-X = 80 NIR spectra of corn samples × 700 wavelengths (1100 to 2498 nm).-Y = 80 reference values of 4 responses : moisture, oil, protein and starch

  and cross-validation errors of the LDA model, as a function of the number of CovSel steps. k = 8 variables were retained. the LDA model when applied on the test set Scores of the LDA calculated on the selected variables.


Marie Chavent
email: .chavent@u-bordeaux2.

Vanessa Kuentz

Benoit Liquet

Jerome Saracco

Package 'ClustOfVar'

Keywords: ClustOfVar Type, Package Version, 0.5

Cluster analysis on a set of variables. Variables can be quantitative, qualitative or a mixture of both. Missing values are replaced by means for quantitative variables and by zeros in the indicator matrix for qualitative variables.

matsim boolean, if 'TRUE', the matrices of similarities between variables in same cluster are calculated.

Value var a list of matrices of squared loadings i.e. for each cluster of variables, the squared loadings on first principal component of PCAmix. For quantitative variables (resp. qualitative), squared loadings are the squared correlations (resp. the correlation ratios) with the first PC (the cluster center).

sim a list of matrices of similarities i.e. for each cluster, similarities between their variables. The similarity between two variables is defined as a square cosine: the square of the Pearson correlation when the two variables are quantitative; the correlation ratio when one variable is quantitative and the other one is qualitative; the square of the canonical correlation between two sets of dummy variables, when the two variables are qualitative. sim is 'NULL' if 'matsim' is 'FALSE'.

cluster a vector of integers indicating the cluster to which each variable is allocated.

wss the within-cluster sum of squares for each cluster: the sum of the correlation ratio (for qualitative variables) and the squared correlation (for quantitative variables) between the variables and the center of the cluster.

E

the pourcentage of homogeneity which is accounted by the partition in k clusters.

size the number of variables in each cluster.

scores a n by k numerical matrix which contains the k cluster centers. The center of a cluster is a synthetic variable: the first principal component calculated by PCAmix. The k columns of scores contain the scores of the n observations units on the first PCs of the k clusters.

Author(s)

Marie Chavent <marie.chavent@u-bordeaux2. hclustvar

Hierarchical clustering of variables

Description

Ascendant hierarchical clustering of a set of variables. Variables can be quantitative, qualitative or a mixture of both. The aggregation criterion is the decrease in homogeneity for the cluster being merged. The homogeneity of a cluster is the sum of the correlation ratio (for qualitative variables) and the squared correlation (for quantitative variables) between the variables and the center of the cluster which is the first principal component of PCAmix. PCAmix is defined for a mixture of qualitative and quantitative variables and includes ordinary principal component analysis (PCA) and multiple correspondence analysis (MCA) as special cases. Missing values are replaced by means for quantitative variables and by zeros in the indicator matrix for qualitative variables.

Usage hclustvar(X.quanti = NULL, X.quali = NULL) Arguments X.quanti a numeric matrix of data, or an object that can be coerced to such a matrix (such as a numeric vector or a data frame with all numeric columns).

X.quali a categorical matrix of data, or an object that can be coerced to such a matrix (such as a character vector, a factor or a data frame with all factor columns).

Value

Description

Iterative relocation algorithm of k-means type which performs a partitionning of a set of variables.

Variables can be quantitative, qualitative or a mixture of both. The center of a cluster of variables is a synthetic variable but is not a 'mean' as for classical k-means. This synthetic variable is the first principal component calculated by PCAmix. PCAmix is defined for a mixture of qualitative and quantitative variables and includes ordinary principal component analysis (PCA) and multiple correspondence analysis (MCA) as special cases. The homogeneity of a cluster of variables is defined as the sum of the correlation ratio (for qualitative variables) and the squared correlation (for quantitative variables) between the variables and the center of the cluster, which is in all cases a numerical variable. Missing values are replaced by means for quantitative variables and by zeros in the indicator matrix for qualitative variables.

kmeansvar Usage kmeansvar(X.quanti = NULL, X.quali = NULL, init, iter.max = 15 , nstart = 1, matsim = FALSE) Arguments X.quanti numeric matrix of data, or an object that can be coerced to such a matrix (such as a numeric vector or a data frame with all numeric columns).

X.quali a categorical matrix of data, or an object that can be coerced to such a matrix (such as a character vector, a factor or a data frame with all factor columns).

init either the number of clusters or an initial partition (a vector of integers indicating the cluster to which each variable is allocated). If init is a number, a random set of (distinct) columns in X.quali and X.quanti is chosen as the initial cluster centers.

iter.max the maximum number of iterations allowed.

nstart if init is a number, nstart corresponds with the number of random sets used in the process. matsim boolean, if 'TRUE', the matrices of similarities between variables in same cluster are calculated.

Value var a list of matrices of squared loadings i.e. for each cluster of variables, the squared loadings on first principal component of PCAmix. For quantitative variables (resp. qualitative), squared loadings are the squared correlations (resp. the correlation ratios) with the first PC (the cluster center).

sim a list of matrices of similarities i.e. for each cluster, similarities between their variables. The similarity between two variables is defined as a square cosine: the square of the Pearson correlation when the two variables are quantitative; the correlation ratio when one variable is quantitative and the other one is qualitative; the square of the canonical correlation between two sets of dummy variables, when the two variables are qualitative. sim is 'NULL' if 'matsim' is 'FALSE'.

cluster a vector of integers indicating the cluster to which each variable is allocated.

wss the within-cluster sum of squares for each cluster: the sum of the correlation ratio (for qualitative variables) and the squared correlation (for quantitative variables) between the variables and the center of the cluster.

E

  into groups of variables Description Cuts a hierarchical tree of variables resulting from hclustvar into several clusters by specifying the desired number of clusters. Usage cutreevar(obj, k = NULL, matsim = FALSE) Arguments obj an object of class 'hclustvar'. k an integer scalar with the desired number of clusters.

  fr>, Vanessa Kuentz, Benoit Liquet, Jerome Saracco

	hclustvar
	Examples
	data(decathlon)
	tree <-hclustvar(decathlon[,1:1 ])
	plot(tree)
	#choice of the number of clusters
	stability(tree,B=4 )
	part <-cutreevar(tree,4)
	print(part)
	summary(part)
	See Also
	hclustvar, summary.clustvar,stability

  If j is positive then the merge was with the cluster formed at the (earlier) stage j of the algorithm. Thus negative entries in merge indicate agglomerations of singletons, and positive entries indicate agglomerations of non-singletons. Marie Chavent <marie.chavent@u-bordeaux2.fr>, Vanessa Kuentz, Benoit Liquet, Jerome Saracco

	height clusmat merge See Also cutreevar, plot.hclustvar, stability,kmeansvar a set of p-1 non-decreasing real values: the values of the aggregation criterion. a p by p matrix with group memberships where each column k corresponds to the elements of the partition in k clusters. a p-1 by 2 matrix. Row i of merge describes the merging of clusters at step i of the clustering. If an element j in the row is negative, then observation -j was Examples #quantitative variables data(decathlon) tree <-hclustvar(decathlon[,1:1 ]) plot(tree) #qualitative variables with missing values data(vnf) tree_NA <-hclustvar(X.quali=vnf) plot(tree_NA) dev.new() vnf2<-na.omit(vnf) tree <-hclustvar(X.quali=vnf2) merged at this stage. Author(s) plot(tree)

#mixture of quantitative and qualitative variables data(wine) X.quanti <-wine[,c(3:29)] X.quali <-wine[,c(1,2)] tree <-hclustvar(X.quanti,X.quali) plot(tree) kmeansvar k-means clustering of variables

  the pourcentage of homogeneity which is accounted by the partition in k clusters. An object of class clustvar generated by the functions cutreevar and kmeansvar.DescriptionEvaluates the stability of partitions obtained from a hierarchy of p variables. This hierarchy is performed with hclustvar and the stability of the partitions of 2 to p-1 clusters is evaluated with a bootstrap approach. The boostrap approch is the following: hclustvar is applied to B boostrap samples of the n rows. The partitions of 2 to p-1 clusters obtained from the B bootstrap hierarchies are compared with the partitions from the initial hierarchy . The mean of the corrected Rand indices is plotted according to the number of clusters. This graphical representation helps in the determination of a suitable numbers of clusters.ValueReturns a list of matrices of squared loadings i.e. for each cluster of variables, the squared loadings on first principal component of PCAmix. For quantitative variables (resp. qualitative), squared loadings are the squared correlations (resp. the correlation ratios) with the first PC (the cluster center). If the partition of variables has been obtained with kmeansvar the number of iteration until convergence is also indicated. Marie Chavent <marie.chavent@u-bordeaux2.fr>, Vanessa Kuentz, Benoit Liquet, Jerome Saracco

	summary.clustab	plot.hclustvar print.hclustvar stability summary.clustvar
	Author(s) Arguments
	Marie Chavent <marie.chavent@u-bordeaux2.fr>, Vanessa Kuentz, Benoit Liquet, Jerome Saracco plot.clustab Plot of an index of stability of partitions of variables x an object of class hclustvar. sub a sub title for the plot. print.clustvar Print a 'clustvar' object stability Stability of partitions from a hierarchy of variables summary.clustab Summary of a 'clustab' object
	See Also summary.clustvar,print.clustvar,stability,cutreevar Author(s) Description Description Description ... further arguments passed to or from other methods. Index
	Examples Plot of the index of stability of the partitions against the number of clusters. Marie Chavent <marie.chavent@u-bordeaux2.fr>, Vanessa Kuentz, Benoit Liquet, Jerome Saracco This is a method for the function print for objects of the class clustvar. This is a method for the function summary for objects of the class clustab. Author(s)
	data(decathlon) #choice of the number of clusters tree <-hclustvar(X.quanti=decathlon[,1:1 ]) stab <-stability(tree,B=6 ) #a random set of variables is chosen as the initial cluster centers, nstart=1 times * Topic cluster rand, 11 Usage ## S3 method for class 'clustab' plot(x, nmin = NULL, nmax = NULL, ...) See Also hclustvar Usage ## S3 method for class 'clustvar' cutreevar, 2 Usage ## S3 method for class 'clustab' stability, 2, 3, 5, 7-9, 11, 12, 13 hclustvar, 4 See Also kmeansvar, 5 summary.clustab, 13 print(x, ...) summary(object, ...) kmeansvar, cutreevar mixedVarSim, 7 summary.clustvar, 3, 7, 13
	part1 <-kmeansvar(X.quanti=decathlon[,1:1 ],init=5,nstart=1 ) summary(part1) #the partition from the hierarchical clustering is chosen as initial partition part_init<-cutreevar(tree,5)$cluster part2<-kmeansvar(X.quanti=decathlon[,1:1 ],init=part_init,matsim=TRUE) summary(part2) part2$sim mixedVarSim Similarity between two variables Examples plot.clustab, 8 Arguments x an object of class clusab. nmin the minimum number of clusters in the plot. nmax the maximum number of clusters in the plot. See Also stability Examples data(decathlon) tree <-hclustvar(X.quanti=decathlon[,1:1 ]) data(wine) X.quanti <-wine[,c(3:29)] Arguments x Arguments plot.hclustvar, Examples rand, 11 Usage object An object of class clustab generated by the function stability. data(decathlon) stability, 12 X.quali <-wine[,c(1,2)] tree <-hclustvar(X.quanti,X.quali) plot(tree) print.clustab Print a 'clustab' object Description This is a method for the function print for objects of the class clustab. Usage ## S3 method for class 'clustab' ... Further arguments to be passed to or from other methods. They are ignored in this function. Author(s) Marie Chavent <marie.chavent@u-bordeaux2.fr>, Vanessa Kuentz, Benoit Liquet, Jerome Saracco See Also cutreevar , kmeansvar print.hclustvar Print a 'hclustvar' object stability(tree, B = 1 , graph = TRUE) Arguments tree an object of class hclustvar. B the number of bootstrap samples. graph boolean, if 'TRUE' a graph is displayed. Value matCR matrix of corrected Rand indices. meanCR vector of mean corrected Rand indices. ... Further arguments passed to or from other methods. Author(s) Marie Chavent <marie.chavent@u-bordeaux2.fr>, Vanessa Kuentz, Benoit Liquet, Jerome Saracco See Also stability summary.clustvar part<-kmeansvar(X.quanti=decathlon[,1:1 ],init=5) summary(part) * Topic multivariate cutreevar, 2 hclustvar, 4 kmeansvar, 5 mixedVarSim, 7 plot.clustab, 8 plot.hclustvar, rand, 11 stability, 12 * Topic package ClustOfVar-package, 2 * Topic print print.clustab, print.clustvar, Summary of a 'hclustvar' object print.hclustvar, 10	CemOA : archive ouverte d'Irstea / Cemagref CemOA : archive ouverte d'Irstea / Cemagref
	stab<-stability(tree,B=2 ) print(x, ...) Author(s)
	plot(stab,nmax=7) ClustOfVar (ClustOfVar-package), 2
	Usage mixedVarSim(X1, X2) plot.hclustvar Arguments x An object of class clustab generated by the function stability. Dendrogram of the hierarchy of variables ... Further arguments to be passed to or from other methods. They are ignored in Description Marie Chavent <marie.chavent@u-bordeaux2.fr>, Vanessa Kuentz, Benoit Liquet, Jerome Saracco ClustOfVar-package, Description cutreevar, 2, 5, 7, 10, This is a method for the function print for objects of the class hclustvar. This is a method for the function summary for objects of the class clustvar. See Also hclustvar, 2, 3, 4, 9-12
	size Arguments X1 X2 Description Dendrogram of the hierarchy of variables resulting from hclustvar. a vector or a factor a vector or a factor this function. Author(s) Usage plot.clustab, hclustvar Usage kmeansvar, 2, 5, 5, 10, ## S3 method for class 'hclustvar' print(x, ...) ## S3 method for class 'clustvar' Examples mixedVarSim, 7 summary(object, ...) Marie Chavent <marie.chavent@u-bordeaux2.fr>, Vanessa Kuentz, Benoit Liquet, Jerome Saracco data(decathlon) the number of variables in each cluster. scores a n by k numerical matrix which contains the k cluster centers. The center of a cluster is a synthetic variable: the first principal component calculated by Author(s) Marie Chavent <marie.chavent@u-bordeaux2.fr>, Vanessa Kuentz, Benoit Liquet, Jerome Saracco Usage ## S3 method for class 'hclustvar' See Also stability Arguments x plot.clustab, 8, 12 tree <-hclustvar(X.quanti=decathlon[,1:1 ]) stab<-stability(tree,B=2 ) Arguments plot.hclustvar, 5, 8 An object of class hclustvar generated by the function hclustvar. ... Further arguments to be passed to or from other methods. They are ignored in plot(stab,nmax=7) object print.clustab, 9 an object of class clustvar. print.clustvar, 7,
	plot(x, sub="", ...) this function. ... further arguments passed to or from other methods. print.hclustvar, 10

PCAmix. The k columns of scores contain the scores of the n observations units on the first PCs of the k clusters.

Description

Returns the similarity between two quantitative variables, two qualitative variables or a quantitative variable and a qualitative variable. The similarity between two variables is defined as a square cosine: the square of the Pearson correlation when the two variables are quantitative; the correlation ratio when one variable is quantitative and the other one is qualitative; the square of the canonical correlation between two sets of dummy variables, when the two variables are qualitative.

... further arguments passed to or from other methods.

Author(s)

Marie Chavent <marie.chavent@u-bordeaux2.fr>, Vanessa Kuentz, Benoit Liquet, Jerome Saracco

CemOA : archive ouverte d'Irstea / Cemagref

CemOA : archive ouverte d'Irstea / Cemagref print.clustab


Hani Abdallah
email: hani.abdallah@teledetection.fr

Jean-Stéphane Bailly

Nicolas Baghdadi

Nicolas Lemarquand

Improving the assessment of ICESat water altimetry accuracy accounting for autocorrelation

1 Given that water resources are scarce and are strained by competing demands, it has become crucial to develop and improve techniques to observe the temporal and spatial variations in the inland water volume. Due to the lack of data and the heterogeneity of water level stations, remote sensing, and especially altimetry from space, appear as complementary techniques for water level monitoring. In addition to spatial resolution and sampling rates in space or time, one of the most relevant criteria for satellite altimetry on inland water is the accuracy of the elevation data. Here, the accuracy of ICESat LIDAR altimetry product is assessed over the Great Lakes in North America.

The accuracy assessment method used in this paper emphasizes on autocorrelation in high temporal frequency ICESat measurements. It also considers uncertainties resulting from both in situ lake level reference data.

abilistic upscaling process was developed. This process is based on several successive ICESat shots averaged in a spatial transect accounting for autocorrelation between successive shots. The method also applies pre-processing of the ICESat data with saturation correction of ICESat waveforms, spatial filtering to avoid measurement disturbance from the land-water transition effects on waveform saturation and data selection to avoid trends in water elevations across space. Initially this paper analyzes 237 collected ICESat transects, consistent with the available hydrometric ground stations for four of the Great Lakes. By adapting a geostatistical framework, a high frequency autocorrelation between successive shot elevation values was observed and then modeled for 45% of the 237 transects. The modeled autocorrelation was therefore used to estimate water elevations at the transect scale and the resulting uncertainty for the 117 transects without trend. This uncertainty was 8 times greater than the usual computed uncertainty, when no temporal correlation is taken into account. This temporal correlation, corresponding to approximately 11 consecutive ICESat shots, could be linked to low transmitted ICESat GLAS energy and to poor weather conditions. Assuming Gaussian uncertainties for both reference data and ICESat data upscaled at the transect scale, we derived GLAS deviations statistics by averaging the results at station and lake scales. An overall bias of -4.6 cm (underestimation) and an overall standard deviation of 11.6 cm were computed for all lakes. Results demonstrated the relevance of taking autocorrelation into account in satellite data uncertainty assesment.

Keywords: GLAS, LIDAR, accuracy, temporal correlation, block kriging, Great Lakes 2

Introduction

Although lakes and rivers correspond to only 0.27% of the global fresh water and 0.007% of the Earth's water budget, they constitute the most accessible inland water resources available for ecosystems and human consumption.

Given that these resources are scarce and are strained by competing demands, it has become crucial to develop and improve techniques to observe the temporal and spatial variations in the water volume of lakes, rivers and wetlands to meet human needs and assess the ongoing impacts of climate change.

For this reason, most countries operate a national network of inland water level stations to collect information for water resource development.

The installation and maintenance of such networks are expensive tasks that only international cooperation programs can initiate and sustain in many countries. The availability and access to data from water level stations are therefore severely limited due to a decline in the number of stations, inadequate monitoring networks, gaps in records, differences in processing and quality control, differences in datum level and differences in data policies [START_REF] Harvey | Expert meeting on hydrological data for global studies[END_REF][START_REF] Chen | Evolutionary artificial neural networks for hydrological systems forecasting[END_REF]. As a consequence, water resource sampling is neither spatially nor temporally homogeneous.

Due to the lack of data and the heterogeneity of water level stations, remote sensing, and especially altimetry from space, appear as complementary techniques for water level monitoring [START_REF] Calmant | Continental surface waters from satellite altimetry[END_REF].

Satellite radar altimetry first appeared in 1974 with Skylab. In 1975, the GEOS-3 radar altimeter was designed to monitor ocean surfaces, followed by the SEASAT (1978) and GEOSAT (1985GEOSAT ( -1989) ) missions. Radar altimetry entered a pre-operational phase in 1992 with the satellites ERS and Topex-Poseidon, followed by Envisat and Jason in 2001. All these missions were originally designed for measuring the level of the ocean through a combination of a radar technique determining the distance from the satellite to the reflecting surface and a satellite positioning technology identifying the precise location (within centimeters) of the satellite. Since 1990, the radar elevation measurements have demonstrated their relevance for inland water observation [START_REF] Bercher | Précision de l'altimetrie satelittaire radar sur les cours d'eau[END_REF]. The application of this technique in the literature has allowed monitoring of the inland seas [START_REF] Aladin | Modern hydro-biological state of the Small Aral sea[END_REF], lakes [START_REF] Birkett | The contribution of TOPEX/POSEIDON to the global monitoring of climatically sensitive lakes[END_REF][START_REF] Birkett | Synergistic remote sensing of Lake Chad: Variability of basin inundation[END_REF][START_REF] Birkett | Surface water dynamics in the Amazon Basin: Application of satellite radar altimetry[END_REF][START_REF] Crétaux | Lake studies from satellite radar altimetry[END_REF] and large rivers [START_REF] Birkett | Contribution of the TOPEX NASA radar altimeter to the global monitoring of large rivers and wetlands[END_REF][START_REF] Mercier | Alimétrie spatiale sur les eaux continentales: apport des missions TOPEX/POSEIDON et ERS1&2 à l'étude des lacs[END_REF][START_REF] Maheu | Water level fluctuations in the Plata Basin (South America) from Topex/Poseidon satellite altimetry[END_REF].

In addition to spatial resolution and sampling rates in space or time, the most relevant criterion for satellite altimetry on inland water is the accuracy of the elevation data; therefore, many studies aimed to quantify these accuracies. [START_REF] Morris | Evaluation of the TOPEX/POSEIDON altimeter system over the Great Lakes[END_REF] found a Topex-Poseidon data RMSE (Root Mean Square Error) of a few centimeters on the Great Lakes of North America. Studies over other world lakes indicated the following RMSE values: i) greater than 10 cm for Lake Chad [START_REF] Birkett | Synergistic remote sensing of Lake Chad: Variability of basin inundation[END_REF]; ii) approximately 5 cm for Issyk Kul Lake, Kyrgyzstan; and iii) 10 cm for Chardarya Lake, Kazakhstan [START_REF] Crétaux | Lake studies from satellite radar altimetry[END_REF]. The accuracy of radar altimetry was also studied for a selection of rivers and floodplains in the Amazon Basin. Substantial errors were observed with RMSE values between 40 cm and 1.1 m [START_REF] Birkett | Surface water dynamics in the Amazon Basin: Application of satellite radar altimetry[END_REF][START_REF] Bercher | Précision de l'altimetrie satelittaire radar sur les cours d'eau[END_REF].

In recent years, LIDAR onboard satellite also appears as a promising tool for accurate, high resolution altimetry. The Geoscience Laser Altimeter System (GLAS) ranging instrument onboard the ICESat (Ice Cloud and land Elevation Satellite) provides elevation data over all Earth surfaces. Its main 4

Author-produced version of the article published in ISPRS Journal of Photogrammetry and Remote Sensing, 2011, 66(6), 833-844. The original publication is available at http://www.isprs.org/ DOI : 10.1016/j.isprsjprs.2011.09.002 strength is its small footprint diameter averaged over 70 m [START_REF] Zwally | ICESat's 28[END_REF] compared to the larger radar footprint size of 250 m to a few km. This relative small footprint size is promising for inland water monitoring, especially in rivers. ICESat was launched on the 13 January 2003 and the mission stopped at 14 august 2010. It included three lasers that transmitted short pulses (4-6 ns) of infrared light (1064 nm) and visible green light (532 nm). GLAS used 1064-nm laser pulses to measure the heights of the surface and dense cloud and 532 nm pulses to measure the vertical distribution of clouds and aerosols. The 1064-nm signal was also separately filtered and digitized at 2 MHz for detection of dense clouds and aerosols at 76.8-m vertical resolution. The three lasers have been operated one at a time, sequentially throughout the mission. To extend mission life, the operational mode included 33-day to 56-day campaigns, several times per year. Each period has been assigned a campaign or operations period identifier, such as Laser 2a, to denote the operating laser (2) and the operations period (a) [START_REF] Schutz | Overview of the ICESat mission[END_REF]. Laser pulses or shots at 40 Hz illuminated sites spaced at 170-meter intervals along track over Earth's surface. The primary mission of GLAS was to measure changes in ice sheet elevation, and secondary objectives included the measurement of cloud and aerosol height profiles, land elevation, vegetation cover and sea ice thickness [START_REF] Zwally | ICESat's 28[END_REF][START_REF] Schutz | Overview of the ICESat mission[END_REF]. Although less devoted to water monitoring compared to radar altimetry, these LIDAR data could also be used to observe the inland water altimetry [START_REF] Urban | A survey of ICE-Sat coastal altimetry applications: continental coast, open ocean island, and inland river[END_REF].

To date, few studies assessed GLAS elevation data accuracy for inland waters. The estimated standard deviation of GLAS elevation data computed on Lake Nasser in southern Egypt range from 3 to 8 cm [START_REF] Chipman | Satellite-based assessment of the dynamics of new lakes in southern Egypt[END_REF]. On Otter Tail County lakes (USA), [START_REF] Bhang | Verification of the vertical error in C-band SRTM DEM using ICESat and Landsat-7, Otter Tail County, MN[END_REF] 5

Author-produced version of the article published in ISPRS Journal of Photogrammetry and Remote Sensing, 2011, 66(6), 833-844. The original publication is available at http://www.isprs.org/ DOI : 10. 1016/j.isprsjprs.2011.09.002 indicates differences between GLAS elevation and hydrologic station elevation ranging from 2 to 35 cm. [START_REF] Urban | A survey of ICE-Sat coastal altimetry applications: continental coast, open ocean island, and inland river[END_REF] reports GLAS elevation RMSE of approximately 3 cm under clear conditions, 8-15 cm under cloudy skies and 25 cm under very cloudy skies on the Tapajos Rivers (Brazil).

Several problems were revealed in previous satellite altimetry accuracy studies. First, detailed descriptions of the assumptions and accuracy statistics computation are lacking. Second, most of the studies compute errors between satellite and ground data considering ground data as 'truth', or exact and free from uncertainty. Consequently, in terms of error distribution estimation, most studies only compute RMSE as an accuracy statistic, but despite its common use, RMSE is inadequate because it mixes the bias (exactness) and dispersion (standard deviation) of satellite measurements, which can have very different consequences on their final utility. Moreover, no study compares a single shot elevation value to the ground reference data; rather, studies compare reference data to an average of successive shots, consistent with the bodies of water under consideration, such as a river section or a continuous lake transect. Therefore, before accuracy statistics computation, an implicit aggregating or upscaling process, which is typically an average based on the arithmetic mean, involves successive shot elevation values that occurs at the scale of the studied water body. As a result from this aggregating process, the satellite 'measurement' is thus compared to the corresponding ground data. The use of an arithmetic average assumes that successive shots are considered as independent [START_REF] Morris | Evaluation of the TOPEX/POSEIDON altimeter system over the Great Lakes[END_REF] and identically distributed, allowing the usual statistics computation on measurement uncertainties, such as the SDOM (Standard Deviation of the Mean). The centimeter differences observed between satellite measurements and ground measurements and the uncertainty of both measurements can make 6

Author-produced version of the article published in ISPRS Journal of Photogrammetry and Remote Sensing, 2011, 66(6), 833-844. The original publication is available at http://www.isprs.org/ DOI : 10.1016/j.isprsjprs.2011.09.002 the error word inappropriate to qualify the observed deviations between satellite and ground measurements. If not, how to estimate the bias (mean) and dispersion of these deviations ? To answer this question, we need to first consider the satellite and ground measurements as random variables and to verify their significant deviations by computing the uncertainty of each measurement, ground measurements and satellite measurements, at water body scale. The uncertainty of satellite measurements at water body scale can be estimated using the SDOM only if the shot elevation values are, for high frequencies, independent in time and space. If not, the correlation between successive shot altimetrical values, denoted further as autocorrelation, must be modeled and estimated. Uncertainty at the water body scale must be computed accordingly. Although infrequent in satellite elevation data, the support effect accounts for the dependence between sequential measurements in upscaling and is a typical practice used on spatial data [START_REF] Chiles | Geostatistics: modeling spatial uncertainty[END_REF] [START_REF] Atkinson | Spatial scale problems and geostatistical solutions: A review[END_REF]. It has already been used for remote sensing data accuracy assessments [START_REF] Crosson | Impacts of spatial scaling errors on soil moisture retrieval accuracy at L-band[END_REF] and even on LIDAR data [START_REF] Bailly | Geostatistical estimation of bathymetric LIDAR errors on rivers[END_REF].

This paper aims to i) explore the high frequencies autocorrelation in GLAS-ICESat water elevation shot data, ii) assess the consequences of this autocorrelation in the uncertainty statistics on water level estimates and iii) assess the consequences of this autocorrelation in deviations to ground station elevation data from the Great Lakes. To ensure a proper comparison, ICESat data pre-processing was first performed. It consisted of i) the exclusion of shots with waveform saturation corresponding to data near land [START_REF] Baghdadi | The relevance of GLAS/ICESat elevation data for the monitoring of river networks[END_REF], ii) the exclusion of shot outliers from atmospheric disturbance, iii) the application of saturation corrections to initial elevations [START_REF] Urban | A survey of ICE-Sat coastal altimetry applications: continental coast, open ocean island, and inland river[END_REF] and iiii) the exclusion of data when lake water surface 7

Author-produced version of the article published in ISPRS Journal of Photogrammetry and Remote Sensing, 2011, 66(6), 833-844. The original publication is available at http://www.isprs.org/ DOI : 10.1016/j.isprsjprs.2011.09.002 can not be considered as flat for accuracy and deviations estimates. The autocorrelation between successive shot elevation values was first tested and modeled for each of the considered sets of shots. A set of shots consist of a GLAS transect within a Great Lake, spatially coherent with a given ground water-level station (a water body) and temporally located in the same satellite track. We thus proposed an explicative model of the autocorrelation presence. In cases of significant autocorrelation and absence of trend in transect shot data, a block kriging process [START_REF] Wackernagel | Multivariate Geostatistics[END_REF] was used to estimate a mean water elevation at the transect scale and its uncertainty, Finally, the ground measurements and satellite measurements and their relative uncertainties at the transect scale were estimated. By considering these two water elevation estimates as random variables, according deviations statistics (bias and standard deviations) were proposed.

Study site and data set description

Study site

Due to the lack of data on one of the five Great Lakes in North America, the assessment of GLAS-ICESat elevation data accuracy was assessed for the following lakes: Superior, Michigan, Erie and Ontario. One of the youngest natural features on the North American continent, the Great Lakes make up the largest surface freshwater system on Earth. Covering more than 94,000 square miles and draining more than twice as much land, these freshwater seas hold an estimated 22.7 trilliards cubic meters of water, approximately one-fifth of the world's surface freshwater and nine-tenths of the U.S. supply.

A network of 52 hydrometric stations monitor the Great Lakes (25 Canadian MEDS station and 27 US NOAA stations) for level measurement (Fig. 1).

8

Author-produced version of the article published in ISPRS Journal of Photogrammetry and Remote Sensing, 2011, 66(6), 833-844. The original publication is available at http://www.isprs.org/ DOI : 10.1016/j.isprsjprs.2011.09.002

The Great Lakes have from 1 to 4 cm tides but no prevalent ocean tides.

Thus, the Great Lakes are an excellent site for accuracy studies.

Here is located Fig. 1 2.2. Data set

ICESat data

The GLAS elevation data used in this study are the Level-2 altimetry products GLA14 and GLA01 (Release 31) that provide surface elevation measurements (land, water, etc.) and corresponding waveforms. GLA14 data also include the laser footprint geolocation with a precision smaller than 4 m [START_REF] Duong | Full waveform analysis: Icesat laser data for land cover classification[END_REF]. GLA01 data include waveforms that have been decomposed into multiple Gaussian distributions [START_REF] Wagner | Gaussian decomposition and calibration of a novel small-footprint full-waveform digitising airborne laser scanner[END_REF][START_REF] Jutzi | Range determination with waveform recording laser systems using a Wiener Filter[END_REF][START_REF] Chen | Assessment of terrain elevation derived from satellite laser altimetry over mountainous forest areas using airborne LIDAR data[END_REF] corresponding to 544 or 1000 samples of received power in volts at 1 ns sampling rate (see [START_REF] Zwally | ICESat's 28[END_REF] for a detailed description of ICESat full waveform data). These data record the period between 20 February 2003 and 11 October 2009. In this study, 237 ICESat transects, with lengths ranging from 5 to 20 km, containing 20,224 shots have been used (Table 1). These 237 selected transects are the set of continuous tracks of shots that are i) in a radius of 20 km from an available hydrometric station to avoid the influence of natural spatial lake level variations and ii) at least 2 km from the shore to avoid measurement disturbance from land-water transition effects on waveform saturation [START_REF] Baghdadi | The relevance of GLAS/ICESat elevation data for the monitoring of river networks[END_REF].

Next we developed a two step GLA14 data cleaning procedure, close to the one proposed by [START_REF] Zhang | Monitoring lake level changes on the Tibetan Plateau using ICESat altimetry data (2003-2009)[END_REF]. In step 1, each of the 237 GLAS 9

Author-produced version of the article published in ISPRS Journal of Photogrammetry and Remote Sensing, 2011, 66(6), 833-844. The original publication is available at http://www.isprs.org/ DOI : 10.1016/j.isprsjprs.2011.09.002

ICESat transects was first visualy inspected on time-elevation bivariate plots (Figs. 5-a and 5-b) in order to remove obvious outliers and to verify that the water level was almost constant. Some parts of the transects were accordingly excluded when the mean water level varied or when observations were 'constant' due to i) cloudy episodes with outlier elevations in transect (indicating a difference in elevation of approximately 1000 m) [START_REF] Wesche | Precise surface topography of the grounded ice ridges at the Ekströmisen, Antarctica, based on several geophysical data sets[END_REF], ii) the proximity of river outlets. Fig. 2 shows the only example of this phenomenon: the extended transect of 14 June 2005 for Station 21 (Holland station) on Lake Michigan is located near an important river outlet (Fig. 2-a). The transect profile shows the 85 km elevation profile track that increases in the middle near the outlet, producing an upper water level of approximately 70 cm. In Step 2, for each transect, the median elevation value of the remaining shots was calculated and the shots which were at least 4 meters far from this median value were removed.

Here is located Fig. 2 2

.2.2. Water level reference data

The reference lake level data in this study were obtained from 27 of the 52 water level stations located throughout the Great Lakes basin. The selected stations are near (less than 20 km) the transect of the GLAS footprints.

Eleven stations are operated by the Marine and Environmental Data Service (MEDS) of Canada's Department of Fisheries and Oceans, and 16 stations are operated by the U.S. National Oceanic and Atmospheric Administration (NOAA). Both NOAA and MEDS stations continuously monitor the surface levels of the four lakes at 6 minutes intervals [START_REF] Braun | ICESat laser altimetry in the Great Lakes[END_REF]) with a 10

Author-produced version of the article published in ISPRS Journal of Photogrammetry and Remote Sensing, 2011, 66(6), 833-844. The original publication is available at http://www.isprs.org/ DOI : 10.1016/j.isprsjprs.2011.09.002 standard deviation of 2 cm [START_REF] Hovis | Effects of hurricane Isabel on water levels[END_REF].

North American Great Lakes are known to have low tides but are sensitive to short-term water-level fluctuations, such as up to 2.44 m in approximately 2 hours, caused by wind or storm surges, which are known as seiches [START_REF] Touchart | Limnologie physique et dynamique : Une géographie des lacs et des étangs[END_REF]. Considering this phenomenon and because the ICESat transect lasts very few seconds and a hydrometric station may have a timing bias of a few minutes, this study analyzed the dispersion of recorded ground measurement (hydrometric station data) in a 6-minute range around transect ICESat times for all stations. This analysis exhibits a constant dispersion for each station of 7 mm for one standard deviation.

In addition, the geoid model G99SSS used in this study to convert the reference data from ellipsoidal height (WGS84) to orthometric height (NAVD88) has a mean error of approximately 2.5 cm [START_REF] Smith | GEOID99 and G99SSS: 1-arc-minute geoid models for the United States[END_REF].

Given that these different sources of uncertainty are independent, they can be summed (in variance). As a result, the computed uncertainty on reference ground data was estimated to 3.3 cm for one standard deviation (σ R ).

Data consistency

GLA14 products provide original elevation data with ellipsoidal heights based on the same ellipsoid used by the Topex satellite. The lake level data collected from hydrometric stations are in the International Great Lakes Datum of 1985 (IGLD85). Consequently, data transformations are required in order to conduct a consistent comparison in the same reference frame.

GLAS heights and lake level heights were therefore converted to the same vertical datum NAVD88. Two successive transformations are required for GLAS heights: the first converts the Topex ellipsoidal heights to WGS84 11

Author-produced version of the article published in ISPRS Journal of Photogrammetry and Remote Sensing, 2011, 66(6), 833-844. The original publication is available at http://www.isprs.org/ DOI : 10.1016/j.isprsjprs.2011.09.002 ellipsoidal heigts using the geoid model EGM96; the second transforms the WGS84 ellipsoidal heights to NAVD88 orthometric heights using the geodic model G99SSS. For hydrometric station elevation data, the IGLD85 lake level datum converts to orthometric height NAVD88.

Methods

GLAS water elevation saturation correction

Some ICESat waveforms exhibit distortions, or saturation events, when the high energy emitted by the laser returns, passes the inadequate automatic gain controls and overloads the detector [START_REF] Sun | ICESat receiver signal dynamic range assessment and correction of range bias due to saturation[END_REF]. This occurs when the amplitude of the energy returned by a number of 1-ns bins is greater than the threshold function of gain [START_REF] Urban | A survey of ICE-Sat coastal altimetry applications: continental coast, open ocean island, and inland river[END_REF][START_REF] Abshire | Geoscience Laser Altimeter System (GLAS) on the ICE-Sat mission: On-orbit measurement performance[END_REF], producing an incorrect elevation. The ICESat science team has developed an analytical saturation correction, also available on the data products. Adding this correction to the elevation is recommended for ice and calm water data, but not over oceans [START_REF] Urban | A survey of ICE-Sat coastal altimetry applications: continental coast, open ocean island, and inland river[END_REF]. Corrections vary between 0 and 1.5 m. For some waveforms, the saturation should be corrected but the method used to calculate the correction coefficient cannot provide a correct value and is marked by '-999.00' flags. Fig. 3 shows some ICESat waveforms obtained across Lake Superior; the saturated waveform is clipped and widened artificially by the returning high energy.

Here is located Fig. 3 Fig. 4 illustrates an example of elevation transects before and after application of the saturation correction. The saturation correction allows a reduction of error bias but does not reduce the error dispersion in the figure .     Here is located Fig. 4 3.2. GLAS water elevation estimation at transect scale To compute the GLAS water elevation estimate ẐG and uncertainty σ G at the transect scale and considering autocorrelation, we used a statistical framework based on the regional variable and random field theories [START_REF] Cressie | Statistics for spatial data[END_REF]. In this framework, a GLAS transect is viewed as a temporal block, or a delimited part of a temporal random field (z(t)). For each GLAS transect, we developed the usual block kriging process (Wackernagel, 1995, p.77) based on an autocorrelation model, allowing both the estimation of the water elevation ẐG and its corresponding uncertainty σ G , at the transect scale.

Autocorrelation test and modeling

A GLAS transect is composed from (t 1 . . . t n ) and (z 1 . . . z n ) series, respectively, where t refers to the laser shot (pulse) time, z to the elevation in the NAVD88 altimetrical system and n is the number of shots in the transect. These series are denoted as z(t) and are considered to be a realization of a stationary temporal random function. We want to infer the main property of z(t), its temporal correlation (autocorrelation), especially the high frequency autocorrelation. When z(t) is stationary, with constant expectation and homogeneous temporal covariance, this inference is possible using the variogram modeling of the field z(t), equivalent to the covariance.

Elevation trend test at transect scale. To ensure a stationary process and reenforce the high frequency autocorrelation test power [START_REF] Armstrong | Problems with universal kriging[END_REF], 13

Author-produced version of the article published in ISPRS Journal of Photogrammetry and Remote Sensing, 2011, 66(6), 833-844. The original publication is available at http://www.isprs.org/ DOI : 10.1016/j.isprsjprs.2011.09.002 a linear trend test is first performed. A linear model z(t) = m(t)+y(t), where m(t) = at + b is the linear regression of z by t, and y(t) are the residuals of this linear regression, is fitted. Using the usual T-test on the coefficient a it is tested if a linear trend is significant.

Next, when the trend tests positively, the initial z(t) field is decomposed in two additive term m(t) + y(t), with m(t) a deterministic trend term and y(t) a term consisting of random residuals. When the trend is negative, the m(t) term disappears and we write z(t) = y(t) for the sake of simplicity. To enable a valid comparison between ground station and GLAS measurements and to avoid natural spatial variations of lake water levels at high spatial frequencies, only transects without linear trends were kept further to compute the GLAS deviations and accuracy.

Autocorrelation significance test. Next, to estimate the autocorrelation of z(t) transect data, an experimental variogram γ(h), with h equal to lag time, is computed on the y(t) data (see Wackernagel (1995, p.35) for the usual experimental variogram formulation). This variogram was computed on 16 regular lag times each 62.5 ms ranging from 0 to 1 s, chosen to obtain numerous shot pairs for an accurate variogram estimation for each lag time h. An autocorrelation is considered significant when the experimental variogram graph shows a clear increasing shape from small time differences to large time differences. When the corresponding experimental variogram graph shows a flat variogram, no autocorrelation for this GLAS transect exists. To test the significance of the autocorrelation calculated from the experimental variogram, we used the standard empirical Mantel test [START_REF] Legendre | Spatial pattern and ecological analysis[END_REF], which simulates data corresponding to the null hypothesis H 0 through data permutations, where H 0 corresponds to the absence of 14

Author-produced version of the article published in ISPRS Journal of Photogrammetry and Remote Sensing, 2011, 66(6), 833-844. The original publication is available at http://www.isprs.org/ DOI : 10.1016/j.isprsjprs.2011.09.002 autocorrelation (Fig. 5-c). At each permutation of the z data, a corresponding experimental variogram is computed. For numerous replications of the permutations, we obtained a quintile distribution of variograms representing H 0 for each lag time h. The H 0 acceptance area is thus represented in the variogram graph through a 95% confidence band of H 0 . In the case of significant autocorrelation, the actual experimental variogram for small values of h is below the 2.5% quintile of H 0 (the lower line of the band). Here, we chose to reject H 0 only when the first variogram point was below the 2.5% quintile of H 0 .

Here is located Fig. 5 Autocorrelation modeling. In the case of a significant experimental variogram (Fig. 5-d), we chose to model it with a single spherical model function given by:

γ(h) =      nu + (si -nu)( 3h 2r -h 3 2r 3 ) if h < r (nu + si) if h ≥ r (1)
Using a weighted least-square estimation [START_REF] Pebesma | Gstat: A program for geostatistical modelling, prediction and simulation[END_REF], the model parameters nugget (nu), range (r) and sill (si) were fitted. This was done automatically for all GLAS transects from initial fitting parameters corresponding to a nugget equal to the experimental variogram at a smaller lag time, a sill equal to the maximum value of the experimental variogram and a range equal to the time duration at which the variogram reaches its maxium. After the fitting process, a variogram model γ(h) was obtained as shown in the example by the fitted line in Fig. 5-d.

Transect water elevation estimation with autocorrelation: block kriging

For a GLAS transect, when autocorrelation is significant and modeled, a mean value ẐG for the transect is estimated from transect data z(t) from Ŷ . The Ŷ estimate is computed through a block kriging process of y(t) (Wackernagel, 1995, p.77) using the variogram model γ(h). Ŷ is therefore a linear combination of y(t) data with kriging weights β i :

Ŷ = n i=1 β i y i (2)
Uncertainty on ẐG is represented in the block kriging estimation standard deviation σ G given by (Wackernagel, 1995, p.77):

σ G = µ BK + γ + 2 n i=1 β i γ(t i , T ) (3)
In this equation, µ BK is the Lagrange multiplier deduced from the block kriging system, γ(t i , T ) is the average variogram computed between the sample point t i and all points of the transect T (block) and γ is the variance within the transect. In practice, both ẐG and σ G are approximated through the average of gridded points within the considered GLAS transect. We chose a gridding step of 10ms here corresponding to 1/4 of the time lag between pulses. Note that σ G only depends on the data temporal sampling level (low or high) in the lake (Fig. 6-b). For the transects that do experience significant correlation the 107 fitted variogram models show:

• a nugget nugget+sill ratio ranging from 0.24 to 0.52 with a mean value of 0.4.

This value means that there is in average 40 % of randomness between two consecutive, i.e. very close, shots.

• a range ranging from 0.27 s to 0.34 s with a mean value of 0.3 s (time the variogram model becomes flat at sill value). This range corresponds to the duration of 11 pulses and a distance along the transect of approximately 1.9 km.

Here is located Fig. 6 4

.2. Autocorrelation explanation

To explore the physical explanation of the significant autocorrelation observed for 45% of the GLAS transects, we attributed to each one of the 237 transects i) the environmental parameters from the ICESat data at the transect time (surface temperature, surface humidity and pressure) and ii) the GLA14 instrument parameters (transmitted energy, received energy, gain, incidence angle, laser campaign, transect direction and transect dates according to the beginning and end of each laser campaign [START_REF] Schutz | Overview of the ICESat mission[END_REF], table 1, p.2)). To explore relationships between the presence of autocorrelation and these available system or environmental parameters, they were all first considered as potential temporal autocorrelation predictors. We used the multivariate classification tree (CART) method of [START_REF] Breiman | Classification And Regression Tree[END_REF] that computes predictor importance that explains autocorrelation by 18

Author-produced version of the article published in ISPRS Journal of Photogrammetry and Remote Sensing, 2011, 66(6), 833-844. The original publication is available at http://www.isprs.org/ DOI : 10.1016/j.isprsjprs.2011.09.002 using indifferrently quantitative or qualitative predictors. The obtained classification tree in Fig. 7 shows separation between the autocorrelated (1) and non-autocorrelated (0) transects based respectively on the transmitted energy (TE), temperature (temp) and pressure (pressure). This tree explain 68% of the autocorrelation. These results show that the transmitted energy plays a major role in autocorrelation among all parameters. We thus draw from this analysis that significant autocorrelation occurs when the GLAS instrument transmits low energy which notably occurs at the end of each of three lasers life times, and when bad weather conditions are marked by either lower temperature or lower pressure. These results are consistent with no systematic autocorrelation observed in space or time that could have come from satellite system behaviour.

Here is located Fig. 7 4.3. Autocorrelation effect on GLAS water elevation uncertainties at transect scale Z G estimates and corresponding σ G uncertainty parameters were computed for the 117 transects where no significant trend has been observed.

When no autocorrelation is observed for the transect, we used the block kriging process (Equation 2), and if not, we used the Standard Deviation

Of the Mean (SDOM). Here is located Fig. 8 4.4. GLAS water elevation deviation distribution

The GLAS deviation, the difference between the GLAS and reference elevation, is considered as a random variable. The distributions of these deviations were computed for all Great Lakes together, for individual lakes and at station scale. For all Great Lakes, the probability of a GLAS water elevation deviation to be lower than 1 cm, 10 cm and 20 cm in absolute values is 7.5%, 66.8% and 91.9%, respectively, which seems close to a N (0, 10cm)

value. The 2.5% and 97.5% quantiles are 29 cm and 16 cm, respectively.

The expectation of this distribution is -4.6 cm.

Examples of distributions are depicted in Fig. 9, for all Great Lakes, at lake scale (for lake Erie and lake Ontario), at station scale (for station Erie, Fermi Power, Rochester and Toronto) and at transect scale (transect number 5,8,26,24,88,85,92,90). The mix of deviation distributions at the lake scale gives the overall Great Lakes deviation distribution. Normality of deviation distributions is clearly low when downscaling. At station scale, there are various shapes of deviation distribution, non-symetric, with positive or negative pseudo-biases. The mix of these deviation distributions at the station scale, with weights equal to transect numbers for each station of a lake, provides the lake deviation distribution.

20

Author-produced version of the article published in ISPRS Journal of Photogrammetry and Remote Sensing, 2011, 66(6), 833-844. The original publication is available at http://www.isprs.org/ DOI : 10.1016/j.isprsjprs.2011.09.002

However, these results show that when aggregating numerous data, i.e.

transects at the Great Lakes scale, and due to the general Law of Large Numbers, the obtained shapes of deviation distribution become nearly Gaussian. This is not true at the station scale where shapes are asymmetrical and complex. In this latter case, attention must be paid in using common accuracy statistics, especially for dispersion parameters.

Here is located Fig. 9 4.5. GLAS accuracy parameters

The right part of Table 1 summarizes the results of the biases and standard deviations statistics at the Great Lakes, lake and station scales. At the Great Lakes scale, the GLAS deviation bias is -4.6 cm. In other words, GLAS is slightly underestimating the water elevation on average. The bias values range from -11.4 cm to -2.1 cm at the lake scale and from -14.5 cm to 3.8 cm at the station scale. The bias found on the stations of the Michigan Lake is the higher in magnitude, reaching up to -14.5 cm for Station Kewaunee. A standard deviation of 11.6 cm was computed at the Great Lakes scale. It varies from 7.9 cm to 12.9 cm at the lake scale and from 3.3 cm to 13.7 cm at the station scale. The highest standard deviation was found at Station Ludington on Lake Michigan. Bias and standard deviation statistics clearly differ from one lake, or one station, to the next. All these standard deviations must be compared in magnitude with the reference data uncertainties (3.3 cm for one standard deviation) and GLAS water elevation estimate uncertainty at the transect scale (Fig. 8). as the GLAS water elevation estimate uncertainty when autocorrelation is significant.

Here is located table 1

Conclusions

This study explored the accuracy of GLAS-ICESat inland water elevation data at various spatial scales with a specific emphasis on uncertainties originating from in situ measurements and impacts of autocorrelation between successive ICESat shots, i.e. at temporal high frequencies. This study also benefited from the pre-processing of the GLAS data with the saturation correction of the GLAS waveforms and spatial filtering to avoid measurement disturbance from land-water transition effects on waveform saturation as previously observed. A set of 237 GLAS transects near the available hydrometric ground stations for four of the Great Lakes of North America was analyzed. A significant autocorrelation between successive shot elevation values was observed for 45% of the transects. This autocorrelation of an approximately 11 pulses duration seems to occur when a combination of a low transmitted energy from the GLAS instrument and bad weather with low temperature and low pressure occurs. The main consequence of this autocorrelation is a drastic increased uncertainty (by 8 times) of the GLAS water elevation at the transect scale.

After removing the 120 transects where a linear trend was observed and assuming Gaussian uncertainties for both reference data and GLAS data upscaled at the transect scale, we derived empirical distributions on GLAS deviations at the Great Lakes, lake and station scales. scale, a bias of -4.6 cm (underestimation) and a standard deviation of 11.6 cm were computed on the various shapes of GLAS deviation distributions.

However, these statistics were highly variable among the stations.

This study indicates that accuracy statistics computation is highly dependent to the assumptions made on satellite data and reference data as well. These assumptions can highly affect the computed accuracy statistics and conclusion. Even if the impact of the temporal correlation in GLAS raw data was partially smoothed by the reference data uncertainty, this can be crucial when reasoning with few data, such as a small number of transects, at the station scale.

The accuracy results of this study confirm that satellite ranging LIDAR provides data with a decimeter accuracy to monitor water level in inland water for the Great Lakes or a wide river section.

  Fig. 8 depicts the distribution of the σ G uncertainty parameters for transects as a function of the pulse number composing the transect and the presence or absence of autocorrelation. It appears clearly on this figure that the presence of autocorrelation multiplies the σ G uncertainty by 8 on average. As expected, accounting for autocorrelation drastically changes the σ G estimation, even if the ẐG values differ only slightly (due to the regular time sampling of GLAS measurements). This usual geostatistical result comes from autocorelation that acts as diminishing the number of data.

Fig 1 :

 1 Fig 1: GLAS-ICESat shots forming transects over the Great Lakes and the location of the 27 hydrometric stations.

Fig 3 :

 3 Fig 3: Example of GLAS waveforms on Lake Superior: (a) Unsaturated waveform ; (b) Waveform with low saturation ; (c) and (d) waveforms with high saturation. Ordinate = energy in volts and abscissae = bins number (bins spacing = 1 ns).

Fig 4 :

 4 Fig 4: Example of a GLAS-ICESat transect compared to reference elevation. Abscissae = shot number and ordinate = elevation (ICESat with and without saturation correction and hydrometric station).

Fig 5 :

 5 Fig 5: (a) and (b): temporal series z(t), and (c) and (d): corresponding experimental variograms obtained on the transect at Port Colborne station from 2 February 2005 (a, c) and Bar Point station from 10 May 2007 (b, d). The dashed lines in (c, d) represent the two H0 envelopes, the straight lines in (a, b) represent the reference elevation, and the curve lines in (c, d) are the spherical model function fitted to the experimental variogram (points on c, d). (a, c): absence of autocorrelation. (b, d): significant autocorrelation. (d) Variogram range = 0.32 s.

Fig 7 :

 7 Fig 7: Classification tree explaining the autocorrelation of the transects. A value of 0 on leaves denotes the absence of autocorrelation, and a value of 1 on leaves denotes autocorrelation. The 0/1 value on a leaf means that it mixes correlated and non-correlated transects in equal proportions.

Fig 8 :

 8 Fig 8: Distribution of the standard deviation σG characterizing the uncertainty of estimated water elevations at the transect scale from GLAS data. This distribution is shown as a function of the shot number composing the transect and presence (right) or absence (left) of autocorrelation.

Fig 9 :

 9 Fig 9: GLAS water elevation deviation distribution and upscaling: For transect scale (down), for station scale, for lake scale and all lake scale (up). Not all distributions for lake, station and transect scales are depicted but only examples. The black wide central line on density plots denote the biases values and the thin lines denotes the 2.5 % and 95.5 % quantiles respectively. The two examples of deviation distribution selected at transect scale for each lake are for a correlated transect (wide on the left) and a not-correlated transect (thin on the right).

  scheme, not on the elevation values z. Consequently, given that the temporal sampling scheme of ICESat data is almost regular, except when outliers are excluded from the initial transect data, σ G primarily depends on the transect Author-produced version of the article published in ISPRS Journal ofPhotogrammetry and Remote Sensing, 2011, 66(6), 833-844. The original publication is available at http://www.isprs.org/ DOI : 10.1016/j.isprsjprs.2011.09.002

	total duration.
	When no autocorrelation is observed for a given GLAS transect (H 0 is
	accepted), γ(h) is a pure nugget effect variogram. In this particular case,
	16

  At the Great Lakes 22Author-produced version of the article published in ISPRS Journal ofPhotogrammetry and Remote Sensing, 2011, 66(6), 833-844. The original publication is available at http://www.isprs.org/ DOI : 10.1016/j.isprsjprs.2011.09.002

Table 1 :

 1 Author-produced version of the article published in ISPRS Journal ofPhotogrammetry and Remote Sensing, 2011, 66(6), 833-844. The original publication is available at http://www.isprs.org/ DOI : 10.1016/j.isprsjprs.2011.09.002

	Lake	Id	Station	Number of	Number of	Bias Standard Deviation
				shots used	transects used	(m)	(m)
	Erie	1	Bar Point	975	8	0.008	0.035
	Erie	2	Fermi Power	1634	16	-0.057	0.102
	Erie	3	Marblehead	198	2	-0.012	0.034
	Erie	4	Cleveland	607	8	0.006	0.055
	Erie	5	Fair Port	240	6	0.038	0.054
	Erie	6	Erie	711	10	-0.058	0.112
	Erie	7	Sturgeon	2148	18	-0.051	0.089
	Erie	8	Port Colborne	2034	16	-0.012	0.093
	Erie	9	Port Dover	287	5	-0.044	0.051
	Erie	10	Port Stanley	782	9	0.028	0.061
	Erie	11	Erieau	651	18	-0.038	0.059
	Erie				Mean	-0.021	0.109
	Ontario	12	Port Weller	1220	14	-0.097	0.092
	Ontario	13	Rochester	979	14	-0.031	0.084
	Ontario	14	Oswego	253	5	-0.140	0.034
	Ontario	15	Kingston	221	5	-0.033	0.038
	Ontario	16	Cobourg	1631	17	-0.028	0.074
	Ontario	17	Toronto	743	8	-0.022	0.084
	Ontario				Mean	-0.048	0.087
	Michigan	18	Green bay	213	3	-0.088	0.124
	Michigan	19	Kewaunee	384	4	-0.145	0.056
	Michigan	20	Calumet	571	7	-0.143	0.128
	Michigan	21	Holland	264	4	-0.022	0.084
	Michigan	22	Ludington	425	6	-0.120	0.137
	Michigan				Mean	-0.114	0.129
	Superior	23	Ontanagon	213	15	-0.048	0.093
	Superior	24	Marquette	384	8	-0.091	0.042
	Superior	25	Ross Port	571	12	-0.042	0.062
	Superior	26	Thunder Bay	264	5	-0.004	0.087
	Superior	27 Grand Marais	425	6	0.010	0.033
	Superior				Mean	-0.042	0.079
	ALL LAKES				Mean	-0.046	0.116

Left : Number of GLAS shots and transects used across each hydrometric station of the four Great Lakes; Right: biases and standard deviation of deviations (GLASreference) resulting from comparison between ICESat elevations and reference elevations 31 CemOA : archive ouverte d'Irstea / Cemagref

CemOA : archive ouverte d'Irstea / Cemagref

Acknowledgements

The authors would like to thank the ICESat team and the National Snow and Ice Data Centre for providing ICESat data and expertise. The authors are very grateful to EADS-Astrium and CNES, which supported this study, and especially to Frederic Fabre from EADS Astrium for his useful advices.

The authors also want to thank Carmen Games-Puertas who built several useful tools for ICESat waveforms data extraction during the training period.

Author-produced version of the article published in ISPRS Journal of Photogrammetry and Remote Sensing, 2011, 66(6), 833-844. The original publication is available at http://www.isprs.org/ DOI : 10.1016/j.isprsjprs.2011.09.002

Author-produced version of the article published in ISPRS Journal of Photogrammetry and Remote Sensing, 2011, 66(6), 833-844. The original publication is available at http://www.isprs.org/ DOI : 10.1016/j.isprsjprs.2011.09.002 kriging weights β i in equation 2 are all equal and equation 3 is simplified to the common SDOM formulation. Consequently ẐG equals the arithmetic mean z of the values z i .

GLAS deviations to reference water elevation at transect scale

For each transect j where no linear trend has been observed, a GLAS estimated value ẐG having σ G Gaussian uncertainty was compared to a reference value Z R having σ R Gaussian uncertainty. As GLAS and reference elevations can be considered independent Gaussian random variables, the distribution of the deviations between GLAS and reference water elevations can be assumed also to be of Gaussian distribution N (µ j , σ j ), with µ j = ẐG -Z R and standard deviation

From these Gaussian deviation distributions obtained for a transect, GLAS deviation statistics (accuracy) were computed at station, lake or Great Lakes scales, by simply averaging the cumulative probability functions of the Gaussian laws N (µ j , σ j ) and computing the resulting empirical deviation distribution (not necessarily Gaussian).

From this empirical distribution, the first (bias) and second (standard deviation) moments were estimated to propose GLAS accuracy statistics.

Results and discussion

Autocorrelation between GLAS shot measurments

A significant high frequency autocorrelation is observed for 107 transects, i.e. 45% of the initial 237 transects. No direct explanation is known to the authors for the occurrence of the autocorrelation as shown in Fig. 6, although corrrelated and non-correlated transects look randomly distributed in space (Fig. 6-a), or in time (Fig. 6