A Hybrid Bio-Optical Transformation for Satellite Bathymetry Modeling Using Sentinel-2 Imagery

The article presents a new hybrid bio-optical transformation (HBT) method for the rapid modelling of bathymetry in coastal areas. The proposed approach exploits free-of-charge multispectral images and their processing by applying limited manpower and resources. The testbed area is a strait between two Greek Islands in the Aegean Sea with many small islets and complex seabed relief. The HBT methodology implements semi-analytical and empirical steps to model sea-water inherent optical properties (IOPs) and apparent optical properties (AOPs) observed by the Sentinel-2A multispectral satellite. The relationships of the calculated IOPs and AOPs are investigated and utilized to classify the study area into sub-regions with similar water optical characteristics, where no environmental observations have previously been collected. The bathymetry model is configured using very few field data (training depths) chosen from existing official nautical charts. The assessment of the HBT indicates the potential for obtaining satellite derived bathymetry with a satisfactory accuracy for depths down to 30 m.


Introduction
Multispectral imagery provides an alternative way to acquire marine environmental data for several operational, as well as scientific, purposes [1][2][3][4], in different time and spatial scales. Satellite imagery can provide users with adequate coverage and often free-of-charge environmental information, like those acquired by the Landsat 8 (OLI) and Sentinel-2 sensors [5]. Currently, the acquisition of spaceborne observations is a common practice, however, the processing of imagery data to extract bathymetry is a more complex procedure, since it requires special treatment for estimating the subsurface irradiance reflectance signal (R(0−)). In general, the range of extracted depths is strongly related to the water clarity, or the so-called ocean color (IOPs, turbidity, etc.) The main issue is that the ocean color depends directly on the optical properties of the water, which in turn vary spatially and temporally [6]. As a result, a special methodology is needed for correcting multispectral images due to atmospheric and water column interactions.
After pre-processing the imagery data with respect to atmospheric interference (i.e., via models such as 6SV, OPERA, etc.) [7,8], and air-water surface interaction, the remaining signal propagating

The Field Data
The training and validation bathymetry used in the Hybrid Bio-optical Transformation (HBT ) model implementation was extracted from the nautical chart No. 423/2 (scale 1:20,000), published in June 1989 by the Hellenic Navy Hydrographic Service. This data originates from hydrographic surveys performed between 1982 and 1987, categorized as CATZOC A2 (Category Zone Confidence A2), with horizontal accuracy of the depicted depths of ±20.0 m and vertical accuracy of ±1.6 m.

Satellite Data
In the present investigation, we used a Sentinel-2A Level-2A cloud-free multispectral image (consisting of 13 spectral channels) acquired on 12/07/2017 (local time 09:00:21). The Level-2A data was downloaded directly from the open-source ESA/Copernicus Science Hub, and is ortho-corrected to the bottom-of-atmosphere (BOA) reflectance. The ESA atmospheric processing consists of combining state-of-the-art techniques for performing atmospheric corrections, including cirrus clouds [13], customized to the Sentinel-2 environment. The former correction is representative of the atmospheric and climate conditions in the study area (i.e., the strait of Paros and Antiparos Islands).
The imagery was referenced onto the UTM cartographic projection and on the WGS84 geodetic reference system, with a spatial resolution of 10 m (for channels 2, 3, 4, and 8 at 490, 560, 665, and 842 nm, respectively); 20 m (for channels 5, 6, 7, 8a, 11, and 12 at 705, 740, 783, 865, 1610, and 2190 nm, respectively); and 60 m (for channels 1,9, and 10 at 443, 945, and 1375 nm, respectively). The 10 m granule channels were used for assessing the Sentinel-2A imagery at its highest spatial analysis. The data analysis, processing, and evaluation were carried out with the use of the ENVI Software, SNAP Desktop/ESA, ArcGIS Pro, IBM SPSS and EXCEL toolboxes.

Satellite Data
In the present investigation, we used a Sentinel-2A Level-2A cloud-free multispectral image (consisting of 13 spectral channels) acquired on 12/07/2017 (local time 09:00:21). The Level-2A data was downloaded directly from the open-source ESA/Copernicus Science Hub, and is ortho-corrected to the bottom-of-atmosphere (BOA) reflectance. The ESA atmospheric processing consists of combining state-of-the-art techniques for performing atmospheric corrections, including cirrus clouds [13], customized to the Sentinel-2 environment. The former correction is representative of the atmospheric and climate conditions in the study area (i.e., the strait of Paros and Antiparos Islands).
The imagery was referenced onto the UTM cartographic projection and on the WGS84 geodetic reference system, with a spatial resolution of 10 m (for channels 2, 3, 4, and 8 at 490, 560, 665, and 842 nm, respectively); 20 m (for channels 5, 6, 7, 8a, 11, and 12 at 705, 740, 783, 865, 1610, and 2190 nm, respectively); and 60 m (for channels 1, 9, and 10 at 443, 945, and 1375 nm, respectively). The 10 m granule channels were used for assessing the Sentinel-2A imagery at its highest spatial analysis. The data analysis, processing, and evaluation were carried out with the use of the ENVI Software, SNAP Desktop/ESA, ArcGIS Pro, IBM SPSS and EXCEL toolboxes.

Methodology-Model Development
Several techniques have been adopted in order to extract bathymetric information from multispectral imagery to date [14,15]. The initial steps for SDB analysis comprise: (a) radiometric, geometric and atmospheric correction, such as the Sentinel-2A Level-2A products evaluated in this work; (b) sub-setting the imagery to the study region, namely the Paros-Antiparos Strait; and (c) imposing a sea-land mask using the NIR band (B8) of Sentinel-2A to emphasize the marine features. Following the above pre-processing phases, the overall methodology implemented in this study is described in the work-flow chart of Figure 2, and analyzed further below. geometric and atmospheric correction, such as the Sentinel-2Α Level-2A products evaluated in this work; (b) sub-setting the imagery to the study region, namely the Paros-Antiparos Strait; and (c) imposing a sea-land mask using the NIR band (B8) of Sentinel-2A to emphasize the marine features. Following the above pre-processing phases, the overall methodology implemented in this study is described in the work-flow chart of Figure 2, and analyzed further below.

Sunglint Correction
The reflectance of the bottom-of-atmosphere (BOA) is a parameter representing how much light is reflected from the sea surface and in which direction. When imagery observations are made above the sea surface, the signal measured by satellites (i.e., the remote sensing signal/reflectance, Rrs) is the sum of the water leaving signal (i.e., water leaving radiance/reflectance) and the sea-surface reflected signal (sunglint).
Sunglint is related to the specular reflection of sun radiation due to wave undulation and seasurface foam. According to Goodman et al. [16], high-resolution images uncorrected for the sunglint effect may introduce an error in the water-depth estimation of up to 30%. Modelling the sunglint error was achieved in this work by applying the NIR band (B8) of Sentinel-2A, according to the algorithm (Equation (1))developed by Kay et al. [17]: where Degl(RBi) is the resulting deglint image band; is the Bottom-of-Atmosphere remote sensing corrected reflectance of Bi band (being either B2, B3, or B4); RNIR is the B8 band of Sentinel 2A; and minNIR is the minimum pixel value of the B8 band.
In order to perform sunglint removal, an area of deep water was selected, with depths more than 50 m. Statistical analysis between B2/B8, B3/B8, and B4/B8 bands resulted in slopes of 0.67 for B2 (blue band), 0.52 for B3 (green band), and 0.31 for B4 (red band), respectively (Figure 3), whereas the minNIR value in this case was 0.055. The processed image after the sunglint error removal is shown in Figure 4.

Sunglint Correction
The reflectance of the bottom-of-atmosphere (BOA) is a parameter representing how much light is reflected from the sea surface and in which direction. When imagery observations are made above the sea surface, the signal measured by satellites (i.e., the remote sensing signal/reflectance, R rs ) is the sum of the water leaving signal (i.e., water leaving radiance/reflectance) and the sea-surface reflected signal (sunglint).
Sunglint is related to the specular reflection of sun radiation due to wave undulation and sea-surface foam. According to Goodman et al. [16], high-resolution images uncorrected for the sunglint effect may introduce an error in the water-depth estimation of up to 30%. Modelling the sunglint error was achieved in this work by applying the NIR band (B8) of Sentinel-2A, according to the algorithm (Equation (1)) developed by Kay et al. [17]: where Degl(R Bi ) is the resulting deglint image band; R BOA Bi is the Bottom-of-Atmosphere remote sensing corrected reflectance of B i band (being either B2, B3, or B4); R NIR is the B8 band of Sentinel 2A; and minNIR is the minimum pixel value of the B8 band.
In order to perform sunglint removal, an area of deep water was selected, with depths more than 50 m. Statistical analysis between B2/B8, B3/B8, and B4/B8 bands resulted in slopes of 0.67 for B2 (blue band), 0.52 for B3 (green band), and 0.31 for B4 (red band), respectively (Figure 3), whereas the minNIR value in this case was 0.055. The processed image after the sunglint error removal is shown in Figure 4.

Water Column Correction
As solar radiation penetrates the sea surface, it interacts with the water column and its constituents. This kind of light propagation reduces its intensity before the radiation reaches the sea bottom [18]. The decreased irradiance energy, usually called light attenuation, is dependent on the degree of water absorption and scattering processes. Lyzenga (1981) [19], introduced an image-based method (water column correction) for compensating for the effect of light attenuation due to absorption and scattering caused by depth variability. In essence, this technique equalizes the spectral signature of the sea bottom by producing a depth-invariant index determined from pairs of spectral bands. The Lyzenga method can be used instead of estimating the seabed reflectance; this technique was implemented in the particular marine area of our study, which consists of the various seabed sediments and depth variability.

Water Column Correction
As solar radiation penetrates the sea surface, it interacts with the water column and its constituents. This kind of light propagation reduces its intensity before the radiation reaches the sea bottom [18]. The decreased irradiance energy, usually called light attenuation, is dependent on the degree of water absorption and scattering processes. Lyzenga (1981) [19], introduced an image-based method (water column correction) for compensating for the effect of light attenuation due to absorption and scattering caused by depth variability. In essence, this technique equalizes the spectral signature of the sea bottom by producing a depth-invariant index determined from pairs of spectral bands. The Lyzenga method can be used instead of estimating the seabed reflectance; this technique was implemented in the particular marine area of our study, which consists of the various seabed sediments and depth variability.
Three areas were selected (indicated with the letters (A), (B), and (C)) to correct the Sentinel-2A image for the water column effect, with depth ranges of (A) 8-15 m; (B) 2-9 m; and (C) 10-15 m ( Figure 5). Subsequently, the sunglint corrected reflectance of Sentinel-2A spectral bands (Degl(R Bi )) were re-modelled using natural logarithms (Equation (2)), in order for their variation to become linear: The ratio of attenuation coefficients K d (Bi) was then calculated between pairs of spectral bands, using data within each marine area of the same seabed sediment, but with different depths. The ratio pairs validated in this study were B2/B3, B2/B4, and B3/B4, with the ratio of each band pair determined by Equations (3) and (4): where σ ii , σ jj in Equation (4) are the variances of the deglinted reflectances of bands B i and B j (i.e., B2 and B3) and σ ij is the covariance between bands B i and B j . were re-modelled using natural logarithms (Equation (2)), in order for their variation to become linear: The ratio of attenuation coefficients Kd(Bi) was then calculated between pairs of spectral bands, using data within each marine area of the same seabed sediment, but with different depths. The ratio pairs validated in this study were B2/B3, B2/B4, and B3/B4, with the ratio of each band pair determined by Equations (3) and (4): where σii, σjj in Equation (4) are the variances of the deglinted reflectances of bands Bi and Bj (i.e., B2 and B3) and σij is the covariance between bands Bi and Bj. The depth-invariant index (DII) is a parameter that represents the image correction due to the water column interaction (i.e., absorption and scattering) and was computed for each band pair using Equation (5): The depth-invariant index (DII) is a parameter that represents the image correction due to the water column interaction (i.e., absorption and scattering) and was computed for each band pair using Equation (5): Figure 6 presents the comparison between the DII and the deglinted image, indicating the spectral contrast of the seabed texture, which is better in the case of the water column corrected image ( Figure 6, right).

Band Ratio Methodology
The next phase of model analysis consisted of identifying the relationship between the bio-optical parameters and sea-bottom reflectance. To examine this, tests were carried out using various spectral band ratios. Previous research [20,21] confirmed that the ratio between two spectral bands allows for more stable depth results when extracting bathymetry in coastal areas with variations in bottom albedo.
In the present study, we adopted the ratio algorithm technique as shown in Equation (6), with an innovation, since initial tests indicated that it is rather difficult to estimate the exact value for the factor n: where, z is the water depth, R i and R j are the reflectances of bands i (blue) and j (green), m 1 is a constant for scaling the ratio to depth, and m o is the offset for z = 0. The n-factor usually takes a fixed value between 500 and 10,000, ensuring that the logarithm remains positive under any seabed texture, and results in a linear correlation between ratio bands and field data. Consequently, all the possible band combinations or band ratio schemes were examined in order to estimate which one provided the best correlation to field data or-in our case-to depths derived from a nautical chart. Initial tests using ratio transformation (Equation (6)) indicated that the best results for Sentinel-2A (BOA reflectance) were obtained when implementing the ratio (Equation (7)) between the sunglint-corrected blue (Degl_B2) and green bands (Degl_B3): where x is the unknown term of the regression equation (Figure 7).
Remote Sens. 2019, 11, x FOR PEER REVIEW 7 of 17 Figure 6 presents the comparison between the DII and the deglinted image, indicating the spectral contrast of the seabed texture, which is better in the case of the water column corrected image ( Figure 6, right).

Band Ratio Methodology
The next phase of model analysis consisted of identifying the relationship between the biooptical parameters and sea-bottom reflectance. To examine this, tests were carried out using various spectral band ratios. Previous research [20][21] confirmed that the ratio between two spectral bands allows for more stable depth results when extracting bathymetry in coastal areas with variations in bottom albedo.
In the present study, we adopted the ratio algorithm technique as shown in Equation (6), with an innovation, since initial tests indicated that it is rather difficult to estimate the exact value for the factor n: where, z is the water depth, Ri and Rj are the reflectances of bands i (blue) and j (green), m1 is a constant for scaling the ratio to depth, and mo is the offset for z = 0. The n-factor usually takes a fixed value between 500 and 10,000, ensuring that the logarithm remains positive under any seabed texture, and results in a linear correlation between ratio bands and field data. Consequently, all the possible band combinations or band ratio schemes were examined in order to estimate which one provided the best correlation to field data or-in our case-to depths derived from a nautical chart. Initial tests using ratio transformation (Equation (6)) indicated that the best results for Sentinel-2A (BOA reflectance) were obtained when implementing the ratio (Equation (7)) between the sunglintcorrected blue (Degl_B2) and green bands (Degl_B3): where x is the unknown term of the regression equation (Figure 7). An initial test was undertaken to assess the robustness of the band correlation between chart depths and the logarithmic ratio of Equation (7). For this reason, a set of training depths was selected, comprising 28 points depicted on the Greek chart no. 423/2, with depths ranging between 0.4 and 41.0 m. The ratio of Equation (7)   The chart depths refer to a variety of seabed textures and depth ranges, since the training depths were chosen arbitrarily from different seabed textures; nevertheless, the above test revealed that there is a good correlation between Sentinel-2A imagery and chart depths. An initial test was undertaken to assess the robustness of the band correlation between chart depths and the logarithmic ratio of Equation (7). For this reason, a set of training depths was selected, comprising 28 points depicted on the Greek chart no. 423/2, with depths ranging between 0.4 and 41.0 m. The ratio of Equation (7) was implemented along with regression analysis and the relevant curve fitting to provide the corresponding determination coefficient values (R 2 ), obtaining values of 0.78 for linear, 0.84 for quadratic, and 0.85 for cubic fit (Figure 7).
The chart depths refer to a variety of seabed textures and depth ranges, since the training depths were chosen arbitrarily from different seabed textures; nevertheless, the above test revealed that there is a good correlation between Sentinel-2A imagery and chart depths.
Based on these results, we further investigated the effects of inherent (IOPs) and apparent optical parameters (AOPs) on depth modelling. In particular, we examined how it is possible to improve the outcome bathymetry of the ratio algorithm from data taken solely by spaceborne sensors, and using the IOPS and AOPs of coastal waters.
Further trials showed that close to the coastline the remote sensing reflectance of B3 (i.e., R rs (B3) green band) is greater than the equivalent reflectance of B2 (blue band) and B4 (red band). This fact revealed the existence of subsurface meadows at depths of approximately 0-8 m, which mainly consist of species such as Posidonia oceanica and Cymodocea nodosa, the latter being characteristic species of the marine flora in the study area [22]. An unsupervised classification methodology was then adopted in order to model bathymetry over various seabed compositions and by implementing as criteria the main IOPs and AOPs (i.e., the total backscattering coefficient, the total absorption coefficient, and the diffuse attenuation coefficient (K d )). The unsupervised classification of the entire study area resulted in five classes, namely, turbid waters, sparse seaweed coverage, dense seaweed coverage, shallow offshore waters without any flora, and deep waters without any flora (Figure 8). The adopted approach in this research took into account the local conditions of the study region with the aim of classifying the image according to its bio-optical water parameters. In order to model these water properties within the study area, the approach of Lee et al. [23] was adopted.

Computation of Bio-Optical Parameters
The adopted approach in this research took into account the local conditions of the study region with the aim of classifying the image according to its bio-optical water parameters. In order to model these water properties within the study area, the approach of Lee et al. [23] was adopted.
Firstly, the back-scattering coefficient (b bp (λ)) related to suspended particles in the water column is computed using Equation (8); its application is based on the fact that since the tested area is entirely coastal, then the water optical properties are wavelength dependent. Hence, the most representative reference wavelength (λ o ) for the areas covered mainly by coastal suspended sediments is the Sentinel-2A band B4 at 665 nm: where α w (B4) = α(665) = 0.428 m −1 is the absorption coefficient of pure seawater given by Pope et al. [24], and R BOA rs (B4) is the corrected Bottom-of-Atmosphere remote sensing reflectance of the Sentinel-2A red band (B4) over water areas.
Knowing the sea surface reflectance R BOA rs (B i ) from the multispectral imagery, it is possible to model the remote sensing reflectance below the sea surface (r rs (λ)) (Equation (9)) for Sentinel-2A bands B2, B3, and B4 (at 10 m resolution): The subsurface reflectance (r rs (λ)) is a parameter of the u-ratio (Equation (10)), which in turn is a function of water IOPs, namely, the back-scattering coefficient (b b (λ)) and the absorption coefficient (α(λ)): The u(B2), u(B3), and u(B4) ratios are then calculated in order to connect the AOPs (i.e., reflectances K d (λ)) with the IOPs (i.e., b b (λ) and α(λ) coefficients) in the equivalent spectral bands of the Sentinel-2A satellite.
Subsequently, the power parameter Y is estimated (Equation (11)), which usually takes values between 0.1 and 2.5: where r rs (B1) and r rs (B4) are the subsurface remote sensing reflectance for bands (B1) and (B4), respectively. The total backscattering coefficient b b (λ) (Equation (12)) can then be computed using the following equation: In Equation (12), the backscattering coefficient (b bp (B4)) of the red band due to suspended particles is determined by Equation (8), the power factor Y is provided by Equation (11), and the backscattering coefficient of pure seawater (b bw (λ)) is obtained by interpolation from Morel measurements (Morel, 1974, [25]. The wavelength (λ) in Equation (9) corresponds to the wavelengths of B2 and B3 Sentinel-2A spectral bands at 490 and 560 nm, respectively.
Next, the total absorption coefficient α(λ) is computed using Equation (10) with Equation (13): where the total backscattering coefficient b b (λ) is given by Equation (12). Finally, the diffuse attenuation coefficient K d (λ) is determined from Equation (14): where m o 1 + 0.005 × θ s , in which θ s is the solar zenith angle; and m 1 , m 2 , and m 3 are model constants with values of 4.18, 0.52, and 10.8, respectively, according to Lee et al. (2002). Equation (14) is then transformed to Equation (15) as follows:

Seabed Classification Results
Experiments shown that the adopted technique (HBT) exploits IOPs and AOPs to perform an unsupervised seabed classification ( Figure 8) according to the relationships of their bio-optical characteristics, as modelled according to Table 1 below. Table 1. Bio-optical characteristics of seabed classes.

Shallow Waters
The statistical behavior of the above-computed bio-optical parameters in the study coastal area is presented in Table 2. The diffuse attenuation coefficient (K d ) of Sentinel-2A band B4 (K d (B4)) reaches minimum values in areas with either a high load of suspended particles or with coverage of sparse flora (meadows), due to increased scattering. Furthermore, in areas with suspended particles, the total absorption coefficient reaches a minimum in band B3 (α(B3)) and maximum in band B1 (α(B1)), respectively. The backscattering coefficients (b b (B i ) have a similar fluctuation in all the classified seabed types, reaching their maximum in band B1 and gradually decreasing to their minimum in band B4. In areas covered by submerged flora, the reflectance of the green band (B3) dominates, since it is the main parameter that characterizes these coastal regions. Offshore waters are purely defined by the combination of the rate of image sunglint correction in bands B1 and B2, and the variation of diffuse attenuation coefficient (K d ) in spectral bands B2 and B4.

Satellite Derived Bathymetry (SDB) Results
In order to examine the robustness of the HBT methodology, statistical tests were undertaken by considering as field observations the depths depicted on the official large-scale hydrographic chart (no. 423/2) that covers the coastal area of interest. A small number (i.e., <10) of map depths was used as training data for extracting seabed-classified SDB by implementing either the band deglint ratios (Equation (7)) or the ratios of the modelled bio-optical properties. Linear and quadratic regression analyses were then applied to the SDB results, according to which of the five seabed classes these results belonged to. There was no need for tide reduction in the SDB outcome, since the regression analysis was based on chart (reduced) data. Table 3 presents the best correlation (after removal of two-three erroneous data points, i.e., errors caused due to remaining glint on some pixels) obtained for each seabed class between, on one hand, the training depths and, on the other hand, the SDB retrieved either from the band deglint or the bio-optical properties ratios. We find that in areas of suspended particles, or dense or sparse flora, SDB exhibits best correlation with the training chart data when implementing the bio-optical properties ratios; over offshore waters, however, the best correlation between the SDB and the training map data is obtained when applying the band-deglint ratios. Overall, the quadratic regression analysis results in a better correlation compared to the linear regression for all seabed classes. The validation set comprised 20 depths that were also depicted on the official nautical chart covering the strait between Paros and Antiparos Islands, within the range from 0.4 to 41.00 m. These depths were arbitrarily chosen, and for all seabed classes contained some of the data used to configure the ratio transformation (Figure 7).
By comparing the coefficients of determination (R 2 ) of each seabed class, we find that the HBT bathymetry seems to be correlated better with the chart data at an average of approximately 13% when a linear regression is applied, and approximately 10% when a quadratic regression analysis is implemented. Depths are calculated with adequate accuracy up to 30 m.
Finally, the bathymetries extracted from each class were merged to create the total bathymetric model for the area of interest (Figure 9). For validation purposes, the HBT results were compared with the bathymetry derived from the ratio transformation. By comparing the coefficients of determination (R 2 ) of each seabed class, we find that the HBT bathymetry seems to be correlated better with the chart data at an average of approximately 13% when a linear regression is applied, and approximately 10% when a quadratic regression analysis is implemented. Depths are calculated with adequate accuracy up to 30 m.
Finally, the bathymetries extracted from each class were merged to create the total bathymetric model for the area of interest (Figure 9). For validation purposes, the HBT results were compared with the bathymetry derived from the ratio transformation.    By comparing the coefficients of determination (R 2 ) of each seabed class, we find that the HBT bathymetry seems to be correlated better with the chart data at an average of approximately 13% when a linear regression is applied, and approximately 10% when a quadratic regression analysis is implemented. Depths are calculated with adequate accuracy up to 30 m.
Finally, the bathymetries extracted from each class were merged to create the total bathymetric model for the area of interest (Figure 9). For validation purposes, the HBT results were compared with the bathymetry derived from the ratio transformation.

Discussion
The hybrid bio-optical transformation (HBT) methodology was presented for extracting SDB from free-of-charge Sentinel-2A multispectral imagery and based on a minimum number of training chart data. HBT utilizes both empirical approaches and semi-analytical algorithms for calculating the main seawater optical parameters. Sentinel-2A data were initially corrected for geometrical distortions, radiometric errors, atmosphere interaction, sunglint noise, and water column interference.
According to Hedley et al. [26] the glint correction of Sentinel-2 data is potentially more challenging than with higher spatial resolution imagery (e.g., WorldView-2) because 10 m pixels may contain a combination of glint at different scales that may be related to surface waves at more or less than 10 m apart. Following this pre-processing, the seawater bio-optical properties were estimated in order to classify the imagery into seabed types of similar optical characteristics. Since there were no previous environmental observations in the region of interest, a semi-analytical model was adopted as the best way to compute IOPs and AOPs.
The HBT method does not require previous knowledge of the bottom albedo for the calculation of both the intrinsic and apparent water optical properties; this seems to be an important issue, since no costly and time-consuming field measurements are needed for the collection of the spectral seabed data. Nevertheless, when missing measurements of IOPs, their uncertainties are known to intrude and propagate into their calculation algorithms. The major sources of errors and uncertainties in deriving IOP products are the uncertainties in the R rs (λ) data, due to sensor noise and atmospheric correction error, along with the assumptions and approximations of the IOP parameterization schemes.
Another issue for consideration is the absorption coefficient α w of pure seawater (Equation (8)), for which we used the values provided by Pope and Fry [24]. Most recent laboratory measurements in the ultra-violet and visible wavelengths by exhibit lower values in the UV spectrum, though Lee et al. [23] derived higher values for seawater in the UV spectrum, possibly due to absorption by dissolved inorganic constituents. These results imply that additional research is required to further bridge the gap among different measuring methods of α w .
Werdell et al. [27] summarize the up-to-date approaches for retrieving marine inherent optical properties from ocean color remote sensing. A variety of techniques for potentially improving IOP retrieval exists, depending also on the sensor, computational capacity, and water types. Several IOP optimization schemes have also been implemented; for example, keeping the power parameter Y in the model of Lee et al. as a free variable seems to improve the errors relative to the standard model for most depths, though retrieved IOPs are not necessarily correct. Hedley et al. [26] proposed an adaptive look-up tables (ALUT) inversion approach for the application of semi-analytical models, which facilitates the uncertainty propagation; it was found that, under suitable conditions, a strong relationship between Sentinel-2A SDB and in situ depths exists, especially when combined with mean filtering techniques. Nevertheless, LUT approaches remain computationally expensive, whereas, in the HBT method, limited resources are required.
After completing the calculations of the bio-optical parameters, the relationships between the main IOPs and AOPs (b b (λ), α(λ), and K d (λ)) were identified, and were used to classify the area of interest in sub-regions of water with similar bio-optical properties. Five water classes were found in the study area by performing an unsupervised classification based on the spectral response of various seabed and water types. This classification was determined exclusively according to the relationships of IOPs and AOPs of each seawater class (Table 1), meaning that the accuracy of classification did not depend on the accuracy of IOPs and AOPs. For this reason, the HBT method can be considered more stable in terms of classification capacity compared to other algorithms, since HBT is applied to water classes of homogeneous physical and spectral bottom characteristics. Following classification, the bathymetry in every seabed class and/or water type was retrieved using the best correlation between the logarithm ratios of band reflectance or the logarithm ratios between the intrinsic and the seawater AOPs for each class. Tests revealed that the HBT method can extract SDB by involving a small number of training data points (less than 10 depths in each class), and both training and validation depths were selected arbitrarily from a nautical chart; this procedure is particularly applicable in areas where insufficient field data exists, or depth soundings remain difficult to collect in situ.
Finally, the HBT SDB was compared with the results from the well-known ratio transformation method [22], using the same data set in both the configuration and validation procedures. SDB extracted using the HBT model yielded better results (RMSE = 2.0 m, R 2 = 0.961) against the chart data, regardless of the seabed type and for depths up to 30 m, compared to the ratio transformation method (RMSE = 5.4 m, R 2 = 0.903) at depths between 15 and 20 m. The HBT bathymetric model was executed via the SNAP Desktop, i.e., ESA's software for processing imagery downloaded from the Copernicus hub, and the ratio algorithm is available on the ENVI software.
Overall, the HBT advantages can be summarized as the following: (i) The SDB from the HBT model is not a function of any empirical factor, such as the n parameter of the ratio transform, or a 0 , a 1 . . . a 4 coefficients of a linear transform; it was confirmed that the use of ratios between bands or between their logarithms is a significant tool for bathymetry modelling, since both the HBT model, as well as the ratio transform, exploit the ratios between the spectral bands for extracting SDB from multispectral imagery. (ii) The HBT model may cover water bio-optical parameters with different flora and bottom spectral reflectance. (iii) The bathymetric information extracted requires limited resources to run and can be executed with appropriate free software; there is no requirement for extensive computations or allocation of expensive resources (computing systems, optimization software), nor the collection of costly field data.
The use of SDB data is especially useful where vessel traffic or the deployment of aids to navigation indicate that charted data is misleading and there have been no recent surveys to update the charts. Nevertheless, SDB is used by many hydrographic offices as a reconnaissance survey only [28]. Hydrographic offices use SDB data for generating depth curves and not to extract smooth-sheet soundings. This categorizes the SDB data as CATZOC C [29]. However, these limitations may be overcome by the employment of new methods, such as HBT and machine learning, providing more reliable SDB data that will meet the relevant IHO specifications for chart production.

Conclusions
The aim of the proposed hybrid bio-optical transformation (HBT) model is the use of minimum manpower, along with low-cost chart data in order to extract SDB from free-of-charge multispectral imagery. The preliminary encouraging results indicate the potential of the HBT method to be further extended into future trials over areas of high importance for navigation safety, but also in other marine applications, such as coastal management, and environmental and climate change monitoring. Given the scientific and commercial interest that SDB has been attracting, the proposed HBT bathymetric model may contribute to the Marine Spatial Data Infrastructure (MSDI), according to the recommendations and guidelines of the IHO (International Hydrographic Organization, [11]) and the EMSA (European Maritime Safety Agency, [29]).
The preliminary results show that the spectral attributes of seabed benthic communities are related directly with the spectral characteristics of physical and optical properties of the supernatants within the water column. In addition, the seabed classification can feasibly be accomplished by analyzing the bio-optical parameters. Hence, it is possible to link temporally and spatially the seawater optical properties with other marine physical features, such as seabed sediment type and/or bathymetry in the area of interest. The conducted tests showed that HBT extracted bathymetry contributes to rapid depth modeling and promises to address the difficulties currently encountered in the rapid collection of bathymetric data. In this framework, there is a need to further test the hybrid bio-optical transform (HBT) model by applying optimization techniques during the model calibration/training phase, in combination with machine learning classification techniques or artificial intelligence techniques.
HBT aims to address the difficulties currently encountered in the rapid collection of bathymetric data, when priority is given to the use of limited human and technical resources, as well as to the short completion, production, and updating of existing nautical charts. In order to develop practical rules of thumb for the updating of nautical charts, according to the relevant specifications of the IHO, there is a clear need to perform further trial tests in other coastal areas, including error analysis, various satellite sensors, and comparisons against in situ data collected by modern systems, such as Lidar, or multibeam sonar. Future plans include such trials and the comparison of their outcome with the international standards of IHO.
Author Contributions: A.K.M. conceived of the presented Bathymetric Transformation Modelling, performed the computations and presents the results; E.O. evaluated the analytical method and supervised the findings of this work; A.P. encouraged A.K.M. to continue the experiments and contributed to the interpretation of the results focusing on the aspects of the Safety of Navigation; S.P. supervised the findings of this work and contributed to the editing of the manuscript; All authors discussed the results and contributed to the final manuscript.