Methods for Satellite Derived Bathymetry from Sentinel-2 images: a comparison

: In recent decades, research has been developed to estimate near-shore bathymetry depth values using satellite imagery. Visible and infrared bands are used to derive elevation profile estimates, so to obtain bathymetric in rapid way without mobilisation of persons or equipment and saving the costs. For consequence, Satellite Derived Bathymetry (SDB) is seen as a valid approach for shallow waters survey: strongly supported by the activity of scholars and researchers, multiple methods are available in the literature. This article aims to investigate and compare different SDB methods for sea depth extraction from Sentinel-2 satellite multispectral images, with particular attention to the accuracy of the results. The experiments are conducted on imagery including Blue, Green, Red and Near Infrared bands, with 10 m resolution, concerning the Bay of Pozzuoli (Italy). After removing the glint, the effects caused by the reflection of sunlight through single scattering from sea surface, three methods are applied: Band Ratio method (BRM), 3rd-degree polynomial regression line method (3DPM), and principal component analysis method (PCAM). 3DPM can be seen as a variant of the BRM where the linear law that interprets the correlation between the band ratio values and the depth values is replaced by the third order function. Models are trained using depth data extracted from an Electronic Navigational Chart (ENC) at 1:7,500 scale, which is also used to verify result accuracy. The experiments demonstrate that the 3DPM is better able to obtain a more precise bathymetric model, confirming the greater adaptability of the 3rd order function to interpretate the variability of the interaction of light with water along the water column.


INTRODUCTION
The seabed is about 70% of the total surface of the planet earth, and 50% of the global continental shelf area (extension of the first 200 m) is still undetected (Duplančić Leder et al., 2020).The reason is that historically the surveys were done through an acquisition campaign on hydrographic vessels equipped with special instrumentation such as echosounder, with both financial and time costs.The measurement of depth is therefore of primary importance because it also allows the study of fundamental phenomena related to coastal erosion (Thompson et al., 2015;Tsiakos and Chalkias, 2023), flooding (Karamouz et al., 2021), hydrodynamic and wave modelling (Belibassakis et al., 2020).Traditional bathymetric surveys are made with singlebeam (SBES) or multi-beam (MBES) (Alcaras, 2022b).In general, these systems sending sound waves into water; the time it takes for the sound to cross the water, i.e., to reach the bottom and return to the transductor, determines the depth of the water (Amoroso and Parente, 2021).The SBES sends a single pulse per measurement, while the MBES sends multiple pulses and reads all scattered data.This allows MBES to have a larger swath width and area coverage (Zirek and Sunar, 2014).These systems are very accurate but as previously mentioned the main disadvantage is that they are expensive both in terms of time and cost.Consequently, even though echo sounders have been used for 5 decades, a large part of the ocean floor remains unexplored: the main reason for this lack of data is that bathymetric surveys carried out by ships are slow and expensive; just think that it would take more than 120 years to systematically map the sea depths using ships (Sandwell et al, 2006).Considering the above-mentioned disadvantages of echosounder use, it has been necessary in recent times to find an alternative that overcomes the problems of traditional bathymetric surveys.The techniques for acquiring the depths of the seas have gradually moved from the surface of the water passing from the air to space.Remote sensing therefore makes it possible to obtain bathymetric data, and one of the most popular instruments in recent times is the Airborne Laser Bathymetry (ALB).These instruments operate in the green wavelength domain, allowing the laser to penetrate the water column and map the seabed (Doneus et al., 2015) (Wang and Philpot, 2007).The acquisition tool is mounted on an aerial platform or drone, permitting to reach even areas that are not easily accessible.In general, it mainly consists of a GNSS receiver (for the position of the aircraft), a laser emitter and an optical receiver.This technology has high accuracy and wide coverage, but the cost of measurement is still high, and the flight conditions of the sensing area may limit its application (Lewicka et al., 2022).Inside remote sensing, specifically in the sector of earth observation from the space, we can distinguish Satellite Derived Bathymetry (SDB) techniques in three main groups: Synthetic Aperture Radar (SAR) techniques, altimeter techniques and optical techniques.SAR techniques used for depth derivation are generally based on the observation of surface hydrodynamic processes, which are sensitive to the bathymetry and therefore reflect the submarine topography (Wiehle and Pleskachevsky, 2018) (Bian et al., 2020).In other words, these techniques are based on the interference between the electromagnetic microwaves signal coming from the satellite with the sea surface roughness that is influenced by the seabed morphology (Cloarec et al., 2016).The obtainable depths range from -10 m to -70 m, therefore the first meters (> -10 m) cannot be achieved due to effects such as wave breaking and coastal waves reflection (Wiehle et al., 2019).This is a huge limitation for shallow water derivation also accompanied by a low geometric resolution in case of free data which can go up to 100 m (Duplančić Leder et al., 2023).The altimeter techniques for SDB are based on altimeter satellite data and the strong correlation between the ocean bottom topography and the shape of the geoid, a gravitational equipotential surface defined by mean sea level.The depth variations of the seafloor can be observed as height variations of mass elements of the density Δρ which is the contrast between the density of the seabed ρsb and seawater ρw (Calmant and Baudry, 1996).Even if in the 1970s a number of studies in the ocean basins had indicated that the gravitational field or geoid is correlated with bathymetry at short to medium wavelengths (30-300 km), Dixon et al. (Dixon et al., 1983), using SEASAT (Evans et al, 2005) altimeter data on a test area in the Musician Seamounts province north of Hawaii, demonstrated that seafloor topography can be predicted from marine gravity anomalies.Smith and Sandwell (Smith and Sandwell, 1994) modelled the bathymetry from marine gravity anomalies integrating data from SEASAT, ERS-1 (Duchossois, 1991), and Geosat (Jensen and Wooldridge, 1987) geodetic mission in the Southern Ocean; they used ship soundings for validation.These techniques present a strong limitation since the map they produce do not have sufficient accuracy and resolution to be used for assessing navigational hazards; nevertheless, they are useful for such diverse applications as locating obstructions to the major ocean currents and shallow seamounts (Jawak et al, 2015).The optical techniques for SDB allow the bathymetry to be obtained for shallow waters (depth up to 20-25 m) and for large areas in a rapid and increasingly precise way, reducing acquisition times and costs (which can be zero) (Apicella et al., 2023) (Tursi et al., 2023).An online questionnaire carried out in 2021, as part of the European Commission's "4S" innovation project on the use of shallow water surveying techniques, highlighted that 28% of 213 stakeholders recognized the ability to map inaccessible places as the main reason to use SDB technology.The survey also showed other interesting statistical data that underline the importance attributed to this new approach to survey shallow waters, e.g., decreased efforts for organization and management of data acquisition operations (18%) (Hartmann et al, 2022).The optical techniques for SDB are based on the principle that radiance is, to varying degrees, attenuated by the water column; the degree of attenuation coefficient depends on different factors, i.e., band wavelength, sea bottom types, and water column properties (Gholamalifard et al., 2013) Even if the use of satellite images is intended to be an alternative to in situ surveying, bathymetry measurements are to be preferred to enhance model training and validation (Çelik, 2023).Nevertheless, the presence of archive data from previous surveys or nautical maps with accuracy adequate for the resolution of the images, can be sufficient for both training and validation of the model (Figliomeni and Parente, 2023).In the literature there are different methods that allow to obtain the depth of the seas from optical satellite images, such as Band Ratio method (BRM) (Stumpf et al., 2003), 3rd-degree polynomial regression line method (3DPM) (Figliomeni and Parente, 2023) and principal component analysis method (PCAM) (Hashim et al., 2020).Since they provide a depth value for each pixel, these methods for SDB return continuous models of the bathymetry and do not require interpolation processes.In this way, another source of error connected to the capability of the algorithm to provide reliable depth values starting from other measurements, is avoided.In this work we compare the above mentioned three SDB methods for sea depth extraction from satellite multispectral images, analysing and discussing with particular attention the accuracy of the results.The experiments are conducted on Sentinel-2 images (Blue, Green, Red and Near Infrared) with 10 m resolution concerning the eastern sector of the Bay of Pozzuoli (Italy).After removing the glint, the effects caused by the reflection of sunlight through single scattering from sea surface, BRM, 3DPM and PCAM are applied.Models are trained using depth data extracted from an Electronic Navigational Chart (ENC) at 1:7,500 scale, which is also used to verify result accuracy.All operations are performed using the free and opensource GIS (Geographic Information System) software named Quantum GIS (QGIS), version 3.22.The document is structured as follows.Section 2 describes the study area, the datasets used and the methodologies for carrying out the SDB.Section 3 presents and discusses the experiment results.Finally, section 4 completes the document with the conclusion.

Study area and dataset
The Bay of Pozzuoli is located in the western area of the city of Naples, in southern Italy, and extends from Capo Miseno to Capo Posillipo.Figure 1 shows the geo-localization of the study area reporting the territorial framework in equirectangular projection and WGS 84 ellipsoidal coordinates (EPSG: 4326) in the upper part; a satellite overview of the Bay of Pozzuoli (true colour RGB composition of Landsat 8 OLI images) in UTM/WGS84 plane coordinates expressed in meters (EPSG: 32632) is shown in the lower part of the figure 1.The Bay of Pozzuoli, included within the Gulf of Naples, is located near the city from which it takes its name, Pozzuoli, the Greek colony of Dicaearchia founded in about 531 BC in Magna Graecia, conquered by the Roman in 195 BC and renamed Puteoli whose roots are in the Latin "puteus" (well or cistern).The geomorphology of this bay has been significantly influenced by the presence of important active volcanoes, such as Monte Nuovo and Solfatara, and by other geological phenomena such as land subsidence (Somma et al., 2016).Furthermore, human activities have greatly altered the morphology of the entire bay.The coast of Bagnoli and the part of the sea in front, in fact, have undergone a significant deposition of marine sediments, given by industrial activities (Di Martino et al., 2020).The coastal area is characterized by rocky cliffs and pebble beaches, while the interior of the bay is home to numerous inlets and lagoons.For the study conducted, images from the Sentinel-2 mission are used, acquired on September 22 nd , 2022.The area investigated has an extension of 2.4 Km x 2.4 Km, it is enclosed within the following UTM/WGS84 -zone 33N coordinates: E1= 423,600 m, E2= 426,000 m, N1= 4,521,300 m, N2= 4,518,900 m.This area is shown in Figure 2

Workflow
The workflow for generating the bathymetric model obtained via satellite image involves several steps.
Once the satellite images have been downloaded from the appropriate portal, we start with a pre-processing phase, in which the water body is separated from the land area.The noise caused by the glint is then removed.When the pre-processing phase is completed, we proceed with the application of the SDB methods.Depth data extracted from the ENC are used to train the model.In this way, Digital Bathymetry Models are obtained, which are finally tested with available bathymetric data to verify the accuracy level.

Preprocessing
To apply the SDB methods, it is necessary to perform image pre-processing.The first thing to do is to separate the sea zone of the study area from the land one.For this reason, the coastline detection in the available dataset is necessary: the method proposed by Alcaras et al. (Alcaras et al., 2022), which automatically allows to separate the water and no-water areas is chosen.This method involves the use of the water index called Normalized Difference Water Index (NDWI), introduced by McFeeters (McFeeters, 1996).In the case of Sentinel-2 this index requires the green band, i.e. band 3, and the NIR band, i.e. band 5, in accordance with the following formula: Subsequently there is the application of the unsupervised classification algorithm called K-means proposed by MacQueen (MacQueen, 1965).K-means allows to describe a process for partitioning an N-dimensional population in k sets with a high degree of similarity within each set and a low degree of similarity between sets (Friedman et al., 2009).Starting from the operator's choice of the number of classes, which in this case are two, k-means permits to obtain the water and no-water thematic map.Finally, through automatic vectorization, two large polygons are obtained which can be chosen as a mask to separate the affected water area (Figliomeni et al., 2023).Once this is done, the removal of the glint is applied, through the use of the methodology proposed by Hedley et al. (Hedley et al., 2007).Glint is the noise generated by the fact that sunlight is partly reflected by the sea surface, obtaining an additional contribution of reflectance values.The satellite sensor therefore acquires this quantity, distorting the final bathymetric model.The method proposed by Hedley involves the selection of sea points where the glint is not present, generally deepwater areas, determining the minimum NIR value.A regression line is established between the NIR values (x-axis) and the values of the band of which you want to correct the glint (y-axis), i.e. the visible bands used.The slope (bi) of this line is calculated.Finally, the following formula is applied to obtain the correct value of the visible band considered: Where; Ri is the initial value to be corrected; bi is the slope of the line; RNIR is the value in the NIR; MinNIR is the minimum value of the NIR band.
Finally, 4 cropped images are obtained in the interested stretch of sea, 3 of which are also corrected by the noise given by the glint.

Band Ratio Method
One of the most used methods in SDB was proposed by Stumpf et al. (Stumpf et al., 2003), called Band Ratio Method (BRM).This method, empirical type, linearizes the relationship between the spectral values, of two bands, and the depth of the seabed with a logarithmic transformation.
The bands used are blue and green, as by making better use of the intrinsic optical properties of water.In fact, in the blue band, clear water has a higher reflectance peak than in the green band; on the other hand, in coastal waters, where organic and inorganic elements are present, the green band has a higher reflectance peak than the blue one.These different spectral responses allow to obtain information for the determination of the bathymetry.Therefore, the depth (z) is calculated according to the following formula: Z= m1 [ln ( n * ρW(λi) ) / ln ( n * ρW(λj) )] -m0 Where: m1 is the constant to scale the depth ratio; n is the fixed constant to make the logarithm positive; m0 is the depth offset where z = 0; ρW is the reflectance of water; λi is blue band; λj is green band.

3rd-degree polynomial method
In 2023, an alternative method to the application of the BRM was proposed (Figliomeni and Parente, 2023), called the 3rd Degree Polynomial Method (3DPM).It exploits the relationship between the blue and green bands but fitting the depth data to the spectral data with a third-order polynomial function.In fact, the propagation of light in water depends on various aspects relating to the interaction of sunlight with the marine environment, and therefore the linear relationship, given by the BRM, is not always able to better interpret the relationship between the spectral values and the real depth.The depth is then determined with the following formula: Where: n is the constant to make the logarithm positive; ρW is the reflectance of water; λi is blue band; λj is green band; a, b, c, and d are the constant coefficient (obtain after applying third polynomial regression).

Principal Component Analysis method
Principal component analysis is used to transform a data set, such as bands, into new ones (called component images) that are uncorrelated with each other.The new bands are organized such that each has lower variance than the previous one.Thus, the first two component images have more than 90% of the variance (Mudrová and Procházka, 2005).
For SDB only the first three bands, corresponding to the blue, green, and red, are transformed with the Principal Component Analysis Method (PCAM).This method preserves the uncertainty caused by depth differences and ignoring small gap due to local shifts in turbidity or bottom shape (Hashim et al., 2020).However, different approaches are available in literature on PCA application for SDB.One of them assumes that the first component resulting from PCA application is the most performing for water depth determination; in other terms, it is the first component that can be calibrated to known water depths (Gholamalifard et al, 2013).Other applications use three principal components.After the transformation, three new bands are obtained and the multiple regression is applied between the of the three resulting images and the depth values, obtaining the following transformation equation: Where: x1, x2, x3 are respectively the first, second and third component image after PCA transformation; a0, a1, a2, a3, are constant value (obtained after applying multiple regression).

Accuracy evaluation
Starting from the three Digital Bathymetry Models (DBM) obtained from the application of the previous described methods, we proceed with the evaluation of their accuracy.To evaluate the effectiveness of the results, the depth data of the ENC are considered as in situ data.In particular, the difference between the true depth values and the respective ones supplied by each DBM is made, thus obtaining the residuals.
In order to make a better evaluation, we provide the statistical values (mean, standard deviation, Root Mean Square Error (RMSE), maximum and minimum) for each model.Particularly, RMSE is calculated using the following formula: Where: N is the length of the dataset (number of considered point values); Vri is real (measured) depth value; Vpi is the predict depth value.

RESULTS AND DISCUSSION
The following figures show the bathymetric models resulting from the application of the previously described methods.
Particularly, the figure 4 shows the BRM derived model, the figure 5 the 3DPM derived model, and the figure 6 the PCAM derived model.Table I shows the statistical parameters that characterize the accuracy evaluation process.In particular, these parameters are calculated for each model considering the residuals between the true depth data extracted from the ENC and the corresponding satellite derived values.The obtained values can be compared with the pixel size, i.e. 10 m, to establish the accuracy of the generated models.In general, the results show that the methods used provide excellent DBM as the RMSE is below the pixel size.Among these methods, the BRM turns out to be the worst, in fact it has the highest RMSE value, i.e., 2.934 m, with a minimum value of -8.252 m, and a maximum value of 10.047 m.The 3DPM is the most performing method with RMSE value of 1.950 m, and maximum (6.736 m) and minimum (-6.842 m) values below the pixel size.Finally, the multiple regression method with PCAM is positioned in a middle position, with RMSE value above 2 metres.

Statistics Values
Figures 7 displays the 3D DBMs generated by the applied methods: BRM (upper), 3DPM (in the middle part), and PCAM (lower).
In order to evaluate the validity of the results, other works present in the literature are examined.Said et al. (Said et al., 2018) tested two SDB methods, the BRM and the Lyzenga method (LM) (Lyzenga, 1981) (Lyzenga, 1985), on Sentinel-2 achieving RMSE equal to 2.347 m for the BRM and 2.474 for LM.

CONCLUSION
The experiments carried out in this work involves the application of three SDB methods on the Sentinel-2 dataset, presenting 10 m resolution and concerning the study area of the Bay of Pozzuoli (Italy), in order to evaluate the most efficient approach to obtain a DBM.In particular, the methods used are BRM, 3DPM, and PCAM: the first two use only two bands (Blue and Green), the third three bands (Blue, Green and Red).
The Near infrared band together with the Green band allow the application of NDWI in order to distinguish water and no-water and select the sea from the context of the image.The methods are trained using depth data from an Electronic Navigational Chart (ENC) at 1:7,500 scale.The extracted depths are also used to evaluate the accuracy of the resulting models.In consideration of the applications, the most performing one among the considered methods is the 3DPM with a RMSE equal to 1.950 m.Subsequently, PCAM is second with RMSE equal to 2.603 m and finally BRM with RMSE equal to 2.934 m.
The experiments demonstrate that the 3DPM is better able to obtain a more precise bathymetric model, confirming the greater adaptability of the 3rd order function to approximate the correlation between the depth values and the variability of the interaction of light with water along water column.
However, the penetration of light along the water column is linked in various ways to the interaction of sunlight with the marine environment, therefore, the analysis of the specific situation is necessary to identify correct empirical modelling.Future developments of this work will focus on the application of SBD methods on a satellite imagery with greater geometric resolution.In addition, broader comparisons will be accomplished by integrating the three methods considered in this work with others that exploit Artificial Intelligence (AI) to determine a more performing model.

measurements.
The International Archives of the Photogrammetry, Remote Sensing and Spatial Information Sciences, 40, 221-227. https://doi.org/10.5194/isprsarchives-XL-7-221-2014 . Sentinel-2 images are acquired by the mission developed by ESA under the Copernicus programme.The mission consists of two identical satellites, in the same orbit but 180° out of phase from each other; the first satellite put into orbit on June 23, 2015, is called Sentinel-2A, while the second is Sentinel-2B put into orbit on June 23, 2017.The sensors on board the sentinel satellites are capable of acquiring 13 multispectral bands between the visible and near infrared.The spatial resolution is variable (10 m, 20 m, 60 m) as it depends on the type of band.The satellite flies over the same point on the earth's surface with the same angle of view every 5 days.The mission mainly supports studies in the fields of climate change, land monitoring, emergency management and security (European Space Agency, 2023).In this work 4 bands are used: Blue, Green Red and Near Infrared (NIR).These have a spatial resolution of 10 m and a maximum correction level called 2A (orthorectification with Bottom of Atmosphere correction).

Figure 1 .
Figure 1.Geo-localization of the study area: the map above is in equirectangular projection and WGS 84 ellipsoidal coordinates (EPSG: 4326); the map below is a satellite overview of the Bay of Pozzuoli (true colour RGB composition of Landsat 8 images) in UTM/WGS84 plane coordinates expressed in meters (EPSG: 32632).Finally, the last data necessary for training and testing the models is obtained from the ENC n.IT50083A (chart named as the Port of Pozzuoli), scale 1:7,500, produced by the Istituto Idrografico della Marina Militare (IIMM), the Italian Hydrographic Office.This chart format is digital, specifically vector, which complies with the standards established by the International Hydrographic Organization (IHO).As ENC it is an example of the database, standardized as to content, structure

Figure 2 .
Figure 2. The eastern sector of the Bay of Pozzuoli as resulting from RGB composition of Sentinel-2 images; the considered area is georeferred in UTM/WGS 84 plane coordinates expressed in meters (EPSG: 32632).ENCs are georeferred in WGS84 geodetic datum.The depths are referred to the Mean Lower Low Water (MLLW), which is determined by averaging each day's lowest tide at a particular location over a recording period (usually 19 years).ENC provides several information, such as shoreline shape and position, seafloor conformation, water depths, aids to navigation, dangers for navigation, anchorages, and other features relevant to navigation (IMO, 2019)(Alcaras et al., 2020).Among the various data provided by the considered nautical chart, only those relating to depth are selected, i.e., the contour lines and the soundings up to -26.50 m.Since the ENC is georeferred in WGS84 (ellipsoidal coordinates), we project the dataset in the Universal Transverse of Mercator (UTM)/WGS84 Zone 33 N (EPSG code: 32633).As shown in the figure3, we group the vertices of contour lines and the soundings in one shape file.

Table 1 .
Statistical values of the residuals generated for each model as difference between the true depth data extracted from the ENC and the corresponding satellite derived values.
(Najar et al., 2022)ot have a direct comparison given that this method is only introduced in 2023 and therefore after Said et al. study, provides greater accuracy.In the experiments conducted by Najar et al.(Najar et al., 2022), Sentinel-2 satellite images and bathymetric survey data are considered: to generate the bathymetric models of two areas, they experiment methods based on deep learning.Accuracy tests report RMSE equal to 7.59 m for the worst model and RMSE equal to 3.26 m for the best one.Also in this case, the results of our applications are very encouraging.
Their results are in line with those we acquire for BRM.Nevertheless, in our experiments the application of