Automatic bathymetry retrieval from SAR images

Bathymetry, the topography of the sea floor, is in high demand due to the increase in offshore constructions like wind parks. It is also an important dataset for climate change modelling, when sea level rises and changes in circulation currents are to be simulated. The retrieval of accurate bathymetry data is a cost-intensive task usually requiring a survey vessel charting the respective area. However, bathymetry can also be retrieved remotely using data from Earth observation satellites. The main point of this study is the development of a processor that allows the automatic derivation of gridded bathymetry information from spaceborne Synthetic Aperture Radar (SAR) data. Observations of sea state modifications in SAR images are used to derive the bathymetry in shelf areas using the shoaling effect, which causes wavelengths to become shorter when reaching shallower waters. The water depth is derived using the dispersion relation for surface water waves, which requires wavelength and wave period as input parameters. While the wavelength can be directly retrieved from the SAR image, for the peak period additional information and procedures are required, e.g. local measurements or complex SAR data. A method for automatically deriving the wave period for swell waves in SAR images was developed and tested in this paper. It uses depth data from public databases as initial values which are compared to derived depths iterating through possible peak periods along the calculation grid; the peak period resulting in a minimal root-mean-square deviation is then used for bathymetry calculation. The bathymetry derived from a TerraSAR-X acquisition of the Channel Islands is presented; the resulting peak wave period of 11.3 s fits well to nearby in situ measurement data.


Introduction
With the rapid development of big data processing techniques, global 3D models of the worldwide topography like recently generated by the TanDEM-X mission [13], are state-of-the-art technologies demanded by a huge community of users. Such 3D models cover the global landmasses with high accuracy and horizontal resolutions of several meters. A follow-on step of these developments is to enable the extension of this 3D elevation models to underwater areas. Oceans cover more than 70% of the globe, however, hidden under the water surface, the bathymetry is far from being studied as good as the elevation above water.
Global bathymetry databases like General Bathymetric Charts of the Oceans (GEBCO) [6] do provide data for all areas worldwide, however, this data is often only interpolated or sometimes based on lead line measurements dating back to the 19th century. Local depth variations like sandbanks and coral reefs on the one hand and temporal morphodynamical evolutions of seabed structures on the other hand can be significant in littoral zones. In the Wadden Sea along the Dutch, German, and Danish North Sea coast, the soft seabed topography can change relatively fast due to storms causing large amounts of sediment transport [9] so that official charts can be quite out of date and sea signs need to be replaced regularly. With the current technological developments, many areas in remote locations, e.g. in Australian and African waters, are now becoming interesting sites for offshore constructions, requiring up-to-date bathymetry data. Accurate knowledge of the water depth is also crucial for many other maritime applications such as coastal protection, dredging, fishing, and ship traffic. In addition, reliable bathymetric information is required for sea level, circulation current, and sea state simulations, which help to predict sea level rises as well as wave climate and flood 1 3 events and, hence, are crucial for human safety. Databases like GEBCO offer bathymetric on a 1∕2 NM grid. However, a series of technical applications require more accurate depth information. For example, obstacles like narrow banks or reefs may be not resolved in this dataset, but since they stop the propagation of waves, they must be implemented into sea state models. Hence, in contrast to its importance, the bathymetry is relatively poorly known on a global scale.
The European Marine Observation Data Network (EMODnet) [4] offers bathymetry data for European waters in medium resolution by collecting survey data from various participating institutions. However, several data gaps remain also within these often sailed waters. As an example, Fig. 1 shows the available data for the Channel Islands area from EMODnet; coloured areas are covered by recent ship surveys and included in the EMODnet database, but in the grey areas, present especially around the islands themselves, only low-resolution and most likely outdated GEBCO data is available. A reason for these data gaps are the high costs of acquiring bathymetric data. The operation of a survey vessel equipped with appropriate sounding equipment and a trained crew is very costly. For shallow coastal waters with minor turbidity, airplanes equipped with LiDAR devices offer an alternative, but these are not available at low costs, either.
With the emerging Copernicus Sentinel satellite fleet, an increasing amount of earth observation (EO) data is available which can be used to derive the bathymetry from space with different techniques. While the topography on land can be accurately measured with radar satellite constellations like TerraSAR-X/TanDEM-X [13], a full measurement of underwater topography cannot be accomplished by a single spaceborne EO technology today. However, different satellite technologies combined can cover many areas of the oceans. For coastal waters, the bathymetry can be derived from optical satellites based on sunlight reflection analysis of the sea bottom, accounting for chemical and physical characteristics of sea water and the sea floor (e.g. [11,23]). This method provides depths in shallow waters up to 20 m in preferably calm weather conditions (no suspension due to wave stress). For mapping very deep areas, altimetry data can be used (e.g. [24]). Altimeter sensors measure the slope of the sea surface, which is modified by underwater terrain structures affecting the local gravity. Since above-sea-level landmasses also influence gravity, this method is typically confined to several kilometers off the shore.
Synthetic Aperture Radar (SAR) sensors can be used for intermediate depths (e.g. [3]). Their technique of depth derivation, which is also used in this work, is based on the observation of hydrodynamic surface processes, which are influenced by the sea bottom and thus reflect the underwater topography. Depending on sea state properties, this SAR data based technique covers the areas between about 70 and 10 m water depth. The acquisition quality of long swell waves depends on the local wind conditions: for wind speeds ≳ 15 m/s , the local wind sea waves interfere with the swell waves and nonlinear distortions and wave breaking effects hinder the accurate estimation of the swell peak wavelength [19]. In shallow areas ≲ 10 m , effects like wave breaking and coastal wave reflection may also hinder accurate automatic derivations. The depths in the range of 20m up to 10m represent the domain where a synergy of data from optical and SAR sources is possible. Thus, the results derived from both sensors can complement each other [18]. Other than with SAR as used in this study, this technique can also be performed with marine radar, e.g. WAMOS [10], and optical data [17].
A different method for bathymetry derivation from SAR is based on investigating changes of the water surface radar reflection [1,8]. Underwater structures with steep slopes cause strong surface currents. These affect the roughness of the surface by altering the small-scale Bragg waves which are responsible for the radar backscatter of the sea surface. Thus, underwater bottom topography becomes visible in SAR images by modulating the image brightness [20]. Using this approach, one can only compute relative spatial changes of the underwater bottom topography. However, a complete accurate mapping using this approach can hardly be applied worldwide.
The use of EO data for a robust and fast generation of derived products, here the bathymetry, requires algorithms suitable to run in automatic data processing chains. This operational data processing is typically part of a Near Real Time (NRT) chain, where no operator interaction is required. Combined with regular global data acquisition as is the case for the Sentinel satellites, the bathymetry can then be retrieved with minimal effort. To accomplish this, we present here an improvement for the SAR bathymetry algorithm by an automatic derivation of the peak wave period, a parameter which previously had to be manually supplied after inspection of the scene.
In the following Sect. 2, the use of TerraSAR-X imagery for bathymetry derivation and the methods applied in this study are described: the automatic retrieval of peak wavelengths and of the wave peak period. Section 3 presents and analyzes the results for a scene in the Channel Islands area for which in situ wave period data was available. A concluding summary is then presented in Sect. 4.

Method
In this paper, data from the TerraSAR-X (TS-X) satellite is used. The X-band SAR satellite TS-X was launched in 2007 and operates from 514 km height at a sun-synchronous orbit with a ground speed of 7 km/s (15 orbits per day). The typical TS-X incidence angle range is between 20 • and 55 • [5]. The TS-X data used for this study are Multi-Look Ground Range Detected (MGD) standard products with a pixel spacing of 1.25 m for StripMap mode (resolution of 3 m).
The detection of sea surface phenomena by SAR is dependent on the availability of wind [16], which is present in oceanic scenarios most of the time. Even very low wind speeds, ≈ 2 m/s , are sufficient to create ripple waves, centimetre-sized distortions of the otherwise smooth water surface [21]. Without these waves in scales similar to the radar wavelength, the radar backscattering does not reflect swell waves. The characteristics of the radar echo averaged over a subscene can directly be related to wind speed [15] and local variation of intensity can be related to wave height by image spectra analysis [19].
When long waves and ripple waves are present, the long waves are visible on the SAR image as regular brightness modulations. Here, the term "long waves" means wave structures which have a length of several pixels in the image-the pixel spacing depends on the satellite and the applied acquisition mode-and visible wave fronts, significantly longer than the wavelength.
The ocean surface gravity waves are constantly moving; hence, the mechanisms of their SAR imaging consist of the iniƟal/modified peak period peak wavelength field dispersion relaƟon database coarse depths iniƟal depth RMSD analysis optimal peak wave period resulƟng bathymetry map min/max for peak period subscenes pre-filtering, FFT peak period iterative selection SAR scene reading, calibraƟon, subscenes selecƟon After initial pre-processing of the scene and splitting into small subscenes, the subscenes are filtered and the FFT spectra are created. From these image spectra, the respective peak wavelength is derived and minimum and maximum peak wave period for the iteration are defined. Each iteration results in an initial depth value for each subscene calculated with the dispersion equation. The root mean square deviation (RMSD) between the initial depth and depth from a database is calculated for every wave period. Finally, the optimal peak wave period with the minimum RMSD is selected and the resulting bathymetry map is created linear transformation of tilt and hydrodynamic modulation, as well as the non-linear effects. This leads, among other effects, to image smearing and to a loss of information beyond the so-called azimuth cut-off wavelength. For TS-X images this cut-off is about 30 m for range traveling swell waves. The non-linear effects like the so-called "velocity bunching" [2] are lower for TS-X than for data from satellites like ENVISAT-ASAR or Sentinel-1 due to the lower TS-X orbit of 514 km , which allows a more stable imaging of the ocean surface (see also [19]).

Bathymetry derived from wave interaction
Radar waves practically do not penetrate the water surface and, hence, cannot directly, i.e. via runtime measurements, determine the height of the ocean floor. Similarly, interferometric observations as done by the TanDEM-X mission are not possible underwater. On the other hand, the advantage is that radar beams penetrate clouds and deliver an unobscured view at the sea surface in every acquisition. SAR bathymetry requires long swell waves because of the shoaling effect, which causes changes in wavelength and wave height when the underlying topography changes [18]. The effect occurs in medium-depth waters, when the waves interact with the seafloor. This requires, approximately, a water depth of less than half of one wavelength. Swell waves do not originate from local wind, but are a result of distant wind fields and can travel several thousand kilometres across the ocean with wavelengths of up to several hundred meters. However, swell waves approaching the coast can be superimposed by waves created under the local wind. A swell system usually has an extension of many kilometres parallel and perpendicular to the wave direction and becomes more homogenous in wavelength with increasing distance from its origin, as recently also shown using TS-X data [7]. Hence, for a SAR acquisition of about 30 × 50 km (TS-X StripMap scene), incoming swell in deep water may typically be considered constant in the wave period and wavelength.
Using a linear approach and neglecting local circulation currents, the connection between peak wavelengths and depth is described by the dispersion relation in the form where p = 2 ∕T p is the peak wave frequency, T p is the peak wave period, L p is the peak wavelength and g is the gravitational acceleration. This can be solved for the depth d to From Eq. (2), it follows that two parameters are required to calculate the water depth: the peak wavelength L p and the peak wave period T p . The determination of these parameters is described in the following sections. Figure 2 gives an overview of the algorithm workflow.

Determination of the wavelength
The wavelength is retrieved via a two-dimensional Fast Fourier Transformation (FFT) of individual subscenes as explained in [18,19]. However, in contrast to [18], not the wave ray technique, but a gridded approach was applied which is more suitable for automatic processing. A requirement for the successful derivation of the wavelength is that the subscene contains several visible wave structures. Subscenes on land or mostly on land, which is decided by a static landmask or an automatically generated landmask for the respective scene [26], are excluded from the further analysis.
Since the goal of the algorithm is automatic operation, it is important that possible errors in the spectrum are accounted for. The sources of these error are in the first place a number of natural features and man-made objects typical for coastal waters: sandbanks and small islands, wave breaking, ships, current boundaries, and also atmospheric fronts. Even internal wave structures impact the image spectra. Such spectral perturbations result in an integrated value which yields a contribution to the total spectral energy that is not connected to the sea state. In order to exclude all these Fig. 3 Overview of the acquisition area at the western exit of the English Channel. The box indicates the location of the TerraSAR-X scene, the circle indicates the position of the Channel Lightvessel errors, the automatic pre-filtering procedure introduced in [19] is applied before the FFT analysis for each analyzed subscene.
Another source of error can occur by the presence of more than one swell system in the scene. For each subscene, only the dominant peak wavelength is retrieved. If the dominant swell field changes in a subscene, strong changes in wavelength and direction occur compared to neighbouring subscenes. These are avoided by applying a filtering routine after initial processing, where in this case the non-dominant peak is identified by limiting the spectrum search to wavelength and direction similar to those of the neighbours.

Determination of the wave period
In contrast to many other maritime parameters, there are only limited ways to derive the peak wave period from SAR acquisitions. Since SAR image are acquired over a certain time span, usually one second or more, it is possible to split a scene into multiple time frames and, by analyzing the differences, observe the wave period [25]. However, this works best for SpotLight scenes where the radar beam is looking at a small area for an extended time, and is less applicable to StripMap or ScanSAR modes. Due to their wide coverage, the latter are more suitable for bathymetry retrieval in extended areas, which is why they are preferable for this type of study. Information about the wave period can be provided by in situ measurements, for example from wave-rider buoys if they are inside or nearby the observed area. Unfortunately, this is rarely the case. Additionally, for automatic operation the respective data would need to be available in a public online database which can be queried, a feature many of the measurement devices do not support.
Another approach requires initial, first-guess information about the depth, which may be derived from sea charts or bathymetric databases. Then, using Eq. (1), the wave period can be calculated in an area with known depth, an approach used in [18]. Brusch et al. [3] used a similar approach, but averaged expected deviations due to the initial depth derivation. Both require having an operator looking at a single area of the scene to set an initial bathymetry value.
In this work, a new approach was used to further automate the bathymetry retrieval algorithm: the derived bathymetry is compared to external bathymetric data on a point-by-point basis. For this, data from freely available sources like GEBCO (available worldwide) or EMODnet (available around Europe only) is used. This external data constitutes an initial depth information. The depth value at the center of each FFT scene is retrieved. Doing this, multiple solutions for the bathymetry are calculated by iterating through a range of possible wave periods. This range depends on the wavelengths of the specific scene: longer wavelengths are connected to longer wave periods and vice versa. Since tanh −1 2 L gT 2 ≤ |1| , a lower limit of the possible wave period exists for each wavelength. For example, a wave of L p = 200 m requires T p ≳ 11.3 s while a wave of L p = 150 m requires only T p ≳ 9.8 s (see also [18], Fig. 14).
An entire satellite scene usually consists of three types of areas relevant for the wave period determination: land and near-land areas with wave breaking and wave reflection effects (1), deep water where the wavelength does not change (2) and intermediate depth water (3). The wave properties in type 1 (land/near-land) cannot be reliably derived since the distortions here are too strong to be filtered. Thus, the land masking applied for the wavelength derivation extends about 1 km into the sea, so these areas are completely filtered out. In type 2, the deep water areas, no interaction between waves and seafloor takes place. The decision if an area is deep water or not is made by comparing the wavelength retrieved for the respective subscene with the database depth for this area. Areas where the water is deeper than half of the wavelength are considered deep water. While the peak period can be derived using the deep water dispersion relation, no information about the bathymetry can be retrieved from it.
The remaining areas, type 3, constitute water with a depth that influences the wavelength. The respective subscenes are used for the calculation of the wave period as follows: First, the longest wavelength found in all the subscenes is used to calculate the minimum wave period. Since the wavelengths get shorter when approaching shallow waters, this ensures the long incoming swell waves are used to determine the lower boundary of the iteration range. In addition, using the longest wavelength ensures the resulting wave period is large enough to calculate the bathymetry for all wavelengths in the non-deep water area. With the minimum wave period calculated as described above, the depth in each of the subscenes is calculated using Eq. (2). This yields a depth value for each subscene, which is combined with the database depth to perform a root mean square deviation (RMSD) analysis: where E i is the deviation for the wave period index i, d is the calculated depth and d ′ the database depth, the index j denotes the individual subscene. After the iteration through all wave periods is complete, the wave period resulting in the least deviation is selected and the final depths are calculated using this wave period.

Results
The presented method for automatic wave period retrieval was successfully applied for a series of locations including northern Mauritius, the Bay of Biskay, Australia and also the Channel Islands, which are selected as a case study for this paper. This region at the southern exit of the English Channel is interesting due to the extended areas of 30-50 m depth and the frequent exposure to swell waves originating from the Atlantic Ocean in the west. This combination of intermediate water depths and long swell waves is well-suited for our applied method of SAR bathymetry. For comparisons of the wave period, in situ measurements from the Channel Lightvessel situated about 56 km north of Jersey is used. The input bathymetry data was taken from the EMODnet bathymetry data portal (accessed on 2017-03-06). Figure 3 shows the location of the scene and of the Channel Lightvessel. The investigated TerraSAR-X scene was acquired over the island of Jersey extending up to the French coast on 2010-03-31 at 06:18 UTC in StripMap mode covering an area of 30 km × 50 km . For analysis, the image was divided into a grid of 24 by 44 equidistant subscenes with an extension of 2.5 km × 2.5 km each. The distance between two subscenes is only 1.25 km ; hence, a 50% overlap is used. This was set to increase the resolution of the retrieved bathymetry, while having a subscene size sufficiently large for wavelength retrieval.
The automatic derivation of the wave period, performed as described in Sect. 2.3, yielded a peak wave period of 11.3 s . The measurement of the Channel Lightvessel gives a wave period of 8 s which was recorded at 07:00 UTC, 42 minutes after the satellite acquisition. This record constitutes the mean wave period averaged over 17.5 min. However, from the SAR image only the dominant swell waves and, hence, the peak swell period is retrieved. For comparing peak and mean wave period, [14] suggested a conversion factor of T peak = 1.41T mean . That means the reported mean wave period of 8 s equals a peak wave period of 11.28 s . Unfortunately, the mean wave period is only reported in full (points). The background shows the EMODnet data used as input for the automatic wave period retrieval. The RMSD for these datasets is 7.1 m seconds, so the possible uncertainty is ± 0.5 s . With this, the peak wave period is in the range of 10.52 s to 11.99 s . Our automatically retrieved value is, therefore, very well within the margin of error provided by the measurement. Figure 4 shows the results of the bathymetry calculation (points) compared to the EMODnet bathymetry (background); Fig. 5 additionally provides a 3D look at the results. The bathymetry retrieved with our algorithm can well reproduce the general features present in this scene, like the deeper areas north-west of Jersey or the very shallow area in the immediate north-east of the island.
The RMSD calculated for this scene and the chosen wave period of 11.3 s is 7.1 m . There are several factors that may influence the resulting accuracy. On the one hand, much of the EMODnet data around Jersey is not yet retrieved from recent survey data but taken from the GEBCO database, as shown in Fig. 1. This data is of relatively coarse resolution and might be rather outdated. Therefore, while EMODnet and GEBCO data certainly provide a sufficiently realistic first guess for bathymetry data, it cannot be ensured that the dataset available for comparison actually depicts the true bathymetry at the acquisition time.
On the other hand, the English Channel and, hence, also the area around the Channel Islands, is subject to strong current flows which were not taken into account in the study. These influence the wavelengths and, consequently, also the calculated bathymetry. The effects of these currents can be implemented by adding the current term Uk j to Eq. (1), where U is the current velocity in j-direction [12]. However, as the selected scene contains areas also behind islands and close to the coast, a single value for the current speed would not be sufficient and a high-resolution current dataset for the acquisition time would be required for this correction. In addition, the investigated scene shows no indications of strong currents which are typically visible in the SAR image by bright or dark lines at the current boundaries, which are present in several other acquisitions of this area. Another deviation is caused by different water levels: since the water height above the sea floor changes during a tidal cycle, the water level at acquisition time is usually different from the lowest astronomical tide (LAT) water level given in EMODnet or often found in nautical charts. Since the area of the Channel Islands is rather open, the tidal deviations can lead to a minor water level gradient across the scene, but will mostly cause an absolute offset for the entire area.
Additionally, in gridded bathymetry data usually the average depth within a grid cell is given. The lowest depth in the respective cell may be different from the average value, but is essential for the wave interaction. For example, small structures like sheet pilings or single rocks are not incorporated into bathymetry, but are important for nautical charts and can be detected with the presented method.
In a simulation study on bathymetry obtained from FFTderived wavelengths, the same method as used here, Shen [22] found errors reaching from 8% to 18% of the depth, depending on the complexity of the underlying bathymetry and the absolute water depth. Table 1 lists the RMSD separated by derived water depths and their respective wavelengths for the retrieved peak wave period of T p = 11.3 s . The deviation increases for larger depths. Parallely, the range of possible wavelengths to yield these depths decreases. Hence, inaccuracies of the wavelength retrieval have a larger effect for longer waves (or deeper waters) than for shorter waves (or shallower waters). Figure 6 shows the effect of an over-or underestimation of the peak wave period on the retrieved depths. By varying the peak period by ± 0.5 s from its determined optimal value of T p = 11.3 s , the RMSE is increased by 2.4 m for T p = 10.8 s and by 1.1 m for T p = 11.8 s . Hence, finding the optimal peak wave period is a crucial step to retrieve the most accurate depths.

Conclusions
A practical approach for exploring the near-coast bathymetry from about −70 m up to −10 m depth using remote sensing data from high-resolution SAR images presented in [18] has been extended by automatization of peak period derivation, enhancing the applicability and accuracy of SAR bathymetry retrieval. The study shows the possibility for an automatic improvement of coarser water depth data, e.g. GEBCObathymetry, using satellite-based SAR. To create a more extended bathymetry dataset, this method can be combined with other spaceborne data sources like bathymetry derived from optical satellite data or altimetry data.
As shown in the presented scene, the algorithm determined a wave period consistent to the period measured by a nearby lightship after the mean wave period was converted to peak wave period. This prooves that the suggested method worked well for the investigated scene and suggests applicability also for the majority of scenes where wave period data is not available, as it is mostly the case. The bathymetry retrieved from the TerraSAR-X scene fitted well to the data from EMODnet for the Channel Islands area. The retrieved RMSD of 7.1 m is considered acceptable regarding the uncertainties contained in the EMODnet dataset itself,  where a large part of the acquisition area is only covered by GEBCO data, and also regarding the current and tide situations in the English Channel. The developments presented here enable the automatic application of SAR bathymetry derivation, a technique possible in many shelf areas worldwide. The identification of scenes containing swell is currently done by manual visual inspection. However, the wavelength retrieval is fully automated and not dependent on specific scene conditions or swell waves [19]. In scenes without swell waves, when only short wind waves are present, no individual wave fronts can be identified in the SAR image and wavelengths cannot be retrieved. When the algorithm is part of an operational service, bathymetry information will only be calculated for scenes and subscenes containing wave patterns with wavelengths suitable for bathymetry retrieval.
The principles of the algorithm also work on optical satellite data, therefore, optical and SAR data may complement each other in areas not frequently covered by clouds. This will help to improve the availability of up-to-date bathymetric data in many coastal regions of the world.