Adrien Poupardin

Déborah Idier

Marcello De Michele

Senior Member, IEEE Daniel Raucoules

Water Depth Inversion From a Single SPOT-5 Dataset

Keywords: Bathymetry, cross correlation, SPOT-5, wavelength, wavelets, waves celerity

Knowing bathymetry at intermediate depth, over large areas, and at a reasonable cost is a key issue. Spaceborne remote sensing techniques must play an essential role in retrieving such bathymetry. In this paper, a method is proposed that aims to address this issue without any in situ measurements by exploiting the characteristics of the SPOT-5 satellite dataset. The proposed method is designed to provide bathymetry from two optical SPOT-5 satellite images separated by a time lag DT of 2.04 s. It relies on the estimation of several clouds of wave celerity and wavelength pairs using wavelet and cross-correlation techniques and on the linear wave dispersion relation. This method has been applied to two SPOT-5 images on a test site characterized by complex bathymetry (Saint-Pierre, La Réunion Island). A comparison of the retrieved bathymetry with in situ bathymetric measurements reveals good morphological agreement. The mean relative error is less than 30% in the 3-80-m water depth range. The methodological choices made during method development are discussed based on additional computations, and guidelines for using the proposed method on other images at other sites are provided.

I. INTRODUCTION

T ODAY, knowledge of nearshore bathymetry is essential for multiple applications. For example, bathymetry provides essential data for the study of submarine morphodynamics. These data are vital for planning sustainable coastal development and the conservation of submarine ecosystems. Moreover, they represent a crucial input for nearshore navigation and submarine resources exploration.

The reasons why spaceborne remote sensing techniques must play an essential role in retrieving nearshore bathymetry are threefold. First, spaceborne imagery makes it possible to access remote areas with wide spatial coverage at high spatial resolution. Second, because spaceborne imagery is acquired on a regular basis, a historical data archive is accessible for A. Poupardin was with the French Geological Survey (BRGM), 45000 Orléans, France. He is now with the CEA, DAM, DIF F., 91297 Arpajon, France (e-mail: adrien.poupardin@cea.fr). D. Idier, M. de Michele, and D. Raucoules are with the French Geological Survey (BRGM), 45000 Orléans, France (e-mail: d.idier@brgm.fr; m.demichele@brgm.fr; d.raucoules@brgm.fr).

Color versions of one or more of the figures in this paper are available online at http://ieeexplore.ieee.org.

Digital Object Identifier 10.1109/TGRS.2015.2499379 most sensors, which enables scientists to access information from the past. Third, the cost of the data is relatively affordable compared with airborne or ground missions.

A number of methods have been proven reliable in retrieving nearshore bathymetry from spaceborne remote sensing imagery. However, each of these presents substantial limitations with respect to ground-based (or airborne) methods, either because wave celerity c cannot be directly measured or because ground calibration measurements are needed. Therefore, previous pioneering techniques in this field made use of either a priori information on wave characteristics or ground measurements to create a reliable nearshore bathymetric reconstruction from space.

Among existing techniques, one prominent approach is the use of multispectral satellite imagery, e.g., see [START_REF] Su | Geographically adaptive inversion model for improving bathymetric retrieval from satellite multispectral imagery[END_REF] and [START_REF] Capo | Decadal morphodynamics evolution of a mixed-energy inlet using multispectral SPOT imagery[END_REF]. This method makes it possible to measure bathymetry by assuming the optical depth properties of the water column in shallow waters. This technique has been proven reliable for retrieving shallow bathymetry, i.e., until 20 m approximatively [START_REF] Su | Geographically adaptive inversion model for improving bathymetric retrieval from satellite multispectral imagery[END_REF], when waters are turbidity free. Even for moderately turbid waters, the method may be used, but the bathymetry may be only estimated until 4 m in the case of the work in [START_REF] Capo | Decadal morphodynamics evolution of a mixed-energy inlet using multispectral SPOT imagery[END_REF]. The main drawback of this method is the need to measure a calibration factor on the ground, which limits this technique to accessible areas. Other techniques to estimate bathymetry from optical (e.g., see [START_REF] Leu | Remotely sensing in detecting the water depths and bed load of shallow waters and their changes[END_REF] and [START_REF] Splinter | Bathymetry estimation from singleframe images of nearshore waves[END_REF]) or radar (e.g., see [START_REF] Pleskachevsky | Estimation of underwater topography using satellite high resolution synthetic aperture radar data[END_REF] and [START_REF] Boccia | Tracking of coastal swell fields in SAR images for sea depth retrieval: Application to ALOS L-band data[END_REF]), airborne or spaceborne imagery make use of wave characteristics. Most of these are based on a dispersion relation, which considers linear wave theory and is particularly well adapted to determine bathymetry at intermediate depths (i.e., λ/20 < h < λ/2) as follows:

h = λ/2π • tanh -1 2πc 2 gλ ( 1
)
where h is the water depth, c is the wave celerity, and λ is the wavelength. Here, the dispersion relation is written for the case where the currents are equal to zero. Although λ can be often directly measured from the spaceborne data themselves, these methods need either to make use of an a priori estimate of c (or wave period T ) from a model or to obtain a groundlevel measurement of c (or T ) from buoys, which is a limitation that results in punctual or sparse availability of measurements. Indeed, compared with ground-based or airborne-based methods such as those described in [START_REF] Holman | cBathy: A robust algorithm for estimating nearshore bathymetry[END_REF] and [START_REF] Plant | Ocean wavenumber estimation from wave-resolving time series imagery[END_REF], one key difference is that such airborne and satellite datasets (see [START_REF] Leu | Remotely sensing in detecting the water depths and bed load of shallow waters and their changes[END_REF]- [START_REF] Boccia | Tracking of coastal swell fields in SAR images for sea depth retrieval: Application to ALOS L-band data[END_REF]) are characterized by single images (i.e., no temporal description over a wave period). As a first pragmatic approach, the wave period is often assumed to be uniform. In the case where no wave period data are available, a constant wave period may be used by considering the offshore wave period computed from the deep-water relation between the wave period and the wavelength [START_REF] Dalrymple | Determining depth from remotely sensed images[END_REF]. It is worth noting that the method by Splinter and Holman [START_REF] Splinter | Bathymetry estimation from singleframe images of nearshore waves[END_REF], which is independent of the wave period, provides good bathymetric estimations. Their method exploits spatial changes in wave direction θ to improve the bathymetry estimation and notably in considering the bathymetry gradients without any celerity information. It represents a good improvement to estimate the bathymetry from a single snapshot, but it is not well adapted to more complex configurations. In contrast, the method described in this paper aims at considering the influence of the wave period that may be useful to study a complex sea characterized by several waves superimposing each other. Hereafter, for the sake of clarity, this method will be called the CWB method for correlation, wavelets, and bathymetry. A pioneering method to retrieve bathymetry from satellite stereo image pairs for retrieving both wavelengths and wave celerities has been already described in [START_REF] Abileah | Methods for Mapping Depth and Surface Current[END_REF]. The CWB method differs from it because it determines depth by finding the best fit of the gravity dispersion equation to the measured power spectrum in the same manner as in [START_REF] Piotrowski | Accuracy of bathymetry and current retrievals from airborne optical time-series imaging of shoaling waves[END_REF]. Furthermore, the method proposed in [START_REF] Abileah | Methods for Mapping Depth and Surface Current[END_REF] has been applied to WorldView-2 and IKONOS satellite stereo image pairs characterized by time lags of 10 and 13 s (e.g., see [START_REF] Abileah | Mapping near shore bathymetry using wave kinematics in a time series of WorldView-2 satellite images[END_REF] for an application with WorldView-2), respectively. Such a time lag is quite close to the wave period (which generally ranges from a few seconds to 18 s for very long swells). Depending on the environment, this might lead to ambiguity in wave celerity or frequency detection. Furthermore, in the case of complex seas, it may be difficult to distinguish the different superimposed waves.

The CWB method uses the SPOT-5 High-Resolution Visible sensor that provides panchromatic and multispectral images acquired with a time lag of DT = 2.04 s more adapted to follow the wave in this displacement. The main objective of this paper is to investigate whether by using the peculiar characteristics of the SPOT-5 dataset it is possible to retrieve bathymetry from space with reasonable accuracy and precision without using any in situ measurements or assumptions about the homogeneity of wave characteristics. To tackle the issue of estimating bathymetry using two images only, a method based on cross correlation and wavelet analysis is proposed that exploits the spatial and temporal characteristics of the SPOT-5-like dataset to extract bathymetry. The proposed method combines the spaceborne direct c measurement method presented in [START_REF] De Michele | Direct measurement of ocean waves velocity field from a single SPOT-5 dataset[END_REF] with an original wavelet-based adaptive λ estimate [START_REF] Poupardin | Water depth inversion from satellite dataset[END_REF] to retrieve a spatially dense series (clouds) of (λ, c) couples that are then used to estimate water depth using dispersion relation [START_REF] Su | Geographically adaptive inversion model for improving bathymetric retrieval from satellite multispectral imagery[END_REF].

A number of methods exist to derive wavelength λ from both optical and radar spaceborne imagery. To measure λ, a dataset is needed with fine enough spatial resolution to isolate the shortest wavelengths. One of the first developments in measuring λ was proposed by Populus et al. [START_REF] Populus | The use of SPOT data for wave analysis[END_REF], who used a two-dimensional discrete Fourier transform applied to SPOT data at a spatial resolution of 10 m. They could distinguish between swells and wind-generated waves, which are locally generated and therefore have a shorter wavelength than swells. However, this method fails to correctly determine wavelengths and directions near the coast because waves evolve rapidly.

To solve this problem, several authors have developed twodimensional continuous wavelet analyses of the sea state [START_REF] Feng | Application of two-dimensional wavelet transform in near-shore X-band radar images[END_REF]- [START_REF] Wu | Ocean remotely sensed image analysis using two-dimensional continuous wavelet transforms[END_REF]. These methods, which are based on ground photos, show good accuracy in characterizing waves, even in coastal areas. Chuang et al. [START_REF] Chuang | Twodimensional continuous wavelet transform of simulated spatial images of waves on a slowly varying topography[END_REF] tested a method on a synthetic wave-field dataset under various wave conditions and seabed slopes. They showed that two-dimensional continuous wavelets are well adapted to characterizing nonstationary and inhomogeneous wave fields. In particular, the method in [START_REF] Chuang | Twodimensional continuous wavelet transform of simulated spatial images of waves on a slowly varying topography[END_REF] shows good precision in detecting wave refraction and reflection, taking shoaling effects into consideration, and in determining wave spectra at any location in the coastal area. Wavelet function parameters were also discussed by Wu et al. [START_REF] Wu | Ocean remotely sensed image analysis using two-dimensional continuous wavelet transforms[END_REF], who used a set of X-band synthetic aperture radar data and panchromatic QuickBird data to extract wave characteristics. Their method showed good results when retrieving ocean wave inhomogeneity and nearshore refraction and for describing local wave information at several locations on the scene, particularly near the shore. In [START_REF] Gorman | Directional wavelet analysis of inhomogeneity in the surface wave field from aerial laser scanning data[END_REF], aerial laser scanning measurements were used along with a continuous wavelet method to analyze the sea state in the spatial domain. Gorman and Hicks [START_REF] Gorman | Directional wavelet analysis of inhomogeneity in the surface wave field from aerial laser scanning data[END_REF] considered wave fields as greatly inhomogeneous and found that the Fourier transform method was not sufficient to retrieve wave characteristics at precise locations.

Variables c or f must be determined to estimate water depth h from dispersion relation [START_REF] Su | Geographically adaptive inversion model for improving bathymetric retrieval from satellite multispectral imagery[END_REF], with f related to c by f = 2πc/λ. Literature contains many methods to measure c or wave frequency f based on optical images. Apart for the study by Abileah [START_REF] Abileah | Methods for Mapping Depth and Surface Current[END_REF], [START_REF] Abileah | Mapping near shore bathymetry using wave kinematics in a time series of WorldView-2 satellite images[END_REF] and de Michele et al. [START_REF] De Michele | Direct measurement of ocean waves velocity field from a single SPOT-5 dataset[END_REF], none of them makes use of spaceborne imagery. As an example, the coastal observation techniques described in [START_REF] Holman | cBathy: A robust algorithm for estimating nearshore bathymetry[END_REF], [START_REF] Plant | Ocean wavenumber estimation from wave-resolving time series imagery[END_REF], and [START_REF] Hampson | Video-Based Nearshore Depth Inversion Using WDM Method[END_REF] used wellresolved airborne or ground-based multitemporal imagery to obtain wave frequency f directly. In [START_REF] Holman | cBathy: A robust algorithm for estimating nearshore bathymetry[END_REF], images were collected at an acquisition frequency of 2 Hz during 1024 s, which is not yet possible from a spaceborne platform. Even with these ground-based experimental configurations, wavelength λ must be deduced using dedicated techniques. For instance, the works in [START_REF] Hampson | Video-Based Nearshore Depth Inversion Using WDM Method[END_REF] and [START_REF] Donelan | Non-stationary analysis of the directional properties of propagating waves[END_REF] developed the wavelet directional method to determine wavelength λ from the phase differences among spectral wavelet transforms of several couples of points in the spatial domain. In this method, the wavelets were applied to a water-surface elevation time series (estimated using ground equipment). In a similar manner, Holman et al. [START_REF] Holman | cBathy: A robust algorithm for estimating nearshore bathymetry[END_REF] performed temporal Fourier transforms on several points in the spatial domain before computing the cross-spectral matrix, providing estimates of the dominant eigenvectors and their associated wavelengths λ. Note that, in the studies described in [START_REF] Holman | cBathy: A robust algorithm for estimating nearshore bathymetry[END_REF] and [START_REF] Plant | Ocean wavenumber estimation from wave-resolving time series imagery[END_REF], a cross correlation is realized to find the time delay DT that maximizes the correlation for a given frequency f , whereas in the CWB method, a cross correlation is realized to find the wave displacement that maximizes the correlation for a given wavelength λ. In this last case, using two images separated by a small temporal lag DT, it is preferable to work with (c, λ) couples rather than (f, λ) couples that are more adapted to time series. Some methods to resolve f and λ are particularly suited for airborne datasets [START_REF] Piotrowski | Accuracy of bathymetry and current retrievals from airborne optical time-series imaging of shoaling waves[END_REF]. In the experiment described in [START_REF] Piotrowski | Accuracy of bathymetry and current retrievals from airborne optical time-series imaging of shoaling waves[END_REF], frames were collected at a frequency of 2 Hz during 480 s and a spatial resolution of 2 m for a spatial domain of 2 km × 2 km. With this optimal configuration, high-spatial-resolution and high-temporal-resolution f /λ spectra can be computed, yielding h with high precision. In their method, Piotrowski and Dugan [START_REF] Piotrowski | Accuracy of bathymetry and current retrievals from airborne optical time-series imaging of shoaling waves[END_REF] used the fact that the largest waves are more strongly impacted by the sea bottom. Then, they considered a wavelength range in the dominant wave direction (i.e., for the long waves impacted by the bottom) and the corresponding frequency range that maximized the spectral energy. Then, they chose the h value that minimized the sum of squares of the distances between the point cloud found determined by the f /λ spectrum and the characteristic dispersion curve given by [START_REF] Su | Geographically adaptive inversion model for improving bathymetric retrieval from satellite multispectral imagery[END_REF].

Ground-based and airborne observation techniques are well adapted to extract wave characteristics and bathymetry at high spatial resolution on well-defined local areas only. Today, because of cost or because of the inaccessibility of the study area, they cannot be routinely used to cover very large surfaces. Moreover, this approach cannot be used to retrieve past bathymetry. Furthermore, ground-based image analysis methods cannot cover a large spatial extent. These limitations were the motivation for developing the spaceborne method proposed here.

This paper proposes a methodology based on publicly available (commercial or not) spaceborne optical imagery to determine h, which only requires one pair of images characterized by a short enough time lag.

This paper is structured as follows. Section II describes the proposed method. Section III presents the application of the method to the selected SPOT-5 dataset to provide bathymetric estimates, which are compared with in situ measurements. In Section IV, the results and the method are discussed, as well as the applicability of the CWB method to other satellite datasets. Finally, conclusions are drawn in Section V.

II. CWB METHOD

In the procedure for estimating bathymetry from a single SPOT-5 dataset, the two main original ideas of the proposed method are: 1) to use wavelet analysis and cross-correlation techniques to estimate wavelengths λ and celerities c; and 2) to consider several clouds of (λ, c) couples to estimate bathymetry. Some of the underlying ideas have been already developed and used within ground-based or airborne approaches. However, the characteristics of spaceborne images (only two images for the temporal dimension, at plurimetric spatial resolution) imply specific issues compared with ground-based or airborne methods so that the present method is original when applied to satellite-based bathymetric estimation. The methodology is illustrated in Fig. 1, and details are given in Section II-A-C. In the general case, the input data consist of two images, possibly with different spatial resolutions (this is the case, for instance, for SPOT-5 panchromatic and multispectral images). In the following paragraphs, the term high-resolution scene (SH) is used to refer to the SPOT-5 image with higher spatial resolution (i.e., the panchromatic image), whereas the term low-resolution scene (SL) is used to refer to the lower resolution SPOT-5 image (i.e., the multispectral image). The corresponding spatial resolutions are designated by DX SH and DX SL , respectively. 3(b)]. The selected ĥ is indicated with a thicker line (here, it is h 1 ). The theoretical curve obtained using the in situ bathymetry is shown by a dotted line.

The principle of the CWB method is, for every subscene of the image where (λ, c) are estimated, to identify several waves (called dominant waves), to estimate several (λ, c) pairs (called a cloud) for each dominant wave, to use the obtained cloud to estimate one water depth per dominant wave, and then to select the "best" ĥ water depth among the various estimated IEEE TRANSACTIONS ON GEOSCIENCE AND REMOTE SENSING depths. This principle is illustrated in Fig. 2 for one point P1 in the spatial domain [indicated by a white cross in Fig. 3(b)].

In this example, four dominant waves were detected, providing four (λ, c) cloud pairs and, therefore, four estimates of the water depth. The thick line indicates the selected ĥ value. Table I presents the parameters of the method, and its application is described in the following.

A. Wavelet Analysis of Higher Resolution Scene SH

A wavelet analysis of higher resolution scene SH was carried out to obtain wavelengths with higher precision than those from a classical Fourier transform method. Indeed, the wavelengths rapidly vary on distances that are close to the wavelengths themselves. The wavelet analysis is well adapted in this case because it may be computed at each pixel, providing a precise value of the wavelength on these pixels. Furthermore, the Fourier analysis considers all the waves and noise included in the entire work subscene. In the case in [START_REF] Boccia | Tracking of coastal swell fields in SAR images for sea depth retrieval: Application to ALOS L-band data[END_REF], smoothing of the wave spectrum was performed to only keep the swell. On the contrary, the wavelet analysis eliminates a part of the noise and of the waves that are not centered on the work subscene, leading to a higher spectrum quality without using any filters, as shown in [START_REF] Wu | Ocean remotely sensed image analysis using two-dimensional continuous wavelet transforms[END_REF], where an example of the comparison between the two analyses is given.

The wavelets used to characterize the image may be scaled, shifted, and rotated, which makes it possible to detect all the waves at a particular location even if they are superimposed on each other. In the study of sea-surface waves, the most common and useful wavelets used in literature are the Morlet wavelets (see [START_REF] Feng | Application of two-dimensional wavelet transform in near-shore X-band radar images[END_REF]- [START_REF] Gorman | Directional wavelet analysis of inhomogeneity in the surface wave field from aerial laser scanning data[END_REF]). Therefore, the Morlet wavelet was selected for use in the proposed method. The associated wavelet functions are defined by

Ψ a, b,θ ( x) = 1 a exp -x -b 2 2a 2 exp 5i( x -b) • r a ( 2
)
where r = (cos θ, sin θ), and x is a two-dimensional spatial vector.

The SH scene (i.e., the scene having higher spatial resolution) is convoluted by the wavelet functions given spectral matrix W as follows:

W a, b,θ ( x) = R 2 I( x)Ψ a, b,θ ( x)dx 2 (3)
where I( x) is the subscene function.

Once wavelet transform matrix W has been computed, it can be expressed in the form W x (a, θ) = W x (λ, θ), with λ related to a by λ = 2π a L S /5 and with L S × L S being the dimension of the subscene function under consideration, i.e., I( x). This means that, for a particular location x, matrix W only depends on the wavelength and the wave direction. Morlet wavelets were applied to subscenes of variable dimension (L S × L S ) so that each subscene contains at least two wave crests and two wave troughs. L S was estimated based on preliminary wavelength estimates (fast Fourier transform or visual inspection). To avoid boundary effects, only the wave spectra in a twice smaller area (dimension L S /2 × L S /2) centered on each subscene are used in the following steps. The wavelet spectra are computed using input subscenes with resolution DX WA . Theoretically, wavelet analysis makes it possible to compute wavelength spectra at the same resolution (DX SH ) as the pixels of the SH scene (e.g., for panchromatic SPOT-5 images). If the computation time must be decreased, DX WA can be reduced, but it still should be small enough to characterize the wavelengths properly. Then, by spatial averaging, mean wavelength spectra are computed at a resolution of DX WO . Then, the subscenes are recovering each other to cover all the scenes with regularly spaced wavelength spectra W x (λ, θ), at resolution DX WO .

Once wavelength spectra W x (λ, θ) have been computed, N dominant wavelengths are numerically extracted to cover a wide range of waves, e.g., several superimposed waves. To achieve this, local energy maxima |W x (λ n , θ n )| are assumed to identify the (λ n , θ n ) polar coordinates (n varying from 1 to N ) associated with each main local wave climate. In this method, N will be always less than or equal to the fixed maximum number of local maxima, which is called N max , to avoid excessive and unnecessary computation time (this point will be further discussed in Section IV). A representation of a wavelength spectrum obtained at point P1 [located in Fig. 3(b)] is given in Fig. 4. In this case, four significant dominant waves were detected. If more than N significant waves can be theoretically produced, the spectral mode search can be limited to the more relevant waves in terms of energy. At this step, one may note that, most of the time, the wind wave is rejected because it does not have sufficient energy. On one hand, it is necessary to robustly localize the local maxima with respect to their neighborhood and to spectral noise. On the other hand, it is clearly desirable to identify wavelength domains (large enough to take advantage of the dispersion phenomenon) around each maxima. It is therefore necessary to focus on spectral peaks with high energy compared with other wavelengths and peaks spread over significant wavelength ranges. Indeed, the dominant waves, i.e., wind waves or swells, are associated to bumps distributed around a local maximum and not just to a peak value of the Dirac type. These bumps are well described by the clouds (λ, c). Some studies had been led to distinguish the wind wave from the swell, e.g., as in [START_REF] Aarnes | Partitioning sequences for the dissection of directional ocean wave spectra: A review[END_REF]. Other methods also use the detection of local maxima to detect the dominant waves in the wave spectrum, as in [START_REF] Boccia | Tracking of coastal swell fields in SAR images for sea depth retrieval: Application to ALOS L-band data[END_REF]. This choice will be further discussed and its relevance illustrated in the application in Section III.

B. Determination of (λ n,m , c n,m ) Couples and Water Depth h n Associated With Each Dominant Wave

Before computing the wave celerity, in the case of images having different spatial resolutions (this is the case for panchromatic and multispectral SPOT-5 images), the lowest resolution scene (SL) should be oversampled to match the higher resolution scene (SH) to correlate images of the same step size DX SH . Oversampling of the SL scene is preferable to downsizing of the SH scene to improve the precision on the wave displacement. This operation is carried out by applying a bicubic interpolation to the SL scene.

In this paper, Fourier-filtered subscenes are not used to compute celerity, as in [START_REF] De Michele | Direct measurement of ocean waves velocity field from a single SPOT-5 dataset[END_REF] or [START_REF] Poupardin | Water depth inversion from satellite dataset[END_REF]. Indeed, the spatial resolution given by the Fourier transform subscenes is too coarse to isolate dominant wavelength λ n as determined by the wavelet analysis, particularly for the largest wavelengths. Once the wavelet spectrum has been obtained and the energy maxima and their corresponding wavelength (λ n ) are identified, the method consists of selecting several wavelength values λ n,m around dominant wavelength λ n , keeping θ = θ n and m varying in the [1, M ] interval but corresponding to sufficient spectral energy to distinguish them from background noise. The next step is to compute the celerities associated with each λ n,m . This approach is inspired by the method proposed in [START_REF] Piotrowski | Accuracy of bathymetry and current retrievals from airborne optical time-series imaging of shoaling waves[END_REF] that considered several (f, λ) couples in the dominant wave direction but adapted to the lack of access to frequencies on the basis of a high-temporal-resolution camera. Only one pair of two images separated by a phase lag (DT) of a few seconds can be accessed. Considering spectrum W x at a point on the grid, M couples (λ n,m , θ n ) with spectral energy values larger than a threshold related to the energy maximum are defined as follows: TC × |W x (λ n , θ n )|. This approach aims to only take into account (λ n,m , θ n ) couples whose energy is still relatively close to one of the n dominant modes (λ n , θ n ).

Next, the celerity can be estimated for each wavelength. At each point on the DX WO grid, for each scene, subscene wavelet transforms were recomputed only considering direction θ n and wavelength λ n,m . Contrary to Section II-A, where resolution DX WA was used to reduce the computation time, to compute the wave displacements, the subscenes are considered on their full resolution DX SH (i.e.,. at the pixel resolution). It does not represent a loss in time because the computation is realized for only one set of wavelet parameters (direction θ n and wavelength λ n,m ). Relations (2) and (3) were used with a = 5λ n,m /2πL and θ = θ n to compute the two spectra, which were expressed in the form W n,m ( x) = W λ n,m , b,θ n ( x), i.e., one spectrum for each subscene (SH, SL). Then, for each subscene, the real part of the spectrum was used to identify the wave shape precisely, from the crest to the trough. A crest is located when the wavelet and the subscenes are in phase, and a trough is located when the wavelet and the subscenes are out of phase. The cross correlation was applied to a spectral window of dimension L S /2 × L S /2, giving cross-correlation matrix CC n,m (i, j), i.e.,

CC n,m (i, j)= S-1 k=0 S-1 l=0 Re W SH n,m (k, l) ×Re W SL n,m (k-i, l-j) ( 4
)
where S is the window size (i.e., the number of pixels), and (i, (λ n,m ,c n,m ) is obtained, which is represented in the (λ, c) plane and completes the characterization of the nth dominant wave regime (λ n , c n ). Then, using a Levenberg-Marquardt (LM) algorithm, the water depth h n that minimizes the sum of squares of the distance between the point cloud and the dispersion curve given by (1) can be estimated. To perform this adjustment, a weight is associated with each point (λ n,m , c n,m ) depending on the quality of the correlation Q n,m . Fig. 2 provides an illustration for a case with four dominant waves.

C. Choice of Water Depth ĥ Among the N Bathymetric Values h n

On each point of the (DX WO × DX WO ) grid, the proposed method is applied to the N dominant waves and therefore provides N depth values h n , with n varying from 1 to N (for an illustration, see Fig. 2). A first natural approach to select the "best" water depth ĥ among the N bathymetric values h n would be to use the error on h n estimated during the aforementioned LM inversion. However, due to the nonlinearity of the dispersion relation, the reliability of the obtained h n is not really characterized by its deviation on h n but rather by the quality of its associated (λ n,m , c n,m ) cloud and how it follows the h-dispersion curves. For instance, for a single metric error on h n , depending on the depth, a cloud with characteristics that are not the best aligned with the h-dispersion curve could be chosen, depending on the h n candidate values. Therefore, an alternative approach is proposed here, which focuses on proper wave detection in the SH and SL subscenes so that the nth dominant wave that maximizes product PW

= |W SH x (λ n , θ n )| × |W SL
x (λ n , θ n )| can be identified. Noteworthy is the fact that the dominant modes must be found on both panchromatic and multispectral scenes. Otherwise, the dominant waves will be rejected by our criterion PW. In other words, a mode is determined on the panchromatic image for having the best possible precision, but it is only used for h-derivation if the mode is also detected using the multispectral image. That points out the importance of the PW criterion for rejecting a posteriori irrelevant modes (i.e., modes that are observed in SH but not in SL) during the process.

Within the method, the dispersion of the (λ n,m , c n,m ) cloud around the fitted dispersion curve is largely due to the lack of precision in celerity estimate c n,m . To avoid using clouds characterized by too large celerity errors, the two following rejection conditions of estimate h n associated with the nth cloud were used: 1) if the standard deviation of the celerity within the nth cloud is larger than a given threshold (the parameter called C c , see Table I), then h n is rejected; and 2) an additional threshold is used so that all estimated h n outside the given range of bathymetric values are rejected (a parameter called DZ, see Table I). The possible causes of threshold violation can be various kinds. First, the displacement estimate can be of poor quality for low-spectral-energy modes, resulting in a high standard deviation on the celerity. It is worth noting that the selection and sorting of the first five modes by decreasing the PW value does not assess the global quality of their estimates but rather their relative importance in the image; therefore, the possibility that the first modes result in poor c n,m estimation cannot be rejected. Second, for a given depth and the same standard deviation on c n,m , the process can result in aberrant h n values for the shortest wavelengths (where the h-curves corresponding to strongly different h values are too close to each other with respect to the amount of precision in c n,m ). However, the main factor used in selecting h n is based on maximizing product PW, with the two thresholds C c and DZ being only secondary factors. The influence of these thresholds on the final estimated bathymetry, as well as the advantages of using such thresholds, will be discussed in Section IV.

III. APPLICATION AND RESULTS

A. Site and Data Description

To investigate the potential of the developed method, it was applied to a SPOT-5 image pair acquired on February 6, 2010 at 10:30 A.M. in the area of Saint-Pierre, which is located on the southwest of La Réunion Island (Fig. 3 shows the panchromatic scene). It occurs at the falling tide, i.e., 1 h and 30 min before the low tide. This experimental area is characterized by a short continental shelf with a width of a few kilometers. On this site and this date, there is no current measurement or hydrodynamic model results of high spatial resolution. Thus, we used the results of the global model Navy Coastal Ocean Model of the National Oceanic and Atmospheric Administration. At the global scale, this model provides local surface currents of about 0.1 and 0.14 m/s on February 6, 2010 at 6 and 12 A.M., respectively. These currents remain weak enough to not be considered in the dispersion relation. It is the same site as that used in [START_REF] De Michele | Direct measurement of ocean waves velocity field from a single SPOT-5 dataset[END_REF]. Furthermore, this site is covered by in situ measurements realized by LIDAR (from Litto3D products, see www.shom.fr for more details) for depths less than 30 m and by a singlebeam echo sounder (SHOM (French Navy Hydrographic and Oceanographic Department) data, see data.shom.fr) for deeper water. These in situ measurements are used as a reference to validate the proposed method. A representation of this measured bathymetry is shown in Fig. 6(a).

SPOT-5 is a satellite equipped with a High-Resolution Geometry (HRG1-2) sensor that acquires data in multispectral mode (HS) at a spatial resolution of DX SL = 10 m and in panchromatic mode (HMA) at a spatial resolution of DX SH = 2.5 m. The temporal shift between the HMA and HS acquisitions is DT = 2.04 s. As shown in [START_REF] De Michele | Direct measurement of ocean waves velocity field from a single SPOT-5 dataset[END_REF], this configuration makes it possible to compute wave celerities in correspondence with the displacement of a cluster of pixels over time.

In conjunction with the HMA band, three multispectral bands (XS1-3) can be used to estimate wave celerity. The XS1 multispectral band (green) with a spectral bandwidth centered on 0.55 μm has been chosen because it is the band whose spectral bandwidth best matches that of the HMA panchromatic band (for the ocean surface). This consideration is important to obtain good image cross-correlation scores [START_REF] Necsoiu | Monitoring migration rates of actives subarctic dune field using optical imagery[END_REF]. A spatial shift was found between the XS1 multispectral band and the HMA panchromatic band. Because it can be assumed that the land-based parts of both bands do not move in the temporal lag DT of 2.04 s, this spatial shift can be measured on land and then compensated for by applying a constant offset to the XS1 multispectral band.

B. Application of the Method

Within the CWB method, several types of parameters must be set, i.e., some related to satellite image characteristics, others from the site and from observed waves, and others using the conventional method. Table I contains all the parameter values associated with the application of the method to the study site and the use of SPOT-5 images.

The subscene sizes (L S ranging from 640 m offshore to 320 m nearshore) were estimated on the basis that the more energetic wavelengths typically range from approximatively 250 m offshore to approximatively 100 m near the coast where the bathymetry is much shallower, enabling the consideration of at least two wave crests and two wave troughs for each subscene. Subscene dimensions and resolutions are shown in Fig. 7. In this application, the two images have different spatial resolutions (DX SH and DX SL ), which means that the computation of the dominant waves was done on the panchromatic image (resolution DX SH of 2.5 m), and the multispectral image of resolution DX SL (i.e., 10 m) was oversampled at the resolution DX SH of the panchromatic images.

Theoretically, wavelet analysis makes it possible to compute wavelength spectra at the same resolution as the panchromatic pixels (i.e., 2.5 m), but to save time, this resolution was reduced to DX WA = 10 m offshore and to DX WA = 5 m near the shore. This resolution allows detecting modes with wavelengths superior to 20 m offshore and to 10 m near the coast. Then, mean wavelength spectra were computed at a resolution DX WO of 20 m. This means that the average took into account 2 × 2 wavelength spectra offshore and 4 × 4 wavelength spectra near the coast. This choice was made as a compromise between the computation time and the physical description of wave characteristics. Indeed, a preliminary analysis (not shown here) based on wavelength variograms showed that the 20-m value was always less than the correlation of the detected wavelength (from 100 to 200 m depending on the spatial areas considered For the maximum number of modes (i.e., clouds) for which bathymetric values were estimated, an N max value of 5 was selected. For the maximum number of (λ, c) pairs per cloud, a value of M = 16 was chosen. For threshold T C related to the energy maxima of the (λ n , c n ) couples, based on some preliminary tests, a value of 0.65 was selected. Fig. 2 provides an example of the clouds obtained for point P1 [location indicated in Fig. 3(b)]. In this case, four dominant modes, and therefore four clouds, were identified. This provides a reasonable number of points to estimate bathymetry h n (see the illustration in Fig. 2). For the last step of water depth selection among the N values, C C was set to 1 m/s, whereas depth range DZ was set to [1-100 m]. The thick line in Fig. 2 indicates the selected ĥ value for location P1, which was used for comparison with the dispersion curve obtained using the in situ water depth (obtained from LIDAR and multibeam measurements).

As a comment on the relevant N max value, analysis of the selected water depth ĥ shows that, most of the time, the first or second mode was selected for the final bathymetry estimation (see Fig. 8), which indicates that there is no need to analyze more than N max = 5 dominant waves. From the CWB selection process (one criterion and two thresholds), long waves were most often selected. Note that the longest waves are most impacted by the bottom [START_REF] Holman | cBathy: A robust algorithm for estimating nearshore bathymetry[END_REF].

C. Results

Bathymetry ĥ, which was estimated from the SPOT-5 dataset, is shown in Fig. 6(b) on an 80 m × 80 m grid, which was obtained by spatially averaging the 20-m-resolution bathymetric grid. It has been also corrected by taking tides into account so that the vertical reference of the obtained bathymetry is the hydrographic zero (SHOM convention, corresponding more or less to the lowest tidal level). Indeed, the sea level was 0.51 m above the reference level (hydrographic zero) when the SPOT-5 satellite took the scenes on February 6, 2010 at 10.30 A.M. In Fig. 6(b), the isobaths correspond to in situ measurements h in situ . It should be remembered that the area after wave breaking is not covered by the CWB method [in black in Fig. 6(b)]. Qualitative good agreement with observations can be seen: The proposed method properly reproduces the main morphologic features. For instance, the short continental shelf is described, with bathymetry reaching a 50-m depth at 1 km from the coast. The estimated bathymetry ĥ decreases faster near the harbor (right part of the scene) than where the coast exhibits a concave curve (left part of the scene). Furthermore, the method succeeded in identifying the submerged headland at the (3.34 × 10 5 m, 76.36 × 10 5 m) Universal Transverse Mercator location [see Fig. 6(b)].

To quantify the quality of the results, maps of the absolute error (AE = ĥh in situ ) and the relative error (RE = ( ĥh in situ )/h in situ ) are shown in Fig. 9(a) and (b), respectively. The relative error is strongly nonuniform, ranging overall from -100% to +100%, with most values ranging from -50 to +50%. Errors very close to 0% can be also observed, and it corresponds to the location where the CWB method works better because the dominant waves may be clearly identified on the scene.

To analyze these errors versus water depth, for each 80 m × 80 m bathymetric cell, ĥ is plotted versus h in situ [see Fig. 10(a)]. First, the scatter decreases with the water depth until approximatively 3 m and then increases with the water depth. Second, with regard to quantile, we have the following: 1) for water depth less than 3 m, the estimated bathymetry overestimates the real bathymetry; 2) for water depth in the range [3-20 m], the estimated bathymetry exhibits fairly good agreement with in situ observations; and 3) for water depth greater than 20 m, the estimated bathymetry underestimates the real bathymetry. These discrepancies may be attributed to a nonoptimal selection of depth h n among the N dominant waves at a given location or to an inaccurate estimate of the celerities. Furthermore, some dominant waves may not have been detected at several locations of the basic panchromatic scene [see Fig. 3(b)] because of low wave amplitudes.

To complete the analysis of errors versus water depth, the average absolute error (ER 1 = 1/K × Σ [h1,h2] |AE|) and the average relative error (ER 2 = 1/K × Σ [h1,h2] |RE|) were computed for several h in situ corresponding to the h in situ quantiles. Such error indicators were used to avoid negative and positive errors that counterbalance each other. Fig. 10(b) shows that the absolute error first decreased with the water depth for the first two quantiles and then increased with the water depth from 1.5 to 22 m, whereas the relative error remained almost constant, ranging from 20% to 30%, except for the two smallest water depth quantiles.

IV. DISCUSSION

To complete the analysis of errors versus water depth, the average absolute error (ER 1 = 1/K × Σ [h1,h2] |AE|) and the average relative error (ER 2 = 1/K × Σ [h1,h2] |RE|) were computed for several h in situ corresponding to the h in situ quantiles. Such error indicators were used to avoid negative and positive errors that counterbalance each other. Fig. 10(b) shows that absolute error.

A. Error Comparison With [3]

Most bathymetric methods using optical satellite images are based on wavelength estimation and assume a uniform wave period. Here, a comparison of the error magnitudes obtained with the CWB method and with the most widely used method will be presented. Differences in images, sites, investigated water depth range, and wave characteristics make impossible any rigorous comparison of the quality of the results with the results obtained by other authors using other techniques. However, to make the various configurations as comparable as possible, the problem can be made partially dimensionless based on wavelength characteristics. This leads to the following dimensionless numbers: R1 = λ/h 1 , R2 = λ/h 2 , and R3 = λ/dx, where λ is a characteristic offshore wavelength, h 1 and h 2 are the lower and upper bounds of the depth range for which the errors are quantified, respectively, and dx is the bathymetric grid resolution. This approach is used to compare the results presented here with those in [START_REF] Leu | Remotely sensing in detecting the water depths and bed load of shallow waters and their changes[END_REF] obtained assuming a uniform wave period. In that study, they obtained a mean relative error ER range. This error is of the same order as that reported in [START_REF] Leu | Remotely sensing in detecting the water depths and bed load of shallow waters and their changes[END_REF]. IEEE TRANSACTIONS ON GEOSCIENCE AND REMOTE SENSING However, it should be noted that the bathymetry in the Saint-Pierre area is particularly difficult to estimate because the bed slope is two to three times steeper than in [START_REF] Leu | Remotely sensing in detecting the water depths and bed load of shallow waters and their changes[END_REF], the bathymetric lines are not parallel, and the seabed is of coral origin and is therefore characterized by small-scale bathymetric heterogeneity.

B. Methodological Choices

The CWB method relies on several methodological choices made during its development. To analyze the outcomes of the proposed method more thoroughly, the added value of the following methodological choices will be illustrated: (C1) using (λ, c) couples rather than λ and assuming a uniform period to estimate bathymetry; (C2) using the real part of wavelet analysis rather than Fourier-filtered scenes to compute c; (C3) considering clouds rather than a single pair (λ, c) to estimate bathymetric values h n ; (C4) choosing the "best" water depth ĥ among the N bathymetric values h n rather than using all couples (λ, c) to estimate a single bathymetric value; and (C5) using c-STD (the celerity standard deviation) and bathymetric range thresholds to reject aberrant values rather than only using energy criteria. To investigate these choices, additional tests were performed on the study site using the same image dataset. Indeed, this test case is a very challenging case: It is characterized by nonlinear bathymetric lines, strong bathymetric heterogeneity, and the overlapping of swells and wind waves.

To assess the added value of choice C1 more thoroughly, the bathymetry was estimated only using wavelength values and assuming a uniform wave period. The wave period is usually obtained from offshore buoy measurements. On the study site, there is no offshore buoy. Therefore, the offshore wave period was estimated based on offshore wavelengths and celerity values (72.4 m and 16.1 m/s, respectively). These values provided a period value of 10.7 s. This value is consistent with wave modeling results obtained in this area for the time of acquisition of the SPOT-5 images. Then, the bathymetry was computed using the dominant mode wavelength. Fig. 11(a) shows that the relative errors of this method are generally larger than those associated with the CWB method [see Fig. 8(b)]. Focusing on specific areas, it is apparent that, for instance, in the northwest, in the 20-25-m depth area, the simplified method overestimated the bathymetry by approximatively 60%, whereas the CWB method provided an error close to zero. In the center of the image, there is an area with an approximate depth of 15-20 m where the CWB method also provided an error close to zero, whereas the simplified method underestimated the bathymetry by approximatively 60%. This comparison illustrates the potential improvement provided by using (λ, c) couples over the most widely used method, which is based on the assumption of a uniform period.

As for the improvement obtained using choice C2, Fig. 11(b) and (c) shows the bathymetric errors for the two following cases: only considering dominant waves and using Fourierfiltered scenes to compute the celerity [see Fig. 11(b)] and only considering dominant waves but using the real part of wavelet analysis to compute the celerity [see Fig. 11(c)]. This comparison shows that, in the first case, water depths are strongly underestimated, particularly offshore, whereas in the second case, the errors are smaller (globally reduced by a factor of two) than when only considering dominant waves and using Fourier-filtered scenes to compute the celerity. This illustrates the added value of using the real part of wavelet analysis for the subscene correlation and, thus, for the celerity computations.

As for the improvements obtained using choice C3, Fig. 11(c As for choice C4, Fig. 11(d) shows the bathymetric errors when considering all clouds together to compute a single water depth. The relative errors are larger than those obtained using the CWB method [see Fig. 9(b)], particularly in water depths less than 40 m. Indeed, at some locations (not shown here), the clouds can be very different. For instance, within a cloud, the celerity may appear to be severely underestimated compared with other clouds. In that case, the cloud method makes it possible to avoid mixing well-resolved and poorly resolved clouds.

Finally, the influence of using thresholds on the c-STD and bathymetric range DZ was investigated (see Table I). Unlike previous comparisons that were made on the 80 m × 80 m grid, this analysis was directly done on the 20 m × 20 m grid. First, using the c-STD threshold implies a bathymetric change (i.e., a change in the selection of the "best" water depth) over 17% of the bathymetric grid, mainly for water depths greater than approximatively 40 m or where the bathymetric slope is steepest. Using the bathymetric range implies changes in ĥ over 19% of the grid, mainly offshore (water depth greater than 40 m) and near the coast (water depth less than 1 m). Using both leads to changes over approximatively 30% of the grid, leading to a significant reduction in relative bathymetric errors, mainly offshore (water depth greater than 40 m) and near the coast (water depth less than 1 m). Focusing on the [1-40 m] bathymetric area, using both thresholds leads to changes over approximatively 5% of the grid, mainly in the steepest areas. Therefore, these thresholds mainly contribute to improving the bathymetry near the coast and offshore, whereas they do not significantly modify it in the intermediate bathymetric range [1-4 m]. For water deeper than 40 m, the waves that are no longer impacted by the bottom must be rejected to retrieve the bathymetry. The threshold on the c-STD is used to avoid using the shortest irrelevant waves. The second threshold based on DZ acts similarly. As the dispersion relation is no longer valid for water depth in the range λ/20 > h > λ/2, the corresponding waves shall not be used to retrieve the bathymetry. In the presence of several superimposed waves, the best option is to consider one (or several) wave(s) that corresponds better to the relation dispersion domain of validity rather than using all the detected waves.

This analysis illustrates the improvements brought about by using (λ, c) pairs to estimate bathymetry (C1), using the real part of wavelet analysis to compute c (C2), considering clouds and estimating one bathymetric value per cloud (C3), choosing the "best" bathymetric value (C4), and using c-STD and bathymetric range thresholds to complement the "best" bathymetry selection obtained only considering spectral energy (C5).

C. Dispersion Relation: Limits and Uncertainty Propagation

One source of uncertainty in the CWB method is related to the linear wave dispersion relation itself. This equation is based on the linear wave theory, which provides a linearized description of the gravity-wave propagation on the surface of a homogeneous fluid layer. One of the main assumptions is that the fluid layer has a uniform mean depth. In the present problem, the water depth is not uniform, meaning that in bathymetric estimation, this equation is being used beyond its theoretical application. However, most of these studies (including this paper) still provide reasonable estimates of water depths.

For further discussion on the limitations of the linear dispersion equation, see [START_REF] Holland | Application of the linear dispersion relation with respect to depth inversion and remotely sensed imagery[END_REF].

A second source of uncertainties is related to the estimation of the wave characteristics used to compute the water depth. Some authors have investigated the sensitivity of the relative error in depth (estimated using the dispersion relation) to relative errors in frequency f and wavenumber k [START_REF] Dalrymple | Determining depth from remotely sensed images[END_REF]. They showed that the sensitivity of the relative error in depth to relative errors in f and k increases with the depth and is greater to the relative errors in f than in k.

In the proposed method, the data used to estimate the bathymetry using this equation are wave celerity c and wavelength λ. Let us investigate the propagation of uncertainty in (λ, c) using the dispersion relation. Fig. 12(a) illustrates theoretical bathymetric variations as a function of (λ, c) values. To investigate uncertainty propagation, coefficients C C and C L relating relative bathymetric errors and relative errors to the wavelength and the celerity, respectively, are numerically computed as follows:

Δh h = C C (λ, c) Δc c + C L (λ, c) Δλ λ . ( 5
)
Fig. 12(b) and (c) shows the patterns of C C and C L , respectively. First, the relative bathymetric errors are more sensitive to the relative celerity error than to the relative wavelength error: |C C | is generally larger than |C L |. C C is characterized by a minimum value of 2, meaning that a relative celerity error of 50% leads to bathymetric errors of at least 100%. This minimum value corresponds to a wavelength that is much greater than the water depth (i.e., h/λ → 0). C C exponentially increases with wave celerity c (and therefore the water depth), decreases with the wavelength, and exponentially increases with ratio h/λ. For instance, at point (λ, c) = 100 m, 10 m/s), C C is equal to 2.86 so that an error of 2 m/s (i.e., 20%) contributes to the relative error in bathymetry by 57%. The minimum value of C L is equal to 0. This can be explained by the asymptote of the dispersion equation, where λ is much larger than h: In this case, water depth h is only controlled by wave celerity c. C L increases with wave celerity and decreases with the wavelength, and its absolute value exponentially increases with ratio h/λ. For instance, at point (λ, c) = (100 m, 10 m/s), C L is equal to -0.48 so that an error of 20 m (i.e., 20%) contributes to the relative error in water depth by 10%.

This investigation highlights that the CWB method is highly sensitive to the quality of the wave celerity estimate. First, even with such high sensitivity, in the [5-40 m] range, the estimated bathymetry is realistic, with a relative error of approximatively 25% [see Fig. 10(b)] compared with the range of the C C factor (generally 2-3 in this bathymetric range). Assuming that the wavelength estimation is perfect and assuming the linear wave theory, this implies that the relative error of the estimated celerity is approximatively 10%. Reducing the relative celerity error to 5% should lead to a significant reduction in the relative bathymetric error (by a factor of two). To reduce these celerity errors, it could be interesting to consider the current effect in the dispersion relation even for areas covered by slight currents. As another perspective on reducing celerity errors, it would be worthwhile to investigate the use of subpixel displacement methods (see [START_REF] De Michele | Direct measurement of ocean waves velocity field from a single SPOT-5 dataset[END_REF], [START_REF] Michel | Deformation due to the 17 August 1999 Izmit, Turkey, earthquake measured from SPOT images[END_REF], and [START_REF] Leprince | Automatic and precise orthorectification, coregistration, and subpixel correlation of satellite images: Application to ground deformation measurements[END_REF]). Still another perspective would involve applying the procedure to data from higher resolution sensors keeping acceptable time delay DT, resulting in better correlation efficiency for a given displacement value.

D. Limitations and Recommendations for Using the Proposed Method to Analyze Satellite Images

As shown in Table I, the CWB method relies on several parameters, with some depending on satellite image characteristics (DT, DX SH , DX SL ), others on the site (i.e., the wave characteristics at the location and date of image acquisition), and still others on the choice of method, e.g., the number of clouds and the number of pairs per cloud. Researchers have no control over the satellite image characteristics, whereas some tuning of site-specific parameters might be possible. In this paper, these parameters were selected mainly based on physical considerations (overall wavelengths and variograms) and computation time. Even without any sensitivity analysis, the method provides reasonable results, considering that it is based on only one image pair. This means that it is possible to obtain reliable bathymetry without any tuning. Of course, when possible, a sensitivity analysis should make it possible to obtain better bathymetry. However, it must be recalled that the initial objective was to estimate bathymetry where no data are available and, therefore, no reference against which to tune parameters.

As guidelines for using this approach in other areas and on other satellite image types, the following recommendations can be made. As a preliminary condition, the waves should be visible on the images (e.g., bathymetry cannot be retrieved from quiet wave conditions or during rain and storm events). Second, the satellite-specific parameters should enable the description of the waves, i.e., not too large DX SL (to enable proper detection of the wavelengths in the wavelet analysis), not too large DT (to avoid ambiguity in the wave displacement estimate), and not too small DT with respect to DX SH (to enable detection of displacement in the cross-correlation step). This means that, physically, for the procedure presented here, the image/sensor characteristics should be such that

[DX SL < λ/2 ; DX SH < 2 c DT ; DT < λ/2c]. ( 6
)
The first inequality defines the minimum wavelength that can be reliably estimated given the sensor characteristics; the two others define the acceptable celerity range. Equation (6) provides the critical limit of the application of the CWB method. For instance, as DX SL becomes small, the results will be reliable, but for DX SL ≥ λ/2, wavelength λ will not be detected at all. Therefore, if conditions (6) are not satisfied, the proposed method cannot be used. If they are satisfied, then the site-specific parameters should be chosen so that the subscenes of dimension L S × L S contain at least two crests and two troughs (which can be determined by visual inspection), and the grid resolution for wavelet analysis should be small enough to describe spatial wavelength variations properly.

As illustrated in Section III, SPOT-5 images satisfy the application constraints of the CWB method for typical wave characteristics encountered around the world (in open sea, waves have wavelengths and celerities of the same order of magnitude as those encountered on the study site). However, many other systems can potentially provide estimates of bathymetry using such an approach. Among others, the Pléiades system provides panchromatic (multispectral) images with a resolution of 0.7 m (2.8 m) characterized by a time lag of 0.16 s (Centre National d'Études Spatiales, personal communication). With such characteristics, conditions [START_REF] Boccia | Tracking of coastal swell fields in SAR images for sea depth retrieval: Application to ALOS L-band data[END_REF] are satisfied because the expected celerities are greater than 3-4 m/s. It is worth noting the fact that using combinations of stereo pairs instead of bands from the same acquisition would not satisfy the condition. For instance, the minimal DT between a Pléiades stereo pair of images is approximatively 15 s; this means that the celerity at a 100-m wavelength cannot be reliably detected if the celerity is greater than 3.3 m/s. The higher resolution of Pléiades (or equivalent, e.g., SPOT 6/7) will provide better precision in the correlation computation and better (λ, c) estimates for the shortest wavelengths (particularly in the 5-20-m range) by using smaller subscene sizes, providing a more complete representation of the dispersion curve.

In an effort to improve the recommendations, an in-depth sensitivity analysis with respect to wave conditions, image quality, and parameter values should be done. It would be worthwhile to do it on a synthetic case using hydrodynamics that can simulate wave dynamics and provides instantaneous images of the wave field. This initiative falls outside the scope of this paper, which was to test the feasibility of estimating bathymetry from single-satellite image pairs, under real conditions, without any specific parameter tuning.

Finally, the CWB method has been developed assuming that currents were equal to 0. Such a hypothesis was justified on the application site that is characterized by surface current weaker than 0.14 m/s. To extend the applicability of the CWB method, it would be worthwhile to add currents as an unknown in the algorithm. The CWB method appears well suited for such integration as it is based on the study of several waves superimposing on each other. Then, as the manner in [START_REF] Piotrowski | Accuracy of bathymetry and current retrievals from airborne optical time-series imaging of shoaling waves[END_REF], once the water depth is estimated with the selected dominant wave and, if it exists, the contribution of the current is in that direction, we could use the other waves to estimate the current in other directions. This would provide better estimates of bathymetry in areas subject to significant currents, thus extending the field of application, but would also yield estimates of the currents themselves.

V. CONCLUSION

This paper has investigated whether a single SPOT-5 dataset can be used to retrieve bathymetry from space with reasonable accuracy and precision and without any in situ measurements. To tackle this issue, a new method (called CWB) was developed based on wave celerity and wavelength estimates. This method jointly uses the following: 1) wavelet analysis and cross correlation to provide (λ, c) clouds (one per detected dominant wavelength); 2) the linear dispersion relation to estimate one bathymetric value per cloud; and 3) criteria to select the "best" estimated bathymetric value at each location. The application of the CWB method to the Saint-Pierre coastal area (La Réunion Island) using one SPOT-5 panchromatic image and the corresponding green multispectral image showed that a single SPOT-5 dataset could be used to retrieve bathymetry from space with reasonable accuracy and precision (an absolute relative error of approximatively 20%-30% for water depths ranging from 3 to 80 m, with local areas providing errors fairly close to 0%).

A sensitivity analysis of the choices made in the CWB method illustrated the improvements brought about by using (λ, c) couples rather than (λ) assuming a uniform period to compute the bathymetry, by using the real part of wavelet analysis rather than Fourier-filtered scenes to compute (c), by considering clouds rather than a single (λ, c) pair to estimate the bathymetric values h n , and by choosing the "best" water depth ĥ among the N bathymetric values h n rather than using all couples (λ, c) to estimate a single bathymetric value. Analysis of the sensitivity of the relative bathymetric errors to the relative errors in the celerity and the wavelength highlighted the strong sensitivity to celerity errors and the potential for improvement by reducing the error in the celerity estimate. Some guidelines and recommendations have been provided for using the proposed approach on other sites and images. Other systems with higher resolution, such as Pléiades, could provide images, making it possible to extend the measured (λ, c) domain toward shorter wavelengths, thus reducing bathymetric errors.

In comparison with more traditional approaches based on LIDAR, multibeam echo sounders, or high-frequency time series of ground-based images, the results are coarse at the moment. However, they can be certainly improved in the near future. Being able to estimate bathymetry from space without any in situ measurements opens new perspectives, such as accessing remote and past bathymetries in many areas around the world, completing existing catalogues (e.g., extending the bathymetric catalogues further offshore or very close to the coast), estimating bathymetry before and between existing in situ acquisitions, and implementing a multitemporal bathymetric dataset.

Manuscript received August 5 ,

 5 2014; revised December 10, 2014, July 2, 2015, and October 19, 2015; accepted November 1, 2015. This work was supported in part by the French National Research Agency (ANR) Carnot Institution, by the French Geological Survey (BRGM), and by the National Center for Scientific Research (CNRS)/National Institute for Earth Sciences and Astronomy (INSU) through Programme National de Télédétection Spatiale (PNTS).

Fig. 1 .

 1 Fig. 1. Flowchart of the CWB methodology for retrieving bathymetry from two satellite images (SH for the higher resolution image and SL for the lower resolution image).

Fig. 2 .

 2 Fig. 2. Point clouds and nonlinear curve fitting corresponding to the n modes [n = 1 (black), n = 2 (blue), n = 3 (red), and n = 4 (green)] at the P1 location [the white cross in Fig.3(b)]. The selected ĥ is indicated with a thicker line (here, it is h 1 ). The theoretical curve obtained using the in situ bathymetry is shown by a dotted line.

Fig. 3 .

 3 Fig. 3. SPOT-5 panchromatic scene. (a) La Réunion Island. (b) Area of Saint-Pierre.

Fig. 4 .

 4 Fig. 4. Normalized wave amplitude spectrum W x at the P1 location [the white cross in Fig. 3(b)]. The spectrum is expressed in terms of the wavenumber in the (kx, ky) plane using the following convention: kx = 2π sin θ/λ and kx = 2π cos θ/λ The dimensionless maximum spectral intensity is W x = 0.1, and the contour intervals are 0.01.

Fig. 5 .

 5 Fig. 5. Windows of 160 m × 160 m centered on (a) panchromatic and (b) multispectral subscenes (in grayscale) and on the (c) dimensionless real part of their respective wavelet transforms and (d) corresponding to the P1 location [the white cross in Fig. 3(b)]. λ 1,1 = 85.4 m and θ = 2.62 rad are used for the wavelet transform.

l=0

  j) ∈ [-(S -1), S -1] are the indexes associated with the displacement in pixels. Correlation computations were carried out with the xcorr2 function of Matlab. CC n,m is the cross-correlation matrix of size (2 × S -1, 2 × S -1). Re(W SH n,m ) (or Re(W SL n,m )) is the real part of the wavelet transform corresponding to wavelength λ n,m and oriented by θ n as applied to the SH (or SL) image. The maximum value of matrix CC n,m gives the displacement d n,m between the two images, corresponding to wavelength λ n,m and ori-ented by θ n . Then, celerity c n,m corresponding to wavelength λ n,m can be deduced by considering the temporal lag DT between the two scenes (c n,m = d n,m /DT). The quality of the correlation Q n,m is given by the maximum of CC n,m normalized by the product of S-1 k=0 S-1 l=0 Re(W SH n,m (k, l)) Re(W SL n,m (k, l)) 2 that represents the autocorrelations of the real parts of both wavelet transforms. Fig. 5 provides an illustration of the center part of a subscene of the image pair and of the real part of its wavelet transforms. Applying this procedure for every index m ∈ [1 : M ], a point cloud

Fig. 6 .

 6 Fig. 6. (a) In situ measured bathymetry (in meters) on an 80 m × 80 m grid realized by LIDAR (from Litto3D products) for depths less than 30 m and by a multibeam echo sounder (SHOM data) for deeper water. The isobaths are represented by black curves and are plotted each from 5 to 50 m. The vertical reference of the bathymetry is the hydrographic zero (SHOM convention). (b) Estimated bathymetry (in meters) on an 80 m × 80 m grid from the SPOT-5 dataset. The black curves represent the isobaths at each 5 m given by the in situ measurements. The black area near the coast is not covered by the method.

Fig. 7 .

 7 Fig. 7. Subscene dimension and resolution depending on the distance from the coast.

Fig. 8 .

 8 Fig. 8. Selection percentage of the modes. The modes are the maxima of spectral energy (with Mode 1 (i.e., n = 1) being the most energetic and Mode 5 (i.e., n = 5) being the least energetic).for variogram analysis). This means that the subscenes overlap to cover the entire scene with regularly spaced wavelength spectra W x (λ, θ) (20-m resolution). A representation of the corresponding wavelength spectrum W x (k x , k y ) for point P1 is given in Fig.4. In this case, the four following dominant waves were detected: 1)λ 1 = 102.5 m, θ 1 = 150 • ; 2) λ 2 = 122.9 m, θ 2 = 128 • ; 3) λ 3 = 49.4 m, θ 3 = 135 • ; and 4) λ 4 = 45.1 m, θ 4 = 165 • .For the maximum number of modes (i.e., clouds) for which bathymetric values were estimated, an N max value of 5 was selected. For the maximum number of (λ, c) pairs per cloud, a value of M = 16 was chosen. For threshold T C related to the energy maxima of the (λ n , c n ) couples, based on some preliminary tests, a value of 0.65 was selected. Fig.2provides an example of the clouds obtained for point P1 [location indicated in Fig.3(b)]. In this case, four dominant modes, and therefore four clouds, were identified. This provides a reasonable number of points to estimate bathymetry h n (see the illustration in Fig.2). For the last step of water depth selection among the N values, C C was set to 1 m/s, whereas depth range DZ was set to [1-100 m]. The thick line in Fig.2indicates the selected ĥ value for location P1, which was used for comparison with the dispersion curve obtained using the in situ water depth (obtained from LIDAR and multibeam measurements).As a comment on the relevant N max value, analysis of the selected water depth ĥ shows that, most of the time, the first or second mode was selected for the final bathymetry estimation (see Fig.8), which indicates that there is no need to analyze more than N max = 5 dominant waves. From the CWB selection process (one criterion and two thresholds), long waves were most often selected. Note that the longest waves are most impacted by the bottom[START_REF] Holman | cBathy: A robust algorithm for estimating nearshore bathymetry[END_REF].

Fig. 9 .

 9 Fig. 9. (a) Absolute and (b) relative errors in bathymetry (in meters and in percentage, respectively) on an 80 m × 80 m grid from the SPOT-5 dataset.The black curves represent the isobaths at each 5-50 m given by the in situ measurements. The black area near the coast is not covered by the method.

Fig. 10 .

 10 Fig. 10. (a) Water depths ĥ estimated from the satellite dataset versus the in situ measured water depths h in situ considering all the points of the 80 m × 80 m grid. The blue crosses represent the 100 percentile ranks computed from h in situ = 2 m. (b) Absolute ( * ) and relative ( ) errors in bathymetry (in meters and in percentage, respectively) for the 15 quantiles of the in situ measured depth.

  2 of 16.2% in the [2 m, 20 m] range for local characteristics such that λ ≈ 60 m, h 1 = 2 m, h 2 = 20 m, and dx = 200 m; therefore, R1 = 30, R2 = 3, and R3 = 0.3. Using these parameters for this study and considering an offshore characteristic value λ ≈ 150 m, the results were h 1 = 2 m, h 2 = 50 m, and dx = 500 m. Estimating ER 2 with these characteristics (h 1 , h 2 , dx) led to an error ER 2 of 14.9% in the [5 m, 50 m]

Fig. 11 .

 11 Fig. 11. Relative errors in bathymetry (in meters) on an 80 m × 80 m grid from the SPOT-5 dataset for the following cases: (a) assuming a uniform wave period of 10.7 s to compute the bathymetry, (b) only using Fourier-filtered scenes and dominant waves to compute the celerity, (c) only using dominant wave characteristics to compute the bathymetry, and (d) using all (λ, c) clouds together to compute a single bathymetric value. The black curves represent the isobaths at each 5-50 m given by the in situ measurements. The black area near the coast is not covered by the method.

Fig. 12 .

 12 Fig. 12. (a) Theoretical water depth (in meters). (b) Coefficient C C . (c) Coefficient C L . (d) Coefficients C C (black) and C L (red) versus water depth divided by the wavelength, with the linear dispersion being applicable for 0.02 < h/λ < 0.5.

  ) shows the bathymetric errors obtained when only considering dominant waves rather than clouds to compute each h n . This figure should be analyzed in comparison with Fig. 9(b), which provides the errors when clouds are used. Using only the dominant waves leads to larger errors, particularly in water depths greater than 40 m and in steep areas.

TABLE I PARAMETERS

 I TO BE ESTIMATED FOR METHOD APPLICATIONS. THE VALUES ARE THOSE CORRESPONDING TO THE SPOT-5 IMAGES IN THE SAINT-PIERRE AREA, FEBRUARY 6, 2010

POUPARDIN et al.: WATER DEPTH INVERSION FROM A SINGLE SPOT-5 DATASET

ACKNOWLEDGMENT

The authors would like to thank the following data providers: the Kalideos team (Centre National d'Études Spatiales) for the SPOT-5 dataset, and the Naval Hydrographic and Oceanographic Service (SHOM) and the French National Geographic Institute for the bathymetric data (Litto3D and the multibeam echo sounder). The authors would also like to thank R. Pedreros for his contribution, as well as the reviewers for their comments.

He is currently with the French Geological Survey (BRGM), Orléans, France. His research interests include interferometric synthetic aperture radar processing and subpixel offset techniques for Earth surface displacement measurements applied to seismotectonics, volcanoes, natural hazards, and ocean swell velocity field from space.

Daniel Raucoules (M'12-SM'15) received the Ph.D. degree from Université Montpellier II, Montpellier, France, in 1997 and the "Habilitation à Diriger des Recherches" degree from the University of Orléans, Orléans, France, in 2011.

Since 1999, he has been a Researcher with the French Geological Survey (BRGM), Orléans, France, in the domain of spaceborne remote sensing tools (mainly radar interferometry and image correlation) applied to the measurement of ground surface deformation. His research interests (carried out in the framework of French and international projects) include methodological developments on techniques and their application to thematics dealing with ground motion, seismotectonic deformation, fluid injection/extraction, and postmining.